Sources Esope | Cast3M

cli111
C CLI111    SOURCE    OF166741  24/12/13    21:15:10     12097          
      SUBROUTINE CLI111(MELEMF,MELEMC,MELECB,MELEFC,INORM,ICHPVO,ICHPSU,
     &     IROC,IVITC,IPC,IGAMC,ICHLIM,ICHRES,ICHRLI)
C************************************************************************
C
C PROJET            :  CASTEM 2000
C
C NOM               :  CLI111
C
C DESCRIPTION       :  Subroutine appellée par CLIM11
C
C LANGAGE           :  FORTRAN 77 + ESOPE 2000 (avec estensions CISI)
C
C AUTEUR            :  A. BECCANTINI, DRN/DMT/SEMT/LTMF
C
C************************************************************************
C
C APPELES (Calcul)  :
C
C************************************************************************
C
C HISTORIQUE (Anomalies et modifications éventuelles)
C
C HISTORIQUE :
C
C************************************************************************
C
      IMPLICIT INTEGER(I-N)
-INC PPARAM
-INC CCOPTIO
-INC SMLMOTS
-INC SMELEME
      POINTEUR MELEFC.MELEME
-INC SMLENTI
      POINTEUR MLEMC.MLENTI, MLEMCB.MLENTI,MLEMF.MLENTI
-INC SMCHPOI
      POINTEUR MPNORM.MPOVAL, MPVOL.MPOVAL, MPSURF.MPOVAL, MPRC.MPOVAL,
     &     MPVC.MPOVAL, MPPC.MPOVAL, MPGAMC.MPOVAL, MPLIM.MPOVAL,
     &     MPRES.MPOVAL,  MPRLI.MPOVAL
 
      INTEGER MELEMF,MELEMC,MELECB,INORM,ICHPVO,ICHPSU, IROC,IVITC,IPC
     &     ,IGAMC,ICHLIM,ICHRES,ICHRLI,ICEL,NFAC,IFAC
     &     ,NGF,NGC,NLF,NLC,NLCB
      REAL*8 VOLU,SURF,RC,PC,UXC,UYC,UZC,GAMC,CNX,CNY,CNZ,CTX,CTY,CTZ
     &     ,CT2X,CT2Y,CT2Z,RF,PF,UXF,UYF,UZF
     &     ,UNC,UNF,UTF,UT2F,SF,ASONC,ASONF
*     &     ,UTC,UT2C,SC
     &     ,USGM1,DSGM1,G1,G3,ASON2,S,UT,UT2,UN,RHO,P,UX,UY,UZ
*     &     ,CACCA,EPS
      CHARACTER*(8) TYPE
C
C
C**** KRIPAD pour la correspondance global/local
C
      CALL KRIPAD(MELEMC,MLEMC)
C     SEGINI MLEMC
      CALL KRIPAD(MELECB,MLEMCB)
C     SEGINI MLEMCB
      CALL KRIPAD(MELEMF,MLEMF)
C     SEGINI MLEMF
C
C**** CHPOINTs de la table DOMAINE
C
      CALL LICHT(INORM,MPNORM,TYPE,ICEL)
      CALL LICHT(ICHPVO,MPVOL,TYPE,ICEL)
      CALL LICHT(ICHPSU,MPSURF,TYPE,ICEL)
C
C**** LICHT active les MPOVALs en *MOD
C
C     SEGACT MPNORM*MOD
C     SEGACT MPOVSU*MOD
C     SEGACT MPOVOL*MOD
C
C
C**** CHPOINTs des variables
C
      CALL LICHT(IROC,MPRC,TYPE,ICEL)
      CALL LICHT(IVITC,MPVC,TYPE,ICEL)
      CALL LICHT(IPC,MPPC,TYPE,ICEL)
      CALL LICHT(IGAMC,MPGAMC,TYPE,ICEL)
      CALL LICHT(ICHLIM,MPLIM,TYPE,ICEL)
      CALL LICHT(ICHRES,MPRES,TYPE,ICEL)
      CALL LICHT(ICHRLI,MPRLI,TYPE,ICEL)
C
C     SEGACT *MOD
C     SEGACT *MOD
C     SEGACT *MOD
C     SEGACT *MOD
C     SEGACT *MOD
C     SEGACT *MOD
C     SEGACT *MOD
C
C
C**** Boucle sur le face pour le calcul des invariants de
C     Riemann et du flux
C
      SEGACT MELEFC
      NFAC=MELEFC.NUM(/2)
      UZC=0.0D0
      UZF=0.0D0
      CNZ=0.0D0
      CTZ=0.0D0
      CT2X=0.0D0
      CT2Y=0.0D0
      CT2Z=0.0D0
      DO IFAC=1,NFAC,1
         NGF=MELEFC.NUM(1,IFAC)
         NGC=MELEFC.NUM(2,IFAC)
         NLF=MLEMF.LECT(NGF)
         NLC=MLEMC.LECT(NGC)
         NLCB=MLEMCB.LECT(NGF)
         VOLU=MPVOL.VPOCHA(NLC,1)
         SURF=MPSURF.VPOCHA(NLF,1)
C        In CASTEM les normales sont sortantes
         CNX=-1*MPNORM.VPOCHA(NLF,1)
         CNY=-1*MPNORM.VPOCHA(NLF,2)
         IF(IDIM.EQ.2)THEN
            CTX=-1.0D0*CNY
            CTY=CNX
         ELSE
            CNZ=-1*MPNORM.VPOCHA(NLF,3)
            CTX=-1*MPNORM.VPOCHA(NLF,4)
            CTY=-1*MPNORM.VPOCHA(NLF,5)
            CTZ=-1*MPNORM.VPOCHA(NLF,6)
            CT2X=-1*MPNORM.VPOCHA(NLF,7)
            CT2Y=-1*MPNORM.VPOCHA(NLF,8)
            CT2Z=-1*MPNORM.VPOCHA(NLF,9)
         ENDIF
C        Variables au centre
         RC=MPRC.VPOCHA(NLC,1)
         PC=MPPC.VPOCHA(NLC,1)
         UXC=MPVC.VPOCHA(NLC,1)
         UYC=MPVC.VPOCHA(NLC,2)
         GAMC=MPGAMC.VPOCHA(NLC,1)
         IF(IDIM.EQ.3)UZC=MPVC.VPOCHA(NLC,3)
C        Variables à la face
         RF=MPLIM.VPOCHA(NLCB,1)
         UXF=MPLIM.VPOCHA(NLCB,2)
         UYF=MPLIM.VPOCHA(NLCB,3)
         IF(IDIM.EQ.3)UZF=MPLIM.VPOCHA(NLCB,4)
         PF=MPLIM.VPOCHA(NLCB,IDIM+2)
C
C******* On calcule UN, UT, UT2, ASON, S
C
         UNC=(UXC*CNX)+(UYC*CNY)+(UZC*CNZ)
*         UTC=(UXC*CTX)+(UYC*CTY)+(UZC*CTZ)
*         UT2C=(UXC*CT2X)+(UYC*CT2Y)+(UZC*CT2Z)
         UNF=(UXF*CNX)+(UYF*CNY)+(UZF*CNZ)
         UTF=(UXF*CTX)+(UYF*CTY)+(UZF*CTZ)
         UT2F=(UXF*CT2X)+(UYF*CT2Y)+(UZF*CT2Z)
C
         ASONC=(GAMC*PC/RC)**0.5D0
         ASONF=(GAMC*PF/RF)**0.5D0
C
*         SC=PC/(RC**GAMC)
         SF=PF/(RF**GAMC)
C
C******* Densite, vitesse, pression sur le bord
C
         USGM1=1.0D0/(GAMC-1.0D0)
         DSGM1=2.0D0*USGM1
         G1=UNC-(DSGM1*ASONC)
         G3=UNF+(DSGM1*ASONF)
         UN=0.5D0*(G1+G3)
         ASON2=(0.5D0*(G3-G1))
         ASON2=ASON2/DSGM1
         ASON2=ASON2*ASON2
         S=SF
         UT=UTF
         UT2=UT2F
         RHO=ASON2/(GAMC*S)
         RHO=RHO**USGM1
         P=RHO*ASON2/GAMC
         UX=(UN*CNX)+(UT*CTX)+(UT2*CT2X)
         UY=(UN*CNY)+(UT*CTY)+(UT2*CT2Y)
         UZ=(UN*CNZ)+(UT*CTZ)+(UT2*CT2Z)
C
         MPRLI.VPOCHA(NLCB,1)=RHO
         MPRLI.VPOCHA(NLCB,2)=UX
         MPRLI.VPOCHA(NLCB,3)=UY
         IF(IDIM.EQ.3) MPRLI.VPOCHA(NLCB,4)=UZ
         MPRLI.VPOCHA(NLCB,IDIM+2)=P
C
C*******************************************************
C******* Test : we compute RHO*UN*SURF/VOLU
C               and its derivative with respect to RHO
CC*******************************************************
CC
C         CACCA=RHO*UN*SURF/VOLU
C         EPS=1.0D-6
C         RC=RC*(1+EPS)
CC
CC******* On calcule UN, UT, UT2, ASON, S
CC
C         UNC=(UXC*CNX)+(UYC*CNY)+(UZC*CNZ)
C*         UTC=(UXC*CTX)+(UYC*CTY)+(UZC*CTZ)
C*         UT2C=(UXC*CT2X)+(UYC*CT2Y)+(UZC*CT2Z)
C         UNF=(UXF*CNX)+(UYF*CNY)+(UZF*CNZ)
C         UTF=(UXF*CTX)+(UYF*CTY)+(UZF*CTZ)
C         UT2F=(UXF*CT2X)+(UYF*CT2Y)+(UZF*CT2Z)
CC
C         ASONC=(GAMC*PC/RC)**0.5D0
C         ASONF=(GAMC*PF/RF)**0.5D0
CC
C*         SC=PC/(RC**GAMC)
C         SF=PF/(RF**GAMC)
CC
CC******* Densite, vitesse, pression sur le bord
CC
C         USGM1=1.0D0/(GAMC-1.0D0)
C         DSGM1=2.0D0*USGM1
C         G1=UNC-(DSGM1*ASONC)
C         G3=UNF+(DSGM1*ASONF)
C         UN=0.5D0*(G1+G3)
C         ASON2=(0.5D0*(G3-G1))
C         ASON2=ASON2/DSGM1
C         ASON2=ASON2*ASON2
C         S=SF
C         UT=UTF
C         UT2=UT2F
C         RHO=ASON2/(GAMC*S)
C         RHO=RHO**USGM1
C         write(*,*) (((RHO*UN*SURF/VOLU) - CACCA)/(RC*EPS))
CC*******************************************************
C*************** FIN TEST ******************************
C*******************************************************
C
C******* Residuum (son SPG a le meme ordre que MELEFC)
C
         MPRES.VPOCHA(IFAC,1)=RHO*UN*SURF/VOLU
         MPRES.VPOCHA(IFAC,2)=(RHO*UN*UX+(P*CNX))*SURF/VOLU
         MPRES.VPOCHA(IFAC,3)=(RHO*UN*UY+(P*CNY))*SURF/VOLU
         IF(IDIM.EQ.3)
     &        MPRES.VPOCHA(IFAC,4)=(RHO*UN*UZ+(P*CNZ))*SURF/VOLU
         MPRES.VPOCHA(IFAC,IDIM+2)=((UN*GAMC*USGM1*P) +
     &        (0.5D0*RHO*UN*(UN*UN+UT*UT+UT2*UT2)))*SURF/VOLU
      ENDDO
C
      SEGDES MELEFC
C
      SEGSUP MLEMC
      SEGSUP MLEMCB
      SEGSUP MLEMF
C
      SEGDES MPNORM
      SEGDES MPVOL
      SEGDES MPSURF
      SEGDES MPRC
      SEGDES MPPC
      SEGDES MPVC
      SEGDES MPGAMC
      SEGDES MPLIM
      SEGDES MPRES
      SEGDES MPRLI
C
 9999 CONTINUE
      RETURN
      END