PROGRAM main c Program to generate four different ensemble deterministic c precipitation forecasts. The forecasts are based on the c Ebert (2001) study. This program uses four ensemble members. c This program is for accumulated precipitation integer gridx, gridy, t, majct, pct, mct, valid parameter (gridx=99, gridy=92, t=5, dm=15) integer is(gridx*gridy), js(gridx*gridy), aindex(gridx*gridy), +findex(gridx*gridy*15), rindex(gridx*gridy) double precision r real prcp(gridx,gridy,t,40), inppre(gridx,gridy,t), a, b, +samp(15), median, major, inprcp(gridx,gridy), +pmat(gridx*gridy*15), pmean(gridx*gridy), pran(gridx*gridy), +enspm(gridx,gridy,t), ranpm(gridx,gridy,t), +old_prcp(gridx,gridy,33,40) integer i,j,p,q,memcount,k,f, rct, m character chr character*25 fname(18), ofname(9) c array of filenames fname(1)='MOD1.dat' fname(2)='MOD2.dat' fname(3)='MOD3.dat' fname(4)='MOD4.dat' fname(5)='MOD5.dat' fname(6)='MOD6.dat' fname(7)='MOD7.dat' fname(8)='MOD8.dat' fname(9)='MOD9.dat' fname(10)='MOD10.dat' fname(11)='MOD11.dat' fname(12)='MOD12.dat' fname(13)='MOD13.dat' fname(14)='MOD14.dat' fname(15)='MOD15.dat' c read in precip file and convert to six hour totals do m=1,dm do k=1,33 open(21,file=fname(m),form='unformatted', +status='old',access='direct',recl=dx*dy*4) read(21) ((old_prcp(i,j,k,m),i=1,gridx),j=1,gridy) enddo close(21) enddo c get 6 hour totals from 3 hour totals do m=1,dm kk1=4 kk2=5 kk3=6 kk4=7 kk5=8 kk6=9 do k=1,t do j=1,gridy do i=1,gridx prcp(i,j,k,m)=old_prcp(i,j,kk1,m)+old_prcp(i,j,kk2,m)+ + old_prcp(i,j,kk3,m)+old_prcp(i,j,kk4,m)+ + old_prcp(i,j,kk5,m)+old_prcp(i,j,kk6,m) enddo enddo kk1=kk1+6 kk2=kk2+6 kk3=kk3+6 kk4=kk4+6 kk5=kk5+6 kk6=kk6+6 enddo enddo c random seed f=time() call sgrnd(f) r=grnd() c now compute the four deterministic forecasts c 1. ensemble mean c 2. ensemble median c 3. majority rules (average nonzero QPFs if at least c half of models forecast precip) c 4. probability matching do k=1,t valid=0 rct=0 do j=1,gridy do i=1,gridx if (prcp(i,j,1,1)<1000 + .AND. prcp(i,j,1,2)<1000 + .AND. prcp(i,j,1,3)<1000 + .AND. prcp(i,j,1,4)<1000 + .AND. prcp(i,j,1,5)<1000 + .AND. prcp(i,j,1,6)<1000 + .AND. prcp(i,j,1,7)<1000 + .AND. prcp(i,j,1,8)<1000 + .AND. prcp(i,j,1,9)<1000 + .AND. prcp(i,j,1,10)<1000 + .AND. prcp(i,j,1,11)<1000 + .AND. prcp(i,j,1,12)<1000 + .AND. prcp(i,j,1,13)<1000 + .AND. prcp(i,j,1,14)<1000 + .AND. prcp(i,j,1,15)<1000) then valid=valid+1 do rr=1,15 samp(rr)=prcp(i,j,k,rr) enddo do q=1,15 if (samp(q)<0) then samp(q)=0.0 end if enddo prcp(i,j,k,17)=(samp(1)+samp(2)+samp(3)+ + samp(4)+samp(5)+samp(6)+samp(7)+samp(8)+samp(9)+samp(10)+ + samp(11)+samp(12)+samp(13)+samp(14)+samp(15))/15.0 pmat(valid*15-14)=samp(1) pmat(valid*15-13)=samp(2) pmat(valid*15-12)=samp(3) pmat(valid*15-11)=samp(4) pmat(valid*15-10)=samp(5) pmat(valid*15-9)=samp(6) pmat(valid*15-8)=samp(7) pmat(valid*15-7)=samp(8) pmat(valid*15-6)=samp(9) pmat(valid*15-5)=samp(10) pmat(valid*15-4)=samp(11) pmat(valid*15-3)=samp(12) pmat(valid*15-2)=samp(13) pmat(valid*15-1)=samp(14) pmat(valid*15)=samp(15) major=0.0 majct=0 call inssort(samp,15) prcp(i,j,k,18)=(samp(4)+samp(5))/2 if (samp(5) > 0.254) then major=samp(1)+samp(2)+samp(3)+samp(4) majct=4 do q=1,4 if (samp(q) > 0.254) then major=major+samp(q) majct=majct+1 endif enddo else major=0.0 majct=1 end if prcp(i,j,k,19)=major/majct is(valid)=i js(valid)=j pmean(valid)=prcp(i,j,k,17) else prcp(i,j,k,17)=9.999E20 prcp(i,j,k,18)=9.999E20 prcp(i,j,k,19)=9.999E20 prcp(i,j,k,20)=9.999E20 ranpm(i,j,k)=9.999E20 end if enddo enddo c sort the data for probability matching call sortrx(valid,pmean,aindex) call sortrx(valid*15,pmat,findex) c random probability match sample do j=1,valid r=grnd() pran(j)=pmat(dint(r*valid*15)+1) if (pran(j)/25.4 > 0.01) then rct=rct+1 end if enddo call sortrx(valid,pran,rindex) c put the sorted data into the prob matching array do p=1,valid prcp(is(aindex(p)),js(aindex(p)),k,20)= + pran(rindex(p)) enddo c if (k .eq. 4) then c open(unit=19,file='pmatch.dat',form='unformatted', c + recl=4,access='direct') c open(unit=18,file='pmean.dat',form='unformatted', c + recl=4,access='direct') c do p=1,valid c write(19,rec=p*4-3) pmat(p*4-3) c write(19,rec=p*4-2) pmat(p*4-2) c write(19,rec=p*4-1) pmat(p*4-1) c write(19,rec=p*4) pmat(p*4) c write(18,rec=p) pmean(p) c enddo c close(19) c close(18) c endif enddo c output prcp(17 and 18) c output results to binary file do k=1,t write(30) ((prcp(i,j,k,20),i=1,gridx),j=1,gridy) c write(31) ((prcp(i,j,k,17),i=1,gridx),j=1,gridy) enddo end c a simple insertion sort subroutine inssort(data,n) integer i,j,n real data(n), index do i=2,n index = data(i) j = i do while ((j > 0) .AND. (data(j-1) > index)) data(j) = data(j-1) j = j - 1 enddo data(j) = index enddo end ************************************************************************ subroutine sgrnd(seed) * implicit integer(a-z) * * Period parameters parameter(N = 624) * dimension mt(0:N-1) * the array for the state vector common /block/mti,mt save /block/ * * setting initial seeds to mt[N] using * the generator Line 25 of Table 1 in * [KNUTH 1981, The Art of Computer Programming * Vol. 2 (2nd Ed.), pp102] * mt(0)= iand(seed,-1) do 1000 mti=1,N-1 mt(mti) = iand(69069 * mt(mti-1),-1) 1000 continue * return end ************************************************************************ double precision function grnd() * implicit integer(a-z) * * Period parameters parameter(N = 624) parameter(N1 = N+1) parameter(M = 397) parameter(MATA = -1727483681) * constant vector a parameter(UMASK = -2147483648) * most significant w-r bits parameter(LMASK = 2147483647) * least significant r bits * Tempering parameters parameter(TMASKB= -1658038656) parameter(TMASKC= -272236544) * dimension mt(0:N-1) * the array for the state vector common /block/mti,mt save /block/ data mti/N1/ * mti==N+1 means mt[N] is not initialized * dimension mag01(0:1) data mag01/0, MATA/ save mag01 * mag01(x) = x * MATA for x=0,1 * TSHFTU(y)=ishft(y,-11) TSHFTS(y)=ishft(y,7) TSHFTT(y)=ishft(y,15) TSHFTL(y)=ishft(y,-18) * if(mti.ge.N) then * generate N words at one time if(mti.eq.N+1) then * if sgrnd() has not been called, call sgrnd(4357) * a default initial seed is used endif * do 1000 kk=0,N-M-1 y=ior(iand(mt(kk),UMASK),iand(mt(kk+1),LMASK)) mt(kk)=ieor(ieor(mt(kk+M),ishft(y,-1)),mag01(iand(y,1))) 1000 continue do 1100 kk=N-M,N-2 y=ior(iand(mt(kk),UMASK),iand(mt(kk+1),LMASK)) mt(kk)=ieor(ieor(mt(kk+(M-N)),ishft(y,-1)),mag01(iand(y,1))) 1100 continue y=ior(iand(mt(N-1),UMASK),iand(mt(0),LMASK)) mt(N-1)=ieor(ieor(mt(M-1),ishft(y,-1)),mag01(iand(y,1))) mti = 0 endif * y=mt(mti) mti=mti+1 y=ieor(y,TSHFTU(y)) y=ieor(y,iand(TSHFTS(y),TMASKB)) y=ieor(y,iand(TSHFTT(y),TMASKC)) y=ieor(y,TSHFTL(y)) * if(y.lt.0) then grnd=(dble(y)+2.0d0**32)/(2.0d0**32-1.0d0) else grnd=dble(y)/(2.0d0**32-1.0d0) endif * return end C From Leonard J. Moss of SLAC: C Here's a hybrid QuickSort I wrote a number of years ago. It's C based on suggestions in Knuth, Volume 3, and performs much better C than a pure QuickSort on short or partially ordered input arrays. SUBROUTINE SORTRX(N,DATA,INDEX) C=================================================================== C C SORTRX -- SORT, Real input, indeX output C C C Input: N INTEGER C DATA REAL C C Output: INDEX INTEGER (DIMENSION N) C C This routine performs an in-memory sort of the first N elements of C array DATA, returning into array INDEX the indices of elements of C DATA arranged in ascending order. Thus, C C DATA(INDEX(1)) will be the smallest number in array DATA; C DATA(INDEX(N)) will be the largest number in DATA. C C The original data is not physically rearranged. The original order C of equal input values is not necessarily preserved. C C=================================================================== C C SORTRX uses a hybrid QuickSort algorithm, based on several C suggestions in Knuth, Volume 3, Section 5.2.2. In particular, the C "pivot key" [my term] for dividing each subsequence is chosen to be C the median of the first, last, and middle values of the subsequence; C and the QuickSort is cut off when a subsequence has 9 or fewer C elements, and a straight insertion sort of the entire array is done C at the end. The result is comparable to a pure insertion sort for C very short arrays, and very fast for very large arrays (of order 12 C micro-sec/element on the 3081K for arrays of 10K elements). It is C also not subject to the poor performance of the pure QuickSort on C partially ordered data. C C Created: 15 Jul 1986 Len Moss C C=================================================================== INTEGER N,INDEX(N) REAL DATA(N) INTEGER LSTK(31),RSTK(31),ISTK INTEGER L,R,I,J,P,INDEXP,INDEXT REAL DATAP C QuickSort Cutoff C C Quit QuickSort-ing when a subsequence contains M or fewer C elements and finish off at end with straight insertion sort. C According to Knuth, V.3, the optimum value of M is around 9. INTEGER M PARAMETER (M=9) C=================================================================== C C Make initial guess for INDEX DO 50 I=1,N INDEX(I)=I 50 CONTINUE C If array is short, skip QuickSort and go directly to C the straight insertion sort. IF (N.LE.M) GOTO 900 C=================================================================== C C QuickSort C C The "Qn:"s correspond roughly to steps in Algorithm Q, C Knuth, V.3, PP.116-117, modified to select the median C of the first, last, and middle elements as the "pivot C key" (in Knuth's notation, "K"). Also modified to leave C data in place and produce an INDEX array. To simplify C comments, let DATA[I]=DATA(INDEX(I)). C Q1: Initialize ISTK=0 L=1 R=N 200 CONTINUE C Q2: Sort the subsequence DATA[L]..DATA[R]. C C At this point, DATA[l] <= DATA[m] <= DATA[r] for all l < L, C r > R, and L <= m <= R. (First time through, there is no C DATA for l < L or r > R.) I=L J=R C Q2.5: Select pivot key C C Let the pivot, P, be the midpoint of this subsequence, C P=(L+R)/2; then rearrange INDEX(L), INDEX(P), and INDEX(R) C so the corresponding DATA values are in increasing order. C The pivot key, DATAP, is then DATA[P]. P=(L+R)/2 INDEXP=INDEX(P) DATAP=DATA(INDEXP) IF (DATA(INDEX(L)) .GT. DATAP) THEN INDEX(P)=INDEX(L) INDEX(L)=INDEXP INDEXP=INDEX(P) DATAP=DATA(INDEXP) ENDIF IF (DATAP .GT. DATA(INDEX(R))) THEN IF (DATA(INDEX(L)) .GT. DATA(INDEX(R))) THEN INDEX(P)=INDEX(L) INDEX(L)=INDEX(R) ELSE INDEX(P)=INDEX(R) ENDIF INDEX(R)=INDEXP INDEXP=INDEX(P) DATAP=DATA(INDEXP) ENDIF C Now we swap values between the right and left sides and/or C move DATAP until all smaller values are on the left and all C larger values are on the right. Neither the left or right C side will be internally ordered yet; however, DATAP will be C in its final position. 300 CONTINUE C Q3: Search for datum on left >= DATAP C C At this point, DATA[L] <= DATAP. We can therefore start scanning C up from L, looking for a value >= DATAP (this scan is guaranteed C to terminate since we initially placed DATAP near the middle of C the subsequence). I=I+1 IF (DATA(INDEX(I)).LT.DATAP) GOTO 300 400 CONTINUE C Q4: Search for datum on right <= DATAP C C At this point, DATA[R] >= DATAP. We can therefore start scanning C down from R, looking for a value <= DATAP (this scan is guaranteed C to terminate since we initially placed DATAP near the middle of C the subsequence). J=J-1 IF (DATA(INDEX(J)).GT.DATAP) GOTO 400 C Q5: Have the two scans collided? IF (I.LT.J) THEN C Q6: No, interchange DATA[I] <--> DATA[J] and continue INDEXT=INDEX(I) INDEX(I)=INDEX(J) INDEX(J)=INDEXT GOTO 300 ELSE C Q7: Yes, select next subsequence to sort C C At this point, I >= J and DATA[l] <= DATA[I] == DATAP <= DATA[r], C for all L <= l < I and J < r <= R. If both subsequences are C more than M elements long, push the longer one on the stack and C go back to QuickSort the shorter; if only one is more than M C elements long, go back and QuickSort it; otherwise, pop a C subsequence off the stack and QuickSort it. IF (R-J .GE. I-L .AND. I-L .GT. M) THEN ISTK=ISTK+1 LSTK(ISTK)=J+1 RSTK(ISTK)=R R=I-1 ELSE IF (I-L .GT. R-J .AND. R-J .GT. M) THEN ISTK=ISTK+1 LSTK(ISTK)=L RSTK(ISTK)=I-1 L=J+1 ELSE IF (R-J .GT. M) THEN L=J+1 ELSE IF (I-L .GT. M) THEN R=I-1 ELSE C Q8: Pop the stack, or terminate QuickSort if empty IF (ISTK.LT.1) GOTO 900 L=LSTK(ISTK) R=RSTK(ISTK) ISTK=ISTK-1 ENDIF GOTO 200 ENDIF 900 CONTINUE C=================================================================== C C Q9: Straight Insertion sort DO 950 I=2,N IF (DATA(INDEX(I-1)) .GT. DATA(INDEX(I))) THEN INDEXP=INDEX(I) DATAP=DATA(INDEXP) P=I-1 920 CONTINUE INDEX(P+1) = INDEX(P) P=P-1 IF (P.GT.0) THEN IF (DATA(INDEX(P)).GT.DATAP) GOTO 920 ENDIF INDEX(P+1) = INDEXP ENDIF 950 CONTINUE C=================================================================== C C All done END