************************************************************************
c     Subroutines in this file are responsible for reading 
c     input files and generating vectors and arrays needed
c     before starting the main loop.
c
c     Initialization of various vectors from supplied data for 
c     initial conditions
      subroutine rdInit(datPath,NumNP,NLayers,wTop,wBot,watType,
     ;                  layer,saltLay,traceLay,NMassMin,NMassFix,
     ;                  orgNFrac,PsiZero,layerMix,ierr,errStr)
      implicit none
      include 'defs.fi'
      character*120 datPath
      integer NumNP,NLayers,watType,layer(NLayers),ierr
      logical layerMix(NLayers)
      double precision wTop,wBot
      double precision hTop,hBot
      double precision traceLay,NMassMin,NMassFix,orgNFrac,saltLay,
     ;                 PsiZero
      dimension saltLay(NLayers),NMassFix(NLayers),traceLay(NLayers),
     ;          orgNFrac(NLayers),PsiZero(NLayers),NMassMin(NLayers)
      character*80 errStr
c
      integer i,lenPath,ldt,n
      character*20 tokens(10),str
      character*120 line
      logical feof
      integer strlen,getTok
      external strlen,getTok
      double precision zero,one
      character*140 fname
      data zero,one /0.0D00, 1.0D00/
c
      ierr=0
      errStr='Init.in'
      lenPath=strlen(datPath)
      fname = datPath(1:lenPath)//'\Init.in'
      open(UNIT=nInitUnit,FILE=fname,ERR=900)
c
      call advance(nInitUnit,line,feof)
      ldt=3
      n = getTok(line,tokens,ldt)
      if (n.ne.ldt) then
        goto 901
      endif
c     watType determines whether data entered are for
c     pressure head - watType = 1 or
c     water content - watType = 2
      read(tokens(1),'(I5)') watType
      if (watType .lt. 1 .or. watType .gt. 2) then
        ierr = 84
        return
      endif
      read(tokens(2),'(F10.4)') wTop
      read(tokens(3),'(F10.4)') wBot
      call advance(nInitUnit,line,feof)
      ldt=7
      n = getTok(line,tokens,ldt)
      if (n.ne.ldt) then
        goto 901
      endif
      do i=NLayers,1,-1
        read(tokens(1),'(I5)') layer(i)
        read(tokens(2),*) traceLay(i)
        read(tokens(3),*) saltLay(i)  ! mg/L
c       NMassFix includes both organic and inorganic N
        read(tokens(4),*) NMassFix(i) !kg N/ha
c       fraction of NMassFix that is organic
        read(tokens(5),*) orgNFrac(i)
c       mass of mineral N applied
        NMassMin(i)=(one-orgNFrac(i))*NMassFix(i)
c       PsiZero is the fraction of organic N in applied fertilizer
c       that is in the slow pool when the fertilizer is applied. 
c       If the applied fertilizer is mixed with existing organic-N
c       the value of Psi will change.
        read(tokens(6),*) PsiZero(i)
        read(tokens(7),'(A)') str
c       the layerMix variable controls which layers will have
c       mixing of the existing and applied organic N when the 
c       organic N is applied
        if (str(1:1).eq.'T'.or.str(1:1).eq.'t') then
          layerMix(i)=.TRUE.
        else
          layerMix(i)=.FALSE.
        endif
        if (i.gt.1) then
          call advance(nInitUnit,line,feof)
          if (feof) goto 903
          n = getTok(line,tokens,ldt)
          if (n.ne.ldt) then
            goto 901
          endif
        endif
      enddo
      close(nInitUnit)
      return
c     File Init.in could not be opened
900   call closer(nInitUnit)
      ierr=1
      return
c     Incorrect number of entries in file
901   ierr=2
      call closer(nInitUnit)
      return
902   ierr=3
      call closer(nInitUnit)
      return
903   ierr=4
      call closer(nInitUnit)
      return
      end
c
c     Read the plant properties from Plant.in
      subroutine rdPlant(datPath,nPlants,cropProp,cropID,maxRoot,
     ;                   ierr,errStr)
      implicit none
      include 'defs.fi'
      character*120 datPath
      integer nPlants,cropID(nPlantD),ierr
      double precision cropProp(nCropPar,*),maxRoot
      character*80 errStr
c 
c   cropProp(1,i) = cropID
c   cropProp(2,i) = xRMax     maximum rooting depth
c   cropProp(3,i) = hShres    threshold pressure head for water uptake
c   cropProp(4,i) = hh5       pressure head for 50% reduction 
c   cropProp(5,i) = sShres    M-H threshold salinity
c   cropProp(6,i) = gg5       osmotic head for 50% reduction
c   cropProp(7,i) = rel_fact  plant growth reduction, entered elsewhere
c   cropProp(8,i) = hPWP      pressure head at permanent wilting
c
c Local variables and arrays
      integer i,lenPath,n,ldt
      double precision MHSlope
      character*20 tokens(11)
      character*140 line,fname
      logical feof
	  double precision zero,one
      data zero,one /0.0D00,1.0D00/
      integer strlen,getTok
      external strlen,getTok
c
      errStr='Plant.in'
      lenPath=strlen(datPath)
c
c Build the first record for the bare soil (fallow) condition
      do i=1,8
        cropProp(i,1)=zero
      enddo
c
c Build the first record for the fallow condition, cropNum=1
      cropProp(1,2)=one
      cropProp(2,2)=5.0D00   !xRMax
      cropProp(3,2)=zero     !hShres
      cropProp(8,2)=-1.5D04  !hPWP
      cropProp(5,2)=2.0D00   !M-H Threshold
      MHSlope=5.0D00   !M-H Slope
	  call vgParams(2,cropProp,MHSlope,ierr)
c
c Add records from data files
      fname = datPath(1:lenPath)//'\Plant.in'
      open(UNIT=nPlntUnit,FILE=fname,ERR=900)
c
      call advance(nPlntUnit,line,feof)
      ldt=1
      n = getTok(line,tokens,ldt)
      if (n.ne.ldt) then
        goto 901
      endif
      read(tokens(1),'(I5)') nPlants
      call advance(nPlntUnit,line,feof)
      ldt=6
      n = getTok(line,tokens,ldt)
      if (n.ne.ldt) then
        goto 901
      endif
      maxRoot=cropProp(2,2)
c     The first two cropID's are reserved so normal crop ID
c     entries start at index 3 in the vector
      do i=3,nPlants+2
        read(tokens(1),'(I5)') cropID(i)
        cropProp(1,i)=dble(cropID(i))
        read(tokens(2),*) cropProp(2,i) !xRMax
        read(tokens(3),*) cropProp(3,i) !hShRes
        read(tokens(4),*) cropProp(8,i) !hPWP
        read(tokens(5),*) cropProp(5,i) !M-H Threshold
        read(tokens(6),*) MHSlope       !needed for gg5
c       get hh5 and gg5
        call vgParams(i,cropProp,MHSlope,ierr)
        if (i.lt.(nPlants+2)) then
          call advance(nPlntUnit,line,feof)
          if (feof) goto 910
          n = getTok(line,tokens,ldt)
          if (n.ne.ldt) then
            goto 901
          endif
        endif
        maxRoot=max(maxRoot,cropProp(2,i))
      enddo
      close(nPlntUnit)
      return
c     File Plant.in could not be opened
900   ierr=1
      return
c     Incorrect number of entries in file
901   ierr=2
      return
c     Incorrect number of lines in file
910   ierr=4
      return
911   ierr=5
      return
      end

c     Calculate the h50 and pi50 values for use in the
c     van Genuchten stress function.
      subroutine vgParams(iRec,cropProp,MHSlope,ierr)
      implicit none
      include "defs.fi"
      integer iRec,ierr
      double precision cropProp(nCropPar,*),MHSlope
c
      double precision zero,one,two
      data zero,one,two /0.0D00,1.0D00,2.0D00/
c     Calculate the hh50 number from the threshold h and PWP h
      if (cropProp(3,iRec).gt.zero .or. 
     ;    cropProp(8,iRec).gt.cropProp(3,iRec)) then
        ierr=83
        return
      endif
      cropProp(4,iRec)=(cropProp(8,iRec)+cropProp(3,iRec))/two
c
c     User has entered MHSlope as a percentage and all M-H slopes 
c     are negative, divide by two to convert to ECsw
      MHSlope=-MHSlope/1.0D02/two
c     Convert ECe to ECe for soil water 
      cropProp(5,iRec)=two*cropProp(5,iRec)
c
c     Calculate the salinity at 50% yield reduction to get the
c     Pi50 coefficient
      cropProp(6,iRec)=(5.0D-01-(one-MHSlope*cropProp(5,iRec)))/MHSlope
c
c     Factor 3.9D02 is used to convert dS/m to osmotic head (cm) using
c     Handbook 60 Chapter 2 composite plot of osmotic pressure vs EC
      cropProp(5,iRec)=cropProp(5,iRec)*(-3.9D02)
      cropProp(6,iRec)=cropProp(6,iRec)*(-3.9D02)
      return
      end
c
c     Read the fundamental variables for the simulation defining the
c     geometry and the type of simulation being performed. Also, if
c     temperature-dependent organic-N decay is being modeled enter 
c     the monthly coefficients for temperature dependence.
      subroutine rdBasic(datPath,NumNP,nLayers,nTab,d,dz,
     ;                   NodeZ,x,depthLay,layer,layMat,lNitrate,
     ;                   NLossRat,TVarCoef,lTCoef,freeD,lPlant,
     ;                   hCritA,ierr,errStr)
      implicit none
      include 'defs.fi'
      integer NumNP,nLayers,nTab,layer,LayMat,ldt,ierr
      character*120 datPath
      logical lNitrate,freeD,lPlant,lTCoef
      double precision d,dz,depthLay,x,nodeZ,NLossRat,TVarCoef,
     ;                 hCritA
      dimension layer(nLayerD),layMat(nLayerD),TVarCoef(nMonths)
      dimension depthLay(NLayerD),x(NumNPD),nodeZ(NumNPD)
      character*80 errStr
c
c     Note TVarCoef are temperature reduction factors for NLossRat
c     The 2 entries for each month are month number and coefficient
      character*20 tokens(10)
      integer i,k,n,lenPath
      character*140 fname
      character*120 line
      character*20 BoolStr,monthStr
      logical found,feof
      integer strlen,getTok
      external strlen,getTok,getMonth
      double precision zero,one,getMonth
      data zero,one /0.0D00,1.0D00/
c
      errStr='Basic.in'
      ierr=0
      lTCoef=.FALSE.
      lenPath=strlen(datPath)
      fname = datPath(1:lenPath)//'\Basic.in'
      open(UNIT=nBasUnit,FILE=fname,ERR=900)
      call advance(nBasUnit,line,feof)
      ldt=4
      n = getTok(line,tokens,ldt)
      if (n.ne.ldt) then
        goto 901
      endif
      read(tokens(1),'(F10.4)') d
      read(tokens(2),'(F10.4)') dz
      read(tokens(3),'(I5)') nLayers
      read(tokens(4),'(I5)') NumNP
      x(1)=zero
c
c     Assume NTab - 100, if needed it can be a data entry
      NTab = 100
      x(NumNP)=zero
      do i=NumNP-1,1,-1
        x(i)=x(i+1)-dz
      enddo
      do i=1,NumNP
        nodeZ(i)=(NumNP-i)*d/(NumNP-1)
      enddo
      call advance(nBasUnit,line,feof)
c     
c     Plant roots, Free drainage, Nitrogen, N Loss Rate,hCritA
      ldt=5
      n = getTok(line,tokens,ldt)
      if (n.ne.ldt .and. n.ne.(ldt-1)) then
        goto 901
      endif
      read(tokens(1),'(A)') BoolStr
      if (BoolStr(1:1).eq.'T' .or. BoolStr(1:1).eq.'t') then
        lPlant=.TRUE.
      else
        lPlant=.FALSE.
      endif
      read(tokens(2),'(A)') BoolStr
      if (BoolStr(1:1).eq.'T' .or. BoolStr(1:1).eq.'t') then
        freeD=.TRUE.
      else
        freeD=.FALSE.
      endif
      read(tokens(3),'(A)') BoolStr
      if (BoolStr(1:1).eq.'T' .or. BoolStr(1:1).eq.'t') then
        lNitrate = .TRUE.
        if (n.ne.ldt) then
          goto 901
        endif
        read(tokens(4),*) NLossRat
      else
        lNitrate = .FALSE.
c       Read placeholder for NLossRat
        read(tokens(4),*) NLossRat
      endif
      read(tokens(5),*) hCritA
      call advance(nBasUnit,line,feof)
c
c Format is LayNum, MatNum, Layer Depth
      ldt=3
c Process the first line of data
      n = getTok(line,tokens,ldt)
      if (n.ne.ldt) then
        goto 901
      endif
c Layers have Layer(1) at the bottom of the profile
c Layer(nLayers) at the top. The depth of layer i
c is a positive number of cm from the soil surface
c to the bottom of layer i.
      do i=nLayers,1,-1
        read(tokens(1),'(I5)') layer(i)
        read(tokens(2),'(I5)') layMat(i)
        read(tokens(3),*) depthLay(i)
        if (i.gt.1) then
          read(nBasUnit,'(A)',ERR=903) line
          n = getTok(line,tokens,ldt)
          if (n.ne.ldt) then
            goto 901
          endif
        endif
      enddo
      if (.not.lNitrate) then
        return
      else
c
c       read(nBasUnit,'(A)',END=800,ERR=903) line
c       Preset the array for coefficient == 1
        do i=1,nMonths
          TVarCoef(i)=one
        enddo
        call advance(nBasUnit,line,feof)
c
c       If the user omitted the section on temperature-dependence
c       of mineralization rate coefficients then return with all
c       coefficients set to unity.
        if (feof) then
          return
        endif
        lTCoef=.TRUE.
        ldt=2
        n = getTok(line,tokens,ldt)
        if (n.ne.ldt) then
          goto 901
        endif
c
c       For those months entered by the user change the
c       array according to user-supplied coefficients
        do i=1,nMonths
          read(tokens(1),'(A)') monthStr
          read(tokens(2),*) TVarCoef(INT(getMonth(monthStr)))
          if (i.lt.nMonths) then
            read(nBasUnit,'(A)',END=800,ERR=902) line
            n = getTok(line,tokens,ldt)
            if (n.ne.ldt) then
              goto 901
            endif
          endif
        enddo
      endif
800   close(nBasUnit)
      return
900        ierr=1
          return
901        ierr=2
          return
902        ierr=3
          return
903        ierr=4
          return
      end

***********************************************************************
c     Read variables related to the computation of time steps and the
c     time period of the run
      subroutine rdTime(datPath,ntOut,dt0,dt,dtMin,dtMax,dtOut,dtOld,
     ;                  dateInit,dateFinal,tOutput,tInit,
     ;                  tFinal,ierr,errStr)
      implicit none
      include 'defs.fi'
      integer ntOut,daysInYr,ierr
      double precision dt0,dtMin,dtMax,dtOut,dt,dtOld
      double precision tOutput,tInit,tFinal
      character*20 dateInit,dateFinal
      dimension tOutput(nTOutD)
      character*120 datPath
      character*80 errStr
c
      integer n,i,ldt
      character*120 line
c
      character*20 tokens(10),dateStr
      integer*4 gettok,j,iPrint,lenPath,specPrnt,getDoy,dateVec(3)
      double precision hours,day,elapsed,convDate,zero
      external gettok,convDate,getDoy
      character*140 fname

      integer strlen
      logical feof
      external strlen
      data zero /0.0D00/
c
      ierr=0
      errStr='Time.in'
      lenPath=strlen(datPath)
      fname = datPath(1:lenPath)//'\Time.in'
      open(UNIT=nTimUnit,FILE=fname,ERR=900)
      call advance(nTimUnit,line,feof)
      ldt=1
      n = getTok(line,tokens,ldt)
      if (n.ne.ldt) then
        goto 901
      endif
      read(tokens(1),'(I10)') specPrnt
      call advance(nTimUnit,line,feof)
      ldt=6
      n = getTok(line,tokens,ldt)
      if (n.ne.ldt) then
        goto 901
      endif
      read(tokens(1),'(E13.6)') dt0
      read(tokens(2),'(E13.6)') dtMin
      read(tokens(3),'(E13.6)') dtMax
      read(tokens(4),'(E13.6)') dtOut
      read(tokens(5),'(A)') dateInit
      read(tokens(6),'(A)') dateFinal
c
c Set the starting dt
      dt = dt0
      dtOld=dt0
c
c Convert the date format for initial and final times to real numbers
      elapsed = convDate(dateInit,dateFinal,ierr)
      call dateParts(dateInit,dateVec(2),dateVec(3),dateVec(1),ierr)
      if (ierr.eq.1) goto 904
      tInit=dble(getDoy(dateVec))
      tFinal=tInit + elapsed
c
c In case of no list of specific print times
      if (specPrnt.eq.0.and.dtOut.gt.zero) then
c       The printTimes vector needs to be built for increment dtOut
        day=tInit
        i=1
        do while (day.lt.tFinal)
          tOutput(i)=day
          day=day+dble(int(dtOut))
          i=i+1
c         Terminate building vector if size is exceeded
          if (i.gt.nTOutD) then
            goto 10
          endif
        enddo
10      continue
        nTOut=i-1
        return !no print times specified so just return
      endif
c
c Read in the printout times, they are entered on
c successive lines of the Time.in file as a date
c in standard US format (mm/dd/yyyy) and a time
c in fractional hours (7.25 is 7:15 AM)
      iPrint=1
      call advance(nTimUnit,line,feof)
      ldt=2
      n = getTok(line,tokens,ldt)
      if (n.ne.ldt) then
        goto 901
      endif
      read(tokens(1),'(A)') dateStr
      read(tokens(2),*) hours
      tOutput(1) = tInit + convDate(dateInit,dateStr,ierr)
     ;             + hours/2.4D01
      if (ierr.gt.0) goto 904
      do while (iPrint.lt.nTOutD) !stop if exceeds vector size 
        read(nTimUnit,'(A)',END=800,ERR=902) line
        n=getTok(line,tokens,ldt)
        if (n.ne.ldt) then
          goto 901
        endif
        iPrint=iPrint+1
        read(tokens(1),'(A)') dateStr
        read(tokens(2),*) hours
        tOutput(iPrint) = tInit + 
     ;    convDate(dateInit,dateStr,ierr) + hours/2.4D01
        if (ierr.gt.0) goto 904
      enddo
800   nTOut = iPrint
      close(nTimUnit)
      return
900		ierr=1
      	return
901		ierr=2
      	return
902		ierr=3
      	return
903		ierr=4
      	return
904		ierr=5 !Error in date specification
      	return
      end
c
c     read soil properties for each designated soil from Soil.in
      subroutine rdSoil(datPath,nSoilPrp,matID,depth,
     ;                  alpha,Ks,nvg,thR,thS,ierr,errStr)
      implicit none
      include 'defs.fi'
      character*120 datPath
      double precision depth(nMatD),alpha(nMatD),
     ;  Ks(nMatD),nvg(nMatD),thR(nMatD),thS(nMatD)
      integer nSoilPrp,NLayers,matID(nMatD),ierr
      character*80 errStr
c
      integer ldt
      character*120 line
      character*20 tokens(10)
      character*140 fname
      integer strlen,lenPath,getTok,i,n
      logical feof
      external strlen,getTok
c
      errStr='Soil.in'
      lenPath=strlen(datPath)
      fname = datPath(1:lenPath)//'\Soil.in'
      open(UNIT=nSoilUnit,FILE=fname,ERR=900)
      call advance(nSoilUnit,line,feof)
      read(line,*,ERR=901) nSoilPrp      !Number of materials
      call advance(nSoilUnit,line,feof)
c
c One ID number and five soil properties per line of Soil.in,
c no comment lines beyond the header
      ldt=6
      n = getTok(line,tokens,ldt)
      if (n.ne.ldt) then
        goto 901
      endif
      do i=1,nSoilPrp
        read(tokens(1),'(I5)') matID(i)
        read(tokens(2),*) alpha(i)
        read(tokens(3),*) Ks(i)
        read(tokens(4),*) nvg(i)
        read(tokens(5),*) thR(i)
        read(tokens(6),*) thS(i)
c     Get the next line if available
        if (i.lt.nSoilPrp) then
          read(nSoilUnit,'(A)',ERR=902,END=903) line
          n = getTok(line,tokens,ldt)
          if (n.ne.ldt) then
            goto 901
          endif
        endif
      enddo
      close(nSoilUnit)
      return
900		ierr=1
      	return
901		ierr=2
      	return
902		ierr=3
      	return
903		ierr=4
      	return
      end
c
c     read soil properties related to chemistry for each designated
c     soil. These include bulk density, molecular diffusivity, and
c     dispersivity. Normally molecular diffusivity will be the same
c     for all soils.
      subroutine rdChmPar(datPath,NMat,matID,chemPar,ierr,errStr)
      implicit none
      include 'defs.fi'
      character*120 datPath
      double precision chemPar(nChemPar,NMat)
      integer NMat,matID(NMatD),ierr
      character*80 errStr
c
      integer ldt,Mater
      character*120 line
      character*20 tokens(10)
      character*140 fname
      integer strlen,lenPath,getTok,i,n
      logical feof
      external strlen,getTok
c
      errStr='ChemPar.in'
      lenPath=strlen(datPath)
      fname = datPath(1:lenPath)//'\ChemPar.in'
      open(UNIT=nChmUnit,FILE=fname,ERR=900)
      call advance(nChmUnit,line,feof)
      ldt=4
      n = getTok(line,tokens,ldt)
      if (n.ne.ldt) then
        goto 901
      endif
      do i=1,NMat
        read(tokens(1),'(I5)') Mater
c       Test material ID number for conformity with Soil.in IDs
        if (matID(i).ne.Mater) then
          goto 904
        endif
        chemPar(1,i)=dble(Mater)        !first entry is material number
        read(tokens(2),*) chemPar(2,i)  !bulk density (g/cm^3)
        read(tokens(3),*) chemPar(3,i)  !Dmolecular (cm^2/d)
        read(tokens(4),*) chemPar(4,i)  !Dispersivity (cm)
        if (i.lt.NMat) then
          read(nChmUnit,'(A)',ERR=902,END=903) line
          n = getTok(line,tokens,ldt)
          if (n.ne.ldt) then
            goto 901
          endif
        endif
      enddo
      close(nChmUnit)
      return
900		ierr=1
      	return
901		ierr=2
      	return
902		ierr=3
      	return
903		ierr=4
      	return
904     ierr=8
        return
      end

      subroutine rdChem(datPath,nSols,irriChem,solNum,ierr,errStr)
      implicit none
      include 'defs.fi'
      character*120 datPath
      double precision irriChem(nIrriC*2,MaxIrriD)
      integer nSols,solNum(MaxIrriD),ierr
      character*80 errStr
c
      integer ldt,n,k
      character*120 line
      character*20 tokens(10)
      character*140 fname
      integer strlen,lenPath,getTok
      logical feof
      external strlen,getTok
c
      errStr='ChemIrrig.in'
      lenPath=strlen(datPath)
      fname = datPath(1:lenPath)//'\ChemIrrig.in'
      open(UNIT=nIrrUnit,FILE=fname,ERR=900)
      call advance(nIrrUnit,line,feof)
      ldt=7
      n = getTok(line,tokens,ldt)
      if (n.ne.ldt) then
        goto 901
      endif
      k=1
      do while (k.le.MaxIrriD)
        read(tokens(1),'(I5)') solNum(k)
        read(tokens(2),*) irriChem(1,k)  !tracer top (arbitrary units, example mg/L)
        read(tokens(3),*) irriChem(2,k)  !salt top (mg/L) convert to mmolc/L
        irriChem(2,k)=irriChem(2,k)/5.6965D01 !Assuming MW for CaCl2 + 2NaCl
        read(tokens(4),*) irriChem(3,k)  !nitrate top convert mg/L to mmol/L using MW of N
        irriChem(3,k)=irriChem(3,k)/1.4D01
        read(tokens(5),*) irriChem(4,k)  !tracer bot
        read(tokens(6),*) irriChem(5,k)  !salt bot
    	irriChem(5,k)=irriChem(5,k)/5.6965D01
        read(tokens(7),*) irriChem(6,k)  !nitrate bot
        irriChem(6,k)=irriChem(6,k)/1.4D01
        if (k.lt.MaxIrriD) then
          read(nIrrUnit,'(A)',ERR=902,END=100) line
          n = getTok(line,tokens,ldt)
          if (n.ne.ldt) then
            goto 901
          endif
        else
          nSols=MaxIrriD
          close(nIrrUnit)
          return
        endif
        k=k+1
      enddo
100   close(nIrrUnit)
      nSols=k
      return
900   ierr=1
      return
901   ierr=2
      return
902   ierr=3
      return
      end
c
c Fertilization, tillage etc involve changes to the plow layer
c the timing and changes are entered in the file Input.in
***********************************************************************
      subroutine rdAdd(datPath,dateInit,tInit,nDiscret,atmFall,
     ;                 mineral,amtAppl,orgFrac,PsiApp,tApplic,
     ;                 tNextEvt,iApplic,ierr,errStr)
      implicit none
      include "defs.fi"
      integer nDiscret,iApplic,ierr
      character*120 datPath
      character*20 dateInit
      character*80 errStr
      double precision tInit,atmFall,mineral,amtAppl(nDisD),
     ;  orgFrac(nDisD),PsiApp(nDisD),tApplic(nDisD),tNextEvt
c
      integer ldt,i,n,nApplic
      character*120 line
      character*20 tokens(10),dateStr
      character*140 fname
      integer strlen,lenPath,getTok
      logical feof
      double precision appTime,convDate
      external strlen,getTok,convDate
c
      ierr=0
      errStr='Input.in'
      lenPath=strlen(datPath)
      appTime=7.0D00/2.4D01 !Assume application times start at 7AM
      fname = datPath(1:lenPath)//'\Input.in'
      open(UNIT=nAddUnit,FILE=fname,ERR=900)
      call advance(nAddUnit,line,feof)
      ldt=3
      n = getTok(line,tokens,ldt)
      if (n.ne.ldt) then
        goto 901
      endif
      read(tokens(1),'(F10.4)') atmFall
      atmFall=atmFall*1.0D-05
      read(tokens(2),'(F10.4)') mineral
      read(tokens(3),'(I5)') nDiscret
c
c Each discrete event has five entries:
c Application ID, amount, mineral fraction, organic fraction
c and time of application
      if (nDiscret.gt.0) then
        ldt=5
        call advance(nAddUnit,line,feof)
        n = getTok(line,tokens,ldt)
        if (n.ne.ldt) then
          goto 901
        endif
        do i=1,nDiscret
          read(tokens(1),'(I10)') nApplic !dummy entry for readability
          read(tokens(2),'(A)') dateStr
          tApplic(i)=tInit + 
     ;      convDate(dateInit,dateStr,ierr) + appTime
          if (ierr.gt.0) then
            ierr=93
          endif
          read(tokens(3),'(F10.4)') amtAppl(i) !kg/ha convert to g N/cm^2 in top layer
          amtAppl(i) = amtAppl(i)*1.D03/1.D08
          read(tokens(4),'(F10.4)') orgFrac(i)
          read(tokens(5),'(F10.4)') PsiApp(i)
          if (i.lt.nDiscret) then
            read(UNIT=nAddUnit,'(A)',END=903,ERR=902) line
            n = getTok(line,tokens,ldt)
            if (n.ne.ldt) then
              goto 901
            endif
          endif
        enddo
c     Set time for first event
        tNextEvt = tApplic(1)
        iApplic = 1
      else
        iApplic = 0
      endif
      close(nAddUnit)
      return
c       open file error
900		ierr=1
      	return
901		ierr=2
      	return
902		ierr=3
      	return
903		ierr=4
      	return
      end
c
c Boundary conditions
************************************************************************
      subroutine rdBounds(datPath,dateInit,tInit,bndCond,nBounds,
     ;                    bndTime,freeD,kodTop,kodBot,lNitrate,
     ;                    ierr,errStr)
      implicit none
      include 'defs.fi'
      integer nBounds,kodTop,kodBot,ierr
      character*120 datPath
      character*20 dateInit
      character*80 errStr
      logical freeD,lNitrate
      double precision tInit,bndCond(nBndVars,*),bndTime(*)
c
      integer ldt,n,i,ibound,chmBCTop,chmBCBot,iSol
      double precision zero,bndEvap,bndIrr,hours,convDate
      character*120 line
      character*140 fname
      character*20 tokens(12),dateStr
      integer strlen,lenPath,getTok
      logical feof
      external strlen,getTok,convDate
      data zero /0.0D00/
c
      ierr=0
      errStr='Bound.in'
      lenPath=strlen(datPath)
      fname = datPath(1:lenPath)//'\Bound.in'
      open(UNIT=nBndUnit,FILE=fname,ERR=900)
      call advance(nBndUnit,line,feof)
      ldt=1
      n=getTok(line,tokens,ldt)
      read(tokens(1),'(I5)',ERR=901) nBounds
      call advance(nBndUnit,line,feof)
      ldt=12
      n=getTok(line,tokens,ldt)
      if (n.ne.ldt) then
        goto 901
      endif
      do i=1,nBounds
        read(tokens(1),*) ibound     !unused variable, readability of data file
c       boundary condition date standard US format mm/dd/yyyy
        read(tokens(2),'(A)') dateStr
        read(tokens(3),*) hours       !time of day in fractional hours
        bndTime(i)=tInit + convDate(dateInit,dateStr,ierr)
     ;             + hours/2.4D01
        if (ierr.eq.1) then
          ierr=92
          return
        endif
        read(tokens(4),*) bndCond(1,i) !Potential evapotranspiration ETo (cm/day)
        read(tokens(5),*) bndCond(2,i) !Irrigation rate (cm/day)
        read(tokens(6),*) bndCond(3,i) !Crop Coefficient
        read(tokens(7),*) bndCond(4,i) !NUptake
c       Convert kg/ha/d to mmol/cm^2/d
        bndCond(4,i)=bndCond(4,i)/1.4D03 !MW NO3-N=14; mg/cm^2->kg/ha=100
        read(tokens(8),*) chmBCTop    !chemical transport BC Top
        read(tokens(9),*) chmBCBot    !chemical transport BC Bot
        bndCond(5,i)=dble(chmBCTop)
        bndCond(6,i)=dble(chmBCBot)
        read(tokens(10),*) iSol  !Solution ID number
        bndCond(7,i)=dble(iSol)
        if (lNitrate) then
          read(tokens(11),*) bndCond(8,i)!user-specified C* for N uptake
        else
          bndCond(8,i)=zero
        endif
        if (freeD) then
          bndCond(9,i)=zero
          kodBot=-1  !bottom condition for free drainage
        else
          read(tokens(12),*) bndCond(9,i)
          kodBot=1   !bottom condition for Dirichlet spec of hBot(t)
        endif
        if (i.lt.nBounds) then
          read(nBndUnit,'(A)',END=903,ERR=902) line
          n = getTok(line,tokens,ldt)
          if (n.ne.ldt) then
              goto 901
          endif
        endif
      enddo
c     We assume atmospheric boundary condition
      kodTop=-1
      close(nBndUnit)
      return
900   ierr=1
      return
901   ierr=2
      return
902   ierr=3
      return
903   ierr=4
      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 rdSched(datPath,iRecs,schedNam,cropSchd,dateInit,
     ;  dateFinal,tInit,tFinal,startNum,ierr,errStr)
      implicit none
      include 'defs.fi'
      integer iRecs,startNum,ierr
      character*120 datPath
      character*20 dateInit,dateFinal,str
      character*80 errStr
      double precision cropSchd,tInit,tFinal
      dimension cropSchd(nSchdPar,*)
      character*20		tokens(10)
      character*20		schedNam(*)
c
c     Local variables
      character*120	Line
      character*10	plntDate,harvDate,RDFDate
      character*140 fname
      double precision currVals(3),plntInit,lastHarv,elapsed,convDate,
     ;                 testDays
      integer	cropNum,nSched,iSched,n,i,lenPath,ldt,dateVec(3)
      integer	strlen,getTok,getDoy
      external	strlen,getTok,getDoy,convDate
c
      logical feof,lBare
      double precision zero,one
      data zero,one /0.0D00, 1.0D00/
c
      errStr='CropSched.in'
      lenPath=strlen(datPath)
      fname = datPath(1:lenPath)//'\CropSched.in'
      open(UNIT=nSchdUnit,File=fname,ERR=900)
      call advance(nSchdUnit,line,feof)
      ldt=1
      n = getTok(line,tokens,ldt)
      if (n.ne.ldt) then
        goto 901
      endif
      read(tokens(1),'(I5)') nSched
      if (nSched.eq.0) then
        return
      endif
c     Find out whether the fallow periods will be represented as bare
c     soil or as the standard cover crop
      call advance(nSchdUnit,line,feof)
      n=getTok(line,tokens,ldt)
      if (n.ne.ldt) then
        goto 901
      endif
      read(tokens(1),'(A)') str
      if (str(1:1).eq.'T'.or.str(1:1).eq.'t') then
        lBare=.TRUE.
      else
        lBare=.FALSE.
      endif
c
      do i=1,120
        Line(i:i)=' '
      end do
c Zero the first column of cropSchd for initialization of time variables
      do i=1,nSchdPar
        cropSchd(i,1)=zero
      enddo
c
c     Advance through comment lines and return the first data line
c
      iSched=0
      iRecs=0
      ldt=6
      call advance(nSchdUnit,line,feof)
c
c     Set the starting crop number to a signal value of -1 to indicate
c     that the starting crop number has not yet been set
      startNum=-1
c
c       Loop through data records for each crop 
20    iSched = iSched + 1 
      if (iSched.le.nSched) then
        n = getTok(line,tokens,ldt)
        if (n .ne. ldt) then
          goto 901
        end if
c
        if (iSched.eq.1) then
c         Assume that the first crop is planted on the same day as
c         specified for the initial time. If this is wrong then a
c         fallow record will be inserted.
          read(tokens(3),'(A10)') plntDate
          call dateParts(plntDate,dateVec(2),dateVec(3),dateVec(1),ierr)
          if (ierr.eq.1) goto 904
          plntInit=dble(getDoy(dateVec))
          if (plntInit.gt.tInit) then
c            Insert a fallow record to start
             iRecs = iRecs + 1
c            subtract 1 from the next crop start to get tEnd
             call mkFallow(cropSchd,iRecs,schedNam,tInit,plntInit-one,
     ;                     lBare)
c            set the starting crop number
             if (lBare) then
               startNum=0
             else 
c              cover crop
               startNum=1
             endif
          endif
        endif
        read(tokens(2),'(I5)') cropNum
c
c       If the starting crop number was not zero or one then set
c       it to the crop number just read because the crop was 
c       planted on the same day as the start of the simulation
        if (startNum.eq.-1) then
          startNum=cropNum
        endif
        read(tokens(3),'(A10)') plntDate
        read(tokens(4),'(A10)') harvDate
        read(tokens(5),'(A10)') RDFDate
c
c       If this is not the first crop then determine whether a record
c       of fallow period needs to be inserted into the crop schedule
        call cropDates(plntDate,harvDate,RDFDate,plntInit,currVals,
     ;                 ierr)
        if (ierr.eq.1) goto 904 !date format error
c
c       Testdays is the number of days from tInit to the planting 
c       of the current crop. If this is larger than plntInit then
c       we need a fallow period to make up the difference 
        testDays=tInit + convDate(dateInit,plntDate,ierr)
        if (ierr.eq.1) goto 904
        if (testDays.gt.plntInit) then
          iRecs = iRecs + 1
c           Insert a fallow record (cropNum==0) up until the day before
c           the planting date of the current data entry
          call mkFallow(cropSchd,iRecs,schedNam,plntInit,testDays-one,
     ;                  lBare)
c           Reset the harvest day to the last day of fallow
          plntInit=testDays
          call cropDates(plntDate,harvDate,RDFDate,plntInit,currVals,
     ;                   ierr)
          if (ierr.eq.1) goto 904
c         else just add the record to the crop schedule
c           the next crop planting date is immediately after the
c           prior harvest, this might apply to a cover crop or similar
        endif
        iRecs=iRecs+1
        read(tokens(1),'(A20)') schedNam(iRecs)
        cropSchd(1,iRecs) = dble(cropNum)
        do i=2,4
          cropSchd(i,iRecs)=currVals(i-1) !Current crop days
        enddo
        do i=5,n-1
c         remaining variables, n-1 occurs because first data entry
c         is the crop name which is a string so use tokens(i+1)
          read(tokens(i+1),*) cropSchd(i,iRecs)
        end do
c
c       next planting date will be one day after the current harvest
c       unless the next data entry indicates a fallow period. 
c       currVals(2) stores the harvest day of the current crop
        plntInit=currVals(2)+one
c
c       Get another line from file, end of file is an error
c       because there are less than nSched records
        if (iSched.lt.nSched) then
          read(nSchdUnit,'(A)',END=903,ERR=902) line
        endif
        goto 20
      else
        close(nSchdUnit)
        return
      endif
c
c     When kRoot=1, lTemp must be .TRUE. and temperature data must
c     be specified, otherwise root growth will not occur
c
900   ierr=1
      return
c     Incorrect number of entries on a line
901   ierr=2
      return
c     File read error
902   ierr=3
      return
c     Incorrect number of lines in file
903   ierr=4
      return
c     Date format problem
904   ierr=5
      return
905   ierr=6
      return
      end

************************************************************************
      subroutine mkFallow(cropSchd,iRecs,schedNam,tStart,tEnd,lBare)
      implicit none
      include "defs.fi"
      double precision cropSchd,tStart,tEnd
      dimension cropSchd(nSchdPar,*)
      integer iRecs
      logical lBare
      character*20		schedNam(*)
      double precision zero,one,two
      data zero,one,two /0.0D00, 1.0D00, 2.0D00/
c
      cropSchd(2,iRecs) = tStart
      cropSchd(3,iRecs) = tEnd
      if (lBare) then
	    schedNam(iRecs)='Bare_Soil'
		cropSchd(1,iRecs)=zero !cropID for bare soil
		cropSchd(4,iRecs)=zero
		cropSchd(5,iRecs)=zero
      else
        schedNam(iRecs)='Fallow'
        cropSchd(1,iRecs) = one !cropID for fallow
        cropSchd(4,iRecs) = (tStart*3.0D00+tEnd)/4.0D00 !RDFDay assumed
c                         to represent 1/4 of the fallow period
        cropSchd(5,iRecs) = two   !assume grass has 2cm start
      endif
      return
      end

************************************************************************
      subroutine bldTFix(nBounds,nTout,nDiscret,bndTime,tOutput,
     ;                   tApplic,tFix,tFinal,nTFix,ipFix,iPrint,ierr)
      implicit none
      include "defs.fi"
c
      character*120 datPath
      integer nBounds,nTout,nDiscret,nTFix,ipFix,iPrint,ierr
      double precision bndTime(nBounds),tOutput(nTout),tApplic(nDiscret)
      double precision tFix(nTFixD)
      double precision tFinal
c
      integer k,ipBound,ipInput,ipPrint
      double precision tLast,tFixCurr,nextBnd,nextOut,nextInput
      double precision zero,x,y,err
      data zero /0.0D00/
      err(x,y)=dabs(1.-x/y)-1.0D-04
      ierr=0
c Build the TFix array considering each of the three events that can cause a
c reset of time-stepping. These are boundary condition record, fertilizer input
c record and printout time record.
      ipBound=0
      ipInput=0
      ipPrint=0
c
c Find the last time among the three records
      if (nDiscret.gt.0) then
        tLast=max(bndTime(nBounds),tOutput(nTout),tApplic(nDiscret))
      else
        tLast=max(bndTime(nBounds),tOutput(nTout))
      endif
      k=1
      tFixCurr=zero
      do while (tFixCurr.lt.tLast)
c   Are there any equalities among the three times, first test for all three
        ipBound=ipBound+1
        ipInput=ipInput+1
        ipPrint=ipPrint+1
        if (ipBound.gt.nBounds) then
          nextBnd=tLast
          ipBound=nBounds+1
        else
          nextBnd=bndTime(ipBound)
        endif
        if (ipInput.gt.nDiscret) then
          nextInput=tLast
          ipInput=nDiscret+1
        else
          nextInput=tApplic(ipInput)
        endif
        if (ipPrint.gt.nTout) then
          nextOut=tLast
          ipPrint=nTout+1
        else
          nextOut=tOutput(ipPrint)
        endif
c Test for equality of all three
        if ((err(nextBnd,nextInput).lt.zero).and.
     ;        (err(nextBnd,nextOut).lt.zero)) then
            tFix(k)=nextBnd
c
c Test for equality of each of the pairs
        elseif (err(nextBnd,nextInput).lt.zero) then
            if (nextBnd.lt.nextOut) then
              tFix(k)=nextBnd
              ipPrint=ipPrint-1
            else 
              tFix(k)=nextOut
              ipBound=ipBound-1
              ipInput=ipInput-1
            endif
        elseif (err(nextBnd,nextOut).lt.zero) then
            if (nextBnd.lt.nextInput) then
              tFix(k)=nextBnd
              ipInput=ipInput-1
            else
              tFix(k)=nextInput
              ipBound=ipBound-1
              ipPrint=ipPrint-1
            endif
        elseif (err(nextInput,nextOut).lt.zero) then
            if (nextInput.lt.nextBnd) then
              tFix(k)=nextInput
              ipBound=ipBound-1
            else
              tFix(k)=nextBnd
              ipInput=ipInput-1
              ipPrint=ipPrint-1
            endif
c
c Most likely case, all three times are different, find the smallest
        else
c           Get the minimum of the three times, find out which it was
c           and decrement the other two
            tFix(k)=min(nextBnd,nextInput,nextOut)
            if (tFix(k).eq.nextBnd) then
              ipInput=ipInput-1
              ipPrint=ipPrint-1
            elseif (tFix(k).eq.nextInput) then
              ipBound=ipBound-1
              ipPrint=ipPrint-1
            else
              ipBound=ipBound-1
              ipInput=ipInput-1
            endif
        endif
        tFixCurr=tFix(k)
        if (k.lt.nTFixD-1) then
          k=k+1
        else
          goto 901
        endif
      enddo !while
c Set the final fixed time from the time arrays
      tFix(k)=tLast
c Add tFinal to the fixed-time array if the last entry in the array
c is less than tFinal
      if ((err(tFix(k),tFinal).lt.zero).and.(tFix(k).lt.tFinal)) then
        k=k+1
        tFix(k)=tFinal
      endif
      nTFix=k
      ipFix=1
      iPrint=1
      return
901   ierr=61
      return
      end

************************************************************************
      subroutine genMat(nMat,MatID,alpha,Ks,nvg,thR,thS,ParD,ierr,
     ;                  errStr)
      implicit none
      include 'defs.fi'
      integer nMat,matID(NMat),ierr
      double precision alpha,Ks,nvg,thR,thS,ParD(NTPar,*)
      dimension alpha(nMat),Ks(nMat),
     ;          nvg(nMat),thR(nMat),thS(nMat)
      character*80 errStr
c
c nTPar==10 is an assumption (defined in defs.fi)
      integer i
c
      ierr=0
      if (nMat.gt.nMatD) then
        ierr=82
        return
      end if
c Build the ParD array of soil hydraulic properties
      do i=1,nMat
        MatID(i)=i  !set material ID number to index
          ParD(1,i) = thR(i)
          ParD(2,i) = thS(i)
          ParD(3,i) = alpha(i)
          ParD(4,i) = nvg(i)
          ParD(5,i) = Ks(i)
          ParD(6,i) = 0.5D00  !Fixed value for exponent
          ParD(7,i) = thS(i)
          ParD(8,i) = thR(i)
          ParD(9,i) = thS(i)
          ParD(10,i) = Ks(i)
      enddo
      return
      end
c
c     Modifying factors for water flow were included in the
c     earlier version of ENVIRO-GRO but had not been implemented.
c     This is included as a record of what potentially might be
c     included
      subroutine setModif(NumNP,Ah,AK,ATh,AThS,AKS)
      implicit none
      integer NumNP
      double precision Ah,AK,ATh,AThS,AKS
      dimension Ah(NumNP),AK(NumNP),ATh(NumNP),AThS(NumNP),
     ;  AKS(NumNP)
c
      integer i
      double precision one
      data one /1.0D00 /
c
c     Modifying factors for water flow calculations all set 
c     to one.
      do i=1,NumNP
        Ah(i) = one
        AK(i) = one
        ATh(i) = one
        AThS(i) = one
        AKS(i) = one
      enddo
      return
      end

************************************************************************
c
c     Set the initial pressure heads for all nodes
      subroutine nodeHead(NumNP,hOld,hTemp,hNew,hTop,hBot)
      implicit none
      integer NumNP
      double precision hOld(NumNP),hNew(NumNP),hTemp(NumNP),hTop,hBot
c
      double precision one
      data one / 1.D00 /
      integer i
c
c     Initialize pressure heads using a simple linear interpolation
c     between the pressure head at the top (hTop) and the bottom (hBot)
      do i=NumNP,1,-1
        hOld(i)=hTop - (hTop-hBot)*
     ;          (dble(NumNP)-dble(i))/(dble(NumNP)-one)
        hTemp(i)=hOld(i)
        hNew(i)=hOld(i)
      enddo
      return
      end

************************************************************************
        subroutine nodeBld(NumNP,NLayers,NMat,layerDep,layMat,MatNum,
     ;                     dz,saltLay,traceLay,NMassFix,Conc,CFix,
     ;                     orgNFrac,Psi,PsiZero,chemPar,prevMass,ierr)
        implicit none
        include "defs.fi"
        integer          NumNP,NLayers,NMat,layMat,MatNum,ierr
        dimension        layMat(NLayers),MatNum(NumNP)
        double precision dz,layerDep(NLayers),chemPar(nChemPar,NMat)
        double precision Conc(NConcD,NumNP),CFix(NFixC,NumNP),
     ;                   Psi(NumNP)
        double precision saltLay(NLayers),traceLay(NLayers),
     ;                   NMassFix(NLayers),PsiZero(NLayers),
     ;                   orgNFrac(NLayers),prevMass(NLayers)
c
        integer i
        double precision zero,one,NConcMin,NConcOrg,top,
     ;                   thick,layThick
        dimension NConcOrg(NLayerD)
        external layThick
        data zero,one / 0.0D00,1.0D00 /
c
c       Convert the salt concentration in mg/L to mmolc/L
        call convSalt(NLayers,saltLay)
		top=zero
		do i=NLayers,1,-1
		  thick = abs(layThick(NLayers,i,top,layerDep))
c         convert from kg/ha -> g/cm^2 multiply by 1e-05
          NConcOrg(i)=(orgNFrac(i)*NMassFix(i)*1.0e-05)/
     ;                (thick*chemPar(2,layMat(i)))  !g/g Soil
          top=layerDep(i)
		enddo
c
c       Build the MatNum vector specifying material ID for each node
        call transInt(NumNP,NLayers,layerDep,dz,MatNum,layMat,ierr)
c
c       Each array contains initialization data by layer so the
c       dimension of the input array is NLayers.
c       The output array contains the same data but it has been
c       copied to form a nodal array with the same values
        call transNod(NumNP,NLayers,layerDep,dz,Conc,1,NConcD,
     ;                traceLay,ierr)
        call transNod(NumNP,NLayers,layerDep,dz,Conc,2,NConcD,
     ;                saltLay,ierr)
        call transVar(NumNP,NLayers,layerDep,dz,Psi,PsiZero,
     ;                ierr)
c
c       Set CFix for organic N fraction
        call transNod(NumNP,NLayers,layerDep,dz,CFix,1,NFixC,
     ;                NConcOrg,ierr)
        call setFracN(NumNP,NLayers,CFix,1,dz,orgNFrac,layerDep,
     ;                ierr)
c
c       Calculate the mass of orgN in each layer
        call layNMass(NumNP,NLayers,NMat,MatNum,dz,chemPar,CFix,
     ;                layerDep,prevMass,ierr)
c
c       cFix(2,:) is initialized to zero in layConc under the assumption
c       that all mineral N initially present will be dissolved into 
c       soil water
        return
        end

************************************************************************
        subroutine transNod(NumNP,NLayers,layerDep,dz,nodalVar,species,
     ;                      ldc,layerVar,ierr)
        implicit none
        integer NumNP,NLayers,species,ldc,ierr
        double precision dz
        double precision layerDep(NLayers),nodalVar(ldc,NumNP),
     ;                   layerVar(NLayers)
c
        integer layNum,layOfNod,j
c
c       Vector concentrations are entered so that the top node has 
c       index NumNP and the bottom node has index 1
        do j=NumNP,1,-1
          layNum = layOfNod(NumNP,layerDep,dz,j,NLayers,ierr)
          nodalVar(species,j) = layerVar(layNum)
        enddo
        return
        end

        subroutine transVar(NumNP,NLayers,layerDep,dz,nodalVar,
     ;                      layerVar,ierr)
        implicit none
        integer NumNP,NLayers,species,ldc,ierr
        double precision dz
        double precision layerDep(NLayers),nodalVar(NumNP),
     ;                   layerVar(NLayers)
c
        integer j,layNum,layOfNod
c
        do j=NumNP,1,-1
          layNum = layOfNod(NumNP,layerDep,dz,j,NLayers,ierr)
          nodalVar(j) = layerVar(layNum)
        enddo
        return
        end

************************************************************************
        subroutine transInt(NumNP,NLayers,layerDep,dz,nodalVar,layerVar,
     ;                      ierr)
        implicit none
        integer NumNP,NLayers,ierr
        double precision dz,layerDep(NLayers)
        integer nodalVar(NumNP),layerVar(NLayers)
c
        integer layNum,layOfNod,j
c
        do j=NumNP,1,-1
          layNum = layOfNod(NumNP,layerDep,dz,j,NLayers,ierr)
          nodalVar(j) = layerVar(layNum)
        enddo
        return
        end

        subroutine setFracN(N,NLayers,CFix,nType,dz,specNFrac,
     ;                      layerDep,ierr)
        implicit none
        include "defs.fi"
        integer N,NLayers,nType,ierr
        double precision CFix,specNFrac,dz,layerDep
        dimension CFix(nFixC,N),specNFrac(NLayers),layerDep(NLayers)
c
        integer j,layOfNod,layNum
c       multiply the total fixed concentration of N both organic
c       and mineral by the fractions of organic and mineral
        do j=N,1,-1
          layNum = layOfNod(N,layerDep,dz,j,NLayers,ierr)
          CFix(nType,j)=CFix(nType,j)*specNFrac(layNum)
        enddo
        return
        end
c
c Convert dates into number of days following the planting day of the
c first crop
        subroutine cropDates(plntDate,harvDate,RDFDate,
     ;                       plntInit,currVals,ierr)
        implicit none
        integer ierr
        character*10 plntDate,harvDate,RDFDate
        double precision plntInit,currVals
        dimension currVals(3)
        double precision convDate
        external convDate
c
        currVals(1)=plntInit !planting date
        currVals(2)=plntInit+convDate(plntDate,harvDate,ierr)
        currVals(3)=plntInit+convDate(plntDate,RDFDate,ierr)
        return
        end
c
c dateStr is assumed to be a standard 10-character string of the
c form mm/dd/yyyy. For other formats modifications are necessary.
************************************************************************
        double precision function convDate(dateStr1,dateStr2,ierr)
        character*(*) dateStr1,dateStr2
        integer       index,ierr
c
c Local variables
        integer mm1,dd1,yyyy1,mm2,dd2,yyyy2
        integer offset1,indexStr,getTok
        integer elapsDays
        external elapsDays
        character*1 test1
c
        ierr=0
        call dateParts(dateStr1,mm1,dd1,yyyy1,ierr)
        call dateParts(dateStr2,mm2,dd2,yyyy2,ierr)
        convDate = dble(elapsDays(yyyy1,mm1,dd1,yyyy2,mm2,dd2,ierr))
        return
        end

************************************************************************
        subroutine dateParts(dateStr,mm,dd,yyyy,ierr)
        implicit none
        character*(*) dateStr
        integer mm,dd,yyyy,ierr
c
        character*20 Tokens
        dimension Tokens(12)
        integer indexStr,getTok,offset1,n
        external indexStr,getTok
c
        character*1 test1
        test1='/'
c       Find the offsets to '/' and replace with blanks
        offset1=indexStr(test1,dateStr)
        if (offset1.eq.0) then
c         Assume the blanks have already been inserted, if
c         they have not that is an error
          n = getTok(dateStr,Tokens,3)
          if (n.ne.3) goto 900
        else
c         convert slashes to blanks
          dateStr(offset1:offset1)=' '
          offset1=indexStr(test1,dateStr)
          dateStr(offset1:offset1)=' '
c         the constant 12 is used as a dimension for Tokens
          n = getTok(dateStr,Tokens,3)
          if (n.ne.3) goto 900
        endif
        read(Tokens(1),'(I5)') mm
        read(Tokens(2) ,'(I5)') dd
        read(Tokens(3) ,'(I5)') yyyy
        return
900     ierr=1
        return
        end

      subroutine setBool(wLayer,TopInf,lTable)
      implicit none
      logical wLayer,topInf,lTable
      wLayer=.TRUE.
      topInf=.TRUE.
      lTable=.TRUE.
      return
      end

      subroutine setFlags(lFixFlag,lPrnFlag)
      implicit none
      logical lFixFlag,lPrnFlag
      lFixFlag=.FALSE.
      lPrnFlag=.FALSE.
      return
      end

      subroutine initLay(NLayers,sLayTopQ,
     ;                    sLayBotQ,sLaySink)
      implicit none
      integer NLayers
      double precision sLayTopQ,sLayBotQ,sLaySink
      dimension sLayTopQ(NLayers),sLayBotQ(NLayers),
     ;          sLaySink(NLayers)
c
      integer i
      double precision zero
      data zero /0.0D00/
c
      do i=NLayers,1,-1
        sLayTopQ(i)=zero
        sLayBotQ(i)=zero
        sLaySink(i)=zero
      enddo
      return
      end

      subroutine layConc(N,NLayers,NMat,z,theta,layerDep,NMassMin,
     ;                   Conc,cFix,chemPar,MatNum,ierr)
      implicit none
      include "defs.fi"
      integer MatNum,NMat,N,NLayers,ierr
      double precision z,theta,layerDep,NMassMin,Conc,cFix,chemPar
      dimension z(N),theta(N),layerDep(NLayers),NMassMin(NLayers),
     ;          Conc(NConcD,N),cFix(nFixC,N),chemPar(nChemPar,NMat),
     ;          MatNum(N)
c
      integer i,j
      double precision massMinN,c,sumVol,dz,zero,two
      dimension c(NLayers)
      data zero,two /0.0D00, 2.0D00/
c
      j = NLayers
      sumVol = zero
      do i=N-1,1,-1
        dz = z(i+1)-z(i)
        sumVol = sumVol + dz*(theta(i+1)+theta(i))/two
        if (abs(z(i)).ge.layerDep(j)) then
c         convert kg/ha to a layer concentration in mmol/L
          c(j)=NMassMin(j)*1.0D06/(sumVol*1.0D-03*1.4D01*1.0D08)
          sumVol=zero
          j=j-1
        endif
c       Initialize cFix for mineral N to zero because all 
c       mineral N has been dissolved upon return from this
c       subroutine
        cFix(2,i)=zero
      enddo
      cFix(2,N)=zero
c
c     Transfer layer concentrations in c to Conc array
      call transNod(N,NLayers,layerDep,dz,Conc,3,NConcD,
     ;              c,ierr)
      return
      end

c     Calculate the mass of plant-available nitrogen in each layer 
c     from the nodal N concentrations in the Conc array
      subroutine layMass(N,NLayers,z,theta,layerDep,Conc,
     ;                   massNLay,ierr)
      implicit none
      include "defs.fi"
      integer N,NLayers,ierr
      double precision z,theta,layerDep,Conc,massNLay
      dimension z(N),theta(N),layerDep(NLayers),
     ;          Conc(nConcD,N),massNLay(NLayers)
c
      integer i,j
      double precision zero,dz
      data zero /0.0D00/
c
      j=N-1
      do i=NLayers,1,-1
        massNLay(i)=zero
        do while (abs(z(j)).le.layerDep(i))
          dz = abs(z(j+1)-z(j))
          massNLay(i) = massNLay(i) + 
     ;           theta(j)*Conc(3,j)*dz*1.D+05*1.4D01/1.D+06
          j=j-1
        enddo
      enddo
      return
      end

      subroutine initNSum(NLayers,dayNSum,tLeachN)
      implicit none
      integer NLayers
      double precision dayNSum,tLeachN
      dimension dayNSum(NLayers)
c
      integer i
      double precision zero
      data zero /0.0D00/
c
      do i=NLayers,1,-1
        dayNSum(i)=zero
      enddo
      tLeachN=zero
      return
      end
          

      subroutine convSalt(NLayers,saltLay)
      implicit none
      integer NLayers
      double precision saltLay
      dimension saltLay(NLayers)
c
      integer i
      double precision MWSalt
c     molecular weight of a mixture of CaCl2,NaCl
c     having equal charge concentration for each component
      data MWSalt /5.6965D01/
c
      do i=NLayers,1,-1
        saltLay(i)=saltLay(i)/MWSalt
      enddo
      return
      end
c
c     Set pressure heads at the top and bottom of the profile
      subroutine setHBnd(NLayers,NMat,iModel,watType,freeD,wTop,
     ;                   wBot,ParD,layMat,bndCond,hBot,hTop)
      implicit none
      include "defs.fi"
      logical freeD
      integer watType,layMat,NLayers,NMat,iModel
      double precision wTop,wBot,ParD,bndCond,hTop,hBot
      dimension ParD(NTPar,NMat),layMat(NLayers),
     ;          bndCond(nBndVars,*)
c
      double precision fh
      external fh
c
      if (watType.eq.1) then
        hTop=wTop
        if (freeD) then
          hBot=wBot
        else
c         boundary condition is specified head
          hBot=bndCond(9,1)
        endif
      elseif (freeD) then
c       The hTop and hBot must be calculated from the entered theta values
        hTop=fH(iModel,wTop,ParD,layMat(NLayers))
        hBot=fH(iModel,wBot,ParD,layMat(1))
      else !bottom BC is hBot specified so override wBot using the BC
        hTop=fH(iModel,wTop,ParD,layMat(NLayers))
        hBot=bndCond(9,1)
      endif
      return
      end