subroutine melt(irtype,wrain,hour) c c c +++PURPOSE+++ c This subroutine calculates the amount of snow melt that c occurs on an hourly basis given the hourly weather and c surface conditions. The amount of melted water is calculated c in units of meters. c c Author(s): Cully Hession and Bruce Lucord, USDA-ARS-NCSRL c Revised by John Witte, UofMn WCES @ USDA-ARS-NCSRL c Date: 03/23/93 cc Verified and tested by Reza Savabi, USDA-ARS, NSERL 317-494-5051 c Modified 11/25/96 by Dennis Flanagan c c +++ARGUMENT DECLARATIONS+++ real wrain integer hour,irtype c c +++ARGUMENT DEFINITIONS+++ c wrain - Amount of daily measured rainfall (m/day). c hour - Hour of the day (hr). c irtype - crop residue type index c c +++PARAMETERS+++ include 'pmxpln.inc' include 'pmxnsl.inc' include 'pmxhil.inc' include 'pmxtil.inc' include 'pmxtls.inc' include 'pntype.inc' include 'pmxres.inc' c c +++COMMON BLOCKS+++ include 'ccrpout.inc' c read: ccr include 'cclim.inc' c read: hradmj,tmax,vwind,hrdewp,hrtemp,rpoth include 'cdecvar.inc' c read: cuthgt(ntype) include 'cstruc.inc' c read: iplane include 'ccover.inc' c read: canhgt,cancov,lanuse include 'cupdate.inc' c read: sdate include 'cwint.inc' c read: snodpy,densg include 'crinpt6.inc' c read: rrough(mxplane) c c +++LOCAL VARIABLES+++ save real amelt1,bmelt1,cmelt1,dmelt1,m1 real height,rough,disp,cldpct,amelt,bmelt,cmelt,windh, 1 dmelt,x,tmptrm c c +++LOCAL DEFINITIONS+++ c height - Height of the canopy which is above the snow layer (m). c rough - Surface roughness of the top layer of the system (m). c disp - Zero-plane displacement of the system's surface (m). c cldpct - Percent cloud cover for the given hour (decimal %). c amelt - First part of WEPP eqn 3.2.1 - represents radiation term. c bmelt - Second part of WEPP 3.2.1 - represents cloud cover term. c tmptrm - Temporary term used in equations. c cmelt - Third partof WEPP 3.2.1 - represents wind term. c windh - Height at which the wind speed is measured, assume 2.0 m. c dmelt - Fourth part of WEPP 3.2.1 - represents temperature term. c x - Temporary variable. c adj - adjustment for wind height and roughness calculations c +++DATA INITIALIZATIONS+++ data windh/2.0/ c c +++END SPECIFICATIONS+++ c -- Calculate the surface conditions -- if(lanuse(iplane).eq.1)then height = cuthgt(irtype) if(height.lt.canhgt(iplane))height = canhgt(iplane) else height = canhgt(iplane) endif c rough = 0.0 if(height .lt. 0.001) rough = 0.0002 if(snodpy(iplane) .lt. 0.01) rough = rrc(iplane) if(height .lt. 0.001) goto 10 height = height - snodpy(iplane) if(height .lt. 0.001) then height = 0.0 rough = 0.0002 else rough = 0.13 * height endif 10 disp = 0.77 * height c c Savabi correction - need to convert delta temperature c from degrees C to degrees F - dcf 11/25/96 c hrtef=hrtemp*(9./5.) c -- Calculating decimal percent cloud cover -- if (rpoth .lt. 0.0001) then cldpct = 1.0 else cldpct=(1.0 - radmj / rpoth) / 0.7 endif if (cldpct .gt. 1.0) cldpct = 1.0 c -- Aug 1989 WEPP manual. First term of 3.2.1 for melt. -- c -- NOTE that 0.0607 is a conversion from MJ to Ly units. -- amelt = 0.0607 * hradmj * (1.0 - cancov(iplane)) c -- NOTE: Canopy may not be the correct variable to use in the c -- WEPP version of winter. -JW c -- Check with Bob to see if cutght(iplane) should be used... c -- Allow some melt when max temp > 0 but mean temperature < 0 -- if(hrtemp .gt. -4.0 .and. hrtemp .le. 0.0) then tmptrm = amelt amelt = tmptrm * (0.36 * hrtemp + 1.0) endif c -- Aug 1989 WEPP manual. Second term of 3.2.1 for melt. -- bmelt =( 0.84 * (1.0 - cldpct))/24. etm(iplane) = etm(iplane) + (.0254 * bmelt) c -- Aug 1989 WEPP manual. Third term of 3.2.1 for melt. c -- This safety check was placed in here to account for case c -- of rough being unset and defaulted to 0.0. if (rough .lt. 0.0001) then cmelt = 0.0 else xtem = windh - disp + rough if(xtem .gt. 0.0)then adj = 1.0 - (1.0 / (log(xtem) / rough)) else adj = 1.0 endif cmelt1 = 0.0188 * vwind * (1.0 - 0.8 * cancov(iplane)) * 1 (0.396 * hrtemp + 1.404 * hrdewp) * adj c c Savabi correction to convert wind velocity from m/sec c to miles/hour. Assume 0.1 tx=0.22 hrtef - dcf 11/25/96 vwmph=(vwind*3600)/1609. hrdtf=hrdewp*(9./5.) cmelt2 = (0.0084/24.) * vwmph*(1.0-0.8*cancov(iplane))* 1 ((0.22*hrtef)+(0.78*hrdtf))*adj endif c -- Aug 1989 WEPP manual. Fourth term of 3.2.1 for melt. c Equation commented out. New equation for dmelt1 from c Savabi is given below. - dcf 11/25/96 c dmelt1 = hrtemp * (0.0382 + 0.014 * cldpct) + hrtemp * c 1 0.000496 * hrrain(hour) rainin=hrrain(hour)*39.37 dmelt2=(hrtef)*(0.025/24.+(0.007/24.*rainin)) c -- Put terms into the following eqn calculates water melt (m). -- c -- NOTE that 0.0254 is a conversion factor from inches to meters. -- wmelt(iplane) = 0.0254 * (amelt - bmelt + cmelt2 + dmelt2) amelt1 = hradmj*23.9 * (1.0 - cancov(iplane)) if(hrtemp .gt. -4.0 .and. hrtemp .le. 0.0) then tmptrm = amelt1 amelt1 = tmptrm * (1-0.7*cldpct)*(0.0025*hradmj*23.9)* 1 (2*0.22*(hrtemp*9./5.+32.)+1) endif bmelt1 = 0.84*(1-cldpct)*(1.0 - cancov(iplane)) if (rough .lt. 0.0001) then cmelt1 = 0.0 else cmelt1 = 0.0084 * vwind*(1608./3600.) * 1 (1.0 - 0.8 * cancov(iplane)) * 1 (0.22*(hrtemp*9./5.+32.) + 0.78 * (hrdewp*9./5.+32.)) endif c c New equation for dmelt1 from Savabi - dcf 11/25/96 dmelt1 = (hrtemp*9./5.+32.)*(0.025+0.007*hrrain(hour)/0.025) m1 = 0.0254*(amelt1 - bmelt1 + cmelt1 + dmelt1) c -- This check converts the water depth to a depth of equal density -- c -- to that of the snowpack. -- c -- The "1000" is the density of water (kg/m^3). c Savabi new line to prevent negative melt values - dcf 11/25/96 if(wmelt(iplane).le.0.0)wmelt(iplane)=0.00 x = wmelt(iplane) * 1000.0 / densg(iplane) if(x .ge. snodpy(iplane)) 1 wmelt(iplane) = snodpy(iplane) * (densg(iplane) / 1000.0) if(wmelt(iplane).lt.0.0) wmelt(iplane) = 0.0 return end