C CLI272 SOURCE OF166741 24/12/13 21:15:41 12097 SUBROUTINE CLI272(NSP,MELEMF,MELEMC,MELECB,MELEFC,MELRES,INORM, & ICHPVO,ICHPSU,LRECP,LRECV, & IROC,IVITC,IPC,IYC,ICHLIM,ILIINC,ILIINP,IJAC,IJACO) C************************************************************************ C C PROJET : CASTEM 2000 C C NOM : CLI112 C C DESCRIPTION : Subroutine appellée par CLIM22 C OPTION: 'INJE' -- Jacobian 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) C INTEGER MELEMF,MELEMC,MELECB,INORM,ICHPVO,ICHPSU, IROC,IVITC,IPC & ,IGAMC,ICHLIM,ICEL,NFAC,IFAC,MELRES,IJACO & ,NGF,NGC,NLF,NLC,NLCB & ,ILIINC,ILIINP,IJAC,II,JJ,K & ,MP, NBEL, NBME, NBSOUS, NKID, NKMT, NMATRI, NP, NRIGE & ,NSP,I, IYC,J, LRECP,LRECV,KV REAL*8 VOLU,SURF,RC,PC,UXC,UYC,GAMC,CNX,CNY & ,PSRF,RHOUF,DECIDP,DUNDPC,DUXDPC,DUYDPC & ,UN,RHO,P,COEF & ,BR1,BOT,TOP CHARACTER*(8) TYPE -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, MPLIM.MPOVAL, MPYC.MPOVAL -INC SMMATRIK POINTEUR CELL.IZAFM C------------------------------------------------------- -INC SMLREEL POINTEUR MLRECP.MLREEL, MLRECV.MLREEL C------------------------------------------------------- C********* Les Jacobians ****************************** C------------------------------------------------------- SEGMENT JACEL REAL*8 JAC(3+NSP,3+NSP) ENDSEGMENT POINTEUR JTL.JACEL, WL.JACEL, JTT.JACEL C------------------------------------------------------------- C******* Les fractionines massiques ************************** C------------------------------------------------------------- SEGMENT FRAMAS REAL*8 YET(NSP) ENDSEGMENT POINTEUR YC.FRAMAS, YF.FRAMAS C------------------------------------------------------- C********** Les CP's and CV's *********************** C------------------------------------------------------- SEGMENT GCONST REAL*8 GC(NSP) ENDSEGMENT POINTEUR CP.GCONST, CV.GCONST C------------------------------------------------------------- C********** Segments for the vectors *********************** C------------------------------------------------------------- SEGMENT VECEL REAL*8 VV(NSP) ENDSEGMENT POINTEUR DGDYC.VECEL C-------------------------------------------------- SEGINI JTL SEGINI JTT SEGINI WL 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(ICHLIM,MPLIM,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--------------------------------- C**** Objet MATRIK C--------------------------------- NRIGE = 7 NMATRI = 1 NKID = 9 NKMT = 7 C--------------------------------- SEGINI MATRIK IJACO = MATRIK MATRIK.IRIGEL(1,1) = MELRES MATRIK.IRIGEL(2,1) = MELRES C--------------------------------- C**** Matrice non symetrique C--------------------------------- MATRIK.IRIGEL(7,1) = 2 C--------------------------------- NBME = (3+NSP)*(3+NSP) NBSOUS = 1 SEGINI IMATRI IF(IJAC.EQ.1)THEN MLMOTS=ILIINC ELSEIF(IJAC.EQ.2)THEN MLMOTS=ILIINP ENDIF SEGACT MLMOTS MATRIK.IRIGEL(4,1) = IMATRI C------------------------------------------- DO 1 J=1,(NSP+3) KV=(J-1)*(3+NSP) IMATRI.LISPRI(KV+1) = MLMOTS.MOTS(1) IMATRI.LISPRI(KV+2) = MLMOTS.MOTS(2) IMATRI.LISPRI(KV+3) = MLMOTS.MOTS(3) IMATRI.LISPRI(KV+4) = MLMOTS.MOTS(4) DO 2 I=1,(NSP-1) IMATRI.LISPRI(KV+4+I) = MLMOTS.MOTS(4+I) 2 CONTINUE 1 CONTINUE C----------------------------------------------- SEGDES MLMOTS MLMOTS=ILIINC SEGACT MLMOTS C----------------------------------------------- DO 3 J=1,(NSP+3) KV=(J-1)*(3+NSP) IMATRI.LISDUA(KV+1) = MLMOTS.MOTS(j) IMATRI.LISDUA(KV+2) = MLMOTS.MOTS(j) IMATRI.LISDUA(KV+3) = MLMOTS.MOTS(j) IMATRI.LISDUA(KV+4) = MLMOTS.MOTS(j) DO 4 I=1,(NSP-1) IMATRI.LISDUA(KV+4+I) = MLMOTS.MOTS(j) 4 CONTINUE 3 CONTINUE C----------------------------------------------- C----------------------------------------------- SEGDES MLMOTS NBEL = NFAC NBSOUS = 1 NP = 1 MP = 1 C----------------------------------------------------------- C----------------------------------------------------------- DO 5 I=1,NBME SEGINI CELL IMATRI.LIZAFM(1,I) = CELL 5 CONTINUE C--------------------------------- C--------------------------------- C**** Fin definition MATRIK C--------------------------------- 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=-1*MPNORM.VPOCHA(NLF,1) CNY=-1*MPNORM.VPOCHA(NLF,2) 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) SEGINI YC SEGACT MPYC DO 100 I=1,(NSP-1) YC.YET(I)=MPYC.VPOCHA(NLC,I) 100 CONTINUE C--------------------------------- C Variables à la face C--------------------------------- RHOUF=MPLIM.VPOCHA(NLCB,1) PSRF=MPLIM.VPOCHA(NLCB,2) SEGINI YF SEGACT MPLIM DO 101 I=1,(NSP-1) YF.YET(I)=MPLIM.VPOCHA(NLCB,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------------------------------------------------------------- SEGINI DGDYC do 41 i=1,(nsp-1) dgdyc.vv(i)=(cp.gc(i)-cp.gc(nsp)- & GAMC*(cv.gc(i)-cv.gc(nsp)))/bot 41 continue c------------------------------------------------------------- C------------------------------------------------ C******* Densite, vitesse, pression sur le bord C------------------------------------------------ P=PC RHO=P/PSRF UN=RHOUF/RHO C------------------------------ C******* Derivatives w.r.t. PC C------------------------------ DUNDPC=-1*RHOUF/(RHO*RHO*PSRF) DUXDPC=CNX*DUNDPC DUYDPC=CNY*DUNDPC DECIDP=UN*DUNDPC C-------------------------------------------------------------- COEF=SURF/VOLU C------------------------------------------------------------- JTL.JAC(1,1) = 0.0D0 JTL.JAC(1,2) = 0.0D0 JTL.JAC(1,3) = 0.0D0 JTL.JAC(1,4) = 0.0D0 DO 107 I=1,(NSP-1) JTL.JAC(1,4+I) = 0.0D0 107 CONTINUE C---------------------------------------- JTL.JAC(2,1) = 0.0D0 JTL.JAC(2,2) = 0.0D0 JTL.JAC(2,3) = 0.0D0 JTL.JAC(2,4) = (RHOUF*DUXDPC+CNX)*COEF DO 108 I=1,(NSP-1) JTL.JAC(2,4+I) = 0.0D0 108 CONTINUE C---------------------------------------- JTL.JAC(3,1) = 0.0D0 JTL.JAC(3,2) = 0.0D0 JTL.JAC(3,3) = 0.0D0 JTL.JAC(3,4) = (RHOUF*DUYDPC+CNY)*COEF DO 109 I=1,(NSP-1) JTL.JAC(3,4+I) = 0.0D0 109 CONTINUE C---------------------------------------- JTL.JAC(4,1) = 0.0D0 JTL.JAC(4,2) = 0.0D0 JTL.JAC(4,3) = 0.0D0 JTL.JAC(4,4) = (RHOUF*DECIDP)*COEF DO 110 I=1,(NSP-1) JTL.JAC(4,4+I) = 0.0D0 110 CONTINUE C---------------------------------------- DO 111 I=1,(NSP-1) JTL.JAC(4+I,1) = 0.0D0 JTL.JAC(4+I,2) = 0.0D0 JTL.JAC(4+I,3) = 0.0D0 JTL.JAC(4+I,4) = 0.0D0 DO 112 J=1,(NSP-1) JTL.JAC(4+I,4+J) = 0.0D0 112 CONTINUE 111 CONTINUE C--------------------------------------------------------- c matrix wl(i,j) represents the derivative of the i-component c of the vector of primitive variables of the left state with c respect to the j-component of the vector of the conservative c variables of the left state. c c Here: (rho, ux, uy, p, Y_1,...,Y_(nsp-1)) - c vector of primitive variables; c (rho, rho ux, rho uy, rho e, rho Y_1,..., rho Y_(nsp-1)) - c vector of conservative variables. c------------------------------------------------------------- wl.jac(1,1)=1.0d0 wl.jac(1,2)=0.0d0 wl.jac(1,3)=0.0d0 wl.jac(1,4)=0.0d0 do 83 i=1,(nsp-1) wl.jac(1,4+i)=0.0d0 83 continue c------------------------------ wl.jac(2,1)=-UXC/RC wl.jac(2,2)=1.0d0/RC wl.jac(2,3)=0.0d0 wl.jac(2,4)=0.0d0 do 84 i=1,(nsp-1) wl.jac(2,4+i)=0.0d0 84 continue c------------------------------ wl.jac(3,1)=-UYC/RC wl.jac(3,2)=0.0d0 wl.jac(3,3)=1.0d0/RC wl.jac(3,4)=0.0d0 do 85 i=1,(nsp-1) wl.jac(3,4+i)=0.0d0 85 continue c------------------------------ br1=0.0d0 do 86 i=1,(nsp-1) br1=br1+dgdyc.vv(i)*yc.yet(i) 86 continue br1=br1*PC/(RC*(GAMC-1.0D0)) wl.jac(4,1)=(GAMC-1.0D0)*(UXC*UXC+UYC*UYC)/2.0d0-br1 wl.jac(4,2)=-UXC*(GAMC-1.0D0) wl.jac(4,3)=-UYC*(GAMC-1.0D0) wl.jac(4,4)=(GAMC-1.0D0) do 87 i=1,(nsp-1) wl.jac(4,4+i)=dgdyc.vv(i)*PC/(RC*(GAMC-1.0D0)) 87 continue c------------------------------ do 88 i=1,(nsp-1) do 89 j=1,4 wl.jac(4+i,j)=0.0d0 if(j.eq.1) wl.jac(4+i,j)=-yc.yet(i)/RC 89 continue c------------ do 890 j=5,(4+nsp-1) wl.jac(4+i,j)=0.0d0 if(4+i.eq.j) then wl.jac(4+i,j)=1.0d0/RC endif 890 continue 88 continue c------------------------------------------------ C------------------------------------------------ do 114 i=1,(3+nsp) do 115 j=1,(3+nsp) jtt.jac(i,j)=0.0d0 do 116 k=1,(3+nsp) jtt.jac(i,j)=jtt.jac(i,j)+jtl.jac(i,k)*wl.jac(k,j) 116 continue 115 continue 114 continue C---------------------------------------------------------------- C******* Jacobian with respect to conservative variables C---------------------------------------------------------------- IF(IJAC.EQ.1)THEN DO 9 II = 1,(3+NSP) DO 15 JJ = 1,(3+NSP) KV = (II-1)*(3+NSP) C---------------------------------- CELL = IMATRI.LIZAFM(1,KV+JJ) CELL.AM(IFAC,1,1) = JTT.JAC(II,JJ) 15 CONTINUE 9 CONTINUE ELSEIF(IJAC.EQ.2)THEN DO 20 II = 1,(3+NSP) DO 25 JJ = 1,(3+NSP) KV = (II-1)*(3+NSP) C---------------------------------- CELL = IMATRI.LIZAFM(1,KV+JJ) CELL.AM(IFAC,1,1) = JTL.JAC(II,JJ) 25 CONTINUE 20 CONTINUE ENDIF c-------------------------------------------------- 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 MPYC SEGDES MPLIM SEGDES YC SEGDES YF SEGDES CP SEGDES CV SEGDES JTL SEGDES JTT SEGDES WL SEGDES DGDYC SEGDES MATRIK DO 80 II=1,NBME CELL = IMATRI.LIZAFM(1,II) SEGDES CELL 80 CONTINUE SEGDES IMATRI C--------------------------------------------- 9999 CONTINUE RETURN END