c$Author: fredfox $
c$Date: 2001-09-27 22:04:42 $
c$Revision: 1.3 $
c$Source: /weru/cvs/weps/weps.src/util/misc/setbds.for,v $

      real function setbds (clay, sand, om)

c     + + + PURPOSE + + +
c     The following function estimates settled soil bulk density from
c     intrinsic properties. see Rawls (1983) Soil Science 135, 123-125.

c     Should eventually be called by MAIN to initialize the values
c     for each subregion (unless soil composition changes).

c     + + + KEYWORDS + + +
c     bulk density, initialization

c     + + + PARAMETERS AND COMMON BLOCKS + + +
c     none

c     + + + ARGUMENT DECLARATIONS + + +
      real clay, sand, om

c     + + + ARGUMENT DEFINITIONS + + +

c     setbd - settled bulk density
c     clay  - fraction of soil clay content
c     sand  - fraction of soil sand content
c     om    - organic matter

c     + + + LOCAL VARIABLES + + +
c
      integer li, lj, hi, hj
      real mi, mj, fi, fj
      real mbdtv (0:10,0:10), mbd, tempsum
      real mbd_hi_hj

c     li  - index into table less than sand content
c     lj  - index into table less than clay content
c     hi  - index into table higher than sand content
c     hj  - index into table higher than clay content
c     mi  - value between indexes for interpolation for sand
c     mj  - value between indexes for interpolation for clay
c     fi  - fraction of distance between grid cells for sand
c     fj  - fraction of distance between grid cells for clay
c     mbdtv (0:10,0:10) - data table of settled bulk density
c                         as a function of sand (across the top)
c                         and clay (down the side)
c     mbd   - mineral bulk density without organic matter
c     mbd_hi_hj - value for mbdtv(hi,hj), if outside triangular
c                 part of table it is reflected from mbdtv(li,lj)
c                 otherwise it is just the real point

c     + + + SUBROUTINES CALLED + + +

c     + + + FUNCTION DECLARATONS + + +

c     + + + DATA INITIALIZATIONS + + +
c     first index in this direction ->
c     second index || goes down 
c                  \/
      data mbdtv /1.48,1.25,1.00,1.06,1.16,1.22,1.30,1.39,1.45,1.51,1.52
     &           ,1.52,1.40,1.19,1.25,1.32,1.40,1.52,1.58,1.63,1.65,0.
     &           ,1.52,1.40,1.25,1.35,1.45,1.53,1.60,1.64,1.72,0.,0.
     &           ,1.52,1.40,1.29,1.41,1.50,1.57,1.63,1.68,0.,0.,0.
     &           ,1.50,1.40,1.35,1.43,1.53,1.61,1.64,0.,0.,0.,0.
     &           ,1.46,1.40,1.40,1.43,1.53,1.62,0.,0.,0.,0.,0.
     &           ,1.45,1.40,1.38,1.42,1.50,0.,0.,0.,0.,0.,0.
     &           ,1.42,1.37,1.33,1.33,0.,0.,0.,0.,0.,0.,0.
     &           ,1.33,1.32,1.20,0.,0.,0.,0.,0.,0.,0.,0.
     &           ,1.23,1.18,0.,0.,0.,0.,0.,0.,0.,0.,0.
     &           ,1.15,0.,0.,0.,0.,0.,0.,0.,0.,0.,0./

c     + + + END SPECIFICATIONS + + +

      tempsum = sand + clay
      if( tempsum.gt.1.0 ) then
          sand = sand / tempsum
          clay = clay / tempsum
          write(*,*) "setbds: sand plus clay fractions greater than 1.0"
          write(*,*) "values adjusted by averaging the difference"
      endif

c      i = nint(sand*100.0/10.0)
c      j = nint(clay*100.0/10.0)

c      mbd = mbdtv(i,j)

      mi = sand*10.0
      li = int(mi)
      fi = mi - li
      hi = min( li+1, 10 )
       
      mj = clay*10.0
      lj = int(mj)
      fj = mj - lj
      hj = min( lj+1, 10 )

c     check for table edge
      if( li + lj .eq. 10 ) then
c         on table edge, no interpolation necessary
          mbd = mbdtv(li,lj)
      else
          if( hi + hj .gt. 10 ) then
c             interpolation on the triangular edge of the table
c             mirror li,lj value to make using grid interpolation possible
              mbd_hi_hj =  mbdtv(li,hj) + mbdtv(hi,lj) - mbdtv(li,lj)
          else
c             interpolation within the table, use grid point
              mbd_hi_hj = mbdtv(hi,hj)
          end if
          mbd = (1-fi) * (1-fj) * mbdtv(li,lj)
     &        + (1-fi) * fj * mbdtv(li,hj)
     &        + fi * (1-fj) * mbdtv(hi,lj)
     &        + fi * fj * mbd_hi_hj
      end if

      setbds = 1.0 / ((om / 0.224)  + (1.0 - om)/ mbd)

      return
      end