c     ENVIRO-GRO text output routines
c
c     Open all standard output files
      subroutine outOpen(datPath,lNitrate,lPlant,ierr)
      implicit none
      include "defs.fi"
      character*120 datPath
      character*149 fname
      logical lNitrate,lPlant
      integer strlen,lenPath,ierr
      external strlen
      lenPath=strlen(datPath)
      fname = datPath(1:lenPath)//'\WatProf.dat'
      open(UNIT=nProfUnit,FILE=fname,ERR=900)
      fname = datPath(1:lenPath)//'\WatBal.dat'
      open(UNIT=nWBalUnit,FILE=fname,ERR=900)
      fname = datPath(1:lenPath)//'\DailyFlux.dat'
      open(UNIT=nDayUnit,FILE=fname,ERR=900)
      fname = datPath(1:lenPath)//'\ChemProf.dat'
      open(UNIT=nChemUnit,FILE=fname,ERR=900)
      fname = datPath(1:lenPath)//'\ChemBal.dat'
      open(UNIT=nCBalUnit,FILE=fname,ERR=900)
      if (lPlant) then
        fname = datPath(1:lenPath)//'\Reduc.dat'
        open(UNIT=nAlfUnit,FILE=fname,ERR=900)
        fname = datPath(1:lenPath)//'\CropFinal.dat'
        open(UNIT=nRelUnit,FILE=fname,ERR=900)
        fname = datPath(1:lenPath)//'\StressProf.dat'
        open(UNIT=nStrsUnit,FILE=fname,ERR=900)
      endif
      if (lNitrate) then
        fname = datPath(1:lenPath)//'\NO3.dat'
        open(UNIT=nNO3Unit,FILE=fname,ERR=900)
        fname = datPath(1:lenPath)//'\OrgNMass.dat'
        open(UNIT=nNRtUnit,File=fname,ERR=900)
      endif
      return
900   ierr=1
      return
901   ierr=2
      return
      end

c     Output to screen
      subroutine sreenOut(N,t,rTop,vTop,vRoot,hTop,hBot,lineCnt,hNew,
     ;                    KodBot,KodTop)
      implicit none
      integer lineCnt,N,KodBot,KodTop
      double precision t,rTop,vTop,vRoot,hTop,hBot,hNew
      dimension hNew(N)
c
      if (lineCnt.eq.20) then
c       print header line 
        write(*,FMT=110)
        write(*,FMT=120)
        lineCnt=1
      else
        lineCnt=lineCnt+1
      endif
      if (KodBot.gt.0.and.KodTop.gt.0) then
        write(*,FMT=100) t,rTop,vTop,vRoot,hTop,hBot
      elseif (KodBot.gt.0.and.KodTop.lt.0) then
        write(*,FMT=100) t,rTop,vTop,vRoot,hNew(N),hBot
      elseif (KodBot.lt.0.and.KodTop.gt.0) then
        write(*,FMT=100) t,rTop,vTop,vRoot,hTop,hNew(1)
      else
        write(*,FMT=100) t,rTop,vTop,vRoot,hNew(N),hNew(1)
      endif
      return
100   format(F8.1,2X,3(F10.1,3X),2(F10.1,4X))
110   format(
     ;'     Time       rTop       vTop       vRoot       hTop',
     ;'       hBot'/)
120   format(
     ;'     (d)       (cm/d)     (cm/d)     (cm/d)       (cm)')
      end
c
c     Create headers for various output files
      subroutine headers(N,CropProp,NLayers,lNitrate,lPlant,
     ;                   x,Conc,theta,ierr)
      implicit none
      include "defs.fi"
      integer N,NLayers,ierr
      logical lNitrate,lPlant
	  double precision x,Conc,theta
      dimension x(N),Conc(NConcD,N),theta(N)
      double precision cropProp,prntVec,zero
      dimension cropProp(nCropPar,*)
      integer prntLen,i
      character*80 errStr
      data prntLen /5/
      dimension prntVec(NPrntD)
      data zero /0.0D00/
c
      do i=1,NPrntD
        prntVec(i)=zero
      enddo
c     Headers are created by the normal output subroutine when
c     the subroutine is called with a TRUE value for the lHead
c     dummy logical variable. The other data are irrelevant for
c     the call.
      call dayPrnt(zero,zero,zero,zero,zero,zero,zero,0,.TRUE.,
     ;             ierr,errStr)
      if (lPlant) then
        call alfaPrn(.TRUE.,zero,cropProp,0,zero,zero,zero,0,ierr)
        call wrtCrop(.TRUE.,zero,0,'',zero,zero,zero,zero,zero,
     ;               zero,zero,lNitrate,ierr)
      endif
      if (lNitrate) then
        call NPrnt(N,zero,zero,zero,.TRUE.,.TRUE.,zero,zero,
     ;           zero,zero,x,zero,Conc,theta,zero,ierr,errStr)
        if (NLayers.eq.1) then
          call orgNPrnt(.TRUE.,4,prntVec)
        else
          call orgNPrnt(.TRUE.,7,prntVec)
        endif
      endif
      return
      end

c     Output profiles of variables related to water flow at fixed
c     times.
      subroutine FixPrnt(N,t,tPrint,x,MatNum,hNew,ThNew,Con,Sink,vNew,
     ;                   ierr,errStr)
      implicit none
      include "defs.fi"
      double precision t,x,tPrint,hNew,ThNew,Con,Sink,vNew
      integer N,MatNum,ierr
      character*80 errStr
      dimension x(N),hNew(N),ThNew(N),Con(N),Sink(N),vNew(N),MatNum(N)
c
      integer i
      errStr='WatProf.dat'
      write(nProfUnit,100,err=901) tPrint
      write(nProfUnit,105,err=901) 
      do i=N,1,-1
        write(nProfUnit,110,err=901) i,x(i),MatNum(i),hNew(i),ThNew(i),
     ;                               Con(i),vNew(i),Sink(i)
      enddo
      return
901   ierr=10
      return
100   format('Time =',f8.2,' (d)')
105   format(
     ;'Node     Depth  Mater   h       Theta',
     ;'     K        v        Sink',/
     ;' -       (cm)     -    (cm)       -  ',
     ;'   (cm/d)   (cm/d)     (1/d)')
110   format(I5,F8.2,2X,I3,2X,F8.1,1X,F8.3,1X,2(1PE10.2),E10.3)
      end
c
c     Output profiles of Sink term, uncompensated alpha and
c     compensated alpha at fixed times.
      subroutine strsPrnt(N,t,tPrint,x,Sink,AlphaHS,AlphaHSR,
     ;                    cropNum,cropName,iCrop,ierr,errStr)
      implicit none
      include "defs.fi"
      double precision t,x,tPrint,Sink,AlphaHS,AlphaHSR
      integer N,MatNum,cropNum,iCrop,ierr
      character*80 errStr
      character*20 cropName
      dimension x(N),Sink(N),AlphaHS(N),AlphaHSR(N)
c
      integer i
c
      errStr='StressProf.dat'
      write(nStrsUnit,100,err=901) tPrint,iCrop,cropNum,cropName
      write(nStrsUnit,105,err=901) 
      do i=N,1,-1
        write(nStrsUnit,110,err=901) i,x(i),Sink(i),AlphaHS(i),
     ;                               AlphaHSR(i)
      enddo
      return
901   ierr=10
      return
100   format('Time = ',f8.2,' (d)',3X,'Crop#',I3,'  CropID=',
     ;       I3,5X,A10)
105   format(
     ;'Node     Depth       Sink         Alpha    AlphaComp'/
     ;' -        (cm)       (1/d)         -           -')
110   format(I5,1X,F10.1,3X,E12.3,2X,F8.2,3X,F8.2)
115   format(//)
      end
c
c     Output water balance for the period between fixed output times.
c     Water balance is the residual of the summed water flow into the profile
c     minus summed water flow out and minus the change in water storage.
c     WatBal handles cases where only a single layer has been defined.
      subroutine watBal(N,z,t,thStart,thNew,accVTop,accVBot,accSink,
     ;                  totInfil,hFlag,ierr,errStr)
      implicit none
      include "defs.fi"
      integer N,ierr
      logical hFlag
      character*80 errStr
      double precision z,t,thStart,thNew,accVTop,accVBot,accSink,
     ;                 totInfil
      dimension z(N),thStart(N),thNew(N)
c
      integer i
      double precision zero,two,delStore,dz,wBalErr
      data zero,two /0.0D00, 2.0D00/
c
      ierr=0
      errStr='WatBal.dat'
c     write the header to the output file
      if (hFlag) then
        write(nWBalUnit,200,err=901)
        write(nWBalUnit,210,err=901)
        return
      endif
      delStore=zero
      do i=N-1,1,-1
        dz=z(i+1)-z(i)
        delStore = delStore + (thNew(i+1)+thNew(i)-
     ;                         thStart(i+1)-thStart(i))*dz/two
      enddo
c     sign convention: when vTop is negative the accumulating flux is added to 
c                      the change in storage so the sign on accVTop is minus
      wBalErr = -accVTop + accVBot - delStore - accSink

      write(nWBalUnit,FMT=100,err=901) t,accVTop,totInfil,accVBot,
     ;      accSink,delStore,wBalErr
c
c     reset the accumulators and thStart
      call veccpr(1,N,ThNew,thStart,ierr)
      call initBal(accVTop,accVBot,accSink,totInfil)
      return
901   ierr=10
      return
100   format(F8.2,2X,4(F8.2,2X),F8.2,2X,E10.2E3)
***********************************************************************
200   format(
     ;'  Time      sumVTop  totInfil   sumVBot   sumSink   Storage',
     ;'   Error')
210   format(
     ;'  (d)        (cm)      (cm)      (cm)       (cm)     (cm)  ',
     ;'   (cm)',/)
      end
c
c     layerBal does the same thing as WatBal but provides water
c     balance for each layer in the profile.
      subroutine layerBal(t,N,NLayers,z,layerDep,thStart,thNew,
     ;                    sLayTopQ,sLayBotQ,sLaySink,ierr,errStr)
      implicit none
      integer N,NLayers,ierr
      double precision z,layerDep,thStart,thNew,sLayTopQ,sLayBotQ,
     ;                 sLaySink,t
      character*80 errStr
      dimension z(N),layerDep(NLayers),sLayTopQ(NLayers),
     ;          sLayBotQ(NLayers),sLaySink(NLayers),thNew(N),
     ;          thStart(N)
c
      integer j,layNum,layOfNod,layLast
      double precision dz,delStore,relErr,wBalErr,watChng
      dimension delStore(NLayers)
      double precision zero,one,two
      external layOfNod
      data zero,one,two /0.0D00,1.0D00,2.0D00/
c
      ierr=0
      errStr='WatBal.dat'
      do j=NLayers,1,-1
        delStore(j)=zero
      enddo
      layLast=layOfNod(N,layerDep,z(N)-z(N-1),N,NLayers,ierr)
c     get the change in storage for each layer
      do j=N-1,1,-1
        dz = z(j+1)-z(j)
        layNum = layOfNod(N,layerDep,dz,j,NLayers,ierr)
        delStore(layNum) = delStore(layNum) + 
     ;      (thNew(j+1)+thNew(j)-thStart(j+1)-thStart(j))*dz/two
      enddo
      
      do j=NLayers,1,-1
c           This section only runs when the balance needs to
c           be calculated for layers above the bottom layer
            wBalErr = -sLayTopQ(j) + sLayBotQ(j) - 
     ;                delStore(j) - sLaySink(j)
c           find the largest component of the water balance: storage
c           flux or sink
            watChng = max(abs(delStore(j)),
     ;                    abs(sLayBotQ(j)-sLayTopQ(j)),
     ;                    abs(sLaySink(j)))
            if (watChng.gt.zero) then
              relErr = 1.D02*abs(wBalErr)/abs(watChng)
            else
              relErr = -one
            endif
            call layBPrnt(t,NLayers,j,sLayBotQ(j),
     ;                    sLayTopQ(j),slaySink(j),delStore(j),
     ;                    wBalErr,relErr,ierr)
      enddo
      return
      end
c
c     Subsidiary subroutine to layerBal for output of data for
c     individual layers
      subroutine layBPrnt(t,NLayers,layNum,sLayBotQ,sLayTopQ,sLaySink,
     ;                    delStore,wBalErr,relErr,ierr)
      implicit none
      include "defs.fi"
      integer NLayers,layNum,ierr
      double precision sLayBotQ,sLayTopQ,sLaySink,delStore,wBalErr,
     ;                 relErr,t
c
      double precision zero
      if (layNum.eq.NLayers) then
        write(nWBalUnit,100,ERR=901) t
        write(nWBalUnit,110,ERR=901)
      endif

      if (relErr.lt.zero) then
        write(nWBalUnit,115,ERR=901) layNum,sLayBotQ,sLayTopQ,
     ;                       sLaySink,delStore,wBalErr
      else
        write(nWBalUnit,120,ERR=901) layNum,sLayBotQ,sLayTopQ,
     ;                       sLaySink,delStore,wBalErr,relErr
      endif
      if (layNum.eq.1) then
        write(nWBalUnit,150,ERR=901)
        write(nWBalUnit,160,ERR=901)
      endif
      return
901   ierr=10
      return
100   format('Time =',F7.2,' (d)')
110   format('Layer# QCumulBot QCumulTop TotUptake delStore',
     ;' ABS_ERR REL_ERR(%)')
115   format(1X,I3,2X,2(F7.2,2X),F7.2,6X,F7.2,2X,E8.1E2,1X,'--N/A--')
120   format(1X,I3,2X,2(F7.2,2X),F7.2,6X,F7.2,2X,E8.1E2,1X,F6.1)
150   format(
     ;'--------------------------------------------------------------')
160   format(/)
      end
c
c     Output data representing each single day of the simulation.
c     The variables represent 24-hour averages and are reinitialized
c     after the output occurs.
      subroutine dayPrnt(t,vTopB,vRootB,rTopB,vBotB,eTopB,sumDt,
     ;                   nSteps,lHead,ierr,errStr)
      implicit none
      include "defs.fi"
      integer nSteps,ierr
      logical lHead
      character*80 errStr
      double precision t,vTopB,vRootB,rTopB,vBotB,sumDt,eTopB
c
      double precision zero
	  integer i
      data zero /0.0D00/
c
      errStr='DailyFlux.dat'
      if (lHead) then
        write(nDayUnit,110,err=901)
      else
        write(nDayUnit,100,err=901) t,sumDt/dble(nSteps),rTopB,vTopB,
     ;                              vRootB,vBotB,eTopB
      endif
c
c     Reset all the mean values to zero to start next day
      vTopB=zero
      vRootB=zero
      rTopB=zero
      vBotB=zero
      nSteps=0
      sumDt=zero
      eTopB=zero
      return
901   ierr=10
      return
100   format(F8.3,1X,E10.2,1X,3(F8.2,1X),2(F8.3,1X))
110   format(
     ;' Time        dt        rTop      vTop    vRoot',
     ;'    vBot    Evap'/
     ;'  (d)       (d)       (cm/d)    (cm/d)   (cm/d)',
     ;'   (cm/d)  (cm/d)'/)
      end


***********************************************************************
c     Update the 24-hour mean values for the current time step
      subroutine updOut(dt,vTop,vBot,vRoot,rTop,vTopB,vRootB,rTopB,
     ;                  eTopB,vBotB,wDayEvap)
      implicit none
      double precision dt,vTop,vBot,vRoot,rTop,vTopB,vRootB,rTopB,
     ;                 vBotB,eTopB,wDayEvap
      double precision zero
      data zero /0.0D00/
c
c     The unit timestep is one day so each variable is 
c     multiplied by the current fraction of a day represented
c     by dt. After accumulation for all time steps within one
c     day the "bar" variables represent 24-hour averages
      vTopB=vTopB + vTop*dt;
      vRootB=vRootB + vRoot*dt;
      rTopB=rTopB + rTop*dt;
      vBotB=vBotB + vBot*dt;
      if (vTop.gt.zero) then
        eTopB=eTopB + vTop*dt
      else
c       sum the wet-day total soil evaporation
        eTopB=eTopB + wDayEvap*dt
      endif
      return
      end
c
c     Profile of chemical concentrations for the three soluble
c     species and the immobile organic N content.
      subroutine chemOut(N,tPrint,x,Conc,CFix,ierr,errStr)
      implicit none
      include "defs.fi"
      double precision t,x,tPrint,Conc,CFix
      integer N,ierr
      character*80 errStr
      dimension x(N),Conc(NConcD,N),CFix(nFixC,N)
c
      integer i
      double precision MWSalt,MWN
      data MWSalt,MWN /5.6965D01,1.4D01/
c
      ierr=0
      errStr='ChemProf.dat'
      write(nChemUnit,100,err=901) tPrint
      write(nChemUnit,105,err=901) 
      do i=N-1,1,-1
        write(nChemUnit,110,err=901) i,x(i),Conc(1,i),
     ;     Conc(2,i)*MWSalt,Conc(3,i)*MWN,CFix(1,i)
      enddo
100   format('Time =',1X,F8.2,'(d)')
105   format(
     ;'Node      Depth        Tracer    Salt     NO3-N      Org_N'/
     ;'          (cm)         (-/L)    (mg/L)    (mg/L)   (g/g Soil)')
110   format(I4,2X,F8.1,1X,3(3X,F8.1),2X,E12.3)
      return
901   ierr=10
      return
      end
c
c     Chemical balance calculation and output
      subroutine chemBal(N,NMat,nConc,volConv,t,x,ThN,Conc,CFix,
     ;                   sumTopQ,sumBotQ,sumRootN,lastM,sConvN,
     ;                   chemPar,MatNum,initFlg,ierr,errStr)
      implicit none
      include "defs.fi"
      logical initFlg
      integer N,NMat,MatNum,nConc,ierr
      character*80 errStr
      double precision t,x,ThN,Conc,sumTopQ,sumBotQ,sumRootN,
     ;                 lastM,volConv,sConvN,cFix,chemPar
      dimension lastM(nConc),x(N),ThN(N),Conc(nConc,N),sumTopQ(nConc),
     ;          sumBotQ(nConc),sumRootN,cFix(nFixC,N),MatNum(N),
     ;          chemPar(nChemPar,NMat)
      double precision zero,two,dx,subOrgN
      data zero,two / 0.0D+00,2.0D+00 /
c
c     Local Arrays
c
      double precision  SubM,errM,relErr,MWN
      dimension         SubM(nConc),errM(nConc),relErr(nConc)
      logical            trdone
      double precision    delM,epserr,rhoBulk,convC,small,eps
      integer i,j
      data MWN,convC,small,eps /1.4D01,1.4D03,1.0D-15,1.0D-30/
c
      ierr=0
      errStr='ChemBal.dat'
      epserr=1.0D-06
c
c         j=1: Tracer in solution
c         j=2: Salt in solution
c         j=3: NO3-N in solution
c
      j=1
      do while (j.le.nConc)
        SubM(j)=zero
        do i=N-1,1,-1
          dx=(x(i+1)-x(i))
c
c         volConv converts concentration in mmol/L to mmol/cm^3
          SubM(j)=SubM(j)+volConv*dx*
     ;            (Conc(j,i+1)*ThN(i+1)+Conc(j,i)*ThN(i))/two

        enddo
        if (abs(SubM(j)).lt.small) then
          SubM(j)=zero
        endif
        j=j+1
      enddo

      if (initFlg) then
        do j=1,nConc
          lastM(j)=SubM(j)
        enddo
        sConvN=zero
        return
      endif
c
c     calculate mass of N as organic N (CFix: g/g soil)
      subOrgN=zero
      do i=N-1,1,-1
        dx=x(i+1)-x(i)
        rhoBulk=chemPar(2,MatNum(i))
        subOrgN = subOrgN + cFix(1,i)*dx*rhoBulk
      enddo
c     convert g/cm^2 to kg/ha
      subOrgN=subOrgN*1.0D05
c
c     calculate the error that has accumulated since the last
c     printout time. Take total column mass at last time, add
c     the cumulative flux in (a negative number), subtract 
c     cumulative flux out and,for NO3-N add the actual cumulative
c     root N uptake.Note when top flux is negative that corresponds
c     to adding salt to the column
      do j=1,nConc
        errM(j) = lastM(j) - sumTopQ(j) + sumBotQ(j)
        errM(j) = subM(j) - errM(j)
        if (j.eq.3) then
          errM(j) = errM(j) + sumRootN - sConvN
        endif
        if (abs(errM(j)).lt.small) then
          errM(j)=zero
        endif
        if (subM(j).gt.zero .and. errM(j).gt.zero) then
          relErr(j) = min(abs(errM(j)*1.0D02/subM(j)),1.0D02)
        else
          relErr(j)=zero
        endif
      enddo
c     Make sure vectors are printable by setting zeros for
c     small values
      do i=1,nConc
        if (abs(sumTopQ(i)).lt.eps) then
          sumTopQ(i)=zero
        endif
        if (abs(sumBotQ(i)).lt.eps) then
          sumBotQ(i)=zero
        endif
      enddo
c
c     convert to units of kg/ha
      write(nCBalUnit,110,err=901) t
      write(nCBalUnit,120,err=901) (convC*sumTopQ(i),i=1,nConc) 
      write(nCBalUnit,130,err=901) (convC*sumBotQ(i),i=1,nConc)
      write(nCBalUnit,140,err=901) (convC*subM(i),i=1,nConc)
      write(nCBalUnit,180,err=901) (convC*errM(i),i=1,nConc)
      write(nCBalUnit,190,err=901) (relErr(i),i=1,nConc)
      write(nCBalUnit,160,err=901) 
c
c     convert sumRootN (mmol/cm^2) to (kg/ha)
c     convert sConvN (mmol/cm^2) to (kg/ha
c     Avoid writing out large negative exponentssoil
      if (abs(subOrgN).lt.small) then
        subOrgN=zero
      endif
      if (abs(sumRootN).lt.small) then
        sumRootN=zero
      endif
      if (abs(sConvN).lt.small) then
        sConvN=zero
      endif
      write(nCBalUnit,170,err=901) subOrgN,sumRootN*MWN*1.D02,
     ;                             sConvN*MWN*1.D02
c
c         Put a lower limit for calculating RelErr, otherwise can get
c         meaningless division by very small numbers. Put 999 as flag.

      write(nCBalUnit,200,err=901)
c
c     Reinitialize all accumulators
      do j=1,nConc
        lastM(j)=subM(j)
        sumTopQ(j)=zero
        sumBotQ(j)=zero
      enddo
      sumRootN=zero
      sConvN=zero
      return
901   ierr=10
      return
110   format(///' Time  [T]:',F12.3/80('-'),/20X,'Tracer',5X,'Salt',
     ;       7X,5HNO3-N,/80('-'))
120   format(' Cum. top flux    ',3(E10.4,1X),'  kg/ha')
130   format(' Cum. lower flux  ',3(E10.4,1X),'    "  ')
140   format(' Mass in solution ',3(E10.4,1X),'    "  ')
150   format(' Mass of Organic N ',E9.2)
160   format(' Mass OrgN        NO3-N Uptake        OrgN->NO3-N'/
     ;       '  (kg/ha)           (kg/ha)             (kg/ha)')
170   format(1X,E10.3,7X,E10.3,11X,E10.3)
180   format(' Absolute error   ',3(E10.3,1X),'    "  ')
190   format(' Relative error[%]',3(F10.3,1X)/)
200   format(80('-'))
      end
c
c     Daily output of potential and actual N uptake and
c     leached N in kg N/ha.
***********************************************************************
      subroutine NPrnt(N,t,sMassN,sPotN,initNO3,hflag,sumDel,sumNPot,
     ;                 leachN,tLeachN,x,xRMax,Conc,theta,volConv,ierr,
     ;                 errStr)
      implicit none
      include "defs.fi"
      integer N,ierr
      logical initNO3,hflag
      character*80 errStr
      double precision t,dt,sMassN,sPotN,sumDel,x,xRMax,Conc,theta,
     ;                 volConv,sumNPot,leachN,tLeachN
      dimension x(N),Conc(NConcD,N),theta(N)
c
      integer i,count
      double precision availN,dx,depth,zero,two,eps,oLeachN,otLeachN,
     ;                 small
      data zero,two,eps,small /0.0D00,2.0D00,1.0D-15,1.0D-30/
c
      errStr='NO3.dat'
      if (hflag) then
        write(nNO3unit,120,ERR=901)
        return
      endif
      oLeachN=leachN
c     Units of tLeachN are kg N/ha
      otLeachN=tLeachN
      if (abs(oLeachN).lt.small) then
        oLeachN=zero
      endif
c     output total leached N is normally a positive number
      if (abs(otLeachN).lt.small) then
        otLeachN=zero
      endif
c     Calculate available N in the depth range from the surface
c     to the maximum rooting depth
      availN=zero
      depth=x(N)-x(N-1)
      i=N-1
      count=1
      do while (depth.le.xRMax)
        dx = x(i+1)-x(i)
        availN=availN + volConv*dx*1.4D01*
     ;      (Conc(3,i)*theta(i))
        i=i-1
        count=count+1
        depth = depth + dx
      enddo
c     convert availN: mmol*1e-03 = mol
c                     mol*MWN = g
c                     g * 1.e-03 = kg
c                     kg/cm^2 *1.0e08 = kg/ha
c     sPotN is maximum potential N uptake, convert uptake numbers
c     and available N from mg/cm^2 to kg/ha by multiplying by 1400
      if (sumDel.lt.eps) sumDel=zero
	  if (sumNPot.lt.eps) sumNPot=zero
      if (sumNPot.gt.zero) then
        write(nNO3unit,100,ERR=901) t,sPotN*1.4D03,sumDel*1.4D03,
     ;                              sumDel/sPotN,oleachN,
     ;                              otLeachN,availN*1.0D02
      else
        write(nNO3unit,100,ERR=901) t,sPotN*1.4D03,sumDel*1.4D03,
     ;                              zero,oleachN,otLeachN,
     ;                              availN*1.D02
      endif
      initNO3=.TRUE.
      return
901   ierr=10
      return
100   format(1X,F8.2,4(1X,E13.3),1X,E13.5,1X,E12.3)
120   format(3X,'Time',5X,'SumPotNup',5X,'SumActualNup',3X,
     ;       'Nuptake/PotN',2X,'Leach N',6X,'Tot Leach N',5X,
     ;       'Avail N',/,
     ;       2X,'(days)',7X,'(kg/ha)',7X,'(kg/ha)',9X,'(kg/ha/d)',
     ;       2X,'(kg/ha)',7X,'(kg/ha)',7X,'(kg/ha)')
      end
c
c     Daily output of variables related to stress effects of salinity
c     drought and N-deficit. sNoComp is the 24-hour average for alpha
c     calculated by the Pang-Letey stress function for salinity and 
c     matric stress. SumEff is compensated value and sumAlfG is 24-hour
c     average of the N-deficit stress. The final entry is the Rf
c     variable representing the stunting effects of all stresses as
c     the season progresses.
      subroutine alfaPrn(hFlag,t,indCrop,crpIndex,cropProp,sNoComp,
     ;                   sumEff,sumAlfG,counter,ierr)
      implicit none
      include "defs.fi"
      logical hFlag
      integer indCrop,crpIndex,counter,ierr
      double precision t,sNoComp,sumEff,sumAlfG,cropProp
      dimension cropProp(nCropPar,*)
c
      double precision zero
      data zero /0.0D00/
c
      if (hFlag) then
        write (nAlfUnit,100,ERR=901)
        write (nAlfUnit,105,ERR=901)
      else
c
        write(nAlfUnit,'(F6.1,$)') t
        write(nAlfUnit,'(I4,$)') indCrop
        write(nAlfUnit,'(4X,F8.3,$)') sNoComp/dble(counter)
        write(nAlfUnit,'(4X,F8.3,$)') sumEff/dble(counter)
        write(nAlfUnit,'(4X,F8.3,$)') sumAlfG/dble(counter)
        write(nAlfUnit,'(3X,F8.1)') cropProp(7,crpIndex)*1.0D02
      endif
      sumEff=zero
      sumAlfG=zero
      sNoComp=zero
      counter=0
      return
901   ierr=1
      return
c100   format(2X,F6.1,2X,I5,4(5X,E12.3))
100   format(
     ;' Time  CropID  alphaNoComp   alphaComp   alphaN    RelFact')
105   format(
     ;' (d)     -         -           -           -          %')
      end
c
c     Output occurs only at the end of a cropping season. Variables
c     are the relative yields calculated for salinity-matric stress
c     and N-deficiency. The cumulative root water uptake is sumRoot (cm)
c     and the cumulative leached N is sLeachN (kg N/ha).
      subroutine wrtCrop(hFlag,t,cropNum,cropName,relY,sumRoot,
     ;                   sPotN,sLeachN,prodNSea,sEPotMax,
     ;                   sumDel,lNitrate,ierr)
      implicit none
      include "defs.fi"
      logical hFlag,lNitrate
      integer cropNum,ierr
      character*20 cropName
      double precision t,relY,sumRoot,sumNPot,sPotN,sLeachN,
     ;                 prodNSea,sEPotMax,sumDel
c
      double precision zero,MWN
      data zero,MWN /0.0D00,1.4D01/
c
      if (hFlag) then
        if (lNitrate) then
          write(nRelUnit,100,ERR=901)
          write(nRelUnit,105,ERR=901)
        else
          write(nRelUnit,110,ERR=901)
          write(nRelUnit,115,ERR=901)
        endif
        return
      endif
      if (lNitrate) then
c       for N uptake convert mg/cm^2 to kg N/ha using factor 1400
        write(nRelUnit,120,ERR=901) t,cropNum,cropName,sumRoot,
     ;                       sPotN*1.4D03,sumDel*1.4D03,
     ;                       prodNSea*MWN*1.0D02,abs(sLeachN),
     ;                       min(sumRoot/sEPotMax*1.D02,1.0D02),
     ;                       min(1.D02*sumDel/sPotN,1.0D02)
c       set production of N sum to zero to start next crop
        prodNSea=zero
      else
        write(nRelUnit,125,ERR=901) t,cropNum,cropName,sumRoot,
     ;                              min(relY*1.D02,1.0D02)
      endif
      return
901   ierr=1
      return
100   format(
     ;'   Time   CropID  Name     SumRoot  PotNup  ActNup NProd'
     ;' NLeach RelYH2O RelNPotUptake')
105   format(
     ;'   (d)     -      -        (cm)    (kg/ha) (kg/ha) (kg/ha)',
     ;' (kg/ha)  (%)       (%)'/)
110   format(
     ;'  Time   CropID  Name       SumRoot   RelYH2O')
115   format(
     ;'   (d)     -      -         (cm)       (%)')
120   format(1X,F8.2,2X,I2,1X,A10,1X,4(F6.1,2X),F6.1,1X,F6.1,1X,
     ;       F6.1)
125   format(1X,F8.2,2X,I2,3X,A10,1X,F6.1,3X,F6.1)
      end
c
c     Daily output of variables related to organic N.
      subroutine orgNOut(N,NLayers,NMat,t,z,cFix,layerDep,chemPar,
     ;                   matNum,prevMass,theta,Conc,prntVec,ierr,
     ;                   errStr)
      implicit none
      include "defs.fi"
      integer N,NLayers,NMat,matNum,ierr
      double precision t,z,cFix,layerDep,chemPar,prevMass,prntVec,
     ;                 theta,Conc
      character*80 errStr
      dimension cFix(nFixC,N),z(N),layerDep(NLayers),
     ;          prntVec(NPrntD),chemPar(nChemPar,NMat),matNum(N),
     ;          prevMass(NLayers),theta(N),Conc(NConcD,N)
c
      double precision dz,rhoBulk,Mass,delMass,massNLay,zero
      dimension Mass(NLayers),delMass(NLayers),massNLay(NLayers)
      integer i,layNum,layOfNod,numPrnt
      external layOfNod
      data zero /0.0D00/
c
      ierr=0
      errStr='OrgNMass.dat'

      dz=z(N)-z(N-1)
      call layNMass(N,NLayers,NMat,matNum,dz,chemPar,CFix,
     ;              layerDep,Mass,ierr)
      do i=NLayers,1,-1
        delMass(i)=Mass(i)-prevMass(i) !rate of N-mineralization (g N/cm^2/d)
        delMass(i)=delMass(i)*1.D05    !convert to kg/ha/d
      enddo
c
c     Get the PAN mass for each layer
      call layMass(N,NLayers,z,theta,layerDep,Conc,
     ;             massNLay,ierr)
c     Don't write the record if an application of N has occurred (delMass>0)
      prntVec(1)=t
      prntVec(2)=Mass(NLayers)*1.D05 !convert to kg/ha
      prntVec(3)=delMass(NLayers)
	  prntVec(4)=massNLay(NLayers) !Daily mean mass PAN per layer
      NumPrnt=4
      if (NLayers.gt.1) then
        prntVec(5)=Mass(NLayers-1)*1.D05 !convert to kg/ha
        prntVec(6)=delMass(NLayers-1)
        prntVec(7)=massNLay(NLayers-1)
        NumPrnt=7
      endif
      if (delMass(NLayers).le.zero) then
        call orgNPrnt(.FALSE.,NumPrnt,prntVec,ierr)
      endif
c
c     End of day reset mass variables
      do i=NLayers,1,-1
        prevMass(i)=Mass(i)
      enddo
      do i=1,NPrntD
        prntVec(i)=zero
      enddo
      return
      end
c
c     Subsidiary subroutine for output of N variables by layer
      subroutine orgNPrnt(hFlag,Nprnt,prntVec,ierr)
      implicit none
      include "defs.fi"
      integer Nprnt,ierr
      logical hFlag
      double precision prntVec(NPrntD)
c
      integer i
c
      if (hFlag) then
        if (NPrnt.eq.4) then
c         Only one layer is present
          write(nNrtUnit,210)
        else
c         Several layers present, print top 2
          write(nNrtUnit,220)
        endif
        return
      endif
      if (Nprnt.eq.4) then
        write(nNrtUnit,110,ERR=900) (prntVec(i),i=1,Nprnt)
c     otherwise just write out the N status of the top two layers
      elseif (Nprnt.eq.7) then
        write(nNrtUnit,120,ERR=900) (prntVec(i),i=1,Nprnt)
      endif
      return
900   ierr=10
      return
110   format(2X,F10.3,3(3X,E12.3))
120   format(2X,F10.3,6(3X,E12.3))
210   format(
     ;'     Time            OrgN             MinRate      PANMass'/
     ;'      (d)           (kg N/ha)         (kg N/ha/d)  (kg N/ha)')
220   format(
     ;'     Time           OrgN(Top)      MinRate       PANMass',
     ;'      OrgN(Next)    MinRate       PANMass',/,
     ;'      (d)           (kg N/ha)      (kg N/ha/d)   (kg N/ha)',
     ;'       (kg N/ha)     (kg N/ha/d)   (kg N/ha)')
      end