C CLI221 SOURCE OF166741 24/12/13 21:15:30 12097 SUBROUTINE CLI221(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 : CLI221 C C DESCRIPTION : Subroutine appellée par CLIM22 C calcul de RESIDU et CLIM at the board C OPTION: 'INRI' 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 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 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(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 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 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 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 1 CONTINUE C SEGDES MELEFC C SEGSUP MLEMC SEGSUP MLEMCB SEGSUP 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 C 9999 CONTINUE RETURN END