C KFM12     SOURCE    OF166741  24/12/13    21:16:05     12097          
      SUBROUTINE KFM12(IRES,IRN,IGN,IRETN,IGAMN,IVCO,
     &     ICHPSU,ICHPDI,ICHPVO,INORM,
     &     MELEMC,MELEMF,MELEFL,MELTFA,DTPS,DCFL,
     &     MELLIM,ICHLIM,NJAC,ICOEF,ICOEV,
     &     IDU,IPROBL)
C************************************************************************
C
C PROJET            :  CASTEM 2000
C
C NOM               :  KFM12
C
C DESCRIPTION       :  Voir aussi KFM1
C
C                      Cas deux dimensions, gaz "calorically perfect"
C                      Méthode de relaxation de type SGS
C
C LANGAGE           :  FORTRAN 77 + ESOPE 2000 (avec estensions CISI)
C
C AUTEUR            :  T. KLOCZKO DEN/DM2S/SFME/LTMF
C
C************************************************************************
C ENTREES
C
C
C  1) Pointeurs des CHPOINTs
C
C     IRES    : CHPOINT 'CENTRE' contenant le residu
C
C     IRN     : CHPOINT 'CENTRE' contenant la masse volumique
C
C     IGN     : CHPOINT 'CENTRE' contenant la q.d.m.
C
C     IRETN   : CHPOINT 'CENTRE' contenant l'energie totale
C
C     IGAMN   : CHPOINT 'CENTRE' contenant le gamma
C
C     IVCO    : CHPOINT 'CENTRE' contenant le cut-off
C
C     ICHLIM  : CHPOINT conditions aux bords
C
C     ICOEV   : CHPOINT 'CENTRE' contenant le (gamma - 1) K / (R \rho)
C               coeff pour calculer le ray. spect. de la matrice
C               visc.
C
C  3) Pointeurs de CHPOINTs de la table DOMAINE
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     INORM   : CHPOINT "CENTRE" contenant le normales aux faces
C
C
C  4) 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     MELEFL  : MELEME 'FACEL' du connectivité FACES -> ELEM
C
C     MELLIM  : MELEME SPG conditions aux bords
C
C  5)
C
C     DCFL = le double de la CFL
C
C     DTPS = pas de temps physique
C
C     NJAC = nombre d'iterations dans le jacobien
C
C
C SORTIES
C
C     IDU  : resultat
C
C     IPROBL : si > 0, il y a un probleme
C
C************************************************************************
C
C HISTORIQUE (Anomalies et modifications éventuelles)
C
C HISTORIQUE : 2004  : creation
C
C
C************************************************************************
C
C
C N.B.: On suppose qu'on a déjà controllé RO, P > 0
C                                         GAMMA \in (1,3)
C                                         Y \in (0,1)
C       Si non il faut le faire!!!
C
C************************************************************************
C
      IMPLICIT INTEGER(I-N)
      INTEGER IRN,IGN,IRETN,IGAMN,ICHPSU,ICHPDI,ICHPVO,INORM
     &     ,MELEMF,IDU,NFAC, NFACE
     &     ,IGEOMF,IGEOMC, NGCEG, NGCED, NGCF, NLCF, NLCEG, NLCED
     &     ,ICEN,NCEN,ICHLIM,MELLIM,NLLIM,IRES,NJAC,ICOMP
     &     ,I1,I2, IBLJAC, IPROBL, IVCO, MELTFA
     &     ,NSOUPO,IFACT,NELT, ISOUP, IELT, IFAC, IFACG, JG
     &     ,ISOUP1, ISOUP2, ICOEF,ICOEV
     &     , ID, IELT1, IELT2, IFAC1, IFAC2
      REAL*8 ROG, RUXG, RUYG, UNG, RETG, GAMG, RECG, RECD, PG, VOLG
     &     , DIAMG, ROD, RUXD, RUYD, UND, RETD, GAMD, PD
     &     , CNX, CNY, DCFL, DTPS
     &     , SURF, SIGMAF
     &     , UNSDT, UNSDTF, UXD, UYD
     &     , FR, FRUX, FRUY, FRET, COEF
     &     , CTX, CTY, RHTD, UTD, UTG, UXG, UYG, QG, CG, UCOG
     &     , ROF, PF, UNF,  GAMF, UCOF, CF, PHI2, COEF1
     &     , QN
C
      CHARACTER*8 TYPE
      LOGICAL LOGINV, LOGJAC
C
C**** LES INCLUDES
C

-INC PPARAM
-INC CCOPTIO
-INC SMCHPOI
      POINTEUR MPOVSU.MPOVAL, MPOVDI.MPOVAL
     &       , MPODU.MPOVAL, MPOVOL.MPOVAL
     &       , MPRN.MPOVAL, MPGN.MPOVAL, MPRETN.MPOVAL, MPGAMN.MPOVAL
     &       , MPNORM.MPOVAL, MPLIM.MPOVAL, MPORES.MPOVAL, MPOCO.MPOVAL
     &       , MPOCO1.MPOVAL, MELEFL.MELEME, MELEMC.MELEME
     &       , MPCOEV.MPOVAL
-INC SMCOORD
-INC SMELEME
-INC SMLMOTS
-INC SMLENTI
      POINTEUR MLELIM.MLENTI, MLESUP.MLENTI
C
C**** POINTEUR vecteurs
C
      SEGMENT VEC
      REAL*8 VALVEC(I1)
      ENDSEGMENT
      POINTEUR SIGMA0.VEC, AN0.VEC, DN0.VEC, D1.VEC
     &     , USDTI.VEC, PHI2V.VEC, SIGMAV.VEC
C
C**** POINTEUR matriciels
C
      SEGMENT MAT
      REAL*8 VAL(I1,I2)
      ENDSEGMENT
      POINTEUR CONS0.MAT, CONS1.MAT, BN0.MAT, BN1.MAT
     &     , CN0.MAT, CN1.MAT, DU.MAT
     &     , RES.MAT
     &     , Q1.MAT, Q2.MAT
C
C
C**** IPROBL: si different de 0, un probleme a été detecté
C
      IPROBL=0
C
C**** Initialisation des MLENTI des conditions aux limites
C
C
      CALL KRIPAD(MELLIM,MLELIM)
C     SEGINI MLELIM
C
C**** Initialisation des MELEMEs
C
C     'CENTRE', 'FACEL'
C
      SEGACT MELEFL
      NFAC = MELEFL.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 table DOMAINE
C
      CALL LICHT(ICHPSU,MPOVSU,TYPE,IGEOMF)
      CALL LICHT(ICHPDI,MPOVDI,TYPE,IGEOMC)
      CALL LICHT(ICHPVO,MPOVOL,TYPE,IGEOMC)
      CALL LICHT(INORM,MPNORM,TYPE,IGEOMC)
C
C**** LICHT active les MPOVALs en *MOD
C
C     i.e.
C
C     SEGACT MPOVSU*MOD
C     SEGACT MPOVOL*MOD
C     SEGACT MPOVDI*MOD
C     SEGACT MPNORM*MOD
C
      CALL LICHT(IDU,MPODU,TYPE,IGEOMC)
C     SEGACT MPODU*MOD
C
C     USDTI initialisé a zero; USDTI = 1 / DT
C
      NCEN=MPOVOL.VPOCHA(/1)
      I1=NCEN
      SEGINI USDTI
C
      CALL LICHT(IRES,MPORES,TYPE,IGEOMC)
      CALL LICHT(IRN,MPRN,TYPE,IGEOMC)
      CALL LICHT(IGN,MPGN,TYPE,IGEOMC)
      CALL LICHT(IRETN,MPRETN,TYPE,IGEOMC)
      CALL LICHT(IGAMN,MPGAMN,TYPE,IGEOMC)
      CALL LICHT(IVCO,MPOCO,TYPE,IGEOMC)
      CALL LICHT(ICHLIM,MPLIM,TYPE,MELLIM)
      CALL LICHT(ICOEV,MPCOEV,TYPE,IGEOMC)
C
C     SEGACT MPORES*MOD
C     SEGACT MPRN*MOD
C     SEGACT MPGN*MOD
C     SEGACT MPRETN*MOD
C     SEGACT MPGAMN*MOD
C     SEGACT MPOCO
C     SEGACT MPLIM*MOD
C     SEGACT MPCOEV*MOD
C
      SEGINI, MPOCO1=MPOCO
C
C**** MPOCO1 is the "corrected cut-off" speed.
C     It is bigger than the given cut-off and the local
C     speed.
C
      IF(IGEOMF .NE. MELEMF)THEN
         WRITE(IOIMP,*) 'Anomalie dedans kfm11.eso'
         WRITE(IOIMP,*) 'Probleme de SPG'
C        21 2
C        Données incompatibles
         CALL ERREUR(21)
         GOTO 9999
      ENDIF
      IF(IGEOMC .NE. MELEMC)THEN
         WRITE(IOIMP,*) 'Anomalie dedans kfm11.eso'
         WRITE(IOIMP,*) 'Probleme de SPG'
C        21 2
C        Données incompatibles
         CALL ERREUR(21)
         GOTO 9999
      ENDIF
C
C**** On initialise GAMMA=I+(1-phi2)/phi2*Q1.(Q2)^t
C
      I1=4
      I2=NCEN
      SEGINI Q1
      SEGINI Q2
C
C**** On initialise le segment de travail D1
C
      I1=4
      SEGINI D1
C
C**** On initialise AN0, DN0, PHI2V
C
      I1=NCEN
      SEGINI AN0
      SEGINI DN0
      SEGINI PHI2V
C
C**** On initialise BN0
C
      I1=4
      I2=NCEN
      SEGINI BN0
C
C**** On initialise BN1
C
      I1=4
      I2=NCEN
      SEGINI BN1
C
C**** On initialise CN0
C
      I1=4
      I2=NCEN
      SEGINI CN0
C
C**** On initialise CN1
C
      I1=4
      I2=NCEN
      SEGINI CN1
C
C**** On initialise SIGMA0, SIGMAV
C
      I1=NFAC
      SEGINI SIGMA0
      SEGINI SIGMAV
C
C**** RES
C
      I1=4
      I2=NCEN
      SEGINI RES
      DO ICEN=1,NCEN,1
         DO ICOMP=1,4,1
            RES.VAL(ICOMP,ICEN)=MPORES.VPOCHA(ICEN,ICOMP)
         ENDDO
      ENDDO
      SEGDES MPORES
C
C***** DU (increment des variables conservatives)
C
      I1=4
      I2=NCEN
      SEGINI DU
C
C**** Les variables conservatives
C
C     CONS1.VAL contient les variables conservatives + gamma
C     i.e. rho(i) = CONS.VAL(1,i)
C          rhoux(i) = CONS.VAL(2,i)
C          rhouy(i) = CONS.VAL(3,i)
C          rhoet(i) = CONS.VAL(4,i)
C          gam(i)   = CONS.VAL(5,i)
C
      I1=5
      I2=NCEN
      SEGINI CONS0
      DO ICEN=1,NCEN,1
         CONS0.VAL(1,ICEN)=MPRN.VPOCHA(ICEN,1)
         CONS0.VAL(2,ICEN)=MPGN.VPOCHA(ICEN,1)
         CONS0.VAL(3,ICEN)=MPGN.VPOCHA(ICEN,2)
         CONS0.VAL(4,ICEN)=MPRETN.VPOCHA(ICEN,1)
         CONS0.VAL(5,ICEN)=MPGAMN.VPOCHA(ICEN,1)
      ENDDO
      SEGDES MPRN
      SEGDES MPRETN
      SEGDES MPGN
      SEGDES MPGAMN
C
      SEGINI, CONS1=CONS0
C
C**** Storage of l-1 values of F at the l Jacobi iteration
C
      IPT1=MELTFA
      SEGACT IPT1
      NSOUPO= IPT1.LISOUS(/1)
      NSOUPO=MAX(NSOUPO,1)
      JG=NSOUPO
      SEGINI MLESUP
      IFACT=0
      IF(NSOUPO .EQ. 1)THEN
         MLESUP.LECT(1)=IPT1
         NFACE=IPT1.NUM(/1)
         NELT=IPT1.NUM(/2)
         IFACT=IFACT+(NFACE*NELT)
      ELSE
         DO ISOUP = 1, NSOUPO, 1
            IPT2=IPT1.LISOUS(ISOUP)
            MLESUP.LECT(ISOUP)=IPT2
            SEGACT IPT2
            NFACE=IPT2.NUM(/1)
            NELT=IPT2.NUM(/2)
            IFACT=IFACT+(NFACE*NELT)
         ENDDO
      ENDIF
C
C***********************************************************************
C**** FIRST JACOBI ITERATION *******************************************
C***********************************************************************
C
C
C**** LOOP on ELTFA to compute:
C
C     AN0_i = 1/(2 \Omega_i) \sum_j \sigma_i,j S_j
C     SIGMA0_j
C     The cut-off in each cell
C     The reciprocal of the dual time step USDTI
C
      ICEN=0
      IFACG=0
      DO ISOUP=1,NSOUPO,1
         IPT2=MLESUP.LECT(ISOUP)
         NFACE=IPT2.NUM(/1)
         NELT=IPT2.NUM(/2)
         DO IELT=1,NELT,1
            ICEN=ICEN+1
            DO IFAC=1,NFACE,1
C              IFACG increase with IFAC,IELT,ISOUP
               IFACG = IFACG+1
               NGCF  = IPT2.NUM(IFAC,IELT)
               NLCF  = MLENT2.LECT(NGCF)
C
C************* NLCF  = numero local du centre de facel
C              NGCF  = numero global 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 = MLENT1.LECT(NGCEG)
               NLCED = MLENT1.LECT(NGCED)
C
               CNX = MPNORM.VPOCHA(NLCF,1)
               CNY = MPNORM.VPOCHA(NLCF,2)
               CTX = -1.0D0 * CNY
               CTY = CNX
C
C************* Left or right?
C
               IF(NLCEG .EQ. ICEN)THEN
                  LOGINV=.FALSE.
               ELSEIF(NLCED .EQ. ICEN)THEN
                  LOGINV=.TRUE.
               ELSE
                  WRITE(IOIMP,*) 'Anomalie dans kfm13.eso'
                  WRITE(IOIMP,*) 'Probleme de SPG'
C                 21 2
C                 Données incompatibles
                  CALL ERREUR(21)
                  GOTO 9999
               ENDIF
C************* If right, inversion
               IF(LOGINV) THEN
                  NLCED=NLCEG
                  NLCEG=ICEN
                  CNX = -1.0D0*CNX
                  CNY = -1.0D0*CNY
                  CTX = -1.0D0*CTX
                  CTY = -1.0D0*CTY
               ENDIF
C
               SURF = MPOVSU.VPOCHA(NLCF,1)
C
C************* Récupération des Etats "gauche" et "droite"
C
C
               ROG  = CONS0.VAL(1,NLCEG)
               RUXG = CONS0.VAL(2,NLCEG)
               RUYG = CONS0.VAL(3,NLCEG)
               UXG  = RUXG / ROG
               UYG  = RUYG / ROG
               UNG  = (UXG * CNX) + (UYG * CNY)
               UTG  = (UXG * CTX) + (UYG * CTY)
               RETG = CONS0.VAL(4,NLCEG)
               GAMG = CONS0.VAL(5,NLCEG)
               RECG = 0.5D0 * ((RUXG * RUXG) + (RUYG*RUYG))
               RECG = RECG / ROG
               PG   = (GAMG - 1.0D0) * (RETG - RECG)
               VOLG = MPOVOL.VPOCHA(NLCEG,1)
               DIAMG = MPOVDI.VPOCHA(NLCEG,1)
C
               ROD  = CONS0.VAL(1,NLCED)
               RUXD = CONS0.VAL(2,NLCED)
               RUYD = CONS0.VAL(3,NLCED)
               UXD  = RUXD / ROD
               UYD  = RUYD / ROD
               RETD = CONS0.VAL(4,NLCED)
               GAMD = CONS0.VAL(5,NLCED)
               RECD = 0.5D0 * ((RUXD * RUXD) + (RUYD*RUYD))
               RECD = RECD / ROD
               PD   = (GAMD - 1.0D0) * (RETD - RECD)
               UND  = (UXD * CNX) + (UYD * CNY)
               UTD  = (UXD * CTX) + (UYD * CTY)
C
               IF(NLCEG .EQ. NLCED)THEN
C           Murs ou condition limite
                  NLLIM=MLELIM.LECT(NGCF)
                  IF(NLLIM .EQ.0)THEN
C              Mur
                     UND  = -1.0D0 * UNG
                     UTD  = UTG
                     UXD  = UND * CNX + UTD * CTX
                     UYD  = UND * CNY + UTD * CTY
                     RUXD = ROD*UXD
                     RUYD = ROD*UYD
                  ELSE
                     ROD  = MPLIM.VPOCHA(NLLIM,1)
                     UXD  = MPLIM.VPOCHA(NLLIM,2)
                     UYD  = MPLIM.VPOCHA(NLLIM,3)
                     RUXD = ROD*UXD
                     RUYD = ROD*UYD
                     PD   = MPLIM.VPOCHA(NLLIM,4)
                     GAMD = GAMG
                     UND  = (UXD * CNX) + (UYD * CNY)
                     UTD  = (UXD * CTX) + (UYD * CTY)
                     RETD = ((1.0D0/(GAMD - 1.0D0))*PD)
     &                    +(0.5D0*ROD*((UXD*UXD)+(UYD*UYD)))
                  ENDIF
               ENDIF
C
C************* We check that density and pressure are positive
C
               IF((ROG .LE. 0.0D0) .OR. (ROD .LE. 0.0D0) .OR.
     &              (PG .LE. 0.0D0) .OR. (PD .LE. 0.0D0))THEN
                  WRITE(IOIMP,*) 'FMM'
                  WRITE(IOIMP,*) 'TEST'
C                 22 0
C**************** Opération malvenue. Résultat douteux
C
                  CALL ERREUR(22)
                  IPROBL = 1
                  GOTO 9998
               ENDIF
C
C************* We compute the interfacial state
C
               ROF  = 0.5D0*(ROG+ROD)
               PF   = 0.5D0*(PG+PD)
               UNF  = 0.5D0*(UNG+UND)
               GAMF = 0.5D0*(GAMG+GAMD)
C
C************* We compute SIGMA and the corrected cut-off speed
C
               UCOF=MAX(MPOCO.VPOCHA(NLCEG,1),MPOCO.VPOCHA(NLCED,1))
               UCOF=MAX(UCOF,(((UNG*UNG)+(UTG*UTG))**0.5D0))
               UCOF=MAX(UCOF,(((UND*UND)+(UTD*UTD))**0.5D0))
               IF(UCOF .GT. MPOCO1.VPOCHA(NLCEG,1))
     &            MPOCO1.VPOCHA(NLCEG,1)= UCOF
C
               QN=ABS(UNF)
               CF=(GAMF*PF/ROF)**0.5D0
               IF(UCOF .GE. CF)THEN
                  PHI2 = 1.0D0
               ELSEIF(QN .GT. CF)THEN
                  PHI2 = 1.0D0
               ELSEIF(QN .LT. UCOF)THEN
                  PHI2 = UCOF / CF
               ELSE
                  PHI2 = QN / CF
               ENDIF
               PHI2=PHI2*PHI2
C
               SIGMAF=(1.0D0-PHI2)*QN
               SIGMAF=SIGMAF*SIGMAF
               SIGMAF=SIGMAF+(4.0D0*PHI2*CF*CF)
               SIGMAF=SIGMAF**0.5D0
               SIGMAF=SIGMAF+((1.0D0+PHI2)*QN)
               SIGMAF=0.5D0*SIGMAF
C
               SIGMA0.VALVEC(NLCF)=SIGMAF
C
C************* We compute AN0
C
               AN0.VALVEC(NLCEG)=AN0.VALVEC(NLCEG)
     &                          +(0.5D0*SIGMAF*SURF/VOLG)
C
C************* We compute the dual time step
C
               UNSDT=USDTI.VALVEC(NLCEG)
               UNSDTF = SIGMAF / DIAMG
               IF(UNSDT .LT. UNSDTF) USDTI.VALVEC(NLCEG)=UNSDTF
C
C************* We compute SIGMAV
C
C              COEF=FG
               COEF=(MCOORD.XCOOR(((NGCEG-1)*3)+1)-
     &                 MCOORD.XCOOR(((NGCF-1)*3)+1))*CNX
               COEF=COEF+
     &              ((MCOORD.XCOOR(((NGCEG-1)*3)+2)-
     &              MCOORD.XCOOR(((NGCF-1)*3)+2))*CNY)
C              COEF1=FD
               COEF1=(MCOORD.XCOOR(((NGCED-1)*3)+1)-
     &              MCOORD.XCOOR(((NGCF-1)*3)+1))*CNX
               COEF1=COEF1+
     &              ((MCOORD.XCOOR(((NGCED-1)*3)+2)-
     &              MCOORD.XCOOR(((NGCF-1)*3)+2))*CNY)
C              COEF=|FG|+|FD|
               COEF=ABS(COEF)+ABS(COEF1)
               SIGMAF=0.5D0*(MPCOEV.VPOCHA(NLCEG,1)+
     &              MPCOEV.VPOCHA(NLCED,1))/COEF
               SIGMAV.VALVEC(NLCF)=SIGMAF
C
C************* We compute the contribution of SIGMAV to
C              the diagonal part of the matrix to invert
C
C              In the given report, DNO is a_i^0, i.e.
C              DN0 = AN0 (previously defined) +
C                    + (\frac{1}{{\Omega_i}}
C                       \sum_j \sigma_{V,i,j}^{0} S_j)
C                    + 1 / \Delta t^n
C
               COEF=SURF/VOLG
               DN0.VALVEC(NLCEG)=DN0.VALVEC(NLCEG)+COEF*SIGMAF
C
C
C************* We compute the flux
C              FLU.VALVEC(1:4,IFAC) flux in IFAC face
C              1 -> mass
C              2 -> momentum (x)
C              3 -> momentum (y)
C              4 -> energy
C
               RHTD = RETD+PD
               FR   = (ROD*UND)
               FRUX = (ROD*UND*UXD)+(PD*CNX)
               FRUY = (ROD*UND*UYD)+(PD*CNY)
               FRET = (UND*RHTD)
C
               COEF = SURF/(2.0D0*VOLG)
C
               BN0.VAL(1,NLCEG) = BN0.VAL(1,NLCEG)
     &                          + FR   * COEF
               BN0.VAL(2,NLCEG) = BN0.VAL(2,NLCEG)
     &                          + FRUX * COEF
               BN0.VAL(3,NLCEG) = BN0.VAL(3,NLCEG)
     &                          + FRUY * COEF
               BN0.VAL(4,NLCEG) = BN0.VAL(4,NLCEG)
     &                          + FRET * COEF
C
C************* Viscous contribution
C
               COEF=SURF*SIGMAV.VALVEC(NLCF)/VOLG
C
               BN0.VAL(1,NLCEG) = BN0.VAL(1,NLCEG)
     &                          - (COEF * ROD)
               BN0.VAL(2,NLCEG) = BN0.VAL(2,NLCEG)
     &                          - (COEF * RUXD)
               BN0.VAL(3,NLCEG) = BN0.VAL(3,NLCEG)
     &                          - (COEF *  RUYD)
               BN0.VAL(4,NLCEG) = BN0.VAL(4,NLCEG)
     &                          - (COEF *  RETD)
C
C************* On calcule CN0
C
C              CN0 in the referred report is given by
C              \sum_j
C              \left\{
C              \frac{S_j}{2 \Omega_i} ~
C              \left(
C              \sigma_{i,j}^{0}
C              \mathbf{U}_{j,R}^{n+1,m}
C              \right)
C              \right\}
C
C
               COEF=SURF*SIGMA0.VALVEC(NLCF)/(2.0D0*VOLG)
C
               CN0.VAL(1,NLCEG) = CN0.VAL(1,NLCEG)
     &                          + (COEF * ROD)
               CN0.VAL(2,NLCEG) = CN0.VAL(2,NLCEG)
     &                          + (COEF * RUXD)
               CN0.VAL(3,NLCEG) = CN0.VAL(3,NLCEG)
     &                          + (COEF *  RUYD)
               CN0.VAL(4,NLCEG) = CN0.VAL(4,NLCEG)
     &                          + (COEF *  RETD)

            ENDDO
C
C********** End of loop on the faces of the element
C
C           N.B.: ICEN = NLCEG
C
C
C           Computation of the low-Mach preconditioning Matrix
C
            CG=(GAMG*PG/ROG)**0.5D0
            UCOG=MPOCO1.VPOCHA(ICEN,1)
            QG=2.0D0*RECG/ROG
            QG=QG**0.5D0
C
            IF(UCOG .GE. CG)THEN
               PHI2 = 1.0D0
            ELSEIF(QG .GT. CG)THEN
               PHI2 = 1.0D0
            ELSEIF(QG .LT. UCOG)THEN
               PHI2 = UCOG / CG
            ELSE
               PHI2 = QG / CG
            ENDIF
            PHI2=PHI2*PHI2
            PHI2V.VALVEC(ICEN)=PHI2
C
            Q2.VAL(1,ICEN)=(0.5D0*QG*QG)
            Q2.VAL(2,ICEN)=-1.0D0*RUXG/ROG
            Q2.VAL(3,ICEN)=-1.0D0*RUYG/ROG
            Q2.VAL(4,ICEN)=1.0D0
C
            COEF1=(GAMG-1.0D0)/(CG*CG)
            COEF=(1.0D0/COEF1)+(0.5D0*QG*QG)
C
            Q1.VAL(1,ICEN)=COEF1
            Q1.VAL(2,ICEN)=COEF1*RUXG/ROG
            Q1.VAL(3,ICEN)=COEF1*RUYG/ROG
            Q1.VAL(4,ICEN)=COEF1*COEF
C
C
C           Computation of the diagonal coefficients
C
            AN0.VALVEC(ICEN)=AN0.VALVEC(ICEN)
C     &                      +(USDTI.VALVEC(ICEN)*2.0D0/DCFL)
C
            DN0.VALVEC(ICEN)=DN0.VALVEC(ICEN)+(1.0D0/DTPS)
     &                      +(USDTI.VALVEC(ICEN)*2.0D0/DCFL)
C
        ENDDO
      ENDDO
C
C***********************************************************************
C**** END FIRST JACOBI ITERATION ***************************************
C***********************************************************************
C
C***********************************************************************
C**** JACOBI LOOP ******************************************************
C***********************************************************************
C
C
      LOGJAC=.TRUE.
C
      DO IBLJAC=1,NJAC,1
C
C
C Si ICOEF=2 (i.e. si on a "LJACOF" (voir kfm1.eso))
C            alors relaxation SGS avec uniquement des balayages "forward"
C
C Si ICOEF=3 (i.e. si on a "LJACOB")
C            alors relaxation SGS avec uniquement des balayages "backward"
C
C Sinon      (i.e. si on a "LJACOFB")
C            alors relaxation SGS avec des balayages "forward" et "backward"
C
         IF(ICOEF .EQ. 2)THEN
            LOGJAC=.TRUE.
         ELSEIF(ICOEF .EQ. 3)THEN
            LOGJAC=.FALSE.
         ELSE
            LOGINV=LOGJAC
            LOGJAC=(IBLJAC / 2 * 2) .EQ. IBLJAC
            IF(LOGJAC)THEN
               LOGJAC= .NOT. LOGINV
            ELSE
               LOGJAC= LOGINV
            ENDIF
         ENDIF
C
         IF(LOGJAC)THEN
            ICEN=0
            IFACG=0
            ID=1
            ISOUP1=1
            ISOUP2=NSOUPO
         ELSE
            ICEN=NCEN+1
            IFACG=IFACT+1
            ID=-1
            ISOUP1=NSOUPO
            ISOUP2=1
         ENDIF
C
         DO ISOUP=ISOUP1,ISOUP2,ID
            IPT2=MLESUP.LECT(ISOUP)
            NFACE=IPT2.NUM(/1)
            NELT=IPT2.NUM(/2)
C
            IF(LOGJAC)THEN
               IELT1=1
               IELT2=NELT
               IFAC1=1
               IFAC2=NFACE
            ELSE
               IELT2=1
               IELT1=NELT
               IFAC2=1
               IFAC1=NFACE
            ENDIF
C
            DO IELT=IELT1,IELT2,ID
               ICEN=ICEN+ID
               DO IFAC=IFAC1,IFAC2,ID
C                 IFACG increase with IFAC,IELT,ISOUP
                  IFACG = IFACG+ID
                  NGCF  = IPT2.NUM(IFAC,IELT)
                  NLCF  = MLENT2.LECT(NGCF)
                  NGCEG = MELEFL.NUM(1,NLCF)
                  NGCED = MELEFL.NUM(3,NLCF)
                  NLCEG = MLENT1.LECT(NGCEG)
                  NLCED = MLENT1.LECT(NGCED)
C
                  CNX = MPNORM.VPOCHA(NLCF,1)
                  CNY = MPNORM.VPOCHA(NLCF,2)
                  CTX = -1.0D0 * CNY
                  CTY = CNX
C
C**************** Left or right?
C
                  IF(NLCEG .EQ. ICEN)THEN
                     LOGINV=.FALSE.
                  ELSEIF(NLCED .EQ. ICEN)THEN
                     LOGINV=.TRUE.
                  ELSE
                     WRITE(IOIMP,*) 'Anomalie dedans kfm11.eso'
                     WRITE(IOIMP,*) 'Probleme de SPG'
C                    21 2
C                    Données incompatibles
                     CALL ERREUR(21)
                     GOTO 9999
                  ENDIF
C**************** If right, inversion
                  IF(LOGINV) THEN
                     NLCED=NLCEG
                     NLCEG=ICEN
                     CNX = -1.0D0*CNX
                     CNY = -1.0D0*CNY
                     CTX = -1.0D0*CTX
                     CTY = -1.0D0*CTY
                  ENDIF
C
                  SURF   = MPOVSU.VPOCHA(NLCF,1)
                  SIGMAF = SIGMA0.VALVEC(NLCF)
C
C**************** Recuperation des Etats "gauche" et "droite"
C
                  ROG  = CONS1.VAL(1,NLCEG)
                  RUXG = CONS1.VAL(2,NLCEG)
                  RUYG = CONS1.VAL(3,NLCEG)
                  ROG  = CONS1.VAL(1,NLCEG)
                  RUXG = CONS1.VAL(2,NLCEG)
                  RUYG = CONS1.VAL(3,NLCEG)
                  UXG  = RUXG / ROG
                  UYG  = RUYG / ROG
                  UNG  = (UXG * CNX) + (UYG * CNY)
                  UTG  = (UXG * CTX) + (UYG * CTY)
                  RETG = CONS1.VAL(4,NLCEG)
                  GAMG = CONS1.VAL(5,NLCEG)
                  RECG = 0.5D0 * ((RUXG * RUXG) + (RUYG*RUYG))
                  RECG = RECG / ROG
                  PG   = (GAMG - 1.0D0) * (RETG - RECG)
                  VOLG = MPOVOL.VPOCHA(NLCEG,1)
C
                  ROD  = CONS1.VAL(1,NLCED)
                  RUXD = CONS1.VAL(2,NLCED)
                  RUYD = CONS1.VAL(3,NLCED)
                  UXD  = RUXD / ROD
                  UYD  = RUYD / ROD
                  RETD = CONS1.VAL(4,NLCED)
                  GAMD = CONS1.VAL(5,NLCED)
                  RECD = 0.5D0 * ((RUXD * RUXD) + (RUYD*RUYD))
                  RECD = RECD / ROD
                  PD   = (GAMD - 1.0D0) * (RETD - RECD)
                  UND  = (UXD * CNX) + (UYD * CNY)
                  UTD  = (UXD * CTX) + (UYD * CTY)
C
                  IF(NLCEG .EQ. NLCED)THEN
C                    Murs ou condition limite
                     NLLIM=MLELIM.LECT(NGCF)
                     IF(NLLIM .EQ.0)THEN
C                       Mur
                        UND  = -1.0D0 * UNG
                        UTD  = UTG
                        UXD  = UND * CNX + UTD * CTX
                        UYD  = UND * CNY + UTD * CTY
                        RUXD = ROD*UXD
                        RUYD = ROD*UYD
                     ELSE
                        ROD  = MPLIM.VPOCHA(NLLIM,1)
                        UXD  = MPLIM.VPOCHA(NLLIM,2)
                        UYD  = MPLIM.VPOCHA(NLLIM,3)
                        RUXD = ROD*UXD
                        RUYD = ROD*UYD
                        PD   = MPLIM.VPOCHA(NLLIM,4)
                        GAMD = GAMG
                        UND  = (UXD * CNX) + (UYD * CNY)
                        UTD  = (UXD * CTX) + (UYD * CTY)
                        RETD = ((1.0D0/(GAMD - 1.0D0))*PD)
     &                       +(0.5D0*ROD*((UXD*UXD)+(UYD*UYD)))
                     ENDIF
                  ENDIF
C
C**************** We check that density and pressure are positive
C
                  IF((ROG .LE. 0.0D0) .OR. (ROD .LE. 0.0D0) .OR.
     &                 (PG .LE. 0.0D0) .OR. (PD .LE. 0.0D0))THEN
                     WRITE(IOIMP,*) 'FMM'
                     WRITE(IOIMP,*) 'ITERATION JAC',IBLJAC
                     WRITE(IOIMP,*) 'ITERATION ICEN',ICEN
                     WRITE(IOIMP,*) ROG,PG
                     WRITE(IOIMP,*) ROD,PD
C
C                    22 0
C******************* Opération malvenue. Résultat douteux
C
                     CALL ERREUR(22)
                     IPROBL = 1
                     GOTO 9998
                  ENDIF
C
C
C**************** We compute the flux
C                 FLU.VALVEC(1:4,IFAC) flux in IFAC face
C                 1 -> mass
C                 2 -> momentum (x)
C                 3 -> momentum (y)
C                 4 -> energy
C
                  RHTD = RETD+PD
                  FR   = (ROD*UND)
                  FRUX = (ROD*UND*UXD)+(PD*CNX)
                  FRUY = (ROD*UND*UYD)+(PD*CNY)
                  FRET = (UND*RHTD)
C
                  COEF = SURF/(2.0D0*VOLG)
C
                  BN1.VAL(1,NLCEG) = BN1.VAL(1,NLCEG)
     &                             + FR   * COEF
                  BN1.VAL(2,NLCEG) = BN1.VAL(2,NLCEG)
     &                             + FRUX * COEF
                  BN1.VAL(3,NLCEG) = BN1.VAL(3,NLCEG)
     &                             + FRUY * COEF
                  BN1.VAL(4,NLCEG) = BN1.VAL(4,NLCEG)
     &                             + FRET * COEF
C
C**************** Viscous contribution
C
                  COEF=SURF*SIGMAV.VALVEC(NLCF)/VOLG
C
                  BN1.VAL(1,NLCEG) = BN1.VAL(1,NLCEG)
     &                             - (COEF * ROD)
                  BN1.VAL(2,NLCEG) = BN1.VAL(2,NLCEG)
     &                             - (COEF * RUXD)
                  BN1.VAL(3,NLCEG) = BN1.VAL(3,NLCEG)
     &                             - (COEF *  RUYD)
                  BN1.VAL(4,NLCEG) = BN1.VAL(4,NLCEG)
     &                             - (COEF *  RETD)
C
C**************** On calcule DIS1
C
C                 DIS1 in the referred report is given by
C                 \sum_j
C                 \left\{
C                 \frac{S_j}{2 \Omega_i} ~
C                 \left(
C                 \sigma_{i,j}^{0}
C                 \mathbf{U}_{j,R}^{n+1,m,l}
C                 \right)
C                 \right\}
C
C
                  COEF=SURF*SIGMA0.VALVEC(NLCF)/(2.0D0*VOLG)
C
                  CN1.VAL(1,NLCEG) = CN1.VAL(1,NLCEG)
     &                             + (COEF * ROD)
                  CN1.VAL(2,NLCEG) = CN1.VAL(2,NLCEG)
     &                             + (COEF * RUXD)
                  CN1.VAL(3,NLCEG) = CN1.VAL(3,NLCEG)
     &                             + (COEF *  RUYD)
                  CN1.VAL(4,NLCEG) = CN1.VAL(4,NLCEG)
     &                             + (COEF *  RETD)

               ENDDO
C
C************* End of the loop on the faces of the element
C
               DO ICOMP=1,4,1
                  D1.VALVEC(ICOMP) = 0.0D0
               ENDDO
C
               COEF=0.0D0
               DO I1=1,4,1
                  COEF=COEF+(Q2.VAL(I1,ICEN)
     &                *(CN1.VAL(I1,ICEN) - CN0.VAL(I1,ICEN)))
               ENDDO
               COEF=COEF*(1.0D0-PHI2V.VALVEC(ICEN))
     &                     /PHI2V.VALVEC(ICEN)
C
               DO I1=1,4,1
                  D1.VALVEC(I1)=(CN1.VAL(I1,ICEN)
     &                         - CN0.VAL(I1,ICEN))
     &                         +(COEF*Q1.VAL(I1,ICEN))
c                  write(6,*)D1.VALVEC(I1)
               ENDDO
C
C************* At this moment D1 = \Gamma \cdot (\Delta CN1)
C
               DO ICOMP=1,4,1
C
C************* First increment
C
                  D1.VALVEC(ICOMP) = (D1.VALVEC(ICOMP)
     &                 +RES.VAL(ICOMP,ICEN)
     &                 +(BN0.VAL(ICOMP,ICEN)-BN1.VAL(ICOMP,ICEN)))
     &                 /(AN0.VALVEC(ICEN)+DN0.VALVEC(ICEN))
               ENDDO
C
C************* At this moment, in the referred report,
C              D1 is (\delta \mathbf{U}')^{l+1}_i in section A3
C
C************* Second increment
C
               COEF=0.0D0
               DO I1=1,4,1
                  COEF=COEF+(Q2.VAL(I1,ICEN)*D1.VALVEC(I1))
               ENDDO
C
C************* COEF is the scalar product between Q2 and the first
C              increment.
C              Let us multiply it by
C              b_i^0 / (a_i^0 + b_i^0)
C
               COEF=COEF*AN0.VALVEC(ICEN)
     &                  *(PHI2V.VALVEC(ICEN)-1.0D0)
               COEF=COEF/(AN0.VALVEC(ICEN)+(DN0.VALVEC(ICEN)
     &                  *PHI2V.VALVEC(ICEN)))
C
               DO I1=1,4,1
                  D1.VALVEC(I1)=D1.VALVEC(I1)+COEF*Q1.VAL(I1,ICEN)
               ENDDO
C
C************* D1 = \delta \mathbf{U}
C
               DO ICOMP=1,4,1
                  DU.VAL(ICOMP,ICEN)    = D1.VALVEC(ICOMP)
C
                  CONS1.VAL(ICOMP,ICEN) = CONS0.VAL(ICOMP,ICEN)
     &                                   + D1.VALVEC(ICOMP)
C
                  BN1.VAL(ICOMP,ICEN)   = 0.0D0
                  CN1.VAL(ICOMP,ICEN)   = 0.0D0
               ENDDO
C
C********** End of the loop on the elements
            ENDDO
C
C******* End of the loop on the Domains
         ENDDO
C
C**** End of the Jacobi loop
      ENDDO
C
 9998 CONTINUE
      DO ICEN=1,NCEN,1
         DO ICOMP=1,4,1
            MPODU.VPOCHA(ICEN,ICOMP) = DU.VAL(ICOMP,ICEN)
         ENDDO
      ENDDO
C
      SEGSUP DU
      SEGSUP RES
      SEGSUP CONS0
      SEGSUP CONS1
      SEGSUP Q1
      SEGSUP Q2
      SEGSUP SIGMA0
      SEGSUP AN0
      SEGSUP AN0
      SEGSUP BN0
      SEGSUP BN1
      SEGSUP CN0
      SEGSUP CN1
      SEGSUP DN0
      SEGSUP PHI2V
      SEGSUP USDTI
C
      SEGDES MPOVSU
      SEGDES MPOVDI
      SEGDES MPOVOL
      SEGDES MPNORM
      SEGDES MPOCO
      SEGDES MPCOEV
C
      SEGSUP MPOCO1
C
      SEGDES MPODU
C
      SEGDES MELEFL
      SEGSUP MLENT1
      SEGSUP MLENT2
C
      SEGSUP MLELIM
      IF(MPLIM .GT. 0) SEGDES MPLIM
C
      DO ISOUP=1,NSOUPO,1
         IPT2=MLESUP.LECT(ISOUP)
         SEGDES IPT2
      ENDDO
      IPT2 = MELTFA
      IF(ISOUP .GT. 1) SEGDES IPT2
      SEGSUP MLESUP
C
 9999 CONTINUE
      RETURN
      END
C