subroutine ranset() 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 include 'crandom.inc' c read: mox, thresh, thres2, vv c write: ranary,g_dimi,g_dimp,g_dsum,g_ssum c 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 real last_r(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/ 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 + + + FUNCTION DECLARATIONS + + + real randn, dstn1 c + + + LOCAL VARIABLES + + + integer i,j real x2sum, g_davg c if(mox.ne.2) then dimi = dim(mox) else if(nt.eq.1) then dimi = 28 else dimi = 29 endif endif c c ------ update the days to end of this month. g_dimi(mox) = g_dimi(mox) + dimi do 999 j=1, 9 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 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 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) g_dimp(mox) = g_dimp(mox) - ldimp endif if((level .gt. thresh(j)) .or. (level2 .gt. thres2(j))) goto 10 999 continue c return end