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

      subroutine param_pot_bc( nlay, bsdblk,
     &                         bsdpart, bhrwcf, bhrwcw,
     &                         bsfcla, bsfom,
     &                         bh0cb, bheaep )

c     + + + PURPOSE + + +
c     
c     This subroutine calculates matric potential parameters from given
c     values of bulk density, gravimetric 1.3 bar water, 15 bar water and
c     clay and organic matter fraction

c     + + + KEYWORDS + + +
c     matric potential 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(*),
     &     bsdpart(*), bhrwcf(*), bhrwcw(*),
     &     bsfcla(*), bsfom(*),
     &     bh0cb(*), bheaep(*)

c     + + + ARGUMENT DEFINITIONS + + +
c     nlay     - number of soil layers to be updated
c     bsdblk   - 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     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     + + + FUNCTION DECLARATIONS + + +
      real volwatadsorb

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

c     + + + LOCAL VARIABLE DEFINITIONS + + +

c     + + + END SPECIFICATIONS + + + 

      do lay=1,nlay
          thetaf = bhrwcf(lay) * bsdblk(lay)        ! field capacity
          thetaw = bhrwcw(lay) * bsdblk(lay)        ! 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 )

c         error check and adjustments to prevent numerical problems
c         with curve fit
          if( thetas.le.thetaf ) then
              write(*,*) 'saturation less than field capacity, layer',
     &                   lay
              stop
          else if( thetaf.le.thetaw ) then
              write(*,*) 'field capacity less than wilting point, layer'
     &                  ,lay
              stop
          end if
          thetar = min( thetar, 0.8 * thetaw )

c         Calculate air entry and b to match saturated, field
c         capacity, permanent wilting point and residual water
c         as calcuated above 
          bh0cb(lay) = -(log(33.333)-log(1500.0))
     &            / (log((thetaf-thetar)/(thetas-thetar))
     &            - log((thetaw-thetar)/(thetas-thetar)))
          bheaep(lay)=-exp(log(1500.0)+bh0cb(lay)
     &             *log((thetaw-thetar)/(thetas-thetar)))

c         error check brooks and corey b value to keep in range and
c         prevent later numberical problems
          if( bh0cb(lay).gt.50.0 ) then
              write(*,*) 'Derived Brooks and Corey b too large, layer',
     &                   lay
              stop
          endif

      end do

      end