scale_atm_dyn_dgm_nonhydro3d_rhot_heve_gpu.F90

Go to the documentation of this file.

!-------------------------------------------------------------------------------
!> module FElib / Fluid dyn solver / Atmosphere / Regional nonhydrostatic model / HEVE
!!
!! @par Description
!!      HEVE DGM scheme for Atmospheric dynamical process which runs on GPU. 
!!      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_nonhydro3d_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
 
  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_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_heve_gpu_init
  public :: atm_dyn_dgm_nonhydro3d_rhot_heve_gpu_final
  public :: atm_dyn_dgm_nonhydro3d_rhot_heve_cal_tend_gpu
  public :: atm_dyn_dgm_nonhydro3d_rhot_heve_cal_tend_gpu2
 
  !-----------------------------------------------------------------------------
  !
  !++ Public parameters & variables
  !
  
  !-----------------------------------------------------------------------------
  !
  !++ Private procedures & variables
  !
  !-------------------
 
contains
  !> Initialize the module for HEVE scheme on GPU
  subroutine atm_dyn_dgm_nonhydro3d_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_nonhydro3d_rhot_heve_gpu_init
 
  !> Finalize the module for HEVE scheme on GPU
  subroutine atm_dyn_dgm_nonhydro3d_rhot_heve_gpu_final()
    implicit none
    !--------------------------------------------
    
    call atm_dyn_dgm_nonhydro3d_common_final()
    return
  end subroutine atm_dyn_dgm_nonhydro3d_rhot_heve_gpu_final
 
  !-------------------------------
 
!OCL SERIAL
  subroutine atm_dyn_dgm_nonhydro3d_rhot_heve_cal_tend_gpu( &
    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_generalvc_gpu
    use scale_element_operation_gpu_driver, only: &
      elementoperationgpudriver, &
      div_kplane => elementoperationgpu_div_kplane, &
      divvar5_kplane => elementoperationgpu_divvar5_kplane, &
      divvar5_z_lift => elementoperationgpu_divvar5_z_lift, &
      vfilterpm1 => elementoperationgpu_vfilterpm1
    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) :: flux2d(elem%nnode_h1d**2,5,2,lmesh%ne)
    real(rp) :: flux2d_(elem%nnode_h1d**2,2,elem%nnode_v,lmesh%ne)   
    real(rp) :: mflxx, mflxy, mflxz
    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) :: u_, v_, w_, pt_
    real(rp) :: cor
    real(rp) :: drho(elem%np,lmesh%ne)
    real(rp) :: rdens_, gsqrth, gsqrtv, rgsqrtv, rgsqrt
    real(rp) :: gsqrt_, gsqrtdpres_, e11, e22, e33
 
    integer :: ke, ke2d
    integer :: p, ph, pz
    integer :: iv
    
    integer :: indexh2dto3d(elem%np)
    integer :: emap3dto2d(lmesh%ne)
 
    type(elementoperationgpudriver) :: element3d_operation_driver
    real(rp) :: tend_tmp(elem%nnode_h1d**2,5,elem%nnode_v,lmesh%ne)
 
    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,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_, dens_hyd, pres_hyd, therm_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)
    !$acc parallel loop gang collapse(2) present(DDENS_,MOMX_,MOMY_,MOMZ_,DRHOT_,DPRES_,DENS_hyd,THERM_hyd,lmesh%Gsqrt,lmesh%GsqrtH,lmesh%Escale)
    do ke = nes, nee
      do pz=1, nnode_v
        ke2d = lmesh%EMap3Dto2D(ke)
 
        !$acc loop vector
        do ph=1, nnode_h1d**2
          p = ph + (pz-1) * nnode_h1d**2
 
          gsqrt_ = lmesh%Gsqrt(p,ke)
          gsqrtv  = gsqrt_ / lmesh%GsqrtH(indexh2dto3d(p),ke2d)
          rgsqrtv = 1.0_rp / gsqrtv
          rgsqrt  = 1.0_rp / gsqrt_
          
          flux2d_(ph,1,pz,ke) = gsqrt_ * momx_(p,ke)
          flux2d_(ph,2,pz,ke) = gsqrt_ * momy_(p,ke)
          fluxz_store(ph,dens_vid,pz,ke) = gsqrt_ * ( &
              momz_(p,ke) * rgsqrtv        &
            + lmesh%GI3(p,ke,1) * momx_(p,ke) &
            + lmesh%GI3(p,ke,2) * momy_(p,ke) )
        end do
        call div_kplane( element3d_operation_driver, flux2d_(:,:,pz,ke), lmesh%Escale(:,ke,1,1), lmesh%Escale(:,ke,2,2), pz, &
          tend_tmp(:,dens_vid,pz,ke) )
  
        !$acc loop vector
        do ph=1, nnode_h1d**2
          p = ph + (pz-1) * nnode_h1d**2
 
          gsqrt_ = lmesh%Gsqrt(p,ke)
          rdens_  = 1.0_rp / (ddens_(p,ke) + dens_hyd(p,ke))
          pt_ = ( therm_hyd(p,ke) + drhot_(p,ke) ) * rdens_
 
          flux2d_(ph,1,pz,ke) = gsqrt_ * momx_(p,ke) * pt_
          flux2d_(ph,2,pz,ke) = gsqrt_ * momy_(p,ke) * pt_
          fluxz_store(ph,rhot_vid,pz,ke)  = fluxz_store(ph,dens_vid,pz,ke) * pt_
        end do
        call div_kplane( element3d_operation_driver, flux2d_(:,:,pz,ke), lmesh%Escale(:,ke,1,1), lmesh%Escale(:,ke,2,2), pz, &
          tend_tmp(:,rhot_vid,pz,ke) )
        
        !$acc loop vector
        do ph=1, nnode_h1d**2
          p = ph + (pz-1) * nnode_h1d**2
 
          gsqrt_ = lmesh%Gsqrt(p,ke)
          gsqrth  = lmesh%GsqrtH(indexh2dto3d(p),ke2d)
          rdens_  = 1.0_rp / (ddens_(p,ke) + dens_hyd(p,ke))
          w_ = momz_(p,ke) * rdens_
 
          flux2d_(ph,1,pz,ke) = gsqrt_ * momx_(p,ke) * w_
          flux2d_(ph,2,pz,ke) = gsqrt_ * momy_(p,ke) * w_
          fluxz_store(ph,momz_vid,pz,ke) = fluxz_store(ph,dens_vid,pz,ke) * w_ + gsqrth * dpres_(p,ke)
        end do
        call div_kplane( element3d_operation_driver, flux2d_(:,:,pz,ke), lmesh%Escale(:,ke,1,1), lmesh%Escale(:,ke,2,2), pz, &
          tend_tmp(:,momz_vid,pz,ke) )  
 
        !$acc loop vector
        do ph=1, nnode_h1d**2
          p = ph + (pz-1) * nnode_h1d**2
 
          gsqrt_ = lmesh%Gsqrt(p,ke)
          gsqrth  = lmesh%GsqrtH(indexh2dto3d(p),ke2d)
          rdens_  = 1.0_rp / (ddens_(p,ke) + dens_hyd(p,ke))
 
          u_ = momx_(p,ke) * rdens_
          gsqrtdpres_ = gsqrt_ * dpres_(p,ke)
 
          flux2d_(ph,1,pz,ke) = gsqrt_ * momx_(p,ke) * u_ + gsqrtdpres_
          flux2d_(ph,2,pz,ke) = gsqrt_ * momy_(p,ke) * u_
          fluxz_store(ph,momx_vid,pz,ke) = fluxz_store(ph,dens_vid,pz,ke) * u_ + gsqrtdpres_ * lmesh%GI3(p,ke,1)
        end do
        call div_kplane( element3d_operation_driver, flux2d_(:,:,pz,ke), lmesh%Escale(:,ke,1,1), lmesh%Escale(:,ke,2,2), pz, &
          tend_tmp(:,momx_vid,pz,ke) )  
 
 
        !$acc loop vector
        do ph=1, nnode_h1d**2
          p = ph + (pz-1) * nnode_h1d**2
 
          gsqrt_ = lmesh%Gsqrt(p,ke)
          gsqrth  = lmesh%GsqrtH(indexh2dto3d(p),ke2d)
          rdens_  = 1.0_rp / (ddens_(p,ke) + dens_hyd(p,ke))
          
          v_ = momy_(p,ke) * rdens_
          gsqrtdpres_ = gsqrt_ * dpres_(p,ke)
 
          flux2d_(ph,1,pz,ke) = gsqrt_ * momx_(p,ke) * v_ + gsqrtdpres_
          flux2d_(ph,2,pz,ke) = gsqrt_ * momy_(p,ke) * v_
          fluxz_store(ph,momy_vid,pz,ke) = fluxz_store(ph,dens_vid,pz,ke) * v_ + gsqrtdpres_ * lmesh%GI3(p,ke,2)
        end do
        call div_kplane( element3d_operation_driver, flux2d_(:,:,pz,ke), lmesh%Escale(:,ke,1,1), lmesh%Escale(:,ke,2,2), pz, &
          tend_tmp(:,momy_vid,pz,ke) )  
 
        ! !$acc loop vector
        ! do ph=1, Nnode_h1D**2
        !   p = ph + (pz-1) * Nnode_h1D**2
 
        !   Gsqrt_ = lmesh%Gsqrt(p,ke)
        !   GsqrtV  = Gsqrt_ / lmesh%GsqrtH(IndexH2Dto3D(p),ke2d)
        !   RGsqrtV = 1.0_RP / GsqrtV
        !   RGsqrt  = 1.0_RP / Gsqrt_
        !   RDENS_  = 1.0_RP / (DDENS_(p,ke) + DENS_hyd(p,ke))
 
        !   !-
        !   mflxX = Gsqrt_ * MOMX_(p,ke)
        !   mflxY = Gsqrt_ * MOMY_(p,ke)
        !   mflxZ = Gsqrt_ * ( &
        !       MOMZ_(p,ke) * RGsqrtV        &
        !     + lmesh%GI3(p,ke,1) * MOMX_(p,ke) &
        !     + lmesh%GI3(p,ke,2) * MOMY_(p,ke) )
 
        !   Flux2D(ph,DENS_VID,1,ke) = mflxX
        !   Flux2D(ph,DENS_VID,2,ke) = mflxY
        !   FluxZ_store(ph,DENS_VID,pz,ke) = mflxZ
          
        !   !-
        !   pt_ = ( THERM_hyd(p,ke) + DRHOT_(p,ke) ) * RDENS_
 
        !   Flux2D(ph,RHOT_VID,1,ke) = mflxX * pt_
        !   Flux2D(ph,RHOT_VID,2,ke) = mflxY * pt_
        !   FluxZ_store(ph,RHOT_VID,pz,ke) = mflxZ * pt_
 
        !   !-
        !   GsqrtDPRES_ = Gsqrt_ * DPRES_(p,ke)
        !   w_ = MOMZ_(p,ke) * RDENS_
        !   Flux2D(ph,MOMZ_VID,1,ke) = mflxX * w_
        !   Flux2D(ph,MOMZ_VID,2,ke) = mflxY * w_
        !   FluxZ_store(ph,MOMZ_VID,pz,ke) = mflxZ * w_ + GsqrtDPRES_ * RGsqrtV
 
        !   !-
        !   u_ = MOMX_(p,ke) * RDENS_
 
        !   Flux2D(ph,MOMX_VID,1,ke) = mflxX * u_ + GsqrtDPRES_
        !   Flux2D(ph,MOMX_VID,2,ke) = mflxY * u_ 
        !   FluxZ_store(ph,MOMX_VID,pz,ke) = mflxZ * u_ + GsqrtDPRES_ * lmesh%GI3(p,ke,1)
 
        !   v_ = MOMY_(p,ke) * RDENS_
        !   Flux2D(ph,MOMY_VID,1,ke) = mflxX * v_
        !   Flux2D(ph,MOMY_VID,2,ke) = mflxY * v_ + GsqrtDPRES_ 
        !   FluxZ_store(ph,MOMY_VID,pz,ke) = mflxZ * v_ + GsqrtDPRES_ * lmesh%GI3(p,ke,2)
        ! end do
        ! call DivVar5_kplane( element3D_operation_driver,              &
        !   Flux2D(:,:,:,ke), lmesh%Escale(:,ke,1,1), lmesh%Escale(:,ke,2,2), pz, &
        !   tend_tmp(:,:,pz,ke) )
      end do
    end do
    !$acc parallel loop gang
    do ke = nes, nee
      call divvar5_z_lift( element3d_operation_driver, &
        fluxz_store(:,:,:,ke),  del_flux(:,:,ke), lmesh%Escale(:,ke,3,3), lmesh%Gsqrt(:,ke), &
        tend_tmp(:,:,:,ke) )
    end do
 
    call vfilterpm1( element3d_operation_driver, ddens_, lmesh%NeA, lmesh%Ne, &
      drho )
 
    !$acc parallel loop gang collapse(2)
    do ke = nes, nee
      !--
      do pz=1, nnode_v
         ke2d = lmesh%EMap3Dto2D(ke)
 
        !$acc loop vector
        do ph=1, nnode_h1d**2
          p = ph + (pz-1) * nnode_h1d**2
 
          dens_dt(p,ke) = - tend_tmp(ph,dens_vid,pz,ke)
          rhot_dt(p,ke) = - tend_tmp(ph,rhot_vid,pz,ke)
 
 
          momz_dt(p,ke) = - tend_tmp(ph,momz_vid,pz,ke) &
                          - grav * drho(p,ke)
          
          cor = coriolis(ph,ke2d)
          momx_dt(p,ke) = - tend_tmp(ph,momx_vid,pz,ke) &
                          - dphyddx(p,ke) &
                          + cor * momy_(p,ke)
          momy_dt(p,ke) = - tend_tmp(ph,momy_vid,pz,ke) &
                          - dphyddy(p,ke) &
                          - cor * momx_(p,ke)
        end do
      end do
    end do
    call prof_rapend('cal_dyn_tend_interior', 3)
    !$acc end data
 
    call element3d_operation_driver%Final()
    return
  end subroutine atm_dyn_dgm_nonhydro3d_rhot_heve_cal_tend_gpu
 
!OCL SERIAL
  subroutine atm_dyn_dgm_nonhydro3d_rhot_heve_cal_tend_gpu2( &
    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_generalvc_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_, dens_hyd, pres_hyd, therm_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 )
    !$acc wait(1)                                                      ! (in)
    call prof_rapend('cal_dyn_tend_bndflux', 3)
 
    call prof_rapstart('cal_dyn_tend_interior', 3)
    call atm_dyn_dgm_nonhydro3d_rhot_heve_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, coriolis, dphyddx, dphyddy,              & ! (in)
      lmesh%Gsqrt,lmesh%GsqrtH,lmesh%GI3(:,:,1),lmesh%GI3(:,:,2),          &
      lmesh%Escale(:,:,1,1), lmesh%Escale(:,:,2,2), lmesh%Escale(:,:,3,3), & ! (in)
      flux2d, fluxz_store, tend_tmp, drho, & ! (in)
      element3d_operation_driver, & ! (in)
      lmesh%EMap3Dto2D, lmesh,elem, lmesh%NeS, lmesh%NeE, lmesh%NeA, lmesh%Ne2DA, elem%Nnode_h1D, elem%Nnode_v )                         ! (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_nonhydro3d_rhot_heve_cal_tend_gpu2
 
  subroutine atm_dyn_dgm_nonhydro3d_rhot_heve_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, coriolis, dphyddx, dphyddy,              & ! (in)
    gsqrt,gsqrth,g13,g23, e11,e22,e33,                            & ! (in)
    flux2d, fluxz_store, tend_tmp, drho,                          & ! (in)
    element3d_operation_driver,                                   & ! (in)
    emap3dto2d, lmesh,elem, nes, nee, nea, ne2da, nnode_h1d, nnode_v )                         ! (in)
    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)  :: coriolis(nnode_h1d**2,ne2da)
    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(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 :: ke, ke2d, ph,pz
 
    real(rp) :: mflxx, mflxy, mflxz
    real(rp) :: u_, v_, w_, pt_
    real(rp) :: cor
    real(rp) :: rdens_, gsqrtv, rgsqrtv, rgsqrt
    real(rp) :: gsqrt_, gsqrtdpres_
 
    !----------------------------
    !!$acc parallel present(DENS_dt,MOMX_dt,MOMY_dt,MOMZ_dt,RHOT_dt,&
    !!$acc DDENS_,MOMX_,MOMY_,MOMZ_,DRHOT_,DPRES_,DENS_hyd,THERM_hyd,&
    !!$acc DPhydDx,DPhydDy,Coriolis, del_flux, tend_tmp,drho, Gsqrt,GsqrtH,G13,G23,E11,E22,E33,Flux2D,FluxZ_store, EMap3Dto2D)
 
    !$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))
 
        !-
        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_ + gsqrtdpres_
        flux2d(ph,momx_vid,pz,ke,2) = mflxy * u_ 
        fluxz_store(ph,momx_vid,pz,ke) = mflxz * u_ + gsqrtdpres_ * g13(ph,pz,ke)
 
        v_ = momy_(ph,pz,ke) * rdens_
        flux2d(ph,momy_vid,pz,ke,1) = mflxx * v_
        flux2d(ph,momy_vid,pz,ke,2) = mflxy * v_ + gsqrtdpres_ 
        fluxz_store(ph,momy_vid,pz,ke) = mflxz * v_ + gsqrtdpres_ * 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,Coriolis, 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)
        
        cor = coriolis(ph,ke2d)
        momx_dt(ph,pz,ke) = - tend_tmp(ph,momx_vid,pz,ke) &
                        - dphyddx(ph,pz,ke) &
                        + cor * momy_(ph,pz,ke)
        momy_dt(ph,pz,ke) = - tend_tmp(ph,momy_vid,pz,ke) &
                        - dphyddy(ph,pz,ke) &
                        - cor * momx_(ph,pz,ke)
      end do
    end do
    end do
    !$acc end parallel
    return
  end subroutine atm_dyn_dgm_nonhydro3d_rhot_heve_cal_tend_interior_gpu
 
end module scale_atm_dyn_dgm_nonhydro3d_rhot_heve_gpu