Télécharger jw3dbm.eso

Retour à la liste

Numérotation des lignes :

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

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