Télécharger la plaquette Cast3M


fr en

Télécharger kgfm12.eso

Retour à la liste

Numérotation des lignes : 

kgfm12

  1. C KGFM12    SOURCE    CB215821  20/11/25    13:31:25     10792          
  2.       SUBROUTINE KGFM12(
  3.      &        LOGCON, NESP,
  4.      &        INDMET,
  5.      &        MLRMGA, MLRMPI, MLRPGA, MLRPPI,
  6.      &        PMIN,
  7.      &        ICHPHI, ICHALP,
  8.      &        IPHI1,IROF1,IVITF1,IPF1,
  9.      &        IFRMAF, IFRALF,
  10.      &        ICHPSU,ICHPDI,ICHPVO,
  11.      &        MELEMC,MELEMF,MELEFE,MELLIM,
  12.      &        ICHRES,
  13.      &        DT,
  14.      &        LOGNC,LOGAN,MESERR)
  15. C************************************************************************
  16. C
  17. C PROJET            :  CASTEM 2000
  18. C
  19. C NOM               :  KGFM12
  20. C
  21. C DESCRIPTION       :  Voir KONV15
  22. C
  23. C                      Cas deux dimensions, GFMP
  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 1) Parameters concerning some gas properties
  34. C
  35. C    GAM1, PINF1, GAM2, PINF2 : gamma and pinf for the "stiffened gases".
  36. C
  37. C 2) Numerical scheme
  38. C
  39. C INDMET  :     1     GODUNOV (exact Riemann solver)
  40. C
  41. C 3) Pointers for the CHPOINTs/MCHAML des variables
  42. C
  43. C ICHPHI   : chpoint CENTRE contenant (phi)
  44. C
  45. C IPHI1   : MCHAML sur "FACEL"  contenant phi
  46. C           ("gauche" et "droite");
  47. C
  48. C IROF1   : MCHAML sur "FACEL"  contenant la masse volumique 1
  49. C           ("gauche" et "droite");
  50. C
  51. C IVITF1  : MCHAML sur "FACEL"  contenant la vitesse 1 dans le repaire
  52. C           local (n,t) et les cosinus directeurs des repaire local;
  53. C
  54. C IPF1    : MCHAML sur "FACEL"  contenant la pression 1;
  55. C
  56. C
  57. C 4) Others
  58. C
  59. C ICHPSU  : CHPOINT "FACE" contenant la surface des faces
  60. C
  61. C ICHPDI  : CHPOINT "CENTRE" contenant le diametre minimum
  62. C           de chaque element
  63. C
  64. C ICHPVO  : CHPOINT "CENTRE" contenant le volume
  65. C           de chaque element
  66. C
  67. C 5) Pointeurs de MELEME de la table DOMAINE
  68. C
  69. C MELEMC  : MELEME 'CENTRE' du SPG des CENTRES
  70. C
  71. C MELEMF  : MELEME 'FACE' du SPG des FACES
  72. C
  73. C MELEFE  : MELEME 'FACEL' du connectivité FACES -> ELEM
  74. C
  75. C MELLIM  : MAILLAGE where the flux (or residu) will not be found
  76. C
  77. C SORTIES (il faudrait dire E/S)
  78. C
  79. C ICHRES  : pointeurs de CHPOINTs "FACE" des flux aux interfaces:
  80. C
  81. C DT      : pas de temps pour le respect de la CFL-like condition
  82. C           DT < DIAMEL /2 /max(Lambda_i)
  83. C           En maillage regulier cette condition garantie la
  84. C           non-interaction des ondes
  85. C
  86. C LOGNC   : (LOGICAL): si .TRUE. la methode de Newton-Rapson, utilisée
  87. C           dans pour la solution du probleme Riemann n'a pas bien
  88. C           marchéee
  89. C
  90. C LOGAN   : (LOGICAL): si .TRUE. une anomalie à été detectée
  91. C
  92. C MESERR  : pour l'ecriture des messages d'erreurs
  93. C
  94. C************************************************************************
  95. C
  96. C HISTORIQUE (Anomalies et modifications éventuelles)
  97. C
  98. C  04/12/2010 - Created
  99. *  25/05/2011 - Evolution in CAST3M
  100. C
  101. C************************************************************************
  102. C
  103. C
  104. C N.B.: On suppose qu'on a déjà controllé RO, P > 0
  105. C                                         Y \in (0,1)
  106. C       Si non il faut le faire!!!
  107. C
  108. C************************************************************************
  109. C
  110. C
  111.       IMPLICIT INTEGER(I-N)
  112.       INTEGER NESP, NESP1, IESP, INDMET, NPAR
  113.      &     , ICHPHI, ICHALP
  114.      &     , IPHI1, IROF1, IVITF1, IPF1
  115.      &     , ICHPSU, ICHPDI, ICHPVO, MELEMC, MELEMF, MELEFE, MELLIM
  116.      &     , IGEOMC, IGEOMF
  117.      &     , ICHRES
  118.      &     , NFAC
  119.      &     , NLCF, NGCEG,  NGCED, NLCEG, NLCED
  120.      &     , NGCF, NLCF1, SPG1, SPG2, NLFL
  121.      &     , NFRA
  122.       REAL*8 GAMMA, PINF, ALP, ALPI, RNUM, DEN, PMIN
  123.       REAL*8 DT, UNSDT, CELLT, CELLTM
  124.      &     , PHIG1C, PHIG1F, ROG1, UNG1, UTG1, PG1
  125.      &     , PHID1C, PHID1F, ROD1, UND1, UTD1, PD1
  126.      &     , SURF,CNX, CNY, CTX , CTY
  127.      &     , CELL, DIAMG, DIAMD, DIAM, VOLG, VOLD
  128.      &     , FLUCEL, FLUCED
  129.      &     , PHIPRO, PINFG, GAMG, PINFD, GAMD
  130.       REAL*8 RHO_B, P_B, U_B, D_B
  131.      &     ,RHO_A, P_A, U_A, D_A
  132.      &     ,RHO_B1, P_B1, U_B1, D_B1
  133.      &     ,RHO_A1, P_A1, U_A1, D_A1
  134.      &     ,ROG1G, PG1G, UNG1G
  135.      &     ,ROD1G, PD1G, UND1G
  136.      &     ,UINT
  137. C
  138.       LOGICAL LOGNC, LOGAN, LOGVAC, LOGDEB, LOGCON
  139.       CHARACTER*(40) MESERR
  140.       CHARACTER*(8) TYPE
  141.       PARAMETER (LOGDEB = .FALSE.)
  142. C
  143. C**** LES INCLUDES
  144. C
  145.  
  146. -INC PPARAM
  147. -INC CCOPTIO
  148. -INC SMCOORD
  149. -INC SMELEME
  150. -INC SMCHAML
  151. C     Normal and tangent
  152.       POINTEUR MELVNX.MELVAL, MELVNY.MELVAL,
  153.      &         MELT1X.MELVAL, MELT1Y.MELVAL
  154.       POINTEUR MELUN1.MELVAL, MELUT1.MELVAL
  155.       POINTEUR MELPH1.MELVAL, MELRO1.MELVAL, MELP1.MELVAL
  156.       POINTEUR MCHAMY.MCHAML, MCHAMA.MCHAML
  157. -INC SMCHPOI
  158.       POINTEUR MPOVSU.MPOVAL, MPOVDI.MPOVAL, MPOVVO.MPOVAL
  159.      &        , MPORES.MPOVAL
  160.      &        , MPPHI1.MPOVAL, MPALPH.MPOVAL
  161. -INC SMLMOTS
  162. -INC SMLENTI
  163.       POINTEUR MLELIM.MLENTI
  164.       INTEGER JG
  165. -INC SMLREEL
  166.       POINTEUR MLRMGA.MLREEL, MLRPGA.MLREEL,
  167.      &     MLRMPI.MLREEL, MLRPPI.MLREEL
  168. C
  169. C**** Segment for the subroutine for the flux
  170. C
  171.       POINTEUR IFLU.MLREEL, WPARG.MLREEL, WPARD.MLREEL
  172. C
  173. C**** Pointer for Y
  174.       SEGMENT FRAMAS
  175.          REAL*8 YET(NFRA)
  176.       ENDSEGMENT
  177.       POINTEUR FRAMAG.FRAMAS, FRAMAD.FRAMAS
  178.      &     , ALPHAG.FRAMAS, ALPHAD.FRAMAS
  179.      &     , ALPHCG.FRAMAS, ALPHCD.FRAMAS
  180. C
  181.       LOGAN = .FALSE.
  182. C
  183. C**** Parameter (phi + fractions massiques)
  184. C
  185.       NPAR = 1 + NESP
  186.       JG = NPAR
  187.       SEGINI WPARG
  188.       SEGINI WPARD
  189. C
  190.       SEGACT MLRMGA
  191.       SEGACT MLRPGA
  192.       SEGACT MLRMPI
  193.       SEGACT MLRPPI
  194. C
  195. C**** KRIPAD pour la correspondance global/local des conditions limits
  196. C
  197.       CALL KRIPAD(MELLIM,MLELIM)
  198. C      SEGINI MLELIM
  199. C
  200. C**** Initialisation des MCHAMLs
  201. C
  202. C     PHI
  203. C
  204.       MCHEL1 = IPHI1
  205.       SEGACT MCHEL1
  206.       MCHAM1 = MCHEL1.ICHAML(1)
  207.       SEGACT MCHAM1
  208.       MELPH1 = MCHAM1.IELVAL(1)
  209.       SEGDES MCHEL1
  210.       SEGDES MCHAM1
  211. C
  212. C     Masse volumique
  213. C
  214.       MCHEL1 = IROF1
  215.       SEGACT MCHEL1
  216.       MCHAM1 = MCHEL1.ICHAML(1)
  217.       SEGACT MCHAM1
  218.       MELRO1 = MCHAM1.IELVAL(1)
  219.       SEGDES MCHEL1
  220.       SEGDES MCHAM1
  221. C
  222. C     Pression
  223. C
  224.       MCHEL1 = IPF1
  225.       SEGACT MCHEL1
  226.       MCHAM1 = MCHEL1.ICHAML(1)
  227.       SEGACT MCHAM1
  228.       MELP1 = MCHAM1.IELVAL(1)
  229.       SEGDES MCHEL1
  230.       SEGDES MCHAM1
  231. C
  232. C     Vitesse et cosinus directeurs du repere (n,t)
  233. C
  234.       MCHEL1 = IVITF1
  235.       SEGACT MCHEL1
  236. C
  237. C     La vitesse a comme SPG MELEFE
  238. C     Le cosinus directeurs ont comme SPG MELEMF
  239. C
  240. C     MCHAM1 -> Cosinus directeurs
  241. C     MCHAM2 -> Vitesse
  242. C
  243.       SPG1 = MCHEL1.IMACHE(1)
  244.       SPG2 = MCHEL1.IMACHE(2)
  245.       IF((SPG1 .EQ. MELEMF)  .AND. (SPG2 .EQ. MELEFE))THEN
  246.          MCHAM1 = MCHEL1.ICHAML(1)
  247.          MCHAM2 = MCHEL1.ICHAML(2)
  248.       ELSEIF((SPG1 .EQ. MELEFE)  .AND. (SPG2 .EQ. MELEMF))THEN
  249.          MCHAM1 = MCHEL1.ICHAML(2)
  250.          MCHAM2 = MCHEL1.ICHAML(1)
  251.       ELSE
  252.          LOGAN = .TRUE.
  253.          CALL ERREUR(251)
  254.          GOTO 9999
  255.       ENDIF
  256.       SEGACT MCHAM1
  257.       MELVNX = MCHAM1.IELVAL(1)
  258.       MELVNY = MCHAM1.IELVAL(2)
  259.       MELT1X = MCHAM1.IELVAL(3)
  260.       MELT1Y = MCHAM1.IELVAL(4)
  261.       SEGDES MCHAM1
  262.       SEGACT MCHAM2
  263.       MELUN1 = MCHAM2.IELVAL(1)
  264.       MELUT1 = MCHAM2.IELVAL(2)
  265.       SEGDES MCHAM2
  266.       SEGDES MCHEL1
  267. C
  268.       IF (NESP .GE. 1) THEN
  269. C
  270. C        Fractions massiques au centre (non-conservative approach)
  271. C
  272.          CALL LICHT(ICHALP, MPALPH, TYPE, IGEOMC)
  273. C        SEGACT MPALPH
  274.  
  275.          NFRA = NESP
  276.          SEGINI FRAMAG
  277.          SEGINI FRAMAD
  278.          SEGINI ALPHAG
  279.          SEGINI ALPHAD
  280.          SEGINI ALPHCG
  281.          SEGINI ALPHCD
  282. C
  283.          MCHEL1 = IFRMAF
  284.          SEGACT MCHEL1
  285.          MCHAMY = MCHEL1.ICHAML(1)
  286.          SEGACT MCHAMY
  287.          SEGDES MCHEL1
  288.          NESP1 =  MCHAMY.IELVAL(/1)
  289.          IF (NESP1 .NE. NESP) THEN
  290.             WRITE(IOIMP,*) 'Mass fractions, dimension = ???'
  291.             CALL ERREUR(21)
  292.             GOTO 9999
  293.          ENDIF
  294. C
  295.          DO IESP = 1, NESP, 1
  296.             MELVA1 = MCHAMY.IELVAL(IESP)
  297.             SEGACT MELVA1
  298.          ENDDO
  299. C
  300.          MCHEL1 = IFRALF
  301.          SEGACT MCHEL1
  302.          MCHAMA = MCHEL1.ICHAML(1)
  303.          SEGACT MCHAMA
  304.          SEGDES MCHEL1
  305.          NESP1 =  MCHAMA.IELVAL(/1)
  306.          IF (NESP1 .NE. NESP) THEN
  307.             WRITE(IOIMP,*) 'Volume fractions dimension = ???'
  308.             CALL ERREUR(21)
  309.             GOTO 9999
  310.          ENDIF
  311. C
  312.          DO IESP = 1, NESP, 1
  313.             MELVA1 = MCHAMA.IELVAL(IESP)
  314.             SEGACT MELVA1
  315. *            write(*,*) melva1
  316.          ENDDO
  317.       ENDIF
  318. C
  319. C      write(*,*) 'Son in kgdm12'
  320. C      write(*,*) 'Ho letto gli MCHAM'
  321. C     Fluxes
  322. C
  323.       JG = IDIM + 3 + NESP
  324.       SEGINI IFLU
  325. C
  326. C**** Initialisation des MELEMEs
  327. C
  328. C     'CENTRE', 'FACEL'
  329. C
  330.       IPT2 = MELEFE
  331.       SEGACT IPT2
  332.       NFAC = IPT2.NUM(/2)
  333. C
  334. C**** KRIPAD pour la correspondance global/local de centre
  335. C
  336.       CALL KRIPAD(MELEMC,MLENT1)
  337. C
  338. C**** MLENTI1 a MCORD.XCOORD(/1)/(IDIM+1) elements
  339. C
  340. C     Si i est le numero global d'un noeud de ICEN,
  341. C     MLENT1.LECT(i) contient sa position, i.e.
  342. C
  343. C     I              = numero global du noeud centre
  344. C     MLENT1.LECT(i) = numero local du noeud centre
  345. C
  346. C     MLENT1 déjà activé, i.e.
  347. C
  348. C     SEGINI MLENT1
  349. C
  350. C
  351. C**** KRIPAD pour la correspondance global/local de 'FACE'
  352. C
  353.       CALL KRIPAD(MELEMF,MLENT2)
  354. C      SEGINI MLENT2
  355. C
  356. C**** CHPOINTs de la table DOMAINE
  357. C
  358.       CALL LICHT(ICHPSU,MPOVSU,TYPE,IGEOMF)
  359.       CALL LICHT(ICHPDI,MPOVDI,TYPE,IGEOMC)
  360.       CALL LICHT(ICHPVO,MPOVVO,TYPE,IGEOMC)
  361. C
  362. C**** LICHT active les MPOVALs en *MOD
  363. C
  364. C     i.e.
  365. C
  366. C     SEGACT MPOVSU*MOD
  367. C     SEGACT MPOVDI*MOD
  368. C     SEGACT MPOVVO*MOD
  369. C
  370. C
  371. C**** Les FLUX aux face
  372. C
  373.       CALL LICHT(ICHRES,MPORES,TYPE,IGEOMC)
  374. C     SEGACT MPORES*MOD
  375. C
  376. C**** Activation des MCHAMLs
  377. C
  378.       SEGACT MELPH1
  379.       SEGACT MELRO1
  380.       SEGACT MELP1
  381.       SEGACT MELUN1
  382.       SEGACT MELUT1
  383.       SEGACT MELVNX
  384.       SEGACT MELVNY
  385.       SEGACT MELT1X
  386.       SEGACT MELT1Y
  387. C
  388. C**** Lecture des MPOVALs ICHPHI, ICHPA2
  389. C
  390.       CALL LICHT(ICHPHI, MPPHI1, TYPE, IGEOMC)
  391. C     SEGACT MPPHI1
  392. C
  393. C**** Initialisation de 1/DT
  394. C
  395.       UNSDT = 0.0D0
  396. C
  397. C**** BOUCLE SUR FACEL pour le calcul du FLUX/RESIDU
  398. C
  399.       DO NLCF = 1, NFAC
  400. C
  401. C******* NLCF  = numero local du centre de facel
  402. C        NGCF  = numero global du centre de facel
  403. C        NLCF1 = numero local du centre de face
  404. C        NGCEG = numero global du centre ELT "gauche"
  405. C        NLCEG = numero local du centre ELT "gauche"
  406. C        NGCED = numero global du centre ELT "droite"
  407. C        NLCED = numero local du centre ELT "droite"
  408. C
  409.          NGCEG = IPT2.NUM(1,NLCF)
  410.          NGCED = IPT2.NUM(3,NLCF)
  411.          NGCF  = IPT2.NUM(2,NLCF)
  412.          NLCF1 = MLENT2.LECT(NGCF)
  413.          NLCEG = MLENT1.LECT(NGCEG)
  414.          NLCED = MLENT1.LECT(NGCED)
  415.          NLFL  = MLELIM.LECT(NGCF)
  416. C
  417.          VOLG = MPOVVO.VPOCHA(NLCEG,1)
  418.          VOLD = MPOVVO.VPOCHA(NLCED,1)
  419. C
  420. C******* NLCF != NLCF1 -> l'auteur (MOI) n'a rien compris.
  421. C
  422.          IF(NLCF .NE. NLCF1)THEN
  423.             MESERR = 'Il ne faut pas jouer avec la console.   '
  424.             LOGAN = .TRUE.
  425.             CALL ERREUR(5)
  426.             GOTO 9999
  427.          ENDIF
  428.          IF(NLFL .EQ. 0)THEN
  429.             CELLTM = 0.0D0
  430. C
  431. C********** Recuperation des Etats "gauche" et "droite"
  432. C
  433.             PHIG1C = MPPHI1.VPOCHA(NLCEG,1)
  434.             PHIG1F = MELPH1.VELCHE(1,NLCF)
  435.             ROG1 = MELRO1.VELCHE(1,NLCF)
  436.             UNG1 = MELUN1.VELCHE(1,NLCF)
  437.             UTG1 = MELUT1.VELCHE(1,NLCF)
  438.             PG1 = MELP1.VELCHE(1,NLCF)
  439.             DO IESP = 1, NESP, 1
  440.                MELVA1 = MCHAMY.IELVAL(IESP)
  441.                FRAMAG.YET(IESP) = MELVA1.VELCHE(1,NLCF)
  442.             ENDDO
  443.             DO IESP = 1, NESP, 1
  444.                MELVA1 = MCHAMA.IELVAL(IESP)
  445.                ALPHAG.YET(IESP) = MELVA1.VELCHE(1,NLCF)
  446. C               IF (NGCF .EQ. 685) THEN
  447. C                  write(*,*) 'IFRALF =', IFRALF
  448. C                  write(*,*) 'MELVA1 =', MELVA1
  449. C                  write(*,*) 'NLCF =', NLCF
  450. C                  write(*,*) 'IESP =', IESP
  451. C                  write(*,*) 'ALPHAG.YET(IESP) =', MELVA1.VELCHE(1,NLCF)
  452. C               ENDIF
  453.             ENDDO
  454.             DO IESP = 1, NESP, 1
  455.                ALPHCG . YET(IESP) = MPALPH.VPOCHA(NLCEG,IESP)
  456.             ENDDO
  457. C
  458.             WPARG.PROG(1) = PHIG1F
  459. C           Fractions massiques traité en forme conservative...
  460.             DO IESP = 1, NESP, 1
  461.                WPARG.PROG (IESP + 1) = FRAMAG.YET(IESP)
  462.             ENDDO
  463. C
  464.             PHID1C = MPPHI1.VPOCHA(NLCED,1)
  465.             PHID1F = MELPH1.VELCHE(3,NLCF)
  466.             ROD1 = MELRO1.VELCHE(3,NLCF)
  467.             UND1 = MELUN1.VELCHE(3,NLCF)
  468.             UTD1 = MELUT1.VELCHE(3,NLCF)
  469.             PD1 = MELP1.VELCHE(3,NLCF)
  470.             DO IESP = 1, NESP, 1
  471.                MELVA1 = MCHAMY.IELVAL(IESP)
  472.                FRAMAD.YET(IESP) = MELVA1.VELCHE(3,NLCF)
  473.             ENDDO
  474.             DO IESP = 1, NESP, 1
  475.                MELVA1 = MCHAMA.IELVAL(IESP)
  476.                ALPHAD.YET(IESP) = MELVA1.VELCHE(3,NLCF)
  477.             ENDDO
  478.             DO IESP = 1, NESP, 1
  479.                ALPHCD . YET(IESP) = MPALPH.VPOCHA(NLCED,IESP)
  480.             ENDDO
  481.             WPARD.PROG(1) = PHID1F
  482.             DO IESP = 1, NESP, 1
  483.                WPARD.PROG (IESP + 1) = FRAMAD.YET(IESP)
  484.             ENDDO
  485. C
  486.             CNX = MELVNX.VELCHE(1,NLCF)
  487.             CNY = MELVNY.VELCHE(1,NLCF)
  488.             CTX = MELT1X.VELCHE(1,NLCF)
  489.             CTY = MELT1Y.VELCHE(1,NLCF)
  490. C
  491.             SURF   = MPOVSU.VPOCHA(NLCF,1)
  492. C
  493.             PHIPRO = PHIG1C * PHID1C
  494. C
  495. C            IF (NGCF .EQ. 685) THEN
  496.             IF (LOGDEB) THEN
  497.                WRITE(*,*)
  498.                WRITE(*,*) 'NGCF =', NGCF
  499.                WRITE(*,*) 'Left '
  500.                WRITE(*,*) 'Phic, phi'
  501.                WRITE(*,*) PHIG1C, PHIG1F
  502.                WRITE(*,*) 'rho, un, ut, pg'
  503.                WRITE(*,*) ROG1, UNG1, UTG1, PG1
  504.                WRITE(*,*) 'Right'
  505.                WRITE(*,*) 'Phic, phi'
  506.                WRITE(*,*) PHID1C, PHID1F
  507.                WRITE(*,*) 'rho, un, ut, pg'
  508.                WRITE(*,*) ROD1, UND1, UTD1, PD1
  509.             ENDIF
  510. C
  511. C********** According to PHIG1C we compute GAMG, PINFG
  512. C
  513.             IF (PHIG1C .LT. 0.0D0) THEN
  514.                MLREE1 = MLRMGA
  515.                MLREE2 = MLRMPI
  516.             ELSE
  517.                MLREE1 = MLRPGA
  518.                MLREE2 = MLRPPI
  519.             ENDIF
  520. C           Left
  521.             ALP = 1.0D0
  522.             DEN = 0.0D0
  523.             RNUM = 0.0D0
  524.             DO IESP = 1, NESP, 1
  525.                GAMMA = MLREE1.PROG(IESP)
  526.                PINF =  MLREE2.PROG(IESP)
  527.                ALPI = ALPHAG.YET(IESP)
  528. C               IF (NGCF .EQ. 685) THEN
  529. C                  write(*,*) 'IESP =', IESP
  530. C                  write(*,*) 'ALPI =', ALPI
  531. C               ENDIF
  532. C
  533.                ALP = ALP - ALPI
  534.                DEN = DEN + (ALPI / (GAMMA - 1.0D0))
  535.                RNUM = RNUM  +
  536.      &              ((ALPI * GAMMA * PINF) / (GAMMA - 1.0D0))
  537. C               IF (NGCF .EQ. 685) THEN
  538. C                  write(*,*) 'IESP =', IESP
  539. C                  write(*,*) 'GAMMA, PINF, ALPI, ALP, DEN, RNUM'
  540. C                  write(*,*) GAMMA, PINF, ALPI, ALP, DEN, RNUM
  541. C               ENDIF
  542.             ENDDO
  543.             GAMMA = MLREE1.PROG(NESP + 1)
  544.             PINF =  MLREE2.PROG(NESP + 1)
  545.             DEN = DEN + (ALP / (GAMMA - 1.0D0))
  546.             RNUM = RNUM  + ((ALP * GAMMA * PINF) / (GAMMA - 1.0D0))
  547. C            IF (NGCF .EQ. 685) THEN
  548. C               write(*,*) 'IESP =', (NESP + 1)
  549. C               write(*,*) 'GAMMA, PINF, ALPI, ALP, DEN, RNUM'
  550. C               write(*,*) GAMMA, PINF, ALPI, ALP, DEN, RNUM
  551. C            ENDIF
  552. C
  553.             PINF = RNUM / DEN
  554.             GAMMA = 1.0D0 / DEN
  555.             GAMMA = GAMMA + 1.0D0
  556.             PINF = PINF / GAMMA
  557. C
  558.             PINFG = PINF
  559.             GAMG = GAMMA
  560. C
  561. C********** According to PHID1C we compute GAMD, PINFD
  562. C
  563.             IF (PHID1C .LT. 0.0D0) THEN
  564.                MLREE1 = MLRMGA
  565.                MLREE2 = MLRMPI
  566.             ELSE
  567.                MLREE1 = MLRPGA
  568.                MLREE2 = MLRPPI
  569.             ENDIF
  570.             ALP = 1.0D0
  571.             DEN = 0.0D0
  572.             RNUM = 0.0D0
  573.             DO IESP = 1, NESP, 1
  574.                GAMMA = MLREE1.PROG(IESP)
  575.                PINF =  MLREE2.PROG(IESP)
  576.                ALPI = ALPHAD.YET(IESP)
  577.                ALP = ALP - ALPI
  578.                DEN = DEN + (ALPI / (GAMMA - 1.0D0))
  579.                RNUM = RNUM  +
  580.      &              ((ALPI * GAMMA * PINF) / (GAMMA - 1.0D0))
  581.             ENDDO
  582.             GAMMA = MLREE1.PROG(NESP + 1)
  583.             PINF =  MLREE2.PROG(NESP + 1)
  584.             DEN = DEN + (ALP / (GAMMA - 1.0D0))
  585.             RNUM = RNUM  + ((ALP * GAMMA * PINF) / (GAMMA - 1.0D0))
  586. C
  587.             PINF = RNUM / DEN
  588.             GAMMA = 1.0D0 / DEN
  589.             GAMMA = GAMMA + 1.0D0
  590.             PINF = PINF / GAMMA
  591. C
  592.             PINFD = PINF
  593.             GAMD = GAMMA
  594. C
  595. C           A priori, there is no vacuum...
  596. C
  597.             LOGVAC = .FALSE.
  598. C
  599. C            IF (NGCF .EQ. 685) THEN
  600. C               WRITE(*,*)
  601. C     &              PINFG, GAMG,
  602. C     &              PINFD, GAMD
  603. C            ENDIF
  604. C
  605.             IF (INDMET .EQ. 1) THEN
  606. C
  607. C                 Exact Riemann solver
  608. C                 Computation of the states on the left and on the
  609. C                 right of the contact discontinuity
  610. C
  611.                CALL FSTIFF(IDIM, NPAR,
  612.      &              PINFG, GAMG,
  613.      &              PINFD, GAMD,
  614.      &              PMIN,
  615.      &              ROG1, PG1, UNG1, UTG1, 0.0D0, WPARG.PROG,
  616.      &              ROD1, PD1, UND1, UTD1, 0.0D0, WPARD.PROG,
  617.      &              RHO_B, P_B, U_B, D_B,
  618.      &              RHO_A, P_A, U_A, D_A,
  619.      &              0.0D0,
  620.      &              UINT,
  621.      &              IFLU.PROG, CELLT,
  622.      &              LOGAN, LOGNC, LOGVAC)
  623.                IF (LOGAN) THEN
  624.                   WRITE(IOIMP,*) 'SUBROUTINE KGFM12'
  625.                   WRITE(IOIMP,*) 'ANOMALY DETECTED'
  626. C
  627.                   WRITE(*,*)
  628.                   WRITE(*,*) 'NGCF =', NGCF
  629.                   WRITE(*,*) 'NLCF =', NLCF
  630.                   WRITE(*,*) 'Left '
  631.                   WRITE(*,*) 'Phic, phi'
  632.                   WRITE(*,*) PHIG1C, PHIG1F
  633.                   WRITE(*,*) 'Pinf, gamma'
  634.                   WRITE(*,*) PINFG, GAMG
  635.                   WRITE(*,*) 'rho, un, ut, pg'
  636.                   WRITE(*,*) ROG1, UNG1, UTG1, PG1
  637.                   WRITE(*,*) 'Right'
  638.                   WRITE(*,*) 'Phic, phi'
  639.                   WRITE(*,*) PHID1C, PHID1F
  640.                   WRITE(*,*) 'Pinf, gamma'
  641.                   WRITE(*,*) PINFD, GAMD
  642.                   WRITE(*,*) 'rho, un, ut, pg'
  643.                   WRITE(*,*) ROD1, UND1, UTD1, PD1
  644. C
  645.                   CALL ERREUR(5)
  646.                   GOTO 9999
  647.                ENDIF
  648.             ENDIF
  649.             CELLTM = MAX(CELLT, CELLTM)
  650. C
  651.             IF ((.NOT. LOGVAC) .AND. (PHIPRO .GT. 0.0D0)) THEN
  652. C
  653. C************* Ecriture du residu
  654. C
  655. C              FLUX(2) = RO Un              RO Un
  656. C              FLUX(3) = RO Un Un + P   ->  RO Un Ux + P CNX
  657. C              FLUX(4) = RO Un Ut       ->  RO Un Uy + P CNY
  658. C              FLUX(5) = RO Un Et           RO Un Et
  659. C              Mass fractions... (5 + IESP)...
  660. C              Volume fractions in non conservative form...
  661. C
  662.                FLUCEL = IFLU.PROG(1) * SURF
  663. C               write(*,*) 'flux', iflu.prog(1)
  664. C               write(*,*) 'nlceg, nlced', nlceg, nlced
  665. C               write(*,*) 'surf', surf
  666.                MPORES.VPOCHA(NLCEG,2) = MPORES.VPOCHA(NLCEG,2) -
  667.      &              (FLUCEL / VOLG)
  668.                IF (NLCED .NE. NLCEG)
  669.      &              MPORES.VPOCHA(NLCED,2) = MPORES.VPOCHA(NLCED,2) +
  670.      &              (FLUCEL / VOLD)
  671.                FLUCEL = (IFLU.PROG(2)*CNX+IFLU.PROG(3)*CTX) * SURF
  672.                MPORES.VPOCHA(NLCEG,3) = MPORES.VPOCHA(NLCEG,3) -
  673.      &              (FLUCEL / VOLG)
  674.                IF (NLCED .NE. NLCEG)
  675.      &              MPORES.VPOCHA(NLCED,3) = MPORES.VPOCHA(NLCED,3) +
  676.      &              (FLUCEL / VOLD)
  677.                FLUCEL = (IFLU.PROG(2)*CNY+IFLU.PROG(3)*CTY) * SURF
  678.                MPORES.VPOCHA(NLCEG,4) = MPORES.VPOCHA(NLCEG,4) -
  679.      &              (FLUCEL / VOLG)
  680.                IF (NLCED .NE. NLCEG)
  681.      &              MPORES.VPOCHA(NLCED,4) = MPORES.VPOCHA(NLCED,4) +
  682.      &              (FLUCEL / VOLD)
  683.                FLUCEL = IFLU.PROG(4) * SURF
  684.                MPORES.VPOCHA(NLCEG,5) = MPORES.VPOCHA(NLCEG,5) -
  685.      &              (FLUCEL / VOLG)
  686.                IF (NLCED .NE. NLCEG)
  687.      &              MPORES.VPOCHA(NLCED,5) = MPORES.VPOCHA(NLCED,5) +
  688.      &              (FLUCEL / VOLD)
  689. C
  690. C************* rho Y in conservative form
  691. C              Components 5 + 1 -> 5 + NESP
  692. C
  693.                DO IESP = 1, NESP, 1
  694.  
  695.                   FLUCEL = IFLU.PROG(5 + IESP) * SURF
  696.                   MPORES.VPOCHA(NLCEG,5 + IESP) =
  697.      &                 MPORES.VPOCHA(NLCEG,5 + IESP) -
  698.      &                 (FLUCEL / VOLG)
  699.                   IF (NLCED .NE. NLCEG)
  700.      &                 MPORES.VPOCHA(NLCED,5 + IESP) =
  701.      &                 MPORES.VPOCHA(NLCED,5 + IESP) +
  702.      &                 (FLUCEL / VOLD)
  703.                ENDDO
  704. C************* Alpha in non conservative form
  705. C              Components 5 + NESP -> 5 + (2 * NESP)
  706.                IF (NESP .GE. 1) THEN
  707.                   IF (UINT .GT. 0.0D0) THEN
  708.                      FRAMAS = ALPHAG
  709.                   ELSE
  710.                      FRAMAS = ALPHAD
  711.                   ENDIF
  712.                   DO IESP = 1, NESP, 1
  713. C                     UINT > 0
  714. C                     FLUCEL = UINT * (PHIG1F - PHIG1C) * SURF
  715. C                     UINT < 0
  716. C                     FLUCEL = UINT * (PHID1F - PHIG1C) * SURF
  717.                      FLUCEL = UINT * (FRAMAS.YET(IESP) -
  718.      &                    ALPHCG.YET(IESP)) * SURF
  719. C                     UINT > 0
  720. C                     FLUCED = UINT * (PHIG1F - PHID1C) * SURF
  721. C                     UINT < 0
  722. C                     FLUCED = UINT * (PHID1F - PHID1C) * SURF
  723.                      FLUCED = UINT * (FRAMAS.YET(IESP) -
  724.      &                    ALPHCD.YET(IESP)) * SURF
  725. C
  726.                      MPORES.VPOCHA(NLCEG,5 + NESP + IESP) =
  727.      &                    MPORES.VPOCHA(NLCEG,5 + NESP + IESP) -
  728.      &                    (FLUCEL / VOLG)
  729.                      IF (NLCED .NE. NLCEG)
  730.      &                    MPORES.VPOCHA(NLCED,5 + NESP + IESP) =
  731.      &                    MPORES.VPOCHA(NLCED,5 + NESP + IESP) +
  732.      &                    (FLUCED / VOLD)
  733.                   ENDDO
  734.                ENDIF
  735.             ENDIF
  736. C
  737. C           Level set
  738. C
  739.             IF (LOGCON) THEN
  740. C
  741. C************* Conservative, variable rho phi...
  742. C
  743.                FLUCEL = IFLU.PROG(5) * SURF
  744. C            write(*,*) 'flux', iflu.prog(1)
  745. C            write(*,*) 'nlceg, nlced', nlceg, nlced
  746. C            write(*,*) 'surf', surf
  747.                MPORES.VPOCHA(NLCEG,1) = MPORES.VPOCHA(NLCEG,1) -
  748.      &              (FLUCEL / VOLG)
  749.                IF (NLCED .NE. NLCEG)
  750.      &              MPORES.VPOCHA(NLCED,1) = MPORES.VPOCHA(NLCED,1) +
  751.      &              (FLUCEL / VOLD)
  752.             ELSE
  753. C
  754. C************ Non conservative, variable phi
  755. C
  756.                IF (UINT .GT. 0.0D0) THEN
  757.                   FLUCEL = UINT * (PHIG1F - PHIG1C) * SURF
  758.                   FLUCED = UINT * (PHIG1F - PHID1C) * SURF
  759.                ELSE
  760.                   FLUCEL = UINT * (PHID1F - PHIG1C) * SURF
  761.                   FLUCED = UINT * (PHID1F - PHID1C) * SURF
  762.                ENDIF
  763. C
  764.                MPORES.VPOCHA(NLCEG,1) = MPORES.VPOCHA(NLCEG,1) -
  765.      &              (FLUCEL / VOLG)
  766.                IF (NLCED .NE. NLCEG)
  767.      &              MPORES.VPOCHA(NLCED,1) = MPORES.VPOCHA(NLCED,1) +
  768.      &              (FLUCED / VOLD)
  769.             ENDIF
  770. C
  771.             IF ((PHIPRO .LE. 0.0D0) .OR. LOGVAC) THEN
  772. *               write(*,*) 'Son qui !!!'
  773. *               stop
  774. C
  775. C           Left flux for everything but not for the level-set
  776. C
  777.                ROD1G = RHO_B
  778.                PD1G  = P_B
  779.                UND1G = U_B
  780.                IF (INDMET .EQ. 1) THEN
  781. C
  782. C              Exact Riemann solver
  783. C
  784.                   CALL FSTIFF(IDIM, NPAR,
  785.      &                 PINFG, GAMG,
  786.      &                 PINFG, GAMG,
  787.      &                 PMIN,
  788.      &                 ROG1, PG1, UNG1, UTG1, 0.0D0, WPARG.PROG,
  789.      &                 ROD1G, PD1G, UND1G, UTG1, 0.0D0, WPARG.PROG,
  790.      &                 RHO_B1, P_B1, U_B1, D_B1,
  791.      &                 RHO_A1, P_A1, U_A1, D_A1,
  792.      &                 0.0D0,
  793.      &                 UINT,
  794.      &                 IFLU.PROG, CELLT,
  795.      &                 LOGAN, LOGNC, LOGVAC)
  796. C                 IF (LOGVAC) THEN
  797. C                    WRITE(IOIMP,*) 'Vacuum of vacuum ???'
  798. C                    WRITE(IOIMP,*) 'SUBROUTINE KGFM12'
  799. C                    WRITE(IOIMP,*) 'ANOMALY DETECTED'
  800. C                    CALL ERREUR(251)
  801. C                    GOTO 9999
  802. C                 ENDIF
  803.                   IF (LOGAN) THEN
  804.                      WRITE(IOIMP,*) 'SUBROUTINE KGFM12'
  805.                      WRITE(IOIMP,*) 'ANOMALY DETECTED'
  806. C
  807.                      WRITE(*,*)
  808.                      WRITE(*,*) 'NGCF =', NGCF
  809.                      WRITE(*,*) 'Left '
  810.                      WRITE(*,*) 'Phic, phi'
  811.                      WRITE(*,*) PHIG1C, PHIG1F
  812.                      WRITE(*,*) 'rho, un, ut, pg'
  813.                      WRITE(*,*) ROG1, UNG1, UTG1, PG1
  814.                      WRITE(*,*) 'Right'
  815.                      WRITE(*,*) 'Phic, phi'
  816.                      WRITE(*,*) PHID1C, PHID1F
  817.                      WRITE(*,*) 'rho, un, ut, pg'
  818.                      WRITE(*,*) ROD1, UND1, UTD1, PD1
  819. C
  820.                      CALL ERREUR(5)
  821.                      GOTO 9999
  822.                   ENDIF
  823.                ENDIF
  824.                CELLTM = MAX(CELLT, CELLTM)
  825. C
  826. C********** Ecriture du residu... gauche only...
  827. C
  828. C           FLUX(2) = RO Un              RO Un
  829. C           FLUX(3) = RO Un Un + P   ->  RO Un Ux + P CNX
  830. C           FLUX(4) = RO Un Ut       ->  RO Un Uy + P CNY
  831. C           FLUX(5) = RO Un Et           RO Un Et
  832. C
  833.                FLUCEL = IFLU.PROG(1) * SURF
  834. C               write(*,*) 'flux', iflu.prog(1)
  835. C               write(*,*) 'nlceg, nlced', nlceg, nlced
  836. C               write(*,*) 'surf', surf
  837.                MPORES.VPOCHA(NLCEG,2) = MPORES.VPOCHA(NLCEG,2) -
  838.      &              (FLUCEL / VOLG)
  839.                FLUCEL = (IFLU.PROG(2)*CNX+IFLU.PROG(3)*CTX) * SURF
  840.                MPORES.VPOCHA(NLCEG,3) = MPORES.VPOCHA(NLCEG,3) -
  841.      &              (FLUCEL / VOLG)
  842.                FLUCEL = (IFLU.PROG(2)*CNY+IFLU.PROG(3)*CTY) * SURF
  843.                MPORES.VPOCHA(NLCEG,4) = MPORES.VPOCHA(NLCEG,4) -
  844.      &              (FLUCEL / VOLG)
  845.                FLUCEL = IFLU.PROG(4) * SURF
  846.                MPORES.VPOCHA(NLCEG,5) = MPORES.VPOCHA(NLCEG,5) -
  847.      &              (FLUCEL / VOLG)
  848. C
  849. C              Level set already done...
  850. C              Level set treated as the other variables (like in the
  851. C              following comment... does not work)
  852. C
  853. C               FLUCEL = IFLU.PROG(5) * SURF
  854. C               MPORES.VPOCHA(NLCEG,1) = MPORES.VPOCHA(NLCEG,1) -
  855. C     &              (FLUCEL / VOLG)
  856. C
  857. C
  858. C************* rho Y in conservative form
  859. C              Components 5 + 1 -> 5 + NESP
  860. C
  861.                DO IESP = 1, NESP, 1
  862.                   FLUCEL = IFLU.PROG(5 + IESP) * SURF
  863.                   MPORES.VPOCHA(NLCEG,5 + IESP) =
  864.      &                 MPORES.VPOCHA(NLCEG,5 + IESP) -
  865.      &                 (FLUCEL / VOLG)
  866.                ENDDO
  867. C************* Alpha in non conservative form
  868. C              Components 5 + NESP -> 5 + (2 * NESP)
  869.                IF (NESP .GE. 1) THEN
  870.                   IF (UINT .GT. 0.0D0) THEN
  871.                      FRAMAS = ALPHAG
  872.                   ELSE
  873.                      FRAMAS = ALPHAD
  874.                   ENDIF
  875.                   DO IESP = 1, NESP, 1
  876. C                     UINT > 0
  877. C                     FLUCEL = UINT * (PHIG1F - PHIG1C) * SURF
  878. C                     UINT < 0
  879. C                     FLUCEL = UINT * (PHID1F - PHIG1C) * SURF
  880.                      FLUCEL = UINT * (FRAMAS.YET(IESP) -
  881.      &                    ALPHCG.YET(IESP)) * SURF
  882. C
  883.                      MPORES.VPOCHA(NLCEG,5 + NESP + IESP) =
  884.      &                    MPORES.VPOCHA(NLCEG,5 + NESP + IESP) -
  885.      &                    (FLUCEL / VOLG)
  886.                   ENDDO
  887.                ENDIF
  888. C
  889. C********** Calcul du pas du temps (CFL)
  890. C
  891. C*********  a) etat a l'interface
  892. C
  893.                DIAMG = MPOVDI.VPOCHA(NLCEG,1)
  894.                DIAMD = MPOVDI.VPOCHA(NLCED,1)
  895.                DIAM = MIN(DIAMG,DIAMD)
  896.                CELL  = CELLTM/DIAM
  897.                IF(CELL .GT. UNSDT)THEN
  898.                   UNSDT  =  CELL
  899.                ENDIF
  900. C
  901. C           Right flux for everything but not for the level-set
  902. C
  903.                ROG1G = RHO_A
  904.                PG1G  = P_A
  905.                UNG1G = U_A
  906. C
  907.                IF (NLCED .NE. NLCEG) THEN
  908.                   IF (INDMET .EQ. 1) THEN
  909. C
  910. C                 Exact Riemann solver
  911. C
  912.                      CALL FSTIFF(IDIM, NPAR,
  913.      &                    PINFD, GAMD,
  914.      &                    PINFD, GAMD,
  915.      &                    PMIN,
  916.      &                    ROG1G, PG1G, UNG1G, UTD1, 0.0D0, WPARD.PROG,
  917.      &                    ROD1, PD1, UND1, UTD1, 0.0D0, WPARD.PROG,
  918.      &                    RHO_B1, P_B1, U_B1, D_B1,
  919.      &                    RHO_A1, P_A1, U_A1, D_A1,
  920.      &                    0.0D0,
  921.      &                    UINT,
  922.      &                    IFLU.PROG, CELLT,
  923.      &                    LOGAN, LOGNC, LOGVAC)
  924. C
  925. C                  IF (LOGVAC) THEN
  926. C                     WRITE(IOIMP,*) 'Vacuum of vacuum ???'
  927. C                     WRITE(IOIMP,*) 'SUBROUTINE KGFM12'
  928. C                     WRITE(IOIMP,*) 'ANOMALY DETECTED'
  929. C                     CALL ERREUR(251)
  930. C                     GOTO 9999
  931. C                  ENDIF
  932. C
  933.                      IF (LOGAN) THEN
  934.                         WRITE(IOIMP,*) 'SUBROUTINE KGFM12'
  935.                         WRITE(IOIMP,*) 'ANOMALY DETECTED'
  936.                         CALL ERREUR(5)
  937.                         GOTO 9999
  938.                      ENDIF
  939. C
  940.                   ENDIF
  941.                   CELLTM = MAX(CELLT, CELLTM)
  942.  
  943. C
  944. C************* Ecriture du flux a droite...
  945. C
  946. C              FLUX(2) = RO Un              RO Un
  947. C              FLUX(3) = RO Un Un + P   ->  RO Un Ux + P CNX
  948. C              FLUX(4) = RO Un Ut       ->  RO Un Uy + P CNY
  949. C              FLUX(5) = RO Un Et           RO Un Et
  950. C
  951.                   FLUCEL = IFLU.PROG(1) * SURF
  952. C               write(*,*) 'flux', iflu.prog(1)
  953. C               write(*,*) 'nlceg, nlced', nlceg, nlced
  954. C               write(*,*) 'surf', surf
  955.                   MPORES.VPOCHA(NLCED,2) = MPORES.VPOCHA(NLCED,2) +
  956.      &                 (FLUCEL / VOLD)
  957.                   FLUCEL = (IFLU.PROG(2)*CNX+IFLU.PROG(3)*CTX) * SURF
  958.  
  959.                   MPORES.VPOCHA(NLCED,3) = MPORES.VPOCHA(NLCED,3) +
  960.      &                 (FLUCEL / VOLD)
  961.                   FLUCEL = (IFLU.PROG(2)*CNY+IFLU.PROG(3)*CTY) * SURF
  962.                   MPORES.VPOCHA(NLCED,4) = MPORES.VPOCHA(NLCED,4) +
  963.      &              (FLUCEL / VOLD)
  964.                   FLUCEL = IFLU.PROG(4) * SURF
  965.                   MPORES.VPOCHA(NLCED,5) = MPORES.VPOCHA(NLCED,5) +
  966.      &                 (FLUCEL / VOLD)
  967. C
  968. C              Level set already done...
  969. C              Level set treated as the other variables (like in the
  970. C              following comment... does not work)
  971. C
  972. C               FLUCEL = IFLU.PROG(5) * SURF
  973. C               MPORES.VPOCHA(NLCED,1) = MPORES.VPOCHA(NLCED,1) +
  974. C     &              (FLUCEL / VOLD)
  975. C
  976. C
  977. C**************** rho Y in conservative form
  978. C                 Components 5 + 1 -> 5 + NESP
  979. C
  980.                   DO IESP = 1, NESP, 1
  981.  
  982.                      FLUCEL = IFLU.PROG(5 + IESP) * SURF
  983.                      MPORES.VPOCHA(NLCED,5 + IESP) =
  984.      &                    MPORES.VPOCHA(NLCED,5 + IESP) +
  985.      &                    (FLUCEL / VOLD)
  986.                   ENDDO
  987. C************* Alpha in non conservative form
  988. C              Components 5 + NESP -> 5 + (2 * NESP)
  989.                   IF (NESP .GE. 1) THEN
  990.                      IF (UINT .GT. 0.0D0) THEN
  991.                         FRAMAS = ALPHAG
  992.                      ELSE
  993.                         FRAMAS = ALPHAD
  994.                      ENDIF
  995.                      DO IESP = 1, NESP, 1
  996.                         FLUCED = UINT * (FRAMAS.YET(IESP) -
  997.      &                       ALPHCD.YET(IESP)) * SURF
  998. C
  999.                         MPORES.VPOCHA(NLCED,5 + NESP + IESP) =
  1000.      &                       MPORES.VPOCHA(NLCED,5 + NESP + IESP) +
  1001.      &                       (FLUCED / VOLD)
  1002.                      ENDDO
  1003.                   ENDIF
  1004.                ENDIF
  1005.             ENDIF
  1006.             DIAMG = MPOVDI.VPOCHA(NLCEG,1)
  1007.             DIAMD = MPOVDI.VPOCHA(NLCED,1)
  1008.             DIAM = MIN(DIAMG,DIAMD)
  1009.             CELL  = CELLTM/DIAM
  1010.             IF(CELL .GT. UNSDT)THEN
  1011.                UNSDT  =  CELL
  1012.             ENDIF
  1013.          ENDIF
  1014. C
  1015. C**** Fin boucle sur FACEL
  1016. C
  1017.       ENDDO
  1018. C
  1019. C**** Pas du temps (condition de non interaction en 1D)
  1020. C
  1021.       DT = 0.5D0 / UNSDT
  1022. C
  1023. C**** Desactivation des segments
  1024. C
  1025.       SEGSUP MLELIM
  1026.       SEGSUP MLENT1
  1027.       SEGSUP MLENT2
  1028. C
  1029.       SEGDES MPOVSU
  1030.       SEGDES MPOVDI
  1031.       SEGDES MPOVVO
  1032.       SEGDES MPORES
  1033.       SEGDES IPT2
  1034. C
  1035. C     Warning. MLRMFR = PGAS . 'MASSFRA'
  1036. C              MLRCHE = PGAS . 'CHEMCOEF'
  1037. C
  1038.       SEGSUP WPARG
  1039.       SEGSUP WPARD
  1040. C
  1041.       SEGSUP IFLU
  1042. C
  1043.       SEGDES MPPHI1
  1044.       SEGDES MELPH1
  1045.       SEGDES MELRO1
  1046.       SEGDES MELP1
  1047.       SEGDES MELUN1
  1048.       SEGDES MELUT1
  1049.       SEGDES MELVNX
  1050.       SEGDES MELVNY
  1051.       SEGDES MELT1X
  1052.       SEGDES MELT1Y
  1053. C
  1054.       SEGDES MLRMGA
  1055.       SEGDES MLRPGA
  1056.       SEGDES MLRMPI
  1057.       SEGDES MLRPPI
  1058. C
  1059.       IF (NESP .GE. 1)THEN
  1060. C
  1061.          SEGDES MPALPH
  1062. C
  1063.          SEGSUP FRAMAG
  1064.          SEGSUP FRAMAD
  1065.          SEGSUP ALPHAG
  1066.          SEGSUP ALPHAD
  1067.          SEGSUP ALPHCG
  1068.          SEGSUP ALPHCD
  1069. C
  1070.          DO IESP = 1, NESP
  1071.             MELVA1 = MCHAMY.IELVAL(IESP)
  1072.             SEGDES MELVA1
  1073.          ENDDO
  1074.          SEGDES MCHAMY
  1075. C
  1076.          DO IESP = 1, NESP
  1077.             MELVA1 = MCHAMA.IELVAL(IESP)
  1078.             SEGDES MELVA1
  1079.          ENDDO
  1080.          SEGDES MCHAMA
  1081.       ENDIF
  1082. C
  1083.  9999 CONTINUE
  1084.       RETURN
  1085.       END
  1086. CC
  1087.  
  1088.  
  1089.  
  1090.  
  1091.  

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