1 : C CONMSP SOURCE PV 09/03/12 21:18:03 6325
2 : SUBROUTINE CONMSP(NSP,JLL,JRR,WVEC_L,WVEC_R,NVECT,TVECT,
3 : & mpyn,lrecp,lrecv,nlcg,nlcd)
4 : C************************************************************************
5 : C
6 : C PROJET : CASTEM 2000
7 : C
8 : C NOM : CONMSP ('convection for multispecies')
9 : C
10 : C DESCRIPTION : Voir KONMSP (appele par KONMSP)
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
25 : c EQUATIONS: 2D Euler equations of gas dynamics - MULTISPECIES GAS
26 : c
27 : c
28 : c REFERENCE: JCP, 129, 364-382 (1996)
29 : c " A Sequel to AUSM: AUSM+ ".
30 : c----------------------------------------------------------------------
31 : c INPUT:
32 : c
33 : c alpha -- parameter of the AUSM+ scheme in the Pressure function;
34 : c ( -3/4 <= alpha <= 3/16 ) (imposed as a parameter)
35 : c
36 : c beta -- parameter of the AUSM+ scheme in the Mach function;
37 : c ( -1/16 <= beta <= 1/2 ) (imposed as a parameter)
38 : c
39 : c nsp -- number of species (total);
40 : c
41 : c wvec_l -- vector of the primitive variables
42 : c (rho,u_x,u_y,p) at the left cell;
43 : c
44 : c wvec_r -- vector of the primitive variables
45 : c (rho,u_x,u_y,p) at the right cell;
46 : c
47 : c nvect -- normal vector to the interface (2 components in 2D);
48 : c
49 : c tvect -- tangential vector to the interface;
50 : c
51 : c mpyn -- pointer to the vectors of the primitive variables
52 : c (Y_1,Y_2,...Y_(nsp-1)) at the left and the right cells;
53 : c
54 : c lrecp -- pointer to the vector of specific heats at constant pressure
55 : c (size of the vector is equal to number of species (nsp));
56 : c
57 : c lrecv -- pointer to the vector of specific heats at constant volume
58 : c (size of the vector is equal to number of species (nsp));
59 : c
60 : c nlcg -- "local" number corresponding to the left cell;
61 : c
62 : c nlcd -- "local" number corresponding to the right cell;
63 : c----------------------------------------------------------------------
64 : c
65 : c OUTPUT:
66 : c
67 : c jll -- jakobian matrix (3+nsp) by (3+nsp) -
68 : c derivatives of the numerical
69 : c flux function with respect to the conservative variables
70 : c from the left cell;
71 : c
72 : c jrr -- jakobian matrix (3+nsp) by (3+nsp) -
73 : c derivatives of the numerical
74 : c flux function with respect to the conservative variables
75 : c from the right cell.
76 : c----------------------------------------------------------------------
77 : IMPLICIT INTEGER(I-N)
78 : integer nsp,jll,jrr,lrecp,lrecv,nlcg,nlcd
79 : real*8 wvec_l(4),wvec_r(4)
80 : real*8 nvect(2),tvect(2)
81 : real*8 alpha,beta
82 : real*8 gal,gar,gm1l,gm1r,temph
83 : real*8 n_x,n_y,t_x,t_y
84 : real*8 un_l, un_r, ut_l, ut_r
85 : real*8 ml,mr,Mplus,Mmin,mmid
86 : real*8 mpl_m, mmin_m,am
87 : real*8 rold_l,uold_l,vold_l,pold_l,eold_l
88 : real*8 rold_r,uold_r,vold_r,pold_r,eold_r
89 : real*8 Pplus,Pmin,pmid
90 : real*8 hr_l,hr_r,det,top,bot
91 : real*8 br1,br2,br3,br4,temp_l,temp_r,brac_l,brac_r
92 : real*8 aleft, arigh
93 : real*8 damr_l,damr_r,damu_l,damu_r
94 : real*8 damv_l,damv_r,damp_l,damp_r
95 : real*8 damg_l,damg_r
96 : real*8 dmlr_l,dmlr_r,dmlu_l,dmlu_r
97 : real*8 dmlv_l,dmlv_r,dmlp_l,dmlp_r
98 : real*8 dmrr_l,dmrr_r,dmru_l,dmru_r
99 : real*8 dmrv_l,dmrv_r,dmrp_l,dmrp_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 dMmr_l,dMmr_r,dMmu_l,dMmu_r
103 : real*8 dMmv_l,dMmv_r,dMmp_l,dMmp_r
104 : real*8 dmir_l,dmir_r,dmiu_l,dmiu_r
105 : real*8 dmiv_l,dmiv_r,dmip_l,dmip_r
106 : real*8 d3mr_l,d3mr_r,d3mu_l,d3mu_r
107 : real*8 d3mv_l,d3mv_r,d3mp_l,d3mp_r
108 : real*8 d2mr_l,d2mr_r,d2mu_l,d2mu_r
109 : real*8 d2mv_l,d2mv_r,d2mp_l,d2mp_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 dPmr_l,dPmr_r,dPmu_l,dPmu_r
113 : real*8 dPmv_l,dPmv_r,dPmp_l,dPmp_r
114 : real*8 dpir_l,dpir_r,dpiu_l,dpiu_r
115 : real*8 dpiv_l,dpiv_r,dpip_l,dpip_r
116 : integer i,j,k
117 : parameter(alpha = 0.1875D0, beta = 0.125D0)
118 : C------------------------------------------------------------
119 : -INC SMCHPOI
120 : POINTEUR MPYN.MPOVAL
121 : C-------------------------------------------------------------
122 : -INC SMLREEL
123 : POINTEUR MLRECP.MLREEL, MLRECV.MLREEL
124 : C-------------------------------------------------------------
125 : C******* Les fractionines massiques **************************
126 : C-------------------------------------------------------------
127 : SEGMENT FRAMAS
128 : REAL*8 YET(NSP)
129 : ENDSEGMENT
130 : POINTEUR YL.FRAMAS, YR.FRAMAS
131 : C-------------------------------------------------------
132 : C********** Les CP's and CV's ***********************
133 : C-------------------------------------------------------
134 : SEGMENT GCONST
135 : REAL*8 GC(NSP)
136 : ENDSEGMENT
137 : POINTEUR CP.GCONST, CV.GCONST
138 : C-------------------------------------------------------------
139 : C******** Segments for the elementary matrixes *************
140 : C-------------------------------------------------------------
141 : SEGMENT JACEL
142 : REAL*8 JAC(3+NSP,3+NSP)
143 : ENDSEGMENT
144 : POINTEUR JTL.JACEL, JTR.JACEL, JL.JACEL, JR.JACEL,
145 : & WL.JACEL, WR.JACEL
146 : c----------------------------------------
147 : SEGINI JTL
148 : SEGINI JTR
149 : SEGINI JL
150 : SEGINI JR
151 : SEGINI WL
152 : SEGINI WR
153 : C-------------------------------------------------------------
154 : C********** Segments for the vectors ***********************
155 : C-------------------------------------------------------------
156 : SEGMENT VECEL
157 : REAL*8 VV(NSP)
158 : ENDSEGMENT
159 : POINTEUR DMLY_L.VECEL, DMLY_R.VECEL,
160 : & dmry_l.vecel, dmry_r.vecel,
161 : & dMpy_l.vecel, dMpy_r.vecel,
162 : & dMmy_l.vecel, dMmy_r.vecel,
163 : & dmiy_l.vecel, dmiy_r.vecel,
164 : & d3my_l.vecel, d3my_r.vecel,
165 : & d2my_l.vecel, d2my_r.vecel,
166 : & dPpy_l.vecel, dPpy_r.vecel,
167 : & dPmy_l.vecel, dPmy_r.vecel,
168 : & dpiy_l.vecel, dpiy_r.vecel,
169 : & dgdyl.vecel, dgdyr.vecel
170 : C----------------------------------------------
171 : SEGINI DMLY_L, DMLY_R,
172 : & dmry_l, dmry_r,
173 : & dMpy_l, dMpy_r,
174 : & dMmy_l, dMmy_r,
175 : & dmiy_l, dmiy_r,
176 : & d3my_l, d3my_r,
177 : & d2my_l, d2my_r,
178 : & dPpy_l, dPpy_r,
179 : & dPmy_l, dPmy_r,
180 : & dpiy_l, dpiy_r,
181 : & dgdyl, dgdyr
182 : C-------------------------------------------------------------
183 : C********** Segments for the flux-vector *******************
184 : C-------------------------------------------------------------
185 : SEGMENT FUNEL
186 : REAL*8 FU(3+NSP)
187 : ENDSEGMENT
188 : POINTEUR f.funel
189 : C-------------------------
190 : SEGINI f
191 : C------------------------------------------------------------
192 : SEGINI YL, YR
193 : SEGACT MPYN
194 : DO 100 I=1,(NSP-1)
195 : YL.YET(I)=MPYN.VPOCHA(NLCG,I)
196 : YR.YET(I)=MPYN.VPOCHA(NLCD,I)
197 : 100 CONTINUE
198 : C----------------------------------------
199 : SEGINI CP, CV
200 : MLRECP = LRECP
201 : MLRECV = LRECV
202 : SEGACT MLRECP, MLRECV
203 : DO 101 I=1,(NSP-1)
204 : CP.GC(I)=MLRECP.PROG(I)
205 : CV.GC(I)=MLRECV.PROG(I)
206 : 101 CONTINUE
207 : CP.GC(NSP)=MLRECP.PROG(NSP)
208 : CV.GC(NSP)=MLRECV.PROG(NSP)
209 : c-------------------------------------------------------------
210 : c Computing GAMMA at the left cell and its derivatives
211 : c with respect to the primitive variables Y_i
212 : c-------------------------------------------------------------
213 : top=0.0D0
214 : bot=0.0D0
215 : do 40 i=1,(nsp-1)
216 : top=top+yl.yet(i)*(cp.gc(i)-cp.gc(nsp))
217 : bot=bot+yl.yet(i)*(cv.gc(i)-cv.gc(nsp))
218 : 40 continue
219 : top=cp.gc(nsp)+top
220 : bot=cv.gc(nsp)+bot
221 : gal=top/bot
222 : gm1l=gal-1.0d0
223 : c-------------------------------------------------------------
224 : do 41 i=1,(nsp-1)
225 : dgdyl.vv(i)=(cp.gc(i)-cp.gc(nsp)-
226 : & gal*(cv.gc(i)-cv.gc(nsp)))/bot
227 : 41 continue
228 : c-------------------------------------------------------------
229 : c Computing GAMMA at the right cell and its derivatives
230 : c with respect to the primitive variables Y_i
231 : c-------------------------------------------------------------
232 : top=0.0D0
233 : bot=0.0D0
234 : do 42 i=1,(nsp-1)
235 : top=top+yr.yet(i)*(cp.gc(i)-cp.gc(nsp))
236 : bot=bot+yr.yet(i)*(cv.gc(i)-cv.gc(nsp))
237 : 42 continue
238 : top=cp.gc(nsp)+top
239 : bot=cv.gc(nsp)+bot
240 : gar=top/bot
241 : gm1r=gar-1.0d0
242 : c-------------------------------------------------------------
243 : do 43 i=1,(nsp-1)
244 : dgdyr.vv(i)=(cp.gc(i)-cp.gc(nsp)-
245 : & gar*(cv.gc(i)-cv.gc(nsp)))/bot
246 : 43 continue
247 : c-------------------------------------------------------------
248 : n_x=nvect(1)
249 : n_y=nvect(2)
250 : t_x=tvect(1)
251 : t_y=tvect(2)
252 : c----------------------------
253 : c----------------------------
254 : rold_l=wvec_l(1)
255 : uold_l=wvec_l(2)
256 : vold_l=wvec_l(3)
257 : pold_l=wvec_l(4)
258 : c-----------------------
259 : rold_r=wvec_r(1)
260 : uold_r=wvec_r(2)
261 : vold_r=wvec_r(3)
262 : pold_r=wvec_r(4)
263 : c------------------------------------------------------------------
264 : c Computation of the specific total energy on the left and right.
265 : c------------------------------------------------------------------
266 : eold_l=(uold_l*uold_l+vold_l*vold_l)/2.0d0
267 : eold_l=eold_l+pold_l/(gm1l*rold_l)
268 : eold_r=(uold_r*uold_r+vold_r*vold_r)/2.0d0
269 : eold_r=eold_r+pold_r/(gm1r*rold_r)
270 : c-------------------------------------------------------------------
271 : c Computation of the speed of sound and its derivatives;
272 : c numerical speed of sound at the interface is taken as an average
273 : c of the speeds of sounds of the neighbouring cells
274 : c-------------------------------------------------------------------
275 : aleft=sqrt(gal*pold_l/rold_l)
276 : arigh=sqrt(gar*pold_r/rold_r)
277 : am=0.5d0*(aleft+arigh)
278 : c-------------------------------------------------------------------
279 : damr_r=-arigh/(4.0d0*rold_r)
280 : damu_r=0.0d0
281 : damv_r=0.0d0
282 : damp_r=gar/(4.0d0*arigh*rold_r)
283 : damg_r=arigh/(4.0d0*gar)
284 : c-----------------------
285 : damr_l=-aleft/(4.0d0*rold_l)
286 : damu_l=0.0d0
287 : damv_l=0.0d0
288 : damp_l=gal/(4.0d0*aleft*rold_l)
289 : damg_l=aleft/(4.0d0*gal)
290 : c-------------------------------------------------------------------
291 : c Computing numerical Mach number and its derivatives,
292 : c see p.370, under (A1).
293 : c-------------------------------------------------------------------
294 : un_l=uold_l*n_x+vold_l*n_y
295 : un_r=uold_r*n_x+vold_r*n_y
296 : ut_l=uold_l*t_x+vold_l*t_y
297 : ut_r=uold_r*t_x+vold_r*t_y
298 : c-------------------------------------------------------------------
299 : ml=un_l/am
300 : mr=un_r/am
301 : c-------------------------------------------------------------------
301 : c Mplus and Mmin are calligraphic lettes M+ and M- from the paper,
302 : c see (19a) and (19b), p.367.
303 : c-------------------------------------------------------------------
304 : if(abs(ml) .ge. 1.0d0) then
305 : Mplus=(ml+abs(ml))/2.0d0
306 : else
307 : Mplus=(ml+1.0d0)*(ml+1.0d0)/4.0d0
308 : Mplus=Mplus+beta*(ml*ml-1.0d0)*(ml*ml-1.0d0)
309 : endif
310 : c-------------------------------------------------------------------
311 : if(abs(mr) .ge. 1.0d0) then
312 : Mmin=(mr-abs(mr))/2.0d0
313 : else
314 : Mmin=-(mr-1.0d0)*(mr-1.0d0)/4.0d0
315 : Mmin=Mmin-beta*(mr*mr-1.0d0)*(mr*mr-1.0d0)
316 : endif
317 : c-------------------------------------------------------------------
318 : c Derivatives of ml and mr with respect to both sets of primitive
319 : c variables: left and right.
320 : c-------------------------------------------------------------------
321 : temp_l=-un_l/(am*am)
322 : temp_r=-un_r/(am*am)
323 : c--------
324 : dmlr_l=temp_l*damr_l
325 : dmlr_r=temp_l*damr_r
326 : dmrr_l=temp_r*damr_l
327 : dmrr_r=temp_r*damr_r
328 : c--------
329 : dmlu_l=n_x/am+temp_l*damu_l
330 : dmlu_r=temp_l*damu_r
331 : dmru_l=temp_r*damu_l
332 : dmru_r=n_x/am+temp_r*damu_r
333 : c--------
334 : dmlv_l=n_y/am+temp_l*damv_l
335 : dmlv_r=temp_l*damv_r
336 : dmrv_l=temp_r*damv_l
337 : dmrv_r=n_y/am+temp_r*damv_r
338 : c--------
339 : dmlp_l=temp_l*damp_l
340 : dmlp_r=temp_l*damp_r
341 : dmrp_l=temp_r*damp_l
342 : dmrp_r=temp_r*damp_r
343 : c--------
344 : do 44 i=1,(nsp-1)
345 : dmly_l.vv(i)=temp_l*damg_l*dgdyl.vv(i)
346 : dmly_r.vv(i)=temp_l*damg_r*dgdyr.vv(i)
347 : dmry_l.vv(i)=temp_r*damg_l*dgdyl.vv(i)
348 : dmry_r.vv(i)=temp_r*damg_r*dgdyr.vv(i)
349 : 44 continue
350 : c-----------------------------------------------------------
351 : c mmid is m_{1/2} (notation as in the paper, see (13),p.366)
352 : c-----------------------------------------------------------
353 : mmid=Mplus+Mmin
354 : c-----------------------------------------------------------
355 : c Computing the derivatives of M+ and M-
356 : c-----------------------------------------------------------
357 : if(ml .ge. 1.0d0) then
358 : dMpr_l=dmlr_l
359 : dMpu_l=dmlu_l
360 : dMpv_l=dmlv_l
361 : dMpp_l=dmlp_l
362 : do 45 i=1,(nsp-1)
363 : dMpy_l.vv(i)=dmly_l.vv(i)
364 : 45 continue
365 : c--------------------
366 : dMpr_r=dmlr_r
367 : dMpu_r=dmlu_r
368 : dMpv_r=dmlv_r
369 : dMpp_r=dmlp_r
370 : do 46 i=1,(nsp-1)
371 : dMpy_r.vv(i)=dmly_r.vv(i)
372 : 46 continue
373 : else
374 : if((ml .gt. -1.0d0) .and. (ml .lt. 1.0d0)) then
375 : temph=(ml+1.0d0)/2.0d0
376 : dMpr_l=(temph+4.0d0*beta*ml*(ml*ml-1.0d0))*dmlr_l
377 : dMpu_l=(temph+4.0d0*beta*ml*(ml*ml-1.0d0))*dmlu_l
378 : dMpv_l=(temph+4.0d0*beta*ml*(ml*ml-1.0d0))*dmlv_l
379 : dMpp_l=(temph+4.0d0*beta*ml*(ml*ml-1.0d0))*dmlp_l
380 : do 47 i=1,(nsp-1)
381 : dMpy_l.vv(i)=(temph+4.0d0*beta*ml*
382 : & (ml*ml-1.0d0))*dmly_l.vv(i)
383 : 47 continue
384 : c--------------------
385 : dMpr_r=(temph+4.0d0*beta*ml*(ml*ml-1.0d0))*dmlr_r
386 : dMpu_r=(temph+4.0d0*beta*ml*(ml*ml-1.0d0))*dmlu_r
387 : dMpv_r=(temph+4.0d0*beta*ml*(ml*ml-1.0d0))*dmlv_r
388 : dMpp_r=(temph+4.0d0*beta*ml*(ml*ml-1.0d0))*dmlp_r
389 : do 48 i=1,(nsp-1)
390 : dMpy_r.vv(i)=(temph+4.0d0*beta*ml*
391 : & (ml*ml-1.0d0))*dmly_r.vv(i)
392 : 48 continue
393 : else
394 : dMpr_l=0.0d0
395 : dMpu_l=0.0d0
396 : dMpv_l=0.0d0
397 : dMpp_l=0.0d0
398 : do 49 i=1,(nsp-1)
399 : dMpy_l.vv(i)=0.0d0
400 : 49 continue
401 : c---------------------
402 : dMpr_r=0.0d0
403 : dMpu_r=0.0d0
404 : dMpv_r=0.0d0
405 : dMpp_r=0.0d0
406 : do 50 i=1,(nsp-1)
407 : dMpy_r.vv(i)=0.0d0
408 : 50 continue
409 : endif
410 : endif
411 : c-----------------------------------------------------------
412 : if(mr .ge. 1.0d0) then
413 : dMmr_l=0.0d0
414 : dMmu_l=0.0d0
415 : dMmv_l=0.0d0
416 : dMmp_l=0.0d0
417 : do 51 i=1,(nsp-1)
418 : dMmy_l.vv(i)=0.0d0
419 : 51 continue
420 : c---------------------
421 : dMmr_r=0.0d0
422 : dMmu_r=0.0d0
423 : dMmv_r=0.0d0
424 : dMmp_r=0.0d0
425 : do 52 i=1,(nsp-1)
426 : dMmy_r.vv(i)=0.0d0
427 : 52 continue
428 : else
429 : if((mr .gt. -1.0d0) .and. (mr .lt. 1.0d0)) then
430 : temph=(-mr+1.0d0)/2.0d0
431 : dMmr_l=(temph-4.0d0*beta*mr*(mr*mr-1.0d0))*dmrr_l
432 : dMmu_l=(temph-4.0d0*beta*mr*(mr*mr-1.0d0))*dmru_l
433 : dMmv_l=(temph-4.0d0*beta*mr*(mr*mr-1.0d0))*dmrv_l
434 : dMmp_l=(temph-4.0d0*beta*mr*(mr*mr-1.0d0))*dmrp_l
435 : do 53 i=1,(nsp-1)
436 : dMmy_l.vv(i)=(temph-4.0d0*beta*mr*
437 : & (mr*mr-1.0d0))*dmry_l.vv(i)
438 : 53 continue
439 : c--------------------
440 : dMmr_r=(temph-4.0d0*beta*mr*(mr*mr-1.0d0))*dmrr_r
441 : dMmu_r=(temph-4.0d0*beta*mr*(mr*mr-1.0d0))*dmru_r
442 : dMmv_r=(temph-4.0d0*beta*mr*(mr*mr-1.0d0))*dmrv_r
443 : dMmp_r=(temph-4.0d0*beta*mr*(mr*mr-1.0d0))*dmrp_r
444 : do 54 i=1,(nsp-1)
445 : dMmy_r.vv(i)=(temph-4.0d0*beta*mr*
446 : & (mr*mr-1.0d0))*dmry_r.vv(i)
447 : 54 continue
448 : else
449 : dMmr_l=dmrr_l
450 : dMmu_l=dmru_l
451 : dMmv_l=dmrv_l
452 : dMmp_l=dmrp_l
453 : do 55 i=1,(nsp-1)
454 : dMmy_l.vv(i)=dmry_l.vv(i)
455 : 55 continue
456 : c--------------------
457 : dMmr_r=dmrr_r
458 : dMmu_r=dmru_r
459 : dMmv_r=dmrv_r
460 : dMmp_r=dmrp_r
461 : do 56 i=1,(nsp-1)
462 : dMmy_r.vv(i)=dmry_r.vv(i)
463 : 56 continue
464 : endif
465 : endif
466 : c-----------------------------------------------------------------
467 : c computing the derivatives of m_{1/2} (notation as in the paper)
468 : c-----------------------------------------------------------------
469 : dmir_l=dMpr_l+dMmr_l
470 : dmir_r=dMpr_r+dMmr_r
471 : c-------------
472 : dmiu_l=dMpu_l+dMmu_l
473 : dmiu_r=dMpu_r+dMmu_r
474 : c-------------
475 : dmiv_l=dMpv_l+dMmv_l
476 : dmiv_r=dMpv_r+dMmv_r
477 : c-------------
478 : dmip_l=dMpp_l+dMmp_l
479 : dmip_r=dMpp_r+dMmp_r
480 : c-------------
481 : do 57 i=1,(nsp-1)
482 : dmiy_l.vv(i)=dMpy_l.vv(i)+dMmy_l.vv(i)
483 : dmiy_r.vv(i)=dMpy_r.vv(i)+dMmy_r.vv(i)
484 : 57 continue
485 : c----------------------------------------------------------------
486 : c computing the main convective variables and their derivatives
487 : c mpl_m is m^{+}_{1/2} (paper's notation) and
488 : c mmin_m is m^{-}_{1/2} (paper's notation), see (A2) on p.370.
489 : c----------------------------------------------------------------
490 : if(mmid .ge. 0.0d0) then
491 : mpl_m = mmid
492 : d2mr_l=dmir_l
493 : d2mu_l=dmiu_l
494 : d2mv_l=dmiv_l
495 : d2mp_l=dmip_l
496 : do 58 i=1,(nsp-1)
497 : d2my_l.vv(i)=dmiy_l.vv(i)
498 : 58 continue
499 : c------------
500 : d2mr_r=dmir_r
501 : d2mu_r=dmiu_r
502 : d2mv_r=dmiv_r
503 : d2mp_r=dmip_r
504 : do 59 i=1,(nsp-1)
505 : d2my_r.vv(i)=dmiy_r.vv(i)
506 : 59 continue
507 : c------------
508 : else
509 : mpl_m = 0.0d0
510 : d2mr_l=0.0d0
511 : d2mu_l=0.0d0
512 : d2mv_l=0.0d0
513 : d2mp_l=0.0d0
514 : do 60 i=1,(nsp-1)
515 : d2my_l.vv(i)=0.0d0
516 : 60 continue
517 : c------------
518 : d2mr_r=0.0d0
519 : d2mu_r=0.0d0
520 : d2mv_r=0.0d0
521 : d2mp_r=0.0d0
522 : do 61 i=1,(nsp-1)
523 : d2my_r.vv(i)=0.0d0
524 : 61 continue
525 : endif
526 : c---------------------------------------------
527 : if(mmid .le. 0.0d0) then
528 : mmin_m = mmid
529 : d3mr_l=dmir_l
530 : d3mu_l=dmiu_l
531 : d3mv_l=dmiv_l
532 : d3mp_l=dmip_l
533 : do 62 i=1,(nsp-1)
534 : d3my_l.vv(i)=dmiy_l.vv(i)
535 : 62 continue
536 : c------------
537 : d3mr_r=dmir_r
538 : d3mu_r=dmiu_r
539 : d3mv_r=dmiv_r
540 : d3mp_r=dmip_r
541 : do 63 i=1,(nsp-1)
542 : d3my_r.vv(i)=dmiy_r.vv(i)
543 : 63 continue
544 : c------------
545 : else
546 : mmin_m = 0.0d0
547 : d3mr_l=0.0d0
548 : d3mu_l=0.0d0
549 : d3mv_l=0.0d0
550 : d3mp_l=0.0d0
551 : do 64 i=1,(nsp-1)
552 : d3my_l.vv(i)=0.0d0
553 : 64 continue
554 : c------------
555 : d3mr_r=0.0d0
556 : d3mu_r=0.0d0
557 : d3mv_r=0.0d0
558 : d3mp_r=0.0d0
559 : do 65 i=1,(nsp-1)
560 : d3my_r.vv(i)=0.0d0
561 : 65 continue
562 : endif
563 : c---------------------------------------------------------------
563 : c Computing the calligraphic P+ and P- with their derivatives,
564 : c see (21a) & (21b) on p.368.
565 : c---------------------------------------------------------------
566 : if(ml .ge. 1.0d0) then
567 : Pplus = 1.0d0
568 : else
569 : if((ml .gt. -1.0d0) .and. (ml .lt. 1.0d0)) then
570 : Pplus=(ml+1.0d0)*(ml+1.0d0)*(2.0d0-ml)/4.0d0
571 : Pplus=Pplus+alpha*ml*(ml*ml-1.0d0)*(ml*ml-1.0d0)
572 : else
573 : Pplus = 0.0d0
574 : endif
575 : endif
576 : c---------------------------------------------------------------
577 : if(mr .ge. 1.0d0) then
578 : Pmin = 0.0d0
579 : else
580 : if((mr .gt. -1.0d0) .and. (mr .lt. 1.0d0)) then
581 : Pmin=(mr-1.0d0)*(mr-1.0d0)*(2.0d0+mr)/4.0d0
582 : Pmin=Pmin-alpha*mr*(mr*mr-1.0d0)*(mr*mr-1.0d0)
583 : else
584 : Pmin = 1.0d0
585 : endif
586 : endif
587 : c---------------------------------------------------------------
588 : brac_l=(ml+1.0d0)*(2.0d0-ml)/2.0d0-(ml+1.0d0)*(ml+1.0d0)/4.0d0
589 : brac_l=brac_l+alpha*(ml*ml-1.0d0)*(ml*ml-1.0d0)
590 : brac_l=brac_l+4.0d0*alpha*ml*ml*(ml*ml-1.0d0)
591 : c--------------
592 : brac_r=(mr-1.0d0)*(2.0d0+mr)/2.0d0+(mr-1.0d0)*(mr-1.0d0)/4.0d0
593 : brac_r=brac_r-alpha*(mr*mr-1.0d0)*(mr*mr-1.0d0)
594 : brac_r=brac_r-4.0d0*alpha*mr*mr*(mr*mr-1.0d0)
595 : c---------------------------------------------------------------
596 : if((ml .gt. -1.0d0) .and. (ml .lt. 1.0d0)) then
597 : dPpr_l=brac_l*dmlr_l
598 : dPpr_r=brac_l*dmlr_r
599 : c------------
600 : dPpu_l=brac_l*dmlu_l
601 : dPpu_r=brac_l*dmlu_r
602 : c------------
603 : dPpv_l=brac_l*dmlv_l
604 : dPpv_r=brac_l*dmlv_r
605 : c------------
606 : dPpp_l=brac_l*dmlp_l
607 : dPpp_r=brac_l*dmlp_r
608 : c------------
609 : do 66 i=1,(nsp-1)
610 : dPpy_l.vv(i)=brac_l*dmly_l.vv(i)
611 : dPpy_r.vv(i)=brac_l*dmly_r.vv(i)
612 : 66 continue
613 : c------------
614 : else
615 : dPpr_l=0.0d0
616 : dPpr_r=0.0d0
617 : c-----------
618 : dPpu_l=0.0d0
619 : dPpu_r=0.0d0
620 : c-----------
621 : dPpv_l=0.0d0
622 : dPpv_r=0.0d0
623 : c-----------
624 : dPpp_l=0.0d0
625 : dPpp_r=0.0d0
626 : c-----------
627 : do 67 i=1,(nsp-1)
628 : dPpy_l.vv(i)=0.0d0
629 : dPpy_r.vv(i)=0.0d0
630 : 67 continue
631 : c-----------
632 : endif
633 : c---------------------------------------------------------------
634 : if((mr .gt. -1.0d0) .and. (mr .lt. 1.0d0)) then
635 : dPmr_l=brac_r*dmrr_l
636 : dPmr_r=brac_r*dmrr_r
637 : c------------
638 : dPmu_l=brac_r*dmru_l
639 : dPmu_r=brac_r*dmru_r
640 : c------------
641 : dPmv_l=brac_r*dmrv_l
642 : dPmv_r=brac_r*dmrv_r
643 : c------------
644 : dPmp_l=brac_r*dmrp_l
645 : dPmp_r=brac_r*dmrp_r
646 : c------------
647 : do 68 i=1,(nsp-1)
648 : dPmy_l.vv(i)=brac_r*dmry_l.vv(i)
649 : dPmy_r.vv(i)=brac_r*dmry_r.vv(i)
650 : 68 continue
651 : c------------
652 : else
653 : dPmr_l=0.0d0
654 : dPmr_r=0.0d0
655 : c-----------
656 : dPmu_l=0.0d0
657 : dPmu_r=0.0d0
658 : c-----------
659 : dPmv_l=0.0d0
660 : dPmv_r=0.0d0
661 : c-----------
662 : dPmp_l=0.0d0
663 : dPmp_r=0.0d0
664 : c-----------
665 : do 69 i=1,(nsp-1)
666 : dPmy_l.vv(i)=0.0d0
667 : dPmy_r.vv(i)=0.0d0
668 : 69 continue
669 : c-----------
670 : endif
671 : c-------------------------------------------------------------------
672 : c computing pmid -- p_{1/2} and its derivatives, see (20b), p.367.
673 : c-------------------------------------------------------------------
674 : pmid=Pplus*pold_l + Pmin*pold_r
675 : dpir_l=dPpr_l*pold_l+dPmr_l*pold_r
676 : dpiu_l=dPpu_l*pold_l+dPmu_l*pold_r
677 : dpiv_l=dPpv_l*pold_l+dPmv_l*pold_r
678 : dpip_l=dPpp_l*pold_l+Pplus+dPmp_l*pold_r
679 : do 70 i=1,(nsp-1)
680 : dpiy_l.vv(i)=dPpy_l.vv(i)*pold_l+dPmy_l.vv(i)*pold_r
681 : 70 continue
682 : c----------------------------
683 : dpir_r=dPpr_r*pold_l+dPmr_r*pold_r
684 : dpiu_r=dPpu_r*pold_l+dPmu_r*pold_r
685 : dpiv_r=dPpv_r*pold_l+dPmv_r*pold_r
686 : dpip_r=dPpp_r*pold_l+Pmin+dPmp_r*pold_r
687 : do 71 i=1,(nsp-1)
688 : dpiy_r.vv(i)=dPpy_r.vv(i)*pold_l+dPmy_r.vv(i)*pold_r
689 : 71 continue
690 : c---------------------------------------------------------------------
691 : c computing JACOBIAN as a derivative of the numerical flux function
692 : c with respect to the primitive variables.
693 : c Notation: jl(i,j) --- is the derivative of the i-component of the
694 : c flux function with respect to the j-component of the
695 : c vector of primitive variables of the left state.
696 : c jr(i,j) --- is the derivative of the i-component of the
697 : c flux function with respect to the j-component of the
698 : c vector of primitive variables of the right state.
699 : c---------------------------------------------------------------------
700 : f.fu(1)=am*(mpl_m*rold_l+mmin_m*rold_r)
701 : c---------------------------------------------------------------------
702 : jl.jac(1,1)=damr_l*f.fu(1)/am+am*(d2mr_l*rold_l+mpl_m)
703 : jl.jac(1,1)=jl.jac(1,1)+am*d3mr_l*rold_r
704 : jl.jac(1,2)=damu_l*f.fu(1)/am+am*(d2mu_l*rold_l+d3mu_l*rold_r)
705 : jl.jac(1,3)=damv_l*f.fu(1)/am+am*(d2mv_l*rold_l+d3mv_l*rold_r)
706 : jl.jac(1,4)=damp_l*f.fu(1)/am+am*(d2mp_l*rold_l+d3mp_l*rold_r)
707 : do 72 i=1,(nsp-1)
708 : jl.jac(1,4+i)=damg_l*dgdyl.vv(i)*f.fu(1)/am+
709 : & am*(d2my_l.vv(i)*rold_l+d3my_l.vv(i)*rold_r)
710 : 72 continue
711 : c------------------------------------
712 : jr.jac(1,1)=damr_r*f.fu(1)/am+am*(d2mr_r*rold_l+mmin_m)
713 : jr.jac(1,1)=jr.jac(1,1)+am*d3mr_r*rold_r
714 : jr.jac(1,2)=damu_r*f.fu(1)/am+am*(d2mu_r*rold_l+d3mu_r*rold_r)
715 : jr.jac(1,3)=damv_r*f.fu(1)/am+am*(d2mv_r*rold_l+d3mv_r*rold_r)
716 : jr.jac(1,4)=damp_r*f.fu(1)/am+am*(d2mp_r*rold_l+d3mp_r*rold_r)
717 : do 73 i=1,(nsp-1)
718 : jr.jac(1,4+i)=damg_r*dgdyr.vv(i)*f.fu(1)/am+
719 : & am*(d2my_r.vv(i)*rold_l+d3my_r.vv(i)*rold_r)
720 : 73 continue
721 : c------------------------------------
722 : c------------------------------------
723 : br1=mpl_m*rold_l*un_l+mmin_m*rold_r*un_r
724 : br2=mpl_m*rold_l*ut_l+mmin_m*rold_r*ut_r
725 : f.fu(2)=n_x*(am*br1+pmid)+am*t_x*br2
726 : c------------------
727 : det=n_x*t_y-n_y*t_x
728 : c---------------------------------------------------------
729 : br3=d2mr_l*rold_l*un_l+mpl_m*un_l+d3mr_l*rold_r*un_r
730 : br4=d2mr_l*rold_l*ut_l+mpl_m*ut_l+d3mr_l*rold_r*ut_r
731 : jl.jac(2,1)=n_x*(damr_l*br1+am*br3+dpir_l)
732 : jl.jac(2,1)=jl.jac(2,1)+t_x*damr_l*br2+am*t_x*br4
733 : c-------------------
734 : br3=d2mu_l*rold_l*un_l+mpl_m*rold_l*t_y/det
735 : br3=br3+d3mu_l*rold_r*un_r
736 : br4=d2mu_l*rold_l*ut_l+mpl_m*rold_l*(-n_y)/det
737 : br4=br4+d3mu_l*rold_r*ut_r
738 : jl.jac(2,2)=n_x*(damu_l*br1+am*br3+dpiu_l)
739 : jl.jac(2,2)=jl.jac(2,2)+t_x*damu_l*br2+am*t_x*br4
740 : c-------------------
741 : br3=d2mv_l*rold_l*un_l+mpl_m*rold_l*(-t_x)/det
742 : br3=br3+d3mv_l*rold_r*un_r
743 : br4=d2mv_l*rold_l*ut_l+mpl_m*rold_l*n_x/det
744 : br4=br4+d3mv_l*rold_r*ut_r
745 : jl.jac(2,3)=n_x*(damv_l*br1+am*br3+dpiv_l)
746 : jl.jac(2,3)=jl.jac(2,3)+t_x*damv_l*br2+am*t_x*br4
747 : c-------------------
748 : br3=d2mp_l*rold_l*un_l+d3mp_l*rold_r*un_r
749 : br4=d2mp_l*rold_l*ut_l+d3mp_l*rold_r*ut_r
750 : jl.jac(2,4)=n_x*(damp_l*br1+am*br3+dpip_l)
751 : jl.jac(2,4)=jl.jac(2,4)+t_x*damp_l*br2+am*t_x*br4
752 : c-------------------
753 : do 74 i=1,(nsp-1)
754 : br3=d2my_l.vv(i)*rold_l*un_l+d3my_l.vv(i)*rold_r*un_r
755 : br4=d2my_l.vv(i)*rold_l*ut_l+d3my_l.vv(i)*rold_r*ut_r
756 : jl.jac(2,4+i)=n_x*(damg_l*dgdyl.vv(i)*br1+
757 : & am*br3+dpiy_l.vv(i))
758 : jl.jac(2,4+i)=jl.jac(2,4+i)+
759 : & t_x*damg_l*dgdyl.vv(i)*br2+am*t_x*br4
760 : 74 continue
761 : c-------------------------------------------------------------
762 : br3=d2mr_r*rold_l*un_l+mmin_m*un_r+d3mr_r*rold_r*un_r
763 : br4=d2mr_r*rold_l*ut_l+mmin_m*ut_r+d3mr_r*rold_r*ut_r
764 : jr.jac(2,1)=n_x*(damr_r*br1+am*br3+dpir_r)
765 : jr.jac(2,1)=jr.jac(2,1)+t_x*damr_r*br2+am*t_x*br4
766 : c-------------------
767 : br3=d2mu_r*rold_l*un_l+mmin_m*rold_r*t_y/det
768 : br3=br3+d3mu_r*rold_r*un_r
769 : br4=d2mu_r*rold_l*ut_l+mmin_m*rold_r*(-n_y/det)
770 : br4=br4+d3mu_r*rold_r*ut_r
771 : jr.jac(2,2)=n_x*(damu_r*br1+am*br3+dpiu_r)
772 : jr.jac(2,2)=jr.jac(2,2)+t_x*damu_r*br2+am*t_x*br4
773 : c-------------------
774 : br3=d2mv_r*rold_l*un_l+mmin_m*rold_r*(-t_x/det)
775 : br3=br3+d3mv_r*rold_r*un_r
776 : br4=d2mv_r*rold_l*ut_l+mmin_m*rold_r*n_x/det
777 : br4=br4+d3mv_r*rold_r*ut_r
778 : jr.jac(2,3)=n_x*(damv_r*br1+am*br3+dpiv_r)
779 : jr.jac(2,3)=jr.jac(2,3)+t_x*damv_r*br2+am*t_x*br4
780 : c-------------------
781 : br3=d2mp_r*rold_l*un_l+d3mp_r*rold_r*un_r
782 : br4=d2mp_r*rold_l*ut_l+d3mp_r*rold_r*ut_r
783 : jr.jac(2,4)=n_x*(damp_r*br1+am*br3+dpip_r)
784 : jr.jac(2,4)=jr.jac(2,4)+t_x*damp_r*br2+am*t_x*br4
785 : c-------------------
786 : do 75 i=1,(nsp-1)
787 : br3=d2my_r.vv(i)*rold_l*un_l+d3my_r.vv(i)*rold_r*un_r
788 : br4=d2my_r.vv(i)*rold_l*ut_l+d3my_r.vv(i)*rold_r*ut_r
789 : jr.jac(2,4+i)=n_x*(damg_r*dgdyr.vv(i)*br1+
790 : & am*br3+dpiy_r.vv(i))
791 : jr.jac(2,4+i)=jr.jac(2,4+i)+
792 : & t_x*damg_r*dgdyr.vv(i)*br2+am*t_x*br4
793 : 75 continue
794 : c-------------------------------------------------------------
795 : c-------------------------------------------------------------
796 : br1=mpl_m*rold_l*un_l+mmin_m*rold_r*un_r
797 : br2=mpl_m*rold_l*ut_l+mmin_m*rold_r*ut_r
798 : f.fu(3)=n_y*(am*br1+pmid)+am*t_y*br2
799 : br3=d2mr_l*rold_l*un_l+mpl_m*un_l+d3mr_l*rold_r*un_r
800 : br4=d2mr_l*rold_l*ut_l+mpl_m*ut_l+d3mr_l*rold_r*ut_r
801 : jl.jac(3,1)=n_y*(damr_l*br1+am*br3+dpir_l)
802 : jl.jac(3,1)=jl.jac(3,1)+t_y*damr_l*br2+am*t_y*br4
803 : c-------------------
804 : br3=d2mu_l*rold_l*un_l+mpl_m*rold_l*t_y/det
805 : br3=br3+d3mu_l*rold_r*un_r
806 : br4=d2mu_l*rold_l*ut_l+mpl_m*rold_l*(-n_y/det)
807 : br4=br4+d3mu_l*rold_r*ut_r
808 : jl.jac(3,2)=n_y*(damu_l*br1+am*br3+dpiu_l)
809 : jl.jac(3,2)=jl.jac(3,2)+t_y*damu_l*br2+am*t_y*br4
810 : c-------------------
811 : br3=d2mv_l*rold_l*un_l+mpl_m*rold_l*(-t_x/det)
812 : br3=br3+d3mv_l*rold_r*un_r
813 : br4=d2mv_l*rold_l*ut_l+mpl_m*rold_l*n_x/det
814 : br4=br4+d3mv_l*rold_r*ut_r
815 : jl.jac(3,3)=n_y*(damv_l*br1+am*br3+dpiv_l)
816 : jl.jac(3,3)=jl.jac(3,3)+t_y*damv_l*br2+am*t_y*br4
817 : c-------------------
818 : br3=d2mp_l*rold_l*un_l+d3mp_l*rold_r*un_r
819 : br4=d2mp_l*rold_l*ut_l+d3mp_l*rold_r*ut_r
820 : jl.jac(3,4)=n_y*(damp_l*br1+am*br3+dpip_l)
821 : jl.jac(3,4)=jl.jac(3,4)+t_y*damp_l*br2+am*t_y*br4
822 : c-------------------
823 : do 76 i=1,(nsp-1)
824 : br3=d2my_l.vv(i)*rold_l*un_l+d3my_l.vv(i)*rold_r*un_r
825 : br4=d2my_l.vv(i)*rold_l*ut_l+d3my_l.vv(i)*rold_r*ut_r
826 : jl.jac(3,4+i)=n_y*(damg_l*dgdyl.vv(i)*br1+
827 : & am*br3+dpiy_l.vv(i))
828 : jl.jac(3,4+i)=jl.jac(3,4+i)+
829 : & t_y*damg_l*dgdyl.vv(i)*br2+am*t_y*br4
830 : 76 continue
831 : c-------------------------------------------------------------
832 : br3=d2mr_r*rold_l*un_l+mmin_m*un_r+d3mr_r*rold_r*un_r
833 : br4=d2mr_r*rold_l*ut_l+mmin_m*ut_r+d3mr_r*rold_r*ut_r
834 : jr.jac(3,1)=n_y*(damr_r*br1+am*br3+dpir_r)
835 : jr.jac(3,1)=jr.jac(3,1)+t_y*damr_r*br2+am*t_y*br4
836 : c-------------------
837 : br3=d2mu_r*rold_l*un_l+mmin_m*rold_r*t_y/det
838 : br3=br3+d3mu_r*rold_r*un_r
839 : br4=d2mu_r*rold_l*ut_l+mmin_m*rold_r*(-n_y/det)
840 : br4=br4+d3mu_r*rold_r*ut_r
841 : jr.jac(3,2)=n_y*(damu_r*br1+am*br3+dpiu_r)
842 : jr.jac(3,2)=jr.jac(3,2)+t_y*damu_r*br2+am*t_y*br4
843 : c-------------------
844 : br3=d2mv_r*rold_l*un_l+mmin_m*rold_r*(-t_x/det)
845 : br3=br3+d3mv_r*rold_r*un_r
846 : br4=d2mv_r*rold_l*ut_l+mmin_m*rold_r*n_x/det
847 : br4=br4+d3mv_r*rold_r*ut_r
848 : jr.jac(3,3)=n_y*(damv_r*br1+am*br3+dpiv_r)
849 : jr.jac(3,3)=jr.jac(3,3)+t_y*damv_r*br2+am*t_y*br4
850 : c-------------------
851 : br3=d2mp_r*rold_l*un_l+d3mp_r*rold_r*un_r
852 : br4=d2mp_r*rold_l*ut_l+d3mp_r*rold_r*ut_r
853 : jr.jac(3,4)=n_y*(damp_r*br1+am*br3+dpip_r)
854 : jr.jac(3,4)=jr.jac(3,4)+t_y*damp_r*br2+am*t_y*br4
855 : c-------------------
856 : do 77 i=1,(nsp-1)
857 : br3=d2my_r.vv(i)*rold_l*un_l+d3my_r.vv(i)*rold_r*un_r
858 : br4=d2my_r.vv(i)*rold_l*ut_l+d3my_r.vv(i)*rold_r*ut_r
859 : jr.jac(3,4+i)=n_y*(damg_r*dgdyr.vv(i)*br1+
860 : & am*br3+dpiy_r.vv(i))
861 : jr.jac(3,4+i)=jr.jac(3,4+i)+
862 : & t_y*damg_r*dgdyr.vv(i)*br2+am*t_y*br4
863 : 77 continue
864 : c-------------------------------------------------------------
865 : c --------------------- f4 ----------------------------
866 : c-------------------------------------------------------------
867 : hr_l=rold_l*(uold_l*uold_l+vold_l*vold_l)/2.0d0+gal*pold_l/gm1l
868 : hr_r=rold_r*(uold_r*uold_r+vold_r*vold_r)/2.0d0+gar*pold_r/gm1r
869 : f.fu(4)=am*(mpl_m*hr_l+mmin_m*hr_r)
870 : c---------------------
871 : br1=d2mr_l*hr_l+mpl_m*(uold_l*uold_l+vold_l*vold_l)/2.0d0
872 : br1=br1+d3mr_l*hr_r
873 : jl.jac(4,1)=damr_l*f.fu(4)/am+am*br1
874 : c---------------------
875 : br1=d2mu_l*hr_l+mpl_m*uold_l*rold_l
876 : br1=br1+d3mu_l*hr_r
877 : jl.jac(4,2)=damu_l*f.fu(4)/am+am*br1
878 : c---------------------
879 : br1=d2mv_l*hr_l+mpl_m*vold_l*rold_l
880 : br1=br1+d3mv_l*hr_r
881 : jl.jac(4,3)=damv_l*f.fu(4)/am+am*br1
882 : c---------------------
883 : br1=d2mp_l*hr_l+mpl_m*gal/gm1l
884 : br1=br1+d3mp_l*hr_r
885 : jl.jac(4,4)=damp_l*f.fu(4)/am+am*br1
886 : c---------------------
887 : do 78 i=1,(nsp-1)
888 : br1=d2my_l.vv(i)*hr_l+mpl_m*(-pold_l)*dgdyl.vv(i)/(gm1l*gm1l)
889 : br1=br1+d3my_l.vv(i)*hr_r
890 : jl.jac(4,4+i)=damg_l*dgdyl.vv(i)*f.fu(4)/am+am*br1
891 : 78 continue
892 : c----------------------------------------------------------
893 : c----------------------------------------------------------
894 : br1=d2mr_r*hr_l+mmin_m*(uold_r*uold_r+vold_r*vold_r)/2.0d0
895 : br1=br1+d3mr_r*hr_r
896 : jr.jac(4,1)=damr_r*f.fu(4)/am+am*br1
897 : c---------------------
898 : br1=d2mu_r*hr_l+mmin_m*uold_r*rold_r
899 : br1=br1+d3mu_r*hr_r
900 : jr.jac(4,2)=damu_r*f.fu(4)/am+am*br1
901 : c---------------------
902 : br1=d2mv_r*hr_l+mmin_m*vold_r*rold_r
903 : br1=br1+d3mv_r*hr_r
904 : jr.jac(4,3)=damv_r*f.fu(4)/am+am*br1
905 : c---------------------
906 : br1=d2mp_r*hr_l+mmin_m*gar/gm1r
907 : br1=br1+d3mp_r*hr_r
908 : jr.jac(4,4)=damp_r*f.fu(4)/am+am*br1
909 : c---------------------
910 : do 79 i=1,(nsp-1)
911 : br1=d2my_r.vv(i)*hr_l+mmin_m*
912 : & (-pold_r)*dgdyr.vv(i)/(gm1r*gm1r)
913 : br1=br1+d3my_r.vv(i)*hr_r
914 : jr.jac(4,4+i)=damg_r*dgdyr.vv(i)*f.fu(4)/am+am*br1
915 : 79 continue
916 : c-------------------------------------------------------------
917 : c ------------------ f5++ ------------------------------------
918 : c-------------------------------------------------------------
919 : do 80 i=1,(nsp-1)
920 : f.fu(4+i)=am*(mpl_m*rold_l*yl.yet(i)+mmin_m*rold_r*yr.yet(i))
921 : c---------------------
922 : jl.jac(4+i,1)=damr_l*f.fu(4+i)/am+
923 : & am*(d2mr_l*rold_l*yl.yet(i)+mpl_m*yl.yet(i))
924 : jl.jac(4+i,1)=jl.jac(4+i,1)+am*d3mr_l*rold_r*yr.yet(i)
925 : jl.jac(4+i,2)=damu_l*f.fu(4+i)/am+am*(d2mu_l*rold_l*yl.yet(i)+
926 : & d3mu_l*rold_r*yr.yet(i))
927 : jl.jac(4+i,3)=damv_l*f.fu(4+i)/am+am*(d2mv_l*rold_l*yl.yet(i)+
928 : & d3mv_l*rold_r*yr.yet(i))
929 : jl.jac(4+i,4)=damp_l*f.fu(4+i)/am+am*(d2mp_l*rold_l*yl.yet(i)+
930 : & d3mp_l*rold_r*yr.yet(i))
931 : do 81 j=5,(4+nsp-1)
932 : if((4+i).eq.j) then
933 : jl.jac(4+i,j)=damg_l*dgdyl.vv(j-4)*f.fu(4+i)/am+
934 : & am*(d2my_l.vv(j-4)*rold_l*yl.yet(i)+mpl_m*rold_l+
935 : & d3my_l.vv(j-4)*rold_r*yr.yet(i))
936 : else
937 : jl.jac(4+i,j)=damg_l*dgdyl.vv(j-4)*f.fu(4+i)/am+
938 : & am*(d2my_l.vv(j-4)*rold_l*yl.yet(i)+
939 : & d3my_l.vv(j-4)*rold_r*yr.yet(i))
940 : endif
941 : 81 continue
942 : c------------------------------------
943 : jr.jac(4+i,1)=damr_r*f.fu(4+i)/am+
944 : & am*(d2mr_r*rold_l*yl.yet(i)+mmin_m*yr.yet(i))
945 : jr.jac(4+i,1)=jr.jac(4+i,1)+am*d3mr_r*rold_r*yr.yet(i)
946 : jr.jac(4+i,2)=damu_r*f.fu(4+i)/am+am*(d2mu_r*rold_l*yl.yet(i)+
947 : & d3mu_r*rold_r*yr.yet(i))
948 : jr.jac(4+i,3)=damv_r*f.fu(4+i)/am+am*(d2mv_r*rold_l*yl.yet(i)+
949 : & d3mv_r*rold_r*yr.yet(i))
950 : jr.jac(4+i,4)=damp_r*f.fu(4+i)/am+am*(d2mp_r*rold_l*yl.yet(i)+
951 : & d3mp_r*rold_r*yr.yet(i))
952 : do 82 j=1,(nsp-1)
953 : if(i.eq.j) then
954 : jr.jac(4+i,4+j)=damg_r*dgdyr.vv(j)*f.fu(4+i)/am+
955 : & am*(d2my_r.vv(j)*rold_l*yl.yet(i)+mmin_m*rold_r+
956 : & d3my_r.vv(j)*rold_r*yr.yet(i))
957 : else
958 : jr.jac(4+i,4+j)=damg_r*dgdyr.vv(j)*f.fu(4+i)/am+
959 : & am*(d2my_r.vv(j)*rold_l*yl.yet(i)+
960 : & d3my_r.vv(j)*rold_r*yr.yet(i))
961 : endif
962 : 82 continue
963 : 80 continue
964 : c-------------------------------------------------------------
965 : c matrix wl(i,j) represents the derivative of the i-component
966 : c of the vector of primitive variables of the left state with
967 : c respect to the j-component of the vector of the conservative
968 : c variables of the left state.
969 : c
970 : c Here: (rho, ux, uy, p, Y_1,...,Y_(nsp-1)) -
971 : c vector of primitive variables;
972 : c (rho, rho ux, rho uy, rho e, rho Y_1,..., rho Y_(nsp-1)) -
973 : c vector of conservative variables.
974 : c-------------------------------------------------------------
975 : wl.jac(1,1)=1.0d0
976 : wl.jac(1,2)=0.0d0
977 : wl.jac(1,3)=0.0d0
978 : wl.jac(1,4)=0.0d0
979 : do 83 i=1,(nsp-1)
980 : wl.jac(1,4+i)=0.0d0
981 : 83 continue
982 : c------------------------------
983 : wl.jac(2,1)=-uold_l/rold_l
984 : wl.jac(2,2)=1.0d0/rold_l
985 : wl.jac(2,3)=0.0d0
986 : wl.jac(2,4)=0.0d0
987 : do 84 i=1,(nsp-1)
988 : wl.jac(2,4+i)=0.0d0
989 : 84 continue
990 : c------------------------------
991 : wl.jac(3,1)=-vold_l/rold_l
992 : wl.jac(3,2)=0.0d0
993 : wl.jac(3,3)=1.0d0/rold_l
994 : wl.jac(3,4)=0.0d0
995 : do 85 i=1,(nsp-1)
996 : wl.jac(3,4+i)=0.0d0
997 : 85 continue
998 : c------------------------------
999 : br1=0.0d0
1000 : do 86 i=1,(nsp-1)
1001 : br1=br1+dgdyl.vv(i)*yl.yet(i)
1002 : 86 continue
1003 : br1=br1*pold_l/(rold_l*gm1l)
1004 : wl.jac(4,1)=gm1l*(uold_l*uold_l+vold_l*vold_l)/2.0d0-br1
1005 : wl.jac(4,2)=-uold_l*gm1l
1006 : wl.jac(4,3)=-vold_l*gm1l
1007 : wl.jac(4,4)=gm1l
1008 : do 87 i=1,(nsp-1)
1009 : wl.jac(4,4+i)=dgdyl.vv(i)*pold_l/(rold_l*gm1l)
1010 : 87 continue
1011 : c------------------------------
1012 : do 88 i=1,(nsp-1)
1013 : do 89 j=1,4
1014 : wl.jac(4+i,j)=0.0d0
1015 : if(j.eq.1) wl.jac(4+i,j)=-yl.yet(i)/rold_l
1016 : 89 continue
1017 : c------------
1018 : do 890 j=5,(4+nsp-1)
1019 : wl.jac(4+i,j)=0.0d0
1020 : if(4+i.eq.j) then
1021 : wl.jac(4+i,j)=1.0d0/rold_l
1022 : endif
1023 : 890 continue
1024 : 88 continue
1025 : c------------------------------
1026 : c------------------------------
1027 : wr.jac(1,1)=1.0d0
1028 : wr.jac(1,2)=0.0d0
1029 : wr.jac(1,3)=0.0d0
1030 : wr.jac(1,4)=0.0d0
1031 : do 90 i=1,(nsp-1)
1032 : wr.jac(1,4+i)=0.0d0
1033 : 90 continue
1034 : c------------------------------
1035 : wr.jac(2,1)=-uold_r/rold_r
1036 : wr.jac(2,2)=1.0d0/rold_r
1037 : wr.jac(2,3)=0.0d0
1038 : wr.jac(2,4)=0.0d0
1039 : do 91 i=1,(nsp-1)
1040 : wr.jac(2,4+i)=0.0d0
1041 : 91 continue
1042 : c------------------------------
1043 : wr.jac(3,1)=-vold_r/rold_r
1044 : wr.jac(3,2)=0.0d0
1045 : wr.jac(3,3)=1.0d0/rold_r
1046 : wr.jac(3,4)=0.0d0
1047 : do 92 i=1,(nsp-1)
1048 : wr.jac(3,4+i)=0.0d0
1049 : 92 continue
1050 : c------------------------------
1051 : br1=0.0d0
1052 : do 860 i=1,(nsp-1)
1053 : br1=br1+dgdyr.vv(i)*yr.yet(i)
1054 : 860 continue
1055 : br1=br1*pold_r/(rold_r*gm1r)
1056 : wr.jac(4,1)=gm1r*(uold_r*uold_r+vold_r*vold_r)/2.0d0-br1
1057 : wr.jac(4,2)=-uold_r*gm1r
1058 : wr.jac(4,3)=-vold_r*gm1r
1059 : wr.jac(4,4)=gm1r
1060 : do 93 i=1,(nsp-1)
1061 : wr.jac(4,4+i)=dgdyr.vv(i)*pold_r/(rold_r*gm1r)
1062 : 93 continue
1063 : c------------------------------
1064 : do 94 i=1,(nsp-1)
1065 : do 95 j=1,4
1066 : wr.jac(4+i,j)=0.0d0
1067 : if(j.eq.1) wr.jac(4+i,j)=-yr.yet(i)/rold_r
1068 : 95 continue
1069 : c----------------------
1070 : do 950 j=5,(3+nsp)
1071 : wr.jac(4+i,j)=0.0d0
1072 : if(4+i.eq.j) wr.jac(4+i,j)=1.0d0/rold_r
1073 : 950 continue
1074 : 94 continue
1075 : c----------------------------------------------
1076 : c----------------------------------------------
1077 : do 1 i=1,(3+nsp)
1078 : do 2 j=1,(3+nsp)
1079 : jtl.jac(i,j)=0.0d0
1080 : jtr.jac(i,j)=0.0d0
1081 : do 3 k=1,(3+nsp)
1082 : jtl.jac(i,j)=jtl.jac(i,j)+jl.jac(i,k)*wl.jac(k,j)
1083 : jtr.jac(i,j)=jtr.jac(i,j)+jr.jac(i,k)*wr.jac(k,j)
1084 : 3 continue
1085 : 2 continue
1086 : 1 continue
1087 : c----------------------------------------------------------------------
1088 : SEGSUP DMLY_L, DMLY_R,
1089 : & dmry_l, dmry_r,
1090 : & dMpy_l, dMpy_r,
1091 : & dMmy_l, dMmy_r,
1092 : & dmiy_l, dmiy_r,
1093 : & d3my_l, d3my_r,
1094 : & d2my_l, d2my_r,
1095 : & dPpy_l, dPpy_r,
1096 : & dPmy_l, dPmy_r,
1097 : & dpiy_l, dpiy_r,
1098 : & dgdyl, dgdyr
1099 : C--------------------------------------------
1100 : SEGSUP f
1101 : C--------------------------------------------
1102 : SEGSUP JL
1103 : SEGSUP JR
1104 : SEGSUP WL
1105 : SEGSUP WR
1106 : C--------------------------------------------
1107 : jll=jtl
1108 : SEGDES JTL
1109 : jrr=jtr
1110 : SEGDES JTR
1111 : SEGDES YL
1112 : SEGDES YR
1113 : SEGDES CP
1114 : SEGDES CV
1115 : SEGDES MLRECP, MLRECV
1116 : C--------------------------------------------
1117 : return
1118 : end
1119 :
1120 :
1121 :
1122 :
1123 :
1124 :
1125 :
1126 :
1127 :
1128 :
1129 :
1130 :
1131 :
1132 :
1133 :
1134 :
1135 :
1136 :
1137 :
1138 :
1139 :
1140 :
1141 :
1142 :
1143 :
1144 :
1145 :
1146 :
1147 :
1148 :
1149 :
1150 :
1151 :
1152 :
© Cast3M 2003 - All rights reserved.
Disclaimer