subroutine winter(wrain,snoflg) c c +++PURPOSE+++ c This subroutine calls the snowmelt, snowdrift, and frost c models from wepp, it will perform the waterbalance until c the frost is out and computes a new depth and density of c of snow on an hourly basis. The equations found in this c program are by Bob Young and George Benoit, etal. c c Authors(s): Cully Hession and Bruce Lucord, USDA-ARS-NCSRL c Revised by John Witte, UofMn WCES @ USDA-ARS-NCSRL c Date: 04/01/93 c c Verified and tested by Reza Savabi, USDA-ARS, NSERL 317-494-5051 c August 1994 c c Modified by Dennis Flanagan, 11/26/96 with Savabi changes. c Modified by Dennis Flanagan, 1/31/97 with Meyer changes. c c Changes: c 1) Made major changes to streamline logic -- primarily that c used to determine WCASE -- also "*1000's" in write statements c -- additionally, various block-if's converted to simple if's, c and divides converted to equivalent multiplys. As a result c Winter went from gobbling 53.6% of CPU time, to 36.7. c 2) Added embedded calling map. c c c Recoded from V-95.1 by Charles R. Meyer c 5/2x/96 - 6/04/96 c c Changes: c 3) The cyclomatic complexity v(g) of Winter was 81. Split out c W_CASE as a function and SNOWD, HR_TMP, STMTIM, RES_DP c as separate subroutines c to give v(g) of 33 (Winter), 11, 17, 12, & 11 respectively. c 4) Eliminated "goto 111" structure. c 5) Removed subroutine W_CASE, as it did not affect any other c part of the code. dcf 2/6/97 c c c WINTER c | c ----------------------------------------------------------- c | | | | | | | c SET_HC SUNMAP The Killer Loop SNOWD HR_TMP STMTIM RES_DP c | (do 500) (add) (add) (add) (add) c CALCHC | c | c ---------------------------------- c | | | | c RADCUR TMPADJ FROST MELT c | | c HRTMP | c | c -------------------------------------------------------------- c | | | | | | | | c CAQDRY CAQWET WMLT TFWMLT CAQOUT AVOID2 CLDFRZ MLTFDP c (gone) (gone) (gone) c c c +++ARGUMENT DECLARATIONS+++ real wrain integer snoflg c c +++ARGUMENT DEFINITIONS+++ c wrain - Daily rainfall amount (m). c snoflg - Flag used for winter output. c c +++PARAMETERS+++ include 'pmxelm.inc' include 'pmxcrp.inc' include 'pmxhil.inc' include 'pmxnsl.inc' include 'pmxpln.inc' include 'pmxpnd.inc' include 'pmxslp.inc' include 'pmxtil.inc' include 'pmxtls.inc' include 'pntype.inc' include 'pmxres.inc' c c +++COMMON BLOCKS+++ include 'cangie.inc' c read: radlat include 'cclim.inc' c read: tave,tmnavg,tmxavg include 'ccrpvr5.inc' c read: diam include 'ccons.inc' c read: bdcons include 'ccrpprm.inc' c read: iresd include 'ccover.inc' c read: lanuse include 'ccrpout.inc' c read: bd include 'chydrol.inc' c read: rain(mxplan) include 'cdiss11.inc' c read: dur include 'cstruc.inc' c read: iplane include 'cupdate.inc' c read: day,year,sdate include 'cwint.inc' c read: wmelt(mxplan),drift,snodpt(mxplan),azm,deglat include 'cparame.inc' c read: sm(mxplan) include 'cslope2.inc' c read: avgslp(mxplan) include 'cwater.inc' c read: solthk(mxplan,mxplan) include 'ccrpvr1.inc' c read: rmogt include 'ctillge.inc' c read: tildep(mxplan,mxplan) include 'cstmflg.inc' c read: jyear c c +++LOCAL VARIABLES+++ integer itflag,isflag,temphr integer hour,numres,snohrs,rnhrs,i,j real tmpvr1,tmpvr2,tmpvr3,tmpvr4,tmpvr5 real tmpvr6,tmpvr7,tmpvr8,tmpvr9 real densgy,denh2o,snodep,densgt,halfdy,langmj, 1 durtn,k(4),randn,wbdtil(mxplan),bdutil(mxplan), 1 smtill(mxplan),smutil(mxplan) real dgutil,porutl(mxplan),portil(mxplan),estrad real driftf,driftg,fdrft,gdrft,smelt,temp24,totmel c save wbdtil,bdutil,porutl save temp24 c c +++LOCAL DEFINITIONS+++ c densgy - Temporary variable representing snow density (Kg/m^3). c denh2o - Density of water (Kg/m^3). c snodep - Temporary variable representing snow depth (m). c snodep - Temporary variable representing snow depth after c adjustments have been made (m). c densgt - Temporty variable representing snow density after c adjustments have been made (Kg/m^3). c halfdy - Length of day from sunrise until noon (hrs). c hour - Hour counter...used in the looping structure (hr). c langmj - Conversion factor used to convert from Langly's to MJ. c temp - Temporary variable. c numres - Number of residue layers. c snohrs - Number of hours which there is snowfall throughout day. c rnhrs - Number of hours which there is rainfall throughout day. cc schrs - Number of hours which there is snow cover during the day. c durtn - Local variable keeping track of storm duration (h). c k - Variables used in random number generator. c randn - Random number used to calculate time of storm. c i - Seed used for random number generator. c hourx - Previous value of hour c smtilx - Previous value of smtill(iplane) c smutix - Previous value of smutil(iplane) c wbdtix - Previous value of wbdtil(iplane) c bdutix - Previous value of bdutil(iplane) c c + + + DATA INITIALIZATIONS + + + C data hourx,smtilx,smutix,wbdtix,bdutix /5 * -1/ c c +++END SPECIFICATIONS+++ c c -- First, we do a little initializing... c denh2o = 1000.0 langmj = 0.04184 halfdy = 0.0 snow = 0.0 snohrs = 0 rnhrs = 0 durtn = 0.0 i = sdate rans=0.0 snodpt(iplane)=snodpy(iplane) do 5 j=1, 24 hrmlt(j,iplane) = 0.0 5 continue c c added following line because variable "temp24" was undefined c dcf 8/16/94 c if(wntflg(iplane) .eq. 0)temp24 = tmin c c cycle(iplane) = 1 c c -- We borrowed CLIGEN'S random number gen. for storm time... k(1)=i k(2)=i k(3)=i k(4)=i k(4)=3*k(4)+k(2) k(3)=3*k(3)+k(1) k(2)=3*k(2) k(1)=3*k(1) i=k(1)/1000 k(1)=k(1)-i*1000 k(2)=k(2)+i i=k(2)/100 k(2)=k(2)-100*i k(3)=k(3)+i i=k(3)/1000 k(3)=k(3)-i*1000 k(4)=k(4)+i i=k(4)/100 k(4)=k(4)-100*i randn=((((k(1)*.001+k(2))*.01+k(3))*.001+k(4))*.01) * 24 wnttim = aint(randn) if (wnttim .lt. 1.0) then wnttim = 1.0 endif c smelt = 0.0 c if (frdp(iplane) .lt. 0.0001) then tfrdp(iplane) = 0.0 tthawd(iplane) = 0.0 frdp(iplane) = 0.0 thdp(iplane) = 0.0 pftill(iplane) = 0.0 pfwutl(iplane) = 0.0 endif c c -- We set the bulk density and soil moisture of the tilled soil c -- equal to that of the first layer (secondary tilled) and second c -- layer (primary tilled). c c -- As for the untilled soil, we compute the average bulk density c -- and soil moisture of the untilled soil from the remaining till c -- layers. c c The following logic is now executed every day in the winter c routines so that water balance is tracked properly. Otherwise, c a dry soil will never freeze - even if a lot of water added to c it as precipitation or snow melt. c tmpvr1 = dg(1,iplane) + dg(2,iplane) wbdtil(iplane) = ((bd(1,iplane) * dg(1,iplane)) + 1 (bd(2,iplane) * dg(2,iplane))) / tmpvr1 smtill(iplane) = (soilw(1,iplane) + soilw(2,iplane)) / tmpvr1 portil(iplane) = (2650.0 - wbdtil(iplane))/2650. bdutil(iplane) = 0.0 smutil(iplane) = 0.0 c c -- NOTE: We ran into problems here whenever we were dealing with c -- a soil profile containing less than 3 layers. So, what we did c -- is simulate an untilled layer using values from the tilled layer. c if (nsl(iplane) .lt. 3) then bdutil(iplane) = wbdtil(iplane) smutil(iplane) = smtill(iplane) porutl(iplane) = portil(iplane) else do 10 i=3, nsl(iplane) bdutil(iplane) = bdutil(iplane) 1 + bdcons(i,iplane) * dg(i,iplane) smutil(iplane) = smutil(iplane) + soilw(i,iplane) 10 continue tilld(iplane) = solthk(2,iplane) c smutil & bdutil are divided by dg(..) for the layer c below the till layer. dgutil = solthk(nsl(iplane),iplane) - solthk(2,iplane) bdutil(iplane)=bdutil(iplane)/dgutil smutil(iplane)=smutil(iplane)/dgutil c c Compute porosity of the untilled layer from the bulk density c value. dcf 11/26/96 c porutl(iplane) = (2650. - bdutil(iplane))/2650. if (frdp(iplane) .lt. tilld(iplane)) then uftild(iplane) = tilld(iplane) - frdp(iplane) else uftild(iplane) = 0.0 endif smtill(iplane)=smtill(iplane)/portil(iplane) smutil(iplane)=smutil(iplane)/porutl(iplane) wftill(iplane) =(tilld(iplane)-uftild(iplane))*smtill(iplane) if(wftill(iplane).gt.0.0)then pftill(iplane) = wftill(iplane) / 1 (tilld(iplane) - uftild(iplane)) else pftill(iplane)=0.0 endif if (frdp(iplane) .gt. tilld(iplane)) then wfutil(iplane) = (frdp(iplane) - tilld(iplane)) 1 * smutil(iplane) pfwutl(iplane) = wfutil(iplane) / 1 (frdp(iplane) - tilld(iplane)) else wfutil(iplane) = 0.0 pfwutl(iplane) = 0.0 endif endif tilld(iplane) = solthk(2,iplane) if (frdp(iplane) .lt. tilld(iplane)) then uftild(iplane) = tilld(iplane) - frdp(iplane) else uftild(iplane) = 0.0 endif c if (pftill(iplane) .le. 0.0) pftill(iplane) = smtill(iplane) if (pfwutl(iplane) .le. 0.0) pfwutl(iplane) = smutil(iplane) c c -- Now we must calculate the total residue depth (m). c call res_dp(lanuse(iplane),numres,iplane,resdep(iplane)) c c c -- The following converts rad to MJoules from Langley's c This will cause problem for ET calculation during c winter, because rad is converted to mJ in evap subroutine c c Note; radmj= radiation in mj/day, radly is ly/day Reza 3/9/94 temp = radly * langmj radmj = temp c c -- We now begin our 24 hour loop for each day... c hour = 1 totmel=0.0 etm(iplane) = 0.0 c c -- Calculate Hydraulic Conductivity values under given conditions. c call set_hc(smtill(iplane),smutil(iplane),hour,wbdtil(iplane), 1 bdutil(iplane)) isflag=0 itflag=0 if(rain(iplane).le.0.0) then if(snodpy(iplane).gt.0.0) then isflag=1 else if((tmax-tmin).le.1) itflag=1 endif endif c c -- Sunmap calculates the radiation on a sloping surface. call sunmap(sdate,radly,estrad,halfdy,dsunmp) hour=0 c c ------------------------------------------------------------ c The 500 hourly loop follows 500 continue hour=hour+1 c c -- Now we call the RADCUR routine to calc hourly radiation. c call radcur(sdate,hour,radlat,dsunmp,cratio) c c -- Now we know the ratio of daily radiation to hourly radiation. c reza, in the next line hradmj is in MJ c call hr_tmp(itflag,jyear,sdate,hour,temp24,radmj,halfdy, 1 estrad,cratio,rpoth,tmin,tmax,hradmj,hrtemp) c call stmtim(rain(iplane),stmdur,hour,wnttim,rans, 1 hrrain(hour),hrsnow(hour),hrtemp, 2 snodpt(iplane),rnhrs,snohrs) c snow = snow + hrsnow(hour) c -------------------------------------------------------------- c SNOW DRIFTING ROUTINES ARE CURRENTLY NOT ACTIVE - dcf 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 Set residue type index to be used to determine height of c residue for drifting computations. c c irtype = iresd(1,iplane) c reza we are not currently simulating snow drift 8/1/94 c reza if((snodpt(iplane).gt.0.0) .or. (hrsnow(hour).gt.0.0)) then c reza call sndrft(irtype,hour,driftf,driftg) c reza, test for water balance cal. during winter c reza fdrft = driftf c reza gdrft = driftg c reza schrs = schrs + 1 c reza else c reza fdrft = 0.0 c reza gdrft = 0.0 c reza endif c driftf = 0.0 driftg = 0.0 fdrft = 0.0 gdrft = 0.0 c -------------------------------------------------------------- c c -- Next, we move on to the FROST subroutine. This routine is c -- called for every hour that the winter simulation is run. c -- The main frost driver is located in the Sfrost.for file. c call frost(hour,smtill(iplane),smutil(iplane), 1 wbdtil(iplane),bdutil(iplane)) c c -- CHECK... Still have to move infiltration capacity calculations c -- out of frost to this point. c c -------------------------- c call snowd(iresd(1,iplane),denh2o,iplane,driftf,driftg, 1 snodep,densgy,densgt,smelt,hour) c c -- Set up the hourly arrays for thaw and melt. These will be c -- passed back to the main WEPP model to account for erosion c -- possibly taking place. c wmelt(iplane) = wmelt(iplane) + hrrain(hour) hrmlt(hour,iplane) = wmelt(iplane) c c Write hourly winter output if(snoflg.eq.1) then tmpvr2=hrsnow(hour)*1000.0 tmpvr3=hrrain(hour)*1000.0 tmpvr4=gdrft*1000.0 tmpvr5=fdrft*1000.0 tmpvr6=wmelt(iplane)*1000.0 tmpvr7=snodpt(iplane)*1000.0 tmpvr8=frdp(iplane)*1000.0 tmpvr9=thdp(iplane)*1000.0 write(42,'(2(i3,1x),i5,1x,9(f6.2,1x),i3,1x,i3)') 1 sdate,hour,year-ibyear,tmpvr2,tmpvr3,tmpvr4,tmpvr5,tmpvr6, 2 tmpvr7,densgt,tmpvr8,tmpvr9,cycle(iplane),iplane endif totmel = totmel + wmelt(iplane) if(hour.eq.24) temp24=hrtemp if(isflag.eq.1.and.snodpy(iplane).le.0.0) then temphr=24-hour endif c if(hour.lt.24) goto 500 c c END OF 500 HOURLY LOOP c ------------------------------------------------------------ wmelt(iplane)=totmel c return end