***********************************************************************
*     MAIN for TSTEROD
***********************************************************************
c
c     +++  PURPOSE +++
c
c     To start a standalone version of the HYDROLOGY submodel
c
c
c     +++ ARGUMENT DECLARATIONS +++
c     
c     + + + GLOBAL COMMON BLOCKS + + +
C Includes
      include 'p1werm.inc'
      include 'c1glob.inc'
      include 'd1glob.inc'
      include 'm1subr.inc'
      include 'm1flag.inc'
      include 's1layd.inc'
      include 's1layr.inc'
      include 's1dbh.inc'
      include 's1phys.inc'
      include 's1surf.inc'
      include 'h1hydro.inc'
      include 'h1temp.inc'
      include 'h1db1.inc'
      include 'h1scs.inc'
      include 'w1wind.inc'
      include 'main/main.inc'
c
c     + + + LOCAL COMMON BLOCKS      
c
c     ++++ ARGUMENT DEFINITIONS +++
c
c     +++ SUBROUTINES CALLED+++
C     
c     ++++ LOCAL VARIABLES +++
	  integer isr
      integer mxdasm
      integer idx
      integer cd, cm, cy
      real grad
	  real temp
      real begtim, endtim
      
c
c     +++ END SPECIFICATIONS +++
c

C use isr so that we don't change all of the parameters in the function calls

      call cpu_time(begtim)
      am0ifl = .true.
      isr = 1
      am0csr = 1
      am0hfl = 2
      am0dfmfl = 0
      am0drmfl = 1
      nsubr = 1

C read in parameters

      write(*,*) 'Enter average slope '
      read(*,*) amrslp(isr)
      write(*,*) 'Enter crop biomass cover fraction '
      read(*,*) acftcv(isr)
      write(*,*) 'Enter crop leaf area index '
      read(*,*) acrlai(isr)
      write(*,*) 'Enter total flat biomass '
      read(*,*) bdmft
      write(*,*) 'Enter root depth '
      read(*,*) aczrtd(isr)
      write(*,*) 'Enter average daily wind speed '
      read(*,*) awudav
      write(*,*) 'Enter minimum air temperature '
      read(*,*) awtdmn
      write(*,*) 'Enter maximum air temperatures '
      read(*,*) awtdmx
      write(*,*) 'Enter number of days to simulate '
      read(*,*) mxdasm
      write(*,*) 'Enter ifc file name '
      read(*,*) sinfil
      write(*,*) 'Enter cli_gen file name '
      read(*,*) clifil
      call fopenk (21, clifil, 'old')
      call fopenk(25, 'h1.out', 'unknown')
      call fopenk(26, 'h2.out', 'unknown')

C read in the ifc file        

      call inpsub

C re-calc the layers
      
      call recalc

C calculate layer depths

      awtdav = (awtdmn + awtdmx) / 2
      amzele = 35.0
      write(*,*) ' awtdav ', awtdav
      aszlyd(1,isr) = aszlyt(1,isr)
      do 20 idx = 2,nslay(isr)
        aszlyd(idx,isr) = aszlyd(idx-1,isr) + aszlyt(idx,isr)
   20 continue

C initialize the hydro variables
      call hydrinit(isr)

C actually call hydro, printing out variables as we go

      daysim = 0
      write(25,1001) daysim, (ahrwca(idx,isr), idx=1,10)
      write(25,1004) daysim, (ahrwc(idx,isr), idx=1,10)
      write(25,1002) daysim, (ahrwcf(idx,isr), idx=1,10)
      write(25,1003) daysim, (ahrwcw(idx,isr), idx=1,10)
      write(25,1007) daysim, (ahrwcs(idx,isr), idx=1,10)

      write(25,1005) 
C *** 1005 format('      arrwca    arrwca    arrwca    arrwca  ',
 1005 format('      awzdpt  ',
     *  '   ahrwc     ahrwc     ahrwc     ahrwc  ')
      write(*,*) ' ah0cnp, ah0cng ', ah0cnp(isr), ah0cng(isr)
      write(25,1006) 0, 0.0, (ahrwc(idx,isr), idx=1,7)
      do 10 daysim=(1+365*2000),mxdasm+365*2000
        awzdpt=0.0
        call caldat (daysim, cd, cm, cy)

        call getcli(cd, cm, cy, awzdpt, awtdmx, awtdmn, grad, awtdpt)
        aweirr = grad * 0.04186

C ***        if (mod(daysim, 25).eq.15) awzdpt=20.0

        call hydro( nslay(isr), amrslp(isr),
     &                  acftcv(isr), acrlai(isr),
     &                  bdmft, aczrtd(isr), ahfwsf(isr),
     &                  aszlyd(1, isr), asdblk(1, isr),
     &                  ahrwc(1, isr), ahrwcs(1, isr),
     &                  ahrwcf(1, isr), ahrwcw(1, isr),
     &                  ah0cb(1,isr), aheaep(1,isr),
     &                  asfsan(1,isr), asfsil(1,isr), asfcla(1,isr),
     &                  ah0cng(isr), ah0cnp(isr),
     &                  ahzper(isr), ahzirr(isr), ahzrun(isr),
     &                  awudav, ahrsk(1, isr),
     &                  ahtsmx(1, isr), ahtsmn(1, isr),
     &                  ahrwc0(1, isr), daysim,
     &                  asfald(isr), asfalw(isr), aszlyt(1,isr),
     *                  awzdpt, awtdmx, awtdmn, ahzwid(isr) )
C
      temp =  sum(ahrwc0(1:24,isr))/(24 *  asdblk(1,isr))


C ***      write(27, 1010) bwzdpt, snwci, snmlt, bhzirr, bhzrun, dinf within the subroutine (inputs and swnci before calculations, calculated values and swnci after calculations). In a full WEPS run, print out the global variables awzdpt, ahzsno, ahzsnd, ahzirr, ahzrun to see if the parameters are passed and updated correctly. Make sure curve numbers are enabled, ie. input file should have curve numbers greater than zero so the runoff portion is tested. 

      write(25,1006) daysim, awzdpt, temp, (ahrwc(idx,isr), idx=1,7)
1006  format(i8, 9f10.4) 
   10 continue
      write(25,1002) daysim, (ahrwcf(idx,isr), idx=1,10)
      write(25,1001) daysim, (ahrwca(idx,isr), idx=1,10)
1001  format(i4, ' ahrwca ', 10f8.4) 
1004  format(i4, ' ahrwc  ', 10f8.4) 
1002  format(i4, ' ahrwcf ', 10f8.4) 
      write(25,1003) daysim, (ahrwcw(idx,isr), idx=1,10)
1003  format(i4, ' ahrwcw ', 10f8.4) 
      write(25,1007) daysim, (ahrwcs(idx,isr), idx=1,10)
1007  format(i4, ' ahrwcs ', 10f8.4) 
 
      call cpu_time(endtim)

      write(*,*) ' times ', begtim, endtim, endtim-begtim
      write(25,*) ' times ', begtim, endtim, endtim-begtim

      close (25)
      end