      PROGRAM EXPLIC
C*******************************************************************************
C                                                                              *
C PROGRAMME PRINCIPAL POUR LA SOLUTION DES ECOULEMENTS NON PERMANENTS          *
C DANS LES CANAUX A SECTION TRAPEZOIDALE, RECTANGULAIRE OU TRIANGULAIRE        *
C                                                                              *
C*******************************************************************************
C
C LISTE DES VARIABLES SIMPLES, DES VECTEURS ET DES MATRICES POUR LE PROGRAMME
C PRINCIPAL ET LES SOUS-PROGRAMMES.
C 
C TYPE NOM   DIMEN.    EXPLICATIONS
C ==== ===== =======   =========================================================
C R*4  B             = LARGEUR DU FOND DU CANAL
C R*4  BETA          = VOIR EQ. 5.24a
C R*4  C0            = CELERITE DES ONDES DE GRAVITE AU TEMPS T = 0
C R*4  DT            = PAS DE TEMPS
C R*4  DTMAX         = PAS DE TEMPS LIMITE SELON LA COND. DE STAB. DE COURANT
C R*4  DX            = PAS D'ESPACE
C CHAR FICH          = NOM DU FICHIER DE SORTIE
C R*4  G             = ACCELERATION GRAVITATIONNELLE (G = 9.81 m/s2)
C R*4  GAMA          = VOIR EQ. 5.24
C R*4  H     (2,101) = PROFOND. AUX NOEUDS POUR LES DEUX DERNIERS PAS DE TEMPS
C R*4  HN            = PROFONDEUR NORMALE
C I*4  I             = COMPTEUR DE NUMERO DES NOEUDS
C I*4  IL            = COMPTEUR DE NUMERO DES LIGNES DU TITRE
C R*4  JF            = PENTE DU FOND
C I*4  K             = COMPTEUR DES BOUCLES
C R*4  LT            = LONGUEUR TOTALE DU CANAL
C R*4  M             = PENTE DES BERGES DU CANAL
C I*4  MAXDIV        = NOMBRE MAX. DE TRONCONS AUTORISES PAR LE PROGRAMME (100)
C R*4  N             = COEFFICIENT DE MANNING
C I*4  NDIV          = NOMBRE DE TRONCONS
C I*4  NN            = NOMBRE DE NOEUDS
C I*4  NNMAX         = NOMBRE MAX. DE NOEUDS AUTORISES PAR LE PROGRAMME (100+1)
C I*4  NOST  (3)     = NUMERO DES NOEUDS AVEC IMPRESSION DES RESULTATS
C I*4  NP            = NUMERO DE PAS DE TEMPS AU TEMPS T
C I*4  NPSTEP        = FREQUENCE (NOMBRE DE PAS) D'ECRITURE DES RESULTATS
C I*4  NUNIT         = NUMERO DE L'UNITE DE SORTIE (POUR LE FICHIER DE SORTIE)
C R*4  P             = PERIMETRE MOUILLE
C R*4  Q     (2,101) = DEBITS AUX NOEUDS POUR LES DEUX DERNIERS PAS DE TEMPS
C R*4  QE            = DEBIT ENTRANT AU TEMPS T
C R*4  Q0            = DEBIT INITIAL (DEBIT D'ECOULEMENT UNIFORME)
C R*4  QMAX          = DEBIT DE POINTE DE LA CRUE
C CHAR REP           = REPONSE ALPHANUMERIQUE
C R*4  RH            = RAYON HYDRAULIQUE
C R*4  S             = SURFACE MOUILLEE
C R*4  T             = DUREE DES CALCULS
C CHAR TITRE (24)    = TITRE A IMPRIMER AU DEBUT DU FICHIER DE SORTIE
C R*4  TMAX          = LIMITE DE TEMPS POUR LES CALCULS
C R*4  TP            = DUREE DE MONTEE DE L'HYDROGRAMME
C R*4  TPD           = DUREE DE DESCENTE DE LA CRUE
C R*4  TTH           = DUREE TOTALE DE LA CRUE (= TP + TPD)
C R*4  U     (2,101) = VITESSES AUX NOEUDS POUR LES DEUX DERNIERS PAS DE TEMPS
C R*4  U0            = VITESSE INITIALE (VITESSE D'ECOULEMENT UNIFORME)
C
C*******************************************************************************
C
C PARAMETRES
      PARAMETER ( NUNIT = 5 )
      PARAMETER ( MAXDIV = 100 , NNMAX = MAXDIV+1 )
      PARAMETER ( G = 9.81 )
C
C DECLARATION DES VARIABLES
      CHARACTER*50  FICH
      REAL          M , JF , N , LT
      DIMENSION     H(2,NNMAX) , U(2,NNMAX) , Q(2,NNMAX) , NOST(3)
C
C APPEL AU SOUS-PROGRAMME "LIRE" POUR LIRE LES DONNEES DU PROBLEME ET
C AUTRES PARAMETRES NECESSAIRES POUR LE PROGRAMME
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 APPEL AU SOUS-PROGRAMME "INIT" POUR CALCULER LES CONDITIONS INITIALES
C A TOUS LES NOEUDS AU TEMPS T = 0
C
      CALL INIT ( NNMAX , NN , H , U , Q , HN , U0 , Q0 )
C
C APPEL AU SOUS-PROGRAMME "CALCUL" POUR CALCULER L'ECOULEMENT NON PERMANENT
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 FIN DU PROGRAMME
C
      STOP
      END
      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                                                                              *
C SOUS-PROGRAMME POUR LIRE INTERACTIVEMENT:                                    *
C 1- LES DONNEES DU PROBLEME A RESOUDRE,                                       *
C 2- LES PARAMETRES NECESSAIRES POUR LE PROGRAMME                              *
C                                                                              *
C*******************************************************************************
C
C DECLARATION DES 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 = 'unknown')
C
C AFFICHER LE TITRE DU PROGRAMME
C
      WRITE(*,*)
      WRITE(7,*)
      WRITE(*,*)
      WRITE(7,*)
      WRITE(*,10)
      WRITE(7,10)
10    FORMAT(//'  PROGRAMA DE CALCULO PARA ESCOAEMENTOS NAO PERMANENTES'
     1       //'           SISTEMA DE UNIDADES : METRICO'//)
C
C LECTURE DES DONNEES CONCERNANT LE CANAL
C
      WRITE(*,20)
      WRITE(7,20)
20    FORMAT(' LEITURA DOS DADOS RELATIVOS AO CANAL :'/
     1       'O CANAL PODE TER UMA SECAO TRAPEZOIDAL,RETANGULAR (Z=0)OU'
     2       /' TRIANGULAR (B= 0 )'//)
    
30    WRITE(*,40)
      WRITE(7,40)
40    FORMAT(' LARGURA DE FUNDO DO CANAL ?                  B (m) = ',$)
      READ(*,*,ERR=30)B
50    WRITE(*,60)
      WRITE(7,60)
60    FORMAT(' INCLINACAO DOS TALUDES ?                     m (-) = ',$)
      READ(*,*,ERR=50)M
70    WRITE(*,80)
      WRITE(7,80)
80    FORMAT(' DECLIVIDADE DE FUNDO DO CANAL ?             Jf (-) = ',$)
      READ(*,*,ERR=70)JF
90    WRITE(*,100)
      WRITE(7,100)
100   FORMAT(' COEFICIENTE DE MANNING ?              n (m^-1/3 s) = ',$)
      READ(*,*,ERR=90)N
110   WRITE(*,120)
      WRITE(7,120)
120   FORMAT(' PROFUNDIDADE NO ESCOAMENTO UNIFORME ?       hn (m) = ',$)
      READ(*,*,ERR=110)HN
130   WRITE(*,140)
      WRITE(7,140)
140   FORMAT(' COMPRIMENTO TOTAL DO CANAL ?                LT (m) = ',$)
      READ(*,*,ERR=130)LT
C
C CALCUL DU DEBIT D'ECOULEMENT UNIFORME
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(//' UTILIZANDO A FORMULA DE MANNING '/ 
     1       ' A VAZAO DO ESCOAMENTO UNIFORME VALE       Qo (m3/s) = ',
     2       F8.3)
C
C LECTURE DES DONNEES CONCERNANT LES CONDITIONS AUX LIMITES
C
      WRITE(*,200)
      WRITE(7,200)
200   FORMAT(///' LEITURA DOS DADOS RELATIVOS AS CONDICOES ',
     1       ' LIMITES (MONTANTE) :'/
     2       ' ADMITE-SE UM HIDROGRAMA TRIANGULAR',
     3       ' DEFINIDO POR TRES PARAMETROS :'/)
210   WRITE(*,220)
      WRITE(7,220)
220   FORMAT(' TEMPO DE ASCENSAO ?          t'' (s) = ',$)
      READ(*,*,ERR=210)TP
230   WRITE(*,240)
      WRITE(7,240)
240   FORMAT(' VAZAO DE PICO ?          QMAX (m3/s) = ',$)
      READ(*,*,ERR=230)QMAX
250   WRITE(*,260)
      WRITE(7,260)
260   FORMAT(' TEMPO DE DEPLECAO ?         t'''' (s) = ',$)
      READ(*,*,ERR=250)TPD
      TTH = TP + TPD
      WRITE(*,270) TTH
      WRITE(7,270) TTH
270   FORMAT(' ATENCAO ! ...'/' PARA t > (t'' + t'''') = ',F10.3,
     1       ' (s), ADMITE-SE QUE Q = Qo'///)
C
C LECTURE DES PAS D'ESPACE ET DE TEMPS
C
      NDIV = MAXDIV
      NN   = NNMAX
      DX   = LT / MAXDIV
300   WRITE(*,310) LT , MAXDIV , DX
      WRITE(7,310) LT , MAXDIV , DX
310   FORMAT(/' LEITURA DOS DADOS RELATIVOS AO PASSO DE ESPACO, Dx :'/
     1       ' O COMPRIMENTO DO CANAL EH                              ',
     2       ' LT (m) = ',F7.2/
     3       ' NUMERO MAXIMO DE TRECHOS AUTORIZADO PARA O PROGRAMA, ',
     4       ' MAXDIV = ',I3/
     5       ' O VALOR MINIMO DE Dx EH DADO POR,               ',
     6       ' LT / MAXDIV = ',F7.2, '(m)'//
     7       ' (Nota importante! Se voce tem necessidade de mais',
     8       ' trechos, modifique'/'  o valor do parametro MAXDIV',
     9       ' no programa principal)'
     A     //' VOCE ESTA DE ACORDO COM ESTE VALOR ? REPONDA',
     B       ' S/CR OU N = ',$)
      READ(*,5)REP
      IF(REP.EQ.'N'.OR.REP.EQ.'n')THEN
320     WRITE(*,330)
        WRITE(7,330)
330     FORMAT(' QUEL EST LE NOMBRE DE TRONCONS,  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(/' LA VALEUR DE DX EST DONC,                        LT',
     1         ' / NDIV = ',F7.2, '(m)'/
     2         ' LE NOMBRE DE NOEUDS EST,                    NN =',
     3         ' NDIV + 1 = ',I3)
      ENDIF
C
      C0    = SQRT(G * (HN*(B+M*HN)) / (B+2*M*HN))
      DTMAX = DX / (U0 + C0)
      DT    = DTMAX
400   WRITE(*,410) Q0 , U0 , C0 , DTMAX
      WRITE(7,410) Q0 , U0 , C0 , DTMAX
410   FORMAT(//' LEITURA DOS DADOS RELATIVOS AO PASSO DE TEMPO, DT :'//
     1      ' A VAZAO NO ESCOAMENTO UNIFORME EH,                 Qo',
     2      ' (m3/s) = ',F8.3/
     3      ' A VELOCIDADE MEDIA DO ESCOAMENTO EH, Uo (m/s) = Q',
     4      ' / S      = ',F8.3/
     5      ' A CELERIDADE DAS ONDAS DE GRAVIDADE EH,   Co (m/s) =',
     6      ' (g*hn)^0.5 = ',F8.3//
     7      ' A CONDICAO DE ESTABILIDADE DE COURANT EXIGE QUE DT < (DX',
     8      ' / (ABS(Uo) + Co)'/
     9      ' O VALOR MAXIMO DE DT EH DADO POR: DTMAX= DX / (ABS(Uo)',
     A      ' + Co) : ',F8.3,' s'//)
420   WRITE(*,430)
      WRITE(7,430)
430   FORMAT(' QUAL EH O VALOR DE DT (s) ? = '$)
      READ(*,*,ERR=420)DT
      IF(DT.LE.0.0.OR.DT.GT.DTMAX) GO TO 400
C
C QUAND FAUT-IL TERMINER LES CALCULS ?
440   WRITE(*,450)
      WRITE(7,450)
450   FORMAT(//' QUAL A DURACAO DOS CALCULOS, TMAX (s) ? = ',$)
      READ(*,*,ERR=440)TMAX
      IF(TMAX.LT.DT) GO TO 440
C
C COMMENT IMPRIMER LES RESULTATS ?
500   WRITE(*,510)NN
      WRITE(7,510)NN
510   FORMAT(//' VOCE TEM ',I3,' ESTACOES AO LONGO DO CANAL;'/
     1         ' VOCE DESEJA IMPRIMIR OS RESULTADOS ?'/
     2         ' ESCREVA OS NUMEROS DE 3 ESTACOES EM FORMATO LIVRE =',$)
      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 FREQUENCE D'ECRITURE DES RESULTATS
550   WRITE(*,560)
      WRITE(7,560)
560   FORMAT(//' COM QUE FREQUENCIA VOCE QUER ESCREVER OS',
     1         ' RESULTADOS ? (NUMERO DE PASSOS DE TEMPO) = ',$)
      READ(*,*,ERR=550)NPSTEP
      IF(NPSTEP.LT.1 .OR. NPSTEP.GT.(TMAX/DT))GO TO 550
C
C FICHIER DE SORTIE
600   WRITE(*,610)
      WRITE(7,610)
610   FORMAT(//' NOME DO ARQUIVO DE SAIDA ? = ',$)
      READ(*,5) FICH
      OPEN( UNIT = NUNIT , FILE = FICH , STATUS = 'NEW' , ERR = 600 )
C
      RETURN
      END
      SUBROUTINE INIT ( NNMAX , NN , H , U , Q , HN , U0 , Q0 )
C*******************************************************************************
C                                                                              *
C SOUS-PROGRAMME POUR L'ATTRIBUTION DES CONDITIONS INITIALES (T=0) AUX NOEUDS. *
C                                                                              *
C*******************************************************************************
C
C DECLARATION DES 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
      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 SOUS-PROGRAMME POUR CALCULER (PAR PAS DE TEMPS) LA PROFONDEUR, LA VITESSSE   *
C ET LE DEBIT A TOUS LES NOEUDS POUR UN ECOULEMENT NON PERMANENT SELON LA      *
C METHODE EXPLICITE (VOIR CHAP. 5.2.3).                                        *
C                                                                              *
C*******************************************************************************
C
C DECLARATION DES VARIABLES
      REAL       M , N , JF
      DIMENSION  H(2,NNMAX) , U(2,NNMAX) , Q(2,NNMAX) , NOST(3)
C
C INITIALISER LE TEMPS
      T = 0
C
C ECRIRE LE TITRE DE LA SORTIE
      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 IMPRIMER LES RESULTATS AU TEMPS T
100   CALL ECRIT ( NNMAX , NN , H , U , Q , T , DT ,
     1             NOST , NPSTEP , NUNIT )
C
C INCREMENTER LE TEMPS SI LES PAS DE TEMPS NE SONT PAS EPUISES
      IF(T+DT.GT.TMAX) GO TO 900
      T = T + DT
C
C CALCULER LE DEBIT ENTRANT EN AMONT AU NOUVEAU TEMPS T
      CALL QENTRE ( TP , QMAX , TPD , T , Q0 , QE )
C
C CALCULER LE NOEUD EN AMONT DU CANAL AU NOUVEAU TEMPS T
      CALL AMONT( NNMAX , NN , H , U , Q ,
     1            B , M , QE , DX , DT )
C
C CALCULER LES POINTS INTERMEDIAIRES AU NOUVEAU TEMPS 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 REMPLACER LES VALEURS AU PAS DE TEMPS PRECEDENT PAR DES VALEURS CALCULEES
      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 ALLER INCREMENTER LE TEMPS ET REPETER LES CALCULS
C SI, LE NOMBRE DE PAS DE TEMPS SPECIFIE N'EST PAS EPUISE.
      GO TO 100
C
900   WRITE(*,910)
      WRITE(NUNIT,910)
910   FORMAT(//' FIM NORMAL DO PROGRAMA')
      RETURN
      END
      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 CE SOUS-PROGRAMME ECRIT LE TITRE ET LES DONNEES DU PROBLEME SUR LE FICHIER   *
C DE SORTIE.                                                                   *
C                                                                              *
C*******************************************************************************
C
C DECLARATION DES VARIABLES
      CHARACTER*131 TITRE(24)
      REAL           M , JF , N , LT
      DIMENSION NOST(3)
C
      DATA TITRE(1)/'                             RESULTADO DOS CALCULOS 
     1 DO ESCOUAMENTO NAO PERMANENTE COM O METODO EXPLICITO'/
      DATA TITRE(2)/'                             ======================
     1==================================================='/
      DATA TITRE(3),TITRE(4)/' ',' '/
      DATA TITRE(5)/'                                      DADOS RELATIV
     1OS AO CANAL E AO ESCOAMENTO UNIFORME'/
      DATA TITRE(6)/'                                      -------------
     1---------------------------------------'/
      DATA TITRE(7)/' LARGURA DO FUNDO DO CANAL,          B (m) =        
     1              COEFICIENTE DE MANNING,     n (m^-1/3 s) = '/
      DATA TITRE(8)/' INCLINACAO DOS TALUDES,            Z (-) =        
     1              PROFUNDIDADE PARA ESCOA. UNIFORME,  hn (m) = '/
      DATA TITRE(9)/' DECLIVIDADE DO FUNDO DO CANAL,    Jf (-) =        
     1              COMPRIMENTO TOTAL DO CANAL,         LT (m) = '/
      DATA TITRE(10)/' '/
      DATA TITRE(11)/'                 VAZAO INICIAL,   Q0 (m3/S) = XXXX
     1.XXX                 VELOCIDADE INICIAL, U0 (m/s) = '/
      DATA TITRE(12)/'                                    CELERIDADE DAS 
     1 ONDAS DE GRAVIDADE C0 (m/s) = '/
      DATA TITRE(13)/' '/
      DATA TITRE(14)/'                                                  
     1   TEMPO DE ASCENSAO,      t'' (s) = '/
      DATA TITRE(15)/' DEFINICAO DO HIDROGRAMA TRIANGULAR :      
     1    VAZAO DE PICO, QMAX (m3/s)    = '/
      DATA TITRE(16)/'                                                  
     1   TEMPO DE DEPLECAO,   t'''' (s) = '/
      DATA TITRE(17)/' '/
      DATA TITRE(18)/' NUMERO DE TRECHOS, NDIV =             COMPRIMENTO
     1 DOS TRECHOS, DX (m) =                 NUMERO DE NOS, NN=NDIV+1
     2 = '/
      DATA TITRE(19)/' PASSO DE TEMPO,     DT (s) =             DURACAO 
     1DOS CALCULOS,  TMAX (s) =              FREQUECIA DE ESCRITA (pas)
     2  = '/
      DATA TITRE(20)/' '/
      DATA TITRE(21)/'                             ESTACAO 1            
     1                ESTACAO 2                            ESTACAO 3'/
      DATA TITRE(22)/'      TEMPO        NO NUMERO =         X =         
     1      NO NUMERO =         X =               NO NUMERO =         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
      SUBROUTINE ECRIT ( NNMAX , NN , H , U , Q , T , DT ,
     1                   NOST , NPSTEP , NUNIT )
C*******************************************************************************
C                                                                              *
C SOUS-PROGRAMME POUR IMPRIMER LES VALEURS DE H , U , ET Q AU TEMPS T.         *
C                                                                              *
C*******************************************************************************
C
C DECLARATION DES 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
      SUBROUTINE QENTRE ( TP , QMAX , TPD , T , Q0 , QE )
C*******************************************************************************
C                                                                              *
C SOUS-PROGRAMME POUR CALCULER LE DEBIT ENTRANT PAR L'AMONT AU NOUVEAU TEMPS T *
C PAR INTERPOLATION LINEAIRE.  LA FORME DE L'HYDROGRAMME EST TRIANGULAIRE.     *
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
      SUBROUTINE AMONT( NNMAX , NN , H , U , Q ,
     1                  B , M , QE , DX , DT )
C*******************************************************************************
C                                                                              *
C SOUS-PROGRAMME POUR CALCULER LA PROFONDEUR , LA VITESSE ET LE DEBIT          *
C AU NOEUD 1, QUI SE TROUVE TOUT EN AMONT DU CANAL, AU NOUVEAU TEMPS.          *
C LES VALEURS POUR LE PAS DE TEMPS PRECEDENT SONT STOCKEES SUR LA LIGNE        *
C I = 1 DES TABLEAUX H , U ET Q. LES VALEURS CALCULEES POUR LE NOUVEAU         *
C PAS DE TEMPS SONT STOCKEES SUR LA LIGNE 2.                                   *
C                                                                              *
C*******************************************************************************
C
C DECLARATION DES VARIABLES
      REAL      M
      DIMENSION H(2,NNMAX) , U(2,NNMAX) , Q(2,NNMAX)
C
C CALCULER D'ABORD LA PROFONDEUR  AU NOEUD "1" AU NOUVEAU TEMPS
      H(2,1) = H(1,1) + (DT/DX) * 
     1         ( U(1,1)*(H(1,1)-H(1,2)) +
     2         (H(1,1)*(B+M*H(1,1))/(B+2*M*H(1,1)))*(U(1,1)-U(1,2)) )
C
C CALCULER ENSUITE LA SURFACE MOUILLEE POUR CETTE PROFONDEUR
      S  = H(2,1) * (B + M * H(2,1))
C
C LE DEBIT  AU NOEUD "1" AU NOUVEAU TEMPS EST EGAL AU DEBIT ENTRANT
      Q(2,1) = QE
C
C VITESSE  AU NOEUD "1" AU NOUVEAU TEMPS
      U(2,1) = Q(2,1) / S
C
      RETURN
      END
      SUBROUTINE INTER( NNMAX , NN , H , U , Q ,
     1                  B , M , JF , N , G , DX , DT , NUNIT )
C*******************************************************************************
C                                                                              *
C SOUS-PROGRAMME POUR CALCULER LA PROFONDEUR , LA VITESSE ET LE DEBIT          *
C AUX NOEUDS INTERMEDIAIRES (2 A NN-1) AU NOUVEAU TEMPS.                       *
C LES VALEURS POUR LE PAS DE TEMPS SOUS-PROG SONT STOCKEES SUR LA LIGNE        *
C I = 1 DES TABLEAUX H , U ET Q. LES VALEURS CALCULEES POUR LE NOUVEAU         *
C PAS DE TEMPS SONT STOCKEES SUR LA LIGNE 2.                                   *
C                                                                              *
C*******************************************************************************
C
C DECLARATION DES VARIABLES
      REAL      M , N , JF
      DIMENSION H(2,NNMAX) , U(2,NNMAX) , Q(2,NNMAX)
C
      DO 300 I = 2 , NN-1
C
C CALCULER D'ABORD LA PROFONDEUR AU NOEUD "NN" AU NOUVEAU TEMPS
      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)*(B+M*H(1,I))/(B+2*M*H(1,I)))*(U(1,I-1)-U(1,I+1)) )
C
C VERIFIER LA CONVERGENCE
      IF(H(2,I).LE.0)THEN
        WRITE(*,100)
100     FORMAT(' ERREUR ! LA PROFONDEUR CALCULEE EST NEGATIVE !'/
     1         ' ESSAYER DE RELANCER LE PROGRAMME EN UTILISANT'/
     2         ' UN PAS DE TEMPS PLUS PETIT'//
     3         ' FIN ANORMALE DU PROGRAMME !'//)
        WRITE(NUNIT,100)
        WRITE(*,*)'TAPER RETURN POUR TERMINER LE PROGRAMME.'
        READ(*,*)
        STOP
      ENDIF
C
C CALCULER ENSUITE LA VITESSE
      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 VERIFIER LA CONDITION DE COURANT
      IF(DX/(ABS(U(2,I))+SQRT((H(2,I)*(B+M*H(2,I))/(B+2*M*H(2,I)))*G))
     1.LT.4*DT)THEN
        WRITE(*,200)
200     FORMAT(' ERREUR ! LA CONDITION DE COURANT NE PEUT PAS'/
     1         ' ETRE RESPECTEE. ESSAYER DE LANCER LE PROGRAMME'/
     2         ' AVEC UN PAS DE TEMPS PLUS PETIT'//
     3         ' FIN ANORMALE DU PROGRAMME !'//)
        WRITE(NUNIT,200)
        WRITE(*,*)'TAPER RETURN POUR TERMINER LE PROGRAMME.'
        READ(*,*)
        STOP
      ENDIF
C
C CALCULER LA SURFACE MOUILLEE
      S  = H(2,I) * (B + M * H(2,I))
C
C DEBIT AU NOEUD I AU NOUVEAU TEMPS
      Q(2,I) = U(2,I) * S
C
300   CONTINUE
C
      RETURN
      END
      SUBROUTINE AVAL( NNMAX , NN , H , U , Q ,
     1                 B , M , JF , N , DX , DT )
C*******************************************************************************
C                                                                              *
C SOUS-PROGRAMME POUR CALCULER LA PROFONDEUR , LA VITESSE ET LE DEBIT          *
C AU NOEUD NN, QUI SE TROUVE TOUT EN AVAL DU CANAL, AU NOUVEAU TEMPS.          *
C LES VALEURS POUR LE PAS DE TEMPS PRECEDENT SONT STOCKEES SUR LA LIGNE        *
C I = 1 DES TABLEAUX H , U ET Q. LES VALEURS CALCULEES POUR LE NOUVEAU         *
C PAS DE TEMPS SONT STOCKEES SUR LA LIGNE 2.                                   *
C                                                                              *
C*******************************************************************************
C
C DECLARATION DES VARIABLES
      REAL      M , N , JF
      DIMENSION H(2,NNMAX) , U(2,NNMAX) , Q(2,NNMAX)
C
C CALCULER D'ABORD LA PROFONDEUR AU NOEUD "NN" AU NOUVEAU TEMPS
      H(2,NN) = H(1,NN) + (DT/DX) * 
     1         ( U(1,NN)*(H(1,NN-1)-H(1,NN)) +
     2   (H(1,NN)*(B+M*H(1,NN))/(B+2*M*H(1,NN)))*(U(1,NN-1)-U(1,NN)) )
C
C CALCULER ENSUITE LA SURFACE MOUILLEE ET LE PERIMETRE
      S   = H(2,NN) * (B + M * H(2,NN))
      P   = B + 2 * H(2,NN)* SQRT(1 + M**2)
      RH  = S / P
C
C CALCULER ENSUITE LA VITESSE AU NOEUD "NN" AU NOUVEAU TEMPS
C A L'AIDE DE LA FORMULE DE MANNING
      U(2,NN) = RH**(2./3.) * SQRT(JF) / N
C
C DEBIT AU NOEUD "NN" AU NOUVEAU TEMPS
      Q(2,NN) = U(2,NN) * S
C
      RETURN
      END

 