1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
| 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 = 'NEW')
C
C AFFICHER LE TITRE DU PROGRAMME
C
WRITE(*,*)
WRITE(7,*)
WRITE(*,*)
WRITE(7,*)
WRITE(*,10)
WRITE(7,10)
10 FORMAT(//' PROGRAMME DE CALCUL POUR ECOULEMENTS NON PERMANENTS'
1 //' SYSTEME D''UNITE : METRIQUE'//)
C
C LECTURE DES DONNEES CONCERNANT LE CANAL
C
WRITE(*,20)
WRITE(7,20)
20 FORMAT(' LECTURE DES DONNEES CONCERNANT LE CANAL :'/
1 ' LE CANAL PEUT AVOIR UNE SECTION TRAPEZOIDALE, RECTANGULA'
2,'IRE (m = 0) OU'/
3' TRIANGULAIRE (B = 0).'//)
30 WRITE(*,40)
WRITE(7,40)
40 FORMAT(' LARGEUR DU FOND DU CANAL ? B (m) = ',$)
READ(*,*,ERR=30)B
50 WRITE(*,60)
WRITE(7,60)
60 FORMAT(' PENTE DES BERGES ? m (-) = ',$)
READ(*,*,ERR=50)M
70 WRITE(*,80)
WRITE(7,80)
80 FORMAT(' PENTE DU FOND DU CANAL ? Jf (-) = ',$)
READ(*,*,ERR=70)JF
90 WRITE(*,100)
WRITE(7,100)
100 FORMAT(' COEFFICIENT DE MANNING ? n (m^-1/3 s) = ',$)
READ(*,*,ERR=90)N
110 WRITE(*,120)
WRITE(7,120)
120 FORMAT(' PROFONDEUR POUR ECOULEMENT UNIFORME ? hn (m) = ',$)
READ(*,*,ERR=110)HN
130 WRITE(*,140)
WRITE(7,140)
140 FORMAT(' LONGUEUR TOTALE DU 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(//' EN UTILISANT LA FORMULE DE MANNING (EQ. 3.16 OU 3.33)'/
1 ' LE DEBIT D''ECOULEMENT UNIFORME EST DONC, Qo (m3/s) = ',
2 F8.3)
C
C LECTURE DES DONNEES CONCERNANT LES CONDITIONS AUX LIMITES
C
WRITE(*,200)
WRITE(7,200)
200 FORMAT(///' LECTURE DES DONNEES CONCERNANT LES CONDITIONS AUX',
1 ' LIMITES (AMONT) :'/
2 ' ON ADMET UN HYDROGRAMME TRIANGULAIRE',
3 ' DEFINI PAR TROIS PARAMETRES :'/)
210 WRITE(*,220)
WRITE(7,220)
220 FORMAT(' TEMPS DE MONTEE ? t'' (s) = ',$)
READ(*,*,ERR=210)TP
230 WRITE(*,240)
WRITE(7,240)
240 FORMAT(' DEBIT DE POINTE ? QMAX (m3/s) = ',$)
READ(*,*,ERR=230)QMAX
250 WRITE(*,260)
WRITE(7,260)
260 FORMAT(' TEMPS DE DESCENTE ? t'''' (s) = ',$)
READ(*,*,ERR=250)TPD
TTH = TP + TPD
WRITE(*,270) TTH
WRITE(7,270) TTH
270 FORMAT(' ATTENTION ! ...'/' POUR t > (t'' + t'''') = ',F10.3,
1 ' (s), ON ADMET 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(/' LECTURE DES DONNEES CONCERNANT LE PAS D''ESPACE, DX :'/
1 ' LA LONGUEUR DU CANAL EST ',
2 ' LT (m) = ',F7.2/
3 ' NOMBRE MAXIMAL DE TRONCONS AUTORISE PAR LE PROGRAMME, ',
4 ' MAXDIV = ',I3/
5 ' LA VALEUR MINIMALE DE DX EST DONC, ',
6 ' LT / MAXDIV = ',F7.2, '(m)'//
7 ' (Note importante! Si vous voulez avoir davantage de',
8 ' troncons, modifiez'/' la valeur du parametre MAXDIV',
9 ' dans le programme principal)'
A //' ETES-VOUS D''ACCORD AVEC CETTE VALEUR ? REPONDRE',
B ' O/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(//' LECTURE DES DONNEES CONCERNANT LE PAS DE TEMPS, DT :'//
1 ' LE DEBIT D''ECOULEMENT UNIFORME EST, Qo',
2 ' (m3/s) = ',F8.3/
3 ' LA VITESSE MOYENNE DE L''ECOULEMENT EST, Uo (m/s) = Q',
4 ' / S = ',F8.3/
5 ' LA CELERITE DES ONDES DE GRAVITE EST, Co (m/s) =',
6 ' (g*hn)^0.5 = ',F8.3//
7 ' LA CONDITION DE STABILITE DE COURANT EXIGE QUE DT < (DX',
8 ' / (ABS(Uo) + Co)'/
9 ' LA VALEUR MAXIMALE DE DT EST DONC: DTMAX= DX / (ABS(Uo)',
A ' + Co) : ',F8.3,' s'//)
420 WRITE(*,430)
WRITE(7,430)
430 FORMAT(' QUELLE EST LA VALEUR 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(//' DONNER LA DUREE DES CALCULS, 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(//' VOUS AVEZ ',I3,' STATIONS LE LONG DU CANAL;'/
1 ' OU VOULEZ-VOUS IMPRIMER LES RESULTATS ?'/
2 ' ECRIRE LES NUMEROS DE 3 STATIONS EN FORMAT LIBRE = ',$)
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(//' A QUELLE FREQUENCE DE PAS VOULEZ-VOUS ECRIRE LES',
1 ' RESULTATS ? = ',$)
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(//' NOM DU FICHIER DE SORTIE ? = ',$)
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(//' FIN NORMALE DU PROGRAMME')
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)/' RESULTAT DES CALCULS D
1''ECOULEMENT NON PERMANENT AVEC LA METHODE EXPLICITE'/
DATA TITRE(2)/' ======================
1==================================================='/
DATA TITRE(3),TITRE(4)/' ',' '/
DATA TITRE(5)/' DONNEES CONCE
1RNANT LE CANAL ET L''ECOULEMENT UNIFORME'/
DATA TITRE(6)/' -------------
1---------------------------------------'/
DATA TITRE(7)/' LARGEUR DU FOND DU CANAL, B (m) =
1 COEFFICIENT DE MANNING, n (m^-1/3 s) = '/
DATA TITRE(8)/' PENTE DES BERGES, m (-) =
1 PROFONDEUR POUR ECOUL. UNIFORME, hn (m) = '/
DATA TITRE(9)/' PENTE DU FOND DU CANAL, Jf (-) =
1 LONGUEUR TOTALE DU CANAL, LT (m) = '/
DATA TITRE(10)/' '/
DATA TITRE(11)/' DEBIT INITIAL, Q0 (M3/S) = XXXX
1.XXX VITESSE INITIALE, U0 (m/s) = '/
DATA TITRE(12)/' CELERITE DES O
1NDES DE GRAVITE C0 (m/s) = '/
DATA TITRE(13)/' '/
DATA TITRE(14)/'
1 TEMPS DE MONTEE, t'' (s) = '/
DATA TITRE(15)/' DEFINITION DE L''HYDROGRAMME TRIANGULAIRE :
1 DEBIT DE POINTE, QMAX (m3/s) = '/
DATA TITRE(16)/'
1 TEMPS DE DESCENTE, t'''' (s) = '/
DATA TITRE(17)/' '/
DATA TITRE(18)/' NOMBRE DE TRONCONS, NDIV = LONGUEUR D
1ES TRONCONS, DX (m) = NOMBRE DE NOEUDS, NN=NDIV+1
2 = '/
DATA TITRE(19)/' PAS DE TEMPS, DT (s) = DUREE DES
1CALCULS, TMAX (s) = FREQUENCE D''ECRITURE (pas)
2 = '/
DATA TITRE(20)/' '/
DATA TITRE(21)/' STATION 1
1 STATION 2 STATION 3'/
DATA TITRE(22)/' TEMPS NOEUD NO = X =
1 NOEUD NO = X = NOEUD 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
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 |
[Fortran 77] Tirer un algorithme d'un code Fortran pour le traduire en C
Répondre avec citation
Partager