subroutine spinew(iendyr,ndiv) c Subroutine to calculate divisional SPIs for different time scalles c Original program written by Ned Guttman,Nov 1997; some of the software c and logic based on software provided by John Kleist, Colorado State Univ. c Extensively modified to run as a subroutine in shared space, Trevor Wallis, c Apr 2016. c ---------------------------------------------------------------------- c CAVEAT c SPI values in central part of the precip distribution can c be considered to be accurate and reasonable if the number c of independent non-zero accumulated precipitation amounts c is greater than about 40; values in the tails of the c distribution can be considered accurate and reasonable c if the number of non-zero accumulated amounts is greater c than about 60. COMMON /COMPATHS/ PATHCPC,PATHDRD,PATHDRDW,cMeanpath,PATHWDATA, : cOdata,cIndata CHARACTER PATHCPC*42,PATHDRD*34,PATHDRDW*40,cMeanpath*36, : PATHWDATA*41,cIndata*38,cOdata*38,cNl*2,cScode*2, : cPcode*2,cAcode*2,CODE*2,cStfl*3 integer nrun(7),num(12),numpos(12) real pp(2400),probne(2400),pcpacc(2400),pzero(12),spi(2400), : prec(200,14),spiout(13),prob(13),pacc(13),pcpn(200,14) double precision par1(12),par2(12),par3(12),xmom1(12), * xmom2(12),xmom3(12),bmssng data nrun /1,2,3,6,9,12,24/ c if (ndiv .eq. 0) then WRITE (9,*) 'subroutine spinew started for states' open(2,file=cIndata//'us-stregs.spi',status='old') else WRITE (9,*) 'subroutine spinew started' open(2,file=cIndata//'usdivs.spi',status='old') endif ibegyr = 1895 nlen = 7 iout = 10 amssng = -9.99 bmssng = -9.99d00 c ams1=-99.99 ams2=-9.999 ams3=-9.99 c Open output files if (ndiv .eq. 0) then cStfl = 'st.' else cStfl = 'dv.' endif do 100,nl=1,7 write(cNl,'(i2.2)') nrun(nl) write(cScode,'(i2.2)') nl+70 write(cPcode,'(i2.2)') nl+80 write(cAcode,'(i2.2)') nl+90 open(nl*10+1,file=cOdata//'spi/scale'//cNl//'_l-moments.'// : cStfl,status='unknown') open(nl*10+2,file=cOdata//'spi/sp'//cNl//cScode//cStfl, : status='old') open(nl*10+3,file=cOdata//'spi/pr'//cNl//cPcode//cStfl, : status='old') open(nl*10+4,file=cOdata//'spi/ac'//cNl//cAcode//cStfl, : status='old') 100 continue c Divisional data starts at 1895 maxyrs=iendyr-1894 MAXY=maxyrs*12 c Step through usdivs.spi ist1=0 idv1=0 1 read (2,345,end=50) CODE,ist,idv 345 format(A2,I3.3,I2.2) if (ist.eq.ist1.and.idv.eq.idv1) go to 11 call GETPOR(ist,idv,1,pcpn) c Process data in=1 do 221 k=1,200 iy=INT(pcpn(k,14)) if (iy.lt.ibegyr) go to 221 do m=1,12 prec(in,m)=pcpn(k,m) enddo in=in+1 221 continue c create 1-dimension array pp() from 2-dimension array prec() ik=0 do 13 j=1,in-1 do k=1,12 pp(k+ik)=prec(j,k) enddo ik=ik+12 13 continue c compute SPI values, probabilities, accumulated precip do 10 i=1,nlen ih1=iout*i+1 ih2=iout*i+2 ih3=iout*i+3 ih4=iout*i+4 ispi= 70+i iprob=80+i ipacc=90+i call spipe3 (nrun(i),pp,par1,par2,par3,pzero,spi,probne, * pcpacc,xmom1,xmom2,xmom3,amssng,ams1,ams2,ams3,bmssng, * num,numpos,maxyrs) c write out monthly arrays of ncounts, lmoments, parameters,etc. do k=1,12 write(ih1,1000)ist,idv,nrun(i),k,xmom1(k),xmom2(k), * xmom3(k),par1(k),par2(k),par3(k),num(k),numpos(k),pzero(k) 1000 format (i3.3,i2.2,2i3,6f13.8,2i5,f6.3) enddo c---------------------------------------------------------------------- c Calculate SPIs, probs, etc and write to output files do im=1,13 spiout(im)=ams1 prob(im) =ams2 pacc(im) =ams3 enddo do 40 j=1,MAXY iy=(j-1)/12+ibegyr im=mod (j-1,12)+1 spiout(im)=spi(j) prob(im)=probne(j) pacc(im)=pcpacc(j) if (im.eq.12.or.j.eq.MAXY) then c write out arrays write (ih2,1001) CODE,ist,idv,ispi,iy,(spiout(im),im=1,12),ams1 write (ih3,1003) CODE,ist,idv,iprob,iy,(prob(im),im=1,12),ams2 write (ih4,1001) CODE,ist,idv,ipacc,iy,(pacc(im),im=1,12),ams3 1001 format (A2,i3.3,i2.2,i2,i4,13f7.2) 1003 format (A2,i3.3,i2.2,i2,i4,13f7.3) c Reset arrays 30 do im=1,13 spiout(im)=ams1 prob(im) =ams2 pacc(im) =ams3 enddo endif 40 continue 10 continue 11 ist1=ist idv1=idv goto 1 50 continue c Close output files close (2) do i=1,nlen do k=1,4 ih=iout*i+k close (ih) enddo enddo WRITE (9,*) 'Subroutine spinew executed' c return end c------------------------------------------------------------------------ SUBROUTINE SPIPE3(nrun,pp,par1,par2,par3,pzero,spi,probne, * pcpacc,xmom1,xmom2,xmom3,amssng,ams1,ams2,ams3,bmssng, * num,numpos,maxyrs) c Subroutine to calculate L-Moments and parameters for Pearson III c distribution c Developed by Ned Guttman, Nov 1997, based on programs by John Kleist, c Colorado Stateniv., that computed spi values from a 2 parameter ML gamma c distribution. Subroutines and functions from Jon Hosking, IBM. c ------------------------------------------------------------------ dimension pp(2400),probne(2400),pcpacc(2400), * pzero(12),spi(2400),num(12),numpos(12) real*8 par1(12),par2(12),par3(12),index(2400), * xmom1(12),xmom2(12),xmom3(12),xmom(3),para(3), * y(2400),x(200),temparr(201),bmssng,uu,vv c MAXY=maxyrs*12 c the first nrun-1 index values will be missing do 10 j=1,nrun-1 index(j) =bmssng probne(j)=ams2 pcpacc(j)=ams3 spi(j) =ams1 10 continue c sum nrun precip values; store them in appropriate index location c if any value is missing, set the sum to missing do 30 j=nrun,MAXY index(j)=0.0D00 do 20 i=0,nrun-1 if(pp(j-i).ne.amssng)then index(j)=index(j)+DBLE(pp(j-i)) pcpacc(j)=index(j) else index(j) =bmssng spi(j) =ams1 probne(j)=ams2 pcpacc(j)=ams3 goto 30 endif 20 continue 30 continue c to maintain seasonality, do everything by month do 50 i=0,11 n=0 nz=0 np=0 if (nrun.le.12) then do 40 j=nrun+i,MAXY,12 if(index(j).ge.0.0D00) then c n-count for all non-missing data, including zeroes n=n+1 if(index(j).gt.0.0D00) then c np-count for all non-missing, non-zero data; temparr is for c non-zero, non-missing precipitation np=np+1 temparr(np)=index(j) else c nz-count for all non-missing, zero data nz=nz+1 endif endif 40 continue elseif (nrun.gt.12) then j=nrun+i c step the sequence of data by the length of nrun to get independent samples if (nrun.gt.12.and.nrun.le.24) nstep=24 if (nrun.gt.24.and.nrun.le.36) nstep=36 if (nrun.gt.36.and.nrun.le.48) nstep=48 if (nrun.gt.48.and.nrun.le.60) nstep=60 if (nrun.gt.60.and.nrun.le.72) nstep=72 41 if (j.gt.MAXY) GOTO 42 c if(index(j).ne.bmssng) then n=n+1 if(index(j).gt.0.0D00) then np=np+1 temparr(np)=index(j) else nz=nz+1 endif j=j+nstep goto 41 else c look for next non-missing nrun>12 value j=j+12 goto 41 endif endif 42 im=mod(nrun+i-1,12)+1 num(im)=n numpos(im)=np if (n.eq.0) go to 50 pzero(im)=float(nz)/float(n) c order the data call sort(temparr,x,np,maxyrs) c fit l-moments for non-zero precip; if 2nd moment=0, routine c fails and all 3 moments set to zero (condition indicates all c data values are equal). also, if not enough non-zero data c (less than 3), routine fails and all moments set to zero. call samlmr(x,np,xmom,3,-0.00D0,0.0D0,ifail) xmom1(im)=xmom(1) xmom2(im)=xmom(2) xmom3(im)=xmom(3) c compute parameters of Pearson type III (3-parameter gamma); the c three parameters are the mean, s.d., skewness of the NON-ZERO data. c The mean for all the data is (1-pzero)*par1 or (1-pzero)*xmom1 c since par1=xmom1 call pelpe3(xmom,para,ifail) c if ifail not = 0, then all 3 parameters are set to zero c if routine fails, output is sent to unit 6 by pelpe3 routine par1(im)=para(1) par2(im)=para(2) par3(im)=para(3) 50 continue c compute the probability (cdfpe3), take into account the c mixed distribution if pzero not equal to zero, truncate c the probability from .001 to .999, transform the probability to c spi (quastn). ifail=1 indicates invalid parameter 2 (negative). c if routines fail, output is sent to unit 6 by called routines. do 60 j=nrun,MAXY im=mod(j-1,12)+1 if (num(im).eq.0) go to 61 para(1)=par1(im) para(2)=par2(im) para(3)=par3(im) c set missing values to amssng 61 if(index(j).lt.0.0D00) then probne(j)=ams2 spi(j) =ams1 pcpacc(j)=ams3 goto 60 endif if(index(j).ge.0.0D00) then uu=index(j) call CDFPE3(uu,para,vv,ifail,bmssng) index(j)=vv c if cdf routine fails, set probne and spi = amssng if(ifail.ne.0) then index(j)=bmssng goto 61 endif y(j)=DBLE(pzero(im))+(1.0D00-DBLE(pzero(im)))*index(j) c force the probabilities and therefore the spi to be bounded by +/- 3.09 if (y(j).gt.0.999D00) y(j)=0.999D00 if (y(j).lt.0.001D00) y(j)=0.001D00 probne(j)=y(j) uu=y(j) call quastn(uu,vv) spi(j)=vv endif 60 continue return end c ---------------------------------------------------------------------------- SUBROUTINE CDFPE3(X,PARA,vv,ifail,bmssng) C* FORTRAN CODE WRITTEN FOR INCLUSION IN IBM RESEARCH REPORT RC20525, C* FORTRAN ROUTINES FOR USE WITH THE METHOD OF L-MOMENTS, VERSION 3 C* J. R. M. HOSKING C* IBM RESEARCH DIVISION C* T. J. WATSON RESEARCH CENTER C* YORKTOWN HEIGHTS C* NEW YORK 10598, U.S. C* VERSION 3 AUGUST 1996 C C DISTRIBUTION FUNCTION OF THE PEARSON TYPE 3 DISTRIBUTION IMPLICIT DOUBLE PRECISION (A-H,O-Z) real*8 PARA(3),vv DATA ZERO/.0D0/,HALF/0.5D0/,ONE/1.0D0/,TWO/2.0D0/,FOUR/4.0D0/ DATA RTHALF/0.70710 67811 86547 524D0/ C SMALL IS USED TO TEST WHETHER SKEWNESS IS EFFECTIVELY ZERO DATA SMALL/1.0D-6/ C ifail=0 vv=ZERO IF(PARA(2).LE.ZERO)GOTO 1000 GAMMA=PARA(3) IF(DABS(GAMMA).LE.SMALL)GOTO 10 ALPHA=FOUR/(GAMMA*GAMMA) Z=TWO*(X-PARA(1))/(PARA(2)*GAMMA)+ALPHA IF(Z.GT.ZERO) then call DLGAMA(ALPHA,DLG) call GAMIND(Z,ALPHA,DLG,GAM) vv=GAM endif IF(GAMMA.LT.ZERO) vv=ONE-vv c RETURN C C ZERO SKEWNESS C 10 Z=(X-PARA(1))/PARA(2) z=Z*RTHALF call DERF(z,DD) vv=HALF+HALF*DD RETURN C 1000 WRITE(6,*) ' ERROR *** ROUTINE CDFPE3 : PARAMETERS INVALID' ifail=1 vv=bmssng RETURN END c c -------------------------------------------------------------------------- SUBROUTINE PELPE3(XMOM,PARA,ifail) C* FORTRAN CODE WRITTEN FOR INCLUSION IN IBM RESEARCH REPORT RC20525, C* FORTRAN ROUTINES FOR USE WITH THE METHOD OF L-MOMENTS, VERSION 3 C* J. R. M. HOSKING C* IBM RESEARCH DIVISION C* T. J. WATSON RESEARCH CENTER C* YORKTOWN HEIGHTS C* NEW YORK 10598, U.S. C* VERSION 3 AUGUST 1996 C C PARAMETER ESTIMATION VIA L-MOMENTS FOR THE PEARSON TYPE 3 DISTRIBUTION C METHOD: RATIONAL APPROXIMATION IS USED TO EXPRESS ALPHA, THE SHAPE C PARAMETER OF THE GAMMA DISTRIBUTION, AS A FUNCTION OF TAU-3. C RELATIVE ACCURACY OF THE APPROXIMATION IS BETTER THAN 3E-5. IMPLICIT DOUBLE PRECISION (A-H,O-Z) real*8 XMOM(3),PARA(3) DATA ZERO/.0D0/,THIRD/0.33333333D0/,HALF/0.5D0/,ONE/1.0D0/, * TWO/2.0D0/ C SMALL IS USED TO TEST WHETHER SKEWNESS IS EFFECTIVELY ZERO DATA SMALL/1D-6/ C CONSTANTS USED IN MINIMAX APPROXIMATIONS DATA C1,C2,C3/ 0.2906D0, 0.1882D0, 0.0442D0/ DATA D1,D2,D3/ 0.36067D0,-0.59567D0, 0.25361D0/ DATA D4,D5,D6/-2.78861D0, 2.56096D0,-0.77045D0/ DATA PI3,ROOTPI/9.4247780D0,1.7724539D0/ C ifail=0 T3=DABS(XMOM(3)) IF(XMOM(2).LE.ZERO.OR.T3.GE.ONE)GOTO 1000 IF(T3.LE.SMALL)GOTO 100 IF(T3.GE.THIRD)GOTO 10 T=PI3*T3*T3 ALPHA=(ONE+C1*T)/(T*(ONE+T*(C2+T*C3))) GOTO 20 10 CONTINUE T=ONE-T3 ALPHA=T*(D1+T*(D2+T*D3))/(ONE+T*(D4+T*(D5+T*D6))) 20 CONTINUE RTALPH=DSQRT(ALPHA) ah=ALPHA+HALF call DLGAMA(ALPHA,dlg1) call DLGAMA(ah,dlg2) BETA=ROOTPI*XMOM(2)*DEXP(dlg1-dlg2) PARA(1)=XMOM(1) PARA(2)=BETA*RTALPH PARA(3)=TWO/RTALPH IF(XMOM(3).LT.ZERO)PARA(3)=-PARA(3) RETURN C ZERO SKEWNESS 100 CONTINUE PARA(1)=XMOM(1) PARA(2)=XMOM(2)*ROOTPI PARA(3)=ZERO RETURN c Failure 1000 DO 1010 I=1,3 PARA(I)=ZERO 1010 continue ifail=1 RETURN END c c -------------------------------------------------------------------------- SUBROUTINE SAMLMR(X,N,XMOM,NMOM,A,B,ifail) C* FORTRAN CODE WRITTEN FOR INCLUSION IN IBM RESEARCH REPORT RC20525, C* FORTRAN ROUTINES FOR USE WITH THE METHOD OF L-MOMENTS, VERSION 3 C* J. R. M. HOSKING C* IBM RESEARCH DIVISION C* T. J. WATSON RESEARCH CENTER C* YORKTOWN HEIGHTS C* NEW YORK 10598, U.S. C* VERSION 3 AUGUST 1996 C C SAMPLE L-MOMENTS OF A DATA ARRAY C FOR UNBIASED ESTIMATES (OF THE LAMBDA'S) SET A=B=ZERO. OTHERWISE, C PLOTTING-POSITION ESTIMATORS ARE USED, BASED ON THE PLOTTING POSITION C (J+A)/(N+B) FOR THE J'TH SMALLEST OF N OBSERVATIONS. FOR EXAMPLE, C A=-0.35D0 AND B=0.0D0 YIELDS THE ESTIMATORS RECOMMENDED BY C HOSKING ET AL. (1985, TECHNOMETRICS) FOR THE GEV DISTRIBUTION. C IMPLICIT DOUBLE PRECISION (A-H,O-Z) real*8 X(N),XMOM(NMOM),SUM(20) DATA ZERO/0.0D0/,ONE/1.0D0/ c ifail=0 IF(NMOM.GT.20.OR.NMOM.GT.N)GOTO 1000 DO 10 J=1,NMOM SUM(J)=ZERO 10 continue IF(A.EQ.ZERO.AND.B.EQ.ZERO)GOTO 50 IF(A.LE.-ONE.OR.A.GE.B)GOTO 1010 C PLOTTING-POSITION ESTIMATES OF PWM'S DO 30 I=1,N PPOS=(I+A)/(N+B) TERM=X(I) SUM(1)=SUM(1)+TERM DO 20 J=2,NMOM TERM=TERM*PPOS SUM(J)=SUM(J)+TERM 20 continue 30 CONTINUE DO 40 J=1,NMOM SUM(J)=SUM(J)/N 40 continue GOTO 100 C UNBIASED ESTIMATES OF PWM'S 50 DO 70 I=1,N Z=I TERM=X(I) SUM(1)=SUM(1)+TERM DO 60 J=2,NMOM Z=Z-ONE TERM=TERM*Z SUM(J)=SUM(J)+TERM 60 continue 70 CONTINUE Y=N Z=N SUM(1)=SUM(1)/Z DO 80 J=2,NMOM Y=Y-ONE Z=Z*Y SUM(J)=SUM(J)/Z 80 continue C L-MOMENTS 100 K=NMOM P0=ONE IF(NMOM-NMOM/2*2.EQ.1)P0=-ONE DO 120 KK=2,NMOM AK=K P0=-P0 P=P0 TEMP=P*SUM(1) DO 110 I=1,K-1 AI=I P=-P*(AK+AI-ONE)*(AK-AI)/(AI*AI) TEMP=TEMP+P*SUM(I+1) 110 continue SUM(K)=TEMP K=K-1 120 continue XMOM(1)=SUM(1) IF(NMOM.EQ.1)RETURN XMOM(2)=SUM(2) IF(SUM(2).EQ.ZERO)GOTO 1020 IF(NMOM.EQ.2)RETURN DO 130 K=3,NMOM XMOM(K)=SUM(K)/SUM(2) 130 continue RETURN C SET RESULTS IF PROCESS HAS FAILED 1000 ifail=1 do i=1,3 xmom(i)=0 enddo WRITE(9,*) 'ERROR ROUTINE SAMLMR : PARAMETER NMOM INVALID' RETURN 1010 WRITE(9,*) 'ERROR ROUTINE SAMLMR : '// : 'PLOTTING-POSITION PARAMETERS INVALID' RETURN 1020 ifail=1 do i=1,3 xmom(i)=0 enddo WRITE(9,*) 'ERROR ROUTINE SAMLMR : ALL DATA VALUES EQUAL' RETURN END c c --------------------------------------------------------------------------- SUBROUTINE DERF(X,DD) C* FORTRAN CODE WRITTEN FOR INCLUSION IN IBM RESEARCH REPORT RC20525, C* FORTRAN ROUTINES FOR USE WITH THE METHOD OF L-MOMENTS, VERSION 3 C* J. R. M. HOSKING C* IBM RESEARCH DIVISION C* T. J. WATSON RESEARCH CENTER C* YORKTOWN HEIGHTS C* NEW YORK 10598, U.S. C* VERSION 3 AUGUST 1996 C C ERROR FUNCTION C BASED ON ALGORITHM 5666, J.F.HART ET AL. (1968) 'COMPUTER C APPROXIMATIONS' C ACCURATE TO 15 DECIMAL PLACES C IMPLICIT DOUBLE PRECISION (A-H, O-Z) DATA ZERO/0D0/,ONE/1D0/,TWO/2D0/,THREE/3D0/,FOUR/4D0/,P65/0.65D0/ C COEFFICIENTS OF RATIONAL-FUNCTION APPROXIMATION DATA P0,P1,P2,P3,P4,P5,P6/ * 0.22020 68679 12376 1D3, 0.22121 35961 69931 1D3, * 0.11207 92914 97870 9D3, 0.33912 86607 83830 0D2, * 0.63739 62203 53165 0D1, 0.70038 30644 43688 1D0, * 0.35262 49659 98910 9D-1/ DATA Q0,Q1,Q2,Q3,Q4,Q5,Q6,Q7/ * 0.44041 37358 24752 2D3, 0.79382 65125 19948 4D3, * 0.63733 36333 78831 1D3, 0.29656 42487 79673 7D3, * 0.86780 73220 29460 8D2, 0.16064 17757 92069 5D2, * 0.17556 67163 18264 2D1, 0.88388 34764 83184 4D-1/ C C1 IS SQRT(2), C2 IS SQRT(2/PI) C BIG IS THE POINT AT WHICH DERF=1 TO MACHINE PRECISION DATA C1/1.4142 13562 37309 5D0/ DATA C2/7.9788 45608 02865 4D-1/ DATA BIG/6.25D0/,CRIT/5.0D0/ C DD=ZERO IF(X.EQ.ZERO)RETURN XX=DABS(X) IF(XX.GT.BIG)GOTO 20 EXPNTL=DEXP(-X*X) ZZ=DABS(X*C1) IF(XX.GT.CRIT)GOTO 10 DD=EXPNTL*((((((P6*ZZ+P5)*ZZ+P4)*ZZ+P3)*ZZ+P2)*ZZ+P1)*ZZ+P0)/ * (((((((Q7*ZZ+Q6)*ZZ+Q5)*ZZ+Q4)*ZZ+Q3)*ZZ+Q2)*ZZ+Q1)*ZZ+Q0) IF(X.GT.ZERO)DD=ONE-TWO*DD IF(X.LT.ZERO)DD=TWO*DD-ONE RETURN C 10 DD=EXPNTL*C2/(ZZ+ONE/(ZZ+TWO/(ZZ+THREE/(ZZ+FOUR/(ZZ+P65))))) IF(X.GT.ZERO)DD=ONE-DD IF(X.LT.ZERO)DD=DD-ONE RETURN C 20 DD=ONE IF(X.LT.ZERO)DD=-ONE RETURN END c c -------------------------------------------------------------------------- SUBROUTINE DLGAMA(X,DLG) C* FORTRAN CODE WRITTEN FOR INCLUSION IN IBM RESEARCH REPORT RC20525, C* FORTRAN ROUTINES FOR USE WITH THE METHOD OF L-MOMENTS, VERSION 3 C* J. R. M. HOSKING C* IBM RESEARCH DIVISION C* T. J. WATSON RESEARCH CENTER C* YORKTOWN HEIGHTS C* NEW YORK 10598, U.S. C* VERSION 3 AUGUST 1996 C C LOGARITHM OF GAMMA FUNCTION C BASED ON ALGORITHM ACM291, COMMUN. ASSOC. COMPUT. MACH. (1966) C IMPLICIT DOUBLE PRECISION (A-H,O-Z) DATA SMALL,CRIT,BIG,TOOBIG/1.0D-7,13D0,1D9,2D36/ C C0 IS 0.5*LOG(2*PI) C C1...C7 ARE THE COEFFTS OF THE ASYMPTOTIC EXPANSION OF DLGAMA DATA C0,C1,C2,C3,C4,C5,C6,C7/ * 0.91893 85332 04672 742D 0, 0.83333 33333 33333 333D-1, * -0.27777 77777 77777 778D-2, 0.79365 07936 50793 651D-3, * -0.59523 80952 38095 238D-3, 0.84175 08417 50841 751D-3, * -0.19175 26917 52691 753D-2, 0.64102 56410 25641 026D-2/ C S1 IS -(EULER'S CONSTANT), S2 IS PI**2/12 DATA S1/-0.57721 56649 01532 861D 0/ DATA S2/ 0.82246 70334 24113 218D 0/ DATA ZERO/0D0/,HALF/0.5D0/,ONE/1D0/,TWO/2D0/ DLG=ZERO IF(X.LE.ZERO)GOTO 1000 IF(X.GT.TOOBIG)GOTO 1000 C USE SMALL-X APPROXIMATION IF X IS NEAR 0, 1 OR 2 IF(DABS(X-TWO).GT.SMALL)GOTO 10 DLG=DLOG(X-ONE) XX=X-TWO GOTO 20 10 IF(DABS(X-ONE).GT.SMALL)GOTO 30 XX=X-ONE 20 DLG=DLG+XX*(S1+XX*S2) RETURN 30 IF(X.GT.SMALL)GOTO 40 DLG=-DLOG(X)+S1*X RETURN C REDUCE TO DLGAMA(X+N) WHERE X+N.GE.CRIT 40 SUM1=ZERO Y=X IF(Y.GE.CRIT)GOTO 60 Z=ONE 50 Z=Z*Y Y=Y+ONE IF(Y.LT.CRIT)GOTO 50 SUM1=SUM1-DLOG(Z) C USE ASYMPTOTIC EXPANSION IF Y.GE.CRIT 60 SUM1=SUM1+(Y-HALF)*DLOG(Y)-Y+C0 SUM2=ZERO IF(Y.GE.BIG)GOTO 70 Z=ONE/(Y*Y) SUM2=((((((C7*Z+C6)*Z+C5)*Z+C4)*Z+C3)*Z+C2)*Z+C1)/Y 70 DLG=SUM1+SUM2 RETURN C 1000 WRITE(9,*) 'ERROR ROUTINE DLGAMA : ARGUMENT OUT OF RANGE' RETURN END c c ------------------------------------------------------------------------ SUBROUTINE GAMIND(X,ALPHA,G,GAM) C* FORTRAN CODE WRITTEN FOR INCLUSION IN IBM RESEARCH REPORT RC20525, C* FORTRAN ROUTINES FOR USE WITH THE METHOD OF L-MOMENTS, VERSION 3 C* J. R. M. HOSKING C* IBM RESEARCH DIVISION C* T. J. WATSON RESEARCH CENTER C* YORKTOWN HEIGHTS C* NEW YORK 10598, U.S. C* VERSION 3 AUGUST 1996 C C THE INCOMPLETE GAMMA INTEGRAL C BASED ON ALGORITHM AS239, APPL. STATIST. (1988) VOL.37 NO.3 C IMPLICIT DOUBLE PRECISION (A-H,O-Z) DATA ZERO/0D0/,HALF/0.5D0/,ONE/1D0/,TWO/2D0/,THREE/3D0/,X13/13D0/, * X36/36D0/,X42/42D0/,X119/119D0/,X1620/1620D0/,X38880/38880D0/, * RTHALF/0.70710 67811 86547 524D0/ C EPS,MAXIT CONTROL THE TEST FOR CONVERGENCE OF THE SERIES AND C CONTINUED-FRACTION EXPANSIONS. C OFL IS A LARGE NUMBER, USED TO RESCALE THE CONTINUED FRACTION. C UFL IS SUCH THAT EXP(UFL) IS JUST .GT. ZERO. C AHILL CONTROLS THE SWITCH TO HILL'S APPROXIMATION. C DATA EPS/1D-12/,MAXIT/100000/,OFL/1D30/,UFL/-180D0/,AHILL/1D4/ GAM=ZERO IF(ALPHA.LE.ZERO)GOTO 1000 IF(X.LT.ZERO)GOTO 1010 IF(X.EQ.ZERO)RETURN C IF(ALPHA.GT.AHILL)GOTO 100 IF(X.GT.ONE.AND.X.GE.ALPHA)GOTO 50 C SERIES EXPANSION SUM=ONE TERM=ONE A=ALPHA DO 10 IT=1,MAXIT A=A+ONE TERM=TERM*X/A SUM=SUM+TERM IF(TERM.LE.EPS)GOTO 20 10 CONTINUE WRITE(9,*) 'WARNING ** ROUTINE GAMIND :'// : ' ITERATION HAS NOT CONVERGED. RESULT MAY BE UNRELIABLE.' 20 ARG=ALPHA*DLOG(X)-X-G+DLOG(SUM/ALPHA) GAM=ZERO IF(ARG.GE.UFL) GAM=DEXP(ARG) RETURN C CONTINUED-FRACTION EXPANSION 50 CONTINUE A=ONE-ALPHA B=A+X+ONE TERM=ZERO PN1=ONE PN2=X PN3=X+ONE PN4=X*B RATIO=PN3/PN4 DO 70 IT=1,MAXIT A=A+ONE B=B+TWO TERM=TERM+ONE AN=A*TERM PN5=B*PN3-AN*PN1 PN6=B*PN4-AN*PN2 IF(PN6.EQ.ZERO)GOTO 60 RN=PN5/PN6 DIFF=DABS(RATIO-RN) IF(DIFF.LE.EPS.AND.DIFF.LE.EPS*RN)GOTO 80 RATIO=RN 60 PN1=PN3 PN2=PN4 PN3=PN5 PN4=PN6 IF(DABS(PN5).LT.OFL)GOTO 70 PN1=PN1/OFL PN2=PN2/OFL PN3=PN3/OFL PN4=PN4/OFL 70 CONTINUE WRITE(9,*) 'WARNING ** ROUTINE GAMIND :'// : ' ITERATION HAS NOT CONVERGED. RESULT MAY BE UNRELIABLE.' 80 ARG=ALPHA*DLOG(X)-X-G+DLOG(RATIO) GAM=ONE IF(ARG.GE.UFL) GAM=ONE-DEXP(ARG) RETURN C ALPHA IS LARGE: USE HILL'S APPROXIMATION (N.L. JOHNSON AND C S. KOTZ, 1970, 'CONTINUOUS UNIVARIATE DISTRIBUTIONS 1', P.180) C THE 'DO 110' LOOP CALCULATES 2*(X-ALPHA-ALPHA*DLOG(X/ALPHA)), C USING POWER-SERIES EXPANSION TO AVOID ROUNDING ERROR 100 CONTINUE R=ONE/DSQRT(ALPHA) Z=(X-ALPHA)*R TERM=Z*Z SUM=HALF*TERM DO 110 I=1,12 TERM=-TERM*Z*R SUM=SUM+TERM/(I+TWO) IF(DABS(TERM).LT.EPS)GOTO 120 110 CONTINUE 120 WW=TWO*SUM W=DSQRT(WW) IF(X.LT.ALPHA)W=-W H1=ONE/THREE H2=-W/X36 H3=(-WW+X13)/X1620 H4=(X42*WW+X119)*W/X38880 Z=(((H4*R+H3)*R+H2)*R+H1)*R+W z=Z*RTHALF call DERF(z,DD) GAM=HALF+HALF*DD RETURN C ERROR MESSAGES to LOG 1000 WRITE(9,*) 'ERROR ROUTINE GAMIND : SHAPE PARAMETER OUT OF RANGE' RETURN 1010 WRITE(9,*) 'ERROR ROUTINE GAMIND : ARGUMENT OF OUT OF RANGE' RETURN END c c ------------------------------------------------------------------------ SUBROUTINE QUASTN(F,QUA) C* FORTRAN CODE WRITTEN FOR INCLUSION IN IBM RESEARCH REPORT RC20525, C* FORTRAN ROUTINES FOR USE WITH THE METHOD OF L-MOMENTS, VERSION 3 C* J. R. M. HOSKING C* IBM RESEARCH DIVISION C* T. J. WATSON RESEARCH CENTER C* YORKTOWN HEIGHTS C* NEW YORK 10598, U.S. C* VERSION 3 AUGUST 1996 C C QUANTILE FUNCTION OF THE STANDARD NORMAL DISTRIBUTION C BASED ON ALGORITHM AS241, APPL. STATIST. (1988) VOL.37 NO.3 C IMPLICIT DOUBLE PRECISION (A-H,O-Z) DATA ZERO/0.0D0/,HALF/0.5D0/,ONE/1.0D0/ DATA SPLIT1/0.425D0/,SPLIT2/5.0D0/,CONST1/0.180625D0/, * CONST2/1.6D0/ C COEFFICIENTS OF RATIONAL-FUNCTION APPROXIMATIONS DATA A0,A1,A2,A3,A4,A5,A6,A7,B1,B2,B3,B4,B5,B6,B7/ * 0.33871 32872 79636 661D 1, * 0.13314 16678 91784 377D 3, 0.19715 90950 30655 144D 4, * 0.13731 69376 55094 611D 5, 0.45921 95393 15498 715D 5, * 0.67265 77092 70087 009D 5, 0.33430 57558 35881 281D 5, * 0.25090 80928 73012 267D 4, 0.42313 33070 16009 113D 2, * 0.68718 70074 92057 908D 3, 0.53941 96021 42475 111D 4, * 0.21213 79430 15865 959D 5, 0.39307 89580 00927 106D 5, * 0.28729 08573 57219 427D 5, 0.52264 95278 85285 456D 4/ DATA C0,C1,C2,C3,C4,C5,C6,C7,D1,D2,D3,D4,D5,D6,D7/ * 0.14234 37110 74968 358D 1, * 0.46303 37846 15654 530D 1, 0.57694 97221 46069 141D 1, * 0.36478 48324 76320 461D 1, 0.12704 58252 45236 838D 1, * 0.24178 07251 77450 612D 0, 0.22723 84498 92691 846D -1, * 0.77454 50142 78341 408D -3, 0.20531 91626 63775 882D 1, * 0.16763 84830 18380 385D 1, 0.68976 73349 85100 005D 0, * 0.14810 39764 27480 075D 0, 0.15198 66656 36164 572D -1, * 0.54759 38084 99534 495D -3, 0.10507 50071 64441 684D -8/ DATA E0,E1,E2,E3,E4,E5,E6,E7,F1,F2,F3,F4,F5,F6,F7/ * 0.66579 04643 50110 378D 1, * 0.54637 84911 16411 437D 1, 0.17848 26539 91729 133D 1, * 0.29656 05718 28504 891D 0, 0.26532 18952 65761 230D -1, * 0.12426 60947 38807 844D -2, 0.27115 55568 74348 758D -4, * 0.20103 34399 29228 813D -6, 0.59983 22065 55887 938D 0, * 0.13692 98809 22735 805D 0, 0.14875 36129 08506 149D -1, * 0.78686 91311 45613 259D -3, 0.18463 18317 51005 468D -4, * 0.14215 11758 31644 589D -6, 0.20442 63103 38993 979D-14/ 101 format(3X,'quastn=',F12.7) c Q=F-HALF IF(DABS(Q).GT.SPLIT1)GOTO 10 R=CONST1-Q*Q QUA=Q*(((((((A7*R+A6)*R+A5)*R+A4)*R+A3)*R+A2)*R+A1)*R+A0) * /(((((((B7*R+B6)*R+B5)*R+B4)*R+B3)*R+B2)*R+B1)*R+ONE) RETURN 10 R=F IF(Q.GE.ZERO)R=ONE-F IF(R.LE.ZERO)GOTO 1000 R=DSQRT(-DLOG(R)) IF(R.GT.SPLIT2)GOTO 20 R=R-CONST2 QUA=(((((((C7*R+C6)*R+C5)*R+C4)*R+C3)*R+C2)*R+C1)*R+C0) * /(((((((D7*R+D6)*R+D5)*R+D4)*R+D3)*R+D2)*R+D1)*R+ONE) GOTO 30 20 R=R-SPLIT2 QUA=(((((((E7*R+E6)*R+E5)*R+E4)*R+E3)*R+E2)*R+E1)*R+E0) * /(((((((F7*R+F6)*R+F5)*R+F4)*R+F3)*R+F2)*R+F1)*R+ONE) 30 IF(Q.LT.ZERO)QUA=-QUA RETURN C 1000 WRITE(9,*) 'ERROR ROUTINE QUASTN : ARGUMENT OF FUNCTION INVALID' QUA=ZERO RETURN END c c --------------------------------------------------------------------------- SUBROUTINE SORT(temparr,X,N,maxyrs) C* FORTRAN CODE WRITTEN FOR INCLUSION IN IBM RESEARCH REPORT RC20525, C* FORTRAN ROUTINES FOR USE WITH THE METHOD OF L-MOMENTS, VERSION 3 C* J. R. M. HOSKING C* IBM RESEARCH DIVISION C* T. J. WATSON RESEARCH CENTER C* YORKTOWN HEIGHTS C* NEW YORK 10598, U.S. C* VERSION 3 AUGUST 1996 C C SORTS THE ARRAY X INTO ASCENDING ORDER C METHOD USED IS SHELL SORT WITH SEQUENCE OF INCREMENTS AS IN C D.F.KNUTH (1969) 'THE ART OF COMPUTER PROGRAMMING', VOL.3, P.95 C IMPLICIT DOUBLE PRECISION (A-H,O-Z) real*8 X(N),temparr(maxyrs+1) IF(N.LE.1)RETURN do i=1,n x(i)=temparr(i) enddo J=4 DO 10 I=1,100 J=3*J+1 IF(J.GE.N)GOTO 20 10 CONTINUE 20 CONTINUE M=(J/3) DO 60 MM=1,100 M=M/3 IF(M.EQ.0)RETURN DO 50 I=M+1,N TEST=X(I) J=I DO 30 JJ=1,100 J=J-M IF(J.LE.0)GOTO 40 IF(TEST.GE.X(J))GOTO 40 X(J+M)=X(J) 30 continue 40 CONTINUE X(J+M)=TEST 50 continue 60 CONTINUE END