scale_history.F90

Go to the documentation of this file.

!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------
module scale_history
  !-----------------------------------------------------------------------------
  !
  !++ Used modules
  !
  use scale_precision
  use scale_stdio
  use scale_prof
  use scale_grid_index
  use scale_land_grid_index
  use scale_urban_grid_index
  !-----------------------------------------------------------------------------
  implicit none
  private
  !-----------------------------------------------------------------------------
  !
  !++ Public procedures
  !
  public :: hist_setup
  public :: hist_switch
  public :: hist_setpres
  public :: hist_reg
  public :: hist_query
  public :: hist_put
  public :: hist_in
  public :: hist_get
  public :: hist_write

  interface hist_put
     module procedure hist_put_0d
     module procedure hist_put_1d
     module procedure hist_put_2d
     module procedure hist_put_3d
  end interface hist_put

  interface hist_in
     module procedure hist_in_0d
     module procedure hist_in_1d
     module procedure hist_in_2d
     module procedure hist_in_3d
  end interface hist_in

  interface hist_get
     module procedure hist_get_1d
     module procedure hist_get_2d
     module procedure hist_get_3d
  end interface hist_get

  !-----------------------------------------------------------------------------
  !
  !++ Public parameters & variables
  !
  !-----------------------------------------------------------------------------
  !
  !++ Private procedures
  !
  private :: hist_put_axes

  !-----------------------------------------------------------------------------
  !
  !++ Private parameters & variables
  !
  integer,                private, parameter   :: I_MODEL = 0
  integer,                private, parameter   :: I_Z     = 1
  integer,                private, parameter   :: I_PRES  = 2

  integer,                private              :: HIST_item_limit
  integer,                private              :: HIST_item_count
  character(len=H_SHORT), private, allocatable :: HIST_item   (:)
  integer,                private, allocatable :: HIST_variant(:)
  integer,                private, allocatable :: HIST_zcoord (:,:)

  integer,                private, parameter   :: HIST_PRES_nlim   = 300
  integer,                private              :: HIST_PRES_nlayer =  -1
  real(RP),               private, allocatable :: HIST_PRES_val(:)

  logical,                private              :: enabled
  integer,                private              :: im,  jm,  km
  integer,                private              :: ims, ime
  integer,                private              :: jms, jme

  logical, private :: HIST_BND = .false.
  !-----------------------------------------------------------------------------
contains
  !-----------------------------------------------------------------------------
  subroutine hist_setup
    use gtool_history, only: &
       historyinit
    use scale_process, only: &
       prc_mpistop,    &
       prc_masterrank, &
       prc_myrank
    use scale_rm_process, only: &
       prc_2drank, &
       prc_num_x, &
       prc_num_y
    use scale_time, only: &
       time_nowdate,     &
       time_nowms,       &
       time_dtsec,       &
       time_startdaysec
    use scale_interpolation, only: &
       interp_setup_pres
    implicit none

    character(len=H_MID) :: history_h_title = 'SCALE-RM HISTORY OUTPUT'
    character(len=H_MID) :: history_t_since

    real(RP)             :: hist_pres(hist_pres_nlim) = 0.0_rp          

    namelist / param_hist / &
       hist_pres_nlayer, &
       hist_pres,        &
       hist_bnd

    integer  :: hist_variant_limit

    real(DP) :: start_daysec
    integer  :: rankidx(2)
    integer  :: ierr
    integer  :: k
    !---------------------------------------------------------------------------

    if( io_l ) write(io_fid_log,*)
    if( io_l ) write(io_fid_log,*) '++++++ Module[HISTORY] / Categ[ATMOS-RM IO] / Origin[SCALElib]'

    !--- read namelist
    rewind(io_fid_conf)
    read(io_fid_conf,nml=param_hist,iostat=ierr)
    if( ierr < 0 ) then !--- missing
       if( io_l ) write(io_fid_log,*) '*** Not found namelist. Default used.'
    elseif( ierr > 0 ) then !--- fatal error
       write(*,*) 'xxx Not appropriate names in namelist PARAM_HIST. Check!'
       call prc_mpistop
    endif
    if( io_nml ) write(io_fid_nml,nml=param_hist)

    ! check pressure coordinate
    if ( hist_pres_nlayer > 0 ) then
       if ( hist_pres_nlayer > 100 ) then
          write(*,*) 'xxx number of layers of pressure is larger the KMAX'
          call prc_mpistop
       end if

       do k = 1, hist_pres_nlayer
          if ( hist_pres(k) <= 0.0_rp ) then
             write(*,*) 'xxx Invalid value found in pressure coordinate! (k,value)=', k, hist_pres(k)
             call prc_mpistop
          elseif ( hist_pres(k+1) >= hist_pres(k) ) then
             write(*,*) 'xxx The value of pressure coordinate must be descending order! ', &
                        '(k,value[k],value[k+1])=', k, hist_pres(k), hist_pres(k+1)
             call prc_mpistop
          endif
       enddo
    else
       if( io_l ) write(io_fid_log,*)
       if( io_l ) write(io_fid_log,*) '*** HIST_PRES_nlayer is not set.'
       if( io_l ) write(io_fid_log,*) '*** Output with pressure coordinate is disabled'
    endif

    call prof_rapstart('FILE_O_NetCDF', 2)

    rankidx(1) = prc_2drank(prc_myrank,1)
    rankidx(2) = prc_2drank(prc_myrank,2)

    start_daysec = time_startdaysec
    if ( time_nowdate(1) > 0 ) then
       write(history_t_since,'(I4.4,5(A1,I2.2))')      time_nowdate(1), &
                                                  '-', time_nowdate(2), &
                                                  '-', time_nowdate(3), &
                                                  ' ', time_nowdate(4), &
                                                  ':', time_nowdate(5), &
                                                  ':', time_nowdate(6)
       start_daysec = time_nowms
    else
       history_t_since = ''
    endif

    if ( hist_bnd ) then
       im  = imaxb
       jm  = jmaxb
       ims = isb
       ime = ieb
       jms = jsb
       jme = jeb
    else
       im  = imax
       jm  = jmax
       ims = is
       ime = ie
       jms = js
       jme = je
    end if

    km = max( lkmax, ukmax, kmax )

    call historyinit( hist_item_limit,                  & ! [OUT]
                      hist_variant_limit,               & ! [OUT]
                      im, jm, km,                       & ! [IN]
                      prc_masterrank,                   & ! [IN]
                      prc_myrank,                       & ! [IN]
                      rankidx(:),                       & ! [IN]
                      history_h_title,                  & ! [IN]
                      h_source,                         & ! [IN]
                      h_institute,                      & ! [IN]
                      time_start     = start_daysec,    & ! [IN]
                      time_interval  = time_dtsec,      & ! [IN]
                      time_since     = history_t_since, & ! [IN]
                      default_zcoord = 'model',         & ! [IN]
                      namelist_fid   = io_fid_conf      ) ! [IN]

    hist_item_count = 0
    if ( hist_item_limit > 0 ) then
       allocate( hist_item(hist_item_limit)                    )
       allocate( hist_variant(hist_item_limit)                    )
       allocate( hist_zcoord(hist_item_limit,hist_variant_limit) )
       hist_item(:)   = ''
       hist_variant(:)   = 0
       hist_zcoord(:,:) = 0
    endif

    if ( hist_pres_nlayer > 0 ) then
       allocate( hist_pres_val(hist_pres_nlayer) )

       do k = 1, hist_pres_nlayer
          hist_pres_val(k) = hist_pres(k) * 100.0_rp ! [hPa->Pa]
       enddo

       call interp_setup_pres( hist_pres_nlayer ) ! [IN]
    endif

    call hist_put_axes

    enabled = .true.

    call prof_rapend  ('FILE_O_NetCDF', 2)

    return
  end subroutine hist_setup

  !-----------------------------------------------------------------------------
  !  only register the axis and coordinate variables into internal buffers
  !  The actual write happens later when calling HIST_write
  subroutine hist_put_axes
    use gtool_history, only: &
       historyputaxis,  &
       historyputassociatedcoordinates
    use scale_const, only: &
       d2r => const_d2r
    use scale_grid, only: &
       grid_cz,    &
       grid_cx,    &
       grid_cy,    &
       grid_fz,    &
       grid_fx,    &
       grid_fy,    &
       grid_cdz,   &
       grid_cdx,   &
       grid_cdy,   &
       grid_fdz,   &
       grid_fdx,   &
       grid_fdy,   &
       grid_cbfz,  &
       grid_cbfx,  &
       grid_cbfy,  &
       grid_fbfz,  &
       grid_fbfx,  &
       grid_fbfy,  &
       grid_cxg,   &
       grid_cyg,   &
       grid_fxg,   &
       grid_fyg,   &
       grid_cdxg,  &
       grid_cdyg,  &
       grid_fdxg,  &
       grid_fdyg,  &
       grid_cbfxg, &
       grid_cbfyg, &
       grid_fbfxg, &
       grid_fbfyg
    use scale_land_grid, only: &
       grid_lcz, &
       grid_lfz, &
       grid_lcdz
    use scale_urban_grid, only: &
       grid_ucz, &
       grid_ufz, &
       grid_ucdz
    use scale_grid_real, only: &
       real_cz,   &
       real_fz,   &
       real_lon,  &
       real_lonx, &
       real_lony, &
       real_lonxy, &
       real_lat,  &
       real_latx, &
       real_laty, &
       real_latxy
    use scale_topography, only: &
       topo_zsfc
    use scale_landuse, only: &
       landuse_frac_land
    use scale_process, only: &
       prc_myrank
    use scale_rm_process, only: &
       prc_2drank, &
       prc_num_x, &
       prc_num_y
    implicit none

    real(RP)         :: axis     (im,jm,kmax)
    character(len=2) :: axis_name(3)

    integer :: k, i, j
    integer :: rankidx(2)
    integer :: start(3)
    integer :: startx, starty, startz
    integer :: xag, yag
    !---------------------------------------------------------------------------

    rankidx(1) = prc_2drank(prc_myrank,1)
    rankidx(2) = prc_2drank(prc_myrank,2)

    ! For parallel I/O, some variables are written by a subset of processes.
    ! 1. Only PRC_myrank 0 writes all z axes
    ! 2. Only south-most processes (rankidx(2) == 0) write x axes
    !         rankidx(1) == 0           writes west HALO
    !         rankidx(1) == PRC_NUM_X-1 writes east HALO
    !         others                    writes without HALO
    ! 3. Only west-most processes (rankidx(1) == 0) write y axes
    !         rankidx(1) == 0           writes south HALO
    !         rankidx(1) == PRC_NUM_Y-1 writes north HALO
    !         others                    writes without HALO

    if ( io_aggregate ) then
       if ( hist_bnd ) then
          startx = isgb   ! global subarray starting index
          starty = jsgb   ! global subarray starting index
          xag    = iagb
          yag    = jagb
       else
          startx = 1 + prc_2drank(prc_myrank,1) * imax    ! no IHALO
          starty = 1 + prc_2drank(prc_myrank,2) * jmax    ! no JHALO
          xag    = imaxg
          yag    = jmaxg
       end if
       startz = 1
       if ( rankidx(2) .GT. 0 ) startx = -1   ! only south-most processes write
       if ( rankidx(1) .GT. 0 ) starty = -1   ! only west-most processes write
       if ( prc_myrank .GT. 0 ) startz = -1   ! only rank 0 writes Z axes
    else
       xag = im
       yag = jm
       startx = 1
       starty = 1
       startz = 1
    end if

    ! for the shared-file I/O method, the axes are global (gsize)
    ! for one-file-per-process I/O method, the axes size is equal to the local buffer size
    call historyputaxis( 'x',   'X',               'm', 'x',   grid_cx(ims:ime), gsize=xag,  start=startx )
    call historyputaxis( 'y',   'Y',               'm', 'y',   grid_cy(jms:jme), gsize=yag,  start=starty )
    call historyputaxis( 'z',   'Z',               'm', 'z',   grid_cz(ks:ke),   gsize=kmax, start=startz )
    call historyputaxis( 'xh',  'X (half level)',  'm', 'xh',  grid_fx(ims:ime), gsize=xag,  start=startx )
    call historyputaxis( 'yh',  'Y (half level)',  'm', 'yh',  grid_fy(jms:jme), gsize=yag,  start=starty )
    call historyputaxis( 'zh',  'Z (half level)',  'm', 'zh',  grid_fz(ks:ke),   gsize=kmax, start=startz )

    if ( hist_pres_nlayer > 0 ) then
       call historyputaxis( 'pressure', 'Pressure', 'hPa', 'pressure',         &
                            hist_pres_val(:)/100.0_rp, down=.true., gsize=hist_pres_nlayer, start=startz )
    endif

    call historyputaxis( 'lz',  'LZ',              'm', 'lz',  grid_lcz(lks:lke), down=.true., gsize=lkmax, start=startz )
    call historyputaxis( 'lzh', 'LZ (half level)', 'm', 'lzh', grid_lfz(lks:lke), down=.true., gsize=lkmax, start=startz )

    call historyputaxis( 'uz',  'UZ',              'm', 'uz',  grid_ucz(uks:uke), down=.true., gsize=ukmax, start=startz )
    call historyputaxis( 'uzh', 'UZ (half level)', 'm', 'uzh', grid_ufz(uks:uke), down=.true., gsize=ukmax, start=startz )

    ! axes below always include halos when written to file regardless of PRC_PERIODIC_X/PRC_PERIODIC_Y
    call historyputaxis( 'CZ',  'Atmos Grid Center Position Z', 'm', 'CZ',  grid_cz,  gsize=ka,    start=startz )
    call historyputaxis( 'FZ',  'Atmos Grid Face Position Z',   'm', 'FZ',  grid_fz,  gsize=ka+1,  start=startz )
    call historyputaxis( 'CDZ',  'Grid Cell length Z', 'm', 'CZ',  grid_cdz,  gsize=ka,    start=startz )
    call historyputaxis( 'FDZ',  'Grid distance Z',    'm', 'FDZ', grid_fdz,  gsize=ka-1,  start=startz )
    if ( io_aggregate ) then
       call historyputaxis( 'CX',  'Atmos Grid Center Position X', 'm', 'CX',  grid_cxg, gsize=iag,   start=startz )
       call historyputaxis( 'CY',  'Atmos Grid Center Position Y', 'm', 'CY',  grid_cyg, gsize=jag,   start=startz )
       call historyputaxis( 'FX',  'Atmos Grid Face Position X',   'm', 'FX',  grid_fxg, gsize=iag+1, start=startz )
       call historyputaxis( 'FY',  'Atmos Grid Face Position Y',   'm', 'FY',  grid_fyg, gsize=jag+1, start=startz )
       call historyputaxis( 'CDX',  'Grid Cell length X', 'm', 'CX',  grid_cdxg, gsize=iag,   start=startz )
       call historyputaxis( 'CDY',  'Grid Cell length Y', 'm', 'CY',  grid_cdyg, gsize=jag,   start=startz )
       call historyputaxis( 'FDX',  'Grid distance X',    'm', 'FDX', grid_fdxg, gsize=iag-1, start=startz )
       call historyputaxis( 'FDY',  'Grid distance Y',    'm', 'FDY', grid_fdyg, gsize=jag-1, start=startz )
    else
       call historyputaxis( 'CX',  'Atmos Grid Center Position X', 'm', 'CX',  grid_cx )
       call historyputaxis( 'CY',  'Atmos Grid Center Position Y', 'm', 'CY',  grid_cy )
       call historyputaxis( 'FX',  'Atmos Grid Face Position X',   'm', 'FX',  grid_fx )
       call historyputaxis( 'FY',  'Atmos Grid Face Position Y',   'm', 'FY',  grid_fy )
       call historyputaxis( 'CDX',  'Grid Cell length X', 'm', 'CX',  grid_cdx )
       call historyputaxis( 'CDY',  'Grid Cell length Y', 'm', 'CY',  grid_cdy )
       call historyputaxis( 'FDX',  'Grid distance X',    'm', 'FDX', grid_fdx )
       call historyputaxis( 'FDY',  'Grid distance Y',    'm', 'FDY', grid_fdy )
    end if

    call historyputaxis( 'LCZ',  'Land Grid Center Position Z', 'm', 'LCZ', grid_lcz, down=.true., gsize=lkmax,   start=startz )
    call historyputaxis( 'LFZ',  'Land Grid Face Position Z',   'm', 'LFZ', grid_lfz, down=.true., gsize=lkmax+1, start=startz )
    call historyputaxis( 'LCDZ', 'Land Grid Cell length Z',     'm', 'LCZ', grid_lcdz,             gsize=lkmax,   start=startz )

    call historyputaxis( 'UCZ',  'Urban Grid Center Position Z', 'm', 'UCZ', grid_ucz, down=.true., gsize=ukmax,   start=startz )
    call historyputaxis( 'UFZ',  'Urban Grid Face Position Z',   'm', 'UFZ', grid_ufz, down=.true., gsize=ukmax+1, start=startz )
    call historyputaxis( 'UCDZ', 'Urban Grid Cell length Z',     'm', 'UCZ', grid_ucdz,             gsize=ukmax,   start=startz )

    call historyputaxis('CBFZ',  'Boundary factor Center Z', '1', 'CZ', grid_cbfz,  gsize=ka,  start=startz )
    call historyputaxis('FBFZ',  'Boundary factor Face Z',   '1', 'CZ', grid_fbfz,  gsize=ka,  start=startz )
    if ( io_aggregate ) then
       call historyputaxis('CBFX',  'Boundary factor Center X', '1', 'CX', grid_cbfxg, gsize=iag, start=startz )
       call historyputaxis('CBFY',  'Boundary factor Center Y', '1', 'CY', grid_cbfyg, gsize=jag, start=startz )
       call historyputaxis('FBFX',  'Boundary factor Face X',   '1', 'CX', grid_fbfxg, gsize=iag, start=startz )
       call historyputaxis('FBFY',  'Boundary factor Face Y',   '1', 'CY', grid_fbfyg, gsize=jag, start=startz )
    else
       call historyputaxis('CBFX',  'Boundary factor Center X', '1', 'CX', grid_cbfx )
       call historyputaxis('CBFY',  'Boundary factor Center Y', '1', 'CY', grid_cbfy )
       call historyputaxis('FBFX',  'Boundary factor Face X',   '1', 'CX', grid_fbfx )
       call historyputaxis('FBFY',  'Boundary factor Face Y',   '1', 'CY', grid_fbfy )
    end if

    ! TODO: skip 8 axes below when IO_AGGREGATE is true, as all axes are now global
    call historyputaxis('CXG',   'Grid Center Position X (global)', 'm', 'CXG', grid_cxg, gsize=iag,   start=startz )
    call historyputaxis('CYG',   'Grid Center Position Y (global)', 'm', 'CYG', grid_cyg, gsize=jag,   start=startz )
    call historyputaxis('FXG',   'Grid Face Position X (global)',   'm', 'FXG', grid_fxg, gsize=iag+1, start=startz )
    call historyputaxis('FYG',   'Grid Face Position Y (global)',   'm', 'FYG', grid_fyg, gsize=jag+1, start=startz )

    call historyputaxis('CBFXG', 'Boundary factor Center X (global)', '1', 'CXG', grid_cbfxg, gsize=iag, start=startz )
    call historyputaxis('CBFYG', 'Boundary factor Center Y (global)', '1', 'CYG', grid_cbfyg, gsize=jag, start=startz )
    call historyputaxis('FBFXG', 'Boundary factor Face X (global)',   '1', 'CXG', grid_fbfxg, gsize=iag, start=startz )
    call historyputaxis('FBFYG', 'Boundary factor Face Y (global)',   '1', 'CYG', grid_fbfyg, gsize=jag, start=startz )

    ! associate coordinates
    if ( io_aggregate ) then
       if ( hist_bnd ) then
          start(1) = isgb   ! global subarray starting index
          start(2) = jsgb   ! global subarray starting index
       else
          start(1) = 1 + prc_2drank(prc_myrank,1) * imax    ! no IHALO
          start(2) = 1 + prc_2drank(prc_myrank,2) * jmax    ! no JHALO
       end if
       start(3) = 1
    else
       start(:) = 1
    end if

    do k = 1, kmax
       axis(1:im,1:jm,k) = real_cz(k+ks-1,ims:ime,jms:jme)
    enddo
    axis_name(1:3) = (/'x ', 'y ', 'z '/)
    call historyputassociatedcoordinates( 'height', 'height above ground level', &
                                          'm', axis_name(1:3), axis(:,:,:), start=start )

    do k = 1, kmax
       axis(1:im,1:jm,k) = real_fz(k+ks-1,ims:ime,jms:jme)
    enddo
    axis_name(1:3) = (/'x ', 'y ', 'zh'/)
    call historyputassociatedcoordinates( 'height_xyw', 'height above ground level (half level xyw)', &
                                          'm' , axis_name(1:3), axis(:,:,:), start=start )

    do k = 1, kmax
    do j = 1, jm
    do i = 1, min(im,ia-ims)
       axis(i,j,k) = ( real_cz(k+ks-1,ims+i-1,jms+j-1) + real_cz(k+ks-1,ims+i,jms+j-1) ) * 0.5_rp
    enddo
    enddo
    enddo
    if ( im == ia-ims+1 ) then
       do k = 1, kmax
       do j = 1, jm
          axis(im,j,k) = real_cz(k+ks-1,ims+im-1,jms+j-1)
       enddo
       enddo
    end if
    axis_name(1:3) = (/'xh', 'y ', 'z '/)
    call historyputassociatedcoordinates( 'height_uyz', 'height above ground level (half level uyz)', &
                                          'm', axis_name(1:3), axis(:,:,:), start=start )

    do k = 1, kmax
    do j = 1, min(jm,ja-jms)
    do i = 1, im
       axis(i,j,k) = ( real_cz(k+ks-1,ims+i-1,jms+j-1) + real_cz(k+ks-1,ims+i-1,jms+j) ) * 0.5_rp
    enddo
    enddo
    enddo
    if ( jm == ja-jms+1 ) then
       do k = 1, kmax
       do i = 1, im
          axis(i,jm,k) = real_cz(k+ks-1,ims+i-1,jms+jm-1)
       enddo
       enddo
    end if
    axis_name(1:3) = (/'x ', 'yh', 'z '/)
    call historyputassociatedcoordinates( 'height_xvz', 'height above ground level (half level xvz)', &
                                          'm', axis_name(1:3), axis(:,:,:), start=start )

    do k = 1, kmax
    do j = 1, min(jm,ja-jms)
    do i = 1, min(im,ia-ims)
       axis(i,j,k) = ( real_cz(k+ks-1,ims+i-1,jms+j-1) + real_cz(k+ks-1,ims+i  ,jms+j-1) &
                     + real_cz(k+ks-1,ims+i-1,jms+j  ) + real_cz(k+ks-1,ims+i  ,jms+j  ) ) * 0.25_rp
    enddo
    enddo
    enddo
    if ( jm == ja-jms+1 ) then
       do k = 1, kmax
       do i = 1, min(im,ia-ims)
          axis(i,jm,k) = ( real_cz(k+ks-1,ims+i-1,jms+jm-1) + real_cz(k+ks-1,ims+i,jms+jm-1) ) * 0.5_rp
       enddo
       enddo
    end if
    if ( im == ia-ims+1 ) then
       do k = 1, kmax
       do j = 1, min(jm,ja-jms)
          axis(im,j,k) = ( real_cz(k+ks-1,ims+im-1,jms+j-1) + real_cz(k+ks-1,ims+im-1,jms+j) ) * 0.5_rp
       enddo
       enddo
    end if
    if ( im == ia-ims+1 .and. jm == ja-jms+1 ) then
       do k = 1, kmax
          axis(im,jm,k) = real_cz(k+ks-1,ims+im-1,jms+jm-1)
       enddo
    end if
    axis_name(1:3) = (/'xh', 'yh', 'z '/)
    call historyputassociatedcoordinates( 'height_uvz', 'height above ground level (half level uvz)', &
                                          'm', axis_name(1:3), axis(:,:,:), start=start )

    do k = 1, kmax
    do j = 1, jm
    do i = 1, min(im,ia-ims)
       axis(i,j,k) = ( real_fz(k+ks-1,ims+i-1,jms+j-1) + real_fz(k+ks-1,ims+i,jms+j-1) ) * 0.5_rp
    enddo
    enddo
    enddo
    if ( im == ia-ims+1 ) then
       do k = 1, kmax
       do j = 1, jm
          axis(im,j,k) = real_fz(k+ks-1,ims+im-1,jms+j-1)
       enddo
       enddo
    end if
    axis_name(1:3) = (/'xh', 'y ', 'zh'/)
    call historyputassociatedcoordinates( 'height_uyw', 'height above ground level (half level uyw)', &
                                          'm', axis_name(1:3), axis(:,:,:), start=start )

    do k = 1, kmax
    do j = 1, min(jm,ja-jms)
    do i = 1, im
       axis(i,j,k) = ( real_fz(k+ks-1,ims+i-1,jms+j-1) + real_fz(k+ks-1,ims+i-1,jms+j) ) * 0.5_rp
    enddo
    enddo
    enddo
    if ( jm == ja-jms+1 ) then
       do k = 1, kmax
       do i = 1, im
          axis(i,jm,k) = real_fz(k+ks-1,ims+i-1,jms+jm-1)
       enddo
       enddo
    end if
    axis_name(1:3) = (/'x ', 'yh', 'zh'/)
    call historyputassociatedcoordinates( 'height_xvw', 'height above ground level (half level xvw)', &
                                          'm', axis_name(1:3), axis(:,:,:), start=start )

    do k = 1, kmax
    do j = 1, min(jm,ja-jms)
    do i = 1, min(im,ia-ims)
       axis(i,j,k) = ( real_fz(k+ks-1,ims+i-1,jms+j-1) + real_fz(k+ks-1,ims+i  ,jms+j-1) &
                     + real_fz(k+ks-1,ims+i-1,jms+j  ) + real_fz(k+ks-1,ims+i  ,jms+j  ) ) * 0.25_rp
    enddo
    enddo
    enddo
    if ( jm == ja-jms+1 ) then
       do k = 1, kmax
       do i = 1, min(im,ia-ims)
          axis(i,jm,k) = ( real_fz(k+ks-1,ims+i-1,jms+jm-1) + real_fz(k+ks-1,ims+i,jms+jm-1) ) * 0.5_rp
       enddo
       enddo
    end if
    if ( im == ia-ims+1 ) then
       do k = 1, kmax
       do j = 1, min(jm,ja-jms)
          axis(im,j,k) = ( real_fz(k+ks-1,ims+im-1,jms+j-1) + real_fz(k+ks-1,ims+im-1,jms+j) ) * 0.5_rp
       enddo
       enddo
    end if
    if ( im == ia-ims+1 .and. jm == ja-jms+1 ) then
       do k = 1, kmax
          axis(im,jm,k) = real_fz(k+ks-1,ims+im-1,jms+jm-1)
       enddo
    end if
    axis_name(1:3) = (/'xh', 'yh', 'zh'/)
    call historyputassociatedcoordinates( 'height_uvw', 'height above ground level (half level uvw)', &
                                          'm', axis_name(1:3), axis(:,:,:), start=start               )

    axis(1:im,1:jm,1) = real_lon(ims:ime,jms:jme) / d2r
    axis_name(1:2) = (/'x ', 'y '/)
    call historyputassociatedcoordinates( 'lon', 'longitude',                                      &
                                          'degrees_east', axis_name(1:2), axis(:,:,1), start=start )

    axis(1:im,1:jm,1) = real_lonx(ims:ime,jms:jme) / d2r
    axis_name(1:2) = (/'xh', 'y '/)
    call historyputassociatedcoordinates( 'lon_uy', 'longitude (half level uy)',                   &
                                          'degrees_east', axis_name(1:2), axis(:,:,1), start=start )

    axis(1:im,1:jm,1) = real_lony(ims:ime,jms:jme) / d2r
    axis_name(1:2) = (/'x ', 'yh'/)
    call historyputassociatedcoordinates( 'lon_xv', 'longitude (half level xv)',                   &
                                          'degrees_east', axis_name(1:2), axis(:,:,1), start=start )

    axis(1:im,1:jm,1) = real_lonxy(ims:ime,jms:jme) / d2r
    axis_name(1:2) = (/'xh', 'yh'/)
    call historyputassociatedcoordinates( 'lon_uv', 'longitude (half level uv)',                   &
                                          'degrees_east', axis_name(1:2), axis(:,:,1), start=start )

    axis(1:im,1:jm,1) = real_lat(ims:ime,jms:jme) / d2r
    axis_name(1:2) = (/'x ', 'y '/)
    call historyputassociatedcoordinates( 'lat', 'latitude',                                        &
                                          'degrees_north', axis_name(1:2), axis(:,:,1), start=start )

    axis(1:im,1:jm,1) = real_latx(ims:ime,jms:jme) / d2r
    axis_name(1:2) = (/'xh', 'y '/)
    call historyputassociatedcoordinates( 'lat_uy', 'latitude (half level uy)',                     &
                                          'degrees_north', axis_name(1:2), axis(:,:,1), start=start )

    axis(1:im,1:jm,1) = real_laty(ims:ime,jms:jme) / d2r
    axis_name(1:2) = (/'x ', 'yh'/)
    call historyputassociatedcoordinates( 'lat_xv', 'latitude (half level xv)',                     &
                                          'degrees_north', axis_name(1:2), axis(:,:,1), start=start )

    axis(1:im,1:jm,1) = real_latxy(ims:ime,jms:jme) / d2r
    axis_name(1:2) = (/'xh', 'yh'/)
    call historyputassociatedcoordinates( 'lat_uv', 'latitude (half level uv)',                     &
                                          'degrees_north', axis_name(1:2), axis(:,:,1), start=start )

    axis(1:im,1:jm,1) = topo_zsfc(ims:ime,jms:jme)
    axis_name(1:2) = (/'x ', 'y '/)
    call historyputassociatedcoordinates( 'topo', 'topography',                         &
                                          'm', axis_name(1:2), axis(:,:,1), start=start )

    axis(1:im,1:jm,1) = landuse_frac_land(ims:ime,jms:jme)
    axis_name(1:2) = (/'x ', 'y '/)
    call historyputassociatedcoordinates( 'lsmask', 'fraction for land-sea mask',       &
                                          '1', axis_name(1:2), axis(:,:,1), start=start )

    return
  end subroutine hist_put_axes

  !-----------------------------------------------------------------------------
  subroutine hist_switch( switch )
    implicit none

    logical, intent(in) :: switch
    !---------------------------------------------------------------------------

    enabled = switch

    return
  end subroutine hist_switch

  !-----------------------------------------------------------------------------
  subroutine hist_setpres( &
       PRES,    &
       SFC_PRES )
    use scale_interpolation, only: &
       interp_update_pres
    implicit none

    real(RP), intent(in) :: pres    (ka,ia,ja) ! pressure in Xi coordinate [Pa]
    real(RP), intent(in) :: sfc_pres(   ia,ja) ! surface pressure          [Pa]
    !---------------------------------------------------------------------------

    if ( hist_pres_nlayer > 0 ) then
       call interp_update_pres( hist_pres_nlayer,     & ! [IN]
                                pres(:,:,:), & ! [IN]
                                sfc_pres(:,:)  , & ! [IN]
                                hist_pres_val(:)      ) ! [IN]
    endif

    return
  end subroutine hist_setpres

  !-----------------------------------------------------------------------------
  subroutine hist_reg( &
       itemid, &
       item,   &
       desc,   &
       unit,   &
       ndim,   &
       xdim,   &
       ydim,   &
       zdim    )
    use mpi, only: &
       mpi_comm_null
    use gtool_history, only: &
       historyaddvariable, &
       historycheck
    use scale_time, only: &
       nowstep => time_nowstep
    use scale_process, only: &
       prc_mpistop,          &
       prc_myrank,           &
       prc_local_comm_world
    use scale_rm_process, only: &
       prc_2drank, &
       prc_num_x, &
       prc_num_y
    implicit none

    integer,          intent(out) :: itemid
    character(len=*), intent(in)  :: item
    character(len=*), intent(in)  :: desc
    character(len=*), intent(in)  :: unit
    integer,          intent(in)  :: ndim
    character(len=*), intent(in), optional :: xdim
    character(len=*), intent(in), optional :: ydim
    character(len=*), intent(in), optional :: zdim

    logical :: flag_half_x
    logical :: flag_half_y
    logical :: flag_half_z

    character(len=H_SHORT) :: dims(3)

    integer                :: nvariant1, nvariant2, nvariant3
    integer                :: v, id
    logical                :: atom

    integer                :: rankidx(2)
    integer                :: start(4), count(4)
    integer                :: comm
    !---------------------------------------------------------------------------

    itemid = -1

    if( .NOT. enabled ) return

    if( hist_item_limit == 0 ) return

    do id = 1, hist_item_count
       if ( item == hist_item(id) ) then ! item exists
          itemid = id
          return
       endif
    enddo

    call prof_rapstart('FILE_O_NetCDF', 2)

    ! Try to add new item

    if ( len_trim(item) >= h_short ) then
       write(*,'(1x,A,I2,A,A)') 'xxx Length of history name should be <= ', h_short-1 ,' chars. STOP', trim(item)
       call prc_mpistop
    endif

    atom = .true.

    rankidx(1) = prc_2drank(prc_myrank,1)
    rankidx(2) = prc_2drank(prc_myrank,2)

    start = 0
    count = 0

    if ( ndim == 1 ) then

       ! check half/full level for vertical
       dims(1) = "z"
       if ( present(zdim) ) then
          if ( zdim == 'half' ) then
             dims(1) = "zh"
          endif
       endif

       ! for shared-file parallel I/O, only rank 0 writes variables with only Z dimension
       start(1) = 1
       count(1) = kmax
       if ( prc_myrank .GT. 0 ) count(1) = 0

    elseif ( ndim == 2 ) then

       ! check half/full level for horizontal
       flag_half_x = .false.
       if ( present(xdim) ) then
          if( xdim == 'half' ) flag_half_x = .true.
       endif

       flag_half_y = .false.
       if ( present(ydim) ) then
          if( ydim == 'half' ) flag_half_y = .true.
       endif

       if    ( flag_half_x .AND. flag_half_y ) then
          dims(1) = 'lon_uv'
          dims(2) = 'lat_uv'
       elseif( flag_half_x ) then
          dims(1) = 'lon_uy'
          dims(2) = 'lat_uy'
       elseif( flag_half_y ) then
          dims(1) = 'lon_xv'
          dims(2) = 'lat_xv'
       else
          dims(1) = 'lon'
          dims(2) = 'lat'
       endif

    elseif ( ndim == 3 ) then

       ! check half/full level for vertical/horizontal
       flag_half_x = .false.

       if ( present(xdim) ) then
          if( xdim == 'half' ) flag_half_x = .true.
       endif

       flag_half_y = .false.
       if ( present(ydim) ) then
          if( ydim == 'half' ) flag_half_y = .true.
       endif

       flag_half_z = .false.
       if ( present(zdim) ) then
          if( zdim == 'half' ) flag_half_z = .true.
       endif

       if    ( flag_half_x .AND. flag_half_y ) then
          dims(1) = 'lon_uv'
          dims(2) = 'lat_uv'
          if ( flag_half_z ) then
             dims(3) = 'height_uvw'
          else
             dims(3) = 'height_uvz'
          endif
       elseif( flag_half_x ) then
          dims(1) = 'lon_uy'
          dims(2) = 'lat_uy'
          if ( flag_half_z ) then
             dims(3) = 'height_uyw'
          else
             dims(3) = 'height_uyz'
          endif
       elseif( flag_half_y ) then
          dims(1) = 'lon_xv'
          dims(2) = 'lat_xv'
          if ( flag_half_z ) then
             dims(3) = 'height_xvw'
          else
             dims(3) = 'height_xvz'
          endif
       else
          dims(1) = 'lon'
          dims(2) = 'lat'
          if ( flag_half_z ) then
             dims(3) = 'height_xyw'
          else
             dims(3) = 'height'
          endif
       endif

       ! start and count will be used by PnetCDF I/O
       start(3) = 1
       count(3) = kmax

       if ( present(zdim) ) then
          if    ( zdim == 'land'      ) then
             dims(3) = 'lz'
             count(3) = lkmax
             atom = .false.
          elseif( zdim == 'landhalf'  ) then
             dims(3) = 'lzh'
             count(3) = lkmax
             atom = .false.
          elseif( zdim == 'urban'     ) then
             dims(3) = 'uz'
             count(3) = ukmax
             atom = .false.
          elseif( zdim == 'urbanhalf' ) then
             dims(3) = 'uzh'
             count(3) = ukmax
             atom = .false.
          endif

       endif

    endif

    if ( ndim >= 2 ) then
       ! start and count will be used by PnetCDF I/O
       if ( hist_bnd ) then
          start(1) = isgb
          start(2) = jsgb
          count(1) = imaxb
          count(2) = jmaxb
       else
          ! for the case the shared-file contains no halos
          start(1) = 1 + prc_2drank(prc_myrank,1) * imax    ! no IHALO
          start(2) = 1 + prc_2drank(prc_myrank,2) * jmax    ! no JHALO
          count(1) = imax
          count(2) = jmax
       end if
    end if


    if ( io_aggregate ) then  ! user input parameter indicates to do PnetCDF I/O
       comm = prc_local_comm_world
    else
       comm = mpi_comm_null
    end if

    if ( atom ) then

       ! model coordinate (terrain following coordinate)
       call historyaddvariable( nvariant1,    & ! [OUT]
                                item,         & ! [IN]
                                dims(1:ndim), & ! [IN]
                                desc,         & ! [IN]
                                unit,         & ! [IN]
                                nowstep,      & ! [IN]
                                'model',      & ! [IN]
                                start=start,  & ! [IN]
                                count=count,  & ! [IN]
                                comm=comm     ) ! [IN]

       ! absolute height coordinate
       dims(3) = 'z'
       call historyaddvariable( nvariant2,    & ! [OUT]
                                item,         & ! [IN]
                                dims(1:ndim), & ! [IN]
                                desc,         & ! [IN]
                                unit,         & ! [IN]
                                nowstep,      & ! [IN]
                                'z',          & ! [IN]
                                start=start,  & ! [IN]
                                count=count,  & ! [IN]
                                comm=comm     ) ! [IN]

       ! pressure coordinate
       if ( hist_pres_nlayer > 0 ) then

          dims(3) = 'pressure'
          call historyaddvariable( nvariant3,    & ! [OUT]
                                   item,         & ! [IN]
                                   dims(1:ndim), & ! [IN]
                                   desc,         & ! [IN]
                                   unit,         & ! [IN]
                                   nowstep,      & ! [IN]
                                   'pressure',   & ! [IN]
                                   start=start,  & ! [IN]
                                   count=count,  & ! [IN]
                                   comm=comm     ) ! [IN]

       else
          nvariant3 = 0
       endif

    else

       call historyaddvariable( nvariant1,    & ! [OUT]
                                item,         & ! [IN]
                                dims(1:ndim), & ! [IN]
                                desc,         & ! [IN]
                                unit,         & ! [IN]
                                nowstep,      & ! [IN]
                                start=start,  & ! [IN]
                                count=count,  & ! [IN]
                                comm=comm     ) ! [IN]

       nvariant2 = 0
       nvariant3 = 0
    endif

    if ( nvariant1 + nvariant2 + nvariant3 > 0 ) then
       hist_item_count   = hist_item_count + 1
       itemid            = hist_item_count
       hist_item(itemid) = item

       do v = 1, nvariant1
          hist_variant(itemid)                      = hist_variant(itemid) + 1
          hist_zcoord(itemid,hist_variant(itemid)) = i_model
       enddo

       do v = 1, nvariant2
          hist_variant(itemid)                      = hist_variant(itemid) + 1
          hist_zcoord(itemid,hist_variant(itemid)) = i_z
       enddo

       do v = 1, nvariant3
          hist_variant(itemid)                      = hist_variant(itemid) + 1
          hist_zcoord(itemid,hist_variant(itemid)) = i_pres
       enddo
    endif

    call prof_rapend  ('FILE_O_NetCDF', 2)

    return
  end subroutine hist_reg

  !-----------------------------------------------------------------------------
  subroutine hist_query( &
       itemid, &
       answer  )
    use gtool_history, only: &
       historyquery
    use scale_time, only: &
       time_nowstep
    implicit none

    integer,  intent(in)  :: itemid
    logical,  intent(out) :: answer
    !---------------------------------------------------------------------------

    answer = .false.

    if( .NOT. enabled ) return

    if( itemid < 0 ) return

    call prof_rapstart('FILE_O_NetCDF', 2)

    call historyquery( hist_item(itemid), & ! [IN]
                       time_nowstep,      & ! [IN]
                       answer             ) ! [OUT]

    call prof_rapend  ('FILE_O_NetCDF', 2)

    return
  end subroutine hist_query

  !-----------------------------------------------------------------------------
  subroutine hist_put_0d( &
       itemid, &
       var     )
    use gtool_history, only: &
       historyput
    use scale_time, only: &
       time_nowstep
    implicit none

    integer,  intent(in)  :: itemid
    real(RP), intent(in)  :: var

    integer  :: n, v, id
    !---------------------------------------------------------------------------

    if( .NOT. enabled ) return

    if( itemid < 0 ) return

    call prof_rapstart('FILE_O_NetCDF', 2)

    id = 0
    do n = 1, itemid-1
       id = id + hist_variant(n)
    enddo

    do v = 1, hist_variant(itemid)
       id = id + 1

       call historyput( id,           & ! [IN]
                        time_nowstep, & ! [IN]
                        var           ) ! [IN]
    enddo

    call prof_rapend  ('FILE_O_NetCDF', 2)

    return
  end subroutine hist_put_0d

  !-----------------------------------------------------------------------------
  subroutine hist_put_1d( &
       itemid, &
       var     )
    use gtool_history, only: &
       historyput
    use scale_time, only: &
       time_nowstep
    implicit none

    integer,  intent(in)  :: itemid
    real(RP), intent(in)  :: var(:)

    real(RP) :: var_trim(KMAX)

    integer  :: n, v, id
    integer  :: k
    !---------------------------------------------------------------------------

    if( .NOT. enabled ) return

    if( itemid < 0 ) return

    call prof_rapstart('FILE_O_NetCDF', 2)

    do k = 1, kmax
       var_trim(k) = var(ks+k-1)
    enddo

    id = 0
    do n = 1, itemid-1
       id = id + hist_variant(n)
    enddo

    do v = 1, hist_variant(itemid)
       id = id + 1

       call historyput( id,           & ! [IN]
                        time_nowstep, & ! [IN]
                        var_trim(:)   ) ! [IN]
    enddo

    call prof_rapend  ('FILE_O_NetCDF', 2)

    return
  end subroutine hist_put_1d

  !-----------------------------------------------------------------------------
  subroutine hist_put_2d( &
       itemid, &
       var,    &
       nohalo  )
    use gtool_file, only: &
       rmiss
    use gtool_history, only: &
       historyput
    use scale_time, only: &
       time_nowstep
    implicit none

    integer,  intent(in)  :: itemid
    real(RP), intent(in)  :: var(:,:)
    logical,  intent(in), optional :: nohalo

    real(RP) :: var_trim(im*jm)
    logical  :: nohalo_

    integer  :: n, v, id
    integer  :: i, j
    !---------------------------------------------------------------------------

    if( .NOT. enabled ) return

    if( itemid < 0 ) return

    call prof_rapstart('FILE_O_NetCDF', 2)

    do j = 1, jm
    do i = 1, im
       var_trim((j-1)*im+i) = var(ims+i-1,jms+j-1)
    enddo
    enddo

    nohalo_ = .false.
    if( present(nohalo) ) nohalo_ = nohalo

    if ( nohalo_ ) then
       ! W halo
       do j = 1, jm
       do i = 1, is-ims
          var_trim((j-1)*im+i) = rmiss
       enddo
       enddo
       ! E halo
       do j = 1, jm
       do i = ie-ims+2, ime-ims+1
          var_trim((j-1)*im+i) = rmiss
       enddo
       enddo
       ! S halo
       do j = 1, js-jms
       do i = 1, im
          var_trim((j-1)*im+i) = rmiss
       enddo
       enddo
       ! N halo
       do j = je-jms+2, jme-jms+1
       do i = 1, im
          var_trim((j-1)*im+i) = rmiss
       enddo
       enddo
    endif

    id = 0
    do n = 1, itemid-1
       id = id + hist_variant(n)
    enddo

    do v = 1, hist_variant(itemid)
       id = id + 1

       call historyput( id,           & ! [IN]
                        time_nowstep, & ! [IN]
                        var_trim(:)   ) ! [IN]
    enddo

    call prof_rapend  ('FILE_O_NetCDF', 2)

    return
  end subroutine hist_put_2d

  !-----------------------------------------------------------------------------
  subroutine hist_put_3d( &
       itemid, &
       var,    &
       xdim,   &
       ydim,   &
       zdim,   &
       nohalo  )
    use gtool_file, only: &
       rmiss
    use gtool_history, only: &
       historyput
    use scale_time, only: &
       time_nowstep
    use scale_interpolation, only: &
       interp_vertical_xi2z, &
       interp_vertical_xi2p, &
       interp_available
    implicit none

    integer,          intent(in)  :: itemid
    real(RP),         intent(in)  :: var(:,:,:)
    character(len=*), intent(in), optional :: xdim
    character(len=*), intent(in), optional :: ydim
    character(len=*), intent(in), optional :: zdim
    logical,          intent(in), optional :: nohalo

    character(len=H_SHORT) :: xd, yd, zd
    integer                :: isize, jsize, ksize
    integer                :: istart, jstart, kstart

    real(RP) :: var_Z(KA              ,IA,JA)
    real(RP) :: var_P(HIST_PRES_nlayer,IA,JA)

    real(RP) :: var_trim(km*im*jm)
    logical  :: nohalo_
    integer  :: s(3)

    integer  :: n, v, id
    integer  :: i, j, k

    intrinsic shape
    !---------------------------------------------------------------------------

    if( .NOT. enabled ) return

    if( itemid < 0 ) return

    call prof_rapstart('FILE_O_NetCDF', 2)

    xd = ''
    yd = ''
    zd = ''
    if( present(xdim) ) xd = xdim
    if( present(ydim) ) yd = ydim
    if( present(zdim) ) zd = zdim

    nohalo_ = .false.
    if( present(nohalo) ) nohalo_ = nohalo

    ! select dimension
    select case( xd )
      case('half')
        isize  = im
        istart = ims
      case default
        isize  = im
        istart = ims
    end select

    select case( yd )
      case('half')
        jsize  = jm
        jstart = jms
      case default
        jsize  = jm
        jstart = jms
    end select

    select case( zd )
      case('land')
        ksize  = lkmax
        kstart = lks
      case('landhalf')
        ksize  = lkmax
        kstart = lks
      case('urban')
        ksize  = ukmax
        kstart = uks
      case('urbanhalf')
        ksize  = ukmax
        kstart = uks
      case('half')
        ksize  = kmax
        kstart = ks
      case default
        ksize  = kmax
        kstart = ks
    end select

    s(:) = shape(var)

    id = 0
    do n = 1, itemid-1
       id = id + hist_variant(n)
    enddo

    do v = 1, hist_variant(itemid)

       if    (       s(1)  == ka                  &
               .AND. ksize == kmax                &
               .AND. hist_zcoord(itemid,v) == i_z &
               .AND. interp_available             ) then ! z*->z interpolation

          call prof_rapstart('FILE_O_interp', 2)
          call interp_vertical_xi2z( var(:,:,:), & ! [IN]
                                     var_z(:,:,:)  ) ! [OUT]
          call prof_rapend  ('FILE_O_interp', 2)

          do k = 1, ksize
          do j = 1, jsize
          do i = 1, isize
             var_trim((k-1)*jsize*isize+(j-1)*isize+i) = var_z(kstart+k-1,istart+i-1,jstart+j-1)
          enddo
          enddo
          enddo

       elseif(       s(1)  == ka                     &
               .AND. ksize == kmax                   &
               .AND. hist_zcoord(itemid,v) == i_pres ) then ! z*->p interpolation

          ksize = hist_pres_nlayer

          call prof_rapstart('FILE_O_interp', 2)
          call interp_vertical_xi2p( hist_pres_nlayer, & ! [IN]
                                     var(:,:,:),     & ! [IN]
                                     var_p(:,:,:)      ) ! [OUT]
          call prof_rapend  ('FILE_O_interp', 2)

          do k = 1, ksize
          do j = 1, jsize
          do i = 1, isize
             var_trim((k-1)*jsize*isize+(j-1)*isize+i) = var_p(k,istart+i-1,jstart+j-1)
          enddo
          enddo
          enddo

       else ! no interpolation

          do k = 1, ksize
          do j = 1, jsize
          do i = 1, isize
             var_trim((k-1)*jsize*isize+(j-1)*isize+i) = var(kstart+k-1,istart+i-1,jstart+j-1)
          enddo
          enddo
          enddo

       endif

       if ( nohalo_ ) then
             ! W halo
          do k = 1, ksize
          do j = 1, jsize
          do i = 1, is-istart
             var_trim((k-1)*jsize*isize+(j-1)*isize+i) = rmiss
          enddo
          enddo
          enddo
          ! E halo
          do k = 1, ksize
          do j = 1, jsize
          do i = ie-istart+2, ime-istart+1
             var_trim((k-1)*jsize*isize+(j-1)*isize+i) = rmiss
          enddo
          enddo
          enddo
          ! S halo
          do k = 1, ksize
          do j = 1, js-jstart
          do i = 1, isize
             var_trim((k-1)*jsize*isize+(j-1)*isize+i) = rmiss
          enddo
          enddo
          enddo
          ! N halo
          do k = 1, ksize
          do j = je-jstart+2, jme-jstart+1
          do i = 1, isize
             var_trim((k-1)*jsize*isize+(j-1)*isize+i) = rmiss
          enddo
          enddo
          enddo
       endif

       id = id + 1

       call historyput( id,                           & ! [IN]
                        time_nowstep,                 & ! [IN]
                        var_trim(1:isize*jsize*ksize) ) ! [IN]
    enddo

    call prof_rapend  ('FILE_O_NetCDF', 2)

    return
  end subroutine hist_put_3d

  !-----------------------------------------------------------------------------
  subroutine hist_in_0d( &
       var,  &
       item, &
       desc, &
       unit )
    implicit none

    real(RP),         intent(in) :: var
    character(len=*), intent(in) :: item
    character(len=*), intent(in) :: desc
    character(len=*), intent(in) :: unit

    integer :: itemid
    logical :: do_put
    !---------------------------------------------------------------------------

    if( .NOT. enabled ) return

    ! Check whether the item has been already registered
    call hist_reg  ( itemid,  & ! [OUT]
                     item,    & ! [IN]
                     desc,    & ! [IN]
                     unit,    & ! [IN]
                     0        ) ! [IN]

    ! Check whether it is time to input the item
    call hist_query( itemid,  & ! [IN]
                     do_put   ) ! [OUT]

    if ( do_put ) then
       call hist_put( itemid, & ! [IN]
                      var     ) ! [IN]
    endif

    return
  end subroutine hist_in_0d

  !-----------------------------------------------------------------------------
  subroutine hist_in_1d( &
       var,  &
       item, &
       desc, &
       unit, &
       zdim  )
    implicit none

    real(RP),         intent(in)  :: var(:)
    character(len=*), intent(in)  :: item
    character(len=*), intent(in)  :: desc
    character(len=*), intent(in)  :: unit
    character(len=*), intent(in), optional :: zdim

    character(len=H_SHORT) :: zd

    integer :: itemid
    logical :: do_put
    !---------------------------------------------------------------------------

    if( .NOT. enabled ) return

    zd = ''
    if( present(zdim) ) zd = zdim

    ! Check whether the item has been already registered
    call hist_reg  ( itemid,  & ! [OUT]
                     item,    & ! [IN]
                     desc,    & ! [IN]
                     unit,    & ! [IN]
                     1,       & ! [IN]
                     zdim=zd  ) ! [IN]

    ! Check whether it is time to input the item
    call hist_query( itemid,  & ! [IN]
                     do_put   ) ! [OUT]

    if ( do_put ) then
       call hist_put( itemid, & ! [IN]
                      var     ) ! [IN]
    endif

    return
  end subroutine hist_in_1d

  !-----------------------------------------------------------------------------
  subroutine hist_in_2d( &
       var,   &
       item,  &
       desc,  &
       unit,  &
       xdim,  &
       ydim,  &
       nohalo )
    implicit none

    real(RP),         intent(in)  :: var(:,:)
    character(len=*), intent(in)  :: item
    character(len=*), intent(in)  :: desc
    character(len=*), intent(in)  :: unit
    character(len=*), intent(in), optional :: xdim
    character(len=*), intent(in), optional :: ydim
    logical,          intent(in), optional :: nohalo

    character(len=H_SHORT) :: xd, yd

    integer :: itemid
    logical :: do_put
    !---------------------------------------------------------------------------

    if( .NOT. enabled ) return

    xd = ''
    yd = ''
    if( present(xdim) ) xd = xdim
    if( present(ydim) ) yd = ydim

    ! Check whether the item has been already registered
    call hist_reg  ( itemid,  & ! [OUT]
                     item,    & ! [IN]
                     desc,    & ! [IN]
                     unit,    & ! [IN]
                     2,       & ! [IN]
                     xdim=xd, & ! [IN]
                     ydim=yd  ) ! [IN]

    ! Check whether it is time to input the item
    call hist_query( itemid,  & ! [IN]
                     do_put   ) ! [OUT]

    if ( do_put ) then
       call hist_put( itemid,       & ! [IN]
                      var,          & ! [IN]
                      nohalo=nohalo ) ! [IN]
    endif

    return
  end subroutine hist_in_2d

  !-----------------------------------------------------------------------------
  subroutine hist_in_3d( &
       var,   &
       item,  &
       desc,  &
       unit,  &
       xdim,  &
       ydim,  &
       zdim,  &
       nohalo )
    implicit none

    real(RP),         intent(in)  :: var(:,:,:)
    character(len=*), intent(in)  :: item
    character(len=*), intent(in)  :: desc
    character(len=*), intent(in)  :: unit
    character(len=*), intent(in), optional :: xdim
    character(len=*), intent(in), optional :: ydim
    character(len=*), intent(in), optional :: zdim
    logical,          intent(in), optional :: nohalo

    character(len=H_SHORT) :: xd, yd, zd

    integer :: itemid
    logical :: do_put
    !---------------------------------------------------------------------------

    if( .NOT. enabled ) return

    xd = ''
    yd = ''
    zd = ''
    if( present(xdim) ) xd = xdim
    if( present(ydim) ) yd = ydim
    if( present(zdim) ) zd = zdim

    ! Check whether the item has been already registered
    call hist_reg  ( itemid,  & ! [OUT]
                     item,    & ! [IN]
                     desc,    & ! [IN]
                     unit,    & ! [IN]
                     3,       & ! [IN]
                     xdim=xd, & ! [IN]
                     ydim=yd, & ! [IN]
                     zdim=zd  ) ! [IN]

    ! Check whether it is time to input the item
    call hist_query( itemid,  & ! [IN]
                     do_put   ) ! [OUT]

    if ( do_put ) then
       call hist_put( itemid,       & ! [IN]
                      var,          & ! [IN]
                      xdim=xd,      & ! [IN]
                      ydim=yd,      & ! [IN]
                      zdim=zd,      & ! [IN]
                      nohalo=nohalo ) ! [IN]
    endif

    return
  end subroutine hist_in_3d

  !-----------------------------------------------------------------------------
  subroutine hist_get_1d( &
       var,          &
       basename,     &
       varname,      &
       step,         &
       allow_missing )
    use gtool_history, only: &
       historyget
    implicit none

    real(RP),         intent(out) :: var(:)
    character(len=*), intent(in)  :: basename
    character(len=*), intent(in)  :: varname
    integer,          intent(in)  :: step
    logical,          intent(in), optional :: allow_missing

    logical :: am
    !---------------------------------------------------------------------------

    call prof_rapstart('FILE_I_NetCDF', 2)

    am = .false.
    if( present(allow_missing) ) am = allow_missing

    call historyget( var(:),          & ! [OUT]
                     basename,        & ! [IN]
                     varname,         & ! [IN]
                     step,            & ! [IN]
                     allow_missing=am ) ! [IN]

    call prof_rapend  ('FILE_I_NetCDF', 2)

    return
  end subroutine hist_get_1d

  !-----------------------------------------------------------------------------
  subroutine hist_get_2d( &
       var,          &
       basename,     &
       varname,      &
       step,         &
       allow_missing )
    use gtool_history, only: &
       historyget
    implicit none

    real(RP),         intent(out) :: var(:,:)
    character(len=*), intent(in)  :: basename
    character(len=*), intent(in)  :: varname
    integer,          intent(in)  :: step
    logical,          intent(in), optional :: allow_missing

    logical :: am
    !---------------------------------------------------------------------------

    call prof_rapstart('FILE_I_NetCDF', 2)

    am = .false.
    if( present(allow_missing) ) am = allow_missing

    call historyget( var(:,:),        & ! [OUT]
                     basename,        & ! [IN]
                     varname,         & ! [IN]
                     step,            & ! [IN]
                     allow_missing=am ) ! [IN]

    call prof_rapend  ('FILE_I_NetCDF', 2)

    return
  end subroutine hist_get_2d

  !-----------------------------------------------------------------------------
  subroutine hist_get_3d( &
       var,          &
       basename,     &
       varname,      &
       step,         &
       allow_missing )
    use gtool_history, only: &
       historyget
    implicit none

    real(RP),         intent(out) :: var(:,:,:)
    character(len=*), intent(in)  :: basename
    character(len=*), intent(in)  :: varname
    integer,          intent(in)  :: step
    logical,          intent(in), optional :: allow_missing

    logical :: am
    !---------------------------------------------------------------------------

    call prof_rapstart('FILE_I_NetCDF', 2)

    am = .false.
    if( present(allow_missing) ) am = allow_missing

    call historyget( var(:,:,:),      & ! [OUT]
                     basename,        & ! [IN]
                     varname,         & ! [IN]
                     step,            & ! [IN]
                     allow_missing=am ) ! [IN]

    call prof_rapend  ('FILE_I_NetCDF', 2)

    return
  end subroutine hist_get_3d

  !-----------------------------------------------------------------------------
  subroutine hist_write
    use gtool_history, only: &
       historywriteall
    use scale_time, only: &
       time_nowstep
    implicit none
    !---------------------------------------------------------------------------

    call prof_rapstart('FILE_O_NetCDF', 2)

    call historywriteall( time_nowstep ) ![IN]

    call prof_rapend  ('FILE_O_NetCDF', 2)

    return
  end subroutine hist_write

end module scale_history