subroutine wmlt(kres,kutill,kuutil,resdp,hour) c c +++PURPOSE+++ c This function adjusts the thaw depth and frost depth in the soil c when there is no top frost present. It calculates the melt c taking place from the top-down based on heat coming down from c above. This function is called from sfrost.for when the hourly c surface temperature is > 0.0. c c Author(s): John Witte, UofMn WCES @ USDA-ARS-NCSRL c Date: 04/06/93 c c Verified and tested by Reza Savabi, USDA-ARS, NSERL 317-494-5051 c August 1994 c c NOTE: The logic of this module almost exactly parallels that c for TFWMLT. This may have implications for optimizing c their code. -- CRM -- 3/18/96. c c Inserted into main WEPP code 2/3/97 by dcf c c +++ARGUMENT DECLARATIONS+++ real kres,kutill,kuutil,resdp integer hour c c +++ARGUMENT DEFINITIONS+++ c kres - Thermal conductivity of the residue layer (W/meter C). c kutill - Thermal cond. of the unfrozen tilled layer (W/meter C). c kuutil - Thermal cond. of the unfrozen untilled layer (W/meter C). c resdp - Depth of the residue cover (m). c hour - Hour of the day being simulated (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: thdp(mxplan),wftill,pfwutl,surtmp c include 'cstruc.inc' c read: iplane include 'cupdate.inc' c read: sdate c +++LOCAL VARIABLES+++ save real drythw,lhfh2o,uftdep,tufutd,numer,denom,tctill,tcutil, 1 tcres,incr, olddep c c +++LOCAL DEFINITIONS+++ c drythw - Heat needed to thaw water in the soil (W/m^3). c lhfh2o - Latent heat of fusion of water (W-hr/m^3). c uftdep - Depth of unfrozen tilled layer (m). c tufutd - Temporary unfrozen untilled depth variable (m). c numer - Temporary variable used to represent numerator in eqn. c denom - Temporary variable used to represent denominator in eqn. c tctill - Thermal conductivity of the unfrozen tilled layer (W/meter C). c tcutil - Thermal cond. of the unfrozen untilled layer (W/meter C). c tcres - Thermal cond. of the residue layer (W/meter C). c incr - Amount of melt taking place this hour (m). c c +++DATA INITIALIZATIONS+++ data lhfh2o/9.3027e04/ c c +++END SPECIFICATIONS+++ c c tctill = 1.0 tcutil = 1.0 incr = 0.0 uftdep = 0.0 tufutd = 0.0 drythw = 0.0 c c c -- Cases 1,3, and 6. c c *** L1 - IF *** if (thdp(iplane) .lt. 0.0001) then thdp(iplane) = 0.0 olddep = thdp(iplane) c c -- Cases 1 and 3. c if (tilld(iplane) .gt. 0.0) then drythw = lhfh2o * pftill(iplane) else c c -- Case 6. c drythw = lhfh2o * pfwutl(iplane) endif c c -- Cases 2,4,5, and 7. c c *** L1 - ELSE *** else olddep = thdp(iplane) c c -- Cases 2 and 4. c if (thdp(iplane) .lt. tilld(iplane)) then uftdep = thdp(iplane) drythw = lhfh2o * pftill(iplane) c c -- Cases 5 and 7. c else uftdep = tilld(iplane) tcutil = kuutil tufutd = thdp(iplane) - tilld(iplane) drythw = lhfh2o * pfwutl(iplane) endif tctill = kutill c c *** L1 - ENDIF *** endif c c c -- Here, we set the value for thermal cond. of residue if present. if (resdp .gt. 0.0001) then tcres = 0.168 else tcres = 1.0 endif c c -- Adjust the thaw depth... c denom = (tcres * (tctill * tufutd + uftdep * tcutil) + 1 resdp * tctill * tcutil) * drythw c if (denom .gt. 0.0) then numer = tcres * tctill * tcutil * surtmp(hour) incr = (numer / denom) c limit the increment to the maximum value of 0.01 - Savabi incr = min(incr,0.01) else incr = 0.0 endif c c ----- depth of thaw from soil surface thdp(iplane) = thdp(iplane) + incr * 0.90909 c c -- NOTE: The reason for 1.1 above and .0909 below is due to c -- ==== the expansion coefficient we used in freezing. We c == now must account for it again. c c ----- Don't let thawing progress below the tilled layer. if ((olddep .lt. tilld(iplane)) .and. 1 (thdp(iplane) .gt. tilld(iplane))) then thdp(iplane) = tilld(iplane) incr = thdp(iplane) - olddep endif c c -- As noted above, the 0.0909 below is due to the expansion c -- coefficient we have used, it is (1.0 - 0.909 = 0.0909). c c ----- depth of freezing (frost) from soil surface frdp(iplane) = frdp(iplane) - (0.0909 * incr) if(frdp(iplane) .lt. 0.0001) frdp(iplane) = 0.0 c c -- Below we see cases 1 - 4. c c *** M1 - ENDIF *** if (thdp(iplane) .lt. tilld(iplane)) then c c -- Cases 1 and 2. c if (frdp(iplane) .lt. tilld(iplane)) then c wftill(iplane) = (frdp(iplane) - thdp(iplane)) * 1 pftill(iplane) wfutil(iplane) = 0.0 c c -- Cases 3 and 4. c else c wftill(iplane) = (tilld(iplane) - thdp(iplane)) * 1 pftill(iplane) wfutil(iplane) = (frdp(iplane) - tilld(iplane)) * 1 pfwutl(iplane) endif c c -- Below we see cases 5,6, and 7. c c *** M1 - ELSE *** else wftill(iplane) = 0.0 wfutil(iplane) = (frdp(iplane) - thdp(iplane)) * 1 pfwutl(iplane) c c *** M1 - ENDIF *** endif c c *** N1 - IF *** if (thdp(iplane) .lt. frdp(iplane)) then if(thdp(iplane).lt.tilld(iplane)) then uftild(iplane) = thdp(iplane) else uftild(iplane) = tilld(iplane) endif c c *** N1 - ELSE *** else frdp(iplane) = 0.0 thdp(iplane) = 0.0 amtfrz(iplane) = 0.0 wftill(iplane) = 0.0 wfutil(iplane) = 0.0 pftill(iplane) = 0.0 pfwutl(iplane) = 0.0 uftild(iplane) = tilld(iplane) c c *** N1 - ENDIF *** endif c c c ----- Increase freeze/thaw cycle (Also in TFWMLT.) if ((olddep .lt. 0.0001).and.(thdp(iplane).gt.0.0001)) 1 cycle(iplane) = cycle(iplane) + 1 c amtfrz(iplane) = wftill(iplane) + wfutil(iplane) c return end