Télécharger pre52f.eso

Retour à la liste

pre52f
C PRE52F    SOURCE    CB215821  20/11/25    13:36:31     10792          
       SUBROUTINE PRE52F(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               :  PRE121
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 APPELES (Outils)  :  KRIPAD, LICHT
C
C APPELES (Calcul)  :  AUCUN
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 LIMITA 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 LIMITA 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 LIMITA 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 LIMITA 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 LIMITA 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 LIMITA 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
C     DELTAT  : REAL*8, encrement en temps pour calculer la prediction
C
C SORTIES
C
C    C SORTIES
C
C     IALPHF      : MCHAML defined en la MELEME of pointers
C                   IFACEL, it contains the void fraction
C                   (a gauche et a droite de chaque face).
C                   Only one component  ('SCAL')
C
C     IUVF        : MCHAML  "FACEL"  vapour velocity and the
C                   director cosines (n,t) in the corresponding face;
C                   in the 2D case 6 composantes:
C                   'UVN' = normal velocity
C                   'UVT' = tangent velocity,
C                   these two velocities are defined in a local
C                   reference framework defined over the MELEME
C                   of pointers IFACE
C                   'NX' = n.x
C                   'NY' = n.y
C                   'TX' = t.x
C                   'TY' = t.y
C
C     IULF        : MCHAML  "FACEL"  liquid velocity
C                   in the 2D case 2 composantes:
C                   'ULN' = normal velocity
C                   'ULT' = tangent velocity
C
C     IPF         : MCHAML "FACEL" pressure
C                   Only one component ('SCAL')
C
C     ITVF        : MCHAML "FACEL"  vapour temperature
C                   Only one component ('SCAL')
C
C     ITVL        : MCHAML "FACEL"  liquid temperature
C                   Only one component ('SCAL')
C
C     IRVF        : MCHAML "FACEL"  vapour density
C                   Only one component ('SCAL')
C
C     IRLF        : MCHAML "FACEL"  liquid density
C                   Only one component ('SCAL')
C
C     LOGAN   : anomalie detectee (changement de la convention dans
C               la table domaine)
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 11.6.98.
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, XC, YC, XD, YD, DELTAT
     &       , DXG, DYG, DXD, DYD
     &       , CNX, CNY, CTX, CTY, ORIENT
     &       , alphag, uvxg, uvyg, uvng, uvtg, ulxg, ulyg, ulng, ultg
     &       , pg, Tvg, Tlg, rvg, rlg
     &       , alphad, uvxd, uvyd, uvnd, uvtd, ulxd, ulyd, ulnd, ultd
     &       , pd, Tvd, Tld, rvd, rld
     &       , VALCEL, DCEL, ALCEL
     &       , alpha, uvx, uvy, ulx, uly, p, Tv, Tl
     &       , rv, rl, drvdp, drvdhv, drldp, drldhl
     &       , hv, hl, cv, cl
     &       , A(8,8), B(8,8), denom
     &       , dVdx(8), dVdy(8), dV(8), LIMITA(8), dVdVp(8,8)
     &       , dVpdV(8,8), P1(8,8), P2(8,8)
     &       , dhvdp, dhvdTv, dhldp, dhldTl
      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,
     &         MLULN.MELVAL, MLULT.MELVAL,
     &         MLVNX.MELVAL, MLVNY.MELVAL, MLVTX.MELVAL, MLVTY.MELVAL,
     &         MLLNX.MELVAL, MLLNY.MELVAL, MLLTX.MELVAL, MLLTY.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
C     VAPOUR PHASE
 
      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) in (x,y) '
      MCHEL1.CONCHE(2) = ' UV in (n,t)    '
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 = 2 -> 4 composantes
C
      N2 = 4
      SEGINI MCHAM1
      MCHEL1.ICHAML(1) = MCHAM1
      MCHAM1.NOMCHE(1) = 'NVX     '
      MCHAM1.NOMCHE(2) = 'NVY     '
      MCHAM1.NOMCHE(3) = 'TVX     '
      MCHAM1.NOMCHE(4) = 'TVY     '
      MCHAM1.TYPCHE(1) = 'REAL*8          '
      MCHAM1.TYPCHE(2) = 'REAL*8          '
      MCHAM1.TYPCHE(3) = 'REAL*8          '
      MCHAM1.TYPCHE(4) = 'REAL*8          '
      SEGINI MLVNX
      SEGINI MLVNY
      SEGINI MLVTX
      SEGINI MLVTY
      MCHAM1.IELVAL(1) = MLVNX
      MCHAM1.IELVAL(2) = MLVNY
      MCHAM1.IELVAL(3) = MLVTX
      MCHAM1.IELVAL(4) = MLVTY
      SEGDES MCHAM1
C
C**** Vitesse
C
      N1EL = NFAC
      N1PTEL = 3
      N2EL = 0
      N2PTEL = 0
C
C**** MCHAML a N2 composantes:
C
C     IDIM = 2 -> 2 composantes
C
      N2 = 2
      SEGINI MCHAM1
      MCHEL1.ICHAML(2) = MCHAM1
      SEGDES MCHEL1
      MCHAM1.NOMCHE(1) = 'UVN     '
      MCHAM1.NOMCHE(2) = 'UVT     '
      MCHAM1.TYPCHE(1) = 'REAL*8         '
      MCHAM1.TYPCHE(2) = 'REAL*8         '
      SEGINI MLUVN
      SEGINI MLUVT
      MCHAM1.IELVAL(1) = MLUVN
      MCHAM1.IELVAL(2) = MLUVT
      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) in (x,y) '
      MCHEL1.CONCHE(2) = ' UL in (n,t)    '
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 = 2 -> 4 composantes
C
      N2 = 4
      SEGINI MCHAM1
      MCHEL1.ICHAML(1) = MCHAM1
      MCHAM1.NOMCHE(1) = 'NLX     '
      MCHAM1.NOMCHE(2) = 'NLY     '
      MCHAM1.NOMCHE(3) = 'TLX     '
      MCHAM1.NOMCHE(4) = 'TLY     '
      MCHAM1.TYPCHE(1) = 'REAL*8          '
      MCHAM1.TYPCHE(2) = 'REAL*8          '
      MCHAM1.TYPCHE(3) = 'REAL*8          '
      MCHAM1.TYPCHE(4) = 'REAL*8          '
      SEGINI MLLNX
      SEGINI MLLNY
      SEGINI MLLTX
      SEGINI MLLTY
      MCHAM1.IELVAL(1) = MLLNX
      MCHAM1.IELVAL(2) = MLLNY
      MCHAM1.IELVAL(3) = MLLTX
      MCHAM1.IELVAL(4) = MLLTY
      SEGDES MCHAM1
C
C**** Vitesse
C
      N1EL = NFAC
      N1PTEL = 3
      N2EL = 0
      N2PTEL = 0
C
C**** MCHAML a N2 composantes:
C
C     IDIM = 2 -> 2 composantes
C
      N2 = 2
      SEGINI MCHAM1
      MCHEL1.ICHAML(2) = MCHAM1
      SEGDES MCHEL1
      MCHAM1.NOMCHE(1) = 'ULN     '
      MCHAM1.NOMCHE(2) = 'ULT     '
      MCHAM1.TYPCHE(1) = 'REAL*8         '
      MCHAM1.TYPCHE(2) = 'REAL*8         '
      SEGINI MLULN
      SEGINI MLULT
      MCHAM1.IELVAL(1) = MLULN
      MCHAM1.IELVAL(2) = MLULT
      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,
C     MLVTX, MLVTY,
C
C     MLLNX, MLLNY,
C     MLLTX, MLLTY
C
C     MLUVN, MLUVT -> vapour velocities
C
C     MLULN, MLULT -> 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     MPALPC  -> 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********* PREDICTION **************************************************
C***********************************************************************
C
Cµµµµµµµµµµµµµµµµµµµµµ AQUI ME QUEDO
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) = MELLUL.VPOCHA(NLCE,1)
               LIMITA(5) = MELLUL.VPOCHA(NLCE,2)
               LIMITA(6) = MELLP.VPOCHA(NLCE,1)
               LIMITA(7) = MELLTV.VPOCHA(NLCE,1)
               LIMITA(8) = MELLTL.VPOCHA(NLCE,1)
            ELSE
C
C************* Gradients non limités
C
               LIMITA(1) = 1.0D0
               LIMITA(2) = 1.0D0
               LIMITA(3) = 1.0D0
               LIMITA(4) = 1.0D0
               LIMITA(5) = 1.0D0
               LIMITA(6) = 1.0D0
               LIMITA(7) = 1.0D0
               LIMITA(8) = 1.0D0
            ENDIF
 
            alpha = MPALPP.VPOCHA(NLCE,1)
            uvx = MPUVP.VPOCHA(NLCE,1)
            uvy = MPUVP.VPOCHA(NLCE,2)
            ulx = MPULP.VPOCHA(NLCE,1)
            uly = MPULP.VPOCHA(NLCE,2)
            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
            cv = 1.d0/SQRT(drvdhv/rv + drvdp)
            drldp  = rl/(p + pil)
            drldhl = -rl/hl
            cl = 1.d0/SQRT(drldhl/rl + drldp)
            dhvdp = 0.d0
            dhvdTv = Cpv
            dhldp = 0.d0
            dhldTl = Cpl
 
            denom = drldhl*(drvdhv*p + rv*(alpha*drvdp*p + rv -
     &              alpha*rv)) + rl*(alpha*rl*rv*cv**(-2) -
     &              (-1 + alpha)*drldp*(drvdhv*p + rv**2))
 
            A(1,1) = (alpha*(drldhl*p + rl**2)*rv*cv**(-2)*ulx -
     &         (-1 + alpha)*rl*cl**(-2)*(drvdhv*p + rv**2)*uvx)
     &         /denom
 
            A(1,2) = -((-1.d0 + alpha)*alpha*rl*cl**(-2)*rv**2)
     $           /denom
            A(1,3) = 0.d0
            A(1,4) = ((-1.d0 + alpha)*alpha*rl**2*rv*cv**(-2))
     $           /denom
            A(1,5) = 0.d0
            A(1,6) = ((-1.d0 + alpha)*alpha*rl*cl**(-2)*rv*cv**(-2)*(ulx
     $           -uvx))/denom
            A(1,7) = 0.d0
            A(1,8) = 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) = 1.d0/rv
            A(2,7) = 0.d0
            A(2,8) = 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(4,1) = p/(-rl + alpha*rl)
            A(4,2) = 0.d0
            A(4,3) = 0.d0
            A(4,4) = ulx
            A(4,5) = 0.d0
            A(4,6) = 1.d0/rl
            A(4,7) = 0.d0
            A(4,8) = 0.d0
 
            A(5,1) = 0.d0
            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) = 0
 
            A(6,1) = -((drldhl*p + rl**2)*(drvdhv*p + rv**2)*(ulx - uvx)
     $           )/denom
            A(6,2) = (alpha*(drldhl*p + rl**2)*rv**2)/denom
            A(6,3) = 0.d0
            A(6,4) = -((-1.d0 + alpha)*rl**2*(drvdhv*p + rv**2))
     $           /denom
            A(6,5) = 0.d0
            A(6,6) = (-((-1 + alpha)*rl*cl**(-2)*(drvdhv*p + rv**2)*ulx)
     $           +alpha*(drldhl*p + rl**2)*rv*cv**(-2)*uvx)/denom
            A(6,7) = 0.d0
            A(6,8) = 0.d0
 
            A(7,1) = ((drldhl*p + rl**2)*(drvdp*p - rv)*(ulx - uvx))
     &         /denom
            A(7,2) = ((drldhl*p + rl*(drldp*(p - alpha*p) + alpha*rl))
     $           *rv)/denom
            A(7,3) = 0.d0
            A(7,4) = ((-1.d0 + alpha)*rl**2*(drvdp*p - rv))/denom
            A(7,5) = 0.d0
            A(7,6) = ((-1.d0 + alpha)*rl*cl**(-2)*(drvdp*p - rv)*(ulx -
     $           uvx))/denom
            A(7,7) = uvx
            A(7,8) = 0.d0
 
            A(8,1) = ((drldp*p - rl)*(drvdhv*p + rv**2)*(ulx - uvx))
     &         /denom
            A(8,2) = (alpha*(-(drldp*p) + rl)*rv**2)/denom
            A(8,3) = 0.d0
            A(8,4) = (rl*(drvdhv*p + rv*(alpha*drvdp*p + rv - alpha*rv))
     $           )/denom
            A(8,5) = 0.d0
            A(8,6) = (alpha*(drldp*p - rl)*rv*cv**(-2)*(ulx - uvx))
     &         /denom
            A(8,7) = 0.d0
            A(8,8) = ulx
 
            B(1,1) = (alpha*(drldhl*p + rl**2)*rv*cv**(-2)*uly -
     &         (-1.d0 + alpha)*rl*cl**(-2)*(drvdhv*p + rv**2)*uvy)
     &         /denom
            B(1,2) = 0.d0
            B(1,3) = -((-1 + alpha)*alpha*rl*cl**(-2)*rv**2)/denom
            B(1,4) = 0.d0
            B(1,5) = ((-1.d0 + alpha)*alpha*rl**2*rv*cv**(-2))
     $           /denom
            B(1,6) = ((-1.d0 + alpha)*alpha*rl*cl**(-2)*rv*cv**(-2)*(uly
     $           -uvy))/denom
            B(1,7) = 0.d0
            B(1,8) = 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(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) = 1.d0/rv
            B(3,7) = 0.d0
            B(3,8) = 0.d0
 
            B(4,1) = 0.d0
            B(4,2) = 0.d0
            B(4,3) = 0.d0
            B(4,4) = uly
            B(4,5) = 0.d0
            B(4,6) = 0.d0
            B(4,7) = 0.d0
            B(4,8) = 0.d0
 
            B(5,1) = p/(-rl + alpha*rl)
            B(5,2) = 0.d0
            B(5,3) = 0.d0
            B(5,4) = 0.d0
            B(5,5) = uly
            B(5,6) = 1.d0/rl
            B(5,7) = 0.d0
            B(5,8) = 0.d0
 
            B(6,1) = -((drldhl*p + rl**2)*(drvdhv*p + rv**2)*(uly - uvy)
     $           )/denom
            B(6,2) = 0.d0
            B(6,3) = (alpha*(drldhl*p + rl**2)*rv**2)/denom
            B(6,4) = 0.d0
            B(6,5) = -((-1.d0 + alpha)*rl**2*(drvdhv*p + rv**2))
     $           /denom
            B(6,6) = (-((-1.d0 + alpha)*rl*cl**(-2)*(drvdhv*p + rv**2)
     $           *uly) +alpha*(drldhl*p + rl**2)*rv*cv**(-2)*uvy)
     $           /denom
            B(6,7) = 0.d0
            B(6,8) = 0.d0
 
            B(7,1) = ((drldhl*p + rl**2)*(drvdp*p - rv)*(uly - uvy))
     &         /denom
            B(7,2) = 0.d0
            B(7,3) = ((drldhl*p + rl*(drldp*(p - alpha*p) + alpha*rl))
     $           *rv)/denom
            B(7,4) = 0.d0
            B(7,5) = ((-1.d0 + alpha)*rl**2*(drvdp*p - rv))/denom
            B(7,6) = ((-1.d0 + alpha)*rl*cl**(-2)*(drvdp*p - rv)*(uly -
     $           uvy))/denom
            B(7,7) = uvy
            B(7,8) = 0.d0
 
            B(8,1) = ((drldp*p - rl)*(drvdhv*p + rv**2)*(uly - uvy))
     &         /denom
            B(8,2) = 0.d0
            B(8,3) = (alpha*(-(drldp*p) + rl)*rv**2)/denom
            B(8,4) = 0.d0
            B(8,5) = (rl*(drvdhv*p + rv*(alpha*drvdp*p + rv - alpha*rv))
     $           )/denom
            B(8,6) = (alpha*(drldp*p - rl)*rv*cv**(-2)*(uly - uvy))
     &         /denom
            B(8,7) = 0.d0
            B(8,8) = uly
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(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(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(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(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(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(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) = dhvdp
        dVpdV(7,7) = dhvdTv
        dVpdV(7,8) = 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) = dhldp
        dVpdV(8,7) = 0.d0
        dVpdV(8,8) = 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(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(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(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(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(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(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) = -dhvdp/dhvdTv
        dVdVp(7,7) = 1.d0/dhvdTv
        dVdVp(7,8) = 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) = dhldp*dhvdp/(dhldTl*dhvdTv)
        dVdVp(8,7) = -dhldp/(dhldTl*dhvdTv)
        dVdVp(8,8) = 1.d0/dhldTl
C
C     dVdVp*A*dVpdV
C
        do l = 1, 8
           do m = 1, 8
              P1(l,m) = 0.d0
           end do
        end do
        do l = 1, 8
           do m = 1, 8
                do n = 1, 8
                   P1(l,m) = P1(l,m) + dVdVp(l,n)*A(n,m)
                end do
           end do
        end do
        do l = 1, 8
           do m = 1, 8
              P2(l,m) = 0.d0
           end do
        end do
        do l = 1, 8
           do m = 1, 8
                do n = 1, 8
                   P2(l,m) = P2(l,m) + P1(l,n)*dVpdV(n,m)
                end do
           end do
        end do
        do l = 1, 8
            do m = 1, 8
                A(l,m) = P2(l,m)
            end do
        end do
C
C     dVdVp*B*dVpdV
C
        do l = 1, 8
           do m = 1, 8
              P1(l,m) = 0.d0
           end do
        end do
        do l = 1, 8
           do m = 1, 8
                do n = 1, 8
                   P1(l,m) = P1(l,m) + dVdVp(l,n)*B(n,m)
                end do
           end do
        end do
        do l = 1, 8
           do m = 1, 8
              P2(l,m) = 0.d0
           end do
        end do
        do l = 1, 8
           do m = 1, 8
                do n = 1, 8
                   P2(l,m) = P2(l,m) + P1(l,n)*dVpdV(n,m)
                end do
           end do
        end do
        do l = 1, 8
            do m = 1, 8
                B(l,m) = P2(l,m)
            end do
        end do
C
C     MEJORAR LO DE LA MULTIPLICACION DE LAS MATRICES DE CAMBIO
C
            dVdx(1) = MPGRAL.VPOCHA(NLCE,1)*LIMITA(1)
            dVdx(2) = MPGRUV.VPOCHA(NLCE,1)*LIMITA(2)
            dVdx(3) = MPGRUV.VPOCHA(NLCE,3)*LIMITA(3)
            dVdx(4) = MPGRUL.VPOCHA(NLCE,1)*LIMITA(4)
            dVdx(5) = MPGRUL.VPOCHA(NLCE,3)*LIMITA(5)
            dVdx(6) = MPGRP.VPOCHA(NLCE,1)*LIMITA(6)
            dVdx(7) = MPGRTV.VPOCHA(NLCE,1)*LIMITA(7)
            dVdx(8) = MPGRTL.VPOCHA(NLCE,1)*LIMITA(8)
 
            dVdy(1) = MPGRAL.VPOCHA(NLCE,2)*LIMITA(1)
            dVdy(2) = MPGRUV.VPOCHA(NLCE,2)*LIMITA(2)
            dVdy(3) = MPGRUV.VPOCHA(NLCE,4)*LIMITA(3)
            dVdy(4) = MPGRUL.VPOCHA(NLCE,2)*LIMITA(4)
            dVdy(5) = MPGRUL.VPOCHA(NLCE,4)*LIMITA(5)
            dVdy(6) = MPGRP.VPOCHA(NLCE,2)*LIMITA(6)
            dVdy(7) = MPGRTV.VPOCHA(NLCE,2)*LIMITA(7)
            dVdy(8) = MPGRTL.VPOCHA(NLCE,2)*LIMITA(8)
 
 
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, 8
               dV(m) = 0.d0
            end do
            do m = 1, 8
               do n = 1, 8
                  dV(m) = dV(m) + A(m,n)*dVdx(n) + B(m,n)*dVdy(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)
            MPULP.VPOCHA(NLCE,1) = ulx - DELTAT*dV(4)
            MPULP.VPOCHA(NLCE,2) = uly - DELTAT*dV(5)
            MPPP.VPOCHA(NLCE,1)  = p  - DELTAT*dV(6)
            MPTVP.VPOCHA(NLCE,1) = Tv - DELTAT*dV(7)
            MPTLP.VPOCHA(NLCE,1) = Tl - DELTAT*dV(8)
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 pre52f.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)
         INDCEL = (NGCF-1)*IDIMP1
         XC = XCOOR(INDCEL + 1)
         YC = XCOOR(INDCEL + 2)
         INDCEL = (NGCED-1)*IDIMP1
         XD = XCOOR(INDCEL+1)
         YD = XCOOR(INDCEL+2)
         DXG = XC - XG
         DYG = YC - YG
         DXD = XC - XD
         DYD = YC - YD
C
C******* On calcule le sign du pruduit scalare
C        (Normales de Castem) * (vecteur "gauche" -> "centre")
C
         CNX = MPNORM.VPOCHA(NLCF,1)
         CNY = MPNORM.VPOCHA(NLCF,2)
         ORIENT = CNX * DXG + CNY * DYG
         ORIENT = SIGN(1.0D0,ORIENT)
         IF(ORIENT .NE. 1.0D0)THEN
            LOGAN = .TRUE.
            MESERR(1:30)=
     &           'PRET , subroutine pre52f.eso. '
            GOTO 9999
         ENDIF
         CNX = CNX * ORIENT
         CNY = CNY * ORIENT
C
C********** Cosinus directeurs de tangent 2D
C
         CTX = -1.0D0 * CNY
         CTY = CNX
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
C
 
            VALCEL = MPALPP.VPOCHA(NLCEG,1)
            ALCEL  = MELLAL.VPOCHA(NLCEG,1)
            DCEL   = MPGRAL.VPOCHA(NLCEG,1)*DXG +
     &               MPGRAL.VPOCHA(NLCEG,2)*DYG
            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
            uvxg   = VALCEL + ALCEL*DCEL
 
            VALCEL = MPUVP.VPOCHA(NLCEG,2)
            ALCEL  = MELLUV.VPOCHA(NLCEG,2)
            DCEL   = MPGRUV.VPOCHA(NLCEG,3)*DXG +
     &               MPGRUV.VPOCHA(NLCEG,4)*DYG
            uvyg   = VALCEL + ALCEL*DCEL
 
            uvng   = uvxg*CNX + uvyg*CNY
            uvtg   = uvxg*CTX + uvyg*CTY
 
            VALCEL = MPULP.VPOCHA(NLCEG,1)
            ALCEL  = MELLUL.VPOCHA(NLCEG,1)
            DCEL   = MPGRUL.VPOCHA(NLCEG,1)*DXG +
     &               MPGRUL.VPOCHA(NLCEG,2)*DYG
            ulxg   = VALCEL + ALCEL*DCEL
 
            VALCEL = MPULP.VPOCHA(NLCEG,2)
            ALCEL  = MELLUL.VPOCHA(NLCEG,2)
            DCEL   = MPGRUL.VPOCHA(NLCEG,3)*DXG +
     &               MPGRUL.VPOCHA(NLCEG,4)*DYG
            ulyg   = VALCEL + ALCEL*DCEL
 
            ulng   = ulxg*CNX + ulyg*CNY
            ultg   = ulxg*CTX + ulyg*CTY
 
            VALCEL = MPPP.VPOCHA(NLCEG, 1)
            ALCEL  = MELLP.VPOCHA(NLCEG, 1)
            DCEL   = MPGRP.VPOCHA(NLCEG, 1)*DXG +
     &               MPGRP.VPOCHA(NLCEG, 2)*DYG
            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
            Tvg    = VALCEL + ALCEL*DCEL
 
            VALCEL = MPTLP.VPOCHA(NLCEG, 1)
            ALCEL  = MELLTV.VPOCHA(NLCEG, 1)
            DCEL   = MPGRTL.VPOCHA(NLCEG, 1)*DXG +
     &               MPGRTL.VPOCHA(NLCEG, 2)*DYG
            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.0D0) .OR. (rvg .LT. 0.0D0)
     &              .OR. (rlg .LT. 0.0D0)
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
               uvtg   = MPULC.VPOCHA(NLCEG,1)*CTX +
     &                  MPUVC.VPOCHA(NLCEG,2)*CTY
               ulng   = MPULC.VPOCHA(NLCEG,1)*CNX +
     &                  MPULC.VPOCHA(NLCEG,2)*CNY
               ultg   = MPULC.VPOCHA(NLCEG,1)*CTX +
     &                  MPULC.VPOCHA(NLCEG,2)*CTY
               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
            ulnd = -1.d0*ulng
            ultd = ultg
            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
            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
            uvxg   = VALCEL + ALCEL*DCEL
 
            VALCEL = MPUVP.VPOCHA(NLCEG,2)
            ALCEL  = MELLUV.VPOCHA(NLCEG,2)
            DCEL   = MPGRUV.VPOCHA(NLCEG,3)*DXG +
     &               MPGRUV.VPOCHA(NLCEG,4)*DYG
            uvyg   = VALCEL + ALCEL*DCEL
 
            uvng   = uvxg*CNX + uvyg*CNY
            uvtg   = uvxg*CTX + uvyg*CTY
 
            VALCEL = MPULP.VPOCHA(NLCEG,1)
            ALCEL  = MELLUL.VPOCHA(NLCEG,1)
            DCEL   = MPGRUL.VPOCHA(NLCEG,1)*DXG +
     &               MPGRUL.VPOCHA(NLCEG,2)*DYG
            ulxg   = VALCEL + ALCEL*DCEL
 
            VALCEL = MPULP.VPOCHA(NLCEG,2)
            ALCEL  = MELLUL.VPOCHA(NLCEG,2)
            DCEL   = MPGRUL.VPOCHA(NLCEG,3)*DXG +
     &               MPGRUL.VPOCHA(NLCEG,4)*DYG
            ulyg   = VALCEL + ALCEL*DCEL
 
            ulng   = ulxg*CNX + ulyg*CNY
            ultg   = ulxg*CTX + ulyg*CTY
 
            VALCEL = MPPP.VPOCHA(NLCEG, 1)
            ALCEL  = MELLP.VPOCHA(NLCEG, 1)
            DCEL   = MPGRP.VPOCHA(NLCEG, 1)*DXG +
     &               MPGRP.VPOCHA(NLCEG, 2)*DYG
            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
            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
            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
               uvtg   = MPULC.VPOCHA(NLCEG,1)*CTX +
     &                  MPUVC.VPOCHA(NLCEG,2)*CTY
               ulng   = MPULC.VPOCHA(NLCEG,1)*CNX +
     &                  MPULC.VPOCHA(NLCEG,2)*CNY
               ultg   = MPULC.VPOCHA(NLCEG,1)*CTX +
     &                  MPULC.VPOCHA(NLCEG,2)*CTY
               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
            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
            uvxd   = VALCEL + ALCEL * DCEL
 
            VALCEL = MPUVP.VPOCHA(NLCED, 2)
            ALCEL  = MELLUV.VPOCHA(NLCED, 2)
            DCEL   = MPGRUV.VPOCHA(NLCED, 3)*DXD +
     &               MPGRUV.VPOCHA(NLCED, 4)*DYD
            uvyd   = VALCEL + ALCEL * DCEL
 
            uvnd   = uvxd*CNX + uvyd*CNY
            uvtd   = uvxd*CTX + uvyd*CTY
 
            VALCEL = MPULP.VPOCHA(NLCED, 1)
            ALCEL  = MELLUL.VPOCHA(NLCED, 1)
            DCEL   = MPGRUL.VPOCHA(NLCED, 1)*DXD +
     &               MPGRUL.VPOCHA(NLCED, 2)*DYD
            ulxd   = VALCEL + ALCEL * DCEL
 
            VALCEL = MPULP.VPOCHA(NLCED,2)
            ALCEL  = MELLUL.VPOCHA(NLCED, 2)
            DCEL   = MPGRUL.VPOCHA(NLCED, 3)*DXD +
     &               MPGRUL.VPOCHA(NLCED, 4)*DYD
            ulyd   = VALCEL + ALCEL * DCEL
 
            ulnd   = ulxd*CNX + ulyd*CNY
            ultd   = ulxd*CTX + ulyd*CTY
 
            VALCEL = MPPP.VPOCHA(NLCED, 1)
            ALCEL  = MELLP.VPOCHA(NLCED, 1)
            DCEL   = MPGRP.VPOCHA(NLCED, 1)*DXD +
     &               MPGRP.VPOCHA(NLCED, 2)*DYD
            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
            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
            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
               uvtd   = MPUVC.VPOCHA(NLCED,1)*CTX +
     &                  MPUVC.VPOCHA(NLCED,2)*CTY
               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
         MLVNX.VELCHE(1,NLCF) = CNX
         MLVNY.VELCHE(1,NLCF) = CNY
         MLVTX.VELCHE(1,NLCF) = CTX
         MLVTY.VELCHE(1,NLCF) = CTY
 
         MLULN.VELCHE(1,NLCF) = ulng
         MLULN.VELCHE(3,NLCF) = ulnd
         MLULT.VELCHE(1,NLCF) = ultg
         MLULT.VELCHE(3,NLCF) = ultd
         MLLNX.VELCHE(1,NLCF) = CNX
         MLLNY.VELCHE(1,NLCF) = CNY
         MLLTX.VELCHE(1,NLCF) = CTX
         MLLTY.VELCHE(1,NLCF) = CTY
C
      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
C
      SEGDES MLALP
      SEGDES MLP
      SEGDES MLTV
      SEGDES MLTL
      SEGDES MLRV
      SEGDES MLRL
      SEGDES MLUVN
      SEGDES MLUVT
      SEGDES MLULN
      SEGDES MLULT
      SEGDES MLVNX
      SEGDES MLVNY
      SEGDES MLVTX
      SEGDES MLVTY
      SEGDES MLLNX
      SEGDES MLLNY
      SEGDES MLLTX
      SEGDES MLLTY
C
C**** Destruction du MELNTI correspondance local/global
C
      SEGSUP MLENT1
C
 9999 CONTINUE
C
      RETURN
      END