subroutine caqout(ksnow,kres,kftill,kfutil,resdp,qout,hour)
c     subroutine caqout(tfrstq,smtill,smutil,ksnow,kres,kftill,kfutil,
c    1                  resdp,qout,hour)
c
c     +++PURPOSE+++
c     This function calculates the layer heat flux based on
c     equation 3.1.1 from the WEPP manual, August 1989.  The actual
c     calculation used here is a reduction of 3.1.1.  When comparing
c     equations 3.1.1 and 3.1.2, the term Zsrf is ultimately found in
c     the numerator and denominator of the calculation and is therefore
c     cancelled out when doing the actual calculation here.
c     Layer heat flux refers to the heat passing down from above
c     through the frozen layer system such as snow, residue, and
c     frozen soil if such layers do exist.
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    ksnow,kres,kftill,kfutil,resdp,qout
      integer hour
c
c     +++ARGUMENT DEFINITIONS+++
c     ksnow  - Thermal Conductivity of the snow layer (W/min C).
c     kres   - Thermal Conductivity of the residue layer (W/min C).
c     kftill - Thermal Conductivity of the frozen tilled layer (W/min C).
c     kfutil - Thermal Conductivity of unfrozen untilled layer (W/min C).
c     resdp  - Depth of the residue layer (m).
c     qout   - Layer heat flux value being calculated here (W/m^2).
c     hour   - Hour of the day being run (hr).
c
c     +++PARAMETERS+++
      include 'pmxtil.inc'
      include 'pmxtls.inc'
      include 'pmxpln.inc'
      include 'pmxhil.inc'
      include 'pmxnsl.inc'
c
c     +++COMMON BLOCKS+++
c
      include  'cwint.inc'
c       read:  snodpt(mxplan),tfrdp(mxplan),frdp(mxplan),surtmp
c
      include  'cstruc.inc'
c       read:  iplane
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,denom,stbldp,hflux,tftdep,tfudep,
     1              tmp,tctill,tcutil,tcsnow,tcres

c
c     +++LOCAL DEFINITIONS+++
c     numer  - Numerator of the reduced fraction.
c     denom  - Denominator of the reduced fraction.
c     stbldp - Depth down to stable soil temperature (m).
c     hflux  - Heat flux being calculated (W/m^2).
c     tftdep - Temporary frozen tilled-layer depth (m).
c     tfudep - Temporary frozen untilled-layer depth (m).
c     tmp    - Temporary variable (deg C).
c     tctill - Thermal cond. of the frozen, tilled soil layer (W/m C).
c     tcutil - Thermal cond. of the frozen, untilled soil layer (W/m C).
c     tcsnow - Thermal cond. of snowpack layer (W/m C).
c     tcres  - Thermal cond. of the residue layer (W/m C).
c
c     +++END SPECIFICATIONS+++

      stbldp = 1.0
c
c -- If the given layers do not exist, then therm cond values are 1.
      tctill = 1.0
      tcutil = 1.0
      tcsnow = 1.0
      tcres  = 1.0
      qout   = 0.0
      numer  = 0.0
      denom  = 0.0
      tftdep = 0.0
      tfudep = 0.0

      if (snodpt(iplane) .gt. 0.0001) then
        tcsnow = ksnow
      endif

      if (resdp .gt. 0.0001) then
        tcres = kres
      endif

c -- If a thawed layer exists...

c
c XXX Changed "gt" to "ge" to be consistent with other limits.  dcf 6/94
c     if (thdp(iplane) .gt. 0.0001) then
      if (thdp(iplane) .ge. 0.0001) then

c -- If a top frost layer exists as well...

c
c XXX   Changed "gt" to "ge" to be consistent with other limits.dcf 6/94
c       if (tfrdp(iplane) .gt. 0.0001) then
        if (tfrdp(iplane) .ge. 0.0001) then

c -- If there is frozen tilled soil present in the system...

          if (tilld(iplane) .gt. 0.0) then
            if (tilld(iplane) .ge. tfrdp(iplane)) then
              tctill = kftill
              tftdep = tfrdp(iplane)
            else
               tctill = kftill
c              tcutl = kfutil
               tftdep = tilld(iplane)
               tfudep = tfrdp(iplane) - tilld(iplane)
            endif
          else
            tcutil = kfutil
            tfudep = tfrdp(iplane)
          endif

        endif

      endif

c -- If no top frost layer exists...

      if (thdp(iplane) .lt. 0.0001) then
c
c       Reset value of THDP to zero.  dcf 6/94
c
        thdp(iplane) = 0.0

c -- If a frost layer exists...

        if (frdp(iplane) .gt. 0.0) then

c -- If no till layer exists...

          if (tilld(iplane) .le. 0.0001) then
            tcutil = kfutil
            tfudep = frdp(iplane)

c -- Otherwise, if there is a tilled layer present...

          else
            if (tilld(iplane) .ge. frdp(iplane)) then
              tctill = kftill
c             tctdep = frdp(iplane)
            else
              tctill = kftill
              tcutil = kfutil
              tftdep = tilld(iplane)
              tfudep = frdp(iplane) - tilld(iplane)
            endif

          endif

        endif

      endif

c -- Now we can get into the calculations...

      tmp = abs(surtmp(hour))

      numer = tcsnow * tcres * tctill * tcutil * tmp

c -- NOTE: This calculation is only made when the adjusted
c -- ===== surface temperature is below zero, this is why
c --       the sign's are switched for adjusted temperature.
c --       See WEPP calculation 3.1.3 for avg thermal cond.

c -- John, check this out...

      denom = ((tcsnow * tcres * tctill * tfudep)   +
     1        (tcsnow * tcres * tftdep * tcutil)   +
     2        (tcsnow * resdp * tctill * tcutil)   +
     3        (snodpt(iplane) * tcres * tctill * tcutil)) * stbldp

      if (denom .gt. 0.0) then
        hflux = numer / denom
      else
        hflux = 0.0
      endif

      qout = hflux

      return

      end