Sources Esope | Cast3M

kfm14
C KFM14     SOURCE    OF166741  24/12/13    21:16:07     12097          
C KFM14     SOURCE    CHAT      05/01/13    00:55:31     5004
      SUBROUTINE KFM14(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               :  KFM14
C
C DESCRIPTION       :  Voir aussi KFM1 et KFM12
C
C                      Cas 3D, gaz "calorically perfect"
C                      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 création
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, RUZG, RETG, GAMG, RHTG, RECG
     &     , ROD, RUXD, RUYD, RUZD, RETD, GAMD, RHTD, RECD
     &     , UXG, UYG, UZG, UNG, UTG, UVG, PG, VOLG, DIAMG, SURF
     &     , UXD, UYD, UZD, UND, UTD, UVD, PD
     &     , CNX, CNY, CNZ, CTX, CTY, CTZ, CVX, CVY, CVZ
     &     , DCFL, DTPS, UNSDT, UNSDTF
     &     , FR, FRUX, FRUY, FRUZ, FRET
     &     , QG, CG, UCOG, ROF, PF, UNF, GAMF, UCOF, CF, SIGMAF
     &     , QN, PHI2, COEF, COEF1
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 dans kfm14.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 dans kfm14.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=5
      I2=NCEN
      SEGINI Q1
      SEGINI Q2
C
C**** On initialise le segment de travail D1
C
      I1=5
      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=5
      I2=NCEN
      SEGINI BN0
C
C**** On initialise BN1
C
      I1=5
      I2=NCEN
      SEGINI BN1
C
C**** On initialise CN0
C
      I1=5
      I2=NCEN
      SEGINI CN0
C
C**** On initialise CN1
C
      I1=5
      I2=NCEN
      SEGINI CN1
C
C**** On initialise SIGMA0, SIGMAV
C
      I1=NFAC
      SEGINI SIGMA0
      SEGINI SIGMAV
C
C**** RES
C
      I1=5
      I2=NCEN
      SEGINI RES
      DO ICEN=1,NCEN,1
         DO ICOMP=1,5,1
            RES.VAL(ICOMP,ICEN)=MPORES.VPOCHA(ICEN,ICOMP)
         ENDDO
      ENDDO
      SEGDES MPORES
C
C***** DU (increment des variables conservatives)
C
      I1=5
      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          rhouz(i) = CONS.VAL(4,i)
C          rhoet(i) = CONS.VAL(5,i)
C          gam(i)   = CONS.VAL(6,i)
C
      I1=6
      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)=MPGN.VPOCHA(ICEN,3)
         CONS0.VAL(5,ICEN)=MPRETN.VPOCHA(ICEN,1)
         CONS0.VAL(6,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,7)
               CNY = MPNORM.VPOCHA(NLCF,8)
               CNZ = MPNORM.VPOCHA(NLCF,9)
               CTX = MPNORM.VPOCHA(NLCF,1)
               CTY = MPNORM.VPOCHA(NLCF,2)
               CTZ = MPNORM.VPOCHA(NLCF,3)
               CVX = MPNORM.VPOCHA(NLCF,4)
               CVY = MPNORM.VPOCHA(NLCF,5)
               CVZ = MPNORM.VPOCHA(NLCF,6)
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
                  CNZ = -1.0D0*CNZ
                  CTX = -1.0D0*CTX
                  CTY = -1.0D0*CTY
                  CTZ = -1.0D0*CTZ
                  CVX = -1.0D0*CVX
                  CVY = -1.0D0*CVY
                  CVZ = -1.0D0*CVZ
               ENDIF
C
               SURF = MPOVSU.VPOCHA(NLCF,1)
C
C************* Recuperation des Etats "gauche" et "droite"
C
C
               ROG  = CONS0.VAL(1,NLCEG)
               RUXG = CONS0.VAL(2,NLCEG)
               RUYG = CONS0.VAL(3,NLCEG)
               RUZG = CONS0.VAL(4,NLCEG)
               RETG = CONS0.VAL(5,NLCEG)
               GAMG = CONS0.VAL(6,NLCEG)
C
               UXG  = RUXG / ROG
               UYG  = RUYG / ROG
               UZG  = RUZG / ROG
C
               UNG  = (UXG * CNX) + (UYG * CNY) + (UZG * CNZ)
               UTG  = (UXG * CTX) + (UYG * CTY) + (UZG * CTZ)
               UVG  = (UXG * CVX) + (UYG * CVY) + (UZG * CVZ)
C
               RECG = 0.5D0 * (RUXG**2 + RUYG**2 + RUZG**2)
               RECG = RECG / ROG
               PG   = (GAMG - 1.0D0) * (RETG - RECG)
C
               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)
               RUZD = CONS0.VAL(4,NLCED)
               RETD = CONS0.VAL(5,NLCED)
               GAMD = CONS0.VAL(6,NLCED)
C
               UXD  = RUXD / ROD
               UYD  = RUYD / ROD
               UZD  = RUZD / ROD
C
               UND  = (UXD * CNX) + (UYD * CNY) + (UZD * CNZ)
               UTD  = (UXD * CTX) + (UYD * CTY) + (UZD * CTZ)
               UVD  = (UXD * CVX) + (UYD * CVY) + (UZD * CVZ)
C
               RECD = 0.5D0 * (RUXD**2 + RUYD**2 + RUZD**2)
               RECD = RECD / ROD
               PD   = (GAMD - 1.0D0) * (RETD - RECD)
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
                     UVD  = UVG
                     UXD  = UND * CNX + UTD * CTX + UVD * CVX
                     UYD  = UND * CNY + UTD * CTY + UVD * CVY
                     UZD  = UND * CNZ + UTD * CTZ + UVD * CVZ
C
                     RUXD = ROD * UXD
                     RUYD = ROD * UYD
                     RUZD = ROD * UZD
C
                  ELSE
                     ROD  = MPLIM.VPOCHA(NLLIM,1)
                     UXD  = MPLIM.VPOCHA(NLLIM,2)
                     UYD  = MPLIM.VPOCHA(NLLIM,3)
                     UZD  = MPLIM.VPOCHA(NLLIM,4)
                     PD   = MPLIM.VPOCHA(NLLIM,5)
C
                     RUXD = ROD*UXD
                     RUYD = ROD*UYD
                     RUZD = ROD * UZD
C
                     GAMD = GAMG
C
                     UND  = (UXD * CNX) + (UYD * CNY) + (UZD * CNZ)
                     UTD  = (UXD * CTX) + (UYD * CTY) + (UZD * CTZ)
                     UVD  = (UXD * CVX) + (UYD * CVY) + (UZD * CVZ)
C
                     RETD = ((1.0D0/(GAMD - 1.0D0))*PD)
     &                    +(0.5D0*ROD*(UXD**2+UYD**2+UZD**2))
                  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**2+UTG**2+UVG**2)**0.5D0))
               UCOF=MAX(UCOF,((UND**2+UTD**2+UVD**2)**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)*4)+1)
     &                 -MCOORD.XCOOR(((NGCF-1)*4)+1)) * CNX
               COEF  = COEF
     &               + (MCOORD.XCOOR(((NGCEG-1)*4)+2)
     &                 -MCOORD.XCOOR(((NGCF-1)*4)+2)) * CNY
                COEF  = COEF
     &               + (MCOORD.XCOOR(((NGCEG-1)*4)+3)
     &                 -MCOORD.XCOOR(((NGCF-1)*4)+3)) * CNZ
C
C              COEF1=FD
               COEF1 = (MCOORD.XCOOR(((NGCED-1)*4)+1)
     &                 -MCOORD.XCOOR(((NGCF-1)*4)+1)) * CNX
               COEF1 = COEF1
     &               + (MCOORD.XCOOR(((NGCED-1)*4)+2)
     &                 -MCOORD.XCOOR(((NGCF-1)*4)+2)) * CNY
                 COEF1 = COEF1
     &               + (MCOORD.XCOOR(((NGCED-1)*4)+3)
     &                 -MCOORD.XCOOR(((NGCF-1)*4)+3)) * CNZ
C
C              COEF=|FG|+|FD|
               COEF   = ABS(COEF)+ABS(COEF1)
               SIGMAF = 0.5D0 * (MPCOEV.VPOCHA(NLCEG,1)
     &                          +MPCOEV.VPOCHA(NLCED,1)) / COEF
C
C
               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)
               FRUZ = (ROD * UND * UZD) + (PD * CNZ)
               FRET = (UND * RHTD)
C
               COEF = 0.5D0 * SURF / 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)
     &                          + FRUZ * COEF
               BN0.VAL(5,NLCEG) = BN0.VAL(5,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 * RUZD)
               BN0.VAL(5,NLCEG) = BN0.VAL(5,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 *  RUZD)
               CN0.VAL(5,NLCEG) = CN0.VAL(5,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 * RUZG / ROG
            Q2.VAL(5,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 * RUZG / ROG
            Q1.VAL(5,ICEN) = COEF1 * COEF
C
C
C           Computation of the diagonal coefficients
C
            AN0.VALVEC(ICEN) = AN0.VALVEC(ICEN)
     &                       + (USDTI.VALVEC(ICEN) * 2.0D0 / DCFL)
C
            DN0.VALVEC(ICEN) = DN0.VALVEC(ICEN) + (1.0D0 / DTPS)
c     &                       + (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            selon la parité de la sous-itération en cours (IBLJAC)...
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,7)
                  CNY = MPNORM.VPOCHA(NLCF,8)
                  CNZ = MPNORM.VPOCHA(NLCF,9)
                  CTX = MPNORM.VPOCHA(NLCF,1)
                  CTY = MPNORM.VPOCHA(NLCF,2)
                  CTZ = MPNORM.VPOCHA(NLCF,3)
                  CVX = MPNORM.VPOCHA(NLCF,4)
                  CVY = MPNORM.VPOCHA(NLCF,5)
                  CVZ = MPNORM.VPOCHA(NLCF,6)
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 kfm14.eso'
                     WRITE(IOIMP,*) 'Problème 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
                     CNZ = -1.0D0*CNZ
                     CTX = -1.0D0*CTX
                     CTY = -1.0D0*CTY
                     CTZ = -1.0D0*CTZ
                     CVX = -1.0D0*CVX
                     CVY = -1.0D0*CVY
                     CVZ = -1.0D0*CVZ
                  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)
                  RUZG = CONS1.VAL(4,NLCEG)
                  RETG = CONS1.VAL(5,NLCEG)
                  GAMG = CONS1.VAL(6,NLCEG)
C
                  UXG  = RUXG / ROG
                  UYG  = RUYG / ROG
                  UZG  = RUZG / ROG
C
                  UNG  = (UXG * CNX) + (UYG * CNY) + (UZG * CNZ)
                  UTG  = (UXG * CTX) + (UYG * CTY) + (UZG * CTZ)
                  UVG  = (UXG * CVX) + (UYG * CVY) + (UZG * CVZ)
C
                  RECG = 0.5D0 * (RUXG**2 + RUYG**2 + RUZG**2)
                  RECG = RECG / ROG
                  PG   = (GAMG - 1.0D0) * (RETG - RECG)
C
                  VOLG = MPOVOL.VPOCHA(NLCEG,1)
C
                  ROD  = CONS1.VAL(1,NLCED)
                  RUXD = CONS1.VAL(2,NLCED)
                  RUYD = CONS1.VAL(3,NLCED)
                  RUZD = CONS1.VAL(4,NLCED)
                  RETD = CONS1.VAL(5,NLCED)
                  GAMD = CONS1.VAL(6,NLCED)
C
                  UXD  = RUXD / ROD
                  UYD  = RUYD / ROD
                  UZD  = RUZD / ROD
C
                  UND  = (UXD * CNX) + (UYD * CNY) + (UZD * CNZ)
                  UTD  = (UXD * CTX) + (UYD * CTY) + (UZD * CTZ)
                  UVD  = (UXD * CVX) + (UYD * CVY) + (UZD * CVZ)
C
                  RECD = 0.5D0 * (RUXD**2 + RUYD**2 + RUZD**2)
                  RECD = RECD / ROD
                  PD   = (GAMD - 1.0D0) * (RETD - RECD)
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
                        UVD  = UVG
                        UXD  = UND * CNX + UTD * CTX + UVD * CVX
                        UYD  = UND * CNY + UTD * CTY + UVD * CVY
                        UZD  = UND * CNZ + UTD * CTZ + UVD * CVZ
C
                        RUXD = ROD*UXD
                        RUYD = ROD*UYD
                        RUZD = ROD * UZD
C
                     ELSE
                        ROD  = MPLIM.VPOCHA(NLLIM,1)
                        UXD  = MPLIM.VPOCHA(NLLIM,2)
                        UYD  = MPLIM.VPOCHA(NLLIM,3)
                        UZD  = MPLIM.VPOCHA(NLLIM,4)
                        PD   = MPLIM.VPOCHA(NLLIM,5)
C
                        RUXD = ROD*UXD
                        RUYD = ROD*UYD
                        RUZD = ROD * UZD
C
                        GAMD = GAMG
C
                        UND  = (UXD * CNX) + (UYD * CNY) + (UZD * CNZ)
                        UTD  = (UXD * CTX) + (UYD * CTY) + (UZD * CTZ)
                        UVD  = (UXD * CVX) + (UYD * CVY) + (UZD * CVZ)
                        RETD = ((1.0D0/(GAMD - 1.0D0))*PD)
     &                       +(0.5D0*ROD*(UXD**2+UYD**2+UZD**2))
                     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)
                    FRUZ = (ROD * UND * UZD) + (PD * CNZ)
                  FRET = (UND * RHTD)
C
                  COEF = 0.5D0 * SURF / 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)
     &                             + FRUZ * COEF
                  BN1.VAL(5,NLCEG) = BN1.VAL(5,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 *  RUZD)
                  BN1.VAL(5,NLCEG) = BN1.VAL(5,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 *  RUZD)
                  CN1.VAL(5,NLCEG) = CN1.VAL(5,NLCEG)
     &                             + (COEF *  RETD)
 
               ENDDO
C
C************* End of the loop on the faces of the element
C
               DO ICOMP=1,5,1
                  D1.VALVEC(ICOMP) = 0.0D0
               ENDDO
C
               COEF=0.0D0
               DO I1=1,5,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,5,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,5,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,5,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,5,1
                  D1.VALVEC(I1) = D1.VALVEC(I1)
     &                          + COEF * Q1.VAL(I1,ICEN)
               ENDDO
C
C************* D1 = \delta \mathbf{U}
C
               DO ICOMP=1,5,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,5,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