subroutine mltfdp(qdry,qout) c c +++PURPOSE+++ c This function calculates the change in the frost depth when c there is no heat flow into or out of the soil from above. c This accounts for melting that occurs from the bottom-up and c this is called whenever there is less heat passed down through c the layer from above than is coming up from the bottom or if c there is any thaw layer present in the frozen layer system. 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 qdry,qout c c +++ARGUMENT DEFINITIONS+++ c qdry - Heat flow from stable soil temperature to the bottom c of the frost layer (W/m^2). c qout - Heat flow across the four frozen layers (W/m^2). 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: frdp(mxplan),wftill,pfwutl,amtfrz, c thdp(mxplan) c include 'cstruc.inc' c read: iplane include 'cupdate.inc' c read: sdate c c +++LOCAL VARIABLES+++ save real tflux,lhfh2o,olddep,incr c c +++LOCAL DEFINITIONS+++ c tflux - Temporary flux value represents available energy (W/m^2). c lhfh2o - Latent heat of fusion of water (W/hr/m^3). c olddep - Original frost depth going into this routine (m). c incr - Calculated amount that the frost depth moves up (m). c c +++DATA INITIALIZATIONS+++ data lhfh2o/9.3027e04/ c c +++END SPECIFICATIONS+++ c c LIMIT TO PREVENT UNDERFLOWS ELSEWHERE - dcf 6/1/94 c c frdp(iplane) = amax1(frdp(iplane), 0.0) if(frdp(iplane) .lt. 0.0001) frdp(iplane) = 0.0 olddep = frdp(iplane) c *** L1 - IF *** if (frdp(iplane) .ge. 0.0001) then c -- Tflux is the net heat available to melt frost, ie. heat from c -- below - cold from above. if (thdp(iplane) .gt. 0.0) then tflux = qdry else tflux = qdry - qout endif c -- Now, if all the frost is in the till layer, then melting occurs c -- in the frozen tilled soil only (1 & 2). if (frdp(iplane) .le. tilld(iplane)) then incr = amax1(tflux / (lhfh2o * pftill(iplane)), 0.0) incr = min(incr,0.01) frdp(iplane) = frdp(iplane) - incr if(frdp(iplane).lt.0.0001) frdp(iplane)=0.0 c frdp(iplane) = amax1(frdp(iplane) - incr , 0.0) c -- Otherwise, the frost extends into the untilled layer. Then c -- melting occurs in the frozen untilled soil only. This is cases c -- 3,4,5,6 and 7. else c c XXX Note - temporary fix of BOMB here - for test data set #2 - c getting here and having a value of PFWUTL = 0.0, then c have a divide by zero. The reason why PFWUTL is zero c needs to be determined and fixed. dcf 2/25/94 c if(pfwutl(iplane).gt.0.0)then incr = amax1(tflux / (lhfh2o * pfwutl(iplane)) , 0.0) frdp(iplane) = amax1(frdp(iplane) - incr , 0.0) endif endif c -- Here, we don't allow melt to pass across till boundary. if ((olddep .gt. tilld(iplane)) .and. 1 (frdp(iplane) .lt. tilld(iplane))) 2 frdp(iplane) = tilld(iplane) c -- If frost is gone. c *** L2 - IF *** if (frdp(iplane) .le. thdp(iplane)) then c -- If the main layer all melts and a top-frost layer exists, we c -- swap the main-frost with the top-frost. c c CHANGE LIMIT TO PREVENT UNDERFLOW dcf 6/1/94 c if (tfrdp(iplane) .gt. 0.0) then if (tfrdp(iplane) .ge. 0.0001) then frdp(iplane) = tfrdp(iplane) thdp(iplane) = tthawd(iplane) tfrdp(iplane) = 0.0 tthawd(iplane) = 0.0 wftill(iplane) = wfttil(iplane) wfutil(iplane) = wftutl(iplane) pftill(iplane) = pfttil(iplane) pfwutl(iplane) = pftutl(iplane) wfttil(iplane) = 0.0 wftutl(iplane) = 0.0 pfttil(iplane) = 0.0 pftutl(iplane) = 0.0 c -- If main-frost is gone, and no top-frost, it all is gone. else frdp(iplane) = 0.0 thdp(iplane) = 0.0 tfrdp(iplane) = 0.0 tthawd(iplane) = 0.0 wftill(iplane) = 0.0 wfutil(iplane) = 0.0 pftill(iplane) = 0.0 pfwutl(iplane) = 0.0 wfttil(iplane) = 0.0 wftutl(iplane) = 0.0 pfttil(iplane) = 0.0 pftutl(iplane) = 0.0 uftild(iplane) = tilld(iplane) endif c *** L2 - ELSE *** else c -- Otherwise, a frost layer still exists. if (frdp(iplane) .le. tilld(iplane)) then c -- Cases 1 and 2. wftill(iplane) = (frdp(iplane) - thdp(iplane)) 1 * pftill(iplane) wfutil(iplane) = 0.0 else c -- Frost extends into the untilled layer. This follows cases 3,4 c -- 5,6, and 7. if (thdp(iplane) .ge. tilld(iplane)) then c -- Cases 5, 6, and 7. wftill(iplane) = 0.0 wfutil(iplane) = frdp(iplane) - thdp(iplane) 1 * pfwutl(iplane) else c -- The frost extends across the till boundary. Cases 3 and 4. wftill(iplane) = (tilld(iplane) - thdp(iplane)) 1 * pftill(iplane) wfutil(iplane) = (frdp(iplane) - tilld(iplane)) 1 * pfwutl(iplane) endif endif c *** L2 - ENDIF *** endif c -- All gone again. c *** L1 - ELSE *** else frdp(iplane) = 0.0 thdp(iplane) = 0.0 tfrdp(iplane) = 0.0 tthawd(iplane) = 0.0 wftill(iplane) = 0.0 wfutil(iplane) = 0.0 pftill(iplane) = 0.0 pfwutl(iplane) = 0.0 wfttil(iplane) = 0.0 wftutl(iplane) = 0.0 pfttil(iplane) = 0.0 pftutl(iplane) = 0.0 uftild(iplane) = tilld(iplane) c *** L1 - ENDIF *** endif amtfrz(iplane) = wftill(iplane) + wfutil(iplane) return end