scale_file_history_cartesC.F90

Go to the documentation of this file.

!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------
#include "scalelib.h"
module scale_file_history_cartesc
  !-----------------------------------------------------------------------------
  !
  !++ Used modules
  !
  use scale_precision
  use scale_io
  use scale_prof
  use scale_atmos_grid_cartesc_index
  use scale_ocean_grid_cartesc_index
  use scale_land_grid_cartesc_index
  use scale_urban_grid_cartesc_index
  use scale_prc, only: &
     prc_abort
  !-----------------------------------------------------------------------------
  implicit none
  private
  !-----------------------------------------------------------------------------
  !
  !++ Public procedures
  !
  public :: file_history_cartesc_setup
  public :: file_history_cartesc_set_pres

  !-----------------------------------------------------------------------------
  !
  !++ Public parameters & variables
  !
  !-----------------------------------------------------------------------------
  !
  !++ Private procedures
  !
  private :: file_history_cartesc_set_dims
  private :: file_history_cartesc_set_axes
  private :: file_history_cartesc_set_axes_attributes

  !-----------------------------------------------------------------------------
  !
  !++ Private parameters & variables
  !
  integer,          parameter :: nzs = 3
  character(len=8), parameter :: zs(nzs) = (/ "model   ", &
                                              "z       ", &
                                              "pressure"  /)

  integer               :: file_history_cartesc_pres_nlayer = 0
  real(RP), allocatable :: file_history_cartesc_pres_val(:)

  integer  :: im,   jm,  km
  integer  :: ims,  ime
  integer  :: jms,  jme
  integer  :: imh,  jmh
  integer  :: imsh, jmsh

  integer  :: file_history_cartesc_startdate(6)
  real(DP) :: file_history_cartesc_startms

  logical  :: file_history_cartesc_boundary = .false.

  !-----------------------------------------------------------------------------
contains
  !-----------------------------------------------------------------------------
  subroutine file_history_cartesc_setup
    use scale_file_h, only: &
       file_hshort
    use scale_file_history, only: &
       file_history_setup, &
       file_history_set_nowdate, &
       file_history_truncate_1d, &
       file_history_truncate_2d, &
       file_history_truncate_3d
    use scale_prc, only: &
       prc_masterrank, &
       prc_myrank
    use scale_prc_cartesc, only: &
       prc_periodic_x, &
       prc_periodic_y, &
       prc_has_w,      &
       prc_has_s
    use scale_time, only: &
       time_nowdate,     &
       time_nowms,       &
       time_nowstep,     &
       time_dtsec,       &
       time_startdaysec
    use scale_calendar, only: &
       calendar_get_name
    use scale_interp_vert, only: &
       interp_vert_alloc_pres
    implicit none

    integer, parameter :: nlayer_max = 300
    real(RP)           :: FILE_HISTORY_CARTESC_PRES(nlayer_max)

    namelist / param_file_history_cartesc / &
       file_history_cartesc_pres_nlayer, &
       file_history_cartesc_pres,        &
       file_history_cartesc_boundary

    character(len=H_MID) :: FILE_HISTORY_CARTESC_H_TITLE = 'SCALE-RM FILE_HISTORY_CARTESC OUTPUT'
    character(len=H_MID) :: FILE_HISTORY_CARTESC_T_SINCE

    character(len=FILE_HSHORT) :: calendar
    real(DP) :: start_daysec
    integer  :: ierr
    integer  :: k
    !---------------------------------------------------------------------------

    log_newline
    log_info("FILE_HISTORY_CARTESC_setup",*) 'Setup'

    file_history_cartesc_pres(:) = 0.0_rp

    !--- read namelist
    rewind(io_fid_conf)
    read(io_fid_conf,nml=param_file_history_cartesc,iostat=ierr)
    if( ierr < 0 ) then !--- missing
       log_info("FILE_HISTORY_CARTESC_setup",*) 'Not found namelist. Default used.'
    elseif( ierr > 0 ) then !--- fatal error
       log_error("FILE_HISTORY_CARTESC_setup",*) 'Not appropriate names in namelist PARAM_FILE_HISTORY_CARTESC. Check!'
       call prc_abort
    endif
    log_nml(param_file_history_cartesc)


    ! check pressure coordinate
    if ( file_history_cartesc_pres_nlayer > 0 ) then
       if ( file_history_cartesc_pres_nlayer > nlayer_max ) then
          log_error("FILE_HISTORY_CARTESC_setup",'(a,i3)') 'FILE_HISTORY_CARTESC_PRES_nlayer must be <= ', nlayer_max
          call prc_abort
       end if
       allocate( file_history_cartesc_pres_val(file_history_cartesc_pres_nlayer) )

       do k = 1, file_history_cartesc_pres_nlayer
          if ( file_history_cartesc_pres(k) <= 0.0_rp ) then
             log_error("FILE_HISTORY_CARTESC_setup",'(a,i3,f7.1)') 'Invalid value found in pressure coordinate! (k,value)=', k, file_history_cartesc_pres(k)
             call prc_abort
          elseif ( file_history_cartesc_pres(k+1) >= file_history_cartesc_pres(k) ) then
             log_error("FILE_HISTORY_CARTESC_setup",'(a,i3,2f7.1)') 'The value of pressure coordinate must be descending order! ', &
                  '(k,value[k],value[k+1])=', k, file_history_cartesc_pres(k), file_history_cartesc_pres(k+1)
             call prc_abort
          endif
          file_history_cartesc_pres_val(k) = file_history_cartesc_pres(k) * 100.0_rp ! [hPa->Pa]
       enddo

       call interp_vert_alloc_pres( file_history_cartesc_pres_nlayer, ia, ja ) ! [IN]
    else
       log_info("FILE_HISTORY_CARTESC_setup",*) 'FILE_HISTORY_CARTESC_PRES_nlayer is not set.'
       log_info("FILE_HISTORY_CARTESC_setup",*) 'Output with pressure coordinate is disabled'
    endif



    file_history_cartesc_startdate(:) = time_nowdate
    file_history_cartesc_startms      = time_nowms

    start_daysec = time_startdaysec
    if ( time_nowdate(1) > 0 ) then
       write(file_history_cartesc_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
       file_history_cartesc_t_since = ''
    endif

    if ( file_history_cartesc_boundary ) then
       ims  = isb
       ime  = ieb
       jms  = jsb
       jme  = jeb

       imsh = ims
       jmsh = jms

       im   = imaxb
       jm   = jmaxb
       imh  = im
       jmh  = jm
    else
       ims  = is
       ime  = ie
       jms  = js
       jme  = je

       if ( prc_has_w .OR. prc_periodic_x ) then
          imsh = ims
       else
          imsh = ims - 1 ! including i = IS-1
       endif
       if ( prc_has_s .OR. prc_periodic_y ) then
          jmsh = jms
       else
          jmsh = jms - 1 ! include j = JS-1
       endif

       im   = ime - ims  + 1
       jm   = jme - jms  + 1
       imh  = ime - imsh + 1
       jmh  = jme - jmsh + 1
    endif

    ! get calendar name
    call calendar_get_name( calendar )

    call file_history_setup( file_history_cartesc_h_title,              & ! [IN]
                             h_source, h_institute,                     & ! [IN]
                             start_daysec, time_dtsec,                  & ! [IN]
                             time_since = file_history_cartesc_t_since, & ! [IN]
                             calendar = calendar,                       & ! [IN]
                             default_zcoord = 'model',                  & ! [IN]
                             myrank = prc_myrank                        ) ! [IN]

    call file_history_set_nowdate( time_nowdate, time_nowms, time_nowstep )

    call file_history_cartesc_set_dims

    call file_history_cartesc_set_axes

    call file_history_cartesc_set_axes_attributes

    file_history_truncate_1d => file_history_cartesc_truncate_1d
    file_history_truncate_2d => file_history_cartesc_truncate_2d
    file_history_truncate_3d => file_history_cartesc_truncate_3d

    return
  end subroutine file_history_cartesc_setup

  !-----------------------------------------------------------------------------
  !-----------------------------------------------------------------------------
  subroutine file_history_cartesc_set_pres( &
       PRES,    &
       PRESH,   &
       SFC_PRES )
    use scale_interp_vert, only: &
       interp_vert_setcoef_pres
    implicit none

    real(RP), intent(in) :: PRES    (:,:,:) ! pressure at the full level [Pa]
    real(RP), intent(in) :: PRESH   (:,:,:) ! pressure at the half level [Pa]
    real(RP), intent(in) :: SFC_PRES(  :,:) ! surface pressure           [Pa]
    !---------------------------------------------------------------------------

    if ( file_history_cartesc_pres_nlayer > 0 ) then
       call interp_vert_setcoef_pres( file_history_cartesc_pres_nlayer, & ! [IN]
                                      ka, ks, ke,                       & ! [IN]
                                      ia, isb, ieb,                     & ! [IN]
                                      ja, jsb, jeb,                     & ! [IN]
                                      pres(:,:,:),             & ! [IN]
                                      presh(:,:,:),             & ! [IN]
                                      sfc_pres(:,:)  ,             & ! [IN]
                                      file_history_cartesc_pres_val(:)  ) ! [IN]
    endif

    return
  end subroutine file_history_cartesc_set_pres

  ! private routines

  !-----------------------------------------------------------------------------
  !-----------------------------------------------------------------------------
  subroutine file_history_cartesc_set_dims
    use scale_prc, only: &
       prc_myrank
    use scale_prc_cartesc, only: &
       prc_2drank, &
       prc_has_w, &
       prc_has_s
    use scale_file_history, only: &
       file_history_set_dim
    use scale_mapprojection, only: &
       mapprojection_get_attributes
    implicit none

    character(len=H_SHORT) :: mapping

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

    integer :: start(3,3), count(3,3)
    integer :: xs, xc, ys, yc
    !---------------------------------------------------------------------------

    ! get start and count for x and y
    if ( file_history_cartesc_boundary ) then
       xs = isgb
       ys = jsgb
       xc = imaxb
       yc = jmaxb
    else
       ! for the case the shared-file contains no halos
       xs = prc_2drank(prc_myrank,1) * imax + 1 ! no IHALO
       xc = imax
       ys = prc_2drank(prc_myrank,2) * jmax + 1 ! no JHALO
       yc = jmax
    end if

    ! get mapping name
    call mapprojection_get_attributes( mapping )


    !  Vertical 1D
    start(1,1) = 1
    dims(1,1) = "z"
    count(1,1) = kmax
    call file_history_set_dim( "Z",   1, 1, dims(:,:), zs(:), start(:,:), count(:,:) ) ! [IN]
    dims(1,1) = "zh"
    count(1,1) = kmax + 1
    call file_history_set_dim( "ZH",  1, 1, dims(:,:), zs(:), start(:,:), count(:,:) ) ! [IN]

    dims(1,1) = "oz"
    count(1,1) = okmax
    call file_history_set_dim( "OZ",  1, 1, dims(:,:), zs(:), start(:,:), count(:,:) ) ! [IN]
    dims(1,1) = "ozh"
    count(1,1) = okmax + 1
    call file_history_set_dim( "OZH", 1, 1, dims(:,:), zs(:), start(:,:), count(:,:) ) ! [IN]

    dims(1,1) = "lz"
    count(1,1) = lkmax
    call file_history_set_dim( "LZ",  1, 1, dims(:,:), zs(:), start(:,:), count(:,:) ) ! [IN]
    dims(1,1) = "lzh"
    count(1,1) = lkmax + 1
    call file_history_set_dim( "LZH", 1, 1, dims(:,:), zs(:), start(:,:), count(:,:) ) ! [IN]

    dims(1,1) = "uz"
    count(1,1) = ukmax
    call file_history_set_dim( "UZ",  1, 1, dims(:,:), zs(:), start(:,:), count(:,:) ) ! [IN]
    dims(1,1) = "uzh"
    count(1,1) = ukmax + 1
    call file_history_set_dim( "UZH", 1, 1, dims(:,:), zs(:), start(:,:), count(:,:) ) ! [IN]

    ! X, Y
    start(1,:) = xs
    start(2,:) = ys
    dims(1,:) = 'lon'
    dims(2,:) = 'lat'
    count(1,:) = xc
    count(2,:) = yc
    call file_history_set_dim( "XY", 2, 1, dims(:,:), zs(:), start(:,:), count(:,:), mapping=mapping, area="cell_area", location="face" ) ! [IN]

    start(3,:) = 1
    dims(3,:) = (/ "height    ", "z         ", "pressure  " /)
    count(3,:) = (/ kmax,   kmax,   file_history_cartesc_pres_nlayer /)
    call file_history_set_dim( "ZXY",  3, nzs, dims(:,:), zs(:), start(:,:), count(:,:), mapping=mapping, & ! [IN]
                               area="cell_area", area_x="cell_area_xyz_x", area_y="cell_area_xyz_y", volume="cell_volume", location="face" ) ! [IN]
    dims(3,:) = (/ "height_xyw", "zh        ", "pressure  " /)
    count(3,:) = (/ kmax+1, kmax+1, file_history_cartesc_pres_nlayer /)
    call file_history_set_dim( "ZHXY", 3, nzs, dims(:,:), zs(:), start(:,:), count(:,:), mapping=mapping, area="cell_area", volume="cell_volume_xyw", location="face" )

    if ( okmax > 0 ) then
       dims(3,1) = "oz"
       count(3,1) = okmax
       call file_history_set_dim( "OXY",  3, 1, dims(:,:), zs(:), start(:,:), count(:,:), mapping=mapping, area="cell_area", volume="cell_volume_xyo", location="face", grid="ocean" ) ! [IN]
       dims(3,1) = "ozh"
       count(3,1) = okmax + 1
       call file_history_set_dim( "OHXY", 3, 1, dims(:,:), zs(:), start(:,:), count(:,:), mapping=mapping, area="cell_area", location="face", grid="ocean" ) ! [IN]
    endif

    if ( lkmax > 0 ) then
       dims(3,1) = "lz"
       count(3,1) = lkmax
       call file_history_set_dim( "LXY",  3, 1, dims(:,:), zs(:), start(:,:), count(:,:), mapping=mapping, area="cell_area", volume="cell_volume_xyl", location="face", grid="land" ) ! [IN]
       dims(3,1) = "lzh"
       count(3,1) = lkmax + 1
       call file_history_set_dim( "LHXY", 3, 1, dims(:,:), zs(:), start(:,:), count(:,:), mapping=mapping, area="cell_area", location="face", grid="land" ) ! [IN]
    endif

    if ( ukmax > 0 ) then
       dims(3,1) = "uz"
       count(3,1) = ukmax
       call file_history_set_dim( "UXY",  3, 1, dims(:,:), zs(:), start(:,:), count(:,:), mapping=mapping, area="cell_area", volume="cell_volume_xyu", grid="urban" ) ! [IN]
       dims(3,1) = "uzh"
       count(3,1) = ukmax + 1
       call file_history_set_dim( "UHXY", 3, 1, dims(:,:), zs(:), start(:,:), count(:,:), mapping=mapping, area="cell_area", location="face", grid="urban" ) ! [IN]
    endif

    ! XH, Y
    dims(1,:) = 'lon_uy'
    dims(2,:) = 'lat_uy'
    if ( prc_has_w ) then
       start(1,:) = xs+1
       count(1,:) = xc
    else
       start(1,:) = xs
       count(1,:) = xc+1
    endif
    call file_history_set_dim( "XHY", 2, 1, dims(:,:), zs(:), start(:,:), count(:,:), mapping=mapping, area="cell_area_uy", location="edge1" ) ! [IN]

    dims(3,:) = (/ "height_uyz", "z         ", "pressure  " /)
    count(3,:) = (/ kmax,   kmax,   file_history_cartesc_pres_nlayer /)
    call file_history_set_dim( "ZXHY",  3, nzs, dims(:,:), zs(:), start(:,:), count(:,:), mapping=mapping, area="cell_area_uy", volume="cell_volume_uyz", location="edge1" ) ! [IN]
    dims(3,:) = (/ "height_uyw", "zh        ", "pressure  " /)
    count(3,:) = (/ kmax+1, kmax+1, file_history_cartesc_pres_nlayer /)
    call file_history_set_dim( "ZHXHY", 3, nzs, dims(:,:), zs(:), start(:,:), count(:,:), mapping=mapping, area="cell_area_uy", location="edge1" ) ! [IN]

    ! X, YH
    dims(1,:) = 'lon_xv'
    dims(2,:) = 'lat_xv'
    start(1,:) = xs
    count(1,:) = xc
    if ( prc_has_s ) then
       start(2,:) = ys+1
       count(2,:) = yc
    else
       start(2,:) = ys
       count(2,:) = yc+1
    endif
    call file_history_set_dim( "XYH", 2, 1, dims(:,:), zs(:), start(:,:), count(:,:), mapping=mapping, area="cell_area_xv", location="edge2" ) ! [IN]

    dims(3,:) = (/ "height_xvz", "z         ", "pressure  " /)
    count(3,:) = (/ kmax,   kmax,   file_history_cartesc_pres_nlayer /)
    call file_history_set_dim( "ZXYH",  3, nzs, dims(:,:), zs(:), start(:,:), count(:,:), mapping=mapping, area="cell_area_xv", volume="cell_volume_xvz", location="edge2" ) ! [IN]
    dims(3,:) = (/ "height_xvw", "zh        ", "pressure  " /)
    count(3,:) = (/ kmax+1, kmax+1, file_history_cartesc_pres_nlayer /)
    call file_history_set_dim( "ZHXYH", 3, nzs, dims(:,:), zs(:), start(:,:), count(:,:), mapping=mapping, area="cell_area_xv", location="edge2" ) ! [IN]

    ! XH, YH
    dims(1,:) = 'lon_uv'
    dims(2,:) = 'lat_uv'
    if ( prc_has_w ) then
       start(1,:) = xs+1
       count(1,:) = xc
    else
       start(1,:) = xs
       count(1,:) = xc+1
    endif
    call file_history_set_dim( "XHYH", 2, 1, dims(:,:), zs(:), start(:,:), count(:,:), mapping=mapping, location="face" ) ! [IN]

    dims(3,:) = (/ "height_uvz", "z         ", "pressure  " /)
    count(3,:) = (/ kmax,   kmax,   file_history_cartesc_pres_nlayer /)
    call file_history_set_dim( "ZXHYH",  3, nzs, dims(:,:), zs(:), start(:,:), count(:,:), mapping=mapping, & ![IN]
                               area="cell_area_uv", area_x="cell_area_uvz_x", area_y="cell_area_uvz_y", location="node" ) ! [IN]
    dims(3,:) = (/ "height_uvw", "zh        ", "pressure  " /)
    count(3,:) = (/ kmax+1, kmax+1, file_history_cartesc_pres_nlayer /)
    call file_history_set_dim( "ZHXHYH", 3, nzs, dims(:,:), zs(:), start(:,:), count(:,:), mapping=mapping, location="node" ) ! [IN]

    return
  end subroutine file_history_cartesc_set_dims

  !-----------------------------------------------------------------------------
  subroutine file_history_cartesc_truncate_1d( &
       src, &
       dim_type, zcoord, fill_halo, &
       dst )
    implicit none

    real(RP),         intent(in) :: src(:)
    character(len=*), intent(in) :: dim_type
    character(len=*), intent(in) :: zcoord
    logical,          intent(in) :: fill_halo ! ignored

    real(DP), intent(out) :: dst(:)

    integer  :: ksize
    integer  :: kstart
    integer  :: k
    !---------------------------------------------------------------------------

    ! select dimension
    select case ( dim_type )
    case ('Z')
        ksize  = kmax
        kstart = ks
    case ('ZH')
       ksize  = kmax+1
       kstart = ks-1
    case ('OZ')
       ksize  = okmax
       kstart = oks
    case ('OZH')
       ksize  = okmax+1
       kstart = oks-1
    case ('LZ')
       ksize  = lkmax
       kstart = lks
    case ('LZH')
       ksize  = lkmax+1
       kstart = lks-1
    case ('UZ')
       ksize  = ukmax
       kstart = uks
    case ('UZH')
       ksize  = ukmax+1
       kstart = uks-1
    case default
       log_error("FILE_HISTORY_CARTESC_truncate_1D",*) 'dim_type is invalid: ', trim(dim_type)
       call prc_abort
    end select

    do k = 1, ksize
       dst(k) = src(kstart+k-1)
    enddo

    return
  end subroutine file_history_cartesc_truncate_1d

  !-----------------------------------------------------------------------------
  subroutine file_history_cartesc_truncate_2d( &
       src, &
       dim_type, zcoord, fill_halo, &
       dst )
    use scale_file_h, only: &
       rmiss => file_rmiss
    implicit none

    real(RP),         intent(in) :: src(:,:)
    character(len=*), intent(in) :: dim_type
    character(len=*), intent(in) :: zcoord    ! ignored
    logical,          intent(in) :: fill_halo

    real(DP), intent(out) :: dst(:)

    integer :: isize, jsize
    integer :: istart, jstart
    integer :: i, j
    !---------------------------------------------------------------------------

    ! select dimension
    select case( dim_type )
    case ( 'XY', 'XYH' )
       isize  = im
       istart = ims
    case ( 'XHY', 'XHYH' )
       isize  = imh
       istart = imsh
    case default
       log_error("FILE_HISTORY_CARTESC_truncate_2D",*) 'dim_type is invalid: ', trim(dim_type)
       call prc_abort
    end select

    select case ( dim_type )
    case ( 'XY', 'XHY' )
       jsize  = jm
       jstart = jms
    case ( 'XYH', 'XHYH' )
       jsize  = jmh
       jstart = jmsh
    case default
       log_error("FILE_HISTORY_CARTESC_truncate_2D",*) 'dim_type is invalid: ', trim(dim_type)
       call prc_abort
    end select

    !$omp parallel do
    do j = 1, jsize
    do i = 1, isize
       dst((j-1)*isize+i) = src(istart+i-1,jstart+j-1)
    enddo
    enddo

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

    end if

    return
  end subroutine file_history_cartesc_truncate_2d

  !-----------------------------------------------------------------------------
  subroutine file_history_cartesc_truncate_3d( &
       src, &
       dim_type, zcoord, fill_halo, &
       dst )
    use scale_file_h, only: &
       rmiss => file_rmiss
    use scale_atmos_grid_cartesc, only: &
       atmos_grid_cartesc_cz, &
       atmos_grid_cartesc_fz
    use scale_atmos_grid_cartesc_real, only: &
       atmos_grid_cartesc_real_cz, &
       atmos_grid_cartesc_real_fz
    use scale_interp_vert, only: &
       interp_vert_xi2z,   &
       interp_vert_xi2p,   &
       interp_vert_xih2zh, &
       interp_vert_xih2p,  &
       interp_available
    implicit none

    real(RP),         intent(in)  :: src(:,:,:)
    character(len=*), intent(in)  :: dim_type
    character(len=*), intent(in)  :: zcoord
    logical,          intent(in)  :: fill_halo
    real(DP),         intent(out) :: dst(:)

    real(RP) :: src_Z(ka,ia,ja)
    real(RP) :: src_P(file_history_cartesc_pres_nlayer,ia,ja)

    integer  :: isize,  jsize,  ksize
    integer  :: istart, jstart, kstart
    integer  :: i, j, k
    !---------------------------------------------------------------------------

    ! select dimension
    if ( index( dim_type, 'XH' ) == 0 ) then
       isize  = im
       istart = ims
    else
       isize  = imh
       istart = imsh
    end if

    if ( index( dim_type, 'YH' ) == 0 ) then
       jsize  = jm
       jstart = jms
    else
       jsize  = jmh
       jstart = jmsh
    end if

    select case( dim_type(1:1) )
    case ( 'Z' )
       ksize  = kmax
       kstart = ks
    case('O')
       ksize  = okmax
       kstart = oks
    case('L')
       ksize  = lkmax
       kstart = lks
    case('U')
       ksize  = ukmax
       kstart = uks
    case default
       log_error("FILE_HISTORY_CARTESC_truncate_3D",*) 'dim_type is invalid: ', trim(dim_type)
       call prc_abort
    end select
    if ( dim_type(2:2) == 'H' ) then
       ksize = ksize + 1
       kstart = kstart - 1
    end if


    if ( ksize == kmax .and. zcoord == "z" .and. interp_available ) then ! z*->z interpolation (full level)

       call prof_rapstart('FILE_O_interp', 2)
       call interp_vert_xi2z( ka, ks, ke, ia, isb, ieb, ja, jsb, jeb, &
                              atmos_grid_cartesc_cz(:),          & ! [IN]
                              atmos_grid_cartesc_real_cz(:,:,:), & ! [IN]
                              src(:,:,:), & ! [IN]
                              src_z(:,:,:)  ) ! [OUT]
       call prof_rapend  ('FILE_O_interp', 2)

       !$omp parallel do
       do k = 1, ksize
       do j = 1, jsize
       do i = 1, isize
          dst((k-1)*jsize*isize+(j-1)*isize+i) = src_z(kstart+k-1,istart+i-1,jstart+j-1)
       enddo
       enddo
       enddo

    else if( ksize == kmax+1 .and. zcoord == "z" .and. interp_available ) then ! z*->z interpolation (half level)


       call prof_rapstart('FILE_O_interp', 2)
       call interp_vert_xih2zh( ka, ks, ke, ia, isb, ieb, ja, jsb, jeb, &
                                atmos_grid_cartesc_fz(:),          & ! [IN]
                                atmos_grid_cartesc_real_fz(:,:,:), & ! [IN]
                                src(:,:,:), & ! [IN]
                                src_z(:,:,:)  ) ! [OUT]
       call prof_rapend  ('FILE_O_interp', 2)

       !$omp parallel do
       do k = 1, ksize
       do j = 1, jsize
       do i = 1, isize
          dst((k-1)*jsize*isize+(j-1)*isize+i) = src_z(kstart+k-1,istart+i-1,jstart+j-1)
       enddo
       enddo
       enddo

    elseif( ksize == kmax .and. zcoord == "pressure" ) then ! z*->p interpolation (full level)
       ksize = file_history_cartesc_pres_nlayer
       if ( ksize == 0 ) then
          log_error("FILE_HISTORY_CARTESC_truncate_3D",*) 'FILE_HISTORY_CARTESC_PRES_nlayer must be set to output variable with the pressure coordinate'
          call prc_abort
       end if

       call prof_rapstart('FILE_O_interp', 2)
       call interp_vert_xi2p( file_history_cartesc_pres_nlayer, & ! [IN]
                              ka,                               & ! [IN]
                              ia, isb, ieb,                     & ! [IN]
                              ja, jsb, jeb,                     & ! [IN]
                              src(:,:,:),                     & ! [IN]
                              src_p(:,:,:)                      ) ! [OUT]
       call prof_rapend  ('FILE_O_interp', 2)

       !$omp parallel do
       do k = 1, ksize
       do j = 1, jsize
       do i = 1, isize
          dst((k-1)*jsize*isize+(j-1)*isize+i) = src_p(k,istart+i-1,jstart+j-1)
       enddo
       enddo
       enddo

    elseif( ksize == kmax+1 .and. zcoord == "pressure" ) then ! z*->p interpolation (half level)
       ksize = file_history_cartesc_pres_nlayer
       if ( ksize == 0 ) then
          log_error("FILE_HISTORY_CARTESC_truncate_3D",*) 'FILE_HISTORY_CARTESC_PRES_nlayer must be set to output variable with the pressure coordinate'
          call prc_abort
       end if

       call prof_rapstart('FILE_O_interp', 2)
       call interp_vert_xih2p( file_history_cartesc_pres_nlayer, & ! [IN]
                               ka,                               & ! [IN]
                               ia, isb, ieb,                     & ! [IN]
                               ja, jsb, jeb,                     & ! [IN]
                               src(:,:,:),                     & ! [IN]
                               src_p(:,:,:)                      ) ! [OUT]
       call prof_rapend  ('FILE_O_interp', 2)

       do k = 1, ksize
       do j = 1, jsize
       do i = 1, isize
          dst((k-1)*jsize*isize+(j-1)*isize+i) = src_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
          dst((k-1)*jsize*isize+(j-1)*isize+i) = src(kstart+k-1,istart+i-1,jstart+j-1)
       enddo
       enddo
       enddo

    endif

    if ( fill_halo ) then
       ! W halo
       do k = 1, ksize
       do j = 1, jsize
       do i = 1, is-istart
          dst((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
          dst((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
          dst((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
          dst((k-1)*jsize*isize+(j-1)*isize+i) = rmiss
       enddo
       enddo
       enddo
    endif

    return
  end subroutine file_history_cartesc_truncate_3d

  !-----------------------------------------------------------------------------
  !  only register the axis and coordinate variables into internal buffers
  !  The actual write happens later when calling FILE_HISTORY_CARTESC_write
  subroutine file_history_cartesc_set_axes
    use scale_file_history, only: &
       file_history_aggregate, &
       file_history_set_axis, &
       file_history_set_associatedcoordinate
    use scale_const, only: &
       undef => const_undef, &
       d2r => const_d2r
    use scale_prc, only: &
       prc_myrank
    use scale_prc_cartesc, only: &
       prc_2drank,     &
       prc_num_x,      &
       prc_num_y,      &
       prc_periodic_x, &
       prc_periodic_y, &
       prc_has_w,      &
       prc_has_s
    use scale_atmos_grid_cartesc, only: &
       atmos_grid_cartesc_cz,    &
       atmos_grid_cartesc_cx,    &
       atmos_grid_cartesc_cy,    &
       atmos_grid_cartesc_fz,    &
       atmos_grid_cartesc_fx,    &
       atmos_grid_cartesc_fy,    &
       atmos_grid_cartesc_cdz,   &
       atmos_grid_cartesc_cdx,   &
       atmos_grid_cartesc_cdy,   &
       atmos_grid_cartesc_fdz,   &
       atmos_grid_cartesc_fdx,   &
       atmos_grid_cartesc_fdy,   &
       atmos_grid_cartesc_cbfz,  &
       atmos_grid_cartesc_cbfx,  &
       atmos_grid_cartesc_cbfy,  &
       atmos_grid_cartesc_fbfz,  &
       atmos_grid_cartesc_fbfx,  &
       atmos_grid_cartesc_fbfy,  &
       atmos_grid_cartesc_cxg,   &
       atmos_grid_cartesc_cyg,   &
       atmos_grid_cartesc_fxg,   &
       atmos_grid_cartesc_fyg,   &
       atmos_grid_cartesc_cdxg,  &
       atmos_grid_cartesc_cdyg,  &
       atmos_grid_cartesc_fdxg,  &
       atmos_grid_cartesc_fdyg,  &
       atmos_grid_cartesc_cbfxg, &
       atmos_grid_cartesc_cbfyg, &
       atmos_grid_cartesc_fbfxg, &
       atmos_grid_cartesc_fbfyg
    use scale_ocean_grid_cartesc, only: &
       ocean_grid_cartesc_cz, &
       ocean_grid_cartesc_fz, &
       ocean_grid_cartesc_cdz
    use scale_land_grid_cartesc, only: &
       land_grid_cartesc_cz, &
       land_grid_cartesc_fz, &
       land_grid_cartesc_cdz
    use scale_urban_grid_cartesc, only: &
       urban_grid_cartesc_cz, &
       urban_grid_cartesc_fz, &
       urban_grid_cartesc_cdz
    use scale_atmos_grid_cartesc_real, only: &
       atmos_grid_cartesc_real_cz,    &
       atmos_grid_cartesc_real_fz,    &
       atmos_grid_cartesc_real_lon,   &
       atmos_grid_cartesc_real_lonuy, &
       atmos_grid_cartesc_real_lonxv, &
       atmos_grid_cartesc_real_lonuv, &
       atmos_grid_cartesc_real_lat,   &
       atmos_grid_cartesc_real_latuy, &
       atmos_grid_cartesc_real_latxv, &
       atmos_grid_cartesc_real_latuv, &
       area   => atmos_grid_cartesc_real_area,   &
       areauy => atmos_grid_cartesc_real_areauy, &
       areaxv => atmos_grid_cartesc_real_areaxv, &
       areazuy_x => atmos_grid_cartesc_real_areazuy_x, &
       areazxv_y => atmos_grid_cartesc_real_areazxv_y, &
       areawuy_x => atmos_grid_cartesc_real_areawuy_x, &
       areawxv_y => atmos_grid_cartesc_real_areawxv_y, &
       areazxy_x => atmos_grid_cartesc_real_areazxy_x, &
       areazuv_y => atmos_grid_cartesc_real_areazuv_y, &
       areazuv_x => atmos_grid_cartesc_real_areazuv_x, &
       areazxy_y => atmos_grid_cartesc_real_areazxy_y, &
       vol    => atmos_grid_cartesc_real_vol,    &
       volwxy => atmos_grid_cartesc_real_volwxy, &
       volzuy => atmos_grid_cartesc_real_volzuy, &
       volzxv => atmos_grid_cartesc_real_volzxv
    use scale_ocean_grid_cartesc_real, only: &
       volo => ocean_grid_cartesc_real_vol
    use scale_land_grid_cartesc_real, only: &
       voll => land_grid_cartesc_real_vol
    use scale_urban_grid_cartesc_real, only: &
       volu => urban_grid_cartesc_real_vol
    use scale_topography, only: &
       topo_zsfc
    use scale_landuse, only: &
       landuse_frac_land
    implicit none

    real(RP)         :: AXIS (imh,jmh,0:kmax)
    real(RP)         :: AXISO(im, jm, okmax)
    real(RP)         :: AXISL(im, jm, lkmax)
    real(RP)         :: AXISU(im, jm, ukmax)
    character(len=2) :: AXIS_name(3)

    integer :: rankidx(2)
    integer :: start(3,4)
    integer :: startX, startY, startZ
    integer :: startXH, startYH
    integer :: XAG, YAG
    integer :: XAGH, YAGH

    real(RP) :: z_bnds(2,ka), zh_bnds(2,0:ka)
    real(RP) :: oz_bnds(2,oka), ozh_bnds(2,0:oka)
    real(RP) :: lz_bnds(2,lka), lzh_bnds(2,0:lka)
    real(RP) :: uz_bnds(2,uka), uzh_bnds(2,0:uka)
    real(RP) :: x_bnds(2,ia), xh_bnds(2,0:ia)
    real(RP) :: y_bnds(2,ja), yh_bnds(2,0:ja)

    real(RP) :: FDXG(0:iag), FDYG(0:jag)
    real(RP) :: FDX(0:ia), FDY(0:ja)

    integer :: k, i, j
    !---------------------------------------------------------------------------

    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 ( file_history_aggregate ) then

       startz = 1

       if ( file_history_cartesc_boundary ) then
          startx  = isgb ! global subarray starting index
          starty  = jsgb ! global subarray starting index
          startxh = startx
          startyh = starty
          xag     = iagb
          yag     = jagb
          xagh    = xag
          yagh    = yag
       else
          startx  = prc_2drank(prc_myrank,1) * imax + 1
          starty  = prc_2drank(prc_myrank,2) * jmax + 1
          startxh = startx
          startyh = starty
          xag     = imaxg
          yag     = jmaxg
          xagh    = xag
          yagh    = yag

          if ( .NOT. prc_periodic_x ) then
             xagh = xagh + 1
             if( prc_has_w ) startxh = startxh + 1
          endif

          if ( .NOT. prc_periodic_y ) then
             yagh = yagh + 1
             if( prc_has_s ) startyh = startyh + 1
          endif
       endif

       ! for shared-file parallel I/O, only a part of rank writes variables
       if ( prc_myrank > 0 ) then ! only rank 0 writes Z axes
          startz = -1
       endif
       if ( rankidx(2) > 0 ) then ! only south-most processes write
          startx  = -1
          startxh = -1
       endif
       if ( rankidx(1) > 0 ) then ! only west-most processes write
          starty  = -1
          startyh = -1
       endif
    else
       startz  = 1
       startx  = 1
       starty  = 1
       startxh = startx
       startyh = starty
       xag     = im
       yag     = jm
       xagh    = imh
       yagh    = jmh
    endif

    ! bounds
    do k = ks, ke
       z_bnds(1,k) = atmos_grid_cartesc_fz(k-1)
       z_bnds(2,k) = atmos_grid_cartesc_fz(k  )
    end do
    do k = ks-1, ke
       zh_bnds(1,k) = atmos_grid_cartesc_cz(k  )
       zh_bnds(2,k) = atmos_grid_cartesc_cz(k+1)
    end do

    if ( oka > 0 ) then
       do k = oks, oke
          oz_bnds(1,k) = ocean_grid_cartesc_fz(k-1)
          oz_bnds(2,k) = ocean_grid_cartesc_fz(k  )
       end do
       ozh_bnds(1,oks-1) = ocean_grid_cartesc_fz(oks-1)
       do k = oks-1, oke-1
          ozh_bnds(2,k  ) = ocean_grid_cartesc_cz(k+1)
          ozh_bnds(1,k+1) = ocean_grid_cartesc_cz(k+1)
       end do
       ozh_bnds(2,oke) = ocean_grid_cartesc_fz(oke)
    end if

    if ( lka > 0 ) then
       do k = lks, lke
          lz_bnds(1,k) = land_grid_cartesc_fz(k-1)
          lz_bnds(2,k) = land_grid_cartesc_fz(k  )
       end do
       lzh_bnds(1,lks-1) = land_grid_cartesc_fz(lks-1)
       do k = lks-1, lke-1
          lzh_bnds(2,k  ) = land_grid_cartesc_cz(k+1)
          lzh_bnds(1,k+1) = land_grid_cartesc_cz(k+1)
       end do
       lzh_bnds(2,lke) = land_grid_cartesc_fz(lke)
    end if

    if ( uka > 0 ) then
       do k = uks, uke
          uz_bnds(1,k) = urban_grid_cartesc_fz(k-1)
          uz_bnds(2,k) = urban_grid_cartesc_fz(k  )
       end do
       uzh_bnds(1,uks-1) = urban_grid_cartesc_fz(uks-1)
       do k = uks-1, uke-1
          uzh_bnds(2,k  ) = urban_grid_cartesc_cz(k+1)
          uzh_bnds(1,k+1) = urban_grid_cartesc_cz(k+1)
       end do
       uzh_bnds(2,uke) = urban_grid_cartesc_fz(uke)
    end if

    do i = ims, ime
       x_bnds(1,i) = atmos_grid_cartesc_fx(i-1)
       x_bnds(2,i) = atmos_grid_cartesc_fx(i  )
    end do
    if ( imsh == 0 ) then
       xh_bnds(1,0) = atmos_grid_cartesc_fx(0)
    else
       xh_bnds(1,imsh) = atmos_grid_cartesc_cx(imsh)
    end if
    do i = imsh, ime-1
       xh_bnds(2,i  ) = atmos_grid_cartesc_cx(i+1)
       xh_bnds(1,i+1) = atmos_grid_cartesc_cx(i+1)
    end do
    if ( ime == ia ) then
       xh_bnds(2,ime) = atmos_grid_cartesc_fx(ia)
    else
       xh_bnds(2,ime) = atmos_grid_cartesc_cx(ime+1)
    end if

    do j = jms, jme
       y_bnds(1,j) = atmos_grid_cartesc_fy(j-1)
       y_bnds(2,j) = atmos_grid_cartesc_fy(j  )
    end do
    if ( jmsh == 0 ) then
       yh_bnds(1,jmsh) = atmos_grid_cartesc_fy(jmsh)
    else
       yh_bnds(1,jmsh) = atmos_grid_cartesc_cy(jmsh)
    end if
    do j = jmsh, jme-1
       yh_bnds(2,j  ) = atmos_grid_cartesc_cy(j+1)
       yh_bnds(1,j+1) = atmos_grid_cartesc_cy(j+1)
    end do
    if ( jme == ja ) then
       yh_bnds(2,jme) = atmos_grid_cartesc_fy(jme)
    else
       yh_bnds(2,jme) = atmos_grid_cartesc_cy(jme+1)
    end if


    fdxg(1:iag-1) = atmos_grid_cartesc_fdxg(:)
    fdxg(0  ) = undef
    fdxg(iag) = undef
    fdyg(1:jag-1) = atmos_grid_cartesc_fdyg(:)
    fdyg(0  ) = undef
    fdyg(jag) = undef

    fdx(1:ia-1) = atmos_grid_cartesc_fdx(:)
    fdx(0 ) = fdxg(is_ing-ihalo-1)
    fdx(ia) = fdxg(ie_ing+ihalo  )
    fdy(1:ja-1) = atmos_grid_cartesc_fdy(:)
    fdy(0 ) = fdyg(js_ing-jhalo-1)
    fdy(ja) = fdyg(je_ing+jhalo  )


    ! 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 file_history_set_axis( 'z',   'Z',               'm', 'z',   atmos_grid_cartesc_cz(ks  :ke), &
                                bounds=z_bnds(:,ks  :ke), gsize=kmax   , start=startz                 )
    call file_history_set_axis( 'zh',  'Z (half level)',  'm', 'zh',  atmos_grid_cartesc_fz(ks-1:ke), &
                                bounds=zh_bnds(:,ks-1:ke), gsize=kmax+1 , start=startz                 )

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

    if ( oka > 0 ) then
       call file_history_set_axis( 'oz',  'OZ',              'm', 'oz',  ocean_grid_cartesc_cz(oks  :oke), &
                                   bounds=oz_bnds(:,oks  :oke), gsize=okmax  , start=startz, down=.true.  )
       call file_history_set_axis( 'ozh', 'OZ (half level)', 'm', 'ozh', ocean_grid_cartesc_fz(oks-1:oke), &
                                   bounds=ozh_bnds(:,oks-1:oke), gsize=okmax+1, start=startz, down=.true.  )
    endif

    if ( lka > 0 ) then
       call file_history_set_axis( 'lz',  'LZ',              'm', 'lz',  land_grid_cartesc_cz(lks  :lke),  &
                                   bounds=lz_bnds(:,lks  :lke), gsize=lkmax  , start=startz, down=.true.  )
       call file_history_set_axis( 'lzh', 'LZ (half level)', 'm', 'lzh', land_grid_cartesc_fz(lks-1:lke),  &
                                   bounds=lzh_bnds(:,lks-1:lke), gsize=lkmax+1, start=startz, down=.true.  )
    endif

    if ( uka > 0 ) then
       call file_history_set_axis( 'uz',  'UZ',              'm', 'uz',  urban_grid_cartesc_cz(uks  :uke), &
                                   bounds=uz_bnds(:,uks  :uke), gsize=ukmax  , start=startz, down=.true.  )
       call file_history_set_axis( 'uzh', 'UZ (half level)', 'm', 'uzh', urban_grid_cartesc_fz(uks-1:uke), &
                                   bounds=uzh_bnds(:,uks-1:uke), gsize=ukmax+1, start=startz, down=.true.  )
    endif

    call file_history_set_axis( 'x',   'X',               'm', 'x',   atmos_grid_cartesc_cx(ims :ime), &
                                bounds=x_bnds(:,ims :ime), gsize=xag , start=startx                    )
    call file_history_set_axis( 'xh',  'X (half level)',  'm', 'xh',  atmos_grid_cartesc_fx(imsh:ime), &
                                bounds=xh_bnds(:,imsh:ime), gsize=xagh, start=startxh                   )

    call file_history_set_axis( 'y',   'Y',               'm', 'y',   atmos_grid_cartesc_cy(jms :jme), &
                                bounds=y_bnds(:,jms :jme), gsize=yag , start=starty                    )
    call file_history_set_axis( 'yh',  'Y (half level)',  'm', 'yh',  atmos_grid_cartesc_fy(jmsh:jme), &
                                bounds=yh_bnds(:,jmsh:jme), gsize=yagh, start=startyh                   )

    ! axes below always include halos when written to file regardless of PRC_PERIODIC_X/PRC_PERIODIC_Y
    call file_history_set_axis( 'CZ',   'Atmos Grid Center Position Z', 'm', 'CZ',  atmos_grid_cartesc_cz,   gsize=ka,      start=startz )
    call file_history_set_axis( 'FZ',   'Atmos Grid Face Position Z',   'm', 'FZ',  atmos_grid_cartesc_fz,   gsize=ka+1,    start=startz )
    call file_history_set_axis( 'CDZ',  'Grid Cell length Z',           'm', 'CZ',  atmos_grid_cartesc_cdz,  gsize=ka,      start=startz )
    call file_history_set_axis( 'FDZ',  'Grid distance Z',              'm', 'FDZ', atmos_grid_cartesc_fdz,  gsize=ka-1,    start=startz )
    call file_history_set_axis( 'CBFZ', 'Boundary factor Center Z',     '1', 'CZ',  atmos_grid_cartesc_cbfz, gsize=ka,      start=startz )
    call file_history_set_axis( 'FBFZ', 'Boundary factor Face Z',       '1', 'FZ',  atmos_grid_cartesc_fbfz, gsize=ka+1,    start=startz )

    if ( okmax > 0 ) then
       call file_history_set_axis( 'OCZ',  'Ocean Grid Center Position Z', 'm', 'OCZ', ocean_grid_cartesc_cz,  gsize=okmax,   start=startz, down=.true. )
       call file_history_set_axis( 'OFZ',  'Ocean Grid Face Position Z',   'm', 'OFZ', ocean_grid_cartesc_fz,  gsize=okmax+1, start=startz, down=.true. )
       call file_history_set_axis( 'OCDZ', 'Ocean Grid Cell length Z',     'm', 'OCZ', ocean_grid_cartesc_cdz, gsize=okmax,   start=startz              )
    endif

    if ( lkmax > 0 ) then
       call file_history_set_axis( 'LCZ',  'Land Grid Center Position Z',  'm', 'LCZ', land_grid_cartesc_cz,  gsize=lkmax,   start=startz, down=.true. )
       call file_history_set_axis( 'LFZ',  'Land Grid Face Position Z',    'm', 'LFZ', land_grid_cartesc_fz,  gsize=lkmax+1, start=startz, down=.true. )
       call file_history_set_axis( 'LCDZ', 'Land Grid Cell length Z',      'm', 'LCZ', land_grid_cartesc_cdz, gsize=lkmax,   start=startz              )
    endif

    if ( ukmax > 0 ) then
       call file_history_set_axis( 'UCZ',  'Urban Grid Center Position Z', 'm', 'UCZ', urban_grid_cartesc_cz,  gsize=ukmax,   start=startz, down=.true. )
       call file_history_set_axis( 'UFZ',  'Urban Grid Face Position Z',   'm', 'UFZ', urban_grid_cartesc_fz,  gsize=ukmax+1, start=startz, down=.true. )
       call file_history_set_axis( 'UCDZ', 'Urban Grid Cell length Z',     'm', 'UCZ', urban_grid_cartesc_cdz, gsize=ukmax,   start=startz              )
    endif

    if ( file_history_aggregate ) then
       call file_history_set_axis( 'CX',   'Atmos Grid Center Position X', 'm', 'CX', atmos_grid_cartesc_cxg,   gsize=iag,   start=startz )
       call file_history_set_axis( 'CY',   'Atmos Grid Center Position Y', 'm', 'CY', atmos_grid_cartesc_cyg,   gsize=jag,   start=startz )
       call file_history_set_axis( 'FX',   'Atmos Grid Face Position X',   'm', 'FX', atmos_grid_cartesc_fxg,   gsize=iag+1, start=startz )
       call file_history_set_axis( 'FY',   'Atmos Grid Face Position Y',   'm', 'FY', atmos_grid_cartesc_fyg,   gsize=jag+1, start=startz )
       call file_history_set_axis( 'CDX',  'Grid Cell length X',           'm', 'CX', atmos_grid_cartesc_cdxg,  gsize=iag,   start=startz )
       call file_history_set_axis( 'CDY',  'Grid Cell length Y',           'm', 'CY', atmos_grid_cartesc_cdyg,  gsize=jag,   start=startz )
       call file_history_set_axis( 'FDX',  'Grid distance X',              'm', 'FX',                    fdxg,  gsize=iag+1, start=startz )
       call file_history_set_axis( 'FDY',  'Grid distance Y',              'm', 'FY',                    fdyg,  gsize=jag+1, start=startz )
       call file_history_set_axis( 'CBFX', 'Boundary factor Center X',     '1', 'CX', atmos_grid_cartesc_cbfxg, gsize=iag,   start=startz )
       call file_history_set_axis( 'CBFY', 'Boundary factor Center Y',     '1', 'CY', atmos_grid_cartesc_cbfyg, gsize=jag,   start=startz )
       call file_history_set_axis( 'FBFX', 'Boundary factor Face X',       '1', 'FX', atmos_grid_cartesc_fbfxg, gsize=iag+1, start=startz )
       call file_history_set_axis( 'FBFY', 'Boundary factor Face Y',       '1', 'FY', atmos_grid_cartesc_fbfyg, gsize=jag+1, start=startz )
    else
       call file_history_set_axis( 'CX',   'Atmos Grid Center Position X', 'm', 'CX', atmos_grid_cartesc_cx   )
       call file_history_set_axis( 'CY',   'Atmos Grid Center Position Y', 'm', 'CY', atmos_grid_cartesc_cy   )
       call file_history_set_axis( 'FX',   'Atmos Grid Face Position X',   'm', 'FX', atmos_grid_cartesc_fx   )
       call file_history_set_axis( 'FY',   'Atmos Grid Face Position Y',   'm', 'FY', atmos_grid_cartesc_fy   )
       call file_history_set_axis( 'CDX',  'Grid Cell length X',           'm', 'CX', atmos_grid_cartesc_cdx  )
       call file_history_set_axis( 'CDY',  'Grid Cell length Y',           'm', 'CY', atmos_grid_cartesc_cdy  )
       call file_history_set_axis( 'FDX',  'Grid distance X',              'm', 'FX',                    fdx  )
       call file_history_set_axis( 'FDY',  'Grid distance Y',              'm', 'FY',                    fdy  )
       call file_history_set_axis( 'CBFX', 'Boundary factor Center X',     '1', 'CX', atmos_grid_cartesc_cbfx )
       call file_history_set_axis( 'CBFY', 'Boundary factor Center Y',     '1', 'CY', atmos_grid_cartesc_cbfy )
       call file_history_set_axis( 'FBFX', 'Boundary factor Face X',       '1', 'FX', atmos_grid_cartesc_fbfx )
       call file_history_set_axis( 'FBFY', 'Boundary factor Face Y',       '1', 'FY', atmos_grid_cartesc_fbfy )
    endif

    call file_history_set_axis('CXG',   'Grid Center Position X (global)',   'm', 'CXG',  atmos_grid_cartesc_cxg,   gsize=iag,   start=startz )
    call file_history_set_axis('CYG',   'Grid Center Position Y (global)',   'm', 'CYG',  atmos_grid_cartesc_cyg,   gsize=jag,   start=startz )
    call file_history_set_axis('FXG',   'Grid Face Position X (global)',     'm', 'FXG',  atmos_grid_cartesc_fxg,   gsize=iag+1, start=startz )
    call file_history_set_axis('FYG',   'Grid Face Position Y (global)',     'm', 'FYG',  atmos_grid_cartesc_fyg,   gsize=jag+1, start=startz )
    call file_history_set_axis('CDXG',  'Grid Cell length X (global)',       'm', 'CXG',  atmos_grid_cartesc_cdxg,  gsize=iag,   start=startz )
    call file_history_set_axis('CDYG',  'Grid Cell length Y (global)',       'm', 'CYG',  atmos_grid_cartesc_cdyg,  gsize=jag,   start=startz )
    call file_history_set_axis('FDXG',  'Grid distance X (global)',          'm', 'FDXG',                    fdxg,  gsize=iag+1, start=startz )
    call file_history_set_axis('FDYG',  'Grid distance Y (global)',          'm', 'FDYG',                    fdyg,  gsize=jag+1, start=startz )
    call file_history_set_axis('CBFXG', 'Boundary factor Center X (global)', '1', 'CXG',  atmos_grid_cartesc_cbfxg, gsize=iag,   start=startz )
    call file_history_set_axis('CBFYG', 'Boundary factor Center Y (global)', '1', 'CYG',  atmos_grid_cartesc_cbfyg, gsize=jag,   start=startz )
    call file_history_set_axis('FBFXG', 'Boundary factor Face X (global)',   '1', 'FXG',  atmos_grid_cartesc_fbfxg, gsize=iag+1, start=startz )
    call file_history_set_axis('FBFYG', 'Boundary factor Face Y (global)',   '1', 'FYG',  atmos_grid_cartesc_fbfyg, gsize=jag+1, start=startz )



    ! associate coordinates
    if ( file_history_aggregate ) then
       if ( file_history_cartesc_boundary ) then
          start(1,:) = isgb   ! global subarray starting index
          start(2,:) = jsgb   ! global subarray starting index
       else
          start(1,:) = prc_2drank(prc_myrank,1) * imax + 1 ! no IHALO
          start(2,:) = prc_2drank(prc_myrank,2) * jmax + 1 ! no JHALO
          if ( (.NOT. prc_periodic_x) .AND. prc_has_w ) then
             start(1,2) = start(1,2) + 1
             start(1,4) = start(1,4) + 1
          endif
          if ( (.NOT. prc_periodic_y) .AND. prc_has_s ) then
             start(2,3) = start(2,3) + 1
             start(2,4) = start(2,4) + 1
          endif
       endif
       start(3,:) = 1
    else
       start(:,:) = 1
    endif

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

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

    do k = 1, kmax
    do j = 1, jm
    do i = 1, min(imh,ia-imsh)
       axis(i,j,k) = ( atmos_grid_cartesc_real_cz(k+ks-1,imsh+i-1,jms+j-1) + atmos_grid_cartesc_real_cz(k+ks-1,imsh+i,jms+j-1) ) * 0.5_rp
    enddo
    enddo
    enddo
    if ( imh == ia-imsh+1 ) then
       do k = 1, kmax
       do j = 1, jm
          axis(imh,j,k) = atmos_grid_cartesc_real_cz(k+ks-1,imsh+imh-1,jms+j-1)
       enddo
       enddo
    endif
    axis_name(1:3) = (/'xh', 'y ', 'z '/)
    call file_history_set_associatedcoordinate( 'height_uyz', 'height above ground level (half level uyz)',    &
                                                'm', axis_name(1:3), axis(1:imh,1:jm,1:kmax), start=start(:,2) )

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

    do k = 1, kmax
    do j = 1, min(jmh,ja-jmsh)
    do i = 1, min(imh,ia-imsh)
       axis(i,j,k) = ( atmos_grid_cartesc_real_cz(k+ks-1,imsh+i-1,jmsh+j-1) + atmos_grid_cartesc_real_cz(k+ks-1,imsh+i  ,jmsh+j-1) &
                     + atmos_grid_cartesc_real_cz(k+ks-1,imsh+i-1,jmsh+j  ) + atmos_grid_cartesc_real_cz(k+ks-1,imsh+i  ,jmsh+j  ) ) * 0.25_rp
    enddo
    enddo
    enddo
    if ( jmh == ja-jmsh+1 ) then
       do k = 1, kmax
       do i = 1, min(imh,ia-imsh)
          axis(i,jmh,k) = ( atmos_grid_cartesc_real_cz(k+ks-1,imsh+i-1,jmsh+jmh-1) + atmos_grid_cartesc_real_cz(k+ks-1,imsh+i,jmsh+jmh-1) ) * 0.5_rp
       enddo
       enddo
    endif
    if ( imh == ia-imsh+1 ) then
       do k = 1, kmax
       do j = 1, min(jmh,ja-jmsh)
          axis(imh,j,k) = ( atmos_grid_cartesc_real_cz(k+ks-1,imsh+imh-1,jmsh+j-1) + atmos_grid_cartesc_real_cz(k+ks-1,imsh+imh-1,jmsh+j) ) * 0.5_rp
       enddo
       enddo
    endif
    if ( imh == ia-imsh+1 .AND. jmh == ja-jmsh+1 ) then
       do k = 1, kmax
          axis(imh,jmh,k) = atmos_grid_cartesc_real_cz(k+ks-1,imsh+imh-1,jmsh+jmh-1)
       enddo
    endif
    axis_name(1:3) = (/'xh', 'yh', 'z '/)
    call file_history_set_associatedcoordinate( 'height_uvz', 'height above ground level (half level uvz)',     &
                                                'm', axis_name(1:3), axis(1:imh,1:jmh,1:kmax), start=start(:,4) )

    do k = 0, kmax
    do j = 1, jm
    do i = 1, min(imh,ia-imsh)
       axis(i,j,k) = ( atmos_grid_cartesc_real_fz(k+ks-1,imsh+i-1,jms+j-1) + atmos_grid_cartesc_real_fz(k+ks-1,imsh+i,jms+j-1) ) * 0.5_rp
    enddo
    enddo
    enddo
    if ( imh == ia-imsh+1 ) then
       do k = 0, kmax
       do j = 1, jm
          axis(imh,j,k) = atmos_grid_cartesc_real_fz(k+ks-1,imsh+imh-1,jms+j-1)
       enddo
       enddo
    endif
    axis_name(1:3) = (/'xh', 'y ', 'zh'/)
    call file_history_set_associatedcoordinate( 'height_uyw', 'height above ground level (half level uyw)',    &
                                                'm', axis_name(1:3), axis(1:imh,1:jm,0:kmax), start=start(:,2) )

    do k = 0, kmax
    do j = 1, min(jmh,ja-jmsh)
    do i = 1, im
       axis(i,j,k) = ( atmos_grid_cartesc_real_fz(k+ks-1,ims+i-1,jmsh+j-1) + atmos_grid_cartesc_real_fz(k+ks-1,ims+i-1,jmsh+j) ) * 0.5_rp
    enddo
    enddo
    enddo
    if ( jmh == ja-jmsh+1 ) then
       do k = 0, kmax
       do i = 1, im
          axis(i,jmh,k) = atmos_grid_cartesc_real_fz(k+ks-1,ims+i-1,jmsh+jmh-1)
       enddo
       enddo
    endif
    axis_name(1:3) = (/'x ', 'yh', 'zh'/)
    call file_history_set_associatedcoordinate( 'height_xvw', 'height above ground level (half level xvw)',    &
                                                'm', axis_name(1:3), axis(1:im,1:jmh,0:kmax), start=start(:,3) )

    do k = 0, kmax
    do j = 1, min(jmh,ja-jmsh)
    do i = 1, min(imh,ia-imsh)
       axis(i,j,k) = ( atmos_grid_cartesc_real_fz(k+ks-1,imsh+i-1,jmsh+j-1) + atmos_grid_cartesc_real_fz(k+ks-1,imsh+i  ,jmsh+j-1) &
                     + atmos_grid_cartesc_real_fz(k+ks-1,imsh+i-1,jmsh+j  ) + atmos_grid_cartesc_real_fz(k+ks-1,imsh+i  ,jmsh+j  ) ) * 0.25_rp
    enddo
    enddo
    enddo
    if ( jmh == ja-jmsh+1 ) then
       do k = 0, kmax
       do i = 1, min(imh,ia-imsh)
          axis(i,jmh,k) = ( atmos_grid_cartesc_real_fz(k+ks-1,imsh+i-1,jmsh+jmh-1) + atmos_grid_cartesc_real_fz(k+ks-1,imsh+i,jmsh+jmh-1) ) * 0.5_rp
       enddo
       enddo
    endif
    if ( imh == ia-imsh+1 ) then
       do k = 0, kmax
       do j = 1, min(jmh,ja-jmsh)
          axis(imh,j,k) = ( atmos_grid_cartesc_real_fz(k+ks-1,imsh+imh-1,jmsh+j-1) + atmos_grid_cartesc_real_fz(k+ks-1,imsh+imh-1,jmsh+j) ) * 0.5_rp
       enddo
       enddo
    endif
    if ( imh == ia-imsh+1 .AND. jm == ja-jmsh+1 ) then
       do k = 0, kmax
          axis(imh,jmh,k) = atmos_grid_cartesc_real_fz(k+ks-1,imsh+imh-1,jmsh+jmh-1)
       enddo
    endif
    axis_name(1:3) = (/'xh', 'yh', 'zh'/)
    call file_history_set_associatedcoordinate( 'height_uvw', 'height above ground level (half level uvw)', 'm', &
                                                axis_name(1:3), axis(1:imh,1:jmh,0:kmax), start=start(:,4)       )

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

    axis(1:imh,1:jm,1) = atmos_grid_cartesc_real_lonuy(imsh:ime,jms:jme) / d2r
    axis_name(1:2) = (/'xh', 'y '/)
    call file_history_set_associatedcoordinate( 'lon_uy', 'longitude (half level uy)', 'degrees_east', &
                                                axis_name(1:2), axis(1:imh,1:jm,1), start=start(:,2)   )

    axis(1:im,1:jmh,1) = atmos_grid_cartesc_real_lonxv(ims:ime,jmsh:jme) / d2r
    axis_name(1:2) = (/'x ', 'yh'/)
    call file_history_set_associatedcoordinate( 'lon_xv', 'longitude (half level xv)', 'degrees_east', &
                                                axis_name(1:2), axis(1:im,1:jmh,1), start=start(:,3)   )

    axis(1:imh,1:jmh,1) = atmos_grid_cartesc_real_lonuv(imsh:ime,jmsh:jme) / d2r
    axis_name(1:2) = (/'xh', 'yh'/)
    call file_history_set_associatedcoordinate( 'lon_uv', 'longitude (half level uv)', 'degrees_east', &
                                                axis_name(1:2), axis(1:imh,1:jmh,1), start=start(:,4)  )

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

    axis(1:imh,1:jm,1) = atmos_grid_cartesc_real_latuy(imsh:ime,jms:jme) / d2r
    axis_name(1:2) = (/'xh', 'y '/)
    call file_history_set_associatedcoordinate( 'lat_uy', 'latitude (half level uy)',  'degrees_north', &
                                                axis_name(1:2), axis(1:imh,1:jm,1), start=start(:,2)    )

    axis(1:im,1:jmh,1) = atmos_grid_cartesc_real_latxv(ims:ime,jmsh:jme) / d2r
    axis_name(1:2) = (/'x ', 'yh'/)
    call file_history_set_associatedcoordinate( 'lat_xv', 'latitude (half level xv)',  'degrees_north', &
                                                axis_name(1:2), axis(1:im,1:jmh,1), start=start(:,3)    )

    axis(1:imh,1:jmh,1) = atmos_grid_cartesc_real_latuv(imsh:ime,jmsh:jme) / d2r
    axis_name(1:2) = (/'xh', 'yh'/)
    call file_history_set_associatedcoordinate( 'lat_uv', 'latitude (half level uv)',  'degrees_north', &
                                                axis_name(1:2), axis(1:imh,1:jmh,1), start=start(:,4)   )

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

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

    axis_name(1:2) = (/'x ','y '/)
    call file_history_set_associatedcoordinate( 'cell_area',    'area of grid cell',                  'm2', axis_name(1:2), &
                                                area(ims:ime,jms:jme),   start=start(:,1)                                   )
    axis_name(1:2) = (/'xh','y '/)
    call file_history_set_associatedcoordinate( 'cell_area_uy', 'area of grid cell (half level uy)',  'm2', axis_name(1:2), &
                                                areauy(imsh:ime,jms:jme), start=start(:,2)                                  )
    axis_name(1:2) = (/'x ','yh'/)
    call file_history_set_associatedcoordinate( 'cell_area_xv', 'area of grid cell (half level xv)',  'm2', axis_name(1:2), &
                                                areaxv(ims:ime,jmsh:jme), start=start(:,3)                                  )

    do k = 1, kmax
    do j = 1, jm
    do i = 1, imh
       axis(i,j,k) = areazuy_x(ks+k-1,imsh+i-1,jms+j-1)
    end do
    end do
    end do
    axis_name = (/'xh', 'y ', 'z '/)
    call file_history_set_associatedcoordinate( 'cell_area_uyz_x', 'area of grid cell face (half level uyz, normal x)', 'm2', &
                                                axis_name(1:3), axis(1:imh,1:jm,1:kmax), start=start(:,2)                     )
    do k = 1, kmax
    do j = 1, jmh
    do i = 1, im
       axis(i,j,k) = areazxv_y(ks+k-1,ims+i-1,jmsh+j-1)
    end do
    end do
    end do
    axis_name = (/'x ', 'yh', 'z '/)
    call file_history_set_associatedcoordinate( 'cell_area_xvz_y', 'area of grid cell face (half level xvz, normal y)', 'm2', &
                                                axis_name(1:3), axis(1:im,1:jmh,1:kmax), start=start(:,3)                     )
    do k = 0, kmax
    do j = 1, jmh
    do i = 1, imh
       axis(i,j,k) = areawuy_x(ks+k-1,imsh+i-1,jmsh+j-1)
    end do
    end do
    end do
    axis_name = (/'xh', 'y ', 'zh'/)
    call file_history_set_associatedcoordinate( 'cell_area_uyw_x', 'area of grid cell face (half level uyw, normal x)', 'm2', &
                                                axis_name(1:3), axis(1:imh,1:jmh,0:kmax), start=start(:,2)                    )
    do k = 0, kmax
    do j = 1, jmh
    do i = 1, im
       axis(i,j,k) = areawxv_y(ks+k-1,ims+i-1,jmsh+j-1)
    end do
    end do
    end do
    axis_name = (/'x ', 'yh', 'zh'/)
    call file_history_set_associatedcoordinate( 'cell_area_xvw_y', 'area of grid cell face (half level xvw, normal y)', 'm2', &
                                                axis_name(1:3), axis(1:im,1:jmh,0:kmax), start=start(:,3)                     )
    do k = 1, kmax
    do j = 1, jm
    do i = 1, im
       axis(i,j,k) = areazxy_x(ks+k-1,ims+i-1,jms+j-1)
    end do
    end do
    end do
    axis_name = (/'x ', 'y ', 'z '/)
    call file_history_set_associatedcoordinate( 'cell_area_xyz_x', 'area of grid cell face (half level xyz, normal x)', 'm2', &
                                                axis_name(1:3), axis(1:im,1:jm,1:kmax), start=start(:,1)                      )
    do k = 1, kmax
    do j = 1, jmh
    do i = 1, imh
       axis(i,j,k) = areazuv_y(ks+k-1,imsh+i-1,jmsh+j-1)
    end do
    end do
    end do
    axis_name = (/'xh', 'yh', 'z '/)
    call file_history_set_associatedcoordinate( 'cell_area_uvz_y', 'area of grid cell face (half level uvz, normal y)', 'm2', &
                                                axis_name(1:3), axis(1:imh,1:jmh,1:kmax), start=start(:,4)                    )
    do k = 1, kmax
    do j = 1, jmh
    do i = 1, imh
       axis(i,j,k) = areazuv_x(ks+k-1,imsh+i-1,jmsh+j-1)
    end do
    end do
    end do
    axis_name = (/'xh', 'yh', 'z '/)
    call file_history_set_associatedcoordinate( 'cell_area_uvz_x', 'area of grid cell face (half level uvz, normal x)', 'm2', &
                                                axis_name(1:3), axis(1:imh,1:jmh,1:kmax), start=start(:,4)                    )
    do k = 1, kmax
    do j = 1, jm
    do i = 1, im
       axis(i,j,k) = areazxy_y(ks+k-1,ims+i-1,jms+j-1)
    end do
    end do
    end do
    axis_name = (/'x ', 'y ', 'z '/)
    call file_history_set_associatedcoordinate( 'cell_area_xyz_y', 'area of grid cell face (half level xyz, normal y)', 'm2', &
                                                axis_name(1:3), axis(1:im,1:jm,1:kmax), start=start(:,1)                      )

    do k = 1, kmax
    do j = 1, jm
    do i = 1, im
       axis(i,j,k) = vol(ks+k-1,ims+i-1,jms+j-1)
    end do
    end do
    end do
    axis_name = (/ 'x ', 'y ', 'z '/)
    call file_history_set_associatedcoordinate( 'cell_volume',     'volume of grid cell',                  'm3', &
                                                axis_name(1:3), axis(1:im,1:jm,1:kmax), start=start(:,1)         )
    do k = 0, kmax
    do j = 1, jm
    do i = 1, im
       axis(i,j,k) = volwxy(ks+k-1,ims+i-1,jms+j-1)
    end do
    end do
    end do
    axis_name = (/'x ', 'y ', 'zh'/)
    call file_history_set_associatedcoordinate( 'cell_volume_xyw', 'volume of grid cell (half level xyw)', 'm3', &
                                                axis_name(1:3), axis(1:im,1:jm,0:kmax), start=start(:,1)         )
    do k = 1, kmax
    do j = 1, jm
    do i = 1, imh
       axis(i,j,k) = volzuy(ks+k-1,imsh+i-1,jms+j-1)
    end do
    end do
    end do
    axis_name = (/'xh', 'y ', 'z '/)
    call file_history_set_associatedcoordinate( 'cell_volume_uyz', 'volume of grid cell (half level uyz)',  'm3', &
                                                axis_name(1:3), axis(1:imh,1:jm,1:kmax), start=start(:,2)         )
    do k = 1, kmax
    do j = 1, jmh
    do i = 1, im
       axis(i,j,k) = volzxv(ks+k-1,ims+i-1,jmsh+j-1)
    end do
    end do
    end do
    axis_name = (/'x ', 'yh', 'z '/)
    call file_history_set_associatedcoordinate( 'cell_volume_xvz', 'volume of grid cell (half level xvz)',  'm3', &
                                                axis_name(1:3), axis(1:im,1:jmh,1:kmax), start=start(:,3)         )

    if ( okmax > 0 ) then
       do k = 1, okmax
       do j = 1, jm
       do i = 1, im
          axiso(i,j,k) = volo(oks+k-1,ims+i-1,jms+j-1)
       end do
       end do
       end do
       axis_name = (/'x ', 'y ', 'oz'/)
       call file_history_set_associatedcoordinate( 'cell_volume_xyo', 'volume of grid cell', 'm3', &
                                                   axis_name(1:3), axiso(:,:,:), start=start(:,1)  )
    endif

    if ( lkmax > 0 ) then
       do k = 1, lkmax
       do j = 1, jm
       do i = 1, im
          axisl(i,j,k) = voll(lks+k-1,ims+i-1,jms+j-1)
       end do
       end do
       end do
       axis_name = (/'x ', 'y ', 'lz'/)
       call file_history_set_associatedcoordinate( 'cell_volume_xyl', 'volume of grid cell', 'm3', &
                                                   axis_name(1:3), axisl(:,:,:), start=start(:,1)  )
    endif

    if ( ukmax > 0 ) then
       do k = 1, ukmax
       do j = 1, jm
       do i = 1, im
          axisu(i,j,k) = volu(uks+k-1,ims+i-1,jms+j-1)
       end do
       end do
       end do
       axis_name = (/'x ', 'y ', 'uz'/)
       call file_history_set_associatedcoordinate( 'cell_volume_xyu', 'volume of grid cell', 'm3', &
                                                   axis_name(1:3), axisu(:,:,:), start=start(:,1)  )
    endif

    return
  end subroutine file_history_cartesc_set_axes

  !-----------------------------------------------------------------------------
  subroutine file_history_cartesc_set_axes_attributes
    use scale_atmos_grid_cartesc, only: &
       atmos_grid_cartesc_name
    use scale_calendar, only: &
       calendar_get_name
    use scale_file_history, only: &
       file_history_aggregate, &
       file_history_set_attribute
    use scale_prc, only: &
       prc_myrank
    use scale_const, only: &
       undef => const_undef
    use scale_prc_cartesc, only: &
       prc_2drank,     &
       prc_num_x,      &
       prc_num_y,      &
       prc_periodic_x, &
       prc_periodic_y, &
       prc_has_w,      &
       prc_has_e,      &
       prc_has_s,      &
       prc_has_n
    use scale_file, only: &
       file_get_cftunits
    use scale_file_cartesc, only: &
       axisattinfo,        &
       mappinginfo
    use scale_mapprojection, only: &
       mapprojection_get_attributes
    implicit none

    character(len=34) :: tunits
    character(len=H_SHORT) :: calendar

    type(axisattinfo) :: ainfo(4) ! x, xh, y, yh
    type(mappinginfo) :: minfo
    !---------------------------------------------------------------------------

    call file_history_set_attribute( "global", "Conventions", "CF-1.6" ) ! [IN]

    call file_history_set_attribute( "global", "grid_name", atmos_grid_cartesc_name ) ! [IN]

    if ( file_history_aggregate ) then
       call file_history_set_attribute( "global", "scale_cartesC_prc_rank_x", (/0/) ) ! [IN]
       call file_history_set_attribute( "global", "scale_cartesC_prc_rank_y", (/0/) ) ! [IN]

       call file_history_set_attribute( "global", "scale_cartesC_prc_num_x", (/1/) ) ! [IN]
       call file_history_set_attribute( "global", "scale_cartesC_prc_num_y", (/1/) ) ! [IN]
    else
       call file_history_set_attribute( "global", "scale_cartesC_prc_rank_x", (/prc_2drank(prc_myrank,1)/) ) ! [IN]
       call file_history_set_attribute( "global", "scale_cartesC_prc_rank_y", (/prc_2drank(prc_myrank,2)/) ) ! [IN]

       call file_history_set_attribute( "global", "scale_cartesC_prc_num_x", (/prc_num_x/) ) ! [IN]
       call file_history_set_attribute( "global", "scale_cartesC_prc_num_y", (/prc_num_y/) ) ! [IN]
    endif

    call file_history_set_attribute( "global", "scale_cartesC_prc_periodic_z", .false.        ) ! [IN]
    call file_history_set_attribute( "global", "scale_cartesC_prc_periodic_x", prc_periodic_x ) ! [IN]
    call file_history_set_attribute( "global", "scale_cartesC_prc_periodic_y", prc_periodic_y ) ! [IN]

    call file_history_set_attribute( "global", "scale_atmos_grid_cartesC_index_imaxg", (/imaxg/) ) ! [IN]
    call file_history_set_attribute( "global", "scale_atmos_grid_cartesC_index_jmaxg", (/jmaxg/) ) ! [IN]

                     call file_history_set_attribute( "global", "scale_atmos_grid_cartesC_index_kmax", (/kmax/)  ) ! [IN]
    if ( okmax > 0 ) call file_history_set_attribute( "global", "scale_ocean_grid_cartesC_index_kmax", (/okmax/) ) ! [IN]
    if ( lkmax > 0 ) call file_history_set_attribute( "global", "scale_land_grid_cartesC_index_kmax",  (/lkmax/) ) ! [IN]
    if ( ukmax > 0 ) call file_history_set_attribute( "global", "scale_urban_grid_cartesC_index_kmax", (/ukmax/) ) ! [IN]

    call file_history_set_attribute( "global", "scale_atmos_grid_cartesC_index_khalo", (/khalo/) ) ! [IN]
    call file_history_set_attribute( "global", "scale_atmos_grid_cartesC_index_ihalo", (/ihalo/) ) ! [IN]
    call file_history_set_attribute( "global", "scale_atmos_grid_cartesC_index_jhalo", (/jhalo/) ) ! [IN]

    call calendar_get_name( calendar )
    if ( calendar /= "" ) call file_history_set_attribute( "global", "calendar", calendar )
    call file_get_cftunits( file_history_cartesc_startdate(:), tunits )
    call file_history_set_attribute( "global", "time_units", tunits )
    call file_history_set_attribute( "global", "time_start", (/file_history_cartesc_startms/) )

    if ( prc_periodic_x ) then
       ainfo(1)%periodic = .true.
       ainfo(2)%periodic = .true.
    else
       ainfo(1)%periodic = .false.
       ainfo(2)%periodic = .false.
    endif

    if ( prc_periodic_y ) then
       ainfo(3)%periodic = .true.
       ainfo(4)%periodic = .true.
    else
       ainfo(3)%periodic = .false.
       ainfo(4)%periodic = .false.
    endif

    ! for x
    if ( prc_periodic_x .OR. .NOT. file_history_cartesc_boundary ) then
       ainfo(1)%size_global (1) = imax * prc_num_x
       ainfo(1)%start_global(1) = is_ing - ihalo
       ainfo(1)%halo_global (1) = 0     ! west side
       ainfo(1)%halo_global (2) = 0     ! east side
       ainfo(1)%halo_local  (1) = 0     ! west side
       ainfo(1)%halo_local  (2) = 0     ! east side
    else
       ainfo(1)%size_global (1) = iag
       ainfo(1)%start_global(1) = isga
       ainfo(1)%halo_global (1) = ihalo ! west side
       ainfo(1)%halo_global (2) = ihalo ! east side
       ainfo(1)%halo_local  (1) = ihalo ! west side
       ainfo(1)%halo_local  (2) = ihalo ! east side
       if ( .not. file_history_aggregate ) then
          if( prc_has_w ) ainfo(1)%halo_local(1) = 0
          if( prc_has_e ) ainfo(1)%halo_local(2) = 0
       end if
    endif

    ! for xh
    ainfo(2) = ainfo(1)
    if ( .NOT. prc_periodic_x .AND. .NOT. file_history_cartesc_boundary ) then
       ainfo(2)%size_global (1) = ainfo(2)%size_global (1) + 1
       ainfo(2)%halo_global (1) = ainfo(2)%halo_global (1) + 1
       if ( prc_has_w .and. ( .not. file_history_aggregate ) ) then
          ainfo(2)%start_global(1) = ainfo(2)%start_global(1) + 1
       else
          ainfo(2)%halo_local  (1) = ainfo(2)%halo_local  (1) + 1
       endif
    endif

    ! for y
    if ( prc_periodic_y .OR. .NOT. file_history_cartesc_boundary ) then
       ainfo(3)%size_global (1) = jmax * prc_num_y
       ainfo(3)%start_global(1) = js_ing - jhalo
       ainfo(3)%halo_global (1) = 0     ! south side
       ainfo(3)%halo_global (2) = 0     ! north side
       ainfo(3)%halo_local  (1) = 0     ! south side
       ainfo(3)%halo_local  (2) = 0     ! north side
    else
       ainfo(3)%size_global (1) = jag
       ainfo(3)%start_global(1) = jsga
       ainfo(3)%halo_global (1) = jhalo ! south side
       ainfo(3)%halo_global (2) = jhalo ! north side
       ainfo(3)%halo_local  (1) = jhalo ! south side
       ainfo(3)%halo_local  (2) = jhalo ! north side
       if ( .not. file_history_aggregate ) then
          if( prc_has_s ) ainfo(3)%halo_local(1) = 0
          if( prc_has_n ) ainfo(3)%halo_local(2) = 0
       end if
    endif

    ! for yh
    ainfo(4) = ainfo(3)
    if ( .NOT. prc_periodic_y .AND. .NOT. file_history_cartesc_boundary ) then
       ainfo(4)%size_global (1) = ainfo(4)%size_global (1) + 1
       ainfo(4)%halo_global (1) = ainfo(4)%halo_global (1) + 1
       if ( prc_has_s .and. ( .not. file_history_aggregate ) ) then
          ainfo(4)%start_global(1) = ainfo(4)%start_global(1) + 1
       else
          ainfo(4)%halo_local  (1) = ainfo(4)%halo_local  (1) + 1
       endif
    endif

    if ( file_history_aggregate ) then
       ainfo(1)%start_global(1) = 1
       ainfo(2)%start_global(1) = 1
       ainfo(3)%start_global(1) = 1
       ainfo(4)%start_global(1) = 1
    end if

    call file_history_set_attribute( "x" , "size_global" , ainfo(1)%size_global (:) )
    call file_history_set_attribute( "x" , "start_global", ainfo(1)%start_global(:) )
    call file_history_set_attribute( "x" , "halo_global" , ainfo(1)%halo_global (:) )
    call file_history_set_attribute( "x" , "halo_local"  , ainfo(1)%halo_local  (:) )
    call file_history_set_attribute( "x" , "periodic"    , ainfo(1)%periodic        )

    call file_history_set_attribute( "xh", "size_global" , ainfo(2)%size_global (:) )
    call file_history_set_attribute( "xh", "start_global", ainfo(2)%start_global(:) )
    call file_history_set_attribute( "xh", "halo_global" , ainfo(2)%halo_global (:) )
    call file_history_set_attribute( "xh", "halo_local"  , ainfo(2)%halo_local  (:) )
    call file_history_set_attribute( "xh", "periodic"    , ainfo(2)%periodic        )

    call file_history_set_attribute( "y" , "size_global" , ainfo(3)%size_global (:) )
    call file_history_set_attribute( "y" , "start_global", ainfo(3)%start_global(:) )
    call file_history_set_attribute( "y" , "halo_global" , ainfo(3)%halo_global (:) )
    call file_history_set_attribute( "y" , "halo_local"  , ainfo(3)%halo_local  (:) )
    call file_history_set_attribute( "y" , "periodic"    , ainfo(3)%periodic        )

    call file_history_set_attribute( "yh", "size_global" , ainfo(4)%size_global (:) )
    call file_history_set_attribute( "yh", "start_global", ainfo(4)%start_global(:) )
    call file_history_set_attribute( "yh", "halo_global" , ainfo(4)%halo_global (:) )
    call file_history_set_attribute( "yh", "halo_local"  , ainfo(4)%halo_local  (:) )
    call file_history_set_attribute( "yh", "periodic"    , ainfo(4)%periodic        )

    ! map projection info
    call mapprojection_get_attributes( minfo%mapping_name,                             & ! [OUT]
                                       minfo%false_easting                        (1), & ! [OUT]
                                       minfo%false_northing                       (1), & ! [OUT]
                                       minfo%longitude_of_central_meridian        (1), & ! [OUT]
                                       minfo%longitude_of_projection_origin       (1), & ! [OUT]
                                       minfo%latitude_of_projection_origin        (1), & ! [OUT]
                                       minfo%straight_vertical_longitude_from_pole(1), & ! [OUT]
                                       minfo%standard_parallel                    (:)  ) ! [OUT]

    if ( minfo%mapping_name /= "" ) then
       call file_history_set_attribute( "x" , "standard_name", "projection_x_coordinate" )
       call file_history_set_attribute( "xh", "standard_name", "projection_x_coordinate" )
       call file_history_set_attribute( "y" , "standard_name", "projection_y_coordinate" )
       call file_history_set_attribute( "yh", "standard_name", "projection_y_coordinate" )

       call file_history_set_attribute( minfo%mapping_name, "grid_mapping_name", minfo%mapping_name, add_variable=.true. )

       if ( minfo%false_easting(1) /= undef ) then
          call file_history_set_attribute( minfo%mapping_name,    & ! [IN]
                                           "false_easting",       & ! [IN]
                                           minfo%false_easting(:) ) ! [IN]
       endif

       if ( minfo%false_northing(1) /= undef ) then
          call file_history_set_attribute( minfo%mapping_name,     & ! [IN]
                                           "false_northing",       & ! [IN]
                                           minfo%false_northing(:) ) ! [IN]
       endif

       if ( minfo%longitude_of_central_meridian(1) /= undef ) then
          call file_history_set_attribute( minfo%mapping_name,                    & ! [IN]
                                           "longitude_of_central_meridian",       & ! [IN]
                                           minfo%longitude_of_central_meridian(:) ) ! [IN]
       endif

       if ( minfo%longitude_of_projection_origin(1) /= undef ) then
          call file_history_set_attribute( minfo%mapping_name,                     & ! [IN]
                                           "longitude_of_projection_origin",       & ! [IN]
                                           minfo%longitude_of_projection_origin(:) ) ! [IN]
       endif

       if ( minfo%latitude_of_projection_origin(1) /= undef ) then
          call file_history_set_attribute( minfo%mapping_name,                    & ! [IN]
                                           "latitude_of_projection_origin",       & ! [IN]
                                           minfo%latitude_of_projection_origin(:) ) ! [IN]
       endif

       if ( minfo%straight_vertical_longitude_from_pole(1) /= undef ) then
          call file_history_set_attribute( minfo%mapping_name,                            & ! [IN]
                                           "straight_vertical_longitude_from_pole",       & ! [IN]
                                           minfo%straight_vertical_longitude_from_pole(:) ) ! [IN]
       endif

       if ( minfo%standard_parallel(1) /= undef ) then
          if ( minfo%standard_parallel(2) /= undef ) then
             call file_history_set_attribute( minfo%mapping_name,          & ! [IN]
                                              "standard_parallel",         & ! [IN]
                                              minfo%standard_parallel(1:2) ) ! [IN]
          else
             call file_history_set_attribute( minfo%mapping_name,          & ! [IN]
                                              "standard_parallel",         & ! [IN]
                                              minfo%standard_parallel(1:1) ) ! [IN]
          endif
       endif
    endif

    ! area and volume
    call file_history_set_attribute( "cell_area",    "standard_name", "area" ) ! [IN]
    call file_history_set_attribute( "cell_area_uy", "standard_name", "area" ) ! [IN]
    call file_history_set_attribute( "cell_area_xv", "standard_name", "area" ) ! [IN]

    call file_history_set_attribute( "cell_area_uyz_x", "standard_name", "area" ) ! [IN]
    call file_history_set_attribute( "cell_area_xvz_y", "standard_name", "area" ) ! [IN]
    call file_history_set_attribute( "cell_area_uyw_x", "standard_name", "area" ) ! [IN]
    call file_history_set_attribute( "cell_area_xvw_y", "standard_name", "area" ) ! [IN]
    call file_history_set_attribute( "cell_area_xyz_x", "standard_name", "area" ) ! [IN]
    call file_history_set_attribute( "cell_area_uvz_y", "standard_name", "area" ) ! [IN]
    call file_history_set_attribute( "cell_area_uvz_x", "standard_name", "area" ) ! [IN]
    call file_history_set_attribute( "cell_area_xyz_y", "standard_name", "area" ) ! [IN]

    call file_history_set_attribute( "cell_volume",     "standard_name", "volume" ) ! [IN]
    call file_history_set_attribute( "cell_volume_xyw", "standard_name", "volume" ) ! [IN]
    call file_history_set_attribute( "cell_volume_uyz", "standard_name", "volume" ) ! [IN]
    call file_history_set_attribute( "cell_volume_xvz", "standard_name", "volume" ) ! [IN]

    if ( okmax > 0 ) then
       call file_history_set_attribute( "cell_volume_xyo", "standard_name", "volume" ) ! [IN]
    endif
    if ( lkmax > 0 ) then
       call file_history_set_attribute( "cell_volume_xyl", "standard_name", "volume" ) ! [IN]
    endif
    if ( ukmax > 0 ) then
       call file_history_set_attribute( "cell_volume_xyu", "standard_name", "volume" ) ! [IN]
    endif

    ! SGRID
    call file_history_set_attribute( "grid", "cf_role",             "grid_topology", add_variable=.true. )
    call file_history_set_attribute( "grid", "topology_dimension",  (/ 2 /) )
    call file_history_set_attribute( "grid", "node_dimensions",     "xh yh" )
    call file_history_set_attribute( "grid", "face_dimensions",     "x: xh (padding: none) y: yh (padding: none)" )
    call file_history_set_attribute( "grid", "node_coordinates",    "lon_uv lat_uv" )
    call file_history_set_attribute( "grid", "face_coordinates",    "lon lat" )
    call file_history_set_attribute( "grid", "edge1_coordinates",   "lon_uy lat_uy" )
    call file_history_set_attribute( "grid", "edge2_coordinates",   "lon_xv lat_xv" )
    call file_history_set_attribute( "grid", "vertical_dimensions", "z: zh (padding: none)" )

    call file_history_set_attribute( "grid_ocean", "cf_role",             "grid_topology", add_variable=.true. )
    call file_history_set_attribute( "grid_ocean", "topology_dimension",  (/ 2 /) )
    call file_history_set_attribute( "grid_ocean", "node_dimensions",     "xh yh" )
    call file_history_set_attribute( "grid_ocean", "face_dimensions",     "x: xh (padding: none) y: yh (padding: none)" )
    call file_history_set_attribute( "grid_ocean", "node_coordinates",    "lon_uv lat_uv" )
    call file_history_set_attribute( "grid_ocean", "face_coordinates",    "lon lat" )
    call file_history_set_attribute( "grid_ocean", "edge1_coordinates",   "lon_uy lat_uy" )
    call file_history_set_attribute( "grid_ocean", "edge2_coordinates",   "lon_xv lat_xv" )
    call file_history_set_attribute( "grid_ocean", "vertical_dimensions", "oz: ozh (padding: none)" )

    call file_history_set_attribute( "grid_land", "cf_role",             "grid_topology", add_variable=.true. )
    call file_history_set_attribute( "grid_land", "topology_dimension",  (/ 2 /) )
    call file_history_set_attribute( "grid_land", "node_dimensions",     "xh yh" )
    call file_history_set_attribute( "grid_land", "face_dimensions",     "x: xh (padding: none) y: yh (padding: none)" )
    call file_history_set_attribute( "grid_land", "node_coordinates",    "lon_uv lat_uv" )
    call file_history_set_attribute( "grid_land", "face_coordinates",    "lon lat" )
    call file_history_set_attribute( "grid_land", "edge1_coordinates",   "lon_uy lat_uy" )
    call file_history_set_attribute( "grid_land", "edge2_coordinates",   "lon_xv lat_xv" )
    call file_history_set_attribute( "grid_land", "vertical_dimensions", "lz: lzh (padding: none)" )

    call file_history_set_attribute( "grid_urban", "cf_role",             "grid_topology", add_variable=.true. )
    call file_history_set_attribute( "grid_urban", "topology_dimension",  (/ 2 /) )
    call file_history_set_attribute( "grid_urban", "node_dimensions",     "xh yh" )
    call file_history_set_attribute( "grid_urban", "face_dimensions",     "x: xh (padding: none) y: yh (padding: none)" )
    call file_history_set_attribute( "grid_urban", "node_coordinates",    "lon_uv lat_uv" )
    call file_history_set_attribute( "grid_urban", "face_coordinates",    "lon lat" )
    call file_history_set_attribute( "grid_urban", "edge1_coordinates",   "lon_uy lat_uy" )
    call file_history_set_attribute( "grid_urban", "edge2_coordinates",   "lon_xv lat_xv" )
    call file_history_set_attribute( "grid_urban", "vertical_dimensions", "uz: uzh (padding: none)" )

    call file_history_set_attribute( "grid_pressure", "cf_role",             "grid_topology", add_variable=.true. )
    call file_history_set_attribute( "grid_pressure", "topology_dimension",  (/ 2 /) )
    call file_history_set_attribute( "grid_pressure", "node_dimensions",     "xh yh" )
    call file_history_set_attribute( "grid_pressure", "face_dimensions",     "x: xh (padding: none) y: yh (padding: none)" )
    call file_history_set_attribute( "grid_pressure", "node_coordinates",    "lon_uv lat_uv" )
    call file_history_set_attribute( "grid_pressure", "face_coordinates",    "lon lat" )
    call file_history_set_attribute( "grid_pressure", "edge1_coordinates",   "lon_uy lat_uy" )
    call file_history_set_attribute( "grid_pressure", "edge2_coordinates",   "lon_xv lat_xv" )
    call file_history_set_attribute( "grid_pressure", "vertical_dimensions", "pressure" )

    call file_history_set_attribute( "grid_z", "cf_role",             "grid_topology", add_variable=.true. )
    call file_history_set_attribute( "grid_z", "topology_dimension",  (/ 2 /) )
    call file_history_set_attribute( "grid_z", "node_dimensions",     "xh yh" )
    call file_history_set_attribute( "grid_z", "face_dimensions",     "x: xh (padding: none) y: yh (padding: none)" )
    call file_history_set_attribute( "grid_z", "node_coordinates",    "lon_uv lat_uv" )
    call file_history_set_attribute( "grid_z", "face_coordinates",    "lon lat" )
    call file_history_set_attribute( "grid_z", "edge1_coordinates",   "lon_uy lat_uy" )
    call file_history_set_attribute( "grid_z", "edge2_coordinates",   "lon_xv lat_xv" )
    call file_history_set_attribute( "grid_z", "vertical_dimensions", "height_xyw: height (padding: none)" )

    call file_history_set_attribute( "grid_model", "cf_role",             "grid_topology", add_variable=.true. )
    call file_history_set_attribute( "grid_model", "topology_dimension",  (/ 2 /) )
    call file_history_set_attribute( "grid_model", "node_dimensions",     "FX FY" )
    call file_history_set_attribute( "grid_model", "face_dimensions",     "CX: FY (padding: none) CY: FY (padding: none)" )
    call file_history_set_attribute( "grid_model", "vertical_dimensions", "CZ: FZ (padding: none)" )

    call file_history_set_attribute( "grid_model_global", "cf_role",             "grid_topology", add_variable=.true. )
    call file_history_set_attribute( "grid_model_global", "topology_dimension",  (/ 2 /) )
    call file_history_set_attribute( "grid_model_global", "node_dimensions",     "FXG FYG" )
    call file_history_set_attribute( "grid_model_global", "face_dimensions",     "CXG: FYG (padding: none) CYG: FYG (padding: none)" )
    call file_history_set_attribute( "grid_model_global", "vertical_dimensions", "CZ: FZ (padding: none)" )

    return
  end subroutine file_history_cartesc_set_axes_attributes

end module scale_file_history_cartesc