1 : C PRE52F SOURCE PV 09/03/12 21:31:10 6325
2 : SUBROUTINE PRE52F(LOGTEM,
3 : & ICEN,IFACE,IFACEL,INORM,
4 : & IALPH, IGRALP, IALALP,
5 : & IUVC, IGRUVC, IALUVC,
6 : & IULC, IGRULC, IALULC,
7 : & IPC, IGRPC, IALPC,
8 : & ITVC, IGRTVC, IALTVC,
9 : & ITLC, IGRTLC, IALTLC,
10 : & IRVC, IRLC,
11 : & DELTAT,
12 : & IALPHF, IUVF, IULF, IPF, ITVF, ITLF,
13 : & IRVF, IRLF,
14 : & LOGAN,LOGNEG,LOGBOR,MESERR,VALER,VAL1,VAL2)
15 : C************************************************************************
16 : C
17 : C PROJET : CASTEM 2000
18 : C
19 : C NOM : PRE121
20 : C
21 : C DESCRIPTION : Voir PRE42F
22 : C
23 : C LANGAGE : FORTRAN 77 + ESOPE 2000 (avec extensions CISI)
24 : C
25 : C AUTEUR : A. BECCANTINI, DRN/DMT/SEMT/TTMF
26 : C Modified for two-fluid flow by
27 : C Jose R. Garcia-Cascales
28 : C
29 : C************************************************************************
30 : C
31 : C APPELES (Outils) : KRIPAD, LICHT
32 : C
33 : C APPELES (Calcul) : AUCUN
34 : C
35 : C************************************************************************
36 : C
37 : C ENTREES
38 : C
39 : C LOGTEM : LOGICAL; si .TRUE. 2em ordre en temps
40 : C sinon 1er ordre en temps
41 : C
42 : C 1) Pointeurs de MELEMEs et de CHPOINTs de la table DOMAINE
43 : C
44 : C ICEN : MELEME de 'POI1' SPG des CENTRES
45 : C
46 : C IFACE : MELEME de 'POI1' SPG des FACES
47 : C
48 : C IFACEL : MELEME de 'SEG3' avec
49 : C CENTRE d'Elt "gauche"
50 : C CENTRE de Face
51 : C CENTRE d'Elt "droite"
52 : C
53 : C N.B. = IFACE.NUM(i,1) = IFACEL.NUM(i,2)
54 : C
55 : C INORM : CHPOINT des cosinus directeurs de normales aux faces
56 : C
57 : C 2) Pointeurs des CHPOINTs
58 : C
59 : C IALPH : CHPOINT "CENTRE" containing void fraction
60 : C
61 : C IGRALP : CHPOINT "CENTRE" contaninig the void fraction gradient
62 : C (2 composantes)
63 : C
64 : C IALALP : CHPOINT "CENTRE" contaning the LIMITA of the
65 : C void fraction gradient
66 : C
67 : C IUVC : CHPOINT "CENTRE" containing the vapour velocity UvX, UVY ;
68 : C
69 : C IGRUVC : CHPOINT "CENTRE" containig the vapour velocity gradient
70 : C (4 composantes)
71 : C
72 : C IALUVC : CHPOINT "CENTRE" contaning the LIMITA of the vapour
73 : C velocity (2 composantes)
74 : C
75 : C IULC : CHPOINT "CENTRE" containing the liquid velocity UvX, UVY ;
76 : C
77 : C IGRULC : CHPOINT "CENTRE" containig the liquid velocity gradient
78 : C (4 composantes)
79 : C
80 : C IALULC : CHPOINT "CENTRE" contaning the LIMITA of the liquid
81 : C velocity (2 composantes)
82 : C
83 : C IPC : CHPOINT "CENTRE" containing the pressure P ;
84 : C
85 : C IGRPC : CHPOINT "CENTRE" containing the gradient of pressure
86 : C (2 composantes)
87 : C
88 : C IALPC : CHPOINT "CENTRE" containig the LIMITA of pressure
89 : C gradient
90 : C
91 : C ITVC : CHPOINT "CENTRE" containing the vapour temperature Tv ;
92 : C
93 : C IGRTVC : CHPOINT "CENTRE" containing the gradient of vapour
94 : C temperature (2 composantes)
95 : C
96 : C IALTVC : CHPOINT "CENTRE" containig the LIMITA of vapour
97 : C temperature gradient
98 : C
99 : C ITLC : CHPOINT "CENTRE" containing the liquid temperature Tl ;
100 : C
101 : C IGRTLC : CHPOINT "CENTRE" containing the gradient of liquid
102 : C temperature (2 composantes)
103 : C
104 : C IALTLC : CHPOINT "CENTRE" containig the LIMITA of liquid
105 : C temperature gradient
106 : C
107 : C IRVC : CHPOINT "CENTRE" containing the vapour density
108 : C
109 : C IRLC : CHPOINT "CENTRE" containing the liquid density
110 : C
111 : C
112 : C DELTAT : REAL*8, encrement en temps pour calculer la prediction
113 : C
114 : C SORTIES
115 : C
116 : C C SORTIES
117 : C
118 : C IALPHF : MCHAML defined en la MELEME of pointers
119 : C IFACEL, it contains the void fraction
120 : C (a gauche et a droite de chaque face).
121 : C Only one component ('SCAL')
122 : C
123 : C IUVF : MCHAML "FACEL" vapour velocity and the
124 : C director cosines (n,t) in the corresponding face;
125 : C in the 2D case 6 composantes:
126 : C 'UVN' = normal velocity
127 : C 'UVT' = tangent velocity,
128 : C these two velocities are defined in a local
129 : C reference framework defined over the MELEME
130 : C of pointers IFACE
131 : C 'NX' = n.x
132 : C 'NY' = n.y
133 : C 'TX' = t.x
134 : C 'TY' = t.y
135 : C
136 : C IULF : MCHAML "FACEL" liquid velocity
137 : C in the 2D case 2 composantes:
138 : C 'ULN' = normal velocity
139 : C 'ULT' = tangent velocity
140 : C
141 : C IPF : MCHAML "FACEL" pressure
142 : C Only one component ('SCAL')
143 : C
144 : C ITVF : MCHAML "FACEL" vapour temperature
145 : C Only one component ('SCAL')
146 : C
147 : C ITVL : MCHAML "FACEL" liquid temperature
148 : C Only one component ('SCAL')
149 : C
150 : C IRVF : MCHAML "FACEL" vapour density
151 : C Only one component ('SCAL')
152 : C
153 : C IRLF : MCHAML "FACEL" liquid density
154 : C Only one component ('SCAL')
155 : C
156 : C LOGAN : anomalie detectee (changement de la convention dans
157 : C la table domaine)
158 : C
159 : C LOGNEG : (LOGICAL): si .TRUE. une pression ou une densité
160 : C negative a été detectée -> en interactif le
161 : C programme s'arrete en GIBIANE
162 : C (erreur stocké en MESERR et VALER)
163 : C
164 : C LOGBOR : (LOGICAL): si .TRUE. un gamma a ete detecte
165 : C dehor 1 et 3 (sa valeur stockée en MESERR et VALER;
166 : C en VAL1 et en VAL2 on stocke 1.0 et 3.0)
167 : C
168 : C MESERR
169 : C VALER
170 : C VAL1,
171 : C VAL2 : pour les messages d'erreur
172 : C
173 : C************************************************************************
174 : C
175 : C HISTORIQUE (Anomalies et modifications éventuelles)
176 : C
177 : C HISTORIQUE : Créée le 11.6.98.
178 : C
179 : C************************************************************************
180 : C
181 : C
182 : C ATTENTION: Cet programme marche si le MAILLAGE est convex;
183 : C si non il faut changer l'argoritme de calcul de
184 : C l'orientation des normales aux faces.
185 : C
186 : C La positivité n'est pas controlle parce que c'est déjà fait
187 : C dans l'operateur PRIM
188 : C
189 : C
190 : C************************************************************************
191 : C
192 : C**** Variables de COOPTIO
193 : C
194 : INTEGER IPLLB, IERPER, IERMAX, IERR, INTERR
195 : & ,IOTER, IOLEC, IOIMP, IOCAR, IOACQ
196 : & ,IOPER, IOSGB, IOGRA, IOSAU, IORES
197 : & ,IECHO, IIMPI, IOSPI
198 : & ,IDIM, IFICLE, IPREFI
199 : C & ,MCOORD
200 : & ,IFOMOD, NIFOUR, IFOUR, NSDPGE, IONIVE
201 : & ,NGMAXY, IZROSF, ISOTYP, IOSCR,LTEXLU
202 : & ,NORINC,NORVAL,NORIND,NORVAD
203 : & ,NUCROU, IPSAUV
204 : C
205 : C**** Les variables
206 : C
207 : INTEGER ICEN, IFACE, IFACEL, INORM,
208 : & IALPH, IGRALP, IALALP,
209 : & IUVC, IGRUVC, IALUVC,
210 : & IULC, IGRULC, IALULC,
211 : & IPC, IGRPC, IALPC,
212 : & ITVC, IGRTVC, IALTVC,
213 : & ITLC, IGRTLC, IALTLC,
214 : & IRVC, IRLC,
215 : & IALPHF, IUVF, IULF, IPF, ITVF, ITLF,
216 : & IRVF, IRLF,
217 : & IGEOM, NFAC, NCEN,
218 : & N1PTEL, N1EL, N2PTEL, N2EL, N2, N1, N3, L1, NLCE,
219 : & NLCF, NGCF, NGCEG, NLCEG, NGCED, NLCED, NGCF1,
220 : & IDIMP1, INDCEL,
221 : & l, m, n
222 : REAL*8 VALER, VAL1, VAL2, XG, YG, XC, YC, XD, YD, DELTAT
223 : & , DXG, DYG, DXD, DYD
224 : & , CNX, CNY, CTX, CTY, ORIENT
225 : & , alphag, uvxg, uvyg, uvng, uvtg, ulxg, ulyg, ulng, ultg
226 : & , pg, Tvg, Tlg, rvg, rlg
227 : & , alphad, uvxd, uvyd, uvnd, uvtd, ulxd, ulyd, ulnd, ultd
228 : & , pd, Tvd, Tld, rvd, rld
229 : & , VALCEL, DCEL, ALCEL
230 : & , alpha, uvx, uvy, ulx, uly, p, Tv, Tl
231 : & , rv, rl, drvdp, drvdhv, drldp, drldhl
232 : & , hv, hl, cv, cl
233 : & , A(8,8), B(8,8), denom
234 : & , dVdx(8), dVdy(8), dV(8), LIMITA(8), dVdVp(8,8)
235 : & , dVpdV(8,8), P1(8,8), P2(8,8)
236 : & , dhvdp, dhvdTv, dhldp, dhldTl
237 : CHARACTER*(40) MESERR
238 : CHARACTER*(8) TYPE
239 : LOGICAL LOGAN,LOGNEG, LOGBOR, LOGTEM, LOGI1, LOGALP
240 : C
241 : C**** perfect and stiffened gas eos
242 : C
243 : real*8 gamv, cpv, gaml, cpl, pil
244 : PARAMETER(gamv = 1.4D0, cpv = 1008.D0,
245 : & gaml = 2.8D0, cpl = 4186.D0, pil = 8.5D8)
246 : C
247 : C**** LOGALP = .TRUE. -> Prediction avec limiteur
248 : C
249 : PARAMETER(LOGALP = .TRUE.)
250 : C
251 : C**** Les Includes
252 : C
253 : -INC SMCOORD
254 : -INC CCOPTIO
255 : -INC SMCHPOI
256 : POINTEUR MPALPC.MPOVAL, MPGRAL.MPOVAL,
257 : & MPUVC.MPOVAL, MPGRUV.MPOVAL,
258 : & MPULC.MPOVAL, MPGRUL.MPOVAL,
259 : & MPPC.MPOVAL, MPGRP.MPOVAL,
260 : & MPTVC.MPOVAL, MPGRTV.MPOVAL,
261 : & MPTLC.MPOVAL, MPGRTL.MPOVAL,
262 : & MPRVC.MPOVAL, MPRLC.MPOVAL,
263 : & MPNORM.MPOVAL,
264 : & MPALPP.MPOVAL, MPUVP.MPOVAL, MPULP.MPOVAL,
265 : & MPPP.MPOVAL, MPTVP.MPOVAL, MPTLP.MPOVAL,
266 : & MPRVP.MPOVAL, MPRLP.MPOVAL
267 : -INC SMCHAML
268 : POINTEUR MLALP.MELVAL, MLP.MELVAL,
269 : & MLTV.MELVAL, MLTL.MELVAL,
270 : & MLRV.MELVAL, MLRL.MELVAL
271 : POINTEUR MLUVN.MELVAL, MLUVT.MELVAL,
272 : & MLULN.MELVAL, MLULT.MELVAL,
273 : & MLVNX.MELVAL, MLVNY.MELVAL, MLVTX.MELVAL, MLVTY.MELVAL,
274 : & MLLNX.MELVAL, MLLNY.MELVAL, MLLTX.MELVAL, MLLTY.MELVAL
275 : C**** Some melemes for the LIMITAs
276 : POINTEUR MELLAL.MPOVAL,
277 : & MELLUV.MPOVAL,
278 : & MELLUL.MPOVAL,
279 : & MELLP.MPOVAL,
280 : & MELLTV.MPOVAL,
281 : & MELLTL.MPOVAL
282 : -INC SMLENTI
283 : -INC SMELEME
284 : C
285 : C
286 : C**** Initialisation des parametres d'erreur déjà faite, i.e.
287 : C
288 : C LOGNEG = .FALSE.
289 : C LOGBOR = .FALSE.
290 : C MESERR = ' '
291 : C MOTERR(1:40) = MESERR(1:40)
292 : C VALER = 0.0D0
293 : C VAL1 = 0.0D0
294 : C VAL2 = 0.0D0
295 : C
296 : C
297 : C**** KRIPAD pour la correspondance global/local de centre
298 : C
298 : CALL KRIPAD(ICEN,MLENT1)
299 : C
300 : C**** MLENTI1 a MCORD.XCOORD(/1)/(IDIM+1) elements
301 : C
302 : C Si i est le numero global d'un noeud de ICEN,
303 : C MLENT1.LECT(i) contient sa position, i.e.
304 : C
305 : C I = numero global du noeud centre
306 : C MLENT1.LECT(i) = numero local du noeud centre
307 : C
308 : C MLENT1 déjà activé, i.e.
309 : C
310 : C SEGACT MLENT1
311 : C
312 : C**** Activation de CHPOINTs
313 : C
314 : C void fraction + grad + limiter
315 : C vapour velocity + grad + limiter
316 : C liquid velocity + grad + limiter
317 : C pressure + grad + limiter
318 : C vapour temperature + grad + limiter
319 : C liquid temperature + grad + limiter
320 : C vapour density
321 : C liquid density
322 : C
323 : C cosinus directeurs des normales aux surface
324 : C
324 : CALL LICHT(IALPH, MPALPC, TYPE, IGEOM)
324 : CALL LICHT(IGRALP, MPGRAL, TYPE, IGEOM)
324 : CALL LICHT(IUVC, MPUVC, TYPE, IGEOM)
324 : CALL LICHT(IGRUVC, MPGRUV, TYPE, IGEOM)
324 : CALL LICHT(IULC, MPULC, TYPE, IGEOM)
324 : CALL LICHT(IGRULC, MPGRUL, TYPE, IGEOM)
324 : CALL LICHT(IPC , MPPC , TYPE, IGEOM)
324 : CALL LICHT(IGRPC, MPGRP, TYPE, IGEOM)
324 : CALL LICHT(ITVC, MPTVC, TYPE, IGEOM)
324 : CALL LICHT(IGRTVC, MPGRTV, TYPE, IGEOM)
324 : CALL LICHT(ITLC, MPTLC, TYPE, IGEOM)
324 : CALL LICHT(IGRTLC, MPGRTL, TYPE, IGEOM)
324 : CALL LICHT(IRVC, MPRVC, TYPE, IGEOM)
324 : CALL LICHT(IRLC, MPRLC, TYPE, IGEOM)
324 : CALL LICHT(INORM, MPNORM, TYPE, IGEOM)
325 : C
326 : C**** Les MPOVALs 'Prediction'
327 : C
328 : IF(LOGTEM)THEN
329 : SEGINI, MPALPP = MPALPC
330 : SEGINI, MPUVP = MPUVC
331 : SEGINI, MPULP = MPULC
332 : SEGINI, MPPP = MPPC
333 : SEGINI, MPTVP = MPTVC
334 : SEGINI, MPTLP = MPTLC
335 : SEGINI, MPRVP = MPRVC
336 : SEGINI, MPRLP = MPRLC
337 : ELSE
338 : MPALPP = MPALPC
339 : MPUVP = MPUVC
340 : MPULP = MPULC
341 : MPPP = MPPC
342 : MPTVP = MPTVC
343 : MPTLP = MPTLC
344 : MPRVP = MPRVC
345 : MPRLP = MPRLC
346 : ENDIF
347 : C
348 : C**** Les Limiteurs
349 : C
349 : CALL LICHT(IALALP, MELLAL, TYPE, IGEOM)
349 : CALL LICHT(IALUVC, MELLUV, TYPE, IGEOM)
349 : CALL LICHT(IALULC, MELLUL, TYPE, IGEOM)
349 : CALL LICHT(IALPC , MELLP, TYPE, IGEOM)
349 : CALL LICHT(IALTVC, MELLTV, TYPE, IGEOM)
349 : CALL LICHT(IALTLC, MELLTL, TYPE, IGEOM)
350 : C
351 : C
352 : C**** MPOVAL sont déjà activés i.e.:
353 : C
354 : C SEGACT MPALPP
355 : C SEGACT MPGRAL
356 : C SEGACT IALALP
357 : C SEGACT MPUVC
358 : C SEGACT IALUVC
359 : C SEGACT MPGRUV
360 : C SEGACT MPULC
361 : C SEGACT MPGRUL
362 : C SEGACT IALULC
363 : C SEGACT MPPC
364 : C SEGACT MPGRP
365 : C SEGACT IALPC
366 : C SEGACT MPTVC
367 : C SEGACT MPGRTV
368 : C SEGACT IALTVC
369 : C SEGACT MPTLC
370 : C SEGACT MPGRTL
371 : C SEGACT IALTLC
372 : C SEGACT MPRVC
373 : C SEGACT MPRLC
374 : C SEGACT MPNORM
375 : C
376 : C**** Le MELEME FACEL
377 : C
378 : IPT1 = IFACEL
379 : IPT2 = IFACE
380 : SEGACT IPT1
381 : SEGACT IPT2
382 : NFAC = IPT1.NUM(/2)
383 : C
384 : C**** Creation de MCHAMLs contenant les etats gauche et droite,
385 : C
386 : C i.e.:
387 : C
388 : C vitesse + cosinus directors du repere local
389 : C densité
390 : C pression
391 : C gamma
392 : C
393 : C**** Cosinus directors du repere local et vitesse
394 : C
395 : C Les cosinus directeurs
396 : C
397 : C VAPOUR PHASE
398 :
399 : N1 = 2
400 : N3 = 6
401 : L1 = 28
402 : SEGINI MCHEL1
403 : IUVF = MCHEL1
404 : MCHEL1.TITCHE = 'UV '
405 : MCHEL1.IMACHE(1) = IFACE
406 : MCHEL1.IMACHE(2) = IFACEL
407 : MCHEL1.CONCHE(1) = ' (n,t) in (x,y) '
408 : MCHEL1.CONCHE(2) = ' UV in (n,t) '
409 : C
410 : C**** Valeurs des cosinus definies par respect au repair global, i.e.
411 : C
412 : MCHEL1.INFCHE(1,1) = 2
413 : MCHEL1.INFCHE(1,3) = NIFOUR
414 : MCHEL1.INFCHE(1,4) = 0
415 : MCHEL1.INFCHE(1,5) = 0
416 : MCHEL1.INFCHE(1,6) = 0
417 : MCHEL1.IFOCHE = IFOUR
418 : C
419 : C**** Valeurs de vitesse definies par respect au repair local, i.e.
420 : C
421 : MCHEL1.INFCHE(2,1) = 1
422 : MCHEL1.INFCHE(2,3) = NIFOUR
423 : MCHEL1.INFCHE(2,4) = 0
424 : MCHEL1.INFCHE(2,5) = 0
425 : MCHEL1.INFCHE(2,6) = 0
426 : C
427 : C**** Le cosinus directeurs
428 : C
429 : N1PTEL = 1
430 : N1EL = NFAC
431 : N2PTEL = 0
432 : N2EL = 0
433 : C
434 : C**** MCHAML a N2 composantes:
435 : C
436 : C cosinus directeurs du repere local (n,t1)
437 : C
438 : C IDIM = 2 -> 4 composantes
439 : C
440 : N2 = 4
441 : SEGINI MCHAM1
442 : MCHEL1.ICHAML(1) = MCHAM1
443 : MCHAM1.NOMCHE(1) = 'NVX '
444 : MCHAM1.NOMCHE(2) = 'NVY '
445 : MCHAM1.NOMCHE(3) = 'TVX '
446 : MCHAM1.NOMCHE(4) = 'TVY '
447 : MCHAM1.TYPCHE(1) = 'REAL*8 '
448 : MCHAM1.TYPCHE(2) = 'REAL*8 '
449 : MCHAM1.TYPCHE(3) = 'REAL*8 '
450 : MCHAM1.TYPCHE(4) = 'REAL*8 '
451 : SEGINI MLVNX
452 : SEGINI MLVNY
453 : SEGINI MLVTX
454 : SEGINI MLVTY
455 : MCHAM1.IELVAL(1) = MLVNX
456 : MCHAM1.IELVAL(2) = MLVNY
457 : MCHAM1.IELVAL(3) = MLVTX
458 : MCHAM1.IELVAL(4) = MLVTY
459 : SEGDES MCHAM1
460 : C
461 : C**** Vitesse
462 : C
463 : N1EL = NFAC
464 : N1PTEL = 3
465 : N2EL = 0
466 : N2PTEL = 0
467 : C
468 : C**** MCHAML a N2 composantes:
469 : C
470 : C IDIM = 2 -> 2 composantes
471 : C
472 : N2 = 2
473 : SEGINI MCHAM1
474 : MCHEL1.ICHAML(2) = MCHAM1
475 : SEGDES MCHEL1
476 : MCHAM1.NOMCHE(1) = 'UVN '
477 : MCHAM1.NOMCHE(2) = 'UVT '
478 : MCHAM1.TYPCHE(1) = 'REAL*8 '
479 : MCHAM1.TYPCHE(2) = 'REAL*8 '
480 : SEGINI MLUVN
481 : SEGINI MLUVT
482 : MCHAM1.IELVAL(1) = MLUVN
483 : MCHAM1.IELVAL(2) = MLUVT
484 : SEGDES MCHAM1
485 : C
486 : C**** Cosinus directors du repere local et vitesse
487 : C
488 : C Les cosinus directeurs
489 : C
490 : C LIQUID PHASE
491 :
492 : N1 = 2
493 : N3 = 6
494 : L1 = 28
495 : SEGINI MCHEL1
496 : IULF = MCHEL1
497 : MCHEL1.TITCHE = 'UL '
498 : MCHEL1.IMACHE(1) = IFACE
499 : MCHEL1.IMACHE(2) = IFACEL
500 : MCHEL1.CONCHE(1) = ' (n,t) in (x,y) '
501 : MCHEL1.CONCHE(2) = ' UL in (n,t) '
502 : C
503 : C**** Valeurs des cosinus definies par respect au repair global, i.e.
504 : C
505 : MCHEL1.INFCHE(1,1) = 2
506 : MCHEL1.INFCHE(1,3) = NIFOUR
507 : MCHEL1.INFCHE(1,4) = 0
508 : MCHEL1.INFCHE(1,5) = 0
509 : MCHEL1.INFCHE(1,6) = 0
510 : MCHEL1.IFOCHE = IFOUR
511 : C
512 : C**** Valeurs de vitesse definies par respect au repair local, i.e.
513 : C
514 : MCHEL1.INFCHE(2,1) = 1
515 : MCHEL1.INFCHE(2,3) = NIFOUR
516 : MCHEL1.INFCHE(2,4) = 0
517 : MCHEL1.INFCHE(2,5) = 0
518 : MCHEL1.INFCHE(2,6) = 0
519 : C
520 : C**** Le cosinus directeurs
521 : C
522 : N1PTEL = 1
523 : N1EL = NFAC
524 : N2PTEL = 0
525 : N2EL = 0
526 : C
527 : C**** MCHAML a N2 composantes:
528 : C
529 : C cosinus directeurs du repere local (n,t1)
530 : C
531 : C IDIM = 2 -> 4 composantes
532 : C
533 : N2 = 4
534 : SEGINI MCHAM1
535 : MCHEL1.ICHAML(1) = MCHAM1
536 : MCHAM1.NOMCHE(1) = 'NLX '
537 : MCHAM1.NOMCHE(2) = 'NLY '
538 : MCHAM1.NOMCHE(3) = 'TLX '
539 : MCHAM1.NOMCHE(4) = 'TLY '
540 : MCHAM1.TYPCHE(1) = 'REAL*8 '
541 : MCHAM1.TYPCHE(2) = 'REAL*8 '
542 : MCHAM1.TYPCHE(3) = 'REAL*8 '
543 : MCHAM1.TYPCHE(4) = 'REAL*8 '
544 : SEGINI MLLNX
545 : SEGINI MLLNY
546 : SEGINI MLLTX
547 : SEGINI MLLTY
548 : MCHAM1.IELVAL(1) = MLLNX
549 : MCHAM1.IELVAL(2) = MLLNY
550 : MCHAM1.IELVAL(3) = MLLTX
551 : MCHAM1.IELVAL(4) = MLLTY
552 : SEGDES MCHAM1
553 : C
554 : C**** Vitesse
555 : C
556 : N1EL = NFAC
557 : N1PTEL = 3
558 : N2EL = 0
559 : N2PTEL = 0
560 : C
561 : C**** MCHAML a N2 composantes:
562 : C
563 : C IDIM = 2 -> 2 composantes
564 : C
565 : N2 = 2
566 : SEGINI MCHAM1
567 : MCHEL1.ICHAML(2) = MCHAM1
568 : SEGDES MCHEL1
569 : MCHAM1.NOMCHE(1) = 'ULN '
570 : MCHAM1.NOMCHE(2) = 'ULT '
571 : MCHAM1.TYPCHE(1) = 'REAL*8 '
572 : MCHAM1.TYPCHE(2) = 'REAL*8 '
573 : SEGINI MLULN
574 : SEGINI MLULT
575 : MCHAM1.IELVAL(1) = MLULN
576 : MCHAM1.IELVAL(2) = MLULT
577 : SEGDES MCHAM1
578 : C
579 : C**** Void fraction
580 : C
581 : N1 = 1
582 : N3 = 6
583 : L1 = 15
584 : SEGINI MCHEL2
585 : IALPHF = MCHEL2
586 : MCHEL2.IMACHE(1) = IFACEL
587 : MCHEL2.TITCHE = 'ALPHA '
588 : MCHEL2.CONCHE(1) = ' '
589 : C
590 : C**** Valeurs independente du repére, i.e.
591 : C
592 : MCHEL2.INFCHE(1,1) = 0
593 : MCHEL2.INFCHE(1,3) = NIFOUR
594 : MCHEL2.INFCHE(1,4) = 0
595 : MCHEL2.INFCHE(1,5) = 0
596 : MCHEL2.INFCHE(1,6) = 0
597 : MCHEL2.IFOCHE = IFOUR
598 : N2 = 1
599 : SEGINI MCHAM1
600 : MCHEL2.ICHAML(1) = MCHAM1
601 : SEGDES MCHEL2
602 : MCHAM1.NOMCHE(1) = 'SCAL '
603 : MCHAM1.TYPCHE(1) = 'REAL*8 '
604 : SEGINI MLALP
605 : MCHAM1.IELVAL(1) = MLALP
606 : SEGDES MCHAM1
607 : C
608 : C**** Pressure
609 : C
610 : MCHEL1 = IALPHF
611 : SEGINI, MCHEL2 = MCHEL1
612 : IPF = MCHEL2
613 : MCHEL2.TITCHE = 'P '
614 : C
615 : C**** MCHAM1 = MCHAML de la alpha
616 : C
617 : SEGINI, MCHAM2 = MCHAM1
618 : MCHEL2.ICHAML(1) = MCHAM2
619 : SEGDES MCHEL2
620 : SEGINI MLP
621 : MCHAM2.IELVAL(1) = MLP
622 : SEGDES MCHAM2
623 : C
624 : C**** Vapour temperature
625 : C
626 : MCHEL1 = IALPHF
627 : SEGINI, MCHEL2 = MCHEL1
628 : ITVF = MCHEL2
629 : MCHEL2.TITCHE = 'TV '
630 : C
631 : C**** MCHAM1 = MCHAML de la alpha
632 : C
633 : SEGINI, MCHAM2 = MCHAM1
634 : MCHEL2.ICHAML(1) = MCHAM2
635 : SEGDES MCHEL2
636 : SEGINI MLTV
637 : MCHAM2.IELVAL(1) = MLTV
638 : SEGDES MCHAM2
639 : C
640 : C**** Liquid temperature
641 : C
642 : MCHEL1 = IALPHF
643 : SEGINI, MCHEL2 = MCHEL1
644 : ITLF = MCHEL2
645 : MCHEL2.TITCHE = 'TL '
646 : C
647 : C**** MCHAM1 = MCHAML de la alpha
648 : C
649 : SEGINI, MCHAM2 = MCHAM1
650 : MCHEL2.ICHAML(1) = MCHAM2
651 : SEGDES MCHEL2
652 : SEGINI MLTL
653 : MCHAM2.IELVAL(1) = MLTL
654 : SEGDES MCHAM2
655 : C
656 : C**** Vapour density
657 : C
658 : MCHEL1 = IALPHF
659 : SEGINI, MCHEL2 = MCHEL1
660 : IRVF = MCHEL2
661 : MCHEL2.TITCHE = 'RV '
662 : C
663 : C**** MCHAM1 = MCHAML de la alpha
664 : C
665 : SEGINI, MCHAM2 = MCHAM1
666 : MCHEL2.ICHAML(1) = MCHAM2
667 : SEGDES MCHEL2
668 : SEGINI MLRV
669 : MCHAM2.IELVAL(1) = MLRV
670 : SEGDES MCHAM2
671 : C
672 : C**** Liquid density
673 : C
674 : MCHEL1 = IALPHF
675 : SEGINI, MCHEL2 = MCHEL1
676 : IRLF = MCHEL2
677 : MCHEL2.TITCHE = 'RL '
678 : C
679 : C**** MCHAM1 = MCHAML de la alpha
680 : C
681 : SEGINI, MCHAM2 = MCHAM1
682 : MCHEL2.ICHAML(1) = MCHAM2
683 : SEGDES MCHEL2
684 : SEGINI MLRL
685 : MCHAM2.IELVAL(1) = MLRL
686 : SEGDES MCHAM2
687 :
688 : C
689 : C**** Recapitulatif: les MELVALs et les MPOVALs actives
690 : C
691 : C MLVNX, MLVNY,
692 : C MLVTX, MLVTY,
693 : C
694 : C MLLNX, MLLNY,
695 : C MLLTX, MLLTY
696 : C
697 : C MLUVN, MLUVT -> vapour velocities
698 : C
699 : C MLULN, MLULT -> liquid velocities
700 : C
701 : C MLALP -> void fraction
702 : C
703 : C MLP -> pressure
704 : C
705 : C MLTV -> vapour temperature
706 : C
707 : C MLTL -> liquid temperature
708 : C
709 : C MLRV -> vapour density
710 : C
711 : C MLRL -> liquid density
712 : C****
713 : C MPALPC -> void fraction
714 : C
715 : C MPUVC -> vapour velocity
716 : C
717 : C MPULC -> liquid velocity
718 : C
719 : C MPPC -> pressure
720 : C
721 : C MPTVC -> vapour temperature
722 : C
723 : C MPTLC -> liquid temperature
724 : C
725 : C MPRVC -> vapour density
726 : C
727 : C MPRLC -> liquid density
728 : C
729 : C MPNORM -> normales aux faces
730 : C
731 : C***********************************************************************
732 : C********* PREDICTION **************************************************
733 : C***********************************************************************
734 : C
735 : Cµµµµµµµµµµµµµµµµµµµµµ AQUI ME QUEDO
736 : C
737 : IF(LOGTEM)THEN
738 : C
739 : IPT3 = ICEN
740 : SEGACT IPT3
741 : NCEN = IPT3.NUM(/2)
742 : DO NLCE = 1, NCEN
743 : IF(LOGALP)THEN
744 : C
745 : C************* Gradients limités
746 : C
747 : LIMITA(1) = MELLAL.VPOCHA(NLCE,1)
748 : LIMITA(2) = MELLUV.VPOCHA(NLCE,1)
749 : LIMITA(3) = MELLUV.VPOCHA(NLCE,2)
750 : LIMITA(4) = MELLUL.VPOCHA(NLCE,1)
751 : LIMITA(5) = MELLUL.VPOCHA(NLCE,2)
752 : LIMITA(6) = MELLP.VPOCHA(NLCE,1)
753 : LIMITA(7) = MELLTV.VPOCHA(NLCE,1)
754 : LIMITA(8) = MELLTL.VPOCHA(NLCE,1)
755 : ELSE
756 : C
757 : C************* Gradients non limités
758 : C
759 : LIMITA(1) = 1.0D0
760 : LIMITA(2) = 1.0D0
761 : LIMITA(3) = 1.0D0
762 : LIMITA(4) = 1.0D0
763 : LIMITA(5) = 1.0D0
764 : LIMITA(6) = 1.0D0
765 : LIMITA(7) = 1.0D0
766 : LIMITA(8) = 1.0D0
767 : ENDIF
768 :
769 : alpha = MPALPP.VPOCHA(NLCE,1)
770 : uvx = MPUVP.VPOCHA(NLCE,1)
771 : uvy = MPUVP.VPOCHA(NLCE,2)
772 : ulx = MPULP.VPOCHA(NLCE,1)
773 : uly = MPULP.VPOCHA(NLCE,2)
774 : p = MPPP.VPOCHA(NLCE,1)
775 : Tv = MPTVP.VPOCHA(NLCE,1)
776 : Tl = MPTLP.VPOCHA(NLCE,1)
777 : rv = MPRVP.VPOCHA(NLCE,1)
778 : rl = MPRLP.VPOCHA(NLCE,1)
779 : C
780 : C************* Value of the enthalpies and the derivatives of
781 : C the densities respect pressure and enthalpies,
782 : C for the moment only for perfect and stiffened gases.
783 : hv = Cpv*Tv
784 : hl = Cpl*Tl
785 : drvdp = rv/p
786 : drvdhv = -rv/hv
787 : cv = 1.d0/dsqrt(drvdhv/rv + drvdp)
788 : drldp = rl/(p + pil)
789 : drldhl = -rl/hl
790 : cl = 1.d0/dsqrt(drldhl/rl + drldp)
791 : dhvdp = 0.d0
792 : dhvdTv = Cpv
793 : dhldp = 0.d0
794 : dhldTl = Cpl
795 :
796 : denom = drldhl*(drvdhv*p + rv*(alpha*drvdp*p + rv -
797 : & alpha*rv)) + rl*(alpha*rl*rv*cv**(-2) -
798 : & (-1 + alpha)*drldp*(drvdhv*p + rv**2))
799 :
800 : A(1,1) = (alpha*(drldhl*p + rl**2)*rv*cv**(-2)*ulx -
801 : & (-1 + alpha)*rl*cl**(-2)*(drvdhv*p + rv**2)*uvx)
802 : & /denom
803 :
804 : A(1,2) = -((-1.d0 + alpha)*alpha*rl*cl**(-2)*rv**2)
805 : $ /denom
806 : A(1,3) = 0.d0
807 : A(1,4) = ((-1.d0 + alpha)*alpha*rl**2*rv*cv**(-2))
808 : $ /denom
809 : A(1,5) = 0.d0
810 : A(1,6) = ((-1.d0 + alpha)*alpha*rl*cl**(-2)*rv*cv**(-2)*(ulx
811 : $ -uvx))/denom
812 : A(1,7) = 0.d0
813 : A(1,8) = 0.d0
814 :
815 : A(2,1) = p/(alpha*rv)
816 : A(2,2) = uvx
817 : A(2,3) = 0.d0
818 : A(2,4) = 0.d0
819 : A(2,5) = 0.d0
820 : A(2,6) = 1.d0/rv
821 : A(2,7) = 0.d0
822 : A(2,8) = 0.d0
823 :
824 : A(3,1) = 0.d0
825 : A(3,2) = 0.d0
826 : A(3,3) = uvx
827 : A(3,4) = 0.d0
828 : A(3,5) = 0.d0
829 : A(3,6) = 0.d0
830 : A(3,7) = 0.d0
831 : A(3,8) = 0.d0
832 :
833 : A(4,1) = p/(-rl + alpha*rl)
834 : A(4,2) = 0.d0
835 : A(4,3) = 0.d0
836 : A(4,4) = ulx
837 : A(4,5) = 0.d0
838 : A(4,6) = 1.d0/rl
839 : A(4,7) = 0.d0
840 : A(4,8) = 0.d0
841 :
842 : A(5,1) = 0.d0
843 : A(5,2) = 0.d0
844 : A(5,3) = 0.d0
845 : A(5,4) = 0.d0
846 : A(5,5) = ulx
847 : A(5,6) = 0.d0
848 : A(5,7) = 0.d0
849 : A(5,8) = 0
850 :
851 : A(6,1) = -((drldhl*p + rl**2)*(drvdhv*p + rv**2)*(ulx - uvx)
852 : $ )/denom
853 : A(6,2) = (alpha*(drldhl*p + rl**2)*rv**2)/denom
854 : A(6,3) = 0.d0
855 : A(6,4) = -((-1.d0 + alpha)*rl**2*(drvdhv*p + rv**2))
856 : $ /denom
857 : A(6,5) = 0.d0
858 : A(6,6) = (-((-1 + alpha)*rl*cl**(-2)*(drvdhv*p + rv**2)*ulx)
859 : $ +alpha*(drldhl*p + rl**2)*rv*cv**(-2)*uvx)/denom
860 : A(6,7) = 0.d0
861 : A(6,8) = 0.d0
862 :
863 : A(7,1) = ((drldhl*p + rl**2)*(drvdp*p - rv)*(ulx - uvx))
864 : & /denom
865 : A(7,2) = ((drldhl*p + rl*(drldp*(p - alpha*p) + alpha*rl))
866 : $ *rv)/denom
867 : A(7,3) = 0.d0
868 : A(7,4) = ((-1.d0 + alpha)*rl**2*(drvdp*p - rv))/denom
869 : A(7,5) = 0.d0
870 : A(7,6) = ((-1.d0 + alpha)*rl*cl**(-2)*(drvdp*p - rv)*(ulx -
871 : $ uvx))/denom
872 : A(7,7) = uvx
873 : A(7,8) = 0.d0
874 :
875 : A(8,1) = ((drldp*p - rl)*(drvdhv*p + rv**2)*(ulx - uvx))
876 : & /denom
877 : A(8,2) = (alpha*(-(drldp*p) + rl)*rv**2)/denom
878 : A(8,3) = 0.d0
879 : A(8,4) = (rl*(drvdhv*p + rv*(alpha*drvdp*p + rv - alpha*rv))
880 : $ )/denom
881 : A(8,5) = 0.d0
882 : A(8,6) = (alpha*(drldp*p - rl)*rv*cv**(-2)*(ulx - uvx))
883 : & /denom
884 : A(8,7) = 0.d0
885 : A(8,8) = ulx
886 :
887 : B(1,1) = (alpha*(drldhl*p + rl**2)*rv*cv**(-2)*uly -
888 : & (-1.d0 + alpha)*rl*cl**(-2)*(drvdhv*p + rv**2)*uvy)
889 : & /denom
890 : B(1,2) = 0.d0
891 : B(1,3) = -((-1 + alpha)*alpha*rl*cl**(-2)*rv**2)/denom
892 : B(1,4) = 0.d0
893 : B(1,5) = ((-1.d0 + alpha)*alpha*rl**2*rv*cv**(-2))
894 : $ /denom
895 : B(1,6) = ((-1.d0 + alpha)*alpha*rl*cl**(-2)*rv*cv**(-2)*(uly
896 : $ -uvy))/denom
897 : B(1,7) = 0.d0
898 : B(1,8) = 0.d0
899 :
900 : B(2,1) = 0.d0
901 : B(2,2) = uvy
902 : B(2,3) = 0.d0
903 : B(2,4) = 0.d0
904 : B(2,5) = 0.d0
905 : B(2,6) = 0.d0
906 : B(2,7) = 0.d0
907 : B(2,8) = 0.d0
908 :
909 : B(3,1) = p/(alpha*rv)
910 : B(3,2) = 0.d0
911 : B(3,3) = uvy
912 : B(3,4) = 0.d0
913 : B(3,5) = 0.d0
914 : B(3,6) = 1.d0/rv
915 : B(3,7) = 0.d0
916 : B(3,8) = 0.d0
917 :
918 : B(4,1) = 0.d0
919 : B(4,2) = 0.d0
920 : B(4,3) = 0.d0
921 : B(4,4) = uly
922 : B(4,5) = 0.d0
923 : B(4,6) = 0.d0
924 : B(4,7) = 0.d0
925 : B(4,8) = 0.d0
926 :
927 : B(5,1) = p/(-rl + alpha*rl)
928 : B(5,2) = 0.d0
929 : B(5,3) = 0.d0
930 : B(5,4) = 0.d0
931 : B(5,5) = uly
932 : B(5,6) = 1.d0/rl
933 : B(5,7) = 0.d0
934 : B(5,8) = 0.d0
935 :
936 : B(6,1) = -((drldhl*p + rl**2)*(drvdhv*p + rv**2)*(uly - uvy)
937 : $ )/denom
938 : B(6,2) = 0.d0
939 : B(6,3) = (alpha*(drldhl*p + rl**2)*rv**2)/denom
940 : B(6,4) = 0.d0
941 : B(6,5) = -((-1.d0 + alpha)*rl**2*(drvdhv*p + rv**2))
942 : $ /denom
943 : B(6,6) = (-((-1.d0 + alpha)*rl*cl**(-2)*(drvdhv*p + rv**2)
944 : $ *uly) +alpha*(drldhl*p + rl**2)*rv*cv**(-2)*uvy)
945 : $ /denom
946 : B(6,7) = 0.d0
947 : B(6,8) = 0.d0
948 :
949 : B(7,1) = ((drldhl*p + rl**2)*(drvdp*p - rv)*(uly - uvy))
950 : & /denom
951 : B(7,2) = 0.d0
952 : B(7,3) = ((drldhl*p + rl*(drldp*(p - alpha*p) + alpha*rl))
953 : $ *rv)/denom
954 : B(7,4) = 0.d0
955 : B(7,5) = ((-1.d0 + alpha)*rl**2*(drvdp*p - rv))/denom
956 : B(7,6) = ((-1.d0 + alpha)*rl*cl**(-2)*(drvdp*p - rv)*(uly -
957 : $ uvy))/denom
958 : B(7,7) = uvy
959 : B(7,8) = 0.d0
960 :
961 : B(8,1) = ((drldp*p - rl)*(drvdhv*p + rv**2)*(uly - uvy))
962 : & /denom
963 : B(8,2) = 0.d0
964 : B(8,3) = (alpha*(-(drldp*p) + rl)*rv**2)/denom
965 : B(8,4) = 0.d0
966 : B(8,5) = (rl*(drvdhv*p + rv*(alpha*drvdp*p + rv - alpha*rv))
967 : $ )/denom
968 : B(8,6) = (alpha*(drldp*p - rl)*rv*cv**(-2)*(uly - uvy))
969 : & /denom
970 : B(8,7) = 0.d0
971 : B(8,8) = uly
972 : C
973 : C dVpdV
974 : C
975 : dVpdV(1,1) = 1.d0
976 : dVpdV(1,2) = 0.d0
977 : dVpdV(1,3) = 0.d0
978 : dVpdV(1,4) = 0.d0
979 : dVpdV(1,5) = 0.d0
980 : dVpdV(1,6) = 0.d0
981 : dVpdV(1,7) = 0.d0
982 : dVpdV(1,8) = 0.d0
983 :
984 : dVpdV(2,1) = 0.d0
985 : dVpdV(2,2) = 1.d0
986 : dVpdV(2,3) = 0.d0
987 : dVpdV(2,4) = 0.d0
988 : dVpdV(2,5) = 0.d0
989 : dVpdV(2,6) = 0.d0
990 : dVpdV(2,7) = 0.d0
991 : dVpdV(2,8) = 0.d0
992 :
993 : dVpdV(3,1) = 0.d0
994 : dVpdV(3,2) = 0.d0
995 : dVpdV(3,3) = 1.d0
996 : dVpdV(3,4) = 0.d0
997 : dVpdV(3,5) = 0.d0
998 : dVpdV(3,6) = 0.d0
999 : dVpdV(3,7) = 0.d0
1000 : dVpdV(3,8) = 0.d0
1001 :
1002 : dVpdV(4,1) = 0.d0
1003 : dVpdV(4,2) = 0.d0
1004 : dVpdV(4,3) = 0.d0
1005 : dVpdV(4,4) = 1.d0
1006 : dVpdV(4,5) = 0.d0
1007 : dVpdV(4,6) = 0.d0
1008 : dVpdV(4,7) = 0.d0
1009 : dVpdV(4,8) = 0.d0
1010 :
1011 : dVpdV(5,1) = 0.d0
1012 : dVpdV(5,2) = 0.d0
1013 : dVpdV(5,3) = 0.d0
1014 : dVpdV(5,4) = 0.d0
1015 : dVpdV(5,5) = 1.d0
1016 : dVpdV(5,6) = 0.d0
1017 : dVpdV(5,7) = 0.d0
1018 : dVpdV(5,8) = 0.d0
1019 :
1020 : dVpdV(6,1) = 0.d0
1021 : dVpdV(6,2) = 0.d0
1022 : dVpdV(6,3) = 0.d0
1023 : dVpdV(6,4) = 0.d0
1024 : dVpdV(6,5) = 0.d0
1025 : dVpdV(6,6) = 1.d0
1026 : dVpdV(6,7) = 0.d0
1027 : dVpdV(6,8) = 0.d0
1028 :
1029 : dVpdV(7,1) = 0.d0
1030 : dVpdV(7,2) = 0.d0
1031 : dVpdV(7,3) = 0.d0
1032 : dVpdV(7,4) = 0.d0
1033 : dVpdV(7,5) = 0.d0
1034 : dVpdV(7,6) = dhvdp
1035 : dVpdV(7,7) = dhvdTv
1036 : dVpdV(7,8) = 0.d0
1037 :
1038 : dVpdV(8,1) = 0.d0
1039 : dVpdV(8,2) = 0.d0
1040 : dVpdV(8,3) = 0.d0
1041 : dVpdV(8,4) = 0.d0
1042 : dVpdV(8,5) = 0.d0
1043 : dVpdV(8,6) = dhldp
1044 : dVpdV(8,7) = 0.d0
1045 : dVpdV(8,8) = dhldTl
1046 : C
1047 : C dVdVp
1048 : C
1049 : dVdVp(1,1) = 1.d0
1050 : dVdVp(1,2) = 0.d0
1051 : dVdVp(1,3) = 0.d0
1052 : dVdVp(1,4) = 0.d0
1053 : dVdVp(1,5) = 0.d0
1054 : dVdVp(1,6) = 0.d0
1055 : dVdVp(1,7) = 0.d0
1056 : dVdVp(1,8) = 0.d0
1057 :
1058 : dVdVp(2,1) = 0.d0
1059 : dVdVp(2,2) = 1.d0
1060 : dVdVp(2,3) = 0.d0
1061 : dVdVp(2,4) = 0.d0
1062 : dVdVp(2,5) = 0.d0
1063 : dVdVp(2,6) = 0.d0
1064 : dVdVp(2,7) = 0.d0
1065 : dVdVp(2,8) = 0.d0
1066 :
1067 : dVdVp(3,1) = 0.d0
1068 : dVdVp(3,2) = 0.d0
1069 : dVdVp(3,3) = 1.d0
1070 : dVdVp(3,4) = 0.d0
1071 : dVdVp(3,5) = 0.d0
1072 : dVdVp(3,6) = 0.d0
1073 : dVdVp(3,7) = 0.d0
1074 : dVdVp(3,8) = 0.d0
1075 :
1076 : dVdVp(4,1) = 0.d0
1077 : dVdVp(4,2) = 0.d0
1078 : dVdVp(4,3) = 0.d0
1079 : dVdVp(4,4) = 1.d0
1080 : dVdVp(4,5) = 0.d0
1081 : dVdVp(4,6) = 0.d0
1082 : dVdVp(4,7) = 0.d0
1083 : dVdVp(4,8) = 0.d0
1084 :
1085 : dVdVp(5,1) = 0.d0
1086 : dVdVp(5,2) = 0.d0
1087 : dVdVp(5,3) = 0.d0
1088 : dVdVp(5,4) = 0.d0
1089 : dVdVp(5,5) = 1.d0
1090 : dVdVp(5,6) = 0.d0
1091 : dVdVp(5,7) = 0.d0
1092 : dVdVp(5,8) = 0.d0
1093 :
1094 : dVdVp(6,1) = 0.d0
1095 : dVdVp(6,2) = 0.d0
1096 : dVdVp(6,3) = 0.d0
1097 : dVdVp(6,4) = 0.d0
1098 : dVdVp(6,5) = 0.d0
1099 : dVdVp(6,6) = 1.d0
1100 : dVdVp(6,7) = 0.d0
1101 : dVdVp(6,8) = 0.d0
1102 :
1103 : dVdVp(7,1) = 0.d0
1104 : dVdVp(7,2) = 0.d0
1105 : dVdVp(7,3) = 0.d0
1106 : dVdVp(7,4) = 0.d0
1107 : dVdVp(7,5) = 0.d0
1108 : dVdVp(7,6) = -dhvdp/dhvdTv
1109 : dVdVp(7,7) = 1.d0/dhvdTv
1110 : dVdVp(7,8) = 0.d0
1111 :
1112 : dVdVp(8,1) = 0.d0
1113 : dVdVp(8,2) = 0.d0
1114 : dVdVp(8,3) = 0.d0
1115 : dVdVp(8,4) = 0.d0
1116 : dVdVp(8,5) = 0.d0
1117 : dVdVp(8,6) = dhldp*dhvdp/(dhldTl*dhvdTv)
1118 : dVdVp(8,7) = -dhldp/(dhldTl*dhvdTv)
1119 : dVdVp(8,8) = 1.d0/dhldTl
1120 : C
1121 : C dVdVp*A*dVpdV
1122 : C
1123 : do l = 1, 8
1124 : do m = 1, 8
1125 : P1(l,m) = 0.d0
1126 : end do
1127 : end do
1128 : do l = 1, 8
1129 : do m = 1, 8
1130 : do n = 1, 8
1131 : P1(l,m) = P1(l,m) + dVdVp(l,n)*A(n,m)
1132 : end do
1133 : end do
1134 : end do
1135 : do l = 1, 8
1136 : do m = 1, 8
1137 : P2(l,m) = 0.d0
1138 : end do
1139 : end do
1140 : do l = 1, 8
1141 : do m = 1, 8
1142 : do n = 1, 8
1143 : P2(l,m) = P2(l,m) + P1(l,n)*dVpdV(n,m)
1144 : end do
1145 : end do
1146 : end do
1147 : do l = 1, 8
1148 : do m = 1, 8
1149 : A(l,m) = P2(l,m)
1150 : end do
1151 : end do
1152 : C
1153 : C dVdVp*B*dVpdV
1154 : C
1155 : do l = 1, 8
1156 : do m = 1, 8
1157 : P1(l,m) = 0.d0
1158 : end do
1159 : end do
1160 : do l = 1, 8
1161 : do m = 1, 8
1162 : do n = 1, 8
1163 : P1(l,m) = P1(l,m) + dVdVp(l,n)*B(n,m)
1164 : end do
1165 : end do
1166 : end do
1167 : do l = 1, 8
1168 : do m = 1, 8
1169 : P2(l,m) = 0.d0
1170 : end do
1171 : end do
1172 : do l = 1, 8
1173 : do m = 1, 8
1174 : do n = 1, 8
1175 : P2(l,m) = P2(l,m) + P1(l,n)*dVpdV(n,m)
1176 : end do
1177 : end do
1178 : end do
1179 : do l = 1, 8
1180 : do m = 1, 8
1181 : B(l,m) = P2(l,m)
1182 : end do
1183 : end do
1184 : C
1185 : C MEJORAR LO DE LA MULTIPLICACION DE LAS MATRICES DE CAMBIO
1186 : C
1187 : dVdx(1) = MPGRAL.VPOCHA(NLCE,1)*LIMITA(1)
1188 : dVdx(2) = MPGRUV.VPOCHA(NLCE,1)*LIMITA(2)
1189 : dVdx(3) = MPGRUV.VPOCHA(NLCE,3)*LIMITA(3)
1190 : dVdx(4) = MPGRUL.VPOCHA(NLCE,1)*LIMITA(4)
1191 : dVdx(5) = MPGRUL.VPOCHA(NLCE,3)*LIMITA(5)
1192 : dVdx(6) = MPGRP.VPOCHA(NLCE,1)*LIMITA(6)
1193 : dVdx(7) = MPGRTV.VPOCHA(NLCE,1)*LIMITA(7)
1194 : dVdx(8) = MPGRTL.VPOCHA(NLCE,1)*LIMITA(8)
1195 :
1196 : dVdy(1) = MPGRAL.VPOCHA(NLCE,2)*LIMITA(1)
1197 : dVdy(2) = MPGRUV.VPOCHA(NLCE,2)*LIMITA(2)
1198 : dVdy(3) = MPGRUV.VPOCHA(NLCE,4)*LIMITA(3)
1199 : dVdy(4) = MPGRUL.VPOCHA(NLCE,2)*LIMITA(4)
1200 : dVdy(5) = MPGRUL.VPOCHA(NLCE,4)*LIMITA(5)
1201 : dVdy(6) = MPGRP.VPOCHA(NLCE,2)*LIMITA(6)
1202 : dVdy(7) = MPGRTV.VPOCHA(NLCE,2)*LIMITA(7)
1203 : dVdy(8) = MPGRTL.VPOCHA(NLCE,2)*LIMITA(8)
1204 :
1205 :
1206 : C No te olvides de afectar de la correspondiente cp las
1207 : C dVdx y dVdy correspondientes a las entalpias, no vayas
1208 : C a meter la pata con eso
1209 : do m = 1, 8
1210 : dV(m) = 0.d0
1211 : end do
1212 : do m = 1, 8
1213 : do n = 1, 8
1214 : dV(m) = dV(m) + A(m,n)*dVdx(n) + B(m,n)*dVdy(n)
1215 : end do
1216 : end do
1217 : C
1218 : C************* Prediction avec/sans gradients limités
1219 : C
1220 : MPALPP.VPOCHA(NLCE,1)= alpha - DELTAT*dV(1)
1221 : MPUVP.VPOCHA(NLCE,1) = uvx - DELTAT*dV(2)
1222 : MPUVP.VPOCHA(NLCE,2) = uvy - DELTAT*dV(3)
1223 : MPULP.VPOCHA(NLCE,1) = ulx - DELTAT*dV(4)
1224 : MPULP.VPOCHA(NLCE,2) = uly - DELTAT*dV(5)
1225 : MPPP.VPOCHA(NLCE,1) = p - DELTAT*dV(6)
1226 : MPTVP.VPOCHA(NLCE,1) = Tv - DELTAT*dV(7)
1227 : MPTLP.VPOCHA(NLCE,1) = Tl - DELTAT*dV(8)
1228 : C
1229 : C************* Updating densities(perfect and stiffened gas eos
1230 : C
1231 : MPRVP.VPOCHA(NLCE,1) = gamv*MPPP.VPOCHA(NLCE,1)/
1232 : & (Cpv*MPTVP.VPOCHA(NLCE,1)*(gamv - 1.d0))
1233 : MPRLP.VPOCHA(NLCE,1) = gaml*(MPPP.VPOCHA(NLCE,1) + pil)/
1234 : & (Cpl* MPTLP.VPOCHA(NLCE,1)*(gaml - 1.d0))
1235 :
1236 : ENDDO
1237 : C
1238 : ENDIF
1239 : C
1240 : C
1241 : C***********************************************************************
1242 : C********* CORRECTION **************************************************
1243 : C***********************************************************************
1244 : C
1245 : C**** Boucle sur le faces
1246 : C
1247 : IDIMP1 = IDIM + 1
1248 : DO NLCF = 1, NFAC
1249 : C
1250 : C******* NLCF = numero local du centre de face
1251 : C NGCF = numero global du centre de face
1252 : C NGCEG = numero global du centre ELT "gauche"
1253 : C NLCEG = numero local du centre ELT "gauche"
1254 : C NGCED = numero global du centre ELT "droite"
1255 : C NLCED = numero local du centre ELT "droite"
1256 : C
1257 : NGCEG = IPT1.NUM(1,NLCF)
1258 : NGCF = IPT1.NUM(2,NLCF)
1259 : NGCED = IPT1.NUM(3,NLCF)
1260 : NLCEG = MLENT1.LECT(NGCEG)
1261 : NLCED = MLENT1.LECT(NGCED)
1262 : C
1263 : C******* TEST: IPT2.NUM(1,NLCF) = IPT1.NUM(2,NLCF)
1264 : C
1265 : NGCF1 = IPT2.NUM(1,NLCF)
1266 : IF( NGCF1 .NE. NGCF) THEN
1267 : LOGAN = .TRUE.
1268 : MESERR(1:40) = 'PRET, subroutine pre52f.eso '
1269 : GOTO 9999
1270 : ENDIF
1271 : C
1272 : C******* Cosinus directeurs des NORMALES aux faces
1273 : C
1274 : C On impose que les normales sont direct "Gauche" -> "Centre"
1275 : C
1276 : INDCEL = (NGCEG-1)*IDIMP1
1277 : XG = XCOOR(INDCEL+1)
1278 : YG = XCOOR(INDCEL+2)
1279 : INDCEL = (NGCF-1)*IDIMP1
1280 : XC = XCOOR(INDCEL + 1)
1281 : YC = XCOOR(INDCEL + 2)
1282 : INDCEL = (NGCED-1)*IDIMP1
1283 : XD = XCOOR(INDCEL+1)
1284 : YD = XCOOR(INDCEL+2)
1285 : DXG = XC - XG
1286 : DYG = YC - YG
1287 : DXD = XC - XD
1288 : DYD = YC - YD
1289 : C
1290 : C******* On calcule le sign du pruduit scalare
1291 : C (Normales de Castem) * (vecteur "gauche" -> "centre")
1292 : C
1293 : CNX = MPNORM.VPOCHA(NLCF,1)
1294 : CNY = MPNORM.VPOCHA(NLCF,2)
1295 : ORIENT = CNX * DXG + CNY * DYG
1296 : ORIENT = SIGN(1.0D0,ORIENT)
1297 : IF(ORIENT .NE. 1.0D0)THEN
1298 : LOGAN = .TRUE.
1299 : MESERR(1:30)=
1300 : & 'PRET , subroutine pre52f.eso. '
1301 : GOTO 9999
1302 : ENDIF
1303 : CNX = CNX * ORIENT
1304 : CNY = CNY * ORIENT
1305 : C
1306 : C********** Cosinus directeurs de tangent 2D
1307 : C
1308 : CTX = -1.0D0 * CNY
1309 : CTY = CNX
1310 : C
1311 : C
1312 : C******* Les autres MELVALs
1313 : C
1314 : C
1315 : C******* N.B.: On suppose qu'on a déjà controlle RO, P > 0
1316 : C GAMMA \in (1,3)
1317 : C Si non il faut le faire, en utilisant LOGBOR,
1318 : C LOGNEG, VALER, VAL1, VAL2
1319 : C
1320 : C
1321 : C
1322 : C******* NGCEG = NGCED -> Mur
1323 : C
1324 : IF( NGCEG .EQ. NGCED)THEN
1325 : C
1326 : C********** Sur le mur on fait de reconstruction sur l'etat gauche
1327 : C
1328 : C
1329 : C********** Etat gauche
1330 : C
1331 :
1332 : VALCEL = MPALPP.VPOCHA(NLCEG,1)
1333 : ALCEL = MELLAL.VPOCHA(NLCEG,1)
1334 : DCEL = MPGRAL.VPOCHA(NLCEG,1)*DXG +
1335 : & MPGRAL.VPOCHA(NLCEG,2)*DYG
1336 : alphag = VALCEL + ALCEL*DCEL
1337 :
1338 : VALCEL = MPUVP.VPOCHA(NLCEG,1)
1339 : ALCEL = MELLUV.VPOCHA(NLCEG,1)
1340 : DCEL = MPGRUV.VPOCHA(NLCEG,1)*DXG +
1341 : & MPGRUV.VPOCHA(NLCEG,2)*DYG
1342 : uvxg = VALCEL + ALCEL*DCEL
1343 :
1344 : VALCEL = MPUVP.VPOCHA(NLCEG,2)
1345 : ALCEL = MELLUV.VPOCHA(NLCEG,2)
1346 : DCEL = MPGRUV.VPOCHA(NLCEG,3)*DXG +
1347 : & MPGRUV.VPOCHA(NLCEG,4)*DYG
1348 : uvyg = VALCEL + ALCEL*DCEL
1349 :
1350 : uvng = uvxg*CNX + uvyg*CNY
1351 : uvtg = uvxg*CTX + uvyg*CTY
1352 :
1353 : VALCEL = MPULP.VPOCHA(NLCEG,1)
1354 : ALCEL = MELLUL.VPOCHA(NLCEG,1)
1355 : DCEL = MPGRUL.VPOCHA(NLCEG,1)*DXG +
1356 : & MPGRUL.VPOCHA(NLCEG,2)*DYG
1357 : ulxg = VALCEL + ALCEL*DCEL
1358 :
1359 : VALCEL = MPULP.VPOCHA(NLCEG,2)
1360 : ALCEL = MELLUL.VPOCHA(NLCEG,2)
1361 : DCEL = MPGRUL.VPOCHA(NLCEG,3)*DXG +
1362 : & MPGRUL.VPOCHA(NLCEG,4)*DYG
1363 : ulyg = VALCEL + ALCEL*DCEL
1364 :
1365 : ulng = ulxg*CNX + ulyg*CNY
1366 : ultg = ulxg*CTX + ulyg*CTY
1367 :
1368 : VALCEL = MPPP.VPOCHA(NLCEG, 1)
1369 : ALCEL = MELLP.VPOCHA(NLCEG, 1)
1370 : DCEL = MPGRP.VPOCHA(NLCEG, 1)*DXG +
1371 : & MPGRP.VPOCHA(NLCEG, 2)*DYG
1372 : pg = VALCEL + ALCEL*DCEL
1373 :
1374 : VALCEL = MPTVP.VPOCHA(NLCEG, 1)
1375 : ALCEL = MELLTV.VPOCHA(NLCEG, 1)
1376 : DCEL = MPGRTV.VPOCHA(NLCEG, 1)*DXG +
1377 : & MPGRTV.VPOCHA(NLCEG, 2)*DYG
1378 : Tvg = VALCEL + ALCEL*DCEL
1379 :
1380 : VALCEL = MPTLP.VPOCHA(NLCEG, 1)
1381 : ALCEL = MELLTV.VPOCHA(NLCEG, 1)
1382 : DCEL = MPGRTL.VPOCHA(NLCEG, 1)*DXG +
1383 : & MPGRTL.VPOCHA(NLCEG, 2)*DYG
1384 : Tlg = VALCEL + ALCEL*DCEL
1385 :
1386 : C********** We calculate densities again, we use for
1387 : C the moment, perfect and stiffened gas eos
1388 : rvg = gamv*pg/(cpv*Tvg*(gamv - 1))
1389 : rlg = gaml*(pg + pil)/(cpl*Tlg*(gaml - 1))
1390 : C
1391 : C********** Si l'on fait pas de prediction, ce n'est pas necessaire de
1392 : C controller la positivite' de la pression et de la densité; elle
1393 : C est déjà garantie par la proprieté LED de limiteur; mais, vu
1394 : C que le limiteur n'est pas calculé ici, mais dans un autre
1395 : C operateur, on le fait
1396 : C
1397 : LOGI1 = (pg .LT. 0.0D0) .OR. (rvg .LT. 0.0D0)
1398 : & .OR. (rlg .LT. 0.0D0)
1399 : C
1400 : IF(LOGI1)THEN
1401 : C
1402 : C************* Premier ordre en espace local
1403 : C
1404 : alphag = MPALPC.VPOCHA(NLCEG,1)
1405 : uvng = MPUVC.VPOCHA(NLCEG,1)*CNX +
1406 : & MPUVC.VPOCHA(NLCEG,2)*CNY
1407 : uvtg = MPULC.VPOCHA(NLCEG,1)*CTX +
1408 : & MPUVC.VPOCHA(NLCEG,2)*CTY
1409 : ulng = MPULC.VPOCHA(NLCEG,1)*CNX +
1410 : & MPULC.VPOCHA(NLCEG,2)*CNY
1411 : ultg = MPULC.VPOCHA(NLCEG,1)*CTX +
1412 : & MPULC.VPOCHA(NLCEG,2)*CTY
1413 : pg = MPPC.VPOCHA(NLCEG,1)
1414 : Tvg = MPTVC.VPOCHA(NLCEG,1)
1415 : Tlg = MPTLC.VPOCHA(NLCEG,1)
1416 : rvg = MPRVC.VPOCHA(NLCEG,1)
1417 : rlg = MPRLC.VPOCHA(NLCEG,1)
1418 :
1419 : ENDIF
1420 : C
1421 : C********** Son etat droite
1422 : C
1423 : alphad = alphag
1424 : uvnd = -1.d0*uvng
1425 : uvtd = uvtg
1426 : ulnd = -1.d0*ulng
1427 : ultd = ultg
1428 : pd = pg
1429 : Tvd = Tvg
1430 : Tld = Tlg
1431 : rvd = rvg
1432 : rld = rlg
1433 : C
1434 : C************* Fin cas mur
1435 : C
1436 : ELSE
1437 : C
1438 : C************* Etat gauche
1439 : C
1440 : VALCEL = MPALPP.VPOCHA(NLCEG,1)
1441 : ALCEL = MELLAL.VPOCHA(NLCEG,1)
1442 : DCEL = MPGRAL.VPOCHA(NLCEG,1)*DXG +
1443 : & MPGRAL.VPOCHA(NLCEG,2)*DYG
1444 : alphag = VALCEL + ALCEL*DCEL
1445 :
1446 : VALCEL = MPUVP.VPOCHA(NLCEG,1)
1447 : ALCEL = MELLUV.VPOCHA(NLCEG,1)
1448 : DCEL = MPGRUV.VPOCHA(NLCEG,1)*DXG +
1449 : & MPGRUV.VPOCHA(NLCEG,2)*DYG
1450 : uvxg = VALCEL + ALCEL*DCEL
1451 :
1452 : VALCEL = MPUVP.VPOCHA(NLCEG,2)
1453 : ALCEL = MELLUV.VPOCHA(NLCEG,2)
1454 : DCEL = MPGRUV.VPOCHA(NLCEG,3)*DXG +
1455 : & MPGRUV.VPOCHA(NLCEG,4)*DYG
1456 : uvyg = VALCEL + ALCEL*DCEL
1457 :
1458 : uvng = uvxg*CNX + uvyg*CNY
1459 : uvtg = uvxg*CTX + uvyg*CTY
1460 :
1461 : VALCEL = MPULP.VPOCHA(NLCEG,1)
1462 : ALCEL = MELLUL.VPOCHA(NLCEG,1)
1463 : DCEL = MPGRUL.VPOCHA(NLCEG,1)*DXG +
1464 : & MPGRUL.VPOCHA(NLCEG,2)*DYG
1465 : ulxg = VALCEL + ALCEL*DCEL
1466 :
1467 : VALCEL = MPULP.VPOCHA(NLCEG,2)
1468 : ALCEL = MELLUL.VPOCHA(NLCEG,2)
1469 : DCEL = MPGRUL.VPOCHA(NLCEG,3)*DXG +
1470 : & MPGRUL.VPOCHA(NLCEG,4)*DYG
1471 : ulyg = VALCEL + ALCEL*DCEL
1472 :
1473 : ulng = ulxg*CNX + ulyg*CNY
1474 : ultg = ulxg*CTX + ulyg*CTY
1475 :
1476 : VALCEL = MPPP.VPOCHA(NLCEG, 1)
1477 : ALCEL = MELLP.VPOCHA(NLCEG, 1)
1478 : DCEL = MPGRP.VPOCHA(NLCEG, 1)*DXG +
1479 : & MPGRP.VPOCHA(NLCEG, 2)*DYG
1480 : pg = VALCEL + ALCEL*DCEL
1481 :
1482 : VALCEL = MPTVP.VPOCHA(NLCEG, 1)
1483 : ALCEL = MELLTV.VPOCHA(NLCEG, 1)
1484 : DCEL = MPGRTV.VPOCHA(NLCEG, 1)*DXG +
1485 : & MPGRTV.VPOCHA(NLCEG, 2)*DYG
1486 : Tvg = VALCEL + ALCEL*DCEL
1487 :
1488 : VALCEL = MPTLP.VPOCHA(NLCEG, 1)
1489 : ALCEL = MELLTL.VPOCHA(NLCEG, 1)
1490 : DCEL = MPGRTL.VPOCHA(NLCEG, 1)*DXG +
1491 : & MPGRTL.VPOCHA(NLCEG, 2)*DYG
1492 : Tlg = VALCEL + ALCEL*DCEL
1493 :
1494 : C********** We calculate densities again, we use for
1495 : C the moment, perfect and stiffened gas eos
1496 : rvg = gamv*pg/(cpv*Tvg*(gamv - 1))
1497 : rlg = gaml*(pg + pil)/(cpl*Tlg*(gaml - 1))
1498 : C
1499 : C********** Positivite
1500 : C
1501 : LOGI1 = (pg .LT. 0.D0) .OR. (rvg .LT. 0.D0)
1502 : & .OR. (rlg .LT. 0.D0)
1503 : C
1504 : IF(LOGI1)THEN
1505 : C
1506 : C************* Premier ordre en espace local
1507 : C
1508 : alphag = MPALPC.VPOCHA(NLCEG,1)
1509 : uvng = MPUVC.VPOCHA(NLCEG,1)*CNX +
1510 : & MPUVC.VPOCHA(NLCEG,2)*CNY
1511 : uvtg = MPULC.VPOCHA(NLCEG,1)*CTX +
1512 : & MPUVC.VPOCHA(NLCEG,2)*CTY
1513 : ulng = MPULC.VPOCHA(NLCEG,1)*CNX +
1514 : & MPULC.VPOCHA(NLCEG,2)*CNY
1515 : ultg = MPULC.VPOCHA(NLCEG,1)*CTX +
1516 : & MPULC.VPOCHA(NLCEG,2)*CTY
1517 : pg = MPPC.VPOCHA(NLCEG,1)
1518 : Tvg = MPTVC.VPOCHA(NLCEG,1)
1519 : Tlg = MPTLC.VPOCHA(NLCEG,1)
1520 : rvg = MPRVC.VPOCHA(NLCEG,1)
1521 : rlg = MPRLC.VPOCHA(NLCEG,1)
1522 :
1523 : ENDIF
1524 : C
1525 : C********** Etat droite
1526 : C
1527 : VALCEL = MPALPP.VPOCHA(NLCED, 1)
1528 : ALCEL = MELLAL.VPOCHA(NLCED, 1)
1529 : DCEL = MPGRAL.VPOCHA(NLCED, 1)*DXD +
1530 : & MPGRAL.VPOCHA(NLCED, 2)*DYD
1531 : alphad = VALCEL + ALCEL * DCEL
1532 :
1533 : VALCEL = MPUVP.VPOCHA(NLCED, 1)
1534 : ALCEL = MELLUV.VPOCHA(NLCED, 1)
1535 : DCEL = MPGRUV.VPOCHA(NLCED, 1)*DXD +
1536 : & MPGRUV.VPOCHA(NLCED, 2)*DYD
1537 : uvxd = VALCEL + ALCEL * DCEL
1538 :
1539 : VALCEL = MPUVP.VPOCHA(NLCED, 2)
1540 : ALCEL = MELLUV.VPOCHA(NLCED, 2)
1541 : DCEL = MPGRUV.VPOCHA(NLCED, 3)*DXD +
1542 : & MPGRUV.VPOCHA(NLCED, 4)*DYD
1543 : uvyd = VALCEL + ALCEL * DCEL
1544 :
1545 : uvnd = uvxd*CNX + uvyd*CNY
1546 : uvtd = uvxd*CTX + uvyd*CTY
1547 :
1548 : VALCEL = MPULP.VPOCHA(NLCED, 1)
1549 : ALCEL = MELLUL.VPOCHA(NLCED, 1)
1550 : DCEL = MPGRUL.VPOCHA(NLCED, 1)*DXD +
1551 : & MPGRUL.VPOCHA(NLCED, 2)*DYD
1552 : ulxd = VALCEL + ALCEL * DCEL
1553 :
1554 : VALCEL = MPULP.VPOCHA(NLCED,2)
1555 : ALCEL = MELLUL.VPOCHA(NLCED, 2)
1556 : DCEL = MPGRUL.VPOCHA(NLCED, 3)*DXD +
1557 : & MPGRUL.VPOCHA(NLCED, 4)*DYD
1558 : ulyd = VALCEL + ALCEL * DCEL
1559 :
1560 : ulnd = ulxd*CNX + ulyd*CNY
1561 : ultd = ulxd*CTX + ulyd*CTY
1562 :
1563 : VALCEL = MPPP.VPOCHA(NLCED, 1)
1564 : ALCEL = MELLP.VPOCHA(NLCED, 1)
1565 : DCEL = MPGRP.VPOCHA(NLCED, 1)*DXD +
1566 : & MPGRP.VPOCHA(NLCED, 2)*DYD
1567 : pd = VALCEL + ALCEL * DCEL
1568 :
1569 : VALCEL = MPTVP.VPOCHA(NLCED, 1)
1570 : ALCEL = MELLTV.VPOCHA(NLCED, 1)
1571 : DCEL = MPGRTV.VPOCHA(NLCED, 1)*DXD +
1572 : & MPGRTV.VPOCHA(NLCED, 2)*DYD
1573 : Tvd = VALCEL + ALCEL * DCEL
1574 :
1575 : VALCEL = MPTLP.VPOCHA(NLCED, 1)
1576 : ALCEL = MELLTL.VPOCHA(NLCED, 1)
1577 : DCEL = MPGRTL.VPOCHA(NLCED,1)*DXD +
1578 : & MPGRTL.VPOCHA(NLCED,2)*DYD
1579 : Tld = VALCEL + ALCEL * DCEL
1580 : C
1581 : C********** We calculate densities again, we use for
1582 : C the moment, perfect and stiffened gas eos
1583 : rvd = gamv*pd/(cpv*Tvd*(gamv - 1))
1584 : rld = gaml*(pd + pil)/(cpl*Tld*(gaml - 1))
1585 : C
1586 : C********** Positivite
1587 : C
1588 : LOGI1 = (pd .LT. 0.0D0) .OR. (rvd .LT. 0.0D0)
1589 : & .OR. (rld .LT. 0.0D0)
1590 : C
1591 : IF(LOGI1)THEN
1592 : C
1593 : C************* Premier ordre en espace local
1594 : C
1595 : alphad = MPALPC.VPOCHA(NLCED,1)
1596 : uvnd = MPUVC.VPOCHA(NLCED,1)*CNX +
1597 : & MPUVC.VPOCHA(NLCED,2)*CNY
1598 : uvtd = MPUVC.VPOCHA(NLCED,1)*CTX +
1599 : & MPUVC.VPOCHA(NLCED,2)*CTY
1600 : pd = MPPC.VPOCHA(NLCED,1)
1601 : Tvd = MPTVC.VPOCHA(NLCED,1)
1602 : Tld = MPTLC.VPOCHA(NLCED,1)
1603 : rvd = MPRVC.VPOCHA(NLCED,1)
1604 : rld = MPRLC.VPOCHA(NLCED,1)
1605 :
1606 : ENDIF
1607 : ENDIF
1608 : C
1609 : C******** Les MELVALs
1610 : C
1611 : MLALP.VELCHE(1,NLCF) = alphag
1612 : MLALP.VELCHE(3,NLCF) = alphad
1613 : MLP.VELCHE(1,NLCF) = pg
1614 : MLP.VELCHE(3,NLCF) = pd
1615 : MLTV.VELCHE(1,NLCF) = Tvg
1616 : MLTV.VELCHE(3,NLCF) = Tvd
1617 : MLTL.VELCHE(1,NLCF) = Tlg
1618 : MLTL.VELCHE(3,NLCF) = Tld
1619 : MLRV.VELCHE(1,NLCF) = rvg
1620 : MLRV.VELCHE(3,NLCF) = rvd
1621 : MLRL.VELCHE(1,NLCF) = rlg
1622 : MLRL.VELCHE(3,NLCF) = rld
1623 :
1624 : MLUVN.VELCHE(1,NLCF) = uvng
1625 : MLUVN.VELCHE(3,NLCF) = uvnd
1626 : MLUVT.VELCHE(1,NLCF) = uvtg
1627 : MLUVT.VELCHE(3,NLCF) = uvtd
1628 : MLVNX.VELCHE(1,NLCF) = CNX
1629 : MLVNY.VELCHE(1,NLCF) = CNY
1630 : MLVTX.VELCHE(1,NLCF) = CTX
1631 : MLVTY.VELCHE(1,NLCF) = CTY
1632 :
1633 : MLULN.VELCHE(1,NLCF) = ulng
1634 : MLULN.VELCHE(3,NLCF) = ulnd
1635 : MLULT.VELCHE(1,NLCF) = ultg
1636 : MLULT.VELCHE(3,NLCF) = ultd
1637 : MLLNX.VELCHE(1,NLCF) = CNX
1638 : MLLNY.VELCHE(1,NLCF) = CNY
1639 : MLLTX.VELCHE(1,NLCF) = CTX
1640 : MLLTY.VELCHE(1,NLCF) = CTY
1641 : C
1642 : ENDDO
1643 : C
1644 : C**** Desactivation des SEGMENTs
1645 : C
1646 : SEGDES IPT1
1647 : SEGDES IPT2
1648 : C
1649 : C**** Le MPOVALs 'Prediction' sont detruits (si existentes)
1650 : C
1651 : IF(LOGTEM)THEN
1652 : SEGSUP MPALPP
1653 : SEGSUP MPUVP
1654 : SEGSUP MPULP
1655 : SEGSUP MPPP
1656 : SEGSUP MPTVP
1657 : SEGSUP MPTLP
1658 : SEGSUP MPRVP
1659 : SEGSUP MPRLP
1660 : ENDIF
1661 : C
1662 : SEGDES MPALPC
1663 : SEGDES MPGRAL
1664 : SEGDES MPUVC
1665 : SEGDES MPGRUV
1666 : SEGDES MPULC
1667 : SEGDES MPGRUL
1668 : SEGDES MPPC
1669 : SEGDES MPGRP
1670 : SEGDES MPTVC
1671 : SEGDES MPGRTV
1672 : SEGDES MPTLC
1673 : SEGDES MPGRTL
1674 : SEGDES MPRVC
1675 : SEGDES MPRLC
1676 : SEGDES MPNORM
1677 : SEGDES MELLAL
1678 : SEGDES MELLUV
1679 : SEGDES MELLUL
1680 : SEGDES MELLP
1681 : SEGDES MELLTV
1682 : SEGDES MELLTL
1683 : C
1684 : SEGDES MLALP
1685 : SEGDES MLP
1686 : SEGDES MLTV
1687 : SEGDES MLTL
1688 : SEGDES MLRV
1689 : SEGDES MLRL
1690 : SEGDES MLUVN
1691 : SEGDES MLUVT
1692 : SEGDES MLULN
1693 : SEGDES MLULT
1694 : SEGDES MLVNX
1695 : SEGDES MLVNY
1696 : SEGDES MLVTX
1697 : SEGDES MLVTY
1698 : SEGDES MLLNX
1699 : SEGDES MLLNY
1700 : SEGDES MLLTX
1701 : SEGDES MLLTY
1702 : C
1703 : C**** Destruction du MELNTI correspondance local/global
1704 : C
1705 : SEGSUP MLENT1
1706 : C
1707 : 9999 CONTINUE
1708 : C
1709 : RETURN
1710 : END
1711 :
1712 :
1713 :
1714 :
1715 :
1716 :
© Cast3M 2003 - All rights reserved.
Disclaimer