scale_atm_dyn_dgm_nonhydro3d_rhot_hevi.F90

Go to the documentation of this file.

!-------------------------------------------------------------------------------
!> module FElib / Fluid dyn solver / Atmosphere / Regional nonhydrostatic model / HEVI
!!
!! @par Description
!!      HEVI DGM scheme for 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
!<
!-------------------------------------------------------------------------------
#include "scaleFElib.h"
module scale_atm_dyn_dgm_nonhydro3d_rhot_hevi
  !-----------------------------------------------------------------------------
  !
  !++ Used modules
  !
  use scale_precision
  use scale_io
  use scale_prc
  use scale_prof
  use scale_const, only: &
    grav => const_grav,  &
    rdry => const_rdry,  &
    cpdry => const_cpdry, &
    cvdry => const_cvdry, &
    pres00 => const_pre00
 
  use scale_sparsemat  
  use scale_element_base, only: &
    elementbase2d, elementbase3d
  use scale_element_modalfilter, only: modalfilter
  use scale_element_hexahedral, only: hexahedralelement
  use scale_element_operation_base, only: elementoperationbase3d
  use scale_localmesh_2d, only: localmesh2d  
  use scale_localmesh_3d, only: localmesh3d
  use scale_mesh_base3d, only: meshbase3d
  use scale_localmeshfield_base, only: localmeshfield3d
  use scale_meshfield_base, only: meshfield3d
 
  use scale_atm_dyn_dgm_nonhydro3d_common, only: &
    atm_dyn_dgm_nonhydro3d_common_init,                       &
    atm_dyn_dgm_nonhydro3d_common_final,                      &
    dens_vid => prgvar_ddens_id, rhot_vid => prgvar_drhot_id, &
    momx_vid => prgvar_momx_id, momy_vid => prgvar_momy_id,   &
    momz_vid => prgvar_momz_id,                               &
    prgvar_num, intrpmat_vpordm1
 
  !-----------------------------------------------------------------------------
  implicit none
  private
  !-----------------------------------------------------------------------------
  !
  !++ Public procedures
  !
  public :: atm_dyn_dgm_nonhydro3d_rhot_hevi_init
  public :: atm_dyn_dgm_nonhydro3d_rhot_hevi_final
  public :: atm_dyn_dgm_nonhydro3d_rhot_hevi_cal_tend
  public :: atm_dyn_dgm_nonhydro3d_rhot_hevi_cal_tend_asis
  public :: atm_dyn_dgm_nonhydro3d_rhot_hevi_cal_vi
  public :: atm_dyn_dgm_nonhydro3d_rhot_hevi_cal_vi_asis
 
  !-----------------------------------------------------------------------------
  !
  !++ Public parameters & variables
  !
  
  !-----------------------------------------------------------------------------
  !
  !++ Private procedures & variables
  !
  !-------------------
 
contains
  subroutine atm_dyn_dgm_nonhydro3d_rhot_hevi_init( mesh )
    use scale_atm_dyn_dgm_nonhydro3d_rhot_hevi_common, only: &
      atm_dyn_dgm_nonhydro3d_rhot_hevi_common_init
    implicit none
    class(meshbase3d), intent(in) :: mesh
    !--------------------------------------------
 
    call atm_dyn_dgm_nonhydro3d_common_init( mesh )
    call atm_dyn_dgm_nonhydro3d_rhot_hevi_common_init()
 
    return
  end subroutine atm_dyn_dgm_nonhydro3d_rhot_hevi_init
 
 
  subroutine atm_dyn_dgm_nonhydro3d_rhot_hevi_final()
    use scale_atm_dyn_dgm_nonhydro3d_rhot_hevi_common, only: &
      atm_dyn_dgm_nonhydro3d_rhot_hevi_common_final
    implicit none
    !--------------------------------------------
 
    call atm_dyn_dgm_nonhydro3d_rhot_hevi_common_final()
    call atm_dyn_dgm_nonhydro3d_common_final()
 
    return
  end subroutine atm_dyn_dgm_nonhydro3d_rhot_hevi_final  
 
  !-------------------------------
 
!OCL SERIAL  
  subroutine atm_dyn_dgm_nonhydro3d_rhot_hevi_cal_tend_asis( &
    DENS_dt, MOMX_dt, MOMY_dt, MOMZ_dt, RHOT_dt,         & ! (out)
    ddens_, momx_, momy_, momz_, drhot_, dpres_,         & ! (in)
    dens_hyd, pres_hyd, pres_hyd_ref, therm_hyd,         & ! (in)
    coriolis, rtot, cvtot, cptot, dphyddx, dphyddy,      & ! (in)
    element3d_operation, dx, dy, dz, sx, sy, sz, lift,   & ! (in)
    lmesh, elem, lmesh2d, elem2d )                         ! (in)
 
    use scale_atm_dyn_dgm_nonhydro3d_rhot_hevi_numflux, only: &
      get_ebnd_flux => atm_dyn_dgm_nonhydro3d_rhot_hevi_numflux_get_generalvc_asis
 
    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)  :: therm_hyd(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) :: gsqrtv(elem%np), rgsqrtv(elem%np)
    real(rp) :: cori(elem%np)
 
    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%GI3(:,:,1), lmesh%GI3(:,:,2),                        & ! (in)    
      lmesh%normal_fn(:,:,1), lmesh%normal_fn(:,:,2), lmesh%normal_fn(:,:,3), & ! (in)
      lmesh%vmapM, lmesh%vmapP,                                               & ! (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
 
    !$omp parallel do private( ke2d, Cori,         &
    !$omp RHOT_, rdens_, u_, v_, w_, wt_,          &
    !$omp DPRES_hyd, GradPhyd_x, GradPhyd_y,       &
    !$omp GsqrtV, RGsqrtV,                         &
    !$omp Fx, Fy, Fz, LiftDelFlx )
    do ke = lmesh%NeS, lmesh%NeE
      !--
      ke2d = lmesh%EMap3Dto2D(ke)
      cori(:) = coriolis(elem%IndexH2Dto3D(:),ke2d)
 
      gsqrtv(:)  = lmesh%Gsqrt(:,ke) / 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_(:) 
 
      !-- 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) + dpres_(:,ke) ), fx)
      call sparsemat_matmul(dy, lmesh%Gsqrt(:,ke) *    v_(:) * momx_(:,ke)                 , fy)
      call sparsemat_matmul(dz, lmesh%Gsqrt(:,ke) * ( wt_(:) * momx_(:,ke)                     &
                                                    + lmesh%GI3(:,ke,1) * 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) &
          - gradphyd_x(:) * rgsqrtv(:)            &
          + cori(:) * momy_(:,ke)
 
      !-- MOMY
      call sparsemat_matmul(dx, lmesh%Gsqrt(:,ke) *    u_(:) * momy_(:,ke)                 , fx)
      call sparsemat_matmul(dy, lmesh%Gsqrt(:,ke) * (  v_(:) * momy_(:,ke) + dpres_(:,ke) ), fy)
      call sparsemat_matmul(dz, lmesh%Gsqrt(:,ke) * ( wt_(:) * momy_(:,ke)                     &
                                                    + lmesh%GI3(:,ke,2) * 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) &
          - gradphyd_y(:) * rgsqrtv(:)            &
          - cori(:) * momx_(:,ke)
 
      !-- 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
    call prof_rapend( 'cal_dyn_tend_interior', 3)
 
    return
  end subroutine atm_dyn_dgm_nonhydro3d_rhot_hevi_cal_tend_asis
 
!OCL SERIAL
  subroutine atm_dyn_dgm_nonhydro3d_rhot_hevi_cal_tend( &
    DENS_dt, MOMX_dt, MOMY_dt, MOMZ_dt, RHOT_dt,         & ! (out)
    ddens_, momx_, momy_, momz_, drhot_, dpres_,         & ! (in)
    dens_hyd, pres_hyd, pres_hyd_ref, therm_hyd,         & ! (in)
    coriolis, rtot, cvtot, cptot, dphyddx, dphyddy,      & ! (in)
    element3d_operation, dx, dy, dz, sx, sy, sz, lift,   & ! (in)
    lmesh, elem, lmesh2d, elem2d )                         ! (in)
 
    use scale_atm_dyn_dgm_nonhydro3d_rhot_hevi_numflux, only: &
      get_ebnd_flux => atm_dyn_dgm_nonhydro3d_rhot_hevi_numflux_get_generalvc
    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)  :: therm_hyd(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) :: flux(elem%np,3,5), dflux(elem%np,4,5)
    real(rp) :: del_flux(elem%nfptot,prgvar_num,lmesh%ne)
    real(rp) :: u_, v_, w_, pt_
    real(rp) :: rhot_(elem%np)
    real(rp) :: rdens_(elem%np), gsqrtv(elem%np), rgsqrtv(elem%np), rgsqrt(elem%np)
    real(rp) :: gsqrt_, gsqrtdpres_, e11, e22, e33
    real(rp) :: cor
 
    integer :: ke, ke2d
    integer :: p
 
    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,                                                               & ! (out)
      ddens_, momx_, momy_, momz_, drhot_, dpres_,                            & ! (in)
      dens_hyd, pres_hyd, therm_hyd, rtot, cvtot, cptot,                      & ! (in)
      lmesh%Gsqrt, lmesh%GI3(:,:,1), lmesh%GI3(:,:,2),                        & ! (in)    
      lmesh%normal_fn(:,:,1), lmesh%normal_fn(:,:,2), lmesh%normal_fn(:,:,3), & ! (in)
      lmesh%vmapM, lmesh%vmapP,                                               & ! (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
 
    !$omp parallel private( ke, ke2d, p,           &
    !$omp u_, v_, w_, pt_,                         &
    !$omp cor, GsqrtDPRES_,                        &
    !$omp RDENS_, RGsqrt, GsqrtV, RGsqrtV, Gsqrt_, &
    !$omp E11, E22, E33, DFlux, Flux    )
 
    !$omp do
    do ke = lmesh%NeS, lmesh%NeE
      !--
      ke2d = lmesh%EMap3Dto2D(ke)
 
      do p=1, elem%Np
        gsqrtv(p)  = lmesh%Gsqrt(p,ke) / lmesh%GsqrtH(elem%IndexH2Dto3D(p),ke2d)
        rgsqrtv(p) = 1.0_rp / gsqrtv(p)
        rgsqrt(p) = 1.0_rp / lmesh%Gsqrt(p,ke)
        rdens_(p) = 1.0_rp / ( ddens_(p,ke) + dens_hyd(p,ke) )
      end do
 
      !--
      do p=1, elem%Np
        gsqrt_ = lmesh%Gsqrt(p,ke)
        flux(p,1,dens_vid) = gsqrt_ * momx_(p,ke)
        flux(p,2,dens_vid) = gsqrt_ * momy_(p,ke)
        flux(p,3,dens_vid) = gsqrt_ * ( &
            momz_(p,ke) * rgsqrtv(p)        &
          + lmesh%GI3(p,ke,1) * momx_(p,ke) &
          + lmesh%GI3(p,ke,2) * momy_(p,ke) )
      end do
      do p=1, elem%Np
        pt_ = ( therm_hyd(p,ke) + drhot_(p,ke) ) * rdens_(p)
        flux(p,1,rhot_vid) = flux(p,1,dens_vid) * pt_
        flux(p,2,rhot_vid) = flux(p,2,dens_vid) * pt_
!        Flux(p,3,RHOT_VID) = 0.0_RP ! for HEVI
 
        w_ = momz_(p,ke) * rdens_(p)
        flux(p,1,momz_vid) = flux(p,1,dens_vid) * w_
        flux(p,2,momz_vid) = flux(p,2,dens_vid) * w_
        flux(p,3,momz_vid) = flux(p,3,dens_vid) * w_
      end do
 
      do p=1, elem%Np
        gsqrtdpres_ = lmesh%Gsqrt(p,ke) * dpres_(p,ke)
        u_ = momx_(p,ke) * rdens_(p)
        v_ = momy_(p,ke) * rdens_(p)
 
        flux(p,1,momx_vid) = flux(p,1,dens_vid) * u_ + gsqrtdpres_
        flux(p,2,momx_vid) = flux(p,2,dens_vid) * u_
        flux(p,3,momx_vid) = flux(p,3,dens_vid) * u_ + gsqrtdpres_ * lmesh%GI3(p,ke,1)
 
        flux(p,1,momy_vid) = flux(p,1,dens_vid) * v_
        flux(p,2,momy_vid) = flux(p,2,dens_vid) * v_ + gsqrtdpres_ 
        flux(p,3,momy_vid) = flux(p,3,dens_vid) * v_ + gsqrtdpres_ * lmesh%GI3(p,ke,2)
      end do
 
      call element3d_operation%Div_var5( &
        flux, del_flux(:,:,ke), & ! (in)
        dflux                   ) ! (out)
 
      do p=1, elem%Np
        e11 = lmesh%Escale(p,ke,1,1)
        e22 = lmesh%Escale(p,ke,2,2)
        e33 = lmesh%Escale(p,ke,3,3)
 
        dens_dt(p,ke) = - ( &
              e11 * dflux(p,1,dens_vid)    &
            + e22 * dflux(p,2,dens_vid)    &
            + dflux(p,4,dens_vid) ) * rgsqrt(p)
          
        rhot_dt(p,ke) = - ( &
              e11 * dflux(p,1,rhot_vid)    &
            + e22 * dflux(p,2,rhot_vid)    &
            + dflux(p,4,rhot_vid) ) * rgsqrt(p)
      end do
 
      do p=1, elem%Np
        e11 = lmesh%Escale(p,ke,1,1)
        e22 = lmesh%Escale(p,ke,2,2)
        e33 = lmesh%Escale(p,ke,3,3)
            
        momz_dt(p,ke) = - ( &
              e11 * dflux(p,1,momz_vid)     &
            + e22 * dflux(p,2,momz_vid)     &
            + e33 * dflux(p,3,momz_vid)     &
            + dflux(p,4,momz_vid) ) * rgsqrt(p)
      end do
 
      !--
      do p=1, elem%Np
        cor = coriolis(elem%IndexH2Dto3D(p),ke2d)
        momx_dt(p,ke) = - dphyddx(p,ke) &
                        + cor * momy_(p,ke)
        momy_dt(p,ke) = - dphyddy(p,ke) &
                        - cor * momx_(p,ke)
      end do
 
      do p=1, elem%Np
        e11 = lmesh%Escale(p,ke,1,1)
        e22 = lmesh%Escale(p,ke,2,2)
        e33 = lmesh%Escale(p,ke,3,3)
  
        momx_dt(p,ke) = momx_dt(p,ke) - ( &
              e11 * dflux(p,1,momx_vid)    &
            + e22 * dflux(p,2,momx_vid)    &
            + e33 * dflux(p,3,momx_vid)    &
            + dflux(p,4,momx_vid) ) * rgsqrt(p)
        
        momy_dt(p,ke) = momy_dt(p,ke) - ( &
              e11 * dflux(p,1,momy_vid)    &
            + e22 * dflux(p,2,momy_vid)    &
            + e33 * dflux(p,3,momy_vid)    &
            + dflux(p,4,momy_vid) ) * rgsqrt(p)
      end do
    end do
    !$omp end parallel
    call prof_rapend('cal_dyn_tend_interior', 3)
    return
  end subroutine atm_dyn_dgm_nonhydro3d_rhot_hevi_cal_tend
 
  !------
 
!OCL SERIAL  
  subroutine atm_dyn_dgm_nonhydro3d_rhot_hevi_cal_vi_asis( &
    DENS_dt, MOMX_dt, MOMY_dt, MOMZ_dt, RHOT_dt,             & ! (out)
    ddens_, momx_, momy_, momz_, drhot_, dens_hyd, pres_hyd, & ! (in)
    ddens0_, momx0_, momy0_, momz0_, drhot0_,                & ! (in)
    rtot, cvtot, cptot,                                      & ! (in)
    element3d_operation, dz, lift,                           & ! (in)
    impl_fac, dt,                                            & ! (in)
    lmesh, elem, lmesh2d, elem2d                             ) ! (in)
 
    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, &
      vi_use_lapack_flag
    use scale_linalgebra, only: linalgebra_solvelineq_bndmat
    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)**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 linalgebra_solvelineq_bndmat( pmatbnd_uv(:,:,ij), b1d_uv(:,:,:,ij,ke_xy), ipiv_uv(:,ij), nz_1d_uv, kl_uv, ku_uv, 2, vi_use_lapack_flag )
 
            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 linalgebra_solvelineq_bndmat( pmatbnd(:,:,ij), b1d(:,:,:,ij,ke_xy), ipiv(:,ij), nz_1d, kl, ku, 1, vi_use_lapack_flag )
            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
  end subroutine atm_dyn_dgm_nonhydro3d_rhot_hevi_cal_vi_asis
 
 
!OCL SERIAL  
  subroutine atm_dyn_dgm_nonhydro3d_rhot_hevi_cal_vi( &
    DENS_dt, MOMX_dt, MOMY_dt, MOMZ_dt, RHOT_dt,             & ! (out)
    ddens_, momx_, momy_, momz_, drhot_, dens_hyd, pres_hyd, & ! (in)
    ddens0_, momx0_, momy0_, momz0_, drhot0_,                & ! (in)
    rtot, cvtot, cptot,                                      & ! (in)
    element3d_operation, dz, lift,                           & ! (in)
    impl_fac, dt,                                            & ! (in)
    lmesh, elem, lmesh2d, elem2d                             ) ! (in)
 
    use scale_atm_dyn_dgm_hevi_common_linalgebra, only: &
      atm_dyn_dgm_hevi_common_linalgebra_get_param
    use scale_atm_dyn_dgm_nonhydro3d_rhot_hevi_common_2, 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_solve => atm_dyn_dgm_nonhydro3d_rhot_hevi_common_solve,           &
      vi_solve_uv => atm_dyn_dgm_nonhydro3d_rhot_hevi_common_solve_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%nex*lmesh%ney,lmesh%nez,prgvar_num)
    real(rp) :: prog_vars0(elem%np,lmesh%nex*lmesh%ney,lmesh%nez,prgvar_num)
    real(rp) :: alph(elem%nfptot,lmesh%ne)
    real(rp) :: gsqrtv(elem%np,lmesh%ne)
 
    integer :: vmapm(elem%nfptot,lmesh%ne)
    integer :: vmapp(elem%nfptot,lmesh%ne)
    integer :: ke_xy, ke_z, ke, ke2d
    integer :: p
    integer :: itr_nlin
    integer :: ij, i, j, im, jm
    logical :: is_converged
 
    real(rp), allocatable :: b_uv(:,:,:,:)
    real(rp), allocatable :: b   (:,:,:)
 
    real(rp) :: dens(elem%np,lmesh%ne)
    real(rp) :: w(elem%np,lmesh%ne)
    real(rp) :: wt(elem%np,lmesh%ne)
    real(rp) :: pot(elem%np,lmesh%ne)
    real(rp) :: dpdrhot(elem%np,lmesh%ne)
    !------------------------------------------------------------------------
 
    call prof_rapstart( 'hevi_cal_vi_prep', 3)
 
    call lmesh%GetVmapZ3D( vmapm, vmapp ) ! (out)
    call atm_dyn_dgm_hevi_common_linalgebra_get_param( im, jm, & ! (out)
      elem%Nnode_h1D ) ! (in)
    
    !-
    !$omp parallel private( ke_xy, ke_z, ke, ke2D, p )
    !$omp do collapse(2)
    do ke_z =1, lmesh%NeZ
    do ke_xy=1, lmesh%NeX * lmesh%NeY
      ke = ke_xy + (ke_z-1)*lmesh%NeX*lmesh%NeY
      ke2d = lmesh%EMap3Dto2D(ke)
 
      prog_vars(:,ke_xy,ke_z,dens_vid) = ddens0_(:,ke)
      prog_vars(:,ke_xy,ke_z,momx_vid) = momx0_(:,ke)
      prog_vars(:,ke_xy,ke_z,momy_vid) = momy0_(:,ke)
      prog_vars(:,ke_xy,ke_z,momz_vid) = momz0_(:,ke)
      prog_vars(:,ke_xy,ke_z,rhot_vid) = drhot0_(:,ke)
 
      do p=1, elem%Np
        gsqrtv(p,ke)  = lmesh%Gsqrt(p,ke) / lmesh%GsqrtH(elem%IndexH2Dto3D(p),ke2d)
      end do
    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)
 
      allocate( b_uv(im*elem%Nnode_v,2,jm,lmesh%Ne) )
      allocate( b(im*3*elem%Nnode_v,jm,lmesh%Ne) )
 
      ! 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_uv', 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, pres_hyd,                           & ! (in)
          rtot, cptot, cvtot, gsqrtv,                   & ! (in)
          impl_fac, dt,                                 & ! (in) 
          lmesh, elem, vmapm, vmapp,                    & ! (in)
          element3d_operation,                          & ! (in)
          im, jm,                                       & ! (in)
          b_uv                                          ) ! (out)
        call prof_rapend( 'hevi_cal_vi_ax_uv', 3)
 
        call vi_solve_uv( prog_vars,        & ! (inout)
          b_uv, alph,                       & ! (in)
          gsqrtv,                           & ! (in)
          impl_fac, dt,                     & ! (in)
          vmapm, vmapp, lmesh, elem, im, jm ) ! (in)
 
        !----
 
        call prof_rapstart( 'hevi_cal_vi_ax', 3)
        call vi_eval_ax( &
          dens_dt(:,:), momz_dt(:,:), rhot_dt(:,:),                    & ! (out) 
          alph(:,:),                                                   & ! (in)
          prog_vars, prog_vars0, ddens_, momx_, momy_, momz_, drhot_,  & ! (in)
          dens_hyd, pres_hyd, rtot, cptot, cvtot, gsqrtv,              & ! (in)
          impl_fac, dt, lmesh, elem, vmapm, vmapp,                     & ! (in)
          element3d_operation,                                         & ! (in)
          im, jm,                                                      & ! (in)
          b, dens, w, wt, pot, dpdrhot                                 ) ! (out)
        call prof_rapend( 'hevi_cal_vi_ax', 3)
 
        call vi_solve( prog_vars,                                      & ! (inout)
          b, dens, w, wt, pot, dpdrhot, alph,                          & ! (in)
          gsqrtv, intrpmat_vpordm1,                                    & ! (in)
          impl_fac, dt,                                                & ! (in)
          vmapm, vmapp, lmesh, elem, im, jm                            ) ! (in)
      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_z=1, lmesh%NeZ
      do ke_xy=1, lmesh%NeX * lmesh%NeY
        ke = ke_xy + (ke_z-1)*lmesh%NeX*lmesh%NeY
        dens_dt(:,ke) = ( prog_vars(:,ke_xy,ke_z,dens_vid) - ddens_(:,ke) ) / impl_fac
        momx_dt(:,ke) = ( prog_vars(:,ke_xy,ke_z,momx_vid) - momx_(:,ke) ) / impl_fac
        momy_dt(:,ke) = ( prog_vars(:,ke_xy,ke_z,momy_vid) - momy_(:,ke) ) / impl_fac
        momz_dt(:,ke) = ( prog_vars(:,ke_xy,ke_z,momz_vid) - momz_(:,ke) ) / impl_fac
        rhot_dt(:,ke) = ( prog_vars(:,ke_xy,ke_z,rhot_vid) - drhot_(:,ke) ) / impl_fac
      end do
      end do
    else
      call vi_eval_ax_uv( &
        momx_dt(:,:), momy_dt(:,:), alph(:,:),        & ! (out)
        prog_vars, prog_vars0,                        & ! (in)
        ddens_, momx_, momy_, momz_, drhot_,          & ! (in)
        dens_hyd, pres_hyd,                           & ! (in)
        rtot, cptot, cvtot, gsqrtv,                   & ! (in)
        impl_fac, dt,                                 & ! (in) 
        lmesh, elem, vmapm, vmapp,                    & ! (in)
        element3d_operation, im, jm                   ) ! (in)
 
      call vi_eval_ax( &
        dens_dt(:,:), momz_dt(:,:), rhot_dt(:,:),                    & ! (out, dummy) 
        alph(:,:),                                                   & ! (in)
        prog_vars, prog_vars0, ddens_, momx_, momy_, momz_, drhot_,  & ! (in)
        dens_hyd, pres_hyd, rtot, cptot, cvtot, gsqrtv,              & ! (in)
        impl_fac, dt, lmesh, elem, vmapm, vmapp,                     & ! (in)
        element3d_operation, im, jm                                  ) ! (in)
    end if
    call prof_rapend( 'hevi_cal_vi_retrun_var', 3)
 
    return
  end subroutine atm_dyn_dgm_nonhydro3d_rhot_hevi_cal_vi  
 
end module scale_atm_dyn_dgm_nonhydro3d_rhot_hevi