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

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