ALE sponge mask_z halo fixes

In the ALE sponge, the `mask_u` and `mask_v` masks are constructed from
`mask_z`, but also require valid halo data on their E/W or N/S
boundaries.

Although a `pass_var` function is called on `mask_z` before computing
these masks, this function will not populate the halos of `mask_z` if
there is no adjacent data, e.g. a non-reentrant boundary or a
land-masked tile.

And even though `mask_z` was initialized to zero, this was undone by the
internal dellocation/reallocation of the array inside of
`horiz_interp_and_extrap_tracer` (although the actual result appears to
be compiler dependent).

There are two major changes in this patch:

* The FMS-based `horiz_interp_and_extrap_tracer` function no longer does
  a deallocate/reallocate of its output arrays, and now simply assumes
  they are unallocated.

  The output arrays are also explicitly declared as intent(out).

  This change clarifies that only the compute domains of `mask_z` and
  associated fields are updated, although it doesn't fully resolve the
  issue described above.

* The ALE sponge code now explicitly initializes the halo values of
  mask_z before interpolating the mask_u and mask_v masks.

  This ensures that `mask_[uv]` boundary values are disabled on
  points where no halo data is available (and hence no halo updates from
  `pass_var`.  When the data is available, sensible values will replace
  these zeros.

These changes prevent anomalous values of mask_z from entering the
halos, and ensuring that `mask_[uv]` contain sensible values.

A similar operation should not be required by the tracer fields, since
the zero-halo values in the mask will correctly disable these values
when no adjacent field is available for halo updates.

src/framework/MOM_horizontal_regridding.F90

@@ -263,12 +263,16 @@ subroutine horiz_interp_and_extrap_tracer_record(filename, varnam,  conversion,
   real,                  intent(in)    :: conversion !< Conversion factor for tracer.
   integer,               intent(in)    :: recnum     !< Record number of tracer to be read.
   type(ocean_grid_type), intent(inout) :: G          !< Grid object
-  real, allocatable, dimension(:,:,:)  :: tr_z       !< pointer to allocatable tracer array on local
+  real, allocatable, dimension(:,:,:), intent(out) :: tr_z
+                                                     !< pointer to allocatable tracer array on local
                                                      !! model grid and input-file vertical levels.
-  real, allocatable, dimension(:,:,:)  :: mask_z     !< pointer to allocatable tracer mask array on
+  real, allocatable, dimension(:,:,:), intent(out) :: mask_z
+                                                     !< pointer to allocatable tracer mask array on
                                                      !! local model grid and input-file vertical levels.
-  real, allocatable,     dimension(:)  :: z_in       !< Cell grid values for input data.
-  real, allocatable,     dimension(:)  :: z_edges_in !< Cell grid edge values for input data.
+  real, allocatable, dimension(:), intent(out) :: z_in
+                                                     !< Cell grid values for input data.
+  real, allocatable, dimension(:), intent(out) :: z_edges_in
+                                                     !< Cell grid edge values for input data.
   real,                  intent(out)   :: missing_value !< The missing value in the returned array.
   logical,               intent(in)    :: reentrant_x !< If true, this grid is reentrant in the x-direction
   logical,               intent(in)    :: tripolar_n !< If true, this is a northern tripolar grid
@@ -329,10 +333,6 @@ subroutine horiz_interp_and_extrap_tracer_record(filename, varnam,  conversion,
   is_ongrid = .false.
   if (present(ongrid)) is_ongrid = ongrid
 
-  if (allocated(tr_z)) deallocate(tr_z)
-  if (allocated(mask_z)) deallocate(mask_z)
-  if (allocated(z_edges_in)) deallocate(z_edges_in)
-
   PI_180 = atan(1.0)/45.
 
   ! Open NetCDF file and if present, extract data and spatial coordinate information
@@ -383,13 +383,6 @@ subroutine horiz_interp_and_extrap_tracer_record(filename, varnam,  conversion,
   rcode = NF90_GET_ATT(ncid, varid, "scale_factor", scale_factor)
   if (rcode /= 0) scale_factor = 1.0
 
-  if (allocated(lon_in)) deallocate(lon_in)
-  if (allocated(lat_in)) deallocate(lat_in)
-  if (allocated(z_in)) deallocate(z_in)
-  if (allocated(z_edges_in)) deallocate(z_edges_in)
-  if (allocated(tr_z)) deallocate(tr_z)
-  if (allocated(mask_z)) deallocate(mask_z)
-
   allocate(lon_in(id), lat_in(jd), z_in(kd), z_edges_in(kd+1))
   allocate(tr_z(isd:ied,jsd:jed,kd), mask_z(isd:ied,jsd:jed,kd))
 
@@ -619,12 +612,16 @@ subroutine horiz_interp_and_extrap_tracer_fms_id(fms_id,  Time, conversion, G, t
   type(time_type),       intent(in)    :: Time       !< A FMS time type
   real,                  intent(in)    :: conversion !< Conversion factor for tracer.
   type(ocean_grid_type), intent(inout) :: G          !< Grid object
-  real, allocatable, dimension(:,:,:)  :: tr_z       !< pointer to allocatable tracer array on local
+  real, allocatable, dimension(:,:,:), intent(out) :: tr_z
+                                                     !< pointer to allocatable tracer array on local
                                                      !! model grid and native vertical levels.
-  real, allocatable, dimension(:,:,:)  :: mask_z     !< pointer to allocatable tracer mask array on
+  real, allocatable, dimension(:,:,:), intent(out) :: mask_z
+                                                     !< pointer to allocatable tracer mask array on
                                                      !! local model grid and native vertical levels.
-  real, allocatable,     dimension(:)  :: z_in       !< Cell grid values for input data.
-  real, allocatable,     dimension(:)  :: z_edges_in !< Cell grid edge values for input data. (Intent out)
+  real, allocatable, dimension(:), intent(out) :: z_in
+                                                     !< Cell grid values for input data.
+  real, allocatable, dimension(:), intent(out) :: z_edges_in
+                                                     !< Cell grid edge values for input data.
   real,                  intent(out)   :: missing_value !< The missing value in the returned array.
   logical,               intent(in)    :: reentrant_x !< If true, this grid is reentrant in the x-direction
   logical,               intent(in)    :: tripolar_n !< If true, this is a northern tripolar grid
@@ -650,8 +647,8 @@ subroutine horiz_interp_and_extrap_tracer_fms_id(fms_id,  Time, conversion, G, t
   integer :: i,j,k
   integer, dimension(4) :: start, count, dims, dim_id
   real, dimension(:,:), allocatable :: x_in, y_in
-  real, dimension(:), allocatable  :: lon_in, lat_in ! The longitude and latitude in the input file
-  real, dimension(:), allocatable  :: lat_inp ! The input file latitudes expanded to the pole
+  real, dimension(:), allocatable :: lon_in, lat_in ! The longitude and latitude in the input file
+  real, dimension(:), allocatable :: lat_inp ! The input file latitudes expanded to the pole
   real :: max_lat, min_lat, pole, max_depth, npole
   real :: roundoff  ! The magnitude of roundoff, usually ~2e-16.
   logical :: add_np
@@ -697,12 +694,6 @@ subroutine horiz_interp_and_extrap_tracer_fms_id(fms_id,  Time, conversion, G, t
   call cpu_clock_begin(id_clock_read)
 
   call get_external_field_info(fms_id, size=fld_sz, axes=axes_data, missing=missing_value)
-  if (allocated(lon_in)) deallocate(lon_in)
-  if (allocated(lat_in)) deallocate(lat_in)
-  if (allocated(z_in)) deallocate(z_in)
-  if (allocated(z_edges_in)) deallocate(z_edges_in)
-  if (allocated(tr_z)) deallocate(tr_z)
-  if (allocated(mask_z)) deallocate(mask_z)
 
   id = fld_sz(1) ; jd  = fld_sz(2) ; kd = fld_sz(3)
 
@@ -899,7 +890,6 @@ subroutine horiz_interp_and_extrap_tracer_fms_id(fms_id,  Time, conversion, G, t
         enddo
       enddo
   endif
-
 end subroutine horiz_interp_and_extrap_tracer_fms_id
 
 !> Create a 2d-mesh of grid coordinates from 1-d arrays.

src/parameterizations/vertical/MOM_ALE_sponge.F90

@@ -903,8 +903,6 @@ subroutine apply_ALE_sponge(h, dt, G, GV, US, CS, Time)
       call MOM_error(FATAL,"apply_ALE_sponge: No time information provided")
     do m=1,CS%fldno
       nz_data = CS%Ref_val(m)%nz_data
-      allocate(sp_val(G%isd:G%ied,G%jsd:G%jed,1:nz_data)); sp_val(:,:,:) = 0.0
-      allocate(mask_z(G%isd:G%ied,G%jsd:G%jed,1:nz_data)); mask_z(:,:,:) = 0.0
       call horiz_interp_and_extrap_tracer(CS%Ref_val(m)%id, Time, 1.0, G, sp_val, mask_z, z_in, &
                      z_edges_in,  missing_value, CS%reentrant_x, CS%tripolar_N, .false., &
                      spongeOnGrid=CS%SpongeDataOngrid, m_to_Z=US%m_to_Z, &
@@ -991,22 +989,20 @@ subroutine apply_ALE_sponge(h, dt, G, GV, US, CS, Time)
          call MOM_error(FATAL,"apply_ALE_sponge: No time information provided")
 
       nz_data = CS%Ref_val_u%nz_data
-      allocate(sp_val(G%isd:G%ied,G%jsd:G%jed,1:nz_data))
-      allocate(sp_val_u(G%isdB:G%iedB,G%jsd:G%jed,1:nz_data))
-      allocate(mask_u(G%isdB:G%iedB,G%jsd:G%jed,1:nz_data))
-      allocate(mask_z(G%isd:G%ied,G%jsd:G%jed,1:nz_data))
-      sp_val(:,:,:) = 0.0
-      sp_val_u(:,:,:) = 0.0
-      mask_u(:,:,:) = 0.0
-      mask_z(:,:,:) = 0.0
       ! Interpolate from the external horizontal grid and in time
       call horiz_interp_and_extrap_tracer(CS%Ref_val_u%id, Time, 1.0, G, sp_val, mask_z, z_in, &
                                           z_edges_in, missing_value, CS%reentrant_x, CS%tripolar_N, .false., &
                                           spongeOnGrid=CS%SpongeDataOngrid, m_to_Z=US%m_to_Z,&
                                           answers_2018=CS%hor_regrid_answers_2018)
 
+      ! Initialize mask_z halos to zero before pass_var, in case of no update
+      mask_z(G%isc-1, G%jsc:G%jec, :) = 0.
+      mask_z(G%iec+1, G%jsc:G%jec, :) = 0.
       call pass_var(sp_val, G%Domain)
       call pass_var(mask_z, G%Domain)
+
+      allocate(sp_val_u(G%isdB:G%iedB,G%jsd:G%jed,1:nz_data))
+      allocate(mask_u(G%isdB:G%iedB,G%jsd:G%jed,1:nz_data))
       do j=G%jsc,G%jec; do I=G%iscB,G%iecB
         sp_val_u(I,j,1:nz_data) = 0.5*(sp_val(i,j,1:nz_data)+sp_val(i+1,j,1:nz_data))
         mask_u(I,j,1:nz_data) = min(mask_z(i,j,1:nz_data),mask_z(i+1,j,1:nz_data))
@@ -1041,20 +1037,19 @@ subroutine apply_ALE_sponge(h, dt, G, GV, US, CS, Time)
       enddo
       deallocate(sp_val, sp_val_u, mask_u, mask_z, hsrc)
       nz_data = CS%Ref_val_v%nz_data
-      allocate(sp_val( G%isd:G%ied,G%jsd:G%jed,1:nz_data))
-      allocate(sp_val_v(G%isd:G%ied,G%jsdB:G%jedB,1:nz_data))
-      allocate(mask_v(G%isd:G%ied,G%jsdB:G%jedB,1:nz_data))
-      allocate(mask_z(G%isd:G%ied,G%jsd:G%jed,1:nz_data))
-      sp_val(:,:,:) = 0.0
-      sp_val_v(:,:,:) = 0.0
-      mask_z(:,:,:) = 0.0
       ! Interpolate from the external horizontal grid and in time
       call horiz_interp_and_extrap_tracer(CS%Ref_val_v%id, Time, 1.0, G, sp_val, mask_z, z_in, &
                                           z_edges_in, missing_value, CS%reentrant_x, CS%tripolar_N, .false., &
                                           spongeOnGrid=CS%SpongeDataOngrid, m_to_Z=US%m_to_Z,&
                                           answers_2018=CS%hor_regrid_answers_2018)
+      ! Initialize mask_z halos to zero before pass_var, in case of no update
+      mask_z(G%isc:G%iec, G%jsc-1, :) = 0.
+      mask_z(G%isc:G%iec, G%jec+1, :) = 0.
       call pass_var(sp_val, G%Domain)
       call pass_var(mask_z, G%Domain)
+
+      allocate(sp_val_v(G%isd:G%ied,G%jsdB:G%jedB,1:nz_data))
+      allocate(mask_v(G%isd:G%ied,G%jsdB:G%jedB,1:nz_data))
       do J=G%jscB,G%jecB; do i=G%isc,G%iec
         sp_val_v(i,J,1:nz_data) = 0.5*(sp_val(i,j,1:nz_data)+sp_val(i,j+1,1:nz_data))
         mask_v(i,J,1:nz_data) = min(mask_z(i,j,1:nz_data),mask_z(i,j+1,1:nz_data))