Sources Esope | Cast3M

kfm13
C KFM13     SOURCE    CB215821  20/11/25    13:31:15     10792          
C KFM13     SOURCE    CHAT      05/01/13    00:55:24     5004
      SUBROUTINE KFM13(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               :  KFM13
C
C DESCRIPTION       :  Voir aussi KFM1
C
C                      Cas 3D, gaz "calorically perfect"
C                      Méthode implicite sans Matrice résolue
C                      par Jacobi par point.
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     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 : 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
C
C**** Variables de COOPTIO
C
C      INTEGER IPLLB, IERPER, IERMAX, IERR, INTERR
C     &        ,IOTER,   IOLEC,  IOIMP,     IOCAR,  IOACQ
C     &        ,IOPER,   IOSGB,  IOGRA,     IOSAU,  IORES
C     &        ,IECHO,   IIMPI,  IOSPI
C     &        ,IDIM
CC     &        ,MCOORD
C     &        ,IFOMOD, NIFOUR, IFOUR, NSDPGE, IONIVE
C     &        ,NGMAXY, IZROSF, ISOTYP, IOSCR,LTEXLU
C     &        ,NORINC,NORVAL,NORIND,NORVAD
C     &        ,NUCROU, IPSAUV, IFICLE, IPREFI, IREFOR, ISAFOR
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, RUZG, RETG, GAMG, RHTG, REC
     &     , ROD, RUXD, RUYD, RUZD, RETD, GAMD, RHTD
     &     , UXG, UYG, UZG, UNG, UTG, UVG, PG, VOLG, DIAMG, SURF
     &     , UXD, UYD, UZD, UND, UTD, UVD, PD, VOLD, DIAMD
     &     , 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
      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, NEWLI.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     SEGACT MPTLI*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 kfm13.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 kfm13.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=5
      I2=NCEN
      SEGINI Q1
      SEGINI Q2
C
C**** On initialise le segment de travail D1
C
      I1=5
      SEGINI D1
      SEGINI NEWLI
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 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)=MPGN.VPOCHA(ICEN,3)
         CONS1.VAL(5,ICEN)=MPRETN.VPOCHA(ICEN,1)
         CONS1.VAL(6,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,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)
 
            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)
            RUZG = CONS2.VAL(4,NLCEG)
            RETG = CONS2.VAL(5,NLCEG)
            GAMG = CONS2.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
            REC  = 0.5D0 * (RUXG**2 + RUYG**2 + RUZG**2)
            REC  = REC / ROG
            PG   = (GAMG - 1.0D0) * (RETG - REC)
C
            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)
            RUZD = CONS2.VAL(4,NLCED)
            RETD = CONS2.VAL(5,NLCED)
            GAMD = CONS2.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
            REC  = 0.5D0 * (RUXD**2 + RUYD**2 + RUZD**2)
            REC  = REC / ROD
            PD   = (GAMD - 1.0D0) * (RETD - REC)
C
            VOLD  = MPOVOL.VPOCHA(NLCED,1)
            DIAMD = MPOVDI.VPOCHA(NLCED,1)
 
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))
                  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
C********** Entahlpy
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**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
               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************* 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)*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
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)
             FRUZ = (ROD * UND * UZD) + (PD * CNZ)
            FRET = (UND * RHTD)
C
C
            COEF = 0.5D0 * SURF / VOLG
C
            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 *  FRUZ)
            BN1.VAL(5,NLCEG) = BN1.VAL(5,NLCEG)
     &                       + (COEF *  FRET)
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 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 = 0.5D0 * SURF * SIGMA0.VALVEC(NLCF) / 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 *  RUZD)
            CN1.VAL(5,NLCEG) = CN1.VAL(5,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)
                 FRUZ = (ROG * UNG * UZG) + (PG * CNZ)
               FRET = (UNG * RHTG)
C
               COEF = 0.5D0 * SURF / 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 *  FRUZ)
               BN1.VAL(5,NLCED) = BN1.VAL(5,NLCED)
     &                          - (COEF *  FRET)
C
C              Viscous contribution
C
               COEF = SURF * SIGMAV.VALVEC(NLCF) / VOLD
C
               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 *  RUZG)
               BN1.VAL(5,NLCED) = BN1.VAL(5,NLCED)
     &                          - (COEF *  RETG)
C
C              Numerical diffusity contribution
C
               COEF = 0.5D0 * SURF * SIGMA0.VALVEC(NLCF) / VOLD
C
               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 *  RUZG)
               CN1.VAL(5,NLCED) = CN1.VAL(5,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)
                 RUZG = CONS2.VAL(4,ICEN)
                 RETG = CONS2.VAL(5,ICEN)
               GAMG = CONS2.VAL(6,ICEN)
C
               UXG  = RUXG / ROG
               UYG  = RUYG / ROG
                 UZG  = RUZG / ROG
C
               REC  = 0.5D0 * (RUXG**2 + RUYG**2 + RUZG**2)
               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 * 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************* Ebd of the computation of AN0 and DN0
C
               AN0.VALVEC(ICEN) = AN0.VALVEC(ICEN)
     &                          + USDTI.VALVEC(ICEN) * 2.0D0 / DCFL
C
               DN0.VALVEC(ICEN) = DN0.VALVEC(ICEN)
     &                          + (1.0D0 / DTPS)
C
C
C************* We save BN1, CN1 into BN0 and CN0
C
               DO ICOMP=1,5,1
                  BN0.VAL(ICOMP,ICEN) = BN1.VAL(ICOMP,ICEN)
               ENDDO
C
               DO ICOMP=1,5,1
                  CN0.VAL(ICOMP,ICEN) = CN1.VAL(ICOMP,ICEN)
               ENDDO
C
            ENDIF
C
C********** IBLJAC => 1
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********** D2 = \delta \mathbf{U}
C
            DO ICOMP=1,5,1
               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
C**** Exit Data
C
      DO ICEN=1,NCEN,1
         DO ICOMP=1,5,1
            MPODU.VPOCHA(ICEN,ICOMP) = DU.VAL(ICOMP,ICEN)
         ENDDO
      ENDDO
C
C
C
      SEGSUP DU
      SEGSUP RES
      SEGSUP CONS1
      SEGSUP CONS2
      SEGSUP Q1
      SEGSUP Q2
      SEGSUP D1
      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
      SEGSUP NEWLI
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