Télécharger pre62f.eso

Retour à la liste

pre62f
C PRE62F    SOURCE    CB215821  20/11/25    13:36:36     10792          
      SUBROUTINE PRE62F(LOGTEM,
     &                  ICEN,IFACE,IFACEL,INORM,
     &                  IALPH, IGRALP, IALALP,
     &                  IUVC, IGRUVC, IALUVC,
     &                  IULC, IGRULC, IALULC,
     &                  IPC,  IGRPC,  IALPC,
     &                  ITVC, IGRTVC, IALTVC,
     &                  ITLC, IGRTLC, IALTLC,
     &                  IRVC, IRLC,
     &                  DELTAT,
     &                  IALPHF, IUVF, IULF, IPF, ITVF, ITLF,
     &                  IRVF, IRLF,
     &                  LOGAN,LOGNEG,LOGBOR,MESERR,VALER,VAL1,VAL2)
C************************************************************************
C
C PROJET            :  CASTEM 2000
C
C NOM               :  PRE62F
C
C DESCRIPTION       :  Voir PRE42F
C
C LANGAGE           :  FORTRAN 77 + ESOPE 2000 (avec extensions CISI)
C
C AUTEUR            :  A. BECCANTINI, DRN/DMT/SEMT/TTMF
C                      Modified for two-fluid flow by
C                      Jose R. Garcia-Cascales
C
C************************************************************************
C
C
C APPELES (Outils)  :  KRIPAD, LICHT
C
C APPELES (Calcul)  :  AUCUN
C
C
C************************************************************************
C
C ENTREES
C
C     LOGTEM  : LOGICAL; si .TRUE. 2em ordre en temps
C                        sinon 1er ordre en temps
C
C     1) Pointeurs de MELEMEs et de CHPOINTs de la table DOMAINE
C
C     ICEN    : MELEME de 'POI1' SPG des CENTRES
C
C     IFACE   : MELEME de 'POI1' SPG des FACES
C
C     IFACEL  : MELEME de 'SEG3' avec
C               CENTRE d'Elt "gauche"
C               CENTRE de Face
C               CENTRE d'Elt "droite"
C
C               N.B. = IFACE.NUM(i,1) = IFACEL.NUM(i,2)
C
C     INORM   : CHPOINT des cosinus directeurs de normales aux faces
C
C     2) Pointeurs des  CHPOINTs
C
C     IALPH   : CHPOINT "CENTRE" containing  void fraction
C
C     IGRALP  : CHPOINT "CENTRE" contaninig the void fraction gradient
C               (2 composantes)
C
C     IALALP  : CHPOINT "CENTRE" contaning the limiter of the
C               void fraction gradient
C
C     IUVC    : CHPOINT "CENTRE" containing the vapour velocity UvX, UVY ;
C
C     IGRUVC  : CHPOINT "CENTRE" containig the vapour velocity gradient
C               (4 composantes)
C
C     IALUVC  : CHPOINT "CENTRE" contaning the limiter of the vapour
C               velocity (2 composantes)
C
C     IULC    : CHPOINT "CENTRE" containing the liquid  velocity UvX, UVY ;
C
C     IGRULC  : CHPOINT "CENTRE" containig the liquid velocity gradient
C               (4 composantes)
C
C     IALULC  : CHPOINT "CENTRE" contaning the limiter of the liquid
C               velocity (2 composantes)
C
C     IPC     : CHPOINT "CENTRE" containing the pressure P ;
C
C     IGRPC   : CHPOINT "CENTRE" containing the gradient of pressure
C               (2 composantes)
C
C     IALPC   : CHPOINT "CENTRE" containig the limiter of pressure
C               gradient
C
C     ITVC     : CHPOINT "CENTRE" containing the vapour temperature Tv ;
C
C     IGRTVC   : CHPOINT "CENTRE" containing the gradient of vapour
C                temperature (2 composantes)
C
C     IALTVC   : CHPOINT "CENTRE" containig the limiter of vapour
C                temperature gradient
C
C     ITLC     : CHPOINT "CENTRE" containing the liquid temperature Tl ;
C
C     IGRTLC   : CHPOINT "CENTRE" containing the gradient of liquid
C                temperature (2 composantes)
C
C     IALTLC   : CHPOINT "CENTRE" containig the limiter of liquid
C                temperature gradient
C
C     IRVC     : CHPOINT "CENTRE" containing the vapour density
C
C     IRLC     : CHPOINT "CENTRE" containing the liquid density
C
C     3)
C
C     DELTAT  : REAL*8, encrement en temps pour calculer la prediction
C
C
C SORTIES
C
C
C     IROF    : MCHAML  defini sur le MELEME de pointeur IFACEL,
C               contenant la masse volumique RHO
C
C     IVITF   : MCHAML  defini sur le MELEME de pointeur IFACEL,
C               contenant la vitesse UN, UT dans le repaire local
C               (n,t) et defini sur le MELEME de pointeur IFACE,
C               contenant les cosinus directeurs du repere local
C
C     IPF     : MCHAML  defini sur le MELEME de pointeur IFACEL,
C               contenant la pression P
C
C     IGAMF   : MCHAML  defini sur le MELEME de pointeur IFACEL,
C               contenant le "gamma" du gaz
C
C     LOGAN   : anomalie detectee
C
C     LOGNEG  : (LOGICAL): si .TRUE. une pression ou une densité
C               negative a été detectée -> en interactif le
C               programme s'arrete  en GIBIANE
C               (erreur stocké en MESERR et VALER)
C
C     LOGBOR  : (LOGICAL): si .TRUE. un gamma a ete detecte
C               dehor 1 et 3 (sa valeur stockée en MESERR et VALER;
C               en VAL1 et en VAL2 on stocke 1.0 et 3.0)
C
C     MESERR
C     VALER
C     VAL1,
C     VAL2    : pour les messages d'erreur
C
C************************************************************************
C
C HISTORIQUE (Anomalies et modifications éventuelles)
C
C HISTORIQUE :  Créée le 02/05/02.
C
C************************************************************************
C
C
C ATTENTION: Cet programme marche si le MAILLAGE est convex;
C            si non il faut changer l'argoritme de calcul de
C            l'orientation des normales aux faces.
C
C            La positivité n'est pas controlle parce que c'est déjà fait
C            dans l'operateur PRIM
C
C
C************************************************************************
C
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8(A-H,O-Z)
 
C
C**** Les variables
C
      INTEGER ICEN, IFACE, IFACEL, INORM,
     &        IALPH, IGRALP, IALALP,
     &        IUVC,  IGRUVC, IALUVC,
     &        IULC,  IGRULC, IALULC,
     &        IPC,   IGRPC,  IALPC,
     &        ITVC,  IGRTVC, IALTVC,
     &        ITLC,  IGRTLC, IALTLC,
     &        IRVC,  IRLC,
     &        IALPHF, IUVF, IULF, IPF, ITVF, ITLF,
     &        IRVF, IRLF,
     &        IGEOM, NFAC, NCEN,
     &        N1PTEL, N1EL, N2PTEL, N2EL, N2, N1, N3, L1, NLCE,
     &        NLCF, NGCF, NGCEG, NLCEG, NGCED, NLCED, NGCF1,
     &        IDIMP1, INDCEL,
     &        l, m, n
      REAL*8 VALER, VAL1, VAL2, XG, YG, ZG, XC, YC, ZC, XD, YD, ZD,
     &        DELTAT
     &       ,DXG, DYG, DZG, DXD, DYD, DZD,ORIENT
     &       , CNX, CNY, CNZ, CTX, CTY, CTZ, CVX, CVY, CVZ
     &       , alphag, uvxg, uvyg, uvzg, uvng, uvtg, uvvg
     &       ,         ulxg, ulyg, ulzg, ulng, ultg, ulvg
     &       , pg, Tvg, Tlg, rvg, rlg
     &       , alphad, uvxd, uvyd, uvzd, uvnd, uvtd, uvvd
     &       ,         ulxd, ulyd, ulzd, ulnd, ultd, ulvd
     &       , pd, Tvd, Tld, rvd, rld
     &       , VALCEL, DCEL, ALCEL
     &       , alpha, uvx, uvy, uvz, ulx, uly, ulz, p, Tv, Tl
     &       , rv, rl, drvdp, drvdhv, drldp, drldhl
     &       , hv, hl
     &       , A(10,10), B(10,10), C(10, 10), denom
     &       , dVdx(10), dVdy(10), dVdz(10), dV(10), LIMITA(10)
     &       , dVdVp(10,10), dVpdV(10,10), P1(10,10)
     &       , dhvdp, dhvdTv, dhldp, dhldTl
C     &       , DROX, DROY, DROZ, DUXX, DUXY, DUXZ, DUYX, DUYY,
C     &         DUYZ, DUZX, DUZY, DUZZ, DPX, DPY, DPZ
C     &       , DRO, DUX, DUY, DUZ, DP, ALPHA(5)
      CHARACTER*(40) MESERR
      CHARACTER*(8) TYPE
      LOGICAL LOGAN,LOGNEG, LOGBOR, LOGTEM, LOGI1, LOGALP
C
C**** perfect and stiffened gas eos
C
      real*8 gamv, cpv, gaml, cpl, pil
      PARAMETER(gamv = 1.4D0, cpv = 1008.D0,
     &          gaml = 2.8D0, cpl = 4186.D0, pil = 8.5D8)
C
C**** LOGALP = .TRUE. -> Prediction avec limiteur
C
      PARAMETER(LOGALP = .TRUE.)
C
C**** Les Includes
C
-INC SMCOORD
 
-INC PPARAM
-INC CCOPTIO
-INC SMCHPOI
      POINTEUR MPALPC.MPOVAL, MPGRAL.MPOVAL,
     &         MPUVC.MPOVAL, MPGRUV.MPOVAL,
     &         MPULC.MPOVAL, MPGRUL.MPOVAL,
     &         MPPC.MPOVAL,  MPGRP.MPOVAL,
     &         MPTVC.MPOVAL, MPGRTV.MPOVAL,
     &         MPTLC.MPOVAL, MPGRTL.MPOVAL,
     &         MPRVC.MPOVAL, MPRLC.MPOVAL,
     &         MPNORM.MPOVAL,
     &         MPALPP.MPOVAL, MPUVP.MPOVAL, MPULP.MPOVAL,
     &         MPPP.MPOVAL,   MPTVP.MPOVAL, MPTLP.MPOVAL,
     &         MPRVP.MPOVAL,  MPRLP.MPOVAL
-INC SMCHAML
      POINTEUR MLALP.MELVAL, MLP.MELVAL,
     &         MLTV.MELVAL, MLTL.MELVAL,
     &         MLRV.MELVAL, MLRL.MELVAL
      POINTEUR MLUVN.MELVAL,  MLUVT.MELVAL,  MLUVV.MELVAL,
     &         MLULN.MELVAL,  MLULT.MELVAL,  MLULV.MELVAL,
     &         MLVNX.MELVAL,  MLVNY.MELVAL,  MLVNZ.MELVAL,
     &         MLVT1X.MELVAL, MLVT1Y.MELVAL, MLVT1Z.MELVAL,
     &         MLVT2X.MELVAL, MLVT2Y.MELVAL, MLVT2Z.MELVAL,
     &         MLLNX.MELVAL,  MLLNY.MELVAL,  MLLNZ.MELVAL,
     &         MLLT1X.MELVAL, MLLT1Y.MELVAL, MLLT1Z.MELVAL,
     &         MLLT2X.MELVAL, MLLT2Y.MELVAL, MLLT2Z.MELVAL
C**** Some melemes for the LIMITAs
      POINTEUR MELLAL.MPOVAL,
     &         MELLUV.MPOVAL,
     &         MELLUL.MPOVAL,
     &         MELLP.MPOVAL,
     &         MELLTV.MPOVAL,
     &         MELLTL.MPOVAL
-INC SMLENTI
-INC SMELEME
C
C
C**** Initialisation des parametres d'erreur déjà faite, i.e.
C
C     LOGNEG  = .FALSE.
C     LOGBOR  = .FALSE.
C     MESERR = '                                        '
C     MOTERR(1:40) = MESERR(1:40)
C     VALER = 0.0D0
C     VAL1 = 0.0D0
C     VAL2 = 0.0D0
C
C
C**** KRIPAD pour la correspondance global/local de centre
C
      CALL KRIPAD(ICEN,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     SEGACT MLENT1
C
C**** Activation de CHPOINTs
C
C     void fraction + grad + limiter
C     vapour velocity  + grad + limiter
C     liquid velocity  + grad + limiter
C     pressure + grad + limiter
C     vapour temperature + grad + limiter
C     liquid temperature + grad + limiter
C     vapour density
C     liquid density
C
C     cosinus directeurs des normales aux surface
C
      CALL LICHT(IALPH, MPALPC, TYPE, IGEOM)
      CALL LICHT(IGRALP, MPGRAL, TYPE, IGEOM)
      CALL LICHT(IUVC, MPUVC, TYPE, IGEOM)
      CALL LICHT(IGRUVC, MPGRUV, TYPE, IGEOM)
      CALL LICHT(IULC, MPULC, TYPE, IGEOM)
      CALL LICHT(IGRULC, MPGRUL, TYPE, IGEOM)
      CALL LICHT(IPC , MPPC , TYPE, IGEOM)
      CALL LICHT(IGRPC, MPGRP, TYPE, IGEOM)
      CALL LICHT(ITVC, MPTVC, TYPE, IGEOM)
      CALL LICHT(IGRTVC, MPGRTV, TYPE, IGEOM)
      CALL LICHT(ITLC, MPTLC, TYPE, IGEOM)
      CALL LICHT(IGRTLC, MPGRTL, TYPE, IGEOM)
      CALL LICHT(IRVC, MPRVC, TYPE, IGEOM)
      CALL LICHT(IRLC, MPRLC, TYPE, IGEOM)
      CALL LICHT(INORM, MPNORM, TYPE, IGEOM)
C
C**** Les MPOVALs 'Prediction'
C
      IF(LOGTEM)THEN
         SEGINI, MPALPP = MPALPC
         SEGINI, MPUVP = MPUVC
         SEGINI, MPULP = MPULC
         SEGINI, MPPP  = MPPC
         SEGINI, MPTVP = MPTVC
         SEGINI, MPTLP = MPTLC
         SEGINI, MPRVP = MPRVC
         SEGINI, MPRLP = MPRLC
      ELSE
         MPALPP = MPALPC
         MPUVP  = MPUVC
         MPULP  = MPULC
         MPPP   = MPPC
         MPTVP  = MPTVC
         MPTLP  = MPTLC
         MPRVP  = MPRVC
         MPRLP  = MPRLC
      ENDIF
C
C**** Les Limiteurs
C
      CALL LICHT(IALALP, MELLAL, TYPE, IGEOM)
      CALL LICHT(IALUVC, MELLUV, TYPE, IGEOM)
      CALL LICHT(IALULC, MELLUL, TYPE, IGEOM)
      CALL LICHT(IALPC , MELLP, TYPE, IGEOM)
      CALL LICHT(IALTVC, MELLTV, TYPE, IGEOM)
      CALL LICHT(IALTLC, MELLTL, TYPE, IGEOM)
C
C
C**** MPOVAL sont déjà activés i.e.:
C
C     SEGACT MPALPP
C     SEGACT MPGRAL
C     SEGACT IALALP
C     SEGACT MPUVC
C     SEGACT IALUVC
C     SEGACT MPGRUV
C     SEGACT MPULC
C     SEGACT MPGRUL
C     SEGACT IALULC
C     SEGACT MPPC
C     SEGACT MPGRP
C     SEGACT IALPC
C     SEGACT MPTVC
C     SEGACT MPGRTV
C     SEGACT IALTVC
C     SEGACT MPTLC
C     SEGACT MPGRTL
C     SEGACT IALTLC
C     SEGACT MPRVC
C     SEGACT MPRLC
C     SEGACT MPNORM
C
C**** Le MELEME FACEL
C
      IPT1 = IFACEL
      IPT2 = IFACE
      SEGACT IPT1
      SEGACT IPT2
      NFAC = IPT1.NUM(/2)
C
C**** Creation de MCHAMLs contenant les etats gauche et droite,
C
C     i.e.:
C
C     vitesse + cosinus directors du repere local
C     densité
C     pression
C     gamma
C
C**** Cosinus directors du repere local et vitesse
C
C     Les cosinus directeurs
C
      N1 = 2
      N3 = 6
      L1 = 28
      SEGINI MCHEL1
      IUVF = MCHEL1
      MCHEL1.TITCHE = 'UV '
      MCHEL1.IMACHE(1) = IFACE
      MCHEL1.IMACHE(2) = IFACEL
      MCHEL1.CONCHE(1) = '(n,t,v)in(x,y,z)'
      MCHEL1.CONCHE(2) = ' UV in (n,t,v)  '
C
C**** Valeurs des cosinus definies par respect au repair global, i.e.
C
      MCHEL1.INFCHE(1,1) = 2
      MCHEL1.INFCHE(1,3) = NIFOUR
      MCHEL1.INFCHE(1,4) = 0
      MCHEL1.INFCHE(1,5) = 0
      MCHEL1.INFCHE(1,6) = 0
      MCHEL1.IFOCHE = IFOUR
C
C**** Valeurs de vitesse definies par respect au repair local, i.e.
C
      MCHEL1.INFCHE(2,1) = 1
      MCHEL1.INFCHE(2,3) = NIFOUR
      MCHEL1.INFCHE(2,4) = 0
      MCHEL1.INFCHE(2,5) = 0
      MCHEL1.INFCHE(2,6) = 0
C
C**** Le cosinus directeurs
C
      N1PTEL = 1
      N1EL = NFAC
      N2PTEL = 0
      N2EL =   0
C
C**** MCHAML a N2 composantes:
C
C     cosinus directeurs du repere local (n,t1,t2)
C
C     IDIM = 3 -> 9 composantes
C
      N2 = 9
      SEGINI MCHAM1
      MCHEL1.ICHAML(1) = MCHAM1
      MCHAM1.NOMCHE(1) = 'NVX     '
      MCHAM1.NOMCHE(2) = 'NVY     '
      MCHAM1.NOMCHE(3) = 'NVZ     '
      MCHAM1.NOMCHE(4) = 'TVX     '
      MCHAM1.NOMCHE(5) = 'TVY     '
      MCHAM1.NOMCHE(6) = 'TVZ     '
      MCHAM1.NOMCHE(7) = 'VVX     '
      MCHAM1.NOMCHE(8) = 'VVY     '
      MCHAM1.NOMCHE(9) = 'VVZ     '
      MCHAM1.TYPCHE(1) = 'REAL*8          '
      MCHAM1.TYPCHE(2) = 'REAL*8          '
      MCHAM1.TYPCHE(3) = 'REAL*8          '
      MCHAM1.TYPCHE(4) = 'REAL*8          '
      MCHAM1.TYPCHE(5) = 'REAL*8          '
      MCHAM1.TYPCHE(6) = 'REAL*8          '
      MCHAM1.TYPCHE(7) = 'REAL*8          '
      MCHAM1.TYPCHE(8) = 'REAL*8          '
      MCHAM1.TYPCHE(9) = 'REAL*8          '
      SEGINI MLVNX
      SEGINI MLVNY
      SEGINI MLVNZ
      SEGINI MLVT1X
      SEGINI MLVT1Y
      SEGINI MLVT1Z
      SEGINI MLVT2X
      SEGINI MLVT2Y
      SEGINI MLVT2Z
      MCHAM1.IELVAL(1) = MLVNX
      MCHAM1.IELVAL(2) = MLVNY
      MCHAM1.IELVAL(3) = MLVNZ
      MCHAM1.IELVAL(4) = MLVT1X
      MCHAM1.IELVAL(5) = MLVT1Y
      MCHAM1.IELVAL(6) = MLVT1Z
      MCHAM1.IELVAL(7) = MLVT2X
      MCHAM1.IELVAL(8) = MLVT2Y
      MCHAM1.IELVAL(9) = MLVT2Z
      SEGDES MCHAM1
C
C**** Vitesse
C
      N1EL = NFAC
      N1PTEL = 3
      N2EL = 0
      N2PTEL = 0
C
C**** MCHAML a N2 composantes:
C
C     IDIM = 3 -> 3 composantes
C
      N2 = 3
      SEGINI MCHAM1
      MCHEL1.ICHAML(2) = MCHAM1
      SEGDES MCHEL1
      MCHAM1.NOMCHE(1) = 'UVN     '
      MCHAM1.NOMCHE(2) = 'UVT     '
      MCHAM1.NOMCHE(3) = 'UVV     '
      MCHAM1.TYPCHE(1) = 'REAL*8          '
      MCHAM1.TYPCHE(2) = 'REAL*8          '
      MCHAM1.TYPCHE(3) = 'REAL*8          '
      SEGINI MLUVN
      SEGINI MLUVT
      SEGINI MLUVV
      MCHAM1.IELVAL(1) = MLUVN
      MCHAM1.IELVAL(2) = MLUVT
      MCHAM1.IELVAL(3) = MLUVV
      SEGDES MCHAM1
C
C**** Cosinus directors du repere local et vitesse
C
C     Les cosinus directeurs
C
C     LIQUID PHASE
 
      N1 = 2
      N3 = 6
      L1 = 28
      SEGINI MCHEL1
      IULF = MCHEL1
      MCHEL1.TITCHE = 'UL '
      MCHEL1.IMACHE(1) = IFACE
      MCHEL1.IMACHE(2) = IFACEL
      MCHEL1.CONCHE(1) = '(n,t,v) in (x,y)'
      MCHEL1.CONCHE(2) = ' UL in (n,t,v)  '
C
C**** Valeurs des cosinus definies par respect au repair global, i.e.
C
      MCHEL1.INFCHE(1,1) = 2
      MCHEL1.INFCHE(1,3) = NIFOUR
      MCHEL1.INFCHE(1,4) = 0
      MCHEL1.INFCHE(1,5) = 0
      MCHEL1.INFCHE(1,6) = 0
      MCHEL1.IFOCHE = IFOUR
C
C**** Valeurs de vitesse definies par respect au repair local, i.e.
C
      MCHEL1.INFCHE(2,1) = 1
      MCHEL1.INFCHE(2,3) = NIFOUR
      MCHEL1.INFCHE(2,4) = 0
      MCHEL1.INFCHE(2,5) = 0
      MCHEL1.INFCHE(2,6) = 0
C
C**** Le cosinus directeurs
C
      N1PTEL = 1
      N1EL = NFAC
      N2PTEL = 0
      N2EL =   0
C
C**** MCHAML a N2 composantes:
C
C     cosinus directeurs du repere local (n,t1)
C
C     IDIM = 3 -> 9 composantes
C
      N2 = 9
      SEGINI MCHAM1
      MCHEL1.ICHAML(1) = MCHAM1
      MCHAM1.NOMCHE(1) = 'NLX     '
      MCHAM1.NOMCHE(2) = 'NLY     '
      MCHAM1.NOMCHE(3) = 'NLZ     '
      MCHAM1.NOMCHE(4) = 'TLX     '
      MCHAM1.NOMCHE(5) = 'TLY     '
      MCHAM1.NOMCHE(6) = 'TLZ     '
      MCHAM1.NOMCHE(7) = 'VLX     '
      MCHAM1.NOMCHE(8) = 'VLY     '
      MCHAM1.NOMCHE(9) = 'VLZ     '
      MCHAM1.TYPCHE(1) = 'REAL*8          '
      MCHAM1.TYPCHE(2) = 'REAL*8          '
      MCHAM1.TYPCHE(3) = 'REAL*8          '
      MCHAM1.TYPCHE(4) = 'REAL*8          '
      MCHAM1.TYPCHE(5) = 'REAL*8          '
      MCHAM1.TYPCHE(6) = 'REAL*8          '
      MCHAM1.TYPCHE(7) = 'REAL*8          '
      MCHAM1.TYPCHE(8) = 'REAL*8          '
      MCHAM1.TYPCHE(9) = 'REAL*8          '
      SEGINI MLLNX
      SEGINI MLLNY
      SEGINI MLLNZ
      SEGINI MLLT1X
      SEGINI MLLT1Y
      SEGINI MLLT1Z
      SEGINI MLLT2X
      SEGINI MLLT2Y
      SEGINI MLLT2Z
      MCHAM1.IELVAL(1) = MLLNX
      MCHAM1.IELVAL(2) = MLLNY
      MCHAM1.IELVAL(3) = MLLNZ
      MCHAM1.IELVAL(4) = MLLT1X
      MCHAM1.IELVAL(5) = MLLT1Y
      MCHAM1.IELVAL(6) = MLLT1Z
      MCHAM1.IELVAL(7) = MLLT2X
      MCHAM1.IELVAL(8) = MLLT2Y
      MCHAM1.IELVAL(9) = MLLT2Z
      SEGDES MCHAM1
C
C**** Vitesse
C
      N1EL = NFAC
      N1PTEL = 3
      N2EL = 0
      N2PTEL = 0
C
C**** MCHAML a N2 composantes:
C
C     IDIM = 3 -> 3 composantes
C
      N2 = 3
      SEGINI MCHAM1
      MCHEL1.ICHAML(2) = MCHAM1
      SEGDES MCHEL1
      MCHAM1.NOMCHE(1) = 'ULN     '
      MCHAM1.NOMCHE(2) = 'ULT     '
      MCHAM1.NOMCHE(3) = 'ULV     '
      MCHAM1.TYPCHE(1) = 'REAL*8         '
      MCHAM1.TYPCHE(2) = 'REAL*8         '
      MCHAM1.TYPCHE(3) = 'REAL*8         '
      SEGINI MLULN
      SEGINI MLULT
      SEGINI MLULV
      MCHAM1.IELVAL(1) = MLULN
      MCHAM1.IELVAL(2) = MLULT
      MCHAM1.IELVAL(3) = MLULV
      SEGDES MCHAM1
C
C**** Void fraction
C
      N1 = 1
      N3 = 6
      L1 = 15
      SEGINI MCHEL2
      IALPHF = MCHEL2
      MCHEL2.IMACHE(1) = IFACEL
      MCHEL2.TITCHE = 'ALPHA           '
      MCHEL2.CONCHE(1) = '                '
C
C**** Valeurs independente du repére, i.e.
C
      MCHEL2.INFCHE(1,1) = 0
      MCHEL2.INFCHE(1,3) = NIFOUR
      MCHEL2.INFCHE(1,4) = 0
      MCHEL2.INFCHE(1,5) = 0
      MCHEL2.INFCHE(1,6) = 0
      MCHEL2.IFOCHE = IFOUR
      N2 = 1
      SEGINI MCHAM1
      MCHEL2.ICHAML(1) = MCHAM1
      SEGDES MCHEL2
      MCHAM1.NOMCHE(1) = 'SCAL    '
      MCHAM1.TYPCHE(1) = 'REAL*8          '
      SEGINI MLALP
      MCHAM1.IELVAL(1) = MLALP
      SEGDES MCHAM1
C
C**** Pressure
C
      MCHEL1 = IALPHF
      SEGINI, MCHEL2 = MCHEL1
      IPF = MCHEL2
      MCHEL2.TITCHE = 'P            '
C
C**** MCHAM1 = MCHAML de la alpha
C
      SEGINI, MCHAM2 = MCHAM1
      MCHEL2.ICHAML(1) = MCHAM2
      SEGDES MCHEL2
      SEGINI MLP
      MCHAM2.IELVAL(1) = MLP
      SEGDES MCHAM2
C
C**** Vapour temperature
C
      MCHEL1 = IALPHF
      SEGINI, MCHEL2 = MCHEL1
      ITVF = MCHEL2
      MCHEL2.TITCHE = 'TV             '
C
C**** MCHAM1 = MCHAML de la alpha
C
      SEGINI, MCHAM2 = MCHAM1
      MCHEL2.ICHAML(1) = MCHAM2
      SEGDES MCHEL2
      SEGINI MLTV
      MCHAM2.IELVAL(1) = MLTV
      SEGDES MCHAM2
C
C**** Liquid temperature
C
      MCHEL1 = IALPHF
      SEGINI, MCHEL2 = MCHEL1
      ITLF = MCHEL2
      MCHEL2.TITCHE = 'TL             '
C
C**** MCHAM1 = MCHAML de la alpha
C
      SEGINI, MCHAM2 = MCHAM1
      MCHEL2.ICHAML(1) = MCHAM2
      SEGDES MCHEL2
      SEGINI MLTL
      MCHAM2.IELVAL(1) = MLTL
      SEGDES MCHAM2
C
C**** Vapour density
C
      MCHEL1 = IALPHF
      SEGINI, MCHEL2 = MCHEL1
      IRVF = MCHEL2
      MCHEL2.TITCHE = 'RV             '
C
C**** MCHAM1 = MCHAML de la alpha
C
      SEGINI, MCHAM2 = MCHAM1
      MCHEL2.ICHAML(1) = MCHAM2
      SEGDES MCHEL2
      SEGINI MLRV
      MCHAM2.IELVAL(1) = MLRV
      SEGDES MCHAM2
C
C**** Liquid density
C
      MCHEL1 = IALPHF
      SEGINI, MCHEL2 = MCHEL1
      IRLF = MCHEL2
      MCHEL2.TITCHE = 'RL             '
C
C**** MCHAM1 = MCHAML de la alpha
C
      SEGINI, MCHAM2 = MCHAM1
      MCHEL2.ICHAML(1) = MCHAM2
      SEGDES MCHEL2
      SEGINI MLRL
      MCHAM2.IELVAL(1) = MLRL
      SEGDES MCHAM2
C
C**** Recapitulatif: les MELVALs et les MPOVALs actives
C
C     MLVNX, MLVNY, MLVNZ
C     MLVTX, MLVTY, MLVTZ
C     MLVVX, MLVVY, MLVVZ
C
C     MLLNX, MLLNY, MLLNZ
C     MLLTX, MLLTY, MLLTZ
C     MLLVX, MLLVY, MLLVZ
C
C     MLUVN, MLUVT, MLUVV -> vapour velocities
C
C     MLULN, MLULT, MLULV -> liquid velocities
C
C     MLALP  -> void fraction
C
C     MLP     -> pressure
C
C     MLTV    -> vapour temperature
C
C     MLTL    -> liquid temperature
C
C     MLRV    -> vapour density
C
C     MLRL    -> liquid density
C****
C     MPALP   -> void fraction
C
C     MPUVC   -> vapour velocity
C
C     MPULC   -> liquid velocity
C
C     MPPC    -> pressure
C
C     MPTVC   -> vapour temperature
C
C     MPTLC   -> liquid temperature
C
C     MPRVC   -> vapour density
C
C     MPRLC   -> liquid density
C
C     MPNORM  -> normales aux faces
C
C
C***********************************************************************
C********* PREDICTION **************************************************
C***********************************************************************
C
      IF(LOGTEM)THEN
C
         IPT3 = ICEN
         SEGACT IPT3
         NCEN = IPT3.NUM(/2)
         DO NLCE = 1, NCEN
            IF(LOGALP)THEN
C
C************* Gradients limités
C
               LIMITA(1) = MELLAL.VPOCHA(NLCE,1)
               LIMITA(2) = MELLUV.VPOCHA(NLCE,1)
               LIMITA(3) = MELLUV.VPOCHA(NLCE,2)
               LIMITA(4) = MELLUV.VPOCHA(NLCE,3)
               LIMITA(5) = MELLUL.VPOCHA(NLCE,1)
               LIMITA(6) = MELLUL.VPOCHA(NLCE,2)
               LIMITA(7) = MELLUL.VPOCHA(NLCE,3)
               LIMITA(8) = MELLP.VPOCHA(NLCE,1)
               LIMITA(9) = MELLTV.VPOCHA(NLCE,1)
               LIMITA(10) = MELLTL.VPOCHA(NLCE,1)
 
            ELSE
C
C************* Gradients non limités
C
               LIMITA(1) = 1.D0
               LIMITA(2) = 1.D0
               LIMITA(3) = 1.D0
               LIMITA(4) = 1.D0
               LIMITA(5) = 1.D0
               LIMITA(6) = 1.D0
               LIMITA(7) = 1.D0
               LIMITA(8) = 1.D0
               LIMITA(9) = 1.D0
               LIMITA(10) = 1.D0
 
            ENDIF
 
            alpha = MPALPP.VPOCHA(NLCE,1)
            uvx = MPUVP.VPOCHA(NLCE,1)
            uvy = MPUVP.VPOCHA(NLCE,2)
            uvz = MPUVP.VPOCHA(NLCE,3)
            ulx = MPULP.VPOCHA(NLCE,1)
            uly = MPULP.VPOCHA(NLCE,2)
            ulz = MPULP.VPOCHA(NLCE,3)
            p   = MPPP.VPOCHA(NLCE,1)
            Tv  = MPTVP.VPOCHA(NLCE,1)
            Tl  = MPTLP.VPOCHA(NLCE,1)
            rv  = MPRVP.VPOCHA(NLCE,1)
            rl  = MPRLP.VPOCHA(NLCE,1)
C
C************* Value of the enthalpies and the derivatives of
C              the densities respect pressure and enthalpies,
C              for the moment only for perfect and stiffened gases.
            hv = cpv*Tv
            hl = cpl*Tl
            drvdp  =  rv/p
            drvdhv = -rv/hv
            drldp  = rl/(p + pil)
            drldhl = -rl/hl
            dhvdp = 0.d0
            dhvdTv = cpv
            dhldp = 0.d0
            dhldTl = cpl
 
      denom = drldhl*(drvdhv*p + rv*(alpha*drvdp*p + rv - alpha*rv)) +
     &     rl*(alpha*rl*(drvdhv + drvdp*rv) - (-1 + alpha)*drldp*(drvdhv
     $     *p + rv**2))
 
       A(1,1) = (alpha*(drldhl*p + rl**2)*(drvdhv + drvdp*rv)*ulx -
     &     (-1.d0 + alpha)*(drldhl + drldp*rl)*(drvdhv*p + rv**2)*uvx)
     $     /denom
       A(1,2) = -(((-1.d0 + alpha)*alpha*(drldhl + drldp*rl)*rv**2)
     $      /denom)
       A(1,3) = 0.d0
       A(1,4) = 0.d0
       A(1,5) = ((-1.d0 + alpha)*alpha*rl**2*(drvdhv + drvdp*rv))/denom
       A(1,6) = 0.d0
       A(1,7) = 0.d0
       A(1,8) = ((-1.d0 + alpha)*alpha*(drldhl + drldp*rl)*(drvdhv +
     $      drvdp*rv)*(ulx - uvx))/denom
       A(1,9) = 0.d0
       A(1,10) = 0.d0
 
       A(2,1) = p/(alpha*rv)
       A(2,2) = uvx
       A(2,3) = 0.d0
       A(2,4) = 0.d0
       A(2,5) = 0.d0
       A(2,6) = 0.d0
       A(2,7) = 0.d0
       A(2,8) = 1.d0/rv
       A(2,9) = 0.d0
       A(2,10) = 0.d0
 
       A(3,1) = 0.d0
       A(3,2) = 0.d0
       A(3,3) = uvx
       A(3,4) = 0.d0
       A(3,5) = 0.d0
       A(3,6) = 0.d0
       A(3,7) = 0.d0
       A(3,8) = 0.d0
       A(3,9) = 0.d0
       A(3,10) = 0.d0
 
       A(4,1) = 0.d0
       A(4,2) = 0.d0
       A(4,3) = 0d0
       A(4,4) = uvx
       A(4,5) = 0.d0
       A(4,6) = 0.d0
       A(4,7) = 0.d0
       A(4,8) = 0.d0
       A(4,9) = 0.d0
       A(4,10) = 0.d0
 
       A(5,1) = p/(-rl + alpha*rl)
       A(5,2) = 0.d0
       A(5,3) = 0.d0
       A(5,4) = 0.d0
       A(5,5) = ulx
       A(5,6) = 0.d0
       A(5,7) = 0.d0
       A(5,8) = 1.d0/rl
       A(5,9) = 0.d0
       A(5,10) = 0.d0
 
       A(6,1) = 0.d0
       A(6,2) = 0.d0
       A(6,3) = 0.d0
       A(6,4) = 0.d0
       A(6,5) = 0.d0
       A(6,6) = ulx
       A(6,7) = 0.d0
       A(6,8) = 0.d0
       A(6,9) = 0.d0
       A(6,10) = 0.d0
 
       A(7,1) = 0.d0
       A(7,2) = 0.d0
       A(7,3) = 0.d0
       A(7,4) = 0.d0
       A(7,5) = 0.d0
       A(7,6) = 0.d0
       A(7,7) = ulx
       A(7,8) = 0.d0
       A(7,9) = 0.d0
       A(7,10) = 0.d0
 
       A(8,1) = -(((drldhl*p + rl**2)*(drvdhv*p + rv**2)*(ulx - uvx))
     $      /denom)
       A(8,2) = (alpha*(drldhl*p + rl**2)*rv**2)/denom
       A(8,3) = 0.d0
       A(8,4) = 0.d0
       A(8,5) = -(((-1.d0 + alpha)*rl**2*(drvdhv*p + rv**2))/denom)
       A(8,6) = 0.d0
       A(8,7) = 0.d0
       A(8,8) = (-((-1.d0 + alpha)*(drldhl + drldp*rl)*(drvdhv*p + rv**2
     $      )*ulx) + alpha*(drldhl*p + rl**2)*(drvdhv + drvdp*rv)*uvx)
     $      /denom
       A(8,9) = 0.d0
       A(8,10) = 0.d0
 
       A(9,1) = ((drldhl*p + rl**2)*(drvdp*p - rv)*(ulx - uvx))/denom
       A(9,2) = ((drldhl*p + rl*(drldp*(p - alpha*p) + alpha*rl))*rv)
     $      /denom
       A(9,3) = 0.d0
       A(9,4) = 0.d0
       A(9,5) = ((-1.d0 + alpha)*rl**2*(drvdp*p - rv))/denom
       A(9,6) = 0.d0
       A(9,7) = 0.d0
       A(9,8) = ((-1.d0 + alpha)*(drldhl + drldp*rl)*(drvdp*p - rv)*(ulx
     $      - uvx))/denom
       A(9,9) = uvx
       A(9,10) = 0.d0
 
       A(10,1) = ((drldp*p - rl)*(drvdhv*p + rv**2)*(ulx - uvx))/denom
       A(10,2) = (alpha*(-(drldp*p) + rl)*rv**2)/denom
       A(10,3) = 0.d0
       A(10,4) = 0.d0
       A(10,5) = (rl*(drvdhv*p + rv*(alpha*drvdp*p + rv - alpha*rv)))
     $      /denom
       A(10,6) = 0.d0
       A(10,7) = 0.d0
       A(10,8) = (alpha*(drldp*p - rl)*(drvdhv + drvdp*rv)*(ulx - uvx))
     $      /denom
       A(10,9) = 0.d0
       A(10,10) = ulx
 
       B(1,1) = (alpha*(drldhl*p + rl**2)*(drvdhv + drvdp*rv)*uly - (-1
     $      .d0 + alpha)*(drldhl + drldp*rl)*(drvdhv*p + rv**2)*uvy)
     $      /denom
       B(1,2) = 0.d0
       B(1,3) = -(((-1.d0 + alpha)*alpha*(drldhl + drldp*rl)*rv**2)
     $      /denom)
       B(1,4) = 0.d0
       B(1,5) = 0.d0
       B(1,6) = ((-1.d0 + alpha)*alpha*rl**2*(drvdhv + drvdp*rv))/denom
       B(1,7) = 0.d0
       B(1,8) = ((-1.d0 + alpha)*alpha*(drldhl + drldp*rl)*(drvdhv +
     $      drvdp*rv)*(uly - uvy))/denom
       B(1,9) = 0.d0
       B(1,10) = 0.d0
 
       B(2,1) = 0.d0
       B(2,2) = uvy
       B(2,3) = 0.d0
       B(2,4) = 0.d0
       B(2,5) = 0.d0
       B(2,6) = 0.d0
       B(2,7) = 0.d0
       B(2,8) = 0.d0
       B(2,9) = 0.d0
       B(2,10) = 0.d0
 
       B(3,1) = p/(alpha*rv)
       B(3,2) = 0.d0
       B(3,3) = uvy
       B(3,4) = 0.d0
       B(3,5) = 0.d0
       B(3,6) = 0.d0
       B(3,7) = 0.d0
       B(3,8) = 1.d0/rv
       B(3,9) = 0.d0
       B(3,10) = 0.d0
 
       B(4,1) = 0.d0
       B(4,2) = 0.d0
       B(4,3) = 0.d0
       B(4,4) = uvy
       B(4,5) = 0.d0
       B(4,6) = 0.d0
       B(4,7) = 0.d0
       B(4,8) = 0.d0
       B(4,9) = 0.d0
       B(4,10) = 0.d0
 
       B(5,1) = 0.d0
       B(5,2) = 0.d0
       B(5,3) = 0.d0
       B(5,4) = 0.d0
       B(5,5) = uly
       B(5,6) = 0.d0
       B(5,7) = 0.d0
       B(5,8) = 0.d0
       B(5,9) = 0.d0
       B(5,10) = 0.d0
 
       B(6,1) = p/(-rl + alpha*rl)
       B(6,2) = 0.d0
       B(6,3) = 0.d0
       B(6,4) = 0.d0
       B(6,5) = 0.d0
       B(6,6) = uly
       B(6,7) = 0.d0
       B(6,8) = 1.d0/rl
       B(6,9) = 0.d0
       B(6,10) = 0.d0
 
       B(7,1) = 0.d0
       B(7,2) = 0.d0
       B(7,3) = 0.d0
       B(7,4) = 0.d0
       B(7,5) = 0.d0
       B(7,6) = 0.d0
       B(7,7) = uly
       B(7,8) = 0.d0
       B(7,9) = 0.d0
       B(7,10) = 0.d0
 
       B(8,1) = -(((drldhl*p + rl**2)*(drvdhv*p + rv**2)*(uly - uvy))
     $      /denom)
       B(8,2) = 0.d0
       B(8,3) = (alpha*(drldhl*p + rl**2)*rv**2)/denom
       B(8,4) = 0.d0
       B(8,5) = 0.d0
       B(8,6) = -(((-1.d0 + alpha)*rl**2*(drvdhv*p + rv**2))/denom)
       B(8,7) = 0.d0
       B(8,8) = (-((-1.d0 + alpha)*(drldhl + drldp*rl)*(drvdhv*p + rv**2
     $      )*uly) + alpha*(drldhl*p + rl**2)*(drvdhv + drvdp*rv)*uvy)
     $      /denom
       B(8,9) = 0.d0
       B(8,10) = 0.d0
 
       B(9,1) = ((drldhl*p + rl**2)*(drvdp*p - rv)*(uly - uvy))/denom
       B(9,2) = 0.d0
       B(9,3) = ((drldhl*p + rl*(drldp*(p - alpha*p) + alpha*rl))*rv)
     $      /denom
       B(9,4) = 0.d0
       B(9,5) = 0.d0
       B(9,6) = ((-1.d0 + alpha)*rl**2*(drvdp*p - rv))/denom
       B(9,7) = 0.d0
       B(9,8) = ((-1.d0 + alpha)*(drldhl + drldp*rl)*(drvdp*p - rv)*(uly
     $      - uvy))/denom
       B(9,9) = uvy
       B(9,10) = 0.d0
 
       B(10,1) = ((drldp*p - rl)*(drvdhv*p + rv**2)*(uly - uvy))/denom
       B(10,2) = 0.d0
       B(10,3) = (alpha*(-(drldp*p) + rl)*rv**2)/denom
       B(10,4) = 0.d0
       B(10,5) = 0.d0
       B(10,6) = (rl*(drvdhv*p + rv*(alpha*drvdp*p + rv - alpha*rv)))
     $      /denom
       B(10,7) = 0.d0
       B(10,8) = (alpha*(drldp*p - rl)*(drvdhv + drvdp*rv)*(uly - uvy))
     $      /denom
       B(10,9) = 0.d0
       B(10,10) = uly
 
       C(1,1) = (alpha*(drldhl*p + rl**2)*(drvdhv + drvdp*rv)*ulz - (-1
     $      .d0 + alpha)*(drldhl + drldp*rl)*(drvdhv*p + rv**2)*uvz)
     $      /denom
       C(1,2) = 0.d0
       C(1,3) = 0.d0
       C(1,4) = -(((-1.d0 + alpha)*alpha*(drldhl + drldp*rl)*rv**2)
     $      /denom)
       C(1,5) = 0.d0
       C(1,6) = 0.d0
       C(1,7) = ((-1.d0 + alpha)*alpha*rl**2*(drvdhv + drvdp*rv))/denom
       C(1,8) = ((-1.d0 + alpha)*alpha*(drldhl + drldp*rl)*(drvdhv +
     $      drvdp*rv)*(ulz - uvz))/denom
       C(1,9) = 0.d0
       C(1,10) = 0.d0
 
       C(2,1) = 0.d0
       C(2,2) = uvz
       C(2,3) = 0.d0
       C(2,4) = 0.d0
       C(2,5) = 0.d0
       C(2,6) = 0.d0
       C(2,7) = 0.d0
       C(2,8) = 0.d0
       C(2,9) = 0.d0
       C(2,10) = 0.d0
 
       C(3,1) = 0.d0
       C(3,2) = 0.d0
       C(3,3) = uvz
       C(3,4) = 0.d0
       C(3,5) = 0.d0
       C(3,6) = 0.d0
       C(3,7) = 0.d0
       C(3,8) = 0.d0
       C(3,9) = 0.d0
       C(3,10) = 0.d0
 
       C(4,1) = p/(alpha*rv)
       C(4,2) = 0.d0
       C(4,3) = 0.d0
       C(4,4) = uvz
       C(4,5) = 0.d0
       C(4,6) = 0.d0
       C(4,7) = 0.d0
       C(4,8) = 1.d0/rv
       C(4,9) = 0.d0
       C(4,10) = 0.d0
 
       C(5,1) = 0.d0
       C(5,2) = 0.d0
       C(5,3) = 0.d0
       C(5,4) = 0.d0
       C(5,5) = ulz
       C(5,6) = 0.d0
       C(5,7) = 0.d0
       C(5,8) = 0.d0
       C(5,9) = 0.d0
       C(5,10) = 0.d0
 
       C(6,1) = 0.d0
       C(6,2) = 0.d0
       C(6,3) = 0.d0
       C(6,4) = 0.d0
       C(6,5) = 0.d0
       C(6,6) = ulz
       C(6,7) = 0.d0
       C(6,8) = 0.d0
       C(6,9) = 0.d0
       C(6,10) = 0.d0
 
       C(7,1) = p/(-rl + alpha*rl)
       C(7,2) = 0.d0
       C(7,3) = 0.d0
       C(7,4) = 0.d0
       C(7,5) = 0.d0
       C(7,6) = 0.d0
       C(7,7) = ulz
       C(7,8) = 1.d0/rl
       C(7,9) = 0.d0
       C(7,10) = 0.d0
 
       C(8,1) = -(((drldhl*p + rl**2)*(drvdhv*p + rv**2)*(ulz - uvz))
     $      /denom)
       C(8,2) = 0.d0
       C(8,3) = 0.d0
       C(8,4) = (alpha*(drldhl*p + rl**2)*rv**2)/denom
       C(8,5) = 0.d0
       C(8,6) = 0.d0
       C(8,7) = -(((-1.d0 + alpha)*rl**2*(drvdhv*p + rv**2))/denom)
       C(8,8) = (-((-1.d0 + alpha)*(drldhl + drldp*rl)*(drvdhv*p + rv**2
     $      )*ulz) + alpha*(drldhl*p + rl**2)*(drvdhv + drvdp*rv)*uvz)
     $      /denom
       C(8,9) = 0.d0
       C(8,10) = 0.d0
 
       C(9,1) = ((drldhl*p + rl**2)*(drvdp*p - rv)*(ulz - uvz))/denom
       C(9,2) = 0.d0
       C(9,3) = 0.d0
       C(9,4) = ((drldhl*p + rl*(drldp*(p - alpha*p) + alpha*rl))*rv)
     $      /denom
       C(9,5) = 0.d0
       C(9,6) = 0.d0
       C(9,7) = ((-1.d0 + alpha)*rl**2*(drvdp*p - rv))/denom
       C(9,8) = ((-1.d0 + alpha)*(drldhl + drldp*rl)*(drvdp*p - rv)*(ulz
     $      - uvz))/denom
       C(9,9) = uvz
       C(9,10) = 0.d0
 
       C(10,1) = ((drldp*p - rl)*(drvdhv*p + rv**2)*(ulz - uvz))/denom
       C(10,2) = 0.d0
       C(10,3) = 0.d0
       C(10,4) = (alpha*(-(drldp*p) + rl)*rv**2)/denom
       C(10,5) = 0.d0
       C(10,6) = 0.d0
       C(10,7) = (rl*(drvdhv*p + rv*(alpha*drvdp*p + rv - alpha*rv)))
     $      /denom
       C(10,8) = (alpha*(drldp*p - rl)*(drvdhv + drvdp*rv)*(ulz - uvz))
     $      /denom
       C(10,9) = 0.d0
       C(10,10) = ulz
C
C       dVpdV
C
        dVpdV(1,1) = 1.d0
        dVpdV(1,2) = 0.d0
        dVpdV(1,3) = 0.d0
        dVpdV(1,4) = 0.d0
        dVpdV(1,5) = 0.d0
        dVpdV(1,6) = 0.d0
        dVpdV(1,7) = 0.d0
        dVpdV(1,8) = 0.d0
        dVpdV(1,9) = 0.d0
        dVpdV(1,10) = 0.d0
 
        dVpdV(2,1) = 0.d0
        dVpdV(2,2) = 1.d0
        dVpdV(2,3) = 0.d0
        dVpdV(2,4) = 0.d0
        dVpdV(2,5) = 0.d0
        dVpdV(2,6) = 0.d0
        dVpdV(2,7) = 0.d0
        dVpdV(2,8) = 0.d0
        dVpdV(2,9) = 0.d0
        dVpdV(2,10) = 0.d0
 
        dVpdV(3,1) = 0.d0
        dVpdV(3,2) = 0.d0
        dVpdV(3,3) = 1.d0
        dVpdV(3,4) = 0.d0
        dVpdV(3,5) = 0.d0
        dVpdV(3,6) = 0.d0
        dVpdV(3,7) = 0.d0
        dVpdV(3,8) = 0.d0
        dVpdV(3,9) = 0.d0
        dVpdV(3,10) = 0.d0
 
        dVpdV(4,1) = 0.d0
        dVpdV(4,2) = 0.d0
        dVpdV(4,3) = 0.d0
        dVpdV(4,4) = 1.d0
        dVpdV(4,5) = 0.d0
        dVpdV(4,6) = 0.d0
        dVpdV(4,7) = 0.d0
        dVpdV(4,8) = 0.d0
        dVpdV(4,9) = 0.d0
        dVpdV(4,10) = 0.d0
 
        dVpdV(5,1) = 0.d0
        dVpdV(5,2) = 0.d0
        dVpdV(5,3) = 0.d0
        dVpdV(5,4) = 0.d0
        dVpdV(5,5) = 1.d0
        dVpdV(5,6) = 0.d0
        dVpdV(5,7) = 0.d0
        dVpdV(5,8) = 0.d0
        dVpdV(5,9) = 0.d0
        dVpdV(5,10) = 0.d0
 
        dVpdV(6,1) = 0.d0
        dVpdV(6,2) = 0.d0
        dVpdV(6,3) = 0.d0
        dVpdV(6,4) = 0.d0
        dVpdV(6,5) = 0.d0
        dVpdV(6,6) = 1.d0
        dVpdV(6,7) = 0.d0
        dVpdV(6,8) = 0.d0
        dVpdV(6,9) = 0.d0
        dVpdV(6,10) = 0.d0
 
        dVpdV(7,1) = 0.d0
        dVpdV(7,2) = 0.d0
        dVpdV(7,3) = 0.d0
        dVpdV(7,4) = 0.d0
        dVpdV(7,5) = 0.d0
        dVpdV(7,6) = 0.d0
        dVpdV(7,7) = 1.d0
        dVpdV(7,8) = 0.d0
        dVpdV(7,9) = 0.d0
        dVpdV(7,10) = 0.d0
 
        dVpdV(8,1) = 0.d0
        dVpdV(8,2) = 0.d0
        dVpdV(8,3) = 0.d0
        dVpdV(8,4) = 0.d0
        dVpdV(8,5) = 0.d0
        dVpdV(8,6) = 0.d0
        dVpdV(8,7) = 0.d0
        dVpdV(8,8) = 1.d0
        dVpdV(8,9) = 0.d0
        dVpdV(8,10) = 0.d0
 
        dVpdV(9,1) = 0.d0
        dVpdV(9,2) = 0.d0
        dVpdV(9,3) = 0.d0
        dVpdV(9,4) = 0.d0
        dVpdV(9,5) = 0.d0
        dVpdV(9,6) = 0.d0
        dVpdV(9,7) = 0.d0
        dVpdV(9,8) = dhvdp
        dVpdV(9,9) = dhvdTv
        dVpdV(9,10) = 0.d0
 
        dVpdV(10,1) = 0.d0
        dVpdV(10,2) = 0.d0
        dVpdV(10,3) = 0.d0
        dVpdV(10,4) = 0.d0
        dVpdV(10,5) = 0.d0
        dVpdV(10,6) = 0.d0
        dVpdV(10,7) = 0.d0
        dVpdV(10,8) = dhldp
        dVpdV(10,9) = 0.d0
        dVpdV(10,10) = dhldTl
C
C       dVdVp
C
        dVdVp(1,1) = 1.d0
        dVdVp(1,2) = 0.d0
        dVdVp(1,3) = 0.d0
        dVdVp(1,4) = 0.d0
        dVdVp(1,5) = 0.d0
        dVdVp(1,6) = 0.d0
        dVdVp(1,7) = 0.d0
        dVdVp(1,8) = 0.d0
        dVdVp(1,9) = 0.d0
        dVdVp(1,10) = 0.d0
 
        dVdVp(2,1) = 0.d0
        dVdVp(2,2) = 1.d0
        dVdVp(2,3) = 0.d0
        dVdVp(2,4) = 0.d0
        dVdVp(2,5) = 0.d0
        dVdVp(2,6) = 0.d0
        dVdVp(2,7) = 0.d0
        dVdVp(2,8) = 0.d0
        dVdVp(2,9) = 0.d0
        dVdVp(2,10) = 0.d0
 
        dVdVp(3,1) = 0.d0
        dVdVp(3,2) = 0.d0
        dVdVp(3,3) = 1.d0
        dVdVp(3,4) = 0.d0
        dVdVp(3,5) = 0.d0
        dVdVp(3,6) = 0.d0
        dVdVp(3,7) = 0.d0
        dVdVp(3,8) = 0.d0
        dVdVp(3,9) = 0.d0
        dVdVp(3,10) = 0.d0
 
        dVdVp(4,1) = 0.d0
        dVdVp(4,2) = 0.d0
        dVdVp(4,3) = 0.d0
        dVdVp(4,4) = 1.d0
        dVdVp(4,5) = 0.d0
        dVdVp(4,6) = 0.d0
        dVdVp(4,7) = 0.d0
        dVdVp(4,8) = 0.d0
        dVdVp(4,9) = 0.d0
        dVdVp(4,10) = 0.d0
 
        dVdVp(5,1) = 0.d0
        dVdVp(5,2) = 0.d0
        dVdVp(5,3) = 0.d0
        dVdVp(5,4) = 0.d0
        dVdVp(5,5) = 1.d0
        dVdVp(5,6) = 0.d0
        dVdVp(5,7) = 0.d0
        dVdVp(5,8) = 0.d0
        dVdVp(5,9) = 0.d0
        dVdVp(5,10) = 0.d0
 
        dVdVp(6,1) = 0.d0
        dVdVp(6,2) = 0.d0
        dVdVp(6,3) = 0.d0
        dVdVp(6,4) = 0.d0
        dVdVp(6,5) = 0.d0
        dVdVp(6,6) = 1.d0
        dVdVp(6,7) = 0.d0
        dVdVp(6,8) = 0.d0
        dVdVp(6,9) = 0.d0
        dVdVp(6,10) = 0.d0
 
        dVdVp(7,1) = 0.d0
        dVdVp(7,2) = 0.d0
        dVdVp(7,3) = 0.d0
        dVdVp(7,4) = 0.d0
        dVdVp(7,5) = 0.d0
        dVdVp(7,6) = 0.d0
        dVdVp(7,7) = 1.d0
        dVdVp(7,8) = 0.d0
        dVdVp(7,9) = 0.d0
        dVdVp(7,10) = 0.d0
 
        dVdVp(8,1) = 0.d0
        dVdVp(8,2) = 0.d0
        dVdVp(8,3) = 0.d0
        dVdVp(8,4) = 0.d0
        dVdVp(8,5) = 0.d0
        dVdVp(8,6) = 0.d0
        dVdVp(8,7) = 0.d0
        dVdVp(8,8) = 1.d0
        dVdVp(8,9) = 0.d0
        dVdVp(8,10) = 0.d0
 
        dVdVp(9,1) = 0.d0
        dVdVp(9,2) = 0.d0
        dVdVp(9,3) = 0.d0
        dVdVp(9,4) = 0.d0
        dVdVp(9,5) = 0.d0
        dVdVp(9,6) = 0.d0
        dVdVp(9,7) = 0.d0
        dVdVp(9,8) = -dhvdp/dhvdTv
        dVdVp(9,9) = 1.d0/dhvdTv
        dVdVp(9,10) = 0.d0
 
        dVdVp(10,1) = 0.d0
        dVdVp(10,2) = 0.d0
        dVdVp(10,3) = 0.d0
        dVdVp(10,4) = 0.d0
        dVdVp(10,5) = 0.d0
        dVdVp(10,6) = 0.d0
        dVdVp(10,7) = 0.d0
        dVdVp(10,8) = dhldp*dhvdp/(dhldTl*dhvdTv)
        dVdVp(10,9) = -dhldp/(dhldTl*dhvdTv)
        dVdVp(10,10) = 1.d0/dhldTl
C
C     dVdVp*A*dVpdV
C
        do l = 1, 10
           do m = 1, 10
              P1(l,m) = 0.d0
           end do
        end do
        do l = 1, 10
           do m = 1, 10
                do n = 1, 10
                   P1(l,m) = P1(l,m) + dVdVp(l,n)*A(n,m)
                end do
           end do
        end do
        do l = 1, 10
           do m = 1, 10
              A(l,m) = 0.d0
           end do
        end do
        do l = 1, 8
            do m = 1, 10
                do n = 1, 10
                   A(l,m) = A(l,m) + P1(l,n)*dVpdV(n,m)
                end do
            end do
        end do
C
C     dVdVp*B*dVpdV
C
        do l = 1, 10
           do m = 1, 10
              P1(l,m) = 0.d0
           end do
        end do
        do l = 1, 10
                do m = 1, 10
                        do n = 1, 10
                                P1(l,m) = P1(l,m) + dVdVp(l,n)*B(n,m)
                end do
            end do
        end do
        do l = 1, 10
           do m = 1, 10
              B(l,m) = 0.d0
           end do
        end do
        do l = 1, 10
                do m = 1, 10
                        do n = 1, 10
                                B(l,m) = B(l,m) + P1(l,n)*dVpdV(n,m)
                end do
            end do
        end do
C
C     dVdVp*C*dVpdV
C
        do l = 1, 10
           do m = 1, 10
              P1(l,m) = 0.d0
           end do
        end do
        do l = 1, 10
           do m = 1, 10
                do n = 1, 10
                   P1(l,m) = P1(l,m) + dVdVp(l,n)*C(n,m)
                end do
           end do
        end do
        do l = 1, 10
           do m = 1, 10
              C(l,m) = 0.d0
           end do
        end do
        do l = 1, 10
           do m = 1, 10
                do n = 1, 10
                   C(l,m) = C(l,m) + P1(l,n)*dVpdV(n,m)
                end do
           end do
        end do
C
            dVdx(1) = MPGRAL.VPOCHA(NLCE,1)*LIMITA(1)
C           duvxdx
            dVdx(2) = MPGRUV.VPOCHA(NLCE,1)*LIMITA(2)
C           duvydx
            dVdx(3) = MPGRUV.VPOCHA(NLCE,4)*LIMITA(3)
C           duvzdx
            dVdx(4) = MPGRUV.VPOCHA(NLCE,7)*LIMITA(4)
C the same with the liquid velocity
            dVdx(5) = MPGRUL.VPOCHA(NLCE,1)*LIMITA(5)
            dVdx(6) = MPGRUL.VPOCHA(NLCE,4)*LIMITA(6)
            dVdx(7) = MPGRUL.VPOCHA(NLCE,7)*LIMITA(7)
            dVdx(8) = MPGRP.VPOCHA(NLCE,1)*LIMITA(8)
            dVdx(9) = MPGRTV.VPOCHA(NLCE,1)*LIMITA(9)
            dVdx(10) = MPGRTL.VPOCHA(NLCE,1)*LIMITA(10)
 
            dVdy(1) = MPGRAL.VPOCHA(NLCE,2)*LIMITA(1)
            dVdy(2) = MPGRUV.VPOCHA(NLCE,2)*LIMITA(2)
            dVdy(3) = MPGRUV.VPOCHA(NLCE,5)*LIMITA(3)
            dVdy(4) = MPGRUV.VPOCHA(NLCE,8)*LIMITA(4)
            dVdy(5) = MPGRUL.VPOCHA(NLCE,2)*LIMITA(5)
            dVdy(6) = MPGRUL.VPOCHA(NLCE,5)*LIMITA(6)
            dVdy(7) = MPGRUL.VPOCHA(NLCE,8)*LIMITA(7)
            dVdy(8) = MPGRP.VPOCHA(NLCE,2)*LIMITA(8)
            dVdy(9) = MPGRTV.VPOCHA(NLCE,2)*LIMITA(9)
            dVdy(10) = MPGRTL.VPOCHA(NLCE,2)*LIMITA(10)
 
            dVdz(1) = MPGRAL.VPOCHA(NLCE,3)*LIMITA(1)
            dVdz(2) = MPGRUV.VPOCHA(NLCE,3)*LIMITA(2)
            dVdz(3) = MPGRUV.VPOCHA(NLCE,6)*LIMITA(3)
            dVdz(4) = MPGRUV.VPOCHA(NLCE,9)*LIMITA(4)
            dVdz(5) = MPGRUL.VPOCHA(NLCE,3)*LIMITA(5)
            dVdz(6) = MPGRUL.VPOCHA(NLCE,6)*LIMITA(6)
            dVdz(7) = MPGRUL.VPOCHA(NLCE,9)*LIMITA(7)
            dVdz(8) = MPGRP.VPOCHA(NLCE,3)*LIMITA(8)
            dVdz(9) = MPGRTV.VPOCHA(NLCE,3)*LIMITA(9)
            dVdz(10) = MPGRTL.VPOCHA(NLCE,3)*LIMITA(10)
 
C No te olvides de afectar de la correspondiente cp las
C dVdx y dVdy correspondientes a las entalpias, no vayas
C a meter la pata con eso
            do m = 1, 10
               dV(m) = 0.d0
            end do
            do m = 1, 10
               do n = 1, 10
                  dV(m) = dV(m) + A(m,n)*dVdx(n) + B(m,n)*dVdy(n)
     &                          + C(m,n)*dVdz(n)
               end do
            end do
C
C************* Prediction avec/sans gradients limités
C
            MPALPP.VPOCHA(NLCE,1)= alpha - DELTAT*dV(1)
            MPUVP.VPOCHA(NLCE,1) = uvx -  DELTAT*dV(2)
            MPUVP.VPOCHA(NLCE,2) = uvy -  DELTAT*dV(3)
            MPUVP.VPOCHA(NLCE,3) = uvz -  DELTAT*dV(4)
            MPULP.VPOCHA(NLCE,1) = ulx -  DELTAT*dV(5)
            MPULP.VPOCHA(NLCE,2) = uly -  DELTAT*dV(6)
            MPULP.VPOCHA(NLCE,3) = ulz -  DELTAT*dV(7)
            MPPP.VPOCHA(NLCE,1)  = p  -   DELTAT*dV(8)
            MPTVP.VPOCHA(NLCE,1) = Tv -   DELTAT*dV(9)
            MPTLP.VPOCHA(NLCE,1) = Tl -   DELTAT*dV(10)
C
C************* If we only want 2nd order in space
C
C            MPALPP.VPOCHA(NLCE,1)= alpha
C            MPUVP.VPOCHA(NLCE,1) = uvx
C            MPUVP.VPOCHA(NLCE,2) = uvy
C            MPUVP.VPOCHA(NLCE,3) = uvz
C            MPULP.VPOCHA(NLCE,1) = ulx
C            MPULP.VPOCHA(NLCE,2) = uly
C            MPULP.VPOCHA(NLCE,3) = ulz
C            MPPP.VPOCHA(NLCE,1)  = p
C            MPTVP.VPOCHA(NLCE,1) = Tv
C            MPTLP.VPOCHA(NLCE,1) = Tl
C
C************* Updating densities(perfect and stiffened gas eos
C
            MPRVP.VPOCHA(NLCE,1) = gamv*MPPP.VPOCHA(NLCE,1)/
     &              (cpv*MPTVP.VPOCHA(NLCE,1)*(gamv - 1.d0))
            MPRLP.VPOCHA(NLCE,1) = gaml*(MPPP.VPOCHA(NLCE,1) + pil)/
     &             (cpl* MPTLP.VPOCHA(NLCE,1)*(gaml - 1.d0))
 
         ENDDO
C
      ENDIF
C
C
C***********************************************************************
C********* CORRECTION **************************************************
C***********************************************************************
C
C**** Boucle sur le faces
C
      IDIMP1 = IDIM + 1
      DO NLCF = 1, NFAC
C
C******* NLCF  = numero local du centre de face
C        NGCF  = numero global 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 = IPT1.NUM(1,NLCF)
         NGCF  = IPT1.NUM(2,NLCF)
         NGCED = IPT1.NUM(3,NLCF)
         NLCEG = MLENT1.LECT(NGCEG)
         NLCED = MLENT1.LECT(NGCED)
C
C******* TEST: IPT2.NUM(1,NLCF) = IPT1.NUM(2,NLCF)
C
         NGCF1 = IPT2.NUM(1,NLCF)
         IF( NGCF1 .NE. NGCF) THEN
            LOGAN = .TRUE.
            MESERR(1:40) = 'PRET, subroutine pre62f.eso             '
            GOTO 9999
         ENDIF
C
C******* Cosinus directeurs des NORMALES aux faces
C
C        On impose que les normales sont direct "Gauche" -> "Centre"
C
         INDCEL = (NGCEG-1)*IDIMP1
         XG = XCOOR(INDCEL+1)
         YG = XCOOR(INDCEL+2)
         ZG = XCOOR(INDCEL+3)
         INDCEL = (NGCF-1)*IDIMP1
         XC = XCOOR(INDCEL + 1)
         YC = XCOOR(INDCEL + 2)
         ZC = XCOOR(INDCEL+3)
         INDCEL = (NGCED-1)*IDIMP1
         XD = XCOOR(INDCEL+1)
         YD = XCOOR(INDCEL+2)
         ZD = XCOOR(INDCEL+3)
         DXG = XC - XG
         DYG = YC - YG
         DZG = ZC - ZG
         DXD = XC - XD
         DYD = YC - YD
         DZD = ZC - ZD
C
C******* On calcule le sign du pruduit scalare
C        (Normales de Castem) * (vecteur "gauche" -> "centre")
C
         CNX = MPNORM.VPOCHA(NLCF,7)
         CNY = MPNORM.VPOCHA(NLCF,8)
         CNZ = MPNORM.VPOCHA(NLCF,9)
         ORIENT = CNX * DXG + CNY * DYG + CNZ * DZG
         ORIENT = SIGN(1.0D0,ORIENT)
         IF(ORIENT .NE. 1.0D0)THEN
            LOGAN = .TRUE.
            MESERR(1:30)=
     &           'PRET , subroutine pre122.eso. '
            GOTO 9999
         ENDIF
         CNX = CNX * ORIENT
         CNY = CNY * ORIENT
         CNZ = CNZ * ORIENT
C
C********** Cosinus directeurs de tangente 1
C
         CTX = MPNORM.VPOCHA(NLCF,1) * ORIENT
         CTY = MPNORM.VPOCHA(NLCF,2) * ORIENT
         CTZ = MPNORM.VPOCHA(NLCF,3) * ORIENT
C
C********** Cosinus directeurs de tangente 2
C
         CVX = MPNORM.VPOCHA(NLCF,4) * ORIENT
         CVY = MPNORM.VPOCHA(NLCF,5) * ORIENT
         CVZ = MPNORM.VPOCHA(NLCF,6) * ORIENT
C
C
C******* Les autres MELVALs
C
C
C******* N.B.: On suppose qu'on a déjà controlle RO, P > 0
C              GAMMA \in (1,3)
C              Si non il faut le faire, en utilisant LOGBOR,
C              LOGNEG, VALER, VAL1, VAL2
C
C
C
C******* NGCEG = NGCED -> Mur
C
         IF( NGCEG .EQ. NGCED)THEN
C
C********** Sur le mur on fait de reconstruction sur l'etat gauche
C
C
C********** Etat gauche
 
            VALCEL = MPALPP.VPOCHA(NLCEG,1)
            ALCEL  = MELLAL.VPOCHA(NLCEG,1)
            DCEL   = MPGRAL.VPOCHA(NLCEG,1)*DXG +
     &               MPGRAL.VPOCHA(NLCEG,2)*DYG +
     &               MPGRAL.VPOCHA(NLCEG,3)*DZG
            alphag = VALCEL + ALCEL*DCEL
 
            VALCEL = MPUVP.VPOCHA(NLCEG,1)
            ALCEL  = MELLUV.VPOCHA(NLCEG,1)
            DCEL   = MPGRUV.VPOCHA(NLCEG,1)*DXG +
     &               MPGRUV.VPOCHA(NLCEG,2)*DYG +
     &               MPGRUV.VPOCHA(NLCEG,3)*DZG
            uvxg   = VALCEL + ALCEL*DCEL
 
            VALCEL = MPUVP.VPOCHA(NLCEG,2)
            ALCEL  = MELLUV.VPOCHA(NLCEG,2)
            DCEL   = MPGRUV.VPOCHA(NLCEG,4)*DXG +
     &               MPGRUV.VPOCHA(NLCEG,5)*DYG +
     &               MPGRUV.VPOCHA(NLCEG,6)*DZG
            uvyg   = VALCEL + ALCEL*DCEL
 
            VALCEL = MPUVP.VPOCHA(NLCEG,3)
            ALCEL  = MELLUV.VPOCHA(NLCEG,3)
            DCEL   = MPGRUV.VPOCHA(NLCEG,7)*DXG +
     &               MPGRUV.VPOCHA(NLCEG,8)*DYG +
     &               MPGRUV.VPOCHA(NLCEG,9)*DZG
            uvzg   = VALCEL + ALCEL*DCEL
 
            uvng   = uvxg*CNX + uvyg*CNY + uvzg*CNZ
            uvtg   = uvxg*CTX + uvyg*CTY + uvzg*CTZ
            uvvg   = uvxg*CVX + uvyg*CVY + uvzg*CVZ
 
            VALCEL = MPULP.VPOCHA(NLCEG,1)
            ALCEL  = MELLUL.VPOCHA(NLCEG,1)
            DCEL   = MPGRUL.VPOCHA(NLCEG,1)*DXG +
     &               MPGRUL.VPOCHA(NLCEG,2)*DYG +
     &               MPGRUL.VPOCHA(NLCEG,3)*DZG
            ulxg   = VALCEL + ALCEL*DCEL
 
            VALCEL = MPULP.VPOCHA(NLCEG,2)
            ALCEL  = MELLUL.VPOCHA(NLCEG,2)
            DCEL   = MPGRUL.VPOCHA(NLCEG,4)*DXG +
     &               MPGRUL.VPOCHA(NLCEG,5)*DYG +
     &               MPGRUL.VPOCHA(NLCEG,6)*DZG
            ulyg   = VALCEL + ALCEL*DCEL
 
            VALCEL = MPULP.VPOCHA(NLCEG,3)
            ALCEL  = MELLUL.VPOCHA(NLCEG,3)
            DCEL   = MPGRUL.VPOCHA(NLCEG,7)*DXG +
     &               MPGRUL.VPOCHA(NLCEG,8)*DYG +
     &               MPGRUL.VPOCHA(NLCEG,9)*DZG
            ulzg   = VALCEL + ALCEL*DCEL
 
            ulng   = ulxg*CNX + ulyg*CNY + ulzg*CNZ
            ultg   = ulxg*CTX + ulyg*CTY + ulzg*CTZ
            ulvg   = ulxg*CVX + ulyg*CVY + ulzg*CVZ
 
            VALCEL = MPPP.VPOCHA(NLCEG, 1)
            ALCEL  = MELLP.VPOCHA(NLCEG, 1)
            DCEL   = MPGRP.VPOCHA(NLCEG, 1)*DXG +
     &               MPGRP.VPOCHA(NLCEG, 2)*DYG +
     &               MPGRP.VPOCHA(NLCEG, 3)*DZG
            pg     = VALCEL + ALCEL*DCEL
 
            VALCEL = MPTVP.VPOCHA(NLCEG, 1)
            ALCEL  = MELLTV.VPOCHA(NLCEG, 1)
            DCEL   = MPGRTV.VPOCHA(NLCEG, 1)*DXG +
     &               MPGRTV.VPOCHA(NLCEG, 2)*DYG +
     &               MPGRTV.VPOCHA(NLCEG, 3)*DZG
            Tvg    = VALCEL + ALCEL*DCEL
 
            VALCEL = MPTLP.VPOCHA(NLCEG, 1)
            ALCEL  = MELLTL.VPOCHA(NLCEG, 1)
            DCEL   = MPGRTL.VPOCHA(NLCEG, 1)*DXG +
     &               MPGRTL.VPOCHA(NLCEG, 2)*DYG +
     &               MPGRTL.VPOCHA(NLCEG, 3)*DZG
            Tlg    = VALCEL + ALCEL*DCEL
 
C********** We calculate densities again, we use for
C           the moment, perfect and stiffened gas eos
            rvg    = gamv*pg/(cpv*Tvg*(gamv - 1))
            rlg    = gaml*(pg + pil)/(cpl*Tlg*(gaml - 1))
C
C********** Si l'on fait pas de prediction, ce n'est pas necessaire de
C           controller la positivite' de la pression et de la densité; elle
C           est déjà garantie par la proprieté LED de limiteur; mais, vu
C           que le limiteur n'est pas calculé ici, mais dans un autre
C           operateur, on le fait
C
            LOGI1 = (pg .LT. 0.D0) .OR. (rvg .LT. 0.D0)
     &              .OR. (rlg .LT. 0.D0)
C
            IF(LOGI1)THEN
C
C************* Premier ordre en espace local
C
               alphag = MPALPC.VPOCHA(NLCEG,1)
               uvng   = MPUVC.VPOCHA(NLCEG,1)*CNX +
     &                  MPUVC.VPOCHA(NLCEG,2)*CNY +
     &                  MPUVC.VPOCHA(NLCEG,3)*CNZ
               uvtg   = MPUVC.VPOCHA(NLCEG,1)*CTX +
     &                  MPUVC.VPOCHA(NLCEG,2)*CTY +
     &                  MPUVC.VPOCHA(NLCEG,3)*CTZ
               uvvg   = MPUVC.VPOCHA(NLCEG,1)*CVX +
     &                  MPUVC.VPOCHA(NLCEG,2)*CVY +
     &                  MPUVC.VPOCHA(NLCEG,3)*CVZ
               ulng   = MPULC.VPOCHA(NLCEG,1)*CNX +
     &                  MPULC.VPOCHA(NLCEG,2)*CNY +
     &                  MPULC.VPOCHA(NLCEG,3)*CNZ
               ultg   = MPULC.VPOCHA(NLCEG,1)*CTX +
     &                  MPULC.VPOCHA(NLCEG,2)*CTY +
     &                  MPULC.VPOCHA(NLCEG,3)*CTZ
               ulvg   = MPULC.VPOCHA(NLCEG,1)*CVX +
     &                  MPULC.VPOCHA(NLCEG,2)*CVY +
     &                  MPULC.VPOCHA(NLCEG,3)*CVZ
               pg     = MPPC.VPOCHA(NLCEG,1)
               Tvg    = MPTVC.VPOCHA(NLCEG,1)
               Tlg    = MPTLC.VPOCHA(NLCEG,1)
               rvg    = MPRVC.VPOCHA(NLCEG,1)
               rlg    = MPRLC.VPOCHA(NLCEG,1)
 
           ENDIF
C
C********** Son etat droite
C
            alphad = alphag
            uvnd = -1.d0*uvng
            uvtd = uvtg
            uvvd = uvvg
            ulnd = -1.d0*ulng
            ultd = ultg
            ulvd = ulvg
            pd   = pg
            Tvd  = Tvg
            Tld  = Tlg
            rvd  = rvg
            rld  = rlg
C
C************* Fin cas mur
C
         ELSE
C
C************* Etat gauche
C
            VALCEL = MPALPP.VPOCHA(NLCEG,1)
            ALCEL  = MELLAL.VPOCHA(NLCEG,1)
            DCEL   = MPGRAL.VPOCHA(NLCEG,1)*DXG +
     &               MPGRAL.VPOCHA(NLCEG,2)*DYG +
     &               MPGRAL.VPOCHA(NLCEG,3)*DZG
            alphag = VALCEL + ALCEL*DCEL
 
            VALCEL = MPUVP.VPOCHA(NLCEG,1)
            ALCEL  = MELLUV.VPOCHA(NLCEG,1)
            DCEL   = MPGRUV.VPOCHA(NLCEG,1)*DXG +
     &               MPGRUV.VPOCHA(NLCEG,2)*DYG +
     &               MPGRUV.VPOCHA(NLCEG,3)*DZG
            uvxg   = VALCEL + ALCEL*DCEL
 
            VALCEL = MPUVP.VPOCHA(NLCEG,2)
            ALCEL  = MELLUV.VPOCHA(NLCEG,2)
            DCEL   = MPGRUV.VPOCHA(NLCEG,4)*DXG +
     &               MPGRUV.VPOCHA(NLCEG,5)*DYG +
     &               MPGRUV.VPOCHA(NLCEG,6)*DZG
            uvyg   = VALCEL + ALCEL*DCEL
 
            VALCEL = MPUVP.VPOCHA(NLCEG,3)
            ALCEL  = MELLUV.VPOCHA(NLCEG,3)
            DCEL   = MPGRUV.VPOCHA(NLCEG,7)*DXG +
     &               MPGRUV.VPOCHA(NLCEG,8)*DYG +
     &               MPGRUV.VPOCHA(NLCEG,9)*DZG
            uvzg   = VALCEL + ALCEL*DCEL
 
            uvng   = uvxg*CNX + uvyg*CNY + uvzg*CNZ
            uvtg   = uvxg*CTX + uvyg*CTY + uvzg*CTZ
            uvvg   = uvxg*CVX + uvyg*CVY + uvzg*CVZ
 
            VALCEL = MPULP.VPOCHA(NLCEG,1)
            ALCEL  = MELLUL.VPOCHA(NLCEG,1)
            DCEL   = MPGRUL.VPOCHA(NLCEG,1)*DXG +
     &               MPGRUL.VPOCHA(NLCEG,2)*DYG +
     &               MPGRUL.VPOCHA(NLCEG,3)*DZG
            ulxg   = VALCEL + ALCEL*DCEL
 
            VALCEL = MPULP.VPOCHA(NLCEG,2)
            ALCEL  = MELLUL.VPOCHA(NLCEG,2)
            DCEL   = MPGRUL.VPOCHA(NLCEG,4)*DXG +
     &               MPGRUL.VPOCHA(NLCEG,5)*DYG +
     &               MPGRUL.VPOCHA(NLCEG,6)*DZG
            ulyg   = VALCEL + ALCEL*DCEL
 
            VALCEL = MPULP.VPOCHA(NLCEG,3)
            ALCEL  = MELLUL.VPOCHA(NLCEG,3)
            DCEL   = MPGRUL.VPOCHA(NLCEG,7)*DXG +
     &               MPGRUL.VPOCHA(NLCEG,8)*DYG +
     &               MPGRUL.VPOCHA(NLCEG,9)*DZG
            ulzg   = VALCEL + ALCEL*DCEL
 
            ulng   = ulxg*CNX + ulyg*CNY + ulzg*CNZ
            ultg   = ulxg*CTX + ulyg*CTY + ulzg*CTZ
            ulvg   = ulxg*CVX + ulyg*CVY + ulzg*CVZ
 
            VALCEL = MPPP.VPOCHA(NLCEG, 1)
            ALCEL  = MELLP.VPOCHA(NLCEG, 1)
            DCEL   = MPGRP.VPOCHA(NLCEG, 1)*DXG +
     &               MPGRP.VPOCHA(NLCEG, 2)*DYG +
     &               MPGRP.VPOCHA(NLCEG, 3)*DZG
            pg     = VALCEL + ALCEL*DCEL
 
            VALCEL = MPTVP.VPOCHA(NLCEG, 1)
            ALCEL  = MELLTV.VPOCHA(NLCEG, 1)
            DCEL   = MPGRTV.VPOCHA(NLCEG, 1)*DXG +
     &               MPGRTV.VPOCHA(NLCEG, 2)*DYG +
     &               MPGRTV.VPOCHA(NLCEG, 3)*DZG
            Tvg    = VALCEL + ALCEL*DCEL
 
            VALCEL = MPTLP.VPOCHA(NLCEG, 1)
            ALCEL  = MELLTL.VPOCHA(NLCEG, 1)
            DCEL   = MPGRTL.VPOCHA(NLCEG, 1)*DXG +
     &               MPGRTL.VPOCHA(NLCEG, 2)*DYG +
     &               MPGRTL.VPOCHA(NLCEG, 3)*DZG
            Tlg    = VALCEL + ALCEL*DCEL
 
C********** We calculate densities again, we use for
C           the moment, perfect and stiffened gas eos
            rvg    = gamv*pg/(cpv*Tvg*(gamv - 1))
            rlg    = gaml*(pg + pil)/(cpl*Tlg*(gaml - 1))
C
C********** Positivite
C
            LOGI1 = (pg .LT. 0.D0) .OR. (rvg .LT. 0.D0)
     &           .OR. (rlg  .LT. 0.d0)
C
            IF(LOGI1)THEN
C
C************* Premier ordre en espace local
C
               alphag = MPALPC.VPOCHA(NLCEG,1)
               uvng   = MPUVC.VPOCHA(NLCEG,1)*CNX +
     &                  MPUVC.VPOCHA(NLCEG,2)*CNY +
     &                  MPUVC.VPOCHA(NLCEG,3)*CNZ
               uvtg   = MPUVC.VPOCHA(NLCEG,1)*CTX +
     &                  MPUVC.VPOCHA(NLCEG,2)*CTY +
     &                  MPUVC.VPOCHA(NLCEG,3)*CTZ
               uvvg   = MPUVC.VPOCHA(NLCEG,1)*CVX +
     &                  MPUVC.VPOCHA(NLCEG,2)*CVY +
     &                  MPUVC.VPOCHA(NLCEG,3)*CVZ
               ulng   = MPULC.VPOCHA(NLCEG,1)*CNX +
     &                  MPULC.VPOCHA(NLCEG,2)*CNY +
     &                  MPULC.VPOCHA(NLCEG,3)*CNZ
               ultg   = MPULC.VPOCHA(NLCEG,1)*CTX +
     &                  MPULC.VPOCHA(NLCEG,2)*CTY +
     &                  MPULC.VPOCHA(NLCEG,3)*CTZ
               ulvg   = MPULC.VPOCHA(NLCEG,1)*CVX +
     &                  MPULC.VPOCHA(NLCEG,2)*CVY +
     &                  MPULC.VPOCHA(NLCEG,3)*CVZ
               pg     = MPPC.VPOCHA(NLCEG,1)
               Tvg    = MPTVC.VPOCHA(NLCEG,1)
               Tlg    = MPTLC.VPOCHA(NLCEG,1)
               rvg    = MPRVC.VPOCHA(NLCEG,1)
               rlg    = MPRLC.VPOCHA(NLCEG,1)
 
            ENDIF
C
C********** Etat droite
C
            VALCEL = MPALPP.VPOCHA(NLCED, 1)
            ALCEL  = MELLAL.VPOCHA(NLCED, 1)
            DCEL   = MPGRAL.VPOCHA(NLCED, 1)*DXD +
     &               MPGRAL.VPOCHA(NLCED, 2)*DYD +
     &               MPGRAL.VPOCHA(NLCED, 3)*DZD
            alphad = VALCEL + ALCEL * DCEL
 
            VALCEL = MPUVP.VPOCHA(NLCED, 1)
            ALCEL  = MELLUV.VPOCHA(NLCED, 1)
            DCEL   = MPGRUV.VPOCHA(NLCED, 1)*DXD +
     &               MPGRUV.VPOCHA(NLCED, 2)*DYD +
     &               MPGRUV.VPOCHA(NLCED, 3)*DZD
            uvxd   = VALCEL + ALCEL * DCEL
 
            VALCEL = MPUVP.VPOCHA(NLCED, 2)
            ALCEL  = MELLUV.VPOCHA(NLCED, 2)
            DCEL   = MPGRUV.VPOCHA(NLCED, 4)*DXD +
     &               MPGRUV.VPOCHA(NLCED, 5)*DYD +
     &               MPGRUV.VPOCHA(NLCED, 6)*DZD
            uvyd   = VALCEL + ALCEL * DCEL
 
            VALCEL = MPUVP.VPOCHA(NLCED, 3)
            ALCEL  = MELLUV.VPOCHA(NLCED, 3)
            DCEL   = MPGRUV.VPOCHA(NLCED, 7)*DXD +
     &               MPGRUV.VPOCHA(NLCED, 8)*DYD +
     &               MPGRUV.VPOCHA(NLCED, 9)*DZD
            uvzd   = VALCEL + ALCEL * DCEL
 
            uvnd   = uvxd*CNX + uvyd*CNY + uvzd*CNZ
            uvtd   = uvxd*CTX + uvyd*CTY + uvzd*CTZ
            uvvd   = uvxd*CVX + uvyd*CVY + uvzd*CVZ
 
            VALCEL = MPULP.VPOCHA(NLCED, 1)
            ALCEL  = MELLUL.VPOCHA(NLCED, 1)
            DCEL   = MPGRUL.VPOCHA(NLCED, 1)*DXD +
     &               MPGRUL.VPOCHA(NLCED, 2)*DYD +
     &               MPGRUL.VPOCHA(NLCED, 3)*DZD
            ulxd   = VALCEL + ALCEL * DCEL
 
            VALCEL = MPULP.VPOCHA(NLCED,2)
            ALCEL  = MELLUL.VPOCHA(NLCED, 2)
            DCEL   = MPGRUL.VPOCHA(NLCED, 4)*DXD +
     &               MPGRUL.VPOCHA(NLCED, 5)*DYD +
     &               MPGRUL.VPOCHA(NLCED, 6)*DZD
            ulyd   = VALCEL + ALCEL * DCEL
 
            VALCEL = MPULP.VPOCHA(NLCED,3)
            ALCEL  = MELLUL.VPOCHA(NLCED, 3)
            DCEL   = MPGRUL.VPOCHA(NLCED, 7)*DXD +
     &               MPGRUL.VPOCHA(NLCED, 8)*DYD +
     &               MPGRUL.VPOCHA(NLCED, 9)*DZD
            ulzd   = VALCEL + ALCEL * DCEL
 
            ulnd   = ulxd*CNX + ulyd*CNY + ulzd*CNZ
            ultd   = ulxd*CTX + ulyd*CTY + ulzd*CTZ
            ulvd   = ulxd*CVX + ulyd*CVY + ulzd*CVZ
 
            VALCEL = MPPP.VPOCHA(NLCED, 1)
            ALCEL  = MELLP.VPOCHA(NLCED, 1)
            DCEL   = MPGRP.VPOCHA(NLCED, 1)*DXD +
     &               MPGRP.VPOCHA(NLCED, 2)*DYD +
     &               MPGRP.VPOCHA(NLCED, 3)*DZD
            pd     = VALCEL + ALCEL * DCEL
 
            VALCEL = MPTVP.VPOCHA(NLCED, 1)
            ALCEL  = MELLTV.VPOCHA(NLCED, 1)
            DCEL   = MPGRTV.VPOCHA(NLCED, 1)*DXD +
     &               MPGRTV.VPOCHA(NLCED, 2)*DYD +
     &               MPGRTV.VPOCHA(NLCED, 3)*DZD
            Tvd    = VALCEL + ALCEL * DCEL
 
            VALCEL = MPTLP.VPOCHA(NLCED, 1)
            ALCEL  = MELLTL.VPOCHA(NLCED, 1)
            DCEL   = MPGRTL.VPOCHA(NLCED,1)*DXD +
     &               MPGRTL.VPOCHA(NLCED,2)*DYD +
     &               MPGRTL.VPOCHA(NLCED,3)*DZD
            Tld    = VALCEL + ALCEL * DCEL
C
C********** We calculate densities again, we use for
C           the moment, perfect and stiffened gas eos
            rvd    = gamv*pd/(cpv*Tvd*(gamv - 1))
            rld    = gaml*(pd + pil)/(cpl*Tld*(gaml - 1))
C
C********** Positivite
C
            LOGI1 = (pd .LT. 0.0D0) .OR. (rvd .LT. 0.0D0)
     &              .OR. (rld .LT. 0.0D0)
C
            IF(LOGI1)THEN
C
C************* Premier ordre en espace local
C
               alphad = MPALPC.VPOCHA(NLCED,1)
               uvnd   = MPUVC.VPOCHA(NLCED,1)*CNX +
     &                  MPUVC.VPOCHA(NLCED,2)*CNY +
     &                  MPUVC.VPOCHA(NLCED,3)*CNZ
               uvtd   = MPUVC.VPOCHA(NLCED,1)*CTX +
     &                  MPUVC.VPOCHA(NLCED,2)*CTY +
     &                  MPUVC.VPOCHA(NLCED,3)*CTZ
               uvvd   = MPUVC.VPOCHA(NLCED,1)*CVX +
     &                  MPUVC.VPOCHA(NLCED,2)*CVY +
     &                  MPUVC.VPOCHA(NLCED,3)*CVZ
               pd     = MPPC.VPOCHA(NLCED,1)
               Tvd    = MPTVC.VPOCHA(NLCED,1)
               Tld    = MPTLC.VPOCHA(NLCED,1)
               rvd    = MPRVC.VPOCHA(NLCED,1)
               rld    = MPRLC.VPOCHA(NLCED,1)
 
            ENDIF
         ENDIF
C
C******** Les MELVALs
C
         MLALP.VELCHE(1,NLCF) = alphag
         MLALP.VELCHE(3,NLCF) = alphad
         MLP.VELCHE(1,NLCF)   = pg
         MLP.VELCHE(3,NLCF)   = pd
         MLTV.VELCHE(1,NLCF)  = Tvg
         MLTV.VELCHE(3,NLCF)  = Tvd
         MLTL.VELCHE(1,NLCF)  = Tlg
         MLTL.VELCHE(3,NLCF)  = Tld
         MLRV.VELCHE(1,NLCF)  = rvg
         MLRV.VELCHE(3,NLCF)  = rvd
         MLRL.VELCHE(1,NLCF)  = rlg
         MLRL.VELCHE(3,NLCF)  = rld
 
         MLUVN.VELCHE(1,NLCF) = uvng
         MLUVN.VELCHE(3,NLCF) = uvnd
         MLUVT.VELCHE(1,NLCF) = uvtg
         MLUVT.VELCHE(3,NLCF) = uvtd
         MLUVV.VELCHE(1,NLCF) = uvvg
         MLUVV.VELCHE(3,NLCF) = uvvd
 
         MLVNX.VELCHE(1,NLCF) = CNX
         MLVNY.VELCHE(1,NLCF) = CNY
         MLVNZ.VELCHE(1,NLCF) = CNZ
         MLVT1X.VELCHE(1,NLCF) = CTX
         MLVT1Y.VELCHE(1,NLCF) = CTY
         MLVT1Z.VELCHE(1,NLCF) = CTZ
         MLVT2X.VELCHE(1,NLCF) = CVX
         MLVT2Y.VELCHE(1,NLCF) = CVY
         MLVT2Z.VELCHE(1,NLCF) = CVZ
 
         MLULN.VELCHE(1,NLCF) = ulng
         MLULN.VELCHE(3,NLCF) = ulnd
         MLULT.VELCHE(1,NLCF) = ultg
         MLULT.VELCHE(3,NLCF) = ultd
         MLULV.VELCHE(1,NLCF) = ulvg
         MLULV.VELCHE(3,NLCF) = ulvd
 
         MLLNX.VELCHE(1,NLCF) = CNX
         MLLNY.VELCHE(1,NLCF) = CNY
         MLLNZ.VELCHE(1,NLCF) = CNZ
         MLLT1X.VELCHE(1,NLCF) = CTX
         MLLT1Y.VELCHE(1,NLCF) = CTY
         MLLT1Z.VELCHE(1,NLCF) = CTZ
         MLVT2X.VELCHE(1,NLCF) = CVX
         MLVT2Y.VELCHE(1,NLCF) = CVY
         MLVT2Z.VELCHE(1,NLCF) = CVZ
 
      ENDDO
C
C**** Desactivation des SEGMENTs
C
      SEGDES IPT1
      SEGDES IPT2
C
C**** Le MPOVALs 'Prediction' sont detruits (si existentes)
C
      IF(LOGTEM)THEN
         SEGSUP MPALPP
         SEGSUP MPUVP
         SEGSUP MPULP
         SEGSUP MPPP
         SEGSUP MPTVP
         SEGSUP MPTLP
         SEGSUP MPRVP
         SEGSUP MPRLP
      ENDIF
C
      SEGDES MPALPC
      SEGDES MPGRAL
      SEGDES MPUVC
      SEGDES MPGRUV
      SEGDES MPULC
      SEGDES MPGRUL
      SEGDES MPPC
      SEGDES MPGRP
      SEGDES MPTVC
      SEGDES MPGRTV
      SEGDES MPTLC
      SEGDES MPGRTL
      SEGDES MPRVC
      SEGDES MPRLC
      SEGDES MPNORM
      SEGDES MELLAL
      SEGDES MELLUV
      SEGDES MELLUL
      SEGDES MELLP
      SEGDES MELLTV
      SEGDES MELLTL
 
      SEGDES MLALP
      SEGDES MLP
      SEGDES MLTV
      SEGDES MLTL
      SEGDES MLRV
      SEGDES MLRL
 
      SEGDES MLUVN
      SEGDES MLUVT
      SEGDES MLUVV
      SEGDES MLULN
      SEGDES MLULT
      SEGDES MLULV
 
      SEGDES MLVNX
      SEGDES MLVNY
      SEGDES MLVNZ
      SEGDES MLVT1X
      SEGDES MLVT1Y
      SEGDES MLVT1Z
      SEGDES MLVT2X
      SEGDES MLVT2Y
      SEGDES MLVT2Z
 
      SEGDES MLLNX
      SEGDES MLLNY
      SEGDES MLLNZ
      SEGDES MLLT1X
      SEGDES MLLT1Y
      SEGDES MLLT1Z
      SEGDES MLLT2X
      SEGDES MLLT2Y
      SEGDES MLLT2Z
C
C**** Destruction du MELNTI correspondance local/global
C
      SEGSUP MLENT1
C
 9999 CONTINUE
C
      RETURN
      END