Télécharger pre221.eso

Retour à la liste

pre221
C PRE221    SOURCE    OF166741  24/10/03    21:15:31     12022          
      SUBROUTINE PRE221(LOGTEM,
     &                  ICEN,IFACE,IFACEL,INORM,
     &                  IROC, IGRROC, IALROC,
     &                  IVITC, IGRVC, IALVC,
     &                  IPC ,IGRPC, IALPC,
     &                  MLRECV, MLRECP, MLMOT1,
     &                  IYC ,IGRYC, IALYC,
     &                  IGAMC,
     &                  DELTAT,
     &                  IROF,IVITF,IPF,IGAMF,IYF,
     &                  LOGAN,LOGNEG,LOGBOR,MESERR,VALER,VAL1,VAL2)
C************************************************************************
C
C PROJET            :  CASTEM 2000
C
C NOM               :  PRE221
C
C DESCRIPTION       :  Voir PRE22
C
C LANGAGE           :  FORTRAN 77 + ESOPE 2000 (avec extensions CISI)
C
C AUTEUR            :  A. BECCANTINI, DRN/DMT/SEMT/TTMF
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) Autres pointeurs
C
C     IROC    : CHPOINT "CENTRE" contenant la masse volumique RHO
C
C     IGRROC  : CHPOINT "CENTRE" contenant le gradient de la
C               masse volumique RHO (2 composantes)
C
C     IALROC  : CHPOINT "CENTRE" contenant le limiteur du gradient
C               de la masse volumique
C
C     IVITC   : CHPOINT "CENTRE" contenant la vitesse UX, UY ;
C
C     IGRVC   : CHPOINT "CENTRE" contenant le gradient de la
C               vitesse (4 composantes)
C
C     IALVC  :  CHPOINT "CENTRE" contenant le limiteur du gradient
C               de la vitesse (2 composantes)
C
C     IPC     : CHPOINT "CENTRE" contenat la pression P;
C
C     IGRPC   : CHPOINT "CENTRE" contenant le gradient de la
C               pression (2 composantes)
C
C     IALPC   : CHPOINT "CENTRE" contenant le limiteur du gradient
C               de la pression
C
C     MLRECV  : pointeur de la LISTREEL des cv
C
C     MLRECP  : pointeur de la LISTREEL des cp
C
C     MLMOT1  : pointeur des noms des gaz qui apparessent dans les
C               equations d'Euler
C
C     IYC     : CHPOINT "CENTRE" contenat les fractions massiques ;
C
C     IGRYC   : CHPOINT "CENTRE" contenant les gradient des fr.mass.;
C
C     IALPC   : CHPOINT "CENTRE" contenant les limiteurs des gradients
C               des fr.mass.;
C
C     IGAMC   : CHPOINT "CENTRE" contenat le "Gamma" du gaz
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     IYF     : MCHAML  defini sur le MELEME de pointeur IFACEL,
C               contenant les  fractions massiques;
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 10.7.98.
C
C************************************************************************
C
C
C ATTENTION: Cet programme marche si le MAILLAGE est convex;
C            si non il faut changer l'algoritme 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
C**** Les variables
C
      IMPLICIT INTEGER(I-N)
      INTEGER ICEN, IFACE, IFACEL, INORM,IROC, IGRROC, IALROC
     &       , IVITC, IGRVC, IALVC
     &       , IPC ,IGRPC, IALPC
     &       , IYC, IGRYC, IALYC
     &       , IGAMC
     &       , IROF, IVITF, IPF, IGAMF, IYF
     &       , IGEOM, NFAC, NCEN
     &       , N1PTEL, N1EL, N2PTEL, N2EL, N2, N1, N3, L1, NLCE
     &       , NLCF, NGCF, NGCEG, NLCEG, NGCED, NLCED, NGCF1
     &       , IDIMP1, INDCEL, I1, NESP
      REAL*8 VALER, VAL1, VAL2, XG, YG, XC, YC, XD, YD, DELTAT
     &        ,DXG, DYG, DXD, DYD
     &       , CNX, CNY, CTX, CTY, ORIENT
     &       , ROG, PG, GAMG, UXG, UYG, UNG, UTG
     &       , ROD, PD, GAMD, UXD, UYD, UND, UTD
     &       , VALCEL, DCEL, ALCEL
     &       , DROX, DROY, DUXX, DUXY, DUYX, DUYY, DPX, DPY
     &       , DRO, DUX, DUY, DP
     &       , ALPHA, CPCEL, CVCEL, DYMAS, SUMY
      CHARACTER*(40) MESERR
      CHARACTER*(8) TYPE
      LOGICAL LOGAN,LOGNEG, LOGBOR, LOGTEM, LOGI1, LOGI2
C
C
C**** Les Includes
C
-INC SMCOORD
 
-INC PPARAM
-INC CCOPTIO
-INC SMCHPOI
      POINTEUR MPROC.MPOVAL,  MPGRR.MPOVAL,
     &         MPVITC.MPOVAL, MPGRV.MPOVAL,
     &         MPPC.MPOVAL,   MPGRP.MPOVAL,
     &         MPYC.MPOVAL,   MPGRY.MPOVAL,
     &         MPGAMC.MPOVAL, MPNORM.MPOVAL,
     &         MPROP.MPOVAL,  MPPP.MPOVAL, MPVITP.MPOVAL,
     &         MPYP.MPOVAL
-INC SMCHAML
      POINTEUR MCHAMY.MCHAML
      POINTEUR MELVNX.MELVAL, MELVNY.MELVAL,
     &         MELT1X.MELVAL, MELT1Y.MELVAL
      POINTEUR MELVUN.MELVAL, MELVUT.MELVAL
      POINTEUR MELRO.MELVAL, MELP.MELVAL,
     &         MELGAM.MELVAL
      POINTEUR MELALR.MPOVAL,
     &         MELALV.MPOVAL,
     &         MELALP.MPOVAL,
     &         MELALY.MPOVAL
-INC SMLENTI
-INC SMELEME
-INC SMLREEL
      POINTEUR MLRECV.MLREEL, MLRECP.MLREEL
-INC SMLMOTS
C
C**** Segments des fractions massiques gauche et droit
C
      SEGMENT FRAMAS
         REAL*8 FRAMG(NESP), FRAMD(NESP)
      ENDSEGMENT
C
      SEGACT MLMOT1
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     AB
C     2010/09/14: LOGI1 computed can destroy the symmetry!
C     LOGI1 = .FALSE.
      LOGI1 = .FALSE.
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     densité  + grad + limiteur
C     vitesse  + grad + limiteur
C     pression + grad + limiteur
C     fract.mass + grad + limiteur
C     gamma
C     cosinus directeurs des normales aux surface
C
      CALL LICHT(IROC  , MPROC  , TYPE, IGEOM)
      CALL LICHT(IGRROC, MPGRR  , TYPE, IGEOM)
      CALL LICHT(IVITC,  MPVITC , TYPE, IGEOM)
      CALL LICHT(IGRVC,  MPGRV  , TYPE, IGEOM)
      CALL LICHT(IPC  ,  MPPC   , TYPE, IGEOM)
      CALL LICHT(IGRPC,  MPGRP  , TYPE, IGEOM)
      CALL LICHT(IYC  ,  MPYC   , TYPE, IGEOM)
      CALL LICHT(IGRYC,  MPGRY  , TYPE, IGEOM)
      CALL LICHT(IGAMC,  MPGAMC , TYPE, IGEOM)
      CALL LICHT(INORM,  MPNORM , TYPE, IGEOM)
C
C**** Les MPOVALs 'Prediction'
C
      IF(LOGTEM)THEN
         SEGINI, MPROP = MPROC
         SEGINI, MPPP = MPPC
         SEGINI, MPVITP = MPVITC
         SEGINI, MPYP = MPYC
      ELSE
         MPROP = MPROC
         MPPP = MPPC
         MPVITP = MPVITC
         MPYP = MPYC
      ENDIF
C
C**** Les Limiteurs
C
      CALL LICHT(IALROC, MELALR , TYPE, IGEOM)
      CALL LICHT(IALVC,  MELALV , TYPE, IGEOM)
      CALL LICHT(IALPC,  MELALP , TYPE, IGEOM)
      CALL LICHT(IALYC,  MELALY , TYPE, IGEOM)
C
C
C**** Les MPOVAL sont déjà activés i.e.:
C
C     SEGACT MPROC
C     SEGACT MPGRR
C     SEGACT MPIALR
C     SEGACT MPVITC
C     SEGACT MPGRV
C     SEGACT MPIALV
C     SEGACT MPPC
C     SEGACT MPGRP
C     SEGACT MPIALP
C     SEGACT MPYC
C     SEGACT MPGRY
C     SEGACT MPIALY
C     SEGACT MPGAMC
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
      IVITF = MCHEL1
      MCHEL1.TITCHE = 'U  '
      MCHEL1.IMACHE(1) = IFACE
      MCHEL1.IMACHE(2) = IFACEL
      MCHEL1.CONCHE(1) = ' (n,t) in (x,y) '
      MCHEL1.CONCHE(2) = ' U 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) = 1
      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) = 1
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) = 'NX      '
      MCHAM1.NOMCHE(2) = 'NY      '
      MCHAM1.NOMCHE(3) = 'TX      '
      MCHAM1.NOMCHE(4) = 'TY      '
      MCHAM1.TYPCHE(1) = 'REAL*8          '
      MCHAM1.TYPCHE(2) = 'REAL*8          '
      MCHAM1.TYPCHE(3) = 'REAL*8          '
      MCHAM1.TYPCHE(4) = 'REAL*8          '
      SEGINI MELVNX
      SEGINI MELVNY
      SEGINI MELT1X
      SEGINI MELT1Y
      MCHAM1.IELVAL(1) = MELVNX
      MCHAM1.IELVAL(2) = MELVNY
      MCHAM1.IELVAL(3) = MELT1X
      MCHAM1.IELVAL(4) = MELT1Y
      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) = 'UN      '
      MCHAM1.NOMCHE(2) = 'UT      '
      MCHAM1.TYPCHE(1) = 'REAL*8          '
      MCHAM1.TYPCHE(2) = 'REAL*8          '
      SEGINI MELVUN
      SEGINI MELVUT
      MCHAM1.IELVAL(1) = MELVUN
      MCHAM1.IELVAL(2) = MELVUT
      SEGDES MCHAM1
C
C**** Densite
C
      N1 = 1
      N3 = 6
      L1 = 15
      SEGINI MCHEL2
      IROF = MCHEL2
      MCHEL2.IMACHE(1) = IFACEL
      MCHEL2.TITCHE = 'RO             '
      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) = 1
      MCHEL2.IFOCHE = IFOUR
      N2 = 1
      SEGINI MCHAM1
      MCHEL2.ICHAML(1) = MCHAM1
      SEGDES MCHEL2
      MCHAM1.NOMCHE(1) = 'SCAL    '
      MCHAM1.TYPCHE(1) = 'REAL*8          '
      SEGINI MELRO
      MCHAM1.IELVAL(1) = MELRO
      SEGDES MCHAM1
C
C**** Pression
C
      MCHEL1 = IROF
      SEGINI, MCHEL2 = MCHEL1
      IPF = MCHEL2
      MCHEL2.TITCHE = 'P            '
C
C**** MCHAM1 = MCHAML de la densite
C
      SEGINI, MCHAM2 = MCHAM1
      MCHEL2.ICHAML(1) = MCHAM2
      SEGDES MCHEL2
      SEGINI MELP
      MCHAM2.IELVAL(1) = MELP
      SEGDES MCHAM2
C
C**** Gamma
C
      MCHEL1 = IROF
      SEGINI, MCHEL2 = MCHEL1
      IGAMF = MCHEL2
      MCHEL2.TITCHE = 'GAMMA        '
C
C**** MCHAM1 = MCHAML de la densite
C
      SEGINI, MCHAM2 = MCHAM1
      MCHEL2.ICHAML(1) = MCHAM2
      SEGDES MCHEL2
      SEGINI MELGAM
      MCHAM2.IELVAL(1) = MELGAM
      SEGDES MCHAM2
C
C**** Le CHAMELEM des fractions massiques
C
      NESP = MPYC.VPOCHA(/2)
      SEGINI FRAMAS
      MCHEL1 = IROF
      SEGINI, MCHEL2 = MCHEL1
      IYF = MCHEL2
      MCHEL2.TITCHE = 'Y            '
      N2 = NESP
      SEGINI MCHAMY
      MCHEL2.ICHAML(1) = MCHAMY
      SEGDES MCHEL2
      N1EL = NFAC
      N1PTEL = 3
      N2EL = 0
      N2PTEL = 0
      DO I1 = 1, NESP
         SEGINI MELVA1
         MCHAMY.IELVAL(I1) = MELVA1
         MCHAMY.NOMCHE(I1) = MLMOT1.MOTS(I1)
         MCHAMY.TYPCHE(I1) = 'REAL*8          '
      ENDDO
C
C
C***********************************************************************
C********* PREDICTION **************************************************
C***********************************************************************
C
C**** Prediction avec gradients limités
C
C
      IF(LOGTEM)THEN
C
         IPT3 = ICEN
         SEGACT IPT3
         NCEN = IPT3.NUM(/2)
         DO NLCE = 1, NCEN
            ROG = MPROP.VPOCHA(NLCE,1)
            UXG = MPVITP.VPOCHA(NLCE,1)
            UYG = MPVITP.VPOCHA(NLCE,2)
            PG = MPPP.VPOCHA(NLCE,1)
            DROX = MPGRR.VPOCHA(NLCE,1)*MELALR.VPOCHA(NLCE,1)
            DROY = MPGRR.VPOCHA(NLCE,2)*MELALR.VPOCHA(NLCE,1)
            DUXX = MPGRV.VPOCHA(NLCE,1)*MELALV.VPOCHA(NLCE,1)
            DUXY = MPGRV.VPOCHA(NLCE,2)*MELALV.VPOCHA(NLCE,1)
            DUYX = MPGRV.VPOCHA(NLCE,3)*MELALV.VPOCHA(NLCE,2)
            DUYY = MPGRV.VPOCHA(NLCE,4)* MELALV.VPOCHA(NLCE,2)
            DPX  = MPGRP.VPOCHA(NLCE,1)* MELALP.VPOCHA(NLCE,1)
            DPY  = MPGRP.VPOCHA(NLCE,2)* MELALP.VPOCHA(NLCE,1)
            GAMG = MPGAMC.VPOCHA(NLCE,1)
            DRO  = UXG * DROX + ROG * ( DUXX + DUYY )
     &           + UYG * DROY
            DUX  = UXG * DUXX + DPX / ROG + UYG * DUXY
            DUY  = UXG * DUYX + UYG * DUYY + DPY / ROG
            DP   = GAMG * PG * (DUXX + DUYY)
     &           + UXG * DPX + UYG * DPY
C
            MPROP.VPOCHA(NLCE,1) = ROG - DELTAT * DRO
            MPVITP.VPOCHA(NLCE,1) = UXG - DELTAT * DUX
            MPVITP.VPOCHA(NLCE,2) = UYG - DELTAT * DUY
            MPPP.VPOCHA(NLCE,1) = PG - DELTAT * DP
            DO I1 = 1, NESP
               INDCEL = 2 * I1 - 1
               ALPHA = MELALY.VPOCHA(NLCE,I1)
               DYMAS = UXG * MPGRY.VPOCHA(NLCE,INDCEL) * ALPHA +
     &              UYG * MPGRY.VPOCHA(NLCE,INDCEL+1) * ALPHA
               MPYP.VPOCHA(NLCE,I1) = MPYC.VPOCHA(NLCE,I1) -
     &              DELTAT * DYMAS
            ENDDO
         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 pre221.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 pre221.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 = MPROP.VPOCHA(NLCEG, 1)
            ALCEL = MELALR.VPOCHA(NLCEG, 1)
            DCEL = MPGRR.VPOCHA(NLCEG, 1)*DXG +
     &             MPGRR.VPOCHA(NLCEG, 2)*DYG
            ROG = VALCEL + ALCEL * DCEL
C
            VALCEL = MPPP.VPOCHA(NLCEG, 1)
            ALCEL =  MELALP.VPOCHA(NLCEG, 1)
            DCEL = MPGRP.VPOCHA(NLCEG, 1)*DXG +
     &             MPGRP.VPOCHA(NLCEG, 2)*DYG
            PG = VALCEL + ALCEL * DCEL
C
            LOGI2 = .FALSE.
            SUMY = 0.0D0
            CPCEL = 0.0D0
            CVCEL = 0.0D0
            DO I1 = 1, NESP
               INDCEL = 2 * I1 - 1
               VALCEL = MPYP.VPOCHA(NLCEG,I1)
               ALCEL =  MELALY.VPOCHA(NLCEG, I1)
               DCEL = MPGRY.VPOCHA(NLCEG, INDCEL)*DXG +
     &              MPGRY.VPOCHA(NLCEG,INDCEL + 1 )*DYG
               ALCEL = VALCEL + ALCEL * DCEL
               SUMY = SUMY + ALCEL
               CPCEL = CPCEL + MLRECP.PROG(I1) * ALCEL
               CVCEL = CVCEL + MLRECV.PROG(I1) * ALCEL
               LOGI2 = LOGI2 .OR. (ALCEL .LT. 0.0D0)
               FRAMAS.FRAMG(I1) = ALCEL
            ENDDO
            LOGI2 = LOGI2 .OR. (SUMY .GT. 1.0D0)
            CPCEL = CPCEL + (1.0D0 - SUMY) * MLRECP.PROG(NESP+1)
            CVCEL = CVCEL + (1.0D0 - SUMY) * MLRECV.PROG(NESP+1)
            GAMG = CPCEL / CVCEL
C
            VALCEL = MPVITP.VPOCHA(NLCEG, 1)
            ALCEL = MELALV.VPOCHA(NLCEG, 1)
            DCEL = MPGRV.VPOCHA(NLCEG, 1)*DXG +
     &             MPGRV.VPOCHA(NLCEG, 2)*DYG
            UXG = VALCEL + ALCEL * DCEL
C
            VALCEL = MPVITP.VPOCHA(NLCEG, 2)
            ALCEL = MELALV.VPOCHA(NLCEG, 2)
            DCEL = MPGRV.VPOCHA(NLCEG, 3)*DXG +
     &             MPGRV.VPOCHA(NLCEG, 4)*DYG
            UYG = VALCEL + ALCEL * DCEL
C
            UNG = UXG * CNX + UYG * CNY
            UTG = UXG * CTX + UYG * CTY
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
C            LOGI1 = (PG .LT. 0.0D0) .OR. (ROG .LT. 0.0D0) .OR. LOGI2
C            2010/09/14: AB this destroys the symmetry!
C            LOGI1 = .FALSE.
C
C
            IF(LOGI1)THEN
C
C************* Premier ordre en espace local
C
               ROG = MPROC.VPOCHA(NLCEG,1)
               ROD = ROG
               PG  = MPPC.VPOCHA(NLCEG,1)
               PD = PG
               UNG = MPVITC.VPOCHA(NLCEG,1)*CNX +
     &              MPVITC.VPOCHA(NLCEG,2)*CNY
               UTG = MPVITC.VPOCHA(NLCEG,1)*CTX +
     &              MPVITC.VPOCHA(NLCEG,2)*CTY
               UND = -1.0D0 * UNG
               UTD = UTG
               DO I1 = 1, NESP
                  FRAMAS.FRAMG(I1) = MPYC.VPOCHA(NLCEG,I1)
                  FRAMAS.FRAMD(I1) = FRAMAS.FRAMG(I1)
               ENDDO
               GAMG = MPGAMC.VPOCHA(NLCEG,1)
               GAMD = GAMG
            ELSE
C
C********** Son etat droite
C
               ROD = ROG
               PD = PG
               UND = -1.0D0 * UNG
               UTD = UTG
               DO I1 = 1, NESP
                  FRAMAS.FRAMD(I1) = FRAMAS.FRAMG(I1)
               ENDDO
               GAMD = GAMG
            ENDIF
C
C************* Fin cas mur
C
         ELSE
C
C************* Etat gauche
C
            VALCEL = MPROP.VPOCHA(NLCEG, 1)
            ALCEL = MELALR.VPOCHA(NLCEG, 1)
            DCEL = MPGRR.VPOCHA(NLCEG, 1)*DXG +
     &             MPGRR.VPOCHA(NLCEG, 2)*DYG
            ROG = VALCEL + ALCEL * DCEL
C
            VALCEL = MPPP.VPOCHA(NLCEG, 1)
            ALCEL =  MELALP.VPOCHA(NLCEG, 1)
            DCEL = MPGRP.VPOCHA(NLCEG, 1)*DXG +
     &             MPGRP.VPOCHA(NLCEG, 2)*DYG
            PG = VALCEL + ALCEL * DCEL
C
            LOGI2 = .FALSE.
            SUMY = 0.0D0
            CPCEL = 0.0D0
            CVCEL = 0.0D0
            DO I1 = 1, NESP
               INDCEL = 2 * I1 - 1
               VALCEL = MPYP.VPOCHA(NLCEG,I1)
               ALCEL =  MELALY.VPOCHA(NLCEG, I1)
               DCEL = MPGRY.VPOCHA(NLCEG, INDCEL)*DXG +
     &              MPGRY.VPOCHA(NLCEG,INDCEL + 1 )*DYG
               ALCEL = VALCEL + ALCEL * DCEL
               SUMY = SUMY + ALCEL
               CPCEL = CPCEL + MLRECP.PROG(I1) * ALCEL
               CVCEL = CVCEL + MLRECV.PROG(I1) * ALCEL
               LOGI2 = LOGI2 .OR. (ALCEL .LT. 0.0D0)
               FRAMAS.FRAMG(I1) = ALCEL
            ENDDO
            LOGI2 = LOGI2 .OR. (SUMY .GT. 1.0D0)
            CPCEL = CPCEL + (1.0D0 - SUMY) * MLRECP.PROG(NESP+1)
            CVCEL = CVCEL + (1.0D0 - SUMY) * MLRECV.PROG(NESP+1)
            GAMG = CPCEL / CVCEL
C
            VALCEL = MPVITP.VPOCHA(NLCEG, 1)
            ALCEL = MELALV.VPOCHA(NLCEG, 1)
            DCEL = MPGRV.VPOCHA(NLCEG, 1)*DXG +
     &             MPGRV.VPOCHA(NLCEG, 2)*DYG
            UXG = VALCEL + ALCEL * DCEL
C
            VALCEL = MPVITP.VPOCHA(NLCEG, 2)
            ALCEL = MELALV.VPOCHA(NLCEG, 2)
            DCEL = MPGRV.VPOCHA(NLCEG, 3)*DXG +
     &             MPGRV.VPOCHA(NLCEG, 4)*DYG
            UYG = VALCEL + ALCEL * DCEL
C
            UNG = UXG * CNX + UYG * CNY
            UTG = UXG * CTX + UYG * CTY
C
C********** Positivite
C
C            LOGI1 = (PG .LT. 0.0D0) .OR. (ROG .LT. 0.0D0) .OR. LOGI2
C            2010/09/14: AB this destroys the symmetry!
C            LOGI1 = .FALSE.
C
            IF(LOGI1)THEN
C
C************* Premier ordre en espace local
C
               ROG = MPROC.VPOCHA(NLCEG,1)
               PG  = MPPC.VPOCHA(NLCEG,1)
               UNG = MPVITC.VPOCHA(NLCEG,1)*CNX +
     &              MPVITC.VPOCHA(NLCEG,2)*CNY
               UTG = MPVITC.VPOCHA(NLCEG,1)*CTX +
     &              MPVITC.VPOCHA(NLCEG,2)*CTY
               DO I1 = 1, NESP
                  FRAMAS.FRAMG(I1) = MPYC.VPOCHA(NLCEG,I1)
               ENDDO
               GAMG = MPGAMC.VPOCHA(NLCEG,1)
            ENDIF
C
C********** Etat droite
C
            VALCEL = MPROP.VPOCHA(NLCED, 1)
            ALCEL = MELALR.VPOCHA(NLCED, 1)
            DCEL = MPGRR.VPOCHA(NLCED, 1)*DXD +
     &             MPGRR.VPOCHA(NLCED, 2)*DYD
            ROD = VALCEL + ALCEL * DCEL
C
            VALCEL = MPPP.VPOCHA(NLCED, 1)
            ALCEL =  MELALP.VPOCHA(NLCED, 1)
            DCEL = MPGRP.VPOCHA(NLCED, 1)*DXD +
     &             MPGRP.VPOCHA(NLCED, 2)*DYD
            PD = VALCEL + ALCEL * DCEL
C
            LOGI2 = .FALSE.
            SUMY = 0.0D0
            CPCEL = 0.0D0
            CVCEL = 0.0D0
            DO I1 = 1, NESP
               INDCEL = 2 * I1 - 1
               VALCEL = MPYP.VPOCHA(NLCED,I1)
               ALCEL =  MELALY.VPOCHA(NLCED, I1)
               DCEL = MPGRY.VPOCHA(NLCED, INDCEL)*DXD +
     &              MPGRY.VPOCHA(NLCED,INDCEL + 1 )*DYD
               ALCEL = VALCEL + ALCEL * DCEL
               SUMY = SUMY + ALCEL
               CPCEL = CPCEL + MLRECP.PROG(I1) * ALCEL
               CVCEL = CVCEL + MLRECV.PROG(I1) * ALCEL
               LOGI2 = LOGI2 .OR. (ALCEL .LT. 0.0D0)
               FRAMAS.FRAMD(I1) = ALCEL
            ENDDO
            LOGI2 = LOGI2 .OR. (SUMY .GT. 1.0D0)
            CPCEL = CPCEL + (1.0D0 - SUMY) * MLRECP.PROG(NESP+1)
            CVCEL = CVCEL + (1.0D0 - SUMY) * MLRECV.PROG(NESP+1)
            GAMD = CPCEL / CVCEL
C
            VALCEL = MPVITP.VPOCHA(NLCED, 1)
            ALCEL = MELALV.VPOCHA(NLCED, 1)
            DCEL = MPGRV.VPOCHA(NLCED, 1)*DXD +
     &             MPGRV.VPOCHA(NLCED, 2)*DYD
            UXD = VALCEL + ALCEL * DCEL
C
            VALCEL = MPVITP.VPOCHA(NLCED, 2)
            ALCEL = MELALV.VPOCHA(NLCED, 2)
            DCEL = MPGRV.VPOCHA(NLCED, 3)*DXD +
     &             MPGRV.VPOCHA(NLCED, 4)*DYD
            UYD = VALCEL + ALCEL * DCEL
C
            UND = UXD * CNX + UYD * CNY
            UTD = UXD * CTX + UYD * CTY
C
C********** Positivite
C
C            LOGI1 = (PD .LT. 0.0D0) .OR. (ROD .LT. 0.0D0) .OR. LOGI2
C            2010/09/14: AB this destroys the symmetry!
C            LOGI1 = .FALSE.
C
            IF(LOGI1)THEN
C
C************* Premier ordre en espace local
C
               ROD = MPROC.VPOCHA(NLCED,1)
               PD  = MPPC.VPOCHA(NLCED,1)
               UND = MPVITC.VPOCHA(NLCED,1)*CNX +
     &              MPVITC.VPOCHA(NLCED,2)*CNY
               UTD = MPVITC.VPOCHA(NLCED,1)*CTX +
     &              MPVITC.VPOCHA(NLCED,2)*CTY
               DO I1 = 1, NESP
                  FRAMAS.FRAMD(I1) = MPYC.VPOCHA(NLCED,I1)
               ENDDO
               GAMD = MPGAMC.VPOCHA(NLCED,1)
            ENDIF
         ENDIF
C
C******** Les MELVALs
C
         MELRO.VELCHE(1,NLCF) = ROG
         MELRO.VELCHE(3,NLCF) = ROD
         MELP.VELCHE(1,NLCF) = PG
         MELP.VELCHE(3,NLCF) = PD
         MELGAM.VELCHE(1,NLCF) = GAMG
         MELGAM.VELCHE(3,NLCF) = GAMD
         MELVUN.VELCHE(1,NLCF) = UNG
         MELVUN.VELCHE(3,NLCF) = UND
         MELVUT.VELCHE(1,NLCF) = UTG
         MELVUT.VELCHE(3,NLCF) = UTD
         MELVNX.VELCHE(1,NLCF) = CNX
         MELVNY.VELCHE(1,NLCF) = CNY
         MELT1X.VELCHE(1,NLCF) = CTX
         MELT1Y.VELCHE(1,NLCF) = CTY
         DO I1 = 1, NESP
            MELVA1 = MCHAMY.IELVAL(I1)
            MELVA1.VELCHE(1,NLCF) = FRAMAS.FRAMG(I1)
            MELVA1.VELCHE(3,NLCF) = FRAMAS.FRAMD(I1)
         ENDDO
      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 MPROP
         SEGSUP MPVITP
         SEGSUP MPPP
         SEGSUP MPYP
      ENDIF
C
      SEGDES MPROC
      SEGDES MPGRR
      SEGDES MELALR
      SEGDES MPVITC
      SEGDES MPGRV
      SEGDES MELALV
      SEGDES MPPC
      SEGDES MPGRP
      SEGDES MELALP
      SEGDES MPYC
      SEGDES MPGRY
      SEGDES MELALY
      SEGDES MPGAMC
      SEGDES MPNORM
C
      SEGDES MELRO
      SEGDES MELP
      SEGDES MELGAM
      SEGDES MELVUN
      SEGDES MELVUT
      SEGDES MELVNX
      SEGDES MELVNY
      SEGDES MELT1X
      SEGDES MELT1Y
C
      SEGDES MPYC
      DO I1 = 1, NESP
         MELVA1 = MCHAMY.IELVAL(I1)
         SEGDES MELVA1
      ENDDO
      SEGDES MCHAMY
      SEGSUP FRAMAS
C
      SEGDES MLMOT1
CC
C**** Destruction du MELNTI correspondance local/global
C
      SEGSUP MLENT1
C
 9999 CONTINUE
C
      RETURN
      END
 
 
 
 
 
 
 
 
 
 
 
 
 
 

© Cast3M 2003 - Tous droits réservés.
Mentions légales