***********************************************************************
c     SetSnkN computes root uptake of NO3-N
      subroutine SetSnkN(N,iRoot,volConv,t,dt,x,tpotN,tpotMax,
     ;                   tPotNMax,Conc,tMassR,theta,Beta,alphaG,
     ;                   cStarN,sumDel,cropProp,crpIndex,
     ;                   ierr,compCnt)
      implicit none
      include "defs.fi"
      integer N,iRoot,crpIndex,ierr,compCnt
      double precision tpotN,x,tPotMax,Conc,theta,Beta,tMassR,
     ;                 volConv,alphaG,dt,t,cStarN,cropProp,
     ;                 tPotNMax,sumDel
      dimension Conc(nConcD,N),x(N),Beta(N),cropProp(nCropPar,*)
c
c     Local variables
      integer i,j,iter,nodeCnt
      double precision Cstar,sumN,dxM,Ci,PMass,cNew,delN,
     ;  mInitN,mFinalN,volSum,rel_fact,sumMass,reqMass,
     ;  diff,sumTemp,Weight,tempMass
      logical nodeMark,deficit,done
      dimension dxM(N),PMass(N),cNew(N),theta(N),Ci(N),Weight(N),
     ;          nodeMark(N),sumMass(N),reqMass(N),tempMass(N)
      double precision zero,one,two,eps,sPMass,half,small,
     ;                 err,stSumN,sumBeta,totReq,totalM,CiplsOld,
     ;                 cTemp,getMassN
      external getMassN
      data zero,one,two,eps /0.0D00, 1.0D00, 2.0D00, 1.0D-15/
      data half,small /5.D-01,1.0D-12/
c
      if (abs(tPotMax).lt.eps.or.abs(tpotN).lt.eps) then
c       No uptake is occurring set reduction factor alphaG=1
c       and return immediately
        alphaG=one
        tMassR=zero
        return
      endif
c     Test the user-entered value of CstarN, if it's greater
c     than zero then the user is providing the value. Otherwise
c     CstarN must be calculated as the ratio of potential
c     nitrate uptake to potential water uptake.
      if (CstarN.gt.zero) then
        Cstar=CstarN
      else
        Cstar=tpotN/tPotMax !mmol/cm^2/d  / cm^3/cm^2/d
        Cstar=Cstar*1.0D03  !mmol/L
      endif
      rel_fact=cropProp(7,crpIndex)
c
c     Total N required is mmol/cm^2/d * d is mmol/cm^2
      tPotN=rel_fact*tPotNMax
      totReq=tPotN*dt
      sumN=zero
      tMassR=zero
      do i=N-1,1,-1
        dxM(i)=x(i+1)-x(i)
        Ci(i)=Conc(3,i) !mmol/L
        ReqMass(i)=zero
        PMass(i)=zero
      enddo
      iRoot=max(iRoot,1)
c
c     Test for deficiency and make sure Beta is non-zero
      deficit=.TRUE.
      sumBeta=zero
      mInitN=getMassN(N,iRoot,Conc,theta,dxM,volConv)
      do i=N-1,iRoot,-1
        cNew(i)=Ci(i)
        if (Ci(i).gt.Cstar) then
          deficit=.FALSE.
          nodeMark(i)=.TRUE.
        else
          nodeMark(i)=.FALSE.
        endif
        sumBeta=sumBeta+Beta(i)
      enddo
c
c     All Beta equal to zero means no uptake: return
      if (sumBeta.le.zero) then
        mFinalN=getMassN(N,iRoot,Conc,theta,dxM,volConv)
        alphaG=one
        return
      endif

      if (.not.deficit) then
c       Get N from nodes
        reqMass(N)=zero
        sumMass(N)=zero
        CiPlsOld=zero
        do i=N-1,iRoot,-1
          sumMass(i)=zero
          if (nodeMark(i)) then
            ReqMass(i)=totReq*Beta(i)
            PMass(i)=volConv*(Ci(i)-Cstar)*dxM(i)*theta(i)
            if (PMass(i).ge.ReqMass(i)) then
              call resCNode(cNew(i),Ci(i),ReqMass(i),theta(i),dxM(i),
     ;                      volConv)
              sumN=sumN+ReqMass(i)
              sumMass(i)=ReqMass(i)
            elseif (PMass(i).gt.zero) then
              cNew(i)=Cstar
              sumN=sumN+PMass(i)
              sumMass(i)=PMass(i)
c           Otherwise this node has a deficit
            endif
          endif
        enddo
        err=abs(one - sumN/(totReq))-small
        if (err.le.zero) then
c         the uptake is complete reset the concentrations
          call resetC(N,iRoot,sumMass,Conc,theta,volConv,dxM)
          tMassR=sumN
        elseif (sumN.gt.totReq) then
c         this is an error because required mass reqMass should never
c         exceed TPotN*dt by a finite amount
          ierr=72
          return
        else
c         there was insufficient mass to meet demand so the remainder
c         is further removed from those nodes with Ci > C*
c         This is the implementation of compensation.
          done=.FALSE.
c         Initialize nodeMark to FALSE, Weight to zeros
          call initMark(N,nodeMark,Weight)
          iter=1
          stSumN=sumN
          diff=totReq - sumN
c         Start iterative loop for compensation
          do while (.not.done)
            sumTemp=zero
c           Identify participating nodes for weighting function
            sumBeta=zero
            do i=N-1,iRoot,-1
              nodeMark(i)=.FALSE.
              if (cNew(i).gt.CStar) then
                nodeMark(i)=.TRUE.
                sumBeta=sumBeta+Beta(i)
              endif
            enddo
            do i=N-1,iRoot,-1
              tempMass(i)=zero
              if (nodeMark(i)) then
                Weight(i)=Beta(i)/sumBeta
              endif
            enddo
            sumMass(N)=zero
            do i=N-1,iRoot,-1
              if (nodeMark(i)) then
C               diff already corrected for deltat so only need
c               to normalize using the weighting function
                ReqMass(i)=diff*Weight(i)
                PMass(i)=
     ;              volConv*dxM(i)*((CNew(i)-CStar)*theta(i)+
     ;              (CNew(i+1)-CStar)*theta(i+1))/two
                if (PMass(i).ge.ReqMass(i)) then
                  sumTemp=sumTemp + reqMass(i)
                  sumMass(i)=sumMass(i) + reqMass(i)
                  tempMass(i)=reqMass(i)
                elseif (PMass(i).gt.zero) then
                  cNew(i)=Cstar
                  sumTemp=sumTemp+PMass(i)
                  sumMass(i)=sumMass(i) + PMass(i)
                  tempMass(i)=PMass(i)
                endif
              endif
            enddo
c
c           If no node had excess N then there is no further
c           compensation and the loop is done, go to N-deficiency
c           calculation for remaining N uptake
            if (sumTemp.le.zero) then
              done=.TRUE.
c             return total mass of N uptake in tMassR
              totalM=zero
              do i=N-1,iRoot,-1
                totalM = totalM + sumMass(i)
              enddo
              call defCon(N,cNew,theta,iRoot,alphaG,tMassR,
     ;                Cstar,volConv,dxM,diff,sumMass,spMass,
     ;                sumN,ierr)
              tMassR=spMass
            else
              sumN=sumN + sumTemp
c             maximum 50 iterations for compensated N uptake
              if (sumN.ge.diff) then
                done=.TRUE.
              elseif (iter.gt.50) then
                ierr=73
c               increment counter for iteration overrun, used
c               to limit the number of lines of error output
                compCnt=compCnt+1
                done=.TRUE.
              else
                iter=iter+1
                diff=totReq - sumN
              endif
            endif
c           Initialize nodemarks and weights and reset the cNew
c           vector of concentrations for the next iteration
            if (.not.done) then
              do i=N-1,iRoot,-1
                Ci(i)=cNew(i)
              enddo
              call initMark(N,nodeMark,Weight)
c             Reset the cNew vector
              call resCNew(N,iRoot,tempMass,cNew,Ci,theta,volConv,
     ;                     dxM)
            else
              tMassR=sumN
            endif
          enddo !end of iteration loop
          call resetC(N,iRoot,sumMass,Conc,theta,volConv,dxM)
        endif
      else
c       deficit conditions, N uptake at all nodes subject
c       to nitrogen deficiency
        do i=N,iRoot,-1
          sumMass(i)=zero
        enddo
        call defCon(N,Ci,theta,iRoot,alphaG,tMassR,Cstar,
     ;              volConv,dxM,totReq,sumMass,spMass,
     ;              sumN,ierr)
c       Reset N-component of Conc Array
        call resetC(N,iRoot,sumMass,Conc,theta,volConv,dxM)
        spMass=zero
        do i=N-1,iRoot,-1
          spMass=spMass+sumMass(i)
        enddo
      endif
c
c     After N uptake recalculate N mass in root zone in mmol/cm^2
      mFinalN=getMassN(N,iRoot,Conc,theta,dxM,volConv)
c     Testing of actual total N uptake
      sumDel=sumDel + mInitN - mFinalN
c     Calculate alphaG for the entire N uptake including both
c     uptake under non-stressed and stressed conditions
      alphaG=tMassR/totReq
      return
      end

      subroutine defCon(N,cNew,theta,iRoot,alphaG,tMassR,
     ;                  Cstar,volConv,dxM,delN,sumMass,
     ;                  sPMass,sumN,ierr)
      implicit none
      integer N,iRoot,ierr
      double precision cNew,theta,alphaG,tMassR,Cstar,
     ;                 volConv,dxM,delN,sumMass,spMass,sumN
      dimension cNew(N),theta(N),dxM(N),sumMass(N)
c     
c         Root-finding operation to obtain cg
c         First, test the minimum cg value to make sure that
c         there is a negative result from getDelta. If the 
c         result is positive then Nuptake is zero and the
c         concentrations remain unchanged.
      integer i
      double precision addMass(N),totMass
      double precision cg,f,getDelta,fmid,rtbis,dx,xacc,Ci,
     ;                 CiPlus,zero,two
      data zero,two /0.0D00,2.0D00/

      cg = 1.0D-02 !minimum value of cg
      f = getDelta(cNew,theta,cg,iRoot,N,
     ;                  Cstar,volConv,dxM,delN)
      if (f.ge.zero) then
        goto 80 !special case all cNew < cg
      endif
      cg=Cstar
      fmid=getDelta(cNew,theta,cg,iRoot,N,
     ;              Cstar,volConv,dxM,delN)
      cg=1.D-02
      f=getDelta(cNew,theta,cg,iRoot,N,
     ;           Cstar,volConv,dxM,delN)
      if (f*fmid.gt.zero) then
        ierr=70
        return
      endif
      if (f.lt.zero) then
        rtbis=zero
        dx=Cstar
      else
        rtbis=Cstar
        dx=-Cstar
      endif
      xacc=1.0D-12
      do while (abs(dx).gt.xacc.and.fmid.ne.zero)
        dx=dx*0.5D00
        cg=rtbis + dx
        fmid=getDelta(cNew,theta,cg,iRoot,N,
     ;                Cstar,volConv,dxM,delN)
        if (fmid.le.zero) then
          rtbis=cg
        endif
      enddo
      alphaG = cg/Cstar
      totMass=zero
      spMass = zero
      if (alphaG.gt.1.D-02) then
        do i=N-1,iRoot,-1
          Ci=cNew(i)
          CiPlus=cNew(i+1)
          addMass(i)=max(Ci-cg,zero)*theta(i)*volConv*dxM(i)
          totMass=totMass+addMass(i)
          sumMass(i)=sumMass(i)+addMass(i)
          spMass=spMass+sumMass(i)
        enddo
      endif
c     sumN may be nonzero on entry due to earlier N uptake
      sumN=sumN + totMass
      tMassR = sumN
      return
80    alphaG=zero
      tMassR=zero
      return
      end

c     For the estimated cg return the difference between the
c     global reduction in delN by the ratio cg/Cstar and the
c     sum of N uptake that results from the choice of cg
      double precision function getDelta(Conc,theta,cg,iRoot,
     ;                          N,Cstar,volConv,dxM,delN)
      implicit none
      integer iRoot,N
      double precision Conc,theta,cg,Cstar,volConv,dxM,delN
      dimension Conc(N),theta(N),dxM(N)
c
      integer i
      double precision sumN,zero,two,Ci,CiPlus
      data zero,two /0.0D00, 2.0D00/
      sumN=zero
      do i=N-1,iRoot,-1
        Ci=Conc(i)
        if (i.eq.N-1) then
          sumN = max(Ci-cg,zero)*theta(i)*volConv*dxM(i)
        else
          CiPlus=Conc(i+1)
          sumN = sumN + max(Ci-cg,zero)*theta(i)*volConv*dxM(i)
        endif
      enddo
      getDelta=cg*delN/Cstar - sumN
      return
      end
c
c     reset the concentration vector after N uptake
      subroutine resetC(N,iRoot,sumMass,Conc,theta,volConv,dxM)
      implicit none
      include 'defs.fi'
      integer N,iRoot
      double precision sumMass,Conc,dxM,theta,volConv
      dimension sumMass(N),Conc(NConcD,N),dxM(N),theta(N)
c
      integer i
      double precision delta,one,eps
      data one,eps /1.0D00,1.0D-15/
c
      do i=N-1,iRoot,-1
        delta=sumMass(i)/(volConv*theta(i)*dxM(i))
        if (abs(Conc(3,i)/delta - one).gt.eps) then
          Conc(3,i)=Conc(3,i) - delta
        endif
      enddo
      return
      end

      subroutine resCNew(N,iRoot,tempMass,CNew,C,theta,volConv,dxM)
      implicit none
      integer N,iRoot
      double precision tempMass(N),C(N),CNew(N),theta(N),volConv,
     ;                 dxM(N)
c
      integer i
      double precision zero
      data zero /0.0D00/
c
      tempMass(N)=zero
      do i=N-1,iRoot,-1
        call resCNode(CNew(i),C(i),tempMass(i),theta(i),dxM(i),
     ;                volConv)
      enddo
      return
      end

      subroutine resCNode(CNew,C,mass,thetaI,dx,volConv)
      implicit none
      double precision C,CNew,mass,thetaI,dx,volConv
c
      double precision one,eps,delta
      data one,eps /1.0D00,1.0D-15/
c
      delta=mass/(volConv*thetaI*dx)
      if (abs(C/delta - one).gt.eps) then
        CNew=C - delta
      else
        CNew=C
      endif
      return
      end

      double precision function getMassN(N,iRoot,Conc,theta,
     ;                                   dxM,volConv)
      implicit none
      integer N,iRoot
      double precision Conc,theta,dxM,volConv
      dimension Conc(3,N),theta(N),dxM(N)
c
      integer i
      double precision zero
      data zero /0.0D00/
c
      getMassN=zero
      do i=N-1,iRoot,-1
        getMassN=getMassN + Conc(3,i)*theta(i)*dxM(i)*volConv
      enddo
	  return
      end

c
c     Calculate the mass of organic-N in each layer
      subroutine layNMass(N,NLayers,NMat,MatNum,dz,chemPar,CFix,
     ;                    layerDep,Mass,ierr)
	  implicit none
      include 'defs.fi'
	  integer N,NLayers,NMat,MatNum,ierr
	  double precision dz,chemPar,CFix,layerDep,Mass
      dimension chemPar(nChemPar,NMat),layerDep(NLayers),
     ;          mass(NLayers),CFix(NFixC,N),MatNum(N)
c
      integer i,layNum,layOfNod
      double precision rhoBulk,zero
      data zero /0.0D00/
c
      do i=NLayers,1,-1
        Mass(i)=zero
      enddo
      do i=N-1,1,-1
        rhoBulk=chemPar(2,MatNum(i)) !g/cm^3
        layNum = layOfNod(N,layerDep,dz,i,NLayers,ierr)
        Mass(layNum)=Mass(layNum) + CFix(1,i)*rhoBulk*dz !grams
      enddo
      return
      end

	  subroutine initMark(N,nodeMark,Weight)
      implicit none
      integer N
      logical nodeMark(N)
      double precision Weight(N)
c
      integer i
      double precision zero
      data zero /0.0D00/
c
      do i=1,N
        nodeMark(i)=.FALSE.
        Weight(i)=zero
      enddo
      return
      end