Télécharger ckon4.eso

Retour à la liste

Numérotation des lignes :

  1. C CKON4 SOURCE PV 09/03/12 21:17:08 6325
  2. SUBROUTINE CKON4(LOGME,INDMET,NORDP1,
  3. & IROF,IVITF,IPF,IFRMAF,ISCAF,PROPHY,
  4. & ICHPSU,ICHPDI,
  5. & MELEMC,MELEMF,MELEFE,
  6. & IZG1,IZG2,IZG3,IZG4,IZG5,DT,DIAMEL,NLCEMI,
  7. & LOGNC,LOGAN,MESERR)
  8. C************************************************************************
  9. C
  10. C PROJET : CASTEM 2000
  11. C
  12. C NOM : CKON4
  13. C
  14. C DESCRIPTION : Voir CKON
  15. C
  16. C Cas trois dimensions, gaz "thermally perfect"
  17. C
  18. C LANGAGE : FORTRAN 77 + ESOPE 2000 (avec estensions CISI)
  19. C
  20. C AUTEUR : R. MOREL, A. BECCANTINI, DRN/DMT/SEMT/TTMF
  21. C
  22. C************************************************************************
  23. C
  24. C
  25. C APPELES (Outils
  26. C CASTEM) : KRIPAD, LICHT
  27. C
  28. C APPELES (Calcul) : FVLHT2, FCOLT2
  29. C
  30. C
  31. C************************************************************************
  32. C
  33. C ENTREES
  34. C
  35. C
  36. C 1) PARAMETRES
  37. C
  38. C LOGME : (LOGICAL); .TRUE. -> MULTI-ESPECES
  39. C .FALSE. -> MONO-ESPECE
  40. C
  41. C INDMET : 1 VLH (van Leer Hanel FVS)
  42. C
  43. C 2 HUSVLH (VLH FVS + Colella-Glaz FDS)
  44. C
  45. C NORDP1 : (ordre des polynoms cv(T)) + 1
  46. C
  47. C 2) Pointeurs des MCHAMLs
  48. C
  49. C IROF : MCHAML sur "FACEL" contenant la masse volumique
  50. C ("gauche" et "droite");
  51. C
  52. C IVITF : MCHAML sur "FACEL" contenant la vitesse dans le repaire
  53. C local (n,t) et les cosinus directeurs des repaire local;
  54. C
  55. C IPF : MCHAML sur "FACEL" contenant la pression;
  56. C
  57. C IFRAMAF : MCHAML sur "FACEL", contenant les fractions massiques
  58. C si LOGME = .TRUE.;
  59. C LOGME = .FALSE. -> IFRAMAF = 0
  60. C
  61. C ISCAF : MCHAML sur "FACEL", contenant les scalires passifs a
  62. C transporter (ou ISCAF = 0)
  63. C
  64. C 3) Pointeur sur la table des proprietes de gaz
  65. C
  66. C PROPHY
  67. C
  68. C 4) Pointeurs de CHPOINTs de la table DOMAINE
  69. C
  70. C ICHPSU : CHPOINT "FACE" contenant la surface des faces
  71. C
  72. C ICHPDI : CHPOINT "CENTRE" contenant le diametre minimum
  73. C de chaque element
  74. C
  75. C
  76. C 5) Pointeurs de MELEME de la table DOMAINE
  77. C
  78. C MELEMC : MELEME 'CENTRE' du SPG des CENTRES
  79. C
  80. C MELEMF : MELEME 'FACE' du SPG des FACES
  81. C
  82. C MELEFE : MELEME 'FACEL' du connectivité FACES -> ELEM
  83. C
  84. C SORTIES (il faudrait dire E/S)
  85. C
  86. C IZGi : pointeurs de CHPOINTs "FACE" dont les valeurs
  87. C se trouvent dans la table KIZX . 'EQEX' . 'KIZG'
  88. C aux indices 'IZGi'
  89. C
  90. C IZG1 : Increment de masse voluique
  91. C
  92. C IZG2 : Increment de quantite de mouvement
  93. C
  94. C IZG3 : Increment de l'energie totale
  95. C
  96. C IZG4 : Increment de les Masse Volumiques des Especies
  97. C (si LOGME = .TRUE.)
  98. C
  99. C IZG5 : Increment de scalaires passifs
  100. C
  101. C DIAMEL : 'minimum' diametre du maillage
  102. C
  103. C NLCEMI : numero local du CENTRE ou le diametre est 'minimum'
  104. C
  105. C DT : pas de temps pour le respect de la CFL-like condition
  106. C DT < DIAMEL /2 /max(Lambda_i)
  107. C En maillage regulier cette condition garantie la
  108. C non-interaction des ondes
  109. C
  110. C
  111. C LOGNC : (LOGICAL): si .TRUE. la methode de Newton-Rapson, utilisée
  112. C dans pour la solution du probleme Riemann exacte ou dans
  113. C l'algorithm HUS, n'a pas bien marchéee; MESERR = 'Goudunov'
  114. C ou 'HUS'.
  115. C
  116. C LOGAN : (LOGICAL): si .TRUE. une anomalie à été detectée
  117. C
  118. C MESERR : pour l'ecriture des messages d'erreurs
  119. C
  120. C************************************************************************
  121. C
  122. C HISTORIQUE (Anomalies et modifications éventuelles)
  123. C
  124. C HISTORIQUE : 22.12.98 Creation
  125. C
  126. C 08.02.00 Ajout du flux numerique de Colella-Glaz
  127. C
  128. C 21.02.00 Ajout de transport de scalaires passifs
  129. C
  130. C************************************************************************
  131. C
  132. C
  133. C N.B.: On suppose qu'on a déjà controllé RO, P > 0
  134. C GAMMA \in (1,3)
  135. C Y \in (0,1)
  136. C Si non il faut le faire!!!
  137. C
  138. C************************************************************************
  139. C
  140. C
  141. C**** Variables de COOPTIO
  142. C
  143. C INTEGER IPLLB, IERPER, IERMAX, IERR, INTERR
  144. C & ,IOTER, IOLEC, IOIMP, IOCAR, IOACQ
  145. C & ,IOPER, IOSGB, IOGRA, IOSAU, IORES
  146. C & ,IECHO, IIMPI, IOSPI
  147. C & ,IDIM
  148. CC & ,MCOORD
  149. C & ,IFOMOD, NIFOUR, IFOUR, NSDPGE, IONIVE
  150. C & ,NGMAXY, IZROSF, ISOTYP, IOSCR,LTEXLU
  151. C & ,NORINC,NORVAL,NORIND,NORVAD
  152. C & ,NUCROU, IPSAUV
  153. C
  154. IMPLICIT INTEGER(I-N)
  155. INTEGER I1, I2
  156. & ,INDMET, NORDP1
  157. & ,IROF,IVITF,IPF,IFRMAF, ISCAF
  158. & ,ICHPSU,ICHPDI,MELEMC,MELEMF,MELEFE
  159. & ,IGEOMC,IGEOMF
  160. & ,IZG1,IZG2,IZG3,IZG4,IZG5,NLCEMI
  161. & ,NESP, NFAC, NSCA, NMA, NMA2
  162. & ,NLCF, NGCEG, NGCED, NLCEG, NLCED
  163. & ,NGCF, NLCF1, SPG1, SPG2
  164. REAL*8 DIAMEL, DT, UNSDT, CELLT
  165. & , ROG, UNG, UTG, UVG, PG, GAMG, TG
  166. & , ROD, UND, UTD, UVD, PD, GAMD, TD
  167. & , SURF,CNX, CNY, CNZ, CTX , CTY, CTZ
  168. & , CVX, CVY, CVZ
  169. & , CELL, DIAMG, DIAMD, DIAM
  170. & , ASON, LAMBDA
  171. & , Y1G, Y1D, RTOTG, RTOTD, CVTOTG, CVTOTD, ETHERG, ETHERD
  172. & , YG, YD, FUNTG, FUNTD, DCV, XST
  173. LOGICAL LOGME, LOGNC, LOGAN
  174. CHARACTER*(40) MESERR
  175. CHARACTER*(8) TYPE
  176. C
  177. C**** Segment des proprietes du gaz
  178. C
  179. SEGMENT PROPHY
  180. REAL*8 ACV(NORDP1,NESP+1), R(NESP+1), H0K(NESP+1)
  181. & ,ACVTOG(NORDP1), ACVTOD(NORDP1)
  182. ENDSEGMENT
  183. C
  184. C**** Les Includes
  185. C
  186. C
  187. C**** LES INCLUDES
  188. C
  189. -INC CCOPTIO
  190. -INC SMCOORD
  191. -INC SMCHAML
  192. POINTEUR MELVNX.MELVAL, MELVNY.MELVAL,MELVNZ.MELVAL,
  193. & MELT1X.MELVAL, MELT1Y.MELVAL,MELT1Z.MELVAL,
  194. & MELT2X.MELVAL, MELT2Y.MELVAL,MELT2Z.MELVAL
  195. POINTEUR MELVUN.MELVAL, MELVUT.MELVAL, MELVUV.MELVAL
  196. POINTEUR MELRO.MELVAL, MELP.MELVAL
  197. -INC SMCHPOI
  198. POINTEUR MPOVSU.MPOVAL, MPOVDI.MPOVAL
  199. & , MPOVG1.MPOVAL, MPOVG2.MPOVAL
  200. & , MPOVG3.MPOVAL, MPOVG4.MPOVAL, MPOVG5.MPOVAL
  201. POINTEUR MCHAMY.MCHAML, MCHAMS.MCHAML
  202. -INC SMELEME
  203. -INC SMLMOTS
  204. -INC SMLENTI
  205. C
  206. C**** Les fractiones massiques.
  207. C
  208. SEGMENT FRAMAS
  209. REAL*8 YET(NMA)
  210. ENDSEGMENT
  211. POINTEUR FRAMAG.FRAMAS, FRAMAD.FRAMAS, SCAG.FRAMAS, SCAD.FRAMAS
  212. C
  213. C**** Les flux aux interface dans le repaire (n,t)
  214. C
  215. SEGMENT IFLUX
  216. REAL*8 FLUX(NMA2)
  217. ENDSEGMENT
  218. POINTEUR IFLU1.IFLUX, IFLU2.IFLUX
  219. C
  220. C
  221. C**** Initialisation des MCHAMLs
  222. C
  223. C**** Masse volumique
  224. C
  225. MCHEL1 = IROF
  226. SEGACT MCHEL1
  227. MCHAM1 = MCHEL1.ICHAML(1)
  228. SEGACT MCHAM1
  229. MELRO = MCHAM1.IELVAL(1)
  230. SEGDES MCHEL1
  231. SEGDES MCHAM1
  232. C
  233. C**** Pression
  234. C
  235. MCHEL1 = IPF
  236. SEGACT MCHEL1
  237. MCHAM1 = MCHEL1.ICHAML(1)
  238. SEGACT MCHAM1
  239. MELP = MCHAM1.IELVAL(1)
  240. SEGDES MCHEL1
  241. SEGDES MCHAM1
  242. C
  243. C**** Vitesse et cosinus directeurs du repere (n,t)
  244. C
  245. MCHEL1 = IVITF
  246. SEGACT MCHEL1
  247. C
  248. C**** La vitesse a comme SPG MELEFE
  249. C Le cosinus directeurs ont comme SPG MELEMF
  250. C
  251. C MCHAM1 -> Cosinus directeurs
  252. C MCHAM2 -> Vitesse
  253. C
  254. SPG1 = MCHEL1.IMACHE(1)
  255. SPG2 = MCHEL1.IMACHE(2)
  256. IF((SPG1 .EQ. MELEMF) .AND. (SPG2 .EQ. MELEFE))THEN
  257. MCHAM1 = MCHEL1.ICHAML(1)
  258. MCHAM2 = MCHEL1.ICHAML(2)
  259. ELSEIF((SPG1 .EQ. MELEFE) .AND. (SPG2 .EQ. MELEMF))THEN
  260. MCHAM1 = MCHEL1.ICHAML(2)
  261. MCHAM2 = MCHEL1.ICHAML(1)
  262. ELSE
  263. LOGAN = .TRUE.
  264. GOTO 9999
  265. ENDIF
  266. SEGACT MCHAM1
  267. MELVNX = MCHAM1.IELVAL(1)
  268. MELVNY = MCHAM1.IELVAL(2)
  269. MELVNZ = MCHAM1.IELVAL(3)
  270. MELT1X = MCHAM1.IELVAL(4)
  271. MELT1Y = MCHAM1.IELVAL(5)
  272. MELT1Z = MCHAM1.IELVAL(6)
  273. MELT2X = MCHAM1.IELVAL(7)
  274. MELT2Y = MCHAM1.IELVAL(8)
  275. MELT2Z = MCHAM1.IELVAL(9)
  276. SEGDES MCHAM1
  277. SEGACT MCHAM2
  278. MELVUN = MCHAM2.IELVAL(1)
  279. MELVUT = MCHAM2.IELVAL(2)
  280. MELVUV = MCHAM2.IELVAL(3)
  281. SEGDES MCHAM2
  282. SEGDES MCHEL1
  283. C
  284. C**** Fractions massiques
  285. C
  286. IF(LOGME)THEN
  287. MCHEL1 = IFRMAF
  288. SEGACT MCHEL1
  289. MCHAMY = MCHEL1.ICHAML(1)
  290. SEGACT MCHAMY
  291. C
  292. C******* Numero d'especes dans les equations d'Euler
  293. C
  294. NESP = MCHAMY.IELVAL(/1)
  295. DO I1 = 1, NESP
  296. MELVA1 = MCHAMY.IELVAL(I1)
  297. SEGACT MELVA1
  298. ENDDO
  299. NMA = NESP
  300. SEGINI FRAMAG
  301. SEGINI FRAMAD
  302. SEGDES MCHEL1
  303. ELSE
  304. C
  305. C******* Definition minimale de YET, necessaire pour transmetre YET aux
  306. C subroutines FORTRAN qui calculent les flux
  307. C
  308. NMA = 1
  309. SEGINI FRAMAG
  310. SEGINI FRAMAD
  311. NESP = 0
  312. ENDIF
  313. C
  314. C**** Scalaires passifs
  315. C
  316. IF(ISCAF .GT. 0)THEN
  317. MCHEL1 = ISCAF
  318. SEGACT MCHEL1
  319. MCHAMS = MCHEL1.ICHAML(1)
  320. SEGACT MCHAMS
  321. NSCA = MCHAMS.IELVAL(/1)
  322. DO I1 = 1, NSCA, 1
  323. MELVA1 = MCHAMS.IELVAL(I1)
  324. SEGACT MELVA1
  325. ENDDO
  326. NMA = NSCA
  327. SEGINI SCAG
  328. SEGINI SCAD
  329. SEGDES MCHEL1
  330. ELSE
  331. C
  332. C******* Definition minimale de YET, necessaire pour transmetre YET aux
  333. C subroutines FORTRAN qui calculent les flux
  334. C
  335. NMA = 1
  336. SEGINI SCAG
  337. SEGINI SCAD
  338. NSCA= 0
  339. ENDIF
  340. C
  341. C
  342. C**** Initialisation des MELEMEs
  343. C
  344. C 'CENTRE', 'FACEL'
  345. C
  346. IPT2 = MELEFE
  347. SEGACT IPT2
  348. NFAC = IPT2.NUM(/2)
  349. C
  350. C**** KRIPAD pour la correspondance global/local de centre
  351. C
  352. CALL KRIPAD(MELEMC,MLENT1)
  353. C
  354. C**** MLENTI1 a MCORD.XCOORD(/1)/(IDIM+1) elements
  355. C
  356. C Si i est le numero global d'un noeud de ICEN,
  357. C MLENT1.LECT(i) contient sa position, i.e.
  358. C
  359. C I = numero global du noeud centre
  360. C MLENT1.LECT(i) = numero local du noeud centre
  361. C
  362. C MLENT1 déjà activé, i.e.
  363. C
  364. C SEGACT MLENT1
  365. C
  366. C
  367. C**** KRIPAD pour la correspondance global/local de 'FACE'
  368. C
  369. CALL KRIPAD(MELEMF,MLENT2)
  370. C
  371. C**** Initialisation de flux
  372. C
  373. NMA2 = 5 + NESP + NSCA
  374. SEGINI IFLU1
  375. SEGINI IFLU2
  376. C
  377. C**** IFLU2 = segment de travail en FLUVLH; c'est plus rapide le definir ici
  378. C
  379. C
  380. C**** CHPOINTs de la table DOMAINE
  381. C
  382. CALL LICHT(ICHPSU,MPOVSU,TYPE,IGEOMF)
  383. CALL LICHT(ICHPDI,MPOVDI,TYPE,IGEOMC)
  384. C
  385. C**** LICHT active les MPOVALs en *MOD
  386. C
  387. C i.e.
  388. C
  389. C SEGACT MPOVSU*MOD
  390. C SEGACT MPOVDI*MOD
  391. C
  392. C
  393. C**** Les FLUX aux face
  394. C
  395. C La densité
  396. C
  397. CALL LICHT(IZG1,MPOVG1,TYPE,IGEOMF)
  398. C
  399. C SEGACT MPOVG1*MOD
  400. C
  401. C**** Les debits
  402. C
  403. CALL LICHT(IZG2,MPOVG2,TYPE,IGEOMF)
  404. C
  405. C SEGACT MPOVG2*MOD
  406. C
  407. C**** L'energie totale volumique
  408. C
  409. CALL LICHT(IZG3,MPOVG3,TYPE,IGEOMF)
  410. C
  411. C SEGACT MPOVG3*MOD
  412. C
  413. C**** Les Fractions Massiques
  414. C
  415. IF(LOGME)THEN
  416. CALL LICHT(IZG4,MPOVG4,TYPE,IGEOMF)
  417. C
  418. C SEGACT MPOVG4*MOD
  419. C
  420. ENDIF
  421. C
  422. C**** Les scalaires passifs
  423. C
  424. IF(ISCAF .GT. 0)THEN
  425. CALL LICHT(IZG5,MPOVG5,TYPE,IGEOMF)
  426. C
  427. C SEGACT MPOVG5*MOD
  428. C
  429. ENDIF
  430. C
  431. C**** Activation des MCHAMLs
  432. C
  433. SEGACT MELRO
  434. SEGACT MELP
  435. SEGACT MELVUN
  436. SEGACT MELVUT
  437. SEGACT MELVUV
  438. SEGACT MELVNX
  439. SEGACT MELVNY
  440. SEGACT MELVNZ
  441. SEGACT MELT1X
  442. SEGACT MELT1Y
  443. SEGACT MELT1Z
  444. SEGACT MELT2X
  445. SEGACT MELT2Y
  446. SEGACT MELT2Z
  447. C
  448. C**** Initialisation de 1/DT
  449. C
  450. UNSDT = 0.0D0
  451. C
  452. C**** BOUCLE SUR FACEL pour le calcul du FLUX
  453. C
  454. DO NLCF = 1, NFAC
  455. C
  456. C******* NLCF = numero local du centre de facel
  457. C NGCF = numero global du centre de facel
  458. C NLCF1 = numero local du centre de face
  459. C NGCEG = numero global du centre ELT "gauche"
  460. C NLCEG = numero local du centre ELT "gauche"
  461. C NGCED = numero global du centre ELT "droite"
  462. C NLCED = numero local du centre ELT "droite"
  463. C
  464. NGCEG = IPT2.NUM(1,NLCF)
  465. NGCED = IPT2.NUM(3,NLCF)
  466. NGCF = IPT2.NUM(2,NLCF)
  467. NLCF1 = MLENT2.LECT(NGCF)
  468. NLCEG = MLENT1.LECT(NGCEG)
  469. NLCED = MLENT1.LECT(NGCED)
  470. C
  471. C******* NLCF != NLCF1 -> l'auteur (MOI) n'a rien compris.
  472. C
  473. IF(NLCF .NE. NLCF1)THEN
  474. MESERR = 'Il ne faut pas jouer avec la console. '
  475. LOGAN = .TRUE.
  476. GOTO 9999
  477. ENDIF
  478. C
  479. C******* Recuperation des Etats "gauche" et "droite"
  480. C
  481. ROG = MELRO.VELCHE(1,NLCF)
  482. UNG = MELVUN.VELCHE(1,NLCF)
  483. UTG = MELVUT.VELCHE(1,NLCF)
  484. UVG = MELVUV.VELCHE(1,NLCF)
  485. PG = MELP.VELCHE(1,NLCF)
  486. C
  487. ROD = MELRO.VELCHE(3,NLCF)
  488. UND = MELVUN.VELCHE(3,NLCF)
  489. UTD = MELVUT.VELCHE(3,NLCF)
  490. UVD = MELVUV.VELCHE(3,NLCF)
  491. PD = MELP.VELCHE(3,NLCF)
  492. C
  493. CNX = MELVNX.VELCHE(1,NLCF)
  494. CNY = MELVNY.VELCHE(1,NLCF)
  495. CNZ = MELVNZ.VELCHE(1,NLCF)
  496. CTX = MELT1X.VELCHE(1,NLCF)
  497. CTY = MELT1Y.VELCHE(1,NLCF)
  498. CTZ = MELT1Z.VELCHE(1,NLCF)
  499. CVX = MELT2X.VELCHE(1,NLCF)
  500. CVY = MELT2Y.VELCHE(1,NLCF)
  501. CVZ = MELT2Z.VELCHE(1,NLCF)
  502. C
  503. C******* Le fractiones massiques, R et ACVTO
  504. C
  505. RTOTG = 0.0D0
  506. RTOTD = 0.0D0
  507. CVTOTG = 0.0D0
  508. CVTOTD = 0.0D0
  509. Y1G = 1.0D0
  510. Y1D = 1.0D0
  511. ETHERG = 0.0D0
  512. ETHERD = 0.0D0
  513. DO I1= 1, NORDP1, 1
  514. PROPHY.ACVTOG(I1) =0.0D0
  515. PROPHY.ACVTOD(I1) =0.0D0
  516. ENDDO
  517. DO I1 = 1, NESP, 1
  518. MELVA1 = MCHAMY.IELVAL(I1)
  519. YG = MELVA1.VELCHE(1,NLCF)
  520. YD = MELVA1.VELCHE(3,NLCF)
  521. Y1G = Y1G - YG
  522. Y1D = Y1D - YD
  523. FRAMAG.YET(I1) = YG
  524. FRAMAD.YET(I1) = YD
  525. RTOTG = RTOTG + YG * PROPHY.R(I1)
  526. RTOTD = RTOTD + YD * PROPHY.R(I1)
  527. DO I2 = 1, NORDP1
  528. PROPHY.ACVTOG(I2) = PROPHY.ACVTOG(I2) +
  529. & YG * PROPHY.ACV(I2,I1)
  530. PROPHY.ACVTOD(I2) = PROPHY.ACVTOD(I2) +
  531. & YD * PROPHY.ACV(I2,I1)
  532. ENDDO
  533. ENDDO
  534. RTOTG = RTOTG + Y1G * PROPHY.R(NESP+1)
  535. RTOTD = RTOTD + Y1D * PROPHY.R(NESP+1)
  536. DO I2 = 1, NORDP1, 1
  537. PROPHY.ACVTOG(I2) = PROPHY.ACVTOG(I2) +
  538. & Y1G * PROPHY.ACV(I2,NESP+1)
  539. PROPHY.ACVTOD(I2) = PROPHY.ACVTOD(I2) +
  540. & Y1D * PROPHY.ACV(I2,NESP+1)
  541. ENDDO
  542. TG = PG / (ROG * RTOTG)
  543. TD = PD / (ROD * RTOTD)
  544. FUNTG = 1.0D0
  545. FUNTD = 1.0D0
  546. CVTOTG = PROPHY.ACVTOG(1)
  547. ETHERG = CVTOTG * TG
  548. CVTOTD = PROPHY.ACVTOD(1)
  549. ETHERD = CVTOTD * TD
  550. DO I1 = 2, NORDP1, 1
  551. FUNTG = FUNTG * TG
  552. FUNTD = FUNTD * TD
  553. DCV = PROPHY.ACVTOG(I1) * FUNTG
  554. CVTOTG = CVTOTG + DCV
  555. ETHERG = ETHERG + DCV * TG / I1
  556. DCV = PROPHY.ACVTOD(I1) * FUNTD
  557. CVTOTD = CVTOTD + DCV
  558. ETHERD = ETHERD + DCV * TD / I1
  559. ENDDO
  560. GAMG = (CVTOTG + RTOTG) /CVTOTG
  561. GAMD = (CVTOTD + RTOTD) /CVTOTD
  562. C
  563. C******* Les scalaires passifs
  564. C
  565. DO I1 = 1, NSCA, 1
  566. MELVA1 = MCHAMS.IELVAL(I1)
  567. SCAG.YET(I1) = MELVA1.VELCHE(1,NLCF)
  568. SCAD.YET(I1) = MELVA1.VELCHE(3,NLCF)
  569. ENDDO
  570. C
  571. C******* On a defini (ROg,ROUNg,ROUTg,Pg,(Yg)), (ROd,ROUNd,ROUTd,Pd,(Yd))
  572. C et on a déjà verifié ROg, ROd, Pg, Pd > 0 et 0<Y_i<1
  573. C
  574. C
  575. C******* Calcul du flux aux interfaces
  576. C
  577. IF(INDMET .EQ. 1)THEN
  578. C
  579. C********** VLH FVS
  580. C
  581. CALL FVLHT2(NESP,NSCA,
  582. & GAMG,ROG,PG,UNG,UTG,UVG,ETHERG,
  583. & GAMD,ROD,PD,UND,UTD,UVD,ETHERD,
  584. & FRAMAG.YET,FRAMAD.YET,
  585. & SCAG.YET,SCAD.YET,
  586. & IFLU1.FLUX,IFLU2.FLUX,
  587. & CELLT)
  588. ELSEIF(INDMET .EQ. 2)THEN
  589. C
  590. C******* Colella-Glaz FDS (avec Entropy-fix)
  591. C
  592. XST=0.0D0
  593. CALL FCOLT2(NESP,NSCA,NORDP1,XST,
  594. & GAMG,RTOTG,PROPHY.ACVTOG,ROG,PG,TG,UNG,UTG,UVG,ETHERG,
  595. & GAMD,RTOTD,PROPHY.ACVTOD,ROD,PD,TD,UND,UTD,UVD,ETHERD,
  596. & FRAMAG.YET,FRAMAD.YET,
  597. & SCAG.YET,SCAD.YET,
  598. & IFLU1.FLUX,CELLT,
  599. & LOGNC,MESERR,LOGAN)
  600. ENDIF
  601. C
  602. IF(LOGAN) GOTO 9999
  603. IF(LOGNC) GOTO 9999
  604. C
  605. C******* Ecriture des flux
  606. C
  607. C FLUX(1) = RO Un RO Un
  608. C FLUX(2) = RO Un Un + P -> RO Un Ux + P CNX
  609. C FLUX(3) = RO Un Ut -> RO Un Uy + P CNY
  610. C FLUX(4) = RO Un Et RO Un Et
  611. C
  612. SURF = MPOVSU.VPOCHA(NLCF,1)
  613. MPOVG1.VPOCHA(NLCF,1) = MPOVG1.VPOCHA(NLCF,1) +
  614. & (IFLU1.FLUX(1) * SURF )
  615. MPOVG2.VPOCHA(NLCF,1) = MPOVG2.VPOCHA(NLCF,1) +
  616. &((IFLU1.FLUX(2)*CNX+IFLU1.FLUX(3)*CTX+IFLU1.FLUX(4)*CVX) * SURF)
  617. MPOVG2.VPOCHA(NLCF,2) = MPOVG2.VPOCHA(NLCF,2) +
  618. &((IFLU1.FLUX(2)*CNY+IFLU1.FLUX(3)*CTY+IFLU1.FLUX(4)*CVY) * SURF)
  619. MPOVG2.VPOCHA(NLCF,3) = MPOVG2.VPOCHA(NLCF,3) +
  620. &((IFLU1.FLUX(2)*CNZ+IFLU1.FLUX(3)*CTZ+IFLU1.FLUX(4)*CVZ) * SURF)
  621. MPOVG3.VPOCHA(NLCF,1) = MPOVG3.VPOCHA(NLCF,1) +
  622. & (IFLU1.FLUX(5) * SURF)
  623. DO I1 = 1, NESP, 1
  624. MPOVG4.VPOCHA(NLCF,I1)=IFLU1.FLUX(5+I1)
  625. & * SURF
  626. ENDDO
  627. DO I1 = 1, NSCA, 1
  628. MPOVG5.VPOCHA(NLCF,I1)=IFLU1.FLUX(5+I1+NESP)
  629. & * SURF
  630. ENDDO
  631. C
  632. C******* Calcul du pas du temps (CFL)
  633. C
  634. C****** a) etat a l'interface
  635. C
  636. DIAMG = MPOVDI.VPOCHA(NLCEG,1)
  637. DIAMD = MPOVDI.VPOCHA(NLCED,1)
  638. DIAM = MIN(DIAMG,DIAMD)
  639. CELL = 1.0D0/DIAM/CELLT
  640. IF(CELL .GT. UNSDT)THEN
  641. UNSDT = CELL
  642. DIAMEL = DIAMG
  643. NLCEMI = NLCEG
  644. ENDIF
  645. C
  646. C****** b) etat gauche
  647. C
  648. ASON = SQRT(GAMG*PG/ROG)
  649. LAMBDA = ABS(UNG) + ASON
  650. CELL = LAMBDA / DIAM
  651. IF(CELL .GT. UNSDT)THEN
  652. UNSDT = CELL
  653. DIAMEL = DIAMG
  654. NLCEMI = NLCEG
  655. ENDIF
  656. C
  657. C****** C) etat droite
  658. C
  659. ASON = SQRT(GAMD*PD/ROD)
  660. LAMBDA = ABS(UND) + ASON
  661. CELL = LAMBDA / DIAM
  662. IF(CELL .GT. UNSDT)THEN
  663. UNSDT = CELL
  664. DIAMEL = DIAMD
  665. NLCEMI = NLCED
  666. ENDIF
  667. C
  668. C
  669. C**** Fin boucle sur FACEL
  670. C
  671. ENDDO
  672. C
  673. C**** Pas du temps (condition de non interaction en 1D)
  674. C
  675. DT = 0.5D0 / UNSDT
  676. C
  677. C**** Desactivation des segments et
  678. C on detruit les MCHAMLs
  679. C
  680. C
  681. C**** SEGSUP FRAMAG
  682. C SEGSUP FRAMAD
  683. C
  684. C meme si LOGME = .FALSE.
  685. C
  686. SEGSUP FRAMAG
  687. SEGSUP FRAMAD
  688. C
  689. SEGSUP MLENT1
  690. SEGDES MLENT2
  691. SEGDES IPT2
  692. C
  693. SEGSUP IFLU1
  694. SEGSUP IFLU2
  695. C
  696. SEGDES MPOVSU
  697. SEGDES MPOVDI
  698. C
  699. SEGDES MPOVG1
  700. SEGDES MPOVG2
  701. SEGDES MPOVG3
  702. C
  703. SEGDES MELRO
  704. SEGDES MELP
  705. SEGDES MELVUN
  706. SEGDES MELVUT
  707. SEGDES MELVNX
  708. SEGDES MELVNY
  709. SEGDES MELT1X
  710. SEGDES MELT1Y
  711. C
  712. IF(LOGME) THEN
  713. DO I1 = 1, NESP
  714. MELVA1 = MCHAMY.IELVAL(I1)
  715. SEGDES MELVA1
  716. ENDDO
  717. SEGDES MPOVG4
  718. C
  719. SEGDES MCHAMY
  720. ENDIF
  721. IF(ISCAF .GT. 0) THEN
  722. DO I1 = 1, NSCA, 1
  723. MELVA1 = MCHAMS.IELVAL(I1)
  724. SEGDES MELVA1
  725. ENDDO
  726. SEGDES MPOVG5
  727. C
  728. SEGDES MCHAMS
  729. ENDIF
  730. C
  731. 9999 CONTINUE
  732. C
  733. RETURN
  734. END
  735. C
  736.  
  737.  
  738.  
  739.  
  740.  
  741.  
  742.  
  743.  
  744.  
  745.  
  746.  
  747.  
  748.  
  749.  
  750.  

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