c     Root water uptake with compensation
      subroutine Sinker(N,NMat,iRoot,crpIndex,dt,dz,hNew,Beta,se0,
     ;                  theta0,cropProp,tPotMax,MatNum,ParD,vRoot,
     ;                  Et_shres,alphaG,alphaHS,alphaHSR,alfaGHS,
     ;                  alfaEff,Sink)
      implicit none
      include "defs.fi"
      integer N,iRoot,crpIndex,MatNum,NMat
      double precision dz,hNew,Beta,tPotMax,vRoot,alphaG,theta0,
     ;                 se0,alphaHS,dt,ParD,Sink,cropProp,alfaGHS,
     ;                 alfaEff,alphaHSR,Et_shres
      dimension hNew(N),Beta(N),se0(N),theta0(N),alphaHS(N),
     ;          MatNum(N),ParD(nTPar,NMat),Sink(N),alphaHSR(N),
     ;          cropProp(nCropPar,*)
c
      double precision hShres,sShres,gg5,hh5,rel_fact,tPot,
     ;                 h,newSink,vadd,diff,theta,FALfa,fh,
     ;                 se,reduc,zero,one,eps
      dimension h(N),newSink(N),theta(N),se(N),reduc(N)
      logical lReduc,lAllNodes,done
      integer i,iter
      external fh
      data zero,one,eps /0.0D00,1.0D00,1.0D-06/
c
c     Zero the reduction function vectors
      do i=N-1,1,-1
        alphaHS(i)=zero
        alphaHSR(i)=zero
      enddo
c     if there is no rooting depth then return immediately
      if (iRoot.eq.0 .or. tPotMax.eq.zero) then
c       Zero the Sink array and global alphas then return
        do i=N-1,1,-1
          Sink(i)=zero
        enddo
        alfaGHS=zero
        alfaEff=zero
        return
      endif
      gg5=cropProp(6,crpIndex) !pi50 coefficient
      hh5=cropProp(4,crpIndex) !h50 coefficient
      rel_fact=cropProp(7,crpIndex)
      tPot=rel_fact*tPotMax
      lReduc=.FALSE.
      lallNodes=.TRUE.
      vRoot=zero
      do i=N-1,iRoot,-1
        se(i)=se0(i)
c
c Calculate the reduction function for the matric head, hRed,
c and salt concentration piRed for nodes with root uptake.
        alphaHS(i)=FAlfa(gg5,hh5,hNew(i),se(i))
        if (alphaHS(i).lt.Et_shres .and.
     ;      hNew(i).ge.cropProp(8,crpIndex)) then
            Sink(i)=Beta(i)*alphaHS(i)*TPot
            alphaHSR(i)=alphaHS(i)
            lReduc=.TRUE.
        elseif (hNew(i).lt.cropProp(8,crpIndex)) then
c           Matric head is less than PWP, uptake is zero
            Sink(i)=zero
            alphaHS(i)=zero
            alphaHSR(i)=zero
            lReduc=.TRUE.
        else
            Sink(i)=Beta(i)*TPot !No reduction if above threshold
            lAllNodes=.FALSE.
            alphaHSR(i)=one
        endif
c       get the new water contents
        theta(i)=theta0(i)-Sink(i)*dt
        vRoot=vRoot+Sink(i)*dz
      enddo
c     Calculate alpha with no compensation
      alfaGHS=vRoot/tPot
      if (.not.lReduc) then !no nodes had alphaHS < threshold
        alfaEff=alfaGHS     !no compenstion occurs
        return
      elseif (lAllNodes) then !All nodes had alphaHS < threshold
        alfaEff=alfaGHS       !no compensation occurs
        return
      else  !perform compensation calculation
        done = .FALSE.
        if (TPot.gt.vRoot) then
          diff = TPot - vRoot
        else
          done=.TRUE. !No iteration if vRoot=TPot
        endif
c
        iter=1
        do while (.not.done)
          done=.TRUE.
          vadd=zero
          do i=N-1,iRoot,-1
            h(i)=fH(IWMODEL,theta(i),ParD,MatNum(i))
c           calculate increase in concentration due to removal of water
            se(i)=se0(i)*theta0(i)/theta(i)
            reduc(i)=FAlfa(gg5,hh5,h(i),se(i))
c           Water is extracted only from those nodes with alphaHS > threshold
            if (reduc(i).gt.Et_shres) then
              newSink(i)=diff * Beta(i) * reduc(i)
              theta(i) = theta(i) - (newSink(i) * dt)
            else
              newSink(i)=zero
            endif
c
c           if any node exceeds threshold then further water uptake
c           can still occur so another iteration is needed
            if (reduc(i).gt.Et_shres) then
              done=.FALSE.
            endif
            vadd=vadd + newSink(i)*dz
          enddo
          vRoot = vRoot + vadd
c
c         if vRoot is less than TPot then continue iterations
          if (.not.done .and. ((one-vRoot/TPot-eps).gt.zero)) then
            diff = TPot - vRoot
          else
            done=.TRUE.
          endif
c
c         Set maximum limit on iterations of 15
          if (iter.ge.15) then
            done=.TRUE.
          else
            iter=iter+1
          endif
        enddo
c       Use the total change in theta to recalculate the sink term,
c       root uptake cannot be greater than potential uptake
        vRoot=zero
        do i=N-1,iRoot,-1
          Sink(i) = max((theta0(i)-theta(i))/dt,zero)
          if (Beta(i).eq.zero .or. TPot.eq.zero) then
            iter=0
            alphaHSR(i)=zero
          else
            alphaHSR(i)=Sink(i)/(Beta(i)*TPot)
          endif
          vRoot=vRoot+Sink(i)*dz
        enddo
        alfaEff=vRoot/TPot
      endif
      return 
      end

************************************************************************
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 SetRel(lNitrate, cropProp, dt, zR, TPotMax,TPotN,
     ;                  TPotNMax,sumNAct, sumNPot,tMassR, vRoot,
     ;                  sumEPot,sumRoot,crpIndex,sEPotMax,ierr)
      implicit none
      include "defs.fi"
      integer crpIndex,ierr
      logical lNitrate
      double precision cropProp,dt,TPotMax,TPotN,vRoot,sumNPot,
     ;                 sumNAct,tMassR,sumEPot,sumRoot,zR,sEPotMax,
     ;                 TPotNMax,diffN,sumDel
      dimension cropProp(nCropPar,*)
c
      double precision rel_fact,alphaW,alphaN
      double precision zero,one,eps
      data zero,one,eps / 0.0D00, 1.0D00, 1.0D-06 /
c
      ierr=0
c
c     Compute maximum possible water uptake
      sEPotMax = sEPotMax + TPotMax*dt
c
c     Compute sum of potential and actual root water uptake
      rel_fact=cropProp(7,crpIndex)
      alphaW = sumRoot/sumEPot
      sumRoot=sumRoot + vRoot*dt
      sumEpot = sumEpot + TPotMax*dt*rel_fact
      if (lNitrate.and.(TPotN.gt.zero)) then
        sumNPot = sumNPot + TPotNMax*dt*rel_fact
        sumNAct = sumNAct + tMassR
        alphaN = sumNAct/sumNPot
c       Only make a change in rel_fact if alphaN differs from
c       alphaW by a finite amount 
        if (abs(alphaN/alphaW - one).gt.eps) then
          rel_fact=min(alphaN,alphaW)
        endif
      elseif (lNitrate .and.(tMassR.le.zero)) then
        rel_fact=sumRoot/sumEPot
        alphaN=zero
      elseif (lNitrate .and.(tMassR.gt.zero)) then
        ierr=73 !fatal
      else !Rel_fact is determined by water uptake
        rel_fact=sumRoot/sumEPot
      endif
c
c     Only recalculate growth factor when rooting depth is
c     greater than 10 cm. Otherwise rel_fact is unstable.
      if (vRoot.eq.zero.or.zR.lt.1.0D01) then
        rel_fact = one
      endif
      cropProp(7,crpIndex)=min(rel_fact,one)
      return
      end

      double precision function FAlfa(gg5,hh5,hRed,piRed)
      implicit none
      double precision gg5,hh5,hRed,piRed
c
c     Local variables
      double precision ps,zero,one,three
      data zero,one,three /0.0D00, 1.0D00, 3.0D00/
c
      ps=min(gg5/hh5*hRed + piRed * (-3.6D01),zero)
      FAlfa=one/(one+(ps/gg5)**three)
      return
      end

      subroutine SetRG(N,x,Beta,xR,t,dz,cropProp,nPlants,crpIndex,
     ;                 cropSchd,cropID,nCrpRecs,iCrop,iRoot,ierr)
      implicit none
      include "defs.fi"
      integer N,nPlants,cropNum,nCrpRecs,iCrop,cropID(nPlants),
     ;        crpIndex,iRoot,ierr
      double precision x(N),Beta(N),xR,t,dz
      double precision cropProp(nCropPar,*),cropSchd(nSchdPar,*)
c
c     cropSchd:
c     1 - cropNum
c     2 - tRMin
c     3 - tRHarv
c     4 - RDFDay
c     5 - xRMin
c
c     cropProp:
c   cropProp(1,i) = cropNum
c   cropProp(2,i) = xRMax     maximum rooting depth
c   cropProp(3,i) = hShres
c   cropProp(4,i) = hh5
c   cropProp(5,i) = sShres
c   cropProp(6,i) = gg5
c   cropProp(7,i) = rel_fact
c   cropProp(8,i) = hPWP
c
c Local Variables
      double precision tRMin,tRHarv,RDFDay,xRMin
      double precision xRMax,RGR,aa,c,SBeta,SSBeta,temp
      double precision rtm,xRMed
      dimension temp(N)
      integer i
      double precision zero,one,two
      data zero,one,two / 0.0D00,1.0D00,2.0D00 /
c
      ierr=0
      cropNum=INT(cropSchd(1,iCrop))
      tRMin=cropSchd(2,iCrop)
      tRHarv=cropSchd(3,iCrop)
      RDFday=cropSchd(4,iCrop)
      xRMin=cropSchd(5,iCrop)
c     Get static crop params
      xRMax=cropProp(2,crpIndex)
      if(t.lt.tRMin.or.t.gt.tRHarv.or.cropNum.eq.0) then
        do i=1,N
          Beta(i)=zero
        enddo
        iRoot=0
        return
      end if
c
c     Test for data errors
      if (tRMin.ge.tRHarv .or. xRMin.gt.xRMax .or.
     ;    tRMin.gt.RDFDay  .or. RDFDay.gt.tRHarv) then
        ierr=31
        return
      endif
c     calculate xRMed using 99% rule
      xRMed=xRMin + 0.99D00*(xRMax-xRMin)
      rtm=RDFDay-tRMin
      RGR = -(1/rtm)*dlog((xRMin*(xRMax-xRMed))/(xRMed*(xRMax-xRMin)))
c
c     xRMin must be at least equal to 0.001
      xRMin=max(xRMin,1.0D-03)
      if(t.ge.tRMin.and.t.le.RDFday) then
        xR=(xRMax*xRMin)/(xRMin+(xRMax-xRMin)*exp(-RGR*(t-tRMin)))
      else if(t.gt.RDFday.and.t.le.tRHarv) then
        xR=xRMax
      endif
      do i=N,1,-1
        if (x(i).le.-xR) then
          iRoot=min(i+1,N-1)
          goto 15
        endif
      enddo
15    continue
c
      SBeta=zero
      do i=N-1,iRoot,-1
        if (x(i).ge.-xR/3.0D00) then
          Beta(i)=one
        elseif (x(i).ge.-xR) then
          Beta(i)=3.0D00*(x(i)+xR)/(two*xR)
        else
          Beta(i)=zero
        endif
        SBeta=SBeta+Beta(i)
      enddo
      do i=N-1,1,-1
        Beta(i)=Beta(i)/SBeta
      enddo
      return
      end !SetRG
c
c     Initialize accumulators for variables at the start
c     of a cropping season
      subroutine initSums(N,sumEPot,sumRoot,cropProp,
     ;                    crpIndex,alphaG,sumEff,sumAlfG,
     ;                    accETp,sNoComp,sumNAct,sumNPot,
     ;                    sPotN,stSeason,sEPotMax,sumDel)
      implicit none
      include "defs.fi"
      integer crpIndex,N
      logical stSeason
      double precision sumEPot,sumRoot,cropProp,alphaG,
     ;                 sumEff,sumAlfG,accETp,sNoComp,sumNAct,
     ;                 sumNPot,sPotN,sEPotMax,sumDel,
     ;                 tLeachN
      dimension cropProp(nCropPar,*)
c
      integer i
      double precision zero,one
      data zero,one /0.0D00, 1.0D00/
c
      stSeason=.True.
      sumEPot = zero
      sumRoot = zero
      accETp = zero
      sEPotMax = zero
      sumDel=zero
c     Initialize at start of crop only
      cropProp(7,crpIndex) = one
c     Set alphaG = 1 because SetSnk will run for the first
c     time before SetSnkN so reduction coefficient should be one
      alphaG=one
      sumEff=zero
      sumAlfG=zero
      sNoComp=zero
      sumNAct=zero
      sumNPot=zero
      sPotN=zero
      return
      end

c     Initialize daily sums for various reduction coefficients
      subroutine initDSum(sumEff,sumAlfG,sNoComp)
      implicit none
      double precision sumEff,sumAlfG,sNoComp
      double precision zero
      data zero /0.0D00/
      sumEff=zero
      sumAlfG=zero
      sNoComp=zero
      return
      end

      integer function getCIndx(nPlants,cropNum,cropProp,ierr,errStr)
      implicit none
      include "defs.fi"
      integer nPlants,cropNum,crpIndex,ierr
      double precision cropProp
      dimension cropProp(nCropPar,*)
      character*80 errStr
c
      integer i
      ierr=0
c
c     Total crop records is nPlants+2 because of the fallow record (cropID=0)
c     and the cover crop record (cropID=1)
      do i=1,nPlants+2
        if (cropNum.eq.INT(cropProp(1,i))) then
          getCIndx=i
          goto 100
        endif
      enddo
      errStr='Plant.in'
      ierr=6
100   return
      end

************************************************************************
c
c     Calculate the N concentration resulting from the decay of organic
c     N fertilizer. Also recalculate Psi, the fraction of organic-N that
c     is in the slow pool. Note CFix(1,:) are the fixed concentrations
c     of total organic N. CFix(2,:) are the inorganic fixed concentrations
c     which may be present due to incorporation of inorganic N in the soil.
      subroutine setSrcN(N,NMat,volConv,x,MatNum,t,dt,CFix,
     ;                    nFixC,Conc,nConc,chemPar,nChemPar,
     ;                    Psi,theta,NLossRat,sConvN,TVarCoef,
     ;                    lTCoef,dateInit,prodNSea,nMonths)
      implicit none
      integer N,iRoot,nFixC,nConc,MatNum,NMat,nChemPar,nMonths
      logical lTCoef
      character*10 dateInit
      double precision volConv,x,t,dt,CFix,Conc,NLossRat,Psi,
     ;                 theta,chemPar,sConvN,TVarCoef,prodNSea
      dimension x(N),Conc(nConc,N),CFix(nFixC,N),theta(N),Psi(N),
     ;          chemPar(nChemPar,NMat),MatNum(N),TVarCoef(nMonths)
c
      integer i
      double precision one,MWN,delN,MassN,dxM,rhoBulk,lambdaF,TCoef,
     ;                 lambdaS,MassNF,MassNS,deltaNF,deltaNS,NTReduc,
     ;                 lambdaS0,eps
      external NTReduc
c     lambdaS is for 1% decay per year, half-life 69 years
      data one,MWN,lambdaS0,eps / 1.0D00,1.4D01,2.753D-05,1.0D-06 /
c
c     Multiply NLossRat by the monthly temperature reduction coefficient
c     if such coefficients were defined in Basic.IN
      if (lTCoef) then
        TCoef = NTReduc(t,dateInit,TVarCoef)
        lambdaF = NLossRat * TCoef
        lambdaS = lambdaS0 * TCoef
      else
        lambdaF = NLossRat
        lambdaS = lambdaS0
      endif
      do i=N-1,1,-1
        dxM = x(i+1)-x(i)
        rhoBulk=chemPar(2,MatNum(i))
c       Calculate decay only for those elements with non-negligible
c       concentration of organic-N
        if (CFix(1,i).gt.eps) then
          MassN=CFix(1,i)*rhoBulk*dxM  !Mass of org N in g in element i
          MassNF=(one-Psi(i))*MassN
          MassNS=Psi(i)*MassN
          deltaNF=MassNF*(one - exp(-lambdaF*dt)) !g
          deltaNS=MassNS*(one - exp(-lambdaS*dt))
          Psi(i) = (MassNs - deltaNs)/
     ;             (MassNF-deltaNF + MassNs - deltaNs)
          delN = deltaNF + deltaNS
c         Conc in solution in mmol/L is mass in g divided by
c         product of MW, volume conversion, theta and dx
          Conc(3,i)=Conc(3,i)+delN*1.D03/(MWN*volConv*theta(i)*dxM)
c         Subtract the amount lost from organic N
          CFix(1,i)=CFix(1,i)-delN/rhoBulk/dxM !g/g Soil
c         sum of mass of N converted from organic to NO3-N (mmol)
          sConvN=sConvN + delN*1.D03/MWN
          prodNSea=prodNSea + delN*1.D03/MWN
        endif
      enddo
c     set top node
      if (CFix(1,N-1).gt.eps) then
        cFix(1,N)=cFix(1,N-1)
        Psi(N)=Psi(N-1)
      endif
      return
      end

      double precision function sumSinkL(N,NLayers,Sink,nodeTop,
     ;                                   nodeBot,dt)
      implicit none
      integer N,NLayers,nodeTop,nodeBot
      double precision Sink,dt
      dimension Sink(N)
c
      double precision zero
      integer j
      data zero /0.0D00/
      sumSinkL=zero
      do j=nodeTop,nodeBot,-1
        sumSinkL=sumSinkL + Sink(j)*dt
      enddo
      return
      end