Fix indexing bugs in get_field_nc (#334)

Fixed two horizontal indexing bugs in `get_field_nc`, where the difference
between the array starting index (always 1 in this subroutine) and the values
in the handle argument were not being taken into account when the array was
being passed in with only its computational domain.

Also initialized the internal `unlim_index` array in `get_netcdf_fields` to fix
a problem with using an uninitialized array that was being flagged when run in
debug mode.

With this commit, the model is once again reproducing the expected answers when
rescaling is applied for vertical distances or time.

Co-authored-by: Marshall Ward <marshall.ward@gmail.com>

src/framework/MOM_io.F90

@@ -2057,23 +2057,24 @@ 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.
+!> 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)
   character(len=*), intent(in) :: filename
     !< Input filename
   character(len=*), intent(in)  :: fieldname
     !< Field variable name
-  real, intent(out) :: values(:,:)
-    !< Field value
+  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
+    !< Rescale factor, omitting this is the same as setting it to 1.
 
   integer :: turns
     ! Number of quarter-turns from input to model grid

src/framework/MOM_io_file.F90

@@ -1681,15 +1681,18 @@ subroutine read_field_chksum_nc(handle, field, chksum, valid_chksum)
 end subroutine read_field_chksum_nc
 
 
-!> Read the values of a netCDF field
+!> Read the values of a netCDF field into an array that might have halos
 subroutine get_field_nc(handle, label, values, rescale)
   class(MOM_netcdf_file), intent(in) :: handle
-    ! Handle of netCDF file to be read
+    !< Handle of netCDF file to be read
   character(len=*), intent(in) :: label
-    ! Field variable name
-  real, intent(out) :: values(:,:)
-    ! Field values read from file
+    !< 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.
   real, optional, intent(in) :: rescale
+    !< A multiplicative rescaling factor for the values that are read.
+    !! Omitting this is the same as setting it to 1.
 
   logical :: data_domain
     ! True if values matches the data domain size
@@ -1701,10 +1704,12 @@ subroutine get_field_nc(handle, label, values, rescale)
     ! Index bounds of compute domain
   integer :: isd, ied, jsd, jed
     ! Index bounds of data domain
+  integer :: iscl, iecl, jscl, jecl
+    ! Local 1-based index bounds of compute domain
   integer :: bounds(2,2)
     ! Index bounds of domain
-  real, allocatable :: values_nc(:,:)
-    ! Local copy of field, used for data-decomposed I/O
+  real, allocatable :: values_c(:,:)
+    ! Field values on the compute domain, used for copying to a data domain
 
   isc = handle%HI%isc
   iec = handle%HI%iec
@@ -1727,36 +1732,45 @@ subroutine get_field_nc(handle, label, values, rescale)
 
   field_nc = handle%fields%get(label)
 
-  if (data_domain) then
-    allocate(values_nc(isc:iec,jsc:jec))
-    values(:,:) = 0.
-  endif
+  if (data_domain) &
+    allocate(values_c(1:iec-isc+1,1:jec-jsc+1))
 
   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)
+      call read_netcdf_field(handle%handle_nc, field_nc, values_c, 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)
+      call read_netcdf_field(handle%handle_nc, field_nc, values_c)
     else
       call read_netcdf_field(handle%handle_nc, field_nc, values)
     endif
   endif
 
-  if (data_domain) &
-    values(isc:iec,jsc:jec) = values_nc(:,:)
+  if (data_domain) then
+    iscl = isc - isd + 1
+    iecl = iec - isd + 1
+    jscl = jsc - jsd + 1
+    jecl = jec - jsd + 1
+
+    values(iscl:iecl,jscl:jecl) = values_c(:,:)
+  else
+    iscl = 1
+    iecl = iec - isc + 1
+    jscl = 1
+    jecl = jec - jsc + 1
+  endif
 
   ! NOTE: It is more efficient to do the rescale in-place while copying
-  ! values_nc(:,:) to values(:,:).  But since rescale is only present for
+  ! values_c(:,:) 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)
+      values(iscl:iecl,jscl:jecl) = rescale * values(iscl:iecl,jscl:jecl)
     endif
   endif
 end subroutine get_field_nc

src/framework/MOM_netcdf.F90

@@ -665,6 +665,9 @@ subroutine get_netcdf_fields(handle, axes, fields)
   call check_netcdf_call(rc, 'get_netcdf_fields', &
       'File "' // trim(handle%filename) // '"')
 
+  ! Initialize unlim_index with an unreachable value (outside [1,ndims])
+  unlim_index = -1
+
   allocate(axes(ndims))
   do i = 1, ndims
     rc = nf90_inquire_dimension(handle%ncid, dimids(i), name=label, len=len)