fms2_io/fms_netcdf_domain_io.F90

   1 !***********************************************************************
   2 !*                   GNU Lesser General Public License
   3 !*
   4 !* This file is part of the GFDL Flexible Modeling System (FMS).
   5 !*
   6 !* FMS is free software: you can redistribute it and/or modify it under
   7 !* the terms of the GNU Lesser General Public License as published by
   8 !* the Free Software Foundation, either version 3 of the License, or (at
   9 !* your option) any later version.
  10 !*
  11 !* FMS is distributed in the hope that it will be useful, but WITHOUT
  12 !* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  13 !* FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
  14 !* for more details.
  15 !*
  16 !* You should have received a copy of the GNU Lesser General Public
  17 !* License along with FMS.  If not, see <http://www.gnu.org/licenses/>.
  18 !***********************************************************************
  19 !> @defgroup fms_netcdf_domain_io_mod fms_netcdf_domain_io_mod
  20 !> @ingroup fms2_io
  21 !> @brief Domain-specific I/O wrappers.
  22
  23 !> @addtogroup fms_netcdf_domain_io_mod
  24 !> @{
  25 module fms_netcdf_domain_io_mod
  26 use netcdf
  27 use mpp_mod
  28 use mpp_domains_mod
  29 use fms_io_utils_mod
  30 use netcdf_io_mod
  31 use platform_mod
  32 implicit none
  33 private
  34
  35
  36 !Module constants.
  37 integer, parameter :: no_domain_decomposed_dimension = 0
  38 integer, parameter, public :: max_num_domain_decomposed_dims = 10
  39 integer, parameter :: variable_not_found = 0
  40 integer, parameter :: default_domain_position = center
  41 character(len=16), parameter :: domain_pos_att = "domain_position"
  42 character(len=16), parameter :: domain_axis_att_name = "domain_axis"
  43 character(len=16), parameter :: x = "x"
  44 character(len=16), parameter :: y = "y"
  45
  46 !> @}
  47
  48 !> @brief Domain variable.
  49 !> @ingroup fms_netcdf_domain_io_mod
  50 type, private :: DomainDimension_t
  51   character(len=nf90_max_name) :: varname !< Variable name.
  52   integer :: pos !< Domain position.
  53 endtype DomainDimension_t
  54
  55
  56 !> @brief netcdf domain file type.
  57 !> @ingroup fms_netcdf_domain_io_mod
  58 type, extends(FmsNetcdfFile_t), public :: FmsNetcdfDomainFile_t
  59   type(domain2d) :: domain !< Two-dimensional domain.
  60   type(DomainDimension_t), dimension(:), allocatable :: xdims !< Dimensions associated
  61                                                               !! with the "x" axis
  62                                                               !! of a 2d domain.
  63   integer :: nx !< Number of "x" dimensions.
  64   type(DomainDimension_t), dimension(:), allocatable :: ydims !< Dimensions associated
  65                                                               !! with the "y" axis
  66                                                               !! of a 2d domain.
  67   integer :: ny !< Number of "y" dimensions.
  68   character(len=256) :: non_mangled_path !< Non-domain-mangled file path.
  69   logical :: adjust_indices !< Flag telling if indices need to be adjusted
  70                             !! for domain-decomposed read.
  71 endtype FmsNetcdfDomainFile_t
  72
  73
  74 public :: open_domain_file
  75 public :: close_domain_file
  76 public :: register_domain_decomposed_dimension
  77 public :: register_domain_variable
  78 public :: register_domain_restart_variable_0d
  79 public :: register_domain_restart_variable_1d
  80 public :: register_domain_restart_variable_2d
  81 public :: register_domain_restart_variable_3d
  82 public :: register_domain_restart_variable_4d
  83 public :: register_domain_restart_variable_5d
  84 public :: domain_read_0d
  85 public :: domain_read_1d
  86 public :: domain_read_2d
  87 public :: domain_read_3d
  88 public :: domain_read_4d
  89 public :: domain_read_5d
  90 public :: domain_write_0d
  91 public :: domain_write_1d
  92 public :: domain_write_2d
  93 public :: domain_write_3d
  94 public :: domain_write_4d
  95 public :: domain_write_5d
  96 public :: save_domain_restart
  97 public :: restore_domain_state
  98 public :: get_compute_domain_dimension_indices
  99 public :: get_global_io_domain_indices
 100 public :: is_dimension_registered
 101 public :: get_mosaic_tile_grid
 102
 103 !> @ingroup fms_netcdf_domain_io_mod
 104 interface compute_global_checksum
 105   module procedure compute_global_checksum_2d
 106   module procedure compute_global_checksum_3d
 107   module procedure compute_global_checksum_4d
 108 end interface compute_global_checksum
 109
 110 !> @addtogroup fms_netcdf_domain_io_mod
 111 !> @{
 112
 113 contains
 114
 115
 116 !> @brief Get the index of a domain decomposed dimension.
 117 !! @return Index of domain decomposed dimension.
 118 function get_domain_decomposed_index(name_, array, size_) &
 119   result(index_)
 120
 121   character(len=*), intent(in) :: name_ !< Name.
 122   type(DomainDimension_t), dimension(:), intent(in) :: array !< Array to search through.
 123   integer, intent(in) :: size_ !< Number of spots to look in.
 124   integer :: index_
 125
 126   integer :: i
 127
 128   index_ = variable_not_found
 129   do i = 1, size_
 130     if (string_compare(array(i)%varname, name_)) then
 131       index_ = i
 132       return
 133     endif
 134   enddo
 135 end function get_domain_decomposed_index
 136
 137
 138 !> @brief Add a domain decomposed dimension to an array.
 139 subroutine append_domain_decomposed_dimension(name_, position_, array, size_)
 140
 141   character(len=*), intent(in) :: name_ !< Variable name.
 142   integer, intent(in) :: position_ !< Domain position.
 143   type(DomainDimension_t), dimension(:), intent(inout) :: array !< Array to search through.
 144   integer, intent(inout) :: size_ !< Number of spots to look in.
 145
 146   integer :: i
 147
 148   do i = 1, size_
 149     if (string_compare(array(i)%varname, name_)) then
 150       call error("variable "//trim(name_)//" already registered.")
 151     endif
 152   enddo
 153   size_ = size_ + 1
 154   if (size_ .gt. size(array)) then
 155     call error("number of domain decomposed variables exceeds limit.")
 156   endif
 157   call string_copy(array(size_)%varname, name_)
 158   array(size_)%pos = position_
 159 end subroutine append_domain_decomposed_dimension
 160
 161
 162 !> @brief Given a domain decomposed dimension, get its domain position.
 163 !! @return Position of the domain decomposed variable.
 164 function get_domain_position(name_, array, size_) &
 165   result(dpos)
 166
 167   character(len=*), intent(in) :: name_ !< Variable name.
 168   type(DomainDimension_t), dimension(:), intent(in) :: array !< Array to search through.
 169   integer, intent(in) :: size_
 170   integer :: dpos
 171
 172   dpos = get_domain_decomposed_index(name_, array, size_)
 173   if (dpos .ne. variable_not_found) then
 174     dpos = array(dpos)%pos
 175   endif
 176 end function get_domain_position
 177
 178
 179 !> @brief Given a variable, get the index of the "x" or "y" domain decomposed
 180 !!        dimension.
 181 !! @return Index of the domain decomposed dimension or else
 182 !!         no_domain_decomposed_dimension.
 183 function get_domain_decomposed_dimension_index(fileobj, variable_name, &
 184                                                xory, broadcast) &
 185   result(index_)
 186
 187   type(FmsNetcdfDomainFile_t), intent(in), target :: fileobj !< File object.
 188   character(len=*), intent(in) :: variable_name !< Variable name.
 189   character(len=*), intent(in) :: xory !< String telling which dimension to
 190                                        !! look for.  Valid values are "x"
 191                                        !! or "y".
 192   logical, intent(in), optional :: broadcast !< Flag controlling whether or
 193                                              !! not the index will be
 194                                              !! broadcast to non "I/O root"
 195                                              !! ranks.  The broadcast will
 196                                              !! be done by default.
 197   integer :: index_
 198
 199   integer :: ndims
 200   character(len=nf90_max_name), dimension(:), allocatable :: dim_names
 201   type(DomainDimension_t), dimension(:), pointer :: p
 202   integer :: n
 203   integer :: i
 204
 205   index_ = no_domain_decomposed_dimension
 206   if (fileobj%is_root) then
 207     ndims = get_variable_num_dimensions(fileobj, variable_name, broadcast=.false.)
 208     allocate(dim_names(ndims))
 209     dim_names(:) = ""
 210     call get_variable_dimension_names(fileobj, variable_name, dim_names, broadcast=.false.)
 211     if (string_compare(xory, x, .true.)) then
 212       p => fileobj%xdims
 213       n = fileobj%nx
 214     elseif (string_compare(xory, y, .true.)) then
 215       p => fileobj%ydims
 216       n = fileobj%ny
 217     else
 218       call error("unrecognized xory flag value.")
 219     endif
 220     do i = 1, size(dim_names)
 221       if (get_domain_decomposed_index(dim_names(i), p, n) .ne. variable_not_found) then
 222         index_ = i
 223         exit
 224       endif
 225     enddo
 226     deallocate(dim_names)
 227   endif
 228   if (present(broadcast)) then
 229     if (.not. broadcast) then
 230       return
 231     endif
 232   endif
 233   call mpp_broadcast(index_, fileobj%io_root, pelist=fileobj%pelist)
 234 end function get_domain_decomposed_dimension_index
 235
 236
 237 !> @brief Determine if a variable is "domain decomposed."
 238 !! @return Flag telling if the variable is "domain decomposed."
 239 function is_variable_domain_decomposed(fileobj, variable_name, broadcast, &
 240                                        xindex, yindex, xpos, ypos) &
 241   result(is_decomposed)
 242
 243   type(FmsNetcdfDomainFile_t), intent(in) :: fileobj !< File object.
 244   character(len=*), intent(in) :: variable_name !< Variable name.
 245   logical, intent(in), optional :: broadcast !< Flag controlling whether or
 246                                              !! not the index will be
 247                                              !! broadcast to non "I/O root"
 248                                              !! ranks.  The broadcast will
 249                                              !! be done by default.
 250   integer, intent(out), optional :: xindex !< The index of the domain
 251                                            !! x dimension.
 252   integer, intent(out), optional :: yindex !< The index of the domain
 253                                            !! y dimension.
 254   integer, intent(out), optional :: xpos !< Domain position of the x dimension.
 255   integer, intent(out), optional :: ypos !< Domain position of the y dimension.
 256   logical :: is_decomposed
 257
 258   integer, dimension(2) :: indices
 259   integer :: ndims
 260   character(len=nf90_max_name), dimension(:), allocatable :: dim_names
 261
 262   indices(1) = get_domain_decomposed_dimension_index(fileobj, variable_name, x, broadcast)
 263   if (present(xindex)) then
 264     xindex = indices(1)
 265   endif
 266   indices(2) = get_domain_decomposed_dimension_index(fileobj, variable_name, y, broadcast)
 267   if (present(yindex)) then
 268     yindex = indices(2)
 269   endif
 270   is_decomposed = (indices(1) .ne. no_domain_decomposed_dimension) .and. &
 271                   (indices(2) .ne. no_domain_decomposed_dimension)
 272   if (is_decomposed) then
 273     if (.not. present(xpos) .and. .not. present(ypos)) then
 274       return
 275     endif
 276     ndims = get_variable_num_dimensions(fileobj, variable_name, broadcast)
 277     allocate(dim_names(ndims))
 278     dim_names(:) = ""
 279     call get_variable_dimension_names(fileobj, variable_name, dim_names, broadcast)
 280     if (present(xpos)) then
 281       xpos = get_domain_position(dim_names(indices(1)), fileobj%xdims, fileobj%nx)
 282     endif
 283     if (present(ypos)) then
 284       ypos = get_domain_position(dim_names(indices(2)), fileobj%ydims, fileobj%ny)
 285     endif
 286     deallocate(dim_names)
 287   else
 288     if (present(xpos)) then
 289       xpos = -1
 290     endif
 291     if (present(ypos)) then
 292       ypos = -1
 293     endif
 294   endif
 295 end function is_variable_domain_decomposed
 296
 297
 298 !> @brief Determine whether a domain-decomposed dimension has been registered to the file object
 299 !! @return Flag telling if the dimension is registered to the file object
 300 function is_dimension_registered(fileobj, dimension_name) &
 301   result(is_registered)
 302
 303   type(FmsNetcdfDomainFile_t), intent(in) :: fileobj !< File object.
 304   character(len=*), intent(in) :: dimension_name !< Dimension name.
 305
 306   ! local
 307   logical :: is_registered
 308   integer :: dpos
 309
 310   dpos = 0
 311   is_registered = .false.
 312   dpos = get_domain_decomposed_index(dimension_name, fileobj%xdims, fileobj%nx)
 313   if (dpos .ne. variable_not_found) then
 314     is_registered = .true.
 315   else
 316     dpos = get_domain_decomposed_index(dimension_name, fileobj%ydims, fileobj%ny)
 317     if (dpos .ne. variable_not_found) is_registered = .true.
 318   endif
 319
 320 end function is_dimension_registered
 321
 322
 323 !> @brief Open a domain netcdf file.
 324 !! @return Flag telling if the open completed successfully.
 325 function open_domain_file(fileobj, path, mode, domain, nc_format, is_restart, dont_add_res_to_filename) &
 326   result(success)
 327
 328   type(FmsNetcdfDomainFile_t),intent(inout) :: fileobj !< File object.
 329   character(len=*), intent(in) :: path !< File path.
 330   character(len=*), intent(in) :: mode !< File mode.  Allowed values
 331                                        !! are "read", "append", "write", or
 332                                        !! "overwrite".
 333   type(domain2d), intent(in) :: domain !< Two-dimensional domain.
 334   character(len=*), intent(in), optional :: nc_format !< Netcdf format that
 335                                                       !! new files are written
 336                                                       !! as.  Allowed values
 337                                                       !! are: "64bit", "classic",
 338                                                       !! or "netcdf4". Defaults to
 339                                                       !! "64bit".
 340   logical, intent(in), optional :: is_restart !< Flag telling if this file
 341                                               !! is a restart file.  Defaults
 342                                               !! to false.
 343   logical, intent(in), optional :: dont_add_res_to_filename !< Flag indicating not to add
 344                                               !! ".res" to the filename
 345   logical :: success
 346
 347   integer, dimension(2) :: io_layout
 348   integer, dimension(1) :: tile_id
 349   character(len=256) :: combined_filepath
 350   type(domain2d), pointer :: io_domain
 351   character(len=256) :: distributed_filepath
 352   integer :: pelist_size
 353   integer, dimension(:), allocatable :: pelist
 354   logical :: success2
 355   type(FmsNetcdfDomainFile_t) :: fileobj2
 356
 357   !Get the path of a "combined" file.
 358   io_layout = mpp_get_io_domain_layout(domain)
 359   tile_id = mpp_get_tile_id(domain)
 360
 361   !< If the number of tiles is greater than 1 or if the current tile is greater
 362   !than 1 add .tileX. to the filename
 363   if (mpp_get_ntile_count(domain) .gt. 1 .or. tile_id(1) > 1) then
 364     call domain_tile_filepath_mangle(combined_filepath, path, tile_id(1))
 365   else
 366     call string_copy(combined_filepath, path)
 367   endif
 368
 369   !Get the path of a "distributed" file.
 370   io_domain => mpp_get_io_domain(domain)
 371   if (.not. associated(io_domain)) then
 372     call error("The domain associated with the file:"//trim(path)//" does not have an io_domain.")
 373   endif
 374   if (io_layout(1)*io_layout(2) .gt. 1) then
 375     tile_id = mpp_get_tile_id(io_domain)
 376     call io_domain_tile_filepath_mangle(distributed_filepath, combined_filepath, tile_id(1))
 377   else
 378     call string_copy(distributed_filepath, combined_filepath)
 379   endif
 380
 381   !Make sure the input domain has an I/O domain and get its pelist.
 382   pelist_size = mpp_get_domain_npes(io_domain)
 383   allocate(pelist(pelist_size))
 384   call mpp_get_pelist(io_domain, pelist)
 385   fileobj%adjust_indices = .true. !Set the default to true
 386
 387   !Open the distibuted files.
 388   success = netcdf_file_open(fileobj, distributed_filepath, mode, nc_format, pelist, &
 389                              is_restart, dont_add_res_to_filename)
 390   if (string_compare(mode, "read", .true.) .or. string_compare(mode, "append", .true.)) then
 391     if (success) then
 392       if (.not. string_compare(distributed_filepath, combined_filepath)) then
 393         success2 = netcdf_file_open(fileobj2, combined_filepath, mode, nc_format, pelist, &
 394                                     is_restart, dont_add_res_to_filename)
 395         if (success2) then
 396           call error("The domain decomposed file:"//trim(path)// &
 397                    & " contains both combined (*.nc) and distributed files (*.nc.XXXX).")
 398         endif
 399       endif
 400     else
 401       success = netcdf_file_open(fileobj, combined_filepath, mode, nc_format, pelist, &
 402                                  is_restart, dont_add_res_to_filename)
 403       !If the file is combined and the layout is not (1,1) set the adjust_indices flag to false
 404       if (success .and. (io_layout(1)*io_layout(2) .gt. 1)) fileobj%adjust_indices = .false.
 405     endif
 406   endif
 407   if (.not. success) then
 408     deallocate(pelist)
 409     return
 410   endif
 411
 412   !Store/initialize necessary properties.
 413   call string_copy(fileobj%non_mangled_path, path)
 414   fileobj%domain = domain
 415   allocate(fileobj%xdims(max_num_domain_decomposed_dims))
 416   fileobj%nx = 0
 417   allocate(fileobj%ydims(max_num_domain_decomposed_dims))
 418   fileobj%ny = 0
 419   call string_copy(fileobj%non_mangled_path, path)
 420
 421   if (string_compare(mode, "write", .true.) .or. string_compare(mode, "overwrite", .true.)) then
 422     !Add global attribute needed by mppnccombine.
 423     call register_global_attribute(fileobj, "NumFilesInSet", io_layout(1)*io_layout(2))
 424   endif
 425 end function open_domain_file
 426
 427
 428 !> @brief Close a domain netcdf file.
 429 subroutine close_domain_file(fileobj)
 430
 431   type(FmsNetcdfDomainFile_t), intent(inout) :: fileobj !< File object.
 432
 433   call netcdf_file_close(fileobj)
 434   deallocate(fileobj%xdims)
 435   fileobj%nx = 0
 436   deallocate(fileobj%ydims)
 437   fileobj%ny = 0
 438 end subroutine close_domain_file
 439
 440
 441 !> @brief Add a dimension to a file associated with a two-dimensional domain.
 442 subroutine register_domain_decomposed_dimension(fileobj, dim_name, xory, domain_position)
 443
 444   type(FmsNetcdfDomainFile_t), intent(inout) :: fileobj !< File object.
 445   character(len=*), intent(in) :: dim_name !< Dimension name.
 446   character(len=*), intent(in) :: xory !< Flag telling if the dimension
 447                                        !! is associated with the "x" or "y"
 448                                        !! axis of the 2d domain.  Allowed
 449                                        !! values are "x" or "y".
 450   integer, intent(in), optional :: domain_position !< Domain position.
 451
 452   integer :: dpos
 453   type(domain2d), pointer :: io_domain
 454   integer :: domain_size
 455   integer :: dim_size
 456
 457   dpos = default_domain_position
 458   if (mpp_domain_is_symmetry(fileobj%domain) .and. present(domain_position)) then
 459     dpos = domain_position
 460   endif
 461   io_domain => mpp_get_io_domain(fileobj%domain)
 462   if (string_compare(xory, x, .true.)) then
 463     if (dpos .ne. center .and. dpos .ne. east) then
 464       call error("Only domain_position=center or domain_position=EAST is supported for x dimensions."// &
 465                   & " Fix your register_axis call for file:"&
 466                   &//trim(fileobj%path)//" and dimension:"//trim(dim_name))
 467     endif
 468     call mpp_get_global_domain(io_domain, xsize=domain_size, position=dpos)
 469     call append_domain_decomposed_dimension(dim_name, dpos, fileobj%xdims, fileobj%nx)
 470   elseif (string_compare(xory, y, .true.)) then
 471     if (dpos .ne. center .and. dpos .ne. north) then
 472       call error("Only domain_position=center or domain_position=NORTH is supported for y dimensions."// &
 473                   & " Fix your register_axis call for file:"&
 474                   &//trim(fileobj%path)//" and dimension:"//trim(dim_name))
 475     endif
 476     call mpp_get_global_domain(io_domain, ysize=domain_size, position=dpos)
 477     call append_domain_decomposed_dimension(dim_name, dpos, fileobj%ydims, fileobj%ny)
 478   else
 479     call error("The register_axis call for file:"//trim(fileobj%path)//" and dimension:"//trim(dim_name)// &
 480                & " has an unrecognized xory flag value:"&
 481                &//trim(xory)//" only 'x' and 'y' are allowed.")
 482   endif
 483   if (fileobj%is_readonly .or. (fileobj%mode_is_append .and. dimension_exists(fileobj, dim_name))) then
 484     call get_dimension_size(fileobj, dim_name, dim_size, broadcast=.true.)
 485     if (dim_size .lt. domain_size) then
 486       call error("dimension "//trim(dim_name)//" in the file "//trim(fileobj%path)//" is smaller than the size of" &
 487                  //" the associated domain "//trim(xory)//" axis.")
 488     endif
 489   else
 490     call netcdf_add_dimension(fileobj, dim_name, domain_size, is_compressed=.false.)
 491   endif
 492 end subroutine register_domain_decomposed_dimension
 493
 494
 495 !> @brief Add a "domain_decomposed" attribute to the axis variables because it is
 496 !!        required by mppnccombine.
 497 !! @internal
 498 subroutine add_domain_attribute(fileobj, variable_name)
 499
 500   type(FmsNetcdfDomainFile_t), intent(inout) :: fileobj !< File object.
 501   character(len=*), intent(in) :: variable_name !< Variable name.
 502
 503   type(domain2d), pointer :: io_domain
 504   integer :: dpos
 505   integer :: sg
 506   integer :: eg
 507   integer :: s
 508   integer :: e
 509   integer, dimension(2) :: io_layout !< Io_layout in the fileobj's domain
 510
 511   !< Don't add the "domain_decomposition" variable attribute if the io_layout is
 512   !! 1,1, to avoid frecheck "failures"
 513   io_layout = mpp_get_io_domain_layout(fileobj%domain)
 514   if (io_layout(1)  .eq. 1 .and. io_layout(2) .eq. 1) return
 515
 516   io_domain => mpp_get_io_domain(fileobj%domain)
 517   dpos = get_domain_decomposed_index(variable_name, fileobj%xdims, fileobj%nx)
 518   if (dpos .ne. variable_not_found) then
 519     dpos = fileobj%xdims(dpos)%pos
 520     call mpp_get_global_domain(fileobj%domain, xbegin=sg, xend=eg, position=dpos)
 521     call mpp_get_global_domain(io_domain, xbegin=s, xend=e, position=dpos)
 522     call register_variable_attribute(fileobj, variable_name, "domain_decomposition", &
 523                                      (/sg, eg, s, e/))
 524   else
 525     dpos = get_domain_decomposed_index(variable_name, fileobj%ydims, fileobj%ny)
 526     if (dpos .ne. variable_not_found) then
 527       dpos = fileobj%ydims(dpos)%pos
 528       call mpp_get_global_domain(fileobj%domain, ybegin=sg, yend=eg, position=dpos)
 529       call mpp_get_global_domain(io_domain, ybegin=s, yend=e, position=dpos)
 530       call register_variable_attribute(fileobj, variable_name, "domain_decomposition", &
 531                                        (/sg, eg, s, e/))
 532     endif
 533   endif
 534 end subroutine add_domain_attribute
 535
 536
 537 !> @brief Add a domain decomposed variable.
 538 subroutine register_domain_variable(fileobj, variable_name, variable_type, dimensions)
 539
 540   type(FmsNetcdfDomainFile_t), intent(inout) :: fileobj !< File object.
 541   character(len=*), intent(in) :: variable_name !< Variable name.
 542   character(len=*), intent(in) :: variable_type !< Variable type.  Allowed
 543                                                 !! values are: "int", "int64",
 544                                                 !! "float", or "double".
 545   character(len=*), dimension(:), intent(in), optional :: dimensions !< Dimension names.
 546
 547   if (.not. fileobj%is_readonly) then
 548     call netcdf_add_variable(fileobj, variable_name, variable_type, dimensions)
 549     if (present(dimensions)) then
 550       if (size(dimensions) .eq. 1) then
 551         call add_domain_attribute(fileobj, variable_name)
 552       endif
 553     endif
 554   endif
 555 end subroutine register_domain_variable
 556
 557
 558 !> @brief Loop through registered restart variables and write them to
 559 !!        a netcdf file.
 560 subroutine save_domain_restart(fileobj, unlim_dim_level)
 561
 562   type(FmsNetcdfDomainFile_t), intent(in) :: fileobj !< File object.
 563   integer, intent(in), optional :: unlim_dim_level !< Unlimited dimension level.
 564
 565   integer :: i
 566   character(len=32) :: chksum
 567   logical :: is_decomposed
 568
 569   if (.not. fileobj%is_restart) then
 570     call error("file "//trim(fileobj%path)// &
 571              & " is not a restart file. You must set is_restart=.true. in your open_file call.")
 572   endif
 573
 574 ! Calculate the variable's checksum and write it to the netcdf file
 575   do i = 1, fileobj%num_restart_vars
 576     if (associated(fileobj%restart_vars(i)%data2d)) then
 577       chksum = compute_global_checksum(fileobj, fileobj%restart_vars(i)%varname, &
 578                                        fileobj%restart_vars(i)%data2d, is_decomposed)
 579       if (is_decomposed) then
 580         call register_variable_attribute(fileobj, fileobj%restart_vars(i)%varname, &
 581                                          "checksum", chksum(1:len(chksum)), str_len=len(chksum))
 582       endif
 583     elseif (associated(fileobj%restart_vars(i)%data3d)) then
 584       chksum = compute_global_checksum(fileobj, fileobj%restart_vars(i)%varname, &
 585                                        fileobj%restart_vars(i)%data3d, is_decomposed)
 586       if (is_decomposed) then
 587         call register_variable_attribute(fileobj, fileobj%restart_vars(i)%varname, &
 588                                          "checksum", chksum(1:len(chksum)), str_len=len(chksum))
 589       endif
 590     elseif (associated(fileobj%restart_vars(i)%data4d)) then
 591       chksum = compute_global_checksum(fileobj, fileobj%restart_vars(i)%varname, &
 592                                        fileobj%restart_vars(i)%data4d, is_decomposed)
 593       if (is_decomposed) then
 594         call register_variable_attribute(fileobj, fileobj%restart_vars(i)%varname, &
 595                                          "checksum", chksum(1:len(chksum)), str_len=len(chksum))
 596       endif
 597     endif
 598   enddo
 599
 600 ! Write the variable's data to the netcdf file
 601   do i = 1, fileobj%num_restart_vars
 602     if (associated(fileobj%restart_vars(i)%data0d)) then
 603       call domain_write_0d(fileobj, fileobj%restart_vars(i)%varname, &
 604                            fileobj%restart_vars(i)%data0d, unlim_dim_level=unlim_dim_level)
 605     elseif (associated(fileobj%restart_vars(i)%data1d)) then
 606       call domain_write_1d(fileobj, fileobj%restart_vars(i)%varname, &
 607                            fileobj%restart_vars(i)%data1d, unlim_dim_level=unlim_dim_level)
 608     elseif (associated(fileobj%restart_vars(i)%data2d)) then
 609       call domain_write_2d(fileobj, fileobj%restart_vars(i)%varname, &
 610                            fileobj%restart_vars(i)%data2d, unlim_dim_level=unlim_dim_level)
 611     elseif (associated(fileobj%restart_vars(i)%data3d)) then
 612       call domain_write_3d(fileobj, fileobj%restart_vars(i)%varname, &
 613                            fileobj%restart_vars(i)%data3d, unlim_dim_level=unlim_dim_level)
 614     elseif (associated(fileobj%restart_vars(i)%data4d)) then
 615       call domain_write_4d(fileobj, fileobj%restart_vars(i)%varname, &
 616                            fileobj%restart_vars(i)%data4d, unlim_dim_level=unlim_dim_level)
 617     else
 618       call error("This routine only accepts data that is scalar, 1d 2d 3d or 4d."//&
 619                  " The data sent in has an unsupported dimensionality")
 620     endif
 621   enddo
 622
 623 end subroutine save_domain_restart
 624
 625
 626 !> @brief Loop through registered restart variables and read them from
 627 !!        a netcdf file.
 628 subroutine restore_domain_state(fileobj, unlim_dim_level, ignore_checksum)
 629
 630   type(FmsNetcdfDomainFile_t), intent(inout) :: fileobj !< File object.
 631   integer, intent(in), optional :: unlim_dim_level !< Unlimited dimension level.
 632   logical, intent(in), optional :: ignore_checksum !< Checksum data integrity flag.
 633
 634   integer :: i
 635   character(len=32) :: chksum_in_file
 636   character(len=32) :: chksum
 637   logical :: chksum_ignore = .FALSE. !< local variable for data integrity checks
 638                                      !! default: .FALSE. - checks enabled
 639   logical :: is_decomposed
 640
 641   if (PRESENT(ignore_checksum)) chksum_ignore = ignore_checksum
 642
 643   if (.not. fileobj%is_restart) then
 644     call error("file "//trim(fileobj%path)// &
 645              & " is not a restart file. You must set is_restart=.true. in your open_file call.")
 646   endif
 647   do i = 1, fileobj%num_restart_vars
 648     if (associated(fileobj%restart_vars(i)%data0d)) then
 649       call domain_read_0d(fileobj, fileobj%restart_vars(i)%varname, &
 650                           fileobj%restart_vars(i)%data0d, unlim_dim_level=unlim_dim_level)
 651     elseif (associated(fileobj%restart_vars(i)%data1d)) then
 652       call domain_read_1d(fileobj, fileobj%restart_vars(i)%varname, &
 653                           fileobj%restart_vars(i)%data1d, unlim_dim_level=unlim_dim_level)
 654     elseif (associated(fileobj%restart_vars(i)%data2d)) then
 655       call domain_read_2d(fileobj, fileobj%restart_vars(i)%varname, &
 656                           fileobj%restart_vars(i)%data2d, unlim_dim_level=unlim_dim_level)
 657       if (.not.chksum_ignore) then
 658         chksum = compute_global_checksum(fileobj, fileobj%restart_vars(i)%varname, &
 659                                          fileobj%restart_vars(i)%data2d, is_decomposed)
 660         if (variable_att_exists(fileobj, fileobj%restart_vars(i)%varname, "checksum") .and. &
 661             is_decomposed) then
 662           call get_variable_attribute(fileobj, fileobj%restart_vars(i)%varname, &
 663                                       "checksum", chksum_in_file)
 664           if (.not. string_compare(trim(adjustl(chksum_in_file)), trim(adjustl(chksum)))) then
 665             call error("The checksum in the file:"//trim(fileobj%path)//" and variable:"// &
 666                       & trim(fileobj%restart_vars(i)%varname)// &
 667                       &" does not match the checksum calculated from the data. file:"//trim(adjustl(chksum_in_file))//&
 668                       &" from data:"//trim(adjustl(chksum)))
 669           endif
 670         endif
 671       endif
 672     elseif (associated(fileobj%restart_vars(i)%data3d)) then
 673       call domain_read_3d(fileobj, fileobj%restart_vars(i)%varname, &
 674                           fileobj%restart_vars(i)%data3d, unlim_dim_level=unlim_dim_level)
 675       if (.not.chksum_ignore) then
 676         chksum = compute_global_checksum(fileobj, fileobj%restart_vars(i)%varname, &
 677                                          fileobj%restart_vars(i)%data3d, is_decomposed)
 678         if (variable_att_exists(fileobj, fileobj%restart_vars(i)%varname, "checksum") .and. &
 679             is_decomposed) then
 680           call get_variable_attribute(fileobj, fileobj%restart_vars(i)%varname, &
 681                                       "checksum", chksum_in_file(1:len(chksum_in_file)))
 682           if (.not. string_compare(trim(adjustl(chksum_in_file)), trim(adjustl(chksum)))) then
 683             call error("The checksum in the file:"//trim(fileobj%path)//" and variable:"// &
 684                       & trim(fileobj%restart_vars(i)%varname)//&
 685                       &" does not match the checksum calculated from the data. file:"//trim(adjustl(chksum_in_file))//&
 686                       &" from data:"//trim(adjustl(chksum)))
 687           endif
 688         endif
 689       endif
 690     elseif (associated(fileobj%restart_vars(i)%data4d)) then
 691       call domain_read_4d(fileobj, fileobj%restart_vars(i)%varname, &
 692                           fileobj%restart_vars(i)%data4d, unlim_dim_level=unlim_dim_level)
 693       if (.not.chksum_ignore) then
 694         chksum = compute_global_checksum(fileobj, fileobj%restart_vars(i)%varname, &
 695                                          fileobj%restart_vars(i)%data4d, is_decomposed)
 696         if (variable_att_exists(fileobj, fileobj%restart_vars(i)%varname, "checksum") .and. &
 697             is_decomposed) then
 698           call get_variable_attribute(fileobj, fileobj%restart_vars(i)%varname, &
 699                                       "checksum", chksum_in_file)
 700           if (.not. string_compare(trim(adjustl(chksum_in_file)), trim(adjustl(chksum)))) then
 701             call error("The checksum in the file:"//trim(fileobj%path)//" and variable:"// &
 702                       & trim(fileobj%restart_vars(i)%varname)//&
 703                       &" does not match the checksum calculated from the data. file:"//trim(adjustl(chksum_in_file))//&
 704                       &" from data:"//trim(adjustl(chksum)))
 705           endif
 706         endif
 707       endif
 708     else
 709       call error("There is no data associated with the variable: "//trim(fileobj%restart_vars(i)%varname)//&
 710                  &" and the file: "//trim(fileobj%path)//". Check your register_restart_variable call")
 711     endif
 712   enddo
 713 end subroutine restore_domain_state
 714
 715
 716 !> @brief Return an array of compute domain indices.
 717 subroutine get_compute_domain_dimension_indices(fileobj, dimname, indices)
 718
 719   type(FmsNetcdfDomainFile_t), intent(in) :: fileobj !< File object.
 720   character(len=*), intent(in) :: dimname !< Name of dimension variable.
 721   integer, dimension(:), allocatable, intent(inout) :: indices !< Compute domain indices.
 722
 723   type(domain2d), pointer :: io_domain
 724   integer :: dpos
 725   integer :: s
 726   integer :: e
 727   integer :: i
 728
 729   io_domain => mpp_get_io_domain(fileobj%domain)
 730   dpos = get_domain_decomposed_index(dimname, fileobj%xdims, fileobj%nx)
 731   if (dpos .ne. variable_not_found) then
 732     dpos = fileobj%xdims(dpos)%pos
 733     call mpp_get_compute_domain(io_domain, xbegin=s, xend=e, position=dpos)
 734   else
 735     dpos = get_domain_decomposed_index(dimname, fileobj%ydims, fileobj%ny)
 736     if (dpos .ne. variable_not_found) then
 737       dpos = fileobj%ydims(dpos)%pos
 738       call mpp_get_compute_domain(io_domain, ybegin=s, yend=e, position=dpos)
 739     else
 740       call error("get_compute_domain_dimension_indices: the input dimension:"//trim(dimname)// &
 741                & " is not domain decomposed.")
 742     endif
 743   endif
 744   if (allocated(indices)) then
 745     deallocate(indices)
 746   endif
 747   allocate(indices(e-s+1))
 748   do i = s, e
 749     indices(i-s+1) = i
 750   enddo
 751 end subroutine get_compute_domain_dimension_indices
 752
 753
 754 !> @brief Utility routine that retrieves domain indices.
 755 !! @internal
 756 subroutine domain_offsets(data_xsize, data_ysize, domain, xpos, ypos, &
 757                           isd, isc, xc_size, jsd, jsc, yc_size, &
 758                           buffer_includes_halos, extra_x_point, &
 759                           extra_y_point, msg)
 760
 761   integer, intent(in) :: data_xsize !< Size of buffer's domain "x" dimension.
 762   integer, intent(in) :: data_ysize !< Size of buffer's domain "y" dimension.
 763   type(domain2d), intent(in) :: domain !< Parent domain.
 764   integer, intent(in) :: xpos !< Variable's domain x dimension position.
 765   integer, intent(in) :: ypos !< Variable's domain y dimension position.
 766   integer, intent(out) :: isd !< Starting index for x dimension of data domain.
 767   integer, intent(out) :: isc !< Starting index for x dimension of compute domain.
 768   integer, intent(out) :: xc_size !< Size of x dimension of compute domain.
 769   integer, intent(out) :: jsd !< Starting index for y dimension of data domain.
 770   integer, intent(out) :: jsc !< Starting index for y dimension of compute domain.
 771   integer, intent(out) :: yc_size !< Size of y dimension of compute domain.
 772   logical, intent(out) :: buffer_includes_halos !< Flag telling if input buffer includes space for halos.
 773   logical, intent(out), optional :: extra_x_point !< Flag indicating if data_array has an extra point in x
 774   logical, intent(out), optional :: extra_y_point !< Flag indicating if data_array has an extra point in y
 775   character(len=*), intent(in), optional :: msg !< Message appended to fatal error
 776
 777   integer :: xd_size
 778   integer :: yd_size
 779   integer :: iec
 780   integer :: xmax
 781   integer :: jec
 782   integer :: ymax
 783   type(domain2d), pointer :: io_domain !< I/O domain variable is decomposed over.
 784
 785   io_domain => mpp_get_io_domain(domain)
 786
 787   call mpp_get_global_domain(domain, xend=xmax, position=xpos)
 788   call mpp_get_global_domain(domain, yend=ymax, position=ypos)
 789   call mpp_get_data_domain(io_domain, xbegin=isd, xsize=xd_size, position=xpos)
 790   call mpp_get_data_domain(io_domain, ybegin=jsd, ysize=yd_size, position=ypos)
 791
 792   call mpp_get_compute_domain(io_domain, xbegin=isc, xend=iec, xsize=xc_size, &
 793                               position=xpos)
 794   ! If the xpos is east and the ending x index is NOT equal to max allowed, set extra_x_point to true
 795   if (present(extra_x_point)) then
 796     if ((xpos .eq. east) .and. (iec .ne. xmax)) then
 797       extra_x_point = .true.
 798     else
 799       extra_x_point = .false.
 800     endif
 801   endif
 802
 803   call mpp_get_compute_domain(io_domain, ybegin=jsc, yend=jec, ysize=yc_size, &
 804                               position=ypos)
 805   ! If the ypost is north and the ending y index is NOT equal to max allowed, set extra_y_point to true
 806   if (present(extra_y_point)) then
 807     if ((ypos .eq. north) .and. (jec .ne. ymax)) then
 808       extra_y_point = .true.
 809     else
 810       extra_y_point = .false.
 811     endif
 812   endif
 813
 814   buffer_includes_halos = (data_xsize .eq. xd_size) .and. (data_ysize .eq. yd_size)
 815   if (.not. buffer_includes_halos .and. data_xsize .ne. xc_size .and. data_ysize &
 816       .ne. yc_size) then
 817      print *, "buffer_includes_halos:", buffer_includes_halos, " data_xsize:", &
 818      data_xsize, " xc_size:", xc_size, " data_ysize:", data_ysize, " yc_size:", &
 819      yc_size
 820     call error(trim(msg)//" The data is not on the compute domain or the data domain")
 821   endif
 822 end subroutine domain_offsets
 823
 824
 825 !> @brief Get starting/ending global indices of the I/O domain for a domain decomposed
 826 !!        file.
 827 subroutine get_global_io_domain_indices(fileobj, dimname, is, ie, indices)
 828
 829   type(FmsNetcdfDomainFile_t), intent(in) :: fileobj !< File object.
 830   character(len=*), intent(in) :: dimname !< Name of dimension variable.
 831   integer, intent(out) :: is !< Staring index of I/O global domain.
 832   integer, intent(out) :: ie !< Ending index of I/O global domain.
 833   integer, dimension(:), allocatable, intent(out), optional :: indices !< Global domain indices
 834
 835   type(domain2d), pointer :: io_domain
 836   integer :: dpos
 837   integer :: i
 838
 839   io_domain => mpp_get_io_domain(fileobj%domain)
 840   dpos = get_domain_decomposed_index(dimname, fileobj%xdims, fileobj%nx)
 841   if (dpos .ne. variable_not_found) then
 842     dpos = fileobj%xdims(dpos)%pos
 843     call mpp_get_global_domain(io_domain, xbegin=is, xend=ie, position=dpos)
 844   else
 845     dpos = get_domain_decomposed_index(dimname, fileobj%ydims, fileobj%ny)
 846     if (dpos .ne. variable_not_found) then
 847       dpos = fileobj%ydims(dpos)%pos
 848       call mpp_get_global_domain(io_domain, ybegin=is, yend=ie, position=dpos)
 849     else
 850       call error("get_global_io_domain_indices: the dimension "//trim(dimname)//" in the file: "//trim(fileobj%path)//&
 851                  &" is not domain decomposed. Check your register_axis call")
 852     endif
 853   endif
 854
 855 ! Allocate indices to the difference between the ending and starting indices and
 856 ! fill indices with the data
 857   if (present(indices)) then
 858     if(allocated(indices)) then
 859       call error("get_global_io_domain_indices: the variable indices should not be allocated.")
 860     endif
 861     allocate(indices(ie-is+1))
 862     do i = is, ie
 863       indices(i-is+1) = i
 864     enddo
 865   endif
 866
 867
 868 end subroutine get_global_io_domain_indices
 869
 870 !> @brief Read a mosaic_file and get the grid filename for the current tile or
 871 !!        for the tile specified
 872 subroutine get_mosaic_tile_grid(grid_file,mosaic_file, domain, tile_count)
 873   character(len=*), intent(out)          :: grid_file !< Filename of the grid file for the
 874                                                       !! current domain tile or for tile
 875                                                       !! specified in tile_count
 876   character(len=*), intent(in)           :: mosaic_file !< Filename that will be read
 877   type(domain2D),   intent(in)           :: domain !< Input domain
 878   integer,          intent(in), optional :: tile_count !< Optional argument indicating
 879                                                        !! the tile you want grid file name for
 880                                                        !! this is for when a pe is in more than
 881                                                        !! tile.
 882   integer                                :: tile !< Current domian tile or tile_count
 883   integer                                :: ntileMe !< Total number of tiles in the domain
 884   integer, dimension(:), allocatable     :: tile_id !< List of tiles in the domain
 885   type(FmsNetcdfFile_t)                  :: fileobj !< Fms2io file object
 886
 887   tile = 1
 888   if(present(tile_count)) tile = tile_count
 889   ntileMe = mpp_get_current_ntile(domain)
 890   allocate(tile_id(ntileMe))
 891   tile_id = mpp_get_tile_id(domain)
 892
 893   if (netcdf_file_open(fileobj, mosaic_file, "read")) then
 894       call netcdf_read_data(fileobj, "gridfiles", grid_file, corner=tile_id(tile))
 895       grid_file = 'INPUT/'//trim(grid_file)
 896       call netcdf_file_close(fileobj)
 897   endif
 898
 899 end subroutine get_mosaic_tile_grid
 900
 901 include "register_domain_restart_variable.inc"
 902 include "domain_read.inc"
 903 include "domain_write.inc"
 904 #include "compute_global_checksum.inc"
 905
 906
 907 end module fms_netcdf_domain_io_mod
 908 !> @}
 909 ! close documentation grouping