Télécharger la plaquette Cast3M


fr en

Télécharger jakwbm.eso

Retour à la liste

Numérotation des lignes : 

jakwbm

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

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