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