subroutine den(
     *  csdblk, csdsblk, csdwblk, cszlyt, csdagd,
     *  chrwc0, chrwc, chrwca, chrwcw,
     *  dcump, chzwid)

c     + + + ARGUMENT DECLARATIONS + + +

      real csdblk, csdsblk, csdwblk, cszlyt, csdagd
      real chrwc0, chrwc, chrwca, chrwcw
      real dcump, chzwid

c     + + + ARGUMENT DEFINITIONS + + +

c     csdblk - present soil bulk density (Mg/m^3)
c     csdsblk - Settled soil bulk density (Mg/m^3)
c     csdwblk - 1/3 bar soil bulk density (Mg/m^3)
c     cszlyt - Soil layer thickness (mm)
c     csdagd - Soil aggregate density (Mg/m^3)
c     chrwc0 - Soil water content from previous day (g/g)
c     chrwc  - Soil water content on present day (g/g)
c     chrwca - Available soil water content (g/g)
c     chrwcw - Wilting point soil water content (g/g)
c     dcump  - daily water application depth (mm)
c     chzwid - water infiltration depth (mm)

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.csdsblk) then
        dsdblk = csdblk + 0.01*(csdsblk - 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 * csdsblk) then
            wsdblk = wsdblk+0.75 * (1-(wsdblk/
     &        (0.97 * csdsblk)))**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 = csdsblk
c
      end