c     file name: extra.for

      real function extra (bszlyd, thetamx, thetamn, theta)

c     + + + PURPOSE + + +
c     This subroutine extrapolates soil water content to the surface
c     from the three uppermost simulation layers.  A numerical
c     solution known as Cramer's rule is used to obtain an estimate
c     of the extrapolated surface soil water content by solving the
c     three simultaneous equations that describe the relationship
c     between soil water content and soil depth for the three
c     uppermost simulation layers.
c     DATE:  09/22/93
c     MODIFIED:  10/06/93

c     + + + KEY WORDS + + +
c     soil, water content

c     + + + ARGUMENT DECLARATIONS + + +
      real bszlyd(*)
      real thetamn, thetamx
      real theta(0:*)
      real thetax

c     + + + ARGUMENT DEFINITIONS + + +
c     bszlyd  - Depth to bottom of soil layer from surface (mm)
c     thetax - Extrapolated surface soil water content (m^3/m^3)

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

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

c     + + + LOCAL VARIABLES + + +
      real d
      real d1

c     + + + LOCAL DEFINITIONS + + +
c     d      - The determinant of the coefficient matrix.
c     d1     - The determinant of the matrix formed by substituting
c              load vector into column 1 of the coefficient matrix.

c     + + + END SPECIFICATIONS + + +


      d = (bszlyd(2)*bszlyd(3)**2) + (bszlyd(3)*bszlyd(1)**2) +
     &    (bszlyd(1)*bszlyd(2)**2) - (bszlyd(2)*bszlyd(1)**2) -
     &    (bszlyd(1)*bszlyd(3)**2) - (bszlyd(3)*bszlyd(2)**2)

      d1 = (theta(1)*bszlyd(2)*bszlyd(3)**2) +
     &     (theta(2)*bszlyd(3)*bszlyd(1)**2) +
     &     (theta(3)*bszlyd(1)*bszlyd(2)**2) -
     &     (theta(3)*bszlyd(2)*bszlyd(1)**2) -
     &     (theta(2)*bszlyd(1)*bszlyd(3)**2) -
     &     (theta(1)*bszlyd(3)*bszlyd(2)**2)

C Check to make sure that "d" is not too close to zero (thetax gets big)
      if (d .lt. 0.0000001) then
        thetax = thetamx
      else 
        thetax= d1/d

        if ( thetax .gt. thetamx )   thetax= thetamx
        if ( thetax .le. 0.0 )         thetax= thetamn
      endif 

      extra = thetax
      end