mod_atmos_phy_mp_vars.F90

Go to the documentation of this file.

!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------
#include "scalelib.h"
module mod_atmos_phy_mp_vars
  !-----------------------------------------------------------------------------
  !
  !++ used modules
  !
  use scale_precision
  use scale_io
  use scale_prof
  use scale_debug
  use scale_atmos_grid_cartesc_index
  use scale_tracer
  !-----------------------------------------------------------------------------
  implicit none
  private
  !-----------------------------------------------------------------------------
  !
  !++ Public procedure
  !
  public :: atmos_phy_mp_vars_setup
  public :: atmos_phy_mp_vars_finalize
  public :: atmos_phy_mp_vars_fillhalo
  public :: atmos_phy_mp_vars_restart_read
  public :: atmos_phy_mp_vars_restart_write
 
  public :: atmos_phy_mp_vars_restart_create
  public :: atmos_phy_mp_vars_restart_open
  public :: atmos_phy_mp_vars_restart_def_var
  public :: atmos_phy_mp_vars_restart_enddef
  public :: atmos_phy_mp_vars_restart_close
 
  public :: atmos_phy_mp_vars_history
 
  public :: atmos_phy_mp_vars_get_diagnostic
  public :: atmos_phy_mp_vars_reset_diagnostics
 
  !-----------------------------------------------------------------------------
  !
  !++ Public parameters & variables
  !
  logical,                public :: atmos_phy_mp_restart_output                = .false.                
 
  character(len=H_LONG),  public :: atmos_phy_mp_restart_in_basename           = ''
  logical,                public :: atmos_phy_mp_restart_in_aggregate
  logical,                public :: atmos_phy_mp_restart_in_postfix_timelabel  = .false.                
  character(len=H_LONG),  public :: atmos_phy_mp_restart_out_basename          = ''
  logical,                public :: atmos_phy_mp_restart_out_aggregate
  logical,                public :: atmos_phy_mp_restart_out_postfix_timelabel = .true.                 
  character(len=H_MID),   public :: atmos_phy_mp_restart_out_title             = 'ATMOS_PHY_MP restart'
  character(len=H_SHORT), public :: atmos_phy_mp_restart_out_dtype             = 'DEFAULT'
 
  real(rp), public :: atmos_phy_mp_cldfrac_thleshold
 
  real(rp), public, allocatable :: atmos_phy_mp_dens_t(:,:,:)    ! tendency DENS [kg/m3/s]
  real(rp), public, allocatable :: atmos_phy_mp_momz_t(:,:,:)    ! tendency MOMZ [kg/m2/s2]
  real(rp), public, allocatable :: atmos_phy_mp_rhou_t(:,:,:)    ! tendency dens*U [kg/m2/s2]
  real(rp), public, allocatable :: atmos_phy_mp_rhov_t(:,:,:)    ! tendency dens*V [kg/m2/s2]
  real(rp), public, allocatable :: atmos_phy_mp_rhot_t(:,:,:)    ! tendency RHOT [K*kg/m3/s]
  real(rp), public, allocatable :: atmos_phy_mp_rhoq_t(:,:,:,:)  ! tendency rho*QTRC [kg/m3/s]
  real(rp), public, allocatable :: atmos_phy_mp_rhoc_t(:,:,:,:)  ! tendency rho*QTRC (charge) [fC/m3/s]
  real(rp), public, allocatable :: atmos_phy_mp_rhoh  (:,:,:)    ! diabatic heating rate [J/kg/s]
 
  real(rp), public, allocatable :: atmos_phy_mp_evaporate(:,:,:) ! number concentration of evaporated cloud [/m3]
  real(rp), public, allocatable :: atmos_phy_mp_sflx_rain(:,:)   ! precipitation flux (liquid) [kg/m2/s]
  real(rp), public, allocatable :: atmos_phy_mp_sflx_snow(:,:)   ! precipitation flux (solid)  [kg/m2/s]
  real(rp), public, allocatable :: atmos_phy_mp_sflx_engi(:,:)   ! internal energy flux        [J/m2/s]
 
  integer, public :: qa_mp = 0
  integer, public :: qs_mp = -1
  integer, public :: qe_mp = -2
 
  !-----------------------------------------------------------------------------
  !
  !++ Private procedure
  !
  !-----------------------------------------------------------------------------
  !
  !++ Private parameters & variables
  !
  integer,                private, parameter :: vmax = 0
 
  character(len=H_SHORT), private            :: var_name(vmax)
  character(len=H_MID),   private            :: var_desc(vmax)
  character(len=H_SHORT), private            :: var_unit(vmax)
  integer,                private            :: var_id(vmax)
  integer,                private            :: restart_fid = -1  ! file ID
 
!  data VAR_NAME / /
!  data VAR_DESC / /
!  data VAR_UNIT / /
 
 
  ! for diagnostics
  real(rp), private, allocatable :: atmos_phy_mp_cldfrac(:,:,:)
  real(rp), private, allocatable :: atmos_phy_mp_re     (:,:,:,:)
  real(rp), private, allocatable :: atmos_phy_mp_qe     (:,:,:,:)
  real(rp), private, allocatable :: atmos_phy_mp_ne     (:,:,:,:)
  logical, private :: diag_cldfrac
  logical, private :: diag_re
  logical, private :: diag_qe
  logical, private :: diag_ne
 
  ! for history
  integer, private              :: hist_cldfrac_id
  logical, private              :: hist_re
  logical, private              :: hist_qe
  logical, private              :: hist_ne
  integer, private, allocatable :: hist_re_id(:)
  integer, private, allocatable :: hist_qe_id(:)
  integer, private, allocatable :: hist_ne_id(:)
 
  !-----------------------------------------------------------------------------
contains
  !-----------------------------------------------------------------------------
  subroutine atmos_phy_mp_vars_setup
    use scale_prc, only: &
       prc_abort
    use scale_const, only: &
       undef => const_undef
    use scale_atmos_hydrometeor, only: &
       n_hyd,    &
       qha,      &
       hyd_name, &
       num_name, &
       hyd_desc
    use scale_file_history, only: &
       file_history_reg
    implicit none
 
    namelist / param_atmos_phy_mp_vars / &
       atmos_phy_mp_restart_in_basename,           &
       atmos_phy_mp_restart_in_aggregate,          &
       atmos_phy_mp_restart_in_postfix_timelabel,  &
       atmos_phy_mp_restart_output,                &
       atmos_phy_mp_restart_out_basename,          &
       atmos_phy_mp_restart_out_aggregate,         &
       atmos_phy_mp_restart_out_postfix_timelabel, &
       atmos_phy_mp_restart_out_title,             &
       atmos_phy_mp_restart_out_dtype
 
    integer :: ierr
    integer :: iv, ih
    !---------------------------------------------------------------------------
 
    log_newline
    log_info("ATMOS_PHY_MP_vars_setup",*) 'Setup'
 
    allocate( atmos_phy_mp_dens_t(ka,ia,ja)    )
    allocate( atmos_phy_mp_momz_t(ka,ia,ja)    )
    allocate( atmos_phy_mp_rhou_t(ka,ia,ja)    )
    allocate( atmos_phy_mp_rhov_t(ka,ia,ja)    )
    allocate( atmos_phy_mp_rhot_t(ka,ia,ja)    )
    allocate( atmos_phy_mp_rhoq_t(ka,ia,ja,qs_mp:qe_mp) )
    allocate( atmos_phy_mp_rhoh(ka,ia,ja)    )
    allocate( atmos_phy_mp_evaporate(ka,ia,ja)    )
    ! tentative approach
 
    atmos_phy_mp_dens_t(:,:,:)   = 0.0_rp
    atmos_phy_mp_momz_t(:,:,:)   = 0.0_rp
    atmos_phy_mp_rhou_t(:,:,:)   = 0.0_rp
    atmos_phy_mp_rhov_t(:,:,:)   = 0.0_rp
    atmos_phy_mp_rhot_t(:,:,:)   = 0.0_rp
    atmos_phy_mp_rhoq_t(:,:,:,:) = 0.0_rp
    atmos_phy_mp_rhoh(:,:,:)   = 0.0_rp
    atmos_phy_mp_evaporate(:,:,:)   = 0.0_rp
    !$acc enter data copyin(ATMOS_PHY_MP_DENS_t, ATMOS_PHY_MP_MOMZ_t, ATMOS_PHY_MP_RHOU_t, ATMOS_PHY_MP_RHOV_t, ATMOS_PHY_MP_RHOT_t, ATMOS_PHY_MP_RHOQ_t, ATMOS_PHY_MP_RHOH, ATMOS_PHY_MP_EVAPORATE)
 
    allocate( atmos_phy_mp_sflx_rain(ia,ja) )
    allocate( atmos_phy_mp_sflx_snow(ia,ja) )
    allocate( atmos_phy_mp_sflx_engi(ia,ja) )
    atmos_phy_mp_sflx_rain(:,:) = undef
    atmos_phy_mp_sflx_snow(:,:) = undef
    atmos_phy_mp_sflx_engi(:,:) = undef
    !$acc enter data create(ATMOS_PHY_MP_SFLX_rain, ATMOS_PHY_MP_SFLX_snow, ATMOS_PHY_MP_SFLX_ENGI)
 
    !--- read namelist
    rewind(io_fid_conf)
    read(io_fid_conf,nml=param_atmos_phy_mp_vars,iostat=ierr)
    if( ierr < 0 ) then !--- missing
       log_info("ATMOS_PHY_MP_vars_setup",*) 'Not found namelist. Default used.'
    elseif( ierr > 0 ) then !--- fatal error
       log_error("ATMOS_PHY_MP_vars_setup",*) 'Not appropriate names in namelist PARAM_ATMOS_PHY_MP_VARS. Check!'
       call prc_abort
    endif
    log_nml(param_atmos_phy_mp_vars)
 
    log_newline
    log_info("ATMOS_PHY_MP_vars_setup",*) '[ATMOS_PHY_MP] prognostic/diagnostic variables'
    log_info_cont('(1x,A,A24,A,A48,A,A12,A)') &
               '      |', 'VARNAME                 ','|', &
               'DESCRIPTION                                     ', '[', 'UNIT        ', ']'
    do iv = 1, vmax
       log_info_cont('(1x,A,I3,A,A24,A,A48,A,A12,A)') &
                  'NO.',iv,'|',var_name(iv),'|',var_desc(iv),'[',var_unit(iv),']'
    enddo
 
    log_newline
    if ( atmos_phy_mp_restart_in_basename /= '' ) then
       log_info("ATMOS_PHY_MP_vars_setup",*) 'Restart input?  : YES, file = ', trim(atmos_phy_mp_restart_in_basename)
       log_info("ATMOS_PHY_MP_vars_setup",*) 'Add timelabel?  : ', atmos_phy_mp_restart_in_postfix_timelabel
    else
       log_info("ATMOS_PHY_MP_vars_setup",*) 'Restart input?  : NO'
    endif
    if (       atmos_phy_mp_restart_output             &
         .AND. atmos_phy_mp_restart_out_basename /= '' ) then
       log_info("ATMOS_PHY_MP_vars_setup",*) 'Restart output? : YES, file = ', trim(atmos_phy_mp_restart_out_basename)
       log_info("ATMOS_PHY_MP_vars_setup",*) 'Add timelabel?  : ', atmos_phy_mp_restart_out_postfix_timelabel
    else
       log_info("ATMOS_PHY_MP_vars_setup",*) 'Restart output? : NO'
       atmos_phy_mp_restart_output = .false.
    endif
 
 
    ! diagnostices
    allocate( atmos_phy_mp_cldfrac(ka,ia,ja) )
    allocate( atmos_phy_mp_re(ka,ia,ja,n_hyd) )
    allocate( atmos_phy_mp_qe(ka,ia,ja,n_hyd) )
    allocate( atmos_phy_mp_ne(ka,ia,ja,n_hyd) )
!OCL XFILL
    atmos_phy_mp_cldfrac(:,:,:) = undef
!OCL XFILL
    atmos_phy_mp_re(:,:,:,:) = undef
!OCL XFILL
    atmos_phy_mp_qe(:,:,:,:) = undef
    atmos_phy_mp_ne(:,:,:,:) = undef
    diag_cldfrac = .false.
    diag_re      = .false.
    diag_qe      = .false.
    diag_ne      = .false.
    !$acc enter data create(ATMOS_PHY_MP_CLDFRAC, ATMOS_PHY_MP_Re, ATMOS_PHY_MP_Qe, ATMOS_PHY_MP_Ne)
 
    ! history
    allocate( hist_re_id(n_hyd) )
    allocate( hist_qe_id(n_hyd) )
    allocate( hist_ne_id(n_hyd) )
 
    call file_history_reg( 'CLDFRAC', 'cloud fraction', '1', hist_cldfrac_id, fill_halo=.true., dim_type='ZXY' )
 
    hist_re = .false.
    do ih = 1, n_hyd
       call file_history_reg( 'Re_'//trim(hyd_name(ih)), 'effective radius of '//trim(hyd_desc(ih)), 'cm', hist_re_id(ih), fill_halo=.true., dim_type='ZXY' )
       if( hist_re_id(ih) > 0 ) hist_re = .true.
    enddo
 
    hist_qe = .false.
    do ih = 1, n_hyd
       call file_history_reg( trim(hyd_name(ih))//'_hyd', 'mass ratio of '//trim(hyd_desc(ih)), 'kg/kg', hist_qe_id(ih), fill_halo=.true., dim_type='ZXY' )
       if( hist_qe_id(ih) > 0 ) hist_qe = .true.
    enddo
 
    hist_ne = .false.
    do ih = 1, n_hyd
       call file_history_reg( trim(num_name(ih))//'_hyd', 'number concentration of '//trim(hyd_desc(ih)), '1/m3', hist_ne_id(ih), fill_halo=.true., dim_type='ZXY' )
       if( hist_ne_id(ih) > 0 ) hist_ne = .true.
    enddo
 
    return
  end subroutine atmos_phy_mp_vars_setup
 
  !-----------------------------------------------------------------------------
  subroutine atmos_phy_mp_vars_finalize
    implicit none
    !---------------------------------------------------------------------------
 
    log_newline
    log_info("ATMOS_PHY_MP_vars_finalize",*) 'Finalize'
 
    !$acc exit data delete(ATMOS_PHY_MP_DENS_t,ATMOS_PHY_MP_MOMZ_t,ATMOS_PHY_MP_RHOU_t,ATMOS_PHY_MP_RHOV_t,ATMOS_PHY_MP_RHOT_t,ATMOS_PHY_MP_RHOH,ATMOS_PHY_MP_EVAPORATE)
    deallocate( atmos_phy_mp_dens_t     )
    deallocate( atmos_phy_mp_momz_t     )
    deallocate( atmos_phy_mp_rhou_t     )
    deallocate( atmos_phy_mp_rhov_t     )
    deallocate( atmos_phy_mp_rhot_t     )
    deallocate( atmos_phy_mp_rhoq_t     )
    deallocate( atmos_phy_mp_rhoh       )
    deallocate( atmos_phy_mp_evaporate  )
 
    !$acc exit data delete(ATMOS_PHY_MP_SFLX_rain,ATMOS_PHY_MP_SFLX_snow,ATMOS_PHY_MP_SFLX_ENGI)
    deallocate( atmos_phy_mp_sflx_rain )
    deallocate( atmos_phy_mp_sflx_snow )
    deallocate( atmos_phy_mp_sflx_engi )
 
    ! diagnostices
    !$acc exit data delete(ATMOS_PHY_MP_CLDFRAC,ATMOS_PHY_MP_Re,ATMOS_PHY_MP_Qe,ATMOS_PHY_MP_Ne)
    deallocate( atmos_phy_mp_cldfrac )
    deallocate( atmos_phy_mp_re      )
    deallocate( atmos_phy_mp_qe      )
    deallocate( atmos_phy_mp_ne      )
 
 
    ! history
    deallocate( hist_re_id )
    deallocate( hist_qe_id )
    deallocate( hist_ne_id )
 
    return
  end subroutine atmos_phy_mp_vars_finalize
 
  !-----------------------------------------------------------------------------
  subroutine atmos_phy_mp_vars_fillhalo
    use scale_comm_cartesc, only: &
       comm_vars8, &
       comm_wait
    implicit none
    !---------------------------------------------------------------------------
 
    call comm_vars8( atmos_phy_mp_sflx_rain(:,:), 1 )
    call comm_vars8( atmos_phy_mp_sflx_snow(:,:), 2 )
    call comm_vars8( atmos_phy_mp_sflx_engi(:,:), 3 )
    call comm_wait ( atmos_phy_mp_sflx_rain(:,:), 1 )
    call comm_wait ( atmos_phy_mp_sflx_snow(:,:), 2 )
    call comm_wait ( atmos_phy_mp_sflx_engi(:,:), 3 )
 
    return
  end subroutine atmos_phy_mp_vars_fillhalo
 
  !-----------------------------------------------------------------------------
  subroutine atmos_phy_mp_vars_restart_open
    use scale_time, only: &
       time_gettimelabel
    use scale_file_cartesc, only: &
       file_cartesc_open
    implicit none
 
    character(len=19)     :: timelabel
    character(len=H_LONG) :: basename
    !---------------------------------------------------------------------------
 
!!$    LOG_NEWLINE
!!$    LOG_INFO("ATMOS_PHY_MP_vars_restart_open",*) 'Open restart file (ATMOS_PHY_MP) '
!!$
!!$    if ( ATMOS_PHY_MP_RESTART_IN_BASENAME /= '' ) then
!!$
!!$       if ( ATMOS_PHY_MP_RESTART_IN_POSTFIX_TIMELABEL ) then
!!$          call TIME_gettimelabel( timelabel )
!!$          basename = trim(ATMOS_PHY_MP_RESTART_IN_BASENAME)//'_'//trim(timelabel)
!!$       else
!!$          basename = trim(ATMOS_PHY_MP_RESTART_IN_BASENAME)
!!$       endif
!!$
!!$       LOG_INFO("ATMOS_PHY_MP_vars_restart_open",*) 'basename: ', trim(basename)
!!$
!!$       call FILE_CARTESC_open( basename, restart_fid, aggregate=ATMOS_PHY_MP_RESTART_IN_AGGREGATE )
!!$    else
!!$       LOG_INFO("ATMOS_PHY_MP_vars_restart_open",*) 'restart file for ATMOS_PHY_MP is not specified.'
!!$    endif
 
    return
  end subroutine atmos_phy_mp_vars_restart_open
 
  !-----------------------------------------------------------------------------
  subroutine atmos_phy_mp_vars_restart_read
    use scale_file, only: &
       file_get_aggregate
    use scale_file_cartesc, only: &
       file_cartesc_read, &
       file_cartesc_flush
    implicit none
 
    !---------------------------------------------------------------------------
 
!!$    if ( restart_fid /= -1 ) then
!!$       LOG_NEWLINE
!!$       LOG_INFO("ATMOS_PHY_MP_vars_restart_read",*) 'Read from restart file (ATMOS_PHY_MP) '
!!$
!!$       call FILE_CARTESC_read( restart_fid, VAR_NAME(1), 'XY', & ! [IN]
!!$                               ATMOS_PHY_MP_hoge(:,:)          ) ! [OUT]
!!$
!!$       if ( FILE_get_AGGREGATE(restart_fid) ) then
!!$          call FILE_CARTESC_flush( restart_fid ) ! X/Y halos have been read from file
!!$          !$acc update device(ATMOS_PHY_MP_hoge)
!!$       else
!!$          call ATMOS_PHY_MP_vars_fillhalo
!!$       end if
!!$
!!$       call ATMOS_PHY_MP_vars_check
!!$
!!$    else
!!$       LOG_INFO("ATMOS_PHY_MP_vars_restart_read",*) 'invalid restart file ID for ATMOS_PHY_MP.'
!!$    endif
 
    return
  end subroutine atmos_phy_mp_vars_restart_read
 
  !-----------------------------------------------------------------------------
  subroutine atmos_phy_mp_vars_restart_create
    use scale_time, only: &
       time_gettimelabel
    use scale_file_cartesc, only: &
       file_cartesc_create
    implicit none
 
    character(len=19)     :: timelabel
    character(len=H_LONG) :: basename
    !---------------------------------------------------------------------------
 
!!$    if ( ATMOS_PHY_MP_RESTART_OUT_BASENAME /= '' ) then
!!$
!!$       LOG_NEWLINE
!!$       LOG_INFO("ATMOS_PHY_MP_vars_restart_create",*) 'Create restart file (ATMOS_PHY_AE) '
!!$
!!$       if ( ATMOS_PHY_MP_RESTART_OUT_POSTFIX_TIMELABEL ) then
!!$          call TIME_gettimelabel( timelabel )
!!$          basename = trim(ATMOS_PHY_MP_RESTART_OUT_BASENAME)//'_'//trim(timelabel)
!!$       else
!!$          basename = trim(ATMOS_PHY_MP_RESTART_OUT_BASENAME)
!!$       endif
!!$
!!$       LOG_INFO("ATMOS_PHY_MP_vars_restart_create",*) 'basename: ', trim(basename)
!!$
!!$       call FILE_CARTESC_create( &
!!$            basename, ATMOS_PHY_MP_RESTART_OUT_TITLE, ATMOS_PHY_MP_RESTART_OUT_DTYPE, & ! [IN]
!!$            restart_fid,                                                              & ! [OUT]
!!$            aggregate=ATMOS_PHY_MP_RESTART_OUT_AGGREGATE                              ) ! [IN]
!!$
!!$    endif
 
    return
  end subroutine atmos_phy_mp_vars_restart_create
 
  !-----------------------------------------------------------------------------
  subroutine atmos_phy_mp_vars_restart_enddef
    use scale_file_cartesc, only: &
       file_cartesc_enddef
    implicit none
 
!!$    if ( restart_fid /= -1 ) then
!!$       call FILE_CARTESC_enddef( restart_fid ) ! [IN]
!!$    endif
 
    return
  end subroutine atmos_phy_mp_vars_restart_enddef
 
  !-----------------------------------------------------------------------------
  subroutine atmos_phy_mp_vars_restart_close
    use scale_file_cartesc, only: &
       file_cartesc_close
    implicit none
    !---------------------------------------------------------------------------
 
!!$    if ( restart_fid /= -1 ) then
!!$       LOG_NEWLINE
!!$       LOG_INFO("ATMOS_PHY_MP_vars_restart_close",*) 'Close restart file (ATMOS_PHY_MP) '
!!$
!!$       call FILE_CARTESC_close( restart_fid ) ! [IN]
!!$
!!$       restart_fid = -1
!!$    endif
 
    return
  end subroutine atmos_phy_mp_vars_restart_close
 
  !-----------------------------------------------------------------------------
  subroutine atmos_phy_mp_vars_restart_def_var
    use scale_file_cartesc, only: &
       file_cartesc_def_var
    implicit none
    integer :: iv
    !---------------------------------------------------------------------------
 
!!$    if ( restart_fid /= -1 ) then
!!$
!!$       do iv = 1, VMAX
!!$          call FILE_CARTESC_def_var( restart_fid, VAR_NAME(iv), VAR_DESC(iv), VAR_UNIT(iv), &
!!$                                     'XY', ATMOS_PHY_MP_RESTART_OUT_DTYPE, &
!!$                                     VAR_ID(iv) )
!!$       end do
!!$
!!$    endif
 
    return
  end subroutine atmos_phy_mp_vars_restart_def_var
 
  !-----------------------------------------------------------------------------
  subroutine atmos_phy_mp_vars_restart_write
    use scale_file_cartesc, only: &
       file_cartesc_write_var
    implicit none
    !---------------------------------------------------------------------------
 
!!$    if ( restart_fid /= -1 ) then
!!$
!!$       call ATMOS_PHY_MP_vars_fillhalo
!!$
!!$       call ATMOS_PHY_MP_vars_check
!!$
!!$       call FILE_CARTESC_write_var( restart_fid, VAR_ID(1), ATMOS_PHY_MP_hoge(:,:), &
!!$                              VAR_NAME(1), 'XY' ) ! [IN]
!!$
!!$    endif
 
    return
  end subroutine atmos_phy_mp_vars_restart_write
 
  !-----------------------------------------------------------------------------
  subroutine atmos_phy_mp_vars_history( &
       DENS, TEMP, QTRC )
    use scale_atmos_hydrometeor, only: &
       n_hyd
    use scale_file_history, only: &
       file_history_query, &
       file_history_put
    use mod_atmos_admin, only: &
       atmos_phy_mp_type
    implicit none
 
    real(rp), intent(in) :: dens(ka,ia,ja)
    real(rp), intent(in) :: temp(ka,ia,ja)
    real(rp), intent(in) :: qtrc(ka,ia,ja,qa)
 
    real(rp) :: work  (ka,ia,ja,n_hyd)
    logical  :: do_put
    integer  :: ih
    !---------------------------------------------------------------------------
 
    !$acc data create(WORK)
 
    if ( hist_cldfrac_id > 0 ) then
       call file_history_query( hist_cldfrac_id, do_put )
 
       if ( do_put ) then
          call atmos_phy_mp_vars_get_diagnostic( &
               dens(:,:,:), temp(:,:,:), qtrc(:,:,:,:), & ! [IN]
               cldfrac=work(:,:,:,1)                    ) ! [OUT]
          call file_history_put( hist_cldfrac_id, work(:,:,:,1) )
       end if
    end if
 
    if ( hist_re ) then
       do ih = 1, n_hyd
          if ( hist_re_id(ih) > 0 ) then
             call file_history_query( hist_re_id(ih), do_put )
             if ( do_put ) then
                call atmos_phy_mp_vars_get_diagnostic( &
                     dens(:,:,:), temp(:,:,:), qtrc(:,:,:,:), & ! [IN]
                     re=work(:,:,:,:)                         ) ! [OUT]
                exit
             end if
          end if
       end do
       if ( do_put ) then
          do ih = 1, n_hyd
             if ( hist_re_id(ih) > 0 ) then
                call file_history_query( hist_re_id(ih), do_put )
                if ( do_put ) then
                   call file_history_put( hist_re_id(ih), work(:,:,:,ih) )
                end if
             end if
          end do
       end if
    end if
 
    if ( hist_qe ) then
       do ih = 1, n_hyd
          if ( hist_qe_id(ih) > 0 ) then
             call file_history_query( hist_qe_id(ih), do_put )
             if ( do_put ) then
                call atmos_phy_mp_vars_get_diagnostic( &
                     dens(:,:,:), temp(:,:,:), qtrc(:,:,:,:), & ! [IN]
                     qe=work(:,:,:,:)                         ) ! [OUT]
                exit
             end if
          end if
       end do
       if ( do_put ) then
          do ih = 1, n_hyd
             if ( hist_qe_id(ih) > 0 ) then
                call file_history_query( hist_qe_id(ih), do_put )
                if( do_put ) then
                   call file_history_put( hist_qe_id(ih), work(:,:,:,ih) )
                end if
             end if
          end do
       end if
    end if
 
    if ( hist_ne ) then
       do ih = 1, n_hyd
          if ( hist_ne_id(ih) > 0 ) then
             call file_history_query( hist_ne_id(ih), do_put )
             if ( do_put ) then
                call atmos_phy_mp_vars_get_diagnostic( &
                     dens(:,:,:), temp(:,:,:), qtrc(:,:,:,:), & ! [IN]
                     ne=work(:,:,:,:)                         ) ! [OUT]
                exit
             end if
          end if
       end do
       if ( do_put ) then
          do ih = 1, n_hyd
             if ( hist_ne_id(ih) > 0 ) then
                call file_history_query( hist_ne_id(ih), do_put )
                if( do_put ) then
                   call file_history_put( hist_ne_id(ih), work(:,:,:,ih) )
                end if
             end if
          end do
       end if
    end if
 
    !$acc end data
 
    return
  end subroutine atmos_phy_mp_vars_history
 
  subroutine atmos_phy_mp_vars_get_diagnostic( &
       DENS, TEMP, QTRC,   &
       CLDFRAC, Re, Qe, Ne )
    use scale_atmos_hydrometeor, only: &
       n_hyd, &
       i_hc,  &
       i_hr,  &
       i_hi,  &
       i_hs,  &
       i_hg,  &
       i_hh
    use scale_atmos_phy_mp_kessler, only: &
       atmos_phy_mp_kessler_qtrc2qhyd, &
       atmos_phy_mp_kessler_effective_radius, &
       atmos_phy_mp_kessler_cloud_fraction
    use scale_atmos_phy_mp_tomita08, only: &
       atmos_phy_mp_tomita08_qtrc2qhyd, &
       atmos_phy_mp_tomita08_effective_radius, &
       atmos_phy_mp_tomita08_cloud_fraction
    use scale_atmos_phy_mp_sn14, only: &
       atmos_phy_mp_sn14_qtrc2qhyd, &
       atmos_phy_mp_sn14_qtrc2nhyd, &
       atmos_phy_mp_sn14_effective_radius, &
       atmos_phy_mp_sn14_cloud_fraction
    use scale_atmos_phy_mp_suzuki10, only: &
       atmos_phy_mp_suzuki10_qtrc2qhyd, &
       atmos_phy_mp_suzuki10_qtrc2nhyd, &
       atmos_phy_mp_suzuki10_effective_radius, &
       atmos_phy_mp_suzuki10_cloud_fraction
    use mod_atmos_admin, only: &
       atmos_phy_mp_type
    implicit none
 
    real(rp), intent(in) :: dens(ka,ia,ja)
    real(rp), intent(in) :: temp(ka,ia,ja)
    real(rp), intent(in) :: qtrc(ka,ia,ja,qa)
    real(rp), intent(out), optional :: cldfrac(ka,ia,ja)
    real(rp), intent(out), optional :: re     (ka,ia,ja,n_hyd)
    real(rp), intent(out), optional :: qe     (ka,ia,ja,n_hyd)
    real(rp), intent(out), optional :: ne     (ka,ia,ja,n_hyd)
 
    integer :: k, i, j, ih
 
    if ( present(cldfrac) ) then
       !$acc data copyout(CLDFRAC)
       if ( .not. diag_cldfrac ) then
          select case ( atmos_phy_mp_type )
          case ( 'KESSLER' )
             !$acc update host(QTRC(:,:,:,QS_MP+1:QE_MP))
             call atmos_phy_mp_kessler_cloud_fraction( &
                  ka, ks, ke, ia, is, ie, ja, js, je, &
                  qtrc(:,:,:,qs_mp+1:qe_mp), atmos_phy_mp_cldfrac_thleshold, & ! [IN]
                  atmos_phy_mp_cldfrac(:,:,:)                                ) ! [OUT]
             !$acc update device(ATMOS_PHY_MP_CLDFRAC)
          case ( 'TOMITA08' )
             call atmos_phy_mp_tomita08_cloud_fraction( &
                  ka, ks, ke, ia, isb, ieb, ja, jsb, jeb, &
                  qtrc(:,:,:,qs_mp+1:qe_mp), atmos_phy_mp_cldfrac_thleshold, & ! [IN]
                  atmos_phy_mp_cldfrac(:,:,:)                                ) ! [OUT]
          case ( 'SN14' )
             !$acc update host(QTRC(:,:,:,QS_MP+1:QE_MP))
             call atmos_phy_mp_sn14_cloud_fraction( &
                  ka, ks, ke, ia, isb, ieb, ja, jsb, jeb, &
                  qtrc(:,:,:,qs_mp+1:qe_mp), atmos_phy_mp_cldfrac_thleshold, & ! [IN]
                  atmos_phy_mp_cldfrac(:,:,:)                                ) ! [OUT]
             !$acc update device(ATMOS_PHY_MP_CLDFRAC)
          case ( 'SUZUKI10' )
             !$acc update host(QTRC(:,:,:,QS_MP+1:QE_MP))
             call atmos_phy_mp_suzuki10_cloud_fraction( &
                  ka, ks, ke, ia, isb, ieb, ja, jsb, jeb, &
                  qtrc(:,:,:,qs_mp+1:qe_mp), atmos_phy_mp_cldfrac_thleshold, & ! [IN]
                  atmos_phy_mp_cldfrac(:,:,:)                                ) ! [OUT]
             !$acc update device(ATMOS_PHY_MP_CLDFRAC)
          case default
!OCL XFILL
             !$acc kernels
             atmos_phy_mp_cldfrac(:,:,:) = 0.0_rp
             !$acc end kernels
          end select
          diag_cldfrac = .true.
       end if
!OCL XFILL
       !$acc kernels
       do j = jsb, jeb
       do i = isb, ieb
       do k = ks, ke
          cldfrac(k,i,j) = atmos_phy_mp_cldfrac(k,i,j)
       end do
       end do
       end do
       !$acc end kernels
 
       !$acc end data
    end if
 
    if ( present(re) ) then
       !$acc data copyout(Re)
       if ( .not. diag_re ) then
          select case ( atmos_phy_mp_type )
          case ( 'KESSLER' )
             !$acc update host(QTRC(:,:,:,QS_MP+1:QE_MP))
             call atmos_phy_mp_kessler_effective_radius( &
                  ka, ks, ke, ia, is, ie, ja, js, je, &
                  dens(:,:,:), temp(:,:,:), qtrc(:,:,:,qs_mp+1:qe_mp), & ! [IN]
                  atmos_phy_mp_re(:,:,:,:)                             ) ! [OUT]
             !$acc update device(ATMOS_PHY_MP_Re)
          case ( 'TOMITA08' )
             call atmos_phy_mp_tomita08_effective_radius( &
                  ka, ks, ke, ia, isb, ieb, ja, jsb, jeb, &
                  dens(:,:,:), temp(:,:,:), qtrc(:,:,:,qs_mp+1:qe_mp), & ! [IN]
                  atmos_phy_mp_re(:,:,:,:)                             ) ! [OUT]
          case ( 'SN14' )
             !$acc update host(QTRC(:,:,:,QS_MP+1:QE_MP))
             call atmos_phy_mp_sn14_effective_radius( &
                  ka, ks, ke, ia, isb, ieb, ja, jsb, jeb, &
                  dens(:,:,:), temp(:,:,:), qtrc(:,:,:,qs_mp+1:qe_mp), & ! [IN]
                  atmos_phy_mp_re(:,:,:,:)                             ) ! [OUT]
             !$acc update device(ATMOS_PHY_MP_Re)
          case ( 'SUZUKI10' )
             !$acc update host(QTRC(:,:,:,QS_MP+1:QE_MP))
             call atmos_phy_mp_suzuki10_effective_radius( &
                  ka, ks, ke, ia, isb, ieb, ja, jsb, jeb, &
                  dens(:,:,:), temp(:,:,:), qtrc(:,:,:,qs_mp+1:qe_mp), & ! [IN]
                  atmos_phy_mp_re(:,:,:,:)                             ) ! [OUT]
             !$acc update device(ATMOS_PHY_MP_Re)
          case default
!OCL XFILL
             !$acc kernels
             atmos_phy_mp_re(:,:,:,:) = 0.0_rp
             !$acc end kernels
          end select
          diag_re = .true.
       end if
!OCL XFILL
       !$acc kernels
       do ih = 1, n_hyd
       do j = jsb, jeb
       do i = isb, ieb
       do k = ks, ke
          re(k,i,j,ih) = atmos_phy_mp_re(k,i,j,ih)
       end do
       end do
       end do
       end do
       !$acc end kernels
 
       !$acc end data
    end if
 
    if ( present(qe) ) then
       !$acc data copyout(Qe)
       if ( .not. diag_qe ) then
          select case ( atmos_phy_mp_type )
          case ( 'KESSLER' )
             !$acc update host(QTRC(:,:,:,QS_MP+1:QE_MP))
             call atmos_phy_mp_kessler_qtrc2qhyd( &
                  ka, ks, ke, ia, is, ie, ja, js, je, &
                  qtrc(:,:,:,qs_mp+1:qe_mp), & ! [IN]
                  atmos_phy_mp_qe(:,:,:,:)   ) ! [OUT]
             !$acc update device(ATMOS_PHY_MP_Qe)
          case ( 'TOMITA08' )
             call atmos_phy_mp_tomita08_qtrc2qhyd( &
                  ka, ks, ke, ia, isb, ieb, ja, jsb, jeb, &
                  qtrc(:,:,:,qs_mp+1:qe_mp), & ! [IN]
                  atmos_phy_mp_qe(:,:,:,:)   ) ! [OUT]
          case ( 'SN14' )
             !$acc update host(QTRC(:,:,:,QS_MP+1:QE_MP))
             call atmos_phy_mp_sn14_qtrc2qhyd( &
                  ka, ks, ke, ia, isb, ieb, ja, jsb, jeb, &
                  qtrc(:,:,:,qs_mp+1:qe_mp), & ! [IN]
                  atmos_phy_mp_qe(:,:,:,:)   ) ! [OUT]
             !$acc update device(ATMOS_PHY_MP_Qe)
          case ( 'SUZUKI10' )
             !$acc update host(QTRC(:,:,:,QS_MP+1:QE_MP))
             call atmos_phy_mp_suzuki10_qtrc2qhyd( &
                  ka, ks, ke, ia, isb, ieb, ja, jsb, jeb, &
                  qtrc(:,:,:,qs_mp+1:qe_mp), & ! [IN]
                  atmos_phy_mp_qe(:,:,:,:)   ) ! [OUT]
             !$acc update device(ATMOS_PHY_MP_Qe)
          case default
!OCL XFILL
             !$acc kernels
             atmos_phy_mp_qe(:,:,:,:) = 0.0_rp
             !$acc end kernels
          end select
          diag_qe = .true.
       end if
!OCL XFILL
       !$acc kernels
       do ih = 1, n_hyd
       do j = jsb, jeb
       do i = isb, ieb
       do k = ks, ke
          qe(k,i,j,ih) = atmos_phy_mp_qe(k,i,j,ih)
       end do
       end do
       end do
       end do
       !$acc end kernels
 
       !$acc end data
    end if
 
    if ( present(ne) ) then
       !$acc data copyout(Ne)
       if ( .not. diag_ne ) then
          select case ( atmos_phy_mp_type )
          case ( 'KESSLER', 'TOMITA08' )
             ! do nothing
          case ( 'SN14' )
             !$acc update host(QTRC(:,:,:,QS_MP+1:QE_MP))
             call atmos_phy_mp_sn14_qtrc2nhyd( &
                  ka, ks, ke, ia, isb, ieb, ja, jsb, jeb, &
                  qtrc(:,:,:,qs_mp+1:qe_mp), & ! [IN]
                  atmos_phy_mp_ne(:,:,:,:)   ) ! [OUT]
             !$acc update device(ATMOS_PHY_MP_Ne)
          case ( 'SUZUKI10' )
             !$acc update host(QTRC(:,:,:,QS_MP+1:QE_MP))
             call atmos_phy_mp_suzuki10_qtrc2nhyd( &
                  ka, ks, ke, ia, isb, ieb, ja, jsb, jeb, &
                  dens(:,:,:), qtrc(:,:,:,qs_mp+1:qe_mp), & ! [IN]
                  atmos_phy_mp_ne(:,:,:,:)                ) ! [OUT]
             !$acc update device(ATMOS_PHY_MP_Ne)
          end select
          diag_ne = .true.
       end if
!OCL XFILL
       !$acc kernels
       do ih = 1, n_hyd
       do j = jsb, jeb
       do i = isb, ieb
       do k = ks, ke
          ne(k,i,j,ih) = atmos_phy_mp_ne(k,i,j,ih)
       end do
       end do
       end do
       end do
       !$acc end kernels
 
       !$acc end data
    end if
 
    return
  end subroutine atmos_phy_mp_vars_get_diagnostic
 
  subroutine atmos_phy_mp_vars_reset_diagnostics
    diag_cldfrac = .false.
    diag_re      = .false.
    diag_qe      = .false.
    diag_ne      = .false.
 
    return
  end subroutine atmos_phy_mp_vars_reset_diagnostics
 
  subroutine atmos_phy_mp_vars_check
    use scale_statistics, only: &
       statistics_total
    use scale_atmos_grid_cartesc_real, only: &
       atmos_grid_cartesc_real_area,    &
       atmos_grid_cartesc_real_totarea, &
       atmos_grid_cartesc_real_vol,     &
       atmos_grid_cartesc_real_totvol
    implicit none
 
!!$    call VALCHECK( IA, IS, IE, JA, JS, JE, &
!!$                   ATMOS_PHY_MP_hoge(:,:),      & ! (in)
!!$                   0.0_RP, 0.0_RP, VAR_NAME(1), & ! (in)
!!$                   __FILE__, __LINE__           ) ! (in)
!!$
!!$    call STATISTICS_total( IA, IS, IE, JA, JS, JE, &
!!$                           ATMOS_PHY_MP_hoge(:,:), VAR_NAME(1), & ! (in)
!!$                           ATMOS_GRID_CARTESC_REAL_AREA(:,:),   & ! (in)
!!$                           ATMOS_GRID_CARTESC_REAL_TOTAREA      ) ! (in)
 
    return
  end subroutine atmos_phy_mp_vars_check
 
 
end module mod_atmos_phy_mp_vars