Télécharger la plaquette Cast3M


fr en

Télécharger copmsp.eso

Retour à la liste

Numérotation des lignes : 

copmsp

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

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