subroutine tfwmlt(kutill,kuutil,resdp,hour) c c +++PURPOSE+++ c This function ajusts the frost and thaw layers when the c adjusted temperature is above freezing and there is a c top frost present. c c NOTE: Top frost melts from top only. Surface temp > 0.0 c ==== 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 NOTE: The logic of this module almost exactly parallels that c for WMLT. This may have implications for optimizing c their code. -- CRM -- 3/18/96. c c Inserted into WEPP model code - 2/3/97 - dcf c c +++ARGUMENT DECLARATIONS+++ real kutill,kuutil,resdp integer hour c c +++ARGUMENT DEFINITIONS+++ c kutill - Thermal cond. of the unfrozen tilled soil 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 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: tthawd(mxplan),wftill,pfwutl,surtmp c include 'cstruc.inc' c read: iplane c c +++LOCAL VARIABLES+++ save real drythw,lhfh2o,uftdep,tufutd,numer,denom,olddep,incr, 1 tctill,tcutil,tcres c c +++LOCAL DEFINITIONS+++ c drythw - Heat needed to thaw water in the soil (W/m^2). 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 olddep - Original top-thaw depth going into routine (m). c incr - Amount calculated top-thaw increment (m). c tctill - Thermal conductivity of the unfrozen tilled layer (W/meter C). c tcutil - Thermal conduct. of the unfrozen untilled layer (W/meter C). c tcres - Thermal conduct. of the residue layer @ soil surface (W/m C). c c +++DATA INITIALIZATIONS+++ data lhfh2o/9.3027e04/ c c +++END SPECIFICATIONS+++ c c -- Variable initialization portion of our program. 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 (tthawd(iplane) .lt. 0.0001) then tthawd(iplane) = 0.0 olddep = tthawd(iplane) c c -- Cases 1 and 3. c if (tilld(iplane) .gt. 0.0) then drythw = lhfh2o * pfttil(iplane) else c c -- Case 6. c drythw = lhfh2o * pftutl(iplane) endif c c -- Cases 2,4,5, and 7. c c *** L1 - ELSE *** else olddep = tthawd(iplane) c c -- Cases 2 and 4. c if (tthawd(iplane) .lt. tilld(iplane)) then uftdep = tthawd(iplane) drythw = lhfh2o * pfttil(iplane) c c -- Cases 5 and 7. c else uftdep = tilld(iplane) tcutil = kuutil tufutd = tthawd(iplane) - tilld(iplane) drythw = lhfh2o * pftutl(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 top 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) else incr = 0.0 endif incr = min(incr,0.01) c tthawd(iplane) = tthawd(iplane) + incr * 0.90909 c c tthawd(iplane) = amax1(tthawd(iplane), 0.0) if(tthawd(iplane) .lt. 0.0001) tthawd(iplane) = 0.0 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 if ((olddep .lt. tilld(iplane)) .and. 1 (tthawd(iplane) .gt. tilld(iplane))) then tthawd(iplane) = tilld(iplane) incr = tilld(iplane) - olddep endif c tfrdp(iplane) = tfrdp(iplane) - (0.0909 * incr) c c tfrdp(iplane) = amax1(tfrdp(iplane), 0.0) if(tfrdp(iplane) .lt. 0.0001) tfrdp(iplane) = 0.0 c c -- Below we see cases 1 - 4. c c *** M1 - IF *** if (tthawd(iplane) .lt. tilld(iplane)) then c c -- Cases 1 and 2. c if (tfrdp(iplane) .lt. tilld(iplane)) then wfttil(iplane) = (tfrdp(iplane) - tthawd(iplane)) * 1 pfttil(iplane) wftutl(iplane) = 0.0 c c -- Cases 3 and 4. c else wfttil(iplane) = (tilld(iplane) - tthawd(iplane)) * 1 pfttil(iplane) wftutl(iplane) = (tfrdp(iplane) - tilld(iplane)) * 1 pftutl(iplane) endif c c -- Below we see cases 5,6, and 7. c c *** M1 - ELSE *** else wfttil(iplane) = 0.0 wftutl(iplane) = (tfrdp(iplane) - tthawd(iplane)) * 1 pftutl(iplane) c c *** M1 - ENDIF *** endif c c -- The following is if the top frost has completely gone out. c c *** N1 - IF *** if (tthawd(iplane) .lt. tfrdp(iplane)) then if (tthawd(iplane) .lt. tilld(iplane)) then uftild(iplane) = tthawd(iplane) else uftild(iplane) = tilld(iplane) endif c c -- Otherwise, if no top-frost exists, we don't care about uftild... c c *** N1 - ELSE *** else tfrdp(iplane) = 0.0 tthawd(iplane) = 0.0 amtfrz(iplane) = 0.0 wfttil(iplane) = 0.0 wftutl(iplane) = 0.0 pfttil(iplane) = 0.0 pftutl(iplane) = 0.0 uftild(iplane) = tilld(iplane) c c *** N1 - ENDIF *** endif c c -- Need to keep track of freeze-thaw cycles... c c ----- Increase freeze/thaw cycle counter (Also in WMLT.) if ((olddep .lt. 0.0001).and.(tthawd(iplane).gt.0.0001)) 1 cycle(iplane) = cycle(iplane) + 1 c amtfrz(iplane) = wftill(iplane) + wfutil(iplane) 1 + wfttil(iplane) + wftutl(iplane) c return end