c$Author: fredfox $
c$Date: 2001-09-27 22:02:44 $
c$Revision: 1.5 $
c$Source: /weru/cvs/weps/weps.src/soil/aggsta.for,v $
      subroutine aggsta(
     *  cseags, cseagmn, cseagmx,
     *  cbhrwc0, chrwcw, chrwca,
     *  cbhrwc,
     *  chtmx0, chtsmn, chtsmx, ck4d,
     *  se0, k4f, se1, k4w)

c     + + + ARGUMENT DECLARATIONS + + +
      real
     *  cseags, cseagmn, cseagmx,
     *  cbhrwc0, chrwcw, chrwca,
     *  cbhrwc,
     *  chtmx0, chtsmn, chtsmx, ck4d,
     *  se0, k4f, se1, k4w

c     + + + LOCAL VARIABLES + + +
      real se, minse, maxse, temp, hrwc0, hrwc1
      parameter (minse = 0.01)
      parameter (maxse = 1.0)

c     + + + LOCAL DEFINITIONS + + +
c     se        - relative aggregate stability with partial update
c     minse     - minimum value allowed for se
c     maxse     - maximum value allowed for se
c     temp      - temporary variable
c     hrwc0     - relative water content on prior day of each layer
c     hrwc1     - relative water content on current day of each layer

c  AGGREGATE STABILITY SECTION:

c     relative agg stability for prior day
      se0 = (cseags - cseagmn)/(cseagmx - cseagmn)

c     relative water content for prior day
      hrwc0 = (cbhrwc0 - chrwcw)/chrwca
      if (hrwc0.lt.0.0) hrwc0 = 0.0

c     relative water content for current day
      hrwc1 = (cbhrwc - chrwcw) / chrwca
      if (hrwc1 .lt. 0.0) hrwc1 = 0.0

c     check for freeze after wet or dry
      if((chtmx0.gt.0.0).and.(chtsmn.lt.0.0)
     &                  .and.(chtsmx.lt.0.0)) then
          go to 70

c     check for thaw process alone (chtsmn .gt. or .lt 0)
      else if((chtmx0.lt.0.0).and.(chtsmx.ge.0.0)) then
          go to 60

c     check for freeze/thaw
      else if((chtmx0 .gt. 0.0).and.(chtsmn.lt.0.0)
     &                         .and.(chtsmx.gt.0.0)) then
          goto 50

c     check for change in water content while frozen
c     all the freeze/thawing cases were caught above
      else if((chtmx0.lt.0).and.(chtsmn.lt.0.0)
     &                     .and.(chtsmx.lt.0.0)) then
          if (hrwc1 .lt. hrwc0) then
c             update stability with drying while frozen
              se1 = se0 - (2.7 - se0) * k4f * (hrwc0 - hrwc1)/
     &             (2.5 - k4f * hrwc1)
          else
c             update stability with wetting while frozen
              se1 = se0
c               write(*,*) 'aggsta: all frozen wc change',hrwc0,hrwc1
          endif
          go to 80
      else
c        temperature is all above
         go to 70
      endif

c     freeze process with prior day water content
   50 se = se0 * (1.0 - k4w * k4f * hrwc0)/
     *     (1.0 - k4w * hrwc0)
      se0 = se + ((2.7 - se)/2.5)*k4f*hrwc0

c     thaw process with prior day water content
   60 if (hrwc0 .gt. 1.0) then   !soil puddling process when saturated
            se0 = max( minse, 0.999-ck4d*hrwc0)
      else                        !thaw then wetting
            se = se0 - (2.7 - se0)*k4f*hrwc0/
     *        (2.5 - k4f * hrwc0)
            se0 = se*(1.0 - ck4d*hrwc0)/
     *        (1.0 - ck4d*k4f*hrwc0)
      endif

c     check for wetting or drying process
   70 if (hrwc1 .gt. hrwc0) then      !wetting process
           se1 = se0*(1.0 - k4w*hrwc1)/(1.0-k4w*hrwc0)
      else                               !drying process
           se1 = se0*(1.0 - ck4d*hrwc1)/
     *       (1.0 - ck4d*hrwc0)
      endif

c     check if freeze after w/d (chtsmn must be .lt. 0.0)
      if ((chtmx0 .gt. 0.0) .and. (chtsmx .lt. 0.0)) then
c          freeze process with today water content
           se = se0*(1.0 - k4w*k4f*hrwc1)/(1.0 - k4w*
     &          hrwc1)
           se1 = se + ((2.7 - se)/2.5)*k4f*hrwc1
      endif

c     calc. new agg. stability                                   (S-28)
   80 if (se1.lt.minse) then
         se1 = minse
      endif
C     if not frozen, don't allow over max
      if ((se1.gt.maxse).and.(chtsmx.gt.0.0)) then
c         write(*,*) "aggsta: se1=",se1," Will be set to",maxse
         se1 = maxse
      endif

c     set resulting aggregate statbility based on range limited se1
      cseags = se1*(cseagmx - cseagmn) + cseagmn

      return      
      end