Sources Esope | Cast3M

kodfl1
C KODFL1    SOURCE    CB215821  20/11/25    13:31:45     10792          
      SUBROUTINE KODFL1(INDK0,NESP,NORD,TMAX,RUNIV,PROPHY,
     &        MLRCHE,MLRMFR,
     &        INDMET,
     &        ICHPA1, ICHPA2,
     &        IALF1,IROF1,IVITF1,IPF1,
     &        IALF2,IROF2,IVITF2,IPF2,
     &        IUINF1, IUINF2,
     &        K0, ICHPK0,IGRALP,EPS,
     &        ICHPSU,ICHPDI,ICHPVO,
     &        MELEMC,MELEMF,MELEFE,MELLIM,
     &        ICHRES,
     &        DT,SURFL,
     &        LOGNC,LOGAN,MESERR)
C************************************************************************
C
C PROJET            :  CASTEM 2000
C
C NOM               :  KODFL1
C
C DESCRIPTION       :  Voir KONV13
C
C                      Cas deux dimensions, DEM
C
C LANGAGE           :  FORTRAN 77 + ESOPE 2000 (avec estensions CISI)
C
C AUTEUR            :  A. BECCANTINI, DEN/DM2S/SFME/LTMF
C
C************************************************************************
C
C ENTREES
C
C 0) INDK0 = 1 => K0 = constant
C    INDK0 = 2 => K0 = variable
C
C 1) Parameters concerning some gas properties
C
C NESP    : nombre d'especes dans le melange.
C
C NORD    : ordre des polynoms du cv_i
C
C TMAX    : maximum temperature for cv expansion
C
C RUNIV   : universal constant for gases
C
C PROPHY  : thermodynamic properties for the gases
C
C MLRCHE  : LISTREEL with the coefS involved in the chemical
C           reaction
C
C MLRMFR  : LISTREEL with the mass fractions before or after
C           the chemical reaction
C
C 2) Numerical scheme
C
C INDMET  :     1        VLH (non-reactive) + SS (reactive)
C
C INDMET  :     2        SS  (non-reactive) + SS (reactive)
C
C INDMET  :     3        AUSM+up (non-reactive) + SS (reactive)
C
C 3) Pointers for the CHPOINTs/MCHAML des variables
C
C ICHPA1   : chpoint CENTRE contenant alpha1
C ICHPA2   : chpoint CENTRE contenant alpha2
C
C IALF1   : MCHAML sur "FACEL"  contenant alpha1
C           ("gauche" et "droite");
C
C IROF1   : MCHAML sur "FACEL"  contenant la masse volumique 1
C           ("gauche" et "droite");
C
C IVITF1  : MCHAML sur "FACEL"  contenant la vitesse 1 dans le repaire
C           local (n,t) et les cosinus directeurs des repaire local;
C
C IPF1    : MCHAML sur "FACEL"  contenant la pression 1;
C
C IALF2   : MCHAML sur "FACEL"  contenant alpha1
C           ("gauche" et "droite");
C
C IROF2   : MCHAML sur "FACEL"  contenant la masse volumique 1
C           ("gauche" et "droite");
C
C IVITF2  : MCHAML sur "FACEL"  contenant la vitesse 1 dans le repaire
C           local (n,t) et les cosinus directeurs des repaire local;
C
C IPF2    : MCHAML sur "FACEL"  contenant la pression 1;
C
C IUINF1   : CHPOINT, one component, "cut-off velocity" 1
C               0 if no Bas MACH
C
C IUINF2   : CHPOINT, one component, "cut-off velocity" 2
C
C
C 4) Others
C
C K0 or
C ICHPK0  : fundamental flame speed
C
C IGRALP  : CHPOINT "FACE" contenant grad(alpha)/alpha
C
C EPS     : "zero"
C
C ICHPSU  : CHPOINT "FACE" contenant la surface des faces
C
C ICHPDI  : CHPOINT "CENTRE" contenant le diametre minimum
C           de chaque element
C
C ICHPVO  : CHPOINT "CENTRE" contenant le volume
C           de chaque element
C
C 5) Pointeurs de MELEME de la table DOMAINE
C
C MELEMC  : MELEME 'CENTRE' du SPG des CENTRES
C
C MELEMF  : MELEME 'FACE' du SPG des FACES
C
C MELEFE  : MELEME 'FACEL' du connectivité FACES -> ELEM
C
C MELLIM  : MAILLAGE where the flux (or residu) will not be found
C
C SORTIES (il faudrait dire E/S)
C
C ICHRES  : pointeurs de CHPOINTs "FACE" des flux aux interfaces:
C
C DT      : pas de temps pour le respect de la CFL-like condition
C           DT &lt; DIAMEL /2 /max(Lambda_i)
C           En maillage regulier cette condition garantie la
C           non-interaction des ondes
C
C LOGNC   : (LOGICAL): si .TRUE. la methode de Newton-Rapson, utilisée
C           dans pour la solution du probleme Riemann n'a pas bien
C           marchéee
C
C LOGAN   : (LOGICAL): si .TRUE. une anomalie à été detectée
C
C MESERR  : pour l'ecriture des messages d'erreurs
C
C************************************************************************
C
C HISTORIQUE (Anomalies et modifications éventuelles)
C
C 07.12.09 Creation
C 21.12.2010 Estension au 3D
C
C************************************************************************
C
C
C N.B.: On suppose qu'on a déjà controllé RO, P > 0
C                                         Y \in (0,1)
C       Si non il faut le faire!!!
C
C************************************************************************
C
C
      IMPLICIT INTEGER(I-N)
      INTEGER I1, I2, ICOMP
     &     , NESP,NORD ,INDMET
     &     , ICHPA1, ICHPA2
     &     , IALF1, IROF1, IVITF1, IPF1
     &     , IALF2, IROF2, IVITF2, IPF2
     &     , IUINF1, IUINF2
     &     , IGRALP, ICHPK0
     &     , ICHPSU, ICHPDI, ICHPVO, MELEMC, MELEMF, MELEFE, MELLIM
     &     , IGEOMC, IGEOMF
     &     , ICHRES
     &     , NFAC
     &     , NLCF, NGCEG,  NGCED, NLCEG, NLCED
     &     , NGCF, NLCF1, SPG1, SPG2, NLFL
     &     , NLHS
      REAL*8 TMAX, RUNIV, K0, K0G, K0D, EPS
      REAL*8 DT, UNSDT, CELLT, CELLTM
     &     , AG1C, ALPG1, ROG1, UNG1, UTG1, UVG1, PG1
     &     , AG2C, ALPG2, ROG2, UNG2, UTG2, UVG2, PG2
     &     , AD1C, ALPD1, ROD1, UND1, UTD1, UVD1, PD1
     &     , AD2C, ALPD2, ROD2, UND2, UTD2, UVD2, PD2
     &     , SURF,CNX, CNY, CNZ, CTX , CTY, CTZ, CVX, CVY, CVZ
     &     , V_INF
     &     , GRALX, GRALY, GRALZ
     &     , CNNA, CTNA, CVNA
     &     , CELL, DIAMG, DIAMD, DIAM, VOLG, VOLD
     &     , FLUCEL, SCON1, SG1, SD1, SCON2, SG2, SD2
     &     , S1_11, S2_22, S1_12, S1_21, S2_12, S2_21
     &     , S1_21l, S1_12l, S2_21l, S2_12l
     &     , S1_21r, S1_12r, S2_21r, S2_12r
     &     , S1_21e, S1_12e, S2_21e, S2_12e
     &     , RESG(12), RESD(12)
     &     , SURFL, SURF0
      LOGICAL LOGNC, LOGAN
      CHARACTER*(40) MESERR
      CHARACTER*(8) TYPE
C
C**** Segment des proprietes du gaz
C
      SEGMENT PROPHY
          REAL*8 XTAB(N1,N2)
C         REAL*8 ACV(NORDP1,NESP+1), W(NESP+1), H0K(NESP+1)
C     &         ,ACVTOG(NORDP1), ACVTOD(NORDP1)
      ENDSEGMENT
C     N1 = NORD + 3
C     N2 = NESP
C
C**** LES INCLUDES
C
 
-INC PPARAM
-INC CCOPTIO
-INC SMCOORD
-INC SMELEME
-INC SMCHAML
C     Normal and tangent
      POINTEUR MELVNX.MELVAL, MELVNY.MELVAL, MELVNZ.MELVAL,
     &     MELT1X.MELVAL, MELT1Y.MELVAL, MELT1Z.MELVAL,
     &     MELT2X.MELVAL, MELT2Y.MELVAL, MELT2Z.MELVAL
      POINTEUR MELUN1.MELVAL, MELUT1.MELVAL, MELUV1.MELVAL
      POINTEUR MELAL1.MELVAL, MELRO1.MELVAL, MELP1.MELVAL
      POINTEUR MELUN2.MELVAL, MELUT2.MELVAL, MELUV2.MELVAL
      POINTEUR MELAL2.MELVAL, MELRO2.MELVAL, MELP2.MELVAL
-INC SMCHPOI
      POINTEUR MPOVSU.MPOVAL, MPOVDI.MPOVAL, MPOVVO.MPOVAL
     &        , MPORES.MPOVAL, MPOALP.MPOVAL
     &        , MPALP1.MPOVAL, MPALP2.MPOVAL
     &        , MPUIN2.MPOVAL, MPUIN1.MPOVAL
     &        , MPOK0.MPOVAL
-INC SMLMOTS
-INC SMLENTI
      POINTEUR MLELIM.MLENTI
      INTEGER JG
-INC SMLREEL
      POINTEUR MLRMFR.MLREEL, MLRCHE.MLREEL
C
C**** Segment for the subroutine prithe
C
      POINTEUR IWST_L.MLREEL, IWST_R.MLREEL, IWWORK.MLREEL
C     1:(4+IDIM+NESP)
      POINTEUR IY.MLREEL
C     1:NESP
      POINTEUR ICVL.MLREEL, ICVR.MLREEL, ICVRSS.MLREEL, ICVWOR.MLREEL
C     0:NORD
C
      POINTEUR IFLU.MLREEL, IFLULA.MLREEL
C
C**** Surface de flamme totale
C
      SURFL = 0.0D0
C
C**** KRIPAD pour la correspondance global/local des conditions limits
C
      CALL KRIPAD(MELLIM,MLELIM)
C      SEGINI MLELIM
C
C**** Initialisation des MCHAMLs
C
C     Alpha
C
      MCHEL1 = IALF1
      SEGACT MCHEL1
      MCHAM1 = MCHEL1.ICHAML(1)
      SEGACT MCHAM1
      MELAL1 = MCHAM1.IELVAL(1)
      SEGDES MCHEL1
      SEGDES MCHAM1
C
C     Masse volumique
C
      MCHEL1 = IROF1
      SEGACT MCHEL1
      MCHAM1 = MCHEL1.ICHAML(1)
      SEGACT MCHAM1
      MELRO1 = MCHAM1.IELVAL(1)
      SEGDES MCHEL1
      SEGDES MCHAM1
C
C     Pression
C
      MCHEL1 = IPF1
      SEGACT MCHEL1
      MCHAM1 = MCHEL1.ICHAML(1)
      SEGACT MCHAM1
      MELP1 = MCHAM1.IELVAL(1)
      SEGDES MCHEL1
      SEGDES MCHAM1
C
C     Vitesse et cosinus directeurs du repere (n,t)
C
      MCHEL1 = IVITF1
      SEGACT MCHEL1
C
C     La vitesse a comme SPG MELEFE
C     Le cosinus directeurs ont comme SPG MELEMF
C
C     MCHAM1 -> Cosinus directeurs
C     MCHAM2 -> Vitesse
C
      SPG1 = MCHEL1.IMACHE(1)
      SPG2 = MCHEL1.IMACHE(2)
      IF((SPG1 .EQ. MELEMF)  .AND. (SPG2 .EQ. MELEFE))THEN
         MCHAM1 = MCHEL1.ICHAML(1)
         MCHAM2 = MCHEL1.ICHAML(2)
      ELSEIF((SPG1 .EQ. MELEFE)  .AND. (SPG2 .EQ. MELEMF))THEN
         MCHAM1 = MCHEL1.ICHAML(2)
         MCHAM2 = MCHEL1.ICHAML(1)
      ELSE
         LOGAN = .TRUE.
         GOTO 9999
      ENDIF
      SEGACT MCHAM1
      SEGACT MCHAM2
      IF (IDIM .EQ. 2) THEN
         MELVNX = MCHAM1.IELVAL(1)
         MELVNY = MCHAM1.IELVAL(2)
         MELT1X = MCHAM1.IELVAL(3)
         MELT1Y = MCHAM1.IELVAL(4)
         MELUN1 = MCHAM2.IELVAL(1)
         MELUT1 = MCHAM2.IELVAL(2)
      ELSE
         MELVNX = MCHAM1.IELVAL(1)
         MELVNY = MCHAM1.IELVAL(2)
         MELVNZ = MCHAM1.IELVAL(3)
         MELT1X = MCHAM1.IELVAL(4)
         MELT1Y = MCHAM1.IELVAL(5)
         MELT1Z = MCHAM1.IELVAL(6)
         MELT2X = MCHAM1.IELVAL(7)
         MELT2Y = MCHAM1.IELVAL(8)
         MELT2Z = MCHAM1.IELVAL(9)
         MELUN1 = MCHAM2.IELVAL(1)
         MELUT1 = MCHAM2.IELVAL(2)
         MELUV1 = MCHAM2.IELVAL(3)
      ENDIF
      SEGDES MCHAM1
      SEGDES MCHAM2
      SEGDES MCHEL1
C
C     Alpha
C
      MCHEL1 = IALF2
      SEGACT MCHEL1
      MCHAM1 = MCHEL1.ICHAML(1)
      SEGACT MCHAM1
      MELAL2 = MCHAM1.IELVAL(1)
      SEGDES MCHEL1
      SEGDES MCHAM1
C
C     Masse volumique
C
      MCHEL1 = IROF2
      SEGACT MCHEL1
      MCHAM1 = MCHEL1.ICHAML(1)
      SEGACT MCHAM1
      MELRO2 = MCHAM1.IELVAL(1)
      SEGDES MCHEL1
      SEGDES MCHAM1
C
C     Pression
C
      MCHEL1 = IPF2
      SEGACT MCHEL1
      MCHAM1 = MCHEL1.ICHAML(1)
      SEGACT MCHAM1
      MELP2 = MCHAM1.IELVAL(1)
      SEGDES MCHEL1
      SEGDES MCHAM1
C
C     Vitesse et cosinus directeurs du repere (n,t)
C
      MCHEL1 = IVITF2
      SEGACT MCHEL1
C
C     La vitesse a comme SPG MELEFE
C     Le cosinus directeurs ont comme SPG MELEMF
C
C     MCHAM1 -> Cosinus directeurs
C     MCHAM2 -> Vitesse
C
      SPG1 = MCHEL1.IMACHE(1)
      SPG2 = MCHEL1.IMACHE(2)
      IF((SPG1 .EQ. MELEMF)  .AND. (SPG2 .EQ. MELEFE))THEN
         MCHAM1 = MCHEL1.ICHAML(1)
         MCHAM2 = MCHEL1.ICHAML(2)
      ELSEIF((SPG1 .EQ. MELEFE)  .AND. (SPG2 .EQ. MELEMF))THEN
         MCHAM1 = MCHEL1.ICHAML(2)
         MCHAM2 = MCHEL1.ICHAML(1)
      ELSE
         LOGAN = .TRUE.
         GOTO 9999
      ENDIF
C      SEGACT MCHAM1
C      MELVNX = MCHAM1.IELVAL(1)
C      MELVNY = MCHAM1.IELVAL(2)
C      MELT1X = MCHAM1.IELVAL(3)
C      MELT1Y = MCHAM1.IELVAL(4)
C      SEGDES MCHAM1
      SEGACT MCHAM2
      MELUN2 = MCHAM2.IELVAL(1)
      MELUT2 = MCHAM2.IELVAL(2)
      IF (IDIM .EQ. 3) MELUV2 = MCHAM2.IELVAL(3)
      SEGDES MCHAM2
      SEGDES MCHEL1
C      write(*,*) 'Son in kodfl1'
C      write(*,*) 'Ho letto gli MCHAM'
C
C**** Segment for fcolre
C
C     Chemical coefficients in MLRCHE.MLREEL.
      SEGACT MLRCHE
      NLHS = 0
      DO I1 = 1, NESP, 1
         IF (MLRCHE.PROG(I1) .GT. 0.0D0) NLHS = NLHS + 1
C         write(*,*) MLRCHE.PROG(I1)
      ENDDO
c      write(*,*) NLHS
C
C     Mass fractions
C     MLRMFR.PROG(1) = YINI_first
C     MLRMFR.PROG(2) = YFIN_first
C     MLRMFR.PROG(3) = YFIN_second
C     ...
C     MLRMFR.PROG(nLHS+1) = YFIN_nLHS
C     MLRMFR.PROG(nLHS+2) = YINI_(nLHS+1)
C
      SEGACT MLRMFR
C
C     cv, molar weicht formation energy in PROPHY.XTAB,
C     already active
C
C     States
      JG = NESP + IDIM + 4
      SEGINI IWWORK
      SEGINI IWST_R
      SEGINI IWST_L
C
C     Mass fractions
      JG = NESP
      SEGINI IY
C
C     Coeff for CV
      JG = NORD + 1
      SEGINI ICVL
      SEGINI ICVR
      SEGINI ICVRSS
      SEGINI ICVWOR
C
C     Fluxes
      JG = NESP + IDIM + 4
      SEGINI IFLU
      SEGINI IFLULA
C
C     Since the mixture here considered is homogeneous, we first fill
C     the part of IWST_L and IWST_R involving the species mass fractions
C     and K0
C
      IWST_L.PROG(4 + IDIM) = MLRMFR.PROG(2)
      IWST_R.PROG(4 + IDIM) = MLRMFR.PROG(2)
      DO I2 = 3, NESP, 1
         I1 = 2 + IDIM + I2
C        DO I1 = 5 + IDIM, 2 + IDIM + NESP
         IWST_L.PROG(I1) = MLRMFR.PROG(I2)
         IWST_R.PROG(I1) = MLRMFR.PROG(I2)
      ENDDO
      IWST_L.PROG(3 + IDIM + NESP) = MLRMFR.PROG(1) - MLRMFR.PROG(2)
      IWST_R.PROG(3 + IDIM + NESP) = MLRMFR.PROG(1) - MLRMFR.PROG(2)
C
C
C**** Initialisation des MELEMEs
C
C     'CENTRE', 'FACEL'
C
      IPT2 = MELEFE
      SEGACT IPT2
      NFAC = IPT2.NUM(/2)
C
C**** KRIPAD pour la correspondance global/local de centre
C
      CALL KRIPAD(MELEMC,MLENT1)
C
C**** MLENTI1 a MCORD.XCOORD(/1)/(IDIM+1) elements
C
C     Si i est le numero global d'un noeud de ICEN,
C     MLENT1.LECT(i) contient sa position, i.e.
C
C     I              = numero global du noeud centre
C     MLENT1.LECT(i) = numero local du noeud centre
C
C     MLENT1 déjà activé, i.e.
C
C     SEGINI MLENT1
C
C
C**** KRIPAD pour la correspondance global/local de 'FACE'
C
      CALL KRIPAD(MELEMF,MLENT2)
C      SEGINI MLENT2
C
C
C**** CHPOINTs de la normale à alpha
C
      CALL LICHT(IGRALP,MPOALP,TYPE,IGEOMF)
      IF (INDK0 .EQ. 2) CALL LICHT(ICHPK0,MPOK0,TYPE,IGEOMF)
C     SEGACT MPOALP
C     SEGDES MPOK0
C
C**** CHPOINTs de la table DOMAINE
C
      CALL LICHT(ICHPSU,MPOVSU,TYPE,IGEOMF)
      CALL LICHT(ICHPDI,MPOVDI,TYPE,IGEOMC)
      CALL LICHT(ICHPVO,MPOVVO,TYPE,IGEOMC)
C
C**** LICHT active les MPOVALs en *MOD
C
C     i.e.
C
C     SEGACT MPOVSU*MOD
C     SEGACT MPOVDI*MOD
C     SEGACT MPOVVO*MOD
C
C
C**** Les FLUX aux face
C
      CALL LICHT(ICHRES,MPORES,TYPE,IGEOMC)
C     SEGACT MPORES*MOD
C
C**** Activation des MCHAMLs
C
      SEGACT MELAL1
      SEGACT MELRO1
      SEGACT MELP1
      SEGACT MELUN1
      SEGACT MELUT1
      SEGACT MELAL2
      SEGACT MELRO2
      SEGACT MELP2
      SEGACT MELUN2
      SEGACT MELUT2
      SEGACT MELVNX
      SEGACT MELVNY
      SEGACT MELT1X
      SEGACT MELT1Y
      IF (IDIM .EQ. 3) THEN
         SEGACT MELUV1
         SEGACT MELUV2
         SEGACT MELVNZ
         SEGACT MELT1Z
         SEGACT MELT2X
         SEGACT MELT2Y
         SEGACT MELT2Z
      ENDIF
C
C**** Lecture des MPOVALs ICHPA1, ICHPA2
C
      CALL LICHT(ICHPA1, MPALP1, TYPE, IGEOMC)
      CALL LICHT(ICHPA2, MPALP2, TYPE, IGEOMC)
C     SEGACT MPALP1
C     SEGACT MPALP2
C
C**** Bas Mach
C
      IF(IUINF1 .NE. 0)THEN
         CALL LICHT(IUINF2,MPUIN2,TYPE,IGEOMC)
         CALL LICHT(IUINF1,MPUIN1,TYPE,IGEOMC)
C        SEGACT MPUIN2*MOD
C        SEGACT MPUIN1*MOD
      ELSE
         MPUIN2=0
         MPUIN1=0
      ENDIF
C
C**** Initialisation de 1/DT et d'autres variables
C
      UNSDT = 0.0D0
      GRALZ = 0.0D0
      CNZ   = 0.0D0
      CTZ   = 0.0D0
      CVX   = 0.0D0
      CVY   = 0.0D0
      CVZ   = 0.0D0
      UVG1  = 0.0D0
      UVD1  = 0.0D0
      UVG2  = 0.0D0
      UVD2  = 0.0D0
C
C**** BOUCLE SUR FACEL pour le calcul du FLUX/RESIDU
C
      DO NLCF = 1, NFAC
C
C******* NLCF  = numero local du centre de facel
C        NGCF  = numero global du centre de facel
C        NLCF1 = numero local du centre de face
C        NGCEG = numero global du centre ELT "gauche"
C        NLCEG = numero local du centre ELT "gauche"
C        NGCED = numero global du centre ELT "droite"
C        NLCED = numero local du centre ELT "droite"
C
         NGCEG = IPT2.NUM(1,NLCF)
         NGCED = IPT2.NUM(3,NLCF)
         NGCF  = IPT2.NUM(2,NLCF)
         NLCF1 = MLENT2.LECT(NGCF)
         NLCEG = MLENT1.LECT(NGCEG)
         NLCED = MLENT1.LECT(NGCED)
         NLFL  = MLELIM.LECT(NGCF)
C
         VOLG = MPOVVO.VPOCHA(NLCEG,1)
         VOLD = MPOVVO.VPOCHA(NLCED,1)
         SURF0 = MPOVSU.VPOCHA(NLCF,1)
C
C******* NLCF != NLCF1 -> l'auteur (MOI) n'a rien compris.
C
         IF(NLCF .NE. NLCF1)THEN
            MESERR = 'Il ne faut pas jouer avec la console.   '
            LOGAN = .TRUE.
            GOTO 9999
         ENDIF
         IF(NLFL .EQ. 0)THEN
            CELLTM = 0.0D0
C
C********** Recuperation des Etats "gauche" et "droite"
C
            AG1C = MPALP1.VPOCHA(NLCEG,1)
            ALPG1 = MELAL1.VELCHE(1,NLCF)
            ROG1 = MELRO1.VELCHE(1,NLCF)
            UNG1 = MELUN1.VELCHE(1,NLCF)
            UTG1 = MELUT1.VELCHE(1,NLCF)
            PG1 = MELP1.VELCHE(1,NLCF)
            AG2C = MPALP2.VPOCHA(NLCEG,1)
            ALPG2 = MELAL2.VELCHE(1,NLCF)
            ROG2 = MELRO2.VELCHE(1,NLCF)
            UNG2 = MELUN2.VELCHE(1,NLCF)
            UTG2 = MELUT2.VELCHE(1,NLCF)
            PG2 = MELP2.VELCHE(1,NLCF)
C
            AD1C = MPALP1.VPOCHA(NLCED,1)
            ALPD1 = MELAL1.VELCHE(3,NLCF)
            ROD1 = MELRO1.VELCHE(3,NLCF)
            UND1 = MELUN1.VELCHE(3,NLCF)
            UTD1 = MELUT1.VELCHE(3,NLCF)
            PD1 = MELP1.VELCHE(3,NLCF)
            AD2C = MPALP2.VPOCHA(NLCED,1)
            ALPD2 = MELAL2.VELCHE(3,NLCF)
            ROD2 = MELRO2.VELCHE(3,NLCF)
            UND2 = MELUN2.VELCHE(3,NLCF)
            UTD2 = MELUT2.VELCHE(3,NLCF)
            PD2 = MELP2.VELCHE(3,NLCF)
C
            CNX = MELVNX.VELCHE(1,NLCF)
            CNY = MELVNY.VELCHE(1,NLCF)
            CTX = MELT1X.VELCHE(1,NLCF)
            CTY = MELT1Y.VELCHE(1,NLCF)
C
            GRALX = MPOALP.VPOCHA(NLCF,1)
            GRALY = MPOALP.VPOCHA(NLCF,2)
C
C            GRALX = CNX
C            GRALY = CNY
C
            IF (IDIM .EQ. 3) THEN
               GRALZ = MPOALP.VPOCHA(NLCF,3)
               CNZ   = MELVNZ.VELCHE(1,NLCF)
               CTZ   = MELT1Z.VELCHE(1,NLCF)
               CVX   = MELT2X.VELCHE(1,NLCF)
               CVY   = MELT2Y.VELCHE(1,NLCF)
               CVZ   = MELT2Z.VELCHE(1,NLCF)
               UVG1 = MELUV1.VELCHE(1,NLCF)
               UVD1 = MELUV1.VELCHE(3,NLCF)
               UVG2 = MELUV2.VELCHE(1,NLCF)
               UVD2 = MELUV2.VELCHE(3,NLCF)
            ENDIF
C
            CNNA = (GRALX * CNX) + (GRALY * CNY) + (GRALZ * CNZ)
            CTNA = (GRALX * CTX) + (GRALY * CTY) + (GRALZ * CTZ)
            CVNA = (GRALX * CVX) + (GRALY * CVY) + (GRALZ * CVZ)
C
C            if ((ngcf .eq. 90) .or. (ngcf .eq. 84) .or. (ngcf .eq. 110))
C     $           then
C               write(*,*) 'ngcf = ', ngcf
C            endif
C
C            write(*,*) cnna
C
C            write(*,*) alpg1, rog1, ung1, utg1, pg1
C            write(*,*) alpg2, rog2, ung2, utg2, pg2
C            write(*,*) alpd1, rod1, und1, utd1, pd1
C            write(*,*) alpd2, rod2, und2, utd2, pd2
C
C
C
C********** Computation of the partition of the surface
C
            S1_11 = min (ALPG1, ALPD1)
            S2_22 = min (ALPG2, ALPD2)
C           S2_22 = 1.0d0 - (max (ALPG1, ALPD1))
            S1_21 = max (0.0d0, (AD1C - AG1C))
            S1_21r = max (0.0d0, (AD1C - ALPD1))
            S1_21e = max (0.0d0, (ALPD1 - ALPG1))
            S1_21l = max (0.0d0, (ALPG1 - AG1C))
            S2_21 = max (0.0d0, (AG2C - AD2C))
            S2_21l = max (0.0d0, (AG2C - ALPG2))
            S2_21e = max (0.0d0, (ALPG2 - ALPD2))
            S2_21r = max (0.0d0, (ALPD2 - AD2C))
C
C         S1_21 = min (S1_21, S2_21)
C         S2_21 = S1_21
            S1_12  = max (0.0d0, (AG1C - AD1C))
            S1_12l = max (0.0d0, (AG1C - ALPG1))
            S1_12e = max (0.0d0, (ALPG1 - ALPD1))
            S1_12r = max (0.0d0, (ALPD1 - AD1C))
            S2_12  = max (0.0d0, (AD2C - AG2C))
            S2_12r = max (0.0d0, (AD2C - ALPD2))
            S2_12e = max (0.0d0, (ALPD2 - ALPG2))
            S2_12l = max (0.0d0, (ALPG2 - AG2C))
C
C            write(*,*) 'NGCF =', NGCF
C            write(*,*) (S1_11 + S2_22 + S1_12 + S1_21)
C            write(*,*) (S1_11 + S2_22 + S2_12 + S2_21)
*
C            write(*,*) (S1_11 + S2_22 + S1_12 + S1_21)
C            write(*,*) (S1_11 + S2_22 + S2_12 + S2_21)
c            CNNA = 1.0D0
c            CTNA = 0.0D0
c            CVNA = 0.0D0
C
C           SURFL = Surface which burns * (CNNA)
*
            SURFL = SURFL + ((ABS ((S1_12 + S1_21) * CNNA)) * SURF0)
C
C           Flame speed
C
            IF (INDK0 .EQ. 1) THEN
               K0G = K0
               K0D = K0
            ELSE
               K0G = MPOK0.VPOCHA(NLCEG,1)
               K0D = MPOK0.VPOCHA(NLCED,1)
            ENDIF
            IWST_L.PROG(4 + IDIM + NESP) = K0G
            IWST_R.PROG(4 + IDIM + NESP) = K0D
C
C********** Riemann problem 1-1
C
            IF (S1_11 .GE. EPS) THEN
               IWST_L.PROG(1) = ROG1
               IWST_L.PROG(2) = UNG1
C               IWST_L.PROG(2) = 0.0D0
               IWST_L.PROG(3) = UTG1
C              IDIM + 2 = 2 in 2D case
               IWST_L.PROG(IDIM + 2)= PG1
               IWST_L.PROG(IDIM + 3) = 0.0D0
C
               IWST_R.PROG(1) = ROD1
               IWST_R.PROG(2) = UND1
C               IWST_R.PROG(2) = 0.0D0
               IWST_R.PROG(3) = UTD1
               IWST_R.PROG(IDIM + 2) = PD1
               IWST_R.PROG(IDIM + 3) = 0.0D0
               IF (IDIM .EQ. 3) THEN
                  IWST_L.PROG(IDIM + 1) = UVG1
                  IWST_R.PROG(IDIM + 1) = UVD1
               ENDIF
C
C               do i2 = 1 , 4 + NESP + IDIM, 1
C                  write(*,*) i2, iwst_l.prog(i2), iwst_r.prog(i2)
C               enddo
C
               IF (INDMET .EQ. 1) THEN
                  CALL FVLHRE(IDIM, NESP, NLHS, NORD,
     &                 MLRCHE.PROG, PROPHY.XTAB, RUNIV,
     &                 TMAX,
     &                 CNNA, CTNA, CVNA,
     &                 IWST_L.PROG,
     &                 IWST_R.PROG,
     &                 IWWORK.PROG, IY.PROG, ICVL.PROG, ICVR.PROG,
     &                 ICVRSS.PROG, ICVWOR.PROG,
     &                 IFLU.PROG, IFLULA.PROG, CELLT,
     &                 LOGAN)
               ELSEIF(INDMET .EQ. 2)THEN
                  CALL FBECRE(IDIM, NESP, NLHS, NORD,
     &                 MLRCHE.PROG, PROPHY.XTAB, RUNIV,
     &                 TMAX,
     &                 CNNA, CTNA, CVNA,
     &                 IWST_L.PROG,
     &                 IWST_R.PROG,
     &                 IWWORK.PROG, IY.PROG, ICVL.PROG, ICVR.PROG,
     &                 ICVRSS.PROG, ICVWOR.PROG,
     &                 IFLU.PROG, IFLULA.PROG, CELLT,
     &                 LOGAN)
               ELSEIF(INDMET .EQ. 3)THEN
C
C********** AUSMP-up
C
                  V_INF=MAX(MPUIN1.VPOCHA(NLCEG,1),
     &                 MPUIN1.VPOCHA(NLCED,1))
                  CALL FAUSRE(IDIM, NESP, NLHS, NORD,
     &                 MLRCHE.PROG, PROPHY.XTAB, RUNIV,
     &                 TMAX,V_INF,
     &                 CNNA, CTNA, CVNA,
     &                 IWST_L.PROG,
     &                 IWST_R.PROG,
     &                 IWWORK.PROG, IY.PROG, ICVL.PROG, ICVR.PROG,
     &                 ICVRSS.PROG, ICVWOR.PROG,
     &                 IFLU.PROG, IFLULA.PROG, CELLT,
     &                 LOGAN)
               ENDIF
 
C                  write(*,*) 'Riemann problem 1-1'
C                  write(*,*) iflula.prog(idim+1), iflu.prog(idim+1)
C
C                subroutine fcolre(ndim, nesp, nLHS, nordpo,
C               &     reacoe, aprop, Runiv,
C               &     Tmaxcv,
C               &     wvect_l,
C               &     wvect_r,
C               &     wwork, y, acv_l, acv_r, acv_rss, acv_work,
C               &     flu, flu_lag, cellt,
C               &     logan)
C
C
               IF (LOGAN) THEN
                  WRITE(IOIMP,*) 'SUBROUTINE KODFL1'
                  WRITE(IOIMP,*) 'RIEMANN PROBLEM 1-1'
                  WRITE(IOIMP,*) 'FACE ', NGCF, ' S11 ', S1_11
                  WRITE(IOIMP,*) '  ICOMP     WVECT_L       WVECT_R'
                  DO I2 = 1, 2 + IDIM
                     WRITE(IOIMP,'(I6,2E16.6)') I2, IWST_L.PROG(I2),
     $                    IWST_R.PROG(I2)
                  ENDDO
                  WRITE(IOIMP,*) 'ANOMALY DETECTED'
                  GOTO 9999
               ENDIF
C
               CELLTM = MAX(CELLT, CELLTM)
C
C************* Ecriture du residu
C
C              FLUX(2) = RO Un              RO Un
C              FLUX(3) = RO Un Un + P   ->  RO Un Ux + P CNX
C              FLUX(4) = RO Un Ut       ->  RO Un Uy + P CNY
C              FLUX(2 + IDIM) = RO Un Et    RO Un Et
C
               SURF   = SURF0 * S1_11
               FLUCEL = IFLU.PROG(1) * SURF
C               write(*,*) 'flux', iflu.prog(1)
C               write(*,*) 'nlceg, nlced', nlceg, nlced
C               write(*,*) 'surf', surf
C
C              Rho un
C
               MPORES.VPOCHA(NLCEG,2) = MPORES.VPOCHA(NLCEG,2) -
     &              (FLUCEL / VOLG)
               IF (NLCED .NE. NLCEG)
     &              MPORES.VPOCHA(NLCED,2) = MPORES.VPOCHA(NLCED,2) +
     &              (FLUCEL / VOLD)
C
C              Rho (Un Ux + P CNX)
C              NB : In 2D (IFLU.PROG(4)*CVX) = RO Un Et * 0 = 0
C                   In 3D (IFLU.PROG(4)*CVX) = RO Un uv * CVX
C
               FLUCEL = ( (IFLU.PROG(2)*CNX) + (IFLU.PROG(3)*CTX) +
     &              (IFLU.PROG(4)*CVX)) * SURF
               MPORES.VPOCHA(NLCEG,3) = MPORES.VPOCHA(NLCEG,3) -
     &              (FLUCEL / VOLG)
               IF (NLCED .NE. NLCEG)
     &              MPORES.VPOCHA(NLCED,3) = MPORES.VPOCHA(NLCED,3) +
     &              (FLUCEL / VOLD)
C
C              Rho (Un Uy + P CNY)
C
               FLUCEL = ( (IFLU.PROG(2)*CNY) + (IFLU.PROG(3)*CTY)
     &              + (IFLU.PROG(4)*CVY) ) * SURF
               MPORES.VPOCHA(NLCEG,4) = MPORES.VPOCHA(NLCEG,4) -
     &              (FLUCEL / VOLG)
               IF (NLCED .NE. NLCEG)
     &              MPORES.VPOCHA(NLCED,4) = MPORES.VPOCHA(NLCED,4) +
     &              (FLUCEL / VOLD)
C
C              Rho (Un Uz + P CNZ)
C
               IF (IDIM .EQ. 3) THEN
                  FLUCEL = ( (IFLU.PROG(2)*CNZ) + (IFLU.PROG(3)*CTZ)
     &              + (IFLU.PROG(4)*CVZ) ) * SURF
                  MPORES.VPOCHA(NLCEG,5) = MPORES.VPOCHA(NLCEG,5) -
     &                 (FLUCEL / VOLG)
                  IF (NLCED .NE. NLCEG)
     &                 MPORES.VPOCHA(NLCED,5) = MPORES.VPOCHA(NLCED,5) +
     &                 (FLUCEL / VOLD)
               ENDIF
C
C              Rho Un et
C
               FLUCEL = IFLU.PROG(2 + IDIM) * SURF
               MPORES.VPOCHA(NLCEG,3 + IDIM) =
     &              MPORES.VPOCHA(NLCEG, 3 + IDIM) -
     &              (FLUCEL / VOLG)
               IF (NLCED .NE. NLCEG)
     &              MPORES.VPOCHA(NLCED, 3 + IDIM) =
     &              MPORES.VPOCHA(NLCED, 3 + IDIM) +
     &              (FLUCEL / VOLD)
C
            ENDIF
C
C********** Riemann problem 2-2
C
            IF (S2_22 .GE. EPS) THEN
               IWST_L.PROG(1) = ROG2
               IWST_L.PROG(2) = UNG2
C               IWST_L.PROG(2) = 0.0D0
               IWST_L.PROG(3) = UTG2
               IWST_L.PROG(IDIM + 2) = PG2
               IWST_L.PROG(IDIM + 3) = 1.0D0
C
               IWST_R.PROG(1) = ROD2
               IWST_R.PROG(2) = UND2
C               IWST_R.PROG(2) = 0.0D0
               IWST_R.PROG(3) = UTD2
               IWST_R.PROG(IDIM + 2) = PD2
               IWST_R.PROG(IDIM + 3) = 1.0D0
               IF (IDIM .EQ. 3) THEN
                  IWST_L.PROG(4) = UVG2
                  IWST_R.PROG(4) = UVD2
               ENDIF
C
C               do i2 = 1 , 4 + NESP + IDIM, 1
C                  write(*,*) i2, iwst_l.prog(i2), iwst_r.prog(i2)
C               enddo
C
               IF (INDMET .EQ. 1) THEN
                  CALL FVLHRE(IDIM, NESP, NLHS, NORD,
     &                 MLRCHE.PROG, PROPHY.XTAB, RUNIV,
     &                 TMAX,
     &                 CNNA, CTNA, CVNA,
     &                 IWST_L.PROG,
     &                 IWST_R.PROG,
     &                 IWWORK.PROG, IY.PROG, ICVL.PROG, ICVR.PROG,
     &                 ICVRSS.PROG, ICVWOR.PROG,
     &                 IFLU.PROG, IFLULA.PROG, CELLT,
     &                 LOGAN)
               ELSEIF(INDMET .EQ. 2)THEN
                  CALL FBECRE(IDIM, NESP, NLHS, NORD,
     &                 MLRCHE.PROG, PROPHY.XTAB, RUNIV,
     &                 TMAX,
     &                 CNNA, CTNA, CVNA,
     &                 IWST_L.PROG,
     &                 IWST_R.PROG,
     &                 IWWORK.PROG, IY.PROG, ICVL.PROG, ICVR.PROG,
     &                 ICVRSS.PROG, ICVWOR.PROG,
     &                 IFLU.PROG, IFLULA.PROG, CELLT,
     &                 LOGAN)
               ELSEIF(INDMET .EQ. 3)THEN
C
C********** AUSMP-up
C
                  V_INF=MAX(MPUIN2.VPOCHA(NLCEG,1),
     &                 MPUIN2.VPOCHA(NLCED,1))
                  CALL FAUSRE(IDIM, NESP, NLHS, NORD,
     &                 MLRCHE.PROG, PROPHY.XTAB, RUNIV,
     &                 TMAX, V_INF,
     &                 CNNA, CTNA, CVNA,
     &                 IWST_L.PROG,
     &                 IWST_R.PROG,
     &                 IWWORK.PROG, IY.PROG, ICVL.PROG, ICVR.PROG,
     &                 ICVRSS.PROG, ICVWOR.PROG,
     &                 IFLU.PROG, IFLULA.PROG, CELLT,
     &                 LOGAN)
               ENDIF
C
C                  write(*,*) 'Riemann problem 2-2'
C                  write(*,*) iflula.prog(idim+1), iflu.prog(idim+1)
C
C                subroutine fcolre(ndim, nesp, nLHS, nordpo,
C               &     reacoe, aprop, Runiv,
C               &     Tmaxcv,
C               &     wvect_l,
C               &     wvect_r,
C               &     wwork, y, acv_l, acv_r, acv_rss, acv_work,
C               &     flu, flu_lag, cellt,
C               &     logan)
C
C
               IF (LOGAN) THEN
                  WRITE(IOIMP,*) 'SUBROUTINE KODFL1'
                  WRITE(IOIMP,*) 'RIEMANN PROBLEM 2-2'
                  WRITE(IOIMP,*) 'FACE ', NGCF, ' S22 ', S2_22
                  WRITE(IOIMP,*) '  ICOMP     WVECT_L       WVECT_R'
                  DO I2 = 1, 2 + IDIM
                     WRITE(IOIMP,'(I6,2E16.6)') I2, IWST_L.PROG(I2),
     $                    IWST_R.PROG(I2)
                  ENDDO
                  WRITE(IOIMP,*) 'ANOMALY DETECTED'
                  GOTO 9999
               ENDIF
C
               CELLTM = MAX(CELLT, CELLTM)
C
C
C************* Ecriture du residu
C
C              FLUX(7) = RO Un              RO Un
C              FLUX(8) = RO Un Un + P   ->  RO Un Ux + P CNX
C              FLUX(9) = RO Un Ut       ->  RO Un Uy + P CNY
C              FLUX(10)= RO Un Et           RO Un Et
C
               SURF   = SURF0 * S2_22
C
C              Rho un
C
               FLUCEL = IFLU.PROG(1) * SURF
               MPORES.VPOCHA(NLCEG, (IDIM + 5)) =
     &              MPORES.VPOCHA(NLCEG, (IDIM + 5)) -
     &              (FLUCEL / VOLG)
               IF (NLCED .NE. NLCEG)
     &              MPORES.VPOCHA(NLCED, (IDIM + 5)) =
     &              MPORES.VPOCHA(NLCED, (IDIM + 5)) +
     &              (FLUCEL / VOLD)
C
C              Rho un ux + P cnx
C
               FLUCEL = ((IFLU.PROG(2)*CNX) + (IFLU.PROG(3)*CTX) +
     &              (IFLU.PROG(4)*CVX)) * SURF
               MPORES.VPOCHA(NLCEG, (IDIM + 6)) =
     &              MPORES.VPOCHA(NLCEG, (IDIM + 6)) -
     &              (FLUCEL / VOLG)
               IF (NLCED .NE. NLCEG)
     &              MPORES.VPOCHA(NLCED, (IDIM + 6)) =
     &              MPORES.VPOCHA(NLCED, (IDIM + 6)) +
     &              (FLUCEL / VOLD)
C
C              Rho un uy + P cny
C
               FLUCEL = ((IFLU.PROG(2)*CNY) + (IFLU.PROG(3)*CTY) +
     &              (IFLU.PROG(4)*CVY)) * SURF
               MPORES.VPOCHA(NLCEG, (IDIM + 7)) =
     &              MPORES.VPOCHA(NLCEG, (IDIM + 7)) -
     &              (FLUCEL / VOLG)
               IF (NLCED .NE. NLCEG)
     &              MPORES.VPOCHA(NLCED, (IDIM + 7)) =
     &              MPORES.VPOCHA(NLCED, (IDIM + 7)) +
     &              (FLUCEL / VOLD)
C
C              Rho un uz + P cnz
C
               IF (IDIM .EQ. 3) THEN
                  FLUCEL = ((IFLU.PROG(2)*CNZ) + (IFLU.PROG(3)*CTZ) +
     &                 (IFLU.PROG(4)*CVZ)) * SURF
                  MPORES.VPOCHA(NLCEG, (IDIM + 8)) =
     &                 MPORES.VPOCHA(NLCEG, (IDIM + 8)) -
     &                 (FLUCEL / VOLG)
                  IF (NLCED .NE. NLCEG)
     &                 MPORES.VPOCHA(NLCED, (IDIM + 8)) =
     &                 MPORES.VPOCHA(NLCED, (IDIM + 8)) +
     &                 (FLUCEL / VOLD)
               ENDIF
C
C              Rho un et
C
               FLUCEL = IFLU.PROG(IDIM + 2) * SURF
               MPORES.VPOCHA(NLCEG, ((2 * IDIM) + 6)) =
     &              MPORES.VPOCHA(NLCEG, ((2 * IDIM) + 6)) -
     &              (FLUCEL / VOLG)
               IF (NLCED .NE. NLCEG)
     &              MPORES.VPOCHA(NLCED, ((2 * IDIM) + 6)) =
     &              MPORES.VPOCHA(NLCED, ((2 * IDIM) + 6)) +
     &              (FLUCEL / VOLD)
            ENDIF
C
C********** Riemann problem 1-2
C           From a numerical point of view, S1_12 and S2_12 are identical.
C           The volume fraction of 1 is larger in l than in r.
C
C           rho, u, p for the phase 1 and 2.
C           NB in case of wall-boundary conditions, we cannot be here
C              since ALPD1 = ALPG1, ALPD2 = ALPG2.
C
            IF ((S1_12 .gt. EPS) .or. (S2_12 .gt. EPS)) THEN
C
               IWST_L.PROG(1) = ROG1
               IWST_L.PROG(2) = UNG1
C               IWST_L.PROG(2) = 0.0D0
               IWST_L.PROG(3) = UTG1
               IWST_L.PROG(IDIM + 2) = PG1
               IWST_L.PROG(IDIM + 3) = 0.0D0
C
               IWST_R.PROG(1) = ROD2
               IWST_R.PROG(2) = UND2
C               IWST_R.PROG(2) = 0.0D0
               IWST_R.PROG(3) = UTD2
               IWST_R.PROG(IDIM + 2) = PD2
               IWST_R.PROG(IDIM + 3) = 1.0D0
               IF (IDIM .EQ. 3) THEN
                  IWST_L.PROG(4) = UVG1
                  IWST_R.PROG(4) = UVD2
               ENDIF
C
C               do i2 = 1 , 4 + NESP + IDIM, 1
C                  write(*,*) i2, iwst_l.prog(i2), iwst_r.prog(i2)
C               enddo
C
               CALL FBECRE(IDIM, NESP, NLHS, NORD,
     &              MLRCHE.PROG, PROPHY.XTAB, RUNIV,
     &              TMAX,
     &              CNNA, CTNA, CVNA,
     &              IWST_L.PROG,
     &              IWST_R.PROG,
     &              IWWORK.PROG, IY.PROG, ICVL.PROG, ICVR.PROG,
     &              ICVRSS.PROG, ICVWOR.PROG,
     &              IFLU.PROG, IFLULA.PROG, CELLT,
     &              LOGAN)
C
C                  write(*,*) 'Riemann problem 1-2'
C                  write(*,*) iflula.prog(idim+1), iflu.prog(idim+1)
C
C                subroutine fcolre(ndim, nesp, nLHS, nordpo,
C               &     reacoe, aprop, Runiv,
C               &     Tmaxcv,
C               &     wvect_l,
C               &     wvect_r,
C               &     wwork, y, acv_l, acv_r, acv_rss, acv_work,
C               &     flu, flu_lag, cellt,
C               &     logan)
C
               IF (LOGAN) THEN
                  WRITE(IOIMP,*) 'SUBROUTINE KODFL1'
                  WRITE(IOIMP,*) 'RIEMANN PROBLEM 1-2'
                  WRITE(IOIMP,*) 'FACE ', NGCF, ' S12 ', S1_12, S2_12
                  WRITE(IOIMP,*) '  ICOMP     WVECT_L       WVECT_R'
                  DO I2 = 1, 2 + IDIM
                     WRITE(IOIMP,'(I6,2E16.6)') I2, IWST_L.PROG(I2),
     $                    IWST_R.PROG(I2)
                  ENDDO
                  WRITE(IOIMP,*) 'ANOMALY DETECTED'
                  GOTO 9999
               ENDIF
C
               CELLTM = MAX(CELLT, CELLTM)
C
C************* Ecriture du residu
C
C              FLUX(7) = RO Un              RO Un
C              FLUX(8) = RO Un Un + P   ->  RO Un Ux + P CNX
C              FLUX(9) = RO Un Ut       ->  RO Un Uy + P CNY
C              FLUX(10)= RO Un Et           RO Un Et
C
C              FLULA(3+IDIM) = - D_resh, with D_resh is the reactive shock
C                                   speed
C
C              We update resu (we are dealing with a reactive
C              Riemann problem)
C
               IF (IFLULA.PROG(IDIM + 3) .LT. 0.0D0) THEN
C
C                 Note that IFLULA.PROG(IDIM + 3) = - D_resh
C
C                              t
C                              |F    / Flag
C                             ---> --->
C                              |   /
C                         1    |  /       2
C                              | /
C                              |/
C                 ------------------------------> x
C                      ifac            ifac+1
C
C                 D_resh > 0
C                 IFLULA.PROG(5) &lt; 0
C
C                 Definition of SCON (conservative surface) and SG and SD
C                 (lagrangian surfaces)
C
c                 NB SCON > 0 if, for a positive conservative flux, gives negative
C                             contribution in ifac and positive contribution
C                             in ifac+1
C                             If follows that SCON \ge 0
C                    SG and SD > 0 if, for a positive lagrangian flux, gives
C                                  negative contribution in ifac and negative
C                                  contributuion in ifac+1
C
 
C                 In first order case :
C                 PHASE 1, conservative contribution in ifac (-F) and ifac+1 (+F)
C                          non-conservative contribution in ifac+1 (-Flag)
C                          negative contribution => SG and SD > 0
C
                  SCON1 = S1_12E
                  SD1 = S1_12E + S1_12R
                  SG1 = S1_12L
C                 In first order case :
C                 PHASE 2, no conservative contribution
C                          non-conservative contribution in ifac+1 (Flag)
C                          positive contribution => SG and SD &lt; 0
C
                  SCON2 = 0.0D0
                  SD2 = -1.0D0 * (S2_12E + S2_12R)
                  SG2 = -1.0D0 * S2_12L
               ELSE
C
C                              t
C                        \Flag | F
C                        ---> --->
C                          \   |
C                      1    \  |       2
C                            \ |
C                             \|
C                 ------------------------------> x
C                      ifac          ifac+1
C
C                 D_resh &lt; 0
C                 IFLULA.PROG(IDIM + 3) > 0
C
C                 In first order case :
C                 PHASE 1, no conservative contribution in ifac and ifac+1
C                          non-conservative contribution in ifac (-Flag)
C                          negative contribution => SG and SD > 0
C
                  SCON1 = 0.0D0
                  SG1 = S1_12E + S1_12L
                  SD1 = S1_12R
C                 In first order case :
C                 PHASE 2, conservative contribution in ifac and ifac+1
C                          non-conservative contribution in ifac (Flag)
C                          positive contribution => SG and SD &lt; 0
C
                  SCON2 = S2_12E
                  SG2 = -1.0D0 * (S2_12L + S2_12E)
                  SD2 = -1.0D0 * S2_12R
               ENDIF
               SURF   = SURF0
               CALL RDEMUP(IDIM, 0,
     &              SCON1, SG1, SD1,
     &              SCON2, SG2, SD2,
     &              SURF,VOLG,VOLD,
     &              CNX,CNY,CNZ,CTX,CTY,CTZ,CVX,CVY,CVZ,
     &              IFLU.PROG,IFLULA.PROG,RESG,RESD)
C               write(*,*) 'RESD'
               DO ICOMP = 1, ((2 * IDIM) + 6)
                  MPORES.VPOCHA(NLCED,ICOMP) =
     &                 MPORES.VPOCHA(NLCED,ICOMP) + RESD(ICOMP)
C                  write(*,*)  ICOMP, RESD(ICOMP)
               ENDDO
C               write(*,*) 'RESG'
               DO ICOMP = 1, ((2 * IDIM) + 6)
                  MPORES.VPOCHA(NLCEG,ICOMP) =
     &                 MPORES.VPOCHA(NLCEG,ICOMP) + RESG(ICOMP)
C                  write(*,*)  ICOMP, RESG(ICOMP)
               ENDDO
            ENDIF
C
C********** Riemann problem 2-1
C           From a numerical point of view, S2_21 and S2_21 are identical.
C           The volume fraction of 2 is larger in l than in r.
C
C           rho, u, p for the phase 1 and 2.
C
C           NB in case of wall-boundary conditions, we cannot be here
C              since ALPD1 = ALPG1, ALPD2 = ALPG2.
C
            IF ((S1_21 .gt. EPS) .or. (S2_21 .gt. EPS)) THEN
C
C
               IWST_L.PROG(1) = ROG2
               IWST_L.PROG(2) = UNG2
CC               IWST_L.PROG(2) = 0.0D0
               IWST_L.PROG(3) = UTG2
               IWST_L.PROG(IDIM + 2) = PG2
               IWST_L.PROG(IDIM + 3) = 1.0D0
CC
               IWST_R.PROG(1) = ROD1
               IWST_R.PROG(2) = UND1
C               IWST_R.PROG(2) = 0.0D0
               IWST_R.PROG(3) = UTD1
               IWST_R.PROG(IDIM + 2) = PD1
               IWST_R.PROG(IDIM + 3) = 0.0D0
C
               IF (IDIM .EQ. 3) THEN
                  IWST_L.PROG(4) = UVG2
                  IWST_R.PROG(4) = UVD1
               ENDIF
C
CC
C               do i2 = 1 , 4 + NESP + IDIM, 1
C                  write(*,*) i2, iwst_l.prog(i2), iwst_r.prog(i2)
C               enddo
C
               CALL FBECRE(IDIM, NESP, NLHS, NORD,
     &              MLRCHE.PROG, PROPHY.XTAB, RUNIV,
     &              TMAX,
     &              CNNA, CTNA, CVNA,
     &              IWST_L.PROG,
     &              IWST_R.PROG,
     &              IWWORK.PROG, IY.PROG, ICVL.PROG, ICVR.PROG,
     &              ICVRSS.PROG, ICVWOR.PROG,
     &              IFLU.PROG, IFLULA.PROG, CELLT,
     &              LOGAN)
C                  write(*,*) 'Riemann problem 1-2'
C                  write(*,*) iflula.prog(idim+1), iflu.prog(idim+1)
C
C                subroutine fcolre(ndim, nesp, nLHS, nordpo,
C               &     reacoe, aprop, Runiv,
C               &     Tmaxcv,
C               &     wvect_l,
C               &     wvect_r,
C               &     wwork, y, acv_l, acv_r, acv_rss, acv_work,
C               &     flu, flu_lag, cellt,
C               &     logan)
C
               IF (LOGAN) THEN
                  WRITE(IOIMP,*) 'SUBROUTINE KODFL1'
                  WRITE(IOIMP,*) 'RIEMANN PROBLEM 2-1'
                  WRITE(IOIMP,*) 'FACE ', NGCF, ' S21 ', S1_21, S2_21
                  WRITE(IOIMP,*) '  ICOMP     WVECT_L       WVECT_R'
                  DO I2 = 1, 2 + IDIM
                     WRITE(IOIMP,'(I6,2E16.6)') I2, IWST_L.PROG(I2),
     $                    IWST_R.PROG(I2)
                  ENDDO
                  WRITE(IOIMP,*) 'ANOMALY DETECTED'
                  GOTO 9999
               ENDIF
C
               CELLTM = MAX(CELLT, CELLTM)
C
C************* Ecriture du residu
C
C              FLUX(7) = RO Un              RO Un
C              FLUX(8) = RO Un Un + P   ->  RO Un Ux + P CNX
C              FLUX(9) = RO Un Ut       ->  RO Un Uy + P CNY
C              FLUX(10)= RO Un Et           RO Un Et
C
C              FLULA(3+IDIM) = - D_resh, with D_resh is the reactive shock
C                                   speed
C
C
C              We update resu (we are dealing with a reactive
C              Riemann problem)
C
               IF (IFLULA.PROG(IDIM + 3) .LT. 0.0D0) THEN
C
C                 Note that IFLULA.PROG(IDIM + 3) = - D_resh
C
C                              t
C                              |F    / Flag
C                             ---> --->
C                              |   /
C                         2    |  /       1
C                              | /
C                              |/
C                 ------------------------------> x
C                      ifac            ifac+1
C
C                 D_resh > 0
C                 IFLULA.PROG(IDIM + 3) &lt; 0
C
C                 In first order case :
C                 PHASE 2, conservative contribution in ifac (-F) and ifac+1 (+F)
C                          non-conservative contribution in ifac+1 (-Flag)
C                          negative contribution => SG and SD > 0
C
                  SCON2 = S2_21E
                  SG2 = S2_21L
                  SD2 = S2_21E + S2_21R
C
C                 PHASE 1, non-conservative contribution in ifac+1 (+Flag)
C                          positive contribution => SG and SD &lt; 0
C
                  SCON1 = 0.0D0
                  SG1 = -1.0D0 * S1_21L
                  SD1 = -1.0D0 * (S1_21E + S1_21R)
C
               ELSE
C
C                              t
C                        \Flag | F
C                        ---> --->
C                          \   |
C                      2    \  |       1
C                            \ |
C                             \|
C                 ------------------------------> x
C                      ifac          ifac+1
C
C                 D_resh &lt; 0
C                 IFLULA.PROG(IDIM + 3) > 0
C
C                 In first order case :
C                 PHASE 2, no conservative contribution
C                          non-conservative contribution in ifac (-Flag)
C                          negative contribution => SG and SD > 0
C
                  SCON2 = 0.0D0
                  SG2 = S2_21E + S2_21L
                  SD2 = S2_21R
C
C                 In first order case :
C                 PHASE 1, conservative contribution in ifac (-F) and
C                          ifac+1 (+F)
C                          non-conservative contribution in ifac (+Flag)
C                          positive contribution => SG and SD > 0
C
C                 ALPHA1
C
                  SCON1 = S1_21E
                  SG1 = -1.0D0 * (S1_21E + S1_21L)
                  SD1 = -1.0D0 * S1_21R
               ENDIF
               SURF   = SURF0
               CALL RDEMUP(IDIM, 0,
     &              SCON1, SG1, SD1,
     &              SCON2, SG2, SD2,
     &              SURF,VOLG,VOLD,
     &              CNX,CNY,CNZ,CTX,CTY,CTZ,CVX,CVY,CVZ,
     &              IFLU.PROG,IFLULA.PROG,RESG,RESD)
C               write(*,*) 'RESD'
               DO ICOMP = 1, ((2 * IDIM) + 6)
                  MPORES.VPOCHA(NLCED,ICOMP) =
     &                 MPORES.VPOCHA(NLCED,ICOMP) + RESD(ICOMP)
C                  write(*,*)  ICOMP, RESD(ICOMP)
               ENDDO
C               write(*,*) 'RESG'
               DO ICOMP = 1, ((2 * IDIM) + 6)
                  MPORES.VPOCHA(NLCEG,ICOMP) =
     &                 MPORES.VPOCHA(NLCEG,ICOMP) + RESG(ICOMP)
C                  write(*,*)  ICOMP, RESG(ICOMP)
               ENDDO
            ENDIF
C
C********** Calcul du pas du temps (CFL)
C
C*********  a) etat a l'interface
C
            DIAMG = MPOVDI.VPOCHA(NLCEG,1)
            DIAMD = MPOVDI.VPOCHA(NLCED,1)
            DIAM = MIN(DIAMG,DIAMD)
            CELL  = CELLTM/DIAM
            IF(CELL .GT. UNSDT)THEN
               UNSDT  =  CELL
            ENDIF
         ENDIF
C
C**** Fin boucle sur FACEL
C
      ENDDO
C
C**** Pas du temps (condition de non interaction en 1D)
C
      DT = 0.5D0 / UNSDT
C
C**** Desactivation des segments
C
      SEGSUP MLELIM
      SEGSUP MLENT1
      SEGSUP MLENT2
C
      SEGDES MPOVSU
      SEGDES MPOVDI
      SEGDES MPOVVO
      SEGDES MPORES
      SEGDES IPT2
C
C     Warning. MLRMFR = PGAS . 'MASSFRA'
C              MLRCHE = PGAS . 'CHEMCOEF'
C
      SEGDES MLRMFR
      SEGDES MLRCHE
      SEGSUP IY
C
      SEGSUP IWWORK
      SEGSUP IWST_R
      SEGSUP IWST_L
      SEGSUP ICVL
      SEGSUP ICVR
      SEGSUP ICVRSS
      SEGSUP ICVWOR
      SEGSUP IFLU
      SEGSUP IFLULA
C
      SEGDES MPALP1
      SEGDES MPALP2
      SEGDES MELAL1
      SEGDES MELRO1
      SEGDES MELP1
      SEGDES MELAL2
      SEGDES MELRO2
      SEGDES MELP2
      SEGDES MELUN2
      SEGDES MELUT2
      IF (IDIM .EQ. 2) THEN
         SEGDES MELUN1
         SEGDES MELUT1
         SEGDES MELVNX
         SEGDES MELVNY
         SEGDES MELT1X
         SEGDES MELT1Y
         SEGDES MELUN2
         SEGDES MELUT2
      ELSE
         SEGDES   MELVNX
         SEGDES   MELVNY
         SEGDES   MELVNZ
         SEGDES   MELT1X
         SEGDES   MELT1Y
         SEGDES   MELT1Z
         SEGDES   MELT2X
         SEGDES   MELT2Y
         SEGDES   MELT2Z
         SEGDES   MELUN1
         SEGDES   MELUT1
         SEGDES   MELUV1
         SEGDES   MELUN2
         SEGDES   MELUT2
         SEGDES   MELUV2
      ENDIF
C
      SEGDES MPOALP
      IF (INDK0 .EQ. 2) SEGDES MPOK0
C
      IF(MPUIN2 .NE. 0)THEN
         SEGDES MPUIN2
         SEGDES MPUIN1
      ENDIF
C
 9999 CONTINUE
      RETURN
      END
CC