      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

