Télécharger la plaquette Cast3M


fr en

Télécharger jakbm.eso

Retour à la liste

Numérotation des lignes : 

jakbm

  1. C JAKBM     SOURCE    CHAT      05/01/13    00:48:48     5004
  2.        SUBROUTINE JAKBM(JTL,JTR,WVEC_L,WVEC_R,NVECT,TVECT,
  3.      &                   ga,v_inf)
  4. C************************************************************************
  5. C
  6. C PROJET            :  CASTEM 2000
  7. C
  8. C NOM               :  JACBM
  9. C
  10. C DESCRIPTION       :  Voir KONJA6
  11. C
  12. C LANGAGE           :  FORTRAN 77 + ESOPE 2000 (avec estensions CISI)
  13. C
  14. C AUTEUR            :  S. KUDRIAKOV, DM2S/SFME/LTMF
  15. C
  16. C************************************************************************
  17. C
  18. c----------------------------------------------------------------------
  19. c GENERAL DESCRIPTION:
  20. c   This subroutine provides the jacobians which are the derivatives
  21. c   of the numerical flux function defined at the cell interface
  22. c   with respect to the conservative variables of the left and right
  23. c   cells (relative to the cell interface).
  24. c   The low-mach number corrections are made for the flux functions
  25. c
  26. c EQUATIONS: 2D Euler equations of gas dynamics
  27. c
  28. c
  29. c REFERENCE: 1) JCP, 129, 364-382 (1996)
  30. c               " A Sequel to AUSM: AUSM+ ".
  31. c                     M.S.Liou
  32. c            2) AIAA Journal, Sept. 1998
  33. c          " Low-Diffusion Flux-Splitting Methods for Flows at All Speeds"
  34. c                J.R.Edwards and M.S.Liou
  35. c----------------------------------------------------------------------
  36. c INPUT:
  37. c
  38. c      alpha   -- parameter of the AUSM+ scheme in the Pressure function;
  39. c                 ( -3/4 <= alpha <= 3/16 ) (imposed as a parameter)
  40. c
  41. c      beta    -- parameter of the AUSM+ scheme in the Mach function;
  42. c                 ( -1/16 <= beta <= 1/2 )  (imposed as a parameter)
  43. c
  44. c      wvec_l  -- vector of the primitive variables (rho,ux,uy,p) at the
  45. c                 left cell;
  46. c
  47. c      wvec_r  -- vector of the primitive variables (rho,ux,uy,p) at the
  48. c                 right cell;
  49. c
  50. c      nvect   -- normal vector to the interface (2 components in 2D);
  51. c
  52. c      tvect   -- tangential vector to the interface;
  53. c
  54. c      ga      -- ratio of the specific heats (assumed constant)
  55. c
  56. c      v_inf   -- parameter for reference velocity, see Ref. 2).
  57. c----------------------------------------------------------------------
  58. c
  59. c OUTPUT:
  60. c
  61. c      jtl     -- jakobian matrix 4 by 4 - derivatives of the numerical
  62. c                 flux function with respect to the conservative variables
  63. c                 from the left cell;
  64. c
  65. c      jtr     -- jakobian matrix 4 by 4 - derivatives of the numerical
  66. c                 flux function with respect to the conservative variables
  67. c                 from the right cell.
  68. c----------------------------------------------------------------------
  69.       IMPLICIT INTEGER(I-N)
  70.        real*8 wvec_l(4),wvec_r(4)
  71.        real*8 nvect(2),tvect(2)
  72.        real*8 jl(4,4),jr(4,4)
  73.        real*8 wl(4,4),wr(4,4)
  74.        real*8 jtl(4,4),jtr(4,4)
  75.        real*8 alpha,beta
  76.        real*8 ga,gm1,temph
  77.        real*8 n_x,n_y,t_x,t_y
  78.        real*8 un_l, un_r, ut_l, ut_r
  79.        real*8 ml,mr,Mplus,Mmin,mmid
  80.        real*8 mpl_m, mmin_m,am
  81.        real*8 rold_l,uold_l,vold_l,pold_l,eold_l
  82.        real*8 rold_r,uold_r,vold_r,pold_r,eold_r
  83.        real*8 Pplus,Pmin
  84.        real*8 hr_l,hr_r,top,bot,bots
  85.        real*8 br1,br3,br4,temp_l,temp_r,brac_l,brac_r
  86.        real*8 aleft, arigh
  87.        real*8 damr_l,damr_r,damu_l,damu_r
  88.        real*8 damv_l,damv_r,damp_l,damp_r
  89.        real*8 dmlr_l,dmlr_r,dmlu_l,dmlu_r
  90.        real*8 dmlv_l,dmlv_r,dmlp_l,dmlp_r
  91.        real*8 dmrr_l,dmrr_r,dmru_l,dmru_r
  92.        real*8 dmrv_l,dmrv_r,dmrp_l,dmrp_r
  93.        real*8 dMpr_l,dMpr_r,dMpu_l,dMpu_r
  94.        real*8 dMpv_l,dMpv_r,dMpp_l,dMpp_r
  95.        real*8 dMmr_l,dMmr_r,dMmu_l,dMmu_r
  96.        real*8 dMmv_l,dMmv_r,dMmp_l,dMmp_r
  97.        real*8 dmir_l,dmir_r,dmiu_l,dmiu_r
  98.        real*8 dmiv_l,dmiv_r,dmip_l,dmip_r
  99.        real*8 d3mr_l,d3mr_r,d3mu_l,d3mu_r
  100.        real*8 d3mv_l,d3mv_r,d3mp_l,d3mp_r
  101.        real*8 d2mr_l,d2mr_r,d2mu_l,d2mu_r
  102.        real*8 d2mv_l,d2mv_r,d2mp_l,d2mp_r
  103.        real*8 dPpr_l,dPpr_r,dPpu_l,dPpu_r
  104.        real*8 dPpv_l,dPpv_r,dPpp_l,dPpp_r
  105.        real*8 dPmr_l,dPmr_r,dPmu_l,dPmu_r
  106.        real*8 dPmv_l,dPmv_r,dPmp_l,dPmp_r
  107.        real*8 dpir_l,dpir_r,dpiu_l,dpiu_r
  108.        real*8 dpiv_l,dpiv_r,dpip_l,dpip_r
  109.        real*8 epsil,qq,amw,Mmin1,Mplus1
  110.        real*8 fmid,mlw,mrw,termp
  111.        real*8 ur_r,ur_l,urm,mhalf,mhalfr
  112.        real*8 durr_l,durr_r,duru_l,duru_r
  113.        real*8 durv_l,durv_r,durp_l,durp_r
  114.        real*8 dmhr_l,dmhr_r,dmhu_l,dmhu_r
  115.        real*8 dmhv_l,dmhv_r,dmhp_l,dmhp_r
  116.        real*8 dmfr_l,dmfr_r,dmfu_l,dmfu_r
  117.        real*8 dmfv_l,dmfv_r,dmfp_l,dmfp_r
  118.        real*8 dfm_mf,dfm_mh
  119.        real*8 dfmr_l,dfmr_r,dfmu_l,dfmu_r
  120.        real*8 dfmv_l,dfmv_r,dfmp_l,dfmp_r
  121.        real*8 m1mr_l,m1mr_r,m1mu_l,m1mu_r
  122.        real*8 m1mv_l,m1mv_r,m1mp_l,m1mp_r
  123.        real*8 m1pr_l,m1pr_r,m1pu_l,m1pu_r
  124.        real*8 m1pv_l,m1pv_r,m1pp_l,m1pp_r
  125.        real*8 tmpr_l,tmpr_r,tmpu_l,tmpu_r
  126.        real*8 tmpv_l,tmpv_r,tmpp_l,tmpp_r
  127.        real*8 coef,canc,v_inf
  128.        real*8 sr_l,sr_r,su_l,su_r,sv_l,sv_r,sp_l,sp_r
  129.        real*8 rum,rumr_l,rumu_l,rumv_l,rump_l
  130.        real*8 rumr_r,rumu_r,rumv_r,rump_r
  131.        integer i,j,k
  132.        parameter(alpha = 0.1875D0, beta = 0.125D0)
  133.        parameter(epsil = 1.0d0)
  134.        canc=1.0d0
  135. c-------------------------------------------------------------
  136.        n_x=nvect(1)
  137.        n_y=nvect(2)
  138.        t_x=tvect(1)
  139.        t_y=tvect(2)
  140. c----------------------------
  141.        gm1=ga-1.0d0
  142. c----------------------------
  143.        rold_l=wvec_l(1)
  144.        uold_l=wvec_l(2)
  145.        vold_l=wvec_l(3)
  146.        pold_l=wvec_l(4)
  147. c-----------------------
  148.        rold_r=wvec_r(1)
  149.        uold_r=wvec_r(2)
  150.        vold_r=wvec_r(3)
  151.        pold_r=wvec_r(4)
  152. c------------------------------------------------------------------
  153. c Computation of the specific total energy on the left and right.
  154. c------------------------------------------------------------------
  155.        eold_l=(uold_l*uold_l+vold_l*vold_l)/2.0d0
  156.        eold_l=eold_l+pold_l/(gm1*rold_l)
  157.        eold_r=(uold_r*uold_r+vold_r*vold_r)/2.0d0
  158.        eold_r=eold_r+pold_r/(gm1*rold_r)
  159. c------------------------------------------------------------------
  160. c   Computing reference velocity and its derivatives
  161. c------------------------------------------------------------------
  162.        aleft=sqrt(ga*pold_l/rold_l)
  163.        arigh=sqrt(ga*pold_r/rold_r)
  164.        qq=sqrt(uold_l*uold_l+vold_l*vold_l)
  165.        if(qq .lt. (epsil*v_inf)) then
  166.           ur_l = epsil*v_inf
  167.           durr_l=0.0d0
  168.           duru_l=0.0d0
  169.           durv_l=0.0d0
  170.           durp_l=0.0d0
  171.        else
  172.           ur_l=qq
  173.           durr_l=0.0d0
  174.           duru_l=uold_l/qq
  175.           durv_l=vold_l/qq
  176.           durp_l=0.0d0
  177.        endif
  178. c------------------------------
  179.        if(ur_l .ge. aleft) then
  180.           ur_l=aleft
  181.           durr_l=-aleft/(2.0d0*rold_l)
  182.           duru_l=0.0d0
  183.           durv_l=0.0d0
  184.           durp_l=ga/(2.0d0*aleft*rold_l)
  185.        endif
  186. c------------------------------------------------------------------
  187.        qq=sqrt(uold_r*uold_r+vold_r*vold_r)
  188.        if(qq .lt. (epsil*v_inf)) then
  189.           ur_r = epsil*v_inf
  190.           durr_r=0.0d0
  191.           duru_r=0.0d0
  192.           durv_r=0.0d0
  193.           durp_r=0.0d0
  194.        else
  195.           ur_r=qq
  196.           durr_r=0.0d0
  197.           duru_r=uold_r/qq
  198.           durv_r=vold_r/qq
  199.           durp_r=0.0d0
  200.        endif
  201. c---------------------------
  202.        if(ur_r .ge. arigh) then
  203.           ur_r=arigh
  204.           durr_r=-arigh/(2.0d0*rold_r)
  205.           duru_r=0.0d0
  206.           durv_r=0.0d0
  207.           durp_r=ga/(2.0d0*arigh*rold_r)
  208.        endif
  209. c------------------------------------------------------------------
  210. c-------------------------------------------------------------------
  211. c  Reference velocity at the interface is taken as an average
  212. c  of the reference velocities of the neighbouring cells
  213. c-------------------------------------------------------------------
  214.        urm=0.5d0*(ur_l+ur_r)
  215.        durr_l=0.5d0*durr_l
  216.        duru_l=0.5d0*duru_l
  217.        durv_l=0.5d0*durv_l
  218.        durp_l=0.5d0*durp_l
  219. c-------------------------
  220.        durr_r=0.5d0*durr_r
  221.        duru_r=0.5d0*duru_r
  222.        durv_r=0.5d0*durv_r
  223.        durp_r=0.5d0*durp_r
  224. c-------------------------------------------------------------------
  225. c Computation of the speed of sound and its derivatives;
  226. c numerical speed of sound at the interface is taken as an average
  227. c of the speeds of sounds of the neighbouring cells
  228. c-------------------------------------------------------------------
  229.        am=0.5d0*(aleft+arigh)
  230. c-------------------------------------------------------------------
  231.        if(abs(urm/am-1.0d0).le.0.000001d0) then
  232.          coef=0.0d0
  233.        else
  234.          coef=1.0d0
  235.        endif
  236. c-------------------------------------------------------------------
  237.        damr_r=-arigh/(4.0d0*rold_r)
  238.        damu_r=0.0d0
  239.        damv_r=0.0d0
  240.        damp_r=ga/(4.0d0*arigh*rold_r)
  241. c-----------------------
  242.        damr_l=-aleft/(4.0d0*rold_l)
  243.        damu_l=0.0d0
  244.        damv_l=0.0d0
  245.        damp_l=ga/(4.0d0*aleft*rold_l)
  246. c-------------------------------------------------------------------
  247. c Computing numerical Mach number and its derivatives,
  248. c see p.370, under (A1).
  249. c-------------------------------------------------------------------
  250.        un_l=uold_l*n_x+vold_l*n_y
  251.        un_r=uold_r*n_x+vold_r*n_y
  252.        ut_l=uold_l*t_x+vold_l*t_y
  253.        ut_r=uold_r*t_x+vold_r*t_y
  254. c-------------------------------------------------------------------
  255.        ml=un_l/am
  256.        mr=un_r/am
  257.        mhalf=0.5d0*(un_l+un_r)/am
  258. c---------------------
  259.        top=0.5d0*(un_l+un_r)/(am*am)
  260.        dmhr_l=-top*damr_l
  261.        dmhu_l=n_x/2.0d0/am-top*damu_l
  262.        dmhv_l=n_y/2.0d0/am-top*damv_l
  263.        dmhp_l=-top*damp_l
  264. c---------------------
  265.        dmhr_r=-top*damr_r
  266.        dmhu_r=n_x/2.0d0/am-top*damu_r
  267.        dmhv_r=n_y/2.0d0/am-top*damv_r
  268.        dmhp_r=-top*damp_r
  269. c--------------------------------
  270.        mhalfr=urm/am
  271. c---------------------
  272.        top=urm/(am*am)
  273.        dmfr_l=durr_l/am-top*damr_l
  274.        dmfu_l=duru_l/am-top*damu_l
  275.        dmfv_l=durv_l/am-top*damv_l
  276.        dmfp_l=durp_l/am-top*damp_l
  277. c---------------------
  278.        dmfr_r=durr_r/am-top*damr_r
  279.        dmfu_r=duru_r/am-top*damu_r
  280.        dmfv_r=durv_r/am-top*damv_r
  281.        dmfp_r=durp_r/am-top*damp_r
  282. c-------------------------------------------------------------------
  283. c Scaling function for the speed of sound and its derivatives
  284. c-------------------------------------------------------------------
  285.        top=(1.0d0-mhalfr*mhalfr)*(1.0d0-mhalfr*mhalfr)
  286.        top=top*mhalf*mhalf+4.0d0*mhalfr*mhalfr
  287.        bot=1.0d0+mhalfr*mhalfr
  288.        if(abs(canc-0.0d0).lt.0.000001d0) then
  289.          fmid=1.0d0
  290.        else
  291.          fmid=sqrt(top)/bot
  292.        endif
  293. c--------------------------
  294.        temph=-4.0d0*(1.0d0-mhalfr*mhalfr)*mhalfr
  295.        temph=temph*mhalf*mhalf+8.0d0*mhalfr
  296.        dfm_mf=temph/(sqrt(top)*2.0d0*bot)
  297.        dfm_mf=dfm_mf-sqrt(top)*2.0d0*mhalfr/(bot*bot)
  298. c--------------------------
  299.        temph=2.0d0*(1.0d0-mhalfr*mhalfr)*(1.0d0-mhalfr*mhalfr)*mhalf
  300.        dfm_mh=temph/(2.0d0*bot*sqrt(top))
  301. c--------------------------
  302.        dfmr_l=dfm_mf*dmfr_l+dfm_mh*dmhr_l
  303.        dfmu_l=dfm_mf*dmfu_l+dfm_mh*dmhu_l
  304.        dfmv_l=dfm_mf*dmfv_l+dfm_mh*dmhv_l
  305.        dfmp_l=dfm_mf*dmfp_l+dfm_mh*dmhp_l
  306. c--------------------------
  307.        dfmr_r=dfm_mf*dmfr_r+dfm_mh*dmhr_r
  308.        dfmu_r=dfm_mf*dmfu_r+dfm_mh*dmhu_r
  309.        dfmv_r=dfm_mf*dmfv_r+dfm_mh*dmhv_r
  310.        dfmp_r=dfm_mf*dmfp_r+dfm_mh*dmhp_r
  311. c--------------------------
  312.        amw=fmid*am
  313.        mlw=un_l/amw
  314.        mrw=un_r/amw
  315. c--------------------------
  316.        damr_l=canc*coef*dfmr_l*am+fmid*damr_l
  317.        damu_l=canc*coef*dfmu_l*am+fmid*damu_l
  318.        damv_l=canc*coef*dfmv_l*am+fmid*damv_l
  319.        damp_l=canc*coef*dfmp_l*am+fmid*damp_l
  320. c--------------------------
  321.        damr_r=canc*coef*dfmr_r*am+fmid*damr_r
  322.        damu_r=canc*coef*dfmu_r*am+fmid*damu_r
  323.        damv_r=canc*coef*dfmv_r*am+fmid*damv_r
  324.        damp_r=canc*coef*dfmp_r*am+fmid*damp_r
  325. c-/-/-/-/-/-/-/-/-/-/-/-/-/-/-/-/-/
  326.        am=amw
  327. c-------------------------------------------------------------------
  328. c  Redefinition of the numerical mach numbers
  329. c-------------------------------------------------------------------
  330.        if(abs(canc-0.0d0).lt.0.000001d0) then
  331.          top=2.0d0
  332.          bot=0.0d0
  333.        else
  334.          top=1.0d0+mhalfr*mhalfr
  335.          bot=1.0d0-mhalfr*mhalfr
  336.        endif
  337.        ml=0.5d0*(top*mlw+bot*mrw)
  338.        mr=0.5d0*(top*mrw+bot*mlw)
  339. c-------------------------------------------------------------------
  340. c  Mplus and Mmin are calligraphic lettes M+ and M- from the paper,
  341. c  see (19a) and (19b), p.367.
  342. c-------------------------------------------------------------------
  343.        if(abs(ml) .ge. 1.0d0) then
  344.          Mplus=(ml+abs(ml))/2.0d0
  345.        else
  346.          Mplus=(ml+1.0d0)*(ml+1.0d0)/4.0d0
  347.          Mplus=Mplus+beta*(ml*ml-1.0d0)*(ml*ml-1.0d0)
  348.        endif
  349.        Mplus1=(ml+abs(ml))/2.0d0
  350. c-------------------------------------------------------------------
  351.        if(abs(mr) .ge. 1.0d0) then
  352.          Mmin=(mr-abs(mr))/2.0d0
  353.        else
  354.          Mmin=-(mr-1.0d0)*(mr-1.0d0)/4.0d0
  355.          Mmin=Mmin-beta*(mr*mr-1.0d0)*(mr*mr-1.0d0)
  356.        endif
  357.        Mmin1=(mr-abs(mr))/2.0d0
  358. c-------------------------------------------------------------------
  359. c Derivatives of ml and mr with respect to both sets of primitive
  360. c variables: left  and right.
  361. c-------------------------------------------------------------------
  362.        temp_l=-un_l/(am*am)
  363.        temp_r=-un_r/(am*am)
  364. c--------
  365.        dmlr_l=temp_l*damr_l
  366.        dmlr_r=temp_l*damr_r
  367.        dmrr_l=temp_r*damr_l
  368.        dmrr_r=temp_r*damr_r
  369. c--------
  370.        dmlu_l=n_x/am+temp_l*damu_l
  371.        dmlu_r=temp_l*damu_r
  372.        dmru_l=temp_r*damu_l
  373.        dmru_r=n_x/am+temp_r*damu_r
  374. c--------
  375.        dmlv_l=n_y/am+temp_l*damv_l
  376.        dmlv_r=temp_l*damv_r
  377.        dmrv_l=temp_r*damv_l
  378.        dmrv_r=n_y/am+temp_r*damv_r
  379. c--------
  380.        dmlp_l=temp_l*damp_l
  381.        dmlp_r=temp_l*damp_r
  382.        dmrp_l=temp_r*damp_l
  383.        dmrp_r=temp_r*damp_r
  384. c-----------------------------
  385.        sr_l=dmlr_l
  386.        sr_r=dmlr_r
  387.        su_l=dmlu_l
  388.        su_r=dmlu_r
  389.        sv_l=dmlv_l
  390.        sv_r=dmlv_r
  391.        sp_l=dmlp_l
  392.        sp_r=dmlp_r
  393. c-----------------------------------------------------------------
  394. c   Redefinition of the derivatives of the ml & mr
  395. c-----------------------------------------------------------------
  396.        temp_l=(mlw-mrw)*mhalfr
  397.        temp_r=-temp_l
  398. c--------
  399.        dmlr_l=0.5d0*(top*dmlr_l+bot*dmrr_l)+canc*coef*temp_l*dmfr_l
  400.        dmlu_l=0.5d0*(top*dmlu_l+bot*dmru_l)+canc*coef*temp_l*dmfu_l
  401.        dmlv_l=0.5d0*(top*dmlv_l+bot*dmrv_l)+canc*coef*temp_l*dmfv_l
  402.        dmlp_l=0.5d0*(top*dmlp_l+bot*dmrp_l)+canc*coef*temp_l*dmfp_l
  403. c--------
  404.        dmlr_r=0.5d0*(top*dmlr_r+bot*dmrr_r)+canc*coef*temp_l*dmfr_r
  405.        dmlu_r=0.5d0*(top*dmlu_r+bot*dmru_r)+canc*coef*temp_l*dmfu_r
  406.        dmlv_r=0.5d0*(top*dmlv_r+bot*dmrv_r)+canc*coef*temp_l*dmfv_r
  407.        dmlp_r=0.5d0*(top*dmlp_r+bot*dmrp_r)+canc*coef*temp_l*dmfp_r
  408. c--------
  409.        dmrr_l=0.5d0*(top*dmrr_l+bot*sr_l)+canc*coef*temp_r*dmfr_l
  410.        dmru_l=0.5d0*(top*dmru_l+bot*su_l)+canc*coef*temp_r*dmfu_l
  411.        dmrv_l=0.5d0*(top*dmrv_l+bot*sv_l)+canc*coef*temp_r*dmfv_l
  412.        dmrp_l=0.5d0*(top*dmrp_l+bot*sp_l)+canc*coef*temp_r*dmfp_l
  413. c--------
  414.        dmrr_r=0.5d0*(top*dmrr_r+bot*sr_r)+canc*coef*temp_r*dmfr_r
  415.        dmru_r=0.5d0*(top*dmru_r+bot*su_r)+canc*coef*temp_r*dmfu_r
  416.        dmrv_r=0.5d0*(top*dmrv_r+bot*sv_r)+canc*coef*temp_r*dmfv_r
  417.        dmrp_r=0.5d0*(top*dmrp_r+bot*sp_r)+canc*coef*temp_r*dmfp_r
  418. c-----------------------------------------------------------
  419. c mmid is m_{1/2} (notation as in the paper, see (13),p.366)
  420. c-----------------------------------------------------------
  421.        mmid=Mplus+Mmin
  422. c-----------------------------------------------------------
  423. c Computing the derivatives of M+ and M-
  424. c-----------------------------------------------------------
  425.        if(ml .ge. 1.0d0) then
  426.          dMpr_l=dmlr_l
  427.          dMpu_l=dmlu_l
  428.          dMpv_l=dmlv_l
  429.          dMpp_l=dmlp_l
  430. c--------------------
  431.          dMpr_r=dmlr_r
  432.          dMpu_r=dmlu_r
  433.          dMpv_r=dmlv_r
  434.          dMpp_r=dmlp_r
  435.        else
  436.         if((ml .gt. -1.0d0) .and. (ml .lt. 1.0d0)) then
  437.          temph=(ml+1.0d0)/2.0d0
  438.          dMpr_l=(temph+4.0d0*beta*ml*(ml*ml-1.0d0))*dmlr_l
  439.          dMpu_l=(temph+4.0d0*beta*ml*(ml*ml-1.0d0))*dmlu_l
  440.          dMpv_l=(temph+4.0d0*beta*ml*(ml*ml-1.0d0))*dmlv_l
  441.          dMpp_l=(temph+4.0d0*beta*ml*(ml*ml-1.0d0))*dmlp_l
  442. c--------------------
  443.            dMpr_r=(temph+4.0d0*beta*ml*(ml*ml-1.0d0))*dmlr_r
  444.            dMpu_r=(temph+4.0d0*beta*ml*(ml*ml-1.0d0))*dmlu_r
  445.            dMpv_r=(temph+4.0d0*beta*ml*(ml*ml-1.0d0))*dmlv_r
  446.            dMpp_r=(temph+4.0d0*beta*ml*(ml*ml-1.0d0))*dmlp_r
  447.          else
  448.            dMpr_l=0.0d0
  449.            dMpu_l=0.0d0
  450.            dMpv_l=0.0d0
  451.            dMpp_l=0.0d0
  452. c---------------------
  453.            dMpr_r=0.0d0
  454.            dMpu_r=0.0d0
  455.            dMpv_r=0.0d0
  456.            dMpp_r=0.0d0
  457.          endif
  458.        endif
  459. c-----------------------------------------------------------
  460. c   addition of low Mach number
  461. c-----------------------------------------------------------
  462.        if(ml .ge. 0.0d0) then
  463.          m1pr_l=dmlr_l
  464.          m1pu_l=dmlu_l
  465.          m1pv_l=dmlv_l
  466.          m1pp_l=dmlp_l
  467. c---------------------
  468.          m1pr_r=dmlr_r
  469.          m1pu_r=dmlu_r
  470.          m1pv_r=dmlv_r
  471.          m1pp_r=dmlp_r
  472.        else
  473.          m1pr_l=0.0d0
  474.          m1pu_l=0.0d0
  475.          m1pv_l=0.0d0
  476.          m1pp_l=0.0d0
  477. c--------------------
  478.          m1pr_r=0.0d0
  479.          m1pu_r=0.0d0
  480.          m1pv_r=0.0d0
  481.          m1pp_r=0.0d0
  482.        endif
  483. c-----------------------------------------------------------
  484.        if(mr .ge. 1.0d0) then
  485.          dMmr_l=0.0d0
  486.          dMmu_l=0.0d0
  487.          dMmv_l=0.0d0
  488.          dMmp_l=0.0d0
  489. c---------------------
  490.          dMmr_r=0.0d0
  491.          dMmu_r=0.0d0
  492.          dMmv_r=0.0d0
  493.          dMmp_r=0.0d0
  494.        else
  495.          if((mr .gt. -1.0d0) .and. (mr .lt. 1.0d0)) then
  496.            temph=(-mr+1.0d0)/2.0d0
  497.            dMmr_l=(temph-4.0d0*beta*mr*(mr*mr-1.0d0))*dmrr_l
  498.            dMmu_l=(temph-4.0d0*beta*mr*(mr*mr-1.0d0))*dmru_l
  499.            dMmv_l=(temph-4.0d0*beta*mr*(mr*mr-1.0d0))*dmrv_l
  500.            dMmp_l=(temph-4.0d0*beta*mr*(mr*mr-1.0d0))*dmrp_l
  501. c--------------------
  502.            dMmr_r=(temph-4.0d0*beta*mr*(mr*mr-1.0d0))*dmrr_r
  503.            dMmu_r=(temph-4.0d0*beta*mr*(mr*mr-1.0d0))*dmru_r
  504.            dMmv_r=(temph-4.0d0*beta*mr*(mr*mr-1.0d0))*dmrv_r
  505.            dMmp_r=(temph-4.0d0*beta*mr*(mr*mr-1.0d0))*dmrp_r
  506.          else
  507.            dMmr_l=dmrr_l
  508.            dMmu_l=dmru_l
  509.            dMmv_l=dmrv_l
  510.            dMmp_l=dmrp_l
  511. c--------------------
  512.            dMmr_r=dmrr_r
  513.            dMmu_r=dmru_r
  514.            dMmv_r=dmrv_r
  515.            dMmp_r=dmrp_r
  516.          endif
  517.        endif
  518. c-----------------------------------------------------------
  519. c   addition of low Mach number
  520. c-----------------------------------------------------------
  521.        if(mr .le. 0.0d0) then
  522.          m1mr_l=dmrr_l
  523.          m1mu_l=dmru_l
  524.          m1mv_l=dmrv_l
  525.          m1mp_l=dmrp_l
  526. c---------------------
  527.          m1mr_r=dmrr_r
  528.          m1mu_r=dmru_r
  529.          m1mv_r=dmrv_r
  530.          m1mp_r=dmrp_r
  531.        else
  532.          m1mr_l=0.0d0
  533.          m1mu_l=0.0d0
  534.          m1mv_l=0.0d0
  535.          m1mp_l=0.0d0
  536. c--------------------
  537.          m1mr_r=0.0d0
  538.          m1mu_r=0.0d0
  539.          m1mv_r=0.0d0
  540.          m1mp_r=0.0d0
  541.        endif
  542. c-----------------------------------------------------------------
  543. c computing the derivatives of m_{1/2} (notation as in the paper)
  544. c-----------------------------------------------------------------
  545.        dmir_l=dMpr_l+dMmr_l
  546.        dmir_r=dMpr_r+dMmr_r
  547. c-------------
  548.        dmiu_l=dMpu_l+dMmu_l
  549.        dmiu_r=dMpu_r+dMmu_r
  550. c-------------
  551.        dmiv_l=dMpv_l+dMmv_l
  552.        dmiv_r=dMpv_r+dMmv_r
  553. c-------------
  554.        dmip_l=dMpp_l+dMmp_l
  555.        dmip_r=dMpp_r+dMmp_r
  556. c----------------------------------------------------------------
  557. c computing the main convective variables and their derivatives
  558. c mpl_m is m^{+}_{1/2} (paper's notation) and
  559. c mmin_m  is m^{-}_{1/2} (paper's notation), see (A2) on p.370.
  560. c----------------------------------------------------------------
  561.        termp=(Mmin1-Mmin+Mplus-Mplus1)*(1.0d0/(mhalfr*mhalfr)-1.0d0)
  562.        termp=termp*(pold_l-pold_r)/(pold_l/rold_l+pold_r/rold_r)
  563. c-------------------------------------------------------------
  564. c    derivatives of the termp
  565. c-------------------------------------------------------------
  566.        top=(Mmin1-Mmin+Mplus-Mplus1)
  567.        bots=1.0d0/(pold_l/rold_l+pold_r/rold_r)
  568.        bot=(pold_l-pold_r)*bots
  569.        temph=1.0d0/(mhalfr*mhalfr)-1.0d0
  570. c---------------------------
  571.        tmpr_l=(m1mr_l-dMmr_l+dMpr_l-m1pr_l)*bot*temph
  572.        tmpr_l=tmpr_l-2.0d0*bot*top*dmfr_l/(mhalfr*mhalfr*mhalfr)
  573.        tmpr_l=tmpr_l+temph*top*bot*bots*(pold_l/rold_l/rold_l)
  574. c---------------------------
  575.        tmpu_l=(m1mu_l-dMmu_l+dMpu_l-m1pu_l)*bot*temph
  576.        tmpu_l=tmpu_l-2.0d0*bot*top*dmfu_l/(mhalfr*mhalfr*mhalfr)
  577. c---------------------------
  578.        tmpv_l=(m1mv_l-dMmv_l+dMpv_l-m1pv_l)*bot*temph
  579.        tmpv_l=tmpv_l-2.0d0*bot*top*dmfv_l/(mhalfr*mhalfr*mhalfr)
  580. c---------------------------
  581.        tmpp_l=(m1mp_l-dMmp_l+dMpp_l-m1pp_l)*bot*temph
  582.        tmpp_l=tmpp_l-2.0d0*bot*top*dmfp_l/(mhalfr*mhalfr*mhalfr)
  583.        tmpp_l=tmpp_l+temph*top*bots*(1.0d0-bot/rold_l)
  584. c------------rrrrrrrr-------
  585. c------------rrrrrrrr-------
  586.        tmpr_r=(m1mr_r-dMmr_r+dMpr_r-m1pr_r)*bot*temph
  587.        tmpr_r=tmpr_r-2.0d0*bot*top*dmfr_r/(mhalfr*mhalfr*mhalfr)
  588.        tmpr_r=tmpr_r+temph*top*bot*bots*(pold_r/rold_r/rold_r)
  589. c---------------------------
  590.        tmpu_r=(m1mu_r-dMmu_r+dMpu_r-m1pu_r)*bot*temph
  591.        tmpu_r=tmpu_r-2.0d0*bot*top*dmfu_r/(mhalfr*mhalfr*mhalfr)
  592. c---------------------------
  593.        tmpv_r=(m1mv_r-dMmv_r+dMpv_r-m1pv_r)*bot*temph
  594.        tmpv_r=tmpv_r-2.0d0*bot*top*dmfv_r/(mhalfr*mhalfr*mhalfr)
  595. c---------------------------
  596.        tmpp_r=(m1mp_r-dMmp_r+dMpp_r-m1pp_r)*bot*temph
  597.        tmpp_r=tmpp_r-2.0d0*bot*top*dmfp_r/(mhalfr*mhalfr*mhalfr)
  598.        tmpp_r=tmpp_r-temph*top*bots*(1.0d0+bot/rold_r)
  599. c-------------------------------------------------------------
  600.        if(mmid .ge. 0.0d0) then
  601.          mpl_m = mmid
  602.          d2mr_l=dmir_l
  603.          d2mu_l=dmiu_l
  604.          d2mv_l=dmiv_l
  605.          d2mp_l=dmip_l
  606. c------------
  607.          d2mr_r=dmir_r
  608.          d2mu_r=dmiu_r
  609.          d2mv_r=dmiv_r
  610.          d2mp_r=dmip_r
  611. c------------
  612.        else
  613.          mpl_m = 0.0d0
  614.          d2mr_l=0.0d0
  615.          d2mu_l=0.0d0
  616.          d2mv_l=0.0d0
  617.          d2mp_l=0.0d0
  618. c------------
  619.          d2mr_r=0.0d0
  620.          d2mu_r=0.0d0
  621.          d2mv_r=0.0d0
  622.          d2mp_r=0.0d0
  623.        endif
  624. c---------------------------------------------
  625. cc       derivatives for the term termm
  626. cc------------------------------------------------------------------
  627.        if(mmid .le. 0.0d0) then
  628.          mmin_m = mmid
  629.          d3mr_l=dmir_l
  630.          d3mu_l=dmiu_l
  631.          d3mv_l=dmiv_l
  632.          d3mp_l=dmip_l
  633. c------------
  634.          d3mr_r=dmir_r
  635.          d3mu_r=dmiu_r
  636.          d3mv_r=dmiv_r
  637.          d3mp_r=dmip_r
  638. c------------
  639.        else
  640.          mmin_m = 0.0d0
  641.          d3mr_l=0.0d0
  642.          d3mu_l=0.0d0
  643.          d3mv_l=0.0d0
  644.          d3mp_l=0.0d0
  645. c------------
  646.          d3mr_r=0.0d0
  647.          d3mu_r=0.0d0
  648.          d3mv_r=0.0d0
  649.          d3mp_r=0.0d0
  650.        endif
  651. c---------------------------------------------------------------
  652. c Computing the calligraphic P+ and P- with their derivatives,
  653. c see (21a) & (21b) on p.368.
  654. c---------------------------------------------------------------
  655.        if(ml .ge. 1.0d0) then
  656.          Pplus = 1.0d0
  657.        else
  658.          if((ml .gt. -1.0d0) .and. (ml .lt. 1.0d0)) then
  659.            Pplus=(ml+1.0d0)*(ml+1.0d0)*(2.0d0-ml)/4.0d0
  660.            Pplus=Pplus+alpha*ml*(ml*ml-1.0d0)*(ml*ml-1.0d0)
  661.          else
  662.            Pplus = 0.0d0
  663.          endif
  664.        endif
  665. c---------------------------------------------------------------
  666.        if(mr .ge. 1.0d0) then
  667.          Pmin = 0.0d0
  668.        else
  669.          if((mr .gt. -1.0d0) .and. (mr .lt. 1.0d0)) then
  670.            Pmin=(mr-1.0d0)*(mr-1.0d0)*(2.0d0+mr)/4.0d0
  671.            Pmin=Pmin-alpha*mr*(mr*mr-1.0d0)*(mr*mr-1.0d0)
  672.          else
  673.            Pmin = 1.0d0
  674.          endif
  675.        endif
  676. c---------------------------------------------------------------
  677.        brac_l=(ml+1.0d0)*(2.0d0-ml)/2.0d0-(ml+1.0d0)*(ml+1.0d0)/4.0d0
  678.        brac_l=brac_l+alpha*(ml*ml-1.0d0)*(ml*ml-1.0d0)
  679.        brac_l=brac_l+4.0d0*alpha*ml*ml*(ml*ml-1.0d0)
  680. c--------------
  681.        brac_r=(mr-1.0d0)*(2.0d0+mr)/2.0d0+(mr-1.0d0)*(mr-1.0d0)/4.0d0
  682.        brac_r=brac_r-alpha*(mr*mr-1.0d0)*(mr*mr-1.0d0)
  683.        brac_r=brac_r-4.0d0*alpha*mr*mr*(mr*mr-1.0d0)
  684. c---------------------------------------------------------------
  685.        if((ml .gt. -1.0d0) .and. (ml .lt. 1.0d0)) then
  686.          dPpr_l=brac_l*dmlr_l
  687.          dPpr_r=brac_l*dmlr_r
  688. c------------
  689.          dPpu_l=brac_l*dmlu_l
  690.          dPpu_r=brac_l*dmlu_r
  691. c------------
  692.          dPpv_l=brac_l*dmlv_l
  693.          dPpv_r=brac_l*dmlv_r
  694. c------------
  695.          dPpp_l=brac_l*dmlp_l
  696.          dPpp_r=brac_l*dmlp_r
  697. c------------
  698.        else
  699.          dPpr_l=0.0d0
  700.          dPpr_r=0.0d0
  701. c-----------
  702.          dPpu_l=0.0d0
  703.          dPpu_r=0.0d0
  704. c-----------
  705.          dPpv_l=0.0d0
  706.          dPpv_r=0.0d0
  707. c-----------
  708.          dPpp_l=0.0d0
  709.          dPpp_r=0.0d0
  710. c-----------
  711.        endif
  712. c---------------------------------------------------------------
  713.        if((mr .gt. -1.0d0) .and. (mr .lt. 1.0d0)) then
  714.          dPmr_l=brac_r*dmrr_l
  715.          dPmr_r=brac_r*dmrr_r
  716. c------------
  717.          dPmu_l=brac_r*dmru_l
  718.          dPmu_r=brac_r*dmru_r
  719. c------------
  720.          dPmv_l=brac_r*dmrv_l
  721.          dPmv_r=brac_r*dmrv_r
  722. c------------
  723.          dPmp_l=brac_r*dmrp_l
  724.          dPmp_r=brac_r*dmrp_r
  725. c------------
  726.        else
  727.          dPmr_l=0.0d0
  728.          dPmr_r=0.0d0
  729. c-----------
  730.          dPmu_l=0.0d0
  731.          dPmu_r=0.0d0
  732. c-----------
  733.          dPmv_l=0.0d0
  734.          dPmv_r=0.0d0
  735. c-----------
  736.          dPmp_l=0.0d0
  737.          dPmp_r=0.0d0
  738. c-----------
  739.        endif
  740. c-------------------------------------------------------------------
  741. c computing pmid -- p_{1/2} and its derivatives, see (20b), p.367.
  742. c-------------------------------------------------------------------
  743.        dpir_l=dPpr_l*pold_l+dPmr_l*pold_r
  744.        dpiu_l=dPpu_l*pold_l+dPmu_l*pold_r
  745.        dpiv_l=dPpv_l*pold_l+dPmv_l*pold_r
  746.        dpip_l=dPpp_l*pold_l+Pplus+dPmp_l*pold_r
  747. c----------------------------
  748.        dpir_r=dPpr_r*pold_l+dPmr_r*pold_r
  749.        dpiu_r=dPpu_r*pold_l+dPmu_r*pold_r
  750.        dpiv_r=dPpv_r*pold_l+dPmv_r*pold_r
  751.        dpip_r=dPpp_r*pold_l+Pmin+dPmp_r*pold_r
  752. c---------------------------------------------------------------------
  753. c  Computation of the mass flux (rho * u)_1/2
  754. c---------------------------------------------------------------
  755.        rum=am*(mpl_m*rold_l+mmin_m*rold_r)+canc*am*termp
  756. c-------------------------------------------------------
  757.        rumr_l=damr_l*(mpl_m*rold_l+mmin_m*rold_r)+
  758.      &        am*(d2mr_l*rold_l+mpl_m+d3mr_l*rold_r)
  759.        rumr_l=rumr_l+canc*coef*(damr_l*termp+am*tmpr_l)
  760.        rumu_l=damu_l*(mpl_m*rold_l+mmin_m*rold_r)+
  761.      &        am*(d2mu_l*rold_l+d3mu_l*rold_r)
  762.        rumu_l=rumu_l+canc*coef*(damu_l*termp+am*tmpu_l)
  763.        rumv_l=damv_l*(mpl_m*rold_l+mmin_m*rold_r)+
  764.      &        am*(d2mv_l*rold_l+d3mv_l*rold_r)
  765.        rumv_l=rumv_l+canc*coef*(damv_l*termp+am*tmpv_l)
  766.        rump_l=damp_l*(mpl_m*rold_l+mmin_m*rold_r)+
  767.      &        am*(d2mp_l*rold_l+d3mp_l*rold_r)
  768.        rump_l=rump_l+canc*coef*(damp_l*termp+am*tmpp_l)
  769. c-------------------------------------------------
  770.        rumr_r=damr_r*(mpl_m*rold_l+mmin_m*rold_r)+
  771.      &        am*(d2mr_r*rold_l+mmin_m+d3mr_r*rold_r)
  772.        rumr_r=rumr_r+canc*coef*(damr_r*termp+am*tmpr_r)
  773.        rumu_r=damu_r*(mpl_m*rold_l+mmin_m*rold_r)+
  774.      &        am*(d2mu_r*rold_l+d3mu_r*rold_r)
  775.        rumu_r=rumu_r+canc*coef*(damu_r*termp+am*tmpu_r)
  776.        rumv_r=damv_r*(mpl_m*rold_l+mmin_m*rold_r)+
  777.      &        am*(d2mv_r*rold_l+d3mv_r*rold_r)
  778.        rumv_r=rumv_r+canc*coef*(damv_r*termp+am*tmpv_r)
  779.        rump_r=damp_r*(mpl_m*rold_l+mmin_m*rold_r)+
  780.      &        am*(d2mp_r*rold_l+d3mp_r*rold_r)
  781.        rump_r=rump_r+canc*coef*(damp_r*termp+am*tmpp_r)
  782. c---------------------------------------------------------------------
  783. c computing JACOBIAN as a derivative of the numerical flux function
  784. c with respect to the primitive variables.
  785. c Notation: jl(i,j) --- is the derivative of the i-component of the
  786. c           flux function with respect to the j-component of the
  787. c           vector of primitive variables of the left state.
  788. c           jr(i,j) --- is the derivative of the i-component of the
  789. c           flux function with respect to the j-component of the
  790. c           vector of primitive variables of the right state.
  791. c---------------------------------------------------------------------
  792.        jl(1,1)=rumr_l
  793.        jl(1,2)=rumu_l
  794.        jl(1,3)=rumv_l
  795.        jl(1,4)=rump_l
  796. c------------------------------------
  797.        jr(1,1)=rumr_r
  798.        jr(1,2)=rumu_r
  799.        jr(1,3)=rumv_r
  800.        jr(1,4)=rump_r
  801. c------------------------------------
  802. c---------------------------------------------------------
  803.        if(rum .ge. 0.0d0) then
  804.          br3=rumr_l*un_l
  805.          br4=rumr_l*ut_l
  806.        else
  807.          br3=rumr_l*un_r
  808.          br4=rumr_l*ut_r
  809.        endif
  810.        jl(2,1)=n_x*(br3+dpir_l)+br4*t_x
  811.        jl(3,1)=n_y*(br3+dpir_l)+br4*t_y
  812. c-------------------
  813.        if(rum .ge. 0.0d0) then
  814.          br3=rumu_l*un_l+rum*n_x
  815.          br4=rumu_l*ut_l+rum*t_x
  816.        else
  817.          br3=rumu_l*un_r
  818.          br4=rumu_l*ut_r
  819.        endif
  820.        jl(2,2)=n_x*(br3+dpiu_l)+br4*t_x
  821.        jl(3,2)=n_y*(br3+dpiu_l)+br4*t_y
  822. c-------------------
  823.        if(rum .ge. 0.0d0) then
  824.          br3=rumv_l*un_l+rum*n_y
  825.          br4=rumv_l*ut_l+rum*t_y
  826.        else
  827.          br3=rumv_l*un_r
  828.          br4=rumv_l*ut_r
  829.        endif
  830.        jl(2,3)=n_x*(br3+dpiv_l)+br4*t_x
  831.        jl(3,3)=n_y*(br3+dpiv_l)+br4*t_y
  832. c-------------------
  833.        if(rum .ge. 0.0d0) then
  834.          br3=rump_l*un_l
  835.          br4=rump_l*ut_l
  836.        else
  837.          br3=rump_l*un_r
  838.          br4=rump_l*ut_r
  839.        endif
  840.        jl(2,4)=n_x*(br3+dpip_l)+br4*t_x
  841.        jl(3,4)=n_y*(br3+dpip_l)+br4*t_y
  842. c-------------------------------------------------------------
  843. c-------------------------------------------------------------
  844.        if(rum .ge. 0.0d0) then
  845.          br3=rumr_r*un_l
  846.          br4=rumr_r*ut_l
  847.        else
  848.          br3=rumr_r*un_r
  849.          br4=rumr_r*ut_r
  850.        endif
  851.        jr(2,1)=n_x*(br3+dpir_r)+br4*t_x
  852.        jr(3,1)=n_y*(br3+dpir_r)+br4*t_y
  853. c-------------------
  854.        if(rum .ge. 0.0d0) then
  855.          br3=rumu_r*un_l
  856.          br4=rumu_r*ut_l
  857.        else
  858.          br3=rumu_r*un_r+rum*n_x
  859.          br4=rumu_r*ut_r+rum*t_x
  860.        endif
  861.        jr(2,2)=n_x*(br3+dpiu_r)+br4*t_x
  862.        jr(3,2)=n_y*(br3+dpiu_r)+br4*t_y
  863. c-------------------
  864.        if(rum .ge. 0.0d0) then
  865.          br3=rumv_r*un_l
  866.          br4=rumv_r*ut_l
  867.        else
  868.          br3=rumv_r*un_r+rum*n_y
  869.          br4=rumv_r*ut_r+rum*t_y
  870.        endif
  871.        jr(2,3)=n_x*(br3+dpiv_r)+br4*t_x
  872.        jr(3,3)=n_y*(br3+dpiv_r)+br4*t_y
  873. c-------------------
  874.        if(rum .ge. 0.0d0) then
  875.          br3=rump_r*un_l
  876.          br4=rump_r*ut_l
  877.        else
  878.          br3=rump_r*un_r
  879.          br4=rump_r*ut_r
  880.        endif
  881.        jr(2,4)=n_x*(br3+dpip_r)+br4*t_x
  882.        jr(3,4)=n_y*(br3+dpip_r)+br4*t_y
  883. c-------------------------------------------------------------
  884. c  ------  f44444444444444444444444444444 ---------
  885. c-------------------------------------------------------------
  886.        hr_l=(uold_l*uold_l+vold_l*vold_l)/2.0d0+ga*pold_l/gm1/rold_l
  887.        hr_r=(uold_r*uold_r+vold_r*vold_r)/2.0d0+ga*pold_r/gm1/rold_r
  888. c-------------------------------------------------------------
  889.        if(rum .ge. 0.0d0) then
  890.          br1=rumr_l*hr_l-rum*ga*pold_l/gm1/rold_l/rold_l
  891.        else
  892.          br1=rumr_l*hr_r
  893.        endif
  894.        jl(4,1)=br1
  895. c---------------------
  896.        if(rum .ge. 0.0d0) then
  897.          br1=rumu_l*hr_l+rum*uold_l
  898.        else
  899.          br1=rumu_l*hr_r
  900.        endif
  901.        jl(4,2)=br1
  902. c---------------------
  903.        if(rum .ge. 0.0d0) then
  904.          br1=rumv_l*hr_l+rum*vold_l
  905.        else
  906.          br1=rumv_l*hr_r
  907.        endif
  908.        jl(4,3)=br1
  909. c---------------------
  910.        if(rum .ge. 0.0d0) then
  911.          br1=rump_l*hr_l+rum*ga/gm1/rold_l
  912.        else
  913.          br1=rump_l*hr_r
  914.        endif
  915.        jl(4,4)=br1
  916. c----------------------------------------------------------
  917.        if(rum .ge. 0.0d0) then
  918.          br1=rumr_r*hr_l
  919.        else
  920.          br1=rumr_r*hr_r-rum*ga*pold_r/gm1/rold_r/rold_r
  921.        endif
  922.        jr(4,1)=br1
  923. c---------------------
  924.        if(rum .ge. 0.0d0) then
  925.          br1=rumu_r*hr_l
  926.        else
  927.          br1=rumu_r*hr_r+rum*uold_r
  928.        endif
  929.        jr(4,2)=br1
  930. c---------------------
  931.        if(rum .ge. 0.0d0) then
  932.          br1=rumv_r*hr_l
  933.        else
  934.          br1=rumv_r*hr_r+rum*vold_r
  935.        endif
  936.        jr(4,3)=br1
  937. c---------------------
  938.        if(rum .ge. 0.0d0) then
  939.          br1=rump_r*hr_l
  940.        else
  941.          br1=rump_r*hr_r+rum*ga/gm1/rold_r
  942.        endif
  943.        jr(4,4)=br1
  944. c-------------------------------------------------------------
  945. c  matrix wl(i,j) represents the derivative of the i-component
  946. c  of the vector of primitive variables of the left state with
  947. c  respect to the j-component of the vector of the conservative
  948. c  variables of the left state.
  949. c
  950. c  Here: (rho, u, v, p) - vector of primitive variables;
  951. c        (rho, rho u, rho v, rho e) -- conservative variables.
  952. c-------------------------------------------------------------
  953.        wl(1,1)=1.0d0
  954.        wl(1,2)=0.0d0
  955.        wl(1,3)=0.0d0
  956.        wl(1,4)=0.0d0
  957. c------------------------------
  958.        wl(2,1)=-uold_l/rold_l
  959.        wl(2,2)=1.0d0/rold_l
  960.        wl(2,3)=0.0d0
  961.        wl(2,4)=0.0d0
  962. c------------------------------
  963.        wl(3,1)=-vold_l/rold_l
  964.        wl(3,2)=0.0d0
  965.        wl(3,3)=1.0d0/rold_l
  966.        wl(3,4)=0.0d0
  967. c------------------------------
  968.        wl(4,1)=gm1*(uold_l*uold_l+vold_l*vold_l)/2.0d0
  969.        wl(4,2)=-uold_l*gm1
  970.        wl(4,3)=-vold_l*gm1
  971.        wl(4,4)=gm1
  972. c------------------------------
  973. c------------------------------
  974.        wr(1,1)=1.0d0
  975.        wr(1,2)=0.0d0
  976.        wr(1,3)=0.0d0
  977.        wr(1,4)=0.0d0
  978. c------------------------------
  979.        wr(2,1)=-uold_r/rold_r
  980.        wr(2,2)=1.0d0/rold_r
  981.        wr(2,3)=0.0d0
  982.        wr(2,4)=0.0d0
  983. c------------------------------
  984.        wr(3,1)=-vold_r/rold_r
  985.        wr(3,2)=0.0d0
  986.        wr(3,3)=1.0d0/rold_r
  987.        wr(3,4)=0.0d0
  988. c------------------------------
  989.        wr(4,1)=gm1*(uold_r*uold_r+vold_r*vold_r)/2.0d0
  990.        wr(4,2)=-uold_r*gm1
  991.        wr(4,3)=-vold_r*gm1
  992.        wr(4,4)=gm1
  993. c----------------------------------------------
  994.        do 1 i=1,4
  995.         do 2 j=1,4
  996.          jtl(i,j)=0.0d0
  997.          jtr(i,j)=0.0d0
  998.          do 3 k=1,4
  999.           jtl(i,j)=jtl(i,j)+jl(i,k)*wl(k,j)
  1000.           jtr(i,j)=jtr(i,j)+jr(i,k)*wr(k,j)
  1001.  3       continue
  1002.  2      continue
  1003.  1     continue
  1004. c----------------------------------------------------------------------
  1005.        return
  1006.        end
  1007.  
  1008.  
  1009.  
  1010.  
  1011.  
  1012.  
  1013.  
  1014.  
  1015.  
  1016.  
  1017.  
  1018.  
  1019.  
  1020.  
  1021.  
  1022.  
  1023.  
  1024.  
  1025.  
  1026.  
  1027.  
  1028.  
  1029.  
  1030.  
  1031.  
  1032.  
  1033.  
  1034.  
  1035.  
  1036.  
  1037.  
  1038.  
  1039.  
  1040.  

© Cast3M 2003 - Tous droits réservés.
Mentions légales