Sources Esope | Cast3M

cl221t
C CL221T    SOURCE    OF166741  24/12/13    21:15:06     12097          
      SUBROUTINE CL221T(NSP,MELEMF,MELEMC,MELECB,MELEFC,INORM,ICHPVO,
     &     ICHPSU,LRECP,LRECV,IROC,IVITC,IPC,IYN,
     &     IKAN,IEPSN,ICHLIM,ICHRES,ICHRLI)
C************************************************************************
C
C PROJET            :  CASTEM 2000
C
C NOM               :  CL221T
C
C DESCRIPTION       :  Subroutine appellée par CLIM22
C                       calcul de RESIDU et CLIM at the board
C                       OPTION: 'INRI' 2D with turbulence (k-\eps)
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 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, MPYN.MPOVAL, MPLIM.MPOVAL,
     &     MPRES.MPOVAL,  MPRLI.MPOVAL, MPKAC.MPOVAL,MPEPSC.MPOVAL
C----------------------------------------
      INTEGER INORM,ICHPVO,ICHPSU, IROC,IVITC,IPC
     &     ,IYN,ICHLIM,ICHRES,ICHRLI,ICEL,NFAC,IFAC,IKAN,IEPSN
     &     ,NGF,NGC,NLF,NLC,NLCB,LRECP,LRECV,I,NSP
      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,ASONC,ASONF,GAMF
     &     ,G1,G3,ASON2,S,UT,UT2,UN,RHO,P,UX,UY,UZ
     &     ,KAC,EPSC,KAF,EPSF
      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**** 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(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
      UZF=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
         KAC=MPKAC.VPOCHA(NLC,1)
         EPSC=MPEPSC.VPOCHA(NLC,1)
C----------------------------
C        Variables à la face
C----------------------------
         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)
         SEGINI YF
         DO 101 I=1,(NSP-1)
           YF.YET(I)=MPLIM.VPOCHA(NLCB,IDIM+2+I)
 101     CONTINUE
         KAF=MPLIM.VPOCHA(NLCB,IDIM+NSP+2)
         EPSF=MPLIM.VPOCHA(NLCB,IDIM+NSP+3)
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
C---------------------------------------
C******* On calcule UN, UT, UT2, ASON, S
C---------------------------------------
         UNC=(UXC*CNX)+(UYC*CNY)+(UZC*CNZ)
         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=(GAMF*PF/RF)**0.5D0
C
         SF=PF/(RF**GAMF)
C-----------------------------------------------
C******* Densite, vitesse, pression sur le bord
C-----------------------------------------------
         G1=UNC-(2.0D0*ASONC)/(GAMC-1.0D0)
         G3=UNF+(2.0D0*ASONF)/(GAMF-1.0D0)
         UN=0.5D0*(G1+G3)
         ASON2=(0.5D0*(G3-G1))
         ASON2=ASON2*(GAMF-1.0D0)/2.0D0
         ASON2=ASON2*ASON2
         S=SF
         UT=UTF
         UT2=UT2F
         RHO=ASON2/(GAMF*S)
         RHO=RHO**(1.0D0/(GAMF-1.0D0))
         P=RHO*ASON2/GAMF
         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
         do 104  i=1,(nsp-1)
           MPRLI.VPOCHA(NLCB,IDIM+2+I)=YF.YET(I)
 104     continue
         MPRLI.VPOCHA(NLCB,IDIM+NSP+2)=KAF
         MPRLI.VPOCHA(NLCB,IDIM+NSP+3)=EPSF
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*GAMF*P/(GAMF-1.0D0)) +
     &        (0.5D0*RHO*UN*(UN*UN+UT*UT+UT2*UT2)))*SURF/VOLU
         do 105  i=1,(nsp-1)
           MPRES.VPOCHA(IFAC,IDIM+2+I)=RHO*YF.YET(I)*UN*SURF/VOLU
 105     continue
         MPRES.VPOCHA(IFAC,IDIM+NSP+2)=RHO*KAF*UN*SURF/VOLU
         MPRES.VPOCHA(IFAC,IDIM+NSP+3)=RHO*EPSF*UN*SURF/VOLU
 1      CONTINUE
C
      SEGDES MELEFC
C
      SEGDES MLEMC
      SEGDES MLEMCB
      SEGDES MLEMF
C
      SEGDES MPNORM
      SEGDES MPVOL
      SEGDES MPSURF
      SEGDES MPRC
      SEGDES MPPC
      SEGDES MPVC
      SEGDES MPYN
      SEGDES MPKAC
      SEGDES MPEPSC
      SEGDES MPLIM
      SEGDES MPRES
      SEGDES MPRLI
      SEGDES MLRECP
      SEGDES MLRECV
      SEGDES YC
      SEGDES YF
C
 9999 CONTINUE
      RETURN
      END