subroutine clgen(ntd,iyear) c c + + + PURPOSE + + + c Simulates daily solar radiation, simulates daily precipitation c and/or maximum and minimum air temperature at the users option, c and calls functions RANDN and DSTN1. c c + + + ARGUMENT DECLARATIONS + + + integer ntd, iyear c c + + + ARGUMENT DEFINITIONS + + + c ntd - days in this year (365 or 366) c c + + + COMMON BLOCKS + + + include 'cbk1.inc' c read: rh c modify: tdp c rh - Avg Monthly Dew Point Temperature. Used to calculate TDP. c tdp - Generated dewpoint temperature (C). c include 'cbk3.inc' c read: ida c ida - Julian Day of Year. Used as a subscript to R. c include 'cbk4.inc' c read: mo, nt c mo - The current month (1=Jan, 2=Feb...). c nt - Set to 1 if IYEAR is not a leap year: otherwise, zero. c include 'cbk7.inc' c read: prw,obmx,obmn,obsl,k1,k2,k3,k4,k5,k9,yls,ylc,pit,nsim,msim, c stdtx,stdtm c modify: rst,v1,v3,v5,v7,v11,ra,tmxg,tmng,rmx,l c prw(1,12) - monthly probability of wet day after wet day c prw(2,12) - monthly probability of wet day after dry day c obmx - Maximum Temperature. c obmn - Minimum Temperature. c k1 ... k9 - Seeds for random number generation. c yls - ??? -- Used to compute CH and YS. sin(ylt/clt) sin(latitude)c ylc - ??? -- Used to compute CH and YC. cos(ylt/clt) cos(latitude)c pit - ??? -- Used to compute SD. Defined as pit=58.13 c nsim - ??? Has value zero or one. Used as a Switch. c msim - ??? Has value zero or one. Used as a Switch. c stdtx - Standard deviation of daily max. temp. for the month. c stdtm - Standard deviation of daily min. temp. for the month. c rst(i,j) - Array Of Monthly precipitation stats. c dim1: month 1..12 c dim2: 1=mean of daily rainfall c mean liquid equivalent precipitation c depth (inches) for a day precipitation c occurs (by month) [=avg total precip c for month / # wet days in month] c 2=std deviation of daily rainfall c standard deviation of the daily precip c value (inches) (by month) c 3=skew coefficient of daily rainfall c c v1 ... v11 - Random numbers used to generate various daily values. c ra - Generated radiation? RADG receives RA's value for output. c tmxg - Generated daily max temp. c tmng - Generated daily min temp. c rmx - Maximum possible solar radiation. c l - Set to either 1 or 2; linked to nsim 0 or 1; selects PRW. c include 'cbk5.inc' c modify: r c r - Daily Precipitation amount (inches of water) c include 'cinterp.inc' c read: lf, rf, o_mo c lf - weighting factor for the midpoint value on this month's end c of the time interval. c rf - weighting factor for the midpoint value on the "other" end. c o_mo - month (on the "other" end) whose average value should be used. c CC include 'ctap.inc' CC include 'ctap2.inc' c write: tap1, tap2, tap3, tap4, tap5, tap6 c include 'crandom.inc' c read: vvx,v2x,v4x,v6x,v8x,v12x c modify: mox,dax c c + + + FUNCTION DECLARATIONS + + + real dstn1, randn, ryf2, fouri2 c + + + LOCAL VARIABLES + + + real rx, tmpv10, tap5, v6, tmpv13, tmpvr9, tmpvr7, tmpvr8 real tmpvr6, tap4, tap2, v4, tap1, v2, v12, xx, tmpvr2 real tmpvr1, xlv, r6, tmpvr3, v8, tap3, h, ch, yc, ys real sd, xi c c + + + END SPECIFICATIONS + + + c xi=ida sd=0.4102*sin((xi-80.25)/pit) ch=-yls*tan(sd)/ylc c if (ch.ge.1.0) then h=0.0 else if (ch.le.-1.0) then h=3.1416 else h=acos(ch) endif c ys=yls*sin(sd) C -- XXX -- Note that YC is also used to store a Y/N user response! C CRM -- 10/21/99 yc=ylc*cos(sd) c ---- max possible solar radiation for this day of the year. rmx=711.0*(h*ys+yc*sin(h)) c c ---- Generate a month's worth of Consecutive Random Numbers for c each parameter. (All at once -- not interleaved! -- CRM) if(mo .ne. mox) then mox = mo dax = 1 c call ranset() call ranset(ntd,iyear) else dax = dax + 1 endif c if (nsim.ne.0) then vv=vvx(dax) if ((prw(l,mo).le.0.0).or.(vv.gt.prw(l,mo))) then r(ida)=0.0 l=2 tap3=0.0 else c -------- A Mutated variant of Equation 2.1.5 c -------- Generated Daily Precip c According to Larry M. Younkin, the amount of precip is c Assumed to follow a Pearson type III c distribution: c c 2s / / g / g \ \3 \ c R = Rbar + --- | | --- | x - --- | + 1 | - 1 | c g \ \ 6 \ 6 / / / c v8=v8x(dax) C if(rst(mo,3).eq.0.0) rst(mo,3)=0.01 C r6=rst(mo,3)/6.0 c ---------(interpolate precip skewness) if(interp.eq.0) then tmpvr3 = rst(mo,3) else if(interp .eq. 1) then tmpvr3 = rst(mo,3)*lf + rst(o_mo,3)*rf else if(interp .eq. 2) then tmpvr3 = fouri2(3) else if(interp .eq. 3) then tmpvr3 = ryf2(mo,dax,ntd,3) endif if(tmpvr3.eq.0.0) tmpvr3=0.01 r6=tmpvr3/6.0 c c -- XXX -- Bandaid: -- CRM 4/26/2000 if(v7 .eq. 0.0) v7 = randn(k5) xlv=(dstn1(v7,v8)-r6)*r6+1.0 xlv=(xlv**3-1.0)*2.0/tmpvr3 tap3=xlv v7=v8 c ---------(interpolate precip std. dev. & mean) if(interp.eq.0) then tmpvr1 = rst(mo,1) tmpvr2 = rst(mo,2) else if(interp .eq. 1) then tmpvr1 = rst(mo,1)*lf + rst(o_mo,1)*rf tmpvr2 = rst(mo,2)*lf + rst(o_mo,2)*rf else if(interp .eq. 2) then tmpvr1 = fouri2(1) tmpvr2 = fouri2(2) else if(interp .eq. 3) then tmpvr1 = ryf2(mo,dax,ntd,1) tmpvr2 = ryf2(mo,dax,ntd,2) endif c r(ida)=xlv*rst(mo,2)+rst(mo,1) r(ida)=xlv*tmpvr2 + tmpvr1 if (r(ida).lt.0.01) r(ida)=0.01 l=1 endif endif c if(msim.eq.0) then xx=1. v12=v12x(dax) tdp=xx*dstn1(v11,v12) v11=v12 else xx=1. v2=v2x(dax) tmxg=xx*dstn1(v1,v2) tap1=tmxg c ------ Shift 2nd constant to 1st position for "reuse". v1=v2 v4=v4x(dax) tmng=xx*dstn1(v3,v4) tap2=tmng v12=v12x(dax) tdp =xx*dstn1(v11,v12) tap4=tdp c ------ Shift 2nd constant to 1st position for "reuse". v3=v4 v11=v12 c c -------Equations 2.1.10 & 2.1.11 C ------(interpolate max & min temp means & SD's) if(interp.eq.0) then tmpvr6 = obmx(mo) tmpvr8 = stdtx(mo) tmpvr7 = obmn(mo) tmpvr9 = stdtm(mo) else if(interp .eq. 1) then tmpvr6 = obmx(mo)*lf + obmx(o_mo)*rf tmpvr8 = stdtx(mo)*lf + stdtx(o_mo)*rf tmpvr7 = obmn(mo)*lf + obmn(o_mo)*rf tmpvr9 = stdtm(mo)*lf + stdtm(o_mo)*rf else if(interp .eq. 2) then tmpvr6 = fouri2(6) tmpvr8 = fouri2(8) tmpvr7 = fouri2(7) tmpvr9 = fouri2(9) else if(interp .eq. 3) then tmpvr6 = ryf2(mo,dax,ntd,6) tmpvr8 = ryf2(mo,dax,ntd,8) tmpvr7 = ryf2(mo,dax,ntd,7) tmpvr9 = ryf2(mo,dax,ntd,9) endif c c ------ Generated Max & Min Daily Temps. C tmxg=obmx(mo)+tmxg*stdtx(mo) C tmng=obmn(mo)+tmng*stdtm(mo) tmxg=tmpvr6 + tmxg*tmpvr8 tmng=tmpvr7 + tmng*tmpvr9 c --- RANGE CHECK on MIN TEMP: if (tmng.gt.tmxg) tmng=tmxg-0.2*abs(tmxg) endif c c TDP now calculated using standard dev. instead of CV. 3/95 c ---- Generated Daily Dew Point Temperature. c ---------(interpolate humidity) if(interp.eq.0) then tmpv13 = rh(mo) else if(interp .eq. 1) then tmpv13 = rh(mo)*lf + rh(o_mo)*rf else if(interp .eq. 2) then tmpv13 = fouri2(13) else if(interp .eq. 3) then tmpv13 = ryf2(mo,dax,ntd,13) endif c -----A mutated version of Equation 2.1.14: c -- XXX -- CRM -- 10/21/99 -- I think there is a mis-placed paren here. c tdp =rh(mo)+(tdp*(stdtx(mo)+stdtm(mo)/2.)) cc I believe it should be as follows: cc tdp =rh(mo)+(tdp*(stdtx(mo)+stdtm(mo))/2.) tdp = tmpv13+tdp*(tmpvr8+tmpvr9)/2.0 CC tdp = tmpv13+(tdp*(tmpvr8+tmpvr9/2.0)) if (tdp.gt.((tmxg+tmng)/2.)) tdp=((tmxg+tmng)/2.)*0.99 c Change to limit - dew point in cold months ADNe if(tdp.lt.-10.) tdp=1.1*tmng c v6=v6x(dax) c ---- Compute Daily Radiation. ra=xx*dstn1(v5,v6) tap5=ra c ---------(interpolate radiation) if(interp.eq.0) then tmpv10 = obsl(mo) else if(interp .eq. 1) then tmpv10 = obsl(mo)*lf + obsl(o_mo)*rf else if(interp .eq. 2) then tmpv10 = fouri2(10) else if(interp .eq. 3) then tmpv10 = ryf2(mo,dax,ntd,10) endif c -----A mutation of Equations 2.1.12 & 2.1.13: (Mis-placed parens, this time in WEPP User Doc ?) C rx=rmx-obsl(mo) C if (obsl(mo).gt.rx) rx=obsl(mo) C ra=obsl(mo)+ra*rx/4.0 rx=rmx-tmpv10 if (tmpv10.gt.rx) rx=tmpv10 ra=tmpv10+ra*rx/4.0 if(ra.ge.rmx) ra=0.90*rmx if (ra.le.0.0) ra=0.05*rmx c ------ Shift 2nd constant to 1st position for "reuse". v5=v6 c return end