subroutine frost(hour,smtill,smutil,bdtill,bdutil) c c +++PURPOSE+++ c This is the main driver for the frost subroutines. Based c on the surface, climate and soil conditions it decides c which subroutines are to be called. It considers energy c flow between the frozen layers, snow depth, melting or c freezing and also calls the infiltration capacity calculations. c c Authors(s): John Witte, UofMn WCES @ USDA-ARS-NCSRL c Date: 04/05/93 c c Verified and tested by Reza Savabi, USDA-ARS, NSERL 317-494-5051 c August 1994 c c Recoded by Charles R. Meyer Winter of '96 c Changes: (1) Eliminated calls to AVOID2, CAQWET, and CAQDRY. c Incorporated their code into FROST. c (2) Added greatly modified diagram from CWINT.INC below c and added many comments related to it. c (3) Added check to see if dummy parameters have changed c since last call to FROST. If not a lot of work can c be avoided. c Incorporated into WEPP code by Dennis Flanagan, 3/97 c c ********************************************************************* c * From CWINT.INC: * c ********************************************************************* c * c * c * The Frozen Layered System: WEPP c * -------------------------- Variable c * Name: c * c * ********************************* - Snow depth snodpt c * c * c * _________________________________ - Residue depth c * XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX - Residue c * ================================= <-- Surface c * ((((((((((((((((((((((((((((((((( c * )))) Soil thawed from above ))))) c * --------------------------------- <-- Top Thaw depth tthawd c * //// Soil frozen from above ///// c * --------------------------------- <-- Top Frost depth tfrdp c * )))) Soil thawed from above ))))) c * --------------------------------- <-- Thaw depth thdp c * \\\\\\\\\\\\ Frozen \\\\\\\\\\\\ c * \\\\\\\\\\\\ Soil \\\\\\\\\\\\\ c * --------------------------------- <-- Frost depth frdp <= 1.0m c * (((((((((( Stable Temp (((((((((( (thawed from below, c * ))))) Always Above Freezing ))))) or entire profile c * has thawed.) c * c * (Note that Top Thaw depth and Top Freeze depth can swap positions.) c * c * c * NOTES: c * ====== c * c * Tilled Layer - Depth from soil surface to the primary tillage c * ------------ depth for the season previous. (tilld = 0.2m) c * c * Untilled Layer - Depth from bottom of the tilled layer to the c * -------------- stable soil depth. c * c * Stable Soil Depth - Depth @ which soil temperature is stable... c * ----------------- this model assumes this depth is 1 meter below c * the lowest 0-degree isotherm depth. If no c * frost is present, this depth is left at 1 m. c * c * Simplifying Assumption c * (to avoid multiple frozen layers below the soil surface): c * If top frost layer gets more than 1/2 thawed, frost "disappears". c * If top thaw layer gets more than 1/2 frozen, thaw "disappears". c * c ********************************************************************* c c c +++ARGUMENT DECLARATIONS+++ integer hour real smtill,smutil,bdtill,bdutil c c +++ARGUMENT DEFINITIONS+++ c hour - The hour of the day that we are calculating. c smtill - Soil moisture in the tilled soil layer (%/100). c smutil - Soil moisture in the untilled soil layer (%/100). c bdtill - Bulk density of the tilled soil layer (Kg/HA). c bdutil - Bulk density of the untilled soil layer (Kg/HA). c c +++PARAMETERS+++ include 'pmxhil.inc' include 'pmxnsl.inc' include 'pmxpln.inc' c c +++COMMON BLOCKS+++ include 'cclim.inc' c read: tave,vwind,hradmj,hrtemp,rpoth,rad c include 'cwint.inc' c read: snodpt(iplane),tfrdp(mxplan),tthawd(mxplan),frdp(mxplan), c thdp(mxplan),densg(mxplan) c c include 'ccrpout.inc' c read: bd c include 'cwater.inc' c read: solthk(mxnsl,mxplan),fctill(mxplan),fcutil(mxplan) c include 'cstruc.inc' c read: iplane c c +++LOCAL VARIABLES+++ c c XXX Use of a global save does not follow WEPP coding conventions. c ONLY variables which need to be reused next time into FROST c should be saved. dcf 5/18/94 c save real qdry,qwet,qout,kutill,kuutil,ksnow,kftill, 1 kfutil,kres c ----- Saved Variables used to avoid unnecessary recalculation: integer sflag real smtilx,smutix,bdtilx,bdutix c c +++LOCAL DEFINITIONS+++ 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 qout - Heat flow across the four frozen layers (W/m^2). c kutill - Thermal conductivity of unfrozen tilled soil (W/hr C). c kuutil - Thermal conductivity of unfrozen untilled soil (W/hr C). c ksnow - Thermal conductivity of the snow pack (W/hr C). c kftill - Thermal conductivity of frozen tilled soil (W/hr C). c kfutil - Thermal conductivity of frozen untilled soil (W/hr C). c kres - Thermal conductivity of residue layer (W/hr C). c sflag - Set (sflag=1) if smtill,smutil,bdtill,bdutil are unchanged c since previous call. c c + + + DATA INITIALIZATIONS + + + data smtilx,smutix,bdtilx,bdutix /4*-1/ c c +++END SPECIFICATIONS+++ c c ----- See if input parameters have changed. If not, some expensive c calculations can be avoided. if(smtill.ne.smtilx) then sflag = 0 smtilx = smtill smutix = smutil bdtilx = bdtill bdutix = bdutil else if(smutil.ne.smutix) then sflag = 0 smutix = smutil bdtilx = bdtill bdutix = bdutil else if(bdtill.ne.bdtilx) then sflag = 0 bdtilx = bdtill bdutix = bdutil else if(bdutil.ne.bdutix) then sflag = 0 bdutix = bdutil else sflag = 1 endif c c ----- If some of the inputs have changed.... c ----- *** K1 IF *** if(sflag.eq.0) then c c ------ We can start by initializing some of these variables... kftill = 1.75 kfutil = 2.1 kres = 0.168 c c ------- Calculate Thermal conductivity of unfrozen tilled soil tmpvr1 = 0.5096 + 7.4493 * smtill - 8.7484 * smtill ** 2 kutill = tmpvr1 * (0.0014139*bdtill - 1.0588) c c c ------- Calculate Thermal conductivity of unfrozen untilled soil tmpvr1 = 0.5096 + 7.4493 * smutil - 8.7484 * smutil ** 2 kuutil = tmpvr1 * (0.0014139*bdutil - 1.0588) c c ----- *** K1 ENDIF *** endif c c ************************************************************* c *** BALANCE the HEAT FLOW from the unfrozen soil layers *** c *** to the frozen layers. Note that the temperature is *** c *** normally assumed constant (at 7 degrees Celsius) at *** c *** depths 1 meter below the freezing front. *** c ************************************************************* c c ****************** c *** DRY LAYERS *** c ****************** c ------ Calculate QDRY, heat flow from [the dry layers] beneath c the frost layer. c c ------- If there is tillage below the frost layer.... c ------- *** K2 IF *** if ((tilld(iplane) .ge. 0.0001) .and. 1 (tilld(iplane) .gt. frdp(iplane))) then c ------ thickness of the layer below the frost layer uftill = tilld(iplane) - frdp(iplane) c ------ Denominator for QDRY denom = ((1.0 - uftill) * kutill) + (uftill * kuutil) c ------ Denominator for QWET, the heat of fusion released by freezing c water at the freezing front. denom2 = ((1.0 - uftill) * usatkt(iplane)) + 1 (uftill * usatku(iplane)) c c Note: The calculations for DENOM & DENOM2 involve 2 soil layers -- c (a) unfrozen tilled layer (top); (b) unfrozen one-meter-thick c untilled temperature transition layer (bottom). The soil c above these 2 layers is frozen, and no longer of interest. c The soil below these is assumed to be at a constant 7 degrees c Celsius, and is also not of interest. c c ------- K uf * [T uf / Z uf] part of eqn 3.8.4 in WEPP manual. if (denom .gt. 0.0) then qdry = kutill * kuutil * 7.0 / denom else qdry = 0.0 endif c c ----- *** K2 ELSE *** else qdry = 7.0 * kuutil denom2 = usatkt(iplane) c ----- *** K2 ENDIF *** endif c c ****************** c *** WET LAYERS *** c ****************** c The EXPO portion of the equation ("ten" ?) here is a multiplier c 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 ----- Calculate QWET, heat of fusion released by freezing water that c is held in place at, or migrates to, the freezing front. c ----- L * K w [P / Z uf] part of eqn 3.8.4. c c ------- If the denominator of the heat flow equation in positive c ------- *** K3 IF *** if (denom2 .gt. 0.0) then c The following line of code was replaced because it only c calculates a PWATER value for the top soil layer (tilled) c and we should actually calculate an average PWATER value c for the tilled and untilled zones. dcf 3/14/97 c pwater = st(1,iplane) / dg(1,iplane) + thetdr(1,iplane) c temp = 0.0 pwatil = 0.0 do 100 i=1,2 temp = st(i,iplane) / dg(i,iplane) + thetdr(i,iplane) pwatil = pwatil + temp 100 continue pwatil = pwatil * 0.5 temp = 0.0 pwatut = 0.0 do 200 i = 3, nsl(iplane), 1 temp = st(i,iplane) / dg(i,iplane) + thetdr(i,iplane) pwatut = pwatut + temp 200 continue if(nsl(iplane)-2 .gt. 0)then pwatut = pwatut / float(nsl(iplane) - 2) else pwatut = pwatil endif c c The following groups of if/else logic were taken directly c from the information in subroutine CAQWET. Corrections c have been made to use the fcu***,por***,slo***,and ten*** c values for the respective tilled and untilled layers, which c was not correctly implemented in v95.7. dcf 3/12/97 c if(thdp(iplane).gt.0.0001) then if(tilld(iplane).gt.0.0001) then portil = 1.0 - bdtill*0.00037735849056 slotil = (fctill(iplane) - portil) * 65.853854492 c ten = 10.0**((100.0*(portil-pwater) / abs(slotil)) - 2.0) ten = 10.0**((100.0*(portil-pwatil) / abs(slotil)) - 2.0) else porutl = 1.0 - bdutil*0.00037735849056 sloutl = (fcutil(iplane) - porutl) * 65.853854492 c ten = 10.0**((100.0*(porutl-pwater) / abs(sloutl)) - 2.0) ten = 10.0**((100.0*(porutl-pwatut) / abs(sloutl)) - 2.0) endif qwet=9.3027e04*usatkt(iplane)*usatku(iplane)*ten / denom2 else if(frdp(iplane).lt.tilld(iplane)) then portil = 1.0 - bdtill*0.00037735849056 slotil = (fctill(iplane) - portil) * 65.853854492 c ten = 10.0**((100.0*(portil-pwater) / abs(slotil)) - 2.0) ten = 10.0**((100.0*(portil-pwatil) / abs(slotil)) - 2.0) else porutl = 1.0 - bdutil*0.00037735849056 sloutl = (fcutil(iplane) - porutl) * 65.853854492 c ten = 10.0**((100.0*(porutl-pwater) / abs(sloutl)) - 2.0) ten = 10.0**((100.0*(porutl-pwatut) / abs(sloutl)) - 2.0) endif endif qwet=9.3027e04*usatkt(iplane)*usatku(iplane)*ten / denom2 c ------- *** K3 ELSE *** else qwet = 0.0 c ------- *** K3 ENDIF *** endif c c c ----- *** L1 IF *** c ----- If the surface temperature is below freezing.... if (surtmp(hour) .le. 0.0) then c ------- Calculate heat flow through the top layers to extend the top c frost depth. c c ------- Calculate Thermal conductivity of the snow pack if (snodpt(iplane) .gt. 0.001) then ksnow = 0.0000029694 * densg(iplane) ** 2 else ksnow = 1.0 endif call caqout(ksnow,kres,kftill,kfutil,resdep,qout,hour) c c ------- *** L2 IF *** c ------- If a top frost layer exists and is frozen to the surface... c (This is cases 2a, 4a, 5a, 5c, and 7a.) if (thdp(iplane) .ge. 0.0001) then c --------- Determine whether to ignore the top frost layer or assume c it if frozen to the surface. c c --------- Prevent the occurence of multiple top frost layers in the soil. c If more than half of the top frost is melted, the top frost is c eliminated, otherwise it is assumed that there is no top thaw. c c --------- If there is no frozen sandwich.... if (tfrdp(iplane) .ge. 0.0001) then if(tthawd(iplane)/tfrdp(iplane) .gt. 0.5)tfrdp(iplane) = 0.0 tthawd(iplane) = 0.0 endif c c ----------- Cold freeze extends top frost depth. call cldfrz(qout,qdry,qwet,smtill,smutil) c c ----------- Melt frost from the bottom up from the stable soil heat. call mltfdp(qdry,qout) c c ------- The soil is frozen to the surface and no top frost layer exists. c ------- *** L2 ELSE *** else thdp(iplane) = 0.0 c --------- If heat from the bottom is greater than heat from above. if (qdry .gt. qout) then c ----------- Melt the frost from the bottom. call mltfdp(qdry,qout) else c ----------- Cold Freeze extends frost depth. call cldfrz(qout,qdry,qwet,smtill,smutil) endif c c ------- *** L2 ENDIF *** endif c c ----- *** L1 ELSE *** else c c ------- A frost layer does exist. c ------- *** M2 IF *** if (frdp(iplane) .ge. 0.0001) then c c --------- No snow if (snodpt(iplane) .lt. 0.00001) then c c ----------- The temperature is above zero and there is no top frost c present in the soil. So, when melt comes from the top c down, it will melt the main frost layer. c c ----------- If there's top-frost, melt it from the top-down... if (tfrdp(iplane) .ge. 0.0001) then call tfwmlt(kutill,kuutil,resdep,hour) c c ----------- No top frost... else tfrdp(iplane) = 0.0 call wmlt(kres,kutill,kuutil,resdep,hour) endif endif c c ------- This is the "no frost" or "do nothing" part... c ------- *** M2 ELSE *** else frdp(iplane) = 0.0 c c ------- *** M2 ENDIF *** endif c c ----- *** L1 ENDIF *** endif c c Check to make sure that computed frost/thaw depths do not c exceed soil thickness. dcf 3/17/97 c if(frdp(iplane).gt.solthk(nsl(iplane),iplane)) 1 frdp(iplane)=solthk(nsl(iplane),iplane) if(thdp(iplane).gt.solthk(nsl(iplane),iplane))then thdp(iplane) = 0.0 frdp(iplane) = 0.0 endif return end