C JBMMW2 SOURCE CB215821 20/11/25 13:30:39 10792 SUBROUTINE JBMMW2(NSP,JLL,WVEC_L,WVEC_R,NVECT,TVECT, & mpyn,lrecp,lrecv,nlcg,v_inf) C************************************************************************ C C PROJET : CASTEM 2000 C C NOM : JACBM C C DESCRIPTION : Voir KONJA6 C C LANGAGE : FORTRAN 77 + ESOPE 2000 (avec estensions CISI) C C AUTEUR : S. KUDRIAKOV, DM2S/SFME/LTMF C C************************************************************************ C c---------------------------------------------------------------------- c GENERAL DESCRIPTION: c This subroutine provides the jacobians which are the derivatives c of the numerical flux function defined at the cell interface c with respect to the conservative variables of the left and right c cells (relative to the cell interface). c The low-mach number corrections are made for the flux functions c c EQUATIONS: 2D Euler equations of gas dynamics c c c REFERENCE: 1) JCP, 129, 364-382 (1996) c " A Sequel to AUSM: AUSM+ ". c M.S.Liou c 2) AIAA Journal, Sept. 1998 c " Low-Diffusion Flux-Splitting Methods for Flows at All Speeds" c J.R.Edwards and M.S.Liou c---------------------------------------------------------------------- c INPUT: c c alpha -- parameter of the AUSM+ scheme in the Pressure function; c ( -3/4 <= alpha <= 3/16 ) (imposed as a parameter) c c beta -- parameter of the AUSM+ scheme in the Mach function; c ( -1/16 <= beta <= 1/2 ) (imposed as a parameter) c c wvec_l -- vector of the primitive variables (rho,ux,uy,p) at the c left cell; c c wvec_r -- vector of the primitive variables (rho,ux,uy,p) at the c right cell; c c nvect -- normal vector to the interface (2 components in 2D); c c tvect -- tangential vector to the interface; c c ga -- ratio of the specific heats (assumed constant) c c v_inf -- parameter for reference velocity, see Ref. 2). c---------------------------------------------------------------------- c c OUTPUT: c c jtl -- jakobian matrix 4 by 4 - derivatives of the numerical c flux function with respect to the conservative variables c from the left cell; c c jtr -- jakobian matrix 4 by 4 - derivatives of the numerical c flux function with respect to the conservative variables c from the right cell. c---------------------------------------------------------------------- IMPLICIT INTEGER(I-N) real*8 wvec_l(4),wvec_r(4) real*8 nvect(2),tvect(2) real*8 alpha,beta,bcoef,tcoef real*8 gal,gar,gm1l,gm1r,temph real*8 n_x,n_y,t_x,t_y real*8 un_l, un_r, ut_l, ut_r real*8 ml,mr,Mplus,Mmin,mmid real*8 mpl_m, mmin_m,am real*8 rold_l,uold_l,vold_l,pold_l,eold_l real*8 rold_r,uold_r,vold_r,pold_r,eold_r real*8 Pplus,Pmin real*8 top,bot,bots real*8 br1,br3,br4,temp_l,temp_r,brac_l,brac_r real*8 aleft, arigh real*8 damr_l,damr_r,damu_l,damu_r real*8 damv_l,damv_r,damp_l,damp_r real*8 damg_l,damg_r real*8 dmlr_l,dmlr_r,dmlu_l,dmlu_r real*8 dmlv_l,dmlv_r,dmlp_l,dmlp_r real*8 dmrr_l,dmrr_r,dmru_l,dmru_r real*8 dmrv_l,dmrv_r,dmrp_l,dmrp_r real*8 dMpr_l,dMpr_r,dMpu_l,dMpu_r real*8 dMpv_l,dMpv_r,dMpp_l,dMpp_r real*8 dMmr_l,dMmr_r,dMmu_l,dMmu_r real*8 dMmv_l,dMmv_r,dMmp_l,dMmp_r real*8 dmir_l,dmir_r,dmiu_l,dmiu_r real*8 dmiv_l,dmiv_r,dmip_l,dmip_r real*8 d3mr_l,d3mr_r,d3mu_l,d3mu_r real*8 d3mv_l,d3mv_r,d3mp_l,d3mp_r real*8 d2mr_l,d2mr_r,d2mu_l,d2mu_r real*8 d2mv_l,d2mv_r,d2mp_l,d2mp_r real*8 dPpr_l,dPpr_r,dPpu_l,dPpu_r real*8 dPpv_l,dPpv_r,dPpp_l,dPpp_r real*8 dPmr_l,dPmr_r,dPmu_l,dPmu_r real*8 dPmv_l,dPmv_r,dPmp_l,dPmp_r real*8 dpir_l,dpir_r,dpiu_l,dpiu_r real*8 dpiv_l,dpiv_r,dpip_l,dpip_r real*8 epsil,qq,amw,Mmin1,Mplus1 real*8 fmid,mlw,mrw,termp real*8 ur_r,ur_l,urm,mhalf,mhalfr real*8 durr_l,durr_r,duru_l,duru_r real*8 durv_l,durv_r,durp_l,durp_r real*8 dmhr_l,dmhr_r,dmhu_l,dmhu_r real*8 dmhv_l,dmhv_r,dmhp_l,dmhp_r real*8 dmfr_l,dmfr_r,dmfu_l,dmfu_r real*8 dmfv_l,dmfv_r,dmfp_l,dmfp_r real*8 dfm_mf,dfm_mh real*8 dfmr_l,dfmr_r,dfmu_l,dfmu_r real*8 dfmv_l,dfmv_r,dfmp_l,dfmp_r real*8 m1mr_l,m1mr_r,m1mu_l,m1mu_r real*8 m1mv_l,m1mv_r,m1mp_l,m1mp_r real*8 m1pr_l,m1pr_r,m1pu_l,m1pu_r real*8 m1pv_l,m1pv_r,m1pp_l,m1pp_r real*8 tmpr_l,tmpr_r,tmpu_l,tmpu_r real*8 tmpv_l,tmpv_r,tmpp_l,tmpp_r real*8 coef,canc,v_inf real*8 sr_l,sr_r,su_l,su_r,sv_l,sv_r,sp_l,sp_r real*8 rum,rumr_l,rumu_l,rumv_l,rump_l real*8 rumr_r,rumu_r,rumv_r,rump_r real*8 dfmg_l,dfmg_r,dmfg_l,dmfg_r real*8 dmhg_l,dmhg_r,durg_l,durg_r integer i,j,k,jll,lrecp,lrecv,nlcd,nlcg,nsp parameter(alpha = 0.1875D0, beta = 0.125D0) parameter(epsil = 1.0d0) canc=1.0d0 c------------------------------------------------------------- -INC SMCHPOI POINTEUR MPYN.MPOVAL C------------------------------------------------------------- -INC SMLREEL POINTEUR MLRECP.MLREEL, MLRECV.MLREEL C------------------------------------------------------------- C******* Les fractionines massiques ************************** C------------------------------------------------------------- SEGMENT FRAMAS REAL*8 YET(NSP) ENDSEGMENT POINTEUR YL.FRAMAS, YR.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 elementary matrixes ************* C------------------------------------------------------------- SEGMENT JACEL REAL*8 JAC(3+NSP,3+NSP) ENDSEGMENT POINTEUR JTL.JACEL, JTR.JACEL, JL.JACEL, JR.JACEL, & WL.JACEL, WR.JACEL c---------------------------------------- SEGINI JTL SEGINI JTR SEGINI JL SEGINI JR SEGINI WL SEGINI WR C------------------------------------------------------------- C********** Segments for the vectors *********************** C------------------------------------------------------------- SEGMENT VECEL REAL*8 VV(NSP) ENDSEGMENT POINTEUR DMLY_L.VECEL, DMLY_R.VECEL, & dmry_l.vecel, dmry_r.vecel, & dMpy_l.vecel, dMpy_r.vecel, & dMmy_l.vecel, dMmy_r.vecel, & dmiy_l.vecel, dmiy_r.vecel, & d3my_l.vecel, d3my_r.vecel, & d2my_l.vecel, d2my_r.vecel, & dPpy_l.vecel, dPpy_r.vecel, & dPmy_l.vecel, dPmy_r.vecel, & dpiy_l.vecel, dpiy_r.vecel, & dgdyl.vecel, dgdyr.vecel, & sy_l.vecel, sy_r.vecel, & m1py_l.vecel, m1py_r.vecel, & m1my_l.vecel, m1my_r.vecel, & tmpy_l.vecel, tmpy_r.vecel, & rumy_l.vecel, rumy_r.vecel C---------------------------------------------- SEGINI DMLY_L, DMLY_R, & dmry_l, dmry_r, & dMpy_l, dMpy_r, & dMmy_l, dMmy_r, & dmiy_l, dmiy_r, & d3my_l, d3my_r, & d2my_l, d2my_r, & dPpy_l, dPpy_r, & dPmy_l, dPmy_r, & dpiy_l, dpiy_r, & dgdyl, dgdyr, & sy_l, sy_r, & m1py_l, m1py_r, & m1my_l, m1my_r, & tmpy_l, tmpy_r, & rumy_l, rumy_r C------------------------------------------------------------- C********** Segments for the flux-vector ******************* C------------------------------------------------------------- SEGMENT FUNEL REAL*8 FU(3+NSP) ENDSEGMENT POINTEUR f.funel C------------------------- SEGINI f C------------------------------------------------------------ SEGINI YL, YR SEGACT MPYN DO 100 I=1,(NSP-1) YL.YET(I)=MPYN.VPOCHA(NLCG,I) YR.YET(I)=MPYN.VPOCHA(NLCG,I) 100 CONTINUE C---------------------------------------- SEGINI CP, CV MLRECP = LRECP MLRECV = LRECV SEGACT MLRECP, MLRECV DO 101 I=1,(NSP-1) CP.GC(I)=MLRECP.PROG(I) CV.GC(I)=MLRECV.PROG(I) 101 CONTINUE CP.GC(NSP)=MLRECP.PROG(NSP) CV.GC(NSP)=MLRECV.PROG(NSP) c------------------------------------------------------------- c Computing GAMMA at the left cell and its derivatives c with respect to the primitive variables Y_i c------------------------------------------------------------- top=0.0D0 bot=0.0D0 do 40 i=1,(nsp-1) top=top+yl.yet(i)*(cp.gc(i)-cp.gc(nsp)) bot=bot+yl.yet(i)*(cv.gc(i)-cv.gc(nsp)) 40 continue top=cp.gc(nsp)+top bot=cv.gc(nsp)+bot gal=top/bot gm1l=gal-1.0d0 c------------------------------------------------------------- do 41 i=1,(nsp-1) dgdyl.vv(i)=(cp.gc(i)-cp.gc(nsp)- & gal*(cv.gc(i)-cv.gc(nsp)))/bot 41 continue c------------------------------------------------------------- c Computing GAMMA at the right cell and its derivatives c with respect to the primitive variables Y_i c------------------------------------------------------------- top=0.0D0 bot=0.0D0 do 42 i=1,(nsp-1) top=top+yr.yet(i)*(cp.gc(i)-cp.gc(nsp)) bot=bot+yr.yet(i)*(cv.gc(i)-cv.gc(nsp)) 42 continue top=cp.gc(nsp)+top bot=cv.gc(nsp)+bot gar=top/bot gm1r=gar-1.0d0 c------------------------------------------------------------- do 43 i=1,(nsp-1) dgdyr.vv(i)=(cp.gc(i)-cp.gc(nsp)- & gar*(cv.gc(i)-cv.gc(nsp)))/bot 43 continue c------------------------------------------------------------- n_x=nvect(1) n_y=nvect(2) t_x=tvect(1) t_y=tvect(2) c---------------------------- rold_l=wvec_l(1) uold_l=wvec_l(2) vold_l=wvec_l(3) pold_l=wvec_l(4) c----------------------- rold_r=wvec_r(1) uold_r=wvec_r(2) vold_r=wvec_r(3) pold_r=wvec_r(4) c------------------------------------------------------------------ c Computation of the specific total energy on the left and right. c------------------------------------------------------------------ eold_l=(uold_l*uold_l+vold_l*vold_l)/2.0d0 eold_l=eold_l+pold_l/(gm1l*rold_l) eold_r=(uold_r*uold_r+vold_r*vold_r)/2.0d0 eold_r=eold_r+pold_r/(gm1r*rold_r) c------------------------------------------------------------------ c Computing reference velocity and its derivatives c------------------------------------------------------------------ aleft=sqrt(gal*pold_l/rold_l) arigh=sqrt(gar*pold_r/rold_r) qq=sqrt(uold_l*uold_l+vold_l*vold_l) if(qq .lt. (epsil*v_inf)) then ur_l = epsil*v_inf durr_l=0.0d0 duru_l=0.0d0 durv_l=0.0d0 durp_l=0.0d0 durg_l=0.0d0 else ur_l=qq durr_l=0.0d0 duru_l=uold_l/qq durv_l=vold_l/qq durp_l=0.0d0 durg_l=0.0d0 endif c------------------------------ if(ur_l .ge. aleft) then ur_l=aleft durr_l=-aleft/(2.0d0*rold_l) duru_l=0.0d0 durv_l=0.0d0 durp_l=gal/(2.0d0*aleft*rold_l) durg_l=aleft/(2.0d0*gal) endif c------------------------------------------------------------------ qq=sqrt(uold_r*uold_r+vold_r*vold_r) if(qq .lt. (epsil*v_inf)) then ur_r = epsil*v_inf durr_r=0.0d0 duru_r=0.0d0 durv_r=0.0d0 durp_r=0.0d0 durg_r=0.0d0 else ur_r=qq durr_r=0.0d0 duru_r=uold_r/qq durv_r=vold_r/qq durp_r=0.0d0 durg_r=0.0d0 endif c--------------------------- if(ur_r .ge. arigh) then ur_r=arigh durr_r=-arigh/(2.0d0*rold_r) duru_r=0.0d0 durv_r=0.0d0 durp_r=gar/(2.0d0*arigh*rold_r) durg_r=arigh/(2.0d0*gar) endif c------------------------------------------------------------------ c------------------------------------------------------------------- c Reference velocity at the interface is taken as an average c of the reference velocities of the neighbouring cells c------------------------------------------------------------------- urm=0.5d0*(ur_l+ur_r) durr_l=0.5d0*durr_l duru_l=0.5d0*duru_l durv_l=0.5d0*durv_l durp_l=0.5d0*durp_l durg_l=0.5d0*durg_l c------------------------- durr_r=0.5d0*durr_r duru_r=0.5d0*duru_r durv_r=0.5d0*durv_r durp_r=0.5d0*durp_r durg_r=0.5d0*durg_r c------------------------------------------------------------------- c Computation of the speed of sound and its derivatives; c numerical speed of sound at the interface is taken as an average c of the speeds of sounds of the neighbouring cells c------------------------------------------------------------------- am=0.5d0*(aleft+arigh) c------------------------------------------------------------------- if(abs(urm/am-1.0d0).le.0.000001d0) then coef=0.0d0 else coef=1.0d0 endif c------------------------------------------------------------------- damr_r=-arigh/(4.0d0*rold_r) damu_r=0.0d0 damv_r=0.0d0 damp_r=gar/(4.0d0*arigh*rold_r) damg_r=arigh/(4.0d0*gar) c----------------------- damr_l=-aleft/(4.0d0*rold_l) damu_l=0.0d0 damv_l=0.0d0 damp_l=gal/(4.0d0*aleft*rold_l) damg_l=aleft/(4.0d0*gal) c------------------------------------------------------------------- c Computing numerical Mach number and its derivatives, c see p.370, under (A1). c------------------------------------------------------------------- un_l=uold_l*n_x+vold_l*n_y un_r=uold_r*n_x+vold_r*n_y ut_l=uold_l*t_x+vold_l*t_y ut_r=uold_r*t_x+vold_r*t_y c------------------------------------------------------------------- ml=un_l/am mr=un_r/am mhalf=0.5d0*(un_l+un_r)/am c--------------------- top=0.5d0*(un_l+un_r)/(am*am) dmhr_l=-top*damr_l dmhu_l=n_x/2.0d0/am-top*damu_l dmhv_l=n_y/2.0d0/am-top*damv_l dmhp_l=-top*damp_l dmhg_l=-top*damg_l c--------------------- dmhr_r=-top*damr_r dmhu_r=n_x/2.0d0/am-top*damu_r dmhv_r=n_y/2.0d0/am-top*damv_r dmhp_r=-top*damp_r dmhg_r=-top*damg_r c-------------------------------- mhalfr=urm/am c--------------------- top=urm/(am*am) dmfr_l=durr_l/am-top*damr_l dmfu_l=duru_l/am-top*damu_l dmfv_l=durv_l/am-top*damv_l dmfp_l=durp_l/am-top*damp_l dmfg_l=durg_l/am-top*damg_l c--------------------- dmfr_r=durr_r/am-top*damr_r dmfu_r=duru_r/am-top*damu_r dmfv_r=durv_r/am-top*damv_r dmfp_r=durp_r/am-top*damp_r dmfg_r=durg_r/am-top*damg_r c------------------------------------------------------------------- c Scaling function for the speed of sound and its derivatives c------------------------------------------------------------------- top=(1.0d0-mhalfr*mhalfr)*(1.0d0-mhalfr*mhalfr) top=top*mhalf*mhalf+4.0d0*mhalfr*mhalfr bot=1.0d0+mhalfr*mhalfr if(abs(canc-0.0d0).lt.0.000001d0) then fmid=1.0d0 else fmid=sqrt(top)/bot endif c-------------------------- temph=-4.0d0*(1.0d0-mhalfr*mhalfr)*mhalfr temph=temph*mhalf*mhalf+8.0d0*mhalfr dfm_mf=temph/(sqrt(top)*2.0d0*bot) dfm_mf=dfm_mf-sqrt(top)*2.0d0*mhalfr/(bot*bot) c-------------------------- temph=2.0d0*(1.0d0-mhalfr*mhalfr)*(1.0d0-mhalfr*mhalfr)*mhalf dfm_mh=temph/(2.0d0*bot*sqrt(top)) c-------------------------- dfmr_l=dfm_mf*dmfr_l+dfm_mh*dmhr_l dfmu_l=dfm_mf*dmfu_l+dfm_mh*dmhu_l dfmv_l=dfm_mf*dmfv_l+dfm_mh*dmhv_l dfmp_l=dfm_mf*dmfp_l+dfm_mh*dmhp_l dfmg_l=dfm_mf*dmfg_l+dfm_mh*dmhg_l c-------------------------- dfmr_r=dfm_mf*dmfr_r+dfm_mh*dmhr_r dfmu_r=dfm_mf*dmfu_r+dfm_mh*dmhu_r dfmv_r=dfm_mf*dmfv_r+dfm_mh*dmhv_r dfmp_r=dfm_mf*dmfp_r+dfm_mh*dmhp_r dfmg_r=dfm_mf*dmfg_r+dfm_mh*dmhg_r c-------------------------- amw=fmid*am mlw=un_l/amw mrw=un_r/amw c-------------------------- damr_l=canc*coef*dfmr_l*am+fmid*damr_l damu_l=canc*coef*dfmu_l*am+fmid*damu_l damv_l=canc*coef*dfmv_l*am+fmid*damv_l damp_l=canc*coef*dfmp_l*am+fmid*damp_l damg_l=canc*coef*dfmg_l*am+fmid*damg_l c-------------------------- damr_r=canc*coef*dfmr_r*am+fmid*damr_r damu_r=canc*coef*dfmu_r*am+fmid*damu_r damv_r=canc*coef*dfmv_r*am+fmid*damv_r damp_r=canc*coef*dfmp_r*am+fmid*damp_r damg_r=canc*coef*dfmg_r*am+fmid*damg_r c-/-/-/-/-/-/-/-/-/-/-/-/-/-/-/-/-/ am=amw c------------------------------------------------------------------- c Redefinition of the numerical mach numbers c------------------------------------------------------------------- if(abs(canc-0.0d0).lt.0.000001d0) then top=2.0d0 bot=0.0d0 else top=1.0d0+mhalfr*mhalfr bot=1.0d0-mhalfr*mhalfr endif ml=0.5d0*(top*mlw+bot*mrw) mr=0.5d0*(top*mrw+bot*mlw) c------------------------------------------------------------------- c Mplus and Mmin are calligraphic lettes M+ and M- from the paper, c see (19a) and (19b), p.367. c------------------------------------------------------------------- if(abs(ml) .ge. 1.0d0) then Mplus=(ml+abs(ml))/2.0d0 else Mplus=(ml+1.0d0)*(ml+1.0d0)/4.0d0 Mplus=Mplus+beta*(ml*ml-1.0d0)*(ml*ml-1.0d0) endif Mplus1=(ml+abs(ml))/2.0d0 c------------------------------------------------------------------- if(abs(mr) .ge. 1.0d0) then Mmin=(mr-abs(mr))/2.0d0 else Mmin=-(mr-1.0d0)*(mr-1.0d0)/4.0d0 Mmin=Mmin-beta*(mr*mr-1.0d0)*(mr*mr-1.0d0) endif Mmin1=(mr-abs(mr))/2.0d0 c------------------------------------------------------------------- c Derivatives of ml and mr with respect to both sets of primitive c variables: left and right. c------------------------------------------------------------------- temp_l=-un_l/(am*am) temp_r=-un_r/(am*am) c-------- dmlr_l=temp_l*damr_l dmlr_r=temp_l*damr_r dmrr_l=temp_r*damr_l dmrr_r=temp_r*damr_r c-------- dmlu_l=n_x/am+temp_l*damu_l dmlu_r=temp_l*damu_r dmru_l=temp_r*damu_l dmru_r=n_x/am+temp_r*damu_r c-------- dmlv_l=n_y/am+temp_l*damv_l dmlv_r=temp_l*damv_r dmrv_l=temp_r*damv_l dmrv_r=n_y/am+temp_r*damv_r c-------- dmlp_l=temp_l*damp_l dmlp_r=temp_l*damp_r dmrp_l=temp_r*damp_l dmrp_r=temp_r*damp_r c-------- do 44 i=1,(nsp-1) dmly_l.vv(i)=temp_l*damg_l*dgdyl.vv(i) dmly_r.vv(i)=temp_l*damg_r*dgdyr.vv(i) dmry_l.vv(i)=temp_r*damg_l*dgdyl.vv(i) dmry_r.vv(i)=temp_r*damg_r*dgdyr.vv(i) 44 continue c----------------------------- sr_l=dmlr_l sr_r=dmlr_r su_l=dmlu_l su_r=dmlu_r sv_l=dmlv_l sv_r=dmlv_r sp_l=dmlp_l sp_r=dmlp_r do 440 i=1,(nsp-1) sy_l.vv(i) = dmly_l.vv(i) sy_r.vv(i) = dmly_r.vv(i) 440 continue c----------------------------------------------------------------- c Redefinition of the derivatives of the ml & mr c----------------------------------------------------------------- temp_l=(mlw-mrw)*mhalfr temp_r=-temp_l c-------- dmlr_l=0.5d0*(top*dmlr_l+bot*dmrr_l)+canc*coef*temp_l*dmfr_l dmlu_l=0.5d0*(top*dmlu_l+bot*dmru_l)+canc*coef*temp_l*dmfu_l dmlv_l=0.5d0*(top*dmlv_l+bot*dmrv_l)+canc*coef*temp_l*dmfv_l dmlp_l=0.5d0*(top*dmlp_l+bot*dmrp_l)+canc*coef*temp_l*dmfp_l c-------- dmlr_r=0.5d0*(top*dmlr_r+bot*dmrr_r)+canc*coef*temp_l*dmfr_r dmlu_r=0.5d0*(top*dmlu_r+bot*dmru_r)+canc*coef*temp_l*dmfu_r dmlv_r=0.5d0*(top*dmlv_r+bot*dmrv_r)+canc*coef*temp_l*dmfv_r dmlp_r=0.5d0*(top*dmlp_r+bot*dmrp_r)+canc*coef*temp_l*dmfp_r c-------- do 441 i=1,(nsp-1) dmly_l.vv(i)=0.5d0*(top*dmly_l.vv(i)+bot*dmry_l.vv(i))+ & canc*coef*temp_l*dmfg_l*dgdyl.vv(i) dmly_r.vv(i)=0.5d0*(top*dmly_r.vv(i)+bot*dmry_r.vv(i))+ & canc*coef*temp_l*dmfg_r*dgdyr.vv(i) 441 continue c-------- dmrr_l=0.5d0*(top*dmrr_l+bot*sr_l)+canc*coef*temp_r*dmfr_l dmru_l=0.5d0*(top*dmru_l+bot*su_l)+canc*coef*temp_r*dmfu_l dmrv_l=0.5d0*(top*dmrv_l+bot*sv_l)+canc*coef*temp_r*dmfv_l dmrp_l=0.5d0*(top*dmrp_l+bot*sp_l)+canc*coef*temp_r*dmfp_l c-------- dmrr_r=0.5d0*(top*dmrr_r+bot*sr_r)+canc*coef*temp_r*dmfr_r dmru_r=0.5d0*(top*dmru_r+bot*su_r)+canc*coef*temp_r*dmfu_r dmrv_r=0.5d0*(top*dmrv_r+bot*sv_r)+canc*coef*temp_r*dmfv_r dmrp_r=0.5d0*(top*dmrp_r+bot*sp_r)+canc*coef*temp_r*dmfp_r do 442 i=1,(nsp-1) dmry_l.vv(i)=0.5d0*(top*dmry_l.vv(i)+bot*sy_l.vv(i))+ & canc*coef*temp_r*dmfg_l*dgdyl.vv(i) dmry_r.vv(i)=0.5d0*(top*dmry_r.vv(i)+bot*sy_r.vv(i))+ & canc*coef*temp_r*dmfg_r*dgdyr.vv(i) 442 continue c----------------------------------------------------------- c mmid is m_{1/2} (notation as in the paper, see (13),p.366) c----------------------------------------------------------- mmid=Mplus+Mmin c----------------------------------------------------------- c Computing the derivatives of M+ and M- c----------------------------------------------------------- if(ml .ge. 1.0d0) then dMpr_l=dmlr_l dMpu_l=dmlu_l dMpv_l=dmlv_l dMpp_l=dmlp_l do 45 i=1,(nsp-1) dMpy_l.vv(i)=dmly_l.vv(i) 45 continue c-------------------- dMpr_r=dmlr_r dMpu_r=dmlu_r dMpv_r=dmlv_r dMpp_r=dmlp_r do 46 i=1,(nsp-1) dMpy_r.vv(i)=dmly_r.vv(i) 46 continue else if((ml .gt. -1.0d0) .and. (ml .lt. 1.0d0)) then temph=(ml+1.0d0)/2.0d0 dMpr_l=(temph+4.0d0*beta*ml*(ml*ml-1.0d0))*dmlr_l dMpu_l=(temph+4.0d0*beta*ml*(ml*ml-1.0d0))*dmlu_l dMpv_l=(temph+4.0d0*beta*ml*(ml*ml-1.0d0))*dmlv_l dMpp_l=(temph+4.0d0*beta*ml*(ml*ml-1.0d0))*dmlp_l do 47 i=1,(nsp-1) dMpy_l.vv(i)=(temph+4.0d0*beta*ml* & (ml*ml-1.0d0))*dmly_l.vv(i) 47 continue c-------------------- dMpr_r=(temph+4.0d0*beta*ml*(ml*ml-1.0d0))*dmlr_r dMpu_r=(temph+4.0d0*beta*ml*(ml*ml-1.0d0))*dmlu_r dMpv_r=(temph+4.0d0*beta*ml*(ml*ml-1.0d0))*dmlv_r dMpp_r=(temph+4.0d0*beta*ml*(ml*ml-1.0d0))*dmlp_r do 48 i=1,(nsp-1) dMpy_r.vv(i)=(temph+4.0d0*beta*ml* & (ml*ml-1.0d0))*dmly_r.vv(i) 48 continue else dMpr_l=0.0d0 dMpu_l=0.0d0 dMpv_l=0.0d0 dMpp_l=0.0d0 do 49 i=1,(nsp-1) dMpy_l.vv(i)=0.0d0 49 continue c--------------------- dMpr_r=0.0d0 dMpu_r=0.0d0 dMpv_r=0.0d0 dMpp_r=0.0d0 do 50 i=1,(nsp-1) dMpy_r.vv(i)=0.0d0 50 continue endif endif c----------------------------------------------------------- c addition of low Mach number c----------------------------------------------------------- if(ml .ge. 0.0d0) then m1pr_l=dmlr_l m1pu_l=dmlu_l m1pv_l=dmlv_l m1pp_l=dmlp_l do 450 i=1,(nsp-1) m1py_l.vv(i)=dmly_l.vv(i) 450 continue c--------------------- m1pr_r=dmlr_r m1pu_r=dmlu_r m1pv_r=dmlv_r m1pp_r=dmlp_r do 460 i=1,(nsp-1) m1py_r.vv(i)=dmly_r.vv(i) 460 continue else m1pr_l=0.0d0 m1pu_l=0.0d0 m1pv_l=0.0d0 m1pp_l=0.0d0 do 451 i=1,(nsp-1) m1py_l.vv(i)=0.0d0 451 continue c-------------------- m1pr_r=0.0d0 m1pu_r=0.0d0 m1pv_r=0.0d0 m1pp_r=0.0d0 do 461 i=1,(nsp-1) m1py_r.vv(i)=0.0d0 461 continue endif c----------------------------------------------------------- if(mr .ge. 1.0d0) then dMmr_l=0.0d0 dMmu_l=0.0d0 dMmv_l=0.0d0 dMmp_l=0.0d0 do 51 i=1,(nsp-1) dMmy_l.vv(i)=0.0d0 51 continue c--------------------- dMmr_r=0.0d0 dMmu_r=0.0d0 dMmv_r=0.0d0 dMmp_r=0.0d0 do 52 i=1,(nsp-1) dMmy_r.vv(i)=0.0d0 52 continue else if((mr .gt. -1.0d0) .and. (mr .lt. 1.0d0)) then temph=(-mr+1.0d0)/2.0d0 dMmr_l=(temph-4.0d0*beta*mr*(mr*mr-1.0d0))*dmrr_l dMmu_l=(temph-4.0d0*beta*mr*(mr*mr-1.0d0))*dmru_l dMmv_l=(temph-4.0d0*beta*mr*(mr*mr-1.0d0))*dmrv_l dMmp_l=(temph-4.0d0*beta*mr*(mr*mr-1.0d0))*dmrp_l do 53 i=1,(nsp-1) dMmy_l.vv(i)=(temph-4.0d0*beta*mr* & (mr*mr-1.0d0))*dmry_l.vv(i) 53 continue c-------------------- dMmr_r=(temph-4.0d0*beta*mr*(mr*mr-1.0d0))*dmrr_r dMmu_r=(temph-4.0d0*beta*mr*(mr*mr-1.0d0))*dmru_r dMmv_r=(temph-4.0d0*beta*mr*(mr*mr-1.0d0))*dmrv_r dMmp_r=(temph-4.0d0*beta*mr*(mr*mr-1.0d0))*dmrp_r do 54 i=1,(nsp-1) dMmy_r.vv(i)=(temph-4.0d0*beta*mr* & (mr*mr-1.0d0))*dmry_r.vv(i) 54 continue else dMmr_l=dmrr_l dMmu_l=dmru_l dMmv_l=dmrv_l dMmp_l=dmrp_l do 55 i=1,(nsp-1) dMmy_l.vv(i)=dmry_l.vv(i) 55 continue c-------------------- dMmr_r=dmrr_r dMmu_r=dmru_r dMmv_r=dmrv_r dMmp_r=dmrp_r do 56 i=1,(nsp-1) dMmy_r.vv(i)=dmry_r.vv(i) 56 continue endif endif c----------------------------------------------------------- c addition of low Mach number c----------------------------------------------------------- if(mr .le. 0.0d0) then m1mr_l=dmrr_l m1mu_l=dmru_l m1mv_l=dmrv_l m1mp_l=dmrp_l do 550 i=1,(nsp-1) m1my_l.vv(i)=dmry_l.vv(i) 550 continue c--------------------- m1mr_r=dmrr_r m1mu_r=dmru_r m1mv_r=dmrv_r m1mp_r=dmrp_r do 560 i=1,(nsp-1) m1my_r.vv(i)=dmry_r.vv(i) 560 continue else m1mr_l=0.0d0 m1mu_l=0.0d0 m1mv_l=0.0d0 m1mp_l=0.0d0 do 551 i=1,(nsp-1) m1my_l.vv(i)=0.0d0 551 continue c-------------------- m1mr_r=0.0d0 m1mu_r=0.0d0 m1mv_r=0.0d0 m1mp_r=0.0d0 do 561 i=1,(nsp-1) m1my_r.vv(i)=0.0d0 561 continue endif c----------------------------------------------------------------- c computing the derivatives of m_{1/2} (notation as in the paper) c----------------------------------------------------------------- dmir_l=dMpr_l+dMmr_l dmir_r=dMpr_r+dMmr_r c------------- dmiu_l=dMpu_l+dMmu_l dmiu_r=dMpu_r+dMmu_r c------------- dmiv_l=dMpv_l+dMmv_l dmiv_r=dMpv_r+dMmv_r c------------- dmip_l=dMpp_l+dMmp_l dmip_r=dMpp_r+dMmp_r c------------- do 57 i=1,(nsp-1) dmiy_l.vv(i)=dMpy_l.vv(i)+dMmy_l.vv(i) dmiy_r.vv(i)=dMpy_r.vv(i)+dMmy_r.vv(i) 57 continue c---------------------------------------------------------------- c computing the main convective variables and their derivatives c mpl_m is m^{+}_{1/2} (paper's notation) and c mmin_m is m^{-}_{1/2} (paper's notation), see (A2) on p.370. c---------------------------------------------------------------- termp=(Mmin1-Mmin+Mplus-Mplus1)*(1.0d0/(mhalfr*mhalfr)-1.0d0) termp=termp*(pold_l-pold_r)/(pold_l/rold_l+pold_r/rold_r) c------------------------------------------------------------- c derivatives of the termp c------------------------------------------------------------- top=(Mmin1-Mmin+Mplus-Mplus1) bots=1.0d0/(pold_l/rold_l+pold_r/rold_r) bot=(pold_l-pold_r)*bots temph=1.0d0/(mhalfr*mhalfr)-1.0d0 c--------------------------- tmpr_l=(m1mr_l-dMmr_l+dMpr_l-m1pr_l)*bot*temph tmpr_l=tmpr_l-2.0d0*bot*top*dmfr_l/(mhalfr*mhalfr*mhalfr) tmpr_l=tmpr_l+temph*top*bot*bots*(pold_l/rold_l/rold_l) c--------------------------- tmpu_l=(m1mu_l-dMmu_l+dMpu_l-m1pu_l)*bot*temph tmpu_l=tmpu_l-2.0d0*bot*top*dmfu_l/(mhalfr*mhalfr*mhalfr) c--------------------------- tmpv_l=(m1mv_l-dMmv_l+dMpv_l-m1pv_l)*bot*temph tmpv_l=tmpv_l-2.0d0*bot*top*dmfv_l/(mhalfr*mhalfr*mhalfr) c--------------------------- tmpp_l=(m1mp_l-dMmp_l+dMpp_l-m1pp_l)*bot*temph tmpp_l=tmpp_l-2.0d0*bot*top*dmfp_l/(mhalfr*mhalfr*mhalfr) tmpp_l=tmpp_l+temph*top*bots*(1.0d0-bot/rold_l) c--------------------------- do 570 i=1,(nsp-1) tmpy_l.vv(i)=(m1my_l.vv(i)-dMmy_l.vv(i)+ & dMpy_l.vv(i)-m1py_l.vv(i))*bot*temph tmpy_l.vv(i)=tmpy_l.vv(i)-2.0d0*bot*top* & dmfg_l*dgdyl.vv(i)/(mhalfr*mhalfr*mhalfr) 570 continue c------------rrrrrrrr------- c------------rrrrrrrr------- tmpr_r=(m1mr_r-dMmr_r+dMpr_r-m1pr_r)*bot*temph tmpr_r=tmpr_r-2.0d0*bot*top*dmfr_r/(mhalfr*mhalfr*mhalfr) tmpr_r=tmpr_r+temph*top*bot*bots*(pold_r/rold_r/rold_r) c--------------------------- tmpu_r=(m1mu_r-dMmu_r+dMpu_r-m1pu_r)*bot*temph tmpu_r=tmpu_r-2.0d0*bot*top*dmfu_r/(mhalfr*mhalfr*mhalfr) c--------------------------- tmpv_r=(m1mv_r-dMmv_r+dMpv_r-m1pv_r)*bot*temph tmpv_r=tmpv_r-2.0d0*bot*top*dmfv_r/(mhalfr*mhalfr*mhalfr) c--------------------------- tmpp_r=(m1mp_r-dMmp_r+dMpp_r-m1pp_r)*bot*temph tmpp_r=tmpp_r-2.0d0*bot*top*dmfp_r/(mhalfr*mhalfr*mhalfr) tmpp_r=tmpp_r-temph*top*bots*(1.0d0+bot/rold_r) c--------------------------- do 571 i=1,(nsp-1) tmpy_r.vv(i)=(m1my_r.vv(i)-dMmy_r.vv(i)+ & dMpy_r.vv(i)-m1py_r.vv(i))*bot*temph tmpy_r.vv(i)=tmpy_r.vv(i)-2.0d0*bot*top* & dmfg_r*dgdyr.vv(i)/(mhalfr*mhalfr*mhalfr) 571 continue c------------------------------------------------------------- if(mmid .ge. 0.0d0) then mpl_m = mmid d2mr_l=dmir_l d2mu_l=dmiu_l d2mv_l=dmiv_l d2mp_l=dmip_l do 58 i=1,(nsp-1) d2my_l.vv(i)=dmiy_l.vv(i) 58 continue c------------ d2mr_r=dmir_r d2mu_r=dmiu_r d2mv_r=dmiv_r d2mp_r=dmip_r do 59 i=1,(nsp-1) d2my_r.vv(i)=dmiy_r.vv(i) 59 continue c------------ else mpl_m = 0.0d0 d2mr_l=0.0d0 d2mu_l=0.0d0 d2mv_l=0.0d0 d2mp_l=0.0d0 do 60 i=1,(nsp-1) d2my_l.vv(i)=0.0d0 60 continue c------------ d2mr_r=0.0d0 d2mu_r=0.0d0 d2mv_r=0.0d0 d2mp_r=0.0d0 do 61 i=1,(nsp-1) d2my_r.vv(i)=0.0d0 61 continue endif c--------------------------------------------- cc derivatives for the term termm cc------------------------------------------------------------------ if(mmid .le. 0.0d0) then mmin_m = mmid d3mr_l=dmir_l d3mu_l=dmiu_l d3mv_l=dmiv_l d3mp_l=dmip_l do 62 i=1,(nsp-1) d3my_l.vv(i)=dmiy_l.vv(i) 62 continue c------------ d3mr_r=dmir_r d3mu_r=dmiu_r d3mv_r=dmiv_r d3mp_r=dmip_r do 63 i=1,(nsp-1) d3my_r.vv(i)=dmiy_r.vv(i) 63 continue c------------ else mmin_m = 0.0d0 d3mr_l=0.0d0 d3mu_l=0.0d0 d3mv_l=0.0d0 d3mp_l=0.0d0 do 64 i=1,(nsp-1) d3my_l.vv(i)=0.0d0 64 continue c------------ d3mr_r=0.0d0 d3mu_r=0.0d0 d3mv_r=0.0d0 d3mp_r=0.0d0 do 65 i=1,(nsp-1) d3my_r.vv(i)=0.0d0 65 continue endif c--------------------------------------------------------------- c Computing the calligraphic P+ and P- with their derivatives, c see (21a) & (21b) on p.368. c--------------------------------------------------------------- if(ml .ge. 1.0d0) then Pplus = 1.0d0 else if((ml .gt. -1.0d0) .and. (ml .lt. 1.0d0)) then Pplus=(ml+1.0d0)*(ml+1.0d0)*(2.0d0-ml)/4.0d0 Pplus=Pplus+alpha*ml*(ml*ml-1.0d0)*(ml*ml-1.0d0) else Pplus = 0.0d0 endif endif c--------------------------------------------------------------- if(mr .ge. 1.0d0) then Pmin = 0.0d0 else if((mr .gt. -1.0d0) .and. (mr .lt. 1.0d0)) then Pmin=(mr-1.0d0)*(mr-1.0d0)*(2.0d0+mr)/4.0d0 Pmin=Pmin-alpha*mr*(mr*mr-1.0d0)*(mr*mr-1.0d0) else Pmin = 1.0d0 endif endif c--------------------------------------------------------------- brac_l=(ml+1.0d0)*(2.0d0-ml)/2.0d0-(ml+1.0d0)*(ml+1.0d0)/4.0d0 brac_l=brac_l+alpha*(ml*ml-1.0d0)*(ml*ml-1.0d0) brac_l=brac_l+4.0d0*alpha*ml*ml*(ml*ml-1.0d0) c-------------- brac_r=(mr-1.0d0)*(2.0d0+mr)/2.0d0+(mr-1.0d0)*(mr-1.0d0)/4.0d0 brac_r=brac_r-alpha*(mr*mr-1.0d0)*(mr*mr-1.0d0) brac_r=brac_r-4.0d0*alpha*mr*mr*(mr*mr-1.0d0) c--------------------------------------------------------------- if((ml .gt. -1.0d0) .and. (ml .lt. 1.0d0)) then dPpr_l=brac_l*dmlr_l dPpr_r=brac_l*dmlr_r c------------ dPpu_l=brac_l*dmlu_l dPpu_r=brac_l*dmlu_r c------------ dPpv_l=brac_l*dmlv_l dPpv_r=brac_l*dmlv_r c------------ dPpp_l=brac_l*dmlp_l dPpp_r=brac_l*dmlp_r c------------ do 66 i=1,(nsp-1) dPpy_l.vv(i)=brac_l*dmly_l.vv(i) dPpy_r.vv(i)=brac_l*dmly_r.vv(i) 66 continue c------------ else dPpr_l=0.0d0 dPpr_r=0.0d0 c----------- dPpu_l=0.0d0 dPpu_r=0.0d0 c----------- dPpv_l=0.0d0 dPpv_r=0.0d0 c----------- dPpp_l=0.0d0 dPpp_r=0.0d0 c----------- do 67 i=1,(nsp-1) dPpy_l.vv(i)=0.0d0 dPpy_r.vv(i)=0.0d0 67 continue c----------- endif c--------------------------------------------------------------- if((mr .gt. -1.0d0) .and. (mr .lt. 1.0d0)) then dPmr_l=brac_r*dmrr_l dPmr_r=brac_r*dmrr_r c------------ dPmu_l=brac_r*dmru_l dPmu_r=brac_r*dmru_r c------------ dPmv_l=brac_r*dmrv_l dPmv_r=brac_r*dmrv_r c------------ dPmp_l=brac_r*dmrp_l dPmp_r=brac_r*dmrp_r c------------ do 68 i=1,(nsp-1) dPmy_l.vv(i)=brac_r*dmry_l.vv(i) dPmy_r.vv(i)=brac_r*dmry_r.vv(i) 68 continue c------------ else dPmr_l=0.0d0 dPmr_r=0.0d0 c----------- dPmu_l=0.0d0 dPmu_r=0.0d0 c----------- dPmv_l=0.0d0 dPmv_r=0.0d0 c----------- dPmp_l=0.0d0 dPmp_r=0.0d0 c----------- do 69 i=1,(nsp-1) dPmy_l.vv(i)=0.0d0 dPmy_r.vv(i)=0.0d0 69 continue c----------- endif c------------------------------------------------------------------- c computing pmid -- p_{1/2} and its derivatives, see (20b), p.367. c------------------------------------------------------------------- dpir_l=dPpr_l*pold_l+dPmr_l*pold_r dpiu_l=dPpu_l*pold_l+dPmu_l*pold_r dpiv_l=dPpv_l*pold_l+dPmv_l*pold_r dpip_l=dPpp_l*pold_l+Pplus+dPmp_l*pold_r do 70 i=1,(nsp-1) dpiy_l.vv(i)=dPpy_l.vv(i)*pold_l+dPmy_l.vv(i)*pold_r 70 continue c---------------------------- dpir_r=dPpr_r*pold_l+dPmr_r*pold_r dpiu_r=dPpu_r*pold_l+dPmu_r*pold_r dpiv_r=dPpv_r*pold_l+dPmv_r*pold_r dpip_r=dPpp_r*pold_l+Pmin+dPmp_r*pold_r do 71 i=1,(nsp-1) dpiy_r.vv(i)=dPpy_r.vv(i)*pold_l+dPmy_r.vv(i)*pold_r 71 continue c--------------------------------------------------------------------- c Computation of the mass flux (rho * u)_1/2 c--------------------------------------------------------------- rum=am*(mpl_m*rold_l+mmin_m*rold_r)+canc*am*termp c------------------------------------------------------- rumr_l=damr_l*(mpl_m*rold_l+mmin_m*rold_r)+ & am*(d2mr_l*rold_l+mpl_m+d3mr_l*rold_r) rumr_l=rumr_l+canc*coef*(damr_l*termp+am*tmpr_l) rumu_l=damu_l*(mpl_m*rold_l+mmin_m*rold_r)+ & am*(d2mu_l*rold_l+d3mu_l*rold_r) rumu_l=rumu_l+canc*coef*(damu_l*termp+am*tmpu_l) rumv_l=damv_l*(mpl_m*rold_l+mmin_m*rold_r)+ & am*(d2mv_l*rold_l+d3mv_l*rold_r) rumv_l=rumv_l+canc*coef*(damv_l*termp+am*tmpv_l) rump_l=damp_l*(mpl_m*rold_l+mmin_m*rold_r)+ & am*(d2mp_l*rold_l+d3mp_l*rold_r) rump_l=rump_l+canc*coef*(damp_l*termp+am*tmpp_l) c------------------ do 710 i=1,(nsp-1) rumy_l.vv(i)=damg_l*dgdyl.vv(i)*(mpl_m*rold_l+ & mmin_m*rold_r)+am*(d2my_l.vv(i)*rold_l+d3my_l.vv(i)*rold_r) rumy_l.vv(i)=rumy_l.vv(i)+canc*coef* & (damg_l*dgdyl.vv(i)*termp+am*tmpy_l.vv(i)) 710 continue c------------------------------------------------- rumr_r=damr_r*(mpl_m*rold_l+mmin_m*rold_r)+ & am*(d2mr_r*rold_l+mmin_m+d3mr_r*rold_r) rumr_r=rumr_r+canc*coef*(damr_r*termp+am*tmpr_r) rumu_r=damu_r*(mpl_m*rold_l+mmin_m*rold_r)+ & am*(d2mu_r*rold_l+d3mu_r*rold_r) rumu_r=rumu_r+canc*coef*(damu_r*termp+am*tmpu_r) rumv_r=damv_r*(mpl_m*rold_l+mmin_m*rold_r)+ & am*(d2mv_r*rold_l+d3mv_r*rold_r) rumv_r=rumv_r+canc*coef*(damv_r*termp+am*tmpv_r) rump_r=damp_r*(mpl_m*rold_l+mmin_m*rold_r)+ & am*(d2mp_r*rold_l+d3mp_r*rold_r) rump_r=rump_r+canc*coef*(damp_r*termp+am*tmpp_r) do 711 i=1,(nsp-1) rumy_r.vv(i)=damg_r*dgdyr.vv(i)*(mpl_m*rold_l+mmin_m*rold_r)+ & am*(d2my_r.vv(i)*rold_l+d3my_r.vv(i)*rold_r) rumy_r.vv(i)=rumy_r.vv(i)+canc*coef* & (damg_r*dgdyr.vv(i)*termp+am*tmpy_r.vv(i)) 711 continue c--------------------------------------------------------------------- c computing JACOBIAN as a derivative of the numerical flux function c with respect to the primitive variables. c Notation: jl(i,j) --- is the derivative of the i-component of the c flux function with respect to the j-component of the c vector of primitive variables of the left state. c jr(i,j) --- is the derivative of the i-component of the c flux function with respect to the j-component of the c vector of primitive variables of the right state. c--------------------------------------------------------------------- jl.jac(1,1)=0.0d0 jl.jac(1,2)=0.0d0 jl.jac(1,3)=0.0d0 jl.jac(1,4)=0.0d0 do 72 i=1,(nsp-1) jl.jac(1,4+i)=0.0d0 72 continue c------------------------------------ jr.jac(1,1)=0.0d0 jr.jac(1,2)=0.0d0 jr.jac(1,3)=0.0d0 jr.jac(1,4)=0.0d0 do 73 i=1,(nsp-1) jr.jac(1,4+i)=0.0d0 73 continue c------------------------------------ c--------------------------------------------------------- if(rum .ge. 0.0d0) then br3=rumr_l*un_l br4=rumr_l*ut_l else br3=rumr_l*un_r br4=rumr_l*ut_r endif jl.jac(2,1)=n_x*(br3+dpir_l)+br4*t_x jl.jac(3,1)=n_y*(br3+dpir_l)+br4*t_y c------------------- if(rum .ge. 0.0d0) then br3=rumu_l*un_l+rum*n_x br4=rumu_l*ut_l+rum*t_x else br3=rumu_l*un_r br4=rumu_l*ut_r endif jl.jac(2,2)=n_x*(br3+dpiu_l)+br4*t_x jl.jac(3,2)=n_y*(br3+dpiu_l)+br4*t_y c------------------- if(rum .ge. 0.0d0) then br3=rumv_l*un_l+rum*n_y br4=rumv_l*ut_l+rum*t_y else br3=rumv_l*un_r br4=rumv_l*ut_r endif jl.jac(2,3)=n_x*(br3+dpiv_l)+br4*t_x jl.jac(3,3)=n_y*(br3+dpiv_l)+br4*t_y c------------------- if(rum .ge. 0.0d0) then br3=rump_l*un_l br4=rump_l*ut_l else br3=rump_l*un_r br4=rump_l*ut_r endif jl.jac(2,4)=n_x*(br3+dpip_l)+br4*t_x jl.jac(3,4)=n_y*(br3+dpip_l)+br4*t_y c------------------------------------------------------------- do 74 i=1,(nsp-1) if(rum .ge. 0.0d0) then br3=rumy_l.vv(i)*un_l br4=rumy_l.vv(i)*ut_l else br3=rumy_l.vv(i)*un_r br4=rumy_l.vv(i)*ut_r endif jl.jac(2,4+i)=n_x*(br3+dpiy_l.vv(i))+br4*t_x jl.jac(3,4+i)=n_y*(br3+dpiy_l.vv(i))+br4*t_y 74 continue c------------------------------------------------------------- if(rum .ge. 0.0d0) then br3=rumr_r*un_l br4=rumr_r*ut_l else br3=rumr_r*un_r br4=rumr_r*ut_r endif jr.jac(2,1)=n_x*(br3+dpir_r)+br4*t_x jr.jac(3,1)=n_y*(br3+dpir_r)+br4*t_y c------------------- if(rum .ge. 0.0d0) then br3=rumu_r*un_l br4=rumu_r*ut_l else br3=rumu_r*un_r+rum*n_x br4=rumu_r*ut_r+rum*t_x endif jr.jac(2,2)=n_x*(br3+dpiu_r)+br4*t_x jr.jac(3,2)=n_y*(br3+dpiu_r)+br4*t_y c------------------- if(rum .ge. 0.0d0) then br3=rumv_r*un_l br4=rumv_r*ut_l else br3=rumv_r*un_r+rum*n_y br4=rumv_r*ut_r+rum*t_y endif jr.jac(2,3)=n_x*(br3+dpiv_r)+br4*t_x jr.jac(3,3)=n_y*(br3+dpiv_r)+br4*t_y c------------------- if(rum .ge. 0.0d0) then br3=rump_r*un_l br4=rump_r*ut_l else br3=rump_r*un_r br4=rump_r*ut_r endif jr.jac(2,4)=n_x*(br3+dpip_r)+br4*t_x jr.jac(3,4)=n_y*(br3+dpip_r)+br4*t_y c-------------------- do 75 i=1,(nsp-1) if(rum .ge. 0.0d0) then br3=rumy_r.vv(i)*un_l br4=rumy_r.vv(i)*ut_l else br3=rumy_r.vv(i)*un_r br4=rumy_r.vv(i)*ut_r endif jr.jac(2,4+i)=n_x*(br3+dpiy_r.vv(i))+br4*t_x jr.jac(3,4+i)=n_y*(br3+dpiy_r.vv(i))+br4*t_y 75 continue c------------------------------------------------------------- c ------ f44444444444444444444444444444 --------- c------------------------------------------------------------- jl.jac(4,1)=0.0d0 jl.jac(4,2)=0.0d0 jl.jac(4,3)=0.0d0 jl.jac(4,4)=0.0d0 c---------------------------------------------------------- jr.jac(4,1)=0.0d0 jr.jac(4,2)=0.0d0 jr.jac(4,3)=0.0d0 jr.jac(4,4)=0.0d0 c--------------------- do 76 i=1,(nsp-1) jl.jac(4,4+i)=0.0d0 jr.jac(4,4+i)=0.0d0 76 continue c------------------------------------------------------------- c------------------------------------------------------------- do 78 i=1,(nsp-1) jl.jac(4+i,1)=0.0d0 jr.jac(4+i,1)=0.0d0 c-------------------------- jl.jac(4+i,2)=0.0d0 jr.jac(4+i,2)=0.0d0 c-------------------------- jl.jac(4+i,3)=0.0d0 jr.jac(4+i,3)=0.0d0 c-------------------------- jl.jac(4+i,4)=0.0d0 jr.jac(4+i,4)=0.0d0 c-------------------------- do 780 j=1,(nsp-1) jl.jac(4+i,4+j)=0.0d0 jr.jac(4+i,4+j)=0.0d0 780 continue 78 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)=-uold_l/rold_l wl.jac(2,2)=1.0d0/rold_l 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)=-vold_l/rold_l wl.jac(3,2)=0.0d0 wl.jac(3,3)=1.0d0/rold_l 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+dgdyl.vv(i)*yl.yet(i) 86 continue br1=br1*pold_l/(rold_l*gm1l) wl.jac(4,1)=gm1l*(uold_l*uold_l+vold_l*vold_l)/2.0d0-br1 wl.jac(4,2)=-uold_l*gm1l wl.jac(4,3)=-vold_l*gm1l wl.jac(4,4)=gm1l do 87 i=1,(nsp-1) wl.jac(4,4+i)=dgdyl.vv(i)*pold_l/(rold_l*gm1l) 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)=-yl.yet(i)/rold_l 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/rold_l endif 890 continue 88 continue c------------------------------ c------------------------------ wr.jac(1,1)=1.0d0 wr.jac(1,2)=0.0d0 wr.jac(1,3)=0.0d0 wr.jac(1,4)=0.0d0 do 90 i=1,(nsp-1) wr.jac(1,4+i)=0.0d0 90 continue c------------------------------ wr.jac(2,1)=-uold_r/rold_r wr.jac(2,2)=1.0d0/rold_r wr.jac(2,3)=0.0d0 wr.jac(2,4)=0.0d0 do 91 i=1,(nsp-1) wr.jac(2,4+i)=0.0d0 91 continue c------------------------------ wr.jac(3,1)=-vold_r/rold_r wr.jac(3,2)=0.0d0 wr.jac(3,3)=1.0d0/rold_r wr.jac(3,4)=0.0d0 do 92 i=1,(nsp-1) wr.jac(3,4+i)=0.0d0 92 continue c------------------------------ br1=0.0d0 do 860 i=1,(nsp-1) br1=br1+dgdyr.vv(i)*yr.yet(i) 860 continue br1=br1*pold_r/(rold_r*gm1r) wr.jac(4,1)=gm1r*(uold_r*uold_r+vold_r*vold_r)/2.0d0-br1 wr.jac(4,2)=-uold_r*gm1r wr.jac(4,3)=-vold_r*gm1r wr.jac(4,4)=gm1r do 93 i=1,(nsp-1) wr.jac(4,4+i)=dgdyr.vv(i)*pold_r/(rold_r*gm1r) 93 continue c------------------------------ do 94 i=1,(nsp-1) do 95 j=1,4 wr.jac(4+i,j)=0.0d0 if(j.eq.1) wr.jac(4+i,j)=-yr.yet(i)/rold_r 95 continue c---------------------- do 950 j=5,(3+nsp) wr.jac(4+i,j)=0.0d0 if(4+i.eq.j) wr.jac(4+i,j)=1.0d0/rold_r 950 continue 94 continue c---------------------------------------------- c---------------------------------------------- do 1 i=1,(3+nsp) do 2 j=1,(3+nsp) jtl.jac(i,j)=0.0d0 jtr.jac(i,j)=0.0d0 do 3 k=1,(3+nsp) jtl.jac(i,j)=jtl.jac(i,j)+jl.jac(i,k)*wl.jac(k,j) jtr.jac(i,j)=jtr.jac(i,j)+jr.jac(i,k)*wr.jac(k,j) 3 continue 2 continue 1 continue c---------------------------------------------------------------------- tcoef=nvect(1)*tvect(2)+tvect(1)*nvect(2) bcoef=nvect(1)*tvect(2)-tvect(1)*nvect(2) c---------------------------------------------------------------------- do 11 i=1,(3+nsp) jtl.jac(i,1)=jtl.jac(i,1)+jtr.jac(i,1) jtl.jac(i,2)=jtl.jac(i,2)+jtr.jac(i,2)*(-tcoef/bcoef)+ & jtr.jac(i,3)*2.0d0*nvect(1)*tvect(1)/bcoef jtl.jac(i,3)=jtl.jac(i,3)+jtr.jac(i,2)* & (-2.0d0*nvect(2)*tvect(2)/bcoef)+ & jtr.jac(i,3)*tcoef/bcoef jtl.jac(i,4)=jtl.jac(i,4)+jtr.jac(i,4) do 777 j=5,(3+nsp) jtl.jac(i,j)=jtl.jac(i,j)+jtr.jac(i,j) 777 continue 11 continue c---------------------------------- c do 11 i=1,(3+nsp) c jtl.jac(i,1)=jtl.jac(i,1)+jtr.jac(i,1) c jtl.jac(i,2)=jtl.jac(i,2)+jtr.jac(i,2)*(-tcoef/bcoef)+ c & jtr.jac(i,3)*2.0d0*nvect(1)*tvect(1)/bcoef c jtl.jac(i,3)=jtl.jac(i,3)+jtr.jac(i,2)* c & (-2.0d0*nvect(2)*tvect(2)/bcoef)+ c & jtr.jac(i,3)*tcoef/bcoef c jtl.jac(i,4)=jl.jac(i,4)+jr.jac(i,4) c do 777 j=5,(3+nsp) c jl.jac(i,j)=jl.jac(i,j)+jr.jac(i,j) c 777 continue c 11 continue c------------------------------------------------------------------- SEGSUP DMLY_L, DMLY_R, & dmry_l, dmry_r, & dMpy_l, dMpy_r, & dMmy_l, dMmy_r, & dmiy_l, dmiy_r, & d3my_l, d3my_r, & d2my_l, d2my_r, & dPpy_l, dPpy_r, & dPmy_l, dPmy_r, & dpiy_l, dpiy_r, & dgdyl, dgdyr C-------------------------------------------- SEGSUP f C-------------------------------------------- c jll = jl SEGDES JL SEGDES JR SEGSUP WL SEGSUP WR C-------------------------------------------- jll=jtl SEGDES JTL SEGDES JTR SEGDES YL SEGDES YR SEGDES CP SEGDES CV SEGDES MLRECP, MLRECV c------------------- return end