67 ddens_, momx_, momy_, momz_, drhot_, dpres, &
68 dens_hyd, pres_hyd, therm_hyd, rtot, cvtot, cptot, &
69 gsqrt, g13, g23, nx, ny, nz, &
70 vmapm, vmapp, lmesh, elem, lmesh2d, elem2d )
78 real(rp),
intent(out) :: del_flux(elem%nfptot,
prgvar_num,lmesh%ne)
79 real(rp),
intent(in) :: ddens_(elem%np*lmesh%nea)
80 real(rp),
intent(in) :: momx_(elem%np*lmesh%nea)
81 real(rp),
intent(in) :: momy_(elem%np*lmesh%nea)
82 real(rp),
intent(in) :: momz_(elem%np*lmesh%nea)
83 real(rp),
intent(in) :: drhot_(elem%np*lmesh%nea)
84 real(rp),
intent(in) :: dpres(elem%np*lmesh%nea)
85 real(rp),
intent(in) :: dens_hyd(elem%np*lmesh%nea)
86 real(rp),
intent(in) :: pres_hyd(elem%np*lmesh%nea)
87 real(rp),
intent(in) :: therm_hyd(elem%np*lmesh%nea)
88 real(rp),
intent(in) :: rtot (elem%np*lmesh%nea)
89 real(rp),
intent(in) :: cvtot(elem%np*lmesh%nea)
90 real(rp),
intent(in) :: cptot(elem%np*lmesh%nea)
91 real(rp),
intent(in) :: gsqrt(elem%np*lmesh%nea)
92 real(rp),
intent(in) :: g13(elem%np*lmesh%nea)
93 real(rp),
intent(in) :: g23(elem%np*lmesh%nea)
94 real(rp),
intent(in) :: nx(elem%nfptot,lmesh%ne)
95 real(rp),
intent(in) :: ny(elem%nfptot,lmesh%ne)
96 real(rp),
intent(in) :: nz(elem%nfptot,lmesh%ne)
97 integer,
intent(in) :: vmapm(elem%nfptot,lmesh%ne)
98 integer,
intent(in) :: vmapp(elem%nfptot,lmesh%ne)
100 integer :: ke, fp, i, ip, im
102 real(rp) :: vel_m, vel_p, alpha
103 real(rp) :: dpres_m, dpres_p
104 real(rp) :: gsqrtdens_m, gsqrtdens_p
105 real(rp) :: gsqrtrhot_m, gsqrtrhot_p
106 real(rp) :: gsqrtddens_m, gsqrtddens_p
107 real(rp) :: gsqrtmomx_m, gsqrtmomx_p
108 real(rp) :: gsqrtmomy_m, gsqrtmomy_p
109 real(rp) :: gsqrtmomz_m, gsqrtmomz_p
110 real(rp) :: gsqrtdrhot_m, gsqrtdrhot_p
111 real(rp) :: gsqrt_m, gsqrt_p
112 real(rp) :: gsqrtv_m, gsqrtv_p
113 real(rp) :: rgsqrtv_m, rgsqrtv_p
114 real(rp) :: g13_m, g13_p
115 real(rp) :: g23_m, g23_p
116 real(rp) :: gnn_m, gnn_p
118 real(rp) :: gamm, rgamm
120 real(rp) :: rovp0, p0ovr
123 real(rp) :: del_flux_tmp_momz, del_flux_tmp_momx, del_flux_tmp_momy
124 real(rp) :: gsqrtdpres_m, gsqrtdpres_p
125 real(rp) :: nx_, ny_, nz_, fscale_
127 integer :: nes, nee, nfptot
131 rgamm = cvdry / cpdry
132 rp0 = 1.0_rp / pres00
134 p0ovr = pres00 / rdry
136 nes = lmesh%NeS; nee = lmesh%NeE
148 im = vmapm(fp,ke); ip = vmapp(fp,ke)
149 ke2d = lmesh%EMap3Dto2D(ke)
151 nx_ = nx(fp,ke); ny_ = ny(fp,ke); nz_ = nz(fp,ke)
159 rgsqrtv_m = 1.0_rp / gsqrtv_m
160 rgsqrtv_p = 1.0_rp / gsqrtv_p
167 gsqrtddens_m = gsqrt_m * ddens_(im)
168 gsqrtddens_p = gsqrt_p * ddens_(ip)
169 gsqrtmomx_m = gsqrt_m * momx_(im)
170 gsqrtmomx_p = gsqrt_p * momx_(ip)
171 gsqrtmomy_m = gsqrt_m * momy_(im)
172 gsqrtmomy_p = gsqrt_p * momy_(ip)
173 gsqrtmomz_m = gsqrt_m * momz_(im)
174 gsqrtmomz_p = gsqrt_p * momz_(ip)
175 gsqrtdrhot_m = gsqrt_m * drhot_(im)
176 gsqrtdrhot_p = gsqrt_p * drhot_(ip)
181 gsqrtdens_p = gsqrtddens_p + gsqrt_p * dens_hyd(ip)
182 gsqrtdens_m = gsqrtddens_m + gsqrt_m * dens_hyd(im)
184 gsqrtrhot_p = gsqrt_p * therm_hyd(ip) + gsqrtdrhot_p
185 gsqrtrhot_m = gsqrt_m * therm_hyd(im) + gsqrtdrhot_m
187 vel_m = ( gsqrtmomx_m * nx_ + gsqrtmomy_m * ny_ &
188 + ( ( gsqrtmomz_m * rgsqrtv_m &
189 + g13_m * gsqrtmomx_m + g23_m * gsqrtmomy_m ) * nz_ ) &
192 vel_p = ( gsqrtmomx_p * nx_ + gsqrtmomy_p * ny_ &
193 + ( ( gsqrtmomz_p * rgsqrtv_p &
194 + g13_p * gsqrtmomx_p + g23_p * gsqrtmomy_p ) * nz_ ) &
200 tmp1 = abs( nx_ ) + abs( ny_ )
202 + ( 1.0_rp * rgsqrtv_m**2 + g13_m**2 + g23_m**2 ) * abs( nz_ )
205 + ( 1.0_rp * rgsqrtv_p**2 + g13_p**2 + g23_p**2 ) * abs( nz_ )
207 alpha = max( sqrt( gnn_m * gamm * ( pres_hyd(im) + dpres_m ) * gsqrt_m / gsqrtdens_m ) + abs(vel_m), &
208 sqrt( gnn_p * gamm * ( pres_hyd(ip) + dpres_p ) * gsqrt_p / gsqrtdens_p ) + abs(vel_p) )
214 fscale_ = 0.5_rp * lmesh%Fscale(fp,ke)
218 del_flux(fp,dens_vid,ke) = fscale_ * ( &
219 gsqrtdens_p * vel_p - gsqrtdens_m * vel_m &
220 - alpha * ( gsqrtddens_p - gsqrtddens_m ) )
222 del_flux(fp,rhot_vid,ke) = fscale_ * ( &
223 gsqrtrhot_p * vel_p - gsqrtrhot_m * vel_m &
224 - alpha * ( gsqrtdrhot_p - gsqrtdrhot_m ) )
227 gsqrtdpres_m = gsqrt_m * dpres_m
228 gsqrtdpres_p = gsqrt_p * dpres_p
230 del_flux(fp,momz_vid,ke) = fscale_ * ( &
231 gsqrtmomz_p * vel_p - gsqrtmomz_m * vel_m &
232 + ( gsqrtdpres_p * rgsqrtv_p &
233 - gsqrtdpres_m * rgsqrtv_m ) * nz_ &
234 - alpha * ( gsqrtmomz_p - gsqrtmomz_m ) )
236 del_flux(fp,momx_vid,ke) = fscale_ * ( &
237 gsqrtmomx_p * vel_p - gsqrtmomx_m * vel_m &
238 + ( nx_ + g13_p * nz_ ) * gsqrtdpres_p &
239 - ( nx_ + g13_m * nz_ ) * gsqrtdpres_m &
240 - alpha * ( gsqrtmomx_p - gsqrtmomx_m ) )
242 del_flux(fp,momy_vid,ke) = fscale_ * ( &
243 gsqrtmomy_p * vel_p - gsqrtmomy_m * vel_m &
244 + ( ny_ + g23_p * nz_ ) * gsqrtdpres_p &
245 - ( ny_ + g23_m * nz_ ) * gsqrtdpres_m &
246 - alpha * ( gsqrtmomy_p - gsqrtmomy_m ) )
256 ddens_, momx_, momy_, momz_, drhot_, dpres, &
257 dens_hyd, pres_hyd, therm_hyd, rtot, cvtot, cptot, &
258 gsqrt, g11, g12, g22, gsqrth, gam, g13, g23, nx, ny, nz, &
259 vmapm, vmapp, im2dto3d, lmesh, elem, lmesh2d, elem2d )
267 real(rp),
intent(out) :: del_flux(elem%nfptot,
prgvar_num,lmesh%ne)
268 real(rp),
intent(in) :: ddens_(elem%np*lmesh%nea)
269 real(rp),
intent(in) :: momx_(elem%np*lmesh%nea)
270 real(rp),
intent(in) :: momy_(elem%np*lmesh%nea)
271 real(rp),
intent(in) :: momz_(elem%np*lmesh%nea)
272 real(rp),
intent(in) :: drhot_(elem%np*lmesh%nea)
273 real(rp),
intent(in) :: dpres(elem%np*lmesh%nea)
274 real(rp),
intent(in) :: dens_hyd(elem%np*lmesh%nea)
275 real(rp),
intent(in) :: pres_hyd(elem%np*lmesh%nea)
276 real(rp),
intent(in) :: therm_hyd(elem%np*lmesh%nea)
277 real(rp),
intent(in) :: rtot (elem%np*lmesh%nea)
278 real(rp),
intent(in) :: cvtot(elem%np*lmesh%nea)
279 real(rp),
intent(in) :: cptot(elem%np*lmesh%nea)
280 real(rp),
intent(in) :: gsqrt(elem%np*lmesh%nea)
281 real(rp),
intent(in) :: g11(elem2d%np,lmesh2d%ne)
282 real(rp),
intent(in) :: g12(elem2d%np,lmesh2d%ne)
283 real(rp),
intent(in) :: g22(elem2d%np,lmesh2d%ne)
284 real(rp),
intent(in) :: gsqrth(elem2d%np,lmesh2d%ne)
285 real(rp),
intent(in) :: gam(elem%np*lmesh%nea)
286 real(rp),
intent(in) :: g13(elem%np*lmesh%nea)
287 real(rp),
intent(in) :: g23(elem%np*lmesh%nea)
288 real(rp),
intent(in) :: nx(elem%nfptot,lmesh%ne)
289 real(rp),
intent(in) :: ny(elem%nfptot,lmesh%ne)
290 real(rp),
intent(in) :: nz(elem%nfptot,lmesh%ne)
291 integer,
intent(in) :: vmapm(elem%nfptot,lmesh%ne)
292 integer,
intent(in) :: vmapp(elem%nfptot,lmesh%ne)
293 integer,
intent(in) :: im2dto3d(elem%nfptot)
295 integer :: ke, fp, ip, im
297 real(rp) :: vel_m, vel_p, alpha
298 real(rp) :: dpres_m, dpres_p
299 real(rp) :: gsqrtdens_m, gsqrtdens_p
300 real(rp) :: gsqrtrhot_m, gsqrtrhot_p
301 real(rp) :: gsqrtddens_m, gsqrtddens_p
302 real(rp) :: gsqrtmomx_m, gsqrtmomx_p
303 real(rp) :: gsqrtmomy_m, gsqrtmomy_p
304 real(rp) :: gsqrtmomz_m, gsqrtmomz_p
305 real(rp) :: gsqrtdrhot_m, gsqrtdrhot_p
306 real(rp) :: gsqrt_m, gsqrt_p
307 real(rp) :: gsqrtv_m, gsqrtv_p
308 real(rp) :: rgsqrtv_m, rgsqrtv_p
309 real(rp) :: g13_m, g13_p
310 real(rp) :: g23_m, g23_p
311 real(rp) :: gnn_m, gnn_p
312 real(rp) :: rgam2_m, rgam2_p
317 real(rp) :: gsqrtdpres_m, gsqrtdpres_p
318 real(rp) :: nx_, ny_, nz_, fscale_
320 integer :: nes, nee, nfptot
321 real(rp) :: g11_, g12_, g22_
322 real(rp) :: gxz_m, gxz_p
323 real(rp) :: gyz_m, gyz_p
324 real(rp) :: g1n_m, g1n_p
325 real(rp) :: g2n_m, g2n_p
330 nes = lmesh%NeS; nee = lmesh%NeE
344 im = vmapm(fp,ke); ip = vmapp(fp,ke)
345 ke2d = lmesh%EMap3Dto2D(ke)
347 nx_ = nx(fp,ke); ny_ = ny(fp,ke); nz_ = nz(fp,ke)
354 rgam2_m = 1.0_rp / gam(im)**2
355 gsqrtv_m = gsqrt_m * rgam2_m / gsqrth(im2dto3d(fp),ke2d)
356 rgsqrtv_m = 1.0_rp / gsqrtv_m
358 rgam2_p = 1.0_rp / gam(ip)**2
359 gsqrtv_p = gsqrt_p * rgam2_p / gsqrth(im2dto3d(fp),ke2d)
360 rgsqrtv_p = 1.0_rp / gsqrtv_p
368 gsqrtddens_m = gsqrt_m * ddens_(im)
369 gsqrtddens_p = gsqrt_p * ddens_(ip)
370 gsqrtmomx_m = gsqrt_m * momx_(im)
371 gsqrtmomx_p = gsqrt_p * momx_(ip)
372 gsqrtmomy_m = gsqrt_m * momy_(im)
373 gsqrtmomy_p = gsqrt_p * momy_(ip)
374 gsqrtmomz_m = gsqrt_m * momz_(im)
375 gsqrtmomz_p = gsqrt_p * momz_(ip)
376 gsqrtdrhot_m = gsqrt_m * drhot_(im)
377 gsqrtdrhot_p = gsqrt_p * drhot_(ip)
382 gsqrtdens_p = gsqrtddens_p + gsqrt_p * dens_hyd(ip)
383 gsqrtdens_m = gsqrtddens_m + gsqrt_m * dens_hyd(im)
385 gsqrtrhot_p = gsqrt_p * therm_hyd(ip) + gsqrtdrhot_p
386 gsqrtrhot_m = gsqrt_m * therm_hyd(im) + gsqrtdrhot_m
388 vel_m = ( gsqrtmomx_m * nx_ + gsqrtmomy_m * ny_ &
389 + ( ( gsqrtmomz_m * rgsqrtv_m &
390 + g13_m * gsqrtmomx_m + g23_m * gsqrtmomy_m ) * nz_ ) &
393 vel_p = ( gsqrtmomx_p * nx_ + gsqrtmomy_p * ny_ &
394 + ( ( gsqrtmomz_p * rgsqrtv_p &
395 + g13_p * gsqrtmomx_p + g23_p * gsqrtmomy_p ) * nz_ ) &
401 g11_ = g11(im2dto3d(fp),ke2d)
402 g12_ = g12(im2dto3d(fp),ke2d)
403 g22_ = g22(im2dto3d(fp),ke2d)
404 tmp1 = abs( g11_ * nx_ ) + abs( g22_ * ny_ )
406 gxz_m = rgam2_m * ( g11_ * g13_m + g12_ * g23_m )
407 gyz_m = rgam2_m * ( g12_ * g13_m + g22_ * g23_m )
408 g1n_m = rgam2_m * ( g11_ * nx_ + g12_ * ny_ )
409 g2n_m = rgam2_m * ( g12_ * nx_ + g22_ * ny_ )
410 gnn_m = rgam2_m * tmp1 &
411 + ( 1.0_rp * rgsqrtv_m**2 + g13_m * gxz_m + g23_m * gyz_m ) * abs( nz_ )
413 gxz_p = rgam2_p * ( g11_ * g13_p + g12_ * g23_p )
414 gyz_p = rgam2_p * ( g12_ * g13_p + g22_ * g23_p )
415 g1n_p = rgam2_p * ( g11_ * nx_ + g12_ * ny_ )
416 g2n_p = rgam2_p * ( g12_ * nx_ + g22_ * ny_ )
417 gnn_p = rgam2_p * tmp1 &
418 + ( 1.0_rp * rgsqrtv_p**2 + g13_p * gxz_p + g23_p * gyz_p ) * abs( nz_ )
420 alpha = max( sqrt( gnn_m * gamm * ( pres_hyd(im) + dpres_m ) * gsqrt_m / gsqrtdens_m ) + abs(vel_m), &
421 sqrt( gnn_p * gamm * ( pres_hyd(ip) + dpres_p ) * gsqrt_p / gsqrtdens_p ) + abs(vel_p) )
424 fscale_ = 0.5_rp * lmesh%Fscale(fp,ke)
428 del_flux(fp,dens_vid,ke) = fscale_ * ( &
429 gsqrtdens_p * vel_p - gsqrtdens_m * vel_m &
430 - alpha * ( gsqrtddens_p - gsqrtddens_m ) )
432 del_flux(fp,rhot_vid,ke) = fscale_ * ( &
433 gsqrtrhot_p * vel_p - gsqrtrhot_m * vel_m &
434 - alpha * ( gsqrtdrhot_p - gsqrtdrhot_m ) )
437 gsqrtdpres_m = gsqrt_m * dpres_m
438 gsqrtdpres_p = gsqrt_p * dpres_p
440 del_flux(fp,momz_vid,ke) = fscale_ * ( &
441 gsqrtmomz_p * vel_p - gsqrtmomz_m * vel_m &
442 + ( gsqrtdpres_p * rgsqrtv_p &
443 - gsqrtdpres_m * rgsqrtv_m ) * nz_ &
444 - alpha * ( gsqrtmomz_p - gsqrtmomz_m ) )
446 del_flux(fp,momx_vid,ke) = fscale_ * ( &
447 gsqrtmomx_p * vel_p - gsqrtmomx_m * vel_m &
448 + ( g1n_p + gxz_p * nz_ ) * gsqrtdpres_p &
449 - ( g1n_m + gxz_m * nz_ ) * gsqrtdpres_m &
450 - alpha * ( gsqrtmomx_p - gsqrtmomx_m ) )
452 del_flux(fp,momy_vid,ke) = fscale_ * ( &
453 gsqrtmomy_p * vel_p - gsqrtmomy_m * vel_m &
454 + ( g2n_p + gyz_p * nz_ ) * gsqrtdpres_p &
455 - ( g2n_m + gyz_m * nz_ ) * gsqrtdpres_m &
456 - alpha * ( gsqrtmomy_p - gsqrtmomy_m ) )