domain_read.inc

!***********************************************************************
!*                             Apache License 2.0
!*
!* This file is part of the GFDL Flexible Modeling System (FMS).
!*
!* Licensed under the Apache License, Version 2.0 (the "License");
!* you may not use this file except in compliance with the License.
!* You may obtain a copy of the License at
!*
!*     http://www.apache.org/licenses/LICENSE-2.0
!*
!* FMS is distributed in the hope that it will be useful, but WITHOUT
!* WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied;
!* without even the implied warranty of MERCHANTABILITY or FITNESS FOR A
!* PARTICULAR PURPOSE. See the License for the specific language
!* governing permissions and limitations under the License.
!***********************************************************************
!> @file
!> @brief Reads domain decomposed variables and scatters data to pe's for the @ref read_data interface
 
!> @addtogroup fms2_io_mod
!> @{
 
!> @brief I/O domain root reads in  a domain decomposed variable at a
!!        specific unlimited dimension level and scatters the data to the
!!        rest of the ranks using its I/O compute domain indices. This
!!        routine may only be used with variables that are "domain
!!        decomposed".
subroutine domain_read_0d(fileobj, variable_name, vdata, unlim_dim_level, corner)
 
  type(fmsnetcdfdomainfile_t), intent(in) :: fileobj !< File object.
  character(len=*), intent(in) :: variable_name !< Variable name.
  class(*), intent(inout) :: vdata !< Data that will
                                   !! be read out
                                   !! to the netcdf file.
  integer, intent(in), optional :: unlim_dim_level !< Level for the unlimited
                                                   !! dimension.
  integer, intent(in), optional :: corner !< Array of starting
                                          !! indices describing
                                          !! where the data
                                          !! will be read to.
 
  call netcdf_read_data(fileobj, variable_name, vdata, &
                          unlim_dim_level=unlim_dim_level, corner=corner, &
                          broadcast=.true.)
 
end subroutine domain_read_0d
 
 
!> @brief I/O domain root reads in  a domain decomposed variable at a
!!        specific unlimited dimension level and scatters the data to the
!!        rest of the ranks using its I/O compute domain indices. This
!!        routine may only be used with variables that are "domain
!!        decomposed".
subroutine domain_read_1d(fileobj, variable_name, vdata, unlim_dim_level, &
                          corner, edge_lengths)
 
  type(fmsnetcdfdomainfile_t), intent(in) :: fileobj !< File object.
  character(len=*), intent(in) :: variable_name !< Variable name.
  class(*), dimension(:), intent(inout) :: vdata !< Data that will
                                                 !! be read out
                                                 !! to the netcdf file.
  integer, intent(in), optional :: unlim_dim_level !< Level for the unlimited
                                                   !! dimension.
  integer, dimension(1), intent(in), optional :: corner !< Array of starting
                                                        !! indices describing
                                                        !! where the data
                                                        !! will be read to.
  integer, dimension(1), intent(in), optional :: edge_lengths !< The number of
                                                              !! elements that
                                                              !! will be read
                                                              !! in each dimension.
 
  call netcdf_read_data(fileobj, variable_name, vdata, &
                          unlim_dim_level=unlim_dim_level, corner=corner, &
                          edge_lengths=edge_lengths, broadcast=.true.)
 
end subroutine domain_read_1d
 
 
!> @brief I/O domain root reads in  a domain decomposed variable at a
!!        specific unlimited dimension level and scatters the data to the
!!        rest of the ranks using its I/O compute domain indices. This
!!        routine may only be used with variables that are "domain
!!        decomposed".
subroutine domain_read_2d(fileobj, variable_name, vdata, unlim_dim_level, &
                          corner, edge_lengths)
  type(fmsnetcdfdomainfile_t),  intent(in)           :: fileobj         !< File object.
  character(len=*),             intent(in)           :: variable_name   !< Variable name.
  class(*), contiguous, target, intent(inout)        :: vdata(:,:)      !< Data that will
                                                                        !! be read out
                                                                        !! to the netcdf file.
  integer,                      intent(in), optional :: unlim_dim_level !< Level for the unlimited
                                                                        !! dimension.
  integer, dimension(2),        intent(in), optional :: corner          !< Array of starting
                                                                        !! indices describing
                                                                        !! where the data
                                                                        !! will be read to.
  integer, dimension(2),        intent(in), optional :: edge_lengths    !< The number of
                                                                        !! elements that
                                                                        !! will be read
                                                                        !! in each dimension.
 
  integer                 :: xdim_index            !< The index of the variable that is the x dimension
  integer                 :: ydim_index            !< The index of the variable that is the y dimension
  integer                 :: xpos                  !< The position of the x axis
  integer                 :: ypos                  !< The position of the y axis
  integer                 :: i                     !< For do loops
  integer                 :: isd                   !< The starting x position of the data io_domain
  integer                 :: isc                   !< The starting x position of the compute io_domain
  integer                 :: xc_size               !< The size of the x compute io_domain
  integer                 :: yc_size               !< The size of the y compute io_domain
  integer                 :: jsd                   !< The ending x position of the data io_domain
  integer                 :: jsc                   !< The ending y position of the compute io_domain
  integer                 :: c(2)                  !< The corners of the data
  integer                 :: e(2)                  !< The number of points (edges)
  logical                 :: buffer_includes_halos !< .True. if vdata includes halo points
  integer                 :: xgbegin               !< Starting x index of global io domain
  integer                 :: xgsize                !< Size of global x io domain
  integer                 :: ygbegin               !< Starting y index of global io domain
  integer                 :: ygsize                !< Size of global y io domain
  type(domain2d), pointer :: io_domain             !< pointer to the io_domain
 
  !< The global data is only allocated by the io root PEs
  integer(kind=i4_kind), dimension(:,:),     allocatable :: buf_i4_kind_pe !< PES section of the data
  integer(kind=i8_kind), dimension(:,:),     allocatable :: buf_i8_kind_pe !< PES section of the data
  real(kind=r4_kind),    dimension(:,:),     allocatable :: buf_r4_kind_pe !< PES section of the data
  real(kind=r8_kind),    dimension(:,:),     allocatable :: buf_r8_kind_pe !< PES section of the data
  integer(kind=i4_kind), dimension(:,:),     allocatable :: buf_i4_kind    !< Global section of the data
  integer(kind=i8_kind), dimension(:,:),     allocatable :: buf_i8_kind    !< Global section of the data
  real(kind=r4_kind),    dimension(:,:),     allocatable :: buf_r4_kind    !< Global section of the data
  real(kind=r8_kind),    dimension(:,:),     allocatable :: buf_r8_kind    !< Global section of the data
  class(*),              dimension(:,:,:,:), pointer     :: vdata_dummy    !< Vdata remapped as 4D
 
  if (.not. is_variable_domain_decomposed(fileobj, variable_name, .true., &
                                          xdim_index, ydim_index, xpos, ypos)) then
    call netcdf_read_data(fileobj, variable_name, vdata, &
                          unlim_dim_level=unlim_dim_level, corner=corner, &
                          edge_lengths=edge_lengths, broadcast=.true.)
    return
  endif
 
  call domain_offsets(size(vdata, xdim_index), size(vdata, ydim_index), fileobj%domain, &
                      xpos, ypos, isd, isc, xc_size, jsd, jsc, yc_size, buffer_includes_halos, &
                      msg="file:"//trim(fileobj%path)//" and variable:"//trim(variable_name))
 
  if (fileobj%use_netcdf_mpi) then
    c = 1
    e = shape(vdata)
 
    if (buffer_includes_halos) then
      c(xdim_index) = isd
      c(ydim_index) = jsd
      e(xdim_index) = xc_size + 2*(isc - isd)
      e(ydim_index) = yc_size + 2*(jsc - jsd)
    else
      c(xdim_index) = isc
      c(ydim_index) = jsc
      e(xdim_index) = xc_size
      e(ydim_index) = yc_size
    endif
 
    call netcdf_read_data(fileobj, variable_name, vdata, unlim_dim_level=unlim_dim_level, corner=c, edge_lengths=e)
    return
  endif
 
  if (xdim_index .ne. 1 .or. ydim_index .ne. 2) then
    ! This is a KLUDGE
    ! mpp_scatter assumes that the variable is (x,y), if that is not the case it remaps the data
    ! to a 4D array and calls domain_read_4d which does not use mpp_scatter yet
    vdata_dummy(1:size(vdata,1),1:size(vdata,2), 1:1, 1:1) => vdata(:,:)
    call domain_read_4d(fileobj, variable_name, vdata_dummy, unlim_dim_level)
    return
  endif
  io_domain => mpp_get_io_domain(fileobj%domain)
  c(:) = 1
  e(:) = shape(vdata)
 
  call mpp_get_global_domain(io_domain, xbegin=xgbegin, xsize=xgsize, position=xpos)
  call mpp_get_global_domain(io_domain, ybegin=ygbegin, ysize=ygsize, position=ypos)
 
  !I/O root reads in the data and scatters it.
  if (fileobj%is_root) then
 
    if (fileobj%adjust_indices) then
      !< If the file is distributed, the file only contains the io global domain
      c(xdim_index) = 1
      c(ydim_index) = 1
    else
      !< If the file is not distributed read, the file contains the global domain,
      !! so you only need to read the global io domain
      c(xdim_index) = xgbegin
      c(ydim_index) = ygbegin
    endif
 
    e(xdim_index) = xgsize
    e(ydim_index) = ygsize
 
    !Read in the global io domain
    select type(vdata)
      type is (integer(kind=i4_kind))
        call allocate_array(buf_i4_kind, e)
        call netcdf_read_data(fileobj, variable_name, buf_i4_kind, &
                              unlim_dim_level=unlim_dim_level, &
                              corner=c, edge_lengths=e, broadcast=.false.)
      type is (integer(kind=i8_kind))
        call allocate_array(buf_i8_kind, e)
        call netcdf_read_data(fileobj, variable_name, buf_i8_kind, &
                              unlim_dim_level=unlim_dim_level, &
                              corner=c, edge_lengths=e, broadcast=.false.)
      type is (real(kind=r4_kind))
        call allocate_array(buf_r4_kind, e)
        call netcdf_read_data(fileobj, variable_name, buf_r4_kind, &
                              unlim_dim_level=unlim_dim_level, &
                              corner=c, edge_lengths=e, broadcast=.false.)
      type is (real(kind=r8_kind))
        call allocate_array(buf_r8_kind, e)
        call netcdf_read_data(fileobj, variable_name, buf_r8_kind, &
                              unlim_dim_level=unlim_dim_level, &
                              corner=c, edge_lengths=e, broadcast=.false.)
      class default
        call error("unsupported variable type: domain_read_2d: file: "//trim(fileobj%path)//" variable:"// &
                  & trim(variable_name))
    end select
 
  endif
 
  c = 1
  e = shape(vdata)
 
  if (buffer_includes_halos) then
    !Adjust if the input buffer has room for halos.
    c(xdim_index) = isc - isd + 1
    c(ydim_index) = jsc - jsd + 1
  else
    c(xdim_index) = 1
    c(ydim_index) = 1
  endif
 
  e(xdim_index) = xc_size
  e(ydim_index) = yc_size
 
  select type(vdata)
  type is (integer(kind=i4_kind))
    call allocate_array(buf_i4_kind_pe, e)
    call mpp_scatter(isc-xgbegin+1, isc+xc_size-xgbegin, jsc-ygbegin+1, jsc+yc_size-ygbegin, fileobj%pelist, &
                     buf_i4_kind_pe, buf_i4_kind, fileobj%is_root)
    call put_array_section(buf_i4_kind_pe, vdata, c, e)
    deallocate(buf_i4_kind_pe)
  type is (integer(kind=i8_kind))
    call allocate_array(buf_i8_kind_pe, e)
    call mpp_scatter(isc-xgbegin+1, isc+xc_size-xgbegin, jsc-ygbegin+1, jsc+yc_size-ygbegin, fileobj%pelist, &
                     buf_i8_kind_pe, buf_i8_kind, fileobj%is_root)
    call put_array_section(buf_i8_kind_pe, vdata, c, e)
    deallocate(buf_i8_kind_pe)
  type is (real(kind=r4_kind))
    call allocate_array(buf_r4_kind_pe, e)
    call mpp_scatter(isc-xgbegin+1, isc+xc_size-xgbegin, jsc-ygbegin+1, jsc+yc_size-ygbegin, fileobj%pelist, &
                     buf_r4_kind_pe, buf_r4_kind, fileobj%is_root)
    call put_array_section(buf_r4_kind_pe, vdata, c, e)
    deallocate(buf_r4_kind_pe)
  type is (real(kind=r8_kind))
    call allocate_array(buf_r8_kind_pe, e)
    call mpp_scatter(isc-xgbegin+1, isc+xc_size-xgbegin, jsc-ygbegin+1, jsc+yc_size-ygbegin, fileobj%pelist, &
                     buf_r8_kind_pe, buf_r8_kind, fileobj%is_root)
    call put_array_section(buf_r8_kind_pe, vdata, c, e)
    deallocate(buf_r8_kind_pe)
  class default
    call error("unsupported variable type: domain_read_2d: file: "//trim(fileobj%path)//" variable:"// &
               & trim(variable_name))
  end select
 
  if (fileobj%is_root) then
    if (allocated(buf_i4_kind)) deallocate(buf_i4_kind)
    if (allocated(buf_i8_kind)) deallocate(buf_i8_kind)
    if (allocated(buf_r4_kind)) deallocate(buf_r4_kind)
    if (allocated(buf_r8_kind)) deallocate(buf_r8_kind)
  endif
end subroutine domain_read_2d
 
 
!> @brief I/O domain root reads in  a domain decomposed variable at a
!!        specific unlimited dimension level and scatters the data to the
!!        rest of the ranks using its I/O compute domain indices. This
!!        routine may only be used with variables that are "domain
!!        decomposed".
subroutine domain_read_3d(fileobj, variable_name, vdata, unlim_dim_level, &
                          corner, edge_lengths)
  type(fmsnetcdfdomainfile_t),  intent(in)           :: fileobj         !< File object.
  character(len=*),             intent(in)           :: variable_name   !< Variable name.
  class(*), contiguous, target, intent(inout)        :: vdata(:,:,:)    !< Data that will
                                                                        !! be read out
                                                                        !! to the netcdf file.
  integer,                      intent(in), optional :: unlim_dim_level !< Level for the unlimited
                                                                        !! dimension.
  integer, dimension(3),        intent(in), optional :: corner          !< Array of starting
                                                                        !! indices describing
                                                                        !! where the data
                                                                        !! will be read to.
  integer, dimension(3),        intent(in), optional :: edge_lengths    !< The number of
                                                                        !! elements that
                                                                        !! will be read
                                                                        !! in each dimension.
 
  integer                 :: xdim_index            !< The index of the variable that is the x dimension
  integer                 :: ydim_index            !< The index of the variable that is the y dimension
  integer                 :: xpos                  !< The position of the x axis
  integer                 :: ypos                  !< The position of the y axis
  integer                 :: i                     !< For do loops
  integer                 :: isd                   !< The starting x position of the data io_domain
  integer                 :: isc                   !< The starting x position of the compute io_domain
  integer                 :: xc_size               !< The size of the x compute io_domain
  integer                 :: yc_size               !< The size of the y compute io_domain
  integer                 :: jsd                   !< The ending x position of the data io_domain
  integer                 :: jsc                   !< The ending y position of the compute io_domain
  integer                 :: c(3)                  !< The corners of the data
  integer                 :: e(3)                  !< The number of points (edges)
  logical                 :: buffer_includes_halos !< .True. if vdata includes halo points
  integer                 :: xgbegin               !< Starting x index of global io domain
  integer                 :: xgsize                !< Size of global x io domain
  integer                 :: ygbegin               !< Starting y index of global io domain
  integer                 :: ygsize                !< Size of global y io domain
  type(domain2d), pointer :: io_domain             !< pointer to the io_domain
 
  !< The global data is only allocated by the io root PEs
  integer(kind=i4_kind), dimension(:,:,:),   allocatable :: buf_i4_kind_pe !< PES section of the data
  integer(kind=i8_kind), dimension(:,:,:),   allocatable :: buf_i8_kind_pe !< PES section of the data
  real(kind=r4_kind),    dimension(:,:,:),   allocatable :: buf_r4_kind_pe !< PES section of the data
  real(kind=r8_kind),    dimension(:,:,:),   allocatable :: buf_r8_kind_pe !< PES section of the data
  integer(kind=i4_kind), dimension(:,:,:),   allocatable :: buf_i4_kind    !< Global section of the data
  integer(kind=i8_kind), dimension(:,:,:),   allocatable :: buf_i8_kind    !< Global section of the data
  real(kind=r4_kind),    dimension(:,:,:),   allocatable :: buf_r4_kind    !< Global section of the data
  real(kind=r8_kind),    dimension(:,:,:),   allocatable :: buf_r8_kind    !< Global section of the data
  class(*),              dimension(:,:,:,:), pointer     :: vdata_dummy    !< Vdata remapped as 4D
 
  if (.not. is_variable_domain_decomposed(fileobj, variable_name, .true., &
                                          xdim_index, ydim_index, xpos, ypos)) then
    call netcdf_read_data(fileobj, variable_name, vdata, &
                          unlim_dim_level=unlim_dim_level, corner=corner, &
                          edge_lengths=edge_lengths, broadcast=.true.)
    return
  endif
 
  call domain_offsets(size(vdata, xdim_index), size(vdata, ydim_index), fileobj%domain, &
                      xpos, ypos, isd, isc, xc_size, jsd, jsc, yc_size, buffer_includes_halos, &
                      msg="file:"//trim(fileobj%path)//" and variable:"//trim(variable_name))
 
  if (fileobj%use_netcdf_mpi) then
    c = 1
    e = shape(vdata)
 
    if (buffer_includes_halos) then
      c(xdim_index) = isd
      c(ydim_index) = jsd
      e(xdim_index) = xc_size + 2*(isc - isd)
      e(ydim_index) = yc_size + 2*(jsc - jsd)
    else
      c(xdim_index) = isc
      c(ydim_index) = jsc
      e(xdim_index) = xc_size
      e(ydim_index) = yc_size
    endif
 
    call netcdf_read_data(fileobj, variable_name, vdata, unlim_dim_level=unlim_dim_level, corner=c, edge_lengths=e)
    return
  endif
 
  if (xdim_index .ne. 1 .or. ydim_index .ne. 2) then
    ! This is a KLUDGE
    ! mpp_scatter assumes that the variable is (x,y), if that is not the case it remaps the data
    ! to a 4D array and calls domain_read_4d which does not use mpp_scatter yet
    vdata_dummy(1:size(vdata,1),1:size(vdata,2), 1:size(vdata,3), 1:1) => vdata(:,:,:)
    call domain_read_4d(fileobj, variable_name, vdata_dummy, unlim_dim_level)
    return
  endif
  io_domain => mpp_get_io_domain(fileobj%domain)
  c(:) = 1
  if (present(corner)) c = corner
 
  e(:) = shape(vdata)
  if (present(edge_lengths)) e = edge_lengths
 
  call mpp_get_global_domain(io_domain, xbegin=xgbegin, xsize=xgsize, position=xpos)
  call mpp_get_global_domain(io_domain, ybegin=ygbegin, ysize=ygsize, position=ypos)
 
  !I/O root reads in the data and scatters it.
  if (fileobj%is_root) then
 
    if (fileobj%adjust_indices) then
      !< If the file is distributed, the file only contains the io global domain
      c(xdim_index) = 1
      c(ydim_index) = 1
    else
      !< If the file is not distributed read, the file contains the global domain,
      !! so you only need to read the global io domain
      c(xdim_index) = xgbegin
      c(ydim_index) = ygbegin
    endif
 
    e(xdim_index) = xgsize
    e(ydim_index) = ygsize
 
    !Read in the global io domain
    select type(vdata)
      type is (integer(kind=i4_kind))
        call allocate_array(buf_i4_kind, e)
        call netcdf_read_data(fileobj, variable_name, buf_i4_kind, &
                              unlim_dim_level=unlim_dim_level, &
                              corner=c, edge_lengths=e, broadcast=.false.)
      type is (integer(kind=i8_kind))
        call allocate_array(buf_i8_kind, e)
        call netcdf_read_data(fileobj, variable_name, buf_i8_kind, &
                              unlim_dim_level=unlim_dim_level, &
                              corner=c, edge_lengths=e, broadcast=.false.)
      type is (real(kind=r4_kind))
        call allocate_array(buf_r4_kind, e)
        call netcdf_read_data(fileobj, variable_name, buf_r4_kind, &
                              unlim_dim_level=unlim_dim_level, &
                              corner=c, edge_lengths=e, broadcast=.false.)
      type is (real(kind=r8_kind))
        call allocate_array(buf_r8_kind, e)
        call netcdf_read_data(fileobj, variable_name, buf_r8_kind, &
                              unlim_dim_level=unlim_dim_level, &
                              corner=c, edge_lengths=e, broadcast=.false.)
      class default
        call error("unsupported variable type: domain_read_2d: file: "//trim(fileobj%path)//" variable:"// &
                  & trim(variable_name))
    end select
 
  endif
 
  c = 1
  e = shape(vdata)
 
  if (buffer_includes_halos) then
    !Adjust if the input buffer has room for halos.
    c(xdim_index) = isc - isd + 1
    c(ydim_index) = jsc - jsd + 1
  else
    c(xdim_index) = 1
    c(ydim_index) = 1
  endif
 
  e(xdim_index) = xc_size
  e(ydim_index) = yc_size
 
  select type(vdata)
  type is (integer(kind=i4_kind))
    call allocate_array(buf_i4_kind_pe, e)
    call mpp_scatter(isc-xgbegin+1, isc+xc_size-xgbegin, jsc-ygbegin+1, jsc+yc_size-ygbegin, e(3), fileobj%pelist, &
                     buf_i4_kind_pe, buf_i4_kind, fileobj%is_root)
    call put_array_section(buf_i4_kind_pe, vdata, c, e)
    deallocate(buf_i4_kind_pe)
  type is (integer(kind=i8_kind))
    call allocate_array(buf_i8_kind_pe, e)
    call mpp_scatter(isc-xgbegin+1, isc+xc_size-xgbegin, jsc-ygbegin+1, jsc+yc_size-ygbegin, e(3), fileobj%pelist, &
                     buf_i8_kind_pe, buf_i8_kind, fileobj%is_root)
    call put_array_section(buf_i8_kind_pe, vdata, c, e)
    deallocate(buf_i8_kind_pe)
  type is (real(kind=r4_kind))
    call allocate_array(buf_r4_kind_pe, e)
    call mpp_scatter(isc-xgbegin+1, isc+xc_size-xgbegin, jsc-ygbegin+1, jsc+yc_size-ygbegin, e(3), fileobj%pelist, &
                     buf_r4_kind_pe, buf_r4_kind, fileobj%is_root)
    call put_array_section(buf_r4_kind_pe, vdata, c, e)
    deallocate(buf_r4_kind_pe)
  type is (real(kind=r8_kind))
    call allocate_array(buf_r8_kind_pe, e)
    call mpp_scatter(isc-xgbegin+1, isc+xc_size-xgbegin, jsc-ygbegin+1, jsc+yc_size-ygbegin, e(3), fileobj%pelist, &
                     buf_r8_kind_pe, buf_r8_kind, fileobj%is_root)
    call put_array_section(buf_r8_kind_pe, vdata, c, e)
    deallocate(buf_r8_kind_pe)
  class default
    call error("unsupported variable type: domain_read_2d: file: "//trim(fileobj%path)//" variable:"// &
               & trim(variable_name))
  end select
 
  if (fileobj%is_root) then
    if (allocated(buf_i4_kind)) deallocate(buf_i4_kind)
    if (allocated(buf_i8_kind)) deallocate(buf_i8_kind)
    if (allocated(buf_r4_kind)) deallocate(buf_r4_kind)
    if (allocated(buf_r8_kind)) deallocate(buf_r8_kind)
  endif
 
end subroutine domain_read_3d
 
 
!> @brief I/O domain root reads in  a domain decomposed variable at a
!!        specific unlimited dimension level and scatters the data to the
!!        rest of the ranks using its I/O compute domain indices. This
!!        routine may only be used with variables that are "domain
!!        decomposed".
subroutine domain_read_4d(fileobj, variable_name, vdata, unlim_dim_level, &
                          corner, edge_lengths)
 
  type(fmsnetcdfdomainfile_t), intent(in) :: fileobj !< File object.
  character(len=*), intent(in) :: variable_name !< Variable name.
  class(*), dimension(:,:,:,:), intent(inout) :: vdata !< Data that will
                                                       !! be read out
                                                       !! to the netcdf file.
  integer, intent(in), optional :: unlim_dim_level !< Level for the unlimited
                                                   !! dimension.
  integer, dimension(4), intent(in), optional :: corner !< Array of starting
                                                        !! indices describing
                                                        !! where the data
                                                        !! will be read to.
  integer, dimension(4), intent(in), optional :: edge_lengths !< The number of
                                                              !! elements that
                                                              !! will be read
                                                              !! in each dimension.
 
  integer :: xdim_index
  integer :: ydim_index
  type(domain2d), pointer :: io_domain
  integer :: xpos
  integer :: ypos
  integer :: i
  integer :: isd
  integer :: isc
  integer :: xc_size
  integer :: jsd
  integer :: jsc
  integer :: yc_size
  integer, dimension(:), allocatable :: pe_isc
  integer, dimension(:), allocatable :: pe_icsize
  integer, dimension(:), allocatable :: pe_jsc
  integer, dimension(:), allocatable :: pe_jcsize
  integer, dimension(4) :: c
  integer, dimension(4) :: e
  integer(kind=i4_kind), dimension(:,:,:,:), allocatable :: buf_i4_kind
  integer(kind=i8_kind), dimension(:,:,:,:), allocatable :: buf_i8_kind
  real(kind=r4_kind), dimension(:,:,:,:), allocatable :: buf_r4_kind
  real(kind=r8_kind), dimension(:,:,:,:), allocatable :: buf_r8_kind
  logical :: buffer_includes_halos
  integer :: xgmin !< Starting x index of global io domain
  integer :: ygmin !< Starting y index of global io domain
 
  if (.not. is_variable_domain_decomposed(fileobj, variable_name, .true., &
                                          xdim_index, ydim_index, xpos, ypos)) then
    call netcdf_read_data(fileobj, variable_name, vdata, &
                          unlim_dim_level=unlim_dim_level, corner=corner, &
                          edge_lengths=edge_lengths, broadcast=.true.)
    return
  endif
 
  call domain_offsets(size(vdata, xdim_index), size(vdata, ydim_index), fileobj%domain, &
                      xpos, ypos, isd, isc, xc_size, jsd, jsc, yc_size, buffer_includes_halos, &
                      msg="file:"//trim(fileobj%path)//" and variable:"//trim(variable_name))
 
  if (fileobj%use_netcdf_mpi) then
    c = 1
    e = shape(vdata)
 
    if (buffer_includes_halos) then
      c(xdim_index) = isd
      c(ydim_index) = jsd
      e(xdim_index) = xc_size + 2*(isc - isd)
      e(ydim_index) = yc_size + 2*(jsc - jsd)
    else
      c(xdim_index) = isc
      c(ydim_index) = jsc
      e(xdim_index) = xc_size
      e(ydim_index) = yc_size
    endif
 
    call netcdf_read_data(fileobj, variable_name, vdata, unlim_dim_level=unlim_dim_level, corner=c, edge_lengths=e)
    return
  endif
 
  io_domain => mpp_get_io_domain(fileobj%domain)
  c(:) = 1
  e(:) = shape(vdata)
  if (present(edge_lengths)) e = edge_lengths
 
  !I/O root reads in the data and scatters it.
  if (fileobj%is_root) then
    allocate(pe_isc(size(fileobj%pelist)))
    allocate(pe_icsize(size(fileobj%pelist)))
    allocate(pe_jsc(size(fileobj%pelist)))
    allocate(pe_jcsize(size(fileobj%pelist)))
    call mpp_get_compute_domains(io_domain, xbegin=pe_isc, xsize=pe_icsize, position=xpos)
    call mpp_get_compute_domains(io_domain, ybegin=pe_jsc, ysize=pe_jcsize, position=ypos)
    call mpp_get_global_domain(io_domain, xbegin=xgmin, position=xpos)
    call mpp_get_global_domain(io_domain, ybegin=ygmin, position=ypos)
    do i = 1, size(fileobj%pelist)
      if (present(corner)) c = corner
      c(xdim_index) = pe_isc(i)
      c(ydim_index) = pe_jsc(i)
      if (fileobj%adjust_indices) then
        c(xdim_index) = c(xdim_index) - xgmin + 1
        c(ydim_index) = c(ydim_index) - ygmin + 1
      endif
      e(xdim_index) = pe_icsize(i)
      e(ydim_index) = pe_jcsize(i)
      select type(vdata)
        type is (integer(kind=i4_kind))
          !Read in the data for fileobj%pelist(i)'s portion of the compute domain.
          call allocate_array(buf_i4_kind, e)
          call netcdf_read_data(fileobj, variable_name, buf_i4_kind, &
                                unlim_dim_level=unlim_dim_level, &
                                corner=c, edge_lengths=e, broadcast=.false.)
          if (i .eq. 1) then
            !Root rank stores data directly.
            c = 1
            if (buffer_includes_halos) then
              !Adjust if the input buffer has room for halos.
              c(xdim_index) = isc - isd + 1
              c(ydim_index) = jsc - jsd + 1
            endif
            call put_array_section(buf_i4_kind, vdata, c, e)
          else
            !Send data to non-root ranks.
            call mpp_send(buf_i4_kind, size(buf_i4_kind), fileobj%pelist(i))
            call mpp_sync_self(check=event_send)
          endif
          deallocate(buf_i4_kind)
        type is (integer(kind=i8_kind))
          !Read in the data for fileobj%pelist(i)'s portion of the compute domain.
          call allocate_array(buf_i8_kind, e)
          call netcdf_read_data(fileobj, variable_name, buf_i8_kind, &
                                unlim_dim_level=unlim_dim_level, &
                                corner=c, edge_lengths=e, broadcast=.false.)
          if (i .eq. 1) then
            !Root rank stores data directly.
            c = 1
            if (buffer_includes_halos) then
              !Adjust if the input buffer has room for halos.
              c(xdim_index) = isc - isd + 1
              c(ydim_index) = jsc - jsd + 1
            endif
            call put_array_section(buf_i8_kind, vdata, c, e)
          else
            !Send data to non-root ranks.
            call mpp_send(buf_i8_kind, size(buf_i8_kind), fileobj%pelist(i))
            call mpp_sync_self(check=event_send)
          endif
          deallocate(buf_i8_kind)
        type is (real(kind=r4_kind))
          !Read in the data for fileobj%pelist(i)'s portion of the compute domain.
          call allocate_array(buf_r4_kind, e)
          call netcdf_read_data(fileobj, variable_name, buf_r4_kind, &
                                unlim_dim_level=unlim_dim_level, &
                                corner=c, edge_lengths=e, broadcast=.false.)
          if (i .eq. 1) then
            !Root rank stores data directly.
            c = 1
            if (buffer_includes_halos) then
              !Adjust if the input buffer has room for halos.
              c(xdim_index) = isc - isd + 1
              c(ydim_index) = jsc - jsd + 1
            endif
            call put_array_section(buf_r4_kind, vdata, c, e)
          else
            !Send data to non-root ranks.
            call mpp_send(buf_r4_kind, size(buf_r4_kind), fileobj%pelist(i))
            call mpp_sync_self(check=event_send)
          endif
          deallocate(buf_r4_kind)
        type is (real(kind=r8_kind))
          !Read in the data for fileobj%pelist(i)'s portion of the compute domain.
          call allocate_array(buf_r8_kind, e)
          call netcdf_read_data(fileobj, variable_name, buf_r8_kind, &
                                unlim_dim_level=unlim_dim_level, &
                                corner=c, edge_lengths=e, broadcast=.false.)
          if (i .eq. 1) then
            !Root rank stores data directly.
            c = 1
            if (buffer_includes_halos) then
              !Adjust if the input buffer has room for halos.
              c(xdim_index) = isc - isd + 1
              c(ydim_index) = jsc - jsd + 1
            endif
            call put_array_section(buf_r8_kind, vdata, c, e)
          else
            !Send data to non-root ranks.
            call mpp_send(buf_r8_kind, size(buf_r8_kind), fileobj%pelist(i))
            call mpp_sync_self(check=event_send)
          endif
          deallocate(buf_r8_kind)
        class default
          call error("unsupported variable type: domain_read_4d: file: "//trim(fileobj%path)//" variable:"// &
                   & trim(variable_name))
      end select
    enddo
    deallocate(pe_isc)
    deallocate(pe_icsize)
    deallocate(pe_jsc)
    deallocate(pe_jcsize)
  else
    c = 1
    if (buffer_includes_halos) then
      c(xdim_index) = isc - isd + 1
      c(ydim_index) = jsc - jsd + 1
    endif
    e(xdim_index) = xc_size
    e(ydim_index) = yc_size
    select type(vdata)
      type is (integer(kind=i4_kind))
        call allocate_array(buf_i4_kind, e)
        call mpp_recv(buf_i4_kind, size(buf_i4_kind), fileobj%io_root, block=.true.)
        call put_array_section(buf_i4_kind, vdata, c, e)
        deallocate(buf_i4_kind)
      type is (integer(kind=i8_kind))
        call allocate_array(buf_i8_kind, e)
        call mpp_recv(buf_i8_kind, size(buf_i8_kind), fileobj%io_root, block=.true.)
        call put_array_section(buf_i8_kind, vdata, c, e)
        deallocate(buf_i8_kind)
      type is (real(kind=r4_kind))
        call allocate_array(buf_r4_kind, e)
        call mpp_recv(buf_r4_kind, size(buf_r4_kind), fileobj%io_root, block=.true.)
        call put_array_section(buf_r4_kind, vdata, c, e)
        deallocate(buf_r4_kind)
      type is (real(kind=r8_kind))
        call allocate_array(buf_r8_kind, e)
        call mpp_recv(buf_r8_kind, size(buf_r8_kind), fileobj%io_root, block=.true.)
        call put_array_section(buf_r8_kind, vdata, c, e)
        deallocate(buf_r8_kind)
      class default
        call error("unsupported variable type: domain_read_4d: file: "//trim(fileobj%path)//" variable:"// &
                 & trim(variable_name))
    end select
  endif
end subroutine domain_read_4d
 
 
!> @brief I/O domain root reads in  a domain decomposed variable at a
!!        specific unlimited dimension level and scatters the data to the
!!        rest of the ranks using its I/O compute domain indices. This
!!        routine may only be used with variables that are "domain
!!        decomposed".
subroutine domain_read_5d(fileobj, variable_name, vdata, unlim_dim_level, &
                          corner, edge_lengths)
 
  type(fmsnetcdfdomainfile_t), intent(in) :: fileobj !< File object.
  character(len=*), intent(in) :: variable_name !< Variable name.
  class(*), dimension(:,:,:,:,:), intent(inout) :: vdata !< Data that will
                                                         !! be read out
                                                         !! to the netcdf file.
  integer, intent(in), optional :: unlim_dim_level !< Level for the unlimited
                                                   !! dimension.
  integer, dimension(5), intent(in), optional :: corner !< Array of starting
                                                        !! indices describing
                                                        !! where the data
                                                        !! will be read to.
  integer, dimension(5), intent(in), optional :: edge_lengths !< The number of
                                                              !! elements that
                                                              !! will be read
                                                              !! in each dimension.
 
  integer :: xdim_index
  integer :: ydim_index
  type(domain2d), pointer :: io_domain
  integer :: xpos
  integer :: ypos
  integer :: i
  integer :: isd
  integer :: isc
  integer :: xc_size
  integer :: jsd
  integer :: jsc
  integer :: yc_size
  integer, dimension(:), allocatable :: pe_isc
  integer, dimension(:), allocatable :: pe_icsize
  integer, dimension(:), allocatable :: pe_jsc
  integer, dimension(:), allocatable :: pe_jcsize
  integer, dimension(5) :: c
  integer, dimension(5) :: e
  integer(kind=i4_kind), dimension(:,:,:,:,:), allocatable :: buf_i4_kind
  integer(kind=i8_kind), dimension(:,:,:,:,:), allocatable :: buf_i8_kind
  real(kind=r4_kind), dimension(:,:,:,:,:), allocatable :: buf_r4_kind
  real(kind=r8_kind), dimension(:,:,:,:,:), allocatable :: buf_r8_kind
  logical :: buffer_includes_halos
  integer :: xgmin !< Starting x index of global io domain
  integer :: ygmin !< Starting y index of global io domain
 
  if (.not. is_variable_domain_decomposed(fileobj, variable_name, .true., &
                                          xdim_index, ydim_index, xpos, ypos)) then
    call netcdf_read_data(fileobj, variable_name, vdata, &
                          unlim_dim_level=unlim_dim_level, corner=corner, &
                          edge_lengths=edge_lengths, broadcast=.true.)
    return
  endif
 
  call domain_offsets(size(vdata, xdim_index), size(vdata, ydim_index), fileobj%domain, &
                      xpos, ypos, isd, isc, xc_size, jsd, jsc, yc_size, buffer_includes_halos, &
                      msg="file:"//trim(fileobj%path)//" and variable:"//trim(variable_name))
 
  if (fileobj%use_netcdf_mpi) then
    c = 1
    e = shape(vdata)
 
    if (buffer_includes_halos) then
      c(xdim_index) = isd
      c(ydim_index) = jsd
      e(xdim_index) = xc_size + 2*(isc - isd)
      e(ydim_index) = yc_size + 2*(jsc - jsd)
    else
      c(xdim_index) = isc
      c(ydim_index) = jsc
      e(xdim_index) = xc_size
      e(ydim_index) = yc_size
    endif
 
    call netcdf_read_data(fileobj, variable_name, vdata, unlim_dim_level=unlim_dim_level, corner=c, edge_lengths=e)
    return
  endif
 
  io_domain => mpp_get_io_domain(fileobj%domain)
  c(:) = 1
  e(:) = shape(vdata)
  if (present(edge_lengths)) e = edge_lengths
 
  !I/O root reads in the data and scatters it.
  if (fileobj%is_root) then
    allocate(pe_isc(size(fileobj%pelist)))
    allocate(pe_icsize(size(fileobj%pelist)))
    allocate(pe_jsc(size(fileobj%pelist)))
    allocate(pe_jcsize(size(fileobj%pelist)))
    call mpp_get_compute_domains(io_domain, xbegin=pe_isc, xsize=pe_icsize, position=xpos)
    call mpp_get_compute_domains(io_domain, ybegin=pe_jsc, ysize=pe_jcsize, position=ypos)
    call mpp_get_global_domain(io_domain, xbegin=xgmin, position=xpos)
    call mpp_get_global_domain(io_domain, ybegin=ygmin, position=ypos)
    do i = 1, size(fileobj%pelist)
      !Calculate the indices of the domain-decomposed chunk relative to its position in the file.
      if (present(corner)) c = corner
      c(xdim_index) = pe_isc(i)
      c(ydim_index) = pe_jsc(i)
      if (fileobj%adjust_indices) then
        c(xdim_index) = c(xdim_index) - xgmin + 1
        c(ydim_index) = c(ydim_index) - ygmin + 1
      endif
      e(xdim_index) = pe_icsize(i)
      e(ydim_index) = pe_jcsize(i)
      select type(vdata)
        type is (integer(kind=i4_kind))
          !Read in the data for fileobj%pelist(i)'s portion of the compute domain.
          call allocate_array(buf_i4_kind, e)
          call netcdf_read_data(fileobj, variable_name, buf_i4_kind, &
                                unlim_dim_level=unlim_dim_level, &
                                corner=c, edge_lengths=e, broadcast=.false.)
          if (i .eq. 1) then
            !Root rank stores data directly.  Re-adjust the indicies relative
            !to the input buffer vdata.
            c = 1
            if (buffer_includes_halos) then
              !Adjust if the input buffer has room for halos.
              c(xdim_index) = isc - isd + 1
              c(ydim_index) = jsc - jsd + 1
            endif
            call put_array_section(buf_i4_kind, vdata, c, e)
          else
            !Send data to non-root ranks.
            call mpp_send(buf_i4_kind, size(buf_i4_kind), fileobj%pelist(i))
            call mpp_sync_self(check=event_send)
          endif
          deallocate(buf_i4_kind)
        type is (integer(kind=i8_kind))
          !Read in the data for fileobj%pelist(i)'s portion of the compute domain.
          call allocate_array(buf_i8_kind, e)
          call netcdf_read_data(fileobj, variable_name, buf_i8_kind, &
                                unlim_dim_level=unlim_dim_level, &
                                corner=c, edge_lengths=e, broadcast=.false.)
          if (i .eq. 1) then
            !Root rank stores data directly.
            c = 1
            if (buffer_includes_halos) then
              !Adjust if the input buffer has room for halos.
              c(xdim_index) = isc - isd + 1
              c(ydim_index) = jsc - jsd + 1
            endif
            call put_array_section(buf_i8_kind, vdata, c, e)
          else
            !Send data to non-root ranks.
            call mpp_send(buf_i8_kind, size(buf_i8_kind), fileobj%pelist(i))
            call mpp_sync_self(check=event_send)
          endif
          deallocate(buf_i8_kind)
        type is (real(kind=r4_kind))
          !Read in the data for fileobj%pelist(i)'s portion of the compute domain.
          call allocate_array(buf_r4_kind, e)
          call netcdf_read_data(fileobj, variable_name, buf_r4_kind, &
                                unlim_dim_level=unlim_dim_level, &
                                corner=c, edge_lengths=e, broadcast=.false.)
          if (i .eq. 1) then
            !Root rank stores data directly.
            c = 1
            if (buffer_includes_halos) then
              !Adjust if the input buffer has room for halos.
              c(xdim_index) = isc - isd + 1
              c(ydim_index) = jsc - jsd + 1
            endif
            call put_array_section(buf_r4_kind, vdata, c, e)
          else
            !Send data to non-root ranks.
            call mpp_send(buf_r4_kind, size(buf_r4_kind), fileobj%pelist(i))
            call mpp_sync_self(check=event_send)
          endif
          deallocate(buf_r4_kind)
        type is (real(kind=r8_kind))
          !Read in the data for fileobj%pelist(i)'s portion of the compute domain.
          call allocate_array(buf_r8_kind, e)
          call netcdf_read_data(fileobj, variable_name, buf_r8_kind, &
                                unlim_dim_level=unlim_dim_level, &
                                corner=c, edge_lengths=e, broadcast=.false.)
          if (i .eq. 1) then
            !Root rank stores data directly.
            c = 1
            if (buffer_includes_halos) then
              !Adjust if the input buffer has room for halos.
              c(xdim_index) = isc - isd + 1
              c(ydim_index) = jsc - jsd + 1
            endif
            call put_array_section(buf_r8_kind, vdata, c, e)
          else
            !Send data to non-root ranks.
            call mpp_send(buf_r8_kind, size(buf_r8_kind), fileobj%pelist(i))
            call mpp_sync_self(check=event_send)
          endif
          deallocate(buf_r8_kind)
        class default
          call error("unsupported variable type: domain_read_5d: file: "//trim(fileobj%path)//" variable:"// &
                   & trim(variable_name))
      end select
    enddo
    deallocate(pe_isc)
    deallocate(pe_icsize)
    deallocate(pe_jsc)
    deallocate(pe_jcsize)
  else
    c = 1
    if (buffer_includes_halos) then
      c(xdim_index) = isc - isd + 1
      c(ydim_index) = jsc - jsd + 1
    endif
    e(xdim_index) = xc_size
    e(ydim_index) = yc_size
    select type(vdata)
      type is (integer(kind=i4_kind))
        call allocate_array(buf_i4_kind, e)
        call mpp_recv(buf_i4_kind, size(buf_i4_kind), fileobj%io_root, block=.true.)
        call put_array_section(buf_i4_kind, vdata, c, e)
        deallocate(buf_i4_kind)
      type is (integer(kind=i8_kind))
        call allocate_array(buf_i8_kind, e)
        call mpp_recv(buf_i8_kind, size(buf_i8_kind), fileobj%io_root, block=.true.)
        call put_array_section(buf_i8_kind, vdata, c, e)
        deallocate(buf_i8_kind)
      type is (real(kind=r4_kind))
        call allocate_array(buf_r4_kind, e)
        call mpp_recv(buf_r4_kind, size(buf_r4_kind), fileobj%io_root, block=.true.)
        call put_array_section(buf_r4_kind, vdata, c, e)
        deallocate(buf_r4_kind)
      type is (real(kind=r8_kind))
        call allocate_array(buf_r8_kind, e)
        call mpp_recv(buf_r8_kind, size(buf_r8_kind), fileobj%io_root, block=.true.)
        call put_array_section(buf_r8_kind, vdata, c, e)
        deallocate(buf_r8_kind)
      class default
        call error("unsupported variable type: domain_read_5d: file: "//trim(fileobj%path)//" variable:"// &
                 & trim(variable_name))
    end select
  endif
end subroutine domain_read_5d
!> @}