C ***************************************************************** C PROGRAM: spi.daily.f (spiprg) C C LANGUAGE: Fortran 77 C C MACHINE: VM C ----------------------------------------------------------------* C PURPOSE: Generate 3, 6, 9, 12, 24-month averaged SPI C ----------------------------------------------------------------* C USAGE: C Define environment variables FORT10, FORT50, FORT60 for Input and Output files C spi.daily.x C ----------------------------------------------------------------* C INPUT FILES: C work/spi/spi_date - target date in yyyymmdd format C work/spi/monthly.600m - 600-month daily precipitation data C OUTPUT FILES: C work/spi/spi3.us.tmp - 3-month averaged SPI C work/spi/spi6.us.tmp - 6-month averaged SPI C work/spi/spi9.us.tmp - 9-month averaged SPI C work/spi/spi12.us.tmp - 12-month averaged SPI C work/spi/spi24.us.tmp - 24-month averaged SPI C ----------------------------------------------------------------* C SUBROUTINES USED: C spigam - Calculate indices assuming incomplete gamma distribution. C gamfit - Estimate incomplete gamma parameters C FUNCTIONS USED: C anvnrm - Return z for the input probability C gamcdf - Compute probability of a<=x using incomplete gamma parameters C gaminv - Compute inverse gamma function C gamser - Functions for the incomplete gamma functions P and Q C gammcf - Evaluate P(a,x) in its continued fraction representation C gammap - Evaluate the incomplete gamma function P(a,x), choosing the most C appropriate representation. C gammaq - Evaluate the incomplete gamma function Q(a,x), choosing the most C appropriate representation. C gammln - ln(gamma) function C ----------------------------------------------------------------* C INPUT VARIABLES: C C LOCAL VARIABLES: C C ----------------------------------------------------------------* C AUTHOR(S): Kingtse Mo C ----------------------------------------------------------------* C DATE 2000-01-01 C ----------------------------------------------------------------* C MODIFICATIONS: C 2018-03-01 Hai-Tien Lee R2O standard compliance C ***************************************************************** program spiprg parameter (nlen=5, maxyrs=50, amssng=-9999.0,jm=60, 1 im=150,ibegyr=600, iendyr=ibegyr+maxyrs) character header*80 dimension len(nlen) real prec(maxyrs*12), beta(12), gamm(12), pzero(12), 1 mask(150,60),uu(150,60) 1 ,spi(maxyrs*12,nlen),p(150,60,ibegyr),a(900),r(150,60) real b(150,60,ibegyr),pp(6) data len / 3, 6, 9, 12, 24/ open(10,file='work/spi/spi_date') open(11,file='work/spi/monthly.600m', c access='direct',recl=150*60*4,form='unformatted') open(20,file='work/spi/spi3.us.tmp',form='unformatted', c access='direct',recl=150*60*4) open(21,file='work/spi/spi6.us.tmp',form='unformatted', c access='direct',recl=150*60*4) open(22,file='work/spi/spi9.us.tmp',form='unformatted', c access='direct',recl=150*60*4) open(23,file='work/spi/spi12.us.tmp',form='unformatted', c access='direct',recl=150*60*4) open(24,file='work/spi/spi24.us.tmp',form='unformatted', c access='direct',recl=150*60*4) c***** c read(10,*) mdate print *,mdate myr=mdate/10000 mmo=(mdate-myr*10000)/100 tmonth=(myr-1950)*12+mmo c Open files. These are defaults under most unix f77's print *,' myr ',myr,' mo ',mmo,' total ',tmonth undef=-999. do j=1,jm do i=1,im uu(i,j)=-999. enddo enddo itest=200 print *,' start to get mask' read(11,rec=itest) r do j=1,jm do i=1,im mask(i,j)=0. if(r(i,j).ge.0.) mask(i,j)=1 enddo enddo c write(31) mask irec=2 do it=1,ibegyr read(11,rec=it+irec) r do jj=1,jm do ii=1,im p(ii,jj,it)=r(ii,jj) enddo enddo print *,it, p(70,25,it) enddo limit=ibegyr print *,'total x points',it,limit c Compute SPI's print *,' computer SPI for 3 month means' do 10 mm=1,5 print *,' start to compute ',mm iunit=19+mm do 6000 jj=1,jm do 6000 ii=1,im if(mask(ii,jj).ne.0.) then do it=1,limit prec(it)=p(ii,jj,it) if(prec(it).gt.70.) prec(it)=70. enddo c print *,prec(70) c do 10 i = 1,3 mi=mm call spigam (len(mm), prec, beta, gamm, pzero, spi(1,mm)) c 10 continue c Skip leading missings ifirst = 0 20 continue ifirst = ifirst + 1 if (spi(ifirst,mi) .eq. amssng) goto 20 c Skip trailing missngs last = maxyrs*12 + 1 30 continue last = last - 1 if (spi(last,mi) .eq. amssng) goto 30 c Ouput SPI's c print *,' output spi ' c do 40 i = ifirst, last do 40 i = ifirst, ibegyr if (spi(i,mi) .eq. amssng) goto 50 iy = (i-1)/12 + 1 imm = mod (i-1,12)+1 c print 1001, iy, imm, spi(i,mi) b(ii,jj,i)=spi(i,mi) 40 continue 50 continue endif 6000 continue print *,' p values' it=limit do j=1,jm do i=1,im r(i,j)=undef if(mask(i,j).eq.1.) then r(i,j)=b(i,j,it) endif enddo enddo print *,r(80,25),it write(iunit,rec=1) r C for output last 35year spi6 if ( len(mm).eq. 6 ) then open(33,file='work/spi/spi6.35y',form='unformatted', c access='direct',recl=150*60*4) i=1 do k=15*12-1,600-12,12 write(33,rec=i)b(:,:,k) i=i+1 end do close(33) end if c enddo 10 continue stop 1000 format(a80) 1001 format(i4,i3,4f7.2) end cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c Read monthly prec data. c c Format: c Header c Data- yyyy mm prec c c Where: c yyyy - year; values > ENDYR and < BEGYR will be skipped. c mm - month [1-12] c prec - precipitation (in 0.01's) c c Special codes: c -9900 = Missing c ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c These functions compute the Standardized Precipitation Index c using an incomplete gamma distribution function to estimate c probabilities. c c Useful references are: c c _Numerical Recipes in C_ by Flannery, Teukolsky and Vetterling c Cambridge University Press, ISBN 0-521-35465-x c c _Handbook of Mathematical Functions_ by Abramowitz and Stegun c Dover, Standard Book Number 486-61272-4 c c Notes for ForTran users: c - There is a companion _Numerical Recipes in Fortran. The c following code was translated from the c code. c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c Calculate indices assuming incomplete gamma distribution. c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc subroutine spigam (nrun, pp, beta, gamm, pzero, index) parameter (nlen=5, maxyrs=50, amssng=-9999.0, 1 ibegyr=600, iendyr=ibegyr+maxyrs) real pp(maxyrs*12), index(maxyrs*12),tmparr(maxyrs+1), 1 beta(12), gamm(12), pzero(12) c c The first nrun-1 index values will be missing. c do 10 j = 1, nrun-1 index(j) = amssng 10 continue c c Sum nrun precip. values; c store them in the appropriate index location. c c If any value is missing; set the sum to missing. c do 30 j = nrun, maxyrs*12 index(j) = 0.0 do 20 i = 0, nrun-1 if(pp(j - i) .ne. amssng) then index(j) = index(j) + pp(j - i) else index(j) = amssng goto 30 endif 20 continue 30 continue c c For nrun<12, the monthly distributions will be substantially c different. So we need to compute gamma parameters for c each month starting with the (nrun-1)th. c do 50 i = 0,11 n = 0 do 40 j = nrun+i, maxyrs*12, 12 if(index(j) .ne. amssng) then n = n + 1 tmparr(n) = index(j) endif 40 continue im = mod (nrun+i-1, 12) + 1 c c Here's where we do the fitting. c call gamfit (tmparr, n, alpha, beta(im), gamm(im), pzero(im)) 50 continue c c Replace precip. sums stored in index with SPI's c do 60 j = nrun, maxyrs*12 im = mod (j-1,12) + 1 if(index(j) .ne. amssng) then c c Get the probability c index(j) = gamcdf(beta(im), gamm(im), pzero(im), index(j)) c c Convert prob. to z value. c index(j) = anvnrm(index(j)) endif 60 continue return end ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c input prob; return z. c c See Abromowitz and Stegun _Handbook of Mathematical Functions_, p. 933 c ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc function anvnrm (prob) data c0, c1, c2 /2.515517, 0.802853, 0.010328/ data d1, d2, d3 /1.432788, 0.189269, 0.001308/ if (prob .gt. 0.5) then sign = 1.0 prob = 1.0 - prob else sign = -1.0 endif if (prob .lt. 0.0) then write(0, *) 'Error in anvnrm(). Prob. not in [0,1.0]' anvnrm = 0.0 return endif if (prob .eq. 0.0) then anvnrm = 1.0e37 * sign return endif t = sqrt(alog (1.0 / (prob * prob))) anvnrm = (sign * (t - ((((c2 * t) + c1) * t) + c0) / 1 ((((((d3 * t) + d2) * t) + d1) * t) + 1.0))) return end cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c Estimate incomplete gamma parameters. c c Input: c datarr - data array c n - size of datarr c c Output: c alpha, beta, gamma - gamma paarameters c pzero - probability of zero. c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc subroutine gamfit (datarr, n, alpha, beta, gamm, pzero) real datarr(*) if (n .le. 0) then write(0, *) 'Error in gamfit - empty data array' stop endif sum = 0.0 sumlog = 0.0 pzero = 0.0 nact = 0 c compute sums do 10 i = 1, n if (datarr(i) .gt. 0.0) then sum = sum + datarr(i) sumlog = sumlog + alog (datarr(i)) nact = nact + 1 else pzero = pzero + 1 endif 10 continue pzero = pzero / n if(nact .ne. 0.0) av = sum / nact c Bogus data array but do something reasonable if(nact .eq. 1) then alpha = 0.0 gamm = 1.0 beta = av return endif c They were all zeroes. if(pzero .eq. 1.0) then alpha = 0.0 gamm = 1.0 beta = av return endif c Use MLE alpha = alog (av) - sumlog / nact gamm = (1.0 + sqrt (1.0 + 4.0 * alpha / 3.0)) / (4.0 * alpha) beta = av / gamm return end cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c Compute probability of a<=x using incomplete gamma parameters. c c Input: c beta, gamma - gamma parameters c pzero - probability of zero. c x - value. c c Return: c Probability a<=x. c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc function gamcdf (beta, gamm, pzero, x) if(x .le. 0.0) then gamcdf = pzero else gamcdf = pzero + (1.0 - pzero) * gammap (gamm, x / beta) endif return end cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c Compute inverse gamma function i.e. return x given p where CDF(x) = p. c c Input: c beta, gamma - gamma parameters c pzero - probability of zero. c prob - probability. c c Return: c x as above. c c Method: c We use a simple binary search to first bracket out initial c guess and then to refine our guesses until two guesses are within c tolerance (eps). Is there a better way to do this? c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc function gaminv (beta, gamm, pzero, prob) data eps /1.0e-7/ c Check if prob < prob of zero if (prob .le. pzero) then gaminv = 0.0 return endif c Otherwise adjust prob prob = (prob - pzero) / (1.0 - pzero) c Make initial guess. Keep doubling estimate until prob is c bracketed. thigh = 2.0*eps 10 continue phigh = gamcdf (beta, gamm, pzero, thigh) if(phigh .ge. prob) goto 20 thigh = thigh*2.0 goto 10 20 continue tlow = thigh / 2.0 c Iterate to find root. niter = 0 30 continue if((thigh - tlow) .le. eps) goto 40 niter = niter + 1 t = (tlow + thigh) / 2.0 p = gamcdf (beta, gamm, pzero, t) if (p .lt. prob) then tlow = t else thigh = t endif goto 30 40 continue gaminv = (tlow + thigh) / 2.0 return end ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c Functions for the incomplete gamma functions P and Q c c 1 /x -t a-1 c P (a, x) = -------- | e t dt, a > 0 c Gamma(x)/ 0 c c Q (a, x) = 1 - P (a, x) c c Reference: Press, Flannery, Teukolsky, and Vetterling, c _Numerical Recipes_, pp. 160-163 c c Thanks to kenny@cs.uiuc.edu c ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c Evaluate P(a,x) by its series representation. c function gamser (a, x) c Maximum number of iterations, and bound on error. parameter (maxitr=100, eps=3.0e-7) data iwarn /0/ gln = gammln (a) if (x .eq. 0.0) then gamser = 0.0 return endif ap = a sum = 1.0 / a del = sum do 10 n = 1, maxitr ap = ap + 1.0 del = del * (x / ap) sum = sum + del if (abs (del) .lt. eps * abs (sum)) goto 20 10 continue iwarn = iwarn + 1 if (iwarn .lt. 20) then write (0, *) 'gamser(',a,x,'): not converging.' write (0, *) 'Approximate value of ',sum,' + /-',del,' used.' endif 20 continue gamser = sum * exp (-x + a * alog (x) - gln) return end c c Evaluate P(a,x) in its continued fraction representation. c function gammcf (a, x) parameter (maxitr=100, eps=3.0e-7) data nwarn / 0 /, g / 0.0 / gln = gammln (a) gold = 0.0 a0 = 1.0 a1 = x b0 = 0.0 b1 = 1.0 fac = 1.0 do 10 n = 1, maxitr an = n ana = an - a a0 = (a1 + a0 * ana) * fac b0 = (b1 + b0 * ana) * fac anf = an * fac a1 = x * a0 + anf * a1 b1 = x * b0 + anf * b1 if (a1 .ne. 0.0) then fac = 1.0 / a1 g = b1 * fac if (abs((g - gold) / g) .lt. eps) goto 20 gold = g endif 10 continue nwarn = nwarn + 1 if (nwarn .lt. 20) then write (0, *) 'gammcf(',a,x,'): not converging.' write (0, *) 'Inaccurate value of ', g, ' +/- ', 1 abs(g - gold), ' used.' endif 20 continue gammcf = g * exp (-x + a * alog (x) - gln) return end c c Evaluate the incomplete gamma function P(a,x), choosing the most c appropriate representation. c function gammap (a, x) if (x .lt. a + 1.0) then gammap = gamser (a, x) else gammap = 1.0 - gammcf (a, x) endif return end c c Evaluate the incomplete gamma function Q(a,x), choosing the most c appropriate representation. c function gammaq (a, x) if (x .lt. a + 1.0) then gammaq = 1.0 - gamser (a, x) else gammaq = gammcf (a, x) endif return end c c For those who don't have a ln(gamma) function. c function gammln(xx) dimension cof(6) data cof /76.18009173, -86.50532033, 24.01409822, -1.231739516, 1 0.120858003e-2, -0.536382e-5/ x = xx - 1.0 tmp = x + 5.5 tmp = tmp - (x+0.5) * alog (tmp) ser = 1.0 do 10 j = 1, 5 x = x + 1.0 ser = ser + cof(j) / x 10 continue gammln = -tmp + alog (2.50662827465 * ser) return end