Sources Esope | Cast3M

cli261
C CLI261    SOURCE    OF166741  24/12/13    21:15:37     12097          
      SUBROUTINE CLI261(NSP,MELEMF,MELEMC,MELECB,MELEFC,INORM,ICHPVO,
     &     ICHPSU,LRECP,LRECV,IROC,IVITC,IPC,IYN,ICHLIM,ICHRES,ICHRLI)
C************************************************************************
C
C PROJET            :  CASTEM 2000
C
C NOM               :  CLI261
C
C DESCRIPTION       :  Subroutine appellée par CLIM22
C                       calcul de RESIDU et CLIM at the board
C                       OPTION: 'INSU' 2D
C
C LANGAGE           :  FORTRAN 77 + ESOPE 2000 (avec estensions CISI)
C
C AUTEUR            :  S.Kudriakov, DEN/DM2S/SFME/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, MPYN.MPOVAL, MPLIM.MPOVAL,
     &     MPRES.MPOVAL,  MPRLI.MPOVAL
C----------------------------------------
      INTEGER MELEMF,MELEMC,MELECB,INORM,ICHPVO,ICHPSU, IROC,IVITC,IPC
     &     ,IYN,ICHLIM,ICHRES,ICHRLI,ICEL,NFAC,IFAC
     &     ,NGF,NGC,NLF,NLC,NLCB,LRECP,LRECV,I,NSP,NESP
      REAL*8 VOLU,SURF,RC,PC,UXC,UYC,UZC,GAMC,CNX,CNY,CNZ,CTX,CTY,CTZ
     &     ,CT2X,CT2Y,CT2Z,RF,PF,UXF,UYF,UZF,TOP,BOT
     &     ,UNC,UNF,UTF,UT2F,SF,GAMF,ECIN,PSRF,HTF,GM1
     &     ,CELLT,UT2C,UTC
      CHARACTER*(8) TYPE
C------------------------------------------------------------
-INC SMLREEL
       POINTEUR MLRECP.MLREEL, MLRECV.MLREEL
C-------------------------------------------------------
C**********  Les CP's and CV's   ***********************
C-------------------------------------------------------
      SEGMENT GCONST
         REAL*8 GC(NSP)
      ENDSEGMENT
      POINTEUR CP.GCONST, CV.GCONST
C-------------------------------------------------------------
C******* Les fractionines massiques **************************
C-------------------------------------------------------------
      SEGMENT FRAMAS
         REAL*8 YET(NSP)
      ENDSEGMENT
      POINTEUR YC.FRAMAS, YF.FRAMAS
C-------------------------------------------------------------
C**********  Segments for the flux-vector  *******************
C-------------------------------------------------------------
       SEGMENT FUNEL
         REAL*8 FU(4+NSP)
       ENDSEGMENT
       POINTEUR flux2D.funel, flux3D.funel
       SEGINI FLUX2D
       SEGINI FLUX3D
C------------------------------------------------------
C**** KRIPAD pour la correspondance global/local
C------------------------------------------------------
      CALL KRIPAD(MELEMC,MLEMC)
      CALL KRIPAD(MELECB,MLEMCB)
      CALL KRIPAD(MELEMF,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**** CHPOINTs des variables
C------------------------------------------------------
      CALL LICHT(IROC,MPRC,TYPE,ICEL)
      CALL LICHT(IVITC,MPVC,TYPE,ICEL)
      CALL LICHT(IPC,MPPC,TYPE,ICEL)
      CALL LICHT(IYN,MPYN,TYPE,ICEL)
      CALL LICHT(ICHLIM,MPLIM,TYPE,ICEL)
      CALL LICHT(ICHRES,MPRES,TYPE,ICEL)
      CALL LICHT(ICHRLI,MPRLI,TYPE,ICEL)
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
      UT2F=0.0D0
      CNZ=0.0D0
      CTZ=0.0D0
      CT2X=0.0D0
      CT2Y=0.0D0
      CT2Z=0.0D0
      DO 1 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----------------------------------------------
C        In CASTEM les normales sont sortantes
C----------------------------------------------
         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----------------------------------------
       SEGINI CP, CV
       MLRECP = LRECP
       MLRECV = LRECV
       SEGACT MLRECP, MLRECV
       DO 10 I=1,(NSP-1)
           CP.GC(I)=MLRECP.PROG(I)
           CV.GC(I)=MLRECV.PROG(I)
 10    CONTINUE
       CP.GC(NSP)=MLRECP.PROG(NSP)
       CV.GC(NSP)=MLRECV.PROG(NSP)
C----------------------------
C        Variables au centre
C----------------------------
         RC=MPRC.VPOCHA(NLC,1)
         PC=MPPC.VPOCHA(NLC,1)
         UXC=MPVC.VPOCHA(NLC,1)
         UYC=MPVC.VPOCHA(NLC,2)
         IF(IDIM.EQ.3)UZC=MPVC.VPOCHA(NLC,3)
         SEGINI YC
         SEGACT MPYN
         DO 100 I=1,(NSP-1)
           YC.YET(I)=MPYN.VPOCHA(NLC,I)
 100     CONTINUE
C----------------------------
C        Variables à la face
C----------------------------
         HTF=MPLIM.VPOCHA(NLCB,1)
         SF=MPLIM.VPOCHA(NLCB,2)
         SEGINI YF
         DO 101 I=1,(NSP-1)
           YF.YET(I)=MPLIM.VPOCHA(NLCB,2+I)
 101     CONTINUE
         UTF=0.0D0
c-------------------------------------------------------------
c  Computing GAMMA at the cell-center
c-------------------------------------------------------------
         top=0.0D0
         bot=0.0D0
         do 102  i=1,(nsp-1)
          top=top+yc.yet(i)*(cp.gc(i)-cp.gc(nsp))
          bot=bot+yc.yet(i)*(cv.gc(i)-cv.gc(nsp))
 102     continue
         top=cp.gc(nsp)+top
         bot=cv.gc(nsp)+bot
         GAMC=top/bot
c-------------------------------------------------------------
c  Computing GAMMA at the face-center
c-------------------------------------------------------------
         top=0.0D0
         bot=0.0D0
         do 103  i=1,(nsp-1)
          top=top+yf.yet(i)*(cp.gc(i)-cp.gc(nsp))
          bot=bot+yf.yet(i)*(cv.gc(i)-cv.gc(nsp))
 103     continue
         top=cp.gc(nsp)+top
         bot=cv.gc(nsp)+bot
         GAMF=top/bot
         GM1=GAMF-1.0D0
C---------------------------------------
C******* On calcule UN, UT, UT2, ASON, S
C---------------------------------------
         UNC=(UXC*CNX)+(UYC*CNY)+(UZC*CNZ)
         UNF=UNC
         UTC=(UXC*CTX)+(UYC*CTY)+(UZC*CTZ)
C----------------------------------
         UXF=UNF*CNX+UTF*CTX+UT2F*CT2X
         UYF=UNF*CNY+UTF*CTY+UT2F*CT2Y
         UZF=UNF*CNZ+UTF*CTZ+UT2F*CT2Z
C----------------------------------
         ECIN=0.5D0*((UXF*UXF)+(UYF*UYF)+(UZF*UZF))
         PSRF=(GM1/GAMF)*(HTF-ECIN)
         RF=PSRF/SF
         RF=RF**(1.0D0/GM1)
         PF=SF*(RF**GAMF)
C-----------------------------------------------
C******* Densite, vitesse, pression sur le bord
C-----------------------------------------------
         MPRLI.VPOCHA(NLCB,1)=RF
         MPRLI.VPOCHA(NLCB,2)=UXF
         MPRLI.VPOCHA(NLCB,3)=UYF
         IF(IDIM.EQ.3) MPRLI.VPOCHA(NLCB,4)=UZF
         MPRLI.VPOCHA(NLCB,IDIM+2)=PF
         do 104  i=1,(nsp-1)
           MPRLI.VPOCHA(NLCB,IDIM+2+I)=YF.YET(I)
 104     continue
C---------------------------------------------------
C******* Probleme de Riemann entre l'etat gauche
C        RF,UNC,UTF,UT2F,PF et l'etat droite
C        RC,UNC,UTC,UT2C,PC
C        On utilise AUSM+
C        Flux dans le repaire normale
C---------------------------------------------------
         NESP=NSP-1
         IF(IDIM.EQ.2)THEN
            CALL FAUSMP(NESP,
     &           GAMF,RF,PF,UNC,UTF,
     &           GAMC,RC,PC,UNC,UTC,
     &           YF.YET,YC.YET,
     &           FLUX2D.FU,
     &           CELLT)
C-------------------------------------------------------
C******* Residuum (son SPG a le meme ordre que MELEFC)
C-------------------------------------------------------
            MPRES.VPOCHA(IFAC,1)=FLUX2D.FU(1)*SURF/VOLU
         MPRES.VPOCHA(IFAC,2)=((FLUX2D.FU(2)*CNX)+(FLUX2D.FU(3)*CTX))
     &           *SURF/VOLU
         MPRES.VPOCHA(IFAC,3)=((FLUX2D.FU(2)*CNY)+(FLUX2D.FU(3)*CTY))
     &           *SURF/VOLU
            MPRES.VPOCHA(IFAC,4)=FLUX2D.FU(4)*SURF/VOLU
            do 105  i=1,(nsp-1)
              MPRES.VPOCHA(IFAC,4+I)=FLUX2D.FU(4+I)*SURF/VOLU
 105       continue
         ELSE
            CALL FAUSM3(NESP,
     &           GAMF,RF,PF,UNC,UTF,UT2F,
     &           GAMC,RC,PC,UNC,UTC,UT2C,
     &           YF.YET,YC.YET,
     &           FLUX3D.FU,
     &           CELLT)
C------------------------------------------------------
C******* Residuum (son SPG a le meme ordre que MELEFC)
C------------------------------------------------------
            MPRES.VPOCHA(IFAC,1)=FLUX3D.FU(1)*SURF/VOLU
         MPRES.VPOCHA(IFAC,2)=((FLUX3D.FU(2)*CNX)+(FLUX3D.FU(3)*CTX)+
     &           (FLUX3D.FU(4)*CT2X))*SURF/VOLU
         MPRES.VPOCHA(IFAC,3)=((FLUX3D.FU(2)*CNY)+(FLUX3D.FU(3)*CTY)+
     &           (FLUX3D.FU(4)*CT2Z))*SURF/VOLU
         MPRES.VPOCHA(IFAC,4)=((FLUX3D.FU(2)*CNZ)+(FLUX3D.FU(3)*CTZ)+
     &           (FLUX3D.FU(4)*CT2Z))*SURF/VOLU
            MPRES.VPOCHA(IFAC,5)=FLUX3D.FU(5)*SURF/VOLU
            do 106  i=1,(nsp-1)
              MPRES.VPOCHA(IFAC,5+I)=FLUX3D.FU(5+I)*SURF/VOLU
 106        continue
         ENDIF
 1      CONTINUE
C
      SEGDES MELEFC
C
c      SEGSUP MLEMC
c      SEGSUP MLEMCB
c      SEGSUP MLEMF
c-------------------------
      SEGDES MLEMC
      SEGDES MLEMCB
      SEGDES MLEMF
C
      SEGDES MPNORM
      SEGDES MPVOL
      SEGDES MPSURF
      SEGDES MPRC
      SEGDES MPPC
      SEGDES MPVC
      SEGDES MPYN
      SEGDES MPLIM
      SEGDES MPRES
      SEGDES MPRLI
      SEGDES MLRECP
      SEGDES MLRECV
      SEGDES YC
      SEGDES YF
      SEGDES FLUX2D
      SEGDES FLUX3D
C
 9999 CONTINUE
      RETURN
      END