Télécharger pre221.eso

Retour à la liste

Numérotation des lignes :

  1. C PRE221 SOURCE BECC 10/10/07 21:15:00 6774
  2. SUBROUTINE PRE221(LOGTEM,
  3. & ICEN,IFACE,IFACEL,INORM,
  4. & IROC, IGRROC, IALROC,
  5. & IVITC, IGRVC, IALVC,
  6. & IPC ,IGRPC, IALPC,
  7. & MLRECV, MLRECP, MLMOT1,
  8. & IYC ,IGRYC, IALYC,
  9. & IGAMC,
  10. & DELTAT,
  11. & IROF,IVITF,IPF,IGAMF,IYF,
  12. & LOGAN,LOGNEG,LOGBOR,MESERR,VALER,VAL1,VAL2)
  13. C************************************************************************
  14. C
  15. C PROJET : CASTEM 2000
  16. C
  17. C NOM : PRE221
  18. C
  19. C DESCRIPTION : Voir PRE22
  20. C
  21. C LANGAGE : FORTRAN 77 + ESOPE 2000 (avec extensions CISI)
  22. C
  23. C AUTEUR : A. BECCANTINI, DRN/DMT/SEMT/TTMF
  24. C
  25. C************************************************************************
  26. C
  27. C
  28. C APPELES (Outils) : KRIPAD, LICHT
  29. C
  30. C APPELES (Calcul) : AUCUN
  31. C
  32. C
  33. C************************************************************************
  34. C
  35. C ENTREES
  36. C
  37. C LOGTEM : LOGICAL; si .TRUE. 2em ordre en temps
  38. C sinon 1er ordre en temps
  39. C
  40. C 1) Pointeurs de MELEMEs et de CHPOINTs de la table DOMAINE
  41. C
  42. C ICEN : MELEME de 'POI1' SPG des CENTRES
  43. C
  44. C IFACE : MELEME de 'POI1' SPG des FACES
  45. C
  46. C IFACEL : MELEME de 'SEG3' avec
  47. C CENTRE d'Elt "gauche"
  48. C CENTRE de Face
  49. C CENTRE d'Elt "droite"
  50. C
  51. C N.B. = IFACE.NUM(i,1) = IFACEL.NUM(i,2)
  52. C
  53. C INORM : CHPOINT des cosinus directeurs de normales aux faces
  54. C
  55. C 2) Autres pointeurs
  56. C
  57. C IROC : CHPOINT "CENTRE" contenant la masse volumique RHO
  58. C
  59. C IGRROC : CHPOINT "CENTRE" contenant le gradient de la
  60. C masse volumique RHO (2 composantes)
  61. C
  62. C IALROC : CHPOINT "CENTRE" contenant le limiteur du gradient
  63. C de la masse volumique
  64. C
  65. C IVITC : CHPOINT "CENTRE" contenant la vitesse UX, UY ;
  66. C
  67. C IGRVC : CHPOINT "CENTRE" contenant le gradient de la
  68. C vitesse (4 composantes)
  69. C
  70. C IALVC : CHPOINT "CENTRE" contenant le limiteur du gradient
  71. C de la vitesse (2 composantes)
  72. C
  73. C IPC : CHPOINT "CENTRE" contenat la pression P;
  74. C
  75. C IGRPC : CHPOINT "CENTRE" contenant le gradient de la
  76. C pression (2 composantes)
  77. C
  78. C IALPC : CHPOINT "CENTRE" contenant le limiteur du gradient
  79. C de la pression
  80. C
  81. C MLRECV : pointeur de la LISTREEL des cv
  82. C
  83. C MLRECP : pointeur de la LISTREEL des cp
  84. C
  85. C MLMOT1 : pointeur des noms des gaz qui apparessent dans les
  86. C equations d'Euler
  87. C
  88. C IYC : CHPOINT "CENTRE" contenat les fractions massiques ;
  89. C
  90. C IGRYC : CHPOINT "CENTRE" contenant les gradient des fr.mass.;
  91. C
  92. C IALPC : CHPOINT "CENTRE" contenant les limiteurs des gradients
  93. C des fr.mass.;
  94. C
  95. C IGAMC : CHPOINT "CENTRE" contenat le "Gamma" du gaz
  96. C
  97. C 3)
  98. C
  99. C DELTAT : REAL*8, encrement en temps pour calculer la prediction
  100. C
  101. C
  102. C SORTIES
  103. C
  104. C
  105. C IROF : MCHAML defini sur le MELEME de pointeur IFACEL,
  106. C contenant la masse volumique RHO
  107. C
  108. C IVITF : MCHAML defini sur le MELEME de pointeur IFACEL,
  109. C contenant la vitesse UN, UT dans le repaire local
  110. C (n,t) et defini sur le MELEME de pointeur IFACE,
  111. C contenant les cosinus directeurs du repere local
  112. C
  113. C IPF : MCHAML defini sur le MELEME de pointeur IFACEL,
  114. C contenant la pression P
  115. C
  116. C IYF : MCHAML defini sur le MELEME de pointeur IFACEL,
  117. C contenant les fractions massiques;
  118. C
  119. C IGAMF : MCHAML defini sur le MELEME de pointeur IFACEL,
  120. C contenant le "gamma" du gaz
  121. C
  122. C LOGAN : anomalie detectee
  123. C
  124. C LOGNEG : (LOGICAL): si .TRUE. une pression ou une densité
  125. C negative a été detectée -> en interactif le
  126. C programme s'arrete en GIBIANE
  127. C (erreur stocké en MESERR et VALER)
  128. C
  129. C LOGBOR : (LOGICAL): si .TRUE. un gamma a ete detecte
  130. C dehor 1 et 3 (sa valeur stockée en MESERR et VALER;
  131. C en VAL1 et en VAL2 on stocke 1.0 et 3.0)
  132. C
  133. C MESERR
  134. C VALER
  135. C VAL1,
  136. C VAL2 : pour les messages d'erreur
  137. C
  138. C************************************************************************
  139. C
  140. C HISTORIQUE (Anomalies et modifications éventuelles)
  141. C
  142. C HISTORIQUE : Créée le 10.7.98.
  143. C
  144. C************************************************************************
  145. C
  146. C
  147. C ATTENTION: Cet programme marche si le MAILLAGE est convex;
  148. C si non il faut changer l'algoritme de calcul de
  149. C l'orientation des normales aux faces.
  150. C
  151. C La positivité n'est pas controlle parce que c'est déjà fait
  152. C dans l'operateur PRIM
  153. C
  154. C
  155. C************************************************************************
  156. C
  157. C**** Variables de COOPTIO
  158. C
  159. C INTEGER IPLLB, IERPER, IERMAX, IERR, INTERR
  160. C & ,IOTER, IOLEC, IOIMP, IOCAR, IOACQ
  161. C & ,IOPER, IOSGB, IOGRA, IOSAU, IORES
  162. C & ,IECHO, IIMPI, IOSPI
  163. C & ,IDIM
  164. CC & ,MCOORD
  165. C & ,IFOMOD, NIFOUR, IFOUR, NSDPGE, IONIVE
  166. C & ,NGMAXY, IZROSF, ISOTYP, IOSCR,LTEXLU
  167. C & ,NORINC,NORVAL,NORIND,NORVAD
  168. C & ,NUCROU, IPSAUV
  169. C
  170. C**** Les variables
  171. C
  172. IMPLICIT INTEGER(I-N)
  173. INTEGER ICEN, IFACE, IFACEL, INORM,IROC, IGRROC, IALROC
  174. & , IVITC, IGRVC, IALVC
  175. & , IPC ,IGRPC, IALPC
  176. & , IYC, IGRYC, IALYC
  177. & , IGAMC
  178. & , IROF, IVITF, IPF, IGAMF, IYF
  179. & , IGEOM, NFAC, NCEN
  180. & , N1PTEL, N1EL, N2PTEL, N2EL, N2, N1, N3, L1, NLCE
  181. & , NLCF, NGCF, NGCEG, NLCEG, NGCED, NLCED, NGCF1
  182. & , IDIMP1, INDCEL, I1, NESP
  183. REAL*8 VALER, VAL1, VAL2, XG, YG, XC, YC, XD, YD, DELTAT
  184. & ,DXG, DYG, DXD, DYD
  185. & , CNX, CNY, CTX, CTY, ORIENT
  186. & , ROG, PG, GAMG, UXG, UYG, UNG, UTG
  187. & , ROD, PD, GAMD, UXD, UYD, UND, UTD
  188. & , VALCEL, DCEL, ALCEL
  189. & , DROX, DROY, DUXX, DUXY, DUYX, DUYY, DPX, DPY
  190. & , DRO, DUX, DUY, DP
  191. & , ALPHA, CPCEL, CVCEL, DYMAS, SUMY
  192. CHARACTER*(40) MESERR
  193. CHARACTER*(8) TYPE
  194. CHARACTER*(4) CARCEL
  195. LOGICAL LOGAN,LOGNEG, LOGBOR, LOGTEM, LOGI1, LOGI2
  196. C
  197. C
  198. C**** Les Includes
  199. C
  200. -INC SMCOORD
  201. -INC CCOPTIO
  202. -INC SMCHPOI
  203. POINTEUR MPROC.MPOVAL, MPGRR.MPOVAL,
  204. & MPVITC.MPOVAL, MPGRV.MPOVAL,
  205. & MPPC.MPOVAL, MPGRP.MPOVAL,
  206. & MPYC.MPOVAL, MPGRY.MPOVAL,
  207. & MPGAMC.MPOVAL, MPNORM.MPOVAL,
  208. & MPROP.MPOVAL, MPPP.MPOVAL, MPVITP.MPOVAL,
  209. & MPYP.MPOVAL
  210. -INC SMCHAML
  211. POINTEUR MCHAMY.MCHAML
  212. POINTEUR MELVNX.MELVAL, MELVNY.MELVAL,
  213. & MELT1X.MELVAL, MELT1Y.MELVAL
  214. POINTEUR MELVUN.MELVAL, MELVUT.MELVAL
  215. POINTEUR MELRO.MELVAL, MELP.MELVAL,
  216. & MELGAM.MELVAL
  217. POINTEUR MELALR.MPOVAL,
  218. & MELALV.MPOVAL,
  219. & MELALP.MPOVAL,
  220. & MELALY.MPOVAL
  221. -INC SMLENTI
  222. -INC SMELEME
  223. -INC SMLREEL
  224. POINTEUR MLRECV.MLREEL, MLRECP.MLREEL
  225. -INC SMLMOTS
  226. C
  227. C**** Segments des fractions massiques gauche et droit
  228. C
  229. SEGMENT FRAMAS
  230. REAL*8 FRAMG(NESP), FRAMD(NESP)
  231. ENDSEGMENT
  232. C
  233. SEGACT MLMOT1
  234. C
  235. C
  236. C**** Initialisation des parametres d'erreur déjà faite, i.e.
  237. C
  238. C LOGNEG = .FALSE.
  239. C LOGBOR = .FALSE.
  240. C MESERR = ' '
  241. C MOTERR(1:40) = MESERR(1:40)
  242. C VALER = 0.0D0
  243. C VAL1 = 0.0D0
  244. C VAL2 = 0.0D0
  245. C
  246. C AB
  247. C 2010/09/14: LOGI1 computed can destroy the symmetry!
  248. C LOGI1 = .FALSE.
  249. LOGI1 = .FALSE.
  250. C
  251. C**** KRIPAD pour la correspondance global/local de centre
  252. C
  253. CALL KRIPAD(ICEN,MLENT1)
  254. C
  255. C**** MLENTI1 a MCORD.XCOORD(/1)/(IDIM+1) elements
  256. C
  257. C Si i est le numero global d'un noeud de ICEN,
  258. C MLENT1.LECT(i) contient sa position, i.e.
  259. C
  260. C I = numero global du noeud centre
  261. C MLENT1.LECT(i) = numero local du noeud centre
  262. C
  263. C MLENT1 déjà activé, i.e.
  264. C
  265. C SEGACT MLENT1
  266. C
  267. C**** Activation de CHPOINTs
  268. C
  269. C densité + grad + limiteur
  270. C vitesse + grad + limiteur
  271. C pression + grad + limiteur
  272. C fract.mass + grad + limiteur
  273. C gamma
  274. C cosinus directeurs des normales aux surface
  275. C
  276. CALL LICHT(IROC , MPROC , TYPE, IGEOM)
  277. CALL LICHT(IGRROC, MPGRR , TYPE, IGEOM)
  278. CALL LICHT(IVITC, MPVITC , TYPE, IGEOM)
  279. CALL LICHT(IGRVC, MPGRV , TYPE, IGEOM)
  280. CALL LICHT(IPC , MPPC , TYPE, IGEOM)
  281. CALL LICHT(IGRPC, MPGRP , TYPE, IGEOM)
  282. CALL LICHT(IYC , MPYC , TYPE, IGEOM)
  283. CALL LICHT(IGRYC, MPGRY , TYPE, IGEOM)
  284. CALL LICHT(IGAMC, MPGAMC , TYPE, IGEOM)
  285. CALL LICHT(INORM, MPNORM , TYPE, IGEOM)
  286. C
  287. C**** Les MPOVALs 'Prediction'
  288. C
  289. IF(LOGTEM)THEN
  290. SEGINI, MPROP = MPROC
  291. SEGINI, MPPP = MPPC
  292. SEGINI, MPVITP = MPVITC
  293. SEGINI, MPYP = MPYC
  294. ELSE
  295. MPROP = MPROC
  296. MPPP = MPPC
  297. MPVITP = MPVITC
  298. MPYP = MPYC
  299. ENDIF
  300. C
  301. C**** Les Limiteurs
  302. C
  303. CALL LICHT(IALROC, MELALR , TYPE, IGEOM)
  304. CALL LICHT(IALVC, MELALV , TYPE, IGEOM)
  305. CALL LICHT(IALPC, MELALP , TYPE, IGEOM)
  306. CALL LICHT(IALYC, MELALY , TYPE, IGEOM)
  307. C
  308. C
  309. C**** Les MPOVAL sont déjà activés i.e.:
  310. C
  311. C SEGACT MPROC
  312. C SEGACT MPGRR
  313. C SEGACT MPIALR
  314. C SEGACT MPVITC
  315. C SEGACT MPGRV
  316. C SEGACT MPIALV
  317. C SEGACT MPPC
  318. C SEGACT MPGRP
  319. C SEGACT MPIALP
  320. C SEGACT MPYC
  321. C SEGACT MPGRY
  322. C SEGACT MPIALY
  323. C SEGACT MPGAMC
  324. C SEGACT MPNORM
  325. C
  326. C**** Le MELEME FACEL
  327. C
  328. IPT1 = IFACEL
  329. IPT2 = IFACE
  330. SEGACT IPT1
  331. SEGACT IPT2
  332. NFAC = IPT1.NUM(/2)
  333. C
  334. C**** Creation de MCHAMLs contenant les etats gauche et droite,
  335. C
  336. C i.e.:
  337. C
  338. C vitesse + cosinus directors du repere local
  339. C densité
  340. C pression
  341. C gamma
  342. C
  343. C**** Cosinus directors du repere local et vitesse
  344. C
  345. C Les cosinus directeurs
  346. C
  347. N1 = 2
  348. N3 = 6
  349. L1 = 28
  350. SEGINI MCHEL1
  351. IVITF = MCHEL1
  352. MCHEL1.TITCHE = 'U '
  353. MCHEL1.IMACHE(1) = IFACE
  354. MCHEL1.IMACHE(2) = IFACEL
  355. MCHEL1.CONCHE(1) = ' (n,t) in (x,y) '
  356. MCHEL1.CONCHE(2) = ' U in (n,t) '
  357. C
  358. C**** Valeurs des cosinus definies par respect au repair global, i.e.
  359. C
  360. MCHEL1.INFCHE(1,1) = 2
  361. MCHEL1.INFCHE(1,3) = NIFOUR
  362. MCHEL1.INFCHE(1,4) = 0
  363. MCHEL1.INFCHE(1,5) = 0
  364. MCHEL1.INFCHE(1,6) = 0
  365. MCHEL1.IFOCHE = IFOUR
  366. C
  367. C**** Valeurs de vitesse definies par respect au repair local, i.e.
  368. C
  369. MCHEL1.INFCHE(2,1) = 1
  370. MCHEL1.INFCHE(2,3) = NIFOUR
  371. MCHEL1.INFCHE(2,4) = 0
  372. MCHEL1.INFCHE(2,5) = 0
  373. MCHEL1.INFCHE(2,6) = 0
  374. C
  375. C**** Le cosinus directeurs
  376. C
  377. N1PTEL = 1
  378. N1EL = NFAC
  379. N2PTEL = 0
  380. N2EL = 0
  381. C
  382. C**** MCHAML a N2 composantes:
  383. C
  384. C cosinus directeurs du repere local (n,t1)
  385. C
  386. C IDIM = 2 -> 4 composantes
  387. C
  388. N2 = 4
  389. SEGINI MCHAM1
  390. MCHEL1.ICHAML(1) = MCHAM1
  391. MCHAM1.NOMCHE(1) = 'NX '
  392. MCHAM1.NOMCHE(2) = 'NY '
  393. MCHAM1.NOMCHE(3) = 'TX '
  394. MCHAM1.NOMCHE(4) = 'TY '
  395. MCHAM1.TYPCHE(1) = 'REAL*8 '
  396. MCHAM1.TYPCHE(2) = 'REAL*8 '
  397. MCHAM1.TYPCHE(3) = 'REAL*8 '
  398. MCHAM1.TYPCHE(4) = 'REAL*8 '
  399. SEGINI MELVNX
  400. SEGINI MELVNY
  401. SEGINI MELT1X
  402. SEGINI MELT1Y
  403. MCHAM1.IELVAL(1) = MELVNX
  404. MCHAM1.IELVAL(2) = MELVNY
  405. MCHAM1.IELVAL(3) = MELT1X
  406. MCHAM1.IELVAL(4) = MELT1Y
  407. SEGDES MCHAM1
  408. C
  409. C**** Vitesse
  410. C
  411. N1EL = NFAC
  412. N1PTEL = 3
  413. N2EL = 0
  414. N2PTEL = 0
  415. C
  416. C**** MCHAML a N2 composantes:
  417. C
  418. C IDIM = 2 -> 2 composantes
  419. C
  420. N2 = 2
  421. SEGINI MCHAM1
  422. MCHEL1.ICHAML(2) = MCHAM1
  423. SEGDES MCHEL1
  424. MCHAM1.NOMCHE(1) = 'UN '
  425. MCHAM1.NOMCHE(2) = 'UT '
  426. MCHAM1.TYPCHE(1) = 'REAL*8 '
  427. MCHAM1.TYPCHE(2) = 'REAL*8 '
  428. SEGINI MELVUN
  429. SEGINI MELVUT
  430. MCHAM1.IELVAL(1) = MELVUN
  431. MCHAM1.IELVAL(2) = MELVUT
  432. SEGDES MCHAM1
  433. C
  434. C**** Densite
  435. C
  436. N1 = 1
  437. N3 = 6
  438. L1 = 15
  439. SEGINI MCHEL2
  440. IROF = MCHEL2
  441. MCHEL2.IMACHE(1) = IFACEL
  442. MCHEL2.TITCHE = 'RO '
  443. MCHEL2.CONCHE(1) = ' '
  444. C
  445. C**** Valeurs independente du repére, i.e.
  446. C
  447. MCHEL2.INFCHE(1,1) = 0
  448. MCHEL2.INFCHE(1,3) = NIFOUR
  449. MCHEL2.INFCHE(1,4) = 0
  450. MCHEL2.INFCHE(1,5) = 0
  451. MCHEL2.INFCHE(1,6) = 0
  452. MCHEL2.IFOCHE = IFOUR
  453. N2 = 1
  454. SEGINI MCHAM1
  455. MCHEL2.ICHAML(1) = MCHAM1
  456. SEGDES MCHEL2
  457. MCHAM1.NOMCHE(1) = 'SCAL '
  458. MCHAM1.TYPCHE(1) = 'REAL*8 '
  459. SEGINI MELRO
  460. MCHAM1.IELVAL(1) = MELRO
  461. SEGDES MCHAM1
  462. C
  463. C**** Pression
  464. C
  465. MCHEL1 = IROF
  466. SEGINI, MCHEL2 = MCHEL1
  467. IPF = MCHEL2
  468. MCHEL2.TITCHE = 'P '
  469. C
  470. C**** MCHAM1 = MCHAML de la densite
  471. C
  472. SEGINI, MCHAM2 = MCHAM1
  473. MCHEL2.ICHAML(1) = MCHAM2
  474. SEGDES MCHEL2
  475. SEGINI MELP
  476. MCHAM2.IELVAL(1) = MELP
  477. SEGDES MCHAM2
  478. C
  479. C**** Gamma
  480. C
  481. MCHEL1 = IROF
  482. SEGINI, MCHEL2 = MCHEL1
  483. IGAMF = MCHEL2
  484. MCHEL2.TITCHE = 'GAMMA '
  485. C
  486. C**** MCHAM1 = MCHAML de la densite
  487. C
  488. SEGINI, MCHAM2 = MCHAM1
  489. MCHEL2.ICHAML(1) = MCHAM2
  490. SEGDES MCHEL2
  491. SEGINI MELGAM
  492. MCHAM2.IELVAL(1) = MELGAM
  493. SEGDES MCHAM2
  494. C
  495. C**** Le CHAMELEM des fractions massiques
  496. C
  497. NESP = MPYC.VPOCHA(/2)
  498. SEGINI FRAMAS
  499. MCHEL1 = IROF
  500. SEGINI, MCHEL2 = MCHEL1
  501. IYF = MCHEL2
  502. MCHEL2.TITCHE = 'Y '
  503. N2 = NESP
  504. SEGINI MCHAMY
  505. MCHEL2.ICHAML(1) = MCHAMY
  506. SEGDES MCHEL2
  507. N1EL = NFAC
  508. N1PTEL = 3
  509. N2EL = 0
  510. N2PTEL = 0
  511. DO I1 = 1, NESP
  512. SEGINI MELVA1
  513. MCHAMY.IELVAL(I1) = MELVA1
  514. C CARCEL = MLMOT1.MOTS(1)
  515. C AB Error. It is corrected by the following line
  516. CARCEL = MLMOT1.MOTS(I1)
  517. TYPE = ' '
  518. TYPE(1:4) = CARCEL(1:4)
  519. MCHAMY.NOMCHE(I1) = TYPE
  520. MCHAMY.TYPCHE(I1) = 'REAL*8 '
  521. ENDDO
  522. C
  523. C
  524. C***********************************************************************
  525. C********* PREDICTION **************************************************
  526. C***********************************************************************
  527. C
  528. C**** Prediction avec gradients limités
  529. C
  530. C
  531. IF(LOGTEM)THEN
  532. C
  533. IPT3 = ICEN
  534. SEGACT IPT3
  535. NCEN = IPT3.NUM(/2)
  536. DO NLCE = 1, NCEN
  537. ROG = MPROP.VPOCHA(NLCE,1)
  538. UXG = MPVITP.VPOCHA(NLCE,1)
  539. UYG = MPVITP.VPOCHA(NLCE,2)
  540. PG = MPPP.VPOCHA(NLCE,1)
  541. DROX = MPGRR.VPOCHA(NLCE,1)*MELALR.VPOCHA(NLCE,1)
  542. DROY = MPGRR.VPOCHA(NLCE,2)*MELALR.VPOCHA(NLCE,1)
  543. DUXX = MPGRV.VPOCHA(NLCE,1)*MELALV.VPOCHA(NLCE,1)
  544. DUXY = MPGRV.VPOCHA(NLCE,2)*MELALV.VPOCHA(NLCE,1)
  545. DUYX = MPGRV.VPOCHA(NLCE,3)*MELALV.VPOCHA(NLCE,2)
  546. DUYY = MPGRV.VPOCHA(NLCE,4)* MELALV.VPOCHA(NLCE,2)
  547. DPX = MPGRP.VPOCHA(NLCE,1)* MELALP.VPOCHA(NLCE,1)
  548. DPY = MPGRP.VPOCHA(NLCE,2)* MELALP.VPOCHA(NLCE,1)
  549. GAMG = MPGAMC.VPOCHA(NLCE,1)
  550. DRO = UXG * DROX + ROG * ( DUXX + DUYY )
  551. & + UYG * DROY
  552. DUX = UXG * DUXX + DPX / ROG + UYG * DUXY
  553. DUY = UXG * DUYX + UYG * DUYY + DPY / ROG
  554. DP = GAMG * PG * (DUXX + DUYY)
  555. & + UXG * DPX + UYG * DPY
  556. C
  557. MPROP.VPOCHA(NLCE,1) = ROG - DELTAT * DRO
  558. MPVITP.VPOCHA(NLCE,1) = UXG - DELTAT * DUX
  559. MPVITP.VPOCHA(NLCE,2) = UYG - DELTAT * DUY
  560. MPPP.VPOCHA(NLCE,1) = PG - DELTAT * DP
  561. DO I1 = 1, NESP
  562. INDCEL = 2 * I1 - 1
  563. ALPHA = MELALY.VPOCHA(NLCE,I1)
  564. DYMAS = UXG * MPGRY.VPOCHA(NLCE,INDCEL) * ALPHA +
  565. & UYG * MPGRY.VPOCHA(NLCE,INDCEL+1) * ALPHA
  566. MPYP.VPOCHA(NLCE,I1) = MPYC.VPOCHA(NLCE,I1) -
  567. & DELTAT * DYMAS
  568. ENDDO
  569. ENDDO
  570. C
  571. ENDIF
  572. C
  573. C
  574. C***********************************************************************
  575. C********* CORRECTION **************************************************
  576. C***********************************************************************
  577. C
  578. C**** Boucle sur le faces
  579. C
  580. IDIMP1 = IDIM + 1
  581. DO NLCF = 1, NFAC
  582. C
  583. C******* NLCF = numero local du centre de face
  584. C NGCF = numero global du centre de face
  585. C NGCEG = numero global du centre ELT "gauche"
  586. C NLCEG = numero local du centre ELT "gauche"
  587. C NGCED = numero global du centre ELT "droite"
  588. C NLCED = numero local du centre ELT "droite"
  589. C
  590. NGCEG = IPT1.NUM(1,NLCF)
  591. NGCF = IPT1.NUM(2,NLCF)
  592. NGCED = IPT1.NUM(3,NLCF)
  593. NLCEG = MLENT1.LECT(NGCEG)
  594. NLCED = MLENT1.LECT(NGCED)
  595. C
  596. C******* TEST: IPT2.NUM(1,NLCF) = IPT1.NUM(2,NLCF)
  597. C
  598. NGCF1 = IPT2.NUM(1,NLCF)
  599. IF( NGCF1 .NE. NGCF) THEN
  600. LOGAN = .TRUE.
  601. MESERR(1:40) = 'PRET, subroutine pre221.eso '
  602. GOTO 9999
  603. ENDIF
  604. C
  605. C******* Cosinus directeurs des NORMALES aux faces
  606. C
  607. C On impose que les normales sont direct "Gauche" -> "Centre"
  608. C
  609. INDCEL = (NGCEG-1)*IDIMP1
  610. XG = XCOOR(INDCEL+1)
  611. YG = XCOOR(INDCEL+2)
  612. INDCEL = (NGCF-1)*IDIMP1
  613. XC = XCOOR(INDCEL + 1)
  614. YC = XCOOR(INDCEL + 2)
  615. INDCEL = (NGCED-1)*IDIMP1
  616. XD = XCOOR(INDCEL+1)
  617. YD = XCOOR(INDCEL+2)
  618. DXG = XC - XG
  619. DYG = YC - YG
  620. DXD = XC - XD
  621. DYD = YC - YD
  622. C
  623. C******* On calcule le sign du pruduit scalare
  624. C (Normales de Castem) * (vecteur "gauche" -> "centre")
  625. C
  626. CNX = MPNORM.VPOCHA(NLCF,1)
  627. CNY = MPNORM.VPOCHA(NLCF,2)
  628. ORIENT = CNX * DXG + CNY * DYG
  629. ORIENT = SIGN(1.0D0,ORIENT)
  630. IF(ORIENT .NE. 1.0D0)THEN
  631. LOGAN = .TRUE.
  632. MESERR(1:30)=
  633. & 'PRET , subroutine pre221.eso. '
  634. GOTO 9999
  635. ENDIF
  636. CNX = CNX * ORIENT
  637. CNY = CNY * ORIENT
  638. C
  639. C********** Cosinus directeurs de tangent 2D
  640. C
  641. CTX = -1.0D0 * CNY
  642. CTY = CNX
  643. C
  644. C
  645. C******* Les autres MELVALs
  646. C
  647. C
  648. C******* N.B.: On suppose qu'on a déjà controlle RO, P > 0
  649. C GAMMA \in (1,3)
  650. C Si non il faut le faire, en utilisant LOGBOR,
  651. C LOGNEG, VALER, VAL1, VAL2
  652. C
  653. C
  654. C
  655. C******* NGCEG = NGCED -> Mur
  656. C
  657. IF( NGCEG .EQ. NGCED)THEN
  658. C
  659. C********** Sur le mur on fait de reconstruction sur l'etat gauche
  660. C
  661. C
  662. C********** Etat gauche
  663. C
  664. VALCEL = MPROP.VPOCHA(NLCEG, 1)
  665. ALCEL = MELALR.VPOCHA(NLCEG, 1)
  666. DCEL = MPGRR.VPOCHA(NLCEG, 1)*DXG +
  667. & MPGRR.VPOCHA(NLCEG, 2)*DYG
  668. ROG = VALCEL + ALCEL * DCEL
  669. C
  670. VALCEL = MPPP.VPOCHA(NLCEG, 1)
  671. ALCEL = MELALP.VPOCHA(NLCEG, 1)
  672. DCEL = MPGRP.VPOCHA(NLCEG, 1)*DXG +
  673. & MPGRP.VPOCHA(NLCEG, 2)*DYG
  674. PG = VALCEL + ALCEL * DCEL
  675. C
  676. LOGI2 = .FALSE.
  677. SUMY = 0.0D0
  678. CPCEL = 0.0D0
  679. CVCEL = 0.0D0
  680. DO I1 = 1, NESP
  681. INDCEL = 2 * I1 - 1
  682. VALCEL = MPYP.VPOCHA(NLCEG,I1)
  683. ALCEL = MELALY.VPOCHA(NLCEG, I1)
  684. DCEL = MPGRY.VPOCHA(NLCEG, INDCEL)*DXG +
  685. & MPGRY.VPOCHA(NLCEG,INDCEL + 1 )*DYG
  686. ALCEL = VALCEL + ALCEL * DCEL
  687. SUMY = SUMY + ALCEL
  688. CPCEL = CPCEL + MLRECP.PROG(I1) * ALCEL
  689. CVCEL = CVCEL + MLRECV.PROG(I1) * ALCEL
  690. LOGI2 = LOGI2 .OR. (ALCEL .LT. 0.0D0)
  691. FRAMAS.FRAMG(I1) = ALCEL
  692. ENDDO
  693. LOGI2 = LOGI2 .OR. (SUMY .GT. 1.0D0)
  694. CPCEL = CPCEL + (1.0D0 - SUMY) * MLRECP.PROG(NESP+1)
  695. CVCEL = CVCEL + (1.0D0 - SUMY) * MLRECV.PROG(NESP+1)
  696. GAMG = CPCEL / CVCEL
  697. C
  698. VALCEL = MPVITP.VPOCHA(NLCEG, 1)
  699. ALCEL = MELALV.VPOCHA(NLCEG, 1)
  700. DCEL = MPGRV.VPOCHA(NLCEG, 1)*DXG +
  701. & MPGRV.VPOCHA(NLCEG, 2)*DYG
  702. UXG = VALCEL + ALCEL * DCEL
  703. C
  704. VALCEL = MPVITP.VPOCHA(NLCEG, 2)
  705. ALCEL = MELALV.VPOCHA(NLCEG, 2)
  706. DCEL = MPGRV.VPOCHA(NLCEG, 3)*DXG +
  707. & MPGRV.VPOCHA(NLCEG, 4)*DYG
  708. UYG = VALCEL + ALCEL * DCEL
  709. C
  710. UNG = UXG * CNX + UYG * CNY
  711. UTG = UXG * CTX + UYG * CTY
  712. C
  713. C********** Si l'on fait pas de prediction, ce n'est pas necessaire de
  714. C controller la positivite' de la pression et de la densité; elle
  715. C est déjà garantie par la proprieté LED de limiteur; mais, vu
  716. C que le limiteur n'est pas calculé ici, mais dans un autre
  717. C operateur, on le fait
  718. C
  719. C LOGI1 = (PG .LT. 0.0D0) .OR. (ROG .LT. 0.0D0) .OR. LOGI2
  720. C 2010/09/14: AB this destroys the symmetry!
  721. C LOGI1 = .FALSE.
  722. C
  723. C
  724. IF(LOGI1)THEN
  725. C
  726. C************* Premier ordre en espace local
  727. C
  728. ROG = MPROC.VPOCHA(NLCEG,1)
  729. ROD = ROG
  730. PG = MPPC.VPOCHA(NLCEG,1)
  731. PD = PG
  732. UNG = MPVITC.VPOCHA(NLCEG,1)*CNX +
  733. & MPVITC.VPOCHA(NLCEG,2)*CNY
  734. UTG = MPVITC.VPOCHA(NLCEG,1)*CTX +
  735. & MPVITC.VPOCHA(NLCEG,2)*CTY
  736. UND = -1.0D0 * UNG
  737. UTD = UTG
  738. DO I1 = 1, NESP
  739. FRAMAS.FRAMG(I1) = MPYC.VPOCHA(NLCEG,I1)
  740. FRAMAS.FRAMD(I1) = FRAMAS.FRAMG(I1)
  741. ENDDO
  742. GAMG = MPGAMC.VPOCHA(NLCEG,1)
  743. GAMD = GAMG
  744. ELSE
  745. C
  746. C********** Son etat droite
  747. C
  748. ROD = ROG
  749. PD = PG
  750. UND = -1.0D0 * UNG
  751. UTD = UTG
  752. DO I1 = 1, NESP
  753. FRAMAS.FRAMD(I1) = FRAMAS.FRAMG(I1)
  754. ENDDO
  755. GAMD = GAMG
  756. ENDIF
  757. C
  758. C************* Fin cas mur
  759. C
  760. ELSE
  761. C
  762. C************* Etat gauche
  763. C
  764. VALCEL = MPROP.VPOCHA(NLCEG, 1)
  765. ALCEL = MELALR.VPOCHA(NLCEG, 1)
  766. DCEL = MPGRR.VPOCHA(NLCEG, 1)*DXG +
  767. & MPGRR.VPOCHA(NLCEG, 2)*DYG
  768. ROG = VALCEL + ALCEL * DCEL
  769. C
  770. VALCEL = MPPP.VPOCHA(NLCEG, 1)
  771. ALCEL = MELALP.VPOCHA(NLCEG, 1)
  772. DCEL = MPGRP.VPOCHA(NLCEG, 1)*DXG +
  773. & MPGRP.VPOCHA(NLCEG, 2)*DYG
  774. PG = VALCEL + ALCEL * DCEL
  775. C
  776. LOGI2 = .FALSE.
  777. SUMY = 0.0D0
  778. CPCEL = 0.0D0
  779. CVCEL = 0.0D0
  780. DO I1 = 1, NESP
  781. INDCEL = 2 * I1 - 1
  782. VALCEL = MPYP.VPOCHA(NLCEG,I1)
  783. ALCEL = MELALY.VPOCHA(NLCEG, I1)
  784. DCEL = MPGRY.VPOCHA(NLCEG, INDCEL)*DXG +
  785. & MPGRY.VPOCHA(NLCEG,INDCEL + 1 )*DYG
  786. ALCEL = VALCEL + ALCEL * DCEL
  787. SUMY = SUMY + ALCEL
  788. CPCEL = CPCEL + MLRECP.PROG(I1) * ALCEL
  789. CVCEL = CVCEL + MLRECV.PROG(I1) * ALCEL
  790. LOGI2 = LOGI2 .OR. (ALCEL .LT. 0.0D0)
  791. FRAMAS.FRAMG(I1) = ALCEL
  792. ENDDO
  793. LOGI2 = LOGI2 .OR. (SUMY .GT. 1.0D0)
  794. CPCEL = CPCEL + (1.0D0 - SUMY) * MLRECP.PROG(NESP+1)
  795. CVCEL = CVCEL + (1.0D0 - SUMY) * MLRECV.PROG(NESP+1)
  796. GAMG = CPCEL / CVCEL
  797. C
  798. VALCEL = MPVITP.VPOCHA(NLCEG, 1)
  799. ALCEL = MELALV.VPOCHA(NLCEG, 1)
  800. DCEL = MPGRV.VPOCHA(NLCEG, 1)*DXG +
  801. & MPGRV.VPOCHA(NLCEG, 2)*DYG
  802. UXG = VALCEL + ALCEL * DCEL
  803. C
  804. VALCEL = MPVITP.VPOCHA(NLCEG, 2)
  805. ALCEL = MELALV.VPOCHA(NLCEG, 2)
  806. DCEL = MPGRV.VPOCHA(NLCEG, 3)*DXG +
  807. & MPGRV.VPOCHA(NLCEG, 4)*DYG
  808. UYG = VALCEL + ALCEL * DCEL
  809. C
  810. UNG = UXG * CNX + UYG * CNY
  811. UTG = UXG * CTX + UYG * CTY
  812. C
  813. C********** Positivite
  814. C
  815. C LOGI1 = (PG .LT. 0.0D0) .OR. (ROG .LT. 0.0D0) .OR. LOGI2
  816. C 2010/09/14: AB this destroys the symmetry!
  817. C LOGI1 = .FALSE.
  818. C
  819. IF(LOGI1)THEN
  820. C
  821. C************* Premier ordre en espace local
  822. C
  823. ROG = MPROC.VPOCHA(NLCEG,1)
  824. PG = MPPC.VPOCHA(NLCEG,1)
  825. UNG = MPVITC.VPOCHA(NLCEG,1)*CNX +
  826. & MPVITC.VPOCHA(NLCEG,2)*CNY
  827. UTG = MPVITC.VPOCHA(NLCEG,1)*CTX +
  828. & MPVITC.VPOCHA(NLCEG,2)*CTY
  829. DO I1 = 1, NESP
  830. FRAMAS.FRAMG(I1) = MPYC.VPOCHA(NLCEG,I1)
  831. ENDDO
  832. GAMG = MPGAMC.VPOCHA(NLCEG,1)
  833. ENDIF
  834. C
  835. C********** Etat droite
  836. C
  837. VALCEL = MPROP.VPOCHA(NLCED, 1)
  838. ALCEL = MELALR.VPOCHA(NLCED, 1)
  839. DCEL = MPGRR.VPOCHA(NLCED, 1)*DXD +
  840. & MPGRR.VPOCHA(NLCED, 2)*DYD
  841. ROD = VALCEL + ALCEL * DCEL
  842. C
  843. VALCEL = MPPP.VPOCHA(NLCED, 1)
  844. ALCEL = MELALP.VPOCHA(NLCED, 1)
  845. DCEL = MPGRP.VPOCHA(NLCED, 1)*DXD +
  846. & MPGRP.VPOCHA(NLCED, 2)*DYD
  847. PD = VALCEL + ALCEL * DCEL
  848. C
  849. LOGI2 = .FALSE.
  850. SUMY = 0.0D0
  851. CPCEL = 0.0D0
  852. CVCEL = 0.0D0
  853. DO I1 = 1, NESP
  854. INDCEL = 2 * I1 - 1
  855. VALCEL = MPYP.VPOCHA(NLCED,I1)
  856. ALCEL = MELALY.VPOCHA(NLCED, I1)
  857. DCEL = MPGRY.VPOCHA(NLCED, INDCEL)*DXD +
  858. & MPGRY.VPOCHA(NLCED,INDCEL + 1 )*DYD
  859. ALCEL = VALCEL + ALCEL * DCEL
  860. SUMY = SUMY + ALCEL
  861. CPCEL = CPCEL + MLRECP.PROG(I1) * ALCEL
  862. CVCEL = CVCEL + MLRECV.PROG(I1) * ALCEL
  863. LOGI2 = LOGI2 .OR. (ALCEL .LT. 0.0D0)
  864. FRAMAS.FRAMD(I1) = ALCEL
  865. ENDDO
  866. LOGI2 = LOGI2 .OR. (SUMY .GT. 1.0D0)
  867. CPCEL = CPCEL + (1.0D0 - SUMY) * MLRECP.PROG(NESP+1)
  868. CVCEL = CVCEL + (1.0D0 - SUMY) * MLRECV.PROG(NESP+1)
  869. GAMD = CPCEL / CVCEL
  870. C
  871. VALCEL = MPVITP.VPOCHA(NLCED, 1)
  872. ALCEL = MELALV.VPOCHA(NLCED, 1)
  873. DCEL = MPGRV.VPOCHA(NLCED, 1)*DXD +
  874. & MPGRV.VPOCHA(NLCED, 2)*DYD
  875. UXD = VALCEL + ALCEL * DCEL
  876. C
  877. VALCEL = MPVITP.VPOCHA(NLCED, 2)
  878. ALCEL = MELALV.VPOCHA(NLCED, 2)
  879. DCEL = MPGRV.VPOCHA(NLCED, 3)*DXD +
  880. & MPGRV.VPOCHA(NLCED, 4)*DYD
  881. UYD = VALCEL + ALCEL * DCEL
  882. C
  883. UND = UXD * CNX + UYD * CNY
  884. UTD = UXD * CTX + UYD * CTY
  885. C
  886. C********** Positivite
  887. C
  888. C LOGI1 = (PD .LT. 0.0D0) .OR. (ROD .LT. 0.0D0) .OR. LOGI2
  889. C 2010/09/14: AB this destroys the symmetry!
  890. C LOGI1 = .FALSE.
  891. C
  892. IF(LOGI1)THEN
  893. C
  894. C************* Premier ordre en espace local
  895. C
  896. ROD = MPROC.VPOCHA(NLCED,1)
  897. PD = MPPC.VPOCHA(NLCED,1)
  898. UND = MPVITC.VPOCHA(NLCED,1)*CNX +
  899. & MPVITC.VPOCHA(NLCED,2)*CNY
  900. UTD = MPVITC.VPOCHA(NLCED,1)*CTX +
  901. & MPVITC.VPOCHA(NLCED,2)*CTY
  902. DO I1 = 1, NESP
  903. FRAMAS.FRAMD(I1) = MPYC.VPOCHA(NLCED,I1)
  904. ENDDO
  905. GAMD = MPGAMC.VPOCHA(NLCED,1)
  906. ENDIF
  907. ENDIF
  908. C
  909. C******** Les MELVALs
  910. C
  911. MELRO.VELCHE(1,NLCF) = ROG
  912. MELRO.VELCHE(3,NLCF) = ROD
  913. MELP.VELCHE(1,NLCF) = PG
  914. MELP.VELCHE(3,NLCF) = PD
  915. MELGAM.VELCHE(1,NLCF) = GAMG
  916. MELGAM.VELCHE(3,NLCF) = GAMD
  917. MELVUN.VELCHE(1,NLCF) = UNG
  918. MELVUN.VELCHE(3,NLCF) = UND
  919. MELVUT.VELCHE(1,NLCF) = UTG
  920. MELVUT.VELCHE(3,NLCF) = UTD
  921. MELVNX.VELCHE(1,NLCF) = CNX
  922. MELVNY.VELCHE(1,NLCF) = CNY
  923. MELT1X.VELCHE(1,NLCF) = CTX
  924. MELT1Y.VELCHE(1,NLCF) = CTY
  925. DO I1 = 1, NESP
  926. MELVA1 = MCHAMY.IELVAL(I1)
  927. MELVA1.VELCHE(1,NLCF) = FRAMAS.FRAMG(I1)
  928. MELVA1.VELCHE(3,NLCF) = FRAMAS.FRAMD(I1)
  929. ENDDO
  930. ENDDO
  931. C
  932. C**** Desactivation des SEGMENTs
  933. C
  934. SEGDES IPT1
  935. SEGDES IPT2
  936. C
  937. C**** Le MPOVALs 'Prediction' sont detruits (si existentes)
  938. C
  939. IF(LOGTEM)THEN
  940. SEGSUP MPROP
  941. SEGSUP MPVITP
  942. SEGSUP MPPP
  943. SEGSUP MPYP
  944. ENDIF
  945. C
  946. SEGDES MPROC
  947. SEGDES MPGRR
  948. SEGDES MELALR
  949. SEGDES MPVITC
  950. SEGDES MPGRV
  951. SEGDES MELALV
  952. SEGDES MPPC
  953. SEGDES MPGRP
  954. SEGDES MELALP
  955. SEGDES MPYC
  956. SEGDES MPGRY
  957. SEGDES MELALY
  958. SEGDES MPGAMC
  959. SEGDES MPNORM
  960. C
  961. SEGDES MELRO
  962. SEGDES MELP
  963. SEGDES MELGAM
  964. SEGDES MELVUN
  965. SEGDES MELVUT
  966. SEGDES MELVNX
  967. SEGDES MELVNY
  968. SEGDES MELT1X
  969. SEGDES MELT1Y
  970. C
  971. SEGDES MPYC
  972. DO I1 = 1, NESP
  973. MELVA1 = MCHAMY.IELVAL(I1)
  974. SEGDES MELVA1
  975. ENDDO
  976. SEGDES MCHAMY
  977. SEGSUP FRAMAS
  978. C
  979. SEGDES MLMOT1
  980. CC
  981. C**** Destruction du MELNTI correspondance local/global
  982. C
  983. SEGSUP MLENT1
  984. C
  985. 9999 CONTINUE
  986. C
  987. RETURN
  988. END
  989.  
  990.  
  991.  
  992.  
  993.  
  994.  
  995.  
  996.  
  997.  
  998.  
  999.  

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