scale_atm_dyn_dgm_nonhydro3d_common.F90

Go to the documentation of this file.

!-------------------------------------------------------------------------------
!> module FElib / Fluid dyn solver / Atmosphere / Nonhydrostatic model / Common
!!
!! @par Description
!!      A common model for atmospheric nonhydrostatic dynamical core 
!!
!! @author Yuta Kawai, Team SCALE
!<
!-------------------------------------------------------------------------------
#include "scaleFElib.h"
module scale_atm_dyn_dgm_nonhydro3d_common
  !-----------------------------------------------------------------------------
  !
  !++ 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_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_model_var_manager, only: &
    modelvarmanager, variableinfo
 
  use scale_element_operation_base, only: elementoperationbase3d
 
  !-----------------------------------------------------------------------------
  implicit none
  private
  !-----------------------------------------------------------------------------
  !
  !++ Public procedures
  !
  public :: atm_dyn_dgm_nonhydro3d_common_init
  public :: atm_dyn_dgm_nonhydro3d_common_final
  public :: atm_dyn_dgm_nonhydro3d_common_setup_variables
  public :: atm_dyn_dgm_nonhydro3d_common_get_varinfo
  public :: atm_dyn_dgm_nonhydro3d_common_calc_pressure
  public :: atm_dyn_dgm_nonhydro3d_common_calc_therm_phyd
  public :: atm_dyn_dgm_nonhydro3d_common_drhot2pres
  public :: atm_dyn_dgm_nonhydro3d_common_entot2pres
  public :: atm_dyn_dgm_nonhydro3d_common_drhot2entot
  public :: atm_dyn_dgm_nonhydro3d_common_calc_rhot_hyd
  public :: atm_dyn_dgm_nonhydro3d_common_calc_phyd_hgrad_lc
 
  !-----------------------------------------------------------------------------
  !
  !++ Public parameters & variables
  !
  
  !-
  integer, public, parameter :: prgvar_ddens_id   = 1
  integer, public, parameter :: prgvar_therm_id   = 2 ! Variable associated with energy equation (DRHOT or ETOT)
  integer, public, parameter :: prgvar_drhot_id   = 2
  integer, public, parameter :: prgvar_etot_id    = 2
  integer, public, parameter :: prgvar_momz_id    = 3
  integer, public, parameter :: prgvar_momx_id    = 4
  integer, public, parameter :: prgvar_momy_id    = 5
  integer, public, parameter :: prgvar_scalar_num = 3
  integer, public, parameter :: prgvar_hvec_num   = 1
  integer, public, parameter :: prgvar_num        = 5
 
  integer, public, parameter :: phytend_dens_id     = 1
  integer, public, parameter :: phytend_momx_id     = 2
  integer, public, parameter :: phytend_momy_id     = 3
  integer, public, parameter :: phytend_momz_id     = 4
  integer, public, parameter :: phytend_rhot_id     = 5
  integer, public, parameter :: phytend_rhoh_id     = 6
  integer, public, parameter :: phytend_num         = 6
 
  !-
  integer, public, parameter :: auxvar_preshydro_id     = 1
  integer, public, parameter :: auxvar_denshydro_id     = 2
  integer, public, parameter :: auxvar_thermhydro_id    = 3
  integer, public, parameter :: auxvar_pres_id          = 4
  integer, public, parameter :: auxvar_pt_id            = 5
  integer, public, parameter :: auxvar_rtot_id          = 6
  integer, public, parameter :: auxvar_cvtot_id         = 7
  integer, public, parameter :: auxvar_cptot_id         = 8
  integer, public, parameter :: auxvar_qdry_id          = 9
  integer, public, parameter :: auxvar_preshydro_ref_id = 10
  integer, public, parameter :: auxvar_num              = 10
 
  
  
  real(rp), public, allocatable :: intrpmat_vpordm1(:,:)
 
  !-----------------------------------------------------------------------------
  !
  !++ Private procedures & variables
  !
  !-------------------
  
contains
  !> Initialize a common module for atmospheric nonhydrostatic dynamical core
  subroutine atm_dyn_dgm_nonhydro3d_common_init( mesh )
    implicit none
    class(meshbase3d), intent(in) :: mesh
 
    type(elementbase3d), pointer :: elem
    !--------------------------------------------
 
    elem => mesh%refElem3D
    allocate( intrpmat_vpordm1(elem%Np,elem%Np) )
    call mesh%refElem3D%Generate_ModalTruncationMat( elem%PolyOrder_h, elem%PolyOrder_v-1, & ! (in)
      intrpmat_vpordm1 ) ! (out)
    return
  end subroutine atm_dyn_dgm_nonhydro3d_common_init
 
 
  !> Finalize a common module for atmospheric nonhydrostatic dynamical core
  subroutine atm_dyn_dgm_nonhydro3d_common_final()
    implicit none
    !--------------------------------------------
    
    deallocate( intrpmat_vpordm1 )
    return
  end subroutine atm_dyn_dgm_nonhydro3d_common_final  
 
  !> Get variable information for atmospheric nonhydrostatic dynamical core
  subroutine atm_dyn_dgm_nonhydro3d_common_get_varinfo( &
    prgvar_info, auxvar_info, phytend_info              )
 
    implicit none
    type(variableinfo), intent(out) :: prgvar_info(prgvar_num)
    type(variableinfo), intent(out) :: auxvar_info(auxvar_num)
    type(variableinfo), intent(out), optional :: phytend_info(phytend_num)
 
    type(variableinfo) :: prgvar_varinfo(prgvar_num)
    DATA prgvar_varinfo / &
      variableinfo( prgvar_ddens_id, 'DDENS', 'deviation of density',         &
                    'kg/m3',  3, 'XYZ',  'air_density'                   ),   &
      variableinfo( prgvar_therm_id, 'THERM', 'THERM',                        &
                    '-', 3, 'XYZ',  ''                                     ), &
      variableinfo( prgvar_momz_id , 'MOMZ', 'momentum z',                    &
                    'kg/m2/s', 3, 'XYZ', 'northward_mass_flux_of_air'      ), &
      variableinfo( prgvar_momx_id , 'MOMX', 'momentum x',                    &
                    'kg/m2/s', 3, 'XYZ', 'upward_mass_flux_of_air'         ), &
      variableinfo( prgvar_momy_id , 'MOMY', 'momentum y',                    &
                    'kg/m2/s', 3, 'XYZ', 'eastward_mass_flux_of_air'       )  /
 
    type(variableinfo) :: auxvar_varinfo(auxvar_num)
    DATA auxvar_varinfo / &
      variableinfo( auxvar_preshydro_id, 'PRES_hyd', 'hydrostatic part of pressure',             &
                      'Pa', 3, 'XYZ', ''                                                      ), &
      variableinfo( auxvar_denshydro_id, 'DENS_hyd', 'hydrostatic part of density',              &
                    'kg/m3', 3, 'XYZ', ''                                                     ), &
      variableinfo( auxvar_thermhydro_id, 'THERM_hyd', 'hydrostatic part of THERM',              &
                    'kg/m3', 3, 'XYZ', ''                                                     ), &
      variableinfo( auxvar_pres_id     ,     'PRES', 'pressure',                                 &
                      'Pa', 3, 'XYZ', 'air_pressure'                                          ), &
      variableinfo( auxvar_pt_id       ,       'PT', 'potential temperature',                    &
                        'K', 3, 'XYZ', 'potential_temperature'                                ), &
      variableinfo( auxvar_rtot_id     ,     'RTOT', 'Total gas constant',                       &
                        'J/kg/K', 3, 'XYZ', ''                                                ), &
      variableinfo( auxvar_cvtot_id    ,    'CVTOT', 'Total heat capacity',                      &
                        'J/kg/K', 3, 'XYZ', ''                                                ), &
      variableinfo( auxvar_cptot_id    ,    'CPTOT', 'Total heat capacity',                      &
                        'J/kg/K', 3, 'XYZ', ''                                                ), &
      variableinfo( auxvar_qdry_id     ,     'QDRY', 'dry air',                                  &
                        'kg/kg', 3, 'XYZ', ''                                                 ), &
      variableinfo( auxvar_preshydro_ref_id, 'PRES_hyd_REF', 'hydrostatic reference pressure',   &
                      'Pa', 3, 'XYZ', ''                                                        )/
 
    type(variableinfo) :: phytend_varinfo(phytend_num)
    DATA phytend_varinfo / &
      variableinfo( phytend_dens_id, 'DENS_tp', 'DENS_tp',                        &
                    'kg/m3/s',  3, 'XYZ',  'tendency of physical process for DENS' ),   &
      variableinfo( phytend_momx_id, 'MOMX_tp', 'MOMX_tp',                        &
                    'kg/m2/s',  3, 'XYZ',  'tendency of physical process for MOMX' ),   &
      variableinfo( phytend_momy_id, 'MOMY_tp', 'MOMY_tp',                        &
                    'kg/m2/s',  3, 'XYZ',  'tendency of physical process for MOMY' ),   &
      variableinfo( phytend_momz_id, 'MOMZ_tp', 'MOMZ_tp',                        &
                    'kg/m2/s',  3, 'XYZ',  'tendency of physical process for MOMZ' ),   &
      variableinfo( phytend_rhot_id, 'RHOT_tp', 'RHOT_tp',                        &
                    'kg/m3.K/s',  3, 'XYZ',  'tendency of physical process for RHOT' ), &
      variableinfo( phytend_rhoh_id,  'RHOH_p',  'RHOH_p',                        &
                    'kg/m3.J/s',  3, 'XYZ',  'heating of physical process for THERM' )   /
 
    !----------------------------------------------------------
 
    prgvar_info(:) = prgvar_varinfo
    auxvar_info(:) = auxvar_varinfo
    if ( present(phytend_info) ) phytend_info(:) = phytend_varinfo
 
    return
  end subroutine atm_dyn_dgm_nonhydro3d_common_get_varinfo
 
  !> Setup variable managers for atmospheric nonhydrostatic dynamical core
  subroutine atm_dyn_dgm_nonhydro3d_common_setup_variables( &
    prgvars, qtrcvars, auxvars, phytends,                             & ! (inout)
    prgvar_manager, qtrcvar_manager, auxvar_manager, phytend_manager, & ! (inout)
    reg_file_hist, do_setup_phytend, phytend_num_tot, mesh3d,         & ! (in)
    prgvar_varinfo )                                                    ! (out)
 
    use scale_atmos_hydrometeor, only: &
      atmos_hydrometeor_dry
    use scale_tracer, only: &
      qa, tracer_name, tracer_desc, tracer_unit    
    use scale_mesh_base3d, only: meshbase3d
    use scale_meshfield_base, only: meshfield3d
    implicit none
    integer, intent(in) :: phytend_num_tot
    type(meshfield3d), intent(inout) :: prgvars(prgvar_num)       !< Array of objects to manage prognostic variables
    type(meshfield3d), intent(inout) :: qtrcvars(0:qa)            !< Array of objects to manage tracer variables
    type(meshfield3d), intent(inout) :: auxvars(auxvar_num)       !< Array of objects to manage auxiliary variables
    type(meshfield3d), intent(inout) :: phytends(phytend_num_tot) !< Array of objects to manage physics tendency variables
    type(modelvarmanager), intent(inout) :: prgvar_manager        !< Object to manage prognostic variables
    type(modelvarmanager), intent(inout) :: qtrcvar_manager       !< Object to manage tracer variables
    type(modelvarmanager), intent(inout) :: auxvar_manager        !< Object to manage auxiliary variables
    type(modelvarmanager), intent(inout) :: phytend_manager       !< Object to manage physics tendency variables
    logical, intent(in) :: reg_file_hist                          !< Flag whether variables are registered for history file output or not
    logical, intent(in) :: do_setup_phytend                       !< Flag whether to setup physics tendency variables or not
    class(meshbase3d), intent(in) :: mesh3d                       !< Object to manage 3D mesh
    type(variableinfo), intent(out) :: prgvar_varinfo(prgvar_num) !< Variable information with prognostic variables
 
    type(variableinfo) :: auxvar_varinfo(auxvar_num)
    type(variableinfo) :: phytend_varinfo(phytend_num)
 
    integer :: iv
    integer :: iq
 
    type(variableinfo) :: qtrc_dry_vinfo_tmp
    type(variableinfo) :: qtrc_dry_tp_vinfo_tmp
    type(variableinfo) :: qtrc_vinfo_tmp
    type(variableinfo) :: qtrc_tp_vinfo_tmp
    !----------------------------------------------------------
 
    call atm_dyn_dgm_nonhydro3d_common_get_varinfo( prgvar_varinfo, auxvar_varinfo, phytend_varinfo ) ! (out)
 
    !- Initialize prognostic variables
 
    do iv = 1, prgvar_num
      call prgvar_manager%Regist(  &
        prgvar_varinfo(iv), mesh3d,                           & ! (in) 
        prgvars(iv),                                          & ! (inout)
        reg_file_hist,  monitor_flag=.true., fill_zero=.true. ) ! (out)
      !$acc update device( prgvars(iv) )
    end do
 
    !- Initialize tracer variables
 
    if ( qa == 0 ) then
      ! Dummy
      qtrc_dry_vinfo_tmp%ndims    = 3
      qtrc_dry_vinfo_tmp%dim_type = 'XYZ'
      qtrc_dry_vinfo_tmp%STDNAME  = ''
 
      qtrc_dry_vinfo_tmp%keyID = 0
      qtrc_dry_vinfo_tmp%NAME  = "QV"
      qtrc_dry_vinfo_tmp%DESC  = "Ratio of Water Vapor mass to total mass (Specific humidity)"
      qtrc_dry_vinfo_tmp%UNIT  = "kg/kg"
      call qtrcvar_manager%Regist( &
        qtrc_dry_vinfo_tmp, mesh3d,                      & ! (in) 
        qtrcvars(0),                                     & ! (inout)
        .false., monitor_flag=.false., fill_zero=.true.  ) ! (in)
    else
      qtrc_vinfo_tmp%ndims    = 3
      qtrc_vinfo_tmp%dim_type = 'XYZ'
      qtrc_vinfo_tmp%STDNAME  = ''
  
      do iq = 1, qa
        qtrc_vinfo_tmp%keyID = iq
        qtrc_vinfo_tmp%NAME  = tracer_name(iq)
        qtrc_vinfo_tmp%DESC  = tracer_desc(iq)
        qtrc_vinfo_tmp%UNIT  = tracer_unit(iq)
        call qtrcvar_manager%Regist( &
          qtrc_vinfo_tmp, mesh3d,                               & ! (in) 
          qtrcvars(iq),                                         & ! (inout)
          reg_file_hist, monitor_flag=.true., fill_zero=.true.  ) ! (in)
      end do
    end if
 
    !- Initialize auxiliary variables
 
    do iv = 1, auxvar_num
      call auxvar_manager%Regist( &
        auxvar_varinfo(iv), mesh3d,       & ! (in) 
        auxvars(iv),                      & ! (inout)
        reg_file_hist, fill_zero=.true.   ) ! (in)
    end do
 
    !- Initialize the tendency of physical processes
 
    if ( do_setup_phytend ) then
 
      do iv = 1, phytend_num    
        call phytend_manager%Regist( &
          phytend_varinfo(iv), mesh3d,     & ! (in) 
          phytends(iv),                    & ! (inout)
          reg_file_hist, fill_zero=.true.  ) ! (in)
      end do
 
      if ( qa == 0 ) then
        ! Dummy
        qtrc_dry_tp_vinfo_tmp%ndims    = 3
        qtrc_dry_tp_vinfo_tmp%dim_type = 'XYZ'
        qtrc_dry_tp_vinfo_tmp%STDNAME  = ''
 
        iv = phytend_num + 1
        qtrc_dry_tp_vinfo_tmp%keyID = iv
        qtrc_dry_tp_vinfo_tmp%NAME  = "QV_tp"
        qtrc_dry_tp_vinfo_tmp%DESC  = "tendency of physical process for QV"
        qtrc_dry_tp_vinfo_tmp%UNIT  = "kg/m3/s"
        call phytend_manager%Regist( &
          qtrc_dry_tp_vinfo_tmp, mesh3d,   & ! (in) 
          phytends(iv),                    & ! (inout)
          .false., fill_zero=.true.        ) ! (in)
      else
        qtrc_tp_vinfo_tmp%ndims    = 3
        qtrc_tp_vinfo_tmp%dim_type = 'XYZ'
        qtrc_tp_vinfo_tmp%STDNAME  = ''
 
        do iq = 1, qa
          iv = phytend_num + iq 
          qtrc_tp_vinfo_tmp%keyID = iv
          qtrc_tp_vinfo_tmp%NAME  = trim(tracer_name(iq))//'_tp'
          qtrc_tp_vinfo_tmp%DESC  = 'tendency of physical process for '//trim(tracer_desc(iq))
          qtrc_tp_vinfo_tmp%UNIT  = trim(tracer_unit(iq))//'/s'
 
          call phytend_manager%Regist( &
            qtrc_tp_vinfo_tmp, mesh3d,       & ! (in) 
            phytends(iv),                    & ! (inout)
            reg_file_hist, fill_zero=.true.  ) ! (in)
        end do    
      end if
 
    end if
 
    return
  end subroutine atm_dyn_dgm_nonhydro3d_common_setup_variables
 
!OCL SERIAL
  subroutine atm_dyn_dgm_nonhydro3d_common_calc_pressure( &
    PRES, DPRES,                                        & ! (inout)
    ddens, momx, momy, momz, therm,                     & ! (in)
    pres_hyd, dens_hyd, therm_hyd, rtot, cvtot, cptot,  & ! (in)
    mesh3d, entot_conserve_scheme_flag                  ) ! (in)
 
    implicit none
    class(meshfield3d), intent(inout) :: pres
    class(meshfield3d), intent(inout) :: dpres
    class(meshfield3d), intent(in) :: ddens
    class(meshfield3d), intent(in) :: momx
    class(meshfield3d), intent(in) :: momy
    class(meshfield3d), intent(in) :: momz
    class(meshfield3d), intent(in) :: therm
    class(meshfield3d), intent(in) :: pres_hyd
    class(meshfield3d), intent(in) :: dens_hyd
    class(meshfield3d), intent(in) :: therm_hyd
    class(meshfield3d), intent(in) :: rtot
    class(meshfield3d), intent(in) :: cvtot
    class(meshfield3d), intent(in) :: cptot
    class(meshbase3d), intent(in), target :: mesh3d
    logical, intent(in) :: entot_conserve_scheme_flag
    
    integer :: n
    class(localmesh3d), pointer :: lcmesh3d
    !---------------------------
 
    do n=1, mesh3d%LOCAL_MESH_NUM
      lcmesh3d => mesh3d%lcmesh_list(n)
 
      if ( entot_conserve_scheme_flag ) then
        call atm_dyn_dgm_nonhydro3d_common_entot2pres( pres%local(n)%val, dpres%local(n)%val,                  &
          ddens%local(n)%val, momx%local(n)%val, momy%local(n)%val, momz%local(n)%val, therm%local(n)%val,     &
          pres_hyd%local(n)%val, dens_hyd%local(n)%val, rtot%local(n)%val, cvtot%local(n)%val,                 &
          lcmesh3d, lcmesh3d%refElem3D )
      else
        call atm_dyn_dgm_nonhydro3d_common_drhot2pres( pres%local(n)%val, dpres%local(n)%val,                                           &
          therm%local(n)%val, pres_hyd%local(n)%val, therm_hyd%local(n)%val, rtot%local(n)%val, cvtot%local(n)%val, cptot%local(n)%val, &
          lcmesh3d, lcmesh3d%refElem3D )
      end if
    end do
 
    return
  end subroutine atm_dyn_dgm_nonhydro3d_common_calc_pressure
 
!OCL SERIAL
  subroutine atm_dyn_dgm_nonhydro3d_common_calc_therm_phyd( &
    THERM_hyd,                                 & ! (inout)
    pres_hyd, dens_hyd,                        & ! (in)
    mesh3d, entot_conserve_scheme_flag         ) ! (in)
 
    implicit none
    class(meshfield3d), intent(inout) :: therm_hyd
    class(meshfield3d), intent(in) :: pres_hyd
    class(meshfield3d), intent(in) :: dens_hyd
    class(meshbase3d), intent(in), target :: mesh3d
    logical, intent(in) :: entot_conserve_scheme_flag
    
    integer :: n
    class(localmesh3d), pointer :: lcmesh3d
    !---------------------------
 
    do n=1, mesh3d%LOCAL_MESH_NUM
      lcmesh3d => mesh3d%lcmesh_list(n)
 
      if ( entot_conserve_scheme_flag ) then
        ! ToDO
      else
        call atm_dyn_dgm_nonhydro3d_common_calc_rhot_hyd( therm_hyd%local(n)%val, &
          pres_hyd%local(n)%val, lcmesh3d, lcmesh3d%refElem3D                     )
      end if
    end do
 
    return
  end subroutine atm_dyn_dgm_nonhydro3d_common_calc_therm_phyd
 
  !> Calculate pressure from the deviation of density-weighted potential temperature (DRHOT)
!OCL SERIAL
  subroutine atm_dyn_dgm_nonhydro3d_common_drhot2pres( PRES, DPRES, &
    DRHOT, PRES_hyd, THERM_hyd, Rtot, CVtot, CPtot,                 &
    lcmesh, elem3D                                                  )
 
    use scale_const, only: &
      rdry => const_rdry,   &
      cpdry => const_cpdry, &
      cvdry => const_cvdry, &
      pres00 => const_pre00
        
    implicit none
    class(localmesh3d), intent(in) :: lcmesh
    class(elementbase3d), intent(in) :: elem3d
    real(rp), intent(out) :: pres(elem3d%np,lcmesh%nea) 
    real(rp), intent(out) :: dpres(elem3d%np,lcmesh%nea)
    real(rp), intent(in) :: drhot(elem3d%np,lcmesh%nea) 
    real(rp), intent(in) :: pres_hyd(elem3d%np,lcmesh%nea)
    real(rp), intent(in) :: therm_hyd(elem3d%np,lcmesh%nea)
    real(rp), intent(in) :: rtot(elem3d%np,lcmesh%nea) 
    real(rp), intent(in) :: cvtot(elem3d%np,lcmesh%nea) 
    real(rp), intent(in) :: cptot(elem3d%np,lcmesh%nea) 
 
    integer :: ke, p
    real(rp) :: rhot(elem3d%np)
    real(rp) :: rhot_
    real(rp) :: rp0
    !---------------------------------------------------------------
 
    rp0 = 1.0_rp / pres00
    !$omp parallel do private( RHOT )
    !$acc parallel loop gang present( PRES, DPRES, DRHOT, PRES_hyd, THERM_hyd, Rtot, CVtot, CPtot, lcmesh,elem3D )
    do ke=lcmesh%NeS, lcmesh%NeE
!      RHOT(:) = PRES00 / Rdry * ( PRES_hyd(:,ke) / PRES00 )**(CvDry/CpDry) + DRHOT(:,ke)
 
#ifdef _OPENACC
      !$acc loop vector
      do p=1, elem3d%Np
        rhot_ = pres00 / rdry * ( pres_hyd(p,ke) / pres00 )**(cvdry/cpdry) + drhot(p,ke)
 
        pres(p,ke) = pres00 * ( rtot(p,ke) * rp0 * rhot_ )**( cptot(p,ke) / cvtot(p,ke) )
        dpres(p,ke) = pres(p,ke) - pres_hyd(p,ke)
      end do
#else
      rhot(:) = therm_hyd(:,ke) + drhot(:,ke)
 
      pres(:,ke) = pres00 * ( rtot(:,ke) * rp0 * rhot(:) )**( cptot(:,ke) / cvtot(:,ke) )
      dpres(:,ke) = pres(:,ke) - pres_hyd(:,ke)
#endif
    end do
 
    return
  end subroutine atm_dyn_dgm_nonhydro3d_common_drhot2pres
 
!OCL SERIAL
  subroutine atm_dyn_dgm_nonhydro3d_common_drhot2entot( EnTot,       &
    DDENS, MOMX, MOMY, MOMZ, DRHOT,                                  &
    DENS_hyd, PRES_hyd, THERM_hyd, Rtot, CVtot, CPtot,               &
    lcmesh, elem3D                                                   )
 
    use scale_const, only: &
      grav => const_grav
        
    implicit none
    class(localmesh3d), intent(in) :: lcmesh
    class(elementbase3d), intent(in) :: elem3d
    real(rp), intent(out) :: entot(elem3d%np,lcmesh%nea) 
    real(rp), intent(in) :: ddens(elem3d%np,lcmesh%nea)
    real(rp), intent(in) :: momx(elem3d%np,lcmesh%nea)
    real(rp), intent(in) :: momy(elem3d%np,lcmesh%nea)
    real(rp), intent(in) :: momz(elem3d%np,lcmesh%nea)
    real(rp), intent(in) :: drhot(elem3d%np,lcmesh%nea) 
    real(rp), intent(in) :: dens_hyd(elem3d%np,lcmesh%nea)
    real(rp), intent(in) :: pres_hyd(elem3d%np,lcmesh%nea)
    real(rp), intent(in) :: therm_hyd(elem3d%np,lcmesh%nea)
    real(rp), intent(in) :: rtot(elem3d%np,lcmesh%nea) 
    real(rp), intent(in) :: cvtot(elem3d%np,lcmesh%nea) 
    real(rp), intent(in) :: cptot(elem3d%np,lcmesh%nea) 
 
    integer :: ke, ke2d
 
    real(rp) :: dens(elem3d%np)
    real(rp) :: mom_u1(elem3d%np), mom_u2(elem3d%np)
 
    real(rp) :: pres(elem3d%np,lcmesh%nea)
    real(rp) :: dpres(elem3d%np,lcmesh%nea)
    !---------------------------------------------------------------
 
    call atm_dyn_dgm_nonhydro3d_common_drhot2pres( pres, dpres, &
      drhot, pres_hyd, therm_hyd, rtot, cvtot, cptot,           &
      lcmesh, elem3d                                            )
    
    !$omp parallel do private( ke2D, DENS, mom_u1, mom_u2 )
    do ke=lcmesh%NeS, lcmesh%NeE
      ke2d = lcmesh%EMap3Dto2D(ke)
 
      dens(:) = dens_hyd(:,ke) + ddens(:,ke)
      mom_u1(:) = lcmesh%G_ij(lcmesh%refElem3D%IndexH2Dto3D,ke2d,1,1) * momx(:,ke) + lcmesh%G_ij(lcmesh%refElem3D%IndexH2Dto3D,ke2d,2,1) * momy(:,ke)
      mom_u2(:) = lcmesh%G_ij(lcmesh%refElem3D%IndexH2Dto3D,ke2d,2,1) * momx(:,ke) + lcmesh%G_ij(lcmesh%refElem3D%IndexH2Dto3D,ke2d,2,2) * momy(:,ke)
 
      entot(:,ke) = pres(:,ke) * cvtot(:,ke) / rtot(:,ke)                                                  &
                  + 0.5_rp * ( momx(:,ke) * mom_u1(:) + momy(:,ke) * mom_u2(:) + momz(:,ke)**2 ) / dens(:) &
                  + grav * dens(:) * lcmesh%zlev(:,ke)
    end do
 
    return
  end subroutine atm_dyn_dgm_nonhydro3d_common_drhot2entot
 
!OCL SERIAL
  subroutine atm_dyn_dgm_nonhydro3d_common_entot2pres( PRES, DPRES, &
    DDENS, MOMX, MOMY, MOMZ, EnTot,         &
    PRES_hyd, DENS_hyd, Rtot, CVtot,        &
    lcmesh, elem3D )
 
    use scale_const, only: &
      grav => const_grav
    
    implicit none
    class(localmesh3d), intent(in) :: lcmesh
    class(elementbase3d), intent(in) :: elem3d
    real(rp), intent(out) :: pres(elem3d%np,lcmesh%nea)
    real(rp), intent(out) :: dpres(elem3d%np,lcmesh%nea)
    real(rp), intent(in) :: ddens(elem3d%np,lcmesh%nea)
    real(rp), intent(in) :: momx(elem3d%np,lcmesh%nea)
    real(rp), intent(in) :: momy(elem3d%np,lcmesh%nea)
    real(rp), intent(in) :: momz(elem3d%np,lcmesh%nea)
    real(rp), intent(in) :: entot(elem3d%np,lcmesh%nea)
    real(rp), intent(in) :: pres_hyd(elem3d%np,lcmesh%nea)
    real(rp), intent(in) :: dens_hyd(elem3d%np,lcmesh%nea)
    real(rp), intent(in) :: rtot(elem3d%np,lcmesh%nea)
    real(rp), intent(in) :: cvtot(elem3d%np,lcmesh%nea)
 
    integer :: ke, ke2d
 
    real(rp) :: dens(elem3d%np)
    real(rp) :: mom_u1(elem3d%np), mom_u2(elem3d%np)
    !---------------------------------------------------------------
 
    !$omp parallel do private( ke2D, DENS, mom_u1, mom_u2 )
    do ke=lcmesh%NeS, lcmesh%NeE
      ke2d = lcmesh%EMap3Dto2D(ke)
 
      dens(:) = dens_hyd(:,ke) + ddens(:,ke)
      mom_u1(:) = lcmesh%G_ij(lcmesh%refElem3D%IndexH2Dto3D,ke2d,1,1) * momx(:,ke) + lcmesh%G_ij(lcmesh%refElem3D%IndexH2Dto3D,ke2d,2,1) * momy(:,ke)
      mom_u2(:) = lcmesh%G_ij(lcmesh%refElem3D%IndexH2Dto3D,ke2d,2,1) * momx(:,ke) + lcmesh%G_ij(lcmesh%refElem3D%IndexH2Dto3D,ke2d,2,2) * momy(:,ke)
 
      pres(:,ke) = (  entot(:,ke) - grav * dens(:) * lcmesh%zlev(:,ke)                                        &
                  - 0.5_rp * ( momx(:,ke) * mom_u1(:) + momy(:,ke) * mom_u2(:) + momz(:,ke)**2 ) / dens(:) &
                ) * rtot(:,ke) / cvtot(:,ke)
      
      dpres(:,ke) = pres(:,ke) - pres_hyd(:,ke)
    end do
    
    return
  end subroutine atm_dyn_dgm_nonhydro3d_common_entot2pres
 
!OCL SERIAL
  subroutine atm_dyn_dgm_nonhydro3d_common_calc_rhot_hyd( RHOT_hyd, &
    PRES_hyd,                                                       &
    lcmesh, elem3D                                                  )
 
    use scale_const, only: &
      rdry => const_rdry,   &
      cpdry => const_cpdry, &
      cvdry => const_cvdry, &
      pres00 => const_pre00
        
    implicit none
    class(localmesh3d), intent(in) :: lcmesh
    class(elementbase3d), intent(in) :: elem3d
    real(rp), intent(out) :: rhot_hyd(elem3d%np,lcmesh%nea) 
    real(rp), intent(in) :: pres_hyd(elem3d%np,lcmesh%nea)
 
    integer :: ke, p
    real(rp) :: rp0
    !---------------------------------------------------------------
 
    rp0 = 1.0_rp / pres00
 
    !$omp parallel do
    !$acc parallel loop gang present( RHOT_hyd, PRES_hyd, lcmesh, elem3D )
    do ke=lcmesh%NeS, lcmesh%NeE
#ifdef _OPENACC
      do p=1, elem3d%Np
        rhot_hyd(p,ke) = pres00 / rdry * ( pres_hyd(p,ke) / pres00 )**(cvdry/cpdry) 
      end do
#else
      rhot_hyd(:,ke) = pres00 / rdry * ( pres_hyd(:,ke) / pres00 )**(cvdry/cpdry) 
#endif
    end do
 
    return
  end subroutine atm_dyn_dgm_nonhydro3d_common_calc_rhot_hyd
 
!> Calculate horizontal graidient of hydrostatic pressure 
!! In this calculation, we assume that PRES_hyd_ref is continuous at element boundaries.
!OCL SERIAL
  subroutine atm_dyn_dgm_nonhydro3d_common_calc_phyd_hgrad_lc( DPhydDx, DPhydDy, &
      PRES_hyd, PRES_hyd_ref,        &
      element3D_operation, lmesh, elem )
    implicit none
    class(localmesh3d), intent(in) :: lmesh
    class(elementbase3d), intent(in) :: elem
    real(rp), intent(out)  :: dphyddx(elem%np,lmesh%nea)
    real(rp), intent(out)  :: dphyddy(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)
    class(elementoperationbase3d), intent(in) :: element3d_operation
 
    integer :: ke, ke2d, p  
 
    real(rp) :: flux(elem%np,3), fz(elem%np), dflux(elem%np,4,2)
    real(rp) :: del_flux_hyd(elem%nfptot,2,lmesh%ne)
    real(rp) :: gsqrtv, rgsqrtv(elem%np)
 
    real(rp) :: e33
    real(rp) :: gradphyd_x, gradphyd_y
    !-----------------------------------------
 
    call get_phyd_hgrad_numflux_generalhvc( del_flux_hyd,                       & ! (out)
      pres_hyd, pres_hyd_ref,                                                   & ! (in)
      lmesh%Gsqrt, 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, lmesh%lcmesh2D, lmesh%lcmesh2D%refElem2D                     ) ! (in)
 
    !$omp parallel do private( ke, ke2D, p, &
    !$omp Flux, Fz, DFlux, GsqrtV, RGsqrtV, E33, &
    !$omp GradPhyd_x, GradPhyd_y )      
    do ke = lmesh%NeS, lmesh%NeE
      ke2d = lmesh%EMap3Dto2D(ke)
 
      do p=1, elem%Np
        gsqrtv  = lmesh%Gsqrt(p,ke) / ( lmesh%gam(p,ke)**2 * lmesh%GsqrtH(elem%IndexH2Dto3D(p),ke2d) )
 
        rgsqrtv(p) = 1.0_rp / gsqrtv
        flux(p,1) = gsqrtv * ( pres_hyd(p,ke) - pres_hyd_ref(p,ke) )
      end do
 
      do p=1, elem%Np        
        flux(p,2) = flux(p,1)
        flux(p,3) = lmesh%GI3(p,ke,1) * flux(p,1)
        fz(p)   = lmesh%GI3(p,ke,2) * flux(p,1)
      end do   
      
      call element3d_operation%Div( flux, del_flux_hyd(:,1,ke), &
        dflux(:,:,1) )
      call element3d_operation%Dz( fz, dflux(:,3,2) )
      call element3d_operation%Lift( del_flux_hyd(:,2,ke), dflux(:,4,2) )
      
      do p=1, elem%Np
        e33 = lmesh%Escale(p,ke,3,3)
 
        gradphyd_x = lmesh%Escale(p,ke,1,1) * dflux(p,1,1) &
                   + e33                    * dflux(p,3,1) &
                   + dflux(p,4,1)
 
        gradphyd_y = lmesh%Escale(p,ke,2,2) * dflux(p,2,1) &
                   + e33                    * dflux(p,3,2) &
                   + dflux(p,4,2)
 
        dphyddx(p,ke) = gradphyd_x * rgsqrtv(p)
        dphyddy(p,ke) = gradphyd_y * rgsqrtv(p)
      end do      
    end do
 
    return
  end subroutine atm_dyn_dgm_nonhydro3d_common_calc_phyd_hgrad_lc
 
!-- private
 
!OCL SERIAL
  subroutine get_phyd_hgrad_numflux_generalhvc( &
    del_flux_hyd,                                        & ! (out)
    pres_hyd, pres_hyd_ref,                              & ! (in)
    gsqrt, gsqrth, gam, g13, g23, nx, ny, nz,            & ! (in)
    vmapm, vmapp, im2dto3d, lmesh, elem, lmesh2d, elem2d ) ! (in)
 
    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) ::  del_flux_hyd(elem%nfptot,2,lmesh%ne)
    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) ::  gsqrth(elem2d%np,lmesh2d%ne)
    real(rp), intent(in) ::  gsqrt(elem%np*lmesh%nea)
    real(rp), intent(in) :: gam(elem%np*lmesh%nea)
    real(rp), intent(in) ::  g13(elem%np*lmesh%nea)
    real(rp), intent(in) ::  g23(elem%np*lmesh%nea)
    real(rp), intent(in) :: nx(elem%nfptot,lmesh%ne)
    real(rp), intent(in) :: ny(elem%nfptot,lmesh%ne)
    real(rp), intent(in) :: nz(elem%nfptot,lmesh%ne)
    integer, intent(in) :: vmapm(elem%nfptot,lmesh%ne)
    integer, intent(in) :: vmapp(elem%nfptot,lmesh%ne)
    integer, intent(in) :: im2dto3d(elem%nfptot)
    
    integer :: ke, fp, i, ip(elem%nfptot), im(elem%nfptot)
    integer :: ke2d
    real(rp) :: dpres_hyd(elem%nfptot,2)
    real(rp) :: gsqrt_(elem%nfptot,2)
    real(rp) :: gsqrtv_(elem%nfptot,2)
    real(rp) :: g13_(elem%nfptot,2)
    real(rp) :: g23_(elem%nfptot,2)
    
    integer, parameter :: in = 1
    integer, parameter :: ex = 2
 
    real(rp) :: tmp1, tmp2
    !------------------------------------------------------------------------
 
    !$omp parallel do private( &
    !$omp ke, iM, iP, ke2D, fp,                    &
    !$omp DPRES_hyd,  Gsqrt_, GsqrtV_, G13_, G23_, &
    !$omp tmp1, tmp2 )
!OCL PREFETCH
    do ke=lmesh%NeS, lmesh%NeE
      im(:) = vmapm(:,ke); ip(:) = vmapp(:,ke)
      ke2d = lmesh%EMap3Dto2D(ke)
 
      gsqrt_(:,in) = gsqrt(im)
      gsqrt_(:,ex) = gsqrt(ip)
      gsqrtv_(:,in) = gsqrt_(:,in) / gam(im)**2 / gsqrth(im2dto3d(:),ke2d)
      gsqrtv_(:,ex) = gsqrt_(:,ex) / gam(ip)**2 / gsqrth(im2dto3d(:),ke2d)
 
      g13_(:,in) = g13(im)
      g13_(:,ex) = g13(ip)
      g23_(:,in) = g23(im)
      g23_(:,ex) = g23(ip)
 
      dpres_hyd(:,in) = pres_hyd(im) - pres_hyd_ref(im)
      dpres_hyd(:,ex) = pres_hyd(ip) - pres_hyd_ref(ip) 
      
      do fp=1, elem%NfpTot
        tmp1 = lmesh%Fscale(fp,ke) * 0.5_rp * gsqrtv_(fp,ex) * dpres_hyd(fp,ex)
        tmp2 = lmesh%Fscale(fp,ke) * 0.5_rp * gsqrtv_(fp,in) * dpres_hyd(fp,in)
 
        del_flux_hyd(fp,1,ke) = &
            ( nx(fp,ke) + g13_(fp,ex) * nz(fp,ke) ) * tmp1 &
          - ( nx(fp,ke) + g13_(fp,in) * nz(fp,ke) ) * tmp2
 
        del_flux_hyd(fp,2,ke) = &
            ( ny(fp,ke) + g23_(fp,ex) * nz(fp,ke) ) * tmp1 &
          - ( ny(fp,ke) + g23_(fp,in) * nz(fp,ke) ) * tmp2
      end do
    end do
 
    return
  end subroutine get_phyd_hgrad_numflux_generalhvc
 
end module scale_atm_dyn_dgm_nonhydro3d_common