www.pudn.com > PREISSMANN.rar > G.FOR
C COMPUTATION OF UNSTEADY, FREE-SURFACE FLOWS BY PREISSMANN
C IMPLICIT SCHEME IN A TRAPEZOIDAL CHANNEL
C CONSTANT FLOW DEPTH ALONG THE CHANNEL IS SPECIFIED AS
C INITIAL CONDITION.
C TRANSIENT CONDITIONS ARE PRODUCED BY THE SUDDEN CLOSURE
C OF A DOWNSTREAM GATE.
C
C ************************* NOTATION ************************
C
C ALPHA = WEIGHTING COEFFICIENT
C AR = STATEMENT FUNCTION FOR FLOW AREA
C B0 = CHANNEL BOTTOM WIDTH
C C = CELERITY
C CENTR = MOMENT OF FLOW AREA
C CMN = MANNING'S COEFFICIENT
C CHL = CHANNEL LENGTH
C G = ACCELERATION OF GRAVITY
C HR = STATEMENT FUNCTION FOR HYDRAULIC RADIUS
C IPRINT = COUNTER FOR PRINTING RESULTS
C MAXITER = MAXIMUM NUMBER OF ITERATIONS
C NSEC = NUMBER OF CHANNEL SECTIONS
C Q0 = INITIAL STEADY STATE DISCHARGE
C S = CHANNEL LATERAL SLOPE
C S0 = CHANNEL BOTTOM SLOPE
C TLAST = TIME FOR TRANSIENT FLOW COMPUTATION
C TOL = TOLERANCE FOR INTERATIONS
C TOP = STATEMENT FUNCTION FOR WATER TOP WIDTH
C V = FLOW VELOCITY
C Y = FLOW DEPTH
C UNITS = ALPHANUMERIC VARIABLE FOR UNITS SYSTEM
C ************************************************************
C
CHARACTER*15 UNITS
DIMENSION Y(60),V(60),C1(60),C2(60),DF(60),EQN(124,125)
AR(D) = (B0+D*S)*D
HR(D) = (B0+D*S)*D/(B0+2*D*SQRT(1+S*S))
TOP(D) = B0+2*D*S
CENTR(D) = D*D*(B0/2+D*S/3)
DCENDY(D)= D*(B0+D*S)
READ(5,*) G,NSEC,TLAST,IPRINT
READ(5,*) CHL,B0,S,CMN,S0,Q0,Y0,ALPHA,TOL,MAXITER
T = 0
IF(G.GT.10) THEN
CMN2 = (CMN*CMN)/2.22
ELSE
CMN2 = CMN*CMN
END IF
C
C STEADY STATE CONDITIONS
C
C = SQRT(G*AR(Y0)/TOP(Y0))
V0 = Q0/AR(Y0)
DX = CHL/NSEC
DT = DX/(V0+C)
DTX2 = 2*DT/DX
YRES = Y0
I = 1
NP1 = NSEC + 1
DO WHILE(I.LE.NP1)
V(I) = V0
Y(I) = Y0
I = I+1
END DO
IF(G.LE.10) UNITS = 'SI'
IF(G.GT.10) UNITS = 'ENGLISH'
WRITE(6,15) UNITS,CHL,B0,S,S0,CMN,Q0,Y0,NP1,TLAST
IFLAG = 0
IP = IPRINT
C
C COMPUTE TRANSIENT CONDITIONS
C
DO WHILE(T.LE.TLAST .AND. IFLAG.EQ.0)
ITER = 0
IF(IPRINT.EQ.IP) THEN
IP = 0
WRITE(6,20) T
WRITE(6,30) (Y(I),I= 1,NP1)
WRITE(6,40) (V(I),I= 1,NP1)
END IF
T = T+DT
C
C GENERATE SYSTEM OF EQUATIONS
C
I = 1
DO WHILE(I.LE.NSEC)
ARI = AR(Y(I))
ARIP1 = AR(Y(I+1))
C1(I)=DTX2*(1-ALPHA)*(ARIP1*V(I+1)-ARI*V(I))-ARI-ARIP1
SF1 = ABS(V(I))*V(I)*CMN2/HR(Y(I))**1.333
SF2 = ABS(V(I+1))*V(I+1)*CMN2/HR(Y(I+1))**1.333
TERM1 = -DT*(1-ALPHA)*(G*ARIP1*(S0-SF2)+G*ARI*(S0-SF1))
TERM2 = -(V(I)*ARI+V(I+1)*ARIP1)
TERM3 = DTX2*(1-ALPHA)*(V(I+1)**2*ARIP1+G*CENTR(Y(I+1))-
1 V(I)**2*ARI-G*CENTR(Y(I)))
C2(I) = TERM1 + TERM2 + TERM3
I = I+1
END DO
NP11 = 2*NP1 + 1
SUM = TOL+10
100 IF(.NOT.(SUM.LE.TOL))THEN
I = 1
DO WHILE(I.LE.2*NP1)
J = 1
DO WHILE(J.LE.NP11)
EQN(I,J) = 0.0
J = J+1
END DO
I = I+1
END DO
ITER = ITER+1
C
C BOUNDARY EQUATIONS
C
EQN(1,1) = 1.0
EQN(1,NP11) = -(Y(1)-YRES)
EQN(2*NP1,2*NP1) = 1.0
EQN(2*NP1,NP11) = -V(NP1)
C
C INTERIOR NODES
C
I = 1
DO WHILE(I.LE.NSEC)
ARI = AR(Y(I))
ARIP1 = AR(Y(I+1))
K = 2*I
EQN(K,NP11)=-(ARI+ARIP1+DTX2*ALPHA*(V(I+1)*ARIP1-V(I)*ARI)+C1(I))
SF1 = ABS(V(I))*V(I)*CMN2/HR(Y(I))**1.333
SF2 = ABS(V(I+1))*V(I+1)*CMN2/HR(Y(I+1))**1.333
TERM1 = DTX2*ALPHA*(V(I+1)**2*ARIP1+G*CENTR(Y(I+1))-V(I)**2*ARI
1 -G*CENTR(Y(I)))
TERM2 = -ALPHA*DT*G*((S0-SF2)*ARIP1+(S0-SF1)*ARI)
EQN(K+1,NP11) = -(V(I)*ARI+V(I+1)*ARIP1+TERM1+TERM2+C2(I))
DAY1 = TOP(Y(I))
DAY2 = TOP(Y(I+1))
EQN(K,K-1) = DAY1*(1-DTX2*ALPHA*V(I))
EQN(K,K) = -DTX2*ALPHA*ARI
EQN(K,K+1) = DAY2*(1+DTX2*ALPHA*V(I+1))
EQN(K,K+2) = DTX2*ALPHA*ARIP1
DCDY1 = DCENDY(Y(I))
DCDY2 = DCENDY(Y(I+1))
DSDV1 = 2*V(I)*CMN2/HR(I)**1.333
DSDV2 = 2*V(I+1)*CMN2/HR(I+1)**1.333
PERI = ARI/HR(Y(I))
TERM1 = 2*SQRT(1+S**2)*ARI - DAY1*PERI
TERM2 = HR(Y(I))**0.333*ARI**2
DSDY1 = 1.333*V(I)*ABS(V(I))*CMN2*TERM1/TERM2
PERIP1 = ARIP1/HR(Y(I+1))
TERM1 = 2*SQRT(1+S**2)*ARIP1-DAY2*PERIP1
TERM2 = HR(Y(I+1))**0.333*ARIP1**2
DSDY2 = 1.333*V(I+1)*ABS(V(I+1))*CMN2*TERM1/TERM2
TERM1 = -DTX2*ALPHA*(V(I)**2*DAY1 + G*DCDY1)
TERM2 = -G*DT*ALPHA*(S0-SF1)*DAY1
EQN(K+1,K-1)=V(I)*DAY1+TERM1+TERM2+G*DT*ALPHA*ARI*DSDY1
EQN(K+1,K)=ARI-DTX2*ALPHA*2*V(I)*ARI+G*DT*ALPHA*ARI*DSDV1
TERM1 = DTX2*ALPHA*(V(I+1)**2*DAY2+G*DCDY2)
TERM2 = -G*DT*G*(S0-SF2)*DAY2
EQN(K+1,K+1) = V(I+1)*DAY2+TERM1+TERM2+ALPHA*DT*G*ARIP1*DSDY2
EQN(K+1,K+2) = ARIP1+DTX2*ALPHA*2*V(I+1)*ARIP1+ALPHA*DT*G
1 *DSDV2*ARIP1
I = I+1
END DO
C
C SOLVE SYSTEM OF EQUATIONS
C
CALL MATSOL(2*NP1,DF,EQN)
I = 1
SUM = 0.0
DO WHILE(I.LE.2*NP1)
SUM = ABS(DF(I))+SUM
IF(MOD(I,2).EQ.1) Y(I/2+1) = Y(I/2+1)+DF(I)
IF(MOD(I,2).EQ.0) V(I/2) = V(I/2) + DF(I)
I = I+1
END DO
C
C CHECK NUMBER OF ITERATIONS
C
IF(ITER.GT.MAXITER) THEN
IFLAG = 1
SUM = TOL
END IF
GO TO 100
END IF
IP = IP+1
END DO
IF(IFLAG.EQ.1) WRITE(6,50)
15 FORMAT(5X,'PREISSMANN SCHEME',
1 /5X,'UNITS = ',A10,/5X,'CHANNEL LENGTH = ',F7.2,
2 /5X,'CHANNEL BOTTOM WIDTH = ',F6.2,/5X,'CHANNEL LATERAL SLOPE =
3 ',F5.2,/5X,'CHANNEL BOTTOM SLOPE = ',F10.5,/5X,'MANNING"S n = '
4,F10.5,/5X,'INITIAL STEADY STATE DISCHARGE = ',F7.2,/5X,
5'UNIFORM FLOW DEPTH = ',F7.2,/5X,'NUMBER OF CHANNEL SECTIONS = ',
6 I3,/5X,'TIME FOR WHICH TRANSIENTS WILL BE COMPUTED = ',F6.2)
20 FORMAT(/2X,'T = ',F8.3)
30 FORMAT(2X,'Y = ',(12F6.2))
40 FORMAT(2X,'V = ',(12F6.3))
50 FORMAT(5X,'MAXIMUM NUMBER OF ITERATIONS EXCEEDED')
STOP
END
C ******************************************************************
C SIMULTANEOUS SOLUTION OF THE SYSTEM OF EQUATIONS
C
SUBROUTINE MATSOL(N,X,A)
DIMENSION X(60),A(124,125),NROW(60)
I = 1
DO WHILE(I.LE.N)
NROW(I) = I
I = I + 1
END DO
I = 1
DO WHILE(I.LE.N-1)
AMAX = A(NROW(I),I)
J = I
IP =I
DO WHILE(J.LE.N)
IF(ABS(A(NROW(J),I)) .GT. AMAX) THEN
AMAX = ABS(A(NROW(J),I))
IP = J
END IF
J = J+1
END DO
IF(ABS(AMAX) .LE. 1E-08) GO TO 100
IF(NROW(I).NE.NROW(IP)) THEN
NC = NROW(I)
NROW(I) = NROW(IP)
NROW(IP) = NC
END IF
J = I+1
DO WHILE(J.LE.N)
COEF = A(NROW(J),I)/A(NROW(I),I)
JJ = I+1
DO WHILE(JJ.LE.N+1)
A(NROW(J),JJ)=A(NROW(J),JJ)-COEF*A(NROW(I),JJ)
JJ = JJ + 1
END DO
J = J+1
END DO
I = I+1
END DO
IF(ABS(A(NROW(N),N)) .LE. 1E-08) GO TO 100
X(N) = A(NROW(N),N+1)/A(NROW(N),N)
I = N-1
DO WHILE(I.GE.1)
SUM = 0.0
J = I+1
DO WHILE(J.LE.N)
SUM = A(NROW(I),J)*X(J) + SUM
J = J+1
END DO
X(I) = (A(NROW(I),N+1) - SUM)/A(NROW(I),I)
I = I-1
END DO
GO TO 200
100 WRITE(*,*) 'SINGULAR MATRIX --> NO UNIQUE SOLUTION EXISTS'
200 RETURN
END