Fms2 io read3d slice (#399)

* Restore functionality for reading slices from 3d volumes in MOM_io

 - The recent MOM_io modifications in support of FMS2_io accidentally
   removed support for reading on-grid data (same horizontal grid as
   model) k-slices. This is needed in some configurations in the model
   state initialization.

* Add FMS1 interfaces

* Additional patches to enable reading ongrid state initialization data

 - read local 3d volume rather than attempting to slice ongrid
   data vertically.

 - Related bugfixes in MOM_io

config_src/infra/FMS1/MOM_coms_infra.F90

@@ -25,7 +25,7 @@ module MOM_coms_infra
 !> Communicate an array, string or scalar from one PE to others
 interface broadcast
   module procedure broadcast_char, broadcast_int32_0D, broadcast_int64_0D, broadcast_int1D
-  module procedure broadcast_real0D, broadcast_real1D, broadcast_real2D
+  module procedure broadcast_real0D, broadcast_real1D, broadcast_real2D, broadcast_real3D
 end interface broadcast
 
 !> Compute a checksum for a field distributed over a PE list.  If no PE list is
@@ -260,6 +260,28 @@ subroutine broadcast_real2D(dat, length, from_PE, PElist, blocking)
 
 end subroutine broadcast_real2D
 
+
+!> Communicate a 3-D array of reals from one PE to others
+subroutine broadcast_real3D(dat, length, from_PE, PElist, blocking)
+  real, dimension(:,:,:), intent(inout) :: dat       !< The data to communicate and destination
+  integer,              intent(in)      :: length    !< The total number of data elements
+  integer,    optional, intent(in)      :: from_PE   !< The source PE, by default the root PE
+  integer,    optional, intent(in)      :: PElist(:) !< The list of participating PEs, by default the
+                                                     !! active PE set as previously set via Set_PElist.
+  logical,    optional, intent(in)      :: blocking  !< If true, barriers are added around the call
+
+  integer :: src_PE   ! The processor that is sending the data
+  logical :: do_block ! If true add synchronizing barriers
+
+  do_block = .false. ; if (present(blocking)) do_block = blocking
+  if (present(from_PE)) then ; src_PE = from_PE ; else ; src_PE = root_PE() ; endif
+
+  if (do_block) call mpp_sync(PElist)
+  call mpp_broadcast(dat, length, src_PE, PElist)
+  if (do_block) call mpp_sync_self(PElist)
+
+end subroutine broadcast_real3D
+
 ! field_chksum wrappers
 
 !> Compute a checksum for a field distributed over a PE list.  If no PE list is

config_src/infra/FMS1/MOM_io_infra.F90

@@ -57,7 +57,7 @@ module MOM_io_infra
 !> Read a data field from a file
 interface read_field
   module procedure read_field_4d
-  module procedure read_field_3d
+  module procedure read_field_3d, read_field_3d_region
   module procedure read_field_2d, read_field_2d_region
   module procedure read_field_1d, read_field_1d_int
   module procedure read_field_0d, read_field_0d_int
@@ -696,6 +696,45 @@ subroutine read_field_3d(filename, fieldname, data, MOM_Domain, &
   endif ; endif
 end subroutine read_field_3d
 
+!> This routine uses the fms_io subroutine read_data to read a region from a distributed or
+!! global 3-D data field named "fieldname" from file "filename".
+subroutine read_field_3d_region(filename, fieldname, data, start, nread, MOM_domain, &
+                                no_domain, 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
+  integer, dimension(:),  intent(in)    :: start     !< The starting index to read in each of 4
+                                                     !! dimensions.  For this 3-d read, the
+                                                     !! 4th values are always 1.
+  integer, dimension(:),  intent(in)    :: nread     !< The number of points to read in each of 4
+                                                     !! dimensions.  For this 3-d read, the
+                                                     !! 4th values are always 1.
+  type(MOM_domain_type), &
+                optional, intent(in)    :: MOM_Domain !< The MOM_Domain that describes the decomposition
+  logical,      optional, intent(in)    :: no_domain !< If present and true, this variable does not
+                                                     !! use domain decomposion.
+  real,         optional, intent(in)    :: scale     !< A scaling factor that the field is multiplied
+                                                     !! by before it is returned.
+
+  if (present(MOM_Domain)) then
+    call read_data(filename, fieldname, data, start, nread, domain=MOM_Domain%mpp_domain, &
+                   no_domain=no_domain)
+  else
+    call read_data(filename, fieldname, data, start, nread, no_domain=no_domain)
+  endif
+
+  if (present(scale)) then ; if (scale /= 1.0) then
+    if (present(MOM_Domain)) then
+      call rescale_comp_data(MOM_Domain, data, scale)
+    else
+      ! Dangerously rescale the whole array
+      data(:,:,:) = scale*data(:,:,:)
+    endif
+  endif ; endif
+end subroutine read_field_3d_region
+
+
 !> This routine uses the fms_io subroutine 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.

config_src/infra/FMS2/MOM_coms_infra.F90

@@ -25,7 +25,7 @@ module MOM_coms_infra
 !> Communicate an array, string or scalar from one PE to others
 interface broadcast
   module procedure broadcast_char, broadcast_int32_0D, broadcast_int64_0D, broadcast_int1D
-  module procedure broadcast_real0D, broadcast_real1D, broadcast_real2D
+  module procedure broadcast_real0D, broadcast_real1D, broadcast_real2D, broadcast_real3D
 end interface broadcast
 
 !> Compute a checksum for a field distributed over a PE list.  If no PE list is
@@ -260,6 +260,27 @@ subroutine broadcast_real2D(dat, length, from_PE, PElist, blocking)
 
 end subroutine broadcast_real2D
 
+!> Communicate a 3-D array of reals from one PE to others
+subroutine broadcast_real3D(dat, length, from_PE, PElist, blocking)
+  real, dimension(:,:,:), intent(inout) :: dat       !< The data to communicate and destination
+  integer,              intent(in)    :: length    !< The total number of data elements
+  integer,    optional, intent(in)    :: from_PE   !< The source PE, by default the root PE
+  integer,    optional, intent(in)    :: PElist(:) !< The list of participating PEs, by default the
+                                                   !! active PE set as previously set via Set_PElist.
+  logical,    optional, intent(in)    :: blocking  !< If true, barriers are added around the call
+
+  integer :: src_PE   ! The processor that is sending the data
+  logical :: do_block ! If true add synchronizing barriers
+
+  do_block = .false. ; if (present(blocking)) do_block = blocking
+  if (present(from_PE)) then ; src_PE = from_PE ; else ; src_PE = root_PE() ; endif
+
+  if (do_block) call mpp_sync(PElist)
+  call mpp_broadcast(dat, length, src_PE, PElist)
+  if (do_block) call mpp_sync_self(PElist)
+
+end subroutine broadcast_real3D
+
 ! field_chksum wrappers
 
 !> Compute a checksum for a field distributed over a PE list.  If no PE list is

config_src/infra/FMS2/MOM_io_infra.F90

@@ -66,7 +66,7 @@ module MOM_io_infra
 !> Read a data field from a file
 interface read_field
   module procedure read_field_4d
-  module procedure read_field_3d
+  module procedure read_field_3d, read_field_3d_region
   module procedure read_field_2d, read_field_2d_region
   module procedure read_field_1d, read_field_1d_int
   module procedure read_field_0d, read_field_0d_int
@@ -1030,6 +1030,74 @@ subroutine read_field_3d(filename, fieldname, data, MOM_Domain, &
 
 end subroutine read_field_3d
 
+!> This routine uses the fms_io subroutine read_data to read a region from a distributed or
+!! global 3-D data field named "fieldname" from file "filename".
+subroutine read_field_3d_region(filename, fieldname, data, start, nread, MOM_domain, &
+                                no_domain, 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
+  integer, dimension(:),  intent(in)    :: start     !< The starting index to read in each of 3
+                                                     !! dimensions.  For this 3-d read, the
+                                                     !! 4th value is always 1.
+  integer, dimension(:),  intent(in)    :: nread     !< The number of points to read in each of 4
+                                                     !! dimensions.  For this 3-d read, the
+                                                     !! 4th values are always 1.
+  type(MOM_domain_type), &
+                optional, intent(in)    :: MOM_Domain !< The MOM_Domain that describes the decomposition
+  logical,      optional, intent(in)    :: no_domain !< If present and true, this variable does not
+                                                     !! use domain decomposion.
+  real,         optional, intent(in)    :: scale     !< A scaling factor that the field is multiplied
+                                                     !! by before it is returned.
+
+  ! Local variables
+  type(FmsNetcdfFile_t)       :: fileObj ! A handle to a non-domain-decomposed file
+  type(FmsNetcdfDomainFile_t) :: fileobj_DD ! A handle to a domain-decomposed file object
+  character(len=96) :: var_to_read ! Name of variable to read from the netcdf file
+  logical :: success               ! True if the file was successfully opened
+
+  if (present(MOM_Domain)) then
+    ! Open the FMS2 file-set.
+    success = fms2_open_file(fileobj_DD, filename, "read", MOM_domain%mpp_domain)
+    if (.not.success) call MOM_error(FATAL, "Failed to open "//trim(filename))
+
+    ! Find the matching case-insensitive variable name in the file and prepare to read it.
+    call prepare_to_read_var(fileobj_DD, fieldname, "read_field_2d_region: ", &
+                             filename, var_to_read)
+
+    ! Read the data.
+    call fms2_read_data(fileobj_DD, var_to_read, data, corner=start(1:3), edge_lengths=nread(1:3))
+
+    ! Close the file-set.
+    if (check_if_open(fileobj_DD)) call fms2_close_file(fileobj_DD)
+  else
+    ! Open the FMS2 file-set.
+    success = fms2_open_file(fileObj, trim(filename), "read")
+    if (.not.success) call MOM_error(FATAL, "Failed to open "//trim(filename))
+
+    ! Find the matching case-insensitive variable name in the file, and determine whether it
+    ! has a time dimension.
+    call find_varname_in_file(fileObj, fieldname, "read_field_2d_region: ", filename, var_to_read)
+
+    ! Read the data.
+    call fms2_read_data(fileobj, var_to_read, data, corner=start(1:3), edge_lengths=nread(1:3))
+
+    ! Close the file-set.
+    if (check_if_open(fileobj)) call fms2_close_file(fileobj)
+  endif
+
+  if (present(scale)) then ; if (scale /= 1.0) then
+    if (present(MOM_Domain)) then
+      call rescale_comp_data(MOM_Domain, data, scale)
+    else
+      ! Dangerously rescale the whole array
+      data(:,:,:) = scale*data(:,:,:)
+    endif
+  endif ; endif
+
+end subroutine read_field_3d_region
+
 !> This routine uses the fms_io subroutine 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.

src/framework/MOM_horizontal_regridding.F90

@@ -309,6 +309,9 @@ subroutine horiz_interp_and_extrap_tracer_record(filename, varnam, recnum, G, tr
   real, dimension(:,:),  allocatable   :: tr_in      !< A 2-d array for holding input data on its
                                                      !! native horizontal grid, with units that change
                                                      !! as the input data is interpreted [a] then [A ~> a]
+  real, dimension(:,:,:),  allocatable   :: tr_in_full  !< A 3-d array for holding input data on the
+                                                     !! model horizontal grid, with units that change
+                                                     !! as the input data is interpreted [a] then [A ~> a]
   real, dimension(:,:),  allocatable   :: tr_inp     !< Native horizontal grid data extended to the poles
                                                      !! with units that change as the input data is
                                                      !! interpreted [a] then [A ~> a]
@@ -448,6 +451,7 @@ subroutine horiz_interp_and_extrap_tracer_record(filename, varnam, recnum, G, tr
 
   if (is_ongrid) then
     allocate(tr_in(is:ie,js:je), source=0.0)
+    allocate(tr_in_full(is:ie,js:je,kd), source=0.0)
     allocate(mask_in(is:ie,js:je), source=0.0)
   else
     call horizontal_interp_init()
@@ -470,14 +474,19 @@ subroutine horiz_interp_and_extrap_tracer_record(filename, varnam, recnum, G, tr
 
   ! Loop through each data level and interpolate to model grid.
   ! After interpolating, fill in points which will be needed to define the layers.
+
+  if (is_ongrid) then
+    start(1) = is+G%HI%idg_offset ; start(2) = js+G%HI%jdg_offset ; start(3) = 1
+    count(1) = ie-is+1 ; count(2) = je-js+1 ; count(3) = kd ; start(4) = 1 ; count(4) = 1
+    call MOM_read_data(trim(filename), trim(varnam), tr_in_full, start, count, G%Domain)
+  endif
+
   do k=1,kd
     mask_in(:,:)  = 0.0
     tr_out(:,:) = 0.0
 
     if (is_ongrid) then
-      start(1) = is+G%HI%idg_offset ; start(2) = js+G%HI%jdg_offset ; start(3) = k
-      count(1) = ie-is+1 ; count(2) = je-js+1 ; count(3) = 1 ; start(4) = 1 ; count(4) = 1
-      call MOM_read_data(trim(filename), trim(varnam), tr_in, start, count, G%Domain)
+      tr_in(is:ie,js:je) = tr_in_full(is:ie,js:je,k)
       do j=js,je
         do i=is,ie
           if (abs(tr_in(i,j)-missing_val_in) > abs(roundoff*missing_val_in)) then
@@ -594,7 +603,10 @@ subroutine horiz_interp_and_extrap_tracer_record(filename, varnam, recnum, G, tr
 
   enddo ! kd
 
-  deallocate(lon_in, lat_in)
+  if (allocated(lat_inp)) deallocate(lat_inp)
+  deallocate(tr_in)
+  if (allocated(tr_inp)) deallocate(tr_inp)
+  if (allocated(tr_in_full)) deallocate(tr_in_full)
 
 end subroutine horiz_interp_and_extrap_tracer_record