C------------------------------------------------------------------------------
C THIS PROGRAM CALCULATES THE MVPMC METHOD INCLUDING THE WATER SURFACE
C SLOPE AND THE VEDERNIKOV NUMBER IN THE HYDRAULIC DIFFUSIVITY
C DEVELOPED 21FEB97 BY V. M. PONCE AND A. LUGO
C------------------------------------------------------------------------------
PROGRAM SYS97_LOOPEDMC
PARAMETER (VERSION=1.11,IVERSIONDATE=971024)
PARAMETER (MR=160,MT=640)
PARAMETER (LU5=5,LU6=6,LU7=7,LU8=8,LU9=9,LU10=10,LU11=11)
DIMENSION Q3(0:MR,0:MT),D3(0:MR,0:MT)
DIMENSION Q4(0:MR,0:MT),D4(0:MR,0:MT)
DIMENSION Q5(0:MR,0:MT),D5(0:MR,0:MT)
DIMENSION Q6(0:MR,0:MT),D6(0:MR,0:MT)
DIMENSION TIME(-1:MT)
DIMENSION WSSL(MR,0:MT),AWSSL(MR,MT)
DIMENSION Q5RAT(MT),D5RAT(MT)
DIMENSION Q5LOWER(MT),Q5UPPER(MT),D5LOWER(MT),D5UPPER(MT)
DIMENSION QBELLY(0:100)
DIMENSION DBELLYLOWER(100),DBELLYUPPER(100),BELLY(100)
COMMON /PARAIN/ QREF,DREF,SLOPE,ALPHA,BETA,REBETA,DT,DX,GR,BEONESQ
COMMON /PARAOUT/ COU,REY
DOUBLE PRECISION Q3,D3,Q4,D4,Q5,D5,Q6,D6,WSSL,AWSSL
DOUBLE PRECISION QTRIAL,DTRIAL,QTOL,DEDIFF,DTOL
DOUBLE PRECISION QREF,DREF
DOUBLE PRECISION SLOPE,ALPHA,BETA,REBETA
DOUBLE PRECISION DT,DX,GR,BEONESQ
DOUBLE PRECISION COU,REY
DOUBLE PRECISION C0,C1,C2
CHARACTER*15 RESOL,PEAK,FLOOD,SLO,DEP,UNIT,DRO,REA,BASE,MANN
OPEN(LU7,FILE='sys97_loopedmc.dat',STATUS='UNKNOWN')
OPEN(LU8,FILE='sys97_loopedmc.res',STATUS='UNKNOWN')
OPEN(LU9,FILE='sys97_loopedmc.aux',STATUS='UNKNOWN')
OPEN(LU10,FILE='sys97_loopedmc.deb',STATUS='UNKNOWN')
OPEN(LU11,FILE='sys97_loopedmc.plo',STATUS='UNKNOWN')
C------------------------------------------------------------------------------
C DATA STATEMENTS
C------------------------------------------------------------------------------
DATA PI /3.1415926/
DATA JPRINT,NPRINT /2*5/
DATA SUMIN23,SUMOUT23,SIMIN43,SUMOUT43 /4*0./
DATA SUMIN24,SUMOUT24,SIMIN44,SUMOUT44 /4*0./
DATA SUMIN25,SUMOUT25,SIMIN45,SUMOUT45 /4*0./
DATA SUMIN26,SUMOUT26,SIMIN46,SUMOUT46 /4*0./
C------------------------------------------------------------------------------
C BATCH DATA INPUT
C------------------------------------------------------------------------------
C IRESOL= RESOLUTION LEVEL, EITHER 1, 2, OR 3 ONLY, INTEGER.
C IN U.S. CUSTOMARY UNITS:
C IRESOL= 1: DX= 25 MILES; DT= 6 HOURS.
C IRESOL= 2: DX= 12.5 MILES; DT= 3 HOURS.
C IRESOL= 3: DX= 6.25 MILES; DT= 1.5 HOURS.
C IN SI UNITS:
C IRESOL= 1: DX= 10 KM; DT= 3 HOURS.
C IRESOL= 2: DX= 5 KM; DT= 1.5 HOURS.
C IRESOL= 3: DX= 2.5 KM; DT= 0.75 HOURS.
C IQPRATIO= RATIO OF INITIAL PEAK TO BASE FLOW, INTEGER.
C PERIOD= PERIOD OF SINUSOIDAL FLOOD WAVE, IN HOURS, REAL.
C ALWE= WEIGHTING FACTOR FOR SLOPE
C ALDE= WEIGHTING FACTOR FOR DEPTH
C INDU= UNIT INDICATOR, 0 FOR U.S. CUSTOMARY; 1 FOR SI.
C------------------------------------------------------------------------------
READ(LU7,*) RESOL,IRESOL
READ(LU7,*) PEAK,IQPRATIO
READ(LU7,*) FLOOD,PERIOD
READ(LU7,*) SLO,ALWE
READ(LU7,*) DEP,DEWE
READ(LU7,*) UNIT,INDU
READ(LU7,*) DRO,DROP
READ(LU7,*) REA,RLENGTH
READ(LU7,*) BASE,QBASE
READ(LU7,*) MANN,XN
C------------------------------------------------------------------------------
C TESTING FOR INCONSISTENT INPUT
C------------------------------------------------------------------------------
IF(IRESOL.LT.1.OR.IRESOL.GT.3) THEN
WRITE(LU6,'(1X,A)')
1'ERROR IN RESOLUTION. PLEASE CHECK INPUT. STOP.'
STOP
ENDIF
IF(INDU.EQ.0) THEN
IF(IRESOL.EQ.1) THEN
DXM= 25.
DTH= 6.
ELSEIF(IRESOL.EQ.2) THEN
DXM= 12.5
DTH= 3.
ELSEIF(IRESOL.EQ.3) THEN
DXM= 6.25
DTH= 1.5
ENDIF
ELSEIF(INDU.EQ.1) THEN
IF(IRESOL.EQ.1) THEN
DXM= 10.
DTH= 3
ELSEIF(IRESOL.EQ.2) THEN
DXM= 5.
DTH= 1.5
ELSEIF(IRESOL.EQ.3) THEN
DXM= 2.5
DTH= 0.75
ENDIF
ENDIF
IF(IQPRATIO.LT.2) THEN
WRITE(LU6,'(1X,A)')
1'ERROR IN PEAK-TO-BASEFLOW RATIO. PLEASE CHECK INPUT. STOP.'
STOP
ELSEIF(IQPRATIO.GT.10) THEN
WRITE(LU6,'(1X,A)')
1'ERROR IN PEAK-TO-BASEFLOW RATIO. PLEASE CHECK INPUT. STOP.'
STOP
ENDIF
IF(ALWE.LT.0.5) THEN
WRITE(LU6,'(1X,A)')
1'ERROR IN SLOPE WEIGHTING. PLEASE CHECK INPUT. STOP.'
STOP
ENDIF
IF(DEWE.LT.0.5) THEN
WRITE(LU6,'(1X,A)')
1'ERROR IN DEPTH WEIGHTING. PLEASE CHECK INPUT. STOP.'
STOP
ENDIF
IF(INDU.LT.0) THEN
WRITE(LU6,'(1X,A)')
1'ERROR IN UNIT INDICATOR. PLEASE CHECK INPUT. STOP.'
STOP
ELSEIF(INDU.GT.1) THEN
WRITE(LU6,'(1X,A)')
1'ERROR IN UNIT INDICATOR. PLEASE CHECK INPUT. STOP.'
STOP
ENDIF
C------------------------------------------------------------------------------
C PROBLEM CONSTANT INPUT
C------------------------------------------------------------------------------
IF(INDU.EQ.0) THEN
GR= 32.17
CMAN= 1.486
XN= 0.02972
BETA= 1.666667
KTTRATIO= 5
QTOL= 0.01
DTOL= 0.001
ELSEIF(INDU.EQ.1) THEN
GR= 9.81
CMAN= 1.
BETA= 1.666667
KTTRATIO= 5
QTOL= 0.01
DTOL= 0.001
ENDIF
C------------------------------------------------------------------------------
WRITE(LU6,'(1X,79(1H-))')
WRITE(LU8,'(1X,79(1H-))')
WRITE(LU6,'(A,/,A,F8.2,/,A,I6)')
1' SYS978: MUSKINGUM-CUNGE ROUTING MODEL WITH VARIABLE SLOPE',
2' VERSION:',VERSION,' DATE: ',IVERSIONDATE
WRITE(LU8,'(A,/,A,F6.2,/,A,I6)')
1' SYS978: MUSKINGUM-CUNGE ROUTING MODEL WITH VARIABLE SLOPE',
2' VERSION ',VERSION,' DATE ',IVERSIONDATE
WRITE(LU6,'(1X,79(1H-))')
WRITE(LU8,'(1X,79(1H-))')
WRITE
*(LU6,'(1X,A,I3,/,1X,A,I3,/,1X,A,F6.0,/,1X,A,F6.2,/,1X,A,F6.2,
*/,1X,A,I3,/,1X,A,F6.2,/,1X,A,F6.0,/,1X,A,F6.0,/,1X,A,F10.5)')
1'RESOLUTION LEVEL = ',IRESOL,
2'PEAK-TO-BASEFLOW RATIO = ',IQPRATIO,
3'FLOOD WAVE PERIOD (HR) = ',PERIOD,
4'SLOPE WEIGHT = ',ALWE,
5'DEPTH WEIGHT = ',DEWE,
6'INDICATOR OF UNITS = ',INDU,
7'DROP (FT/MILE OR M/KM) = ',DROP,
8'REACH LENGTH (MI OR KM) = ',RLENGTH,
9'BASEFLOW (CFS OR M3/S) = ',QBASE,
*'MANNING N = ',XN
WRITE
*(LU8,'(1X,A,I3,/,1X,A,I3,/,1X,A,F6.0,/,1X,A,F6.2,/,1X,A,F6.2,
*/,1X,A,I3,/,1X,A,F6.2,/,1X,A,F6.0,/,1X,A,F6.0,/,1X,A,F10.5)')
1'RESOLUTION LEVEL = ',IRESOL,
2'PEAK-TO-BASEFLOW RATIO = ',IQPRATIO,
3'FLOOD WAVE PERIOD (HR) = ',PERIOD,
4'SLOPE WEIGHT = ',ALWE,
5'DEPTH WEIGHT = ',DEWE,
6'INDICATOR OF UNITS = ',INDU,
7'DROP (FT/MI OR M/KM) = ',DROP,
8'REACH LENGTH (MI OR KM) = ',RLENGTH,
9'BASEFLOW (CFS OR M3/S) = ',QBASE,
*'MANNING N = ',XN
WRITE
*(LU11,'(1X,A,I3,/,1X,A,I3,/,1X,A,F6.0,/,1X,A,F6.2,/,1X,A,F6.2,
*/,1X,A,I3,/,1X,A,F6.2,/,1X,A,F6.0,/,1X,A,F6.0,/,1X,A,F10.5)')
1'RESOLUTION LEVEL = ',IRESOL,
2'PEAK-TO-BASEFLOW RATIO = ',IQPRATIO,
3'FLOOD WAVE PERIOD (HR) = ',PERIOD,
4'SLOPE WEIGHT = ',ALWE,
5'DEPTH WEIGHT = ',DEWE,
6'INDICATOR OF UNITS = ',INDU,
7'DROP (FT/MI OR M/KM) = ',DROP,
8'REACH LENGTH (MI OR KM) = ',RLENGTH,
9'BASEFLOW (CFS OR M3/S) = ',QBASE,
*'MANNING N = ',XN
C------------------------------------------------------------------------------
C SETTING INITIAL AND BOUNDARY CONDITIONS
C------------------------------------------------------------------------------
BEONESQ= (BETA - 1.)**2
QPEAK= QBASE*IQPRATIO
TTIME= PERIOD*KTTRATIO
ALWER= 1. - ALWE
DEWER= 1. - DEWE
NR= RLENGTH/DXM
NT= TTIME/DTH
IF(NR.GT.MR) THEN
WRITE(LU6,'(1X,A)')
1'NR GREATER THAN MR. PLEASE CHECK ARRAY SIZE. STOP.'
STOP
ENDIF
IF(NT.GT.MT) THEN
WRITE(LU6,'(1X,A)')
1'NT GREATER THAN MT. PLEASE CHECK ARRAY SIZE. STOP.'
STOP
ENDIF
IF(INDU.EQ.0) THEN
DX= DXM*5280.
SLOPE0= DROP/5280.
ELSE
DX= DXM*1000.
SLOPE0= DROP/1000.
ENDIF
DT= DTH*3600.
TIME(-1)= -DTH
ALPHA0= (CMAN/XN)*SQRT(SLOPE0)
REBETA= 1./BETA
DBASE= (QBASE/ALPHA0)**REBETA
C------------------------------------------------------------------------------
C CALCULATION OF INITIAL CONDITION (BASEFLOW)
C------------------------------------------------------------------------------
DO 10 J= 0,NR
Q3(J,0)= QBASE
D3(J,0)= DBASE
Q4(J,0)= QBASE
D4(J,0)= DBASE
Q5(J,0)= QBASE
D5(J,0)= DBASE
Q6(J,0)= QBASE
D6(J,0)= DBASE
10 ENDDO
DO 15 J= 1,NR
WSSL(J,0)= SLOPE0
15 ENDDO
C------------------------------------------------------------------------------
C CALCULATION OF BOUNDARY CONDITION (INFLOW HYDROGRAPH)
C------------------------------------------------------------------------------
QA= 0.5*(QPEAK + QBASE)
QD= 0.5*(QPEAK - QBASE)
DO 20 N= 1,NT
Q3(0,N)= QBASE
D3(0,N)= DBASE
Q4(0,N)= QBASE
D4(0,N)= DBASE
Q5(0,N)= QBASE
D5(0,N)= DBASE
Q6(0,N)= QBASE
D6(0,N)= DBASE
20 ENDDO
T= -DTH
N= -1
DO 25 WHILE (T.LT.PERIOD)
T= T + DTH
N= N + 1
Q3(0,N)= QA - QD*COS(2*PI*T/PERIOD)
D3(0,N)= (Q3(0,N)/ALPHA0)**REBETA
Q4(0,N)= Q3(0,N)
D4(0,N)= D3(0,N)
Q5(0,N)= Q3(0,N)
D5(0,N)= D3(0,N)
Q6(0,N)= Q3(0,N)
D6(0,N)= D3(0,N)
25 ENDDO
C------------------------------------------------------------------------------
WRITE(LU9,*) 'IRESOL, IQPRATIO, PERIOD ',IRESOL,IQPRATIO,PERIOD
WRITE(LU9,*) 'INDU ',INDU
C------------------------------------------------------------------------------
C MAIN LOOP
C------------------------------------------------------------------------------
DO 60 J= 1,NR
DO 50 N= 1,NT
C------------------------------------------------------------------------------
C THREE-POINT COMPUTATION (FIRST GUESS OF FOUR-POINT)
C------------------------------------------------------------------------------
SLOPE= SLOPE0
ALPHA= ALPHA0
QREF= 0.3333333*(Q3(J-1,N) + Q3(J-1,N-1) + Q3(J,N-1))
DREF= (QREF/ALPHA)**REBETA
CALL PARAME
CALL COEFF(C0,C1,C2)
Q3(J,N)= C0*Q3(J-1,N) + C1*Q3(J-1,N-1) + C2*Q3(J,N-1)
D3(J,N)= (Q3(J,N)/ALPHA)**REBETA
QTRIAL= Q3(J,N)
DTRIAL= D3(J,N)
IFLAG= 0
ITRIAL= 0
C------------------------------------------------------------------------------
C FOUR-POINT COMPUTATION (BASED ON BED SLOPE)
C------------------------------------------------------------------------------
DO 30 WHILE (IFLAG.EQ.0)
ITRIAL= ITRIAL + 1
DREF=
10.5*(DEWER*(D4(J-1,N-1) + D4(J,N-1)) + DEWE*(D4(J-1,N) + DTRIAL))
QREF= ALPHA*DREF**BETA
CALL PARAME
CALL COEFF(C0,C1,C2)
Q4(J,N)= C0*Q4(J-1,N) + C1*Q4(J-1,N-1) + C2*Q4(J,N-1)
D4(J,N)= (Q4(J,N)/ALPHA)**REBETA
IF(ABS(Q4(J,N)-QTRIAL).LT.QTOL) THEN
IFLAG= 1
ENDIF
QTRIAL= Q4(J,N)
DTRIAL= D4(J,N)
WRITE(LU10,'(A,3I5,2F9.3,F10.6)')
1'J,N,ITRIAL,QTRIAL,DTRIAL,SLOPE= ',
2J,N,ITRIAL,QTRIAL,DTRIAL,SLOPE
30 ENDDO
WRITE (LU9,'(A,2I5,2F10.2)') 'J,N,COU,REY',J,N,COU,REY
C------------------------------------------------------------------------------
C FIVE-POINT COMPUTATION (BASED ON WATER SURFACE SLOPE)
C------------------------------------------------------------------------------
JFLAG= 0
JTRIAL= 0
DO 40 WHILE (JFLAG.EQ.0)
JTRIAL= JTRIAL + 1
DEDIFF= DTRIAL - D5(J-1,N)
WSSL(J,N)= SLOPE0 - DEDIFF/DX
AWSSL(J,N) = ALWE*WSSL(J,N) + ALWER*WSSL(J,N-1)
SLOPE= AWSSL(J,N)
IF(SLOPE.LT.0.) THEN
WRITE(LU6,*) ' 5-POINT: SLOPE IS LESS THAN ZERO. STOP. J,N ',J,N
STOP
ENDIF
ALPHA= (CMAN/XN)*SQRT(SLOPE)
DREF=
10.5*(DEWER*(D5(J-1,N-1) + D5(J,N-1)) + DEWE*(D5(J-1,N) + DTRIAL))
QREF= ALPHA*DREF**BETA
CALL PARAME
CALL COEFF(C0,C1,C2)
Q5(J,N)= C0*Q5(J-1,N) + C1*Q5(J-1,N-1) + C2*Q5(J,N-1)
D5(J,N)= (Q5(J,N)/ALPHA)**REBETA
IF(ABS(Q5(J,N)-QTRIAL).LT.QTOL) THEN
JFLAG= 1
ENDIF
QTRIAL= Q5(J,N)
DTRIAL= D5(J,N)
WRITE(LU10,'(A,3I5,2F9.3,1X,F10.6)')
1'J,N,JTRIAL,QTRIAL,DTRIAL,SLOPE= ',
2J,N,JTRIAL,QTRIAL,DTRIAL,SLOPE
40 ENDDO
IF(MOD(J,JPRINT).EQ.0) THEN
IF(MOD(N,NPRINT).EQ.0) THEN
C WRITE(LU6,'(1X,A,4I5)') 'STILL WORKING, J,N ',J,N
ENDIF
ENDIF
50 ENDDO
60 ENDDO
C------------------------------------------------------------------------------
C PRINTING RESULTS OF FOUR AND FIVE-POINT METHODS
C------------------------------------------------------------------------------
WRITE(LU6,'(1X,79(1H-))')
WRITE(LU8,'(1X,79(1H-))')
WRITE(LU11,'(1X,79(1H-))')
WRITE(LU6,*)
1' STEP TIME INFLOW OUTFLOW4 DEPTH4 OUTFLOW5
2 DEPTH5'
WRITE(LU8,*)
1' STEP TIME INFLOW OUTFLOW4 DEPTH4 OUTFLOW5
2 DEPTH5'
WRITE(LU11,*)
1' STEP TIME INFLOW OUTFLOW4 DEPTH4 OUTFLOW5
2 DEPTH5'
WRITE(LU6,'(1X,79(1H-))')
WRITE(LU8,'(1X,79(1H-))')
WRITE(LU11,'(1X,79(1H-))')
DO 65 N= 0,NT
TIME(N) = TIME(N-1) + DTH
WRITE(LU6,'(I6,F8.2,5F12.3)')
1N,TIME(N),Q4(0,N),Q4(NR,N),D4(NR,N),Q5(NR,N),D5(NR,N)
WRITE(LU8,'(I6,F8.2,5F12.3)')
1N,TIME(N),Q4(0,N),Q4(NR,N),D4(NR,N),Q5(NR,N),D5(NR,N)
WRITE(LU11,'(I5,1H,,F8.2,5(1H,,F11.3))')
1N,TIME(N),Q4(0,N),Q4(NR,N),D4(NR,N),Q5(NR,N),D5(NR,N)
65 ENDDO
WRITE(LU6,'(1X,79(1H-))')
WRITE(LU8,'(1X,79(1H-))')
WRITE(LU11,'(1X,79(1H-))')
NCOUNT= 0
DO 70 N= 1,NT-1,2
SUMIN43= SUMIN43 + Q3(0,N) - QBASE
SUMOUT43= SUMOUT43 + Q3(NR,N) - QBASE
SUMIN44= SUMIN43
SUMOUT44= SUMOUT44 + Q4(NR,N) - QBASE
SUMIN45= SUMIN43
SUMOUT45= SUMOUT45 + Q5(NR,N) - QBASE
IF((Q5(NR,N) - QBASE).GT.0.001) THEN
NCOUNT= NCOUNT + 1
ENDIF
70 ENDDO
DO 75 N= 2,NT-2,2
SUMIN23= SUMIN23 + Q3(0,N) - QBASE
SUMOUT23= SUMOUT23 + Q3(NR,N) - QBASE
SUMIN24= SUMIN23
SUMOUT24= SUMOUT24 + Q4(NR,N) - QBASE
SUMIN25= SUMIN23
SUMOUT25= SUMOUT25 + Q5(NR,N) - QBASE
IF((Q5(NR,N) - QBASE).GT.0.001) THEN
NCOUNT= NCOUNT + 1
ENDIF
75 ENDDO
VOLIN3=
1(DT/3)*(Q3(0,0) - 2*QBASE + 4*SUMIN43 + 2*SUMIN23 + Q3(0,NT))
VOLOUT3=
1(DT/3)*(Q3(NR,0) - 2*QBASE + 4*SUMOUT43 + 2*SUMOUT23 + Q3(NR,NT))
VOLIN4= VOLIN3
VOLOUT4=
1(DT/3)*(Q4(NR,0) - 2*QBASE + 4*SUMOUT44 + 2*SUMOUT24 + Q4(NR,NT))
VOLIN5= VOLIN3
VOLOUT5=
1(DT/3)*(Q5(NR,0) - 2*QBASE + 4*SUMOUT45 + 2*SUMOUT25 + Q5(NR,NT))
PERC3= 100.*(VOLOUT3/VOLIN3)
PERC4= 100.*(VOLOUT4/VOLIN4)
PERC5= 100.*(VOLOUT5/VOLIN5)
WRITE(LU6,'(1X,A,F10.3)')
1'PERCENTAGE MASS CONSERVATION (3-POINT METHOD)= ',PERC3
WRITE(LU8,'(1X,A,F10.3)')
1'PERCENTAGE MASS CONSERVATION (3-POINT METHOD)= ',PERC3
WRITE(LU6,'(1X,A,F10.3)')
1'PERCENTAGE MASS CONSERVATION (4-POINT METHOD)= ',PERC4
WRITE(LU8,'(1X,A,F10.3)')
1'PERCENTAGE MASS CONSERVATION (4-POINT METHOD)= ',PERC4
WRITE(LU6,'(1X,A,F10.3)')
1'PERCENTAGE MASS CONSERVATION (5-POINT METHOD)= ',PERC5
WRITE(LU8,'(1X,A,F10.3)')
1'PERCENTAGE MASS CONSERVATION (5-POINT METHOD)= ',PERC5
C------------------------------------------------------------------------------
C COMPUTATION OF SIZE OF BELLY IN LOOPED RATING, 5-POINT METHOD
C------------------------------------------------------------------------------
NK= 0
DO 200 N= 1,NT
Q5TEMP= Q5(NR,N) - QBASE
IF(Q5TEMP.GE.0.001) THEN
NK= NK + 1
Q5RAT(NK)= Q5(NR,N)
D5RAT(NK)= D5(NR,N)
ENDIF
200 ENDDO
C WRITE(LU6,*) 'NK =',NK
C WRITE(LU8,*) 'NK =',NK
Q5RATMAX= QBASE
MFLAG= 0
DO 210 WHILE(MFLAG.EQ.0)
DO 205 L= 1,NK
IF(Q5RAT(L).GT.Q5RATMAX) THEN
Q5RATMAX= Q5RAT(L)
NPL= L
ELSE
MFLAG= 1
ENDIF
205 ENDDO
210 ENDDO
C WRITE(LU6,*) 'Q5RATMAX =',Q5RATMAXC
C WRITE(LU8,*) 'Q5RATMAX =',Q5RATMAX
NPU= NK - NPL + 1
DO 220 L= 1,NPL
Q5LOWER(L)= Q5RAT(L)
D5LOWER(L)= D5RAT(L)
220 ENDDO
DO 230 K= 1,NPU
L= NK - K + 1
Q5UPPER(K) = Q5RAT(L)
D5UPPER(K) = D5RAT(L)
230 ENDDO
C WRITE(LU6,*) NPL,(Q5LOWER(L),L=1,NPL)
C WRITE(LU6,*) NPU,(Q5UPPER(L),L=1,NPU)
C WRITE(LU6,*) NPL,(D5LOWER(L),L=1,NPL)
C WRITE(LU6,*) NPU,(D5UPPER(L),L=1,NPU)
C WRITE(LU8,*) NPL,(Q5LOWER(L),L=1,NPL)
C WRITE(LU8,*) NPU,(Q5UPPER(L),L=1,NPU)
C WRITE(LU8,*) NPL,(D5LOWER(L),L=1,NPL)
C WRITE(LU8,*) NPU,(D5UPPER(L),L=1,NPU)
QINTERVAL= 0.01*(Q5RATMAX - QBASE)
QBELLY(0)= QBASE
DO 240 L= 1,100
QBELLY(L)= QBELLY(L-1) + QINTERVAL
XL= QBELLY(L)
CALL PWLIN(Q5LOWER,D5LOWER,NPL,XL,YL)
DBELLYLOWER(L)= YL
CALL PWLIN(Q5UPPER,D5UPPER,NPU,XL,YL)
DBELLYUPPER(L)= YL
BELLY(L)= DBELLYUPPER(L) - DBELLYLOWER(L)
240 ENDDO
C WRITE(LU6,*) 'QBELLY(L)= ',(QBELLY(L),L=1,100)
C WRITE(LU8,*) 'QBELLY(L)= ',(QBELLY(L),L=1,100)
C WRITE(LU6,*) 'DBELLYUPPER(L)= ',(DBELLYUPPER(L),L=1,100)
C WRITE(LU8,*) 'DBELLYUPPER(L)= ',(DBELLYUPPER(L),L=1,100)
C WRITE(LU6,*) 'DBELLYLOWER(L)= ',(DBELLYLOWER(L),L=1,100)
C WRITE(LU8,*) 'DBELLYLOWER(L)= ',(DBELLYLOWER(L),L=1,100)
C WRITE(LU6,*) 'BELLY(L)= ',(BELLY(L),L=1,100)
C WRITE(LU8,*) 'BELLY(L)= ',(BELLY(L),L=1,100)
BELLYMAX= 0.
DO 250 L= 1,100
IF(BELLY(L).GT.BELLYMAX) THEN
BELLYMAX= BELLY(L)
LMAX= L
ENDIF
250 ENDDO
DEPTHBELLY= 0.5*(DBELLYUPPER(LMAX) + DBELLYLOWER(LMAX))
BELLYMAXRELA= BELLYMAX/DEPTHBELLY
DIMQBELLY= (QBELLY(LMAX) - QBASE)/(QPEAK - QBASE)
DIMQBELLY1= (QBELLY(LMAX) - QBASE)/(Q5RATMAX - QBASE)
WRITE(LU6,'(1X,A,F10.3)')
1'MAXIMUM DEPTH DIFFERENCE IN THE LOOP IS (M) = ',BELLYMAX
WRITE(LU8,'(1X,A,F10.3)')
1'MAXIMUM DEPTH DIFFERENCE IN THE LOOP IS (M) = ',BELLYMAX
WRITE(LU6,'(1X,A,F10.3)')
1'MAXIMUM RELATIVE DIFFERENCE IN THE LOOP IS = ',BELLYMAXRELA
WRITE(LU8,'(1X,A,F10.3)')
1'MAXIMUM RELATIVE DIFFERENCE IN THE LOOP IS = ',BELLYMAXRELA
WRITE(LU6,'(1X,A,F10.3)')
1'DIMENSIONLESS DISCHARGE AT THE MAXIMUM LOOP (OVER QPEAK) IS = ',
2DIMQBELLY
WRITE(LU8,'(1X,A,F10.3)')
1'DIMENSIONLESS DISCHARGE AT THE MAXIMUM LOOP (OVER QPEAK) IS = ',
2DIMQBELLY
WRITE(LU6,'(1X,A,F10.3)')
1'DIMENSIONLESS DISCHARGE AT THE MAXIMUM LOOP (OVER Q5RAT) IS = ',
2DIMQBELLY1
WRITE(LU8,'(1X,A,F10.3)')
1'DIMENSIONLESS DISCHARGE AT THE MAXIMUM LOOP (OVER Q5RAT) IS = ',
2DIMQBELLY1
WRITE(LU6,'(1X,79(1H-))')
WRITE(LU8,'(1X,79(1H-))')
C------------------------------------------------------------------------------
C SIX-POINT COMPUTATION (TO CORRECT MASS CONSERVATION)
C------------------------------------------------------------------------------
NDEN= 0.5*(NT/KTTRATIO + NCOUNT)
CORRECTION= (100. - (PERC5 - 100.)/NDEN)/100.
DO 90 J= 1,NR
DO 86 N= 1,NT
JFLAG= 0
JTRIAL= 0
QTRIAL= Q4(J,N)
DTRIAL= D4(J,N)
DO 84 WHILE (JFLAG.EQ.0)
JTRIAL= JTRIAL + 1
DEDIFF= DTRIAL - D6(J-1,N)
WSSL(J,N)= SLOPE0 - DEDIFF/DX
AWSSL(J,N) = ALWE*WSSL(J,N) + ALWER*WSSL(J,N-1)
SLOPE= AWSSL(J,N)
ALPHA= (CMAN/XN)*SQRT(SLOPE)
DREF=
10.5*(DEWER*(D6(J-1,N-1) + D6(J,N-1)) + DEWE*(D6(J-1,N) + DTRIAL))
QREF= ALPHA*DREF**BETA
CALL PARAME
CALL COEFF(C0,C1,C2)
Q6(J,N)= C0*Q6(J-1,N) + C1*Q6(J-1,N-1) + C2*Q6(J,N-1)
D6(J,N)= (Q6(J,N)/ALPHA)**REBETA
IF(ABS(Q6(J,N)-QTRIAL).LT.QTOL) THEN
JFLAG= 1
Q6(J,N)= CORRECTION*(Q6(J,N) - QBASE) + QBASE
D6(J,N)= (Q6(J,N)/ALPHA)**REBETA
ENDIF
QTRIAL= Q6(J,N)
DTRIAL= D6(J,N)
84 ENDDO
IF(MOD(J,JPRINT).EQ.0) THEN
IF(MOD(N,NPRINT).EQ.0) THEN
C WRITE(LU6,'(A,4I5)') ' STILL WORKING, J,N ',J,N
ENDIF
ENDIF
86 ENDDO
90 ENDDO
C------------------------------------------------------------------------------
C PRINTING RESULTS OF FOUR AND SIX-POINT METHODS
C------------------------------------------------------------------------------
WRITE(LU6,'(1X,79(1H-))')
WRITE(LU8,'(1X,79(1H-))')
WRITE(LU11,'(1X,79(1H-))')
WRITE(LU6,*)
1' STEP TIME INFLOW OUTFLOW4 DEPTH4 OUTFLOW6
2 DEPTH6'
WRITE(LU8,*)
1' STEP TIME INFLOW OUTFLOW4 DEPTH4 OUTFLOW6
2 DEPTH6'
WRITE(LU11,*)
1' STEP TIME INFLOW OUTFLOW4 DEPTH4 OUTFLOW6
2 DEPTH6'
WRITE(LU6,'(1X,79(1H-))')
WRITE(LU8,'(1X,79(1H-))')
WRITE(LU11,'(1X,79(1H-))')
DO 92 N= 0,NT
TIME(N) = TIME(N-1) + DTH
WRITE(LU6,'(I6,F8.0,5F12.3)')
1N,TIME(N),Q4(0,N),Q4(NR,N),D4(NR,N),Q6(NR,N),D6(NR,N)
WRITE(LU8,'(I6,F8.0,5F12.3)')
1N,TIME(N),Q4(0,N),Q4(NR,N),D4(NR,N),Q6(NR,N),D6(NR,N)
WRITE(LU11,'(I6,1H,,F7.0,5(1H,,F11.3))')
1N,TIME(N),Q4(0,N),Q4(NR,N),D4(NR,N),Q6(NR,N),D6(NR,N)
92 ENDDO
DO 94 N= 1,NT-1,2
SUMIN46= SUMIN43
SUMOUT46= SUMOUT46 + Q6(NR,N) - QBASE
94 ENDDO
DO 96 N= 2,NT-2,2
SUMIN26= SUMIN23
SUMOUT26= SUMOUT26 + Q6(NR,N) - QBASE
96 ENDDO
VOLIN6= VOLIN3
VOLOUT6=
1(DT/3)*(Q6(NR,0) - 2*QBASE + 4*SUMOUT46 + 2*SUMOUT26 + Q6(NR,NT))
PERC6= 100.*(VOLOUT6/VOLIN6)
WRITE(LU6,'(1X,79(1H-))')
WRITE(LU8,'(1X,79(1H-))')
WRITE(LU6,'(1X,A,F10.3)')
1'PERCENTAGE MASS CONSERVATION (6-POINT METHOD)= ',PERC6
WRITE(LU8,'(1X,A,F10.3)')
1'PERCENTAGE MASS CONSERVATION (6-POINT METHOD)= ',PERC6
C------------------------------------------------------------------------------
C COMPUTATION OF SIZE OF BELLY IN LOOPED RATING, 6-POINT METHOD
C------------------------------------------------------------------------------
NK= 0
DO 300 N= 1,NT
Q5TEMP= Q6(NR,N) - QBASE
IF(Q5TEMP.GE.0.001) THEN
NK= NK + 1
Q5RAT(NK)= Q6(NR,N)
D5RAT(NK)= D6(NR,N)
ENDIF
300 ENDDO
C WRITE(LU6,*) 'NK =',NK
C WRITE(LU8,*) 'NK =',NK
Q5RATMAX= QBASE
MFLAG= 0
DO 310 WHILE(MFLAG.EQ.0)
DO 305 L= 1,NK
IF(Q5RAT(L).GT.Q5RATMAX) THEN
Q5RATMAX= Q5RAT(L)
NPL= L
ELSE
MFLAG= 1
ENDIF
305 ENDDO
310 ENDDO
C WRITE(LU6,*) 'Q5RATMAX =',Q5RATMAXC
C WRITE(LU8,*) 'Q5RATMAX =',Q5RATMAX
NPU= NK - NPL + 1
DO 320 L= 1,NPL
Q5LOWER(L)= Q5RAT(L)
D5LOWER(L)= D5RAT(L)
320 ENDDO
DO 330 K= 1,NPU
L= NK - K + 1
Q5UPPER(K) = Q5RAT(L)
D5UPPER(K) = D5RAT(L)
330 ENDDO
C WRITE(LU6,*) NPL,(Q5LOWER(L),L=1,NPL)
C WRITE(LU6,*) NPU,(Q5UPPER(L),L=1,NPU)
C WRITE(LU6,*) NPL,(D5LOWER(L),L=1,NPL)
C WRITE(LU6,*) NPU,(D5UPPER(L),L=1,NPU)
C WRITE(LU8,*) NPL,(Q5LOWER(L),L=1,NPL)
C WRITE(LU8,*) NPU,(Q5UPPER(L),L=1,NPU)
C WRITE(LU8,*) NPL,(D5LOWER(L),L=1,NPL)
C WRITE(LU8,*) NPU,(D5UPPER(L),L=1,NPU)
QINTERVAL= 0.01*(Q5RATMAX - QBASE)
QBELLY(0)= QBASE
DO 340 L= 1,100
QBELLY(L)= QBELLY(L-1) + QINTERVAL
XL= QBELLY(L)
CALL PWLIN(Q5LOWER,D5LOWER,NPL,XL,YL)
DBELLYLOWER(L)= YL
CALL PWLIN(Q5UPPER,D5UPPER,NPU,XL,YL)
DBELLYUPPER(L)= YL
BELLY(L)= DBELLYUPPER(L) - DBELLYLOWER(L)
340 ENDDO
C WRITE(LU6,*) 'QBELLY(L)= ',(QBELLY(L),L=1,100)
C WRITE(LU8,*) 'QBELLY(L)= ',(QBELLY(L),L=1,100)
C WRITE(LU6,*) 'DBELLYUPPER(L)= ',(DBELLYUPPER(L),L=1,100)
C WRITE(LU8,*) 'DBELLYUPPER(L)= ',(DBELLYUPPER(L),L=1,100)
C WRITE(LU6,*) 'DBELLYLOWER(L)= ',(DBELLYLOWER(L),L=1,100)
C WRITE(LU8,*) 'DBELLYLOWER(L)= ',(DBELLYLOWER(L),L=1,100)
C WRITE(LU6,*) 'BELLY(L)= ',(BELLY(L),L=1,100)
C WRITE(LU8,*) 'BELLY(L)= ',(BELLY(L),L=1,100)
BELLYMAX= 0.
DO 350 L= 1,100
IF(BELLY(L).GT.BELLYMAX) THEN
BELLYMAX= BELLY(L)
LMAX= L
ENDIF
350 ENDDO
DEPTHBELLY= 0.5*(DBELLYUPPER(LMAX) + DBELLYLOWER(LMAX))
BELLYMAXRELA= BELLYMAX/DEPTHBELLY
DIMQBELLY= (QBELLY(LMAX) - QBASE)/(QPEAK - QBASE)
DIMQBELLY1= (QBELLY(LMAX) - QBASE)/(Q5RATMAX - QBASE)
WRITE(LU6,'(1X,A,F10.3)')
1'MAXIMUM DEPTH DIFFERENCE IN THE LOOP IS (M) = ',BELLYMAX
WRITE(LU8,'(1X,A,F10.3)')
1'MAXIMUM DEPTH DIFFERENCE IN THE LOOP IS (M) = ',BELLYMAX
WRITE(LU6,'(1X,A,F10.3)')
1'MAXIMUM RELATIVE DIFFERENCE IN THE LOOP IS = ',BELLYMAXRELA
WRITE(LU8,'(1X,A,F10.3)')
1'MAXIMUM RELATIVE DIFFERENCE IN THE LOOP IS = ',BELLYMAXRELA
WRITE(LU6,'(1X,A,F10.3)')
1'DIMENSIONLESS DISCHARGE AT THE MAXIMUM LOOP (OVER QPEAK) IS = ',
2DIMQBELLY
WRITE(LU8,'(1X,A,F10.3)')
1'DIMENSIONLESS DISCHARGE AT THE MAXIMUM LOOP (OVER QPEAK) IS = ',
2DIMQBELLY
WRITE(LU6,'(1X,A,F10.3)')
1'DIMENSIONLESS DISCHARGE AT THE MAXIMUM LOOP (OVER Q6RAT) IS = ',
2DIMQBELLY1
WRITE(LU8,'(1X,A,F10.3)')
1'DIMENSIONLESS DISCHARGE AT THE MAXIMUM LOOP (OVER Q6RAT) IS = ',
2DIMQBELLY1
WRITE(LU6,'(1X,79(1H-))')
WRITE(LU8,'(1X,79(1H-))')
WRITE(LU11,'(1X,79(1H-))')
C------------------------------------------------------------------------------
END
C------------------------------------------------------------------------------
C SUBROUTINE PARAME
C------------------------------------------------------------------------------
SUBROUTINE PARAME
COMMON /PARAIN/ QREF,DREF,SLOPE,ALPHA,BETA,REBETA,DT,DX,GR,BEONESQ
COMMON /PARAOUT/ COU,REY
DOUBLE PRECISION COU,REY
DOUBLE PRECISION QREF,DREF
DOUBLE PRECISION SLOPE,ALPHA,BETA,REBETA
DOUBLE PRECISION DT,DX,GR,BEONESQ
CEL= BETA*QREF*(QREF/ALPHA)**-REBETA
COU= CEL*DT/DX
FROUDESQ= QREF**2/(GR*DREF**3)
VEDERSQ= BEONESQ*FROUDESQ
REY= (QREF/(SLOPE*CEL*DX))*(1. - VEDERSQ)
RETURN
C------------------------------------------------------------------------------
END
C------------------------------------------------------------------------------
C SUBROUTINE COEFF
C------------------------------------------------------------------------------
SUBROUTINE COEFF(C0,C1,C2)
C------------------------------------------------------------------------------
COMMON /PARAOUT/ COU,REY
DOUBLE PRECISION COU,REY
DOUBLE PRECISION C0,C1,C2
DOUBLE PRECISION C3
C3= 1./(1. + COU + REY)
C0= C3*(-1. + COU + REY)
C1= C3*( 1. + COU - REY)
C2= C3*( 1. - COU + REY)
RETURN
C------------------------------------------------------------------------------
END
C------------------------------------------------------------------------------
C SUBROUTINE PWLIN(X,Y,N,XL,YL)
C------------------------------------------------------------------------------
SUBROUTINE PWLIN(X,Y,N,XL,YL)
DIMENSION X(N),Y(N)
IF(XL.LT.X(1)) THEN
YL= ((Y(2)-Y(1))*XL+(Y(1)*X(2)-Y(2)*X(1)))/(X(2)-X(1))
ELSEIF(XL.GE.X(N)) THEN
K= N-1
YL= ((Y(N)-Y(K))*XL+(Y(K)*X(N)-Y(N)*X(K)))/(X(N)-X(K))
ELSE
DO 1 J= 2,N
IF(XL.GE.X(J-1).AND.XL.LT.X(J)) M=J
1 CONTINUE
L= M-1
YL= ((Y(M)-Y(L))*XL+(Y(L)*X(M)-Y(M)*X(L)))/(X(M)-X(L))
ENDIF
RETURN
C------------------------------------------------------------------------------
END
C------------------------------------------------------------------------------