Télécharger ylap2e.eso

Retour à la liste

ylap2e
C YLAP2E    SOURCE    CB215821  20/11/25    13:44:20     10792          
      SUBROUTINE YLAP2E(ICOGRV,MPGRVF,MPROC,MPVITC,
     $     MPVOLU,MPNORM,MPSURF,MELEFL,
     $     KRFACE,KRCENT,
     $     LCLIMV,KRVIMP,MPVIMP,
     $     LCLITO,KRTOIM,
     $     NOMINC,
     $     MU,
     $     IJACO,
     $     IMPR,IRET)
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8 (A-H,O-Z)
C***********************************************************************
C NOM         : YLAP2E
C DESCRIPTION : Calcul de la matrice jacobienne du résidu du laplacien
C               VF 3D.
C               Ici, on ne calcule que des contributions 'compliquées' à
C               la matrice jacobienne faisant intervenir les
C               coefficients pour le calcul des gradients de vitesse
C               (ICOGRV)
C               (contributions à (d Res_{\rho e_t} / d var)
C               var prenant successivement les valeurs :
C               \rho,   \rho u,    \rho v,     \rho w,    \rho e_t )
C               Les contributions sont plus "compliquées" à calculer que
C               les simples (cf. ylap1c) car on a à dériver des produits
C               de fonctions de la vitesse.
C               ylap2d calcule aussi une partie des contributions
C               'compliquées'.
C
C NOTE PERSO  : On pourrait programmer ça plus lisiblement en stockant
C               les contributions dans un tableau de pointeurs (2
C               indices, c'est possible ?) et en effectuant des produits
C               matriciels (coeff. x matrices de dérivées).
C               On n'y coupera pas si on veut regrouper 2D et 3D...
C               On ne va pas le faire.
C
C LANGAGE     : ESOPE
C AUTEUR      : Stéphane GOUNAND (CEA/DEN/DM2S/SFME/LTMF)
C               mél : gounand@semt2.smts.cea.fr
C***********************************************************************
C APPELES (UTIL)   : AJMTK : ajoute un objet de type MATSIM (non
C                            standard) à un objet de type MATRIK.
C APPELE PAR       : YLAP2A : Calcul de la matrice jacobienne du
C                    résidu du laplacien VF 3D.
C***********************************************************************
C ENTREES            : ICOGRV (type MCHELM) : coefficients pour le
C                        calcul du gradient de la vitesse aux
C                        interfaces.
C                      MPGRVF (type MPOVAL) : gradient de la vitesse
C                        aux interfaces.
C                      MPROC (type MPOVAL) : masse volumique par
C                                            élément.
C                      MPVITC (type MPOVAL) : vitesse par élément.
C                      MPVOLU (type MPOVAL) : volume des éléments.
C                      MPNORM (type MPOVAL) : normale aux faces.
C                      MPSURF (type MPOVAL) : surface des faces.
C                      MELEFL (type MELEME) : connectivités face-(centre
C                                             gauche, centre droit).
C                      KRFACE (type MLENTI) : tableau de repérage dans
C                        le maillage des faces des éléments.
C                      KRCENT (type MLENTI) : tableau de repérage dans
C                        le maillage des centres des éléments.
C                      LCLIMV (type logique) : .TRUE. => CL de Dirichlet
C                        sur la vitesse.
C                      KRVIMP (type MLENTI) : tableau de repérage dans
C                        maillage des CL de Dirichlet sur la
C                        vitesse.
C                      MPVIMP (type MPOVAL) : valeurs des CL de
C                        Dirichlet sur la vitesse.
C                      LCLITO (type logique) : .TRUE. => CL de Dirichlet
C                        sur le tenseur des contraintes.
C                      KRTOIM (type MLENTI) : tableau de repérage dans
C                        maillage des CL de Dirichlet sur le tenseur des
C                        contraintes
C                      NOMINC (type MLMOTS) : noms des inconnues.
C                      MU (type réel) : viscosité dynamique (SI).
C ENTREES/SORTIES    : IJACO (type MATRIK) : matrice jacobienne du
C                        résidu du laplacien VF 3D.
C SORTIES            : -
C CODE RETOUR (IRET) : = 0 si tout s'est bien passé
C***********************************************************************
C VERSION    : v1, 28/08/2001, version initiale
C HISTORIQUE : v1, 28/08/2001, création
C HISTORIQUE :
C HISTORIQUE :
C***********************************************************************
C Prière de PRENDRE LE TEMPS de compléter les commentaires
C en cas de modification de ce sous-programme afin de faciliter
C la maintenance !
C***********************************************************************
 
-INC PPARAM
-INC CCOPTIO
-INC SMCOORD
-INC SMCHPOI
      POINTEUR MPGRVF.MPOVAL
      POINTEUR MPROC.MPOVAL ,MPVITC.MPOVAL
      POINTEUR MPSURF.MPOVAL,MPNORM.MPOVAL,MPVOLU.MPOVAL
      POINTEUR MPVIMP.MPOVAL
-INC SMCHAML
      POINTEUR ICOGRV.MCHELM,JCOGRV.MCHAML
      POINTEUR KDUNDX.MELVAL,KDUNDY.MELVAL,KDUNDZ.MELVAL
-INC SMELEME
      POINTEUR MELEFL.MELEME
      POINTEUR MCOGRV.MELEME
-INC SMLENTI
      POINTEUR KRVIMP.MLENTI,KRTOIM.MLENTI
      POINTEUR KRCENT.MLENTI,KRFACE.MLENTI
-INC SMLMOTS
      POINTEUR NOMINC.MLMOTS
-INC SMMATRIK
      POINTEUR IJACO.MATRIK
*
* Objet matrice élémentaire simplifié
*
      SEGMENT GMATSI
      INTEGER POIPR1(NPP1,NEL1)
      INTEGER POIDU1(1,NEL1)
      INTEGER POIPR2(NPP2,NEL2)
      INTEGER POIDU2(2,NEL2)
      POINTEUR LMATSI(0).MATSIM
      ENDSEGMENT
*  Contributions compliquées 1 de la part du gradient de
*  vitesse (CCGRV1)
      POINTEUR CCGRV1.GMATSI
      SEGMENT MATSIM
      CHARACTER*8 NOMPRI,NOMDUA
      REAL*8 VALMA1(1,NPP1,NEL1)
      REAL*8 VALMA2(2,NPP2,NEL2)
      ENDSEGMENT
      POINTEUR RETRHO.MATSIM
      POINTEUR RETROU.MATSIM
      POINTEUR RETROV.MATSIM
      POINTEUR RETROW.MATSIM
*
      REAL*8 MU
*
      INTEGER IMPR,IRET
*
      LOGICAL LCLIMV,LCLITO
      LOGICAL LMUR
      LOGICAL LCTR1A,LCTR1B
*
      INTEGER IELEM,IPD,IPP,ISOUCH,IEL1,IEL2
      INTEGER NELEM,NPD,NPP,NSOUCH,NEL1,NEL2,NPP1,NPP2
      INTEGER NGCDRO,NGCGAU,NGFACE,NPPRIM,NPDUAL
      INTEGER NLCENP,NLCEND,NLFACE,NLCLTO,NLCLV
      INTEGER NLCED,NLCEG,NLFVI
      INTEGER NPTEL
      INTEGER ICOORX
*
      REAL*8 ALPHA,UMALPH
      REAL*8 DRD,DRG
      REAL*8 DUXXF,DUXYF,DUXZF
      REAL*8 DUYXF,DUYYF,DUYZF
      REAL*8 DUZXF,DUZYF,DUZZF
      REAL*8 UXD,UXF,UXG
      REAL*8 UYD,UYF,UYG
      REAL*8 UZD,UZF,UZG
      REAL*8 XD,XF,XG,XFMXD,XFMXG
      REAL*8 YD,YF,YG,YFMYD,YFMYG
      REAL*8 ZD,ZF,ZG,ZFMZD,ZFMZG
      REAL*8 ALPHAX,ALPHAY,ALPHAZ
      REAL*8 CNX,CNY,CNZ
      REAL*8 SIGNOR,SURFFA,VOLUEL
      REAL*8 RHOP,UP,VP,WP
      REAL*8 FACTOR,FACTDU,FACTDV,FACTDW
      REAL*8 DUDRHO,DUDROU
      REAL*8 DVDRHO,DVDROV
      REAL*8 DWDRHO,DWDROW
*
* Executable statements
*
      IF (IMPR.GT.2) WRITE(IOIMP,*) 'Entrée dans ylap2e.eso'
*   On calcule les contributions à (d Res_{\rho e_t} / d var) ; var
*   prenant successivement les valeurs :
*      \rho,   \rho u,    \rho v,      \rho w,     \rho e_t
*   On dérive les termes :  (\tens{\tau(\grad{u})} \prod \vect{u})
*                            \pscal \vect{n}
*   ce qui donne deux contributions.
*   Contribution 1 :
*     ( (d \tens{\tau} / d var) \prod \vect{u}) \pscal \vect{n}
*   Contribution 2 :
*     ( \tens{\tau}  \prod (d \vect{u} / d var)) \pscal \vect{n}
*   Note :
*     pas de contribution à (d Res_{\rho e_t} / d \rho e_t).
*   Les noms de matrices élémentaires (type MATSIM) associées sont :
*   RETRHO, RETROU, RETROV, RETROW
      IF (LCLIMV) THEN
         SEGACT KRVIMP
         SEGACT MPVIMP
      ENDIF
      IF (LCLITO) THEN
         SEGACT KRTOIM
      ENDIF
      SEGACT NOMINC
      SEGACT KRCENT
      SEGACT KRFACE
      SEGACT MELEFL
      SEGACT MPSURF
      SEGACT MPNORM
      SEGACT MPVOLU
      SEGACT MPROC
      SEGACT MPVITC
      SEGACT MPGRVF
      SEGACT ICOGRV
      NSOUCH=ICOGRV.IMACHE(/1)
      DO 1 ISOUCH=1,NSOUCH
         MCOGRV=ICOGRV.IMACHE(ISOUCH)
         JCOGRV=ICOGRV.ICHAML(ISOUCH)
         SEGACT JCOGRV
         KDUNDX=JCOGRV.IELVAL(1)
         KDUNDY=JCOGRV.IELVAL(2)
         KDUNDZ=JCOGRV.IELVAL(3)
         SEGDES JCOGRV
         SEGACT KDUNDX
         SEGACT KDUNDY
         SEGACT KDUNDZ
         SEGACT MCOGRV
         NELEM=MCOGRV.NUM(/2)
         NPTEL=MCOGRV.NUM(/1)
         NPP1=NPTEL-1
         NPP2=NPTEL-1
         NEL1=NELEM
         NEL2=NELEM
         IEL1=1
         IEL2=1
         SEGINI RETRHO
         SEGINI RETROU
         SEGINI RETROV
         SEGINI RETROW
         SEGINI CCGRV1
         RETRHO.NOMPRI(1:4)=NOMINC.MOTS(1)
         RETRHO.NOMPRI(5:8)='    '
         RETRHO.NOMDUA(1:4)=NOMINC.MOTS(5)
         RETRHO.NOMDUA(5:8)='    '
         RETROU.NOMPRI(1:4)=NOMINC.MOTS(2)
         RETROU.NOMPRI(5:8)='    '
         RETROU.NOMDUA(1:4)=NOMINC.MOTS(5)
         RETROU.NOMDUA(5:8)='    '
         RETROV.NOMPRI(1:4)=NOMINC.MOTS(3)
         RETROV.NOMPRI(5:8)='    '
         RETROV.NOMDUA(1:4)=NOMINC.MOTS(5)
         RETROV.NOMDUA(5:8)='    '
         RETROW.NOMPRI(1:4)=NOMINC.MOTS(4)
         RETROW.NOMPRI(5:8)='    '
         RETROW.NOMDUA(1:4)=NOMINC.MOTS(5)
         RETROW.NOMDUA(5:8)='    '
         DO 12 IELEM=1,NELEM
*     Le premier point du support de ICOGRV est un point FACE
            NGFACE=MCOGRV.NUM(1,IELEM)
            NLFACE=KRFACE.LECT(NGFACE)
            IF (NLFACE.EQ.0) THEN
               WRITE(IOIMP,*) 'Erreur de programmation n°1'
               GOTO 9999
            ENDIF
*     On calcule la contribution 1 à la matrice jacobienne IJACO de la
*     face NGFACE
*           (points duaux     : centres à gauche et à droite de la face)
*           (points primaux   : une partie (bicoz conditions aux limites)
*            de ceux du stencil pour le calcul du gradient
*            à la face, ils doivent être des points centres)
*       Si le tenseur des contraintes sur la face est imposé par les
*       conditions aux limites, la contribution 1 de la face à IJACO est
*       nulle.
            LCTR1A=.TRUE.
            IF (LCLITO) THEN
               NLCLTO=KRTOIM.LECT(NGFACE)
               IF (NLCLTO.NE.0) THEN
                  LCTR1A=.FALSE.
               ENDIF
            ENDIF
            IF (LCTR1A) THEN
               NGCGAU=MELEFL.NUM(1,NLFACE)
               NGCDRO=MELEFL.NUM(3,NLFACE)
               LMUR=(NGCGAU.EQ.NGCDRO)
*     On calcule tout d'abord une interpolation de la vitesse sur la
*     face NGFACE de la même manière que dans ylap13.eso.
               IF (LCLIMV) THEN
                  NLFVI=KRVIMP.LECT(NGFACE)
               ELSE
                  NLFVI=0
               ENDIF
*       La vitesse est imposée sur la face, rien à calculer
               IF (NLFVI.NE.0) THEN
                  UXF=MPVIMP.VPOCHA(NLFVI,1)
                  UYF=MPVIMP.VPOCHA(NLFVI,2)
                  UZF=MPVIMP.VPOCHA(NLFVI,3)
               ELSE
                  NLCEG=KRCENT.LECT(NGCGAU)
                  NLCED=KRCENT.LECT(NGCDRO)
*         Cas non au bord
                  IF (.NOT.LMUR) THEN
*           Paramètres géométriques
                     ICOORX = ((IDIM + 1) * (NGFACE - 1))+1
                     XF = MCOORD.XCOOR(ICOORX)
                     YF = MCOORD.XCOOR(ICOORX+1)
                     ZF = MCOORD.XCOOR(ICOORX+2)
                     ICOORX = ((IDIM + 1) * (NGCGAU - 1))+1
                     XG = MCOORD.XCOOR(ICOORX)
                     YG = MCOORD.XCOOR(ICOORX+1)
                     ZG = MCOORD.XCOOR(ICOORX+2)
                     XFMXG = XF - XG
                     YFMYG = YF - YG
                     ZFMZG = ZF - ZG
                     DRG=SQRT((XFMXG*XFMXG)+(YFMYG*YFMYG)+(ZFMZG*ZFMZG))
                     ICOORX = ((IDIM + 1) * (NGCDRO - 1))+1
                     XD = MCOORD.XCOOR(ICOORX)
                     YD = MCOORD.XCOOR(ICOORX+1)
                     ZD = MCOORD.XCOOR(ICOORX+2)
                     XFMXD = XF - XD
                     YFMYD = YF - YD
                     ZFMZD = ZF - ZD
                     DRD=SQRT((XFMXD*XFMXD)+(YFMYD*YFMYD)+(ZFMZD*ZFMZD))
                     ALPHA=DRG/(DRG+DRD)
                     UMALPH= 1.0D0 - ALPHA
                     DUXXF = MPGRVF.VPOCHA(NLFACE,1)
                     DUXYF = MPGRVF.VPOCHA(NLFACE,2)
                     DUXZF = MPGRVF.VPOCHA(NLFACE,3)
                     DUYXF = MPGRVF.VPOCHA(NLFACE,4)
                     DUYYF = MPGRVF.VPOCHA(NLFACE,5)
                     DUYZF = MPGRVF.VPOCHA(NLFACE,6)
                     DUZXF = MPGRVF.VPOCHA(NLFACE,7)
                     DUZYF = MPGRVF.VPOCHA(NLFACE,8)
                     DUZZF = MPGRVF.VPOCHA(NLFACE,9)
*            Calcul de la vitesse
                     UXG = MPVITC.VPOCHA(NLCEG,1)
                     UYG = MPVITC.VPOCHA(NLCEG,2)
                     UZG = MPVITC.VPOCHA(NLCEG,3)
                     UXD = MPVITC.VPOCHA(NLCED,1)
                     UYD = MPVITC.VPOCHA(NLCED,2)
                     UZD = MPVITC.VPOCHA(NLCED,3)
                     UXF = UMALPH * UXG + ALPHA * UXD
                     UYF = UMALPH * UYG + ALPHA * UYD
                     UZF = UMALPH * UZG + ALPHA * UZD
*            Correction de la vitesse lineaire exacte
                     UXF = UXF +
     &              (DUXXF * ((XFMXG * UMALPH)+ (XFMXD * ALPHA))) +
     &              (DUXYF * ((YFMYG * UMALPH)+ (YFMYD * ALPHA))) +
     &              (DUXZF * ((ZFMZG * UMALPH)+ (ZFMZD * ALPHA)))
                     UYF = UYF +
     &              (DUYXF * ((XFMXG * UMALPH)+ (XFMXD * ALPHA))) +
     &              (DUYYF * ((YFMYG * UMALPH)+ (YFMYD * ALPHA))) +
     &              (DUYZF * ((ZFMZG * UMALPH)+ (ZFMZD * ALPHA)))
                     UZF = UZF +
     &              (DUZXF * ((XFMXG * UMALPH)+ (XFMXD * ALPHA))) +
     &              (DUZYF * ((YFMYG * UMALPH)+ (YFMYD * ALPHA))) +
     &              (DUZZF * ((ZFMZG * UMALPH)+ (ZFMZD * ALPHA)))
 
*         Cas au bord
                  ELSE
*          Parametres geometriques
                     ICOORX = ((IDIM + 1) * (NGFACE - 1))+1
                     XF = MCOORD.XCOOR(ICOORX)
                     YF = MCOORD.XCOOR(ICOORX+1)
                     ZF = MCOORD.XCOOR(ICOORX+2)
                     ICOORX = ((IDIM + 1) * (NGCGAU - 1))+1
                     XG = MCOORD.XCOOR(ICOORX)
                     YG = MCOORD.XCOOR(ICOORX+1)
                     ZG = MCOORD.XCOOR(ICOORX+2)
                     XFMXG = XF - XG
                     YFMYG = YF - YG
                     ZFMZG = ZF - ZG
                     DUXXF = MPGRVF.VPOCHA(NLFACE,1)
                     DUXYF = MPGRVF.VPOCHA(NLFACE,2)
                     DUXZF = MPGRVF.VPOCHA(NLFACE,3)
                     DUYXF = MPGRVF.VPOCHA(NLFACE,4)
                     DUYYF = MPGRVF.VPOCHA(NLFACE,5)
                     DUYZF = MPGRVF.VPOCHA(NLFACE,6)
                     DUZXF = MPGRVF.VPOCHA(NLFACE,7)
                     DUZYF = MPGRVF.VPOCHA(NLFACE,8)
                     DUZZF = MPGRVF.VPOCHA(NLFACE,9)
*            Calcul de la vitesse
                     UXF = MPVITC.VPOCHA(NLCEG,1)
                     UYF = MPVITC.VPOCHA(NLCEG,2)
                     UZF = MPVITC.VPOCHA(NLCEG,3)
*            Correction de la vitesse lineaire exacte
                     UXF = UXF +
     &                    (DUXXF * XFMXG ) +
     &                    (DUXYF * YFMYG ) +
     &                    (DUXZF * ZFMZG )
                     UYF = UYF +
     &                    (DUYXF * XFMXG ) +
     &                    (DUYYF * YFMYG ) +
     &                    (DUYZF * ZFMZG )
                     UZF = UZF +
     &                    (DUZXF * XFMXG ) +
     &                    (DUZYF * YFMYG ) +
     &                    (DUZZF * ZFMZG )
                  ENDIF
               ENDIF
*     On distingue le cas où la face est un bord du maillage (mur)
*                  du cas où la face est interne au maillage
               IF (.NOT.LMUR) THEN
                  NPD=2
               ELSE
                  NPD=1
               ENDIF
               NPP=NPTEL-1
*     IPD=1 : point à gauche du point NGFACE
*     IPD=2 : point à droite du point NGFACE
               DO 122 IPD=1,NPD
                  NPDUAL=MELEFL.NUM((2*IPD)-1,NLFACE)
                  IF (.NOT.LMUR) THEN
                     CCGRV1.POIDU2(IPD,IEL2)=NPDUAL
                  ELSE
                     CCGRV1.POIDU1(IPD,IEL1)=NPDUAL
                  ENDIF
                  NLCEND=KRCENT.LECT(NPDUAL)
                  IF (NLCEND.EQ.0) THEN
                     WRITE(IOIMP,*) 'Erreur grave n°1'
                     GOTO 9999
                  ENDIF
                  DO 124 IPP=1,NPP
                     NPPRIM=MCOGRV.NUM(IPP+1,IELEM)
                     NLCENP=KRCENT.LECT(NPPRIM)
                     IF(NLCENP .EQ. 0)THEN
*     Lorsque une contribution est nulle, on fixe artificiellement le
*     point primal égal au point dual.
                        IF (.NOT.LMUR) THEN
                           CCGRV1.POIPR2(IPP,IEL2)=NPDUAL
                           RETRHO.VALMA2(IPD,IPP,IEL2)=0.D0
                           RETROU.VALMA2(IPD,IPP,IEL2)=0.D0
                           RETROV.VALMA2(IPD,IPP,IEL2)=0.D0
                           RETROW.VALMA2(IPD,IPP,IEL2)=0.D0
                        ELSE
                           CCGRV1.POIPR1(IPP,IEL1)=NPDUAL
                           RETRHO.VALMA1(IPD,IPP,IEL1)=0.D0
                           RETROU.VALMA1(IPD,IPP,IEL1)=0.D0
                           RETROV.VALMA1(IPD,IPP,IEL1)=0.D0
                           RETROW.VALMA1(IPD,IPP,IEL1)=0.D0
                        ENDIF
                     ELSE
*     Les contributions 1 valent :
*     (d Res_{\rho e_t})_d / (d var)_p =
*     +/-1 (normale sortante, rentrante)  (1/V_d) * (S_f) * \mu
*     * [ [ [ (n_x * (( 4/3 * u_f * \alpha_x) + (v_f * \alpha_y)
*                                             + (w_f * \alpha_z)))
*           + (n_y * ((u_f * \alpha_y) + (-2/3 * v_f * \alpha_x)))
*           + (n_z * ((u_f * \alpha_z) + (-2/3 * w_f * \alpha_x)))
*           ]
*         * ((du)_p / (d var)_p)
*         ]
*       + [ [ (n_x * ((-2/3 * u_f * \alpha_y) + (v_f * \alpha_x)))
*           + (n_y * ((u_f * \alpha_x) + ( 4/3 * v_f * \alpha_y)
*                                      +        (w_f * \alpha_z)))
*           + (n_z * ((v_f * \alpha_z) + (-2/3 * w_f * \alpha_y)))
*           ]
*         * ((dv)_p / (d var)_p)
*         ]
*       + [ [ (n_x * ((-2/3 * u_f * \alpha_z) + (w_f * \alpha_x)))
*           + (n_y * ((-2/3 * v_f * \alpha_z) + (w_f * \alpha_y)))
*           + (n_z * ((u_f * \alpha_x) +        (v_f * \alpha_y)
*                                      + ( 4/3 * w_f * \alpha_z)))
*           ]
*         * ((dw)_p / (d var)_p)
*         ]
*       ]
*
*       normale sortante pour IPD=1, rentrante pour IPD=2
                        SIGNOR=(-1.D0)**(IPD+1)
                        VOLUEL=MPVOLU.VPOCHA(NLCEND,1)
                        SURFFA=MPSURF.VPOCHA(NLFACE,1)
                        CNX   =MPNORM.VPOCHA(NLFACE,1)
                        CNY   =MPNORM.VPOCHA(NLFACE,2)
                        CNZ   =MPNORM.VPOCHA(NLFACE,3)
                        ALPHAX=KDUNDX.VELCHE(IPP+1,IELEM)
                        ALPHAY=KDUNDY.VELCHE(IPP+1,IELEM)
                        ALPHAZ=KDUNDZ.VELCHE(IPP+1,IELEM)
                        RHOP  =MPROC.VPOCHA(NLCENP,1)
                        UP    =MPVITC.VPOCHA(NLCENP,1)
                        VP    =MPVITC.VPOCHA(NLCENP,2)
                        WP    =MPVITC.VPOCHA(NLCENP,3)
                        FACTOR=SIGNOR*(1.D0/VOLUEL)*SURFFA*MU
                        FACTDU=(CNX*((( 4.D0/3.D0)*UXF*ALPHAX)
     $                              +(UYF*ALPHAY)
     $                              +(UZF*ALPHAZ)))
     $                       + (CNY*((UXF*ALPHAY)
     $                              +((-2.D0/3.D0)*UYF*ALPHAX)))
     $                       + (CNZ*((UXF*ALPHAZ)
     $                              +((-2.D0/3.D0)*UZF*ALPHAX)))
                        FACTDV=(CNX*(((-2.D0/3.D0)*UXF*ALPHAY)
     $                              +(UYF*ALPHAX)))
     $                       + (CNY*((UXF*ALPHAX)
     $                              +(( 4.D0/3.D0)*UYF*ALPHAY)
     $                              +(UZF*ALPHAZ)))
     $                       + (CNZ*((UYF*ALPHAZ)
     $                              +((-2.D0/3.D0)*UZF*ALPHAY)))
                        FACTDW=(CNX*(((-2.D0/3.D0)*UXF*ALPHAZ)
     $                              +(UZF*ALPHAX)))
     $                       + (CNY*(((-2.D0/3.D0)*UYF*ALPHAZ)
     $                              +(UZF*ALPHAY)))
     $                       + (CNZ*((UXF*ALPHAX)
     $                              +(UYF*ALPHAY)
     $                              +(( 4.D0/3.D0)*UZF*ALPHAZ)))
                        DUDRHO=-UP /RHOP
                        DUDROU=1.D0/RHOP
                        DVDRHO=-VP /RHOP
                        DVDROV=1.D0/RHOP
                        DWDRHO=-WP /RHOP
                        DWDROW=1.D0/RHOP
                        IF (.NOT.LMUR) THEN
                           CCGRV1.POIPR2(IPP,IEL2)=NPPRIM
                           RETRHO.VALMA2(IPD,IPP,IEL2)=
     $                          FACTOR*( (FACTDU*DUDRHO)
     $                                  +(FACTDV*DVDRHO)
     $                                  +(FACTDW*DWDRHO))
                           RETROU.VALMA2(IPD,IPP,IEL2)=
     $                          FACTOR*  (FACTDU*DUDROU)
                           RETROV.VALMA2(IPD,IPP,IEL2)=
     $                          FACTOR*  (FACTDV*DVDROV)
                           RETROW.VALMA2(IPD,IPP,IEL2)=
     $                          FACTOR*  (FACTDW*DWDROW)
                        ELSE
                           CCGRV1.POIPR1(IPP,IEL1)=NPPRIM
                           RETRHO.VALMA1(IPD,IPP,IEL1)=
     $                          FACTOR*( (FACTDU*DUDRHO)
     $                                  +(FACTDV*DVDRHO)
     $                                  +(FACTDW*DWDRHO))
                           RETROU.VALMA1(IPD,IPP,IEL1)=
     $                          FACTOR*  (FACTDU*DUDROU)
                           RETROV.VALMA1(IPD,IPP,IEL1)=
     $                          FACTOR*  (FACTDV*DVDROV)
                           RETROW.VALMA1(IPD,IPP,IEL1)=
     $                          FACTOR*  (FACTDW*DWDROW)
                        ENDIF
                     ENDIF
 124              CONTINUE
 122           CONTINUE
               IF (.NOT.LMUR) THEN
                  IEL2=IEL2+1
               ELSE
                  IEL1=IEL1+1
               ENDIF
            ENDIF
 12      CONTINUE
         NPP1=NPTEL-1
         NPP2=NPTEL-1
         NEL1=IEL1-1
         NEL2=IEL2-1
         SEGADJ RETRHO
         SEGADJ RETROU
         SEGADJ RETROV
         SEGADJ RETROW
         SEGADJ CCGRV1
         CCGRV1.LMATSI(**)=RETRHO
         CCGRV1.LMATSI(**)=RETROU
         CCGRV1.LMATSI(**)=RETROV
         CCGRV1.LMATSI(**)=RETROW
* On accumule les matrices résultantes dans IJACO
         CALL AJMTK(CCGRV1,IJACO,IMPR,IRET)
         IF (IRET.NE.0) GOTO 9999
         SEGSUP RETRHO
         SEGSUP RETROU
         SEGSUP RETROV
         SEGSUP RETROW
         SEGSUP CCGRV1
*
         SEGDES MCOGRV
         SEGDES KDUNDZ
         SEGDES KDUNDY
         SEGDES KDUNDX
 1    CONTINUE
      SEGDES ICOGRV
      SEGDES MPGRVF
      SEGDES MPVITC
      SEGDES MPROC
      SEGDES MPVOLU
      SEGDES MPNORM
      SEGDES MPSURF
      SEGDES MELEFL
      SEGDES KRFACE
      SEGDES KRCENT
      SEGDES NOMINC
      IF (LCLITO) THEN
         SEGDES KRTOIM
      ENDIF
      IF (LCLIMV) THEN
         SEGDES KRVIMP
         SEGDES MPVIMP
      ENDIF
*
* Normal termination
*
      IRET=0
      RETURN
*
* Format handling
*
*
* Error handling
*
 9999 CONTINUE
      IRET=1
      WRITE(IOIMP,*) 'An error was detected in subroutine ylap2e'
      RETURN
*
* End of subroutine YLAP2E
*
      END
 
 
 
 
 
 
 
 
 
 
 
 
 
 

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