c Subroutines related to time including setting the time step
c setting boundary conditions, crop variables and fertilization
c as needed by the various schedules.
c
c Determine the next time step, dt, and update the time,t
      subroutine setTime(NTFix,ipFix,t,tOld,tFix,dt,dtMin,dtMax,
     ;                   dtOld,Iter,tOutput,nTOut,iPrint,itDay,
     ;                   lPrint,lFixFlag,lPrnFlag,dayFlag,ierr,
     ;                   errStr)
      implicit none
      double precision t,tFix,dt,dtMin,dtMax,tOld,dtOld,
     ;                 tOutput
      integer ipFix,NTFix,Iter,iPrint,nTout,itDay,ierr
      logical lFixFlag,lPrnFlag,lPrint,dayFlag
      dimension tFix(NTFix),tOutput(nTOut)
      character*80 errStr
c
c     local variables
      double precision zero,one,err,errSmall,x,y,dtMul,dtMul1,dtMul2
      data zero,one,dtMul1,dtMul2 /0.0D00,1.0D00,0.5D00,1.3D00/
      err(x,y)=dabs(1.D00-x/y)-1.0D-03
      errSmall(x,y)=dabs(1.D00-x/y)-1.0D-05
      if (iPrint.lt.nTOut.and.(tOutput(iPrint).gt.tOutput(iPrint+1)))
     ;then
        errStr='Time.in'
        ierr=90
        return
      endif
c
c     Set the dt multiplier
      if (Iter.gt.7) then
        dtMul=dtMul1
      elseif (Iter.le.3) then
        dtMul=dtMul2
      else
        dtMul=1.0D00
      endif
c
c     When time approaches a printout set the trigger flag
      if (err(t,tOutput(iPrint)).lt.zero) then
        lPrnFlag=.true.
      endif
      dtOld=dt
c Determine whether adding dt from watflow to the current time
c will be greater than t

      if ((err(t+dt,tFix(ipFix)).lt.zero).or.
     ;   (t+dt.gt.tFix(ipFix))) then
c if the projected time and the next fixed time are almost identical
c or the projected time is greater than the next fixed time
c then set dt to be one half the current delta t to the next fixed time. 
c This will decrease time steps to dtMin in advance of a change in boundary
c conditions.
        lFixFlag=.true.
        if ((tFix(ipFix)-t)/1.3D00.gt.dtMin) then
          dt=(tFix(ipFix)-t)/1.3D00
        else
          dt=dtMin
        endif
      else !Use the dtMul factor
        dt=min(dtMul*dt,dtMax)
      endif
      dt=min(dt,dtMax)
c Compute next time

      if (t.lt.tFix(ipFix).and.tFix(ipFix)-t.lt.dt) then
        dt=tFix(ipFix)-t
      endif
      tOld=t
      t=t+dt
c
c Test for the new time exceeding the fixed time and the flag indicating
c that the increment of ipFix has not yet occurred.

      if (t.gt.tFix(ipFix).and.lFixFlag) then
        ipFix=ipFix+1
        lFixFlag=.false.
      endif
c
c As soon as time is equal to or exceeds the printout time and the
c trigger flag is set then set the print flag
      if (t.ge.tOutput(iPrint).and.lPrnFlag) then
        lPrint=.TRUE.
        lPrnFlag=.FALSE.
      endif
c
c Daily values are printed immediately after midnight using dayFlag
      if (t.gt.dble(itday)+one) then
        itday=int(t)
        dayFlag=.TRUE.
      endif
      return
      end
c
c SetBC transfers values from boundary condition 2-D array to variables
      subroutine setBC(nBounds,iBound,nSols,bndCond,irriChem,t,tFinal,
     ;                 dtMax,bndTime,nextTime,rTop,tPotMax,tpotN,
     ;                 tPotNMax,KodTop,KodBot,kTOld,kBOld,chmBCTop,
     ;                 chmBCBot,rBot,cTop,cBot,cropNum,cStarN,hBot,
     ;                 wDayEvap,lNitrate,freeD,initFlag,ierr)
      implicit none
      include "defs.fi"
      logical initFlag,lNitrate,freeD
      integer iBound,nBounds,nSols,KodTop,KodBot,kTOld,kBOld,cropNum,
     ;        chmBCTop,chmBCBot,ierr
      double precision bndCond,irriChem,bndTime,rTop,t,tFinal,
     ;       tPotMax,tpotN,nextTime,rBot,cTop,cBot,dtMax,cStarN,hBot,
     ;       wDayEvap,tPotNMax
      dimension bndCond(nBndVars,*),bndTime(nBounds),
     ;          irriChem(nIrriC*2,nSols),cTop(nIrriC),
     ;          cBot(nIrriC)
c
      integer iSol,i
      double precision zero,one,rSoil,PotETo,bndIrr,Kc,NUptake,x,y,err
      data zero,one /0.0D00,1.0D00/
      err(x,y)=dabs(one-x/y)-1.0D-02
c
c     Test for out-of-sequence boundary condition
      if (iBound.lt.nBounds.and.(bndTime(iBound).ge.bndTime(iBound+1)))
     ;then
        ierr=91
        return
      endif

      if (t.lt.bndTime(1)) then
        if (bndTime(1)-t.lt.dtMax) then
c         time to first BC change is very short so skip the first one
          nextTime=bndTime(2)
        else
          nextTime=bndTime(1)
        endif
c     return the time of the next scheduled BC change
      elseif (iBound.lt.nBounds) then
        nextTime=bndTime(iBound+1)
      else
        nextTime=tFinal
      endif

      if (initFlag) then
        wDayEvap=zero  !first day rTop is zero
        rBot=zero
        kodTop = -4
        if (freeD) then
          kodBot = -1  !free drainage
          chmBCBot=0
        else
          kodBot = 1
          chmBCBot=1  !Dirichlet condition
          !start hBot at the first hBot in the BC file to avoid
          !sudden rapid flow due to jump discontinuity
          hBot=bndCond(9,1)
        endif
        kTOld = kodTop
        kBOld = kodBot
        rTop = 0
        chmBCTop=-1
        do i=1,nIrriC
          cTop(i)=0
          cBot(i)=0
        enddo
        cStarN=zero !Initialize critical nitrogen concentration
      endif
c
c       If the user has specified the first record for a time
c       later than the start of the simulation then initialize
c       the boundary conditions to zeros.
      if (t.lt.bndTime(1)) then
        wDayEvap=zero
        TPotMax=zero
        NUptake=zero
        chmBCTop=-1
        chmBCBot=0
        tPotNMax=zero
      else
        PotETo=bndCond(1,iBound)   !Potential ET (cm/day)
        bndIrr=bndCond(2,iBound)   !Irrigation rate (cm/day)
        Kc=bndCond(3,iBound)       !Crop Coefficient
        NUptake=bndCond(4,iBound)  !NUptake
        chmBCTop=int(bndCond(5,iBound))
        chmBCBot=int(bndCond(6,iBound))
        iSol=int(bndCond(7,iBound))
        if (iSol.gt.0) then
          if (iSol.gt.nSols) then
            ierr=80
            return
          endif
        endif
        if (lNitrate) then
c         Set user-supplied cStarN which must be zero if the N uptake
c         routine is using TPotN/TPot as the cStar value
          cStarN=bndCond(8,iBound)
        else
          cStarN=zero
        endif
c
c       When free drainage is false the pressure head at the bottom node
c       is specified as hBot
        if (KodBot.ge.1) then
          hBot=bndCond(9,iBound)
        endif
c       Calculate the soil surface boundary condition based on existence
c       of a crop or bare soil surface in case of evaporation. Normally
c       the user should specify Kc=1.0 for soil surface evaporation meaning
c       that the actual evaporation will be equal to potential.
        wDayEvap=zero
        if (cropNum.eq.0) then
          rTop=Kc*PotETo !fallow condition, soil surface evaporation
          if (bndIrr.gt.1.0D-02) then
            rTop=rTop-bndIrr !subtract finite irrigation rate
            wDayEvap=Kc*PotETo !soil surface evap when water was applied
          endif
          tpotN=zero
          tpotMax=zero
          tPotNMax=zero
        else
c         All the water flux to the atmosphere is transpiration
          rTop=min(-bndIrr,zero)
          TPotMax=Kc*PotETo
          tpotN=NUptake
          tPotNMax=NUptake
        endif
c
c       Set the concentration vectors for upper and lower boundaries
        if (iSol.gt.0) then
          do i=1,nIrriC
            cTop(i)=irriChem(i,iSol)
            cBot(i)=irriChem(i+3,iSol)
          enddo
        else !assume all concentrations are zero
          do i=1,nIrriC
            cTop(i)=0
            cBot(i)=0
          enddo
        endif
      endif
      return
      end
c
c Crop Schedule variables:
c cropSchd(1,iRecs) = cropNum
c cropSchd(2,iRecs) = tRMin
c cropSchd(3,iRecs) = tRHarv
c cropSchd(4,iRecs) = RDFDay
c cropSchd(5,iRecs) = xRMin
      subroutine advCrop(cropSchd,cropNum,iCrop,scdNam,cropName,
     ;                   nPlants,cropProp,crpIndex,t,relY,sumRoot,
     ;                   sumDel,sPotN,sLeachN,prodNSea,Et_Shres,
     ;                   initFlag,nCrpRecs,cropDone,lNitrate,
     ;                   prevCrop,sEPotMax,ierr,errStr)
      implicit none
      include "defs.fi"
      integer cropNum,iCrop,nCrpRecs,crpIndex,nPlants,prevCrop
      logical initFlag,cropDone,lNitrate
      double precision cropSchd,t,cropProp,relY,sumRoot,sumDel,
     ;                 sPotN,sLeachN,Et_shres,prodNSea,sEPotMax
      dimension cropSchd(nSchdPar,*),cropProp(nCropPar,*)
      character*20 scdNam(maxScd),cropName
      character*80 errStr
      integer ierr
c     local variables
      double precision hShres,sShres,gg5,hh5,one
      logical lChange
c
      integer getCIndx
      external getCIndx
      data one /1.0D00/
      ierr = 0
      cropDone=.FALSE.
      lChange=.FALSE.
c
c       There must be either a crop or fallow if lPlant is true
c       so set the first crop if initFlag is set
      if (initFlag) then
        cropNum=INT(cropSchd(1,iCrop)) !Crop identification number
        cropName=scdNam(iCrop)
        initFlag=.FALSE.
        lChange=.TRUE.
      elseif (iCrop.ge.nCrpRecs .and. t.ge.cropSchd(3,iCrop)) then
c       Dont write final crop records for bare soil or fallow
        if (cropNum.gt.1) then
          call wrtCrop(.FALSE.,t,cropNum,cropName,relY,sumRoot,
     ;                 sPotN,sLeachN,prodNSea,sEPotMax,
     ;                 sumDel,lNitrate,ierr)
        endif
        cropDone=.TRUE.
        prevCrop=cropNum
c       the last crop is finished
        return
      elseif (iCrop.lt.nCrpRecs .and. t.ge.cropSchd(2,iCrop+1)) then
c       write out data for the end of the previous crop and advance to
c       the next one
        if (cropNum.gt.0) then
          call wrtCrop(.FALSE.,t,cropNum,cropName,relY,sumRoot,
     ;                 sPotN,sLeachN,prodNSea,sEPotMax,sumDel,
     ;                 lNitrate,ierr)
        endif
        prevCrop=cropNum
        iCrop=iCrop+1
        cropNum=INT(cropSchd(1,iCrop))
        cropName=scdNam(iCrop)
        lChange=.TRUE.
      endif
c     Find the correct index for this crop in the cropProp 
c     array and use the index to obtain the stress parameter
c     values for calculating the threshold fractional stress
c     that causes water uptake reduction
      if (lChange) then
        crpIndex=getCIndx(nPlants,cropNum,cropProp,ierr,errStr)
        if (ierr.gt.0) call errHand(nErrLog,ierr,errStr)
c
c       Calculate the threshold for sink term.
        hShres=cropProp(3,crpIndex) !(cm)
        sShres=cropProp(5,crpIndex) !units for osmotic head (cm)
        gg5=cropProp(6,crpIndex)    !(cm)
        hh5=cropProp(4,crpIndex)    !(cm)
        Et_shres = one/(one + (hShres/hh5 + sShres/gg5)**3)
      endif
      return
      end
c
c     This subroutine handles a change in crop and returns the
c     correct rTop,tPotMax and TPotN given the existing boundary
c     conditions.
      subroutine bcAccom(bndCond,iBound,cropProp,cropNum,rTop,
     ;                   tPotMax,tPotN,tPotNMax,tPot,nPlants,
     ;                   ierr,errStr)
      implicit none
      include "defs.fi"
      double precision bndCond(nBndVars,*),rTop,tPotMax,TPotN,tPot,
     ;                 cropProp(nCropPar,*),tPotNMax
      integer iBound,cropNum,nPlants,ierr
      character*80 errStr
c
      integer crpIndex,getCIndx
      double precision Nuptake,PotETo,Kc,Irrig,zero,one
      data zero,one / 0.0D00,1.0D00 /
c     Accomodate the existing boundary conditions, if there is
c     irrigation then bndCond(2,iBound) < 0
      PotETo=bndCond(1,iBound)
      Kc=bndCond(3,iBound)
      Irrig=bndCond(2,iBound)
      NUptake=bndCond(4,iBound)  !NUptake
      if (Irrig.gt.zero) then
        if (cropNum.eq.0) then
c         Include fallow soil surface potential evap in rTop
          rTop = Kc*PotETo - Irrig
          TPotN = zero
          tPotMax = zero
          tPot = zero
        else
c         Don't change the irrigation regime
c         but set the irrigation B.C. and potential transpiration
          rTop = -Irrig
          tPotMax=Kc*PotETo
c         Initialize tPot to tPotMax at the start of a crop
          tPot = tPotMax
          TPotN = Nuptake
          TPotNMax = Nuptake
        endif
      elseif (cropNum.eq.0) then
c       For fallow use soil surface evap as rTop
        rTop = Kc*PotETo
        tPotMax = zero
        tPot = zero
        TPotN = zero
      elseif (cropNum.gt.0) then
c       If there is a crop then there is no soil surface 
c       evaporation and all of Kc*ETo is potential transpiration
        rTop = zero
        TPotMax = Kc*PotETo
        tPot= TPotMax
        TPotN = NUptake
        TPotNMax = Nuptake
      else
        ierr = 1
      endif
	  if (cropNum.gt.0) then
        crpIndex=getCIndx(nPlants,cropNum,cropProp,ierr,errStr)
        if (ierr.gt.0) call errHand(nErrLog,ierr,errStr)
c       Set rel_fact at the start of the crop to one
        cropProp(7,crpIndex)=one
      endif
      return
      end
c
c     If mineral N was added then dissolve it into Conc array, if 
c     organic N was added then it updates the cFix array
      subroutine addFert(N,NLayers,nDiscret,NMat,z,iApplic,t,tApplic,
     ;                   tNextEvt,dtMax,amtAppl,minFrac,orgFrac,CFix,
     ;                   layerMix,tFinal,layerDep,MatNum,chemPar,
     ;                   Conc,theta,Psi,PsiApp,volConv,mStart,ierr)
      implicit none
      include "defs.fi"
      integer N,NLayers,iApplic,nDiscret,MatNum,NMat,ierr
      logical layerMix(NLayers)
      double precision tApplic,tNextEvt,amtAppl,minFrac,orgFrac,CFix,t,
     ;                 dtMax,layerDep,chemPar,tFinal,z,Conc,theta,
     ;                 Psi,PsiApp,volConv,mStart,tempMass,sumMass
      dimension tApplic(nDiscret),amtAppl(nDiscret),minFrac(nDiscret),
     ;          orgFrac(nDiscret),MatNum(N),cFix(nFixC,N),
     ;          layerDep(NLayers),chemPar(nChemPar,NMat),z(N),
     ;          Conc(nConcD,N),theta(N),Psi(N),PsiApp(nDiscret)
c
      integer i,M,layNum,layOfNod,lastNum
      double precision dz,massOrgN,massMinN,thick,layThick,top,
     ;                 massApp,massNs,massNf,zero,one,two,eps
      dimension massMinN(N)
      data zero,one,two,eps /0.0D00,1.0D00,2.0D00,1.D-04/
c
      if (iApplic.gt.nDiscret) then
        ierr=81
        return
      endif
      if (t.lt.tApplic(1)) then
        if (tApplic(1)-t.lt.dtMax) then
c         time to first fertilization is very short so skip the first one
          tNextEvt=tApplic(2)
        else
          tNextEvt=tApplic(1)
        endif
c     return the time of the next scheduled fertilizer addition
      elseif (iApplic.lt.nDiscret) then
        tNextEvt=tApplic(iApplic+1)
      else
        tNextEvt=tFinal
      endif
c
      top=zero
      dz=z(N)-z(N-1)
      lastNum=layOfNod(N,layerDep,dz,N-1,NLayers,ierr)
      tempMass=zero
      sumMass=zero
      do i=N-1,1,-1
        MassMinN(i)=zero
        dz=z(i+1)-z(i)
        M = MatNum(i)
        layNum = layOfNod(N,layerDep,dz,i,NLayers,ierr)
        thick = abs(layThick(NLayers,layNum,top,layerDep))
        if (layNum.ne.lastNum) then
          top=layerDep(lastNum)
          lastNum=layNum
        endif
c       do the mixing for this node if the layer is being mixed
        if (layerMix(layNum)) then
c           get mass of N already present, g/g soil*volume of element*rhoBulk
            massOrgN = cFix(1,i)*dz*chemPar(2,M)
            massMinN(i) = cFix(2,i)*dz*chemPar(2,M)
c
c           Recalculate Psi(z) based on Psi of the existing orgN in 
c           the element and the orgN added with PsiApp of the applied
c           material.
            if (orgFrac(iApplic).gt.eps) then
              massApp = amtAppl(iApplic)*orgFrac(iApplic)*(dz/thick)
              massNs = massOrgN*Psi(i) + massApp*PsiApp(iApplic)
              massNf = massOrgN*(one - Psi(i)) +
     ;               massApp*(one - PsiApp(iApplic))
              Psi(i) = massNs/(massNs + massNf)
            else
c             Assume that no organic-N was applied. Psi(i) remains 
c             unaltered because no organic-N was applied so there was
c             no update to the Ns and Nf fractions of any existing
c             organic-N in CFix
              massApp=zero
              orgFrac(iApplic)=zero
            endif
            massOrgN = massOrgN + massApp
c           additional mass amtAppl is g/cm^2 to be applied in
c           layer layNum. The total applied to the layer is divided among
c           the elements according to a fractional element thickness by
c           layer thickness
c
c           massOrgN in g in the element, amtAppl in g/cm2 in the top layer
            massMinN(i)=(Conc(3,i)*theta(i))*dz*1.0D-06*1.4D01
            tempMass=tempMass+massMinN(i)*1.0D05
            massMinN(i) = massMinN(i) + amtAppl(iApplic)*
     ;                 (one - orgFrac(iApplic))*(dz/thick)
            sumMass=sumMass+massMinN(i)*1.0D05
c           recompute concentrations mmol
            cFix(1,i) = massOrgN/(dz*chemPar(2,M)) !cFix in g/g soil
c
c           transfer the cFix mineral mass by dissolving it in the
c           soil water. (mmol/L)*cm*(cm^3/cm^3)*(L/cm^3)*(mol/mmol)*(g/mol).
c           Starting massMinN(g) is the N mass in element i (g)
            cFix(2,i)=zero
        endif
      enddo
      mStart=zero
      do i=N-1,1,-1
        dz=z(i+1)-z(i)
        M = MatNum(i)
        layNum = layOfNod(N,layerDep,dz,i,NLayers,ierr)
c       Only recalculate concentrations for the layers that have 
c       fertilizer mixed in. If fertilizer is applied on top of
c       soil it should be dissolved in irrigation water.
        if (layerMix(layNum)) then
            Conc(3,i)=massMinN(i)/(dz*theta(i)*1.4D01*1.0D-06)
        endif
c       mstart is for debugging
        if (layerMix(layNum)) then
          mStart=mStart + dz*1.4D01*1.0D-06*1.0D05*
     ;        (Conc(3,i)*theta(i))
        endif
      enddo
      Conc(3,N)=Conc(3,N-1)
c     Top Node
      cFix(1,N)=cFix(1,N-1)
      cFix(2,N)=cFix(2,N-1)
      return
      end
