Télécharger pre62f.eso

Retour à la liste

Numérotation des lignes :

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

© Cast3M 2003 - Tous droits réservés.
Mentions légales