1 : C JAKBM SOURCE CHAT 05/01/13 00:48:48 5004
2 : SUBROUTINE JAKBM(JTL,JTR,WVEC_L,WVEC_R,NVECT,TVECT,
3 : & ga,v_inf)
4 : C************************************************************************
5 : C
6 : C PROJET : CASTEM 2000
7 : C
8 : C NOM : JACBM
9 : C
10 : C DESCRIPTION : Voir KONJA6
11 : C
12 : C LANGAGE : FORTRAN 77 + ESOPE 2000 (avec estensions CISI)
13 : C
14 : C AUTEUR : S. KUDRIAKOV, DM2S/SFME/LTMF
15 : C
16 : C************************************************************************
17 : C
18 : c----------------------------------------------------------------------
19 : c GENERAL DESCRIPTION:
20 : c This subroutine provides the jacobians which are the derivatives
21 : c of the numerical flux function defined at the cell interface
22 : c with respect to the conservative variables of the left and right
23 : c cells (relative to the cell interface).
24 : c The low-mach number corrections are made for the flux functions
25 : c
26 : c EQUATIONS: 2D Euler equations of gas dynamics
27 : c
28 : c
29 : c REFERENCE: 1) JCP, 129, 364-382 (1996)
30 : c " A Sequel to AUSM: AUSM+ ".
31 : c M.S.Liou
32 : c 2) AIAA Journal, Sept. 1998
33 : c " Low-Diffusion Flux-Splitting Methods for Flows at All Speeds"
34 : c J.R.Edwards and M.S.Liou
35 : c----------------------------------------------------------------------
36 : c INPUT:
37 : c
38 : c alpha -- parameter of the AUSM+ scheme in the Pressure function;
39 : c ( -3/4 <= alpha <= 3/16 ) (imposed as a parameter)
40 : c
41 : c beta -- parameter of the AUSM+ scheme in the Mach function;
42 : c ( -1/16 <= beta <= 1/2 ) (imposed as a parameter)
43 : c
44 : c wvec_l -- vector of the primitive variables (rho,ux,uy,p) at the
45 : c left cell;
46 : c
47 : c wvec_r -- vector of the primitive variables (rho,ux,uy,p) at the
48 : c right cell;
49 : c
50 : c nvect -- normal vector to the interface (2 components in 2D);
51 : c
52 : c tvect -- tangential vector to the interface;
53 : c
54 : c ga -- ratio of the specific heats (assumed constant)
55 : c
56 : c v_inf -- parameter for reference velocity, see Ref. 2).
57 : c----------------------------------------------------------------------
58 : c
59 : c OUTPUT:
60 : c
61 : c jtl -- jakobian matrix 4 by 4 - derivatives of the numerical
62 : c flux function with respect to the conservative variables
63 : c from the left cell;
64 : c
65 : c jtr -- jakobian matrix 4 by 4 - derivatives of the numerical
66 : c flux function with respect to the conservative variables
67 : c from the right cell.
68 : c----------------------------------------------------------------------
69 : IMPLICIT INTEGER(I-N)
70 : real*8 wvec_l(4),wvec_r(4)
71 : real*8 nvect(2),tvect(2)
72 : real*8 jl(4,4),jr(4,4)
73 : real*8 wl(4,4),wr(4,4)
74 : real*8 jtl(4,4),jtr(4,4)
75 : real*8 alpha,beta
76 : real*8 ga,gm1,temph
77 : real*8 n_x,n_y,t_x,t_y
78 : real*8 un_l, un_r, ut_l, ut_r
79 : real*8 ml,mr,Mplus,Mmin,mmid
80 : real*8 mpl_m, mmin_m,am
81 : real*8 rold_l,uold_l,vold_l,pold_l,eold_l
82 : real*8 rold_r,uold_r,vold_r,pold_r,eold_r
83 : real*8 Pplus,Pmin
84 : real*8 hr_l,hr_r,top,bot,bots
85 : real*8 br1,br3,br4,temp_l,temp_r,brac_l,brac_r
86 : real*8 aleft, arigh
87 : real*8 damr_l,damr_r,damu_l,damu_r
88 : real*8 damv_l,damv_r,damp_l,damp_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 dmrr_l,dmrr_r,dmru_l,dmru_r
92 : real*8 dmrv_l,dmrv_r,dmrp_l,dmrp_r
93 : real*8 dMpr_l,dMpr_r,dMpu_l,dMpu_r
94 : real*8 dMpv_l,dMpv_r,dMpp_l,dMpp_r
95 : real*8 dMmr_l,dMmr_r,dMmu_l,dMmu_r
96 : real*8 dMmv_l,dMmv_r,dMmp_l,dMmp_r
97 : real*8 dmir_l,dmir_r,dmiu_l,dmiu_r
98 : real*8 dmiv_l,dmiv_r,dmip_l,dmip_r
99 : real*8 d3mr_l,d3mr_r,d3mu_l,d3mu_r
100 : real*8 d3mv_l,d3mv_r,d3mp_l,d3mp_r
101 : real*8 d2mr_l,d2mr_r,d2mu_l,d2mu_r
102 : real*8 d2mv_l,d2mv_r,d2mp_l,d2mp_r
103 : real*8 dPpr_l,dPpr_r,dPpu_l,dPpu_r
104 : real*8 dPpv_l,dPpv_r,dPpp_l,dPpp_r
105 : real*8 dPmr_l,dPmr_r,dPmu_l,dPmu_r
106 : real*8 dPmv_l,dPmv_r,dPmp_l,dPmp_r
107 : real*8 dpir_l,dpir_r,dpiu_l,dpiu_r
108 : real*8 dpiv_l,dpiv_r,dpip_l,dpip_r
109 : real*8 epsil,qq,amw,Mmin1,Mplus1
110 : real*8 fmid,mlw,mrw,termp
111 : real*8 ur_r,ur_l,urm,mhalf,mhalfr
112 : real*8 durr_l,durr_r,duru_l,duru_r
113 : real*8 durv_l,durv_r,durp_l,durp_r
114 : real*8 dmhr_l,dmhr_r,dmhu_l,dmhu_r
115 : real*8 dmhv_l,dmhv_r,dmhp_l,dmhp_r
116 : real*8 dmfr_l,dmfr_r,dmfu_l,dmfu_r
117 : real*8 dmfv_l,dmfv_r,dmfp_l,dmfp_r
118 : real*8 dfm_mf,dfm_mh
119 : real*8 dfmr_l,dfmr_r,dfmu_l,dfmu_r
120 : real*8 dfmv_l,dfmv_r,dfmp_l,dfmp_r
121 : real*8 m1mr_l,m1mr_r,m1mu_l,m1mu_r
122 : real*8 m1mv_l,m1mv_r,m1mp_l,m1mp_r
123 : real*8 m1pr_l,m1pr_r,m1pu_l,m1pu_r
124 : real*8 m1pv_l,m1pv_r,m1pp_l,m1pp_r
125 : real*8 tmpr_l,tmpr_r,tmpu_l,tmpu_r
126 : real*8 tmpv_l,tmpv_r,tmpp_l,tmpp_r
127 : real*8 coef,canc,v_inf
128 : real*8 sr_l,sr_r,su_l,su_r,sv_l,sv_r,sp_l,sp_r
129 : real*8 rum,rumr_l,rumu_l,rumv_l,rump_l
130 : real*8 rumr_r,rumu_r,rumv_r,rump_r
131 : integer i,j,k
132 : parameter(alpha = 0.1875D0, beta = 0.125D0)
133 : parameter(epsil = 1.0d0)
134 : canc=1.0d0
135 : c-------------------------------------------------------------
136 : n_x=nvect(1)
137 : n_y=nvect(2)
138 : t_x=tvect(1)
139 : t_y=tvect(2)
140 : c----------------------------
141 : gm1=ga-1.0d0
142 : c----------------------------
143 : rold_l=wvec_l(1)
144 : uold_l=wvec_l(2)
145 : vold_l=wvec_l(3)
146 : pold_l=wvec_l(4)
147 : c-----------------------
148 : rold_r=wvec_r(1)
149 : uold_r=wvec_r(2)
150 : vold_r=wvec_r(3)
151 : pold_r=wvec_r(4)
152 : c------------------------------------------------------------------
153 : c Computation of the specific total energy on the left and right.
154 : c------------------------------------------------------------------
155 : eold_l=(uold_l*uold_l+vold_l*vold_l)/2.0d0
156 : eold_l=eold_l+pold_l/(gm1*rold_l)
157 : eold_r=(uold_r*uold_r+vold_r*vold_r)/2.0d0
158 : eold_r=eold_r+pold_r/(gm1*rold_r)
159 : c------------------------------------------------------------------
160 : c Computing reference velocity and its derivatives
161 : c------------------------------------------------------------------
162 : aleft=sqrt(ga*pold_l/rold_l)
163 : arigh=sqrt(ga*pold_r/rold_r)
164 : qq=sqrt(uold_l*uold_l+vold_l*vold_l)
165 : if(qq .lt. (epsil*v_inf)) then
166 : ur_l = epsil*v_inf
167 : durr_l=0.0d0
168 : duru_l=0.0d0
169 : durv_l=0.0d0
170 : durp_l=0.0d0
171 : else
172 : ur_l=qq
173 : durr_l=0.0d0
174 : duru_l=uold_l/qq
175 : durv_l=vold_l/qq
176 : durp_l=0.0d0
177 : endif
178 : c------------------------------
179 : if(ur_l .ge. aleft) then
180 : ur_l=aleft
181 : durr_l=-aleft/(2.0d0*rold_l)
182 : duru_l=0.0d0
183 : durv_l=0.0d0
184 : durp_l=ga/(2.0d0*aleft*rold_l)
185 : endif
186 : c------------------------------------------------------------------
187 : qq=sqrt(uold_r*uold_r+vold_r*vold_r)
188 : if(qq .lt. (epsil*v_inf)) then
189 : ur_r = epsil*v_inf
190 : durr_r=0.0d0
191 : duru_r=0.0d0
192 : durv_r=0.0d0
193 : durp_r=0.0d0
194 : else
195 : ur_r=qq
196 : durr_r=0.0d0
197 : duru_r=uold_r/qq
198 : durv_r=vold_r/qq
199 : durp_r=0.0d0
200 : endif
201 : c---------------------------
202 : if(ur_r .ge. arigh) then
203 : ur_r=arigh
204 : durr_r=-arigh/(2.0d0*rold_r)
205 : duru_r=0.0d0
206 : durv_r=0.0d0
207 : durp_r=ga/(2.0d0*arigh*rold_r)
208 : endif
209 : c------------------------------------------------------------------
210 : c-------------------------------------------------------------------
211 : c Reference velocity at the interface is taken as an average
212 : c of the reference velocities of the neighbouring cells
213 : c-------------------------------------------------------------------
214 : urm=0.5d0*(ur_l+ur_r)
215 : durr_l=0.5d0*durr_l
216 : duru_l=0.5d0*duru_l
217 : durv_l=0.5d0*durv_l
218 : durp_l=0.5d0*durp_l
219 : c-------------------------
220 : durr_r=0.5d0*durr_r
221 : duru_r=0.5d0*duru_r
222 : durv_r=0.5d0*durv_r
223 : durp_r=0.5d0*durp_r
224 : c-------------------------------------------------------------------
225 : c Computation of the speed of sound and its derivatives;
226 : c numerical speed of sound at the interface is taken as an average
227 : c of the speeds of sounds of the neighbouring cells
228 : c-------------------------------------------------------------------
229 : am=0.5d0*(aleft+arigh)
230 : c-------------------------------------------------------------------
231 : if(abs(urm/am-1.0d0).le.0.000001d0) then
232 : coef=0.0d0
233 : else
234 : coef=1.0d0
235 : endif
236 : c-------------------------------------------------------------------
237 : damr_r=-arigh/(4.0d0*rold_r)
238 : damu_r=0.0d0
239 : damv_r=0.0d0
240 : damp_r=ga/(4.0d0*arigh*rold_r)
241 : c-----------------------
242 : damr_l=-aleft/(4.0d0*rold_l)
243 : damu_l=0.0d0
244 : damv_l=0.0d0
245 : damp_l=ga/(4.0d0*aleft*rold_l)
246 : c-------------------------------------------------------------------
247 : c Computing numerical Mach number and its derivatives,
248 : c see p.370, under (A1).
249 : c-------------------------------------------------------------------
250 : un_l=uold_l*n_x+vold_l*n_y
251 : un_r=uold_r*n_x+vold_r*n_y
252 : ut_l=uold_l*t_x+vold_l*t_y
253 : ut_r=uold_r*t_x+vold_r*t_y
254 : c-------------------------------------------------------------------
255 : ml=un_l/am
256 : mr=un_r/am
257 : mhalf=0.5d0*(un_l+un_r)/am
258 : c---------------------
259 : top=0.5d0*(un_l+un_r)/(am*am)
260 : dmhr_l=-top*damr_l
261 : dmhu_l=n_x/2.0d0/am-top*damu_l
262 : dmhv_l=n_y/2.0d0/am-top*damv_l
263 : dmhp_l=-top*damp_l
264 : c---------------------
265 : dmhr_r=-top*damr_r
266 : dmhu_r=n_x/2.0d0/am-top*damu_r
267 : dmhv_r=n_y/2.0d0/am-top*damv_r
268 : dmhp_r=-top*damp_r
269 : c--------------------------------
270 : mhalfr=urm/am
271 : c---------------------
272 : top=urm/(am*am)
273 : dmfr_l=durr_l/am-top*damr_l
274 : dmfu_l=duru_l/am-top*damu_l
275 : dmfv_l=durv_l/am-top*damv_l
276 : dmfp_l=durp_l/am-top*damp_l
277 : c---------------------
278 : dmfr_r=durr_r/am-top*damr_r
279 : dmfu_r=duru_r/am-top*damu_r
280 : dmfv_r=durv_r/am-top*damv_r
281 : dmfp_r=durp_r/am-top*damp_r
282 : c-------------------------------------------------------------------
283 : c Scaling function for the speed of sound and its derivatives
284 : c-------------------------------------------------------------------
285 : top=(1.0d0-mhalfr*mhalfr)*(1.0d0-mhalfr*mhalfr)
286 : top=top*mhalf*mhalf+4.0d0*mhalfr*mhalfr
287 : bot=1.0d0+mhalfr*mhalfr
288 : if(abs(canc-0.0d0).lt.0.000001d0) then
289 : fmid=1.0d0
290 : else
291 : fmid=sqrt(top)/bot
292 : endif
293 : c--------------------------
294 : temph=-4.0d0*(1.0d0-mhalfr*mhalfr)*mhalfr
295 : temph=temph*mhalf*mhalf+8.0d0*mhalfr
296 : dfm_mf=temph/(sqrt(top)*2.0d0*bot)
297 : dfm_mf=dfm_mf-sqrt(top)*2.0d0*mhalfr/(bot*bot)
298 : c--------------------------
299 : temph=2.0d0*(1.0d0-mhalfr*mhalfr)*(1.0d0-mhalfr*mhalfr)*mhalf
300 : dfm_mh=temph/(2.0d0*bot*sqrt(top))
301 : c--------------------------
302 : dfmr_l=dfm_mf*dmfr_l+dfm_mh*dmhr_l
303 : dfmu_l=dfm_mf*dmfu_l+dfm_mh*dmhu_l
304 : dfmv_l=dfm_mf*dmfv_l+dfm_mh*dmhv_l
305 : dfmp_l=dfm_mf*dmfp_l+dfm_mh*dmhp_l
306 : c--------------------------
307 : dfmr_r=dfm_mf*dmfr_r+dfm_mh*dmhr_r
308 : dfmu_r=dfm_mf*dmfu_r+dfm_mh*dmhu_r
309 : dfmv_r=dfm_mf*dmfv_r+dfm_mh*dmhv_r
310 : dfmp_r=dfm_mf*dmfp_r+dfm_mh*dmhp_r
311 : c--------------------------
312 : amw=fmid*am
313 : mlw=un_l/amw
314 : mrw=un_r/amw
315 : c--------------------------
316 : damr_l=canc*coef*dfmr_l*am+fmid*damr_l
317 : damu_l=canc*coef*dfmu_l*am+fmid*damu_l
318 : damv_l=canc*coef*dfmv_l*am+fmid*damv_l
319 : damp_l=canc*coef*dfmp_l*am+fmid*damp_l
320 : c--------------------------
321 : damr_r=canc*coef*dfmr_r*am+fmid*damr_r
322 : damu_r=canc*coef*dfmu_r*am+fmid*damu_r
323 : damv_r=canc*coef*dfmv_r*am+fmid*damv_r
324 : damp_r=canc*coef*dfmp_r*am+fmid*damp_r
325 : c-/-/-/-/-/-/-/-/-/-/-/-/-/-/-/-/-/
326 : am=amw
327 : c-------------------------------------------------------------------
328 : c Redefinition of the numerical mach numbers
329 : c-------------------------------------------------------------------
330 : if(abs(canc-0.0d0).lt.0.000001d0) then
331 : top=2.0d0
332 : bot=0.0d0
333 : else
334 : top=1.0d0+mhalfr*mhalfr
335 : bot=1.0d0-mhalfr*mhalfr
336 : endif
337 : ml=0.5d0*(top*mlw+bot*mrw)
338 : mr=0.5d0*(top*mrw+bot*mlw)
339 : c-------------------------------------------------------------------
339 : c Mplus and Mmin are calligraphic lettes M+ and M- from the paper,
340 : c see (19a) and (19b), p.367.
341 : c-------------------------------------------------------------------
342 : if(abs(ml) .ge. 1.0d0) then
343 : Mplus=(ml+abs(ml))/2.0d0
344 : else
345 : Mplus=(ml+1.0d0)*(ml+1.0d0)/4.0d0
346 : Mplus=Mplus+beta*(ml*ml-1.0d0)*(ml*ml-1.0d0)
347 : endif
348 : Mplus1=(ml+abs(ml))/2.0d0
349 : c-------------------------------------------------------------------
350 : if(abs(mr) .ge. 1.0d0) then
351 : Mmin=(mr-abs(mr))/2.0d0
352 : else
353 : Mmin=-(mr-1.0d0)*(mr-1.0d0)/4.0d0
354 : Mmin=Mmin-beta*(mr*mr-1.0d0)*(mr*mr-1.0d0)
355 : endif
356 : Mmin1=(mr-abs(mr))/2.0d0
357 : c-------------------------------------------------------------------
358 : c Derivatives of ml and mr with respect to both sets of primitive
359 : c variables: left and right.
360 : c-------------------------------------------------------------------
361 : temp_l=-un_l/(am*am)
362 : temp_r=-un_r/(am*am)
363 : c--------
364 : dmlr_l=temp_l*damr_l
365 : dmlr_r=temp_l*damr_r
366 : dmrr_l=temp_r*damr_l
367 : dmrr_r=temp_r*damr_r
368 : c--------
369 : dmlu_l=n_x/am+temp_l*damu_l
370 : dmlu_r=temp_l*damu_r
371 : dmru_l=temp_r*damu_l
372 : dmru_r=n_x/am+temp_r*damu_r
373 : c--------
374 : dmlv_l=n_y/am+temp_l*damv_l
375 : dmlv_r=temp_l*damv_r
376 : dmrv_l=temp_r*damv_l
377 : dmrv_r=n_y/am+temp_r*damv_r
378 : c--------
379 : dmlp_l=temp_l*damp_l
380 : dmlp_r=temp_l*damp_r
381 : dmrp_l=temp_r*damp_l
382 : dmrp_r=temp_r*damp_r
383 : c-----------------------------
384 : sr_l=dmlr_l
385 : sr_r=dmlr_r
386 : su_l=dmlu_l
387 : su_r=dmlu_r
388 : sv_l=dmlv_l
389 : sv_r=dmlv_r
390 : sp_l=dmlp_l
391 : sp_r=dmlp_r
392 : c-----------------------------------------------------------------
393 : c Redefinition of the derivatives of the ml & mr
394 : c-----------------------------------------------------------------
395 : temp_l=(mlw-mrw)*mhalfr
396 : temp_r=-temp_l
397 : c--------
398 : dmlr_l=0.5d0*(top*dmlr_l+bot*dmrr_l)+canc*coef*temp_l*dmfr_l
399 : dmlu_l=0.5d0*(top*dmlu_l+bot*dmru_l)+canc*coef*temp_l*dmfu_l
400 : dmlv_l=0.5d0*(top*dmlv_l+bot*dmrv_l)+canc*coef*temp_l*dmfv_l
401 : dmlp_l=0.5d0*(top*dmlp_l+bot*dmrp_l)+canc*coef*temp_l*dmfp_l
402 : c--------
403 : dmlr_r=0.5d0*(top*dmlr_r+bot*dmrr_r)+canc*coef*temp_l*dmfr_r
404 : dmlu_r=0.5d0*(top*dmlu_r+bot*dmru_r)+canc*coef*temp_l*dmfu_r
405 : dmlv_r=0.5d0*(top*dmlv_r+bot*dmrv_r)+canc*coef*temp_l*dmfv_r
406 : dmlp_r=0.5d0*(top*dmlp_r+bot*dmrp_r)+canc*coef*temp_l*dmfp_r
407 : c--------
408 : dmrr_l=0.5d0*(top*dmrr_l+bot*sr_l)+canc*coef*temp_r*dmfr_l
409 : dmru_l=0.5d0*(top*dmru_l+bot*su_l)+canc*coef*temp_r*dmfu_l
410 : dmrv_l=0.5d0*(top*dmrv_l+bot*sv_l)+canc*coef*temp_r*dmfv_l
411 : dmrp_l=0.5d0*(top*dmrp_l+bot*sp_l)+canc*coef*temp_r*dmfp_l
412 : c--------
413 : dmrr_r=0.5d0*(top*dmrr_r+bot*sr_r)+canc*coef*temp_r*dmfr_r
414 : dmru_r=0.5d0*(top*dmru_r+bot*su_r)+canc*coef*temp_r*dmfu_r
415 : dmrv_r=0.5d0*(top*dmrv_r+bot*sv_r)+canc*coef*temp_r*dmfv_r
416 : dmrp_r=0.5d0*(top*dmrp_r+bot*sp_r)+canc*coef*temp_r*dmfp_r
417 : c-----------------------------------------------------------
418 : c mmid is m_{1/2} (notation as in the paper, see (13),p.366)
419 : c-----------------------------------------------------------
420 : mmid=Mplus+Mmin
421 : c-----------------------------------------------------------
422 : c Computing the derivatives of M+ and M-
423 : c-----------------------------------------------------------
424 : if(ml .ge. 1.0d0) then
425 : dMpr_l=dmlr_l
426 : dMpu_l=dmlu_l
427 : dMpv_l=dmlv_l
428 : dMpp_l=dmlp_l
429 : c--------------------
430 : dMpr_r=dmlr_r
431 : dMpu_r=dmlu_r
432 : dMpv_r=dmlv_r
433 : dMpp_r=dmlp_r
434 : else
435 : if((ml .gt. -1.0d0) .and. (ml .lt. 1.0d0)) then
436 : temph=(ml+1.0d0)/2.0d0
437 : dMpr_l=(temph+4.0d0*beta*ml*(ml*ml-1.0d0))*dmlr_l
438 : dMpu_l=(temph+4.0d0*beta*ml*(ml*ml-1.0d0))*dmlu_l
439 : dMpv_l=(temph+4.0d0*beta*ml*(ml*ml-1.0d0))*dmlv_l
440 : dMpp_l=(temph+4.0d0*beta*ml*(ml*ml-1.0d0))*dmlp_l
441 : c--------------------
442 : dMpr_r=(temph+4.0d0*beta*ml*(ml*ml-1.0d0))*dmlr_r
443 : dMpu_r=(temph+4.0d0*beta*ml*(ml*ml-1.0d0))*dmlu_r
444 : dMpv_r=(temph+4.0d0*beta*ml*(ml*ml-1.0d0))*dmlv_r
445 : dMpp_r=(temph+4.0d0*beta*ml*(ml*ml-1.0d0))*dmlp_r
446 : else
447 : dMpr_l=0.0d0
448 : dMpu_l=0.0d0
449 : dMpv_l=0.0d0
450 : dMpp_l=0.0d0
451 : c---------------------
452 : dMpr_r=0.0d0
453 : dMpu_r=0.0d0
454 : dMpv_r=0.0d0
455 : dMpp_r=0.0d0
456 : endif
457 : endif
458 : c-----------------------------------------------------------
459 : c addition of low Mach number
460 : c-----------------------------------------------------------
461 : if(ml .ge. 0.0d0) then
462 : m1pr_l=dmlr_l
463 : m1pu_l=dmlu_l
464 : m1pv_l=dmlv_l
465 : m1pp_l=dmlp_l
466 : c---------------------
467 : m1pr_r=dmlr_r
468 : m1pu_r=dmlu_r
469 : m1pv_r=dmlv_r
470 : m1pp_r=dmlp_r
471 : else
472 : m1pr_l=0.0d0
473 : m1pu_l=0.0d0
474 : m1pv_l=0.0d0
475 : m1pp_l=0.0d0
476 : c--------------------
477 : m1pr_r=0.0d0
478 : m1pu_r=0.0d0
479 : m1pv_r=0.0d0
480 : m1pp_r=0.0d0
481 : endif
482 : c-----------------------------------------------------------
483 : if(mr .ge. 1.0d0) then
484 : dMmr_l=0.0d0
485 : dMmu_l=0.0d0
486 : dMmv_l=0.0d0
487 : dMmp_l=0.0d0
488 : c---------------------
489 : dMmr_r=0.0d0
490 : dMmu_r=0.0d0
491 : dMmv_r=0.0d0
492 : dMmp_r=0.0d0
493 : else
494 : if((mr .gt. -1.0d0) .and. (mr .lt. 1.0d0)) then
495 : temph=(-mr+1.0d0)/2.0d0
496 : dMmr_l=(temph-4.0d0*beta*mr*(mr*mr-1.0d0))*dmrr_l
497 : dMmu_l=(temph-4.0d0*beta*mr*(mr*mr-1.0d0))*dmru_l
498 : dMmv_l=(temph-4.0d0*beta*mr*(mr*mr-1.0d0))*dmrv_l
499 : dMmp_l=(temph-4.0d0*beta*mr*(mr*mr-1.0d0))*dmrp_l
500 : c--------------------
501 : dMmr_r=(temph-4.0d0*beta*mr*(mr*mr-1.0d0))*dmrr_r
502 : dMmu_r=(temph-4.0d0*beta*mr*(mr*mr-1.0d0))*dmru_r
503 : dMmv_r=(temph-4.0d0*beta*mr*(mr*mr-1.0d0))*dmrv_r
504 : dMmp_r=(temph-4.0d0*beta*mr*(mr*mr-1.0d0))*dmrp_r
505 : else
506 : dMmr_l=dmrr_l
507 : dMmu_l=dmru_l
508 : dMmv_l=dmrv_l
509 : dMmp_l=dmrp_l
510 : c--------------------
511 : dMmr_r=dmrr_r
512 : dMmu_r=dmru_r
513 : dMmv_r=dmrv_r
514 : dMmp_r=dmrp_r
515 : endif
516 : endif
517 : c-----------------------------------------------------------
518 : c addition of low Mach number
519 : c-----------------------------------------------------------
520 : if(mr .le. 0.0d0) then
521 : m1mr_l=dmrr_l
522 : m1mu_l=dmru_l
523 : m1mv_l=dmrv_l
524 : m1mp_l=dmrp_l
525 : c---------------------
526 : m1mr_r=dmrr_r
527 : m1mu_r=dmru_r
528 : m1mv_r=dmrv_r
529 : m1mp_r=dmrp_r
530 : else
531 : m1mr_l=0.0d0
532 : m1mu_l=0.0d0
533 : m1mv_l=0.0d0
534 : m1mp_l=0.0d0
535 : c--------------------
536 : m1mr_r=0.0d0
537 : m1mu_r=0.0d0
538 : m1mv_r=0.0d0
539 : m1mp_r=0.0d0
540 : endif
541 : c-----------------------------------------------------------------
542 : c computing the derivatives of m_{1/2} (notation as in the paper)
543 : c-----------------------------------------------------------------
544 : dmir_l=dMpr_l+dMmr_l
545 : dmir_r=dMpr_r+dMmr_r
546 : c-------------
547 : dmiu_l=dMpu_l+dMmu_l
548 : dmiu_r=dMpu_r+dMmu_r
549 : c-------------
550 : dmiv_l=dMpv_l+dMmv_l
551 : dmiv_r=dMpv_r+dMmv_r
552 : c-------------
553 : dmip_l=dMpp_l+dMmp_l
554 : dmip_r=dMpp_r+dMmp_r
555 : c----------------------------------------------------------------
556 : c computing the main convective variables and their derivatives
557 : c mpl_m is m^{+}_{1/2} (paper's notation) and
558 : c mmin_m is m^{-}_{1/2} (paper's notation), see (A2) on p.370.
559 : c----------------------------------------------------------------
560 : termp=(Mmin1-Mmin+Mplus-Mplus1)*(1.0d0/(mhalfr*mhalfr)-1.0d0)
561 : termp=termp*(pold_l-pold_r)/(pold_l/rold_l+pold_r/rold_r)
562 : c-------------------------------------------------------------
563 : c derivatives of the termp
564 : c-------------------------------------------------------------
565 : top=(Mmin1-Mmin+Mplus-Mplus1)
566 : bots=1.0d0/(pold_l/rold_l+pold_r/rold_r)
567 : bot=(pold_l-pold_r)*bots
568 : temph=1.0d0/(mhalfr*mhalfr)-1.0d0
569 : c---------------------------
570 : tmpr_l=(m1mr_l-dMmr_l+dMpr_l-m1pr_l)*bot*temph
571 : tmpr_l=tmpr_l-2.0d0*bot*top*dmfr_l/(mhalfr*mhalfr*mhalfr)
572 : tmpr_l=tmpr_l+temph*top*bot*bots*(pold_l/rold_l/rold_l)
573 : c---------------------------
574 : tmpu_l=(m1mu_l-dMmu_l+dMpu_l-m1pu_l)*bot*temph
575 : tmpu_l=tmpu_l-2.0d0*bot*top*dmfu_l/(mhalfr*mhalfr*mhalfr)
576 : c---------------------------
577 : tmpv_l=(m1mv_l-dMmv_l+dMpv_l-m1pv_l)*bot*temph
578 : tmpv_l=tmpv_l-2.0d0*bot*top*dmfv_l/(mhalfr*mhalfr*mhalfr)
579 : c---------------------------
580 : tmpp_l=(m1mp_l-dMmp_l+dMpp_l-m1pp_l)*bot*temph
581 : tmpp_l=tmpp_l-2.0d0*bot*top*dmfp_l/(mhalfr*mhalfr*mhalfr)
582 : tmpp_l=tmpp_l+temph*top*bots*(1.0d0-bot/rold_l)
583 : c------------rrrrrrrr-------
584 : c------------rrrrrrrr-------
585 : tmpr_r=(m1mr_r-dMmr_r+dMpr_r-m1pr_r)*bot*temph
586 : tmpr_r=tmpr_r-2.0d0*bot*top*dmfr_r/(mhalfr*mhalfr*mhalfr)
587 : tmpr_r=tmpr_r+temph*top*bot*bots*(pold_r/rold_r/rold_r)
588 : c---------------------------
589 : tmpu_r=(m1mu_r-dMmu_r+dMpu_r-m1pu_r)*bot*temph
590 : tmpu_r=tmpu_r-2.0d0*bot*top*dmfu_r/(mhalfr*mhalfr*mhalfr)
591 : c---------------------------
592 : tmpv_r=(m1mv_r-dMmv_r+dMpv_r-m1pv_r)*bot*temph
593 : tmpv_r=tmpv_r-2.0d0*bot*top*dmfv_r/(mhalfr*mhalfr*mhalfr)
594 : c---------------------------
595 : tmpp_r=(m1mp_r-dMmp_r+dMpp_r-m1pp_r)*bot*temph
596 : tmpp_r=tmpp_r-2.0d0*bot*top*dmfp_r/(mhalfr*mhalfr*mhalfr)
597 : tmpp_r=tmpp_r-temph*top*bots*(1.0d0+bot/rold_r)
598 : c-------------------------------------------------------------
599 : if(mmid .ge. 0.0d0) then
600 : mpl_m = mmid
601 : d2mr_l=dmir_l
602 : d2mu_l=dmiu_l
603 : d2mv_l=dmiv_l
604 : d2mp_l=dmip_l
605 : c------------
606 : d2mr_r=dmir_r
607 : d2mu_r=dmiu_r
608 : d2mv_r=dmiv_r
609 : d2mp_r=dmip_r
610 : c------------
611 : else
612 : mpl_m = 0.0d0
613 : d2mr_l=0.0d0
614 : d2mu_l=0.0d0
615 : d2mv_l=0.0d0
616 : d2mp_l=0.0d0
617 : c------------
618 : d2mr_r=0.0d0
619 : d2mu_r=0.0d0
620 : d2mv_r=0.0d0
621 : d2mp_r=0.0d0
622 : endif
623 : c---------------------------------------------
624 : cc derivatives for the term termm
625 : cc------------------------------------------------------------------
626 : if(mmid .le. 0.0d0) then
627 : mmin_m = mmid
628 : d3mr_l=dmir_l
629 : d3mu_l=dmiu_l
630 : d3mv_l=dmiv_l
631 : d3mp_l=dmip_l
632 : c------------
633 : d3mr_r=dmir_r
634 : d3mu_r=dmiu_r
635 : d3mv_r=dmiv_r
636 : d3mp_r=dmip_r
637 : c------------
638 : else
639 : mmin_m = 0.0d0
640 : d3mr_l=0.0d0
641 : d3mu_l=0.0d0
642 : d3mv_l=0.0d0
643 : d3mp_l=0.0d0
644 : c------------
645 : d3mr_r=0.0d0
646 : d3mu_r=0.0d0
647 : d3mv_r=0.0d0
648 : d3mp_r=0.0d0
649 : endif
650 : c---------------------------------------------------------------
650 : c Computing the calligraphic P+ and P- with their derivatives,
651 : c see (21a) & (21b) on p.368.
652 : c---------------------------------------------------------------
653 : if(ml .ge. 1.0d0) then
654 : Pplus = 1.0d0
655 : else
656 : if((ml .gt. -1.0d0) .and. (ml .lt. 1.0d0)) then
657 : Pplus=(ml+1.0d0)*(ml+1.0d0)*(2.0d0-ml)/4.0d0
658 : Pplus=Pplus+alpha*ml*(ml*ml-1.0d0)*(ml*ml-1.0d0)
659 : else
660 : Pplus = 0.0d0
661 : endif
662 : endif
663 : c---------------------------------------------------------------
664 : if(mr .ge. 1.0d0) then
665 : Pmin = 0.0d0
666 : else
667 : if((mr .gt. -1.0d0) .and. (mr .lt. 1.0d0)) then
668 : Pmin=(mr-1.0d0)*(mr-1.0d0)*(2.0d0+mr)/4.0d0
669 : Pmin=Pmin-alpha*mr*(mr*mr-1.0d0)*(mr*mr-1.0d0)
670 : else
671 : Pmin = 1.0d0
672 : endif
673 : endif
674 : c---------------------------------------------------------------
675 : brac_l=(ml+1.0d0)*(2.0d0-ml)/2.0d0-(ml+1.0d0)*(ml+1.0d0)/4.0d0
676 : brac_l=brac_l+alpha*(ml*ml-1.0d0)*(ml*ml-1.0d0)
677 : brac_l=brac_l+4.0d0*alpha*ml*ml*(ml*ml-1.0d0)
678 : c--------------
679 : brac_r=(mr-1.0d0)*(2.0d0+mr)/2.0d0+(mr-1.0d0)*(mr-1.0d0)/4.0d0
680 : brac_r=brac_r-alpha*(mr*mr-1.0d0)*(mr*mr-1.0d0)
681 : brac_r=brac_r-4.0d0*alpha*mr*mr*(mr*mr-1.0d0)
682 : c---------------------------------------------------------------
683 : if((ml .gt. -1.0d0) .and. (ml .lt. 1.0d0)) then
684 : dPpr_l=brac_l*dmlr_l
685 : dPpr_r=brac_l*dmlr_r
686 : c------------
687 : dPpu_l=brac_l*dmlu_l
688 : dPpu_r=brac_l*dmlu_r
689 : c------------
690 : dPpv_l=brac_l*dmlv_l
691 : dPpv_r=brac_l*dmlv_r
692 : c------------
693 : dPpp_l=brac_l*dmlp_l
694 : dPpp_r=brac_l*dmlp_r
695 : c------------
696 : else
697 : dPpr_l=0.0d0
698 : dPpr_r=0.0d0
699 : c-----------
700 : dPpu_l=0.0d0
701 : dPpu_r=0.0d0
702 : c-----------
703 : dPpv_l=0.0d0
704 : dPpv_r=0.0d0
705 : c-----------
706 : dPpp_l=0.0d0
707 : dPpp_r=0.0d0
708 : c-----------
709 : endif
710 : c---------------------------------------------------------------
711 : if((mr .gt. -1.0d0) .and. (mr .lt. 1.0d0)) then
712 : dPmr_l=brac_r*dmrr_l
713 : dPmr_r=brac_r*dmrr_r
714 : c------------
715 : dPmu_l=brac_r*dmru_l
716 : dPmu_r=brac_r*dmru_r
717 : c------------
718 : dPmv_l=brac_r*dmrv_l
719 : dPmv_r=brac_r*dmrv_r
720 : c------------
721 : dPmp_l=brac_r*dmrp_l
722 : dPmp_r=brac_r*dmrp_r
723 : c------------
724 : else
725 : dPmr_l=0.0d0
726 : dPmr_r=0.0d0
727 : c-----------
728 : dPmu_l=0.0d0
729 : dPmu_r=0.0d0
730 : c-----------
731 : dPmv_l=0.0d0
732 : dPmv_r=0.0d0
733 : c-----------
734 : dPmp_l=0.0d0
735 : dPmp_r=0.0d0
736 : c-----------
737 : endif
738 : c-------------------------------------------------------------------
739 : c computing pmid -- p_{1/2} and its derivatives, see (20b), p.367.
740 : c-------------------------------------------------------------------
741 : dpir_l=dPpr_l*pold_l+dPmr_l*pold_r
742 : dpiu_l=dPpu_l*pold_l+dPmu_l*pold_r
743 : dpiv_l=dPpv_l*pold_l+dPmv_l*pold_r
744 : dpip_l=dPpp_l*pold_l+Pplus+dPmp_l*pold_r
745 : c----------------------------
746 : dpir_r=dPpr_r*pold_l+dPmr_r*pold_r
747 : dpiu_r=dPpu_r*pold_l+dPmu_r*pold_r
748 : dpiv_r=dPpv_r*pold_l+dPmv_r*pold_r
749 : dpip_r=dPpp_r*pold_l+Pmin+dPmp_r*pold_r
750 : c---------------------------------------------------------------------
751 : c Computation of the mass flux (rho * u)_1/2
752 : c---------------------------------------------------------------
753 : rum=am*(mpl_m*rold_l+mmin_m*rold_r)+canc*am*termp
754 : c-------------------------------------------------------
755 : rumr_l=damr_l*(mpl_m*rold_l+mmin_m*rold_r)+
756 : & am*(d2mr_l*rold_l+mpl_m+d3mr_l*rold_r)
757 : rumr_l=rumr_l+canc*coef*(damr_l*termp+am*tmpr_l)
758 : rumu_l=damu_l*(mpl_m*rold_l+mmin_m*rold_r)+
759 : & am*(d2mu_l*rold_l+d3mu_l*rold_r)
760 : rumu_l=rumu_l+canc*coef*(damu_l*termp+am*tmpu_l)
761 : rumv_l=damv_l*(mpl_m*rold_l+mmin_m*rold_r)+
762 : & am*(d2mv_l*rold_l+d3mv_l*rold_r)
763 : rumv_l=rumv_l+canc*coef*(damv_l*termp+am*tmpv_l)
764 : rump_l=damp_l*(mpl_m*rold_l+mmin_m*rold_r)+
765 : & am*(d2mp_l*rold_l+d3mp_l*rold_r)
766 : rump_l=rump_l+canc*coef*(damp_l*termp+am*tmpp_l)
767 : c-------------------------------------------------
768 : rumr_r=damr_r*(mpl_m*rold_l+mmin_m*rold_r)+
769 : & am*(d2mr_r*rold_l+mmin_m+d3mr_r*rold_r)
770 : rumr_r=rumr_r+canc*coef*(damr_r*termp+am*tmpr_r)
771 : rumu_r=damu_r*(mpl_m*rold_l+mmin_m*rold_r)+
772 : & am*(d2mu_r*rold_l+d3mu_r*rold_r)
773 : rumu_r=rumu_r+canc*coef*(damu_r*termp+am*tmpu_r)
774 : rumv_r=damv_r*(mpl_m*rold_l+mmin_m*rold_r)+
775 : & am*(d2mv_r*rold_l+d3mv_r*rold_r)
776 : rumv_r=rumv_r+canc*coef*(damv_r*termp+am*tmpv_r)
777 : rump_r=damp_r*(mpl_m*rold_l+mmin_m*rold_r)+
778 : & am*(d2mp_r*rold_l+d3mp_r*rold_r)
779 : rump_r=rump_r+canc*coef*(damp_r*termp+am*tmpp_r)
780 : c---------------------------------------------------------------------
781 : c computing JACOBIAN as a derivative of the numerical flux function
782 : c with respect to the primitive variables.
783 : c Notation: jl(i,j) --- is the derivative of the i-component of the
784 : c flux function with respect to the j-component of the
785 : c vector of primitive variables of the left state.
786 : c jr(i,j) --- is the derivative of the i-component of the
787 : c flux function with respect to the j-component of the
788 : c vector of primitive variables of the right state.
789 : c---------------------------------------------------------------------
790 : jl(1,1)=rumr_l
791 : jl(1,2)=rumu_l
792 : jl(1,3)=rumv_l
793 : jl(1,4)=rump_l
794 : c------------------------------------
795 : jr(1,1)=rumr_r
796 : jr(1,2)=rumu_r
797 : jr(1,3)=rumv_r
798 : jr(1,4)=rump_r
799 : c------------------------------------
800 : c---------------------------------------------------------
801 : if(rum .ge. 0.0d0) then
802 : br3=rumr_l*un_l
803 : br4=rumr_l*ut_l
804 : else
805 : br3=rumr_l*un_r
806 : br4=rumr_l*ut_r
807 : endif
808 : jl(2,1)=n_x*(br3+dpir_l)+br4*t_x
809 : jl(3,1)=n_y*(br3+dpir_l)+br4*t_y
810 : c-------------------
811 : if(rum .ge. 0.0d0) then
812 : br3=rumu_l*un_l+rum*n_x
813 : br4=rumu_l*ut_l+rum*t_x
814 : else
815 : br3=rumu_l*un_r
816 : br4=rumu_l*ut_r
817 : endif
818 : jl(2,2)=n_x*(br3+dpiu_l)+br4*t_x
819 : jl(3,2)=n_y*(br3+dpiu_l)+br4*t_y
820 : c-------------------
821 : if(rum .ge. 0.0d0) then
822 : br3=rumv_l*un_l+rum*n_y
823 : br4=rumv_l*ut_l+rum*t_y
824 : else
825 : br3=rumv_l*un_r
826 : br4=rumv_l*ut_r
827 : endif
828 : jl(2,3)=n_x*(br3+dpiv_l)+br4*t_x
829 : jl(3,3)=n_y*(br3+dpiv_l)+br4*t_y
830 : c-------------------
831 : if(rum .ge. 0.0d0) then
832 : br3=rump_l*un_l
833 : br4=rump_l*ut_l
834 : else
835 : br3=rump_l*un_r
836 : br4=rump_l*ut_r
837 : endif
838 : jl(2,4)=n_x*(br3+dpip_l)+br4*t_x
839 : jl(3,4)=n_y*(br3+dpip_l)+br4*t_y
840 : c-------------------------------------------------------------
841 : c-------------------------------------------------------------
842 : if(rum .ge. 0.0d0) then
843 : br3=rumr_r*un_l
844 : br4=rumr_r*ut_l
845 : else
846 : br3=rumr_r*un_r
847 : br4=rumr_r*ut_r
848 : endif
849 : jr(2,1)=n_x*(br3+dpir_r)+br4*t_x
850 : jr(3,1)=n_y*(br3+dpir_r)+br4*t_y
851 : c-------------------
852 : if(rum .ge. 0.0d0) then
853 : br3=rumu_r*un_l
854 : br4=rumu_r*ut_l
855 : else
856 : br3=rumu_r*un_r+rum*n_x
857 : br4=rumu_r*ut_r+rum*t_x
858 : endif
859 : jr(2,2)=n_x*(br3+dpiu_r)+br4*t_x
860 : jr(3,2)=n_y*(br3+dpiu_r)+br4*t_y
861 : c-------------------
862 : if(rum .ge. 0.0d0) then
863 : br3=rumv_r*un_l
864 : br4=rumv_r*ut_l
865 : else
866 : br3=rumv_r*un_r+rum*n_y
867 : br4=rumv_r*ut_r+rum*t_y
868 : endif
869 : jr(2,3)=n_x*(br3+dpiv_r)+br4*t_x
870 : jr(3,3)=n_y*(br3+dpiv_r)+br4*t_y
871 : c-------------------
872 : if(rum .ge. 0.0d0) then
873 : br3=rump_r*un_l
874 : br4=rump_r*ut_l
875 : else
876 : br3=rump_r*un_r
877 : br4=rump_r*ut_r
878 : endif
879 : jr(2,4)=n_x*(br3+dpip_r)+br4*t_x
880 : jr(3,4)=n_y*(br3+dpip_r)+br4*t_y
881 : c-------------------------------------------------------------
882 : c ------ f44444444444444444444444444444 ---------
883 : c-------------------------------------------------------------
884 : hr_l=(uold_l*uold_l+vold_l*vold_l)/2.0d0+ga*pold_l/gm1/rold_l
885 : hr_r=(uold_r*uold_r+vold_r*vold_r)/2.0d0+ga*pold_r/gm1/rold_r
886 : c-------------------------------------------------------------
887 : if(rum .ge. 0.0d0) then
888 : br1=rumr_l*hr_l-rum*ga*pold_l/gm1/rold_l/rold_l
889 : else
890 : br1=rumr_l*hr_r
891 : endif
892 : jl(4,1)=br1
893 : c---------------------
894 : if(rum .ge. 0.0d0) then
895 : br1=rumu_l*hr_l+rum*uold_l
896 : else
897 : br1=rumu_l*hr_r
898 : endif
899 : jl(4,2)=br1
900 : c---------------------
901 : if(rum .ge. 0.0d0) then
902 : br1=rumv_l*hr_l+rum*vold_l
903 : else
904 : br1=rumv_l*hr_r
905 : endif
906 : jl(4,3)=br1
907 : c---------------------
908 : if(rum .ge. 0.0d0) then
909 : br1=rump_l*hr_l+rum*ga/gm1/rold_l
910 : else
911 : br1=rump_l*hr_r
912 : endif
913 : jl(4,4)=br1
914 : c----------------------------------------------------------
915 : if(rum .ge. 0.0d0) then
916 : br1=rumr_r*hr_l
917 : else
918 : br1=rumr_r*hr_r-rum*ga*pold_r/gm1/rold_r/rold_r
919 : endif
920 : jr(4,1)=br1
921 : c---------------------
922 : if(rum .ge. 0.0d0) then
923 : br1=rumu_r*hr_l
924 : else
925 : br1=rumu_r*hr_r+rum*uold_r
926 : endif
927 : jr(4,2)=br1
928 : c---------------------
929 : if(rum .ge. 0.0d0) then
930 : br1=rumv_r*hr_l
931 : else
932 : br1=rumv_r*hr_r+rum*vold_r
933 : endif
934 : jr(4,3)=br1
935 : c---------------------
936 : if(rum .ge. 0.0d0) then
937 : br1=rump_r*hr_l
938 : else
939 : br1=rump_r*hr_r+rum*ga/gm1/rold_r
940 : endif
941 : jr(4,4)=br1
942 : c-------------------------------------------------------------
943 : c matrix wl(i,j) represents the derivative of the i-component
944 : c of the vector of primitive variables of the left state with
945 : c respect to the j-component of the vector of the conservative
946 : c variables of the left state.
947 : c
948 : c Here: (rho, u, v, p) - vector of primitive variables;
949 : c (rho, rho u, rho v, rho e) -- conservative variables.
950 : c-------------------------------------------------------------
951 : wl(1,1)=1.0d0
952 : wl(1,2)=0.0d0
953 : wl(1,3)=0.0d0
954 : wl(1,4)=0.0d0
955 : c------------------------------
956 : wl(2,1)=-uold_l/rold_l
957 : wl(2,2)=1.0d0/rold_l
958 : wl(2,3)=0.0d0
959 : wl(2,4)=0.0d0
960 : c------------------------------
961 : wl(3,1)=-vold_l/rold_l
962 : wl(3,2)=0.0d0
963 : wl(3,3)=1.0d0/rold_l
964 : wl(3,4)=0.0d0
965 : c------------------------------
966 : wl(4,1)=gm1*(uold_l*uold_l+vold_l*vold_l)/2.0d0
967 : wl(4,2)=-uold_l*gm1
968 : wl(4,3)=-vold_l*gm1
969 : wl(4,4)=gm1
970 : c------------------------------
971 : c------------------------------
972 : wr(1,1)=1.0d0
973 : wr(1,2)=0.0d0
974 : wr(1,3)=0.0d0
975 : wr(1,4)=0.0d0
976 : c------------------------------
977 : wr(2,1)=-uold_r/rold_r
978 : wr(2,2)=1.0d0/rold_r
979 : wr(2,3)=0.0d0
980 : wr(2,4)=0.0d0
981 : c------------------------------
982 : wr(3,1)=-vold_r/rold_r
983 : wr(3,2)=0.0d0
984 : wr(3,3)=1.0d0/rold_r
985 : wr(3,4)=0.0d0
986 : c------------------------------
987 : wr(4,1)=gm1*(uold_r*uold_r+vold_r*vold_r)/2.0d0
988 : wr(4,2)=-uold_r*gm1
989 : wr(4,3)=-vold_r*gm1
990 : wr(4,4)=gm1
991 : c----------------------------------------------
992 : do 1 i=1,4
993 : do 2 j=1,4
994 : jtl(i,j)=0.0d0
995 : jtr(i,j)=0.0d0
996 : do 3 k=1,4
997 : jtl(i,j)=jtl(i,j)+jl(i,k)*wl(k,j)
998 : jtr(i,j)=jtr(i,j)+jr(i,k)*wr(k,j)
999 : 3 continue
1000 : 2 continue
1001 : 1 continue
1002 : c----------------------------------------------------------------------
1003 : return
1004 : end
1005 :
1006 :
1007 :
1008 :
1009 :
1010 :
1011 :
1012 :
1013 :
1014 :
1015 :
1016 :
1017 :
1018 :
1019 :
1020 :
1021 :
1022 :
1023 :
1024 :
1025 :
1026 :
1027 :
1028 :
1029 :
1030 :
1031 :
1032 :
1033 :
1034 :
1035 :
1036 :
1037 :
© Cast3M 2003 - All rights reserved.
Disclaimer