C KFM11     SOURCE    OF166741  24/12/13    21:16:04     12097          
      SUBROUTINE KFM11(IRES,IRN,IGN,IRETN,IGAMN,IVCO,
     &     ICHPSU,ICHPDI,ICHPVO,INORM,
     &     MELEMC,MELEMF,MELEFE,DTPS,DCFL,
     &     MELLIM,ICHLIM,NJAC,ICOEV,
     &     IDU,IPROBL)
C************************************************************************
C
C PROJET            :  CASTEM 2000
C
C NOM               :  KFM11
C
C DESCRIPTION       :  Voir aussi KFM1
C
C                      Cas deux dimensions, gaz "calorically perfect"
C                      Méthode de ralaxation de type Jacobi par point
C
C LANGAGE           :  FORTRAN 77 + ESOPE 2000 (avec estensions CISI)
C
C AUTEUR            :  A. BECCANTINI, DRN/DMT/SEMT/TTMF
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     MELEFE  : 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 : Avril 2002  : creation
C              Janvier 2003: implementation de condition aux limites
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
     &     ,MELEMC,MELEMF,MELEFE,IDU,NFAC
     &     ,IGEOMF,IGEOMC, NGCEG, NGCED, NGCF, NLCF, NLCF1, NLCEG, NLCED
     &     ,ICEN,NCEN,ICHLIM,MELLIM,NLLIM,IRES,NJAC,ICOMP
     &     ,I1,I2, IBLJAC, IPROBL, IVCO, ICOEV
      REAL*8 ROG, RUXG, RUYG, UNG, RETG, GAMG, REC, PG, VOLG, DIAMG
     &     , ROD, RUXD, RUYD, UND, RETD, GAMD, PD, VOLD, DIAMD
     &     , CNX, CNY, DCFL, DTPS
     &     , SURF, SIGMAF
     &     , UNSDT, UNSDTF, UXD, UYD
     &     , FR, FRUX, FRUY, FRET, COEF
     &     , CTX, CTY, RHTD, RHTG, UTD, UTG, UXG, UYG, QG, CG, UCOG
     &     , ROF,PF,UNF,  GAMF, UCOF, CF, PHI2, COEF1
     &     , QN
      CHARACTER*8 TYPE
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, MPCOEV.MPOVAL
-INC SMCOORD
-INC SMELEME
-INC SMLMOTS
-INC SMLENTI
      POINTEUR MLELIM.MLENTI
C
C**** POINTEUR vecteurs
C
      SEGMENT VEC
      REAL*8 VALVEC(I1)
      ENDSEGMENT
      POINTEUR D1.VEC,SIGMA0.VEC, AN0.VEC, DN0.VEC
     &     ,USDTI.VEC, PHI2V.VEC, SIGMAV.VEC
C     &     ,D2.VEC
C
C**** POINTEUR matriciels
C
      SEGMENT MAT
      REAL*8 VAL(I1,I2)
      ENDSEGMENT
      POINTEUR CONS1.MAT, CONS2.MAT, BN0.MAT, BN1.MAT
     &     ,CN0.MAT, CN1.MAT, DU.MAT
     &     , RES.MAT
     &     , Q1.MAT, Q2.MAT
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
      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 nbpts 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+Q1.(Q2)^t
C
      I1=4
      I2=NCEN
      SEGINI Q1
      SEGINI Q2
C
C**** On initialise le segment de travail D1 et D2
C
      I1=4
      SEGINI D1
C      SEGINI D2
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 CONS1
      DO ICEN=1,NCEN,1
         CONS1.VAL(1,ICEN)=MPRN.VPOCHA(ICEN,1)
         CONS1.VAL(2,ICEN)=MPGN.VPOCHA(ICEN,1)
         CONS1.VAL(3,ICEN)=MPGN.VPOCHA(ICEN,2)
         CONS1.VAL(4,ICEN)=MPRETN.VPOCHA(ICEN,1)
         CONS1.VAL(5,ICEN)=MPGAMN.VPOCHA(ICEN,1)
      ENDDO
      SEGDES MPRN
      SEGDES MPRETN
      SEGDES MPGN
      SEGDES MPGAMN
C
      SEGINI, CONS2=CONS1
C
C**** BOUCLE JACOBI
C
      DO IBLJAC=1,NJAC,1
C
C**** BOUCLE SUR FACEL pour le calcul de AN0, BN0
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)
C
C******* NLCF != NLCF1 -> l'auteur (MOI) n'a rien compris.
C
            IF(NLCF .NE. NLCF1)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
            CNX = MPNORM.VPOCHA(NLCF,1)
            CNY = MPNORM.VPOCHA(NLCF,2)
            CTX = -1.0D0 * CNY
            CTY = CNX
            SURF= MPOVSU.VPOCHA(NLCF,1)
C
C******* Recuperation des Etats "gauche" et "droite"
C
C
            ROG = CONS2.VAL(1,NLCEG)
            RUXG = CONS2.VAL(2,NLCEG)
            RUYG = CONS2.VAL(3,NLCEG)
            UXG = RUXG / ROG
            UYG = RUYG / ROG
            UNG = (UXG * CNX) + (UYG * CNY)
            UTG = (UXG * CTX) + (UYG * CTY)
            RETG = CONS2.VAL(4,NLCEG)
            GAMG = CONS2.VAL(5,NLCEG)
            REC = 0.5D0 * ((RUXG * RUXG) + (RUYG*RUYG))
            REC = REC / ROG
            PG = (GAMG - 1.0D0) * (RETG - REC)
            VOLG = MPOVOL.VPOCHA(NLCEG,1)
            DIAMG = MPOVDI.VPOCHA(NLCEG,1)
C
            ROD = CONS2.VAL(1,NLCED)
            RUXD = CONS2.VAL(2,NLCED)
            RUYD = CONS2.VAL(3,NLCED)
            UXD = RUXD / ROD
            UYD = RUYD / ROD
            RETD = CONS2.VAL(4,NLCED)
            GAMD = CONS2.VAL(5,NLCED)
            REC = 0.5D0 * ((RUXD * RUXD) + (RUYD*RUYD))
            REC = REC / ROD
            PD = (GAMD - 1.0D0) * (RETD - REC)
            VOLD = MPOVOL.VPOCHA(NLCED,1)
            DIAMD = MPOVDI.VPOCHA(NLCED,1)
            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)))
                  VOLD = VOLG
               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'
C            22 0
C*********** Opération malvenue. Résultat douteux
C
               CALL ERREUR(22)
               IPROBL = 1
               GOTO 9998
            ENDIF
C
            RHTG=RETG+PG
            RHTD=RETD+PD
C
C********** FIRST JACOBI ITERATION
C
            IF(IBLJAC .EQ. 1)THEN
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              Here SIGMA is the spectral of the preconditioned Roe
C              matrix (|\Gamma^{-1} A| in the referred report)
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
               IF(UCOF .GT. MPOCO1.VPOCHA(NLCED,1)) MPOCO1.VPOCHA(NLCED
     $              ,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
c              write(*,*) qn,sigmaf,phi2
c
c              write(*,*) 'UCOF',UCOF
c              write(*,*) 'PHI2',PHI2
c              write(*,*) 'JACOBI,SIGMAF',IBLJAC,SIGMAF
c
C
C************* We store SIGMA into SIGMA0
C
               SIGMA0.VALVEC(NLCF)=SIGMAF
C
C************* We compute the dual time step
C
               UNSDT=USDTI.VALVEC(NLCEG)
               UNSDTF = SIGMAF / DIAMG
               IF(UNSDT .LT. UNSDTF) USDTI.VALVEC(NLCEG)=UNSDTF
               UNSDT=USDTI.VALVEC(NLCED)
               UNSDTF = SIGMAF / DIAMD
               IF(UNSDT .LT. UNSDTF) USDTI.VALVEC(NLCED)=UNSDTF
C
C************* We compute AN0
C              ANO in the referred report is equal to
C              \left(
C              \frac{1}{\Delta \tau_i ^0}
C              + \frac{1}{2{\Omega_i}}
C              \sum_j \sigma_{i,j}^{0} S_j
C              \right)
C
               AN0.VALVEC(NLCEG)=AN0.VALVEC(NLCEG) +
     &              (0.5D0*SIGMAF*SURF/VOLG)
               IF(NLCEG .NE. NLCED)
     &              AN0.VALVEC(NLCED)=AN0.VALVEC(NLCED) +
     &              (0.5D0*SIGMAF*SURF/VOLD)
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              write(*,*) 'sigmaf =',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
               IF(NLCEG .NE. NLCED)THEN
                  COEF=SURF/VOLD
                  DN0.VALVEC(NLCED)=DN0.VALVEC(NLCED)+COEF*SIGMAF
               ENDIF
            ENDIF
C
C********** ALL JACOBI ITERATIONS
C
C
C********** On calcule BN1
C
C           BN1 in the referred report is given by
C           \sum_j
C           \left\{
C           \frac{S_j}{2 \Omega_i} ~
C           \left(
C           \mathcal{F}(\mathbf{U}_{j,R}^{\ell})
C           -
C           2 \sigma_{V,i,j}^{0}
C           \mathbf{U}_{j,R}^{\ell}
C           _right)
C           \right\}
C
C           Central contribution
C
            FR=(ROD*UND)
            FRUX=(ROD*UND*UXD)+(PD*CNX)
            FRUY=(ROD*UND*UYD)+(PD*CNY)
            FRET=(UND*RHTD)
C
            COEF=SURF/(2.0D0*VOLG)
            BN1.VAL(1,NLCEG) = BN1.VAL(1,NLCEG) +
     &           (COEF * FR)
            BN1.VAL(2,NLCEG) = BN1.VAL(2,NLCEG) +
     &           (COEF * FRUX)
            BN1.VAL(3,NLCEG) = BN1.VAL(3,NLCEG) +
     &           (COEF *  FRUY)
            BN1.VAL(4,NLCEG) = BN1.VAL(4,NLCEG) +
     &           (COEF *  FRET)
C
C           Viscous contribution
C
            COEF=SURF*SIGMAV.VALVEC(NLCF)/VOLG
            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 CN1
C
C           CN1 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}^{\ell}
C           \right)
C           \right\}
C
C
            COEF=SURF*SIGMA0.VALVEC(NLCF)/(2.0D0*VOLG)
            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)
C
            IF(NLCEG .NE. NLCED)THEN
C
C********** Domain interieur
C
C              Central contribution
C
               FR=(ROG*UNG)
               FRUX=(ROG*UNG*UXG)+(PG*CNX)
               FRUY=(ROG*UNG*UYG)+(PG*CNY)
               FRET=(UNG*RHTG)
               COEF=SURF/(2.0D0*VOLD)
C
               BN1.VAL(1,NLCED) = BN1.VAL(1,NLCED) -
     &              (COEF * FR)
               BN1.VAL(2,NLCED) = BN1.VAL(2,NLCED) -
     &              (COEF * FRUX)
               BN1.VAL(3,NLCED) = BN1.VAL(3,NLCED) -
     &              (COEF *  FRUY)
               BN1.VAL(4,NLCED) = BN1.VAL(4,NLCED) -
     &              (COEF *  FRET)
C
C              Viscous contribution
C
               COEF=SURF*SIGMAV.VALVEC(NLCF)/VOLD
               BN1.VAL(1,NLCED) = BN1.VAL(1,NLCED) -
     &              (COEF * ROG)
               BN1.VAL(2,NLCED) = BN1.VAL(2,NLCED) -
     &              (COEF * RUXG)
               BN1.VAL(3,NLCED) = BN1.VAL(3,NLCED) -
     &              (COEF *  RUYG)
               BN1.VAL(4,NLCED) = BN1.VAL(4,NLCED) -
     &              (COEF *  RETG)
C
C              Numerical diffusity contribution
C
               COEF=SURF*SIGMA0.VALVEC(NLCF)/(2.0D0*VOLD)
               CN1.VAL(1,NLCED) = CN1.VAL(1,NLCED) +
     &              (COEF * ROG)
               CN1.VAL(2,NLCED) = CN1.VAL(2,NLCED) +
     &              (COEF * RUXG)
               CN1.VAL(3,NLCED) = CN1.VAL(3,NLCED) +
     &              (COEF *  RUYG)
               CN1.VAL(4,NLCED) = CN1.VAL(4,NLCED) +
     &              (COEF *  RETG)

            ENDIF
C
C********** Fin boucle sur les faces
C
         ENDDO

C
C********** Loop on the centre
C           Computation of GAMMA
C                          DMAT
C
         DO ICEN=1,NCEN,1
C
            IF(IBLJAC .EQ. 1)THEN
C
               ROG = CONS2.VAL(1,ICEN)
               RUXG = CONS2.VAL(2,ICEN)
               RUYG = CONS2.VAL(3,ICEN)
               UXG=RUXG/ROG
               UYG=RUYG/ROG
               RETG = CONS2.VAL(4,ICEN)
               GAMG = CONS2.VAL(5,ICEN)
               REC = 0.5D0 * ((RUXG * RUXG) + (RUYG*RUYG))
               REC = REC / ROG
               PG = (GAMG - 1.0D0) * (RETG - REC)
C
C************* Note that the Cut-off is the corrected one
C
               CG=(GAMG*PG/ROG)**0.5D0
               UCOG=MPOCO1.VPOCHA(ICEN,1)
               QG=2.0D0*REC/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
C************* We compute GAMMA
C              GAMMA=I+((1-PHI2)/PHI2*Q1.(Q2)^t)
C
C              In the referred report, Q2 is the line vector
C
C              ( \frac{q^2}{2} , -u ,  -v , 1)
C
C              and Q1 is the column vector
C
C              \frac{\gamma - 1}{c^2}
C              \left(
C              \begin{array}{c}
C                    1 \\
C                    u \\
C                    v \\
C                    H \\
C              \end{array}
C
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              COEF=HT
               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
               AN0.VALVEC(ICEN)=USDTI.VALVEC(ICEN)*2.0D0/DCFL
     &              +AN0.VALVEC(ICEN)
               DN0.VALVEC(ICEN)=DN0.VALVEC(ICEN)+(1.0D0/DTPS)
C
C
C************* We save BN1, CN1 into BN0 and CN0
C
               DO ICOMP=1,4,1
                  BN0.VAL(ICOMP,ICEN) = BN1.VAL(ICOMP,ICEN)
               ENDDO
C
               DO ICOMP=1,4,1
                  CN0.VAL(ICOMP,ICEN) = CN1.VAL(ICOMP,ICEN)
               ENDDO
C
            ENDIF
C
C********** IBLJAC => 1
C
            DO ICOMP=1,4,1
               D1.VALVEC(ICOMP) = 0.0D0
C               D2.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))
            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
C               D2.VALVEC(I1)=D1.VALVEC(I1)+COEF*Q1.VAL(I1,ICEN)
               D1.VALVEC(I1)=D1.VALVEC(I1)+COEF*Q1.VAL(I1,ICEN)
            ENDDO
C
C********** D2 = \delta \mathbf{U}
C
            DO ICOMP=1,4,1
C               DU.VAL(ICOMP,ICEN)=D2.VALVEC(ICOMP)
C               CONS2.VAL(ICOMP,ICEN) = CONS1.VAL(ICOMP,ICEN) +
C     &              D2.VALVEC(ICOMP)
               DU.VAL(ICOMP,ICEN)=D1.VALVEC(ICOMP)
               CONS2.VAL(ICOMP,ICEN) = CONS1.VAL(ICOMP,ICEN) +
     &              D1.VALVEC(ICOMP)
               BN1.VAL(ICOMP,ICEN)=0.0D0
               CN1.VAL(ICOMP,ICEN)=0.0D0
            ENDDO
         ENDDO
C
C**** Fin boucle JACOBI
C
      ENDDO
C
 9998 CONTINUE
c      IF(IPROBL .EQ. 0)THEN
         DO ICEN=1,NCEN,1
            DO ICOMP=1,4,1
               MPODU.VPOCHA(ICEN,ICOMP) = DU.VAL(ICOMP,ICEN)
            ENDDO
         ENDDO
c      ELSE
c         DO ICEN=1,NCEN,1
c            DO ICOMP=1,4,1
c               MPODU.VPOCHA(ICEN,ICOMP) = 0.0D0
c            ENDDO
c         ENDDO
c      ENDIF
C
      SEGSUP DU
      SEGSUP RES
      SEGSUP CONS1
      SEGSUP CONS2
      SEGSUP Q1
      SEGSUP Q2
      SEGSUP D1
C      SEGSUP D2
      SEGSUP SIGMA0
      SEGSUP SIGMAV
      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 IPT2
      SEGSUP MLENT1
      SEGSUP MLENT2
C
      SEGSUP MLELIM
      IF(MPLIM .GT. 0) SEGDES MPLIM
C
 9999 CONTINUE
      RETURN
      END
C