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