c$Author: fredfox $
c$Date: 2001-11-21 22:38:05 $
c$Revision: 1.1 $
c$Source: /weru/cvs/weps/weps.src/hydro/param_blkden_adj.for,v $

      subroutine param_blkden_adj( nlay, bsdblk, bsdblk0,
     &                         bsdpart, bhrwcf, bhrwcw, bhrwca,
     &                         bsfcla, bsfom,
     &                         bh0cb, bheaep, bhrsk )

c     + + + PURPOSE + + +
c     
c     This subroutine adjusts the air entry potential and saturated hydraulic
c     for changes in bulk density

c     + + + KEYWORDS + + +
c     bulk density adjustment hydraulic parameters

      include 'hydro/clayomprop.inc'

c     + + + ACCESSED COMMON BLOCK VARIABLE DEFINITIONS + + +
c     claygrav80rh - gravimetric water content of clay at 80 percent relative humidity
c     orggrav80rh  - gravimetric water content of soil organic matter at 80 percent relative humidity

c     + + + ARGUMENT DECLARATIONS + + +
      integer nlay
      real bsdblk(*), bsdblk0(*),
     &     bsdpart(*), bhrwcf(*), bhrwcw(*), bhrwca(*),
     &     bsfcla(*), bsfom(*),
     &     bh0cb(*), bheaep(*), bhrsk(*)

c     + + + ARGUMENT DEFINITIONS + + +
c     nlay     - number of soil layers to be updated
c     bsdblk   - bulk density (Mg/m^3)
c     bsdblk0  - previous day bulk density (Mg/m^3)
c     bsdpart  - particle density (Mg/m^3)
c     bhrwcf   - gravimetric 1/3 bar water
c     bhrwcw   - gravimetric 15 bar water
c     bhrwca   - gravimetric plant available water
c     bsfcla   - fraction of soil mineral portion which is clay
c     bsfomf   - fraction of total soil mass which is organic matter
c     bh0cb    - Brooks and Corey pore size interation exponent b
c     bheaep   - Brooks and Corey air entry potential
c     bhrsk    - Saturated hydraulic conductivity (m/s)

c     + + + FUNCTION DECLARATIONS + + +
      real volwatadsorb
      real volwat_matpot_bc

c     + + + LOCAL VARIABLES + + +
      integer lay
      real thetas, thetaf, thetaw, thetar
      real temp, temp1

c     + + + LOCAL VARIABLE DEFINITIONS + + +

c     + + + END SPECIFICATIONS + + + 

      do lay=1,nlay
c         adjust air entry potential from Campbell (1985) pg 46
          bheaep(lay) = bheaep(lay)
     &                * (bsdblk(lay)/bsdblk0(lay))**(0.67*bh0cb(lay))
c         adjust saturated hydraulic conductivity based on Campbell (1985) pg 54
          bhrsk(lay) = bhrsk(lay)
     &               * (bsdblk0(lay)/bsdblk(lay))**(1.34*bh0cb(lay))
c         update previous day bulk density since changes done
          bsdblk0(lay) = bsdblk(lay)

c             reverse calculation of field capacity and wilting point
          thetas = 1 - bsdblk(lay) / bsdpart(lay)   ! saturation

c         use theta corresponding to 80% relhum in soil for thetar
          temp = bsdblk(lay)*1000.0  !convert Mg/m^3 to kg/m^3
          thetar = volwatadsorb( temp, bsfcla(lay), bsfom(lay),
     &                           claygrav80rh, orggrav80rh )

          temp = -33.33
          temp1 = 1.0 / bh0cb(lay)
          thetaf = volwat_matpot_bc(temp, thetar, thetas,
     &                              bheaep(lay), temp1)
          temp = -1500.0
          thetaw = volwat_matpot_bc(temp, thetar, thetas,
     &                              bheaep(lay), temp1)

c         update gravimetric values for these properties
          bhrwcf(lay) = thetaf / bsdblk(lay)        ! field capacity
          bhrwcw(lay) = thetaw / bsdblk(lay)        ! wilting point
          bhrwca(lay) = bhrwcf(lay) - bhrwcw(lay)   ! plant available capacity
      end do

      end