Télécharger la plaquette Cast3M


fr en

Télécharger jw3dbm.eso

Retour à la liste

Numérotation des lignes : 

jw3dbm

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

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