18#include "scaleFElib.h"
68 private :: atm_dyn_dgm_trcadvect3d_heve_cal_alphdens_dyn
69 private :: atm_dyn_dgm_trcadvect3d_heve_get_netoutwardflux_generalhvc
70 private :: atm_dyn_dgm_trcadvect3d_heve_get_delflux_generalhvc
79 type(
sparsemat),
intent(inout) :: faceintmat
83 real(rp),
allocatable :: intweight_lgl1dpts_h(:)
84 real(rp),
allocatable :: intweight_lgl1dpts_v(:)
85 real(rp),
allocatable :: intweight_h(:)
86 real(rp),
allocatable :: intweight_v(:)
92 real(rp),
allocatable :: intweight(:,:)
95 lcmesh => mesh%lcmesh_list(1)
96 elem => lcmesh%refElem3D
97 allocate( intweight(elem%Nfaces,elem%NfpTot) )
98 intweight(:,:) = 0.0_rp
100 allocate( intweight_lgl1dpts_h(elem%Nnode_h1D) )
101 allocate( intweight_lgl1dpts_v(elem%Nnode_v) )
102 allocate( intweight_h(elem%Nnode_h1D*elem%Nnode_v) )
103 allocate( intweight_v(elem%Nnode_h1D**2) )
108 do f=1, elem%Nfaces_h
110 do i=1, elem%Nnode_h1D
111 l = i + (k-1)*elem%Nnode_h1D
112 intweight_h(l) = intweight_lgl1dpts_h(i) * intweight_lgl1dpts_v(k)
116 is = (f-1)*elem%Nfp_h + 1
117 ie = is + elem%Nfp_h - 1
118 intweight(f,is:ie) = intweight_h(:)
121 do f=1, elem%Nfaces_v
122 do j=1, elem%Nnode_h1D
123 do i=1, elem%Nnode_h1D
124 l = i + (j-1)*elem%Nnode_h1D
125 intweight_v(l) = intweight_lgl1dpts_h(i) * intweight_lgl1dpts_h(j)
129 is = elem%Nfaces_h*elem%Nfp_h + (f-1)*elem%Nfp_v + 1
130 ie = is + elem%Nfp_v - 1
131 intweight(elem%Nfaces_h+f,is:ie) = intweight_v(:)
134 call faceintmat%Init( intweight )
151 qtrc_, momx_, momy_, momz_, &
152 alphdens_m, alphdens_p, fct_coef, &
154 element3d_operation, faceintmat, &
155 lmesh, elem, lmesh2d, elem2d )
162 real(rp),
intent(out) :: qtrc_dt(elem%np,lmesh%nea)
163 real(rp),
intent(in) :: qtrc_(elem%np,lmesh%nea)
164 real(rp),
intent(in) :: momx_(elem%np,lmesh%nea)
165 real(rp),
intent(in) :: momy_(elem%np,lmesh%nea)
166 real(rp),
intent(in) :: momz_(elem%np,lmesh%nea)
167 real(rp),
intent(in) :: alphdens_m(elem%nfptot,lmesh%ne)
168 real(rp),
intent(in) :: alphdens_p(elem%nfptot,lmesh%ne)
169 real(rp),
intent(in) :: fct_coef(elem%np,lmesh%nea)
170 real(rp),
intent(in) :: rhoq_tp(elem%np,lmesh%nea)
172 type(
sparsemat),
intent(in) :: faceintmat
174 real(rp) :: flux(elem%np,3), dflux(elem%np,4)
175 real(rp) :: del_flux(elem%nfptot,lmesh%ne)
178 real(rp) :: rgsqrt(elem%np), rgsqrtv(elem%np)
182 real(rp) :: q0, q1, vol
186 call prof_rapstart(
'cal_trcadv_tend_bndflux', 3)
187 call atm_dyn_dgm_trcadvect3d_heve_get_delflux_generalhvc( &
189 qtrc_, momx_, momy_, momz_, alphdens_m, alphdens_p, fct_coef, &
190 lmesh%Gsqrt, lmesh%GsqrtH, lmesh%GI3(:,:,1), lmesh%GI3(:,:,2), &
191 lmesh%normal_fn(:,:,1), lmesh%normal_fn(:,:,2), lmesh%normal_fn(:,:,3), &
192 lmesh%J(:,:), lmesh%Fscale(:,:), &
193 lmesh%vmapM, lmesh%vmapP, elem%IndexH2Dto3D_bnd, faceintmat, &
194 lmesh, elem, lmesh2d, elem2d )
195 call prof_rapend(
'cal_trcadv_tend_bndflux', 3)
197 call prof_rapstart(
'cal_trcadv_tend_interior', 3)
200 do ke=lmesh%NeS, lmesh%NeE
201 ke2d = lmesh%EMap3Dto2D(ke)
204 rgsqrt(p) = 1.0_rp / lmesh%Gsqrt(p,ke)
205 rgsqrtv(p) = lmesh%GsqrtH(elem%IndexH2Dto3D(p),ke2d) * rgsqrt(p)
209 gsqrt_ = lmesh%Gsqrt(p,ke)
210 flux(p,1) = gsqrt_ * momx_(p,ke) * qtrc_(p,ke)
211 flux(p,2) = gsqrt_ * momy_(p,ke) * qtrc_(p,ke)
212 flux(p,3) = gsqrt_ * ( &
213 momz_(p,ke) * rgsqrtv(p) &
214 + lmesh%GI3(p,ke,1) * momx_(p,ke) &
215 + lmesh%GI3(p,ke,2) * momy_(p,ke) ) * qtrc_(p,ke)
217 call element3d_operation%Div( &
218 flux, del_flux(:,ke), &
221 qtrc_dt(:,ke) = - ( &
222 lmesh%Escale(:,ke,1,1) * dflux(:,1) &
223 + lmesh%Escale(:,ke,2,2) * dflux(:,2) &
224 + lmesh%Escale(:,ke,3,3) * dflux(:,3) &
225 + dflux(:,4) ) / lmesh%Gsqrt(:,ke) &
228 call prof_rapend(
'cal_trcadv_tend_interior', 3)
236 qtrc_, momx_, momy_, momz_, rhoq_tp_, alphdens_m, alphdens_p, &
237 dens_hyd, ddens_, ddens0_, rk_c_ssm1, dt, &
238 faceintmat, lmesh, elem, lmesh2d, elem2d, &
246 real(rp),
intent(out) :: fct_coef(elem%np,lmesh%nea)
247 real(rp),
intent(in) :: qtrc_(elem%np,lmesh%nea)
248 real(rp),
intent(in) :: momx_(elem%np,lmesh%nea)
249 real(rp),
intent(in) :: momy_(elem%np,lmesh%nea)
250 real(rp),
intent(in) :: momz_(elem%np,lmesh%nea)
251 real(rp),
intent(in) :: rhoq_tp_(elem%np,lmesh%nea)
252 real(rp),
intent(in) :: alphdens_m(elem%nfptot,lmesh%ne)
253 real(rp),
intent(in) :: alphdens_p(elem%nfptot,lmesh%ne)
254 real(rp),
intent(in) :: dens_hyd(elem%np,lmesh%nea)
255 real(rp),
intent(in) :: ddens_ (elem%np,lmesh%nea)
256 real(rp),
intent(in) :: ddens0_(elem%np,lmesh%nea)
257 real(rp),
intent(in) :: rk_c_ssm1
258 real(rp),
intent(in) :: dt
259 type(
sparsemat),
intent(in) :: faceintmat
260 logical,
intent(in),
optional :: disable_limiter
262 real(rp) :: netoutwardflux(lmesh%ne)
263 real(rp) :: momwt_(elem%np)
267 real(rp) :: dens_ssm1(elem%np)
270 if (
present(disable_limiter) )
then
271 if ( disable_limiter )
then
273 do ke=lmesh%NeS, lmesh%NeE
274 fct_coef(:,ke) = 1.0_rp
280 call prof_rapstart(
'cal_trcadv_fct_coef_bndflux', 3)
281 call atm_dyn_dgm_trcadvect3d_heve_get_netoutwardflux_generalhvc( &
283 qtrc_, momx_, momy_, momz_, alphdens_m, alphdens_p, &
284 lmesh%Gsqrt, lmesh%GsqrtH, lmesh%GI3(:,:,1), lmesh%GI3(:,:,2), &
285 lmesh%normal_fn(:,:,1), lmesh%normal_fn(:,:,2), lmesh%normal_fn(:,:,3), &
286 lmesh%J(:,:), lmesh%Fscale(:,:), lmesh%vmapM, lmesh%vmapP, elem%IndexH2Dto3D_bnd, &
288 lmesh, elem, lmesh2d, elem2d )
289 call prof_rapend(
'cal_trcadv_fct_coef_bndflux', 3)
291 call prof_rapstart(
'cal_trcadv_fct_coef', 3)
294 do ke=lmesh%NeS, lmesh%NeE
296 dens_ssm1(:) = dens_hyd(:,ke) &
297 + ( 1.0_rp - rk_c_ssm1 ) * ddens0_(:,ke) + rk_c_ssm1 * ddens_(:,ke)
298 q = sum( lmesh%Gsqrt(:,ke) * lmesh%J(:,ke) * elem%IntWeight_lgl(:) * ( dens_ssm1(:) * qtrc_(:,ke) / dt + rhoq_tp_(:,ke) ) )
300 fct_coef(:,ke) = max( 0.0_rp, min( 1.0_rp, q / ( netoutwardflux(ke) + 1.0e-10_rp ) ) )
303 call prof_rapend(
'cal_trcadv_fct_coef', 3)
312 dens_hyd, ddens_, lmesh, elem, lmesh2d, elem2d )
319 real(rp),
intent(inout) :: qtrc_(elem%np,lmesh%nea)
320 real(rp),
intent(in) :: dens_hyd(elem%np,lmesh%nea)
321 real(rp),
intent(in) :: ddens_ (elem%np,lmesh%nea)
325 real(rp) :: q(elem%np)
326 real(rp) :: dens(elem%np)
330 do ke=lmesh%NeS, lmesh%NeE
331 dens(:) = dens_hyd(:,ke) + ddens_(:,ke)
332 q(:) = max( 0.0_rp, qtrc_(:,ke) )
334 q0 = sum( lmesh%Gsqrt(:,ke) * lmesh%J(:,ke) * elem%IntWeight_lgl(:) * dens(:) * qtrc_(:,ke) )
335 q1 = sum( lmesh%Gsqrt(:,ke) * lmesh%J(:,ke) * elem%IntWeight_lgl(:) * dens(:) * q(:) )
336 qtrc_(:,ke) = q0 / ( q1 + 1.0e-32_rp ) * q(:)
344 MFLX_x_tavg, MFLX_y_tavg, MFLX_z_tavg, alph_dens_M, alph_dens_P, &
345 DDENS, MOMX, MOMY, MOMZ, DPRES, DENS_hyd, PRES_hyd, &
346 Rtot, CVtot, CPtot, &
347 lmesh, elem, rkstage, tavg_weight_h, tavg_weight_v, is_hevi )
352 real(rp),
intent(inout) :: mflx_x_tavg(elem%np,lmesh%nea)
353 real(rp),
intent(inout) :: mflx_y_tavg(elem%np,lmesh%nea)
354 real(rp),
intent(inout) :: mflx_z_tavg(elem%np,lmesh%nea)
355 real(rp),
intent(inout) :: alph_dens_m(elem%nfptot,lmesh%ne)
356 real(rp),
intent(inout) :: alph_dens_p(elem%nfptot,lmesh%ne)
357 real(rp),
intent(in) :: ddens(elem%np,lmesh%nea)
358 real(rp),
intent(in) :: momx(elem%np,lmesh%nea)
359 real(rp),
intent(in) :: momy(elem%np,lmesh%nea)
360 real(rp),
intent(in) :: momz(elem%np,lmesh%nea)
361 real(rp),
intent(in) :: dpres(elem%np,lmesh%nea)
362 real(rp),
intent(in) :: dens_hyd(elem%np,lmesh%nea)
363 real(rp),
intent(in) :: pres_hyd(elem%np,lmesh%nea)
364 real(rp),
intent(in) :: rtot(elem%np,lmesh%nea)
365 real(rp),
intent(in) :: cvtot(elem%np,lmesh%nea)
366 real(rp),
intent(in) :: cptot(elem%np,lmesh%nea)
367 integer,
intent(in) :: rkstage
368 real(rp),
intent(in) :: tavg_weight_h
369 real(rp),
intent(in) :: tavg_weight_v
370 logical,
intent(in) :: is_hevi
376 do ke=lmesh%NeS, lmesh%NeE
377 if (rkstage == 1)
then
378 mflx_x_tavg(:,ke) = tavg_weight_h * momx(:,ke)
379 mflx_y_tavg(:,ke) = tavg_weight_h * momy(:,ke)
380 mflx_z_tavg(:,ke) = tavg_weight_v * momz(:,ke)
381 alph_dens_m(:,ke) = 0.0_rp
382 alph_dens_p(:,ke) = 0.0_rp
384 mflx_x_tavg(:,ke) = mflx_x_tavg(:,ke) + tavg_weight_h * momx(:,ke)
385 mflx_y_tavg(:,ke) = mflx_y_tavg(:,ke) + tavg_weight_h * momy(:,ke)
386 mflx_z_tavg(:,ke) = mflx_z_tavg(:,ke) + tavg_weight_v * momz(:,ke)
390 call atm_dyn_dgm_trcadvect3d_heve_cal_alphdens_dyn( &
391 alph_dens_m, alph_dens_p, &
392 ddens, momx, momy, momz, dpres, dens_hyd, pres_hyd, &
393 rtot, cvtot, cptot, &
394 lmesh%Gsqrt, lmesh%GIJ(:,:,1,1), lmesh%GIJ(:,:,1,2), lmesh%GIJ(:,:,2,2), &
395 lmesh%normal_fn(:,:,1), lmesh%normal_fn(:,:,2), lmesh%normal_fn(:,:,3), &
396 lmesh%vmapM, lmesh%vmapP, elem%IndexH2Dto3D_bnd, &
397 lmesh, lmesh%lcmesh2D, lmesh%refElem3D, lmesh%lcmesh2D%refElem2D, &
398 tavg_weight_h, tavg_weight_v, is_hevi )
405 DDENS_, MOMX_, MOMY_, MOMZ_, DENS_hyd, &
406 Gsqrt, nx, ny, nz, vmapM, vmapP, lmesh, elem )
408 use scale_const,
only: &
409 grav => const_grav, &
410 rdry => const_rdry, &
411 cpdry => const_cpdry, &
412 cvdry => const_cvdry, &
413 pres00 => const_pre00
419 real(rp),
intent(inout) :: alph_dens_m(elem%nfptot*lmesh%ne)
420 real(rp),
intent(inout) :: alph_dens_p(elem%nfptot*lmesh%ne)
421 real(rp),
intent(in) :: ddens_(elem%np*lmesh%nea)
422 real(rp),
intent(in) :: momx_(elem%np*lmesh%nea)
423 real(rp),
intent(in) :: momy_(elem%np*lmesh%nea)
424 real(rp),
intent(in) :: momz_(elem%np*lmesh%nea)
425 real(rp),
intent(in) :: dens_hyd(elem%np*lmesh%nea)
426 real(rp),
intent(in) :: gsqrt(elem%np*lmesh%nea)
427 real(rp),
intent(in) :: nx(elem%nfptot*lmesh%ne)
428 real(rp),
intent(in) :: ny(elem%nfptot*lmesh%ne)
429 real(rp),
intent(in) :: nz(elem%nfptot*lmesh%ne)
430 integer,
intent(in) :: vmapm(elem%nfptot*lmesh%ne)
431 integer,
intent(in) :: vmapp(elem%nfptot*lmesh%ne)
434 real(rp) :: velp, velm, alpha, densm, densp
440 do i=1, elem%NfpTot*lmesh%Ne
441 im = vmapm(i); ip = vmapp(i)
443 densm = ddens_(im) + dens_hyd(im)
444 densp = ddens_(ip) + dens_hyd(ip)
446 velm = ( momx_(im) * nx(i) + momy_(im) * ny(i) + momz_(im) * nz(i) ) / densm
447 velp = ( momx_(ip) * nx(i) + momy_(ip) * ny(i) + momz_(ip) * nz(i) ) / densp
449 alpha = max( abs(velm), abs(velp) )
450 alph_dens_m(i) = alpha * densm * gsqrt(im)
451 alph_dens_p(i) = alpha * densp * gsqrt(ip)
460 subroutine atm_dyn_dgm_trcadvect3d_heve_cal_alphdens_dyn( alph_dens_M, alph_dens_P, &
461 DDENS_, MOMX_, MOMY_, MOMZ_, DPRES_, DENS_hyd, PRES_hyd, &
462 Rtot, CVtot, CPtot, &
463 Gsqrt, G11, G12, G22, nx, ny, nz, vmapM, vmapP, iM2Dto3D, &
464 lmesh, lmesh2D, elem, elem2D, &
465 tavg_weight_h, tavg_weight_v, is_hevi )
467 use scale_const,
only: &
468 grav => const_grav, &
469 rdry => const_rdry, &
470 cpdry => const_cpdry, &
471 cvdry => const_cvdry, &
472 pres00 => const_pre00
480 real(rp),
intent(inout) :: alph_dens_m(elem%nfptot,lmesh%ne)
481 real(rp),
intent(inout) :: alph_dens_p(elem%nfptot,lmesh%ne)
482 real(rp),
intent(in) :: ddens_(elem%np*lmesh%nea)
483 real(rp),
intent(in) :: momx_(elem%np*lmesh%nea)
484 real(rp),
intent(in) :: momy_(elem%np*lmesh%nea)
485 real(rp),
intent(in) :: momz_(elem%np*lmesh%nea)
486 real(rp),
intent(in) :: dpres_(elem%np*lmesh%nea)
487 real(rp),
intent(in) :: dens_hyd(elem%np*lmesh%nea)
488 real(rp),
intent(in) :: pres_hyd(elem%np*lmesh%nea)
489 real(rp),
intent(in) :: rtot(elem%np*lmesh%nea)
490 real(rp),
intent(in) :: cvtot(elem%np*lmesh%nea)
491 real(rp),
intent(in) :: cptot(elem%np*lmesh%nea)
492 real(rp),
intent(in) :: gsqrt(elem%np*lmesh%nea)
493 real(rp),
intent(in) :: g11(elem2d%np,lmesh2d%ne)
494 real(rp),
intent(in) :: g12(elem2d%np,lmesh2d%ne)
495 real(rp),
intent(in) :: g22(elem2d%np,lmesh2d%ne)
496 real(rp),
intent(in) :: nx(elem%nfptot,lmesh%ne)
497 real(rp),
intent(in) :: ny(elem%nfptot,lmesh%ne)
498 real(rp),
intent(in) :: nz(elem%nfptot,lmesh%ne)
499 integer,
intent(in) :: vmapm(elem%nfptot,lmesh%ne)
500 integer,
intent(in) :: vmapp(elem%nfptot,lmesh%ne)
501 integer,
intent(in) :: im2dto3d(elem%nfptot)
502 real(rp),
intent(in) :: tavg_weight_h
503 real(rp),
intent(in) :: tavg_weight_v
504 logical,
intent(in) :: is_hevi
506 integer :: ke, ip(elem%nfptot), im(elem%nfptot)
508 real(rp) :: velp(elem%nfptot), velm(elem%nfptot), alpha(elem%nfptot)
509 real(rp) :: densm(elem%nfptot), densp(elem%nfptot)
510 real(rp) :: gamm, rgamm
511 real(rp) :: tavg_weight(elem%nfptot)
513 real(rp) :: gnn_m(elem%nfptot), gnn_p(elem%nfptot)
522 do ke=lmesh%NeS, lmesh%NeE
523 im(:) = vmapm(:,ke); ip(:) = vmapp(:,ke)
524 ke2d = lmesh%EMap3Dto2D(ke)
526 densm(:) = ddens_(im) + dens_hyd(im)
527 densp(:) = ddens_(ip) + dens_hyd(ip)
530 gnn_m(:) = g11(im2dto3d(:),ke2d) * abs( nx(:,ke) ) + g22(im2dto3d(:),ke2d) * abs( ny(:,ke) )
531 gnn_p(:) = g11(im2dto3d(:),ke2d) * abs( nx(:,ke) ) + g22(im2dto3d(:),ke2d) * abs( ny(:,ke) )
533 gnn_m(:) = g11(im2dto3d(:),ke2d) * abs( nx(:,ke) ) + g22(im2dto3d(:),ke2d) * abs( ny(:,ke) ) &
535 gnn_p(:) = g11(im2dto3d(:),ke2d) * abs( nx(:,ke) ) + g22(im2dto3d(:),ke2d) * abs( ny(:,ke) ) &
539 velm(:) = ( momx_(im)*nx(:,ke) + momy_(im)*ny(:,ke) + momz_(im)*nz(:,ke) ) / densm(:)
540 velp(:) = ( momx_(ip)*nx(:,ke) + momy_(ip)*ny(:,ke) + momz_(ip)*nz(:,ke) ) / densp(:)
542 alpha(:) = max( sqrt( gnn_m(:) * gamm * ( pres_hyd(im) + dpres_(im) ) / densm(:) ) + abs(velm(:)), &
543 sqrt( gnn_p(:) * gamm * ( pres_hyd(ip) + dpres_(ip) ) / densp(:) ) + abs(velp(:)) )
544 tavg_weight = tavg_weight_h * ( abs(nx(:,ke)) + abs(ny(:,ke)) ) + tavg_weight_v * abs(nz(:,ke))
546 alph_dens_m(:,ke) = alph_dens_m(:,ke) + tavg_weight * alpha(:) * densm(:) * gsqrt(im)
547 alph_dens_p(:,ke) = alph_dens_p(:,ke) + tavg_weight * alpha(:) * densp(:) * gsqrt(ip)
551 end subroutine atm_dyn_dgm_trcadvect3d_heve_cal_alphdens_dyn
554 subroutine atm_dyn_dgm_trcadvect3d_heve_get_delflux_generalhvc( &
556 qtrc_, momx_, momy_, momz_, alphdens_m, alphdens_p, fct_coef, &
557 gsqrt, gsqrth, g13, g23, nx, ny, nz, j, fscale, &
558 vmapm, vmapp, im2dto3d, faceintmat, &
559 lmesh, elem, lmesh2d, elem2d )
567 real(rp),
intent(out) :: del_flux(elem%nfptot,lmesh%ne)
568 real(rp),
intent(in) :: qtrc_(elem%np*lmesh%nea)
569 real(rp),
intent(in) :: momx_(elem%np*lmesh%nea)
570 real(rp),
intent(in) :: momy_(elem%np*lmesh%nea)
571 real(rp),
intent(in) :: momz_(elem%np*lmesh%nea)
572 real(rp),
intent(in) :: alphdens_m(elem%nfptot,lmesh%ne)
573 real(rp),
intent(in) :: alphdens_p(elem%nfptot,lmesh%ne)
574 real(rp),
intent(in) :: fct_coef(elem%np*lmesh%nea)
575 real(rp),
intent(in) :: gsqrt(elem%np*lmesh%nea)
576 real(rp),
intent(in) :: gsqrth(elem2d%np,lmesh2d%ne)
577 real(rp),
intent(in) :: g13(elem%np*lmesh%nea)
578 real(rp),
intent(in) :: g23(elem%np*lmesh%nea)
579 real(rp),
intent(in) :: nx(elem%nfptot,lmesh%ne)
580 real(rp),
intent(in) :: ny(elem%nfptot,lmesh%ne)
581 real(rp),
intent(in) :: nz(elem%nfptot,lmesh%ne)
582 real(rp),
intent(in) :: j(elem%np*lmesh%ne)
583 real(rp),
intent(in) :: fscale(elem%nfptot,lmesh%ne)
584 integer,
intent(in) :: vmapm(elem%nfptot,lmesh%ne)
585 integer,
intent(in) :: vmapp(elem%nfptot,lmesh%ne)
586 integer,
intent(in) :: im2dto3d(elem%nfptot)
587 type(
sparsemat),
intent(in) :: faceintmat
589 integer :: ke, i, ip(elem%nfptot), im(elem%nfptot)
591 real(rp) :: momflxp(elem%nfptot), momflxm(elem%nfptot), alpha(elem%nfptot)
592 real(rp) :: qtrc_p(elem%nfptot), qtrc_m(elem%nfptot)
593 real(rp) :: gsqrtmomx_p(elem%nfptot), gsqrtmomx_m(elem%nfptot)
594 real(rp) :: gsqrtmomy_p(elem%nfptot), gsqrtmomy_m(elem%nfptot)
595 real(rp) :: gsqrtmomz_p(elem%nfptot), gsqrtmomz_m(elem%nfptot)
596 real(rp) :: gsqrt_p(elem%nfptot), gsqrt_m(elem%nfptot)
597 real(rp) :: gsqrtv_p(elem%nfptot), gsqrtv_m(elem%nfptot)
598 real(rp) :: g13_m(elem%nfptot), g13_p(elem%nfptot)
599 real(rp) :: g23_m(elem%nfptot), g23_p(elem%nfptot)
601 real(rp) :: r_m(elem%nfptot), r_p(elem%nfptot)
602 real(rp) :: numflux(elem%nfptot)
603 real(rp) :: outward_flux_tmp(elem%nfaces)
604 integer :: f, p, fp, is, ie
613 do ke=lmesh%NeS, lmesh%NeE
614 im(:) = vmapm(:,ke); ip(:) = vmapp(:,ke)
615 ke2d = lmesh%EMap3Dto2D(ke)
617 gsqrt_m(:) = gsqrt(im)
618 gsqrt_p(:) = gsqrt(ip)
619 gsqrtv_m(:) = gsqrt_m(:) / gsqrth(im2dto3d(:),ke2d)
620 gsqrtv_p(:) = gsqrt_p(:) / gsqrth(im2dto3d(:),ke2d)
627 r_m(:) = fct_coef(im)
628 r_p(:) = fct_coef(ip)
630 qtrc_m(:) = qtrc_(im)
631 qtrc_p(:) = qtrc_(ip)
632 gsqrtmomx_m(:) = gsqrt_m(:) * momx_(im)
633 gsqrtmomx_p(:) = gsqrt_p(:) * momx_(ip)
634 gsqrtmomy_m(:) = gsqrt_m(:) * momy_(im)
635 gsqrtmomy_p(:) = gsqrt_p(:) * momy_(ip)
636 gsqrtmomz_m(:) = gsqrt_m(:) * momz_(im)
637 gsqrtmomz_p(:) = gsqrt_p(:) * momz_(ip)
639 momflxm(:) = ( gsqrtmomx_m(:) * nx(:,ke) + gsqrtmomy_m(:) * ny(:,ke) &
640 + ( ( gsqrtmomz_m(:) / gsqrtv_m(:) &
641 + g13_m(:) * gsqrtmomx_m(:) + g23_m(:) * gsqrtmomy_m(:) ) * nz(:,ke) ) &
643 momflxp(:) = ( gsqrtmomx_p(:) * nx(:,ke) + gsqrtmomy_p(:) * ny(:,ke) &
644 + ( ( gsqrtmomz_p(:) / gsqrtv_p(:) &
645 + g13_p(:) * gsqrtmomx_p(:) + g23_p(:) * gsqrtmomy_p(:) ) * nz(:,ke) ) &
648 numflux(:) = 0.5_rp * ( ( qtrc_p(:) * momflxp(:) + qtrc_m(:) * momflxm(:) ) &
649 - alphdens_p(:,ke) * qtrc_p(:) + alphdens_m(:,ke) * qtrc_m(:) )
654 call sparsemat_matmul( faceintmat, j(im) * fscale(:,ke) * numflux(:), outward_flux_tmp )
655 do f=1, elem%Nfaces_h
657 fp = p + (f-1)*elem%Nfp_h
658 del_flux(fp,ke) = lmesh%Fscale(fp,ke) * &
659 ( numflux(fp) * 0.5_rp * ( r_p(fp) + r_m(fp) - ( r_p(fp) - r_m(fp) ) * sign( 1.0_rp, outward_flux_tmp(f) ) ) &
660 - qtrc_m(fp) * momflxm(fp) )
663 do f=1, elem%Nfaces_v
665 fp = p + (f-1)*elem%Nfp_v + elem%Nfaces_h * elem%Nfp_h
666 del_flux(fp,ke) = lmesh%Fscale(fp,ke) * &
667 ( numflux(fp) * 0.5_rp * ( r_p(fp) + r_m(fp) - ( r_p(fp) - r_m(fp) ) * sign( 1.0_rp, outward_flux_tmp(elem%Nfaces_h+f) ) ) &
668 - qtrc_m(fp) * momflxm(fp) )
674 end subroutine atm_dyn_dgm_trcadvect3d_heve_get_delflux_generalhvc
678 subroutine atm_dyn_dgm_trcadvect3d_heve_get_netoutwardflux_generalhvc( &
679 net_outward_flux, & ! (out)
680 qtrc_, momx_, momy_, momz_, alphdens_m, alphdens_p, &
681 gsqrt, gsqrth, g13, g23, nx, ny, nz, j, fscale, &
682 vmapm, vmapp, im2dto3d, faceintmat, &
683 lmesh, elem, lmesh2d, elem2d )
691 real(rp),
intent(out) :: net_outward_flux(lmesh%ne)
692 real(rp),
intent(in) :: qtrc_(elem%np*lmesh%nea)
693 real(rp),
intent(in) :: momx_(elem%np*lmesh%nea)
694 real(rp),
intent(in) :: momy_(elem%np*lmesh%nea)
695 real(rp),
intent(in) :: momz_(elem%np*lmesh%nea)
696 real(rp),
intent(in) :: alphdens_m(elem%nfptot,lmesh%ne)
697 real(rp),
intent(in) :: alphdens_p(elem%nfptot,lmesh%ne)
698 real(rp),
intent(in) :: gsqrt(elem%np*lmesh%nea)
699 real(rp),
intent(in) :: gsqrth(elem2d%np,lmesh2d%ne)
700 real(rp),
intent(in) :: g13(elem%np*lmesh%nea)
701 real(rp),
intent(in) :: g23(elem%np*lmesh%nea)
702 real(rp),
intent(in) :: nx(elem%nfptot,lmesh%ne)
703 real(rp),
intent(in) :: ny(elem%nfptot,lmesh%ne)
704 real(rp),
intent(in) :: nz(elem%nfptot,lmesh%ne)
705 real(rp),
intent(in) :: j(elem%np*lmesh%ne)
706 real(rp),
intent(in) :: fscale(elem%nfptot,lmesh%ne)
707 integer,
intent(in) :: vmapm(elem%nfptot,lmesh%ne)
708 integer,
intent(in) :: vmapp(elem%nfptot,lmesh%ne)
709 integer,
intent(in) :: im2dto3d(elem%nfptot)
710 type(
sparsemat),
intent(in) :: faceintmat
712 integer :: ke, i, ip(elem%nfptot), im(elem%nfptot)
714 real(rp) :: momflxp(elem%nfptot), momflxm(elem%nfptot), alpha(elem%nfptot)
715 real(rp) :: qtrc_p(elem%nfptot), qtrc_m(elem%nfptot)
716 real(rp) :: gsqrtmomx_p(elem%nfptot), gsqrtmomx_m(elem%nfptot)
717 real(rp) :: gsqrtmomy_p(elem%nfptot), gsqrtmomy_m(elem%nfptot)
718 real(rp) :: gsqrtmomz_p(elem%nfptot), gsqrtmomz_m(elem%nfptot)
719 real(rp) :: gsqrt_p(elem%nfptot), gsqrt_m(elem%nfptot)
720 real(rp) :: gsqrtv_p(elem%nfptot), gsqrtv_m(elem%nfptot)
721 real(rp) :: g13_m(elem%nfptot), g13_p(elem%nfptot)
722 real(rp) :: g23_m(elem%nfptot), g23_p(elem%nfptot)
724 real(rp) :: numflux(elem%nfptot)
725 real(rp) :: outward_flux_tmp(elem%nfaces)
733 do ke=lmesh%NeS, lmesh%NeE
734 im(:) = vmapm(:,ke); ip(:) = vmapp(:,ke)
735 ke2d = lmesh%EMap3Dto2D(ke)
737 gsqrt_m(:) = gsqrt(im)
738 gsqrt_p(:) = gsqrt(ip)
739 gsqrtv_m(:) = gsqrt_m(:) / gsqrth(im2dto3d(:),ke2d)
740 gsqrtv_p(:) = gsqrt_p(:) / gsqrth(im2dto3d(:),ke2d)
747 qtrc_m(:) = qtrc_(im)
748 qtrc_p(:) = qtrc_(ip)
749 gsqrtmomx_m(:) = gsqrt_m(:) * momx_(im)
750 gsqrtmomx_p(:) = gsqrt_p(:) * momx_(ip)
751 gsqrtmomy_m(:) = gsqrt_m(:) * momy_(im)
752 gsqrtmomy_p(:) = gsqrt_p(:) * momy_(ip)
753 gsqrtmomz_m(:) = gsqrt_m(:) * momz_(im)
754 gsqrtmomz_p(:) = gsqrt_p(:) * momz_(ip)
756 momflxm(:) = ( gsqrtmomx_m(:) * nx(:,ke) + gsqrtmomy_m(:) * ny(:,ke) &
757 + ( ( gsqrtmomz_m(:) / gsqrtv_m(:) &
758 + g13_m(:) * gsqrtmomx_m(:) + g23_m(:) * gsqrtmomy_m(:) ) * nz(:,ke) ) &
760 momflxp(:) = ( gsqrtmomx_p(:) * nx(:,ke) + gsqrtmomy_p(:) * ny(:,ke) &
761 + ( ( gsqrtmomz_p(:) / gsqrtv_p(:) &
762 + g13_p(:) * gsqrtmomx_p(:) + g23_p(:) * gsqrtmomy_p(:) ) * nz(:,ke) ) &
765 alpha(:) = max( abs(momflxm(:)), abs(momflxp(:)) )
767 numflux(:) = 0.5_rp * ( ( qtrc_p(:) * momflxp(:) + qtrc_m(:) * momflxm(:) ) &
768 - alphdens_p(:,ke) * qtrc_p(:) + alphdens_m(:,ke) * qtrc_m(:) )
772 call sparsemat_matmul( faceintmat, j(im) * fscale(:,ke) * numflux(:), outward_flux_tmp )
773 net_outward_flux(ke) = sum( max( 0.0_rp, outward_flux_tmp(:) ) )
777 end subroutine atm_dyn_dgm_trcadvect3d_heve_get_netoutwardflux_generalhvc
module FElib / Fluid dyn solver / Atmosphere / Tracer advection
subroutine, public atm_dyn_dgm_trcadvect3d_heve_final()
subroutine, public atm_dyn_dgm_trcadvect3d_tmar(qtrc_, dens_hyd, ddens_, lmesh, elem, lmesh2d, elem2d)
Second Step of limiter in which nonlinear truncation and mass aware rescaling (TMAR)
subroutine, public atm_dyn_dgm_trcadvect3d_save_massflux(mflx_x_tavg, mflx_y_tavg, mflx_z_tavg, alph_dens_m, alph_dens_p, ddens, momx, momy, momz, dpres, dens_hyd, pres_hyd, rtot, cvtot, cptot, lmesh, elem, rkstage, tavg_weight_h, tavg_weight_v, is_hevi)
subroutine, public atm_dyn_dgm_trcadvect3d_heve_calc_fct_coef(fct_coef, qtrc_, momx_, momy_, momz_, rhoq_tp_, alphdens_m, alphdens_p, dens_hyd, ddens_, ddens0_, rk_c_ssm1, dt, faceintmat, lmesh, elem, lmesh2d, elem2d, disable_limiter)
subroutine, public atm_dyn_dgm_trcadvect3d_heve_cal_tend(qtrc_dt, qtrc_, momx_, momy_, momz_, alphdens_m, alphdens_p, fct_coef, rhoq_tp, element3d_operation, faceintmat, lmesh, elem, lmesh2d, elem2d)
subroutine, public atm_dyn_dgm_trcadvect3d_heve_cal_alphdens_advtest(alph_dens_m, alph_dens_p, ddens_, momx_, momy_, momz_, dens_hyd, gsqrt, nx, ny, nz, vmapm, vmapp, lmesh, elem)
subroutine, public atm_dyn_dgm_trcadvect3d_heve_init(mesh, faceintmat)
module FElib / Element / Base
module FElib / Element / hexahedron
module FElib / Element / Operation / Base
module FElib / Mesh / Local 2D
module FElib / Mesh / Local 3D
module FElib / Data / base
module FElib / Mesh / Base 2D
module FElib / Mesh / Base 3D
module FElib / Data / base
Module common / Polynomial.
real(rp) function, dimension(nord+1), public polynomial_gengausslobattoptintweight(nord)
A function to calculate the Gauss-Lobbato weights.
Module common / sparsemat.
Derived type representing a 2D reference element.
Derived type representing a 3D reference element.
Derived type representing a hexahedral element.
Base type for elementwise operations.
Derived type representing a local mesh for 2D domain.
Derived type to manage a local 3D computational domain.
Derived type representing a field with 3D local mesh.
Derived type to manage a computational mesh (base type for 2D domain)
Derived type to manage a computational mesh (base type for 3D domain)
Derived type representing a field with 3D mesh.
Derived type to manage a sparse matrix.