Télécharger jbmms2.eso

Retour à la liste

Numérotation des lignes :

  1. C JBMMS2 SOURCE CHAT 05/01/13 00:48:55 5004
  2. SUBROUTINE JBMMS2(NSP,JLL,JRR,WVEC_L,WVEC_R,NVECT,TVECT,
  3. & mpyn,lrecp,lrecv,nlcg,nlcd,v_inf)
  4. C************************************************************************
  5. C
  6. C PROJET : CASTEM 2000
  7. C
  8. C NOM : JACBM
  9. C
  10. C DESCRIPTION : Voir KONJA6
  11. C
  12. C LANGAGE : FORTRAN 77 + ESOPE 2000 (avec estensions CISI)
  13. C
  14. C AUTEUR : S. KUDRIAKOV, DM2S/SFME/LTMF
  15. C
  16. C************************************************************************
  17. C
  18. c----------------------------------------------------------------------
  19. c GENERAL DESCRIPTION:
  20. c This subroutine provides the 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 The low-mach number corrections are made for the flux functions
  25. c
  26. c EQUATIONS: 2D Euler equations of gas dynamics
  27. c
  28. c
  29. c REFERENCE: 1) JCP, 129, 364-382 (1996)
  30. c " A Sequel to AUSM: AUSM+ ".
  31. c M.S.Liou
  32. c 2) AIAA Journal, Sept. 1998
  33. c " Low-Diffusion Flux-Splitting Methods for Flows at All Speeds"
  34. c J.R.Edwards and M.S.Liou
  35. c----------------------------------------------------------------------
  36. c INPUT:
  37. c
  38. c alpha -- parameter of the AUSM+ scheme in the Pressure function;
  39. c ( -3/4 <= alpha <= 3/16 ) (imposed as a parameter)
  40. c
  41. c beta -- parameter of the AUSM+ scheme in the Mach function;
  42. c ( -1/16 <= beta <= 1/2 ) (imposed as a parameter)
  43. c
  44. c wvec_l -- vector of the primitive variables (rho,ux,uy,p) at the
  45. c left cell;
  46. c
  47. c wvec_r -- vector of the primitive variables (rho,ux,uy,p) at the
  48. c right cell;
  49. c
  50. c nvect -- normal vector to the interface (2 components in 2D);
  51. c
  52. c tvect -- tangential vector to the interface;
  53. c
  54. c ga -- ratio of the specific heats (assumed constant)
  55. c
  56. c v_inf -- parameter for reference velocity, see Ref. 2).
  57. c----------------------------------------------------------------------
  58. c
  59. c OUTPUT:
  60. c
  61. c jtl -- jakobian matrix 4 by 4 - derivatives of the numerical
  62. c flux function with respect to the conservative variables
  63. c from the left cell;
  64. c
  65. c jtr -- jakobian matrix 4 by 4 - derivatives of the numerical
  66. c flux function with respect to the conservative variables
  67. c from the right cell.
  68. c----------------------------------------------------------------------
  69. IMPLICIT INTEGER(I-N)
  70. real*8 wvec_l(4),wvec_r(4)
  71. real*8 nvect(2),tvect(2)
  72. real*8 alpha,beta
  73. real*8 gal,gar,gm1l,gm1r,temph
  74. real*8 n_x,n_y,t_x,t_y
  75. real*8 un_l, un_r, ut_l, ut_r
  76. real*8 ml,mr,Mplus,Mmin,mmid
  77. real*8 mpl_m, mmin_m,am
  78. real*8 rold_l,uold_l,vold_l,pold_l,eold_l
  79. real*8 rold_r,uold_r,vold_r,pold_r,eold_r
  80. real*8 Pplus,Pmin
  81. real*8 hr_l,hr_r,top,bot,bots
  82. real*8 br1,br2,br3,br4,temp_l,temp_r,brac_l,brac_r
  83. real*8 aleft, arigh
  84. real*8 damr_l,damr_r,damu_l,damu_r
  85. real*8 damv_l,damv_r,damp_l,damp_r
  86. real*8 damg_l,damg_r
  87. real*8 dmlr_l,dmlr_r,dmlu_l,dmlu_r
  88. real*8 dmlv_l,dmlv_r,dmlp_l,dmlp_r
  89. real*8 dmrr_l,dmrr_r,dmru_l,dmru_r
  90. real*8 dmrv_l,dmrv_r,dmrp_l,dmrp_r
  91. real*8 dMpr_l,dMpr_r,dMpu_l,dMpu_r
  92. real*8 dMpv_l,dMpv_r,dMpp_l,dMpp_r
  93. real*8 dMmr_l,dMmr_r,dMmu_l,dMmu_r
  94. real*8 dMmv_l,dMmv_r,dMmp_l,dMmp_r
  95. real*8 dmir_l,dmir_r,dmiu_l,dmiu_r
  96. real*8 dmiv_l,dmiv_r,dmip_l,dmip_r
  97. real*8 d3mr_l,d3mr_r,d3mu_l,d3mu_r
  98. real*8 d3mv_l,d3mv_r,d3mp_l,d3mp_r
  99. real*8 d2mr_l,d2mr_r,d2mu_l,d2mu_r
  100. real*8 d2mv_l,d2mv_r,d2mp_l,d2mp_r
  101. real*8 dPpr_l,dPpr_r,dPpu_l,dPpu_r
  102. real*8 dPpv_l,dPpv_r,dPpp_l,dPpp_r
  103. real*8 dPmr_l,dPmr_r,dPmu_l,dPmu_r
  104. real*8 dPmv_l,dPmv_r,dPmp_l,dPmp_r
  105. real*8 dpir_l,dpir_r,dpiu_l,dpiu_r
  106. real*8 dpiv_l,dpiv_r,dpip_l,dpip_r
  107. real*8 epsil,qq,amw,Mmin1,Mplus1
  108. real*8 fmid,mlw,mrw,termp
  109. real*8 ur_r,ur_l,urm,mhalf,mhalfr
  110. real*8 durr_l,durr_r,duru_l,duru_r
  111. real*8 durv_l,durv_r,durp_l,durp_r
  112. real*8 dmhr_l,dmhr_r,dmhu_l,dmhu_r
  113. real*8 dmhv_l,dmhv_r,dmhp_l,dmhp_r
  114. real*8 dmfr_l,dmfr_r,dmfu_l,dmfu_r
  115. real*8 dmfv_l,dmfv_r,dmfp_l,dmfp_r
  116. real*8 dfm_mf,dfm_mh
  117. real*8 dfmr_l,dfmr_r,dfmu_l,dfmu_r
  118. real*8 dfmv_l,dfmv_r,dfmp_l,dfmp_r
  119. real*8 m1mr_l,m1mr_r,m1mu_l,m1mu_r
  120. real*8 m1mv_l,m1mv_r,m1mp_l,m1mp_r
  121. real*8 m1pr_l,m1pr_r,m1pu_l,m1pu_r
  122. real*8 m1pv_l,m1pv_r,m1pp_l,m1pp_r
  123. real*8 tmpr_l,tmpr_r,tmpu_l,tmpu_r
  124. real*8 tmpv_l,tmpv_r,tmpp_l,tmpp_r
  125. real*8 coef,canc,v_inf
  126. real*8 sr_l,sr_r,su_l,su_r,sv_l,sv_r,sp_l,sp_r
  127. real*8 rum,rumr_l,rumu_l,rumv_l,rump_l
  128. real*8 rumr_r,rumu_r,rumv_r,rump_r
  129. real*8 dfmg_l,dfmg_r,dmfg_l,dmfg_r
  130. real*8 dmhg_l,dmhg_r,durg_l,durg_r
  131. integer i,j,k,jll,jrr,lrecp,lrecv,nlcd,nlcg,nsp
  132. parameter(alpha = 0.1875D0, beta = 0.125D0)
  133. parameter(epsil = 1.0d0)
  134. canc=1.0d0
  135. c-------------------------------------------------------------
  136. -INC SMCHPOI
  137. POINTEUR MPYN.MPOVAL
  138. C-------------------------------------------------------------
  139. -INC SMLREEL
  140. POINTEUR MLRECP.MLREEL, MLRECV.MLREEL
  141. C-------------------------------------------------------------
  142. C******* Les fractionines massiques **************************
  143. C-------------------------------------------------------------
  144. SEGMENT FRAMAS
  145. REAL*8 YET(NSP)
  146. ENDSEGMENT
  147. POINTEUR YL.FRAMAS, YR.FRAMAS
  148. C-------------------------------------------------------
  149. C********** Les CP's and CV's ***********************
  150. C-------------------------------------------------------
  151. SEGMENT GCONST
  152. REAL*8 GC(NSP)
  153. ENDSEGMENT
  154. POINTEUR CP.GCONST, CV.GCONST
  155. C-------------------------------------------------------------
  156. C******** Segments for the elementary matrixes *************
  157. C-------------------------------------------------------------
  158. SEGMENT JACEL
  159. REAL*8 JAC(3+NSP,3+NSP)
  160. ENDSEGMENT
  161. POINTEUR JTL.JACEL, JTR.JACEL, JL.JACEL, JR.JACEL,
  162. & WL.JACEL, WR.JACEL
  163. c----------------------------------------
  164. SEGINI JTL
  165. SEGINI JTR
  166. SEGINI JL
  167. SEGINI JR
  168. SEGINI WL
  169. SEGINI WR
  170. C-------------------------------------------------------------
  171. C********** Segments for the vectors ***********************
  172. C-------------------------------------------------------------
  173. SEGMENT VECEL
  174. REAL*8 VV(NSP)
  175. ENDSEGMENT
  176. POINTEUR DMLY_L.VECEL, DMLY_R.VECEL,
  177. & dmry_l.vecel, dmry_r.vecel,
  178. & dMpy_l.vecel, dMpy_r.vecel,
  179. & dMmy_l.vecel, dMmy_r.vecel,
  180. & dmiy_l.vecel, dmiy_r.vecel,
  181. & d3my_l.vecel, d3my_r.vecel,
  182. & d2my_l.vecel, d2my_r.vecel,
  183. & dPpy_l.vecel, dPpy_r.vecel,
  184. & dPmy_l.vecel, dPmy_r.vecel,
  185. & dpiy_l.vecel, dpiy_r.vecel,
  186. & dgdyl.vecel, dgdyr.vecel,
  187. & sy_l.vecel, sy_r.vecel,
  188. & m1py_l.vecel, m1py_r.vecel,
  189. & m1my_l.vecel, m1my_r.vecel,
  190. & tmpy_l.vecel, tmpy_r.vecel,
  191. & rumy_l.vecel, rumy_r.vecel
  192. C----------------------------------------------
  193. SEGINI DMLY_L, DMLY_R,
  194. & dmry_l, dmry_r,
  195. & dMpy_l, dMpy_r,
  196. & dMmy_l, dMmy_r,
  197. & dmiy_l, dmiy_r,
  198. & d3my_l, d3my_r,
  199. & d2my_l, d2my_r,
  200. & dPpy_l, dPpy_r,
  201. & dPmy_l, dPmy_r,
  202. & dpiy_l, dpiy_r,
  203. & dgdyl, dgdyr,
  204. & sy_l, sy_r,
  205. & m1py_l, m1py_r,
  206. & m1my_l, m1my_r,
  207. & tmpy_l, tmpy_r,
  208. & rumy_l, rumy_r
  209. C-------------------------------------------------------------
  210. C********** Segments for the flux-vector *******************
  211. C-------------------------------------------------------------
  212. SEGMENT FUNEL
  213. REAL*8 FU(3+NSP)
  214. ENDSEGMENT
  215. POINTEUR f.funel
  216. C-------------------------
  217. SEGINI f
  218. C------------------------------------------------------------
  219. SEGINI YL, YR
  220. SEGACT MPYN
  221. DO 100 I=1,(NSP-1)
  222. YL.YET(I)=MPYN.VPOCHA(NLCG,I)
  223. YR.YET(I)=MPYN.VPOCHA(NLCD,I)
  224. 100 CONTINUE
  225. C----------------------------------------
  226. SEGINI CP, CV
  227. MLRECP = LRECP
  228. MLRECV = LRECV
  229. SEGACT MLRECP, MLRECV
  230. DO 101 I=1,(NSP-1)
  231. CP.GC(I)=MLRECP.PROG(I)
  232. CV.GC(I)=MLRECV.PROG(I)
  233. 101 CONTINUE
  234. CP.GC(NSP)=MLRECP.PROG(NSP)
  235. CV.GC(NSP)=MLRECV.PROG(NSP)
  236. c-------------------------------------------------------------
  237. c Computing GAMMA at the left cell and its derivatives
  238. c with respect to the primitive variables Y_i
  239. c-------------------------------------------------------------
  240. top=0.0D0
  241. bot=0.0D0
  242. do 40 i=1,(nsp-1)
  243. top=top+yl.yet(i)*(cp.gc(i)-cp.gc(nsp))
  244. bot=bot+yl.yet(i)*(cv.gc(i)-cv.gc(nsp))
  245. 40 continue
  246. top=cp.gc(nsp)+top
  247. bot=cv.gc(nsp)+bot
  248. gal=top/bot
  249. gm1l=gal-1.0d0
  250. c-------------------------------------------------------------
  251. do 41 i=1,(nsp-1)
  252. dgdyl.vv(i)=(cp.gc(i)-cp.gc(nsp)-
  253. & gal*(cv.gc(i)-cv.gc(nsp)))/bot
  254. 41 continue
  255. c-------------------------------------------------------------
  256. c Computing GAMMA at the right cell and its derivatives
  257. c with respect to the primitive variables Y_i
  258. c-------------------------------------------------------------
  259. top=0.0D0
  260. bot=0.0D0
  261. do 42 i=1,(nsp-1)
  262. top=top+yr.yet(i)*(cp.gc(i)-cp.gc(nsp))
  263. bot=bot+yr.yet(i)*(cv.gc(i)-cv.gc(nsp))
  264. 42 continue
  265. top=cp.gc(nsp)+top
  266. bot=cv.gc(nsp)+bot
  267. gar=top/bot
  268. gm1r=gar-1.0d0
  269. c-------------------------------------------------------------
  270. do 43 i=1,(nsp-1)
  271. dgdyr.vv(i)=(cp.gc(i)-cp.gc(nsp)-
  272. & gar*(cv.gc(i)-cv.gc(nsp)))/bot
  273. 43 continue
  274. c-------------------------------------------------------------
  275. n_x=nvect(1)
  276. n_y=nvect(2)
  277. t_x=tvect(1)
  278. t_y=tvect(2)
  279. c----------------------------
  280. rold_l=wvec_l(1)
  281. uold_l=wvec_l(2)
  282. vold_l=wvec_l(3)
  283. pold_l=wvec_l(4)
  284. c-----------------------
  285. rold_r=wvec_r(1)
  286. uold_r=wvec_r(2)
  287. vold_r=wvec_r(3)
  288. pold_r=wvec_r(4)
  289. c------------------------------------------------------------------
  290. c Computation of the specific total energy on the left and right.
  291. c------------------------------------------------------------------
  292. eold_l=(uold_l*uold_l+vold_l*vold_l)/2.0d0
  293. eold_l=eold_l+pold_l/(gm1l*rold_l)
  294. eold_r=(uold_r*uold_r+vold_r*vold_r)/2.0d0
  295. eold_r=eold_r+pold_r/(gm1r*rold_r)
  296. c------------------------------------------------------------------
  297. c Computing reference velocity and its derivatives
  298. c------------------------------------------------------------------
  299. aleft=sqrt(gal*pold_l/rold_l)
  300. arigh=sqrt(gar*pold_r/rold_r)
  301. qq=sqrt(uold_l*uold_l+vold_l*vold_l)
  302. if(qq .lt. (epsil*v_inf)) then
  303. ur_l = epsil*v_inf
  304. durr_l=0.0d0
  305. duru_l=0.0d0
  306. durv_l=0.0d0
  307. durp_l=0.0d0
  308. durg_l=0.0d0
  309. else
  310. ur_l=qq
  311. durr_l=0.0d0
  312. duru_l=uold_l/qq
  313. durv_l=vold_l/qq
  314. durp_l=0.0d0
  315. durg_l=0.0d0
  316. endif
  317. c------------------------------
  318. if(ur_l .ge. aleft) then
  319. ur_l=aleft
  320. durr_l=-aleft/(2.0d0*rold_l)
  321. duru_l=0.0d0
  322. durv_l=0.0d0
  323. durp_l=gal/(2.0d0*aleft*rold_l)
  324. durg_l=aleft/(2.0d0*gal)
  325. endif
  326. c------------------------------------------------------------------
  327. qq=sqrt(uold_r*uold_r+vold_r*vold_r)
  328. if(qq .lt. (epsil*v_inf)) then
  329. ur_r = epsil*v_inf
  330. durr_r=0.0d0
  331. duru_r=0.0d0
  332. durv_r=0.0d0
  333. durp_r=0.0d0
  334. durg_r=0.0d0
  335. else
  336. ur_r=qq
  337. durr_r=0.0d0
  338. duru_r=uold_r/qq
  339. durv_r=vold_r/qq
  340. durp_r=0.0d0
  341. durg_r=0.0d0
  342. endif
  343. c---------------------------
  344. if(ur_r .ge. arigh) then
  345. ur_r=arigh
  346. durr_r=-arigh/(2.0d0*rold_r)
  347. duru_r=0.0d0
  348. durv_r=0.0d0
  349. durp_r=gar/(2.0d0*arigh*rold_r)
  350. durg_r=arigh/(2.0d0*gar)
  351. endif
  352. c------------------------------------------------------------------
  353. c-------------------------------------------------------------------
  354. c Reference velocity at the interface is taken as an average
  355. c of the reference velocities of the neighbouring cells
  356. c-------------------------------------------------------------------
  357. urm=0.5d0*(ur_l+ur_r)
  358. durr_l=0.5d0*durr_l
  359. duru_l=0.5d0*duru_l
  360. durv_l=0.5d0*durv_l
  361. durp_l=0.5d0*durp_l
  362. durg_l=0.5d0*durg_l
  363. c-------------------------
  364. durr_r=0.5d0*durr_r
  365. duru_r=0.5d0*duru_r
  366. durv_r=0.5d0*durv_r
  367. durp_r=0.5d0*durp_r
  368. durg_r=0.5d0*durg_r
  369. c-------------------------------------------------------------------
  370. c Computation of the speed of sound and its derivatives;
  371. c numerical speed of sound at the interface is taken as an average
  372. c of the speeds of sounds of the neighbouring cells
  373. c-------------------------------------------------------------------
  374. am=0.5d0*(aleft+arigh)
  375. c-------------------------------------------------------------------
  376. if(abs(urm/am-1.0d0).le.0.000001d0) then
  377. coef=0.0d0
  378. else
  379. coef=1.0d0
  380. endif
  381. c-------------------------------------------------------------------
  382. damr_r=-arigh/(4.0d0*rold_r)
  383. damu_r=0.0d0
  384. damv_r=0.0d0
  385. damp_r=gar/(4.0d0*arigh*rold_r)
  386. damg_r=arigh/(4.0d0*gar)
  387. c-----------------------
  388. damr_l=-aleft/(4.0d0*rold_l)
  389. damu_l=0.0d0
  390. damv_l=0.0d0
  391. damp_l=gal/(4.0d0*aleft*rold_l)
  392. damg_l=aleft/(4.0d0*gal)
  393. c-------------------------------------------------------------------
  394. c Computing numerical Mach number and its derivatives,
  395. c see p.370, under (A1).
  396. c-------------------------------------------------------------------
  397. un_l=uold_l*n_x+vold_l*n_y
  398. un_r=uold_r*n_x+vold_r*n_y
  399. ut_l=uold_l*t_x+vold_l*t_y
  400. ut_r=uold_r*t_x+vold_r*t_y
  401. c-------------------------------------------------------------------
  402. ml=un_l/am
  403. mr=un_r/am
  404. mhalf=0.5d0*(un_l+un_r)/am
  405. c---------------------
  406. top=0.5d0*(un_l+un_r)/(am*am)
  407. dmhr_l=-top*damr_l
  408. dmhu_l=n_x/2.0d0/am-top*damu_l
  409. dmhv_l=n_y/2.0d0/am-top*damv_l
  410. dmhp_l=-top*damp_l
  411. dmhg_l=-top*damg_l
  412. c---------------------
  413. dmhr_r=-top*damr_r
  414. dmhu_r=n_x/2.0d0/am-top*damu_r
  415. dmhv_r=n_y/2.0d0/am-top*damv_r
  416. dmhp_r=-top*damp_r
  417. dmhg_r=-top*damg_r
  418. c--------------------------------
  419. mhalfr=urm/am
  420. c---------------------
  421. top=urm/(am*am)
  422. dmfr_l=durr_l/am-top*damr_l
  423. dmfu_l=duru_l/am-top*damu_l
  424. dmfv_l=durv_l/am-top*damv_l
  425. dmfp_l=durp_l/am-top*damp_l
  426. dmfg_l=durg_l/am-top*damg_l
  427. c---------------------
  428. dmfr_r=durr_r/am-top*damr_r
  429. dmfu_r=duru_r/am-top*damu_r
  430. dmfv_r=durv_r/am-top*damv_r
  431. dmfp_r=durp_r/am-top*damp_r
  432. dmfg_r=durg_r/am-top*damg_r
  433. c-------------------------------------------------------------------
  434. c Scaling function for the speed of sound and its derivatives
  435. c-------------------------------------------------------------------
  436. top=(1.0d0-mhalfr*mhalfr)*(1.0d0-mhalfr*mhalfr)
  437. top=top*mhalf*mhalf+4.0d0*mhalfr*mhalfr
  438. bot=1.0d0+mhalfr*mhalfr
  439. if(abs(canc-0.0d0).lt.0.000001d0) then
  440. fmid=1.0d0
  441. else
  442. fmid=sqrt(top)/bot
  443. endif
  444. c--------------------------
  445. temph=-4.0d0*(1.0d0-mhalfr*mhalfr)*mhalfr
  446. temph=temph*mhalf*mhalf+8.0d0*mhalfr
  447. dfm_mf=temph/(sqrt(top)*2.0d0*bot)
  448. dfm_mf=dfm_mf-sqrt(top)*2.0d0*mhalfr/(bot*bot)
  449. c--------------------------
  450. temph=2.0d0*(1.0d0-mhalfr*mhalfr)*(1.0d0-mhalfr*mhalfr)*mhalf
  451. dfm_mh=temph/(2.0d0*bot*sqrt(top))
  452. c--------------------------
  453. dfmr_l=dfm_mf*dmfr_l+dfm_mh*dmhr_l
  454. dfmu_l=dfm_mf*dmfu_l+dfm_mh*dmhu_l
  455. dfmv_l=dfm_mf*dmfv_l+dfm_mh*dmhv_l
  456. dfmp_l=dfm_mf*dmfp_l+dfm_mh*dmhp_l
  457. dfmg_l=dfm_mf*dmfg_l+dfm_mh*dmhg_l
  458. c--------------------------
  459. dfmr_r=dfm_mf*dmfr_r+dfm_mh*dmhr_r
  460. dfmu_r=dfm_mf*dmfu_r+dfm_mh*dmhu_r
  461. dfmv_r=dfm_mf*dmfv_r+dfm_mh*dmhv_r
  462. dfmp_r=dfm_mf*dmfp_r+dfm_mh*dmhp_r
  463. dfmg_r=dfm_mf*dmfg_r+dfm_mh*dmhg_r
  464. c--------------------------
  465. amw=fmid*am
  466. mlw=un_l/amw
  467. mrw=un_r/amw
  468. c--------------------------
  469. damr_l=canc*coef*dfmr_l*am+fmid*damr_l
  470. damu_l=canc*coef*dfmu_l*am+fmid*damu_l
  471. damv_l=canc*coef*dfmv_l*am+fmid*damv_l
  472. damp_l=canc*coef*dfmp_l*am+fmid*damp_l
  473. damg_l=canc*coef*dfmg_l*am+fmid*damg_l
  474. c--------------------------
  475. damr_r=canc*coef*dfmr_r*am+fmid*damr_r
  476. damu_r=canc*coef*dfmu_r*am+fmid*damu_r
  477. damv_r=canc*coef*dfmv_r*am+fmid*damv_r
  478. damp_r=canc*coef*dfmp_r*am+fmid*damp_r
  479. damg_r=canc*coef*dfmg_r*am+fmid*damg_r
  480. c-/-/-/-/-/-/-/-/-/-/-/-/-/-/-/-/-/
  481. am=amw
  482. c-------------------------------------------------------------------
  483. c Redefinition of the numerical mach numbers
  484. c-------------------------------------------------------------------
  485. if(abs(canc-0.0d0).lt.0.000001d0) then
  486. top=2.0d0
  487. bot=0.0d0
  488. else
  489. top=1.0d0+mhalfr*mhalfr
  490. bot=1.0d0-mhalfr*mhalfr
  491. endif
  492. ml=0.5d0*(top*mlw+bot*mrw)
  493. mr=0.5d0*(top*mrw+bot*mlw)
  494. c-------------------------------------------------------------------
  495. c Mplus and Mmin are calligraphic lettes M+ and M- from the paper,
  496. c see (19a) and (19b), p.367.
  497. c-------------------------------------------------------------------
  498. if(abs(ml) .ge. 1.0d0) then
  499. Mplus=(ml+abs(ml))/2.0d0
  500. else
  501. Mplus=(ml+1.0d0)*(ml+1.0d0)/4.0d0
  502. Mplus=Mplus+beta*(ml*ml-1.0d0)*(ml*ml-1.0d0)
  503. endif
  504. Mplus1=(ml+abs(ml))/2.0d0
  505. c-------------------------------------------------------------------
  506. if(abs(mr) .ge. 1.0d0) then
  507. Mmin=(mr-abs(mr))/2.0d0
  508. else
  509. Mmin=-(mr-1.0d0)*(mr-1.0d0)/4.0d0
  510. Mmin=Mmin-beta*(mr*mr-1.0d0)*(mr*mr-1.0d0)
  511. endif
  512. Mmin1=(mr-abs(mr))/2.0d0
  513. c-------------------------------------------------------------------
  514. c Derivatives of ml and mr with respect to both sets of primitive
  515. c variables: left and right.
  516. c-------------------------------------------------------------------
  517. temp_l=-un_l/(am*am)
  518. temp_r=-un_r/(am*am)
  519. c--------
  520. dmlr_l=temp_l*damr_l
  521. dmlr_r=temp_l*damr_r
  522. dmrr_l=temp_r*damr_l
  523. dmrr_r=temp_r*damr_r
  524. c--------
  525. dmlu_l=n_x/am+temp_l*damu_l
  526. dmlu_r=temp_l*damu_r
  527. dmru_l=temp_r*damu_l
  528. dmru_r=n_x/am+temp_r*damu_r
  529. c--------
  530. dmlv_l=n_y/am+temp_l*damv_l
  531. dmlv_r=temp_l*damv_r
  532. dmrv_l=temp_r*damv_l
  533. dmrv_r=n_y/am+temp_r*damv_r
  534. c--------
  535. dmlp_l=temp_l*damp_l
    dmlp_r=temp_l*damp_r
  536. dmrp_l=temp_r*damp_l
  537. dmrp_r=temp_r*damp_r
  538. c--------
  539. do 44 i=1,(nsp-1)
  540. dmly_l.vv(i)=temp_l*damg_l*dgdyl.vv(i)
  541. dmly_r.vv(i)=temp_l*damg_r*dgdyr.vv(i)
  542. dmry_l.vv(i)=temp_r*damg_l*dgdyl.vv(i)
  543. dmry_r.vv(i)=temp_r*damg_r*dgdyr.vv(i)
  544. 44 continue
  545. c-----------------------------
  546. sr_l=dmlr_l
  547. sr_r=dmlr_r
  548. su_l=dmlu_l
  549. su_r=dmlu_r
  550. sv_l=dmlv_l
  551. sv_r=dmlv_r
  552. sp_l=dmlp_l
  553. sp_r=dmlp_r
  554. do 440 i=1,(nsp-1)
  555. sy_l.vv(i) = dmly_l.vv(i)
  556. sy_r.vv(i) = dmly_r.vv(i)
  557. 440 continue
  558. c-----------------------------------------------------------------
  559. c Redefinition of the derivatives of the ml & mr
  560. c-----------------------------------------------------------------
  561. temp_l=(mlw-mrw)*mhalfr
  562. temp_r=-temp_l
  563. c--------
  564. dmlr_l=0.5d0*(top*dmlr_l+bot*dmrr_l)+canc*coef*temp_l*dmfr_l
  565. dmlu_l=0.5d0*(top*dmlu_l+bot*dmru_l)+canc*coef*temp_l*dmfu_l
  566. dmlv_l=0.5d0*(top*dmlv_l+bot*dmrv_l)+canc*coef*temp_l*dmfv_l
  567. dmlp_l=0.5d0*(top*dmlp_l+bot*dmrp_l)+canc*coef*temp_l*dmfp_l
  568. c--------
  569. dmlr_r=0.5d0*(top*dmlr_r+bot*dmrr_r)+canc*coef*temp_l*dmfr_r
  570. dmlu_r=0.5d0*(top*dmlu_r+bot*dmru_r)+canc*coef*temp_l*dmfu_r
  571. dmlv_r=0.5d0*(top*dmlv_r+bot*dmrv_r)+canc*coef*temp_l*dmfv_r
  572. dmlp_r=0.5d0*(top*dmlp_r+bot*dmrp_r)+canc*coef*temp_l*dmfp_r
  573. c--------
  574. do 441 i=1,(nsp-1)
  575. dmly_l.vv(i)=0.5d0*(top*dmly_l.vv(i)+bot*dmry_l.vv(i))+
  576. & canc*coef*temp_l*dmfg_l*dgdyl.vv(i)
  577. dmly_r.vv(i)=0.5d0*(top*dmly_r.vv(i)+bot*dmry_r.vv(i))+
  578. & canc*coef*temp_l*dmfg_r*dgdyr.vv(i)
  579. 441 continue
  580. c--------
  581. dmrr_l=0.5d0*(top*dmrr_l+bot*sr_l)+canc*coef*temp_r*dmfr_l
  582. dmru_l=0.5d0*(top*dmru_l+bot*su_l)+canc*coef*temp_r*dmfu_l
  583. dmrv_l=0.5d0*(top*dmrv_l+bot*sv_l)+canc*coef*temp_r*dmfv_l
  584. dmrp_l=0.5d0*(top*dmrp_l+bot*sp_l)+canc*coef*temp_r*dmfp_l
  585. c--------
  586. dmrr_r=0.5d0*(top*dmrr_r+bot*sr_r)+canc*coef*temp_r*dmfr_r
  587. dmru_r=0.5d0*(top*dmru_r+bot*su_r)+canc*coef*temp_r*dmfu_r
  588. dmrv_r=0.5d0*(top*dmrv_r+bot*sv_r)+canc*coef*temp_r*dmfv_r
  589. dmrp_r=0.5d0*(top*dmrp_r+bot*sp_r)+canc*coef*temp_r*dmfp_r
  590. do 442 i=1,(nsp-1)
  591. dmry_l.vv(i)=0.5d0*(top*dmry_l.vv(i)+bot*sy_l.vv(i))+
  592. & canc*coef*temp_r*dmfg_l*dgdyl.vv(i)
  593. dmry_r.vv(i)=0.5d0*(top*dmry_r.vv(i)+bot*sy_r.vv(i))+
  594. & canc*coef*temp_r*dmfg_r*dgdyr.vv(i)
  595. 442 continue
  596. c-----------------------------------------------------------
  597. c mmid is m_{1/2} (notation as in the paper, see (13),p.366)
  598. c-----------------------------------------------------------
  599. mmid=Mplus+Mmin
  600. c-----------------------------------------------------------
  601. c Computing the derivatives of M+ and M-
  602. c-----------------------------------------------------------
  603. if(ml .ge. 1.0d0) then
  604. dMpr_l=dmlr_l
  605. dMpu_l=dmlu_l
  606. dMpv_l=dmlv_l
  607. dMpp_l=dmlp_l
  608. do 45 i=1,(nsp-1)
  609. dMpy_l.vv(i)=dmly_l.vv(i)
  610. 45 continue
  611. c--------------------
  612. dMpr_r=dmlr_r
  613. dMpu_r=dmlu_r
  614. dMpv_r=dmlv_r
  615. dMpp_r=dmlp_r
  616. do 46 i=1,(nsp-1)
  617. dMpy_r.vv(i)=dmly_r.vv(i)
  618. 46 continue
  619. else
  620. if((ml .gt. -1.0d0) .and. (ml .lt. 1.0d0)) then
  621. temph=(ml+1.0d0)/2.0d0
  622. dMpr_l=(temph+4.0d0*beta*ml*(ml*ml-1.0d0))*dmlr_l
  623. dMpu_l=(temph+4.0d0*beta*ml*(ml*ml-1.0d0))*dmlu_l
  624. dMpv_l=(temph+4.0d0*beta*ml*(ml*ml-1.0d0))*dmlv_l
  625. dMpp_l=(temph+4.0d0*beta*ml*(ml*ml-1.0d0))*dmlp_l
  626. do 47 i=1,(nsp-1)
  627. dMpy_l.vv(i)=(temph+4.0d0*beta*ml*
  628. & (ml*ml-1.0d0))*dmly_l.vv(i)
  629. 47 continue
  630. c--------------------
  631. dMpr_r=(temph+4.0d0*beta*ml*(ml*ml-1.0d0))*dmlr_r
  632. dMpu_r=(temph+4.0d0*beta*ml*(ml*ml-1.0d0))*dmlu_r
  633. dMpv_r=(temph+4.0d0*beta*ml*(ml*ml-1.0d0))*dmlv_r
  634. dMpp_r=(temph+4.0d0*beta*ml*(ml*ml-1.0d0))*dmlp_r
  635. do 48 i=1,(nsp-1)
  636. dMpy_r.vv(i)=(temph+4.0d0*beta*ml*
  637. & (ml*ml-1.0d0))*dmly_r.vv(i)
  638. 48 continue
  639. else
  640. dMpr_l=0.0d0
  641. dMpu_l=0.0d0
  642. dMpv_l=0.0d0
  643. dMpp_l=0.0d0
  644. do 49 i=1,(nsp-1)
  645. dMpy_l.vv(i)=0.0d0
  646. 49 continue
  647. c---------------------
  648. dMpr_r=0.0d0
  649. dMpu_r=0.0d0
  650. dMpv_r=0.0d0
  651. dMpp_r=0.0d0
  652. do 50 i=1,(nsp-1)
  653. dMpy_r.vv(i)=0.0d0
  654. 50 continue
  655. endif
  656. endif
  657. c-----------------------------------------------------------
  658. c addition of low Mach number
  659. c-----------------------------------------------------------
  660. if(ml .ge. 0.0d0) then
  661. m1pr_l=dmlr_l
  662. m1pu_l=dmlu_l
  663. m1pv_l=dmlv_l
  664. m1pp_l=dmlp_l
  665. do 450 i=1,(nsp-1)
  666. m1py_l.vv(i)=dmly_l.vv(i)
  667. 450 continue
  668. c---------------------
  669. m1pr_r=dmlr_r
  670. m1pu_r=dmlu_r
  671. m1pv_r=dmlv_r
  672. m1pp_r=dmlp_r
  673. do 460 i=1,(nsp-1)
  674. m1py_r.vv(i)=dmly_r.vv(i)
  675. 460 continue
  676. else
  677. m1pr_l=0.0d0
  678. m1pu_l=0.0d0
  679. m1pv_l=0.0d0
  680. m1pp_l=0.0d0
  681. do 451 i=1,(nsp-1)
  682. m1py_l.vv(i)=0.0d0
  683. 451 continue
  684. c--------------------
  685. m1pr_r=0.0d0
  686. m1pu_r=0.0d0
  687. m1pv_r=0.0d0
  688. m1pp_r=0.0d0
  689. do 461 i=1,(nsp-1)
  690. m1py_r.vv(i)=0.0d0
  691. 461 continue
  692. endif
  693. c-----------------------------------------------------------
  694. if(mr .ge. 1.0d0) then
  695. dMmr_l=0.0d0
  696. dMmu_l=0.0d0
  697. dMmv_l=0.0d0
  698. dMmp_l=0.0d0
  699. do 51 i=1,(nsp-1)
  700. dMmy_l.vv(i)=0.0d0
  701. 51 continue
  702. c---------------------
  703. dMmr_r=0.0d0
  704. dMmu_r=0.0d0
  705. dMmv_r=0.0d0
  706. dMmp_r=0.0d0
  707. do 52 i=1,(nsp-1)
  708. dMmy_r.vv(i)=0.0d0
  709. 52 continue
  710. else
  711. if((mr .gt. -1.0d0) .and. (mr .lt. 1.0d0)) then
  712. temph=(-mr+1.0d0)/2.0d0
  713. dMmr_l=(temph-4.0d0*beta*mr*(mr*mr-1.0d0))*dmrr_l
  714. dMmu_l=(temph-4.0d0*beta*mr*(mr*mr-1.0d0))*dmru_l
  715. dMmv_l=(temph-4.0d0*beta*mr*(mr*mr-1.0d0))*dmrv_l
  716. dMmp_l=(temph-4.0d0*beta*mr*(mr*mr-1.0d0))*dmrp_l
  717. do 53 i=1,(nsp-1)
  718. dMmy_l.vv(i)=(temph-4.0d0*beta*mr*
  719. & (mr*mr-1.0d0))*dmry_l.vv(i)
  720. 53 continue
  721. c--------------------
  722. dMmr_r=(temph-4.0d0*beta*mr*(mr*mr-1.0d0))*dmrr_r
  723. dMmu_r=(temph-4.0d0*beta*mr*(mr*mr-1.0d0))*dmru_r
  724. dMmv_r=(temph-4.0d0*beta*mr*(mr*mr-1.0d0))*dmrv_r
  725. dMmp_r=(temph-4.0d0*beta*mr*(mr*mr-1.0d0))*dmrp_r
  726. do 54 i=1,(nsp-1)
  727. dMmy_r.vv(i)=(temph-4.0d0*beta*mr*
  728. & (mr*mr-1.0d0))*dmry_r.vv(i)
  729. 54 continue
  730. else
  731. dMmr_l=dmrr_l
  732. dMmu_l=dmru_l
  733. dMmv_l=dmrv_l
  734. dMmp_l=dmrp_l
  735. do 55 i=1,(nsp-1)
  736. dMmy_l.vv(i)=dmry_l.vv(i)
  737. 55 continue
  738. c--------------------
  739. dMmr_r=dmrr_r
  740. dMmu_r=dmru_r
  741. dMmv_r=dmrv_r
  742. dMmp_r=dmrp_r
  743. do 56 i=1,(nsp-1)
  744. dMmy_r.vv(i)=dmry_r.vv(i)
  745. 56 continue
  746. endif
  747. endif
  748. c-----------------------------------------------------------
  749. c addition of low Mach number
  750. c-----------------------------------------------------------
  751. if(mr .le. 0.0d0) then
  752. m1mr_l=dmrr_l
  753. m1mu_l=dmru_l
  754. m1mv_l=dmrv_l
  755. m1mp_l=dmrp_l
  756. do 550 i=1,(nsp-1)
  757. m1my_l.vv(i)=dmry_l.vv(i)
  758. 550 continue
  759. c---------------------
  760. m1mr_r=dmrr_r
  761. m1mu_r=dmru_r
  762. m1mv_r=dmrv_r
  763. m1mp_r=dmrp_r
  764. do 560 i=1,(nsp-1)
  765. m1my_r.vv(i)=dmry_r.vv(i)
  766. 560 continue
  767. else
  768. m1mr_l=0.0d0
  769. m1mu_l=0.0d0
  770. m1mv_l=0.0d0
  771. m1mp_l=0.0d0
  772. do 551 i=1,(nsp-1)
  773. m1my_l.vv(i)=0.0d0
  774. 551 continue
  775. c--------------------
  776. m1mr_r=0.0d0
  777. m1mu_r=0.0d0
  778. m1mv_r=0.0d0
  779. m1mp_r=0.0d0
  780. do 561 i=1,(nsp-1)
  781. m1my_r.vv(i)=0.0d0
  782. 561 continue
  783. endif
  784. c-----------------------------------------------------------------
  785. c computing the derivatives of m_{1/2} (notation as in the paper)
  786. c-----------------------------------------------------------------
  787. dmir_l=dMpr_l+dMmr_l
  788. dmir_r=dMpr_r+dMmr_r
  789. c-------------
  790. dmiu_l=dMpu_l+dMmu_l
  791. dmiu_r=dMpu_r+dMmu_r
  792. c-------------
  793. dmiv_l=dMpv_l+dMmv_l
  794. dmiv_r=dMpv_r+dMmv_r
  795. c-------------
  796. dmip_l=dMpp_l+dMmp_l
  797. dmip_r=dMpp_r+dMmp_r
  798. c-------------
  799. do 57 i=1,(nsp-1)
  800. dmiy_l.vv(i)=dMpy_l.vv(i)+dMmy_l.vv(i)
  801. dmiy_r.vv(i)=dMpy_r.vv(i)+dMmy_r.vv(i)
  802. 57 continue
  803. c----------------------------------------------------------------
  804. c computing the main convective variables and their derivatives
  805. c mpl_m is m^{+}_{1/2} (paper's notation) and
  806. c mmin_m is m^{-}_{1/2} (paper's notation), see (A2) on p.370.
  807. c----------------------------------------------------------------
  808. termp=(Mmin1-Mmin+Mplus-Mplus1)*(1.0d0/(mhalfr*mhalfr)-1.0d0)
  809. termp=termp*(pold_l-pold_r)/(pold_l/rold_l+pold_r/rold_r)
  810. c-------------------------------------------------------------
  811. c derivatives of the termp
  812. c-------------------------------------------------------------
  813. top=(Mmin1-Mmin+Mplus-Mplus1)
  814. bots=1.0d0/(pold_l/rold_l+pold_r/rold_r)
  815. bot=(pold_l-pold_r)*bots
  816. temph=1.0d0/(mhalfr*mhalfr)-1.0d0
  817. c---------------------------
  818. tmpr_l=(m1mr_l-dMmr_l+dMpr_l-m1pr_l)*bot*temph
  819. tmpr_l=tmpr_l-2.0d0*bot*top*dmfr_l/(mhalfr*mhalfr*mhalfr)
  820. tmpr_l=tmpr_l+temph*top*bot*bots*(pold_l/rold_l/rold_l)
  821. c---------------------------
  822. tmpu_l=(m1mu_l-dMmu_l+dMpu_l-m1pu_l)*bot*temph
  823. tmpu_l=tmpu_l-2.0d0*bot*top*dmfu_l/(mhalfr*mhalfr*mhalfr)
  824. c---------------------------
  825. tmpv_l=(m1mv_l-dMmv_l+dMpv_l-m1pv_l)*bot*temph
  826. tmpv_l=tmpv_l-2.0d0*bot*top*dmfv_l/(mhalfr*mhalfr*mhalfr)
  827. c---------------------------
  828. tmpp_l=(m1mp_l-dMmp_l+dMpp_l-m1pp_l)*bot*temph
  829. tmpp_l=tmpp_l-2.0d0*bot*top*dmfp_l/(mhalfr*mhalfr*mhalfr)
  830. tmpp_l=tmpp_l+temph*top*bots*(1.0d0-bot/rold_l)
  831. c---------------------------
  832. do 570 i=1,(nsp-1)
  833. tmpy_l.vv(i)=(m1my_l.vv(i)-dMmy_l.vv(i)+
  834. & dMpy_l.vv(i)-m1py_l.vv(i))*bot*temph
  835. tmpy_l.vv(i)=tmpy_l.vv(i)-2.0d0*bot*top*
  836. & dmfg_l*dgdyl.vv(i)/(mhalfr*mhalfr*mhalfr)
  837. 570 continue
  838. c------------rrrrrrrr-------
  839. c------------rrrrrrrr-------
  840. tmpr_r=(m1mr_r-dMmr_r+dMpr_r-m1pr_r)*bot*temph
  841. tmpr_r=tmpr_r-2.0d0*bot*top*dmfr_r/(mhalfr*mhalfr*mhalfr)
  842. tmpr_r=tmpr_r+temph*top*bot*bots*(pold_r/rold_r/rold_r)
  843. c---------------------------
  844. tmpu_r=(m1mu_r-dMmu_r+dMpu_r-m1pu_r)*bot*temph
  845. tmpu_r=tmpu_r-2.0d0*bot*top*dmfu_r/(mhalfr*mhalfr*mhalfr)
  846. c---------------------------
  847. tmpv_r=(m1mv_r-dMmv_r+dMpv_r-m1pv_r)*bot*temph
  848. tmpv_r=tmpv_r-2.0d0*bot*top*dmfv_r/(mhalfr*mhalfr*mhalfr)
  849. c---------------------------
  850. tmpp_r=(m1mp_r-dMmp_r+dMpp_r-m1pp_r)*bot*temph
  851. tmpp_r=tmpp_r-2.0d0*bot*top*dmfp_r/(mhalfr*mhalfr*mhalfr)
  852. tmpp_r=tmpp_r-temph*top*bots*(1.0d0+bot/rold_r)
  853. c---------------------------
  854. do 571 i=1,(nsp-1)
  855. tmpy_r.vv(i)=(m1my_r.vv(i)-dMmy_r.vv(i)+
  856. & dMpy_r.vv(i)-m1py_r.vv(i))*bot*temph
  857. tmpy_r.vv(i)=tmpy_r.vv(i)-2.0d0*bot*top*
  858. & dmfg_r*dgdyr.vv(i)/(mhalfr*mhalfr*mhalfr)
  859. 571 continue
  860. c-------------------------------------------------------------
  861. if(mmid .ge. 0.0d0) then
  862. mpl_m = mmid
  863. d2mr_l=dmir_l
  864. d2mu_l=dmiu_l
  865. d2mv_l=dmiv_l
  866. d2mp_l=dmip_l
  867. do 58 i=1,(nsp-1)
  868. d2my_l.vv(i)=dmiy_l.vv(i)
  869. 58 continue
  870. c------------
  871. d2mr_r=dmir_r
  872. d2mu_r=dmiu_r
  873. d2mv_r=dmiv_r
  874. d2mp_r=dmip_r
  875. do 59 i=1,(nsp-1)
  876. d2my_r.vv(i)=dmiy_r.vv(i)
  877. 59 continue
  878. c------------
  879. else
  880. mpl_m = 0.0d0
  881. d2mr_l=0.0d0
  882. d2mu_l=0.0d0
  883. d2mv_l=0.0d0
  884. d2mp_l=0.0d0
  885. do 60 i=1,(nsp-1)
  886. d2my_l.vv(i)=0.0d0
  887. 60 continue
  888. c------------
  889. d2mr_r=0.0d0
  890. d2mu_r=0.0d0
  891. d2mv_r=0.0d0
  892. d2mp_r=0.0d0
  893. do 61 i=1,(nsp-1)
  894. d2my_r.vv(i)=0.0d0
  895. 61 continue
  896. endif
  897. c---------------------------------------------
  898. cc derivatives for the term termm
  899. cc------------------------------------------------------------------
  900. if(mmid .le. 0.0d0) then
  901. mmin_m = mmid
  902. d3mr_l=dmir_l
  903. d3mu_l=dmiu_l
  904. d3mv_l=dmiv_l
  905. d3mp_l=dmip_l
  906. do 62 i=1,(nsp-1)
  907. d3my_l.vv(i)=dmiy_l.vv(i)
  908. 62 continue
  909. c------------
  910. d3mr_r=dmir_r
  911. d3mu_r=dmiu_r
  912. d3mv_r=dmiv_r
  913. d3mp_r=dmip_r
  914. do 63 i=1,(nsp-1)
  915. d3my_r.vv(i)=dmiy_r.vv(i)
  916. 63 continue
  917. c------------
  918. else
  919. mmin_m = 0.0d0
  920. d3mr_l=0.0d0
  921. d3mu_l=0.0d0
  922. d3mv_l=0.0d0
  923. d3mp_l=0.0d0
  924. do 64 i=1,(nsp-1)
  925. d3my_l.vv(i)=0.0d0
  926. 64 continue
  927. c------------
  928. d3mr_r=0.0d0
  929. d3mu_r=0.0d0
  930. d3mv_r=0.0d0
  931. d3mp_r=0.0d0
  932. do 65 i=1,(nsp-1)
  933. d3my_r.vv(i)=0.0d0
  934. 65 continue
  935. endif
  936. c---------------------------------------------------------------
  937. c Computing the calligraphic P+ and P- with their derivatives,
  938. c see (21a) & (21b) on p.368.
  939. c---------------------------------------------------------------
  940. if(ml .ge. 1.0d0) then
  941. Pplus = 1.0d0
  942. else
  943. if((ml .gt. -1.0d0) .and. (ml .lt. 1.0d0)) then
  944. Pplus=(ml+1.0d0)*(ml+1.0d0)*(2.0d0-ml)/4.0d0
  945. Pplus=Pplus+alpha*ml*(ml*ml-1.0d0)*(ml*ml-1.0d0)
  946. else
  947. Pplus = 0.0d0
  948. endif
  949. endif
  950. c---------------------------------------------------------------
  951. if(mr .ge. 1.0d0) then
  952. Pmin = 0.0d0
  953. else
  954. if((mr .gt. -1.0d0) .and. (mr .lt. 1.0d0)) then
  955. Pmin=(mr-1.0d0)*(mr-1.0d0)*(2.0d0+mr)/4.0d0
  956. Pmin=Pmin-alpha*mr*(mr*mr-1.0d0)*(mr*mr-1.0d0)
  957. else
  958. Pmin = 1.0d0
  959. endif
  960. endif
  961. c---------------------------------------------------------------
  962. brac_l=(ml+1.0d0)*(2.0d0-ml)/2.0d0-(ml+1.0d0)*(ml+1.0d0)/4.0d0
  963. brac_l=brac_l+alpha*(ml*ml-1.0d0)*(ml*ml-1.0d0)
  964. brac_l=brac_l+4.0d0*alpha*ml*ml*(ml*ml-1.0d0)
  965. c--------------
  966. brac_r=(mr-1.0d0)*(2.0d0+mr)/2.0d0+(mr-1.0d0)*(mr-1.0d0)/4.0d0
  967. brac_r=brac_r-alpha*(mr*mr-1.0d0)*(mr*mr-1.0d0)
  968. brac_r=brac_r-4.0d0*alpha*mr*mr*(mr*mr-1.0d0)
  969. c---------------------------------------------------------------
  970. if((ml .gt. -1.0d0) .and. (ml .lt. 1.0d0)) then
  971. dPpr_l=brac_l*dmlr_l
  972. dPpr_r=brac_l*dmlr_r
  973. c------------
  974. dPpu_l=brac_l*dmlu_l
  975. dPpu_r=brac_l*dmlu_r
  976. c------------
  977. dPpv_l=brac_l*dmlv_l
  978. dPpv_r=brac_l*dmlv_r
  979. c------------
  980. dPpp_l=brac_l*dmlp_l
  981. dPpp_r=brac_l*dmlp_r
  982. c------------
  983. do 66 i=1,(nsp-1)
  984. dPpy_l.vv(i)=brac_l*dmly_l.vv(i)
  985. dPpy_r.vv(i)=brac_l*dmly_r.vv(i)
  986. 66 continue
  987. c------------
  988. else
  989. dPpr_l=0.0d0
  990. dPpr_r=0.0d0
  991. c-----------
  992. dPpu_l=0.0d0
  993. dPpu_r=0.0d0
  994. c-----------
  995. dPpv_l=0.0d0
  996. dPpv_r=0.0d0
  997. c-----------
  998. dPpp_l=0.0d0
  999. dPpp_r=0.0d0
  1000. c-----------
  1001. do 67 i=1,(nsp-1)
  1002. dPpy_l.vv(i)=0.0d0
  1003. dPpy_r.vv(i)=0.0d0
  1004. 67 continue
  1005. c-----------
  1006. endif
  1007. c---------------------------------------------------------------
  1008. if((mr .gt. -1.0d0) .and. (mr .lt. 1.0d0)) then
  1009. dPmr_l=brac_r*dmrr_l
  1010. dPmr_r=brac_r*dmrr_r
  1011. c------------
  1012. dPmu_l=brac_r*dmru_l
  1013. dPmu_r=brac_r*dmru_r
  1014. c------------
  1015. dPmv_l=brac_r*dmrv_l
  1016. dPmv_r=brac_r*dmrv_r
  1017. c------------
  1018. dPmp_l=brac_r*dmrp_l
  1019. dPmp_r=brac_r*dmrp_r
  1020. c------------
  1021. do 68 i=1,(nsp-1)
  1022. dPmy_l.vv(i)=brac_r*dmry_l.vv(i)
  1023. dPmy_r.vv(i)=brac_r*dmry_r.vv(i)
  1024. 68 continue
  1025. c------------
  1026. else
  1027. dPmr_l=0.0d0
  1028. dPmr_r=0.0d0
  1029. c-----------
  1030. dPmu_l=0.0d0
  1031. dPmu_r=0.0d0
  1032. c-----------
  1033. dPmv_l=0.0d0
  1034. dPmv_r=0.0d0
  1035. c-----------
  1036. dPmp_l=0.0d0
  1037. dPmp_r=0.0d0
  1038. c-----------
  1039. do 69 i=1,(nsp-1)
  1040. dPmy_l.vv(i)=0.0d0
  1041. dPmy_r.vv(i)=0.0d0
  1042. 69 continue
  1043. c-----------
  1044. endif
  1045. c-------------------------------------------------------------------
  1046. c computing pmid -- p_{1/2} and its derivatives, see (20b), p.367.
  1047. c-------------------------------------------------------------------
  1048. dpir_l=dPpr_l*pold_l+dPmr_l*pold_r
  1049. dpiu_l=dPpu_l*pold_l+dPmu_l*pold_r
  1050. dpiv_l=dPpv_l*pold_l+dPmv_l*pold_r
  1051. dpip_l=dPpp_l*pold_l+Pplus+dPmp_l*pold_r
  1052. do 70 i=1,(nsp-1)
  1053. dpiy_l.vv(i)=dPpy_l.vv(i)*pold_l+dPmy_l.vv(i)*pold_r
  1054. 70 continue
  1055. c----------------------------
  1056. dpir_r=dPpr_r*pold_l+dPmr_r*pold_r
  1057. dpiu_r=dPpu_r*pold_l+dPmu_r*pold_r
  1058. dpiv_r=dPpv_r*pold_l+dPmv_r*pold_r
  1059. dpip_r=dPpp_r*pold_l+Pmin+dPmp_r*pold_r
  1060. do 71 i=1,(nsp-1)
  1061. dpiy_r.vv(i)=dPpy_r.vv(i)*pold_l+dPmy_r.vv(i)*pold_r
  1062. 71 continue
  1063. c---------------------------------------------------------------------
  1064. c Computation of the mass flux (rho * u)_1/2
  1065. c---------------------------------------------------------------
  1066. rum=am*(mpl_m*rold_l+mmin_m*rold_r)+canc*am*termp
  1067. c-------------------------------------------------------
  1068. rumr_l=damr_l*(mpl_m*rold_l+mmin_m*rold_r)+
  1069. & am*(d2mr_l*rold_l+mpl_m+d3mr_l*rold_r)
  1070. rumr_l=rumr_l+canc*coef*(damr_l*termp+am*tmpr_l)
  1071. rumu_l=damu_l*(mpl_m*rold_l+mmin_m*rold_r)+
  1072. & am*(d2mu_l*rold_l+d3mu_l*rold_r)
  1073. rumu_l=rumu_l+canc*coef*(damu_l*termp+am*tmpu_l)
  1074. rumv_l=damv_l*(mpl_m*rold_l+mmin_m*rold_r)+
  1075. & am*(d2mv_l*rold_l+d3mv_l*rold_r)
  1076. rumv_l=rumv_l+canc*coef*(damv_l*termp+am*tmpv_l)
  1077. rump_l=damp_l*(mpl_m*rold_l+mmin_m*rold_r)+
  1078. & am*(d2mp_l*rold_l+d3mp_l*rold_r)
  1079. rump_l=rump_l+canc*coef*(damp_l*termp+am*tmpp_l)
  1080. c------------------
  1081. do 710 i=1,(nsp-1)
  1082. rumy_l.vv(i)=damg_l*dgdyl.vv(i)*(mpl_m*rold_l+
  1083. & mmin_m*rold_r)+am*(d2my_l.vv(i)*rold_l+d3my_l.vv(i)*rold_r)
  1084. rumy_l.vv(i)=rumy_l.vv(i)+canc*coef*
  1085. & (damg_l*dgdyl.vv(i)*termp+am*tmpy_l.vv(i))
  1086. 710 continue
  1087. c-------------------------------------------------
  1088. rumr_r=damr_r*(mpl_m*rold_l+mmin_m*rold_r)+
  1089. & am*(d2mr_r*rold_l+mmin_m+d3mr_r*rold_r)
  1090. rumr_r=rumr_r+canc*coef*(damr_r*termp+am*tmpr_r)
  1091. rumu_r=damu_r*(mpl_m*rold_l+mmin_m*rold_r)+
  1092. & am*(d2mu_r*rold_l+d3mu_r*rold_r)
  1093. rumu_r=rumu_r+canc*coef*(damu_r*termp+am*tmpu_r)
  1094. rumv_r=damv_r*(mpl_m*rold_l+mmin_m*rold_r)+
  1095. & am*(d2mv_r*rold_l+d3mv_r*rold_r)
  1096. rumv_r=rumv_r+canc*coef*(damv_r*termp+am*tmpv_r)
  1097. rump_r=damp_r*(mpl_m*rold_l+mmin_m*rold_r)+
  1098. & am*(d2mp_r*rold_l+d3mp_r*rold_r)
  1099. rump_r=rump_r+canc*coef*(damp_r*termp+am*tmpp_r)
  1100. do 711 i=1,(nsp-1)
  1101. rumy_r.vv(i)=damg_r*dgdyr.vv(i)*(mpl_m*rold_l+mmin_m*rold_r)+
  1102. & am*(d2my_r.vv(i)*rold_l+d3my_r.vv(i)*rold_r)
  1103. rumy_r.vv(i)=rumy_r.vv(i)+canc*coef*
  1104. & (damg_r*dgdyr.vv(i)*termp+am*tmpy_r.vv(i))
  1105. 711 continue
  1106. c---------------------------------------------------------------------
  1107. c computing JACOBIAN as a derivative of the numerical flux function
  1108. c with respect to the primitive variables.
  1109. c Notation: jl(i,j) --- is the derivative of the i-component of the
  1110. c flux function with respect to the j-component of the
  1111. c vector of primitive variables of the left state.
  1112. c jr(i,j) --- is the derivative of the i-component of the
  1113. c flux function with respect to the j-component of the
  1114. c vector of primitive variables of the right state.
  1115. c---------------------------------------------------------------------
  1116. jl.jac(1,1)=rumr_l
  1117. jl.jac(1,2)=rumu_l
  1118. jl.jac(1,3)=rumv_l
  1119. jl.jac(1,4)=rump_l
  1120. do 72 i=1,(nsp-1)
  1121. jl.jac(1,4+i)=rumy_l.vv(i)
  1122. 72 continue
  1123. c------------------------------------
  1124. jr.jac(1,1)=rumr_r
  1125. jr.jac(1,2)=rumu_r
  1126. jr.jac(1,3)=rumv_r
  1127. jr.jac(1,4)=rump_r
  1128. do 73 i=1,(nsp-1)
  1129. jr.jac(1,4+i)=rumy_r.vv(i)
  1130. 73 continue
  1131. c------------------------------------
  1132. c---------------------------------------------------------
  1133. if(rum .ge. 0.0d0) then
  1134. br3=rumr_l*un_l
  1135. br4=rumr_l*ut_l
  1136. else
  1137. br3=rumr_l*un_r
  1138. br4=rumr_l*ut_r
  1139. endif
  1140. jl.jac(2,1)=n_x*(br3+dpir_l)+br4*t_x
  1141. jl.jac(3,1)=n_y*(br3+dpir_l)+br4*t_y
  1142. c-------------------
  1143. if(rum .ge. 0.0d0) then
  1144. br3=rumu_l*un_l+rum*n_x
  1145. br4=rumu_l*ut_l+rum*t_x
  1146. else
  1147. br3=rumu_l*un_r
  1148. br4=rumu_l*ut_r
  1149. endif
  1150. jl.jac(2,2)=n_x*(br3+dpiu_l)+br4*t_x
  1151. jl.jac(3,2)=n_y*(br3+dpiu_l)+br4*t_y
  1152. c-------------------
  1153. if(rum .ge. 0.0d0) then
  1154. br3=rumv_l*un_l+rum*n_y
  1155. br4=rumv_l*ut_l+rum*t_y
  1156. else
  1157. br3=rumv_l*un_r
  1158. br4=rumv_l*ut_r
  1159. endif
  1160. jl.jac(2,3)=n_x*(br3+dpiv_l)+br4*t_x
  1161. jl.jac(3,3)=n_y*(br3+dpiv_l)+br4*t_y
  1162. c-------------------
  1163. if(rum .ge. 0.0d0) then
  1164. br3=rump_l*un_l
  1165. br4=rump_l*ut_l
  1166. else
  1167. br3=rump_l*un_r
  1168. br4=rump_l*ut_r
  1169. endif
  1170. jl.jac(2,4)=n_x*(br3+dpip_l)+br4*t_x
  1171. jl.jac(3,4)=n_y*(br3+dpip_l)+br4*t_y
  1172. c-------------------------------------------------------------
  1173. do 74 i=1,(nsp-1)
  1174. if(rum .ge. 0.0d0) then
  1175. br3=rumy_l.vv(i)*un_l
  1176. br4=rumy_l.vv(i)*ut_l
  1177. else
  1178. br3=rumy_l.vv(i)*un_r
  1179. br4=rumy_l.vv(i)*ut_r
  1180. endif
  1181. jl.jac(2,4+i)=n_x*(br3+dpiy_l.vv(i))+br4*t_x
  1182. jl.jac(3,4+i)=n_y*(br3+dpiy_l.vv(i))+br4*t_y
  1183. 74 continue
  1184. c-------------------------------------------------------------
  1185. if(rum .ge. 0.0d0) then
  1186. br3=rumr_r*un_l
  1187. br4=rumr_r*ut_l
  1188. else
  1189. br3=rumr_r*un_r
  1190. br4=rumr_r*ut_r
  1191. endif
  1192. jr.jac(2,1)=n_x*(br3+dpir_r)+br4*t_x
  1193. jr.jac(3,1)=n_y*(br3+dpir_r)+br4*t_y
  1194. c-------------------
  1195. if(rum .ge. 0.0d0) then
  1196. br3=rumu_r*un_l
  1197. br4=rumu_r*ut_l
  1198. else
  1199. br3=rumu_r*un_r+rum*n_x
  1200. br4=rumu_r*ut_r+rum*t_x
  1201. endif
  1202. jr.jac(2,2)=n_x*(br3+dpiu_r)+br4*t_x
  1203. jr.jac(3,2)=n_y*(br3+dpiu_r)+br4*t_y
  1204. c-------------------
  1205. if(rum .ge. 0.0d0) then
  1206. br3=rumv_r*un_l
  1207. br4=rumv_r*ut_l
  1208. else
  1209. br3=rumv_r*un_r+rum*n_y
  1210. br4=rumv_r*ut_r+rum*t_y
  1211. endif
  1212. jr.jac(2,3)=n_x*(br3+dpiv_r)+br4*t_x
  1213. jr.jac(3,3)=n_y*(br3+dpiv_r)+br4*t_y
  1214. c-------------------
  1215. if(rum .ge. 0.0d0) then
  1216. br3=rump_r*un_l
  1217. br4=rump_r*ut_l
  1218. else
  1219. br3=rump_r*un_r
  1220. br4=rump_r*ut_r
  1221. endif
  1222. jr.jac(2,4)=n_x*(br3+dpip_r)+br4*t_x
  1223. jr.jac(3,4)=n_y*(br3+dpip_r)+br4*t_y
  1224. c--------------------
  1225. do 75 i=1,(nsp-1)
  1226. if(rum .ge. 0.0d0) then
  1227. br3=rumy_r.vv(i)*un_l
  1228. br4=rumy_r.vv(i)*ut_l
  1229. else
  1230. br3=rumy_r.vv(i)*un_r
  1231. br4=rumy_r.vv(i)*ut_r
  1232. endif
  1233. jr.jac(2,4+i)=n_x*(br3+dpiy_r.vv(i))+br4*t_x
  1234. jr.jac(3,4+i)=n_y*(br3+dpiy_r.vv(i))+br4*t_y
  1235. 75 continue
  1236. c-------------------------------------------------------------
  1237. c ------ f44444444444444444444444444444 ---------
  1238. c-------------------------------------------------------------
  1239. hr_l=(uold_l*uold_l+vold_l*vold_l)/2.0d0+
  1240. & gal*pold_l/gm1l/rold_l
  1241. hr_r=(uold_r*uold_r+vold_r*vold_r)/2.0d0+
  1242. & gar*pold_r/gm1r/rold_r
  1243. c-------------------------------------------------------------
  1244. if(rum .ge. 0.0d0) then
  1245. br1=rumr_l*hr_l-rum*gal*pold_l/gm1l/rold_l/rold_l
  1246. else
  1247. br1=rumr_l*hr_r
  1248. endif
  1249. jl.jac(4,1)=br1
  1250. c---------------------
  1251. if(rum .ge. 0.0d0) then
  1252. br1=rumu_l*hr_l+rum*uold_l
  1253. else
  1254. br1=rumu_l*hr_r
  1255. endif
  1256. jl.jac(4,2)=br1
  1257. c---------------------
  1258. if(rum .ge. 0.0d0) then
  1259. br1=rumv_l*hr_l+rum*vold_l
  1260. else
  1261. br1=rumv_l*hr_r
  1262. endif
  1263. jl.jac(4,3)=br1
  1264. c---------------------
  1265. if(rum .ge. 0.0d0) then
  1266. br1=rump_l*hr_l+rum*gal/gm1l/rold_l
  1267. else
  1268. br1=rump_l*hr_r
  1269. endif
  1270. jl.jac(4,4)=br1
  1271. c---------------------
  1272. do 76 i=1,(nsp-1)
  1273. if(rum .ge. 0.0d0) then
  1274. br1=rumy_l.vv(i)*hr_l-
  1275. & rum*(pold_l*dgdyl.vv(i))/(rold_l*gm1l*gm1l)
  1276. else
  1277. br1=rumy_l.vv(i)*hr_r
  1278. endif
  1279. jl.jac(4,4+i)=br1
  1280. 76 continue
  1281. c----------------------------------------------------------
  1282. if(rum .ge. 0.0d0) then
  1283. br1=rumr_r*hr_l
  1284. else
  1285. br1=rumr_r*hr_r-rum*gar*pold_r/gm1r/rold_r/rold_r
  1286. endif
  1287. jr.jac(4,1)=br1
  1288. c---------------------
  1289. if(rum .ge. 0.0d0) then
  1290. br1=rumu_r*hr_l
  1291. else
  1292. br1=rumu_r*hr_r+rum*uold_r
  1293. endif
  1294. jr.jac(4,2)=br1
  1295. c---------------------
  1296. if(rum .ge. 0.0d0) then
  1297. br1=rumv_r*hr_l
  1298. else
  1299. br1=rumv_r*hr_r+rum*vold_r
  1300. endif
  1301. jr.jac(4,3)=br1
  1302. c---------------------
  1303. if(rum .ge. 0.0d0) then
  1304. br1=rump_r*hr_l
  1305. else
  1306. br1=rump_r*hr_r+rum*gar/gm1r/rold_r
  1307. endif
  1308. jr.jac(4,4)=br1
  1309. c---------------------
  1310. do 77 i=1,(nsp-1)
  1311. if(rum .ge. 0.0d0) then
  1312. br1=rumy_r.vv(i)*hr_l
  1313. else
  1314. br1=rumy_r.vv(i)*hr_r-
  1315. & rum*(pold_r*dgdyr.vv(i))/(rold_r*gm1r*gm1r)
  1316. endif
  1317. jr.jac(4,4+i)=br1
  1318. 77 continue
  1319. c-------------------------------------------------------------
  1320. c-------------------------------------------------------------
  1321. do 78 i=1,(nsp-1)
  1322. if(rum .ge. 0.0d0) then
  1323. br1=rumr_l*yl.yet(i)
  1324. br2=rumr_r*yl.yet(i)
  1325. else
  1326. br1=rumr_l*yr.yet(i)
  1327. br2=rumr_r*yr.yet(i)
  1328. endif
  1329. jl.jac(4+i,1)=br1
  1330. jr.jac(4+i,1)=br2
  1331. c-------------------------
  1332. if(rum .ge. 0.0d0) then
  1333. br1=rumu_l*yl.yet(i)
  1334. br2=rumu_r*yl.yet(i)
  1335. else
  1336. br1=rumu_l*yr.yet(i)
  1337. br2=rumu_r*yr.yet(i)
  1338. endif
  1339. jl.jac(4+i,2)=br1
  1340. jr.jac(4+i,2)=br2
  1341. c-------------------------
  1342. if(rum .ge. 0.0d0) then
  1343. br1=rumv_l*yl.yet(i)
  1344. br2=rumv_r*yl.yet(i)
  1345. else
  1346. br1=rumv_l*yr.yet(i)
  1347. br2=rumv_r*yr.yet(i)
  1348. endif
  1349. jl.jac(4+i,3)=br1
  1350. jr.jac(4+i,3)=br2
  1351. c-------------------------
  1352. if(rum .ge. 0.0d0) then
  1353. br1=rump_l*yl.yet(i)
  1354. br2=rump_r*yl.yet(i)
  1355. else
  1356. br1=rump_l*yr.yet(i)
  1357. br2=rump_r*yr.yet(i)
  1358. endif
  1359. jl.jac(4+i,4)=br1
  1360. jr.jac(4+i,4)=br2
  1361. c--------------------------
  1362. do 780 j=1,(nsp-1)
  1363. if(rum .ge. 0.0d0) then
  1364. if(i .eq. j) then
  1365. br1=rumy_l.vv(j)*yl.yet(i)+rum
  1366. br2=rumy_r.vv(j)*yl.yet(i)
  1367. else
  1368. br1=rumy_l.vv(j)*yl.yet(i)
  1369. br2=rumy_r.vv(j)*yl.yet(i)
  1370. endif
  1371. c--------------
  1372. else
  1373. if(i .eq. j) then
  1374. br1=rumy_l.vv(j)*yr.yet(i)
  1375. br2=rumy_r.vv(j)*yr.yet(i)+rum
  1376. else
  1377. br1=rumy_l.vv(j)*yr.yet(i)
  1378. br2=rumy_r.vv(j)*yr.yet(i)
  1379. endif
  1380. endif
  1381. jl.jac(4+i,4+j)=br1
  1382. jr.jac(4+i,4+j)=br2
  1383. 780 continue
  1384. 78 continue
  1385. c-------------------------------------------------------------
  1386. c matrix wl(i,j) represents the derivative of the i-component
  1387. c of the vector of primitive variables of the left state with
  1388. c respect to the j-component of the vector of the conservative
  1389. c variables of the left state.
  1390. c
  1391. c Here: (rho, ux, uy, p, Y_1,...,Y_(nsp-1)) -
  1392. c vector of primitive variables;
  1393. c (rho, rho ux, rho uy, rho e, rho Y_1,..., rho Y_(nsp-1)) -
  1394. c vector of conservative variables.
  1395. c-------------------------------------------------------------
  1396. wl.jac(1,1)=1.0d0
  1397. wl.jac(1,2)=0.0d0
  1398. wl.jac(1,3)=0.0d0
  1399. wl.jac(1,4)=0.0d0
  1400. do 83 i=1,(nsp-1)
  1401. wl.jac(1,4+i)=0.0d0
  1402. 83 continue
  1403. c------------------------------
  1404. wl.jac(2,1)=-uold_l/rold_l
  1405. wl.jac(2,2)=1.0d0/rold_l
  1406. wl.jac(2,3)=0.0d0
  1407. wl.jac(2,4)=0.0d0
  1408. do 84 i=1,(nsp-1)
  1409. wl.jac(2,4+i)=0.0d0
  1410. 84 continue
  1411. c------------------------------
  1412. wl.jac(3,1)=-vold_l/rold_l
  1413. wl.jac(3,2)=0.0d0
  1414. wl.jac(3,3)=1.0d0/rold_l
  1415. wl.jac(3,4)=0.0d0
  1416. do 85 i=1,(nsp-1)
  1417. wl.jac(3,4+i)=0.0d0
  1418. 85 continue
  1419. c------------------------------
  1420. br1=0.0d0
  1421. do 86 i=1,(nsp-1)
  1422. br1=br1+dgdyl.vv(i)*yl.yet(i)
  1423. 86 continue
  1424. br1=br1*pold_l/(rold_l*gm1l)
  1425. wl.jac(4,1)=gm1l*(uold_l*uold_l+vold_l*vold_l)/2.0d0-br1
  1426. wl.jac(4,2)=-uold_l*gm1l
  1427. wl.jac(4,3)=-vold_l*gm1l
  1428. wl.jac(4,4)=gm1l
  1429. do 87 i=1,(nsp-1)
  1430. wl.jac(4,4+i)=dgdyl.vv(i)*pold_l/(rold_l*gm1l)
  1431. 87 continue
  1432. c------------------------------
  1433. do 88 i=1,(nsp-1)
  1434. do 89 j=1,4
  1435. wl.jac(4+i,j)=0.0d0
  1436. if(j.eq.1) wl.jac(4+i,j)=-yl.yet(i)/rold_l
  1437. 89 continue
  1438. c------------
  1439. do 890 j=5,(4+nsp-1)
  1440. wl.jac(4+i,j)=0.0d0
  1441. if(4+i.eq.j) then
  1442. wl.jac(4+i,j)=1.0d0/rold_l
  1443. endif
  1444. 890 continue
  1445. 88 continue
  1446. c------------------------------
  1447. c------------------------------
  1448. wr.jac(1,1)=1.0d0
  1449. wr.jac(1,2)=0.0d0
  1450. wr.jac(1,3)=0.0d0
  1451. wr.jac(1,4)=0.0d0
  1452. do 90 i=1,(nsp-1)
  1453. wr.jac(1,4+i)=0.0d0
  1454. 90 continue
  1455. c------------------------------
  1456. wr.jac(2,1)=-uold_r/rold_r
  1457. wr.jac(2,2)=1.0d0/rold_r
  1458. wr.jac(2,3)=0.0d0
  1459. wr.jac(2,4)=0.0d0
  1460. do 91 i=1,(nsp-1)
  1461. wr.jac(2,4+i)=0.0d0
  1462. 91 continue
  1463. c------------------------------
  1464. wr.jac(3,1)=-vold_r/rold_r
  1465. wr.jac(3,2)=0.0d0
  1466. wr.jac(3,3)=1.0d0/rold_r
  1467. wr.jac(3,4)=0.0d0
  1468. do 92 i=1,(nsp-1)
  1469. wr.jac(3,4+i)=0.0d0
  1470. 92 continue
  1471. c------------------------------
  1472. br1=0.0d0
  1473. do 860 i=1,(nsp-1)
  1474. br1=br1+dgdyr.vv(i)*yr.yet(i)
  1475. 860 continue
  1476. br1=br1*pold_r/(rold_r*gm1r)
  1477. wr.jac(4,1)=gm1r*(uold_r*uold_r+vold_r*vold_r)/2.0d0-br1
  1478. wr.jac(4,2)=-uold_r*gm1r
  1479. wr.jac(4,3)=-vold_r*gm1r
  1480. wr.jac(4,4)=gm1r
  1481. do 93 i=1,(nsp-1)
  1482. wr.jac(4,4+i)=dgdyr.vv(i)*pold_r/(rold_r*gm1r)
  1483. 93 continue
  1484. c------------------------------
  1485. do 94 i=1,(nsp-1)
  1486. do 95 j=1,4
  1487. wr.jac(4+i,j)=0.0d0
  1488. if(j.eq.1) wr.jac(4+i,j)=-yr.yet(i)/rold_r
  1489. 95 continue
  1490. c----------------------
  1491. do 950 j=5,(3+nsp)
  1492. wr.jac(4+i,j)=0.0d0
  1493. if(4+i.eq.j) wr.jac(4+i,j)=1.0d0/rold_r
  1494. 950 continue
  1495. 94 continue
  1496. c----------------------------------------------
  1497. c----------------------------------------------
  1498. do 1 i=1,(3+nsp)
  1499. do 2 j=1,(3+nsp)
  1500. jtl.jac(i,j)=0.0d0
  1501. jtr.jac(i,j)=0.0d0
  1502. do 3 k=1,(3+nsp)
  1503. jtl.jac(i,j)=jtl.jac(i,j)+jl.jac(i,k)*wl.jac(k,j)
  1504. jtr.jac(i,j)=jtr.jac(i,j)+jr.jac(i,k)*wr.jac(k,j)
  1505. 3 continue
  1506. 2 continue
  1507. 1 continue
  1508. c-------------------------------------------------------------------
  1509. SEGSUP DMLY_L, DMLY_R,
  1510. & dmry_l, dmry_r,
  1511. & dMpy_l, dMpy_r,
  1512. & dMmy_l, dMmy_r,
  1513. & dmiy_l, dmiy_r,
  1514. & d3my_l, d3my_r,
  1515. & d2my_l, d2my_r,
  1516. & dPpy_l, dPpy_r,
  1517. & dPmy_l, dPmy_r,
  1518. & dpiy_l, dpiy_r,
  1519. & dgdyl, dgdyr
  1520. C--------------------------------------------
  1521. SEGSUP f
  1522. C--------------------------------------------
  1523. c jll = jl
  1524. c jrr = jr
  1525. SEGDES JL
  1526. SEGDES JR
  1527. SEGSUP WL
  1528. SEGSUP WR
  1529. C--------------------------------------------
  1530. jll=jtl
  1531. SEGDES JTL
  1532. jrr=jtr
  1533. SEGDES JTR
  1534. SEGDES YL
  1535. SEGDES YR
  1536. SEGDES CP
  1537. SEGDES CV
  1538. SEGDES MLRECP, MLRECV
  1539. c-------------------
  1540. return
  1541. end
  1542.  
  1543.  
  1544.  
  1545.  
  1546.  
  1547.  
  1548.  
  1549.  
  1550.  
  1551.  
  1552.  
  1553.  
  1554.  
  1555.  
  1556.  
  1557.  
  1558.  
  1559.  
  1560.  
  1561.  
  1562.  
  1563.  
  1564.  
  1565.  
  1566.  
  1567.  
  1568.  
  1569.  
  1570.  
  1571.  
  1572.  
  1573.  
  1574.  
  1575.  
  1576.  

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