C FLACR3    SOURCE    OF166741  24/12/13    21:15:52     12097          
      SUBROUTINE FLACR3(EPSCSI,EPSI,DELTAT,MCEN,MFACEL,IRC,IYC,IYINIT,
     &     IYFINA,IVCAR,IDX,MLRMAS,MLRH0K,MLRECO,ICHRET,ICHRY)
C************************************************************************
C
C PROJET            :  CASTEM 2000
C
C NOM               :  FLACR3
C
C DESCRIPTION       :  CREBCOM: modele non-homogene
C                      voir FLACR2
C
C LANGAGE           :  FORTRAN 77 + ESOPE 2000 (avec extensions CISI)
C
C AUTEUR            :  A. BECCANTINI, DM2S/SFME/LTMF
C
C************************************************************************
C
C INPUTS
C
C EPSCSI   : (REAL*8), parametre pour controler si CSI=DY/(YF - YI)
C            \in(0,1). Si  CSI<EPSCSI ou CSI>1+EPSCSI un message de
C            warnin est donné.
C
C EPSI     : parametre du critère CREBCOM (REAL*8)
C
C DELTAT   : pas de temps (REAL*8)
C
C MCEN    : pointeur du MELEME contenant les centres des ELTs
C
C MFACEL  : pointeur du MELEME contenant la correspondence
C           CENTRE-FACE-CENTRE
C
C IRC     : pointeur du CHPOINT  contenant la masse volumique.
C
C IYC     : pointeur du CHPOINT contenant les fractions massiques
C
C IYINIT  : pointeur du CHPOINT contenant la fraction massique initiale
C           de la premiere composante de IYC;
C
C IYFINA  : pointeur du CHPOINT contenant la fraction massique finale de
C           la premiere composante de IYC ;
C
C IVCAR   : pointeur du CHPOINT contenant la vitesse caractéristique de
C           la flamme
C
C IDX     : pointeur du CHPOINT contenant la diemnsion de la maille
C
C MLRMAS  : pointeur du LISTREEL contenant les masses molaires
C
C MLRH0K  : pointeur du LISTREEL contenant les energies des formation à
C           0K
C
C MLRCOE  : pointeur du LISTREEL contenant les coeff. stoch. de la
C           reaction chimique
C
C OUTPUTS :
C
C ICHRET  : pointeur de l'increment de l'energie totale par unité de
C           volume
C
C ICHRY   : pointeur de l'increment des densités massiques
C
C************************************************************************
C
C HISTORIQUE (Anomalies et modifications éventuelles)
C
C HISTORIQUE :
C
C
C************************************************************************
C
      IMPLICIT INTEGER(I-N)

-INC PPARAM
-INC CCOPTIO
-INC SMLREEL
      POINTEUR MLRECO.MLREEL, MLRMAS.MLREEL, MLRH0K.MLREEL
-INC SMELEME
      POINTEUR MELCEN.MELEME, MELEFL.MELEME
-INC SMCHPOI
      INTEGER N, NC
      POINTEUR MPOVRO.MPOVAL, MPOVAY.MPOVAL,MPOVYI.MPOVAL,
     &     MPOVYF.MPOVAL, MPOVC.MPOVAL, MPODRE.MPOVAL, MPODRY.MPOVAL,
     &     MPOCSI.MPOVAL,MPOCRI.MPOVAL,MPOVDX.MPOVAL
-INC SMLENTI
      POINTEUR MLECEN.MLENTI
C
C**** Les variables
C
      INTEGER MCEN,MFACEL,IRC,IYC,IYINIT,IYFINA,IVCAR,ICHRET,ICHRY
     $     ,IGEOM,ICEN,NCEN,NFAC,NLCF,NGCEG,NGCED,NLCEG,NLCED
     &     ,IESP,NESP,IDX
      REAL*8 EPSI, DELTAT, Y1F, Y1I, Y1, VCSIG, VCSID, VCSI2G, VCSI2D
     &     ,EPS12, DY, DYMAX, DCSI, DCSI1, DY1, RHO, DY2, YMAX, DYMAX1
     &     ,ERRTOL, EPSCSI, CSIG
      CHARACTER*8 TYPE
      PARAMETER(ERRTOL=1.0D-6)
C
      NESP=MLRH0K.PROG(/1)
      MELCEN=MCEN
      MELEFL=MFACEL
C
C**** KRIPAD pour la correspondance global/local de centre
C
      CALL KRIPAD(MELCEN,MLECEN)
      IF(IERR .NE. 0)GOTO 9999
C     SEGINI MLECEN
      SEGACT MELCEN
      NCEN=MELCEN.NUM(/2)
      SEGACT MELEFL
      NFAC=MELEFL.NUM(/2)
C
C**** Lectures de CHPOINT
C
      CALL LICHT(IRC,MPOVRO,TYPE,IGEOM)
C     SEGACT MPOVRO*MOD
      CALL LICHT(IYC,MPOVAY,TYPE,IGEOM)
C     SEGACT MPOVAY*MOD
      CALL LICHT(IYINIT,MPOVYI,TYPE,IGEOM)
C     SEGACT MPOVYI*MOD
      CALL LICHT(IYFINA,MPOVYF,TYPE,IGEOM)
C     SEGACT MPOVYF*MOD
      CALL LICHT(IVCAR,MPOVC,TYPE,IGEOM)
C     SEGACT MPOVYF*MOD
      CALL LICHT(IDX,MPOVDX,TYPE,IGEOM)
C     SEGACT MPOVDX*MOD
      CALL LICHT(ICHRET,MPODRE,TYPE,IGEOM)
C     SEGACT MPODRE*MOD
      CALL LICHT(ICHRY,MPODRY,TYPE,IGEOM)
C     SEGACT MPODRY*MOD
C
C**** Creation du MPOVAL qui contient csi
C     Creation du MPOVAL du critere
C
      N=NCEN
      NC=1
      SEGINI MPOCSI
      SEGINI MPOCRI
C
C**** Calcul de DYMAX pour la premiere espece
C     Controle de la densité positive
C     Controle des fractions massiques (Yi>0, sum_i Y_i < 1)
C
      DYMAX = 0.0D0
      YMAX=0.0D0
      DO ICEN=1,NCEN,1
         Y1F = MPOVYF.VPOCHA(ICEN,1)
         Y1I = MPOVYI.VPOCHA(ICEN,1)
         DY = Y1F - Y1I
         DYMAX=MAX(ABS(DY),DYMAX)
         Y1=MPOVAY.VPOCHA(ICEN,1)
         IF(Y1 .LT. 0.0D0)THEN
            WRITE(IOIMP,*) 'CHPO2 < 0 ???'
            CALL ERREUR(21)
            GOTO 9999
         ELSE
            YMAX=Y1
         ENDIF
         RHO=MPOVRO.VPOCHA(ICEN,1)
         IF(RHO .LT. 0)THEN
            WRITE(IOIMP,*) 'CHPO1 < 0 ???'
            CALL ERREUR(21)
            GOTO 9999
         ENDIF
         RHO=MPOVYI.VPOCHA(ICEN,1)
         IF((RHO .LT. 0) .OR. (RHO .GT. 1))THEN
            WRITE(IOIMP,*) 'CHPO3 < 0 ou CHPO3 > 1???'
            CALL ERREUR(21)
            GOTO 9999
         ENDIF
         RHO=MPOVYF.VPOCHA(ICEN,1)
         IF((RHO .LT. 0.0D0) .OR. (RHO .GT. 1.0D0))THEN
            WRITE(IOIMP,*) 'CHPO4 < 0 ou CHPO4 > 1???'
            CALL ERREUR(21)
            GOTO 9999
         ENDIF
         RHO=MPOVC.VPOCHA(ICEN,1)
         IF((RHO .LT. 0))THEN
            WRITE(IOIMP,*) 'CHPO5 < 0'
            CALL ERREUR(21)
            GOTO 9999
         ENDIF
         DO IESP=2,NESP,1
            Y1=MPOVAY.VPOCHA(ICEN,IESP)
            IF(Y1 .LT. 0.0D0)THEN
               WRITE(IOIMP,*) 'CHPO2 < 0 ???'
               CALL ERREUR(21)
               GOTO 9999
            ELSE
               YMAX=YMAX+Y1
            ENDIF
         ENDDO
         IF(YMAX .GT. 1.0D0)THEN
            WRITE(IOIMP,*) 'sum CHPO2 > 1 ???'
            CALL ERREUR(21)
            GOTO 9999
         ENDIF
      ENDDO
C
C**** Calcul de CSI \in (0,1)
C
      DO ICEN=1,NCEN,1
         Y1F = MPOVYF.VPOCHA(ICEN,1)
         Y1I = MPOVYI.VPOCHA(ICEN,1)
         DY = Y1F - Y1I
         IF(ABS(DY) .LE. (ERRTOL*DYMAX))THEN
C
C********** Pas de combustion en cette region
C
            MPOCSI.VPOCHA(ICEN,1)=0.0D0
C
         ELSE
            Y1 = MPOVAY.VPOCHA(ICEN,1)
            VCSIG = (Y1 - Y1I) / DY
C
C********** On n'accepte pas csi > 1.1 or csi < -0.1
C
            IF((VCSIG .GT. (1.0D0+EPSCSI)) .OR.
     &           (VCSIG .LT. (-1*EPSCSI)))THEN
               WRITE(IOIMP,*) 'Progress variable = ???'
C              21 2
C              Données incompatibles
               CALL ERREUR(21)
               GOTO 9999
            ELSEIF(VCSIG .GT. 1.0D0)THEN
               VCSIG = 1.0D0
            ELSEIF(VCSIG .LT. 0.0D0)THEN
               VCSIG = 0.0D0
            ENDIF
            MPOCSI.VPOCHA(ICEN,1)= VCSIG
         ENDIF
      ENDDO
C
C**** Le critere CREBCOM
C
      DO NLCF = 1, NFAC, 1
C
C******* NLCF  = numero local du centre de facel
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 = MELEFL.NUM(1,NLCF)
         NGCED = MELEFL.NUM(3,NLCF)
         NLCEG = MLECEN.LECT(NGCEG)
         NLCED = MLECEN.LECT(NGCED)
C
         VCSIG=MPOCSI.VPOCHA(NLCEG,1)
         VCSID=MPOCSI.VPOCHA(NLCED,1)
         VCSI2G=VCSIG*VCSIG
         VCSI2D=VCSID*VCSID
C
         IF(NLCEG .EQ. NLCED)THEN
C
C********** Murs
C
            MPOCRI.VPOCHA(NLCEG,1)=MPOCRI.VPOCHA(NLCEG,1) + (0.5D0 *
     &           VCSI2D)
C
         ELSE
C
            MPOCRI.VPOCHA(NLCEG,1)=MPOCRI.VPOCHA(NLCEG,1) +
     &           (VCSI2D - (0.5D0 * VCSI2G))
            MPOCRI.VPOCHA(NLCED,1)=MPOCRI.VPOCHA(NLCED,1) +
     &           (VCSI2G - (0.5D0 * VCSI2D))
C
         ENDIF
      ENDDO
C
C**** Calcul des increments de l'energie et des densités
C
      EPS12 = EPSI * EPSI
      DO ICEN = 1, NCEN, 1
         CSIG = MPOCSI.VPOCHA(ICEN,1)
         VCSIG =  MPOCRI.VPOCHA(ICEN,1)
C
C******* In 2D, contribution of the ideal upper and lower cells
C
         IF(IDIM .EQ. 2) VCSIG = VCSIG + (CSIG * CSIG)
         IF(VCSIG .GE. (EPS12*(1.0D0+ERRTOL)))THEN
C
C********** Il y a combustion
C
            DCSI=MPOVC.VPOCHA(ICEN,1)*DELTAT/MPOVDX.VPOCHA(ICEN,1)
            DCSI1=1.0D0-CSIG
            DCSI=MIN(DCSI,DCSI1)
            DY1=(MPOVYF.VPOCHA(ICEN,1)-MPOVYI.VPOCHA(ICEN,1))*DCSI
C
C********** On force la positivité de MPOVAY.VPOCHA(ICEN,1)
C
            IF(DY1 .GT. 0.0D0)THEN
               DY2 = 1.0D0 - MPOVAY.VPOCHA(ICEN,1)
            ELSE
               DY2 = MPOVAY.VPOCHA(ICEN,1)
            ENDIF
            DY1=SIGN(MIN(ABS(DY2),ABS(DY1)),DY1)
C
            RHO=MPOVRO.VPOCHA(ICEN,1)
            MPODRY.VPOCHA(ICEN,1)=DY1*RHO
            MPODRE.VPOCHA(ICEN,1)=-1.0D0*DY1*RHO*MLRH0K.PROG(1)
            DY1=DY1/(MLRMAS.PROG(1)*MLRECO.PROG(1))
            DO IESP=2,NESP,1
               DY=DY1*(MLRMAS.PROG(IESP)*MLRECO.PROG(IESP))
               DYMAX1=ABS(DYMAX/(MLRMAS.PROG(1)*MLRECO.PROG(1))*
     &              (MLRMAS.PROG(IESP)*MLRECO.PROG(IESP)))
               IF(DY .GT. 0)THEN
                  DY2=1.0D0 - MPOVAY.VPOCHA(ICEN,IESP)
               ELSE
                  DY2=MPOVAY.VPOCHA(ICEN,IESP)
               ENDIF
C
C************* On force la positivité de MPOVAY.VPOCHA(ICEN,IESP)
C
               IF((ABS(DY)-ABS(DY2)).GE.(0.5D0*DYMAX1))THEN
                  WRITE(IOIMP,*) 'CHPO2, CHPO3, CHPO4 = ???'
C                 Données incompatibles
                  CALL ERREUR(21)
                  GOTO 9999
               ELSE
                  DY=SIGN(MIN(ABS(DY2),ABS(DY)),DY)
               ENDIF
               MPODRY.VPOCHA(ICEN,IESP)=DY*RHO
               MPODRE.VPOCHA(ICEN,1)=MPODRE.VPOCHA(ICEN,1)-(DY*RHO
     $              *MLRH0K.PROG(IESP))
            ENDDO
         ENDIF
      ENDDO
C
      SEGDES MELCEN
      SEGDES MELEFL
      SEGSUP MLECEN
C
      SEGDES MPOVRO
      SEGDES MPOVAY
      SEGDES MPOVYI
      SEGDES MPOVYF
      SEGDES MPOVC
      SEGDES MPOVDX
      SEGDES MPODRE
      SEGDES MPODRY
      SEGSUP MPOCSI
      SEGSUP MPOCRI
C
 9999 RETURN
      END