mod_atmos_vars.F90

Go to the documentation of this file.

!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------
#include "scaleFElib.h"
module mod_atmos_vars
  !-----------------------------------------------------------------------------
  !
  !++ Used modules
  !
  use scale_precision
  use scale_io
  use scale_prc
  use scale_debug
  use scale_tracer, only: &
    qa, tracer_name, tracer_desc, tracer_unit
  use scale_atmos_hydrometeor, only: &
    atmos_hydrometeor_dry
 
  use scale_element_base, only: &
    elementbase, elementbase2d, elementbase3d
  use scale_mesh_base, only: meshbase
  use scale_mesh_base2d, only: &
    meshbase2d
  use scale_mesh_base3d, only: &
    meshbase3d,                              &
    dimtype_xyz  => meshbase3d_dimtypeid_xyz
 
  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_meshfield_base, only: &
    meshfieldbase, meshfield2d, meshfield3d, meshfield3dlist
  
  use scale_file_restart_meshfield, only: &
    file_restart_meshfield_component
  
  use scale_meshfieldcomm_base, only: meshfieldcontainer
 
  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_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
    
  !-----------------------------------------------------------------------------
  implicit none
  private
 
  !-----------------------------------------------------------------------------
  !
  !++ Public type & procedures
  !
  type, public :: atmosvars
    class(atmosmesh), pointer :: mesh
 
    !- prognostic variables
    type(meshfield3d), allocatable :: prog_vars(:)
    type(modelvarmanager) :: progvars_manager
    integer :: prog_vars_commid
 
    !- tracer variables
    type(meshfield3d), allocatable :: qtrc_vars(:)
    type(modelvarmanager) :: qtrcvars_manager
    integer :: qtrc_vars_commid 
 
    !- auxiliary variables    
    type(meshfield3d), allocatable :: aux_vars(:)
    type(modelvarmanager) :: auxvars_manager 
    integer :: aux_vars_commid
 
    !- auxiliary variables (2D)
    type(meshfield2d), allocatable :: aux_vars2d(:)
    type(modelvarmanager) :: auxvars2d_manager 
    
    !-
    type(modelvarmanager), pointer :: ptr_mp_auxvars2d_manager
 
    logical :: moist
    type(meshfield3d), pointer :: qv
    type(meshfield3d) :: zero
 
    !- Tendency with physics
    type(meshfield3d), allocatable :: phy_tend(:)
    type(modelvarmanager) :: phytends_manager 
    integer :: phytends_commid
    integer :: phytend_num_tot
 
    !--
    integer, allocatable :: diagvars2d_histid(:)
    integer, allocatable :: diagvars3d_histid(:)
 
    type(file_restart_meshfield_component) :: restart_file
    
    logical :: check_range
    logical :: check_total
 
  contains
    procedure :: init => atmosvars_init
    procedure :: final => atmosvars_final
    procedure :: calc_diagnostics => atmosvars_calculatediagnostics
    procedure :: calc_diagvar => atmosvars_calcdiagvar
    procedure :: calc_diagvar2d => atmosvars_calcdiagvar2d
    procedure :: history => atmosvars_history
    procedure :: check   => atmosvars_check
    procedure :: monitor => atmosvars_monitor
    procedure :: read_restart_file => atmosvar_read_restart_file
    procedure :: write_restart_file => atmosvar_write_restart_file
    procedure :: regist_physvar_manager => atmosvars_regist_physvar_manager
  end type atmosvars
 
  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
  !
  real(rp), parameter :: progvars_check_min(prgvar_num) = (/ -1.0_rp, -100.0_rp, -200.0_rp, -200.0_rp, -200.0_rp /)
  real(rp), parameter :: progvars_check_max(prgvar_num) = (/  1.0_rp,   100.0_rp, 200.0_rp,  200.0_rp,  200.0_rp /)
              
 
  ! 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', ''  )                    /
  
                    
 
  ! Diagnostic variables
 
  integer, public, parameter :: atmos_diagvars_dens_id   = 1
  integer, public, parameter :: atmos_diagvars_u_id      = 2
  integer, public, parameter :: atmos_diagvars_v_id      = 3
  integer, public, parameter :: atmos_diagvars_w_id      = 4
  integer, public, parameter :: atmos_diagvars_t_id      = 5  
  integer, public, parameter :: atmos_diagvars_umet_id   = 6
  integer, public, parameter :: atmos_diagvars_vmet_id   = 7
  integer, public, parameter :: atmos_diagvars_qdry_id   = 8
  integer, public, parameter :: atmos_diagvars_rh_id     = 9
  integer, public, parameter :: atmos_diagvars_engt_id   = 10
  integer, public, parameter :: atmos_diagvars_engp_id   = 11
  integer, public, parameter :: atmos_diagvars_engk_id   = 12
  integer, public, parameter :: atmos_diagvars_engi_id   = 13
  integer, public, parameter :: atmos_diagvars3d_num     = 13
 
  type(variableinfo), public :: atmos_diagvars3d_vinfo(atmos_diagvars3d_num)
  DATA atmos_diagvars3d_vinfo / &
    variableinfo( atmos_diagvars_dens_id, 'DENS',  'density'              , 'kg/m3', 3, 'XYZ', 'air_density'          ), &
    variableinfo( atmos_diagvars_u_id   , 'U'   ,  'velocity u'           , 'm/s'  , 3, 'XYZ', 'x_wind'               ), &
    variableinfo( atmos_diagvars_v_id   , 'V'   ,  'velocity v'           , 'm/s'  , 3, 'XYZ', 'y_wind'               ), &  
    variableinfo( atmos_diagvars_w_id   , 'W'   ,  'velocity w'           , 'm/s'  , 3, 'XYZ', 'upward_air_velocity'  ), &
    variableinfo( atmos_diagvars_t_id   , 'T'   ,  'temperature'          , 'K'    , 3, 'XYZ', 'air_temperature'      ), &
    variableinfo( atmos_diagvars_umet_id, 'Umet',  'eastward velocity'    , 'm/s'  , 3, 'XYZ', 'x_wind'               ), &
    variableinfo( atmos_diagvars_vmet_id, 'Vmet',  'northward velocity'   , 'm/s'  , 3, 'XYZ', 'y_wind'               ), &  
    variableinfo( atmos_diagvars_qdry_id, 'QDRY',  'dry air'              , 'kg/kg', 3, 'XYZ', ''                     ), &
    variableinfo( atmos_diagvars_rh_id  ,   'RH', 'relative humidity(liq)',     '%', 3, 'XYZ', 'relative_humidity'    ), &
    variableinfo( atmos_diagvars_engt_id, 'ENGT',  'total energy'         , 'J/m3' , 3, 'XYZ', ''                     ), &
    variableinfo( atmos_diagvars_engp_id, 'ENGP',  'potential energy'     , 'J/m3' , 3, 'XYZ', ''                     ), &
    variableinfo( atmos_diagvars_engk_id, 'ENGK',  'kinetic energy'       , 'J/m3' , 3, 'XYZ', ''                     ), &
    variableinfo( atmos_diagvars_engi_id, 'ENGI',   'internal energy'     , 'J/m3' , 3, 'XYZ', ''                     )  /
 
 
  integer, public, parameter :: atmos_diagvars_rain_id   = 1
  integer, public, parameter :: atmos_diagvars_snow_id   = 2
  integer, public, parameter :: atmos_diagvars2d_num     = 2
  
  type(variableinfo), public :: atmos_diagvars2d_vinfo(atmos_diagvars2d_num)
  DATA atmos_diagvars2d_vinfo / &
    variableinfo( atmos_diagvars_rain_id, 'RAIN', 'surface rain flux'              , 'kg/m2/s', 2, 'XY', 'rainfall_flux'       ), &
    variableinfo( atmos_diagvars_snow_id, 'SNOW', 'surface snow flux'              , 'kg/m2/s', 2, 'XY', 'snowfall_flux'       )  /
 
  !-----------------------------------------------------------------------------
  !
  !++ Private procedures & variables
  !
  !-------------------
 
  private :: vars_calc_diagnosevar_lc
  private :: vars_calc_diagnosevar2d_lc
 
  ! for monitor
 
  integer, private, parameter   :: im_qdry         =  1
  integer, private, parameter   :: im_qtot         =  2
  integer, private, parameter   :: im_engt         =  3
  integer, private, parameter   :: im_engp         =  4
  integer, private, parameter   :: im_engk         =  5
  integer, private, parameter   :: im_engi         =  6
  integer, private, parameter   :: dvm_nmax        =  6
  integer, private              :: dv_monit_id(dvm_nmax)
 
contains
 
!OCL SERIAL
  subroutine atmosvars_init( this, atm_mesh )
    use scale_atmos_hydrometeor, only: &
      atmos_hydrometeor_dry
    use scale_file_monitor_meshfield, only:    &
      monitor_reg => file_monitor_meshfield_reg
 
    use scale_atm_dyn_dgm_nonhydro3d_common, only: &
      prgvar_scalar_num, prgvar_hvec_num,          &
      atm_dyn_dgm_nonhydro3d_common_setup_variables
    implicit none
 
    class(atmosvars), target, intent(inout) :: this
    class(atmosmesh), target, intent(inout) :: atm_mesh
 
    integer :: n
    integer :: iv
    integer :: iq
    logical :: reg_file_hist
 
    type(modelvarmanager) :: diagvar_manager               ! dummy
    type(meshfield2d) :: diag_vars2d(atmos_diagvars2d_num) ! dummy
    type(meshfield3d) :: diag_vars3d(atmos_diagvars3d_num) ! dummy
 
    logical :: check_range    = .false.
    logical :: check_total    = .false.
 
    namelist / param_atmos_vars / &
      check_range, &
      check_total
 
    character(len=H_LONG) :: in_basename           = ''
    logical :: in_postfix_timelabel                = .false.   
    character(len=H_LONG) :: out_basename          = ''
    logical :: out_postfix_timelabel               = .true.    
    character(len=H_MID) :: out_title              = ''
    character(len=H_SHORT) :: out_dtype            = 'DEFAULT'
 
    namelist / param_atmos_vars_restart / &
      in_basename,           &
      in_postfix_timelabel,  &
      out_basename,          &
      out_postfix_timelabel, &
      out_title,             &
      out_dtype    
    
    integer :: ierr
    logical :: is_specified
 
    class(meshbase3d), pointer :: mesh3d
    class(meshbase2d), pointer :: mesh2d
 
    type(variableinfo) :: prgvar_info(prgvar_num)
    !--------------------------------------------------
 
    log_info('AtmosVars_Init',*)
 
    !- read namelist
    rewind(io_fid_conf)
    read(io_fid_conf,nml=param_atmos_vars,iostat=ierr)
    if( ierr < 0 ) then !--- missing
       log_info("ATMOS_vars_setup",*) 'Not found namelist. Default used.'
    elseif( ierr > 0 ) then !--- fatal error
       log_error("ATMOS_vars_setup",*) 'Invalid names in namelist PARAM_ATMOS_VARS. Check!'
       call prc_abort
    endif
    log_nml(param_atmos_vars)
 
    !- Set the pointer of mesh
    this%mesh => atm_mesh
    mesh3d => atm_mesh%ptr_mesh
    call mesh3d%GetMesh2D( mesh2d )
 
    !- Initialize variables associated with dynamical core 
    !  (prognostic variables, tracer variables, 3D auxiliary variables, and tendencies of physical processes)
 
    call this%PROGVARS_manager%Init()
    call this%QTRCVARS_manager%Init()
    call this%AUXVARS_manager%Init()
    call this%PHYTENDS_manager%Init()
 
    allocate( this%PROG_VARS(prgvar_num) )
    allocate( this%QTRC_VARS(0:qa) )
    allocate( this%AUX_VARS(auxvar_num) )
 
    this%PHYTEND_NUM_TOT = phytend_num1 + max(1,qa)
    allocate( this%PHY_TEND(this%PHYTEND_NUM_TOT) )
 
    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)
      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 ( .not. atmos_hydrometeor_dry ) 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
 
    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
 
 
    !- Initialize 2D auxiliary variables
    call this%AUXVARS2D_manager%Init()
    allocate( this%AUX_VARS2D(atmos_auxvars2d_num) )
    
    reg_file_hist = .true.
    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, fill_zero=.true.      ) ! (in)
    end do
  
 
    !- Initialize diagnostic variables for output
 
    ! 3D
    call diagvar_manager%Init()
    allocate( this%DIAGVARS3D_HISTID(atmos_diagvars3d_num) )
 
    reg_file_hist = .true.
    do iv = 1, atmos_diagvars3d_num
      call diagvar_manager%Regist( &
        atmos_diagvars3d_vinfo(iv), atm_mesh%ptr_mesh, & ! (in) 
        diag_vars3d(iv), reg_file_hist                   ) ! (out)
      
      this%DIAGVARS3D_HISTID(iv) = diag_vars3d(iv)%hist_id
    end do
    call diagvar_manager%Final()
 
    ! 2D
    call diagvar_manager%Init()
    allocate( this%DIAGVARS2D_HISTID(atmos_diagvars2d_num) )
 
    reg_file_hist = .true.
    do iv = 1, atmos_diagvars2d_num
      call diagvar_manager%Regist( &
        atmos_diagvars2d_vinfo(iv), mesh2d,  & ! (in) 
        diag_vars2d(iv), reg_file_hist       ) ! (out)
      
      this%DIAGVARS2D_HISTID(iv) = diag_vars2d(iv)%hist_id
    end do
    call diagvar_manager%Final()
 
    !-- Setup information for input/output restart files. 
 
    is_specified = .true.
    !- read namelist
    rewind(io_fid_conf)
    read(io_fid_conf,nml=param_atmos_vars_restart,iostat=ierr)
    if( ierr < 0 ) then !--- missing
       log_info("AtmosVars_Init",*) 'Not found namelist PARAM_ATMOS_VARS_RESTART. Default used.'
       is_specified = .false.
    elseif( ierr > 0 ) then !--- fatal error
       log_error("AtmosVars_Init",*) 'Not appropriate names in namelist PARAM_ATMOS_VARS_RESTART. Check!'
       call prc_abort
    endif
    log_nml(param_atmos_vars_restart)
 
    if (is_specified) then
      call atm_mesh%Setup_restartfile( this%restart_file,           &
        in_basename, in_postfix_timelabel,                          &
        out_basename, out_postfix_timelabel, out_dtype, out_title,  &
        prgvar_num + auxvar_num                                     )
    else
      call atm_mesh%Setup_restartfile( this%restart_file, &
        prgvar_num + auxvar_num                           )
    end if
 
    !-----< monitor output setup >-----
    
    call monitor_reg( 'QTOT',         'water mass',            'kg', & ! (in)
                      dv_monit_id(im_qtot),                          & ! (out)
                      dim_type='ATM3D', is_tendency=.false.          ) ! (in)
    
    call monitor_reg( 'ENGT',         'total     energy',       'J', & ! (in)
                      dv_monit_id(im_engt),                          & ! (out)
                      dim_type='ATM3D', is_tendency=.false.          ) ! (in)
    call monitor_reg( 'ENGP',         'potential energy',       'J', & ! (in)
                      dv_monit_id(im_engp),                          & ! (out)
                      dim_type='ATM3D', is_tendency=.false.          ) ! (in)
    call monitor_reg( 'ENGK',         'kinetic   energy',       'J', & ! (in)
                      dv_monit_id(im_engk),                          & ! (out)
                      dim_type='ATM3D', is_tendency=.false.          ) ! (in)
    call monitor_reg( 'ENGI',         'internal  energy',       'J', & ! (in)
                      dv_monit_id(im_engi),                          & ! (out)
                      dim_type='ATM3D', is_tendency=.false.          ) ! (in)
 
 
    !-----< check the range of values >-----
 
    this%check_range = check_range
    this%check_total = check_total
    log_info("ATMOS_vars_setup",*) 'Check value range of variables?     : ', check_range
    log_info("ATMOS_vars_setup",*) 'Check total value of variables?     : ', check_total
 
    return
  end subroutine atmosvars_init
 
!OCL SERIAL
  subroutine atmosvars_final( this )
    implicit none
    class(atmosvars), intent(inout) :: this
 
    !--------------------------------------------------
 
    log_info('AtmosVars_Final',*)
 
    call this%restart_file%Final()
 
    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 )
 
    call this%AUXVARS2D_manager%Final()
    deallocate( this%AUX_VARS2D )
 
    call this%PHYTENDS_manager%Final()
    deallocate( this%PHY_TEND )
 
    deallocate( this%DIAGVARS3D_HISTID )
 
    return
  end subroutine atmosvars_final
 
!OCL SERIAL
  subroutine atmosvars_regist_physvar_manager( this, &
    mp_AUXVARS2D_manager )
    implicit none
 
    class(atmosvars), target, intent(inout) :: this
    type(modelvarmanager), intent(in), target :: mp_auxvars2d_manager
    !----------------------------------------------
 
    this%ptr_MP_AUXVARS2D_manager => mp_auxvars2d_manager
 
    return
  end subroutine atmosvars_regist_physvar_manager
 
!OCL SERIAL
  subroutine atmosvars_history( this )
    use scale_file_history_meshfield, only: file_history_meshfield_put
    implicit none
    class(atmosvars), intent(inout) :: this
  
    integer :: v
    integer :: hst_id
 
    type(meshfield3d) :: tmp_field3d
    class(meshbase3d), pointer :: mesh3d
 
    type(meshfield2d) :: tmp_field2d
    class(meshbase2d), pointer :: mesh2d
    !-------------------------------------------------------------------------
 
    mesh3d => this%PROG_VARS(1)%mesh
    call mesh3d%GetMesh2D(mesh2d)
 
    do v = 1, prgvar_num
      hst_id = this%PROG_VARS(v)%hist_id
      if ( hst_id > 0 ) call file_history_meshfield_put( hst_id, this%PROG_VARS(v) )
    end do
 
    do v = 1, qa
      hst_id = this%QTRC_VARS(v)%hist_id
      if ( hst_id > 0 ) call file_history_meshfield_put( hst_id, this%QTRC_VARS(v) )
    end do
 
    call this%Calc_diagnostics()
    
    do v = 1, auxvar_num
      hst_id = this%AUX_VARS(v)%hist_id
      if ( hst_id > 0 ) call file_history_meshfield_put( hst_id, this%AUX_VARS(v) )
    end do
    do v = 1, atmos_auxvars2d_num
      hst_id = this%AUX_VARS2D(v)%hist_id
      if ( hst_id > 0 ) call file_history_meshfield_put( hst_id, this%AUX_VARS2D(v) )
    end do
 
    do v = 1, this%PHYTEND_NUM_TOT
      hst_id = this%PHY_TEND(v)%hist_id
      if ( hst_id > 0 ) call file_history_meshfield_put( hst_id, this%PHY_TEND(v) )
    end do
 
    !- Output diagnostic variables
 
    ! 3D
    call tmp_field3d%Init( "tmp_field", "", mesh3d)
    do v = 1, atmos_diagvars3d_num
      hst_id = this%DIAGVARS3D_HISTID(v)
      if ( hst_id > 0 ) then
        call atmosvars_calcdiagvar( this, atmos_diagvars3d_vinfo(v)%NAME, tmp_field3d )
        call file_history_meshfield_put( hst_id, tmp_field3d )
      end if
    end do
    call tmp_field3d%Final()
 
    ! 2D
    call tmp_field2d%Init( "tmp_field", "", mesh2d)
    do v = 1, atmos_diagvars2d_num
      hst_id = this%DIAGVARS2D_HISTID(v)
      if ( hst_id > 0 ) then
        call atmosvars_calcdiagvar2d( this, atmos_diagvars2d_vinfo(v)%NAME, tmp_field2d )
        call file_history_meshfield_put( hst_id, tmp_field2d )
      end if
    end do
    call tmp_field2d%Final()
 
    return
  end subroutine atmosvars_history
 
!OCL SERIAL
  subroutine atmosvar_read_restart_file( this, atmos_mesh, dyncore  )
 
    use scale_meshfieldcomm_cubedom3d, only: meshfieldcommcubedom3d
    use scale_meshfieldcomm_base, only: meshfieldcontainer
    use scale_atm_dyn_dgm_driver_nonhydro3d, only: atmdyndgmdriver_nonhydro3d
 
    use scale_atm_dyn_dgm_nonhydro3d_common, only: &
      auxvar_denshydro_id    
    implicit none
    
    class(atmosvars), intent(inout), target :: this
    class(atmosmesh), intent(in) :: atmos_mesh
    class(atmdyndgmdriver_nonhydro3d), intent(inout) :: dyncore
 
    integer :: iv
 
    integer :: domid
    integer :: ke
    class(meshbase3d), pointer :: mesh3d
    class(localmesh3d), pointer :: lcmesh3d
    class(meshfield3d), pointer :: phyd_ref
    !---------------------------------------
 
    log_newline
    log_info("ATMOSVar_read_restart_file",*) 'Open restart file (ATMOS) '
        
    !- Open restart file
    call this%restart_file%Open()
 
    !- Read restart file
    
    do iv=1, prgvar_num
      call this%restart_file%Read_var( dimtype_xyz, this%PROG_VARS(iv)%varname, &
        this%PROG_VARS(iv)                                                      )
    end do
    do iv=1, auxvar_denshydro_id
      call this%restart_file%Read_var( dimtype_xyz, this%AUX_VARS(iv)%varname, &
        this%AUX_VARS(iv)                                                      )
    end do
    do iv=1, qa
      call this%restart_file%Read_var( dimtype_xyz, this%QTRC_VARS(iv)%varname, &
        this%QTRC_VARS(iv)                                                      )
    end do
 
    !- Close restart file
    log_info("ATMOSVar_read_restart_file",*) 'Close restart file (ATMOS) '
    call this%restart_file%Close()
 
    !-- Diagnostic pressure
    call vars_calc_specific_heat( this )
    
    call dyncore%calc_pressure( this%AUX_VARS(auxvar_pres_id), &
      this%PROGVARS_manager, this%AUXVARS_manager              )
 
    ! Set reference value of hydrostatic pressure
    phyd_ref => this%AUX_VARS(auxvar_preshydro_ref_id)
    mesh3d => phyd_ref%mesh
    do domid=1, mesh3d%LOCAL_MESH_NUM
      lcmesh3d => mesh3d%lcmesh_list(domid)
      !$omp parallel do
      do ke=lcmesh3d%NeS, lcmesh3d%NeE
        phyd_ref%local(domid)%val(:,ke) = 0.0_rp
      end do
    end do
 
    !-- Check read data
    call this%Check( force = .true. )
 
    !-- Calculate diagnostic variables
    call this%Calc_diagnostics()   
 
    !-- Communicate halo data of hydrostatic & diagnostic variables
    call this%AUXVARS_manager%MeshFieldComm_Exchange()
 
    !-- Set horizontal gradient of hydrostatic pressure
    call dyncore%update_phyd_hgrad( this%AUX_VARS(auxvar_preshydro_id), phyd_ref, &
      mesh3d, atmos_mesh%element3D_operation )
 
    return
  end subroutine atmosvar_read_restart_file
 
!OCL SERIAL
  subroutine atmosvar_write_restart_file( this )
    use scale_atm_dyn_dgm_nonhydro3d_common, only: &
      atm_dyn_dgm_nonhydro3d_common_get_varinfo
    
    implicit none
    class(atmosvars), intent(inout) :: this
 
    type(variableinfo) :: prgvar_info(prgvar_num)
    type(variableinfo) :: auxvar_info(auxvar_num)
 
    integer :: iv, rf_vid 
    !---------------------------------------
    
    log_newline
    log_info("ATMOSVar_write_restart_file",*) 'Create restart file (ATMOS) '
 
    !- Check data which will be written to restart file
    call this%Check( force = .true. )
 
    !- Create restart file
    call this%restart_file%Create()
    call prc_mpibarrier()
 
    !- Define variables
 
    call atm_dyn_dgm_nonhydro3d_common_get_varinfo( prgvar_info, auxvar_info )
 
    do iv=1, prgvar_num
      rf_vid = iv
      call this%restart_file%Def_var( this%PROG_VARS(iv),  &
        prgvar_info(iv)%DESC, rf_vid, dimtype_xyz          )
    end do
    do iv=1, auxvar_denshydro_id
      rf_vid = prgvar_num + iv
      call this%restart_file%Def_var( this%AUX_VARS(iv),   &
        auxvar_info(iv)%DESC, rf_vid, dimtype_xyz          )
    end do
    do iv=1, qa
      rf_vid = rf_vid + 1
      call this%restart_file%Def_var( this%QTRC_VARS(iv), &
        tracer_desc(iv), rf_vid, dimtype_xyz              )    
    end do
 
    call this%restart_file%End_def()
 
    !- Write restart file
    do iv=1, prgvar_num
      rf_vid = iv
      call this%restart_file%Write_var(rf_vid, this%PROG_VARS(iv) )
    end do
    do iv=1, auxvar_denshydro_id
      rf_vid = prgvar_num + iv
      call this%restart_file%Write_var(rf_vid, this%AUX_VARS(iv) )
    end do
    do iv=1, qa
      rf_vid = rf_vid + 1
      call this%restart_file%Write_var(rf_vid, this%QTRC_VARS(iv) )
    end do
 
    !- Close restart file
    log_info("ATMOSVar_write_restart_file",*) 'Close restart file (ATMOS) '
    call this%restart_file%Close()
 
    return
  end subroutine atmosvar_write_restart_file
 
!OCL SERIAL
  subroutine atmosvars_check( this, force )
 
    use scale_meshfield_statistics, only: &
      meshfield_statistics_total,         &
      meshfield_statistics_detail
    
    implicit none
    class(atmosvars), intent(inout) :: this
    logical, intent(in), optional :: force
 
    integer :: iv
    integer :: iv_diag
    integer :: n
    logical  :: check
 
    class(meshbase3d), pointer :: mesh3d
    class(localmeshbase), pointer :: lcmesh
    class(localmeshfieldbase), pointer :: lcfield
    type(elementbase), pointer :: elem
    character(len=H_MID) :: varname
 
    type(meshfield3d) :: vel_fields(3)
    type(meshfield3d) :: work
    !--------------------------------------------------------------------------
 
    if ( present(force) ) then
      check = force
    else
      check = this%check_range
    end if
 
    if (check) then
      do iv=1, prgvar_num
        if ( iv == prgvar_therm_id ) cycle
        
        mesh3d => this%PROG_VARS(iv)%mesh
        do n=1, mesh3d%LOCAL_MESH_NUM
          lcmesh => mesh3d%lcmesh_list(n)
          elem => lcmesh%refElem
          call this%PROG_VARS(iv)%GetLocalMeshField(n, lcfield)
          write(varname,'(a,i3.3,a)') this%PROG_VARS(iv)%varname//'(domID=', n, ')' 
          call valcheck( elem%Np, 1, elem%Np, lcmesh%NeA, lcmesh%NeS, lcmesh%NeE, lcfield%val(:,:), &
            progvars_check_min(iv), progvars_check_max(iv), trim(varname), __file__, __line__       )
        end do
      end do
 
      mesh3d => this%PROG_VARS(1)%mesh
      do iv=1, 3
        iv_diag = atmos_diagvars_u_id + iv - 1
        call vel_fields(iv)%Init( atmos_diagvars3d_vinfo(iv_diag)%NAME, "",  mesh3d )
        call atmosvars_calcdiagvar( this, vel_fields(iv)%varname, vel_fields(iv) )
      end do
      call meshfield_statistics_detail( vel_fields(:) )
      do iv=1, 3
        call vel_fields(iv)%Final()
      end do
    end if
 
    if ( present(force) ) then
      check = force
    else
      check = this%check_total
    end if
    if (check) then
      mesh3d => this%PROG_VARS(1)%mesh
      call work%Init("tmp", "", mesh3d)
      call work%Final()
    end if
 
    return
  end subroutine atmosvars_check
 
!OCL SERIAL
  subroutine atmosvars_monitor( this )
    use scale_file_monitor_meshfield, only: &
      file_monitor_meshfield_put
 
    implicit none
    class(atmosvars), intent(inout) :: this
 
    integer :: iv
    class(meshbase3d), pointer :: mesh3d
    type(meshfield3d) :: work
 
    integer :: n
    integer :: ke
    class(localmesh3d), pointer :: lcmesh
    !--------------------------------------------------------------------------
 
    mesh3d => this%PROG_VARS(1)%mesh
    call work%Init("tmp", "", mesh3d)
 
    do iv=1, prgvar_num
      call file_monitor_meshfield_put( this%PROG_VARS(iv)%monitor_id, this%PROG_VARS(iv) )
    end do
  
    do iv=1, qa
      if ( this%QTRC_VARS(iv)%monitor_id > 0 ) then
        do n=1, mesh3d%LOCAL_MESH_NUM
          lcmesh => mesh3d%lcmesh_list(n)
          !$omp parallel do
          do ke=lcmesh%NeS, lcmesh%NeE
            work%local(n)%val(:,ke) = ( this%AUX_VARS(auxvar_denshydro_id)%local(n)%val(:,ke) &
                                      + this%PROG_VARS(prgvar_ddens_id  )%local(n)%val(:,ke) &
                                      ) * this%QTRC_VARS(iv)%local(n)%val(:,ke)
          end do
        end do
        call file_monitor_meshfield_put( this%QTRC_VARS(iv)%monitor_id, work )
      end if
    end do
    if ( dv_monit_id(im_qtot) > 0 ) then
      do n=1, mesh3d%LOCAL_MESH_NUM
        lcmesh => mesh3d%lcmesh_list(n)
        !$omp parallel do
        do ke=lcmesh%NeS, lcmesh%NeE
          work%local(n)%val(:,ke) = ( this%AUX_VARS(auxvar_denshydro_id)%local(n)%val(:,ke) &
                                    + this%PROG_VARS(prgvar_ddens_id   )%local(n)%val(:,ke) &
                                    ) * ( 1.0_rp - this%AUX_VARS(auxvar_qdry_id)%local(n)%val(:,ke) )
        end do
      end do
      call file_monitor_meshfield_put( dv_monit_id(im_qtot), work )
    end if
 
    !##### Energy Budget #####
 
    if ( dv_monit_id(im_engt) > 0 ) then
      call atmosvars_calcdiagvar( this, 'ENGT', work )
      call file_monitor_meshfield_put( dv_monit_id(im_engt), work )
    end if
    if ( dv_monit_id(im_engp) > 0 ) then
      call atmosvars_calcdiagvar( this, 'ENGP', work )
      call file_monitor_meshfield_put( dv_monit_id(im_engp), work )
    end if
    if ( dv_monit_id(im_engk) > 0 ) then
      call atmosvars_calcdiagvar( this, 'ENGK', work )
      call file_monitor_meshfield_put( dv_monit_id(im_engk), work )
    end if
    if ( dv_monit_id(im_engi) > 0 ) then
      call atmosvars_calcdiagvar( this, 'ENGI', work )
      call file_monitor_meshfield_put( dv_monit_id(im_engi), work )
    end if
 
    call work%Final()
 
    return
  end subroutine atmosvars_monitor
 
  !----  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) :: 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)
   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
    class(localmeshbase), pointer :: lcmesh
    !-------------------------------------------------------
 
    call phytends_list%Get(phytend_num1 + qtrcid, field)
    call field%GetLocalMeshField(domid, rhoq_tp)
 
    return
  end subroutine atmosvars_getlocalmeshqtrcphytend  
 
  !-----------------------------------------------------------------------------
!OCL SERIAL  
  subroutine atmosvars_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(atmosvars), intent(inout), target :: this
 
    class(localmesh3d), pointer :: lcmesh3d
    integer :: n
    integer :: varid
    integer :: ke
 
    class(meshfield3d), pointer :: field
    class(elementbase3d), pointer :: elem3d
 
    type(localmeshfieldbaselist) :: qtrc(qa)
    !-------------------------------------------------------
 
    ! Calculate specific heat
    call vars_calc_specific_heat( this )
    
    ! Calculate diagnostic variables
    do varid=auxvar_denshydro_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
 
    return
  end subroutine atmosvars_calculatediagnostics
 
!OCL SERIAL
  subroutine atmosvars_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(atmosvars), intent(inout) :: this
    character(*), intent(in) :: field_name
    type(meshfield3d), intent(inout) :: field_work
 
    class(localmesh3d), pointer :: lcmesh3d
    class(elementbase3d), pointer :: elem3d
    integer :: n
    integer :: ke
    integer :: iq
 
    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)
        !$omp parallel do
        do ke=lcmesh3d%NeS, lcmesh3d%NeE
          field_work%local(n)%val(:,ke) = field_work_uvmet(uvmet_i)%local(n)%val(:,ke) &
            / ( this%AUX_VARS (auxvar_denshydro_id)%local(n)%val(:,ke)                 &
              + this%PROG_VARS(prgvar_ddens_id   )%local(n)%val(:,ke)                  )
        end do
        call field_work_uvmet(1)%Final()
        call field_work_uvmet(2)%Final()
      end if
    end do
 
    return
  end subroutine atmosvars_calcdiagvar
 
!OCL SERIAL
  subroutine atmosvars_calcdiagvar2d( this, field_name, field_work ) 
    implicit none
    class(atmosvars), intent(inout) :: this
    character(*), intent(in) :: field_name
    type(meshfield2d), intent(inout), target :: field_work
 
    class(localmesh2d), pointer :: lcmesh2d
    class(elementbase2d), pointer :: elem2d
 
    integer :: n
    !--------------------------------------------------
 
    field_work%varname = field_name
 
    do n=1, field_work%mesh%LOCAL_MESH_NUM
      lcmesh2d => field_work%mesh%lcmesh_list(n)
      call vars_calc_diagnosevar2d_lc( field_name, field_work%local(n)%val,  &
        this%ptr_MP_AUXVARS2D_manager,                                       &
        field_work%mesh, lcmesh2d, lcmesh2d%refElem2D                        )
    end do
 
    return
  end subroutine atmosvars_calcdiagvar2d
 
!-- private -----------------------------------------------------------------------
    
!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_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 :: iq
    real(rp) :: dens(elem%np), temp(elem%np)
    real(rp) :: mom_u1(elem%np), mom_u2(elem%np), g_11(elem%np), g_12(elem%np), g_22(elem%np)
    real(rp) :: psat(elem%np)
 
    integer :: iq_qv
    !-------------------------------------------------------------------------
 
    select case(trim(field_name))
    case('DENS')
      !$omp parallel do
      do ke=1, lcmesh%Ne
        var_out(:,ke) = ddens_(:,ke) + dens_hyd(:,ke)
      end do
    
    case('U')
      !$omp parallel do private (DENS)
      do ke=1, lcmesh%Ne
        dens(:) = ddens_(:,ke) + dens_hyd(:,ke)
        var_out(:,ke) = momx_(:,ke) / dens(:)
      end do
    
    case('V')
      !$omp parallel do private (DENS)
      do ke=1, lcmesh%Ne
        dens(:) = ddens_(:,ke) + dens_hyd(:,ke)
        var_out(:,ke) = momy_(:,ke) / dens(:)
      end do  
          
    case('W')
      !$omp parallel do private (DENS)
      do ke=1, lcmesh%Ne
        dens(:) = ddens_(:,ke) + dens_hyd(:,ke)
        var_out(:,ke) = momz_(:,ke) / dens(:)
      end do
    
    case ( 'PRES' )
    case('PRES_diff')  
      !$omp parallel do
      do ke=1, lcmesh%Ne
        var_out(:,ke) = pres_(:,ke) - pres_hyd(:,ke)
      end do
    
    case('T')
      !$omp parallel do
      do ke=1, lcmesh%Ne
        var_out(:,ke) = pres_(:,ke) / (rtot(:,ke) * (ddens_(:,ke) + dens_hyd(:,ke)) )
      end do
    
    case('T_diff')
      !$omp parallel do
      do ke=1, lcmesh%Ne
        var_out(:,ke) = pres_(:,ke) / ( rtot(:,ke) * (ddens_(:,ke) + dens_hyd(:,ke)) ) &
                      - pres_hyd(:,ke) / ( rdry * dens_hyd(:,ke) )
      end do
    
    case('PT')
      !$omp parallel do private( DENS )
      do ke=1, lcmesh%Ne
        dens(:) = ddens_(:,ke) + dens_hyd(:,ke)
        var_out(:,ke) = pres_(:,ke) / (rtot(:,ke) * dens(:) ) * ( pres00 / pres_(:,ke) )**( rtot(:,ke) / cptot(:,ke) )
      end do 
    
    case('PT_diff')
      !$omp parallel do private( DENS )
      do ke=1, lcmesh%Ne
        dens(:) = ddens_(:,ke) + dens_hyd(:,ke)
        var_out(:,ke) = pres_(:,ke) / (rtot(:,ke) * dens(:) ) * ( pres00 / pres_(:,ke) )**( rtot(:,ke) / cptot(:,ke) ) &
                      - pres00/rdry * (pres_hyd(:,ke)/pres00)**(cvdry/cpdry) / dens_hyd(:,ke)
      end do 
    
    case( 'RH', 'RHL' )
      if ( atmos_hydrometeor_dry ) then
        var_out(:,ke) = 0.0_rp
      else
        call tracer_inq_id( "QV", iq_qv )
 
        !$omp parallel do private (TEMP, PSAT)
        do ke=1, lcmesh%Ne
          temp(:) = pres_(:,ke) / (rtot(:,ke) * (ddens_(:,ke) + dens_hyd(:,ke)) )
 
          call atmos_saturation_psat_liq( &
            elem%Np, 1, elem%Np, temp(:),     & ! (in)
            psat(:)                           ) ! (out)
 
          var_out(:,ke) = ( ddens_(:,ke) + dens_hyd(:,ke) ) * qtrc(iq_qv)%ptr%val(:,ke) &
                        / psat(:) * rvap * temp(:) * 100.0_rp
        end do 
      end if
    
    case('ENGK')
      !$omp parallel do private (ke2D, DENS, mom_u1, mom_u2, G_11, G_12, G_22)
      do ke=1, lcmesh%Ne
        ke2d = lcmesh%EMap3Dto2D(ke)
 
        dens(:) = ddens_(:,ke) + dens_hyd(:,ke)
        g_11(:) = lcmesh%G_ij(elem%IndexH2Dto3D,ke2d,1,1)
        g_12(:) = lcmesh%G_ij(elem%IndexH2Dto3D,ke2d,1,2)
        g_22(:) = lcmesh%G_ij(elem%IndexH2Dto3D,ke2d,2,2)
 
        mom_u1(:) = g_11(:) * momx_(:,ke) + g_12(:) * momy_(:,ke)
        mom_u2(:) = g_12(:) * momx_(:,ke) + g_22(:) * momy_(:,ke)
 
        var_out(:,ke) = 0.5_rp * ( momx_(:,ke) * mom_u1(:) + momy_(:,ke) * mom_u2(:) + momz_(:,ke)**2 ) / dens(:)
      end do
    
    case('ENGP')
      !$omp parallel do private (DENS)
      do ke=1, lcmesh%Ne
        dens(:) = ddens_(:,ke) + dens_hyd(:,ke)
        var_out(:,ke) = dens(:) * grav * lcmesh%zlev(:,ke)
      end do
    
    case('ENGI')
      !$omp parallel do private (DENS, iq)
      do ke=1, lcmesh%Ne
        dens(:) = ddens_(:,ke) + dens_hyd(:,ke)
        var_out(:,ke) = qdry_(:,ke) * pres_(:,ke) / rtot(:,ke) * cvdry
        do iq = 1, qa
          var_out(:,ke) = var_out(:,ke) &
            + qtrc(iq)%ptr%val(:,ke) * ( pres_(:,ke) / rtot(:,ke) * tracer_cv(iq) + dens(:) * tracer_engi0(iq) )
        end do
      end do
    
    case('ENGT')
      !$omp parallel do private (ke2D, DENS, mom_u1, mom_u2, iq, G_11, G_12, G_22)
      do ke=1, lcmesh%Ne
        ke2d = lcmesh%EMap3Dto2D(ke)
 
        dens(:) = ddens_(:,ke) + dens_hyd(:,ke)
 
        g_11(:) = lcmesh%G_ij(elem%IndexH2Dto3D,ke2d,1,1)
        g_12(:) = lcmesh%G_ij(elem%IndexH2Dto3D,ke2d,1,2)
        g_22(:) = lcmesh%G_ij(elem%IndexH2Dto3D,ke2d,2,2)
        mom_u1(:) = g_11(:) * momx_(:,ke) + g_12(:) * momy_(:,ke)
        mom_u2(:) = g_12(:) * momx_(:,ke) + g_22(:) * momy_(:,ke)
 
        ! ENGI
        var_out(:,ke) = qdry_(:,ke) * pres_(:,ke) / rtot(:,ke) * cvdry
        do iq = 1, qa
          var_out(:,ke) = var_out(:,ke) &
            + qtrc(iq)%ptr%val(:,ke) * ( pres_(:,ke) / rtot(:,ke) * tracer_cv(iq) + dens(:) * tracer_engi0(iq) )
        end do
        ! ENGT
        var_out(:,ke) = &
            0.5_rp * ( momx_(:,ke) * mom_u1(:) + momy_(:,ke) * mom_u2(:) + momz_(:,ke)**2 ) / dens(:) & ! ENGK       
          + var_out(:,ke)                                                                             & ! ENGI
          + dens(:) * grav * lcmesh%pos_en(:,ke,3)                                                      ! ENGP
      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
 
!OCL SERIAL
  subroutine vars_calc_diagnosevar2d_lc( field_name, & ! (in)
    var_out,                                         & ! (out)
    mp_auxvars2d, mesh2d, lcmesh, elem )               ! (in)
 
    use mod_atmos_phy_mp_vars, only: &
      atmosphympvars_getlocalmeshfields_sfcflx
 
    implicit none
    class(localmesh2d), intent(in) :: lcmesh
    class(elementbase2d), intent(in) :: elem
    character(*), intent(in) :: field_name
    real(rp), intent(out) :: var_out(elem%np,lcmesh%nea)
    class(modelvarmanager), intent(inout) :: mp_auxvars2d
    class(meshbase2d), intent(in) :: mesh2d
 
    integer :: ke
 
    class(localmeshfieldbase), pointer :: sflx_rain_mp, sflx_snow_mp, sflx_engi_mp
    !-------------------------------------------------------------------------
 
    select case(trim(field_name))
    case('RAIN', 'SNOW')
      call atmosphympvars_getlocalmeshfields_sfcflx( &
        lcmesh%lcdomID, mesh2d, mp_auxvars2d,        &
        sflx_rain_mp, sflx_snow_mp, sflx_engi_mp     )
    end select
 
    select case(trim(field_name))
    case('RAIN')
      !$omp parallel do
      do ke=lcmesh%NeS, lcmesh%NeE
        var_out(:,ke) = sflx_rain_mp%val(:,ke)
      end do
    case('SNOW')
      !$omp parallel do
      do ke=lcmesh%NeS, lcmesh%NeE
        var_out(:,ke) = sflx_snow_mp%val(:,ke)
      end do
    case default
      log_error("AtmosVars_calc_diagnoseVar2D_lc",*) 'The name of diagnostic variable is not suported. Check!', field_name
      call prc_abort
    end select
 
    return
  end subroutine vars_calc_diagnosevar2d_lc
 
!OCL SERIAL  
  subroutine vars_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(atmosvars), intent(inout), target :: this
 
    class(localmesh3d), pointer :: lcmesh3d
    integer :: n
    integer :: varid
    integer :: ke
    integer :: iq
 
    class(elementbase3d), pointer :: elem3d
 
    real(rp), allocatable :: q_tmp(:,:)
    !-------------------------------------------------------
 
    ! Calculate specific heat
    do n=1, this%AUX_VARS(1)%mesh%LOCAL_MESH_NUM
      lcmesh3d => this%AUX_VARS(1)%mesh%lcmesh_list(n)
      elem3d => lcmesh3d%refElem3D
      allocate( q_tmp(elem3d%Np,qa) )
 
      !$omp parallel do private(ke, iq, q_tmp)
      do ke = lcmesh3d%NeS, lcmesh3d%NeE
        do iq = 1, qa
          q_tmp(:,iq) = this%QTRC_VARS(iq)%local(n)%val(:,ke)
        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)
          this%AUX_VARS(auxvar_qdry_id )%local(n)%val(:,ke),            & ! (out)
          this%AUX_VARS(auxvar_rtot_id )%local(n)%val(:,ke),            & ! (out)
          this%AUX_VARS(auxvar_cvtot_id)%local(n)%val(:,ke),            & ! (out)
          this%AUX_VARS(auxvar_cptot_id)%local(n)%val(:,ke)             ) ! (out)
      end do
      deallocate(q_tmp)
    end do
  end subroutine vars_calc_specific_heat
 
end module mod_atmos_vars