Sources Esope | Cast3M

cl251t
C CL251T    SOURCE    OF166741  24/12/13    21:15:07     12097          
      SUBROUTINE CL251T(NSP,MELEMF,MELEMC,MELECB,MELEFC,INORM,ICHPVO,
     &    ICHPSU,LRECP,LRECV,IROC,IVITC,IPC,IYC,
     &    IKAN,IEPSN,IK0N,ICHLIM,ICHRES,ICHRLI)
C************************************************************************
C
C PROJET            :  CASTEM 2000
C
C NOM               :  CLI251T
C
C DESCRIPTION       :  Subroutine appellée par CLIM22
C                      Outlet B.C. (known pressure)
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
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
 
-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,
     &     MPK0C.MPOVAL
C
      INTEGER INORM,ICHPVO,ICHPSU, IROC,IVITC,IPC
     &     ,IYC,ICHLIM,ICHRES,ICHRLI,ICEL,NFAC,IFAC,IKAN, IEPSN
     &     ,IK0N,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,K0C
      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)
      IF (IK0N .GT. 0) THEN
        CALL LICHT(IK0N,MPK0C,TYPE,ICEL)
      ENDIF
      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)
         IF (IK0N .GT. 0) THEN
           K0C=MPK0C.VPOCHA(NLC,1)
         ENDIF
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=MPLIM.VPOCHA(NLCB,1)
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
         IF (IK0N .GT. 0) THEN
           MPRLI.VPOCHA(NLCB,IDIM+NSP+4)=K0C
         ENDIF
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
             IF (IK0N .GT. 0) THEN
               MPRES.VPOCHA(IFAC,6+NSP)=-1*K0C*FLUX2D.FU(1)*SURF/VOLU
             ENDIF
         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)*CT2Y))*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
            IF (IK0N .GT. 0) THEN
              MPRES.VPOCHA(IFAC,7+NSP)=-1*K0C*FLUX3D.FU(1)*SURF/VOLU
            ENDIF
         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
      IF (IK0N .GT. 0) THEN
        SEGDES MPK0C
      ENDIF
      SEGDES MPLIM
      SEGDES MPRES
      SEGDES MPRLI
      SEGDES YC
      SEGDES FLUX2D
      SEGDES FLUX3D
C
 9999 CONTINUE
      RETURN
      END