fms_netcdf_domain_io.F90

!***********************************************************************
!*                             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.
!***********************************************************************
!> @defgroup fms_netcdf_domain_io_mod fms_netcdf_domain_io_mod
!> @ingroup fms2_io
!> @brief Domain-specific I/O wrappers.
 
!> @addtogroup fms_netcdf_domain_io_mod
!> @{
module fms_netcdf_domain_io_mod
use netcdf
use mpp_mod
use mpp_domains_mod
use fms_io_utils_mod
use netcdf_io_mod
use platform_mod
implicit none
private
 
 
!Module constants.
integer, parameter :: no_domain_decomposed_dimension = 0
integer, parameter, public :: max_num_domain_decomposed_dims = 10
integer, parameter :: variable_not_found = 0
integer, parameter :: default_domain_position = center
character(len=16), parameter :: domain_pos_att = "domain_position"
character(len=16), parameter :: domain_axis_att_name = "domain_axis"
character(len=16), parameter :: x = "x"
character(len=16), parameter :: y = "y"
 
!> @}
 
!> @brief Domain variable.
!> @ingroup fms_netcdf_domain_io_mod
type, private :: domaindimension_t
  character(len=nf90_max_name) :: varname !< Variable name.
  integer :: pos !< Domain position.
endtype domaindimension_t
 
 
!> @brief netcdf domain file type.
!> @ingroup fms_netcdf_domain_io_mod
type, extends(fmsnetcdffile_t), public :: fmsnetcdfdomainfile_t
  type(domain2d) :: domain !< Two-dimensional domain.
  type(domaindimension_t), dimension(:), allocatable :: xdims !< Dimensions associated
                                                              !! with the "x" axis
                                                              !! of a 2d domain.
  integer :: nx !< Number of "x" dimensions.
  type(domaindimension_t), dimension(:), allocatable :: ydims !< Dimensions associated
                                                              !! with the "y" axis
                                                              !! of a 2d domain.
  integer :: ny !< Number of "y" dimensions.
  character(len=FMS_PATH_LEN) :: non_mangled_path !< Non-domain-mangled file path.
  logical :: adjust_indices !< Flag telling if indices need to be adjusted
                            !! for domain-decomposed read.
endtype fmsnetcdfdomainfile_t
 
 
public :: open_domain_file
public :: close_domain_file
public :: register_domain_decomposed_dimension
public :: register_domain_variable
public :: register_domain_restart_variable_0d
public :: register_domain_restart_variable_1d
public :: register_domain_restart_variable_2d
public :: register_domain_restart_variable_3d
public :: register_domain_restart_variable_4d
public :: register_domain_restart_variable_5d
public :: domain_read_0d
public :: domain_read_1d
public :: domain_read_2d
public :: domain_read_3d
public :: domain_read_4d
public :: domain_read_5d
public :: domain_write_0d
public :: domain_write_1d
public :: domain_write_2d
public :: domain_write_3d
public :: domain_write_4d
public :: domain_write_5d
public :: save_domain_restart
public :: restore_domain_state
public :: get_compute_domain_dimension_indices
public :: get_global_io_domain_indices
public :: is_dimension_registered
public :: get_mosaic_tile_grid
 
!> @ingroup fms_netcdf_domain_io_mod
interface compute_global_checksum
  module procedure compute_global_checksum_2d
  module procedure compute_global_checksum_3d
  module procedure compute_global_checksum_4d
end interface compute_global_checksum
 
!> @addtogroup fms_netcdf_domain_io_mod
!> @{
 
contains
 
 
!> @brief Get the index of a domain decomposed dimension.
!! @return Index of domain decomposed dimension.
function get_domain_decomposed_index(name_, array, size_) &
  result(index_)
 
  character(len=*), intent(in) :: name_ !< Name.
  type(domaindimension_t), dimension(:), intent(in) :: array !< Array to search through.
  integer, intent(in) :: size_ !< Number of spots to look in.
  integer :: index_
 
  integer :: i
 
  index_ = variable_not_found
  do i = 1, size_
    if (string_compare(array(i)%varname, name_)) then
      index_ = i
      return
    endif
  enddo
end function get_domain_decomposed_index
 
 
!> @brief Add a domain decomposed dimension to an array.
subroutine append_domain_decomposed_dimension(name_, position_, array, size_)
 
  character(len=*), intent(in) :: name_ !< Variable name.
  integer, intent(in) :: position_ !< Domain position.
  type(domaindimension_t), dimension(:), intent(inout) :: array !< Array to search through.
  integer, intent(inout) :: size_ !< Number of spots to look in.
 
  integer :: i
 
  do i = 1, size_
    if (string_compare(array(i)%varname, name_)) then
      call error("variable "//trim(name_)//" already registered.")
    endif
  enddo
  size_ = size_ + 1
  if (size_ .gt. size(array)) then
    call error("number of domain decomposed variables exceeds limit.")
  endif
  call string_copy(array(size_)%varname, name_)
  array(size_)%pos = position_
end subroutine append_domain_decomposed_dimension
 
 
!> @brief Given a domain decomposed dimension, get its domain position.
!! @return Position of the domain decomposed variable.
function get_domain_position(name_, array, size_) &
  result(dpos)
 
  character(len=*), intent(in) :: name_ !< Variable name.
  type(domaindimension_t), dimension(:), intent(in) :: array !< Array to search through.
  integer, intent(in) :: size_
  integer :: dpos
 
  dpos = get_domain_decomposed_index(name_, array, size_)
  if (dpos .ne. variable_not_found) then
    dpos = array(dpos)%pos
  endif
end function get_domain_position
 
 
!> @brief Given a variable, get the index of the "x" or "y" domain decomposed
!!        dimension.
!! @return Index of the domain decomposed dimension or else
!!         no_domain_decomposed_dimension.
function get_domain_decomposed_dimension_index(fileobj, variable_name, &
                                               xory, broadcast) &
  result(index_)
 
  type(fmsnetcdfdomainfile_t), intent(in), target :: fileobj !< File object.
  character(len=*), intent(in) :: variable_name !< Variable name.
  character(len=*), intent(in) :: xory !< String telling which dimension to
                                       !! look for.  Valid values are "x"
                                       !! or "y".
  logical, intent(in), optional :: broadcast !< Flag controlling whether or
                                             !! not the index will be
                                             !! broadcast to non "I/O root"
                                             !! ranks.  The broadcast will
                                             !! be done by default.
  integer :: index_
 
  integer :: ndims
  character(len=nf90_max_name), dimension(:), allocatable :: dim_names
  type(domaindimension_t), dimension(:), pointer :: p
  integer :: n
  integer :: i
 
  index_ = no_domain_decomposed_dimension
  if (fileobj%is_root) then
    ndims = get_variable_num_dimensions(fileobj, variable_name, broadcast=.false.)
    allocate(dim_names(ndims))
    dim_names(:) = ""
    call get_variable_dimension_names(fileobj, variable_name, dim_names, broadcast=.false.)
    if (string_compare(xory, x, .true.)) then
      p => fileobj%xdims
      n = fileobj%nx
    elseif (string_compare(xory, y, .true.)) then
      p => fileobj%ydims
      n = fileobj%ny
    else
      call error("unrecognized xory flag value.")
    endif
    do i = 1, size(dim_names)
      if (get_domain_decomposed_index(dim_names(i), p, n) .ne. variable_not_found) then
        index_ = i
        exit
      endif
    enddo
    deallocate(dim_names)
  endif
  if (present(broadcast)) then
    if (.not. broadcast) then
      return
    endif
  endif
  call mpp_broadcast(index_, fileobj%io_root, pelist=fileobj%pelist)
end function get_domain_decomposed_dimension_index
 
 
!> @brief Determine if a variable is "domain decomposed."
!! @return Flag telling if the variable is "domain decomposed."
function is_variable_domain_decomposed(fileobj, variable_name, broadcast, &
                                       xindex, yindex, xpos, ypos) &
  result(is_decomposed)
 
  type(fmsnetcdfdomainfile_t), intent(in) :: fileobj !< File object.
  character(len=*), intent(in) :: variable_name !< Variable name.
  logical, intent(in), optional :: broadcast !< Flag controlling whether or
                                             !! not the index will be
                                             !! broadcast to non "I/O root"
                                             !! ranks.  The broadcast will
                                             !! be done by default.
  integer, intent(out), optional :: xindex !< The index of the domain
                                           !! x dimension.
  integer, intent(out), optional :: yindex !< The index of the domain
                                           !! y dimension.
  integer, intent(out), optional :: xpos !< Domain position of the x dimension.
  integer, intent(out), optional :: ypos !< Domain position of the y dimension.
  logical :: is_decomposed
 
  integer, dimension(2) :: indices
  integer :: ndims
  character(len=nf90_max_name), dimension(:), allocatable :: dim_names
 
  indices(1) = get_domain_decomposed_dimension_index(fileobj, variable_name, x, broadcast)
  if (present(xindex)) then
    xindex = indices(1)
  endif
  indices(2) = get_domain_decomposed_dimension_index(fileobj, variable_name, y, broadcast)
  if (present(yindex)) then
    yindex = indices(2)
  endif
  is_decomposed = (indices(1) .ne. no_domain_decomposed_dimension) .and. &
                  (indices(2) .ne. no_domain_decomposed_dimension)
  if (is_decomposed) then
    if (.not. present(xpos) .and. .not. present(ypos)) then
      return
    endif
    ndims = get_variable_num_dimensions(fileobj, variable_name, broadcast)
    allocate(dim_names(ndims))
    dim_names(:) = ""
    call get_variable_dimension_names(fileobj, variable_name, dim_names, broadcast)
    if (present(xpos)) then
      xpos = get_domain_position(dim_names(indices(1)), fileobj%xdims, fileobj%nx)
    endif
    if (present(ypos)) then
      ypos = get_domain_position(dim_names(indices(2)), fileobj%ydims, fileobj%ny)
    endif
    deallocate(dim_names)
  else
    if (present(xpos)) then
      xpos = -1
    endif
    if (present(ypos)) then
      ypos = -1
    endif
  endif
end function is_variable_domain_decomposed
 
 
!> @brief Determine whether a domain-decomposed dimension has been registered to the file object
!! @return Flag telling if the dimension is registered to the file object
function is_dimension_registered(fileobj, dimension_name) &
  result(is_registered)
 
  type(fmsnetcdfdomainfile_t), intent(in) :: fileobj !< File object.
  character(len=*), intent(in) :: dimension_name !< Dimension name.
 
  ! local
  logical :: is_registered
  integer :: dpos
 
  dpos = 0
  is_registered = .false.
  dpos = get_domain_decomposed_index(dimension_name, fileobj%xdims, fileobj%nx)
  if (dpos .ne. variable_not_found) then
    is_registered = .true.
  else
    dpos = get_domain_decomposed_index(dimension_name, fileobj%ydims, fileobj%ny)
    if (dpos .ne. variable_not_found) is_registered = .true.
  endif
 
end function is_dimension_registered
 
!> @brief Open a NetCDF-4 file in parallel write mode
!! @return True on success
function open_collective_netcdf_file(fileobj, path, mode, domain, is_restart, dont_add_res_to_filename) &
  result(success)
 
  type(fmsnetcdfdomainfile_t),intent(inout) :: fileobj !< File object.
  character(len=*), intent(in) :: path !< File path.
  character(len=*), intent(in) :: mode !< File mode.  Allowed values
                                       !! are "read", "append", "write", or
                                       !! "overwrite".
  type(domain2d), intent(in) :: domain !< Two-dimensional domain.
  logical, intent(in), optional :: is_restart !< Flag telling if this file
                                              !! is a restart file.  Defaults
                                              !! to false.
  logical, intent(in), optional :: dont_add_res_to_filename !< Flag indicating not to add
                                              !! ".res" to the filename
 
  integer :: nc_format_param
  integer :: tile_id(1)
  integer :: err
  character(len=FMS_PATH_LEN) :: combined_filepath
  character(len=FMS_PATH_LEN) :: full_path
  logical :: is_res
  logical :: dont_add_res
  integer :: success
  integer :: domain_size
 
  success = .true.
  call string_copy(fileobj%non_mangled_path, path)
 
  !TODO Lots of duplicate code between this and netcdf_file_open
 
  !! Determine the name of your file !!
 
  !! If the number of tiles is greater than 1 or if the current tile is greater
  !than 1 add .tileX. to the filename
  tile_id = mpp_get_tile_id(domain)
  if (mpp_get_ntile_count(domain) .gt. 1 .or. tile_id(1) > 1) then
    call domain_tile_filepath_mangle(combined_filepath, path, tile_id(1))
  else
    call string_copy(combined_filepath, path)
  endif
 
  !< Only add ".res" to the file path if is_restart is set to true
  !! and dont_add_res_to_filename is set to false.
  is_res = .false.
  if (present(is_restart)) then
    is_res = is_restart
  endif
  fileobj%is_restart = is_res
 
  dont_add_res = .false.
  if (present(dont_add_res_to_filename)) then
    dont_add_res = dont_add_res_to_filename
  endif
 
  if (is_res .and. .not. dont_add_res) then
    call restart_filepath_mangle(full_path, trim(combined_filepath))
  else
    call string_copy(full_path, trim(combined_filepath))
  endif
 
  call string_copy(fileobj%path, trim(full_path))
 
  nc_format_param = ior(nf90_netcdf4, nf90_mpiio)
  fileobj%is_netcdf4 = .true.
 
  fileobj%domain = domain
  call mpp_get_global_domain(fileobj%domain, xsize=domain_size)
 
  if (string_compare(mode, "read", .true.) .or. string_compare(mode, "append", .true.)) &
    call error("The use_netcdf_mpi = .true. option for reads is currently not supported")
 
  if (string_compare(mode, "write", .true.)) then
    err = nf90_create(trim(fileobj%path), ior(nf90_noclobber, nc_format_param), fileobj%ncid, &
                      comm = mpp_get_domain_tile_commid(fileobj%domain), info = mpp_info_null)
  elseif (string_compare(mode,"overwrite",.true.)) then
    err = nf90_create(trim(fileobj%path), ior(nf90_clobber, nc_format_param), fileobj%ncid, &
                      comm = mpp_get_domain_tile_commid(fileobj%domain), info = mpp_info_null)
  endif
  call check_netcdf_code(err, "open_collective_netcdf_file:"//trim(fileobj%path))
 
  allocate(fileobj%xdims(max_num_domain_decomposed_dims))
  fileobj%nx = 0
  allocate(fileobj%ydims(max_num_domain_decomposed_dims))
  fileobj%ny = 0
 
  ! Every rank is the root PE of its own pelist. This forces all ranks to hit any NetCDF calls,
  ! which are usually inside `if (fileobj%is_root)` blocks.
  allocate(fileobj%pelist(1))
  fileobj%pelist(1) = mpp_pe()
  fileobj%io_root = mpp_pe()
  fileobj%is_root = .true.
 
  fileobj%use_collective = .false. !TODO
  fileobj%is_diskless = .false.
 
  if (fileobj%is_restart) then
    allocate(fileobj%restart_vars(max_num_restart_vars))
    fileobj%num_restart_vars = 0
  endif
 
  fileobj%is_readonly = string_compare(mode, "read", .true.)
  fileobj%mode_is_append = string_compare(mode, "append", .true.)
  allocate(fileobj%compressed_dims(max_num_compressed_dims))
  fileobj%num_compressed_dims = 0
  ! Set the is_open flag to true for this file object.
  if (.not.allocated(fileobj%is_open)) allocate(fileobj%is_open)
  fileobj%is_open = .true.
 
  fileobj%bc_dimensions%xlen = 0
  fileobj%bc_dimensions%ylen = 0
  fileobj%bc_dimensions%zlen = 0
  fileobj%bc_dimensions%cur_dim_len = 0
 
end function open_collective_netcdf_file
 
!> @brief Open a domain netcdf file.
!! @return Flag telling if the open completed successfully.
function open_domain_file(fileobj, path, mode, domain, nc_format, is_restart, dont_add_res_to_filename, &
                          use_netcdf_mpi) result(success)
 
  type(fmsnetcdfdomainfile_t),intent(inout) :: fileobj !< File object.
  character(len=*), intent(in) :: path !< File path.
  character(len=*), intent(in) :: mode !< File mode.  Allowed values
                                       !! are "read", "append", "write", or
                                       !! "overwrite".
  type(domain2d), intent(in) :: domain !< Two-dimensional domain.
  character(len=*), intent(in), optional :: nc_format !< Netcdf format that
                                                      !! new files are written
                                                      !! as.  Allowed values
                                                      !! are: "64bit", "classic",
                                                      !! or "netcdf4". Defaults to
                                                      !! "64bit".
  logical, intent(in), optional :: is_restart !< Flag telling if this file
                                              !! is a restart file.  Defaults
                                              !! to false.
  logical, intent(in), optional :: dont_add_res_to_filename !< Flag indicating not to add
                                              !! ".res" to the filename
  logical, intent(in), optional :: use_netcdf_mpi !< Flag telling if this file should be using netcdf4 collective
                                                  !! reads and writes. Defaults to false.
                                                  !! nc_format is automatically set to netcdf4
  logical :: success
 
  integer, dimension(2) :: io_layout
  integer, dimension(1) :: tile_id
  character(len=FMS_PATH_LEN) :: combined_filepath
  type(domain2d), pointer :: io_domain
  character(len=FMS_PATH_LEN) :: distributed_filepath
  integer :: pelist_size
  integer, dimension(:), allocatable :: pelist
  logical :: success2
  type(fmsnetcdfdomainfile_t) :: fileobj2
 
  io_domain => mpp_get_io_domain(domain)
 
  fileobj%use_netcdf_mpi = .false.
  if (present(use_netcdf_mpi)) fileobj%use_netcdf_mpi = use_netcdf_mpi
 
  if (fileobj%use_netcdf_mpi) then
#ifdef NO_NC_PARALLEL4
    call mpp_error(fatal, "NetCDF was not built with HDF5 parallel I/O features, so parallel writes are not supported. &
                          &Please turn parallel writes off for the file: " // trim(path))
#endif
 
    if (associated(io_domain)) then
      call mpp_error(note, "NetCDF MPI is enabled: ignoring I/O domain. Only one output file will be produced.")
    endif
 
    success = open_collective_netcdf_file(fileobj, path, mode, domain, is_restart, dont_add_res_to_filename)
    return
  endif
 
  !Get the path of a "combined" file.
  io_layout = mpp_get_io_domain_layout(domain)
  tile_id = mpp_get_tile_id(domain)
 
  !< If the number of tiles is greater than 1 or if the current tile is greater
  !than 1 add .tileX. to the filename
  if (mpp_get_ntile_count(domain) .gt. 1 .or. tile_id(1) > 1) then
    call domain_tile_filepath_mangle(combined_filepath, path, tile_id(1))
  else
    call string_copy(combined_filepath, path)
  endif
 
  !Get the path of a "distributed" file.
  if (.not. associated(io_domain)) then
    call error("The domain associated with the file:"//trim(path)//" does not have an io_domain.")
  endif
  if (io_layout(1)*io_layout(2) .gt. 1) then
    tile_id = mpp_get_tile_id(io_domain)
    call io_domain_tile_filepath_mangle(distributed_filepath, combined_filepath, tile_id(1))
  else
    call string_copy(distributed_filepath, combined_filepath)
  endif
 
  !Make sure the input domain has an I/O domain and get its pelist.
  pelist_size = mpp_get_domain_npes(io_domain)
  allocate(pelist(pelist_size))
  call mpp_get_pelist(io_domain, pelist)
  fileobj%adjust_indices = .true. !Set the default to true
 
  !Open the distibuted files.
  success = netcdf_file_open(fileobj, distributed_filepath, mode, nc_format, pelist, &
                             is_restart, dont_add_res_to_filename)
  if (string_compare(mode, "read", .true.) .or. string_compare(mode, "append", .true.)) then
    if (success) then
      if (.not. string_compare(distributed_filepath, combined_filepath)) then
        success2 = netcdf_file_open(fileobj2, combined_filepath, mode, nc_format, pelist, &
                                    is_restart, dont_add_res_to_filename)
        if (success2) then
          call error("The domain decomposed file:"//trim(path)// &
                   & " contains both combined (*.nc) and distributed files (*.nc.XXXX).")
        endif
      endif
    else
      success = netcdf_file_open(fileobj, combined_filepath, mode, nc_format, pelist, &
                                 is_restart, dont_add_res_to_filename)
      !If the file is combined and the layout is not (1,1) set the adjust_indices flag to false
      if (success .and. (io_layout(1)*io_layout(2) .gt. 1)) fileobj%adjust_indices = .false.
    endif
  endif
  if (.not. success) then
    deallocate(pelist)
    return
  endif
 
  !Store/initialize necessary properties.
  call string_copy(fileobj%non_mangled_path, path)
  fileobj%domain = domain
  allocate(fileobj%xdims(max_num_domain_decomposed_dims))
  fileobj%nx = 0
  allocate(fileobj%ydims(max_num_domain_decomposed_dims))
  fileobj%ny = 0
  call string_copy(fileobj%non_mangled_path, path)
 
  if (string_compare(mode, "write", .true.) .or. string_compare(mode, "overwrite", .true.)) then
    !Add global attribute needed by mppnccombine.
    call register_global_attribute(fileobj, "NumFilesInSet", io_layout(1)*io_layout(2))
  endif
end function open_domain_file
 
 
!> @brief Close a domain netcdf file.
subroutine close_domain_file(fileobj)
 
  type(fmsnetcdfdomainfile_t), intent(inout) :: fileobj !< File object.
 
  call netcdf_file_close(fileobj)
  deallocate(fileobj%xdims)
  fileobj%nx = 0
  deallocate(fileobj%ydims)
  fileobj%ny = 0
end subroutine close_domain_file
 
 
!> @brief Add a dimension to a file associated with a two-dimensional domain.
subroutine register_domain_decomposed_dimension(fileobj, dim_name, xory, domain_position)
 
  type(fmsnetcdfdomainfile_t), target, intent(inout) :: fileobj !< File object.
  character(len=*), intent(in) :: dim_name !< Dimension name.
  character(len=*), intent(in) :: xory !< Flag telling if the dimension
                                       !! is associated with the "x" or "y"
                                       !! axis of the 2d domain.  Allowed
                                       !! values are "x" or "y".
  integer, intent(in), optional :: domain_position !< Domain position.
 
  integer :: dpos
  type(domain2d), pointer :: io_domain
  integer :: domain_size
  integer :: dim_size
 
  dpos = default_domain_position
  if (mpp_domain_is_symmetry(fileobj%domain) .and. present(domain_position)) then
    dpos = domain_position
  endif
 
  ! If using NetCDF MPI, the IO domain is ignored, so use the domain to determine the correct size of each
  ! domain-decomposed dimension.
  if (fileobj%use_netcdf_mpi) then
    io_domain => fileobj%domain
  else
    io_domain => mpp_get_io_domain(fileobj%domain)
  endif
 
  if (string_compare(xory, x, .true.)) then
    if (dpos .ne. center .and. dpos .ne. east) then
      call error("Only domain_position=center or domain_position=EAST is supported for x dimensions."// &
                  & " Fix your register_axis call for file:"&
                  &//trim(fileobj%path)//" and dimension:"//trim(dim_name))
    endif
    call mpp_get_global_domain(io_domain, xsize=domain_size, position=dpos)
    call append_domain_decomposed_dimension(dim_name, dpos, fileobj%xdims, fileobj%nx)
  elseif (string_compare(xory, y, .true.)) then
    if (dpos .ne. center .and. dpos .ne. north) then
      call error("Only domain_position=center or domain_position=NORTH is supported for y dimensions."// &
                  & " Fix your register_axis call for file:"&
                  &//trim(fileobj%path)//" and dimension:"//trim(dim_name))
    endif
    call mpp_get_global_domain(io_domain, ysize=domain_size, position=dpos)
    call append_domain_decomposed_dimension(dim_name, dpos, fileobj%ydims, fileobj%ny)
  else
    call error("The register_axis call for file:"//trim(fileobj%path)//" and dimension:"//trim(dim_name)// &
               & " has an unrecognized xory flag value:"&
               &//trim(xory)//" only 'x' and 'y' are allowed.")
  endif
  if (fileobj%is_readonly .or. (fileobj%mode_is_append .and. dimension_exists(fileobj, dim_name))) then
    call get_dimension_size(fileobj, dim_name, dim_size, broadcast=.true.)
    if (dim_size .lt. domain_size) then
      call error("dimension "//trim(dim_name)//" in the file "//trim(fileobj%path)//" is smaller than the size of" &
                 //" the associated domain "//trim(xory)//" axis.")
    endif
  else
    call netcdf_add_dimension(fileobj, dim_name, domain_size, is_compressed=.false.)
  endif
end subroutine register_domain_decomposed_dimension
 
 
!> @brief Add a "domain_decomposed" attribute to the axis variables because it is
!!        required by mppnccombine.
!! @internal
subroutine add_domain_attribute(fileobj, variable_name)
 
  type(fmsnetcdfdomainfile_t), intent(inout) :: fileobj !< File object.
  character(len=*), intent(in) :: variable_name !< Variable name.
 
  type(domain2d), pointer :: io_domain
  integer :: dpos
  integer :: sg
  integer :: eg
  integer :: s
  integer :: e
  integer, dimension(2) :: io_layout !< Io_layout in the fileobj's domain
 
  !< Don't add the "domain_decomposition" variable attribute if the io_layout is
  !! 1,1, or if using mpi netcdf for writes, to avoid frecheck "failures"
  io_layout = mpp_get_io_domain_layout(fileobj%domain)
  if (io_layout(1)  .eq. 1 .and. io_layout(2) .eq. 1) return
  if (fileobj%use_netcdf_mpi) return
 
  io_domain => mpp_get_io_domain(fileobj%domain)
  dpos = get_domain_decomposed_index(variable_name, fileobj%xdims, fileobj%nx)
  if (dpos .ne. variable_not_found) then
    dpos = fileobj%xdims(dpos)%pos
    call mpp_get_global_domain(fileobj%domain, xbegin=sg, xend=eg, position=dpos)
    call mpp_get_global_domain(io_domain, xbegin=s, xend=e, position=dpos)
    call register_variable_attribute(fileobj, variable_name, "domain_decomposition", &
                                     (/sg, eg, s, e/))
  else
    dpos = get_domain_decomposed_index(variable_name, fileobj%ydims, fileobj%ny)
    if (dpos .ne. variable_not_found) then
      dpos = fileobj%ydims(dpos)%pos
      call mpp_get_global_domain(fileobj%domain, ybegin=sg, yend=eg, position=dpos)
      call mpp_get_global_domain(io_domain, ybegin=s, yend=e, position=dpos)
      call register_variable_attribute(fileobj, variable_name, "domain_decomposition", &
                                       (/sg, eg, s, e/))
    endif
  endif
end subroutine add_domain_attribute
 
 
!> @brief Add a domain decomposed variable.
subroutine register_domain_variable(fileobj, variable_name, variable_type, dimensions, chunksizes)
 
  type(fmsnetcdfdomainfile_t), intent(inout) :: fileobj !< File object.
  character(len=*), intent(in) :: variable_name !< Variable name.
  character(len=*), intent(in) :: variable_type !< Variable type.  Allowed
                                                !! values are: "int", "int64",
                                                !! "float", or "double".
  character(len=*), dimension(:), intent(in), optional :: dimensions !< Dimension names.
  integer, intent(in), optional :: chunksizes(:) !< netcdf chunksize to use for this variable (netcdf4 only)
 
  if (.not. fileobj%is_readonly) then
    call netcdf_add_variable(fileobj, variable_name, variable_type, dimensions, chunksizes)
    if (present(dimensions)) then
      if (size(dimensions) .eq. 1) then
        call add_domain_attribute(fileobj, variable_name)
      endif
    endif
  endif
end subroutine register_domain_variable
 
 
!> @brief Loop through registered restart variables and write them to
!!        a netcdf file.
subroutine save_domain_restart(fileobj, unlim_dim_level)
 
  type(fmsnetcdfdomainfile_t), intent(in) :: fileobj !< File object.
  integer, intent(in), optional :: unlim_dim_level !< Unlimited dimension level.
 
  integer :: i
  character(len=32) :: chksum
  logical :: is_decomposed
 
  if (.not. fileobj%is_restart) then
    call error("file "//trim(fileobj%path)// &
             & " is not a restart file. You must set is_restart=.true. in your open_file call.")
  endif
 
! Calculate the variable's checksum and write it to the netcdf file
  do i = 1, fileobj%num_restart_vars
    if (associated(fileobj%restart_vars(i)%data2d)) then
      chksum = compute_global_checksum(fileobj, fileobj%restart_vars(i)%varname, &
                                       fileobj%restart_vars(i)%data2d, is_decomposed)
      if (is_decomposed) then
        call register_variable_attribute(fileobj, fileobj%restart_vars(i)%varname, &
                                         "checksum", chksum(1:len(chksum)), str_len=len(chksum))
      endif
    elseif (associated(fileobj%restart_vars(i)%data3d)) then
      chksum = compute_global_checksum(fileobj, fileobj%restart_vars(i)%varname, &
                                       fileobj%restart_vars(i)%data3d, is_decomposed)
      if (is_decomposed) then
        call register_variable_attribute(fileobj, fileobj%restart_vars(i)%varname, &
                                         "checksum", chksum(1:len(chksum)), str_len=len(chksum))
      endif
    elseif (associated(fileobj%restart_vars(i)%data4d)) then
      chksum = compute_global_checksum(fileobj, fileobj%restart_vars(i)%varname, &
                                       fileobj%restart_vars(i)%data4d, is_decomposed)
      if (is_decomposed) then
        call register_variable_attribute(fileobj, fileobj%restart_vars(i)%varname, &
                                         "checksum", chksum(1:len(chksum)), str_len=len(chksum))
      endif
    endif
  enddo
 
! Write the variable's data to the netcdf file
  do i = 1, fileobj%num_restart_vars
    if (associated(fileobj%restart_vars(i)%data0d)) then
      call domain_write_0d(fileobj, fileobj%restart_vars(i)%varname, &
                           fileobj%restart_vars(i)%data0d, unlim_dim_level=unlim_dim_level)
    elseif (associated(fileobj%restart_vars(i)%data1d)) then
      call domain_write_1d(fileobj, fileobj%restart_vars(i)%varname, &
                           fileobj%restart_vars(i)%data1d, unlim_dim_level=unlim_dim_level)
    elseif (associated(fileobj%restart_vars(i)%data2d)) then
      call domain_write_2d(fileobj, fileobj%restart_vars(i)%varname, &
                           fileobj%restart_vars(i)%data2d, unlim_dim_level=unlim_dim_level)
    elseif (associated(fileobj%restart_vars(i)%data3d)) then
      call domain_write_3d(fileobj, fileobj%restart_vars(i)%varname, &
                           fileobj%restart_vars(i)%data3d, unlim_dim_level=unlim_dim_level)
    elseif (associated(fileobj%restart_vars(i)%data4d)) then
      call domain_write_4d(fileobj, fileobj%restart_vars(i)%varname, &
                           fileobj%restart_vars(i)%data4d, unlim_dim_level=unlim_dim_level)
    else
      call error("This routine only accepts data that is scalar, 1d 2d 3d or 4d."//&
                 " The data sent in has an unsupported dimensionality")
    endif
  enddo
 
end subroutine save_domain_restart
 
 
!> @brief Loop through registered restart variables and read them from
!!        a netcdf file.
subroutine restore_domain_state(fileobj, unlim_dim_level, ignore_checksum)
 
  type(fmsnetcdfdomainfile_t), intent(inout) :: fileobj !< File object.
  integer, intent(in), optional :: unlim_dim_level !< Unlimited dimension level.
  logical, intent(in), optional :: ignore_checksum !< Checksum data integrity flag.
 
  integer :: i
  character(len=32) :: chksum_in_file
  character(len=32) :: chksum
  logical :: chksum_ignore = .false. !< local variable for data integrity checks
                                     !! default: .FALSE. - checks enabled
  logical :: is_decomposed
 
  if (PRESENT(ignore_checksum)) chksum_ignore = ignore_checksum
 
  if (.not. fileobj%is_restart) then
    call error("file "//trim(fileobj%path)// &
             & " is not a restart file. You must set is_restart=.true. in your open_file call.")
  endif
  do i = 1, fileobj%num_restart_vars
    if (associated(fileobj%restart_vars(i)%data0d)) then
      call domain_read_0d(fileobj, fileobj%restart_vars(i)%varname, &
                          fileobj%restart_vars(i)%data0d, unlim_dim_level=unlim_dim_level)
    elseif (associated(fileobj%restart_vars(i)%data1d)) then
      call domain_read_1d(fileobj, fileobj%restart_vars(i)%varname, &
                          fileobj%restart_vars(i)%data1d, unlim_dim_level=unlim_dim_level)
    elseif (associated(fileobj%restart_vars(i)%data2d)) then
      call domain_read_2d(fileobj, fileobj%restart_vars(i)%varname, &
                          fileobj%restart_vars(i)%data2d, unlim_dim_level=unlim_dim_level)
      if (.not.chksum_ignore) then
        chksum = compute_global_checksum(fileobj, fileobj%restart_vars(i)%varname, &
                                         fileobj%restart_vars(i)%data2d, is_decomposed)
        if (variable_att_exists(fileobj, fileobj%restart_vars(i)%varname, "checksum") .and. &
            is_decomposed) then
          call get_variable_attribute(fileobj, fileobj%restart_vars(i)%varname, &
                                      "checksum", chksum_in_file)
          if (.not. string_compare(trim(adjustl(chksum_in_file)), trim(adjustl(chksum)))) then
            call error("The checksum in the file:"//trim(fileobj%path)//" and variable:"// &
                      & trim(fileobj%restart_vars(i)%varname)// &
                      &" does not match the checksum calculated from the data. file:"//trim(adjustl(chksum_in_file))//&
                      &" from data:"//trim(adjustl(chksum)))
          endif
        endif
      endif
    elseif (associated(fileobj%restart_vars(i)%data3d)) then
      call domain_read_3d(fileobj, fileobj%restart_vars(i)%varname, &
                          fileobj%restart_vars(i)%data3d, unlim_dim_level=unlim_dim_level)
      if (.not.chksum_ignore) then
        chksum = compute_global_checksum(fileobj, fileobj%restart_vars(i)%varname, &
                                         fileobj%restart_vars(i)%data3d, is_decomposed)
        if (variable_att_exists(fileobj, fileobj%restart_vars(i)%varname, "checksum") .and. &
            is_decomposed) then
          call get_variable_attribute(fileobj, fileobj%restart_vars(i)%varname, &
                                      "checksum", chksum_in_file(1:len(chksum_in_file)))
          if (.not. string_compare(trim(adjustl(chksum_in_file)), trim(adjustl(chksum)))) then
            call error("The checksum in the file:"//trim(fileobj%path)//" and variable:"// &
                      & trim(fileobj%restart_vars(i)%varname)//&
                      &" does not match the checksum calculated from the data. file:"//trim(adjustl(chksum_in_file))//&
                      &" from data:"//trim(adjustl(chksum)))
          endif
        endif
      endif
    elseif (associated(fileobj%restart_vars(i)%data4d)) then
      call domain_read_4d(fileobj, fileobj%restart_vars(i)%varname, &
                          fileobj%restart_vars(i)%data4d, unlim_dim_level=unlim_dim_level)
      if (.not.chksum_ignore) then
        chksum = compute_global_checksum(fileobj, fileobj%restart_vars(i)%varname, &
                                         fileobj%restart_vars(i)%data4d, is_decomposed)
        if (variable_att_exists(fileobj, fileobj%restart_vars(i)%varname, "checksum") .and. &
            is_decomposed) then
          call get_variable_attribute(fileobj, fileobj%restart_vars(i)%varname, &
                                      "checksum", chksum_in_file)
          if (.not. string_compare(trim(adjustl(chksum_in_file)), trim(adjustl(chksum)))) then
            call error("The checksum in the file:"//trim(fileobj%path)//" and variable:"// &
                      & trim(fileobj%restart_vars(i)%varname)//&
                      &" does not match the checksum calculated from the data. file:"//trim(adjustl(chksum_in_file))//&
                      &" from data:"//trim(adjustl(chksum)))
          endif
        endif
      endif
    else
      call error("There is no data associated with the variable: "//trim(fileobj%restart_vars(i)%varname)//&
                 &" and the file: "//trim(fileobj%path)//". Check your register_restart_variable call")
    endif
  enddo
end subroutine restore_domain_state
 
 
!> @brief Return an array of compute domain indices.
subroutine get_compute_domain_dimension_indices(fileobj, dimname, indices)
 
  type(fmsnetcdfdomainfile_t), intent(in) :: fileobj !< File object.
  character(len=*), intent(in) :: dimname !< Name of dimension variable.
  integer, dimension(:), allocatable, intent(inout) :: indices !< Compute domain indices.
 
  type(domain2d), pointer :: io_domain
  integer :: dpos
  integer :: s
  integer :: e
  integer :: i
 
  io_domain => mpp_get_io_domain(fileobj%domain)
  dpos = get_domain_decomposed_index(dimname, fileobj%xdims, fileobj%nx)
  if (dpos .ne. variable_not_found) then
    dpos = fileobj%xdims(dpos)%pos
    call mpp_get_compute_domain(io_domain, xbegin=s, xend=e, position=dpos)
  else
    dpos = get_domain_decomposed_index(dimname, fileobj%ydims, fileobj%ny)
    if (dpos .ne. variable_not_found) then
      dpos = fileobj%ydims(dpos)%pos
      call mpp_get_compute_domain(io_domain, ybegin=s, yend=e, position=dpos)
    else
      call error("get_compute_domain_dimension_indices: the input dimension:"//trim(dimname)// &
               & " is not domain decomposed.")
    endif
  endif
  if (allocated(indices)) then
    deallocate(indices)
  endif
  allocate(indices(e-s+1))
  do i = s, e
    indices(i-s+1) = i
  enddo
end subroutine get_compute_domain_dimension_indices
 
 
!> @brief Utility routine that retrieves domain indices.
!! @internal
subroutine domain_offsets(data_xsize, data_ysize, domain, xpos, ypos, &
                          isd, isc, xc_size, jsd, jsc, yc_size, &
                          buffer_includes_halos, extra_x_point, &
                          extra_y_point, msg)
 
  integer, intent(in) :: data_xsize !< Size of buffer's domain "x" dimension.
  integer, intent(in) :: data_ysize !< Size of buffer's domain "y" dimension.
  type(domain2d), intent(in) :: domain !< Parent domain.
  integer, intent(in) :: xpos !< Variable's domain x dimension position.
  integer, intent(in) :: ypos !< Variable's domain y dimension position.
  integer, intent(out) :: isd !< Starting index for x dimension of data domain.
  integer, intent(out) :: isc !< Starting index for x dimension of compute domain.
  integer, intent(out) :: xc_size !< Size of x dimension of compute domain.
  integer, intent(out) :: jsd !< Starting index for y dimension of data domain.
  integer, intent(out) :: jsc !< Starting index for y dimension of compute domain.
  integer, intent(out) :: yc_size !< Size of y dimension of compute domain.
  logical, intent(out) :: buffer_includes_halos !< Flag telling if input buffer includes space for halos.
  logical, intent(out), optional :: extra_x_point !< Flag indicating if data_array has an extra point in x
  logical, intent(out), optional :: extra_y_point !< Flag indicating if data_array has an extra point in y
  character(len=*), intent(in), optional :: msg !< Message appended to fatal error
 
  integer :: xd_size
  integer :: yd_size
  integer :: iec
  integer :: xmax
  integer :: jec
  integer :: ymax
  type(domain2d), pointer :: io_domain !< I/O domain variable is decomposed over.
 
  io_domain => mpp_get_io_domain(domain)
 
  call mpp_get_global_domain(domain, xend=xmax, position=xpos)
  call mpp_get_global_domain(domain, yend=ymax, position=ypos)
  call mpp_get_data_domain(io_domain, xbegin=isd, xsize=xd_size, position=xpos)
  call mpp_get_data_domain(io_domain, ybegin=jsd, ysize=yd_size, position=ypos)
 
  call mpp_get_compute_domain(io_domain, xbegin=isc, xend=iec, xsize=xc_size, &
                              position=xpos)
  ! If the xpos is east and the ending x index is NOT equal to max allowed, set extra_x_point to true
  if (present(extra_x_point)) then
    if ((xpos .eq. east) .and. (iec .ne. xmax)) then
      extra_x_point = .true.
    else
      extra_x_point = .false.
    endif
  endif
 
  call mpp_get_compute_domain(io_domain, ybegin=jsc, yend=jec, ysize=yc_size, &
                              position=ypos)
  ! If the ypost is north and the ending y index is NOT equal to max allowed, set extra_y_point to true
  if (present(extra_y_point)) then
    if ((ypos .eq. north) .and. (jec .ne. ymax)) then
      extra_y_point = .true.
    else
      extra_y_point = .false.
    endif
  endif
 
  buffer_includes_halos = (data_xsize .eq. xd_size) .and. (data_ysize .eq. yd_size)
  if (.not. buffer_includes_halos .and. data_xsize .ne. xc_size .and. data_ysize &
      .ne. yc_size) then
     print *, "buffer_includes_halos:", buffer_includes_halos, " data_xsize:", &
     data_xsize, " xc_size:", xc_size, " data_ysize:", data_ysize, " yc_size:", &
     yc_size
    call error(trim(msg)//" The data is not on the compute domain or the data domain")
  endif
end subroutine domain_offsets
 
 
!> @brief Get starting/ending global indices of the I/O domain for a domain decomposed
!!        file.
subroutine get_global_io_domain_indices(fileobj, dimname, is, ie, indices)
 
  type(fmsnetcdfdomainfile_t), intent(in) :: fileobj !< File object.
  character(len=*), intent(in) :: dimname !< Name of dimension variable.
  integer, intent(out) :: is !< Starting index of I/O global domain.
  integer, intent(out) :: ie !< Ending index of I/O global domain.
  integer, dimension(:), allocatable, intent(out), optional :: indices !< Global domain indices
 
  type(domain2d), pointer :: io_domain
  integer :: dpos
  integer :: i
 
  io_domain => mpp_get_io_domain(fileobj%domain)
  dpos = get_domain_decomposed_index(dimname, fileobj%xdims, fileobj%nx)
  if (dpos .ne. variable_not_found) then
    dpos = fileobj%xdims(dpos)%pos
    call mpp_get_global_domain(io_domain, xbegin=is, xend=ie, position=dpos)
  else
    dpos = get_domain_decomposed_index(dimname, fileobj%ydims, fileobj%ny)
    if (dpos .ne. variable_not_found) then
      dpos = fileobj%ydims(dpos)%pos
      call mpp_get_global_domain(io_domain, ybegin=is, yend=ie, position=dpos)
    else
      call error("get_global_io_domain_indices: the dimension "//trim(dimname)//" in the file: "//trim(fileobj%path)//&
                 &" is not domain decomposed. Check your register_axis call")
    endif
  endif
 
! Allocate indices to the difference between the ending and starting indices and
! fill indices with the data
  if (present(indices)) then
    if(allocated(indices)) then
      call error("get_global_io_domain_indices: the variable indices should not be allocated.")
    endif
    allocate(indices(ie-is+1))
    do i = is, ie
      indices(i-is+1) = i
    enddo
  endif
 
 
end subroutine get_global_io_domain_indices
 
!> @brief Read a mosaic_file and get the grid filename for the current tile or
!!        for the tile specified
subroutine get_mosaic_tile_grid(grid_file,mosaic_file, domain, tile_count)
  character(len=*), intent(out)          :: grid_file !< Filename of the grid file for the
                                                      !! current domain tile or for tile
                                                      !! specified in tile_count
  character(len=*), intent(in)           :: mosaic_file !< Filename that will be read
  type(domain2d),   intent(in)           :: domain !< Input domain
  integer,          intent(in), optional :: tile_count !< Optional argument indicating
                                                       !! the tile you want grid file name for
                                                       !! this is for when a pe is in more than
                                                       !! tile.
  integer                                :: tile !< Current domian tile or tile_count
  integer                                :: ntileme !< Total number of tiles in the domain
  integer, dimension(:), allocatable     :: tile_id !< List of tiles in the domain
  type(fmsnetcdffile_t)                  :: fileobj !< Fms2io file object
 
  tile = 1
  if(present(tile_count)) tile = tile_count
  ntileme = mpp_get_current_ntile(domain)
  allocate(tile_id(ntileme))
  tile_id = mpp_get_tile_id(domain)
 
  if (netcdf_file_open(fileobj, mosaic_file, "read")) then
      call netcdf_read_data(fileobj, "gridfiles", grid_file, corner=tile_id(tile))
      grid_file = 'INPUT/'//trim(grid_file)
      call netcdf_file_close(fileobj)
  endif
 
end subroutine get_mosaic_tile_grid
 
include "register_domain_restart_variable.inc"
include "domain_read.inc"
include "domain_write.inc"
#include "compute_global_checksum.inc"
 
 
end module fms_netcdf_domain_io_mod
!> @}
! close documentation grouping