c$Author: fredfox $
c$Date: 2001-09-27 20:37:45 $
c$Revision: 1.2 $
c$Source: /weru/cvs/weps/weps.src/hydro/matricpot_bc.for,v $
      subroutine matricpot_bc(theta, thetar, thetas, airentry, lambda,
     &                        thetaw, theta80rh, soiltemp,
     &                        matricpot, soilrh )
c     returns: matricpot, soilrh
c     returns the matric potential in meters of water as defined by the 
c     Brooks and Corey function down to wilting point. Below wilting point,
c     the calculation is done based on clay and organic matter adsorption
c     isotherms. Coincidentally, the soil relative humidity is used to find 
c     the matric potential from the clay isotherms and is also returned for
c     potentials in the wetter range.

c*** Argument declarations ***
      real*4  theta, thetar, thetas, airentry, lambda,
     &        thetaw, theta80rh, soiltemp,
     &        matricpot, soilrh
c     theta      - present volumetric water content
c     thetar     - volumetric water content where hydraulic conductivity becomes zero
c     thetas     - saturated volumetric water content
c     airentry   - Brooks/Corey air entry potential (1/pressure) (modify to set units returned)
c     lambda     - Brooks/Corey pore size interaction parameter 
c     thetaw     - volumetric water content at wilting (15 bar or 1.5 MPa)
c     theta80rh  - volumetric water content at %80 relative humidity (300 bar or 30 MPa)
c     soiltemp   - soil temperature (C)
c     matricpot  - matric potential (meters of water)
c     soilrh     - relative humidity of soil air (fraction)

c*** function declarations ***
      real*4 soilrelhum

c*** Include files ***
      include 'hydro/vapprop.inc'

c*** Local variable declarations ***
      real*4  satrat
c     satrat     - conductivity relative saturation ratio

      if( theta .ge. thetaw ) then
          satrat = (theta-thetar)/(thetas-thetar)
          if( satrat .le. 0.0 ) then
              write(*,*) 'matricpot_bc: thetar= ',thetar,
     &                   ' thetaw= ', thetaw
              write(*,*) 'residual water content is greater than wilting
     & point'
              stop
          else if( satrat .ge. 1.0 ) then
              matricpot = airentry
          else
              matricpot = airentry*satrat**(-1.0/lambda)
          end if
          soilrh = soilrelhum(theta, thetaw, theta80rh, soiltemp,
     &                        matricpot)
      else
          soilrh = soilrelhum(theta, thetaw, theta80rh, soiltemp,
     &                        matricpot)
          matricpot = rgas*(soiltemp+zerokelvin)*(log(soilrh))
     &                 / (molewater * gravconst)
      end if
      return
      end