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-------------------------------------------------------------------
418 : c Mplus and Mmin are calligraphic lettes M+ and M- from the paper,
419 : c see (19a) and (19b), p.367. of the Ref.1)
420 : c-------------------------------------------------------------------
421 : if(abs(ml) .ge. 1.0d0) then
422 : Mplus=(ml+abs(ml))/2.0d0
423 : else
424 : Mplus=(ml+1.0d0)*(ml+1.0d0)/4.0d0
425 : Mplus=Mplus+beta*(ml*ml-1.0d0)*(ml*ml-1.0d0)
426 : endif
427 : Mplus1=(ml+abs(ml))/2.0d0
428 : c-------------------------------------------------------------------
429 : if(abs(mr) .ge. 1.0d0) then
430 : Mmin=(mr-abs(mr))/2.0d0
431 : else
432 : Mmin=-(mr-1.0d0)*(mr-1.0d0)/4.0d0
433 : Mmin=Mmin-beta*(mr*mr-1.0d0)*(mr*mr-1.0d0)
434 : endif
435 : Mmin1=(mr-abs(mr))/2.0d0
436 : c-------------------------------------------------------------------
437 : c---------------------------------------
438 : c Derivatives of ml and mr with respect
439 : c to both sets of primitive variables:
440 : c left and right
441 : c--------------------------------------
442 : temp_l=-un_l/(am*am)
443 : temp_r=-un_r/(am*am)
444 : c--------
445 : dmlr_l=temp_l*damr_l
446 : dmlr_r=temp_l*damr_r
447 : dmrr_l=temp_r*damr_l
448 : dmrr_r=temp_r*damr_r
449 : c--------
450 : dmlu_l=n_x/am+temp_l*damu_l
451 : dmlu_r=temp_l*damu_r
452 : dmru_l=temp_r*damu_l
453 : dmru_r=n_x/am+temp_r*damu_r
454 : c--------
455 : dmlv_l=n_y/am+temp_l*damv_l
456 : dmlv_r=temp_l*damv_r
457 : dmrv_l=temp_r*damv_l
458 : dmrv_r=n_y/am+temp_r*damv_r
459 : c--------
460 : dmlw_l=n_z/am+temp_l*damw_l
461 : dmlw_r=temp_l*damw_r
462 : dmrw_l=temp_r*damw_l
463 : dmrw_r=n_z/am+temp_r*damw_r
464 : c--------
465 : dmlp_l=temp_l*damp_l
466 : dmlp_r=temp_l*damp_r
467 : dmrp_l=temp_r*damp_l
468 : dmrp_r=temp_r*damp_r
469 : c-----------------------------
470 : sr_l=dmlr_l
471 : sr_r=dmlr_r
472 : su_l=dmlu_l
473 : su_r=dmlu_r
474 : sv_l=dmlv_l
475 : sv_r=dmlv_r
476 : sw_l=dmlw_l
477 : sw_r=dmlw_r
478 : sp_l=dmlp_l
479 : sp_r=dmlp_r
480 : c-----------------------------------------------------------------
481 : c Redefinition of the derivatives of the ml & mr
482 : c-----------------------------------------------------------------
483 : temp_l=(mlw-mrw)*mhalfr
484 : temp_r=-temp_l
485 : c--------
486 : dmlr_l=0.5d0*(top*dmlr_l+bot*dmrr_l)+canc*coef*temp_l*dmfr_l
487 : dmlu_l=0.5d0*(top*dmlu_l+bot*dmru_l)+canc*coef*temp_l*dmfu_l
488 : dmlv_l=0.5d0*(top*dmlv_l+bot*dmrv_l)+canc*coef*temp_l*dmfv_l
489 : dmlw_l=0.5d0*(top*dmlw_l+bot*dmrw_l)+canc*coef*temp_l*dmfw_l
490 : dmlp_l=0.5d0*(top*dmlp_l+bot*dmrp_l)+canc*coef*temp_l*dmfp_l
491 : c--------
492 : dmlr_r=0.5d0*(top*dmlr_r+bot*dmrr_r)+canc*coef*temp_l*dmfr_r
493 : dmlu_r=0.5d0*(top*dmlu_r+bot*dmru_r)+canc*coef*temp_l*dmfu_r
494 : dmlv_r=0.5d0*(top*dmlv_r+bot*dmrv_r)+canc*coef*temp_l*dmfv_r
495 : dmlw_r=0.5d0*(top*dmlw_r+bot*dmrw_r)+canc*coef*temp_l*dmfw_r
496 : dmlp_r=0.5d0*(top*dmlp_r+bot*dmrp_r)+canc*coef*temp_l*dmfp_r
497 : c--------
498 : dmrr_l=0.5d0*(top*dmrr_l+bot*sr_l)+canc*coef*temp_r*dmfr_l
499 : dmru_l=0.5d0*(top*dmru_l+bot*su_l)+canc*coef*temp_r*dmfu_l
500 : dmrv_l=0.5d0*(top*dmrv_l+bot*sv_l)+canc*coef*temp_r*dmfv_l
501 : dmrw_l=0.5d0*(top*dmrw_l+bot*sw_l)+canc*coef*temp_r*dmfw_l
502 : dmrp_l=0.5d0*(top*dmrp_l+bot*sp_l)+canc*coef*temp_r*dmfp_l
503 : c--------
504 : dmrr_r=0.5d0*(top*dmrr_r+bot*sr_r)+canc*coef*temp_r*dmfr_r
505 : dmru_r=0.5d0*(top*dmru_r+bot*su_r)+canc*coef*temp_r*dmfu_r
506 : dmrv_r=0.5d0*(top*dmrv_r+bot*sv_r)+canc*coef*temp_r*dmfv_r
507 : dmrw_r=0.5d0*(top*dmrw_r+bot*sw_r)+canc*coef*temp_r*dmfw_r
508 : dmrp_r=0.5d0*(top*dmrp_r+bot*sp_r)+canc*coef*temp_r*dmfp_r
509 : c-----------------------------------------------------------
510 : c mmid is m_{1/2} (notation as in the paper)
511 : c see Eq.(13),p.366 of the Ref.1)
512 : c----------------------------------------------------------
513 : mmid=Mplus+Mmin
514 : c-----------------------------------------------------------
515 : c Computing the derivatives of M+ and M-
516 : c-----------------------------------------------------------
517 : if(ml .ge. 1.0d0) then
518 : dMpr_l=dmlr_l
519 : dMpu_l=dmlu_l
520 : dMpv_l=dmlv_l
521 : dMpw_l=dmlw_l
522 : dMpp_l=dmlp_l
523 : c--------------------
524 : dMpr_r=dmlr_r
525 : dMpu_r=dmlu_r
526 : dMpv_r=dmlv_r
527 : dMpw_r=dmlw_r
528 : dMpp_r=dmlp_r
529 : else
530 : if((ml .gt. -1.0d0) .and. (ml .lt. 1.0d0)) then
531 : temph=(ml+1.0d0)/2.0d0
532 : dMpr_l=(temph+4.0d0*beta*ml*(ml*ml-1.0d0))*dmlr_l
533 : dMpu_l=(temph+4.0d0*beta*ml*(ml*ml-1.0d0))*dmlu_l
534 : dMpv_l=(temph+4.0d0*beta*ml*(ml*ml-1.0d0))*dmlv_l
535 : dMpw_l=(temph+4.0d0*beta*ml*(ml*ml-1.0d0))*dmlw_l
536 : dMpp_l=(temph+4.0d0*beta*ml*(ml*ml-1.0d0))*dmlp_l
537 : c--------------------
538 : dMpr_r=(temph+4.0d0*beta*ml*(ml*ml-1.0d0))*dmlr_r
539 : dMpu_r=(temph+4.0d0*beta*ml*(ml*ml-1.0d0))*dmlu_r
540 : dMpv_r=(temph+4.0d0*beta*ml*(ml*ml-1.0d0))*dmlv_r
541 : dMpw_r=(temph+4.0d0*beta*ml*(ml*ml-1.0d0))*dmlw_r
542 : dMpp_r=(temph+4.0d0*beta*ml*(ml*ml-1.0d0))*dmlp_r
543 : else
544 : dMpr_l=0.0d0
545 : dMpu_l=0.0d0
546 : dMpv_l=0.0d0
547 : dMpw_l=0.0d0
548 : dMpp_l=0.0d0
549 : c---------------------
550 : dMpr_r=0.0d0
551 : dMpu_r=0.0d0
552 : dMpv_r=0.0d0
553 : dMpw_r=0.0d0
554 : dMpp_r=0.0d0
555 : endif
556 : endif
557 : c-----------------------------------------------------------
558 : c addition of low Mach number terms
559 : c-----------------------------------------------------------
560 : if(ml .ge. 0.0d0) then
561 : m1pr_l=dmlr_l
562 : m1pu_l=dmlu_l
563 : m1pv_l=dmlv_l
564 : m1pw_l=dmlw_l
565 : m1pp_l=dmlp_l
566 : c---------------------
567 : m1pr_r=dmlr_r
568 : m1pu_r=dmlu_r
569 : m1pv_r=dmlv_r
570 : m1pw_r=dmlw_r
571 : m1pp_r=dmlp_r
572 : else
573 : m1pr_l=0.0d0
574 : m1pu_l=0.0d0
575 : m1pv_l=0.0d0
576 : m1pw_l=0.0d0
577 : m1pp_l=0.0d0
578 : c--------------------
579 : m1pr_r=0.0d0
580 : m1pu_r=0.0d0
581 : m1pv_r=0.0d0
582 : m1pw_r=0.0d0
583 : m1pp_r=0.0d0
584 : endif
585 : c-----------------------------------------------------------
586 : c-----------------------------------------------------------
587 : if(mr .ge. 1.0d0) then
588 : dMmr_l=0.0d0
589 : dMmu_l=0.0d0
590 : dMmv_l=0.0d0
591 : dMmw_l=0.0d0
592 : dMmp_l=0.0d0
593 : c---------------------
594 : dMmr_r=0.0d0
595 : dMmu_r=0.0d0
596 : dMmv_r=0.0d0
597 : dMmw_r=0.0d0
598 : dMmp_r=0.0d0
599 : else
600 : if((mr .gt. -1.0d0) .and. (mr .lt. 1.0d0)) then
601 : temph=(-mr+1.0d0)/2.0d0
602 : dMmr_l=(temph-4.0d0*beta*mr*(mr*mr-1.0d0))*dmrr_l
603 : dMmu_l=(temph-4.0d0*beta*mr*(mr*mr-1.0d0))*dmru_l
604 : dMmv_l=(temph-4.0d0*beta*mr*(mr*mr-1.0d0))*dmrv_l
605 : dMmw_l=(temph-4.0d0*beta*mr*(mr*mr-1.0d0))*dmrw_l
606 : dMmp_l=(temph-4.0d0*beta*mr*(mr*mr-1.0d0))*dmrp_l
607 : c--------------------
608 : dMmr_r=(temph-4.0d0*beta*mr*(mr*mr-1.0d0))*dmrr_r
609 : dMmu_r=(temph-4.0d0*beta*mr*(mr*mr-1.0d0))*dmru_r
610 : dMmv_r=(temph-4.0d0*beta*mr*(mr*mr-1.0d0))*dmrv_r
611 : dMmw_r=(temph-4.0d0*beta*mr*(mr*mr-1.0d0))*dmrw_r
612 : dMmp_r=(temph-4.0d0*beta*mr*(mr*mr-1.0d0))*dmrp_r
613 : else
614 : dMmr_l=dmrr_l
615 : dMmu_l=dmru_l
616 : dMmv_l=dmrv_l
617 : dMmw_l=dmrw_l
618 : dMmp_l=dmrp_l
619 : c--------------------
620 : dMmr_r=dmrr_r
621 : dMmu_r=dmru_r
622 : dMmv_r=dmrv_r
623 : dMmw_r=dmrw_r
624 : dMmp_r=dmrp_r
625 : endif
626 : endif
627 : c-----------------------------------------------------------
628 : c addition of low Mach number terms
629 : c-----------------------------------------------------------
630 : if(mr .le. 0.0d0) then
631 : m1mr_l=dmrr_l
632 : m1mu_l=dmru_l
633 : m1mv_l=dmrv_l
634 : m1mw_l=dmrw_l
635 : m1mp_l=dmrp_l
636 : c---------------------
637 : m1mr_r=dmrr_r
638 : m1mu_r=dmru_r
639 : m1mv_r=dmrv_r
640 : m1mw_r=dmrw_r
641 : m1mp_r=dmrp_r
642 : else
643 : m1mr_l=0.0d0
644 : m1mu_l=0.0d0
645 : m1mv_l=0.0d0
646 : m1mw_l=0.0d0
647 : m1mp_l=0.0d0
648 : c--------------------
649 : m1mr_r=0.0d0
650 : m1mu_r=0.0d0
651 : m1mv_r=0.0d0
652 : m1mw_r=0.0d0
653 : m1mp_r=0.0d0
654 : endif
655 : c-----------------------------------------------------------------
656 : c computing the derivatives of m_{1/2} (notation as in the paper)
657 : c-----------------------------------------------------------------
658 : dmir_l=dMpr_l+dMmr_l
659 : dmir_r=dMpr_r+dMmr_r
660 : c-------------
661 : dmiu_l=dMpu_l+dMmu_l
662 : dmiu_r=dMpu_r+dMmu_r
663 : c-------------
664 : dmiv_l=dMpv_l+dMmv_l
665 : dmiv_r=dMpv_r+dMmv_r
666 : c-------------
667 : dmiw_l=dMpw_l+dMmw_l
668 : dmiw_r=dMpw_r+dMmw_r
669 : c-------------
670 : dmip_l=dMpp_l+dMmp_l
671 : dmip_r=dMpp_r+dMmp_r
672 : c-------------------------------------------------------------
673 : c computing the main convective variables and their derivatives
674 : c mpl_m is m^{+}_{1/2} (paper's notation) and
675 : c mmin_m is m^{-}_{1/2} (paper's notation),
676 : c see Eq.(A2) on p.370 of the Ref.1)
677 : c----------------------------------------------------------------
678 : termp=(Mmin1-Mmin+Mplus-Mplus1)*(1.0d0/(mhalfr*mhalfr)-1.0d0)
679 : termp=termp*(pold_l-pold_r)/(pold_l/rold_l+pold_r/rold_r)
680 : c-------------------------------------------------------------
681 : c derivatives of the termp
682 : c-------------------------------------------------------------
683 : top=(Mmin1-Mmin+Mplus-Mplus1)
684 : bots=1.0d0/(pold_l/rold_l+pold_r/rold_r)
685 : bot=(pold_l-pold_r)*bots
686 : temph=1.0d0/(mhalfr*mhalfr)-1.0d0
687 : c---------------------------
688 : tmpr_l=(m1mr_l-dMmr_l+dMpr_l-m1pr_l)*bot*temph
689 : tmpr_l=tmpr_l-2.0d0*bot*top*dmfr_l/(mhalfr*mhalfr*mhalfr)
690 : tmpr_l=tmpr_l+temph*top*bot*bots*(pold_l/rold_l/rold_l)
691 : c---------------------------
692 : tmpu_l=(m1mu_l-dMmu_l+dMpu_l-m1pu_l)*bot*temph
693 : tmpu_l=tmpu_l-2.0d0*bot*top*dmfu_l/(mhalfr*mhalfr*mhalfr)
694 : c---------------------------
695 : tmpv_l=(m1mv_l-dMmv_l+dMpv_l-m1pv_l)*bot*temph
696 : tmpv_l=tmpv_l-2.0d0*bot*top*dmfv_l/(mhalfr*mhalfr*mhalfr)
697 : c---------------------------
698 : tmpw_l=(m1mw_l-dMmw_l+dMpw_l-m1pw_l)*bot*temph
699 : tmpw_l=tmpw_l-2.0d0*bot*top*dmfw_l/(mhalfr*mhalfr*mhalfr)
700 : c---------------------------
701 : tmpp_l=(m1mp_l-dMmp_l+dMpp_l-m1pp_l)*bot*temph
702 : tmpp_l=tmpp_l-2.0d0*bot*top*dmfp_l/(mhalfr*mhalfr*mhalfr)
703 : tmpp_l=tmpp_l+temph*top*bots*(1.0d0-bot/rold_l)
704 : c------------rrrrrrrr-------
705 : c------------rrrrrrrr-------
706 : tmpr_r=(m1mr_r-dMmr_r+dMpr_r-m1pr_r)*bot*temph
707 : tmpr_r=tmpr_r-2.0d0*bot*top*dmfr_r/(mhalfr*mhalfr*mhalfr)
708 : tmpr_r=tmpr_r+temph*top*bot*bots*(pold_r/rold_r/rold_r)
709 : c---------------------------
710 : tmpu_r=(m1mu_r-dMmu_r+dMpu_r-m1pu_r)*bot*temph
711 : tmpu_r=tmpu_r-2.0d0*bot*top*dmfu_r/(mhalfr*mhalfr*mhalfr)
712 : c---------------------------
713 : tmpv_r=(m1mv_r-dMmv_r+dMpv_r-m1pv_r)*bot*temph
714 : tmpv_r=tmpv_r-2.0d0*bot*top*dmfv_r/(mhalfr*mhalfr*mhalfr)
715 : c---------------------------
716 : tmpw_r=(m1mw_r-dMmw_r+dMpw_r-m1pw_r)*bot*temph
717 : tmpw_r=tmpw_r-2.0d0*bot*top*dmfw_r/(mhalfr*mhalfr*mhalfr)
718 : c---------------------------
719 : tmpp_r=(m1mp_r-dMmp_r+dMpp_r-m1pp_r)*bot*temph
720 : tmpp_r=tmpp_r-2.0d0*bot*top*dmfp_r/(mhalfr*mhalfr*mhalfr)
721 : tmpp_r=tmpp_r-temph*top*bots*(1.0d0+bot/rold_r)
722 : c-------------------------------------------------------------
723 : if(mmid .ge. 0.0d0) then
724 : mpl_m = mmid
725 : d2mr_l=dmir_l
726 : d2mu_l=dmiu_l
727 : d2mv_l=dmiv_l
728 : d2mw_l=dmiw_l
729 : d2mp_l=dmip_l
730 : c------------
731 : d2mr_r=dmir_r
732 : d2mu_r=dmiu_r
733 : d2mv_r=dmiv_r
734 : d2mw_r=dmiw_r
735 : d2mp_r=dmip_r
736 : c------------
737 : else
738 : mpl_m = 0.0d0
739 : d2mr_l=0.0d0
740 : d2mu_l=0.0d0
741 : d2mv_l=0.0d0
742 : d2mw_l=0.0d0
743 : d2mp_l=0.0d0
744 : c------------
745 : d2mr_r=0.0d0
746 : d2mu_r=0.0d0
747 : d2mv_r=0.0d0
748 : d2mw_r=0.0d0
749 : d2mp_r=0.0d0
750 : endif
751 : c------------------------------------------------------------------
752 : if(mmid .le. 0.0d0) then
753 : mmin_m = mmid
754 : d3mr_l=dmir_l
755 : d3mu_l=dmiu_l
756 : d3mv_l=dmiv_l
757 : d3mw_l=dmiw_l
758 : d3mp_l=dmip_l
759 : c------------
760 : d3mr_r=dmir_r
761 : d3mu_r=dmiu_r
762 : d3mv_r=dmiv_r
763 : d3mw_r=dmiw_r
764 : d3mp_r=dmip_r
765 : c------------
766 : else
767 : mmin_m = 0.0d0
768 : d3mr_l=0.0d0
769 : d3mu_l=0.0d0
770 : d3mv_l=0.0d0
771 : d3mw_l=0.0d0
772 : d3mp_l=0.0d0
773 : c------------
774 : d3mr_r=0.0d0
775 : d3mu_r=0.0d0
776 : d3mv_r=0.0d0
777 : d3mw_r=0.0d0
778 : d3mp_r=0.0d0
779 : endif
780 : c---------------------------------------------------------------
780 : c Computing the calligraphic P+ and P- with their derivatives
781 : c---------------------------------------------------------------
782 : if(ml .ge. 1.0d0) then
783 : Pplus = 1.0d0
784 : else
785 : if((ml .gt. -1.0d0) .and. (ml .lt. 1.0d0)) then
786 : Pplus=(ml+1.0d0)*(ml+1.0d0)*(2.0d0-ml)/4.0d0
787 : Pplus=Pplus+alpha*ml*(ml*ml-1.0d0)*(ml*ml-1.0d0)
788 : else
789 : Pplus = 0.0d0
790 : endif
791 : endif
792 : c---------------------------------------------------------------
793 : if(mr .ge. 1.0d0) then
794 : Pmin = 0.0d0
795 : else
796 : if((mr .gt. -1.0d0) .and. (mr .lt. 1.0d0)) then
797 : Pmin=(mr-1.0d0)*(mr-1.0d0)*(2.0d0+mr)/4.0d0
798 : Pmin=Pmin-alpha*mr*(mr*mr-1.0d0)*(mr*mr-1.0d0)
799 : else
800 : Pmin = 1.0d0
801 : endif
802 : endif
803 : c---------------------------------------------------------------
804 : brac_l=(ml+1.0d0)*(2.0d0-ml)/2.0d0-(ml+1.0d0)*(ml+1.0d0)/4.0d0
805 : brac_l=brac_l+alpha*(ml*ml-1.0d0)*(ml*ml-1.0d0)
806 : brac_l=brac_l+4.0d0*alpha*ml*ml*(ml*ml-1.0d0)
807 : c--------------
808 : brac_r=(mr-1.0d0)*(2.0d0+mr)/2.0d0+(mr-1.0d0)*(mr-1.0d0)/4.0d0
809 : brac_r=brac_r-alpha*(mr*mr-1.0d0)*(mr*mr-1.0d0)
810 : brac_r=brac_r-4.0d0*alpha*mr*mr*(mr*mr-1.0d0)
811 : c---------------------------------------------------------------
812 : if((ml .gt. -1.0d0) .and. (ml .lt. 1.0d0)) then
813 : dPpr_l=brac_l*dmlr_l
814 : dPpr_r=brac_l*dmlr_r
815 : c------------
816 : dPpu_l=brac_l*dmlu_l
817 : dPpu_r=brac_l*dmlu_r
818 : c------------
819 : dPpv_l=brac_l*dmlv_l
820 : dPpv_r=brac_l*dmlv_r
821 : c------------
822 : dPpw_l=brac_l*dmlw_l
823 : dPpw_r=brac_l*dmlw_r
824 : c------------
825 : dPpp_l=brac_l*dmlp_l
826 : dPpp_r=brac_l*dmlp_r
827 : c------------
828 : else
829 : dPpr_l=0.0d0
830 : dPpr_r=0.0d0
831 : c-----------
832 : dPpu_l=0.0d0
833 : dPpu_r=0.0d0
834 : c-----------
835 : dPpv_l=0.0d0
836 : dPpv_r=0.0d0
837 : c-----------
838 : dPpw_l=0.0d0
839 : dPpw_r=0.0d0
840 : c-----------
841 : dPpp_l=0.0d0
842 : dPpp_r=0.0d0
843 : c-----------
844 : endif
845 : c---------------------------------------------------------------
846 : if((mr .gt. -1.0d0) .and. (mr .lt. 1.0d0)) then
847 : dPmr_l=brac_r*dmrr_l
848 : dPmr_r=brac_r*dmrr_r
849 : c------------
850 : dPmu_l=brac_r*dmru_l
851 : dPmu_r=brac_r*dmru_r
852 : c------------
853 : dPmv_l=brac_r*dmrv_l
854 : dPmv_r=brac_r*dmrv_r
855 : c------------
856 : dPmw_l=brac_r*dmrw_l
857 : dPmw_r=brac_r*dmrw_r
858 : c------------
859 : dPmp_l=brac_r*dmrp_l
860 : dPmp_r=brac_r*dmrp_r
861 : c------------
862 : else
863 : dPmr_l=0.0d0
864 : dPmr_r=0.0d0
865 : c-----------
866 : dPmu_l=0.0d0
867 : dPmu_r=0.0d0
868 : c-----------
869 : dPmv_l=0.0d0
870 : dPmv_r=0.0d0
871 : c-----------
872 : dPmw_l=0.0d0
873 : dPmw_r=0.0d0
874 : c-----------
875 : dPmp_l=0.0d0
876 : dPmp_r=0.0d0
877 : c-----------
878 : endif
879 : c---------------------------------------------------------------------
880 : c computing pmid -- p_{1/2} and its derivatives
881 : c---------------------------------------------------------------------
882 : dpir_l=dPpr_l*pold_l+dPmr_l*pold_r
883 : dpiu_l=dPpu_l*pold_l+dPmu_l*pold_r
884 : dpiv_l=dPpv_l*pold_l+dPmv_l*pold_r
885 : dpiw_l=dPpw_l*pold_l+dPmw_l*pold_r
886 : dpip_l=dPpp_l*pold_l+Pplus+dPmp_l*pold_r
887 : c----------------------------
888 : dpir_r=dPpr_r*pold_l+dPmr_r*pold_r
889 : dpiu_r=dPpu_r*pold_l+dPmu_r*pold_r
890 : dpiv_r=dPpv_r*pold_l+dPmv_r*pold_r
891 : dpiw_r=dPpw_r*pold_l+dPmw_r*pold_r
892 : dpip_r=dPpp_r*pold_l+Pmin+dPmp_r*pold_r
893 : c---------------------------------------------------------------------
894 : c !!!!!!!!!!!!!!!!!! JACOBIAN !!!!!!!!!!!!!!!!!!!!!!!!!!!!
895 : c---------------------------------------------------------------------
896 : rum=am*(mpl_m*rold_l+mmin_m*rold_r)+canc*am*termp
897 : c---------------------------------------------------------------------
898 : rumr_l=damr_l*(mpl_m*rold_l+mmin_m*rold_r)+
899 : & am*(d2mr_l*rold_l+mpl_m+d3mr_l*rold_r)
900 : rumr_l=rumr_l+canc*coef*(damr_l*termp+am*tmpr_l)
901 : rumu_l=damu_l*(mpl_m*rold_l+mmin_m*rold_r)+
902 : & am*(d2mu_l*rold_l+d3mu_l*rold_r)
903 : rumu_l=rumu_l+canc*coef*(damu_l*termp+am*tmpu_l)
904 : rumv_l=damv_l*(mpl_m*rold_l+mmin_m*rold_r)+
905 : & am*(d2mv_l*rold_l+d3mv_l*rold_r)
906 : rumv_l=rumv_l+canc*coef*(damv_l*termp+am*tmpv_l)
907 : rumw_l=damw_l*(mpl_m*rold_l+mmin_m*rold_r)+
908 : & am*(d2mw_l*rold_l+d3mw_l*rold_r)
909 : rumw_l=rumw_l+canc*coef*(damw_l*termp+am*tmpw_l)
910 : rump_l=damp_l*(mpl_m*rold_l+mmin_m*rold_r)+
911 : & am*(d2mp_l*rold_l+d3mp_l*rold_r)
912 : rump_l=rump_l+canc*coef*(damp_l*termp+am*tmpp_l)
913 : c-------------------------------------------------
914 : rumr_r=damr_r*(mpl_m*rold_l+mmin_m*rold_r)+
915 : & am*(d2mr_r*rold_l+mmin_m+d3mr_r*rold_r)
916 : rumr_r=rumr_r+canc*coef*(damr_r*termp+am*tmpr_r)
917 : rumu_r=damu_r*(mpl_m*rold_l+mmin_m*rold_r)+
918 : & am*(d2mu_r*rold_l+d3mu_r*rold_r)
919 : rumu_r=rumu_r+canc*coef*(damu_r*termp+am*tmpu_r)
920 : rumv_r=damv_r*(mpl_m*rold_l+mmin_m*rold_r)+
921 : & am*(d2mv_r*rold_l+d3mv_r*rold_r)
922 : rumv_r=rumv_r+canc*coef*(damv_r*termp+am*tmpv_r)
923 : rumw_r=damw_r*(mpl_m*rold_l+mmin_m*rold_r)+
924 : & am*(d2mw_r*rold_l+d3mw_r*rold_r)
925 : rumw_r=rumw_r+canc*coef*(damw_r*termp+am*tmpw_r)
926 : rump_r=damp_r*(mpl_m*rold_l+mmin_m*rold_r)+
927 : & am*(d2mp_r*rold_l+d3mp_r*rold_r)
928 : rump_r=rump_r+canc*coef*(damp_r*termp+am*tmpp_r)
929 : c---------------------------------------------------------------------
930 : c---------------------------------------------------------------------
931 : jl(1,1)=rumr_l
932 : jl(1,2)=rumu_l
933 : jl(1,3)=rumv_l
934 : jl(1,4)=rumw_l
935 : jl(1,5)=rump_l
936 : c------------------------------------
937 : jr(1,1)=rumr_r
938 : jr(1,2)=rumu_r
939 : jr(1,3)=rumv_r
940 : jr(1,4)=rumw_r
941 : jr(1,5)=rump_r
942 : c------------------------------------
943 : ccccc f222222222222 -------------
944 : c derivatives with respect to left set of
945 : c primitive variables
946 : c------------------------------------
947 : if(rum .ge. 0.0d0) then
948 : br1=rum*un_l
949 : br2=rum*ut1_l
950 : br3=rum*ut2_l
951 : else
952 : br1=rum*un_r
953 : br2=rum*ut1_r
954 : br3=rum*ut2_r
955 : endif
956 : c------------------------------------
957 : if(rum .ge. 0.0d0) then
958 : br1=rumr_l*un_l
959 : br2=rumr_l*ut1_l
960 : br3=rumr_l*ut2_l
961 : else
962 : br1=rumr_l*un_r
963 : br2=rumr_l*ut1_r
964 : br3=rumr_l*ut2_r
965 : endif
966 : jl(2,1)=n_x*(br1+dpir_l)+br2*t1_x+br3*t2_x
967 : jl(3,1)=n_y*(br1+dpir_l)+br2*t1_y+br3*t2_y
968 : jl(4,1)=n_z*(br1+dpir_l)+br2*t1_z+br3*t2_z
969 : c-------------------------------------------------------------------------
970 : if(rum .ge. 0.0d0) then
971 : br1=rumu_l*un_l+rum*n_x
972 : br2=rumu_l*ut1_l+rum*t1_x
973 : br3=rumu_l*ut2_l+rum*t2_x
974 : else
975 : br1=rumu_l*un_r
976 : br2=rumu_l*ut1_r
977 : br3=rumu_l*ut2_r
978 : endif
979 : jl(2,2)=n_x*(br1+dpiu_l)+br2*t1_x+br3*t2_x
980 : jl(3,2)=n_y*(br1+dpiu_l)+br2*t1_y+br3*t2_y
981 : jl(4,2)=n_z*(br1+dpiu_l)+br2*t1_z+br3*t2_z
982 : c-------------------------------------------------------------------------
983 : if(rum .ge. 0.0d0) then
984 : br1=rumv_l*un_l+rum*n_y
985 : br2=rumv_l*ut1_l+rum*t1_y
986 : br3=rumv_l*ut2_l+rum*t2_y
987 : else
988 : br1=rumv_l*un_r
989 : br2=rumv_l*ut1_r
990 : br3=rumv_l*ut2_r
991 : endif
992 : jl(2,3)=n_x*(br1+dpiv_l)+br2*t1_x+br3*t2_x
993 : jl(3,3)=n_y*(br1+dpiv_l)+br2*t1_y+br3*t2_y
994 : jl(4,3)=n_z*(br1+dpiv_l)+br2*t1_z+br3*t2_z
995 : c-------------------------------------------------------------------------
996 : if(rum .ge. 0.0d0) then
997 : br1=rumw_l*un_l+rum*n_z
998 : br2=rumw_l*ut1_l+rum*t1_z
999 : br3=rumw_l*ut2_l+rum*t2_z
1000 : else
1001 : br1=rumw_l*un_r
1002 : br2=rumw_l*ut1_r
1003 : br3=rumw_l*ut2_r
1004 : endif
1005 : jl(2,4)=n_x*(br1+dpiw_l)+br2*t1_x+br3*t2_x
1006 : jl(3,4)=n_y*(br1+dpiw_l)+br2*t1_y+br3*t2_y
1007 : jl(4,4)=n_z*(br1+dpiw_l)+br2*t1_z+br3*t2_z
1008 : c-------------------------------------------------------------------------
1009 : if(rum .ge. 0.0d0) then
1010 : br1=rump_l*un_l
1011 : br2=rump_l*ut1_l
1012 : br3=rump_l*ut2_l
1013 : else
1014 : br1=rump_l*un_r
1015 : br2=rump_l*ut1_r
1016 : br3=rump_l*ut2_r
1017 : endif
1018 : jl(2,5)=n_x*(br1+dpip_l)+br2*t1_x+br3*t2_x
1019 : jl(3,5)=n_y*(br1+dpip_l)+br2*t1_y+br3*t2_y
1020 : jl(4,5)=n_z*(br1+dpip_l)+br2*t1_z+br3*t2_z
1021 : c-------------------------------------------------------------------------
1022 : c rrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrr
1023 : c-------------------------------------------------------------------------
1024 : if(rum .ge. 0.0d0) then
1025 : br1=rumr_r*un_l
1026 : br2=rumr_r*ut1_l
1027 : br3=rumr_r*ut2_l
1028 : else
1029 : br1=rumr_r*un_r
1030 : br2=rumr_r*ut1_r
1031 : br3=rumr_r*ut2_r
1032 : endif
1033 : jr(2,1)=n_x*(br1+dpir_r)+br2*t1_x+br3*t2_x
1034 : jr(3,1)=n_y*(br1+dpir_r)+br2*t1_y+br3*t2_y
1035 : jr(4,1)=n_z*(br1+dpir_r)+br2*t1_z+br3*t2_z
1036 : c------------------------------------------------------------------------
1037 : if(rum .ge. 0.0d0) then
1038 : br1=rumu_r*un_l
1039 : br2=rumu_r*ut1_l
1040 : br3=rumu_r*ut2_l
1041 : else
1042 : br1=rumu_r*un_r+rum*n_x
1043 : br2=rumu_r*ut1_r+rum*t1_x
1044 : br3=rumu_r*ut2_r+rum*t2_x
1045 : endif
1046 : jr(2,2)=n_x*(br1+dpiu_r)+br2*t1_x+br3*t2_x
1047 : jr(3,2)=n_y*(br1+dpiu_r)+br2*t1_y+br3*t2_y
1048 : jr(4,2)=n_z*(br1+dpiu_r)+br2*t1_z+br3*t2_z
1049 : c------------------------------------------------------------------------
1050 : if(rum .ge. 0.0d0) then
1051 : br1=rumv_r*un_l
1052 : br2=rumv_r*ut1_l
1053 : br3=rumv_r*ut2_l
1054 : else
1055 : br1=rumv_r*un_r+rum*n_y
1056 : br2=rumv_r*ut1_r+rum*t1_y
1057 : br3=rumv_r*ut2_r+rum*t2_y
1058 : endif
1059 : jr(2,3)=n_x*(br1+dpiv_r)+br2*t1_x+br3*t2_x
1060 : jr(3,3)=n_y*(br1+dpiv_r)+br2*t1_y+br3*t2_y
1061 : jr(4,3)=n_z*(br1+dpiv_r)+br2*t1_z+br3*t2_z
1062 : c-----------------------------------------------------------------------
1063 : if(rum .ge. 0.0d0) then
1064 : br1=rumw_r*un_l
1065 : br2=rumw_r*ut1_l
1066 : br3=rumw_r*ut2_l
1067 : else
1068 : br1=rumw_r*un_r+rum*n_z
1069 : br2=rumw_r*ut1_r+rum*t1_z
1070 : br3=rumw_r*ut2_r+rum*t2_z
1071 : endif
1072 : jr(2,4)=n_x*(br1+dpiw_r)+br2*t1_x+br3*t2_x
1073 : jr(3,4)=n_y*(br1+dpiw_r)+br2*t1_y+br3*t2_y
1074 : jr(4,4)=n_z*(br1+dpiw_r)+br2*t1_z+br3*t2_z
1075 : c-----------------------------------------------------------------------
1076 : if(rum .ge. 0.0d0) then
1077 : br1=rump_r*un_l
1078 : br2=rump_r*ut1_l
1079 : br3=rump_r*ut2_l
1080 : else
1081 : br1=rump_r*un_r
1082 : br2=rump_r*ut1_r
1083 : br3=rump_r*ut2_r
1084 : endif
1085 : jr(2,5)=n_x*(br1+dpip_r)+br2*t1_x+br3*t2_x
1086 : jr(3,5)=n_y*(br1+dpip_r)+br2*t1_y+br3*t2_y
1087 : jr(4,5)=n_z*(br1+dpip_r)+br2*t1_z+br3*t2_z
1088 : c-------------------------------------------------------
1089 : c-------------------------------------------------------
1090 : hr_l=(eold_l+pold_l/rold_l)
1091 : hr_r=(eold_r+pold_r/rold_r)
1092 : c-------------------------------------------------
1093 : if(rum .ge. 0.0d0) then
1094 : br1=rumr_l*hr_l-rum*ga*pold_l/gm1/rold_l/rold_l
1095 : else
1096 : br1=rumr_l*hr_r
1097 : endif
1098 : jl(5,1)=br1
1099 : c-------------------------------------------------
1100 : if(rum .ge. 0.0d0) then
1101 : br1=rumu_l*hr_l+rum*uold_l
1102 : else
1103 : br1=rumu_l*hr_r
1104 : endif
1105 : jl(5,2)=br1
1106 : c-------------------------------------------------
1107 : if(rum .ge. 0.0d0) then
1108 : br1=rumv_l*hr_l+rum*vold_l
1109 : else
1110 : br1=rumv_l*hr_r
1111 : endif
1112 : jl(5,3)=br1
1113 : c-------------------------------------------------
1114 : if(rum .ge. 0.0d0) then
1115 : br1=rumw_l*hr_l+rum*wold_l
1116 : else
1117 : br1=rumw_l*hr_r
1118 : endif
1119 : jl(5,4)=br1
1120 : c-------------------------------------------------
1121 : if(rum .ge. 0.0d0) then
1122 : br1=rump_l*hr_l+rum*ga/gm1/rold_l
1123 : else
1124 : br1=rump_l*hr_r
1125 : endif
1126 : jl(5,5)=br1
1127 : c-------------------------------------------------
1128 : c-------------------------------------------------
1129 : if(rum .ge. 0.0d0) then
1130 : br1=rumr_r*hr_l
1131 : else
1132 : br1=rumr_r*hr_r-rum*ga*pold_r/gm1/rold_r/rold_r
1133 : endif
1134 : jr(5,1)=br1
1135 : c---------------------
1136 : if(rum .ge. 0.0d0) then
1137 : br1=rumu_r*hr_l
1138 : else
1139 : br1=rumu_r*hr_r+rum*uold_r
1140 : endif
1141 : jr(5,2)=br1
1142 : c---------------------
1143 : if(rum .ge. 0.0d0) then
1144 : br1=rumv_r*hr_l
1145 : else
1146 : br1=rumv_r*hr_r+rum*vold_r
1147 : endif
1148 : jr(5,3)=br1
1149 : c---------------------
1150 : if(rum .ge. 0.0d0) then
1151 : br1=rumw_r*hr_l
1152 : else
1153 : br1=rumw_r*hr_r+rum*wold_r
1154 : endif
1155 : jr(5,4)=br1
1156 : c---------------------
1157 : if(rum .ge. 0.0d0) then
1158 : br1=rump_r*hr_l
1159 : else
1160 : br1=rump_r*hr_r+rum*ga/gm1/rold_r
1161 : endif
1162 : jr(5,5)=br1
1163 : c---------------------------------------------------
1164 : c----------------------------------------------------------
1165 : c matrixes wwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwww
1166 : c----------------------------------------------------------
1167 : wl(1,1)=1.0d0
1168 : wl(1,2)=0.0d0
1169 : wl(1,3)=0.0d0
1170 : wl(1,4)=0.0d0
1171 : wl(1,5)=0.0d0
1172 : c------------------------------
1173 : wl(2,1)=-uold_l/rold_l
1174 : wl(2,2)=1.0d0/rold_l
1175 : wl(2,3)=0.0d0
1176 : wl(2,4)=0.0d0
1177 : wl(2,5)=0.0d0
1178 : c------------------------------
1179 : wl(3,1)=-vold_l/rold_l
1180 : wl(3,2)=0.0d0
1181 : wl(3,3)=1.0d0/rold_l
1182 : wl(3,4)=0.0d0
1183 : wl(3,5)=0.0d0
1184 : c------------------------------
1185 : wl(4,1)=-wold_l/rold_l
1186 : wl(4,2)=0.0d0
1187 : wl(4,3)=0.0d0
1188 : wl(4,4)=1.0d0/rold_l
1189 : wl(4,5)=0.0d0
1190 : c------------------------------
1191 : wl(5,1)=gm1*(uold_l*uold_l+vold_l*vold_l+wold_l*wold_l)/2.0d0
1192 : wl(5,2)=-uold_l*gm1
1193 : wl(5,3)=-vold_l*gm1
1194 : wl(5,4)=-wold_l*gm1
1195 : wl(5,5)=gm1
1196 : c------------------------------
1197 : c------------------------------
1198 : wr(1,1)=1.0d0
1199 : wr(1,2)=0.0d0
1200 : wr(1,3)=0.0d0
1201 : wr(1,4)=0.0d0
1202 : wr(1,5)=0.0d0
1203 : c------------------------------
1204 : wr(2,1)=-uold_r/rold_r
1205 : wr(2,2)=1.0d0/rold_r
1206 : wr(2,3)=0.0d0
1207 : wr(2,4)=0.0d0
1208 : wr(2,5)=0.0d0
1209 : c------------------------------
1210 : wr(3,1)=-vold_r/rold_r
1211 : wr(3,2)=0.0d0
1212 : wr(3,3)=1.0d0/rold_r
1213 : wr(3,4)=0.0d0
1214 : wr(3,5)=0.0d0
1215 : c------------------------------
1216 : wr(4,1)=-wold_r/rold_r
1217 : wr(4,2)=0.0d0
1218 : wr(4,3)=0.0d0
1219 : wr(4,4)=1.0d0/rold_r
1220 : wr(4,5)=0.0d0
1221 : c------------------------------
1222 : wr(5,1)=gm1*(uold_r*uold_r+vold_r*vold_r+wold_r*wold_r)/2.0d0
1223 : wr(5,2)=-uold_r*gm1
1224 : wr(5,3)=-vold_r*gm1
1225 : wr(5,4)=-wold_r*gm1
1226 : wr(5,5)=gm1
1227 : c----------------------------------------------
1228 : c----------------------------------------------
1229 : do 1 i=1,5
1230 : do 2 j=1,5
1231 : jtl(i,j)=0.0d0
1232 : jtr(i,j)=0.0d0
1233 : do 3 k=1,5
1234 : jtl(i,j)=jtl(i,j)+jl(i,k)*wl(k,j)
1235 : jtr(i,j)=jtr(i,j)+jr(i,k)*wr(k,j)
1236 : 3 continue
1237 : 2 continue
1238 : 1 continue
1239 : c----------------------------------------------
1240 : return
1241 : end
1242 :
1243 :
1244 :
1245 :
1246 :
1247 :
1248 :
1249 :
1250 :
1251 :
1252 :
1253 :
1254 :
1255 :
1256 :
1257 :
1258 :
1259 :
1260 :
1261 :
© Cast3M 2003 - All rights reserved.
Disclaimer