mod_atmos_vars_container.F90

Go to the documentation of this file.

!-------------------------------------------------------------------------------
!> module Atmosphere / Variables
!!
!! @par Description
!!          Container for atmospheric variables
!!
!! @author Yuta Kawai, Team SCALE
!!
!<
!-------------------------------------------------------------------------------
#include "scaleFElib.h"
module mod_atmos_vars_container
  !-----------------------------------------------------------------------------
  !
  !++ Used modules
  !
  use scale_precision
  use scale_io
  use scale_prc
  use scale_tracer, only: qa
 
  use scale_element_base, only: &
    elementbase, elementbase2d, elementbase3d
  use scale_localmesh_base, only: localmeshbase
  use scale_localmesh_2d, only: localmesh2d
  use scale_localmesh_3d, only: localmesh3d
  use scale_localmeshfield_base, only: &
    localmeshfieldbase, localmeshfieldbaselist
  use scale_mesh_base, only: meshbase
  use scale_mesh_base2d, only: meshbase2d
  use scale_mesh_base3d, only: meshbase3d
 
  use scale_meshfield_base, only: &
    meshfieldbase, meshfield2d, meshfield3d
  use scale_model_var_manager, only: &
    modelvarmanager, variableinfo
 
  use scale_atm_dyn_dgm_nonhydro3d_common, only: &
    prgvar_num, auxvar_num, phytend_num1 => phytend_num,                                        &
    prgvar_ddens_id, prgvar_therm_id, prgvar_momz_id, prgvar_momx_id, prgvar_momy_id,           &
    auxvar_denshydro_id, auxvar_preshydro_id, auxvar_thermhydro_id, auxvar_preshydro_ref_id,    &
    auxvar_rtot_id, auxvar_cptot_id, auxvar_cvtot_id,                                           &
    auxvar_pres_id, auxvar_pt_id, auxvar_qdry_id,                                               &
    phytend_dens_id, phytend_momx_id, phytend_momy_id, phytend_momz_id, phytend_rhot_id,        &
    phytend_rhoh_id
 
  use mod_atmos_mesh, only: atmosmesh  
  use mod_atmos_phy_preproc, only: &
    physpreprocbase,        &
    physpreprocnone,        &
    physpreprocmodalfilter, &
    physpreprocglobalfilter
  !-----------------------------------------------------------------------------
  implicit none
  private
 
  !-----------------------------------------------------------------------------
  !
  !++ Public type & procedures
  !
 
  !> Derived type to manage a set of variables (prognostic variables, tracer variables, auxiliary variables, and tendencies of physical processes) with atmospheric component
  type, public :: atmosvarscontainer
    class(atmosmesh), pointer :: mesh              !< Pointer to an object to manage mesh for atmospheric component
 
    !- prognostic variables
    type(meshfield3d), allocatable :: prog_vars(:) !< Array of 3D prognostic variables
    type(modelvarmanager) :: progvars_manager      !< Object to manage 3D prognostic variables
    integer :: prog_vars_commid
 
    !- tracer variables
    type(meshfield3d), allocatable :: qtrc_vars(:) !< Array of 3D tracer variables
    type(modelvarmanager) :: qtrcvars_manager      !< Object to manage 3D tracer variables
    integer :: qtrc_vars_commid 
 
    !- auxiliary variables    
    type(meshfield3d), allocatable :: aux_vars(:)  !< Array of 3D auxiliary variables
    type(modelvarmanager) :: auxvars_manager       !< Object to manage 3D auxiliary variables
    integer :: aux_vars_commid
 
    !- auxiliary variables (2D)
    type(meshfield2d), allocatable :: aux_vars2d(:) !< Array of 2D auxiliary variables
    type(modelvarmanager) :: auxvars2d_manager      !< Object to manage 3D auxiliary variables
 
    !- Tendency with physics
    type(meshfield3d), allocatable :: phy_tend(:)  !< Array of tendency variables with physics
    type(modelvarmanager) :: phytends_manager      !< Object to manage tendency variables with physics
    integer :: phytends_commid
    integer :: phytend_num_tot                     !< Total number of tendency variables with physics
 
    !-
    integer :: container_type
    integer :: phy_preproc_operation_type
    class(physpreprocbase), allocatable :: phy_preproc
 
  contains
    procedure :: init => atmosvarscontainer_init
    procedure :: final => atmosvarscontainer_final
    procedure :: calc_diagnostics => atmosvarscontainer_calculatediagnostics
    procedure :: calc_diagvar => atmosvarscontainer_calcdiagvar
    procedure :: calc_specificheat => atmosvarscontainer_calc_specific_heat
    procedure :: preproc_operation_for_phys => atmosvarscontainer_physics_preoperation
 
    procedure, private :: setup_phys_preoperation
  end type atmosvarscontainer
 
  integer, parameter, public :: atm_vars_container_primary_id = 1
 
  public :: atmosvars_getlocalmeshprgvar
  public :: atmosvars_getlocalmeshprgvars
  public :: atmosvars_getlocalmeshsfcvar
  public :: atmosvars_getlocalmeshqtrcvar  
  public :: atmosvars_getlocalmeshqtrcvarlist
  public :: atmosvars_getlocalmeshqtrc_qv
  public :: atmosvars_getlocalmeshphyauxvars
  public :: atmosvars_getlocalmeshphytends
  public :: atmosvars_getlocalmeshqtrcphytend
 
  !-----------------------------------------------------------------------------
  !
  !++ Public parameters & variables
  !
  !-----------------------------------------------------------------------------
  !
  !++ Private procedures & variables
  !
  private :: vars_calc_diagnosevar_lc
 
  ! Surface variables
  
  integer, public, parameter :: atmos_auxvars2d_prec_id      = 1
  integer, public, parameter :: atmos_auxvars2d_prec_engi_id = 2
  integer, public, parameter :: atmos_auxvars2d_num          = 2
 
  type(variableinfo), public :: atmos_auxvars2d_vinfo(atmos_auxvars2d_num)
  DATA atmos_auxvars2d_vinfo / &
    variableinfo( atmos_auxvars2d_prec_id     ,      'PREC', 'surface precipitaion flux'        , 'kg/m2/s', 2, 'XY', 'precipitation_flux'  ), &
    variableinfo( atmos_auxvars2d_prec_engi_id, 'PREC_ENGI', 'internal energy of precipitation' ,    'J/m2', 2, 'XY', ''  )                    /
 
  !
  integer, parameter :: phy_preoperation_typeid_none         = 0
  integer, parameter :: phy_preoperation_typeid_modalfilter  = 1
  integer, parameter :: phy_preoperation_typeid_globalfilter = 2
 
contains
 
  !> Initalize an object to manage a set of variables with atmospheric component
!OCL SERIAL
  subroutine atmosvarscontainer_init( this, &
    container_type, phy_preproc_file_basename,     &
    atm_mesh )
    use scale_const, only: &
       undef => const_undef    
    use scale_tracer, only: &
      tracer_name, tracer_desc, tracer_unit
    use scale_atm_dyn_dgm_nonhydro3d_common, only: &
      prgvar_scalar_num, prgvar_hvec_num,          &
      atm_dyn_dgm_nonhydro3d_common_setup_variables
    
    implicit none
    class(atmosvarscontainer), intent(inout) :: this
    integer, intent(in) :: container_type
    character(len=*), intent(in) :: phy_preproc_file_basename
    class(atmosmesh), target, intent(inout) :: atm_mesh
 
    integer :: iv
    integer :: idom
 
    type(variableinfo) :: prgvar_info(prgvar_num)
    logical :: do_setup_phytend
    logical :: do_setup_auxvar2d
    logical :: reg_file_hist
 
    class(meshbase3d), pointer :: mesh3d
    class(meshbase2d), pointer :: mesh2d
    !---------------------------------------------------------
    
    if ( container_type == atm_vars_container_primary_id ) then
      do_setup_auxvar2d = .true.
      do_setup_phytend = .true.
      reg_file_hist = .true.
    else
      do_setup_auxvar2d = .false.
      do_setup_phytend = .false.
      reg_file_hist = .false.
    end if
 
    this%mesh => atm_mesh
    mesh3d => atm_mesh%ptr_mesh
    call mesh3d%GetMesh2D( mesh2d )
 
    !-
    call this%PROGVARS_manager%Init()
    call this%QTRCVARS_manager%Init()
    call this%AUXVARS_manager%Init()
 
    allocate( this%PROG_VARS(prgvar_num) )
    allocate( this%QTRC_VARS(0:qa) )
    allocate( this%AUX_VARS(auxvar_num) )
    !$acc enter data create( this%PROG_VARS, this%QTRC_VARS, this%AUX_VARS )    
 
    if ( do_setup_phytend ) then
      call this%PHYTENDS_manager%Init()
 
      this%PHYTEND_NUM_TOT = phytend_num1 + max(1,qa)
      allocate( this%PHY_TEND(this%PHYTEND_NUM_TOT) )
      !$acc enter data create( this%PHY_TEND )
    end if
 
    call atm_dyn_dgm_nonhydro3d_common_setup_variables( &
      this%PROG_VARS, this%QTRC_VARS, this%AUX_VARS, this%PHY_TEND,                              & ! (inout)
      this%PROGVARS_manager, this%QTRCVARS_manager, this%AUXVARS_manager, this%PHYTENDS_manager, & ! (inout)
      reg_file_hist, do_setup_phytend, this%PHYTEND_NUM_TOT, mesh3d,                             & ! (in)
      prgvar_info ) ! (out)
 
    ! Setup communicator
    
    call atm_mesh%Create_communicator( &
      prgvar_scalar_num, prgvar_hvec_num, 0,              & ! (in)
      this%PROGVARS_manager,                              & ! (inout)
      this%PROG_VARS,                                     & ! (in)
      this%PROG_VARS_commID                               ) ! (out)
        
 
    if ( qa > 0 ) then
      call atm_mesh%Create_communicator( &
        qa, 0, 0,                        & ! (in)
        this%QTRCVARS_manager,           & ! (inout)
        this%QTRC_VARS(1:qa),            & ! (in)
        this%QTRC_VARS_commID            ) ! (out)
    end if
 
    call atm_mesh%Create_communicator( &
      auxvar_num, 0, 0,                & ! (in)
      this%AUXVARS_manager,            & ! (inout)
      this%AUX_VARS,                   & ! (in)
      this%AUX_VARS_commID             ) ! (out)
 
    ! Output list of prognostic variables
      
    if ( container_type == atm_vars_container_primary_id ) then
      log_newline
      log_info("ATMOS_vars_setup",*) 'List of prognostic variables (ATMOS) '
      log_info_cont('(1x,A,A24,A,A48,A,A12,A)') &
                '      |', 'VARNAME                 ','|', &
                'DESCRIPTION                                     ', '[', 'UNIT        ', ']'
      do iv = 1, prgvar_num
        log_info_cont('(1x,A,I3,A,A24,A,A48,A,A12,A)') &
        'NO.',iv,'|',prgvar_info(iv)%NAME,'|', prgvar_info(iv)%DESC,'[', prgvar_info(iv)%UNIT,']'
      end do
      do iv = 1, qa
        log_info_cont('(1x,A,I3,A,A24,A,A48,A,A12,A)') &
        'NO.',prgvar_num+iv,'|',tracer_name(iv),'|', tracer_desc(iv),'[', tracer_unit(iv),']'
      end do
      log_newline
    end if
    
    !- Initialize 2D auxiliary variables
    if ( do_setup_auxvar2d ) then
      call this%AUXVARS2D_manager%Init()
      allocate( this%AUX_VARS2D(atmos_auxvars2d_num) )
      !$acc enter data create( this%AUX_VARS2D )
      
      do iv = 1, atmos_auxvars2d_num
        call this%AUXVARS2D_manager%Regist(    &
          atmos_auxvars2d_vinfo(iv), mesh2d,   & ! (in) 
          this%AUX_VARS2D(iv),                 & ! (inout)
          reg_file_hist                        ) ! (in)
        do idom=1, mesh2d%LOCAL_MESH_NUM
          this%AUX_VARS2D(iv)%local(idom)%val(:,:) = undef
        end do
      end do
    end if
 
    !- Setup preoperation before physics
    this%container_type = container_type
    call this%Setup_phys_preoperation( phy_preproc_file_basename, mesh3d, mesh3d%refElem3D )
 
    return
  end subroutine atmosvarscontainer_init
 
  !> Finalize an object to manage a set of variables with atmospheric component
!OCL SERIAL
  subroutine atmosvarscontainer_final( this )
    implicit none
    class(atmosvarscontainer), intent(inout) :: this
    !----------------------------------------
 
    !$acc exit data delete( this%PROG_VARS, this%QTRC_VARS, this%AUX_VARS )
 
    call this%PROGVARS_manager%Final()
    deallocate( this%PROG_VARS )
    
    call this%QTRCVARS_manager%Final()
    deallocate( this%QTRC_VARS )
 
    call this%AUXVARS_manager%Final()
    deallocate( this%AUX_VARS )
 
    if ( this%container_type == atm_vars_container_primary_id ) then
      !$acc exit data delete( this%PHY_TEND, this%AUX_VARS2D )
      
      call this%PHYTENDS_manager%Final()
      deallocate( this%PHY_TEND )
 
      call this%AUXVARS2D_manager%Final()
      deallocate( this%AUX_VARS2D )
    else
      ! Finalize preoperation before physics
      call this%phy_preproc%Final()
    end if
 
    return
  end subroutine atmosvarscontainer_final
 
  !> Calculate diagnostic variables from prognostic variables and auxiliary variables
!OCL SERIAL  
  subroutine atmosvarscontainer_calculatediagnostics( this )    
    use scale_const, only: &
      rdry => const_rdry,      &
      cpdry => const_cpdry,    &
      cvdry => const_cvdry,    &
      pres00 => const_pre00
    use scale_tracer, only: &
      tracer_mass, tracer_r, tracer_cv, tracer_cp    
    use scale_atmos_thermodyn, only: &
      atmos_thermodyn_specific_heat
    implicit none
    class(atmosvarscontainer), intent(inout), target :: this
 
    class(localmesh3d), pointer :: lcmesh3d
    integer :: n
    integer :: varid
 
    class(meshfield3d), pointer :: field
    class(elementbase3d), pointer :: elem3d
 
    type(localmeshfieldbaselist) :: qtrc(qa)
    !-------------------------------------------------------
 
    ! Calculate specific heat
    call this%Calc_SpecificHeat()
    
    ! Calculate diagnostic variables
    do varid=auxvar_thermhydro_id+1, auxvar_pt_id
      field => this%AUX_VARS(varid)
      do n=1, field%mesh%LOCAL_MESH_NUM
        call atmosvars_getlocalmeshqtrcvarlist( n, &
          field%mesh, this%QTRCVARS_manager,       &
          1, qtrc, lcmesh3d )
 
        elem3d => lcmesh3d%refElem3D
 
        call vars_calc_diagnosevar_lc( &
          field%varname, field%local(n)%val,                       &
          this%PROG_VARS(prgvar_ddens_id)%local(n)%val,            &
          this%PROG_VARS(prgvar_momx_id)%local(n)%val,             &
          this%PROG_VARS(prgvar_momy_id)%local(n)%val,             &
          this%PROG_VARS(prgvar_momz_id)%local(n)%val,             &
          this%AUX_VARS(auxvar_pres_id)%local(n)%val,              &
          this%AUX_VARS(auxvar_qdry_id)%local(n)%val,              &
          qtrc,                                                    &
          this%AUX_VARS(auxvar_denshydro_id)%local(n)%val,         & 
          this%AUX_VARS(auxvar_preshydro_id)%local(n)%val,         &
          this%AUX_VARS(auxvar_rtot_id )%local(n)%val,             & 
          this%AUX_VARS(auxvar_cvtot_id)%local(n)%val,             & 
          this%AUX_VARS(auxvar_cptot_id)%local(n)%val,             & 
          lcmesh3d, lcmesh3d%refElem3D )
      end do
    end do
    !$acc wait(1)
 
    return
  end subroutine atmosvarscontainer_calculatediagnostics
 
  !> Calcualte a diagnostic variable from prognostic variables and auxiliary variables
!OCL SERIAL
  subroutine atmosvarscontainer_calcdiagvar( this, field_name, field_work ) 
    use scale_const, only: &
      rdry => const_rdry,      &
      cpdry => const_cpdry,    &
      cvdry => const_cvdry,    &
      pres00 => const_pre00
    use scale_tracer, only: &
      tracer_mass, tracer_r, tracer_cv, tracer_cp    
    use scale_atmos_thermodyn, only: &
      atmos_thermodyn_specific_heat
 
    implicit none
    class(atmosvarscontainer), intent(inout) :: this
    character(*), intent(in) :: field_name         !< Name of the diagnostic variable to be calculated
    type(meshfield3d), intent(inout) :: field_work !< A work space for the diagnostic variable to be calculated
 
    class(localmesh3d), pointer :: lcmesh3d
    integer :: n
 
    type(meshfield3d) :: field_work_uvmet(2)
    logical :: is_uvmet
    integer :: uvmet_i
 
    type(localmeshfieldbaselist) :: qtrc(qa)
    !--------------------------------------------------
 
    is_uvmet = .false.
    if ( field_name == 'Umet' ) then
      is_uvmet = .true.; uvmet_i = 1
    else if ( field_name == 'Vmet' ) then      
      is_uvmet = .true.; uvmet_i = 2
    end if
 
    field_work%varname = field_name
 
    do n=1, field_work%mesh%LOCAL_MESH_NUM
      call atmosvars_getlocalmeshqtrcvarlist( n, &
        field_work%mesh, this%QTRCVARS_manager,  &
        1, qtrc, lcmesh3d )
      
      if ( .not. is_uvmet ) then
        call vars_calc_diagnosevar_lc( field_name, field_work%local(n)%val,  &
          this%PROG_VARS(prgvar_ddens_id)%local(n)%val,               &
          this%PROG_VARS(prgvar_momx_id)%local(n)%val,                &
          this%PROG_VARS(prgvar_momy_id)%local(n)%val,                &
          this%PROG_VARS(prgvar_momz_id)%local(n)%val,                &
          this%AUX_VARS(auxvar_pres_id)%local(n)%val,                 &
          this%AUX_VARS(auxvar_qdry_id)%local(n)%val,                 &
          qtrc,                                                       &
          this%AUX_VARS(auxvar_denshydro_id)%local(n)%val,            & 
          this%AUX_VARS(auxvar_preshydro_id)%local(n)%val,            &
          this%AUX_VARS(auxvar_rtot_id )%local(n)%val,                & 
          this%AUX_VARS(auxvar_cvtot_id)%local(n)%val,                & 
          this%AUX_VARS(auxvar_cptot_id)%local(n)%val,                & 
          lcmesh3d, lcmesh3d%refElem3D )
      else
        call field_work_uvmet(1)%Init( 'Umet', '', field_work%mesh )
        call field_work_uvmet(2)%Init( 'Vmet', '', field_work%mesh )
        
        call this%mesh%Calc_UVmet( &
          this%PROG_VARS(prgvar_momx_id), this%PROG_VARS(prgvar_momy_id), & ! (in)
          field_work_uvmet(1), field_work_uvmet(2)                        ) ! (inout)
        call momtovel( field_work%local(n)%val,            & ! (out)
          field_work_uvmet(uvmet_i)%local(n)%val,          & ! (in)
          this%PROG_VARS(prgvar_ddens_id)%local(n)%val,    & ! (in)
          this%AUX_VARS(auxvar_denshydro_id)%local(n)%val, & ! (in)
          lcmesh3d, lcmesh3d%refElem3D ) ! (in)
        
        call field_work_uvmet(1)%Final()
        call field_work_uvmet(2)%Final()
      end if
    end do
    !$acc wait(1)
    return
  contains
    subroutine momtovel( vel, &
      mom, ddens, dens_hyd, lmesh, elem )
      implicit none
      class(localmesh3d), intent(in) :: lmesh
      class(elementbase3d), intent(in) :: elem
      real(RP), intent(inout) :: vel(elem%Np,lmesh%NeA)
      real(RP), intent(in) :: mom(elem%Np,lmesh%NeA)
      real(RP), intent(in) :: ddens(elem%Np,lmesh%NeA)
      real(RP), intent(in) :: dens_hyd(elem%Np,lmesh%NeA)
 
      integer :: ke, p
      !-----------------------------------------------
 
      !$omp parallel do
      !$acc parallel loop collapse(2) present(vel, mom, ddens, dens_hyd) async(1)
      do ke =lmesh%NeS, lmesh%NeE
      do p = 1, elem%Np
        vel(p,ke) = mom(p,ke) / ( ddens(p,ke) + dens_hyd(p,ke) )
      end do
      end do
 
      return
    end subroutine momtovel
  end subroutine atmosvarscontainer_calcdiagvar
 
  !> Calculate specific heat coefficients which are weighted by the mass of dry air and tracers
!OCL SERIAL  
  subroutine atmosvarscontainer_calc_specific_heat( this )
    use scale_const, only: &
      rdry => const_rdry,      &
      cpdry => const_cpdry,    &
      cvdry => const_cvdry,    &
      pres00 => const_pre00
    use scale_tracer, only: &
      tracer_mass, tracer_r, tracer_cv, tracer_cp    
    use scale_atmos_thermodyn, only: &
      atmos_thermodyn_specific_heat
    implicit none
    class(atmosvarscontainer), intent(inout), target :: this
 
    class(meshbase3d), pointer :: mesh3D
    class(localmesh3d), pointer :: lcmesh3D
    integer :: n
    integer :: ke, p
    integer :: iq
 
    class(elementbase3d), pointer :: elem3D
 
    real(RP), allocatable :: q_tmp(:,:)
    class(localmeshfieldbase), pointer :: Qdry, Rtot, CVtot, CPtot
    type(localmeshfieldbaselist) :: QTRC(QA)    
    !-------------------------------------------------------
 
    mesh3d => this%AUX_VARS(1)%mesh
 
    ! Calculate specific heat
    do n=1, mesh3d%LOCAL_MESH_NUM
      lcmesh3d => mesh3d%lcmesh_list(n)
      elem3d => lcmesh3d%refElem3D
      !$acc enter data attach(lcmesh3D, elem3D)
 
      allocate( q_tmp(elem3d%Np,qa) )
 
      call atmosvars_getlocalmeshqtrcvarlist( n, &
        mesh3d, this%QTRCVARS_manager,  &
        1, qtrc, lcmesh3d )
#ifdef _OPENACC
      do iq=1, qa
        !$acc enter data attach(QTRC(iq)%ptr)
      end do
#endif        
 
      qdry => this%AUX_VARS(auxvar_qdry_id )%local(n)
      rtot => this%AUX_VARS(auxvar_rtot_id )%local(n)
      cvtot => this%AUX_VARS(auxvar_cvtot_id)%local(n)
      cptot => this%AUX_VARS(auxvar_cptot_id)%local(n)
      !$acc enter data attach(Qdry, Rtot, CVtot, CPtot)
 
      !$omp parallel do private(ke, iq, q_tmp)
      !$acc parallel loop gang private(q_tmp) present(Qdry%val, Rtot%val, CVtot%val, CPtot%val, TRACER_MASS, TRACER_R, TRACER_CV, TRACER_CP, lcmesh3D,elem3D)
      do ke = lcmesh3d%NeS, lcmesh3d%NeE
        do iq = 1, qa
        !$acc loop vector
        do p=1, elem3d%Np
          q_tmp(p,iq) = qtrc(iq)%ptr%val(p,ke)
        end do
        end do
 
        call atmos_thermodyn_specific_heat( &
          elem3d%Np, 1, elem3d%Np, qa,                                       & ! (in)
          q_tmp, tracer_mass, tracer_r, tracer_cv, tracer_cp,                & ! (in)
          qdry%val(:,ke), rtot%val(:,ke), cvtot%val(:,ke), cptot%val(:,ke)   ) ! (out)
      end do
 
      deallocate(q_tmp)
    end do
    return
  end subroutine atmosvarscontainer_calc_specific_heat
 
  !> Perform preprocessing operation for variables used in physics
!OCL SERIAL  
  subroutine atmosvarscontainer_physics_preoperation( this, &
    container_ori, dyncore )
    use scale_file_history_meshfield, only: &
      file_history_meshfield_in
    use scale_atm_dyn_dgm_driver_nonhydro3d, only: atmdyndgmdriver_nonhydro3d
 
    use scale_atmos_hydrometeor, only: &
      qla, qia
    use scale_atm_phy_mp_dgm_common, only: &
      atm_phy_mp_dgm_common_negative_fixer    
    implicit none
    class(atmosvarscontainer), intent(inout), target :: this
    class(atmosvarscontainer), intent(in), target :: container_ori
    class(atmdyndgmdriver_nonhydro3d), intent(in) :: dyncore
 
    integer :: n
    integer :: iv, iq
 
    character(len=H_SHORT) :: varname_ori
    character(len=H_SHORT) :: typeid_s
 
    class(meshbase3d), pointer :: mesh3D
    class(localmesh3d), pointer :: lcmesh3D
    integer :: trcid_list(QA)
    type(localmeshfieldbaselist) :: lc_qtrc(QA)
    !----------------------------------------------
 
    if ( this%container_type == atm_vars_container_primary_id ) return
 
    log_info("AtmosVarsContainer_physics_preoperation",*) "..."
    
    call this%phy_preproc%Operate( this%PROG_VARS, this%QTRC_VARS, this%AUX_VARS,  & ! (inout)
      container_ori%PROG_VARS, container_ori%QTRC_VARS, container_ori%AUX_VARS ) ! (in)
 
    call this%Calc_SpecificHeat()
 
    call dyncore%update_therm_hyd( this%AUXVARS_manager )
    call dyncore%calc_pressure( this%AUX_VARS(auxvar_pres_id), &
      this%PROGVARS_manager, this%AUXVARS_manager              )
 
    !- Apply negative fixer for prognostic variables and auxiliary variables which are used in physics
    do iq=1, qa 
      trcid_list(iq) = iq
    end do
 
    mesh3d => this%AUX_VARS(1)%mesh
    do n=1, mesh3d%LOCAL_MESH_NUM
      call this%QTRCVARS_manager%GetLocalMeshFieldList( trcid_list, n, lc_qtrc )
 
      lcmesh3d => mesh3d%lcmesh_list(n)
      call atm_phy_mp_dgm_common_negative_fixer( &
        lc_qtrc, this%PROG_VARS(prgvar_ddens_id)%local(n)%val, this%AUX_VARS(auxvar_pres_id)%local(n)%val, &
        this%AUX_VARS(auxvar_cvtot_id)%local(n)%val, this%AUX_VARS(auxvar_cptot_id)%local(n)%val, &
        this%AUX_VARS(auxvar_rtot_id)%local(n)%val, &        
        this%AUX_VARS(auxvar_denshydro_id)%local(n)%val, this%AUX_VARS(auxvar_preshydro_id)%local(n)%val, &
        1.0_rp, lcmesh3d, lcmesh3d%refElem3D, qa, qla, qia )
    end do
 
    !-
    call this%Calc_diagnostics()
 
    !- Tentative : Register the variables for file history
 
    write(typeid_s,'(I2.2)') this%container_type
    varname_ori = this%PROG_VARS(prgvar_ddens_id)%varname
    this%PROG_VARS(prgvar_ddens_id)%varname = trim(varname_ori) // "_"//trim(typeid_s)
    call file_history_meshfield_in( this%PROG_VARS(prgvar_ddens_id), "" )
    this%PROG_VARS(prgvar_ddens_id)%varname = varname_ori
 
    varname_ori = this%AUX_VARS(auxvar_pres_id)%varname
    this%AUX_VARS(auxvar_pres_id)%varname = trim(varname_ori) // "_"//trim(typeid_s)
    call file_history_meshfield_in( this%AUX_VARS(auxvar_pres_id), "" )
    this%AUX_VARS(auxvar_pres_id)%varname = varname_ori
 
    do iv=1, 3
      varname_ori = this%QTRC_VARS(iv)%varname
      this%QTRC_VARS(iv)%varname = trim(varname_ori) // "_"//trim(typeid_s)
      call file_history_meshfield_in( this%QTRC_VARS(iv), "" )
      this%QTRC_VARS(iv)%varname = varname_ori
    end do
 
    return
  end subroutine atmosvarscontainer_physics_preoperation
 
  !----  Getter ---------------------------------------------------------------------------
 
!OCL SERIAL
  subroutine atmosvars_getlocalmeshprgvar( domID, mesh, prgvars_list, auxvars_list, &
     varid,                                                                         &
     var, DENS_hyd, PRES_hyd, lcmesh3D                                              )
   
    implicit none
 
    integer, intent(in) :: domid
    class(meshbase), intent(in) :: mesh
    class(modelvarmanager), intent(inout) :: prgvars_list
    class(modelvarmanager), intent(inout) :: auxvars_list
    integer, intent(in) :: varid
    class(localmeshfieldbase), pointer, intent(out) :: var
    class(localmeshfieldbase), pointer, intent(out), optional :: dens_hyd, pres_hyd
    class(localmesh3d), pointer, intent(out), optional :: lcmesh3d
 
    class(meshfieldbase), pointer :: field
    class(localmeshbase), pointer :: lcmesh
    !-------------------------------------------------------
 
    !--
    call prgvars_list%Get(varid, field)
    call field%GetLocalMeshField(domid, var)
 
    if (present(dens_hyd)) then
      call auxvars_list%Get(auxvar_denshydro_id, field)
      call field%GetLocalMeshField(domid, dens_hyd)
    end if
    if (present(pres_hyd)) then
      call auxvars_list%Get(auxvar_preshydro_id, field)
      call field%GetLocalMeshField(domid, pres_hyd)
    end if
 
    if (present(lcmesh3d)) then
      call mesh%GetLocalMesh( domid, lcmesh )
      nullify( lcmesh3d )
 
      select type(lcmesh)
      type is (localmesh3d)
        if (present(lcmesh3d)) lcmesh3d => lcmesh
      end select
    end if
 
    return
  end subroutine atmosvars_getlocalmeshprgvar
 
!OCL SERIAL
  subroutine atmosvars_getlocalmeshprgvars( domID, mesh, prgvars_list, auxvars_list, &
    DDENS, MOMX, MOMY, MOMZ, THERM,                                                  &
    DENS_hyd, PRES_hyd, Rtot, CVtot, CPtot,                                          &
    lcmesh3D                                                                         )
    implicit none
    integer, intent(in) :: domid
    class(meshbase), intent(in) :: mesh
    class(modelvarmanager), intent(inout) :: prgvars_list
    class(modelvarmanager), intent(inout) :: auxvars_list
    class(localmeshfieldbase), pointer, intent(out) :: ddens, momx, momy, momz, therm
    class(localmeshfieldbase), pointer, intent(out) :: dens_hyd, pres_hyd
    class(localmeshfieldbase), pointer, intent(out) :: rtot, cvtot, cptot
    class(localmesh3d), pointer, intent(out), optional :: lcmesh3d
 
    class(meshfieldbase), pointer :: field
    class(localmeshbase), pointer :: lcmesh
    !-------------------------------------------------------
 
    !--
    call prgvars_list%Get(prgvar_ddens_id, field)
    call field%GetLocalMeshField(domid, ddens)
 
    call prgvars_list%Get(prgvar_momx_id, field)
    call field%GetLocalMeshField(domid, momx)
    
    call prgvars_list%Get(prgvar_momy_id, field)
    call field%GetLocalMeshField(domid, momy)
 
    call prgvars_list%Get(prgvar_momz_id, field)
    call field%GetLocalMeshField(domid, momz)
 
    call prgvars_list%Get(prgvar_therm_id, field)
    call field%GetLocalMeshField(domid, therm)
  
    !--
    call auxvars_list%Get(auxvar_denshydro_id, field)
    call field%GetLocalMeshField(domid, dens_hyd)
 
    call auxvars_list%Get(auxvar_preshydro_id, field)
    call field%GetLocalMeshField(domid, pres_hyd)
 
    call auxvars_list%Get(auxvar_rtot_id, field)
    call field%GetLocalMeshField(domid, rtot)
 
    call auxvars_list%Get(auxvar_cvtot_id, field)
    call field%GetLocalMeshField(domid, cvtot)
 
    call auxvars_list%Get(auxvar_cptot_id, field)
    call field%GetLocalMeshField(domid, cptot)
 
    !---
    
    if ( present(lcmesh3d) ) then
      call mesh%GetLocalMesh( domid, lcmesh )
      nullify( lcmesh3d )
 
      select type(lcmesh)
      type is (localmesh3d)
        if (present(lcmesh3d)) lcmesh3d => lcmesh
      end select
    end if
 
    return
  end subroutine atmosvars_getlocalmeshprgvars
 
!OCL SERIAL
  subroutine atmosvars_getlocalmeshsfcvar( domID, mesh, auxvars2D_list, &
    PREC, PREC_ENGI, lcmesh2D                                           )
 
   implicit none
   integer, intent(in) :: domid
   class(meshbase), intent(in) :: mesh
   class(modelvarmanager), intent(inout) :: auxvars2d_list
   class(localmeshfieldbase), pointer, intent(out) :: prec, prec_engi
   class(localmesh2d), pointer, intent(out), optional :: lcmesh2d
 
   class(meshfieldbase), pointer :: field
   class(localmeshbase), pointer :: lcmesh
   !-------------------------------------------------------
 
   !--
   call auxvars2d_list%Get(atmos_auxvars2d_prec_id, field)
   call field%GetLocalMeshField(domid, prec)
 
   call auxvars2d_list%Get(atmos_auxvars2d_prec_engi_id, field)
   call field%GetLocalMeshField(domid, prec_engi)
 
   if (present(lcmesh2d)) then
     call mesh%GetLocalMesh( domid, lcmesh )
     nullify( lcmesh2d )
 
     select type(lcmesh)
     type is (localmesh2d)
       if (present(lcmesh2d)) lcmesh2d => lcmesh
     end select
   end if
 
   return
 end subroutine atmosvars_getlocalmeshsfcvar
 
!OCL SERIAL
  subroutine atmosvars_getlocalmeshqtrcvar( domID, mesh, trcvars_list,  &
    varid,                                                              &
    var, lcmesh3D                                                       )
 
   implicit none
   integer, intent(in) :: domid
   class(meshbase), intent(in) :: mesh
   class(modelvarmanager), intent(inout) :: trcvars_list
   integer, intent(in) :: varid
   class(localmeshfieldbase), pointer, intent(out) :: var
   class(localmesh3d), pointer, intent(out), optional :: lcmesh3d
 
   class(meshfieldbase), pointer :: field
   class(localmeshbase), pointer :: lcmesh
   !-------------------------------------------------------
 
   !--
   call trcvars_list%Get(varid, field)
   call field%GetLocalMeshField(domid, var)
 
   if (present(lcmesh3d)) then
     call mesh%GetLocalMesh( domid, lcmesh )
     nullify( lcmesh3d )
 
     select type(lcmesh)
     type is (localmesh3d)
       if (present(lcmesh3d)) lcmesh3d => lcmesh
     end select
   end if
 
   return
  end subroutine atmosvars_getlocalmeshqtrcvar
 
!OCL SERIAL
  subroutine atmosvars_getlocalmeshqtrc_qv( domID, mesh, trcvars_list, forcing_list, &
    var, var_tp, lcmesh3D                                               )
 
    use scale_atmos_hydrometeor, only: &
      atmos_hydrometeor_dry, &
      i_qv    
   implicit none
   integer, intent(in) :: domid
   class(meshbase), intent(in) :: mesh
   class(modelvarmanager), intent(inout) :: trcvars_list
   class(modelvarmanager), intent(inout) :: forcing_list
   class(localmeshfieldbase), pointer, intent(out) :: var
   class(localmeshfieldbase), pointer, intent(out), optional :: var_tp
   class(localmesh3d), pointer, intent(out), optional :: lcmesh3d
 
   class(meshfieldbase), pointer :: field
   class(localmeshbase), pointer :: lcmesh
 
   integer :: iq, tend_iq
   !-------------------------------------------------------
 
   !--
 
   if ( atmos_hydrometeor_dry ) then
     iq = 0; tend_iq = phytend_num1+1
   else
     iq = i_qv; tend_iq = phytend_num1 + i_qv
   end if
 
   call trcvars_list%Get(iq, field)
   call field%GetLocalMeshField(domid, var)
 
   call forcing_list%Get(tend_iq, field)
   if (present(var_tp)) call field%GetLocalMeshField(domid, var_tp)
 
   if (present(lcmesh3d)) then
     call mesh%GetLocalMesh( domid, lcmesh )
     nullify( lcmesh3d )
 
     select type(lcmesh)
     type is (localmesh3d)
       if (present(lcmesh3d)) lcmesh3d => lcmesh
     end select
   end if
 
   return
  end subroutine atmosvars_getlocalmeshqtrc_qv
 
!OCL SERIAL
  subroutine atmosvars_getlocalmeshqtrcvarlist( domID, mesh, trcvars_list,  &
    varid_s,                                                                &
    var_list, lcmesh3D                                                      )
 
    implicit none
    integer, intent(in) :: domid
    class(meshbase), intent(in) :: mesh
    class(modelvarmanager), intent(inout) :: trcvars_list
    integer, intent(in) :: varid_s
    type(localmeshfieldbaselist), intent(out) :: var_list(:)
    class(localmesh3d), pointer, intent(out), optional :: lcmesh3d
 
    class(meshfieldbase), pointer :: field
    class(localmeshbase), pointer :: lcmesh
 
    integer :: iq
    !-------------------------------------------------------
 
    !--
    do iq = varid_s, varid_s + size(var_list) - 1
      call trcvars_list%Get(iq, field)
      call field%GetLocalMeshField(domid, var_list(iq-varid_s+1)%ptr)
    end do
    if (present(lcmesh3d)) then
      call mesh%GetLocalMesh( domid, lcmesh )
      nullify( lcmesh3d )
 
      select type(lcmesh)
      type is (localmesh3d)
        if (present(lcmesh3d)) lcmesh3d => lcmesh
      end select
    end if
 
    return
  end subroutine atmosvars_getlocalmeshqtrcvarlist
 
!OCL SERIAL
  subroutine atmosvars_getlocalmeshphyauxvars( domID, mesh, phyauxvars_list, &
    PRES, PT,                                                                &
    lcmesh3D                                                                 )
 
    implicit none
    integer, intent(in) :: domid
    class(meshbase), intent(in) :: mesh
    class(modelvarmanager), intent(inout) :: phyauxvars_list
    class(localmeshfieldbase), pointer, intent(out) :: pres, pt
    class(localmesh3d), pointer, intent(out), optional :: lcmesh3d
 
    class(meshfieldbase), pointer :: field
    class(localmeshbase), pointer :: lcmesh
    !-------------------------------------------------------
 
    !--    
    call phyauxvars_list%Get(auxvar_pres_id, field)
    call field%GetLocalMeshField(domid, pres)
 
    call phyauxvars_list%Get(auxvar_pt_id, field)
    call field%GetLocalMeshField(domid, pt)
 
    !---
    
    if (present(lcmesh3d)) then
      call mesh%GetLocalMesh( domid, lcmesh )
      nullify( lcmesh3d )
 
      select type(lcmesh)
      type is (localmesh3d)
        if (present(lcmesh3d)) lcmesh3d => lcmesh
      end select
    end if
 
    return
  end subroutine atmosvars_getlocalmeshphyauxvars
 
!OCL SERIAL
  subroutine atmosvars_getlocalmeshphytends( domID, mesh, phytends_list,  &
    DENS_tp, MOMX_tp, MOMY_tp, MOMZ_tp, RHOT_tp, RHOH_p,                  &
    RHOQ_tp,                                                              &
    lcmesh3D                                                              )
 
    implicit none
    integer, intent(in) :: domid
    class(meshbase), intent(in) :: mesh
    class(modelvarmanager), intent(inout) :: phytends_list
    class(localmeshfieldbase), pointer, intent(out) :: dens_tp, momx_tp, momy_tp, momz_tp, rhot_tp
    class(localmeshfieldbase), pointer, intent(out) :: rhoh_p
    type(localmeshfieldbaselist), intent(inout), optional :: rhoq_tp(qa)
    class(localmesh3d), pointer, intent(out), optional :: lcmesh3d
 
    class(meshfieldbase), pointer :: field
    class(localmeshbase), pointer :: lcmesh
 
    integer :: iq
    !-------------------------------------------------------
 
    !--
    call phytends_list%Get(phytend_dens_id, field)
    call field%GetLocalMeshField(domid, dens_tp)
 
    call phytends_list%Get(phytend_momx_id, field)
    call field%GetLocalMeshField(domid, momx_tp)
    
    call phytends_list%Get(phytend_momy_id, field)
    call field%GetLocalMeshField(domid, momy_tp)
 
    call phytends_list%Get(phytend_momz_id, field)
    call field%GetLocalMeshField(domid, momz_tp)
 
    call phytends_list%Get(phytend_rhot_id, field)
    call field%GetLocalMeshField(domid, rhot_tp)
 
    call phytends_list%Get(phytend_rhoh_id, field)
    call field%GetLocalMeshField(domid, rhoh_p)
 
    if ( present(rhoq_tp) ) then
      do iq = 1, qa
        call phytends_list%Get(phytend_num1+iq, field)
        call field%GetLocalMeshField(domid, rhoq_tp(iq)%ptr)  
      end do
    end if
 
    !---
    if ( present(lcmesh3d) ) then
      call mesh%GetLocalMesh( domid, lcmesh )
      nullify( lcmesh3d )
 
      select type(lcmesh)
      type is (localmesh3d)
        if ( present(lcmesh3d) ) lcmesh3d => lcmesh
      end select
    end if
 
    return
  end subroutine atmosvars_getlocalmeshphytends
 
!OCL SERIAL
  subroutine atmosvars_getlocalmeshqtrcphytend( domID, mesh, phytends_list,  &
    qtrcid,                                                                  &
    RHOQ_tp                                                                  )
 
    implicit none
    integer, intent(in) :: domid
    class(meshbase), intent(in) :: mesh
    class(modelvarmanager), intent(inout) :: phytends_list
    integer, intent(in) :: qtrcid
    class(localmeshfieldbase), pointer, intent(out) :: rhoq_tp
 
    class(meshfieldbase), pointer :: field
    !-------------------------------------------------------
 
    call phytends_list%Get(phytend_num1 + qtrcid, field)
    call field%GetLocalMeshField(domid, rhoq_tp)
 
    return
  end subroutine atmosvars_getlocalmeshqtrcphytend  
  
!-- private -----------------------------------------------------------------------
!OCL SERIAL
  subroutine setup_phys_preoperation( this, phy_preproc_file_basename, mesh3D, elem3D )
    use scale_prc, only: prc_ismaster
    implicit none
    class(atmosvarscontainer), intent(inout) :: this
    character(len=*), intent(in) :: phy_preproc_file_basename
    class(meshbase3d), intent(in) :: mesh3d
    class(elementbase3d), intent(in) :: elem3d
 
    integer :: fid
    character(len=H_SHORT) :: phy_preproc_operation_type = 'None' ! ModalFilter / GlobalFilter
    real(rp) :: mf_alpha_h
    integer :: mf_order_h
    real(rp) :: mf_alpha_v
    integer :: mf_order_v
 
    character(len=H_SHORT) :: glfilteroptrtype
    character(len=H_SHORT) :: glfiltershape
    real(rp) :: glfilterwidthfac
    integer :: nnode_h1d_reconst
    namelist / param_atmos_vars_container / &
      phy_preproc_operation_type, &
      mf_alpha_h, mf_order_h, &
      mf_alpha_v, mf_order_v, &
      glfilteroptrtype, glfiltershape, glfilterwidthfac, nnode_h1d_reconst
 
    character(len=H_LONG) :: cnfname
    integer :: ierr
 
    class(physpreprocbase), pointer :: phys_pp_ptr
    !----------------------------------------
 
    log_info("ATMOS_vars_container/setup_phys_preoperation",*) "Setting up physical pre-operation for container type ", this%container_type
    
    if ( this%container_type == atm_vars_container_primary_id ) then
      this%phy_preproc_operation_type = phy_preoperation_typeid_none
      return
    end if
 
    !--
    mf_alpha_h = 0.0_rp; mf_order_h = 16
    mf_alpha_v = 0.0_rp; mf_order_v = 16
 
    glfilteroptrtype = ''
    glfiltershape    = 'GAUSSIAN'
    glfilterwidthfac = 1.0_rp
 
    nnode_h1d_reconst = -1
 
    write(cnfname,'(A,I2.2,A)') trim(phy_preproc_file_basename), this%container_type, ".conf"
    fid = io_cnf_open(trim(cnfname), prc_ismaster)
 
    rewind(fid)
    read(fid,nml=param_atmos_vars_container,iostat=ierr)
    if( ierr < 0 ) then !--- missing
      log_info("ATMOS_vars_container/setup_phys_preoperation",*) 'Not found namelist. Default used.'
    elseif( ierr > 0 ) then !--- fatal error
      log_error("ATMOS_vars_container/setup_phys_preoperation",*) 'Invalid names in namelist PARAM_ATMOS_VARS_CONTAINER. Check!'
      call prc_abort
    endif
    log_nml(param_atmos_vars_container)
    close(fid)
 
    !--
    select case( phy_preproc_operation_type ) 
    case( 'None' )
      this%phy_preproc_operation_type = phy_preoperation_typeid_none
      allocate( physpreprocnone :: this%phy_preproc )
    case( 'ModalFilter' )
      this%phy_preproc_operation_type = phy_preoperation_typeid_modalfilter 
      allocate( physpreprocmodalfilter :: this%phy_preproc )
    case( 'GlobalFilter' )
      this%phy_preproc_operation_type = phy_preoperation_typeid_globalfilter
      allocate( physpreprocglobalfilter :: this%phy_preproc )
    case default
      log_error("ATMOS_vars_container/setup_phys_preoperation",*) 'Unsupported PHY_PREPROC_OPERATION_TYPE is specified. Check!', phy_preproc_operation_type
      call prc_abort
    end select
 
    select type( phys_pp_ptr => this%phy_preproc )
    class is (physpreprocnone)
      call phys_pp_ptr%Init( mesh3d )
    class is (physpreprocmodalfilter)
      call phys_pp_ptr%Init( mf_alpha_h, mf_order_h, mf_alpha_v, mf_order_v, mesh3d )
    class is (physpreprocglobalfilter)
      if ( nnode_h1d_reconst < 0 ) nnode_h1d_reconst = elem3d%Nnode_h1D
      call phys_pp_ptr%Init( glfilteroptrtype, &
        glfiltershape, glfilterwidthfac, &
        nnode_h1d_reconst,               &
        mesh3d )
    end select
 
    return
  end subroutine setup_phys_preoperation
 
!OCL SERIAL
  subroutine vars_calc_diagnosevar_lc( field_name, var_out,  &
    DDENS_, MOMX_, MOMY_, MOMZ_, PRES_, QDRY_, QTRC,         &
    DENS_hyd, PRES_hyd, Rtot, CVtot, CPTot,                  &
    lcmesh, elem )
 
    use scale_const, only: &
      grav => const_grav,   &
      rdry => const_rdry,   &
      rvap => const_rvap,   &
      cpdry => const_cpdry, &
      cvdry => const_cvdry, &
      pres00 => const_pre00
    use scale_tracer, only: &
      tracer_inq_id,          &
      tracer_cv, tracer_engi0
    use scale_atmos_hydrometeor, only: &
      atmos_hydrometeor_dry      
    use scale_atmos_saturation, only: &
      atmos_saturation_psat_liq
    implicit none
 
    class(localmesh3d), intent(in) :: lcmesh
    class(elementbase3d), intent(in) :: elem
    character(*), intent(in) :: field_name
    real(rp), intent(out) :: var_out(elem%np,lcmesh%nea)
    real(rp), intent(in) :: ddens_(elem%np,lcmesh%nea)
    real(rp), intent(in) :: momx_(elem%np,lcmesh%nea)
    real(rp), intent(in) :: momy_(elem%np,lcmesh%nea)
    real(rp), intent(in) :: momz_(elem%np,lcmesh%nea)
    real(rp), intent(in) :: pres_(elem%np,lcmesh%nea)
    real(rp), intent(in) :: qdry_(elem%np,lcmesh%nea)
    type(localmeshfieldbaselist), intent(in) :: qtrc(qa)
    real(rp), intent(in) :: dens_hyd(elem%np,lcmesh%nea)
    real(rp), intent(in) :: pres_hyd(elem%np,lcmesh%nea)
    real(rp), intent(in) :: rtot (elem%np,lcmesh%nea)
    real(rp), intent(in) :: cvtot(elem%np,lcmesh%nea)
    real(rp), intent(in) :: cptot(elem%np,lcmesh%nea)
 
    integer :: ke, ke2d
    integer :: p
    integer :: iq
    real(rp) :: dens
    real(rp) :: mom_u1, mom_u2, g_11, g_12, g_22
    real(rp) :: temp(elem%np), psat(elem%np)
 
    integer :: iq_qv
    integer :: indexh2dto3d(elem%np)
 
    integer :: ne, np
    !-------------------------------------------------------------------------
 
    ne = lcmesh%Ne; np = elem%Np
 
    select case(trim(field_name))
    case('DENS')
      !$omp parallel do
      !$acc parallel loop collapse(2) present(DDENS_, DENS_hyd, var_out) async(1)
      do ke=1, ne
      do p=1, np
        var_out(p,ke) = ddens_(p,ke) + dens_hyd(p,ke)
      end do
      end do
    
    case('U')
      !$omp parallel do private (DENS)
      !$acc parallel loop collapse(2) present(DDENS_, DENS_hyd, MOMX_, var_out) async(1)
      do ke=1, ne
      do p=1, np
        dens = ddens_(p,ke) + dens_hyd(p,ke)
        var_out(p,ke) = momx_(p,ke) / dens
      end do
      end do
    
    case('V')
      !$omp parallel do private (DENS)
      !$acc parallel loop collapse(2) present(DDENS_, DENS_hyd, MOMY_, var_out) async(1)
      do ke=1, ne
      do p=1, np
        dens = ddens_(p,ke) + dens_hyd(p,ke)
        var_out(p,ke) = momy_(p,ke) / dens
      end do
      end do  
          
    case('W')
      !$omp parallel do private (DENS)
      !$acc parallel loop collapse(2) present(DDENS_, DENS_hyd, MOMZ_, var_out) async(1)
      do ke=1, ne
      do p=1, np
        dens = ddens_(p,ke) + dens_hyd(p,ke)
        var_out(p,ke) = momz_(p,ke) / dens
      end do
      end do
    
    case ( 'PRES' )
    case('PRES_diff')  
      !$omp parallel do
      !$acc parallel loop collapse(2) present(PRES_, PRES_hyd, var_out) async(1)
      do ke=1, ne
      do p=1, np
        var_out(p,ke) = pres_(p,ke) - pres_hyd(p,ke)
      end do
      end do
    
    case('T')
      !$omp parallel do
      !$acc parallel loop collapse(2) present(PRES_, Rtot, DDENS_, DENS_hyd, var_out) async(1)
      do ke=1, ne
      do p=1, np
        var_out(p,ke) = pres_(p,ke) / (rtot(p,ke) * (ddens_(p,ke) + dens_hyd(p,ke)) )
      end do
      end do
    
    case('T_diff')
      !$omp parallel do
      !$acc parallel loop collapse(2) present(PRES_, Rtot, DDENS_, DENS_hyd, PRES_hyd, var_out) async(1)
      do ke=1, ne
      do p=1, np
        var_out(p,ke) = pres_(p,ke) / ( rtot(p,ke) * (ddens_(p,ke) + dens_hyd(p,ke)) ) &
                      - pres_hyd(p,ke) / ( rdry * dens_hyd(p,ke) )
      end do
      end do
    
    case('PT')
      !$omp parallel do private(DENS)
      !$acc parallel loop collapse(2) present(PRES_, Rtot, CVtot, CPtot, DDENS_, DENS_hyd, var_out) async(1)
      do ke=1, ne
      do p=1, np
        dens = ddens_(p,ke) + dens_hyd(p,ke)
        var_out(p,ke) = pres_(p,ke) / (rtot(p,ke) * dens ) * ( pres00 / pres_(p,ke) )**( rtot(p,ke) / cptot(p,ke) )
      end do 
      end do
    
    case('PT_diff')
      !$omp parallel do private( DENS )
      !$acc parallel loop collapse(2) present(PRES_, Rtot, CVtot, CPtot, DDENS_, DENS_hyd, PRES_hyd, var_out) async(1)
      do ke=1, ne
      do p=1, np
        dens = ddens_(p,ke) + dens_hyd(p,ke)
        var_out(p,ke) = pres_(p,ke) / (rtot(p,ke) * dens ) * ( pres00 / pres_(p,ke) )**( rtot(p,ke) / cptot(p,ke) ) &
                      - pres00/rdry * (pres_hyd(p,ke)/pres00)**(cvdry/cpdry) / dens_hyd(p,ke)
      end do 
      end do
    
    case( 'RH', 'RHL' )
      if ( atmos_hydrometeor_dry ) then
        !$omp parallel do
        !$acc parallel loop collapse(2) present(var_out) async(1)
        do ke=1, ne
        do p=1, np
          var_out(p,ke) = 0.0_rp
        end do
        end do
      else
        call tracer_inq_id( "QV", iq_qv )
 
        !$omp parallel do private (TEMP, PSAT)
        !$acc parallel loop gang private(TEMP, PSAT) present(PRES_, Rtot, DDENS_, DENS_hyd, var_out) async(1)
        do ke=1, ne
          !$acc loop vector
          do p=1, np
            temp(p) = pres_(p,ke) / (rtot(p,ke) * (ddens_(p,ke) + dens_hyd(p,ke)) )
#ifdef _OPENACC
            call atmos_saturation_psat_liq( temp(p), & ! (in)
              psat(p)                                ) ! (out)
#endif
          end do
 
#ifndef _OPENACC
          call atmos_saturation_psat_liq( &
            np, 1, np, temp(:),     & ! (in)
            psat(:)                           ) ! (out)
#endif
 
          !$acc loop vector
          do p=1, np
            var_out(p,ke) = ( ddens_(p,ke) + dens_hyd(p,ke) ) * qtrc(iq_qv)%ptr%val(p,ke) &
                          / psat(p) * rvap * temp(p) * 100.0_rp
          end do
        end do 
      end if
    
    case('ENGK')
      indexh2dto3d(:) = elem%IndexH2Dto3D(:)
 
      !$omp parallel do private (ke2D, DENS, mom_u1, mom_u2, G_11, G_12, G_22)
      !$acc parallel loop collapse(2) present(DDENS_, DENS_hyd, MOMX_, MOMY_, MOMZ_, var_out, lcmesh,elem) copyin(IndexH2Dto3D) async(1)
      do ke=1, ne
      do p=1, np
        ke2d = lcmesh%EMap3Dto2D(ke)
        
        dens = ddens_(p,ke) + dens_hyd(p,ke)
        g_11 = lcmesh%G_ij(indexh2dto3d(p),ke2d,1,1)
        g_12 = lcmesh%G_ij(indexh2dto3d(p),ke2d,1,2)
        g_22 = lcmesh%G_ij(indexh2dto3d(p),ke2d,2,2)
 
        mom_u1 = g_11 * momx_(p,ke) + g_12 * momy_(p,ke)
        mom_u2 = g_12 * momx_(p,ke) + g_22 * momy_(p,ke)
 
        var_out(p,ke) = 0.5_rp * ( momx_(p,ke) * mom_u1 + momy_(p,ke) * mom_u2 + momz_(p,ke)**2 ) / dens
      end do
      end do
    
    case('ENGP')
      !$omp parallel do private (DENS)
      !$acc parallel loop collapse(2) present(DDENS_, DENS_hyd, var_out, lcmesh) async(1)
      do ke=1, ne
      do p=1, np
        dens = ddens_(p,ke) + dens_hyd(p,ke)
        var_out(p,ke) = dens * grav * lcmesh%zlev(p,ke)
      end do
      end do
    
    case('ENGI')
      !$omp parallel do private (DENS, iq)
      !$acc parallel loop collapse(2) present(DDENS_, DENS_hyd, PRES_, Rtot, QDRY_, var_out) async(1)
      do ke=1, ne
      do p=1, np
        dens = ddens_(p,ke) + dens_hyd(p,ke)
        var_out(p,ke) = qdry_(p,ke) * pres_(p,ke) / rtot(p,ke) * cvdry
        !$acc loop seq
        do iq = 1, qa
          var_out(p,ke) = var_out(p,ke) &
            + qtrc(iq)%ptr%val(p,ke) * ( pres_(p,ke) / rtot(p,ke) * tracer_cv(iq) + dens * tracer_engi0(iq) )
        end do
      end do
      end do
    
    case('ENGT')
      indexh2dto3d(:) = elem%IndexH2Dto3D(:)
 
      !$omp parallel do private (ke2D, DENS, mom_u1, mom_u2, iq, G_11, G_12, G_22)
      !$acc parallel loop collapse(2) present(DDENS_, DENS_hyd, MOMX_, MOMY_, MOMZ_, PRES_, Rtot, QDRY_, var_out, lcmesh, elem) copyin(IndexH2Dto3D) async(1)
      do ke=1, ne
      do p=1, np
        ke2d = lcmesh%EMap3Dto2D(ke)
 
        dens = ddens_(p,ke) + dens_hyd(p,ke)
 
        g_11 = lcmesh%G_ij(indexh2dto3d(p),ke2d,1,1)
        g_12 = lcmesh%G_ij(indexh2dto3d(p),ke2d,1,2)
        g_22 = lcmesh%G_ij(indexh2dto3d(p),ke2d,2,2)
        mom_u1 = g_11 * momx_(p,ke) + g_12 * momy_(p,ke)
        mom_u2 = g_12 * momx_(p,ke) + g_22 * momy_(p,ke)
 
        ! ENGI
        var_out(p,ke) = qdry_(p,ke) * pres_(p,ke) / rtot(p,ke) * cvdry
        !$acc loop seq
        do iq = 1, qa
          var_out(p,ke) = var_out(p,ke) &
            + qtrc(iq)%ptr%val(p,ke) * ( pres_(p,ke) / rtot(p,ke) * tracer_cv(iq) + dens * tracer_engi0(iq) )
        end do
        ! ENGT
        var_out(p,ke) = &
            0.5_rp * ( momx_(p,ke) * mom_u1 + momy_(p,ke) * mom_u2 + momz_(p,ke)**2 ) / dens & ! ENGK       
          + var_out(p,ke)                                                                    & ! ENGI
          + dens * grav * lcmesh%pos_en(p,ke,3)                                                ! ENGP
      end do
      end do
    
    case default
      log_error("AtmosVars_calc_diagnoseVar_lc",*) 'The name of diagnostic variable is not suported. Check!', field_name
      call prc_abort
    
    end select
 
    return
  end subroutine vars_calc_diagnosevar_lc
end module mod_atmos_vars_container