scale_atm_dyn_dgm_globalnonhydro3d_rhot_hevi Module Reference

module FElib / Fluid dyn solver / Atmosphere / Global nonhydrostatic model / HEVI More...

Functions/Subroutines
subroutine, public	atm_dyn_dgm_globalnonhydro3d_rhot_hevi_init (mesh)
subroutine, public	atm_dyn_dgm_globalnonhydro3d_rhot_hevi_final ()
subroutine, public	atm_dyn_dgm_globalnonhydro3d_rhot_hevi_cal_tend (dens_dt, momx_dt, momy_dt, momz_dt, rhot_dt, ddens_, momx_, momy_, momz_, drhot_, dpres_, dens_hyd, pres_hyd, pres_hyd_ref, coriolis, rtot, cvtot, cptot, dphyddx, dphyddy, element3d_operation, dx, dy, dz, sx, sy, sz, lift, lmesh, elem, lmesh2d, elem2d)
subroutine, public	atm_dyn_dgm_globalnonhydro3d_rhot_hevi_cal_vi (dens_dt, momx_dt, momy_dt, momz_dt, rhot_dt, ddens_, momx_, momy_, momz_, drhot_, dens_hyd, pres_hyd, ddens0_, momx0_, momy0_, momz0_, drhot0_, rtot, cvtot, cptot, element3d_operation, dz, lift, impl_fac, dt, lmesh, elem, lmesh2d, elem2d)

Detailed Description

module FElib / Fluid dyn solver / Atmosphere / Global nonhydrostatic model / HEVI

Description

HEVI DGM scheme for Global Atmospheric Dynamical process. The governing equations is a fully compressible nonhydrostatic equations, which consist of mass, momentum, and thermodynamics (density * potential temperature conservation) equations.

Author

Yuta Kawai, Team SCALE

Function/Subroutine Documentation

atm_dyn_dgm_globalnonhydro3d_rhot_hevi_init()

subroutine, public scale_atm_dyn_dgm_globalnonhydro3d_rhot_hevi::atm_dyn_dgm_globalnonhydro3d_rhot_hevi_init

(

class(meshbase3d), intent(in)

mesh

)

Definition at line 75 of file scale_atm_dyn_dgm_globalnonhydro3d_rhot_hevi.F90.

    implicit none
    class(MeshBase3D), intent(in) :: mesh
    !--------------------------------------------
 
    call atm_dyn_dgm_nonhydro3d_common_init( mesh )
    return

References scale_atm_dyn_dgm_nonhydro3d_common::atm_dyn_dgm_nonhydro3d_common_init().

atm_dyn_dgm_globalnonhydro3d_rhot_hevi_final()

subroutine, public scale_atm_dyn_dgm_globalnonhydro3d_rhot_hevi::atm_dyn_dgm_globalnonhydro3d_rhot_hevi_final

Definition at line 85 of file scale_atm_dyn_dgm_globalnonhydro3d_rhot_hevi.F90.

    implicit none
    !--------------------------------------------
    
    call atm_dyn_dgm_nonhydro3d_common_final()    
    return

References scale_atm_dyn_dgm_nonhydro3d_common::atm_dyn_dgm_nonhydro3d_common_final().

atm_dyn_dgm_globalnonhydro3d_rhot_hevi_cal_tend()

subroutine, public scale_atm_dyn_dgm_globalnonhydro3d_rhot_hevi::atm_dyn_dgm_globalnonhydro3d_rhot_hevi_cal_tend	(	real(rp), dimension(elem%np,lmesh%nea), intent(out)	dens_dt,
		real(rp), dimension(elem%np,lmesh%nea), intent(out)	momx_dt,
		real(rp), dimension(elem%np,lmesh%nea), intent(out)	momy_dt,
		real(rp), dimension(elem%np,lmesh%nea), intent(out)	momz_dt,
		real(rp), dimension(elem%np,lmesh%nea), intent(out)	rhot_dt,
		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)	drhot_,
		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)	pres_hyd_ref,
		real(rp), dimension(elem2d%np,lmesh2d%nea), intent(in)	coriolis,
		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,
		real(rp), dimension(elem%np,lmesh%nea), intent(in)	dphyddx,
		real(rp), dimension(elem%np,lmesh%nea), intent(in)	dphyddy,
		class(elementoperationbase3d), intent(in)	element3d_operation,
		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,
		class(localmesh3d), intent(in)	lmesh,
		class(elementbase3d), intent(in)	elem,
		class(localmesh2d), intent(in)	lmesh2d,
		class(elementbase2d), intent(in)	elem2d )

Definition at line 96 of file scale_atm_dyn_dgm_globalnonhydro3d_rhot_hevi.F90.

 
    use scale_atm_dyn_dgm_nonhydro3d_rhot_hevi_numflux, only: &
      get_ebnd_flux => atm_dyn_dgm_nonhydro3d_rhot_hevi_numflux_get_generalhvc
    use scale_const, only: &
      ohm => const_ohm
    implicit none
 
    class(LocalMesh3D), intent(in) :: lmesh
    class(ElementBase3D), intent(in) :: elem
    class(LocalMesh2D), intent(in) :: lmesh2D
    class(ElementBase2D), intent(in) :: elem2D
    class(ElementOperationBase3D), intent(in) :: element3D_operation
    type(SparseMat), intent(in) :: Dx, Dy, Dz, Sx, Sy, Sz, Lift
    real(RP), intent(out) :: DENS_dt(elem%Np,lmesh%NeA)
    real(RP), intent(out) :: MOMX_dt(elem%Np,lmesh%NeA)
    real(RP), intent(out) :: MOMY_dt(elem%Np,lmesh%NeA)
    real(RP), intent(out) :: MOMZ_dt(elem%Np,lmesh%NeA)
    real(RP), intent(out) :: RHOT_dt(elem%Np,lmesh%NeA)
    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)  :: DRHOT_(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)  :: PRES_hyd_ref(elem%Np,lmesh%NeA)
    real(RP), intent(in)  :: CORIOLIS(elem2D%Np,lmesh2D%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)
    real(RP), intent(in) :: DPhydDx(elem%Np,lmesh%NeA)
    real(RP), intent(in) :: DPhydDy(elem%Np,lmesh%NeA)
    
    real(RP) :: Fx(elem%Np), Fy(elem%Np), Fz(elem%Np), LiftDelFlx(elem%Np)
    real(RP) :: DPRES_hyd(elem%Np), GradPhyd_x(elem%Np), GradPhyd_y(elem%Np)
    real(RP) :: del_flux(elem%NfpTot,lmesh%Ne,PRGVAR_NUM)
    real(RP) :: del_flux_hyd(elem%NfpTot,lmesh%Ne,2)
    real(RP) :: RHOT_(elem%Np)
    real(RP) :: rdens_(elem%Np), u_(elem%Np), v_(elem%Np), w_(elem%Np), wt_(elem%Np)
 
    real(RP) :: G11(elem%Np), G12(elem%Np), G22(elem%Np), Rgam2(elem%Np)
    real(RP) :: GsqrtV(elem%Np), RGsqrtV(elem%Np)
    real(RP) :: X2D(elem2D%Np,lmesh2D%Ne), Y2D(elem2D%Np,lmesh2D%Ne)
    real(RP) :: X(elem%Np), Y(elem%Np), twoOVdel2(elem%Np)
    real(RP) :: CORI(elem%Np,2)
    logical :: is_panel1to4
    real(RP) :: s
 
    integer :: ke, ke2d
 
    real(RP) :: gamm, rgamm    
    real(RP) :: rP0
    real(RP) :: RovP0, P0ovR
    !------------------------------------------------------------------------
 
    call prof_rapstart('cal_dyn_tend_bndflux', 3)
    call get_ebnd_flux( &
      del_flux, del_flux_hyd,                                                  & ! (out)
      ddens_, momx_, momy_, momz_, drhot_, 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%GsqrtH, lmesh%GI3(:,:,1), lmesh%GI3(:,:,2),                        & ! (in)    
      lmesh%normal_fn(:,:,1), lmesh%normal_fn(:,:,2), lmesh%normal_fn(:,:,3),  & ! (in)
      lmesh%vmapM, lmesh%vmapP, elem%IndexH2Dto3D_bnd,                         & ! (in)
      lmesh, elem, lmesh2d, elem2d                                             ) ! (in)
    call prof_rapend('cal_dyn_tend_bndflux', 3)
    
    !-----
    call prof_rapstart('cal_dyn_tend_interior', 3)
    gamm  = cpdry / cvdry
    rgamm = cvdry / cpdry
    rp0   = 1.0_rp / pres00
    rovp0 = rdry * rp0
    p0ovr = pres00 / rdry
 
    s = 1.0_rp 
    is_panel1to4 = .true.
    if ( lmesh%panelID == 5 ) then
      is_panel1to4 = .false.
    else if ( lmesh%panelID == 6 ) then
      is_panel1to4 = .false.
      s = - 1.0_rp
    end if
 
    !$omp parallel private(                        &
    !$omp RHOT_, rdens_, u_, v_, w_, wt_,          &
    !$omp Fx, Fy, Fz, LiftDelFlx,                  &
    !$omp DPRES_hyd, GradPhyd_x, GradPhyd_y,       &
    !$omp G11, G12, G22, Rgam2, GsqrtV, RGsqrtV,   &
    !$omp X, Y, twoOVdel2,                         &
    !$omp CORI, ke, ke2D                           )
 
    !$omp do
    do ke2d = lmesh2d%NeS, lmesh2d%NeE
      x2d(:,ke2d) = tan(lmesh2d%pos_en(:,ke2d,1))
      y2d(:,ke2d) = tan(lmesh2d%pos_en(:,ke2d,2))
    end do
 
    !$omp do
    do ke = lmesh%NeS, lmesh%NeE
      !--
      ke2d = lmesh%EMap3Dto2D(ke)
      rgam2(:) = 1.0_rp / lmesh%gam(:,ke)**2
      g11(:) = lmesh%GIJ(elem%IndexH2Dto3D,ke2d,1,1) * rgam2(:)
      g12(:) = lmesh%GIJ(elem%IndexH2Dto3D,ke2d,1,2) * rgam2(:)
      g22(:) = lmesh%GIJ(elem%IndexH2Dto3D,ke2d,2,2) * rgam2(:)
      gsqrtv(:)  = lmesh%Gsqrt(:,ke) * rgam2(:) / lmesh%GsqrtH(elem%IndexH2Dto3D,ke2d)
      rgsqrtv(:) = 1.0_rp / gsqrtv(:)
 
      !--
      rhot_(:) = p0ovr * ( pres_hyd(:,ke) * rp0 )**rgamm + drhot_(:,ke)
      ! DPRES_(:) = PRES00 * ( Rtot(:,ke) * rP0 * RHOT_(:) )**( CPtot(:,ke) / CVtot(:,ke) ) &
      !           - PRES_hyd(:,ke)
      
      rdens_(:) = 1.0_rp / ( ddens_(:,ke) + dens_hyd(:,ke) )
      u_(:) = momx_(:,ke) * rdens_(:)
      v_(:) = momy_(:,ke) * rdens_(:)
      w_(:) = momz_(:,ke) * rdens_(:)
      wt_(:) = w_(:) * rgsqrtv(:) + lmesh%GI3(:,ke,1) * u_(:) + lmesh%GI3(:,ke,2) * v_(:) 
 
      x(:) = x2d(elem%IndexH2Dto3D,ke2d)
      y(:) = y2d(elem%IndexH2Dto3D,ke2d)
      twoovdel2(:) = 2.0_rp / ( 1.0_rp + x(:)**2 + y(:)**2 )
 
      cori(:,1) = s * ohm * twoovdel2(:) * ( - x(:) * y(:)        * momx_(:,ke) + (1.0_rp + y(:)**2) * momy_(:,ke) )
      cori(:,2) = s * ohm * twoovdel2(:) * ( - (1.0_rp + x(:)**2) * momx_(:,ke) +  x(:) * y(:)       * momy_(:,ke) )
      if ( is_panel1to4 ) then
        cori(:,1) = s * y(:) * cori(:,1)
        cori(:,2) = s * y(:) * cori(:,2)
      end if
      
      !-- Gradient hydrostatic pressure
 
      dpres_hyd(:) = pres_hyd(:,ke) - pres_hyd_ref(:,ke)
 
      call sparsemat_matmul(dx, gsqrtv(:) * dpres_hyd(:), fx)
      call sparsemat_matmul(dz, gsqrtv(:) * lmesh%GI3(:,ke,1) * dpres_hyd(:), fz)
      call sparsemat_matmul(lift, lmesh%Fscale(:,ke) * del_flux_hyd(:,ke,1), liftdelflx)
      gradphyd_x(:) = lmesh%Escale(:,ke,1,1) * fx(:) &
                    + lmesh%Escale(:,ke,3,3) * fz(:) &
                    + liftdelflx(:)
 
      call sparsemat_matmul(dy, gsqrtv(:) * dpres_hyd(:), fy)
      call sparsemat_matmul(dz, gsqrtv(:) * lmesh%GI3(:,ke,2) * dpres_hyd(:), fz)
      call sparsemat_matmul(lift, lmesh%Fscale(:,ke) * del_flux_hyd(:,ke,2), liftdelflx)
      gradphyd_y(:) = lmesh%Escale(:,ke,2,2) * fy(:) &
                    + lmesh%Escale(:,ke,3,3) * fz(:) &
                    + liftdelflx(:)
      
      !-- DENS
      call sparsemat_matmul(dx, lmesh%Gsqrt(:,ke) * momx_(:,ke), fx)
      call sparsemat_matmul(dy, lmesh%Gsqrt(:,ke) * momy_(:,ke), fy)
      call sparsemat_matmul(lift, lmesh%Fscale(:,ke) * del_flux(:,ke,dens_vid), liftdelflx)
 
      dens_dt(:,ke) = - ( &
            lmesh%Escale(:,ke,1,1) * fx(:)     &
          + lmesh%Escale(:,ke,2,2) * fy(:)     &
          + liftdelflx(:) ) / lmesh%Gsqrt(:,ke)
      
      !-- MOMX
      call sparsemat_matmul(dx, lmesh%Gsqrt(:,ke) * ( u_(:) * momx_(:,ke) + g11(:) * dpres_(:,ke) ), fx)
      call sparsemat_matmul(dy, lmesh%Gsqrt(:,ke) * ( v_(:) * momx_(:,ke) + g12(:) * dpres_(:,ke) ), fy)
      call sparsemat_matmul(dz, lmesh%Gsqrt(:,ke) * ( wt_(:) * momx_(:,ke)                              &
                                                    + ( lmesh%GI3(:,ke,1) * g11(:) + lmesh%GI3(:,ke,2) * g12(:) ) * dpres_(:,ke) ), fz)
      call sparsemat_matmul(lift, lmesh%Fscale(:,ke) * del_flux(:,ke,momx_vid), liftdelflx)
 
      momx_dt(:,ke) = &
          - ( lmesh%Escale(:,ke,1,1) * fx(:)                                     &
            + lmesh%Escale(:,ke,2,2) * fy(:)                                     &
            + lmesh%Escale(:,ke,3,3) * fz(:)                                     &
            + liftdelflx(:)                   ) / lmesh%Gsqrt(:,ke)              &
          + twoovdel2(:) * y(:) *                                                &
            ( - x(:) * y(:) * u_(:) + (1.0_rp + y(:)**2) * v_(:) ) * momx_(:,ke) &
          - ( g11(:) * gradphyd_x(:) + g12(:) * gradphyd_y(:) ) * rgsqrtv(:)     &
          + cori(:,1)              
      
      !-- MOMY
      call sparsemat_matmul(dx, lmesh%Gsqrt(:,ke) * ( u_(:) * momy_(:,ke) + g12(:) * dpres_(:,ke) ), fx)
      call sparsemat_matmul(dy, lmesh%Gsqrt(:,ke) * ( v_(:) * momy_(:,ke) + g22(:) * dpres_(:,ke) ), fy)
      call sparsemat_matmul(dz, lmesh%Gsqrt(:,ke) * ( wt_(:) * momy_(:,ke)                              &
                                                    + ( lmesh%GI3(:,ke,1) * g12(:) + lmesh%GI3(:,ke,2) * g22(:) ) * dpres_(:,ke) ), fz)
      call sparsemat_matmul(lift, lmesh%Fscale(:,ke) * del_flux(:,ke,momy_vid), liftdelflx)
 
      momy_dt(:,ke) = &
            - ( lmesh%Escale(:,ke,1,1) * fx(:)                                   &
              + lmesh%Escale(:,ke,2,2) * fy(:)                                   &
              + lmesh%Escale(:,ke,3,3) * fz(:)                                   &
              + liftdelflx(:)                  ) / lmesh%Gsqrt(:,ke)             &
            + twoovdel2(:) * x(:) *                                              &
              ( (1.0_rp + x(:)**2) * u_(:) - x(:) * y(:) * v_(:) ) * momy_(:,ke) &
            - ( g12(:) * gradphyd_x(:) + g22(:) * gradphyd_y(:) ) * rgsqrtv(:)   &
            + cori(:,2) 
      
      !-- MOMZ
      call sparsemat_matmul(dx, lmesh%Gsqrt(:,ke) * u_(:) * momz_(:,ke), fx)
      call sparsemat_matmul(dy, lmesh%Gsqrt(:,ke) * v_(:) * momz_(:,ke), fy)
      call sparsemat_matmul(dz, lmesh%Gsqrt(:,ke) * wt_(:) * momz_(:,ke), fz)
      call sparsemat_matmul(lift, lmesh%Fscale(:,ke) * del_flux(:,ke,momz_vid), liftdelflx)
      
      momz_dt(:,ke) = &
          - ( lmesh%Escale(:,ke,1,1) * fx(:)       &
            + lmesh%Escale(:,ke,2,2) * fy(:)       &
            + lmesh%Escale(:,ke,3,3) * fz(:)       &
            + liftdelflx(:) ) / lmesh%Gsqrt(:,ke)  
 
      !-- RHOT
      call sparsemat_matmul(dx, lmesh%Gsqrt(:,ke) * u_(:) * rhot_(:), fx)
      call sparsemat_matmul(dy, lmesh%Gsqrt(:,ke) * v_(:) * rhot_(:), fy)
      call sparsemat_matmul(lift, lmesh%Fscale(:,ke) * del_flux(:,ke,rhot_vid), liftdelflx)
      
      rhot_dt(:,ke) = &
          - ( lmesh%Escale(:,ke,1,1) * fx(:)      &
            + lmesh%Escale(:,ke,2,2) * fy(:)      &
            + liftdelflx(:) ) / lmesh%Gsqrt(:,ke)
    end do
    !$omp end parallel
    call prof_rapend('cal_dyn_tend_interior', 3)
 
    return

References scale_atm_dyn_dgm_nonhydro3d_rhot_hevi_numflux::atm_dyn_dgm_nonhydro3d_rhot_hevi_numflux_get_generalhvc().

atm_dyn_dgm_globalnonhydro3d_rhot_hevi_cal_vi()

subroutine, public scale_atm_dyn_dgm_globalnonhydro3d_rhot_hevi::atm_dyn_dgm_globalnonhydro3d_rhot_hevi_cal_vi	(	real(rp), dimension(elem%np,lmesh%nea), intent(out)	dens_dt,
		real(rp), dimension(elem%np,lmesh%nea), intent(out)	momx_dt,
		real(rp), dimension(elem%np,lmesh%nea), intent(out)	momy_dt,
		real(rp), dimension(elem%np,lmesh%nea), intent(out)	momz_dt,
		real(rp), dimension(elem%np,lmesh%nea), intent(out)	rhot_dt,
		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)	drhot_,
		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)	ddens0_,
		real(rp), dimension(elem%np,lmesh%nea), intent(in)	momx0_,
		real(rp), dimension(elem%np,lmesh%nea), intent(in)	momy0_,
		real(rp), dimension(elem%np,lmesh%nea), intent(in)	momz0_,
		real(rp), dimension(elem%np,lmesh%nea), intent(in)	drhot0_,
		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(elementoperationbase3d), intent(in)	element3d_operation,
		class(sparsemat), intent(in)	dz,
		class(sparsemat), intent(in)	lift,
		real(rp), intent(in)	impl_fac,
		real(rp), intent(in)	dt,
		class(localmesh3d), intent(in)	lmesh,
		class(elementbase3d), intent(in)	elem,
		class(localmesh2d), intent(in)	lmesh2d,
		class(elementbase2d), intent(in)	elem2d )

Definition at line 325 of file scale_atm_dyn_dgm_globalnonhydro3d_rhot_hevi.F90.

 
    use scale_atm_dyn_dgm_nonhydro3d_rhot_hevi_common, only: &
      vi_eval_ax => atm_dyn_dgm_nonhydro3d_rhot_hevi_common_eval_ax,                  &
      vi_eval_ax_uv => atm_dyn_dgm_nonhydro3d_rhot_hevi_common_eval_ax_uv,                  &
      vi_construct_matbnd => atm_dyn_dgm_nonhydro3d_rhot_hevi_common_construct_matbnd, &
      vi_construct_matbnd_uv => atm_dyn_dgm_nonhydro3d_rhot_hevi_common_construct_matbnd_uv
  
    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) :: DENS_dt(elem%Np,lmesh%NeA)
    real(RP), intent(out) :: MOMX_dt(elem%Np,lmesh%NeA)
    real(RP), intent(out) :: MOMY_dt(elem%Np,lmesh%NeA)
    real(RP), intent(out) :: MOMZ_dt(elem%Np,lmesh%NeA)
    real(RP), intent(out) :: RHOT_dt(elem%Np,lmesh%NeA)
    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)  :: DRHOT_(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)  :: DDENS0_(elem%Np,lmesh%NeA)
    real(RP), intent(in)  :: MOMX0_(elem%Np,lmesh%NeA)
    real(RP), intent(in)  :: MOMY0_(elem%Np,lmesh%NeA)
    real(RP), intent(in)  :: MOMZ0_(elem%Np,lmesh%NeA)
    real(RP), intent(in)  :: DRHOT0_(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)
    class(ElementOperationBase3D), intent(in) :: element3D_operation
    class(SparseMat), intent(in) :: Dz, Lift
    real(RP), intent(in) :: impl_fac
    real(RP), intent(in) :: dt
 
    real(RP) :: PROG_VARS (elem%Np,lmesh%NeZ,PRGVAR_NUM,lmesh%NeX*lmesh%NeY)
    real(RP) :: PROG_VARS0(elem%Np,lmesh%NeZ,PRGVAR_NUM,lmesh%NeX*lmesh%NeY)
    real(RP) :: b1D(3,elem%Nnode_v,lmesh%NeZ,elem%Nnode_h1D**2,lmesh%NeX*lmesh%NeY)
    integer :: ipiv(elem%Nnode_v*3*lmesh%NeZ,elem%Nnode_h1D**2)
    real(RP) :: b1D_uv(elem%Nnode_v,lmesh%NeZ,2,elem%Nnode_h1D**2,lmesh%NeX*lmesh%NeY)
    integer :: ipiv_uv(elem%Nnode_v*1*lmesh%NeZ,elem%Nnode_h1D**2)
    real(RP) :: alph(elem%NfpTot,lmesh%NeZ,lmesh%NeX*lmesh%NeY)
    real(RP) :: Rtot_z(elem%Np,lmesh%NeZ,lmesh%NeX*lmesh%NeY)
    real(RP) :: CPtot_ov_CVtot(elem%Np,lmesh%NeZ,lmesh%NeX*lmesh%NeY)
    real(RP) :: DENS_hyd_z(elem%Np,lmesh%NeZ,lmesh%NeX*lmesh%NeY)
    real(RP) :: PRES_hyd_z(elem%Np,lmesh%NeZ,lmesh%NeX*lmesh%NeY)
    real(RP) :: GnnM_z(elem%Np,lmesh%NeZ,lmesh%NeX*lmesh%NeY)
    real(RP) :: G13_z(elem%Np,lmesh%NeZ,lmesh%NeX*lmesh%NeY)
    real(RP) :: G23_z(elem%Np,lmesh%NeZ,lmesh%NeX*lmesh%NeY)
    real(RP) :: GsqrtV_z(elem%Np,lmesh%NeZ,lmesh%NeX*lmesh%NeY)
    real(RP) :: nz(elem%NfpTot,lmesh%NeZ,lmesh%NeX*lmesh%NeY)
    integer :: vmapM(elem%NfpTot,lmesh%NeZ)
    integer :: vmapP(elem%NfpTot,lmesh%NeZ)
    integer :: ColMask(elem%Nnode_v)
    integer :: ke_xy, ke_z, ke, ke2D, v
    integer :: itr_nlin
    integer :: kl, ku, nz_1D
    integer :: kl_uv, ku_uv, nz_1D_uv
    integer :: ij
    logical :: is_converged
 
    real(RP), allocatable :: PmatBnd(:,:,:)
    real(RP), allocatable :: PmatBnd_uv(:,:,:)
    integer :: info, info_uv
    !------------------------------------------------------------------------
 
    call prof_rapstart( 'hevi_cal_vi_prep', 3)
 
    nz_1d = elem%Nnode_v * 3 * lmesh%NeZ
    kl = ( elem%Nnode_v + 1 ) * 3 - 1
    ku = kl
    nz_1d_uv = elem%Nnode_v * 1 * lmesh%NeZ
    kl_uv = elem%Nnode_v
    ku_uv = kl_uv
    allocate( pmatbnd(2*kl+ku+1,nz_1d,elem%Nnode_h1D**2) )
    allocate( pmatbnd_uv(2*kl_uv+ku_uv+1,nz_1d_uv,elem%Nnode_h1D**2) )
 
    call lmesh%GetVmapZ1D( vmapm, vmapp ) ! (out)
 
    !-
    
    !$omp parallel private( ke_xy, ke_z, ke, ke2D )
    !$omp do collapse(2)
    do ke_xy=1, lmesh%NeX*lmesh%NeY
    do ke_z=1, lmesh%NeZ
      ke = ke_xy + (ke_z-1)*lmesh%NeX*lmesh%NeY
      ke2d = lmesh%EMap3Dto2D(ke)
 
      prog_vars(:,ke_z,dens_vid,ke_xy) = ddens0_(:,ke)
      prog_vars(:,ke_z,momx_vid,ke_xy) = momx0_(:,ke)
      prog_vars(:,ke_z,momy_vid,ke_xy) = momy0_(:,ke)
      prog_vars(:,ke_z,momz_vid,ke_xy) = momz0_(:,ke)
      prog_vars(:,ke_z,rhot_vid,ke_xy) = drhot0_(:,ke)
 
      dens_hyd_z(:,ke_z,ke_xy) = dens_hyd(:,ke)
      pres_hyd_z(:,ke_z,ke_xy) = pres_hyd(:,ke)
 
      rtot_z(:,ke_z,ke_xy) = rtot(:,ke)
      cptot_ov_cvtot(:,ke_z,ke_xy) = cptot(:,ke) / cvtot(:,ke)
      
      nz(:,ke_z,ke_xy) = lmesh%normal_fn(:,ke,3)
      g13_z(:,ke_z,ke_xy) = lmesh%GI3(:,ke,1)
      g23_z(:,ke_z,ke_xy) = lmesh%GI3(:,ke,2)
      gsqrtv_z(:,ke_z,ke_xy) = lmesh%Gsqrt(:,ke) / lmesh%GsqrtH(elem%IndexH2Dto3D,ke2d)
 
      gnnm_z(:,ke_z,ke_xy) = ( &
          1.0_rp / gsqrtv_z(:,ke_z,ke_xy)**2                                                    &
        + g13_z(:,ke_z,ke_xy) * ( lmesh%GIJ(elem%IndexH2Dto3D,ke2d,1,1) * g13_z(:,ke_z,ke_xy)   &
                                + lmesh%GIJ(elem%IndexH2Dto3D,ke2d,1,2) * g23_z(:,ke_z,ke_xy) ) &
        + g23_z(:,ke_z,ke_xy) * ( lmesh%GIJ(elem%IndexH2Dto3D,ke2d,1,2) * g13_z(:,ke_z,ke_xy)   &
                                + lmesh%GIJ(elem%IndexH2Dto3D,ke2d,2,2) * g23_z(:,ke_z,ke_xy) ) )
    end do
    end do
    !$omp end do
    !$omp workshare
    prog_vars0(:,:,:,:) = prog_vars(:,:,:,:)
    !$omp end workshare
    !$omp end parallel
      
    call prof_rapend( 'hevi_cal_vi_prep', 3)
 
    !--
 
    if ( abs(impl_fac) > 0.0_rp ) then
      call prof_rapstart( 'hevi_cal_vi_itr', 3)
 
      ! G = (q^n+1 - q^n*) + impl_fac * A(q^n+1) = 0
      ! dG/dq^n+1 del[q] = - G(q^n*)
      do itr_nlin = 1, 1
        call prof_rapstart( 'hevi_cal_vi_ax', 3)
 
        call vi_eval_ax_uv( &
          momx_dt(:,:), momy_dt(:,:), alph(:,:,:),        & ! (out)
          prog_vars, prog_vars0,                          & ! (in)
          ddens_, momx_, momy_, momz_, drhot_,            & ! (in)
          dens_hyd_z, pres_hyd_z,                         & ! (in)
          rtot_z, cptot_ov_cvtot,                         & ! (in)
          dz, lift, intrpmat_vpordm1,                     & ! (in)
          gnnm_z, g13_z, g23_z, gsqrtv_z,                 & ! (in)
          impl_fac, dt,                                   & ! (in) 
          lmesh, elem, nz, vmapm, vmapp,                  & ! (in)
          b1d_uv(:,:,:,:,:)                               ) ! (out)
 
        call prof_rapend( 'hevi_cal_vi_ax', 3)
 
        do ke_xy=1, lmesh%NeX * lmesh%NeY
          call prof_rapstart( 'hevi_cal_vi_matbnd', 3)
 
          call vi_construct_matbnd_uv( pmatbnd_uv(:,:,:),              & ! (out)
            kl_uv, ku_uv, nz_1d_uv,                                    & ! (in)
            prog_vars(:,:,:,ke_xy),                                    & ! (in)
            dens_hyd_z(:,:,ke_xy), pres_hyd_z(:,:,ke_xy),              & ! (in)
            g13_z(:,:,ke_xy), g23_z(:,:,ke_xy), gsqrtv_z(:,:,ke_xy),   & ! (in)
            alph(:,:,ke_xy),                                           & ! (in)
            rtot_z(:,:,ke_xy), cptot_ov_cvtot(:,:,ke_xy),              & ! (in)
            dz, lift, intrpmat_vpordm1,                                & ! (in)
            impl_fac, dt,                                              & ! (in)
            lmesh, elem, nz(:,:,ke_xy), vmapm, vmapp, ke_xy, 1         ) ! (in)
 
          call prof_rapend( 'hevi_cal_vi_matbnd', 3)
          
          call prof_rapstart( 'hevi_cal_vi_lin', 3)
          !$omp parallel private(ij, v, ke_z, info_uv, ColMask)
          !$omp do
          do ij=1, elem%Nnode_h1D**2
            call dgbsv( nz_1d_uv, kl_uv, ku_uv, 2, pmatbnd_uv(:,:,ij), 2*kl_uv+ku_uv+1, ipiv_uv(:,ij), b1d_uv(:,:,:,ij,ke_xy), nz_1d_uv, info_uv)
 
            colmask(:) = elem%Colmask(:,ij)
            do ke_z=1, lmesh%NeZ
              prog_vars(colmask(:),ke_z,momx_vid,ke_xy) = prog_vars(colmask(:),ke_z,momx_vid,ke_xy) + b1d_uv(:,ke_z,1,ij,ke_xy)
              prog_vars(colmask(:),ke_z,momy_vid,ke_xy) = prog_vars(colmask(:),ke_z,momy_vid,ke_xy) + b1d_uv(:,ke_z,2,ij,ke_xy)
            end do
          end do ! for ij
         !$omp end do
         !$omp end parallel
          call prof_rapend( 'hevi_cal_vi_lin', 3)
        end do ! for ke_xy
 
        call prof_rapstart( 'hevi_cal_vi_ax', 3)
        call vi_eval_ax( &
          dens_dt(:,:), momz_dt(:,:), rhot_dt(:,:),                             & ! (out, dummy) 
          alph(:,:,:),                                                          & ! (in)
          prog_vars, prog_vars0,                                                & ! (in)
          ddens_, momx_, momy_, momz_, drhot_,                                  & ! (in)
          dens_hyd_z, pres_hyd_z,                                               & ! (in)
          rtot_z, cptot_ov_cvtot,                                               & ! (in)
          dz, lift, intrpmat_vpordm1,                                           & ! (in)
          gnnm_z, g13_z, g23_z, gsqrtv_z,                                       & ! (in)
          impl_fac, dt,                                                         & ! (in) 
          lmesh, elem, nz, vmapm, vmapp,                                        & ! (in)
          b1d(:,:,:,:,:)                                                        ) ! (out)
        call prof_rapend( 'hevi_cal_vi_ax', 3)
 
        do ke_xy=1, lmesh%NeX * lmesh%NeY
          call prof_rapstart( 'hevi_cal_vi_matbnd', 3)
          call vi_construct_matbnd( pmatbnd(:,:,:),                    & ! (out)
            kl, ku, nz_1d,                                             & ! (in)
            prog_vars(:,:,:,ke_xy),                                    & ! (in)
            dens_hyd_z(:,:,ke_xy), pres_hyd_z(:,:,ke_xy),              & ! (in)
            g13_z(:,:,ke_xy), g23_z(:,:,ke_xy), gsqrtv_z(:,:,ke_xy),   & ! (in)
            alph(:,:,ke_xy),                                           & ! (in)
            rtot_z(:,:,ke_xy), cptot_ov_cvtot(:,:,ke_xy),              & ! (in)
            dz, lift, intrpmat_vpordm1,                                & ! (in)
            impl_fac, dt,                                              & ! (in)
            lmesh, elem, nz(:,:,ke_xy), vmapm, vmapp, ke_xy, 1         ) ! (in)
 
          call prof_rapend( 'hevi_cal_vi_matbnd', 3)
          
          call prof_rapstart( 'hevi_cal_vi_lin', 3)
          !$omp parallel private(ij, v, ke_z, info, ColMask)
          !$omp do
          do ij=1, elem%Nnode_h1D**2
            call dgbsv( nz_1d, kl, ku, 1, pmatbnd(:,:,ij), 2*kl+ku+1, ipiv(:,ij), b1d(:,:,:,ij,ke_xy), nz_1d, info)
 
            colmask(:) = elem%Colmask(:,ij)
            do ke_z=1, lmesh%NeZ
              prog_vars(colmask(:),ke_z,dens_vid,ke_xy) = prog_vars(colmask(:),ke_z,dens_vid,ke_xy) + b1d(1,:,ke_z,ij,ke_xy)
              prog_vars(colmask(:),ke_z,momz_vid,ke_xy) = prog_vars(colmask(:),ke_z,momz_vid,ke_xy) + b1d(2,:,ke_z,ij,ke_xy)
              prog_vars(colmask(:),ke_z,rhot_vid,ke_xy) = prog_vars(colmask(:),ke_z,rhot_vid,ke_xy) + b1d(3,:,ke_z,ij,ke_xy)
            end do
          end do ! for ij
         !$omp end do
         !$omp end parallel          
        call prof_rapend( 'hevi_cal_vi_lin', 3)
 
        end do ! for ke_xy
      end do ! itr nlin
 
      call prof_rapend( 'hevi_cal_vi_itr', 3)
    end if
 
    call prof_rapstart( 'hevi_cal_vi_retrun_var', 3)
    if ( abs(impl_fac) > 0.0_rp) then
      !$omp parallel do collapse(2) private(ke_xy, ke_z, ke)
      do ke_xy=1, lmesh%NeX * lmesh%NeY
      do ke_z=1, lmesh%NeZ
        ke = ke_xy + (ke_z-1)*lmesh%NeX*lmesh%NeY
        dens_dt(:,ke) = ( prog_vars(:,ke_z,dens_vid,ke_xy) - ddens_(:,ke) ) / impl_fac
        momx_dt(:,ke) = ( prog_vars(:,ke_z,momx_vid,ke_xy) - momx_(:,ke) ) / impl_fac
        momy_dt(:,ke) = ( prog_vars(:,ke_z,momy_vid,ke_xy) - momy_(:,ke) ) / impl_fac
        momz_dt(:,ke) = ( prog_vars(:,ke_z,momz_vid,ke_xy) - momz_(:,ke) ) / impl_fac
        rhot_dt(:,ke) = ( prog_vars(:,ke_z,rhot_vid,ke_xy) - drhot_(:,ke) ) / impl_fac
      end do
      end do
    else
      call vi_eval_ax_uv( & 
        momx_dt(:,:), momy_dt(:,:),                                           & ! (out) 
        alph(:,:,:),                                                          & ! (out, dummy)
        prog_vars, prog_vars0,                                                & ! (in)
        ddens_, momx_, momy_, momz_, drhot_,                                  & ! (in)
        dens_hyd_z, pres_hyd_z,                                               & ! (in)
        rtot_z, cptot_ov_cvtot,                                               & ! (in)
        dz, lift, intrpmat_vpordm1,                                           & ! (in)
        gnnm_z, g13_z, g23_z, gsqrtv_z,                                       & ! (in)
        impl_fac, dt,                                                         & ! (in) 
        lmesh, elem, nz, vmapm, vmapp                                         ) ! (in)
 
      call vi_eval_ax( & 
        dens_dt(:,:), momz_dt(:,:), rhot_dt(:,:), & ! (out) 
        alph(:,:,:),                                                          & ! (in, dummy)
        prog_vars, prog_vars0,                                                & ! (in)
        ddens_, momx_, momy_, momz_, drhot_,                                  & ! (in)
        dens_hyd_z, pres_hyd_z,                                               & ! (in)
        rtot_z, cptot_ov_cvtot,                                               & ! (in)
        dz, lift, intrpmat_vpordm1,                                           & ! (in)
        gnnm_z, g13_z, g23_z, gsqrtv_z,                                       & ! (in)
        impl_fac, dt,                                                         & ! (in) 
        lmesh, elem, nz, vmapm, vmapp                                         ) ! (in)
    end if
    call prof_rapend( 'hevi_cal_vi_retrun_var', 3)
 
    return

References scale_atm_dyn_dgm_nonhydro3d_rhot_hevi_common::atm_dyn_dgm_nonhydro3d_rhot_hevi_common_construct_matbnd(), scale_atm_dyn_dgm_nonhydro3d_rhot_hevi_common::atm_dyn_dgm_nonhydro3d_rhot_hevi_common_construct_matbnd_uv(), scale_atm_dyn_dgm_nonhydro3d_rhot_hevi_common::atm_dyn_dgm_nonhydro3d_rhot_hevi_common_eval_ax(), scale_atm_dyn_dgm_nonhydro3d_rhot_hevi_common::atm_dyn_dgm_nonhydro3d_rhot_hevi_common_eval_ax_uv(), and scale_atm_dyn_dgm_nonhydro3d_common::intrpmat_vpordm1.