C CLI121 SOURCE OF166741 24/12/13 21:15:12 12097 SUBROUTINE CLI121(MELEMF,MELEMC,MELECB,MELEFC,INORM,ICHPVO,ICHPSU, & IROC,IVITC,IPC,IGAMC,ICHLIM,ICHRES,ICHRLI) C************************************************************************ C C PROJET : CASTEM 2000 C C NOM : CLI121 C C DESCRIPTION : Subroutine appellée par CLIM11 C C LANGAGE : FORTRAN 77 + ESOPE 2000 (avec estensions CISI) C C AUTEUR : A. BECCANTINI, DRN/DMT/SEMT/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, MPGAMC.MPOVAL, MPLIM.MPOVAL, & MPRES.MPOVAL, MPRLI.MPOVAL INTEGER MELEMF,MELEMC,MELECB,INORM,ICHPVO,ICHPSU, IROC,IVITC,IPC & ,IGAMC,ICHLIM,ICHRES,ICHRLI,ICEL,NFAC,IFAC & ,NGF,NGC,NLF,NLC,NLCB REAL*8 VOLU,SURF,RC,PC,UXC,UYC,UZC,GAMC,CNX,CNY,CNZ,CTX,CTY,CTZ & ,CT2X,CT2Y,CT2Z,RF,PF,UXF,UYF,UZF & ,UNC,UNF,ASONC,ASONF * & ,UTF,UT2F,SF & ,UTC,UT2C,SC & ,USGM1,DSGM1,G1,G3,ASON2,S,UT,UT2,UN,RHO,P,UX,UY,UZ * & ,CACCA,EPS CHARACTER*(8) TYPE C C C**** KRIPAD pour la correspondance global/local C CALL KRIPAD(MELEMC,MLEMC) C SEGINI MLEMC CALL KRIPAD(MELECB,MLEMCB) C SEGINI MLEMCB CALL KRIPAD(MELEMF,MLEMF) C SEGINI 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**** LICHT active les MPOVALs en *MOD C C SEGACT MPNORM*MOD C SEGACT MPOVSU*MOD C SEGACT MPOVOL*MOD C 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(IGAMC,MPGAMC,TYPE,ICEL) CALL LICHT(ICHLIM,MPLIM,TYPE,ICEL) CALL LICHT(ICHRES,MPRES,TYPE,ICEL) CALL LICHT(ICHRLI,MPRLI,TYPE,ICEL) C C SEGACT *MOD C SEGACT *MOD C SEGACT *MOD C SEGACT *MOD C SEGACT *MOD C SEGACT *MOD C SEGACT *MOD C 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 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 Variables au centre RC=MPRC.VPOCHA(NLC,1) PC=MPPC.VPOCHA(NLC,1) UXC=MPVC.VPOCHA(NLC,1) UYC=MPVC.VPOCHA(NLC,2) GAMC=MPGAMC.VPOCHA(NLC,1) IF(IDIM.EQ.3)UZC=MPVC.VPOCHA(NLC,3) C Variables à la face 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) C C******* On calcule UN, UT, UT2, ASON, S C UNC=(UXC*CNX)+(UYC*CNY)+(UZC*CNZ) UTC=(UXC*CTX)+(UYC*CTY)+(UZC*CTZ) UT2C=(UXC*CT2X)+(UYC*CT2Y)+(UZC*CT2Z) 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=(GAMC*PF/RF)**0.5D0 C SC=PC/(RC**GAMC) * SF=PF/(RF**GAMC) C C******* Densite, vitesse, pression sur le bord C USGM1=1.0D0/(GAMC-1.0D0) DSGM1=2.0D0*USGM1 G1=UNF-(DSGM1*ASONF) G3=UNC+(DSGM1*ASONC) UN=0.5D0*(G1+G3) ASON2=(0.5D0*(G3-G1)) ASON2=ASON2/DSGM1 ASON2=ASON2*ASON2 S=SC UT=UTC UT2=UT2C RHO=ASON2/(GAMC*S) RHO=RHO**USGM1 P=RHO*ASON2/GAMC 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 C C******************************************************* C******* Test : we compute RHO*UN*SURF/VOLU C and its derivative with respect to RHO CC******************************************************* CC C CACCA=RHO*UN*SURF/VOLU C EPS=1.0D-6 C RC=RC*(1+EPS) CC CC******* On calcule UN, UT, UT2, ASON, S CC C UNC=(UXC*CNX)+(UYC*CNY)+(UZC*CNZ) C* UTC=(UXC*CTX)+(UYC*CTY)+(UZC*CTZ) C* UT2C=(UXC*CT2X)+(UYC*CT2Y)+(UZC*CT2Z) C UNF=(UXF*CNX)+(UYF*CNY)+(UZF*CNZ) C UTF=(UXF*CTX)+(UYF*CTY)+(UZF*CTZ) C UT2F=(UXF*CT2X)+(UYF*CT2Y)+(UZF*CT2Z) CC C ASONC=(GAMC*PC/RC)**0.5D0 C ASONF=(GAMC*PF/RF)**0.5D0 CC C* SC=PC/(RC**GAMC) C SF=PF/(RF**GAMC) CC CC******* Densite, vitesse, pression sur le bord CC C USGM1=1.0D0/(GAMC-1.0D0) C DSGM1=2.0D0*USGM1 C G1=UNC-(DSGM1*ASONC) C G3=UNF+(DSGM1*ASONF) C UN=0.5D0*(G1+G3) C ASON2=(0.5D0*(G3-G1)) C ASON2=ASON2/DSGM1 C ASON2=ASON2*ASON2 C S=SF C UT=UTF C UT2=UT2F C RHO=ASON2/(GAMC*S) C RHO=RHO**USGM1 C write(*,*) (((RHO*UN*SURF/VOLU) - CACCA)/(RC*EPS)) CC******************************************************* C*************** FIN TEST ****************************** C******************************************************* C C******* Residuum (son SPG a le meme ordre que MELEFC) C MPRES.VPOCHA(IFAC,1)=-1*RHO*UN*SURF/VOLU MPRES.VPOCHA(IFAC,2)=-1*(RHO*UN*UX+(P*CNX))*SURF/VOLU MPRES.VPOCHA(IFAC,3)=-1*(RHO*UN*UY+(P*CNY))*SURF/VOLU IF(IDIM.EQ.3) & MPRES.VPOCHA(IFAC,4)=-1*(RHO*UN*UZ+(P*CNZ))*SURF/VOLU MPRES.VPOCHA(IFAC,IDIM+2)=-1*((UN*GAMC*USGM1*P) + & (0.5D0*RHO*UN*(UN*UN+UT*UT+UT2*UT2)))*SURF/VOLU ENDDO C SEGDES MELEFC C SEGSUP MLEMC SEGSUP MLEMCB SEGSUP MLEMF C SEGDES MPNORM SEGDES MPVOL SEGDES MPSURF SEGDES MPRC SEGDES MPPC SEGDES MPVC SEGDES MPGAMC SEGDES MPLIM SEGDES MPRES SEGDES MPRLI C 9999 CONTINUE RETURN END