c     file name: hinit.for

      subroutine   hinit (layrsn, bsdblk, bhrwc, bhrwcs, bhrwcf,
     &                    bhrwcw, bh0cb, bheaep, bhrsk, bsfsan, bsfsil,
     *                    bsfcla, bszlyd, bszlyt)

c     + + + PURPOSE + + +
c     This subroutine controls the initialization of the HYDROLOGY
c     sumbodel of WEPS.  The program initializes the depth variables
c     of the soil simulation layers and converts the soil water
c     content variables from mass basis to volume basis.
c     DATE:  09/16/93
c     MODIFIED:  12/13/93
c     MODIFIED:  07/28/95
c     MODIFIED:  07/29/95
*     This change was done to determine the average of soil
*     properties from the first simulation layer (10 mm) and the second
*     uppermost simulation layer (40 mm).  The average for the new
*     uppermost simulation layer for the HYDROLOGY submodel (50 mm thick)
*     Using 50 mm as the thickness of the uppermost simulation layer for the
*     HYDROLOGY submodel will increase the speed of simulation and reduce the
*     the potential for errors.
*
*
c     + + + KEYWORDS + + +
c     initialization, hydrology

c     + + + ARGUMENT DECLARATIONS + + +
      integer layrsn
      real bsdblk(*)
      real bhrwc(*)
      real bhrwcs(*)
      real bhrwcf(*)
      real bhrwcw(*)
      real bh0cb(*)
      real bheaep(*)
      real bhrsk(*)
      real bsfsan(*)
      real bsfsil(*)
      real bsfcla(*)
      real bszlyd(*)
      real bszlyt(*)

c     + + + ARGUMENT DEFINITIONS + + +
c     layrsn - Number of soil layers used in simulation
c     bsdblk     - Soil bulk density (Mg/m^3)
c     bhrwc     - Soil water content (mg/mg)
c     bhrwcs    - Soil water content at saturation (mg/mg)
c     bhrwcf   - Soil water content at field capacity (mg/mg)
c     bhrwcw   - Soil water content at wilting point (mg/mg)
c     bh0cb     - Power of campbell's water release curve model (unitless)
c     bheaep   - Soil air entry potential (j/kg)
c     bsfsan   - Sand fractions
c     bhrsk   - Saturated hydraulic conductivity (m/s)
c     bsfsil   - Silt fractions
c     bsfcla   - Clay fractions

c     + + + COMMON BLOCKS + + +
*$noereference
      include 'p1werm.inc'
      include 'p1unconv.inc'

c     + + + LOCAL COMMON BLOCKS + + +
      include 'hydro/htheta.inc'
C ***      include 'hydro/hlayrs.inc'
      include 'hydro/hpsd.inc'
C ***      include 'hydro/htemp.inc'
*$reference

c     + + + LOCAL VARIABLES + + +
      integer k
      real potes(mnsz)

c     + + + LOCAL DEFINITION + + +
c     potes  - Air entry potential at a std. bsdblk of 1.3 Mg/m^3

c     + + + SUBROUTINES CALLED + + +
c     psd
c     extra

c     + + + FUNCTION DECLARATIONS + + +
      real waterk
      real calctht0

c     + + + END SPECIFICATIONS + + +
c     initialize various soil layer references

c     convert soil water contents from mass basis to volume basis

      do 105 k=1,layrsn
          theta(k) = bhrwc(k) * bsdblk(k)
C ***      write(*,*) ' bhrwc(k), bsdblk(k) ', k, bhrwc(k), bsdblk(k)p
          thetaf(k) = bhrwcf(k) * bsdblk(k)        ! field capacity
          thetaw(k) = bhrwcw(k) * bsdblk(k)        ! wilting point
C *** new way of calculating thetas per ffox
          thetas(k) = 1 - bsdblk(k) / 2.65
C ***          thetas(k) = bhrwcs(k) * bsdblk(k)        ! saturation
          thetes(k) = thetas(k) * 0.883            ! reduced saturation content
          thetar(k) = thetaw(k) * 0.369            ! realllllyyyy dry
          wc(k) = theta(k) * bszlyt(k)             ! layer water content
 105  continue

      swci = sum(wc(1:layrsn))

      theta(0) = calctht0(bszlyd, thetes, thetar, theta,
     *  thetaw, thetaf(1) - thetaw(1), 0.0_8, 0.0_8)                        !H-64,65,66

c      write(*,*) 'calctht0: hinit - theta(0) after',theta(0)

c     call subroutine psd and calculate soil hydraulic parameters
      do 110 k=1,layrsn
          if ( bh0cb(k) .eq. -99.9 ) then
              call psd( bsfsan(k), bsfsil(k), bsfcla(k), gmd(k),
     *          gsd(k) )
              potes(k) = -0.2 * gmd(k)**(-0.5)         !H-77
              bh0cb(k) = -2. * potes(k) + 0.2 * gsd(k)  !H-78
              if ( bheaep(k) .eq. -99.90 ) then
                 bheaep(k) = potes(k)*(bsdblk(k)/1.3)**(0.67*bh0cb(k))    !H-79
              end if
          end if
          if ( bhrsk(k) .eq. -99.90 ) then
             bhrsk(k) = waterk(bsdblk(k), bh0cb(k), bsfcla(k),
     *         bsfsil(k))
          end if
C ***       write (*,*) ' theta(k), etc. ', k, theta(k),thetas(k),-bh0cb(k)
110    continue

      return
      end
*
*