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