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
  536. dmlp_r=temp_l*damp_r
  537. dmrp_l=temp_r*damp_l
  538. dmrp_r=temp_r*damp_r
  539. c--------
  540. do 44 i=1,(nsp-1)
  541. dmly_l.vv(i)=temp_l*damg_l*dgdyl.vv(i)
  542. dmly_r.vv(i)=temp_l*damg_r*dgdyr.vv(i)
  543. dmry_l.vv(i)=temp_r*damg_l*dgdyl.vv(i)
  544. dmry_r.vv(i)=temp_r*damg_r*dgdyr.vv(i)
  545. 44 continue
  546. c-----------------------------
  547. sr_l=dmlr_l
  548. sr_r=dmlr_r
  549. su_l=dmlu_l
  550. su_r=dmlu_r
  551. sv_l=dmlv_l
  552. sv_r=dmlv_r
  553. sp_l=dmlp_l
  554. sp_r=dmlp_r
  555. do 440 i=1,(nsp-1)
  556. sy_l.vv(i) = dmly_l.vv(i)
  557. sy_r.vv(i) = dmly_r.vv(i)
  558. 440 continue
  559. c-----------------------------------------------------------------
  560. c Redefinition of the derivatives of the ml & mr
  561. c-----------------------------------------------------------------
  562. temp_l=(mlw-mrw)*mhalfr
  563. temp_r=-temp_l
  564. c--------
  565. dmlr_l=0.5d0*(top*dmlr_l+bot*dmrr_l)+canc*coef*temp_l*dmfr_l
  566. dmlu_l=0.5d0*(top*dmlu_l+bot*dmru_l)+canc*coef*temp_l*dmfu_l
  567. dmlv_l=0.5d0*(top*dmlv_l+bot*dmrv_l)+canc*coef*temp_l*dmfv_l
  568. dmlp_l=0.5d0*(top*dmlp_l+bot*dmrp_l)+canc*coef*temp_l*dmfp_l
  569. c--------
  570. dmlr_r=0.5d0*(top*dmlr_r+bot*dmrr_r)+canc*coef*temp_l*dmfr_r
  571. dmlu_r=0.5d0*(top*dmlu_r+bot*dmru_r)+canc*coef*temp_l*dmfu_r
  572. dmlv_r=0.5d0*(top*dmlv_r+bot*dmrv_r)+canc*coef*temp_l*dmfv_r
  573. dmlp_r=0.5d0*(top*dmlp_r+bot*dmrp_r)+canc*coef*temp_l*dmfp_r
  574. c--------
  575. do 441 i=1,(nsp-1)
  576. dmly_l.vv(i)=0.5d0*(top*dmly_l.vv(i)+bot*dmry_l.vv(i))+
  577. & canc*coef*temp_l*dmfg_l*dgdyl.vv(i)
  578. dmly_r.vv(i)=0.5d0*(top*dmly_r.vv(i)+bot*dmry_r.vv(i))+
  579. & canc*coef*temp_l*dmfg_r*dgdyr.vv(i)
  580. 441 continue
  581. c--------
  582. dmrr_l=0.5d0*(top*dmrr_l+bot*sr_l)+canc*coef*temp_r*dmfr_l
  583. dmru_l=0.5d0*(top*dmru_l+bot*su_l)+canc*coef*temp_r*dmfu_l
  584. dmrv_l=0.5d0*(top*dmrv_l+bot*sv_l)+canc*coef*temp_r*dmfv_l
  585. dmrp_l=0.5d0*(top*dmrp_l+bot*sp_l)+canc*coef*temp_r*dmfp_l
  586. c--------
  587. dmrr_r=0.5d0*(top*dmrr_r+bot*sr_r)+canc*coef*temp_r*dmfr_r
  588. dmru_r=0.5d0*(top*dmru_r+bot*su_r)+canc*coef*temp_r*dmfu_r
  589. dmrv_r=0.5d0*(top*dmrv_r+bot*sv_r)+canc*coef*temp_r*dmfv_r
  590. dmrp_r=0.5d0*(top*dmrp_r+bot*sp_r)+canc*coef*temp_r*dmfp_r
  591. do 442 i=1,(nsp-1)
  592. dmry_l.vv(i)=0.5d0*(top*dmry_l.vv(i)+bot*sy_l.vv(i))+
  593. & canc*coef*temp_r*dmfg_l*dgdyl.vv(i)
  594. dmry_r.vv(i)=0.5d0*(top*dmry_r.vv(i)+bot*sy_r.vv(i))+
  595. & canc*coef*temp_r*dmfg_r*dgdyr.vv(i)
  596. 442 continue
  597. c-----------------------------------------------------------
  598. c mmid is m_{1/2} (notation as in the paper, see (13),p.366)
  599. c-----------------------------------------------------------
  600. mmid=Mplus+Mmin
  601. c-----------------------------------------------------------
  602. c Computing the derivatives of M+ and M-
  603. c-----------------------------------------------------------
  604. if(ml .ge. 1.0d0) then
  605. dMpr_l=dmlr_l
  606. dMpu_l=dmlu_l
  607. dMpv_l=dmlv_l
  608. dMpp_l=dmlp_l
  609. do 45 i=1,(nsp-1)
  610. dMpy_l.vv(i)=dmly_l.vv(i)
  611. 45 continue
  612. c--------------------
  613. dMpr_r=dmlr_r
  614. dMpu_r=dmlu_r
  615. dMpv_r=dmlv_r
  616. dMpp_r=dmlp_r
  617. do 46 i=1,(nsp-1)
  618. dMpy_r.vv(i)=dmly_r.vv(i)
  619. 46 continue
  620. else
  621. if((ml .gt. -1.0d0) .and. (ml .lt. 1.0d0)) then
  622. temph=(ml+1.0d0)/2.0d0
  623. dMpr_l=(temph+4.0d0*beta*ml*(ml*ml-1.0d0))*dmlr_l
  624. dMpu_l=(temph+4.0d0*beta*ml*(ml*ml-1.0d0))*dmlu_l
  625. dMpv_l=(temph+4.0d0*beta*ml*(ml*ml-1.0d0))*dmlv_l
  626. dMpp_l=(temph+4.0d0*beta*ml*(ml*ml-1.0d0))*dmlp_l
  627. do 47 i=1,(nsp-1)
  628. dMpy_l.vv(i)=(temph+4.0d0*beta*ml*
  629. & (ml*ml-1.0d0))*dmly_l.vv(i)
  630. 47 continue
  631. c--------------------
  632. dMpr_r=(temph+4.0d0*beta*ml*(ml*ml-1.0d0))*dmlr_r
  633. dMpu_r=(temph+4.0d0*beta*ml*(ml*ml-1.0d0))*dmlu_r
  634. dMpv_r=(temph+4.0d0*beta*ml*(ml*ml-1.0d0))*dmlv_r
  635. dMpp_r=(temph+4.0d0*beta*ml*(ml*ml-1.0d0))*dmlp_r
  636. do 48 i=1,(nsp-1)
  637. dMpy_r.vv(i)=(temph+4.0d0*beta*ml*
  638. & (ml*ml-1.0d0))*dmly_r.vv(i)
  639. 48 continue
  640. else
  641. dMpr_l=0.0d0
  642. dMpu_l=0.0d0
  643. dMpv_l=0.0d0
  644. dMpp_l=0.0d0
  645. do 49 i=1,(nsp-1)
  646. dMpy_l.vv(i)=0.0d0
  647. 49 continue
  648. c---------------------
  649. dMpr_r=0.0d0
  650. dMpu_r=0.0d0
  651. dMpv_r=0.0d0
  652. dMpp_r=0.0d0
  653. do 50 i=1,(nsp-1)
  654. dMpy_r.vv(i)=0.0d0
  655. 50 continue
  656. endif
  657. endif
  658. c-----------------------------------------------------------
  659. c addition of low Mach number
  660. c-----------------------------------------------------------
  661. if(ml .ge. 0.0d0) then
  662. m1pr_l=dmlr_l
  663. m1pu_l=dmlu_l
  664. m1pv_l=dmlv_l
  665. m1pp_l=dmlp_l
  666. do 450 i=1,(nsp-1)
  667. m1py_l.vv(i)=dmly_l.vv(i)
  668. 450 continue
  669. c---------------------
  670. m1pr_r=dmlr_r
  671. m1pu_r=dmlu_r
  672. m1pv_r=dmlv_r
  673. m1pp_r=dmlp_r
  674. do 460 i=1,(nsp-1)
  675. m1py_r.vv(i)=dmly_r.vv(i)
  676. 460 continue
  677. else
  678. m1pr_l=0.0d0
  679. m1pu_l=0.0d0
  680. m1pv_l=0.0d0
  681. m1pp_l=0.0d0
  682. do 451 i=1,(nsp-1)
  683. m1py_l.vv(i)=0.0d0
  684. 451 continue
  685. c--------------------
  686. m1pr_r=0.0d0
  687. m1pu_r=0.0d0
  688. m1pv_r=0.0d0
  689. m1pp_r=0.0d0
  690. do 461 i=1,(nsp-1)
  691. m1py_r.vv(i)=0.0d0
  692. 461 continue
  693. endif
  694. c-----------------------------------------------------------
  695. if(mr .ge. 1.0d0) then
  696. dMmr_l=0.0d0
  697. dMmu_l=0.0d0
  698. dMmv_l=0.0d0
  699. dMmp_l=0.0d0
  700. do 51 i=1,(nsp-1)
  701. dMmy_l.vv(i)=0.0d0
  702. 51 continue
  703. c---------------------
  704. dMmr_r=0.0d0
  705. dMmu_r=0.0d0
  706. dMmv_r=0.0d0
  707. dMmp_r=0.0d0
  708. do 52 i=1,(nsp-1)
  709. dMmy_r.vv(i)=0.0d0
  710. 52 continue
  711. else
  712. if((mr .gt. -1.0d0) .and. (mr .lt. 1.0d0)) then
  713. temph=(-mr+1.0d0)/2.0d0
  714. dMmr_l=(temph-4.0d0*beta*mr*(mr*mr-1.0d0))*dmrr_l
  715. dMmu_l=(temph-4.0d0*beta*mr*(mr*mr-1.0d0))*dmru_l
  716. dMmv_l=(temph-4.0d0*beta*mr*(mr*mr-1.0d0))*dmrv_l
  717. dMmp_l=(temph-4.0d0*beta*mr*(mr*mr-1.0d0))*dmrp_l
  718. do 53 i=1,(nsp-1)
  719. dMmy_l.vv(i)=(temph-4.0d0*beta*mr*
  720. & (mr*mr-1.0d0))*dmry_l.vv(i)
  721. 53 continue
  722. c--------------------
  723. dMmr_r=(temph-4.0d0*beta*mr*(mr*mr-1.0d0))*dmrr_r
  724. dMmu_r=(temph-4.0d0*beta*mr*(mr*mr-1.0d0))*dmru_r
  725. dMmv_r=(temph-4.0d0*beta*mr*(mr*mr-1.0d0))*dmrv_r
  726. dMmp_r=(temph-4.0d0*beta*mr*(mr*mr-1.0d0))*dmrp_r
  727. do 54 i=1,(nsp-1)
  728. dMmy_r.vv(i)=(temph-4.0d0*beta*mr*
  729. & (mr*mr-1.0d0))*dmry_r.vv(i)
  730. 54 continue
  731. else
  732. dMmr_l=dmrr_l
  733. dMmu_l=dmru_l
  734. dMmv_l=dmrv_l
  735. dMmp_l=dmrp_l
  736. do 55 i=1,(nsp-1)
  737. dMmy_l.vv(i)=dmry_l.vv(i)
  738. 55 continue
  739. c--------------------
  740. dMmr_r=dmrr_r
  741. dMmu_r=dmru_r
  742. dMmv_r=dmrv_r
  743. dMmp_r=dmrp_r
  744. do 56 i=1,(nsp-1)
  745. dMmy_r.vv(i)=dmry_r.vv(i)
  746. 56 continue
  747. endif
  748. endif
  749. c-----------------------------------------------------------
  750. c addition of low Mach number
  751. c-----------------------------------------------------------
  752. if(mr .le. 0.0d0) then
  753. m1mr_l=dmrr_l
  754. m1mu_l=dmru_l
  755. m1mv_l=dmrv_l
  756. m1mp_l=dmrp_l
  757. do 550 i=1,(nsp-1)
  758. m1my_l.vv(i)=dmry_l.vv(i)
  759. 550 continue
  760. c---------------------
  761. m1mr_r=dmrr_r
  762. m1mu_r=dmru_r
  763. m1mv_r=dmrv_r
  764. m1mp_r=dmrp_r
  765. do 560 i=1,(nsp-1)
  766. m1my_r.vv(i)=dmry_r.vv(i)
  767. 560 continue
  768. else
  769. m1mr_l=0.0d0
  770. m1mu_l=0.0d0
  771. m1mv_l=0.0d0
  772. m1mp_l=0.0d0
  773. do 551 i=1,(nsp-1)
  774. m1my_l.vv(i)=0.0d0
  775. 551 continue
  776. c--------------------
  777. m1mr_r=0.0d0
  778. m1mu_r=0.0d0
  779. m1mv_r=0.0d0
  780. m1mp_r=0.0d0
  781. do 561 i=1,(nsp-1)
  782. m1my_r.vv(i)=0.0d0
  783. 561 continue
  784. endif
  785. c-----------------------------------------------------------------
  786. c computing the derivatives of m_{1/2} (notation as in the paper)
  787. c-----------------------------------------------------------------
  788. dmir_l=dMpr_l+dMmr_l
  789. dmir_r=dMpr_r+dMmr_r
  790. c-------------
  791. dmiu_l=dMpu_l+dMmu_l
  792. dmiu_r=dMpu_r+dMmu_r
  793. c-------------
  794. dmiv_l=dMpv_l+dMmv_l
  795. dmiv_r=dMpv_r+dMmv_r
  796. c-------------
  797. dmip_l=dMpp_l+dMmp_l
  798. dmip_r=dMpp_r+dMmp_r
  799. c-------------
  800. do 57 i=1,(nsp-1)
  801. dmiy_l.vv(i)=dMpy_l.vv(i)+dMmy_l.vv(i)
  802. dmiy_r.vv(i)=dMpy_r.vv(i)+dMmy_r.vv(i)
  803. 57 continue
  804. c----------------------------------------------------------------
  805. c computing the main convective variables and their derivatives
  806. c mpl_m is m^{+}_{1/2} (paper's notation) and
  807. c mmin_m is m^{-}_{1/2} (paper's notation), see (A2) on p.370.
  808. c----------------------------------------------------------------
  809. termp=(Mmin1-Mmin+Mplus-Mplus1)*(1.0d0/(mhalfr*mhalfr)-1.0d0)
  810. termp=termp*(pold_l-pold_r)/(pold_l/rold_l+pold_r/rold_r)
  811. c-------------------------------------------------------------
  812. c derivatives of the termp
  813. c-------------------------------------------------------------
  814. top=(Mmin1-Mmin+Mplus-Mplus1)
  815. bots=1.0d0/(pold_l/rold_l+pold_r/rold_r)
  816. bot=(pold_l-pold_r)*bots
  817. temph=1.0d0/(mhalfr*mhalfr)-1.0d0
  818. c---------------------------
  819. tmpr_l=(m1mr_l-dMmr_l+dMpr_l-m1pr_l)*bot*temph
  820. tmpr_l=tmpr_l-2.0d0*bot*top*dmfr_l/(mhalfr*mhalfr*mhalfr)
  821. tmpr_l=tmpr_l+temph*top*bot*bots*(pold_l/rold_l/rold_l)
  822. c---------------------------
  823. tmpu_l=(m1mu_l-dMmu_l+dMpu_l-m1pu_l)*bot*temph
  824. tmpu_l=tmpu_l-2.0d0*bot*top*dmfu_l/(mhalfr*mhalfr*mhalfr)
  825. c---------------------------
  826. tmpv_l=(m1mv_l-dMmv_l+dMpv_l-m1pv_l)*bot*temph
  827. tmpv_l=tmpv_l-2.0d0*bot*top*dmfv_l/(mhalfr*mhalfr*mhalfr)
  828. c---------------------------
  829. tmpp_l=(m1mp_l-dMmp_l+dMpp_l-m1pp_l)*bot*temph
  830. tmpp_l=tmpp_l-2.0d0*bot*top*dmfp_l/(mhalfr*mhalfr*mhalfr)
  831. tmpp_l=tmpp_l+temph*top*bots*(1.0d0-bot/rold_l)
  832. c---------------------------
  833. do 570 i=1,(nsp-1)
  834. tmpy_l.vv(i)=(m1my_l.vv(i)-dMmy_l.vv(i)+
  835. & dMpy_l.vv(i)-m1py_l.vv(i))*bot*temph
  836. tmpy_l.vv(i)=tmpy_l.vv(i)-2.0d0*bot*top*
  837. & dmfg_l*dgdyl.vv(i)/(mhalfr*mhalfr*mhalfr)
  838. 570 continue
  839. c------------rrrrrrrr-------
  840. c------------rrrrrrrr-------
  841. tmpr_r=(m1mr_r-dMmr_r+dMpr_r-m1pr_r)*bot*temph
  842. tmpr_r=tmpr_r-2.0d0*bot*top*dmfr_r/(mhalfr*mhalfr*mhalfr)
  843. tmpr_r=tmpr_r+temph*top*bot*bots*(pold_r/rold_r/rold_r)
  844. c---------------------------
  845. tmpu_r=(m1mu_r-dMmu_r+dMpu_r-m1pu_r)*bot*temph
  846. tmpu_r=tmpu_r-2.0d0*bot*top*dmfu_r/(mhalfr*mhalfr*mhalfr)
  847. c---------------------------
  848. tmpv_r=(m1mv_r-dMmv_r+dMpv_r-m1pv_r)*bot*temph
  849. tmpv_r=tmpv_r-2.0d0*bot*top*dmfv_r/(mhalfr*mhalfr*mhalfr)
  850. c---------------------------
  851. tmpp_r=(m1mp_r-dMmp_r+dMpp_r-m1pp_r)*bot*temph
  852. tmpp_r=tmpp_r-2.0d0*bot*top*dmfp_r/(mhalfr*mhalfr*mhalfr)
  853. tmpp_r=tmpp_r-temph*top*bots*(1.0d0+bot/rold_r)
  854. c---------------------------
  855. do 571 i=1,(nsp-1)
  856. tmpy_r.vv(i)=(m1my_r.vv(i)-dMmy_r.vv(i)+
  857. & dMpy_r.vv(i)-m1py_r.vv(i))*bot*temph
  858. tmpy_r.vv(i)=tmpy_r.vv(i)-2.0d0*bot*top*
  859. & dmfg_r*dgdyr.vv(i)/(mhalfr*mhalfr*mhalfr)
  860. 571 continue
  861. c-------------------------------------------------------------
  862. if(mmid .ge. 0.0d0) then
  863. mpl_m = mmid
  864. d2mr_l=dmir_l
  865. d2mu_l=dmiu_l
  866. d2mv_l=dmiv_l
  867. d2mp_l=dmip_l
  868. do 58 i=1,(nsp-1)
  869. d2my_l.vv(i)=dmiy_l.vv(i)
  870. 58 continue
  871. c------------
  872. d2mr_r=dmir_r
  873. d2mu_r=dmiu_r
  874. d2mv_r=dmiv_r
  875. d2mp_r=dmip_r
  876. do 59 i=1,(nsp-1)
  877. d2my_r.vv(i)=dmiy_r.vv(i)
  878. 59 continue
  879. c------------
  880. else
  881. mpl_m = 0.0d0
  882. d2mr_l=0.0d0
  883. d2mu_l=0.0d0
  884. d2mv_l=0.0d0
  885. d2mp_l=0.0d0
  886. do 60 i=1,(nsp-1)
  887. d2my_l.vv(i)=0.0d0
  888. 60 continue
  889. c------------
  890. d2mr_r=0.0d0
  891. d2mu_r=0.0d0
  892. d2mv_r=0.0d0
  893. d2mp_r=0.0d0
  894. do 61 i=1,(nsp-1)
  895. d2my_r.vv(i)=0.0d0
  896. 61 continue
  897. endif
  898. c---------------------------------------------
  899. cc derivatives for the term termm
  900. cc------------------------------------------------------------------
  901. if(mmid .le. 0.0d0) then
  902. mmin_m = mmid
  903. d3mr_l=dmir_l
  904. d3mu_l=dmiu_l
  905. d3mv_l=dmiv_l
  906. d3mp_l=dmip_l
  907. do 62 i=1,(nsp-1)
  908. d3my_l.vv(i)=dmiy_l.vv(i)
  909. 62 continue
  910. c------------
  911. d3mr_r=dmir_r
  912. d3mu_r=dmiu_r
  913. d3mv_r=dmiv_r
  914. d3mp_r=dmip_r
  915. do 63 i=1,(nsp-1)
  916. d3my_r.vv(i)=dmiy_r.vv(i)
  917. 63 continue
  918. c------------
  919. else
  920. mmin_m = 0.0d0
  921. d3mr_l=0.0d0
  922. d3mu_l=0.0d0
  923. d3mv_l=0.0d0
  924. d3mp_l=0.0d0
  925. do 64 i=1,(nsp-1)
  926. d3my_l.vv(i)=0.0d0
  927. 64 continue
  928. c------------
  929. d3mr_r=0.0d0
  930. d3mu_r=0.0d0
  931. d3mv_r=0.0d0
  932. d3mp_r=0.0d0
  933. do 65 i=1,(nsp-1)
  934. d3my_r.vv(i)=0.0d0
  935. 65 continue
  936. endif
  937. c---------------------------------------------------------------
  938. c Computing the calligraphic P+ and P- with their derivatives,
  939. c see (21a) & (21b) on p.368.
  940. c---------------------------------------------------------------
  941. if(ml .ge. 1.0d0) then
  942. Pplus = 1.0d0
  943. else
  944. if((ml .gt. -1.0d0) .and. (ml .lt. 1.0d0)) then
  945. Pplus=(ml+1.0d0)*(ml+1.0d0)*(2.0d0-ml)/4.0d0
  946. Pplus=Pplus+alpha*ml*(ml*ml-1.0d0)*(ml*ml-1.0d0)
  947. else
  948. Pplus = 0.0d0
  949. endif
  950. endif
  951. c---------------------------------------------------------------
  952. if(mr .ge. 1.0d0) then
  953. Pmin = 0.0d0
  954. else
  955. if((mr .gt. -1.0d0) .and. (mr .lt. 1.0d0)) then
  956. Pmin=(mr-1.0d0)*(mr-1.0d0)*(2.0d0+mr)/4.0d0
  957. Pmin=Pmin-alpha*mr*(mr*mr-1.0d0)*(mr*mr-1.0d0)
  958. else
  959. Pmin = 1.0d0
  960. endif
  961. endif
  962. c---------------------------------------------------------------
  963. brac_l=(ml+1.0d0)*(2.0d0-ml)/2.0d0-(ml+1.0d0)*(ml+1.0d0)/4.0d0
  964. brac_l=brac_l+alpha*(ml*ml-1.0d0)*(ml*ml-1.0d0)
  965. brac_l=brac_l+4.0d0*alpha*ml*ml*(ml*ml-1.0d0)
  966. c--------------
  967. brac_r=(mr-1.0d0)*(2.0d0+mr)/2.0d0+(mr-1.0d0)*(mr-1.0d0)/4.0d0
  968. brac_r=brac_r-alpha*(mr*mr-1.0d0)*(mr*mr-1.0d0)
  969. brac_r=brac_r-4.0d0*alpha*mr*mr*(mr*mr-1.0d0)
  970. c---------------------------------------------------------------
  971. if((ml .gt. -1.0d0) .and. (ml .lt. 1.0d0)) then
  972. dPpr_l=brac_l*dmlr_l
  973. dPpr_r=brac_l*dmlr_r
  974. c------------
  975. dPpu_l=brac_l*dmlu_l
  976. dPpu_r=brac_l*dmlu_r
  977. c------------
  978. dPpv_l=brac_l*dmlv_l
  979. dPpv_r=brac_l*dmlv_r
  980. c------------
  981. dPpp_l=brac_l*dmlp_l
  982. dPpp_r=brac_l*dmlp_r
  983. c------------
  984. do 66 i=1,(nsp-1)
  985. dPpy_l.vv(i)=brac_l*dmly_l.vv(i)
  986. dPpy_r.vv(i)=brac_l*dmly_r.vv(i)
  987. 66 continue
  988. c------------
  989. else
  990. dPpr_l=0.0d0
  991. dPpr_r=0.0d0
  992. c-----------
  993. dPpu_l=0.0d0
  994. dPpu_r=0.0d0
  995. c-----------
  996. dPpv_l=0.0d0
  997. dPpv_r=0.0d0
  998. c-----------
  999. dPpp_l=0.0d0
  1000. dPpp_r=0.0d0
  1001. c-----------
  1002. do 67 i=1,(nsp-1)
  1003. dPpy_l.vv(i)=0.0d0
  1004. dPpy_r.vv(i)=0.0d0
  1005. 67 continue
  1006. c-----------
  1007. endif
  1008. c---------------------------------------------------------------
  1009. if((mr .gt. -1.0d0) .and. (mr .lt. 1.0d0)) then
  1010. dPmr_l=brac_r*dmrr_l
  1011. dPmr_r=brac_r*dmrr_r
  1012. c------------
  1013. dPmu_l=brac_r*dmru_l
  1014. dPmu_r=brac_r*dmru_r
  1015. c------------
  1016. dPmv_l=brac_r*dmrv_l
  1017. dPmv_r=brac_r*dmrv_r
  1018. c------------
  1019. dPmp_l=brac_r*dmrp_l
  1020. dPmp_r=brac_r*dmrp_r
  1021. c------------
  1022. do 68 i=1,(nsp-1)
  1023. dPmy_l.vv(i)=brac_r*dmry_l.vv(i)
  1024. dPmy_r.vv(i)=brac_r*dmry_r.vv(i)
  1025. 68 continue
  1026. c------------
  1027. else
  1028. dPmr_l=0.0d0
  1029. dPmr_r=0.0d0
  1030. c-----------
  1031. dPmu_l=0.0d0
  1032. dPmu_r=0.0d0
  1033. c-----------
  1034. dPmv_l=0.0d0
  1035. dPmv_r=0.0d0
  1036. c-----------
  1037. dPmp_l=0.0d0
  1038. dPmp_r=0.0d0
  1039. c-----------
  1040. do 69 i=1,(nsp-1)
  1041. dPmy_l.vv(i)=0.0d0
  1042. dPmy_r.vv(i)=0.0d0
  1043. 69 continue
  1044. c-----------
  1045. endif
  1046. c-------------------------------------------------------------------
  1047. c computing pmid -- p_{1/2} and its derivatives, see (20b), p.367.
  1048. c-------------------------------------------------------------------
  1049. dpir_l=dPpr_l*pold_l+dPmr_l*pold_r
  1050. dpiu_l=dPpu_l*pold_l+dPmu_l*pold_r
  1051. dpiv_l=dPpv_l*pold_l+dPmv_l*pold_r
  1052. dpip_l=dPpp_l*pold_l+Pplus+dPmp_l*pold_r
  1053. do 70 i=1,(nsp-1)
  1054. dpiy_l.vv(i)=dPpy_l.vv(i)*pold_l+dPmy_l.vv(i)*pold_r
  1055. 70 continue
  1056. c----------------------------
  1057. dpir_r=dPpr_r*pold_l+dPmr_r*pold_r
  1058. dpiu_r=dPpu_r*pold_l+dPmu_r*pold_r
  1059. dpiv_r=dPpv_r*pold_l+dPmv_r*pold_r
  1060. dpip_r=dPpp_r*pold_l+Pmin+dPmp_r*pold_r
  1061. do 71 i=1,(nsp-1)
  1062. dpiy_r.vv(i)=dPpy_r.vv(i)*pold_l+dPmy_r.vv(i)*pold_r
  1063. 71 continue
  1064. c---------------------------------------------------------------------
  1065. c Computation of the mass flux (rho * u)_1/2
  1066. c---------------------------------------------------------------
  1067. rum=am*(mpl_m*rold_l+mmin_m*rold_r)+canc*am*termp
  1068. c-------------------------------------------------------
  1069. rumr_l=damr_l*(mpl_m*rold_l+mmin_m*rold_r)+
  1070. & am*(d2mr_l*rold_l+mpl_m+d3mr_l*rold_r)
  1071. rumr_l=rumr_l+canc*coef*(damr_l*termp+am*tmpr_l)
  1072. rumu_l=damu_l*(mpl_m*rold_l+mmin_m*rold_r)+
  1073. & am*(d2mu_l*rold_l+d3mu_l*rold_r)
  1074. rumu_l=rumu_l+canc*coef*(damu_l*termp+am*tmpu_l)
  1075. rumv_l=damv_l*(mpl_m*rold_l+mmin_m*rold_r)+
  1076. & am*(d2mv_l*rold_l+d3mv_l*rold_r)
  1077. rumv_l=rumv_l+canc*coef*(damv_l*termp+am*tmpv_l)
  1078. rump_l=damp_l*(mpl_m*rold_l+mmin_m*rold_r)+
  1079. & am*(d2mp_l*rold_l+d3mp_l*rold_r)
  1080. rump_l=rump_l+canc*coef*(damp_l*termp+am*tmpp_l)
  1081. c------------------
  1082. do 710 i=1,(nsp-1)
  1083. rumy_l.vv(i)=damg_l*dgdyl.vv(i)*(mpl_m*rold_l+
  1084. & mmin_m*rold_r)+am*(d2my_l.vv(i)*rold_l+d3my_l.vv(i)*rold_r)
  1085. rumy_l.vv(i)=rumy_l.vv(i)+canc*coef*
  1086. & (damg_l*dgdyl.vv(i)*termp+am*tmpy_l.vv(i))
  1087. 710 continue
  1088. c-------------------------------------------------
  1089. rumr_r=damr_r*(mpl_m*rold_l+mmin_m*rold_r)+
  1090. & am*(d2mr_r*rold_l+mmin_m+d3mr_r*rold_r)
  1091. rumr_r=rumr_r+canc*coef*(damr_r*termp+am*tmpr_r)
  1092. rumu_r=damu_r*(mpl_m*rold_l+mmin_m*rold_r)+
  1093. & am*(d2mu_r*rold_l+d3mu_r*rold_r)
  1094. rumu_r=rumu_r+canc*coef*(damu_r*termp+am*tmpu_r)
  1095. rumv_r=damv_r*(mpl_m*rold_l+mmin_m*rold_r)+
  1096. & am*(d2mv_r*rold_l+d3mv_r*rold_r)
  1097. rumv_r=rumv_r+canc*coef*(damv_r*termp+am*tmpv_r)
  1098. rump_r=damp_r*(mpl_m*rold_l+mmin_m*rold_r)+
  1099. & am*(d2mp_r*rold_l+d3mp_r*rold_r)
  1100. rump_r=rump_r+canc*coef*(damp_r*termp+am*tmpp_r)
  1101. do 711 i=1,(nsp-1)
  1102. rumy_r.vv(i)=damg_r*dgdyr.vv(i)*(mpl_m*rold_l+mmin_m*rold_r)+
  1103. & am*(d2my_r.vv(i)*rold_l+d3my_r.vv(i)*rold_r)
  1104. rumy_r.vv(i)=rumy_r.vv(i)+canc*coef*
  1105. & (damg_r*dgdyr.vv(i)*termp+am*tmpy_r.vv(i))
  1106. 711 continue
  1107. c---------------------------------------------------------------------
  1108. c computing JACOBIAN as a derivative of the numerical flux function
  1109. c with respect to the primitive variables.
  1110. c Notation: jl(i,j) --- is the derivative of the i-component of the
  1111. c flux function with respect to the j-component of the
  1112. c vector of primitive variables of the left state.
  1113. c jr(i,j) --- is the derivative of the i-component of the
  1114. c flux function with respect to the j-component of the
  1115. c vector of primitive variables of the right state.
  1116. c---------------------------------------------------------------------
  1117. jl.jac(1,1)=rumr_l
  1118. jl.jac(1,2)=rumu_l
  1119. jl.jac(1,3)=rumv_l
  1120. jl.jac(1,4)=rump_l
  1121. do 72 i=1,(nsp-1)
  1122. jl.jac(1,4+i)=rumy_l.vv(i)
  1123. 72 continue
  1124. c------------------------------------
  1125. jr.jac(1,1)=rumr_r
  1126. jr.jac(1,2)=rumu_r
  1127. jr.jac(1,3)=rumv_r
  1128. jr.jac(1,4)=rump_r
  1129. do 73 i=1,(nsp-1)
  1130. jr.jac(1,4+i)=rumy_r.vv(i)
  1131. 73 continue
  1132. c------------------------------------
  1133. c---------------------------------------------------------
  1134. if(rum .ge. 0.0d0) then
  1135. br3=rumr_l*un_l
  1136. br4=rumr_l*ut_l
  1137. else
  1138. br3=rumr_l*un_r
  1139. br4=rumr_l*ut_r
  1140. endif
  1141. jl.jac(2,1)=n_x*(br3+dpir_l)+br4*t_x
  1142. jl.jac(3,1)=n_y*(br3+dpir_l)+br4*t_y
  1143. c-------------------
  1144. if(rum .ge. 0.0d0) then
  1145. br3=rumu_l*un_l+rum*n_x
  1146. br4=rumu_l*ut_l+rum*t_x
  1147. else
  1148. br3=rumu_l*un_r
  1149. br4=rumu_l*ut_r
  1150. endif
  1151. jl.jac(2,2)=n_x*(br3+dpiu_l)+br4*t_x
  1152. jl.jac(3,2)=n_y*(br3+dpiu_l)+br4*t_y
  1153. c-------------------
  1154. if(rum .ge. 0.0d0) then
  1155. br3=rumv_l*un_l+rum*n_y
  1156. br4=rumv_l*ut_l+rum*t_y
  1157. else
  1158. br3=rumv_l*un_r
  1159. br4=rumv_l*ut_r
  1160. endif
  1161. jl.jac(2,3)=n_x*(br3+dpiv_l)+br4*t_x
  1162. jl.jac(3,3)=n_y*(br3+dpiv_l)+br4*t_y
  1163. c-------------------
  1164. if(rum .ge. 0.0d0) then
  1165. br3=rump_l*un_l
  1166. br4=rump_l*ut_l
  1167. else
  1168. br3=rump_l*un_r
  1169. br4=rump_l*ut_r
  1170. endif
  1171. jl.jac(2,4)=n_x*(br3+dpip_l)+br4*t_x
  1172. jl.jac(3,4)=n_y*(br3+dpip_l)+br4*t_y
  1173. c-------------------------------------------------------------
  1174. do 74 i=1,(nsp-1)
  1175. if(rum .ge. 0.0d0) then
  1176. br3=rumy_l.vv(i)*un_l
  1177. br4=rumy_l.vv(i)*ut_l
  1178. else
  1179. br3=rumy_l.vv(i)*un_r
  1180. br4=rumy_l.vv(i)*ut_r
  1181. endif
  1182. jl.jac(2,4+i)=n_x*(br3+dpiy_l.vv(i))+br4*t_x
  1183. jl.jac(3,4+i)=n_y*(br3+dpiy_l.vv(i))+br4*t_y
  1184. 74 continue
  1185. c-------------------------------------------------------------
  1186. if(rum .ge. 0.0d0) then
  1187. br3=rumr_r*un_l
  1188. br4=rumr_r*ut_l
  1189. else
  1190. br3=rumr_r*un_r
  1191. br4=rumr_r*ut_r
  1192. endif
  1193. jr.jac(2,1)=n_x*(br3+dpir_r)+br4*t_x
  1194. jr.jac(3,1)=n_y*(br3+dpir_r)+br4*t_y
  1195. c-------------------
  1196. if(rum .ge. 0.0d0) then
  1197. br3=rumu_r*un_l
  1198. br4=rumu_r*ut_l
  1199. else
  1200. br3=rumu_r*un_r+rum*n_x
  1201. br4=rumu_r*ut_r+rum*t_x
  1202. endif
  1203. jr.jac(2,2)=n_x*(br3+dpiu_r)+br4*t_x
  1204. jr.jac(3,2)=n_y*(br3+dpiu_r)+br4*t_y
  1205. c-------------------
  1206. if(rum .ge. 0.0d0) then
  1207. br3=rumv_r*un_l
  1208. br4=rumv_r*ut_l
  1209. else
  1210. br3=rumv_r*un_r+rum*n_y
  1211. br4=rumv_r*ut_r+rum*t_y
  1212. endif
  1213. jr.jac(2,3)=n_x*(br3+dpiv_r)+br4*t_x
  1214. jr.jac(3,3)=n_y*(br3+dpiv_r)+br4*t_y
  1215. c-------------------
  1216. if(rum .ge. 0.0d0) then
  1217. br3=rump_r*un_l
  1218. br4=rump_r*ut_l
  1219. else
  1220. br3=rump_r*un_r
  1221. br4=rump_r*ut_r
  1222. endif
  1223. jr.jac(2,4)=n_x*(br3+dpip_r)+br4*t_x
  1224. jr.jac(3,4)=n_y*(br3+dpip_r)+br4*t_y
  1225. c--------------------
  1226. do 75 i=1,(nsp-1)
  1227. if(rum .ge. 0.0d0) then
  1228. br3=rumy_r.vv(i)*un_l
  1229. br4=rumy_r.vv(i)*ut_l
  1230. else
  1231. br3=rumy_r.vv(i)*un_r
  1232. br4=rumy_r.vv(i)*ut_r
  1233. endif
  1234. jr.jac(2,4+i)=n_x*(br3+dpiy_r.vv(i))+br4*t_x
  1235. jr.jac(3,4+i)=n_y*(br3+dpiy_r.vv(i))+br4*t_y
  1236. 75 continue
  1237. c-------------------------------------------------------------
  1238. c ------ f44444444444444444444444444444 ---------
  1239. c-------------------------------------------------------------
  1240. hr_l=(uold_l*uold_l+vold_l*vold_l)/2.0d0+
  1241. & gal*pold_l/gm1l/rold_l
  1242. hr_r=(uold_r*uold_r+vold_r*vold_r)/2.0d0+
  1243. & gar*pold_r/gm1r/rold_r
  1244. c-------------------------------------------------------------
  1245. if(rum .ge. 0.0d0) then
  1246. br1=rumr_l*hr_l-rum*gal*pold_l/gm1l/rold_l/rold_l
  1247. else
  1248. br1=rumr_l*hr_r
  1249. endif
  1250. jl.jac(4,1)=br1
  1251. c---------------------
  1252. if(rum .ge. 0.0d0) then
  1253. br1=rumu_l*hr_l+rum*uold_l
  1254. else
  1255. br1=rumu_l*hr_r
  1256. endif
  1257. jl.jac(4,2)=br1
  1258. c---------------------
  1259. if(rum .ge. 0.0d0) then
  1260. br1=rumv_l*hr_l+rum*vold_l
  1261. else
  1262. br1=rumv_l*hr_r
  1263. endif
  1264. jl.jac(4,3)=br1
  1265. c---------------------
  1266. if(rum .ge. 0.0d0) then
  1267. br1=rump_l*hr_l+rum*gal/gm1l/rold_l
  1268. else
  1269. br1=rump_l*hr_r
  1270. endif
  1271. jl.jac(4,4)=br1
  1272. c---------------------
  1273. do 76 i=1,(nsp-1)
  1274. if(rum .ge. 0.0d0) then
  1275. br1=rumy_l.vv(i)*hr_l-
  1276. & rum*(pold_l*dgdyl.vv(i))/(rold_l*gm1l*gm1l)
  1277. else
  1278. br1=rumy_l.vv(i)*hr_r
  1279. endif
  1280. jl.jac(4,4+i)=br1
  1281. 76 continue
  1282. c----------------------------------------------------------
  1283. if(rum .ge. 0.0d0) then
  1284. br1=rumr_r*hr_l
  1285. else
  1286. br1=rumr_r*hr_r-rum*gar*pold_r/gm1r/rold_r/rold_r
  1287. endif
  1288. jr.jac(4,1)=br1
  1289. c---------------------
  1290. if(rum .ge. 0.0d0) then
  1291. br1=rumu_r*hr_l
  1292. else
  1293. br1=rumu_r*hr_r+rum*uold_r
  1294. endif
  1295. jr.jac(4,2)=br1
  1296. c---------------------
  1297. if(rum .ge. 0.0d0) then
  1298. br1=rumv_r*hr_l
  1299. else
  1300. br1=rumv_r*hr_r+rum*vold_r
  1301. endif
  1302. jr.jac(4,3)=br1
  1303. c---------------------
  1304. if(rum .ge. 0.0d0) then
  1305. br1=rump_r*hr_l
  1306. else
  1307. br1=rump_r*hr_r+rum*gar/gm1r/rold_r
  1308. endif
  1309. jr.jac(4,4)=br1
  1310. c---------------------
  1311. do 77 i=1,(nsp-1)
  1312. if(rum .ge. 0.0d0) then
  1313. br1=rumy_r.vv(i)*hr_l
  1314. else
  1315. br1=rumy_r.vv(i)*hr_r-
  1316. & rum*(pold_r*dgdyr.vv(i))/(rold_r*gm1r*gm1r)
  1317. endif
  1318. jr.jac(4,4+i)=br1
  1319. 77 continue
  1320. c-------------------------------------------------------------
  1321. c-------------------------------------------------------------
  1322. do 78 i=1,(nsp-1)
  1323. if(rum .ge. 0.0d0) then
  1324. br1=rumr_l*yl.yet(i)
  1325. br2=rumr_r*yl.yet(i)
  1326. else
  1327. br1=rumr_l*yr.yet(i)
  1328. br2=rumr_r*yr.yet(i)
  1329. endif
  1330. jl.jac(4+i,1)=br1
  1331. jr.jac(4+i,1)=br2
  1332. c-------------------------
  1333. if(rum .ge. 0.0d0) then
  1334. br1=rumu_l*yl.yet(i)
  1335. br2=rumu_r*yl.yet(i)
  1336. else
  1337. br1=rumu_l*yr.yet(i)
  1338. br2=rumu_r*yr.yet(i)
  1339. endif
  1340. jl.jac(4+i,2)=br1
  1341. jr.jac(4+i,2)=br2
  1342. c-------------------------
  1343. if(rum .ge. 0.0d0) then
  1344. br1=rumv_l*yl.yet(i)
  1345. br2=rumv_r*yl.yet(i)
  1346. else
  1347. br1=rumv_l*yr.yet(i)
  1348. br2=rumv_r*yr.yet(i)
  1349. endif
  1350. jl.jac(4+i,3)=br1
  1351. jr.jac(4+i,3)=br2
  1352. c-------------------------
  1353. if(rum .ge. 0.0d0) then
  1354. br1=rump_l*yl.yet(i)
  1355. br2=rump_r*yl.yet(i)
  1356. else
  1357. br1=rump_l*yr.yet(i)
  1358. br2=rump_r*yr.yet(i)
  1359. endif
  1360. jl.jac(4+i,4)=br1
  1361. jr.jac(4+i,4)=br2
  1362. c--------------------------
  1363. do 780 j=1,(nsp-1)
  1364. if(rum .ge. 0.0d0) then
  1365. if(i .eq. j) then
  1366. br1=rumy_l.vv(j)*yl.yet(i)+rum
  1367. br2=rumy_r.vv(j)*yl.yet(i)
  1368. else
  1369. br1=rumy_l.vv(j)*yl.yet(i)
  1370. br2=rumy_r.vv(j)*yl.yet(i)
  1371. endif
  1372. c--------------
  1373. else
  1374. if(i .eq. j) then
  1375. br1=rumy_l.vv(j)*yr.yet(i)
  1376. br2=rumy_r.vv(j)*yr.yet(i)+rum
  1377. else
  1378. br1=rumy_l.vv(j)*yr.yet(i)
  1379. br2=rumy_r.vv(j)*yr.yet(i)
  1380. endif
  1381. endif
  1382. jl.jac(4+i,4+j)=br1
  1383. jr.jac(4+i,4+j)=br2
  1384. 780 continue
  1385. 78 continue
  1386. c-------------------------------------------------------------
  1387. c matrix wl(i,j) represents the derivative of the i-component
  1388. c of the vector of primitive variables of the left state with
  1389. c respect to the j-component of the vector of the conservative
  1390. c variables of the left state.
  1391. c
  1392. c Here: (rho, ux, uy, p, Y_1,...,Y_(nsp-1)) -
  1393. c vector of primitive variables;
  1394. c (rho, rho ux, rho uy, rho e, rho Y_1,..., rho Y_(nsp-1)) -
  1395. c vector of conservative variables.
  1396. c-------------------------------------------------------------
  1397. wl.jac(1,1)=1.0d0
  1398. wl.jac(1,2)=0.0d0
  1399. wl.jac(1,3)=0.0d0
  1400. wl.jac(1,4)=0.0d0
  1401. do 83 i=1,(nsp-1)
  1402. wl.jac(1,4+i)=0.0d0
  1403. 83 continue
  1404. c------------------------------
  1405. wl.jac(2,1)=-uold_l/rold_l
  1406. wl.jac(2,2)=1.0d0/rold_l
  1407. wl.jac(2,3)=0.0d0
  1408. wl.jac(2,4)=0.0d0
  1409. do 84 i=1,(nsp-1)
  1410. wl.jac(2,4+i)=0.0d0
  1411. 84 continue
  1412. c------------------------------
  1413. wl.jac(3,1)=-vold_l/rold_l
  1414. wl.jac(3,2)=0.0d0
  1415. wl.jac(3,3)=1.0d0/rold_l
  1416. wl.jac(3,4)=0.0d0
  1417. do 85 i=1,(nsp-1)
  1418. wl.jac(3,4+i)=0.0d0
  1419. 85 continue
  1420. c------------------------------
  1421. br1=0.0d0
  1422. do 86 i=1,(nsp-1)
  1423. br1=br1+dgdyl.vv(i)*yl.yet(i)
  1424. 86 continue
  1425. br1=br1*pold_l/(rold_l*gm1l)
  1426. wl.jac(4,1)=gm1l*(uold_l*uold_l+vold_l*vold_l)/2.0d0-br1
  1427. wl.jac(4,2)=-uold_l*gm1l
  1428. wl.jac(4,3)=-vold_l*gm1l
  1429. wl.jac(4,4)=gm1l
  1430. do 87 i=1,(nsp-1)
  1431. wl.jac(4,4+i)=dgdyl.vv(i)*pold_l/(rold_l*gm1l)
  1432. 87 continue
  1433. c------------------------------
  1434. do 88 i=1,(nsp-1)
  1435. do 89 j=1,4
  1436. wl.jac(4+i,j)=0.0d0
  1437. if(j.eq.1) wl.jac(4+i,j)=-yl.yet(i)/rold_l
  1438. 89 continue
  1439. c------------
  1440. do 890 j=5,(4+nsp-1)
  1441. wl.jac(4+i,j)=0.0d0
  1442. if(4+i.eq.j) then
  1443. wl.jac(4+i,j)=1.0d0/rold_l
  1444. endif
  1445. 890 continue
  1446. 88 continue
  1447. c------------------------------
  1448. c------------------------------
  1449. wr.jac(1,1)=1.0d0
  1450. wr.jac(1,2)=0.0d0
  1451. wr.jac(1,3)=0.0d0
  1452. wr.jac(1,4)=0.0d0
  1453. do 90 i=1,(nsp-1)
  1454. wr.jac(1,4+i)=0.0d0
  1455. 90 continue
  1456. c------------------------------
  1457. wr.jac(2,1)=-uold_r/rold_r
  1458. wr.jac(2,2)=1.0d0/rold_r
  1459. wr.jac(2,3)=0.0d0
  1460. wr.jac(2,4)=0.0d0
  1461. do 91 i=1,(nsp-1)
  1462. wr.jac(2,4+i)=0.0d0
  1463. 91 continue
  1464. c------------------------------
  1465. wr.jac(3,1)=-vold_r/rold_r
  1466. wr.jac(3,2)=0.0d0
  1467. wr.jac(3,3)=1.0d0/rold_r
  1468. wr.jac(3,4)=0.0d0
  1469. do 92 i=1,(nsp-1)
  1470. wr.jac(3,4+i)=0.0d0
  1471. 92 continue
  1472. c------------------------------
  1473. br1=0.0d0
  1474. do 860 i=1,(nsp-1)
  1475. br1=br1+dgdyr.vv(i)*yr.yet(i)
  1476. 860 continue
  1477. br1=br1*pold_r/(rold_r*gm1r)
  1478. wr.jac(4,1)=gm1r*(uold_r*uold_r+vold_r*vold_r)/2.0d0-br1
  1479. wr.jac(4,2)=-uold_r*gm1r
  1480. wr.jac(4,3)=-vold_r*gm1r
  1481. wr.jac(4,4)=gm1r
  1482. do 93 i=1,(nsp-1)
  1483. wr.jac(4,4+i)=dgdyr.vv(i)*pold_r/(rold_r*gm1r)
  1484. 93 continue
  1485. c------------------------------
  1486. do 94 i=1,(nsp-1)
  1487. do 95 j=1,4
  1488. wr.jac(4+i,j)=0.0d0
  1489. if(j.eq.1) wr.jac(4+i,j)=-yr.yet(i)/rold_r
  1490. 95 continue
  1491. c----------------------
  1492. do 950 j=5,(3+nsp)
  1493. wr.jac(4+i,j)=0.0d0
  1494. if(4+i.eq.j) wr.jac(4+i,j)=1.0d0/rold_r
  1495. 950 continue
  1496. 94 continue
  1497. c----------------------------------------------
  1498. c----------------------------------------------
  1499. do 1 i=1,(3+nsp)
  1500. do 2 j=1,(3+nsp)
  1501. jtl.jac(i,j)=0.0d0
  1502. jtr.jac(i,j)=0.0d0
  1503. do 3 k=1,(3+nsp)
  1504. jtl.jac(i,j)=jtl.jac(i,j)+jl.jac(i,k)*wl.jac(k,j)
  1505. jtr.jac(i,j)=jtr.jac(i,j)+jr.jac(i,k)*wr.jac(k,j)
  1506. 3 continue
  1507. 2 continue
  1508. 1 continue
  1509. c-------------------------------------------------------------------
  1510. SEGSUP DMLY_L, DMLY_R,
  1511. & dmry_l, dmry_r,
  1512. & dMpy_l, dMpy_r,
  1513. & dMmy_l, dMmy_r,
  1514. & dmiy_l, dmiy_r,
  1515. & d3my_l, d3my_r,
  1516. & d2my_l, d2my_r,
  1517. & dPpy_l, dPpy_r,
  1518. & dPmy_l, dPmy_r,
  1519. & dpiy_l, dpiy_r,
  1520. & dgdyl, dgdyr
  1521. C--------------------------------------------
  1522. SEGSUP f
  1523. C--------------------------------------------
  1524. c jll = jl
  1525. c jrr = jr
  1526. SEGDES JL
  1527. SEGDES JR
  1528. SEGSUP WL
  1529. SEGSUP WR
  1530. C--------------------------------------------
  1531. jll=jtl
  1532. SEGDES JTL
  1533. jrr=jtr
  1534. SEGDES JTR
  1535. SEGDES YL
  1536. SEGDES YR
  1537. SEGDES CP
  1538. SEGDES CV
  1539. SEGDES MLRECP, MLRECV
  1540. c-------------------
  1541. return
  1542. end
  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.  
  1577.  

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