scale_atm_dyn_dgm_trcadvect3d_heve Module Reference

module FElib / Fluid dyn solver / Atmosphere / Tracer advection More...

Functions/Subroutines
subroutine, public	atm_dyn_dgm_trcadvect3d_heve_init (mesh, faceintmat)
subroutine, public	atm_dyn_dgm_trcadvect3d_heve_final ()
subroutine, public	atm_dyn_dgm_trcadvect3d_heve_cal_tend (qtrc_dt, qtrc_, momx_, momy_, momz_, alphdens_m, alphdens_p, fct_coef, rhoq_tp, dx, dy, dz, sx, sy, sz, lift, faceintmat, lmesh, elem, lmesh2d, elem2d)
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, dx, dy, dz, sx, sy, sz, lift, faceintmat, lmesh, elem, lmesh2d, elem2d, disable_limiter)
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)
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)

Detailed Description

module FElib / Fluid dyn solver / Atmosphere / Tracer advection

Description

HEVE DGM scheme for tracer advection.

To preserve nonnegativity, a limiter proposed by Light and Durran (2016, MWR) is used: we apply FCT for the lowest mode and the truncation and mass aware rescaling (TMAR).

Reference

• Light and Durran 2016: Preserving Nonnegativity in Discontinuous Galerkin Approximations to Scalar Transport via Truncation and Mass Aware Rescaling (TMAR). Monthly Weather Review, 144(12), 4771–4786.

Author

Yuta Kawai, Team SCALE

Function/Subroutine Documentation

atm_dyn_dgm_trcadvect3d_heve_init()

subroutine, public scale_atm_dyn_dgm_trcadvect3d_heve::atm_dyn_dgm_trcadvect3d_heve_init	(	class(meshbase3d), intent(in), target	mesh,
		type(sparsemat), intent(inout)	faceintmat )

Definition at line 74 of file scale_atm_dyn_dgm_trcadvect3d_heve.F90.

    use scale_polynominal, only: polynominal_gengausslobattoptintweight
    implicit none
    class(MeshBase3D), intent(in), target :: mesh
    type(SparseMat), intent(inout) :: FaceIntMat
 
    class(LocalMesh3D), pointer :: lcmesh
    class(ElementBase3D), pointer :: elem
    real(RP), allocatable :: intWeight_lgl1DPts_h(:)
    real(RP), allocatable :: intWeight_lgl1DPts_v(:)   
    real(RP), allocatable :: intWeight_h(:) 
    real(RP), allocatable :: intWeight_v(:)  
    
    integer :: f
    integer :: i, j, k, l
    integer :: is, ie
 
    real(RP), allocatable :: IntWeight(:,:)
    !-------------------------------------------------------------
 
    lcmesh => mesh%lcmesh_list(1)
    elem => lcmesh%refElem3D
    allocate( intweight(elem%Nfaces,elem%NfpTot) )
    intweight(:,:) = 0.0_rp
 
    allocate( intweight_lgl1dpts_h(elem%Nnode_h1D) )
    allocate( intweight_lgl1dpts_v(elem%Nnode_v) )
    allocate( intweight_h(elem%Nnode_h1D*elem%Nnode_v) )
    allocate( intweight_v(elem%Nnode_h1D**2) )
 
    intweight_lgl1dpts_h(:) = polynominal_gengausslobattoptintweight(elem%PolyOrder_h)
    intweight_lgl1dpts_v(:) = polynominal_gengausslobattoptintweight(elem%PolyOrder_v)
 
    do f=1, elem%Nfaces_h
      do k=1, elem%Nnode_v
      do i=1, elem%Nnode_h1D
        l = i + (k-1)*elem%Nnode_h1D
        intweight_h(l) = intweight_lgl1dpts_h(i) * intweight_lgl1dpts_v(k)
      end do
      end do
 
      is = (f-1)*elem%Nfp_h + 1
      ie = is + elem%Nfp_h - 1
      intweight(f,is:ie) = intweight_h(:)
    end do
 
    do f=1, elem%Nfaces_v
      do j=1, elem%Nnode_h1D
      do i=1, elem%Nnode_h1D
        l = i + (j-1)*elem%Nnode_h1D
        intweight_v(l) = intweight_lgl1dpts_h(i) * intweight_lgl1dpts_h(j)
      end do
      end do
 
      is = elem%Nfaces_h*elem%Nfp_h + (f-1)*elem%Nfp_v + 1
      ie = is + elem%Nfp_v - 1
      intweight(elem%Nfaces_h+f,is:ie) = intweight_v(:)
    end do
 
    call faceintmat%Init( intweight )
 
    return

References scale_polynominal::polynominal_gengausslobattoptintweight().

atm_dyn_dgm_trcadvect3d_heve_final()

subroutine, public scale_atm_dyn_dgm_trcadvect3d_heve::atm_dyn_dgm_trcadvect3d_heve_final

Definition at line 139 of file scale_atm_dyn_dgm_trcadvect3d_heve.F90.

    implicit none
    !--------------------------------------------
 
    return    

atm_dyn_dgm_trcadvect3d_heve_cal_tend()

subroutine, public scale_atm_dyn_dgm_trcadvect3d_heve::atm_dyn_dgm_trcadvect3d_heve_cal_tend	(	real(rp), dimension(elem%np,lmesh%nea), intent(out)	qtrc_dt,
		real(rp), dimension(elem%np,lmesh%nea), intent(in)	qtrc_,
		real(rp), dimension(elem%np,lmesh%nea), intent(in)	momx_,
		real(rp), dimension(elem%np,lmesh%nea), intent(in)	momy_,
		real(rp), dimension(elem%np,lmesh%nea), intent(in)	momz_,
		real(rp), dimension(elem%nfptot,lmesh%ne), intent(in)	alphdens_m,
		real(rp), dimension(elem%nfptot,lmesh%ne), intent(in)	alphdens_p,
		real(rp), dimension(elem%np,lmesh%nea), intent(in)	fct_coef,
		real(rp), dimension(elem%np,lmesh%nea), intent(in)	rhoq_tp,
		type(sparsemat), intent(in)	dx,
		type(sparsemat), intent(in)	dy,
		type(sparsemat), intent(in)	dz,
		type(sparsemat), intent(in)	sx,
		type(sparsemat), intent(in)	sy,
		type(sparsemat), intent(in)	sz,
		type(sparsemat), intent(in)	lift,
		type(sparsemat), intent(in)	faceintmat,
		class(localmesh3d), intent(in)	lmesh,
		class(elementbase3d), intent(in)	elem,
		class(localmesh2d), intent(in)	lmesh2d,
		class(elementbase2d), intent(in)	elem2d )

Definition at line 148 of file scale_atm_dyn_dgm_trcadvect3d_heve.F90.

 
    class(LocalMesh3D), intent(in) :: lmesh
    class(ElementBase3D), intent(in) :: elem
    class(LocalMesh2D), intent(in) :: lmesh2D
    class(ElementBase2D), intent(in) :: elem2D
    type(SparseMat), intent(in) :: Dx, Dy, Dz, Sx, Sy, Sz, Lift, FaceIntMat
    
    real(RP), intent(out) :: QTRC_dt(elem%Np,lmesh%NeA)
    real(RP), intent(in) :: QTRC_(elem%Np,lmesh%NeA)
    real(RP), intent(in) :: MOMX_(elem%Np,lmesh%NeA)
    real(RP), intent(in) :: MOMY_(elem%Np,lmesh%NeA)    
    real(RP), intent(in) :: MOMZ_(elem%Np,lmesh%NeA)    
    real(RP), intent(in) :: alphDENS_M(elem%NfpTot,lmesh%Ne)
    real(RP), intent(in) :: alphDENS_P(elem%NfpTot,lmesh%Ne)
    real(RP), intent(in) :: fct_coef(elem%Np,lmesh%NeA)
    real(RP), intent(in) :: RHOQ_tp(elem%Np,lmesh%NeA)
 
    real(RP) :: Fx(elem%Np), Fy(elem%Np), Fz(elem%Np), LiftDelFlx(elem%Np)    
    real(RP) :: del_flux(elem%NfpTot,lmesh%Ne)
    real(RP) :: momwt_(elem%Np)    
    
    integer :: ke, ke2d 
 
    real(RP) :: Q0, Q1, vol
    integer :: p
    !---------------------------------------------------------------------------
    
    call prof_rapstart('cal_trcadv_tend_bndflux', 3)
    call atm_dyn_dgm_trcadvect3d_heve_get_delflux_generalhvc( &
      del_flux,                                                                & ! (out)
      qtrc_, momx_, momy_, momz_, alphdens_m, alphdens_p, fct_coef,            & ! (in)
      lmesh%Gsqrt, lmesh%GsqrtH, lmesh%GI3(:,:,1), lmesh%GI3(:,:,2),           & ! (in)    
      lmesh%normal_fn(:,:,1), lmesh%normal_fn(:,:,2), lmesh%normal_fn(:,:,3),  & ! (in)
      lmesh%J(:,:), lmesh%Fscale(:,:),                                         & ! (in) 
      lmesh%vmapM, lmesh%vmapP, elem%IndexH2Dto3D_bnd, faceintmat,             & ! (in)
      lmesh, elem, lmesh2d, elem2d                                             ) ! (in)
    call prof_rapend('cal_trcadv_tend_bndflux', 3)
 
    call prof_rapstart('cal_trcadv_tend_interior', 3)
    !$omp parallel do private( &
    !$omp momwt_, Fx, Fy, Fz, LiftDelFlx,  &    
    !$omp ke, ke2d                         )  
    do ke=lmesh%NeS, lmesh%NeE
      ke2d = lmesh%EMap3Dto2D(ke)
      
      momwt_(:) = momz_(:,ke) * lmesh%GsqrtH(elem%IndexH2Dto3D,ke2d) / lmesh%Gsqrt(:,ke) &
                + lmesh%GI3(:,ke,1) * momx_(:,ke) + lmesh%GI3(:,ke,2) * momy_(:,ke)
 
      call sparsemat_matmul(dx, lmesh%Gsqrt(:,ke) * momx_(:,ke) * qtrc_(:,ke), fx)
      call sparsemat_matmul(dy, lmesh%Gsqrt(:,ke) * momy_(:,ke) * qtrc_(:,ke), fy)
      call sparsemat_matmul(dz, lmesh%Gsqrt(:,ke) * momwt_(:)   * qtrc_(:,ke), fz)
      call sparsemat_matmul(lift, lmesh%Fscale(:,ke) * del_flux(:,ke), liftdelflx)
 
      qtrc_dt(:,ke) = ( &
          - lmesh%Escale(:,ke,1,1) * fx(:)     &
          - lmesh%Escale(:,ke,2,2) * fy(:)     &
          - lmesh%Escale(:,ke,3,3) * fz(:)     &
          - liftdelflx(:)                   )  &
          / lmesh%Gsqrt(:,ke)                  &
          + rhoq_tp(:,ke)
    end do
    call prof_rapend('cal_trcadv_tend_interior', 3)
 
    return

atm_dyn_dgm_trcadvect3d_heve_calc_fct_coef()

subroutine, public scale_atm_dyn_dgm_trcadvect3d_heve::atm_dyn_dgm_trcadvect3d_heve_calc_fct_coef	(	real(rp), dimension(elem%np,lmesh%nea), intent(out)	fct_coef,
		real(rp), dimension(elem%np,lmesh%nea), intent(in)	qtrc_,
		real(rp), dimension(elem%np,lmesh%nea), intent(in)	momx_,
		real(rp), dimension(elem%np,lmesh%nea), intent(in)	momy_,
		real(rp), dimension(elem%np,lmesh%nea), intent(in)	momz_,
		real(rp), dimension(elem%np,lmesh%nea), intent(in)	rhoq_tp_,
		real(rp), dimension(elem%nfptot,lmesh%ne), intent(in)	alphdens_m,
		real(rp), dimension(elem%nfptot,lmesh%ne), intent(in)	alphdens_p,
		real(rp), dimension(elem%np,lmesh%nea), intent(in)	dens_hyd,
		real(rp), dimension (elem%np,lmesh%nea), intent(in)	ddens_,
		real(rp), dimension(elem%np,lmesh%nea), intent(in)	ddens0_,
		real(rp), intent(in)	rk_c_ssm1,
		real(rp), intent(in)	dt,
		type(sparsemat), intent(in)	dx,
		type(sparsemat), intent(in)	dy,
		type(sparsemat), intent(in)	dz,
		type(sparsemat), intent(in)	sx,
		type(sparsemat), intent(in)	sy,
		type(sparsemat), intent(in)	sz,
		type(sparsemat), intent(in)	lift,
		type(sparsemat), intent(in)	faceintmat,
		class(localmesh3d), intent(in)	lmesh,
		class(elementbase3d), intent(in)	elem,
		class(localmesh2d), intent(in)	lmesh2d,
		class(elementbase2d), intent(in)	elem2d,
		logical, intent(in), optional	disable_limiter )

Definition at line 222 of file scale_atm_dyn_dgm_trcadvect3d_heve.F90.

 
    implicit none
    class(LocalMesh3D), intent(in) :: lmesh
    class(ElementBase3D), intent(in) :: elem
    class(LocalMesh2D), intent(in) :: lmesh2D
    class(ElementBase2D), intent(in) :: elem2D    
    real(RP), intent(out) :: fct_coef(elem%Np,lmesh%NeA)
    real(RP), intent(in) :: QTRC_(elem%Np,lmesh%NeA)
    real(RP), intent(in) :: MOMX_(elem%Np,lmesh%NeA)
    real(RP), intent(in) :: MOMY_(elem%Np,lmesh%NeA)    
    real(RP), intent(in) :: MOMZ_(elem%Np,lmesh%NeA)
    real(RP), intent(in) :: RHOQ_tp_(elem%Np,lmesh%NeA)
    real(RP), intent(in) :: alphDENS_M(elem%NfpTot,lmesh%Ne)
    real(RP), intent(in) :: alphDENS_P(elem%NfpTot,lmesh%Ne)
    real(RP), intent(in) :: DENS_hyd(elem%Np,lmesh%NeA)     
    real(RP), intent(in) :: DDENS_ (elem%Np,lmesh%NeA)
    real(RP), intent(in) :: DDENS0_(elem%Np,lmesh%NeA)
    real(RP), intent(in) :: rk_c_ssm1
    real(RP), intent(in) :: dt
    type(SparseMat), intent(in) :: Dx, Dy, Dz, Sx, Sy, Sz, Lift, FaceIntMat
    logical, intent(in), optional :: disable_limiter
 
    real(RP) :: netOutwardFlux(lmesh%Ne)
    real(RP) :: momwt_(elem%Np)    
    
    integer :: ke
    real(RP) :: Q
    real(RP) :: dens_ssm1(elem%Np)
    !---------------------------------------------------------------------------
 
    if ( present(disable_limiter) ) then
      if ( disable_limiter ) then
        !$omp parallel do
        do ke=lmesh%NeS, lmesh%NeE
          fct_coef(:,ke) = 1.0_rp
        end do
        return
      end if 
    end if
 
    call prof_rapstart('cal_trcadv_fct_coef_bndflux', 3)
    call atm_dyn_dgm_trcadvect3d_heve_get_netoutwardflux_generalhvc( &
      netoutwardflux,                                                                   & ! (out)
      qtrc_, momx_, momy_, momz_, alphdens_m, alphdens_p,                               & ! (in)
      lmesh%Gsqrt, lmesh%GsqrtH, lmesh%GI3(:,:,1), lmesh%GI3(:,:,2),                    & ! (in)    
      lmesh%normal_fn(:,:,1), lmesh%normal_fn(:,:,2), lmesh%normal_fn(:,:,3),           & ! (in)
      lmesh%J(:,:), lmesh%Fscale(:,:), lmesh%vmapM, lmesh%vmapP, elem%IndexH2Dto3D_bnd, & ! (in)
      faceintmat,                                                                       & ! (in)
      lmesh, elem, lmesh2d, elem2d                                                      ) ! (in)
    call prof_rapend('cal_trcadv_fct_coef_bndflux', 3)
 
    call prof_rapstart('cal_trcadv_fct_coef', 3)
 
    !$omp parallel do private( ke, Q, dens_ssm1 )  
    do ke=lmesh%NeS, lmesh%NeE
 
      dens_ssm1(:) = dens_hyd(:,ke) &
                    + ( 1.0_rp - rk_c_ssm1 ) * ddens0_(:,ke) + rk_c_ssm1 * ddens_(:,ke)
      q = sum( lmesh%Gsqrt(:,ke) * lmesh%J(:,ke) * elem%IntWeight_lgl(:) * ( dens_ssm1(:) *  qtrc_(:,ke) / dt + rhoq_tp_(:,ke) ) )      
      
      fct_coef(:,ke) = max( 0.0_rp, min( 1.0_rp, q / ( netoutwardflux(ke) + 1.0e-10_rp ) ) )
    end do
 
    call prof_rapend('cal_trcadv_fct_coef', 3)
 
    return

atm_dyn_dgm_trcadvect3d_tmar()

subroutine, public scale_atm_dyn_dgm_trcadvect3d_heve::atm_dyn_dgm_trcadvect3d_tmar	(	real(rp), dimension(elem%np,lmesh%nea), intent(inout)	qtrc_,
		real(rp), dimension(elem%np,lmesh%nea), intent(in)	dens_hyd,
		real(rp), dimension (elem%np,lmesh%nea), intent(in)	ddens_,
		class(localmesh3d), intent(in)	lmesh,
		class(elementbase3d), intent(in)	elem,
		class(localmesh2d), intent(in)	lmesh2d,
		class(elementbase2d), intent(in)	elem2d )

Second Step of limiter in which nonlinear truncation and mass aware rescaling (TMAR)

Definition at line 300 of file scale_atm_dyn_dgm_trcadvect3d_heve.F90.

    
    implicit none
    class(LocalMesh3D), intent(in) :: lmesh
    class(ElementBase3D), intent(in) :: elem
    class(LocalMesh2D), intent(in) :: lmesh2D
    class(ElementBase2D), intent(in) :: elem2D
    real(RP), intent(inout) :: QTRC_(elem%Np,lmesh%NeA)
    real(RP), intent(in) :: DENS_hyd(elem%Np,lmesh%NeA)     
    real(RP), intent(in) :: DDENS_  (elem%Np,lmesh%NeA)
 
    integer :: ke
    real(RP) :: Q0, Q1
    real(RP) :: Q(elem%Np)
    real(RP) :: dens(elem%Np)
    !--------------------------------------------------------
 
    !$omp parallel do private( dens, Q, Q0, Q1 )
    do ke=lmesh%NeS, lmesh%NeE
      dens(:) = dens_hyd(:,ke) + ddens_(:,ke)
      q(:) = max( 0.0_rp, qtrc_(:,ke) )
 
      q0 = sum( lmesh%Gsqrt(:,ke) * lmesh%J(:,ke) * elem%IntWeight_lgl(:) * dens(:) * qtrc_(:,ke) )
      q1 = sum( lmesh%Gsqrt(:,ke) * lmesh%J(:,ke) * elem%IntWeight_lgl(:) * dens(:) * q(:)        )
      qtrc_(:,ke) = q0 / ( q1 + 1.0e-32_rp ) * q(:)
    end do
 
    return

atm_dyn_dgm_trcadvect3d_save_massflux()

subroutine, public scale_atm_dyn_dgm_trcadvect3d_heve::atm_dyn_dgm_trcadvect3d_save_massflux	(	real(rp), dimension(elem%np,lmesh%nea), intent(inout)	mflx_x_tavg,
		real(rp), dimension(elem%np,lmesh%nea), intent(inout)	mflx_y_tavg,
		real(rp), dimension(elem%np,lmesh%nea), intent(inout)	mflx_z_tavg,
		real(rp), dimension(elem%nfptot,lmesh%ne), intent(inout)	alph_dens_m,
		real(rp), dimension(elem%nfptot,lmesh%ne), intent(inout)	alph_dens_p,
		real(rp), dimension(elem%np,lmesh%nea), intent(in)	ddens,
		real(rp), dimension(elem%np,lmesh%nea), intent(in)	momx,
		real(rp), dimension(elem%np,lmesh%nea), intent(in)	momy,
		real(rp), dimension(elem%np,lmesh%nea), intent(in)	momz,
		real(rp), dimension(elem%np,lmesh%nea), intent(in)	dpres,
		real(rp), dimension(elem%np,lmesh%nea), intent(in)	dens_hyd,
		real(rp), dimension(elem%np,lmesh%nea), intent(in)	pres_hyd,
		real(rp), dimension(elem%np,lmesh%nea), intent(in)	rtot,
		real(rp), dimension(elem%np,lmesh%nea), intent(in)	cvtot,
		real(rp), dimension(elem%np,lmesh%nea), intent(in)	cptot,
		class(localmesh3d), intent(in)	lmesh,
		class(elementbase3d), intent(in)	elem,
		integer, intent(in)	rkstage,
		real(rp), intent(in)	tavg_weight_h,
		real(rp), intent(in)	tavg_weight_v )

Definition at line 332 of file scale_atm_dyn_dgm_trcadvect3d_heve.F90.

   
    implicit none
    class(LocalMesh3D), intent(in) :: lmesh
    class(ElementBase3D), intent(in) :: elem
    real(RP), intent(inout) ::  MFLX_x_tavg(elem%Np,lmesh%NeA)
    real(RP), intent(inout) ::  MFLX_y_tavg(elem%Np,lmesh%NeA)
    real(RP), intent(inout) ::  MFLX_z_tavg(elem%Np,lmesh%NeA)
    real(RP), intent(inout) ::  alph_dens_M(elem%NfpTot,lmesh%Ne)
    real(RP), intent(inout) ::  alph_dens_P(elem%NfpTot,lmesh%Ne)
    real(RP), intent(in) ::  DDENS(elem%Np,lmesh%NeA)
    real(RP), intent(in) ::  MOMX(elem%Np,lmesh%NeA)  
    real(RP), intent(in) ::  MOMY(elem%Np,lmesh%NeA)  
    real(RP), intent(in) ::  MOMZ(elem%Np,lmesh%NeA)  
    real(RP), intent(in) ::  DPRES(elem%Np,lmesh%NeA)
    real(RP), intent(in) ::  DENS_hyd(elem%Np,lmesh%NeA)
    real(RP), intent(in) ::  PRES_hyd(elem%Np,lmesh%NeA)
    real(RP), intent(in) ::  Rtot(elem%Np,lmesh%NeA)
    real(RP), intent(in) ::  CVtot(elem%Np,lmesh%NeA)
    real(RP), intent(in) ::  CPtot(elem%Np,lmesh%NeA)
    integer, intent(in) :: rkstage
    real(RP), intent(in) :: tavg_weight_h
    real(RP), intent(in) :: tavg_weight_v
 
    integer :: ke
    !--------------------------------------------------------------
 
    !$omp parallel do private(ke)
    do ke=lmesh%NeS, lmesh%NeE
      if (rkstage == 1) then
        mflx_x_tavg(:,ke) = tavg_weight_h * momx(:,ke)
        mflx_y_tavg(:,ke) = tavg_weight_h * momy(:,ke)
        mflx_z_tavg(:,ke) = tavg_weight_v * momz(:,ke)
        alph_dens_m(:,ke) = 0.0_rp
        alph_dens_p(:,ke) = 0.0_rp
      else
        mflx_x_tavg(:,ke) = mflx_x_tavg(:,ke) + tavg_weight_h * momx(:,ke)
        mflx_y_tavg(:,ke) = mflx_y_tavg(:,ke) + tavg_weight_h * momy(:,ke)
        mflx_z_tavg(:,ke) = mflx_z_tavg(:,ke) + tavg_weight_v * momz(:,ke)
      end if  
    end do
 
    call atm_dyn_dgm_trcadvect3d_heve_cal_alphdens_dyn( &
      alph_dens_m, alph_dens_p,                                               & ! (out)
      ddens, momx, momy, momz, dpres, dens_hyd, pres_hyd,                     & ! (in)
      rtot, cvtot, cptot,                                                     & ! (in)
      lmesh%Gsqrt, lmesh%GIJ(:,:,1,1), lmesh%GIJ(:,:,1,2), lmesh%GIJ(:,:,2,2), & ! (in)
      lmesh%normal_fn(:,:,1), lmesh%normal_fn(:,:,2), lmesh%normal_fn(:,:,3), & ! (in)
      lmesh%vmapM, lmesh%vmapP, elem%IndexH2Dto3D_bnd,                        & ! (in)
      lmesh, lmesh%lcmesh2D, lmesh%refElem3D, lmesh%lcmesh2D%refElem2D,       & ! (in)
      tavg_weight_h, tavg_weight_v                                            ) ! (in)
    
    return

atm_dyn_dgm_trcadvect3d_heve_cal_alphdens_advtest()

subroutine, public scale_atm_dyn_dgm_trcadvect3d_heve::atm_dyn_dgm_trcadvect3d_heve_cal_alphdens_advtest	(	real(rp), dimension(elem%nfptot*lmesh%ne), intent(inout)	alph_dens_m,
		real(rp), dimension(elem%nfptot*lmesh%ne), intent(inout)	alph_dens_p,
		real(rp), dimension(elem%np*lmesh%nea), intent(in)	ddens_,
		real(rp), dimension(elem%np*lmesh%nea), intent(in)	momx_,
		real(rp), dimension(elem%np*lmesh%nea), intent(in)	momy_,
		real(rp), dimension(elem%np*lmesh%nea), intent(in)	momz_,
		real(rp), dimension(elem%np*lmesh%nea), intent(in)	dens_hyd,
		real(rp), dimension(elem%np*lmesh%nea), intent(in)	gsqrt,
		real(rp), dimension(elem%nfptot*lmesh%ne), intent(in)	nx,
		real(rp), dimension(elem%nfptot*lmesh%ne), intent(in)	ny,
		real(rp), dimension(elem%nfptot*lmesh%ne), intent(in)	nz,
		integer, dimension(elem%nfptot*lmesh%ne), intent(in)	vmapm,
		integer, dimension(elem%nfptot*lmesh%ne), intent(in)	vmapp,
		class(localmesh3d), intent(in)	lmesh,
		class(elementbase3d), intent(in)	elem )

Definition at line 392 of file scale_atm_dyn_dgm_trcadvect3d_heve.F90.

 
    use scale_const, only: &
     grav => const_grav,  &
     rdry => const_rdry,  &
     cpdry => const_cpdry, &
     cvdry => const_cvdry, &
     pres00 => const_pre00
   
    implicit none
 
    class(LocalMesh3D), intent(in) :: lmesh
    class(ElementBase3D), intent(in) :: elem  
    real(RP), intent(inout) ::  alph_dens_M(elem%NfpTot*lmesh%Ne)
    real(RP), intent(inout) ::  alph_dens_P(elem%NfpTot*lmesh%Ne)
    real(RP), intent(in) ::  DDENS_(elem%Np*lmesh%NeA)
    real(RP), intent(in) ::  MOMX_(elem%Np*lmesh%NeA)  
    real(RP), intent(in) ::  MOMY_(elem%Np*lmesh%NeA)  
    real(RP), intent(in) ::  MOMZ_(elem%Np*lmesh%NeA)  
    real(RP), intent(in) ::  DENS_hyd(elem%Np*lmesh%NeA)
    real(RP), intent(in) :: Gsqrt(elem%Np*lmesh%NeA)
    real(RP), intent(in) :: nx(elem%NfpTot*lmesh%Ne)
    real(RP), intent(in) :: ny(elem%NfpTot*lmesh%Ne)
    real(RP), intent(in) :: nz(elem%NfpTot*lmesh%Ne)
    integer, intent(in) :: vmapM(elem%NfpTot*lmesh%Ne)
    integer, intent(in) :: vmapP(elem%NfpTot*lmesh%Ne)
    
    integer :: i, iP, iM
    real(RP) :: VelP, VelM, alpha, densM, densP
    !------------------------------------------------------------------------
 
 
    !$omp parallel do private( &
    !$omp iM, iP, VelP, VelM, alpha, densM, densP )
    do i=1, elem%NfpTot*lmesh%Ne
      im = vmapm(i); ip = vmapp(i)
  
      densm = ddens_(im) + dens_hyd(im)
      densp = ddens_(ip) + dens_hyd(ip)
 
      velm = ( momx_(im) * nx(i) + momy_(im) * ny(i) + momz_(im) * nz(i) ) / densm
      velp = ( momx_(ip) * nx(i) + momy_(ip) * ny(i) + momz_(ip) * nz(i) ) / densp
 
      alpha = max( abs(velm), abs(velp)  )
      alph_dens_m(i) = alpha * densm * gsqrt(im)
      alph_dens_p(i) = alpha * densp * gsqrt(ip)
    end do
 
    return