MOM_io.F90

!> This module contains I/O framework code
module mom_io
 
! This file is part of MOM6. See LICENSE.md for the license.
 
 
use mom_error_handler,    only : mom_error, note, fatal, warning
use mom_domains,          only : mom_domain_type, agrid, bgrid_ne, cgrid_ne
use mom_domains,          only : get_simple_array_i_ind, get_simple_array_j_ind
use mom_file_parser,      only : log_version, param_file_type
use mom_grid,             only : ocean_grid_type
use mom_dyn_horgrid,      only : dyn_horgrid_type
use mom_string_functions, only : lowercase, slasher
use mom_verticalgrid,     only : verticalgrid_type
 
use ensemble_manager_mod, only : get_ensemble_id
use fms_mod,              only : write_version_number, open_namelist_file, check_nml_error
use fms_io_mod,           only : file_exist, field_size, read_data
use fms_io_mod,           only : field_exists => field_exist, io_infra_end=>fms_io_exit
use fms_io_mod,           only : get_filename_appendix => get_filename_appendix
use mpp_domains_mod,      only : domain1d, domain2d, mpp_get_domain_components
use mpp_domains_mod,      only : center, corner, north_face=>north, east_face=>east
use mpp_io_mod,           only : open_file => mpp_open, close_file => mpp_close
use mpp_io_mod,           only : mpp_write_meta, write_field => mpp_write, mpp_get_info
use mpp_io_mod,           only : mpp_get_atts, mpp_get_axes, get_axis_data=>mpp_get_axis_data, axistype
use mpp_io_mod,           only : mpp_get_fields, fieldtype, axistype, flush_file => mpp_flush
use mpp_io_mod,           only : append_file=>mpp_append, ascii_file=>mpp_ascii
use mpp_io_mod,           only : multiple=>mpp_multi, netcdf_file=>mpp_netcdf
use mpp_io_mod,           only : overwrite_file=>mpp_overwr, readonly_file=>mpp_rdonly
use mpp_io_mod,           only : single_file=>mpp_single, writeonly_file=>mpp_wronly
use mpp_io_mod,           only : mpp_append, mpp_multi, mpp_overwr, mpp_netcdf, mpp_rdonly
use mpp_io_mod,           only : get_file_info=>mpp_get_info, get_file_atts=>mpp_get_atts
use mpp_io_mod,           only : get_file_fields=>mpp_get_fields, get_file_times=>mpp_get_times
use mpp_io_mod,           only : io_infra_init=>mpp_io_init
 
use netcdf
 
implicit none ; private
 
public :: close_file, create_file, field_exists, field_size, fieldtype, get_filename_appendix
public :: file_exists, flush_file, get_file_info, get_file_atts, get_file_fields
public :: get_file_times, open_file, read_axis_data, read_data
public :: num_timelevels, mom_read_data, mom_read_vector, ensembler
public :: reopen_file, slasher, write_field, write_version_number, mom_io_init
public :: open_namelist_file, check_nml_error, io_infra_init, io_infra_end
public :: append_file, ascii_file, multiple, netcdf_file, overwrite_file
public :: readonly_file, single_file, writeonly_file
public :: center, corner, north_face, east_face
public :: var_desc, modify_vardesc, query_vardesc, cmor_long_std
public :: get_axis_data
 
!> Type for describing a variable, typically a tracer
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
end type vardesc
 
!> Indicate whether a file exists, perhaps with domain decomposition
interface file_exists
  module procedure fms_file_exists
  module procedure mom_file_exists
end interface
 
!> Read a data field from a file
interface mom_read_data
  module procedure mom_read_data_4d
  module procedure mom_read_data_3d
  module procedure mom_read_data_2d
  module procedure mom_read_data_1d
end interface
 
!> Read a pair of data fields representing the two components of a vector from a file
interface mom_read_vector
  module procedure mom_read_vector_3d
  module procedure mom_read_vector_2d
end interface
 
contains
 
!> Routine creates a new NetCDF file.  It also sets up
!! structures that describe this file and variables that will
!! later be written to this file. Type for describing a variable, typically a tracer
subroutine create_file(unit, filename, vars, novars, fields, threading, timeunit, G, dG, GV, checksums)
  integer,               intent(out)   :: unit       !< unit id of an open file or -1 on a
                                                     !! nonwriting PE with single file output
  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=8), optional,      intent(in)    :: checksums(:,:)  !< checksums of vars
 
  logical        :: use_lath, use_lonh, use_latq, use_lonq, use_time
  logical        :: use_layer, use_int, use_periodic
  logical        :: one_file, domain_set
  type(axistype) :: axis_lath, axis_latq, axis_lonh, axis_lonq
  type(axistype) :: axis_layer, axis_int, axis_time, axis_periodic
  type(axistype) :: axes(4)
  type(mom_domain_type), pointer :: domain => null()
  type(domain1d) :: x_domain, y_domain
  integer        :: numaxes, pack, thread, k
  integer        :: isg, ieg, jsg, jeg, isgb, iegb, jsgb, jegb
  integer        :: var_periods, num_periods=0
  real, dimension(:), allocatable :: period_val
  real, pointer, dimension(:) :: &
    gridlatt => null(), & ! The latitude or longitude of T or B points for
    gridlatb => null(), & ! the purpose of labeling the output axes.
    gridlont => null(), gridlonb => null()
  character(len=40) :: time_units, x_axis_units, y_axis_units
  character(len=8)  :: t_grid, t_grid_read
 
  use_lath  = .false. ; use_lonh     = .false.
  use_latq  = .false. ; use_lonq     = .false.
  use_time  = .false. ; use_periodic = .false.
  use_layer = .false. ; use_int      = .false.
 
  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
    call open_file(unit, filename, mpp_overwr, mpp_netcdf, threading=thread)
  else
    call open_file(unit, filename, mpp_overwr, mpp_netcdf, domain=domain%mpp_domain)
  endif
 
! Define the coordinates.
  do k=1,novars
    select case (vars(k)%hor_grid)
      case ('h') ; use_lath = .true. ; use_lonh = .true.
      case ('q') ; use_latq = .true. ; use_lonq = .true.
      case ('u') ; use_lath = .true. ; use_lonq = .true.
      case ('v') ; use_latq = .true. ; use_lonh = .true.
      case ('T')  ; use_lath = .true. ; use_lonh = .true.
      case ('Bu') ; use_latq = .true. ; use_lonq = .true.
      case ('Cu') ; use_lath = .true. ; use_lonq = .true.
      case ('Cv') ; use_latq = .true. ; use_lonh = .true.
      case ('1') ! Do nothing.
      case default
        call mom_error(warning, "MOM_io create_file: "//trim(vars(k)%name)//&
                        " has unrecognized hor_grid "//trim(vars(k)%hor_grid))
    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
  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 mpp_get_domain_components(domain%mpp_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.")
 
! Specify all optional arguments to mpp_write_meta: name, units, longname, cartesian, calendar, sense,
! domain, data, min). Otherwise if optional arguments are added to mpp_write_meta the compiler may
! (and in case of GNU does) get confused and crash.
  if (use_lath) &
    call mpp_write_meta(unit, axis_lath, name="lath", units=y_axis_units, longname="Latitude", &
                   cartesian='Y', domain = y_domain, data=gridlatt(jsg:jeg))
 
  if (use_lonh) &
    call mpp_write_meta(unit, axis_lonh, name="lonh", units=x_axis_units, longname="Longitude", &
                   cartesian='X', domain = x_domain, data=gridlont(isg:ieg))
 
  if (use_latq) &
    call mpp_write_meta(unit, axis_latq, name="latq", units=y_axis_units, longname="Latitude", &
                   cartesian='Y', domain = y_domain, data=gridlatb(jsgb:jegb))
 
  if (use_lonq) &
    call mpp_write_meta(unit, axis_lonq, name="lonq", units=x_axis_units, longname="Longitude", &
                   cartesian='X', domain = x_domain, data=gridlonb(isgb:iegb))
 
  if (use_layer) &
    call mpp_write_meta(unit, axis_layer, name="Layer", units=trim(gv%zAxisUnits), &
          longname="Layer "//trim(gv%zAxisLongName), cartesian='Z', &
          sense=1, data=gv%sLayer(1:gv%ke))
 
  if (use_int) &
    call mpp_write_meta(unit, axis_int, name="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
 
    call mpp_write_meta(unit, axis_time, name="Time", units=time_units, longname="Time", cartesian='T')
  else
    call mpp_write_meta(unit, axis_time, name="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(period_val(num_periods))
    do k=1,num_periods ; period_val(k) = real(k) ; enddo
    call mpp_write_meta(unit, axis_periodic, name="Period", units="nondimensional", &
          longname="Periods for cyclical varaiables", cartesian= 't', data=period_val)
    deallocate(period_val)
  endif
 
  do k=1,novars
    numaxes = 0
    select case (vars(k)%hor_grid)
      case ('h')  ; numaxes = 2 ; axes(1) = axis_lonh ; axes(2) = axis_lath
      case ('q')  ; numaxes = 2 ; axes(1) = axis_lonq ; axes(2) = axis_latq
      case ('u')  ; numaxes = 2 ; axes(1) = axis_lonq ; axes(2) = axis_lath
      case ('v')  ; numaxes = 2 ; axes(1) = axis_lonh ; axes(2) = axis_latq
      case ('T')  ; numaxes = 2 ; axes(1) = axis_lonh ; axes(2) = axis_lath
      case ('Bu') ; numaxes = 2 ; axes(1) = axis_lonq ; axes(2) = axis_latq
      case ('Cu') ; numaxes = 2 ; axes(1) = axis_lonq ; axes(2) = axis_lath
      case ('Cv') ; numaxes = 2 ; axes(1) = axis_lonh ; axes(2) = axis_latq
      case ('1') ! Do nothing.
      case default
        call mom_error(warning, "MOM_io create_file: "//trim(vars(k)%name)//&
                        " has unrecognized hor_grid "//trim(vars(k)%hor_grid))
    end select
    select case (vars(k)%z_grid)
      case ('L') ; numaxes = numaxes+1 ; axes(numaxes) = axis_layer
      case ('i') ; numaxes = numaxes+1 ; axes(numaxes) = axis_int
      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') ; 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
       call mpp_write_meta(unit, fields(k), axes(1:numaxes), vars(k)%name, vars(k)%units, &
           vars(k)%longname, pack = pack, checksum=checksums(k,:))
    else
       call mpp_write_meta(unit, fields(k), axes(1:numaxes), vars(k)%name, vars(k)%units, &
           vars(k)%longname, pack = pack)
    endif
  enddo
 
  if (use_lath) call write_field(unit, axis_lath)
  if (use_latq) call write_field(unit, axis_latq)
  if (use_lonh) call write_field(unit, axis_lonh)
  if (use_lonq) call write_field(unit, axis_lonq)
  if (use_layer) call write_field(unit, axis_layer)
  if (use_int) call write_field(unit, axis_int)
  if (use_periodic) call write_field(unit, axis_periodic)
 
end subroutine create_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_file(unit, filename, vars, novars, fields, threading, timeunit, G, dG, GV)
  integer,               intent(out)   :: unit       !< unit id of an open file or -1 on a
                                                     !! nonwriting PE with single file output
  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(mom_domain_type), pointer :: domain => null()
  character(len=200) :: check_name, mesg
  integer :: length, ndim, nvar, natt, ntime, thread
  logical :: exists, one_file, domain_set
 
  thread = single_file
  if (PRESENT(threading)) thread = threading
 
  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_file(unit, filename, vars, novars, fields, threading, timeunit, &
                     g=g, dg=dg, gv=gv)
  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 open_file(unit, filename, mpp_append, mpp_netcdf, threading=thread)
    else
      call open_file(unit, filename, mpp_append, mpp_netcdf, domain=domain%mpp_domain)
    endif
    if (unit < 0) return
 
    call mpp_get_info(unit, ndim, nvar, natt, ntime)
 
    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_file(unit, filename, vars, novars, fields, threading, timeunit, g=g, gv=gv)
    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 mpp_get_fields(unit,fields(1:nvar))
 
    ! Check the field names...
!    do i=1,nvar
!      call mpp_get_field_atts(fields(i),name)
!      !if (trim(name) /= trim(vars%name) then
!      !write (mesg,'("Reopening file ",a," variable ",a," is called ",a,".")',&
!      !    filename,vars%name,name)
!      !call MOM_error(NOTE,"MOM_io: "//mesg)
!    enddo
  endif
 
end subroutine reopen_file
 
!> Read the data associated with a named axis in a file
subroutine read_axis_data(filename, axis_name, var)
  character(len=*),   intent(in)  :: filename  !< Name of the file to read
  character(len=*),   intent(in)  :: axis_name !< Name of the axis to read
  real, dimension(:), intent(out) :: var       !< The axis location data
 
  integer :: i,len,unit, ndim, nvar, natt, ntime
  logical :: axis_found
  type(axistype), allocatable :: axes(:)
  type(axistype) :: time_axis
  character(len=32) :: name, units
 
  call open_file(unit, trim(filename), action=mpp_rdonly, form=mpp_netcdf, &
                 threading=mpp_multi, fileset=single_file)
 
!Find the number of variables (nvar) in this file
  call mpp_get_info(unit, ndim, nvar, natt, ntime)
! -------------------------------------------------------------------
! Allocate space for the number of axes in the data file.
! -------------------------------------------------------------------
  allocate(axes(ndim))
  call mpp_get_axes(unit, axes, time_axis)
 
  axis_found = .false.
  do i = 1, ndim
    call mpp_get_atts(axes(i), name=name,len=len,units=units)
    if (name == axis_name) then
      axis_found = .true.
      call get_axis_data(axes(i),var)
      exit
    endif
  enddo
 
  if (.not.axis_found) call mom_error(fatal, "MOM_io read_axis_data: "//&
    "Unable to find axis "//trim(axis_name)//" in file "//trim(filename))
 
  deallocate(axes)
 
end subroutine read_axis_data
 
!> This function determines how many time levels a variable has.
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
 
  logical :: found
  character(len=200) :: msg
  character(len=nf90_max_name) :: name
  integer :: ncid, nvars, status, varid, ndims, n
  integer, allocatable :: varids(:)
  integer, dimension(nf90_max_var_dims) :: dimids
 
  n_time = -1
  found = .false.
 
  status = nf90_open(filename, nf90_nowrite, ncid)
  if (status /= nf90_noerr) then
    call mom_error(warning,"num_timelevels: "//&
        " Difficulties opening "//trim(filename)//" - "//&
        trim(nf90_strerror(status)))
    return
  endif
 
  status = nf90_inquire(ncid, nvariables=nvars)
  if (status /= nf90_noerr) then
    call mom_error(warning,"num_timelevels: "//&
        " Difficulties getting the number of variables in file "//&
        trim(filename)//" - "//trim(nf90_strerror(status)))
    return
  endif
 
  if (nvars < 1) then
    call mom_error(warning,"num_timelevels: "//&
        " There appear not to be any variables in "//trim(filename))
    return
  endif
 
 
  allocate(varids(nvars))
 
  status = nf90_inq_varids(ncid, nvars, varids)
  if (status /= nf90_noerr) then
    call mom_error(warning,"num_timelevels: "//&
        " Difficulties getting the variable IDs in file "//&
        trim(filename)//" - "//trim(nf90_strerror(status)))
    deallocate(varids) ; return
  endif
 
  do n = 1,nvars
    status = nf90_inquire_variable(ncid, varids(n), name=name)
    if (status /= nf90_noerr) then
      call mom_error(warning,"num_timelevels: "//&
          " Difficulties getting a variable name in file "//&
          trim(filename)//" - "//trim(nf90_strerror(status)))
    endif
 
    if (trim(lowercase(name)) == trim(lowercase(varname))) then
      if (found) then
        call mom_error(warning,"num_timelevels: "//&
          " Two variables match the case-insensitive name "//trim(varname)//&
          " in file "//trim(filename)//" - "//trim(nf90_strerror(status)))
      else
        varid = varids(n) ; found = .true.
      endif
    endif
  enddo
 
  deallocate(varids)
 
  if (.not.found) then
    call mom_error(warning,"num_timelevels: "//&
        " variable "//trim(varname)//" was not found in file "//&
        trim(filename))
    return
  endif
 
  status = nf90_inquire_variable(ncid, varid, ndims = ndims)
  if (status /= nf90_noerr) then
    call mom_error(warning,"num_timelevels: "//&
      trim(nf90_strerror(status))//" Getting number of dimensions of "//&
      trim(varname)//" in "//trim(filename))
    return
  endif
 
  if (present(min_dims)) then
    if (ndims < min_dims-1) then
      write(msg, '(I3)') min_dims
      call mom_error(warning, "num_timelevels: variable "//trim(varname)//&
        " in file "//trim(filename)//" has fewer than min_dims = "//trim(msg)//&
        " dimensions.")
    elseif (ndims == min_dims - 1) then
      n_time = 0 ; return
    endif
  endif
 
  status = nf90_inquire_variable(ncid, varid, dimids = dimids(1:ndims))
  if (status /= nf90_noerr) then
    call mom_error(warning,"num_timelevels: "//&
      trim(nf90_strerror(status))//" Getting last dimension ID for "//&
      trim(varname)//" in "//trim(filename))
    return
  endif
 
  status = nf90_inquire_dimension(ncid, dimids(ndims), len=n_time)
  if (status /= nf90_noerr) call mom_error(warning,"num_timelevels: "//&
      trim(nf90_strerror(status))//" Getting number of time levels of "//&
      trim(varname)//" in "//trim(filename))
 
  return
 
end function num_timelevels
 
 
!> 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) 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           !< variable horizonal staggering
  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
  character(len=*), optional, intent(in) :: caller             !< calling routine?
  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%z_grid = 'L' ; vd%t_grid = 's'
 
  vd%cmor_field_name  =  ""
  vd%cmor_units       =  ""
  vd%cmor_longname    =  ""
  vd%conversion       =  1.0
 
  call modify_vardesc(vd, units=units, longname=longname, hor_grid=hor_grid, &
                      z_grid=z_grid, t_grid=t_grid,                          &
                      cmor_field_name=cmor_field_name,cmor_units=cmor_units, &
                      cmor_longname=cmor_longname, conversion=conversion, caller=cllr)
 
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)
  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        !< horizonal 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      !< for unit conversions, such as needed
                                                               !! to convert from intensive to extensive
  character(len=*), optional, intent(in)    :: caller          !< calling routine?
 
  character(len=120) :: cllr
  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(hor_grid))  call safe_string_copy(hor_grid, vd%hor_grid, &
                               "vd%hor_grid of "//trim(vd%name), cllr)
  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)
 
end subroutine modify_vardesc
 
!> 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)
  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           !< horiz staggering of variable
  character(len=*), optional, intent(out) :: z_grid             !< vert 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
  character(len=*), optional, intent(in)  :: caller             !< calling routine?
 
 
  character(len=120) :: cllr
  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)
 
end subroutine query_vardesc
 
 
!> 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, ie
 
  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, '(I10)') ens_no ; ens_num_char = adjustl(ens_num_char)
  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
 
 
!> Returns true if the named file or its domain-decomposed variant exists.
function mom_file_exists(filename, MOM_Domain)
  character(len=*),       intent(in) :: filename   !< The name of the file being inquired about
  type(mom_domain_type),  intent(in) :: mom_domain !< The MOM_Domain that describes the decomposition
 
! This function uses the fms_io function file_exist to determine whether
! a named file (or its decomposed variant) exists.
 
  logical :: mom_file_exists
 
  mom_file_exists = file_exist(filename, mom_domain%mpp_domain)
 
end function mom_file_exists
 
!> Returns true if the named file or its domain-decomposed variant exists.
function fms_file_exists(filename, domain, no_domain)
  character(len=*), intent(in)         :: filename  !< The name of the file being inquired about
  type(domain2d), optional, intent(in) :: domain    !< The mpp domain2d that describes the decomposition
  logical,        optional, intent(in) :: no_domain !< This file does not use domain decomposition
! This function uses the fms_io function file_exist to determine whether
! a named file (or its decomposed variant) exists.
 
  logical :: fms_file_exists
 
  fms_file_exists = file_exist(filename, domain, no_domain)
 
end function fms_file_exists
 
!> This function uses the fms_io function read_data to read 1-D
!! data field named "fieldname" from file "filename".
subroutine mom_read_data_1d(filename, fieldname, data, timelevel, scale)
  character(len=*),       intent(in)    :: filename  !< The name of the file to read
  character(len=*),       intent(in)    :: fieldname !< The variable name of the data in the file
  real, dimension(:),     intent(inout) :: data      !< The 1-dimensional array into which the data
  integer,      optional, intent(in)    :: timelevel !< The time level in the file to read
  real,         optional, intent(in)    :: scale     !< A scaling factor that the field is multiplied
                                                     !! by before they are returned.
 
  call read_data(filename, fieldname, data, timelevel=timelevel, no_domain=.true.)
 
  if (present(scale)) then ; if (scale /= 1.0) then
    data(:) = scale*data(:)
  endif ; endif
 
end subroutine mom_read_data_1d
 
!> This function uses the fms_io function read_data to read a distributed
!! 2-D data field named "fieldname" from file "filename".  Valid values for
!! "position" include CORNER, CENTER, EAST_FACE and NORTH_FACE.
subroutine mom_read_data_2d(filename, fieldname, data, MOM_Domain, &
                            timelevel, position, scale)
  character(len=*),       intent(in)    :: filename  !< The name of the file to read
  character(len=*),       intent(in)    :: fieldname !< The variable name of the data in the file
  real, dimension(:,:),   intent(inout) :: data      !< The 2-dimensional array into which the data
                                                     !! should be read
  type(mom_domain_type),  intent(in)    :: MOM_Domain !< The MOM_Domain that describes the decomposition
  integer,      optional, intent(in)    :: timelevel !< The time level in the file to read
  integer,      optional, intent(in)    :: position  !< A flag indicating where this data is located
  real,         optional, intent(in)    :: scale     !< A scaling factor that the field is multiplied
                                                     !! by before they are returned.
 
  integer :: is, ie, js, je
 
  call read_data(filename, fieldname, data, mom_domain%mpp_domain, &
                 timelevel=timelevel, position=position)
 
  if (present(scale)) then ; if (scale /= 1.0) then
    call get_simple_array_i_ind(mom_domain, size(data,1), is, ie)
    call get_simple_array_j_ind(mom_domain, size(data,2), js, je)
    data(is:ie,js:je) = scale*data(is:ie,js:je)
  endif ; endif
 
end subroutine mom_read_data_2d
 
!> This function uses the fms_io function read_data to read a distributed
!! 3-D data field named "fieldname" from file "filename".  Valid values for
!! "position" include CORNER, CENTER, EAST_FACE and NORTH_FACE.
subroutine mom_read_data_3d(filename, fieldname, data, MOM_Domain, &
                            timelevel, position, scale)
  character(len=*),       intent(in)    :: filename  !< The name of the file to read
  character(len=*),       intent(in)    :: fieldname !< The variable name of the data in the file
  real, dimension(:,:,:), intent(inout) :: data      !< The 3-dimensional array into which the data
                                                     !! should be read
  type(mom_domain_type),  intent(in)    :: MOM_Domain !< The MOM_Domain that describes the decomposition
  integer,      optional, intent(in)    :: timelevel !< The time level in the file to read
  integer,      optional, intent(in)    :: position  !< A flag indicating where this data is located
  real,         optional, intent(in)    :: scale     !< A scaling factor that the field is multiplied
                                                     !! by before they are returned.
 
  integer :: is, ie, js, je
 
  call read_data(filename, fieldname, data, mom_domain%mpp_domain, &
                 timelevel=timelevel, position=position)
 
  if (present(scale)) then ; if (scale /= 1.0) then
    call get_simple_array_i_ind(mom_domain, size(data,1), is, ie)
    call get_simple_array_j_ind(mom_domain, size(data,2), js, je)
    data(is:ie,js:je,:) = scale*data(is:ie,js:je,:)
  endif ; endif
 
end subroutine mom_read_data_3d
 
!> This function uses the fms_io function read_data to read a distributed
!! 4-D data field named "fieldname" from file "filename".  Valid values for
!! "position" include CORNER, CENTER, EAST_FACE and NORTH_FACE.
subroutine mom_read_data_4d(filename, fieldname, data, MOM_Domain, &
                            timelevel, position, scale)
  character(len=*),       intent(in)    :: filename  !< The name of the file to read
  character(len=*),       intent(in)    :: fieldname !< The variable name of the data in the file
  real, dimension(:,:,:,:), intent(inout) :: data    !< The 4-dimensional array into which the data
                                                     !! should be read
  type(mom_domain_type),  intent(in)    :: MOM_Domain !< The MOM_Domain that describes the decomposition
  integer,      optional, intent(in)    :: timelevel !< The time level in the file to read
  integer,      optional, intent(in)    :: position  !< A flag indicating where this data is located
  real,         optional, intent(in)    :: scale     !< A scaling factor that the field is multiplied
                                                     !! by before they are returned.
 
  integer :: is, ie, js, je
 
  call read_data(filename, fieldname, data, mom_domain%mpp_domain, &
                 timelevel=timelevel, position=position)
 
  if (present(scale)) then ; if (scale /= 1.0) then
    call get_simple_array_i_ind(mom_domain, size(data,1), is, ie)
    call get_simple_array_j_ind(mom_domain, size(data,2), js, je)
    data(is:ie,js:je,:,:) = scale*data(is:ie,js:je,:,:)
  endif ; endif
 
end subroutine mom_read_data_4d
 
 
!> This function uses the fms_io function read_data to read a pair of distributed
!! 2-D data fields with names given by "[uv]_fieldname" from file "filename".  Valid values for
!! "stagger" include CGRID_NE, BGRID_NE, and AGRID.
subroutine mom_read_vector_2d(filename, u_fieldname, v_fieldname, u_data, v_data, MOM_Domain, &
                              timelevel, stagger, scalar_pair, scale)
  character(len=*),       intent(in)    :: filename  !< The name of the file to read
  character(len=*),       intent(in)    :: u_fieldname !< The variable name of the u data in the file
  character(len=*),       intent(in)    :: v_fieldname !< The variable name of the v data in the file
  real, dimension(:,:),   intent(inout) :: u_data    !< The 2 dimensional array into which the
                                                     !! u-component of the data should be read
  real, dimension(:,:),   intent(inout) :: v_data    !< The 2 dimensional array into which the
                                                     !! v-component of the data should be read
  type(mom_domain_type),  intent(in)    :: MOM_Domain !< The MOM_Domain that describes the decomposition
  integer,      optional, intent(in)    :: timelevel !< The time level in the file to read
  integer,      optional, intent(in)    :: stagger   !< A flag indicating where this vector is discretized
  logical,      optional, intent(in)    :: scalar_pair !< If true, a pair of scalars are to be read.cretized
  real,         optional, intent(in)    :: scale     !< A scaling factor that the fields are multiplied
                                                     !! by before they are returned.
  integer :: is, ie, js, je
  integer :: u_pos, v_pos
 
  u_pos = east_face ; v_pos = north_face
  if (present(stagger)) then
    if (stagger == cgrid_ne) then ; u_pos = east_face ; v_pos = north_face
    elseif (stagger == bgrid_ne) then ; u_pos = corner ; v_pos = corner
    elseif (stagger == agrid) then ; u_pos = center ; v_pos = center ; endif
  endif
 
  call read_data(filename, u_fieldname, u_data, mom_domain%mpp_domain, &
                 timelevel=timelevel, position=u_pos)
  call read_data(filename, v_fieldname, v_data, mom_domain%mpp_domain, &
                 timelevel=timelevel, position=v_pos)
 
  if (present(scale)) then ; if (scale /= 1.0) then
    call get_simple_array_i_ind(mom_domain, size(u_data,1), is, ie)
    call get_simple_array_j_ind(mom_domain, size(u_data,2), js, je)
    u_data(is:ie,js:je) = scale*u_data(is:ie,js:je)
    call get_simple_array_i_ind(mom_domain, size(v_data,1), is, ie)
    call get_simple_array_j_ind(mom_domain, size(v_data,2), js, je)
    v_data(is:ie,js:je) = scale*v_data(is:ie,js:je)
  endif ; endif
 
end subroutine mom_read_vector_2d
 
 
!> This function uses the fms_io function read_data to read a pair of distributed
!! 3-D data fields with names given by "[uv]_fieldname" from file "filename".  Valid values for
!! "stagger" include CGRID_NE, BGRID_NE, and AGRID.
subroutine mom_read_vector_3d(filename, u_fieldname, v_fieldname, u_data, v_data, MOM_Domain, &
                              timelevel, stagger, scalar_pair, scale)
  character(len=*),       intent(in)    :: filename  !< The name of the file to read
  character(len=*),       intent(in)    :: u_fieldname !< The variable name of the u data in the file
  character(len=*),       intent(in)    :: v_fieldname !< The variable name of the v data in the file
  real, dimension(:,:,:), intent(inout) :: u_data    !< The 3 dimensional array into which the
                                                     !! u-component of the data should be read
  real, dimension(:,:,:), intent(inout) :: v_data    !< The 3 dimensional array into which the
                                                     !! v-component of the data should be read
  type(mom_domain_type),  intent(in)    :: MOM_Domain !< The MOM_Domain that describes the decomposition
  integer,      optional, intent(in)    :: timelevel !< The time level in the file to read
  integer,      optional, intent(in)    :: stagger   !< A flag indicating where this vector is discretized
  logical,      optional, intent(in)    :: scalar_pair !< If true, a pair of scalars are to be read.cretized
  real,         optional, intent(in)    :: scale     !< A scaling factor that the fields are multiplied
                                                     !! by before they are returned.
 
  integer :: is, ie, js, je
  integer :: u_pos, v_pos
 
  u_pos = east_face ; v_pos = north_face
  if (present(stagger)) then
    if (stagger == cgrid_ne) then ; u_pos = east_face ; v_pos = north_face
    elseif (stagger == bgrid_ne) then ; u_pos = corner ; v_pos = corner
    elseif (stagger == agrid) then ; u_pos = center ; v_pos = center ; endif
  endif
 
  call read_data(filename, u_fieldname, u_data, mom_domain%mpp_domain, &
                 timelevel=timelevel, position=u_pos)
  call read_data(filename, v_fieldname, v_data, mom_domain%mpp_domain, &
                 timelevel=timelevel, position=v_pos)
 
  if (present(scale)) then ; if (scale /= 1.0) then
    call get_simple_array_i_ind(mom_domain, size(u_data,1), is, ie)
    call get_simple_array_j_ind(mom_domain, size(u_data,2), js, je)
    u_data(is:ie,js:je,:) = scale*u_data(is:ie,js:je,:)
    call get_simple_array_i_ind(mom_domain, size(v_data,1), is, ie)
    call get_simple_array_j_ind(mom_domain, size(v_data,2), js, je)
    v_data(is:ie,js:je,:) = scale*v_data(is:ie,js:je,:)
  endif ; endif
 
end subroutine mom_read_vector_3d
 
 
!> 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
 
 
!> \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