Télécharger la plaquette Cast3M


fr en

Télécharger cms3d.eso

Retour à la liste

Numérotation des lignes : 

cms3d

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

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