1 : C JBMMW2 SOURCE CHAT 05/01/13 00:49:08 5004 2 : SUBROUTINE JBMMW2(NSP,JLL,WVEC_L,WVEC_R,NVECT,TVECT, 3 : & mpyn,lrecp,lrecv,nlcg,v_inf) 4 : C************************************************************************ 5 : C 6 : C PROJET : CASTEM 2000 7 : C 8 : C NOM : JACBM 9 : C 10 : C DESCRIPTION : Voir KONJA6 11 : C 12 : C LANGAGE : FORTRAN 77 + ESOPE 2000 (avec estensions CISI) 13 : C 14 : C AUTEUR : S. KUDRIAKOV, DM2S/SFME/LTMF 15 : C 16 : C************************************************************************ 17 : C 18 : c---------------------------------------------------------------------- 19 : c GENERAL DESCRIPTION: 20 : c This subroutine provides the jacobians which are the derivatives 21 : c of the numerical flux function defined at the cell interface 22 : c with respect to the conservative variables of the left and right 23 : c cells (relative to the cell interface). 24 : c The low-mach number corrections are made for the flux functions 25 : c 26 : c EQUATIONS: 2D Euler equations of gas dynamics 27 : c 28 : c 29 : c REFERENCE: 1) JCP, 129, 364-382 (1996) 30 : c " A Sequel to AUSM: AUSM+ ". 31 : c M.S.Liou 32 : c 2) AIAA Journal, Sept. 1998 33 : c " Low-Diffusion Flux-Splitting Methods for Flows at All Speeds" 34 : c J.R.Edwards and M.S.Liou 35 : c---------------------------------------------------------------------- 36 : c INPUT: 37 : c 38 : c alpha -- parameter of the AUSM+ scheme in the Pressure function; 39 : c ( -3/4 <= alpha <= 3/16 ) (imposed as a parameter) 40 : c 41 : c beta -- parameter of the AUSM+ scheme in the Mach function; 42 : c ( -1/16 <= beta <= 1/2 ) (imposed as a parameter) 43 : c 44 : c wvec_l -- vector of the primitive variables (rho,ux,uy,p) at the 45 : c left cell; 46 : c 47 : c wvec_r -- vector of the primitive variables (rho,ux,uy,p) at the 48 : c right cell; 49 : c 50 : c nvect -- normal vector to the interface (2 components in 2D); 51 : c 52 : c tvect -- tangential vector to the interface; 53 : c 54 : c ga -- ratio of the specific heats (assumed constant) 55 : c 56 : c v_inf -- parameter for reference velocity, see Ref. 2). 57 : c---------------------------------------------------------------------- 58 : c 59 : c OUTPUT: 60 : c 61 : c jtl -- jakobian matrix 4 by 4 - derivatives of the numerical 62 : c flux function with respect to the conservative variables 63 : c from the left cell; 64 : c 65 : c jtr -- jakobian matrix 4 by 4 - derivatives of the numerical 66 : c flux function with respect to the conservative variables 67 : c from the right cell. 68 : c---------------------------------------------------------------------- 69 : IMPLICIT INTEGER(I-N) 70 : real*8 wvec_l(4),wvec_r(4) 71 : real*8 nvect(2),tvect(2) 72 : real*8 alpha,beta,bcoef,tcoef 73 : real*8 gal,gar,gm1l,gm1r,temph 74 : real*8 n_x,n_y,t_x,t_y 75 : real*8 un_l, un_r, ut_l, ut_r 76 : real*8 ml,mr,Mplus,Mmin,mmid 77 : real*8 mpl_m, mmin_m,am 78 : real*8 rold_l,uold_l,vold_l,pold_l,eold_l 79 : real*8 rold_r,uold_r,vold_r,pold_r,eold_r 80 : real*8 Pplus,Pmin 81 : real*8 top,bot,bots 82 : real*8 br1,br3,br4,temp_l,temp_r,brac_l,brac_r 83 : real*8 aleft, arigh 84 : real*8 damr_l,damr_r,damu_l,damu_r 85 : real*8 damv_l,damv_r,damp_l,damp_r 86 : real*8 damg_l,damg_r 87 : real*8 dmlr_l,dmlr_r,dmlu_l,dmlu_r 88 : real*8 dmlv_l,dmlv_r,dmlp_l,dmlp_r 89 : real*8 dmrr_l,dmrr_r,dmru_l,dmru_r 90 : real*8 dmrv_l,dmrv_r,dmrp_l,dmrp_r 91 : real*8 dMpr_l,dMpr_r,dMpu_l,dMpu_r 92 : real*8 dMpv_l,dMpv_r,dMpp_l,dMpp_r 93 : real*8 dMmr_l,dMmr_r,dMmu_l,dMmu_r 94 : real*8 dMmv_l,dMmv_r,dMmp_l,dMmp_r 95 : real*8 dmir_l,dmir_r,dmiu_l,dmiu_r 96 : real*8 dmiv_l,dmiv_r,dmip_l,dmip_r 97 : real*8 d3mr_l,d3mr_r,d3mu_l,d3mu_r 98 : real*8 d3mv_l,d3mv_r,d3mp_l,d3mp_r 99 : real*8 d2mr_l,d2mr_r,d2mu_l,d2mu_r 100 : real*8 d2mv_l,d2mv_r,d2mp_l,d2mp_r 101 : real*8 dPpr_l,dPpr_r,dPpu_l,dPpu_r 102 : real*8 dPpv_l,dPpv_r,dPpp_l,dPpp_r 103 : real*8 dPmr_l,dPmr_r,dPmu_l,dPmu_r 104 : real*8 dPmv_l,dPmv_r,dPmp_l,dPmp_r 105 : real*8 dpir_l,dpir_r,dpiu_l,dpiu_r 106 : real*8 dpiv_l,dpiv_r,dpip_l,dpip_r 107 : real*8 epsil,qq,amw,Mmin1,Mplus1 108 : real*8 fmid,mlw,mrw,termp 109 : real*8 ur_r,ur_l,urm,mhalf,mhalfr 110 : real*8 durr_l,durr_r,duru_l,duru_r 111 : real*8 durv_l,durv_r,durp_l,durp_r 112 : real*8 dmhr_l,dmhr_r,dmhu_l,dmhu_r 113 : real*8 dmhv_l,dmhv_r,dmhp_l,dmhp_r 114 : real*8 dmfr_l,dmfr_r,dmfu_l,dmfu_r 115 : real*8 dmfv_l,dmfv_r,dmfp_l,dmfp_r 116 : real*8 dfm_mf,dfm_mh 117 : real*8 dfmr_l,dfmr_r,dfmu_l,dfmu_r 118 : real*8 dfmv_l,dfmv_r,dfmp_l,dfmp_r 119 : real*8 m1mr_l,m1mr_r,m1mu_l,m1mu_r 120 : real*8 m1mv_l,m1mv_r,m1mp_l,m1mp_r 121 : real*8 m1pr_l,m1pr_r,m1pu_l,m1pu_r 122 : real*8 m1pv_l,m1pv_r,m1pp_l,m1pp_r 123 : real*8 tmpr_l,tmpr_r,tmpu_l,tmpu_r 124 : real*8 tmpv_l,tmpv_r,tmpp_l,tmpp_r 125 : real*8 coef,canc,v_inf 126 : real*8 sr_l,sr_r,su_l,su_r,sv_l,sv_r,sp_l,sp_r 127 : real*8 rum,rumr_l,rumu_l,rumv_l,rump_l 128 : real*8 rumr_r,rumu_r,rumv_r,rump_r 129 : real*8 dfmg_l,dfmg_r,dmfg_l,dmfg_r 130 : real*8 dmhg_l,dmhg_r,durg_l,durg_r 131 : integer i,j,k,jll,lrecp,lrecv,nlcd,nlcg,nsp 132 : parameter(alpha = 0.1875D0, beta = 0.125D0) 133 : parameter(epsil = 1.0d0) 134 : canc=1.0d0 135 : c------------------------------------------------------------- 136 : -INC SMCHPOI 137 : POINTEUR MPYN.MPOVAL 138 : C------------------------------------------------------------- 139 : -INC SMLREEL 140 : POINTEUR MLRECP.MLREEL, MLRECV.MLREEL 141 : C------------------------------------------------------------- 142 : C******* Les fractionines massiques ************************** 143 : C------------------------------------------------------------- 144 : SEGMENT FRAMAS 145 : REAL*8 YET(NSP) 146 : ENDSEGMENT 147 : POINTEUR YL.FRAMAS, YR.FRAMAS 148 : C------------------------------------------------------- 149 : C********** Les CP's and CV's *********************** 150 : C------------------------------------------------------- 151 : SEGMENT GCONST 152 : REAL*8 GC(NSP) 153 : ENDSEGMENT 154 : POINTEUR CP.GCONST, CV.GCONST 155 : C------------------------------------------------------------- 156 : C******** Segments for the elementary matrixes ************* 157 : C------------------------------------------------------------- 158 : SEGMENT JACEL 159 : REAL*8 JAC(3+NSP,3+NSP) 160 : ENDSEGMENT 161 : POINTEUR JTL.JACEL, JTR.JACEL, JL.JACEL, JR.JACEL, 162 : & WL.JACEL, WR.JACEL 163 : c---------------------------------------- 164 : SEGINI JTL 165 : SEGINI JTR 166 : SEGINI JL 167 : SEGINI JR 168 : SEGINI WL 169 : SEGINI WR 170 : C------------------------------------------------------------- 171 : C********** Segments for the vectors *********************** 172 : C------------------------------------------------------------- 173 : SEGMENT VECEL 174 : REAL*8 VV(NSP) 175 : ENDSEGMENT 176 : POINTEUR DMLY_L.VECEL, DMLY_R.VECEL, 177 : & dmry_l.vecel, dmry_r.vecel, 178 : & dMpy_l.vecel, dMpy_r.vecel, 179 : & dMmy_l.vecel, dMmy_r.vecel, 180 : & dmiy_l.vecel, dmiy_r.vecel, 181 : & d3my_l.vecel, d3my_r.vecel, 182 : & d2my_l.vecel, d2my_r.vecel, 183 : & dPpy_l.vecel, dPpy_r.vecel, 184 : & dPmy_l.vecel, dPmy_r.vecel, 185 : & dpiy_l.vecel, dpiy_r.vecel, 186 : & dgdyl.vecel, dgdyr.vecel, 187 : & sy_l.vecel, sy_r.vecel, 188 : & m1py_l.vecel, m1py_r.vecel, 189 : & m1my_l.vecel, m1my_r.vecel, 190 : & tmpy_l.vecel, tmpy_r.vecel, 191 : & rumy_l.vecel, rumy_r.vecel 192 : C---------------------------------------------- 193 : SEGINI DMLY_L, DMLY_R, 194 : & dmry_l, dmry_r, 195 : & dMpy_l, dMpy_r, 196 : & dMmy_l, dMmy_r, 197 : & dmiy_l, dmiy_r, 198 : & d3my_l, d3my_r, 199 : & d2my_l, d2my_r, 200 : & dPpy_l, dPpy_r, 201 : & dPmy_l, dPmy_r, 202 : & dpiy_l, dpiy_r, 203 : & dgdyl, dgdyr, 204 : & sy_l, sy_r, 205 : & m1py_l, m1py_r, 206 : & m1my_l, m1my_r, 207 : & tmpy_l, tmpy_r, 208 : & rumy_l, rumy_r 209 : C------------------------------------------------------------- 210 : C********** Segments for the flux-vector ******************* 211 : C------------------------------------------------------------- 212 : SEGMENT FUNEL 213 : REAL*8 FU(3+NSP) 214 : ENDSEGMENT 215 : POINTEUR f.funel 216 : C------------------------- 217 : SEGINI f 218 : C------------------------------------------------------------ 219 : SEGINI YL, YR 220 : SEGACT MPYN 221 : DO 100 I=1,(NSP-1) 222 : YL.YET(I)=MPYN.VPOCHA(NLCG,I) 223 : YR.YET(I)=MPYN.VPOCHA(NLCG,I) 224 : 100 CONTINUE 225 : C---------------------------------------- 226 : SEGINI CP, CV 227 : MLRECP = LRECP 228 : MLRECV = LRECV 229 : SEGACT MLRECP, MLRECV 230 : DO 101 I=1,(NSP-1) 231 : CP.GC(I)=MLRECP.PROG(I) 232 : CV.GC(I)=MLRECV.PROG(I) 233 : 101 CONTINUE 234 : CP.GC(NSP)=MLRECP.PROG(NSP) 235 : CV.GC(NSP)=MLRECV.PROG(NSP) 236 : c------------------------------------------------------------- 237 : c Computing GAMMA at the left cell and its derivatives 238 : c with respect to the primitive variables Y_i 239 : c------------------------------------------------------------- 240 : top=0.0D0 241 : bot=0.0D0 242 : do 40 i=1,(nsp-1) 243 : top=top+yl.yet(i)*(cp.gc(i)-cp.gc(nsp)) 244 : bot=bot+yl.yet(i)*(cv.gc(i)-cv.gc(nsp)) 245 : 40 continue 246 : top=cp.gc(nsp)+top 247 : bot=cv.gc(nsp)+bot 248 : gal=top/bot 249 : gm1l=gal-1.0d0 250 : c------------------------------------------------------------- 251 : do 41 i=1,(nsp-1) 252 : dgdyl.vv(i)=(cp.gc(i)-cp.gc(nsp)- 253 : & gal*(cv.gc(i)-cv.gc(nsp)))/bot 254 : 41 continue 255 : c------------------------------------------------------------- 256 : c Computing GAMMA at the right cell and its derivatives 257 : c with respect to the primitive variables Y_i 258 : c------------------------------------------------------------- 259 : top=0.0D0 260 : bot=0.0D0 261 : do 42 i=1,(nsp-1) 262 : top=top+yr.yet(i)*(cp.gc(i)-cp.gc(nsp)) 263 : bot=bot+yr.yet(i)*(cv.gc(i)-cv.gc(nsp)) 264 : 42 continue 265 : top=cp.gc(nsp)+top 266 : bot=cv.gc(nsp)+bot 267 : gar=top/bot 268 : gm1r=gar-1.0d0 269 : c------------------------------------------------------------- 270 : do 43 i=1,(nsp-1) 271 : dgdyr.vv(i)=(cp.gc(i)-cp.gc(nsp)- 272 : & gar*(cv.gc(i)-cv.gc(nsp)))/bot 273 : 43 continue 274 : c------------------------------------------------------------- 275 : n_x=nvect(1) 276 : n_y=nvect(2) 277 : t_x=tvect(1) 278 : t_y=tvect(2) 279 : c---------------------------- 280 : rold_l=wvec_l(1) 281 : uold_l=wvec_l(2) 282 : vold_l=wvec_l(3) 283 : pold_l=wvec_l(4) 284 : c----------------------- 285 : rold_r=wvec_r(1) 286 : uold_r=wvec_r(2) 287 : vold_r=wvec_r(3) 288 : pold_r=wvec_r(4) 289 : c------------------------------------------------------------------ 290 : c Computation of the specific total energy on the left and right. 291 : c------------------------------------------------------------------ 292 : eold_l=(uold_l*uold_l+vold_l*vold_l)/2.0d0 293 : eold_l=eold_l+pold_l/(gm1l*rold_l) 294 : eold_r=(uold_r*uold_r+vold_r*vold_r)/2.0d0 295 : eold_r=eold_r+pold_r/(gm1r*rold_r) 296 : c------------------------------------------------------------------ 297 : c Computing reference velocity and its derivatives 298 : c------------------------------------------------------------------ 299 : aleft=sqrt(gal*pold_l/rold_l) 300 : arigh=sqrt(gar*pold_r/rold_r) 301 : qq=sqrt(uold_l*uold_l+vold_l*vold_l) 302 : if(qq .lt. (epsil*v_inf)) then 303 : ur_l = epsil*v_inf 304 : durr_l=0.0d0 305 : duru_l=0.0d0 306 : durv_l=0.0d0 307 : durp_l=0.0d0 308 : durg_l=0.0d0 309 : else 310 : ur_l=qq 311 : durr_l=0.0d0 312 : duru_l=uold_l/qq 313 : durv_l=vold_l/qq 314 : durp_l=0.0d0 315 : durg_l=0.0d0 316 : endif 317 : c------------------------------ 318 : if(ur_l .ge. aleft) then 319 : ur_l=aleft 320 : durr_l=-aleft/(2.0d0*rold_l) 321 : duru_l=0.0d0 322 : durv_l=0.0d0 323 : durp_l=gal/(2.0d0*aleft*rold_l) 324 : durg_l=aleft/(2.0d0*gal) 325 : endif 326 : c------------------------------------------------------------------ 327 : qq=sqrt(uold_r*uold_r+vold_r*vold_r) 328 : if(qq .lt. (epsil*v_inf)) then 329 : ur_r = epsil*v_inf 330 : durr_r=0.0d0 331 : duru_r=0.0d0 332 : durv_r=0.0d0 333 : durp_r=0.0d0 334 : durg_r=0.0d0 335 : else 336 : ur_r=qq 337 : durr_r=0.0d0 338 : duru_r=uold_r/qq 339 : durv_r=vold_r/qq 340 : durp_r=0.0d0 341 : durg_r=0.0d0 342 : endif 343 : c--------------------------- 344 : if(ur_r .ge. arigh) then 345 : ur_r=arigh 346 : durr_r=-arigh/(2.0d0*rold_r) 347 : duru_r=0.0d0 348 : durv_r=0.0d0 349 : durp_r=gar/(2.0d0*arigh*rold_r) 350 : durg_r=arigh/(2.0d0*gar) 351 : endif 352 : c------------------------------------------------------------------ 353 : c------------------------------------------------------------------- 354 : c Reference velocity at the interface is taken as an average 355 : c of the reference velocities of the neighbouring cells 356 : c------------------------------------------------------------------- 357 : urm=0.5d0*(ur_l+ur_r) 358 : durr_l=0.5d0*durr_l 359 : duru_l=0.5d0*duru_l 360 : durv_l=0.5d0*durv_l 361 : durp_l=0.5d0*durp_l 362 : durg_l=0.5d0*durg_l 363 : c------------------------- 364 : durr_r=0.5d0*durr_r 365 : duru_r=0.5d0*duru_r 366 : durv_r=0.5d0*durv_r 367 : durp_r=0.5d0*durp_r 368 : durg_r=0.5d0*durg_r 369 : c------------------------------------------------------------------- 370 : c Computation of the speed of sound and its derivatives; 371 : c numerical speed of sound at the interface is taken as an average 372 : c of the speeds of sounds of the neighbouring cells 373 : c------------------------------------------------------------------- 374 : am=0.5d0*(aleft+arigh) 375 : c------------------------------------------------------------------- 376 : if(abs(urm/am-1.0d0).le.0.000001d0) then 377 : coef=0.0d0 378 : else 379 : coef=1.0d0 380 : endif 381 : c------------------------------------------------------------------- 382 : damr_r=-arigh/(4.0d0*rold_r) 383 : damu_r=0.0d0 384 : damv_r=0.0d0 385 : damp_r=gar/(4.0d0*arigh*rold_r) 386 : damg_r=arigh/(4.0d0*gar) 387 : c----------------------- 388 : damr_l=-aleft/(4.0d0*rold_l) 389 : damu_l=0.0d0 390 : damv_l=0.0d0 391 : damp_l=gal/(4.0d0*aleft*rold_l) 392 : damg_l=aleft/(4.0d0*gal) 393 : c------------------------------------------------------------------- 394 : c Computing numerical Mach number and its derivatives, 395 : c see p.370, under (A1). 396 : c------------------------------------------------------------------- 397 : un_l=uold_l*n_x+vold_l*n_y 398 : un_r=uold_r*n_x+vold_r*n_y 399 : ut_l=uold_l*t_x+vold_l*t_y 400 : ut_r=uold_r*t_x+vold_r*t_y 401 : c------------------------------------------------------------------- 402 : ml=un_l/am 403 : mr=un_r/am 404 : mhalf=0.5d0*(un_l+un_r)/am 405 : c--------------------- 406 : top=0.5d0*(un_l+un_r)/(am*am) 407 : dmhr_l=-top*damr_l 408 : dmhu_l=n_x/2.0d0/am-top*damu_l 409 : dmhv_l=n_y/2.0d0/am-top*damv_l 410 : dmhp_l=-top*damp_l 411 : dmhg_l=-top*damg_l 412 : c--------------------- 413 : dmhr_r=-top*damr_r 414 : dmhu_r=n_x/2.0d0/am-top*damu_r 415 : dmhv_r=n_y/2.0d0/am-top*damv_r 416 : dmhp_r=-top*damp_r 417 : dmhg_r=-top*damg_r 418 : c-------------------------------- 419 : mhalfr=urm/am 420 : c--------------------- 421 : top=urm/(am*am) 422 : dmfr_l=durr_l/am-top*damr_l 423 : dmfu_l=duru_l/am-top*damu_l 424 : dmfv_l=durv_l/am-top*damv_l 425 : dmfp_l=durp_l/am-top*damp_l 426 : dmfg_l=durg_l/am-top*damg_l 427 : c--------------------- 428 : dmfr_r=durr_r/am-top*damr_r 429 : dmfu_r=duru_r/am-top*damu_r 430 : dmfv_r=durv_r/am-top*damv_r 431 : dmfp_r=durp_r/am-top*damp_r 432 : dmfg_r=durg_r/am-top*damg_r 433 : c------------------------------------------------------------------- 434 : c Scaling function for the speed of sound and its derivatives 435 : c------------------------------------------------------------------- 436 : top=(1.0d0-mhalfr*mhalfr)*(1.0d0-mhalfr*mhalfr) 437 : top=top*mhalf*mhalf+4.0d0*mhalfr*mhalfr 438 : bot=1.0d0+mhalfr*mhalfr 439 : if(abs(canc-0.0d0).lt.0.000001d0) then 440 : fmid=1.0d0 441 : else 442 : fmid=sqrt(top)/bot 443 : endif 444 : c-------------------------- 445 : temph=-4.0d0*(1.0d0-mhalfr*mhalfr)*mhalfr 446 : temph=temph*mhalf*mhalf+8.0d0*mhalfr 447 : dfm_mf=temph/(sqrt(top)*2.0d0*bot) 448 : dfm_mf=dfm_mf-sqrt(top)*2.0d0*mhalfr/(bot*bot) 449 : c-------------------------- 450 : temph=2.0d0*(1.0d0-mhalfr*mhalfr)*(1.0d0-mhalfr*mhalfr)*mhalf 451 : dfm_mh=temph/(2.0d0*bot*sqrt(top)) 452 : c-------------------------- 453 : dfmr_l=dfm_mf*dmfr_l+dfm_mh*dmhr_l 454 : dfmu_l=dfm_mf*dmfu_l+dfm_mh*dmhu_l 455 : dfmv_l=dfm_mf*dmfv_l+dfm_mh*dmhv_l 456 : dfmp_l=dfm_mf*dmfp_l+dfm_mh*dmhp_l 457 : dfmg_l=dfm_mf*dmfg_l+dfm_mh*dmhg_l 458 : c-------------------------- 459 : dfmr_r=dfm_mf*dmfr_r+dfm_mh*dmhr_r 460 : dfmu_r=dfm_mf*dmfu_r+dfm_mh*dmhu_r 461 : dfmv_r=dfm_mf*dmfv_r+dfm_mh*dmhv_r 462 : dfmp_r=dfm_mf*dmfp_r+dfm_mh*dmhp_r 463 : dfmg_r=dfm_mf*dmfg_r+dfm_mh*dmhg_r 464 : c-------------------------- 465 : amw=fmid*am 466 : mlw=un_l/amw 467 : mrw=un_r/amw 468 : c-------------------------- 469 : damr_l=canc*coef*dfmr_l*am+fmid*damr_l 470 : damu_l=canc*coef*dfmu_l*am+fmid*damu_l 471 : damv_l=canc*coef*dfmv_l*am+fmid*damv_l 472 : damp_l=canc*coef*dfmp_l*am+fmid*damp_l 473 : damg_l=canc*coef*dfmg_l*am+fmid*damg_l 474 : c-------------------------- 475 : damr_r=canc*coef*dfmr_r*am+fmid*damr_r 476 : damu_r=canc*coef*dfmu_r*am+fmid*damu_r 477 : damv_r=canc*coef*dfmv_r*am+fmid*damv_r 478 : damp_r=canc*coef*dfmp_r*am+fmid*damp_r 479 : damg_r=canc*coef*dfmg_r*am+fmid*damg_r 480 : c-/-/-/-/-/-/-/-/-/-/-/-/-/-/-/-/-/ 481 : am=amw 482 : c------------------------------------------------------------------- 483 : c Redefinition of the numerical mach numbers 484 : c------------------------------------------------------------------- 485 : if(abs(canc-0.0d0).lt.0.000001d0) then 486 : top=2.0d0 487 : bot=0.0d0 488 : else 489 : top=1.0d0+mhalfr*mhalfr 490 : bot=1.0d0-mhalfr*mhalfr 491 : endif 492 : ml=0.5d0*(top*mlw+bot*mrw) 493 : mr=0.5d0*(top*mrw+bot*mlw) 494 : c------------------------------------------------------------------- 494 : c Mplus and Mmin are calligraphic lettes M+ and M- from the paper, 495 : c see (19a) and (19b), p.367. 496 : c------------------------------------------------------------------- 497 : if(abs(ml) .ge. 1.0d0) then 498 : Mplus=(ml+abs(ml))/2.0d0 499 : else 500 : Mplus=(ml+1.0d0)*(ml+1.0d0)/4.0d0 501 : Mplus=Mplus+beta*(ml*ml-1.0d0)*(ml*ml-1.0d0) 502 : endif 503 : Mplus1=(ml+abs(ml))/2.0d0 504 : c------------------------------------------------------------------- 505 : if(abs(mr) .ge. 1.0d0) then 506 : Mmin=(mr-abs(mr))/2.0d0 507 : else 508 : Mmin=-(mr-1.0d0)*(mr-1.0d0)/4.0d0 509 : Mmin=Mmin-beta*(mr*mr-1.0d0)*(mr*mr-1.0d0) 510 : endif 511 : Mmin1=(mr-abs(mr))/2.0d0 512 : c------------------------------------------------------------------- 513 : c Derivatives of ml and mr with respect to both sets of primitive 514 : c variables: left and right. 515 : c------------------------------------------------------------------- 516 : temp_l=-un_l/(am*am) 517 : temp_r=-un_r/(am*am) 518 : c-------- 519 : dmlr_l=temp_l*damr_l 520 : dmlr_r=temp_l*damr_r 521 : dmrr_l=temp_r*damr_l 522 : dmrr_r=temp_r*damr_r 523 : c-------- 524 : dmlu_l=n_x/am+temp_l*damu_l 525 : dmlu_r=temp_l*damu_r 526 : dmru_l=temp_r*damu_l 527 : dmru_r=n_x/am+temp_r*damu_r 528 : c-------- 529 : dmlv_l=n_y/am+temp_l*damv_l 530 : dmlv_r=temp_l*damv_r 531 : dmrv_l=temp_r*damv_l 532 : dmrv_r=n_y/am+temp_r*damv_r 533 : c-------- 534 : dmlp_l=temp_l*damp_l 535 : dmlp_r=temp_l*damp_r 536 : dmrp_l=temp_r*damp_l 537 : dmrp_r=temp_r*damp_r 538 : c-------- 539 : do 44 i=1,(nsp-1) 540 : dmly_l.vv(i)=temp_l*damg_l*dgdyl.vv(i) 541 : dmly_r.vv(i)=temp_l*damg_r*dgdyr.vv(i) 542 : dmry_l.vv(i)=temp_r*damg_l*dgdyl.vv(i) 543 : dmry_r.vv(i)=temp_r*damg_r*dgdyr.vv(i) 544 : 44 continue 545 : c----------------------------- 546 : sr_l=dmlr_l 547 : sr_r=dmlr_r 548 : su_l=dmlu_l 549 : su_r=dmlu_r 550 : sv_l=dmlv_l 551 : sv_r=dmlv_r 552 : sp_l=dmlp_l 553 : sp_r=dmlp_r 554 : do 440 i=1,(nsp-1) 555 : sy_l.vv(i) = dmly_l.vv(i) 556 : sy_r.vv(i) = dmly_r.vv(i) 557 : 440 continue 558 : c----------------------------------------------------------------- 559 : c Redefinition of the derivatives of the ml & mr 560 : c----------------------------------------------------------------- 561 : temp_l=(mlw-mrw)*mhalfr 562 : temp_r=-temp_l 563 : c-------- 564 : dmlr_l=0.5d0*(top*dmlr_l+bot*dmrr_l)+canc*coef*temp_l*dmfr_l 565 : dmlu_l=0.5d0*(top*dmlu_l+bot*dmru_l)+canc*coef*temp_l*dmfu_l 566 : dmlv_l=0.5d0*(top*dmlv_l+bot*dmrv_l)+canc*coef*temp_l*dmfv_l 567 : dmlp_l=0.5d0*(top*dmlp_l+bot*dmrp_l)+canc*coef*temp_l*dmfp_l 568 : c-------- 569 : dmlr_r=0.5d0*(top*dmlr_r+bot*dmrr_r)+canc*coef*temp_l*dmfr_r 570 : dmlu_r=0.5d0*(top*dmlu_r+bot*dmru_r)+canc*coef*temp_l*dmfu_r 571 : dmlv_r=0.5d0*(top*dmlv_r+bot*dmrv_r)+canc*coef*temp_l*dmfv_r 572 : dmlp_r=0.5d0*(top*dmlp_r+bot*dmrp_r)+canc*coef*temp_l*dmfp_r 573 : c-------- 574 : do 441 i=1,(nsp-1) 575 : dmly_l.vv(i)=0.5d0*(top*dmly_l.vv(i)+bot*dmry_l.vv(i))+ 576 : & canc*coef*temp_l*dmfg_l*dgdyl.vv(i) 577 : dmly_r.vv(i)=0.5d0*(top*dmly_r.vv(i)+bot*dmry_r.vv(i))+ 578 : & canc*coef*temp_l*dmfg_r*dgdyr.vv(i) 579 : 441 continue 580 : c-------- 581 : dmrr_l=0.5d0*(top*dmrr_l+bot*sr_l)+canc*coef*temp_r*dmfr_l 582 : dmru_l=0.5d0*(top*dmru_l+bot*su_l)+canc*coef*temp_r*dmfu_l 583 : dmrv_l=0.5d0*(top*dmrv_l+bot*sv_l)+canc*coef*temp_r*dmfv_l 584 : dmrp_l=0.5d0*(top*dmrp_l+bot*sp_l)+canc*coef*temp_r*dmfp_l 585 : c-------- 586 : dmrr_r=0.5d0*(top*dmrr_r+bot*sr_r)+canc*coef*temp_r*dmfr_r 587 : dmru_r=0.5d0*(top*dmru_r+bot*su_r)+canc*coef*temp_r*dmfu_r 588 : dmrv_r=0.5d0*(top*dmrv_r+bot*sv_r)+canc*coef*temp_r*dmfv_r 589 : dmrp_r=0.5d0*(top*dmrp_r+bot*sp_r)+canc*coef*temp_r*dmfp_r 590 : do 442 i=1,(nsp-1) 591 : dmry_l.vv(i)=0.5d0*(top*dmry_l.vv(i)+bot*sy_l.vv(i))+ 592 : & canc*coef*temp_r*dmfg_l*dgdyl.vv(i) 593 : dmry_r.vv(i)=0.5d0*(top*dmry_r.vv(i)+bot*sy_r.vv(i))+ 594 : & canc*coef*temp_r*dmfg_r*dgdyr.vv(i) 595 : 442 continue 596 : c----------------------------------------------------------- 597 : c mmid is m_{1/2} (notation as in the paper, see (13),p.366) 598 : c----------------------------------------------------------- 599 : mmid=Mplus+Mmin 600 : c----------------------------------------------------------- 601 : c Computing the derivatives of M+ and M- 602 : c----------------------------------------------------------- 603 : if(ml .ge. 1.0d0) then 604 : dMpr_l=dmlr_l 605 : dMpu_l=dmlu_l 606 : dMpv_l=dmlv_l 607 : dMpp_l=dmlp_l 608 : do 45 i=1,(nsp-1) 609 : dMpy_l.vv(i)=dmly_l.vv(i) 610 : 45 continue 611 : c-------------------- 612 : dMpr_r=dmlr_r 613 : dMpu_r=dmlu_r 614 : dMpv_r=dmlv_r 615 : dMpp_r=dmlp_r 616 : do 46 i=1,(nsp-1) 617 : dMpy_r.vv(i)=dmly_r.vv(i) 618 : 46 continue 619 : else 620 : if((ml .gt. -1.0d0) .and. (ml .lt. 1.0d0)) then 621 : temph=(ml+1.0d0)/2.0d0 622 : dMpr_l=(temph+4.0d0*beta*ml*(ml*ml-1.0d0))*dmlr_l 623 : dMpu_l=(temph+4.0d0*beta*ml*(ml*ml-1.0d0))*dmlu_l 624 : dMpv_l=(temph+4.0d0*beta*ml*(ml*ml-1.0d0))*dmlv_l 625 : dMpp_l=(temph+4.0d0*beta*ml*(ml*ml-1.0d0))*dmlp_l 626 : do 47 i=1,(nsp-1) 627 : dMpy_l.vv(i)=(temph+4.0d0*beta*ml* 628 : & (ml*ml-1.0d0))*dmly_l.vv(i) 629 : 47 continue 630 : c-------------------- 631 : dMpr_r=(temph+4.0d0*beta*ml*(ml*ml-1.0d0))*dmlr_r 632 : dMpu_r=(temph+4.0d0*beta*ml*(ml*ml-1.0d0))*dmlu_r 633 : dMpv_r=(temph+4.0d0*beta*ml*(ml*ml-1.0d0))*dmlv_r 634 : dMpp_r=(temph+4.0d0*beta*ml*(ml*ml-1.0d0))*dmlp_r 635 : do 48 i=1,(nsp-1) 636 : dMpy_r.vv(i)=(temph+4.0d0*beta*ml* 637 : & (ml*ml-1.0d0))*dmly_r.vv(i) 638 : 48 continue 639 : else 640 : dMpr_l=0.0d0 641 : dMpu_l=0.0d0 642 : dMpv_l=0.0d0 643 : dMpp_l=0.0d0 644 : do 49 i=1,(nsp-1) 645 : dMpy_l.vv(i)=0.0d0 646 : 49 continue 647 : c--------------------- 648 : dMpr_r=0.0d0 649 : dMpu_r=0.0d0 650 : dMpv_r=0.0d0 651 : dMpp_r=0.0d0 652 : do 50 i=1,(nsp-1) 653 : dMpy_r.vv(i)=0.0d0 654 : 50 continue 655 : endif 656 : endif 657 : c----------------------------------------------------------- 658 : c addition of low Mach number 659 : c----------------------------------------------------------- 660 : if(ml .ge. 0.0d0) then 661 : m1pr_l=dmlr_l 662 : m1pu_l=dmlu_l 663 : m1pv_l=dmlv_l 664 : m1pp_l=dmlp_l 665 : do 450 i=1,(nsp-1) 666 : m1py_l.vv(i)=dmly_l.vv(i) 667 : 450 continue 668 : c--------------------- 669 : m1pr_r=dmlr_r 670 : m1pu_r=dmlu_r 671 : m1pv_r=dmlv_r 672 : m1pp_r=dmlp_r 673 : do 460 i=1,(nsp-1) 674 : m1py_r.vv(i)=dmly_r.vv(i) 675 : 460 continue 676 : else 677 : m1pr_l=0.0d0 678 : m1pu_l=0.0d0 679 : m1pv_l=0.0d0 680 : m1pp_l=0.0d0 681 : do 451 i=1,(nsp-1) 682 : m1py_l.vv(i)=0.0d0 683 : 451 continue 684 : c-------------------- 685 : m1pr_r=0.0d0 686 : m1pu_r=0.0d0 687 : m1pv_r=0.0d0 688 : m1pp_r=0.0d0 689 : do 461 i=1,(nsp-1) 690 : m1py_r.vv(i)=0.0d0 691 : 461 continue 692 : endif 693 : c----------------------------------------------------------- 694 : if(mr .ge. 1.0d0) then 695 : dMmr_l=0.0d0 696 : dMmu_l=0.0d0 697 : dMmv_l=0.0d0 698 : dMmp_l=0.0d0 699 : do 51 i=1,(nsp-1) 700 : dMmy_l.vv(i)=0.0d0 701 : 51 continue 702 : c--------------------- 703 : dMmr_r=0.0d0 704 : dMmu_r=0.0d0 705 : dMmv_r=0.0d0 706 : dMmp_r=0.0d0 707 : do 52 i=1,(nsp-1) 708 : dMmy_r.vv(i)=0.0d0 709 : 52 continue 710 : else 711 : if((mr .gt. -1.0d0) .and. (mr .lt. 1.0d0)) then 712 : temph=(-mr+1.0d0)/2.0d0 713 : dMmr_l=(temph-4.0d0*beta*mr*(mr*mr-1.0d0))*dmrr_l 714 : dMmu_l=(temph-4.0d0*beta*mr*(mr*mr-1.0d0))*dmru_l 715 : dMmv_l=(temph-4.0d0*beta*mr*(mr*mr-1.0d0))*dmrv_l 716 : dMmp_l=(temph-4.0d0*beta*mr*(mr*mr-1.0d0))*dmrp_l 717 : do 53 i=1,(nsp-1) 718 : dMmy_l.vv(i)=(temph-4.0d0*beta*mr* 719 : & (mr*mr-1.0d0))*dmry_l.vv(i) 720 : 53 continue 721 : c-------------------- 722 : dMmr_r=(temph-4.0d0*beta*mr*(mr*mr-1.0d0))*dmrr_r 723 : dMmu_r=(temph-4.0d0*beta*mr*(mr*mr-1.0d0))*dmru_r 724 : dMmv_r=(temph-4.0d0*beta*mr*(mr*mr-1.0d0))*dmrv_r 725 : dMmp_r=(temph-4.0d0*beta*mr*(mr*mr-1.0d0))*dmrp_r 726 : do 54 i=1,(nsp-1) 727 : dMmy_r.vv(i)=(temph-4.0d0*beta*mr* 728 : & (mr*mr-1.0d0))*dmry_r.vv(i) 729 : 54 continue 730 : else 731 : dMmr_l=dmrr_l 732 : dMmu_l=dmru_l 733 : dMmv_l=dmrv_l 734 : dMmp_l=dmrp_l 735 : do 55 i=1,(nsp-1) 736 : dMmy_l.vv(i)=dmry_l.vv(i) 737 : 55 continue 738 : c-------------------- 739 : dMmr_r=dmrr_r 740 : dMmu_r=dmru_r 741 : dMmv_r=dmrv_r 742 : dMmp_r=dmrp_r 743 : do 56 i=1,(nsp-1) 744 : dMmy_r.vv(i)=dmry_r.vv(i) 745 : 56 continue 746 : endif 747 : endif 748 : c----------------------------------------------------------- 749 : c addition of low Mach number 750 : c----------------------------------------------------------- 751 : if(mr .le. 0.0d0) then 752 : m1mr_l=dmrr_l 753 : m1mu_l=dmru_l 754 : m1mv_l=dmrv_l 755 : m1mp_l=dmrp_l 756 : do 550 i=1,(nsp-1) 757 : m1my_l.vv(i)=dmry_l.vv(i) 758 : 550 continue 759 : c--------------------- 760 : m1mr_r=dmrr_r 761 : m1mu_r=dmru_r 762 : m1mv_r=dmrv_r 763 : m1mp_r=dmrp_r 764 : do 560 i=1,(nsp-1) 765 : m1my_r.vv(i)=dmry_r.vv(i) 766 : 560 continue 767 : else 768 : m1mr_l=0.0d0 769 : m1mu_l=0.0d0 770 : m1mv_l=0.0d0 771 : m1mp_l=0.0d0 772 : do 551 i=1,(nsp-1) 773 : m1my_l.vv(i)=0.0d0 774 : 551 continue 775 : c-------------------- 776 : m1mr_r=0.0d0 777 : m1mu_r=0.0d0 778 : m1mv_r=0.0d0 779 : m1mp_r=0.0d0 780 : do 561 i=1,(nsp-1) 781 : m1my_r.vv(i)=0.0d0 782 : 561 continue 783 : endif 784 : c----------------------------------------------------------------- 785 : c computing the derivatives of m_{1/2} (notation as in the paper) 786 : c----------------------------------------------------------------- 787 : dmir_l=dMpr_l+dMmr_l 788 : dmir_r=dMpr_r+dMmr_r 789 : c------------- 790 : dmiu_l=dMpu_l+dMmu_l 791 : dmiu_r=dMpu_r+dMmu_r 792 : c------------- 793 : dmiv_l=dMpv_l+dMmv_l 794 : dmiv_r=dMpv_r+dMmv_r 795 : c------------- 796 : dmip_l=dMpp_l+dMmp_l 797 : dmip_r=dMpp_r+dMmp_r 798 : c------------- 799 : do 57 i=1,(nsp-1) 800 : dmiy_l.vv(i)=dMpy_l.vv(i)+dMmy_l.vv(i) 801 : dmiy_r.vv(i)=dMpy_r.vv(i)+dMmy_r.vv(i) 802 : 57 continue 803 : c---------------------------------------------------------------- 804 : c computing the main convective variables and their derivatives 805 : c mpl_m is m^{+}_{1/2} (paper's notation) and 806 : c mmin_m is m^{-}_{1/2} (paper's notation), see (A2) on p.370. 807 : c---------------------------------------------------------------- 808 : termp=(Mmin1-Mmin+Mplus-Mplus1)*(1.0d0/(mhalfr*mhalfr)-1.0d0) 809 : termp=termp*(pold_l-pold_r)/(pold_l/rold_l+pold_r/rold_r) 810 : c------------------------------------------------------------- 811 : c derivatives of the termp 812 : c------------------------------------------------------------- 813 : top=(Mmin1-Mmin+Mplus-Mplus1) 814 : bots=1.0d0/(pold_l/rold_l+pold_r/rold_r) 815 : bot=(pold_l-pold_r)*bots 816 : temph=1.0d0/(mhalfr*mhalfr)-1.0d0 817 : c--------------------------- 818 : tmpr_l=(m1mr_l-dMmr_l+dMpr_l-m1pr_l)*bot*temph 819 : tmpr_l=tmpr_l-2.0d0*bot*top*dmfr_l/(mhalfr*mhalfr*mhalfr) 820 : tmpr_l=tmpr_l+temph*top*bot*bots*(pold_l/rold_l/rold_l) 821 : c--------------------------- 822 : tmpu_l=(m1mu_l-dMmu_l+dMpu_l-m1pu_l)*bot*temph 823 : tmpu_l=tmpu_l-2.0d0*bot*top*dmfu_l/(mhalfr*mhalfr*mhalfr) 824 : c--------------------------- 825 : tmpv_l=(m1mv_l-dMmv_l+dMpv_l-m1pv_l)*bot*temph 826 : tmpv_l=tmpv_l-2.0d0*bot*top*dmfv_l/(mhalfr*mhalfr*mhalfr) 827 : c--------------------------- 828 : tmpp_l=(m1mp_l-dMmp_l+dMpp_l-m1pp_l)*bot*temph 829 : tmpp_l=tmpp_l-2.0d0*bot*top*dmfp_l/(mhalfr*mhalfr*mhalfr) 830 : tmpp_l=tmpp_l+temph*top*bots*(1.0d0-bot/rold_l) 831 : c--------------------------- 832 : do 570 i=1,(nsp-1) 833 : tmpy_l.vv(i)=(m1my_l.vv(i)-dMmy_l.vv(i)+ 834 : & dMpy_l.vv(i)-m1py_l.vv(i))*bot*temph 835 : tmpy_l.vv(i)=tmpy_l.vv(i)-2.0d0*bot*top* 836 : & dmfg_l*dgdyl.vv(i)/(mhalfr*mhalfr*mhalfr) 837 : 570 continue 838 : c------------rrrrrrrr------- 839 : c------------rrrrrrrr------- 840 : tmpr_r=(m1mr_r-dMmr_r+dMpr_r-m1pr_r)*bot*temph 841 : tmpr_r=tmpr_r-2.0d0*bot*top*dmfr_r/(mhalfr*mhalfr*mhalfr) 842 : tmpr_r=tmpr_r+temph*top*bot*bots*(pold_r/rold_r/rold_r) 843 : c--------------------------- 844 : tmpu_r=(m1mu_r-dMmu_r+dMpu_r-m1pu_r)*bot*temph 845 : tmpu_r=tmpu_r-2.0d0*bot*top*dmfu_r/(mhalfr*mhalfr*mhalfr) 846 : c--------------------------- 847 : tmpv_r=(m1mv_r-dMmv_r+dMpv_r-m1pv_r)*bot*temph 848 : tmpv_r=tmpv_r-2.0d0*bot*top*dmfv_r/(mhalfr*mhalfr*mhalfr) 849 : c--------------------------- 850 : tmpp_r=(m1mp_r-dMmp_r+dMpp_r-m1pp_r)*bot*temph 851 : tmpp_r=tmpp_r-2.0d0*bot*top*dmfp_r/(mhalfr*mhalfr*mhalfr) 852 : tmpp_r=tmpp_r-temph*top*bots*(1.0d0+bot/rold_r) 853 : c--------------------------- 854 : do 571 i=1,(nsp-1) 855 : tmpy_r.vv(i)=(m1my_r.vv(i)-dMmy_r.vv(i)+ 856 : & dMpy_r.vv(i)-m1py_r.vv(i))*bot*temph 857 : tmpy_r.vv(i)=tmpy_r.vv(i)-2.0d0*bot*top* 858 : & dmfg_r*dgdyr.vv(i)/(mhalfr*mhalfr*mhalfr) 859 : 571 continue 860 : c------------------------------------------------------------- 861 : if(mmid .ge. 0.0d0) then 862 : mpl_m = mmid 863 : d2mr_l=dmir_l 864 : d2mu_l=dmiu_l 865 : d2mv_l=dmiv_l 866 : d2mp_l=dmip_l 867 : do 58 i=1,(nsp-1) 868 : d2my_l.vv(i)=dmiy_l.vv(i) 869 : 58 continue 870 : c------------ 871 : d2mr_r=dmir_r 872 : d2mu_r=dmiu_r 873 : d2mv_r=dmiv_r 874 : d2mp_r=dmip_r 875 : do 59 i=1,(nsp-1) 876 : d2my_r.vv(i)=dmiy_r.vv(i) 877 : 59 continue 878 : c------------ 879 : else 880 : mpl_m = 0.0d0 881 : d2mr_l=0.0d0 882 : d2mu_l=0.0d0 883 : d2mv_l=0.0d0 884 : d2mp_l=0.0d0 885 : do 60 i=1,(nsp-1) 886 : d2my_l.vv(i)=0.0d0 887 : 60 continue 888 : c------------ 889 : d2mr_r=0.0d0 890 : d2mu_r=0.0d0 891 : d2mv_r=0.0d0 892 : d2mp_r=0.0d0 893 : do 61 i=1,(nsp-1) 894 : d2my_r.vv(i)=0.0d0 895 : 61 continue 896 : endif 897 : c--------------------------------------------- 898 : cc derivatives for the term termm 899 : cc------------------------------------------------------------------ 900 : if(mmid .le. 0.0d0) then 901 : mmin_m = mmid 902 : d3mr_l=dmir_l 903 : d3mu_l=dmiu_l 904 : d3mv_l=dmiv_l 905 : d3mp_l=dmip_l 906 : do 62 i=1,(nsp-1) 907 : d3my_l.vv(i)=dmiy_l.vv(i) 908 : 62 continue 909 : c------------ 910 : d3mr_r=dmir_r 911 : d3mu_r=dmiu_r 912 : d3mv_r=dmiv_r 913 : d3mp_r=dmip_r 914 : do 63 i=1,(nsp-1) 915 : d3my_r.vv(i)=dmiy_r.vv(i) 916 : 63 continue 917 : c------------ 918 : else 919 : mmin_m = 0.0d0 920 : d3mr_l=0.0d0 921 : d3mu_l=0.0d0 922 : d3mv_l=0.0d0 923 : d3mp_l=0.0d0 924 : do 64 i=1,(nsp-1) 925 : d3my_l.vv(i)=0.0d0 926 : 64 continue 927 : c------------ 928 : d3mr_r=0.0d0 929 : d3mu_r=0.0d0 930 : d3mv_r=0.0d0 931 : d3mp_r=0.0d0 932 : do 65 i=1,(nsp-1) 933 : d3my_r.vv(i)=0.0d0 934 : 65 continue 935 : endif 936 : c--------------------------------------------------------------- 936 : c Computing the calligraphic P+ and P- with their derivatives, 937 : c see (21a) & (21b) on p.368. 938 : c--------------------------------------------------------------- 939 : if(ml .ge. 1.0d0) then 940 : Pplus = 1.0d0 941 : else 942 : if((ml .gt. -1.0d0) .and. (ml .lt. 1.0d0)) then 943 : Pplus=(ml+1.0d0)*(ml+1.0d0)*(2.0d0-ml)/4.0d0 944 : Pplus=Pplus+alpha*ml*(ml*ml-1.0d0)*(ml*ml-1.0d0) 945 : else 946 : Pplus = 0.0d0 947 : endif 948 : endif 949 : c--------------------------------------------------------------- 950 : if(mr .ge. 1.0d0) then 951 : Pmin = 0.0d0 952 : else 953 : if((mr .gt. -1.0d0) .and. (mr .lt. 1.0d0)) then 954 : Pmin=(mr-1.0d0)*(mr-1.0d0)*(2.0d0+mr)/4.0d0 955 : Pmin=Pmin-alpha*mr*(mr*mr-1.0d0)*(mr*mr-1.0d0) 956 : else 957 : Pmin = 1.0d0 958 : endif 959 : endif 960 : c--------------------------------------------------------------- 961 : brac_l=(ml+1.0d0)*(2.0d0-ml)/2.0d0-(ml+1.0d0)*(ml+1.0d0)/4.0d0 962 : brac_l=brac_l+alpha*(ml*ml-1.0d0)*(ml*ml-1.0d0) 963 : brac_l=brac_l+4.0d0*alpha*ml*ml*(ml*ml-1.0d0) 964 : c-------------- 965 : brac_r=(mr-1.0d0)*(2.0d0+mr)/2.0d0+(mr-1.0d0)*(mr-1.0d0)/4.0d0 966 : brac_r=brac_r-alpha*(mr*mr-1.0d0)*(mr*mr-1.0d0) 967 : brac_r=brac_r-4.0d0*alpha*mr*mr*(mr*mr-1.0d0) 968 : c--------------------------------------------------------------- 969 : if((ml .gt. -1.0d0) .and. (ml .lt. 1.0d0)) then 970 : dPpr_l=brac_l*dmlr_l 971 : dPpr_r=brac_l*dmlr_r 972 : c------------ 973 : dPpu_l=brac_l*dmlu_l 974 : dPpu_r=brac_l*dmlu_r 975 : c------------ 976 : dPpv_l=brac_l*dmlv_l 977 : dPpv_r=brac_l*dmlv_r 978 : c------------ 979 : dPpp_l=brac_l*dmlp_l 980 : dPpp_r=brac_l*dmlp_r 981 : c------------ 982 : do 66 i=1,(nsp-1) 983 : dPpy_l.vv(i)=brac_l*dmly_l.vv(i) 984 : dPpy_r.vv(i)=brac_l*dmly_r.vv(i) 985 : 66 continue 986 : c------------ 987 : else 988 : dPpr_l=0.0d0 989 : dPpr_r=0.0d0 990 : c----------- 991 : dPpu_l=0.0d0 992 : dPpu_r=0.0d0 993 : c----------- 994 : dPpv_l=0.0d0 995 : dPpv_r=0.0d0 996 : c----------- 997 : dPpp_l=0.0d0 998 : dPpp_r=0.0d0 999 : c----------- 1000 : do 67 i=1,(nsp-1) 1001 : dPpy_l.vv(i)=0.0d0 1002 : dPpy_r.vv(i)=0.0d0 1003 : 67 continue 1004 : c----------- 1005 : endif 1006 : c--------------------------------------------------------------- 1007 : if((mr .gt. -1.0d0) .and. (mr .lt. 1.0d0)) then 1008 : dPmr_l=brac_r*dmrr_l 1009 : dPmr_r=brac_r*dmrr_r 1010 : c------------ 1011 : dPmu_l=brac_r*dmru_l 1012 : dPmu_r=brac_r*dmru_r 1013 : c------------ 1014 : dPmv_l=brac_r*dmrv_l 1015 : dPmv_r=brac_r*dmrv_r 1016 : c------------ 1017 : dPmp_l=brac_r*dmrp_l 1018 : dPmp_r=brac_r*dmrp_r 1019 : c------------ 1020 : do 68 i=1,(nsp-1) 1021 : dPmy_l.vv(i)=brac_r*dmry_l.vv(i) 1022 : dPmy_r.vv(i)=brac_r*dmry_r.vv(i) 1023 : 68 continue 1024 : c------------ 1025 : else 1026 : dPmr_l=0.0d0 1027 : dPmr_r=0.0d0 1028 : c----------- 1029 : dPmu_l=0.0d0 1030 : dPmu_r=0.0d0 1031 : c----------- 1032 : dPmv_l=0.0d0 1033 : dPmv_r=0.0d0 1034 : c----------- 1035 : dPmp_l=0.0d0 1036 : dPmp_r=0.0d0 1037 : c----------- 1038 : do 69 i=1,(nsp-1) 1039 : dPmy_l.vv(i)=0.0d0 1040 : dPmy_r.vv(i)=0.0d0 1041 : 69 continue 1042 : c----------- 1043 : endif 1044 : c------------------------------------------------------------------- 1045 : c computing pmid -- p_{1/2} and its derivatives, see (20b), p.367. 1046 : c------------------------------------------------------------------- 1047 : dpir_l=dPpr_l*pold_l+dPmr_l*pold_r 1048 : dpiu_l=dPpu_l*pold_l+dPmu_l*pold_r 1049 : dpiv_l=dPpv_l*pold_l+dPmv_l*pold_r 1050 : dpip_l=dPpp_l*pold_l+Pplus+dPmp_l*pold_r 1051 : do 70 i=1,(nsp-1) 1052 : dpiy_l.vv(i)=dPpy_l.vv(i)*pold_l+dPmy_l.vv(i)*pold_r 1053 : 70 continue 1054 : c---------------------------- 1055 : dpir_r=dPpr_r*pold_l+dPmr_r*pold_r 1056 : dpiu_r=dPpu_r*pold_l+dPmu_r*pold_r 1057 : dpiv_r=dPpv_r*pold_l+dPmv_r*pold_r 1058 : dpip_r=dPpp_r*pold_l+Pmin+dPmp_r*pold_r 1059 : do 71 i=1,(nsp-1) 1060 : dpiy_r.vv(i)=dPpy_r.vv(i)*pold_l+dPmy_r.vv(i)*pold_r 1061 : 71 continue 1062 : c--------------------------------------------------------------------- 1063 : c Computation of the mass flux (rho * u)_1/2 1064 : c--------------------------------------------------------------- 1065 : rum=am*(mpl_m*rold_l+mmin_m*rold_r)+canc*am*termp 1066 : c------------------------------------------------------- 1067 : rumr_l=damr_l*(mpl_m*rold_l+mmin_m*rold_r)+ 1068 : & am*(d2mr_l*rold_l+mpl_m+d3mr_l*rold_r) 1069 : rumr_l=rumr_l+canc*coef*(damr_l*termp+am*tmpr_l) 1070 : rumu_l=damu_l*(mpl_m*rold_l+mmin_m*rold_r)+ 1071 : & am*(d2mu_l*rold_l+d3mu_l*rold_r) 1072 : rumu_l=rumu_l+canc*coef*(damu_l*termp+am*tmpu_l) 1073 : rumv_l=damv_l*(mpl_m*rold_l+mmin_m*rold_r)+ 1074 : & am*(d2mv_l*rold_l+d3mv_l*rold_r) 1075 : rumv_l=rumv_l+canc*coef*(damv_l*termp+am*tmpv_l) 1076 : rump_l=damp_l*(mpl_m*rold_l+mmin_m*rold_r)+ 1077 : & am*(d2mp_l*rold_l+d3mp_l*rold_r) 1078 : rump_l=rump_l+canc*coef*(damp_l*termp+am*tmpp_l) 1079 : c------------------ 1080 : do 710 i=1,(nsp-1) 1081 : rumy_l.vv(i)=damg_l*dgdyl.vv(i)*(mpl_m*rold_l+ 1082 : & mmin_m*rold_r)+am*(d2my_l.vv(i)*rold_l+d3my_l.vv(i)*rold_r) 1083 : rumy_l.vv(i)=rumy_l.vv(i)+canc*coef* 1084 : & (damg_l*dgdyl.vv(i)*termp+am*tmpy_l.vv(i)) 1085 : 710 continue 1086 : c------------------------------------------------- 1087 : rumr_r=damr_r*(mpl_m*rold_l+mmin_m*rold_r)+ 1088 : & am*(d2mr_r*rold_l+mmin_m+d3mr_r*rold_r) 1089 : rumr_r=rumr_r+canc*coef*(damr_r*termp+am*tmpr_r) 1090 : rumu_r=damu_r*(mpl_m*rold_l+mmin_m*rold_r)+ 1091 : & am*(d2mu_r*rold_l+d3mu_r*rold_r) 1092 : rumu_r=rumu_r+canc*coef*(damu_r*termp+am*tmpu_r) 1093 : rumv_r=damv_r*(mpl_m*rold_l+mmin_m*rold_r)+ 1094 : & am*(d2mv_r*rold_l+d3mv_r*rold_r) 1095 : rumv_r=rumv_r+canc*coef*(damv_r*termp+am*tmpv_r) 1096 : rump_r=damp_r*(mpl_m*rold_l+mmin_m*rold_r)+ 1097 : & am*(d2mp_r*rold_l+d3mp_r*rold_r) 1098 : rump_r=rump_r+canc*coef*(damp_r*termp+am*tmpp_r) 1099 : do 711 i=1,(nsp-1) 1100 : rumy_r.vv(i)=damg_r*dgdyr.vv(i)*(mpl_m*rold_l+mmin_m*rold_r)+ 1101 : & am*(d2my_r.vv(i)*rold_l+d3my_r.vv(i)*rold_r) 1102 : rumy_r.vv(i)=rumy_r.vv(i)+canc*coef* 1103 : & (damg_r*dgdyr.vv(i)*termp+am*tmpy_r.vv(i)) 1104 : 711 continue 1105 : c--------------------------------------------------------------------- 1106 : c computing JACOBIAN as a derivative of the numerical flux function 1107 : c with respect to the primitive variables. 1108 : c Notation: jl(i,j) --- is the derivative of the i-component of the 1109 : c flux function with respect to the j-component of the 1110 : c vector of primitive variables of the left state. 1111 : c jr(i,j) --- is the derivative of the i-component of the 1112 : c flux function with respect to the j-component of the 1113 : c vector of primitive variables of the right state. 1114 : c--------------------------------------------------------------------- 1115 : jl.jac(1,1)=0.0d0 1116 : jl.jac(1,2)=0.0d0 1117 : jl.jac(1,3)=0.0d0 1118 : jl.jac(1,4)=0.0d0 1119 : do 72 i=1,(nsp-1) 1120 : jl.jac(1,4+i)=0.0d0 1121 : 72 continue 1122 : c------------------------------------ 1123 : jr.jac(1,1)=0.0d0 1124 : jr.jac(1,2)=0.0d0 1125 : jr.jac(1,3)=0.0d0 1126 : jr.jac(1,4)=0.0d0 1127 : do 73 i=1,(nsp-1) 1128 : jr.jac(1,4+i)=0.0d0 1129 : 73 continue 1130 : c------------------------------------ 1131 : c--------------------------------------------------------- 1132 : if(rum .ge. 0.0d0) then 1133 : br3=rumr_l*un_l 1134 : br4=rumr_l*ut_l 1135 : else 1136 : br3=rumr_l*un_r 1137 : br4=rumr_l*ut_r 1138 : endif 1139 : jl.jac(2,1)=n_x*(br3+dpir_l)+br4*t_x 1140 : jl.jac(3,1)=n_y*(br3+dpir_l)+br4*t_y 1141 : c------------------- 1142 : if(rum .ge. 0.0d0) then 1143 : br3=rumu_l*un_l+rum*n_x 1144 : br4=rumu_l*ut_l+rum*t_x 1145 : else 1146 : br3=rumu_l*un_r 1147 : br4=rumu_l*ut_r 1148 : endif 1149 : jl.jac(2,2)=n_x*(br3+dpiu_l)+br4*t_x 1150 : jl.jac(3,2)=n_y*(br3+dpiu_l)+br4*t_y 1151 : c------------------- 1152 : if(rum .ge. 0.0d0) then 1153 : br3=rumv_l*un_l+rum*n_y 1154 : br4=rumv_l*ut_l+rum*t_y 1155 : else 1156 : br3=rumv_l*un_r 1157 : br4=rumv_l*ut_r 1158 : endif 1159 : jl.jac(2,3)=n_x*(br3+dpiv_l)+br4*t_x 1160 : jl.jac(3,3)=n_y*(br3+dpiv_l)+br4*t_y 1161 : c------------------- 1162 : if(rum .ge. 0.0d0) then 1163 : br3=rump_l*un_l 1164 : br4=rump_l*ut_l 1165 : else 1166 : br3=rump_l*un_r 1167 : br4=rump_l*ut_r 1168 : endif 1169 : jl.jac(2,4)=n_x*(br3+dpip_l)+br4*t_x 1170 : jl.jac(3,4)=n_y*(br3+dpip_l)+br4*t_y 1171 : c------------------------------------------------------------- 1172 : do 74 i=1,(nsp-1) 1173 : if(rum .ge. 0.0d0) then 1174 : br3=rumy_l.vv(i)*un_l 1175 : br4=rumy_l.vv(i)*ut_l 1176 : else 1177 : br3=rumy_l.vv(i)*un_r 1178 : br4=rumy_l.vv(i)*ut_r 1179 : endif 1180 : jl.jac(2,4+i)=n_x*(br3+dpiy_l.vv(i))+br4*t_x 1181 : jl.jac(3,4+i)=n_y*(br3+dpiy_l.vv(i))+br4*t_y 1182 : 74 continue 1183 : c------------------------------------------------------------- 1184 : if(rum .ge. 0.0d0) then 1185 : br3=rumr_r*un_l 1186 : br4=rumr_r*ut_l 1187 : else 1188 : br3=rumr_r*un_r 1189 : br4=rumr_r*ut_r 1190 : endif 1191 : jr.jac(2,1)=n_x*(br3+dpir_r)+br4*t_x 1192 : jr.jac(3,1)=n_y*(br3+dpir_r)+br4*t_y 1193 : c------------------- 1194 : if(rum .ge. 0.0d0) then 1195 : br3=rumu_r*un_l 1196 : br4=rumu_r*ut_l 1197 : else 1198 : br3=rumu_r*un_r+rum*n_x 1199 : br4=rumu_r*ut_r+rum*t_x 1200 : endif 1201 : jr.jac(2,2)=n_x*(br3+dpiu_r)+br4*t_x 1202 : jr.jac(3,2)=n_y*(br3+dpiu_r)+br4*t_y 1203 : c------------------- 1204 : if(rum .ge. 0.0d0) then 1205 : br3=rumv_r*un_l 1206 : br4=rumv_r*ut_l 1207 : else 1208 : br3=rumv_r*un_r+rum*n_y 1209 : br4=rumv_r*ut_r+rum*t_y 1210 : endif 1211 : jr.jac(2,3)=n_x*(br3+dpiv_r)+br4*t_x 1212 : jr.jac(3,3)=n_y*(br3+dpiv_r)+br4*t_y 1213 : c------------------- 1214 : if(rum .ge. 0.0d0) then 1215 : br3=rump_r*un_l 1216 : br4=rump_r*ut_l 1217 : else 1218 : br3=rump_r*un_r 1219 : br4=rump_r*ut_r 1220 : endif 1221 : jr.jac(2,4)=n_x*(br3+dpip_r)+br4*t_x 1222 : jr.jac(3,4)=n_y*(br3+dpip_r)+br4*t_y 1223 : c-------------------- 1224 : do 75 i=1,(nsp-1) 1225 : if(rum .ge. 0.0d0) then 1226 : br3=rumy_r.vv(i)*un_l 1227 : br4=rumy_r.vv(i)*ut_l 1228 : else 1229 : br3=rumy_r.vv(i)*un_r 1230 : br4=rumy_r.vv(i)*ut_r 1231 : endif 1232 : jr.jac(2,4+i)=n_x*(br3+dpiy_r.vv(i))+br4*t_x 1233 : jr.jac(3,4+i)=n_y*(br3+dpiy_r.vv(i))+br4*t_y 1234 : 75 continue 1235 : c------------------------------------------------------------- 1236 : c ------ f44444444444444444444444444444 --------- 1237 : c------------------------------------------------------------- 1238 : jl.jac(4,1)=0.0d0 1239 : jl.jac(4,2)=0.0d0 1240 : jl.jac(4,3)=0.0d0 1241 : jl.jac(4,4)=0.0d0 1242 : c---------------------------------------------------------- 1243 : jr.jac(4,1)=0.0d0 1244 : jr.jac(4,2)=0.0d0 1245 : jr.jac(4,3)=0.0d0 1246 : jr.jac(4,4)=0.0d0 1247 : c--------------------- 1248 : do 76 i=1,(nsp-1) 1249 : jl.jac(4,4+i)=0.0d0 1250 : jr.jac(4,4+i)=0.0d0 1251 : 76 continue 1252 : c------------------------------------------------------------- 1253 : c------------------------------------------------------------- 1254 : do 78 i=1,(nsp-1) 1255 : jl.jac(4+i,1)=0.0d0 1256 : jr.jac(4+i,1)=0.0d0 1257 : c-------------------------- 1258 : jl.jac(4+i,2)=0.0d0 1259 : jr.jac(4+i,2)=0.0d0 1260 : c-------------------------- 1261 : jl.jac(4+i,3)=0.0d0 1262 : jr.jac(4+i,3)=0.0d0 1263 : c-------------------------- 1264 : jl.jac(4+i,4)=0.0d0 1265 : jr.jac(4+i,4)=0.0d0 1266 : c-------------------------- 1267 : do 780 j=1,(nsp-1) 1268 : jl.jac(4+i,4+j)=0.0d0 1269 : jr.jac(4+i,4+j)=0.0d0 1270 : 780 continue 1271 : 78 continue 1272 : c------------------------------------------------------------- 1273 : c matrix wl(i,j) represents the derivative of the i-component 1274 : c of the vector of primitive variables of the left state with 1275 : c respect to the j-component of the vector of the conservative 1276 : c variables of the left state. 1277 : c 1278 : c Here: (rho, ux, uy, p, Y_1,...,Y_(nsp-1)) - 1279 : c vector of primitive variables; 1280 : c (rho, rho ux, rho uy, rho e, rho Y_1,..., rho Y_(nsp-1)) - 1281 : c vector of conservative variables. 1282 : c------------------------------------------------------------- 1283 : wl.jac(1,1)=1.0d0 1284 : wl.jac(1,2)=0.0d0 1285 : wl.jac(1,3)=0.0d0 1286 : wl.jac(1,4)=0.0d0 1287 : do 83 i=1,(nsp-1) 1288 : wl.jac(1,4+i)=0.0d0 1289 : 83 continue 1290 : c------------------------------ 1291 : wl.jac(2,1)=-uold_l/rold_l 1292 : wl.jac(2,2)=1.0d0/rold_l 1293 : wl.jac(2,3)=0.0d0 1294 : wl.jac(2,4)=0.0d0 1295 : do 84 i=1,(nsp-1) 1296 : wl.jac(2,4+i)=0.0d0 1297 : 84 continue 1298 : c------------------------------ 1299 : wl.jac(3,1)=-vold_l/rold_l 1300 : wl.jac(3,2)=0.0d0 1301 : wl.jac(3,3)=1.0d0/rold_l 1302 : wl.jac(3,4)=0.0d0 1303 : do 85 i=1,(nsp-1) 1304 : wl.jac(3,4+i)=0.0d0 1305 : 85 continue 1306 : c------------------------------ 1307 : br1=0.0d0 1308 : do 86 i=1,(nsp-1) 1309 : br1=br1+dgdyl.vv(i)*yl.yet(i) 1310 : 86 continue 1311 : br1=br1*pold_l/(rold_l*gm1l) 1312 : wl.jac(4,1)=gm1l*(uold_l*uold_l+vold_l*vold_l)/2.0d0-br1 1313 : wl.jac(4,2)=-uold_l*gm1l 1314 : wl.jac(4,3)=-vold_l*gm1l 1315 : wl.jac(4,4)=gm1l 1316 : do 87 i=1,(nsp-1) 1317 : wl.jac(4,4+i)=dgdyl.vv(i)*pold_l/(rold_l*gm1l) 1318 : 87 continue 1319 : c------------------------------ 1320 : do 88 i=1,(nsp-1) 1321 : do 89 j=1,4 1322 : wl.jac(4+i,j)=0.0d0 1323 : if(j.eq.1) wl.jac(4+i,j)=-yl.yet(i)/rold_l 1324 : 89 continue 1325 : c------------ 1326 : do 890 j=5,(4+nsp-1) 1327 : wl.jac(4+i,j)=0.0d0 1328 : if(4+i.eq.j) then 1329 : wl.jac(4+i,j)=1.0d0/rold_l 1330 : endif 1331 : 890 continue 1332 : 88 continue 1333 : c------------------------------ 1334 : c------------------------------ 1335 : wr.jac(1,1)=1.0d0 1336 : wr.jac(1,2)=0.0d0 1337 : wr.jac(1,3)=0.0d0 1338 : wr.jac(1,4)=0.0d0 1339 : do 90 i=1,(nsp-1) 1340 : wr.jac(1,4+i)=0.0d0 1341 : 90 continue 1342 : c------------------------------ 1343 : wr.jac(2,1)=-uold_r/rold_r 1344 : wr.jac(2,2)=1.0d0/rold_r 1345 : wr.jac(2,3)=0.0d0 1346 : wr.jac(2,4)=0.0d0 1347 : do 91 i=1,(nsp-1) 1348 : wr.jac(2,4+i)=0.0d0 1349 : 91 continue 1350 : c------------------------------ 1351 : wr.jac(3,1)=-vold_r/rold_r 1352 : wr.jac(3,2)=0.0d0 1353 : wr.jac(3,3)=1.0d0/rold_r 1354 : wr.jac(3,4)=0.0d0 1355 : do 92 i=1,(nsp-1) 1356 : wr.jac(3,4+i)=0.0d0 1357 : 92 continue 1358 : c------------------------------ 1359 : br1=0.0d0 1360 : do 860 i=1,(nsp-1) 1361 : br1=br1+dgdyr.vv(i)*yr.yet(i) 1362 : 860 continue 1363 : br1=br1*pold_r/(rold_r*gm1r) 1364 : wr.jac(4,1)=gm1r*(uold_r*uold_r+vold_r*vold_r)/2.0d0-br1 1365 : wr.jac(4,2)=-uold_r*gm1r 1366 : wr.jac(4,3)=-vold_r*gm1r 1367 : wr.jac(4,4)=gm1r 1368 : do 93 i=1,(nsp-1) 1369 : wr.jac(4,4+i)=dgdyr.vv(i)*pold_r/(rold_r*gm1r) 1370 : 93 continue 1371 : c------------------------------ 1372 : do 94 i=1,(nsp-1) 1373 : do 95 j=1,4 1374 : wr.jac(4+i,j)=0.0d0 1375 : if(j.eq.1) wr.jac(4+i,j)=-yr.yet(i)/rold_r 1376 : 95 continue 1377 : c---------------------- 1378 : do 950 j=5,(3+nsp) 1379 : wr.jac(4+i,j)=0.0d0 1380 : if(4+i.eq.j) wr.jac(4+i,j)=1.0d0/rold_r 1381 : 950 continue 1382 : 94 continue 1383 : c---------------------------------------------- 1384 : c---------------------------------------------- 1385 : do 1 i=1,(3+nsp) 1386 : do 2 j=1,(3+nsp) 1387 : jtl.jac(i,j)=0.0d0 1388 : jtr.jac(i,j)=0.0d0 1389 : do 3 k=1,(3+nsp) 1390 : jtl.jac(i,j)=jtl.jac(i,j)+jl.jac(i,k)*wl.jac(k,j) 1391 : jtr.jac(i,j)=jtr.jac(i,j)+jr.jac(i,k)*wr.jac(k,j) 1392 : 3 continue 1393 : 2 continue 1394 : 1 continue 1395 : c---------------------------------------------------------------------- 1396 : tcoef=nvect(1)*tvect(2)+tvect(1)*nvect(2) 1397 : bcoef=nvect(1)*tvect(2)-tvect(1)*nvect(2) 1398 : c---------------------------------------------------------------------- 1399 : do 11 i=1,(3+nsp) 1400 : jtl.jac(i,1)=jtl.jac(i,1)+jtr.jac(i,1) 1401 : jtl.jac(i,2)=jtl.jac(i,2)+jtr.jac(i,2)*(-tcoef/bcoef)+ 1402 : & jtr.jac(i,3)*2.0d0*nvect(1)*tvect(1)/bcoef 1403 : jtl.jac(i,3)=jtl.jac(i,3)+jtr.jac(i,2)* 1404 : & (-2.0d0*nvect(2)*tvect(2)/bcoef)+ 1405 : & jtr.jac(i,3)*tcoef/bcoef 1406 : jtl.jac(i,4)=jtl.jac(i,4)+jtr.jac(i,4) 1407 : do 777 j=5,(3+nsp) 1408 : jtl.jac(i,j)=jtl.jac(i,j)+jtr.jac(i,j) 1409 : 777 continue 1410 : 11 continue 1411 : c---------------------------------- 1412 : c do 11 i=1,(3+nsp) 1413 : c jtl.jac(i,1)=jtl.jac(i,1)+jtr.jac(i,1) 1414 : c jtl.jac(i,2)=jtl.jac(i,2)+jtr.jac(i,2)*(-tcoef/bcoef)+ 1415 : c & jtr.jac(i,3)*2.0d0*nvect(1)*tvect(1)/bcoef 1416 : c jtl.jac(i,3)=jtl.jac(i,3)+jtr.jac(i,2)* 1417 : c & (-2.0d0*nvect(2)*tvect(2)/bcoef)+ 1418 : c & jtr.jac(i,3)*tcoef/bcoef 1419 : c jtl.jac(i,4)=jl.jac(i,4)+jr.jac(i,4) 1420 : c do 777 j=5,(3+nsp) 1421 : c jl.jac(i,j)=jl.jac(i,j)+jr.jac(i,j) 1422 : c 777 continue 1423 : c 11 continue 1424 : c------------------------------------------------------------------- 1425 : SEGSUP DMLY_L, DMLY_R, 1426 : & dmry_l, dmry_r, 1427 : & dMpy_l, dMpy_r, 1428 : & dMmy_l, dMmy_r, 1429 : & dmiy_l, dmiy_r, 1430 : & d3my_l, d3my_r, 1431 : & d2my_l, d2my_r, 1432 : & dPpy_l, dPpy_r, 1433 : & dPmy_l, dPmy_r, 1434 : & dpiy_l, dpiy_r, 1435 : & dgdyl, dgdyr 1436 : C-------------------------------------------- 1437 : SEGSUP f 1438 : C-------------------------------------------- 1439 : c jll = jl 1440 : SEGDES JL 1441 : SEGDES JR 1442 : SEGSUP WL 1443 : SEGSUP WR 1444 : C-------------------------------------------- 1445 : jll=jtl 1446 : SEGDES JTL 1447 : SEGDES JTR 1448 : SEGDES YL 1449 : SEGDES YR 1450 : SEGDES CP 1451 : SEGDES CV 1452 : SEGDES MLRECP, MLRECV 1453 : c------------------- 1454 : return 1455 : end 1456 : 1457 : 1458 : 1459 : 1460 : 1461 : 1462 : 1463 : 1464 : 1465 : 1466 : 1467 : 1468 : 1469 : 1470 : 1471 : 1472 : 1473 : 1474 : 1475 : 1476 : 1477 : 1478 : 1479 : 1480 : 1481 : 1482 : 1483 : 1484 : 1485 : 1486 : 1487 : 1488 : 1489 :
© Cast3M 2003 - All rights reserved.
Disclaimer