Direct netCDF data reads

This patch introduces `read_netCDF_data`, a new method for reading
netCDF datasets using the native netCDF I/O interface.  It is designed
to resemble the existing `MOM_read_data` function.

Motivation
----------

Legacy input files may contain content which is not supported by the
newest framework I/O (FMS).  In order to retain support for these input
files, particularly over a wider range of compilers, this patch provides
an alternative method for reading these files.

Interface
---------

As with `MOM_read_data`, the function is provided with a netCDF filepath
and a variable name, and returns the values to a provided variable.

The `global_file` and `file_may_be_4d` flags have been dropped, since
they are related to specific FMS2 compatibility issues.  (Global vs
domain-decomposed reads is controlled by the presence of a `MOM_domain`)

Limited domain-decomposed I/O is supported, providing parallel I/O over
a single file, to the extent supported by the filesystem.
Parallelization over multiple files, as in FMS I/O, is not supported.
Each FMS PE (MPI rank) reads its own segment, as defined by its
MOM_domain.

Output can be saved to either compute or data domains; as in FMS, the
appropriate placement is inferred from the size of the output array.

Support is currently limited to time-independent 2D arrays with
center-cell indexing.  That is, the `position` and `timelevel` arguments
are not yet supported.  The subroutines raise an error if these are
provided, as an indication that they may support them in the future.

Implementation
--------------

Internally, the function opens a `MOM_netcdf_file`, generates its
field/axis manifest, and reads the field contents.  As with
`MOM_read_data`, an internal rotation may be applied.  The file is
closed upon completion.

(This behavior is designed to emulate the existing `MOM_read_data`; in a
future implementation, we may want to use a persistent file handle which
reduces the number of I/O operations.)

Opening a `MOM_netcdf_file` now supports a `MOM_domain` argument, which
is used to determine the index bounds of its local segment of the global
domain.  This is used to compute appropriate bounds for the native
netCDF IO calls.

As part of these changes, the `get_file_fields` function has been
separated into itself and a new function, `update_file_contents`, which
populates the internal axis and field metadata list.

Usage
-----

The following fields have been moved to the native netCDF IO:

* `tideamp` (tidal mixing, FMS surface forcing)
* `gustiness` (solo and FMS surface forcing)
* `h2` (roughness in tidal mixing)

This only comprises the fields which must be handled natively in order
for the GFDL regression suite to pass with the PGI compiler; more files
could be moved to native I/O in the future.

Bugfixes
--------

Some bugfixes to the netCDF I/O are also included:

* `filename` attribute is now only written in an writeable state

* Previously, `get_file_fields` (and now `update_file_contents`) assumed
  that every axis had an equivalent variable, which could lead to
  potential errors if an axis had no equvalent field, such as index
  bounds.

  We now count the number of variables with matching dimension names,
  tagged as axes, rather than assuming that every axis has a variable,
  and exclude them from the field list.

* Not a bugfix, but `hor_index_init` was modified so that `param_file`
  is now an optional input.  This function is used in `MOM_netcdf_file`,
  where `param_file` is not available.  The argument is only used to
  call `log_param`.

Previous usage of these functions was restricted to writing output with
well-defined content, so were unaffected by these issues.

config_src/drivers/FMS_cap/MOM_surface_forcing_gfdl.F90

@@ -28,6 +28,7 @@ module MOM_surface_forcing_gfdl
 use MOM_interpolate,      only : init_external_field, time_interp_external
 use MOM_interpolate,      only : time_interp_external_init
 use MOM_io,               only : slasher, write_version_number, MOM_read_data
+use MOM_io,               only : read_netCDF_data
 use MOM_io,               only : stdout_if_root
 use MOM_restart,          only : register_restart_field, restart_init, MOM_restart_CS
 use MOM_restart,          only : restart_init_end, save_restart, restore_state
@@ -1507,7 +1508,10 @@ subroutine surface_forcing_init(Time, G, US, param_file, diag, CS, wind_stagger)
 
   if (CS%read_TIDEAMP) then
     TideAmp_file = trim(CS%inputdir) // trim(TideAmp_file)
-    call MOM_read_data(TideAmp_file,'tideamp',CS%TKE_tidal,G%domain,timelevel=1, scale=US%m_to_Z*US%T_to_s)
+    ! NOTE: There are certain cases where FMS is unable to read this file, so
+    ! we use read_netCDF_data in place of MOM_read_data.
+    call read_netCDF_data(TideAmp_file, 'tideamp', CS%TKE_tidal, G%Domain, &
+        rescale=US%m_to_Z*US%T_to_s)
     do j=jsd, jed; do i=isd, ied
       utide = CS%TKE_tidal(i,j)
       CS%TKE_tidal(i,j) = G%mask2dT(i,j)*CS%Rho0*CS%cd_tides*(utide*utide*utide)
@@ -1537,8 +1541,10 @@ subroutine surface_forcing_init(Time, G, US, param_file, diag, CS, wind_stagger)
 
     call safe_alloc_ptr(CS%gust,isd,ied,jsd,jed)
     gust_file = trim(CS%inputdir) // trim(gust_file)
-    call MOM_read_data(gust_file, 'gustiness', CS%gust, G%domain, timelevel=1, &
-               scale=US%kg_m3_to_R*US%m_s_to_L_T**2*US%L_to_Z) ! units in file should be Pa
+    ! NOTE: There are certain cases where FMS is unable to read this file, so
+    ! we use read_netCDF_data in place of MOM_read_data.
+    call read_netCDF_data(gust_file, 'gustiness', CS%gust, G%Domain, &
+        rescale=US%kg_m3_to_R*US%m_s_to_L_T**2*US%L_to_Z) ! units in file should be Pa
   endif
   call get_param(param_file, mdl, "DEFAULT_ANSWER_DATE", default_answer_date, &
                  "This sets the default value for the various _ANSWER_DATE parameters.", &

config_src/drivers/solo_driver/MOM_surface_forcing.F90

@@ -30,6 +30,7 @@ module MOM_surface_forcing
 use MOM_grid,                only : ocean_grid_type
 use MOM_get_input,           only : Get_MOM_Input, directories
 use MOM_io,                  only : file_exists, MOM_read_data, MOM_read_vector, slasher
+use MOM_io,                  only : read_netCDF_data
 use MOM_io,                  only : EAST_FACE, NORTH_FACE, num_timelevels
 use MOM_restart,             only : register_restart_field, restart_init, MOM_restart_CS
 use MOM_restart,             only : restart_init_end, save_restart, restore_state
@@ -1866,8 +1867,10 @@ subroutine surface_forcing_init(Time, G, US, param_file, diag, CS, tracer_flow_C
                  "variable gustiness.", fail_if_missing=.true.)
     call safe_alloc_ptr(CS%gust,G%isd,G%ied,G%jsd,G%jed)
     filename = trim(CS%inputdir) // trim(gust_file)
-    call MOM_read_data(filename,'gustiness',CS%gust,G%domain, timelevel=1, &
-                   scale=Pa_to_RLZ_T2) ! units in file should be Pa
+    ! NOTE: There are certain cases where FMS is unable to read this file, so
+    ! we use read_netCDF_data in place of MOM_read_data.
+    call read_netCDF_data(filename, 'gustiness', CS%gust, G%Domain, &
+        rescale=Pa_to_RLZ_T2) ! units in file should be Pa
   endif
 
 !  All parameter settings are now known.

src/framework/MOM_hor_index.F90

@@ -63,7 +63,7 @@ module MOM_hor_index
 subroutine hor_index_init(Domain, HI, param_file, local_indexing, index_offset)
   type(MOM_domain_type),  intent(in)    :: Domain     !< The MOM domain from which to extract information.
   type(hor_index_type),   intent(inout) :: HI         !< A horizontal index type to populate with data
-  type(param_file_type),  intent(in)    :: param_file !< Parameter file handle
+  type(param_file_type), optional, intent(in) :: param_file !< Parameter file handle
   logical, optional,      intent(in)    :: local_indexing !< If true, all tracer data domains start at 1
   integer, optional,      intent(in)    :: index_offset   !< A fixed additional offset to all indices
 
@@ -80,8 +80,9 @@ subroutine hor_index_init(Domain, HI, param_file, local_indexing, index_offset)
   call get_global_shape(Domain, HI%niglobal, HI%njglobal)
 
   ! Read all relevant parameters and write them to the model log.
-  call log_version(param_file, "MOM_hor_index", version, &
-                   "Sets the horizontal array index types.", all_default=.true.)
+  if (present(param_file)) &
+    call log_version(param_file, "MOM_hor_index", version, &
+                     "Sets the horizontal array index types.", all_default=.true.)
 
   HI%IscB = HI%isc ; HI%JscB = HI%jsc
   HI%IsdB = HI%isd ; HI%JsdB = HI%jsd

src/framework/MOM_io.F90

@@ -61,6 +61,7 @@ module MOM_io
 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
@@ -108,6 +109,15 @@ module MOM_io
   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
@@ -2033,6 +2043,59 @@ subroutine MOM_read_data_2d(filename, fieldname, data, MOM_Domain, &
 end subroutine MOM_read_data_2d
 
 
+!> Read a 2d array from file using native netCDF I/O.
+subroutine read_netCDF_data_2d(filename, fieldname, values, MOM_Domain, &
+                            timelevel, position, rescale)
+  character(len=*), intent(in) :: filename
+    !< Input filename
+  character(len=*), intent(in)  :: fieldname
+    !< Field variable name
+  real, intent(out) :: values(:,:)
+    !< Field value
+  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
+
+  integer :: turns
+    ! 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()
+
+  turns = MOM_domain%turns
+  if (turns == 0) then
+    call handle%read(fieldname, values, rescale=rescale)
+  else
+    call allocate_rotated_array(values, [1,1], -turns, values_in)
+    call handle%read(fieldname, values_in, rescale=rescale)
+    call rotate_array(values_in, turns, 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)

src/framework/MOM_io_file.F90

@@ -15,6 +15,9 @@ module MOM_io_file
 use MOM_io_infra, only : get_field_atts
 use MOM_io_infra, only : read_field_chksum
 
+use MOM_hor_index, only : hor_index_type
+use MOM_hor_index, only : hor_index_init
+
 use MOM_netcdf, only : netcdf_file_type
 use MOM_netcdf, only : netcdf_axis
 use MOM_netcdf, only : netcdf_field
@@ -28,6 +31,7 @@ module MOM_io_file
 use MOM_netcdf, only : write_netcdf_attribute
 use MOM_netcdf, only : get_netcdf_size
 use MOM_netcdf, only : get_netcdf_fields
+use MOM_netcdf, only : read_netcdf_field
 
 use MOM_error_handler, only : MOM_error, FATAL
 use MOM_error_handler, only : is_root_PE
@@ -43,9 +47,9 @@ module MOM_io_file
 
 ! Internal types
 
-! NOTE: MOM_axis and MOM_field do not represent the actual axes and
-! fields stored in the file.  They are only very thin wrappers to the keys (as
-! strings) used to reference the associated object inside the MOM_file.
+! NOTE: MOM_axis and MOM_field do not contain the actual axes and fields stored
+! in the file.  They are very thin wrappers to the keys (as strings) used to
+! reference the associated object inside of the MOM_file.
 
 !> Handle for axis in MOM file
 type :: MOM_axis
@@ -316,6 +320,10 @@ module MOM_io_file
   type(field_list_nc) :: fields
   !> True if the file has been opened
   logical :: is_open = .false.
+  !> True if I/O content is domain-decomposed
+  logical :: domain_decomposed = .false.
+  !> True if I/O content is domain-decomposed
+  type(hor_index_type) :: HI
 
   contains
 
@@ -356,6 +364,12 @@ module MOM_io_file
   procedure :: get_field_atts => get_field_atts_nc
   !> Get checksum from a netCDF field
   procedure :: read_field_chksum => read_field_chksum_nc
+
+  ! NOTE: These are currently exclusive to netCDF I/O but could be generalized
+  !> Read the values of a netCDF field
+  procedure :: read => get_field_nc
+  !> Update the axes and fields descriptors of a MOM netCDF file
+  procedure :: update => update_file_contents_nc
 end type MOM_netcdf_file
 
 
@@ -1281,11 +1295,16 @@ subroutine open_file_nc(handle, filename, action, MOM_domain, threading, fileset
   integer, intent(in), optional :: threading
   integer, intent(in), optional :: fileset
 
-  if (.not. is_root_PE()) return
+  if (.not. present(MOM_domain) .and. .not. is_root_PE()) return
 
   call open_netcdf_file(handle%handle_nc, filename, action)
-
   handle%is_open = .true.
+
+  if (present(MOM_domain)) then
+    handle%domain_decomposed = .true.
+    call hor_index_init(MOM_domain, handle%HI)
+  endif
+
   call handle%axes%init()
   call handle%fields%init()
 end subroutine open_file_nc
@@ -1295,7 +1314,7 @@ end subroutine open_file_nc
 subroutine close_file_nc(handle)
   class(MOM_netcdf_file), intent(inout) :: handle
 
-  if (.not. is_root_PE()) return
+  if (.not. handle%domain_decomposed .and. .not. is_root_PE()) return
 
   handle%is_open = .false.
   call close_netcdf_file(handle%handle_nc)
@@ -1575,31 +1594,56 @@ subroutine get_file_info_nc(handle, ndim, nvar, ntime)
 end subroutine get_file_info_nc
 
 
-!> Return the field metadata associated with a MOM netCDF file
-subroutine get_file_fields_nc(handle, fields)
+!> Update the axes and fields descriptors of a MOM netCDF file
+subroutine update_file_contents_nc(handle)
   class(MOM_netcdf_file), intent(inout) :: handle
     !< Handle for a file that is open for I/O
-  type(MOM_field), intent(inout) :: fields(:)
-    !< Field-type descriptions of all of the variables in a file.
 
   type(netcdf_axis), allocatable :: axes_nc(:)
+    ! netCDF axis descriptors
   type(netcdf_field), allocatable :: fields_nc(:)
+    ! netCDF field descriptors
   integer :: i
+    ! Index counter
 
-  if (.not. is_root_PE()) return
+  if (.not. handle%domain_decomposed .and. .not. is_root_PE()) return
 
   call get_netcdf_fields(handle%handle_nc, axes_nc, fields_nc)
-  if (size(fields) /= size(fields_nc)) &
-    call MOM_error(FATAL, 'Number of fields in file does not match field(:).')
 
   do i = 1, size(axes_nc)
     call handle%axes%append(axes_nc(i), axes_nc(i)%label)
   enddo
 
-  do i = 1, size(fields)
-    fields(i)%label = trim(fields_nc(i)%label)
+  do i = 1, size(fields_nc)
     call handle%fields%append(fields_nc(i), fields_nc(i)%label)
   enddo
+end subroutine update_file_contents_nc
+
+
+!> Return the field descriptors of a MOM netCDF file
+subroutine get_file_fields_nc(handle, fields)
+  class(MOM_netcdf_file), intent(inout) :: handle
+    !< Handle for a file that is open for I/O
+  type(MOM_field), intent(inout) :: fields(:)
+    !< Field-type descriptions of all of the variables in a file.
+
+  type(field_node_nc), pointer :: node => null()
+    ! Current field list node
+  integer :: n
+    ! Field counter
+
+  if (.not. is_root_PE()) return
+
+  ! Generate the manifest of axes and fields
+  call handle%update()
+
+  n = 0
+  node => handle%fields%head
+  do while (associated(node%next))
+    n = n + 1
+    fields(n)%label = trim(node%label)
+    node => node%next
+  enddo
 end subroutine get_file_fields_nc
 
 
@@ -1637,6 +1681,87 @@ subroutine read_field_chksum_nc(handle, field, chksum, valid_chksum)
 end subroutine read_field_chksum_nc
 
 
+!> Read the values of a netCDF field
+subroutine get_field_nc(handle, label, values, rescale)
+  class(MOM_netcdf_file), intent(in) :: handle
+    ! Handle of netCDF file to be read
+  character(len=*), intent(in) :: label
+    ! Field variable name
+  real, intent(out) :: values(:,:)
+    ! Field values read from file
+  real, optional, intent(in) :: rescale
+
+  logical :: data_domain
+    ! True if values matches the data domain size
+  logical :: compute_domain
+    ! True if values matches the compute domain size
+  type(netcdf_field) :: field_nc
+    ! netCDF field associated with label
+  integer :: isc, iec, jsc, jec
+    ! Index bounds of compute domain
+  integer :: isd, ied, jsd, jed
+    ! Index bounds of data domain
+  integer :: bounds(2,2)
+    ! Index bounds of domain
+  real, allocatable :: values_nc(:,:)
+    ! Local copy of field, used for data-decomposed I/O
+
+  isc = handle%HI%isc
+  iec = handle%HI%iec
+  jsc = handle%HI%jsc
+  jec = handle%HI%jec
+
+  isd = handle%HI%isd
+  ied = handle%HI%ied
+  jsd = handle%HI%jsd
+  jed = handle%HI%jed
+
+  data_domain = all(shape(values) == [ied-isd+1, jed-jsd+1])
+  compute_domain = all(shape(values) == [iec-isc+1, jec-jsc+1])
+
+  ! NOTE: Data on face and vertex points is not yet supported.  This is a
+  ! temporary check to detect such cases, but may be removed in the future.
+  if (.not. (compute_domain .or. data_domain)) &
+    call MOM_error(FATAL, 'get_field_nc: Only compute and data domains ' // &
+        'are currently supported.')
+
+  field_nc = handle%fields%get(label)
+
+  if (data_domain) then
+    allocate(values_nc(isc:iec,jsc:jec))
+    values(:,:) = 0.
+  endif
+
+  if (handle%domain_decomposed) then
+    bounds(1,:) = [isc, jsc] + [handle%HI%idg_offset, handle%HI%jdg_offset]
+    bounds(2,:) = [iec, jec] + [handle%HI%idg_offset, handle%HI%jdg_offset]
+    if (data_domain) then
+      call read_netcdf_field(handle%handle_nc, field_nc, values_nc, bounds=bounds)
+    else
+      call read_netcdf_field(handle%handle_nc, field_nc, values, bounds=bounds)
+    endif
+  else
+    if (data_domain) then
+      call read_netcdf_field(handle%handle_nc, field_nc, values_nc)
+    else
+      call read_netcdf_field(handle%handle_nc, field_nc, values)
+    endif
+  endif
+
+  if (data_domain) &
+    values(isc:iec,jsc:jec) = values_nc(:,:)
+
+  ! NOTE: It is more efficient to do the rescale in-place while copying
+  ! values_nc(:,:) to values(:,:).  But since rescale is only present for
+  ! debugging, we can probably disregard this impact on performance.
+  if (present(rescale)) then
+    if (rescale /= 1.0) then
+      values(isc:iec,jsc:jec) = rescale * values(isc:iec,jsc:jec)
+    endif
+  endif
+end subroutine get_field_nc
+
+
 !> \namespace MOM_IO_file
 !!
 !! This file defines the MOM_file classes used to inferface with the internal