Télécharger la plaquette Cast3M


fr en

Télécharger copms2.eso

Retour à la liste

Numérotation des lignes : 

copms2

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

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