Télécharger jbmmw2.eso

Retour à la liste

jbmmw2
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