1 : C CKON4 SOURCE PV 09/03/12 21:17:08 6325
2 : SUBROUTINE CKON4(LOGME,INDMET,NORDP1,
3 : & IROF,IVITF,IPF,IFRMAF,ISCAF,PROPHY,
4 : & ICHPSU,ICHPDI,
5 : & MELEMC,MELEMF,MELEFE,
6 : & IZG1,IZG2,IZG3,IZG4,IZG5,DT,DIAMEL,NLCEMI,
7 : & LOGNC,LOGAN,MESERR)
8 : C************************************************************************
9 : C
10 : C PROJET : CASTEM 2000
11 : C
12 : C NOM : CKON4
13 : C
14 : C DESCRIPTION : Voir CKON
15 : C
16 : C Cas trois dimensions, gaz "thermally perfect"
17 : C
18 : C LANGAGE : FORTRAN 77 + ESOPE 2000 (avec estensions CISI)
19 : C
20 : C AUTEUR : R. MOREL, A. BECCANTINI, DRN/DMT/SEMT/TTMF
21 : C
22 : C************************************************************************
23 : C
24 : C
25 : C APPELES (Outils
26 : C CASTEM) : KRIPAD, LICHT
27 : C
28 : C APPELES (Calcul) : FVLHT2, FCOLT2
29 : C
30 : C
31 : C************************************************************************
32 : C
33 : C ENTREES
34 : C
35 : C
36 : C 1) PARAMETRES
37 : C
38 : C LOGME : (LOGICAL); .TRUE. -> MULTI-ESPECES
39 : C .FALSE. -> MONO-ESPECE
40 : C
41 : C INDMET : 1 VLH (van Leer Hanel FVS)
42 : C
43 : C 2 HUSVLH (VLH FVS + Colella-Glaz FDS)
44 : C
45 : C NORDP1 : (ordre des polynoms cv(T)) + 1
46 : C
47 : C 2) Pointeurs des MCHAMLs
48 : C
49 : C IROF : MCHAML sur "FACEL" contenant la masse volumique
50 : C ("gauche" et "droite");
51 : C
52 : C IVITF : MCHAML sur "FACEL" contenant la vitesse dans le repaire
53 : C local (n,t) et les cosinus directeurs des repaire local;
54 : C
55 : C IPF : MCHAML sur "FACEL" contenant la pression;
56 : C
57 : C IFRAMAF : MCHAML sur "FACEL", contenant les fractions massiques
58 : C si LOGME = .TRUE.;
59 : C LOGME = .FALSE. -> IFRAMAF = 0
60 : C
61 : C ISCAF : MCHAML sur "FACEL", contenant les scalires passifs a
62 : C transporter (ou ISCAF = 0)
63 : C
64 : C 3) Pointeur sur la table des proprietes de gaz
65 : C
66 : C PROPHY
67 : C
68 : C 4) Pointeurs de CHPOINTs de la table DOMAINE
69 : C
70 : C ICHPSU : CHPOINT "FACE" contenant la surface des faces
71 : C
72 : C ICHPDI : CHPOINT "CENTRE" contenant le diametre minimum
73 : C de chaque element
74 : C
75 : C
76 : C 5) Pointeurs de MELEME de la table DOMAINE
77 : C
78 : C MELEMC : MELEME 'CENTRE' du SPG des CENTRES
79 : C
80 : C MELEMF : MELEME 'FACE' du SPG des FACES
81 : C
82 : C MELEFE : MELEME 'FACEL' du connectivité FACES -> ELEM
83 : C
84 : C SORTIES (il faudrait dire E/S)
85 : C
86 : C IZGi : pointeurs de CHPOINTs "FACE" dont les valeurs
87 : C se trouvent dans la table KIZX . 'EQEX' . 'KIZG'
88 : C aux indices 'IZGi'
89 : C
90 : C IZG1 : Increment de masse voluique
91 : C
92 : C IZG2 : Increment de quantite de mouvement
93 : C
94 : C IZG3 : Increment de l'energie totale
95 : C
96 : C IZG4 : Increment de les Masse Volumiques des Especies
97 : C (si LOGME = .TRUE.)
98 : C
99 : C IZG5 : Increment de scalaires passifs
100 : C
101 : C DIAMEL : 'minimum' diametre du maillage
102 : C
103 : C NLCEMI : numero local du CENTRE ou le diametre est 'minimum'
104 : C
105 : C DT : pas de temps pour le respect de la CFL-like condition
106 : C DT < DIAMEL /2 /max(Lambda_i)
107 : C En maillage regulier cette condition garantie la
108 : C non-interaction des ondes
109 : C
110 : C
111 : C LOGNC : (LOGICAL): si .TRUE. la methode de Newton-Rapson, utilisée
112 : C dans pour la solution du probleme Riemann exacte ou dans
113 : C l'algorithm HUS, n'a pas bien marchéee; MESERR = 'Goudunov'
114 : C ou 'HUS'.
115 : C
116 : C LOGAN : (LOGICAL): si .TRUE. une anomalie à été detectée
117 : C
118 : C MESERR : pour l'ecriture des messages d'erreurs
119 : C
120 : C************************************************************************
121 : C
122 : C HISTORIQUE (Anomalies et modifications éventuelles)
123 : C
124 : C HISTORIQUE : 22.12.98 Creation
125 : C
126 : C 08.02.00 Ajout du flux numerique de Colella-Glaz
127 : C
128 : C 21.02.00 Ajout de transport de scalaires passifs
129 : C
130 : C************************************************************************
131 : C
132 : C
133 : C N.B.: On suppose qu'on a déjà controllé RO, P > 0
134 : C GAMMA \in (1,3)
135 : C Y \in (0,1)
136 : C Si non il faut le faire!!!
137 : C
138 : C************************************************************************
139 : C
140 : C
141 : C**** Variables de COOPTIO
142 : C
143 : C INTEGER IPLLB, IERPER, IERMAX, IERR, INTERR
144 : C & ,IOTER, IOLEC, IOIMP, IOCAR, IOACQ
145 : C & ,IOPER, IOSGB, IOGRA, IOSAU, IORES
146 : C & ,IECHO, IIMPI, IOSPI
147 : C & ,IDIM
148 : CC & ,MCOORD
149 : C & ,IFOMOD, NIFOUR, IFOUR, NSDPGE, IONIVE
150 : C & ,NGMAXY, IZROSF, ISOTYP, IOSCR,LTEXLU
151 : C & ,NORINC,NORVAL,NORIND,NORVAD
152 : C & ,NUCROU, IPSAUV
153 : C
154 : IMPLICIT INTEGER(I-N)
155 : INTEGER I1, I2
156 : & ,INDMET, NORDP1
157 : & ,IROF,IVITF,IPF,IFRMAF, ISCAF
158 : & ,ICHPSU,ICHPDI,MELEMC,MELEMF,MELEFE
159 : & ,IGEOMC,IGEOMF
160 : & ,IZG1,IZG2,IZG3,IZG4,IZG5,NLCEMI
161 : & ,NESP, NFAC, NSCA, NMA, NMA2
162 : & ,NLCF, NGCEG, NGCED, NLCEG, NLCED
163 : & ,NGCF, NLCF1, SPG1, SPG2
164 : REAL*8 DIAMEL, DT, UNSDT, CELLT
165 : & , ROG, UNG, UTG, UVG, PG, GAMG, TG
166 : & , ROD, UND, UTD, UVD, PD, GAMD, TD
167 : & , SURF,CNX, CNY, CNZ, CTX , CTY, CTZ
168 : & , CVX, CVY, CVZ
169 : & , CELL, DIAMG, DIAMD, DIAM
170 : & , ASON, LAMBDA
171 : & , Y1G, Y1D, RTOTG, RTOTD, CVTOTG, CVTOTD, ETHERG, ETHERD
172 : & , YG, YD, FUNTG, FUNTD, DCV, XST
173 : LOGICAL LOGME, LOGNC, LOGAN
174 : CHARACTER*(40) MESERR
175 : CHARACTER*(8) TYPE
176 : C
177 : C**** Segment des proprietes du gaz
178 : C
179 : SEGMENT PROPHY
180 : REAL*8 ACV(NORDP1,NESP+1), R(NESP+1), H0K(NESP+1)
181 : & ,ACVTOG(NORDP1), ACVTOD(NORDP1)
182 : ENDSEGMENT
183 : C
184 : C**** Les Includes
185 : C
186 : C
187 : C**** LES INCLUDES
188 : C
189 : -INC CCOPTIO
190 : -INC SMCOORD
191 : -INC SMCHAML
192 : POINTEUR MELVNX.MELVAL, MELVNY.MELVAL,MELVNZ.MELVAL,
193 : & MELT1X.MELVAL, MELT1Y.MELVAL,MELT1Z.MELVAL,
194 : & MELT2X.MELVAL, MELT2Y.MELVAL,MELT2Z.MELVAL
195 : POINTEUR MELVUN.MELVAL, MELVUT.MELVAL, MELVUV.MELVAL
196 : POINTEUR MELRO.MELVAL, MELP.MELVAL
197 : -INC SMCHPOI
198 : POINTEUR MPOVSU.MPOVAL, MPOVDI.MPOVAL
199 : & , MPOVG1.MPOVAL, MPOVG2.MPOVAL
200 : & , MPOVG3.MPOVAL, MPOVG4.MPOVAL, MPOVG5.MPOVAL
201 : POINTEUR MCHAMY.MCHAML, MCHAMS.MCHAML
202 : -INC SMELEME
203 : -INC SMLMOTS
204 : -INC SMLENTI
205 : C
206 : C**** Les fractiones massiques.
207 : C
208 : SEGMENT FRAMAS
209 : REAL*8 YET(NMA)
210 : ENDSEGMENT
211 : POINTEUR FRAMAG.FRAMAS, FRAMAD.FRAMAS, SCAG.FRAMAS, SCAD.FRAMAS
212 : C
213 : C**** Les flux aux interface dans le repaire (n,t)
214 : C
215 : SEGMENT IFLUX
216 : REAL*8 FLUX(NMA2)
217 : ENDSEGMENT
218 : POINTEUR IFLU1.IFLUX, IFLU2.IFLUX
219 : C
220 : C
221 : C**** Initialisation des MCHAMLs
222 : C
223 : C**** Masse volumique
224 : C
225 : MCHEL1 = IROF
226 : SEGACT MCHEL1
227 : MCHAM1 = MCHEL1.ICHAML(1)
228 : SEGACT MCHAM1
229 : MELRO = MCHAM1.IELVAL(1)
230 : SEGDES MCHEL1
231 : SEGDES MCHAM1
232 : C
233 : C**** Pression
234 : C
235 : MCHEL1 = IPF
236 : SEGACT MCHEL1
237 : MCHAM1 = MCHEL1.ICHAML(1)
238 : SEGACT MCHAM1
239 : MELP = MCHAM1.IELVAL(1)
240 : SEGDES MCHEL1
241 : SEGDES MCHAM1
242 : C
243 : C**** Vitesse et cosinus directeurs du repere (n,t)
244 : C
245 : MCHEL1 = IVITF
246 : SEGACT MCHEL1
247 : C
248 : C**** La vitesse a comme SPG MELEFE
249 : C Le cosinus directeurs ont comme SPG MELEMF
250 : C
251 : C MCHAM1 -> Cosinus directeurs
252 : C MCHAM2 -> Vitesse
253 : C
254 : SPG1 = MCHEL1.IMACHE(1)
255 : SPG2 = MCHEL1.IMACHE(2)
256 : IF((SPG1 .EQ. MELEMF) .AND. (SPG2 .EQ. MELEFE))THEN
257 : MCHAM1 = MCHEL1.ICHAML(1)
258 : MCHAM2 = MCHEL1.ICHAML(2)
259 : ELSEIF((SPG1 .EQ. MELEFE) .AND. (SPG2 .EQ. MELEMF))THEN
260 : MCHAM1 = MCHEL1.ICHAML(2)
261 : MCHAM2 = MCHEL1.ICHAML(1)
262 : ELSE
263 : LOGAN = .TRUE.
264 : GOTO 9999
265 : ENDIF
266 : SEGACT MCHAM1
267 : MELVNX = MCHAM1.IELVAL(1)
268 : MELVNY = MCHAM1.IELVAL(2)
269 : MELVNZ = MCHAM1.IELVAL(3)
270 : MELT1X = MCHAM1.IELVAL(4)
271 : MELT1Y = MCHAM1.IELVAL(5)
272 : MELT1Z = MCHAM1.IELVAL(6)
273 : MELT2X = MCHAM1.IELVAL(7)
274 : MELT2Y = MCHAM1.IELVAL(8)
275 : MELT2Z = MCHAM1.IELVAL(9)
276 : SEGDES MCHAM1
277 : SEGACT MCHAM2
278 : MELVUN = MCHAM2.IELVAL(1)
279 : MELVUT = MCHAM2.IELVAL(2)
280 : MELVUV = MCHAM2.IELVAL(3)
281 : SEGDES MCHAM2
282 : SEGDES MCHEL1
283 : C
284 : C**** Fractions massiques
285 : C
286 : IF(LOGME)THEN
287 : MCHEL1 = IFRMAF
288 : SEGACT MCHEL1
289 : MCHAMY = MCHEL1.ICHAML(1)
290 : SEGACT MCHAMY
291 : C
292 : C******* Numero d'especes dans les equations d'Euler
293 : C
294 : NESP = MCHAMY.IELVAL(/1)
295 : DO I1 = 1, NESP
296 : MELVA1 = MCHAMY.IELVAL(I1)
297 : SEGACT MELVA1
298 : ENDDO
299 : NMA = NESP
300 : SEGINI FRAMAG
301 : SEGINI FRAMAD
302 : SEGDES MCHEL1
303 : ELSE
304 : C
305 : C******* Definition minimale de YET, necessaire pour transmetre YET aux
306 : C subroutines FORTRAN qui calculent les flux
307 : C
308 : NMA = 1
309 : SEGINI FRAMAG
310 : SEGINI FRAMAD
311 : NESP = 0
312 : ENDIF
313 : C
314 : C**** Scalaires passifs
315 : C
316 : IF(ISCAF .GT. 0)THEN
317 : MCHEL1 = ISCAF
318 : SEGACT MCHEL1
319 : MCHAMS = MCHEL1.ICHAML(1)
320 : SEGACT MCHAMS
321 : NSCA = MCHAMS.IELVAL(/1)
322 : DO I1 = 1, NSCA, 1
323 : MELVA1 = MCHAMS.IELVAL(I1)
324 : SEGACT MELVA1
325 : ENDDO
326 : NMA = NSCA
327 : SEGINI SCAG
328 : SEGINI SCAD
329 : SEGDES MCHEL1
330 : ELSE
331 : C
332 : C******* Definition minimale de YET, necessaire pour transmetre YET aux
333 : C subroutines FORTRAN qui calculent les flux
334 : C
335 : NMA = 1
336 : SEGINI SCAG
337 : SEGINI SCAD
338 : NSCA= 0
339 : ENDIF
340 : C
341 : C
342 : C**** Initialisation des MELEMEs
343 : C
344 : C 'CENTRE', 'FACEL'
345 : C
346 : IPT2 = MELEFE
347 : SEGACT IPT2
348 : NFAC = IPT2.NUM(/2)
349 : C
350 : C**** KRIPAD pour la correspondance global/local de centre
351 : C
351 : CALL KRIPAD(MELEMC,MLENT1)
352 : C
353 : C**** MLENTI1 a MCORD.XCOORD(/1)/(IDIM+1) elements
354 : C
355 : C Si i est le numero global d'un noeud de ICEN,
356 : C MLENT1.LECT(i) contient sa position, i.e.
357 : C
358 : C I = numero global du noeud centre
359 : C MLENT1.LECT(i) = numero local du noeud centre
360 : C
361 : C MLENT1 déjà activé, i.e.
362 : C
363 : C SEGACT MLENT1
364 : C
365 : C
366 : C**** KRIPAD pour la correspondance global/local de 'FACE'
367 : C
367 : CALL KRIPAD(MELEMF,MLENT2)
368 : C
369 : C**** Initialisation de flux
370 : C
371 : NMA2 = 5 + NESP + NSCA
372 : SEGINI IFLU1
373 : SEGINI IFLU2
374 : C
375 : C**** IFLU2 = segment de travail en FLUVLH; c'est plus rapide le definir ici
376 : C
377 : C
378 : C**** CHPOINTs de la table DOMAINE
379 : C
379 : CALL LICHT(ICHPSU,MPOVSU,TYPE,IGEOMF)
379 : CALL LICHT(ICHPDI,MPOVDI,TYPE,IGEOMC)
380 : C
381 : C**** LICHT active les MPOVALs en *MOD
382 : C
383 : C i.e.
384 : C
385 : C SEGACT MPOVSU*MOD
386 : C SEGACT MPOVDI*MOD
387 : C
388 : C
389 : C**** Les FLUX aux face
390 : C
391 : C La densité
392 : C
392 : CALL LICHT(IZG1,MPOVG1,TYPE,IGEOMF)
393 : C
394 : C SEGACT MPOVG1*MOD
395 : C
396 : C**** Les debits
397 : C
397 : CALL LICHT(IZG2,MPOVG2,TYPE,IGEOMF)
398 : C
399 : C SEGACT MPOVG2*MOD
400 : C
401 : C**** L'energie totale volumique
402 : C
402 : CALL LICHT(IZG3,MPOVG3,TYPE,IGEOMF)
403 : C
404 : C SEGACT MPOVG3*MOD
405 : C
406 : C**** Les Fractions Massiques
407 : C
408 : IF(LOGME)THEN
408 : CALL LICHT(IZG4,MPOVG4,TYPE,IGEOMF)
409 : C
410 : C SEGACT MPOVG4*MOD
411 : C
412 : ENDIF
413 : C
414 : C**** Les scalaires passifs
415 : C
416 : IF(ISCAF .GT. 0)THEN
416 : CALL LICHT(IZG5,MPOVG5,TYPE,IGEOMF)
417 : C
418 : C SEGACT MPOVG5*MOD
419 : C
420 : ENDIF
421 : C
422 : C**** Activation des MCHAMLs
423 : C
424 : SEGACT MELRO
425 : SEGACT MELP
426 : SEGACT MELVUN
427 : SEGACT MELVUT
428 : SEGACT MELVUV
429 : SEGACT MELVNX
430 : SEGACT MELVNY
431 : SEGACT MELVNZ
432 : SEGACT MELT1X
433 : SEGACT MELT1Y
434 : SEGACT MELT1Z
435 : SEGACT MELT2X
436 : SEGACT MELT2Y
437 : SEGACT MELT2Z
438 : C
439 : C**** Initialisation de 1/DT
440 : C
441 : UNSDT = 0.0D0
442 : C
443 : C**** BOUCLE SUR FACEL pour le calcul du FLUX
444 : C
445 : DO NLCF = 1, NFAC
446 : C
447 : C******* NLCF = numero local du centre de facel
448 : C NGCF = numero global du centre de facel
449 : C NLCF1 = numero local du centre de face
450 : C NGCEG = numero global du centre ELT "gauche"
451 : C NLCEG = numero local du centre ELT "gauche"
452 : C NGCED = numero global du centre ELT "droite"
453 : C NLCED = numero local du centre ELT "droite"
454 : C
455 : NGCEG = IPT2.NUM(1,NLCF)
456 : NGCED = IPT2.NUM(3,NLCF)
457 : NGCF = IPT2.NUM(2,NLCF)
458 : NLCF1 = MLENT2.LECT(NGCF)
459 : NLCEG = MLENT1.LECT(NGCEG)
460 : NLCED = MLENT1.LECT(NGCED)
461 : C
462 : C******* NLCF != NLCF1 -> l'auteur (MOI) n'a rien compris.
463 : C
464 : IF(NLCF .NE. NLCF1)THEN
465 : MESERR = 'Il ne faut pas jouer avec la console. '
466 : LOGAN = .TRUE.
467 : GOTO 9999
468 : ENDIF
469 : C
470 : C******* Recuperation des Etats "gauche" et "droite"
471 : C
472 : ROG = MELRO.VELCHE(1,NLCF)
473 : UNG = MELVUN.VELCHE(1,NLCF)
474 : UTG = MELVUT.VELCHE(1,NLCF)
475 : UVG = MELVUV.VELCHE(1,NLCF)
476 : PG = MELP.VELCHE(1,NLCF)
477 : C
478 : ROD = MELRO.VELCHE(3,NLCF)
479 : UND = MELVUN.VELCHE(3,NLCF)
480 : UTD = MELVUT.VELCHE(3,NLCF)
481 : UVD = MELVUV.VELCHE(3,NLCF)
482 : PD = MELP.VELCHE(3,NLCF)
483 : C
484 : CNX = MELVNX.VELCHE(1,NLCF)
485 : CNY = MELVNY.VELCHE(1,NLCF)
486 : CNZ = MELVNZ.VELCHE(1,NLCF)
487 : CTX = MELT1X.VELCHE(1,NLCF)
488 : CTY = MELT1Y.VELCHE(1,NLCF)
489 : CTZ = MELT1Z.VELCHE(1,NLCF)
490 : CVX = MELT2X.VELCHE(1,NLCF)
491 : CVY = MELT2Y.VELCHE(1,NLCF)
492 : CVZ = MELT2Z.VELCHE(1,NLCF)
493 : C
494 : C******* Le fractiones massiques, R et ACVTO
495 : C
496 : RTOTG = 0.0D0
497 : RTOTD = 0.0D0
498 : CVTOTG = 0.0D0
499 : CVTOTD = 0.0D0
500 : Y1G = 1.0D0
501 : Y1D = 1.0D0
502 : ETHERG = 0.0D0
503 : ETHERD = 0.0D0
504 : DO I1= 1, NORDP1, 1
505 : PROPHY.ACVTOG(I1) =0.0D0
506 : PROPHY.ACVTOD(I1) =0.0D0
507 : ENDDO
508 : DO I1 = 1, NESP, 1
509 : MELVA1 = MCHAMY.IELVAL(I1)
510 : YG = MELVA1.VELCHE(1,NLCF)
511 : YD = MELVA1.VELCHE(3,NLCF)
512 : Y1G = Y1G - YG
513 : Y1D = Y1D - YD
514 : FRAMAG.YET(I1) = YG
515 : FRAMAD.YET(I1) = YD
516 : RTOTG = RTOTG + YG * PROPHY.R(I1)
517 : RTOTD = RTOTD + YD * PROPHY.R(I1)
518 : DO I2 = 1, NORDP1
519 : PROPHY.ACVTOG(I2) = PROPHY.ACVTOG(I2) +
520 : & YG * PROPHY.ACV(I2,I1)
521 : PROPHY.ACVTOD(I2) = PROPHY.ACVTOD(I2) +
522 : & YD * PROPHY.ACV(I2,I1)
523 : ENDDO
524 : ENDDO
525 : RTOTG = RTOTG + Y1G * PROPHY.R(NESP+1)
526 : RTOTD = RTOTD + Y1D * PROPHY.R(NESP+1)
527 : DO I2 = 1, NORDP1, 1
528 : PROPHY.ACVTOG(I2) = PROPHY.ACVTOG(I2) +
529 : & Y1G * PROPHY.ACV(I2,NESP+1)
530 : PROPHY.ACVTOD(I2) = PROPHY.ACVTOD(I2) +
531 : & Y1D * PROPHY.ACV(I2,NESP+1)
532 : ENDDO
533 : TG = PG / (ROG * RTOTG)
534 : TD = PD / (ROD * RTOTD)
535 : FUNTG = 1.0D0
536 : FUNTD = 1.0D0
537 : CVTOTG = PROPHY.ACVTOG(1)
538 : ETHERG = CVTOTG * TG
539 : CVTOTD = PROPHY.ACVTOD(1)
540 : ETHERD = CVTOTD * TD
541 : DO I1 = 2, NORDP1, 1
542 : FUNTG = FUNTG * TG
543 : FUNTD = FUNTD * TD
544 : DCV = PROPHY.ACVTOG(I1) * FUNTG
545 : CVTOTG = CVTOTG + DCV
546 : ETHERG = ETHERG + DCV * TG / I1
547 : DCV = PROPHY.ACVTOD(I1) * FUNTD
548 : CVTOTD = CVTOTD + DCV
549 : ETHERD = ETHERD + DCV * TD / I1
550 : ENDDO
551 : GAMG = (CVTOTG + RTOTG) /CVTOTG
552 : GAMD = (CVTOTD + RTOTD) /CVTOTD
553 : C
554 : C******* Les scalaires passifs
555 : C
556 : DO I1 = 1, NSCA, 1
557 : MELVA1 = MCHAMS.IELVAL(I1)
558 : SCAG.YET(I1) = MELVA1.VELCHE(1,NLCF)
559 : SCAD.YET(I1) = MELVA1.VELCHE(3,NLCF)
560 : ENDDO
561 : C
562 : C******* On a defini (ROg,ROUNg,ROUTg,Pg,(Yg)), (ROd,ROUNd,ROUTd,Pd,(Yd))
563 : C et on a déjà verifié ROg, ROd, Pg, Pd > 0 et 0FVLHT2(NESP,NSCA,
572 : & GAMG,ROG,PG,UNG,UTG,UVG,ETHERG,
573 : & GAMD,ROD,PD,UND,UTD,UVD,ETHERD,
574 : & FRAMAG.YET,FRAMAD.YET,
575 : & SCAG.YET,SCAD.YET,
576 : & IFLU1.FLUX,IFLU2.FLUX,
577 : & CELLT)
578 : ELSEIF(INDMET .EQ. 2)THEN
579 : C
580 : C******* Colella-Glaz FDS (avec Entropy-fix)
581 : C
582 : XST=0.0D0
582 : CALL FCOLT2(NESP,NSCA,NORDP1,XST,
583 : & GAMG,RTOTG,PROPHY.ACVTOG,ROG,PG,TG,UNG,UTG,UVG,ETHERG,
584 : & GAMD,RTOTD,PROPHY.ACVTOD,ROD,PD,TD,UND,UTD,UVD,ETHERD,
585 : & FRAMAG.YET,FRAMAD.YET,
586 : & SCAG.YET,SCAD.YET,
587 : & IFLU1.FLUX,CELLT,
588 : & LOGNC,MESERR,LOGAN)
589 : ENDIF
590 : C
591 : IF(LOGAN) GOTO 9999
592 : IF(LOGNC) GOTO 9999
593 : C
594 : C******* Ecriture des flux
595 : C
596 : C FLUX(1) = RO Un RO Un
597 : C FLUX(2) = RO Un Un + P -> RO Un Ux + P CNX
598 : C FLUX(3) = RO Un Ut -> RO Un Uy + P CNY
599 : C FLUX(4) = RO Un Et RO Un Et
600 : C
601 : SURF = MPOVSU.VPOCHA(NLCF,1)
602 : MPOVG1.VPOCHA(NLCF,1) = MPOVG1.VPOCHA(NLCF,1) +
603 : & (IFLU1.FLUX(1) * SURF )
604 : MPOVG2.VPOCHA(NLCF,1) = MPOVG2.VPOCHA(NLCF,1) +
605 : &((IFLU1.FLUX(2)*CNX+IFLU1.FLUX(3)*CTX+IFLU1.FLUX(4)*CVX) * SURF)
606 : MPOVG2.VPOCHA(NLCF,2) = MPOVG2.VPOCHA(NLCF,2) +
607 : &((IFLU1.FLUX(2)*CNY+IFLU1.FLUX(3)*CTY+IFLU1.FLUX(4)*CVY) * SURF)
608 : MPOVG2.VPOCHA(NLCF,3) = MPOVG2.VPOCHA(NLCF,3) +
609 : &((IFLU1.FLUX(2)*CNZ+IFLU1.FLUX(3)*CTZ+IFLU1.FLUX(4)*CVZ) * SURF)
610 : MPOVG3.VPOCHA(NLCF,1) = MPOVG3.VPOCHA(NLCF,1) +
611 : & (IFLU1.FLUX(5) * SURF)
612 : DO I1 = 1, NESP, 1
613 : MPOVG4.VPOCHA(NLCF,I1)=IFLU1.FLUX(5+I1)
614 : & * SURF
615 : ENDDO
616 : DO I1 = 1, NSCA, 1
617 : MPOVG5.VPOCHA(NLCF,I1)=IFLU1.FLUX(5+I1+NESP)
618 : & * SURF
619 : ENDDO
620 : C
621 : C******* Calcul du pas du temps (CFL)
622 : C
623 : C****** a) etat a l'interface
624 : C
625 : DIAMG = MPOVDI.VPOCHA(NLCEG,1)
626 : DIAMD = MPOVDI.VPOCHA(NLCED,1)
627 : DIAM = MIN(DIAMG,DIAMD)
628 : CELL = 1.0D0/DIAM/CELLT
629 : IF(CELL .GT. UNSDT)THEN
630 : UNSDT = CELL
631 : DIAMEL = DIAMG
632 : NLCEMI = NLCEG
633 : ENDIF
634 : C
635 : C****** b) etat gauche
636 : C
637 : ASON = SQRT(GAMG*PG/ROG)
638 : LAMBDA = ABS(UNG) + ASON
639 : CELL = LAMBDA / DIAM
640 : IF(CELL .GT. UNSDT)THEN
641 : UNSDT = CELL
642 : DIAMEL = DIAMG
643 : NLCEMI = NLCEG
644 : ENDIF
645 : C
646 : C****** C) etat droite
647 : C
648 : ASON = SQRT(GAMD*PD/ROD)
649 : LAMBDA = ABS(UND) + ASON
650 : CELL = LAMBDA / DIAM
651 : IF(CELL .GT. UNSDT)THEN
652 : UNSDT = CELL
653 : DIAMEL = DIAMD
654 : NLCEMI = NLCED
655 : ENDIF
656 : C
657 : C
658 : C**** Fin boucle sur FACEL
659 : C
660 : ENDDO
661 : C
662 : C**** Pas du temps (condition de non interaction en 1D)
663 : C
664 : DT = 0.5D0 / UNSDT
665 : C
666 : C**** Desactivation des segments et
667 : C on detruit les MCHAMLs
668 : C
669 : C
670 : C**** SEGSUP FRAMAG
671 : C SEGSUP FRAMAD
672 : C
673 : C meme si LOGME = .FALSE.
674 : C
675 : SEGSUP FRAMAG
676 : SEGSUP FRAMAD
677 : C
678 : SEGSUP MLENT1
679 : SEGDES MLENT2
680 : SEGDES IPT2
681 : C
682 : SEGSUP IFLU1
683 : SEGSUP IFLU2
684 : C
685 : SEGDES MPOVSU
686 : SEGDES MPOVDI
687 : C
688 : SEGDES MPOVG1
689 : SEGDES MPOVG2
690 : SEGDES MPOVG3
691 : C
692 : SEGDES MELRO
693 : SEGDES MELP
694 : SEGDES MELVUN
695 : SEGDES MELVUT
696 : SEGDES MELVNX
697 : SEGDES MELVNY
698 : SEGDES MELT1X
699 : SEGDES MELT1Y
700 : C
701 : IF(LOGME) THEN
702 : DO I1 = 1, NESP
703 : MELVA1 = MCHAMY.IELVAL(I1)
704 : SEGDES MELVA1
705 : ENDDO
706 : SEGDES MPOVG4
707 : C
708 : SEGDES MCHAMY
709 : ENDIF
710 : IF(ISCAF .GT. 0) THEN
711 : DO I1 = 1, NSCA, 1
712 : MELVA1 = MCHAMS.IELVAL(I1)
713 : SEGDES MELVA1
714 : ENDDO
715 : SEGDES MPOVG5
716 : C
717 : SEGDES MCHAMS
718 : ENDIF
719 : C
720 : 9999 CONTINUE
721 : C
722 : RETURN
723 : END
724 : C
725 :
726 :
727 :
728 :
729 :
730 :
731 :
732 :
733 :
734 :
735 :
736 :
737 :
738 :
© Cast3M 2003 - All rights reserved.
Disclaimer