Télécharger pre52f.eso

Retour à la liste

Numérotation des lignes :

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

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