mod_atmos_vars.F90

Go to the documentation of this file.

!-------------------------------------------------------------------------------
!> module Atmosphere / Variables
!!
!! @par Description
!!          Module to manage variables with atmospheric component
!!
!! @author Yuta Kawai, Team SCALE
!!
!<
!-------------------------------------------------------------------------------
#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
 
  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_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_vars_container, only: &
    atmosvarscontainer, &
    atmosvars_getlocalmeshprgvar, atmosvars_getlocalmeshprgvars,       &
    atmosvars_getlocalmeshsfcvar,                                      &
    atmosvars_getlocalmeshqtrcvar, atmosvars_getlocalmeshqtrcvarlist,  &
    atmosvars_getlocalmeshqtrc_qv,                                     &
    atmosvars_getlocalmeshphyauxvars,                                  &
    atmosvars_getlocalmeshphytends, atmosvars_getlocalmeshqtrcphytend, &
    atmos_auxvars2d_num,                                               &
    atm_vars_container_primary_id
    
  !-----------------------------------------------------------------------------
  implicit none
  private
 
  !-----------------------------------------------------------------------------
  !
  !++ Public type & procedures
  !
  integer, parameter :: ATM_VARS_CONTAINER_LIST_MAX = 16
 
  !> Derived type to manage variables with atmospheric component
  type, public :: atmosvars
    class(atmosmesh), pointer :: mesh                !< Pointer to an object to manage a primary mesh for atmospheric component
 
    type(atmosvarscontainer) :: container_list(atm_vars_container_list_max) !< Array of containers to manage variables. The first one is the primary container.
    logical :: init_containers_item_flag(atm_vars_container_list_max)       !< Flag whether the item of container_list is initialized or not. 
    
    type(atmosvarscontainer), pointer :: container              !< Pointer to a primary container to manage variables
    type(modelvarmanager), pointer :: progvars_manager          !< Pointer to manage prognostic variables with primary container
    type(modelvarmanager), pointer :: qtrcvars_manager          !< Pointer to manage tracer variables with primary container
    type(modelvarmanager), pointer :: auxvars_manager           !< Pointer to manage auxiliary variables with primary container
    type(modelvarmanager), pointer :: auxvars2d_manager         !< Pointer to manage 2D auxiliary variables with primary container
    type(modelvarmanager), pointer :: phytends_manager          !< Pointer to manage tendency variables with physics
   
    !-
    type(modelvarmanager), pointer :: ptr_mp_auxvars2d_manager
    type(modelvarmanager), pointer :: ptr_cp_auxvars2d_manager
 
    logical :: moist
    type(meshfield3d), pointer :: qv
    type(meshfield3d) :: zero
 
    !--
    integer, allocatable :: diagvars2d_histid(:)
    integer, allocatable :: diagvars3d_histid(:)
 
    type(file_restart_meshfield_component) :: restart_file !< Object to manage restart file for atmospheric component
    
    character(len=H_MID) :: phy_preproc_file_basename !< Basename of configuration file for preprocesses before physics
 
    logical :: check_range
    logical :: check_total
 
  contains
    procedure :: init => atmosvars_init
    procedure :: final => atmosvars_final
    procedure :: setup_container => atmosvars_setup_container   
    procedure :: get_container => atmosvars_get_container 
    procedure :: calc_diagnostics => atmosvars_calculatediagnostics
    procedure :: calc_diagvar => atmosvars_calcdiagvar
    procedure :: calc_diagvar2d => atmosvars_calcdiagvar2d
    procedure :: preprocoperationforphys => atmosvars_preprocoperationforphys
    procedure :: history => atmosvars_history
    procedure :: check   => atmosvars_check
    procedure :: monitor => atmosvars_monitor
    procedure :: read_restart_file => atmosvars_read_restart_file
    procedure :: write_restart_file_prep => atmosvars_write_restart_file_prep
    procedure :: write_restart_file => atmosvars_write_restart_file
    procedure :: write_restart_file_post => atmosvars_write_restart_file_post
    procedure :: regist_physvar_manager => atmosvars_regist_physvar_manager
  end type atmosvars
 
  ! Cascade
  public :: atmosvarscontainer
  public :: atm_vars_container_primary_id
  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 /)
                      
 
  ! 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_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
 
!> Setup an object to manage variables with atmospheric component
!!
!! @param model_mesh Object to manage computational mesh of atmospheric model 
!!
!OCL SERIAL
  subroutine atmosvars_init( this, atm_mesh )
    use scale_atm_dyn_dgm_nonhydro3d_common, only: &
      prgvar_scalar_num, prgvar_hvec_num    
    use scale_file_monitor_meshfield, only:    &
      monitor_reg => file_monitor_meshfield_reg
    implicit none
 
    class(atmosvars), target, intent(inout) :: this
    class(atmosmesh), target, intent(inout) :: atm_mesh
 
    integer :: iv
    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.  !< Flag whether the range of values is checked
    logical :: check_total    = .false.
 
    character(len=H_MID) :: phy_preproc_file_basename = 'phy_preoperation' !< Basename of configuration file for preprocesses before physics
 
    namelist / param_atmos_vars / &
      phy_preproc_file_basename, &
      check_range, &
      check_total
 
    character(len=H_LONG) :: in_basename           = ''        !< Basename of the input  file
    logical :: in_postfix_timelabel                = .false.   !< Add timelabel to the basename of input  file?
    character(len=H_LONG) :: out_basename          = ''        !< Basename of the output file
    logical :: out_postfix_timelabel               = .true.    !< Add timelabel to the basename of output file?
    character(len=H_MID) :: out_title              = ''        !< Title    of the output file
    character(len=H_SHORT) :: out_dtype            = 'DEFAULT' !< REAL4 or REAL8  
 
    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
    !--------------------------------------------------
 
    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)
 
    this%phy_preproc_file_basename = phy_preproc_file_basename
 
    !- 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)
 
    this%init_containers_item_flag(:) = .false.
    call this%Setup_container( 1, atm_mesh )  
 
    !- 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
 
    !-- Set the pointer of 2D auxiliary variable manager with MP and CP components to output precipitation fluxes
 
    nullify( this%ptr_MP_AUXVARS2D_manager )
    nullify( this%ptr_CP_AUXVARS2D_manager )
 
    return
  end subroutine atmosvars_init
 
  !> Setup a container to manage variables with atmospheric component
  !!
!OCL SERIAL
  subroutine atmosvars_setup_container( this, container_type, atm_mesh )
    implicit none
    class(atmosvars), intent(inout), target :: this
    integer, intent(in) :: container_type                   !< Type of container to setup. 1: primary container, >=2: other containers
    class(atmosmesh), target, intent(inout) :: atm_mesh
    !---------------------------------------------
    
    if ( container_type < atm_vars_container_primary_id .or. container_type > atm_vars_container_list_max ) then
      log_error("ATMOS_vars_setup_container",*) 'container_type is out of range. Check!', container_type
      call prc_abort
    end if
 
    if ( .not. this%init_containers_item_flag(container_type) ) then
      log_info("ATMOS_vars_setup_container",*) 'container_type: ', container_type
 
      this%init_containers_item_flag(container_type) = .true.
      call this%container_list(container_type)%Init( container_type, this%phy_preproc_file_basename, atm_mesh )
 
      if ( container_type == atm_vars_container_primary_id ) then
        this%container => this%container_list(container_type)
        this%PROGVARS_manager => this%container%PROGVARS_manager
        this%QTRCVARS_manager => this%container%QTRCVARS_manager
        this%AUXVARS_manager  => this%container%AUXVARS_manager
        this%AUXVARS2D_manager  => this%container%AUXVARS2D_manager
        this%PHYTENDS_manager  => this%container%PHYTENDS_manager
      end if
    end if
    return    
  end subroutine atmosvars_setup_container
 
  !> Get a container to manage variables with atmospheric component
  !!
!OCL SERIAL
  subroutine atmosvars_get_container( this, container_type, &
    container )
    implicit none
    class(atmosvars), intent(inout), target :: this
    integer, intent(in) :: container_type            !< Type of container to get. 1: primary container, >=2: other containers
    class(atmosvarscontainer), pointer :: container  !< Pointer to a container to manage variables with atmospheric component
    !---------------------------------------------
 
    if ( container_type < atm_vars_container_primary_id .or. container_type > atm_vars_container_list_max ) then
      log_error("ATMOS_vars_get_container",*) 'container_type is out of range. Check!', container_type
      call prc_abort
    end if
 
    if ( this%init_containers_item_flag(container_type) ) then
      container => this%container_list(container_type)
    else
      log_error("ATMOS_vars_get_container",*) 'Not initialized container. Check!'
      call prc_abort
    end if
    return
  end subroutine atmosvars_get_container
 
  !> Finalize an object to manage variables with atmospheric component
  !!
!OCL SERIAL
  subroutine atmosvars_final( this )
    implicit none
    class(atmosvars), intent(inout) :: this
 
    integer :: i    
    !--------------------------------------------------
 
    log_info('AtmosVars_Final',*)
 
    call this%restart_file%Final()
 
    do i=1, atm_vars_container_list_max
      if ( this%init_containers_item_flag(i) ) call this%container_list(i)%Final()
    end do
 
    deallocate( this%DIAGVARS3D_HISTID )
 
    return
  end subroutine atmosvars_final
 
  !> Set the pointer of 2D ModelVarManager with cloud microphysics or cumlus parameterization
  !! to output surface variables with precipitation fluxes
!OCL SERIAL
  subroutine atmosvars_regist_physvar_manager( this, &
    mp_AUXVARS2D_manager, cp_AUXVARS2D_manager )
    implicit none
 
    class(atmosvars), target, intent(inout) :: this
    type(modelvarmanager), intent(in), target, optional:: mp_auxvars2d_manager
    type(modelvarmanager), intent(in), target, optional:: cp_auxvars2d_manager
    !----------------------------------------------
 
    if ( present(mp_auxvars2d_manager) ) then
      this%ptr_MP_AUXVARS2D_manager => mp_auxvars2d_manager
    end if
    if ( present(cp_auxvars2d_manager) ) then
      this%ptr_CP_AUXVARS2D_manager => cp_auxvars2d_manager
    end if
    return
  end subroutine atmosvars_regist_physvar_manager
 
  !> Put data with atmospheric variables to history file
  !!
!OCL SERIAL
  subroutine atmosvars_history( this )
    use scale_file_history_meshfield, only: file_history_meshfield_put
    implicit none
    class(atmosvars), intent(inout), target :: this
  
    integer :: v
    integer :: hst_id
 
    type(meshfield3d) :: tmp_field3d
    class(meshbase3d), pointer :: mesh3d
 
    type(meshfield2d) :: tmp_field2d
    class(meshbase2d), pointer :: mesh2d
    !-------------------------------------------------------------------------
 
    mesh3d => this%container%PROG_VARS(1)%mesh
    call mesh3d%GetMesh2D(mesh2d)
 
    !-
    call this%container%Calc_diagnostics()
 
    do v = 1, prgvar_num
      hst_id = this%container%PROG_VARS(v)%hist_id
      if ( hst_id > 0 ) call file_history_meshfield_put( hst_id, this%container%PROG_VARS(v) )
    end do
 
    do v = 1, qa
      hst_id = this%container%QTRC_VARS(v)%hist_id
      if ( hst_id > 0 ) call file_history_meshfield_put( hst_id, this%container%QTRC_VARS(v) )
    end do
 
    do v = 1, auxvar_num
      hst_id = this%container%AUX_VARS(v)%hist_id
      if ( hst_id > 0 ) call file_history_meshfield_put( hst_id, this%container%AUX_VARS(v) )
    end do
    do v = 1, atmos_auxvars2d_num
      hst_id = this%container%AUX_VARS2D(v)%hist_id
      if ( hst_id > 0 ) call file_history_meshfield_put( hst_id, this%container%AUX_VARS2D(v) )
    end do
 
    do v = 1, this%container%PHYTEND_NUM_TOT
      hst_id = this%container%PHY_TEND(v)%hist_id
      if ( hst_id > 0 ) call file_history_meshfield_put( hst_id, this%container%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
 
!> Read data with atmospheric variables from restart file
!!
!OCL SERIAL
  subroutine atmosvars_read_restart_file( this, atmos_mesh, dyncore  )
    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("AtmosVars_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%container%PROG_VARS(iv)%varname, &
        this%container%PROG_VARS(iv)                                                      )
    end do
    do iv=1, auxvar_denshydro_id
      call this%restart_file%Read_var( dimtype_xyz, this%container%AUX_VARS(iv)%varname, &
        this%container%AUX_VARS(iv)                                                      )
    end do
    do iv=1, qa
      call this%restart_file%Read_var( dimtype_xyz, this%container%QTRC_VARS(iv)%varname, &
        this%container%QTRC_VARS(iv)                                                      )
    end do
 
    !- Close restart file
    log_info("AtmosVars_Read_restart_file",*) 'Close restart file (ATMOS) '
    call this%restart_file%Close()
 
    !-- Prepare diagnostic variables
 
    ! Calculate specific heat
    call this%container%Calc_SpecificHeat()
 
    ! Set a basic state of thermodynamics variable
    call dyncore%update_therm_hyd( this%container%AUXVARS_manager )
 
    ! Calculate pressure
    call dyncore%calc_pressure( this%container%AUX_VARS(auxvar_pres_id), &
      this%container%PROGVARS_manager, this%container%AUXVARS_manager    )
 
    ! Set reference value of hydrostatic pressure
    phyd_ref => this%container%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
      !$acc update device(Phyd_ref%local(domid)%val)
    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%container%AUXVARS_manager%MeshFieldComm_Exchange()
 
    !-- Set horizontal gradient of hydrostatic pressure
    call dyncore%update_phyd_hgrad( this%container%AUX_VARS(auxvar_preshydro_id), phyd_ref, &
      mesh3d, atmos_mesh%element3D_operation )
 
    return
  end subroutine atmosvars_read_restart_file
 
!> Write data with atmospheric variables to restart file
!!
!OCL SERIAL
  subroutine atmosvars_write_restart_file_prep( this )
    use scale_tracer, only: &
      tracer_desc    
    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("AtmosVars_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%container%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%container%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%container%QTRC_VARS(iv), &
        tracer_desc(iv), rf_vid, dimtype_xyz                        )    
    end do
    return
  end subroutine atmosvars_write_restart_file_prep
 
!> Write data with atmospheric variables to restart file
!!
!OCL SERIAL
  subroutine atmosvars_write_restart_file( this )
    implicit none
    class(atmosvars), intent(inout) :: this
 
    integer :: iv, rf_vid 
    !---------------------------------------
 
    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%container%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%container%AUX_VARS(iv) )
    end do
    do iv=1, qa
      rf_vid = rf_vid + 1
      call this%restart_file%Write_var(rf_vid, this%container%QTRC_VARS(iv) )
    end do
    return
  end subroutine atmosvars_write_restart_file
 
  !> Close restart file with atmospheric variables
  !!
!OCL SERIAL
  subroutine atmosvars_write_restart_file_post( this )
    implicit none
    class(atmosvars), intent(inout) :: this
    !---------------------------------------
    log_info("AtmosVars_Write_restart_file_post",*) 'Close restart file (ATMOS) '
 
    call this%restart_file%Close()
    return
  end subroutine atmosvars_write_restart_file_post
 
!> Check the range of values with atmospheric variables
!!
!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), pointer :: var    
    type(meshfield3d) :: vel_fields(3)
    type(meshfield3d) :: work
 
    integer :: ldomid
    !--------------------------------------------------------------------------
 
    if ( present(force) ) then
      check = force
    else
      check = this%check_range
    end if
 
    if (check) then
      mesh3d => this%container%PROG_VARS(1)%mesh
 
      do iv=1, prgvar_num
        if ( iv == prgvar_therm_id ) cycle
        
        var => this%container%PROG_VARS(iv)        
        do n=1, mesh3d%LOCAL_MESH_NUM
          lcmesh => mesh3d%lcmesh_list(n)
          elem => lcmesh%refElem
 
          call var%GetLocalMeshField(n, lcfield)
 
          write(varname,'(a,i3.3,a)') var%varname//'(domID=', n, ')' 
 
          ! Note : *acc update host* is called in VALCHECK of SCALE library. 
          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
 
      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) )
#ifdef _OPENACC        
        do ldomid=1, mesh3d%LOCAL_MESH_NUM
          !$acc update host(vel_fields(iv)%local(ldomID)%val)
        end do
#endif        
      end do
      call meshfield_statistics_detail( vel_fields )
      do iv=1, 3
        call vel_fields(iv)%Final()
      end do
 
      ! Check total
      call work%Init("tmp", "", mesh3d)
      call work%Final()
    end if
 
    return
  end subroutine atmosvars_check
 
!> Put the stastics with atmospheric variables
!!
!OCL SERIAL
  subroutine atmosvars_monitor( this )
    implicit none
    class(atmosvars), intent(inout) :: this
    !--------------------------------------------------------------------------
    call atmosvars_monitor_core( this%container )
    return
  end subroutine atmosvars_monitor
!OCL SERIAL
  subroutine atmosvars_monitor_core( this )  
    use scale_file_monitor_meshfield, only: &
      file_monitor_meshfield_put
    implicit none
    class(atmosvarscontainer), 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 this%Calc_diagVar( 'ENGT', work )
      call file_monitor_meshfield_put( dv_monit_id(im_engt), work )
    end if
    if ( dv_monit_id(im_engp) > 0 ) then
      call this%Calc_diagVar( 'ENGP', work )
      call file_monitor_meshfield_put( dv_monit_id(im_engp), work )
    end if
    if ( dv_monit_id(im_engk) > 0 ) then
      call this%Calc_diagVar( 'ENGK', work )
      call file_monitor_meshfield_put( dv_monit_id(im_engk), work )
    end if
    if ( dv_monit_id(im_engi) > 0 ) then
      call this%Calc_diagVar( 'ENGI', work )
      call file_monitor_meshfield_put( dv_monit_id(im_engi), work )
    end if
 
    call work%Final()
 
    return
  end subroutine atmosvars_monitor_core
 
  !> Preprocess operation for physical processes
  !!
!OCL SERIAL
  subroutine atmosvars_preprocoperationforphys( this, dyncore )
    use scale_atm_dyn_dgm_driver_nonhydro3d, only: atmdyndgmdriver_nonhydro3d
    implicit none
    class(atmosvars), intent(inout) :: this
    class(atmdyndgmdriver_nonhydro3d), intent(in) :: dyncore
 
    integer :: ic
    !--------------------------------------------------------------------------
 
    do ic=atm_vars_container_primary_id+1, atm_vars_container_list_max
      if ( this%init_containers_item_flag(ic) ) then
        call this%container_list(ic)%Preproc_operation_for_phys( this%container, dyncore )
      end if
    end do
    return
  end subroutine atmosvars_preprocoperationforphys
 
  !-----------------------------------------------------------------------------
  !> Calculate diagnostic variables
!OCL SERIAL  
  subroutine atmosvars_calculatediagnostics( this )
    implicit none
    class(atmosvars), intent(inout), target :: this
    !-------------------------------------------------------
    call this%container%Calc_diagnostics()
    return
  end subroutine atmosvars_calculatediagnostics
 
!OCL SERIAL
  subroutine atmosvars_calcdiagvar( this, field_name, field_work ) 
    implicit none
    class(atmosvars), intent(inout) :: this
    character(*), intent(in) :: field_name
    type(meshfield3d), intent(inout) :: field_work
    !--------------------------------------------------
    call this%container%Calc_diagVar( field_name, field_work )
    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
 
    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, this%ptr_CP_AUXVARS2D_manager,        &
        field_work%mesh, lcmesh2d, lcmesh2d%refElem2D                        )
    end do
    !$acc wait(1)
    return
  end subroutine atmosvars_calcdiagvar2d
 
!--- private -----
 
!OCL SERIAL
  subroutine vars_calc_diagnosevar2d_lc( field_name, & ! (in)
    var_out,                                         & ! (out)
    mp_auxvars2d, cp_auxvars2d, mesh2d, lcmesh, elem )  ! (in)
 
    use mod_atmos_phy_mp_vars, only: &
      atmosphympvars_getlocalmeshfields_sfcflx
    use mod_atmos_phy_cp_vars, only: &
      atmosphycpvars_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)
    type(modelvarmanager), intent(inout), pointer :: mp_auxvars2d
    type(modelvarmanager), intent(inout), pointer :: cp_auxvars2d
    class(meshbase2d), intent(in) :: mesh2d
 
    integer :: ke, p
 
    logical :: sw_mp, sw_cp
    class(localmeshfieldbase), pointer :: sflx_rain_mp, sflx_snow_mp, sflx_engi_mp
    class(localmeshfieldbase), pointer :: sflx_rain_cp, sflx_snow_cp, sflx_engi_cp
    !-------------------------------------------------------------------------
 
    sw_mp = associated(mp_auxvars2d)
    sw_cp = associated(cp_auxvars2d)
 
    select case(trim(field_name))
    case('RAIN', 'SNOW')
      if ( sw_mp ) then
        call atmosphympvars_getlocalmeshfields_sfcflx( &
          lcmesh%lcdomID, mesh2d, mp_auxvars2d,        &
          sflx_rain_mp, sflx_snow_mp, sflx_engi_mp     )
      end if
      if ( sw_cp ) then
        call atmosphycpvars_getlocalmeshfields_sfcflx( &
          lcmesh%lcdomID, mesh2d, cp_auxvars2d,        &
          sflx_rain_cp, sflx_snow_cp, sflx_engi_cp     )
      end if
    end select
 
    select case(trim(field_name))
    case('RAIN')
      !$omp parallel 
      !$acc parallel present(var_out) async(1)
      !$omp do
      !$acc loop collapse(2) 
      do ke=lcmesh%NeS, lcmesh%NeE
      do p=1, elem%Np
        var_out(p,ke) = 0.0_rp
      end do
      end do
      if ( sw_mp ) then
        !$omp do
        !$acc loop collapse(2)
        do ke=lcmesh%NeS, lcmesh%NeE
        do p=1, elem%Np
          var_out(p,ke) = var_out(p,ke) + sflx_rain_mp%val(p,ke)
        end do
        end do
      end if
      if ( sw_cp ) then
        !$omp do
        !$acc loop collapse(2)
        do ke=lcmesh%NeS, lcmesh%NeE
        do p=1, elem%Np
          var_out(p,ke) = var_out(p,ke) + sflx_rain_cp%val(p,ke)
        end do
        end do
      end if
      !$omp end parallel
      !$acc end parallel
    case('SNOW')
      !$omp parallel 
      !$acc parallel present(var_out) async(1)
      !$omp do
      !$acc loop collapse(2) 
      do ke=lcmesh%NeS, lcmesh%NeE
      do p=1, elem%Np
        var_out(p,ke) = 0.0_rp
      end do
      end do
      if ( sw_mp ) then
        !$omp do
        !$acc loop collapse(2)
        do ke=lcmesh%NeS, lcmesh%NeE
        do p=1, elem%Np
          var_out(p,ke) = var_out(p,ke) + sflx_snow_mp%val(p,ke)
        end do
        end do
      end if
      if ( sw_cp ) then
        !$omp do
        !$acc loop collapse(2)
        do ke=lcmesh%NeS, lcmesh%NeE
        do p=1, elem%Np
          var_out(p,ke) = var_out(p,ke) + sflx_snow_cp%val(p,ke)
        end do
        end do
      end if
      !$omp end parallel
      !$acc end parallel
    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  
end module mod_atmos_vars