MOM_io.F90

! This file is part of MOM6, the Modular Ocean Model version 6.

! See the LICENSE file for licensing information.

! SPDX-License-Identifier: Apache-2.0

!> This module contains I/O framework code

module mom_io

use mom_array_transform,  only : allocate_rotated_array, rotate_array

use mom_array_transform,  only : rotate_array_pair, rotate_vector

use mom_domains,          only : mom_domain_type, domain1d, broadcast, get_domain_components

use mom_domains,          only : rescale_comp_data, num_pes, agrid, bgrid_ne, cgrid_ne

use mom_dyn_horgrid,      only : dyn_horgrid_type

use mom_ensemble_manager, only : get_ensemble_id

use mom_error_handler,    only : mom_error, note, fatal, warning, is_root_pe

use mom_file_parser,      only : log_version, param_file_type

use mom_grid,             only : ocean_grid_type

use mom_io_infra,         only : read_field, read_vector

use mom_io_infra,         only : read_data => read_field ! Deprecated

use mom_io_infra,         only : read_field_chksum

use mom_io_infra,         only : file_exists

use mom_io_infra,         only : open_ascii_file, close_file, file_is_open

use mom_io_infra,         only : get_field_size, field_exists, get_field_atts

use mom_io_infra,         only : get_axis_data, get_filename_suffix

use mom_io_infra,         only : write_version

use mom_io_infra,         only : mom_namelist_file, check_namelist_error, io_infra_init, io_infra_end

use mom_io_infra,         only : append_file, ascii_file, multiple, netcdf_file, overwrite_file

use mom_io_infra,         only : readonly_file, single_file, writeonly_file

use mom_io_infra,         only : center, corner, north_face, east_face

use mom_io_file,          only : mom_file, mom_infra_file, mom_netcdf_file

use mom_io_file,          only : mom_axis, mom_field

use mom_string_functions, only : lowercase, slasher

use mom_verticalgrid,     only : verticalgrid_type

use iso_fortran_env,      only : int32, int64, stdout_iso=>output_unit, stderr_iso=>error_unit

use netcdf,               only : nf90_open, nf90_inq_varid, nf90_inq_varids, nf90_inquire, nf90_close

use netcdf,               only : nf90_inquire_variable, nf90_get_var, nf90_get_att, nf90_inquire_attribute

use netcdf,               only : nf90_strerror, nf90_inquire_dimension

use netcdf,               only : nf90_nowrite, nf90_noerr, nf90_global, nf90_enotatt, nf90_char

! The following are not used in MOM6, but may be used by externals (e.g. SIS2).

use mom_io_infra, only : axistype   ! still used but soon to be nuked

use mom_io_infra, only : fieldtype

use mom_io_infra, only : file_type

use mom_io_infra, only : get_file_info

use mom_io_infra, only : get_file_fields

use mom_io_infra, only : get_file_times

use mom_io_infra, only : open_file

use mom_io_infra, only : write_field

implicit none ; private

! These interfaces are actually implemented in this file.

public :: create_mom_file, reopen_mom_file, cmor_long_std, ensembler, mom_io_init

public :: mom_field

public :: mom_write_field, var_desc, modify_vardesc, query_vardesc, position_from_horgrid

public :: open_namelist_file, check_namelist_error, check_nml_error

public :: get_var_sizes, verify_variable_units, num_timelevels, read_variable, read_attribute

public :: open_file_to_read, close_file_to_read

! The following are simple pass throughs of routines from MOM_io_infra or other modules.

public :: file_exists, open_ascii_file, close_file

public :: mom_file, mom_infra_file, mom_netcdf_file

public :: field_exists, get_filename_appendix

public :: insert_ensemble_appendix

public :: fieldtype, field_size, get_field_atts

public :: axistype, get_axis_data

public :: mom_read_data, mom_read_vector, read_field_chksum

public :: read_netcdf_data

public :: slasher, write_version_number

public :: io_infra_init, io_infra_end

public :: stdout_if_root

public :: get_var_axes_info

public :: get_axis_info

! This is used to set up information descibing non-domain-decomposed axes.

public :: axis_info, set_axis_info, delete_axis_info

! This is used to set up global file attributes

public :: attribute_info, set_attribute_info, delete_attribute_info

! This API is here just to support potential use by non-FMS drivers, and should not persist.

public :: read_data

!> These encoding constants are used to indicate the file format

public :: ascii_file, netcdf_file

!> These encoding constants are used to indicate whether the file is domain decomposed

public :: multiple, single_file

!> These encoding constants are used to indicate the access mode for a file

public :: append_file, overwrite_file, readonly_file, writeonly_file

!> These encoding constants are used to indicate the discretization position of a variable

public :: center, corner, north_face, east_face

! The following are not used in MOM6, but may be used by externals (e.g. SIS2).

public :: create_file

public :: reopen_file

public :: file_type

public :: open_file

public :: get_file_info

public :: get_file_fields

public :: get_file_times

!> Read a field from file using the infrastructure I/O.

interface mom_read_data

  module procedure mom_read_data_0d

  module procedure mom_read_data_0d_int

  module procedure mom_read_data_1d

  module procedure mom_read_data_1d_int

  module procedure mom_read_data_2d

  module procedure mom_read_data_2d_region

  module procedure mom_read_data_3d

  module procedure mom_read_data_3d_region

  module procedure mom_read_data_4d

end interface mom_read_data

!> Read a vector from file using the infrastructure I/O.

interface mom_read_vector

  module procedure mom_read_vector_2d

  module procedure mom_read_vector_3d

end interface mom_read_vector

!> Read a field using native netCDF I/O

!!

!! This function is primarily used for unstructured data which may contain

!! content that cannot be parsed by infrastructure I/O.

interface read_netcdf_data

  ! NOTE: Only 2D I/O is currently used; this should be expanded as needed.

  module procedure read_netcdf_data_2d

end interface read_netcdf_data

!> Write a registered field to an output file, potentially with rotation

interface mom_write_field

  module procedure mom_write_field_legacy_4d

  module procedure mom_write_field_legacy_3d

  module procedure mom_write_field_legacy_2d

  module procedure mom_write_field_legacy_1d

  module procedure mom_write_field_legacy_0d

  module procedure mom_write_field_4d

  module procedure mom_write_field_3d

  module procedure mom_write_field_2d

  module procedure mom_write_field_1d

  module procedure mom_write_field_0d

end interface mom_write_field

!> Read an entire named variable from a named netCDF file using netCDF calls directly, rather

!! than any infrastructure routines and broadcast it from the root PE to the other PEs.

interface read_variable

  module procedure read_variable_0d, read_variable_0d_int

  module procedure read_variable_1d, read_variable_1d_int

  module procedure read_variable_2d, read_variable_3d

end interface read_variable

!> Read a global or variable attribute from a named netCDF file using netCDF calls

!! directly, in some cases reading from the root PE before broadcasting to the other PEs.

interface read_attribute

  module procedure read_attribute_str, read_attribute_real

  module procedure read_attribute_int32, read_attribute_int64

end interface read_attribute

!> Type that stores information that can be used to create a non-decomposed axis.

type :: axis_info

  character(len=32)  :: name = ""       !< The name of this axis for use in files

  character(len=256) :: longname = ""   !< A longer name describing this axis

  character(len=48)  :: units = ""      !< The units of the axis labels

  character(len=8)   :: cartesian = "N" !< A variable indicating which direction

                                        !! this axis corresponds with. Valid values

                                        !! include 'X', 'Y', 'Z', 'T', and 'N' for none.

  integer            :: sense = 0       !< This is 1 for axes whose values increase upward, or -1

                                        !! if they increase downward.  The default, 0, is ignored.

  integer            :: ax_size = 0     !< The number of elements in this axis

  real, allocatable, dimension(:) :: ax_data !< The values of the data on the axis [arbitrary]

end type axis_info

!> Type for describing a 3-d variable for output

type, public :: vardesc

  character(len=64)  :: name               !< Variable name in a NetCDF file

  character(len=48)  :: units              !< Physical dimensions of the variable

  character(len=240) :: longname           !< Long name of the variable

  character(len=8)   :: hor_grid           !< Horizontal grid:  u, v, h, q, Cu, Cv, T, Bu, or 1

  character(len=8)   :: z_grid             !< Vertical grid:  L, i, or 1

  character(len=8)   :: t_grid             !< Time description: s, p, or 1

  character(len=64)  :: cmor_field_name    !< CMOR name

  character(len=64)  :: cmor_units         !< CMOR physical dimensions of the variable

  character(len=240) :: cmor_longname      !< CMOR long name of the variable

  real               :: conversion         !< for unit conversions, such as needed to convert

                                           !! from intensive to extensive [various] or [a A-1 ~> 1]

                                           !! to undo internal dimensional rescaling

  character(len=32)  :: dim_names(5)       !< The names in the file of the axes for this variable

  integer            :: position = -1      !< An integer encoding the horizontal position, it may

                                           !! CENTER, CORNER, EAST_FACE, NORTH_FACE, or 0.

  type(axis_info) :: extra_axes(5)         !< dimensions other than space-time

end type vardesc

!> Type that stores for a global file attribute

type :: attribute_info ; private

  character(len=:), allocatable :: name    !< The name of this attribute

  character(len=:), allocatable :: att_val !< The values of this attribute

end type attribute_info

integer, public :: stdout = stdout_iso  !< standard output unit

integer, public :: stderr = stderr_iso  !< standard output unit

! A note on unit descriptions in comments: MOM6 uses units that can be rescaled for dimensional

! consistency testing. These are noted in comments with units like Z, H, L, and T, along with

! their mks counterparts with notation like "a velocity [Z T-1 ~> m s-1]".  If the units

! vary with the Boussinesq approximation, the Boussinesq variant is given first.

! The functions in this module work with variables with arbitrary units, in which case the

! arbitrary rescaled units are indicated with [A ~> a], while the unscaled units are just [a].

contains

!> `create_MOM_file` wrapper for the legacy file handle, `file_type`.

!! NOTE: This function may be removed in a future release.

subroutine create_file(IO_handle, filename, vars, novars, fields, threading, &

    timeunit, G, dG, GV, checksums, extra_axes, global_atts)

  type(file_type), intent(inout) :: io_handle

    !< Handle for a files or fileset that is to be opened or reopened for

    !! writing

  character(len=*), intent(in) :: filename

    !< full path to the file to create

  type(vardesc), intent(in) :: vars(:)

    !< structures describing fields written to filename

  integer, intent(in) :: novars

    !< number of fields written to filename

  type(fieldtype), intent(inout) :: fields(:)

    !< array of fieldtypes for each variable

  integer, optional, intent(in) :: threading

    !< SINGLE_FILE or MULTIPLE

  real, optional, intent(in) :: timeunit

    !< length of the units for time [s]. The default value is 86400.0, for 1

    !! day.

  type(ocean_grid_type), optional, intent(in) :: g

    !< ocean horizontal grid structure; G or dG is required if the new file

    !! uses any horizontal grid axes.

  type(dyn_horgrid_type), optional, intent(in) :: dg

    !< dynamic horizontal grid structure; G or dG is required if the new file

    !! uses any horizontal grid axes.

  type(verticalgrid_type), optional, intent(in) :: gv

    !< ocean vertical grid structure, which is ! required if the new file uses

    !! any vertical grid axes.

  integer(kind=int64), optional, intent(in) :: checksums(:,:)

    !< checksums of vars

  type(axis_info), optional, intent(in) :: extra_axes(:)

    !< Types with information about some axes that might be used in this file

  type(attribute_info), optional, intent(in) :: global_atts(:)

    !< Global attributes to write to this file

  type(mom_infra_file) :: new_file

  type(mom_field) :: new_fields(novars)

  new_file%handle_infra = io_handle

  call create_mom_file(new_file, filename, vars, novars, new_fields, &

      threading=threading, timeunit=timeunit, g=g, dg=dg, gv=gv, &

      checksums=checksums, extra_axes=extra_axes, global_atts=global_atts)

  io_handle = new_file%handle_infra

  call new_file%get_file_fieldtypes(fields(:novars))

end subroutine create_file

!! Create a new netCDF file and register the MOM_fields to be written.

subroutine create_mom_file(IO_handle, filename, vars, novars, fields, &

    threading, timeunit, G, dG, GV, checksums, extra_axes, global_atts)

  class(mom_file),       intent(inout) :: io_handle  !< Handle for a files or fileset that is to be

                                                     !! opened or reopened for writing

  character(len=*),      intent(in)    :: filename   !< full path to the file to create

  type(vardesc),         intent(in)    :: vars(:)    !< structures describing fields written to filename

  integer,               intent(in)    :: novars     !< number of fields written to filename

  type(mom_field),       intent(inout) :: fields(:)  !< array of fieldtypes for each variable

  integer, optional,     intent(in)    :: threading  !< SINGLE_FILE or MULTIPLE

  real, optional,        intent(in)    :: timeunit   !< length of the units for time [s]. The

                                                     !! default value is 86400.0, for 1 day.

  type(ocean_grid_type),   optional, intent(in) :: g !< ocean horizontal grid structure; G or dG

                                                     !! is required if the new file uses any

                                                     !! horizontal grid axes.

  type(dyn_horgrid_type),  optional, intent(in) :: dg !< dynamic horizontal grid structure; G or dG

                                                     !! is required if the new file uses any

                                                     !! horizontal grid axes.

  type(verticalgrid_type), optional, intent(in) :: gv !< ocean vertical grid structure, which is

                                                     !! required if the new file uses any

                                                     !! vertical grid axes.

  integer(kind=int64),     optional, intent(in) :: checksums(:,:) !< checksums of vars

  type(axis_info),         dimension(:), &

                           optional, intent(in) :: extra_axes !< Types with information about

                                                     !! some axes that might be used in this file

  type(attribute_info),    optional, intent(in) :: global_atts(:) !< Global attributes to

                                                     !! write to this file

  logical        :: use_lath, use_lonh, use_latq, use_lonq, use_time

  logical        :: use_layer, use_int, use_periodic

  logical        :: one_file, domain_set, dim_found

  logical, dimension(:), allocatable :: use_extra_axis

  type(mom_axis) :: axis_lath, axis_latq, axis_lonh, axis_lonq

  type(mom_axis) :: axis_layer, axis_int, axis_time, axis_periodic

  type(mom_axis), dimension(:), allocatable :: more_axes ! Axes generated from extra_axes

  type(mom_axis) :: axes(5)     ! The axes of a variable

  type(mom_domain_type), pointer :: domain => null()

  type(domain1d) :: x_domain, y_domain

  integer        :: position, numaxes, pack, thread, k, n, m

  integer        :: num_extra_dims ! The number of extra possible dimensions from extra_axes

  integer        :: isg, ieg, jsg, jeg, isgb, iegb, jsgb, jegb

  integer        :: var_periods, num_periods=0

  real, dimension(:), allocatable :: axis_val ! Axis label values [various]

  real, pointer, dimension(:) :: &

    gridlatt => null(), & ! The latitude of T or B points for the purpose of labeling

    gridlatb => null(), & ! the output axes, often in units of [degrees_N] or [km] or [m].

    gridlont => null(), & ! The longitude of T or B points for the purpose of labeling

    gridlonb => null()    ! the output axes, often in units of [degrees_E] or [km] or [m].

  character(len=40) :: time_units, x_axis_units, y_axis_units

  character(len=8)  :: t_grid, t_grid_read

  character(len=64) :: ax_name(5)  ! The axis names of a variable

  use_lath  = .false. ; use_lonh     = .false.

  use_latq  = .false. ; use_lonq     = .false.

  use_time  = .false. ; use_periodic = .false.

  use_layer = .false. ; use_int      = .false.

  num_extra_dims = 0

  if (present(extra_axes)) then

    num_extra_dims = size(extra_axes)

    if (num_extra_dims > 0) then

      allocate(use_extra_axis(num_extra_dims)) ; use_extra_axis = .false.

      allocate(more_axes(num_extra_dims))

    endif

  endif

  thread = single_file

  if (PRESENT(threading)) thread = threading

  domain_set = .false.

  if (present(g)) then

    domain_set = .true. ; domain => g%Domain

    gridlatt => g%gridLatT ; gridlatb => g%gridLatB

    gridlont => g%gridLonT ; gridlonb => g%gridLonB

    x_axis_units = g%x_axis_units ; y_axis_units = g%y_axis_units

    isg = g%isg ; ieg = g%ieg ; jsg = g%jsg ; jeg = g%jeg

    isgb = g%IsgB ; iegb = g%IegB ; jsgb = g%JsgB ; jegb = g%JegB

  elseif (present(dg)) then

    domain_set = .true. ; domain => dg%Domain

    gridlatt => dg%gridLatT ; gridlatb => dg%gridLatB

    gridlont => dg%gridLonT ; gridlonb => dg%gridLonB

    x_axis_units = dg%x_axis_units ; y_axis_units = dg%y_axis_units

    isg = dg%isg ; ieg = dg%ieg ; jsg = dg%jsg ; jeg = dg%jeg

    isgb = dg%IsgB ; iegb = dg%IegB ; jsgb = dg%JsgB ; jegb = dg%JegB

  endif

  one_file = .true.

  if (domain_set) one_file = (thread == single_file)

  if (one_file) then

    if (domain_set) then

      call io_handle%open(filename, action=overwrite_file, &

          mom_domain=domain, threading=thread, fileset=single_file)

    else

      call io_handle%open(filename, action=overwrite_file, threading=thread, &

          fileset=single_file)

    endif

  else

    call io_handle%open(filename, action=overwrite_file, mom_domain=domain, &

        threading=thread, fileset=thread)

  endif

! Define the coordinates.

  do k=1,novars

    position = vars(k)%position

    if (position == -1) position = position_from_horgrid(vars(k)%hor_grid)

    select case (position)

      case (center) ; use_lath = .true. ; use_lonh = .true.

      case (corner) ; use_latq = .true. ; use_lonq = .true.

      case (east_face) ; use_lath = .true. ; use_lonq = .true.

      case (north_face) ; use_latq = .true. ; use_lonh = .true.

      case (0) ! Do nothing.

      case default

        call mom_error(warning, "MOM_io create_file: "//trim(vars(k)%name)//" has an unrecognized value of postion")

    end select

    select case (vars(k)%z_grid)

      case ('L') ; use_layer = .true.

      case ('i') ; use_int = .true.

      case ('1') ! Do nothing.

      case default

        call mom_error(fatal, "MOM_io create_file: "//trim(vars(k)%name)//&

                        " has unrecognized z_grid "//trim(vars(k)%z_grid))

    end select

    t_grid = adjustl(vars(k)%t_grid)

    select case (t_grid(1:1))

      case ('s', 'a', 'm') ; use_time = .true.

      case ('p') ; use_periodic = .true.

        if (len_trim(t_grid(2:8)) <= 0) call mom_error(fatal, &

          "MOM_io create_file: No periodic axis length was specified in "//&

          trim(vars(k)%t_grid) // " in the periodic axes of variable "//&

          trim(vars(k)%name)//" in file "//trim(filename))

        var_periods = -9999999

        t_grid_read = adjustl(t_grid(2:8))

        read(t_grid_read,*) var_periods

        if (var_periods == -9999999) call mom_error(fatal, &

          "MOM_io create_file: Failed to read the number of periods from "//&

          trim(vars(k)%t_grid) // " in the periodic axes of variable "//&

          trim(vars(k)%name)//" in file "//trim(filename))

        if (var_periods < 1) call mom_error(fatal, "MOM_io create_file: "//&

           "variable "//trim(vars(k)%name)//" in file "//trim(filename)//&

           " uses a periodic time axis, and must have a positive "//&

           "value for the number of periods in "//vars(k)%t_grid )

        if ((num_periods > 0) .and. (var_periods /= num_periods)) &

          call mom_error(fatal, "MOM_io create_file: "//&

            "Only one value of the number of periods can be used in the "//&

            "create_file call for file "//trim(filename)//".  The second is "//&

            "variable "//trim(vars(k)%name)//" with t_grid "//vars(k)%t_grid )

        num_periods = var_periods

      case ('1') ! Do nothing.

      case default

        call mom_error(warning, "MOM_io create_file: "//trim(vars(k)%name)//&

                        " has unrecognized t_grid "//trim(vars(k)%t_grid))

    end select

    do n=1,5 ; if (len_trim(vars(k)%dim_names(n)) > 0) then

      dim_found = .false.

      do m=1,num_extra_dims

        if (lowercase(trim(vars(k)%dim_names(n))) == lowercase(trim(extra_axes(m)%name))) then

          use_extra_axis(m) = .true.

          dim_found = .true.

          exit

        endif

      enddo

      if (.not.dim_found) call mom_error(fatal, "Unable to find a match for dimension "//&

          trim(vars(k)%dim_names(n))//" for variable "//trim(vars(k)%name)//" in file "//trim(filename))

    endif ; enddo

  enddo

  if ((use_lath .or. use_lonh .or. use_latq .or. use_lonq)) then

    if (.not.domain_set) call mom_error(fatal, "create_file: "//&

      "An ocean_grid_type or dyn_horgrid_type is required to create a file with a horizontal coordinate.")

    call get_domain_components(domain, x_domain, y_domain)

  endif

  if ((use_layer .or. use_int) .and. .not.present(gv)) call mom_error(fatal, &

    "create_file: A vertical grid type is required to create a file with a vertical coordinate.")

  if (use_lath) &

    axis_lath = io_handle%register_axis("lath", units=y_axis_units, longname="Latitude", &

                        cartesian='Y', domain=y_domain, data=gridlatt(jsg:jeg))

  if (use_lonh) &

    axis_lonh = io_handle%register_axis("lonh", units=x_axis_units, longname="Longitude", &

                        cartesian='X', domain=x_domain, data=gridlont(isg:ieg))

  if (use_latq) &

    axis_latq = io_handle%register_axis("latq", units=y_axis_units, longname="Latitude", &

                        cartesian='Y', domain=y_domain, data=gridlatb(jsgb:jegb), edge_axis=.true.)

  if (use_lonq) &

    axis_lonq = io_handle%register_axis("lonq", units=x_axis_units, longname="Longitude", &

                        cartesian='X', domain=x_domain, data=gridlonb(isgb:iegb), edge_axis=.true.)

  if (use_layer) &

    axis_layer = io_handle%register_axis("Layer", units=trim(gv%zAxisUnits), &

                        longname="Layer "//trim(gv%zAxisLongName), cartesian='Z', &

                        sense=1, data=gv%sLayer(1:gv%ke))

  if (use_int) &

    axis_int = io_handle%register_axis("Interface", units=trim(gv%zAxisUnits), &

                        longname="Interface "//trim(gv%zAxisLongName), cartesian='Z', &

                        sense=1, data=gv%sInterface(1:gv%ke+1))

  if (use_time) then ; if (present(timeunit)) then

    ! Set appropriate units, depending on the value.

    if (timeunit < 0.0) then

      time_units = "days" ! The default value.

    elseif ((timeunit >= 0.99) .and. (timeunit < 1.01)) then

      time_units = "seconds"

    elseif ((timeunit >= 3599.0) .and. (timeunit < 3601.0)) then

      time_units = "hours"

    elseif ((timeunit >= 86399.0) .and. (timeunit < 86401.0)) then

      time_units = "days"

    elseif ((timeunit >= 3.0e7) .and. (timeunit < 3.2e7)) then

      time_units = "years"

    else

      write(time_units,'(es8.2," s")') timeunit

    endif

    axis_time = io_handle%register_axis("Time", units=time_units, longname="Time", cartesian='T')

  else

    axis_time = io_handle%register_axis("Time", units="days", longname="Time", cartesian='T')

  endif ; endif

  if (use_periodic) then

    if (num_periods <= 1) call mom_error(fatal, "MOM_io create_file: "//&

      "num_periods for file "//trim(filename)//" must be at least 1.")

    ! Define a periodic axis with unit labels.

    allocate(axis_val(num_periods))

    do k=1,num_periods ; axis_val(k) = real(k) ; enddo

    axis_periodic = io_handle%register_axis("Period", units="nondimensional", &

        longname="Periods for cyclical variables", cartesian='T', data=axis_val)

    deallocate(axis_val)

  endif

  do m=1,num_extra_dims ; if (use_extra_axis(m)) then

    if (allocated(extra_axes(m)%ax_data)) then

      more_axes(m) = io_handle%register_axis(extra_axes(m)%name, units=extra_axes(m)%units, &

                          longname=extra_axes(m)%longname, cartesian=extra_axes(m)%cartesian, &

                          sense=extra_axes(m)%sense, data=extra_axes(m)%ax_data)

    elseif (trim(extra_axes(m)%cartesian) == "T") then

      more_axes(m) = io_handle%register_axis(extra_axes(m)%name, units=extra_axes(m)%units, &

                          longname=extra_axes(m)%longname, cartesian=extra_axes(m)%cartesian)

    else

      ! FMS requires that non-time axes have variables that label their values, even if they are trivial.

      allocate (axis_val(extra_axes(m)%ax_size))

      do k=1,extra_axes(m)%ax_size ; axis_val(k) = real(k) ; enddo

      more_axes(m) = io_handle%register_axis(extra_axes(m)%name, units=extra_axes(m)%units, &

                          longname=extra_axes(m)%longname, cartesian=extra_axes(m)%cartesian, &

                          sense=extra_axes(m)%sense, data=axis_val)

      deallocate(axis_val)

    endif

  endif ; enddo

  do k=1,novars

    numaxes = 0

    position = vars(k)%position

    if (position == -1) position = position_from_horgrid(vars(k)%hor_grid)

    select case (position)

      case (center)

        numaxes = 2 ; axes(1) = axis_lonh ; axes(2) = axis_lath ; ax_name(1) = "lonh" ; ax_name(2) = "lath"

      case (corner)

        numaxes = 2 ; axes(1) = axis_lonq ; axes(2) = axis_latq ; ax_name(1) = "lonq" ; ax_name(2) = "latq"

      case (east_face)

        numaxes = 2 ; axes(1) = axis_lonq ; axes(2) = axis_lath ; ax_name(1) = "lonq" ; ax_name(2) = "lath"

      case (north_face)

        numaxes = 2 ; axes(1) = axis_lonh ; axes(2) = axis_latq ; ax_name(1) = "lonh" ; ax_name(2) = "latq"

      case (0) ! Do nothing.

      case default

        call mom_error(warning, "MOM_io create_file: "//trim(vars(k)%name)//&

                        " has unrecognized position, hor_grid = "//trim(vars(k)%hor_grid))

    end select

    select case (vars(k)%z_grid)

      case ('L') ; numaxes = numaxes+1 ; axes(numaxes) = axis_layer ; ax_name(numaxes) = "Layer"

      case ('i') ; numaxes = numaxes+1 ; axes(numaxes) = axis_int ; ax_name(numaxes) = "Interface"

      case ('1') ! Do nothing.

      case default

        call mom_error(fatal, "MOM_io create_file: "//trim(vars(k)%name)//&

                        " has unrecognized z_grid "//trim(vars(k)%z_grid))

    end select

    do n=1,numaxes

      if ( (len_trim(vars(k)%dim_names(n)) > 0) .and. (trim(ax_name(n)) /= trim(vars(k)%dim_names(n))) ) &

        call mom_error(warning, "MOM_io create_file: dimension "//trim(ax_name(n))//&

               " of variable "//trim(vars(k)%name)//" in "//trim(filename)//&

               " is being set inconsistently as "//trim(vars(k)%dim_names(n)))

    enddo

    do n=numaxes+1,5 ; if (len_trim(vars(k)%dim_names(n)) > 0) then

      dim_found = .false.

      do m=1,num_extra_dims

        if (lowercase(trim(vars(k)%dim_names(n))) == lowercase(trim(extra_axes(m)%name))) then

          numaxes = numaxes+1 ; axes(numaxes) = more_axes(m)

          exit

        endif

      enddo

    endif ; enddo

    t_grid = adjustl(vars(k)%t_grid)

    select case (t_grid(1:1))

      case ('s', 'a', 'm') ; numaxes = numaxes+1 ; axes(numaxes) = axis_time

      case ('p')           ; numaxes = numaxes+1 ; axes(numaxes) = axis_periodic

      case ('1') ! Do nothing.

      case default

        call mom_error(warning, "MOM_io create_file: "//trim(vars(k)%name)//&

                        " has unrecognized t_grid "//trim(vars(k)%t_grid))

    end select

    pack = 1

    if (present(checksums)) then

      fields(k) = io_handle%register_field(axes(1:numaxes), vars(k)%name, vars(k)%units, &

                          vars(k)%longname, pack=pack, checksum=checksums(k,:), conversion=vars(k)%conversion)

    else

      fields(k) = io_handle%register_field(axes(1:numaxes), vars(k)%name, vars(k)%units, &

                          vars(k)%longname, pack=pack, conversion=vars(k)%conversion)

    endif

  enddo

  if (present(global_atts)) then

    do n=1,size(global_atts)

      if (allocated(global_atts(n)%name) .and. allocated(global_atts(n)%att_val)) &

        call io_handle%write_attribute(global_atts(n)%name, global_atts(n)%att_val)

    enddo

  endif

  ! Now write the variables with the axis label values

  if (use_lath) call io_handle%write_field(axis_lath)

  if (use_latq) call io_handle%write_field(axis_latq)

  if (use_lonh) call io_handle%write_field(axis_lonh)

  if (use_lonq) call io_handle%write_field(axis_lonq)

  if (use_layer) call io_handle%write_field(axis_layer)

  if (use_int) call io_handle%write_field(axis_int)

  if (use_periodic) call io_handle%write_field(axis_periodic)

  do m=1,num_extra_dims ; if (use_extra_axis(m)) then

    call io_handle%write_field(more_axes(m))

  endif ; enddo

  if (num_extra_dims > 0) then

    deallocate(use_extra_axis, more_axes)

  endif

end subroutine create_mom_file

!> `reopen_MOM_file` wrapper for the legacy file handle, `file_type`.

!! NOTE: This function may be removed in a future release.

subroutine reopen_file(IO_handle, filename, vars, novars, fields, threading, &

    timeunit, G, dG, GV, extra_axes, global_atts)

  type(file_type), intent(inout) :: io_handle

    !< Handle for a file or fileset that is to be opened or reopened for

    !! writing

  character(len=*), intent(in) :: filename

    !< full path to the file to create

  type(vardesc), intent(in) :: vars(:)

    !< structures describing fields written to filename

  integer, intent(in) :: novars

    !< number of fields written to filename

  type(fieldtype), intent(inout) :: fields(:)

    !< array of fieldtypes for each variable

  integer, optional, intent(in) :: threading

    !< SINGLE_FILE or MULTIPLE

  real, optional, intent(in) :: timeunit

    !< length of the units for time [s]. The default value is 86400.0, for 1

    !! day.

  type(ocean_grid_type), optional, intent(in) :: g

    !< ocean horizontal grid structure; G or dG is required if a new file uses

    !! any horizontal grid axes.

  type(dyn_horgrid_type), optional, intent(in) :: dg

    !< dynamic horizontal grid structure; G or dG is required if a new file

    !! uses any horizontal grid axes.

  type(verticalgrid_type), optional, intent(in) :: gv

    !< ocean vertical grid structure, which is required if a new file uses any

    !! vertical grid axes.

  type(axis_info), optional, intent(in) :: extra_axes(:)

    !< Types with information about some axes that might be used in this file

  type(attribute_info), optional, intent(in) :: global_atts(:)

    !< Global attributes to write to this file

  type(mom_infra_file) :: mfile

    !< Wrapper to MOM file

  type(mom_field), allocatable :: mfields(:)

    !< Wrapper to MOM fields

  integer :: i

  mfile%handle_infra = io_handle

  allocate(mfields(size(fields)))

  call reopen_mom_file(mfile, filename, vars, novars, mfields, &

      threading=threading, timeunit=timeunit, g=g, dg=dg, gv=gv, &

      extra_axes=extra_axes, global_atts=global_atts)

  io_handle = mfile%handle_infra

  call get_file_fields(io_handle, fields)

end subroutine reopen_file

!> This routine opens an existing NetCDF file for output.  If it

!! does not find the file, a new file is created.  It also sets up

!! structures that describe this file and the variables that will

!! later be written to this file.

subroutine reopen_mom_file(IO_handle, filename, vars, novars, fields, &

    threading, timeunit, G, dG, GV, extra_axes, global_atts)

  class(mom_file),       intent(inout) :: io_handle  !< Handle for a file or fileset that is to be

                                                     !! opened or reopened for writing

  character(len=*),      intent(in)    :: filename   !< full path to the file to create

  type(vardesc),         intent(in)    :: vars(:)    !< structures describing fields written to filename

  integer,               intent(in)    :: novars     !< number of fields written to filename

  type(mom_field),       intent(inout) :: fields(:)  !< array of fieldtypes for each variable

  integer, optional,     intent(in)    :: threading  !< SINGLE_FILE or MULTIPLE

  real, optional,        intent(in)    :: timeunit   !< length of the units for time [s]. The

                                                     !! default value is 86400.0, for 1 day.

  type(ocean_grid_type),   optional, intent(in) :: g !< ocean horizontal grid structure; G or dG

                                                     !! is required if a new file uses any

                                                     !! horizontal grid axes.

  type(dyn_horgrid_type),  optional, intent(in) :: dg !< dynamic horizontal grid structure; G or dG

                                                     !! is required if a new file uses any

                                                     !! horizontal grid axes.

  type(verticalgrid_type), optional, intent(in) :: gv !< ocean vertical grid structure, which is

                                                     !! required if a new file uses any

                                                     !! vertical grid axes.

  type(axis_info),         optional, intent(in) :: extra_axes(:)  !< Types with information about

                                                     !! some axes that might be used in this file

  type(attribute_info),    optional, intent(in) :: global_atts(:) !< Global attributes to

                                                     !! write to this file

  type(mom_domain_type), pointer :: domain => null()

  character(len=200) :: check_name, mesg

  integer :: length, nvar, thread

  logical :: exists, one_file, domain_set

  thread = single_file

  if (PRESENT(threading)) thread = threading

  ! For single-file IO, only the root PE is required to set up the fields.

  ! This permits calls by either the root PE or all PEs

  if (.not. is_root_pe() .and. thread == single_file) return

  ! For multiple IO domains, we would need additional functionality:

  ! * Identify ranks as IO PEs

  ! * Determine the filename of

  ! Neither of these tasks should be handed by MOM6, so we cannot safely use

  ! this function.  A framework-specific `inquire()` function is needed.

  ! Until it exists, we will disable this function.

  if (thread == multiple) &

    call mom_error(fatal, 'reopen_MOM_file does not yet support files with ' &

      // 'multiple I/O domains.')

  check_name = filename

  length = len(trim(check_name))

  if (check_name(length-2:length) /= ".nc") check_name = trim(check_name)//".nc"

  if (thread /= single_file) check_name = trim(check_name)//".0000"

  inquire(file=check_name,exist=exists)

  if (.not.exists) then

    call create_mom_file(io_handle, filename, vars, novars, fields, &

        threading, timeunit, g=g, dg=dg, gv=gv, extra_axes=extra_axes, &

        global_atts=global_atts)

  else

    domain_set = .false.

    if (present(g)) then

      domain_set = .true. ; domain => g%Domain

    elseif (present(dg)) then

      domain_set = .true. ; domain => dg%Domain

    endif

    one_file = .true.

    if (domain_set) one_file = (thread == single_file)

    if (one_file) then

      call io_handle%open(filename, append_file, threading=thread)

    else

      call io_handle%open(filename, append_file, mom_domain=domain)

    endif

    if (.not. io_handle%file_is_open()) return

    call io_handle%get_file_info(nvar=nvar)

    if (nvar == -1) then

      write (mesg,*) "Reopening file ",trim(filename)," apparently had ",nvar,&

                     " variables. Clobbering and creating file with ",novars," instead."

      call mom_error(warning,"MOM_io: "//mesg)

      call create_mom_file(io_handle, filename, vars, novars, fields, &

          threading, timeunit, g=g, dg=dg, gv=gv, extra_axes=extra_axes, &

          global_atts=global_atts)

    elseif (nvar /= novars) then

      write (mesg,*) "Reopening file ",trim(filename)," with ",novars,&

                     " variables instead of ",nvar,"."

      call mom_error(fatal,"MOM_io: "//mesg)

    endif

    if (nvar > 0) call io_handle%get_file_fields(fields(1:nvar))

  endif

end subroutine reopen_mom_file

!> Return the index of sdtout if called from the root PE, or 0 for other PEs.

integer function stdout_if_root()

  stdout_if_root = 0

  if (is_root_pe()) stdout_if_root = stdout

end function stdout_if_root

!> This function determines how many time levels a variable has in a file.

function num_timelevels(filename, varname, min_dims) result(n_time)

  character(len=*),  intent(in) :: filename   !< name of the file to read

  character(len=*),  intent(in) :: varname    !< variable whose number of time levels

                                              !! are to be returned

  integer, optional, intent(in) :: min_dims   !< The minimum number of dimensions a variable must have

                                              !! if it has a time dimension.  If the variable has 1 less

                                              !! dimension than this, then 0 is returned.

  integer :: n_time                           !< number of time levels varname has in filename

  character(len=256) :: msg

  integer :: ndims

  integer :: sizes(8)

  n_time = -1

  ! To do almost the same via MOM_io_infra calls, we could do the following:

  !   found = field_exists(filename, varname)

  !   if (found) then

  !     call open_file(ncid, filename, action=READONLY_FILE, form=NETCDF_FILE, threading=MULTIPLE)

  !     call get_file_info(ncid, ntime=n_time)

  !   endif

  ! However, this does not handle the case where the time axis for the variable is not the record

  ! axis and min_dims is not used.

  call get_var_sizes(filename, varname, ndims, sizes, match_case=.false., caller="num_timelevels")

  if (ndims > 0) n_time = sizes(ndims)

  if (present(min_dims)) then

    if (ndims < min_dims-1) then

      write(msg, '(I0)') min_dims

      call mom_error(warning, "num_timelevels: variable "//trim(varname)//" in file "//&

          trim(filename)//" has fewer than min_dims = "//trim(msg)//" dimensions.")

      n_time = -1

    elseif (ndims == min_dims - 1) then

      n_time = 0

    endif

  endif

end function num_timelevels

!> get_var_sizes returns the number and size of dimensions associate with a variable in a file.

!! Usually only the root PE does the read, and then the information is broadcast

subroutine get_var_sizes(filename, varname, ndims, sizes, match_case, caller, all_read, dim_names, ncid_in)

  character(len=*),      intent(in)  :: filename   !< Name of the file to read, used here in messages

  character(len=*),      intent(in)  :: varname    !< The variable name, used here for messages

  integer,               intent(out) :: ndims      !< The number of dimensions to the variable

  integer, dimension(:), intent(out) :: sizes      !< The dimension sizes, or 0 for extra values

  logical,     optional, intent(in)  :: match_case !< If false, allow for variables name matches to be

                                                   !! case insensitive, but take a perfect match if

                                                   !! found.  The default is true.

  character(len=*), optional, intent(in) :: caller !< The name of a calling routine for use in error messages

  logical,     optional, intent(in)  :: all_read   !< If present and true, all PEs that call this

                                                   !! routine actually do the read, otherwise only

                                                   !! root PE reads and then it broadcasts the results.

  character(len=*), dimension(:), &

               optional, intent(out) :: dim_names  !< The names of the dimensions for this variable

  integer,     optional, intent(in)  :: ncid_in    !< The netCDF ID of an open file.  If absent, the

                                                   !! file is opened and closed within this routine.

  logical :: do_read, do_broadcast

  integer, allocatable :: size_msg(:)  ! An array combining the number of dimensions and the sizes.

  integer :: n, nval

  do_read = is_root_pe()

  if (present(all_read)) do_read = all_read .or. do_read

  do_broadcast = .true. ; if (present(all_read)) do_broadcast = .not.all_read

  if (do_read) call read_var_sizes(filename, varname, ndims, sizes, match_case, caller, dim_names, ncid_in)

  if (do_broadcast) then

    ! Distribute the sizes from the root PE.

    nval = size(sizes) + 1

    allocate(size_msg(nval))

    size_msg(1) = ndims

    do n=2,nval ; size_msg(n) = sizes(n-1) ; enddo

    call broadcast(size_msg, nval, blocking=.true.)

    ndims = size_msg(1)

    do n=2,nval ;  sizes(n-1) = size_msg(n) ; enddo

    deallocate(size_msg)

    if (present(dim_names) .and. (ndims > 0)) then

      nval = min(ndims, size(dim_names))

      call broadcast(dim_names(1:nval), len(dim_names(1)), blocking=.true.)

    endif

  endif

end subroutine get_var_sizes

!> read_var_sizes returns the number and size of dimensions associated with a variable in a file.

!! If the variable is not in the file the returned sizes are all 0 and ndims is -1.

!! Every processor for which this is called does the reading.

subroutine read_var_sizes(filename, varname, ndims, sizes, match_case, caller, dim_names, ncid_in)

  character(len=*),      intent(in)  :: filename   !< Name of the file to read, used here in messages

  character(len=*),      intent(in)  :: varname    !< The variable name, used here for messages

  integer,               intent(out) :: ndims      !< The number of dimensions to the variable

  integer, dimension(:), intent(out) :: sizes      !< The dimension sizes, or 0 for extra values

  logical,     optional, intent(in)  :: match_case !< If false, allow for variables name matches to be

                                                   !! case insensitive, but take a perfect match if

                                                   !! found.  The default is true.

  character(len=*), &

               optional, intent(in)  :: caller     !< The name of a calling routine for use in error messages

  character(len=*), dimension(:), &

               optional, intent(out) :: dim_names  !< The names of the dimensions for this variable

  integer,     optional, intent(in)  :: ncid_in    !< The netCDF ID of an open file.  If absent, the

                                                   !! file is opened and closed within this routine.

  character(len=256) :: hdr, dimname

  integer, allocatable :: dimids(:)

  integer :: varid, ncid, n, status

  logical :: success, found

  hdr = "get_var_size: " ; if (present(caller)) hdr = trim(hdr)//": "

  sizes(:) = 0 ; ndims = -1

  if (present(ncid_in)) then

    ncid = ncid_in

  else

    call open_file_to_read(filename, ncid, success=success)

    if (.not.success) then

      call mom_error(warning, "Unsuccessfully attempted to open file "//trim(filename))

      return

    endif

  endif

  ! Get the dimension sizes of the variable varname.

  call get_varid(varname, ncid, filename, varid, match_case=match_case, found=found)

  if (.not.found) then

    call mom_error(warning, "Could not find variable "//trim(varname)//" in file "//trim(filename))

    return

  endif

  status = nf90_inquire_variable(ncid, varid, ndims=ndims)

  if (status /= nf90_noerr) then

    call mom_error(warning, trim(hdr) // trim(nf90_strerror(status)) //&

      " Getting number of dimensions of "//trim(varname)//" in "//trim(filename))

    return

  endif

  if (ndims < 1) return

  allocate(dimids(ndims))

  status = nf90_inquire_variable(ncid, varid, dimids=dimids(1:ndims))

  if (status /= nf90_noerr) then

    call mom_error(warning, trim(hdr) // trim(nf90_strerror(status)) //&

      " Getting dimension IDs for "//trim(varname)//" in "//trim(filename))

    deallocate(dimids) ; return

  endif

  do n = 1, min(ndims,size(sizes))

    status = nf90_inquire_dimension(ncid, dimids(n), name=dimname, len=sizes(n))

    if (status /= nf90_noerr) call mom_error(warning, trim(hdr) // trim(nf90_strerror(status)) //&

        " Getting dimension length for "//trim(varname)//" in "//trim(filename))

    if (present(dim_names)) then

      if (n <= size(dim_names)) dim_names(n) = trim(dimname)

    endif

  enddo

  deallocate(dimids)

  if (.not.present(ncid_in)) call close_file_to_read(ncid, filename)

end subroutine read_var_sizes

!> Read a real scalar variable from a netCDF file with the root PE, and broadcast the

!! results to all the other PEs.

subroutine read_variable_0d(filename, varname, var, ncid_in, scale)

  character(len=*),  intent(in)    :: filename !< The name of the file to read

  character(len=*),  intent(in)    :: varname  !< The variable name of the data in the file

  real,              intent(inout) :: var      !< The scalar into which to read the data in arbitrary units [A ~> a]

  integer, optional, intent(in)    :: ncid_in  !< The netCDF ID of an open file.  If absent, the

                                               !! file is opened and closed within this routine

  real,    optional, intent(in)    :: scale    !< A scaling factor that the variable is multiplied by

                                               !! before it is returned to convert from the units in the file

                                               !! to the internal units for this variable [A a-1 ~> 1]

  integer :: varid, ncid, rc

  character(len=256) :: hdr

  hdr = "read_variable_0d"

  if (is_root_pe()) then

    if (present(ncid_in)) then

      ncid = ncid_in

    else

      call open_file_to_read(filename, ncid)

    endif

    call get_varid(varname, ncid, filename, varid, match_case=.false.)

    if (varid < 0) call mom_error(fatal, "Unable to get netCDF varid for "//trim(varname)//&

                                         " in "//trim(filename))

    rc = nf90_get_var(ncid, varid, var)

    if (rc /= nf90_noerr) call mom_error(fatal, trim(hdr) // trim(nf90_strerror(rc)) //&

          " Difficulties reading "//trim(varname)//" from "//trim(filename))

    if (.not.present(ncid_in)) call close_file_to_read(ncid, filename)

    if (present(scale)) var = scale * var

  endif

  call broadcast(var, blocking=.true.)

end subroutine read_variable_0d

!> Read a 1-d real variable from a netCDF file with the root PE, and broadcast the

!! results to all the other PEs.

subroutine read_variable_1d(filename, varname, var, ncid_in, scale)

  character(len=*),   intent(in)    :: filename !< The name of the file to read

  character(len=*),   intent(in)    :: varname  !< The variable name of the data in the file

  real, dimension(:), intent(inout) :: var      !< The 1-d array into which to read the data in arbitrary units [A ~> a]

  integer,  optional, intent(in)    :: ncid_in  !< The netCDF ID of an open file.  If absent, the

                                                !! file is opened and closed within this routine

  real,     optional, intent(in)    :: scale    !< A scaling factor that the variable is multiplied by

                                                !! before it is returned to convert from the units in the file

                                                !! to the internal units for this variable [A a-1 ~> 1]

  integer :: varid, ncid, rc

  character(len=256) :: hdr

  hdr = "read_variable_1d"

  if (is_root_pe()) then

    if (present(ncid_in)) then

      ncid = ncid_in

    else

      call open_file_to_read(filename, ncid)

    endif

    call get_varid(varname, ncid, filename, varid, match_case=.false.)

    if (varid < 0) call mom_error(fatal, "Unable to get netCDF varid for "//trim(varname)//&

                                         " in "//trim(filename))

    rc = nf90_get_var(ncid, varid, var)

    if (rc /= nf90_noerr) call mom_error(fatal, trim(hdr) // trim(nf90_strerror(rc)) //&

          " Difficulties reading "//trim(varname)//" from "//trim(filename))

    if (.not.present(ncid_in)) call close_file_to_read(ncid, filename)

    if (present(scale)) then ; if (scale /= 1.0) then

      var(:) = scale * var(:)

    endif ; endif

  endif

  call broadcast(var, size(var), blocking=.true.)

end subroutine read_variable_1d

!> Read a integer scalar variable from a netCDF file with the root PE, and broadcast the

!! results to all the other PEs.

subroutine read_variable_0d_int(filename, varname, var, ncid_in)

  character(len=*),  intent(in)    :: filename !< The name of the file to read

  character(len=*),  intent(in)    :: varname  !< The variable name of the data in the file

  integer,           intent(inout) :: var      !< The scalar into which to read the data

  integer, optional, intent(in)    :: ncid_in  !< The netCDF ID of an open file.  If absent, the

                                               !! file is opened and closed within this routine.

  integer :: varid, ncid, rc

  character(len=256) :: hdr

  hdr = "read_variable_0d_int"

  if (is_root_pe()) then

    if (present(ncid_in)) then

      ncid = ncid_in

    else

      call open_file_to_read(filename, ncid)

    endif

    call get_varid(varname, ncid, filename, varid, match_case=.false.)

    if (varid < 0) call mom_error(fatal, "Unable to get netCDF varid for "//trim(varname)//&

                                         " in "//trim(filename))

    rc = nf90_get_var(ncid, varid, var)

    if (rc /= nf90_noerr) call mom_error(fatal, trim(hdr) // trim(nf90_strerror(rc)) //&

          " Difficulties reading "//trim(varname)//" from "//trim(filename))

    if (.not.present(ncid_in)) call close_file_to_read(ncid, filename)

  endif

  call broadcast(var, blocking=.true.)

end subroutine read_variable_0d_int

!> Read a 1-d integer variable from a netCDF file with the root PE, and broadcast the

!! results to all the other PEs.

subroutine read_variable_1d_int(filename, varname, var, ncid_in)

  character(len=*),      intent(in)    :: filename !< The name of the file to read

  character(len=*),      intent(in)    :: varname  !< The variable name of the data in the file

  integer, dimension(:), intent(inout) :: var      !< The 1-d array into which to read the data

  integer, optional,     intent(in)    :: ncid_in  !< The netCDF ID of an open file.  If absent, the

                                                   !! file is opened and closed within this routine.

  integer :: varid, ncid, rc

  character(len=256) :: hdr

  hdr = "read_variable_1d_int"

  if (is_root_pe()) then

    if (present(ncid_in)) then

      ncid = ncid_in

    else

      call open_file_to_read(filename, ncid)

    endif

    call get_varid(varname, ncid, filename, varid, match_case=.false.)

    if (varid < 0) call mom_error(fatal, "Unable to get netCDF varid for "//trim(varname)//&

                                         " in "//trim(filename))

    rc = nf90_get_var(ncid, varid, var)

    if (rc /= nf90_noerr) call mom_error(fatal, trim(hdr) // trim(nf90_strerror(rc)) //&

          " Difficulties reading "//trim(varname)//" from "//trim(filename))

    if (.not.present(ncid_in)) call close_file_to_read(ncid, filename)

  endif

  call broadcast(var, size(var), blocking=.true.)

end subroutine read_variable_1d_int

!> Read a 2d array from a netCDF input file and save to a variable.

!!

!! Start and nread lenths may exceed var rank.  This allows for reading slices

!! of larger arrays.

!!

!! Previous versions of the model required a time axis on IO fields.  This

!! constraint was dropped in later versions.  As a result, versions both with

!! and without a time axis now exist.  In order to support all such versions,

!! we use a reshaped version of start and nread in order to read the variable

!! as it exists in the file.

!!

!! Certain constraints are still applied to start and nread in order to ensure

!! that varname is a valid 2d array, or contains valid 2d slices.

!!

!! I/O occurs only on the root PE, and data is broadcast to other ranks.

!! Due to potentially large memory communication and storage, this subroutine

!! should only be used when domain-decomposition is unavaialable.

subroutine read_variable_2d(filename, varname, var, start, nread, ncid_in)

  character(len=*), intent(in) :: filename  !< Name of file to be read

  character(len=*), intent(in) :: varname   !< Name of variable to be read

  real, intent(out)            :: var(:,:)  !< Output array of variable [arbitrary]

  integer, optional, intent(in) :: start(:) !< Starting index on each axis.

  integer, optional, intent(in) :: nread(:) !< Number of values to be read along each axis

  integer, optional, intent(in) :: ncid_in  !< netCDF ID of an opened file.

              !! If absent, the file is opened and closed within this routine.

  integer :: ncid, varid

  integer :: field_ndims, dim_len

  integer, allocatable :: field_dimids(:), field_shape(:)

  integer, allocatable :: field_start(:), field_nread(:)

  integer :: i, rc

  character(len=*), parameter :: hdr = "read_variable_2d: "

  ! Validate shape of start and nread

  if (present(start)) then

    if (size(start) < 2) &

      call mom_error(fatal, hdr // trim(nf90_strerror(rc)) &

        // " start must have at least two dimensions.")

  endif

  if (present(nread)) then

    if (size(nread) < 2) &

      call mom_error(fatal, hdr // trim(nf90_strerror(rc)) &

        // " nread must have at least two dimensions.")

    if (any(nread(3:) > 1)) &

      call mom_error(fatal, hdr // trim(nf90_strerror(rc)) &

        // " nread may only read a single level in higher dimensions.")

  endif

  ! Since start and nread may be reshaped, we cannot rely on netCDF to ensure

  ! that their lengths are equivalent, and must do it here.

  if (present(start) .and. present(nread)) then

    if (size(start) /= size(nread)) &

      call mom_error(fatal, hdr // trim(nf90_strerror(rc)) &

        // " start and nread must have the same length.")

  endif

  ! Open and read `varname` from `filename`

  if (is_root_pe()) then

    if (present(ncid_in)) then

      ncid = ncid_in

    else

      call open_file_to_read(filename, ncid)

    endif

    call get_varid(varname, ncid, filename, varid, match_case=.false.)

    if (varid < 0) call mom_error(fatal, "Unable to get netCDF varid for "//trim(varname)//&

                                         " in "//trim(filename))

    ! Query for the dimensionality of the input field

    rc = nf90_inquire_variable(ncid, varid, ndims=field_ndims)

    if (rc /= nf90_noerr) call mom_error(fatal, hdr // trim(nf90_strerror(rc)) //&

          ": Difficulties reading "//trim(varname)//" from "//trim(filename))

    ! Confirm that field is at least 2d

    if (field_ndims < 2) &

      call mom_error(fatal, hdr // trim(nf90_strerror(rc)) // " " // &

          trim(varname) // " from " // trim(filename) // " is not a 2D field.")

    ! If start and nread are present, then reshape them to match field dims

    if (present(start) .or. present(nread)) then

      allocate(field_shape(field_ndims))

      allocate(field_dimids(field_ndims))

      rc = nf90_inquire_variable(ncid, varid, dimids=field_dimids)

      if (rc /= nf90_noerr) call mom_error(fatal, hdr // trim(nf90_strerror(rc)) //&

            ": Difficulties reading "//trim(varname)//" from "//trim(filename))

      do i = 1, field_ndims

        rc = nf90_inquire_dimension(ncid, field_dimids(i), len=dim_len)

        if (rc /= nf90_noerr) &

          call mom_error(fatal, hdr // trim(nf90_strerror(rc)) &

              // ": Difficulties reading dimensions from " // trim(filename))

        field_shape(i) = dim_len

      enddo

      ! Reshape start(:) and nreads(:) in case ranks differ

      allocate(field_start(field_ndims))

      field_start(:) = 1

      if (present(start)) then

        dim_len = min(size(start), size(field_start))

        field_start(:dim_len) = start(:dim_len)

      endif

      allocate(field_nread(field_ndims))

      field_nread(:2) = field_shape(:2)

      field_nread(3:) = 1

      if (present(nread)) field_nread(:2) = nread(:2)

      rc = nf90_get_var(ncid, varid, var, field_start, field_nread)

      deallocate(field_start)

      deallocate(field_nread)

      deallocate(field_shape)

      deallocate(field_dimids)

    else

      rc = nf90_get_var(ncid, varid, var)

    endif

    if (rc /= nf90_noerr) call mom_error(fatal, hdr // trim(nf90_strerror(rc)) //&

          " Difficulties reading "//trim(varname)//" from "//trim(filename))

    if (.not.present(ncid_in)) call close_file_to_read(ncid, filename)

  endif

  call broadcast(var, size(var), blocking=.true.)

end subroutine read_variable_2d

subroutine read_variable_3d(filename, varname, var, start, nread, ncid_in)

  character(len=*), intent(in) :: filename  !< Name of file to be read

  character(len=*), intent(in) :: varname   !< Name of variable to be read

  real, intent(out)            :: var(:,:,:)  !< Output array of variable [arbitrary]

  integer, optional, intent(in) :: start(:) !< Starting index on each axis.

  integer, optional, intent(in) :: nread(:) !< Number of values to be read along each axis

  integer, optional, intent(in) :: ncid_in  !< netCDF ID of an opened file.

              !! If absent, the file is opened and closed within this routine.

  integer :: ncid, varid

  integer :: field_ndims, dim_len

  integer, allocatable :: field_dimids(:), field_shape(:)

  integer, allocatable :: field_start(:), field_nread(:)

  integer :: i, rc

  character(len=*), parameter :: hdr = "read_variable_3d: "

  ! Validate shape of start and nread

  if (present(start)) then

    if (size(start) < 2) &

      call mom_error(fatal, hdr // trim(nf90_strerror(rc)) &

        // " start must have at least two dimensions.")

  endif

  if (present(nread)) then

    if (size(nread) < 2) &

      call mom_error(fatal, hdr // trim(nf90_strerror(rc)) &

        // " nread must have at least two dimensions.")

    if (any(nread(3:) > 1)) &

      call mom_error(fatal, hdr // trim(nf90_strerror(rc)) &

        // " nread may only read a single level in higher dimensions.")

  endif

  ! Since start and nread may be reshaped, we cannot rely on netCDF to ensure

  ! that their lengths are equivalent, and must do it here.

  if (present(start) .and. present(nread)) then

    if (size(start) /= size(nread)) &

      call mom_error(fatal, hdr // trim(nf90_strerror(rc)) &

        // " start and nread must have the same length.")

  endif

  ! Open and read `varname` from `filename`

  if (is_root_pe()) then

    if (present(ncid_in)) then

      ncid = ncid_in

    else

      call open_file_to_read(filename, ncid)

    endif

    call get_varid(varname, ncid, filename, varid, match_case=.false.)

    if (varid < 0) call mom_error(fatal, "Unable to get netCDF varid for "//trim(varname)//&

                                         " in "//trim(filename))

    ! Query for the dimensionality of the input field

    rc = nf90_inquire_variable(ncid, varid, ndims=field_ndims)

    if (rc /= nf90_noerr) call mom_error(fatal, hdr // trim(nf90_strerror(rc)) //&

          ": Difficulties reading "//trim(varname)//" from "//trim(filename))

    ! Confirm that field is at least 2d

    if (field_ndims < 2) &

      call mom_error(fatal, hdr // trim(nf90_strerror(rc)) // " " // &

          trim(varname) // " from " // trim(filename) // " is not a 2D field.")

    ! If start and nread are present, then reshape them to match field dims

    if (present(start) .or. present(nread)) then

      allocate(field_shape(field_ndims))

      allocate(field_dimids(field_ndims))

      rc = nf90_inquire_variable(ncid, varid, dimids=field_dimids)

      if (rc /= nf90_noerr) call mom_error(fatal, hdr // trim(nf90_strerror(rc)) //&

            ": Difficulties reading "//trim(varname)//" from "//trim(filename))

      do i = 1, field_ndims

        rc = nf90_inquire_dimension(ncid, field_dimids(i), len=dim_len)

        if (rc /= nf90_noerr) &

          call mom_error(fatal, hdr // trim(nf90_strerror(rc)) &

              // ": Difficulties reading dimensions from " // trim(filename))

        field_shape(i) = dim_len

      enddo

      ! Reshape start(:) and nreads(:) in case ranks differ

      allocate(field_start(field_ndims))

      field_start(:) = 1

      if (present(start)) then

        dim_len = min(size(start), size(field_start))

        field_start(:dim_len) = start(:dim_len)

      endif

      allocate(field_nread(field_ndims))

      field_nread(:3) = field_shape(:3)

      !field_nread(3:) = 1

      if (present(nread)) field_nread(:3) = nread(:3)

      rc = nf90_get_var(ncid, varid, var, field_start, field_nread)

      deallocate(field_start)

      deallocate(field_nread)

      deallocate(field_shape)

      deallocate(field_dimids)

    else

      rc = nf90_get_var(ncid, varid, var)

    endif

    if (rc /= nf90_noerr) call mom_error(fatal, hdr // trim(nf90_strerror(rc)) //&

          " Difficulties reading "//trim(varname)//" from "//trim(filename))

    if (.not.present(ncid_in)) call close_file_to_read(ncid, filename)

  endif

  call broadcast(var, size(var), blocking=.true.)

end subroutine read_variable_3d

!> Read a character-string global or variable attribute

subroutine read_attribute_str(filename, attname, att_val, varname, found, all_read, ncid_in)

  character(len=*),           intent(in)  :: filename !< Name of the file to read

  character(len=*),           intent(in)  :: attname  !< Name of the attribute to read

  character(:), allocatable,  intent(out) :: att_val  !< The value of the attribute

  character(len=*), optional, intent(in)  :: varname  !< The name of the variable whose attribute will

                                                      !! be read. If missing, read a global attribute.

  logical,          optional, intent(out) :: found    !< Returns true if the attribute is found

  logical,          optional, intent(in)  :: all_read !< If present and true, all PEs that call this

                                                      !! routine actually do the read, otherwise only

                                                      !! root PE reads and then broadcasts the results.

  integer,          optional, intent(in)  :: ncid_in  !< The netCDF ID of an open file.  If absent, the

                                                      !! file is opened and closed within this routine.

  logical :: do_read, do_broadcast

  integer :: rc, ncid, varid, att_type, att_len, info(2)

  character(len=256) :: hdr, att_str

  character(len=:), dimension(:), allocatable :: tmp_str

  hdr = "read_attribute_str"

  att_len = 0

  do_read = is_root_pe() ; if (present(all_read)) do_read = all_read .or. do_read

  do_broadcast = .true. ; if (present(all_read)) do_broadcast = .not.all_read

  if (do_read) then

    if (present(ncid_in)) then

      ncid = ncid_in

    else

      call open_file_to_read(filename, ncid, success=found)

      if (present(found)) then ; if (.not.found) do_read = .false. ; endif

    endif

  endif

  if (do_read) then

    rc = nf90_enotatt ; att_len = 0

    if (present(varname)) then  ! Read a variable attribute

      call get_varid(varname, ncid, filename, varid, match_case=.false., found=found)

      att_str = "att "//trim(attname)//" for "//trim(varname)//" from "//trim(filename)

    else   ! Read a global attribute

      varid = nf90_global

      att_str = "global att "//trim(attname)//" from "//trim(filename)

    endif

    if ((varid > 0) .or. (varid == nf90_global)) then ! The named variable does exist, and found would be true.

      rc = nf90_inquire_attribute(ncid, varid, attname, xtype=att_type, len=att_len)

      if ((.not. present(found)) .or. (rc /= nf90_enotatt)) then

        if ((rc /= nf90_noerr) .and. (rc /= nf90_enotatt)) &

          call mom_error(fatal, trim(hdr) // trim(nf90_strerror(rc)) //" Error getting info for "//trim(att_str))

        if (att_type /= nf90_char) &

          call mom_error(fatal, trim(hdr)//": Attribute data type is not a char for "//trim(att_str))

  !      if (att_len > len(att_val)) &

  !        call MOM_error(FATAL, trim(hdr)//": Insufficiently long string passed in to read "//trim(att_str))

        allocate(character(att_len) :: att_val)

        if (rc == nf90_noerr) then

          rc = nf90_get_att(ncid, varid, attname, att_val)

          if ((rc /= nf90_noerr) .and. (rc /= nf90_enotatt)) &

            call mom_error(fatal, trim(hdr) // trim(nf90_strerror(rc)) //" Difficulties reading "//trim(att_str))

        endif

      endif

    endif

    if (present(found)) found = (rc == nf90_noerr)

    if (.not.present(ncid_in)) call close_file_to_read(ncid, filename)

  endif

  if (do_broadcast) then

    ! Communicate the string length

    info(1) = att_len ; info(2) = 0 ; if (do_read .and. found) info(2) = 1

    call broadcast(info, 2, blocking=.true.)

    if (present(found)) then

      found = (info(2) /= 0)

      if (.not. found) return

    endif

    att_len = info(1)

    if (att_len > 0) then

      ! These extra copies are here because broadcast only supports arrays of strings.

      allocate(character(att_len) :: tmp_str(1))

      if (.not.do_read) allocate(character(att_len) :: att_val)

      if (do_read) tmp_str(1) = att_val

      call broadcast(tmp_str, att_len, blocking=.true.)

      att_val = tmp_str(1)

    elseif (.not.allocated(att_val)) then

      allocate(character(4) :: att_val) ; att_val = ''

    endif

  elseif (.not.allocated(att_val)) then

    allocate(character(4) :: att_val) ; att_val = ''

  endif

end subroutine read_attribute_str

!> Read a 32-bit integer global or variable attribute

subroutine read_attribute_int32(filename, attname, att_val, varname, found, all_read, ncid_in)

  character(len=*),           intent(in)  :: filename !< Name of the file to read

  character(len=*),           intent(in)  :: attname  !< Name of the attribute to read

  integer(kind=int32),        intent(out) :: att_val  !< The value of the attribute

  character(len=*), optional, intent(in)  :: varname  !< The name of the variable whose attribute will

                                                      !! be read. If missing, read a global attribute.

  logical,          optional, intent(out) :: found    !< Returns true if the attribute is found

  logical,          optional, intent(in)  :: all_read !< If present and true, all PEs that call this

                                                      !! routine actually do the read, otherwise only

                                                      !! root PE reads and then broadcasts the results.

  integer,        optional, intent(in)    :: ncid_in  !< The netCDF ID of an open file.  If absent, the

                                                      !! file is opened and closed within this routine.

  logical :: do_read, do_broadcast

  integer :: rc, ncid, varid, is_found

  character(len=256) :: hdr

  hdr = "read_attribute_int32"

  att_val = 0

  do_read = is_root_pe() ; if (present(all_read)) do_read = all_read .or. do_read

  do_broadcast = .true. ; if (present(all_read)) do_broadcast = .not.all_read

  if (do_read) then

    if (present(ncid_in)) then

      ncid = ncid_in

    else

      call open_file_to_read(filename, ncid, success=found)

      if (present(found)) then ; if (.not.found) do_read = .false. ; endif

    endif

  endif

  if (do_read) then

    rc = nf90_enotatt

    if (present(varname)) then  ! Read a variable attribute

      call get_varid(varname, ncid, filename, varid, match_case=.false., found=found)

      if (varid >= 0) then ! The named variable does exist, and found would be true.

        rc = nf90_get_att(ncid, varid, attname, att_val)

        if ((rc /= nf90_noerr) .and. (rc /= nf90_enotatt)) &

          call mom_error(fatal, trim(hdr) // trim(nf90_strerror(rc)) //" Difficulties reading att "//&

                trim(attname)//" for "//trim(varname)//" from "//trim(filename))

      endif

    else  ! Read a global attribute

      rc = nf90_get_att(ncid, nf90_global, attname, att_val)

      if ((rc /= nf90_noerr) .and. (rc /= nf90_enotatt)) &

        call mom_error(fatal, trim(hdr) // trim(nf90_strerror(rc)) //&

                " Difficulties reading global att "//trim(attname)//" from "//trim(filename))

    endif

    if (present(found)) found = (rc == nf90_noerr)

    if (.not.present(ncid_in)) call close_file_to_read(ncid, filename)

  endif

  if (do_broadcast) then

    if (present(found)) then

      is_found = 0 ; if (is_root_pe() .and. found) is_found = 1

      call broadcast(is_found, blocking=.false.)

    endif

    call broadcast(att_val, blocking=.true.)

    if (present(found)) found = (is_found /= 0)

  endif

end subroutine read_attribute_int32

!> Read a 64-bit integer global or variable attribute

subroutine read_attribute_int64(filename, attname, att_val, varname, found, all_read, ncid_in)

  character(len=*),           intent(in)  :: filename !< Name of the file to read

  character(len=*),           intent(in)  :: attname  !< Name of the attribute to read

  integer(kind=int64),        intent(out) :: att_val  !< The value of the attribute

  character(len=*), optional, intent(in)  :: varname  !< The name of the variable whose attribute will

                                                      !! be read. If missing, read a global attribute.

  logical,          optional, intent(out) :: found    !< Returns true if the attribute is found

  logical,          optional, intent(in)  :: all_read !< If present and true, all PEs that call this

                                                      !! routine actually do the read, otherwise only

                                                      !! root PE reads and then broadcasts the results.

  integer,        optional, intent(in)    :: ncid_in  !< The netCDF ID of an open file.  If absent, the

                                                      !! file is opened and closed within this routine.

  logical :: do_read, do_broadcast

  integer :: rc, ncid, varid, is_found

  character(len=256) :: hdr

  hdr = "read_attribute_int64"

  att_val = 0

  do_read = is_root_pe() ; if (present(all_read)) do_read = all_read .or. do_read

  do_broadcast = .true. ; if (present(all_read)) do_broadcast = .not.all_read

  if (do_read) then

    if (present(ncid_in)) then

      ncid = ncid_in

    else

      call open_file_to_read(filename, ncid, success=found)

      if (present(found)) then ; if (.not.found) do_read = .false. ; endif

    endif

  endif

  if (do_read) then

    rc = nf90_enotatt

    if (present(varname)) then  ! Read a variable attribute

      call get_varid(varname, ncid, filename, varid, match_case=.false., found=found)

      if (varid >= 0) then ! The named variable does exist, and found would be true.

        rc = nf90_get_att(ncid, varid, attname, att_val)

        if ((rc /= nf90_noerr) .and. (rc /= nf90_enotatt)) &

          call mom_error(fatal, trim(hdr) // trim(nf90_strerror(rc)) //" Difficulties reading att "//&

                trim(attname)//" for "//trim(varname)//" from "//trim(filename))

      endif

    else  ! Read a global attribute

      rc = nf90_get_att(ncid, nf90_global, attname, att_val)

      if ((rc /= nf90_noerr) .and. (rc /= nf90_enotatt)) &

        call mom_error(fatal, trim(hdr) // trim(nf90_strerror(rc)) //&

                " Difficulties reading global att "//trim(attname)//" from "//trim(filename))

    endif

    if (present(found)) found = (rc == nf90_noerr)

    rc = nf90_close(ncid)

  endif

  if (do_broadcast) then

    if (present(found)) then

      is_found = 0 ; if (is_root_pe() .and. found) is_found = 1

      call broadcast(is_found, blocking=.false.)

    endif

    call broadcast(att_val, blocking=.true.)

    if (present(found)) found = (is_found /= 0)

  endif

end subroutine read_attribute_int64

!> Read a real global or variable attribute

subroutine read_attribute_real(filename, attname, att_val, varname, found, all_read, ncid_in)

  character(len=*),           intent(in)  :: filename !< Name of the file to read

  character(len=*),           intent(in)  :: attname  !< Name of the attribute to read

  real,                       intent(out) :: att_val  !< The value of the attribute [arbitrary]

  character(len=*), optional, intent(in)  :: varname  !< The name of the variable whose attribute will

                                                      !! be read. If missing, read a global attribute.

  logical,          optional, intent(out) :: found    !< Returns true if the attribute is found

  logical,          optional, intent(in)  :: all_read !< If present and true, all PEs that call this

                                                      !! routine actually do the read, otherwise only

                                                      !! root PE reads and then broadcasts the results.

  integer,          optional, intent(in)  :: ncid_in  !< The netCDF ID of an open file.  If absent, the

                                                      !! file is opened and closed within this routine.

  logical :: do_read, do_broadcast

  integer :: rc, ncid, varid, is_found

  character(len=256) :: hdr

  hdr = "read_attribute_real"

  att_val = 0.0

  do_read = is_root_pe() ; if (present(all_read)) do_read = all_read .or. do_read

  do_broadcast = .true. ; if (present(all_read)) do_broadcast = .not.all_read

  if (do_read) then

    if (present(ncid_in)) then

      ncid = ncid_in

    else

      call open_file_to_read(filename, ncid, success=found)

      if (present(found)) then ; if (.not.found) do_read = .false. ; endif

    endif

  endif

  if (do_read) then

    rc = nf90_enotatt

    if (present(varname)) then  ! Read a variable attribute

      call get_varid(varname, ncid, filename, varid, match_case=.false., found=found)

      if (varid >= 0) then ! The named variable does exist, and found would be true.

        rc = nf90_get_att(ncid, varid, attname, att_val)

        if ((rc /= nf90_noerr) .and. (rc /= nf90_enotatt)) &

          call mom_error(fatal, trim(hdr) // trim(nf90_strerror(rc)) //" Difficulties reading att "//&

                trim(attname)//" for "//trim(varname)//" from "//trim(filename))

      endif

    else  ! Read a global attribute

      rc = nf90_get_att(ncid, nf90_global, attname, att_val)

      if ((rc /= nf90_noerr) .and. (rc /= nf90_enotatt)) &

        call mom_error(fatal, trim(hdr) // trim(nf90_strerror(rc)) //&

                " Difficulties reading global att "//trim(attname)//" from "//trim(filename))

    endif

    if (present(found)) found = (rc == nf90_noerr)

    if (.not.present(ncid_in)) call close_file_to_read(ncid, filename)

  endif

  if (do_broadcast) then

    if (present(found)) then

      is_found = 0 ; if (is_root_pe() .and. found) is_found = 1

      call broadcast(is_found, blocking=.false.)

    endif

    call broadcast(att_val, blocking=.true.)

    if (present(found)) found = (is_found /= 0)

  endif

end subroutine read_attribute_real

!> Open a netcdf file for reading, with error handling

subroutine open_file_to_read(filename, ncid, success)

  character(len=*),  intent(in)  :: filename   !< path and name of the file to open for reading

  integer,           intent(out) :: ncid       !< The netcdf handle for the file

  logical, optional, intent(out) :: success    !< Returns true if the file was opened, or if this

                                               !! argument is not present, failure is fatal error.

  ! Local variables

  integer rc

  rc = nf90_open(trim(filename), nf90_nowrite, ncid)

  if (present(success)) then

    success = (rc == nf90_noerr)

  elseif (rc /= nf90_noerr) then

    call mom_error(fatal, "Difficulties opening "//trim(filename)//" - "//trim(nf90_strerror(rc)) )

  endif

end subroutine open_file_to_read

!> Close a netcdf file that had been opened for reading, with error handling

subroutine close_file_to_read(ncid, filename)

  integer,                    intent(inout) :: ncid       !< The netcdf handle for the file to close

  character(len=*), optional, intent(in)    :: filename   !< path and name of the file to close

  integer :: rc

  if (ncid >= 0) then

    rc = nf90_close(ncid)

    if (present(filename) .and. (rc /= nf90_noerr)) then

      call mom_error(warning, "Difficulties closing "//trim(filename)//": "//trim(nf90_strerror(rc)))

    elseif (rc /= nf90_noerr) then

      call mom_error(warning, "Difficulties closing file: "//trim(nf90_strerror(rc)))

    endif

  endif

  ncid = -1

end subroutine close_file_to_read

!> get_varid finds the netcdf handle for the potentially case-insensitive variable name in a file

subroutine get_varid(varname, ncid, filename, varid, match_case, found)

  character(len=*),  intent(in)  :: varname    !< The name of the variable that is being sought

  integer,           intent(in)  :: ncid       !< The open netcdf handle for the file

  character(len=*),  intent(in)  :: filename   !< name of the file to read, used here in messages

  integer,           intent(out) :: varid      !< The netcdf handle for the variable

  logical, optional, intent(in)  :: match_case !< If false, allow for variables name matches to be

                                               !! case insensitive, but take a perfect match if

                                               !! found.  The default is true.

  logical, optional, intent(out) :: found      !< Returns true if the attribute is found

  logical :: var_found, insensitive

  character(len=256) :: name

  integer, allocatable :: varids(:)

  integer :: nvars, status, n

  varid = -1

  var_found = .false.

  insensitive = .false. ; if (present(match_case)) insensitive = .not.match_case

  if (insensitive) then

    ! This code ounddoes a case-insensitive search for a variable in the file.

    status = nf90_inquire(ncid, nvariables=nvars)

    if (present(found) .and. ((status /= nf90_noerr) .or. (nvars < 1))) then

      found = .false. ; return

    elseif (status /= nf90_noerr) then

      call mom_error(fatal, "get_varid:  Difficulties getting the number of variables in file "//&

          trim(filename)//" - "//trim(nf90_strerror(status)))

    elseif (nvars < 1) then

      call mom_error(fatal, "get_varid: There appear not to be any variables in "//trim(filename))

    endif

    allocate(varids(nvars))

    status = nf90_inq_varids(ncid, nvars, varids)

    if (status /= nf90_noerr) then

      call mom_error(warning, "get_varid: Difficulties getting the variable IDs in file "//&

          trim(filename)//" - "//trim(nf90_strerror(status)))

      nvars = -1  ! Full error handling will occur after the do-loop.

    endif

    do n = 1,nvars

      status = nf90_inquire_variable(ncid, varids(n), name=name)

      if (status /= nf90_noerr) then

        call mom_error(warning, "get_varid:  Difficulties getting a variable name in file "//&

            trim(filename)//" - "//trim(nf90_strerror(status)))

      endif

      if (trim(lowercase(name)) == trim(lowercase(varname))) then

        if (var_found) then

          call mom_error(warning, "get_varid: Two variables match the case-insensitive name "//&

                  trim(varname)//" in file "//trim(filename))

          ! Replace the first variable if the second one is a case-sensitive match

          if (trim(name) == trim(varname)) varid = varids(n)

        else

          varid = varids(n) ; var_found = .true.

        endif

      endif

    enddo

    if (present(found)) found = var_found

    if ((.not.var_found) .and. .not.present(found)) call mom_error(fatal, &

        "get_varid: variable "//trim(varname)//" was not found in file "//trim(filename))

    deallocate(varids)

  else

    status = nf90_inq_varid(ncid, trim(varname), varid)

    if (present(found)) found = (status == nf90_noerr)

    if ((status /= nf90_noerr) .and. .not.present(found)) then

      call mom_error(fatal, "get_varid: Difficulties getting a variable id for "//&

          trim(varname)//" in file "//trim(filename)//" - "//trim(nf90_strerror(status)))

    endif

  endif

end subroutine get_varid

!> Verify that a file contains a named variable with the expected units.

subroutine verify_variable_units(filename, varname, expected_units, msg, ierr, alt_units)

  character(len=*),           intent(in)    :: filename  !< File name

  character(len=*),           intent(in)    :: varname   !< Variable name

  character(len=*),           intent(in)    :: expected_units !< Expected units of variable

  character(len=*),           intent(inout) :: msg       !< Message to use for errors

  logical,                    intent(out)   :: ierr      !< True if an error occurs

  character(len=*), optional, intent(in)    :: alt_units !< Alterate acceptable units of variable

  ! Local variables

  character (len=200) :: units

  logical :: units_correct, success

  integer :: i, ncid, status, vid

  if (.not.is_root_pe()) then ! Only the root PE should do the verification.

    ierr = .false. ; msg = '' ; return

  endif

  ierr = .true.

  call open_file_to_read(filename, ncid, success)

  if (.not.success) then

    msg = 'File not found: '//trim(filename)

    return

  endif

  status = nf90_inq_varid(ncid, trim(varname), vid)

  if (status /= nf90_noerr) then

    msg = 'Var not found: '//trim(varname)

  else

    status = nf90_get_att(ncid, vid, "units", units)

    if (status /= nf90_noerr) then

      msg = 'Attribute not found: units'

    else

      ! NF90_GET_ATT can return attributes with null characters, which TRIM will not truncate.

      ! This loop replaces any null characters with a space so that the subsequent check

      ! between the read units and the expected units will pass

      do i=1,len_trim(units)

        if (units(i:i) == char(0)) units(i:i) = " "

      enddo

      units_correct = (trim(units) == trim(expected_units))

      if (present(alt_units)) then

        units_correct = units_correct .or. (trim(units) == trim(alt_units))

      endif

      if (units_correct) then

        ierr = .false.

        msg = ''

      else

        msg = 'Units incorrect: '//trim(units)//' /= '//trim(expected_units)

      endif

    endif

  endif

  status = nf90_close(ncid)

end subroutine verify_variable_units

!> Returns a vardesc type whose elements have been filled with the provided

!! fields.  The argument name is required, while the others are optional and

!! have default values that are empty strings or are appropriate for a 3-d

!! tracer field at the tracer cell centers.

function var_desc(name, units, longname, hor_grid, z_grid, t_grid, cmor_field_name, &

                  cmor_units, cmor_longname, conversion, caller, position, dim_names, &

                  extra_axes, fixed) result(vd)

  character(len=*),           intent(in) :: name            !< variable name

  character(len=*), optional, intent(in) :: units           !< variable units

  character(len=*), optional, intent(in) :: longname        !< variable long name

  character(len=*), optional, intent(in) :: hor_grid        !< A character string indicating the horizontal

                                                            !! position of this variable

  character(len=*), optional, intent(in) :: z_grid          !< variable vertical staggering

  character(len=*), optional, intent(in) :: t_grid          !< time description: s, p, or 1

  character(len=*), optional, intent(in) :: cmor_field_name !< CMOR name

  character(len=*), optional, intent(in) :: cmor_units      !< CMOR physical dimensions of variable

  character(len=*), optional, intent(in) :: cmor_longname   !< CMOR long name

  real            , optional, intent(in) :: conversion      !< for unit conversions, such as needed to

                                                            !! convert from intensive to extensive

                                                            !! [various] or [a A-1 ~> 1]

  character(len=*), optional, intent(in) :: caller          !< The calling routine for error messages

  integer,          optional, intent(in) :: position        !< A coded integer indicating the horizontal position

                                                            !! of this variable if it has such dimensions.

                                                            !! Valid values include CORNER, CENTER, EAST_FACE

                                                            !! NORTH_FACE, and 0 for no horizontal dimensions.

  character(len=*), dimension(:), &

                    optional, intent(in) :: dim_names       !< The names of the dimensions of this variable

  type(axis_info),  dimension(:), &

                    optional, intent(in) :: extra_axes      !< dimensions other than space-time

  logical,          optional, intent(in) :: fixed           !< If true, this does not evolve with time

  type(vardesc)                          :: vd              !< vardesc type that is created

  character(len=120) :: cllr

  cllr = "var_desc"

  if (present(caller)) cllr = trim(caller)

  call safe_string_copy(name, vd%name, "vd%name", cllr)

  vd%longname = "" ; vd%units = ""

  vd%hor_grid = 'h' ; vd%position = center ; vd%z_grid = 'L' ; vd%t_grid = 's'

  if (present(dim_names)) vd%z_grid = '1' ! In this case the names are used to set the non-horizontal axes

  if (present(fixed)) then ; if (fixed) vd%t_grid = '1' ; endif

  vd%cmor_field_name  =  ""

  vd%cmor_units       =  ""

  vd%cmor_longname    =  ""

  vd%conversion       =  1.0

  vd%dim_names(:)     =  ""

  call modify_vardesc(vd, units=units, longname=longname, hor_grid=hor_grid, &

                      z_grid=z_grid, t_grid=t_grid, position=position, dim_names=dim_names, &

                      cmor_field_name=cmor_field_name, cmor_units=cmor_units, &

                      cmor_longname=cmor_longname, conversion=conversion, caller=cllr, &

                      extra_axes=extra_axes)

end function var_desc

!> This routine modifies the named elements of a vardesc type.

!! All arguments are optional, except the vardesc type to be modified.

subroutine modify_vardesc(vd, name, units, longname, hor_grid, z_grid, t_grid, &

                 cmor_field_name, cmor_units, cmor_longname, conversion, caller, position, dim_names, &

                 extra_axes)

  type(vardesc),              intent(inout) :: vd              !< vardesc type that is modified

  character(len=*), optional, intent(in)    :: name            !< name of variable

  character(len=*), optional, intent(in)    :: units           !< units of variable

  character(len=*), optional, intent(in)    :: longname        !< long name of variable

  character(len=*), optional, intent(in)    :: hor_grid        !< horizontal staggering of variable

  character(len=*), optional, intent(in)    :: z_grid          !< vertical staggering of variable

  character(len=*), optional, intent(in)    :: t_grid          !< time description: s, p, or 1

  character(len=*), optional, intent(in)    :: cmor_field_name !< CMOR name

  character(len=*), optional, intent(in)    :: cmor_units      !< CMOR physical dimensions of variable

  character(len=*), optional, intent(in)    :: cmor_longname   !< CMOR long name

  real            , optional, intent(in)    :: conversion      !< A multiplicative factor for unit conversions,

                                                               !! such as needed to convert from intensive to

                                                               !! extensive or dimensional consistency testing

                                                               !! [various] or [a A-1 ~> 1]

  character(len=*), optional, intent(in)    :: caller          !< The calling routine for error messages

  integer,          optional, intent(in)    :: position        !< A coded integer indicating the horizontal position

                                                               !! of this variable if it has such dimensions.

                                                               !! Valid values include CORNER, CENTER, EAST_FACE

                                                               !! NORTH_FACE, and 0 for no horizontal dimensions.

  character(len=*), dimension(:), &

                    optional, intent(in)    :: dim_names       !< The names of the dimensions of this variable

  type(axis_info),  dimension(:), &

                    optional, intent(in)    :: extra_axes      !< dimensions other than space-time

  character(len=120) :: cllr

  integer :: n

  cllr = "mod_vardesc" ; if (present(caller)) cllr = trim(caller)

  if (present(name))      call safe_string_copy(name, vd%name, "vd%name", cllr)

  if (present(longname))  call safe_string_copy(longname, vd%longname, &

                               "vd%longname of "//trim(vd%name), cllr)

  if (present(units))     call safe_string_copy(units, vd%units,       &

                               "vd%units of "//trim(vd%name), cllr)

  if (present(position)) then

    vd%position = position

    select case (position)

      case (center)     ; vd%hor_grid = 'T'

      case (corner)     ; vd%hor_grid = 'Bu'

      case (east_face)  ; vd%hor_grid = 'Cu'

      case (north_face) ; vd%hor_grid = 'Cv'

      case (0)          ; vd%hor_grid = '1'

      case default

        call mom_error(fatal, "modify_vardesc: "//trim(vd%name)//" has unrecognized position argument")

    end select

  endif

  if (present(hor_grid)) then

    call safe_string_copy(hor_grid, vd%hor_grid, "vd%hor_grid of "//trim(vd%name), cllr)

    vd%position = position_from_horgrid(vd%hor_grid)

    if (present(caller) .and. (vd%position == -1)) then

      call mom_error(fatal, "modify_vardesc called by "//trim(caller)//": "//trim(vd%name)//&

                   " has an unrecognized hor_grid argument "//trim(vd%hor_grid))

    elseif (vd%position == -1) then

      call mom_error(fatal, "modify_vardesc called with bad hor_grid argument "//trim(vd%hor_grid))

    endif

  endif

  if (present(z_grid))    call safe_string_copy(z_grid, vd%z_grid,     &

                               "vd%z_grid of "//trim(vd%name), cllr)

  if (present(t_grid))    call safe_string_copy(t_grid, vd%t_grid,     &

                               "vd%t_grid of "//trim(vd%name), cllr)

  if (present(cmor_field_name)) call safe_string_copy(cmor_field_name, vd%cmor_field_name, &

                                     "vd%cmor_field_name of "//trim(vd%name), cllr)

  if (present(cmor_units))      call safe_string_copy(cmor_units, vd%cmor_units, &

                                     "vd%cmor_units of "//trim(vd%name), cllr)

  if (present(cmor_longname))   call safe_string_copy(cmor_longname, vd%cmor_longname, &

                                     "vd%cmor_longname of "//trim(vd%name), cllr)

  if (present(conversion)) vd%conversion = conversion

  if (present(dim_names)) then

    do n=1,min(5,size(dim_names)) ; if (len_trim(dim_names(n)) > 0) then

      call safe_string_copy(dim_names(n), vd%dim_names(n), "vd%dim_names of "//trim(vd%name), cllr)

    endif ; enddo

  endif

  if (present(extra_axes)) then

    do n=1,size(extra_axes) ; if (len_trim(extra_axes(n)%name) > 0) then

      vd%extra_axes(n) = extra_axes(n)

    endif ; enddo

  endif

end subroutine modify_vardesc

integer function position_from_horgrid(hor_grid)

  character(len=*), intent(in)    :: hor_grid        !< horizontal staggering of variable

  select case (trim(hor_grid))

    case ('h')   ; position_from_horgrid = center

    case ('q')   ; position_from_horgrid = corner

    case ('u')   ; position_from_horgrid = east_face

    case ('v')   ; position_from_horgrid = north_face

    case ('T')   ; position_from_horgrid = center

    case ('Bu')  ; position_from_horgrid = corner

    case ('Cu')  ; position_from_horgrid = east_face

    case ('Cv')  ; position_from_horgrid = north_face

    case ('1')   ; position_from_horgrid = 0

    case default ; position_from_horgrid = -1 ! This is a bad-value flag.

  end select

end function position_from_horgrid

!> Store information that can be used to create an axis in a subsequent call to create_file.

subroutine set_axis_info(axis, name, units, longname, ax_size, ax_data, cartesian, sense)

  type(axis_info),              intent(inout) :: axis  !< A type with information about a named axis

  character(len=*),             intent(in)    :: name  !< The name of this axis for use in files

  character(len=*),   optional, intent(in)    :: units !< The units of the axis labels

  character(len=*),   optional, intent(in)    :: longname  !< Long name of the axis variable

  integer,            optional, intent(in)    :: ax_size !< The number of elements in this axis

  real, dimension(:), optional, intent(in)    :: ax_data !< The values of the data on the axis [arbitrary]

  character(len=*),   optional, intent(in)    :: cartesian !< A variable indicating which direction this axis

                                                       !! axis corresponds with. Valid values

                                                       !! include 'X', 'Y', 'Z', 'T', and 'N' (the default) for none.

  integer,            optional, intent(in)    :: sense !< This is 1 for axes whose values increase upward, or -1

                                                       !! if they increase downward.  The default, 0, is ignored.

  call safe_string_copy(name, axis%name, "axis%name of "//trim(name), "set_axis_info")

  ! Set the default values.

  axis%longname = trim(axis%name) ; axis%units = "" ; axis%cartesian = "N" ; axis%sense = 0

  if (present(longname))  call safe_string_copy(longname, axis%longname, &

                                                "axis%longname of "//trim(name), "set_axis_info")

  if (present(units))     call safe_string_copy(units, axis%units, &

                                                "axis%units of "//trim(name), "set_axis_info")

  if (present(cartesian)) call safe_string_copy(cartesian, axis%cartesian, &

                                                "axis%cartesian of "//trim(name), "set_axis_info")

  if (present(sense)) axis%sense = sense

  if (.not.(present(ax_size) .or. present(ax_data)) ) then

    call mom_error(fatal, "set_axis_info called for "//trim(name)//&

      "without either an ax_size or an ax_data argument.")

  elseif (present(ax_size) .and. present(ax_data)) then

    if (size(ax_data) /= ax_size) call mom_error(fatal, "set_axis_info called for "//trim(name)//&

      "with an inconsistent value of ax_size and size of ax_data.")

  endif

  if (present(ax_size)) then

    axis%ax_size = ax_size

  else

    axis%ax_size = size(ax_data)

  endif

  if (present(ax_data)) then

    allocate(axis%ax_data(axis%ax_size)) ; axis%ax_data(:) = ax_data(:)

  endif

end subroutine set_axis_info

!> Delete the information in an array of axis_info types and deallocate memory in them.

subroutine delete_axis_info(axes)

  type(axis_info), dimension(:), intent(inout) :: axes  !< An array with information about named axes

  integer :: n

  do n=1,size(axes)

    axes(n)%name = "" ; axes(n)%longname = "" ; axes(n)%units = "" ; axes(n)%cartesian = "N"

    axes(n)%sense = 0 ; axes(n)%ax_size = 0

    if (allocated(axes(n)%ax_data)) deallocate(axes(n)%ax_data)

  enddo

end subroutine delete_axis_info

!> Retrieve the information from an axis_info type.

subroutine get_axis_info(axis,name,longname,units,cartesian,ax_size,ax_data)

  type(axis_info), intent(in) :: axis                               !< An axis type

  character(len=*), intent(out), optional    :: name                !< The axis name.

  character(len=*), intent(out), optional    :: longname            !< The axis longname.

  character(len=*), intent(out), optional    :: units               !< The axis units.

  character(len=*), intent(out), optional    :: cartesian           !< The cartesian attribute

                                                                    !! of the axis [X,Y,Z,T].

  integer,          intent(out), optional   :: ax_size              !< The size of the axis.

  real, optional, allocatable, dimension(:), intent(out) :: ax_data !< The axis label data [arbitrary]

  if (present(ax_data)) then

    if (allocated(ax_data)) deallocate(ax_data)

    allocate(ax_data(axis%ax_size))

    ax_data(:) = axis%ax_data

  endif

  if (present(name)) name = axis%name

  if (present(longname)) longname = axis%longname

  if (present(units)) units = axis%units

  if (present(cartesian)) cartesian = axis%cartesian

  if (present(ax_size)) ax_size = axis%ax_size

end subroutine get_axis_info

!> Store information that can be used to create an attribute in a subsequent call to create_file.

subroutine set_attribute_info(attribute, name, str_value)

  type(attribute_info), intent(inout) :: attribute !< A type with information about a named attribute

  character(len=*),     intent(in)    :: name      !< The name of this attribute for use in files

  character(len=*),     intent(in)    :: str_value !< The value of this attribute

  attribute%name = trim(name)

  attribute%att_val = trim(str_value)

end subroutine set_attribute_info

!> Delete the information in an array of attribute_info types and deallocate memory in them.

subroutine delete_attribute_info(atts)

  type(attribute_info), dimension(:), intent(inout) :: atts  !< An array of global attributes

  integer :: n

  do n=1,size(atts)

    if (allocated(atts(n)%name)) deallocate(atts(n)%name)

    if (allocated(atts(n)%att_val)) deallocate(atts(n)%att_val)

  enddo

end subroutine delete_attribute_info

!> This function returns the CMOR standard name given a CMOR longname, based on

!! the standard pattern of character conversions.

function cmor_long_std(longname) result(std_name)

  character(len=*), intent(in) :: longname  !< The CMOR longname being converted

  character(len=len(longname)) :: std_name  !< The CMOR standard name generated from longname

  integer :: k

  std_name = lowercase(longname)

  do k=1, len_trim(std_name)

    if (std_name(k:k) == ' ') std_name(k:k) = '_'

  enddo

end function cmor_long_std

!> This routine queries vardesc

subroutine query_vardesc(vd, name, units, longname, hor_grid, z_grid, t_grid, &

                         cmor_field_name, cmor_units, cmor_longname, conversion, caller, &

                         extra_axes, position, dim_names)

  type(vardesc),              intent(in)  :: vd                 !< vardesc type that is queried

  character(len=*), optional, intent(out) :: name               !< name of variable

  character(len=*), optional, intent(out) :: units              !< units of variable

  character(len=*), optional, intent(out) :: longname           !< long name of variable

  character(len=*), optional, intent(out) :: hor_grid           !< horizontal staggering of variable

  character(len=*), optional, intent(out) :: z_grid             !< verticle staggering of variable

  character(len=*), optional, intent(out) :: t_grid             !< time description: s, p, or 1

  character(len=*), optional, intent(out) :: cmor_field_name    !< CMOR name

  character(len=*), optional, intent(out) :: cmor_units         !< CMOR physical dimensions of variable

  character(len=*), optional, intent(out) :: cmor_longname      !< CMOR long name

  real            , optional, intent(out) :: conversion         !< for unit conversions, such as needed to

                                                                !! convert from intensive to extensive

                                                                !! [various] or [a A-1 ~> 1]

  character(len=*), optional, intent(in)  :: caller             !< calling routine?

  type(axis_info),  dimension(5), &

                    optional, intent(out) :: extra_axes      !< dimensions other than space-time

  integer,          optional, intent(out) :: position        !< A coded integer indicating the horizontal position

                                                            !! of this variable if it has such dimensions.

                                                            !! Valid values include CORNER, CENTER, EAST_FACE

                                                            !! NORTH_FACE, and 0 for no horizontal dimensions.

  character(len=*), dimension(:), &

                    optional, intent(out) :: dim_names       !< The names of the dimensions of this variable

  integer :: n

  integer, parameter :: nmax_extraaxes = 5

  character(len=120) :: cllr, varname

  cllr = "mod_vardesc"

  if (present(caller)) cllr = trim(caller)

  if (present(name))      call safe_string_copy(vd%name, name,         &

                               "vd%name of "//trim(vd%name), cllr)

  if (present(longname))  call safe_string_copy(vd%longname, longname, &

                               "vd%longname of "//trim(vd%name), cllr)

  if (present(units))     call safe_string_copy(vd%units, units,       &

                               "vd%units of "//trim(vd%name), cllr)

  if (present(hor_grid))  call safe_string_copy(vd%hor_grid, hor_grid, &

                               "vd%hor_grid of "//trim(vd%name), cllr)

  if (present(z_grid))    call safe_string_copy(vd%z_grid, z_grid,     &

                               "vd%z_grid of "//trim(vd%name), cllr)

  if (present(t_grid))    call safe_string_copy(vd%t_grid, t_grid,     &

                               "vd%t_grid of "//trim(vd%name), cllr)

  if (present(cmor_field_name)) call safe_string_copy(vd%cmor_field_name, cmor_field_name, &

                                     "vd%cmor_field_name of "//trim(vd%name), cllr)

  if (present(cmor_units))      call safe_string_copy(vd%cmor_units, cmor_units,          &

                                     "vd%cmor_units of "//trim(vd%name), cllr)

  if (present(cmor_longname))   call safe_string_copy(vd%cmor_longname, cmor_longname, &

                                     "vd%cmor_longname of "//trim(vd%name), cllr)

  if (present(conversion)) conversion = vd%conversion

  if (present(position)) then

    position = vd%position

    if (position == -1) position = position_from_horgrid(vd%hor_grid)

  endif

  if (present(dim_names)) then

    do n=1,min(5,size(dim_names))

      call safe_string_copy(vd%dim_names(n), dim_names(n), "vd%dim_names of "//trim(vd%name), cllr)

    enddo

  endif

  if (present(extra_axes)) then

    ! save_restart expects 5 extra axes (can be empty)

    do n=1, nmax_extraaxes

      if (vd%extra_axes(n)%ax_size>=1) then

        extra_axes(n) = vd%extra_axes(n)

      else

        ! return an empty axis

        write(varname,"('dummy',i1.1)") n

        call set_axis_info(extra_axes(n), name=trim(varname), ax_size=1)

      endif

    enddo

  endif

end subroutine query_vardesc

!> Read a scalar from file using infrastructure I/O.

subroutine mom_read_data_0d(filename, fieldname, data, timelevel, scale, MOM_Domain, &

                            global_file, file_may_be_4d)

  character(len=*), intent(in)  :: filename     !< Input filename

  character(len=*), intent(in)  :: fieldname    !< Field variable name

  real, intent(inout)           :: data         !< Field value in arbitrary units [A ~> a]

  integer, optional, intent(in) :: timelevel    !< Time level to read in file

  real, optional, intent(in)    :: scale        !< A scaling factor that the variable is multiplied by

                                                !! before it is returned to convert from the units in the file

                                                !! to the internal units for this variable [A a-1 ~> 1]

  type(mom_domain_type), optional, intent(in) :: MOM_Domain !< Model domain decomposition

  logical, optional, intent(in) :: global_file    !< If true, read from a single file

  logical, optional, intent(in) :: file_may_be_4d !< If true, fields may be stored

                                                  !! as 4d arrays in the file.

  call read_field(filename, fieldname, data, &

                  timelevel=timelevel, scale=scale, mom_domain=mom_domain, &

                  global_file=global_file, file_may_be_4d=file_may_be_4d)

end subroutine mom_read_data_0d

!> Read a scalar integer from file using infrastructure I/O.

subroutine mom_read_data_0d_int(filename, fieldname, data, timelevel)

  character(len=*), intent(in) :: filename    !< Input filename

  character(len=*), intent(in) :: fieldname   !< Field variable name

  integer, intent(inout) :: data              !< Field value

  integer, optional, intent(in) :: timelevel  !< Time level to read in file

  call read_field(filename, fieldname, data, timelevel=timelevel)

end subroutine mom_read_data_0d_int

!> Read a 1d array from file using infrastructure I/O.

subroutine mom_read_data_1d(filename, fieldname, data, timelevel, scale, MOM_Domain, &

                            global_file, file_may_be_4d)

  character(len=*), intent(in)  :: filename   !< Input filename

  character(len=*), intent(in)  :: fieldname  !< Field variable name

  real, dimension(:), intent(inout) :: data   !< Field value in arbitrary units [A ~> a]

  integer, optional, intent(in) :: timelevel  !< Time level to read in file

  real, optional, intent(in)    :: scale      !< A scaling factor that the variable is multiplied by

                                              !! before it is returned to convert from the units in the file

                                              !! to the internal units for this variable [A a-1 ~> 1]

  type(mom_domain_type), optional, intent(in) :: MOM_Domain !< Model domain decomposition

  logical, optional, intent(in) :: global_file    !< If true, read from a single file

  logical, optional, intent(in) :: file_may_be_4d !< If true, fields may be stored

                                                  !! as 4d arrays in the file.

  call read_field(filename, fieldname, data, &

                  timelevel=timelevel, scale=scale, mom_domain=mom_domain, &

                  global_file=global_file, file_may_be_4d=file_may_be_4d)

end subroutine mom_read_data_1d

!> Read a 1d integer array from file using infrastructure I/O.

subroutine mom_read_data_1d_int(filename, fieldname, data, timelevel)

  character(len=*), intent(in) :: filename    !< Input filename

  character(len=*), intent(in) :: fieldname   !< Field variable name

  integer, dimension(:), intent(inout) :: data  !< Field value

  integer, optional, intent(in) :: timelevel  !< Time level to read in file

  call read_field(filename, fieldname, data, timelevel=timelevel)

end subroutine mom_read_data_1d_int

!> Read a 2d array from file using infrastructure I/O.

subroutine mom_read_data_2d(filename, fieldname, data, MOM_Domain, timelevel, position, &

                            scale, global_file, file_may_be_4d, turns)

  character(len=*), intent(in)  :: filename  !< Input filename

  character(len=*), intent(in)  :: fieldname !< Field variable name

  real, dimension(:,:), intent(inout) :: data   !< Field value in arbitrary units [A ~> a]

  type(mom_domain_type), target, &

                     intent(in) :: MOM_Domain !< Model domain decomposition

  integer, optional, intent(in) :: timelevel !< Time level to read in file

  integer, optional, intent(in) :: position  !< Grid positioning flag

  real, optional, intent(in)    :: scale     !< A scaling factor that the variable is multiplied by

                                             !! before it is returned to convert from the units in the file

                                             !! to the internal units for this variable [A a-1 ~> 1]

  logical, optional, intent(in) :: global_file    !< If true, read from a single file

  logical, optional, intent(in) :: file_may_be_4d !< If true, fields may be stored

                                                  !! as 4d arrays in the file.

  integer, optional, intent(in) :: turns        !< Number of quarter-turns to rotate the data.  If absent

                                                !! the number of turns is taken from MOM_Domain.

  ! Local variables

  integer :: qturns   ! Number of quarter-turns from input to model grid

  real, allocatable :: data_in(:,:)  ! Field array on the input grid in arbitrary units [A ~> a]

  type(mom_domain_type), pointer :: domain_ptr => null()  ! Pointer to the unrotated domain for reading

  qturns = mom_domain%turns ; if (present(turns)) qturns = modulo(turns, 4)

  domain_ptr => mom_domain

  if (associated(mom_domain%domain_in) .and. (qturns /= 0)) domain_ptr => mom_domain%domain_in

  if (qturns == 0) then

    call read_field(filename, fieldname, data, mom_domain, &

                    timelevel=timelevel, position=position, scale=scale, &

                    global_file=global_file, file_may_be_4d=file_may_be_4d)

  else

    call allocate_rotated_array(data, [1,1], -qturns, data_in)

    call rotate_array(data, -qturns, data_in)

    call read_field(filename, fieldname, data_in, domain_ptr, &

                    timelevel=timelevel, position=position, scale=scale, &

                    global_file=global_file, file_may_be_4d=file_may_be_4d)

    call rotate_array(data_in, qturns, data)

    deallocate(data_in)

  endif

end subroutine mom_read_data_2d

!> Read a 2d array (which might have halos) from a file using native netCDF I/O.

subroutine read_netcdf_data_2d(filename, fieldname, values, MOM_Domain, &

                            timelevel, position, rescale, turns)

  character(len=*), intent(in) :: filename

    !< Input filename

  character(len=*), intent(in)  :: fieldname

    !< Field variable name

  real, intent(inout) :: values(:,:)

    !< Field values read from the file.  It would be intent(out) but for the

    !! need to preserve any initialized values in the halo regions.

  type(mom_domain_type), intent(in) :: MOM_Domain

    !< Model domain decomposition

  integer, optional, intent(in) :: timelevel

    !< Time level to read in file

  integer, optional, intent(in) :: position

    !< Grid positioning flag

  real, optional, intent(in) :: rescale

    !< Rescale factor, omitting this is the same as setting it to 1.

  integer, optional, intent(in) :: turns

    !< Number of quarter-turns to rotate the data.  If absent the number of turns is taken

    !! from MOM_Domain.

  integer :: qturns

    ! Number of quarter-turns from input to model grid

  real, allocatable :: values_in(:,:)

    ! Field array on the unrotated input grid

  type(mom_netcdf_file) :: handle

    ! netCDF file handle

  ! General-purpose IO will require the following arguments, but they are not

  ! yet implemented, so we raise an error if they are present.

  ! Fields are currently assumed on cell centers, and position is unsupported

  if (present(position)) &

    call mom_error(fatal, 'read_netCDF_data: position is not yet supported.')

  ! Timelevels are not yet supported

  if (present(timelevel)) &

    call mom_error(fatal, 'read_netCDF_data: timelevel is not yet supported.')

  call handle%open(filename, action=readonly_file, mom_domain=mom_domain)

  call handle%update()

  qturns = mom_domain%turns ; if (present(turns)) qturns = modulo(turns, 4)

  if (qturns == 0) then

    call handle%read(fieldname, values, rescale=rescale)

  else

    call allocate_rotated_array(values, [1,1], -qturns, values_in)

    call rotate_array(values, -qturns, values_in)

    call handle%read(fieldname, values_in, rescale=rescale)

    call rotate_array(values_in, qturns, values)

    deallocate(values_in)

  endif

  call handle%close()

end subroutine read_netcdf_data_2d

!> Read a 2d region array from file using infrastructure I/O.

subroutine mom_read_data_2d_region(filename, fieldname, data, start, nread, MOM_domain, &

                                   no_domain, scale, turns)

  character(len=*), intent(in)  :: filename   !< Input filename

  character(len=*), intent(in)  :: fieldname  !< Field variable name

  real, dimension(:,:), intent(inout) :: data !< Field value in arbitrary units [A ~> a]

  integer, dimension(:), intent(in) :: start  !< Starting index for each axis.

                                              !! In 2d, start(3:4) must be 1.

  integer, dimension(:), intent(in) :: nread  !< Number of values to read along each axis.

                                              !! In 2d, nread(3:4) must be 1.

  type(mom_domain_type), optional, intent(in) :: MOM_Domain !< Model domain decomposition

  logical, optional, intent(in) :: no_domain  !< If true, field does not use

                                              !! domain decomposion.

  real, optional, intent(in)    :: scale      !< A scaling factor that the variable is multiplied by

                                              !! before it is returned to convert from the units in the file

                                              !! to the internal units for this variable [A a-1 ~> 1]

  integer, optional, intent(in) :: turns      !< Number of quarter turns from

                                              !! input to model grid

  integer :: qturns                   ! Number of quarter turns

  real, allocatable :: data_in(:,:)   ! Field array on the input grid in arbitrary units [A ~> a]

  qturns = 0

  if (present(turns)) qturns = modulo(turns, 4)

  if (qturns == 0) then

    call read_field(filename, fieldname, data, start, nread, &

                    mom_domain=mom_domain, no_domain=no_domain, scale=scale)

  else

    call allocate_rotated_array(data, [1,1], -qturns, data_in)

    call rotate_array(data, -qturns, data_in)

    if (associated(mom_domain%domain_in)) then

      call read_field(filename, fieldname, data_in, start, nread, &

                      mom_domain=mom_domain%domain_in, no_domain=no_domain, scale=scale)

    else

      call read_field(filename, fieldname, data_in, start, nread, &

                      mom_domain=mom_domain, no_domain=no_domain, scale=scale)

    endif

    call rotate_array(data_in, qturns, data)

    deallocate(data_in)

  endif

end subroutine mom_read_data_2d_region

!> Read a 3d array from file using infrastructure I/O.

subroutine mom_read_data_3d(filename, fieldname, data, MOM_Domain, timelevel, position, &

                            scale, global_file, file_may_be_4d, turns)

  character(len=*), intent(in)  :: filename     !< Input filename

  character(len=*), intent(in)  :: fieldname    !< Field variable name

  real, dimension(:,:,:), intent(inout) :: data !< Field value in arbitrary units [A ~> a]

  type(mom_domain_type), target, &

                     intent(in) :: MOM_Domain   !< Model domain decomposition

  integer, optional, intent(in) :: timelevel    !< Time level to read in file

  integer, optional, intent(in) :: position     !< Grid positioning flag

  real, optional, intent(in)    :: scale        !< A scaling factor that the variable is multiplied by

                                                !! before it is returned to convert from the units in the file

                                                !! to the internal units for this variable [A a-1 ~> 1]

  logical, optional, intent(in) :: global_file  !< If true, read from a single file

  logical, optional, intent(in) :: file_may_be_4d !< If true, fields may be stored

                                                  !! as 4d arrays in the file.

  integer, optional, intent(in) :: turns        !< Number of quarter-turns to rotate the data.  If absent

                                                !! the number of turns is taken from MOM_Domain.

  ! Local variables

  integer :: qturns   ! Number of quarter-turns from input to model grid

  real, allocatable :: data_in(:,:,:)  ! Field array on the input grid in arbitrary units [A ~> a]

  type(mom_domain_type), pointer :: domain_ptr => null()  ! Pointer to the unrotated domain for reading

  domain_ptr => mom_domain

  if (associated(mom_domain%domain_in) .and. (qturns /= 0)) domain_ptr => mom_domain%domain_in

  qturns = mom_domain%turns ; if (present(turns)) qturns = modulo(turns, 4)

  if (qturns == 0) then

    call read_field(filename, fieldname, data, mom_domain, &

                    timelevel=timelevel, position=position, scale=scale, &

                    global_file=global_file, file_may_be_4d=file_may_be_4d)

  else

    call allocate_rotated_array(data, [1,1,1], -qturns, data_in)

    call rotate_array(data, -qturns, data_in)

    call read_field(filename, fieldname, data_in, domain_ptr, &

                    timelevel=timelevel, position=position, scale=scale, &

                    global_file=global_file, file_may_be_4d=file_may_be_4d)

    call rotate_array(data_in, qturns, data)

    deallocate(data_in)

  endif

end subroutine mom_read_data_3d

!> Read a 3d region array from file using infrastructure I/O.

subroutine mom_read_data_3d_region(filename, fieldname, data, start, nread, MOM_domain, &

                                   no_domain, scale, turns)

  character(len=*), intent(in)  :: filename   !< Input filename

  character(len=*), intent(in)  :: fieldname  !< Field variable name

  real, dimension(:,:,:), intent(inout) :: data !< Field value in arbitrary units [A ~> a]

  integer, dimension(:), intent(in) :: start  !< Starting index for each axis.

  integer, dimension(:), intent(in) :: nread  !< Number of values to read along each axis.

  type(mom_domain_type), optional, intent(in) :: MOM_Domain !< Model domain decomposition

  logical, optional, intent(in) :: no_domain  !< If true, field does not use

                                              !! domain decomposion.

  real, optional, intent(in)    :: scale      !< A scaling factor that the variable is multiplied by

                                              !! before it is returned to convert from the units in the file

                                              !! to the internal units for this variable [A a-1 ~> 1]

  integer, optional, intent(in) :: turns      !< Number of quarter turns from

                                              !! input to model grid

  integer :: qturns                   ! Number of quarter turns

  real, allocatable :: data_in(:,:,:)   ! Field array on the input grid in arbitrary units [A ~> a]

  qturns = 0

  if (present(turns)) qturns = modulo(turns, 4)

  if (qturns == 0) then

    call read_field(filename, fieldname, data, start, nread, &

                    mom_domain=mom_domain, no_domain=no_domain, scale=scale)

  else

    call allocate_rotated_array(data, [1,1,1], -qturns, data_in)

    call rotate_array(data, -qturns, data_in)

    if (associated(mom_domain%domain_in)) then

      call read_field(filename, fieldname, data_in, start, nread, &

                      mom_domain=mom_domain%domain_in, no_domain=no_domain, scale=scale)

    else

      call read_field(filename, fieldname, data_in, start, nread, &

                      mom_domain=mom_domain, no_domain=no_domain, scale=scale)

    endif

    call rotate_array(data_in, qturns, data)

    deallocate(data_in)

  endif

end subroutine mom_read_data_3d_region

!> Read a 4d array from file using infrastructure I/O.

subroutine mom_read_data_4d(filename, fieldname, data, MOM_Domain, &

                            timelevel, position, scale, global_file, turns)

  character(len=*), intent(in) :: filename      !< Input filename

  character(len=*), intent(in) :: fieldname     !< Field variable name

  real, dimension(:,:,:,:), intent(inout) :: data !< Field value in arbitrary units [A ~> a]

  type(mom_domain_type), target, &

                     intent(in) :: MOM_Domain   !< Model domain decomposition

  integer, optional, intent(in) :: timelevel    !< Time level to read in file

  integer, optional, intent(in) :: position     !< Grid positioning flag

  real,    optional, intent(in) :: scale        !< A scaling factor that the variable is multiplied by

                                                !! before it is returned to convert from the units in the file

                                                !! to the internal units for this variable [A a-1 ~> 1]

  logical, optional, intent(in) :: global_file  !< If true, read from a single file

  integer, optional, intent(in) :: turns        !< Number of quarter-turns to rotate the data.  If absent

                                                !! the number of turns is taken from MOM_Domain.

  ! Local variables

  integer :: qturns   ! Number of quarter-turns from input to model grid

  real, allocatable :: data_in(:,:,:,:)  ! Field array on the input grid in arbitrary units [A ~> a]

  type(mom_domain_type), pointer :: domain_ptr => null()  ! Pointer to the unrotated domain for reading

  qturns = mom_domain%turns ; if (present(turns)) qturns = modulo(turns, 4)

  domain_ptr => mom_domain

  if (associated(mom_domain%domain_in) .and. (qturns /= 0)) domain_ptr => mom_domain%domain_in

  if (qturns == 0) then

    call read_field(filename, fieldname, data, mom_domain, &

        timelevel=timelevel, position=position, scale=scale, &

        global_file=global_file)

  else

    ! Read field along the input grid and rotate to the model grid

    call allocate_rotated_array(data, [1,1,1,1], -qturns, data_in)

    call rotate_array(data, -qturns, data_in)

    call read_field(filename, fieldname, data_in, domain_ptr, timelevel=timelevel, &

                    position=position, scale=scale, global_file=global_file)

    call rotate_array(data_in, qturns, data)

    deallocate(data_in)

  endif

end subroutine mom_read_data_4d

!> Read a 2d vector tuple from file using infrastructure I/O.

subroutine mom_read_vector_2d(filename, u_fieldname, v_fieldname, u_data, v_data, MOM_Domain, &

                              timelevel, stagger, scalar_pair, scale, turns)

  character(len=*), intent(in) :: filename      !< Input filename

  character(len=*), intent(in) :: u_fieldname   !< Field variable name in u

  character(len=*), intent(in) :: v_fieldname   !< Field variable name in v

  real, dimension(:,:), intent(inout) :: u_data !< Field value at u points in arbitrary units [A ~> a]

  real, dimension(:,:), intent(inout) :: v_data !< Field value at v points in arbitrary units  [A ~> a]

  type(mom_domain_type), target, &

                     intent(in) :: MOM_Domain   !< Model domain decomposition

  integer, optional, intent(in) :: timelevel    !< Time level to read in file

  integer, optional, intent(in) :: stagger      !< Grid staggering flag

  logical, optional, intent(in) :: scalar_pair  !< True if tuple is not a vector

  real,    optional, intent(in) :: scale        !< A scaling factor that the vector is multiplied by

                                                !! before it is returned to convert from the units in the file

                                                !! to the internal units for this variable [A a-1 ~> 1]

  integer, optional, intent(in) :: turns        !< Number of quarter-turns to rotate the data.  If absent

                                                !! the number of turns is taken from MOM_Domain.

  ! Local variables

  integer :: qturns ! Number of quarter-turns from input to model grid

  real, allocatable :: u_data_in(:,:), v_data_in(:,:)   ! [uv] on the input grid in arbitrary units [A ~> a]

  type(mom_domain_type), pointer :: domain_ptr => null()  ! Pointer to the unrotated domain for reading

  qturns = mom_domain%turns ; if (present(turns)) qturns = modulo(turns, 4)

  domain_ptr => mom_domain

  if (associated(mom_domain%domain_in) .and. (qturns /= 0)) domain_ptr => mom_domain%domain_in

  if (qturns == 0) then

    call read_vector(filename, u_fieldname, v_fieldname, &

                     u_data, v_data, mom_domain, timelevel=timelevel, stagger=stagger, &

                     scalar_pair=scalar_pair, scale=scale)

  else

    call allocate_rotated_array(u_data, [1,1], -qturns, u_data_in)

    call allocate_rotated_array(v_data, [1,1], -qturns, v_data_in)

    if (scalar_pair) then

      call rotate_array_pair(u_data, v_data, -qturns, u_data_in, v_data_in)

    else

      call rotate_vector(u_data, v_data, -qturns, u_data_in, v_data_in)

    endif

    call read_vector(filename, u_fieldname, v_fieldname, u_data_in, v_data_in, &

                     domain_ptr, timelevel=timelevel, &

                     stagger=stagger, scalar_pair=scalar_pair, scale=scale)

    if (scalar_pair) then

      call rotate_array_pair(u_data_in, v_data_in, qturns, u_data, v_data)

    else

      call rotate_vector(u_data_in, v_data_in, qturns, u_data, v_data)

    endif

    deallocate(v_data_in)

    deallocate(u_data_in)

  endif

end subroutine mom_read_vector_2d

!> Read a 3d vector tuple from file using infrastructure I/O.

subroutine mom_read_vector_3d(filename, u_fieldname, v_fieldname, u_data, v_data, MOM_Domain, &

                              timelevel, stagger, scalar_pair, scale, turns)

  character(len=*), intent(in) :: filename      !< Input filename

  character(len=*), intent(in) :: u_fieldname   !< Field variable name in u

  character(len=*), intent(in) :: v_fieldname   !< Field variable name in v

  real, dimension(:,:,:), intent(inout) :: u_data !< Field value in u in arbitrary units [A ~> a]

  real, dimension(:,:,:), intent(inout) :: v_data !< Field value in v in arbitrary units [A ~> a]

  type(mom_domain_type), target, &

                     intent(in) :: MOM_Domain   !< Model domain decomposition

  integer, optional, intent(in) :: timelevel    !< Time level to read in file

  integer, optional, intent(in) :: stagger      !< Grid staggering flag

  logical, optional, intent(in) :: scalar_pair  !< True if tuple is not a vector

  real,    optional, intent(in) :: scale        !< A scaling factor that the vector is multiplied by

                                                !! before it is returned to convert from the units in the file

                                                !! to the internal units for this variable [A a-1 ~> 1]

  integer, optional, intent(in) :: turns        !< Number of quarter-turns to rotate the data.  If absent

                                                !! the number of turns is taken from MOM_Domain.

  ! Local variables

  integer :: qturns ! Number of quarter-turns from input to model grid

  real, allocatable :: u_data_in(:,:,:), v_data_in(:,:,:) ! [uv] on the input grid in arbitrary units [A ~> a]

  type(mom_domain_type), pointer :: domain_ptr => null()  ! Pointer to the unrotated domain for reading

  qturns = mom_domain%turns ; if (present(turns)) qturns = modulo(turns, 4)

  domain_ptr => mom_domain

  if (associated(mom_domain%domain_in) .and. (qturns /= 0)) domain_ptr => mom_domain%domain_in

  if (qturns == 0) then

    call read_vector(filename, u_fieldname, v_fieldname, &

                     u_data, v_data, mom_domain, timelevel=timelevel, stagger=stagger, &

                     scalar_pair=scalar_pair, scale=scale)

  else

    call allocate_rotated_array(u_data, [1,1,1], -qturns, u_data_in)

    call allocate_rotated_array(v_data, [1,1,1], -qturns, v_data_in)

    if (scalar_pair) then

      call rotate_array_pair(u_data, v_data, -qturns, u_data_in, v_data_in)

    else

      call rotate_vector(u_data, v_data, -qturns, u_data_in, v_data_in)

    endif

    call read_vector(filename, u_fieldname, v_fieldname, u_data_in, v_data_in, &

                     domain_ptr, timelevel=timelevel, &

                     stagger=stagger, scalar_pair=scalar_pair, scale=scale)

    if (scalar_pair) then

      call rotate_array_pair(u_data_in, v_data_in, qturns, u_data, v_data)

    else

      call rotate_vector(u_data_in, v_data_in, qturns, u_data, v_data)

    endif

    deallocate(v_data_in)

    deallocate(u_data_in)

  endif

end subroutine mom_read_vector_3d

!> Write a 4d field to an output file, potentially with rotation

subroutine mom_write_field_legacy_4d(IO_handle, field_md, MOM_domain, field, tstamp, tile_count, &

                              fill_value, turns, scale, unscale, zero_zeros)

  type(file_type),          intent(inout) :: IO_handle  !< Handle for a file that is open for writing

  type(fieldtype),          intent(in)    :: field_md   !< Field type with metadata

  type(mom_domain_type),    intent(in)    :: MOM_domain !< The MOM_Domain that describes the decomposition

  real, dimension(:,:,:,:), intent(inout) :: field      !< Unrotated field to write in arbitrary units [A ~> a]

  real,           optional, intent(in)    :: tstamp     !< Model timestamp, often in [days]

  integer,        optional, intent(in)    :: tile_count !< PEs per tile (default: 1)

  real,           optional, intent(in)    :: fill_value !< Missing data fill value in the units used in the file [a]

  integer,        optional, intent(in)    :: turns      !< Number of quarter-turns to rotate the data

  real,           optional, intent(in)    :: scale      !< A scaling factor that the field is multiplied by before

                                                        !! it is written [a A-1 ~> 1], for example to convert it

                                                        !! from its internal units to the desired units for output

  real,           optional, intent(in)    :: unscale    !< A scaling factor that the field is multiplied by before

                                                        !! it is written [a A-1 ~> 1], for example to convert it

                                                        !! from its internal units to the desired units for output.

                                                        !! Here scale and unscale are synonymous, but unscale

                                                        !! takes precedence if both are present.

  logical,        optional, intent(in)    :: zero_zeros !< If present and true, convert negative zeros

                                                        !! into ordinary signless zeros.

  ! Local variables

  real, allocatable :: field_rot(:,:,:,:)  ! A rotated version of field, with the same units [a] or

                                           ! rescaled [A ~> a] then [a]

  real :: scale_fac ! A scaling factor to use before writing the array [a A-1 ~> 1]

  integer :: qturns ! The number of quarter turns through which to rotate field

  qturns = 0 ; if (present(turns)) qturns = modulo(turns, 4)

  scale_fac = 1.0 ; if (present(scale)) scale_fac = scale

  if (present(unscale)) scale_fac = unscale

  if ((qturns == 0) .and. (scale_fac == 1.0) .and. .not.present(zero_zeros)) then

    call write_field(io_handle, field_md, mom_domain, field, tstamp=tstamp, &

                         tile_count=tile_count, fill_value=fill_value)

  else

    call allocate_rotated_array(field, [1,1,1,1], qturns, field_rot)

    call rotate_array(field, qturns, field_rot)

    call rescale_comp_data(mom_domain, field_rot, scale_fac, zero_zeros)

    call write_field(io_handle, field_md, mom_domain, field_rot, tstamp=tstamp, &

                         tile_count=tile_count, fill_value=fill_value)

    deallocate(field_rot)

  endif

end subroutine mom_write_field_legacy_4d

!> Write a 3d field to an output file, potentially with rotation

subroutine mom_write_field_legacy_3d(IO_handle, field_md, MOM_domain, field, tstamp, tile_count, &

                              fill_value, turns, scale, unscale, zero_zeros)

  type(file_type),        intent(inout) :: IO_handle  !< Handle for a file that is open for writing

  type(fieldtype),        intent(in)    :: field_md   !< Field type with metadata

  type(mom_domain_type),  intent(in)    :: MOM_domain !< The MOM_Domain that describes the decomposition

  real, dimension(:,:,:), intent(inout) :: field      !< Unrotated field to write in arbitrary units [A ~> a]

  real,         optional, intent(in)    :: tstamp     !< Model timestamp, often in [days]

  integer,      optional, intent(in)    :: tile_count !< PEs per tile (default: 1)

  real,         optional, intent(in)    :: fill_value !< Missing data fill value in the units used in the file [a]

  integer,      optional, intent(in)    :: turns      !< Number of quarter-turns to rotate the data

  real,         optional, intent(in)    :: scale      !< A scaling factor that the field is multiplied by before

                                                      !! it is written [a A-1 ~> 1], for example to convert it

                                                      !! from its internal units to the desired units for output

  real,         optional, intent(in)    :: unscale    !< A scaling factor that the field is multiplied by before

                                                      !! it is written [a A-1 ~> 1], for example to convert it

                                                      !! from its internal units to the desired units for output.

                                                      !! Here scale and unscale are synonymous, but unscale

                                                      !! takes precedence if both are present.

  logical,      optional, intent(in)    :: zero_zeros !< If present and true, convert negative zeros

                                                      !! into ordinary signless zeros.

  ! Local variables

  real, allocatable :: field_rot(:,:,:)  ! A rotated version of field, with the same units [a] or

                                         ! rescaled [A ~> a] then [a]

  real :: scale_fac ! A scaling factor to use before writing the array [a A-1 ~> 1]

  integer :: qturns ! The number of quarter turns through which to rotate field

  qturns = 0 ; if (present(turns)) qturns = modulo(turns, 4)

  scale_fac = 1.0 ; if (present(scale)) scale_fac = scale

  if (present(unscale)) scale_fac = unscale

  if ((qturns == 0) .and. (scale_fac == 1.0) .and. .not.present(zero_zeros)) then

    call write_field(io_handle, field_md, mom_domain, field, tstamp=tstamp, &

                         tile_count=tile_count, fill_value=fill_value)

  else

    call allocate_rotated_array(field, [1,1,1], qturns, field_rot)

    call rotate_array(field, qturns, field_rot)

    call rescale_comp_data(mom_domain, field_rot, scale_fac, zero_zeros)

    call write_field(io_handle, field_md, mom_domain, field_rot, tstamp=tstamp, &

                         tile_count=tile_count, fill_value=fill_value)

    deallocate(field_rot)

  endif

end subroutine mom_write_field_legacy_3d

!> Write a 2d field to an output file, potentially with rotation

subroutine mom_write_field_legacy_2d(IO_handle, field_md, MOM_domain, field, tstamp, tile_count, &

                              fill_value, turns, scale, unscale, zero_zeros)

  type(file_type),        intent(inout) :: IO_handle  !< Handle for a file that is open for writing

  type(fieldtype),        intent(in)    :: field_md   !< Field type with metadata

  type(mom_domain_type),  intent(in)    :: MOM_domain !< The MOM_Domain that describes the decomposition

  real, dimension(:,:),   intent(inout) :: field      !< Unrotated field to write in arbitrary units [A ~> a]

  real,         optional, intent(in)    :: tstamp     !< Model timestamp, often in [days]

  integer,      optional, intent(in)    :: tile_count !< PEs per tile (default: 1)

  real,         optional, intent(in)    :: fill_value !< Missing data fill value

  integer,      optional, intent(in)    :: turns      !< Number of quarter-turns to rotate the data

  real,         optional, intent(in)    :: scale      !< A scaling factor that the field is multiplied by before

                                                      !! it is written [a A-1 ~> 1], for example to convert it

                                                      !! from its internal units to the desired units for output

  real,         optional, intent(in)    :: unscale    !< A scaling factor that the field is multiplied by before

                                                      !! it is written [a A-1 ~> 1], for example to convert it

                                                      !! from its internal units to the desired units for output.

                                                      !! Here scale and unscale are synonymous, but unscale

                                                      !! takes precedence if both are present.

  logical,      optional, intent(in)    :: zero_zeros !< If present and true, convert negative zeros

                                                      !! into ordinary signless zeros.

  ! Local variables

  real, allocatable :: field_rot(:,:)  ! A rotated version of field, with the same units [a] or

                                       ! rescaled [A ~> a] then [a]

  real :: scale_fac ! A scaling factor to use before writing the array [a A-1 ~> 1]

  integer :: qturns ! The number of quarter turns through which to rotate field

  qturns = 0 ; if (present(turns)) qturns = modulo(turns, 4)

  scale_fac = 1.0 ; if (present(scale)) scale_fac = scale

  if (present(unscale)) scale_fac = unscale

  if ((qturns == 0) .and. (scale_fac == 1.0) .and. .not.present(zero_zeros)) then

    call write_field(io_handle, field_md, mom_domain, field, tstamp=tstamp, &

                         tile_count=tile_count, fill_value=fill_value)

  else

    call allocate_rotated_array(field, [1,1], qturns, field_rot)

    call rotate_array(field, qturns, field_rot)

    call rescale_comp_data(mom_domain, field_rot, scale_fac, zero_zeros)

    call write_field(io_handle, field_md, mom_domain, field_rot, tstamp=tstamp, &

                         tile_count=tile_count, fill_value=fill_value)

    deallocate(field_rot)

  endif

end subroutine mom_write_field_legacy_2d

!> Write a 1d field to an output file

subroutine mom_write_field_legacy_1d(IO_handle, field_md, field, tstamp, fill_value, scale, unscale, zero_zeros)

  type(file_type),        intent(inout) :: IO_handle  !< Handle for a file that is open for writing

  type(fieldtype),        intent(in)    :: field_md   !< Field type with metadata

  real, dimension(:),     intent(in)    :: field      !< Field to write in arbitrary units [A ~> a]

  real,         optional, intent(in)    :: tstamp     !< Model timestamp, often in [days]

  real,         optional, intent(in)    :: fill_value !< Missing data fill value [a]

  real,         optional, intent(in)    :: scale      !< A scaling factor that the field is multiplied by before

                                                      !! it is written [a A-1 ~> 1], for example to convert it

                                                      !! from its internal units to the desired units for output

  real,         optional, intent(in)    :: unscale    !< A scaling factor that the field is multiplied by before

                                                      !! it is written [a A-1 ~> 1], for example to convert it

                                                      !! from its internal units to the desired units for output.

                                                      !! Here scale and unscale are synonymous, but unscale

                                                      !! takes precedence if both are present.

  logical,      optional, intent(in)    :: zero_zeros !< If present and true, convert negative zeros

                                                      !! into ordinary signless zeros.

  ! Local variables

  real, dimension(:), allocatable :: array ! A rescaled copy of field [a]

  real :: scale_fac ! A scaling factor to use before writing the array [a A-1 ~> 1]

  logical :: design_zeros ! If true, convert negative zeros into ordinary signless zeros.

  integer :: i

  scale_fac = 1.0 ; if (present(scale)) scale_fac = scale

  if (present(unscale)) scale_fac = unscale

  design_zeros = .false. ; if (present(zero_zeros)) design_zeros = zero_zeros

  if ((scale_fac == 1.0) .and. (.not.design_zeros)) then

    call write_field(io_handle, field_md, field, tstamp=tstamp)

  else

    allocate(array(size(field)))

    array(:) = scale_fac * field(:)

    if (present(fill_value)) then

      do i=1,size(field) ; if (field(i) == fill_value) array(i) = fill_value ; enddo

    endif

    if (design_zeros) then ! Convert negative zeros into zeros

      do i=1,size(field) ; if (array(i) == 0.0) array(i) = 0.0 ; enddo

    endif

    call write_field(io_handle, field_md, array, tstamp=tstamp)

    deallocate(array)

  endif

end subroutine mom_write_field_legacy_1d

!> Write a 0d field to an output file

subroutine mom_write_field_legacy_0d(IO_handle, field_md, field, tstamp, fill_value, scale, unscale, zero_zeros)

  type(file_type),        intent(inout) :: IO_handle  !< Handle for a file that is open for writing

  type(fieldtype),        intent(in)    :: field_md   !< Field type with metadata

  real,                   intent(in)    :: field      !< Field to write in arbitrary units [A ~> a]

  real,         optional, intent(in)    :: tstamp     !< Model timestamp, often in [days]

  real,         optional, intent(in)    :: fill_value !< Missing data fill value [a]

  real,         optional, intent(in)    :: scale      !< A scaling factor that the field is multiplied by before

                                                      !! it is written [a A-1 ~> 1], for example to convert it

                                                      !! from its internal units to the desired units for output

  real,         optional, intent(in)    :: unscale    !< A scaling factor that the field is multiplied by before

                                                      !! it is written [a A-1 ~> 1], for example to convert it

                                                      !! from its internal units to the desired units for output.

                                                      !! Here scale and unscale are synonymous, but unscale

                                                      !! takes precedence if both are present.

  logical,      optional, intent(in)    :: zero_zeros !< If present and true, convert negative zeros

                                                      !! into ordinary signless zeros.

  ! Local variables

  real :: scale_fac  ! A scaling factor to use before writing the field [a A-1 ~> 1]

  real :: scaled_val ! A rescaled copy of field [a]

  scale_fac = 1.0 ; if (present(scale)) scale_fac = scale

  if (present(unscale)) scale_fac = unscale

  scaled_val = field * scale_fac

  if (present(fill_value)) then ; if (field == fill_value) scaled_val = fill_value ; endif

  if (present(zero_zeros)) then ; if (zero_zeros .and. (scaled_val == 0.0)) scaled_val = 0.0 ; endif

  call write_field(io_handle, field_md, scaled_val, tstamp=tstamp)

end subroutine mom_write_field_legacy_0d

!> Write a 4d field to an output file, potentially with rotation

subroutine mom_write_field_4d(IO_handle, field_md, MOM_domain, field, tstamp, tile_count, &

                              fill_value, turns, scale, unscale, zero_zeros)

  class(mom_file),          intent(inout) :: IO_handle  !< Handle for a file that is open for writing

  type(mom_field),          intent(in)    :: field_md   !< Field type with metadata

  type(mom_domain_type),    intent(in)    :: MOM_domain !< The MOM_Domain that describes the decomposition

  real, dimension(:,:,:,:), intent(inout) :: field      !< Unrotated field to write in arbitrary units [A ~> a]

  real,           optional, intent(in)    :: tstamp     !< Model timestamp, often in [days]

  integer,        optional, intent(in)    :: tile_count !< PEs per tile (default: 1)

  real,           optional, intent(in)    :: fill_value !< Missing data fill value [a]

  integer,        optional, intent(in)    :: turns      !< Number of quarter-turns to rotate the data

  real,           optional, intent(in)    :: scale      !< A scaling factor that the field is multiplied by before

                                                        !! it is written [a A-1 ~> 1], for example to convert it

                                                        !! from its internal units to the desired units for output

  real,           optional, intent(in)    :: unscale    !< A scaling factor that the field is multiplied by before

                                                        !! it is written [a A-1 ~> 1], for example to convert it

                                                        !! from its internal units to the desired units for output.

                                                        !! Here scale and unscale are synonymous, but unscale

                                                        !! takes precedence if both are present.

  logical,        optional, intent(in)    :: zero_zeros !< If present and true, convert negative zeros

                                                        !! into ordinary signless zeros.

  ! Local variables

  real, allocatable :: field_rot(:,:,:,:)  ! A rotated version of field, with the same units or rescaled [a]

  real :: scale_fac ! A scaling factor to use before writing the array [a A-1 ~> 1]

  integer :: qturns ! The number of quarter turns through which to rotate field

  qturns = 0 ; if (present(turns)) qturns = modulo(turns, 4)

  scale_fac = 1.0 ; if (present(scale)) scale_fac = scale

  if (present(unscale)) scale_fac = unscale

  if ((qturns == 0) .and. (scale_fac == 1.0) .and. .not.present(zero_zeros)) then

    call io_handle%write_field(field_md, mom_domain, field, tstamp=tstamp, &

                         tile_count=tile_count, fill_value=fill_value)

  else

    call allocate_rotated_array(field, [1,1,1,1], qturns, field_rot)

    call rotate_array(field, qturns, field_rot)

    call rescale_comp_data(mom_domain, field_rot, scale_fac, zero_zeros)

    call io_handle%write_field(field_md, mom_domain, field_rot, tstamp=tstamp, &

                         tile_count=tile_count, fill_value=fill_value)

    deallocate(field_rot)

  endif

end subroutine mom_write_field_4d

!> Write a 3d field to an output file, potentially with rotation

subroutine mom_write_field_3d(IO_handle, field_md, MOM_domain, field, tstamp, tile_count, &

                              fill_value, turns, scale, unscale, zero_zeros)

  class(mom_file),        intent(inout) :: IO_handle  !< Handle for a file that is open for writing

  type(mom_field),        intent(in)    :: field_md   !< Field type with metadata

  type(mom_domain_type),  intent(in)    :: MOM_domain !< The MOM_Domain that describes the decomposition

  real, dimension(:,:,:), intent(inout) :: field      !< Unrotated field to write in arbitrary units [A ~> a]

  real,         optional, intent(in)    :: tstamp     !< Model timestamp, often in [days]

  integer,      optional, intent(in)    :: tile_count !< PEs per tile (default: 1)

  real,         optional, intent(in)    :: fill_value !< Missing data fill value [a]

  integer,      optional, intent(in)    :: turns      !< Number of quarter-turns to rotate the data

  real,         optional, intent(in)    :: scale      !< A scaling factor that the field is multiplied by before

                                                      !! it is written [a A-1 ~> 1], for example to convert it

                                                      !! from its internal units to the desired units for output

  real,         optional, intent(in)    :: unscale    !< A scaling factor that the field is multiplied by before

                                                      !! it is written [a A-1 ~> 1], for example to convert it

                                                      !! from its internal units to the desired units for output.

                                                      !! Here scale and unscale are synonymous, but unscale

                                                      !! takes precedence if both are present.

  logical,      optional, intent(in)    :: zero_zeros !< If present and true, convert negative zeros

                                                      !! into ordinary signless zeros.

  ! Local variables

  real, allocatable :: field_rot(:,:,:)  ! A rotated version of field, with the same units or rescaled [a]

  real :: scale_fac ! A scaling factor to use before writing the array [a A-1 ~> 1]

  integer :: qturns ! The number of quarter turns through which to rotate field

  qturns = 0 ; if (present(turns)) qturns = modulo(turns, 4)

  scale_fac = 1.0 ; if (present(scale)) scale_fac = scale

  if (present(unscale)) scale_fac = unscale

  if ((qturns == 0) .and. (scale_fac == 1.0) .and. .not.present(zero_zeros)) then

    call io_handle%write_field(field_md, mom_domain, field, tstamp=tstamp, &

                         tile_count=tile_count, fill_value=fill_value)

  else

    call allocate_rotated_array(field, [1,1,1], qturns, field_rot)

    call rotate_array(field, qturns, field_rot)

    call rescale_comp_data(mom_domain, field_rot, scale_fac, zero_zeros)

    call io_handle%write_field(field_md, mom_domain, field_rot, tstamp=tstamp, &

                         tile_count=tile_count, fill_value=fill_value)

    deallocate(field_rot)

  endif

end subroutine mom_write_field_3d

!> Write a 2d field to an output file, potentially with rotation

subroutine mom_write_field_2d(IO_handle, field_md, MOM_domain, field, tstamp, tile_count, &

                              fill_value, turns, scale, unscale, zero_zeros)

  class(mom_file),        intent(inout) :: IO_handle  !< Handle for a file that is open for writing

  type(mom_field),        intent(in)    :: field_md   !< Field type with metadata

  type(mom_domain_type),  intent(in)    :: MOM_domain !< The MOM_Domain that describes the decomposition

  real, dimension(:,:),   intent(inout) :: field      !< Unrotated field to write in arbitrary units [A ~> a]

  real,         optional, intent(in)    :: tstamp     !< Model timestamp, often in [days]

  integer,      optional, intent(in)    :: tile_count !< PEs per tile (default: 1)

  real,         optional, intent(in)    :: fill_value !< Missing data fill value [a]

  integer,      optional, intent(in)    :: turns      !< Number of quarter-turns to rotate the data

  real,         optional, intent(in)    :: scale      !< A scaling factor that the field is multiplied by before

                                                      !! it is written [a A-1 ~> 1], for example to convert it

                                                      !! from its internal units to the desired units for output

  real,         optional, intent(in)    :: unscale    !< A scaling factor that the field is multiplied by before

                                                      !! it is written [a A-1 ~> 1], for example to convert it

                                                      !! from its internal units to the desired units for output.

                                                      !! Here scale and unscale are synonymous, but unscale

                                                      !! takes precedence if both are present.

  logical,      optional, intent(in)    :: zero_zeros !< If present and true, convert negative zeros

                                                      !! into ordinary signless zeros.

  ! Local variables

  real, allocatable :: field_rot(:,:)  ! A rotated version of field, with the same units or rescaled [a]

  real :: scale_fac ! A scaling factor to use before writing the array [a A-1 ~> 1]

  integer :: qturns ! The number of quarter turns through which to rotate field

  qturns = 0 ; if (present(turns)) qturns = modulo(turns, 4)

  scale_fac = 1.0 ; if (present(scale)) scale_fac = scale

  if (present(unscale)) scale_fac = unscale

  if ((qturns == 0) .and. (scale_fac == 1.0) .and. .not.present(zero_zeros)) then

    call io_handle%write_field(field_md, mom_domain, field, tstamp=tstamp, &

                         tile_count=tile_count, fill_value=fill_value)

  else

    call allocate_rotated_array(field, [1,1], qturns, field_rot)

    call rotate_array(field, qturns, field_rot)

    call rescale_comp_data(mom_domain, field_rot, scale_fac, zero_zeros)

    call io_handle%write_field(field_md, mom_domain, field_rot, tstamp=tstamp, &

                         tile_count=tile_count, fill_value=fill_value)

    deallocate(field_rot)

  endif

end subroutine mom_write_field_2d

!> Write a 1d field to an output file

subroutine mom_write_field_1d(IO_handle, field_md, field, tstamp, fill_value, scale, unscale, zero_zeros)

  class(mom_file),        intent(inout) :: IO_handle  !< Handle for a file that is open for writing

  type(mom_field),        intent(in)    :: field_md   !< Field type with metadata

  real, dimension(:),     intent(in)    :: field      !< Field to write in arbitrary units [A ~> a]

  real,         optional, intent(in)    :: tstamp     !< Model timestamp, often in [days]

  real,         optional, intent(in)    :: fill_value !< Missing data fill value [a]

  real,         optional, intent(in)    :: scale      !< A scaling factor that the field is multiplied by before

                                                      !! it is written [a A-1 ~> 1], for example to convert it

                                                      !! from its internal units to the desired units for output

  real,         optional, intent(in)    :: unscale    !< A scaling factor that the field is multiplied by before

                                                      !! it is written [a A-1 ~> 1], for example to convert it

                                                      !! from its internal units to the desired units for output.

                                                      !! Here scale and unscale are synonymous, but unscale

                                                      !! takes precedence if both are present.

  logical,      optional, intent(in)    :: zero_zeros !< If present and true, convert negative zeros

                                                      !! into ordinary signless zeros.

  ! Local variables

  real, dimension(:), allocatable :: array ! A rescaled copy of field in arbtrary unscaled units [a]

  real :: scale_fac ! A scaling factor to use before writing the array [a A-1 ~> 1]

  logical :: design_zeros ! If true, convert negative zeros into ordinary signless zeros.

  integer :: i

  scale_fac = 1.0 ; if (present(scale)) scale_fac = scale

  if (present(unscale)) scale_fac = unscale

  design_zeros = .false. ; if (present(zero_zeros)) design_zeros = zero_zeros

  if ((scale_fac == 1.0) .and. (.not.design_zeros)) then

    call io_handle%write_field(field_md, field, tstamp=tstamp)

  else

    allocate(array(size(field)))

    array(:) = scale_fac * field(:)

    if (present(fill_value)) then

      do i=1,size(field) ; if (field(i) == fill_value) array(i) = fill_value ; enddo

    endif

    if (design_zeros) then ! Convert negative zeros into zeros

      do i=1,size(field) ; if (array(i) == 0.0) array(i) = 0.0 ; enddo

    endif

    call io_handle%write_field(field_md, array, tstamp=tstamp)

    deallocate(array)

  endif

end subroutine mom_write_field_1d

!> Write a 0d field to an output file

subroutine mom_write_field_0d(IO_handle, field_md, field, tstamp, fill_value, scale, unscale, zero_zeros)

  class(mom_file),        intent(inout) :: IO_handle  !< Handle for a file that is open for writing

  type(mom_field),        intent(in)    :: field_md   !< Field type with metadata

  real,                   intent(in)    :: field      !< Field to write in arbitrary units [A ~> a]

  real,         optional, intent(in)    :: tstamp     !< Model timestamp, often in [days]

  real,         optional, intent(in)    :: fill_value !< Missing data fill value [a]

  real,         optional, intent(in)    :: scale      !< A scaling factor that the field is multiplied by before

                                                      !! it is written [a A-1 ~> 1], for example to convert it

                                                      !! from its internal units to the desired units for output

  real,         optional, intent(in)    :: unscale    !< A scaling factor that the field is multiplied by before

                                                      !! it is written [a A-1 ~> 1], for example to convert it

                                                      !! from its internal units to the desired units for output.

                                                      !! Here scale and unscale are synonymous, but unscale

                                                      !! takes precedence if both are present.

  logical,      optional, intent(in)    :: zero_zeros !< If present and true, convert negative zeros

                                                      !! into ordinary signless zeros.

  ! Local variables

  real :: scale_fac  ! A scaling factor to use before writing the field [a A-1 ~> 1]

  real :: scaled_val ! A rescaled copy of field in arbtrary unscaled units [a]

  scale_fac = 1.0 ; if (present(scale)) scale_fac = scale

  if (present(unscale)) scale_fac = unscale

  scaled_val = field * scale_fac

  if (present(fill_value)) then ; if (field == fill_value) scaled_val = fill_value ; endif

  if (present(zero_zeros)) then ; if (zero_zeros .and. (scaled_val == 0.0)) scaled_val = 0.0 ; endif

  call io_handle%write_field(field_md, scaled_val, tstamp=tstamp)

end subroutine mom_write_field_0d

!> Insert the ensemble appendix into a filename. If provided, the appendix is inserted after

!! the last occurrence of the insert_after substring in the filename.

subroutine insert_ensemble_appendix(filename, insert_after)

  character(len=*),           intent(inout) :: filename     !< The filename to which the appendix is inserted

  character(len=*), optional, intent(in)    :: insert_after !< The string after which the appendix is inserted.

                                                            !! If not provided or found, the appendix is inserted

                                                            !! at the end of the filename.

  ! Local variables

  character(len=32) :: filename_appendix            ! ensemble appendix to be inserted into the filename

  character(len=:), allocatable :: filename_tr      ! trimmed filename

  character(len=:), allocatable :: insert_after_tr  ! trimmed insert_after

  integer :: pos ! The filename string index after which the appendix is to be inserted

  call get_filename_appendix(filename_appendix)

  if (len_trim(filename_appendix) == 0) return

  filename_tr = trim(adjustl(filename))

  pos = len(filename_tr)

  ! If insert_after is provided, find the last occurrence of insert_after in the filename and set pos accordingly.

  if (present(insert_after)) then

    insert_after_tr = trim(adjustl(insert_after))

    pos = index(filename_tr, insert_after_tr, back=.true.)

    if (pos == 0) then

      call mom_error(fatal, "insert_ensemble_appendix: The string " // insert_after_tr // &

                     " was not found in the filename " // filename_tr)

    endif

    pos = pos + len(insert_after_tr) - 1

  endif

  ! Insert the ensemble appendix into the filename. If the appendix is to be added to

  ! the end of the filename, do so before the .nc extension if it exists.

  if (pos>3 .and. pos == len(filename_tr)) then

    if (filename_tr(pos-2:pos) == ".nc") then

      pos = pos - 3 ! Position before the .nc extension

    endif

  endif

  filename = filename_tr(1:pos) // trim(filename_appendix) // filename_tr(pos+1:)

end subroutine insert_ensemble_appendix

!> Given filename and fieldname, this subroutine returns the size of the field in the file

subroutine field_size(filename, fieldname, sizes, field_found, no_domain, ndims, ncid_in)

  character(len=*),      intent(in)    :: filename  !< The name of the file to read

  character(len=*),      intent(in)    :: fieldname !< The name of the variable whose sizes are returned

  integer, dimension(:), intent(inout) :: sizes     !< The sizes of the variable in each dimension

  logical,     optional, intent(out)   :: field_found !< This indicates whether the field was found in

                                                    !! the input file.  Without this argument, there

                                                    !! is a fatal error if the field is not found.

  logical,     optional, intent(in)    :: no_domain !< If present and true, do not check for file

                                                    !! names with an appended tile number.  If

                                                    !! ndims is present, the default changes to true.

  integer,     optional, intent(out)   :: ndims     !< The number of dimensions to the variable

  integer,     optional, intent(in)    :: ncid_in   !< The netCDF ID of an open file.  If absent, the

                                                    !! file is opened and closed within this routine.

  if (present(ndims)) then

    if (present(no_domain)) then ; if (.not.no_domain) call mom_error(fatal, &

          "field_size does not support the ndims argument when no_domain is present and false.")

    endif

    call get_var_sizes(filename, fieldname, ndims, sizes, match_case=.false., ncid_in=ncid_in)

    if (present(field_found)) field_found = (ndims >= 0)

    if ((ndims < 0) .and. .not.present(field_found)) then

      call mom_error(fatal, "Variable "//trim(fieldname)//" not found in "//trim(filename) )

    endif

  else

    call get_field_size(filename, fieldname, sizes, field_found=field_found, no_domain=no_domain)

  endif

end subroutine field_size

!> Copies a string

subroutine safe_string_copy(str1, str2, fieldnm, caller)

  character(len=*),           intent(in)  :: str1    !< The string being copied

  character(len=*),           intent(out) :: str2    !< The string being copied into

  character(len=*), optional, intent(in)  :: fieldnm !< The name of the field for error messages

  character(len=*), optional, intent(in)  :: caller  !< The calling routine for error messages

  if (len(trim(str1)) > len(str2)) then

    if (present(fieldnm) .and. present(caller)) then

      call mom_error(fatal, trim(caller)//" attempted to copy the overly long string "//&

                     trim(str1)//" into "//trim(fieldnm))

    else

      call mom_error(fatal, "safe_string_copy: The string "//trim(str1)//&

                     " is longer than its intended target.")

    endif

  endif

  str2 = trim(str1)

end subroutine safe_string_copy

!> Returns a name with "%#E" or "%E" replaced with the ensemble member number.

function ensembler(name, ens_no_in) result(en_nm)

  character(len=*),  intent(in) :: name       !< The name to be modified

  integer, optional, intent(in) :: ens_no_in  !< The number of the current ensemble member

  character(len=len(name)) :: en_nm  !< The name encoded with the ensemble number

  ! This function replaces "%#E" or "%E" with the ensemble number anywhere it

  ! occurs in name, with %E using 4 or 6 digits (depending on the ensemble size)

  ! and %#E using # digits, where # is a number from 1 to 9.

  character(len=len(name)) :: tmp

  character(10) :: ens_num_char

  character(3)  :: code_str

  integer :: ens_no

  integer :: n, is

  en_nm = trim(name)

  if (index(name,"%") == 0) return

  if (present(ens_no_in)) then

    ens_no = ens_no_in

  else

    ens_no = get_ensemble_id()

  endif

  write(ens_num_char, '(I0)') ens_no

  do

    is = index(en_nm,"%E")

    if (is == 0) exit

    if (len(en_nm) < len(trim(en_nm)) + len(trim(ens_num_char)) - 2) &

      call mom_error(fatal, "MOM_io ensembler: name "//trim(name)// &

      " is not long enough for %E expansion for ens_no "//trim(ens_num_char))

    tmp = en_nm(1:is-1)//trim(ens_num_char)//trim(en_nm(is+2:))

    en_nm = tmp

  enddo

  if (index(name,"%") == 0) return

  write(ens_num_char, '(I10.10)') ens_no

  do n=1,9 ; do

    write(code_str, '("%",I1,"E")') n

    is = index(en_nm,code_str)

    if (is == 0) exit

    if (ens_no < 10**n) then

      if (len(en_nm) < len(trim(en_nm)) + n-3) call mom_error(fatal, &

        "MOM_io ensembler: name "//trim(name)//" is not long enough for %E expansion.")

      tmp = en_nm(1:is-1)//trim(ens_num_char(11-n:10))//trim(en_nm(is+3:))

    else

      call mom_error(fatal, "MOM_io ensembler: Ensemble number is too large "//&

          "to be encoded with "//code_str//" in "//trim(name))

    endif

    en_nm = tmp

  enddo ; enddo

end function ensembler

!> Provide a string to append to filenames, to differentiate ensemble members, for example.

subroutine get_filename_appendix(suffix)

  character(len=*), intent(out) :: suffix !< A string to append to filenames

  call get_filename_suffix(suffix)

end subroutine get_filename_appendix

!> Write a file version number to the log file or other output file

subroutine write_version_number(version, tag, unit)

  character(len=*),           intent(in) :: version !< A string that contains the routine name and version

  character(len=*), optional, intent(in) :: tag  !< A tag name to add to the message

  integer,          optional, intent(in) :: unit !< An alternate unit number for output

  call write_version(version, tag, unit)

end subroutine write_version_number

!> Open a single namelist file that is potentially readable by all PEs.

function open_namelist_file(file) result(unit)

  character(len=*), optional, intent(in) :: file !< The file to open, by default "input.nml"

  integer                                :: unit !< The opened unit number of the namelist file

  unit = mom_namelist_file(file)

end function open_namelist_file

!> Checks the iostat argument that is returned after reading a namelist variable and writes a

!! message if there is an error.

function check_nml_error(IOstat, nml_name) result(ierr)

  integer,          intent(in) :: iostat   !< An I/O status field from a namelist read call

  character(len=*), intent(in) :: nml_name !< The name of the namelist

  integer :: ierr    !< A copy of IOstat that is returned to preserve legacy function behavior

  call check_namelist_error(iostat, nml_name)

  ierr = iostat

end function check_nml_error

!> Initialize the MOM_io module

subroutine mom_io_init(param_file)

  type(param_file_type), intent(in) :: param_file  !< structure indicating the open file to

                                                   !! parse for model parameter values.

  ! This include declares and sets the variable "version".

# include "version_variable.h"

  character(len=40)  :: mdl = "MOM_io" ! This module's name.

  call log_version(param_file, mdl, version)

end subroutine mom_io_init

!> Returns the dimension variable information for a netCDF variable

subroutine get_var_axes_info(filename, fieldname, axes_info)

  character(len=*), intent(in) ::            filename  !< A filename from which to read

  character(len=*), intent(in) ::            fieldname !< The name of the field to read

  type(axis_info), dimension(4), intent(inout) :: axes_info !< A returned array of field axis information

  !! local variables

  integer ::  rcode

  logical :: success

  integer ::  ncid, varid, ndims

  integer :: id, jd, kd

  integer, dimension(4) :: dims, dim_id

  character(len=128)  :: dim_name(4)

  integer, dimension(1) :: start, count

  !! cartesian axis data

  real, allocatable, dimension(:) :: x ! x-axis labels, often [degrees_E] or [km] or [m]

  real, allocatable, dimension(:) :: y ! y-axis labels, often [degrees_N] or [km] or [m]

  real, allocatable, dimension(:) :: z ! vertical axis labels [various], often [m] or [kg m-3]

  call open_file_to_read(filename, ncid, success=success)

  rcode = nf90_inq_varid(ncid, trim(fieldname), varid)

  if (rcode /= 0) call mom_error(fatal,"error finding variable "//trim(fieldname)//&

                                 " in file "//trim(filename)//" in hinterp_extrap")

  rcode = nf90_inquire_variable(ncid, varid, ndims=ndims, dimids=dims)

  if (rcode /= 0) call mom_error(fatal, "Error inquiring about the dimensions of "//trim(fieldname)//&

                                 " in file "//trim(filename)//" in hinterp_extrap")

  if (ndims < 3) call mom_error(fatal,"Variable "//trim(fieldname)//" in file "//trim(filename)// &

                                " has too few dimensions to be read as a 3-d array.")

  rcode = nf90_inquire_dimension(ncid, dims(1), dim_name(1), len=id)

  if (rcode /= 0) call mom_error(fatal,"error reading dimension 1 data for "// &

                trim(fieldname)//" in file "// trim(filename)//" in hinterp_extrap")

  rcode = nf90_inq_varid(ncid, dim_name(1), dim_id(1))

  if (rcode /= 0) call mom_error(fatal,"error finding variable "//trim(dim_name(1))//&

                                 " in file "//trim(filename)//" in hinterp_extrap")

  rcode = nf90_inquire_dimension(ncid, dims(2), dim_name(2), len=jd)

  if (rcode /= 0) call mom_error(fatal,"error reading dimension 2 data for "// &

                trim(fieldname)//" in file "// trim(filename)//" in hinterp_extrap")

  rcode = nf90_inq_varid(ncid, dim_name(2), dim_id(2))

  if (rcode /= 0) call mom_error(fatal,"error finding variable "//trim(dim_name(2))//&

                                 " in file "//trim(filename)//" in hinterp_extrap")

  rcode = nf90_inquire_dimension(ncid, dims(3), dim_name(3), len=kd)

  if (rcode /= 0) call mom_error(fatal,"error reading dimension 3 data for "// &

                trim(fieldname)//" in file "// trim(filename)//" in hinterp_extrap")

  rcode = nf90_inq_varid(ncid, dim_name(3), dim_id(3))

  if (rcode /= 0) call mom_error(fatal,"error finding variable "//trim(dim_name(3))//&

                                 " in file "//trim(filename)//" in hinterp_extrap")

  allocate(x(id), y(jd), z(kd))

  start = 1 ; count = 1 ; count(1) = id

  rcode = nf90_get_var(ncid, dim_id(1), x, start, count)

  if (rcode /= 0) call mom_error(fatal,"error reading dimension 1 values for var_name "// &

                trim(fieldname)//",dim_name "//trim(dim_name(1))//" in file "// trim(filename)//" in hinterp_extrap")

  start = 1 ; count = 1 ; count(1) = jd

  rcode = nf90_get_var(ncid, dim_id(2), y, start, count)

  if (rcode /= 0) call mom_error(fatal,"error reading dimension 2 values for var_name "// &

                trim(fieldname)//",dim_name "//trim(dim_name(2))//" in file "// trim(filename)//" in  hinterp_extrap")

  start = 1 ; count = 1 ; count(1) = kd

  rcode = nf90_get_var(ncid, dim_id(3), z, start, count)

  if (rcode /= 0) call mom_error(fatal,"error reading dimension 3 values for var_name "// &

                trim(fieldname//",dim_name "//trim(dim_name(3)))//" in file "// trim(filename)//" in  hinterp_extrap")

  call set_axis_info(axes_info(1), name=trim(dim_name(1)), ax_size=id, ax_data=x,cartesian='X')

  call set_axis_info(axes_info(2), name=trim(dim_name(2)), ax_size=jd, ax_data=y,cartesian='Y')

  call set_axis_info(axes_info(3), name=trim(dim_name(3)), ax_size=kd, ax_data=z,cartesian='Z')

  call close_file_to_read(ncid, filename)

  deallocate(x,y,z)

end subroutine get_var_axes_info

!> \namespace mom_io

!!

!!   This file contains a number of subroutines that manipulate

!!  NetCDF files and handle input and output of fields.  These

!!  subroutines, along with their purpose, are:

!!

!!   * create_file: create a new file and set up structures that are

!!       needed for subsequent output and write out the coordinates.

!!   * reopen_file: reopen an existing file for writing and set up

!!       structures that are needed for subsequent output.

!!   * open_input_file: open the indicated file for reading only.

!!   * close_file: close an open file.

!!   * synch_file: flush the buffers, completing all pending output.

!!

!!   * write_field: write a field to an open file.

!!   * write_time: write a value of the time axis to an open file.

!!   * read_data: read a variable from an open file.

!!   * read_time: read a time from an open file.

!!

!!   * name_output_file: provide a name for an output file based on a

!!       name root and the time of the output.

!!   * find_input_file: find a file that has been previously written by

!!       MOM and named by name_output_file and open it for reading.

!!

!!   * handle_error: write an error code and quit.

end module mom_io