subroutine ranset(ntd,iyear) c ---- Purpose: Generate Array of Random Numbers, an entire month at a c time, for each parameter. Ensure that for the run to c date, the numbers will produce standard normal deviates c with a mean at the specified confidence level ("thresh") c and a varience at the specified confidence level ("thres2"). c C. R. Meyer -- 4/6-24/2000 and 5/16-17/2000. c c Modified: Was using NT to determine whether the current year is a c Leap Year; however, NT only gets set for options 4 & 7, c and only for the *initial year*. Modified to use NTD. c C. R. Meyer -- 3/7/2001 c c + + + ARGUMENT DECLARATIONS + + + integer ntd,iyear c c + + + ARGUMENT DEFINITIONS + + + c ntd - days in this year (365 or 366) c include 'crandom.inc' c read: mox, thresh, thres2, vv c write: ranary,g_dimi,g_dimp,g_dsum,g_ssum c RanAry elements: c ranary(1,1) -- precip prob c ranary(1,2) -- max temp c ranary(1,3) -- min temp c ranary(1,4) -- radiation c ranary(1,5) -- precip amount c ranary(1,6) -- wind dir c ranary(1,7) -- wind vel c ranary(1,8) -- TDP c ranary(1,9) -- TP c C -- XXX -- NT is no longer used. Probably can eliminate this include: C CRM -- 3/20/01. include 'cbk4.inc' c read: nt c nt - Set to 1 if IYEAR is not a leap year: otherwise, zero c include 'cbk7.inc' c modify: k1-k6, k8-k10, v1-v12 c c ---- dim -- Days in each Month. integer dim(12), dimi, ldimp, ell, ellx character*15 params(9) real last_r(9), lst_rx(9) data dim/31,28,31,30,31,30,31,31,30,31,30,31/ c ---- ell -- ell=1 ==> Precip; ell=2 ==> No Precip. data ell/2/ data last_r/9* -1.0/ data params/"Prob. of Precip", "Max. Temp.", "Min. Temp.", & "Radiation", "Precip. Amt.", "Wind Dir.", "Wind Vel.", & "Temp. Dew Pt.", "Time to Peak"/ save ell, last_r c c real ransum,ranavg,level real ransum,level,level2 c ransum -- sum of this month's random std norm deviates c x2sum -- sum of S^2 (ie, Xi^2) for this month's random std norm deviates c ranavg -- avg of this month's random std norm deviates c level -- level at which we're confident of a difference in the means. c level2 -- level at which we're confident of a difference in the S^2's. c g_dsum -- Sum of std norm deviates, for this parameter, for this month, c from the beginning of the simulation. c g_ssum -- Sum of S^2 of std norm deviates, for this parameter, for this c month, from the beginning of the simulation. c g_dimi -- Total days for this month of the year, from the beginning c of simulation to end of current month. c g_dimp -- Total days with precip for this month of the year, from the c beginning of simulation to end of current month. c g_davg -- avg to date for this parm's random std norm deviates, for c this month of the year. c ldimp -- Local variable for days in the month with precip. c last_r -- Value of last random deviate generated for this parm. c ellx -- Saves the value of "ell" in case we re-do preip amts. c lst_rx -- Saves the value of "last_r(j)" in case of a re-do c + + + FUNCTION DECLARATIONS + + + real randn, dstn1 c + + + LOCAL VARIABLES + + + integer i,j integer iredo real x2sum, g_davg c if(mox.ne.2) then dimi = dim(mox) else C if(nt.eq.1) then if(ntd.eq.365) then dimi = 28 else dimi = 29 endif endif c c ------ save the value of "ell" case we re-do precip amount. ellx = ell c ------ prepare to count number of re-do's. iredo = 0 c c ------ update the days to end of this month. g_dimi(mox) = g_dimi(mox) + dimi do 999 j=1, 9 c ------ save "last_r(j)" in case we have a re-do. lst_rx(j) = last_r(j) c c ------ Sum up the standard normal deviates for this month. 10 continue ldimp = 0 ransum=0.0 x2sum =0.0 c c ------ BEGIN random value generation for the current parm. c (Each parm uses its own instance of the random number generator.) do 998 i=1, dimi if(j .eq. 1) then ranary(i,j) = randn(k1) else if(j .eq. 2) then ranary(i,j) = randn(k2) else if(j .eq. 3) then ranary(i,j) = randn(k3) else if(j .eq. 4) then ranary(i,j) = randn(k4) else if(j .eq. 5) then c ---------- we had precip today. if(ranary(i,1) .le. prw(ell,mox)) then ranary(i,j) = randn(k5) ell = 1 ldimp = ldimp + 1 c ---------- we didnt have precip today. else ranary(i,j) = 0.0 ell = 2 endif else if(j .eq. 6) then ranary(i,j) = randn(k6) else if(j .eq. 7) then ranary(i,j) = randn(k8) else if(j .eq. 8) then ranary(i,j) = randn(k9) else if(j .eq. 9) then ranary(i,j) = randn(k10) endif c ------ END random value generation for the current parm. c c ------ BEGIN summing of Random Normal Deviates for the current parm. & month c -------- its not the 1st day of the run. if(last_r(j) .ne. -1.0) then if(j .ne. 5 .or. ell .eq. 1) then ransum = ransum + dstn1(last_r(j),ranary(i,j)) x2sum = x2sum + dstn1(last_r(j),ranary(i,j))**2 last_r(j) = ranary(i,j) endif c c -------- it is the 1st day of the run. else if(j .eq. 1) then if(vv .eq. 0.0) vv = randn(k1) last_r(j) = vv else if(j .eq. 2) then last_r(j) = v1 else if(j .eq. 3) then last_r(j) = v3 else if(j .eq. 4) then last_r(j) = v5 c ---------- we have precip. else if(j .eq. 5 .and. ell .eq. 1) then if(v7 .eq. 0.0) v7 = randn(k5) last_r(j) = v7 C -- XXX -- Huh? -- CRM -- 3/7/01: C ldimp = 1 else if(j .eq. 6) then if(fx .eq. 0.0) fx = randn(k6) last_r(j) = fx else if(j .eq. 7) then last_r(j) = v9 else if(j .eq. 8) then last_r(j) = v11 else if(j .eq. 9) then if(z .eq. 0.0) z = randn(k10) last_r(j) = z endif if(j .ne. 5 .or. ell .eq. 1) then ransum = ransum + dstn1(last_r(j),ranary(i,j)) x2sum = x2sum + dstn1(last_r(j),ranary(i,j))**2 last_r(j) = ranary(i,j) endif endif c ------ END summing of Random Normal Deviates for the current parm. & month c 998 continue c c ------ Update the Grand Sum of the standard normal deviates for this month. g_dsum(j,mox) = g_dsum(j,mox) + ransum g_ssum(j,mox) = g_ssum(j,mox) + x2sum c ------ Compute the Grand Average of the Std Norm Dev's for this month. if(j .ne. 5) then g_davg = g_dsum(j,mox) / float(g_dimi(mox)) call conflm(g_davg,g_dimi(mox),0.0,1.0,level) call confls(real(g_ssum(j,mox)),g_dimi(mox),level2) else g_dimp(mox) = g_dimp(mox) + ldimp if(g_dimp(mox) .gt. 0) then g_davg = g_dsum(j,mox) / float(g_dimp(mox)) else g_davg = 0.0 endif call conflm(g_davg,g_dimp(mox),0.0,1.0,level) call confls(real(g_ssum(j,mox)),g_dimp(mox),level2) endif c ------ If the result is throws us outside the confidence limits, RE-DO. if((level .gt. thresh(j)) .or. (level2 .gt. thres2(j))) then c------- Count the number of "re-do's" and exit after the prescribed number. iredo = iredo + 1 if(iredo .gt. 9995) then if(j .ne. 5) then write(*,*) "ParmNr Yr Mon Mlevel SDlevel S2sum DiM:", & j,iyear,mox,level,level2,g_ssum(j,mox),g_dimi(mox) else write(*,*) "ParmNr Yr Mon Mlevel SDlevel S2sum DiMP:", & j,iyear,mox,level,level2,g_ssum(j,mox),g_dimp(mox) endif if(iredo .eq. 10000) then write(*, '("*** ERROR *** Could not produce desired", & " level of quality in",/20x, "<< ", a, " >>", & " random deviates.")') params(j) endif endif CC if(j.ne.5) then CC write(*,*) mox,j,level,level2,g_ssum(j,mox),g_dimi(mox), CC 1 thresh(j),thres2(j) CC else CC write(*,*) mox,j,level,level2,g_ssum(j,mox),g_dimp(mox), CC 1 thresh(j),thres2(j) CC endif CC write(*,*) " ********** REDO ***********" g_dsum(j,mox) = g_dsum(j,mox) - ransum g_ssum(j,mox) = g_ssum(j,mox) - x2sum if(j .eq. 5) then g_dimp(mox) = g_dimp(mox) - ldimp ell = ellx endif last_r(j) = lst_rx(j) endif c if((level .gt. thresh(j)) .or. (level2 .gt. thres2(j))) goto 10 if(((level .gt. thresh(j)) .or. (level2 .gt. thres2(j))) & .and. iredo.lt.10000) goto 10 999 continue c return end