subroutine snowd(irtype,denh2o,iplane,driftf,driftg,snodep, 1 densgy,densgt,smelt,hour) c c +++ PURPOSE +++ c This subroutine handles snowdrift and snow depth calculations. c It was originally incorporated into WINTER, but was split out c by Charles R. Meyer 6/4/96 to reduce code complexity. c c +++ ARGUMENT DECLARATIONS +++ integer irtype,iplane,hour real denh2o,driftf,driftg,snodep,densgy,densgt,smelt c c +++ARGUMENT DEFINITIONS+++ c irtype - c iplane - c hour - c denh2o - c driftf - c driftg - c snodep - c densgy - c densgt - c smelt - c c +++ PARAMETERS +++ include 'pmxpln.inc' include 'pmxnsl.inc' include 'pmxhil.inc' c c +++ COMMON BLOCKS +++ include 'cwint.inc' include 'cclim.inc' c c -------------------------- c c -- SNOW DRIFT calculations are found in the sndrft.for file. c -- It's purpose is to calculate the amount of drifting and/or c -- scouring of new and old snow. c c -------------------------- c c -- We know must calculate the depth and density of the snowpack c -- layer given the weather conditions for the past hour. c snodep = snodpy(iplane) snodpt(iplane) = snodep densgy = densg(iplane) densgt = densgy c c -- Daily snow settling factor. c if (hour .eq. 1) wdayct(iplane) = wdayct(iplane) + 1 c if (hrsnow(hour) .gt. 0.0) wdayct(iplane) = 1 c c -- When there is no snow on the ground. c if (snodep .le. 0.0) then c c -- If not snowing, then nothing happens. c if (hrsnow(hour) .le. 0.0) then densgt = 0.0 densgy = 0.0 snodep = 0.0 snodpt(iplane) = 0.0 wmelt(iplane) = 0.0 c c -- Otherwise, it is snowing and possibly drifting. c else snodep = hrsnow(hour) + driftf snodpt(iplane) = snodep densgt = 100.0 densgy = densgt wmelt(iplane)=0.0 c -- Now we have to make sure that all of the snow didn't blow away. c if (snodpt(iplane) .le. 0.0) then snodpt(iplane) = 0.0 densgt = 0.0 densgy = densgt endif c endif c c -- Otherwise, we look at the case where there was existing snow c -- on the ground. c else c c -- If there is no melting occurring. c if (hrtemp .lt. -4.0) then c c -- If it is not snowing...but is possibly drifting. c if (hrsnow(hour) .le. 0.0) then snodep = snodpt(iplane) + driftg densgt = densgy c c -- If all snow on the ground has blown away. c if (snodep .le. 0.0) then snodep = 0.0 snodpt(iplane) = 0.0 densgt = 0.0 densgy = 0.0 endif c c -- Now we must account for the daily snow settling factor. C ----- Equation 3.7.1 CRM -- 10/27/95 setf = ((exp(-(float(wdayct(iplane)) * 2.0)))*0.0416667)+1.0 densgt = densgt * setf c c -- We set an upper cap for snow density. c if (densgt .gt. 522) densgt = 522 c C ----- Equation 3.7.2 CRM -- 10/27/95 snodep = snodep * (densgy / densgt) c c --------- It is snowing. Now we check if it is possibly drifting. else snodep = snodpt(iplane) + hrsnow(hour)+driftf+driftg C ----- Equation 3.7.3 CRM -- 10/27/95 C - XXX -- Shouldn't "snodpt" be replaced by the snow depth *yesterday* ? densgt = ((densgy * (snodpt(iplane) + driftg)) + 1 (100 * (hrsnow(hour) + driftf))) / snodep endif c c -- This ELSE is if temp > -4...ie. melt is occurring. c else snodep = snodpt(iplane) + hrsnow(hour) + driftf + driftg C ----- Equation 3.7.3 CRM -- 10/27/95 C - XXX -- Shouldn't "snodpt" be replaced by the snow depth *yesterday* ? densgt = ((densgy * (snodpt(iplane) + driftg)) + 1 (100 * (hrsnow(hour) + driftf))) / snodep c if (snodep .gt. 0.0) then c c -- Note that melt is calculated in terms of meters of water melted. c call melt(irtype,wrain,hour) c if (wmelt(iplane) .gt. 0.0) then densgy = densg(iplane) smelt = (wmelt(iplane) * denh2o) / densgy endif c c -- Nothing to melt. c else smelt = 0.0 wmelt(iplane) = 0.0 endif c c -- Now we must check the new depths. c snodep = snodpt(iplane) - smelt c c -- If all of the snow has been melted, we must recalculate wmelt. c if (snodep .le. 0.0) then snodep = 0.0 densgt = 0.0 c c -- In the below equation, we take the depth of snow that was left, c -- and convert it into an equivalent depth of melted water. This c -- is the maximum amount of melt that could have taken place for c -- the hour. The 1000 divisor is the density of water (Kg/m^3). c wmelt(iplane) = (snodpt(iplane) * densg(iplane)) * 0.001 c -- Otherwise, not all of the snowpack melted. It just became more c -- dense. c else densgt = densg(iplane) * (snodpt(iplane) / snodep) c c -- Up until the density of the snowpack reaches 350 Kg/m^3 we don't c -- allow any water-runoff to occur. The melting just accounts for c -- increasing the snowpack density. However, after 350 Kg/m^3 we c -- calculate an amount of water available for runoff/infiltration. c if (densgt .le. 350) then wmelt(iplane) = 0.0 c else wmelt(iplane) = ((densgt - 350) * snodep) * 0.001 densgt = 350 endif c endif c c -- Now we add in new snow after the melt and calculate our new c -- snow density. c if (hrsnow(hour) .gt. 0.0) then snodpt(iplane) = snodep densg(iplane) = densgt snodep = snodpt(iplane) + hrsnow(hour) + driftf + driftg C ----- Equation 3.7.5 CRM -- 10/27/95 C - XXX -- Note: This equation differs from the on in the User Doc. densgt = ((densg(iplane)) * (snodpt(iplane) + driftg) + 1 (100 * (hrsnow(hour) + driftf))) / snodep endif endif c c -- Set the global variables and move on to next hour/day. c if (densgt .gt. 522) densgt = 522 if (snodep .le. 0.0) densgt = 0.0 c endif snodpt(iplane) = snodep snodpy(iplane) = snodpt(iplane) densg(iplane) = densgt c return end