subroutine caqwet(smtill,smutil,qwet,bdtill,bdutil)
c
c     +++PURPOSE+++
c     This function calculates the heat of fusion released by freezing H2O
c     that is held in place or migrates to the freezing front based on
c     equation 3.1.4 from the WEPP manual August 1989.  The calculation
c     here is L * K w [P/Z uf].  This equation has been modified,
c     the Z uf terms cancel from the numerator, and denominator.  The
c     EXPO portion of the equation here is a multiplier because...
c       "When water in the soil freezes, the water in the unfrozen soil
c     below is depleted, therefore the hydraulic conductivity of the
c     soil is depleted" George Benoit, July 1992.
c
c     Author(s):  John Witte, UofMn WCES @ USDA-ARS-NCSRL
c     Date: 04/06/93

c     Verified and tested by Reza Savabi, USDA-ARS, NSERL 317-494-5051
c                  August 1994
c
c     +++ARGUMENT DECLARATIONS+++
      real    smtill,smutil
c
c     +++ARGUMENT DEFINITIONS+++
c     smtill - Soil moisture in the tilled layer (Dec %).
c     smutil - Soil moisture in the untilled layer (Dec %).
c     qwet   - Flux required to freeze water (W/m^2).
c     bdtill - Bulk density of the tilled soil layer (Kg/m^3).
c     bdutil - Bulk density of the untilled soil layer (Kg/m^3).
c
c     +++PARAMETERS+++
      include  'pmxhil.inc'
      include  'pmxpln.inc'
      include  'pmxnsl.inc'
c
c     +++COMMON BLOCKS+++
      include  'cwint.inc'
c       read:  amtfrz
c
      include  'cstruc.inc'
c       read:  iplane
c
      include  'cwater.inc'
c       read:  thetfc,st,dg
c
c     +++LOCAL VARIABLES+++
c
c XXX Use of the global SAVE does not follow the WEPP coding conventions
c     This needs to be fixed so that only the local variables which
c     need to have their values retained are saved.  dcf  5/18/94
      save
      real    numer,lhfh2o,hctill,hcutil,tentil,tenutl,denom,ufutdp,
     1        uftdep,expo,wflux,slmst,portil,porutl,slotil,sloutl,ten,
     2        pwater,bdtill,bdutil,qwet
c
c     +++LOCAL DEFINITIONS+++
c     pwater - Total soil water content in topmost soil layer (m/m).
c     numer  - Temporary variable for numerator part of 3.1.4.
c     lhfh2o - Latent heat of fusion of water (W/hr/m^3).
c     hctill - Hydraulic conductivity of water in tilled soil (m/hr).
c     hcutil - Hydraulic conductivity of water in untilld soil (m/hr).
c     tentil - Tension at the tilled layer (m).
c     tenutl - Tension at the untilled layer (m).
c     denom  - Temporary variable for denominator part of 3.1.4.
c     ufutdp - Depth of the unfrozen untilled layer (m).
c     uftdep - Depth of the unfrozen tilled layer (m).
c     expo   - Soil moisture adjustment (unitless).
c     wflux  - Temporary variable for water flux freeze (W/m^2).
c     portil - Porosity of the tilled soil layer.
c     porutl - Porosity of the untilled soil layer.
c     slotil - Slope of the retension curve - tilled.
c     sloutl - Slope of the retension curve - untilled.
c     slmst  - Soil moisture (%).
c     ten    - Applicable soil-water tension at the frost front (m).

c
c     +++DATA INITIALIZATIONS+++
      data   lhfh2o/9.3027e04/
c
c     +++END SPECIFICATIONS+++

      hctill = 1.0
      hcutil = 1.0
      numer  = 0.0
      denom  = 0.0
      expo   = 0.0
      wflux  = 0.0
      slmst  = 0.0
      qwet   = 0.0
      ufutdp = 0.0
      uftdep = 0.0

c -- If the entire tilled layer is frozen, the hctill value remains
c -- at 1.0 in order to be overlooked in the numerator part of the eq.

      pwater = (st(1,iplane) + thetdr(1,iplane) * dg(1,iplane)) /
     1                dg(1,iplane)

      portil = 1.0 - ((bdtill / 1000) / 2.65)
      porutl = 1.0 - ((bdutil / 1000) / 2.65)

      slotil = 100. * (thetfc(1,iplane) - portil) / 1.518514
      sloutl = 100. * (thetfc(1,iplane) - portil) / 1.518514

      tentil = (10.**((100.*(portil - pwater)) / abs(slotil))) / 100.
      tenutl = (10.**((100.*(porutl - pwater)) / abs(sloutl))) / 100.

      if (thdp(iplane) .gt. 0.0001) then
        if (tilld(iplane) .gt. 0.0001) then
          ten = tentil
        else
          ten = tenutl
        endif
      else
        if (frdp(iplane) .lt. tilld(iplane)) then
          ten = tentil
        else
          ten = tenutl
        endif
      endif

c -- If the entire tilled layer is frozen, the hctill value remains
c -- at 1.0 in order to be overlooked in the numerator part of the
c -- equation.

      if ((frdp(iplane) .lt. tilld(iplane)) .and.
     1                                 (tilld(iplane) .gt. 0.0001)) then
        uftdep = tilld(iplane) - frdp(iplane)
        ufutdp = (frdp(iplane) + 1.0) - tilld(iplane)
        hctill = usatkt(iplane)
        hcutil = usatku(iplane)
      else
        ufutdp = 1.0
        hcutil = usatku(iplane)
      endif

      numer = lhfh2o * hctill * hcutil * ten

c -- Note that in above eqn, tenabv-tenblw is equivalent to P in 3.1.4.

      denom = (hctill * ufutdp) + (hcutil * uftdep)

c -- Add C's version weighted calculation here...

      slmst = (smtill + smutil) / 2

      expo = exp(- amtfrz(iplane) / slmst)

c -- Turned this adjustment factor off.

      expo = 1.0

      if (denom .gt. 0.0) then
        wflux = (numer / denom) * expo
      else
        wflux = 0.0
      endif

      qwet = wflux

      return

      end