scale_atm_dyn_dgm_globalnonhydro3d_rhot_heve_gpu.F90

Go to the documentation of this file.

!-------------------------------------------------------------------------------
!> module FElib / Fluid dyn solver / Atmosphere / Global nonhydrostatic model / HEVE
!!
!! @par Description
!!      HEVE 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, Xuanzhengbo Ren, and Team SCALE
!<
!-------------------------------------------------------------------------------
#include "scaleFElib.h"
module scale_atm_dyn_dgm_globalnonhydro3d_rhot_heve_gpu
  !-----------------------------------------------------------------------------
  !
  !++ 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, &
    rplanet => const_radius
 
  use scale_sparsemat  
  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_base2d, only: meshbase2d
  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_globalnonhydro3d_rhot_heve_gpu_init
  public :: atm_dyn_dgm_globalnonhydro3d_rhot_heve_gpu_final
  public :: atm_dyn_dgm_globalnonhydro3d_rhot_heve_gpu_cal_tend_shallow_atm
  public :: atm_dyn_dgm_globalnonhydro3d_rhot_heve_gpu_cal_tend_deep_atm
 
  !-----------------------------------------------------------------------------
  !
  !++ Public parameters & variables
  !
  
  !-----------------------------------------------------------------------------
  !
  !++ Private procedures & variables
  !
  !-------------------
 
contains
!OCL SERIAL
  subroutine atm_dyn_dgm_globalnonhydro3d_rhot_heve_gpu_init( mesh )
    implicit none
    class(meshbase3d), intent(in) :: mesh
    !--------------------------------------------
 
    call atm_dyn_dgm_nonhydro3d_common_init( mesh )
 
    return
  end subroutine atm_dyn_dgm_globalnonhydro3d_rhot_heve_gpu_init
 
!OCL SERIAL
  subroutine atm_dyn_dgm_globalnonhydro3d_rhot_heve_gpu_final()
    implicit none
    !--------------------------------------------
    
    call atm_dyn_dgm_nonhydro3d_common_final()    
    return
  end subroutine atm_dyn_dgm_globalnonhydro3d_rhot_heve_gpu_final
 
  !-------------------------------
 
  !OCL SERIAL
  subroutine atm_dyn_dgm_globalnonhydro3d_rhot_heve_gpu_cal_tend_shallow_atm( &
    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_heve_numflux_gpu, only: &
      get_ebnd_flux => atm_dyn_dgm_nonhydro3d_rhot_heve_numflux_get_generalhvc_gpu
    use scale_element_operation_gpu_driver, only: &
      elementoperationgpudriver
    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)
 
    type(elementoperationgpudriver) :: element3d_operation_driver
    real(rp) :: tend_tmp(elem%nnode_h1d**2,5,elem%nnode_v,lmesh%ne)
    real(rp) :: flux2d(elem%nnode_h1d**2,5,elem%nnode_v,lmesh%ne,2)
    real(rp) :: fluxz_store(elem%nnode_h1d**2,5,elem%nnode_v,lmesh%ne)
    real(rp) :: del_flux(elem%nfptot,prgvar_num,lmesh%ne)
    real(rp) :: drho(elem%np,lmesh%ne)
 
    integer :: indexh2dto3d(elem%np)
 
    integer :: nes, nee, nnode_h1d, nnode_v
    !------------------------------------------------------------------------
 
    indexh2dto3d(:) = elem%IndexH2Dto3D(:)
    call element3d_operation_driver%Init( element3d_operation )
    nes = lmesh%NeS
    nee = lmesh%NeE
    nnode_h1d = elem%Nnode_h1D
    nnode_v = elem%Nnode_v
 
    !$acc data present(DDENS_,MOMX_,MOMY_,MOMZ_,DRHOT_,DPRES_,    &
    !$acc              DENS_hyd,PRES_hyd,THERM_hyd,               &
    !$acc              CORIOLIS,Rtot,CVtot,CPtot,DPhydDx,DPhydDy, &
    !$acc              DENS_dt,MOMX_dt,MOMY_dt,MOMZ_dt,RHOT_dt, lmesh,elem) &
    !$acc      create(del_flux,Flux2D,FluxZ_store,tend_tmp,drho) copyin(IndexH2Dto3D)
 
    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%GIJ(:,:,1,1), lmesh%GIJ(:,:,1,2), lmesh%GIJ(:,:,2,2), & ! (in)
      lmesh%GsqrtH, lmesh%gam, 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)
    !$acc wait(1)
    call prof_rapend('cal_dyn_tend_bndflux', 3)
 
    !-----
    call prof_rapstart('cal_dyn_tend_interior', 3)
    call cal_tend_interior_gpu( &
      dens_dt, momx_dt, momy_dt, momz_dt, rhot_dt,                         & ! (out)
      ddens_, momx_, momy_, momz_, drhot_, dpres_, del_flux,               & ! (in)
      dens_hyd, therm_hyd, dphyddx, dphyddy,                               & ! (in)
      lmesh%Gsqrt,lmesh%GsqrtH,lmesh%GI3(:,:,1),lmesh%GI3(:,:,2),          & ! (in)
      lmesh%Escale(:,:,1,1), lmesh%Escale(:,:,2,2), lmesh%Escale(:,:,3,3), & ! (in)
      lmesh2d%pos_en(:,:,1), lmesh2d%pos_en(:,:,2),                        & ! (in)
      flux2d, fluxz_store, tend_tmp, drho,                                 & ! (in)
      element3d_operation_driver, lmesh%EMap3Dto2D, indexh2dto3d,          & ! (in)
      lmesh,elem, lmesh%NeS, lmesh%NeE, lmesh%NeA, lmesh%Ne2DA, elem%Nnode_h1D, elem%Nnode_v ) ! (in)                     ! (in)    
    !$acc wait(1)
    call prof_rapend('cal_dyn_tend_interior', 3)
    
    !$acc end data
 
    call element3d_operation_driver%Final()
    return
  end subroutine atm_dyn_dgm_globalnonhydro3d_rhot_heve_gpu_cal_tend_shallow_atm
 
  subroutine cal_tend_interior_gpu( &
    DENS_dt, MOMX_dt, MOMY_dt, MOMZ_dt, RHOT_dt,                  & ! (out)
    ddens_, momx_, momy_, momz_, drhot_, dpres_, del_flux,        & ! (in)
    dens_hyd, therm_hyd, dphyddx, dphyddy,                        & ! (in)
    gsqrt,gsqrth,g13,g23, e11,e22,e33, alph2d, beta2d,            & ! (in)
    flux2d, fluxz_store, tend_tmp, drho,                          & ! (in)
    element3d_operation_driver,                                   & ! (in)
    emap3dto2d, indexh2dto3d, lmesh,elem, nes, nee, nea, ne2da, nnode_h1d, nnode_v ) ! (in)
    use scale_const, only: &
      ohm => const_ohm
    use scale_element_operation_gpu_driver, only: &
      elementoperationgpudriver, &
      div_kplane => elementoperationgpu_div_kplane, &
      divvar5_kplane => elementoperationgpu_divvar5_kplane2, &
      divvar5_z_lift => elementoperationgpu_divvar5_z_lift, &
      divvar5 => elementoperationgpu_divvar5, &
      vfilterpm1 => elementoperationgpu_vfilterpm1
    implicit none
    class(localmesh3d), intent(in) :: lmesh
    class(elementbase3d), intent(in) :: elem
    integer, intent(in) :: nes, nee, nea, ne2da
    integer, intent(in) :: nnode_h1d, nnode_v
    real(rp), intent(out) :: dens_dt(nnode_h1d**2,nnode_v,nea)
    real(rp), intent(out) :: momx_dt(nnode_h1d**2,nnode_v,nea)
    real(rp), intent(out) :: momy_dt(nnode_h1d**2,nnode_v,nea)
    real(rp), intent(out) :: momz_dt(nnode_h1d**2,nnode_v,nea)
    real(rp), intent(out) :: rhot_dt(nnode_h1d**2,nnode_v,nea)
    real(rp), intent(in)  :: ddens_(nnode_h1d**2,nnode_v,nea)
    real(rp), intent(in)  :: momx_(nnode_h1d**2,nnode_v,nea)
    real(rp), intent(in)  :: momy_(nnode_h1d**2,nnode_v,nea)
    real(rp), intent(in)  :: momz_(nnode_h1d**2,nnode_v,nea)
    real(rp), intent(in)  :: drhot_(nnode_h1d**2,nnode_v,nea)
    real(rp), intent(in)  :: dpres_(nnode_h1d**2,nnode_v,nea)
    real(rp), intent(in)  :: del_flux(elem%nfptot,prgvar_num,lmesh%ne)
    real(rp), intent(in)  :: dens_hyd(nnode_h1d**2,nnode_v,nea)
    real(rp), intent(in)  :: therm_hyd(nnode_h1d**2,nnode_v,nea)
    real(rp), intent(in)  :: dphyddx(nnode_h1d**2,nnode_v,nea)
    real(rp), intent(in)  :: dphyddy(nnode_h1d**2,nnode_v,nea)
    real(rp), intent(in)  :: gsqrt(nnode_h1d**2,nnode_v,nea)
    real(rp), intent(in)  :: gsqrth(nnode_h1d**2,lmesh%ne2d)
    real(rp), intent(in) :: g13(nnode_h1d**2,nnode_v,nea)
    real(rp), intent(in) :: g23(nnode_h1d**2,nnode_v,nea)
    real(rp), intent(in)  :: e11(nnode_h1d**2,nnode_v,lmesh%ne)
    real(rp), intent(in)  :: e22(nnode_h1d**2,nnode_v,lmesh%ne)
    real(rp), intent(in)  :: e33(nnode_h1d**2,nnode_v,lmesh%ne)
    real(rp), intent(in)  :: alph2d(nnode_h1d**2,lmesh%ne2d)
    real(rp), intent(in)  :: beta2d(nnode_h1d**2,lmesh%ne2d)
    real(rp), intent(out) :: flux2d(elem%nnode_h1d**2,5,elem%nnode_v,lmesh%ne,2)
    real(rp), intent(out) :: fluxz_store(elem%nnode_h1d**2,5,elem%nnode_v,lmesh%ne)
    real(rp), intent(out) :: tend_tmp(elem%nnode_h1d**2,5,elem%nnode_v,lmesh%ne)
    real(rp), intent(out) :: drho(nnode_h1d**2,nnode_v,lmesh%ne)
    type(elementoperationgpudriver), intent(in) :: element3d_operation_driver
    integer, intent(in) :: emap3dto2d(lmesh%ne)
    integer, intent(in) :: indexh2dto3d(elem%nnode_h1d**2,nnode_v)
 
    integer :: ke, ke2d
    integer :: ph, pz
 
    real(rp) :: mflxx, mflxy, mflxz
    real(rp) :: u_, v_, w_, pt_
    real(rp) :: cor_x, cor_y
    real(rp) :: rdens_, gsqrtv, rgsqrtv, rgsqrt
    real(rp) :: gsqrt_, gsqrtdpres_
 
    real(rp) :: g11, g12, g22
    real(rp) :: x, y
    real(rp) :: twoovdel2
 
    real(rp) :: s
    logical :: is_panel1to4
    !----------------------------
 
    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
 
    !$acc parallel loop gang collapse(2) &
    !$acc   present(DDENS_,MOMX_,MOMY_,MOMZ_,DRHOT_,DPRES_,DENS_hyd,THERM_hyd,Flux2D,FluxZ_store, EMap3Dto2D) async(1)
    do ke = nes, nee
    do pz=1, nnode_v
      ke2d = emap3dto2d(ke)
 
      !$acc loop vector
      do ph=1, nnode_h1d**2
        gsqrt_ = gsqrt(ph,pz,ke)
        gsqrtv  = gsqrt_ / gsqrth(ph,ke2d)
        rgsqrtv = 1.0_rp / gsqrtv
        rgsqrt  = 1.0_rp / gsqrt_
        rdens_  = 1.0_rp / ( ddens_(ph,pz,ke) + dens_hyd(ph,pz,ke) )
 
        g11 = lmesh%GIJ(ph,ke2d,1,1)
        g12 = lmesh%GIJ(ph,ke2d,1,2)
        g22 = lmesh%GIJ(ph,ke2d,2,2)
 
        !-
        mflxx = gsqrt_ * momx_(ph,pz,ke)
        mflxy = gsqrt_ * momy_(ph,pz,ke)
        mflxz = gsqrt_ * ( &
            momz_(ph,pz,ke) * rgsqrtv        &
          + g13(ph,pz,ke) * momx_(ph,pz,ke) &
          + g23(ph,pz,ke) * momy_(ph,pz,ke) )
 
        flux2d(ph,dens_vid,pz,ke,1) = mflxx
        flux2d(ph,dens_vid,pz,ke,2) = mflxy
        fluxz_store(ph,dens_vid,pz,ke) = mflxz
        
        !-
        pt_ = ( therm_hyd(ph,pz,ke) + drhot_(ph,pz,ke) ) * rdens_
 
        flux2d(ph,rhot_vid,pz,ke,1) = mflxx * pt_
        flux2d(ph,rhot_vid,pz,ke,2) = mflxy * pt_
        fluxz_store(ph,rhot_vid,pz,ke) = mflxz * pt_
 
        !-
        gsqrtdpres_ = gsqrt_ * dpres_(ph,pz,ke)
        w_ = momz_(ph,pz,ke) * rdens_
        flux2d(ph,momz_vid,pz,ke,1) = mflxx * w_
        flux2d(ph,momz_vid,pz,ke,2) = mflxy * w_
        fluxz_store(ph,momz_vid,pz,ke) = mflxz * w_ + gsqrtdpres_ * rgsqrtv
 
        !-        
        u_ = momx_(ph,pz,ke) * rdens_
        flux2d(ph,momx_vid,pz,ke,1) = mflxx * u_ + g11 * gsqrtdpres_
        flux2d(ph,momx_vid,pz,ke,2) = mflxy * u_ + g12 * gsqrtdpres_
        fluxz_store(ph,momx_vid,pz,ke) = mflxz * u_ + gsqrtdpres_ * ( g11 * g13(ph,pz,ke) + g12 * g23(ph,pz,ke) )
 
        v_ = momy_(ph,pz,ke) * rdens_
        flux2d(ph,momy_vid,pz,ke,1) = mflxx * v_ + g12 * gsqrtdpres_
        flux2d(ph,momy_vid,pz,ke,2) = mflxy * v_ + g22 * gsqrtdpres_
        fluxz_store(ph,momy_vid,pz,ke) = mflxz * v_ + gsqrtdpres_ * ( g12 * g13(ph,pz,ke) + g22 * g23(ph,pz,ke) )
      end do
    end do
    end do
 
    call divvar5( element3d_operation_driver, flux2d(:,:,:,:,1), flux2d(:,:,:,:,2), fluxz_store, del_flux, &
      e11, e22, e33, gsqrt, lmesh%Ne, tend_tmp )
    
    call vfilterpm1( element3d_operation_driver, ddens_, lmesh%NeA, lmesh%Ne, &
      drho )
 
 
    !$acc parallel loop gang collapse(2) &
    !$acc   present(DENS_dt,MOMX_dt,MOMY_dt,MOMZ_dt,RHOT_dt, DPhydDx,DPhydDy, tend_tmp, drho, EMap3Dto2D) async(1)
    do ke = nes, nee
    do pz=1, nnode_v
      ke2d = emap3dto2d(ke)
      !$acc loop vector
      do ph=1, nnode_h1d**2
        dens_dt(ph,pz,ke) = - tend_tmp(ph,dens_vid,pz,ke)
        rhot_dt(ph,pz,ke) = - tend_tmp(ph,rhot_vid,pz,ke)
 
 
        momz_dt(ph,pz,ke) = - tend_tmp(ph,momz_vid,pz,ke) &
                            - grav * drho(ph,pz,ke)
      end do
      !$acc loop vector
      do ph=1, nnode_h1d**2
        !-        
        x = tan(alph2d(ph,ke2d))
        y = tan(beta2d(ph,ke2d))
        twoovdel2 = 2.0_rp / ( 1.0_rp + x**2 + y**2 )
 
        cor_x = s * ohm * twoovdel2 * ( - x * y             * momx_(ph,pz,ke) + ( 1.0_rp + y**2 ) * momy_(ph,pz,ke) )
        cor_y = s * ohm * twoovdel2 * ( - ( 1.0_rp + x**2 ) * momx_(ph,pz,ke) + x * y             * momy_(ph,pz,ke) )
        if ( is_panel1to4 ) then
          cor_x = s * y * cor_x
          cor_y = s * y * cor_y
        end if
 
        g11 = lmesh%GIJ(ph,ke2d,1,1)
        g12 = lmesh%GIJ(ph,ke2d,1,2)
        g22 = lmesh%GIJ(ph,ke2d,2,2)
 
        rdens_  = 1.0_rp / ( ddens_(ph,pz,ke) + dens_hyd(ph,pz,ke) )
        u_ = momx_(ph,pz,ke) * rdens_
        v_ = momy_(ph,pz,ke) * rdens_
 
 
        momx_dt(ph,pz,ke) = - tend_tmp(ph,momx_vid,pz,ke) &
                          - ( g11 * dphyddx(ph,pz,ke) + g12 * dphyddy(ph,pz,ke) )     &
                          - twoovdel2 * y *                                           &
                            ( x * y * u_ + ( 1.0_rp + y**2 ) * v_ ) * momx_(ph,pz,ke) &
                          + cor_x
        momy_dt(ph,pz,ke) = - tend_tmp(ph,momy_vid,pz,ke) &
                          - ( g12 * dphyddx(ph,pz,ke) + g22 * dphyddy(ph,pz,ke) )       &
                          - twoovdel2 * y *                                             &
                            ( - ( 1.0_rp + x**2 ) * u_ + x * y * v_ ) * momy_(ph,pz,ke) &
                          - cor_y
      end do
    end do
    end do
    !$acc end parallel
    return
  end subroutine cal_tend_interior_gpu
 
!OCL SERIAL
  subroutine atm_dyn_dgm_globalnonhydro3d_rhot_heve_gpu_cal_tend_deep_atm( &
    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_heve_numflux, only: &
      get_ebnd_flux => atm_dyn_dgm_nonhydro3d_rhot_heve_numflux_get_generalhvc_asis
    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) :: 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), drho(elem%np)
 
    real(rp) :: g11(elem%np), g12(elem%np), g22(elem%np)
    real(rp) :: gsqrtv(elem%np), rgsqrtv(elem%np), rgam2(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) :: om1(elem%np), om2(elem%np), om3(elem%np), del(elem%np), r(elem%np)
    logical :: is_panel1to4
    real(rp) :: s
 
    integer :: ke, ke2d
    integer :: p
    
    real(rp) :: rgamm    
    real(rp) :: rp0
    real(rp) :: 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%gam, 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)
    rgamm = cvdry / cpdry
    rp0   = 1.0_rp / pres00
    p0ovr = pres00 / rdry
 
    s = 2.0_rp * ohm
    is_panel1to4 = .true.
    if ( lmesh%panelID == 5 ) then
      is_panel1to4 = .false.
    else if ( lmesh%panelID == 6 ) then
      is_panel1to4 = .false.
      s = - s
    end if
 
    !$omp parallel private(                           &
    !$omp RHOT_, rdens_, u_, v_, w_, wt_,             &
    !$omp Fx, Fy, Fz, LiftDelFlx,                     &
    !$omp drho, DPRES_hyd, GradPhyd_x, GradPhyd_y,    &
    !$omp G11, G12, G22, Rgam2, GsqrtV, RGsqrtV,      &
    !$omp X, Y, twoOVdel2,                            &
    !$omp OM1, OM2, OM3, DEL, R, 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)
      del(:) = sqrt( 1.0_rp + x(:)**2 + y(:)**2 )
      twoovdel2(:) = 2.0_rp / ( 1.0_rp + x(:)**2 + y(:)**2 )
 
      r(:) = rplanet * lmesh%gam(:,ke)
 
      !-  pnl=1~4: OM1:                     0, OM2: s del / (r (1+Y^2)) ,   s OM3 : Y /del
      !-  pnl=5,6: OM1: - s X del/(r (1+X^2)), OM2: -  s Y del/(r(1+Y^2)), OM3 : s/del
      if ( is_panel1to4 ) then
        om1(:) = 0.0_rp
        om2(:) = s * del(:) / ( r(:) * ( 1.0_rp + y(:)**2 ) )
        om3(:) = s * y(:) / del(:)
      else
        om1(:) = - s * x(:) * del(:) / ( r(:) * ( 1.0_rp + x(:)**2 ) )
        om2(:) = - s * y(:) * del(:) / ( r(:) * ( 1.0_rp + y(:)**2 ) )
        om3(:) =   s / del(:)
      end if
 
      drho(:) = matmul(intrpmat_vpordm1, ddens_(:,ke))
 
      !-- 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(dz, lmesh%Gsqrt(:,ke) * ( ddens_(:,ke) + dens_hyd(:,ke) ) * wt_(:), fz)
      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(:)     &
          + lmesh%Escale(:,ke,3,3) * fz(:)     &
          + 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(:) *                                                           & !-> metric terms
            ( x(:) * y(:) * u_(:) - ( 1.0_rp + y(:)**2 ) * v_(:) ) * momx_(:,ke)            & !
          - 2.0_rp * u_(:) * momz_(:,ke) / r(:)                                             & !<-
          - ( g11(:) * gradphyd_x(:) + g12(:) * gradphyd_y(:) ) * rgsqrtv(:)                & !-> gradient hydrostatic pressure
          - lmesh%Gsqrt(:,ke) * (  g11(:) * ( om2(:) * momz_(:,ke) - om3(:) * momy_(:,ke) ) & !-> Coriolis term
                                 - g12(:) * ( om1(:) * momz_(:,ke) - om3(:) * momx_(:,ke) ) )
 
      !-- 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(:) *                                                           & !-> metric terms
              ( - (1.0_rp + x(:)**2) * u_(:) + x(:) * y(:) * v_(:) ) * momy_(:,ke)            & !
            - 2.0_rp * v_(:) * momz_(:,ke) / r(:)                                             & !<-
            - ( g12(:) * gradphyd_x(:) + g22(:) * gradphyd_y(:) ) * rgsqrtv(:)                & !-> gradient hydrostatic pressure
            - lmesh%Gsqrt(:,ke) * (  g12(:) * ( om2(:) * momz_(:,ke) - om3(:) * momy_(:,ke) ) & !-> Coriolis term
                                   - g22(:) * ( om1(:) * momz_(:,ke) - om3(:) * 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) + rgsqrtv(:) * dpres_(:,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)                                                               &
          - 0.25_rp * r(:) * twoovdel2(:)**2 * ( 1.0_rp * x(:)**2 ) * ( 1.0_rp * y(:)**2 )                      & !-> metric terms
            * ( - ( 1.0_rp + x(:)**2 ) * momx_(:,ke) * u_(:)                                                    & !
                + 2.0_rp * x(:) * y(:) * momx_(:,ke) * v_(:)                                                    & !
                - ( 1.0_rp + y(:)**2 ) * momy_(:,ke) * v_(:)  )                                                 & !<-
          + 2.0_rp * dpres_(:,ke) / r(:)                                                                        & !-> metric term with gradient of pressure deviaition
          - lmesh%Gsqrt(:,ke) * ( om1(:) * momy_(:,ke) - om2(:) * momx_(:,ke) )                                 & !-> Coriolis term
          - grav * rgam2(:) * drho(:)                                                                             !-> buoyancy term
 
      !-- RHOT
      call sparsemat_matmul(dx, lmesh%Gsqrt(:,ke) * u_(:) * rhot_(:), fx)
      call sparsemat_matmul(dy, lmesh%Gsqrt(:,ke) * v_(:) * rhot_(:), fy)
      call sparsemat_matmul(dz, lmesh%Gsqrt(:,ke) * wt_(:) * rhot_(:), fz)
      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(:)      &
            + lmesh%Escale(:,ke,3,3) * fz(:)      &
            + liftdelflx(:) ) / lmesh%Gsqrt(:,ke) 
    end do
    !$omp end do
    !$omp end parallel
    call prof_rapend('cal_dyn_tend_interior', 3)
 
    return
  end subroutine atm_dyn_dgm_globalnonhydro3d_rhot_heve_gpu_cal_tend_deep_atm
 
end module scale_atm_dyn_dgm_globalnonhydro3d_rhot_heve_gpu