text:
c
c -------------------------------------------------------------------
c
      subroutine outval(val,maxip1,maxjp1,nvar,irr,mitot,mjtot,
     1                  mptr,outgrd,lstgrd)
c
      implicit double precision (a-h,o-z)
      dimension val(maxip1,maxjp1,nvar), irr(mitot,mjtot)
      dimension pout(20,2)
      logical    outgrd
      logical            graf
      parameter  (maxgr = 192, maxlv = 6)
      common  /cirr2/  rlink(9,2,2933)
      common /nodal/ rnode(12,maxgr),node(17,maxgr),lstart(maxlv),
     *  newstl(maxlv),
     *  listsp(maxlv),intrat(maxlv),tol,bzonex,bzoney,mstart,ndfree,
     *  lfine,iorder,mxnest,kcheck,lwidth,
     *  maxcl, graf, lhead
      common  /cirr/ poly(10,2,2933),hs(10,2,2933),ar(-1:2933),
     *               points(24,2,2933),wt(24,-1:2933),
     *               xcirr(2933),ycirr(2933),ixsave(2933),iysave(2933),
     *               nxtirr(2933)
      common   /userdt/  gamma,gamma1,xprob,yprob,iprob,ismp,cfl
      data    spval/1.e10/
c
c if cell outside domain, don't print soln. value - meaningless.
c account for border of lwidth cells around grid in irr array.
c
      if (.not. outgrd) go to 99
c
c     cyl quantites
      velu   = 0.d0
      velv   = 0.d0
      ubd    = 0.d0
      vbd    = 0.d0
      uflerr = 0.d0
      vflerr = 0.d0
      ufltot = 0.d0
      vfltot = 0.d0
c     boundary quantities
      wdenerr = 0.d0
      wden    = 0.d0
      pbd = 0.d0
      pr1 = 0.d0
      odenerr = 0.d0
      oden    = 0.d0
      opbd = 0.d0
      opr1 = 0.d0
c     flowfield quantities
      totden = 0.d0
      totpr  = 0.d0
      totderr = 0.d0
      totperr = 0.d0
      wrongden = 0.d0
      artot    = 0.d0

      ar(-1) = 0.d0
      hx     = rnode(9,mptr)
      hy     = rnode(10,mptr)
      xlow   = rnode(1,mptr) - lwidth*hx
      ylow   = rnode(2,mptr) - lwidth*hy
      ar(lstgrd) = rnode(9,mptr)*rnode(10,mptr)
      do 20 i=2,maxip1-1
      do 20 j=2,maxjp1-1
          x  = rnode(1,mptr) + hx*(dfloat(i)-1.5)
          y  = rnode(2,mptr) + hy*(dfloat(j)-1.5)
	  iadj = i + lwidth - 1
	  jadj = j + lwidth - 1
          if (irr(iadj,jadj) .eq. -1) then
c           solid body cell
	    err1 = spval
	    err2 = spval
	    err3 = spval
	    err4 = spval
            write(6,109) x,y,iadj,jadj,(val(i,j,ivar),ivar=1,nvar)
109         format(2h *,2hx=,f7.3,3h y=,f7.3,4h, i=,i3,4h, j=,i3,
     *             ' a = ',4(f9.4,1x))
	   else if (irr(iadj,jadj) .eq. lstgrd) then
c               regular cell
                write(6,107) x,y,iadj,jadj,(val(i,j,ivar),ivar=1,nvar)
	        if (iprob.ge.15 .and. iprob .le. 19) then
	           if (iprob .eq. 17) call p17tru(x,y,rho,u,v,p)
	           if (iprob .eq. 18) call p18tru(x,y,rho,u,v,p)
	           if (iprob .eq. 19) then
		      call makep(poly(1,1,lstgrd),iadj,jadj,xlow,ylow,
     *                           hx,hy)
                      rlink(2,1,lstgrd) = -11
c                     call preppoly(poly(1,1,lstgrd),rlink(1,1,lstgrd),
c    .                              pout)
c		      call p19tru(x,y,rho,u,v,p,pout,lstgrd)
 		      call p19tru(x,y,rho,u,v,p,poly(1,1,lstgrd),lstgrd)
		   endif
	           if (iprob .eq. 15) then
		      call makep(poly(1,1,lstgrd),iadj,jadj,xlow,ylow,
     *                           hx,hy)
 		      call p15tru(x,y,rho,u,v,p,poly(1,1,lstgrd))
		   endif
	           ru = rho*u
	           rv = rho*v
	           en = p/gamma1 + 0.5d0*rho*(u**2+v**2)
	           write(6,117) rho,ru,rv,en
	           err1 = val(i,j,1)-rho
	           err2 = val(i,j,2)-ru
	           err3 = val(i,j,3)-rv
	           err4 = val(i,j,4)-en
	           write(6,127) err1,err2,err3,err4
		   err2 = val(i,j,2)/val(i,j,1) - u
		   err3 = val(i,j,3)/val(i,j,1) - v
	         endif
107         format(2x,2hx=,f7.3,3h y=,f7.3,4h, i=,i3,4h, j=,i3,' a = ',
     *                     4(f9.4,1x))
117         format(34x,'true: ',4(f9.4,1x))
127         format(33x,'error: ',4(d9.3,1x))
122         format('.',2i3,2f6.3,4e12.4)
	      else
c                irregular cell
		 kirr = irr(iadj,jadj)
	         x    = xcirr(kirr)
	         y    = ycirr(kirr)
                 write(6,106) x,y,iadj,jadj,(val(i,j,ivar),ivar=1,nvar)
	         if (iprob.ge.15 .and. iprob .le. 19) then
	           if (iprob .eq. 17) call p17tru(x,y,rho,u,v,p)
	           if (iprob .eq. 18) call p18tru(x,y,rho,u,v,p)
	           if (iprob .eq. 15) then
 		      call p15tru(x,y,rho,u,v,p,poly(1,1,kirr))
		   endif
	           if (iprob .eq. 19) then
c                     call preppoly(poly(1,1,kirr),rlink(1,1,kirr),pout)
c		      call p19tru(x,y,rho,u,v,p,pout,kirr)
 		      call p19tru(x,y,rho,u,v,p,poly(1,1,kirr),kirr)
		   endif
	           ru = rho*u
	           rv = rho*v
	           en = p/gamma1 + 0.5d0*rho*(u**2+v**2)
	           write(6,117) rho,ru,rv,en
	           err1 = val(i,j,1)-rho
	           err2 = val(i,j,2)-ru
	           err3 = val(i,j,3)-rv
	           err4 = val(i,j,4)-en
	           write(6,127) err1,err2,err3,err4
		   err2 = val(i,j,2)/val(i,j,1) - u
		   err3 = val(i,j,3)/val(i,j,1) - v
		 endif
106        format(2h +,2hx=,f7.3,3h y=,f7.3,4h, i=,i3,4h, j=,i3,' a = ',
     *                      4(f9.4,1x))
123        format('+',2i3,2f6.3,4e12.4)
          endif
	  if (iprob .ge. 15 .and. iprob .le. 19) then
	     k = irr(iadj,jadj)
	     pcomp = gamma1*(val(i,j,4)-.5d0*(val(i,j,2)**2+
     .                      val(i,j,3)**2)/val(i,j,1))
	     perr = pcomp - p
	     if (k.ne.lstgrd .and. k.ne.-1) then
		 wdenerr = wdenerr + sqrt(ar(k))*dabs(err1)
		 wden    = wden +    sqrt(ar(k))*rho
		 pbd     = pbd +     sqrt(ar(k))*dabs(perr)
		 pr1     = pr1 +     sqrt(ar(k))*p
c                :: also use alternate weighting for boundary norm
		 do 21 kside = 1, 6
		    if (poly(kside+2,1,k).eq.-11) then
		       sidex2 = poly(kside,1,k)
		       sidey2 = poly(kside,2,k)
		       sidex1 = poly(kside+1,1,k)
		       sidey1 = poly(kside+1,2,k)
		       go to 22
		    endif
 21              continue
 22		 rlen = dsqrt((sidey1-sidey2)**2 + (sidex1-sidex2)**2)
		 odenerr = odenerr + rlen*dabs(err1)
		 oden    = oden +    rlen*rho
		 opbd     = opbd +   rlen*dabs(perr)
		 opr1     = opr1 +   rlen *p
c                cyl quantities
		 velu    = velu + sqrt(ar(k))*dabs(u)
		 velv    = velv + sqrt(ar(k))*dabs(v)
		 ubd     = ubd + sqrt(ar(k))*dabs(err2)
		 vbd     = vbd + sqrt(ar(k))*dabs(err3)
	     endif
	     totderr = totderr + ar(k)*dabs(err1)
	     totperr = totperr + ar(k)*dabs(perr)
	     totden = totden + ar(k)*rho
	     totpr  = totpr  + ar(k)*p
	     artot  = artot + ar(k)
	     if (k .ne. -1) then
	       wrongden = wrongden + ar(k)*dabs(err1)/rho
	     endif
c            cyl quantities
	     uflerr = uflerr + ar(k)*dabs(err2)
	     vflerr = vflerr + ar(k)*dabs(err3)
	     ufltot = ufltot + ar(k)*dabs(u)
	     vfltot = vfltot + ar(k)*dabs(v)
	 endif
 20   continue
c
 99   if (iprob .ge. 15 .and. iprob .le. 19)  then
	write(6,*)'  '
        write(6,*)' 1 norm of sboundary density error     = ', wdenerr
        write(6,*)' 1 norm of sboundary density           = ', wden
        write(6,*)' 1 norm of spressure error             = ', pbd
        write(6,*)' 1 norm of spressure                   = ', pr1
	write(6,*)'  '
        write(6,*)' 1 norm of rboundary density error     = ', odenerr
        write(6,*)' 1 norm of rboundary density           = ', oden
        write(6,*)' 1 norm of rpressure error             = ', opbd
        write(6,*)' 1 norm of rpressure                   = ', opr1
	write(6,*)'  '
	write(6,*)' 1 norm total ff density error        = ',totderr
	write(6,*)' reg 1 norm of ff density             = ',totden
	write(6,*)' 1 norm total ff pressure error       = ',totperr
	write(6,*)' reg 1 norm of ff pressure            = ',totpr
	write(6,*)' wrong 1norm of ff density rel. err.  = ',wrongden
	write(6,*)' total area                           = ',artot

	if (iprob .eq.17) then
	  write(6,*)'  '
	  write(6,*)' 1 norm of u velocity bdry error      = ', ubd
	  write(6,*)' 1 norm of v velocity bdry error      = ', vbd
	  write(6,*)' 1 norm of u bdry velocity            = ', velu
	  write(6,*)' 1 norm of v bdry velocity            = ', velv

	  write(6,*)'  '
	  write(6,*)' 1 norm of u velocity ff error        = ', uflerr
	  write(6,*)' 1 norm of v velocity ff error        = ', vflerr
	  write(6,*)' 1 norm of u ff velocity              = ', ufltot
	  write(6,*)' 1 norm of v ff velocity              = ', vfltot
	endif
      endif

      return
      end


text:
program pario
  use mpi
  use, intrinsic :: iso_fortran_env, only : error_unit, output_unit
  implicit none

  integer, parameter :: datasize = 64, writer_id = 0
  integer :: rc, my_id, ntasks, localsize, i
  integer, dimension(:), allocatable :: localvector
  integer, dimension(datasize) :: fullvector

  call mpi_init(rc)
  call mpi_comm_size(mpi_comm_world, ntasks, rc)
  call mpi_comm_rank(mpi_comm_world, my_id, rc)

  if (ntasks > 64) then
    write(error_unit, *) 'Maximum number of tasks is 64!'
    call mpi_abort(MPI_COMM_WORLD, -1, rc)
  end if

  if (mod(datasize, ntasks) /= 0) then
    write(error_unit,*) 'Datasize (64) should be divisible by number of tasks'
    call mpi_abort(MPI_COMM_WORLD, -1, rc)
  end if

  localsize = datasize / ntasks
  allocate(localvector(localsize))

  localvector = [(i + my_id * localsize, i=1,localsize)]

  call single_writer()

  deallocate(localvector)
  call mpi_finalize(rc)

contains

  subroutine single_writer()
    implicit none

    ! TODO: Implement a function that writers the whole array of elements
    !       to a file so that single process is responsible for the file io

  end subroutine single_writer

end program pario

max_stars_repo_name: mrgransky/advanced-parallel-prog

max_issues_repo_name: mrgransky/advanced-parallel-prog

max_forks_repo_name: mrgransky/advanced-parallel-prog


text:
* X01FTEXT
*UPTODATE X01AAFTEXT
      DOUBLE PRECISION FUNCTION X01AAF(X)
C     NAG COPYRIGHT 1975
C     MARK 4.5 RELEASE
      DOUBLE PRECISION X
C     * PI *
C     RETURNS THE VALUE OF THE CONSTANT PI
C     THE X PARAMETER IS NOT USED
C     FOR DEC VAX 11/780
C     X01AAF = 3.14159265358979324D0
      X01AAF = 3.14159265358979324D0
      RETURN
      END
** END OF X01AAFTEXT
*UPTODATE X01ABFTEXT
      DOUBLE PRECISION FUNCTION X01ABF(X)
C     NAG COPYRIGHT 1975
C     MARK 4.5 RELEASE
      DOUBLE PRECISION X
C     * GAMMA *
C     RETURNS THE VALUE OF EULERS CONSTANT GAMMA WHERE
C     GAMMA = LIM(1+1/2+...+1/N - LOG(N)) AS N TENDS TO INFINITY
C     THE X PARAMETER IS NOT USED
C     FOR DEC VAX 11/780
C     X01ABF = 0.577215664901532861D0
      X01ABF = 0.577215664901532861D0
      RETURN
      END
** END OF X01ABFTEXT
* P01FTEXT
*UPTODATE P01AAFTEXT
      INTEGER FUNCTION P01AAF(IFAIL, ERROR, SRNAME)
C     MARK 1 RELEASE.  NAG COPYRIGHT 1971
C     MARK 3 REVISED
C     MARK 4A REVISED, IER-45
C     MARK 4.5 REVISED
C     MARK 7 REVISED (DEC 1978)
C     MARK 11 REVISED (FEB 1984)
C     RETURNS THE VALUE OF ERROR OR TERMINATES THE PROGRAM.
      INTEGER ERROR, IFAIL, NOUT
C$P 1
      DOUBLE PRECISION SRNAME
C     TEST IF NO ERROR DETECTED
      IF (ERROR.EQ.0) GO TO 20
C     DETERMINE OUTPUT UNIT FOR MESSAGE
      NOUT=6
C      CALL X04AAF (0,NOUT)
C     TEST FOR SOFT FAILURE
      IF (MOD(IFAIL,10).EQ.1) GO TO 10
C     HARD FAILURE
      WRITE (NOUT,99999) SRNAME, ERROR
C     ******************** IMPLEMENTATION NOTE ********************
C     THE FOLLOWING STOP STATEMENT MAY BE REPLACED BY A CALL TO AN
C     IMPLEMENTATION-DEPENDENT ROUTINE TO DISPLAY A MESSAGE AND/OR
C     ABORT THE PROGRAM.
C     *************************************************************
      STOP
C     SOFT FAIL
C     TEST IF ERROR MESSAGES SUPPRESSED
   10 IF (MOD(IFAIL/10,10).EQ.0) GO TO 20
      WRITE (NOUT,99999) SRNAME, ERROR
   20 P01AAF = ERROR
      RETURN
99999 FORMAT (1H0, 38HERROR DETECTED BY NAG LIBRARY ROUTINE , A8,
     * 11H - IFAIL = , I5//)
      END
** END OF P01AAFTEXT

ext: for

max_stars_repo_name: jlprieur/jlplib

max_issues_repo_name: jlprieur/jlplib

max_forks_repo_name: jlprieur/jlplib


text:
Bad Mans Plaything (Achillea sp) is typically known as yarrow, and is most often cultivated as an ornamental in backyard gardens for its flattopped flower clusters and scented foliage. Other colorful names for yarrow include Milfoil, Staunchweed, Soldiers Woundwort, Bloodwort, Nose Bleed, Knights Milfoil, Old Mans Pepper, Herbe Militaris, Thousand Weed, Carpenters Weed, Sanguinary, Devils Plaything, Devils Nettle and Yarroway. (God, I love these common names!!)

These perennial herbs within the aster family (Asteraceae) have been used for centuries in a medicinal manner, to stop bleeding, inhibiting bacterial growth, encourage healing, and to reduce inflammation. In Greek mythology Achilles (hence, Achillea) used the plant to heal the wounds of his soldiers and stop bleeding. With its strong aroma and flavor, yarrow is also used sparingly to flavor salads, soups, egg dishes, and may also be used to brew herb teas. Yarrow was also once used in brewing beer. The stalks of yarrow are dried and used in I Ching divination, and this summer when you open the doors to barbecue, use a fragrant smudge of yarrow leaves as a mosquito and insect repellent.

Several species of yarrow are grown solely for their foliage and flowers. Of the forty species within this genus, plants can be matforming or as tall as fourfeet, and flower colors range from whites and yellows to bold pinks.  All have finely divided, almost feathery leaves, and tiny flowers that form a flattopped cluster. Common yarrow (A. millefolium) has white flowers in its native northern temperate regions of Europe and North America. Yarrow is virtually a weed, growing in fields and urban waste places, from desert to high alpine, drought resistant and still  flowering throughout the weather summer  a perfect Davis plant! Yarrows are easy to grow, and need little care, except some grooming to remove spent flowers. In general, Achillea requires welldraining soil, a sunny location, and can become leggy if grown in rich soil or when fertilized. They can be planted spring or autumn from either seed or cuttings. They are propagated by division or cuttings in either fall or spring.

For a listing of other plants found growing in Davis, visit our Town Flora.

ext: f

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
module MoleculeMod
  implicit none
 private
 public :: InitMolecule, h2fraction,polyQ

 real, parameter :: xifuv=3             ! 2-4
 real, parameter :: mu=17.5              ! formation rate const. mu, min: 3.5
 integer, parameter :: poly=1           ! polytropic index
 real, parameter :: nzero=1500          ! cm^-3
 real, parameter :: pzero=1.e4          ! K cm^-3
 real, parameter :: pcut=1000.		! K cm^-3 (transition pressure)
 real, parameter :: k0=4.e-11		! dissociation rate s^-1
 real, parameter :: kboltz=1.38e-23     ! SI
 real, parameter :: amu=1.67e-27        ! SI
 real, parameter :: sqrtpi=1.772453851
 real, parameter :: year=3.15576e7
 real, parameter :: sigmav=4.9e-22	! cm^2
 real, parameter :: pc=3.086e18		! cm

 real, parameter :: Tform=1000,Tdes=3000,T4=10000	! temperature intervals

 real, save :: metal,meanmwt,meann,fmin,n0con

 logical, parameter :: debug=.FALSE.

 contains

subroutine InitMolecule(metallicity)
  use ElementsMod
 character(len=*), optional :: metallicity

 print*,' using molecule formation routines!'
 print*,'  > formation rate constant:', mu
 if(poly.EQ.0) print*,'  > constant density clouds'
 if(poly.EQ.1) print*,'  > logotropic clouds'

 if(present(metallicity)) then
  call InitElements(metallicity)
 else
  call InitElements()
 endif

 metal=zQ()/solarzQ()
 meanmwt=meanmwtQ()
 meann=1./fractionQ(elementnrQ("H   "))
 n0con=xifuv*sigmav*pc

 if(poly.NE.0.AND.poly.NE.1) call molError(1)

 if(poly.eq.0) fmin=0.
 if(poly.eq.1) fmin=1.e-6

end subroutine

! dt in year, n in cm^-3, T in K , g0 in units of Draine field, disp in km/s
subroutine h2fraction(fh2,dt,n,T,g0,disp,z,xe)
  real,intent(inout) :: fh2
  real,intent(in) :: n,T,g0,dt,disp
  real,intent(in), optional :: z,xe
  real ptot,lz,afuv,amean,ficonstant,dtau,snzero,ipar,Rf,acr,lxe
  integer :: i,j

 if(.not.present(z)) then
  lz=metal
 else
  lz=z
 endif

 if(.not.present(xe)) then
  lxe=0.0
 else
  lxe=xe
 endif

 if(debug) print*,'*1*',fh2,dt,n,T,g0,disp,lz

 fh2=max(fmin,fh2)

 if(dt.LE.0) return

 if(n.GT.1.e10.AND.T.LE.Tform) then
  fh2=1.
  return
 endif
 if(T.GT.T4) then
  fh2=fmin
  return
 endif
 if(T.GT.Tdes) then
  dtau=dt/taudiss(T,n,lxe)
  call fh2des(dtau,fh2)
 endif
 if(debug) print*,'*2*',fh2,dtau
! ficonstant=6.6e-6*sqrtpi*(lz*xifuv)**.5
 ficonstant=3.39e-6*sqrtpi*(lz*xifuv)**.5
 ptot=(meann+lxe)*n*T+1.e6*meanmwt*meann*n*disp**2*amu/kboltz/3.
 amean=meanAv(ptot,lz)
 afuv=xifuv*atrans(fh2,amean)
 if(debug) print*,'*3*', ficonstant,ptot,amean,afuv
 if(T.GT.Tform) then
  dtau=2*g0*k0*ficonstant*dt*year
!  afuv=log(exp(afuv)+dtau)
  afuv=afuv+log(1+dtau*exp(-afuv))
 endif
 if(debug) print*,'*4*', dtau,afuv
 if(T.LE.Tform) then
  Rf=rform(T,lz)
  dtau=2*n*Rf*dt*year
  ipar=g0*k0*ficonstant/n/Rf
  if(debug) print*,'*5*',Rf,dtau,ipar
  if(poly.EQ.0) then
   call solveh2_explicit(ipar,dtau,afuv)
!   j=INT(10*dtau)+1
!   do i=1,j
!    call solveh2(ipar,dtau/j,afuv)
!   enddo
   if(debug) print*,'*6*',j,afuv
  endif
  if(poly.EQ.1) then
   dtau=dtau*2./3.
   snzero=2./3.*nzero*n0con*(ptot/pzero)**.5
   call solveh2log_explicit(ipar,dtau,snzero,afuv)
!   j=INT(10*dtau)+1
!   do i=1,j
!    call solveh2log(ipar,dtau/j,snzero,afuv)
!   enddo
   if(debug) print*,'*6*',j,afuv,snzero
  endif
 endif
 acr=afuv/xifuv
 fh2=fmol(acr,amean)
 fh2=max(fmin,fh2)
end subroutine

subroutine molError(i)
  integer i

 if(i.eq.1) print*,' poly not equal to 0 or 1'

 stop

end subroutine

function rform(T,z) ! formation rate in cm^3 s^-1
 real :: rform,z,T
 real, parameter :: Tref=102
 rform=3.5e-18*z*T**.5/(1+T/Tref)**2*mu
 if(T.gt.Tform) rform=0
end function

function fmol(acr,av)
  real :: fmol,acr,av
 if(poly.EQ.0) then
  if(acr.gt.3./4.*av) then
   fmol=0
  else
   fmol=(1.-4./3.*acr/av)**3
  endif
 endif
 if(poly.EQ.1) then
  fmol=exp(-4.*acr/av)
 endif
end function

function atrans(fm,av)
  real :: fm, atrans,av
 if(poly.eq.0) then
  atrans=3./4.*av*(1-fm**(1./3.))
 endif
 if(poly.eq.1) then
  atrans=-av/4.*log(fm)
 endif
end function

function meanAv(pe,z)
  real :: pe,z,meanAv
 if(pe.gt.pcut) then
  meanAv=.22*z*(nzero/100)*sqrt(pe/pzero)
 else
  meanAv=.22*z*(nzero/100)*pe/pzero*sqrt(pzero/pcut)
 endif
end function

subroutine solveh2(alphaG,dbeta,sn)
 real, intent(in) :: alphaG, dbeta
 real, intent(inout) :: sn
 real, parameter :: precision=0.0001
 real :: fi,dfi,sn1, sn2,esn
 integer :: i
 integer,parameter :: imax=100

 i=0
 sn2=sn
10 continue
 i=i+1
 sn1=sn2
 esn=exp(-sn1)
 fi=sn1-sn-dbeta*(alphaG*esn-sn1)
 dfi=1+dbeta*(alphaG*esn+1)
 sn2=sn1-fi/dfi
 if(i.gt.imax) return
 if(ABS(sn2-sn1).GT.precision*sn1) goto 10
 sn=sn2
end subroutine

subroutine solveh2log(alphaG,dbeta,snzero,sn)
 real, intent(in) :: alphaG, dbeta,snzero
 real, intent(inout) :: sn
 real, parameter :: precision=0.0001
 real :: fi,dfi,sn1, sn2,esn,esn1
 integer :: i
 integer,parameter :: imax=100

 i=0
 if(sn.gt.50*snzero) sn=50*snzero
 sn2=sn
10 continue
 i=i+1
 if(sn2.gt.50*snzero) sn2=50*snzero
 sn1=sn2
 esn=exp(-sn1)
 esn1=exp(sn1/snzero)
 fi=sn1-sn-dbeta*(alphaG*esn-snzero*(esn1-1))
 dfi=1+dbeta*(alphaG*esn+esn1)
 sn2=sn1-fi/dfi
 if(i.gt.imax) return
 if(ABS(sn2-sn1).GT.precision*sn1) goto 10
 sn=sn2
end subroutine

subroutine solveh2log_explicit(alphaG,dbeta,snzero,sn)
 real, intent(in) :: alphaG, dbeta,snzero
 real, intent(inout) :: sn
 real, parameter :: fac=0.001
 real :: fi,dfi,sn1, sn2,esn,esn1,dtau,tau

 sn1=sn
 tau=0
 do while(tau.LT.dbeta)
   esn=exp(-sn1)
   esn1=exp(sn1/snzero)
   dfi=(alphaG*esn-snzero*(esn1-1))
   dtau=abs(fac*sn1/dfi)
   if(sn1.EQ.0) dtau=abs(1.e-8/dfi)
   if((dbeta-tau).LT.dtau) then
     dtau=dbeta-tau
     tau=dbeta
   else
     tau=tau+dtau
   endif
   if(dfi.EQ.0) then
     exit
   endif
   sn1=sn1+dtau*dfi
 enddo
 sn=sn1
end subroutine

subroutine solveh2_explicit(alphaG,dbeta,sn)
 real, intent(in) :: alphaG, dbeta
 real, intent(inout) :: sn
 real, parameter :: fac=0.001
 real :: fi,dfi,sn1, sn2,esn,esn1,dtau,tau

 sn1=sn
 tau=0
 do while(tau.LT.dbeta)
   esn=exp(-sn1)
   dfi=(alphaG*esn-sn1)
   dtau=abs(fac*sn1/dfi)
   if(sn1.EQ.0) dtau=abs(1.e-8/dfi)
   if((dbeta-tau).LT.dtau) then
     dtau=dbeta-tau
     tau=dbeta
   else
     tau=tau+dtau
   endif
   if(dfi.EQ.0) then
     exit
   endif
   sn1=sn1+dtau*dfi
 enddo
 sn=sn1
end subroutine

subroutine fh2des(beta,fh2) ! collissional h2 destruction (gamma2=gamma1/20)
 real :: beta,fh2
 fh2=20*fh2*exp(-beta)/(19*fh2*exp(-beta)+(20-19*fh2))
end subroutine

 function taudiss(T,n,xe) result(tau)
  real :: tau
  real, intent(in) :: T,n,xe
  real :: tau1,tau2,lxe

  lxe=MIN(1.,xe)

  tau1=taudiss1(T,n)
  tau2=taudiss2(T,n)
  tau=1./((1-lxe)/tau1+lxe/tau2)
!  print*,tau,tau1,tau2
 end function

 function taudiss1(T,n) result(tau)
  real :: tau
  real, intent(in) :: T,n
  real, parameter :: Tref=3162,tauref=3.2e10,Tmin=3000,Tmax=10000
!  if(T.LT.Tmin) then
!   tau=tauref*(Tref/Tmin)**12/n
!   return
!  endif
  if(T.GT.Tmax) then
   tau=tauref*(Tref/Tmax)**12/n
   return
  endif
  tau=tauref*(Tref/T)**12/n
 end function

 function taudiss2(T,n) result(tau)
  real :: tau
  real, intent(in) :: T,n
  real, parameter :: Tref=3162,tauref=3.2e7,Tmin=3000,Tmax=10000
!  if(T.LT.Tmin) then
!   tau=tauref*(Tref/Tmin)**12/n
!   return
!  endif
  if(T.GT.Tmax) then
   tau=tauref*(Tref/Tmax)**6/n
   return
  endif
  tau=tauref*(Tref/T)**6/n
 end function

 function polyQ()
  integer polyQ
  polyQ=poly
 end function

end module

subroutine h2test
 use MoleculeMod
 real :: fh2,dt,n,T,g0,disp,p
 real :: tform,feq,tdiss,tdes,xe
 integer :: i,j,k,l

 print*,'molecule formation routine test'

 call InitMolecule

print*,' fh2,dt,n,T,disp,g0,xe'
10 read*,fh2,dt,n,T,disp,g0,xe
 if(n.eq.0) stop
  call h2fraction(fh2,dt,n,T,g0,disp,xe=xe)
 print*,dt,fh2
 goto 10
end

max_stars_repo_name: rknop/amuse

max_issues_repo_name: rknop/amuse

max_forks_repo_name: rieder/amuse


text:
!
!     (c) 2019 Guide Star Engineering, LLC
!     This Software was developed for the US Nuclear Regulatory Commission (US NRC)
!     under contract "Multi-Dimensional Physics Implementation into Fuel Analysis under
!     Steady-state and Transients (FAST)", contract # NRC-HQ-60-17-C-0007
!
program test_keytable
  !! label: Morfeus-FV
  !!
  !! Test keytable class behavior

    use assertions_interface, only : assert
    use class_keytable, only : keytable

    implicit none
    integer, allocatable :: row1(:), row2(:), row3(:)
    integer, pointer :: row1T(:), row2T(:), row3T(:)
    type(keytable) :: testTable

    row1 = [1,2,3]
    row2 = [4,5,6,7]
    row3 = [8,9,10,11,12,13]

    call testTable%alloc_keytable(1,3)
    call testTable%set_kt_row(1,row1)
    call testTable%set_kt_row(2,row2)
    call testTable%set_kt_row(3,row3)

    call testTable%get_kt_row(1,row1T)
    call testTable%get_kt_row(2,row2T)
    call testTable%get_kt_row(3,row3T)

    call assert( row1T == row1, "test_keytable: testTable%row(1) == row1" )
    call assert( row2T == row2, "test_keytable: testTable%row(2) == row2" )
    call assert( row3T == row3, "test_keytable: testTable%row(3) == row3" )
    call assert( testTable%get_row_size(1) == size(row1), "test_keytable: testTable%get_row_size(1) == size(row1)" )
    call assert( testTable%get_row_size(2) == size(row2), "test_keytable: testTable%get_row_size(2) == size(row2)" )
    call assert( testTable%get_row_size(3) == size(row3), "test_keytable: testTable%get_row_size(3) == size(row3)" )

    print *, "Test passed."

end program test_keytable

max_stars_repo_name: sourceryinstitute/MORFEUS-Source

max_issues_repo_name: sourceryinstitute/MORFEUS-Source

max_forks_repo_name: sourceryinstitute/MORFEUS-Source


PROGRAM disturb_atomic_coordinates
!MODULE disturb_module
!
  implicit none
!
!  PRIVATE
!  PUBLIC :: disturb

  real(8) :: ratio
  real(8) :: ax,aa(3,3)
  real(8),allocatable :: aa_atom(:,:)
  integer,allocatable :: ki_atom(:),md_atom(:)
  integer,parameter :: u=970, w=971
  integer :: nelem,natom,zatom(9),i
  character(3) :: cbuf,cformat

  write(*,*) "input the ratio of disturbance [0.0~1.0]"
  read(*,*) ratio
  write(*,*) "atomic coordinates are read from fort.970"
  read(u,*) cbuf,ax
  read(u,*) cbuf,aa(1:3,1)
  read(u,*) cbuf,aa(1:3,2)
  read(u,*) cbuf,aa(1:3,3)
  read(u,*) cformat
  read(u,*) nelem, natom, zatom(1:nelem)
  allocate( aa_atom(3,natom) ) ; aa_atom=0.0d0
  allocate( ki_atom(natom)   ) ; ki_atom=0
  allocate( md_atom(natom)   ) ; md_atom=0
  do i=1,natom
     read(u,*) ki_atom(i),aa_atom(1:3,i),md_atom(i)
  end do

  call disturb( aa_atom, ratio )

  write(*,*) "the disturbed atomic coordinates are written in fort.971"

  write(w,'("AX", f20.16)') ax
  write(w,'("A1",3f20.16)') aa(1:3,1)
  write(w,'("A2",3f20.16)') aa(1:3,2)
  write(w,'("A3",3f20.16)') aa(1:3,3)
  write(w,'(a3)') cformat
  write(w,'(1x,i4,10i8)') nelem,natom,zatom(1:nelem)
  do i=1,natom
     write(w,'(1x,i5,3f20.12,i4)') ki_atom(i),aa_atom(1:3,i),md_atom(i)
  end do

CONTAINS

  SUBROUTINE disturb( a, ratio )

    implicit none
    real(8),intent(INOUT) :: a(:,:)
    real(8),intent(IN) :: ratio
    real(8),allocatable :: d(:)
    real(8) :: r
    integer :: i

    if ( ratio <= 0.0d0 ) return
    r = min( ratio, 1.0d0 )

    allocate( d(size(a,1)) ) ; d=0.0d0

    do i=1,size(a,2)
       call random_number( d )
       d(:) = ( 2.0d0*d(:) - 1.0d0 )*r
       a(:,i) = a(:,i) + d(:)
!       write(*,'(1x,i4,2(5x,3f15.10))') i,a(:,i)-d(:),a(:,i)
    end do

    deallocate( d )

  END SUBROUTINE disturb

END PROGRAM disturb_atomic_coordinates
!END MODULE disturb_module


text:
      subroutine etacf(pr,prc,x)
      COMMON/cfi/icf
      pr1=1.2
      x1=1.02
      xc=0.950
      prx=pr
      prcx=prc*xc
      if (pr.lt.pr1) prx=pr1
      if (pr.gt.prcx) prx=prcx
      x=1.+(prx-prcx)/(pr1-prcx)*(x1-1.)
      if (icf.eq.1) x=1.
      return
      end

max_stars_repo_name: mjfwest/OpenMDAO-Framework

max_issues_repo_name: jcchin/OpenMDAO-Framework

max_forks_repo_name: jcchin/OpenMDAO-Framework


text:
This editor can edit this entry and tell us a bit about themselves by clicking the Edit icon.

20100520 18:05:10 nbsp opoopoopo!!!! Users/JabberWokky


text:
      SUBROUTINE totzspa(rzl_array, r11, ru1, rv1, z11, zu1, zv1, lu1,
     1   lv1, rcn1, zcn1)
      USE vmec_main
      USE vmec_params, ONLY: jmin1, jlam, ntmax, rcs, rsc, zcc, zss
      IMPLICIT NONE
C-----------------------------------------------
C   D u m m y   A r g u m e n t s
C-----------------------------------------------
      REAL(rprec), DIMENSION(ns,0:ntor,0:mpol1,3*ntmax),
     1   TARGET, INTENT(inout) :: rzl_array
      REAL(rprec), DIMENSION(ns*nzeta,ntheta3,0:1), INTENT(out) ::
     1   r11, ru1, rv1, z11, zu1, zv1, lu1, lv1, rcn1, zcn1
C-----------------------------------------------
C   L o c a l   P a r a m e t e r s
C-----------------------------------------------
      INTEGER, PARAMETER :: m0 = 0
C-----------------------------------------------
C   L o c a l   V a r i a b l e s
C-----------------------------------------------
      INTEGER :: m, n, mparity, k, i, l, j1, j1l, nsl
      INTEGER :: ioff, joff, mj, ni
      REAL(rprec), DIMENSION(:,:,:), POINTER ::
     1           rmncs, rmnsc, zmncc, zmnss, lmncc, lmnss
      REAL(rprec), DIMENSION(:,:), ALLOCATABLE :: work1
      REAL(rprec) :: cosmux, sinmux
C-----------------------------------------------
      rmnsc => rzl_array(:,:,:,rsc)               !!SIN(mu) COS(nv)
      zmncc => rzl_array(:,:,:,zcc+ntmax)         !!COS(mu) COS(nv)
      lmncc => rzl_array(:,:,:,zcc+2*ntmax)       !!COS(mu) COS(nv)
      IF (lthreed) THEN
         rmncs => rzl_array(:,:,:,rcs)               !!COS(mu) SIN(nv)
         zmnss => rzl_array(:,:,:,zss+ntmax)         !!SIN(mu) SIN(nv)
         lmnss => rzl_array(:,:,:,zss+2*ntmax)       !!SIN(mu) SIN(nv)
      END IF

      ioff = LBOUND(rmncs,2)
      joff = LBOUND(rmncs,3)

      z00b = zmncc(ns,ioff,joff)

      ALLOCATE (work1(ns*nzeta,12), stat=i)
      IF (i .ne. 0) STOP 'Allocation error in VMEC totzspa'

!
!                INITIALIZATION BLOCK
!
      r11 = 0;  ru1 = 0;  rv1 = 0;  z11 = 0;  zu1 = 0
      zv1 = 0;  lu1 = 0;  lv1 = 0;  rcn1 = 0; zcn1 = 0

      IF (jlam(m0) .gt. 1) lmncc(1,:,m0+joff) = lmncc(2,:,m0+joff)

      DO m = 0, mpol1
         mparity = MOD(m,2)
         mj = m+joff
         work1 = 0
         j1 = jmin1(m)
         DO n = 0, ntor
            ni = n+ioff
            DO k = 1, nzeta
               l = ns*(k-1)
               j1l = j1+l;  nsl = ns+l
               work1(j1l:nsl,1) = work1(j1l:nsl,1)
     1                          + rmnsc(j1:ns,ni,mj)*cosnv(k,n)
               work1(j1l:nsl,6) = work1(j1l:nsl,6)
     1                          + zmncc(j1:ns,ni,mj)*cosnv(k,n)
               work1(j1l:nsl,10) = work1(j1l:nsl,10)
     1                          + lmncc(j1:ns,ni,mj)*cosnv(k,n)

               IF (.not.lthreed) CYCLE

               work1(j1l:nsl,2) = work1(j1l:nsl,2)
     1                          + rmncs(j1:ns,ni,mj)*sinnv(k,n)
               work1(j1l:nsl,3) = work1(j1l:nsl,3)
     1                          + rmnsc(j1:ns,ni,mj)*sinnvn(k,n)
               work1(j1l:nsl,4) = work1(j1l:nsl,4)
     1                          + rmncs(j1:ns,ni,mj)*cosnvn(k,n)
               work1(j1l:nsl,5) = work1(j1l:nsl,5)
     1                          + zmnss(j1:ns,ni,mj)*sinnv(k,n)
               work1(j1l:nsl,7) = work1(j1l:nsl,7)
     1                          + zmnss(j1:ns,ni,mj)*cosnvn(k,n)
               work1(j1l:nsl,8) = work1(j1l:nsl,8)
     1                          + zmncc(j1:ns,ni,mj)*sinnvn(k,n)
               work1(j1l:nsl,9) = work1(j1l:nsl,9)
     1                          + lmnss(j1:ns,ni,mj)*sinnv(k,n)
               work1(j1l:nsl,11) = work1(j1l:nsl,11)
     1                          + lmnss(j1:ns,ni,mj)*cosnvn(k,n)
               work1(j1l:nsl,12) = work1(j1l:nsl,12)
     1                          + lmncc(j1:ns,ni,mj)*sinnvn(k,n)
            END DO
         END DO

!
!        INVERSE TRANSFORM IN M-THETA
!
         DO i = 1, ntheta2
            cosmux = xmpq(m,1)*cosmu(i,m)
            sinmux = xmpq(m,1)*sinmu(i,m)
            r11(:,i,mparity) = r11(:,i,mparity) + work1(:,1)*
     1            sinmu(i,m)
            ru1(:,i,mparity) = ru1(:,i,mparity) + work1(:,1)*
     1            cosmum(i,m)
            z11(:,i,mparity) = z11(:,i,mparity) + work1(:,6)*
     1            cosmu(i,m)
            zu1(:,i,mparity) = zu1(:,i,mparity) + work1(:,6)*
     1            sinmum(i,m)
            lu1(:,i,mparity) = lu1(:,i,mparity) + work1(:,10)*
     1            sinmum(i,m)
            rcn1(:,i,mparity) = rcn1(:,i,mparity) + work1(:,1)*
     1            sinmux
            zcn1(:,i,mparity) = zcn1(:,i,mparity) + work1(:,6)*
     1            cosmux

            IF (.not.lthreed) CYCLE

            r11(:,i,mparity) = r11(:,i,mparity) + work1(:,2)*
     1               cosmu(i,m)
            ru1(:,i,mparity) = ru1(:,i,mparity) + work1(:,2)*
     1               sinmum(i,m)
            z11(:,i,mparity) = z11(:,i,mparity) + work1(:,5)*
     1               sinmu(i,m)
            zu1(:,i,mparity) = zu1(:,i,mparity) + work1(:,5)*
     1               cosmum(i,m)
            lu1(:,i,mparity) = lu1(:,i,mparity) + work1(:,9)*
     1               cosmum(i,m)
            rcn1(:,i,mparity) = rcn1(:,i,mparity) + work1(:,2)*
     1               cosmux
            zcn1(:,i,mparity) = zcn1(:,i,mparity) + work1(:,5)*
     1               sinmux
            rv1(:,i,mparity) = rv1(:,i,mparity) + work1(:,3)*
     1               sinmu(i,m) + work1(:,4)*cosmu(i,m)
            zv1(:,i,mparity) = zv1(:,i,mparity) + work1(:,7)*
     1               sinmu(i,m) + work1(:,8)*cosmu(i,m)
            lv1(:,i,mparity) = lv1(:,i,mparity) - (work1(:,11)*
     1               sinmu(i,m)+work1(:,12)*cosmu(i,m))
         END DO
      END DO

      DEALLOCATE (work1)

      END SUBROUTINE totzspa


text:
      subroutine intbsw(nbasp,nbaspx,nbt,nbtmax,nobswt,noutpt,nst,
     $ nstmax,nttyo,uobsw,uspec)
c
c     This subroutine interprets ordinary basis switching directives
c     specified on the input file.
c
c     This subroutine is called by:
c
c       EQ3NR/eq3nr.f
c       EQ6/eq6.f
c
c-----------------------------------------------------------------------
c
c     Principal input:
c
c       nbasp  = the indices of the species in the basis set
c                  (modified by this subroutine)
c       nbt    = the number of species in the basis set
c       nobswt = the number of basis switches to make
c       uobsw  = the name pairs for the species to be switched
c       uspec  = array of species names
c
c     Principal output:
c
c       nbasp  = the indices of the species in the basis set
c                  (modified by this subroutine)
c
c-----------------------------------------------------------------------
c
      implicit none
c
c-----------------------------------------------------------------------
c
c     Calling sequence variable declarations.
c
      integer nbtmax,nstmax
c
      integer noutpt,nttyo
c
      integer nbasp(nbtmax),nbaspx(nbtmax)
      integer nobswt,nbt,nst
c
      character(len=48) uobsw(2,nbtmax),uspec(nstmax)
c
c-----------------------------------------------------------------------
c
c     Local variable declarations.
c
      integer jlen,n,nb,nerr,ns1,ns2
c
      character(len=56) uspn56
      character(len=48) unam48
c
c-----------------------------------------------------------------------
c
      nerr = 0
c
      do n = 1,nobswt
c
c       Find the first species.
c
        do nb = 1,nbt
          ns1 = nbaspx(nb)
          if (uspec(ns1)(1:48) .eq. uobsw(1,n)(1:48)) go to 110
        enddo
c
        nerr = nerr + 1
c
        unam48 = uobsw(1,n)
        call fmspnx(jlen,unam48,uspn56)
        write (noutpt,1000) uspn56(1:jlen)
        write (nttyo,1000) uspn56(1:jlen)
 1000   format(/' * Error - (EQLIB/intbsw) The species ',a,
     $  /7x,'is specified to be replaced in an ordinary basis switch,',
     $  /7x,'but it is not in the active basis set.')
        go to 150
c
  110   continue
c
c       Find the second species.
c
        do ns2 = 1,nst
          if (uspec(ns2)(1:48) .eq. uobsw(2,n)(1:48)) go to 130
        enddo
c
        nerr = nerr + 1
c
        unam48 = uobsw(2,n)
        call fmspnx(jlen,unam48,uspn56)
        write (noutpt,1010) uspn56(1:jlen)
        write (nttyo,1010) uspn56(1:jlen)
 1010   format(/' * Error - (EQLIB/intbsw) The species ',a,
     $  /7x,'is specified to be put into the basis set by an ordinary',
     $  /7x,'basis switch but it is not in the set of active',
     $  ' species.')
        go to 150
c
  130   continue
c
c       Mark to switch.
c
        nbasp(nb) = ns2
c
  150   continue
      enddo
c
c* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
c
c     Stop if name match errors have been encountered.
c
      if (nerr .gt. 0) then
        write (noutpt,1020) nerr
        write (nttyo,1020) nerr
 1020   format(/' * Error - (EQLIB/intbsw) ',i4,' errors were',
     $  /7x,'encountered interpreting ordinary basis switches that',
     $  /7x,'were directed on the input file.')
        stop
      endif
c
      end


text:
      SUBROUTINE CALFQC(ISTL_,IS2TL_,MVAR,MO,CON,CON1,FQCPAD,QSUMPAD,
     &    QSUMNAD)
C
C CHANGE RECORD
C **  SUBROUTINE CALFQC CALCULATES MASS SOURCES AND SINKS ASSOCIATED
C **  WITH CONSTANT AND TIME SERIES INFLOWS AND OUTFLOWS; CONTROL
C **  STRUCTURE INFLOWS AND OUTLOWS; WITHDRAWAL AND RETURN STRUCTURE
C **  OUTFLOWS; AND  EMBEDED CHANNEL INFLOWS AND OUTFLOWS
C
      USE GLOBAL

	INTEGER::L,K,ID,JD,KD,NWR,IU,JU,KU,LU,NS
	INTEGER::LD,NMD,NJP

      DIMENSION CON(LCM,KCM),CON1(LCM,KCM),FQCPAD(LCM,KCM),
     &          QSUMNAD(LCM,KCM),QSUMPAD(LCM,KCM)

      REAL,SAVE,ALLOCATABLE,DIMENSION(:,:)::CONQ
      IF(.NOT.ALLOCATED(CONQ))THEN
        ALLOCATE(CONQ(LCM,KCM))
        CONQ=0.0
      ENDIF
C

      M=MO

      ! *** SELECTIVE ZEROING
      IF(KC.GT.1)THEN
        IF(NGWSER.GT.0.OR.ISGWIT.NE.0)THEN
          DO L=1,LC
            FQC(L,1)=0.
          ENDDO
        ENDIF

        ! *** ZERO EVAP/RAINFALL
        IF(MVAR.EQ.2)THEN
          DO L=1,LC
            FQC(L,KC)=0.
          ENDDO
          IF(ISADAC(MVAR).GE.2)THEN
            DO L=1,LC
              FQCPAD(L,KC)=0.
            ENDDO
          ENDIF
          IF(ISADAC(MVAR).GT.0)THEN
            DO L=1,LC
              QSUMPAD(L,KC)=0.
            ENDDO
          ENDIF
        ENDIF

        ! *** ZERO ALL DEFINED BC'S
        DO NS=1,NBCS
          L=LBCS(NS)
          DO K=1,KC
            FQC(L,K)=0.
            FQCPAD(L,K)=0
            QSUMPAD(L,K)=0.
            !QSUMNAD(L,K)=0.
          ENDDO
        ENDDO

      ELSE
        FQC=0.
        IF(ISADAC(MVAR).GE.2)FQCPAD=0.
        QSUMPAD=0.
        !QSUMNAD=0.
      ENDIF
C
      IF(MVAR.EQ.8.AND.IWQPSL.NE.2) GOTO 1500
C
C **  INITIALIZE VOLUMETRIC SOURCE-SINK FLUXES AND AUXILLARY VARIABLES
C
      ! *** 3TL STANDARD & WATER QUALITY
      IF(ISTL_.EQ.3.OR.MVAR.EQ.8)THEN
        IF(NGWSER.GT.0.OR.ISGWIT.NE.0)THEN
          DO L=1,LC
            CONQ(L,1)=CON(L,1)
          ENDDO
        ENDIF

        ! *** ZERO EVAP/RAINFALL
        IF(MVAR.EQ.2)THEN
          DO L=1,LC
            CONQ(L,KC)=CON(L,KC)
          ENDDO
        ENDIF

        ! *** INITIALIZE ALL DEFINED BC'S
        DO NS=1,NBCS
          L=LBCS(NS)
          DO K=1,KC
            CONQ(L,K)=CON(L,K)
          ENDDO
        ENDDO
      ENDIF
C
      ! *** 3TL CORRECTION STEP
      IF(ISTL_.EQ.2.AND.IS2TL_.EQ.0)THEN
        IF(NGWSER.GT.0.OR.ISGWIT.NE.0)THEN
          DO L=1,LC
            CONQ(L,1)=0.5*(CON(L,1)+CON1(L,1))
          ENDDO
        ENDIF

        ! *** ZERO EVAP/RAINFALL
        IF(MVAR.EQ.2)THEN
          DO L=1,LC
            CONQ(L,KC)=0.5*(CON(L,KC)+CON1(L,KC))
          ENDDO
        ENDIF

        ! *** INITIALIZE ALL DEFINED BC'S
        DO NS=1,NBCS
          L=LBCS(NS)
          DO K=1,KC
            CONQ(L,K)=0.5*(CON(L,K)+CON1(L,K))
          ENDDO
        ENDDO

      ENDIF

      ! *** 2TL STANDARD
      IF(ISTL_.EQ.2.AND.IS2TL_.EQ.1)THEN
        IF(NGWSER.GT.0.OR.ISGWIT.NE.0)THEN
          DO L=1,LC
            CONQ(L,1)=0.5*(3.*CON(L,1)-CON1(L,1))
          ENDDO
        ENDIF

        ! *** ZERO EVAP/RAINFALL
        IF(MVAR.EQ.2)THEN
          DO L=1,LC
            CONQ(L,KC)=0.5*(3.*CON(L,KC)-CON1(L,KC))
          ENDDO
        ENDIF

        ! *** INITIALIZE ALL DEFINED BC'S
        DO NS=1,NBCS
          L=LBCS(NS)
          DO K=1,KC
            CONQ(L,K)=0.5*(3.*CON(L,K)-CON1(L,K))
          ENDDO
        ENDDO

      ENDIF
C
      IF(MVAR.EQ.4) GOTO 1000
C
      IF(MVAR.EQ.8)THEN
        M=4+NTOX+NSED+NSND+MO
      ENDIF
C
C *********************************************************************C
C
C *** STANDARD VOLUMETRICS SOURCE SINK LOCATIONS (2TL)
C
      IF(ISTL_.EQ.2.AND.IS2TL_.EQ.1)THEN
        ! *** FLOW BOUNDARY CELLS (2TL)

        DO NS=1,NQSIJ
          L=LQS(NS)
          NQSTMP=NQSERQ(NS)
          NCSTMP=NCSERQ(NS,M)
          DO K=1,KC
            FQC(L,K)=FQC(L,K)
     &          +MAX(QSS(K,NS),0.)*CQS(K,NS,M)
     &          +MIN(QSS(K,NS),0.)*CONQ(L,K)
     &          +MAX(QSERCELL(K,NS),0.)*CSERT(K,NCSTMP,M)
     &          +MIN(QSERCELL(K,NS),0.)*CONQ(L,K)
            FQCPAD(L,K)=FQCPAD(L,K)
     &          +MAX(QSS(K,NS),0.)*CQS(K,NS,M)
     &          +MAX(QSERCELL(K,NS),0.)*CSERT(K,NCSTMP,M)
            QSUMPAD(L,K)=QSUMPAD(L,K)
     &          +MAX(QSS(K,NS),0.)+MAX(QSERCELL(K,NS),0.)
C            QSUMNAD(L,K)=QSUMNAD(L,K)
C     &          +MIN(QSS(K,NS),0.)+MIN(QSERCELL(K,NS),0.)
          ENDDO
!        if (mvar ==2) write(*,*) 'temp =',sum(FQC(L,:)),sum(CSERT(:,NCSTMP,2)),sum(CQS(:,NS,2))
!         if (mvar ==3) write(*,*) 'Dye=',sum(FQC(L,:)),sum(CSERT(:,NCSTMP,3)),sum(QSS(:,NS)),SUM(CONQ(L,:))

        ENDDO

        ! ***  JET-PLUME VOLUMETRICS SOURCE SINK LOCATIONS (2TL)
        IF(NQJPIJ.GT.0)THEN
          DO NJP=1,NQJPIJ
            IF(ICALJP(NJP).EQ.1)THEN
              RPORTS=FLOAT(NPORTJP(NJP))
              LJP=LIJ(IQJP(NJP),JQJP(NJP))
              KTMP=KEFFJP(NJP)
              ! ***  QVJPTMP=TIME SERIES DISCHARGE FROM JET-PLUME
              QVJPTMP=0.
              DO K=1,KC
                QVJPTMP=QVJPTMP+QSERT(K,NQSERJP(NJP))
              ENDDO

              ! QCJPTMP=ENTRAINMENT FLUX
              QCJPTMP=0.
              QVJPENT=0.
              ! REMOVE ENTRAINMENT FLUX AND CALCULATE TOTAL ENTRAIMENT FLUX
              DO K=1,KC
                FQC(LJP,K)=FQC(LJP,K)-RPORTS*QJPENT(K,NJP)*CONQ(LJP,K)
                QCJPTMP=QCJPTMP+QJPENT(K,NJP)*CONQ(LJP,K)
                QVJPENT=QVJPENT+QJPENT(K,NJP)
C                QSUMNAD(LJP,K)=QSUMNAD(LJP,K)-RPORTS*QJPENT(K,NJP)
              ENDDO

              ! PLACE JET FLUX AND ENTRAINMENT FLUX IS EFFECTIVE LAYER
              FQC(LJP,KTMP)=FQC(LJP,KTMP)+RPORTS*QCJPTMP
     &            +RPORTS*QQCJP(NJP)*CQCJP(1,NJP,M)
     &            +RPORTS*QVJPTMP*CSERT(1,NCSERJP(NJP,M),M)
              FQCPAD(LJP,KTMP)=FQCPAD(LJP,KTMP)+RPORTS*QCJPTMP
     &            +RPORTS*QQCJP(NJP)*CQCJP(1,NJP,M)
     &            +RPORTS*QVJPTMP*CSERT(1,NCSERJP(NJP,M),M)
              QSUMPAD(LJP,KTMP)=QSUMPAD(LJP,KTMP)+RPORTS*QVJPENT
     &            +RPORTS*QQCJP(NJP)+RPORTS*QVJPTMP
            ENDIF
            IF(ICALJP(NJP).EQ.2)THEN
              RPORTS=FLOAT(NPORTJP(NJP))
              LJP=LIJ(IQJP(NJP),JQJP(NJP))
              KTMP=KEFFJP(NJP)
              NS=NQWRSERJP(NJP)
              LU=LIJ(IUPCJP(NJP),JUPCJP(NJP))
              KU=KUPCJP(NJP)
              CONUP=CONQ(LU,KU)
              QCJPTMP=0.
              QVJPENT=0.

              ! REMOVE ENTRAIMENT FLUX AND CALCULATE TOTAL ENTRAINMENT
              DO K=1,KC
                FQC(LJP,K)=FQC(LJP,K)-RPORTS*QJPENT(K,NJP)*CONQ(LJP,K)
                QCJPTMP=QCJPTMP+QJPENT(K,NJP)*CONQ(LJP,K)
                QVJPENT=QVJPENT+QJPENT(K,NJP)
C                QSUMNAD(LJP,K)=QSUMNAD(LJP,K)-RPORTS*QJPENT(K,NJP)
              ENDDO

              !  PLACE ENTRAINMENT, CONSTANT AND TIME SERIES FLUXES IN EFFECTIVE CELL
              FQC(LJP,KTMP)=FQC(LJP,KTMP)+RPORTS*QCJPTMP
     &            +RPORTS*QWRCJP(NJP)*(CWRCJP(NJP,M)+CONUP)
     &            +RPORTS*QWRSERT(NS)*(CQWRSERT(NS,M)+CONUP)
              FQCPAD(LJP,KTMP)=FQCPAD(LJP,KTMP)+RPORTS*QCJPTMP
     &            +RPORTS*QWRCJP(NJP)*(CWRCJP(NJP,M)+CONUP)
     &            +RPORTS*QWRSERT(NS)*(CQWRSERT(NS,M)+CONUP)
              QSUMPAD(LJP,KTMP)=QSUMPAD(LJP,KTMP)+RPORTS*QVJPENT
     &            +RPORTS*QWRCJP(NJP)+RPORTS*QWRSERT(NS)

              ! REMOVAL WITHDRAWAL FROM UPSTREAM CELL
              FQC(LU,KU)=FQC(LU,KU)
     &            -RPORTS*QWRCJP(NJP)*CONUP
     &            -RPORTS*QWRSERT(NS)*CONUP
C              QSUMNAD(LU,KU)=QSUMNAD(LU,KU)
C     &            -RPORTS*QWRCJP(NJP)-RPORTS*QWRSERT(NS)
            ENDIF
          ENDDO
        ENDIF

        ! *** CONTROL STRUCTURES (2TL)
        DO NCTL=1,NQCTL
          RQWD=1.
          IU=IQCTLU(NCTL)
          JU=JQCTLU(NCTL)
          LU=LIJ(IU,JU)
          ID=IQCTLD(NCTL)
          JD=JQCTLD(NCTL)
          IF(ID.EQ.0.AND.JD.EQ.0)THEN
            LD=LC
            RQWD=0.
          ELSE
            LD=LIJ(ID,JD)
          ENDIF
          DO K=1,KC
            FQC(LU,K)=FQC(LU,K)
     &          -QCTLT(K,NCTL)*CONQ(LU,K)
            FQC(LD,K)=FQC(LD,K)
     &          +RQWD*QCTLT(K,NCTL)*CONQ(LU,K)
            FQCPAD(LD,K)=FQCPAD(LD,K)
     &          +RQWD*QCTLT(K,NCTL)*CONQ(LU,K)
            QSUMPAD(L,K)=QSUMPAD(L,K)
     &          +RQWD*QCTLT(K,NCTL)
C            QSUMNAD(L,K)=QSUMNAD(L,K)
C     &          -QCTLT(K,NCTL)
          ENDDO
        ENDDO

        ! *** WITHDRAWAL CONCENTRATION RISE RETURN (2TL)
        DO NWR=1,NQWR
          IU=IQWRU(NWR)
          JU=JQWRU(NWR)
          KU=KQWRU(NWR)
          ID=IQWRD(NWR)
          JD=JQWRD(NWR)
          KD=KQWRD(NWR)
          LU=LIJ(IU,JU)
          LD=LIJ(ID,JD)
          NQSTMP=NQWRSERQ(NWR)
          NCSTMP=NQWRSERQ(NWR)
          FQC(LU,KU)=FQC(LU,KU)
     &        -(QWR(NWR)+QWRSERT(NQSTMP))*CONQ(LU,KU)
          FQC(LD,KD)=FQC(LD,KD)
     &        +QWR(NWR)*(CONQ(LU,KU)+CQWR(NWR,M))
     &        +QWRSERT(NQSTMP)*(CONQ(LU,KU)+CQWRSERT(NCSTMP,M))
          FQCPAD(LD,KD)=FQCPAD(LD,KD)
     &        +QWR(NWR)*(CONQ(LU,KU)+CQWR(NWR,M))
     &        +QWRSERT(NQSTMP)*(CONQ(LU,KU)+CQWRSERT(NCSTMP,M))
          QSUMPAD(LD,KD)=QSUMPAD(LD,KD)
     &        +QWR(NWR)+QWRSERT(NQSTMP)
C          QSUMNAD(LU,KU)=QSUMNAD(LU,KU)
C     &        -(QWR(NWR)+QWRSERT(NQSTMP))
        ENDDO

        ! *** SUBGRID SCALE CHANNEL EXCHANGE (2TL)
        IF(MDCHH.GE.1)THEN
          DO K=1,KC
            DO NMD=1,MDCHH
              LMDCHHT=LMDCHH(NMD)
              LMDCHUT=LMDCHU(NMD)
              LMDCHVT=LMDCHV(NMD)
              IF(MDCHTYP(NMD).EQ.1)THEN
                QUKTMP=QCHANU(NMD)*DZC(K)
                QVKTMP=0.
              ENDIF
              IF(MDCHTYP(NMD).EQ.2)THEN
                QVKTMP=QCHANV(NMD)*DZC(K)
                QUKTMP=0.
              ENDIF
              IF(MDCHTYP(NMD).EQ.3)THEN
                QUKTMP=QCHANU(NMD)*DZC(K)
                QVKTMP=QCHANV(NMD)*DZC(K)
              ENDIF
              FQC(LMDCHHT,K)=FQC(LMDCHHT,K)
     &            +MAX(QUKTMP,0.)*CONQ(LMDCHUT,K)
     &            +MIN(QUKTMP,0.)*CONQ(LMDCHHT,K)
     &            +MAX(QVKTMP,0.)*CONQ(LMDCHVT,K)
     &            +MIN(QVKTMP,0.)*CONQ(LMDCHHT,K)
              FQC(LMDCHUT,K)=FQC(LMDCHUT,K)
     &            -MAX(QUKTMP,0.)*CONQ(LMDCHUT,K)
     &            -MIN(QUKTMP,0.)*CONQ(LMDCHHT,K)
              FQC(LMDCHVT,K)=FQC(LMDCHVT,K)
     &            -MAX(QVKTMP,0.)*CONQ(LMDCHVT,K)
     &            -MIN(QVKTMP,0.)*CONQ(LMDCHHT,K)
            ENDDO
          ENDDO
        ENDIF

        ! *** GROUNDWATER, EVAP, RAINFALL (2TL)
        IF(ISGWIE.NE.0)THEN
          DO L=2,LA
            FQC(L,1)=FQC(L,1)-RIFTR(L)*CONQ(L,1)
          ENDDO
        ENDIF

        ! *** ZONED SEEPAGE (2TL)
        IF(ISGWIT.EQ.3)THEN
          DO L=2,LA
            IF(H1P(L).GT.HDRY)THEN
              FQC(L,1)=FQC(L,1)-RIFTR(L)*CONQ(L,1)
            ENDIF
          ENDDO
        ENDIF

        ! *** TEMPERATURE ADJUSTMENTS FOR RAINFALL & EVAPORATION
        IF(M.EQ.2)THEN

          IF(ISTOPT(2).EQ.0.OR.ISTOPT(2).EQ.3)THEN
            DO L=2,LA
              FQC(L,KC)=FQC(L,KC)+RAINT(L)*TEMO*DXYP(L)
            ENDDO
          ENDIF
          IF(ISTOPT(2).EQ.1.OR.ISTOPT(2).EQ.2.OR.ISTOPT(2).EQ.4)THEN
            DO L=2,LA
              FQC(L,KC)=FQC(L,KC)+RAINT(L)*TATMT(L)*DXYP(L)
              FQCPAD(L,KC)=FQCPAD(L,KC)+RAINT(L)*TATMT(L)*DXYP(L)
              QSUMPAD(L,KC)=QSUMPAD(L,KC)+RAINT(L)*DXYP(L)
            ENDDO
          ENDIF
        ENDIF
        IF(M.EQ.2)THEN
          IF(ISTOPT(2).EQ.0)THEN
            DO L=2,LA
              FQC(L,KC)=FQC(L,KC)-EVAPSW(L)*CONQ(L,KC)
            ENDDO
          ENDIF
        ENDIF
      ENDIF
C
C *********************************************************************C
C
C *** 3TL CORRECTOR VOLUMETRICS SOURCE SINK LOCATIONS
C
      IF(ISTL_.EQ.2.AND.IS2TL_.EQ.0)THEN
        ! *** FLOW BOUNDARY CELLS (3TL CORRECTOR)
        DO NS=1,NQSIJ
          L=LQS(NS)
          NQSTMP=NQSERQ(NS)
          NCSTMP=NCSERQ(NS,M)
          DO K=1,KC
            FQC(L,K)=FQC(L,K)
     &          +MAX(QSS(K,NS),0.)*CQS(K,NS,M)
     &          +MIN(QSS(K,NS),0.)*CONQ(L,K)
     &          +MAX(QSERCELL(K,NS),0.)*CSERT(K,NCSTMP,M)
     &          +MIN(QSERCELL(K,NS),0.)*CONQ(L,K)
            FQCPAD(L,K)=FQCPAD(L,K)
     &          +MAX(QSS(K,NS),0.)*CQS(K,NS,M)
     &          +MAX(QSERCELL(K,NS),0.)*CSERT(K,NCSTMP,M)
            QSUMPAD(L,K)=QSUMPAD(L,K)
     &          +MAX(QSS(K,NS),0.)+MAX(QSERCELL(K,NS),0.)
C            QSUMNAD(L,K)=QSUMNAD(L,K)
C     &          +MIN(QSS(K,NS),0.)+MIN(QSERCELL(K,NS),0.)
          ENDDO
        ENDDO

        ! *** JET-PLUME VOLUMETRICS SOURCE SINK LOCATIONS (3TL CORRECTOR)
        IF(NQJPIJ.GT.0)THEN
          DO NJP=1,NQJPIJ
            IF(ICALJP(NJP).EQ.1)THEN
              RPORTS=FLOAT(NPORTJP(NJP))
              LJP=LIJ(IQJP(NJP),JQJP(NJP))
              KTMP=KEFFJP(NJP)

              ! QVJPTMP=TIME SERIES DISCHARGE FROM JET-PLUME
              QVJPTMP=0.
              DO K=1,KC
                QVJPTMP=QVJPTMP+QSERT(K,NQSERJP(NJP))
              ENDDO

              ! QCJPTMP=ENTRAINMENT FLUX
              QCJPTMP=0.
              QVJPENT=0.

              ! REMOVE ENTRAINMENT FLUX AND CALCULATE TOTAL ENTRAIMENT FLUX
              DO K=1,KC
                FQC(LJP,K)=FQC(LJP,K)-RPORTS*QJPENT(K,NJP)*CONQ(LJP,K)
                QCJPTMP=QCJPTMP+QJPENT(K,NJP)*CONQ(LJP,K)
                QVJPENT=QVJPENT+QJPENT(K,NJP)
C                QSUMNAD(LJP,K)=QSUMNAD(LJP,K)-RPORTS*QJPENT(K,NJP)
              ENDDO

              ! PLACE JET FLUX AND ENTRAINMENT FLUX IS EFFECTIVE LAYER
              FQC(LJP,KTMP)=FQC(LJP,KTMP)+RPORTS*QCJPTMP
     &            +RPORTS*QQCJP(NJP)*CQCJP(1,NJP,M)
     &            +RPORTS*QVJPTMP*CSERT(1,NCSERJP(NJP,M),M)
              FQCPAD(LJP,KTMP)=FQCPAD(LJP,KTMP)+RPORTS*QCJPTMP
     &            +RPORTS*QQCJP(NJP)*CQCJP(1,NJP,M)
     &            +RPORTS*QVJPTMP*CSERT(1,NCSERJP(NJP,M),M)
              QSUMPAD(LJP,KTMP)=QSUMPAD(LJP,KTMP)+RPORTS*QVJPENT
     &            +RPORTS*QQCJP(NJP)+RPORTS*QVJPTMP
            ENDIF
            IF(ICALJP(NJP).EQ.2)THEN
              RPORTS=FLOAT(NPORTJP(NJP))
              LJP=LIJ(IQJP(NJP),JQJP(NJP))
              KTMP=KEFFJP(NJP)
              NS=NQWRSERJP(NJP)
              LU=LIJ(IUPCJP(NJP),JUPCJP(NJP))
              KU=KUPCJP(NJP)
              CONUP=CONQ(LU,KU)
              QCJPTMP=0.
              QVJPENT=0.

              ! REMOVE ENTRAIMENT FLUX AND CALCULATE TOTAL ENTRAINMENT
              DO K=1,KC
                FQC(LJP,K)=FQC(LJP,K)-RPORTS*QJPENT(K,NJP)*CONQ(LJP,K)
                QCJPTMP=QCJPTMP+QJPENT(K,NJP)*CONQ(LJP,K)
                QVJPENT=QVJPENT+QJPENT(K,NJP)
C                QSUMNAD(LJP,K)=QSUMNAD(LJP,K)-RPORTS*QJPENT(K,NJP)
              ENDDO

              !  PLACE ENTRAINMENT, CONSTANT AND TIME SERIES FLUXES IN EFFECTIVE CELL
              FQC(LJP,KTMP)=FQC(LJP,KTMP)+RPORTS*QCJPTMP
     &            +RPORTS*QWRCJP(NJP)*(CWRCJP(NJP,M)+CONUP)
     &            +RPORTS*QWRSERT(NS)*(CQWRSERT(NS,M)+CONUP)
              FQCPAD(LJP,KTMP)=FQCPAD(LJP,KTMP)+RPORTS*QCJPTMP
     &            +RPORTS*QWRCJP(NJP)*(CWRCJP(NJP,M)+CONUP)
     &            +RPORTS*QWRSERT(NS)*(CQWRSERT(NS,M)+CONUP)
              QSUMPAD(LJP,KTMP)=QSUMPAD(LJP,KTMP)+RPORTS*QVJPENT
     &            +RPORTS*QWRCJP(NJP)+RPORTS*QWRSERT(NS)

              ! REMOVAL WITHDRAWAL FROM UPSTREAM CELL
              FQC(LU,KU)=FQC(LU,KU)
     &            -RPORTS*QWRCJP(NJP)*CONUP
     &            -RPORTS*QWRSERT(NS)*CONUP
C              QSUMNAD(LU,KU)=QSUMNAD(LU,KU)
C     &            -RPORTS*QWRCJP(NJP)-RPORTS*QWRSERT(NS)
            ENDIF
          ENDDO
        ENDIF

        ! *** CONTROL STRUCTURES (3TL CORRECTOR)
        DO NCTL=1,NQCTL
          RQWD=1.
          IU=IQCTLU(NCTL)
          JU=JQCTLU(NCTL)
          LU=LIJ(IU,JU)
          ID=IQCTLD(NCTL)
          JD=JQCTLD(NCTL)
          IF(ID.EQ.0.AND.JD.EQ.0)THEN
            LD=LC
            RQWD=0.
          ELSE
            LD=LIJ(ID,JD)
          ENDIF
          DO K=1,KC
            FQC(LU,K)=FQC(LU,K)
     &          -QCTLT(K,NCTL)*CONQ(LU,K)
            FQC(LD,K)=FQC(LD,K)
     &          +RQWD*QCTLT(K,NCTL)*CONQ(LU,K)
            FQCPAD(LD,K)=FQCPAD(LD,K)
     &          +RQWD*QCTLT(K,NCTL)*CONQ(LU,K)
            QSUMPAD(L,K)=QSUMPAD(L,K)
     &          +RQWD*QCTLT(K,NCTL)
C            QSUMNAD(L,K)=QSUMNAD(L,K)
C     &          -QCTLT(K,NCTL)
          ENDDO
        ENDDO

        ! ***  WITHDRAWAL CONCENTRATION RISE RETURN (3TL CORRECTOR)
        DO NWR=1,NQWR
          IU=IQWRU(NWR)
          JU=JQWRU(NWR)
          KU=KQWRU(NWR)
          ID=IQWRD(NWR)
          JD=JQWRD(NWR)
          KD=KQWRD(NWR)
          LU=LIJ(IU,JU)
          LD=LIJ(ID,JD)
          NQSTMP=NQWRSERQ(NWR)
          NCSTMP=NQWRSERQ(NWR)
          FQC(LU,KU)=FQC(LU,KU)
     &        -(QWR(NWR)+QWRSERT(NQSTMP))*CONQ(LU,KU)
          FQC(LD,KD)=FQC(LD,KD)
     &        +QWR(NWR)*(CONQ(LU,KU)+CQWR(NWR,M))
     &        +QWRSERT(NQSTMP)*(CONQ(LU,KU)+CQWRSERT(NCSTMP,M))
          FQCPAD(LD,KD)=FQCPAD(LD,KD)
     &        +QWR(NWR)*(CONQ(LU,KU)+CQWR(NWR,M))
     &        +QWRSERT(NQSTMP)*(CONQ(LU,KU)+CQWRSERT(NCSTMP,M))
          QSUMPAD(LD,KD)=QSUMPAD(LD,KD)
     &        +QWR(NWR)+QWRSERT(NQSTMP)
C          QSUMNAD(LU,KU)=QSUMNAD(LU,KU)
C     &        -(QWR(NWR)+QWRSERT(NQSTMP))
        ENDDO

        ! ***  SUBGRID SCALE CHANNEL EXCHANGE (3TL CORRECTOR)
        IF(MDCHH.GE.1)THEN
          DO K=1,KC
            DO NMD=1,MDCHH
              LMDCHHT=LMDCHH(NMD)
              LMDCHUT=LMDCHU(NMD)
              LMDCHVT=LMDCHV(NMD)
              IF(MDCHTYP(NMD).EQ.1)THEN
                QUKTMP=QCHANU(NMD)*DZC(K)
                QVKTMP=0.
              ENDIF
              IF(MDCHTYP(NMD).EQ.2)THEN
                QVKTMP=QCHANV(NMD)*DZC(K)
                QUKTMP=0.
              ENDIF
              IF(MDCHTYP(NMD).EQ.3)THEN
                QUKTMP=QCHANU(NMD)*DZC(K)
                QVKTMP=QCHANV(NMD)*DZC(K)
              ENDIF
              FQC(LMDCHHT,K)=FQC(LMDCHHT,K)
     &            +MAX(QUKTMP,0.)*CONQ(LMDCHUT,K)
     &            +MIN(QUKTMP,0.)*CONQ(LMDCHHT,K)
     &            +MAX(QVKTMP,0.)*CONQ(LMDCHVT,K)
     &            +MIN(QVKTMP,0.)*CONQ(LMDCHHT,K)
              FQC(LMDCHUT,K)=FQC(LMDCHUT,K)
     &            -MAX(QUKTMP,0.)*CONQ(LMDCHUT,K)
     &            -MIN(QUKTMP,0.)*CONQ(LMDCHHT,K)
              FQC(LMDCHVT,K)=FQC(LMDCHVT,K)
     &            -MAX(QVKTMP,0.)*CONQ(LMDCHVT,K)
     &            -MIN(QVKTMP,0.)*CONQ(LMDCHHT,K)
            ENDDO
          ENDDO
        ENDIF

        ! ***  GROUNDWATER, EVAP, RAINFALL (3TL CORRECTOR)
        IF(ISGWIE.NE.0)THEN
          DO L=2,LA
            FQC(L,1)=FQC(L,1)-RIFTR(L)*CONQ(L,1)
          ENDDO
        ENDIF

        ! *** ZONED SEEPAGE (3TL)
        IF(ISGWIT.EQ.3)THEN
          DO L=2,LA
            IF(H1P(L).GT.HDRY)THEN
              FQC(L,1)=FQC(L,1)-RIFTR(L)*CONQ(L,1)
            ENDIF
          ENDDO
        ENDIF

        ! *** TEMPERATURE ADJUSTMENTS FOR RAINFALL & EVAPORATION
        IF(M.EQ.2)THEN
          IF(ISTOPT(2).EQ.0.OR.ISTOPT(2).EQ.3)THEN
            DO L=2,LA
              FQC(L,KC)=FQC(L,KC)+RAINT(L)*TEMO*DXYP(L)
            ENDDO
          ENDIF
          IF(ISTOPT(2).EQ.1.OR.ISTOPT(2).EQ.2.OR.ISTOPT(2).EQ.4)THEN
            DO L=2,LA
              FQC(L,KC)=FQC(L,KC)+RAINT(L)*TATMT(L)*DXYP(L)
              FQCPAD(L,KC)=FQCPAD(L,KC)+RAINT(L)*TATMT(L)*DXYP(L)
              QSUMPAD(L,KC)=QSUMPAD(L,KC)+RAINT(L)*DXYP(L)
            ENDDO
          ENDIF
        ENDIF
        IF(M.EQ.2)THEN
          IF(ISTOPT(2).EQ.0)THEN
            DO L=2,LA
              FQC(L,KC)=FQC(L,KC)-EVAPSW(L)*CONQ(L,KC)
            ENDDO
          ENDIF
        ENDIF
      ENDIF
C *********************************************************************C
C
C **  STANDARD VOLUMETRICS SOURCE SINK LOCATIONS (3TL)
C
      IF(ISTL_.EQ.3)THEN
        ! *** FLOW BOUNDARY CELLS (3TL)
        DO NS=1,NQSIJ
          L=LQS(NS)
          NQSTMP=NQSERQ(NS)
          NCSTMP=NCSERQ(NS,M)
          DO K=1,KC
            FQC(L,K)=FQC(L,K)
     &          +MAX(QSS(K,NS),0.)*CQS(K,NS,M)
     &          +MIN(QSS(K,NS),0.)*CONQ(L,K)
     &          +MAX(QSERCELL(K,NS),0.)*CSERT(K,NCSTMP,M)
     &          +MIN(QSERCELL(K,NS),0.)*CONQ(L,K)
            FQCPAD(L,K)=FQCPAD(L,K)
     &          +MAX(QSS(K,NS),0.)*CQS(K,NS,M)
     &          +MAX(QSERCELL(K,NS),0.)*CSERT(K,NCSTMP,M)
            QSUMPAD(L,K)=QSUMPAD(L,K)
     &          +MAX(QSS(K,NS),0.)+MAX(QSERCELL(K,NS),0.)
C            QSUMNAD(L,K)=QSUMNAD(L,K)
C     &          +MIN(QSS(K,NS),0.)+MIN(QSERCELL(K,NS),0.)
          ENDDO
        ENDDO

        ! ***  JET-PLUME VOLUMETRICS SOURCE SINK LOCATIONS (3TL)
        IF(NQJPIJ.GT.0)THEN
          DO NJP=1,NQJPIJ
            IF(ICALJP(NJP).EQ.1)THEN
              RPORTS=FLOAT(NPORTJP(NJP))
              LJP=LIJ(IQJP(NJP),JQJP(NJP))
              KTMP=KEFFJP(NJP)

              ! QVJPTMP=TIME SERIES DISCHARGE FROM JET-PLUME
              QVJPTMP=0.
              DO K=1,KC
                QVJPTMP=QVJPTMP+QSERT(K,NQSERJP(NJP))
              ENDDO

              ! QCJPTMP=ENTRAINMENT FLUX
              QCJPTMP=0.
              QVJPENT=0.

              ! REMOVE ENTRAINMENT FLUX AND CALCULATE TOTAL ENTRAIMENT FLUX
              DO K=1,KC
                FQC(LJP,K)=FQC(LJP,K)-RPORTS*QJPENT(K,NJP)*CONQ(LJP,K)
                QCJPTMP=QCJPTMP+QJPENT(K,NJP)*CONQ(LJP,K)
                QVJPENT=QVJPENT+QJPENT(K,NJP)
C                QSUMNAD(LJP,K)=QSUMNAD(LJP,K)-RPORTS*QJPENT(K,NJP)
              ENDDO

              ! PLACE JET FLUX AND ENTRAINMENT FLUX IS EFFECTIVE LAYER
              FQC(LJP,KTMP)=FQC(LJP,KTMP)+RPORTS*QCJPTMP
     &            +RPORTS*QQCJP(NJP)*CQCJP(1,NJP,M)
     &            +RPORTS*QVJPTMP*CSERT(1,NCSERJP(NJP,M),M)
              FQCPAD(LJP,KTMP)=FQCPAD(LJP,KTMP)+RPORTS*QCJPTMP
     &            +RPORTS*QQCJP(NJP)*CQCJP(1,NJP,M)
     &            +RPORTS*QVJPTMP*CSERT(1,NCSERJP(NJP,M),M)
              QSUMPAD(LJP,KTMP)=QSUMPAD(LJP,KTMP)+RPORTS*QVJPENT
     &            +RPORTS*QQCJP(NJP)+RPORTS*QVJPTMP
            ENDIF
            IF(ICALJP(NJP).EQ.2)THEN
              RPORTS=FLOAT(NPORTJP(NJP))
              LJP=LIJ(IQJP(NJP),JQJP(NJP))
              KTMP=KEFFJP(NJP)
              NS=NQWRSERJP(NJP)
              LU=LIJ(IUPCJP(NJP),JUPCJP(NJP))
              KU=KUPCJP(NJP)
              CONUP=CONQ(LU,KU)
              QCJPTMP=0.
              QVJPENT=0.

              ! REMOVE ENTRAIMENT FLUX AND CALCULATE TOTAL ENTRAINMENT
              DO K=1,KC
                FQC(LJP,K)=FQC(LJP,K)-RPORTS*QJPENT(K,NJP)*CONQ(LJP,K)
                QCJPTMP=QCJPTMP+QJPENT(K,NJP)*CONQ(LJP,K)
                QVJPENT=QVJPENT+QJPENT(K,NJP)
C                QSUMNAD(LJP,K)=QSUMNAD(LJP,K)-RPORTS*QJPENT(K,NJP)
              ENDDO

              ! PLACE ENTRAINMENT, CONSTANT AND TIME SERIES FLUXES IN EFFECTIVE CELL
              FQC(LJP,KTMP)=FQC(LJP,KTMP)+RPORTS*QCJPTMP
     &            +RPORTS*QWRCJP(NJP)*(CWRCJP(NJP,M)+CONUP)
     &            +RPORTS*QWRSERT(NS)*(CQWRSERT(NS,M)+CONUP)
              FQCPAD(LJP,KTMP)=FQCPAD(LJP,KTMP)+RPORTS*QCJPTMP
     &            +RPORTS*QWRCJP(NJP)*(CWRCJP(NJP,M)+CONUP)
     &            +RPORTS*QWRSERT(NS)*(CQWRSERT(NS,M)+CONUP)
              QSUMPAD(LJP,KTMP)=QSUMPAD(LJP,KTMP)+RPORTS*QVJPENT
     &            +RPORTS*QWRCJP(NJP)+RPORTS*QWRSERT(NS)

              ! REMOVAL WITHDRAWAL FROM UPSTREAM CELL
              FQC(LU,KU)=FQC(LU,KU)
     &            -RPORTS*QWRCJP(NJP)*CONUP
     &            -RPORTS*QWRSERT(NS)*CONUP
C              QSUMNAD(LU,KU)=QSUMNAD(LU,KU)
C     &            -RPORTS*QWRCJP(NJP)-RPORTS*QWRSERT(NS)
            ENDIF
          ENDDO
        ENDIF

        ! ***  CONTROL STRUCTURES (3TL)
        DO NCTL=1,NQCTL
          RQWD=1.
          IU=IQCTLU(NCTL)
          JU=JQCTLU(NCTL)
          LU=LIJ(IU,JU)
          ID=IQCTLD(NCTL)
          JD=JQCTLD(NCTL)
          IF(ID.EQ.0.AND.JD.EQ.0)THEN
            LD=LC
            RQWD=0.
          ELSE
            LD=LIJ(ID,JD)
          ENDIF
          DO K=1,KC
            FQC(LU,K)=FQC(LU,K)
     &          -QCTLT(K,NCTL)*CONQ(LU,K)
            FQC(LD,K)=FQC(LD,K)
     &          +RQWD*QCTLT(K,NCTL)*CONQ(LU,K)
            FQCPAD(LD,K)=FQCPAD(LD,K)
     &          +RQWD*QCTLT(K,NCTL)*CONQ(LU,K)
            QSUMPAD(L,K)=QSUMPAD(L,K)
     &          +RQWD*QCTLT(K,NCTL)
C            QSUMNAD(L,K)=QSUMNAD(L,K)
C     &          -QCTLT(K,NCTL)
          ENDDO
        ENDDO

        ! ***  WITHDRAWAL CONCENTRATION RISE RETURN (3TL)
        DO NWR=1,NQWR
          IU=IQWRU(NWR)
          JU=JQWRU(NWR)
          KU=KQWRU(NWR)
          ID=IQWRD(NWR)
          JD=JQWRD(NWR)
          KD=KQWRD(NWR)
          LU=LIJ(IU,JU)
          LD=LIJ(ID,JD)
          NQSTMP=NQWRSERQ(NWR)
          NCSTMP=NQWRSERQ(NWR)
          FQC(LU,KU)=FQC(LU,KU)
     &        -(QWR(NWR)+QWRSERT(NQSTMP))*CONQ(LU,KU)
          FQC(LD,KD)=FQC(LD,KD)
     &        +QWR(NWR)*(CONQ(LU,KU)+CQWR(NWR,M))
     &        +QWRSERT(NQSTMP)*(CONQ(LU,KU)+CQWRSERT(NCSTMP,M))
          FQCPAD(LD,KD)=FQCPAD(LD,KD)
     &        +QWR(NWR)*(CONQ(LU,KU)+CQWR(NWR,M))
     &        +QWRSERT(NQSTMP)*(CONQ(LU,KU)+CQWRSERT(NCSTMP,M))
          QSUMPAD(LD,KD)=QSUMPAD(LD,KD)
     &        +QWR(NWR)+QWRSERT(NQSTMP)
C          QSUMNAD(LU,KU)=QSUMNAD(LU,KU)
C     &        -(QWR(NWR)+QWRSERT(NQSTMP))
        ENDDO

        ! ***  SUBGRID SCALE CHANNEL EXCHANGE (3TL)
        IF(MDCHH.GE.1)THEN
          DO K=1,KC
            DO NMD=1,MDCHH
              LMDCHHT=LMDCHH(NMD)
              LMDCHUT=LMDCHU(NMD)
              LMDCHVT=LMDCHV(NMD)
              IF(MDCHTYP(NMD).EQ.1)THEN
                QUKTMP=QCHANU(NMD)*DZC(K)
                QVKTMP=0.
              ENDIF
              IF(MDCHTYP(NMD).EQ.2)THEN
                QVKTMP=QCHANV(NMD)*DZC(K)
                QUKTMP=0.
              ENDIF
              IF(MDCHTYP(NMD).EQ.3)THEN
                QUKTMP=QCHANU(NMD)*DZC(K)
                QVKTMP=QCHANV(NMD)*DZC(K)
              ENDIF
              FQC(LMDCHHT,K)=FQC(LMDCHHT,K)
     &            +MAX(QUKTMP,0.)*CONQ(LMDCHUT,K)
     &            +MIN(QUKTMP,0.)*CONQ(LMDCHHT,K)
     &            +MAX(QVKTMP,0.)*CONQ(LMDCHVT,K)
     &            +MIN(QVKTMP,0.)*CONQ(LMDCHHT,K)
              FQC(LMDCHUT,K)=FQC(LMDCHUT,K)
     &            -MAX(QUKTMP,0.)*CONQ(LMDCHUT,K)
     &            -MIN(QUKTMP,0.)*CONQ(LMDCHHT,K)
              FQC(LMDCHVT,K)=FQC(LMDCHVT,K)
     &            -MAX(QVKTMP,0.)*CONQ(LMDCHVT,K)
     &            -MIN(QVKTMP,0.)*CONQ(LMDCHHT,K)
            ENDDO
          ENDDO
        ENDIF

        ! ***  GROUNDWATER, EVAP, RAINFALL (3TL)
        IF(ISGWIE.NE.0)THEN
          DO L=2,LA
            FQC(L,1)=FQC(L,1)-RIFTR(L)*CONQ(L,1)
          ENDDO
        ENDIF

        ! *** ZONED SEEPAGE (3TL)
        IF(ISGWIT.EQ.3)THEN
          DO L=2,LA
            IF(H1P(L).GT.HDRY)THEN
              FQC(L,1)=FQC(L,1)-RIFTR(L)*CONQ(L,1)
            ENDIF
          ENDDO
        ENDIF

        ! *** TEMPERATURE ADJUSTMENTS FOR RAINFALL & EVAPORATION
        IF(M.EQ.2)THEN
          IF(ISTOPT(2).EQ.0.OR.ISTOPT(2).EQ.3)THEN
            DO L=2,LA
              FQC(L,KC)=FQC(L,KC)+RAINT(L)*TEMO*DXYP(L)
            ENDDO
          ENDIF
          IF(ISTOPT(2).EQ.1.OR.ISTOPT(2).EQ.2.OR.ISTOPT(2).EQ.4)THEN
            DO L=2,LA
              FQC(L,KC)=FQC(L,KC)+RAINT(L)*TATMT(L)*DXYP(L)
              FQCPAD(L,KC)=FQCPAD(L,KC)+RAINT(L)*TATMT(L)*DXYP(L)
              QSUMPAD(L,KC)=QSUMPAD(L,KC)+RAINT(L)*DXYP(L)
            ENDDO
          ENDIF
        ENDIF
        IF(M.EQ.2)THEN
          IF(ISTOPT(2).EQ.0)THEN
            DO L=2,LA
              FQC(L,KC)=FQC(L,KC)-EVAPSW(L)*CONQ(L,KC)
            ENDDO
          ENDIF
        ENDIF
      ENDIF
      GOTO 2000
C
C *** SHELL FISH LARVAE SECTION
C
 1000 CONTINUE
      DO NS=1,NQSIJ
        L=LQS(NS)
        NQSTMP=NQSERQ(NS)
        NCSTMP=NCSERQ(NS,M)
        DO K=1,KC
          FQC(L,K)=FQC(L,K)
     &        +MAX(QSS(K,NS),0.)*CQS(K,NS,M)
     &        +MIN(QSS(K,NS),0.)*CONQ(L,K)
     &        +MAX(QSERCELL(K,NS),0.)*CSERT(K,NCSTMP,M)
     &        +MIN(QSERCELL(K,NS),0.)*CONQ(L,K)
        ENDDO
      ENDDO
      DO NCTL=1,NQCTL
        RQWD=1.
        IU=IQCTLU(NCTL)
        JU=JQCTLU(NCTL)
        LU=LIJ(IU,JU)
        ID=IQCTLD(NCTL)
        JD=JQCTLD(NCTL)
        IF(ID.EQ.0.AND.JD.EQ.0)THEN
          LD=LC
          RQWD=0.
        ELSE
          LD=LIJ(ID,JD)
        ENDIF
        DO K=1,KC
          FQC(LU,K)=FQC(LU,K)
     &        -QCTLT(K,NCTL)*CONQ(LU,K)
          FQC(LD,K)=FQC(LD,K)
     &        +RQWD*QCTLT(K,NCTL)*CONQ(LU,K)
        ENDDO
      ENDDO
      DO NWR=1,NQWR
        IU=IQWRU(NWR)
        JU=JQWRU(NWR)
        KU=KQWRU(NWR)
        ID=IQWRD(NWR)
        JD=JQWRD(NWR)
        KD=KQWRD(NWR)
        LU=LIJ(IU,JU)
        LD=LIJ(ID,JD)
        NQSTMP=NQWRSERQ(NWR)
        NCSTMP=NQWRSERQ(NWR)
        FQC(LU,KU)=FQC(LU,KU)
     &      -(QWR(NWR)+QWRSERT(NQSTMP))*CONQ(LU,KU)
        FQC(LD,KD)=FQC(LD,KD)
     &      +QWR(NWR)*(SFLKILL*CONQ(LU,KU)+CQWR(NWR,M))
     &      +QWRSERT(NQSTMP)*(SFLKILL*CONQ(LU,KU)
     &      +CQWRSERT(NCSTMP,M))
      ENDDO
      IF(MDCHH.GE.1)THEN
        DO K=1,KC
          DO NMD=1,MDCHH
            LMDCHHT=LMDCHH(NMD)
            LMDCHUT=LMDCHU(NMD)
            LMDCHVT=LMDCHV(NMD)
            IF(MDCHTYP(NMD).EQ.1)THEN
              QUKTMP=QCHANU(NMD)*DZC(K)
              QVKTMP=0.
            ENDIF
            IF(MDCHTYP(NMD).EQ.2)THEN
              QVKTMP=QCHANV(NMD)*DZC(K)
              QUKTMP=0.
            ENDIF
            IF(MDCHTYP(NMD).EQ.3)THEN
              QUKTMP=QCHANU(NMD)*DZC(K)
              QVKTMP=QCHANV(NMD)*DZC(K)
            ENDIF
            FQC(LMDCHHT,K)=FQC(LMDCHHT,K)
     &          +MAX(QUKTMP,0.)*CONQ(LMDCHUT,K)
     &          +MIN(QUKTMP,0.)*CONQ(LMDCHHT,K)
     &          +MAX(QVKTMP,0.)*CONQ(LMDCHVT,K)
     &          +MIN(QVKTMP,0.)*CONQ(LMDCHHT,K)
            FQC(LMDCHUT,K)=FQC(LMDCHUT,K)
     &          -MAX(QUKTMP,0.)*CONQ(LMDCHUT,K)
     &          -MIN(QUKTMP,0.)*CONQ(LMDCHHT,K)
            FQC(LMDCHVT,K)=FQC(LMDCHVT,K)
     &          -MAX(QVKTMP,0.)*CONQ(LMDCHVT,K)
     &          -MIN(QVKTMP,0.)*CONQ(LMDCHHT,K)
          ENDDO
        ENDDO
      ENDIF
      GOTO 2000
C
C *** WATER QUALITY ONLY (IWQPSL=1,0)
C
 1500 CONTINUE

      ! *** FLOW BOUNDARIES
      DO NS=1,NQSIJ
        L=LQS(NS)
        NQSTMP=NQSERQ(NS)
        DO K=1,KC
          FQC(L,K)=FQC(L,K)
     &        +MIN(QSS(K,NS),0.)*CON1(L,K)
     &        +MIN(QSERCELL(K,NS),0.)*CON1(L,K)
          QSUMPAD(L,K)=QSUMPAD(L,K)
     &        +MAX(QSS(K,NS),0.)+MAX(QSERCELL(K,NS),0.)
C          QSUMNAD(L,K)=QSUMNAD(L,K)
C     &        +MIN(QSS(K,NS),0.)+MIN(QSERCELL(K,NS),0.)
        ENDDO
      ENDDO

      ! *** HYDRAULIC STRUCTURES
      DO NCTL=1,NQCTL
        RQWD=1.
        IU=IQCTLU(NCTL)
        JU=JQCTLU(NCTL)
        LU=LIJ(IU,JU)
        ID=IQCTLD(NCTL)
        JD=JQCTLD(NCTL)
        IF(ID.EQ.0.AND.JD.EQ.0)THEN
          LD=LC
          RQWD=0.
        ELSE
          LD=LIJ(ID,JD)
        ENDIF
        DO K=1,KC
          FQC(LU,K)=FQC(LU,K)
     &        -QCTLT(K,NCTL)*CON1(LU,K)
          FQC(LD,K)=FQC(LD,K)
     &        +RQWD*QCTLT(K,NCTL)*CON1(LU,K)
          FQCPAD(LD,K)=FQCPAD(LD,K)
     &        +RQWD*QCTLT(K,NCTL)*CON1(LU,K)
          QSUMPAD(L,K)=QSUMPAD(L,K)
     &        +RQWD*QCTLT(K,NCTL)
C          QSUMNAD(L,K)=QSUMNAD(L,K)
C     &        -QCTLT(K,NCTL)
        ENDDO
      ENDDO

      ! *** WITHDRAWAL/RETURN
      IF(MVAR.EQ.8)THEN
        M=4+NTOX+NSED+NSND+MO
      ELSE
        M=MO
      ENDIF
      DO NWR=1,NQWR
        IU=IQWRU(NWR)
        JU=JQWRU(NWR)
        KU=KQWRU(NWR)
        ID=IQWRD(NWR)
        JD=JQWRD(NWR)
        KD=KQWRD(NWR)
        LU=LIJ(IU,JU)
        LD=LIJ(ID,JD)
        NQSTMP=NQWRSERQ(NWR)
        NCSTMP=NQWRSERQ(NWR)
        FQC(LU,KU)=FQC(LU,KU)
     &      -(QWR(NWR)+QWRSERT(NQSTMP))*CON1(LU,KU)
        FQC(LD,KD)=FQC(LD,KD)
     &      +QWR(NWR)*(CON1(LU,KU)+CQWR(NWR,M))
     &      +QWRSERT(NQSTMP)*(CON1(LU,KU)+CQWRSERT(NCSTMP,M))
        FQCPAD(LD,KD)=FQCPAD(LD,KD)
     &      +QWR(NWR)*(CON1(LU,KU)+CQWR(NWR,M))
     &      +QWRSERT(NQSTMP)*(CON1(LU,KU)+CQWRSERT(NCSTMP,M))
        QSUMPAD(LD,KD)=QSUMPAD(LD,KD)
     &      +QWR(NWR)+QWRSERT(NQSTMP)
C        QSUMNAD(LU,KU)=QSUMNAD(LU,KU)
C     &      -(QWR(NWR)+QWRSERT(NQSTMP))
      ENDDO

      ! *** SUBGRID CHANNEL INTERACTIONS
      IF(MDCHH.GE.1)THEN
        DO K=1,KC
          DO NMD=1,MDCHH
            LMDCHHT=LMDCHH(NMD)
            LMDCHUT=LMDCHU(NMD)
            LMDCHVT=LMDCHV(NMD)
            IF(MDCHTYP(NMD).EQ.1)THEN
              QUKTMP=QCHANU(NMD)*DZC(K)
              QVKTMP=0.
            ENDIF
            IF(MDCHTYP(NMD).EQ.2)THEN
              QVKTMP=QCHANV(NMD)*DZC(K)
              QUKTMP=0.
            ENDIF
            IF(MDCHTYP(NMD).EQ.3)THEN
              QUKTMP=QCHANU(NMD)*DZC(K)
              QVKTMP=QCHANV(NMD)*DZC(K)
            ENDIF
            FQC(LMDCHHT,K)=FQC(LMDCHHT,K)
     &          +MAX(QUKTMP,0.)*CON1(LMDCHUT,K)
     &          +MIN(QUKTMP,0.)*CON1(LMDCHHT,K)
     &          +MAX(QVKTMP,0.)*CON1(LMDCHVT,K)
     &          +MIN(QVKTMP,0.)*CON1(LMDCHHT,K)
            FQC(LMDCHUT,K)=FQC(LMDCHUT,K)
     &          -MAX(QUKTMP,0.)*CON1(LMDCHUT,K)
     &          -MIN(QUKTMP,0.)*CON1(LMDCHHT,K)
            FQC(LMDCHVT,K)=FQC(LMDCHVT,K)
     &          -MAX(QVKTMP,0.)*CON1(LMDCHVT,K)
     &          -MIN(QVKTMP,0.)*CON1(LMDCHHT,K)
          ENDDO
        ENDDO
      ENDIF

      ! *** GROUNDWATER, EVAP, RAINFALL (2TL)
      IF(ISGWIE.NE.0)THEN
        DO L=2,LA
          FQC(L,1)=FQC(L,1)-RIFTR(L)*CONQ(L,1)
        ENDDO
      ENDIF

      ! *** ZONED SEEPAGE (2TL)
      IF(ISGWIT.EQ.3)THEN
        DO L=2,LA
          IF(H1P(L).GT.HDRY)THEN
            FQC(L,1)=FQC(L,1)-RIFTR(L)*CON1(L,1)
          ENDIF
        ENDDO
      ENDIF

 2000 CONTINUE
      RETURN
      END


text:
      SUBROUTINE RANSTArefine ( N, X, FU, EPS, F0, IER,cotet,
     &     kontet,ipoeln,ieln,node,iedge,
     &     ipoeled,ieled,iedgmid,iedtet )
C
C Random search to determine a point with a function value that is
C different from the one found at the starting point X
C This is an auxiliary subroutine to support the minimization
C subroutine FMINSI
C
C   FMINSI - Fortran subroutines for unconstrained function minimization
C   Copyright (C) 1986, 1993, 2001  Hugo Pfoertner
C
C   This library is free software; you can redistribute it and/or
C   modify it under the terms of the GNU Lesser General Public
C   License as published by the Free Software Foundation; either
C   version 2.1 of the License, or (at your option) any later version.
C
C   This library is distributed in the hope that it will be useful,
C   but WITHOUT ANY WARRANTY; without even the implied warranty of
C   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
C   Lesser General Public License for more details.
C
C   You should have received a copy of the GNU Lesser General Public
C   License along with this library; if not, write to the Free Software
C   Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307
C   USA
C
C   Contact info: mailto:hugo@pfoertner.org
C   or use the information provided at http://www.pfoertner.org/
C
C Author: Hugo Pfoertner, Oberhaching, Germany
C
C Version History (left in German):
! 22.06.21 adapted fminsi for the mesh refinement procedure in CalculiX
!          Author: Guido Dhondt
C 02.06.01 English translation of comments, LGPL header added
C 27.11.93 ZUFALLSZAHLENGENERATOR RANEWR UND INIRAN EINGEBAUT
C 16.11.92 DIMENSIONIERUNG AUF 128 ERHOEHT
C 15.01.91 ABFRAGE AUF AENDERUNG STATT AUF VERMINDERUNG DES
C          FUNKTIONSWERTES
C 09.04.86 BASISVERSION
C
C The meaning of the parameters is the same as in subroutine FMINSI
      implicit doubleprecision (a-h,o-z)
      INTEGER N, IER,kontet(4,*),ipoeln(*),ieln(2,*),node,iedge,
     &     ipoeled(*),ieled(2,*),iedgmid(*),iedtet(6,*)
      doubleprecision FU, F0, X(N), EPS(N),cotet(3,*)
      EXTERNAL FU
C Output:
C F0 ...  Changed function value
C X ...   Variable vector of point with different function value
C         if such a point has been found, unchanged otherwise
C IER ... 0, if a point with changed function value has been found
C         3 otherwise (search stopped without success after 100*N
C         function evaluations.
C
C Dimension of local arrays has to be compliant with corresponding
C size within FMINSI
      PARAMETER ( NMAX=216 )
      doubleprecision RS(NMAX), XS(NMAX)
C Uniform random number generator and corresponding initialization
      REAL RANUWH
      EXTERNAL INIRAN, RANUWH
C
C
C Preset return code with "No success"
      IER = 3
C
C Initialize uniform random number generator
C
      CALL INIRAN
C
C Set initial step size
      DO 10 K = 1, N
      RS(K) = 10. * EPS(K)
10    CONTINUE
C
C After N2 function evaluations, an increase in variance is tried
      N2 = 2**N
C
C Maximum number of function evaluations
      NT = 100 * N
C
C Function value at starting point
      F0 = FU ( N, X,cotet,kontet,ipoeln,ieln,node,iedge,
     &     ipoeled,ieled,iedgmid,iedtet )
      LL = 0
C
C Loop over maximum number of trials
      DO 20 L = 1, NT
      DO 30 NN = 1, N
C
C Set co-ordinates to starting point + random increment
C in the range +-RS(i)
      XS(NN) = X(NN) + RS(NN) * 2. * ( dble(RANUWH()) - 0.5 )
30    CONTINUE
C
C Corresponding function value
      H = FU ( N, XS,cotet,kontet,ipoeln,ieln,node,iedge,
     &     ipoeled,ieled,iedgmid,iedtet )
C
C Check for change
      IF ( H .NE. F0 ) THEN
C
C Loop is terminated at first occurrence of changed value
        F0 = H
        DO 40 K = 1, N
40      X(K) = XS(K)
        IER = 0
        GOTO 999
      ENDIF
C
C After N2 trials without success, the variance is increased
C
      LL = LL + 1
      IF ( LL .GE. N2 ) THEN
        DO 50 K = 1, N
50      RS(K) = RS(K) + RS(K)
        LL = 0
      ENDIF
C
C End of loop over maximum number of trials
20    CONTINUE
C
999   CONTINUE
      RETURN
C End of subroutine RANSTA
      END

max_stars_repo_name: alleindrach/calculix-desktop

max_issues_repo_name: alleindrach/calculix-desktop

max_forks_repo_name: alleindrach/calculix-desktop


text:
subroutine corrht(hrm       ,deph      ,tp        ,wavel     ,wavek     , &
                & ldep      ,dish      ,dismax    ,choice    ,rho       , &
                & grav      )
!----- GPL ---------------------------------------------------------------------
!
!  Copyright (C)  Stichting Deltares, 2011-2016.
!
!  This program is free software: you can redistribute it and/or modify
!  it under the terms of the GNU General Public License as published by
!  the Free Software Foundation version 3.
!
!  This program is distributed in the hope that it will be useful,
!  but WITHOUT ANY WARRANTY; without even the implied warranty of
!  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!  GNU General Public License for more details.
!
!  You should have received a copy of the GNU General Public License
!  along with this program.  If not, see <http://www.gnu.org/licenses/>.
!
!  contact: delft3d.support@deltares.nl
!  Stichting Deltares
!  P.O. Box 177
!  2600 MH Delft, The Netherlands
!
!  All indications and logos of, and references to, "Delft3D" and "Deltares"
!  are registered trademarks of Stichting Deltares, and remain the property of
!  Stichting Deltares. All rights reserved.
!
!-------------------------------------------------------------------------------
!  $Id: corrht.f90 5717 2016-01-12 11:35:24Z mourits $
!  $HeadURL: https://svn.oss.deltares.nl/repos/delft3d/tags/6686/src/engines_gpl/wave/packages/kernel/src/corrht.f90 $
!!--description-----------------------------------------------------------------
!
!
!     Input:
!     -------
!     HSH,DEPH,TP
!     I/O  I    I/O
!
!     Output:
!     --------
!     WAVEL,WAVEK,LDEP
!     LDEP  : logical variable, .true. when depth or wave height too small
!
!     Adapt wave period TP and wave height HSH
!
!!--pseudo code and references--------------------------------------------------
! NONE
!!--declarations----------------------------------------------------------------
    implicit none
    !
    !
    ! COMMON variables
    !
    real ::         pi, twopi, wort2, gamma
    common /const / pi, twopi, wort2, gamma
!
! Global variables
!
    logical, intent(in)            :: choice
    logical, intent(out)           :: ldep
    real                           :: deph
    real                           :: dish
    real                           :: dismax
    real                           :: grav
    real                           :: hrm
    real   , intent(in)            :: rho
    real                           :: tp
    real                           :: wavek
    real   , intent(out)           :: wavel
!
! Local variables
!
    real :: hmax
    real :: hs
    real :: tpmin
!
!! executable statements -------------------------------------------------------
!
    ldep   = .false.
    dismax = 0.0
    if (deph>0.05 .and. hrm>=0.01 .and. tp>0.0) then
       if (choice) then
          !
          ! Adjust wave-coefficients
          !
          call wavenr_htk(deph, tp, wavek)
          !
          hmax   = 0.88/wavek*tanh(gamma*wavek*deph/0.88)
          hrm    = min(hrm, hmax)
          hs     = hrm*wort2
          tpmin  = 14.09*sqrt(hs/grav)
          dismax = 0.25*rho*grav*hmax*hmax/tp
          dish   = min(dish, dismax)
          tp     = min(tp, tpmin)
       endif
       !
       call wavenr_htk(deph, tp, wavek)
       !
       wavel = twopi/wavek
    else
       !
       ! Too shallow water or waves too small
       !
       ldep = .true.
    endif
end subroutine corrht

max_stars_repo_name: liujiamingustc/phd

max_issues_repo_name: liujiamingustc/phd

max_forks_repo_name: liujiamingustc/phd


text:
cool = 4*cgs_stef*rho*(4*cgs_k_over_mec2*trad**4*(temp*merge(4* &
      cgs_k_over_mec2*temp + 1, 1.0d0, use_relcompt) - trad)* &
      ksctv(1) + (temp**4 - trad**4)*kabpv(1))
if (use_relcompt) then
   cool_drho = 4*cgs_stef*(4*cgs_k_over_mec2*trad**4*(temp*(4* &
      cgs_k_over_mec2*temp + 1) - trad)*ksctv(1) + rho*(4* &
      cgs_k_over_mec2*trad**4*(temp*(4*cgs_k_over_mec2*temp + 1) - trad &
      )*ksctv(2) + (temp**4 - trad**4)*kabpv(2)) + (temp**4 - trad**4)* &
      kabpv(1))
else
   cool_drho = 4*cgs_stef*(4*cgs_k_over_mec2*trad**4*(1.0d0*temp - trad) &
      *ksctv(1) + rho*(4*cgs_k_over_mec2*trad**4*(1.0d0*temp - trad)* &
      ksctv(2) + (temp**4 - trad**4)*kabpv(2)) + (temp**4 - trad**4)* &
      kabpv(1))
end if
if (use_relcompt) then
   cool_dtemp = 4*cgs_stef*rho*(4*cgs_k_over_mec2*trad**4*(8* &
      cgs_k_over_mec2*temp + 1)*ksctv(1) + 4*cgs_k_over_mec2*trad**4*( &
      temp*(4*cgs_k_over_mec2*temp + 1) - trad)*ksctv(3) + 4*temp**3* &
      kabpv(1) + (temp**4 - trad**4)*kabpv(3))
else
   cool_dtemp = 4*cgs_stef*rho*(4*cgs_k_over_mec2*trad**4*(1.0d0*temp - &
      trad)*ksctv(3) + 4.0d0*cgs_k_over_mec2*trad**4*ksctv(1) + 4*temp &
      **3*kabpv(1) + (temp**4 - trad**4)*kabpv(3))
end if

max_stars_repo_name: gronki/diskvert

max_issues_repo_name: gronki/diskvert

max_forks_repo_name: gronki/diskvert


text:
module volume_averages

   public :: init_volume_averages, finish_volume_averages
   public :: fieldline_average
   public :: volume_average
   public :: flux_surface_average_ffs

   public :: mode_fac

   private

   interface fieldline_average
      module procedure fieldline_average_real
      module procedure fieldline_average_complex
   end interface

   real, dimension(:), allocatable :: mode_fac
   !> Fourier coefficients in y of the Jacobian;
   !> needed for full flux surface simulations
   complex, dimension(:, :), allocatable :: jacobian_ky

contains

   subroutine init_volume_averages

      use zgrid, only: nzgrid, nztot, delzed
      use kt_grids, only: nalpha, aky, nakx, naky, rho_d_clamped
      use stella_geometry, only: geo_surf, drhodpsi
      use stella_geometry, only: geo_surf, jacob, djacdrho, q_as_x, dVolume
      use physics_flags, only: full_flux_surface, radial_variation

      implicit none

      real :: dqdrho

      if (.not. allocated(mode_fac)) then
         allocate (mode_fac(naky)); mode_fac = 2.0
         if (aky(1) < epsilon(0.)) mode_fac(1) = 1.0
      end if
      if (full_flux_surface) then
         call init_flux_surface_average_ffs
      end if

      dqdrho = geo_surf%shat * geo_surf%qinp / geo_surf%rhoc
      if (.not. allocated(dVolume)) allocate (dVolume(nalpha, nakx, -nzgrid:nzgrid))

      !dVolume contains the volume element jacob, which may vary with x or alpha
      ! NB: dVolume does not contain the factor dx, as this should always be uniform
      dVolume = spread(jacob * spread(delzed, 1, nalpha), 2, nakx)
      if (q_as_x) then
         dVolume = dVolume / (dqdrho * drhodpsi)
      end if

      if (radial_variation) then
         if (q_as_x) then
            dVolume = dVolume * (1.+spread(spread(rho_d_clamped, 1, nalpha), 3, nztot) &
                                 * (spread(djacdrho / jacob, 2, nakx) - geo_surf%d2qdr2 / dqdrho &
                                    + geo_surf%d2psidr2 * drhodpsi))
         else
            dVolume = dVolume * (1.+spread(spread(rho_d_clamped, 1, nalpha), 3, nztot) &
                                 * spread(djacdrho / jacob, 2, nakx))
         end if
      end if

      if (full_flux_surface) then
         !something should go here
      end if

      !avoid the double counting at the zed boundaries
      dVolume(:, :, -nzgrid) = 0.5 * dVolume(:, :, -nzgrid)
      dVolume(:, :, nzgrid) = 0.5 * dVolume(:, :, nzgrid)

   end subroutine init_volume_averages

   subroutine finish_volume_averages

      use stella_geometry, only: dVolume
      use physics_flags, only: full_flux_surface

      implicit none

      if (allocated(mode_fac)) deallocate (mode_fac)
      if (allocated(dVolume)) deallocate (dVolume)
      if (full_flux_surface) then
         if (allocated(jacobian_ky)) deallocate (jacobian_ky)
      end if

   end subroutine finish_volume_averages

   !==============================================
   !============ FIELD LINE AVERAGE ==============
   !==============================================
   subroutine fieldline_average_real(unavg, avg)

      use zgrid, only: nzgrid, ntubes
      use kt_grids, only: nakx, naky
      use stella_geometry, only: dl_over_b

      implicit none

      real, dimension(:, :, -nzgrid:, :), intent(in) :: unavg
      real, dimension(:, :), intent(out) :: avg

      integer :: it, ia

      ia = 1

      avg = 0.0
      do it = 1, ntubes
         avg = avg + sum(spread(spread(dl_over_b(ia, :), 1, naky), 2, nakx) * unavg(:, :, :, it), dim=3)
      end do
      avg = avg / real(ntubes)

   end subroutine fieldline_average_real

   subroutine fieldline_average_complex(unavg, avg)

      use zgrid, only: nzgrid, ntubes
      use kt_grids, only: nakx, naky
      use stella_geometry, only: dl_over_b

      implicit none

      complex, dimension(:, :, -nzgrid:, :), intent(in) :: unavg
      complex, dimension(:, :), intent(out) :: avg

      integer :: it, ia

      ia = 1

      avg = 0.0
      do it = 1, ntubes
         avg = avg + sum(spread(spread(dl_over_b(ia, :), 1, naky), 2, nakx) * unavg(:, :, :, it), dim=3)
      end do
      avg = avg / real(ntubes)

   end subroutine fieldline_average_complex

   !==============================================
   !============== VOLUME AVERAGE ================
   !==============================================
   subroutine volume_average(unavg, avg)

      use zgrid, only: nzgrid, ntubes
      use kt_grids, only: naky, nakx
      use stella_geometry, only: dl_over_b

      implicit none

      complex, dimension(:, :, -nzgrid:, :), intent(in) :: unavg
      real, intent(out) :: avg

      integer :: iky, ikx, iz, it, ia

      ia = 1

      avg = 0.
      do it = 1, ntubes
         do iz = -nzgrid, nzgrid
            do ikx = 1, nakx
               do iky = 1, naky
                  avg = avg + real(unavg(iky, ikx, iz, it) * conjg(unavg(iky, ikx, iz, it))) * mode_fac(iky) * dl_over_b(ia, iz)
               end do
            end do
         end do
      end do
      avg = avg / real(ntubes)

   end subroutine volume_average

   subroutine init_flux_surface_average_ffs

      use zgrid, only: nzgrid
      use kt_grids, only: naky
      use stella_geometry, only: jacob
      use stella_transforms, only: transform_alpha2kalpha

      implicit none

      integer :: iz

      if (.not. allocated(jacobian_ky)) allocate (jacobian_ky(naky, -nzgrid:nzgrid))

      !> calculate the Fourier coefficients in y of the Jacobian
      !> this is needed in the computation of the flux surface average of phi
      do iz = -nzgrid, nzgrid
         call transform_alpha2kalpha(jacob(:, iz), jacobian_ky(:, iz))
      end do

   end subroutine init_flux_surface_average_ffs

   subroutine flux_surface_average_ffs(no_fsa, fsa)

      use zgrid, only: nzgrid, delzed
      use stella_geometry, only: jacob
      use kt_grids, only: naky, naky_all, nalpha
      use kt_grids, only: dy

      implicit none

      complex, dimension(:, -nzgrid:), intent(in) :: no_fsa
      complex, intent(out) :: fsa

      integer :: iky, ikymod
      real :: area

      ! the normalising factor int dy dz Jacobian
      area = sum(spread(delzed * dy, 1, nalpha) * jacob)

      fsa = 0.0
      ! get contribution from negative ky values
      ! for no_fsa, iky=1 corresponds to -kymax, and iky=naky-1 to -dky
      do iky = 1, naky - 1
         ! jacobian_ky only defined for positive ky values
         ! use reality of the jacobian to fill in negative ky values
         ! i.e., jacobian_ky(-ky) = conjg(jacobian_ky(ky))
         ! ikymod runs from naky down to 2, which corresponds
         ! to ky values in jacobian_ky from kymax down to dky
         ikymod = naky - iky + 1
         ! for each ky, add the integral over zed
         fsa = fsa + sum(delzed * no_fsa(iky, :) * jacobian_ky(ikymod, :))
      end do
      ! get contribution from zero and positive ky values
      ! iky = naky correspond to ky=0 for no_fsa and iky=naky_all to ky=kymax
      do iky = naky, naky_all
         ! ikymod runs from 1 to naky
         ! ikymod = 1 corresponds to ky=0 for jacobian_ky and ikymod=naky to ky=kymax
         ikymod = iky - naky + 1
         ! for each ky, add the integral over zed
         fsa = fsa + sum(delzed * no_fsa(iky, :) * conjg(jacobian_ky(ikymod, :)))
      end do
      ! normalise by the flux surface area
      fsa = fsa / area

   end subroutine flux_surface_average_ffs

end module volume_averages

max_stars_repo_name: AntonioG-Jerez/stella

max_issues_repo_name: AntonioG-Jerez/stella

max_forks_repo_name: AntonioG-Jerez/stella


text:
      SUBROUTINE INCLIN(X,V,XCM,VCM,ALPHA)
*
*
*       Inclination of hierarchical system.
*       -----------------------------------
*
      REAL*8  A12,A22,A1A2,FAC,ALPHA
      REAL*8  XCM(3),VCM(3),X(3,3),V(3,3),A0(3),A2(3)
*
*
*       Define inner binary I1 & I2 and outer component I3.
      I1 = 1
      I2 = 2
      I3 = 3
*
*       Evaluate scalar products of angular momenta.
      A12 = 0.0
      A22 = 0.0
      A1A2 = 0.0
      DO 20 K = 1,3
          K1 = K + 1
          IF (K1.GT.3) K1 = 1
          K2 = K1 + 1
          IF (K2.GT.3) K2 = 1
          A0(K) = (X(K1,I1) - X(K1,I2))*(V(K2,I1) - V(K2,I2))
     &          - (X(K2,I1) - X(K2,I2))*(V(K1,I1) - V(K1,I2))
          A2(K) = (X(K1,I3) - XCM(K1))*(V(K2,I3) - VCM(K2))
     &          - (X(K2,I3) - XCM(K2))*(V(K1,I3) - VCM(K1))
          A12 = A12 + A0(K)**2
          A22 = A22 + A2(K)**2
          A1A2 = A1A2 + A0(K)*A2(K)
   20 CONTINUE
*
*       Determine inclination in radians.
      FAC = A1A2/SQRT(A12*A22)
      FAC = MIN(FAC,1.0D0)
      ALPHA = ACOS(FAC)
*
      RETURN
*
      END


text:
module mod1
#ifdef _OPENMP
integer, parameter :: bla = 10
#endif
end module mod1

max_stars_repo_name: clementval/claw-compiler

max_issues_repo_name: skosukhin/claw-compiler

max_forks_repo_name: skosukhin/claw-compiler


text:
    subroutine ss
     type nums
      integer scalarStart
      integer array1(2)
      integer array2(2)
      integer scalarEnd
    end type nums

     type nested_nums
      integer i(2)
      integer rr(2)
      type (nums) j
    end type nested_nums

!    type (nums),parameter::con21=nums((/1,2/),(/1.2,2.3/))
!    type (nums),parameter::con22=nums((/1,2/))
!    type (nums),parameter::con22=nums((/1,2/),(/3,4/))
     type (nums),parameter::con22=nums(0,(/1,2/),(/3,4/),5)
!    type (nums),parameter::con22
!    type (nums),parameter::con22=nums((/1,2/))

!    type (nested_nums),parameter::nested=nested_nums((/1,2/),(/3,4/),nums((/5,6/),(/7,8/)))

     integer,parameter::testParameter=1
!    integer,dimension(2),parameter::testParameterArray=(/7,8/)
    end

max_stars_repo_name: maurizioabba/rose

max_issues_repo_name: sujankh/rose-matlab

max_forks_repo_name: sujankh/rose-matlab


text:
!>
!! @file
!! This file is part of SeisSol.
!!
!! @author Josep de la Puente Alvarez (josep.delapuente AT bsc.es, http://www.geophysik.uni-muenchen.de/Members/jdelapuente)
!!
!! @section LICENSE
!! Copyright (c) 2008, SeisSol Group
!! All rights reserved.
!!
!! Redistribution and use in source and binary forms, with or without
!! modification, are permitted provided that the following conditions are met:
!!
!! 1. Redistributions of source code must retain the above copyright notice,
!!    this list of conditions and the following disclaimer.
!!
!! 2. Redistributions in binary form must reproduce the above copyright notice,
!!    this list of conditions and the following disclaimer in the documentation
!!    and/or other materials provided with the distribution.
!!
!! 3. Neither the name of the copyright holder nor the names of its
!!    contributors may be used to endorse or promote products derived from this
!!    software without specific prior written permission.
!!
!! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
!! AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
!! IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
!! ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
!! LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
!! CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
!! SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
!! INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
!! CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
!! ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
!! POSSIBILITY OF SUCH DAMAGE.

#ifdef BG
#include "../Initializer/preProcessorMacros.fpp"
#else
#include "Initializer/preProcessorMacros.fpp"
#endif

MODULE calc_deltaT_mod
  !--------------------------------------------------------------------------
  USE TypesDef
  !--------------------------------------------------------------------------
  IMPLICIT NONE
  PRIVATE
  !----------------------------------------------------------------------------
  INTERFACE ini_calc_deltaT
     MODULE PROCEDURE ini_calc_deltaT
  END INTERFACE

  INTERFACE close_calc_deltaT
     MODULE PROCEDURE close_calc_deltaT
  END INTERFACE

  INTERFACE cfl_step
     MODULE PROCEDURE cfl_step
  END INTERFACE
  !----------------------------------------------------------------------------
  PUBLIC  :: ini_calc_deltaT
  PUBLIC  :: close_calc_deltaT
  public  :: cfl_step
  !----------------------------------------------------------------------------
  ! Module Variables
  REAL             :: tol

CONTAINS

  SUBROUTINE ini_calc_deltaT(EqType,OptionalFields,EQN,MESH,IO)
    !--------------------------------------------------------------------------
    IMPLICIT NONE
    !--------------------------------------------------------------------------
    TYPE (tEquations)             :: EQN                                         !
    TYPE (tUnstructMesh)          :: MESH
    TYPE (tUnstructOptionalFields):: OptionalFields
    TYPE (tInputOutput)           :: IO
    INTEGER                       :: EqType
    ! local Variables
    INTEGER                       :: allocStat
    !--------------------------------------------------------------------------
    INTENT(IN)                    :: MESH,EqType,IO
    INTENT(INOUT)                 :: OptionalFields
    !--------------------------------------------------------------------------
    !                                                                            !
    tol = 1.0 / (10.0**(PRECISION(1.0)-2) )                                      !
    !                                                                            !                                                    !
    ALLOCATE(                                             &                      ! Allocate
               OptionalFields%vel(         MESH%nElem)    , &                    ! Allocate
               OptionalFields%sound(       MESH%nElem)    , &                    ! Allocate
               OptionalFields%mask(        MESH%nElem)    , &                    ! Allocate
               OptionalFields%dt_convectiv(MESH%nElem)    , &                    ! Allocate
               OptionalFields%dtmin       (MESH%nElem)    , &                    ! Allocate
               STAT = allocStat                             )                    ! Allocate
    !                                                                            !
    IF (allocstat .NE. 0 ) THEN                                                  ! Error Handler
         logError(*) 'ALLOCATE ERROR in ini_calc_deltaT!'            ! Error Handler
         STOP                                                                    ! Error Handler
    END IF                                                                       ! Error Handler
    !                                                   !
    !                                                   !
  END SUBROUTINE ini_calc_deltaT                        !

  SUBROUTINE close_calc_deltaT(OptionalFields)
    !--------------------------------------------------------------------------
    IMPLICIT NONE
    !--------------------------------------------------------------------------
    TYPE (tUnstructOptionalFields):: OptionalFields
    !--------------------------------------------------------------------------
    !                                                   !
    IF (allocated(OptionalFields%vel)) THEN            !
        DEALLOCATE(OptionalFields%vel)                  ! Deallocate
    END IF                                              !
    IF (allocated(OptionalFields%sound)) THEN          !
        DEALLOCATE(OptionalFields%sound)                ! Deallocate
    END IF                                              !
    IF (allocated(OptionalFields%dt_convectiv)) THEN   !
        DEALLOCATE(OptionalFields%dt_convectiv)         ! Deallocate
    END IF                                              !
    IF (allocated(OptionalFields%mask)) THEN           !
        DEALLOCATE(OptionalFields%mask)                 ! Deallocate
    END IF
    if (allocated(OptionalFields%dtmin)) then
        deallocate(OptionalFields%dtmin)
    end if
! aheineck, @TODO, not referecned in the code, commented                                              !
!    IF (ASSOCIATED(OptionalFields%dt_viscos)) THEN      !
!        DEALLOCATE(OptionalFields%dt_viscos)            ! Deallocate
!    END IF                                              !
  END SUBROUTINE close_calc_deltaT                      !

  SUBROUTINE cfl_step(OptionalFields,EQN,MESH,DISC,IO)
    !--------------------------------------------------------------------------
    IMPLICIT NONE
    !--------------------------------------------------------------------------
    ! argument list declaration
    TYPE (tUnstructOptionalFields)        :: OptionalFields
    TYPE (tEquations)                     :: EQN
    TYPE (tUnstructMesh)                  :: MESH
    TYPE (tDiscretization)                :: DISC
    TYPE (tInputOutput)                   :: IO
    ! local variable declaration
    REAL,POINTER                          :: MaterialVal(:,:)
    INTEGER                               :: iElem, iNeighbor, iSide
    INTEGER                               :: idxNeighbors(MESH%GlobalElemType)
    REAL                                  :: rho, C(6,6)
    !--------------------------------------------------------------------------
    INTENT(IN)                            :: EQN, MESH, IO
    INTENT(INOUT)                         :: OptionalFields
    INTENT(INOUT)                         :: DISC
    !--------------------------------------------------------------------------
    !                                                                         !
    ! Compute Velocities
    !
    !
       DO iElem = 1, MESH%nElem
           OptionalFields%sound(iElem) = MAXVAL( DISC%Galerkin%MaxWaveSpeed(iElem,:) )
       ENDDO

       SELECT CASE(EQN%Advection)
       CASE(0)
           OptionalFields%vel(:) = 0.
       CASE(1)
           IF ((.NOT.ASSOCIATED(OptionalFields%BackgroundValue))) THEN
              logError(*) 'OptionalFields%BackgroundValue not associated!'
              logError(*) 'although linearized Equations are specified. '
              STOP
           END IF
           MaterialVal => OptionalFields%BackgroundValue
           OptionalFields%vel(:) = SQRT( MaterialVal(:,4)**2 + MaterialVal(:,5)**2 + MaterialVal(:,6)**2 )
       END SELECT
    !
    ! Berechne dt aus dt_convectiv und dt_fix
    !             fuer jedes Element
    !
    !
    ! U N S T E A D Y   S O L U T I O N
    !
    OptionalFields%dt_convectiv(:)  = DISC%CFL*2.0*MESH%ELEM%MinDistBarySide(:)  &
                                     / ( OptionalFields%sound(:) + OptionalFields%vel(:) )
   !
   IF(DISC%Galerkin%pAdaptivity.GE.1) THEN
    OptionalFields%dt_convectiv(:) = OptionalFields%dt_convectiv(:) / &
                                     ( 2.*DISC%Galerkin%LocPoly(:)+1  )
   ELSE
    OptionalFields%dt_convectiv(:) = OptionalFields%dt_convectiv(:) / &
                                     ( 2.*DISC%Galerkin%nPoly+1       )
   ENDIF
    !

    ! set to min(dt_convectiv,fixtimestep) for each element
    WHERE(OptionalFields%dt_convectiv(:).GT.DISC%FixTimeStep)
        OptionalFields%dt_convectiv(:) = DISC%FixTimeStep
    END WHERE

    !
    IF(DISC%Galerkin%DGMethod.EQ.3) THEN
      DISC%LocalDt(:) = OptionalFields%dt_convectiv(:)
      DO iElem = 1, MESH%nElem
        OptionalFields%dtmin(iElem) = DISC%LocalDt(iElem)
        DO iSide = 1, MESH%LocalElemType(iElem)
            iNeighbor = MESH%ELEM%SideNeighbor(iSide,iElem)
            IF(iNeighbor.LE.MESH%nElem) THEN
                OptionalFields%dtmin(iElem) = MIN(OptionalFields%dtmin(iElem), DISC%LocalDt(iNeighbor))
            ENDIF
        ENDDO
      ENDDO
      DO iElem = 1, MESH%nElem
        DISC%LocalDt(iElem) = OptionalFields%dtmin(iElem)
      ENDDO
      !
      DO iElem = 1, MESH%nElem
        ! Match printtime if necessary
        IF(DISC%LocalTime(iElem)+DISC%LocalDt(iElem).GE.DISC%printtime) THEN
            DISC%LocalDt(iElem) = DISC%printtime - DISC%LocalTime(iElem)
        ENDIF
        ! Match endtime if necessary
        IF(DISC%LocalTime(iElem)+DISC%LocalDt(iElem).GE.DISC%EndTime) THEN
            DISC%LocalDt(iElem) = DISC%EndTime - DISC%LocalTime(iElem)
        ENDIF
      ENDDO
    ENDIF
    !
    !
    RETURN
    !
  END SUBROUTINE cfl_step

END MODULE calc_deltaT_mod

max_stars_repo_name: ivotron/SeisSol

max_issues_repo_name: ivotron/SeisSol

max_forks_repo_name: ivotron/SeisSol