PROGRAM EXPLIC C******************************************************************************* C * C MAIN PROGRAM FOR THE SOLUTION OF UNSTEADY FLOWS IN A CHANNEL OF TRAPEZOIDAL, * C RECTANGULAR OR TRAINGULAR CROSS SECTION * C * C PREPARED BY PROF. MUSTAFA ALTINAKAR, UNIVERSITY OF MISSISSIPPI * C MODIFIED BY PROF. FABIAN BOMBARDELLI, UC DAVIS * C * C JANUARY, 2009 * C******************************************************************************* C C LIST OF SIMPLE VARIABLES, OF VECTORS AND MATRICES FOR THE MAIN PROGRAM C AND SUB-ROUTINES C C TYPE NAME DIM. EXPLANATION C ==== ===== ======= ========================================================= C R*4 B = CHANNEL WIDTH C R*4 BETA = SEE EQN. IN THE TEXT OF THE PROJECT C R*4 C0 = CELERITY OF GRAVITY WAVES AT TIME T = 0 C R*4 DT = TIME STEP C R*4 DTMAX = TIME STEP ACCORDING TO THE COURANT STABILITY CRITERION C R*4 DX = SPATIAL STEP C CHAR FICH = NAME OF THE OUTPUT FILE C R*4 G = ACCELERATION OF GRAVITY (G = 9.81 m/s2) C R*4 GAMA = SEE EQN. IN THE TEXT OF THE PROJECT C R*4 H (2,101) = WATER DEPTH IN THE NODES FOR THE PREVIOUS TIME STEP C R*4 HN = NORMAL DEPTH C I*4 I = COUNTER OF NUMBER OF NODES C I*4 IL = COUNTER OF NUMBER OF LINES OF TITLE C R*4 JF = BOTTOM SLOPE OF THE CHANNEL C I*4 K = COUNTER OF LOOPS C R*4 LT = TOTAL LENGTH OF CHANNEL C R*4 M = SLOPE OF SIDES OF CROSS SECTION C I*4 MAXDIV = MAX. NUMBER OF SEGMENTS ALLOWED BY THE PROGRAM (100) C R*4 N = MANNING COEFFICIENT C I*4 NDIV = NUMBER OF SEGMENT C I*4 NN = NUMBER OF NODES C I*4 NNMAX = MAX. NUMBER OF NODES ALLOWED BY THE PROGRAM (100+1) C I*4 NOST (3) = NUMBER OF NODES WITH OUTPUT OF RESULTS C I*4 NP = NUMBER OF TIME STEPS AT TIME T C I*4 NPSTEP = FREQUENCY OF OUTPUT (NUMBER OF STEPS) OF OUTPUT OF RESULTS C I*4 NUNIT = NUMBER OF THE UNITY OF OUTPUT C R*4 P = WETTED PERIMETER C R*4 Q (2,101) = DISCHARGE AT NODES AT THE PREVIOUS TIME STEP C R*4 QE = INCOMING DISCHARGE AT TIME T C R*4 Q0 = INITIAL DISCHARGE (UNIFORM FLOW) C R*4 QMAX = MAX. DISCHARGE OF THE FLOOD C CHAR REP = ALPHANUMERIC RESPONSE C R*4 RH = HYDRAULIC RADIUS C R*4 S = WETTED SURFACE C R*4 T = DURATION OF COMPUTATIONS C CHAR TITRE (24) = TITLE OF THE OUTPUT FILE C R*4 TMAX = TIME LIMIT FOR THE COMPUTATIONS C R*4 TP = DURATION OF THE PEAK OF THE FLOOD C R*4 TPD = DURATION OF THE RECESS OF THE FLOOD C R*4 TTH = TOTAL DURATION OF THE FLOOD (= TP + TPD) C R*4 U (2,101) = VELOCITY AT THE NODES AT THE PREVIOUS TIME STEP C R*4 U0 = INITIAL VELOCITY (UNIFORM FLOW) C C******************************************************************************* C C PARAMETERS PARAMETER ( NUNIT = 5 ) PARAMETER ( MAXDIV = 100 , NNMAX = MAXDIV+1 ) PARAMETER ( G = 9.81 ) C C DECLARATION OF VARIABLES CHARACTER*50 FICH REAL M , JF , N , LT DIMENSION H(2,NNMAX) , U(2,NNMAX) , Q(2,NNMAX) , NOST(3) C C CALLS SUB-ROUTINE "LIRE" TO READ THE DATA OF THE PROBLEM C AND OTHER NEEDED PARAMETERS C CALL LIRE( MAXDIV , NNMAX , G , 1 B , M , JF , N , HN , LT , Q0 , U0 , C0 , 2 TP , QMAX , TPD , 3 NDIV , NN , DX , DT , TMAX , 4 NOST , NPSTEP , NUNIT , FICH ) C C CALLS SUB-ROUTINE "INIT" TO CALCULATE THE INITIAL CONDITIONS AT ALL C NODES AT TIME T=0 C CALL INIT ( NNMAX , NN , H , U , Q , HN , U0 , Q0 ) C C CALLS SUB-ROUTINE "CALCUL" TO CALCULATE THE UNSTEADY FLOW C CALL CALCUL ( NNMAX , NN , H , U , Q , HN , LT , 1 B , M , JF , N , Q0 , U0 , C0 , G , 2 TP , QMAX , TPD , DX , DT , TMAX , 3 NOST , NPSTEP , NUNIT ) C C END OF PROGRAM C STOP END CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC SUBROUTINE LIRE( MAXDIV , NNMAX , G , 1 B , M , JF , N , HN , LT , Q0 , U0 , C0 , 2 TP , QMAX , TPD , 3 NDIV , NN , DX , DT , TMAX , 4 NOST , NPSTEP , NUNIT , FICH ) C******************************************************************************* C READS THE DATA FOR THE PROGRAM * C * C******************************************************************************* C C DECLARATION OF VARIABLES CHARACTER*1 REP CHARACTER*50 FICH REAL M , JF , N , LT DIMENSION NOST(3) C 5 FORMAT(A) OPEN( UNIT = 7 , FILE = 'DIALOG.DAT' , STATUS = 'NEW') C C PROGRAMM TITLE C WRITE(*,*) WRITE(7,*) WRITE(*,*) WRITE(7,*) WRITE(*,10) WRITE(7,10) 10 FORMAT(//' PROGRAM TO CALCULATE UNSTEADY FLOWS IN CHANNELS' 1 //' SYSTEM OF UNITS : METRIC'//) C C READING DATA REGARDING THE CHANNEL C WRITE(*,20) WRITE(7,20) 20 FORMAT(' DATA CONCERNING THE CHANNEL :'/ 1 ' THE CHANNEL MAY BE TRAPEZOIDAL, RECTANGULAR' 2,'(m = 0) OR'/ 3' TRIANGULAR (B = 0).'//) 30 WRITE(*,40) WRITE(7,40) 40 FORMAT(' WIDTH OF THE CHANNEL ? B (m) = ',$) READ(*,*,ERR=30)B 50 WRITE(*,60) WRITE(7,60) 60 FORMAT(' SLOPE OF SIDES OF CROSS SECTION ? m (-) = ',$) READ(*,*,ERR=50)M 70 WRITE(*,80) WRITE(7,80) 80 FORMAT(' CHANNEL SLOPE ? Jf (-) = ',$) READ(*,*,ERR=70)JF 90 WRITE(*,100) WRITE(7,100) 100 FORMAT(' MANNING COEFFICIENT ? n (m^-1/3 s) = ',$) READ(*,*,ERR=90)N 110 WRITE(*,120) WRITE(7,120) 120 FORMAT(' UNIFORM FLOW DEPTH ? hn (m) = ',$) READ(*,*,ERR=110)HN 130 WRITE(*,140) WRITE(7,140) 140 FORMAT(' TOTAL LENGTH OF CHANNEL ? LT (m) = ',$) READ(*,*,ERR=130)LT C C CALCULATE THE DISCHARGE FOR THE UNIFORM FLOW C S = HN * (B + M * HN) P = B + 2 * HN * SQRT(1 + M**2) RH = S / P U0 = RH**(2./3.) * SQRT(JF) / N Q0 = U0 * S WRITE(*,150) Q0 WRITE(7,150) Q0 150 FORMAT(//' UTILIZING MANNINGS FORMULA '/ 1 ' THE DISCHARGE FOR THE UNIFORM FLOW IS, Qo (m3/s) = ', 2 F8.3) C C READING OF THE DATA CONCERNING THE CONDITIONS AT BOUNDARIES C WRITE(*,200) WRITE(7,200) 200 FORMAT(///' READING OF THE DATA CONCERNING THE CONDITIONS', 1 ' AT BOUNDARIES :'/ 2 ' IT ADMITS A TRIANGULAR HYDROGRAPH', 3 ' DEFINED BY THREE PARAMETERS :'/) 210 WRITE(*,220) WRITE(7,220) 220 FORMAT(' TIME TO THE PEAK ? t'' (s) = ',$) READ(*,*,ERR=210)TP 230 WRITE(*,240) WRITE(7,240) 240 FORMAT(' PEAK DISCHARGE ? QMAX (m3/s) = ',$) READ(*,*,ERR=230)QMAX 250 WRITE(*,260) WRITE(7,260) 260 FORMAT(' TIME OF DECAY ? t'''' (s) = ',$) READ(*,*,ERR=250)TPD TTH = TP + TPD WRITE(*,270) TTH WRITE(7,270) TTH 270 FORMAT(' ATTENTION ! ...'/' FOR t > (t'' + t'''') = ',F10.3, 1 ' (s), IT IS ASSUMED THAT Q = Qo'///) C C READING OF SPATIAL AND TIME DATA C NDIV = MAXDIV NN = NNMAX DX = LT / MAXDIV 300 WRITE(*,310) LT , MAXDIV , DX WRITE(7,310) LT , MAXDIV , DX 310 FORMAT(/' READING OF DATA CONCERNING THE SPATIAL STEP, DX :'/ 1 ' THE CHANNEL LENGTH IS ', 2 ' LT (m) = ',F7.2/ 3 ' MAXIMUM NUMBER OF STEPS ALLOWED BY THE PROGRAM, ', 4 ' MAXDIV = ',I3/ 5 ' THE MINIMUM VALUE OF DX IS, ', 6 ' LT / MAXDIV = ',F7.2, '(m)'// 7 ' (Important note! In you want to have more steps', 8 ' modify'/' the value of the parameter MAXDIV', 9 ' on the main program)' A //' DO YOU AGREE WITH THIS VALUE ? ANSWER', B ' O/CR OR N = ',$) READ(*,5)REP IF(REP.EQ.'N'.OR.REP.EQ.'n')THEN 320 WRITE(*,330) WRITE(7,330) 330 FORMAT(' WHAT IS THE NUMBER OF STEPS, NDIV ? = '$) READ(*,*,ERR=320)NDIV IF(NDIV.LE.0.OR.NDIV.GT.MAXDIV) GO TO 300 DX = LT / NDIV NN = NDIV + 1 WRITE(*,340) DX , NN WRITE(7,340) DX , NN 340 FORMAT(/' THE VALUE OF DX IS, THEN, LT', 1 ' / NDIV = ',F7.2, '(m)'/ 2 ' THE NUMBER OF NODES IS, NN =', 3 ' NDIV + 1 = ',I3) ENDIF C C0 = SQRT(G * HN) DTMAX = DX / (U0 + C0) DT = DTMAX 400 WRITE(*,410) Q0 , U0 , C0 , DTMAX WRITE(7,410) Q0 , U0 , C0 , DTMAX 410 FORMAT(//' READING OF DATA CONCERNING THE TIME STEP, DT :'// 1 ' THE DISCHARGE OF UNIFORM FLOW IS, Qo', 2 ' (m3/s) = ',F8.3/ 3 ' THE FLOW AVERAGED VELOCITY VALUE IS, Uo (m/s) = Q', 4 ' / S = ',F8.3/ 5 ' THE CELERITY OF A GRAVITY WAVE IS, Co (m/s) =', 6 ' (g*hn)^0.5 = ',F8.3// 7 ' THE COURANT STABILITY CONDITION SAYS THAT DT < (DX', 8 ' / (ABS(Uo) + Co)'/ 9 ' THE MAXIMUM VALUE OF DT IS: DTMAX= DX / (ABS(Uo)', A ' + Co) : ',F8.3,' s'//) 420 WRITE(*,430) WRITE(7,430) 430 FORMAT(' WHICH IS THE VALUE OF DT (s) ? = '$) READ(*,*,ERR=420)DT IF(DT.LE.0.0.OR.DT.GT.DTMAX) GO TO 400 C C WHEN WILL THE COMPUTATIONS FINISH ? 440 WRITE(*,450) WRITE(7,450) 450 FORMAT(//' INPUT THE COMPUTATION TIME, TMAX (s) ? = ',$) READ(*,*,ERR=440)TMAX IF(TMAX.LT.DT) GO TO 440 C C HOW TO PRINT THE RESULTS ? 500 WRITE(*,510)NN WRITE(7,510)NN 510 FORMAT(//' YOU HAVE ',I3,' STATIONS ALONG THE CHANNEL;'/ 1 ' WHERE DO YOU WANT TO HAVE THE RESULTS ?'/ 2 ' WRITE THE NUMBERS OF 3 STATIONS IN FREE FOMAT = ',$) READ(*,*,ERR=500) ( NOST(K) , K= 1,3) DO 520 K=1,3 IF(NOST(K).LT.1 .OR. NOST(K).GT.NN)GO TO 500 520 CONTINUE C C FREQUENCY OF WRITING OF RESULTS 550 WRITE(*,560) WRITE(7,560) 560 FORMAT(//' AT WHICH FREQUENCY DO YOU WANT TO WRITE YOUR', 1 ' RESULTS ? = ',$) READ(*,*,ERR=550)NPSTEP IF(NPSTEP.LT.1 .OR. NPSTEP.GT.(TMAX/DT))GO TO 550 C C OUTPUT FILE 600 WRITE(*,610) WRITE(7,610) 610 FORMAT(//' NAME OF THE OUTPUT FILE ? = ',$) READ(*,5) FICH OPEN( UNIT = NUNIT , FILE = FICH , STATUS = 'NEW' , ERR = 600 ) C RETURN END CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC SUBROUTINE INIT ( NNMAX , NN , H , U , Q , HN , U0 , Q0 ) C******************************************************************************* C C SUB-ROUTINE FOR THE SETTING OF INITIAL CONDITIONS (T=0) AT NODES * C * C******************************************************************************* C C DECLARATION OF VARIABLES DIMENSION H(2,NNMAX) , U(2,NNMAX) , Q(2,NNMAX) C DO 10 I = 1,NN H(1,I) = HN U(1,I) = U0 Q(1,I) = Q0 10 CONTINUE C RETURN END CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC SUBROUTINE CALCUL ( NNMAX , NN , H , U , Q , HN , LT , 1 B , M , JF , N , Q0 , U0 , C0 , G , 2 TP , QMAX , TPD , DX , DT , TMAX , 3 NOST , NPSTEP , NUNIT ) C******************************************************************************* C * C SUBROUTINE FOR THE CALCULATION (PER TIME STEP) OF THE DEPTH, THE VELOCITY * C AND THE DISCHARGE AT EVERY NODE FOR AN UNSTEADY FLOW ACCORDING TO THE * C EXPLICIT METHOD (SEE THEORY) C * C******************************************************************************* C C DECLARATION OF VARIABLES REAL M , N , JF , LT DIMENSION H(2,NNMAX) , U(2,NNMAX) , Q(2,NNMAX) , NOST(3) C C INITIALIZATION OF TIMES T = 0 C C WRITE TITLES IN THE OUTPUT CALL TITRES ( B , N , M , HN , JF , LT , Q0 , U0 , C0 , 1 TP , QMAX , TPD , NDIV , DX , NN , DT , 2 TMAX , NPSTEP , NOST , NUNIT ) C C PRINT RESULTS AT TIME T 100 CALL ECRIT ( NNMAX , NN , H , U , Q , T , DT , 1 NOST , NPSTEP , NUNIT ) C C INCREASE TIME IF IT IS APPROPRIATE IF(T+DT.GT.TMAX) GO TO 900 T = T + DT C C COMPUTE THE INCOMING DISCHARGE AT THE NEW TIME T CALL QENTRE ( TP , QMAX , TPD , T , Q0 , QE ) C C CALCULATE THE NODES FOR THE NEW TIME CALL AMONT( NNMAX , NN , H , U , Q , 1 B , M , QE , DX , DT ) C C COMPUTE THE INTERMEDIATE POINTS AT THE NEW TIME T CALL INTER( NNMAX , NN , H , U , Q , 1 B , M , JF , N , G , DX , DT , NUNIT ) C C CALCULER LE NOEUD EN AVAL DU CANAL AU NOUVEAU TEMPS T CALL AVAL( NNMAX , NN , H , U , Q , 1 B , M , JF , N , DX , DT ) C C REPLACE THE VALUES OF THE PREVIOUS TIME STEP BY THE NEW ONES DO 200 I = 1, NN H(1,I) = H(2,I) U(1,I) = U(2,I) Q(1,I) = Q(2,I) 200 CONTINUE C C ALLOW INCREMENTS OF TIME AND REPEAT THE COMPUTATIONS GO TO 100 C 900 WRITE(*,910) WRITE(NUNIT,910) 910 FORMAT(//' NORMAL TERMINATION OF PROGRAM') RETURN END CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC SUBROUTINE TITRES ( B , N , M , HN , JF , LT , Q0 , U0 , C0 , 1 TP , QMAX , TPD , NDIV , DX , NN , DT , 2 TMAX , NPSTEP , NOST , NUNIT ) C******************************************************************************* C * C THIS SUB-ROUTINE WRITES THE TITLE AND THE DATA OF THE PROBLEM IN THE OUTPUT * C FILE C * C******************************************************************************* C C DECLARATION OF VARIABLES CHARACTER*131 TITRE(24) REAL M , JF , N , LT DIMENSION NOST(3) C DATA TITRE(1)/' RESULTS OF COMPUTATIONS 1 FOR THE UNSTEADY FLOW WITH THE EXPLICIT METHOD'/ DATA TITRE(2)/' ====================== 1================================================'/ DATA TITRE(3),TITRE(4)/' ',' '/ DATA TITRE(5)/' DATA CONCE 1RNING THE CHANNEL AND THE UNIFORM FLOW'/ DATA TITRE(6)/' ------------- 1---------------------------------------'/ DATA TITRE(7)/' CHANNEL WIDTH, B (m) = 1 MANNING COEFFICIENT, n = '/ DATA TITRE(8)/' SLOPES OF SIDES, m (-) = 1 DEPTH FOR UNIFORM FLOW, hn (m) = '/ DATA TITRE(9)/' SLOPE OF THE CHANNEL BOTTOM, Jf (-) = 1 TOTAL CHANNEL LENGTH, LT (m) = '/ DATA TITRE(10)/' '/ DATA TITRE(11)/' INITIAL DISCHARGE, Q0 (M3/S) = XXXX 1.XXX INITIAL VELOCITY, U0 (m/s) = '/ DATA TITRE(12)/' CELERITY OF 1 GRAVITY WAVE C0 (m/s) = '/ DATA TITRE(13)/' '/ DATA TITRE(14)/' 1 TIME TO PEAK, t'' (s) = '/ DATA TITRE(15)/' DEFINITION OF THE TRIANGULAR HYDROGRAPH : 1 PEAK DISCHARGE, QMAX (m3/s) = '/ DATA TITRE(16)/' 1 TIME OF DECAY, t'''' (s) = '/ DATA TITRE(17)/' '/ DATA TITRE(18)/' NUMBER OF STEPS, NDIV = LENGTH OF 1 STEPS, DX (m) = NUMBER OF NODES, NN=NDIV+1 2 = '/ DATA TITRE(19)/' TIME STEP, DT (s) = DURATION OF 1 COMPUTATIONS, TMAX (s) = FREQUENCY OF OUTPUT (STEPS) 2 = '/ DATA TITRE(20)/' '/ DATA TITRE(21)/' STATION 1 1 STATION 2 STATION 3'/ DATA TITRE(22)/' TIME NODE NO = X = 1 NODE NO = X = NODE NO = X = 2'/ DATA TITRE(23)/' T Q H U 1 Q H U Q H 2 U'/ DATA TITRE(24)/' (s) (m3/s) (m) (m/s) 1 (m3/s) (m) (m/s) (m3/s) (m) 2(m/s)'/ C WRITE(TITRE(7)(45:51),'(F7.3)')B WRITE(TITRE(7)(109:114),'(F6.4)')N WRITE(TITRE(8)(45:51),'(F7.3)')M WRITE(TITRE(8)(109:114),'(F6.3)')HN WRITE(TITRE(9)(45:52),'(F8.5)')JF WRITE(TITRE(9)(109:116),'(F8.3)')LT WRITE(TITRE(11)(47:54),'(F8.3)')Q0 WRITE(TITRE(11)(101:108),'(F8.3)')U0 WRITE(TITRE(12)(78:85),'(F8.3)')C0 WRITE(TITRE(14)(85:93),'(F9.2)')TP WRITE(TITRE(15)(85:92),'(F8.3)')QMAX WRITE(TITRE(16)(85:93),'(F9.2)')TPD WRITE(TITRE(18)(29:31),'(I3)')NN-1 WRITE(TITRE(18)(73:80),'(F8.3)')DX WRITE(TITRE(18)(120:122),'(I3)')NN WRITE(TITRE(19)(29:35),'(F7.3)')DT WRITE(TITRE(19)(73:83),'(F11.3)')TMAX WRITE(TITRE(19)(120:124),'(I5)')NPSTEP C WRITE(TITRE(22)(31:33),'(I3)')NOST(1) WRITE(TITRE(22)(43:50),'(F8.3)')DX*(NOST(1)-1) WRITE(TITRE(22)(68:70),'(I3)')NOST(2) WRITE(TITRE(22)(80:87),'(F8.3)')DX*(NOST(2)-1) WRITE(TITRE(22)(105:107),'(I3)')NOST(3) WRITE(TITRE(22)(117:124),'(F8.3)')DX*(NOST(3)-1) C DO 200 IL = 1 , 24 WRITE(NUNIT,100)TITRE(IL) 100 FORMAT(A) 200 CONTINUE RETURN END CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC SUBROUTINE ECRIT ( NNMAX , NN , H , U , Q , T , DT , 1 NOST , NPSTEP , NUNIT ) C******************************************************************************* C * C SUB-ROUTINE TO OUTPUT THE VALUES OF H , U , AND Q AT TIME T. * C * C******************************************************************************* C C DECLARATION OF VARIABLES DIMENSION H(2,NNMAX) , U(2,NNMAX) , Q(2,NNMAX) , NOST(3) C NP = T / DT IF( MOD(NP,NPSTEP).NE.0 ) GO TO 100 WRITE(NUNIT,10) T , Q(1,NOST(1)) , H(1,NOST(1)) , U(1,NOST(1)) , 1 Q(1,NOST(2)) , H(1,NOST(2)) , U(1,NOST(2)) , 2 Q(1,NOST(3)) , H(1,NOST(3)) , U(1,NOST(3)) 10 FORMAT(3X,F10.2,3(5X,F10.3,1X,F10.3,1X,F10.3)) C 100 RETURN END CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC SUBROUTINE QENTRE ( TP , QMAX , TPD , T , Q0 , QE ) C******************************************************************************* C * C SUB-ROUTINE TO CALCULATE THE DISCHARGE AT THE UPSTREAM BOUNDARY AT THE NEW * C TIME T USING LINEAR INTERPOLATION. THE SHAPE OF THE HYDROGRAPH IS TRIANGULAR * C * C******************************************************************************* IF(T.LE.TP)THEN C C PARTIE MONTANTE DE L'HYDROGRAMME QE = Q0 + T * (QMAX - Q0) / TP RETURN ENDIF C IF(T.LT.TP+TPD)THEN C C PARTIE DESCENDANTE DE L'HYDROGRAMME QE = QMAX - (T - TP) * (QMAX - Q0) / TPD RETURN ENDIF C C ON A T >= (TP + TPD), L'HYDROGRAMME EST TERMINE. QE = Q0 C RETURN END CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC SUBROUTINE AMONT( NNMAX , NN , H , U , Q , 1 B , M , QE , DX , DT ) C******************************************************************************* C * C SUB-ROUTINE TO CALCULATE THE DEPTH, THE VELOCITY AND THE DISCHARGE * C AT NODE 1, FOR THE NEW TIME C * C******************************************************************************* C C DECLARATION OF VARIABLES REAL M DIMENSION H(2,NNMAX) , U(2,NNMAX) , Q(2,NNMAX) C H(2,1) = H(1,1) + (DT/DX) * 1 ( U(1,1)*(H(1,1)-H(1,2)) + H(1,1)*(U(1,1)-U(1,2)) ) C S = H(2,1) * (B + M * H(2,1)) C Q(2,1) = QE C U(2,1) = Q(2,1) / S C RETURN END CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC SUBROUTINE INTER( NNMAX , NN , H , U , Q , 1 B , M , JF , N , G , DX , DT , NUNIT ) C******************************************************************************* C * C SUB-ROUTINE TO CALCULATE THE DEPTH, THE VELOCITY AND THE DISCHARGE * C AT INTERMADIATE NODES, FOR THE NEW TIME * C * C******************************************************************************* C C DECLARATION OF VARIABLES REAL M , N , JF DIMENSION H(2,NNMAX) , U(2,NNMAX) , Q(2,NNMAX) C DO 300 I = 2 , NN-1 C H(2,I) = H(1,I) + 0.5 * (DT/DX) * 1 ( U(1,I)*(H(1,I-1)-H(1,I+1)) + 2 H(1,I)*(U(1,I-1)-U(1,I+1)) ) C C VERIFY THE CONVERGENCE IF(H(2,I).LE.0)THEN WRITE(*,100) 100 FORMAT(' ERROR ! THE COMPUTED DEPTH IS NEGATIVE !'/ 1 ' TRY TO RESTART THE PROGRAM USING'/ 2 ' A SMALLER TIME STEP'// 3 ' END OF PROGRAM !'//) WRITE(NUNIT,100) WRITE(*,*)'PRESS RETURN TO FINISH THE PROGRAM.' READ(*,*) STOP ENDIF C BETA = U(1,I) + 1 0.5 * (DT/DX) * U(1,I)*(U(1,I-1)-U(1,I+1)) + 2 0.5 * G * (DT/DX) * (H(1,I-1)-H(1,I+1)) + 3 G * DT * JF RH = (B + M * H(2,I)) * H(2,I) / (B + 2 * H(2,I) * SQRT(1 + M**2)) GAMA = RH**(4./3.) / (N**2 * G * DT) U(2,I) = 0.5 * (SQRT(GAMA**2 + 4.0*GAMA*BETA) - GAMA) C C VERIFY THE COURANT CONDITION IF(DX/(ABS(U(2,I))+SQRT(H(2,I)*G)).LT.4*DT)THEN WRITE(*,200) 200 FORMAT(' ERROR ! THE COURANT CONDITION IS NOT MET'/ 1 ' TRY TO RESTART THE PROGRAM USING'/ 2 ' A SMALLER TIME STEP'// 3 ' END OF PROGRAM !'//) WRITE(NUNIT,200) WRITE(*,*)'PRESS RETURN TO FINISH THE PROGRAM.' READ(*,*) STOP ENDIF C S = H(2,I) * (B + M * H(2,I)) C Q(2,I) = U(2,I) * S C 300 CONTINUE C RETURN END CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC SUBROUTINE AVAL( NNMAX , NN , H , U , Q , 1 B , M , JF , N , DX , DT ) C******************************************************************************* C * C * C******************************************************************************* C C DECLARATION OF VARIABLES REAL M , N , JF DIMENSION H(2,NNMAX) , U(2,NNMAX) , Q(2,NNMAX) C H(2,NN) = H(1,NN) + (DT/DX) * 1 ( U(1,NN)*(H(1,NN-1)-H(1,NN)) + 2 H(1,NN)*(U(1,NN-1)-U(1,NN)) ) C S = H(2,NN) * (B + M * H(2,NN)) P = B + 2 * H(2,NN)* SQRT(1 + M**2) RH = S / P C U(2,NN) = RH**(2./3.) * SQRT(JF) / N C Q(2,NN) = U(2,NN) * S C RETURN END