subroutine cldfrz(qout,qdry,qwet,smtill,smutil) c c +++PURPOSE+++ c This function is responsible for extending the frost depth c if qwet > qout - qdry (or) b > Qf - a. c c The purpose of this program is to extend the frost depth with c excess heat flow when the heat flow through the frozen layer c system is greater than heat flow from the thermal conductivity c of the unfrozen soil. The excess heat flow is balanced by the c heat of fusion released by freezing water. 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 Reworked by Charles R. Meyer, with the McCabe Tool, March 1996 c c Inserted into WEPP code 2/3/97, dcf c c c c OVERVIEW OF MODULE LOGIC STRUCTURE: c c c | c | c (thaw depth = 0) | (thaw depth > 0) c --------------------------------- c | [L1] | c | | c | | c (till depth > | (till depth < (till depth > | (till depth < c frost depth) | frost depth) frost depth) | frost depth) c ----------------- ----------------- c | [L2] | | [M2] | c | | | | c | | | | c ----------------- ----------------- c | | c | | c --------------------------------- c | c | c | c c c c +++ARGUMENT DECLARATIONS+++ real qout,qdry,qwet,smtill,smutil c c +++ARGUMENT DEFINITIONS+++ c qout - Heat flow across the four frozen layers (W/m^2). c qdry - Heat flow from stable soil temperature to the bottom c of the frost layer (W/m^2). c qwet - Heat flow required to freeze H2O in the soil (W/m^2). c smtill - Soil moisture in the tilled layer (Dec %). c smutil - Soil moisture in the untilled layer (Dec %). 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: tfrdp(mxplan),amtfrz,wftill,pfwtil,wfutil,pfwut c include 'cstruc.inc' c read: iplane include 'cupdate.inc' c read: sdate 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 htreq,slmst,ceh2o,lhfh2o,olddep,otdep,incr c c +++LOCAL DEFINITIONS+++ c htreq - Heat (energy) left over, required to freeze soil (W/m^2). c slmst - Soil moisture in the soil layer (Dec %). c ceh2o - Coefficient of expansion for water (unitless). c lhfh2o - Latent heat of fusion of water (W/hr/m^3). c olddep - Frost depth @ beginning of the hour (m). c otdep - Top frost depth @ beginning of the hour (m). c incr - Depth that the frost will be extended on that hour (m). c c +++DATA INITIALIZATIONS+++ data ceh2o/1.1/ data lhfh2o/9.3027e04/ c +++END SPECIFICATIONS+++ c c LIMITS TO PREVENT UNDERFLOW ERRORS ELSEWHERE - dcf 6/1/94 c if(frdp(iplane) .lt. 0.0001) frdp(iplane) = 0.0 if(tfrdp(iplane) .lt. 0.0001) tfrdp(iplane) = 0.0 c c -- As usual, we being by doing some initialization... slmst = 0.0 olddep = frdp(iplane) otdep = tfrdp(iplane) if (tfrdp(iplane) .lt. thdp(iplane)) qdry = 0.0 if (frdp(iplane) .lt. tilld(iplane)) then slmst = smtill else slmst = smutil endif c -- Soil is frozen to the surface. c *** L1 - IF *** if (thdp(iplane) .lt. 0.0001) then tmpvr1 = qout - qdry if (qwet .gt. tmpvr1) then htreq = tmpvr1 else htreq = qwet endif c *** L2 - IF *** c -- Note: We have already ensured that FRDP>=0, so just need FRDP<=TILLD c Checking this change - new code commented out - caused differences c dcf 3/6/97 c if(tilld(iplane) .gt. frdp(iplane)) then if ((tilld(iplane) .gt. 0.0 ) .and. 1 (frdp(iplane) .le. tilld(iplane))) then c -- Frost extension in the tilled soil layer. c -- NOTE: In the following calculations, we can use either smtill or c -- slmst. Currently, we are using smtill for readability. The c -- same is true for the untilled portion below. if (slmst .ge. 0.0001) then incr = (ceh2o * htreq) / (lhfh2o * smtill) incr = min(incr,0.01) frdp(iplane) = frdp(iplane) + incr endif c c LIMIT TO PREVENT UNDERFLOW dcf 6/1/94 c if(frdp(iplane) .lt. 0.0001)frdp(iplane) = 0.0 c -- The following check is put here to see if the frost crossed the c -- till boundary. We allow the frost to go up to, but not cross c -- the till layer in a single hour. if ((olddep .lt. tilld(iplane)) .and. 1 (frdp(iplane) .gt. tilld(iplane)))then frdp(iplane) = tilld(iplane) incr = tilld(iplane) - olddep htreq = (incr * lhfh2o * slmst) / ceh2o endif wftill(iplane) = wftill(iplane) + htreq / lhfh2o wfutil(iplane) = 0.0 if(frdp(iplane) .gt. 0.0) then pftill(iplane) = wftill(iplane) / frdp(iplane) else pftill(iplane)=0.0 endif pfwutl(iplane) = 0.0 c -- Otherwise, there may or may not be a tilled layer or c -- a frost layer but frost extends into untilled soil. c *** L2 - ELSE *** else incr = (ceh2o * htreq) / (lhfh2o * smutil) frdp(iplane) = frdp(iplane) + incr c c LIMIT TO PREVENT UNDERFLOW dcf 6/1/94 if(frdp(iplane) .lt. 0.0001)frdp(iplane) = 0.0 c wfutil(iplane) = wfutil(iplane) + htreq / lhfh2o c c XXX This equation may be incorrect if FRDP is less than TILLD c or if FRDP is set back to zero. dcf 6/1/94 c pfwutl(iplane) = wfutil(iplane) / 1 (frdp(iplane) - tilld(iplane)) c *** L2 - ENDIF *** endif amtfrz(iplane) = wftill(iplane) + wfutil(iplane) c -- Otherwise, there is a thawed layer. c *** L1 - ELSE *** else if (tilld(iplane) .gt. 0.0) then slmst = smtill else slmst = smutil endif htreq = qout if (slmst .gt. 0.0) then incr = (ceh2o * htreq) / (lhfh2o * slmst) tfrdp(iplane) = tfrdp(iplane) + incr endif c *** M2 - IF *** c -- All top frost is in the tilled layer and extends into c -- the tilled soil. if(tilld(iplane) .gt. tfrdp(iplane)) then c c LIMIT TO PREVENT UNDERFLOW dcf 6/1/94 if(tfrdp(iplane) .lt. 0.0001)tfrdp(iplane) = 0.0 c -- The following check put here to see if the top_frost crossed the c -- tillage boundary. We allow the frost to go up to, but not cross c -- the till layer in a single hour. if ((otdep .lt. tilld(iplane)) .and. 1 (tfrdp(iplane) .gt. tilld(iplane))) then tfrdp(iplane) = tilld(iplane) incr = tilld(iplane) - otdep htreq = incr * lhfh2o * smtill / ceh2o endif c -- This check is put in place to make sure that the frost layers c -- do not overlap. It is done in the same way as for the check c -- for the till layer. if ((otdep .lt. thdp(iplane)) .and. 1 (tfrdp(iplane) .gt. thdp(iplane))) then tfrdp(iplane) = thdp(iplane) incr = thdp(iplane) - otdep htreq = incr * lhfh2o * smtill / ceh2o endif wfttil(iplane) = wfttil(iplane) + htreq / lhfh2o wftutl(iplane) = 0.0 c c XXX Added if/else to prevent divide by zero in the c case that the TFRDP value is so small that it has c been reset to zero. dcf 6/1/94 c if(tfrdp(iplane) .gt. 0.0)then pfttil(iplane) = wfttil(iplane) / tfrdp(iplane) else pfttil(iplane) = 0.0 endif pftutl(iplane) = 0.0 c -- Otherwise, the frost extends into the untilled soil. c *** M2 - ELSE *** else c LIMIT TO PREVENT UNDERFLOW dcf 6/1/94 if(tfrdp(iplane) .lt. 0.0001)tfrdp(iplane) = 0.0 c -- This check is put in place to make sure that the frost c -- layers do not overlap. It is done in the same way as for c -- the check for the till layer. if ((otdep .lt. thdp(iplane)) .and. 1 (tfrdp(iplane) .gt. thdp(iplane))) then tfrdp(iplane) = thdp(iplane) incr = thdp(iplane) - otdep htreq = incr * lhfh2o * smtill / ceh2o endif wftutl(iplane) = wftutl(iplane) + htreq / lhfh2o c c Added following if/else to prevent new divide by zero c dcf 2/3/97 c if(tfrdp(iplane).gt.tilld(iplane))then pftutl(iplane) = wftutl(iplane) / (tfrdp(iplane) 1 - tilld(iplane)) else pftutl(iplane) = 0.0 endif c *** M2 - ENDIF *** endif amtfrz(iplane) = wftill(iplane) + wfutil(iplane) + 1 wfttil(iplane) + wftutl(iplane) c *** L1 - ENDIF *** endif c -- One last case to check is if the 2 frost layers meet... if (tfrdp(iplane) .ge. thdp(iplane)) then tfrdp(iplane) = 0.0 tthawd(iplane) = 0.0 thdp(iplane) = 0.0 wfttil(iplane) = 0.0 wftutl(iplane) = 0.0 pfttil(iplane) = 0.0 pftutl(iplane) = 0.0 endif return end