Sources Esope | Cast3M

konjr1
C KONJR1    SOURCE    OF166741  24/12/13    21:16:42     12097          
      SUBROUTINE KONJR1(ILINC,IRN,IUN,IPN,IGAMN,INORM,ICHPVO
     $     ,ICHPSU,IUINF,IUPRI,MELEMC,MELEFE,MELLIM,IMAT)
C
C************************************************************************
C
C PROJET            :  CASTEM 2000
C
C NOM               :  KONJR1
C
C DESCRIPTION       :  Voir KON12
C                      Calcul du jacobien du résidu pour la méthode
C                      RUSANOLM (CENTERED + diff. num)
C
C                      Cas deux dimensions, gaz "calorically perfect"
C
C LANGAGE           :  FORTRAN 77 + ESOPE 2000 (avec estensions CISI)
C
C AUTEUR            :  A. BECCANTINI, DRN/DMT/SEMT/LTMF
C
C************************************************************************
C
C
C APPELES (Outils
C          CASTEM)  :  KRIPAD, LICHT, ERREUR
C
C APPELES (Calcul)  :  CENTJ0, CENTJ2, FRBM1
C
C************************************************************************
C
C ENTREES
C
C     ILINC  : liste des inconnues (pointeur d'un objet de type LISTMOTS)
C
C  1) Pointeurs des CHPOINT
C
C     IRN     : CHPOINT CENTRE  contenant la masse volumique ;
C
C     IUN     : CHPOINT CENTRE  contenant la vitesse ;
C
C     IPN     : CHPOINT CENTRE  contenant la pression ;
C
C     IGAMN   : CHPOINT CENTRE  contenant le gamma ;
C
C     INORM   : CHPOINT FACE    contenant les normales aux faces ;
C
C     ICHPVO  : CHPOINT VOLUME  contenant le volume
C
C     ICHPSU  : CHPOINT FACE    contenant la surface des faces
C
C
C  2) Pointeurs de MELEME de la table DOMAINE
C
C     MELEMC  : MELEME 'CENTRE' du SPG des CENTRES
C
C     MELEFE  : MELEME 'FACEL' du connectivité Faces -> Elts
C
C     MELLIM  : MELEME sur lequel on ne veut pas calculer
C               la contribution à la matrice jacobienne
C
C  3) Cas Bas Mach
C
C     IUINF   : CHPOINT, one component, "cut-off velocity"
C               0 if there is no Bas MACH
C
C     IUPRI   : CHPOINT, one component, second "cut-off velocity"
C               0 if there is no Bas MACH
C
C
C SORTIES
C
C     IMAT    : pointeur de la MATRIK du jacobien du residu
C
C************************************************************************
C
C HISTORIQUE (Anomalies et modifications éventuelles)
C
C HISTORIQUE :
C
C************************************************************************
C
C
C N.B.: On suppose qu'on a déjà controllé RO, P > 0
C                                         GAMMA \in (1,3)
C       Si non il faut le faire!!!
C
C************************************************************************
C
      IMPLICIT INTEGER(I-N)
      INTEGER ILINC, IRN,IUN,IPN,IGAMN,INORM,ICHPVO,ICHPSU, IUINF, IUPRI
     &        , IMAT, IGEOMC, IGEOMF
     &        , NFAC, NBSOUS, NBREF, NBELEM, NBNN, NRIGE, NMATRI, NKID
     &        , NKMT, NBME, NBEL,  MP, NP
     &        , IFAC, NGCF, NLCF, NGCG, NGCD, NLCG, NLCD, NLFL, I1
      REAL*8 ROG, PG, UXG, UYG, GAMG, VOLG
     &       , ROD, PD, UXD, UYD, GAMD, VOLD
     &       , SURF, CNX, CNY, CTX, CTY, FUNCEL
     &       , DFRO(4), DFRET(4), DFRUN(4), DFRUT(4), VINF
     &       , UNG, UTG, UND, UTD, ROF, PF, UNF, UTF, GAMF, AMAT(4,4)
     &       , RLAMMA, BMAT(4,4)
      CHARACTER*8 TYPE
C
C**** LES INCLUDES
C
 
-INC PPARAM
-INC CCOPTIO
-INC SMCHPOI
-INC SMMATRIK
-INC SMELEME
-INC SMLMOTS
-INC SMLENTI
      POINTEUR MPRN.MPOVAL, MPUN.MPOVAL, MPPN.MPOVAL, MPGAMN.MPOVAL,
     &         MPNORM.MPOVAL, MPVOLU.MPOVAL, MPOVSU.MPOVAL,
     &         MPUPRI.MPOVAL, MPUINF.MPOVAL
      POINTEUR MELEMC.MELEME, MELEMF.MELEME, MELEFE.MELEME,
     &         MELEDU.MELEME, MELLIM.MELEME
      POINTEUR MLENTC.MLENTI, MLENTF.MLENTI, MLELIM.MLENTI
      POINTEUR RR.IZAFM, RUX.IZAFM, RUY.IZAFM, RRET.IZAFM,
     &         UXR.IZAFM, UXUX.IZAFM, UXUY.IZAFM, UXRET.IZAFM,
     &         UYR.IZAFM, UYUX.IZAFM, UYUY.IZAFM, UYRET.IZAFM,
     &         RETR.IZAFM, RETUX.IZAFM, RETUY.IZAFM, RETRET.IZAFM
      POINTEUR MLMINC.MLMOTS
C
C**** KRIPAD pour la correspondance global/local des conditions limits
C
      CALL KRIPAD(MELLIM,MLELIM)
c     SEGACT MELLIM
C
C**** KRIPAD pour la correspondance global/local des centres
C
      CALL KRIPAD(MELEMC,MLENTC)
C
C     SEGACT MLENTC
      SEGACT MELEMC
C
      SEGACT MELEFE
C
      CALL LICHT(ICHPSU,MPOVSU,TYPE,IGEOMF)
      CALL LICHT(INORM,MPNORM,TYPE,IGEOMF)
      CALL LICHT(ICHPVO,MPVOLU,TYPE,IGEOMC)
C
C**** LICHT active les MPOVALs en *MOD
C
C     i.e.
C
C     SEGACT MPOVSU*MOD
C     SEGACT MPOVNO*MOD
C     SEGACT MPVOLU*MOD
C
      MELEMF = IGEOMF
      CALL KRIPAD(MELEMF,MLENTF)
C
C     SEGACT MLENTF
      SEGACT MELEMF
C
      CALL LICHT(IRN,MPRN,TYPE,IGEOMC)
      CALL LICHT(IPN,MPPN,TYPE,IGEOMC)
      CALL LICHT(IUN,MPUN,TYPE,IGEOMC)
      CALL LICHT(IGAMN,MPGAMN,TYPE,IGEOMC)
C
C     SEGACT MPRN*MOD
C     SEGACT MPPN*MOD
C     SEGACT MPUN*MOD
C     SEGACT MPGAMN*MOD
C
C**************************************************************
C    Bas Mach
C**************************************************************
      CALL LICHT(IUPRI,MPUPRI,TYPE,IGEOMC)
      CALL LICHT(IUINF,MPUINF,TYPE,IGEOMC)
C
      NFAC = MELEFE.NUM(/2)
C
C**** Maillage des inconnues primales
C
      NBSOUS = 0
      NBREF = 0
      NBELEM =  NFAC
      NBNN = 2
C
      SEGINI MELEDU
C     MELEPR = MELEDU
C
C**** MELEDU = 'SEG2'
C
      MELEDU.ITYPEL =  2
C
      NRIGE = 7
      NMATRI = 1
      NKID = 9
      NKMT = 7
C
      SEGINI MATRIK
      IMAT = MATRIK
      MATRIK.IRIGEL(1,1) = MELEDU
      MATRIK.IRIGEL(2,1) = MELEDU
C
C**** Matrice non symetrique
C
      MATRIK.IRIGEL(7,1) = 2
C
      NBME = 16
      NBSOUS = 1
      SEGINI IMATRI
      MLMINC = ILINC
      SEGACT MLMINC
      MATRIK.IRIGEL(4,1) = IMATRI
C
      IMATRI.LISPRI(1) = MLMINC.MOTS(1)
      IMATRI.LISPRI(2) = MLMINC.MOTS(2)
      IMATRI.LISPRI(3) = MLMINC.MOTS(3)
      IMATRI.LISPRI(4) = MLMINC.MOTS(4)
      IMATRI.LISPRI(5) = MLMINC.MOTS(1)
      IMATRI.LISPRI(6) = MLMINC.MOTS(2)
      IMATRI.LISPRI(7) = MLMINC.MOTS(3)
      IMATRI.LISPRI(8) = MLMINC.MOTS(4)
      IMATRI.LISPRI(9) = MLMINC.MOTS(1)
      IMATRI.LISPRI(10) = MLMINC.MOTS(2)
      IMATRI.LISPRI(11) = MLMINC.MOTS(3)
      IMATRI.LISPRI(12) = MLMINC.MOTS(4)
      IMATRI.LISPRI(13) = MLMINC.MOTS(1)
      IMATRI.LISPRI(14) = MLMINC.MOTS(2)
      IMATRI.LISPRI(15) = MLMINC.MOTS(3)
      IMATRI.LISPRI(16) = MLMINC.MOTS(4)
C
      IMATRI.LISDUA(1) = MLMINC.MOTS(1)
      IMATRI.LISDUA(2) = MLMINC.MOTS(1)
      IMATRI.LISDUA(3) = MLMINC.MOTS(1)
      IMATRI.LISDUA(4) = MLMINC.MOTS(1)
      IMATRI.LISDUA(5) = MLMINC.MOTS(2)
      IMATRI.LISDUA(6) = MLMINC.MOTS(2)
      IMATRI.LISDUA(7) = MLMINC.MOTS(2)
      IMATRI.LISDUA(8) = MLMINC.MOTS(2)
      IMATRI.LISDUA(9) = MLMINC.MOTS(3)
      IMATRI.LISDUA(10) = MLMINC.MOTS(3)
      IMATRI.LISDUA(11) = MLMINC.MOTS(3)
      IMATRI.LISDUA(12) = MLMINC.MOTS(3)
      IMATRI.LISDUA(13) = MLMINC.MOTS(4)
      IMATRI.LISDUA(14) = MLMINC.MOTS(4)
      IMATRI.LISDUA(15) = MLMINC.MOTS(4)
      IMATRI.LISDUA(16) = MLMINC.MOTS(4)
C
      NBEL = NBELEM
      NBSOUS = 1
      NP = 2
      MP = 2
      SEGINI   RR  , RUX , RUY , RRET ,
     &         UXR , UXUX , UXUY , UXRET ,
     &         UYR , UYUX , UYUY , UYRET ,
     &         RETR , RETUX , RETUY , RETRET
C
C**** Duale = IMATRI.LISDUA(1) = 'RN'
C     Primale = IMATRI.LISPRI(1) = 'RN'
C     -> IMATRI.LIZAFM(1,1) = RR
C
C     Duale = IMATRI.LISDUA(2) = 'RN'
C     Primale = IMATRI.LISPRI(1) = 'RUXN'
C     -> IMATRI.LIZAFM(1,2) = RUX
C     ...
C
      IMATRI.LIZAFM(1,1) = RR
      IMATRI.LIZAFM(1,2) = RUX
      IMATRI.LIZAFM(1,3) = RUY
      IMATRI.LIZAFM(1,4) = RRET
      IMATRI.LIZAFM(1,5) = UXR
      IMATRI.LIZAFM(1,6) = UXUX
      IMATRI.LIZAFM(1,7) = UXUY
      IMATRI.LIZAFM(1,8) = UXRET
      IMATRI.LIZAFM(1,9) = UYR
      IMATRI.LIZAFM(1,10) = UYUX
      IMATRI.LIZAFM(1,11) = UYUY
      IMATRI.LIZAFM(1,12) = UYRET
      IMATRI.LIZAFM(1,13) = RETR
      IMATRI.LIZAFM(1,14) = RETUX
      IMATRI.LIZAFM(1,15) = RETUY
      IMATRI.LIZAFM(1,16) = RETRET
C
      DO IFAC = 1, NFAC, 1
         NGCF = MELEFE.NUM(2,IFAC)
         NLCF = MLENTF.LECT(NGCF)
         IF(NLCF .NE. IFAC)THEN
            WRITE(IOIMP,*) 'Il ne faut pas jouer avec la table domaine'
            CALL ERREUR(5)
            GOTO 9999
         ENDIF
         NGCG = MELEFE.NUM(1,IFAC)
         NGCD = MELEFE.NUM(3,IFAC)
         NLFL = MLELIM.LECT(NGCF)
         IF(NLFL .NE. 0)THEN
C
C********** The point belongs on BC -> No contribution to jacobian!
C
            MELEDU.NUM(1,IFAC) = NGCG
            MELEDU.NUM(2,IFAC) = NGCD
 
         ELSEIF(NGCG .NE. NGCD)THEN
C
C********** Les MELEMEs
C
            MELEDU.NUM(1,IFAC) = NGCG
            MELEDU.NUM(2,IFAC) = NGCD
C
C********** Les etats G et D
C
            NLCG = MLENTC.LECT(NGCG)
            NLCD = MLENTC.LECT(NGCD)
C
            ROG = MPRN.VPOCHA(NLCG,1)
            PG = MPPN.VPOCHA(NLCG,1)
            UXG = MPUN.VPOCHA(NLCG,1)
            UYG = MPUN.VPOCHA(NLCG,2)
            GAMG = MPGAMN.VPOCHA(NLCG,1)
            VOLG = MPVOLU.VPOCHA(NLCG,1)
C
            ROD = MPRN.VPOCHA(NLCD,1)
            PD = MPPN.VPOCHA(NLCD,1)
            UXD = MPUN.VPOCHA(NLCD,1)
            UYD = MPUN.VPOCHA(NLCD,2)
            GAMD = MPGAMN.VPOCHA(NLCD,1)
            VOLD = MPVOLU.VPOCHA(NLCD,1)
C
C********** La normale G->D
C           La tangente
C
            SURF = MPOVSU.VPOCHA(NLCF,1)
            CNX = MPNORM.VPOCHA(NLCF,1)
            CNY = MPNORM.VPOCHA(NLCF,2)
            CTX = -1.0D0 * CNY
            CTY = CNX
C
C********** Cut off
C
            VINF=MPUPRI.VPOCHA(NLCG,1)
            VINF=MAX(VINF,MPUPRI.VPOCHA(NLCD,1))
            VINF=MAX(VINF,MPUINF.VPOCHA(NLCG,1))
            VINF=MAX(VINF,MPUINF.VPOCHA(NLCD,1))
C
C********** L'etat moyen au centre
C
            UNG=UXG*CNX+UYG*CNY
            UTG=UXG*CTX+UYG*CTY
            UND=UXD*CNX+UYD*CNY
            UTD=UXD*CTX+UYD*CTY
            UNF=0.5D0*(UNG+UND)
            UTF=0.5D0*(UTG+UTD)
            ROF=0.5D0*(ROG+ROD)
            PF=0.5D0*(PG+PD)
            GAMF=0.5D0*(GAMG+GAMD)
C
C********** We include in the cut-off UNG,UND,UTG,UTD in order to
C           avoid low diffusivity in stagnation regions
C
            VINF=MAX(VINF,((UNG*UNG+UTG*UTG)**0.5D0))
            VINF=MAX(VINF,((UND*UND+UTD*UTD)**0.5D0))
C
C********** Calcule de la matrice de diffusivité pour le variables
C           conservatives.
C           Dans le repaire (n,t) local
C           F(UL,Ur) = 0.5 (F(UL)+F(UR)+AMAT*RLAMMA*(UL-UR))
C
            CALL FRBM1(VINF, ROF,UNF,UTF,PF,GAMF,AMAT, RLAMMA)
C
C           Ici U=(\rho,\rho un, \rho ut, \rho et)
C           Il faut calculer la contribution sur F dans le repaire
C           (n,t) par rapport aux variables
C           U=(\rho,\rho ux, \rho uy, \rho et)
C
            DO I1 = 1,4,1
               BMAT(I1,1) = AMAT(I1,1)
               BMAT(I1,2) = AMAT(I1,2) * CNX + AMAT(I1,3) * CTX
               BMAT(I1,3) = AMAT(I1,2) * CNY + AMAT(I1,3) * CTY
               BMAT(I1,4) = AMAT(I1,4)
            ENDDO
C
C********** La contribution de Gauche
C
C           Partie centrée
C
            CALL CENTJ0(ROG,UXG,UYG,PG,GAMG,CNX,CNY,CTX,CTY,
     &           DFRO,DFRUN,DFRUT,DFRET)
C
C********** Centrée + Decentrée
C
C           DFRO      = derivative of the numerical approx. of (\rho un)
C                       in the frame (x,y) (centered part)
C           BMAT(1,*) = derivative of the numerical approx. of (\rho un)
C                       in the frame (x,y) (diff. part)
C
            DO I1=1,4,1
               DFRO(I1) = DFRO(I1) + 0.5D0 * RLAMMA * BMAT(1,I1)
               DFRUN(I1) = DFRUN(I1) +  0.5D0 * RLAMMA * BMAT(2,I1)
               DFRUT(I1) = DFRUT(I1) +  0.5D0 * RLAMMA * BMAT(3,I1)
               DFRET(I1) = DFRET(I1) +  0.5D0 * RLAMMA * BMAT(4,I1)
            ENDDO
C
C
C********** AB.AM(IFAC,IPRIM,IDUAL)
C           A = nom de l'inconnu duale (Ro,rUX,rUY,RET)
C           B = nom de l'inconnu primale (Ro,rUX,rUY,RET)
C           IPRIM = 1, 2  -> G, D
C           IDUAL = 1, 2  -> G, D
C           i.e.
C           A_IDUAL = AB.AM(IFAC,IPRIM,IDUAL) * B_IPRIM + ...
C
C
C********** Dual RN
C
            FUNCEL = SURF * DFRO(1)
            RR.AM(IFAC,1,1) = -1.0D0 * FUNCEL / VOLG
            RR.AM(IFAC,1,2) = FUNCEL / VOLD
C
            FUNCEL = SURF * DFRO(2)
            RUX.AM(IFAC,1,1) = -1.0D0 * FUNCEL / VOLG
            RUX.AM(IFAC,1,2) = FUNCEL / VOLD
C
            FUNCEL = SURF * DFRO(3)
            RUY.AM(IFAC,1,1) = -1.0D0 * FUNCEL / VOLG
            RUY.AM(IFAC,1,2) = FUNCEL / VOLD
C
            FUNCEL = SURF * DFRO(4)
            RRET.AM(IFAC,1,1) = -1.0D0 * FUNCEL / VOLG
            RRET.AM(IFAC,1,2) = FUNCEL / VOLD
C
C********** Dual RUXN
C
            FUNCEL = SURF * (DFRUN(1) * CNX + DFRUT(1) * CTX)
            UXR.AM(IFAC,1,1) = -1.0D0 * FUNCEL / VOLG
            UXR.AM(IFAC,1,2) = FUNCEL / VOLD
C
            FUNCEL = SURF * (DFRUN(2) * CNX + DFRUT(2) * CTX)
            UXUX.AM(IFAC,1,1) = -1.0D0 * FUNCEL / VOLG
            UXUX.AM(IFAC,1,2) = FUNCEL / VOLD
C
            FUNCEL = SURF * (DFRUN(3) * CNX + DFRUT(3) * CTX)
            UXUY.AM(IFAC,1,1) = -1.0D0 * FUNCEL / VOLG
            UXUY.AM(IFAC,1,2) = FUNCEL / VOLD
C
            FUNCEL = SURF * (DFRUN(4) * CNX + DFRUT(4) * CTX)
            UXRET.AM(IFAC,1,1) = -1.0D0 * FUNCEL / VOLG
            UXRET.AM(IFAC,1,2) = FUNCEL / VOLD
C
C********** Dual RUYN
C
            FUNCEL = SURF * (DFRUN(1) * CNY + DFRUT(1) * CTY)
            UYR.AM(IFAC,1,1) = -1.0D0 * FUNCEL / VOLG
            UYR.AM(IFAC,1,2) = FUNCEL / VOLD
C
            FUNCEL = SURF * (DFRUN(2) * CNY + DFRUT(2) * CTY)
            UYUX.AM(IFAC,1,1) = -1.0D0 * FUNCEL / VOLG
            UYUX.AM(IFAC,1,2) = FUNCEL / VOLD
C
            FUNCEL = SURF * (DFRUN(3) * CNY + DFRUT(3) * CTY)
            UYUY.AM(IFAC,1,1) = -1.0D0 * FUNCEL / VOLG
            UYUY.AM(IFAC,1,2) = FUNCEL / VOLD
C
            FUNCEL = SURF * (DFRUN(4) * CNY + DFRUT(4) * CTY)
            UYRET.AM(IFAC,1,1) = -1.0D0 * FUNCEL / VOLG
            UYRET.AM(IFAC,1,2) = FUNCEL / VOLD
C
C********** Dual RETN
C
            FUNCEL = SURF * DFRET(1)
            RETR.AM(IFAC,1,1) = -1.0D0 * FUNCEL / VOLG
            RETR.AM(IFAC,1,2) = FUNCEL / VOLD
C
            FUNCEL = SURF * DFRET(2)
            RETUX.AM(IFAC,1,1) = -1.0D0 * FUNCEL / VOLG
            RETUX.AM(IFAC,1,2) = FUNCEL / VOLD
C
            FUNCEL = SURF * DFRET(3)
            RETUY.AM(IFAC,1,1) = -1.0D0 * FUNCEL / VOLG
            RETUY.AM(IFAC,1,2) = FUNCEL / VOLD
C
            FUNCEL = SURF * DFRET(4)
            RETRET.AM(IFAC,1,1) = -1.0D0 * FUNCEL / VOLG
            RETRET.AM(IFAC,1,2) = FUNCEL / VOLD
C
C
C********** La contribution de D
C           Je change l'orientation de normales
C           Dans ce cas F(UD,UG) = 0.5 * (F(UD)+F(UD)
C                                + RLAMMA * AMAT * (UD - UG))
C
            CNX = -1.0D0 * CNX
            CNY = -1.0D0 * CNY
            CTX = -1.0D0 * CTX
            CTY = -1.0D0 * CTY
C
C********** L'etat moyen au centre
C           Meme que avant; mais UNF, UTF change de signe
C
            UNF=-1.0D0*UNF
            UTF=-1.0D0*UTF
C
C********** On recalcule de la matrice de diffusivité pour le variables
C           conservatives (pas d'optimisation illisible!!!)
C
            CALL FRBM1(VINF, ROF,UNF,UTF,PF,GAMF,AMAT, RLAMMA)
C
C           Ici U=(\rho,\rho un, \rho ut, \rho et)
C           Il faut calculer la contribution sur F dans le repaire
C           (n,t) par rapport aux variables
C           U=(\rho,\rho ux, \rho uy, \rho et)
C
            DO I1 = 1,4,1
               BMAT(I1,1) = AMAT(I1,1)
               BMAT(I1,2) = AMAT(I1,2) * CNX + AMAT(I1,3) * CTX
               BMAT(I1,3) = AMAT(I1,2) * CNY + AMAT(I1,3) * CTY
               BMAT(I1,4) = AMAT(I1,4)
            ENDDO
C
            CALL CENTJ0(ROD,UXD,UYD,PD,GAMD,CNX,CNY,CTX,CTY,
     &           DFRO,DFRUN,DFRUT,DFRET)
C
C********** Centrée + Decentrée
C
            DO I1=1,4,1
               DFRO(I1) = DFRO(I1) + 0.5D0 * RLAMMA * BMAT(1,I1)
               DFRUN(I1) = DFRUN(I1) +  0.5D0 * RLAMMA * BMAT(2,I1)
               DFRUT(I1) = DFRUT(I1) +  0.5D0 * RLAMMA * BMAT(3,I1)
               DFRET(I1) = DFRET(I1) +  0.5D0 * RLAMMA * BMAT(4,I1)
            ENDDO
C
C
C********** Dual RN
C
            FUNCEL = SURF * DFRO(1)
            RR.AM(IFAC,2,2) = -1.0D0 * FUNCEL / VOLD
            RR.AM(IFAC,2,1) = FUNCEL / VOLG
C
            FUNCEL = SURF * DFRO(2)
            RUX.AM(IFAC,2,2) = -1.0D0 * FUNCEL / VOLD
            RUX.AM(IFAC,2,1) = FUNCEL / VOLG
C
            FUNCEL = SURF * DFRO(3)
            RUY.AM(IFAC,2,2) = -1.0D0 * FUNCEL / VOLD
            RUY.AM(IFAC,2,1) = FUNCEL / VOLG
C
            FUNCEL = SURF * DFRO(4)
            RRET.AM(IFAC,2,2) = -1.0D0 * FUNCEL / VOLD
            RRET.AM(IFAC,2,1) = FUNCEL / VOLG
C
C********** Dual RUXN
C
            FUNCEL = SURF * (DFRUN(1) * CNX + DFRUT(1) * CTX)
            UXR.AM(IFAC,2,2) = -1.0D0 * FUNCEL / VOLD
            UXR.AM(IFAC,2,1) = FUNCEL / VOLG
C
            FUNCEL = SURF * (DFRUN(2) * CNX + DFRUT(2) * CTX)
            UXUX.AM(IFAC,2,2) = -1.0D0 * FUNCEL / VOLD
            UXUX.AM(IFAC,2,1) = FUNCEL / VOLG
C
            FUNCEL = SURF * (DFRUN(3) * CNX + DFRUT(3) * CTX)
            UXUY.AM(IFAC,2,2) = -1.0D0 * FUNCEL / VOLD
            UXUY.AM(IFAC,2,1) = FUNCEL / VOLG
C
            FUNCEL = SURF * (DFRUN(4) * CNX + DFRUT(4) * CTX)
            UXRET.AM(IFAC,2,2) = -1.0D0 * FUNCEL / VOLD
            UXRET.AM(IFAC,2,1) = FUNCEL / VOLG
C
C********** Dual RUYN
C
            FUNCEL = SURF * (DFRUN(1) * CNY + DFRUT(1) * CTY)
            UYR.AM(IFAC,2,2) = -1.0D0 * FUNCEL / VOLD
            UYR.AM(IFAC,2,1) = FUNCEL / VOLG
C
            FUNCEL = SURF * (DFRUN(2) * CNY + DFRUT(2) * CTY)
            UYUX.AM(IFAC,2,2) = -1.0D0 * FUNCEL / VOLD
            UYUX.AM(IFAC,2,1) = FUNCEL / VOLG
C
            FUNCEL = SURF * (DFRUN(3) * CNY + DFRUT(3) * CTY)
            UYUY.AM(IFAC,2,2) = -1.0D0 * FUNCEL / VOLD
            UYUY.AM(IFAC,2,1) = FUNCEL / VOLG
C
            FUNCEL = SURF * (DFRUN(4) * CNY + DFRUT(4) * CTY)
            UYRET.AM(IFAC,2,2) = -1.0D0 * FUNCEL / VOLD
            UYRET.AM(IFAC,2,1) = FUNCEL / VOLG
C
C********** Dual RETN
C
            FUNCEL = SURF * DFRET(1)
            RETR.AM(IFAC,2,2) = -1.0D0 * FUNCEL / VOLD
            RETR.AM(IFAC,2,1) = FUNCEL / VOLG
C
            FUNCEL = SURF * DFRET(2)
            RETUX.AM(IFAC,2,2) = -1.0D0 * FUNCEL / VOLD
            RETUX.AM(IFAC,2,1) = FUNCEL / VOLG
C
            FUNCEL = SURF * DFRET(3)
            RETUY.AM(IFAC,2,2) = -1.0D0 * FUNCEL / VOLD
            RETUY.AM(IFAC,2,1) = FUNCEL / VOLG
C
            FUNCEL = SURF * DFRET(4)
            RETRET.AM(IFAC,2,2) = -1.0D0 * FUNCEL / VOLD
            RETRET.AM(IFAC,2,1) = FUNCEL / VOLG
C
         ELSE
C
C********** Murs (NGCG = NGCD)
C
C
C********** Les MELEMEs
C
            MELEDU.NUM(1,IFAC) = NGCG
            MELEDU.NUM(2,IFAC) = NGCD
            NLCG = MLENTC.LECT(NGCG)
C
            ROG = MPRN.VPOCHA(NLCG,1)
            PG = MPPN.VPOCHA(NLCG,1)
            UXG = MPUN.VPOCHA(NLCG,1)
            UYG = MPUN.VPOCHA(NLCG,2)
            GAMG = MPGAMN.VPOCHA(NLCG,1)
            VOLG = MPVOLU.VPOCHA(NLCG,1)
C
C********** Cut off
C
            VINF=MPUPRI.VPOCHA(NLCG,1)
            VINF=MAX(VINF,MPUINF.VPOCHA(NLCG,1))
C
C********** La normale sortante
C
            SURF = MPOVSU.VPOCHA(NLCF,1)
            CNX = MPNORM.VPOCHA(NLCF,1)
            CNY = MPNORM.VPOCHA(NLCF,2)
C
C********** L'etat moyen au centre
C
            UNF=0.0D0
            UTF=0.5D0*(UXG*CTX+UYG*CTY)
            ROF=ROG
            PF=PG
            GAMF=GAMG
C
C********** We include in the cut-off UG
C
            VINF=MAX(VINF,((UXG*UXG+UYG*UYG)**0.5D0))
C
C********** Calcule de la matrice de diffusivité pour le variables
C           conservatives.
C           Dans le repaire (n,t) local
C           F(UL,Ur) = 0.5 (F(UL)+F(UR)+AMAT*RLAMMA*(UL-UR))
C
            CALL FRBM1(VINF, ROF,UNF,UTF,PF,GAMF,AMAT, RLAMMA)
C
C           Ici U=(\rho,\rho un, \rho ut, \rho et)
C           Il faut calculer la contribution sur F dans le repaire
C           (n,t) par rapport aux variables
C           U=(\rho,\rho ux, \rho uy, \rho et)
C           Dans le cas murs, le seul flux normal donne une contribution
C           non nulle
C           En plus AMAT(1,2)=AMAT(3,2)=AMAT(4,2)=0
C           AMAT(2,2)=1.0D0
C
C            write(*,*) 'AMAT(*,2)=',(AMAT(I1,2),I1=1,4)
C
C********** Contribution centrée
C
            CALL CENTJ2(ROG,UXG,UYG,GAMG,CNX,CNY,
     &           DFRUN)
C
C********** Centrée + Decentrée
C
            DFRUN(2) = DFRUN(2) +  (RLAMMA * CNX)
            DFRUN(3) = DFRUN(3) +  (RLAMMA * CNY)
C
C********** Dual RN
C
            RR.AM(IFAC,1,1) = 0.0D0
            RR.AM(IFAC,1,2) = 0.0D0
C
            RUX.AM(IFAC,1,1) = 0.0D0
            RUX.AM(IFAC,1,2) = 0.0D0
C
            RUY.AM(IFAC,1,1) = 0.0D0
            RUY.AM(IFAC,1,2) = 0.0D0
C
            RRET.AM(IFAC,1,1) = 0.0D0
            RRET.AM(IFAC,1,2) = 0.0D0
C
C********** Dual RUXN
C
            FUNCEL = SURF * DFRUN(1) * CNX
            UXR.AM(IFAC,1,1) = -1.0D0 * FUNCEL / VOLG
            UXR.AM(IFAC,1,2) = 0.0D0
C
            FUNCEL = SURF * DFRUN(2) * CNX
            UXUX.AM(IFAC,1,1) = -1.0D0 * FUNCEL / VOLG
            UXUX.AM(IFAC,1,2) = 0.0D0
C
            FUNCEL = SURF * DFRUN(3) * CNX
            UXUY.AM(IFAC,1,1) = -1.0D0 * FUNCEL / VOLG
            UXUY.AM(IFAC,1,2) = 0.0D0
C
            FUNCEL = SURF * DFRUN(4) * CNX
            UXRET.AM(IFAC,1,1) = -1.0D0 * FUNCEL / VOLG
            UXRET.AM(IFAC,1,2) = 0.0D0
C
C********** Dual RUYN
C
            FUNCEL = SURF * DFRUN(1) * CNY
            UYR.AM(IFAC,1,1) = -1.0D0 * FUNCEL / VOLG
            UYR.AM(IFAC,1,2) = 0.0D0
C
            FUNCEL = SURF * DFRUN(2) * CNY
            UYUX.AM(IFAC,1,1) = -1.0D0 * FUNCEL / VOLG
            UYUX.AM(IFAC,1,2) = 0.0D0
C
            FUNCEL = SURF * DFRUN(3) * CNY
            UYUY.AM(IFAC,1,1) = -1.0D0 * FUNCEL / VOLG
            UYUY.AM(IFAC,1,2) = 0.0D0
C
            FUNCEL = SURF * DFRUN(4) * CNY
            UYRET.AM(IFAC,1,1) = -1.0D0 * FUNCEL / VOLG
            UYRET.AM(IFAC,1,2) = 0.0D0
C
C********** Dual RETN
C
            RETR.AM(IFAC,1,1) = 0.0D0
            RETR.AM(IFAC,1,2) = 0.0D0
C
            RETUX.AM(IFAC,1,1) = 0.0D0
            RETUX.AM(IFAC,1,2) = 0.0D0
C
            RETUY.AM(IFAC,1,1) = 0.0D0
            RETUY.AM(IFAC,1,2) = 0.0D0
C
            RETRET.AM(IFAC,1,1) = 0.0D0
            RETRET.AM(IFAC,1,2) = 0.0D0
C
C********** Dual RN
C
            RR.AM(IFAC,2,2) = 0.0D0
            RR.AM(IFAC,2,1) = 0.0D0
C
            RUX.AM(IFAC,2,2) = 0.0D0
            RUX.AM(IFAC,2,1) = 0.0D0
C
            RUY.AM(IFAC,2,2) = 0.0D0
            RUY.AM(IFAC,2,1) = 0.0D0
C
            RRET.AM(IFAC,2,2) = 0.0D0
            RRET.AM(IFAC,2,1) = 0.0D0
C
C********** Dual RUXN
C
            UXR.AM(IFAC,2,2) = 0.0D0
            UXR.AM(IFAC,2,1) = 0.0D0
C
            UXUX.AM(IFAC,2,2) = 0.0D0
            UXUX.AM(IFAC,2,1) = 0.0D0
C
            UXUY.AM(IFAC,2,2) = 0.0D0
            UXUY.AM(IFAC,2,1) = 0.0D0
C
            UXRET.AM(IFAC,2,2) = 0.0D0
            UXRET.AM(IFAC,2,1) = 0.0D0
C
C********** Dual RUYN
C
            UYR.AM(IFAC,2,2) = 0.0D0
            UYR.AM(IFAC,2,1) = 0.0D0
C
            UYUX.AM(IFAC,2,2) = 0.0D0
            UYUX.AM(IFAC,2,1) = 0.0D0
C
            UYUY.AM(IFAC,2,2) = 0.0D0
            UYUY.AM(IFAC,2,1) = 0.0D0
C
            UYRET.AM(IFAC,2,2) = 0.0D0
            UYRET.AM(IFAC,2,1) = 0.0D0
C
C********** Dual RETN
C
            RETR.AM(IFAC,2,2) = 0.0D0
            RETR.AM(IFAC,2,1) = 0.0D0
C
            RETUX.AM(IFAC,2,2) = 0.0D0
            RETUX.AM(IFAC,2,1) = 0.0D0
C
            RETUY.AM(IFAC,2,2) = 0.0D0
            RETUY.AM(IFAC,2,1) = 0.0D0
C
            RETRET.AM(IFAC,2,2) = 0.0D0
            RETRET.AM(IFAC,2,1) = 0.0D0
C
         ENDIF
      ENDDO
C
      SEGDES MELEMC
      SEGDES MELEFE
      SEGDES MELEMF
C
      SEGDES MPOVSU
      SEGDES MPVOLU
      SEGDES MPNORM
C
      SEGDES MPRN
      SEGDES MPPN
      SEGDES MPUN
      SEGDES MPGAMN
C
      SEGDES MELEDU
      SEGDES MATRIK
      SEGDES IMATRI
C
      SEGDES   RR  , RUX , RUY , RRET ,
     &         UXR , UXUX , UXUY , UXRET ,
     &         UYR , UYUX , UYUY , UYRET ,
     &         RETR , RETUX , RETUY , RETRET
 
      SEGSUP MLENTC
      SEGSUP MLENTF
      SEGDES MLMINC
      SEGSUP MLELIM
C
      IF(MELLIM .NE.0) SEGDES MELLIM
C
      SEGDES MPUPRI
      SEGDES MPUINF
C
 9999 CONTINUE
      RETURN
      END