C CLI281 SOURCE CB215821 20/11/25 13:20:56 10792 SUBROUTINE CLI281(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 : CLI281 C C DESCRIPTION : Subroutine appellée par CLIM22 C calcul de RESIDU et CLIM at the board C OPTION: 'INOU' 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) IMPLICIT REAL*8(A-H,O-Z) -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 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,TEMP 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) c write(*,*) 'hell', NFAC UZC=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=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) 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 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 We decide whether the boundary conditions are 'inlet' or C 'outlet' by looking at the UNC (UNC > 0.0 -- 'outlet') C------------------------------------------------------------- UNC=(UXC*CNX)+(UYC*CNY)+(UZC*CNZ) c write(*,*) 'before', UNC, NLC,NGC IF(UNC .GE. 0.0D0) THEN C-------------------------------------- C 'OUTP' -- outlet with given pressure C-------------------------------------- C----------------------------------------- C Variables à la face C----------------------------------------- PF=MPLIM.VPOCHA(NLCB,NSP+2) C--------------------------------------- C******* On calcule UN, UT, UT2 C--------------------------------------- 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 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 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 ENDIF TEMP=-1000.0 ELSE c write(*,*) 'insuuuuuuuu', NLCB C------------------------------------------------------- C 'INSU' -- 'inlet' b.c. with given Ht and S (entropy) C------------------------------------------------------- UZF=0.0D0 UT2F=0.0D0 C---------------------------- C Variables à la face C---------------------------- HTF=MPLIM.VPOCHA(NLCB,1) SF=MPLIM.VPOCHA(NLCB,2) SEGINI YF TEMP=1000.0 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 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 changing direction of the normal vector C----------------------------------------- CNX=-1.0D0*CNX CNY=-1.0D0*CNY IF(IDIM.EQ.2)THEN CTX=-1.0D0*CNY CTY=CNX ELSE CNZ=-1.0D0*CNZ CTX=-1.0D0*CTX CTY=-1.0D0*CTY CTZ=-1.0D0*CTZ CT2X=-1.0D0*CT2X CT2Y=-1.0D0*CT2Y CT2Z=-1.0D0*CT2Z ENDIF 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 1040 i=1,(nsp-1) MPRLI.VPOCHA(NLCB,IDIM+2+I)=YF.YET(I) 1040 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 1050 i=1,(nsp-1) MPRES.VPOCHA(IFAC,4+I)=FLUX2D.FU(4+I)*SURF/VOLU 1050 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 1060 i=1,(nsp-1) MPRES.VPOCHA(IFAC,5+I)=FLUX3D.FU(5+I)*SURF/VOLU 1060 continue ENDIF ENDIF 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 MPLIM SEGDES MPRES SEGDES MPRLI SEGDES MLRECP SEGDES MLRECV SEGDES YC IF(TEMP .GT. 0.0D0) SEGDES YF SEGDES FLUX2D SEGDES FLUX3D C RETURN END