text:
!--------------------------------------------------------------------------------
! Copyright (c) 2016 Peter Grünberg Institut, Forschungszentrum Jülich, Germany
! This file is part of FLEUR and available as free software under the conditions
! of the MIT license as expressed in the LICENSE file in more detail.
!--------------------------------------------------------------------------------

MODULE m_od_vvac
CONTAINS
  SUBROUTINE od_vvac(&
       &     stars,vacuum,cell,&
       &     psq,rht,&
       &     vz)

    !     subroutine which calculates the non warped part of the
    !     vacuum potential (m=0,gz=0)
    !                               Y. Mokrousov
    !     the potential in this subroutine can be defined in two
    !     equivalent ways, which nevertheless give a bit defferent
    !     results, 2nd one seems to be more precise

    USE m_qsf
    USE m_od_cylbes
    USE m_types
    USE m_constants
    IMPLICIT NONE
    TYPE(t_vacuum),INTENT(IN)   :: vacuum
    TYPE(t_stars),INTENT(IN)    :: stars
    TYPE(t_cell),INTENT(IN)     :: cell

    COMPLEX, INTENT (IN) :: psq(stars%ng3)
    REAL,    INTENT (IN) :: rht(:,:) !(vacuum%nmzd,2)
    REAL,    INTENT (OUT) :: vz(:,:) !(vacuum%nmzd,2)

    COMPLEX  rhobar
    INTEGER  k1,k2,irec3,irec2,i,j,ivac,imz,imz1
    REAL     g2 ,a(vacuum%nmzd)
    REAL     fJ,z,zp,phi
    REAL     rht1(vacuum%nmzd)
    REAL     f2(vacuum%nmzd),f22(vacuum%nmzd)

    INTRINSIC cmplx

    DO i = 1,vacuum%nmz
       f2(i) = 0.
       f22(i) = 0.
       DO ivac = 1,vacuum%nvac
          vz(i,ivac) = 0.
       END DO
    END DO

    rhobar = -psq(1)

    DO  k1 = -stars%mx1,stars%mx1
       DO  k2 = -stars%mx2,stars%mx2
          irec3 = stars%ig(k1,k2,0)
          IF (irec3.NE.0) THEN
             irec2 = stars%ig2(irec3)
             IF (irec2.NE.1) THEN
                g2 = stars%sk2(irec2)
                phi = stars%phi2(irec2)
                CALL od_cylbes(1,cell%z1*g2,fJ)
                rhobar = rhobar - 2.*psq(irec3)*CMPLX(fJ/(g2*cell%z1),0.0)

             END IF
          END IF
       ENDDO
    ENDDO
    !----> 1st equivalent way

    DO  i=1,vacuum%nmz
       rht1(i) = fpi_const*(cell%z1+(i-1)*vacuum%delz)*rht(i,1)
    ENDDO
    CALL qsf(vacuum%delz,rht1(1),f2(1),vacuum%nmz,1)

    DO  i = 1,vacuum%nmz
       f2(i) = tpi_const*cell%z1*cell%z1*rhobar-f2(i)
    ENDDO

    DO  i = 1,vacuum%nmz
       DO  j = 1,vacuum%nmz
          IF (j.LT.i) THEN
             f22(j) = 0.0
          ELSE
             f22(j) = f2(j)/(cell%z1+vacuum%delz*(j-1))
          END IF
       ENDDO
       CALL qsf(vacuum%delz,f22(1),a,vacuum%nmz,0)
       DO  ivac =1,vacuum%nvac
          vz(i,ivac) = -a(1)
       ENDDO
    ENDDO
    !----> 2nd equivalent way (via the Green function)

    DO imz = 1,vacuum%nmz
       z = cell%z1 + (imz-1)*vacuum%delz
       DO imz1 = 1,vacuum%nmz
          zp = cell%z1 +  (imz1-1)*vacuum%delz
          IF (imz1.LE.imz) THEN
             rht1(imz1) = fpi_const*LOG(z)*zp*rht(imz1,1)
          ELSE
             rht1(imz1) = fpi_const*LOG(zp)*zp*rht(imz1,1)
          END IF

       END DO
       CALL qsf(vacuum%delz,rht1,a,vacuum%nmz,0)
       vz(imz,1) = tpi_const*LOG(z)*(cell%z1*cell%z1)*rhobar - a(1)
    END DO

    RETURN
  END SUBROUTINE od_vvac
END MODULE m_od_vvac

max_stars_repo_name: MRedies/FLEUR

max_issues_repo_name: MRedies/FLEUR

max_forks_repo_name: MRedies/FLEUR


text:
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++!
!                          Futility Development Group                          !
!                             All rights reserved.                             !
!                                                                              !
! Futility is a jointly-maintained, open-source project between the University !
! of Michigan and Oak Ridge National Laboratory.  The copyright and license    !
! can be found in LICENSE.txt in the head directory of this repository.        !
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++!
PROGRAM testFileType_Log
#include "UnitTest.h"
USE ISO_FORTRAN_ENV
USE UnitTest
USE ExceptionHandler
USE FileType_Fortran
USE FileType_Log

IMPLICIT NONE

CHARACTER(LEN=256) :: string
TYPE(ExceptionHandlerType),TARGET :: e
TYPE(FortranFileType) :: testFile
TYPE(LogFileType) :: testLogFile

CREATE_TEST('FILETYPE_LOG')

CALL e%setStopOnError(.FALSE.)
CALL e%setQuietMode(.TRUE.)

REGISTER_SUBTEST('LogFileType',testLogFileType)

FINALIZE_TEST()
!
!===============================================================================
CONTAINS
!
!-------------------------------------------------------------------------------
SUBROUTINE testLogFileType()
  LOGICAL :: bool
  CALL testLogFile%e%addSurrogate(e)
  ASSERT(.NOT.(testLogFile%isEcho()),'%isEcho()')
  CALL testLogFile%setEcho(.TRUE.)
  ASSERT(testLogFile%isEcho(),'%setEcho(...)')
  CALL testLogFile%initialize(UNIT=66,FILE='./test.log')
  ASSERT(TRIM(testLogFile%getFileName()) == 'test','%initialize(...)')
  CALL testLogFile%fopen()
  CALL testLogFile%message('Passed: CALL testLogFile%message(...)',.TRUE.,.TRUE.)
  CALL testLogFile%message('Passed: CALL testLogFile%message(...)',.FALSE.,.TRUE.)
  CALL testLogFile%clear(.FALSE.)
  ASSERT(testLogFile%getUnitNo() == -1,'%clear')
  CALL testFile%initialize(UNIT=66,FILE='./test.log',STATUS='OLD',ACTION='READ')
  CALL testFile%fopen()
  READ(66,'(a)') string
  bool=(TRIM(string(12:LEN(string))) == ' Passed: CALL testLogFile%message(...)')
  ASSERT(bool,'%message(...)')
  CALL testFile%clear(.TRUE.)
ENDSUBROUTINE testLogFileType
!
ENDPROGRAM testFileType_Log

max_stars_repo_name: wgurecky/Futility

max_issues_repo_name: wgurecky/Futility

max_forks_repo_name: wgurecky/Futility


text: #redirect Marrone Bio Innovations

ext: f

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
subroutine AverageToOneDimension(a,AvgDime,array,wopos,istep,ttime,Mesh,Memor)
   use typre
   use Mod_Element
   use Mod_Memor
   use Mod_Mesh
   use MPI
   use Mod_Postpr
   implicit none
   class(PostprFile) :: a
   integer(ip) :: AvgDime
   real(rp) :: array(:,:)
   character(*), intent(in)  :: wopos
   integer(ip) :: istep
   real(rp) :: ttime
   class(FemMesh)   :: Mesh
   type(MemoryMan) :: Memor

   integer(ip) :: ndofn

   real(rp) :: range(2,3),subcoord
   integer(ip) :: nsubdivisions
   real(rp)    :: sublength

   integer(ip) :: nelem,ielem,isub,maxnode,minnode,maxpoin,minpoin,minsub,maxsub,insub,ipoin
   real(rp)    :: maxcoord,mincoord,dvol,alpha

   class(FiniteElement), pointer :: e => NULL()
   real(rp), allocatable :: AuxSubArray(:,:)
   real(rp),   pointer   :: coord(:) => NULL()

   integer(ip) :: npoin,gnpoin

   integer(ip) :: gminmaxsub(2)
   integer(ip), allocatable :: rGminmaxsub(:,:)
   real(rp),    allocatable :: rbuffer(:,:)

   !MPI
   integer, parameter :: mtag1 = 1, mtag2 = 2,mtag3 = 3
   integer status(MPI_STATUS_SIZE),irequest1,irequest2
   integer :: ierr
   integer(ip) :: irank

   integer(ip) :: npoinLocal,local_nodes,inode,ndime

   !Open the graf file if I did not do it
   call OpenGrafFileTrue(a)

   ndofn = size(array,1)

   !GetDimensions
   call Mesh%GetGRange(range)
   call Mesh%GetGNpoin(gnpoin)
   call Mesh%GetNpoinLocal(npoinLocal)
   call Mesh%GetNdime(ndime)
   call Mesh%GetNelem(nelem)

   !We divide the AvgDime in how many parts? Depending on the number of points
   !We want 10 subdivisions per point, just to make sure everything is fine
   nsubdivisions = 10.0_rp*(gnpoin**(1.0_rp/ndime))
   sublength = (range(2,Avgdime)-range(1,Avgdime))/nsubdivisions

   !Memor allocations
   call Memor%alloc(ndofn+1,nsubdivisions+1,AuxSubArray,'AuxSubArray','AverageToOneDimension')
   call Mesh%ElementAlloc(e,Memor,'ForceClosedRule','vmass')
   if (Mesh%MPIrank == Mesh%MPIroot) then
      call Memor%alloc(2,Mesh%MPIsize,rGminmaxsub,'rGminmaxsub','AverageToOneDimension')
      call Memor%alloc(ndofn+1,nsubdivisions+1,rbuffer,'rbuffer','AverageToOneDimension')
   endif

   !Minimum and maximum subdivision contribution for every processor, we send only this range
   !of subdivisions
   gminmaxsub(1) = nsubdivisions+2
   gminmaxsub(2) = -1

   do ielem = 1,nelem
      call Mesh%ElementLoad(ielem,e)
      call e%elmdel

      dvol = e%detjm

      !Now I need to add the corresponding part to each subelement contribution
      maxcoord = maxval(e%elcod(Avgdime,1:e%pnode))
      maxnode  = maxloc(e%elcod(Avgdime,1:e%pnode),1)
      maxpoin  = e%lnods(maxnode)

      mincoord = minval(e%elcod(Avgdime,1:e%pnode))
      minnode  = minloc(e%elcod(Avgdime,1:e%pnode),1)
      minpoin  = e%lnods(minnode)

      !I need to add to all the subdivisions
      minsub = ceiling((mincoord-1e-8-range(1,Avgdime))/sublength)+1
      maxsub = floor((maxcoord+1e-8-range(1,Avgdime))/sublength)+1
      if (maxsub > nsubdivisions+1) maxsub = nsubdivisions+1
      insub = maxsub-minsub+1

      gminmaxsub(1) = min(gminmaxsub(1),minsub)
      gminmaxsub(2) = max(gminmaxsub(2),maxsub)

      !Count how many of lcal nodes and add only the proportional part
      local_nodes = 0
      do inode = 1,e%pnode
         if (e%lnods(inode) <= npoinLocal) then
            local_nodes = local_nodes+1
         endif
      enddo

      dvol = dvol/insub*(1.0_rp*local_nodes)/(1.0_rp*e%pnode)

      do isub = minsub,maxsub
         subcoord = range(1,Avgdime)+(isub-1)*sublength

         alpha = (subcoord-mincoord)/(maxcoord-mincoord)

         !Mass
         AuxSubArray(1,isub) = AuxSubArray(1,isub) + dvol
         !RHS
         AuxSubArray(2:1+ndofn,isub) = AuxSubArray(2:1+ndofn,isub)+dvol*(array(:,maxpoin)*alpha+array(:,minpoin)*(1-alpha))
      enddo

   enddo

   !We gather how many results root needs to receive
   call MPI_Gather( gminmaxsub, 2, MPI_INTEGER4, rGminmaxsub,2, MPI_INTEGER4, Mesh%MPIroot, Mesh%MPIcomm,ierr);

   !Now we send the results root needs and root receives
   call MPI_ISEND(AuxSubArray(1:ndofn+1,gminmaxsub(1):gminmaxsub(2)),(ndofn+1)*(gminmaxsub(2)-gminmaxsub(1)+1), MPI_REAL8, Mesh%MPIroot, mtag1, Mesh%MPIcomm,irequest1, ierr)
   !And receive
   if (Mesh%MPIrank == Mesh%MPIroot) then
      do irank = 0,Mesh%MPIsize-1
         call MPI_RECV(rbuffer,(ndofn+1)*(rGminmaxsub(2,irank+1)-rGminmaxsub(1,irank+1)+1),MPI_REAL8,irank,mtag1,Mesh%MPIcomm,status,ierr)
         AuxSubArray(:,rGminmaxsub(1,irank+1):rGminmaxsub(2,irank+1)) = AuxSubArray(:,rGminmaxsub(1,irank+1):rGminmaxsub(2,irank+1)) + &
                                                                  rbuffer(:,1:rGminmaxsub(2,irank+1)-rGminmaxsub(1,irank+1)+1)
      enddo
   endif
   !We all wait for the send to be complete
   call MPI_WAIT(irequest1, status, ierr)

   if (Mesh%MPIrank == Mesh%MPIroot) then
      !We divide by the mass to have the mean value
      do isub = 1,nsubdivisions+1
         AuxSubArray(2:1+ndofn,isub) = AuxSubArray(2:1+ndofn,isub)/AuxSubArray(1,isub)
      enddo

      write(a%lun_outpu_graf,101) adjustl(trim(wopos)),istep,ttime, adjustl(trim(wopos)),istep,ttime
      do isub = 1,nsubdivisions+1
         write(a%lun_outpu_graf,*) (isub-1)*sublength, AuxSubArray(2:1+ndofn,isub)
      enddo
   endif

   !Memory deallocations
   call Memor%dealloc(ndofn+1,nsubdivisions+1,AuxSubArray,'AuxSubArray','AverageToOneDimension')
   call Mesh%ElementDeAlloc(e,Memor,'ForceClosedRule','vmass')
   if (Mesh%MPIrank == Mesh%MPIroot) then
      call Memor%dealloc(2,Mesh%MPIsize,rGminmaxsub,'rGminmaxsub','AverageToOneDimension')
      call Memor%dealloc(ndofn+1,nsubdivisions+1,rbuffer,'rbuffer','AverageToOneDimension')
   endif

!    call MPI_BCAST(AuxSubArray, (ndofn+1)*(nsubdivisions+1), MPI_REAL8, Mesh%MPIroot, a%MPIcomm, ierr)
!    !Now we put the Averaged value in each of the nodes
!    call Mesh%GetNpoin(npoin)
!
!    do ipoin = 1,npoin
!       call Mesh%GetPointCoord(ipoin,coord)
!       subcoord = coord(AvgDime)
!
!       isub = (subcoord-range(1,AvgDime))/(range(2,AvgDime)-range(1,AvgDime))*nsubdivisions+1
!       if (isub > nsubdivisions) isub = nsubdivisions
!       alpha = (subcoord-(range(1,AvgDime)+(isub-1)*sublength))/(sublength)
!
!       AvgArray(1:ndofn,ipoin) = AuxSubArray(2:ndofn+1,isub+1)*alpha + AuxSubArray(2:ndofn+1,isub)*(1-alpha)
!    enddo
!
!    !Finally we communicate the ghosts
!    call Mesh%ArrayCommunicator%GhostCommunicate(ndofn,AvgArray(1:ndofn,1:npoin))

   101 format ('# Graph: ''',a22, i6,f12.6,''' ',/,&
         '#'                                                    ,/,&
         '# X: ''Distance'' Y: ''',a22, i6,f12.6,''' '     ,/,&
         '#'                                                    ,/,&
         '# Units: {} {} ')

end subroutine

subroutine AverageToOneDimension1D(a,AvgDime,array,wopos,istep,ttime,Mesh,Memor)
   use typre
   use Mod_Memor
   use Mod_Mesh
   use Mod_Postpr
   implicit none
   class(PostprFile) :: a
   integer(ip) :: AvgDime
   real(rp), target :: array(:)
   character(*), intent(in) :: wopos
   integer(ip) :: istep
   real(rp) :: ttime
   type(FemMesh)   :: Mesh
   type(MemoryMan) :: Memor

   interface
      subroutine AverageToOneDimension(a,AvgDime,array,wopos,istep,ttime,Mesh,Memor)
         use typre
         use Mod_Memor
         use Mod_Mesh
         use Mod_Postpr
         implicit none
         class(PostprFile) :: a
         integer(ip) :: AvgDime
         real(rp) :: array(:,:)
         character(*), intent(in)  :: wopos
         integer(ip) :: istep
         real(rp) :: ttime
         type(FemMesh)   :: Mesh
         type(MemoryMan) :: Memor
      end subroutine

   end interface

   real(rp), pointer :: aux_array(:,:) => NULL()
   integer(ip) :: array_size

   array_size = size(array)

   aux_array(1:1,1:array_size) => array

   call AverageToOneDimension(a,AvgDime,aux_array,wopos,istep,ttime,Mesh,Memor)

end subroutine

max_stars_repo_name: ciaid-colombia/InsFEM

max_issues_repo_name: ciaid-colombia/InsFEM

max_forks_repo_name: ciaid-colombia/InsFEM


! { dg-do run }
! PR libfortran/31210
program test
  implicit none
  integer :: l = 0
  character(len=20) :: s

  write(s,'(A,I1)') foo(), 0
  if (trim(s) /= "0") STOP 1

contains

  function foo()
    character(len=l) :: foo
    foo = "XXXX"
  end function

end program test


! Check that the status of a failed image is retrieved correctly.
! Image two is to fail, all others to continue.

program image_fail_and_status_test_1
  use iso_fortran_env , only : STAT_FAILED_IMAGE, STAT_STOPPED_IMAGE
  implicit none
  integer :: i, stat

  associate(np => num_images(), me => this_image())
    if (np < 3) error stop "I need at least 3 images to function."
    do i= 1, np
      if (image_status(i) /= 0) error stop "image_status(i) should not fail"
    end do

    ! Need to sync here or above image_status might catch a fail already.
    sync all
    if (me == 2) fail image
    sync all (STAT=stat)
    ! Check that all images returning from the sync report the failure of an image
    print *,"sync all (STAT=", stat, ")"
    if (stat /= STAT_FAILED_IMAGE) error stop "Expected sync all (STAT == STAT_FAILED_IMAGE)."

    do i= 1, np
      stat = image_status(i)
      if (i /= 2 .AND. stat /= 0 .AND. stat /= STAT_STOPPED_IMAGE) error stop "image_status(i) should not fail"
      if (i == 2 .AND. stat /= STAT_FAILED_IMAGE) error stop "image_status(2) should report fail"
    end do

    if (me == 1) print *,"Test passed."
  end associate

end program image_fail_and_status_test_1


text:
/*****************************************************************************
 * wb_test_master_test_base.f
 *****************************************************************************/

+UVM_TESTNAME=wb_test_master_test_base


text:
      subroutine rngint(ncrop,nowcrp)
c
c     NOTE: BECAUSE "GDAY" IS USED TO DETERMINE WHICH EQUATION IS USED
c           FOR DEFAULT VALUES OF "TLIVE", "TILLAGE" WHICH READS GDAY
c           *MUST* ALWAYS BE CALLED BEFORE "RNGINT".
c
c   ***************************************************************
c   * NOTE: ALL SOURCE FILES THAT REFERENCE "CRINPT1.INC" MUST BE *
c   *       COMPILED WITH THE *NEW* CRINPT1.INC TO WORK WITH THIS *
c   *       FILE.                                                 *
c   ***************************************************************
c
c   ***************************************************************
c   * NOTE: ALL SOURCE FILES THAT REFERENCE "CRINPT3.INC" MUST BE *
c   *       COMPILED WITH THE *NEW* CRINPT3.INC TO WORK WITH THIS *
c   *       FILE.                                                 *
c   ***************************************************************
c
c   ***************************************************************
c   * NOTE: ALL SOURCE FILES THAT REFERENCE "CROUT.INC" MUST BE   *
c   *       COMPILED WITH THE *NEW* CROUT.INC TO WORK WITH THIS   *
c   *       FILE.                                                 *
c   ***************************************************************
c
c     Note: Need definitions of PSCDST & PSCD2S in include file CRINPT2.INC.
c           CRM -- 4/09/93.
c
c     + + + PURPOSE + + +
c     Initializes variables used in RANGE, when land use is rangeland.
c     Called by OFE from CONTIN.
c
c     Called from CONTIN
c     Author(s): Weltz, Meyer
c     Reference in User Guide: Chapter 8
c
c     Changes:
c          1) Deleted references to common block UPDATE.
c          2) Deleted local variables SHEEP & NOWRES.
c          3) Added SPAI, & TPAI to CRINPT3.INC permit single initial
c             calculation of SPAI & TPAI in RNGINT, rather than daily
c             computation in PATRIB.  This required a change in the
c             logic: a loop to initialize each plant type.  Deleted
c             the appropriate code in PATRIB.
c          4) Deleted generic SAVE of all variables.
c          5) Dimensioned rooty(MXNSL) to rooty(MXNSL,MXPLAN) in CROUT.INC.
c             Re: 5/6/93 conversation w/Weltz. -- CRM
c          6) Subscripts were reversed:
c                 if(gday(i,1,iplane) .gt. 0) then
c             changed to:
c                 if(gday(1,i,iplane) .gt. 0) then
c          7) Did a major overhaul relative to the variable PLIVE.  In
c             its original state, PLIVE would not handle more than one
c             OFE correctly since it is modified if burning or herbicide
c             are used.  It also would not supply a default value for a
c             plant that was grazed on one OFE and not grazed on another.
c             Since modules need to be able to swap between plant types
c             on an OFE during a simulation, PLIVE(NTYPE) was redimen-
c             sioned to PLIVE(NTYPE,IPLANE).  This permits an OFE to be
c             altered without affecting the other OFE's.  CRINPT1.INC
c             was altered to accommodate this.  -- CRM -- 5/13/93.
c          8) PROOT was altered analogous to PLIVE. -- CRM -- 5/14/93.
c          9) ROOT(MXPLAN) was removed from common CROUT.INC and made
c             a local variable only - since it is computed and used now
c             only in this routine.  dcf  5/19/93
c         10) statement setting rdf(plant) to 0.43 moved to top of
c             code - otherwise it is undefined at first point of use.
c         11) variable rdf was a constant value of 0.43, and was only used
c             in this subroutine.  Thus it was removed from common CRINPT1.INC
c             and uses in RNGINT were changed to values of 0.43.  dcf  5/21/93
c
c     Version: This module recoded from WEPP Version 92.25.
c     Date recoded: 04/07/93 - 4/09/93.
c                   05/06/93 - 5/13/93.
c     Recoded by: Charles R. Meyer.
c
c     + + + KEYWORDS + + +
c
c     + + + PARAMETERS + + +
      include 'pntype.inc'
      include 'pmxnsl.inc'
      include 'ptilty.inc'
      include 'pmxres.inc'
      include 'pmxpln.inc'
      include 'pmxgrz.inc'
      include 'pmxcrp.inc'
      include 'pmxtls.inc'
      include 'pmxtil.inc'
      include 'pmxhil.inc'
c
c     + + + ARGUMENT DECLARATIONS + + +
      integer ncrop, nowcrp
c
c     + + + ARGUMENT DEFINITIONS + + +
c     ncrop  - number of input crops (rangeland plant communities)
c     nowcrp - current crop
c
c     + + + COMMON BLOCKS + + +
      include 'ccrpgro.inc'
c      modify: ytn
c       write: daymin
c
      include 'ccrpvr1.inc'
c       write: rtm(mxres,mxplan)
c
      include 'ccrpvr2.inc'
c     modify: vdmt(mxplan), cf(ntype)
c
      include 'ccrpvr3.inc'
c      modify: bbb(ntype)
c
      include 'ccover.inc'
c       write: canhgt(mxplan),inrcov(mxplan),rilcov(mxplan)
c      modify: cancov(mxplan),gcover(mxplan)
c
      include 'ccrpout.inc'
c       write: lai(mxplan), rtd(mxplan)
c      modify: rtmass(mxplan),
c
      include 'ccrpprm.inc'
c        read: itype(mxcrop,mxplan)
c
      include 'cends.inc'
c      modify: rspace(mxplan) (?)
c
      include 'crinpt1.inc'
c      modify: ffk(ntype),plive(ntype,mxplan),proot(ntype,mxplan),
c              yield(mxplan)
c       write: pyield(mxplan)
c
c     include 'crinpt1a.inc'
c
      include 'crinpt2.inc'
c
      include 'crinpt3.inc'
c      modify: spai(ntype), tpai(ntype), ghgt(ntype), basden
c
      include 'crinpt4.inc'
c      modify: dwood(mxplan)
c       write: ptlive(mxplan),first(ntype),xlive(mxplan)
c
      include 'crinpt5.inc'
c        read: tdiam(ntype), tcoeff(ntype), tpop(ntype), thgt(ntype)
c
      include 'crout.inc'
c      modify: tlive(mxplan),slive
c       write: droot(mxnsl),rooty(mxnsl,mxplan)
c
      include 'crinpt6.inc'
c      modify: bcover(mxplan), rufcov(mxplan)
c
      include 'cstruc.inc'
c        read: iplane
c
      include 'cwater.inc'
c       write: watstr(mxplan)
c
      include 'cwint.inc'
c        read: deglat
c
c     + + + LOCAL VARIABLES + + +
      real gpai, totpai, tarea, prevrt, ch, h
      real tmpvr2, tmpvr3, tmpvr4, tmpvr5, tmpvr6, root(mxnsl)
      integer plant, i, j
c
c     + + + LOCAL DEFINITIONS + + +
c     gpai   - total projected area from the side, intercepted for grass &
c              herbaceous plants, per 100 meter transect; ie, vertical area
c              "blocked" by plants, measured perpendicular to the slope.
c              Used for snowdrift calculations.
c     totpai - total projected area from the side, for all plants combined
c     tarea  - area of the vertical rectangle circumscribed around the
c              largest plant type present (trees, shrubs, & herbs).
c     prevrt - root mass in the previous layer (one directly above)
c     ch     -
c     h      -
c
c     + + + END SPECIFICATIONS + + +
c
c
      plant = itype(nowcrp,iplane)
c
c     ---- root distribution coefficient for mass by depth
c     rdf(plant) = 0.43
c     NOTE - problems with use of RDF in code below - and since it is
c     always set to a constant value of 0.43 - change from use
c     of an array to a numerical constant in the appropriate
c     equations.    dcf  5/21/93
c
c     ---- Julian date peak standing crop occurs, for first growth peak
      pscdst(plant) = pscday(plant)
c     ---- Julian date peak standing crop occurs, for second growth peak
      pscd2s(plant) = scday2(plant)
c
c     From a user-supplied value of PLIVE for each plant, generate
c     an initial value of PLIVE for each plant and each OFE.  This
c     value may be modified later by RHERB or RBURN, and the plant
c     type may change on this OFE during the simulation.  Also, if
c     the user entered a zero for the current plant, estimate a
c     starting value of PLIVE which depends on whether the OFE is
c     grazed in the first year or not.
c
c     *** L1 IF ***
c     ---- Execute this code the first time it is encountered (OFE #1).
      if (iplane.eq.1) then
c       *** BEGIN L2 DO-LOOP ***
c       Initialize PLIVE values for each plant by OFE.
        do 30 i = 1, ncrop
c         Next line added because NCROP is not readily available to RNGINT.
          if (plive(i,1).ge.0.0) then
c
c           *** L3 IF ***
c           -------- user supplied an initial value
            if (plive(i,1).ne.0.0) then
              do 10 j = 1, mxplan
c               Next line added because NELMT is not readily available to RNGINT.
                if (plive(i,1).gt.0.0) then
c                 ------------ Set PLIVE for each OFE to user-input value.
                  plive(i,j) = plive(i,1)
c                 ------------ maximum potential root biomass
c                 (WEPP Equation 8.5.25)
c                 proot(i,j) = root10(i) / (10.0**rdf(i))
c                 NOTE - assuming that value for rdf=0.43 - using this to simplify equation
c                 dcf 5/21/93
c
                  proot(i,j) = root10(i) / 2.6915348
                end if
   10         continue
c
c           *** L3 ELSE ***
c           -------- model must generate initial value...
            else
c             ---------- set PLIVE for this plant for *all* OFE's
              do 20 j = 1, mxplan
c               Next line added because NELMT is not readily available to RNGINT.
                if (plive(i,1).gt.0.0) then
c                 ---------- check for grazing on this OFE the first year.
                  if (gday(1,i,j).gt.0) then
c                   -------------- grazed areas
c                   (WEPP Equation 8.5.6)
                    plive(i,j) = 0.3569 + 0.0036 * pptg(i)
c                 ------------ non-grazed areas
                  else
c                   (WEPP Equation 8.5.7)
                    plive(i,j) = 0.7723 + 0.0030 * pptg(i)
                  end if
c                 ------------ maximum potential root biomass
c                 (WEPP Equation 8.5.25)
c                 proot(i,j) = root10(i) / (10.0**rdf(i))
c                 NOTE - assuming that value for rdf=0.43 - using this to simplify equation
c                 dcf 5/21/93
c
                  proot(i,j) = root10(i) / 2.6915348
                end if
   20         continue
c           *** L3 ENDIF ***
            end if
          end if
c
c       *** END L2 DO-LOOP ***
   30   continue
c     *** L1 ENDIF ***
      end if
c
c
c     ---- Assume rooting depth (RTD) equals depth to bottom of last soil layer
      rtd(iplane) = solthk(nsl(iplane),iplane)
c
c
c     *******************************************
c     * Calculate root distribution by depth. *
c     *******************************************
c
c     ---- maximum potential root biomass
c     (WEPP Equation 8.5.25)
c     proot(plant) = root10(plant) / (10.0**rdf(plant))
c
      prevrt = 0.0
      do 40 i = 1, nsl(iplane)
c       ------ cumulative total roots to bottom of current layer
c       (WEPP Equation 8.5.24)
        rtmass(iplane) = rootf(plant) * proot(plant,iplane) * (100.0*
     1      solthk(i,iplane)) ** 0.43
c
c       1                         (100.0*solthk(i,iplane))**rdf(plant)
c       NOTE - since RDF is only used in SR RNGINT - removed from common,
c       and since value is always 0.43 for all plant types - use
c       this.  Problems in do 8 loop with values of rdf(i)
c       dcf 5/21/93
c
c       ------ root mass in current soil layer today
        root(i) = rtmass(iplane) - prevrt
c       ------ save root mass to bottom of this layer (top of next layer)
        prevrt = rtmass(iplane)
   40 continue
c
      if (rtmass(iplane).gt.0.0) then
        do 50 i = 1, nsl(iplane)
c         -------- fraction of rootmass in each soil layer (constant throughout sim'n.)
c         (WEPP Equation 8.5.26)
          droot(i) = root(i) / rtmass(iplane)
c         -------- initialize root mass yesterday
          rooty(i,iplane) = root(i)
   50   continue
      else
c       ------ ensure non-negative rootmass
        rtmass(iplane) = 0.0
      end if
c
c     ---- fraction of growth permitted by water stress
      watstr(iplane) = 1.0
c     ---- new plant growth today (min_fraction_live_biomass * max_amt_live_biomass)
      slive = rgcmin(plant) * plive(plant,iplane)
c     ---- total live plant material today
      tlive(iplane) = slive
c     ---- previous day's total live plant material
      ptlive(iplane) = tlive(iplane)
c     ---- total above ground plant production for a simulation year
      yield(iplane) = slive
c     ---- daily net primary plant production
      pyield(iplane) = yield(iplane)
c     ---- previous year's date of initial growth
      first(plant) = strrgc(plant)
c     ----  standing dead biomass left after burning or herbicide
      dwood(iplane) = wood(plant) * rmagt(iplane)
c     ---- evergreen leafy component
      xlive(iplane) = slive
c
c     *************************************************************
c     * NOTE: New equations to calculate  canhgt, cancov, rcov,   *
c     *       rilcov, lai for first day of simulation, should be  *
c     *       added to this section when plant growth subroutine  *
c     *       is updated.                                         *
c     *************************************************************
c
c     ---- total above ground dry weight
      vdmt(iplane) = slive + rmagt(iplane)
c     ---- leaf area index
      lai(iplane) = tlive(iplane) * aleaf(plant)
c
c     --- If this is the first OFE, compute SPAI & TPAI for ALL plant groups.
      if (iplane.eq.1) then
        do 60, i = 1, ncrop
c         -------- projected area for shrubs (remains constant throughout simulation)
c         (WEPP Equation 8.5.17)
          if (sdiam(i).gt.0.0) then
            tmpvr4 = sdiam(i) * scoeff(i) * spop(i)
c           ---------- prevent horizontal measure from exceeding 100 m total
            if (tmpvr4.gt.100.0) tmpvr4 = 100.0
            spai(i) = shgt(i) * tmpvr4
          else
            spai(i) = 0.0
          end if
c
c         -------- projected area for trees (remains constant throughout simulation)
c         (WEPP Equation 8.5.17)
          if (tdiam(i).gt.0.0) then
            tmpvr5 = tdiam(i) * tcoeff(i) * tpop(i)
            if (tmpvr5.gt.100.0) tmpvr5 = 100.0
            tpai(i) = thgt(i) * tmpvr5
          else
            tpai(i) = 0.0
          end if
   60   continue
      end if
c
c     ---- projected area of herbaceous plants (changes daily)
c     (WEPP Equation 8.5.17)
      if (gdiam(plant).gt.0.0) then
        tmpvr6 = gdiam(plant) * gcoeff(plant) * gpop(plant)
        if (tmpvr6.gt.100.0) tmpvr6 = 100.0
        gpai = ghgt(plant) * tmpvr6
      else
        gpai = 0.0
      end if
c
c     **********************************************************************
c     * NOTE: We assume NO OVERLAP between the grass, shrubs, and trees.   *
c     *       Since the current purpose of TOTPAI is solely to compute     *
c     *       a weighted average for canopy height, adding the complexity  *
c     *       of considering overlap did not seem worth the cost.          *
c     *       CRM -- 4/7/93.                                               *
c     **********************************************************************
c     ---- total projected area for vegetation
c     (WEPP Equation 8.5.16)
      totpai = spai(plant) + gpai + tpai(plant)
c
c     ---- vertical area required for largest plant type, using 100 meter transect.
c     (WEPP Equation 8.5.18)
      tarea = max(thgt(plant),shgt(plant),ghgt(plant)) * 100.0
c
c     ---- basal density of plants; ie, the fraction of the "vertical rectangle"
c     (re: TAREA) covered by plants
c     (WEPP Equation 8.5.15)
      if (tarea.gt.0.0) then
        basden = totpai / tarea
      else
        basden = 0.0
      end if
c
c     ---- canopy height, a weighted mean.
c     (WEPP Equation 8.5.13)
      if (basden.gt.0.0) then
        canhgt(iplane) = (ghgt(plant)*gpai+shgt(plant)*spai(plant)+
     1      thgt(plant)*tpai(plant)) / totpai
      else
        canhgt(iplane) = 0.0
      end if
c
c     ... calculate rill spacing
c     (WEPP Equation 8.5.19)
c     rspace(iplane)=100.0/(spop(plant)+gpop(plant)+tpop(plant)+1)
c     if(rspace(iplane) .gt. 5.0) rspace(iplane)=5.0
c     if(rspace(iplane) .lt. 0.5) rspace(iplane)=0.5
c
c     ******************************************************************
c     * NOTE: To calculate grass height requires height coefficient    *
c     *       "bbb", which is estimated from the following regression  *
c     *       equations.                                               *
c     ******************************************************************
c     ---- vegetative dry matter at maturity
      vdmmax(iplane) = plive(plant,iplane) + rmagt(iplane) -
     1    dwood(iplane)
c     ----  parameter value for canopy height equation
      if (bbb(plant).ge.0.0) then
        if (bbb(plant).lt.1.5) then
          bbb(plant) = 2.0 * (vdmmax(iplane)**(-1.0)) + 0.3031455e-06
        else if (bbb(plant).lt.2.5) then
c         Original Code:
c         bbb(plant)=4.0*vdmmax(iplane)**(-1.0)+.23155e-05
          bbb(plant) = 4.0 / vdmmax(iplane) + 0.23155e-05
        else
c         Original Code:
c         bbb(plant)=6.0*vdmmax(iplane)**(-1.0)-.1694262e-04
          bbb(plant) = 6.0 / vdmmax(iplane) - 0.1694262e-04
        end if
      end if
c
c     ---- calculate residue cover from residue mass
c     (WEPP Equation 8.5.20)
c     XXX - commented out by Kidwell/Weltz  3/95   dcf
c     tmpvr2 = 1.085 * (1.0+1.069*exp(-12.61*rmogt(1,iplane))) - 0.5838
c     cf(plant) = 10.0 ** (tmpvr2)
c     (WEPP Equation 8.5.20)
c     rescov(iplane) = 1.0 - exp(-cf(plant)*rmogt(1,iplane))
c
c     **Added by Kidwell 5/25/95**
      if (rescof(iplane) .le. 0.0001) then
        if (rescov(iplane).le.0.0 .or. rmogt(1,iplane).le.0.0) then
                cf(plant) = -63.9
                rescov(iplane) = 0.0
                rmogt(1,iplane) = 0.0
        else
                tmpvr2 = log(1.0 - rescov(iplane))
                cf(plant) = tmpvr2/rmogt(1,iplane)
        end if
      else if (rescof(iplane) .gt. 0.0) then
        if(rescov(iplane).le.0.0 .or. rmogt(1,iplane).le.0.0) then
                rescov(iplane) = 0.0
                rmogt(1,iplane) = 0.0
        else
                cf(plant) = rescof(iplane)
        end if
      end if
c
c     ---- calculate canopy cover from plant mass corresponding to 100% cover (COLD)
c     (WEPP Equation 8.5.22)
c     Original Code:
c     ffk(plant) =21.39148-54.90758*cold(plant)+61.11016*cold(plant)
c     1  **2-30.44471*cold(plant)**3+5.561994*cold(plant)**4
c     XXX - commented out next 3 lines (Kidwell/Weltz) 3/95  dcf
c     tmpvr3 = cold(plant) ** 2
c     ffk(plant) = 21.39148 - 54.90758 * cold(plant) + 61.11016 * tmpvr3
c    1    - 30.44471 * cold(plant) * tmpvr3 + 5.561994 * tmpvr3 ** 2
c     -- XXX -- Need a second equation number in the User Doc.
c     (WEPP Equation 8.5.22)
c     XXX Next 2 lines commented out by Kidwell/Weltz  3/95  dcf
c     cancov(iplane) = 1.0 - exp(-ffk(plant)*vdmt(iplane))
c     if (cancov(iplane).lt.0.0) cancov(iplane) = 0.0
c     **Added by Kidwell 5/25/95**
      if (cancof(iplane) .le. 0.0001) then
        if (cancov(iplane).le.0.0 .or. vdmt(iplane).le.0.0) then
                ffk(plant) = -31.5
                cancov(iplane) = 0.0
                vdmt(iplane) = 0.0
        else
                tmpvr3 = log(1.0 - cancov(iplane))
                ffk(plant) = tmpvr3/vdmt(iplane)
        end if
      else if (cancof(iplane) .gt. 0.0) then
        if(cancov(iplane).le.0.0 .or. vdmt(iplane).le.0.0) then
                cancov(iplane) = 0.0
                vdmt(iplane) = 0.0
        else
                ffk(plant) = cancof(iplane)
        end if
      end if

c
c     -- XXX -- Can't find this equation in the User Doc. -- CRM -- 4/07/93.
c     ---- basal plant cover (ground area occupied by stalks, stems, trunks)
c     XXX Next 5 lines commented out by Kidwell/Weltz  3/95   dcf
c     if (wood(plant).gt.0.0) then
c       bcover(iplane) = 0.3354242 * cancov(iplane)
c     else
c       bcover(iplane) = 0.429118 * cancov(iplane)
c     end if
c
c     -- XXX -- Can't find this equation in the User Doc. -- CRM -- 4/07/93.
c     ---- calculate hydraulic roughness cover from basal cover and moss & lichens
c     XXX Next 2 lines changed by Kidwell/Weltz  3/95   dcf
c     rufcov(iplane) = bcover(iplane) + crypto(iplane)
c     if (rufcov(iplane).lt.0.0) rufcov(iplane) = 0.0
      rufcov(iplane) = bascov(iplane)*fbasr(iplane) +
     1                 crycov(iplane)*fcryr(iplane)
c
c     ---- total ground cover (residue, rocks & gravel, stems and moss & lichens)
c     XXX Next line changed by Kidwell/Weltz  3/95   dcf
c     gcover(iplane) = rescov(iplane) + wcf(iplane) + rufcov(iplane)
      gcover(iplane) = rescov(iplane) + rokcov(iplane) + bascov(iplane)
     1                 + crycov(iplane)
      if (gcover(iplane).gt.1.0) gcover(iplane) = 1.0
c
c     -- XXX -- Need an equation number in the User Doc.
c     Following changed by Mary Kidwell 3/95  dcf
c     rilcov(iplane) = rescov(iplane) + wcf(iplane)
c
c     -- XXX -- Need an equation number in the User Doc.
c     inrcov(iplane) = gcover(iplane)
c
c     --- total rill cover
c
      rilcov(iplane) = rescov(iplane)*fresr(iplane) +
     1                 rokcov(iplane)*frokr(iplane) +
     1                 bascov(iplane)*fbasr(iplane) +
     1                 crycov(iplane)*fcryr(iplane)
c
c     ---- total interrill cover
c
      inrcov(iplane) = rescov(iplane)*fresi(iplane) +
     1                 rokcov(iplane)*froki(iplane) +
     1                 bascov(iplane)*fbasr(iplane) +
     1                 crycov(iplane)*fcryi(iplane)
c
c     ---- initial value of dead root mass for day 1 and when management changes
      rtm(1,iplane) = rtmass(iplane)
c
c     -- XXX -- Is DAYMIN needed here?  It is computed in WATBAL before calling
c     RANGE. -- CRM -- 5/14/93.
c     calculate minimum day length for the input latitude
c
c     original code 10/1/91
c     ytn = tan(deglat / 58.09)
c
c     code changed because of problem with Lahey compiler and
c     tangent function when running on AT&T 6300 machine  dcf
c
c     ---- function of latitude - used in day length calculation (EPIC 2.194)
      ytn = (sin(deglat/58.09)) / (cos(deglat/58.09))
      ch = 0.4349 * abs(ytn)
      if (ch.ge.1.0) then
        h = 0.0
      else
        h = acos(ch)
      end if
c     ---- minimum day length
      daymin = 7.72 * h
c
      return
      end

ext: for

max_stars_repo_name: jarad/dep

max_issues_repo_name: jarad/dep

max_forks_repo_name: akrherz/idep


text:
module Mod_nsm_FreeSurfaceStabilization
   use typre
   use Mod_PointerSetter
   use Mod_nsm_BaseElmope
   use Mod_nsm_InterpolateGradients
   use Mod_nsm_ComputeAdvectionVelocity
   implicit none
   private
   public SetPointersFSS, SetPointersFSS_Endite

   logical :: isFirst2Layers
   real(rp), allocatable :: projected_elGrpreAndForce(:,:), projected_elgrvel(:,:,:), projected_GrPreAndForce(:), projected_grvel(:,:)
   real(rp), allocatable :: elgrvel(:,:,:), elgrpre(:,:),elmass(:)
   real(rp), allocatable :: vmass(:)
   integer(ip) :: ielem2

   !For elmopes
   type, extends(PointerSetter) :: SPFreeSurfaceStabilization
contains
      procedure :: SpecificSet => SpecificSetFSS
   end type
   type(SPFreeSurfaceStabilization) :: SetPointersFSS

   !For endite
   type, extends(PointerSetter) :: SPFSS_Endite
contains
      procedure :: SpecificSet => SpecificSetFSS_Endite
   end type

   type(SPFSS_Endite) :: SetPointersFSS_Endite

contains
   subroutine SpecificSetFSS(d)
      implicit none
      class(SPFreeSurfaceStabilization) :: d

      if (a%kfl_StabilizeFreeSurface == 1) then
         call ConcatenateProcedures(ProcHook_Initializations,Allocations)
         call ConcatenateProcedures(ProcHook_Finalizations,Deallocations)
         call PrependProcedure(ProcHook_PreGauss,StabilizationMatrices)
         call PrependProcedure(ProcHook_PreGauss,Gathers)
         call PrependProcedure(ProcHook_PreGauss,CheckIsFirst2Layers)
      endif
   end subroutine

   subroutine SpecificSetFSS_Endite(d)
      implicit none
      class(SPFSS_Endite) :: d

      if (a%kfl_StabilizeFreeSurface == 1) then

         !We will need the interpolated gradients
         call SetPointersInterpolateGradients%Set
         call ConcatenateProcedures(ProcHook_Initializations,Allocations_Endite)
         !We prepend so that we make sure it is done with the full gauss points
         call PrependProcedure(ProcHook_PreGauss,InGaussOperations_Endite)
         call PrependProcedure(ProcHook_PreGauss,ElmatsToZero_Endite)
         call PrependProcedure(ProcHook_PreGauss,CheckIsFirst2Layers)
         call ConcatenateProcedures(ProcHook_AssemblyEndite,Assembly_Endite)
         call ConcatenateProcedures(ProcHook_Finalizations,Finalizations_Endite)
         call ConcatenateProcedures(ProcHook_Finalizations,Deallocations_Endite)
      endif
   end subroutine

   !------------------------------------------------------------------------
   !Free surface stabilization
   subroutine Allocations
      implicit none

      call a%Mesh%GetNdime(ndime)
      call a%Memor%alloc(ndime,ndime,e%mnode,projected_elgrvel,'projected_elgrvel','alloc_freeSurfStabilization')
      call a%Memor%alloc(ndime,e%mnode,projected_elGrpreAndForce,'projected_elGrpreAndForce','alloc_freeSurfStabilization')
      call a%Memor%alloc(ndime,ndime,projected_grvel,'projected_grvel','alloc_freeSurfStabilization')
      call a%Memor%alloc(ndime,projected_GrPreAndForce,'projected_GrPreAndForce','alloc_freeSurfStabilization')
   end subroutine

   subroutine DeAllocations
      implicit none

      call a%Memor%dealloc(ndime,ndime,size(projected_elgrvel,3),projected_elgrvel,'projected_elgrvel','alloc_freeSurfStabilization')
      call a%Memor%dealloc(ndime,size(projected_elGrpreAndForce,2),projected_elGrpreAndForce,'projected_elGrpreAndForce','alloc_freeSurfStabilization')
      call a%Memor%dealloc(ndime,ndime,projected_grvel,'projected_grvel','alloc_freeSurfStabilization')
      call a%Memor%dealloc(ndime,projected_GrPreAndForce,'projected_GrPreAndForce','alloc_freeSurfStabilization')
   end subroutine

   subroutine CheckIsFirst2Layers
      implicit none
      integer(ip) :: inode,ipoin,poinStatus

      isFirst2Layers = .false.
      do inode = 1,e%pnode
         ipoin = e%lnods(inode)
         call a%CutMesh%GetPointType(ipoin,poinStatus)
         if (poinStatus == 1) then
            !This means the element belongs to either the cut layer or the adjacent one inside
            !I want to add stabilization in these layers
            isFirst2Layers = .true.
            exit
         endif
      enddo
   end subroutine

   subroutine Gathers
      implicit none

      if (isFirst2Layers) then
         call e%gather(ndime*ndime,projected_elgrvel,a%StabilizeFreeSurface_VelocityGradient)
         call e%gather(ndime,projected_elGrpreAndForce,a%StabilizeFreeSurface_PressureGradientAndForce)
      endif
   end subroutine

   subroutine StabilizationMatrices
      implicit none
      real(rp) :: velgraddiff(ndime,ndime), pregraddiff(ndime)
      integer(ip) :: idime,inode,jnode
      real(rp)  :: ghosttimom
      logical :: dopressure

      if (isFirst2Layers) then
         dopressure = .false.
         if (size(elmat,1) > e%ndime) dopressure = .true.
         gauss_points : do igaus=1,e%pgaus
            e%igaus = igaus
            dvol = e%weigp(e%igaus)*e%detjm

            !Interpolate
            call InterpolateGpVelocities
            !Hook
            call ProcHook_Interpolates
            !Physical Properties
            call a%GetPhysicalParameters(imat,acden,acvis)
            !Hook
            call ProcHook_PhysicalProp

            !Advection velocity
            call ProcPointer_ComputeAdvectionVelocity
            if (a%kfl_advec == 1) gpadv = gpvel(:,1)
            if (a%kfl_advco == 1) gpadv(1:e%ndime) = a%advco(1:e%ndime)

            !Advection velocity norm
            call vecnor(gpadv,e%ndime,gpvno,2)
            !Compute a·grad(V)
            call ComputeAGradV(e,gpadv,AGradV)

            !Compute the stability parameters
            call ProcHook_ComputeTaus
            call e%interpg(ndime*ndime,projected_elgrvel,projected_grvel)
            call e%interpg(ndime,projected_elGrpreAndForce,projected_GrPreAndForce)

            !Compute Elext
            elext=0.0_rp
            !Compute vector of external forces
            call ProcPointer_ExternalForces

            !Stabilize the velocity array
            !first the right hand side
            do inode = 1,e%pnode
               !Stabilize velocity array
               do idime = 1,ndime
                  elrhs(idime,inode) = elrhs(idime,inode) + a%StabilizeFreeSurface_Param*tidiv*dot_product(e%cartd(:,inode),projected_grvel(:,idime))*e%weigp(e%igaus)*e%detjm
               enddo
            enddo
            !Second the left hand side
            do inode = 1,e%pnode
               do jnode = 1,e%pnode
                  !Stabilize the velocity array
                  do idime = 1,ndime
                     elmat(idime,inode,idime,jnode) = elmat(idime,inode,idime,jnode) + a%StabilizeFreeSurface_Param*tidiv*dot_product(e%cartd(:,inode),e%cartd(:,jnode))*e%weigp(e%igaus)*e%detjm
                  enddo
               enddo
            enddo

            !Stabilize pressure array
            if (dopressure) then
               do inode = 1,e%pnode
                   elrhs(ndime+1,inode) = elrhs(ndime+1,inode) + a%StabilizeFreeSurface_Param*timom*dot_product(e%cartd(:,inode),(projected_GrPreAndForce(:)+elext(1:e%ndime)))*e%weigp(e%igaus)*e%detjm
               enddo
               !Second the left hand side
               do inode = 1,e%pnode
                  do jnode = 1,e%pnode
                     elmat(ndime+1,inode,ndime+1,jnode) = elmat(ndime+1,inode,ndime+1,jnode) + a%StabilizeFreeSurface_Param*timom*dot_product(e%cartd(:,inode),e%cartd(:,jnode))*e%weigp(e%igaus)*e%detjm
                  enddo
               enddo
            endif
         enddo gauss_points
      endif

   end subroutine

   !------------------------------------------------------------------
   !For endite
   subroutine Allocations_Endite
      implicit none
      integer(ip) :: npoin

      a%StabilizeFreeSurface_VelocityGradient = 0.0_rp
      a%StabilizeFreeSurface_PressureGradientAndForce = 0.0_rp
      call a%Mesh%GetNpoin(npoin)
      call a%Mesh%GetNdime(ndime)
      call a%Memor%alloc(ndime,ndime,e%mnode,elgrvel,'elgrvel','Mod_nsm_FreeSurfaceStabilization')
      call a%Memor%alloc(ndime,e%mnode,elgrpre,'elgrpre','Mod_nsm_FreeSurfaceStabilization')
      call a%Memor%alloc(e%mnode,elmass,'elmass','Mod_nsm_FreeSurfaceStabilization')
      call a%Memor%alloc(npoin,vmass,'auxvmass','Mod_nsm_FreeSurfaceStabilization')

   end subroutine

   subroutine Deallocations_Endite
      implicit none
      integer(ip) :: npoin

      call a%Mesh%GetNpoin(npoin)
      call a%Memor%dealloc(ndime,ndime,e%mnode,elgrvel,'elgrvel','Mod_nsm_FreeSurfaceStabilization')
      call a%Memor%dealloc(ndime,e%mnode,elgrpre,'elgrpre','Mod_nsm_FreeSurfaceStabilization')
      call a%Memor%dealloc(e%mnode,elmass,'elmass','Mod_nsm_FreeSurfaceStabilization')
      call a%Memor%dealloc(npoin,vmass,'auxvmass','Mod_nsm_FreeSurfaceStabilization')

   end subroutine

   subroutine ElmatsToZero_Endite
      implicit none

      !isFirst2Layers = .true.
      ielem2 = ielem
      if (isFirst2Layers) then
         elgrvel = 0.0_rp
         elgrpre = 0.0_rp
         elmass = 0.0_rp
      endif

   end subroutine

   subroutine InGaussOperations_Endite
      implicit none
      integer(ip) :: inode

      if (isFirst2Layers) then

         gauss_points : do igaus=1,e%pgaus
            e%igaus = igaus
            dvol = e%weigp(e%igaus)*e%detjm

            !Interpolate
            call InterpolateGpVelocities
            !Hook
            call ProcHook_Interpolates
            !Compute Elext, Temporal Derivatives
            elext=0.0_rp
            !Compute vector of external forces
            call ProcPointer_ExternalForces
            do inode = 1,e%pnode
               elgrvel(:,:,inode) = elgrvel(:,:,inode) + e%shape(inode,e%igaus)*grvel*dvol
               elgrpre(:,inode)   = elgrpre(:,inode)   + e%shape(inode,e%igaus)*(grpre(1,:)- elext(1:e%ndime))*dvol
               elmass(inode) = elmass(inode) + e%shape(inode,e%igaus)*dvol
            enddo
         enddo gauss_points
      endif

   end subroutine

   subroutine Assembly_Endite
      implicit none
      if (isFirst2Layers) then
         call a%Mesh%AssemblyToArray(e,ndime*ndime,elgrvel,a%StabilizeFreeSurface_VelocityGradient)
         call a%Mesh%AssemblyToArray(e,ndime,elgrpre,a%StabilizeFreeSurface_PressureGradientAndForce)
         call a%Mesh%AssemblyToArray(e,1_ip,elmass,vmass)
      endif

   end subroutine

   subroutine Finalizations_Endite
      use Mod_int2str
      implicit none
      interface
         subroutine AfterProjectOperations(a,ndofn,array)
            use typre
            use Mod_Memor
            use Mod_Mesh
            implicit none
            class(FemMesh) :: a
            integer(ip) :: ndofn
            real(rp)    :: array(:,:)
         end subroutine
      end interface

      integer(ip) :: ipoin,npoin,idime
      real(rp), pointer :: PointerArray(:,:) => NULL()
      character(150) :: auxname

      auxname = 'GraPress'//trim(adjustl(int2str(a%itera)))
      call a%Mesh%GetNpoin(npoin)
      do ipoin = 1,npoin
         if (vmass(ipoin) /= 0.0_rp) then
            a%StabilizeFreeSurface_VelocityGradient(:,:,ipoin) = a%StabilizeFreeSurface_VelocityGradient(:,:,ipoin)/vmass(ipoin)
            a%StabilizeFreeSurface_PressureGradientAndForce(:,ipoin) = a%StabilizeFreeSurface_PressureGradientAndForce(:,ipoin)/vmass(ipoin)
         endif
      enddo

      PointerArray(1:ndime*ndime,1:npoin) => a%StabilizeFreeSurface_VelocityGradient
      call AfterProjectOperations(a%Mesh,ndime*ndime,PointerArray)
      call AfterProjectOperations(a%Mesh,ndime      ,a%StabilizeFreeSurface_PressureGradientAndForce)

   end subroutine

end module


! Copyright 2005-2019 ECMWF.
!
! This software is licensed under the terms of the Apache Licence Version 2.0
! which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
!
! In applying this licence, ECMWF does not waive the privileges and immunities granted to it by
! virtue of its status as an intergovernmental organisation nor does it submit to any jurisdiction.
!
!
!
program new_from_file
use eccodes
  implicit none
  integer           :: ifile
  integer           :: iret
  integer           :: count1=0

  ! Message identifier.
  integer            :: igrib

  ifile=5

  call codes_open_file(ifile,'../../data/collection.grib1','r')

  ! Loop on all the messages in a file.
  call codes_grib_new_from_file(ifile,igrib, iret)

  do while (iret==CODES_SUCCESS)
    count1=count1+1
    print *, "===== Message #",count1
    call codes_grib_new_from_file(ifile,igrib, iret)

  end do
  if (iret /= CODES_END_OF_FILE) then
    call codes_check(iret,'new_from_file','')
  endif

  call codes_close_file(ifile)

end program


text:
*     ***************************************
*     *										*
*     *      fcoord_to_real					*
*     *										*
*     ***************************************
*
*     This routine converts from fractional coordinates to real
*     coordinates.
*
*     Entry - nion
*             a(3,3): lattice vectors
*             ion: fractional coordinates
*    Exit -
*             ion: real coordinates
*

      subroutine fcoord_to_real(nion,ion)
      implicit none
      integer nion
      real*8 ion(3,*)

c     **** local variables ****
      integer i,j
      real*8 tion(3)
      real*8 a(3,3)

*     *** external functions ****
      real*8   lattice_unita
      external lattice_unita

      do j=1,3
      do i=1,3
        a(i,j) = lattice_unita(i,j)
      end do
      end do

      do i=1,nion
         tion(1) = ion(1,i)
         tion(2) = ion(2,i)
         tion(3) = ion(3,i)

         ion(1,i) = a(1,1)*tion(1)
     >            + a(1,2)*tion(2)
     >            + a(1,3)*tion(3)
         ion(2,i) = a(2,1)*tion(1)
     >            + a(2,2)*tion(2)
     >            + a(2,3)*tion(3)
         ion(3,i) = a(3,1)*tion(1)
     >            + a(3,2)*tion(2)
     >            + a(3,3)*tion(3)
      end do

      return
      end

*     ***************************************
*     *										*
*     *      fcoord_to_frac 				*
*     *										*
*     ***************************************
*
*     This routine converts from real coordinates to fractional
*     coordinates.
*
*     Entry - nion
*             a(3,3): lattice vectors
*             ion: real coordinates
*     Exit -
*             ion: fractional coordinates

      subroutine fcoord_to_frac(nion,ion)
      implicit none
      integer nion
      real*8 ion(3,*)

*     **** local variables ****
      integer i,j
      real*8 a(3,3),b(3,3),volume
      real*8 tion(3)

*     *** external functions ****
      real*8   lattice_unita
      external lattice_unita

      do j=1,3
      do i=1,3
        a(i,j) = lattice_unita(i,j)
      end do
      end do

      b(1,1) = a(2,2)*a(3,3) - a(3,2)*a(2,3)
      b(2,1) = a(3,2)*a(1,3) - a(1,2)*a(3,3)
      b(3,1) = a(1,2)*a(2,3) - a(2,2)*a(1,3)
      b(1,2) = a(2,3)*a(3,1) - a(3,3)*a(2,1)
      b(2,2) = a(3,3)*a(1,1) - a(1,3)*a(3,1)
      b(3,2) = a(1,3)*a(2,1) - a(2,3)*a(1,1)
      b(1,3) = a(2,1)*a(3,2) - a(3,1)*a(2,2)
      b(2,3) = a(3,1)*a(1,2) - a(1,1)*a(3,2)
      b(3,3) = a(1,1)*a(2,2) - a(2,1)*a(1,2)
      volume = a(1,1)*b(1,1)
     >       + a(2,1)*b(2,1)
     >       + a(3,1)*b(3,1)

      volume = 1.0d0/volume
      call dscal(9,volume,b,1)

      do i=1,nion
         tion(1) = ion(1,i)
         tion(2) = ion(2,i)
         tion(3) = ion(3,i)
         ion(1,i) = b(1,1)*tion(1)
     >            + b(2,1)*tion(2)
     >            + b(3,1)*tion(3)

         ion(2,i) = b(1,2)*tion(1)
     >            + b(2,2)*tion(2)
     >            + b(3,2)*tion(3)

         ion(3,i) = b(1,3)*tion(1)
     >            + b(2,3)*tion(2)
     >            + b(3,3)*tion(3)

      end do

      return
      end

c $Id$

max_stars_repo_name: dinisAbranches/nwchem

max_issues_repo_name: dinisAbranches/nwchem

max_forks_repo_name: dinisAbranches/nwchem


text:
Teknolust is a fulllength, feature film released in 2002 that was written and directed by Lynn HershmanLeeson who was faculty in the TechnoCultural Studies Department, but was in the Art Department at the time of production.

Basic plot: Nerdy science woman (Dr. Rosetta Stone) invents computer programs that are women who can exist outside of the computer. They are powered by semen, so they fellate men wearing condoms and make semen tea to survive. Fellated men become stricken with virus. Women become autonomous, chaos ensues.

Cool things about movie: lead actress Tilda Swinton (http://www.imdb.com/name/nm0842770/ IMBD) plays four roles and is totally hot as three bots and manages to be nerdy and attractive as science lady. Features EU II/Kemper Hall a few times for a total of a good 10 seconds (the signature floor tiles, a room in the basement, and a transition scene of people going up the spiral staircase) as well as a lookalike of the little green electric car you might have seen around town.

The movie was shot in http://www.imdb.com/title/tt0270688/technical HD, and is breathtakingly colorful and fun. In February 2003, it was screened at the Mondavi Center; Lynn HershmanLeeson and Tilda Swinton were part of a panel discussion prior to the screening.

More at http://www.imdb.com/title/tt0270688/ IMDB and there was a website Xhttp://www.agentruby.com/indexflash.html websiteX run by Hershman (required flash) but it now seems to either be porn or down (Title line saidLive nude webcams online.agentruby.com).

The DVD is out, so you can buy it or rent it (I came across it at a Hollywood Video in Palo Alto, CA) Users/PaulIvanov

As of the end of Winter 2005, Professor Hershman is no longer a member of the UCD Faculty.  Due to complaints by several students about the way her Technocultural Studies 100 class was being run, she was hassled by the administration and eventually quit out of frustration with the university bureaucracy.  It was sad to see how it went down last quarter, being in her last course after all. Users/StephenHo

    What were the complaints? Users/JosephBleckman  20051229

I saw it last night, it was awesome, the best movie I have ever seen. Users/NickSchmalenberger 20051229

20081003 00:04:58 nbsp This movie is really amazingly awesome. Users/NickSchmalenberger


text:
!> Tuning of HV data
subroutine thermopack_getHVparam(i,j,alpha_ij,alpha_ji,aGE_ij,aGE_ji,bGE_ij,bGE_ji,cGE_ij,cGE_ji)
  use thermopack_var
  use cubic_eos, only: cb_eos
  implicit none
  integer, intent(in) :: i,j
  real, intent(out) :: alpha_ij,alpha_ji,aGE_ij,aGE_ji,bGE_ij,bGE_ji,cGE_ij,cGE_ji
  class(base_eos_param), pointer :: act_eos_ptr
  act_eos_ptr => get_active_eos()
  !
  select type(p_eos => act_eos_ptr)
  class is (cb_eos)
    if (.not. allocated(p_eos%mixGE%alpha)) then
      call stoperror('p_eos%mixGEalpha not allocated')
    endif
    if (.not. allocated(p_eos%mixGE%aGE)) then
      call stoperror('p_eos%mixGE%aGE not allocated')
    endif
    if (.not. allocated(p_eos%mixGE%bGE)) then
      call stoperror('p_eos%mixGE%bGE not allocated')
    endif
    if (.not. allocated(p_eos%mixGE%cGE)) then
      call stoperror('p_eos%mixGE%cGE not allocated')
    endif
    alpha_ij = p_eos%mixGE%alpha(i,j)
    alpha_ji = p_eos%mixGE%alpha(j,i)
    aGE_ij = p_eos%mixGE%aGE(i,j)
    aGE_ji = p_eos%mixGE%aGE(j,i)
    bGE_ij = p_eos%mixGE%bGE(i,j)
    bGE_ji = p_eos%mixGE%bGE(j,i)
    cGE_ij = p_eos%mixGE%cGE(i,j)
    cGE_ji = p_eos%mixGE%cGE(j,i)
  class default
    print *,"thermopack_getHVparam: Wrong model - no HV parameters"
  end select
end subroutine thermopack_getHVparam

!> Tuning of HV data
subroutine thermopack_setHVparam(i,j,alpha_ij,alpha_ji,aGE_ij,aGE_ji,bGE_ij,bGE_ji,cGE_ij,cGE_ji)
  use thermopack_var
  use cubic_eos, only: cb_eos
  implicit none
  integer, intent(in) :: i,j
  real, intent(in) :: alpha_ij,alpha_ji,aGE_ij,aGE_ji,bGE_ij,bGE_ji,cGE_ij,cGE_ji
  !
  class(base_eos_param), pointer :: act_eos_ptr
  act_eos_ptr => get_active_eos()
  !
  select type(p_eos => act_eos_ptr)
  class is (cb_eos)
    if (.not. allocated(p_eos%mixGE%alpha)) then
      call stoperror('p_eos%mixGE%alpha not allocated')
    endif
    if (.not. allocated(p_eos%mixGE%aGE)) then
      call stoperror('p_eos%mixGE%aGE not allocated')
    endif
    if (.not. allocated(p_eos%mixGE%bGE)) then
      call stoperror('p_eos%mixGE%bGE not allocated')
    endif
    if (.not. allocated(p_eos%mixGE%cGE)) then
      call stoperror('p_eos%mixGE%cGE not allocated')
    endif
    p_eos%mixGE%alpha(i,j) = alpha_ij
    p_eos%mixGE%alpha(j,i) = alpha_ji
    p_eos%mixGE%aGE(i,j) = aGE_ij
    p_eos%mixGE%aGE(j,i) = aGE_ji
    p_eos%mixGE%bGE(i,j) = bGE_ij
    p_eos%mixGE%bGE(j,i) = bGE_ji
    p_eos%mixGE%cGE(i,j) = cGE_ij
    p_eos%mixGE%cGE(j,i) = cGE_ji
    p_eos%mixGE%correlation(i,j) = 1
    p_eos%mixGE%correlation(j,i) = 1
  class default
    print *,"thermopack_setHVparam: Wrong model - no HV parameters"
  end select
end subroutine thermopack_setHVparam

!> Tuning of WS data
subroutine thermopack_setWSparam(i,j, alpha_ij,alpha_ji, k_ij,k_ji, tau_ij,tau_ji)
  use thermopack_var
  use cubic_eos, only: cb_eos, fraction
  implicit none
  integer, intent(in) :: i,j
  real, intent(in) :: alpha_ij,alpha_ji, k_ij,k_ji, tau_ij,tau_ji
  !
  type(fraction):: frac, fracji
  class(base_eos_param), pointer :: act_eos_ptr
  act_eos_ptr => get_active_eos()
  !
  select type(p_eos => act_eos_ptr)
  class is (cb_eos)
     if (.not. allocated(p_eos%mixWS%f_kij)) then
        call stoperror('p_eos%mixWS%f_kij not allocated')
     endif
     if (.not. allocated(p_eos%mixWS%f_tauij)) then
        call stoperror('p_eos%mixWS%f_tauij not allocated')
     endif
     if (.not. allocated(p_eos%mixWS%alphaij)) then
        call stoperror('p_eos%mixWS%alphaij not allocated')
     endif

     frac%pNum = 0.0
     frac%pDen = 0.0
     frac%pDen(1) = 1.0

     ! Set ij parameters
     p_eos%mixWS%alphaij(i,j) = alpha_ij

     frac%pnum(1) = k_ij
     p_eos%mixWS%f_kij(i,j) = frac

     frac%pnum(1) = tau_ij
     p_eos%mixWS%f_tauij(i,j) = frac

     ! Set ji parameters
     p_eos%mixWS%alphaij(j,i) = alpha_ji

     frac%pnum(1) = k_ji
     p_eos%mixWS%f_kij(j,i) = frac

     frac%pnum(1) = tau_ji
     p_eos%mixWS%f_tauij(j,i) = frac

  class default
     print *,"thermopack_setWSparam: Wrong model - no WS parameters"
  end select
end subroutine thermopack_setWSparam

! Tuning of WS data
subroutine thermopack_getWSparam(i,j,alpha_ij,alpha_ji, k_ij,k_ji, tau_ij,tau_ji)
  use thermopack_var
  use cubic_eos, only: cb_eos
  implicit none
  integer, intent(in) :: i,j
  real, intent(out) :: alpha_ij,alpha_ji, k_ij,k_ji, tau_ij,tau_ji
  !
  class(base_eos_param), pointer :: act_eos_ptr
  act_eos_ptr => get_active_eos()
  !
  select type(p_eos => act_eos_ptr)
  class is (cb_eos)
     if (.not. allocated(p_eos%mixWS%f_kij)) then
        call stoperror('p_eos%mixWS%f_kij not allocated')
     endif
     if (.not. allocated(p_eos%mixWS%f_tauij)) then
        call stoperror('p_eos%mixWS%f_tauij not allocated')
     endif
     if (.not. allocated(p_eos%mixWS%alphaij)) then
        call stoperror('p_eos%mixWS%alphaij not allocated')
     endif

     k_ij = p_eos%mixWS%f_kij(i,j)%pNum(1)
     k_ji = p_eos%mixWS%f_kij(j,i)%pNum(1)

     tau_ij = p_eos%mixWS%f_tauij(i,j)%pNum(1)
     tau_ji = p_eos%mixWS%f_tauij(j,i)%pNum(1)

     alpha_ij = p_eos%mixWS%alphaij(i,j)
     alpha_ji = p_eos%mixWS%alphaij(j,i)
  class default
     print *,"thermopack_getWSparam: Wrong model - no WS parameters"
  end select
end subroutine thermopack_getWSparam

!> Tuning of TWU alpha correlation
subroutine thermopack_getTWUparam(i,c_1,c_2,c_3)
  use thermopack_var
  use cubic_eos, only: cb_eos
  implicit none
  integer, intent(in) :: i
  real, intent(out) :: c_1,c_2,c_3
  class(base_eos_param), pointer :: act_eos_ptr
  !
  act_eos_ptr => get_active_eos()
  !
  select type(p_eos => act_eos_ptr)
  class is (cb_eos)
    c_1 = p_eos%single(i)%alphaParams(1)
    c_2 = p_eos%single(i)%alphaParams(2)
    c_3 = p_eos%single(i)%alphaParams(3)
  class default
    print *,"thermopack_getTWUparam: Wrong model - no TWU parameters"
  end select
end subroutine thermopack_getTWUparam

!> Tuning of TWU alpha correlation. If eoslibinit wasn't initilized with
!> TWU, TWU becomes the new alpha correlation.
subroutine thermopack_setTWUparam(i,c_1,c_2,c_3)
  use cbAlpha, only: setSingleAlphaCorr
  use thermopack_var
  use cubic_eos, only: cb_eos, cbAlphaTwuIdx
  implicit none
  integer, intent(in) :: i
  real, intent(in) :: c_1,c_2,c_3
  class(base_eos_param), pointer :: act_eos_ptr
  !
  act_eos_ptr => get_active_eos()
  !
  select type(p_eos => act_eos_ptr)
  class is (cb_eos)
    call setSingleAlphaCorr(i, p_eos, alphaIdx=cbAlphaTwuIdx, alphaParams=(/c_1,c_2,c_3/))
  class default
    print *,"thermopack_setTWUparam: Wrong model - no TWU parameters"
  end select
end subroutine thermopack_setTWUparam

!> Retrieve TWU parameters
subroutine thermopack_getMCparam(i,c_1,c_2,c_3)
  use thermopack_var
  use cubic_eos, only: cb_eos
  implicit none
  integer, intent(in) :: i
  real, intent(out) :: c_1,c_2,c_3
  !
  class(base_eos_param), pointer :: act_eos_ptr
  !
  act_eos_ptr => get_active_eos()
  !
  select type(p_eos => act_eos_ptr)
  class is (cb_eos)
    c_1 = p_eos%single(i)%alphaParams(1)
    c_2 = p_eos%single(i)%alphaParams(2)
    c_3 = p_eos%single(i)%alphaParams(3)
  class default
    print *,"thermopack_getMCparam: Wrong model - no MC parameters"
  end select
end subroutine thermopack_getMCparam

!> Tuning of MC alpha correlation. If eoslibinit wasn't initilized with
!> MC, MC becomes the new alpha correlation.
subroutine thermopack_setMCparam(i,c_1,c_2,c_3)
  use cbAlpha, only: setSingleAlphaCorr
  use thermopack_var
  use cubic_eos, only: cb_eos, cbAlphaMCIdx
  implicit none
  integer, intent(in) :: i
  real, intent(in) :: c_1,c_2,c_3
  class(base_eos_param), pointer :: act_eos_ptr
  !
  act_eos_ptr => get_active_eos()
  !
  select type(p_eos => act_eos_ptr)
  class is (cb_eos)
    call setSingleAlphaCorr(i, p_eos, alphaIdx=cbAlphaMCIdx, alphaParams=(/c_1,c_2,c_3/))
  class default
    print *,"thermopack_setMCparam: Wrong model - no MC parameters"
  end select
end subroutine thermopack_setMCparam

subroutine thermopack_setClassicFitparam(i,c_1)
  use cbAlpha, only: setSingleAlphaCorr
  use thermopack_var
  use cubic_eos, only: cb_eos
  implicit none
  integer, intent(in) :: i
  real, intent(in) :: c_1
  ! Locals
  integer :: alphaIdx
  class(base_eos_param), pointer :: act_eos_ptr
  !
  act_eos_ptr => get_active_eos()
  !
  select type(p_eos => act_eos_ptr)
  class is (cb_eos)
    alphaIdx = p_eos%single(i)%alphaMethod
    call setSingleAlphaCorr(i, p_eos, alphaIdx=alphaIdx, alphaParams=(/c_1/))
  class default
    print *,"thermopack_setClassicFitparam: Wrong model - no ClassicFit parameters"
  end select
end subroutine thermopack_setClassicFitparam

!> Tuning of vdW kij interaction parameters
subroutine thermopack_getkij(i,j,kij)
  use thermopack_var
  use cubic_eos, only: cb_eos
  implicit none
  integer, intent(in) :: i,j
  real, intent(out) :: kij
  class(base_eos_param), pointer :: act_eos_ptr
  !
  act_eos_ptr => get_active_eos()
  !
  select type(p_eos => act_eos_ptr)
  class is (cb_eos)
    kij = p_eos%kij(i,j)
  class default
    print *,"thermopack_getkij: Wrong model - not cubic"
  end select
end subroutine thermopack_getkij

!> Tuning of vdW kij interaction parameters
subroutine thermopack_setkij(i,j,kij)
  use thermopack_var
  use cubic_eos, only: cb_eos
  implicit none
  integer, intent(in) :: i,j
  real, intent(in) :: kij
  class(base_eos_param), pointer :: act_eos_ptr
  !
  act_eos_ptr => get_active_eos()
  !
  select type(p_eos => act_eos_ptr)
  class is (cb_eos)
    p_eos%kij(i,j) = kij
  class default
    print *,"thermopack_setkij: Wrong model - not cubic"
  end select
end subroutine thermopack_setkij

!> Tuning of kij interaction parameters
subroutine thermopack_setkijandji(i,j,kij)
  use thermopack_var
  use eos_parameters, only: single_eos
  use cubic_eos, only: cb_eos, cpa_eos
  use saftvrmie_containers, only: set_saftvrmie_eps_kij, saftvrmie_eos
  use saft_interface, only: pc_saft_set_kij
  use pc_saft_nonassoc, only: PCSAFT_eos
  implicit none
  integer, intent(in) :: i,j
  real, intent(in) :: kij
  type(thermo_model), pointer :: act_mod_ptr
  class(base_eos_param), pointer :: act_eos_ptr
  !
  act_mod_ptr => get_active_thermo_model()
  act_eos_ptr => get_active_eos()
  select type (p_eos => act_eos_ptr)
  class is (cb_eos)
    p_eos%kij(i,j) = kij
    p_eos%kij(j,i) = kij
  type is(cpa_eos)
    call stoperror("Not able to set binary kij for CPA eos.")
  type is(PCSAFT_eos)
    call pc_saft_set_kij(i,j,kij)
  type is(single_eos)
    call stoperror("Not able to set binary kij for single comp. eos.")
  type is(saftvrmie_eos)
    call set_saftvrmie_eps_kij(i,j,kij)
  class default
    call stoperror("Not able to set binary kij for selected model")
  end select

end subroutine thermopack_setkijandji

!> Tuning of lij interaction parameters
subroutine thermopack_setlijandji(i,j,lij)
  use thermopack_var
  use eos_parameters, only: single_eos
  use cubic_eos, only: cb_eos, cpa_eos
  use saftvrmie_containers, only: set_saftvrmie_sigma_lij, saftvrmie_eos
  use saft_interface, only: pc_saft_set_kij
  use pc_saft_nonassoc, only: PCSAFT_eos
  implicit none
  integer, intent(in) :: i,j
  real, intent(in) :: lij
  type(thermo_model), pointer :: act_mod_ptr
  class(base_eos_param), pointer :: act_eos_ptr
  !
  act_mod_ptr => get_active_thermo_model()
  act_eos_ptr => get_active_eos()
  select type (p_eos => act_eos_ptr)
  class is (cb_eos)
    call stoperror("Not able to set binary lij for Cubic eos.")
  type is(cpa_eos)
    call stoperror("Not able to set binary lij for CPA eos.")
  type is(PCSAFT_eos)
    call stoperror("Not able to set binary lij for PC-SAFT eos.")
  type is(single_eos)
    call stoperror("Not able to set binary lij for single comp. eos.")
  type is(saftvrmie_eos)
    call set_saftvrmie_sigma_lij(i,j,lij)
  class default
    call stoperror("Not able to set binary lij for selected model")
  end select
end subroutine thermopack_setlijandji

!> Get volume translation parameters
subroutine thermopack_get_volume_shift_parameters(i,ciA,ciB,ciC,ci_type)
  use thermopack_var
  implicit none
  integer, intent(in) :: i
  real, intent(out) :: ciA,ciB,ciC
  integer, intent(out) :: ci_type
  ! Locals
  type(thermo_model), pointer :: act_mod_ptr
  !
  act_mod_ptr => get_active_thermo_model()
  ciA = act_mod_ptr%comps(i)%p_comp%cid%ciA
  ciB = act_mod_ptr%comps(i)%p_comp%cid%ciB
  ciC = act_mod_ptr%comps(i)%p_comp%cid%ciC
  ci_type = act_mod_ptr%comps(i)%p_comp%cid%c_type
end subroutine thermopack_get_volume_shift_parameters

!> Tuning of volume translation parameters
subroutine thermopack_set_volume_shift_parameters(i,ciA,ciB,ciC,ci_type)
  use thermopack_var
  implicit none
  integer, intent(in) :: i
  real, intent(in) :: ciA,ciB,ciC
  integer, intent(in) :: ci_type
  ! Locals
  type(thermo_model), pointer :: act_mod_ptr
  !
  act_mod_ptr => get_active_thermo_model()
  act_mod_ptr%comps(i)%p_comp%cid%ciA = ciA
  act_mod_ptr%comps(i)%p_comp%cid%ciB = ciB
  act_mod_ptr%comps(i)%p_comp%cid%ciC = ciC
  act_mod_ptr%comps(i)%p_comp%cid%c_type = ci_type
end subroutine thermopack_set_volume_shift_parameters

max_stars_repo_name: SINTEF/Thermopack

max_issues_repo_name: SINTEF/Thermopack

max_forks_repo_name: SINTEF/Thermopack


!    Copyright (C) 2012 The SPEED FOUNDATION
!    Author: Ilario Mazzieri
!
!    This file is part of SPEED.
!
!    SPEED is free software; you can redistribute it and/or modify it
!    under the terms of the GNU Affero General Public License as
!    published by the Free Software Foundation, either version 3 of the
!    License, or (at your option) any later version.
!
!    SPEED is distributed in the hope that it will be useful, but
!    WITHOUT ANY WARRANTY; without even the implied warranty of
!    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
!    Affero General Public License for more details.
!
!    You should have received a copy of the GNU Affero General Public License
!    along with SPEED.  If not, see <http://www.gnu.org/licenses/>.

!> @brief Makes a pointer for spectral connectivity vector.
!! @author Ilario Mazzieri
!> @date September, 2013
!> @version 1.0
!> @param[in] nb_node  GLL nodes
!> @param[in] spx_con spectral connectivity vector
!> @param[in] node_wgt node_wgt(i) contains the multiplicity of the GLL node i
!!                  (e.g. if the GLL node i is shared by 4 elements node_wgt(i) = 4
!> @param[in] nnz  number of GLL nodes + number of GLL nodes including repetitions + 1
!> @param[out] node_pointer as explained in MAKE_GRID_NODES.f90

      subroutine MAKE_SPX_NODES(nb_node,spx_con_nnz,spx_con,node_wgt,nnz,node_pointer)

      implicit none

      integer*4 :: nb_node,spx_con_nnz,nnz
      integer*4 :: i,j,ie,nb_elem,nn3

      integer*4, dimension(:), allocatable :: i4count
      integer*4, dimension(0:spx_con_nnz) :: spx_con
      integer*4, dimension(nb_node) :: node_wgt
      integer*4, dimension(0:nnz) :: node_pointer

      allocate(i4count(nb_node))

      node_pointer(0) = nb_node + 1
      do j = 1, nb_node
         node_pointer(j) = node_pointer(j -1) + node_wgt(j)
      enddo

      do j = 1, nb_node
         i4count(j) = node_pointer(j -1)
      enddo

      nb_elem = spx_con(0) - 1
      do ie = 1, nb_elem

         nn3 = spx_con(ie) - spx_con(ie -1) -1
         do i = 1, nn3

            j = spx_con(spx_con(ie -1) +i)
            node_pointer(i4count(j)) = ie
            i4count(j) = i4count(j) + 1
         enddo
      enddo

      deallocate(i4count)

      return

      end subroutine MAKE_SPX_NODES


text:
c
c-----------------------------------------------------------------------
c
c
      subroutine p18(xupw,yupw,n)
      implicit double precision(a-h,o-z)
      dimension xupw(n),yupw(n)
      common/p18c/ delta,pmu0,theta1
c
c     # find the upper wall of the rarefaction channel for problem 18
c
      pi2 = datan(1.d0)*2.d0
      delta = dsqrt(2.4d0/0.4d0)
      amach = 1.30255d0
      q0 = 1.d0
      a0 = q0/amach
      p0 = 9.04545d0
      rho0 = 1.4d0*p0/a0**2
      dmu0 = dasin(1.d0/amach)
      pmu0 = delta*datan(delta*dtan(dmu0)) - dmu0
      h = 1.7d0/dfloat(n)
      xw0 = .1d0
      yw0 = .302d0
      width0 = .2d0
      theta0 = 0.d0
      x0 = xw0 + width0/dtan(theta0+dmu0)
      y0 = yw0 + width0
      xupw(1) = 0.d0
      yupw(1) = y0
      xupw(2) = x0
      yupw(2) = y0
c
      do 10 i=1,n-2
	 xw = xw0 + i*h
	 yw = p18yw(xw)
	 theta = p18th(xw)
	 dmu = p18mu(theta)
	 aratio = dsqrt((1.d0+0.4d0/(2.d0*dsin(dmu0)**2)) /
     &        (1.d0+0.4d0/(2.d0*dsin(dmu)**2)))
	 q = a0*aratio/dsin(dmu)
	 aa = a0*aratio
	 p = p0*((1.d0+0.4d0/(2.d0*dsin(dmu0)**2)) /
     &        (1.d0+0.4d0/(2.d0*dsin(dmu)**2)))**(1.4d0/0.4d0)
	 rho = 1.4d0*p/aa**2
	 width = rho0*width0*q0 / (rho*q*dcos(pi2-dmu))
	 x = xw + width*dcos(theta+dmu)
	 y = yw + width*dsin(theta+dmu)
	 xupw(i+2) = x
	 yupw(i+2) = y
  10     continue
      return
      end


text:
C    Copyright(C) 2014-2017 National Technology & Engineering Solutions of
C    Sandia, LLC (NTESS).  Under the terms of Contract DE-NA0003525 with
C    NTESS, the U.S. Government retains certain rights in this software.
C
C    Redistribution and use in source and binary forms, with or without
C    modification, are permitted provided that the following conditions are
C    met:
C
C    * Redistributions of source code must retain the above copyright
C       notice, this list of conditions and the following disclaimer.
C
C    * Redistributions in binary form must reproduce the above
C      copyright notice, this list of conditions and the following
C      disclaimer in the documentation and/or other materials provided
C      with the distribution.
C
C    * Neither the name of NTESS nor the names of its
C      contributors may be used to endorse or promote products derived
C      from this software without specific prior written permission.
C
C    THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
C    "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
C    LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
C    A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
C    OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
C    SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
C    LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
C    DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
C    THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
C    (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
C    OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
C

C $Id: nps.f,v 1.1 1990/11/30 11:12:52 gdsjaar Exp $
C $Log: nps.f,v $
C Revision 1.1  1990/11/30 11:12:52  gdsjaar
C Initial revision
C
C
CC* FILE: [.QMESH]NPS.FOR
CC* MODIFIED BY: TED BLACKER
CC* MODIFICATION DATE: 7/6/90
CC* MODIFICATION: COMPLETED HEADER INFORMATION
C
      SUBROUTINE NPS (ML, MS, MNNPS, NS, ISLIST, NINT, IFLINE, NLPS,
     &   ILLIST, LINKL, LINKS, NNPS, ERR)
C***********************************************************************
C
C  SUBROUTINE NPS = GIVES A LIST OF THE NUMBER OF PERIMETER NODES
C                   ON EACH OF A REGION'S SIDES
C
C***********************************************************************
C
C  SUBROUTINE CALLED BY:
C     QMESH = GENERATES QUAD ELEMENTS
C
C***********************************************************************
C
C  VARIABLES USED:
C     NNPS  = ARRAY OF NUMBER OF NODES PER SIDE
C     CCW   = .TRUE. IF THE SIDE IS ORIENTED CCW
C     KS    = COUNTER OF THE NUMBER OF SIDES
C
C***********************************************************************
C
      DIMENSION NNPS (MNNPS), ISLIST (NS), LINKL (2, ML), LINKS (2, MS)
      DIMENSION NLPS (MS), IFLINE (MS), ILLIST (MS * 3), NINT (ML)
C
      LOGICAL ERR, ADDLNK
C
      ERR = .TRUE.
      ADDLNK = .FALSE.
C
      KS = 0
      DO 110 I = 1, NS
         IF (ISLIST (I) .LT. 0) THEN
            KS = KS + 1
            IL = IABS (ISLIST (I))
            CALL LTSORT (ML, LINKL, IL, IPNTR, ADDLNK)
            NNPS (KS) = IABS (NINT (IPNTR)) + 1
         ELSEIF (ISLIST (I) .GT. 0) THEN
            CALL LTSORT (MS, LINKS, ISLIST (I), IPNTR, ADDLNK)
            J1 = IFLINE (IPNTR)
            J2 = J1 + NLPS (IPNTR) - 1
            KS = KS + 1
            NNPS (KS) = 0
            DO 100 J = J1, J2
               IL = ILLIST (J)
               CALL LTSORT (ML, LINKL, IL, IPNTR, ADDLNK)
               NNPS (KS) = NNPS (KS) + IABS (NINT (IPNTR))
  100       CONTINUE
            NNPS (KS) = NNPS (KS) + 1
         ELSE
            RETURN
         ENDIF
  110 CONTINUE
      ERR = .FALSE.
C
      RETURN
C
      END

max_stars_repo_name: mathstuf/seacas

max_stars_repo_licenses: BSD-3-Clause NetCDF Zlib MIT

max_issues_repo_name: mathstuf/seacas

max_issues_repo_licenses: BSD-3-Clause NetCDF Zlib MIT

max_forks_repo_name: mathstuf/seacas

max_forks_repo_licenses: BSD-3-Clause NetCDF Zlib MIT


text:
      subroutine fpback(a,z,n,k,c,nest)
c  subroutine fpback calculates the solution of the system of
c  equations a*c = z with a a n x n upper triangular matrix
c  of bandwidth k.
c  ..
c  ..scalar arguments..
      integer n,k,nest
c  ..array arguments..
      real a(nest,k),z(n),c(n)
c  ..local scalars..
      real store
      integer i,i1,j,k1,l,m
c  ..
      k1 = k-1
      c(n) = z(n)/a(n,1)
      i = n-1
      if(i.eq.0) go to 30
      do 20 j=2,n
        store = z(i)
        i1 = k1
        if(j.le.k1) i1 = j-1
        m = i
        do 10 l=1,i1
          m = m+1
          store = store-c(m)*a(i,l+1)
  10    continue
        c(i) = store/a(i,1)
        i = i-1
  20  continue
  30  return
      end

max_stars_repo_name: Astro-mh/MathematicaDierckxSurfit

max_issues_repo_name: Astro-mh/MathematicaDierckxSurfit

max_forks_repo_name: Astro-mh/MathematicaDierckxSurfit


text:
SUBROUTINE pmsh_ensi(argv,nlen,step,loads)
!
! This subroutine outputs pressure in the Ensight gold format for
! visualization in ParaView. ParaView also requires the output of subroutine
! mesh_ensi for the geometry.
!
  IMPLICIT none

  INTEGER,PARAMETER             :: iwp=SELECTED_REAL_KIND(15)
  INTEGER,   INTENT(IN)         :: nlen,step
  INTEGER                       :: i
  REAL(iwp), INTENT(IN)         :: loads(:)
  CHARACTER(LEN=15), INTENT(IN) :: argv
  CHARACTER(LEN=5)              :: ch

  WRITE(ch,'(I5.5)') step ! convert integer to string using internal file

  OPEN(17,FILE=argv(1:nlen)//'.ensi.PRESSURE-'//ch)

  WRITE(17,'(A)') "Alya Ensight Gold --- Scalar per-node variable file"
  WRITE(17,'(A/A/A)') "part", "      1","coordinates"

  DO i=1,UBOUND(loads,1)
    WRITE(17,'(E12.5)') loads(i)
  END DO

  CLOSE(17)

  RETURN

END SUBROUTINE pmsh_ensi

max_stars_repo_name: leemargetts/PoreFEM

max_issues_repo_name: leemargetts/PoreFEM

max_forks_repo_name: leemargetts/PoreFEM