FMS2: get_field_sizes with FMS1 emulation

This patch modifies the `get_field_sizes` function to emulate the
FMS1 behavior when the requested `sizes(:)` is larger than the field's
number of dimensions.

In FMS1, `field_sizes` would always expect a `sizes` array of at least
size 4, and in the format of (nx, ny, nz, nt), and would lean on
netCDF attributes like `cartesian_axes` when assigning the values.

FMS2 lacks similar functionality, due its more general approach.
Previously, this was emulated in a way favorable to certain situations,
but raised issues in others.

This patch attempts to resolve the issue by using the
`categorize_axes()` function's ability to identify axes, either by
attributes, presumed names, or the `unlimited` property.

config_src/infra/FMS2/MOM_io_infra.F90

@@ -609,11 +609,14 @@ subroutine get_field_size(filename, fieldname, sizes, field_found, no_domain)
                                                     !! names with an appended tile number
   ! Local variables
   type(FmsNetcdfFile_t) :: fileobj_read ! A handle to a non-domain-decomposed file for obtaining information
-                                              ! about the exiting time axis entries in append mode.
+                                        ! about the exiting time axis entries in append mode.
   logical :: success         ! If true, the file was opened successfully
   logical :: field_exists    ! True if filename exists and field_name is in filename
   integer :: i, ndims
-
+  character(len=512), allocatable :: dimnames(:)  ! Field dimension names
+  logical, allocatable :: is_x(:), is_y(:), is_t(:)     ! True if index matches axis type
+  integer :: size_indices(4)        ! Mapping of size index to FMS1 convention
+  integer :: idx, swap
 
   if (FMS2_reads) then
     field_exists = .false.
@@ -626,12 +629,53 @@ subroutine get_field_size(filename, fieldname, sizes, field_found, no_domain)
           if (ndims > size(sizes)) call MOM_error(FATAL, &
             "get_field_size called with too few sizes for "//trim(fieldname)//" in "//trim(filename))
           call get_variable_size(fileobj_read, fieldname, sizes(1:ndims))
+
           do i=ndims+1,size(sizes) ; sizes(i) = 0 ; enddo
-          ! This preserves previous behavior when reading time-varying data without
-          ! a vertical extent.
-          if (size(sizes)==ndims+1) then
-            sizes(ndims+1)=sizes(ndims)
-            sizes(ndims)=1
+
+          ! If sizes exceeds ndims, then we fallback to the FMS1 convention
+          ! where sizes has at least 4 dimension, and try to position values.
+          if (size(sizes) > ndims)  then
+            ! Assume FMS1 positioning rules: (nx, ny, nz, nt, ...)
+            if (size(sizes) < 4) &
+              call MOM_error(FATAL, "If sizes(:) exceeds field dimensions, "&
+                  &"then its length must be at least 4.")
+
+            ! Fall back to the FMS1 default values of 1 (from mpp field%size)
+            sizes(ndims+1:) = 1
+
+            ! Gather the field dimension names
+            allocate(dimnames(ndims))
+            dimnames(:) = ""
+            call get_variable_dimension_names(fileObj_read, trim(fieldname), &
+                                              dimnames)
+
+            ! Test the dimensions against standard (x,y,t) names and attributes
+            allocate(is_x(ndims), is_y(ndims), is_t(ndims))
+            is_x(:) = .false.
+            is_y(:) = .false.
+            is_t(:) = .false.
+            call categorize_axes(fileObj_read, filename, ndims, dimnames, &
+                                 is_x, is_y, is_t)
+
+            ! Currently no z-test is supported, so disable assignment with 0
+            size_indices = [ &
+                findloc(is_x, .true.), &
+                findloc(is_y, .true.), &
+                0, &
+                findloc(is_t, .true.) &
+            ]
+
+            do i = 1, size(size_indices)
+              idx = size_indices(i)
+              if (idx > 0) then
+                swap = sizes(i)
+                sizes(i) = sizes(idx)
+                sizes(idx) = swap
+              endif
+            enddo
+
+            deallocate(is_x, is_y, is_t)
+            deallocate(dimnames)
           endif
         endif
       endif
@@ -1483,7 +1527,7 @@ end subroutine MOM_register_variable_axes
 !> Determine whether a variable's axes are associated with x-, y- or time-dimensions.  Other
 !! unlimited dimensions are also labeled as time axes for these purposes.
 subroutine categorize_axes(fileObj, filename, ndims, dim_names, is_x, is_y, is_t)
-  type(FmsNetcdfDomainFile_t), intent(in)  :: fileObj  !< Handle to an open FMS2 netCDF file object
+  class(FmsNetcdfFile_t), intent(in) :: fileObj    !< Handle to an open FMS2 netCDF file object
   character(len=*),            intent(in)  :: filename !< The name of the file to read
   integer,                     intent(in)  :: ndims    !< The number of dimensions associated with a variable
   character(len=*), dimension(ndims), intent(in) :: dim_names !< Names of the dimensions associated with a variable
@@ -1530,8 +1574,10 @@ subroutine categorize_axes(fileObj, filename, ndims, dim_names, is_x, is_y, is_t
   if (.not.(x_found .and. y_found)) then
     ! Look for hints from CF-compliant axis units for uncharacterized axes
     do i=1,ndims ; if (.not.(is_x(i) .or. is_y(i) .or. is_t(i))) then
-      call get_variable_units(fileobj, trim(dim_names(i)), units)
-      call categorize_axis_from_units(units, is_x(i), is_y(i))
+      if (variable_exists(fileobj, trim(dim_names(i)))) then
+        call get_variable_units(fileobj, trim(dim_names(i)), units)
+        call categorize_axis_from_units(units, is_x(i), is_y(i))
+      endif
       if (is_x(i)) x_found = .true.
       if (is_y(i)) y_found = .true.
     endif ; enddo