Télécharger kodfl1.eso

Retour à la liste

Numérotation des lignes :

  1. C KODFL1 SOURCE BECC 11/06/08 21:15:36 7000
  2. SUBROUTINE KODFL1(INDK0,NESP,NORD,TMAX,RUNIV,PROPHY,
  3. & MLRCHE,MLRMFR,
  4. & INDMET,
  5. & ICHPA1, ICHPA2,
  6. & IALF1,IROF1,IVITF1,IPF1,
  7. & IALF2,IROF2,IVITF2,IPF2,
  8. & IUINF1, IUINF2,
  9. & K0, ICHPK0,IGRALP,EPS,
  10. & ICHPSU,ICHPDI,ICHPVO,
  11. & MELEMC,MELEMF,MELEFE,MELLIM,
  12. & ICHRES,
  13. & DT,SURFL,
  14. & LOGNC,LOGAN,MESERR)
  15. C************************************************************************
  16. C
  17. C PROJET : CASTEM 2000
  18. C
  19. C NOM : KODFL1
  20. C
  21. C DESCRIPTION : Voir KONV13
  22. C
  23. C Cas deux dimensions, DEM
  24. C
  25. C LANGAGE : FORTRAN 77 + ESOPE 2000 (avec estensions CISI)
  26. C
  27. C AUTEUR : A. BECCANTINI, DEN/DM2S/SFME/LTMF
  28. C
  29. C************************************************************************
  30. C
  31. C ENTREES
  32. C
  33. C 0) INDK0 = 1 => K0 = constant
  34. C INDK0 = 2 => K0 = variable
  35. C
  36. C 1) Parameters concerning some gas properties
  37. C
  38. C NESP : nombre d'especes dans le melange.
  39. C
  40. C NORD : ordre des polynoms du cv_i
  41. C
  42. C TMAX : maximum temperature for cv expansion
  43. C
  44. C RUNIV : universal constant for gases
  45. C
  46. C PROPHY : thermodynamic properties for the gases
  47. C
  48. C MLRCHE : LISTREEL with the coefS involved in the chemical
  49. C reaction
  50. C
  51. C MLRMFR : LISTREEL with the mass fractions before or after
  52. C the chemical reaction
  53. C
  54. C 2) Numerical scheme
  55. C
  56. C INDMET : 1 VLH (non-reactive) + SS (reactive)
  57. C
  58. C INDMET : 2 SS (non-reactive) + SS (reactive)
  59. C
  60. C INDMET : 3 AUSM+up (non-reactive) + SS (reactive)
  61. C
  62. C 3) Pointers for the CHPOINTs/MCHAML des variables
  63. C
  64. C ICHPA1 : chpoint CENTRE contenant alpha1
  65. C ICHPA2 : chpoint CENTRE contenant alpha2
  66. C
  67. C IALF1 : MCHAML sur "FACEL" contenant alpha1
  68. C ("gauche" et "droite");
  69. C
  70. C IROF1 : MCHAML sur "FACEL" contenant la masse volumique 1
  71. C ("gauche" et "droite");
  72. C
  73. C IVITF1 : MCHAML sur "FACEL" contenant la vitesse 1 dans le repaire
  74. C local (n,t) et les cosinus directeurs des repaire local;
  75. C
  76. C IPF1 : MCHAML sur "FACEL" contenant la pression 1;
  77. C
  78. C IALF2 : MCHAML sur "FACEL" contenant alpha1
  79. C ("gauche" et "droite");
  80. C
  81. C IROF2 : MCHAML sur "FACEL" contenant la masse volumique 1
  82. C ("gauche" et "droite");
  83. C
  84. C IVITF2 : MCHAML sur "FACEL" contenant la vitesse 1 dans le repaire
  85. C local (n,t) et les cosinus directeurs des repaire local;
  86. C
  87. C IPF2 : MCHAML sur "FACEL" contenant la pression 1;
  88. C
  89. C IUINF1 : CHPOINT, one component, "cut-off velocity" 1
  90. C 0 if no Bas MACH
  91. C
  92. C IUINF2 : CHPOINT, one component, "cut-off velocity" 2
  93. C
  94. C
  95. C 4) Others
  96. C
  97. C K0 or
  98. C ICHPK0 : fundamental flame speed
  99. C
  100. C IGRALP : CHPOINT "FACE" contenant grad(alpha)/alpha
  101. C
  102. C EPS : "zero"
  103. C
  104. C ICHPSU : CHPOINT "FACE" contenant la surface des faces
  105. C
  106. C ICHPDI : CHPOINT "CENTRE" contenant le diametre minimum
  107. C de chaque element
  108. C
  109. C ICHPVO : CHPOINT "CENTRE" contenant le volume
  110. C de chaque element
  111. C
  112. C 5) Pointeurs de MELEME de la table DOMAINE
  113. C
  114. C MELEMC : MELEME 'CENTRE' du SPG des CENTRES
  115. C
  116. C MELEMF : MELEME 'FACE' du SPG des FACES
  117. C
  118. C MELEFE : MELEME 'FACEL' du connectivité FACES -> ELEM
  119. C
  120. C MELLIM : MAILLAGE where the flux (or residu) will not be found
  121. C
  122. C SORTIES (il faudrait dire E/S)
  123. C
  124. C ICHRES : pointeurs de CHPOINTs "FACE" des flux aux interfaces:
  125. C
  126. C DT : pas de temps pour le respect de la CFL-like condition
  127. C DT < DIAMEL /2 /max(Lambda_i)
  128. C En maillage regulier cette condition garantie la
  129. C non-interaction des ondes
  130. C
  131. C LOGNC : (LOGICAL): si .TRUE. la methode de Newton-Rapson, utilisée
  132. C dans pour la solution du probleme Riemann n'a pas bien
  133. C marchéee
  134. C
  135. C LOGAN : (LOGICAL): si .TRUE. une anomalie à été detectée
  136. C
  137. C MESERR : pour l'ecriture des messages d'erreurs
  138. C
  139. C************************************************************************
  140. C
  141. C HISTORIQUE (Anomalies et modifications éventuelles)
  142. C
  143. C 07.12.09 Creation
  144. C 21.12.2010 Estension au 3D
  145. C
  146. C************************************************************************
  147. C
  148. C
  149. C N.B.: On suppose qu'on a déjà controllé RO, P > 0
  150. C Y \in (0,1)
  151. C Si non il faut le faire!!!
  152. C
  153. C************************************************************************
  154. C
  155. C
  156. IMPLICIT INTEGER(I-N)
  157. INTEGER I1, I2, ICOMP
  158. & , NESP,NORD ,INDMET
  159. & , ICHPA1, ICHPA2
  160. & , IALF1, IROF1, IVITF1, IPF1
  161. & , IALF2, IROF2, IVITF2, IPF2
  162. & , IUINF1, IUINF2
  163. & , IGRALP, ICHPK0
  164. & , ICHPSU, ICHPDI, ICHPVO, MELEMC, MELEMF, MELEFE, MELLIM
  165. & , IGEOMC, IGEOMF
  166. & , ICHRES
  167. & , NFAC
  168. & , NLCF, NGCEG, NGCED, NLCEG, NLCED
  169. & , NGCF, NLCF1, SPG1, SPG2, NLFL
  170. & , NLHS
  171. REAL*8 TMAX, RUNIV, K0, K0G, K0D, EPS
  172. REAL*8 DT, UNSDT, CELLT, CELLTM
  173. & , AG1C, ALPG1, ROG1, UNG1, UTG1, UVG1, PG1
  174. & , AG2C, ALPG2, ROG2, UNG2, UTG2, UVG2, PG2
  175. & , AD1C, ALPD1, ROD1, UND1, UTD1, UVD1, PD1
  176. & , AD2C, ALPD2, ROD2, UND2, UTD2, UVD2, PD2
  177. & , SURF,CNX, CNY, CNZ, CTX , CTY, CTZ, CVX, CVY, CVZ
  178. & , V_INF
  179. & , GRALX, GRALY, GRALZ
  180. & , CNNA, CTNA, CVNA
  181. & , CELL, DIAMG, DIAMD, DIAM, VOLG, VOLD
  182. & , FLUCEL, SCON1, SG1, SD1, SCON2, SG2, SD2
  183. & , S1_11, S2_22, S1_12, S1_21, S2_12, S2_21
  184. & , S1_21l, S1_12l, S2_21l, S2_12l
  185. & , S1_21r, S1_12r, S2_21r, S2_12r
  186. & , S1_21e, S1_12e, S2_21e, S2_12e
  187. & , RESG(12), RESD(12)
  188. & , SURFL, SURF0
  189. LOGICAL LOGNC, LOGAN
  190. CHARACTER*(40) MESERR
  191. CHARACTER*(8) TYPE
  192. C
  193. C**** Segment des proprietes du gaz
  194. C
  195. SEGMENT PROPHY
  196. REAL*8 XTAB(N1,N2)
  197. C REAL*8 ACV(NORDP1,NESP+1), W(NESP+1), H0K(NESP+1)
  198. C & ,ACVTOG(NORDP1), ACVTOD(NORDP1)
  199. ENDSEGMENT
  200. C N1 = NORD + 3
  201. C N2 = NESP
  202. C
  203. C**** LES INCLUDES
  204. C
  205. -INC CCOPTIO
  206. -INC SMCOORD
  207. -INC SMELEME
  208. -INC SMCHAML
  209. C Normal and tangent
  210. POINTEUR MELVNX.MELVAL, MELVNY.MELVAL, MELVNZ.MELVAL,
  211. & MELT1X.MELVAL, MELT1Y.MELVAL, MELT1Z.MELVAL,
  212. & MELT2X.MELVAL, MELT2Y.MELVAL, MELT2Z.MELVAL
  213. POINTEUR MELUN1.MELVAL, MELUT1.MELVAL, MELUV1.MELVAL
  214. POINTEUR MELAL1.MELVAL, MELRO1.MELVAL, MELP1.MELVAL
  215. POINTEUR MELUN2.MELVAL, MELUT2.MELVAL, MELUV2.MELVAL
  216. POINTEUR MELAL2.MELVAL, MELRO2.MELVAL, MELP2.MELVAL
  217. -INC SMCHPOI
  218. POINTEUR MPOVSU.MPOVAL, MPOVDI.MPOVAL, MPOVVO.MPOVAL
  219. & , MPORES.MPOVAL, MPOALP.MPOVAL
  220. & , MPALP1.MPOVAL, MPALP2.MPOVAL
  221. & , MPUIN2.MPOVAL, MPUIN1.MPOVAL
  222. & , MPOK0.MPOVAL
  223. -INC SMLMOTS
  224. -INC SMLENTI
  225. POINTEUR MLELIM.MLENTI
  226. INTEGER JG
  227. -INC SMLREEL
  228. POINTEUR MLRMFR.MLREEL, MLRCHE.MLREEL
  229. C
  230. C**** Segment for the subroutine prithe
  231. C
  232. POINTEUR IWST_L.MLREEL, IWST_R.MLREEL, IWWORK.MLREEL
  233. C 1:(4+IDIM+NESP)
  234. POINTEUR IY.MLREEL
  235. C 1:NESP
  236. POINTEUR ICVL.MLREEL, ICVR.MLREEL, ICVRSS.MLREEL, ICVWOR.MLREEL
  237. C 0:NORD
  238. C
  239. POINTEUR IFLU.MLREEL, IFLULA.MLREEL
  240. C
  241. C**** Surface de flamme totale
  242. C
  243. SURFL = 0.0D0
  244. C
  245. C**** KRIPAD pour la correspondance global/local des conditions limits
  246. C
  247. CALL KRIPAD(MELLIM,MLELIM)
  248. C SEGINI MLELIM
  249. C
  250. C**** Initialisation des MCHAMLs
  251. C
  252. C Alpha
  253. C
  254. MCHEL1 = IALF1
  255. SEGACT MCHEL1
  256. MCHAM1 = MCHEL1.ICHAML(1)
  257. SEGACT MCHAM1
  258. MELAL1 = MCHAM1.IELVAL(1)
  259. SEGDES MCHEL1
  260. SEGDES MCHAM1
  261. C
  262. C Masse volumique
  263. C
  264. MCHEL1 = IROF1
  265. SEGACT MCHEL1
  266. MCHAM1 = MCHEL1.ICHAML(1)
  267. SEGACT MCHAM1
  268. MELRO1 = MCHAM1.IELVAL(1)
  269. SEGDES MCHEL1
  270. SEGDES MCHAM1
  271. C
  272. C Pression
  273. C
  274. MCHEL1 = IPF1
  275. SEGACT MCHEL1
  276. MCHAM1 = MCHEL1.ICHAML(1)
  277. SEGACT MCHAM1
  278. MELP1 = MCHAM1.IELVAL(1)
  279. SEGDES MCHEL1
  280. SEGDES MCHAM1
  281. C
  282. C Vitesse et cosinus directeurs du repere (n,t)
  283. C
  284. MCHEL1 = IVITF1
  285. SEGACT MCHEL1
  286. C
  287. C La vitesse a comme SPG MELEFE
  288. C Le cosinus directeurs ont comme SPG MELEMF
  289. C
  290. C MCHAM1 -> Cosinus directeurs
  291. C MCHAM2 -> Vitesse
  292. C
  293. SPG1 = MCHEL1.IMACHE(1)
  294. SPG2 = MCHEL1.IMACHE(2)
  295. IF((SPG1 .EQ. MELEMF) .AND. (SPG2 .EQ. MELEFE))THEN
  296. MCHAM1 = MCHEL1.ICHAML(1)
  297. MCHAM2 = MCHEL1.ICHAML(2)
  298. ELSEIF((SPG1 .EQ. MELEFE) .AND. (SPG2 .EQ. MELEMF))THEN
  299. MCHAM1 = MCHEL1.ICHAML(2)
  300. MCHAM2 = MCHEL1.ICHAML(1)
  301. ELSE
  302. LOGAN = .TRUE.
  303. GOTO 9999
  304. ENDIF
  305. SEGACT MCHAM1
  306. SEGACT MCHAM2
  307. IF (IDIM .EQ. 2) THEN
  308. MELVNX = MCHAM1.IELVAL(1)
  309. MELVNY = MCHAM1.IELVAL(2)
  310. MELT1X = MCHAM1.IELVAL(3)
  311. MELT1Y = MCHAM1.IELVAL(4)
  312. MELUN1 = MCHAM2.IELVAL(1)
  313. MELUT1 = MCHAM2.IELVAL(2)
  314. ELSE
  315. MELVNX = MCHAM1.IELVAL(1)
  316. MELVNY = MCHAM1.IELVAL(2)
  317. MELVNZ = MCHAM1.IELVAL(3)
  318. MELT1X = MCHAM1.IELVAL(4)
  319. MELT1Y = MCHAM1.IELVAL(5)
  320. MELT1Z = MCHAM1.IELVAL(6)
  321. MELT2X = MCHAM1.IELVAL(7)
  322. MELT2Y = MCHAM1.IELVAL(8)
  323. MELT2Z = MCHAM1.IELVAL(9)
  324. MELUN1 = MCHAM2.IELVAL(1)
  325. MELUT1 = MCHAM2.IELVAL(2)
  326. MELUV1 = MCHAM2.IELVAL(3)
  327. ENDIF
  328. SEGDES MCHAM1
  329. SEGDES MCHAM2
  330. SEGDES MCHEL1
  331. C
  332. C Alpha
  333. C
  334. MCHEL1 = IALF2
  335. SEGACT MCHEL1
  336. MCHAM1 = MCHEL1.ICHAML(1)
  337. SEGACT MCHAM1
  338. MELAL2 = MCHAM1.IELVAL(1)
  339. SEGDES MCHEL1
  340. SEGDES MCHAM1
  341. C
  342. C Masse volumique
  343. C
  344. MCHEL1 = IROF2
  345. SEGACT MCHEL1
  346. MCHAM1 = MCHEL1.ICHAML(1)
  347. SEGACT MCHAM1
  348. MELRO2 = MCHAM1.IELVAL(1)
  349. SEGDES MCHEL1
  350. SEGDES MCHAM1
  351. C
  352. C Pression
  353. C
  354. MCHEL1 = IPF2
  355. SEGACT MCHEL1
  356. MCHAM1 = MCHEL1.ICHAML(1)
  357. SEGACT MCHAM1
  358. MELP2 = MCHAM1.IELVAL(1)
  359. SEGDES MCHEL1
  360. SEGDES MCHAM1
  361. C
  362. C Vitesse et cosinus directeurs du repere (n,t)
  363. C
  364. MCHEL1 = IVITF2
  365. SEGACT MCHEL1
  366. C
  367. C La vitesse a comme SPG MELEFE
  368. C Le cosinus directeurs ont comme SPG MELEMF
  369. C
  370. C MCHAM1 -> Cosinus directeurs
  371. C MCHAM2 -> Vitesse
  372. C
  373. SPG1 = MCHEL1.IMACHE(1)
  374. SPG2 = MCHEL1.IMACHE(2)
  375. IF((SPG1 .EQ. MELEMF) .AND. (SPG2 .EQ. MELEFE))THEN
  376. MCHAM1 = MCHEL1.ICHAML(1)
  377. MCHAM2 = MCHEL1.ICHAML(2)
  378. ELSEIF((SPG1 .EQ. MELEFE) .AND. (SPG2 .EQ. MELEMF))THEN
  379. MCHAM1 = MCHEL1.ICHAML(2)
  380. MCHAM2 = MCHEL1.ICHAML(1)
  381. ELSE
  382. LOGAN = .TRUE.
  383. GOTO 9999
  384. ENDIF
  385. C SEGACT MCHAM1
  386. C MELVNX = MCHAM1.IELVAL(1)
  387. C MELVNY = MCHAM1.IELVAL(2)
  388. C MELT1X = MCHAM1.IELVAL(3)
  389. C MELT1Y = MCHAM1.IELVAL(4)
  390. C SEGDES MCHAM1
  391. SEGACT MCHAM2
  392. MELUN2 = MCHAM2.IELVAL(1)
  393. MELUT2 = MCHAM2.IELVAL(2)
  394. IF (IDIM .EQ. 3) MELUV2 = MCHAM2.IELVAL(3)
  395. SEGDES MCHAM2
  396. SEGDES MCHEL1
  397. C write(*,*) 'Son in kodfl1'
  398. C write(*,*) 'Ho letto gli MCHAM'
  399. C
  400. C**** Segment for fcolre
  401. C
  402. C Chemical coefficients in MLRCHE.MLREEL.
  403. SEGACT MLRCHE
  404. NLHS = 0
  405. DO I1 = 1, NESP, 1
  406. IF (MLRCHE.PROG(I1) .GT. 0.0D0) NLHS = NLHS + 1
  407. C write(*,*) MLRCHE.PROG(I1)
  408. ENDDO
  409. c write(*,*) NLHS
  410. C
  411. C Mass fractions
  412. C MLRMFR.PROG(1) = YINI_first
  413. C MLRMFR.PROG(2) = YFIN_first
  414. C MLRMFR.PROG(3) = YFIN_second
  415. C ...
  416. C MLRMFR.PROG(nLHS+1) = YFIN_nLHS
  417. C MLRMFR.PROG(nLHS+2) = YINI_(nLHS+1)
  418. C
  419. SEGACT MLRMFR
  420. C
  421. C cv, molar weicht formation energy in PROPHY.XTAB,
  422. C already active
  423. C
  424. C States
  425. JG = NESP + IDIM + 4
  426. SEGINI IWWORK
  427. SEGINI IWST_R
  428. SEGINI IWST_L
  429. C
  430. C Mass fractions
  431. JG = NESP
  432. SEGINI IY
  433. C
  434. C Coeff for CV
  435. JG = NORD + 1
  436. SEGINI ICVL
  437. SEGINI ICVR
  438. SEGINI ICVRSS
  439. SEGINI ICVWOR
  440. C
  441. C Fluxes
  442. JG = NESP + IDIM + 4
  443. SEGINI IFLU
  444. SEGINI IFLULA
  445. C
  446. C Since the mixture here considered is homogeneous, we first fill
  447. C the part of IWST_L and IWST_R involving the species mass fractions
  448. C and K0
  449. C
  450. IWST_L.PROG(4 + IDIM) = MLRMFR.PROG(2)
  451. IWST_R.PROG(4 + IDIM) = MLRMFR.PROG(2)
  452. DO I2 = 3, NESP, 1
  453. I1 = 2 + IDIM + I2
  454. C DO I1 = 5 + IDIM, 2 + IDIM + NESP
  455. IWST_L.PROG(I1) = MLRMFR.PROG(I2)
  456. IWST_R.PROG(I1) = MLRMFR.PROG(I2)
  457. ENDDO
  458. IWST_L.PROG(3 + IDIM + NESP) = MLRMFR.PROG(1) - MLRMFR.PROG(2)
  459. IWST_R.PROG(3 + IDIM + NESP) = MLRMFR.PROG(1) - MLRMFR.PROG(2)
  460. C
  461. C
  462. C**** Initialisation des MELEMEs
  463. C
  464. C 'CENTRE', 'FACEL'
  465. C
  466. IPT2 = MELEFE
  467. SEGACT IPT2
  468. NFAC = IPT2.NUM(/2)
  469. C
  470. C**** KRIPAD pour la correspondance global/local de centre
  471. C
  472. CALL KRIPAD(MELEMC,MLENT1)
  473. C
  474. C**** MLENTI1 a MCORD.XCOORD(/1)/(IDIM+1) elements
  475. C
  476. C Si i est le numero global d'un noeud de ICEN,
  477. C MLENT1.LECT(i) contient sa position, i.e.
  478. C
  479. C I = numero global du noeud centre
  480. C MLENT1.LECT(i) = numero local du noeud centre
  481. C
  482. C MLENT1 déjà activé, i.e.
  483. C
  484. C SEGINI MLENT1
  485. C
  486. C
  487. C**** KRIPAD pour la correspondance global/local de 'FACE'
  488. C
  489. CALL KRIPAD(MELEMF,MLENT2)
  490. C SEGINI MLENT2
  491. C
  492. C
  493. C**** CHPOINTs de la normale à alpha
  494. C
  495. CALL LICHT(IGRALP,MPOALP,TYPE,IGEOMF)
  496. IF (INDK0 .EQ. 2) CALL LICHT(ICHPK0,MPOK0,TYPE,IGEOMF)
  497. C SEGACT MPOALP
  498. C SEGDES MPOK0
  499. C
  500. C**** CHPOINTs de la table DOMAINE
  501. C
  502. CALL LICHT(ICHPSU,MPOVSU,TYPE,IGEOMF)
  503. CALL LICHT(ICHPDI,MPOVDI,TYPE,IGEOMC)
  504. CALL LICHT(ICHPVO,MPOVVO,TYPE,IGEOMC)
  505. C
  506. C**** LICHT active les MPOVALs en *MOD
  507. C
  508. C i.e.
  509. C
  510. C SEGACT MPOVSU*MOD
  511. C SEGACT MPOVDI*MOD
  512. C SEGACT MPOVVO*MOD
  513. C
  514. C
  515. C**** Les FLUX aux face
  516. C
  517. CALL LICHT(ICHRES,MPORES,TYPE,IGEOMC)
  518. C SEGACT MPORES*MOD
  519. C
  520. C**** Activation des MCHAMLs
  521. C
  522. SEGACT MELAL1
  523. SEGACT MELRO1
  524. SEGACT MELP1
  525. SEGACT MELUN1
  526. SEGACT MELUT1
  527. SEGACT MELAL2
  528. SEGACT MELRO2
  529. SEGACT MELP2
  530. SEGACT MELUN2
  531. SEGACT MELUT2
  532. SEGACT MELVNX
  533. SEGACT MELVNY
  534. SEGACT MELT1X
  535. SEGACT MELT1Y
  536. IF (IDIM .EQ. 3) THEN
  537. SEGACT MELUV1
  538. SEGACT MELUV2
  539. SEGACT MELVNZ
  540. SEGACT MELT1Z
  541. SEGACT MELT2X
  542. SEGACT MELT2Y
  543. SEGACT MELT2Z
  544. ENDIF
  545. C
  546. C**** Lecture des MPOVALs ICHPA1, ICHPA2
  547. C
  548. CALL LICHT(ICHPA1, MPALP1, TYPE, IGEOMC)
  549. CALL LICHT(ICHPA2, MPALP2, TYPE, IGEOMC)
  550. C SEGACT MPALP1
  551. C SEGACT MPALP2
  552. C
  553. C**** Bas Mach
  554. C
  555. IF(IUINF1 .NE. 0)THEN
  556. CALL LICHT(IUINF2,MPUIN2,TYPE,IGEOMC)
  557. CALL LICHT(IUINF1,MPUIN1,TYPE,IGEOMC)
  558. C SEGACT MPUIN2*MOD
  559. C SEGACT MPUIN1*MOD
  560. ELSE
  561. MPUIN2=0
  562. MPUIN1=0
  563. ENDIF
  564. C
  565. C**** Initialisation de 1/DT et d'autres variables
  566. C
  567. UNSDT = 0.0D0
  568. GRALZ = 0.0D0
  569. CNZ = 0.0D0
  570. CTZ = 0.0D0
  571. CVX = 0.0D0
  572. CVY = 0.0D0
  573. CVZ = 0.0D0
  574. UVG1 = 0.0D0
  575. UVD1 = 0.0D0
  576. UVG2 = 0.0D0
  577. UVD2 = 0.0D0
  578. C
  579. C**** BOUCLE SUR FACEL pour le calcul du FLUX/RESIDU
  580. C
  581. DO NLCF = 1, NFAC
  582. C
  583. C******* NLCF = numero local du centre de facel
  584. C NGCF = numero global du centre de facel
  585. C NLCF1 = numero local du centre de face
  586. C NGCEG = numero global du centre ELT "gauche"
  587. C NLCEG = numero local du centre ELT "gauche"
  588. C NGCED = numero global du centre ELT "droite"
  589. C NLCED = numero local du centre ELT "droite"
  590. C
  591. NGCEG = IPT2.NUM(1,NLCF)
  592. NGCED = IPT2.NUM(3,NLCF)
  593. NGCF = IPT2.NUM(2,NLCF)
  594. NLCF1 = MLENT2.LECT(NGCF)
  595. NLCEG = MLENT1.LECT(NGCEG)
  596. NLCED = MLENT1.LECT(NGCED)
  597. NLFL = MLELIM.LECT(NGCF)
  598. C
  599. VOLG = MPOVVO.VPOCHA(NLCEG,1)
  600. VOLD = MPOVVO.VPOCHA(NLCED,1)
  601. SURF0 = MPOVSU.VPOCHA(NLCF,1)
  602. C
  603. C******* NLCF != NLCF1 -> l'auteur (MOI) n'a rien compris.
  604. C
  605. IF(NLCF .NE. NLCF1)THEN
  606. MESERR = 'Il ne faut pas jouer avec la console. '
  607. LOGAN = .TRUE.
  608. GOTO 9999
  609. ENDIF
  610. IF(NLFL .EQ. 0)THEN
  611. CELLTM = 0.0D0
  612. C
  613. C********** Recuperation des Etats "gauche" et "droite"
  614. C
  615. AG1C = MPALP1.VPOCHA(NLCEG,1)
  616. ALPG1 = MELAL1.VELCHE(1,NLCF)
  617. ROG1 = MELRO1.VELCHE(1,NLCF)
  618. UNG1 = MELUN1.VELCHE(1,NLCF)
  619. UTG1 = MELUT1.VELCHE(1,NLCF)
  620. PG1 = MELP1.VELCHE(1,NLCF)
  621. AG2C = MPALP2.VPOCHA(NLCEG,1)
  622. ALPG2 = MELAL2.VELCHE(1,NLCF)
  623. ROG2 = MELRO2.VELCHE(1,NLCF)
  624. UNG2 = MELUN2.VELCHE(1,NLCF)
  625. UTG2 = MELUT2.VELCHE(1,NLCF)
  626. PG2 = MELP2.VELCHE(1,NLCF)
  627. C
  628. AD1C = MPALP1.VPOCHA(NLCED,1)
  629. ALPD1 = MELAL1.VELCHE(3,NLCF)
  630. ROD1 = MELRO1.VELCHE(3,NLCF)
  631. UND1 = MELUN1.VELCHE(3,NLCF)
  632. UTD1 = MELUT1.VELCHE(3,NLCF)
  633. PD1 = MELP1.VELCHE(3,NLCF)
  634. AD2C = MPALP2.VPOCHA(NLCED,1)
  635. ALPD2 = MELAL2.VELCHE(3,NLCF)
  636. ROD2 = MELRO2.VELCHE(3,NLCF)
  637. UND2 = MELUN2.VELCHE(3,NLCF)
  638. UTD2 = MELUT2.VELCHE(3,NLCF)
  639. PD2 = MELP2.VELCHE(3,NLCF)
  640. C
  641. CNX = MELVNX.VELCHE(1,NLCF)
  642. CNY = MELVNY.VELCHE(1,NLCF)
  643. CTX = MELT1X.VELCHE(1,NLCF)
  644. CTY = MELT1Y.VELCHE(1,NLCF)
  645. C
  646. GRALX = MPOALP.VPOCHA(NLCF,1)
  647. GRALY = MPOALP.VPOCHA(NLCF,2)
  648. C
  649. C GRALX = CNX
  650. C GRALY = CNY
  651. C
  652. IF (IDIM .EQ. 3) THEN
  653. GRALZ = MPOALP.VPOCHA(NLCF,3)
  654. CNZ = MELVNZ.VELCHE(1,NLCF)
  655. CTZ = MELT1Z.VELCHE(1,NLCF)
  656. CVX = MELT2X.VELCHE(1,NLCF)
  657. CVY = MELT2Y.VELCHE(1,NLCF)
  658. CVZ = MELT2Z.VELCHE(1,NLCF)
  659. UVG1 = MELUV1.VELCHE(1,NLCF)
  660. UVD1 = MELUV1.VELCHE(3,NLCF)
  661. UVG2 = MELUV2.VELCHE(1,NLCF)
  662. UVD2 = MELUV2.VELCHE(3,NLCF)
  663. ENDIF
  664. C
  665. CNNA = (GRALX * CNX) + (GRALY * CNY) + (GRALZ * CNZ)
  666. CTNA = (GRALX * CTX) + (GRALY * CTY) + (GRALZ * CTZ)
  667. CVNA = (GRALX * CVX) + (GRALY * CVY) + (GRALZ * CVZ)
  668. C
  669. C if ((ngcf .eq. 90) .or. (ngcf .eq. 84) .or. (ngcf .eq. 110))
  670. C $ then
  671. C write(*,*) 'ngcf = ', ngcf
  672. C endif
  673. C
  674. C write(*,*) cnna
  675. C
  676. C write(*,*) alpg1, rog1, ung1, utg1, pg1
  677. C write(*,*) alpg2, rog2, ung2, utg2, pg2
  678. C write(*,*) alpd1, rod1, und1, utd1, pd1
  679. C write(*,*) alpd2, rod2, und2, utd2, pd2
  680. C
  681. C
  682. C
  683. C********** Computation of the partition of the surface
  684. C
  685. S1_11 = min (ALPG1, ALPD1)
  686. S2_22 = min (ALPG2, ALPD2)
  687. C S2_22 = 1.0d0 - (max (ALPG1, ALPD1))
  688. S1_21 = max (0.0d0, (AD1C - AG1C))
  689. S1_21r = max (0.0d0, (AD1C - ALPD1))
  690. S1_21e = max (0.0d0, (ALPD1 - ALPG1))
  691. S1_21l = max (0.0d0, (ALPG1 - AG1C))
  692. S2_21 = max (0.0d0, (AG2C - AD2C))
  693. S2_21l = max (0.0d0, (AG2C - ALPG2))
  694. S2_21e = max (0.0d0, (ALPG2 - ALPD2))
  695. S2_21r = max (0.0d0, (ALPD2 - AD2C))
  696. C
  697. C S1_21 = min (S1_21, S2_21)
  698. C S2_21 = S1_21
  699. S1_12 = max (0.0d0, (AG1C - AD1C))
  700. S1_12l = max (0.0d0, (AG1C - ALPG1))
  701. S1_12e = max (0.0d0, (ALPG1 - ALPD1))
  702. S1_12r = max (0.0d0, (ALPD1 - AD1C))
  703. S2_12 = max (0.0d0, (AD2C - AG2C))
  704. S2_12r = max (0.0d0, (AD2C - ALPD2))
  705. S2_12e = max (0.0d0, (ALPD2 - ALPG2))
  706. S2_12l = max (0.0d0, (ALPG2 - AG2C))
  707. C
  708. C write(*,*) 'NGCF =', NGCF
  709. C write(*,*) (S1_11 + S2_22 + S1_12 + S1_21)
  710. C write(*,*) (S1_11 + S2_22 + S2_12 + S2_21)
  711. *
  712. C write(*,*) (S1_11 + S2_22 + S1_12 + S1_21)
  713. C write(*,*) (S1_11 + S2_22 + S2_12 + S2_21)
  714. c CNNA = 1.0D0
  715. c CTNA = 0.0D0
  716. c CVNA = 0.0D0
  717. C
  718. C SURFL = Surface which burns * (CNNA)
  719. *
  720. SURFL = SURFL + ((ABS ((S1_12 + S1_21) * CNNA)) * SURF0)
  721. C
  722. C Flame speed
  723. C
  724. IF (INDK0 .EQ. 1) THEN
  725. K0G = K0
  726. K0D = K0
  727. ELSE
  728. K0G = MPOK0.VPOCHA(NLCEG,1)
  729. K0D = MPOK0.VPOCHA(NLCED,1)
  730. ENDIF
  731. IWST_L.PROG(4 + IDIM + NESP) = K0G
  732. IWST_R.PROG(4 + IDIM + NESP) = K0D
  733. C
  734. C********** Riemann problem 1-1
  735. C
  736. IF (S1_11 .GE. EPS) THEN
  737. IWST_L.PROG(1) = ROG1
  738. IWST_L.PROG(2) = UNG1
  739. C IWST_L.PROG(2) = 0.0D0
  740. IWST_L.PROG(3) = UTG1
  741. C IDIM + 2 = 2 in 2D case
  742. IWST_L.PROG(IDIM + 2)= PG1
  743. IWST_L.PROG(IDIM + 3) = 0.0D0
  744. C
  745. IWST_R.PROG(1) = ROD1
  746. IWST_R.PROG(2) = UND1
  747. C IWST_R.PROG(2) = 0.0D0
  748. IWST_R.PROG(3) = UTD1
  749. IWST_R.PROG(IDIM + 2) = PD1
  750. IWST_R.PROG(IDIM + 3) = 0.0D0
  751. IF (IDIM .EQ. 3) THEN
  752. IWST_L.PROG(IDIM + 1) = UVG1
  753. IWST_R.PROG(IDIM + 1) = UVD1
  754. ENDIF
  755. C
  756. C do i2 = 1 , 4 + NESP + IDIM, 1
  757. C write(*,*) i2, iwst_l.prog(i2), iwst_r.prog(i2)
  758. C enddo
  759. C
  760. IF (INDMET .EQ. 1) THEN
  761. CALL FVLHRE(IDIM, NESP, NLHS, NORD,
  762. & MLRCHE.PROG, PROPHY.XTAB, RUNIV,
  763. & TMAX,
  764. & CNNA, CTNA, CVNA,
  765. & IWST_L.PROG,
  766. & IWST_R.PROG,
  767. & IWWORK.PROG, IY.PROG, ICVL.PROG, ICVR.PROG,
  768. & ICVRSS.PROG, ICVWOR.PROG,
  769. & IFLU.PROG, IFLULA.PROG, CELLT,
  770. & LOGAN)
  771. ELSEIF(INDMET .EQ. 2)THEN
  772. CALL FBECRE(IDIM, NESP, NLHS, NORD,
  773. & MLRCHE.PROG, PROPHY.XTAB, RUNIV,
  774. & TMAX,
  775. & CNNA, CTNA, CVNA,
  776. & IWST_L.PROG,
  777. & IWST_R.PROG,
  778. & IWWORK.PROG, IY.PROG, ICVL.PROG, ICVR.PROG,
  779. & ICVRSS.PROG, ICVWOR.PROG,
  780. & IFLU.PROG, IFLULA.PROG, CELLT,
  781. & LOGAN)
  782. ELSEIF(INDMET .EQ. 3)THEN
  783. C
  784. C********** AUSMP-up
  785. C
  786. V_INF=MAX(MPUIN1.VPOCHA(NLCEG,1),
  787. & MPUIN1.VPOCHA(NLCED,1))
  788. CALL FAUSRE(IDIM, NESP, NLHS, NORD,
  789. & MLRCHE.PROG, PROPHY.XTAB, RUNIV,
  790. & TMAX,V_INF,
  791. & CNNA, CTNA, CVNA,
  792. & IWST_L.PROG,
  793. & IWST_R.PROG,
  794. & IWWORK.PROG, IY.PROG, ICVL.PROG, ICVR.PROG,
  795. & ICVRSS.PROG, ICVWOR.PROG,
  796. & IFLU.PROG, IFLULA.PROG, CELLT,
  797. & LOGAN)
  798. ENDIF
  799.  
  800. C write(*,*) 'Riemann problem 1-1'
  801. C write(*,*) iflula.prog(idim+1), iflu.prog(idim+1)
  802. C
  803. C subroutine fcolre(ndim, nesp, nLHS, nordpo,
  804. C & reacoe, aprop, Runiv,
  805. C & Tmaxcv,
  806. C & wvect_l,
  807. C & wvect_r,
  808. C & wwork, y, acv_l, acv_r, acv_rss, acv_work,
  809. C & flu, flu_lag, cellt,
  810. C & logan)
  811. C
  812. C
  813. IF (LOGAN) THEN
  814. WRITE(IOIMP,*) 'SUBROUTINE KODFL1'
  815. WRITE(IOIMP,*) 'RIEMANN PROBLEM 1-1'
  816. WRITE(IOIMP,*) 'FACE ', NGCF, ' S11 ', S1_11
  817. WRITE(IOIMP,*) ' ICOMP WVECT_L WVECT_R'
  818. DO I2 = 1, 2 + IDIM
  819. WRITE(IOIMP,'(I6,2E16.6)') I2, IWST_L.PROG(I2),
  820. $ IWST_R.PROG(I2)
  821. ENDDO
  822. WRITE(IOIMP,*) 'ANOMALY DETECTED'
  823. GOTO 9999
  824. ENDIF
  825. C
  826. CELLTM = MAX(CELLT, CELLTM)
  827. C
  828. C************* Ecriture du residu
  829. C
  830. C FLUX(2) = RO Un RO Un
  831. C FLUX(3) = RO Un Un + P -> RO Un Ux + P CNX
  832. C FLUX(4) = RO Un Ut -> RO Un Uy + P CNY
  833. C FLUX(2 + IDIM) = RO Un Et RO Un Et
  834. C
  835. SURF = SURF0 * S1_11
  836. FLUCEL = IFLU.PROG(1) * SURF
  837. C write(*,*) 'flux', iflu.prog(1)
  838. C write(*,*) 'nlceg, nlced', nlceg, nlced
  839. C write(*,*) 'surf', surf
  840. C
  841. C Rho un
  842. C
  843. MPORES.VPOCHA(NLCEG,2) = MPORES.VPOCHA(NLCEG,2) -
  844. & (FLUCEL / VOLG)
  845. IF (NLCED .NE. NLCEG)
  846. & MPORES.VPOCHA(NLCED,2) = MPORES.VPOCHA(NLCED,2) +
  847. & (FLUCEL / VOLD)
  848. C
  849. C Rho (Un Ux + P CNX)
  850. C NB : In 2D (IFLU.PROG(4)*CVX) = RO Un Et * 0 = 0
  851. C In 3D (IFLU.PROG(4)*CVX) = RO Un uv * CVX
  852. C
  853. FLUCEL = ( (IFLU.PROG(2)*CNX) + (IFLU.PROG(3)*CTX) +
  854. & (IFLU.PROG(4)*CVX)) * SURF
  855. MPORES.VPOCHA(NLCEG,3) = MPORES.VPOCHA(NLCEG,3) -
  856. & (FLUCEL / VOLG)
  857. IF (NLCED .NE. NLCEG)
  858. & MPORES.VPOCHA(NLCED,3) = MPORES.VPOCHA(NLCED,3) +
  859. & (FLUCEL / VOLD)
  860. C
  861. C Rho (Un Uy + P CNY)
  862. C
  863. FLUCEL = ( (IFLU.PROG(2)*CNY) + (IFLU.PROG(3)*CTY)
  864. & + (IFLU.PROG(4)*CVY) ) * SURF
  865. MPORES.VPOCHA(NLCEG,4) = MPORES.VPOCHA(NLCEG,4) -
  866. & (FLUCEL / VOLG)
  867. IF (NLCED .NE. NLCEG)
  868. & MPORES.VPOCHA(NLCED,4) = MPORES.VPOCHA(NLCED,4) +
  869. & (FLUCEL / VOLD)
  870. C
  871. C Rho (Un Uz + P CNZ)
  872. C
  873. IF (IDIM .EQ. 3) THEN
  874. FLUCEL = ( (IFLU.PROG(2)*CNZ) + (IFLU.PROG(3)*CTZ)
  875. & + (IFLU.PROG(4)*CVZ) ) * SURF
  876. MPORES.VPOCHA(NLCEG,5) = MPORES.VPOCHA(NLCEG,5) -
  877. & (FLUCEL / VOLG)
  878. IF (NLCED .NE. NLCEG)
  879. & MPORES.VPOCHA(NLCED,5) = MPORES.VPOCHA(NLCED,5) +
  880. & (FLUCEL / VOLD)
  881. ENDIF
  882. C
  883. C Rho Un et
  884. C
  885. FLUCEL = IFLU.PROG(2 + IDIM) * SURF
  886. MPORES.VPOCHA(NLCEG,3 + IDIM) =
  887. & MPORES.VPOCHA(NLCEG, 3 + IDIM) -
  888. & (FLUCEL / VOLG)
  889. IF (NLCED .NE. NLCEG)
  890. & MPORES.VPOCHA(NLCED, 3 + IDIM) =
  891. & MPORES.VPOCHA(NLCED, 3 + IDIM) +
  892. & (FLUCEL / VOLD)
  893. C
  894. ENDIF
  895. C
  896. C********** Riemann problem 2-2
  897. C
  898. IF (S2_22 .GE. EPS) THEN
  899. IWST_L.PROG(1) = ROG2
  900. IWST_L.PROG(2) = UNG2
  901. C IWST_L.PROG(2) = 0.0D0
  902. IWST_L.PROG(3) = UTG2
  903. IWST_L.PROG(IDIM + 2) = PG2
  904. IWST_L.PROG(IDIM + 3) = 1.0D0
  905. C
  906. IWST_R.PROG(1) = ROD2
  907. IWST_R.PROG(2) = UND2
  908. C IWST_R.PROG(2) = 0.0D0
  909. IWST_R.PROG(3) = UTD2
  910. IWST_R.PROG(IDIM + 2) = PD2
  911. IWST_R.PROG(IDIM + 3) = 1.0D0
  912. IF (IDIM .EQ. 3) THEN
  913. IWST_L.PROG(4) = UVG2
  914. IWST_R.PROG(4) = UVD2
  915. ENDIF
  916. C
  917. C do i2 = 1 , 4 + NESP + IDIM, 1
  918. C write(*,*) i2, iwst_l.prog(i2), iwst_r.prog(i2)
  919. C enddo
  920. C
  921. IF (INDMET .EQ. 1) THEN
  922. CALL FVLHRE(IDIM, NESP, NLHS, NORD,
  923. & MLRCHE.PROG, PROPHY.XTAB, RUNIV,
  924. & TMAX,
  925. & CNNA, CTNA, CVNA,
  926. & IWST_L.PROG,
  927. & IWST_R.PROG,
  928. & IWWORK.PROG, IY.PROG, ICVL.PROG, ICVR.PROG,
  929. & ICVRSS.PROG, ICVWOR.PROG,
  930. & IFLU.PROG, IFLULA.PROG, CELLT,
  931. & LOGAN)
  932. ELSEIF(INDMET .EQ. 2)THEN
  933. CALL FBECRE(IDIM, NESP, NLHS, NORD,
  934. & MLRCHE.PROG, PROPHY.XTAB, RUNIV,
  935. & TMAX,
  936. & CNNA, CTNA, CVNA,
  937. & IWST_L.PROG,
  938. & IWST_R.PROG,
  939. & IWWORK.PROG, IY.PROG, ICVL.PROG, ICVR.PROG,
  940. & ICVRSS.PROG, ICVWOR.PROG,
  941. & IFLU.PROG, IFLULA.PROG, CELLT,
  942. & LOGAN)
  943. ELSEIF(INDMET .EQ. 3)THEN
  944. C
  945. C********** AUSMP-up
  946. C
  947. V_INF=MAX(MPUIN2.VPOCHA(NLCEG,1),
  948. & MPUIN2.VPOCHA(NLCED,1))
  949. CALL FAUSRE(IDIM, NESP, NLHS, NORD,
  950. & MLRCHE.PROG, PROPHY.XTAB, RUNIV,
  951. & TMAX, V_INF,
  952. & CNNA, CTNA, CVNA,
  953. & IWST_L.PROG,
  954. & IWST_R.PROG,
  955. & IWWORK.PROG, IY.PROG, ICVL.PROG, ICVR.PROG,
  956. & ICVRSS.PROG, ICVWOR.PROG,
  957. & IFLU.PROG, IFLULA.PROG, CELLT,
  958. & LOGAN)
  959. ENDIF
  960. C
  961. C write(*,*) 'Riemann problem 2-2'
  962. C write(*,*) iflula.prog(idim+1), iflu.prog(idim+1)
  963. C
  964. C subroutine fcolre(ndim, nesp, nLHS, nordpo,
  965. C & reacoe, aprop, Runiv,
  966. C & Tmaxcv,
  967. C & wvect_l,
  968. C & wvect_r,
  969. C & wwork, y, acv_l, acv_r, acv_rss, acv_work,
  970. C & flu, flu_lag, cellt,
  971. C & logan)
  972. C
  973. C
  974. IF (LOGAN) THEN
  975. WRITE(IOIMP,*) 'SUBROUTINE KODFL1'
  976. WRITE(IOIMP,*) 'RIEMANN PROBLEM 2-2'
  977. WRITE(IOIMP,*) 'FACE ', NGCF, ' S22 ', S2_22
  978. WRITE(IOIMP,*) ' ICOMP WVECT_L WVECT_R'
  979. DO I2 = 1, 2 + IDIM
  980. WRITE(IOIMP,'(I6,2E16.6)') I2, IWST_L.PROG(I2),
  981. $ IWST_R.PROG(I2)
  982. ENDDO
  983. WRITE(IOIMP,*) 'ANOMALY DETECTED'
  984. GOTO 9999
  985. ENDIF
  986. C
  987. CELLTM = MAX(CELLT, CELLTM)
  988. C
  989. C
  990. C************* Ecriture du residu
  991. C
  992. C FLUX(7) = RO Un RO Un
  993. C FLUX(8) = RO Un Un + P -> RO Un Ux + P CNX
  994. C FLUX(9) = RO Un Ut -> RO Un Uy + P CNY
  995. C FLUX(10)= RO Un Et RO Un Et
  996. C
  997. SURF = SURF0 * S2_22
  998. C
  999. C Rho un
  1000. C
  1001. FLUCEL = IFLU.PROG(1) * SURF
  1002. MPORES.VPOCHA(NLCEG, (IDIM + 5)) =
  1003. & MPORES.VPOCHA(NLCEG, (IDIM + 5)) -
  1004. & (FLUCEL / VOLG)
  1005. IF (NLCED .NE. NLCEG)
  1006. & MPORES.VPOCHA(NLCED, (IDIM + 5)) =
  1007. & MPORES.VPOCHA(NLCED, (IDIM + 5)) +
  1008. & (FLUCEL / VOLD)
  1009. C
  1010. C Rho un ux + P cnx
  1011. C
  1012. FLUCEL = ((IFLU.PROG(2)*CNX) + (IFLU.PROG(3)*CTX) +
  1013. & (IFLU.PROG(4)*CVX)) * SURF
  1014. MPORES.VPOCHA(NLCEG, (IDIM + 6)) =
  1015. & MPORES.VPOCHA(NLCEG, (IDIM + 6)) -
  1016. & (FLUCEL / VOLG)
  1017. IF (NLCED .NE. NLCEG)
  1018. & MPORES.VPOCHA(NLCED, (IDIM + 6)) =
  1019. & MPORES.VPOCHA(NLCED, (IDIM + 6)) +
  1020. & (FLUCEL / VOLD)
  1021. C
  1022. C Rho un uy + P cny
  1023. C
  1024. FLUCEL = ((IFLU.PROG(2)*CNY) + (IFLU.PROG(3)*CTY) +
  1025. & (IFLU.PROG(4)*CVY)) * SURF
  1026. MPORES.VPOCHA(NLCEG, (IDIM + 7)) =
  1027. & MPORES.VPOCHA(NLCEG, (IDIM + 7)) -
  1028. & (FLUCEL / VOLG)
  1029. IF (NLCED .NE. NLCEG)
  1030. & MPORES.VPOCHA(NLCED, (IDIM + 7)) =
  1031. & MPORES.VPOCHA(NLCED, (IDIM + 7)) +
  1032. & (FLUCEL / VOLD)
  1033. C
  1034. C Rho un uz + P cnz
  1035. C
  1036. IF (IDIM .EQ. 3) THEN
  1037. FLUCEL = ((IFLU.PROG(2)*CNZ) + (IFLU.PROG(3)*CTZ) +
  1038. & (IFLU.PROG(4)*CVZ)) * SURF
  1039. MPORES.VPOCHA(NLCEG, (IDIM + 8)) =
  1040. & MPORES.VPOCHA(NLCEG, (IDIM + 8)) -
  1041. & (FLUCEL / VOLG)
  1042. IF (NLCED .NE. NLCEG)
  1043. & MPORES.VPOCHA(NLCED, (IDIM + 8)) =
  1044. & MPORES.VPOCHA(NLCED, (IDIM + 8)) +
  1045. & (FLUCEL / VOLD)
  1046. ENDIF
  1047. C
  1048. C Rho un et
  1049. C
  1050. FLUCEL = IFLU.PROG(IDIM + 2) * SURF
  1051. MPORES.VPOCHA(NLCEG, ((2 * IDIM) + 6)) =
  1052. & MPORES.VPOCHA(NLCEG, ((2 * IDIM) + 6)) -
  1053. & (FLUCEL / VOLG)
  1054. IF (NLCED .NE. NLCEG)
  1055. & MPORES.VPOCHA(NLCED, ((2 * IDIM) + 6)) =
  1056. & MPORES.VPOCHA(NLCED, ((2 * IDIM) + 6)) +
  1057. & (FLUCEL / VOLD)
  1058. ENDIF
  1059. C
  1060. C********** Riemann problem 1-2
  1061. C From a numerical point of view, S1_12 and S2_12 are identical.
  1062. C The volume fraction of 1 is larger in l than in r.
  1063. C
  1064. C rho, u, p for the phase 1 and 2.
  1065. C NB in case of wall-boundary conditions, we cannot be here
  1066. C since ALPD1 = ALPG1, ALPD2 = ALPG2.
  1067. C
  1068. IF ((S1_12 .gt. EPS) .or. (S2_12 .gt. EPS)) THEN
  1069. C
  1070. IWST_L.PROG(1) = ROG1
  1071. IWST_L.PROG(2) = UNG1
  1072. C IWST_L.PROG(2) = 0.0D0
  1073. IWST_L.PROG(3) = UTG1
  1074. IWST_L.PROG(IDIM + 2) = PG1
  1075. IWST_L.PROG(IDIM + 3) = 0.0D0
  1076. C
  1077. IWST_R.PROG(1) = ROD2
  1078. IWST_R.PROG(2) = UND2
  1079. C IWST_R.PROG(2) = 0.0D0
  1080. IWST_R.PROG(3) = UTD2
  1081. IWST_R.PROG(IDIM + 2) = PD2
  1082. IWST_R.PROG(IDIM + 3) = 1.0D0
  1083. IF (IDIM .EQ. 3) THEN
  1084. IWST_L.PROG(4) = UVG1
  1085. IWST_R.PROG(4) = UVD2
  1086. ENDIF
  1087. C
  1088. C do i2 = 1 , 4 + NESP + IDIM, 1
  1089. C write(*,*) i2, iwst_l.prog(i2), iwst_r.prog(i2)
  1090. C enddo
  1091. C
  1092. CALL FBECRE(IDIM, NESP, NLHS, NORD,
  1093. & MLRCHE.PROG, PROPHY.XTAB, RUNIV,
  1094. & TMAX,
  1095. & CNNA, CTNA, CVNA,
  1096. & IWST_L.PROG,
  1097. & IWST_R.PROG,
  1098. & IWWORK.PROG, IY.PROG, ICVL.PROG, ICVR.PROG,
  1099. & ICVRSS.PROG, ICVWOR.PROG,
  1100. & IFLU.PROG, IFLULA.PROG, CELLT,
  1101. & LOGAN)
  1102. C
  1103. C write(*,*) 'Riemann problem 1-2'
  1104. C write(*,*) iflula.prog(idim+1), iflu.prog(idim+1)
  1105. C
  1106. C subroutine fcolre(ndim, nesp, nLHS, nordpo,
  1107. C & reacoe, aprop, Runiv,
  1108. C & Tmaxcv,
  1109. C & wvect_l,
  1110. C & wvect_r,
  1111. C & wwork, y, acv_l, acv_r, acv_rss, acv_work,
  1112. C & flu, flu_lag, cellt,
  1113. C & logan)
  1114. C
  1115. IF (LOGAN) THEN
  1116. WRITE(IOIMP,*) 'SUBROUTINE KODFL1'
  1117. WRITE(IOIMP,*) 'RIEMANN PROBLEM 1-2'
  1118. WRITE(IOIMP,*) 'FACE ', NGCF, ' S12 ', S1_12, S2_12
  1119. WRITE(IOIMP,*) ' ICOMP WVECT_L WVECT_R'
  1120. DO I2 = 1, 2 + IDIM
  1121. WRITE(IOIMP,'(I6,2E16.6)') I2, IWST_L.PROG(I2),
  1122. $ IWST_R.PROG(I2)
  1123. ENDDO
  1124. WRITE(IOIMP,*) 'ANOMALY DETECTED'
  1125. GOTO 9999
  1126. ENDIF
  1127. C
  1128. CELLTM = MAX(CELLT, CELLTM)
  1129. C
  1130. C************* Ecriture du residu
  1131. C
  1132. C FLUX(7) = RO Un RO Un
  1133. C FLUX(8) = RO Un Un + P -> RO Un Ux + P CNX
  1134. C FLUX(9) = RO Un Ut -> RO Un Uy + P CNY
  1135. C FLUX(10)= RO Un Et RO Un Et
  1136. C
  1137. C FLULA(3+IDIM) = - D_resh, with D_resh is the reactive shock
  1138. C speed
  1139. C
  1140. C We update resu (we are dealing with a reactive
  1141. C Riemann problem)
  1142. C
  1143. IF (IFLULA.PROG(IDIM + 3) .LT. 0.0D0) THEN
  1144. C
  1145. C Note that IFLULA.PROG(IDIM + 3) = - D_resh
  1146. C
  1147. C t
  1148. C |F / Flag
  1149. C ---> --->
  1150. C | /
  1151. C 1 | / 2
  1152. C | /
  1153. C |/
  1154. C ------------------------------> x
  1155. C ifac ifac+1
  1156. C
  1157. C D_resh > 0
  1158. C IFLULA.PROG(5) < 0
  1159. C
  1160. C Definition of SCON (conservative surface) and SG and SD
  1161. C (lagrangian surfaces)
  1162. C
  1163. c NB SCON > 0 if, for a positive conservative flux, gives negative
  1164. C contribution in ifac and positive contribution
  1165. C in ifac+1
  1166. C If follows that SCON \ge 0
  1167. C SG and SD > 0 if, for a positive lagrangian flux, gives
  1168. C negative contribution in ifac and negative
  1169. C contributuion in ifac+1
  1170. C
  1171.  
  1172. C In first order case :
  1173. C PHASE 1, conservative contribution in ifac (-F) and ifac+1 (+F)
  1174. C non-conservative contribution in ifac+1 (-Flag)
  1175. C negative contribution => SG and SD > 0
  1176. C
  1177. SCON1 = S1_12E
  1178. SD1 = S1_12E + S1_12R
  1179. SG1 = S1_12L
  1180. C In first order case :
  1181. C PHASE 2, no conservative contribution
  1182. C non-conservative contribution in ifac+1 (Flag)
  1183. C positive contribution => SG and SD < 0
  1184. C
  1185. SCON2 = 0.0D0
  1186. SD2 = -1.0D0 * (S2_12E + S2_12R)
  1187. SG2 = -1.0D0 * S2_12L
  1188. ELSE
  1189. C
  1190. C t
  1191. C \Flag | F
  1192. C ---> --->
  1193. C \ |
  1194. C 1 \ | 2
  1195. C \ |
  1196. C \|
  1197. C ------------------------------> x
  1198. C ifac ifac+1
  1199. C
  1200. C D_resh < 0
  1201. C IFLULA.PROG(IDIM + 3) > 0
  1202. C
  1203. C In first order case :
  1204. C PHASE 1, no conservative contribution in ifac and ifac+1
  1205. C non-conservative contribution in ifac (-Flag)
  1206. C negative contribution => SG and SD > 0
  1207. C
  1208. SCON1 = 0.0D0
  1209. SG1 = S1_12E + S1_12L
  1210. SD1 = S1_12R
  1211. C In first order case :
  1212. C PHASE 2, conservative contribution in ifac and ifac+1
  1213. C non-conservative contribution in ifac (Flag)
  1214. C positive contribution => SG and SD < 0
  1215. C
  1216. SCON2 = S2_12E
  1217. SG2 = -1.0D0 * (S2_12L + S2_12E)
  1218. SD2 = -1.0D0 * S2_12R
  1219. ENDIF
  1220. SURF = SURF0
  1221. CALL RDEMUP(IDIM, 0,
  1222. & SCON1, SG1, SD1,
  1223. & SCON2, SG2, SD2,
  1224. & SURF,VOLG,VOLD,
  1225. & CNX,CNY,CNZ,CTX,CTY,CTZ,CVX,CVY,CVZ,
  1226. & IFLU.PROG,IFLULA.PROG,RESG,RESD)
  1227. C write(*,*) 'RESD'
  1228. DO ICOMP = 1, ((2 * IDIM) + 6)
  1229. MPORES.VPOCHA(NLCED,ICOMP) =
  1230. & MPORES.VPOCHA(NLCED,ICOMP) + RESD(ICOMP)
  1231. C write(*,*) ICOMP, RESD(ICOMP)
  1232. ENDDO
  1233. C write(*,*) 'RESG'
  1234. DO ICOMP = 1, ((2 * IDIM) + 6)
  1235. MPORES.VPOCHA(NLCEG,ICOMP) =
  1236. & MPORES.VPOCHA(NLCEG,ICOMP) + RESG(ICOMP)
  1237. C write(*,*) ICOMP, RESG(ICOMP)
  1238. ENDDO
  1239. ENDIF
  1240. C
  1241. C********** Riemann problem 2-1
  1242. C From a numerical point of view, S2_21 and S2_21 are identical.
  1243. C The volume fraction of 2 is larger in l than in r.
  1244. C
  1245. C rho, u, p for the phase 1 and 2.
  1246. C
  1247. C NB in case of wall-boundary conditions, we cannot be here
  1248. C since ALPD1 = ALPG1, ALPD2 = ALPG2.
  1249. C
  1250. IF ((S1_21 .gt. EPS) .or. (S2_21 .gt. EPS)) THEN
  1251. C
  1252. C
  1253. IWST_L.PROG(1) = ROG2
  1254. IWST_L.PROG(2) = UNG2
  1255. CC IWST_L.PROG(2) = 0.0D0
  1256. IWST_L.PROG(3) = UTG2
  1257. IWST_L.PROG(IDIM + 2) = PG2
  1258. IWST_L.PROG(IDIM + 3) = 1.0D0
  1259. CC
  1260. IWST_R.PROG(1) = ROD1
  1261. IWST_R.PROG(2) = UND1
  1262. C IWST_R.PROG(2) = 0.0D0
  1263. IWST_R.PROG(3) = UTD1
  1264. IWST_R.PROG(IDIM + 2) = PD1
  1265. IWST_R.PROG(IDIM + 3) = 0.0D0
  1266. C
  1267. IF (IDIM .EQ. 3) THEN
  1268. IWST_L.PROG(4) = UVG2
  1269. IWST_R.PROG(4) = UVD1
  1270. ENDIF
  1271. C
  1272. CC
  1273. C do i2 = 1 , 4 + NESP + IDIM, 1
  1274. C write(*,*) i2, iwst_l.prog(i2), iwst_r.prog(i2)
  1275. C enddo
  1276. C
  1277. CALL FBECRE(IDIM, NESP, NLHS, NORD,
  1278. & MLRCHE.PROG, PROPHY.XTAB, RUNIV,
  1279. & TMAX,
  1280. & CNNA, CTNA, CVNA,
  1281. & IWST_L.PROG,
  1282. & IWST_R.PROG,
  1283. & IWWORK.PROG, IY.PROG, ICVL.PROG, ICVR.PROG,
  1284. & ICVRSS.PROG, ICVWOR.PROG,
  1285. & IFLU.PROG, IFLULA.PROG, CELLT,
  1286. & LOGAN)
  1287. C write(*,*) 'Riemann problem 1-2'
  1288. C write(*,*) iflula.prog(idim+1), iflu.prog(idim+1)
  1289. C
  1290. C subroutine fcolre(ndim, nesp, nLHS, nordpo,
  1291. C & reacoe, aprop, Runiv,
  1292. C & Tmaxcv,
  1293. C & wvect_l,
  1294. C & wvect_r,
  1295. C & wwork, y, acv_l, acv_r, acv_rss, acv_work,
  1296. C & flu, flu_lag, cellt,
  1297. C & logan)
  1298. C
  1299. IF (LOGAN) THEN
  1300. WRITE(IOIMP,*) 'SUBROUTINE KODFL1'
  1301. WRITE(IOIMP,*) 'RIEMANN PROBLEM 2-1'
  1302. WRITE(IOIMP,*) 'FACE ', NGCF, ' S21 ', S1_21, S2_21
  1303. WRITE(IOIMP,*) ' ICOMP WVECT_L WVECT_R'
  1304. DO I2 = 1, 2 + IDIM
  1305. WRITE(IOIMP,'(I6,2E16.6)') I2, IWST_L.PROG(I2),
  1306. $ IWST_R.PROG(I2)
  1307. ENDDO
  1308. WRITE(IOIMP,*) 'ANOMALY DETECTED'
  1309. GOTO 9999
  1310. ENDIF
  1311. C
  1312. CELLTM = MAX(CELLT, CELLTM)
  1313. C
  1314. C************* Ecriture du residu
  1315. C
  1316. C FLUX(7) = RO Un RO Un
  1317. C FLUX(8) = RO Un Un + P -> RO Un Ux + P CNX
  1318. C FLUX(9) = RO Un Ut -> RO Un Uy + P CNY
  1319. C FLUX(10)= RO Un Et RO Un Et
  1320. C
  1321. C FLULA(3+IDIM) = - D_resh, with D_resh is the reactive shock
  1322. C speed
  1323. C
  1324. C
  1325. C We update resu (we are dealing with a reactive
  1326. C Riemann problem)
  1327. C
  1328. IF (IFLULA.PROG(IDIM + 3) .LT. 0.0D0) THEN
  1329. C
  1330. C Note that IFLULA.PROG(IDIM + 3) = - D_resh
  1331. C
  1332. C t
  1333. C |F / Flag
  1334. C ---> --->
  1335. C | /
  1336. C 2 | / 1
  1337. C | /
  1338. C |/
  1339. C ------------------------------> x
  1340. C ifac ifac+1
  1341. C
  1342. C D_resh > 0
  1343. C IFLULA.PROG(IDIM + 3) < 0
  1344. C
  1345. C In first order case :
  1346. C PHASE 2, conservative contribution in ifac (-F) and ifac+1 (+F)
  1347. C non-conservative contribution in ifac+1 (-Flag)
  1348. C negative contribution => SG and SD > 0
  1349. C
  1350. SCON2 = S2_21E
  1351. SG2 = S2_21L
  1352. SD2 = S2_21E + S2_21R
  1353. C
  1354. C PHASE 1, non-conservative contribution in ifac+1 (+Flag)
  1355. C positive contribution => SG and SD < 0
  1356. C
  1357. SCON1 = 0.0D0
  1358. SG1 = -1.0D0 * S1_21L
  1359. SD1 = -1.0D0 * (S1_21E + S1_21R)
  1360. C
  1361. ELSE
  1362. C
  1363. C t
  1364. C \Flag | F
  1365. C ---> --->
  1366. C \ |
  1367. C 2 \ | 1
  1368. C \ |
  1369. C \|
  1370. C ------------------------------> x
  1371. C ifac ifac+1
  1372. C
  1373. C D_resh < 0
  1374. C IFLULA.PROG(IDIM + 3) > 0
  1375. C
  1376. C In first order case :
  1377. C PHASE 2, no conservative contribution
  1378. C non-conservative contribution in ifac (-Flag)
  1379. C negative contribution => SG and SD > 0
  1380. C
  1381. SCON2 = 0.0D0
  1382. SG2 = S2_21E + S2_21L
  1383. SD2 = S2_21R
  1384. C
  1385. C In first order case :
  1386. C PHASE 1, conservative contribution in ifac (-F) and
  1387. C ifac+1 (+F)
  1388. C non-conservative contribution in ifac (+Flag)
  1389. C positive contribution => SG and SD > 0
  1390. C
  1391. C ALPHA1
  1392. C
  1393. SCON1 = S1_21E
  1394. SG1 = -1.0D0 * (S1_21E + S1_21L)
  1395. SD1 = -1.0D0 * S1_21R
  1396. ENDIF
  1397. SURF = SURF0
  1398. CALL RDEMUP(IDIM, 0,
  1399. & SCON1, SG1, SD1,
  1400. & SCON2, SG2, SD2,
  1401. & SURF,VOLG,VOLD,
  1402. & CNX,CNY,CNZ,CTX,CTY,CTZ,CVX,CVY,CVZ,
  1403. & IFLU.PROG,IFLULA.PROG,RESG,RESD)
  1404. C write(*,*) 'RESD'
  1405. DO ICOMP = 1, ((2 * IDIM) + 6)
  1406. MPORES.VPOCHA(NLCED,ICOMP) =
  1407. & MPORES.VPOCHA(NLCED,ICOMP) + RESD(ICOMP)
  1408. C write(*,*) ICOMP, RESD(ICOMP)
  1409. ENDDO
  1410. C write(*,*) 'RESG'
  1411. DO ICOMP = 1, ((2 * IDIM) + 6)
  1412. MPORES.VPOCHA(NLCEG,ICOMP) =
  1413. & MPORES.VPOCHA(NLCEG,ICOMP) + RESG(ICOMP)
  1414. C write(*,*) ICOMP, RESG(ICOMP)
  1415. ENDDO
  1416. ENDIF
  1417. C
  1418. C********** Calcul du pas du temps (CFL)
  1419. C
  1420. C********* a) etat a l'interface
  1421. C
  1422. DIAMG = MPOVDI.VPOCHA(NLCEG,1)
  1423. DIAMD = MPOVDI.VPOCHA(NLCED,1)
  1424. DIAM = MIN(DIAMG,DIAMD)
  1425. CELL = CELLTM/DIAM
  1426. IF(CELL .GT. UNSDT)THEN
  1427. UNSDT = CELL
  1428. ENDIF
  1429. ENDIF
  1430. C
  1431. C**** Fin boucle sur FACEL
  1432. C
  1433. ENDDO
  1434. C
  1435. C**** Pas du temps (condition de non interaction en 1D)
  1436. C
  1437. DT = 0.5D0 / UNSDT
  1438. C
  1439. C**** Desactivation des segments
  1440. C
  1441. SEGSUP MLELIM
  1442. SEGSUP MLENT1
  1443. SEGSUP MLENT2
  1444. C
  1445. SEGDES MPOVSU
  1446. SEGDES MPOVDI
  1447. SEGDES MPOVVO
  1448. SEGDES MPORES
  1449. SEGDES IPT2
  1450. C
  1451. C Warning. MLRMFR = PGAS . 'MASSFRA'
  1452. C MLRCHE = PGAS . 'CHEMCOEF'
  1453. C
  1454. SEGDES MLRMFR
  1455. SEGDES MLRCHE
  1456. SEGSUP IY
  1457. C
  1458. SEGSUP IWWORK
  1459. SEGSUP IWST_R
  1460. SEGSUP IWST_L
  1461. SEGSUP ICVL
  1462. SEGSUP ICVR
  1463. SEGSUP ICVRSS
  1464. SEGSUP ICVWOR
  1465. SEGSUP IFLU
  1466. SEGSUP IFLULA
  1467. C
  1468. SEGDES MPALP1
  1469. SEGDES MPALP2
  1470. SEGDES MELAL1
  1471. SEGDES MELRO1
  1472. SEGDES MELP1
  1473. SEGDES MELAL2
  1474. SEGDES MELRO2
  1475. SEGDES MELP2
  1476. SEGDES MELUN2
  1477. SEGDES MELUT2
  1478. IF (IDIM .EQ. 2) THEN
  1479. SEGDES MELUN1
  1480. SEGDES MELUT1
  1481. SEGDES MELVNX
  1482. SEGDES MELVNY
  1483. SEGDES MELT1X
  1484. SEGDES MELT1Y
  1485. SEGDES MELUN2
  1486. SEGDES MELUT2
  1487. ELSE
  1488. SEGDES MELVNX
  1489. SEGDES MELVNY
  1490. SEGDES MELVNZ
  1491. SEGDES MELT1X
  1492. SEGDES MELT1Y
  1493. SEGDES MELT1Z
  1494. SEGDES MELT2X
  1495. SEGDES MELT2Y
  1496. SEGDES MELT2Z
  1497. SEGDES MELUN1
  1498. SEGDES MELUT1
  1499. SEGDES MELUV1
  1500. SEGDES MELUN2
  1501. SEGDES MELUT2
  1502. SEGDES MELUV2
  1503. ENDIF
  1504. C
  1505. SEGDES MPOALP
  1506. IF (INDK0 .EQ. 2) SEGDES MPOK0
  1507. C
  1508. IF(MPUIN2 .NE. 0)THEN
  1509. SEGDES MPUIN2
  1510. SEGDES MPUIN1
  1511. ENDIF
  1512. C
  1513. 9999 CONTINUE
  1514. RETURN
  1515. END
  1516. CC
  1517.  
  1518.  
  1519.  
  1520.  
  1521.  

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