Télécharger ja3dbm.eso

Retour à la liste

Numérotation des lignes :

  1. C JA3DBM SOURCE CHAT 05/01/13 00:47:49 5004
  2. SUBROUTINE JA3DBM(JTL,JTR,WVEC_L,WVEC_R,NVECT,TVEC1,TVEC2,
  3. & ga,v_inf)
  4. C************************************************************************
  5. C
  6. C PROJET : CASTEM 2000
  7. C
  8. C NOM : JA3DBM (Bas-Mach jacobian for 3D)
  9. C
  10. C DESCRIPTION : Voir KONJA2
  11. C
  12. C LANGAGE : FORTRAN 77 + ESOPE 2000 (avec estensions CISI)
  13. C
  14. C AUTEUR : S. KUDRIAKOV, A. BECCANTINI DM2S/SFME/LTMF
  15. C
  16. C************************************************************************
  17. C
  18. c----------------------------------------------------------------------
  19. c GENERAL DESCRIPTION:
  20. c This subroutine provides the jacobians, which are
  21. c the derivatives of the numerical flux function,
  22. c with respect to the right and left conservative variables
  23. c The low-mach number corrections are made.
  24. c
  25. c EQUATIONS: 3D Euler equations of gas dynamics
  26. c
  27. c
  28. c REFERENCE: 1) JCP, 129, 364-382 (1996)
  29. c " A Sequel to AUSM: AUSM+ ";
  30. c M.S.Liou
  31. c 2) AIAA Journal, Sept. 1998
  32. c "Low-Diffusion Flux-Splitting Methods for Flows at All Speeds"
  33. c J.R.Edwards and M.S.Liou
  34. c----------------------------------------------------------------------
  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 == (gamg,rhog,pg,ung,utg,uvg) --
  45. c vector of the primitive variables
  46. c at the left cell;
  47. c
  48. c wvec_r == (gamd,rhod,pd,und,utd,uvd) --
  49. c vector of the primitive variables
  50. c at the right cell;
  51. c
  52. c nvect, tvec1, tvec2 --- normal and tangent vectors
  53. c to the cell interface
  54. c
  55. c ga ---- ratio of the specific heats of the gas
  56. c
  57. c v_inf -- parameter for choosing the reference velocity
  58. c when the magnitude of the physical velocity
  59. c is close to zero
  60.  
  61. c----------------------------------------------------------------------
  62. c
  63. c OUTPUT:
  64. c
  65. c jtl -- jakobian matrix 5 by 5 - derivatives of the numerical
  66. c flux function with respect to the conservative variables
  67. c from the left cell;
  68. c
  69. c jtr -- jakobian matrix 5 by 5 - derivatives of the numerical
  70. c flux function with respect to the conservative variables
  71. c from the right cell.
  72. c----------------------------------------------------------------------
  73. IMPLICIT INTEGER(I-N)
  74. real*8 wvec_l(5),wvec_r(5)
  75. real*8 nvect(3),tvec1(3),tvec2(3)
  76. real*8 jl(5,5),jr(5,5)
  77. real*8 wl(5,5),wr(5,5)
  78. real*8 jtl(5,5),jtr(5,5)
  79. real*8 alpha,beta,upr_l,upr_r,v_inf
  80. real*8 ga, gm1,temph
  81. real*8 n_x,n_y,n_z,t1_x,t1_y,t1_z,t2_x,t2_y,t2_z
  82. real*8 un_l,un_r,ut1_l,ut1_r,ut2_l,ut2_r
  83. real*8 ml,mr,Mplus,Mmin,mmid
  84. real*8 mpl_m, mmin_m,am
  85. real*8 rold_l,uold_l,vold_l,wold_l,pold_l,eold_l
  86. real*8 rold_r,uold_r,vold_r,wold_r,pold_r,eold_r
  87. real*8 Pplus,Pmin
  88. real*8 hr_l,hr_r
  89. real*8 br1,br2,br3,temp_l,temp_r,brac_l,brac_r
  90. real*8 aleft, arigh
  91. real*8 damr_l,damr_r,damu_l,damu_r
  92. real*8 damv_l,damv_r,damp_l,damp_r
  93. real*8 damw_l,damw_r
  94. real*8 dmlr_l,dmlr_r,dmlu_l,dmlu_r
  95. real*8 dmlv_l,dmlv_r,dmlp_l,dmlp_r
  96. real*8 dmlw_l,dmlw_r
  97. real*8 dmrr_l,dmrr_r,dmru_l,dmru_r
  98. real*8 dmrv_l,dmrv_r,dmrp_l,dmrp_r
  99. real*8 dmrw_l,dmrw_r
  100. real*8 dMpr_l,dMpr_r,dMpu_l,dMpu_r
  101. real*8 dMpv_l,dMpv_r,dMpp_l,dMpp_r
  102. real*8 dMpw_l,dMpw_r
  103. real*8 dMmr_l,dMmr_r,dMmu_l,dMmu_r
  104. real*8 dMmv_l,dMmv_r,dMmp_l,dMmp_r
  105. real*8 dMmw_l,dMmw_r
  106. real*8 dmir_l,dmir_r,dmiu_l,dmiu_r
  107. real*8 dmiv_l,dmiv_r,dmip_l,dmip_r
  108. real*8 dmiw_l,dmiw_r
  109. real*8 d3mr_l,d3mr_r,d3mu_l,d3mu_r
  110. real*8 d3mv_l,d3mv_r,d3mp_l,d3mp_r
  111. real*8 d3mw_l,d3mw_r
  112. real*8 d2mr_l,d2mr_r,d2mu_l,d2mu_r
  113. real*8 d2mv_l,d2mv_r,d2mp_l,d2mp_r
  114. real*8 d2mw_l,d2mw_r
  115. real*8 dPpr_l,dPpr_r,dPpu_l,dPpu_r
  116. real*8 dPpv_l,dPpv_r,dPpp_l,dPpp_r
  117. real*8 dPpw_l,dPpw_r
  118. real*8 dPmr_l,dPmr_r,dPmu_l,dPmu_r
  119. real*8 dPmv_l,dPmv_r,dPmp_l,dPmp_r
  120. real*8 dPmw_l,dPmw_r
  121. real*8 dpir_l,dpir_r,dpiu_l,dpiu_r
  122. real*8 dpiv_l,dpiv_r,dpip_l,dpip_r
  123. real*8 dpiw_l,dpiw_r
  124. real*8 epsil,qq,amw,Mmin1,Mplus1
  125. real*8 fmid,mlw,mrw,termp
  126. real*8 ur_r,ur_l,urm,mhalf,mhalfr
  127. real*8 durr_l,durr_r,duru_l,duru_r
  128. real*8 durv_l,durv_r,durp_l,durp_r
  129. real*8 durw_l,durw_r
  130. real*8 dmhr_l,dmhr_r,dmhu_l,dmhu_r
  131. real*8 dmhv_l,dmhv_r,dmhp_l,dmhp_r
  132. real*8 dmhw_l,dmhw_r
  133. real*8 dmfr_l,dmfr_r,dmfu_l,dmfu_r
  134. real*8 dmfv_l,dmfv_r,dmfp_l,dmfp_r
  135. real*8 dmfw_l,dmfw_r
  136. real*8 dfm_mf,dfm_mh
  137. real*8 dfmr_l,dfmr_r,dfmu_l,dfmu_r
  138. real*8 dfmv_l,dfmv_r,dfmp_l,dfmp_r
  139. real*8 dfmw_l,dfmw_r
  140. real*8 m1mr_l,m1mr_r,m1mu_l,m1mu_r
  141. real*8 m1mv_l,m1mv_r,m1mp_l,m1mp_r
  142. real*8 m1mw_l,m1mw_r
  143. real*8 m1pr_l,m1pr_r,m1pu_l,m1pu_r
  144. real*8 m1pv_l,m1pv_r,m1pp_l,m1pp_r
  145. real*8 m1pw_l,m1pw_r
  146. real*8 tmpr_l,tmpr_r,tmpu_l,tmpu_r
  147. real*8 tmpv_l,tmpv_r,tmpp_l,tmpp_r
  148. real*8 tmpw_l,tmpw_r,rum
  149. real*8 rumr_l,rumu_l,rumv_l,rumw_l,rump_l
  150. real*8 rumr_r,rumu_r,rumv_r,rumw_r,rump_r
  151. real*8 sw_l,sw_r,top,bot,bots,coef,canc
  152. real*8 sr_l,sr_r,su_l,su_r,sv_l,sv_r,sp_l,sp_r
  153. integer i,j,k
  154. parameter(alpha = 0.1875D0, beta = 0.125D0)
  155. parameter(epsil = 1.0d0)
  156. canc=1.0d0
  157. upr_l=0.0d0
  158. upr_r=0.0d0
  159. c-----------------------
  160. gm1=ga-1.0d0
  161. c-----------------------
  162. n_x=nvect(1)
  163. n_y=nvect(2)
  164. n_z=nvect(3)
  165. c-------------------
  166. t1_x=tvec1(1)
  167. t1_y=tvec1(2)
  168. t1_z=tvec1(3)
  169. c-------------------
  170. t2_x=tvec2(1)
  171. t2_y=tvec2(2)
  172. t2_z=tvec2(3)
  173. c-------------------------------------
  174. rold_l=wvec_l(1)
  175. uold_l=wvec_l(2)
  176. vold_l=wvec_l(3)
  177. wold_l=wvec_l(4)
  178. pold_l=wvec_l(5)
  179. c-----------------------
  180. rold_r=wvec_r(1)
  181. uold_r=wvec_r(2)
  182. vold_r=wvec_r(3)
  183. wold_r=wvec_r(4)
  184. pold_r=wvec_r(5)
  185. c-----------------------
  186. eold_l=(uold_l*uold_l+vold_l*vold_l+wold_l*wold_l)/2.0d0
  187. eold_l=eold_l+pold_l/(gm1*rold_l)
  188. eold_r=(uold_r*uold_r+vold_r*vold_r+wold_r*wold_r)/2.0d0
  189. eold_r=eold_r+pold_r/(gm1*rold_r)
  190. c------------------------------------------------------------------
  191. c Computing reference velocity and its derivatives
  192. c see Eq.(2) of the Ref.(2).
  193. c------------------------------------------------------------------
  194. aleft=sqrt(ga*pold_l/rold_l)
  195. arigh=sqrt(ga*pold_r/rold_r)
  196. qq=sqrt(uold_l*uold_l+vold_l*vold_l+wold_l*wold_l)
  197. if(qq .lt. (epsil*v_inf)) then
  198. ur_l = epsil*v_inf
  199. durr_l=0.0d0
  200. duru_l=0.0d0
  201. durv_l=0.0d0
  202. durw_l=0.0d0
  203. durp_l=0.0d0
  204. else
  205. ur_l=qq
  206. durr_l=0.0d0
  207. duru_l=uold_l/qq
  208. durv_l=vold_l/qq
  209. durw_l=wold_l/qq
  210. durp_l=0.0d0
  211. endif
  212. c-------------------------------------------------------------------
  213. if(ur_l .ge. aleft) then
  214. ur_l=aleft
  215. durr_l=-aleft/(2.0d0*rold_l)
  216. duru_l=0.0d0
  217. durv_l=0.0d0
  218. durw_l=0.0d0
  219. durp_l=ga/(2.0d0*aleft*rold_l)
  220. endif
  221. c-------------------------------------------------------------------
  222. if(ur_l .lt. upr_l) then
  223. ur_l=upr_l
  224. durr_l=0.0d0
  225. duru_l=0.0d0
  226. durv_l=0.0d0
  227. durw_l=0.0d0
  228. durp_l=0.0d0
  229. endif
  230. c-------------------------------------------------------------------
  231. qq=sqrt(uold_r*uold_r+vold_r*vold_r+wold_r*wold_r)
  232. if(qq .lt. (epsil*v_inf)) then
  233. ur_r = epsil*v_inf
  234. durr_r=0.0d0
  235. duru_r=0.0d0
  236. durv_r=0.0d0
  237. durw_r=0.0d0
  238. durp_r=0.0d0
  239. else
  240. ur_r=qq
  241. durr_r=0.0d0
  242. duru_r=uold_r/qq
  243. durv_r=vold_r/qq
  244. durw_r=wold_r/qq
  245. durp_r=0.0d0
  246. endif
  247. c------------------------------------------------------------------
  248. if(ur_r .ge. arigh) then
  249. ur_r=arigh
  250. durr_r=-arigh/(2.0d0*rold_r)
  251. duru_r=0.0d0
  252. durv_r=0.0d0
  253. durw_r=0.0d0
  254. durp_r=ga/(2.0d0*arigh*rold_r)
  255. endif
  256. c------------------------------------------------------------------
  257. if(ur_r .lt. upr_r) then
  258. ur_r=upr_r
  259. durr_r=0.0d0
  260. duru_r=0.0d0
  261. durv_r=0.0d0
  262. durw_r=0.0d0
  263. durp_r=0.0d0
  264. endif
  265. c------------------------------------------------------------------
  266. c Reference velocity at the interface is taken as an average
  267. c of the reference velocities of the neighbouring cells
  268. c-------------------------------------------------------------------
  269. urm=0.5d0*(ur_l+ur_r)
  270. durr_l=0.5d0*durr_l
  271. duru_l=0.5d0*duru_l
  272. durv_l=0.5d0*durv_l
  273. durw_l=0.5d0*durw_l
  274. durp_l=0.5d0*durp_l
  275. c-------------------------
  276. durr_r=0.5d0*durr_r
  277. duru_r=0.5d0*duru_r
  278. durv_r=0.5d0*durv_r
  279. durw_r=0.5d0*durw_r
  280. durp_r=0.5d0*durp_r
  281. c-------------------------------------------------------------------
  282. c Computation of the speed of sound and its derivatives;
  283. c numerical speed of sound at the interface is taken as an average
  284. c of the speeds of sounds of the neighbouring cells
  285. c-------------------------------------------------------------------
  286. am=0.5d0*(aleft+arigh)
  287. c-------------------------------------------------------------------
  288. if(abs(urm/am-1.0d0).le.0.000001d0) then
  289. coef=0.0d0
  290. else
  291. coef=1.0d0
  292. endif
  293. c-------------------------------------------------------------------
  294. c Computation of the speed of sound and its derivatives;
  295. c numerical speed of sound at the interface is taken as an average
  296. c of the speeds of sounds of the neighbouring cells
  297. c--------------------------------------------------------------------
  298. damr_r=-arigh/(4.0d0*rold_r)
  299. damu_r=0.0d0
  300. damv_r=0.0d0
  301. damw_r=0.0d0
  302. damp_r=ga/(4.0d0*arigh*rold_r)
  303. c-----------------------
  304. damr_l=-aleft/(4.0d0*rold_l)
  305. damu_l=0.0d0
  306. damv_l=0.0d0
  307. damw_l=0.0d0
  308. damp_l=ga/(4.0d0*aleft*rold_l)
  309. c----------------------------------------------------------------------
  310. c computing numerical Mach number and its derivatives
  311. c----------------------------------------------------------------------
  312. un_l=uold_l*n_x+vold_l*n_y+wold_l*n_z
  313. un_r=uold_r*n_x+vold_r*n_y+wold_r*n_z
  314. c----------
  315. ut1_l=uold_l*t1_x+vold_l*t1_y+wold_l*t1_z
  316. ut1_r=uold_r*t1_x+vold_r*t1_y+wold_r*t1_z
  317. c----------
  318. ut2_l=uold_l*t2_x+vold_l*t2_y+wold_l*t2_z
  319. ut2_r=uold_r*t2_x+vold_r*t2_y+wold_r*t2_z
  320. c----------------------------------------
  321. ml=un_l/am
  322. mr=un_r/am
  323. mhalf=0.5d0*(un_l+un_r)/am
  324. c---------------------
  325. top=0.5d0*(un_l+un_r)/(am*am)
  326. dmhr_l=-top*damr_l
  327. dmhu_l=n_x/2.0d0/am-top*damu_l
  328. dmhv_l=n_y/2.0d0/am-top*damv_l
  329. dmhw_l=n_z/2.0d0/am-top*damw_l
  330. dmhp_l=-top*damp_l
  331. c---------------------
  332. dmhr_r=-top*damr_r
  333. dmhu_r=n_x/2.0d0/am-top*damu_r
  334. dmhv_r=n_y/2.0d0/am-top*damv_r
  335. dmhw_r=n_z/2.0d0/am-top*damw_r
  336. dmhp_r=-top*damp_r
  337. c--------------------------------
  338. mhalfr=urm/am
  339. c--------------------------------
  340. top=urm/(am*am)
  341. dmfr_l=durr_l/am-top*damr_l
  342. dmfu_l=duru_l/am-top*damu_l
  343. dmfv_l=durv_l/am-top*damv_l
  344. dmfw_l=durw_l/am-top*damw_l
  345. dmfp_l=durp_l/am-top*damp_l
  346. c---------------------
  347. dmfr_r=durr_r/am-top*damr_r
  348. dmfu_r=duru_r/am-top*damu_r
  349. dmfv_r=durv_r/am-top*damv_r
  350. dmfw_r=durw_r/am-top*damw_r
  351. dmfp_r=durp_r/am-top*damp_r
  352. c-------------------------------------------------------------------
  353. c Scaling function for the speed of sound and its derivatives
  354. c see Eq.(32) of the Ref.2)
  355. c-------------------------------------------------------------------
  356. top=(1.0d0-mhalfr*mhalfr)*(1.0d0-mhalfr*mhalfr)
  357. top=top*mhalf*mhalf+4.0d0*mhalfr*mhalfr
  358. bot=1.0d0+mhalfr*mhalfr
  359. if(abs(canc-0.0d0).lt.0.000001d0) then
  360. fmid=1.0d0
  361. else
  362. fmid=sqrt(top)/bot
  363. endif
  364. c--------------------------
  365. temph=-4.0d0*(1.0d0-mhalfr*mhalfr)*mhalfr
  366. temph=temph*mhalf*mhalf+8.0d0*mhalfr
  367. dfm_mf=temph/(sqrt(top)*2.0d0*bot)
  368. dfm_mf=dfm_mf-sqrt(top)*2.0d0*mhalfr/(bot*bot)
  369. c--------------------------
  370. temph=2.0d0*(1.0d0-mhalfr*mhalfr)*(1.0d0-mhalfr*mhalfr)*mhalf
  371. dfm_mh=temph/(2.0d0*bot*sqrt(top))
  372. c--------------------------
  373. dfmr_l=dfm_mf*dmfr_l+dfm_mh*dmhr_l
  374. dfmu_l=dfm_mf*dmfu_l+dfm_mh*dmhu_l
  375. dfmv_l=dfm_mf*dmfv_l+dfm_mh*dmhv_l
  376. dfmw_l=dfm_mf*dmfw_l+dfm_mh*dmhw_l
  377. dfmp_l=dfm_mf*dmfp_l+dfm_mh*dmhp_l
  378. c--------------------------
  379. dfmr_r=dfm_mf*dmfr_r+dfm_mh*dmhr_r
  380. dfmu_r=dfm_mf*dmfu_r+dfm_mh*dmhu_r
  381. dfmv_r=dfm_mf*dmfv_r+dfm_mh*dmhv_r
  382. dfmw_r=dfm_mf*dmfw_r+dfm_mh*dmhw_r
  383. dfmp_r=dfm_mf*dmfp_r+dfm_mh*dmhp_r
  384. c------------------------------------------------------------------
  385. c 'New' speed of sound 'amw' defined as a product of the scaling
  386. c function 'fmid' and the 'Old' speed of sound 'am'; see (31) of Ref.2)
  387. c------------------------------------------------------------------
  388. amw=fmid*am
  389. mlw=un_l/amw
  390. mrw=un_r/amw
  391. c--------------------------
  392. damr_l=canc*coef*dfmr_l*am+fmid*damr_l
  393. damu_l=canc*coef*dfmu_l*am+fmid*damu_l
  394. damv_l=canc*coef*dfmv_l*am+fmid*damv_l
  395. damw_l=canc*coef*dfmw_l*am+fmid*damw_l
  396. damp_l=canc*coef*dfmp_l*am+fmid*damp_l
  397. c--------------------------
  398. damr_r=canc*coef*dfmr_r*am+fmid*damr_r
  399. damu_r=canc*coef*dfmu_r*am+fmid*damu_r
  400. damv_r=canc*coef*dfmv_r*am+fmid*damv_r
  401. damw_r=canc*coef*dfmw_r*am+fmid*damw_r
  402. damp_r=canc*coef*dfmp_r*am+fmid*damp_r
  403. c-/-/-/-/-/-/-/-/-/-/-/-/-/-/-/-/-/
  404. am=amw
  405. c-------------------------------------------------------------------
  406. c Redefinition of the numerical mach numbers;
  407. c see Eqs. (33) and (34) of the Ref.2)
  408. c-------------------------------------------------------------------
  409. if(abs(canc-0.0d0).lt.0.000001d0) then
  410. top=2.0d0
  411. bot=0.0d0
  412. else
  413. top=1.0d0+mhalfr*mhalfr
  414. bot=1.0d0-mhalfr*mhalfr
  415. endif
  416. ml=0.5d0*(top*mlw+bot*mrw)
  417. mr=0.5d0*(top*mrw+bot*mlw)
  418. c-------------------------------------------------------------------
  419. c Mplus and Mmin are calligraphic lettes M+ and M- from the paper,
  420. c see (19a) and (19b), p.367. of the Ref.1)
  421. c-------------------------------------------------------------------
  422. if(abs(ml) .ge. 1.0d0) then
  423. Mplus=(ml+abs(ml))/2.0d0
  424. else
  425. Mplus=(ml+1.0d0)*(ml+1.0d0)/4.0d0
  426. Mplus=Mplus+beta*(ml*ml-1.0d0)*(ml*ml-1.0d0)
  427. endif
  428. Mplus1=(ml+abs(ml))/2.0d0
  429. c-------------------------------------------------------------------
  430. if(abs(mr) .ge. 1.0d0) then
  431. Mmin=(mr-abs(mr))/2.0d0
  432. else
  433. Mmin=-(mr-1.0d0)*(mr-1.0d0)/4.0d0
  434. Mmin=Mmin-beta*(mr*mr-1.0d0)*(mr*mr-1.0d0)
  435. endif
  436. Mmin1=(mr-abs(mr))/2.0d0
  437. c-------------------------------------------------------------------
  438. c---------------------------------------
  439. c Derivatives of ml and mr with respect
  440. c to both sets of primitive variables:
  441. c left and right
  442. c--------------------------------------
  443. temp_l=-un_l/(am*am)
  444. temp_r=-un_r/(am*am)
  445. c--------
  446. dmlr_l=temp_l*damr_l
  447. dmlr_r=temp_l*damr_r
  448. dmrr_l=temp_r*damr_l
  449. dmrr_r=temp_r*damr_r
  450. c--------
  451. dmlu_l=n_x/am+temp_l*damu_l
  452. dmlu_r=temp_l*damu_r
  453. dmru_l=temp_r*damu_l
  454. dmru_r=n_x/am+temp_r*damu_r
  455. c--------
  456. dmlv_l=n_y/am+temp_l*damv_l
  457. dmlv_r=temp_l*damv_r
  458. dmrv_l=temp_r*damv_l
  459. dmrv_r=n_y/am+temp_r*damv_r
  460. c--------
  461. dmlw_l=n_z/am+temp_l*damw_l
  462. dmlw_r=temp_l*damw_r
  463. dmrw_l=temp_r*damw_l
  464. dmrw_r=n_z/am+temp_r*damw_r
  465. c--------
  466. dmlp_l=temp_l*damp_l
  467. dmlp_r=temp_l*damp_r
  468. dmrp_l=temp_r*damp_l
  469. dmrp_r=temp_r*damp_r
  470. c-----------------------------
  471. sr_l=dmlr_l
  472. sr_r=dmlr_r
  473. su_l=dmlu_l
  474. su_r=dmlu_r
  475. sv_l=dmlv_l
  476. sv_r=dmlv_r
  477. sw_l=dmlw_l
  478. sw_r=dmlw_r
  479. sp_l=dmlp_l
  480. sp_r=dmlp_r
  481. c-----------------------------------------------------------------
  482. c Redefinition of the derivatives of the ml & mr
  483. c-----------------------------------------------------------------
  484. temp_l=(mlw-mrw)*mhalfr
  485. temp_r=-temp_l
  486. c--------
  487. dmlr_l=0.5d0*(top*dmlr_l+bot*dmrr_l)+canc*coef*temp_l*dmfr_l
  488. dmlu_l=0.5d0*(top*dmlu_l+bot*dmru_l)+canc*coef*temp_l*dmfu_l
  489. dmlv_l=0.5d0*(top*dmlv_l+bot*dmrv_l)+canc*coef*temp_l*dmfv_l
  490. dmlw_l=0.5d0*(top*dmlw_l+bot*dmrw_l)+canc*coef*temp_l*dmfw_l
  491. dmlp_l=0.5d0*(top*dmlp_l+bot*dmrp_l)+canc*coef*temp_l*dmfp_l
  492. c--------
  493. dmlr_r=0.5d0*(top*dmlr_r+bot*dmrr_r)+canc*coef*temp_l*dmfr_r
  494. dmlu_r=0.5d0*(top*dmlu_r+bot*dmru_r)+canc*coef*temp_l*dmfu_r
  495. dmlv_r=0.5d0*(top*dmlv_r+bot*dmrv_r)+canc*coef*temp_l*dmfv_r
  496. dmlw_r=0.5d0*(top*dmlw_r+bot*dmrw_r)+canc*coef*temp_l*dmfw_r
  497. dmlp_r=0.5d0*(top*dmlp_r+bot*dmrp_r)+canc*coef*temp_l*dmfp_r
  498. c--------
  499. dmrr_l=0.5d0*(top*dmrr_l+bot*sr_l)+canc*coef*temp_r*dmfr_l
  500. dmru_l=0.5d0*(top*dmru_l+bot*su_l)+canc*coef*temp_r*dmfu_l
  501. dmrv_l=0.5d0*(top*dmrv_l+bot*sv_l)+canc*coef*temp_r*dmfv_l
  502. dmrw_l=0.5d0*(top*dmrw_l+bot*sw_l)+canc*coef*temp_r*dmfw_l
  503. dmrp_l=0.5d0*(top*dmrp_l+bot*sp_l)+canc*coef*temp_r*dmfp_l
  504. c--------
  505. dmrr_r=0.5d0*(top*dmrr_r+bot*sr_r)+canc*coef*temp_r*dmfr_r
  506. dmru_r=0.5d0*(top*dmru_r+bot*su_r)+canc*coef*temp_r*dmfu_r
  507. dmrv_r=0.5d0*(top*dmrv_r+bot*sv_r)+canc*coef*temp_r*dmfv_r
  508. dmrw_r=0.5d0*(top*dmrw_r+bot*sw_r)+canc*coef*temp_r*dmfw_r
  509. dmrp_r=0.5d0*(top*dmrp_r+bot*sp_r)+canc*coef*temp_r*dmfp_r
  510. c-----------------------------------------------------------
  511. c mmid is m_{1/2} (notation as in the paper)
  512. c see Eq.(13),p.366 of the Ref.1)
  513. c----------------------------------------------------------
  514. mmid=Mplus+Mmin
  515. c-----------------------------------------------------------
  516. c Computing the derivatives of M+ and M-
  517. c-----------------------------------------------------------
  518. if(ml .ge. 1.0d0) then
  519. dMpr_l=dmlr_l
  520. dMpu_l=dmlu_l
  521. dMpv_l=dmlv_l
  522. dMpw_l=dmlw_l
  523. dMpp_l=dmlp_l
  524. c--------------------
  525. dMpr_r=dmlr_r
  526. dMpu_r=dmlu_r
  527. dMpv_r=dmlv_r
  528. dMpw_r=dmlw_r
  529. dMpp_r=dmlp_r
  530. else
  531. if((ml .gt. -1.0d0) .and. (ml .lt. 1.0d0)) then
  532. temph=(ml+1.0d0)/2.0d0
  533. dMpr_l=(temph+4.0d0*beta*ml*(ml*ml-1.0d0))*dmlr_l
  534. dMpu_l=(temph+4.0d0*beta*ml*(ml*ml-1.0d0))*dmlu_l
  535. dMpv_l=(temph+4.0d0*beta*ml*(ml*ml-1.0d0))*dmlv_l
  536. dMpw_l=(temph+4.0d0*beta*ml*(ml*ml-1.0d0))*dmlw_l
  537. dMpp_l=(temph+4.0d0*beta*ml*(ml*ml-1.0d0))*dmlp_l
  538. c--------------------
  539. dMpr_r=(temph+4.0d0*beta*ml*(ml*ml-1.0d0))*dmlr_r
  540. dMpu_r=(temph+4.0d0*beta*ml*(ml*ml-1.0d0))*dmlu_r
  541. dMpv_r=(temph+4.0d0*beta*ml*(ml*ml-1.0d0))*dmlv_r
  542. dMpw_r=(temph+4.0d0*beta*ml*(ml*ml-1.0d0))*dmlw_r
  543. dMpp_r=(temph+4.0d0*beta*ml*(ml*ml-1.0d0))*dmlp_r
  544. else
  545. dMpr_l=0.0d0
  546. dMpu_l=0.0d0
  547. dMpv_l=0.0d0
  548. dMpw_l=0.0d0
  549. dMpp_l=0.0d0
  550. c---------------------
  551. dMpr_r=0.0d0
  552. dMpu_r=0.0d0
  553. dMpv_r=0.0d0
  554. dMpw_r=0.0d0
  555. dMpp_r=0.0d0
  556. endif
  557. endif
  558. c-----------------------------------------------------------
  559. c addition of low Mach number terms
  560. c-----------------------------------------------------------
  561. if(ml .ge. 0.0d0) then
  562. m1pr_l=dmlr_l
  563. m1pu_l=dmlu_l
  564. m1pv_l=dmlv_l
  565. m1pw_l=dmlw_l
  566. m1pp_l=dmlp_l
  567. c---------------------
  568. m1pr_r=dmlr_r
  569. m1pu_r=dmlu_r
  570. m1pv_r=dmlv_r
  571. m1pw_r=dmlw_r
  572. m1pp_r=dmlp_r
  573. else
  574. m1pr_l=0.0d0
  575. m1pu_l=0.0d0
  576. m1pv_l=0.0d0
  577. m1pw_l=0.0d0
  578. m1pp_l=0.0d0
  579. c--------------------
  580. m1pr_r=0.0d0
  581. m1pu_r=0.0d0
  582. m1pv_r=0.0d0
  583. m1pw_r=0.0d0
  584. m1pp_r=0.0d0
  585. endif
  586. c-----------------------------------------------------------
  587. c-----------------------------------------------------------
  588. if(mr .ge. 1.0d0) then
  589. dMmr_l=0.0d0
  590. dMmu_l=0.0d0
  591. dMmv_l=0.0d0
  592. dMmw_l=0.0d0
  593. dMmp_l=0.0d0
  594. c---------------------
  595. dMmr_r=0.0d0
  596. dMmu_r=0.0d0
  597. dMmv_r=0.0d0
  598. dMmw_r=0.0d0
  599. dMmp_r=0.0d0
  600. else
  601. if((mr .gt. -1.0d0) .and. (mr .lt. 1.0d0)) then
  602. temph=(-mr+1.0d0)/2.0d0
  603. dMmr_l=(temph-4.0d0*beta*mr*(mr*mr-1.0d0))*dmrr_l
  604. dMmu_l=(temph-4.0d0*beta*mr*(mr*mr-1.0d0))*dmru_l
  605. dMmv_l=(temph-4.0d0*beta*mr*(mr*mr-1.0d0))*dmrv_l
  606. dMmw_l=(temph-4.0d0*beta*mr*(mr*mr-1.0d0))*dmrw_l
  607. dMmp_l=(temph-4.0d0*beta*mr*(mr*mr-1.0d0))*dmrp_l
  608. c--------------------
  609. dMmr_r=(temph-4.0d0*beta*mr*(mr*mr-1.0d0))*dmrr_r
  610. dMmu_r=(temph-4.0d0*beta*mr*(mr*mr-1.0d0))*dmru_r
  611. dMmv_r=(temph-4.0d0*beta*mr*(mr*mr-1.0d0))*dmrv_r
  612. dMmw_r=(temph-4.0d0*beta*mr*(mr*mr-1.0d0))*dmrw_r
  613. dMmp_r=(temph-4.0d0*beta*mr*(mr*mr-1.0d0))*dmrp_r
  614. else
  615. dMmr_l=dmrr_l
  616. dMmu_l=dmru_l
  617. dMmv_l=dmrv_l
  618. dMmw_l=dmrw_l
  619. dMmp_l=dmrp_l
  620. c--------------------
  621. dMmr_r=dmrr_r
  622. dMmu_r=dmru_r
  623. dMmv_r=dmrv_r
  624. dMmw_r=dmrw_r
  625. dMmp_r=dmrp_r
  626. endif
  627. endif
  628. c-----------------------------------------------------------
  629. c addition of low Mach number terms
  630. c-----------------------------------------------------------
  631. if(mr .le. 0.0d0) then
  632. m1mr_l=dmrr_l
  633. m1mu_l=dmru_l
  634. m1mv_l=dmrv_l
  635. m1mw_l=dmrw_l
  636. m1mp_l=dmrp_l
  637. c---------------------
  638. m1mr_r=dmrr_r
  639. m1mu_r=dmru_r
  640. m1mv_r=dmrv_r
  641. m1mw_r=dmrw_r
  642. m1mp_r=dmrp_r
  643. else
  644. m1mr_l=0.0d0
  645. m1mu_l=0.0d0
  646. m1mv_l=0.0d0
  647. m1mw_l=0.0d0
  648. m1mp_l=0.0d0
  649. c--------------------
  650. m1mr_r=0.0d0
  651. m1mu_r=0.0d0
  652. m1mv_r=0.0d0
  653. m1mw_r=0.0d0
  654. m1mp_r=0.0d0
  655. endif
  656. c-----------------------------------------------------------------
  657. c computing the derivatives of m_{1/2} (notation as in the paper)
  658. c-----------------------------------------------------------------
  659. dmir_l=dMpr_l+dMmr_l
  660. dmir_r=dMpr_r+dMmr_r
  661. c-------------
  662. dmiu_l=dMpu_l+dMmu_l
  663. dmiu_r=dMpu_r+dMmu_r
  664. c-------------
  665. dmiv_l=dMpv_l+dMmv_l
  666. dmiv_r=dMpv_r+dMmv_r
  667. c-------------
  668. dmiw_l=dMpw_l+dMmw_l
  669. dmiw_r=dMpw_r+dMmw_r
  670. c-------------
  671. dmip_l=dMpp_l+dMmp_l
  672. dmip_r=dMpp_r+dMmp_r
  673. c-------------------------------------------------------------
  674. c computing the main convective variables and their derivatives
  675. c mpl_m is m^{+}_{1/2} (paper's notation) and
  676. c mmin_m is m^{-}_{1/2} (paper's notation),
  677. c see Eq.(A2) on p.370 of the Ref.1)
  678. c----------------------------------------------------------------
  679. termp=(Mmin1-Mmin+Mplus-Mplus1)*(1.0d0/(mhalfr*mhalfr)-1.0d0)
  680. termp=termp*(pold_l-pold_r)/(pold_l/rold_l+pold_r/rold_r)
  681. c-------------------------------------------------------------
  682. c derivatives of the termp
  683. c-------------------------------------------------------------
  684. top=(Mmin1-Mmin+Mplus-Mplus1)
  685. bots=1.0d0/(pold_l/rold_l+pold_r/rold_r)
  686. bot=(pold_l-pold_r)*bots
  687. temph=1.0d0/(mhalfr*mhalfr)-1.0d0
  688. c---------------------------
  689. tmpr_l=(m1mr_l-dMmr_l+dMpr_l-m1pr_l)*bot*temph
  690. tmpr_l=tmpr_l-2.0d0*bot*top*dmfr_l/(mhalfr*mhalfr*mhalfr)
  691. tmpr_l=tmpr_l+temph*top*bot*bots*(pold_l/rold_l/rold_l)
  692. c---------------------------
  693. tmpu_l=(m1mu_l-dMmu_l+dMpu_l-m1pu_l)*bot*temph
  694. tmpu_l=tmpu_l-2.0d0*bot*top*dmfu_l/(mhalfr*mhalfr*mhalfr)
  695. c---------------------------
  696. tmpv_l=(m1mv_l-dMmv_l+dMpv_l-m1pv_l)*bot*temph
  697. tmpv_l=tmpv_l-2.0d0*bot*top*dmfv_l/(mhalfr*mhalfr*mhalfr)
  698. c---------------------------
  699. tmpw_l=(m1mw_l-dMmw_l+dMpw_l-m1pw_l)*bot*temph
  700. tmpw_l=tmpw_l-2.0d0*bot*top*dmfw_l/(mhalfr*mhalfr*mhalfr)
  701. c---------------------------
  702. tmpp_l=(m1mp_l-dMmp_l+dMpp_l-m1pp_l)*bot*temph
  703. tmpp_l=tmpp_l-2.0d0*bot*top*dmfp_l/(mhalfr*mhalfr*mhalfr)
  704. tmpp_l=tmpp_l+temph*top*bots*(1.0d0-bot/rold_l)
  705. c------------rrrrrrrr-------
  706. c------------rrrrrrrr-------
  707. tmpr_r=(m1mr_r-dMmr_r+dMpr_r-m1pr_r)*bot*temph
  708. tmpr_r=tmpr_r-2.0d0*bot*top*dmfr_r/(mhalfr*mhalfr*mhalfr)
  709. tmpr_r=tmpr_r+temph*top*bot*bots*(pold_r/rold_r/rold_r)
  710. c---------------------------
  711. tmpu_r=(m1mu_r-dMmu_r+dMpu_r-m1pu_r)*bot*temph
  712. tmpu_r=tmpu_r-2.0d0*bot*top*dmfu_r/(mhalfr*mhalfr*mhalfr)
  713. c---------------------------
  714. tmpv_r=(m1mv_r-dMmv_r+dMpv_r-m1pv_r)*bot*temph
  715. tmpv_r=tmpv_r-2.0d0*bot*top*dmfv_r/(mhalfr*mhalfr*mhalfr)
  716. c---------------------------
  717. tmpw_r=(m1mw_r-dMmw_r+dMpw_r-m1pw_r)*bot*temph
  718. tmpw_r=tmpw_r-2.0d0*bot*top*dmfw_r/(mhalfr*mhalfr*mhalfr)
  719. c---------------------------
  720. tmpp_r=(m1mp_r-dMmp_r+dMpp_r-m1pp_r)*bot*temph
  721. tmpp_r=tmpp_r-2.0d0*bot*top*dmfp_r/(mhalfr*mhalfr*mhalfr)
  722. tmpp_r=tmpp_r-temph*top*bots*(1.0d0+bot/rold_r)
  723. c-------------------------------------------------------------
  724. if(mmid .ge. 0.0d0) then
  725. mpl_m = mmid
  726. d2mr_l=dmir_l
  727. d2mu_l=dmiu_l
  728. d2mv_l=dmiv_l
  729. d2mw_l=dmiw_l
  730. d2mp_l=dmip_l
  731. c------------
  732. d2mr_r=dmir_r
  733. d2mu_r=dmiu_r
  734. d2mv_r=dmiv_r
  735. d2mw_r=dmiw_r
  736. d2mp_r=dmip_r
  737. c------------
  738. else
  739. mpl_m = 0.0d0
  740. d2mr_l=0.0d0
  741. d2mu_l=0.0d0
  742. d2mv_l=0.0d0
  743. d2mw_l=0.0d0
  744. d2mp_l=0.0d0
  745. c------------
  746. d2mr_r=0.0d0
  747. d2mu_r=0.0d0
  748. d2mv_r=0.0d0
  749. d2mw_r=0.0d0
  750. d2mp_r=0.0d0
  751. endif
  752. c------------------------------------------------------------------
  753. if(mmid .le. 0.0d0) then
  754. mmin_m = mmid
  755. d3mr_l=dmir_l
  756. d3mu_l=dmiu_l
  757. d3mv_l=dmiv_l
  758. d3mw_l=dmiw_l
  759. d3mp_l=dmip_l
  760. c------------
  761. d3mr_r=dmir_r
  762. d3mu_r=dmiu_r
  763. d3mv_r=dmiv_r
  764. d3mw_r=dmiw_r
  765. d3mp_r=dmip_r
  766. c------------
  767. else
  768. mmin_m = 0.0d0
  769. d3mr_l=0.0d0
  770. d3mu_l=0.0d0
  771. d3mv_l=0.0d0
  772. d3mw_l=0.0d0
  773. d3mp_l=0.0d0
  774. c------------
  775. d3mr_r=0.0d0
  776. d3mu_r=0.0d0
  777. d3mv_r=0.0d0
  778. d3mw_r=0.0d0
  779. d3mp_r=0.0d0
  780. endif
  781. c---------------------------------------------------------------
  782. c Computing the calligraphic P+ and P- with their derivatives
  783. c---------------------------------------------------------------
  784. if(ml .ge. 1.0d0) then
  785. Pplus = 1.0d0
  786. else
  787. if((ml .gt. -1.0d0) .and. (ml .lt. 1.0d0)) then
  788. Pplus=(ml+1.0d0)*(ml+1.0d0)*(2.0d0-ml)/4.0d0
  789. Pplus=Pplus+alpha*ml*(ml*ml-1.0d0)*(ml*ml-1.0d0)
  790. else
  791. Pplus = 0.0d0
  792. endif
  793. endif
  794. c---------------------------------------------------------------
  795. if(mr .ge. 1.0d0) then
  796. Pmin = 0.0d0
  797. else
  798. if((mr .gt. -1.0d0) .and. (mr .lt. 1.0d0)) then
  799. Pmin=(mr-1.0d0)*(mr-1.0d0)*(2.0d0+mr)/4.0d0
  800. Pmin=Pmin-alpha*mr*(mr*mr-1.0d0)*(mr*mr-1.0d0)
  801. else
  802. Pmin = 1.0d0
  803. endif
  804. endif
  805. c---------------------------------------------------------------
  806. brac_l=(ml+1.0d0)*(2.0d0-ml)/2.0d0-(ml+1.0d0)*(ml+1.0d0)/4.0d0
  807. brac_l=brac_l+alpha*(ml*ml-1.0d0)*(ml*ml-1.0d0)
  808. brac_l=brac_l+4.0d0*alpha*ml*ml*(ml*ml-1.0d0)
  809. c--------------
  810. brac_r=(mr-1.0d0)*(2.0d0+mr)/2.0d0+(mr-1.0d0)*(mr-1.0d0)/4.0d0
  811. brac_r=brac_r-alpha*(mr*mr-1.0d0)*(mr*mr-1.0d0)
  812. brac_r=brac_r-4.0d0*alpha*mr*mr*(mr*mr-1.0d0)
  813. c---------------------------------------------------------------
  814. if((ml .gt. -1.0d0) .and. (ml .lt. 1.0d0)) then
  815. dPpr_l=brac_l*dmlr_l
  816. dPpr_r=brac_l*dmlr_r
  817. c------------
  818. dPpu_l=brac_l*dmlu_l
  819. dPpu_r=brac_l*dmlu_r
  820. c------------
  821. dPpv_l=brac_l*dmlv_l
  822. dPpv_r=brac_l*dmlv_r
  823. c------------
  824. dPpw_l=brac_l*dmlw_l
  825. dPpw_r=brac_l*dmlw_r
  826. c------------
  827. dPpp_l=brac_l*dmlp_l
  828. dPpp_r=brac_l*dmlp_r
  829. c------------
  830. else
  831. dPpr_l=0.0d0
  832. dPpr_r=0.0d0
  833. c-----------
  834. dPpu_l=0.0d0
  835. dPpu_r=0.0d0
  836. c-----------
  837. dPpv_l=0.0d0
  838. dPpv_r=0.0d0
  839. c-----------
  840. dPpw_l=0.0d0
  841. dPpw_r=0.0d0
  842. c-----------
  843. dPpp_l=0.0d0
  844. dPpp_r=0.0d0
  845. c-----------
  846. endif
  847. c---------------------------------------------------------------
  848. if((mr .gt. -1.0d0) .and. (mr .lt. 1.0d0)) then
  849. dPmr_l=brac_r*dmrr_l
  850. dPmr_r=brac_r*dmrr_r
  851. c------------
  852. dPmu_l=brac_r*dmru_l
  853. dPmu_r=brac_r*dmru_r
  854. c------------
  855. dPmv_l=brac_r*dmrv_l
  856. dPmv_r=brac_r*dmrv_r
  857. c------------
  858. dPmw_l=brac_r*dmrw_l
  859. dPmw_r=brac_r*dmrw_r
  860. c------------
  861. dPmp_l=brac_r*dmrp_l
  862. dPmp_r=brac_r*dmrp_r
  863. c------------
  864. else
  865. dPmr_l=0.0d0
  866. dPmr_r=0.0d0
  867. c-----------
  868. dPmu_l=0.0d0
  869. dPmu_r=0.0d0
  870. c-----------
  871. dPmv_l=0.0d0
  872. dPmv_r=0.0d0
  873. c-----------
  874. dPmw_l=0.0d0
  875. dPmw_r=0.0d0
  876. c-----------
  877. dPmp_l=0.0d0
  878. dPmp_r=0.0d0
  879. c-----------
  880. endif
  881. c---------------------------------------------------------------------
  882. c computing pmid -- p_{1/2} and its derivatives
  883. c---------------------------------------------------------------------
  884. dpir_l=dPpr_l*pold_l+dPmr_l*pold_r
  885. dpiu_l=dPpu_l*pold_l+dPmu_l*pold_r
  886. dpiv_l=dPpv_l*pold_l+dPmv_l*pold_r
  887. dpiw_l=dPpw_l*pold_l+dPmw_l*pold_r
  888. dpip_l=dPpp_l*pold_l+Pplus+dPmp_l*pold_r
  889. c----------------------------
  890. dpir_r=dPpr_r*pold_l+dPmr_r*pold_r
  891. dpiu_r=dPpu_r*pold_l+dPmu_r*pold_r
  892. dpiv_r=dPpv_r*pold_l+dPmv_r*pold_r
  893. dpiw_r=dPpw_r*pold_l+dPmw_r*pold_r
  894. dpip_r=dPpp_r*pold_l+Pmin+dPmp_r*pold_r
  895. c---------------------------------------------------------------------
  896. c !!!!!!!!!!!!!!!!!! JACOBIAN !!!!!!!!!!!!!!!!!!!!!!!!!!!!
  897. c---------------------------------------------------------------------
  898. rum=am*(mpl_m*rold_l+mmin_m*rold_r)+canc*am*termp
  899. c---------------------------------------------------------------------
  900. rumr_l=damr_l*(mpl_m*rold_l+mmin_m*rold_r)+
  901. & am*(d2mr_l*rold_l+mpl_m+d3mr_l*rold_r)
  902. rumr_l=rumr_l+canc*coef*(damr_l*termp+am*tmpr_l)
  903. rumu_l=damu_l*(mpl_m*rold_l+mmin_m*rold_r)+
  904. & am*(d2mu_l*rold_l+d3mu_l*rold_r)
  905. rumu_l=rumu_l+canc*coef*(damu_l*termp+am*tmpu_l)
  906. rumv_l=damv_l*(mpl_m*rold_l+mmin_m*rold_r)+
  907. & am*(d2mv_l*rold_l+d3mv_l*rold_r)
  908. rumv_l=rumv_l+canc*coef*(damv_l*termp+am*tmpv_l)
  909. rumw_l=damw_l*(mpl_m*rold_l+mmin_m*rold_r)+
  910. & am*(d2mw_l*rold_l+d3mw_l*rold_r)
  911. rumw_l=rumw_l+canc*coef*(damw_l*termp+am*tmpw_l)
  912. rump_l=damp_l*(mpl_m*rold_l+mmin_m*rold_r)+
  913. & am*(d2mp_l*rold_l+d3mp_l*rold_r)
  914. rump_l=rump_l+canc*coef*(damp_l*termp+am*tmpp_l)
  915. c-------------------------------------------------
  916. rumr_r=damr_r*(mpl_m*rold_l+mmin_m*rold_r)+
  917. & am*(d2mr_r*rold_l+mmin_m+d3mr_r*rold_r)
  918. rumr_r=rumr_r+canc*coef*(damr_r*termp+am*tmpr_r)
  919. rumu_r=damu_r*(mpl_m*rold_l+mmin_m*rold_r)+
  920. & am*(d2mu_r*rold_l+d3mu_r*rold_r)
  921. rumu_r=rumu_r+canc*coef*(damu_r*termp+am*tmpu_r)
  922. rumv_r=damv_r*(mpl_m*rold_l+mmin_m*rold_r)+
  923. & am*(d2mv_r*rold_l+d3mv_r*rold_r)
  924. rumv_r=rumv_r+canc*coef*(damv_r*termp+am*tmpv_r)
  925. rumw_r=damw_r*(mpl_m*rold_l+mmin_m*rold_r)+
  926. & am*(d2mw_r*rold_l+d3mw_r*rold_r)
  927. rumw_r=rumw_r+canc*coef*(damw_r*termp+am*tmpw_r)
  928. rump_r=damp_r*(mpl_m*rold_l+mmin_m*rold_r)+
  929. & am*(d2mp_r*rold_l+d3mp_r*rold_r)
  930. rump_r=rump_r+canc*coef*(damp_r*termp+am*tmpp_r)
  931. c---------------------------------------------------------------------
  932. c---------------------------------------------------------------------
  933. jl(1,1)=rumr_l
  934. jl(1,2)=rumu_l
  935. jl(1,3)=rumv_l
  936. jl(1,4)=rumw_l
  937. jl(1,5)=rump_l
  938. c------------------------------------
  939. jr(1,1)=rumr_r
  940. jr(1,2)=rumu_r
  941. jr(1,3)=rumv_r
  942. jr(1,4)=rumw_r
  943. jr(1,5)=rump_r
  944. c------------------------------------
  945. ccccc f222222222222 -------------
  946. c derivatives with respect to left set of
  947. c primitive variables
  948. c------------------------------------
  949. if(rum .ge. 0.0d0) then
  950. br1=rum*un_l
  951. br2=rum*ut1_l
  952. br3=rum*ut2_l
  953. else
  954. br1=rum*un_r
  955. br2=rum*ut1_r
  956. br3=rum*ut2_r
  957. endif
  958. c------------------------------------
  959. if(rum .ge. 0.0d0) then
  960. br1=rumr_l*un_l
  961. br2=rumr_l*ut1_l
  962. br3=rumr_l*ut2_l
  963. else
  964. br1=rumr_l*un_r
  965. br2=rumr_l*ut1_r
  966. br3=rumr_l*ut2_r
  967. endif
  968. jl(2,1)=n_x*(br1+dpir_l)+br2*t1_x+br3*t2_x
  969. jl(3,1)=n_y*(br1+dpir_l)+br2*t1_y+br3*t2_y
  970. jl(4,1)=n_z*(br1+dpir_l)+br2*t1_z+br3*t2_z
  971. c-------------------------------------------------------------------------
  972. if(rum .ge. 0.0d0) then
  973. br1=rumu_l*un_l+rum*n_x
  974. br2=rumu_l*ut1_l+rum*t1_x
  975. br3=rumu_l*ut2_l+rum*t2_x
  976. else
  977. br1=rumu_l*un_r
  978. br2=rumu_l*ut1_r
  979. br3=rumu_l*ut2_r
  980. endif
  981. jl(2,2)=n_x*(br1+dpiu_l)+br2*t1_x+br3*t2_x
  982. jl(3,2)=n_y*(br1+dpiu_l)+br2*t1_y+br3*t2_y
  983. jl(4,2)=n_z*(br1+dpiu_l)+br2*t1_z+br3*t2_z
  984. c-------------------------------------------------------------------------
  985. if(rum .ge. 0.0d0) then
  986. br1=rumv_l*un_l+rum*n_y
  987. br2=rumv_l*ut1_l+rum*t1_y
  988. br3=rumv_l*ut2_l+rum*t2_y
  989. else
  990. br1=rumv_l*un_r
  991. br2=rumv_l*ut1_r
  992. br3=rumv_l*ut2_r
  993. endif
  994. jl(2,3)=n_x*(br1+dpiv_l)+br2*t1_x+br3*t2_x
  995. jl(3,3)=n_y*(br1+dpiv_l)+br2*t1_y+br3*t2_y
  996. jl(4,3)=n_z*(br1+dpiv_l)+br2*t1_z+br3*t2_z
  997. c-------------------------------------------------------------------------
  998. if(rum .ge. 0.0d0) then
  999. br1=rumw_l*un_l+rum*n_z
  1000. br2=rumw_l*ut1_l+rum*t1_z
  1001. br3=rumw_l*ut2_l+rum*t2_z
  1002. else
  1003. br1=rumw_l*un_r
  1004. br2=rumw_l*ut1_r
  1005. br3=rumw_l*ut2_r
  1006. endif
  1007. jl(2,4)=n_x*(br1+dpiw_l)+br2*t1_x+br3*t2_x
  1008. jl(3,4)=n_y*(br1+dpiw_l)+br2*t1_y+br3*t2_y
  1009. jl(4,4)=n_z*(br1+dpiw_l)+br2*t1_z+br3*t2_z
  1010. c-------------------------------------------------------------------------
  1011. if(rum .ge. 0.0d0) then
  1012. br1=rump_l*un_l
  1013. br2=rump_l*ut1_l
  1014. br3=rump_l*ut2_l
  1015. else
  1016. br1=rump_l*un_r
  1017. br2=rump_l*ut1_r
  1018. br3=rump_l*ut2_r
  1019. endif
  1020. jl(2,5)=n_x*(br1+dpip_l)+br2*t1_x+br3*t2_x
  1021. jl(3,5)=n_y*(br1+dpip_l)+br2*t1_y+br3*t2_y
  1022. jl(4,5)=n_z*(br1+dpip_l)+br2*t1_z+br3*t2_z
  1023. c-------------------------------------------------------------------------
  1024. c rrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrr
  1025. c-------------------------------------------------------------------------
  1026. if(rum .ge. 0.0d0) then
  1027. br1=rumr_r*un_l
  1028. br2=rumr_r*ut1_l
  1029. br3=rumr_r*ut2_l
  1030. else
  1031. br1=rumr_r*un_r
  1032. br2=rumr_r*ut1_r
  1033. br3=rumr_r*ut2_r
  1034. endif
  1035. jr(2,1)=n_x*(br1+dpir_r)+br2*t1_x+br3*t2_x
  1036. jr(3,1)=n_y*(br1+dpir_r)+br2*t1_y+br3*t2_y
  1037. jr(4,1)=n_z*(br1+dpir_r)+br2*t1_z+br3*t2_z
  1038. c------------------------------------------------------------------------
  1039. if(rum .ge. 0.0d0) then
  1040. br1=rumu_r*un_l
  1041. br2=rumu_r*ut1_l
  1042. br3=rumu_r*ut2_l
  1043. else
  1044. br1=rumu_r*un_r+rum*n_x
  1045. br2=rumu_r*ut1_r+rum*t1_x
  1046. br3=rumu_r*ut2_r+rum*t2_x
  1047. endif
  1048. jr(2,2)=n_x*(br1+dpiu_r)+br2*t1_x+br3*t2_x
  1049. jr(3,2)=n_y*(br1+dpiu_r)+br2*t1_y+br3*t2_y
  1050. jr(4,2)=n_z*(br1+dpiu_r)+br2*t1_z+br3*t2_z
  1051. c------------------------------------------------------------------------
  1052. if(rum .ge. 0.0d0) then
  1053. br1=rumv_r*un_l
  1054. br2=rumv_r*ut1_l
  1055. br3=rumv_r*ut2_l
  1056. else
  1057. br1=rumv_r*un_r+rum*n_y
  1058. br2=rumv_r*ut1_r+rum*t1_y
  1059. br3=rumv_r*ut2_r+rum*t2_y
  1060. endif
  1061. jr(2,3)=n_x*(br1+dpiv_r)+br2*t1_x+br3*t2_x
  1062. jr(3,3)=n_y*(br1+dpiv_r)+br2*t1_y+br3*t2_y
  1063. jr(4,3)=n_z*(br1+dpiv_r)+br2*t1_z+br3*t2_z
  1064. c-----------------------------------------------------------------------
  1065. if(rum .ge. 0.0d0) then
  1066. br1=rumw_r*un_l
  1067. br2=rumw_r*ut1_l
  1068. br3=rumw_r*ut2_l
  1069. else
  1070. br1=rumw_r*un_r+rum*n_z
  1071. br2=rumw_r*ut1_r+rum*t1_z
  1072. br3=rumw_r*ut2_r+rum*t2_z
  1073. endif
  1074. jr(2,4)=n_x*(br1+dpiw_r)+br2*t1_x+br3*t2_x
  1075. jr(3,4)=n_y*(br1+dpiw_r)+br2*t1_y+br3*t2_y
  1076. jr(4,4)=n_z*(br1+dpiw_r)+br2*t1_z+br3*t2_z
  1077. c-----------------------------------------------------------------------
  1078. if(rum .ge. 0.0d0) then
  1079. br1=rump_r*un_l
  1080. br2=rump_r*ut1_l
  1081. br3=rump_r*ut2_l
  1082. else
  1083. br1=rump_r*un_r
  1084. br2=rump_r*ut1_r
  1085. br3=rump_r*ut2_r
  1086. endif
  1087. jr(2,5)=n_x*(br1+dpip_r)+br2*t1_x+br3*t2_x
  1088. jr(3,5)=n_y*(br1+dpip_r)+br2*t1_y+br3*t2_y
  1089. jr(4,5)=n_z*(br1+dpip_r)+br2*t1_z+br3*t2_z
  1090. c-------------------------------------------------------
  1091. c-------------------------------------------------------
  1092. hr_l=(eold_l+pold_l/rold_l)
  1093. hr_r=(eold_r+pold_r/rold_r)
  1094. c-------------------------------------------------
  1095. if(rum .ge. 0.0d0) then
  1096. br1=rumr_l*hr_l-rum*ga*pold_l/gm1/rold_l/rold_l
  1097. else
  1098. br1=rumr_l*hr_r
  1099. endif
  1100. jl(5,1)=br1
  1101. c-------------------------------------------------
  1102. if(rum .ge. 0.0d0) then
  1103. br1=rumu_l*hr_l+rum*uold_l
  1104. else
  1105. br1=rumu_l*hr_r
  1106. endif
  1107. jl(5,2)=br1
  1108. c-------------------------------------------------
  1109. if(rum .ge. 0.0d0) then
  1110. br1=rumv_l*hr_l+rum*vold_l
  1111. else
  1112. br1=rumv_l*hr_r
  1113. endif
  1114. jl(5,3)=br1
  1115. c-------------------------------------------------
  1116. if(rum .ge. 0.0d0) then
  1117. br1=rumw_l*hr_l+rum*wold_l
  1118. else
  1119. br1=rumw_l*hr_r
  1120. endif
  1121. jl(5,4)=br1
  1122. c-------------------------------------------------
  1123. if(rum .ge. 0.0d0) then
  1124. br1=rump_l*hr_l+rum*ga/gm1/rold_l
  1125. else
  1126. br1=rump_l*hr_r
  1127. endif
  1128. jl(5,5)=br1
  1129. c-------------------------------------------------
  1130. c-------------------------------------------------
  1131. if(rum .ge. 0.0d0) then
  1132. br1=rumr_r*hr_l
  1133. else
  1134. br1=rumr_r*hr_r-rum*ga*pold_r/gm1/rold_r/rold_r
  1135. endif
  1136. jr(5,1)=br1
  1137. c---------------------
  1138. if(rum .ge. 0.0d0) then
  1139. br1=rumu_r*hr_l
  1140. else
  1141. br1=rumu_r*hr_r+rum*uold_r
  1142. endif
  1143. jr(5,2)=br1
  1144. c---------------------
  1145. if(rum .ge. 0.0d0) then
  1146. br1=rumv_r*hr_l
  1147. else
  1148. br1=rumv_r*hr_r+rum*vold_r
  1149. endif
  1150. jr(5,3)=br1
  1151. c---------------------
  1152. if(rum .ge. 0.0d0) then
  1153. br1=rumw_r*hr_l
  1154. else
  1155. br1=rumw_r*hr_r+rum*wold_r
  1156. endif
  1157. jr(5,4)=br1
  1158. c---------------------
  1159. if(rum .ge. 0.0d0) then
  1160. br1=rump_r*hr_l
  1161. else
  1162. br1=rump_r*hr_r+rum*ga/gm1/rold_r
  1163. endif
  1164. jr(5,5)=br1
  1165. c---------------------------------------------------
  1166. c----------------------------------------------------------
  1167. c matrixes wwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwww
  1168. c----------------------------------------------------------
  1169. wl(1,1)=1.0d0
  1170. wl(1,2)=0.0d0
  1171. wl(1,3)=0.0d0
  1172. wl(1,4)=0.0d0
  1173. wl(1,5)=0.0d0
  1174. c------------------------------
  1175. wl(2,1)=-uold_l/rold_l
  1176. wl(2,2)=1.0d0/rold_l
  1177. wl(2,3)=0.0d0
  1178. wl(2,4)=0.0d0
  1179. wl(2,5)=0.0d0
  1180. c------------------------------
  1181. wl(3,1)=-vold_l/rold_l
  1182. wl(3,2)=0.0d0
  1183. wl(3,3)=1.0d0/rold_l
  1184. wl(3,4)=0.0d0
  1185. wl(3,5)=0.0d0
  1186. c------------------------------
  1187. wl(4,1)=-wold_l/rold_l
  1188. wl(4,2)=0.0d0
  1189. wl(4,3)=0.0d0
  1190. wl(4,4)=1.0d0/rold_l
  1191. wl(4,5)=0.0d0
  1192. c------------------------------
  1193. wl(5,1)=gm1*(uold_l*uold_l+vold_l*vold_l+wold_l*wold_l)/2.0d0
  1194. wl(5,2)=-uold_l*gm1
  1195. wl(5,3)=-vold_l*gm1
  1196. wl(5,4)=-wold_l*gm1
  1197. wl(5,5)=gm1
  1198. c------------------------------
  1199. c------------------------------
  1200. wr(1,1)=1.0d0
  1201. wr(1,2)=0.0d0
  1202. wr(1,3)=0.0d0
  1203. wr(1,4)=0.0d0
  1204. wr(1,5)=0.0d0
  1205. c------------------------------
  1206. wr(2,1)=-uold_r/rold_r
  1207. wr(2,2)=1.0d0/rold_r
  1208. wr(2,3)=0.0d0
  1209. wr(2,4)=0.0d0
  1210. wr(2,5)=0.0d0
  1211. c------------------------------
  1212. wr(3,1)=-vold_r/rold_r
  1213. wr(3,2)=0.0d0
  1214. wr(3,3)=1.0d0/rold_r
  1215. wr(3,4)=0.0d0
  1216. wr(3,5)=0.0d0
  1217. c------------------------------
  1218. wr(4,1)=-wold_r/rold_r
  1219. wr(4,2)=0.0d0
  1220. wr(4,3)=0.0d0
  1221. wr(4,4)=1.0d0/rold_r
  1222. wr(4,5)=0.0d0
  1223. c------------------------------
  1224. wr(5,1)=gm1*(uold_r*uold_r+vold_r*vold_r+wold_r*wold_r)/2.0d0
  1225. wr(5,2)=-uold_r*gm1
  1226. wr(5,3)=-vold_r*gm1
  1227. wr(5,4)=-wold_r*gm1
  1228. wr(5,5)=gm1
  1229. c----------------------------------------------
  1230. c----------------------------------------------
  1231. do 1 i=1,5
  1232. do 2 j=1,5
  1233. jtl(i,j)=0.0d0
  1234. jtr(i,j)=0.0d0
  1235. do 3 k=1,5
  1236. jtl(i,j)=jtl(i,j)+jl(i,k)*wl(k,j)
  1237. jtr(i,j)=jtr(i,j)+jr(i,k)*wr(k,j)
  1238. 3 continue
  1239. 2 continue
  1240. 1 continue
  1241. c----------------------------------------------
  1242. return
  1243. end
  1244.  
  1245.  
  1246.  
  1247.  
  1248.  
  1249.  
  1250.  
  1251.  
  1252.  
  1253.  
  1254.  
  1255.  
  1256.  
  1257.  
  1258.  
  1259.  
  1260.  
  1261.  
  1262.  
  1263.  
  1264.  

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