From 456d4a90076011deb4fa560b17fb0571aebe8a22 Mon Sep 17 00:00:00 2001 From: Marshall Ward Date: Thu, 2 Sep 2021 00:35:34 -0400 Subject: [PATCH 1/2] 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 | 46 ++++++++----------- .../vertical/MOM_ALE_sponge.F90 | 29 +++++------- 2 files changed, 30 insertions(+), 45 deletions(-) diff --git a/src/framework/MOM_horizontal_regridding.F90 b/src/framework/MOM_horizontal_regridding.F90 index d1a4b7f45d..ecfe226860 100644 --- a/src/framework/MOM_horizontal_regridding.F90 +++ b/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. diff --git a/src/parameterizations/vertical/MOM_ALE_sponge.F90 b/src/parameterizations/vertical/MOM_ALE_sponge.F90 index 419b012387..9eeb8867db 100644 --- a/src/parameterizations/vertical/MOM_ALE_sponge.F90 +++ b/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)) From b47b4463523e936753e5fa5389f2468915746ea9 Mon Sep 17 00:00:00 2001 From: Marshall Ward Date: Sat, 4 Sep 2021 16:58:02 -0400 Subject: [PATCH 2/2] ALE sponge: Only update fields on uv over masks Currently, fields from `horiz_interp_and_extrap_tracer` are interpolated onto [uv] points under the assumption that halos have sensible values. This is generally true due to the preceding `pass_var` call, but it may not be valid if the halo is not updated, e.g. along a boundary or land-masked tile. The mask is designed to be set to zero when this happens, but it could still result in an assignment of an uninitialized value to the fields (`sp_val_[uv])`. Although the value is not used, it can produce compiler warnings and errors under strict debugging builds. This patch moves the calculation of `sp_val_[uv]` so that it is only conditionally computed when the mask is set to 1. --- .../vertical/MOM_ALE_sponge.F90 | 22 ++++++++++--------- 1 file changed, 12 insertions(+), 10 deletions(-) diff --git a/src/parameterizations/vertical/MOM_ALE_sponge.F90 b/src/parameterizations/vertical/MOM_ALE_sponge.F90 index 604ac03626..0962fd2ec8 100644 --- a/src/parameterizations/vertical/MOM_ALE_sponge.F90 +++ b/src/parameterizations/vertical/MOM_ALE_sponge.F90 @@ -863,8 +863,6 @@ subroutine apply_ALE_sponge(h, dt, G, GV, US, CS, Time) real, dimension(SZK_(GV)) :: tmp_val1 ! data values remapped to model grid real, dimension(SZK_(GV)) :: h_col ! A column of thicknesses at h, u or v points [H ~> m or kg m-2] real, allocatable, dimension(:,:,:) :: sp_val ! A temporary array for fields - real, allocatable, dimension(:,:,:) :: sp_val_u ! A temporary array for fields - real, allocatable, dimension(:,:,:) :: sp_val_v ! A temporary array for fields real, allocatable, dimension(:,:,:) :: mask_z ! A temporary array for field mask at h pts real, allocatable, dimension(:,:,:) :: mask_u ! A temporary array for field mask at u pts real, allocatable, dimension(:,:,:) :: mask_v ! A temporary array for field mask at v pts @@ -883,6 +881,8 @@ subroutine apply_ALE_sponge(h, dt, G, GV, US, CS, Time) real :: Idt ! The inverse of the timestep [T-1 ~> s-1] real :: h_neglect, h_neglect_edge ! Negligible thicknesses [H ~> m or kg m-2] real :: zTopOfCell, zBottomOfCell ! Interface heights (positive upward) in the input dataset [Z ~> m]. + real :: sp_val_u ! Interpolation of sp_val to u-points + real :: sp_val_v ! Interpolation of sp_val to v-points integer :: nPoints is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke @@ -1001,10 +1001,8 @@ subroutine apply_ALE_sponge(h, dt, G, GV, US, CS, Time) 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)) enddo ; enddo @@ -1014,7 +1012,10 @@ subroutine apply_ALE_sponge(h, dt, G, GV, US, CS, Time) ! Therefore we use c as per C code and increment the index where necessary. i = CS%col_i_u(c) ; j = CS%col_j_u(c) if (mask_u(i,j,1) == 1.0) then - CS%Ref_val_u%p(1:nz_data,c) = sp_val_u(i,j,1:nz_data) + do k=1,nz_data + sp_val_u = 0.5 * (sp_val(i,j,k) + sp_val(i+1,j,k)) + CS%Ref_val_u%p(k,c) = sp_val_u + enddo else CS%Ref_val_u%p(1:nz_data,c) = 0.0 endif @@ -1035,7 +1036,7 @@ subroutine apply_ALE_sponge(h, dt, G, GV, US, CS, Time) hsrc(nz_data) = hsrc(nz_data) + ( zTopOfCell + G%bathyT(i,j) ) CS%Ref_val_u%h(1:nz_data,c) = GV%Z_to_H*hsrc(1:nz_data) enddo - deallocate(sp_val, sp_val_u, mask_u, mask_z, hsrc) + deallocate(sp_val, mask_u, mask_z, hsrc) nz_data = CS%Ref_val_v%nz_data ! 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, & @@ -1048,10 +1049,8 @@ subroutine apply_ALE_sponge(h, dt, G, GV, US, CS, Time) 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)) enddo ; enddo !call pass_var(mask_z,G%Domain) @@ -1061,7 +1060,10 @@ subroutine apply_ALE_sponge(h, dt, G, GV, US, CS, Time) ! Therefore we use c as per C code and increment the index where necessary. i = CS%col_i_v(c) ; j = CS%col_j_v(c) if (mask_v(i,j,1) == 1.0) then - CS%Ref_val_v%p(1:nz_data,c) = sp_val_v(i,j,1:nz_data) + do k=1,nz_data + sp_val_v = 0.5 * (sp_val(i,j,k) + sp_val(i,j+1,k)) + CS%Ref_val_v%p(k,c) = sp_val_v + enddo else CS%Ref_val_v%p(1:nz_data,c) = 0.0 endif @@ -1082,7 +1084,7 @@ subroutine apply_ALE_sponge(h, dt, G, GV, US, CS, Time) hsrc(nz_data) = hsrc(nz_data) + ( zTopOfCell + G%bathyT(i,j) ) CS%Ref_val_v%h(1:nz_data,c) = GV%Z_to_H*hsrc(1:nz_data) enddo - deallocate(sp_val, sp_val_v, mask_v, mask_z, hsrc) + deallocate(sp_val, mask_v, mask_z, hsrc) endif call pass_var(h,G%Domain)