subroutine aggsta(
     *  cseags, cseagmn, cseagmx,
     *  chrwc0, cbhrwc0, chrwcw, chrwca,
     *  chrwc1, cbhrwc,
     *  chtmx0, chtsmn, chtsmx, ck4d,
     *  se0, k4f, se1, k4w, se)

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

c     + + + LOCAL VARIABLES + + + 

c     + + + LOCAL DEFINITIONS + + +

c  AGGREGATE STABILITY SECTION:

c        calc. relative agg stability and water content for prior day
         se0 = (cseags - cseagmn)/(cseagmx - cseagmn)
         chrwc0 = (cbhrwc0 - chrwcw)/chrwca
         if (chrwc0.lt.0.0) chrwc0 = 0.00001

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

c        check soil temp to determine process

c        check if unfrozen then wet or dry
         if ((chtmx0 .gt. 0.0) .and. (chtsmn .gt. 0.0)) 
     *     go to 70

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

c        check for thaw process alone
         if ((chtmx0 .lt. 0.0) .and. (chtsmx .ge. 0.0)
     &      .and. (chtsmn .ge. 0.0)) go to 60

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

c        check for drying while frozen
         elseif ((chtsmn .lt. 0.0) .and. (chtmx0.lt.0)
     &         .and. (cbhrwc .lt. cbhrwc0)) then
c            update stability with drying while frozen
             se1 = se0 - (2.7 - se0) * k4f * (chrwc0 - chrwc1)/
     &             (2.5 - k4f * chrwc1)
             go to 80
         else
             go to 80
         endif

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

c     thaw process with prior day water content
   60    continue
         if (chrwc0 .gt. 1.0) then
c     soil puddling process when saturated
            se0 = 0.74/chrwc0
         else
            se = se0 - (2.7 - se0)*k4f*chrwc0/
     *        (2.5 - k4f * chrwc0)
            se0 = se*(1.0 - ck4d*chrwc0)/
     *        (1.0 - ck4d*k4f*chrwc0)
        endif
        go to 70

c     check if wetting process
   70   if (cbhrwc .gt. cbhrwc0) then
c     wetting process
C ***             write(*,*) 'soil a2: se0, k4w, chrwc0, chrwc1 ',
C ***     *         se0, k4w, chrwc0, chrwc1
           se1 = se0*(1.0 - k4w*chrwc1)/(1.0-k4w*chrwc0)
        else
c     drying process
C ***             write(*,*) 'soil a3: se0, k4d, chrwc0, chrwc1 ',
C ***     *         se0, k4d, chrwc0, chrwc1
           se1 = se0*(1.0 - ck4d*chrwc1)/
     *       (1.0 - ck4d*chrwc0)
        endif

c     check if freeze after w/d
        if ((chtmx0 .gt. 0.0) .and. (chtsmx .lt. 0.0)
     &     .and. (chtsmn .lt. 0.0)) then
c     freeze process with today water content
           se = se0*(1.0 - k4w*k4f*chrwc1)/(1.0 - k4w*
     &          chrwc1)
C ***             write(*,*) 'soil a4: se, k4f, chrwc1 ',
C ***     *         se, k4f, chrwc1
           se1 = se + ((2.7 - se)/2.5)*k4f*chrwc1
        endif

c     calc. new agg. stability                                     (S-28)
   80   continue
C ***       write(*,*) 'soil: se1, cseagmx, cseagmn ',
C ***     *   se1, cseagmx, cseagmn
       if (se1.lt.0.0) then
c         write(*,*) ' debugging fix se1'
         se1 = 0.0
       endif

        cseags = se1*(cseagmx - cseagmn) + cseagmn
C
C if not frozen, don't allow over max
        if ((chtsmx.gt.0.0) .and. (cseags.gt.cseagmx))
     *    cseags = cseagmx
C
C do not allow to go under minimum
        if (cseags .lt. cseagmn) cseags = cseagmn
      end