Télécharger la plaquette Cast3M


fr en

Télécharger conmsp.eso

Retour à la liste

Numérotation des lignes : 

conmsp

  1. C CONMSP    SOURCE    CB215821  20/11/25    13:22:01     10792          
  2.        SUBROUTINE CONMSP(NSP,JLL,JRR,WVEC_L,WVEC_R,NVECT,TVECT,
  3.      &                   mpyn,lrecp,lrecv,nlcg,nlcd)
  4. C************************************************************************
  5. C
  6. C PROJET            :  CASTEM 2000
  7. C
  8. C NOM               :  CONMSP ('convection for multispecies')
  9. C
  10. C DESCRIPTION       :  Voir KONMSP (appele par KONMSP)
  11. C
  12. C LANGAGE           :  FORTRAN 77 + ESOPE 2000 (avec estensions CISI)
  13. C
  14. C AUTEUR            :  S. KUDRIAKOV, DM2S/SFME/LTMF
  15. C
  16. C************************************************************************
  17. C
  18. c----------------------------------------------------------------------
  19. c GENERAL DESCRIPTION:
  20. c   This subroutine provides the jacobians which are the derivatives
  21. c   of the numerical flux function defined at the cell interface
  22. c   with respect to the conservative variables of the left and right
  23. c   cells (relative to the cell interface).
  24. c
  25. c EQUATIONS: 2D Euler equations of gas dynamics - MULTISPECIES GAS
  26. c
  27. c
  28. c REFERENCE: JCP, 129, 364-382 (1996)
  29. c            " A Sequel to AUSM: AUSM+ ".
  30. c----------------------------------------------------------------------
  31. c INPUT:
  32. c
  33. c      alpha   -- parameter of the AUSM+ scheme in the Pressure function;
  34. c                 ( -3/4 <= alpha <= 3/16 ) (imposed as a parameter)
  35. c
  36. c      beta    -- parameter of the AUSM+ scheme in the Mach function;
  37. c                 ( -1/16 <= beta <= 1/2 )  (imposed as a parameter)
  38. c
  39. c      nsp     -- number of species (total);
  40. c
  41. c      wvec_l  -- vector of the primitive variables
  42. c                 (rho,u_x,u_y,p) at the left cell;
  43. c
  44. c      wvec_r  -- vector of the primitive variables
  45. c                 (rho,u_x,u_y,p) at the right cell;
  46. c
  47. c      nvect   -- normal vector to the interface (2 components in 2D);
  48. c
  49. c      tvect   -- tangential vector to the interface;
  50. c
  51. c      mpyn    -- pointer to the vectors of the primitive variables
  52. c                 (Y_1,Y_2,...Y_(nsp-1)) at the left and the right cells;
  53. c
  54. c      lrecp   -- pointer to the vector of specific heats at constant pressure
  55. c                 (size of the vector is equal to number of species (nsp));
  56. c
  57. c      lrecv   -- pointer to the vector of specific heats at constant volume
  58. c                 (size of the vector is equal to number of species (nsp));
  59. c
  60. c      nlcg    -- "local" number corresponding to the left cell;
  61. c
  62. c      nlcd    -- "local" number corresponding to the right cell;
  63. c----------------------------------------------------------------------
  64. c
  65. c OUTPUT:
  66. c
  67. c      jll     -- jakobian matrix (3+nsp) by (3+nsp) -
  68. c                 derivatives of the numerical
  69. c                 flux function with respect to the conservative variables
  70. c                 from the left cell;
  71. c
  72. c      jrr     -- jakobian matrix (3+nsp) by (3+nsp) -
  73. c                 derivatives of the numerical
  74. c                 flux function with respect to the conservative variables
  75. c                 from the right cell.
  76. c----------------------------------------------------------------------
  77.       IMPLICIT INTEGER(I-N)
  78.        integer nsp,jll,jrr,lrecp,lrecv,nlcg,nlcd
  79.        real*8 wvec_l(4),wvec_r(4)
  80.        real*8 nvect(2),tvect(2)
  81.        real*8 alpha,beta
  82.        real*8 gal,gar,gm1l,gm1r,temph
  83.        real*8 n_x,n_y,t_x,t_y
  84.        real*8 un_l, un_r, ut_l, ut_r
  85.        real*8 ml,mr,Mplus,Mmin,mmid
  86.        real*8 mpl_m, mmin_m,am
  87.        real*8 rold_l,uold_l,vold_l,pold_l,eold_l
  88.        real*8 rold_r,uold_r,vold_r,pold_r,eold_r
  89.        real*8 Pplus,Pmin,pmid
  90.        real*8 hr_l,hr_r,det,top,bot
  91.        real*8 br1,br2,br3,br4,temp_l,temp_r,brac_l,brac_r
  92.        real*8 aleft, arigh
  93.        real*8 damr_l,damr_r,damu_l,damu_r
  94.        real*8 damv_l,damv_r,damp_l,damp_r
  95.        real*8 damg_l,damg_r
  96.        real*8 dmlr_l,dmlr_r,dmlu_l,dmlu_r
  97.        real*8 dmlv_l,dmlv_r,dmlp_l,dmlp_r
  98.        real*8 dmrr_l,dmrr_r,dmru_l,dmru_r
  99.        real*8 dmrv_l,dmrv_r,dmrp_l,dmrp_r
  100.        real*8 dMpr_l,dMpr_r,dMpu_l,dMpu_r
  101.        real*8 dMpv_l,dMpv_r,dMpp_l,dMpp_r
  102.        real*8 dMmr_l,dMmr_r,dMmu_l,dMmu_r
  103.        real*8 dMmv_l,dMmv_r,dMmp_l,dMmp_r
  104.        real*8 dmir_l,dmir_r,dmiu_l,dmiu_r
  105.        real*8 dmiv_l,dmiv_r,dmip_l,dmip_r
  106.        real*8 d3mr_l,d3mr_r,d3mu_l,d3mu_r
  107.        real*8 d3mv_l,d3mv_r,d3mp_l,d3mp_r
  108.        real*8 d2mr_l,d2mr_r,d2mu_l,d2mu_r
  109.        real*8 d2mv_l,d2mv_r,d2mp_l,d2mp_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 dPmr_l,dPmr_r,dPmu_l,dPmu_r
  113.        real*8 dPmv_l,dPmv_r,dPmp_l,dPmp_r
  114.        real*8 dpir_l,dpir_r,dpiu_l,dpiu_r
  115.        real*8 dpiv_l,dpiv_r,dpip_l,dpip_r
  116.        integer i,j,k
  117.        parameter(alpha = 0.1875D0, beta = 0.125D0)
  118. C------------------------------------------------------------
  119. -INC SMCHPOI
  120.        POINTEUR MPYN.MPOVAL
  121. C-------------------------------------------------------------
  122. -INC SMLREEL
  123.        POINTEUR MLRECP.MLREEL, MLRECV.MLREEL
  124. C-------------------------------------------------------------
  125. C******* Les fractionines massiques **************************
  126. C-------------------------------------------------------------
  127.        SEGMENT FRAMAS
  128.          REAL*8 YET(NSP)
  129.       ENDSEGMENT
  130.       POINTEUR YL.FRAMAS, YR.FRAMAS
  131. C-------------------------------------------------------
  132. C**********  Les CP's and CV's   ***********************
  133. C-------------------------------------------------------
  134.       SEGMENT GCONST
  135.          REAL*8 GC(NSP)
  136.       ENDSEGMENT
  137.       POINTEUR CP.GCONST, CV.GCONST
  138. C-------------------------------------------------------------
  139. C********  Segments for the elementary matrixes  *************
  140. C-------------------------------------------------------------
  141.        SEGMENT JACEL
  142.          REAL*8 JAC(3+NSP,3+NSP)
  143.        ENDSEGMENT
  144.        POINTEUR JTL.JACEL, JTR.JACEL, JL.JACEL, JR.JACEL,
  145.      &          WL.JACEL, WR.JACEL
  146. c----------------------------------------
  147.        SEGINI JTL
  148.        SEGINI JTR
  149.        SEGINI JL
  150.        SEGINI JR
  151.        SEGINI WL
  152.        SEGINI WR
  153. C-------------------------------------------------------------
  154. C**********  Segments for the vectors  ***********************
  155. C-------------------------------------------------------------
  156.        SEGMENT VECEL
  157.          REAL*8 VV(NSP)
  158.        ENDSEGMENT
  159.        POINTEUR DMLY_L.VECEL, DMLY_R.VECEL,
  160.      &          dmry_l.vecel, dmry_r.vecel,
  161.      &          dMpy_l.vecel, dMpy_r.vecel,
  162.      &          dMmy_l.vecel, dMmy_r.vecel,
  163.      &          dmiy_l.vecel, dmiy_r.vecel,
  164.      &          d3my_l.vecel, d3my_r.vecel,
  165.      &          d2my_l.vecel, d2my_r.vecel,
  166.      &          dPpy_l.vecel, dPpy_r.vecel,
  167.      &          dPmy_l.vecel, dPmy_r.vecel,
  168.      &          dpiy_l.vecel, dpiy_r.vecel,
  169.      &          dgdyl.vecel, dgdyr.vecel
  170. C----------------------------------------------
  171.        SEGINI DMLY_L, DMLY_R,
  172.      &          dmry_l, dmry_r,
  173.      &          dMpy_l, dMpy_r,
  174.      &          dMmy_l, dMmy_r,
  175.      &          dmiy_l, dmiy_r,
  176.      &          d3my_l, d3my_r,
  177.      &          d2my_l, d2my_r,
  178.      &          dPpy_l, dPpy_r,
  179.      &          dPmy_l, dPmy_r,
  180.      &          dpiy_l, dpiy_r,
  181.      &          dgdyl, dgdyr
  182. C-------------------------------------------------------------
  183. C**********  Segments for the flux-vector  *******************
  184. C-------------------------------------------------------------
  185.        SEGMENT FUNEL
  186.          REAL*8 FU(3+NSP)
  187.        ENDSEGMENT
  188.        POINTEUR f.funel
  189. C-------------------------
  190.        SEGINI f
  191. C------------------------------------------------------------
  192.        SEGINI YL, YR
  193.        SEGACT MPYN
  194.        DO 100 I=1,(NSP-1)
  195.            YL.YET(I)=MPYN.VPOCHA(NLCG,I)
  196.            YR.YET(I)=MPYN.VPOCHA(NLCD,I)
  197.  100   CONTINUE
  198. C----------------------------------------
  199.        SEGINI CP, CV
  200.        MLRECP = LRECP
  201.        MLRECV = LRECV
  202.        SEGACT MLRECP, MLRECV
  203.        DO 101 I=1,(NSP-1)
  204.            CP.GC(I)=MLRECP.PROG(I)
  205.            CV.GC(I)=MLRECV.PROG(I)
  206.  101   CONTINUE
  207.        CP.GC(NSP)=MLRECP.PROG(NSP)
  208.        CV.GC(NSP)=MLRECV.PROG(NSP)
  209. c-------------------------------------------------------------
  210. c  Computing GAMMA at the left cell and its derivatives
  211. c  with respect to the primitive variables Y_i
  212. c-------------------------------------------------------------
  213.        top=0.0D0
  214.        bot=0.0D0
  215.        do 40 i=1,(nsp-1)
  216.         top=top+yl.yet(i)*(cp.gc(i)-cp.gc(nsp))
  217.         bot=bot+yl.yet(i)*(cv.gc(i)-cv.gc(nsp))
  218.  40    continue
  219.        top=cp.gc(nsp)+top
  220.        bot=cv.gc(nsp)+bot
  221.        gal=top/bot
  222.        gm1l=gal-1.0d0
  223. c-------------------------------------------------------------
  224.        do 41 i=1,(nsp-1)
  225.         dgdyl.vv(i)=(cp.gc(i)-cp.gc(nsp)-
  226.      &    gal*(cv.gc(i)-cv.gc(nsp)))/bot
  227.  41    continue
  228. c-------------------------------------------------------------
  229. c  Computing GAMMA at the right cell and its derivatives
  230. c  with respect to the primitive variables Y_i
  231. c-------------------------------------------------------------
  232.        top=0.0D0
  233.        bot=0.0D0
  234.        do 42 i=1,(nsp-1)
  235.         top=top+yr.yet(i)*(cp.gc(i)-cp.gc(nsp))
  236.         bot=bot+yr.yet(i)*(cv.gc(i)-cv.gc(nsp))
  237.  42    continue
  238.        top=cp.gc(nsp)+top
  239.        bot=cv.gc(nsp)+bot
  240.        gar=top/bot
  241.        gm1r=gar-1.0d0
  242. c-------------------------------------------------------------
  243.        do 43 i=1,(nsp-1)
  244.         dgdyr.vv(i)=(cp.gc(i)-cp.gc(nsp)-
  245.      &     gar*(cv.gc(i)-cv.gc(nsp)))/bot
  246.  43    continue
  247. c-------------------------------------------------------------
  248.        n_x=nvect(1)
  249.        n_y=nvect(2)
  250.        t_x=tvect(1)
  251.        t_y=tvect(2)
  252. c----------------------------
  253. c----------------------------
  254.        rold_l=wvec_l(1)
  255.        uold_l=wvec_l(2)
  256.        vold_l=wvec_l(3)
  257.        pold_l=wvec_l(4)
  258. c-----------------------
  259.        rold_r=wvec_r(1)
  260.        uold_r=wvec_r(2)
  261.        vold_r=wvec_r(3)
  262.        pold_r=wvec_r(4)
  263. c------------------------------------------------------------------
  264. c Computation of the specific total energy on the left and right.
  265. c------------------------------------------------------------------
  266.        eold_l=(uold_l*uold_l+vold_l*vold_l)/2.0d0
  267.        eold_l=eold_l+pold_l/(gm1l*rold_l)
  268.        eold_r=(uold_r*uold_r+vold_r*vold_r)/2.0d0
  269.        eold_r=eold_r+pold_r/(gm1r*rold_r)
  270. c-------------------------------------------------------------------
  271. c Computation of the speed of sound and its derivatives;
  272. c numerical speed of sound at the interface is taken as an average
  273. c of the speeds of sounds of the neighbouring cells
  274. c-------------------------------------------------------------------
  275.        aleft=sqrt(gal*pold_l/rold_l)
  276.        arigh=sqrt(gar*pold_r/rold_r)
  277.        am=0.5d0*(aleft+arigh)
  278. c-------------------------------------------------------------------
  279.          damr_r=-arigh/(4.0d0*rold_r)
  280.          damu_r=0.0d0
  281.          damv_r=0.0d0
  282.          damp_r=gar/(4.0d0*arigh*rold_r)
  283.          damg_r=arigh/(4.0d0*gar)
  284. c-----------------------
  285.          damr_l=-aleft/(4.0d0*rold_l)
  286.          damu_l=0.0d0
  287.          damv_l=0.0d0
  288.          damp_l=gal/(4.0d0*aleft*rold_l)
  289.          damg_l=aleft/(4.0d0*gal)
  290. c-------------------------------------------------------------------
  291. c Computing numerical Mach number and its derivatives,
  292. c see p.370, under (A1).
  293. c-------------------------------------------------------------------
  294.        un_l=uold_l*n_x+vold_l*n_y
  295.        un_r=uold_r*n_x+vold_r*n_y
  296.        ut_l=uold_l*t_x+vold_l*t_y
  297.        ut_r=uold_r*t_x+vold_r*t_y
  298. c-------------------------------------------------------------------
  299.        ml=un_l/am
  300.        mr=un_r/am
  301. c-------------------------------------------------------------------
  302. c  Mplus and Mmin are calligraphic lettes M+ and M- from the paper,
  303. c  see (19a) and (19b), p.367.
  304. c-------------------------------------------------------------------
  305.        if(abs(ml) .ge. 1.0d0) then
  306.          Mplus=(ml+abs(ml))/2.0d0
  307.        else
  308.          Mplus=(ml+1.0d0)*(ml+1.0d0)/4.0d0
  309.          Mplus=Mplus+beta*(ml*ml-1.0d0)*(ml*ml-1.0d0)
  310.        endif
  311. c-------------------------------------------------------------------
  312.        if(abs(mr) .ge. 1.0d0) then
  313.          Mmin=(mr-abs(mr))/2.0d0
  314.        else
  315.          Mmin=-(mr-1.0d0)*(mr-1.0d0)/4.0d0
  316.          Mmin=Mmin-beta*(mr*mr-1.0d0)*(mr*mr-1.0d0)
  317.        endif
  318. c-------------------------------------------------------------------
  319. c Derivatives of ml and mr with respect to both sets of primitive
  320. c variables: left  and right.
  321. c-------------------------------------------------------------------
  322.        temp_l=-un_l/(am*am)
  323.        temp_r=-un_r/(am*am)
  324. c--------
  325.        dmlr_l=temp_l*damr_l
  326.        dmlr_r=temp_l*damr_r
  327.        dmrr_l=temp_r*damr_l
  328.        dmrr_r=temp_r*damr_r
  329. c--------
  330.        dmlu_l=n_x/am+temp_l*damu_l
  331.        dmlu_r=temp_l*damu_r
  332.        dmru_l=temp_r*damu_l
  333.        dmru_r=n_x/am+temp_r*damu_r
  334. c--------
  335.        dmlv_l=n_y/am+temp_l*damv_l
  336.        dmlv_r=temp_l*damv_r
  337.        dmrv_l=temp_r*damv_l
  338.        dmrv_r=n_y/am+temp_r*damv_r
  339. c--------
  340.        dmlp_l=temp_l*damp_l
  341.        dmlp_r=temp_l*damp_r
  342.        dmrp_l=temp_r*damp_l
  343.        dmrp_r=temp_r*damp_r
  344. c--------
  345.        do 44 i=1,(nsp-1)
  346.          dmly_l.vv(i)=temp_l*damg_l*dgdyl.vv(i)
  347.          dmly_r.vv(i)=temp_l*damg_r*dgdyr.vv(i)
  348.          dmry_l.vv(i)=temp_r*damg_l*dgdyl.vv(i)
  349.          dmry_r.vv(i)=temp_r*damg_r*dgdyr.vv(i)
  350.  44    continue
  351. c-----------------------------------------------------------
  352. c mmid is m_{1/2} (notation as in the paper, see (13),p.366)
  353. c-----------------------------------------------------------
  354.        mmid=Mplus+Mmin
  355. c-----------------------------------------------------------
  356. c Computing the derivatives of M+ and M-
  357. c-----------------------------------------------------------
  358.        if(ml .ge. 1.0d0) then
  359.          dMpr_l=dmlr_l
  360.          dMpu_l=dmlu_l
  361.          dMpv_l=dmlv_l
  362.          dMpp_l=dmlp_l
  363.          do 45 i=1,(nsp-1)
  364.            dMpy_l.vv(i)=dmly_l.vv(i)
  365.  45      continue
  366. c--------------------
  367.          dMpr_r=dmlr_r
  368.          dMpu_r=dmlu_r
  369.          dMpv_r=dmlv_r
  370.          dMpp_r=dmlp_r
  371.          do 46 i=1,(nsp-1)
  372.            dMpy_r.vv(i)=dmly_r.vv(i)
  373.  46      continue
  374.        else
  375.         if((ml .gt. -1.0d0) .and. (ml .lt. 1.0d0)) then
  376.          temph=(ml+1.0d0)/2.0d0
  377.          dMpr_l=(temph+4.0d0*beta*ml*(ml*ml-1.0d0))*dmlr_l
  378.          dMpu_l=(temph+4.0d0*beta*ml*(ml*ml-1.0d0))*dmlu_l
  379.          dMpv_l=(temph+4.0d0*beta*ml*(ml*ml-1.0d0))*dmlv_l
  380.          dMpp_l=(temph+4.0d0*beta*ml*(ml*ml-1.0d0))*dmlp_l
  381.          do 47 i=1,(nsp-1)
  382.          dMpy_l.vv(i)=(temph+4.0d0*beta*ml*
  383.      &      (ml*ml-1.0d0))*dmly_l.vv(i)
  384.  47      continue
  385. c--------------------
  386.            dMpr_r=(temph+4.0d0*beta*ml*(ml*ml-1.0d0))*dmlr_r
  387.            dMpu_r=(temph+4.0d0*beta*ml*(ml*ml-1.0d0))*dmlu_r
  388.            dMpv_r=(temph+4.0d0*beta*ml*(ml*ml-1.0d0))*dmlv_r
  389.            dMpp_r=(temph+4.0d0*beta*ml*(ml*ml-1.0d0))*dmlp_r
  390.            do 48 i=1,(nsp-1)
  391.              dMpy_r.vv(i)=(temph+4.0d0*beta*ml*
  392.      &         (ml*ml-1.0d0))*dmly_r.vv(i)
  393.  48        continue
  394.          else
  395.            dMpr_l=0.0d0
  396.            dMpu_l=0.0d0
  397.            dMpv_l=0.0d0
  398.            dMpp_l=0.0d0
  399.            do 49 i=1,(nsp-1)
  400.              dMpy_l.vv(i)=0.0d0
  401.  49        continue
  402. c---------------------
  403.            dMpr_r=0.0d0
  404.            dMpu_r=0.0d0
  405.            dMpv_r=0.0d0
  406.            dMpp_r=0.0d0
  407.            do 50 i=1,(nsp-1)
  408.              dMpy_r.vv(i)=0.0d0
  409.  50        continue
  410.          endif
  411.        endif
  412. c-----------------------------------------------------------
  413.        if(mr .ge. 1.0d0) then
  414.          dMmr_l=0.0d0
  415.          dMmu_l=0.0d0
  416.          dMmv_l=0.0d0
  417.          dMmp_l=0.0d0
  418.          do 51 i=1,(nsp-1)
  419.            dMmy_l.vv(i)=0.0d0
  420.  51      continue
  421. c---------------------
  422.          dMmr_r=0.0d0
  423.          dMmu_r=0.0d0
  424.          dMmv_r=0.0d0
  425.          dMmp_r=0.0d0
  426.          do 52 i=1,(nsp-1)
  427.            dMmy_r.vv(i)=0.0d0
  428.  52      continue
  429.        else
  430.          if((mr .gt. -1.0d0) .and. (mr .lt. 1.0d0)) then
  431.            temph=(-mr+1.0d0)/2.0d0
  432.            dMmr_l=(temph-4.0d0*beta*mr*(mr*mr-1.0d0))*dmrr_l
  433.            dMmu_l=(temph-4.0d0*beta*mr*(mr*mr-1.0d0))*dmru_l
  434.            dMmv_l=(temph-4.0d0*beta*mr*(mr*mr-1.0d0))*dmrv_l
  435.            dMmp_l=(temph-4.0d0*beta*mr*(mr*mr-1.0d0))*dmrp_l
  436.            do 53 i=1,(nsp-1)
  437.              dMmy_l.vv(i)=(temph-4.0d0*beta*mr*
  438.      &        (mr*mr-1.0d0))*dmry_l.vv(i)
  439.  53        continue
  440. c--------------------
  441.            dMmr_r=(temph-4.0d0*beta*mr*(mr*mr-1.0d0))*dmrr_r
  442.            dMmu_r=(temph-4.0d0*beta*mr*(mr*mr-1.0d0))*dmru_r
  443.            dMmv_r=(temph-4.0d0*beta*mr*(mr*mr-1.0d0))*dmrv_r
  444.            dMmp_r=(temph-4.0d0*beta*mr*(mr*mr-1.0d0))*dmrp_r
  445.            do 54 i=1,(nsp-1)
  446.              dMmy_r.vv(i)=(temph-4.0d0*beta*mr*
  447.      &          (mr*mr-1.0d0))*dmry_r.vv(i)
  448.  54        continue
  449.          else
  450.            dMmr_l=dmrr_l
  451.            dMmu_l=dmru_l
  452.            dMmv_l=dmrv_l
  453.            dMmp_l=dmrp_l
  454.            do 55 i=1,(nsp-1)
  455.              dMmy_l.vv(i)=dmry_l.vv(i)
  456.  55        continue
  457. c--------------------
  458.            dMmr_r=dmrr_r
  459.            dMmu_r=dmru_r
  460.            dMmv_r=dmrv_r
  461.            dMmp_r=dmrp_r
  462.            do 56 i=1,(nsp-1)
  463.              dMmy_r.vv(i)=dmry_r.vv(i)
  464.  56        continue
  465.          endif
  466.        endif
  467. c-----------------------------------------------------------------
  468. c computing the derivatives of m_{1/2} (notation as in the paper)
  469. c-----------------------------------------------------------------
  470.        dmir_l=dMpr_l+dMmr_l
  471.        dmir_r=dMpr_r+dMmr_r
  472. c-------------
  473.        dmiu_l=dMpu_l+dMmu_l
  474.        dmiu_r=dMpu_r+dMmu_r
  475. c-------------
  476.        dmiv_l=dMpv_l+dMmv_l
  477.        dmiv_r=dMpv_r+dMmv_r
  478. c-------------
  479.        dmip_l=dMpp_l+dMmp_l
  480.        dmip_r=dMpp_r+dMmp_r
  481. c-------------
  482.        do 57 i=1,(nsp-1)
  483.          dmiy_l.vv(i)=dMpy_l.vv(i)+dMmy_l.vv(i)
  484.          dmiy_r.vv(i)=dMpy_r.vv(i)+dMmy_r.vv(i)
  485.  57    continue
  486. c----------------------------------------------------------------
  487. c computing the main convective variables and their derivatives
  488. c mpl_m is m^{+}_{1/2} (paper's notation) and
  489. c mmin_m  is m^{-}_{1/2} (paper's notation), see (A2) on p.370.
  490. c----------------------------------------------------------------
  491.        if(mmid .ge. 0.0d0) then
  492.          mpl_m = mmid
  493.          d2mr_l=dmir_l
  494.          d2mu_l=dmiu_l
  495.          d2mv_l=dmiv_l
  496.          d2mp_l=dmip_l
  497.          do 58 i=1,(nsp-1)
  498.            d2my_l.vv(i)=dmiy_l.vv(i)
  499.  58      continue
  500. c------------
  501.          d2mr_r=dmir_r
  502.          d2mu_r=dmiu_r
  503.          d2mv_r=dmiv_r
  504.          d2mp_r=dmip_r
  505.          do 59 i=1,(nsp-1)
  506.            d2my_r.vv(i)=dmiy_r.vv(i)
  507.  59      continue
  508. c------------
  509.        else
  510.          mpl_m = 0.0d0
  511.          d2mr_l=0.0d0
  512.          d2mu_l=0.0d0
  513.          d2mv_l=0.0d0
  514.          d2mp_l=0.0d0
  515.          do 60 i=1,(nsp-1)
  516.            d2my_l.vv(i)=0.0d0
  517.  60      continue
  518. c------------
  519.          d2mr_r=0.0d0
  520.          d2mu_r=0.0d0
  521.          d2mv_r=0.0d0
  522.          d2mp_r=0.0d0
  523.          do 61 i=1,(nsp-1)
  524.            d2my_r.vv(i)=0.0d0
  525.  61      continue
  526.        endif
  527. c---------------------------------------------
  528.        if(mmid .le. 0.0d0) then
  529.          mmin_m = mmid
  530.          d3mr_l=dmir_l
  531.          d3mu_l=dmiu_l
  532.          d3mv_l=dmiv_l
  533.          d3mp_l=dmip_l
  534.          do 62 i=1,(nsp-1)
  535.            d3my_l.vv(i)=dmiy_l.vv(i)
  536.  62      continue
  537. c------------
  538.          d3mr_r=dmir_r
  539.          d3mu_r=dmiu_r
  540.          d3mv_r=dmiv_r
  541.          d3mp_r=dmip_r
  542.          do 63 i=1,(nsp-1)
  543.            d3my_r.vv(i)=dmiy_r.vv(i)
  544.  63      continue
  545. c------------
  546.        else
  547.          mmin_m = 0.0d0
  548.          d3mr_l=0.0d0
  549.          d3mu_l=0.0d0
  550.          d3mv_l=0.0d0
  551.          d3mp_l=0.0d0
  552.          do 64 i=1,(nsp-1)
  553.            d3my_l.vv(i)=0.0d0
  554.  64      continue
  555. c------------
  556.          d3mr_r=0.0d0
  557.          d3mu_r=0.0d0
  558.          d3mv_r=0.0d0
  559.          d3mp_r=0.0d0
  560.          do 65 i=1,(nsp-1)
  561.            d3my_r.vv(i)=0.0d0
  562.  65      continue
  563.        endif
  564. c---------------------------------------------------------------
  565. c Computing the calligraphic P+ and P- with their derivatives,
  566. c see (21a) & (21b) on p.368.
  567. c---------------------------------------------------------------
  568.        if(ml .ge. 1.0d0) then
  569.          Pplus = 1.0d0
  570.        else
  571.          if((ml .gt. -1.0d0) .and. (ml .lt. 1.0d0)) then
  572.            Pplus=(ml+1.0d0)*(ml+1.0d0)*(2.0d0-ml)/4.0d0
  573.            Pplus=Pplus+alpha*ml*(ml*ml-1.0d0)*(ml*ml-1.0d0)
  574.          else
  575.            Pplus = 0.0d0
  576.          endif
  577.        endif
  578. c---------------------------------------------------------------
  579.        if(mr .ge. 1.0d0) then
  580.          Pmin = 0.0d0
  581.        else
  582.          if((mr .gt. -1.0d0) .and. (mr .lt. 1.0d0)) then
  583.            Pmin=(mr-1.0d0)*(mr-1.0d0)*(2.0d0+mr)/4.0d0
  584.            Pmin=Pmin-alpha*mr*(mr*mr-1.0d0)*(mr*mr-1.0d0)
  585.          else
  586.            Pmin = 1.0d0
  587.          endif
  588.        endif
  589. c---------------------------------------------------------------
  590.        brac_l=(ml+1.0d0)*(2.0d0-ml)/2.0d0-(ml+1.0d0)*(ml+1.0d0)/4.0d0
  591.        brac_l=brac_l+alpha*(ml*ml-1.0d0)*(ml*ml-1.0d0)
  592.        brac_l=brac_l+4.0d0*alpha*ml*ml*(ml*ml-1.0d0)
  593. c--------------
  594.        brac_r=(mr-1.0d0)*(2.0d0+mr)/2.0d0+(mr-1.0d0)*(mr-1.0d0)/4.0d0
  595.        brac_r=brac_r-alpha*(mr*mr-1.0d0)*(mr*mr-1.0d0)
  596.        brac_r=brac_r-4.0d0*alpha*mr*mr*(mr*mr-1.0d0)
  597. c---------------------------------------------------------------
  598.        if((ml .gt. -1.0d0) .and. (ml .lt. 1.0d0)) then
  599.          dPpr_l=brac_l*dmlr_l
  600.          dPpr_r=brac_l*dmlr_r
  601. c------------
  602.          dPpu_l=brac_l*dmlu_l
  603.          dPpu_r=brac_l*dmlu_r
  604. c------------
  605.          dPpv_l=brac_l*dmlv_l
  606.          dPpv_r=brac_l*dmlv_r
  607. c------------
  608.          dPpp_l=brac_l*dmlp_l
  609.          dPpp_r=brac_l*dmlp_r
  610. c------------
  611.          do 66 i=1,(nsp-1)
  612.            dPpy_l.vv(i)=brac_l*dmly_l.vv(i)
  613.            dPpy_r.vv(i)=brac_l*dmly_r.vv(i)
  614.  66      continue
  615. c------------
  616.        else
  617.          dPpr_l=0.0d0
  618.          dPpr_r=0.0d0
  619. c-----------
  620.          dPpu_l=0.0d0
  621.          dPpu_r=0.0d0
  622. c-----------
  623.          dPpv_l=0.0d0
  624.          dPpv_r=0.0d0
  625. c-----------
  626.          dPpp_l=0.0d0
  627.          dPpp_r=0.0d0
  628. c-----------
  629.          do 67 i=1,(nsp-1)
  630.            dPpy_l.vv(i)=0.0d0
  631.            dPpy_r.vv(i)=0.0d0
  632.  67      continue
  633. c-----------
  634.        endif
  635. c---------------------------------------------------------------
  636.        if((mr .gt. -1.0d0) .and. (mr .lt. 1.0d0)) then
  637.          dPmr_l=brac_r*dmrr_l
  638.          dPmr_r=brac_r*dmrr_r
  639. c------------
  640.          dPmu_l=brac_r*dmru_l
  641.          dPmu_r=brac_r*dmru_r
  642. c------------
  643.          dPmv_l=brac_r*dmrv_l
  644.          dPmv_r=brac_r*dmrv_r
  645. c------------
  646.          dPmp_l=brac_r*dmrp_l
  647.          dPmp_r=brac_r*dmrp_r
  648. c------------
  649.          do 68 i=1,(nsp-1)
  650.            dPmy_l.vv(i)=brac_r*dmry_l.vv(i)
  651.            dPmy_r.vv(i)=brac_r*dmry_r.vv(i)
  652.  68      continue
  653. c------------
  654.        else
  655.          dPmr_l=0.0d0
  656.          dPmr_r=0.0d0
  657. c-----------
  658.          dPmu_l=0.0d0
  659.          dPmu_r=0.0d0
  660. c-----------
  661.          dPmv_l=0.0d0
  662.          dPmv_r=0.0d0
  663. c-----------
  664.          dPmp_l=0.0d0
  665.          dPmp_r=0.0d0
  666. c-----------
  667.          do 69 i=1,(nsp-1)
  668.            dPmy_l.vv(i)=0.0d0
  669.            dPmy_r.vv(i)=0.0d0
  670.  69      continue
  671. c-----------
  672.        endif
  673. c-------------------------------------------------------------------
  674. c computing pmid -- p_{1/2} and its derivatives, see (20b), p.367.
  675. c-------------------------------------------------------------------
  676.        pmid=Pplus*pold_l + Pmin*pold_r
  677.        dpir_l=dPpr_l*pold_l+dPmr_l*pold_r
  678.        dpiu_l=dPpu_l*pold_l+dPmu_l*pold_r
  679.        dpiv_l=dPpv_l*pold_l+dPmv_l*pold_r
  680.        dpip_l=dPpp_l*pold_l+Pplus+dPmp_l*pold_r
  681.        do 70 i=1,(nsp-1)
  682.          dpiy_l.vv(i)=dPpy_l.vv(i)*pold_l+dPmy_l.vv(i)*pold_r
  683.  70    continue
  684. c----------------------------
  685.        dpir_r=dPpr_r*pold_l+dPmr_r*pold_r
  686.        dpiu_r=dPpu_r*pold_l+dPmu_r*pold_r
  687.        dpiv_r=dPpv_r*pold_l+dPmv_r*pold_r
  688.        dpip_r=dPpp_r*pold_l+Pmin+dPmp_r*pold_r
  689.        do 71 i=1,(nsp-1)
  690.          dpiy_r.vv(i)=dPpy_r.vv(i)*pold_l+dPmy_r.vv(i)*pold_r
  691.  71    continue
  692. c---------------------------------------------------------------------
  693. c computing JACOBIAN as a derivative of the numerical flux function
  694. c with respect to the primitive variables.
  695. c Notation: jl(i,j) --- is the derivative of the i-component of the
  696. c           flux function with respect to the j-component of the
  697. c           vector of primitive variables of the left state.
  698. c           jr(i,j) --- is the derivative of the i-component of the
  699. c           flux function with respect to the j-component of the
  700. c           vector of primitive variables of the right state.
  701. c---------------------------------------------------------------------
  702.        f.fu(1)=am*(mpl_m*rold_l+mmin_m*rold_r)
  703. c---------------------------------------------------------------------
  704.        jl.jac(1,1)=damr_l*f.fu(1)/am+am*(d2mr_l*rold_l+mpl_m)
  705.        jl.jac(1,1)=jl.jac(1,1)+am*d3mr_l*rold_r
  706.        jl.jac(1,2)=damu_l*f.fu(1)/am+am*(d2mu_l*rold_l+d3mu_l*rold_r)
  707.        jl.jac(1,3)=damv_l*f.fu(1)/am+am*(d2mv_l*rold_l+d3mv_l*rold_r)
  708.        jl.jac(1,4)=damp_l*f.fu(1)/am+am*(d2mp_l*rold_l+d3mp_l*rold_r)
  709.        do 72 i=1,(nsp-1)
  710.         jl.jac(1,4+i)=damg_l*dgdyl.vv(i)*f.fu(1)/am+
  711.      &     am*(d2my_l.vv(i)*rold_l+d3my_l.vv(i)*rold_r)
  712.  72    continue
  713. c------------------------------------
  714.        jr.jac(1,1)=damr_r*f.fu(1)/am+am*(d2mr_r*rold_l+mmin_m)
  715.        jr.jac(1,1)=jr.jac(1,1)+am*d3mr_r*rold_r
  716.        jr.jac(1,2)=damu_r*f.fu(1)/am+am*(d2mu_r*rold_l+d3mu_r*rold_r)
  717.        jr.jac(1,3)=damv_r*f.fu(1)/am+am*(d2mv_r*rold_l+d3mv_r*rold_r)
  718.        jr.jac(1,4)=damp_r*f.fu(1)/am+am*(d2mp_r*rold_l+d3mp_r*rold_r)
  719.        do 73 i=1,(nsp-1)
  720.         jr.jac(1,4+i)=damg_r*dgdyr.vv(i)*f.fu(1)/am+
  721.      &     am*(d2my_r.vv(i)*rold_l+d3my_r.vv(i)*rold_r)
  722.  73    continue
  723. c------------------------------------
  724. c------------------------------------
  725.        br1=mpl_m*rold_l*un_l+mmin_m*rold_r*un_r
  726.        br2=mpl_m*rold_l*ut_l+mmin_m*rold_r*ut_r
  727.        f.fu(2)=n_x*(am*br1+pmid)+am*t_x*br2
  728. c------------------
  729.        det=n_x*t_y-n_y*t_x
  730. c---------------------------------------------------------
  731.        br3=d2mr_l*rold_l*un_l+mpl_m*un_l+d3mr_l*rold_r*un_r
  732.        br4=d2mr_l*rold_l*ut_l+mpl_m*ut_l+d3mr_l*rold_r*ut_r
  733.        jl.jac(2,1)=n_x*(damr_l*br1+am*br3+dpir_l)
  734.        jl.jac(2,1)=jl.jac(2,1)+t_x*damr_l*br2+am*t_x*br4
  735. c-------------------
  736.        br3=d2mu_l*rold_l*un_l+mpl_m*rold_l*t_y/det
  737.        br3=br3+d3mu_l*rold_r*un_r
  738.        br4=d2mu_l*rold_l*ut_l+mpl_m*rold_l*(-n_y)/det
  739.        br4=br4+d3mu_l*rold_r*ut_r
  740.        jl.jac(2,2)=n_x*(damu_l*br1+am*br3+dpiu_l)
  741.        jl.jac(2,2)=jl.jac(2,2)+t_x*damu_l*br2+am*t_x*br4
  742. c-------------------
  743.        br3=d2mv_l*rold_l*un_l+mpl_m*rold_l*(-t_x)/det
  744.        br3=br3+d3mv_l*rold_r*un_r
  745.        br4=d2mv_l*rold_l*ut_l+mpl_m*rold_l*n_x/det
  746.        br4=br4+d3mv_l*rold_r*ut_r
  747.        jl.jac(2,3)=n_x*(damv_l*br1+am*br3+dpiv_l)
  748.        jl.jac(2,3)=jl.jac(2,3)+t_x*damv_l*br2+am*t_x*br4
  749. c-------------------
  750.        br3=d2mp_l*rold_l*un_l+d3mp_l*rold_r*un_r
  751.        br4=d2mp_l*rold_l*ut_l+d3mp_l*rold_r*ut_r
  752.        jl.jac(2,4)=n_x*(damp_l*br1+am*br3+dpip_l)
  753.        jl.jac(2,4)=jl.jac(2,4)+t_x*damp_l*br2+am*t_x*br4
  754. c-------------------
  755.        do 74 i=1,(nsp-1)
  756.          br3=d2my_l.vv(i)*rold_l*un_l+d3my_l.vv(i)*rold_r*un_r
  757.          br4=d2my_l.vv(i)*rold_l*ut_l+d3my_l.vv(i)*rold_r*ut_r
  758.          jl.jac(2,4+i)=n_x*(damg_l*dgdyl.vv(i)*br1+
  759.      &       am*br3+dpiy_l.vv(i))
  760.          jl.jac(2,4+i)=jl.jac(2,4+i)+
  761.      &       t_x*damg_l*dgdyl.vv(i)*br2+am*t_x*br4
  762.  74    continue
  763. c-------------------------------------------------------------
  764.        br3=d2mr_r*rold_l*un_l+mmin_m*un_r+d3mr_r*rold_r*un_r
  765.        br4=d2mr_r*rold_l*ut_l+mmin_m*ut_r+d3mr_r*rold_r*ut_r
  766.        jr.jac(2,1)=n_x*(damr_r*br1+am*br3+dpir_r)
  767.        jr.jac(2,1)=jr.jac(2,1)+t_x*damr_r*br2+am*t_x*br4
  768. c-------------------
  769.        br3=d2mu_r*rold_l*un_l+mmin_m*rold_r*t_y/det
  770.        br3=br3+d3mu_r*rold_r*un_r
  771.        br4=d2mu_r*rold_l*ut_l+mmin_m*rold_r*(-n_y/det)
  772.        br4=br4+d3mu_r*rold_r*ut_r
  773.        jr.jac(2,2)=n_x*(damu_r*br1+am*br3+dpiu_r)
  774.        jr.jac(2,2)=jr.jac(2,2)+t_x*damu_r*br2+am*t_x*br4
  775. c-------------------
  776.        br3=d2mv_r*rold_l*un_l+mmin_m*rold_r*(-t_x/det)
  777.        br3=br3+d3mv_r*rold_r*un_r
  778.        br4=d2mv_r*rold_l*ut_l+mmin_m*rold_r*n_x/det
  779.        br4=br4+d3mv_r*rold_r*ut_r
  780.        jr.jac(2,3)=n_x*(damv_r*br1+am*br3+dpiv_r)
  781.        jr.jac(2,3)=jr.jac(2,3)+t_x*damv_r*br2+am*t_x*br4
  782. c-------------------
  783.        br3=d2mp_r*rold_l*un_l+d3mp_r*rold_r*un_r
  784.        br4=d2mp_r*rold_l*ut_l+d3mp_r*rold_r*ut_r
  785.        jr.jac(2,4)=n_x*(damp_r*br1+am*br3+dpip_r)
  786.        jr.jac(2,4)=jr.jac(2,4)+t_x*damp_r*br2+am*t_x*br4
  787. c-------------------
  788.        do 75 i=1,(nsp-1)
  789.          br3=d2my_r.vv(i)*rold_l*un_l+d3my_r.vv(i)*rold_r*un_r
  790.          br4=d2my_r.vv(i)*rold_l*ut_l+d3my_r.vv(i)*rold_r*ut_r
  791.          jr.jac(2,4+i)=n_x*(damg_r*dgdyr.vv(i)*br1+
  792.      &       am*br3+dpiy_r.vv(i))
  793.          jr.jac(2,4+i)=jr.jac(2,4+i)+
  794.      &      t_x*damg_r*dgdyr.vv(i)*br2+am*t_x*br4
  795.  75    continue
  796. c-------------------------------------------------------------
  797. c-------------------------------------------------------------
  798.        br1=mpl_m*rold_l*un_l+mmin_m*rold_r*un_r
  799.        br2=mpl_m*rold_l*ut_l+mmin_m*rold_r*ut_r
  800.        f.fu(3)=n_y*(am*br1+pmid)+am*t_y*br2
  801.        br3=d2mr_l*rold_l*un_l+mpl_m*un_l+d3mr_l*rold_r*un_r
  802.        br4=d2mr_l*rold_l*ut_l+mpl_m*ut_l+d3mr_l*rold_r*ut_r
  803.        jl.jac(3,1)=n_y*(damr_l*br1+am*br3+dpir_l)
  804.        jl.jac(3,1)=jl.jac(3,1)+t_y*damr_l*br2+am*t_y*br4
  805. c-------------------
  806.        br3=d2mu_l*rold_l*un_l+mpl_m*rold_l*t_y/det
  807.        br3=br3+d3mu_l*rold_r*un_r
  808.        br4=d2mu_l*rold_l*ut_l+mpl_m*rold_l*(-n_y/det)
  809.        br4=br4+d3mu_l*rold_r*ut_r
  810.        jl.jac(3,2)=n_y*(damu_l*br1+am*br3+dpiu_l)
  811.        jl.jac(3,2)=jl.jac(3,2)+t_y*damu_l*br2+am*t_y*br4
  812. c-------------------
  813.        br3=d2mv_l*rold_l*un_l+mpl_m*rold_l*(-t_x/det)
  814.        br3=br3+d3mv_l*rold_r*un_r
  815.        br4=d2mv_l*rold_l*ut_l+mpl_m*rold_l*n_x/det
  816.        br4=br4+d3mv_l*rold_r*ut_r
  817.        jl.jac(3,3)=n_y*(damv_l*br1+am*br3+dpiv_l)
  818.        jl.jac(3,3)=jl.jac(3,3)+t_y*damv_l*br2+am*t_y*br4
  819. c-------------------
  820.        br3=d2mp_l*rold_l*un_l+d3mp_l*rold_r*un_r
  821.        br4=d2mp_l*rold_l*ut_l+d3mp_l*rold_r*ut_r
  822.        jl.jac(3,4)=n_y*(damp_l*br1+am*br3+dpip_l)
  823.        jl.jac(3,4)=jl.jac(3,4)+t_y*damp_l*br2+am*t_y*br4
  824. c-------------------
  825.        do 76 i=1,(nsp-1)
  826.          br3=d2my_l.vv(i)*rold_l*un_l+d3my_l.vv(i)*rold_r*un_r
  827.          br4=d2my_l.vv(i)*rold_l*ut_l+d3my_l.vv(i)*rold_r*ut_r
  828.          jl.jac(3,4+i)=n_y*(damg_l*dgdyl.vv(i)*br1+
  829.      &       am*br3+dpiy_l.vv(i))
  830.          jl.jac(3,4+i)=jl.jac(3,4+i)+
  831.      &     t_y*damg_l*dgdyl.vv(i)*br2+am*t_y*br4
  832.  76    continue
  833. c-------------------------------------------------------------
  834.        br3=d2mr_r*rold_l*un_l+mmin_m*un_r+d3mr_r*rold_r*un_r
  835.        br4=d2mr_r*rold_l*ut_l+mmin_m*ut_r+d3mr_r*rold_r*ut_r
  836.        jr.jac(3,1)=n_y*(damr_r*br1+am*br3+dpir_r)
  837.        jr.jac(3,1)=jr.jac(3,1)+t_y*damr_r*br2+am*t_y*br4
  838. c-------------------
  839.        br3=d2mu_r*rold_l*un_l+mmin_m*rold_r*t_y/det
  840.        br3=br3+d3mu_r*rold_r*un_r
  841.        br4=d2mu_r*rold_l*ut_l+mmin_m*rold_r*(-n_y/det)
  842.        br4=br4+d3mu_r*rold_r*ut_r
  843.        jr.jac(3,2)=n_y*(damu_r*br1+am*br3+dpiu_r)
  844.        jr.jac(3,2)=jr.jac(3,2)+t_y*damu_r*br2+am*t_y*br4
  845. c-------------------
  846.        br3=d2mv_r*rold_l*un_l+mmin_m*rold_r*(-t_x/det)
  847.        br3=br3+d3mv_r*rold_r*un_r
  848.        br4=d2mv_r*rold_l*ut_l+mmin_m*rold_r*n_x/det
  849.        br4=br4+d3mv_r*rold_r*ut_r
  850.        jr.jac(3,3)=n_y*(damv_r*br1+am*br3+dpiv_r)
  851.        jr.jac(3,3)=jr.jac(3,3)+t_y*damv_r*br2+am*t_y*br4
  852. c-------------------
  853.        br3=d2mp_r*rold_l*un_l+d3mp_r*rold_r*un_r
  854.        br4=d2mp_r*rold_l*ut_l+d3mp_r*rold_r*ut_r
  855.        jr.jac(3,4)=n_y*(damp_r*br1+am*br3+dpip_r)
  856.        jr.jac(3,4)=jr.jac(3,4)+t_y*damp_r*br2+am*t_y*br4
  857. c-------------------
  858.        do 77 i=1,(nsp-1)
  859.          br3=d2my_r.vv(i)*rold_l*un_l+d3my_r.vv(i)*rold_r*un_r
  860.          br4=d2my_r.vv(i)*rold_l*ut_l+d3my_r.vv(i)*rold_r*ut_r
  861.          jr.jac(3,4+i)=n_y*(damg_r*dgdyr.vv(i)*br1+
  862.      &       am*br3+dpiy_r.vv(i))
  863.          jr.jac(3,4+i)=jr.jac(3,4+i)+
  864.      &      t_y*damg_r*dgdyr.vv(i)*br2+am*t_y*br4
  865.  77    continue
  866. c-------------------------------------------------------------
  867. c  ---------------------  f4 ----------------------------
  868. c-------------------------------------------------------------
  869.        hr_l=rold_l*(uold_l*uold_l+vold_l*vold_l)/2.0d0+gal*pold_l/gm1l
  870.        hr_r=rold_r*(uold_r*uold_r+vold_r*vold_r)/2.0d0+gar*pold_r/gm1r
  871.        f.fu(4)=am*(mpl_m*hr_l+mmin_m*hr_r)
  872. c---------------------
  873.        br1=d2mr_l*hr_l+mpl_m*(uold_l*uold_l+vold_l*vold_l)/2.0d0
  874.        br1=br1+d3mr_l*hr_r
  875.        jl.jac(4,1)=damr_l*f.fu(4)/am+am*br1
  876. c---------------------
  877.        br1=d2mu_l*hr_l+mpl_m*uold_l*rold_l
  878.        br1=br1+d3mu_l*hr_r
  879.        jl.jac(4,2)=damu_l*f.fu(4)/am+am*br1
  880. c---------------------
  881.        br1=d2mv_l*hr_l+mpl_m*vold_l*rold_l
  882.        br1=br1+d3mv_l*hr_r
  883.        jl.jac(4,3)=damv_l*f.fu(4)/am+am*br1
  884. c---------------------
  885.        br1=d2mp_l*hr_l+mpl_m*gal/gm1l
  886.        br1=br1+d3mp_l*hr_r
  887.        jl.jac(4,4)=damp_l*f.fu(4)/am+am*br1
  888. c---------------------
  889.        do 78 i=1,(nsp-1)
  890.          br1=d2my_l.vv(i)*hr_l+mpl_m*(-pold_l)*dgdyl.vv(i)/(gm1l*gm1l)
  891.          br1=br1+d3my_l.vv(i)*hr_r
  892.          jl.jac(4,4+i)=damg_l*dgdyl.vv(i)*f.fu(4)/am+am*br1
  893.  78    continue
  894. c----------------------------------------------------------
  895. c----------------------------------------------------------
  896.        br1=d2mr_r*hr_l+mmin_m*(uold_r*uold_r+vold_r*vold_r)/2.0d0
  897.        br1=br1+d3mr_r*hr_r
  898.        jr.jac(4,1)=damr_r*f.fu(4)/am+am*br1
  899. c---------------------
  900.        br1=d2mu_r*hr_l+mmin_m*uold_r*rold_r
  901.        br1=br1+d3mu_r*hr_r
  902.        jr.jac(4,2)=damu_r*f.fu(4)/am+am*br1
  903. c---------------------
  904.        br1=d2mv_r*hr_l+mmin_m*vold_r*rold_r
  905.        br1=br1+d3mv_r*hr_r
  906.        jr.jac(4,3)=damv_r*f.fu(4)/am+am*br1
  907. c---------------------
  908.        br1=d2mp_r*hr_l+mmin_m*gar/gm1r
  909.        br1=br1+d3mp_r*hr_r
  910.        jr.jac(4,4)=damp_r*f.fu(4)/am+am*br1
  911. c---------------------
  912.        do 79 i=1,(nsp-1)
  913.          br1=d2my_r.vv(i)*hr_l+mmin_m*
  914.      &     (-pold_r)*dgdyr.vv(i)/(gm1r*gm1r)
  915.          br1=br1+d3my_r.vv(i)*hr_r
  916.          jr.jac(4,4+i)=damg_r*dgdyr.vv(i)*f.fu(4)/am+am*br1
  917.  79    continue
  918. c-------------------------------------------------------------
  919. c ------------------ f5++ ------------------------------------
  920. c-------------------------------------------------------------
  921.        do 80 i=1,(nsp-1)
  922.        f.fu(4+i)=am*(mpl_m*rold_l*yl.yet(i)+mmin_m*rold_r*yr.yet(i))
  923. c---------------------
  924.        jl.jac(4+i,1)=damr_l*f.fu(4+i)/am+
  925.      &       am*(d2mr_l*rold_l*yl.yet(i)+mpl_m*yl.yet(i))
  926.        jl.jac(4+i,1)=jl.jac(4+i,1)+am*d3mr_l*rold_r*yr.yet(i)
  927.        jl.jac(4+i,2)=damu_l*f.fu(4+i)/am+am*(d2mu_l*rold_l*yl.yet(i)+
  928.      &           d3mu_l*rold_r*yr.yet(i))
  929.        jl.jac(4+i,3)=damv_l*f.fu(4+i)/am+am*(d2mv_l*rold_l*yl.yet(i)+
  930.      &           d3mv_l*rold_r*yr.yet(i))
  931.        jl.jac(4+i,4)=damp_l*f.fu(4+i)/am+am*(d2mp_l*rold_l*yl.yet(i)+
  932.      &           d3mp_l*rold_r*yr.yet(i))
  933.        do 81 j=5,(4+nsp-1)
  934.         if((4+i).eq.j) then
  935.           jl.jac(4+i,j)=damg_l*dgdyl.vv(j-4)*f.fu(4+i)/am+
  936.      &    am*(d2my_l.vv(j-4)*rold_l*yl.yet(i)+mpl_m*rold_l+
  937.      &    d3my_l.vv(j-4)*rold_r*yr.yet(i))
  938.         else
  939.           jl.jac(4+i,j)=damg_l*dgdyl.vv(j-4)*f.fu(4+i)/am+
  940.      &    am*(d2my_l.vv(j-4)*rold_l*yl.yet(i)+
  941.      &    d3my_l.vv(j-4)*rold_r*yr.yet(i))
  942.         endif
  943.  81    continue
  944. c------------------------------------
  945.        jr.jac(4+i,1)=damr_r*f.fu(4+i)/am+
  946.      &       am*(d2mr_r*rold_l*yl.yet(i)+mmin_m*yr.yet(i))
  947.        jr.jac(4+i,1)=jr.jac(4+i,1)+am*d3mr_r*rold_r*yr.yet(i)
  948.        jr.jac(4+i,2)=damu_r*f.fu(4+i)/am+am*(d2mu_r*rold_l*yl.yet(i)+
  949.      &           d3mu_r*rold_r*yr.yet(i))
  950.        jr.jac(4+i,3)=damv_r*f.fu(4+i)/am+am*(d2mv_r*rold_l*yl.yet(i)+
  951.      &           d3mv_r*rold_r*yr.yet(i))
  952.        jr.jac(4+i,4)=damp_r*f.fu(4+i)/am+am*(d2mp_r*rold_l*yl.yet(i)+
  953.      &           d3mp_r*rold_r*yr.yet(i))
  954.        do 82 j=1,(nsp-1)
  955.         if(i.eq.j) then
  956.           jr.jac(4+i,4+j)=damg_r*dgdyr.vv(j)*f.fu(4+i)/am+
  957.      &    am*(d2my_r.vv(j)*rold_l*yl.yet(i)+mmin_m*rold_r+
  958.      &    d3my_r.vv(j)*rold_r*yr.yet(i))
  959.         else
  960.           jr.jac(4+i,4+j)=damg_r*dgdyr.vv(j)*f.fu(4+i)/am+
  961.      &    am*(d2my_r.vv(j)*rold_l*yl.yet(i)+
  962.      &    d3my_r.vv(j)*rold_r*yr.yet(i))
  963.         endif
  964.  82    continue
  965.  80    continue
  966. c-------------------------------------------------------------
  967. c  matrix wl(i,j) represents the derivative of the i-component
  968. c  of the vector of primitive variables of the left state with
  969. c  respect to the j-component of the vector of the conservative
  970. c  variables of the left state.
  971. c
  972. c  Here: (rho, ux, uy, p, Y_1,...,Y_(nsp-1)) -
  973. c                               vector of primitive variables;
  974. c  (rho, rho ux, rho uy, rho e, rho Y_1,..., rho Y_(nsp-1)) -
  975. c                            vector of conservative variables.
  976. c-------------------------------------------------------------
  977.        wl.jac(1,1)=1.0d0
  978.        wl.jac(1,2)=0.0d0
  979.        wl.jac(1,3)=0.0d0
  980.        wl.jac(1,4)=0.0d0
  981.        do 83 i=1,(nsp-1)
  982.          wl.jac(1,4+i)=0.0d0
  983.  83    continue
  984. c------------------------------
  985.        wl.jac(2,1)=-uold_l/rold_l
  986.        wl.jac(2,2)=1.0d0/rold_l
  987.        wl.jac(2,3)=0.0d0
  988.        wl.jac(2,4)=0.0d0
  989.        do 84 i=1,(nsp-1)
  990.          wl.jac(2,4+i)=0.0d0
  991.  84    continue
  992. c------------------------------
  993.        wl.jac(3,1)=-vold_l/rold_l
  994.        wl.jac(3,2)=0.0d0
  995.        wl.jac(3,3)=1.0d0/rold_l
  996.        wl.jac(3,4)=0.0d0
  997.        do 85 i=1,(nsp-1)
  998.          wl.jac(3,4+i)=0.0d0
  999.  85    continue
  1000. c------------------------------
  1001.        br1=0.0d0
  1002.        do 86 i=1,(nsp-1)
  1003.          br1=br1+dgdyl.vv(i)*yl.yet(i)
  1004.  86    continue
  1005.        br1=br1*pold_l/(rold_l*gm1l)
  1006.        wl.jac(4,1)=gm1l*(uold_l*uold_l+vold_l*vold_l)/2.0d0-br1
  1007.        wl.jac(4,2)=-uold_l*gm1l
  1008.        wl.jac(4,3)=-vold_l*gm1l
  1009.        wl.jac(4,4)=gm1l
  1010.        do 87 i=1,(nsp-1)
  1011.          wl.jac(4,4+i)=dgdyl.vv(i)*pold_l/(rold_l*gm1l)
  1012.  87    continue
  1013. c------------------------------
  1014.        do 88 i=1,(nsp-1)
  1015.         do 89 j=1,4
  1016.           wl.jac(4+i,j)=0.0d0
  1017.           if(j.eq.1) wl.jac(4+i,j)=-yl.yet(i)/rold_l
  1018.  89     continue
  1019. c------------
  1020.         do 890 j=5,(4+nsp-1)
  1021.           wl.jac(4+i,j)=0.0d0
  1022.           if(4+i.eq.j)  then
  1023.             wl.jac(4+i,j)=1.0d0/rold_l
  1024.           endif
  1025.  890    continue
  1026.  88    continue
  1027. c------------------------------
  1028. c------------------------------
  1029.        wr.jac(1,1)=1.0d0
  1030.        wr.jac(1,2)=0.0d0
  1031.        wr.jac(1,3)=0.0d0
  1032.        wr.jac(1,4)=0.0d0
  1033.        do 90 i=1,(nsp-1)
  1034.          wr.jac(1,4+i)=0.0d0
  1035.  90    continue
  1036. c------------------------------
  1037.        wr.jac(2,1)=-uold_r/rold_r
  1038.        wr.jac(2,2)=1.0d0/rold_r
  1039.        wr.jac(2,3)=0.0d0
  1040.        wr.jac(2,4)=0.0d0
  1041.        do 91 i=1,(nsp-1)
  1042.          wr.jac(2,4+i)=0.0d0
  1043.  91    continue
  1044. c------------------------------
  1045.        wr.jac(3,1)=-vold_r/rold_r
  1046.        wr.jac(3,2)=0.0d0
  1047.        wr.jac(3,3)=1.0d0/rold_r
  1048.        wr.jac(3,4)=0.0d0
  1049.        do 92 i=1,(nsp-1)
  1050.          wr.jac(3,4+i)=0.0d0
  1051.  92    continue
  1052. c------------------------------
  1053.        br1=0.0d0
  1054.        do 860 i=1,(nsp-1)
  1055.          br1=br1+dgdyr.vv(i)*yr.yet(i)
  1056.  860   continue
  1057.        br1=br1*pold_r/(rold_r*gm1r)
  1058.        wr.jac(4,1)=gm1r*(uold_r*uold_r+vold_r*vold_r)/2.0d0-br1
  1059.        wr.jac(4,2)=-uold_r*gm1r
  1060.        wr.jac(4,3)=-vold_r*gm1r
  1061.        wr.jac(4,4)=gm1r
  1062.        do 93 i=1,(nsp-1)
  1063.          wr.jac(4,4+i)=dgdyr.vv(i)*pold_r/(rold_r*gm1r)
  1064.  93    continue
  1065. c------------------------------
  1066.        do 94 i=1,(nsp-1)
  1067.         do 95 j=1,4
  1068.           wr.jac(4+i,j)=0.0d0
  1069.           if(j.eq.1) wr.jac(4+i,j)=-yr.yet(i)/rold_r
  1070.  95     continue
  1071. c----------------------
  1072.         do 950 j=5,(3+nsp)
  1073.           wr.jac(4+i,j)=0.0d0
  1074.           if(4+i.eq.j) wr.jac(4+i,j)=1.0d0/rold_r
  1075.  950    continue
  1076.  94    continue
  1077. c----------------------------------------------
  1078. c----------------------------------------------
  1079.        do 1 i=1,(3+nsp)
  1080.         do 2 j=1,(3+nsp)
  1081.          jtl.jac(i,j)=0.0d0
  1082.          jtr.jac(i,j)=0.0d0
  1083.          do 3 k=1,(3+nsp)
  1084.           jtl.jac(i,j)=jtl.jac(i,j)+jl.jac(i,k)*wl.jac(k,j)
  1085.           jtr.jac(i,j)=jtr.jac(i,j)+jr.jac(i,k)*wr.jac(k,j)
  1086.  3       continue
  1087.  2      continue
  1088.  1     continue
  1089. c----------------------------------------------------------------------
  1090.        SEGSUP DMLY_L, DMLY_R,
  1091.      &          dmry_l, dmry_r,
  1092.      &          dMpy_l, dMpy_r,
  1093.      &          dMmy_l, dMmy_r,
  1094.      &          dmiy_l, dmiy_r,
  1095.      &          d3my_l, d3my_r,
  1096.      &          d2my_l, d2my_r,
  1097.      &          dPpy_l, dPpy_r,
  1098.      &          dPmy_l, dPmy_r,
  1099.      &          dpiy_l, dpiy_r,
  1100.      &          dgdyl, dgdyr
  1101. C--------------------------------------------
  1102.        SEGSUP f
  1103. C--------------------------------------------
  1104.        SEGSUP JL
  1105.        SEGSUP JR
  1106.        SEGSUP WL
  1107.        SEGSUP WR
  1108. C--------------------------------------------
  1109.        jll=jtl
  1110.        SEGDES JTL
  1111.        jrr=jtr
  1112.        SEGDES JTR
  1113.        SEGDES YL
  1114.        SEGDES YR
  1115.        SEGDES CP
  1116.        SEGDES CV
  1117.        SEGDES MLRECP, MLRECV
  1118. C--------------------------------------------
  1119.        return
  1120.        end
  1121.  
  1122.  
  1123.  
  1124.  
  1125.  
  1126.  
  1127.  
  1128.  
  1129.  
  1130.  
  1131.  
  1132.  
  1133.  
  1134.  
  1135.  
  1136.  
  1137.  
  1138.  
  1139.  
  1140.  
  1141.  
  1142.  
  1143.  
  1144.  
  1145.  
  1146.  
  1147.  
  1148.  
  1149.  
  1150.  
  1151.  
  1152.  
  1153.  
  1154.  
  1155.  
  1156.  

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