1 : C PRE62F SOURCE PV 09/03/12 21:31:14 6325 2 : SUBROUTINE PRE62F(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 : PRE62F 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 32 : C APPELES (Outils) : KRIPAD, LICHT 33 : C 34 : C APPELES (Calcul) : AUCUN 35 : C 36 : C 37 : C************************************************************************ 38 : C 39 : C ENTREES 40 : C 41 : C LOGTEM : LOGICAL; si .TRUE. 2em ordre en temps 42 : C sinon 1er ordre en temps 43 : C 44 : C 1) Pointeurs de MELEMEs et de CHPOINTs de la table DOMAINE 45 : C 46 : C ICEN : MELEME de 'POI1' SPG des CENTRES 47 : C 48 : C IFACE : MELEME de 'POI1' SPG des FACES 49 : C 50 : C IFACEL : MELEME de 'SEG3' avec 51 : C CENTRE d'Elt "gauche" 52 : C CENTRE de Face 53 : C CENTRE d'Elt "droite" 54 : C 55 : C N.B. = IFACE.NUM(i,1) = IFACEL.NUM(i,2) 56 : C 57 : C INORM : CHPOINT des cosinus directeurs de normales aux faces 58 : C 59 : C 2) Pointeurs des CHPOINTs 60 : C 61 : C IALPH : CHPOINT "CENTRE" containing void fraction 62 : C 63 : C IGRALP : CHPOINT "CENTRE" contaninig the void fraction gradient 64 : C (2 composantes) 65 : C 66 : C IALALP : CHPOINT "CENTRE" contaning the limiter of the 67 : C void fraction gradient 68 : C 69 : C IUVC : CHPOINT "CENTRE" containing the vapour velocity UvX, UVY ; 70 : C 71 : C IGRUVC : CHPOINT "CENTRE" containig the vapour velocity gradient 72 : C (4 composantes) 73 : C 74 : C IALUVC : CHPOINT "CENTRE" contaning the limiter of the vapour 75 : C velocity (2 composantes) 76 : C 77 : C IULC : CHPOINT "CENTRE" containing the liquid velocity UvX, UVY ; 78 : C 79 : C IGRULC : CHPOINT "CENTRE" containig the liquid velocity gradient 80 : C (4 composantes) 81 : C 82 : C IALULC : CHPOINT "CENTRE" contaning the limiter of the liquid 83 : C velocity (2 composantes) 84 : C 85 : C IPC : CHPOINT "CENTRE" containing the pressure P ; 86 : C 87 : C IGRPC : CHPOINT "CENTRE" containing the gradient of pressure 88 : C (2 composantes) 89 : C 90 : C IALPC : CHPOINT "CENTRE" containig the limiter of pressure 91 : C gradient 92 : C 93 : C ITVC : CHPOINT "CENTRE" containing the vapour temperature Tv ; 94 : C 95 : C IGRTVC : CHPOINT "CENTRE" containing the gradient of vapour 96 : C temperature (2 composantes) 97 : C 98 : C IALTVC : CHPOINT "CENTRE" containig the limiter of vapour 99 : C temperature gradient 100 : C 101 : C ITLC : CHPOINT "CENTRE" containing the liquid temperature Tl ; 102 : C 103 : C IGRTLC : CHPOINT "CENTRE" containing the gradient of liquid 104 : C temperature (2 composantes) 105 : C 106 : C IALTLC : CHPOINT "CENTRE" containig the limiter of liquid 107 : C temperature gradient 108 : C 109 : C IRVC : CHPOINT "CENTRE" containing the vapour density 110 : C 111 : C IRLC : CHPOINT "CENTRE" containing the liquid density 112 : C 113 : C 3) 114 : C 115 : C DELTAT : REAL*8, encrement en temps pour calculer la prediction 116 : C 117 : C 118 : C SORTIES 119 : C 120 : C 121 : C IROF : MCHAML defini sur le MELEME de pointeur IFACEL, 122 : C contenant la masse volumique RHO 123 : C 124 : C IVITF : MCHAML defini sur le MELEME de pointeur IFACEL, 125 : C contenant la vitesse UN, UT dans le repaire local 126 : C (n,t) et defini sur le MELEME de pointeur IFACE, 127 : C contenant les cosinus directeurs du repere local 128 : C 129 : C IPF : MCHAML defini sur le MELEME de pointeur IFACEL, 130 : C contenant la pression P 131 : C 132 : C IGAMF : MCHAML defini sur le MELEME de pointeur IFACEL, 133 : C contenant le "gamma" du gaz 134 : C 135 : C LOGAN : anomalie detectee 136 : C 137 : C LOGNEG : (LOGICAL): si .TRUE. une pression ou une densité 138 : C negative a été detectée -> en interactif le 139 : C programme s'arrete en GIBIANE 140 : C (erreur stocké en MESERR et VALER) 141 : C 142 : C LOGBOR : (LOGICAL): si .TRUE. un gamma a ete detecte 143 : C dehor 1 et 3 (sa valeur stockée en MESERR et VALER; 144 : C en VAL1 et en VAL2 on stocke 1.0 et 3.0) 145 : C 146 : C MESERR 147 : C VALER 148 : C VAL1, 149 : C VAL2 : pour les messages d'erreur 150 : C 151 : C************************************************************************ 152 : C 153 : C HISTORIQUE (Anomalies et modifications éventuelles) 154 : C 155 : C HISTORIQUE : Créée le 02/05/02. 156 : C 157 : C************************************************************************ 158 : C 159 : C 160 : C ATTENTION: Cet programme marche si le MAILLAGE est convex; 161 : C si non il faut changer l'argoritme de calcul de 162 : C l'orientation des normales aux faces. 163 : C 164 : C La positivité n'est pas controlle parce que c'est déjà fait 165 : C dans l'operateur PRIM 166 : C 167 : C 168 : C************************************************************************ 169 : C 170 : C**** Variables de COOPTIO 171 : C 172 : INTEGER IPLLB, IERPER, IERMAX, IERR, INTERR 173 : & ,IOTER, IOLEC, IOIMP, IOCAR, IOACQ 174 : & ,IOPER, IOSGB, IOGRA, IOSAU, IORES 175 : & ,IECHO, IIMPI, IOSPI 176 : & ,IDIM, IFICLE, IPREFI 177 : C & ,MCOORD 178 : & ,IFOMOD, NIFOUR, IFOUR, NSDPGE, IONIVE 179 : & ,NGMAXY, IZROSF, ISOTYP, IOSCR,LTEXLU 180 : & ,NORINC,NORVAL,NORIND,NORVAD 181 : & ,NUCROU, IPSAUV 182 : C 183 : C**** Les variables 184 : C 185 : INTEGER ICEN, IFACE, IFACEL, INORM, 186 : & IALPH, IGRALP, IALALP, 187 : & IUVC, IGRUVC, IALUVC, 188 : & IULC, IGRULC, IALULC, 189 : & IPC, IGRPC, IALPC, 190 : & ITVC, IGRTVC, IALTVC, 191 : & ITLC, IGRTLC, IALTLC, 192 : & IRVC, IRLC, 193 : & IALPHF, IUVF, IULF, IPF, ITVF, ITLF, 194 : & IRVF, IRLF, 195 : & IGEOM, NFAC, NCEN, 196 : & N1PTEL, N1EL, N2PTEL, N2EL, N2, N1, N3, L1, NLCE, 197 : & NLCF, NGCF, NGCEG, NLCEG, NGCED, NLCED, NGCF1, 198 : & IDIMP1, INDCEL, 199 : & l, m, n 200 : REAL*8 VALER, VAL1, VAL2, XG, YG, ZG, XC, YC, ZC, XD, YD, ZD, 201 : & DELTAT 202 : & ,DXG, DYG, DZG, DXD, DYD, DZD,ORIENT 203 : & , CNX, CNY, CNZ, CTX, CTY, CTZ, CVX, CVY, CVZ 204 : & , alphag, uvxg, uvyg, uvzg, uvng, uvtg, uvvg 205 : & , ulxg, ulyg, ulzg, ulng, ultg, ulvg 206 : & , pg, Tvg, Tlg, rvg, rlg 207 : & , alphad, uvxd, uvyd, uvzd, uvnd, uvtd, uvvd 208 : & , ulxd, ulyd, ulzd, ulnd, ultd, ulvd 209 : & , pd, Tvd, Tld, rvd, rld 210 : & , VALCEL, DCEL, ALCEL 211 : & , alpha, uvx, uvy, uvz, ulx, uly, ulz, p, Tv, Tl 212 : & , rv, rl, drvdp, drvdhv, drldp, drldhl 213 : & , hv, hl 214 : & , A(10,10), B(10,10), C(10, 10), denom 215 : & , dVdx(10), dVdy(10), dVdz(10), dV(10), LIMITA(10) 216 : & , dVdVp(10,10), dVpdV(10,10), P1(10,10) 217 : & , dhvdp, dhvdTv, dhldp, dhldTl 218 : C & , DROX, DROY, DROZ, DUXX, DUXY, DUXZ, DUYX, DUYY, 219 : C & DUYZ, DUZX, DUZY, DUZZ, DPX, DPY, DPZ 220 : C & , DRO, DUX, DUY, DUZ, DP, ALPHA(5) 221 : CHARACTER*(40) MESERR 222 : CHARACTER*(8) TYPE 223 : LOGICAL LOGAN,LOGNEG, LOGBOR, LOGTEM, LOGI1, LOGALP 224 : C 225 : C**** perfect and stiffened gas eos 226 : C 227 : real*8 gamv, cpv, gaml, cpl, pil 228 : PARAMETER(gamv = 1.4D0, cpv = 1008.D0, 229 : & gaml = 2.8D0, cpl = 4186.D0, pil = 8.5D8) 230 : C 231 : C**** LOGALP = .TRUE. -> Prediction avec limiteur 232 : C 233 : PARAMETER(LOGALP = .TRUE.) 234 : C 235 : C**** Les Includes 236 : C 237 : -INC SMCOORD 238 : -INC CCOPTIO 239 : -INC SMCHPOI 240 : POINTEUR MPALPC.MPOVAL, MPGRAL.MPOVAL, 241 : & MPUVC.MPOVAL, MPGRUV.MPOVAL, 242 : & MPULC.MPOVAL, MPGRUL.MPOVAL, 243 : & MPPC.MPOVAL, MPGRP.MPOVAL, 244 : & MPTVC.MPOVAL, MPGRTV.MPOVAL, 245 : & MPTLC.MPOVAL, MPGRTL.MPOVAL, 246 : & MPRVC.MPOVAL, MPRLC.MPOVAL, 247 : & MPNORM.MPOVAL, 248 : & MPALPP.MPOVAL, MPUVP.MPOVAL, MPULP.MPOVAL, 249 : & MPPP.MPOVAL, MPTVP.MPOVAL, MPTLP.MPOVAL, 250 : & MPRVP.MPOVAL, MPRLP.MPOVAL 251 : -INC SMCHAML 252 : POINTEUR MLALP.MELVAL, MLP.MELVAL, 253 : & MLTV.MELVAL, MLTL.MELVAL, 254 : & MLRV.MELVAL, MLRL.MELVAL 255 : POINTEUR MLUVN.MELVAL, MLUVT.MELVAL, MLUVV.MELVAL, 256 : & MLULN.MELVAL, MLULT.MELVAL, MLULV.MELVAL, 257 : & MLVNX.MELVAL, MLVNY.MELVAL, MLVNZ.MELVAL, 258 : & MLVT1X.MELVAL, MLVT1Y.MELVAL, MLVT1Z.MELVAL, 259 : & MLVT2X.MELVAL, MLVT2Y.MELVAL, MLVT2Z.MELVAL, 260 : & MLLNX.MELVAL, MLLNY.MELVAL, MLLNZ.MELVAL, 261 : & MLLT1X.MELVAL, MLLT1Y.MELVAL, MLLT1Z.MELVAL, 262 : & MLLT2X.MELVAL, MLLT2Y.MELVAL, MLLT2Z.MELVAL 263 : C**** Some melemes for the LIMITAs 264 : POINTEUR MELLAL.MPOVAL, 265 : & MELLUV.MPOVAL, 266 : & MELLUL.MPOVAL, 267 : & MELLP.MPOVAL, 268 : & MELLTV.MPOVAL, 269 : & MELLTL.MPOVAL 270 : -INC SMLENTI 271 : -INC SMELEME 272 : C 273 : C 274 : C**** Initialisation des parametres d'erreur déjà faite, i.e. 275 : C 276 : C LOGNEG = .FALSE. 277 : C LOGBOR = .FALSE. 278 : C MESERR = ' ' 279 : C MOTERR(1:40) = MESERR(1:40) 280 : C VALER = 0.0D0 281 : C VAL1 = 0.0D0 282 : C VAL2 = 0.0D0 283 : C 284 : C 285 : C**** KRIPAD pour la correspondance global/local de centre 286 : C 286 : CALL KRIPAD(ICEN,MLENT1) 287 : C 288 : C**** MLENTI1 a MCORD.XCOORD(/1)/(IDIM+1) elements 289 : C 290 : C Si i est le numero global d'un noeud de ICEN, 291 : C MLENT1.LECT(i) contient sa position, i.e. 292 : C 293 : C I = numero global du noeud centre 294 : C MLENT1.LECT(i) = numero local du noeud centre 295 : C 296 : C MLENT1 déjà activé, i.e. 297 : C 298 : C SEGACT MLENT1 299 : C 300 : C**** Activation de CHPOINTs 301 : C 302 : C void fraction + grad + limiter 303 : C vapour velocity + grad + limiter 304 : C liquid velocity + grad + limiter 305 : C pressure + grad + limiter 306 : C vapour temperature + grad + limiter 307 : C liquid temperature + grad + limiter 308 : C vapour density 309 : C liquid density 310 : C 311 : C cosinus directeurs des normales aux surface 312 : C 312 : CALL LICHT(IALPH, MPALPC, TYPE, IGEOM) 312 : CALL LICHT(IGRALP, MPGRAL, TYPE, IGEOM) 312 : CALL LICHT(IUVC, MPUVC, TYPE, IGEOM) 312 : CALL LICHT(IGRUVC, MPGRUV, TYPE, IGEOM) 312 : CALL LICHT(IULC, MPULC, TYPE, IGEOM) 312 : CALL LICHT(IGRULC, MPGRUL, TYPE, IGEOM) 312 : CALL LICHT(IPC , MPPC , TYPE, IGEOM) 312 : CALL LICHT(IGRPC, MPGRP, TYPE, IGEOM) 312 : CALL LICHT(ITVC, MPTVC, TYPE, IGEOM) 312 : CALL LICHT(IGRTVC, MPGRTV, TYPE, IGEOM) 312 : CALL LICHT(ITLC, MPTLC, TYPE, IGEOM) 312 : CALL LICHT(IGRTLC, MPGRTL, TYPE, IGEOM) 312 : CALL LICHT(IRVC, MPRVC, TYPE, IGEOM) 312 : CALL LICHT(IRLC, MPRLC, TYPE, IGEOM) 312 : CALL LICHT(INORM, MPNORM, TYPE, IGEOM) 313 : C 314 : C**** Les MPOVALs 'Prediction' 315 : C 316 : IF(LOGTEM)THEN 317 : SEGINI, MPALPP = MPALPC 318 : SEGINI, MPUVP = MPUVC 319 : SEGINI, MPULP = MPULC 320 : SEGINI, MPPP = MPPC 321 : SEGINI, MPTVP = MPTVC 322 : SEGINI, MPTLP = MPTLC 323 : SEGINI, MPRVP = MPRVC 324 : SEGINI, MPRLP = MPRLC 325 : ELSE 326 : MPALPP = MPALPC 327 : MPUVP = MPUVC 328 : MPULP = MPULC 329 : MPPP = MPPC 330 : MPTVP = MPTVC 331 : MPTLP = MPTLC 332 : MPRVP = MPRVC 333 : MPRLP = MPRLC 334 : ENDIF 335 : C 336 : C**** Les Limiteurs 337 : C 337 : CALL LICHT(IALALP, MELLAL, TYPE, IGEOM) 337 : CALL LICHT(IALUVC, MELLUV, TYPE, IGEOM) 337 : CALL LICHT(IALULC, MELLUL, TYPE, IGEOM) 337 : CALL LICHT(IALPC , MELLP, TYPE, IGEOM) 337 : CALL LICHT(IALTVC, MELLTV, TYPE, IGEOM) 337 : CALL LICHT(IALTLC, MELLTL, TYPE, IGEOM) 338 : C 339 : C 340 : C**** MPOVAL sont déjà activés i.e.: 341 : C 342 : C SEGACT MPALPP 343 : C SEGACT MPGRAL 344 : C SEGACT IALALP 345 : C SEGACT MPUVC 346 : C SEGACT IALUVC 347 : C SEGACT MPGRUV 348 : C SEGACT MPULC 349 : C SEGACT MPGRUL 350 : C SEGACT IALULC 351 : C SEGACT MPPC 352 : C SEGACT MPGRP 353 : C SEGACT IALPC 354 : C SEGACT MPTVC 355 : C SEGACT MPGRTV 356 : C SEGACT IALTVC 357 : C SEGACT MPTLC 358 : C SEGACT MPGRTL 359 : C SEGACT IALTLC 360 : C SEGACT MPRVC 361 : C SEGACT MPRLC 362 : C SEGACT MPNORM 363 : C 364 : C**** Le MELEME FACEL 365 : C 366 : IPT1 = IFACEL 367 : IPT2 = IFACE 368 : SEGACT IPT1 369 : SEGACT IPT2 370 : NFAC = IPT1.NUM(/2) 371 : C 372 : C**** Creation de MCHAMLs contenant les etats gauche et droite, 373 : C 374 : C i.e.: 375 : C 376 : C vitesse + cosinus directors du repere local 377 : C densité 378 : C pression 379 : C gamma 380 : C 381 : C**** Cosinus directors du repere local et vitesse 382 : C 383 : C Les cosinus directeurs 384 : C 385 : N1 = 2 386 : N3 = 6 387 : L1 = 28 388 : SEGINI MCHEL1 389 : IUVF = MCHEL1 390 : MCHEL1.TITCHE = 'UV ' 391 : MCHEL1.IMACHE(1) = IFACE 392 : MCHEL1.IMACHE(2) = IFACEL 393 : MCHEL1.CONCHE(1) = '(n,t,v)in(x,y,z)' 394 : MCHEL1.CONCHE(2) = ' UV in (n,t,v) ' 395 : C 396 : C**** Valeurs des cosinus definies par respect au repair global, i.e. 397 : C 398 : MCHEL1.INFCHE(1,1) = 2 399 : MCHEL1.INFCHE(1,3) = NIFOUR 400 : MCHEL1.INFCHE(1,4) = 0 401 : MCHEL1.INFCHE(1,5) = 0 402 : MCHEL1.INFCHE(1,6) = 0 403 : MCHEL1.IFOCHE = IFOUR 404 : C 405 : C**** Valeurs de vitesse definies par respect au repair local, i.e. 406 : C 407 : MCHEL1.INFCHE(2,1) = 1 408 : MCHEL1.INFCHE(2,3) = NIFOUR 409 : MCHEL1.INFCHE(2,4) = 0 410 : MCHEL1.INFCHE(2,5) = 0 411 : MCHEL1.INFCHE(2,6) = 0 412 : C 413 : C**** Le cosinus directeurs 414 : C 415 : N1PTEL = 1 416 : N1EL = NFAC 417 : N2PTEL = 0 418 : N2EL = 0 419 : C 420 : C**** MCHAML a N2 composantes: 421 : C 422 : C cosinus directeurs du repere local (n,t1,t2) 423 : C 424 : C IDIM = 3 -> 9 composantes 425 : C 426 : N2 = 9 427 : SEGINI MCHAM1 428 : MCHEL1.ICHAML(1) = MCHAM1 429 : MCHAM1.NOMCHE(1) = 'NVX ' 430 : MCHAM1.NOMCHE(2) = 'NVY ' 431 : MCHAM1.NOMCHE(3) = 'NVZ ' 432 : MCHAM1.NOMCHE(4) = 'TVX ' 433 : MCHAM1.NOMCHE(5) = 'TVY ' 434 : MCHAM1.NOMCHE(6) = 'TVZ ' 435 : MCHAM1.NOMCHE(7) = 'VVX ' 436 : MCHAM1.NOMCHE(8) = 'VVY ' 437 : MCHAM1.NOMCHE(9) = 'VVZ ' 438 : MCHAM1.TYPCHE(1) = 'REAL*8 ' 439 : MCHAM1.TYPCHE(2) = 'REAL*8 ' 440 : MCHAM1.TYPCHE(3) = 'REAL*8 ' 441 : MCHAM1.TYPCHE(4) = 'REAL*8 ' 442 : MCHAM1.TYPCHE(5) = 'REAL*8 ' 443 : MCHAM1.TYPCHE(6) = 'REAL*8 ' 444 : MCHAM1.TYPCHE(7) = 'REAL*8 ' 445 : MCHAM1.TYPCHE(8) = 'REAL*8 ' 446 : MCHAM1.TYPCHE(9) = 'REAL*8 ' 447 : SEGINI MLVNX 448 : SEGINI MLVNY 449 : SEGINI MLVNZ 450 : SEGINI MLVT1X 451 : SEGINI MLVT1Y 452 : SEGINI MLVT1Z 453 : SEGINI MLVT2X 454 : SEGINI MLVT2Y 455 : SEGINI MLVT2Z 456 : MCHAM1.IELVAL(1) = MLVNX 457 : MCHAM1.IELVAL(2) = MLVNY 458 : MCHAM1.IELVAL(3) = MLVNZ 459 : MCHAM1.IELVAL(4) = MLVT1X 460 : MCHAM1.IELVAL(5) = MLVT1Y 461 : MCHAM1.IELVAL(6) = MLVT1Z 462 : MCHAM1.IELVAL(7) = MLVT2X 463 : MCHAM1.IELVAL(8) = MLVT2Y 464 : MCHAM1.IELVAL(9) = MLVT2Z 465 : SEGDES MCHAM1 466 : C 467 : C**** Vitesse 468 : C 469 : N1EL = NFAC 470 : N1PTEL = 3 471 : N2EL = 0 472 : N2PTEL = 0 473 : C 474 : C**** MCHAML a N2 composantes: 475 : C 476 : C IDIM = 3 -> 3 composantes 477 : C 478 : N2 = 3 479 : SEGINI MCHAM1 480 : MCHEL1.ICHAML(2) = MCHAM1 481 : SEGDES MCHEL1 482 : MCHAM1.NOMCHE(1) = 'UVN ' 483 : MCHAM1.NOMCHE(2) = 'UVT ' 484 : MCHAM1.NOMCHE(3) = 'UVV ' 485 : MCHAM1.TYPCHE(1) = 'REAL*8 ' 486 : MCHAM1.TYPCHE(2) = 'REAL*8 ' 487 : MCHAM1.TYPCHE(3) = 'REAL*8 ' 488 : SEGINI MLUVN 489 : SEGINI MLUVT 490 : SEGINI MLUVV 491 : MCHAM1.IELVAL(1) = MLUVN 492 : MCHAM1.IELVAL(2) = MLUVT 493 : MCHAM1.IELVAL(3) = MLUVV 494 : SEGDES MCHAM1 495 : C 496 : C**** Cosinus directors du repere local et vitesse 497 : C 498 : C Les cosinus directeurs 499 : C 500 : C LIQUID PHASE 501 : 502 : N1 = 2 503 : N3 = 6 504 : L1 = 28 505 : SEGINI MCHEL1 506 : IULF = MCHEL1 507 : MCHEL1.TITCHE = 'UL ' 508 : MCHEL1.IMACHE(1) = IFACE 509 : MCHEL1.IMACHE(2) = IFACEL 510 : MCHEL1.CONCHE(1) = '(n,t,v) in (x,y)' 511 : MCHEL1.CONCHE(2) = ' UL in (n,t,v) ' 512 : C 513 : C**** Valeurs des cosinus definies par respect au repair global, i.e. 514 : C 515 : MCHEL1.INFCHE(1,1) = 2 516 : MCHEL1.INFCHE(1,3) = NIFOUR 517 : MCHEL1.INFCHE(1,4) = 0 518 : MCHEL1.INFCHE(1,5) = 0 519 : MCHEL1.INFCHE(1,6) = 0 520 : MCHEL1.IFOCHE = IFOUR 521 : C 522 : C**** Valeurs de vitesse definies par respect au repair local, i.e. 523 : C 524 : MCHEL1.INFCHE(2,1) = 1 525 : MCHEL1.INFCHE(2,3) = NIFOUR 526 : MCHEL1.INFCHE(2,4) = 0 527 : MCHEL1.INFCHE(2,5) = 0 528 : MCHEL1.INFCHE(2,6) = 0 529 : C 530 : C**** Le cosinus directeurs 531 : C 532 : N1PTEL = 1 533 : N1EL = NFAC 534 : N2PTEL = 0 535 : N2EL = 0 536 : C 537 : C**** MCHAML a N2 composantes: 538 : C 539 : C cosinus directeurs du repere local (n,t1) 540 : C 541 : C IDIM = 3 -> 9 composantes 542 : C 543 : N2 = 9 544 : SEGINI MCHAM1 545 : MCHEL1.ICHAML(1) = MCHAM1 546 : MCHAM1.NOMCHE(1) = 'NLX ' 547 : MCHAM1.NOMCHE(2) = 'NLY ' 548 : MCHAM1.NOMCHE(3) = 'NLZ ' 549 : MCHAM1.NOMCHE(4) = 'TLX ' 550 : MCHAM1.NOMCHE(5) = 'TLY ' 551 : MCHAM1.NOMCHE(6) = 'TLZ ' 552 : MCHAM1.NOMCHE(7) = 'VLX ' 553 : MCHAM1.NOMCHE(8) = 'VLY ' 554 : MCHAM1.NOMCHE(9) = 'VLZ ' 555 : MCHAM1.TYPCHE(1) = 'REAL*8 ' 556 : MCHAM1.TYPCHE(2) = 'REAL*8 ' 557 : MCHAM1.TYPCHE(3) = 'REAL*8 ' 558 : MCHAM1.TYPCHE(4) = 'REAL*8 ' 559 : MCHAM1.TYPCHE(5) = 'REAL*8 ' 560 : MCHAM1.TYPCHE(6) = 'REAL*8 ' 561 : MCHAM1.TYPCHE(7) = 'REAL*8 ' 562 : MCHAM1.TYPCHE(8) = 'REAL*8 ' 563 : MCHAM1.TYPCHE(9) = 'REAL*8 ' 564 : SEGINI MLLNX 565 : SEGINI MLLNY 566 : SEGINI MLLNZ 567 : SEGINI MLLT1X 568 : SEGINI MLLT1Y 569 : SEGINI MLLT1Z 570 : SEGINI MLLT2X 571 : SEGINI MLLT2Y 572 : SEGINI MLLT2Z 573 : MCHAM1.IELVAL(1) = MLLNX 574 : MCHAM1.IELVAL(2) = MLLNY 575 : MCHAM1.IELVAL(3) = MLLNZ 576 : MCHAM1.IELVAL(4) = MLLT1X 577 : MCHAM1.IELVAL(5) = MLLT1Y 578 : MCHAM1.IELVAL(6) = MLLT1Z 579 : MCHAM1.IELVAL(7) = MLLT2X 580 : MCHAM1.IELVAL(8) = MLLT2Y 581 : MCHAM1.IELVAL(9) = MLLT2Z 582 : SEGDES MCHAM1 583 : C 584 : C**** Vitesse 585 : C 586 : N1EL = NFAC 587 : N1PTEL = 3 588 : N2EL = 0 589 : N2PTEL = 0 590 : C 591 : C**** MCHAML a N2 composantes: 592 : C 593 : C IDIM = 3 -> 3 composantes 594 : C 595 : N2 = 3 596 : SEGINI MCHAM1 597 : MCHEL1.ICHAML(2) = MCHAM1 598 : SEGDES MCHEL1 599 : MCHAM1.NOMCHE(1) = 'ULN ' 600 : MCHAM1.NOMCHE(2) = 'ULT ' 601 : MCHAM1.NOMCHE(3) = 'ULV ' 602 : MCHAM1.TYPCHE(1) = 'REAL*8 ' 603 : MCHAM1.TYPCHE(2) = 'REAL*8 ' 604 : MCHAM1.TYPCHE(3) = 'REAL*8 ' 605 : SEGINI MLULN 606 : SEGINI MLULT 607 : SEGINI MLULV 608 : MCHAM1.IELVAL(1) = MLULN 609 : MCHAM1.IELVAL(2) = MLULT 610 : MCHAM1.IELVAL(3) = MLULV 611 : SEGDES MCHAM1 612 : C 613 : C**** Void fraction 614 : C 615 : N1 = 1 616 : N3 = 6 617 : L1 = 15 618 : SEGINI MCHEL2 619 : IALPHF = MCHEL2 620 : MCHEL2.IMACHE(1) = IFACEL 621 : MCHEL2.TITCHE = 'ALPHA ' 622 : MCHEL2.CONCHE(1) = ' ' 623 : C 624 : C**** Valeurs independente du repére, i.e. 625 : C 626 : MCHEL2.INFCHE(1,1) = 0 627 : MCHEL2.INFCHE(1,3) = NIFOUR 628 : MCHEL2.INFCHE(1,4) = 0 629 : MCHEL2.INFCHE(1,5) = 0 630 : MCHEL2.INFCHE(1,6) = 0 631 : MCHEL2.IFOCHE = IFOUR 632 : N2 = 1 633 : SEGINI MCHAM1 634 : MCHEL2.ICHAML(1) = MCHAM1 635 : SEGDES MCHEL2 636 : MCHAM1.NOMCHE(1) = 'SCAL ' 637 : MCHAM1.TYPCHE(1) = 'REAL*8 ' 638 : SEGINI MLALP 639 : MCHAM1.IELVAL(1) = MLALP 640 : SEGDES MCHAM1 641 : C 642 : C**** Pressure 643 : C 644 : MCHEL1 = IALPHF 645 : SEGINI, MCHEL2 = MCHEL1 646 : IPF = MCHEL2 647 : MCHEL2.TITCHE = 'P ' 648 : C 649 : C**** MCHAM1 = MCHAML de la alpha 650 : C 651 : SEGINI, MCHAM2 = MCHAM1 652 : MCHEL2.ICHAML(1) = MCHAM2 653 : SEGDES MCHEL2 654 : SEGINI MLP 655 : MCHAM2.IELVAL(1) = MLP 656 : SEGDES MCHAM2 657 : C 658 : C**** Vapour temperature 659 : C 660 : MCHEL1 = IALPHF 661 : SEGINI, MCHEL2 = MCHEL1 662 : ITVF = MCHEL2 663 : MCHEL2.TITCHE = 'TV ' 664 : C 665 : C**** MCHAM1 = MCHAML de la alpha 666 : C 667 : SEGINI, MCHAM2 = MCHAM1 668 : MCHEL2.ICHAML(1) = MCHAM2 669 : SEGDES MCHEL2 670 : SEGINI MLTV 671 : MCHAM2.IELVAL(1) = MLTV 672 : SEGDES MCHAM2 673 : C 674 : C**** Liquid temperature 675 : C 676 : MCHEL1 = IALPHF 677 : SEGINI, MCHEL2 = MCHEL1 678 : ITLF = MCHEL2 679 : MCHEL2.TITCHE = 'TL ' 680 : C 681 : C**** MCHAM1 = MCHAML de la alpha 682 : C 683 : SEGINI, MCHAM2 = MCHAM1 684 : MCHEL2.ICHAML(1) = MCHAM2 685 : SEGDES MCHEL2 686 : SEGINI MLTL 687 : MCHAM2.IELVAL(1) = MLTL 688 : SEGDES MCHAM2 689 : C 690 : C**** Vapour density 691 : C 692 : MCHEL1 = IALPHF 693 : SEGINI, MCHEL2 = MCHEL1 694 : IRVF = MCHEL2 695 : MCHEL2.TITCHE = 'RV ' 696 : C 697 : C**** MCHAM1 = MCHAML de la alpha 698 : C 699 : SEGINI, MCHAM2 = MCHAM1 700 : MCHEL2.ICHAML(1) = MCHAM2 701 : SEGDES MCHEL2 702 : SEGINI MLRV 703 : MCHAM2.IELVAL(1) = MLRV 704 : SEGDES MCHAM2 705 : C 706 : C**** Liquid density 707 : C 708 : MCHEL1 = IALPHF 709 : SEGINI, MCHEL2 = MCHEL1 710 : IRLF = MCHEL2 711 : MCHEL2.TITCHE = 'RL ' 712 : C 713 : C**** MCHAM1 = MCHAML de la alpha 714 : C 715 : SEGINI, MCHAM2 = MCHAM1 716 : MCHEL2.ICHAML(1) = MCHAM2 717 : SEGDES MCHEL2 718 : SEGINI MLRL 719 : MCHAM2.IELVAL(1) = MLRL 720 : SEGDES MCHAM2 721 : C 722 : C**** Recapitulatif: les MELVALs et les MPOVALs actives 723 : C 724 : C MLVNX, MLVNY, MLVNZ 725 : C MLVTX, MLVTY, MLVTZ 726 : C MLVVX, MLVVY, MLVVZ 727 : C 728 : C MLLNX, MLLNY, MLLNZ 729 : C MLLTX, MLLTY, MLLTZ 730 : C MLLVX, MLLVY, MLLVZ 731 : C 732 : C MLUVN, MLUVT, MLUVV -> vapour velocities 733 : C 734 : C MLULN, MLULT, MLULV -> liquid velocities 735 : C 736 : C MLALP -> void fraction 737 : C 738 : C MLP -> pressure 739 : C 740 : C MLTV -> vapour temperature 741 : C 742 : C MLTL -> liquid temperature 743 : C 744 : C MLRV -> vapour density 745 : C 746 : C MLRL -> liquid density 747 : C**** 748 : C MPALP -> void fraction 749 : C 750 : C MPUVC -> vapour velocity 751 : C 752 : C MPULC -> liquid velocity 753 : C 754 : C MPPC -> pressure 755 : C 756 : C MPTVC -> vapour temperature 757 : C 758 : C MPTLC -> liquid temperature 759 : C 760 : C MPRVC -> vapour density 761 : C 762 : C MPRLC -> liquid density 763 : C 764 : C MPNORM -> normales aux faces 765 : C 766 : C 767 : C*********************************************************************** 768 : C********* PREDICTION ************************************************** 769 : C*********************************************************************** 770 : C 771 : IF(LOGTEM)THEN 772 : C 773 : IPT3 = ICEN 774 : SEGACT IPT3 775 : NCEN = IPT3.NUM(/2) 776 : DO NLCE = 1, NCEN 777 : IF(LOGALP)THEN 778 : C 779 : C************* Gradients limités 780 : C 781 : LIMITA(1) = MELLAL.VPOCHA(NLCE,1) 782 : LIMITA(2) = MELLUV.VPOCHA(NLCE,1) 783 : LIMITA(3) = MELLUV.VPOCHA(NLCE,2) 784 : LIMITA(4) = MELLUV.VPOCHA(NLCE,3) 785 : LIMITA(5) = MELLUL.VPOCHA(NLCE,1) 786 : LIMITA(6) = MELLUL.VPOCHA(NLCE,2) 787 : LIMITA(7) = MELLUL.VPOCHA(NLCE,3) 788 : LIMITA(8) = MELLP.VPOCHA(NLCE,1) 789 : LIMITA(9) = MELLTV.VPOCHA(NLCE,1) 790 : LIMITA(10) = MELLTL.VPOCHA(NLCE,1) 791 : 792 : ELSE 793 : C 794 : C************* Gradients non limités 795 : C 796 : LIMITA(1) = 1.D0 797 : LIMITA(2) = 1.D0 798 : LIMITA(3) = 1.D0 799 : LIMITA(4) = 1.D0 800 : LIMITA(5) = 1.D0 801 : LIMITA(6) = 1.D0 802 : LIMITA(7) = 1.D0 803 : LIMITA(8) = 1.D0 804 : LIMITA(9) = 1.D0 805 : LIMITA(10) = 1.D0 806 : 807 : ENDIF 808 : 809 : alpha = MPALPP.VPOCHA(NLCE,1) 810 : uvx = MPUVP.VPOCHA(NLCE,1) 811 : uvy = MPUVP.VPOCHA(NLCE,2) 812 : uvz = MPUVP.VPOCHA(NLCE,3) 813 : ulx = MPULP.VPOCHA(NLCE,1) 814 : uly = MPULP.VPOCHA(NLCE,2) 815 : ulz = MPULP.VPOCHA(NLCE,3) 816 : p = MPPP.VPOCHA(NLCE,1) 817 : Tv = MPTVP.VPOCHA(NLCE,1) 818 : Tl = MPTLP.VPOCHA(NLCE,1) 819 : rv = MPRVP.VPOCHA(NLCE,1) 820 : rl = MPRLP.VPOCHA(NLCE,1) 821 : C 822 : C************* Value of the enthalpies and the derivatives of 823 : C the densities respect pressure and enthalpies, 824 : C for the moment only for perfect and stiffened gases. 825 : hv = cpv*Tv 826 : hl = cpl*Tl 827 : drvdp = rv/p 828 : drvdhv = -rv/hv 829 : drldp = rl/(p + pil) 830 : drldhl = -rl/hl 831 : dhvdp = 0.d0 832 : dhvdTv = cpv 833 : dhldp = 0.d0 834 : dhldTl = cpl 835 : 836 : denom = drldhl*(drvdhv*p + rv*(alpha*drvdp*p + rv - alpha*rv)) + 837 : & rl*(alpha*rl*(drvdhv + drvdp*rv) - (-1 + alpha)*drldp*(drvdhv 838 : $ *p + rv**2)) 839 : 840 : A(1,1) = (alpha*(drldhl*p + rl**2)*(drvdhv + drvdp*rv)*ulx - 841 : & (-1.d0 + alpha)*(drldhl + drldp*rl)*(drvdhv*p + rv**2)*uvx) 842 : $ /denom 843 : A(1,2) = -(((-1.d0 + alpha)*alpha*(drldhl + drldp*rl)*rv**2) 844 : $ /denom) 845 : A(1,3) = 0.d0 846 : A(1,4) = 0.d0 847 : A(1,5) = ((-1.d0 + alpha)*alpha*rl**2*(drvdhv + drvdp*rv))/denom 848 : A(1,6) = 0.d0 849 : A(1,7) = 0.d0 850 : A(1,8) = ((-1.d0 + alpha)*alpha*(drldhl + drldp*rl)*(drvdhv + 851 : $ drvdp*rv)*(ulx - uvx))/denom 852 : A(1,9) = 0.d0 853 : A(1,10) = 0.d0 854 : 855 : A(2,1) = p/(alpha*rv) 856 : A(2,2) = uvx 857 : A(2,3) = 0.d0 858 : A(2,4) = 0.d0 859 : A(2,5) = 0.d0 860 : A(2,6) = 0.d0 861 : A(2,7) = 0.d0 862 : A(2,8) = 1.d0/rv 863 : A(2,9) = 0.d0 864 : A(2,10) = 0.d0 865 : 866 : A(3,1) = 0.d0 867 : A(3,2) = 0.d0 868 : A(3,3) = uvx 869 : A(3,4) = 0.d0 870 : A(3,5) = 0.d0 871 : A(3,6) = 0.d0 872 : A(3,7) = 0.d0 873 : A(3,8) = 0.d0 874 : A(3,9) = 0.d0 875 : A(3,10) = 0.d0 876 : 877 : A(4,1) = 0.d0 878 : A(4,2) = 0.d0 879 : A(4,3) = 0d0 880 : A(4,4) = uvx 881 : A(4,5) = 0.d0 882 : A(4,6) = 0.d0 883 : A(4,7) = 0.d0 884 : A(4,8) = 0.d0 885 : A(4,9) = 0.d0 886 : A(4,10) = 0.d0 887 : 888 : A(5,1) = p/(-rl + alpha*rl) 889 : A(5,2) = 0.d0 890 : A(5,3) = 0.d0 891 : A(5,4) = 0.d0 892 : A(5,5) = ulx 893 : A(5,6) = 0.d0 894 : A(5,7) = 0.d0 895 : A(5,8) = 1.d0/rl 896 : A(5,9) = 0.d0 897 : A(5,10) = 0.d0 898 : 899 : A(6,1) = 0.d0 900 : A(6,2) = 0.d0 901 : A(6,3) = 0.d0 902 : A(6,4) = 0.d0 903 : A(6,5) = 0.d0 904 : A(6,6) = ulx 905 : A(6,7) = 0.d0 906 : A(6,8) = 0.d0 907 : A(6,9) = 0.d0 908 : A(6,10) = 0.d0 909 : 910 : A(7,1) = 0.d0 911 : A(7,2) = 0.d0 912 : A(7,3) = 0.d0 913 : A(7,4) = 0.d0 914 : A(7,5) = 0.d0 915 : A(7,6) = 0.d0 916 : A(7,7) = ulx 917 : A(7,8) = 0.d0 918 : A(7,9) = 0.d0 919 : A(7,10) = 0.d0 920 : 921 : A(8,1) = -(((drldhl*p + rl**2)*(drvdhv*p + rv**2)*(ulx - uvx)) 922 : $ /denom) 923 : A(8,2) = (alpha*(drldhl*p + rl**2)*rv**2)/denom 924 : A(8,3) = 0.d0 925 : A(8,4) = 0.d0 926 : A(8,5) = -(((-1.d0 + alpha)*rl**2*(drvdhv*p + rv**2))/denom) 927 : A(8,6) = 0.d0 928 : A(8,7) = 0.d0 929 : A(8,8) = (-((-1.d0 + alpha)*(drldhl + drldp*rl)*(drvdhv*p + rv**2 930 : $ )*ulx) + alpha*(drldhl*p + rl**2)*(drvdhv + drvdp*rv)*uvx) 931 : $ /denom 932 : A(8,9) = 0.d0 933 : A(8,10) = 0.d0 934 : 935 : A(9,1) = ((drldhl*p + rl**2)*(drvdp*p - rv)*(ulx - uvx))/denom 936 : A(9,2) = ((drldhl*p + rl*(drldp*(p - alpha*p) + alpha*rl))*rv) 937 : $ /denom 938 : A(9,3) = 0.d0 939 : A(9,4) = 0.d0 940 : A(9,5) = ((-1.d0 + alpha)*rl**2*(drvdp*p - rv))/denom 941 : A(9,6) = 0.d0 942 : A(9,7) = 0.d0 943 : A(9,8) = ((-1.d0 + alpha)*(drldhl + drldp*rl)*(drvdp*p - rv)*(ulx 944 : $ - uvx))/denom 945 : A(9,9) = uvx 946 : A(9,10) = 0.d0 947 : 948 : A(10,1) = ((drldp*p - rl)*(drvdhv*p + rv**2)*(ulx - uvx))/denom 949 : A(10,2) = (alpha*(-(drldp*p) + rl)*rv**2)/denom 950 : A(10,3) = 0.d0 951 : A(10,4) = 0.d0 952 : A(10,5) = (rl*(drvdhv*p + rv*(alpha*drvdp*p + rv - alpha*rv))) 953 : $ /denom 954 : A(10,6) = 0.d0 955 : A(10,7) = 0.d0 956 : A(10,8) = (alpha*(drldp*p - rl)*(drvdhv + drvdp*rv)*(ulx - uvx)) 957 : $ /denom 958 : A(10,9) = 0.d0 959 : A(10,10) = ulx 960 : 961 : B(1,1) = (alpha*(drldhl*p + rl**2)*(drvdhv + drvdp*rv)*uly - (-1 962 : $ .d0 + alpha)*(drldhl + drldp*rl)*(drvdhv*p + rv**2)*uvy) 963 : $ /denom 964 : B(1,2) = 0.d0 965 : B(1,3) = -(((-1.d0 + alpha)*alpha*(drldhl + drldp*rl)*rv**2) 966 : $ /denom) 967 : B(1,4) = 0.d0 968 : B(1,5) = 0.d0 969 : B(1,6) = ((-1.d0 + alpha)*alpha*rl**2*(drvdhv + drvdp*rv))/denom 970 : B(1,7) = 0.d0 971 : B(1,8) = ((-1.d0 + alpha)*alpha*(drldhl + drldp*rl)*(drvdhv + 972 : $ drvdp*rv)*(uly - uvy))/denom 973 : B(1,9) = 0.d0 974 : B(1,10) = 0.d0 975 : 976 : B(2,1) = 0.d0 977 : B(2,2) = uvy 978 : B(2,3) = 0.d0 979 : B(2,4) = 0.d0 980 : B(2,5) = 0.d0 981 : B(2,6) = 0.d0 982 : B(2,7) = 0.d0 983 : B(2,8) = 0.d0 984 : B(2,9) = 0.d0 985 : B(2,10) = 0.d0 986 : 987 : B(3,1) = p/(alpha*rv) 988 : B(3,2) = 0.d0 989 : B(3,3) = uvy 990 : B(3,4) = 0.d0 991 : B(3,5) = 0.d0 992 : B(3,6) = 0.d0 993 : B(3,7) = 0.d0 994 : B(3,8) = 1.d0/rv 995 : B(3,9) = 0.d0 996 : B(3,10) = 0.d0 997 : 998 : B(4,1) = 0.d0 999 : B(4,2) = 0.d0 1000 : B(4,3) = 0.d0 1001 : B(4,4) = uvy 1002 : B(4,5) = 0.d0 1003 : B(4,6) = 0.d0 1004 : B(4,7) = 0.d0 1005 : B(4,8) = 0.d0 1006 : B(4,9) = 0.d0 1007 : B(4,10) = 0.d0 1008 : 1009 : B(5,1) = 0.d0 1010 : B(5,2) = 0.d0 1011 : B(5,3) = 0.d0 1012 : B(5,4) = 0.d0 1013 : B(5,5) = uly 1014 : B(5,6) = 0.d0 1015 : B(5,7) = 0.d0 1016 : B(5,8) = 0.d0 1017 : B(5,9) = 0.d0 1018 : B(5,10) = 0.d0 1019 : 1020 : B(6,1) = p/(-rl + alpha*rl) 1021 : B(6,2) = 0.d0 1022 : B(6,3) = 0.d0 1023 : B(6,4) = 0.d0 1024 : B(6,5) = 0.d0 1025 : B(6,6) = uly 1026 : B(6,7) = 0.d0 1027 : B(6,8) = 1.d0/rl 1028 : B(6,9) = 0.d0 1029 : B(6,10) = 0.d0 1030 : 1031 : B(7,1) = 0.d0 1032 : B(7,2) = 0.d0 1033 : B(7,3) = 0.d0 1034 : B(7,4) = 0.d0 1035 : B(7,5) = 0.d0 1036 : B(7,6) = 0.d0 1037 : B(7,7) = uly 1038 : B(7,8) = 0.d0 1039 : B(7,9) = 0.d0 1040 : B(7,10) = 0.d0 1041 : 1042 : B(8,1) = -(((drldhl*p + rl**2)*(drvdhv*p + rv**2)*(uly - uvy)) 1043 : $ /denom) 1044 : B(8,2) = 0.d0 1045 : B(8,3) = (alpha*(drldhl*p + rl**2)*rv**2)/denom 1046 : B(8,4) = 0.d0 1047 : B(8,5) = 0.d0 1048 : B(8,6) = -(((-1.d0 + alpha)*rl**2*(drvdhv*p + rv**2))/denom) 1049 : B(8,7) = 0.d0 1050 : B(8,8) = (-((-1.d0 + alpha)*(drldhl + drldp*rl)*(drvdhv*p + rv**2 1051 : $ )*uly) + alpha*(drldhl*p + rl**2)*(drvdhv + drvdp*rv)*uvy) 1052 : $ /denom 1053 : B(8,9) = 0.d0 1054 : B(8,10) = 0.d0 1055 : 1056 : B(9,1) = ((drldhl*p + rl**2)*(drvdp*p - rv)*(uly - uvy))/denom 1057 : B(9,2) = 0.d0 1058 : B(9,3) = ((drldhl*p + rl*(drldp*(p - alpha*p) + alpha*rl))*rv) 1059 : $ /denom 1060 : B(9,4) = 0.d0 1061 : B(9,5) = 0.d0 1062 : B(9,6) = ((-1.d0 + alpha)*rl**2*(drvdp*p - rv))/denom 1063 : B(9,7) = 0.d0 1064 : B(9,8) = ((-1.d0 + alpha)*(drldhl + drldp*rl)*(drvdp*p - rv)*(uly 1065 : $ - uvy))/denom 1066 : B(9,9) = uvy 1067 : B(9,10) = 0.d0 1068 : 1069 : B(10,1) = ((drldp*p - rl)*(drvdhv*p + rv**2)*(uly - uvy))/denom 1070 : B(10,2) = 0.d0 1071 : B(10,3) = (alpha*(-(drldp*p) + rl)*rv**2)/denom 1072 : B(10,4) = 0.d0 1073 : B(10,5) = 0.d0 1074 : B(10,6) = (rl*(drvdhv*p + rv*(alpha*drvdp*p + rv - alpha*rv))) 1075 : $ /denom 1076 : B(10,7) = 0.d0 1077 : B(10,8) = (alpha*(drldp*p - rl)*(drvdhv + drvdp*rv)*(uly - uvy)) 1078 : $ /denom 1079 : B(10,9) = 0.d0 1080 : B(10,10) = uly 1081 : 1082 : C(1,1) = (alpha*(drldhl*p + rl**2)*(drvdhv + drvdp*rv)*ulz - (-1 1083 : $ .d0 + alpha)*(drldhl + drldp*rl)*(drvdhv*p + rv**2)*uvz) 1084 : $ /denom 1085 : C(1,2) = 0.d0 1086 : C(1,3) = 0.d0 1087 : C(1,4) = -(((-1.d0 + alpha)*alpha*(drldhl + drldp*rl)*rv**2) 1088 : $ /denom) 1089 : C(1,5) = 0.d0 1090 : C(1,6) = 0.d0 1091 : C(1,7) = ((-1.d0 + alpha)*alpha*rl**2*(drvdhv + drvdp*rv))/denom 1092 : C(1,8) = ((-1.d0 + alpha)*alpha*(drldhl + drldp*rl)*(drvdhv + 1093 : $ drvdp*rv)*(ulz - uvz))/denom 1094 : C(1,9) = 0.d0 1095 : C(1,10) = 0.d0 1096 : 1097 : C(2,1) = 0.d0 1098 : C(2,2) = uvz 1099 : C(2,3) = 0.d0 1100 : C(2,4) = 0.d0 1101 : C(2,5) = 0.d0 1102 : C(2,6) = 0.d0 1103 : C(2,7) = 0.d0 1104 : C(2,8) = 0.d0 1105 : C(2,9) = 0.d0 1106 : C(2,10) = 0.d0 1107 : 1108 : C(3,1) = 0.d0 1109 : C(3,2) = 0.d0 1110 : C(3,3) = uvz 1111 : C(3,4) = 0.d0 1112 : C(3,5) = 0.d0 1113 : C(3,6) = 0.d0 1114 : C(3,7) = 0.d0 1115 : C(3,8) = 0.d0 1116 : C(3,9) = 0.d0 1117 : C(3,10) = 0.d0 1118 : 1119 : C(4,1) = p/(alpha*rv) 1120 : C(4,2) = 0.d0 1121 : C(4,3) = 0.d0 1122 : C(4,4) = uvz 1123 : C(4,5) = 0.d0 1124 : C(4,6) = 0.d0 1125 : C(4,7) = 0.d0 1126 : C(4,8) = 1.d0/rv 1127 : C(4,9) = 0.d0 1128 : C(4,10) = 0.d0 1129 : 1130 : C(5,1) = 0.d0 1131 : C(5,2) = 0.d0 1132 : C(5,3) = 0.d0 1133 : C(5,4) = 0.d0 1134 : C(5,5) = ulz 1135 : C(5,6) = 0.d0 1136 : C(5,7) = 0.d0 1137 : C(5,8) = 0.d0 1138 : C(5,9) = 0.d0 1139 : C(5,10) = 0.d0 1140 : 1141 : C(6,1) = 0.d0 1142 : C(6,2) = 0.d0 1143 : C(6,3) = 0.d0 1144 : C(6,4) = 0.d0 1145 : C(6,5) = 0.d0 1146 : C(6,6) = ulz 1147 : C(6,7) = 0.d0 1148 : C(6,8) = 0.d0 1149 : C(6,9) = 0.d0 1150 : C(6,10) = 0.d0 1151 : 1152 : C(7,1) = p/(-rl + alpha*rl) 1153 : C(7,2) = 0.d0 1154 : C(7,3) = 0.d0 1155 : C(7,4) = 0.d0 1156 : C(7,5) = 0.d0 1157 : C(7,6) = 0.d0 1158 : C(7,7) = ulz 1159 : C(7,8) = 1.d0/rl 1160 : C(7,9) = 0.d0 1161 : C(7,10) = 0.d0 1162 : 1163 : C(8,1) = -(((drldhl*p + rl**2)*(drvdhv*p + rv**2)*(ulz - uvz)) 1164 : $ /denom) 1165 : C(8,2) = 0.d0 1166 : C(8,3) = 0.d0 1167 : C(8,4) = (alpha*(drldhl*p + rl**2)*rv**2)/denom 1168 : C(8,5) = 0.d0 1169 : C(8,6) = 0.d0 1170 : C(8,7) = -(((-1.d0 + alpha)*rl**2*(drvdhv*p + rv**2))/denom) 1171 : C(8,8) = (-((-1.d0 + alpha)*(drldhl + drldp*rl)*(drvdhv*p + rv**2 1172 : $ )*ulz) + alpha*(drldhl*p + rl**2)*(drvdhv + drvdp*rv)*uvz) 1173 : $ /denom 1174 : C(8,9) = 0.d0 1175 : C(8,10) = 0.d0 1176 : 1177 : C(9,1) = ((drldhl*p + rl**2)*(drvdp*p - rv)*(ulz - uvz))/denom 1178 : C(9,2) = 0.d0 1179 : C(9,3) = 0.d0 1180 : C(9,4) = ((drldhl*p + rl*(drldp*(p - alpha*p) + alpha*rl))*rv) 1181 : $ /denom 1182 : C(9,5) = 0.d0 1183 : C(9,6) = 0.d0 1184 : C(9,7) = ((-1.d0 + alpha)*rl**2*(drvdp*p - rv))/denom 1185 : C(9,8) = ((-1.d0 + alpha)*(drldhl + drldp*rl)*(drvdp*p - rv)*(ulz 1186 : $ - uvz))/denom 1187 : C(9,9) = uvz 1188 : C(9,10) = 0.d0 1189 : 1190 : C(10,1) = ((drldp*p - rl)*(drvdhv*p + rv**2)*(ulz - uvz))/denom 1191 : C(10,2) = 0.d0 1192 : C(10,3) = 0.d0 1193 : C(10,4) = (alpha*(-(drldp*p) + rl)*rv**2)/denom 1194 : C(10,5) = 0.d0 1195 : C(10,6) = 0.d0 1196 : C(10,7) = (rl*(drvdhv*p + rv*(alpha*drvdp*p + rv - alpha*rv))) 1197 : $ /denom 1198 : C(10,8) = (alpha*(drldp*p - rl)*(drvdhv + drvdp*rv)*(ulz - uvz)) 1199 : $ /denom 1200 : C(10,9) = 0.d0 1201 : C(10,10) = ulz 1202 : C 1203 : C dVpdV 1204 : C 1205 : dVpdV(1,1) = 1.d0 1206 : dVpdV(1,2) = 0.d0 1207 : dVpdV(1,3) = 0.d0 1208 : dVpdV(1,4) = 0.d0 1209 : dVpdV(1,5) = 0.d0 1210 : dVpdV(1,6) = 0.d0 1211 : dVpdV(1,7) = 0.d0 1212 : dVpdV(1,8) = 0.d0 1213 : dVpdV(1,9) = 0.d0 1214 : dVpdV(1,10) = 0.d0 1215 : 1216 : dVpdV(2,1) = 0.d0 1217 : dVpdV(2,2) = 1.d0 1218 : dVpdV(2,3) = 0.d0 1219 : dVpdV(2,4) = 0.d0 1220 : dVpdV(2,5) = 0.d0 1221 : dVpdV(2,6) = 0.d0 1222 : dVpdV(2,7) = 0.d0 1223 : dVpdV(2,8) = 0.d0 1224 : dVpdV(2,9) = 0.d0 1225 : dVpdV(2,10) = 0.d0 1226 : 1227 : dVpdV(3,1) = 0.d0 1228 : dVpdV(3,2) = 0.d0 1229 : dVpdV(3,3) = 1.d0 1230 : dVpdV(3,4) = 0.d0 1231 : dVpdV(3,5) = 0.d0 1232 : dVpdV(3,6) = 0.d0 1233 : dVpdV(3,7) = 0.d0 1234 : dVpdV(3,8) = 0.d0 1235 : dVpdV(3,9) = 0.d0 1236 : dVpdV(3,10) = 0.d0 1237 : 1238 : dVpdV(4,1) = 0.d0 1239 : dVpdV(4,2) = 0.d0 1240 : dVpdV(4,3) = 0.d0 1241 : dVpdV(4,4) = 1.d0 1242 : dVpdV(4,5) = 0.d0 1243 : dVpdV(4,6) = 0.d0 1244 : dVpdV(4,7) = 0.d0 1245 : dVpdV(4,8) = 0.d0 1246 : dVpdV(4,9) = 0.d0 1247 : dVpdV(4,10) = 0.d0 1248 : 1249 : dVpdV(5,1) = 0.d0 1250 : dVpdV(5,2) = 0.d0 1251 : dVpdV(5,3) = 0.d0 1252 : dVpdV(5,4) = 0.d0 1253 : dVpdV(5,5) = 1.d0 1254 : dVpdV(5,6) = 0.d0 1255 : dVpdV(5,7) = 0.d0 1256 : dVpdV(5,8) = 0.d0 1257 : dVpdV(5,9) = 0.d0 1258 : dVpdV(5,10) = 0.d0 1259 : 1260 : dVpdV(6,1) = 0.d0 1261 : dVpdV(6,2) = 0.d0 1262 : dVpdV(6,3) = 0.d0 1263 : dVpdV(6,4) = 0.d0 1264 : dVpdV(6,5) = 0.d0 1265 : dVpdV(6,6) = 1.d0 1266 : dVpdV(6,7) = 0.d0 1267 : dVpdV(6,8) = 0.d0 1268 : dVpdV(6,9) = 0.d0 1269 : dVpdV(6,10) = 0.d0 1270 : 1271 : dVpdV(7,1) = 0.d0 1272 : dVpdV(7,2) = 0.d0 1273 : dVpdV(7,3) = 0.d0 1274 : dVpdV(7,4) = 0.d0 1275 : dVpdV(7,5) = 0.d0 1276 : dVpdV(7,6) = 0.d0 1277 : dVpdV(7,7) = 1.d0 1278 : dVpdV(7,8) = 0.d0 1279 : dVpdV(7,9) = 0.d0 1280 : dVpdV(7,10) = 0.d0 1281 : 1282 : dVpdV(8,1) = 0.d0 1283 : dVpdV(8,2) = 0.d0 1284 : dVpdV(8,3) = 0.d0 1285 : dVpdV(8,4) = 0.d0 1286 : dVpdV(8,5) = 0.d0 1287 : dVpdV(8,6) = 0.d0 1288 : dVpdV(8,7) = 0.d0 1289 : dVpdV(8,8) = 1.d0 1290 : dVpdV(8,9) = 0.d0 1291 : dVpdV(8,10) = 0.d0 1292 : 1293 : dVpdV(9,1) = 0.d0 1294 : dVpdV(9,2) = 0.d0 1295 : dVpdV(9,3) = 0.d0 1296 : dVpdV(9,4) = 0.d0 1297 : dVpdV(9,5) = 0.d0 1298 : dVpdV(9,6) = 0.d0 1299 : dVpdV(9,7) = 0.d0 1300 : dVpdV(9,8) = dhvdp 1301 : dVpdV(9,9) = dhvdTv 1302 : dVpdV(9,10) = 0.d0 1303 : 1304 : dVpdV(10,1) = 0.d0 1305 : dVpdV(10,2) = 0.d0 1306 : dVpdV(10,3) = 0.d0 1307 : dVpdV(10,4) = 0.d0 1308 : dVpdV(10,5) = 0.d0 1309 : dVpdV(10,6) = 0.d0 1310 : dVpdV(10,7) = 0.d0 1311 : dVpdV(10,8) = dhldp 1312 : dVpdV(10,9) = 0.d0 1313 : dVpdV(10,10) = dhldTl 1314 : C 1315 : C dVdVp 1316 : C 1317 : dVdVp(1,1) = 1.d0 1318 : dVdVp(1,2) = 0.d0 1319 : dVdVp(1,3) = 0.d0 1320 : dVdVp(1,4) = 0.d0 1321 : dVdVp(1,5) = 0.d0 1322 : dVdVp(1,6) = 0.d0 1323 : dVdVp(1,7) = 0.d0 1324 : dVdVp(1,8) = 0.d0 1325 : dVdVp(1,9) = 0.d0 1326 : dVdVp(1,10) = 0.d0 1327 : 1328 : dVdVp(2,1) = 0.d0 1329 : dVdVp(2,2) = 1.d0 1330 : dVdVp(2,3) = 0.d0 1331 : dVdVp(2,4) = 0.d0 1332 : dVdVp(2,5) = 0.d0 1333 : dVdVp(2,6) = 0.d0 1334 : dVdVp(2,7) = 0.d0 1335 : dVdVp(2,8) = 0.d0 1336 : dVdVp(2,9) = 0.d0 1337 : dVdVp(2,10) = 0.d0 1338 : 1339 : dVdVp(3,1) = 0.d0 1340 : dVdVp(3,2) = 0.d0 1341 : dVdVp(3,3) = 1.d0 1342 : dVdVp(3,4) = 0.d0 1343 : dVdVp(3,5) = 0.d0 1344 : dVdVp(3,6) = 0.d0 1345 : dVdVp(3,7) = 0.d0 1346 : dVdVp(3,8) = 0.d0 1347 : dVdVp(3,9) = 0.d0 1348 : dVdVp(3,10) = 0.d0 1349 : 1350 : dVdVp(4,1) = 0.d0 1351 : dVdVp(4,2) = 0.d0 1352 : dVdVp(4,3) = 0.d0 1353 : dVdVp(4,4) = 1.d0 1354 : dVdVp(4,5) = 0.d0 1355 : dVdVp(4,6) = 0.d0 1356 : dVdVp(4,7) = 0.d0 1357 : dVdVp(4,8) = 0.d0 1358 : dVdVp(4,9) = 0.d0 1359 : dVdVp(4,10) = 0.d0 1360 : 1361 : dVdVp(5,1) = 0.d0 1362 : dVdVp(5,2) = 0.d0 1363 : dVdVp(5,3) = 0.d0 1364 : dVdVp(5,4) = 0.d0 1365 : dVdVp(5,5) = 1.d0 1366 : dVdVp(5,6) = 0.d0 1367 : dVdVp(5,7) = 0.d0 1368 : dVdVp(5,8) = 0.d0 1369 : dVdVp(5,9) = 0.d0 1370 : dVdVp(5,10) = 0.d0 1371 : 1372 : dVdVp(6,1) = 0.d0 1373 : dVdVp(6,2) = 0.d0 1374 : dVdVp(6,3) = 0.d0 1375 : dVdVp(6,4) = 0.d0 1376 : dVdVp(6,5) = 0.d0 1377 : dVdVp(6,6) = 1.d0 1378 : dVdVp(6,7) = 0.d0 1379 : dVdVp(6,8) = 0.d0 1380 : dVdVp(6,9) = 0.d0 1381 : dVdVp(6,10) = 0.d0 1382 : 1383 : dVdVp(7,1) = 0.d0 1384 : dVdVp(7,2) = 0.d0 1385 : dVdVp(7,3) = 0.d0 1386 : dVdVp(7,4) = 0.d0 1387 : dVdVp(7,5) = 0.d0 1388 : dVdVp(7,6) = 0.d0 1389 : dVdVp(7,7) = 1.d0 1390 : dVdVp(7,8) = 0.d0 1391 : dVdVp(7,9) = 0.d0 1392 : dVdVp(7,10) = 0.d0 1393 : 1394 : dVdVp(8,1) = 0.d0 1395 : dVdVp(8,2) = 0.d0 1396 : dVdVp(8,3) = 0.d0 1397 : dVdVp(8,4) = 0.d0 1398 : dVdVp(8,5) = 0.d0 1399 : dVdVp(8,6) = 0.d0 1400 : dVdVp(8,7) = 0.d0 1401 : dVdVp(8,8) = 1.d0 1402 : dVdVp(8,9) = 0.d0 1403 : dVdVp(8,10) = 0.d0 1404 : 1405 : dVdVp(9,1) = 0.d0 1406 : dVdVp(9,2) = 0.d0 1407 : dVdVp(9,3) = 0.d0 1408 : dVdVp(9,4) = 0.d0 1409 : dVdVp(9,5) = 0.d0 1410 : dVdVp(9,6) = 0.d0 1411 : dVdVp(9,7) = 0.d0 1412 : dVdVp(9,8) = -dhvdp/dhvdTv 1413 : dVdVp(9,9) = 1.d0/dhvdTv 1414 : dVdVp(9,10) = 0.d0 1415 : 1416 : dVdVp(10,1) = 0.d0 1417 : dVdVp(10,2) = 0.d0 1418 : dVdVp(10,3) = 0.d0 1419 : dVdVp(10,4) = 0.d0 1420 : dVdVp(10,5) = 0.d0 1421 : dVdVp(10,6) = 0.d0 1422 : dVdVp(10,7) = 0.d0 1423 : dVdVp(10,8) = dhldp*dhvdp/(dhldTl*dhvdTv) 1424 : dVdVp(10,9) = -dhldp/(dhldTl*dhvdTv) 1425 : dVdVp(10,10) = 1.d0/dhldTl 1426 : C 1427 : C dVdVp*A*dVpdV 1428 : C 1429 : do l = 1, 10 1430 : do m = 1, 10 1431 : P1(l,m) = 0.d0 1432 : end do 1433 : end do 1434 : do l = 1, 10 1435 : do m = 1, 10 1436 : do n = 1, 10 1437 : P1(l,m) = P1(l,m) + dVdVp(l,n)*A(n,m) 1438 : end do 1439 : end do 1440 : end do 1441 : do l = 1, 10 1442 : do m = 1, 10 1443 : A(l,m) = 0.d0 1444 : end do 1445 : end do 1446 : do l = 1, 8 1447 : do m = 1, 10 1448 : do n = 1, 10 1449 : A(l,m) = A(l,m) + P1(l,n)*dVpdV(n,m) 1450 : end do 1451 : end do 1452 : end do 1453 : C 1454 : C dVdVp*B*dVpdV 1455 : C 1456 : do l = 1, 10 1457 : do m = 1, 10 1458 : P1(l,m) = 0.d0 1459 : end do 1460 : end do 1461 : do l = 1, 10 1462 : do m = 1, 10 1463 : do n = 1, 10 1464 : P1(l,m) = P1(l,m) + dVdVp(l,n)*B(n,m) 1465 : end do 1466 : end do 1467 : end do 1468 : do l = 1, 10 1469 : do m = 1, 10 1470 : B(l,m) = 0.d0 1471 : end do 1472 : end do 1473 : do l = 1, 10 1474 : do m = 1, 10 1475 : do n = 1, 10 1476 : B(l,m) = B(l,m) + P1(l,n)*dVpdV(n,m) 1477 : end do 1478 : end do 1479 : end do 1480 : C 1481 : C dVdVp*C*dVpdV 1482 : C 1483 : do l = 1, 10 1484 : do m = 1, 10 1485 : P1(l,m) = 0.d0 1486 : end do 1487 : end do 1488 : do l = 1, 10 1489 : do m = 1, 10 1490 : do n = 1, 10 1491 : P1(l,m) = P1(l,m) + dVdVp(l,n)*C(n,m) 1492 : end do 1493 : end do 1494 : end do 1495 : do l = 1, 10 1496 : do m = 1, 10 1497 : C(l,m) = 0.d0 1498 : end do 1499 : end do 1500 : do l = 1, 10 1501 : do m = 1, 10 1502 : do n = 1, 10 1503 : C(l,m) = C(l,m) + P1(l,n)*dVpdV(n,m) 1504 : end do 1505 : end do 1506 : end do 1507 : C 1508 : dVdx(1) = MPGRAL.VPOCHA(NLCE,1)*LIMITA(1) 1509 : C duvxdx 1510 : dVdx(2) = MPGRUV.VPOCHA(NLCE,1)*LIMITA(2) 1511 : C duvydx 1512 : dVdx(3) = MPGRUV.VPOCHA(NLCE,4)*LIMITA(3) 1513 : C duvzdx 1514 : dVdx(4) = MPGRUV.VPOCHA(NLCE,7)*LIMITA(4) 1515 : C the same with the liquid velocity 1516 : dVdx(5) = MPGRUL.VPOCHA(NLCE,1)*LIMITA(5) 1517 : dVdx(6) = MPGRUL.VPOCHA(NLCE,4)*LIMITA(6) 1518 : dVdx(7) = MPGRUL.VPOCHA(NLCE,7)*LIMITA(7) 1519 : dVdx(8) = MPGRP.VPOCHA(NLCE,1)*LIMITA(8) 1520 : dVdx(9) = MPGRTV.VPOCHA(NLCE,1)*LIMITA(9) 1521 : dVdx(10) = MPGRTL.VPOCHA(NLCE,1)*LIMITA(10) 1522 : 1523 : dVdy(1) = MPGRAL.VPOCHA(NLCE,2)*LIMITA(1) 1524 : dVdy(2) = MPGRUV.VPOCHA(NLCE,2)*LIMITA(2) 1525 : dVdy(3) = MPGRUV.VPOCHA(NLCE,5)*LIMITA(3) 1526 : dVdy(4) = MPGRUV.VPOCHA(NLCE,8)*LIMITA(4) 1527 : dVdy(5) = MPGRUL.VPOCHA(NLCE,2)*LIMITA(5) 1528 : dVdy(6) = MPGRUL.VPOCHA(NLCE,5)*LIMITA(6) 1529 : dVdy(7) = MPGRUL.VPOCHA(NLCE,8)*LIMITA(7) 1530 : dVdy(8) = MPGRP.VPOCHA(NLCE,2)*LIMITA(8) 1531 : dVdy(9) = MPGRTV.VPOCHA(NLCE,2)*LIMITA(9) 1532 : dVdy(10) = MPGRTL.VPOCHA(NLCE,2)*LIMITA(10) 1533 : 1534 : dVdz(1) = MPGRAL.VPOCHA(NLCE,3)*LIMITA(1) 1535 : dVdz(2) = MPGRUV.VPOCHA(NLCE,3)*LIMITA(2) 1536 : dVdz(3) = MPGRUV.VPOCHA(NLCE,6)*LIMITA(3) 1537 : dVdz(4) = MPGRUV.VPOCHA(NLCE,9)*LIMITA(4) 1538 : dVdz(5) = MPGRUL.VPOCHA(NLCE,3)*LIMITA(5) 1539 : dVdz(6) = MPGRUL.VPOCHA(NLCE,6)*LIMITA(6) 1540 : dVdz(7) = MPGRUL.VPOCHA(NLCE,9)*LIMITA(7) 1541 : dVdz(8) = MPGRP.VPOCHA(NLCE,3)*LIMITA(8) 1542 : dVdz(9) = MPGRTV.VPOCHA(NLCE,3)*LIMITA(9) 1543 : dVdz(10) = MPGRTL.VPOCHA(NLCE,3)*LIMITA(10) 1544 : 1545 : C No te olvides de afectar de la correspondiente cp las 1546 : C dVdx y dVdy correspondientes a las entalpias, no vayas 1547 : C a meter la pata con eso 1548 : do m = 1, 10 1549 : dV(m) = 0.d0 1550 : end do 1551 : do m = 1, 10 1552 : do n = 1, 10 1553 : dV(m) = dV(m) + A(m,n)*dVdx(n) + B(m,n)*dVdy(n) 1554 : & + C(m,n)*dVdz(n) 1555 : end do 1556 : end do 1557 : C 1558 : C************* Prediction avec/sans gradients limités 1559 : C 1560 : MPALPP.VPOCHA(NLCE,1)= alpha - DELTAT*dV(1) 1561 : MPUVP.VPOCHA(NLCE,1) = uvx - DELTAT*dV(2) 1562 : MPUVP.VPOCHA(NLCE,2) = uvy - DELTAT*dV(3) 1563 : MPUVP.VPOCHA(NLCE,3) = uvz - DELTAT*dV(4) 1564 : MPULP.VPOCHA(NLCE,1) = ulx - DELTAT*dV(5) 1565 : MPULP.VPOCHA(NLCE,2) = uly - DELTAT*dV(6) 1566 : MPULP.VPOCHA(NLCE,3) = ulz - DELTAT*dV(7) 1567 : MPPP.VPOCHA(NLCE,1) = p - DELTAT*dV(8) 1568 : MPTVP.VPOCHA(NLCE,1) = Tv - DELTAT*dV(9) 1569 : MPTLP.VPOCHA(NLCE,1) = Tl - DELTAT*dV(10) 1570 : C 1571 : C************* If we only want 2nd order in space 1572 : C 1573 : C MPALPP.VPOCHA(NLCE,1)= alpha 1574 : C MPUVP.VPOCHA(NLCE,1) = uvx 1575 : C MPUVP.VPOCHA(NLCE,2) = uvy 1576 : C MPUVP.VPOCHA(NLCE,3) = uvz 1577 : C MPULP.VPOCHA(NLCE,1) = ulx 1578 : C MPULP.VPOCHA(NLCE,2) = uly 1579 : C MPULP.VPOCHA(NLCE,3) = ulz 1580 : C MPPP.VPOCHA(NLCE,1) = p 1581 : C MPTVP.VPOCHA(NLCE,1) = Tv 1582 : C MPTLP.VPOCHA(NLCE,1) = Tl 1583 : C 1584 : C************* Updating densities(perfect and stiffened gas eos 1585 : C 1586 : MPRVP.VPOCHA(NLCE,1) = gamv*MPPP.VPOCHA(NLCE,1)/ 1587 : & (cpv*MPTVP.VPOCHA(NLCE,1)*(gamv - 1.d0)) 1588 : MPRLP.VPOCHA(NLCE,1) = gaml*(MPPP.VPOCHA(NLCE,1) + pil)/ 1589 : & (cpl* MPTLP.VPOCHA(NLCE,1)*(gaml - 1.d0)) 1590 : 1591 : ENDDO 1592 : C 1593 : ENDIF 1594 : C 1595 : C 1596 : C*********************************************************************** 1597 : C********* CORRECTION ************************************************** 1598 : C*********************************************************************** 1599 : C 1600 : C**** Boucle sur le faces 1601 : C 1602 : IDIMP1 = IDIM + 1 1603 : DO NLCF = 1, NFAC 1604 : C 1605 : C******* NLCF = numero local du centre de face 1606 : C NGCF = numero global du centre de face 1607 : C NGCEG = numero global du centre ELT "gauche" 1608 : C NLCEG = numero local du centre ELT "gauche" 1609 : C NGCED = numero global du centre ELT "droite" 1610 : C NLCED = numero local du centre ELT "droite" 1611 : C 1612 : NGCEG = IPT1.NUM(1,NLCF) 1613 : NGCF = IPT1.NUM(2,NLCF) 1614 : NGCED = IPT1.NUM(3,NLCF) 1615 : NLCEG = MLENT1.LECT(NGCEG) 1616 : NLCED = MLENT1.LECT(NGCED) 1617 : C 1618 : C******* TEST: IPT2.NUM(1,NLCF) = IPT1.NUM(2,NLCF) 1619 : C 1620 : NGCF1 = IPT2.NUM(1,NLCF) 1621 : IF( NGCF1 .NE. NGCF) THEN 1622 : LOGAN = .TRUE. 1623 : MESERR(1:40) = 'PRET, subroutine pre62f.eso ' 1624 : GOTO 9999 1625 : ENDIF 1626 : C 1627 : C******* Cosinus directeurs des NORMALES aux faces 1628 : C 1629 : C On impose que les normales sont direct "Gauche" -> "Centre" 1630 : C 1631 : INDCEL = (NGCEG-1)*IDIMP1 1632 : XG = XCOOR(INDCEL+1) 1633 : YG = XCOOR(INDCEL+2) 1634 : ZG = XCOOR(INDCEL+3) 1635 : INDCEL = (NGCF-1)*IDIMP1 1636 : XC = XCOOR(INDCEL + 1) 1637 : YC = XCOOR(INDCEL + 2) 1638 : ZC = XCOOR(INDCEL+3) 1639 : INDCEL = (NGCED-1)*IDIMP1 1640 : XD = XCOOR(INDCEL+1) 1641 : YD = XCOOR(INDCEL+2) 1642 : ZD = XCOOR(INDCEL+3) 1643 : DXG = XC - XG 1644 : DYG = YC - YG 1645 : DZG = ZC - ZG 1646 : DXD = XC - XD 1647 : DYD = YC - YD 1648 : DZD = ZC - ZD 1649 : C 1650 : C******* On calcule le sign du pruduit scalare 1651 : C (Normales de Castem) * (vecteur "gauche" -> "centre") 1652 : C 1653 : CNX = MPNORM.VPOCHA(NLCF,7) 1654 : CNY = MPNORM.VPOCHA(NLCF,8) 1655 : CNZ = MPNORM.VPOCHA(NLCF,9) 1656 : ORIENT = CNX * DXG + CNY * DYG + CNZ * DZG 1657 : ORIENT = SIGN(1.0D0,ORIENT) 1658 : IF(ORIENT .NE. 1.0D0)THEN 1659 : LOGAN = .TRUE. 1660 : MESERR(1:30)= 1661 : & 'PRET , subroutine pre122.eso. ' 1662 : GOTO 9999 1663 : ENDIF 1664 : CNX = CNX * ORIENT 1665 : CNY = CNY * ORIENT 1666 : CNZ = CNZ * ORIENT 1667 : C 1668 : C********** Cosinus directeurs de tangente 1 1669 : C 1670 : CTX = MPNORM.VPOCHA(NLCF,1) * ORIENT 1671 : CTY = MPNORM.VPOCHA(NLCF,2) * ORIENT 1672 : CTZ = MPNORM.VPOCHA(NLCF,3) * ORIENT 1673 : C 1674 : C********** Cosinus directeurs de tangente 2 1675 : C 1676 : CVX = MPNORM.VPOCHA(NLCF,4) * ORIENT 1677 : CVY = MPNORM.VPOCHA(NLCF,5) * ORIENT 1678 : CVZ = MPNORM.VPOCHA(NLCF,6) * ORIENT 1679 : C 1680 : C 1681 : C******* Les autres MELVALs 1682 : C 1683 : C 1684 : C******* N.B.: On suppose qu'on a déjà controlle RO, P > 0 1685 : C GAMMA \in (1,3) 1686 : C Si non il faut le faire, en utilisant LOGBOR, 1687 : C LOGNEG, VALER, VAL1, VAL2 1688 : C 1689 : C 1690 : C 1691 : C******* NGCEG = NGCED -> Mur 1692 : C 1693 : IF( NGCEG .EQ. NGCED)THEN 1694 : C 1695 : C********** Sur le mur on fait de reconstruction sur l'etat gauche 1696 : C 1697 : C 1698 : C********** Etat gauche 1699 : 1700 : VALCEL = MPALPP.VPOCHA(NLCEG,1) 1701 : ALCEL = MELLAL.VPOCHA(NLCEG,1) 1702 : DCEL = MPGRAL.VPOCHA(NLCEG,1)*DXG + 1703 : & MPGRAL.VPOCHA(NLCEG,2)*DYG + 1704 : & MPGRAL.VPOCHA(NLCEG,3)*DZG 1705 : alphag = VALCEL + ALCEL*DCEL 1706 : 1707 : VALCEL = MPUVP.VPOCHA(NLCEG,1) 1708 : ALCEL = MELLUV.VPOCHA(NLCEG,1) 1709 : DCEL = MPGRUV.VPOCHA(NLCEG,1)*DXG + 1710 : & MPGRUV.VPOCHA(NLCEG,2)*DYG + 1711 : & MPGRUV.VPOCHA(NLCEG,3)*DZG 1712 : uvxg = VALCEL + ALCEL*DCEL 1713 : 1714 : VALCEL = MPUVP.VPOCHA(NLCEG,2) 1715 : ALCEL = MELLUV.VPOCHA(NLCEG,2) 1716 : DCEL = MPGRUV.VPOCHA(NLCEG,4)*DXG + 1717 : & MPGRUV.VPOCHA(NLCEG,5)*DYG + 1718 : & MPGRUV.VPOCHA(NLCEG,6)*DZG 1719 : uvyg = VALCEL + ALCEL*DCEL 1720 : 1721 : VALCEL = MPUVP.VPOCHA(NLCEG,3) 1722 : ALCEL = MELLUV.VPOCHA(NLCEG,3) 1723 : DCEL = MPGRUV.VPOCHA(NLCEG,7)*DXG + 1724 : & MPGRUV.VPOCHA(NLCEG,8)*DYG + 1725 : & MPGRUV.VPOCHA(NLCEG,9)*DZG 1726 : uvzg = VALCEL + ALCEL*DCEL 1727 : 1728 : uvng = uvxg*CNX + uvyg*CNY + uvzg*CNZ 1729 : uvtg = uvxg*CTX + uvyg*CTY + uvzg*CTZ 1730 : uvvg = uvxg*CVX + uvyg*CVY + uvzg*CVZ 1731 : 1732 : VALCEL = MPULP.VPOCHA(NLCEG,1) 1733 : ALCEL = MELLUL.VPOCHA(NLCEG,1) 1734 : DCEL = MPGRUL.VPOCHA(NLCEG,1)*DXG + 1735 : & MPGRUL.VPOCHA(NLCEG,2)*DYG + 1736 : & MPGRUL.VPOCHA(NLCEG,3)*DZG 1737 : ulxg = VALCEL + ALCEL*DCEL 1738 : 1739 : VALCEL = MPULP.VPOCHA(NLCEG,2) 1740 : ALCEL = MELLUL.VPOCHA(NLCEG,2) 1741 : DCEL = MPGRUL.VPOCHA(NLCEG,4)*DXG + 1742 : & MPGRUL.VPOCHA(NLCEG,5)*DYG + 1743 : & MPGRUL.VPOCHA(NLCEG,6)*DZG 1744 : ulyg = VALCEL + ALCEL*DCEL 1745 : 1746 : VALCEL = MPULP.VPOCHA(NLCEG,3) 1747 : ALCEL = MELLUL.VPOCHA(NLCEG,3) 1748 : DCEL = MPGRUL.VPOCHA(NLCEG,7)*DXG + 1749 : & MPGRUL.VPOCHA(NLCEG,8)*DYG + 1750 : & MPGRUL.VPOCHA(NLCEG,9)*DZG 1751 : ulzg = VALCEL + ALCEL*DCEL 1752 : 1753 : ulng = ulxg*CNX + ulyg*CNY + ulzg*CNZ 1754 : ultg = ulxg*CTX + ulyg*CTY + ulzg*CTZ 1755 : ulvg = ulxg*CVX + ulyg*CVY + ulzg*CVZ 1756 : 1757 : VALCEL = MPPP.VPOCHA(NLCEG, 1) 1758 : ALCEL = MELLP.VPOCHA(NLCEG, 1) 1759 : DCEL = MPGRP.VPOCHA(NLCEG, 1)*DXG + 1760 : & MPGRP.VPOCHA(NLCEG, 2)*DYG + 1761 : & MPGRP.VPOCHA(NLCEG, 3)*DZG 1762 : pg = VALCEL + ALCEL*DCEL 1763 : 1764 : VALCEL = MPTVP.VPOCHA(NLCEG, 1) 1765 : ALCEL = MELLTV.VPOCHA(NLCEG, 1) 1766 : DCEL = MPGRTV.VPOCHA(NLCEG, 1)*DXG + 1767 : & MPGRTV.VPOCHA(NLCEG, 2)*DYG + 1768 : & MPGRTV.VPOCHA(NLCEG, 3)*DZG 1769 : Tvg = VALCEL + ALCEL*DCEL 1770 : 1771 : VALCEL = MPTLP.VPOCHA(NLCEG, 1) 1772 : ALCEL = MELLTL.VPOCHA(NLCEG, 1) 1773 : DCEL = MPGRTL.VPOCHA(NLCEG, 1)*DXG + 1774 : & MPGRTL.VPOCHA(NLCEG, 2)*DYG + 1775 : & MPGRTL.VPOCHA(NLCEG, 3)*DZG 1776 : Tlg = VALCEL + ALCEL*DCEL 1777 : 1778 : C********** We calculate densities again, we use for 1779 : C the moment, perfect and stiffened gas eos 1780 : rvg = gamv*pg/(cpv*Tvg*(gamv - 1)) 1781 : rlg = gaml*(pg + pil)/(cpl*Tlg*(gaml - 1)) 1782 : C 1783 : C********** Si l'on fait pas de prediction, ce n'est pas necessaire de 1784 : C controller la positivite' de la pression et de la densité; elle 1785 : C est déjà garantie par la proprieté LED de limiteur; mais, vu 1786 : C que le limiteur n'est pas calculé ici, mais dans un autre 1787 : C operateur, on le fait 1788 : C 1789 : LOGI1 = (pg .LT. 0.D0) .OR. (rvg .LT. 0.D0) 1790 : & .OR. (rlg .LT. 0.D0) 1791 : C 1792 : IF(LOGI1)THEN 1793 : C 1794 : C************* Premier ordre en espace local 1795 : C 1796 : alphag = MPALPC.VPOCHA(NLCEG,1) 1797 : uvng = MPUVC.VPOCHA(NLCEG,1)*CNX + 1798 : & MPUVC.VPOCHA(NLCEG,2)*CNY + 1799 : & MPUVC.VPOCHA(NLCEG,3)*CNZ 1800 : uvtg = MPUVC.VPOCHA(NLCEG,1)*CTX + 1801 : & MPUVC.VPOCHA(NLCEG,2)*CTY + 1802 : & MPUVC.VPOCHA(NLCEG,3)*CTZ 1803 : uvvg = MPUVC.VPOCHA(NLCEG,1)*CVX + 1804 : & MPUVC.VPOCHA(NLCEG,2)*CVY + 1805 : & MPUVC.VPOCHA(NLCEG,3)*CVZ 1806 : ulng = MPULC.VPOCHA(NLCEG,1)*CNX + 1807 : & MPULC.VPOCHA(NLCEG,2)*CNY + 1808 : & MPULC.VPOCHA(NLCEG,3)*CNZ 1809 : ultg = MPULC.VPOCHA(NLCEG,1)*CTX + 1810 : & MPULC.VPOCHA(NLCEG,2)*CTY + 1811 : & MPULC.VPOCHA(NLCEG,3)*CTZ 1812 : ulvg = MPULC.VPOCHA(NLCEG,1)*CVX + 1813 : & MPULC.VPOCHA(NLCEG,2)*CVY + 1814 : & MPULC.VPOCHA(NLCEG,3)*CVZ 1815 : pg = MPPC.VPOCHA(NLCEG,1) 1816 : Tvg = MPTVC.VPOCHA(NLCEG,1) 1817 : Tlg = MPTLC.VPOCHA(NLCEG,1) 1818 : rvg = MPRVC.VPOCHA(NLCEG,1) 1819 : rlg = MPRLC.VPOCHA(NLCEG,1) 1820 : 1821 : ENDIF 1822 : C 1823 : C********** Son etat droite 1824 : C 1825 : alphad = alphag 1826 : uvnd = -1.d0*uvng 1827 : uvtd = uvtg 1828 : uvvd = uvvg 1829 : ulnd = -1.d0*ulng 1830 : ultd = ultg 1831 : ulvd = ulvg 1832 : pd = pg 1833 : Tvd = Tvg 1834 : Tld = Tlg 1835 : rvd = rvg 1836 : rld = rlg 1837 : C 1838 : C************* Fin cas mur 1839 : C 1840 : ELSE 1841 : C 1842 : C************* Etat gauche 1843 : C 1844 : VALCEL = MPALPP.VPOCHA(NLCEG,1) 1845 : ALCEL = MELLAL.VPOCHA(NLCEG,1) 1846 : DCEL = MPGRAL.VPOCHA(NLCEG,1)*DXG + 1847 : & MPGRAL.VPOCHA(NLCEG,2)*DYG + 1848 : & MPGRAL.VPOCHA(NLCEG,3)*DZG 1849 : alphag = VALCEL + ALCEL*DCEL 1850 : 1851 : VALCEL = MPUVP.VPOCHA(NLCEG,1) 1852 : ALCEL = MELLUV.VPOCHA(NLCEG,1) 1853 : DCEL = MPGRUV.VPOCHA(NLCEG,1)*DXG + 1854 : & MPGRUV.VPOCHA(NLCEG,2)*DYG + 1855 : & MPGRUV.VPOCHA(NLCEG,3)*DZG 1856 : uvxg = VALCEL + ALCEL*DCEL 1857 : 1858 : VALCEL = MPUVP.VPOCHA(NLCEG,2) 1859 : ALCEL = MELLUV.VPOCHA(NLCEG,2) 1860 : DCEL = MPGRUV.VPOCHA(NLCEG,4)*DXG + 1861 : & MPGRUV.VPOCHA(NLCEG,5)*DYG + 1862 : & MPGRUV.VPOCHA(NLCEG,6)*DZG 1863 : uvyg = VALCEL + ALCEL*DCEL 1864 : 1865 : VALCEL = MPUVP.VPOCHA(NLCEG,3) 1866 : ALCEL = MELLUV.VPOCHA(NLCEG,3) 1867 : DCEL = MPGRUV.VPOCHA(NLCEG,7)*DXG + 1868 : & MPGRUV.VPOCHA(NLCEG,8)*DYG + 1869 : & MPGRUV.VPOCHA(NLCEG,9)*DZG 1870 : uvzg = VALCEL + ALCEL*DCEL 1871 : 1872 : uvng = uvxg*CNX + uvyg*CNY + uvzg*CNZ 1873 : uvtg = uvxg*CTX + uvyg*CTY + uvzg*CTZ 1874 : uvvg = uvxg*CVX + uvyg*CVY + uvzg*CVZ 1875 : 1876 : VALCEL = MPULP.VPOCHA(NLCEG,1) 1877 : ALCEL = MELLUL.VPOCHA(NLCEG,1) 1878 : DCEL = MPGRUL.VPOCHA(NLCEG,1)*DXG + 1879 : & MPGRUL.VPOCHA(NLCEG,2)*DYG + 1880 : & MPGRUL.VPOCHA(NLCEG,3)*DZG 1881 : ulxg = VALCEL + ALCEL*DCEL 1882 : 1883 : VALCEL = MPULP.VPOCHA(NLCEG,2) 1884 : ALCEL = MELLUL.VPOCHA(NLCEG,2) 1885 : DCEL = MPGRUL.VPOCHA(NLCEG,4)*DXG + 1886 : & MPGRUL.VPOCHA(NLCEG,5)*DYG + 1887 : & MPGRUL.VPOCHA(NLCEG,6)*DZG 1888 : ulyg = VALCEL + ALCEL*DCEL 1889 : 1890 : VALCEL = MPULP.VPOCHA(NLCEG,3) 1891 : ALCEL = MELLUL.VPOCHA(NLCEG,3) 1892 : DCEL = MPGRUL.VPOCHA(NLCEG,7)*DXG + 1893 : & MPGRUL.VPOCHA(NLCEG,8)*DYG + 1894 : & MPGRUL.VPOCHA(NLCEG,9)*DZG 1895 : ulzg = VALCEL + ALCEL*DCEL 1896 : 1897 : ulng = ulxg*CNX + ulyg*CNY + ulzg*CNZ 1898 : ultg = ulxg*CTX + ulyg*CTY + ulzg*CTZ 1899 : ulvg = ulxg*CVX + ulyg*CVY + ulzg*CVZ 1900 : 1901 : VALCEL = MPPP.VPOCHA(NLCEG, 1) 1902 : ALCEL = MELLP.VPOCHA(NLCEG, 1) 1903 : DCEL = MPGRP.VPOCHA(NLCEG, 1)*DXG + 1904 : & MPGRP.VPOCHA(NLCEG, 2)*DYG + 1905 : & MPGRP.VPOCHA(NLCEG, 3)*DZG 1906 : pg = VALCEL + ALCEL*DCEL 1907 : 1908 : VALCEL = MPTVP.VPOCHA(NLCEG, 1) 1909 : ALCEL = MELLTV.VPOCHA(NLCEG, 1) 1910 : DCEL = MPGRTV.VPOCHA(NLCEG, 1)*DXG + 1911 : & MPGRTV.VPOCHA(NLCEG, 2)*DYG + 1912 : & MPGRTV.VPOCHA(NLCEG, 3)*DZG 1913 : Tvg = VALCEL + ALCEL*DCEL 1914 : 1915 : VALCEL = MPTLP.VPOCHA(NLCEG, 1) 1916 : ALCEL = MELLTL.VPOCHA(NLCEG, 1) 1917 : DCEL = MPGRTL.VPOCHA(NLCEG, 1)*DXG + 1918 : & MPGRTL.VPOCHA(NLCEG, 2)*DYG + 1919 : & MPGRTL.VPOCHA(NLCEG, 3)*DZG 1920 : Tlg = VALCEL + ALCEL*DCEL 1921 : 1922 : C********** We calculate densities again, we use for 1923 : C the moment, perfect and stiffened gas eos 1924 : rvg = gamv*pg/(cpv*Tvg*(gamv - 1)) 1925 : rlg = gaml*(pg + pil)/(cpl*Tlg*(gaml - 1)) 1926 : C 1927 : C********** Positivite 1928 : C 1929 : LOGI1 = (pg .LT. 0.D0) .OR. (rvg .LT. 0.D0) 1930 : & .OR. (rlg .LT. 0.d0) 1931 : C 1932 : IF(LOGI1)THEN 1933 : C 1934 : C************* Premier ordre en espace local 1935 : C 1936 : alphag = MPALPC.VPOCHA(NLCEG,1) 1937 : uvng = MPUVC.VPOCHA(NLCEG,1)*CNX + 1938 : & MPUVC.VPOCHA(NLCEG,2)*CNY + 1939 : & MPUVC.VPOCHA(NLCEG,3)*CNZ 1940 : uvtg = MPUVC.VPOCHA(NLCEG,1)*CTX + 1941 : & MPUVC.VPOCHA(NLCEG,2)*CTY + 1942 : & MPUVC.VPOCHA(NLCEG,3)*CTZ 1943 : uvvg = MPUVC.VPOCHA(NLCEG,1)*CVX + 1944 : & MPUVC.VPOCHA(NLCEG,2)*CVY + 1945 : & MPUVC.VPOCHA(NLCEG,3)*CVZ 1946 : ulng = MPULC.VPOCHA(NLCEG,1)*CNX + 1947 : & MPULC.VPOCHA(NLCEG,2)*CNY + 1948 : & MPULC.VPOCHA(NLCEG,3)*CNZ 1949 : ultg = MPULC.VPOCHA(NLCEG,1)*CTX + 1950 : & MPULC.VPOCHA(NLCEG,2)*CTY + 1951 : & MPULC.VPOCHA(NLCEG,3)*CTZ 1952 : ulvg = MPULC.VPOCHA(NLCEG,1)*CVX + 1953 : & MPULC.VPOCHA(NLCEG,2)*CVY + 1954 : & MPULC.VPOCHA(NLCEG,3)*CVZ 1955 : pg = MPPC.VPOCHA(NLCEG,1) 1956 : Tvg = MPTVC.VPOCHA(NLCEG,1) 1957 : Tlg = MPTLC.VPOCHA(NLCEG,1) 1958 : rvg = MPRVC.VPOCHA(NLCEG,1) 1959 : rlg = MPRLC.VPOCHA(NLCEG,1) 1960 : 1961 : ENDIF 1962 : C 1963 : C********** Etat droite 1964 : C 1965 : VALCEL = MPALPP.VPOCHA(NLCED, 1) 1966 : ALCEL = MELLAL.VPOCHA(NLCED, 1) 1967 : DCEL = MPGRAL.VPOCHA(NLCED, 1)*DXD + 1968 : & MPGRAL.VPOCHA(NLCED, 2)*DYD + 1969 : & MPGRAL.VPOCHA(NLCED, 3)*DZD 1970 : alphad = VALCEL + ALCEL * DCEL 1971 : 1972 : VALCEL = MPUVP.VPOCHA(NLCED, 1) 1973 : ALCEL = MELLUV.VPOCHA(NLCED, 1) 1974 : DCEL = MPGRUV.VPOCHA(NLCED, 1)*DXD + 1975 : & MPGRUV.VPOCHA(NLCED, 2)*DYD + 1976 : & MPGRUV.VPOCHA(NLCED, 3)*DZD 1977 : uvxd = VALCEL + ALCEL * DCEL 1978 : 1979 : VALCEL = MPUVP.VPOCHA(NLCED, 2) 1980 : ALCEL = MELLUV.VPOCHA(NLCED, 2) 1981 : DCEL = MPGRUV.VPOCHA(NLCED, 4)*DXD + 1982 : & MPGRUV.VPOCHA(NLCED, 5)*DYD + 1983 : & MPGRUV.VPOCHA(NLCED, 6)*DZD 1984 : uvyd = VALCEL + ALCEL * DCEL 1985 : 1986 : VALCEL = MPUVP.VPOCHA(NLCED, 3) 1987 : ALCEL = MELLUV.VPOCHA(NLCED, 3) 1988 : DCEL = MPGRUV.VPOCHA(NLCED, 7)*DXD + 1989 : & MPGRUV.VPOCHA(NLCED, 8)*DYD + 1990 : & MPGRUV.VPOCHA(NLCED, 9)*DZD 1991 : uvzd = VALCEL + ALCEL * DCEL 1992 : 1993 : uvnd = uvxd*CNX + uvyd*CNY + uvzd*CNZ 1994 : uvtd = uvxd*CTX + uvyd*CTY + uvzd*CTZ 1995 : uvvd = uvxd*CVX + uvyd*CVY + uvzd*CVZ 1996 : 1997 : VALCEL = MPULP.VPOCHA(NLCED, 1) 1998 : ALCEL = MELLUL.VPOCHA(NLCED, 1) 1999 : DCEL = MPGRUL.VPOCHA(NLCED, 1)*DXD + 2000 : & MPGRUL.VPOCHA(NLCED, 2)*DYD + 2001 : & MPGRUL.VPOCHA(NLCED, 3)*DZD 2002 : ulxd = VALCEL + ALCEL * DCEL 2003 : 2004 : VALCEL = MPULP.VPOCHA(NLCED,2) 2005 : ALCEL = MELLUL.VPOCHA(NLCED, 2) 2006 : DCEL = MPGRUL.VPOCHA(NLCED, 4)*DXD + 2007 : & MPGRUL.VPOCHA(NLCED, 5)*DYD + 2008 : & MPGRUL.VPOCHA(NLCED, 6)*DZD 2009 : ulyd = VALCEL + ALCEL * DCEL 2010 : 2011 : VALCEL = MPULP.VPOCHA(NLCED,3) 2012 : ALCEL = MELLUL.VPOCHA(NLCED, 3) 2013 : DCEL = MPGRUL.VPOCHA(NLCED, 7)*DXD + 2014 : & MPGRUL.VPOCHA(NLCED, 8)*DYD + 2015 : & MPGRUL.VPOCHA(NLCED, 9)*DZD 2016 : ulzd = VALCEL + ALCEL * DCEL 2017 : 2018 : ulnd = ulxd*CNX + ulyd*CNY + ulzd*CNZ 2019 : ultd = ulxd*CTX + ulyd*CTY + ulzd*CTZ 2020 : ulvd = ulxd*CVX + ulyd*CVY + ulzd*CVZ 2021 : 2022 : VALCEL = MPPP.VPOCHA(NLCED, 1) 2023 : ALCEL = MELLP.VPOCHA(NLCED, 1) 2024 : DCEL = MPGRP.VPOCHA(NLCED, 1)*DXD + 2025 : & MPGRP.VPOCHA(NLCED, 2)*DYD + 2026 : & MPGRP.VPOCHA(NLCED, 3)*DZD 2027 : pd = VALCEL + ALCEL * DCEL 2028 : 2029 : VALCEL = MPTVP.VPOCHA(NLCED, 1) 2030 : ALCEL = MELLTV.VPOCHA(NLCED, 1) 2031 : DCEL = MPGRTV.VPOCHA(NLCED, 1)*DXD + 2032 : & MPGRTV.VPOCHA(NLCED, 2)*DYD + 2033 : & MPGRTV.VPOCHA(NLCED, 3)*DZD 2034 : Tvd = VALCEL + ALCEL * DCEL 2035 : 2036 : VALCEL = MPTLP.VPOCHA(NLCED, 1) 2037 : ALCEL = MELLTL.VPOCHA(NLCED, 1) 2038 : DCEL = MPGRTL.VPOCHA(NLCED,1)*DXD + 2039 : & MPGRTL.VPOCHA(NLCED,2)*DYD + 2040 : & MPGRTL.VPOCHA(NLCED,3)*DZD 2041 : Tld = VALCEL + ALCEL * DCEL 2042 : C 2043 : C********** We calculate densities again, we use for 2044 : C the moment, perfect and stiffened gas eos 2045 : rvd = gamv*pd/(cpv*Tvd*(gamv - 1)) 2046 : rld = gaml*(pd + pil)/(cpl*Tld*(gaml - 1)) 2047 : C 2048 : C********** Positivite 2049 : C 2050 : LOGI1 = (pd .LT. 0.0D0) .OR. (rvd .LT. 0.0D0) 2051 : & .OR. (rld .LT. 0.0D0) 2052 : C 2053 : IF(LOGI1)THEN 2054 : C 2055 : C************* Premier ordre en espace local 2056 : C 2057 : alphad = MPALPC.VPOCHA(NLCED,1) 2058 : uvnd = MPUVC.VPOCHA(NLCED,1)*CNX + 2059 : & MPUVC.VPOCHA(NLCED,2)*CNY + 2060 : & MPUVC.VPOCHA(NLCED,3)*CNZ 2061 : uvtd = MPUVC.VPOCHA(NLCED,1)*CTX + 2062 : & MPUVC.VPOCHA(NLCED,2)*CTY + 2063 : & MPUVC.VPOCHA(NLCED,3)*CTZ 2064 : uvvd = MPUVC.VPOCHA(NLCED,1)*CVX + 2065 : & MPUVC.VPOCHA(NLCED,2)*CVY + 2066 : & MPUVC.VPOCHA(NLCED,3)*CVZ 2067 : pd = MPPC.VPOCHA(NLCED,1) 2068 : Tvd = MPTVC.VPOCHA(NLCED,1) 2069 : Tld = MPTLC.VPOCHA(NLCED,1) 2070 : rvd = MPRVC.VPOCHA(NLCED,1) 2071 : rld = MPRLC.VPOCHA(NLCED,1) 2072 : 2073 : ENDIF 2074 : ENDIF 2075 : C 2076 : C******** Les MELVALs 2077 : C 2078 : MLALP.VELCHE(1,NLCF) = alphag 2079 : MLALP.VELCHE(3,NLCF) = alphad 2080 : MLP.VELCHE(1,NLCF) = pg 2081 : MLP.VELCHE(3,NLCF) = pd 2082 : MLTV.VELCHE(1,NLCF) = Tvg 2083 : MLTV.VELCHE(3,NLCF) = Tvd 2084 : MLTL.VELCHE(1,NLCF) = Tlg 2085 : MLTL.VELCHE(3,NLCF) = Tld 2086 : MLRV.VELCHE(1,NLCF) = rvg 2087 : MLRV.VELCHE(3,NLCF) = rvd 2088 : MLRL.VELCHE(1,NLCF) = rlg 2089 : MLRL.VELCHE(3,NLCF) = rld 2090 : 2091 : MLUVN.VELCHE(1,NLCF) = uvng 2092 : MLUVN.VELCHE(3,NLCF) = uvnd 2093 : MLUVT.VELCHE(1,NLCF) = uvtg 2094 : MLUVT.VELCHE(3,NLCF) = uvtd 2095 : MLUVV.VELCHE(1,NLCF) = uvvg 2096 : MLUVV.VELCHE(3,NLCF) = uvvd 2097 : 2098 : MLVNX.VELCHE(1,NLCF) = CNX 2099 : MLVNY.VELCHE(1,NLCF) = CNY 2100 : MLVNZ.VELCHE(1,NLCF) = CNZ 2101 : MLVT1X.VELCHE(1,NLCF) = CTX 2102 : MLVT1Y.VELCHE(1,NLCF) = CTY 2103 : MLVT1Z.VELCHE(1,NLCF) = CTZ 2104 : MLVT2X.VELCHE(1,NLCF) = CVX 2105 : MLVT2Y.VELCHE(1,NLCF) = CVY 2106 : MLVT2Z.VELCHE(1,NLCF) = CVZ 2107 : 2108 : MLULN.VELCHE(1,NLCF) = ulng 2109 : MLULN.VELCHE(3,NLCF) = ulnd 2110 : MLULT.VELCHE(1,NLCF) = ultg 2111 : MLULT.VELCHE(3,NLCF) = ultd 2112 : MLULV.VELCHE(1,NLCF) = ulvg 2113 : MLULV.VELCHE(3,NLCF) = ulvd 2114 : 2115 : MLLNX.VELCHE(1,NLCF) = CNX 2116 : MLLNY.VELCHE(1,NLCF) = CNY 2117 : MLLNZ.VELCHE(1,NLCF) = CNZ 2118 : MLLT1X.VELCHE(1,NLCF) = CTX 2119 : MLLT1Y.VELCHE(1,NLCF) = CTY 2120 : MLLT1Z.VELCHE(1,NLCF) = CTZ 2121 : MLVT2X.VELCHE(1,NLCF) = CVX 2122 : MLVT2Y.VELCHE(1,NLCF) = CVY 2123 : MLVT2Z.VELCHE(1,NLCF) = CVZ 2124 : 2125 : ENDDO 2126 : C 2127 : C**** Desactivation des SEGMENTs 2128 : C 2129 : SEGDES IPT1 2130 : SEGDES IPT2 2131 : C 2132 : C**** Le MPOVALs 'Prediction' sont detruits (si existentes) 2133 : C 2134 : IF(LOGTEM)THEN 2135 : SEGSUP MPALPP 2136 : SEGSUP MPUVP 2137 : SEGSUP MPULP 2138 : SEGSUP MPPP 2139 : SEGSUP MPTVP 2140 : SEGSUP MPTLP 2141 : SEGSUP MPRVP 2142 : SEGSUP MPRLP 2143 : ENDIF 2144 : C 2145 : SEGDES MPALPC 2146 : SEGDES MPGRAL 2147 : SEGDES MPUVC 2148 : SEGDES MPGRUV 2149 : SEGDES MPULC 2150 : SEGDES MPGRUL 2151 : SEGDES MPPC 2152 : SEGDES MPGRP 2153 : SEGDES MPTVC 2154 : SEGDES MPGRTV 2155 : SEGDES MPTLC 2156 : SEGDES MPGRTL 2157 : SEGDES MPRVC 2158 : SEGDES MPRLC 2159 : SEGDES MPNORM 2160 : SEGDES MELLAL 2161 : SEGDES MELLUV 2162 : SEGDES MELLUL 2163 : SEGDES MELLP 2164 : SEGDES MELLTV 2165 : SEGDES MELLTL 2166 : 2167 : SEGDES MLALP 2168 : SEGDES MLP 2169 : SEGDES MLTV 2170 : SEGDES MLTL 2171 : SEGDES MLRV 2172 : SEGDES MLRL 2173 : 2174 : SEGDES MLUVN 2175 : SEGDES MLUVT 2176 : SEGDES MLUVV 2177 : SEGDES MLULN 2178 : SEGDES MLULT 2179 : SEGDES MLULV 2180 : 2181 : SEGDES MLVNX 2182 : SEGDES MLVNY 2183 : SEGDES MLVNZ 2184 : SEGDES MLVT1X 2185 : SEGDES MLVT1Y 2186 : SEGDES MLVT1Z 2187 : SEGDES MLVT2X 2188 : SEGDES MLVT2Y 2189 : SEGDES MLVT2Z 2190 : 2191 : SEGDES MLLNX 2192 : SEGDES MLLNY 2193 : SEGDES MLLNZ 2194 : SEGDES MLLT1X 2195 : SEGDES MLLT1Y 2196 : SEGDES MLLT1Z 2197 : SEGDES MLLT2X 2198 : SEGDES MLLT2Y 2199 : SEGDES MLLT2Z 2200 : C 2201 : C**** Destruction du MELNTI correspondance local/global 2202 : C 2203 : SEGSUP MLENT1 2204 : C 2205 : 9999 CONTINUE 2206 : C 2207 : RETURN 2208 : END 2209 : 2210 : 2211 : 2212 : 2213 : 2214 : 2215 :
© Cast3M 2003 - All rights reserved.
Disclaimer