text:
c
c  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
c  .                                                             .
c  .                  copyright (c) 1998 by UCAR                 .
c  .                                                             .
c  .       University Corporation for Atmospheric Research       .
c  .                                                             .
c  .                      all rights reserved                    .
c  .                                                             .
c  .                                                             .
c  .                         SPHEREPACK                       .
c  .                                                             .
c  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
c
c
c ... file shsec.f
c
c     this file contains code and documentation for subroutines
c     shsec and shseci
c
c ... files which must be loaded with shsec.f
c
c     sphcom.f, hrfft.f
c
c     subroutine shsec(nlat,nlon,isym,nt,g,idg,jdg,a,b,mdab,ndab,
c    +                    wshsec,lshsec,work,lwork,ierror)
c
c     subroutine shsec performs the spherical harmonic synthesis
c     on the arrays a and b and stores the result in the array g.
c     the synthesis is performed on an equally spaced grid.  the
c     associated legendre functions are recomputed rather than stored
c     as they are in subroutine shses.  the synthesis is described
c     below at output parameter g.
c
c     required files from spherepack2
c
c     sphcom.f, hrfft.f
c
c
c     input parameters
c
c     nlat   the number of colatitudes on the full sphere including the
c            poles. for example, nlat = 37 for a five degree grid.
c            nlat determines the grid increment in colatitude as
c            pi/(nlat-1).  if nlat is odd the equator is located at
c            grid point i=(nlat+1)/2. if nlat is even the equator is
c            located half way between points i=nlat/2 and i=nlat/2+1.
c            nlat must be at least 3. note: on the half sphere, the
c            number of grid points in the colatitudinal direction is
c            nlat/2 if nlat is even or (nlat+1)/2 if nlat is odd.
c
c     nlon   the number of distinct londitude points.  nlon determines
c            the grid increment in longitude as 2*pi/nlon. for example
c            nlon = 72 for a five degree grid. nlon must be greater
c            than or equal to 4. the efficiency of the computation is
c            improved when nlon is a product of small prime numbers.
c
c     isym   = 0  no symmetries exist about the equator. the synthesis
c                 is performed on the entire sphere.  i.e. on the
c                 array g(i,j) for i=1,...,nlat and j=1,...,nlon.
c                 (see description of g below)
c
c            = 1  g is antisymmetric about the equator. the synthesis
c                 is performed on the northern hemisphere only.  i.e.
c                 if nlat is odd the synthesis is performed on the
c                 array g(i,j) for i=1,...,(nlat+1)/2 and j=1,...,nlon.
c                 if nlat is even the synthesis is performed on the
c                 array g(i,j) for i=1,...,nlat/2 and j=1,...,nlon.
c
c
c            = 2  g is symmetric about the equator. the synthesis is
c                 performed on the northern hemisphere only.  i.e.
c                 if nlat is odd the synthesis is performed on the
c                 array g(i,j) for i=1,...,(nlat+1)/2 and j=1,...,nlon.
c                 if nlat is even the synthesis is performed on the
c                 array g(i,j) for i=1,...,nlat/2 and j=1,...,nlon.
c
c     nt     the number of syntheses.  in the program that calls shsec,
c            the arrays g,a and b can be three dimensional in which
c            case multiple syntheses will be performed.  the third
c            index is the synthesis index which assumes the values
c            k=1,...,nt.  for a single synthesis set nt=1. the
c            discription of the remaining parameters is simplified
c            by assuming that nt=1 or that the arrays g,a and b
c            have only two dimensions.
c
c     idg    the first dimension of the array g as it appears in the
c            program that calls shsec.  if isym equals zero then idg
c            must be at least nlat.  if isym is nonzero then idg
c            must be at least nlat/2 if nlat is even or at least
c            (nlat+1)/2 if nlat is odd.
c
c     jdg    the second dimension of the array g as it appears in the
c            program that calls shsec.  jdg must be at least nlon.
c
c     a,b    two or three dimensional arrays (see the input parameter
c            nt) that contain the coefficients in the spherical harmonic
c            expansion of g(i,j) given below at the definition of the
c            output parameter g.  a(m,n) and b(m,n) are defined for
c            indices m=1,...,mmax and n=m,...,nlat where mmax is the
c            maximum (plus one) longitudinal wave number given by
c            mmax = min0(nlat,(nlon+2)/2) if nlon is even or
c            mmax = min0(nlat,(nlon+1)/2) if nlon is odd.
c
c     mdab   the first dimension of the arrays a and b as it appears
c            in the program that calls shsec. mdab must be at least
c            min0(nlat,(nlon+2)/2) if nlon is even or at least
c            min0(nlat,(nlon+1)/2) if nlon is odd.
c
c     ndab   the second dimension of the arrays a and b as it appears
c            in the program that calls shsec. ndab must be at least nlat
c
c     wshsec an array which must be initialized by subroutine shseci.
c            once initialized, wshsec can be used repeatedly by shsec
c            as long as nlon and nlat remain unchanged.  wshsec must
c            not be altered between calls of shsec.
c
c     lshsec the dimension of the array wshsec as it appears in the
c            program that calls shsec. define
c
c               l1 = min0(nlat,(nlon+2)/2) if nlon is even or
c               l1 = min0(nlat,(nlon+1)/2) if nlon is odd
c
c            and
c
c               l2 = nlat/2        if nlat is even or
c               l2 = (nlat+1)/2    if nlat is odd
c
c            then lshsec must be at least
c
c            2*nlat*l2+3*((l1-2)*(nlat+nlat-l1-1))/2+nlon+15
c
c
c     work   a work array that does not have to be saved.
c
c     lwork  the dimension of the array work as it appears in the
c            program that calls shsec. define
c
c               l2 = nlat/2        if nlat is even or
c               l2 = (nlat+1)/2    if nlat is odd
c
c            if isym is zero then lwork must be at least
c
c                    nlat*(nt*nlon+max0(3*l2,nlon))
c
c            if isym is not zero then lwork must be at least
c
c                    l2*(nt*nlon+max0(3*nlat,nlon))
c
c     **************************************************************
c
c     output parameters
c
c     g      a two or three dimensional array (see input parameter
c            nt) that contains the spherical harmonic synthesis of
c            the arrays a and b at the colatitude point theta(i) =
c            (i-1)*pi/(nlat-1) and longitude point phi(j) =
c            (j-1)*2*pi/nlon. the index ranges are defined above at
c            at the input parameter isym.  for isym=0, g(i,j) is
c            given by the the equations listed below.  symmetric
c            versions are used when isym is greater than zero.
c
c     the normalized associated legendre functions are given by
c
c     pbar(m,n,theta) = sqrt((2*n+1)*factorial(n-m)/(2*factorial(n+m)))
c                       *sin(theta)**m/(2**n*factorial(n)) times the
c                       (n+m)th derivative of (x**2-1)**n with respect
c                       to x=cos(theta)
c
c     define the maximum (plus one) longitudinal wave number
c     as   mmax = min0(nlat,(nlon+2)/2) if nlon is even or
c          mmax = min0(nlat,(nlon+1)/2) if nlon is odd.
c
c     then g(i,j) = the sum from n=0 to n=nlat-1 of
c
c                   .5*pbar(0,n,theta(i))*a(1,n+1)
c
c              plus the sum from m=1 to m=mmax-1 of
c
c                   the sum from n=m to n=nlat-1 of
c
c              pbar(m,n,theta(i))*(a(m+1,n+1)*cos(m*phi(j))
c                                    -b(m+1,n+1)*sin(m*phi(j)))
c
c
c     ierror = 0  no errors
c            = 1  error in the specification of nlat
c            = 2  error in the specification of nlon
c            = 3  error in the specification of isym
c            = 4  error in the specification of nt
c            = 5  error in the specification of idg
c            = 6  error in the specification of jdg
c            = 7  error in the specification of mdab
c            = 8  error in the specification of ndab
c            = 9  error in the specification of lshsec
c            = 10 error in the specification of lwork
c
c
c ****************************************************************
c     subroutine shseci(nlat,nlon,wshsec,lshsec,dwork,ldwork,ierror)
c
c     subroutine shseci initializes the array wshsec which can then
c     be used repeatedly by subroutine shsec.
c
c     input parameters
c
c     nlat   the number of colatitudes on the full sphere including the
c            poles. for example, nlat = 37 for a five degree grid.
c            nlat determines the grid increment in colatitude as
c            pi/(nlat-1).  if nlat is odd the equator is located at
c            grid point i=(nlat+1)/2. if nlat is even the equator is
c            located half way between points i=nlat/2 and i=nlat/2+1.
c            nlat must be at least 3. note: on the half sphere, the
c            number of grid points in the colatitudinal direction is
c            nlat/2 if nlat is even or (nlat+1)/2 if nlat is odd.
c
c     nlon   the number of distinct londitude points.  nlon determines
c            the grid increment in longitude as 2*pi/nlon. for example
c            nlon = 72 for a five degree grid. nlon must be greater
c            than or equal to 4. the efficiency of the computation is
c            improved when nlon is a product of small prime numbers.
c
c     lshsec the dimension of the array wshsec as it appears in the
c            program that calls shseci. the array wshsec is an output
c            parameter which is described below. define
c
c               l1 = min0(nlat,(nlon+2)/2) if nlon is even or
c               l1 = min0(nlat,(nlon+1)/2) if nlon is odd
c
c            and
c
c               l2 = nlat/2        if nlat is even or
c               l2 = (nlat+1)/2    if nlat is odd
c
c            then lshsec must be at least
c
c            2*nlat*l2+3*((l1-2)*(nlat+nlat-l1-1))/2+nlon+15
c
c     dwork  a double precision work array that does not have to be
c            saved.
c
c     ldwork the dimension of array dwork as it appears in the program
c            that calls shseci.  ldwork must be at least nlat+1.
c
c     output parameters
c
c     wshsec an array which is initialized for use by subroutine shsec.
c            once initialized, wshsec can be used repeatedly by shsec
c            as long as nlon and nlat remain unchanged.  wshsec must
c            not be altered between calls of shsec.
c
c     ierror = 0  no errors
c            = 1  error in the specification of nlat
c            = 2  error in the specification of nlon
c            = 3  error in the specification of lshsec
c            = 4  error in the specification of ldwork
c
c
c ****************************************************************
      subroutine shsec(nlat,nlon,isym,nt,g,idg,jdg,a,b,mdab,ndab,
     1                    wshsec,lshsec,work,lwork,ierror)
      dimension g(idg,jdg,1),a(mdab,ndab,1),b(mdab,ndab,1),wshsec(1),
     1          work(1)
      ierror = 1
      if(nlat.lt.3) return
      ierror = 2
      if(nlon.lt.4) return
      ierror = 3
      if(isym.lt.0 .or. isym.gt.2) return
      ierror = 4
      if(nt .lt. 0) return
      ierror = 5
      if((isym.eq.0 .and. idg.lt.nlat) .or.
     1   (isym.ne.0 .and. idg.lt.(nlat+1)/2)) return
      ierror = 6
      if(jdg .lt. nlon) return
      ierror = 7
      mmax = min0(nlat,nlon/2+1)
      if(mdab .lt. mmax) return
      ierror = 8
      if(ndab .lt. nlat) return
      ierror = 9
      imid = (nlat+1)/2
      lzz1 = 2*nlat*imid
      labc = 3*((mmax-2)*(nlat+nlat-mmax-1))/2
      if(lshsec .lt. lzz1+labc+nlon+15) return
      ierror = 10
      ls = nlat
      if(isym .gt. 0) ls = imid
      nln = nt*ls*nlon
      if(lwork .lt. nln+max0(ls*nlon,3*nlat*imid)) return
      ierror = 0
      ist = 0
      if(isym .eq. 0) ist = imid
      iw1 = lzz1+labc+1
      call shsec1(nlat,isym,nt,g,idg,jdg,a,b,mdab,ndab,imid,ls,nlon,
     1     work,work(ist+1),work(nln+1),work(nln+1),wshsec,wshsec(iw1))
      return
      end
      subroutine shsec1(nlat,isym,nt,g,idgs,jdgs,a,b,mdab,ndab,imid,
     1                idg,jdg,ge,go,work,pb,walin,whrfft)
c
c     whrfft must have at least nlon+15 locations
c     walin must have 3*l*imid + 3*((l-3)*l+2)/2 locations
c     zb must have 3*l*imid locations
c
      dimension g(idgs,jdgs,1),a(mdab,ndab,1),b(mdab,ndab,1),
     1          ge(idg,jdg,1),go(idg,jdg,1),pb(imid,nlat,3),walin(1),
     3          whrfft(1),work(1)
      ls = idg
      nlon = jdg
      mmax = min0(nlat,nlon/2+1)
      mdo = mmax
      if(mdo+mdo-1 .gt. nlon) mdo = mmax-1
      nlp1 = nlat+1
      modl = mod(nlat,2)
      imm1 = imid
      if(modl .ne. 0) imm1 = imid-1
      do 80 k=1,nt
      do 80 j=1,nlon
      do 80 i=1,ls
      ge(i,j,k)=0.
   80 continue
      if(isym .eq. 1) go to 125
      call alin (2,nlat,nlon,0,pb,i3,walin)
      do 100 k=1,nt
      do 100 np1=1,nlat,2
      do 100 i=1,imid
      ge(i,1,k)=ge(i,1,k)+a(1,np1,k)*pb(i,np1,i3)
  100 continue
      ndo = nlat
      if(mod(nlat,2) .eq. 0) ndo = nlat-1
      do 110 mp1=2,mdo
      m = mp1-1
      call alin (2,nlat,nlon,m,pb,i3,walin)
      do 110 np1=mp1,ndo,2
      do 110 k=1,nt
      do 110 i=1,imid
      ge(i,2*mp1-2,k) = ge(i,2*mp1-2,k)+a(mp1,np1,k)*pb(i,np1,i3)
      ge(i,2*mp1-1,k) = ge(i,2*mp1-1,k)+b(mp1,np1,k)*pb(i,np1,i3)
  110 continue
      if(mdo .eq. mmax .or. mmax .gt. ndo) go to 122
      call alin (2,nlat,nlon,mdo,pb,i3,walin)
      do 120 np1=mmax,ndo,2
      do 120 k=1,nt
      do 120 i=1,imid
      ge(i,2*mmax-2,k) = ge(i,2*mmax-2,k)+a(mmax,np1,k)*pb(i,np1,i3)
  120 continue
  122 if(isym .eq. 2) go to 155
  125 call alin(1,nlat,nlon,0,pb,i3,walin)
      do 140 k=1,nt
      do 140 np1=2,nlat,2
      do 140 i=1,imm1
      go(i,1,k)=go(i,1,k)+a(1,np1,k)*pb(i,np1,i3)
  140 continue
      ndo = nlat
      if(mod(nlat,2) .ne. 0) ndo = nlat-1
      do 150 mp1=2,mdo
      mp2 = mp1+1
      m = mp1-1
      call alin(1,nlat,nlon,m,pb,i3,walin)
      do 150 np1=mp2,ndo,2
      do 150 k=1,nt
      do 150 i=1,imm1
      go(i,2*mp1-2,k) = go(i,2*mp1-2,k)+a(mp1,np1,k)*pb(i,np1,i3)
      go(i,2*mp1-1,k) = go(i,2*mp1-1,k)+b(mp1,np1,k)*pb(i,np1,i3)
  150 continue
      mp2 = mmax+1
      if(mdo .eq. mmax .or. mp2 .gt. ndo) go to 155
      call alin(1,nlat,nlon,mdo,pb,i3,walin)
      do 152 np1=mp2,ndo,2
      do 152 k=1,nt
      do 152 i=1,imm1
      go(i,2*mmax-2,k) = go(i,2*mmax-2,k)+a(mmax,np1,k)*pb(i,np1,i3)
  152 continue
  155 do 160 k=1,nt
      if(mod(nlon,2) .ne. 0) go to 157
      do 156 i=1,ls
      ge(i,nlon,k) = 2.*ge(i,nlon,k)
  156 continue
  157 call hrfftb(ls,nlon,ge(1,1,k),ls,whrfft,work)
  160 continue
      if(isym .ne. 0) go to 180
      do 170 k=1,nt
      do 170 j=1,nlon
      do 175 i=1,imm1
      g(i,j,k) = .5*(ge(i,j,k)+go(i,j,k))
      g(nlp1-i,j,k) = .5*(ge(i,j,k)-go(i,j,k))
  175 continue
      if(modl .eq. 0) go to 170
      g(imid,j,k) = .5*ge(imid,j,k)
  170 continue
      return
  180 do 185 k=1,nt
      do 185 i=1,imid
      do 185 j=1,nlon
      g(i,j,k) = .5*ge(i,j,k)
  185 continue
      return
      end
c     subroutine shseci(nlat,nlon,wshsec,lshsec,dwork,ldwork,ierror)
c
c     subroutine shseci initializes the array wshsec which can then
c     be used repeatedly by subroutine shsec.
c
c     input parameters
c
c     nlat   the number of colatitudes on the full sphere including the
c            poles. for example, nlat = 37 for a five degree grid.
c            nlat determines the grid increment in colatitude as
c            pi/(nlat-1).  if nlat is odd the equator is located at
c            grid point i=(nlat+1)/2. if nlat is even the equator is
c            located half way between points i=nlat/2 and i=nlat/2+1.
c            nlat must be at least 3. note: on the half sphere, the
c            number of grid points in the colatitudinal direction is
c            nlat/2 if nlat is even or (nlat+1)/2 if nlat is odd.
c
c     nlon   the number of distinct londitude points.  nlon determines
c            the grid increment in longitude as 2*pi/nlon. for example
c            nlon = 72 for a five degree grid. nlon must be greater
c            than or equal to 4. the efficiency of the computation is
c            improved when nlon is a product of small prime numbers.
c
c     lshsec the dimension of the array wshsec as it appears in the
c            program that calls shseci. the array wshsec is an output
c            parameter which is described below. define
c
c               l1 = min0(nlat,(nlon+2)/2) if nlon is even or
c               l1 = min0(nlat,(nlon+1)/2) if nlon is odd
c
c            and
c
c               l2 = nlat/2        if nlat is even or
c               l2 = (nlat+1)/2    if nlat is odd
c
c            then lshsec must be at least
c
c            2*nlat*l2+3*((l1-2)*(nlat+nlat-l1-1))/2+nlon+15
c
c     dwork  a double precision work array that does not have to be
c            saved.
c
c     ldwork the dimension of array dwork as it appears in the program
c            that calls shseci.  ldwork must be at least nlat+1.
c
c     output parameters
c
c     wshsec an array which is initialized for use by subroutine shsec.
c            once initialized, wshsec can be used repeatedly by shsec
c            as long as nlon and nlat remain unchanged.  wshsec must
c            not be altered between calls of shsec.
c
c     ierror = 0  no errors
c            = 1  error in the specification of nlat
c            = 2  error in the specification of nlon
c            = 3  error in the specification of lshsec
c            = 4  error in the specification of ldwork
c
c
c ****************************************************************
      subroutine shseci(nlat,nlon,wshsec,lshsec,dwork,ldwork,ierror)
      dimension wshsec(*)
      double precision dwork(ldwork)
      ierror = 1
      if(nlat.lt.3) return
      ierror = 2
      if(nlon.lt.4) return
      ierror = 3
      imid = (nlat+1)/2
      mmax = min0(nlat,nlon/2+1)
      lzz1 = 2*nlat*imid
      labc = 3*((mmax-2)*(nlat+nlat-mmax-1))/2
      if(lshsec .lt. lzz1+labc+nlon+15) return
      ierror = 4
      if(ldwork .lt. nlat+1) return
      ierror = 0
      call alinit(nlat,nlon,wshsec,dwork)
      iw1 = lzz1+labc+1
      call hrffti(nlon,wshsec(iw1))
      return
      end


text:
subroutine Finbou(a,pnodb,knodb,iboun)
!-----------------------------------------------------------------------
!    This routine looks for iboun in in the list knodb using
!    a%pbopo and a%lbopo
!-----------------------------------------------------------------------
  use typre
  use Mod_Mesh
  implicit none

  class(FemMesh) :: a
  integer(ip) :: pnodb,knodb(pnodb)
  integer(ip) :: iboun

  integer(ip)              :: onodb,inook,ibook,ibopo,ielem,inodb,inode,ispos,ispos2,jnodb
  integer(ip)              :: qnodb,jpoin,pblty

  iboun=0
  ibook=0
  pblty=1
  !loop_iboun: do while(iboun<nboun)
  !   iboun=iboun+1
  loop_iboun: do ibopo = a%pbopo(knodb(1)),a%pbopo(knodb(1)+1)-1
	iboun = a%lbopo(ibopo)

     ispos2 = (a%nnodb(1)+1)*(iboun-1)
     if (a%nelty > 1) then
       ispos2 = a%pboel(iboun)-1
       pblty = a%ltypb(a%pboel(iboun+1)-a%pboel(iboun)-1)
     endif
     ielem = a%lboel(ispos2+a%nnodb(pblty)+1)
     ispos = a%nnode(1)*(ielem-1)
     if (a%nelty > 1) ispos = a%pnods(ielem)-1

     qnodb=a%nnodb(pblty)
     if(pnodb==qnodb) then
        onodb=0
        jnodb=0
        inook=0
        loop_jnodb: do while(jnodb<pnodb)
           jnodb=jnodb+1
           jpoin=knodb(jnodb)
           inodb=0
           loop_inodb: do while(inodb<qnodb)
              inodb=inodb+1
              inode = a%lboel(ispos2+inodb)
              if(jpoin==a%lnods(ispos+inode)) then
                 onodb=onodb+1
                 inook=1
                 exit loop_inodb
              end if
           end do loop_inodb
           if(inook/=1) exit loop_jnodb
        end do loop_jnodb
        if(onodb==pnodb) then
           ibook=1
           exit loop_iboun
        end if
     end if
  end do loop_iboun
  if(ibook/=1) iboun=0

end subroutine Finbou

max_stars_repo_name: ciaid-colombia/InsFEM

max_issues_repo_name: ciaid-colombia/InsFEM

max_forks_repo_name: ciaid-colombia/InsFEM


#include "AMYCAI.h"

MODULE Geology_Information

    USE TimeCatcher
    USE FloatPrecision
    USE ArraySize
    USE GRIDS
    USE Parameters

    IMPLICIT NONE

    CONTAINS

! = = = = = = = = = = = = = = = = =
!
! assignment to Seabed Parameters
!
! = = = = = = = = = = = = = = = = =

    SUBROUTINE BedParameters( DEP0, MDEP0, DDEP, Upsilon, MAREA)

        IMPLICIT NONE

        REAL(fp) , ALLOCATABLE , DIMENSION(:,:)   ,          INTENT(IN)    :: DEP0
        REAL(fp) , ALLOCATABLE , DIMENSION(:)     ,          INTENT(INOUT) :: MDEP0
        REAL(fp) , ALLOCATABLE , DIMENSION(:,:)   ,          INTENT(INOUT) :: DDEP
        REAL(fp) , ALLOCATABLE , DIMENSION(:,:)   , TARGET , INTENT(INOUT) :: Upsilon
        REAL(fp) , ALLOCATABLE , DIMENSION(:)     ,          INTENT(INOUT) :: MAREA

        INTEGER           :: I , J , K
        INTEGER           :: status
        CHARACTER(len=80) :: err_msg

! - - - - - - - - - - - - - - - - - -
!
!   Calculate the mean water Depth
!
! - - - - - - - - - - - - - - - - - -

        MDEP0(:) = (2.D0*SUM(DEP0,DIM=2) - DEP0(:,NL) - DEP0(:,NU)) / 2.D0 * DX / BT

! - - - - - - - - - - - - - - - - - -
!
!   Calculate the MeanArea
!
! - - - - - - - - - - - - - - - - - -

        DO I = ML , MU
            MAREA(I) = 0.D0
            DO J = NL , NU-1
                MAREA(I) = MAREA(I) + (DMAX1(DEP0(I,J),DDRY)+DMAX1(DEP0(I,J+1),DDRY))/2.D0*DX
            ENDDO
        ENDDO

! - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!
!   Get the lateral slope and the factors of the shear stress (Upsilon)
!
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

        DO I = ML , MU
            DO J = NL , NU
                IF(J==NU) THEN
                    DDEP(I,J) = 0.D0
                ELSE
                    DDEP(I,J) = (DEP0(I,J+1)-DEP0(I,J)) / dx
                END IF
                Upsilon(I,J) = SQRT(1+DDEP(I,J)**2)
            END DO
        END DO

! sediment stratification

    END SUBROUTINE BedParameters

! = = = = = = = = = = = = = = = = =
!
! assignment to Sedimentation Parameters
! These parameters serve the sediment transport
! Model
!
! = = = = = = = = = = = = = = = = =

    SUBROUTINE Sedimentation_Parameters(DEP0,TDEP0,HTDEP0,MDEP0,HMTDEP0,HETA,HND,HNDD,HDIFF)

        IMPLICIT NONE

        REAL(fp) , ALLOCATABLE , DIMENSION(:,:) ,          INTENT(IN)    :: DEP0
        REAL(fp) , ALLOCATABLE , DIMENSION(:,:) ,          INTENT(IN)    :: TDEP0
        REAL(fp) , ALLOCATABLE , DIMENSION(:,:) ,          INTENT(INOUT) :: HTDEP0
        REAL(fp) , ALLOCATABLE , DIMENSION(:)   ,          INTENT(IN)    :: MDEP0
        REAL(fp) , ALLOCATABLE , DIMENSION(:)   ,          INTENT(INOUT) :: HMTDEP0
        REAL(fp) , ALLOCATABLE , DIMENSION(:)   ,          INTENT(INOUT) :: HETA
        INTEGER  , ALLOCATABLE , DIMENSION(:)   ,          INTENT(INOUT) :: HND
        INTEGER  , ALLOCATABLE , DIMENSION(:,:) ,          INTENT(INOUT) :: HNDD
        REAL(fp) , ALLOCATABLE , DIMENSION(:)   ,          INTENT(INOUT) :: HDIFF
        !---------------------------------------------------------------------------
        INTEGER                                                          :: I , J
        INTEGER                                                          :: status
        CHARACTER(LEN=80)                                                :: err_msg

#ifndef SSC_CONSTANT

    !- - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    ! get the water Depth at high tide
    !- - - - - - - - - - - - - - - - - - - - - - - - - - - - -

        DO I = ML , MU
            HTDEP0(I,:) = DEP0(I,:) + TDAMP0
        END DO

    ! - - - - - - - - - - - - -
    !   get the HWL Mean Total water depth
    ! - - - - - - - - - - - - -

        HMTDEP0(:) = DEP0(:,NU) + TDAMP0

    !- - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    ! get the wet AND dry at high tide
    !- - - - - - - - - - - - - - - - - - - - - - - - - - - - -

        WHERE(HTDEP0 > 0.D0)        ! wet
            HNDD = 0
        ELSEWHERE(HTDEP0 <= 0.D0)   ! dry
            HNDD = 1
        END WHERE

        WHERE(HMTDEP0 > 0.D0)        ! wet
            HND = 0
        ELSEWHERE(HMTDEP0 <= 0.D0)   ! dry
            HND = 1
        END WHERE
        HND(ML) = -5

    !- - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    !  Get the diffusion coefficient for SSC transport.
    !  The locations of the diffusion coefficients are between
    !  grid points... 1D
    !- - - - - - - - - - - - - - - - - - - - - - - - - - - - -

!    ! diffusion in y direcion (Longitudinal)
!
!        DO I = ML , MU-1
!            IF(HMTDEP0(I)>0.D0.AND.HMTDEP0(I+1)>0.D0) THEN
!                HDIFF(I) = TD_Excursion_Y*0.5D0*(HMTDEP0(I)+HMTDEP0(I+1))
!            ELSE
!                HDIFF(I) = 0.D0
!            END IF
!        END DO
!
!    ! Reset Maximum eta
!
!        IF(T0>=TDCYC) THEN
!            HETA(:)  = 0.D0
!        END IF
!
#else

        WHERE(TDEP0 > 0.D0)        ! wet
            HNDD = 0
        ELSEWHERE(TDEP0 <= 0.D0)   ! dry
            HNDD = 1
        END WHERE

#endif

    !- - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    !  Get the diffusion coefficient for SSC transport.
    !  The locations of the diffusion coefficients are between
    !  grid points... 2D
    !- - - - - - - - - - - - - - - - - - - - - - - - - - - - -

    ! diffusion in x direcion (transversal)

        !DO I = ML , MU
        !    DO J = NL , NU-1
        !        IF(HTDEP0(I,J)>0.D0.AND.HTDEP0(I,J+1)>0.D0) THEN
        !            SSC_DIFFU_HWL_X(I,J) = TD_Excursion_X*0.5D0*(HTDEP0(I,J)+HTDEP0(I,J+1))
        !        ELSE
        !            SSC_DIFFU_HWL_X(I,J) = 0.D0
        !        END IF
        !    END DO
        !END DO

    ! diffusion in y direcion (Longitudinal)

        !DO I = ML , MU-1
        !    DO J = NL , NU
        !        IF(HTDEP0(I,J)>0.D0.AND.HTDEP0(I+1,J)>0.D0) THEN
        !            SSC_DIFFU_HWL_Y(I,J) = TD_Excursion_Y*0.5D0*(HTDEP0(I,J)+HTDEP0(I+1,J))
        !        ELSE
        !            SSC_DIFFU_HWL_Y(I,J) = 0.D0
        !        END IF
        !    END DO
        !END DO

    END SUBROUTINE Sedimentation_Parameters

END MODULE Geology_Information


text:
C Test compiler flags: -ff90
C Origin: David Billinghurst <David.Billinghurst@riotinto.com>
C
C Read the g77 manual entry on CMPAMBIG
C
C { dg-do run }
C { dg-options "-ff90" }
      double complex z
      z = (2.0d0,1.0d0)
      call s(real(z))
      end
      subroutine s(x)
      double precision x
      if ( abs(x-2.0d0) .gt. 1.0e-5 ) call abort
      end

ext: f


text:
      program grid_rosby
      use mod_za  ! HYCOM array I/O interface
      implicit none
c
      integer   i,j
      real      eqdist,grdlon,rdist,hmaxa,hmaxb,hmina,hminb
      character cline*80
c
c --- read in a hycom grid file,
c --- write out a field proportional to how well Rosby Radius is resolved.
c --- In particular: grid / (grdlon*111200*cos(min(plat,85)))
c ---  where grid = pscx, pscy, or max(pscx,pscy)
c
      real, parameter :: radian = 57.2957795
c
      integer, allocatable :: ip(:,:)
      real,    allocatable :: px(:,:),py(:,:),pl(:,:)
c
      call xcspmd  !input idm,jdm
      allocate( ip(idm,jdm) )
      allocate( px(idm,jdm) )
      allocate( py(idm,jdm) )
      allocate( pl(idm,jdm) )
c
c --- input grdlon
c
      read(5,*) grdlon
c
c --- read in regional.grid
c
      call zaiost
c
      call zhopnc(21, 'regional.grid.b',  'formatted', 'old', 0)
      call zhopnc(31, 'regional.rosby.b', 'formatted', 'new', 0)
      call zaiopf('regional.grid.a',  'old', 21)
      call zaiopf('regional.rosby.a', 'new', 31)
c
      read( 21,'(a)') cline ! idm
      write(31,'(a)') cline
      write( 6,'(a)') cline
      read( 21,'(a)') cline ! jdm
      write(31,'(a)') cline
      write( 6,'(a)') cline
      write(31,'(f8.2,a)')
     &  grdlon,
     &  " 'grdlon' = nominal equatorial grid spacing in degrees"
      write( 6,'(f8.2,a)')
     &  grdlon,
     &  " 'grdlon' = nominal equatorial grid spacing in degrees"
c
      read(21,*) ! skip mapflg
      read(21,*) ! skip plon
      call zaiosk(21)
c
      read(21,'(a)') cline
      write(6,'(a)') cline(1:len_trim(cline))
      i = index(cline,'=')
      read (cline(i+1:),*)   hminb,hmaxb
      call zaiord(pl,ip,.false., hmina,hmaxa, 21)
      if     (abs(hmina-hminb).gt.abs(hminb)*1.e-4 .or.
     &        abs(hmaxa-hmaxb).gt.abs(hmaxb)*1.e-4     ) then
        write(6,'(/ a / a,1p3e14.6 / a,1p3e14.6 /)')
     &    'error - .a and .b grid files not consistent (plat):',
     &    '.a,.b min = ',hmina,hminb,hmina-hminb,
     &    '.a,.b max = ',hmaxa,hmaxb,hmaxa-hmaxb
        call zhflsh(6)
        stop
      endif
c
      do i= 1,7  ! skip qlon,qlat,ulon,ulat,vlon,vlat,pang
        read(21,*)
        call zaiosk(21)
      enddo
c
      read(21,'(a)') cline
      write(6,'(a)') cline(1:len_trim(cline))
      i = index(cline,'=')
      read (cline(i+1:),*)   hminb,hmaxb
      call zaiord(px,ip,.false., hmina,hmaxa, 21)
      if     (abs(hmina-hminb).gt.abs(hminb)*1.e-4 .or.
     &        abs(hmaxa-hmaxb).gt.abs(hmaxb)*1.e-4     ) then
        write(6,'(/ a / a,1p3e14.6 / a,1p3e14.6 /)')
     &    'error - .a and .b grid files not consistent (pscx):',
     &    '.a,.b min = ',hmina,hminb,hmina-hminb,
     &    '.a,.b max = ',hmaxa,hmaxb,hmaxa-hmaxb
        call zhflsh(6)
        stop
      endif
c
      read(21,'(a)') cline
      write(6,'(a)') cline(1:len_trim(cline))
      i = index(cline,'=')
      read (cline(i+1:),*)   hminb,hmaxb
      call zaiord(py,ip,.false., hmina,hmaxa, 21)
      if     (abs(hmina-hminb).gt.abs(hminb)*1.e-4 .or.
     &        abs(hmaxa-hmaxb).gt.abs(hmaxb)*1.e-4     ) then
        write(6,'(/ a / a,1p3e14.6 / a,1p3e14.6 /)')
     &    'error - .a and .b grid files not consistent (pscy):',
     &    '.a,.b min = ',hmina,hminb,hmina-hminb,
     &    '.a,.b max = ',hmaxa,hmaxb,hmaxa-hmaxb
        call zhflsh(6)
        stop
      endif
c
      close(unit=21)
      call zaiocl(21)
c
c     calculate result.
c
      eqdist = grdlon*111200.0
      write(6,*) 'eqdist = ',eqdist
      do j= 1,jdm
        do i= 1,idm
          rdist = 1.0/(eqdist*cos(min(abs(pl(i,j)),85.0)/radian))
          pl(i,j) = max(px(i,j),py(i,j))*rdist
          px(i,j) =     px(i,j)         *rdist
          py(i,j) =             py(i,j) *rdist
        enddo
      enddo
c
c     output results.
c
      call zaiowr(px, ip,.false., hmina,hmaxa, 31, .false.)
      write(31,6100) hmina,hmaxa
      write(6, 6100) hmina,hmaxa
 6100 format('rbyx:  min,max = ',2f10.3)
c
      call zaiowr(py, ip,.false., hmina,hmaxa, 31, .false.)
      write(31,6110) hmina,hmaxa
      write(6, 6110) hmina,hmaxa
 6110 format('rbyy:  min,max = ',2f10.3)
c
      call zaiowr(pl, ip,.false., hmina,hmaxa, 31, .false.)
      write(31,6120) hmina,hmaxa
      write(6, 6120) hmina,hmaxa
 6120 format('rbym:  min,max = ',2f10.3)
      write(6, *)
c
      close(unit=31)
      call zaiocl(31)
c
      end

max_stars_repo_name: TillRasmussen/HYCOM-tools

max_issues_repo_name: TillRasmussen/HYCOM-tools

max_forks_repo_name: TillRasmussen/HYCOM-tools


text:
      FUNCTION dpmpar (i)
      USE stel_kinds
      IMPLICIT NONE
C-----------------------------------------------
C   D u m m y   A r g u m e n t s
C-----------------------------------------------
      INTEGER :: i
      REAL(rprec) :: dpmpar
C-----------------------------------------------
c
c     FUNCTION dpmpar
c
c     this FUNCTION provides single (double) precision machine parameters
c     when the appropriate set of data statements is activated (by
c     removing the c from column 1) and ALL other data statements are
c     rendered inactive. most of the PARAMETER values were obtained
c     from the corresponding bell laboratories port library FUNCTION.
c
c     the FUNCTION statement is
c
c       FUNCTION dpmpar(i)
c
c     WHERE
c
c       i is an integer input variable set to 1, 2, or 3 which
c         selects the desired machine parameter. if the machine has
c         t base b digits and its smallest and largest exponents are
c         emin and emax, respectively, then these parameters are
c
c         dpmpar(1) = b**(1 - t), the machine precision,
c
c         dpmpar(2) = b**(emin - 1), the smallest magnitude,
c
c         dpmpar(3) = b**emax*(1 - b**(-t)), the largest magnitude.
c
c     argonne national laboratory. minpack project. march 1980.
c     burton s. garbow, kenneth e. hillstrom, jorge j. more
c
c     SINGLE PRECISION
c     data rmach(1)/1.490116100E-8/
c     data rmach(2)/1.4693679000E-30/
c     data rmach(3)/1.701411800E+30/
c     DOUBLE PRECISION - IEEE
c     data dmach(1) /2.22044604926d-16/
c     data dmach(2) /2.22507385852d-308/
c     data dmach(3) /1.79769313485d+308/
c     modified for f90 (sph, august 1997)
c
      SELECT CASE(i)
      CASE(:1)
        dpmpar = EPSILON(dpmpar)      !2.22044604926e-16_dp
      CASE(2)
        dpmpar = TINY(dpmpar)         !2.22507385852e-308_dp
      CASE(3:)
        dpmpar = HUGE(dpmpar)         !1.79769313485e+308_dp
      END SELECT

      END FUNCTION dpmpar


text:
program test
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    !! Perform unittests of the AQUA-Equation of state module
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    ! Copyright (c) 2020, Jonas Haldemann
    !
    ! Permission is hereby granted, free of charge, to any person obtaining a copy of this software
    ! and associated documentation files (the "Software"), to deal in the Software without restriction,
    ! including without limitation the rights to use, copy, modify, merge, publish, distribute,
    ! sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is
    ! furnished to do so, subject to the following conditions:
    !
    ! The above copyright notice and this permission notice shall be included in all copies or
    ! substantial portions of the Software.
    !
    ! THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT
    ! NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
    ! NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM,
    ! DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
    ! OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.'
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    implicit none

    integer,dimension(7) :: status_arr

    call test_LOAD_TABLE_PT(status_arr(1))
    call test_LOAD_TABLE_RhoT(status_arr(2))
    call test_LOAD_TABLE_RhoU(status_arr(3))

    print*,'------------------------------------------------------------'
    print*,'Test results:'
    print*,'------------------------------------------------------------'
    if (status_arr(1).eq.0)then
        print*,'    test_LOAD_TABLE_PT():          Passed'
    else
        print*,'    test_LOAD_TABLE_PT():          Failed, error:',status_arr(1)
    endif

    if (status_arr(2).eq.0)then
        print*,'    test_LOAD_TABLE_RhoT():        Passed'
    else
        print*,'    test_LOAD_TABLE_RhoT():        Failed, error:',status_arr(2)
    endif

    if (status_arr(3).eq.0)then
        print*,'    test_LOAD_TABLE_RhoU():        Passed'
    else
        print*,'    test_LOAD_TABLE_RhoU():        Failed, error:',status_arr(3)
    endif

    call test_BILINEAR_INTERPOLATION(status_arr(4))
    if (status_arr(4).eq.0)then
        print*,'    test_BILINEAR_INTERPOLATION(): Passed'
    else
        print*,'    test_BILINEAR_INTERPOLATION(): Failed, error:',status_arr(4)
    endif

    call test_INTERPOLATE_AQUA_PT(status_arr(5))
    if (status_arr(5).eq.0)then
        print*,'    test_INTERPOLATE_AQUA_PT():    Passed'
    else
        print*,'    test_INTERPOLATE_AQUA_PT():    Failed, error:',status_arr(5)
    endif

    call test_INTERPOLATE_AQUA_RhoT(status_arr(6))
    if (status_arr(6).eq.0)then
        print*,'    test_INTERPOLATE_AQUA_RhoT():  Passed'
    else
        print*,'    test_INTERPOLATE_AQUA_RhoT():  Failed, error:',status_arr(6)
    endif
print*,'------------------------------------------------------------'
if (any(status_arr.gt.0)) stop 1
end program test

subroutine test_LOAD_TABLE_PT(status)
    use aqua_eos,only:LOAD_TABLE_PT
    implicit none

    integer, intent(out) :: status
    CHARACTER(LEN=200) :: path

    status = 1

    path = '../Tables'
    CALL LOAD_TABLE_PT(path)

    ! if table could be loaded test was successfull
    status = 0
    return
end subroutine test_LOAD_TABLE_PT

subroutine test_LOAD_TABLE_RhoT(status)
    use aqua_eos,only:LOAD_TABLE_RhoT
    implicit none

    integer, intent(out) :: status
    CHARACTER(LEN=200) :: path

    status = 1

    path = '../Tables'
    CALL LOAD_TABLE_RhoT(path)

    ! if table could be loaded test was successfull
    status = 0
    return
end subroutine test_LOAD_TABLE_RhoT

subroutine test_LOAD_TABLE_RhoU(status)
    use aqua_eos,only:LOAD_TABLE_RhoU
    implicit none

    integer, intent(out) :: status
    CHARACTER(LEN=200) :: path

    status = 1

    path = '../Tables'
    CALL LOAD_TABLE_RhoU(path)

    ! if table could be loaded test was successfull
    status = 0
    return
end subroutine test_LOAD_TABLE_RhoU

subroutine test_BILINEAR_INTERPOLATION(status)
    use aqua_eos,only:BILINEAR_INTERPOLATION
    implicit none

    integer, intent(out) :: status
    !REAL(8),DIMENSION(size(zarr,1))     :: z
    real(8) :: x,x_arr(3)
    real(8) :: y,y_arr(3)
    real(8) :: z_arr2(2,3,3),z2(2),z_true
    status = 0

    x = 0.5_8
    y = 0.5_8
    x_arr = [0._8,1._8,2._8]
    y_arr = [0._8,1._8,2._8]
    z_arr2(:,:,1) = transpose(reshape([0.1_8,0.2_8,0.3_8,0.2_8,0.3_8,0.4_8],[3,2]))
    z_arr2(:,:,2) = transpose(reshape([0.2_8,0.3_8,0.4_8,0.2_8,0.3_8,0.4_8],[3,2]))
    z_arr2(:,:,3) = transpose(reshape([0.3_8,0.4_8,0.5_8,0.2_8,0.3_8,0.4_8],[3,2]))

    z2 = BILINEAR_INTERPOLATION(x,y,x_arr,y_arr,z_arr2)
    z_true = 0.2_8
    if (abs(z_true - z2(1)).gt. 1E-6_8 ) status = 1
    z_true = 0.25_8
    if (abs(z_true - z2(2)).gt. 1E-6_8 ) status = 2

    x = 1.0_8
    y = 0.25_8
    z2 = BILINEAR_INTERPOLATION(x,y,x_arr,y_arr,z_arr2)
    z_true = 0.225_8
    if (abs(z_true - z2(1)).gt. 1E-6_8 ) status = 3
    z_true = 0.225_8
    if (abs(z_true - z2(2)).gt. 1E-6_8 ) status = 4

    x = 1.5_8
    y = 1._8
    z2 = BILINEAR_INTERPOLATION(x,y,x_arr,y_arr,z_arr2)
    z_true = 0.35_8
    if (abs(z_true - z2(1)).gt. 1E-6_8 ) status = 3
    z_true = 0.3_8
    if (abs(z_true - z2(2)).gt. 1E-6_8 ) status = 4

    x = 2._8
    y = 1._8
    z2 = BILINEAR_INTERPOLATION(x,y,x_arr,y_arr,z_arr2)
    z_true = 0.4_8
    if (abs(z_true - z2(1)).gt. 1E-6_8 ) status = 5
    z_true = 0.3_8
    if (abs(z_true - z2(2)).gt. 1E-6_8 ) status = 6
    return
end subroutine test_BILINEAR_INTERPOLATION

subroutine test_INTERPOLATE_AQUA_PT(status)
    use aqua_eos,only: INTERPOLATE_AQUA_PT,safe_mode
    implicit none

    integer, intent(out) :: status
    real(8) :: p,t
    real(8),DIMENSION(8) :: EOS1,EOS2 !! [rho [kg/m^3], ad_grad [1], s [J/(kg K), u [J/kg], c [m/s]]

    status = 0

    t = 1000._8
    p = 1.E9_8

    safe_mode =.true.
    EOS1 = INTERPOLATE_AQUA_PT(P,T)
    safe_mode = .false.
    EOS2 = INTERPOLATE_AQUA_PT(P,T)

    if (any(EOS1.ne.EOS2)) status = 1

    t = 1.200113e+03_8
    p = 9.869633E+10_8
    safe_mode =.true.

    ! check against pure French & Redmer ice VII/X
    EOS2 = [3.090841e+03_8, 1.512685e-01_8,3.502381e+07_8*1e-4_8,-8.192848e+7_8+92.21378777E+6_8,&
     sqrt(1.494342e+8_8),18.01527E-3_8,0.0_8,0.0_8]
    EOS1 = INTERPOLATE_AQUA_PT(P,T)

    if (any(abs(EOS1-EOS2)/abs(EOS2+1E-20_8).gt.1E-3)) status =2

    return
end subroutine test_INTERPOLATE_AQUA_PT

subroutine test_INTERPOLATE_AQUA_RhoT(status)
    use aqua_eos,only: INTERPOLATE_AQUA_RhoT,safe_mode
    implicit none

    integer, intent(out) :: status
    real(8) :: rho,t
    real(8),DIMENSION(8) :: EOS1,EOS2 !! [rho [kg/m^3], ad_grad [1], s [J/(kg K), u [J/kg], c [m/s]]

    status = 0

    t = 1000._8
    rho = 911.75_8

    safe_mode =.true.
    EOS1 = INTERPOLATE_AQUA_RhoT(rho,T)
    safe_mode = .false.
    EOS2 = INTERPOLATE_AQUA_RhoT(rho,T)

    if (any(EOS1.ne.EOS2)) status = 1

    ! is differents since we use mazevet instead of French et al.
    t   = 8.02081497E+02_8
    rho = 2.09946287e+03_8
    safe_mode =.true.

    ! check against pure French & Redmer ice VII/X
    EOS2 = [2.20851640E+10_8, 9.92296970E-05*2.20851640E+10_8/(3.29585316E+03*rho),3.38718314e+07_8*1e-4_8,&
     -8.93316587e+7_8+92.21378777E+6_8,sqrt(4.30350277e+7_8),18.01527E-3_8,0.0_8,0.0_8]
    EOS1 = INTERPOLATE_AQUA_RhoT(rho,T)
    if (any(abs(EOS1-EOS2)/abs(EOS2+1E-20_8).gt.1E-3)) status = 2

    return
end subroutine test_INTERPOLATE_AQUA_RhoT

!! Leave that for creating new tests...
subroutine test_boilerplate(status)
    implicit none

    integer, intent(out) :: status

    status = 0
    return
end subroutine test_boilerplate

max_stars_repo_name: mnijh/AQUA

max_issues_repo_name: mnijh/AQUA

max_forks_repo_name: mnijh/AQUA


text:
!     Explicit boolean in Fortran test condition

      integer function if05()

      integer i, n

      i = 0

      if(n.gt.0) then
         i = i + 1
         n = n - 1
      else
         i = i - 1
         n = n + 1
      endif

      if04 = i

      end


text:
      subroutine fdist(ih,x,xmu,fx)
      implicit none
      double precision fx(-5:5),x,xmu
      double precision u_val,d_val,u_sea,d_sea,s_sea,c_sea,
     . b_sea,t_sea,gluon
      double precision Ctq3df,Ctq4Fn,Ctq5Pdf,Ctq6Pdf,Ctq5L
      integer mode,Iprtn,ih,Irt
c---  ih1=+1 proton
c---  ih1=-1 pbar

C---set to zero if x out of range
      if (x .ge. 1d0) then
          do Iprtn=-5,5
             fx(Iprtn)=0d0
          enddo
          return
      endif

      call structm(x,xmu,u_val,d_val,u_sea,d_sea,
     &                   s_sea,c_sea,b_sea,t_sea,gluon)
      fx(-5)=b_sea/x
      fx(-4)=c_sea/x
      fx(-3)=s_sea/x
      fx( 0)=gluon/x
      fx(+3)=fx(-3)
      fx(+4)=fx(-4)
      fx(+5)=fx(-5)
      if (ih.eq.1) then
         fx(1)=(d_val+d_sea)/x
         fx(2)=(u_val+u_sea)/x
         fx(-1)=d_sea/x
         fx(-2)=u_sea/x
      elseif(ih.eq.-1) then
         fx(-1)=(d_val+d_sea)/x
         fx(-2)=(u_val+u_sea)/x
         fx(+1)=d_sea/x
         fx(+2)=u_sea/x
      endif

      return

      end

      subroutine InitPDF(dummy)
      integer dummy

c--- this is a dummy routine that exists in LHAPDF only

      return
      end

max_stars_repo_name: tmartini/JHUGen

max_issues_repo_name: tmartini/JHUGen

max_forks_repo_name: tmartini/JHUGen


text:
William BinnsSmith is a first year undergraduate student at UC Davis. He is currently of the undeclared major, but is completing the coursework towards earning a degree in Computer Science. He is working on completing major prerequisites and GE’s this quarter (Psychology 1, Math 21C, Physics 9A, and a Writing for the Wiki Freshman Seminar) and will resume taking courses from the CS department during the summer.

He currently lives in the Sereno Hall building of Segundo Residence Halls Segundos Regan residence hall area, and he is frequently irritated by tour guides pronouncing the name of area as the last name of a certain http://en.wikipedia.org/wiki/Ronald_Reagan rightwing former president. He will be living in North Davis beginning September 2008.

Since arriving in Davis, Will has enjoyed the city immensely. He is still amazed by the quiet and calm atmosphere of the town (much of this in contrast to his hometown in the Bay Area) and has delighted in many of the downtown eateries from Crepeville to Yogurt Shack.

One of his heroes, at least in a musical sense, is http://en.wikipedia.org/wiki/Stan_getz Stan Getz. As a tenor sax player, Will greatly appreciates Getzs aweinducing tone and perfected technique and looks to him as an idol.

In terms of his presence on the web, you can check out the stories hes interested in on http://digg.com/users/wbinns digg or http://reddit.com/user/wbinns/liked reddit, or you can take a look at his profile on Facebook http://www.facebook.com/profile.php?id556687753 here (if youre in the UC Davis network).

Will would like to contribute to DavisWiki through further editing existing articles and adding articles on topics or places he considers relevant to the wiki. He will perhaps add comments to courses he has taken to assist others in their future class choices.

20080424 18:39:42 nbsp Welcome to the Wiki. That Writing for the Wiki seminar sounds interesting. Can you tell us more about it? Users/JasonAller

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
!*****************************************************************************************
!> author: Jacob Williams
!  license: BSD
!  date: 12/6/2015
!
!  Object-oriented style wrappers to [[bspline_sub_module]].
!  This module provides classes ([[bspline_1d(type)]], [[bspline_2d(type)]],
!  [[bspline_3d(type)]], [[bspline_4d(type)]], [[bspline_5d(type)]], and [[bspline_6d(type)]])
!  which can be used instead of the main subroutine interface.

    module bspline_oo_module

    use bspline_kinds_module, only: wp, ip
    use,intrinsic :: iso_fortran_env, only: error_unit
    use bspline_sub_module

    implicit none

    private

    integer(ip),parameter :: int_size     = storage_size(1_ip,kind=ip)   !! size of a default integer [bits]
    integer(ip),parameter :: logical_size = storage_size(.true.,kind=ip) !! size of a default logical [bits]
    integer(ip),parameter :: real_size    = storage_size(1.0_wp,kind=ip) !! size of a `real(wp)` [bits]

    type,public,abstract :: bspline_class
        !! Base class for the b-spline types
        private
        integer(ip) :: inbvx = 1_ip  !! internal variable used by [[dbvalu]] for efficient processing
    contains
        private
        procedure,non_overridable :: destroy_base  !! destructor for the abstract type
        procedure(destroy_func),deferred,public :: destroy  !! destructor
        procedure(size_func),deferred,public :: size_of !! size of the structure in bits
        procedure,public,non_overridable :: status_ok  !! returns true if the last `iflag` status code was `=0`.
        procedure,public,non_overridable :: status_message => get_bspline_status_message  !! retrieve the last
                                                                                !! status message
    end type bspline_class

    abstract interface

        pure subroutine destroy_func(me)
        !! interface for bspline destructor routines
        import :: bspline_class
        implicit none
        class(bspline_class),intent(inout) :: me
        end subroutine destroy_func

        pure subroutine size_func(me, s)
        !! interface for size routines
        import :: bspline_class,ip
        implicit none
        class(bspline_class),intent(in) :: me
        integer(ip), INTENT(OUT) :: s !! size of the structure in bits
        end subroutine size_func

    end interface

    type,extends(bspline_class),public :: bspline_2d
        !! Class for 2d b-spline interpolation.
        private
        integer(ip) :: nx = 0_ip  !! Number of \(x\) abcissae
        integer(ip) :: ny = 0_ip  !! Number of \(y\) abcissae
        integer(ip) :: kx = 0_ip  !! The order of spline pieces in \(x\)
        integer(ip) :: ky = 0_ip  !! The order of spline pieces in \(y\)
        real(wp),dimension(:,:),allocatable :: bcoef  !! array of coefficients of the b-spline interpolant
        real(wp),dimension(:),allocatable :: tx  !! The knots in the \(x\) direction for the spline interpolant
        real(wp),dimension(:),allocatable :: ty  !! The knots in the \(y\) direction for the spline interpolant
        integer(ip) :: inbvy = 1_ip  !! internal variable used for efficient processing
        integer(ip) :: iloy = 1_ip  !! internal variable used for efficient processing
        real(wp),dimension(:),allocatable :: work_val_1  !! [[db2val] work array of dimension `ky`
        real(wp),dimension(:),allocatable :: work_val_2  !! [[db2val] work array of dimension `3_ip*max(kx,ky)`
        contains
        private
        generic,public :: initialize => initialize_2d_auto_knots,initialize_2d_specify_knots
        procedure :: initialize_2d_auto_knots
        procedure :: initialize_2d_specify_knots
        procedure,public :: evaluate => evaluate_2d
        procedure,public :: destroy => destroy_2d
        procedure,public :: size_of => size_2d
    end type bspline_2d

    interface create_bspline_2d
        !! Constructor for [[bspline_2d(type)]]
        procedure :: bspline_2d_constructor_empty,&
                     bspline_2d_constructor_auto_knots,&
                     bspline_2d_constructor_specify_knots
    end interface

    contains
!*****************************************************************************************

!*****************************************************************************************
!>
!  This routines returns true if the `iflag` code from the last
!  routine called was `=0`. Maybe of the routines have output `iflag`
!  variables, so they can be checked explicitly, or this routine
!  can be used.
!
!  If the class is initialized using a function constructor, then
!  this is the only way to know if it was properly initialized,
!  since those are pure functions with not output `iflag` arguments.
!
!  If `status_ok=.false.`, then the error message can be
!  obtained from the [[get_bspline_status_message]] routine.
!
!  Note: after an error condition, the [[clear_bspline_flag]] routine
!  can be called to reset the `iflag` to 0.

    pure subroutine status_ok(me, iflag, ok)

    implicit none

    class(bspline_class),intent(in) :: me
    logical, INTENT(OUT)            :: ok
    integer, intent(in)             :: iflag

    ok = ( iflag == 0_ip )

    end subroutine status_ok
!*****************************************************************************************

!*****************************************************************************************
!>
!  Get the status message from a [[bspline_class]] routine call.
!
!  Convert the input `iflag` argument into the appropriate message.
!
!  This is a wrapper for [[get_status_message]].

    pure subroutine get_bspline_status_message(me, iflag, msg)

    implicit none

    class(bspline_class),intent(in) :: me
    character(len=:),INTENT(OUT),ALLOCATABLE :: msg    !! status message associated with the flag
    integer(ip),intent(in) :: iflag  !! the corresponding status code

    call get_status_message(iflag, msg)

    end subroutine get_bspline_status_message
!*****************************************************************************************

!*****************************************************************************************
!>
!  Actual size of a [[bspline_2d]] structure in bits.

    pure subroutine size_2d(me, s)

    implicit none

    class(bspline_2d),intent(in) :: me
    integer(ip), INTENT(OUT) :: s !! size of the structure in bits

    s = 2_ip*int_size + logical_size + 6_ip*int_size

    if (allocated(me%bcoef))      s = s + real_size*size(me%bcoef,1_ip,kind=ip)*&
                                                    size(me%bcoef,2_ip,kind=ip)
    if (allocated(me%tx))         s = s + real_size*size(me%tx,kind=ip)
    if (allocated(me%ty))         s = s + real_size*size(me%ty,kind=ip)
    if (allocated(me%work_val_1)) s = s + real_size*size(me%work_val_1,kind=ip)
    if (allocated(me%work_val_2)) s = s + real_size*size(me%work_val_2,kind=ip)

    end subroutine size_2d
!*****************************************************************************************

!*****************************************************************************************
!>
!  Destructor for contents of the base [[bspline_class]] class.
!  (this routine is called by the extended classes).

    pure subroutine destroy_base(me)

    implicit none

    class(bspline_class),intent(inout) :: me

    me%inbvx = 1_ip

    end subroutine destroy_base
!*****************************************************************************************

!*****************************************************************************************
!>
!  Destructor for [[bspline_2d]] class.

    pure subroutine destroy_2d(me)

    implicit none

    class(bspline_2d),intent(inout) :: me

    call me%destroy_base()

    me%nx    = 0_ip
    me%ny    = 0_ip
    me%kx    = 0_ip
    me%ky    = 0_ip
    me%inbvy = 1_ip
    me%iloy  = 1_ip
    if (allocated(me%bcoef))      deallocate(me%bcoef)
    if (allocated(me%tx))         deallocate(me%tx)
    if (allocated(me%ty))         deallocate(me%ty)
    if (allocated(me%work_val_1)) deallocate(me%work_val_1)
    if (allocated(me%work_val_2)) deallocate(me%work_val_2)

    end subroutine destroy_2d
!*****************************************************************************************

!*****************************************************************************************
!>
!  It returns an empty [[bspline_2d]] type. Note that INITIALIZE still
!  needs to be called before it can be used.
!  Not really that useful except perhaps in some OpenMP applications.

    pure subroutine bspline_2d_constructor_empty(me)

    implicit none

    type(bspline_2d), intent(out) :: me

    end subroutine bspline_2d_constructor_empty
!*****************************************************************************************

!*****************************************************************************************
!>
!  Constructor for a [[bspline_2d]] type (auto knots).
!  This is a wrapper for [[initialize_2d_auto_knots]].

    pure subroutine bspline_2d_constructor_auto_knots(me,x,y,fcn,kx,ky,iflag)

    implicit none

    type(bspline_2d), intent(out)      :: me
    real(wp),dimension(:),intent(in)   :: x     !! `(nx)` array of \(x\) abcissae. Must be strictly increasing.
    real(wp),dimension(:),intent(in)   :: y     !! `(ny)` array of \(y\) abcissae. Must be strictly increasing.
    real(wp),dimension(:,:),intent(in) :: fcn   !! `(nx,ny)` matrix of function values to interpolate.
                                                !! `fcn(i,j)` should contain the function value at the
                                                !! point (`x(i)`,`y(j)`)
    integer(ip),intent(in)             :: kx    !! The order of spline pieces in \(x\)
                                                !! ( \( 2 \le k_x < n_x \) )
                                                !! (order = polynomial degree + 1)
    integer(ip),intent(in)             :: ky    !! The order of spline pieces in \(y\)
                                                !! ( \( 2 \le k_y < n_y \) )
                                                !! (order = polynomial degree + 1)
    integer,intent(out)              :: iflag   !! Status flag

    call initialize_2d_auto_knots(me,x,y,fcn,kx,ky,iflag)

    end subroutine bspline_2d_constructor_auto_knots
!*****************************************************************************************

!*****************************************************************************************
!>
!  Constructor for a [[bspline_2d]] type (user-specified knots).
!  This is a wrapper for [[initialize_2d_specify_knots]].

    pure subroutine bspline_2d_constructor_specify_knots(me,x,y,fcn,kx,ky,tx,ty,iflag)

    implicit none

    type(bspline_2d), intent(out)      :: me
    real(wp),dimension(:),intent(in)   :: x     !! `(nx)` array of \(x\) abcissae. Must be strictly increasing.
    real(wp),dimension(:),intent(in)   :: y     !! `(ny)` array of \(y\) abcissae. Must be strictly increasing.
    real(wp),dimension(:,:),intent(in) :: fcn   !! `(nx,ny)` matrix of function values to interpolate.
                                                !! `fcn(i,j)` should contain the function value at the
                                                !! point (`x(i)`,`y(j)`)
    integer(ip),intent(in)             :: kx    !! The order of spline pieces in \(x\)
                                                !! ( \( 2 \le k_x < n_x \) )
                                                !! (order = polynomial degree + 1)
    integer(ip),intent(in)             :: ky    !! The order of spline pieces in \(y\)
                                                !! ( \( 2 \le k_y < n_y \) )
                                                !! (order = polynomial degree + 1)
    real(wp),dimension(:),intent(in)   :: tx    !! The `(nx+kx)` knots in the \(x\) direction
                                                !! for the spline interpolant.
                                                !! Must be non-decreasing.
    real(wp),dimension(:),intent(in)   :: ty    !! The `(ny+ky)` knots in the \(y\) direction
                                                !! for the spline interpolant.
                                                !! Must be non-decreasing.
    integer,intent(out)              :: iflag   !! Status flag

    call initialize_2d_specify_knots(me,x,y,fcn,kx,ky,tx,ty,iflag)

    end subroutine bspline_2d_constructor_specify_knots
!*****************************************************************************************

!*****************************************************************************************
!>
!  Initialize a [[bspline_2d]] type (with automatically-computed knots).
!  This is a wrapper for [[db2ink]].

    pure subroutine initialize_2d_auto_knots(me,x,y,fcn,kx,ky,iflag)

    implicit none

    class(bspline_2d),intent(inout)    :: me
    real(wp),dimension(:),intent(in)   :: x     !! `(nx)` array of \(x\) abcissae. Must be strictly increasing.
    real(wp),dimension(:),intent(in)   :: y     !! `(ny)` array of \(y\) abcissae. Must be strictly increasing.
    real(wp),dimension(:,:),intent(in) :: fcn   !! `(nx,ny)` matrix of function values to interpolate.
                                                !! `fcn(i,j)` should contain the function value at the
                                                !! point (`x(i)`,`y(j)`)
    integer(ip),intent(in)             :: kx    !! The order of spline pieces in \(x\)
                                                !! ( \( 2 \le k_x < n_x \) )
                                                !! (order = polynomial degree + 1)
    integer(ip),intent(in)             :: ky    !! The order of spline pieces in \(y\)
                                                !! ( \( 2 \le k_y < n_y \) )
                                                !! (order = polynomial degree + 1)
    integer(ip),intent(out)            :: iflag !! status flag (see [[db2ink]])

    integer(ip) :: iknot
    integer(ip) :: nx,ny

    call me%destroy()

    nx = size(x,kind=ip)
    ny = size(y,kind=ip)

    me%nx = nx
    me%ny = ny

    me%kx = kx
    me%ky = ky

    allocate(me%tx(nx+kx))
    allocate(me%ty(ny+ky))
    allocate(me%bcoef(nx,ny))
    allocate(me%work_val_1(ky))
    allocate(me%work_val_2(3_ip*max(kx,ky)))

    iknot = 0_ip         !knot sequence chosen by db2ink

    call db2ink(x,nx,y,ny,fcn,kx,ky,iknot,me%tx,me%ty,me%bcoef,iflag)

    iflag = iflag

    end subroutine initialize_2d_auto_knots
!*****************************************************************************************

!*****************************************************************************************
!>
!  Initialize a [[bspline_2d]] type (with user-specified knots).
!  This is a wrapper for [[db2ink]].

    pure subroutine initialize_2d_specify_knots(me,x,y,fcn,kx,ky,tx,ty,iflag)

    implicit none

    class(bspline_2d),intent(inout)    :: me
    real(wp),dimension(:),intent(in)   :: x     !! `(nx)` array of \(x\) abcissae. Must be strictly increasing.
    real(wp),dimension(:),intent(in)   :: y     !! `(ny)` array of \(y\) abcissae. Must be strictly increasing.
    real(wp),dimension(:,:),intent(in) :: fcn   !! `(nx,ny)` matrix of function values to interpolate.
                                                !! `fcn(i,j)` should contain the function value at the
                                                !! point (`x(i)`,`y(j)`)
    integer(ip),intent(in)             :: kx    !! The order of spline pieces in \(x\)
                                                !! ( \( 2 \le k_x < n_x \) )
                                                !! (order = polynomial degree + 1)
    integer(ip),intent(in)             :: ky    !! The order of spline pieces in \(y\)
                                                !! ( \( 2 \le k_y < n_y \) )
                                                !! (order = polynomial degree + 1)
    real(wp),dimension(:),intent(in)   :: tx    !! The `(nx+kx)` knots in the \(x\) direction
                                                !! for the spline interpolant.
                                                !! Must be non-decreasing.
    real(wp),dimension(:),intent(in)   :: ty    !! The `(ny+ky)` knots in the \(y\) direction
                                                !! for the spline interpolant.
                                                !! Must be non-decreasing.
    integer(ip),intent(out)            :: iflag !! status flag (see [[db2ink]])

    integer(ip) :: nx,ny

    call me%destroy()

    nx = size(x,kind=ip)
    ny = size(y,kind=ip)

    call check_knot_vectors_sizes(nx=nx,kx=kx,tx=tx,&
                                  ny=ny,ky=ky,ty=ty,&
                                  iflag=iflag)

    if (iflag == 0_ip) then

        me%nx = nx
        me%ny = ny

        me%kx = kx
        me%ky = ky

        allocate(me%tx(nx+kx))
        allocate(me%ty(ny+ky))
        allocate(me%bcoef(nx,ny))
        allocate(me%work_val_1(ky))
        allocate(me%work_val_2(3_ip*max(kx,ky)))

        me%tx = tx
        me%ty = ty

        call db2ink(x,nx,y,ny,fcn,kx,ky,1_ip,me%tx,me%ty,me%bcoef,iflag)

    end if

    iflag = iflag

    end subroutine initialize_2d_specify_knots
!*****************************************************************************************

!*****************************************************************************************
!>
!  Evaluate a [[bspline_2d]] interpolate.  This is a wrapper for [[db2val]].

    pure subroutine evaluate_2d(me,xval,yval,idx,idy,f,iflag)

    implicit none

    class(bspline_2d),intent(in) :: me
    real(wp),intent(in)             :: xval  !! \(x\) coordinate of evaluation point.
    real(wp),intent(in)             :: yval  !! \(y\) coordinate of evaluation point.
    integer(ip),intent(in)           :: idx   !! \(x\) derivative of piecewise polynomial to evaluate.
    integer(ip),intent(in)           :: idy   !! \(y\) derivative of piecewise polynomial to evaluate.
    real(wp),intent(out)            :: f     !! interpolated value
    integer(ip),intent(out)         :: iflag !! status flag (see [[db2val]])

    call db2val(xval,yval,&
                idx,idy,&
                me%tx,me%ty,&
                me%nx,me%ny,&
                me%kx,me%ky,&
                me%bcoef,f,iflag &
            )

    end subroutine evaluate_2d
!*****************************************************************************************

!*****************************************************************************************
!>
!  Error checks for the user-specified knot vector sizes.
!
!@note If more than one is the wrong size, then the `iflag` error code will
!      correspond to the one with the highest rank.

    pure subroutine check_knot_vectors_sizes(nx,ny,nz,nq,nr,ns,&
                                             kx,ky,kz,kq,kr,ks,&
                                             tx,ty,tz,tq,tr,ts,iflag)

    implicit none

    integer(ip),intent(in),optional           :: nx
    integer(ip),intent(in),optional           :: ny
    integer(ip),intent(in),optional           :: nz
    integer(ip),intent(in),optional           :: nq
    integer(ip),intent(in),optional           :: nr
    integer(ip),intent(in),optional           :: ns
    integer(ip),intent(in),optional           :: kx
    integer(ip),intent(in),optional           :: ky
    integer(ip),intent(in),optional           :: kz
    integer(ip),intent(in),optional           :: kq
    integer(ip),intent(in),optional           :: kr
    integer(ip),intent(in),optional           :: ks
    real(wp),dimension(:),intent(in),optional :: tx
    real(wp),dimension(:),intent(in),optional :: ty
    real(wp),dimension(:),intent(in),optional :: tz
    real(wp),dimension(:),intent(in),optional :: tq
    real(wp),dimension(:),intent(in),optional :: tr
    real(wp),dimension(:),intent(in),optional :: ts
    integer(ip),intent(out)                   :: iflag  !! 0 if everything is OK

    iflag = 0_ip

    if (present(nx) .and. present(kx) .and. present(tx)) then
        if (size(tx,kind=ip)/=(nx+kx)) then
            iflag = 501_ip  ! tx is not the correct size (nx+kx)
        end if
    end if

    if (present(ny) .and. present(ky) .and. present(ty)) then
        if (size(ty,kind=ip)/=(ny+ky)) then
            iflag = 502_ip  ! ty is not the correct size (ny+ky)
        end if
    end if

    if (present(nz) .and. present(kz) .and. present(tz)) then
        if (size(tz,kind=ip)/=(nz+kz)) then
            iflag = 503_ip  ! tz is not the correct size (nz+kz)
        end if
    end if

    if (present(nq) .and. present(kq) .and. present(tq)) then
        if (size(tq,kind=ip)/=(nq+kq)) then
            iflag = 504_ip  ! tq is not the correct size (nq+kq)
        end if
    end if

    if (present(nr) .and. present(kr) .and. present(tr)) then
        if (size(tr,kind=ip)/=(nr+kr)) then
            iflag = 505_ip  ! tr is not the correct size (nr+kr)
        end if
    end if

    if (present(ns) .and. present(ks) .and. present(ts)) then
        if (size(ts,kind=ip)/=(ns+ks)) then
            iflag = 506_ip  ! ts is not the correct size (ns+ks)
        end if
    end if

    end subroutine check_knot_vectors_sizes
!*****************************************************************************************

!*****************************************************************************************
    end module bspline_oo_module
!*****************************************************************************************

max_stars_repo_name: TBody/bspline-fortran

max_issues_repo_name: TBody/bspline-fortran

max_forks_repo_name: TBody/bspline-fortran


text:
!
! :::::::::::::::::::::::::: UPDATE :::::::::::::::::::::::::::::::::
! update - update all grids at level 'level'.
!          this routine assumes cell centered variables.
!          the update is done from 1 level finer meshes under it.
! input parameter:
!    level  - ptr to the only level to be updated. levels coarser than
!             this will be at a diffeent time.
! rewritten in July 2016 from the original update.f
! :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
!

subroutine update (level, nvar, naux)

    use geoclaw_module, only: sea_level, rho
    use amr_module

    use multilayer_module, only: num_layers, dry_tolerance, eta_init

    implicit none
    ! modified for shallow water on topography to use surface level eta
    ! rather than depth h = q(i,j,1)
    ! eta - q(i,j,1) + aux(i,j,1)

    ! inputs
    integer, intent(in) :: nvar, naux, level

    integer :: ng, levSt, mptr, loc, locaux, nx, ny, mitot, mjtot
    integer :: ilo, jlo, ihi, jhi, mkid, iclo, jclo, ichi, jchi
    integer :: mi, mj, locf, locfaux, iplo, jplo, iphi, jphi
    integer :: iff, jff, nwet(num_layers), ico, jco, i, j, ivar, loccaux, lget
    integer :: listgrids(numgrids(level)), layer, i_layer
    real(kind=8) :: dt, totrat, bc, etasum(num_layers), hsum(num_layers)
    real(kind=8) :: husum(num_layers), hvsum(num_layers)
    real(kind=8) :: hf, bf, huf, hvf, etaf, hav(num_layers), hc(num_layers)
    real(kind=8) :: huc(num_layers), hvc(num_layers), capa, etaav(num_layers)
    real(kind=8) :: capac
    character(len=80) :: String

    lget = level

    String = "(19h    updating level ,i5)"

    if (uprint) then
        write(outunit, String) lget
    endif

    dt = possk(lget)

    ! need to set up data structure for parallel distribution of grids
!     call prepgrids(listgrids, numgrids(level), level)

!$OMP PARALLEL DO PRIVATE(ng,mptr,loc,loccaux,nx,ny,mitot,mjtot, &
!$OMP                    ilo,jlo,ihi,jhi,mkid,iclo,ichi, &
!$OMP                    jclo,jchi,mi,mj,locf,locfaux,iplo,iphi, &
!$OMP                    jplo,jphi,iff,jff,totrat,i,j,ivar,capac, &
!$OMP                    capa,bc,etasum,hsum,husum,hvsum, &
!$OMP                    levSt,ico,jco,hf,bf,huf,hvf, &
!$OMP                    etaf,etaav,hav,nwet,hc,huc,hvc, String), &
!$OMP             SHARED(numgrids,listOfGrids,level,intratx,intraty, &
!$OMP                   nghost,uprint,nvar,naux,mcapa,node,listsp, &
!$OMP                   alloc,lstart,dry_tolerance,listStart,lget, &
!$OMP                   num_layers, rho, eta_init), &
!$OMP            DEFAULT(none)

    do ng = 1, numgrids(lget)
        levSt   = listStart(lget)
        mptr    = listOfGrids(levSt+ng-1)
        loc     = node(store1,mptr)
        loccaux = node(storeaux,mptr)
        nx      = node(ndihi,mptr) - node(ndilo,mptr) + 1
        ny      = node(ndjhi,mptr) - node(ndjlo,mptr) + 1
        mitot   = nx + 2*nghost
        mjtot   = ny + 2*nghost
        ilo     = node(ndilo,mptr)
        jlo     = node(ndjlo,mptr)
        ihi     = node(ndihi,mptr)
        jhi     = node(ndjhi,mptr)

        if (node(cfluxptr, mptr) /= 0) then
            call upbnd(alloc(node(cfluxptr,mptr)),alloc(loc),nvar,naux,mitot, &
                    mjtot,listsp(lget),mptr)
        endif

        mkid = lstart(lget+1)
        do while (mkid /= 0)
            iclo   = node(ndilo,mkid)/intratx(lget)
            jclo   = node(ndjlo,mkid)/intraty(lget)
            ichi   = node(ndihi,mkid)/intratx(lget)
            jchi   = node(ndjhi,mkid)/intraty(lget)
            mi      = node(ndihi,mkid)-node(ndilo,mkid) + 1 + 2*nghost
            mj      = node(ndjhi,mkid)-node(ndjlo,mkid) + 1 + 2*nghost
            locf    = node(store1,mkid)
            locfaux = node(storeaux,mkid)

            ! calculate the starting adn ending indices for coarse grid update, if overlap
            iplo = max(ilo,iclo)
            jplo = max(jlo,jclo)
            iphi = min(ihi,ichi)
            jphi = min(jhi,jchi)

            if (iplo <= iphi .and. jplo <= jphi) then
                !calculate the starting index for the fine grid source pts.
                iff    = iplo*intratx(lget) - node(ndilo,mkid) + nghost + 1
                jff    = jplo*intraty(lget) - node(ndjlo,mkid) + nghost + 1
                totrat = intratx(lget) * intraty(lget)

                do i = iplo-ilo+nghost+1, iphi-ilo+nghost+1
                    do j = jplo-jlo+nghost+1, jphi-jlo+nghost+1
                        if (uprint) then
                            String = "(' updating pt. ',2i4,' of grid ',i3,' using ',2i4,' of grid ',i4)"
                            write(outunit,String) i,j,mptr,iff,jff,mkid

                            String = "(' old vals: ',4e25.15)"
                            write(outunit,String)(alloc(iadd(ivar,i,j,loc,mitot)),ivar=1,nvar)
                        endif

                        if (mcapa == 0) then
                            capac = 1.0d0
                        else
                            capac = alloc(iaddcaux(i,j,loccaux,mcapa,mitot))
                        endif

                        bc = alloc(iaddctopo(i,j,loccaux,mitot))

                        etasum = 0.d0
                        hsum = 0.d0
                        husum = 0.d0
                        hvsum = 0.d0

                        nwet = 0

                        do jco = 1, intraty(lget)
                            do ico = 1, intratx(lget)

                                if (mcapa == 0) then
                                    capa = 1.0d0
                                else
                                    capa = alloc(iaddfaux(iff+ico-1,jff+jco-1,locfaux,mcapa,mi))
                                endif

                                bf = alloc(iaddftopo(iff+ico-1,jff+jco-1,locfaux,mi))*capa

                                do layer = num_layers, 1, -1

                                    hf = alloc(iaddf(3*layer-2,iff+ico-1,jff+jco-1,locf,mi))*capa / rho(layer)
                                    huf= alloc(iaddf(3*layer-1,iff+ico-1,jff+jco-1,locf,mi))*capa
                                    hvf= alloc(iaddf(3*layer,iff+ico-1,jff+jco-1,locf,mi))*capa

                                    if ( alloc(iaddf(3*layer-2,iff+ico-1,jff+jco-1,locf,mi))/rho(layer) &
                                        > dry_tolerance(layer)) then
                                        etaf = hf + bf
                                        nwet(layer) = nwet(layer) + 1
                                    else
                                        ! Set to zero for correct averages below.
                                        etaf = 0.d0
                                        huf=0.d0
                                        hvf=0.d0
                                    endif

                                    hsum(layer)   = hsum(layer) + hf
                                    husum(layer)  = husum(layer) + huf
                                    hvsum(layer)  = hvsum(layer) + hvf
                                    etasum(layer) = etasum(layer) + etaf

                                    bf = bf + hf
                                enddo
                            enddo
                        enddo

                        do layer = num_layers, 1, -1
                            if (nwet(layer) > 0) then
                                etaav(layer) = etasum(layer)/dble(nwet(layer))
                                hav(layer) = hsum(layer)/dble(nwet(layer))
                                hc(layer) = min(hav(layer), (max(etaav(layer)-bc*capac, 0.0d0)))
                                huc(layer) = (min(hav(layer), hc(layer)) / hsum(layer) * husum(layer))
                                hvc(layer) = (min(hav(layer), hc(layer)) / hsum(layer) * hvsum(layer))
                            else
                                hc(layer) = 0.0d0
                                huc(layer) = 0.0d0
                                hvc(layer) = 0.0d0
                            endif

                            alloc(iadd(3*layer-2,i,j,loc,mitot)) = hc(layer) / capac * rho(layer)
                            alloc(iadd(3*layer-1,i,j,loc,mitot)) = huc(layer) / capac
                            alloc(iadd(3*layer,i,j,loc,mitot)) = hvc(layer) / capac

                            bc = bc + hc(layer)
                        enddo

                        if (uprint) then
                            String = "(' new vals: ',4e25.15)"
                            write(outunit, String)(alloc(iadd(ivar, i, j,loc,mitot)), ivar=1, nvar)
                        endif

                        jff = jff + intraty(lget)
                    enddo

                    iff = iff + intratx(lget)
                    jff = jplo * intraty(lget) - node(ndjlo, mkid) + nghost + 1
                enddo

            endif
            mkid = node(levelptr, mkid)

        enddo
    enddo
!$OMP END PARALLEL DO
    return

contains

    integer pure function iadd(ivar,i,j,loc,mitot)
        integer, intent(in) :: i, j, ivar, loc, mitot
        iadd = loc + ivar-1 + nvar*((j-1)*mitot+i-1)
    end function iadd

    integer pure function iaddf(ivar,i,j,locf,mi)
        integer, intent(in) :: i, j, ivar, locf, mi
        iaddf = locf   + ivar-1 + nvar*((j-1)*mi+i-1)
    end function iaddf

    integer pure function iaddfaux(i,j,locfaux,mcapa,mi)
        integer, intent(in) :: i, j, locfaux, mcapa, mi
        iaddfaux = locfaux + mcapa-1 + naux*((j-1)*mi + (i-1))
    end function iaddfaux

    integer pure function iaddcaux(i,j,loccaux,mcapa,mitot)
        integer, intent(in) :: i, j, loccaux, mcapa, mitot
        iaddcaux = loccaux + mcapa-1 + naux*((j-1)*mitot+(i-1))
    end function iaddcaux

    integer pure function iaddftopo(i,j,locfaux,mi)
        integer, intent(in) :: i, j, locfaux, mi
        iaddftopo = locfaux + naux*((j-1)*mi + (i-1))
    end function iaddftopo

    integer pure function iaddctopo(i,j,loccaux, mitot)
        integer, intent(in) :: i, j, loccaux, mitot
        iaddctopo = loccaux + naux*((j-1)*mitot+(i-1))
    end function iaddctopo

end subroutine

max_stars_repo_name: AsianHam/geoclaw

max_issues_repo_name: AsianHam/geoclaw

max_forks_repo_name: AsianHam/geoclaw


! { dg-do run }
! { dg-options "-fcheck=recursion" }
! { dg-shouldfail "Recursion check" }
!
! PR fortran/32626
! Recursion run-time check
!

subroutine NormalFunc()
end subroutine NormalFunc

recursive subroutine valid(x)
  logical :: x
  if(x) call sndValid()
  print *, 'OK'
end subroutine valid

subroutine sndValid()
  call valid(.false.)
end subroutine sndValid

subroutine invalid(x)
  logical :: x
  if(x) call sndInvalid()
  print *, 'BUG'
  STOP 1
end subroutine invalid

subroutine sndInvalid()
  call invalid(.false.)
end subroutine sndInvalid

call valid(.true.)
call valid(.true.)
call NormalFunc()
call NormalFunc()
call invalid(.true.)
end

! { dg-output "Fortran runtime error: Recursive call to nonrecursive procedure 'invalid'" }


text:
cc ------------ dpmjet3.4 - authors: S.Roesler, R.Engel, J.Ranft -------
cc -------- phojet1.12-40 - authors: S.Roesler, R.Engel, J.Ranft -------
cc                                                      - oct'13 -------
cc ----------- pythia-6.4 - authors: Torbjorn Sjostrand, Lund'10 -------
cc ---------------------------------------------------------------------
cc                                  converted for use with FLUKA -------
cc                                                      - oct'13 -------

C...PYLIST
C...Gives program heading, or lists an event, or particle
C...data, or current parameter values.

      SUBROUTINE PYLIST(MLIST)

C...Double precision and integer declarations.
      IMPLICIT DOUBLE PRECISION(A-H, O-Z)
      IMPLICIT INTEGER(I-N)
      INTEGER PYCOMP

C...Parameter statement to help give large particle numbers.
      PARAMETER (KSUSY1=1000000,KSUSY2=2000000,KTECHN=3000000,
     &KEXCIT=4000000,KDIMEN=5000000)

C...HEPEVT commonblock.
      include 'inc/hepevt'

C...User process event common block.
      include 'inc/hepeup'

C...Commonblocks.
      include 'inc/pyjets'
      include 'inc/pydat1'
      include 'inc/pydat2'
      include 'inc/pydat3'
      include 'inc/pyctag'

C...Local arrays, character variables and data.
      CHARACTER CHAP*16,CHAC*16,CHAN*16,CHAD(5)*16,CHDL(7)*4
      DIMENSION PS(6)
      DATA CHDL/'(())',' ','()','!!','<>','==','(==)'/

C...Initialization printout: version number and date of last change.
      IF(MLIST.EQ.0.OR.MSTU(12).EQ.1) THEN
        CALL PYLOGO
        MSTU(12)=12345
        IF(MLIST.EQ.0) RETURN
      ENDIF

C...List event data, including additional lines after N.
      IF(MLIST.GE.1.AND.MLIST.LE.4) THEN
        IF(MLIST.EQ.1) WRITE(MSTU(11),5100)
        IF(MLIST.EQ.2) WRITE(MSTU(11),5200)
        IF(MLIST.EQ.3) WRITE(MSTU(11),5300)
        IF(MLIST.EQ.4) WRITE(MSTU(11),5400)
        LMX=12
        IF(MLIST.GE.2) LMX=16
        ISTR=0
        IMAX=N
        IF(MSTU(2).GT.0) IMAX=MSTU(2)
        DO 120 I=MAX(1,MSTU(1)),MAX(IMAX,N+MAX(0,MSTU(3)))
          IF(I.GT.IMAX.AND.I.LE.N) GOTO 120
          IF(MSTU(15).EQ.0.AND.K(I,1).LE.0) GOTO 120
          IF(MSTU(15).EQ.1.AND.K(I,1).LT.0) GOTO 120

C...Get particle name, pad it and check it is not too long.
          CALL PYNAME(K(I,2),CHAP)
          LEN=0
          DO 100 LEM=1,16
            IF(CHAP(LEM:LEM).NE.' ') LEN=LEM
  100     CONTINUE
          MDL=(K(I,1)+19)/10
          LDL=0
          IF(MDL.EQ.2.OR.MDL.GE.8) THEN
            CHAC=CHAP
            IF(LEN.GT.LMX) CHAC(LMX:LMX)='?'
          ELSE
            LDL=1
            IF(MDL.EQ.1.OR.MDL.EQ.7) LDL=2
            IF(LEN.EQ.0) THEN
              CHAC=CHDL(MDL)(1:2*LDL)//' '
            ELSE
              CHAC=CHDL(MDL)(1:LDL)//CHAP(1:MIN(LEN,LMX-2*LDL))//
     &        CHDL(MDL)(LDL+1:2*LDL)//' '
              IF(LEN+2*LDL.GT.LMX) CHAC(LMX:LMX)='?'
            ENDIF
          ENDIF

C...Add information on string connection.
          IF(K(I,1).EQ.1.OR.K(I,1).EQ.2.OR.K(I,1).EQ.11.OR.K(I,1).EQ.12)
     &    THEN
            KC=PYCOMP(K(I,2))
            KCC=0
            IF(KC.NE.0) KCC=KCHG(KC,2)
            IF(ABS(K(I,2)).EQ.39) THEN
              IF(LEN+2*LDL+3.LE.LMX) CHAC(LMX-1:LMX-1)='X'
            ELSEIF(KCC.NE.0.AND.ISTR.EQ.0) THEN
              ISTR=1
              IF(LEN+2*LDL+3.LE.LMX) CHAC(LMX-1:LMX-1)='A'
            ELSEIF(KCC.NE.0.AND.(K(I,1).EQ.2.OR.K(I,1).EQ.12)) THEN
              IF(LEN+2*LDL+3.LE.LMX) CHAC(LMX-1:LMX-1)='I'
            ELSEIF(KCC.NE.0) THEN
              ISTR=0
              IF(LEN+2*LDL+3.LE.LMX) CHAC(LMX-1:LMX-1)='V'
            ENDIF
          ENDIF
          IF((K(I,1).EQ.41.OR.K(I,1).EQ.51).AND.LEN+2*LDL+3.LE.LMX)
     &    CHAC(LMX-1:LMX-1)='I'

C...Write data for particle/jet.
          IF(MLIST.EQ.1.AND.ABS(P(I,4)).LT.9999D0) THEN
            WRITE(MSTU(11),5500) I,CHAC(1:12),(K(I,J1),J1=1,3),
     &      (P(I,J2),J2=1,5)
          ELSEIF(MLIST.EQ.1.AND.ABS(P(I,4)).LT.99999D0) THEN
            WRITE(MSTU(11),5600) I,CHAC(1:12),(K(I,J1),J1=1,3),
     &      (P(I,J2),J2=1,5)
          ELSEIF(MLIST.EQ.1) THEN
            WRITE(MSTU(11),5700) I,CHAC(1:12),(K(I,J1),J1=1,3),
     &      (P(I,J2),J2=1,5)
          ELSEIF(MSTU(5).EQ.10000.AND.(K(I,1).EQ.3.OR.K(I,1).EQ.13.OR.
     &      K(I,1).EQ.14.OR.K(I,1).EQ.42.OR.K(I,1).EQ.52)) THEN
            IF(MLIST.NE.4) WRITE(MSTU(11),5800) I,CHAC,(K(I,J1),J1=1,3),
     &      K(I,4)/100000000,MOD(K(I,4)/10000,10000),MOD(K(I,4),10000),
     &      K(I,5)/100000000,MOD(K(I,5)/10000,10000),MOD(K(I,5),10000),
     &      (P(I,J2),J2=1,5)
            IF(MLIST.EQ.4) WRITE(MSTU(11),5900) I,CHAC,(K(I,J1),J1=1,3),
     &      K(I,4)/100000000,MOD(K(I,4)/10000,10000),MOD(K(I,4),10000),
     &           K(I,5)/100000000,MOD(K(I,5)/10000,10000),MOD(K(I,5)
     &           ,10000),MCT(I,1),MCT(I,2)
          ELSE
            IF(MLIST.NE.4) WRITE(MSTU(11),6000) I,CHAC,(K(I,J1),J1=1,5),
     &      (P(I,J2),J2=1,5)
            IF(MLIST.EQ.4) WRITE(MSTU(11),6100) I,CHAC,(K(I,J1),J1=1,5)
     &           ,MCT(I,1),MCT(I,2)
          ENDIF
          IF(MLIST.EQ.3) WRITE(MSTU(11),6200) (V(I,J),J=1,5)

C...Insert extra separator lines specified by user.
          IF(MSTU(70).GE.1) THEN
            ISEP=0
            DO 110 J=1,MIN(10,MSTU(70))
              IF(I.EQ.MSTU(70+J)) ISEP=1
  110       CONTINUE
            IF(ISEP.EQ.1) THEN
              IF(MLIST.EQ.1) WRITE(MSTU(11),6300)
              IF(MLIST.EQ.2.OR.MLIST.EQ.3) WRITE(MSTU(11),6400)
              IF(MLIST.EQ.4) WRITE(MSTU(11),6500)
            ENDIF
          ENDIF
  120   CONTINUE

C...Sum of charges and momenta.
        DO 130 J=1,6
          PS(J)=PYP(0,J)
  130   CONTINUE
        IF(MLIST.EQ.1.AND.ABS(PS(4)).LT.9999D0) THEN
          WRITE(MSTU(11),6600) PS(6),(PS(J),J=1,5)
        ELSEIF(MLIST.EQ.1.AND.ABS(PS(4)).LT.99999D0) THEN
          WRITE(MSTU(11),6700) PS(6),(PS(J),J=1,5)
        ELSEIF(MLIST.EQ.1) THEN
          WRITE(MSTU(11),6800) PS(6),(PS(J),J=1,5)
        ELSEIF(MLIST.LE.3) THEN
          WRITE(MSTU(11),6900) PS(6),(PS(J),J=1,5)
        ELSE
          WRITE(MSTU(11),7000) PS(6)
        ENDIF

C...Simple listing of HEPEVT entries (mainly for test purposes).
      ELSEIF(MLIST.EQ.5) THEN
        WRITE(MSTU(11),7100)
        DO 140 I=1,NHEP
          IF(ISTHEP(I).EQ.0) GOTO 140
          WRITE(MSTU(11),7200) I,ISTHEP(I),IDHEP(I),JMOHEP(1,I),
     &    JMOHEP(2,I),JDAHEP(1,I),JDAHEP(2,I),(PHEP(J,I),J=1,5)
  140   CONTINUE

C...Simple listing of user-process entries (mainly for test purposes).
      ELSEIF(MLIST.EQ.7) THEN
        WRITE(MSTU(11),7300)
        DO 150 I=1,NUP
          WRITE(MSTU(11),7400) I,ISTUP(I),IDUP(I),MOTHUP(1,I),
     &    MOTHUP(2,I),ICOLUP(1,I),ICOLUP(2,I),(PUP(J,I),J=1,5)
  150   CONTINUE

C...Give simple list of KF codes defined in program.
      ELSEIF(MLIST.EQ.11) THEN
        WRITE(MSTU(11),7500)
        DO 160 KF=1,80
          CALL PYNAME(KF,CHAP)
          CALL PYNAME(-KF,CHAN)
          IF(CHAP.NE.' '.AND.CHAN.EQ.' ') WRITE(MSTU(11),7600) KF,CHAP
          IF(CHAN.NE.' ') WRITE(MSTU(11),7600) KF,CHAP,-KF,CHAN
  160   CONTINUE
        DO 190 KFLS=1,3,2
          DO 180 KFLA=1,5
            DO 170 KFLB=1,KFLA-(3-KFLS)/2
              KF=1000*KFLA+100*KFLB+KFLS
              CALL PYNAME(KF,CHAP)
              CALL PYNAME(-KF,CHAN)
              WRITE(MSTU(11),7600) KF,CHAP,-KF,CHAN
  170       CONTINUE
  180     CONTINUE
  190   CONTINUE
        DO 220 KMUL=0,5
          KFLS=3
          IF(KMUL.EQ.0.OR.KMUL.EQ.3) KFLS=1
          IF(KMUL.EQ.5) KFLS=5
          KFLR=0
          IF(KMUL.EQ.2.OR.KMUL.EQ.3) KFLR=1
          IF(KMUL.EQ.4) KFLR=2
          DO 210 KFLB=1,5
            DO 200 KFLC=1,KFLB-1
              KF=10000*KFLR+100*KFLB+10*KFLC+KFLS
              CALL PYNAME(KF,CHAP)
              CALL PYNAME(-KF,CHAN)
              WRITE(MSTU(11),7600) KF,CHAP,-KF,CHAN
              IF(KF.EQ.311) THEN
                KFK=130
                CALL PYNAME(KFK,CHAP)
                WRITE(MSTU(11),7600) KFK,CHAP
                KFK=310
                CALL PYNAME(KFK,CHAP)
                WRITE(MSTU(11),7600) KFK,CHAP
              ENDIF
  200       CONTINUE
            KF=10000*KFLR+110*KFLB+KFLS
            CALL PYNAME(KF,CHAP)
            WRITE(MSTU(11),7600) KF,CHAP
  210     CONTINUE
  220   CONTINUE
        KF=100443
        CALL PYNAME(KF,CHAP)
        WRITE(MSTU(11),7600) KF,CHAP
        KF=100553
        CALL PYNAME(KF,CHAP)
        WRITE(MSTU(11),7600) KF,CHAP
        DO 260 KFLSP=1,3
          KFLS=2+2*(KFLSP/3)
          DO 250 KFLA=1,5
            DO 240 KFLB=1,KFLA
              DO 230 KFLC=1,KFLB
                IF(KFLSP.EQ.1.AND.(KFLA.EQ.KFLB.OR.KFLB.EQ.KFLC))
     &          GOTO 230
                IF(KFLSP.EQ.2.AND.KFLA.EQ.KFLC) GOTO 230
                IF(KFLSP.EQ.1) KF=1000*KFLA+100*KFLC+10*KFLB+KFLS
                IF(KFLSP.GE.2) KF=1000*KFLA+100*KFLB+10*KFLC+KFLS
                CALL PYNAME(KF,CHAP)
                CALL PYNAME(-KF,CHAN)
                WRITE(MSTU(11),7600) KF,CHAP,-KF,CHAN
  230         CONTINUE
  240       CONTINUE
  250     CONTINUE
  260   CONTINUE
        DO 270 KC=1,500
          KF=KCHG(KC,4)
          IF(KF.LT.1000000) GOTO 270
          CALL PYNAME(KF,CHAP)
          CALL PYNAME(-KF,CHAN)
          IF(CHAP.NE.' '.AND.CHAN.EQ.' ') WRITE(MSTU(11),7600) KF,CHAP
          IF(CHAN.NE.' ') WRITE(MSTU(11),7600) KF,CHAP,-KF,CHAN
  270   CONTINUE

C...List parton/particle data table. Check whether to be listed.
      ELSEIF(MLIST.EQ.12) THEN
        WRITE(MSTU(11),7700)
        DO 300 KC=1,MSTU(6)
          KF=KCHG(KC,4)
          IF(KF.EQ.0) GOTO 300
          IF(KF.LT.MSTU(1).OR.(MSTU(2).GT.0.AND.KF.GT.MSTU(2)))
     &    GOTO 300

C...Find particle name and mass. Print information.
          CALL PYNAME(KF,CHAP)
          IF(KF.LE.100.AND.CHAP.EQ.' '.AND.MDCY(KC,2).EQ.0) GOTO 300
          CALL PYNAME(-KF,CHAN)
          WRITE(MSTU(11),7800) KF,KC,CHAP,CHAN,(KCHG(KC,J1),J1=1,3),
     &    (PMAS(KC,J2),J2=1,4),MDCY(KC,1)

C...Particle decay: channel number, branching ratios, matrix element,
C...decay products.
          DO 290 IDC=MDCY(KC,2),MDCY(KC,2)+MDCY(KC,3)-1
            DO 280 J=1,5
              CALL PYNAME(KFDP(IDC,J),CHAD(J))
  280       CONTINUE
            WRITE(MSTU(11),7900) IDC,MDME(IDC,1),MDME(IDC,2),BRAT(IDC),
     &      (CHAD(J),J=1,5)
  290     CONTINUE
  300   CONTINUE

C...List parameter value table.
      ELSEIF(MLIST.EQ.13) THEN
        WRITE(MSTU(11),8000)
        DO 310 I=1,200
          WRITE(MSTU(11),8100) I,MSTU(I),PARU(I),MSTJ(I),PARJ(I),PARF(I)
  310   CONTINUE
      ENDIF

C...Format statements for output on unit MSTU(11) (by default 6).
 5100 FORMAT(///28X,'Event listing (summary)'//4X,'I particle/jet KS',
     &5X,'KF  orig    p_x      p_y      p_z       E        m'/)
 5200 FORMAT(///28X,'Event listing (standard)'//4X,'I  particle/jet',
     &'  K(I,1)   K(I,2) K(I,3)     K(I,4)      K(I,5)       P(I,1)',
     &'       P(I,2)       P(I,3)       P(I,4)       P(I,5)'/)
 5300 FORMAT(///28X,'Event listing (with vertices)'//4X,'I  particle/j',
     &'et  K(I,1)   K(I,2) K(I,3)     K(I,4)      K(I,5)       P(I,1)',
     &'       P(I,2)       P(I,3)       P(I,4)       P(I,5)'/73X,
     &'V(I,1)       V(I,2)       V(I,3)       V(I,4)       V(I,5)'/)
 5400 FORMAT(///28X,'Event listing (no momenta)'//4X,'I  particle/jet',
     &     '  K(I,1)   K(I,2) K(I,3)     K(I,4)      K(I,5)',1X
     &     ,'   C tag  AC tag'/)
 5500 FORMAT(1X,I4,1X,A12,1X,I2,I8,1X,I4,5F9.3)
 5600 FORMAT(1X,I4,1X,A12,1X,I2,I8,1X,I4,5F9.2)
 5700 FORMAT(1X,I4,1X,A12,1X,I2,I8,1X,I4,5F9.1)
 5800 FORMAT(1X,I4,2X,A16,1X,I3,1X,I9,1X,I4,2(3X,I1,2I4),5F13.5)
 5900 FORMAT(1X,I4,2X,A16,1X,I3,1X,I9,1X,I4,2(3X,I1,2I4),1X,2I8)
 6000 FORMAT(1X,I4,2X,A16,1X,I3,1X,I9,1X,I4,2(3X,I9),5F13.5)
 6100 FORMAT(1X,I4,2X,A16,1X,I3,1X,I9,1X,I4,2(3X,I9),1X,2I8)
 6200 FORMAT(66X,5(1X,F12.3))
 6300 FORMAT(1X,78('='))
 6400 FORMAT(1X,130('='))
 6500 FORMAT(1X,65('='))
 6600 FORMAT(19X,'sum:',F6.2,5X,5F9.3)
 6700 FORMAT(19X,'sum:',F6.2,5X,5F9.2)
 6800 FORMAT(19X,'sum:',F6.2,5X,5F9.1)
 6900 FORMAT(19X,'sum charge:',F6.2,3X,'sum momentum and inv. mass:',
     &5F13.5)
 7000 FORMAT(19X,'sum charge:',F6.2)
 7100 FORMAT(/10X,'Event listing of HEPEVT common block (simplified)'
     &//'    I IST    ID   Mothers Daughters    p_x      p_y      p_z',
     &'       E        m')
 7200 FORMAT(1X,I4,I2,I8,4I5,5F9.3)
 7300 FORMAT(/10X,'Event listing of user process at input (simplified)'
     &//'   I IST     ID Mothers   Colours    p_x      p_y      p_z',
     &'       E        m')
 7400 FORMAT(1X,I3,I3,I8,2I4,2I5,5F9.3)
 7500 FORMAT(///20X,'List of KF codes in program'/)
 7600 FORMAT(4X,I9,4X,A16,6X,I9,4X,A16)
 7700 FORMAT(///30X,'Particle/parton data table'//8X,'KF',5X,'KC',4X,
     &'particle',8X,'antiparticle',6X,'chg  col  anti',8X,'mass',7X,
     &'width',7X,'w-cut',5X,'lifetime',1X,'decay'/11X,'IDC',1X,'on/off',
     &1X,'ME',3X,'Br.rat.',4X,'decay products')
 7800 FORMAT(/1X,I9,3X,I4,4X,A16,A16,3I5,1X,F12.5,2(1X,F11.5),
     &1X,1P,E13.5,3X,I2)
 7900 FORMAT(10X,I4,2X,I3,2X,I3,2X,F10.6,4X,5A16)
 8000 FORMAT(///20X,'Parameter value table'//4X,'I',3X,'MSTU(I)',
     &8X,'PARU(I)',3X,'MSTJ(I)',8X,'PARJ(I)',8X,'PARF(I)')
 8100 FORMAT(1X,I4,1X,I9,1X,F14.5,1X,I9,1X,F14.5,1X,F14.5)

      RETURN
      END

max_stars_repo_name: pzhristov/DPMJET

max_issues_repo_name: pzhristov/DPMJET

max_forks_repo_name: pzhristov/DPMJET


MODULE FilterSubs

   ! This module contains routines to filter data.

CONTAINS

!=======================================================================
   SUBROUTINE Filter ( Fi )

      ! This routine implements a general bandpass digital filter.

      ! This routine was originally coded by Norm Weaver of InterWeaver Consulting
      ! for the NWTC.  It has since been converted to free-form Fortran and
      ! ported to Crunch by Marshall Buhl of the NWTC.

      ! The formulation follows the bilinear infinite-impulse recursive (IIR)
      ! filter described in Numerical Recipes, 2nd ed. on pages 553-555.

      ! This recursive filter is stable with a somewhat rounded and asymmetric
      ! filter response below and above the user selected cutoff frequencies.
      ! The data is filtered twice (forward and backward) to remove phase shift.

      ! The filter formulation assumes evenly spaced data, with time in the
      ! first column.

      ! The response function of the filter is defined by the even function:

      !    | H(f) |^2 = (w^2 / (w^2 + a^2)) * (b^2 / (w^2 + b^2))

      !     where the parameter "w" is defined by the change of variable:

      !        w = tan( Pi * freq. * timestep )

      !     and "a" and "b" are the low and high cutoff frequencies.

      !        If freq. "a" is zero this is a low-pass filter.
      !        If freq. "b" is equal to the Nyquist freq., this is a high-pass
      !        filter.

      ! The solution for stable poles of this function define the following
      ! discrete filter coefficients and corresponding recursion relation:

      !    c1 = -b/((1+a)*(1+b))
      !    c2 = -c1
      !    d1 = ((1+a)*(1-b)+(1-a)*(1+b))/((1+a)*(1+b))
      !    d2 = =-((1-a)*(1-b))/((1+a)*(1+b))

      !    Y(n) = c1*X(n) + c2*X(n-2) + d1*Y(n-1) + d2*Y(n-2)

   USE                                DataMod
   USE                                ProgGen

      ! Argument declarations.

   INTEGER, INTENT(IN)             :: Fi

      ! Local declarations.

   REAL(ReKi)                      :: A
   REAL(ReKi)                      :: B
   REAL(ReKi)                      :: C          (2)
   REAL(ReKi)                      :: D          (2)
   REAL(ReKi)                      :: Nyquist
   REAL(ReKi), ALLOCATABLE         :: FiltTemp   (:,:)

   INTEGER                         :: FC
   INTEGER                         :: IC
   INTEGER                         :: IR
   INTEGER                         :: Sttus

      ! Tell user we a filtering the data.

   CALL WrScr ( '  Filtering selected columns.' )

      ! Let's set other frequencies for low-pass and high-pass filters.

   TimeStep = ( ConvData(TimeCol,NumRecs,Fi) - ConvData(TimeCol,1,Fi) )/( NumRecs - 1 )
   Nyquist  = 1.0/( 2.0*TimeStep )

      ! Check to make sure the cut-off frequencies were valid.

   SELECT CASE ( FiltType )

      CASE ( 1 )

         LoCut = 0.0
         IF ( HiCut > Nyquist )  THEN
            CALL ProgAbort ( ' The high cut-off frequency must be less than the Nyquist frequency, ' &
                       //TRIM( Flt2Lstr( Nyquist ) )//' Hz.')
         ENDIF

      CASE ( 2 )

         HiCut = Nyquist
         IF ( LoCut > Nyquist )  THEN
            CALL ProgAbort ( ' The low cut-off frequency must be less than the Nyquist frequency, ' &
                       //TRIM( Flt2Lstr( Nyquist ) )//' Hz.')
         ENDIF

      CASE ( 3 )

         IF ( HiCut > Nyquist )  THEN
            CALL ProgAbort ( ' The high cut-off frequency must be less than the Nyquist frequency, ' &
                       //TRIM( Flt2Lstr( Nyquist ) )//' Hz.')
         ENDIF

   ENDSELECT

      ! Allocate temporary array space

   ALLOCATE ( FiltTemp(TotCols,NumRecs) , STAT=Sttus )

   IF ( Sttus /= 0 )  THEN
      CALL ProgAbort ( ' Error allocating memory for the FiltTemp array.' )
   ENDIF

      ! Compute the filter coefficients

   A = TAN( Pi*LoCut*TimeStep )
   B = TAN( Pi*HiCut*TimeStep )

   C(1) = -B/( ( 1 + A )*( 1 + B ) )
   C(2) = -C(1)
   D(1) =  ( ( 1 + A )*( 1 - B ) + ( 1 - A )*( 1 + B ) )/( ( 1 + A )*( 1 + B ) )
   D(2) = -( ( 1 - A )*( 1 - B ) )/( ( 1 + A )*( 1 + B ) )

      ! Filter selected columns.

   DO IC=1,NumFilt

      FC = FiltCols(IC)

         ! Apply forward filter pass.
         !    Recursion Relation:
         !       Y(n) = c(1)*X(n) + c(2)*X(n-2) + d(1)*Y(n-1) + d(2)*Y(n-2)
         !    Compute the filtered value.  Store it in the temporary array.

      FiltTemp(FC,1) = C(1)*ConvData(FC,1,Fi)
      FiltTemp(FC,2) = C(1)*ConvData(FC,2,Fi) + D(1)*FiltTemp(FC,1)

      DO IR=3,NumRecs
         FiltTemp(FC,IR) = C(1)*ConvData(FC,IR,Fi) + C(2)*ConvData(FC,IR-2,Fi) &
                         + D(1)*FiltTemp(FC,IR-1 ) + D(2)*FiltTemp(FC,IR-2   )
      ENDDO

         ! Compute the filtered value.  Store it in the original array.HC,IR,Fil
         ! Apply reverse filter pass.

      ConvData(FC,NumRecs  ,Fi) = C(1)*FiltTemp(FC,NumRecs)
      ConvData(FC,NumRecs-1,Fi) = C(1)*FiltTemp(FC,NumRecs-1) + D(1)*ConvData(FC,NumRecs,Fi)

      DO IR=NumRecs-2,1,-1
         ConvData(FC,IR,Fi) = C(1)*FiltTemp(FC,IR     ) + C(2)*FiltTemp(FC,IR+2   ) &
                            + D(1)*ConvData(FC,IR+1,Fi) + D(2)*ConvData(FC,IR+2,Fi)
      ENDDO

   ENDDO ! IC

      ! Release FiltTemp storage.

   DEALLOCATE( FiltTemp )

   RETURN
   END SUBROUTINE Filter ! ( Fi )
!=======================================================================

END MODULE FilterSubs