subroutine den(
     *  csdblk, csdsbk, cszlyt, csdagd,
     *  dcump, chzwid)

c     + + + ARGUMENT DECLARATIONS + + +

      real csdblk, csdsbk, cszlyt, csdagd
      real dcump, chzwid

c     + + + LOCAL VARIABLES + + + 

      real bsdbk0
      integer j, nj
      real wsdblk, dsdblk
      real wszlyt

c     + + + LOCAL DEFINITIONS + + +

c   bsdbk0    - bulk density prior to update by SOIL, Mg/m^3
C   wsdblk    - bulk density of wet soil
C   dsdblk    - bulk density of dry soil
C   wszlyt    - depth of wetness in this layer

c  DENSITY SECTION:
c     update bulk density

c     store initial value of layer density
      bsdbk0 = csdblk

c     daily update density for other forces -- long term
c     only if current bulk density is less than settled bd
      if (csdblk.lt.csdsbk) then
        dsdblk = csdblk + 0.01*(csdsbk - csdblk)
      else
        dsdblk = csdblk
      endif

C if water has infiltrated into the layer
      if (chzwid .gt. 0) then
c     update for water additions in 5 mm increments
        nj =nint(dcump/5.)
        wsdblk = csdblk
        do j = 1,nj
          if (wsdblk .lt. 0.97 * csdsbk) then
            wsdblk = wsdblk+0.75 * (1-(wsdblk/
     &        (0.97 * csdsbk)))**1.5
          else
            exit
          endif
        enddo
      else 
C set wsdblk to something so that it won't bomb out as uninit'd later on
        wsdblk = dsdblk
      endif

C get weighted average of wet and dry densities
      wszlyt = min(cszlyt, chzwid)
      csdblk = (wsdblk*wszlyt + dsdblk*(cszlyt-wszlyt)) /
     *  cszlyt

C reduce wet depth by this layer's thickness or wet depth
      chzwid = chzwid - wszlyt
   
c     update layer thickness
      cszlyt = cszlyt * bsdbk0 / csdblk
      
c     update aggregate density
      csdagd = csdsbk
c
      end