Sources Esope | Cast3M

cl291t
C CL291T    SOURCE    OF166741  24/12/13    21:15:09     12097          
      SUBROUTINE CL291T(NSP,MELEMF,MELEMC,MELECB,MELEFC,INORM,ICHPVO,
     &    ICHPSU,LRECP,LRECV,IROC,IVITC,IPC,IYC,
     &    IKAN,IEPSN,ICHLIM,ICHRES,ICHRLI)
C************************************************************************
C
C PROJET            :  CASTEM 2000
C
C NOM               :  CLI291T
C
C DESCRIPTION       :  Subroutine appellée par CLIM22
C                      Supersonic Outlet B.C.
C
C LANGAGE           :  FORTRAN 77 + ESOPE 2000 (avec estensions CISI)
C
C AUTEUR            :  S. Kudriakov, DEN/DM2S/SFME/LTMF
C
C************************************************************************
C ENTREES :   NSP (type ENTIER) : full number of species
C                                 ('ESPEULE' + 1) ;
C             MELEMF (type MELEME) : maillage des faces des
C                                             éléments.
C             MELEMC (type MELEME) : maillage des centres des
C                                             éléments.
C             MELECB (type MELEME) : maillage des centres des éléments
C                                    de la frontière
C             MELEFC (type MELEME) : connectivités face-(centre
C                                             gauche, centre droit).
C             INORM  (type MCHPOI) : normales at the faces
C             ICHPVO (type MCHPOI) : volumes of the cells
C             ICHPSU (type MCHPOI) : surfaces of the cell-interfaces
C             LRECP  (type MLREEL) : list of CP's of the species
C             LRECV  (type MLREEL) : list of CV's of the species
C             IROC   (type MCHPOI) : densityes at the centres of the domain
C             IVITC  (type MCHPOI) : velocities at the centres of the domain
C             IPC    (type MCHPOI) : pressure at the centres of the domain
C             IYN    (type MCHPOI) : mass fraction of the species at the
C                                    centres of the domain
C             IKAN   (type MCHPOI) : turbulent kinetic energy, k,
C                                    at the centres of the domain
C             IEPSN  (type MCHPOI) : rate of dissipated turb. energy
C                                    at the centres of the domain
C             ICHLIM (type MCHPOI) : boundary conditions at the centers
C                                    of the boundary (CHPVID here)
C-----------------------------------------------------
C SORTIES:    ICHRES (type MCHPOI) : the contribution to the residuum
C                                    due to the boundary conditions
C                                    given at the centres of the cells
C                                    next to the boundary
C             ICHRLI (type MCHPOI) : the values at the boundary faces
C                                    found by the procedure.
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,MELEMF.MELEME,MELEMC.MELEME,MELECB.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, MPYC.MPOVAL, MPLIM.MPOVAL,
     &     MPRES.MPOVAL,  MPRLI.MPOVAL,MPKAC.MPOVAL,MPEPSC.MPOVAL
C
      INTEGER INORM,ICHPVO,ICHPSU, IROC,IVITC,IPC
     &     ,IYC,ICHLIM,ICHRES,ICHRLI,ICEL,NFAC,IFAC,IKAN, IEPSN
     &     ,NGF,NGC,NLF,NLC,NLCB,LRECP,LRECV,NSP,NESP,I
      REAL*8 VOLU,SURF,GAMC,CNX,CNY,CNZ,CTX,CTY,CTZ
     &     ,CT2X,CT2Y,CT2Z,RC,PC,UXC,UYC,UZC,PF,TOP,BOT
     &     ,UNC,UTC,UT2C,CELLT,EPSC,KAC
      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
C-------------------------------------------------------------
C**********  Segments for the flux-vector  *******************
C-------------------------------------------------------------
       SEGMENT FUNEL
         REAL*8 FU(IDIM+1+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(IYC,MPYC,TYPE,ICEL)
      CALL LICHT(IKAN,MPKAC,TYPE,ICEL)
      CALL LICHT(IEPSN,MPEPSC,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
      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=MPNORM.VPOCHA(NLF,1)
         CNY=MPNORM.VPOCHA(NLF,2)
         IF(IDIM.EQ.2)THEN
            CTX=-1.0D0*CNY
            CTY=CNX
         ELSE
            CNZ=MPNORM.VPOCHA(NLF,3)
            CTX=MPNORM.VPOCHA(NLF,4)
            CTY=MPNORM.VPOCHA(NLF,5)
            CTZ=MPNORM.VPOCHA(NLF,6)
            CT2X=MPNORM.VPOCHA(NLF,7)
            CT2Y=MPNORM.VPOCHA(NLF,8)
            CT2Z=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)
         UXC=MPVC.VPOCHA(NLC,1)
         UYC=MPVC.VPOCHA(NLC,2)
         IF(IDIM.EQ.3)UZC=MPVC.VPOCHA(NLC,3)
         PC=MPPC.VPOCHA(NLC,1)
         SEGINI YC
         SEGACT MPYC
         DO 101 I=1,(NSP-1)
           YC.YET(I)=MPYC.VPOCHA(NLC,I)
 101     CONTINUE
         KAC=MPKAC.VPOCHA(NLC,1)
         EPSC=MPEPSC.VPOCHA(NLC,1)
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        Variables à la face
C-----------------------------------------
         PF = PC
C---------------------------------------
C******* On calcule UN, UT, UT2
C---------------------------------------
         UNC=(UXC*CNX)+(UYC*CNY)+(UZC*CNZ)
         UTC=(UXC*CTX)+(UYC*CTY)+(UZC*CTZ)
         UT2C=(UXC*CT2X)+(UYC*CT2Y)+(UZC*CT2Z)
C-----------------------------------------------
C******* Densite, vitesse, pression sur le bord
C-----------------------------------------------
         MPRLI.VPOCHA(NLCB,1)=RC
         MPRLI.VPOCHA(NLCB,2)=UXC
         MPRLI.VPOCHA(NLCB,3)=UYC
         IF(IDIM.EQ.3) MPRLI.VPOCHA(NLCB,4)=UZC
         MPRLI.VPOCHA(NLCB,IDIM+2)=PF
         do 104  i=1,(nsp-1)
           MPRLI.VPOCHA(NLCB,IDIM+2+I)=YC.YET(I)
 104     continue
         MPRLI.VPOCHA(NLCB,IDIM+NSP+2)=KAC
         MPRLI.VPOCHA(NLCB,IDIM+NSP+3)=EPSC
C---------------------------------------------------
C******* Probleme de Riemann entre l'etat gauche
C        RC,UNC,UTC,UT2C,PC et l'etat droite
C        RC,UNC,UTC,UT2C,PF
C        On utilise AUSM+
C        Flux dans le repaire normale
C---------------------------------------------------
         NESP=NSP-1
         IF(IDIM.EQ.2)THEN
            CALL FAUSMP(NESP,
     &           GAMC,RC,PC,UNC,UTC,
     &           GAMC,RC,PF,UNC,UTC,
     &           YC.YET,YC.YET,
     &           FLUX2D.FU,
     &           CELLT)
C-------------------------------------------------------
C******* Residuum (son SPG a le meme ordre que MELEFC)
C-------------------------------------------------------
            MPRES.VPOCHA(IFAC,1)=-1*FLUX2D.FU(1)*SURF/VOLU
             MPRES.VPOCHA(IFAC,2)=-1*((FLUX2D.FU(2)*CNX)+
     &          (FLUX2D.FU(3)*CTX))*SURF/VOLU
             MPRES.VPOCHA(IFAC,3)=-1*((FLUX2D.FU(2)*CNY)+
     &          (FLUX2D.FU(3)*CTY))*SURF/VOLU
             MPRES.VPOCHA(IFAC,4)=-1*FLUX2D.FU(4)*SURF/VOLU
             do 105  i=1,(nsp-1)
              MPRES.VPOCHA(IFAC,4+I)=-1*FLUX2D.FU(4+I)*SURF/VOLU
 105       continue
             MPRES.VPOCHA(IFAC,4+NSP)=-1*KAC*FLUX2D.FU(1)*SURF/VOLU
             MPRES.VPOCHA(IFAC,5+NSP)=-1*EPSC*FLUX2D.FU(1)*SURF/VOLU
         ELSE
            CALL FAUSM3(NESP,
     &           GAMC,RC,PC,UNC,UTC,UT2C,
     &           GAMC,RC,PF,UNC,UTC,UT2C,
     &           YC.YET,YC.YET,
     &           FLUX3D.FU,
     &           CELLT)
C------------------------------------------------------
C******* Residuum (son SPG a le meme ordre que MELEFC)
C------------------------------------------------------
            MPRES.VPOCHA(IFAC,1)=-1*FLUX3D.FU(1)*SURF/VOLU
            MPRES.VPOCHA(IFAC,2)=-1*((FLUX3D.FU(2)*CNX)+
     &       (FLUX3D.FU(3)*CTX)+(FLUX3D.FU(4)*CT2X))*SURF/VOLU
            MPRES.VPOCHA(IFAC,3)=-1*((FLUX3D.FU(2)*CNY)+
     &       (FLUX3D.FU(3)*CTY)+(FLUX3D.FU(4)*CT2Z))*SURF/VOLU
            MPRES.VPOCHA(IFAC,4)=-1*((FLUX3D.FU(2)*CNZ)+
     &       (FLUX3D.FU(3)*CTZ)+(FLUX3D.FU(4)*CT2Z))*SURF/VOLU
            MPRES.VPOCHA(IFAC,5)=-1*FLUX3D.FU(5)*SURF/VOLU
            do 106  i=1,(nsp-1)
              MPRES.VPOCHA(IFAC,5+I)=-1*FLUX3D.FU(5+I)*SURF/VOLU
 106        continue
            MPRES.VPOCHA(IFAC,5+NSP)=-1*KAC*FLUX3D.FU(1)*SURF/VOLU
            MPRES.VPOCHA(IFAC,6+NSP)=-1*EPSC*FLUX3D.FU(1)*SURF/VOLU
         ENDIF
      ENDDO
C
      SEGDES MELEFC
C
      SEGDES MLEMC
      SEGDES MLEMCB
      SEGDES MLEMF
C
      SEGDES MPNORM
      SEGDES MPVOL
      SEGDES MPSURF
      SEGDES MPRC
      SEGDES MPPC
      SEGDES MPVC
      SEGDES MPYC
      SEGDES MPKAC
      SEGDES MPEPSC
c      SEGDES MPLIM
      SEGDES MPRES
      SEGDES MPRLI
      SEGDES YC
      SEGDES FLUX2D
      SEGDES FLUX3D
C
 9999 CONTINUE
      RETURN
      END