scale_atm_dyn_dgm_nonhydro3d_rhot_hevi_splitform.F90

Go to the documentation of this file.

!-------------------------------------------------------------------------------
!> module Atmosphere / Dynamics HEVI 
!!
!! @par Description
!!      HEVI DGM scheme for Atmospheric dynamical process. 
!!      To improve the numerical instability due to the aliasing errors, 
!!      the split form based on Gassner et al. (2016, JCP) is used for advection terms. 
!!
!! @author Yuta Kawai, Team SCALE
!<
!-------------------------------------------------------------------------------
#include "scaleFElib.h"
module scale_atm_dyn_dgm_nonhydro3d_rhot_hevi_splitform
  !-----------------------------------------------------------------------------
  !
  !++ 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_modalfilter, only: modalfilter
  use scale_element_base, only: &
    elementbase2d, elementbase3d
  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_splitform_init
  public :: atm_dyn_dgm_nonhydro3d_rhot_hevi_splitform_final
  public :: atm_dyn_dgm_nonhydro3d_rhot_hevi_splitform_cal_tend
  public :: atm_dyn_dgm_nonhydro3d_rhot_hevi_splitform_cal_vi
 
  !-----------------------------------------------------------------------------
  !
  !++ Public parameters & variables
  !
  
  !-----------------------------------------------------------------------------
  !
  !++ Private procedures & variables
  !
  !-------------------
 
  real(RP), private, allocatable :: DxT1D_(:,:)
  real(RP), private, allocatable :: DyT1D_(:,:)
  real(RP), private, allocatable :: DzT1D_(:,:)
 
  private :: dx_ab, dy_ab, dz_ab
  private :: dx_abc, dy_abc, dz_abc
 
contains
  subroutine atm_dyn_dgm_nonhydro3d_rhot_hevi_splitform_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
 
    integer :: p1, p2, p3, p_
    type(elementbase3d), pointer :: elem
    !--------------------------------------------
 
    call atm_dyn_dgm_nonhydro3d_common_init( mesh )
    call atm_dyn_dgm_nonhydro3d_rhot_hevi_common_init()
 
    !-
    elem => mesh%refElem3D
 
    allocate( dxt1d_(elem%Nnode_h1D,elem%Nnode_h1D) )
    allocate( dyt1d_(elem%Nnode_h1D,elem%Nnode_h1D) )
    allocate( dzt1d_(elem%Nnode_v,elem%Nnode_v) )
 
    do p1=1, elem%Nnode_h1D
      dxt1d_(:,p1) = elem%Dx1(p1,1:elem%Nnode_h1D)
    end do
 
    do p2=1, elem%Nnode_h1D
    do p_=1, elem%Nnode_h1D
      dyt1d_(p_,p2) = elem%Dx2(1+(p2-1)*elem%Nnode_h1D,1+(p_-1)*elem%Nnode_h1D)
    end do
    end do
 
    do p3=1, elem%Nnode_v
      dzt1d_(:,p3) = elem%Dx3(elem%Colmask(p3,1),elem%Colmask(:,1))
    end do
 
    return
  end subroutine atm_dyn_dgm_nonhydro3d_rhot_hevi_splitform_init
 
 
  subroutine atm_dyn_dgm_nonhydro3d_rhot_hevi_splitform_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()
 
    deallocate( dxt1d_, dyt1d_, dzt1d_ )
 
    return
  end subroutine atm_dyn_dgm_nonhydro3d_rhot_hevi_splitform_final  
 
  !-------------------------------
 
  subroutine atm_dyn_dgm_nonhydro3d_rhot_hevi_splitform_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_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) :: fx_sp(elem%np), fy_sp(elem%np), fz_sp(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) :: gsqrtdens_(elem%np), rdens_(elem%np), rhot_hyd(elem%np), rhot_(elem%np)
    real(rp) :: u_(elem%np), v_(elem%np), w_(elem%np), wt_(elem%np), pot_(elem%np)
    real(rp) :: cori(elem%np)
    real(rp) :: gsqrtv(elem%np), rgsqrtv(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_, GsqrtDens_, rdens_, u_, v_, w_, wt_, pot_,         &
    !$omp DPRES_hyd, GradPhyd_x, GradPhyd_y,                        &
    !$omp GsqrtV, RGsqrtV,                                          &
    !$omp Fx, Fy, Fz, Fx_sp, Fy_sp, Fz_sp, 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)
 
      gsqrtdens_(:) = lmesh%Gsqrt(:,ke) * ( ddens_(:,ke) + dens_hyd(:,ke) )
      rdens_(:) = 1.0_rp / gsqrtdens_(:)
      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_(:) 
      pot_(:) = rhot_(:) * rdens_(:)
 
      ke2d = lmesh%EMap3Dto2D(ke)
      cori(:) = coriolis(elem%IndexH2Dto3D(:),ke2d)
 
      !-- 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 dx_ab( dxt1d_, gsqrtdens_(:), u_, elem%Nnode_h1D, elem%Nnode_v, fx_sp )
      call dy_ab( dyt1d_, gsqrtdens_(:), v_, elem%Nnode_h1D, elem%Nnode_v, fy_sp )
      call dz_ab( dzt1d_, gsqrtdens_(:), wt_(:) - w_(:) * rgsqrtv(:), elem%Nnode_h1D, elem%Nnode_v, fz_sp )
      call sparsemat_matmul(lift, lmesh%Fscale(:,ke) * del_flux(:,ke,dens_vid), liftdelflx)
 
      dens_dt(:,ke) = - ( &
            lmesh%Escale(:,ke,1,1) * fx_sp(:) &
          + lmesh%Escale(:,ke,2,2) * fy_sp(:) &
          + lmesh%Escale(:,ke,3,3) * fz_sp(:) &
          + liftdelflx(:) )
      
      !-- MOMX
      call dx_abc( dxt1d_, gsqrtdens_, u_,  u_, elem%Nnode_h1D, elem%Nnode_v, fx_sp )
      call dy_abc( dyt1d_, gsqrtdens_, u_,  v_, elem%Nnode_h1D, elem%Nnode_v, fy_sp )
      call dz_abc( dzt1d_, gsqrtdens_, u_, wt_, elem%Nnode_h1D, elem%Nnode_v, fz_sp )
      call sparsemat_matmul(dx, lmesh%Gsqrt(:,ke) * dpres_(:,ke)              , fx)
      call sparsemat_matmul(lift, lmesh%Fscale(:,ke) * del_flux(:,ke,momx_vid), liftdelflx)
 
      momx_dt(:,ke) = &
          - ( lmesh%Escale(:,ke,1,1) * ( fx_sp(:) + fx(:) )  &
            + lmesh%Escale(:,ke,2,2) * fy_sp(:)              &
            + lmesh%Escale(:,ke,3,3) * fz_sp(:)              &
            + liftdelflx(:) ) / lmesh%Gsqrt(:,ke)            &
          - gradphyd_x(:) * rgsqrtv(:)                       &
          + cori(:) * momy_(:,ke)
 
      !-- MOMY
      call dx_abc( dxt1d_, gsqrtdens_, v_,  u_, elem%Nnode_h1D, elem%Nnode_v, fx_sp )
      call dy_abc( dyt1d_, gsqrtdens_, v_,  v_, elem%Nnode_h1D, elem%Nnode_v, fy_sp )  
      call dz_abc( dzt1d_, gsqrtdens_, v_, wt_, elem%Nnode_h1D, elem%Nnode_v, fz_sp )
      call sparsemat_matmul(dy, lmesh%Gsqrt(:,ke) * dpres_(:,ke)              , fy)
      call sparsemat_matmul(lift, lmesh%Fscale(:,ke) * del_flux(:,ke,momy_vid), liftdelflx)
 
      momy_dt(:,ke) = &
          - ( lmesh%Escale(:,ke,1,1) * fx_sp(:)             &
            + lmesh%Escale(:,ke,2,2) * ( fy_sp(:) + fy(:) ) &
            + lmesh%Escale(:,ke,3,3) * fz_sp(:)             &
            + liftdelflx(:) ) / lmesh%Gsqrt(:,ke)           &
          - gradphyd_y(:) * rgsqrtv(:)                      &
          - cori(:) * momx_(:,ke)
 
      !-- MOMZ
      call dx_abc( dxt1d_, gsqrtdens_, w_,  u_, elem%Nnode_h1D, elem%Nnode_v, fx_sp )
      call dy_abc( dyt1d_, gsqrtdens_, w_,  v_, elem%Nnode_h1D, elem%Nnode_v, fy_sp )
      call dz_abc( dzt1d_, gsqrtdens_, w_, wt_, elem%Nnode_h1D, elem%Nnode_v, fz_sp )
      call sparsemat_matmul(lift, lmesh%Fscale(:,ke) * del_flux(:,ke,momz_vid), liftdelflx)
      
      momz_dt(:,ke) = - ( &
            lmesh%Escale(:,ke,1,1) * fx_sp(:)   &
          + lmesh%Escale(:,ke,2,2) * fy_sp(:)   &
          + lmesh%Escale(:,ke,3,3) * fz_sp(:)   &
          + liftdelflx(:)                       )   
 
      !-- RHOT
      call dx_abc( dxt1d_, gsqrtdens_, pot_,  u_, elem%Nnode_h1D, elem%Nnode_v, fx_sp )
      call dy_abc( dyt1d_, gsqrtdens_, pot_,  v_, elem%Nnode_h1D, elem%Nnode_v, fy_sp )        
      call dz_abc( dzt1d_, gsqrtdens_, pot_,  wt_(:) - w_(:) * rgsqrtv(:), elem%Nnode_h1D, elem%Nnode_v, fz_sp )        
      call sparsemat_matmul(lift, lmesh%Fscale(:,ke) * del_flux(:,ke,rhot_vid), liftdelflx)
      
      rhot_dt(:,ke) =  - (  &
            lmesh%Escale(:,ke,1,1) * fx_sp(:) &
          + lmesh%Escale(:,ke,2,2) * fy_sp(:) &
          + lmesh%Escale(:,ke,3,3) * fz_sp(:) &
          + liftdelflx(:)                     ) 
    end do
    call prof_rapend( 'cal_dyn_tend_interior', 3)
 
    return
  end subroutine atm_dyn_dgm_nonhydro3d_rhot_hevi_splitform_cal_tend
 
  !------
 
  subroutine atm_dyn_dgm_nonhydro3d_rhot_hevi_splitform_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_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(elem%nnode_v,3,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, info
    logical :: is_converged
 
    real(rp), allocatable :: pmatbnd(:,:,:)
    real(rp), allocatable :: pmatbnd_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, 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)
            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 dgbsv( nz_1D, kl, ku, 1, PmatBnd(:,:,ij), 2*kl+ku+1, ipiv(:,ij), b1D(:,:,:,ij,ke_xy), nz_1D, info)
            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_splitform_cal_vi
 
  !-----------------------------------------
 
  subroutine dx_ab(DxT1D, a, b, Nnode_h1D, Nnode_v, fx)
    integer, intent(in) :: nnode_h1d, nnode_v
    real(rp), intent(in) :: dxt1d(nnode_h1d,nnode_h1d)
    real(rp), intent(in) :: a(nnode_h1d,nnode_h1d,nnode_v)
    real(rp), intent(in) :: b(nnode_h1d,nnode_h1d,nnode_v)
    real(rp), intent(out) :: fx(nnode_h1d,nnode_h1d,nnode_v)
 
    integer :: i, j, k, p
    integer :: i2, j2, k2, p2
    !-------------------------------------------------
 
    do k=1, nnode_v
    do j=1, nnode_h1d    
    do i=1, nnode_h1d
      fx(i,j,k) = 0.5_rp * sum( dxt1d(:,i) * (a(i,j,k) + a(:,j,k)) * (b(i,j,k) + b(:,j,k)) )
    end do
    end do
    end do
 
    return
  end subroutine dx_ab
 
  subroutine dy_ab(DyT1D, a, b, Nnode_h1D, Nnode_v, fy)
    integer, intent(in) :: nnode_h1d, nnode_v
    real(rp), intent(in) :: dyt1d(nnode_h1d,nnode_h1d)
    real(rp), intent(in) :: a(nnode_h1d,nnode_h1d,nnode_v)
    real(rp), intent(in) :: b(nnode_h1d,nnode_h1d,nnode_v)
    real(rp), intent(out) :: fy(nnode_h1d,nnode_h1d,nnode_v)
 
    integer :: i, j, k, p
    integer :: i2, j2, k2, p2
    !-------------------------------------------------
 
    do k=1, nnode_v
    do j=1, nnode_h1d    
    do i=1, nnode_h1d
      fy(i,j,k) = 0.5_rp * sum( dyt1d(:,j) * (a(i,j,k) + a(i,:,k)) * (b(i,j,k) + b(i,:,k)) )
    end do
    end do
    end do
 
    return
  end subroutine dy_ab
  
  subroutine dz_ab(DzT1D, a, b, Nnode_h1D, Nnode_v, fz)
    integer, intent(in) :: nnode_h1d, nnode_v
    real(rp), intent(in) :: dzt1d(nnode_v,nnode_v)
    real(rp), intent(in) :: a(nnode_h1d,nnode_h1d,nnode_v)
    real(rp), intent(in) :: b(nnode_h1d,nnode_h1d,nnode_v)
    real(rp), intent(out) :: fz(nnode_h1d,nnode_h1d,nnode_v)
 
    integer :: i, j, k, p
    integer :: i2, j2, k2, p2
    !-------------------------------------------------
 
    do k=1, nnode_v
    do j=1, nnode_h1d    
    do i=1, nnode_h1d
      fz(i,j,k) = 0.5_rp * sum( dzt1d(:,k) * (a(i,j,k) + a(i,j,:)) * (b(i,j,k) + b(i,j,:)) )
    end do
    end do
    end do
 
    return
  end subroutine dz_ab
 
  subroutine dx_abc(DxT1D, a, b, c, Nnode_h1D, Nnode_v, fx)
    integer, intent(in) :: nnode_h1d, nnode_v
    real(rp), intent(in) :: dxt1d(nnode_h1d,nnode_h1d)
    real(rp), intent(in) :: a(nnode_h1d,nnode_h1d,nnode_v)
    real(rp), intent(in) :: b(nnode_h1d,nnode_h1d,nnode_v)
    real(rp), intent(in) :: c(nnode_h1d,nnode_h1d,nnode_v)
    real(rp), intent(out) :: fx(nnode_h1d,nnode_h1d,nnode_v)
 
    integer :: i, j, k, p
    integer :: i2, j2, k2, p2
    !-------------------------------------------------
 
    do k=1, nnode_v
    do j=1, nnode_h1d    
    do i=1, nnode_h1d
      fx(i,j,k) = 0.25_rp * sum( dxt1d(:,i) * (a(i,j,k) + a(:,j,k)) * (b(i,j,k) + b(:,j,k)) * (c(i,j,k) + c(:,j,k)) )
    end do
    end do
    end do
 
    return
  end subroutine dx_abc
 
  subroutine dy_abc(DyT1D, a, b, c, Nnode_h1D, Nnode_v, fy)
    integer, intent(in) :: nnode_h1d, nnode_v
    real(rp), intent(in) :: dyt1d(nnode_h1d,nnode_h1d)
    real(rp), intent(in) :: a(nnode_h1d,nnode_h1d,nnode_v)
    real(rp), intent(in) :: b(nnode_h1d,nnode_h1d,nnode_v)
    real(rp), intent(in) :: c(nnode_h1d,nnode_h1d,nnode_v)
    real(rp), intent(out) :: fy(nnode_h1d,nnode_h1d,nnode_v)
 
    integer :: i, j, k, p
    integer :: i2, j2, k2, p2
    !-------------------------------------------------
 
    do k=1, nnode_v
    do j=1, nnode_h1d    
    do i=1, nnode_h1d
      fy(i,j,k) = 0.25_rp * sum( dyt1d(:,j) * (a(i,j,k) + a(i,:,k)) * (b(i,j,k) + b(i,:,k)) * (c(i,j,k) + c(i,:,k)) )
    end do
    end do
    end do
 
    return
  end subroutine dy_abc
 
  subroutine dz_abc(DzT1D, a, b, c, Nnode_h1D, Nnode_v, fz)
    integer, intent(in) :: nnode_h1d, nnode_v
    real(rp), intent(in) :: dzt1d(nnode_v,nnode_v)
    real(rp), intent(in) :: a(nnode_h1d,nnode_h1d,nnode_v)
    real(rp), intent(in) :: b(nnode_h1d,nnode_h1d,nnode_v)
    real(rp), intent(in) :: c(nnode_h1d,nnode_h1d,nnode_v)
    real(rp), intent(out) :: fz(nnode_h1d,nnode_h1d,nnode_v)
 
    integer :: i, j, k, p
    integer :: i2, j2, k2, p2
    !-------------------------------------------------
 
    do k=1, nnode_v
    do j=1, nnode_h1d    
    do i=1, nnode_h1d
      fz(i,j,k) = 0.25_rp * sum( dzt1d(:,k) * (a(i,j,k) + a(i,j,:)) * (b(i,j,k) + b(i,j,:)) * (c(i,j,k) + c(i,j,:)) )
    end do
    end do
    end do
 
    return
  end subroutine dz_abc
 
end module scale_atm_dyn_dgm_nonhydro3d_rhot_hevi_splitform