C Last change: GSL 30 Sep 2003 12:49 pm c 2-segment linear fit added to HA3 for both SRE & MRE c Program to estimate flood frequency in Tennessee c NALL=453, NHA1=211, NHA2=115, NHA3=65, NHA4=62 REAL rho, ss, & sum, years, yhat INTEGER i, iout, istep, j, jpeak, nest, ireg INTEGER ltd, ltm, lts, lnd, lnm, lns INCLUDE 'tnff.cmn' REAL dist(490), da(490), cf(490), region(490),s(3), & peak(490,7), sd(490), rhoc(490,490), mcon(490,490), & id(490), sl(490), pf(490), glat(490), glng(490), & map(490) REAL d, c, WLS INTEGER indx(490) c CHARACTER*8 ylabsv(7) CHARACTER*32 vout,details CHARACTER*3 method c CHARACTER*1 iansw c CHARACTER*14 reg(4) COMMON /LTLN/ltd, ltm, lts, lnd, lnm, lns EXTERNAL WLS c Standard deviations N=453/CDA, SLOPE, CF2^M DATA s / 0.9039, 0.5311, 0.01880 / c DATA reg/' HA 1 ',' HA 2 ', c & ' HA 3 ',' HA 4 '/ c OPEN (8,file='cgrid.krg') OPEN (10,file='v203inp.txt') OPEN (12,file='v203m2') OPEN (13,file='v203m1') c PRINT *, ' This program predicts T-year flood-peak discharges' PRINT *, ' for unregulated streams in Tennessee using:' PRINT *, ' ' PRINT *, ' Single-Variable Regional-Regression Equations (SRE),' PRINT *, ' Multivariable Regional-Regression Equations (MRE),' PRINT *, ' and the Region-of-Influence Regression Method (ROI).' PRINT *, ' ' PRINT *, ' See the report "Flood-Frequency Prediction Methods for' PRINT *, ' Unregulated Streams of Tennessee, 2000" by George S.' PRINT *, ' Law and Gary D. Tasker. (U.S. Geological Survey Water-' PRINT *, ' Resources Investigations Report 03-4176.)' PRINT *, ' ' PRINT *, ' -------------------------------------------------' PRINT *, ' * No warranty, expressed or implied, is made by *' PRINT *, ' * the USGS as to the accuracy and functioning *' PRINT *, ' * of the program and related program material. *' PRINT *, ' * TDOT VERSION 2.0.3--30SEP2003 *' PRINT *, ' -------------------------------------------------' PRINT *, ' ' PRINT *, ' ' PRINT *, ' Output will be to two user-named output files.' PRINT *, ' ENTER name of output file for flood-peak discharges: ' READ (*,9005) vout OPEN (16,file=vout) 10 CONTINUE DO 3007 i=1,4 DO 3006 j=1,4 jwarn(i,j)=0 3006 continue 3007 continue PRINT *, ' ENTER indentifer for ungaged site: ' READ (*,9006) Siteid print *, ' Is drainage area for site in one Hydrologic' print *, ' Area or two Hydrologic Areas?' print *, ' ENTER 1 or 2: ' read (*,*) isplit if (isplit.eq.1)then xsplit=1.0 zsplit=0.0 psplit=100.0 qsplit=0.0 ireg2=0 end if do 5001 jsplit=1,isplit if (isplit.eq.1)then PRINT *, ' ENTER hydrologic area in which' PRINT *, ' site is located. (1,2,3, or 4): ' read(*,*) ireg end if if (isplit.eq.2.and.jsplit.eq.1)then PRINT *, ' ENTER first hydrologic area in which' PRINT *, ' site is located. (1,2,3, or 4): ' read(*,*) ireg PRINT *, ' ENTER percent of total basin in first' PRINT *, ' hydrologic area: Percent = ' read(*,*) psplit qsplit=100.0 - psplit xsplit=psplit/100.0 zsplit=1.0-xsplit PRINT *, ' ENTER second hydrologic area in which' PRINT *, ' site is located. (1,2,3, or 4): ' read(*,*) ireg2 end if if(isplit.eq.1.or.jsplit.eq.1)then PRINT *, ' ENTER watershed characteristics for site.' PRINT *, ' ' PRINT *, ' Suggested ranges for contributing drainage area:' PRINT *, ' Hydrologic area 1: 0.20 to 9000 square miles' PRINT *, ' Hydrologic area 2: 0.47 to 2557 square miles' PRINT *, ' Hydrologic area 3: 0.17 to 2048 square miles' PRINT *, ' Hydrologic area 4: 0.76 to 2308 square miles' PRINT *, ' ' PRINT *, ' ENTER contributing drainage area, in sq. miles: ' READ (*,*) Area PRINT *, ' ENTER site latitude (dd mm ss): ' READ (*,*) altd, altm, alts PRINT *, ' ENTER site longitude (dd mm ss): ' READ (*,*) alnd, alnm, alns ltd=altd ltm=altm lts=alts lnd=alnd lnm=alnm lns=alns c PRINT *, ' Use simple regression equations, SRE,' c PRINT *, ' or multiple regression equations, MRE,' c PRINT *, ' or region of influence, ROI ' c PRINT *, ' or all three methods, ALL ' c PRINT *, ' (ENTER SRE, MRE, ROI, or ALL)' c READ (*,*) method method='ALL' end if IF (method.EQ.'SRE' .OR. method.EQ.'sre'.OR.method.EQ.'ALL' & .OR.method.eq.'all') THEN if(isplit.eq.2.and.jsplit.eq.2)then CALL SRE(ireg2,jsplit) CALL REGOUT(2,1,ireg,ireg2) print *, ' Press ENTER to continue' read (*,221)blank 221 format(a1) else if(isplit.eq.1)then CALL SRE(ireg,jsplit) CALL REGOUT(1,1,ireg,ireg2) print *, ' Press ENTER to continue' read (*,221)blank else CALL SRE(ireg,jsplit) end if end if IF (method.EQ.'SRE' .OR. method.EQ.'sre')then if(isplit.eq.2.and.jsplit.eq.1)then go to 5001 else STOP end if END IF ENDIF if(isplit.eq.1.or.jsplit.eq.1)then PRINT *, ' For the MRE and ROI methods channel slope must be' PRINT *, ' given. The program determines the climate factor' PRINT *, ' and physiographic-region factor from the latitude, ' PRINT *, ' longitude, and hydrologic area of the ungaged site.' PRINT *, ' ' PRINT *, ' Suggested ranges for channel slope:' PRINT *, ' Hydrologic area 1: 3.29 to 950 feet/mile' PRINT *, ' Hydrologic area 2: 1.90 to 343 feet/mile' PRINT *, ' Hydrologic area 3: 2.12 to 132 feet/mile' PRINT *, ' Hydrologic area 4: 0.89 to 63 feet/mile' PRINT *, ' ' PRINT *, ' ENTER channel slope, in feet/mile: ' Read (*,*) slp c CALL CFX(icall,altd,altm,alts,alnd,alnm,alns,Cf2,Cf25,Cf100) icall=1 Rewind(8) end if IF (method.EQ.'MRE' .OR. method.EQ.'mre'.or.method.eq.'ALL' & .or.method.eq.'all')then if(isplit.eq.1.and.ireg.eq.4)then go to 5001 end if if(isplit.eq.2.and.jsplit.eq.1)then CALL MRE(ireg,jsplit) go to 5001 else if(isplit.eq.1)then CALL MRE(ireg,jsplit) CALL REGOUT(1,2,ireg,ireg2) IF (method.EQ.'MRE' .OR. method.EQ.'mre')STOP end if if(isplit.eq.2.and.jsplit.eq.2)then CALL MRE(ireg2,jsplit) CALL REGOUT(2,2,ireg,ireg2) IF (method.EQ.'MRE' .OR. method.EQ.'mre')STOP end if end if END IF 5001 continue c c ROI method c c reset jwarn c DO 3008 i=1,4 DO 3009 j=1,4 jwarn(i,j)=0 3009 continue 3008 continue PRINT *, ' ENTER name of output file for data and regression' PRINT *, ' details for the ROI method: ' READ (*,9005)details OPEN (18,file=details) 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.87816 Xlab(1) = 'LOG(CDA)' Xlab(2) = 'LOG(CS) ' Xlab(3) = 'LOG(PF) ' Xlab(4) = 'LOG(CF)' c Xlab(5) = 'log(EL) ' Xlab(6) = 'constant' c c read in estimation data c nest=0 do 1 i=1,500 read(10,200,end=99)id(i),sd(i), da(i),sl(i), + pf(i),map(i),cf (i),glat(i),glng(i),region(i), + peak(i,1),peak(i,2),peak(i,3), + peak(i,4),peak(i,5),peak(i,6),peak(i,7) c c Transform to logs c c da(i)=alog10(da(i)) c cf(i)=alog10(cf(i)) c do 42 jj=1,8 c peak(i,jj)=alog10(peak(i,jj)) c 42 continue read(13,201)(mcon(i,j),j=1,i) 200 format(2x,f8.0,4f12.5,f12.5,3f12.5,8f12.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 c Nsites = 60 icall=0 40 Continue c Physiographic region factor equations N=453, NHA1=211, NHA2=115, NHA3=65, NHA4=62 d = ALOG10(Area) c = ALOG10(Cf2) slope=ALOG10(slp) c c Compute distances c DO 50 i = 1, nest dist(i) = SQRT(((d-da(i))/s(1))**2+((c-cf(i))/s(3)) & **2+((slope-sl(i))/s(2))**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 if(ireg.eq.1)p=-.2130 + .0626*d if(ireg.eq.2)p= .0168 + .0353*d if(ireg.eq.3)p= .2319 - .0242*d if(ireg.eq.4)p= .3044 - .1541*d if(isplit.eq.2)then if(ireg2.eq.1)p2sav=-.2130 + .0626*d if(ireg2.eq.2)p2sav= .0168 + .0353*d if(ireg2.eq.3)p2sav= .2319 - .0242*d if(ireg2.eq.4)p2sav= .3044 - .1541*d end if c ysav=0.0 DO 170 jpeak = 1, 7 Ne = 5 Uv(1) = 1.0 Uv(2) = d Uv(3) = slope Uv(4) = p Uv(5) = c Uv2(1)=1.0 Uv2(2)=d Uv2(3)=slope Uv2(4)=p if(isplit.eq.2)Uv2(4)=p2sav Uv2(5)=c c Uv(6) = ele sig=0.0 c*** c write(*,7001)(Uv(iii),iii=1,5) c write(*,7002)(Uv2(iii),iii=1,5) c7001 format(' Uv ',5f10.4) c7002 format(' Uv2',5f10.4) c*** c c if(jpeak.eq.1) write(18,2034)siteid 2034 format(' DATA FOR TDOT VERSION 2.0.3 REGION-OF-INFLUENCE (ROI)', &' METHOD FOR TENNESSEE:',/, & ' SITE ID: ',a65,/,' ID HA ', &'LATITUDE LONGITUDE MAP NO. LOG(CDA) LOG(CS) ', &' LOG(PF) LOG(CF)') Do 71 i=1,Nsites X(i,1) = 1.0 Xt(1,i) = 1.0 X(i,2) = da(indx(i)) X(i,3) = sl(indx(i)) X(i,4) = pf(indx(i)) X(i,5) = cf(indx(i)) Xt(2,i) = X(i,2) Xt(3,i) = X(i,3) Xt(4,i) = X(i,4) Xt(5,i) = x(i,5) Stano(i) = id(indx(i)) sig=sig+sd(indx(i)) if(jpeak.eq.1)Write (18,2033)id(indx(i)),region(indx(i)), & glat(INDX(i)),glng(indx(i)),map(indx(i)),da(indx(i)), & sl(indx(i)), pf(indx(i)), & cf(indx(i)) 2033 format(f10.0,f5.0,2f10.5,f10.0,4f10.5) 71 CONTINUE sig=sig/nsites c DO 170 jpeak = 1, 7 DO 60 j = 1, 4 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 c iflag=0 1000 continue 110 CALL SECANT(WLS) istep = istep + 1 c c compute t and do step-backward drop of all variables with t<1.67 c if(Ne.eq.2)go to 9999 tsave = 2.0 isave = 100 DO 1100 i = 2, Ne sdbeta = SQRT(Xtxinv(i,i)) tbeta = ABS(Bhat(i,1)/sdbeta) IF ( tbeta.LT.1.67.AND. tbeta.LT.tsave ) THEN tsave = tbeta isave = i ENDIF 1100 CONTINUE IF ( isave.EQ.Ne ) THEN Ne = Ne - 1 GOTO 1000 ENDIF IF ( isave.LT.Ne ) THEN Ne = Ne - 1 DO 1150 j = isave, Ne Ylab(j-1) = Ylab(j) 1150 CONTINUE DO 1200 i = 1, Nsites DO 1160 j = isave, Ne X(i,j) = X(i,j+1) Xt(j,i) = X(i,j) 1160 CONTINUE 1200 CONTINUE DO 1250 j = isave, Ne Uv(j) = Uv(j+1) Uv2(j)=Uv2(j+1) 1250 CONTINUE GOTO 1000 ENDIF c c c output final model c 9999 continue c c check if pf is in equation c IOUT=99 CALL OUTPUT(iout,jpeak,yhat,ireg,ireg2,isplit) c c check to see if predicted value is greater than previous prediction c c if(yhat.lt.ysav)then c print *, ' CAUTION: Predicted T-year flow is smaller' c print *, ' than T-year flow with lower recurrence ' c print *, ' interval. See output. ' c end if c ysav=yhat 170 CONTINUE if(isplit.eq.1)then call regout(1,3,ireg,ireg2) STOP end if if(isplit.eq.2)then call regout(2,3,ireg,ireg2) STOP end if c5001 CONTINUE 9005 FORMAT (a32) 9006 FORMAT (a65) 9010 FORMAT (f10.0,13F10.5) 9015 FORMAT (20F4.0) END block data pname include 'tnff.cmn' DATA pklab/' 2 ',' 5 ',' 10 ',' 25 ', &' 50 ',' 100 ',' 500 ',' 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(490), q INTEGER i, INDX(490), 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 'tnff.cmn' REAL GAMA2, det, xtx(10,10), work(10,90), work2(10,90), & work1(1,90) 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,90,det,Wt,Cov) c c compute weighted xtx c CALL MLTPLY(work,Xt,Wt,Ne,Nsites,Nsites,10,10,90) CALL MLTPLY(xtx,work,X,Ne,Nsites,Ne,10,10,90) 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,90) CALL MLTPLY(Bhat,work2,Y,Ne,Nsites,1,10,10,90) c c compute sum of errors c CALL MLTPLY(E,X,Bhat,Nsites,Ne,1,90,90,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,90) CALL MLTPLY(C1,work1,E,1,Nsites,1,1,1,90) 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 'tnff.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 search 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 search 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(90,90), a(90,90) 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(90,90), 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,ireg,ireg2,isplit) REAL cookd, delres, errmod, hatdig, hmax, pred, & press, pru, pv4, resid, samerr, sdbeta, sepu, & smod, hat(90,90), hats(90,90) REAL ssam, stdres, STUTP, tbeta, tv4, & varres, vpu, work1(90,1), work3(10,90), work2(10,1), xo(1,10) & , xot(10,1), ccc(1,1) INTEGER i, IOUT, IPK, iu, l, ndf, isplit INCLUDE 'tnff.cmn' C C pru = 0.0 pru2=0.0 DO 30 iu = 1, Ne pru = pru + Bhat(iu,1)*Uv(iu) pru2= pru2+Bhat(iu,1)*Uv2(iu) 30 CONTINUE y1(ipk,1)=pru y1(ipk,2)=pru2 IF (IOUT.GT.0) THEN WRITE (18,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 (18,9015) c WRITE (*,9015) sdbeta = SQRT(Xtxinv(1,1)) tbeta = Bhat(1,1)/sdbeta c WRITE (*,9020) Bhat(1,1), sdbeta, tbeta WRITE (18,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 (18,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,90,90,10) C C WRITE PREDICTED VALUES ETC. C IF (IOUT.EQ.99) THEN WRITE (18,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,90,90,10) CALL MLTPLY(hat,hats,Wt,Nsites,Nsites,Nsites,90,90,90) 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)) c test1 = 2.*Ne/Nsites c test2 = 4./Nsites c IF ( hatdig.GT.test1 .OR. cookd.GT.test2 ) THEN WRITE (18,9035) Stano(i), Y(i,1), pred, stdres, hatdig, & cookd c WRITE (*,9035) Stano(i), Y(i,1), pred, stdres, hatdig, c & cookd 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)) c arhoc = Atscov/Atse WRITE (18,9040) ssam, errmod, press, hmax c WRITE (*,9040) ssam, errmod, press, Atse,Arhoc, hmax END IF END IF c c Write out prediction for ungaged site c DO 40 l = 1, Ne xo(1,l) = Uv(l) xot(l,1) = Uv(l) 40 CONTINUE 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 v1(ipk,1)=vpu sepu1=sepu DO 42 l = 1, Ne xo(1,l) = Uv2(l) xot(l,1) = Uv2(l) 42 CONTINUE CALL MLTPLY(work2,Xtxinv,xot,Ne,Ne,1,10,10,10) CALL MLTPLY(ccc,xo,work2,1,Ne,1,1,1,10) eqyrs = sig**2*(1.+AK(IPK)**2/2.)/vpu write (18,9041)ireg,sepu,eqyrs sepu = ccc(1,1) c write (18,9041)sepu vpu = Gamasq + sepu v1(ipk,2)=vpu sepu2=sepu eqyrs = sig**2*(1.+AK(IPK)**2/2.)/vpu if(isplit.eq.2)write(18,9042)ireg2,sepu2,eqyrs IF ( sepu1.GT.hmax ) jwarn(ireg,4)=1 IF ( sepu2.GT.hmax ) jwarn(ireg2,4)=1 RETURN 9005 FORMAT (//,' SUMMARY OF REGION-OF-INFLUENCE (ROI) REGRESSION FOR:', & /,1x,a65,/,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 ',t16,'LOG(OBS)',t27,'LOG(PRED)',t39,' STD RES', & ,t53,'LEVERAGE',t64,' COOKS D',/) 9035 FORMAT (1X,F10.0,10F12.5) 9040 FORMAT (//,' AVERAGE SAMPLING ERROR VARIANCE ',F10.4,/, & ' MODEL ERROR VARIANCE ',F10.4,/, & ' PRESS/N ',F10.4,/, & ' MAXIMUM SAMPLING ERROR VARIANCE ',f10.4) 9041 FORMAT (' HA',i2,' SITE',/, & ' SAMPLING ERROR VARIANCE ',f10.4,/, & ' EQUIVALENT YEARS ',f10.2) 9042 FORMAT (' HA',i2,' SITE',/, & ' SAMPLING ERROR VARIANCE ',f10.4,/, & ' EQUIVALENT YEARS ',f10.2) 9045 FORMAT (//,t10,' *********************************************',/, & ' For ',a30,/,' cda=',f10.2,': cs =',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,/,' cda=',f10.2,' : Cf=',f7.2, & /,a8,'YEAR PEAK',/, & ' PREDICTED(cms) - SE (%) + SE (%) 90% PRED INT' & ,/,f12.1,f12.1,f12.1,f12.1,f12.1,/,t10, & ' *********************************************') 9055 FORMAT (2x,a8,f12.1,f12.1,f12.1,f12.1,f12.1) 9060 FORMAT (2x,a8,f12.1,f12.1,f12.1,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 EYEARS(sig,zp,skew,vpi,eqyrs) REAL sig, zp, skew, vpi, whkp, hard, eqyrs c c compute equivalent years c c Wilson-Hilferty approximation (Handbook of Hydrology, eqn 18.2.29) if(skew.gt.0.0)then whkp=(2./skew)*(1.+(skew*zp/6.0)-skew**2/36.)**3 & -2.0/skew else whkp=zp end if c Handbook of hydrology eqn 18.4.11 hard=1.0+skew*whkp+0.5*(1.+.75*skew**2)*whkp**2 eqyrs = sig**2*(hard)/vpi RETURN END c c subroutine sre(ireg,jsplit) include 'tnff.cmn' integer ireg if (ireg.eq.1)call sre1(jsplit) if (ireg.eq.2)call sre2(jsplit) if (ireg.eq.3.and.area.lt.30.2)call sre3a(jsplit) if (ireg.eq.3.and.area.ge.30.2)call sre3b(jsplit) if (ireg.eq.4)call sre4(jsplit) return end subroutine mre(ireg,jsplit) include 'tnff.cmn' integer ireg if (ireg.eq.1)call mre1(jsplit) if (ireg.eq.2)call mre2(jsplit) if (ireg.eq.3.and.area.lt.30.2)call mre3a(jsplit) if (ireg.eq.3.and.area.ge.30.2)call mre3b(jsplit) if (ireg.eq.4)call mre4(jsplit) return end c These are the hydrologic area models for simple and multiple c regression as determined by GLSNET for NALL=453, NHA1=211, c NHA2=115, HHA3ALL=65:NHA3LOW=44:NHA3HGH=45 (OVERLAP), and NHA4=62. c **GSLaw 28APR2003** SUBROUTINE SRE1(jsplit) INCLUDE 'tnff.cmn' INTEGER i, ip, j REAL yhat REAL bs, v, vt, xtxi, temp, temp2, vmodel, & vpi, xt1, xt2 DIMENSION bs(2,7), v(1,2), vt(2,1), xt1(2,2,4), xtxi(2,2), & temp(1,2), temp2(1,1), vmodel(7), & xt2(2,2,3) INTEGER ltd,ltm,lts,lnd,lnm,lns COMMON /LTLN/ltd, ltm, lts, lnd, lnm, lns c SIMPLE REGRESSION METHOD FOR HYDROLOGIC AREA 1 c DATA stut/1.65/ c DATA bs/2.07650,0.75542, 2.29505,0.73968,2.41101,0.73116, & 2.53442,0.72180,2.61364,0.71566,2.68453,0.71010, & 2.82740,0.69869/ c DATA vmodel/0.31377E-01,0.30349E-01,0.31267E-01,0.33949E-01, & 0.36830E-01,0.40320E-01,0.50443E-01/ c DATA xt1/0.12582E-02, -0.46690E-03, & -0.46690E-03, 0.23477E-03, & 0.13882E-02, -0.49753E-03, & -0.49753E-03, 0.24364E-03, & 0.15629E-02, -0.54914E-03, & -0.54914E-03, 0.26466E-03, & 0.18443E-02, -0.63743E-03, & -0.63743E-03, 0.30301E-03/ DATA xt2/0.20857E-02, -0.71579E-03, & -0.71579E-03, 0.33814E-03, & 0.23480E-02, -0.80261E-03, & -0.80261E-03, 0.37770E-03, & 0.30281E-02, -0.10330E-02, & -0.10330E-02, 0.48458E-03/ c v(1,1) = 1.0 v(1,2) = ALOG10(area) vt(1,1) = 1.0 vt(2,1) = v(1,2) DO 30 ip = 1, 7 yhat = bs(1,ip) + bs(2,ip)*v(1,2) 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) if(jsplit.eq.1)then y1(ip,1)=yhat v1(ip,1)=vpi else y1(ip,2)=yhat v1(ip,2)=vpi end if 30 CONTINUE jwarn(1,1)=0 IF (area.lt..20.or.area.gt.21400.)THEN jwarn(1,1)=1 END IF RETURN END c SUBROUTINE SRE2(jsplit) INCLUDE 'tnff.cmn' INTEGER i, ip, j REAL yhat REAL bs, v, vt, xtxi, temp, temp2, vmodel, & vpi, xt1, xt2 DIMENSION bs(2,7), v(1,2), vt(2,1), xt1(2,2,4), xtxi(2,2), & temp(1,2), temp2(1,1), vmodel(7), & xt2(2,2,3) INTEGER ltd,ltm,lts,lnd,lnm,lns COMMON /LTLN/ltd, ltm, lts, lnd, lnm, lns c c SIMPLE REGRESSION METHOD FOR HYDROLOGIC AREA 2^M c DATA stut/1.66/ c DATA bs/2.30928,0.72677,2.53099,0.71634,2.64273,0.71231, & 2.75811,0.70885,2.83070,0.70691,2.89480,0.70530, & 3.02130,0.70234/ c DATA vmodel/0.17720E-01,0.15767E-01,0.16611E-01,0.18933E-01, & 0.21270E-01,0.23988E-01,0.31507E-01/ c DATA xt1/0.15347E-02,-0.52817E-03, & -0.52817E-03, 0.27021E-03, & 0.15708E-02,-0.51955E-03, & -0.51955E-03, 0.25962E-03, & 0.17846E-02,-0.58243E-03, & -0.58243E-03, 0.28809E-03, & 0.21585E-02,-0.69984E-03, & -0.69984E-03, 0.34395E-03/ DATA xt2/0.24827E-02,-0.80465E-03, & -0.80465E-03, 0.39478E-03, & 0.28234E-02,-0.91962E-03, & -0.91962E-03, 0.45115E-03, & 0.37198E-02,-0.12179E-02, & -0.12179E-02, 0.59926E-03/ c v(1,1) = 1.0 v(1,2) = ALOG10(area) vt(1,1) = 1.0 vt(2,1) = v(1,2) DO 30 ip = 1, 7 yhat = bs(1,ip) + bs(2,ip)*v(1,2) c 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) if(jsplit.eq.1)then y1(ip,1)=yhat v1(ip,1)=vpi else y1(ip,2)=yhat v1(ip,2)=vpi end if 30 CONTINUE IF (area.lt..47.or.area.gt.2557.)THEN jwarn(2,1)=1 END IF RETURN END c SUBROUTINE SRE3a(jsplit) INCLUDE 'tnff.cmn' INTEGER i, ip, j REAL yhat REAL bs, v, vt, xtxi, temp, temp2, vmodel, & vpi, xt1, xt2 INTEGER ltd,ltm,lts,lnd,lnm,lns COMMON /LTLN/ltd, ltm, lts, lnd, lnm, lns DIMENSION bs(2,7) DIMENSION vt(2,1),v(1,2) DIMENSION xt1(2,2,4), xtxi(2,2) DIMENSION temp(1,2), temp2(1,1), vmodel(7) DIMENSION xt2(2,2,3) c REVISED--------REVISED--------REVISED^M c SIMPLE REGRESSION METHOD FOR HYDROLOGIC AREA 3^M c 2-SEGMENT FIT LOW END DATA (DA<1.48 LOGS) (stut/1.67/^M c REVISED 28APR2003 GSLAW DATA bs/2.44649,0.78902,2.65552,0.76920,2.75926,0.76062, & 2.86515,0.75271,2.93071,0.74826,2.98768,0.74470, & 3.09705,0.73870/ c DATA vmodel/0.19593E-01,0.19321E-01,0.19686E-01,0.20591E-01, & 0.21559E-01,0.22758E-01,0.26396E-01/ c DATA xt1/0.25104E-02,-0.15590E-02, & -0.15590E-02, 0.15052E-02, & 0.25852E-02,-0.16006E-02, & -0.16006E-02, 0.15412E-02, & 0.28200E-02,-0.17292E-02, & -0.17292E-02, 0.16492E-02, & 0.32065E-02,-0.19431E-02, & -0.19431E-02, 0.18320E-02/ DATA xt2/0.35326E-02,-0.21264E-02, & -0.21264E-02, 0.19915E-02, & 0.38808E-02,-0.23250E-02, & -0.23250E-02, 0.21669E-02, & 0.47630E-02,-0.28401E-02, & -0.28401E-02, 0.26316E-02/ c v(1,1) = 1.0 v(1,2) = ALOG10(area) vt(1,1) = 1.0 vt(2,1) = v(1,2) DO 30 ip = 1, 7 yhat = bs(1,ip) + bs(2,ip)*v(1,2) c 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,w,xtxi,1,2,2,1,1,2) CALL MLTPLY(temp2,temp,vt,1,2,1,1,1,2) vpi = vmodel(ip) + temp2(1,1) if(jsplit.eq.1)then y1(ip,1)=yhat v1(ip,1)=vpi else y1(ip,2)=yhat v1(ip,2)=vpi end if 30 CONTINUE IF (area.lt..170.or.area.gt.2048.)THEN jwarn(3,1)=1 END IF RETURN END c SUBROUTINE SRE3b(jsplit) INCLUDE 'tnff.cmn' INTEGER i, ip, j REAL yhat REAL bs, v, vt, xtxi, temp, temp2, vmodel, & vpi, xt1, xt2 INTEGER ltd,ltm,lts,lnd,lnm,lns COMMON /LTLN/ltd, ltm, lts, lnd, lnm, lns DIMENSION bs(2,7) DIMENSION vt(2,1),v(1,2) DIMENSION xt1(2,2,4), xtxi(2,2) DIMENSION temp(1,2), temp2(1,1), vmodel(7) DIMENSION xt2(2,2,3) c REVISED--------REVISED--------REVISED^M c SIMPLE REGRESSION METHOD FOR HYDROLOGIC AREA 3^M c 2-SEGMENT FIT HIGH END DATA (DA>=1.48 LOGS) stut/1.67/^M c REVISED 28APR03 GSLAW DATA bs/2.83162,0.52708,3.01727,0.52335,3.10863,0.52319, & 3.20036,0.52477,3.25634,0.52685,3.30453,0.52946, & 3.39615,0.53688/ c DATA vmodel/0.12658E-01,0.13152E-01,0.14636E-01,0.17420E-01, & 0.20033E-01,0.23034E-01,0.31387E-01/ c DATA xt1/0.39573E-02,-0.17034E-02, & -0.17034E-02, 0.86375E-03, & 0.42430E-02,-0.18198E-02, & -0.18198E-02, 0.92296E-03, & 0.48656E-02,-0.20793E-02, & -0.20793E-02, 0.10529E-02, & 0.59361E-02,-0.25317E-02, & -0.25317E-02, 0.12808E-02/ DATA xt2/0.68823E-02,-0.29357E-02, & -0.29357E-02, 0.14856E-02, & 0.79270E-02,-0.33853E-02, & -0.33853E-02, 0.17144E-02, & 0.10687E-01,-0.45859E-02, & -0.45859E-02, 0.23290E-02/ c v(1,1) = 1.0 v(1,2) = ALOG10(area) vt(1,1) = 1.0 vt(2,1) = v(1,2) DO 30 ip = 1, 7 yhat = bs(1,ip) + bs(2,ip)*v(1,2) c 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,w,xtxi,1,2,2,1,1,2) CALL MLTPLY(temp2,temp,vt,1,2,1,1,1,2) vpi = vmodel(ip) + temp2(1,1) if(jsplit.eq.1)then y1(ip,1)=yhat v1(ip,1)=vpi else y1(ip,2)=yhat v1(ip,2)=vpi end if 30 CONTINUE IF (area.lt..170.or.area.gt.2048.)THEN jwarn(3,1)=1 END IF RETURN END c SUBROUTINE SRE4(jsplit) INCLUDE 'tnff.cmn' INTEGER i, ip, j REAL yhat REAL bs, v, vt, xtxi, temp, temp2, vmodel, & vpi, xt1, xt2 DIMENSION bs(2,7), v(1,2), vt(2,1), xt1(2,2,4), xtxi(2,2), & temp(1,2), temp2(1,1), vmodel(7), & xt2(2,2,3) INTEGER ltd,ltm,lts,lnd,lnm,lns COMMON /LTLN/ltd, ltm, lts, lnd, lnm, lns c^M c SIMPLE REGRESSION METHOD FOR HYDROLOGIC AREA 4^M c DATA stut/1.67/^M c^M DATA bs/2.63942,0.52708,2.79077,0.54491,2.86648,0.55400, & 2.94365,0.56373,2.99173,0.56995,3.03402,0.57543, & 3.11752,0.58598/ c DATA vmodel/0.25116E-01,0.23343E-01,0.24079E-01,0.26570E-01, & 0.29244E-01,0.32424E-01,0.41309E-01/ c DATA xt1/0.30529E-02,-0.11750E-02, & -0.11750E-02, 0.60037E-03, & 0.29696E-02,-0.11310E-02, & -0.11310E-02, 0.57658E-03, & 0.31910E-02,-0.12055E-02, & -0.12055E-02, 0.61393E-03, & 0.36651E-02,-0.13753E-02, & -0.13753E-02, 0.69968E-03/ DATA xt2/0.41136E-02,-0.15397E-02, & -0.15397E-02, 0.78275E-03, & 0.46181E-02,-0.17267E-02, & -0.17267E-02, 0.87726E-03, & 0.59513E-02,-0.22273E-02, & -0.22273E-02, 0.11300E-02/ c c v(1,1) = 1.0 v(1,2) = ALOG10(area) vt(1,1) = 1.0 vt(2,1) = v(1,2) DO 30 ip = 1, 7 yhat = bs(1,ip) + bs(2,ip)*v(1,2) c 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) if(jsplit.eq.1)then y1(ip,1)=yhat v1(ip,1)=vpi else y1(ip,2)=yhat v1(ip,2)=vpi end if 30 CONTINUE IF (area.lt..760.or.area.gt.2308.)THEN jwarn(4,1)=1 END IF RETURN END c SUBROUTINE MRE1(jsplit) INCLUDE 'tnff.cmn' INTEGER i, ip, j REAL yhat REAL bs, v, vt, xtxi, temp, temp2, vmodel, & vpi, xt1, xt2, m, xm, c , xmt, temp3, xc DIMENSION bs(4,7), v(1,4), vt(4,1), xt1(4,4,4), xtxi(4,4), & temp(1,4), temp2(1,1), vmodel(7),xc(1,1), & xt2(4,4,3), xm(1,3),xmt(3,1), m(3,3),temp3(1,3) c MINIMUM COVERING ELLIPSOID FOR HYDROLOGIC AREA 1^M INTEGER ltd,ltm,lts,lnd,lnm,lns COMMON /LTLN/ltd, ltm, lts, lnd, lnm, lns DATA m /1.7900,1.8012,-0.2842, & 1.8012,5.6536,-4.6668, & -0.2842,-4.6668,2489.8/ DATA xm/1.8131, 1.6496, .35897/ DATA c /8.46/ c MULTIPLE REGRESSION METHOD FOR HYDROLOGIC AREA 1^M c DATA stut/1.65/ c DATA bs/0.23497,0.79845,0.11165,4.58118, & 0.53234,0.78339,0.11443,4.33021, & 0.72718,0.77544,0.11623,4.08716, & 0.95415,0.76645,0.11723,3.77827, & 1.10860,0.76030,0.11715,3.56024, & 1.25219,0.75452,0.11656,3.35371, & 1.55710,0.74211,0.11397,2.90416/ c DATA vmodel/0.26069E-01,0.25700E-01,0.27146E-01,0.30393E-01, & 0.33622E-01,0.37415E-01,0.48149E-01/ c DATA xt1/0.93197E-01, -0.15441E-02,-0.29363E-02,-0.24614, & -0.15441E-02, 0.31751E-03, 0.32177E-03, 0.12442E-02, & -0.29363E-02, 0.32177E-03, 0.89690E-03, 0.28041E-02, & -0.24614 , 0.12442E-02, 0.28041E-02, 0.68971 , & 0.10328 ,-0.15417E-02,-0.29005E-02,-0.27460 , & -0.15417E-02, 0.33931E-03, 0.34147E-03, 0.98932E-03, & -0.29005E-02, 0.34147E-03, 0.93011E-03, 0.24319E-02, & -0.27460 , 0.98932E-03, 0.24319E-02, 0.77322 , & 0.11915 ,-0.16555E-02,-0.31044E-02,-0.31815 , & -0.16555E-02, 0.38000E-03, 0.38225E-03, 0.88139E-03, & -0.31044E-02, 0.38225E-03, 0.10263E-02, 0.23663E-02, & -0.31815 , 0.88139E-03, 0.23663E-02, 0.89838 , & 0.14548 ,-0.18852E-02,-0.35352E-02,-0.38993 , & -0.18852E-02, 0.45003E-03, 0.45336E-03, 0.80319E-03, & -0.35352E-02, 0.45336E-03, 0.12022E-02, 0.24392E-02, & -0.38993 , 0.80319E-03, 0.24392E-02, 1.1036 / DATA xt2/0.16832,-0.21028E-02,-0.39500E-02,-0.45198 , & -0.21028E-02, 0.51230E-03, 0.51692E-03, 0.77694E-03, & -0.39500E-02, 0.51692E-03, 0.13629E-02, 0.25782E-02, & -0.45198 , 0.77694E-03, 0.25782E-02, 1.2806 , & 0.19323 ,-0.23516E-02,-0.44278E-02,-0.51954 , & -0.23516E-02, 0.58141E-03, 0.58765E-03, 0.77182E-03, & -0.44278E-02, 0.58765E-03, 0.15436E-02, 0.27724E-02, & -0.51954 , 0.77182E-03, 0.27724E-02, 1.4731 , & 0.25810 ,-0.30348E-02,-0.57482E-02,-0.69499 , & -0.30348E-02, 0.76553E-03, 0.77671E-03, 0.82519E-03, & -0.57482E-02, 0.77671E-03, 0.20320E-02, 0.33973E-02, & -0.69499 , 0.82519E-03, 0.33973E-02, 1.9721 / c v(1,1) = 1.0 v(1,2) = ALOG10(area) v(1,3) = ALOG10(slp) v(1,4) = ALOG10(cf2) vt(1,1) = 1.0 vt(2,1) = v(1,2) vt(3,1) = v(1,3) vt(4,1) = v(1,4) DO 30 ip = 1, 7 yhat = bs(1,ip)+bs(2,ip)*v(1,2)+bs(3,ip)*v(1,3) & +bs(4,ip)*v(1,4) c yhat = 10**yhat c c Compute CI c DO 20 i = 1, 4 DO 10 j = 1, 4 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,4,4,1,1,4) CALL MLTPLY(temp2,temp,vt,1,4,1,1,1,4) vpi = vmodel(ip) + temp2(1,1) if(jsplit.eq.1)then y2(ip,1)=yhat v2(ip,1)=vpi else y2(ip,2)=yhat v2(ip,2)=vpi end if 30 CONTINUE IF (v(1,2).lt.-.69898.or.v(1,2).gt.4.3305)THEN jwarn(1,1)=1 END IF IF (v(1,3).lt..51718.or.v(1,3).gt.2.9778)THEN jwarn(1,2)=1 END IF IF (v(1,4).lt..31296.or.v(1,4).gt..38307)THEN jwarn(1,3)=1 END IF c c Test for extrapolation beyond MCE c xm(1,1)=v(1,2)-xm(1,1) xm(1,2)=v(1,3)-xm(1,2) xm(1,3)=v(1,4)-xm(1,3) xmt(1,1)=xm(1,1) xmt(2,1)=xm(1,2) xmt(3,1)=xm(1,3) CALL MLTPLY(temp3,xm,m,1,3,3,1,1,3) CALL MLTPLY(xc,temp3,xmt,1,3,1,1,1,3) if (xc(1,1).gt.c)then jwarn(1,4)=1 end if c RETURN END c SUBROUTINE MRE2(jsplit) INCLUDE 'tnff.cmn' INTEGER i, ip, j REAL yhat REAL bs, v, vt, xtxi, temp, temp2, vmodel, & vpi, xt1, xt2, m, xm, c , xmt, temp3, xc DIMENSION bs(3,7), v(1,3), vt(3,1), xt1(3,3,4), xtxi(3,3), & temp(1,3), temp2(1,1), vmodel(7),xc(1,1), & xt2(3,3,3), xm(1,2),xmt(2,1), m(2,2),temp3(1,2) c MINIMUM COVERING ELLIPSOID FOR HYDROLOGIC AREA 2^M INTEGER ltd,ltm,lts,lnd,lnm,lns COMMON /LTLN/ltd, ltm, lts, lnd, lnm, lns DATA m /1.4407,3.2352, & 3.2352,10.039/ DATA xm/.37285,1.5854/ DATA c /9.67/ c MULTIPLE REGRESSION METHOD FOR HYDROLOGIC AREA 2^M c DATA stut/1.66/ c DATA bs/2.02441,0.78657,0.15134, & 2.23146,0.77896,0.15835, & 2.33934,0.77555,0.16034, & 2.45547,0.77170,0.16006, & 2.53106,0.76900,0.15860, & 2.59930,0.76641,0.15654, & 2.73768,0.76081,0.15050/ c DATA vmodel/0.15891E-01,0.13787E-01,0.14618E-01,0.16996E-01, & 0.19405E-01,0.22208E-01,0.29955E-01/ c DATA xt1/0.76949E-02,-0.17969E-02,-0.33375E-02, & -0.17969E-02, 0.52148E-03, 0.69904E-03, & -0.33375E-02, 0.69904E-03, 0.17829E-02, & 0.72572E-02,-0.16767E-02,-0.30935E-02, & -0.16767E-02, 0.48396E-03, 0.64266E-03, & -0.30935E-02, 0.64266E-03, 0.16531E-02, & 0.80998E-02,-0.18632E-02,-0.34368E-02, & -0.18632E-02, 0.53648E-03, 0.71057E-03, & -0.34368E-02, 0.71057E-03, 0.18389E-02, & 0.98122E-02,-0.22513E-02,-0.41592E-02, & -0.22513E-02, 0.64721E-03, 0.85646E-03, & -0.41592E-02, 0.85646E-03, 0.22284E-02/ DATA xt2/0.11389E-01,-0.26114E-02,-0.48323E-02, & -0.26114E-02, 0.75046E-03, 0.99322E-03, & -0.48323E-02, 0.99322E-03, 0.25905E-02, & 0.13146E-01,-0.30143E-02,-0.55864E-02, & -0.30143E-02, 0.86615E-03, 0.11469E-02, & -0.55864E-02, 0.11469E-02, 0.29957E-02, & 0.17778E-01,-0.40807E-02,-0.75882E-02, & -0.40807E-02, 0.11732E-02, 0.15561E-02, & -0.75882E-02, 0.15561E-02, 0.40695E-02/ c v(1,1) = 1.0 v(1,2) = ALOG10(area) v(1,3) = ALOG10(slp) vt(1,1) = 1.0 vt(2,1) = v(1,2) vt(3,1) = v(1,3) DO 30 ip = 1, 7 yhat = bs(1,ip) + bs(2,ip)*v(1,2) +bs(3,ip)*v(1,3) c yhat = 10**yhat c c Compute CI c DO 20 i = 1, 3 DO 10 j = 1, 3 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,3,3,1,1,3) CALL MLTPLY(temp2,temp,vt,1,3,1,1,1,3) vpi = vmodel(ip) + temp2(1,1) if(jsplit.eq.1)then y2(ip,1)=yhat v2(ip,1)=vpi else y2(ip,2)=yhat v2(ip,2)=vpi end if 30 CONTINUE IF (v(1,2).lt.-.328.or.v(1,2).gt.3.4078)THEN jwarn(2,1)=1 END IF IF (v(1,3).lt..2787.or.v(1,3).gt.2.5353)THEN jwarn(2,2)=1 END IF c c Test for extrapolation beyond MCE c xm(1,1)=v(1,2)-xm(1,1) xm(1,2)=v(1,3)-xm(1,2) xmt(1,1)=xm(1,1) xmt(2,1)=xm(1,2) CALL MLTPLY(temp3,xm,m,1,2,2,1,1,2) CALL MLTPLY(xc,temp3,xmt,1,2,1,1,1,2) if (xc(1,1).gt.c)then jwarn(2,4)=1 end if c RETURN END c SUBROUTINE MRE3a(jsplit) INCLUDE 'tnff.cmn' INTEGER i, ip, j REAL yhat REAL bs, v, vt, xtxi, temp, temp2, vmodel, & vpi, xt1, xt2, m, xm, c , xmt, temp3, xc DIMENSION bs(3,7), v(1,3), vt(3,1), xt1(3,3,4), xtxi(3,3), & temp(1,3), temp2(1,1), vmodel(7),xc(1,1), & xt2(3,3,3), xm(1,2),xmt(2,1), m(2,2),temp3(1,2) c REVISED--------REVISED--------REVISED---GSLAW 28APR03 c MRE HA3 2-SEGMENT LINEAR FIT LOW END DATA (DA<1.48 LOGS) c MINIMUM COVERING ELLIPSOID FOR HYDROLOGIC AREA 3 LOW END INTEGER ltd,ltm,lts,lnd,lnm,lns COMMON /LTLN/ltd, ltm, lts, lnd, lnm, lns DATA m /6.0611,8.8985, & 8.8985,28.738/ DATA xm/.62432,1.6318/ DATA c /7.40/ c MULTIPLE REGRESSION METHOD FOR HYDROLOGIC AREA 3 LOW END^M c DATA stut/1.67/ c DATA bs/2.32341,0.81489,0.06316, & 2.51731,0.79835,0.07085, & 2.60696,0.79284,0.07804, & 2.69650,0.78852,0.08638, & 2.75190,0.78632,0.09156, & 2.80056,0.78460,0.09578, & 2.89713,0.78149,0.10228/ c DATA vmodel/0.20005E-01,0.19713E-01,0.20079E-01,0.21001E-01, & 0.21991E-01,0.23221E-01,0.26967E-01/ c DATA xt1/0.71914E-01,-0.16257E-01,-0.35509E-01, & -0.16257E-01, 0.46364E-02, 0.75099E-02, & -0.35509E-01, 0.75099E-02, 0.18177E-01, & 0.74185E-01,-0.16825E-01,-0.36619E-01, & -0.16825E-01, 0.47963E-02, 0.77765E-02, & -0.36619E-01, 0.77765E-02, 0.18738E-01, & 0.79558E-01,-0.18099E-01,-0.39232E-01, & -0.18099E-01, 0.51595E-02, 0.83586E-02, & -0.39232E-01, 0.83586E-02, 0.20066E-01, & 0.88674E-01,-0.20240E-01,-0.43668E-01, & -0.20240E-01, 0.57686E-02, 0.93373E-02, & -0.43668E-01, 0.93373E-02, 0.22321E-01/ DATA xt2/0.96701E-01,-0.22115E-01,-0.47579E-01, & -0.22115E-01, 0.63009E-02, 0.10196E-01, & -0.47579E-01, 0.10196E-01, 0.24309E-01, & 0.10559 ,-0.24185E-01,-0.51915E-01, & -0.24185E-01, 0.68880E-02, 0.11145E-01, & -0.51915E-01, 0.11145E-01, 0.26512E-01, & 0.12932 ,-0.29694E-01,-0.63520E-01, & -0.29694E-01, 0.84499E-02, 0.13679E-01, & -0.63520E-01, 0.13679E-01, 0.32410E-01/ c v(1,1) = 1.0 v(1,2) = ALOG10(area) v(1,3) = ALOG10(slp) vt(1,1) = 1.0 vt(2,1) = v(1,2) vt(3,1) = v(1,3) DO 30 ip = 1, 7 yhat = bs(1,ip) + bs(2,ip)*v(1,2) +bs(3,ip)*v(1,3) c yhat = 10**yhat c c Compute CI c DO 20 i = 1, 3 DO 10 j = 1, 3 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,3,3,1,1,3) CALL MLTPLY(temp2,temp,vt,1,3,1,1,1,3) vpi = vmodel(ip) + temp2(1,1) if(jsplit.eq.1)then y2(ip,1)=yhat v2(ip,1)=vpi else y2(ip,2)=yhat v2(ip,2)=vpi end if 30 CONTINUE IF (v(1,2).lt.-.76956.or.v(1,2).gt.3.31133)THEN jwarn(3,1)=1 END IF IF (v(1,3).lt..3263.or.v(1,3).gt.2.1206)THEN jwarn(3,2)=1 END IF c c Test for extrapolation beyond MCE c xm(1,1)=v(1,2)-xm(1,1) xm(1,2)=v(1,3)-xm(1,2) xmt(1,1)=xm(1,1) xmt(2,1)=xm(1,2) CALL MLTPLY(temp3,xm,m,1,2,2,1,1,2) CALL MLTPLY(xc,temp3,xmt,1,2,1,1,1,2) if (xc(1,1).gt.c)then jwarn(3,4)=1 end if RETURN END c SUBROUTINE MRE3b(jsplit) INCLUDE 'tnff.cmn' INTEGER i, ip, j REAL yhat REAL bs, v, vt, xtxi, temp, temp2, vmodel, & vpi, xt1, xt2, m, xm, c , xmt, temp3, xc DIMENSION bs(3,7), v(1,3), vt(3,1), xt1(3,3,4), xtxi(3,3), & temp(1,3), temp2(1,1), vmodel(7),xc(1,1), & xt2(3,3,3), xm(1,2),xmt(2,1), m(2,2),temp3(1,2) c REVISED--------REVISED--------REVISED---28APR03 c MRE HA3 2-SEGMENT LINEAR FIT HIGH END (DA>=1.48 LOGS) c MINIMUM COVERING ELLIPSOID FOR HYDROLOGIC AREA 3 HIGH END INTEGER ltd,ltm,lts,lnd,lnm,lns COMMON /LTLN/ltd, ltm, lts, lnd, lnm, lns DATA m /17.924,29.903, & 29.903,60.673/ DATA xm/2.0642,0.98761/ DATA c /5.83/ c MULTIPLE REGRESSION METHOD FOR HYDROLOGIC AREA 3 HIGH END^M c DATA stut/1.67/ c DATA bs/2.61127,0.58440,0.10209, & 2.88502,0.55795,0.06105, & 2.99114,0.55403,0.05410, & 3.07910,0.55664,0.05582, & 3.12378,0.56167,0.06111, & 3.15662,0.56827,0.06832, & 3.20319,0.58735,0.08951/ c DATA vmodel/0.12806E-01,0.13449E-01,0.15002E-01,0.17878E-01, & 0.20564E-01,0.23643E-01,0.32193E-01/ c DATA xt1/0.99052E-01,-0.26308E-01,-0.44057E-01, & -0.26308E-01, 0.72314E-02, 0.11398E-01, & -0.44057E-01, 0.11398E-01, 0.20414E-01, & 0.10611 ,-0.28179E-01,-0.47290E-01, & -0.28179E-01, 0.77470E-02, 0.12232E-01, & -0.47290E-01, 0.12232E-01, 0.21963E-01, & 0.12127 ,-0.32190E-01,-0.54068E-01, & -0.32190E-01, 0.88466E-02, 0.13981E-01, & -0.54068E-01, 0.13981E-01, 0.25128E-01, & 0.14789 ,-0.39242E-01,-0.65973E-01, & -0.39242E-01, 0.10780E-01, 0.17055E-01, & -0.65973E-01, 0.17055E-01, 0.30678E-01/ DATA xt2/0.17185 ,-0.45592E-01,-0.76698E-01, & -0.45592E-01, 0.12522E-01, 0.19826E-01, & -0.76698E-01, 0.19826E-01, 0.35680E-01, & 0.19865 ,-0.52699E-01,-0.88712E-01, & -0.52699E-01, 0.14472E-01, 0.22929E-01, & -0.88712E-01, 0.22929E-01, 0.41285E-01, & 0.27077 ,-0.71828E-01,-0.12109 , & -0.71828E-01, 0.19723E-01, 0.31297E-01, & -0.12109 , 0.31297E-01, 0.56410E-01/ c v(1,1) = 1.0 v(1,2) = ALOG10(area) v(1,3) = ALOG10(slp) vt(1,1) = 1.0 vt(2,1) = v(1,2) vt(3,1) = v(1,3) DO 30 ip = 1, 7 yhat = bs(1,ip) + bs(2,ip)*v(1,2) +bs(3,ip)*v(1,3) c yhat = 10**yhat c c Compute CI c DO 20 i = 1, 3 DO 10 j = 1, 3 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,3,3,1,1,3) CALL MLTPLY(temp2,temp,vt,1,3,1,1,1,3) vpi = vmodel(ip) + temp2(1,1) if(jsplit.eq.1)then y2(ip,1)=yhat v2(ip,1)=vpi else y2(ip,2)=yhat v2(ip,2)=vpi end if 30 CONTINUE IF (v(1,2).lt.-.76956.or.v(1,2).gt.3.31133)THEN jwarn(3,1)=1 END IF IF (v(1,3).lt..3263.or.v(1,3).gt.2.1206)THEN jwarn(3,2)=1 END IF c c Test for extrapolation beyond MCE c xm(1,1)=v(1,2)-xm(1,1) xm(1,2)=v(1,3)-xm(1,2) xmt(1,1)=xm(1,1) xmt(2,1)=xm(1,2) CALL MLTPLY(temp3,xm,m,1,2,2,1,1,2) CALL MLTPLY(xc,temp3,xmt,1,2,1,1,1,2) if (xc(1,1).gt.c)then jwarn(3,4)=1 end if RETURN END c SUBROUTINE MRE4(jsplit) INCLUDE 'tnff.cmn' INTEGER i, ip, j REAL yhat REAL bs, v, vt, xtxi, temp, temp2, vmodel, & vpi, xt1, xt2 DIMENSION bs(2,7), v(1,2), vt(2,1), xt1(2,2,4), xtxi(2,2), & temp(1,2), temp2(1,1), vmodel(7), & xt2(2,2,3) INTEGER ltd,ltm,lts,lnd,lnm,lns COMMON /LTLN/ltd, ltm, lts, lnd, lnm, lns c c MULTIPLE REGRESSION METHOD FOR HYDROLOGIC AREA 4^M c DATA stut/1.67/ c DATA bs/2.63942,0.52708,2.79077,0.54491,2.86648,0.55400, & 2.94365,0.56373,2.99173,0.56995,3.03402,0.57543, & 3.11752,0.58598/ c DATA vmodel/0.25116E-01,0.23343E-01,0.24079E-01,0.26570E-01, & 0.29244E-01,0.32424E-01,0.41309E-01/ c DATA xt1/0.30529E-02,-0.11750E-02, & -0.11750E-02, 0.60037E-03, & 0.29696E-02,-0.11310E-02, & -0.11310E-02, 0.57658E-03, & 0.31910E-02,-0.12055E-02, & -0.12055E-02, 0.61393E-03, & 0.36651E-02,-0.13753E-02, & -0.13753E-02, 0.69968E-03/ DATA xt2/0.41136E-02,-0.15397E-02, & -0.15397E-02, 0.78275E-03, & 0.46181E-02,-0.17267E-02, & -0.17267E-02, 0.87726E-03, & 0.59513E-02,-0.22273E-02, & -0.22273E-02, 0.11300E-02/ c c v(1,1) = 1.0 v(1,2) = ALOG10(area) vt(1,1) = 1.0 vt(2,1) = v(1,2) DO 30 ip = 1, 7 yhat = bs(1,ip) + bs(2,ip)*v(1,2) c 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) if(jsplit.eq.1)then y2(ip,1)=yhat v2(ip,1)=vpi else y2(ip,2)=yhat v2(ip,2)=vpi end if 30 CONTINUE IF (area.lt..760.or.area.gt.2308.)THEN jwarn(4,1)=1 END IF RETURN END c subroutine regout(isplit,imethod,ireg,ireg2) include 'tnff.cmn' INTEGER ltd, ltm, lts, lnd, lnm, lns COMMON /LTLN/ltd, ltm, lts, lnd, lnm, lns real ysav, sp, yold character*6 hregion character*1 adj dimension hregion(4), adj(7), tsav(7) dimension ysav(7),sp(7),yold(7) data hregion/' HA 1 ',' HA 2 ',' HA 3 ',' HA 4 '/ data stut/1.66/ if(isplit.eq.1) then if(imethod.eq.1)then write(*,9011)siteid,hregion(ireg),ltd,ltm,lts,lnd,lnm,lns, & area write(16,9011)siteid,hregion(ireg),ltd,ltm,lts,lnd,lnm,lns, & area end if if(imethod.eq.2)then if(ireg.eq.1)then write(*,9021)siteid,hregion(ireg),ltd,ltm,lts,lnd,lnm,lns, & area,slp,cf2 write(16,9021)siteid,hregion(ireg),ltd,ltm,lts,lnd,lnm,lns, & area,slp,cf2 end if if(ireg.eq.2.or.ireg.eq.3)then write(*,9023)siteid,hregion(ireg),ltd,ltm,lts,lnd,lnm,lns, & area,slp write(16,9023)siteid,hregion(ireg),ltd,ltm,lts,lnd,lnm,lns, & area,slp end if if(ireg.eq.4)then write(*,9011)siteid,hregion(ireg),ltd,ltm,lts,lnd,lnm,lns, & area write(16,9011)siteid,hregion(ireg),ltd,ltm,lts,lnd,lnm,lns, & area end if end if if (imethod.eq.3)then write(*,9031)siteid,hregion(ireg),ltd,ltm,lts,lnd,lnm,lns, & area,slp,cf2,p write(16,9031)siteid,hregion(ireg),ltd,ltm,lts,lnd,lnm,lns, & area,slp,cf2,p end if end if if(isplit.eq.2)then if(imethod.eq.1) then write(*,9012)siteid,hregion(ireg),psplit,hregion(ireg2), & qsplit,ltd,ltm,lts,lnd,lnm,lns, & area write(16,9012)siteid,hregion(ireg),psplit,hregion(ireg2), & qsplit,ltd,ltm,lts,lnd,lnm,lns, & area end if if(imethod.eq.2)then if(ireg.ne.1.and.ireg2.ne.1)then write(*,9024)siteid,hregion(ireg),psplit,hregion(ireg2), & qsplit,ltd,ltm,lts,lnd,lnm,lns, & area,slp write(16,9024)siteid,hregion(ireg),psplit,hregion(ireg2), & qsplit,ltd,ltm,lts,lnd,lnm,lns, & area,slp else write(*,9022)siteid,hregion(ireg),psplit,hregion(ireg2), & qsplit,ltd,ltm,lts,lnd,lnm,lns, & area,slp,cf2 write(16,9022)siteid,hregion(ireg),psplit,hregion(ireg2), & qsplit,ltd,ltm,lts,lnd,lnm,lns, & area,slp,cf2 end if end if if(imethod.eq.3)then write(*,9032)siteid,hregion(ireg),psplit,hregion(ireg2), & qsplit,ltd,ltm,lts,lnd,lnm,lns, & area,slp,cf2,p,hregion(ireg),p2sav,hregion(ireg2) write(16,9032)siteid,hregion(ireg),psplit,hregion(ireg2), & qsplit,ltd,ltm,lts,lnd,lnm,lns, & area,slp,cf2,p,hregion(ireg),p2sav,hregion(ireg2) end if end if do 30 ip=1,7 if(isplit.eq.1)then if(imethod.eq.1.or.imethod.eq.3)then ysav(ip)=10**y1(ip,1) sp(ip)=sqrt(v1(ip,1)) else ysav(ip)=10**y2(ip,1) sp(ip)=sqrt(v2(ip,1)) end if end if if(isplit.eq.2)then if(imethod.eq.1.or.imethod.eq.3)then ysav(ip)=10**(y1(ip,1)*xsplit+y1(ip,2)*zsplit) sp(ip)=sqrt(v1(ip,1)*xsplit+v1(ip,2)*zsplit) else ysav(ip)=10**(y2(ip,1)*xsplit+y2(ip,2)*zsplit) sp(ip)=sqrt(v2(ip,1)*xsplit+v2(ip,2)*zsplit) end if end if 30 continue c c check for inconsistent result c adj(1)=' ' ihave=0 do 31 ip=2,7 adj(ip)=' ' if(ysav(ip-1).gt.ysav(ip))then if(ip.eq.2) then dely1=alog10(ysav(3))-alog10(ysav(1)) dely2=alog10(ysav(3))-alog10(ysav(2)) adj(1)='*' adj(2)='*' delx1=ak(3)-ak(1) delx2=ak(3)-ak(2) slp1=dely1/delx1 slp2=dely2/delx2 slp1=amax1(slp1,0.0) slp2=amax1(slp2,0.0) slpx=(slp1+slp2)/2.0 adjy1=delx1*slpx adjy2=delx2*slpx yold(1)=ysav(1) yold(2)=ysav(2) tsav(1)=1 tsav(2)=2 ihave=2 ysav(1)=alog10(ysav(3))-adjy1 ysav(2)=ysav(1)+adjy2 ysav(1)=10**ysav(1) ysav(2)=10**ysav(2) end if if (ip.lt.7.and.ip.gt.2)then ihave=ihave+1 adj(ip)='*' dely=alog10(ysav(ip+1))-alog10(ysav(ip-1)) delx=ak(ip+1)-ak(ip-1) adjy=(ak(ip)-ak(ip-1))*dely/delx yold(ip)=ysav(ip) ysav(ip)=alog10(ysav(ip-1))+adjy ysav(ip)=10**ysav(ip) tsav(ihave)=ip end if if (ip.eq.7)then dely1=alog10(ysav(7))-alog10(ysav(5)) dely2=alog10(ysav(6))-alog10(ysav(5)) ihave=ihave+2 adj(6)='*' adj(7)='*' delx1=ak(7)-ak(5) delx2=ak(6)-ak(5) slp1=dely1/delx1 slp2=dely2/delx2 slp1=amax1(slp1,0.0) slp2=amax1(slp2,0.0) slpx=(slp1+slp2)/2.0 adjy1=delx1*slpx adjy2=delx2*slpx yold(6)=ysav(6) yold(7)=ysav(7) ysav(7)=alog10(ysav(5))+adjy1 ysav(6)=alog10(ysav(5))+adjy2 ysav(7)=10**ysav(7) ysav(6)=10**ysav(6) tsav(ihave-1)=6 tsav(ihave)=7 end if c tsave=ip c yold(ip)=ysav(ip) c if(ihave.eq.1)then c ysav(ip)=alog10(ysav(ip-1))+adjy c ysav(ip)=10**ysav(ip) c end if c if(ihave.eq.2)then c tsv2=ip-1 c yold(ip-1)=ysav(ip-1) c ysav(ip)=alog10(ysav(ip-2))+adjy1 c ysav(ip-1)=alog10(ysav(ip-2))+adjy2 c ysav(ip)=10**ysav(ip) c ysav(ip-1)=10**ysav(ip-1) c end if end if 31 continue do 32 ip=1,7 yhat=ysav(ip) asep=sp(ip) 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,adj(ip),sepl,sepu,cl,cu write(16,9000)pklab(ip),yhat,adj(ip),sepl,sepu,cl,cu 9000 format(2x,a8,f12.1,a1,2f11.1,2f12.1) 32 CONTINUE if(ihave.gt.0)then write(16,8002) write(*,8002) do 321 ih=1,ihave write(16,8001)pklab(tsav(ih)),yold(tsav(ih)) write(*,8001)pklab(tsav(ih)),yold(tsav(ih)) 321 continue end if c c check warnings c do 41 jsplit=1,isplit if(jsplit.eq.1)then jreg=ireg else jreg=ireg2 end if if(jwarn(jreg,1).eq.1)then write(*,9001)hregion(jreg) write(16,9001)hregion(jreg) end if if(jwarn(jreg,2).eq.1)then write(*,9002)hregion(jreg) write(16,9002)hregion(jreg) end if if(jwarn(jreg,3).eq.1)then write(*,9003)hregion(jreg) write(16,9003)hregion(jreg) end if if(jwarn(jreg,4).eq.1)then write(*,9004)hregion(jreg) write(16,9004)hregion(jreg) end if 41 continue return 8001 FORMAT (2x,a8,f12.1) 8002 FORMAT (/,' * Discharge adjusted for consistency with other', &' peak estimates.',/,' UNADJUSTED DISCHARGES', & /,' RI Discharge (cfs)') 9001 FORMAT (/,' WARNING-- In hydrologic area',a6, & ' Drainage area out of range') 9002 FORMAT (/,' WARNING-- In hydrologic area',a6, & ' Channel slope out of range') 9003 FORMAT (/,' WARNING-- In hydrologic area',a6, & ' Climate factor out of range') 9004 FORMAT (/,' WARNING-- In hydrologic area',a6, & ' Extrapolation beyond observed data') 9023 FORMAT (//,' MULTIVARIABLE REGIONAL-REGRESSION EQUATION (MRE)', & ' METHOD FOR TENNESSEE',/, & ' Flood frequency estimates for:',/,1x,a65,/, & ' Hydrologic Area:',a6,' LAT: ',3i4,3x,' LNG: ',3i4,/, & ' Explanatory variables:',/, & ' Contributing drainage area:',f8.2,' square miles',/, & ' Channel slope:',f7.2,' ft/mi',/, & ' RI DISCHARGE - SE (%) + SE (%) 90%', & ' PRED. INTERVAL',/,t12,' (cfs)') 9021 FORMAT (//,' MULTIVARIABLE REGIONAL-REGRESSION EQUATION (MRE)', & ' METHOD FOR TENNESSEE',/, & ' Flood frequency estimates for:',/,1x,a65,/, & ' Hydrologic Area:',a6,' LAT: ',3i4,3x,' LNG: ',3i4,/, & ' Explanatory variables:',/, & ' Contributing drainage area:',f8.2,' square miles',/, & ' Channel slope:',f7.2,' ft/mi',/, & ' Climate factor:',f7.2,/, & ' RI DISCHARGE - SE (%) + SE (%) 90%', & ' PRED. INTERVAL',/,t12,' (cfs)') 9011 FORMAT (/,' TDOT Version 2.0.3',//, & ' SINGLE-VARIABLE REGIONAL-REGRESSION EQUATION (SRE)', & ' METHOD FOR TENNESSEE',/, & ' Flood frequency estimates for:',/,1x,a65,/, & ' Hydrologic Area:',a6,' LAT: ',3i4,3x,' LNG: ',3i4,/, & ' Explanatory variable:',/, & ' Contributing drainage area:',f8.2,' square miles',/, & ' RI DISCHARGE - SE (%) + SE (%) 90%', & ' PRED. INTERVAL',/,t12,' (cfs)') 9012 FORMAT (/,' TDOT Version 2.0.3',//, & ' SINGLE-VARIABLE REGIONAL-REGRESSION EQUATION (SRE)', & ' METHOD FOR TENNESSEE',/, & ' Flood frequency estimates for:',/,1x,a65,/, & ' Hydrologic Areas (percent):',a6,'(',f5.1,')', & a6,'(',f5.1,')',/,' LAT: ',3i4,3x,' LNG: ',3i4,/, & ' Explanatory variable:',/, & ' Contributing drainage area:',f8.2,' square miles',/, & ' RI DISCHARGE - SE (%) + SE (%) 90%', & ' PRED. INTERVAL',/,t12,' (cfs)') 9031 FORMAT (//,' REGION-OF-INFLUENCE (ROI) METHOD FOR TENNESSEE',/, & ' Flood frequency estimates for:',/,1x,a65,/, & ' Hydrologic Area:',a6,' LAT: ',3i4,3x,' LNG: ',3i4,/, & ' Explanatory variables:',/, & ' Contributing drainage area:',f8.2,' square miles',/, & ' Channel slope:',f7.2,' ft/mi',/, & ' Climate factor:',f7.2,/, & ' Log(Physiographic factor):',f9.4,/, & ' RI DISCHARGE - SE (%) + SE (%) 90%', & ' PRED. INTERVAL',/,t12,' (cfs)') 9032 FORMAT (//,' REGION-OF-INFLUENCE (ROI) METHOD FOR TENNESSEE',/, & ' Flood frequency estimates for:',/,1x,a65,/, & ' Hydrologic Areas (percent):',a6,'(',f5.1,')', & a6,'(',f5.1,')',/,' LAT: ',3i4,3x,' LNG: ',3i4,/, & ' Explanatory variables:',/, & ' Contributing drainage area:',f8.2,' square miles',/, & ' Channel slope:',f7.2,' ft/mi',/, & ' Climate factor:',f7.2,/, & ' Log(Physiographic Factor):',f8.3,'(',a6,')',f10.3,'(',a6,')', & /,' RI DISCHARGE - SE (%) + SE (%) 90%', & ' PRED. INTERVAL',/,t12,' (cfs)') 9022 FORMAT (//,' MULTIVARIABLE REGIONAL-REGRESSION EQUATION (MRE)', & ' METHOD FOR TENNESSEE',/, & ' Flood frequency estimates for:',/,1x,a65,/, & ' Hydrologic Areas (percent):',a6,'(',f5.1,')', & a6,'(',f5.1,')',/,' LAT: ',3i4,3x,' LNG: ',3i4,/, & ' Explanatory variables:',/, & ' Contributing drainage area:',f8.2,' square miles',/, & ' Channel slope:',f7.2,' ft/mi',/, & ' Climate factor:',f7.2,/, & ' RI DISCHARGE - SE (%) + SE (%) 90%', & ' PRED. INTERVAL',/,t12,' (cfs)') 9024 FORMAT (//,' MULTIVARIABLE REGIONAL-REGRESSION EQUATION (MRE)', & ' METHOD FOR TENNESSEE',/, & ' Flood frequency estimates for:',/,1x,a65,/, & ' Hydrologic Areas (percent):',a6,'(',f5.1,')', & a6,'(',f5.1,')',/,' LAT: ',3i4,3x,' LNG: ',3i4,/, & ' Explanatory variables:',/, & ' Contributing drainage area:',f8.2,' square miles',/, & ' Channel slope:',f7.2,' ft/mi',/, & ' RI DISCHARGE - SE (%) + SE (%) 90%', & ' PRED. INTERVAL',/,t12,' (cfs)') end