c c Program to estimate flood frequency in North Carolina c REAL rho, ss, & sum, years, yhat INTEGER i, staid, iout, istep, j, jpeak, nest, ireg INCLUDE 'nc.cmn' REAL dist(303), da(303), cf(303), region(303),s(6),lat(303), & lng(303), peak(303,8), sd(303), rhoc(303,303), mcon(303,303), & id(303) REAL d, c, r, WLS INTEGER indx(303) c CHARACTER*8 ylabsv(7) CHARACTER*32 vout CHARACTER*3 method c CHARACTER*1 iansw CHARACTER*14 reg(4) EXTERNAL WLS DATA reg/' Blue Ridge ',' Piedmont ', & ' Coastal Plain',' Sand Hills '/ c PRINT *, ' This program computes estimates of T-year floods' PRINT *, ' for ungaged sites in North Carolina based on ' PRINT *, ' the Regional Regression Equation (RRE) method or' PRINT *, ' the Region Of Influence (ROI) method. (see the ' PRINT *, ' report " Estimating the magnitude and frequency ' PRINT *, ' of floods in rural basins of North Carolina" by ' PRINT *, ' B.F. Pope, G.D. Tasker, and J.C. Robbins USGS ' PRINT *, ' Water Resources Investigations Report 01-4207). ' PRINT *, ' ' PRINT *, '* No warranty, expressed or implied, is made by the' PRINT *, '* USGS as to the accuracy and functioning of the ' PRINT *, '* program and related program material.' PRINT *, ' ++++++++++++++++++++++++++++++++++++++++++++++++++' PRINT *, ' Revised 9-26-01 ' c PRINT *, ' ' PRINT *, ' Estimated T-year peaks for the ungaged site will' PRINT *, ' be printed on screen along with estimates of' PRINT *, ' accuracy. The estimates will be based on regional' PRINT *, ' regression equations based on data from 3 different' PRINT *, ' geographic regions (RRE method) or on' PRINT *, ' regression equations based on data for 30 gaged sites' PRINT *, ' with similar basin characteristics (ROI method).' PRINT *, ' Output will be to a user named output file.' PRINT *, ' ENTER name of output file ' READ (*,9005) vout OPEN (16,file=vout) 10 CONTINUE PRINT *, ' ENTER indentifer for ungaged site' READ (*,9005) Siteid PRINT *, ' ENTER hydrologic area in' PRINT *, ' which site is located. (1,2,or 3)' PRINT *, ' 1 = Blue Ridge/Piedmont' PRINT *, ' 2 = Coastal Plain' PRINT *, ' 3 = Sand Hills' read(*,*) jreg if(jreg.eq.1)ireg=1 if(jreg.eq.2)ireg=3 if(jreg.eq.3)ireg=4 PRINT *, ' ENTER watershed characteristics for site' PRINT *, ' Drainage area (sq. mi.) =' READ (*,*) Area c PRINT *, ' Use regional regression equations or region of' PRINT *, ' influence method, or both methods?' PRINT *, ' (ENTER RRE or ROI or BOTH)' READ (*,*) method IF (method.EQ.'RRE' .OR. method.EQ.'rre'.OR.method.EQ.'BOT' & .OR.method.eq.'bot') THEN CALL RRE(ireg) IF (method.EQ.'RRE' .OR. method.EQ.'rre')STOP END IF if(ireg.eq.4)then PRINT *, ' ROI method does not apply to Sand Hills' STOP end if c OPEN (8,file='cgrid.krg') OPEN (10,file='roi.in') OPEN (12,file='roi.rho') OPEN (13,file='roi.rec') AK(1) = 0.0 AK(2) = 0.84162 AK(3) = 1.28155 AK(4) = 1.75069 AK(5) = 2.05375 AK(6) = 2.32635 Ak(7) = 2.57500 AK(8) = 2.87816 Xlab(1) = 'log(DA) ' Xlab(2) = 'log(CF) ' c Xlab(3) = 'log(PR) ' c Xlab(4) = 'log(SH) ' c Xlab(5) = 'log(EL) ' Xlab(6) = 'constant' c c read in estimation data c nest=0 do 1 i=1,350 read(10,200,end=99)id(i),sd(i), lat(i), lng(i), + da(i), cf(i), region(i), + peak(i,1),peak(i,2),peak(i,3), + peak(i,4),peak(i,5),peak(i,6),peak(i,7), peak(i,8) c c Transform to logs c da(i)=alog10(da(i)) cf(i)=alog10(cf(i)) do 42 jj=1,8 peak(i,jj)=alog10(peak(i,jj)) 42 continue read(13,201)(mcon(i,j),j=1,i) 200 format(f8.0,13x,e13.5,13x,5e13.5,13x,8e13.5) 201 format(20f4.0) read(12,201)(rhoc(i,j),j=1,i) nest=nest+1 do 1 j=1,i mcon(j,i)=mcon(i,j) rhoc(j,i)=rhoc(i,j) 1 continue c c Initialize and read in ungaged site information c 99 continue s(1)=.9342 s(2)=.0090 s(3)=.5717 s(4)=.2537 s(5)=.0010 c Nsites = 30 icall=0 40 Continue regx=0.0 if(ireg.eq.3)regx=1.0 PRINT *, ' For the ROI method Climate Factor is ' PRINT *, ' needed. The program determines this from' PRINT *, ' the latitude and longitude of the site' PRINT *, 'Enter site latitude (dd mm ss)' READ (*,*) altd, altm, alts PRINT *, 'Enter site longitude (dd mm ss)' READ (*,*) alnd, alnm, alns CALL CFX(icall,altd,altm,alts,alnd,alnm,alns,Cf2,Cf25,Cf100) icall=1 Rewind(8) c WRITE (*,9025) Siteid, reg(ireg),Area, Cf25 d = ALOG10(Area) c = ALOG10(Cf25) r = regx c c Compute distances c DO 50 i = 1, nest dist(i) = SQRT(((d-da(i))/s(1))**2+((c-cf(i))/s(2)) & **2+((r-region(i))/s(5))**2) 50 CONTINUE c c Rank distances c CALL INDEXX(nest,dist,indx) c c Select closest nsites stations for regression c loop through dependent variables and do regression c ysav=0.0 c DO 170 jpeak = 1, 8 Ne = 2 Uv(1) = 1.0 Uv(2) = d c Uv(3) = c c Uv(4) = p c Uv(5) = sha c Uv(6) = ele sig=0.0 write(16,2034)siteid 2034 format(' Data for Region of Influence method',/, & ' SITE ID: ',a32,/,' ID REGION', &' LAT LNG ', &' LOG(DA) LOG(CF) LOG(P2) LOG(P5)', &' LOG(P10) LOG(P25) LOG(50) LOG(100)', &' LOG(P200) LOG(P500)') Do 71 i=1,Nsites X(i,1) = 1.0 Xt(1,i) = 1.0 X(i,2) = da(indx(i)) c X(i,3) = cf(indx(i)) Xt(2,i) = X(i,2) c Xt(3,i) = X(i,3) c Xt(4,i) = X(i,4) Stano(i) = id(indx(i)) sig=sig+sd(indx(i)) Write (16,2033)id(indx(i)),region(indx(i)),lat(indx(i)), & lng(indx(i)), da(indx(i)), & cf(indx(i)),(peak(indx(i),jj),jj=1,8) 2033 format(f10.0,f5.0,2f7.2,10f10.5) 71 CONTINUE sig=sig/nsites DO 170 jpeak = 1, 8 DO 60 j = 1, 2 Ylab(j) = Xlab(j) 60 CONTINUE DO 70 i = 1, Nsites Y(i,1) = peak(indx(i),jpeak) Stano(i) = id(indx(i)) 70 CONTINUE sum = 0.0 ss = 0.0 c c compute regional average standard deviation c and time sampling error, sta(i), for each site c DO 80 i = 1, Nsites c sig = sig + sd(indx(i)) sum = sum + Y(i,1) ss = ss + Y(i,1)**2 80 CONTINUE c sig = sig/Nsites Yvar = (ss-sum**2/Nsites)/(Nsites-1.) Atse = 0.0 Atscov = 0.0 DO 100 i = 1, Nsites DO 90 j = 1, i years = mcon(indx(i),indx(j)) & /(mcon(indx(i),indx(i))*mcon(indx(j),indx(j))) rho = rhoc(indx(i),indx(j)) Cov(i,j) = rho*sig**2*(1.+rho*.5*AK(jpeak)**2)*years Cov(j,i) = Cov(i,j) IF (i.EQ.j) THEN Sta(i) = Cov(i,i) Atse = Atse + Sta(i) ELSE Atscov = Atscov + Cov(i,j) ENDIF 90 CONTINUE 100 CONTINUE c c do regression c istep = 0 iflag=0 110 CALL SECANT(WLS) c c check to see if predicted value is greater than previous prediction c if(yhat.lt.ysav)then print *, ' CAUTION: Predicted T-year flow is smaller' print *, ' than T-year flow with lower recurrence ' print *, ' interval. See output. ' end if c c output final model c iout = 99 CALL OUTPUT(iout,jpeak,yhat,staid,iflag) ysav=yhat 170 CONTINUE STOP 9005 FORMAT (a32) 9010 FORMAT (f10.0,13F10.5) 9015 FORMAT (20F4.0) 9016 FORMAT (' Estimates adjusted for proximity to station',i8) 9017 FORMAT (' Station',i8,' unknown: NO ADJUSTMENT MADE') 9018 FORMAT (' Difference in drainage area for Station',i8, $' too great: NO ADUSTMENT MADE') 9020 FORMAT (//,' Region of Influence Method',/, & ' Flood frequency estimates for',/,1x,a32,/,' area=', & f10.2,' : slope=',f7.2,' precip',f6.1,' elev',f5.0, & ' shape',f7.2,/, & ' RI DISCHARGE SE(%) EQ. YRS. 90%', & ' PRED. INTERVAL',/,t12,' (cms)') 9025 FORMAT (//,' REGION OF INFLUENCE METHOD',/, & ' Flood frequency estimates for',/,1x,a32,/, & ' Region:',a14,/, & ' Drainage area:',f7.2,' square miles', & /,' Climate Factor:',f7.2,//, & ' RI DISCHARGE - SE (%) + SE (%) 90%', & ' PRED. INTERVAL',/,t12,' (cfs)') END block data pname include 'nc.cmn' DATA pklab/' 2 ',' 5 ',' 10 ',' 25 ', &' 50 ',' 100 ',' 200 ',' 500 '/ end c Subroutine INDEXX indexs an array ARRIN of length N, outputs the c array INDX such that ARRIN(INDX(J)) is in ascending order for c J=1,2,..,N. The input quantities ARRIN and N are not changed c (ref. Numerical Recipes, p. 233) c SUBROUTINE INDEXX(N,ARRIN,INDX) REAL ARRIN(303), q INTEGER i, INDX(303), indxt, ir, j, l, N C DO 10 j = 1, N INDX(j) = j 10 CONTINUE l = N/2 + 1 ir = N 20 IF (l.GT.1) THEN l = l - 1 indxt = INDX(l) q = ARRIN(indxt) ELSE indxt = INDX(ir) q = ARRIN(indxt) INDX(ir) = INDX(1) ir = ir - 1 IF (ir.EQ.1) THEN INDX(1) = indxt RETURN ENDIF ENDIF i = l j = l + l 30 IF (j.LE.ir) THEN IF (j.LT.ir) THEN IF (ARRIN(INDX(j)).LT.ARRIN(INDX(j+1))) j = j + 1 ENDIF IF (q.LT.ARRIN(INDX(j))) THEN INDX(i) = INDX(j) i = j j = j + j ELSE j = ir + 1 ENDIF GOTO 30 ENDIF INDX(i) = indxt GOTO 20 END c c function computes regression coefficients and weighted SSE c REAL FUNCTION WLS(GAMA2) INCLUDE 'nc.cmn' REAL GAMA2, det, xtx(10,10), work(10,50), work2(10,50), & work1(1,50) INTEGER k c c c compute weighting matrix, wt c DO 10 k = 1, Nsites Cov(k,k) = GAMA2 + Sta(k) 10 CONTINUE CALL INVERT(Nsites,50,det,Wt,Cov) c c compute weighted xtx c CALL MLTPLY(work,Xt,Wt,Ne,Nsites,Nsites,10,10,50) CALL MLTPLY(xtx,work,X,Ne,Nsites,Ne,10,10,50) c c compute regression coefficients c CALL INVERT(Ne,10,det,Xtxinv,xtx) CALL MLTPLY(work,Xtxinv,Xt,Ne,Ne,Nsites,10,10,10) CALL MLTPLY(work2,work,Wt,Ne,Nsites,Nsites,10,10,50) CALL MLTPLY(Bhat,work2,Y,Ne,Nsites,1,10,10,50) c c compute sum of errors c CALL MLTPLY(E,X,Bhat,Nsites,Ne,1,50,50,10) DO 20 k = 1, Nsites E(k,1) = Y(k,1) - E(k,1) Et(1,k) = E(k,1) 20 CONTINUE CALL MLTPLY(work1,Et,Wt,1,Nsites,Nsites,1,1,50) CALL MLTPLY(C1,work1,E,1,Nsites,1,1,1,50) c WLS = (Nsites-Ne)/C1(1,1) - 1. RETURN END c c SUBROUTINE SECANT(FX) REAL f1, f2, f3, fnew, FX, x1, x2, x3, xnew INTEGER i, j INCLUDE 'nc.cmn' EXTERNAL FX x1 = 0.0 x3 = 0.0 x2 = Yvar*2. f2 = FX(x2) f1 = FX(x1) IF (f1.LT.0.) THEN c c midpoint serch for good starting point c DO 10 j = 1, 3 xnew = (x1+x2)/2. fnew = FX(xnew) IF (fnew.LT.0.) THEN x1 = xnew fnew = fnew ELSE x2 = xnew f2 = fnew ENDIF 10 CONTINUE c c search for gama sq using secant serch c DO 20 i = 1, 30 x3 = x1 - f1*(x2-x1)/(f2-f1) IF (x3.LT.0.) THEN x3 = AMIN1(x2,x1)/2. f3 = FX(x3) ELSE f3 = FX(x3) ENDIF IF (ABS(f3).LT..0001) GOTO 30 IF (ABS(f1).LT.ABS(f2)) THEN x2 = x3 f2 = f3 ELSE x1 = x3 f1 = f3 ENDIF 20 CONTINUE ENDIF 30 Gamasq = x3 RETURN END C================================================== STUTP ======= FUNCTION STUTP(X,N) REAL a, b, GAUSCF, rhpi, STUTP, t, X, y, z INTEGER j, N, nn C C STUDENT T PROBABILITY C STUTP = PROB( STUDENT T WITH N DEG FR .LT. X ) C C NOTE - PROB(ABS(T).GT.X) = 2.*STUTP(-X,N) (FOR X .GT. 0.) C C SUBPGM USED - GAUSCF C C REF - G.W. HILL, ACM ALGOR 395, OCTOBER 1970. C C USGS - WK 12/79. C C DATA rhpi/.63661977/ C STUTP = .5 IF (N.LT.1) RETURN C nn = N z = 1. t = X**2 y = t/nn b = 1.0 + y C IF (.NOT.(nn.GE.20.AND.t.LT.nn.OR.nn.GT.200)) THEN C ( OR IF NN NON-INTEGER) C IF (nn.LT.20 .AND. t.LT.4.) THEN C C -- NESTED SUMMATION OF COSINE SERIES y = SQRT(y) a = y IF (nn.EQ.1) a = 0. ELSE C C -- TAIL SERIES FOR LARGE T a = SQRT(b) y = a*nn j = 0 10 j = j + 2 IF (a.EQ.z) THEN nn = nn + 2 z = 0. y = 0. a = -a ELSE z = a y = y*(j-1)/(b*j) a = a + y/(nn+j) GOTO 10 ENDIF ENDIF 20 nn = nn - 2 IF (nn.LE.1) THEN IF (nn.EQ.0) a = a/SQRT(b) IF (nn.NE.0) a = (ATAN(y)+a/b)*rhpi STUTP = 0.5*(z-a) IF (X.GT.0.) STUTP = 1. - STUTP RETURN ELSE a = (nn-1)/(b*nn)*a + y GOTO 20 ENDIF ENDIF C C -- ASYMPTOTIC SERIES FOR LARGE OR NONINTEGER N IF (y.GT.1E-6) y = ALOG(b) a = nn - 0.5 b = 48.*a**2 y = a*y y = (((((-0.4*y-3.3)*y-24.)*y-85.5)/(0.8*y**2+100.+b)+y+3.)/b+1.) & *SQRT(y) STUTP = GAUSCF(-y) IF (X.GT.0.) STUTP = 1. - STUTP RETURN C END SUBROUTINE INVERT(N,NDIM,DET,COVINV,COV) c IMPLICIT REAL*8 (A-H,O-Z) INTEGER i, im, j, k, N, NDIM REAL DET, detl, sum, temp, COVINV(NDIM,NDIM), COV(NDIM,NDIM), & b(50,50), a(50,50) C -------------------------------------------------------------- C COV IS AN N*N MATRIX C SUBROUTINE COMPUTES DETERMINANT OF COV THEN REPLACES COV WITH ITS C INVERSE C B IS THE LOWER TRIANGULAR DECOMPOSITION OF COV C -------------------------------------------------------------- c COMMON /PSTAR/ PINV(75,75) C IF (N.EQ.2) THEN DET = COV(1,1)*COV(2,2) - COV(1,2)**2 temp = COV(1,1)/DET COVINV(1,1) = COV(2,2)/DET COVINV(2,2) = temp COVINV(1,2) = -COV(1,2)/DET COVINV(2,1) = COVINV(1,2) ELSE C c WRITE(*,9)((B(I,J),J=1,N), I=1,N) CALL DECOMP(N,NDIM,COV,b) detl = b(1,1) DO 10 i = 2, N detl = detl*b(i,i) 10 CONTINUE DET = detl**2 c IF ( DET.EQ.0.0 ) THEN c WRITE (6,9010) C STOP c ENDIF C a(1,1) = 1./b(1,1) a(2,2) = 1./b(2,2) a(2,1) = -b(2,1)*a(1,1)*a(2,2) C DO 40 i = 3, N a(i,i) = 1./b(i,i) im = i - 1 DO 30 k = 1, im sum = 0. DO 20 j = k, im sum = sum + b(i,j)*a(j,k) 20 CONTINUE a(i,k) = -sum*a(i,i) 30 CONTINUE 40 CONTINUE C DO 70 i = 1, N DO 60 j = 1, i sum = 0. DO 50 k = i, N sum = sum + a(k,i)*a(k,j) 50 CONTINUE COVINV(i,j) = sum COVINV(j,i) = sum 60 CONTINUE 70 CONTINUE ENDIF RETURN 9005 FORMAT (1X,' PROCESSING STOPPED -- SINGULAR MATRIX') END SUBROUTINE DECOMP(N,NDIM,XLAM,B) c IMPLICIT REAL*8 (A-H,O-Z) INTEGER is, ism, js, jsm, ks, N, NDIM REAL XLAM(NDIM,NDIM), B(50,50), bh, bn C -------------------------------------------------------------- C CHOLESKY DECOMPOSITION BB-TRANSPOSE = XLAM C -------------------------------------------------------------- IF (XLAM(1,1).LE.0. .OR. XLAM(2,2).LE.0.) WRITE (*,9005) N, NDIM, & XLAM(1,1), XLAM(2,1), XLAM(2,2), XLAM(1,2) B(1,1) = SQRT(XLAM(1,1)) B(1,2) = 0. B(2,1) = XLAM(2,1)/B(1,1) C WRITE(1,9010)B(2,2) B(2,2) = SQRT(XLAM(2,2)-B(2,1)**2) CC WRITE(6,9015) NDIM, B(1,1),B(2,1),B(2,2) IF (N.LE.2) RETURN C DO 30 is = 3, N B(is,1) = XLAM(is,1)/B(1,1) bn = XLAM(is,is) - B(is,1)**2 ism = is - 1 DO 20 js = 2, ism jsm = js - 1 bh = XLAM(is,js) DO 10 ks = 1, jsm bh = bh - B(is,ks)*B(js,ks) 10 CONTINUE B(is,js) = bh/B(js,js) bn = bn - B(is,js)**2 20 CONTINUE IF (bn.LE.0.) WRITE (1,9010) bn B(is,is) = SQRT(AMAX1(bn,0.00)) 30 CONTINUE RETURN 9005 FORMAT (' IN DECOMP/ N, NDIM,XLAM 1-1,2-1,2-2,1-2 = ',2I5,/, & 4G13.4,/,' COVARIANCE MATRIX NOT POSITIVE DEFINITE') 9010 FORMAT (1X,' COVARIANCE MATRIX NOT POSITIVE DEFINITE BN=',F10.3) END C###gausex.spg processed by SPAG 3.023 at 12:04 on 8 Feb 1995 C================================================== C================================================== GAUSEX ====== FUNCTION GAUSEX(EXPROB) REAL ax, c0, c1, c2, CUMPRB, d, d1, d2, d3, EXPROB, GAUSEX, p, pr, & t, xlim, XX C C GAUSSIAN PROBABILITY FUNCTIONS W.KIRBY JUNE 71 C GAUSEX=VALUE EXCEEDED WITH PROB EXPROB C GAUSAB=VALUE (NOT EXCEEDED) WITH PROBCUMPROB C GAUSCF=CUMULATIVE PROBABILITY FUNCTION C GAUSDY=DENSITY FUNCTION C SUBPGMS USED -- NONE C C GAUSCF MODIFIED 740906 WK -- REPLACED ERF FCN REF BY RATIONAL APPRX N C ALSO REMOVED DOUBLE PRECISION FROM GAUSEX AND GAUSAB. C 76-05-04 WK -- TRAP UNDERFLOWS IN EXP IN GUASCF AND DY. C C DATA xlim/18.3/ DATA c0, c1, c2/2.51551700, .8028530000, .0103280000/ DATA d1, d2, d3/1.432788000, .1892690000, .0013080000/ C p = EXPROB 10 IF (p.GE.1.0) THEN GAUSEX = -10. ELSEIF (p.GT.0.) THEN pr = p IF (p.GT..5) pr = 1.00 - pr t = SQRT(-2.00*ALOG(pr)) GAUSEX = t - (c0+t*(c1+t*c2))/(1.0+t*(d1+t*(d2+t*d3))) IF (p.GT..5) GAUSEX = -GAUSEX ELSE GAUSEX = +10. ENDIF RETURN C ENTRY GAUSAB(CUMPRB) GAUSAB = 0. p = 1. - CUMPRB GOTO 10 C ENTRY GAUSCF(XX) ax = ABS(XX) GAUSCF = 1. IF (ax.LE.xlim) THEN t = 1.0/(1.0+.2316419*ax) d = 0.3989423*EXP(-XX*XX*.5) GAUSCF = 1. - & d*t*((((1.330274*t-1.821256)*t+1.781478)*t-0.3565638) & *t+0.3193815) ENDIF IF (XX.LT.0) GAUSCF = 1. - GAUSCF RETURN C ENTRY GAUSDY(XX) GAUSDY = 0. IF (ABS(XX).GT.xlim) RETURN GAUSDY = .3989423*EXP(-.500*XX*XX) RETURN END c c SUBROUTINE MLTPLY(PROD,X,Y,K1,K2,K3,N1,N2,N3) c IMPLICIT REAL*8 (A-H,O-Z) INTEGER i, j, k, K1, K2, K3, N1, N2, N3 REAL PROD(N1,K3), X(N2,K2), Y(N3,K3), sum C -------------------------------------------------------------- C X IS K1*K2 MATRIX C Y IS K2*K3 MATRIX C PROD = X*Y IS A K1*K3 MATRIX C -------------------------------------------------------------- DO 30 i = 1, K1 DO 20 k = 1, K3 sum = 0. DO 10 j = 1, K2 sum = sum + X(i,j)*Y(j,k) 10 CONTINUE PROD(i,k) = sum 20 CONTINUE 30 CONTINUE RETURN END c subroutine outputs results to screen and file c SUBROUTINE OUTPUT(IOUT,IPK,PRU,STAID,IFLAG) REAL AK, arhoc, sig REAL eqyrs, cl90, cookd, cu90, delres, errmod, hatdig, hmax, pred, & press, pru, prx, pv4, resid, samerr, sdbeta, sepc, sepu, & smod, hat(50,50), hats(50,50) REAL ssam, stdres, STUTP, tbeta, test1, test2, tstat, tv4, & varres, vpu, work1(50,1), work3(10,50), work2(10,1), xo(1,10) & , xot(10,1), ccc(1,1) INTEGER i, IOUT, IPK, iu, l, ndf,staid INCLUDE 'nc.cmn' C C pru = 0.0 DO 30 iu = 1, Ne pru = pru + Bhat(iu,1)*Uv(iu) 30 CONTINUE IF (IOUT.GT.0) THEN WRITE (16,9005) Siteid, Pklab(IPK) c WRITE (*,9005) Pklab(IPK) C C WRITE BETAS C c WRITE (*,9010) IOUT IF (IOUT.LT.99) WRITE (16,9010) IOUT WRITE (16,9015) c WRITE (*,9015) sdbeta = SQRT(Xtxinv(1,1)) tbeta = Bhat(1,1)/sdbeta c WRITE (*,9020) Bhat(1,1), sdbeta, tbeta WRITE (16,9020) Bhat(1,1), sdbeta, tbeta DO 10 i = 2, Ne sdbeta = SQRT(Xtxinv(i,i)) tbeta = Bhat(i,1)/sdbeta ndf = Nsites - Ne tv4 = ABS(tbeta) pv4 = 2.0*STUTP(-tv4,ndf) IF (pv4.LT..0001) pv4 = .0001 WRITE (16,9025) Ylab(i-1), Bhat(i,1), sdbeta, tbeta, pv4 c WRITE (*,9025) Ylab(i-1), Bhat(i,1), sdbeta, tbeta, pv4 10 CONTINUE CALL MLTPLY(work1,X,Bhat,Nsites,Ne,1,50,50,10) C C WRITE PREDICTED VALUES ETC. C IF (IOUT.EQ.99) THEN C WRITE (16,9030) c WRITE (*,9030) smod = 0. ssam = 0. press = 0.0 hmax = 0.0 CALL MLTPLY(work3,Xtxinv,Xt,Ne,Ne,Nsites,10,10,10) CALL MLTPLY(hats,X,work3,Nsites,Ne,Nsites,50,50,10) CALL MLTPLY(hat,hats,Wt,Nsites,Nsites,Nsites,50,50,50) DO 20 i = 1, Nsites pred = work1(i,1) resid = Y(i,1) - work1(i,1) samerr = hats(i,i) IF (samerr.GT.hmax) hmax = samerr hatdig = hat(i,i) delres = resid/(1.0-hatdig) errmod = Gamasq varres = Gamasq + Sta(i) - samerr stdres = resid/SQRT(varres) cookd = stdres**2*samerr/(Ne*(Gamasq+Sta(i)-samerr)) test1 = 2.*Ne/Nsites test2 = 4./Nsites c IF ( hatdig.GT.test1 .OR. cookd.GT.test2 ) THEN C WRITE (16,9035) Stano(i), Y(i,1), pred, stdres, hatdig, C & cookd c WRITE (*,9035) Stano(i), Y(i,1), pred, stdres, hatdig, c & cookd c ENDIF ssam = ssam + samerr press = press + delres**2 20 CONTINUE C C WRITE AVG SAMPLING ERROR AND MODEL ERROR C ssam = ssam/Nsites smod = smod/Nsites press = press/Nsites Atse = Atse/Nsites Atscov = Atscov/(Nsites/2.0*(Nsites-1)) arhoc = Atscov/Atse C WRITE (16,9040) ssam, errmod, press, Atse, arhoc, hmax c WRITE (*,9040) ssam, errmod, press, Atse,Arhoc, hmax c c Write out prediction for ungaged site c DO 40 l = 1, Ne c sum = 0.0 c DO 90 j = 1, Ne c sum = sum + Xtxinv(l,j)*Uv(j) c90 CONTINUE c work2(l,1) = sum xo(1,l) = Uv(l) xot(l,1) = Uv(l) 40 CONTINUE c sum = 0.0 c DO 120 j = 1, Ne c sum = sum + Uv(j)*work2(j,1) CALL MLTPLY(work2,Xtxinv,xot,Ne,Ne,1,10,10,10) CALL MLTPLY(ccc,xo,work2,1,Ne,1,1,1,10) sepu = ccc(1,1) vpu = Gamasq + sepu eqyrs = sig**2*(1.+AK(IPK)**2/2.)/vpu c c make adjustments c prx = 10**pru c IF (iflag.eq.3)THEN c CALL GADJ (area,staid,eqyrs,ipk,radj,yadj,iflag) c prx=prx*radj c vpu=vpu*eqyrs/(eqyrs+yadj) c eqyrs=eqyrs+yadj c END IF c c convert to cfs and percent error c sepc = 100.*SQRT(EXP(vpu*5.302)-1.) asep=sqrt(vpu) splus=100.*(10**(asep)-1.0) sminu=100.*(10**(-asep)-1.0) tstat = SQRT(vpu)*1.65 cu90 = 10**(tstat+pru) cl90 = 10**(pru-tstat) call round(prx) call round(cl90) call round(cu90) c IF (units.EQ.'M' .OR. units.EQ.'m') THEN c WRITE (16,9050) Siteid, Area, Slope, Precip, Elev, Shape, c & Pklab(IPK), prx, sepc, eqyrs, cl90, cu90 c WRITE (*,9060) Pklab(IPK), prx, sepc, eqyrs, cl90, cu90 c ELSE WRITE (16,9045) Siteid, Area, CF25, & Pklab(IPK), prx, sminu,splus, cl90, cu90 WRITE (*,9055) Pklab(IPK), prx, sminu,splus, cl90, cu90 c ENDIF IF (sepu.GT.hmax) WRITE (16,9065) IF ( sepu.GT.hmax ) then print *, ' WARNING: Prediction is an extrapolation beyond' print *, ' observed data. Check for errors in input basin' print *, ' characteristics. If no errors use results with' print *, ' caution.' END IF c ENDFILE (16) ENDIF ENDIF RETURN 9005 FORMAT (//,' SUMMARY OF REGION OF INFLUENCE REGRESSION FOR',/, & 1x,a32,/,1x,a8,'YR-PEAK') 9010 FORMAT (/,' **** STEP',i2,'****') 9015 FORMAT (/,' REGRESSION COEFFICIENTS',/, & ' VARIABLE COEFFICIENT STANDARD ERROR T FOR H0:BETA=0' & ,' PROB>|T|',/) 9020 FORMAT (' CONSTANT',5X,3F15.5) 9025 FORMAT (2X,A8,5X,3F15.5,f15.4) 9030 FORMAT (//,' Residuals and influence statistics ',/, & ' id ',t15,' obs',t27,' pred',t39,' std res',t51, & ' leverage',t62,' cook D',/) 9035 FORMAT (1X,F10.0,10F12.5) 9040 FORMAT (//,' AVERAGE SAMPLING ERROR VARIANCE',F10.4,/, & ' AVERAGE MODEL ERROR VARIANCE ',F10.4,/, & ' PRESS/N ',F10.4,/, & ' AVERAGE SQ. TIME-SAMPLING ERROR',f10.4,/, & ' Average Cross Correlation ',f10.4,/, & ' H MAX ',f10.4) 9045 FORMAT (//,t10,' *********************************************',/, & ' For ',a30,/,' area=',f10.2,' : cf=',f7.2, & /,a8,'YEAR PEAK',/, & ' PREDICTED(cfs) - SE (%) + SE (%) 90% PRED INT' & ,/,f12.0,f12.0,f12.0,f12.0,f12.0,/,t10, & ' *********************************************') 9050 FORMAT (//,t10,' *********************************************',/, & ' For ',a30,/,' area=',f10.2,' : Cf=',f7.2, & /,a8,'YEAR PEAK',/, & ' PREDICTED(cms) - SE (%) + SE (%) 90% PRED INT' & ,/,f12.1,f12.0,f12.2,f12.1,f12.1,/,t10, & ' *********************************************') 9055 FORMAT (a8,f12.0,f12.0,f12.0,f12.0,f12.0) 9060 FORMAT (a8,f12.1,f12.0,f12.2,f12.1,f12.1) 9065 FORMAT (/, & ' WARNING: Prediction is outside range of observed data ' & ,/) END c c SUBROUTINE ROUND(X) REAL div, test, X INTEGER i, ix c c Subroutine rounds peaks to three significant figures c for writing in table. c DO 10 i = 3, 8 test = 10**i div = 10**(i-2) IF (X.GT.test) THEN X = X/div ix = X + .5 X = div*ix ENDIF 10 CONTINUE RETURN END c c SUBROUTINE CFX(ICALL,ALTD,ALTM,ALTS,ALND,ALNM,ALNS,C2,C25,C100) REAL ALND, ALNM, ALNS, ALTD, ALTM, ALTS, C100, C100a, C2, C25, & C25a, C2a, sx, sy INTEGER i, ix, iy, j, npts c subroutine computes climate factors from lat and long c C COMPUTES X,Y COORDINATES(KILOMETERS) BASED ON LAT AND LONG, C USING ALBERS EQUAL-AREA CONIC PROJECTION--WITH EARTH AS ELLIPSOID. C EQUATIONS FROM BULLETIN 1532, PAGES 96-99. C A = a, equatorial radius of ellipsoid C E2 = e*2, square of eccentricity of ellipsoid (0.006768658). C PHI1,PHI2 = standard parallels of latitude (29.5,45.5). C PHI0 = middle latitude, or latitude chosen as the origin of y-coordinate ( c 38.0). C LAM0 = central meridian of the map, or longitude chosen as the origin C of x-coordinate (96.0). C PHI = latitude of desired y-coordinate. C LAM = longitude of desired x-coordinate. C XOFF = offset to x-coordinate to make all positive for KRIGING. C YOFF = offset to y-coordinate to make all positive for KRIGING. COMMON /CT / C2a(30,26), C25a(30,26), C100a(30,26) DOUBLE PRECISION q, q0, q1, q2, phi, phi0, phi1, phi2, e, e2, m1, & m2, n, rho, rho0, a, c, theta, lam, lam0, x, y, & xoff, yoff, dp, mp, sp, dm, mm, sm, one, two, tpi DATA a/6378206.4/, e2/0.006768658/, tpi/6.2831853/ DATA phi1/29.5/, phi2/45.5/, phi0/38.0/, lam0/96.0/ c OPEN (8,file='cgrid.krg') xoff = 1000.0D0 yoff = 1500.0D0 one = 1.0D0 two = 2.0D0 e = DSQRT(e2) phi0 = phi0*tpi/360.D0 phi1 = phi1*tpi/360.D0 phi2 = phi2*tpi/360.D0 q0 = (one-e2)*(DSIN(phi0)/(one-e2*DSIN(phi0)*DSIN(phi0)) & -(one/(two*e))*DLOG((one-e*DSIN(phi0))/(one+e*DSIN(phi0)))) q1 = (one-e2)*(DSIN(phi1)/(one-e2*DSIN(phi1)*DSIN(phi1)) & -(one/(two*e))*DLOG((one-e*DSIN(phi1))/(one+e*DSIN(phi1)))) q2 = (one-e2)*(DSIN(phi2)/(one-e2*DSIN(phi2)*DSIN(phi2)) & -(one/(two*e))*DLOG((one-e*DSIN(phi2))/(one+e*DSIN(phi2)))) m1 = DCOS(phi1)/DSQRT(one-e2*DSIN(phi1)*DSIN(phi1)) m2 = DCOS(phi2)/DSQRT(one-e2*DSIN(phi2)*DSIN(phi2)) n = (m1*m1-m2*m2)/(q2-q1) c = m1*m1 + n*q1 rho0 = a*DSQRT(c-n*q0)/n C ZERO ARRAYS THEN READ IN AVAILABLE KRIGED CLIMATE FACTORS if(icall.eq.1)go to 201 DO 100 i = 1, 30 DO 50 j = 1, 26 C2a(i,j) = 0.0 C25a(i,j) = 0.0 C100a(i,j) = 0.0 50 CONTINUE 100 CONTINUE READ (8,9005) npts DO 200 i = 1, npts READ (8,9010) ix, iy, C2, C25, C100 C2a(ix,iy) = C2 C25a(ix,iy) = C25 C100a(ix,iy) = C100 200 CONTINUE 201 CONTINUE dp = DBLE(ALTD) mp = DBLE(ALTM) sp = DBLE(ALTS) dm = DBLE(ALND) mm = DBLE(ALNM) sm = DBLE(ALNS) mp = mp + sp/60.D0 phi = (dp+mp/60.D0)*tpi/360.D0 mm = mm + sm/60.D0 lam = (dm+mm/60.D0) q = (one-e2)*(DSIN(phi)/(one-e2*DSIN(phi)*DSIN(phi))-(one/(two*e)) & *DLOG((one-e*DSIN(phi))/(one+e*DSIN(phi)))) theta = n*(lam0-lam)*tpi/360.D0 rho = a*DSQRT(c-n*q)/n x = rho*DSIN(theta)/1000.D0 + xoff y = (rho0-rho*DCOS(theta))/1000.D0 + yoff sx = SNGL(x) sy = SNGL(y) CALL WGT(sx,sy,C2,C25,C100) RETURN 9005 FORMAT (I4) 9010 FORMAT (I8,I9,3F15.0) 9015 FORMAT (t9,2F11.3) 9020 FORMAT (A20) 9025 FORMAT (2X,A8,6X,3F2.0,F3.0,2F2.0) 9030 FORMAT (2F8.1,3F8.5) END SUBROUTINE WGT(X,Y,C2,C25,C100) REAL C100, C100a, C2, C25, C25a, C2a, cp, r1, r2, r3, r4, rx, ry, & wr1, wr2, wr3, wr4, X, xr, Y, yr INTEGER ix, ixp, iy, iyp COMMON /CT / C2a(30,26), C25a(30,26), C100a(30,26) wr1 = 0.0 wr2 = 0.0 wr3 = 0.0 wr4 = 0.0 C ESTIMATE CLIMATE FACTOR FROM INTERPOLATION OF KRIGED VALUES AT 4 POINTS C SURROUNDING X,Y LOCATION OF INTEREST. HOWEVER, USE KRIGED VALUE IF C DISTANCE TO NEAREST POINT IS 10KM OR LESS. ix = X/100.0 rx = ix xr = X - 100.0*rx iy = Y/100.0 ry = iy yr = Y - 100.0*ry ix = ix + 1 iy = iy + 1 r1 = SQRT(xr**2+yr**2) IF ( r1.LE.10. ) THEN wr1 = 1.0 GOTO 100 ENDIF r2 = SQRT(xr**2+(100.0-yr)**2) IF ( r2.LE.10. ) THEN wr2 = 1.0 GOTO 100 ENDIF r3 = SQRT((100.0-xr)**2+(100.0-yr)**2) IF ( r3.LE.10. ) THEN wr3 = 1.0 GOTO 100 ENDIF r4 = SQRT((100.0-xr)**2+yr**2) IF ( r4.LE.10. ) THEN wr4 = 1.0 GOTO 100 ENDIF cp = 1.0/(1.0/r1+1.0/r2+1.0/r3+1.0/r4) wr1 = cp/r1 wr2 = cp/r2 wr3 = cp/r3 wr4 = cp/r4 100 CONTINUE ixp = ix + 1 iyp = iy + 1 C2 = wr1*C2a(ix,iy) + wr2*C2a(ix,iyp) + wr3*C2a(ixp,iyp) & + wr4*C2a(ixp,iy) C25 = wr1*C25a(ix,iy) + wr2*C25a(ix,iyp) + wr3*C25a(ixp,iyp) & + wr4*C25a(ixp,iy) C100 = wr1*C100a(ix,iy) + wr2*C100a(ix,iyp) + wr3*C100a(ixp,iyp) & + wr4*C100a(ixp,iy) RETURN END c subroutine rre(ireg) include 'nc.cmn' integer ireg if (ireg.eq.1.or.ireg.eq.2)call regbrp if (ireg.eq.3)call regcp if (ireg.eq.4)call regsh return end c SUBROUTINE REGBRP INCLUDE 'nc.cmn' INTEGER i, ip, j REAL yhat, sepl, sepu, cl, cu,stut REAL bs, v, vt, xtxi, temp, temp2, vmodel, & t, vpi, xt1, xt2 DIMENSION bs(2,8), v(1,2), vt(2,1), xt1(2,2,4), xtxi(2,2), & temp(1,2), temp2(1,1), vmodel(8), & xt2(2,2,4) c c DATA stut/1.65/ c DATA bs/2.13053,0.70220,2.38355,0.67675, & 2.52355,0.66195,2.67772,0.64544, & 2.77948,0.63456,2.87223,0.62473, & 2.95800,0.61572,3.06300,0.60484/ c DATA vmodel/0.28865E-01,0.28798E-01, & 0.29724E-01,0.31790E-01, & 0.33842E-01,0.36231E-01, & 0.38904E-01,0.42812E-01/ c DATA xt1/ & 0.14072E-02,-0.49612E-03, & -0.49612E-03, 0.24350E-03, & 0.16431E-02,-0.55517E-03, & -0.55517E-03, 0.26322E-03, & 0.18985E-02,-0.62424E-03, & -0.62424E-03, 0.28912E-03, & 0.22833E-02,-0.73140E-03, & -0.73140E-03, 0.33094E-03/ DATA xt2/ & 0.25999E-02,-0.82124E-03, & -0.82124E-03, 0.36687E-03, & 0.29342E-02,-0.91725E-03, & -0.91725E-03, 0.40581E-03, & 0.32839E-02,-0.10186E-02, & -0.10186E-02, 0.44737E-03, & 0.37671E-02,-0.11600E-02, & -0.11600E-02, 0.50586E-03/ c v(1,1) = 1.0 v(1,2) = ALOG10(area) vt(1,1) = 1.0 vt(2,1) = v(1,2) write(*,9025)siteid, area write(16,9025)siteid, area DO 30 ip = 1, 8 yhat = bs(1,ip) + bs(2,ip)*v(1,2) yhat = 10**yhat c c Compute CI c DO 20 i = 1, 2 DO 10 j = 1, 2 IF (ip.LE.4) THEN xtxi(i,j) = xt1(i,j,ip) ELSE xtxi(i,j) = xt2(i,j,ip-4) ENDIF 10 CONTINUE 20 CONTINUE CALL MLTPLY(temp,v,xtxi,1,2,2,1,1,2) CALL MLTPLY(temp2,temp,vt,1,2,1,1,1,2) vpi = vmodel(ip) + temp2(1,1) asep=sqrt(vpi) sepl = 100.*(10**(-asep) -1.0) sepu = 100.*(10**(asep) -1.0) t = 10**(stut*asep) cu = yhat*t cl = yhat/t call round(yhat) call round(cl) call round(cu) write(*,9000)pklab(ip),yhat,sepl,sepu,cl,cu write(16,9000)pklab(ip),yhat,sepl,sepu,cl,cu 9000 format(2x,a8,5f12.1) 30 CONTINUE IF (area.lt..14.or.area.gt.8386.)THEN WRITE(*,9101) WRITE(16,9101) END IF 9101 FORMAT (' WARNING - Drainage area out of range & of observed data') 9025 FORMAT (//,' REGIONAL REGRESSION METHOD',/, & ' Flood frequency estimates for',/,1x,a32,/, & ' Region: Blue Ridge/Piedmont',/, & ' Drainage area:',f7.2,' square miles',/, & ' RI DISCHARGE - SE (%) + SE (%) 90%', & ' PRED. INTERVAL',/,t12,' (cfs)') RETURN END SUBROUTINE REGCP INCLUDE 'nc.cmn' INTEGER i, ip, j REAL yhat, sepl, sepu, cl, cu,stut REAL bs, v, vt, xtxi, temp, temp2, vmodel, & t, vpi, xt1, xt2 DIMENSION bs(2,8), v(1,2), vt(2,1), xt1(2,2,4), xtxi(2,2), & temp(1,2), temp2(1,1), vmodel(8), & xt2(2,2,4) c c DATA stut/1.65/ c DATA bs/1.81081,0.67301,2.11209,0.63489, & 2.27330,0.61504,2.44894,0.59335, & 2.56462,0.57899,2.67018,0.56584, & 2.76806,0.55364,2.88828,0.53865/ c DATA vmodel/0.23562E-01,0.20959E-01, & 0.21143E-01,0.22706E-01, & 0.24530E-01,0.26755E-01, & 0.29313E-01,0.33116E-01/ c DATA xt1/ & 0.31893E-02,-0.98276E-03, & -0.98276E-03, 0.43777E-03, & 0.37147E-02,-0.10840E-02, & -0.10840E-02, 0.45031E-03, & 0.43917E-02,-0.12495E-02, & -0.12495E-02, 0.50084E-03, & 0.54496E-02,-0.15204E-02, & -0.15204E-02, 0.59196E-03/ DATA xt2/ & 0.63333E-02,-0.17517E-02, & -0.17517E-02, 0.67303E-03, & 0.72726E-02,-0.20005E-02, & -0.20005E-02, 0.76200E-03, & 0.82596E-02,-0.22640E-02, & -0.22640E-02, 0.85763E-03, & 0.96272E-02,-0.26319E-02, & -0.26319E-02, 0.99272E-03/ c v(1,1) = 1.0 v(1,2) = ALOG10(area) vt(1,1) = 1.0 vt(2,1) = v(1,2) write(*,9025)siteid, area write(16,9025)siteid, area DO 30 ip = 1, 8 yhat = bs(1,ip) + bs(2,ip)*v(1,2) yhat = 10**yhat c c Compute CI c DO 20 i = 1, 2 DO 10 j = 1, 2 IF (ip.LE.4) THEN xtxi(i,j) = xt1(i,j,ip) ELSE xtxi(i,j) = xt2(i,j,ip-4) ENDIF 10 CONTINUE 20 CONTINUE CALL MLTPLY(temp,v,xtxi,1,2,2,1,1,2) CALL MLTPLY(temp2,temp,vt,1,2,1,1,1,2) vpi = vmodel(ip) + temp2(1,1) asep=sqrt(vpi) sepl = 100.*(10**(-asep) -1.0) sepu = 100.*(10**(asep) -1.0) t = 10**(stut*asep) cu = yhat*t cl = yhat/t call round(yhat) call round(cl) call round(cu) write(*,9000)pklab(ip),yhat,sepl,sepu,cl,cu write(16,9000)pklab(ip),yhat,sepl,sepu,cl,cu 9000 format(2x,a8,5f12.1) 30 CONTINUE IF (area.lt..35.or.area.gt.8671.)THEN WRITE(*,9101) WRITE(16,9101) END IF 9101 FORMAT (' WARNING - Drainage area out of range & of observed data') 9025 FORMAT (//,' REGIONAL REGRESSION METHOD',/, & ' Flood frequency estimates for',/,1x,a32,/, & ' Region: Coastal Plain',/, & ' Drainage area:',f7.2,' square miles',/, & ' RI DISCHARGE - SE (%) + SE (%) 90%', & ' PRED. INTERVAL',/,t12,' (cfs)') RETURN end SUBROUTINE REGSH INCLUDE 'nc.cmn' INTEGER i, ip, j REAL yhat, sepl, sepu, cl, cu,stut REAL bs, v, vt, xtxi, temp, temp2, vmodel, & t, vpi, xt1, xt2 DIMENSION bs(2,8), v(1,2), vt(2,1), xt1(2,2,4), xtxi(2,2), & temp(1,2), temp2(1,1), vmodel(8), & xt2(2,2,4) c c DATA stut/1.78/ c DATA bs/1.52489,0.71201,1.74461,0.70091, & 1.86260,0.69679,1.99183,0.69300, & 2.07760,0.69062,2.15643,0.68842, & 2.23006,0.68632,2.32124,0.68363/ c DATA vmodel/0.21935E-01,0.26212E-01, & 0.29412E-01,0.34075E-01, & 0.37981E-01,0.42210E-01, & 0.46753E-01,0.53218E-01/ c DATA xt1/ & 0.10200E-01,-0.44648E-02, & -0.44648E-02, 0.25105E-02, & 0.12971E-01,-0.56042E-02, & -0.56042E-02, 0.31270E-02, & 0.15456E-01,-0.65943E-02, & -0.65943E-02, 0.36522E-02, & 0.19118E-01,-0.80470E-02, & -0.80470E-02, 0.44209E-02/ DATA xt2/ & 0.22136E-01,-0.92451E-02, & -0.92451E-02, 0.50554E-02, & 0.25348E-01,-0.10523E-01, & -0.10523E-01, 0.57330E-02, & 0.28740E-01,-0.11875E-01, & -0.11875E-01, 0.64516E-02, & 0.33482E-01,-0.13772E-01, & -0.13772E-01, 0.74616E-02/ c v(1,1) = 1.0 v(1,2) = ALOG10(area) vt(1,1) = 1.0 vt(2,1) = v(1,2) write(*,9025)siteid, area write(16,9025)siteid, area DO 30 ip = 1, 8 yhat = bs(1,ip) + bs(2,ip)*v(1,2) yhat = 10**yhat c c Compute CI c DO 20 i = 1, 2 DO 10 j = 1, 2 IF (ip.LE.4) THEN xtxi(i,j) = xt1(i,j,ip) ELSE xtxi(i,j) = xt2(i,j,ip-4) ENDIF 10 CONTINUE 20 CONTINUE CALL MLTPLY(temp,v,xtxi,1,2,2,1,1,2) CALL MLTPLY(temp2,temp,vt,1,2,1,1,1,2) vpi = vmodel(ip) + temp2(1,1) asep=sqrt(vpi) sepl = 100.*(10**(-asep) -1.0) sepu = 100.*(10**(asep) -1.0) t = 10**(stut*asep) cu = yhat*t cl = yhat/t call round(yhat) call round(cl) call round(cu) write(*,9000)pklab(ip),yhat,sepl,sepu,cl,cu write(16,9000)pklab(ip),yhat,sepl,sepu,cl,cu 9000 format(2x,a8,5f12.1) 30 CONTINUE IF (area.lt.2.19.or.area.gt.1228.)THEN WRITE(*,9101) WRITE(16,9101) END IF 9101 FORMAT (' WARNING - Drainage area out of range & of observed data') 9025 FORMAT (//,' REGIONAL REGRESSION METHOD',/, & ' Flood frequency estimates for',/,1x,a32,/, & ' Region: Sand Hills',/, & ' Drainage area:',f7.2,' square miles',/, & ' RI DISCHARGE - SE (%) + SE (%) 90%', & ' PRED. INTERVAL',/,t12,' (cfs)') RETURN end