Télécharger la plaquette Cast3M


fr en

Télécharger ja3dbm.eso

Retour à la liste

Numérotation des lignes : 

ja3dbm

  1. C JA3DBM    SOURCE    CHAT      05/01/13    00:47:49     5004
  2.       SUBROUTINE JA3DBM(JTL,JTR,WVEC_L,WVEC_R,NVECT,TVEC1,TVEC2,
  3.      &                   ga,v_inf)
  4. C************************************************************************
  5. C
  6. C PROJET            :  CASTEM 2000
  7. C
  8. C NOM               :  JA3DBM (Bas-Mach jacobian for 3D)
  9. C
  10. C DESCRIPTION       :  Voir KONJA2
  11. C
  12. C LANGAGE           :  FORTRAN 77 + ESOPE 2000 (avec estensions CISI)
  13. C
  14. C AUTEUR            :  S. KUDRIAKOV, A. BECCANTINI DM2S/SFME/LTMF
  15. C
  16. C************************************************************************
  17. C
  18. c----------------------------------------------------------------------
  19. c GENERAL DESCRIPTION:
  20. c   This subroutine provides the jacobians, which are
  21. c   the derivatives of the numerical flux function,
  22. c   with respect to the right and left conservative variables
  23. c   The low-mach number corrections are made.
  24. c
  25. c EQUATIONS: 3D Euler equations of gas dynamics
  26. c
  27. c
  28. c REFERENCE: 1) JCP, 129, 364-382 (1996)
  29. c             " A Sequel to AUSM: AUSM+ ";
  30. c                M.S.Liou
  31. c            2) AIAA Journal, Sept. 1998
  32. c          "Low-Diffusion Flux-Splitting Methods for Flows at All Speeds"
  33. c             J.R.Edwards and M.S.Liou
  34. c----------------------------------------------------------------------
  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 == (gamg,rhog,pg,ung,utg,uvg)  --
  45. c                                vector of the primitive variables
  46. c                                     at the left cell;
  47. c
  48. c      wvec_r == (gamd,rhod,pd,und,utd,uvd)  --
  49. c                                vector of the primitive variables
  50. c                                     at the right cell;
  51. c
  52. c      nvect, tvec1, tvec2  --- normal and tangent vectors
  53. c                               to the cell interface
  54. c
  55. c      ga    ----  ratio of the specific heats of the gas
  56. c
  57. c      v_inf  -- parameter for choosing the reference velocity
  58. c                 when the magnitude of the physical velocity
  59. c                 is close to zero
  60.  
  61. c----------------------------------------------------------------------
  62. c
  63. c OUTPUT:
  64. c
  65. c      jtl     -- jakobian matrix 5 by 5 - derivatives of the numerical
  66. c                 flux function with respect to the conservative variables
  67. c                 from the left cell;
  68. c
  69. c      jtr     -- jakobian matrix 5 by 5 - derivatives of the numerical
  70. c                 flux function with respect to the conservative variables
  71. c                 from the right cell.
  72. c----------------------------------------------------------------------
  73.       IMPLICIT INTEGER(I-N)
  74.        real*8 wvec_l(5),wvec_r(5)
  75.        real*8 nvect(3),tvec1(3),tvec2(3)
  76.        real*8 jl(5,5),jr(5,5)
  77.        real*8 wl(5,5),wr(5,5)
  78.        real*8 jtl(5,5),jtr(5,5)
  79.        real*8 alpha,beta,upr_l,upr_r,v_inf
  80.        real*8 ga, gm1,temph
  81.        real*8 n_x,n_y,n_z,t1_x,t1_y,t1_z,t2_x,t2_y,t2_z
  82.        real*8 un_l,un_r,ut1_l,ut1_r,ut2_l,ut2_r
  83.        real*8 ml,mr,Mplus,Mmin,mmid
  84.        real*8 mpl_m, mmin_m,am
  85.        real*8 rold_l,uold_l,vold_l,wold_l,pold_l,eold_l
  86.        real*8 rold_r,uold_r,vold_r,wold_r,pold_r,eold_r
  87.        real*8 Pplus,Pmin
  88.        real*8 hr_l,hr_r
  89.        real*8 br1,br2,br3,temp_l,temp_r,brac_l,brac_r
  90.        real*8 aleft, arigh
  91.        real*8 damr_l,damr_r,damu_l,damu_r
  92.        real*8 damv_l,damv_r,damp_l,damp_r
  93.        real*8 damw_l,damw_r
  94.        real*8 dmlr_l,dmlr_r,dmlu_l,dmlu_r
  95.        real*8 dmlv_l,dmlv_r,dmlp_l,dmlp_r
  96.        real*8 dmlw_l,dmlw_r
  97.        real*8 dmrr_l,dmrr_r,dmru_l,dmru_r
  98.        real*8 dmrv_l,dmrv_r,dmrp_l,dmrp_r
  99.        real*8 dmrw_l,dmrw_r
  100.        real*8 dMpr_l,dMpr_r,dMpu_l,dMpu_r
  101.        real*8 dMpv_l,dMpv_r,dMpp_l,dMpp_r
  102.        real*8 dMpw_l,dMpw_r
  103.        real*8 dMmr_l,dMmr_r,dMmu_l,dMmu_r
  104.        real*8 dMmv_l,dMmv_r,dMmp_l,dMmp_r
  105.        real*8 dMmw_l,dMmw_r
  106.        real*8 dmir_l,dmir_r,dmiu_l,dmiu_r
  107.        real*8 dmiv_l,dmiv_r,dmip_l,dmip_r
  108.        real*8 dmiw_l,dmiw_r
  109.        real*8 d3mr_l,d3mr_r,d3mu_l,d3mu_r
  110.        real*8 d3mv_l,d3mv_r,d3mp_l,d3mp_r
  111.        real*8 d3mw_l,d3mw_r
  112.        real*8 d2mr_l,d2mr_r,d2mu_l,d2mu_r
  113.        real*8 d2mv_l,d2mv_r,d2mp_l,d2mp_r
  114.        real*8 d2mw_l,d2mw_r
  115.        real*8 dPpr_l,dPpr_r,dPpu_l,dPpu_r
  116.        real*8 dPpv_l,dPpv_r,dPpp_l,dPpp_r
  117.        real*8 dPpw_l,dPpw_r
  118.        real*8 dPmr_l,dPmr_r,dPmu_l,dPmu_r
  119.        real*8 dPmv_l,dPmv_r,dPmp_l,dPmp_r
  120.        real*8 dPmw_l,dPmw_r
  121.        real*8 dpir_l,dpir_r,dpiu_l,dpiu_r
  122.        real*8 dpiv_l,dpiv_r,dpip_l,dpip_r
  123.        real*8 dpiw_l,dpiw_r
  124.        real*8 epsil,qq,amw,Mmin1,Mplus1
  125.        real*8 fmid,mlw,mrw,termp
  126.        real*8 ur_r,ur_l,urm,mhalf,mhalfr
  127.        real*8 durr_l,durr_r,duru_l,duru_r
  128.        real*8 durv_l,durv_r,durp_l,durp_r
  129.        real*8 durw_l,durw_r
  130.        real*8 dmhr_l,dmhr_r,dmhu_l,dmhu_r
  131.        real*8 dmhv_l,dmhv_r,dmhp_l,dmhp_r
  132.        real*8 dmhw_l,dmhw_r
  133.        real*8 dmfr_l,dmfr_r,dmfu_l,dmfu_r
  134.        real*8 dmfv_l,dmfv_r,dmfp_l,dmfp_r
  135.        real*8 dmfw_l,dmfw_r
  136.        real*8 dfm_mf,dfm_mh
  137.        real*8 dfmr_l,dfmr_r,dfmu_l,dfmu_r
  138.        real*8 dfmv_l,dfmv_r,dfmp_l,dfmp_r
  139.        real*8 dfmw_l,dfmw_r
  140.        real*8 m1mr_l,m1mr_r,m1mu_l,m1mu_r
  141.        real*8 m1mv_l,m1mv_r,m1mp_l,m1mp_r
  142.        real*8 m1mw_l,m1mw_r
  143.        real*8 m1pr_l,m1pr_r,m1pu_l,m1pu_r
  144.        real*8 m1pv_l,m1pv_r,m1pp_l,m1pp_r
  145.        real*8 m1pw_l,m1pw_r
  146.        real*8 tmpr_l,tmpr_r,tmpu_l,tmpu_r
  147.        real*8 tmpv_l,tmpv_r,tmpp_l,tmpp_r
  148.        real*8 tmpw_l,tmpw_r,rum
  149.        real*8 rumr_l,rumu_l,rumv_l,rumw_l,rump_l
  150.        real*8 rumr_r,rumu_r,rumv_r,rumw_r,rump_r
  151.        real*8 sw_l,sw_r,top,bot,bots,coef,canc
  152.        real*8 sr_l,sr_r,su_l,su_r,sv_l,sv_r,sp_l,sp_r
  153.        integer i,j,k
  154.        parameter(alpha = 0.1875D0, beta = 0.125D0)
  155.        parameter(epsil = 1.0d0)
  156.        canc=1.0d0
  157.        upr_l=0.0d0
  158.        upr_r=0.0d0
  159. c-----------------------
  160.        gm1=ga-1.0d0
  161. c-----------------------
  162.        n_x=nvect(1)
  163.        n_y=nvect(2)
  164.        n_z=nvect(3)
  165. c-------------------
  166.        t1_x=tvec1(1)
  167.        t1_y=tvec1(2)
  168.        t1_z=tvec1(3)
  169. c-------------------
  170.        t2_x=tvec2(1)
  171.        t2_y=tvec2(2)
  172.        t2_z=tvec2(3)
  173. c-------------------------------------
  174.        rold_l=wvec_l(1)
  175.        uold_l=wvec_l(2)
  176.        vold_l=wvec_l(3)
  177.        wold_l=wvec_l(4)
  178.        pold_l=wvec_l(5)
  179. c-----------------------
  180.        rold_r=wvec_r(1)
  181.        uold_r=wvec_r(2)
  182.        vold_r=wvec_r(3)
  183.        wold_r=wvec_r(4)
  184.        pold_r=wvec_r(5)
  185. c-----------------------
  186.        eold_l=(uold_l*uold_l+vold_l*vold_l+wold_l*wold_l)/2.0d0
  187.        eold_l=eold_l+pold_l/(gm1*rold_l)
  188.        eold_r=(uold_r*uold_r+vold_r*vold_r+wold_r*wold_r)/2.0d0
  189.        eold_r=eold_r+pold_r/(gm1*rold_r)
  190. c------------------------------------------------------------------
  191. c   Computing reference velocity and its derivatives
  192. c   see Eq.(2) of the Ref.(2).
  193. c------------------------------------------------------------------
  194.        aleft=sqrt(ga*pold_l/rold_l)
  195.        arigh=sqrt(ga*pold_r/rold_r)
  196.        qq=sqrt(uold_l*uold_l+vold_l*vold_l+wold_l*wold_l)
  197.        if(qq .lt. (epsil*v_inf)) then
  198.           ur_l = epsil*v_inf
  199.           durr_l=0.0d0
  200.           duru_l=0.0d0
  201.           durv_l=0.0d0
  202.           durw_l=0.0d0
  203.           durp_l=0.0d0
  204.        else
  205.           ur_l=qq
  206.           durr_l=0.0d0
  207.           duru_l=uold_l/qq
  208.           durv_l=vold_l/qq
  209.           durw_l=wold_l/qq
  210.           durp_l=0.0d0
  211.        endif
  212. c-------------------------------------------------------------------
  213.        if(ur_l .ge. aleft) then
  214.             ur_l=aleft
  215.             durr_l=-aleft/(2.0d0*rold_l)
  216.             duru_l=0.0d0
  217.             durv_l=0.0d0
  218.             durw_l=0.0d0
  219.             durp_l=ga/(2.0d0*aleft*rold_l)
  220.        endif
  221. c-------------------------------------------------------------------
  222.        if(ur_l .lt. upr_l) then
  223.           ur_l=upr_l
  224.           durr_l=0.0d0
  225.           duru_l=0.0d0
  226.           durv_l=0.0d0
  227.           durw_l=0.0d0
  228.           durp_l=0.0d0
  229.        endif
  230. c-------------------------------------------------------------------
  231.        qq=sqrt(uold_r*uold_r+vold_r*vold_r+wold_r*wold_r)
  232.        if(qq .lt. (epsil*v_inf)) then
  233.           ur_r = epsil*v_inf
  234.           durr_r=0.0d0
  235.           duru_r=0.0d0
  236.           durv_r=0.0d0
  237.           durw_r=0.0d0
  238.           durp_r=0.0d0
  239.        else
  240.           ur_r=qq
  241.           durr_r=0.0d0
  242.           duru_r=uold_r/qq
  243.           durv_r=vold_r/qq
  244.           durw_r=wold_r/qq
  245.           durp_r=0.0d0
  246.        endif
  247. c------------------------------------------------------------------
  248.        if(ur_r .ge. arigh) then
  249.             ur_r=arigh
  250.             durr_r=-arigh/(2.0d0*rold_r)
  251.             duru_r=0.0d0
  252.             durv_r=0.0d0
  253.             durw_r=0.0d0
  254.             durp_r=ga/(2.0d0*arigh*rold_r)
  255.        endif
  256. c------------------------------------------------------------------
  257.        if(ur_r .lt. upr_r) then
  258.           ur_r=upr_r
  259.           durr_r=0.0d0
  260.           duru_r=0.0d0
  261.           durv_r=0.0d0
  262.           durw_r=0.0d0
  263.           durp_r=0.0d0
  264.        endif
  265. c------------------------------------------------------------------
  266. c  Reference velocity at the interface is taken as an average
  267. c  of the reference velocities of the neighbouring cells
  268. c-------------------------------------------------------------------
  269.        urm=0.5d0*(ur_l+ur_r)
  270.        durr_l=0.5d0*durr_l
  271.        duru_l=0.5d0*duru_l
  272.        durv_l=0.5d0*durv_l
  273.        durw_l=0.5d0*durw_l
  274.        durp_l=0.5d0*durp_l
  275. c-------------------------
  276.        durr_r=0.5d0*durr_r
  277.        duru_r=0.5d0*duru_r
  278.        durv_r=0.5d0*durv_r
  279.        durw_r=0.5d0*durw_r
  280.        durp_r=0.5d0*durp_r
  281. c-------------------------------------------------------------------
  282. c Computation of the speed of sound and its derivatives;
  283. c numerical speed of sound at the interface is taken as an average
  284. c of the speeds of sounds of the neighbouring cells
  285. c-------------------------------------------------------------------
  286.        am=0.5d0*(aleft+arigh)
  287. c-------------------------------------------------------------------
  288.        if(abs(urm/am-1.0d0).le.0.000001d0) then
  289.          coef=0.0d0
  290.        else
  291.          coef=1.0d0
  292.        endif
  293. c-------------------------------------------------------------------
  294. c Computation of the speed of sound and its derivatives;
  295. c numerical speed of sound at the interface is taken as an average
  296. c of the speeds of sounds of the neighbouring cells
  297. c--------------------------------------------------------------------
  298.        damr_r=-arigh/(4.0d0*rold_r)
  299.        damu_r=0.0d0
  300.        damv_r=0.0d0
  301.        damw_r=0.0d0
  302.        damp_r=ga/(4.0d0*arigh*rold_r)
  303. c-----------------------
  304.        damr_l=-aleft/(4.0d0*rold_l)
  305.        damu_l=0.0d0
  306.        damv_l=0.0d0
  307.        damw_l=0.0d0
  308.        damp_l=ga/(4.0d0*aleft*rold_l)
  309. c----------------------------------------------------------------------
  310. c  computing numerical Mach number and its derivatives
  311. c----------------------------------------------------------------------
  312.        un_l=uold_l*n_x+vold_l*n_y+wold_l*n_z
  313.        un_r=uold_r*n_x+vold_r*n_y+wold_r*n_z
  314. c----------
  315.        ut1_l=uold_l*t1_x+vold_l*t1_y+wold_l*t1_z
  316.        ut1_r=uold_r*t1_x+vold_r*t1_y+wold_r*t1_z
  317. c----------
  318.        ut2_l=uold_l*t2_x+vold_l*t2_y+wold_l*t2_z
  319.        ut2_r=uold_r*t2_x+vold_r*t2_y+wold_r*t2_z
  320. c----------------------------------------
  321.        ml=un_l/am
  322.        mr=un_r/am
  323.        mhalf=0.5d0*(un_l+un_r)/am
  324. c---------------------
  325.        top=0.5d0*(un_l+un_r)/(am*am)
  326.        dmhr_l=-top*damr_l
  327.        dmhu_l=n_x/2.0d0/am-top*damu_l
  328.        dmhv_l=n_y/2.0d0/am-top*damv_l
  329.        dmhw_l=n_z/2.0d0/am-top*damw_l
  330.        dmhp_l=-top*damp_l
  331. c---------------------
  332.        dmhr_r=-top*damr_r
  333.        dmhu_r=n_x/2.0d0/am-top*damu_r
  334.        dmhv_r=n_y/2.0d0/am-top*damv_r
  335.        dmhw_r=n_z/2.0d0/am-top*damw_r
  336.        dmhp_r=-top*damp_r
  337. c--------------------------------
  338.        mhalfr=urm/am
  339. c--------------------------------
  340.        top=urm/(am*am)
  341.        dmfr_l=durr_l/am-top*damr_l
  342.        dmfu_l=duru_l/am-top*damu_l
  343.        dmfv_l=durv_l/am-top*damv_l
  344.        dmfw_l=durw_l/am-top*damw_l
  345.        dmfp_l=durp_l/am-top*damp_l
  346. c---------------------
  347.        dmfr_r=durr_r/am-top*damr_r
  348.        dmfu_r=duru_r/am-top*damu_r
  349.        dmfv_r=durv_r/am-top*damv_r
  350.        dmfw_r=durw_r/am-top*damw_r
  351.        dmfp_r=durp_r/am-top*damp_r
  352. c-------------------------------------------------------------------
  353. c Scaling function for the speed of sound and its derivatives
  354. c see Eq.(32) of the Ref.2)
  355. c-------------------------------------------------------------------
  356.        top=(1.0d0-mhalfr*mhalfr)*(1.0d0-mhalfr*mhalfr)
  357.        top=top*mhalf*mhalf+4.0d0*mhalfr*mhalfr
  358.        bot=1.0d0+mhalfr*mhalfr
  359.        if(abs(canc-0.0d0).lt.0.000001d0) then
  360.          fmid=1.0d0
  361.        else
  362.          fmid=sqrt(top)/bot
  363.        endif
  364. c--------------------------
  365.        temph=-4.0d0*(1.0d0-mhalfr*mhalfr)*mhalfr
  366.        temph=temph*mhalf*mhalf+8.0d0*mhalfr
  367.        dfm_mf=temph/(sqrt(top)*2.0d0*bot)
  368.        dfm_mf=dfm_mf-sqrt(top)*2.0d0*mhalfr/(bot*bot)
  369. c--------------------------
  370.        temph=2.0d0*(1.0d0-mhalfr*mhalfr)*(1.0d0-mhalfr*mhalfr)*mhalf
  371.        dfm_mh=temph/(2.0d0*bot*sqrt(top))
  372. c--------------------------
  373.        dfmr_l=dfm_mf*dmfr_l+dfm_mh*dmhr_l
  374.        dfmu_l=dfm_mf*dmfu_l+dfm_mh*dmhu_l
  375.        dfmv_l=dfm_mf*dmfv_l+dfm_mh*dmhv_l
  376.        dfmw_l=dfm_mf*dmfw_l+dfm_mh*dmhw_l
  377.        dfmp_l=dfm_mf*dmfp_l+dfm_mh*dmhp_l
  378. c--------------------------
  379.        dfmr_r=dfm_mf*dmfr_r+dfm_mh*dmhr_r
  380.        dfmu_r=dfm_mf*dmfu_r+dfm_mh*dmhu_r
  381.        dfmv_r=dfm_mf*dmfv_r+dfm_mh*dmhv_r
  382.        dfmw_r=dfm_mf*dmfw_r+dfm_mh*dmhw_r
  383.        dfmp_r=dfm_mf*dmfp_r+dfm_mh*dmhp_r
  384. c------------------------------------------------------------------
  385. c 'New' speed of sound 'amw' defined as a product of the scaling
  386. c  function 'fmid' and the 'Old' speed of sound 'am'; see (31) of Ref.2)
  387. c------------------------------------------------------------------
  388.        amw=fmid*am
  389.        mlw=un_l/amw
  390.        mrw=un_r/amw
  391. c--------------------------
  392.        damr_l=canc*coef*dfmr_l*am+fmid*damr_l
  393.        damu_l=canc*coef*dfmu_l*am+fmid*damu_l
  394.        damv_l=canc*coef*dfmv_l*am+fmid*damv_l
  395.        damw_l=canc*coef*dfmw_l*am+fmid*damw_l
  396.        damp_l=canc*coef*dfmp_l*am+fmid*damp_l
  397. c--------------------------
  398.        damr_r=canc*coef*dfmr_r*am+fmid*damr_r
  399.        damu_r=canc*coef*dfmu_r*am+fmid*damu_r
  400.        damv_r=canc*coef*dfmv_r*am+fmid*damv_r
  401.        damw_r=canc*coef*dfmw_r*am+fmid*damw_r
  402.        damp_r=canc*coef*dfmp_r*am+fmid*damp_r
  403. c-/-/-/-/-/-/-/-/-/-/-/-/-/-/-/-/-/
  404.        am=amw
  405. c-------------------------------------------------------------------
  406. c  Redefinition of the numerical mach numbers;
  407. c  see Eqs. (33) and (34) of the Ref.2)
  408. c-------------------------------------------------------------------
  409.        if(abs(canc-0.0d0).lt.0.000001d0) then
  410.          top=2.0d0
  411.          bot=0.0d0
  412.        else
  413.          top=1.0d0+mhalfr*mhalfr
  414.          bot=1.0d0-mhalfr*mhalfr
  415.        endif
  416.        ml=0.5d0*(top*mlw+bot*mrw)
  417.        mr=0.5d0*(top*mrw+bot*mlw)
  418. c-------------------------------------------------------------------
  419. c  Mplus and Mmin are calligraphic lettes M+ and M- from the paper,
  420. c  see (19a) and (19b), p.367. of the Ref.1)
  421. c-------------------------------------------------------------------
  422.        if(abs(ml) .ge. 1.0d0) then
  423.          Mplus=(ml+abs(ml))/2.0d0
  424.        else
  425.          Mplus=(ml+1.0d0)*(ml+1.0d0)/4.0d0
  426.          Mplus=Mplus+beta*(ml*ml-1.0d0)*(ml*ml-1.0d0)
  427.        endif
  428.        Mplus1=(ml+abs(ml))/2.0d0
  429. c-------------------------------------------------------------------
  430.        if(abs(mr) .ge. 1.0d0) then
  431.          Mmin=(mr-abs(mr))/2.0d0
  432.        else
  433.          Mmin=-(mr-1.0d0)*(mr-1.0d0)/4.0d0
  434.          Mmin=Mmin-beta*(mr*mr-1.0d0)*(mr*mr-1.0d0)
  435.        endif
  436.        Mmin1=(mr-abs(mr))/2.0d0
  437. c-------------------------------------------------------------------
  438. c---------------------------------------
  439. c Derivatives of ml and mr with respect
  440. c to both sets of primitive variables:
  441. c left  and right
  442. c--------------------------------------
  443.        temp_l=-un_l/(am*am)
  444.        temp_r=-un_r/(am*am)
  445. c--------
  446.        dmlr_l=temp_l*damr_l
  447.        dmlr_r=temp_l*damr_r
  448.        dmrr_l=temp_r*damr_l
  449.        dmrr_r=temp_r*damr_r
  450. c--------
  451.        dmlu_l=n_x/am+temp_l*damu_l
  452.        dmlu_r=temp_l*damu_r
  453.        dmru_l=temp_r*damu_l
  454.        dmru_r=n_x/am+temp_r*damu_r
  455. c--------
  456.        dmlv_l=n_y/am+temp_l*damv_l
  457.        dmlv_r=temp_l*damv_r
  458.        dmrv_l=temp_r*damv_l
  459.        dmrv_r=n_y/am+temp_r*damv_r
  460. c--------
  461.        dmlw_l=n_z/am+temp_l*damw_l
  462.        dmlw_r=temp_l*damw_r
  463.        dmrw_l=temp_r*damw_l
  464.        dmrw_r=n_z/am+temp_r*damw_r
  465. c--------
  466.        dmlp_l=temp_l*damp_l
  467.        dmlp_r=temp_l*damp_r
  468.        dmrp_l=temp_r*damp_l
  469.        dmrp_r=temp_r*damp_r
  470. c-----------------------------
  471.        sr_l=dmlr_l
  472.        sr_r=dmlr_r
  473.        su_l=dmlu_l
  474.        su_r=dmlu_r
  475.        sv_l=dmlv_l
  476.        sv_r=dmlv_r
  477.        sw_l=dmlw_l
  478.        sw_r=dmlw_r
  479.        sp_l=dmlp_l
  480.        sp_r=dmlp_r
  481. c-----------------------------------------------------------------
  482. c   Redefinition of the derivatives of the ml & mr
  483. c-----------------------------------------------------------------
  484.        temp_l=(mlw-mrw)*mhalfr
  485.        temp_r=-temp_l
  486. c--------
  487.        dmlr_l=0.5d0*(top*dmlr_l+bot*dmrr_l)+canc*coef*temp_l*dmfr_l
  488.        dmlu_l=0.5d0*(top*dmlu_l+bot*dmru_l)+canc*coef*temp_l*dmfu_l
  489.        dmlv_l=0.5d0*(top*dmlv_l+bot*dmrv_l)+canc*coef*temp_l*dmfv_l
  490.        dmlw_l=0.5d0*(top*dmlw_l+bot*dmrw_l)+canc*coef*temp_l*dmfw_l
  491.        dmlp_l=0.5d0*(top*dmlp_l+bot*dmrp_l)+canc*coef*temp_l*dmfp_l
  492. c--------
  493.        dmlr_r=0.5d0*(top*dmlr_r+bot*dmrr_r)+canc*coef*temp_l*dmfr_r
  494.        dmlu_r=0.5d0*(top*dmlu_r+bot*dmru_r)+canc*coef*temp_l*dmfu_r
  495.        dmlv_r=0.5d0*(top*dmlv_r+bot*dmrv_r)+canc*coef*temp_l*dmfv_r
  496.        dmlw_r=0.5d0*(top*dmlw_r+bot*dmrw_r)+canc*coef*temp_l*dmfw_r
  497.        dmlp_r=0.5d0*(top*dmlp_r+bot*dmrp_r)+canc*coef*temp_l*dmfp_r
  498. c--------
  499.        dmrr_l=0.5d0*(top*dmrr_l+bot*sr_l)+canc*coef*temp_r*dmfr_l
  500.        dmru_l=0.5d0*(top*dmru_l+bot*su_l)+canc*coef*temp_r*dmfu_l
  501.        dmrv_l=0.5d0*(top*dmrv_l+bot*sv_l)+canc*coef*temp_r*dmfv_l
  502.        dmrw_l=0.5d0*(top*dmrw_l+bot*sw_l)+canc*coef*temp_r*dmfw_l
  503.        dmrp_l=0.5d0*(top*dmrp_l+bot*sp_l)+canc*coef*temp_r*dmfp_l
  504. c--------
  505.        dmrr_r=0.5d0*(top*dmrr_r+bot*sr_r)+canc*coef*temp_r*dmfr_r
  506.        dmru_r=0.5d0*(top*dmru_r+bot*su_r)+canc*coef*temp_r*dmfu_r
  507.        dmrv_r=0.5d0*(top*dmrv_r+bot*sv_r)+canc*coef*temp_r*dmfv_r
  508.        dmrw_r=0.5d0*(top*dmrw_r+bot*sw_r)+canc*coef*temp_r*dmfw_r
  509.        dmrp_r=0.5d0*(top*dmrp_r+bot*sp_r)+canc*coef*temp_r*dmfp_r
  510. c-----------------------------------------------------------
  511. c mmid is m_{1/2} (notation as in the paper)
  512. c see Eq.(13),p.366 of the Ref.1)
  513. c----------------------------------------------------------
  514.        mmid=Mplus+Mmin
  515. c-----------------------------------------------------------
  516. c Computing the derivatives of M+ and M-
  517. c-----------------------------------------------------------
  518.        if(ml .ge. 1.0d0) then
  519.          dMpr_l=dmlr_l
  520.          dMpu_l=dmlu_l
  521.          dMpv_l=dmlv_l
  522.          dMpw_l=dmlw_l
  523.          dMpp_l=dmlp_l
  524. c--------------------
  525.          dMpr_r=dmlr_r
  526.          dMpu_r=dmlu_r
  527.          dMpv_r=dmlv_r
  528.          dMpw_r=dmlw_r
  529.          dMpp_r=dmlp_r
  530.        else
  531.         if((ml .gt. -1.0d0) .and. (ml .lt. 1.0d0)) then
  532.          temph=(ml+1.0d0)/2.0d0
  533.          dMpr_l=(temph+4.0d0*beta*ml*(ml*ml-1.0d0))*dmlr_l
  534.          dMpu_l=(temph+4.0d0*beta*ml*(ml*ml-1.0d0))*dmlu_l
  535.          dMpv_l=(temph+4.0d0*beta*ml*(ml*ml-1.0d0))*dmlv_l
  536.          dMpw_l=(temph+4.0d0*beta*ml*(ml*ml-1.0d0))*dmlw_l
  537.          dMpp_l=(temph+4.0d0*beta*ml*(ml*ml-1.0d0))*dmlp_l
  538. c--------------------
  539.            dMpr_r=(temph+4.0d0*beta*ml*(ml*ml-1.0d0))*dmlr_r
  540.            dMpu_r=(temph+4.0d0*beta*ml*(ml*ml-1.0d0))*dmlu_r
  541.            dMpv_r=(temph+4.0d0*beta*ml*(ml*ml-1.0d0))*dmlv_r
  542.            dMpw_r=(temph+4.0d0*beta*ml*(ml*ml-1.0d0))*dmlw_r
  543.            dMpp_r=(temph+4.0d0*beta*ml*(ml*ml-1.0d0))*dmlp_r
  544.          else
  545.            dMpr_l=0.0d0
  546.            dMpu_l=0.0d0
  547.            dMpv_l=0.0d0
  548.            dMpw_l=0.0d0
  549.            dMpp_l=0.0d0
  550. c---------------------
  551.            dMpr_r=0.0d0
  552.            dMpu_r=0.0d0
  553.            dMpv_r=0.0d0
  554.            dMpw_r=0.0d0
  555.            dMpp_r=0.0d0
  556.          endif
  557.        endif
  558. c-----------------------------------------------------------
  559. c   addition of low Mach number terms
  560. c-----------------------------------------------------------
  561.        if(ml .ge. 0.0d0) then
  562.          m1pr_l=dmlr_l
  563.          m1pu_l=dmlu_l
  564.          m1pv_l=dmlv_l
  565.          m1pw_l=dmlw_l
  566.          m1pp_l=dmlp_l
  567. c---------------------
  568.          m1pr_r=dmlr_r
  569.          m1pu_r=dmlu_r
  570.          m1pv_r=dmlv_r
  571.          m1pw_r=dmlw_r
  572.          m1pp_r=dmlp_r
  573.        else
  574.          m1pr_l=0.0d0
  575.          m1pu_l=0.0d0
  576.          m1pv_l=0.0d0
  577.          m1pw_l=0.0d0
  578.          m1pp_l=0.0d0
  579. c--------------------
  580.          m1pr_r=0.0d0
  581.          m1pu_r=0.0d0
  582.          m1pv_r=0.0d0
  583.          m1pw_r=0.0d0
  584.          m1pp_r=0.0d0
  585.        endif
  586. c-----------------------------------------------------------
  587. c-----------------------------------------------------------
  588.        if(mr .ge. 1.0d0) then
  589.          dMmr_l=0.0d0
  590.          dMmu_l=0.0d0
  591.          dMmv_l=0.0d0
  592.          dMmw_l=0.0d0
  593.          dMmp_l=0.0d0
  594. c---------------------
  595.          dMmr_r=0.0d0
  596.          dMmu_r=0.0d0
  597.          dMmv_r=0.0d0
  598.          dMmw_r=0.0d0
  599.          dMmp_r=0.0d0
  600.        else
  601.          if((mr .gt. -1.0d0) .and. (mr .lt. 1.0d0)) then
  602.            temph=(-mr+1.0d0)/2.0d0
  603.            dMmr_l=(temph-4.0d0*beta*mr*(mr*mr-1.0d0))*dmrr_l
  604.            dMmu_l=(temph-4.0d0*beta*mr*(mr*mr-1.0d0))*dmru_l
  605.            dMmv_l=(temph-4.0d0*beta*mr*(mr*mr-1.0d0))*dmrv_l
  606.            dMmw_l=(temph-4.0d0*beta*mr*(mr*mr-1.0d0))*dmrw_l
  607.            dMmp_l=(temph-4.0d0*beta*mr*(mr*mr-1.0d0))*dmrp_l
  608. c--------------------
  609.            dMmr_r=(temph-4.0d0*beta*mr*(mr*mr-1.0d0))*dmrr_r
  610.            dMmu_r=(temph-4.0d0*beta*mr*(mr*mr-1.0d0))*dmru_r
  611.            dMmv_r=(temph-4.0d0*beta*mr*(mr*mr-1.0d0))*dmrv_r
  612.            dMmw_r=(temph-4.0d0*beta*mr*(mr*mr-1.0d0))*dmrw_r
  613.            dMmp_r=(temph-4.0d0*beta*mr*(mr*mr-1.0d0))*dmrp_r
  614.          else
  615.            dMmr_l=dmrr_l
  616.            dMmu_l=dmru_l
  617.            dMmv_l=dmrv_l
  618.            dMmw_l=dmrw_l
  619.            dMmp_l=dmrp_l
  620. c--------------------
  621.            dMmr_r=dmrr_r
  622.            dMmu_r=dmru_r
  623.            dMmv_r=dmrv_r
  624.            dMmw_r=dmrw_r
  625.            dMmp_r=dmrp_r
  626.          endif
  627.        endif
  628. c-----------------------------------------------------------
  629. c   addition of low Mach number terms
  630. c-----------------------------------------------------------
  631.        if(mr .le. 0.0d0) then
  632.          m1mr_l=dmrr_l
  633.          m1mu_l=dmru_l
  634.          m1mv_l=dmrv_l
  635.          m1mw_l=dmrw_l
  636.          m1mp_l=dmrp_l
  637. c---------------------
  638.          m1mr_r=dmrr_r
  639.          m1mu_r=dmru_r
  640.          m1mv_r=dmrv_r
  641.          m1mw_r=dmrw_r
  642.          m1mp_r=dmrp_r
  643.        else
  644.          m1mr_l=0.0d0
  645.          m1mu_l=0.0d0
  646.          m1mv_l=0.0d0
  647.          m1mw_l=0.0d0
  648.          m1mp_l=0.0d0
  649. c--------------------
  650.          m1mr_r=0.0d0
  651.          m1mu_r=0.0d0
  652.          m1mv_r=0.0d0
  653.          m1mw_r=0.0d0
  654.          m1mp_r=0.0d0
  655.        endif
  656. c-----------------------------------------------------------------
  657. c computing the derivatives of m_{1/2} (notation as in the paper)
  658. c-----------------------------------------------------------------
  659.        dmir_l=dMpr_l+dMmr_l
  660.        dmir_r=dMpr_r+dMmr_r
  661. c-------------
  662.        dmiu_l=dMpu_l+dMmu_l
  663.        dmiu_r=dMpu_r+dMmu_r
  664. c-------------
  665.        dmiv_l=dMpv_l+dMmv_l
  666.        dmiv_r=dMpv_r+dMmv_r
  667. c-------------
  668.        dmiw_l=dMpw_l+dMmw_l
  669.        dmiw_r=dMpw_r+dMmw_r
  670. c-------------
  671.        dmip_l=dMpp_l+dMmp_l
  672.        dmip_r=dMpp_r+dMmp_r
  673. c-------------------------------------------------------------
  674. c computing the main convective variables and their derivatives
  675. c mpl_m is m^{+}_{1/2} (paper's notation) and
  676. c mmin_m  is m^{-}_{1/2} (paper's notation),
  677. c see Eq.(A2) on p.370 of the Ref.1)
  678. c----------------------------------------------------------------
  679.        termp=(Mmin1-Mmin+Mplus-Mplus1)*(1.0d0/(mhalfr*mhalfr)-1.0d0)
  680.        termp=termp*(pold_l-pold_r)/(pold_l/rold_l+pold_r/rold_r)
  681. c-------------------------------------------------------------
  682. c    derivatives of the termp
  683. c-------------------------------------------------------------
  684.        top=(Mmin1-Mmin+Mplus-Mplus1)
  685.        bots=1.0d0/(pold_l/rold_l+pold_r/rold_r)
  686.        bot=(pold_l-pold_r)*bots
  687.        temph=1.0d0/(mhalfr*mhalfr)-1.0d0
  688. c---------------------------
  689.        tmpr_l=(m1mr_l-dMmr_l+dMpr_l-m1pr_l)*bot*temph
  690.        tmpr_l=tmpr_l-2.0d0*bot*top*dmfr_l/(mhalfr*mhalfr*mhalfr)
  691.        tmpr_l=tmpr_l+temph*top*bot*bots*(pold_l/rold_l/rold_l)
  692. c---------------------------
  693.        tmpu_l=(m1mu_l-dMmu_l+dMpu_l-m1pu_l)*bot*temph
  694.        tmpu_l=tmpu_l-2.0d0*bot*top*dmfu_l/(mhalfr*mhalfr*mhalfr)
  695. c---------------------------
  696.        tmpv_l=(m1mv_l-dMmv_l+dMpv_l-m1pv_l)*bot*temph
  697.        tmpv_l=tmpv_l-2.0d0*bot*top*dmfv_l/(mhalfr*mhalfr*mhalfr)
  698. c---------------------------
  699.        tmpw_l=(m1mw_l-dMmw_l+dMpw_l-m1pw_l)*bot*temph
  700.        tmpw_l=tmpw_l-2.0d0*bot*top*dmfw_l/(mhalfr*mhalfr*mhalfr)
  701. c---------------------------
  702.        tmpp_l=(m1mp_l-dMmp_l+dMpp_l-m1pp_l)*bot*temph
  703.        tmpp_l=tmpp_l-2.0d0*bot*top*dmfp_l/(mhalfr*mhalfr*mhalfr)
  704.        tmpp_l=tmpp_l+temph*top*bots*(1.0d0-bot/rold_l)
  705. c------------rrrrrrrr-------
  706. c------------rrrrrrrr-------
  707.        tmpr_r=(m1mr_r-dMmr_r+dMpr_r-m1pr_r)*bot*temph
  708.        tmpr_r=tmpr_r-2.0d0*bot*top*dmfr_r/(mhalfr*mhalfr*mhalfr)
  709.        tmpr_r=tmpr_r+temph*top*bot*bots*(pold_r/rold_r/rold_r)
  710. c---------------------------
  711.        tmpu_r=(m1mu_r-dMmu_r+dMpu_r-m1pu_r)*bot*temph
  712.        tmpu_r=tmpu_r-2.0d0*bot*top*dmfu_r/(mhalfr*mhalfr*mhalfr)
  713. c---------------------------
  714.        tmpv_r=(m1mv_r-dMmv_r+dMpv_r-m1pv_r)*bot*temph
  715.        tmpv_r=tmpv_r-2.0d0*bot*top*dmfv_r/(mhalfr*mhalfr*mhalfr)
  716. c---------------------------
  717.        tmpw_r=(m1mw_r-dMmw_r+dMpw_r-m1pw_r)*bot*temph
  718.        tmpw_r=tmpw_r-2.0d0*bot*top*dmfw_r/(mhalfr*mhalfr*mhalfr)
  719. c---------------------------
  720.        tmpp_r=(m1mp_r-dMmp_r+dMpp_r-m1pp_r)*bot*temph
  721.        tmpp_r=tmpp_r-2.0d0*bot*top*dmfp_r/(mhalfr*mhalfr*mhalfr)
  722.        tmpp_r=tmpp_r-temph*top*bots*(1.0d0+bot/rold_r)
  723. c-------------------------------------------------------------
  724.        if(mmid .ge. 0.0d0) then
  725.          mpl_m = mmid
  726.          d2mr_l=dmir_l
  727.          d2mu_l=dmiu_l
  728.          d2mv_l=dmiv_l
  729.          d2mw_l=dmiw_l
  730.          d2mp_l=dmip_l
  731. c------------
  732.          d2mr_r=dmir_r
  733.          d2mu_r=dmiu_r
  734.          d2mv_r=dmiv_r
  735.          d2mw_r=dmiw_r
  736.          d2mp_r=dmip_r
  737. c------------
  738.        else
  739.          mpl_m = 0.0d0
  740.          d2mr_l=0.0d0
  741.          d2mu_l=0.0d0
  742.          d2mv_l=0.0d0
  743.          d2mw_l=0.0d0
  744.          d2mp_l=0.0d0
  745. c------------
  746.          d2mr_r=0.0d0
  747.          d2mu_r=0.0d0
  748.          d2mv_r=0.0d0
  749.          d2mw_r=0.0d0
  750.          d2mp_r=0.0d0
  751.        endif
  752. c------------------------------------------------------------------
  753.        if(mmid .le. 0.0d0) then
  754.          mmin_m = mmid
  755.          d3mr_l=dmir_l
  756.          d3mu_l=dmiu_l
  757.          d3mv_l=dmiv_l
  758.          d3mw_l=dmiw_l
  759.          d3mp_l=dmip_l
  760. c------------
  761.          d3mr_r=dmir_r
  762.          d3mu_r=dmiu_r
  763.          d3mv_r=dmiv_r
  764.          d3mw_r=dmiw_r
  765.          d3mp_r=dmip_r
  766. c------------
  767.        else
  768.          mmin_m = 0.0d0
  769.          d3mr_l=0.0d0
  770.          d3mu_l=0.0d0
  771.          d3mv_l=0.0d0
  772.          d3mw_l=0.0d0
  773.          d3mp_l=0.0d0
  774. c------------
  775.          d3mr_r=0.0d0
  776.          d3mu_r=0.0d0
  777.          d3mv_r=0.0d0
  778.          d3mw_r=0.0d0
  779.          d3mp_r=0.0d0
  780.        endif
  781. c---------------------------------------------------------------
  782. c Computing the calligraphic P+ and P- with their derivatives
  783. c---------------------------------------------------------------
  784.        if(ml .ge. 1.0d0) then
  785.          Pplus = 1.0d0
  786.        else
  787.          if((ml .gt. -1.0d0) .and. (ml .lt. 1.0d0)) then
  788.            Pplus=(ml+1.0d0)*(ml+1.0d0)*(2.0d0-ml)/4.0d0
  789.            Pplus=Pplus+alpha*ml*(ml*ml-1.0d0)*(ml*ml-1.0d0)
  790.          else
  791.            Pplus = 0.0d0
  792.          endif
  793.        endif
  794. c---------------------------------------------------------------
  795.        if(mr .ge. 1.0d0) then
  796.          Pmin = 0.0d0
  797.        else
  798.          if((mr .gt. -1.0d0) .and. (mr .lt. 1.0d0)) then
  799.            Pmin=(mr-1.0d0)*(mr-1.0d0)*(2.0d0+mr)/4.0d0
  800.            Pmin=Pmin-alpha*mr*(mr*mr-1.0d0)*(mr*mr-1.0d0)
  801.          else
  802.            Pmin = 1.0d0
  803.          endif
  804.        endif
  805. c---------------------------------------------------------------
  806.        brac_l=(ml+1.0d0)*(2.0d0-ml)/2.0d0-(ml+1.0d0)*(ml+1.0d0)/4.0d0
  807.        brac_l=brac_l+alpha*(ml*ml-1.0d0)*(ml*ml-1.0d0)
  808.        brac_l=brac_l+4.0d0*alpha*ml*ml*(ml*ml-1.0d0)
  809. c--------------
  810.        brac_r=(mr-1.0d0)*(2.0d0+mr)/2.0d0+(mr-1.0d0)*(mr-1.0d0)/4.0d0
  811.        brac_r=brac_r-alpha*(mr*mr-1.0d0)*(mr*mr-1.0d0)
  812.        brac_r=brac_r-4.0d0*alpha*mr*mr*(mr*mr-1.0d0)
  813. c---------------------------------------------------------------
  814.        if((ml .gt. -1.0d0) .and. (ml .lt. 1.0d0)) then
  815.          dPpr_l=brac_l*dmlr_l
  816.          dPpr_r=brac_l*dmlr_r
  817. c------------
  818.          dPpu_l=brac_l*dmlu_l
  819.          dPpu_r=brac_l*dmlu_r
  820. c------------
  821.          dPpv_l=brac_l*dmlv_l
  822.          dPpv_r=brac_l*dmlv_r
  823. c------------
  824.          dPpw_l=brac_l*dmlw_l
  825.          dPpw_r=brac_l*dmlw_r
  826. c------------
  827.          dPpp_l=brac_l*dmlp_l
  828.          dPpp_r=brac_l*dmlp_r
  829. c------------
  830.        else
  831.          dPpr_l=0.0d0
  832.          dPpr_r=0.0d0
  833. c-----------
  834.          dPpu_l=0.0d0
  835.          dPpu_r=0.0d0
  836. c-----------
  837.          dPpv_l=0.0d0
  838.          dPpv_r=0.0d0
  839. c-----------
  840.          dPpw_l=0.0d0
  841.          dPpw_r=0.0d0
  842. c-----------
  843.          dPpp_l=0.0d0
  844.          dPpp_r=0.0d0
  845. c-----------
  846.        endif
  847. c---------------------------------------------------------------
  848.        if((mr .gt. -1.0d0) .and. (mr .lt. 1.0d0)) then
  849.          dPmr_l=brac_r*dmrr_l
  850.          dPmr_r=brac_r*dmrr_r
  851. c------------
  852.          dPmu_l=brac_r*dmru_l
  853.          dPmu_r=brac_r*dmru_r
  854. c------------
  855.          dPmv_l=brac_r*dmrv_l
  856.          dPmv_r=brac_r*dmrv_r
  857. c------------
  858.          dPmw_l=brac_r*dmrw_l
  859.          dPmw_r=brac_r*dmrw_r
  860. c------------
  861.          dPmp_l=brac_r*dmrp_l
  862.          dPmp_r=brac_r*dmrp_r
  863. c------------
  864.        else
  865.          dPmr_l=0.0d0
  866.          dPmr_r=0.0d0
  867. c-----------
  868.          dPmu_l=0.0d0
  869.          dPmu_r=0.0d0
  870. c-----------
  871.          dPmv_l=0.0d0
  872.          dPmv_r=0.0d0
  873. c-----------
  874.          dPmw_l=0.0d0
  875.          dPmw_r=0.0d0
  876. c-----------
  877.          dPmp_l=0.0d0
  878.          dPmp_r=0.0d0
  879. c-----------
  880.        endif
  881. c---------------------------------------------------------------------
  882. c computing pmid -- p_{1/2} and its derivatives
  883. c---------------------------------------------------------------------
  884.        dpir_l=dPpr_l*pold_l+dPmr_l*pold_r
  885.        dpiu_l=dPpu_l*pold_l+dPmu_l*pold_r
  886.        dpiv_l=dPpv_l*pold_l+dPmv_l*pold_r
  887.        dpiw_l=dPpw_l*pold_l+dPmw_l*pold_r
  888.        dpip_l=dPpp_l*pold_l+Pplus+dPmp_l*pold_r
  889. c----------------------------
  890.        dpir_r=dPpr_r*pold_l+dPmr_r*pold_r
  891.        dpiu_r=dPpu_r*pold_l+dPmu_r*pold_r
  892.        dpiv_r=dPpv_r*pold_l+dPmv_r*pold_r
  893.        dpiw_r=dPpw_r*pold_l+dPmw_r*pold_r
  894.        dpip_r=dPpp_r*pold_l+Pmin+dPmp_r*pold_r
  895. c---------------------------------------------------------------------
  896. c  !!!!!!!!!!!!!!!!!!  JACOBIAN  !!!!!!!!!!!!!!!!!!!!!!!!!!!!
  897. c---------------------------------------------------------------------
  898.        rum=am*(mpl_m*rold_l+mmin_m*rold_r)+canc*am*termp
  899. c---------------------------------------------------------------------
  900.        rumr_l=damr_l*(mpl_m*rold_l+mmin_m*rold_r)+
  901.      &        am*(d2mr_l*rold_l+mpl_m+d3mr_l*rold_r)
  902.        rumr_l=rumr_l+canc*coef*(damr_l*termp+am*tmpr_l)
  903.        rumu_l=damu_l*(mpl_m*rold_l+mmin_m*rold_r)+
  904.      &        am*(d2mu_l*rold_l+d3mu_l*rold_r)
  905.        rumu_l=rumu_l+canc*coef*(damu_l*termp+am*tmpu_l)
  906.        rumv_l=damv_l*(mpl_m*rold_l+mmin_m*rold_r)+
  907.      &        am*(d2mv_l*rold_l+d3mv_l*rold_r)
  908.        rumv_l=rumv_l+canc*coef*(damv_l*termp+am*tmpv_l)
  909.        rumw_l=damw_l*(mpl_m*rold_l+mmin_m*rold_r)+
  910.      &        am*(d2mw_l*rold_l+d3mw_l*rold_r)
  911.        rumw_l=rumw_l+canc*coef*(damw_l*termp+am*tmpw_l)
  912.        rump_l=damp_l*(mpl_m*rold_l+mmin_m*rold_r)+
  913.      &        am*(d2mp_l*rold_l+d3mp_l*rold_r)
  914.        rump_l=rump_l+canc*coef*(damp_l*termp+am*tmpp_l)
  915. c-------------------------------------------------
  916.        rumr_r=damr_r*(mpl_m*rold_l+mmin_m*rold_r)+
  917.      &        am*(d2mr_r*rold_l+mmin_m+d3mr_r*rold_r)
  918.        rumr_r=rumr_r+canc*coef*(damr_r*termp+am*tmpr_r)
  919.        rumu_r=damu_r*(mpl_m*rold_l+mmin_m*rold_r)+
  920.      &        am*(d2mu_r*rold_l+d3mu_r*rold_r)
  921.        rumu_r=rumu_r+canc*coef*(damu_r*termp+am*tmpu_r)
  922.        rumv_r=damv_r*(mpl_m*rold_l+mmin_m*rold_r)+
  923.      &        am*(d2mv_r*rold_l+d3mv_r*rold_r)
  924.        rumv_r=rumv_r+canc*coef*(damv_r*termp+am*tmpv_r)
  925.        rumw_r=damw_r*(mpl_m*rold_l+mmin_m*rold_r)+
  926.      &        am*(d2mw_r*rold_l+d3mw_r*rold_r)
  927.        rumw_r=rumw_r+canc*coef*(damw_r*termp+am*tmpw_r)
  928.        rump_r=damp_r*(mpl_m*rold_l+mmin_m*rold_r)+
  929.      &        am*(d2mp_r*rold_l+d3mp_r*rold_r)
  930.        rump_r=rump_r+canc*coef*(damp_r*termp+am*tmpp_r)
  931. c---------------------------------------------------------------------
  932. c---------------------------------------------------------------------
  933.        jl(1,1)=rumr_l
  934.        jl(1,2)=rumu_l
  935.        jl(1,3)=rumv_l
  936.        jl(1,4)=rumw_l
  937.        jl(1,5)=rump_l
  938. c------------------------------------
  939.        jr(1,1)=rumr_r
  940.        jr(1,2)=rumu_r
  941.        jr(1,3)=rumv_r
  942.        jr(1,4)=rumw_r
  943.        jr(1,5)=rump_r
  944. c------------------------------------
  945. ccccc    f222222222222  -------------
  946. c derivatives with respect to left set of
  947. c primitive variables
  948. c------------------------------------
  949.        if(rum .ge. 0.0d0) then
  950.          br1=rum*un_l
  951.          br2=rum*ut1_l
  952.          br3=rum*ut2_l
  953.        else
  954.          br1=rum*un_r
  955.          br2=rum*ut1_r
  956.          br3=rum*ut2_r
  957.        endif
  958. c------------------------------------
  959.        if(rum .ge. 0.0d0) then
  960.          br1=rumr_l*un_l
  961.          br2=rumr_l*ut1_l
  962.          br3=rumr_l*ut2_l
  963.        else
  964.          br1=rumr_l*un_r
  965.          br2=rumr_l*ut1_r
  966.          br3=rumr_l*ut2_r
  967.        endif
  968.        jl(2,1)=n_x*(br1+dpir_l)+br2*t1_x+br3*t2_x
  969.        jl(3,1)=n_y*(br1+dpir_l)+br2*t1_y+br3*t2_y
  970.        jl(4,1)=n_z*(br1+dpir_l)+br2*t1_z+br3*t2_z
  971. c-------------------------------------------------------------------------
  972.        if(rum .ge. 0.0d0) then
  973.          br1=rumu_l*un_l+rum*n_x
  974.          br2=rumu_l*ut1_l+rum*t1_x
  975.          br3=rumu_l*ut2_l+rum*t2_x
  976.        else
  977.          br1=rumu_l*un_r
  978.          br2=rumu_l*ut1_r
  979.          br3=rumu_l*ut2_r
  980.        endif
  981.        jl(2,2)=n_x*(br1+dpiu_l)+br2*t1_x+br3*t2_x
  982.        jl(3,2)=n_y*(br1+dpiu_l)+br2*t1_y+br3*t2_y
  983.        jl(4,2)=n_z*(br1+dpiu_l)+br2*t1_z+br3*t2_z
  984. c-------------------------------------------------------------------------
  985.        if(rum .ge. 0.0d0) then
  986.          br1=rumv_l*un_l+rum*n_y
  987.          br2=rumv_l*ut1_l+rum*t1_y
  988.          br3=rumv_l*ut2_l+rum*t2_y
  989.        else
  990.          br1=rumv_l*un_r
  991.          br2=rumv_l*ut1_r
  992.          br3=rumv_l*ut2_r
  993.        endif
  994.        jl(2,3)=n_x*(br1+dpiv_l)+br2*t1_x+br3*t2_x
  995.        jl(3,3)=n_y*(br1+dpiv_l)+br2*t1_y+br3*t2_y
  996.        jl(4,3)=n_z*(br1+dpiv_l)+br2*t1_z+br3*t2_z
  997. c-------------------------------------------------------------------------
  998.        if(rum .ge. 0.0d0) then
  999.          br1=rumw_l*un_l+rum*n_z
  1000.          br2=rumw_l*ut1_l+rum*t1_z
  1001.          br3=rumw_l*ut2_l+rum*t2_z
  1002.        else
  1003.          br1=rumw_l*un_r
  1004.          br2=rumw_l*ut1_r
  1005.          br3=rumw_l*ut2_r
  1006.        endif
  1007.        jl(2,4)=n_x*(br1+dpiw_l)+br2*t1_x+br3*t2_x
  1008.        jl(3,4)=n_y*(br1+dpiw_l)+br2*t1_y+br3*t2_y
  1009.        jl(4,4)=n_z*(br1+dpiw_l)+br2*t1_z+br3*t2_z
  1010. c-------------------------------------------------------------------------
  1011.        if(rum .ge. 0.0d0) then
  1012.          br1=rump_l*un_l
  1013.          br2=rump_l*ut1_l
  1014.          br3=rump_l*ut2_l
  1015.        else
  1016.          br1=rump_l*un_r
  1017.          br2=rump_l*ut1_r
  1018.          br3=rump_l*ut2_r
  1019.        endif
  1020.        jl(2,5)=n_x*(br1+dpip_l)+br2*t1_x+br3*t2_x
  1021.        jl(3,5)=n_y*(br1+dpip_l)+br2*t1_y+br3*t2_y
  1022.        jl(4,5)=n_z*(br1+dpip_l)+br2*t1_z+br3*t2_z
  1023. c-------------------------------------------------------------------------
  1024. c        rrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrr
  1025. c-------------------------------------------------------------------------
  1026.        if(rum .ge. 0.0d0) then
  1027.          br1=rumr_r*un_l
  1028.          br2=rumr_r*ut1_l
  1029.          br3=rumr_r*ut2_l
  1030.        else
  1031.          br1=rumr_r*un_r
  1032.          br2=rumr_r*ut1_r
  1033.          br3=rumr_r*ut2_r
  1034.        endif
  1035.        jr(2,1)=n_x*(br1+dpir_r)+br2*t1_x+br3*t2_x
  1036.        jr(3,1)=n_y*(br1+dpir_r)+br2*t1_y+br3*t2_y
  1037.        jr(4,1)=n_z*(br1+dpir_r)+br2*t1_z+br3*t2_z
  1038. c------------------------------------------------------------------------
  1039.        if(rum .ge. 0.0d0) then
  1040.          br1=rumu_r*un_l
  1041.          br2=rumu_r*ut1_l
  1042.          br3=rumu_r*ut2_l
  1043.        else
  1044.          br1=rumu_r*un_r+rum*n_x
  1045.          br2=rumu_r*ut1_r+rum*t1_x
  1046.          br3=rumu_r*ut2_r+rum*t2_x
  1047.        endif
  1048.        jr(2,2)=n_x*(br1+dpiu_r)+br2*t1_x+br3*t2_x
  1049.        jr(3,2)=n_y*(br1+dpiu_r)+br2*t1_y+br3*t2_y
  1050.        jr(4,2)=n_z*(br1+dpiu_r)+br2*t1_z+br3*t2_z
  1051. c------------------------------------------------------------------------
  1052.        if(rum .ge. 0.0d0) then
  1053.          br1=rumv_r*un_l
  1054.          br2=rumv_r*ut1_l
  1055.          br3=rumv_r*ut2_l
  1056.        else
  1057.          br1=rumv_r*un_r+rum*n_y
  1058.          br2=rumv_r*ut1_r+rum*t1_y
  1059.          br3=rumv_r*ut2_r+rum*t2_y
  1060.        endif
  1061.        jr(2,3)=n_x*(br1+dpiv_r)+br2*t1_x+br3*t2_x
  1062.        jr(3,3)=n_y*(br1+dpiv_r)+br2*t1_y+br3*t2_y
  1063.        jr(4,3)=n_z*(br1+dpiv_r)+br2*t1_z+br3*t2_z
  1064. c-----------------------------------------------------------------------
  1065.        if(rum .ge. 0.0d0) then
  1066.          br1=rumw_r*un_l
  1067.          br2=rumw_r*ut1_l
  1068.          br3=rumw_r*ut2_l
  1069.        else
  1070.          br1=rumw_r*un_r+rum*n_z
  1071.          br2=rumw_r*ut1_r+rum*t1_z
  1072.          br3=rumw_r*ut2_r+rum*t2_z
  1073.        endif
  1074.        jr(2,4)=n_x*(br1+dpiw_r)+br2*t1_x+br3*t2_x
  1075.        jr(3,4)=n_y*(br1+dpiw_r)+br2*t1_y+br3*t2_y
  1076.        jr(4,4)=n_z*(br1+dpiw_r)+br2*t1_z+br3*t2_z
  1077. c-----------------------------------------------------------------------
  1078.        if(rum .ge. 0.0d0) then
  1079.          br1=rump_r*un_l
  1080.          br2=rump_r*ut1_l
  1081.          br3=rump_r*ut2_l
  1082.        else
  1083.          br1=rump_r*un_r
  1084.          br2=rump_r*ut1_r
  1085.          br3=rump_r*ut2_r
  1086.        endif
  1087.        jr(2,5)=n_x*(br1+dpip_r)+br2*t1_x+br3*t2_x
  1088.        jr(3,5)=n_y*(br1+dpip_r)+br2*t1_y+br3*t2_y
  1089.        jr(4,5)=n_z*(br1+dpip_r)+br2*t1_z+br3*t2_z
  1090. c-------------------------------------------------------
  1091. c-------------------------------------------------------
  1092.        hr_l=(eold_l+pold_l/rold_l)
  1093.        hr_r=(eold_r+pold_r/rold_r)
  1094. c-------------------------------------------------
  1095.        if(rum .ge. 0.0d0) then
  1096.          br1=rumr_l*hr_l-rum*ga*pold_l/gm1/rold_l/rold_l
  1097.        else
  1098.          br1=rumr_l*hr_r
  1099.        endif
  1100.        jl(5,1)=br1
  1101. c-------------------------------------------------
  1102.        if(rum .ge. 0.0d0) then
  1103.          br1=rumu_l*hr_l+rum*uold_l
  1104.        else
  1105.          br1=rumu_l*hr_r
  1106.        endif
  1107.        jl(5,2)=br1
  1108. c-------------------------------------------------
  1109.        if(rum .ge. 0.0d0) then
  1110.          br1=rumv_l*hr_l+rum*vold_l
  1111.        else
  1112.          br1=rumv_l*hr_r
  1113.        endif
  1114.        jl(5,3)=br1
  1115. c-------------------------------------------------
  1116.        if(rum .ge. 0.0d0) then
  1117.          br1=rumw_l*hr_l+rum*wold_l
  1118.        else
  1119.          br1=rumw_l*hr_r
  1120.        endif
  1121.        jl(5,4)=br1
  1122. c-------------------------------------------------
  1123.        if(rum .ge. 0.0d0) then
  1124.          br1=rump_l*hr_l+rum*ga/gm1/rold_l
  1125.        else
  1126.          br1=rump_l*hr_r
  1127.        endif
  1128.        jl(5,5)=br1
  1129. c-------------------------------------------------
  1130. c-------------------------------------------------
  1131.        if(rum .ge. 0.0d0) then
  1132.          br1=rumr_r*hr_l
  1133.        else
  1134.          br1=rumr_r*hr_r-rum*ga*pold_r/gm1/rold_r/rold_r
  1135.        endif
  1136.        jr(5,1)=br1
  1137. c---------------------
  1138.        if(rum .ge. 0.0d0) then
  1139.          br1=rumu_r*hr_l
  1140.        else
  1141.          br1=rumu_r*hr_r+rum*uold_r
  1142.        endif
  1143.        jr(5,2)=br1
  1144. c---------------------
  1145.        if(rum .ge. 0.0d0) then
  1146.          br1=rumv_r*hr_l
  1147.        else
  1148.          br1=rumv_r*hr_r+rum*vold_r
  1149.        endif
  1150.        jr(5,3)=br1
  1151. c---------------------
  1152.        if(rum .ge. 0.0d0) then
  1153.          br1=rumw_r*hr_l
  1154.        else
  1155.          br1=rumw_r*hr_r+rum*wold_r
  1156.        endif
  1157.        jr(5,4)=br1
  1158. c---------------------
  1159.        if(rum .ge. 0.0d0) then
  1160.          br1=rump_r*hr_l
  1161.        else
  1162.          br1=rump_r*hr_r+rum*ga/gm1/rold_r
  1163.        endif
  1164.        jr(5,5)=br1
  1165. c---------------------------------------------------
  1166. c----------------------------------------------------------
  1167. c  matrixes wwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwww
  1168. c----------------------------------------------------------
  1169.        wl(1,1)=1.0d0
  1170.        wl(1,2)=0.0d0
  1171.        wl(1,3)=0.0d0
  1172.        wl(1,4)=0.0d0
  1173.        wl(1,5)=0.0d0
  1174. c------------------------------
  1175.        wl(2,1)=-uold_l/rold_l
  1176.        wl(2,2)=1.0d0/rold_l
  1177.        wl(2,3)=0.0d0
  1178.        wl(2,4)=0.0d0
  1179.        wl(2,5)=0.0d0
  1180. c------------------------------
  1181.        wl(3,1)=-vold_l/rold_l
  1182.        wl(3,2)=0.0d0
  1183.        wl(3,3)=1.0d0/rold_l
  1184.        wl(3,4)=0.0d0
  1185.        wl(3,5)=0.0d0
  1186. c------------------------------
  1187.        wl(4,1)=-wold_l/rold_l
  1188.        wl(4,2)=0.0d0
  1189.        wl(4,3)=0.0d0
  1190.        wl(4,4)=1.0d0/rold_l
  1191.        wl(4,5)=0.0d0
  1192. c------------------------------
  1193.        wl(5,1)=gm1*(uold_l*uold_l+vold_l*vold_l+wold_l*wold_l)/2.0d0
  1194.        wl(5,2)=-uold_l*gm1
  1195.        wl(5,3)=-vold_l*gm1
  1196.        wl(5,4)=-wold_l*gm1
  1197.        wl(5,5)=gm1
  1198. c------------------------------
  1199. c------------------------------
  1200.        wr(1,1)=1.0d0
  1201.        wr(1,2)=0.0d0
  1202.        wr(1,3)=0.0d0
  1203.        wr(1,4)=0.0d0
  1204.        wr(1,5)=0.0d0
  1205. c------------------------------
  1206.        wr(2,1)=-uold_r/rold_r
  1207.        wr(2,2)=1.0d0/rold_r
  1208.        wr(2,3)=0.0d0
  1209.        wr(2,4)=0.0d0
  1210.        wr(2,5)=0.0d0
  1211. c------------------------------
  1212.        wr(3,1)=-vold_r/rold_r
  1213.        wr(3,2)=0.0d0
  1214.        wr(3,3)=1.0d0/rold_r
  1215.        wr(3,4)=0.0d0
  1216.        wr(3,5)=0.0d0
  1217. c------------------------------
  1218.        wr(4,1)=-wold_r/rold_r
  1219.        wr(4,2)=0.0d0
  1220.        wr(4,3)=0.0d0
  1221.        wr(4,4)=1.0d0/rold_r
  1222.        wr(4,5)=0.0d0
  1223. c------------------------------
  1224.        wr(5,1)=gm1*(uold_r*uold_r+vold_r*vold_r+wold_r*wold_r)/2.0d0
  1225.        wr(5,2)=-uold_r*gm1
  1226.        wr(5,3)=-vold_r*gm1
  1227.        wr(5,4)=-wold_r*gm1
  1228.        wr(5,5)=gm1
  1229. c----------------------------------------------
  1230. c----------------------------------------------
  1231.        do 1 i=1,5
  1232.         do 2 j=1,5
  1233.          jtl(i,j)=0.0d0
  1234.          jtr(i,j)=0.0d0
  1235.          do 3 k=1,5
  1236.           jtl(i,j)=jtl(i,j)+jl(i,k)*wl(k,j)
  1237.           jtr(i,j)=jtr(i,j)+jr(i,k)*wr(k,j)
  1238.  3       continue
  1239.  2      continue
  1240.  1     continue
  1241. c----------------------------------------------
  1242.        return
  1243.        end
  1244.  
  1245.  
  1246.  
  1247.  
  1248.  
  1249.  
  1250.  
  1251.  
  1252.  
  1253.  
  1254.  
  1255.  
  1256.  
  1257.  
  1258.  
  1259.  
  1260.  
  1261.  
  1262.  
  1263.  
  1264.  

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