*******************************************************************************
c Enviro-Gro Program is a subroutine of the console application
      subroutine E_G(cropSchd,cropProp)
      implicit none
c     Fixed parameter values
      include "defs.fi"
c     Declarations
      integer NumNP,NLayers,ldt,nSoilPrp,nTOut,Iter,ierr
c     Integer vectors
      integer layer(NLayerD),MatNum(NumNPD),LayMat(NLayerD)
c     Time variables
      logical dayFlag,pFlag
      integer ipFix,nTFix,itCum,itDay,nSteps
      double precision dt0,dt,dtMin,dtMax,dtOut,tInit,tFinal,t,
     ;                 tOld,dtOld,sumDt
c     Time vectors
      double precision tOutput(nTOutD),tFix(nTFixD)
c     Geometry variables
      double precision d,dz,depth
c     Geometry arrays
      double precision nodeZ(NumNPD),z(NumNPD),
     ;                 depthLay(NLayerD)
c     Water Flow vectors
      double precision hNew(NumNPD),hTemp(NumNPD),hOld(NumNPD),
     ;  Con(NumNPD),Cap(NumNPD),P(NumNPD),R(NumNPD),S(NumNPD),
     ;  ThNew(NumNPD),ThOld(NumNPD),ThStart(NumNPD),ConO(NumNPD),
     ;  ConR(NumNPD),AThS(NumNPD),ThRR(NumNPD),AKS(NumNPD),Sink(NumNPD),
     ;  Ah(NumNPD), AK(NumNPD), ATh(NumNPD),vOld(NumNPD),vNew(NumNPD),
     ;  Beta(NumNPD),Psi(NumNPD)
c     Soil properties
      integer nTab,nMat,MatID(NMatD)
      double precision thr(NMatD),ths(NMatD),hSat(NMatD),conSat(NMatD),
     ;                 alpha(NMatD),nvg(NMatD),Ks(NMatD),
     ;                 ParD(NTPar,NMatD)
      double precision hTab(NTabD)
      double precision conTab(NTabD,NMatD),capTab(NTabD,NMatD),
     ;                 Thetab(NTabD,NMatD),chemPar(nChemPar,NMatD)
c     Water flow variables
      double precision hTop,hBot,vTop,vBot,accVTop,accVBot,accSink,
     ;                 accETp,totInfil,hCritA,wDayEvap
      integer          watType
      double precision wTop,wBot,sumEPot,sumRoot,rel_fact
      double precision tPot,vRoot
      logical lFixFlag,lPlant,initFlg1,initFlg2,lTable,freeD,wLayer,
     ;        topInf,convgF
      double precision sumTPot,sumTPMax,sEPotMax
c     Boundary conditions
      double precision bndCond(nBndVars,nBndsD),bndTime(nBndsD),
     ;                 rTop,rBot,Transp,tNextBC
      integer          iBound,nBounds,kodTop,kodBot,ktOld,kbOld,
     ;                 chmBCTop,chmBCBot
c     Fertilizer additions
      double precision atmFall,mineral,amtAppl(nDisD),orgFrac(nDisD),
     ;                 minFrac(nDisD),tApplic(nDisD),tNextEvt,
     ;                 TVarCoef(nMonths),prevMass(NLayerD),
     ;                 PsiApp(nDisD)
      integer nDiscret,iApplic
      logical layerMix(NLayerD),lTCoef
c
c     Concentrations
      integer          maxCh,nConc,solNum(MaxIrriD),nSols
      logical          lNitrate
      double precision traceLay(NLayerD),saltLay(NLayerD),
     ;                 NMassFix(NLayerD),NMassMin(NLayerD),cTop(nIrriC),
     ;                 cBot(nIrriC),CFix(NFixC,NumNPD),DOld(NumNPD),
     ;                 DNew(NumNPD),Conc(nConcD,NumNPD),
     ;                 RConc(nConcD,NumNPD),TConc(nConcD,NumNPD),
     ;                 irriChem(nIrriC*2,MaxIrriD)
      double precision orgNFrac(NLayerD),PsiZero(NLayerD)
c     Transportable species are:
c                      Tracer - Conc(1,:)
c                      Salt -   Conc(2,:)
c                      Nitrate  Conc(3,:)
c     Vectors allocated for use by subroutines
      double precision Bvec(NumNPD),Dvec(NumNPD),Evec(NumNPD),
     ;                 Fvec(NumNPD),Tracer(NumNPD)
c     Chemical transport and N uptake
      integer          compCnt
      double precision Peclet,Courant,cvTop(nIrriC),cvBot(nIrriC),
     ;                 dtMaxC,sumTopQ(nIrriC),sumBotQ(nIrriC),sumRootN,
     ;                 lastM(nIrriC),tMassR,sMassN,sPotN,tPotN,mStart,
     ;                 alphaG,tPotMax,tPotNMax,sConvN,prodNSea,tLeachN,
     ;                 cStarN,sumNAct,sumNPot,leachN,sLeachN,sumDel
c
c     Chemical variables
      double precision NLossRat
c
c     Crop Schedule
      character*20 scdNam(maxScd),dateInit,dateFinal
      double precision cropSchd(nSchdPar,maxScd)
      integer nCrpRecs,cropNum
      character*80 errStr
      logical stSeason,cropDone
      data stSeason,cropDone /.True.,.False. /
c     Crop properties
      double precision cropProp(nCropPar,maxCrop)
      integer cropID(maxCrop),iCrop,currCrop,nPlants,crpIndex,
     ;        getCIndx
      character*20 cropName
c
c     Root Growth
      integer iRoot
      double precision zR,maxRoot
c
c     Stress functions
      double precision alphaHS(NumNPD),alphaHSR(NumNPD)
      double precision sNoComp,alfaEff,sumAlfG,sumEff,alfaGHS,
     ;                 Et_shres
c
c     Output
      character*149 fname
      integer strlen,lenPath,cntCall
      external strlen
      logical lPrint,lPrnFlag,initNO3
      integer iPrint,lineCnt,cntBare,prevCrop,i
      double precision vTopB,vRootB,rTopB,vBotB,startT,eTopB
      double precision sLayTopQ,sLayBotQ,sLaySink,prntVec,
     ;                 dayNSum
      dimension sLayTopQ(NLayerD),sLayBotQ(NLayerD),
     ;          sLaySink(NLayerD),prntVec(NPrntD),dayNSum(NLayerD)
c     External functions
      double precision FH,NTReduc
      external FH,getCIndx,NTReduc
c     Constants
      double precision zero,volConv
      data zero,volConv /0.0D00, 1.0D-03/
c     Initialization
      data nConc,vTopB,vRootB,rTopB, vBotB
     ;     /3,  0.0D00,0.0D00,0.0D00,0.0D00/
      data initNO3,pFlag /.TRUE.,.TRUE. /
c
c Data directory path is stored in dataPath.txt, open this file
c which must be in the same directory as the data, and read
c in the path to the directory of input and output files.
c If this file is not present further execution of the program
c is halted.
      character*120 datPath
      open(UNIT=10,FILE='dataPath.txt',STATUS='OLD',ERR=900)
      read(UNIT=10,*,ERR=900) datPath
      close(UNIT=10)
c
c Open error log file and write header
      lenPath=strlen(datPath)
      fname = datPath(1:lenPath)//'\ErrLog.dat'
      open(UNIT=nErrLog,FILE=fname,ERR=901)
      write(nErrLog,1010,ERR=901)
c
      call setFlags(lFixFlag,lPrnFlag)
c Read input data
      call rdTime(datPath,ntOut,dt0,dt,dtMin,dtMax,dtOut,dtOld,
     ;            dateInit,dateFinal,tOutput,tInit,tFinal,ierr,
     ;            errStr)
      if (ierr.gt.0) call errHand(nErrLog,ierr,errStr)
      call rdBasic(datPath,NumNP,NLayers,nTab,d,dz,
     ;             nodeZ,z,depthLay,layer,layMat,lNitrate,
     ;             NLossRat,TVarCoef,lTCoef,freeD,lPlant,
     ;             hCritA,ierr,errStr)
      if (ierr.gt.0) call errHand(nErrLog,ierr,errStr)
c
c     Open the output files and write headers
      call outOpen(datPath,lNitrate,lPlant,ierr)
      if (ierr.gt.0) call errHand(nErrLog,ierr,errStr)
      call headers(NumNP,cropProp,NLayers,lNitrate,lPlant,z,
     ;             Conc,ThNew,ierr)
c     write watBal header if there is only one layer
      if (NLayers.eq.1) then
        call watBal(0,z,zero,thStart,thNew,zero,zero,zero,
     ;              zero,.TRUE.,ierr,errStr)
      endif
c
c     read soil properties
      call rdSoil(datPath,nMat,MatID,depth,alpha,Ks,
     ;            nvg,thR,thS,ierr,errStr)
      if (ierr.gt.0) call errHand(nErrLog,ierr,errStr)
c
c     read soil chemical parameters
      call rdChmPar(datPath,NMat,MatID,chemPar,ierr,errStr)
      if (ierr.gt.0) call errHand(nErrLog,ierr,errStr)
c
c     using the vectors of each soil parameter build the array
c     of soil hydraulic properties (ParD)
      call genMat(nMat,MatID,alpha,Ks,nvg,thR,thS,ParD,ierr,
     ;            errStr)
      if (ierr.gt.0) call errHand(nErrLog,ierr,errStr)
c
c     read initialization data for each layer
      call rdInit(datPath,NumNP,NLayers,wTop,wBot,watType,
     ;            layer,saltLay,traceLay,NMassMin,NMassFix,
     ;            orgNFrac,PsiZero,layerMix,ierr,errStr)
      if (ierr.gt.0) call errHand(nErrLog,ierr,errStr)
c
c     Nodal concentrations in Conc are transportable species, in
c     CFix species are not transportable they are only reactive
      call nodeBld(NumNP,NLayers,NMat,depthLay,layMat,MatNum,dz,
     ;             saltLay,traceLay,NMassFix,Conc,CFix,orgNFrac,
     ;             Psi,PsiZero,chemPar,prevMass,ierr)
      if (ierr.gt.0) call errHand(nErrLog,ierr,errStr)
c
c     These boolean variables are all assigned true pending further
c     potential developments of the program.
      call setBool(wLayer,topInf,lTable)
c
c     read the boundary conditions for the entire run
      call rdBounds(datPath,dateInit,tInit,bndCond,nBounds,bndTime,
     ;              freeD,kodTop,kodBot,lNitrate,ierr,errStr)
      if (ierr.gt.0) call errHand(nErrLog,ierr,errStr)
      call setHBnd(NLayers,NMat,IWMODEL,watType,freeD,wTop,wBot,
     ;             ParD,layMat,bndCond,hBot,hTop)
c
c     Set the initial heads by linear interpolation
      call nodeHead(NumNP,hOld,hTemp,hNew,hTop,hBot)
c
      call setModif(NumNP,Ah,AK,ATh,AThS,AKS)
c     Use the simplified V-G functions for calculation of water retention
c     IWMODEL=1
      call setTab(NTab,NMat,NumNP,hTab,IWMODEL,ParD,hSat,thr,ths,conSat,
     ;            conTab,capTab,theTab,ATh)
c
      Iter=0
c
c     Initial calculation of hydraulic conductivity, capacitance and
c     theta for all nodes based on pressure head (hNew)
      call setMat(NumNP,NTab,NMat,Iter,hTab,ConTab,CapTab,hNew,MatNum,
     ;            ParD,Con,Cap,ConSat,Ah,AK,ATh,hSat,hTemp,TheTab,
     ;            ThOld,thr,ths,ConO,AThS,ThRR,ConR,AKS,IWMODEL,
     ;            lTable)
c
c     Convert initial dissolved N specified in kg/ha to Conc (mmol/L)
      call layConc(NumNP,NLayers,NMat,z,ThOld,depthLay,NMassMin,Conc,
     ;             cFix,chemPar,MatNum,ierr)
c
c     vector copy for theta variables
      call veccpr(1,NumNP,ThOld,ThNew,ierr)
      call veccpr(1,NumNP,ThOld,ThStart,ierr)
c
c     intitialize accumulators for water balance between fixed
c     output times
      call initBal(accVTop,accVBot,accSink,totInfil)
      if (ierr.gt.0) call errHand(nErrLog,ierr,errStr)
c
c     Get the irrigation water chemical composition
      call rdChem(datPath,nSols,irriChem,solNum,ierr,errStr)
      if (ierr.gt.0) call errHand(nErrLog,ierr,errStr)
c
c     Read the plant data including the plant-specific properties
c     and the crop schedule
      if (lPlant) then
        call rdPlant(datPath,nPlants,cropProp,cropID,maxRoot,ierr,
     ;               errStr)
        if (ierr.gt.0) call errHand(nErrLog,ierr,errStr)
c
c     Get the crop schedule and set tInit as DOY of the starting year,
c     tFinal as tInit plus the elapsed days between dateInit and dateFinaL
        call rdSched(datPath,nCrpRecs,scdNam,cropSchd,dateInit,
     ;               dateFinal,tInit,tFinal,cropNum,ierr,errStr)
        if (ierr.gt.0) call errHand(nErrLog,ierr,errStr)
      endif
c
c     Read fertilization event data for organic-N and/or conventional
c     nitrate.
      if (lNitrate) then
        call rdAdd(datPath,dateInit,tInit,nDiscret,atmFall,
     ;           mineral,amtAppl,orgFrac,PsiApp,tApplic,tNextEvt,
     ;           iApplic,ierr,errStr)

        if (ierr.gt.0) call errHand(nErrLog,ierr,errStr)
        call initNSum(NLayers,dayNSum,tLeachN)
      else
        nDiscret=0
      endif
c
c     Build the tFix array by merging tOutput,bndTime and tApplic vectors
      call bldTFix(nBounds,nTout,nDiscret,bndTime,tOutput,tApplic,
     ;             tFix,tFinal,nTFix,ipFix,iPrint,ierr)
      if (ierr.gt.0) call errHand(nErrLog,ierr,errStr)
c
c     Initialize variables prior to main loop
      t = tInit
      tOld = tInit
      sumDt = zero
      nSteps = 0
      itDay = int(t)
      lineCnt=20
      iCrop=1
      initFlg1=.true.
      initFlg2=.true.
      compCnt=0
      cntCall=0
      if (lPlant) then
        crpIndex=getCIndx(nPlants,cropNum,cropProp,ierr,errStr)
        if (ierr.gt.0) call errHand(nErrLog,ierr,errStr)
c
c     Set root growth and Beta for the first scheduled crop
        call SetRG(NumNP,z,Beta,zR,t,dz,cropProp,nPlants,crpIndex,
     ;           cropSchd,cropID,nCrpRecs,iCrop,iRoot,ierr)
        if (ierr.gt.0) call errHand(nErrLog,ierr,errStr)
        cropNum=INT(cropSchd(1,iCrop))
      else
        crpIndex=1 !need to initialize anyway
        cropNum=0
      endif
      call initSums(NumNP,sumEPot,sumRoot,cropProp,iCrop,
     ;              alphaG,sumEff,sumAlfG,accETp,sNoComp,
     ;              sumNAct,sumNPot,sPotN,stSeason,sEPotMax,
     ;              sumDel)

      if (NLayers.ne.1) then
        call initLay(NLayers,sLayTopQ,sLayBotQ,sLaySink)
      endif
c
c     Initialize chemical balance calculation
      call chemBal(NumNP,NMat,nConcD,volConv,t,z,ThNew,Conc,CFix,
     ;             sumTopQ,sumBotQ,sumRootN,lastM,sConvN,
     ;             chemPar,MatNum,.TRUE.,ierr,errStr)
      if (ierr.gt.0) call errHand(nErrLog,ierr,errStr)
c
c     Boundary condition counter
      iBound=1
      call setBC(nBounds,iBound,nSols,bndCond,irriChem,t,tFinal,dtMax,
     ;           bndTime,tNextBC,rTop,tPotMax,tPotN,tPotNMax,KodTop,
     ;           KodBot,kTOld,kBOld,chmBCTop,chmBCBot,rBot,cTop,cBot,
     ;           cropNum,cStarN,hBot,wDayEvap,lNitrate,freeD,.TRUE.,
     ;           ierr)
c
c     Bare soil evap period counter
      cntBare=0
      prodNSea=zero
c
c     Special case where the crop starts growing immediately at t=tInit
      if (lPlant .and. cropNum.gt.0) then
        call bcAccom(bndCond,iBound,cropProp,cropNum,rTop,tPotMax,tPotN,
     ;               tPotNMax,tpot,nPlants,ierr,errStr)
      endif

c     Start main time-based loop
      do while (t.le.tFinal .and. .not.cropDone)
c
      sumDt=sumDt + dt
      nSteps=nSteps + 1
      if (lPlant) then
        currCrop=iCrop
        call advCrop(cropSchd,cropNum,iCrop,scdNam,cropName,nPlants,
     ;               cropProp,crpIndex,t,sumRoot/accETp,sumRoot,
     ;               sumDel,sPotN,sLeachN,prodNSea,Et_shres,initFlg1,
     ;               nCrpRecs,cropDone,lNitrate,prevCrop,sEPotMax,
     ;               ierr,errStr)
c
c       Reset the sums of potential ET and root water uptake for new crop
c       and apply boundary conditions for new crop
        if (iCrop.gt.currCrop) then
          call initSums(NumNP,sumEPot,sumRoot,cropProp,crpIndex,
     ;                  alphaG,sumEff,sumAlfG,accETp,sNoComp,sumNAct,
     ;                  sumNPot,sPotN,stSeason,sEPotMax,sumDel)
          call initBal(accVTop,accVBot,accSink,totInfil)
c
          call bcAccom(bndCond,iBound,cropProp,cropNum,rTop,tPotMax,
     ;                 tPotN,tPotNMax,tPot,nPlants,ierr,errStr)
c         calculate root growth if crop present
          if (cropNum.gt.0) then
            call setRG(NumNP,z,Beta,zR,t,dz,cropProp,nPlants,crpIndex,
     ;           cropSchd,cropID,nCrpRecs,iCrop,iRoot,ierr)
          else
            vRoot=zero
          endif
c
c         Set seasonal N production to zero
          if (prevCrop.eq.0) then
            prodNSea=zero
          endif
        endif
      endif
c
c     Calculate new pressure heads and water contents
      call watFlow(NumNP,NTab,NMat,hTab,ConTab,CapTab,hNew,
     ;             hOld,MatNum,ParD,Con,Cap,ConSat,Ah,AK,ATh,
     ;             hSat,hTemp,KodTop,KodBot,rTop,rBot,dt,
     ;             z,Sink,P,R,S,FreeD,hTop,hBot,WLayer,Iter,
     ;             ItCum,TopInf,KTOld,KBOld,dtMin,t,tOld,
     ;             ConvgF,TheTab,ThNew,ThOld,thr,ths,vTop,
     ;             AThS,ThRR,ConO,ConR,AKS,hCritA,IWMODEL,
     ;             lTable,ierr)
      if (ierr.gt.0) then
        call errHand(nErrLog,ierr,errStr)
      elseif (Iter.eq.20) then
c       If number of iterations is 20 then there was failure to
c       converge, write this to the error log file
        write(nErrLog,970) t,hNew(NumNP),hNew(NumNP-1),
     ;                     hNew(NumNP-2)
      endif
c
c Recalculate the velocity as a function of depth
      call W_Flux(NumNP,dt,vTop,vBot,vNew,ThNew,ThOld,z,Sink)
c
c Calculate new solute concentrations
      call TrnsPort(NumNP,NMat,volConv,t,dt,z,MatNum,ChemPar,ThOld,
     ;              ThNew,vOld,vNew,Dold,DNew,cTop,cBot,iRoot,
     ;              chmBCTop,chmBCBot,Conc,TConc,CFix,NLossRat,
     ;              Bvec,Dvec,Evec,Fvec,cvTop,cvBot,Peclet,
     ;              Courant,dtMaxC,RConc,Sink,Psi,MaxCh,tPotMax,
     ;              tPotN,tPotNMax,tMassR,alphaG,Beta,sConvN,
     ;              prodNSea,TVarCoef,lTCoef,dateInit,cStarN,
     ;              sumDel,cropProp,crpIndex,nConc,cropNum,
     ;              compCnt,lNitrate,ierr)
      if (ierr.gt.0) call errHand(nErrLog,ierr,errStr)
c          
c     Plant root water uptake
      if (lPlant .and. cropNum.gt.0) then
c       Root growth and root density function (Beta) update
        call setRG(NumNP,z,Beta,zR,t,dz,cropProp,nPlants,crpIndex,
     ;             cropSchd,cropID,nCrpRecs,iCrop,iRoot,ierr)
        if (ierr.gt.0) call errHand(nErrLog,ierr,errStr)
c       Root water uptake, update of Sink vector
        call sinker(NumNP,NMat,iRoot,crpIndex,dt,dz,hNew,Beta,
     ;              Conc(2,:),ThNew,cropProp,tPotMax,MatNum,ParD,
     ;              vRoot,Et_shres,alphaG,alphaHS,alphaHSR,
     ;              alfaGHS,alfaEff,Sink)
      endif
c
c     Accumulate water variables for daily and seasonal output
      call W_Cumul(NumNP,NLayers,cropNum,dt,dz,vTop,vBot,vNew,
     ;             vRoot,accVTop,totInfil,accVBot,accSink,alfaEff,
     ;             alphaG,sumEff,sumAlfG,tPotMax,accETp,sNoComp,
     ;             alfaGHS,sLayTopQ,sLayBotQ,sLaySink,nodeZ,Beta,
     ;             cropProp,crpIndex,Sink,depthLay,cntCall,
     ;             lPlant,ierr)
      if (ierr.gt.0) call errHand(nErrLog,ierr,errStr)
c     Chemical variables
      call C_Cumul(NumNP,NLayers,nConcD,z,volConv,dt,vTop,vBot,
     ;             cTop,Conc,thNew,sumTopQ,sumBotQ,tMassR,tPotN,
     ;             tPotNMax,sumRootN,sMassN,sPotN,initFlg2,
     ;             initNO3,depthLay,cvBot,leachN,sLeachN,tLeachN,
     ;             stSeason,freeD,lPlant,lNitrate)

c     Set flag for initializing sumBotQ and sumTopQ to false, reset
c     the relative growth factor variable stored in cropProp
      initFlg2=.FALSE.
      if (lPlant .and. cropNum.gt.0) then
        call setRel(lNitrate,cropProp,dt,zR,TPotMax,tPotN,tPotNMax,
     ;              sumNAct,sumNPot,tMassR,vRoot,sumEPot,sumRoot,
     ;              crpIndex,sEPotMax,ierr)
        if (ierr.gt.0) call errHand(nErrLog,ierr,errStr)
      endif
c
c     Printout various profiles at tOutput times
      if (lPrint) then
        call fixPrnt(NumNP,t,tOutput(iPrint),z,MatNum,hNew,
     ;               ThNew,Con,Sink,vNew,ierr,errStr)
        if (ierr.gt.0) call errHand(nErrLog,ierr,errStr)

        call chemOut(NumNP,tOutput(iPrint),z,Conc,CFix,ierr,errStr)
        if (ierr.gt.0) call errHand(nErrLog,ierr,errStr)
c
c       Use the simple watBal function for a single layer model
        if (NLayers.eq.1) then
          call watBal(NumNP,z,t,thStart,thNew,accVTop,accVBot,accSink,
     ;                totInfil,.FALSE.,ierr,errStr)
          if (ierr.gt.0) call errHand(nErrLog,ierr,errStr)

          call veccpr(1,NumNP,ThNew,ThStart,ierr)
        else
c         Printout the water balance for each layer in the model
          call layerBal(t,NumNP,NLayers,z,depthLay,thStart,thNew,
     ;                  sLayTopQ,sLayBotQ,sLaySink,ierr,errStr)
          if (ierr.gt.0) call errHand(nErrLog,ierr,errStr)
c
c         reinitialize the layer sums and thStart
          call initLay(NLayers,sLayTopQ,sLayBotQ,sLaySink)
          call veccpr(1,NumNP,ThNew,thStart,ierr)
        endif

        call chemBal(NumNP,NMat,nConcD,volConv,t,z,ThNew,Conc,CFix,
     ;               sumTopQ,sumBotQ,sumRootN,lastM,sConvN,
     ;               chemPar,MatNum,.FALSE.,ierr,errStr)
        if (ierr.gt.0) call errHand(nErrLog,ierr,errStr)

        if (lPlant .and. cropNum.gt.0) then
          call strsPrnt(NumNP,t,tOutput(iPrint),z,Sink,AlphaHS,
     ;                AlphaHSR,cropNum,cropName,iCrop,ierr,errStr)
          if (ierr.gt.0) call errHand(nErrLog,ierr,errStr)
        endif
        lPrint = .FALSE.
        iPrint = iPrint + 1
      endif

      call updOut(dt,vTop,vBot,vRoot,rTop,vTopB,vRootB,rTopB,eTopB,
     ;            vBotB,wDayEvap)
c
c     Daily printout to screen and files
      if (dayFlag) then
        startT=dble(int(t))
        call sreenOut(NumNP,t,rTop,vTop,vRoot,hTop,hBot,lineCnt,hNew,
     ;                KodBot,KodTop)
        call dayPrnt(t,vTopB,vRootB,rTopB,vBotB,eTopB,sumDt,nSteps,
     ;               .FALSE.,ierr,errStr)
        if (ierr.gt.0) call errHand(nErrLog,ierr,errStr)

        if (lNitrate) then
c         Set Rel_fact if N-stress is greatest
          call NPrnt(NumNP,t,sMassN,sPotN,initNO3,.FALSE.,sumDel,
     ;               sumNPot,leachN,tLeachN,z,maxRoot,Conc,ThNew,
     ;               volConv,ierr,errStr)
          if (ierr.gt.0) call errHand(nErrLog,ierr,errStr)

          call orgNOut(NumNP,NLayers,NMat,t,z,cFix,depthLay,
     ;                 chemPar,matNum,prevMass,ThNew,Conc,prntVec,
     ;                 ierr,errStr)
          if (ierr.gt.0) call errHand(nErrLog,ierr,errStr)

        endif
        if (lPlant .and. cropNum.gt.0) then
          call alfaPrn(.FALSE.,t,cropNum,crpIndex,cropProp,sNoComp,
     ;                 sumEff,sumAlfG,cntCall,ierr)
          call initDSum(sumEff,sumAlfG,sNoComp)
        endif
        dayFlag=.FALSE.
      endif !dayFlag
c
c Calculate dt and the time for the next time step
      call setTime(NTFix,ipFix,t,tOld,tFix,dt,dtMin,dtMax,dtOld,
     ;             Iter,tOutput,nTOut,iPrint,itDay,lPrint,lFixFlag,
     ;             lPrnFlag,dayFlag,ierr,errStr)
      if (ierr.gt.0) call errHand(nErrLog,ierr,errStr)
c
c Check if the time exceeds the next boundary condition time, if so
c increment the boundary condition counter and get new conditions
      if (t.ge.tNextBC) then
        iBound = iBound + 1
        call setBC(nBounds,iBound,nSols,bndCond,irriChem,t,tFinal,dtMax,
     ;             bndTime,tNextBC,rTop,tPotMax,tPotN,tPotNMax,KodTop,
     ;             KodBot,kTOld,kBOld,chmBCTop,chmBCBot,rBot,cTop,
     ;             cBot,cropNum,cStarN,hBot,wDayEvap,lNitrate,freeD,
     ;             .FALSE.,ierr)
        if (ierr.gt.0) call errHand(nErrLog,ierr,errStr)
      endif
c
c Check for time exceeding next fertilizing event, if so rebuild the
c fertilizer profile in layers
      if (t.ge.tNextEvt .and. lNitrate) then
        call addFert(NumNp,NLayers,nDiscret,NMat,z,iApplic,t,tApplic,
     ;               tNextEvt,dtMax,amtAppl,minFrac,orgFrac,CFix,
     ;               layerMix,tFinal,depthLay,MatNum,chemPar,Conc,
     ;               ThNew,Psi,PsiApp,volConv,mStart,ierr)
        iApplic=iApplic + 1
      endif
c
c     Update water flow vectors
      call copyWat(NumNP,hOld,hNew,hTemp,thOld,thNew,vOld,vNew,
     ;            dt,dtOld)
c
c End of the main loop
      end do
c If the main loop ended without writing the final crop record write it now
      if (.not.cropDone .and. nCrpRecs.gt.0 .and. cropNum.gt.0) then
        call wrtCrop(.FALSE.,cropSchd(3,nCrpRecs),
     ;               int(cropSchd(1,nCrpRecs)),scdNam(nCrpRecs),
     ;               sumRoot/AccETp,sumRoot,sPotN,
     ;               sLeachN,prodNSea,sEPotMax,sumDel,lNitrate,ierr)
      endif
c
c Close output files
      call closer(nProfUnit)
      call closer(nWBalUnit)
      call closer(nDayUnit)
      call closer(nChemUnit)
      call closer(nCBalUnit)
      call closer(nNO3Unit)
      return
900   continue
      write(*,*) 
     ;'dataPath.txt file must be in same directory as input files'
      stop
901   write(*,*) 'File system problem: Unable to open Error Log file'
      stop
970   format('Iter=20, Time, Pressure head top 3 nodes','/',
     ;'      t=',F12.5,3X,3(E15.6,3X))
1010  format('Error Log File')
      end