From 5d3b3397f6334f36abb9c7934a62c3a9938146bc Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Tue, 10 Jun 2025 13:43:47 -0400 Subject: [PATCH 1/2] *Fix PE-layout reproducibility with interior OBCs Corrected the previous PE-count and layout dependency of answers for cases with open boundaries are not at the natural edges of the domain by storing the previous barotropic velocities in ubt_prev and vbt_prev over a wider range of loop indices, and by increasing the stencil size used in btstep_timeloop in such cases. Because the required loop range no longer matches the range over which the new velocities are set in bt_loop_update_u() and bt_loop_update_v(), the code setting these arrays has been moved out of these routines and into the main body of btstep_timeloop(). This commit will change answers (and make them invariant to the PE count and layout) in some cases, such as the ESMG-configs Kelvin_wave/rotate45_BT and Kelvin_wave/rotate_BT test cases, where there are southern open boundaries appearing in the northern halo regions on any PEs, for example. In all the Kelvin_wave test cases, the answers now agree with the previous single-PE answers. Most regional configurations with OBCs only apply them at the edges of the domain and hence will not be impacted, and no answers are changed in cases without open boundary conditions. --- src/core/MOM_barotropic.F90 | 43 ++++++++++++++++++++++++++----------- 1 file changed, 31 insertions(+), 12 deletions(-) diff --git a/src/core/MOM_barotropic.F90 b/src/core/MOM_barotropic.F90 index 4688e646c5..46ffdf871b 100644 --- a/src/core/MOM_barotropic.F90 +++ b/src/core/MOM_barotropic.F90 @@ -4,6 +4,7 @@ module MOM_barotropic ! This file is part of MOM6. See LICENSE.md for the license. use MOM_checksums, only : chksum0 +use MOM_coms, only : any_across_PEs use MOM_cpu_clock, only : cpu_clock_id, cpu_clock_begin, cpu_clock_end, CLOCK_ROUTINE use MOM_debugging, only : hchksum, uvchksum use MOM_diag_mediator, only : post_data, query_averaging_enabled, register_diag_field @@ -92,6 +93,9 @@ module MOM_barotropic integer, allocatable :: v_OBC_type(:,:) !< An integer encoding the type and direction of v-point OBCs logical :: u_OBCs_on_PE !< True if this PE has an open boundary at any u-points. logical :: v_OBCs_on_PE !< True if this PE has an open boundary at any v-points. + logical :: wide_stencil !< If true, a wider stencil must be used in the barotropic iteration because + !! of the presence of OBCs in the opposite halo (e.g., southern OBCs in a + !! northern halo) on any processors anywhere in the domain. This is uncommon. !>@{ Index ranges on the local PE for the open boundary conditions in various directions integer :: Is_u_W_obc, Ie_u_W_obc, js_u_W_obc, je_u_W_obc integer :: Is_u_E_obc, Ie_u_E_obc, js_u_E_obc, je_u_E_obc @@ -2406,6 +2410,8 @@ subroutine btstep_timeloop(eta, ubt, vbt, uhbt0, Datu, BTCL_u, vhbt0, Datv, BTCL stencil = 1 if ((.not.use_BT_cont) .and. CS%Nonlinear_continuity .and. (CS%Nonlin_cont_update_period > 0)) & stencil = 2 + if (CS%BT_OBC%wide_stencil) stencil = 2 + num_cycles = 1 if (CS%use_wide_halos) & num_cycles = min((is-CS%isdw) / stencil, (js-CS%jsdw) / stencil) @@ -2508,6 +2514,14 @@ subroutine btstep_timeloop(eta, ubt, vbt, uhbt0, Datu, BTCL_u, vhbt0, Datv, BTCL jsv = jsv+stencil ; jev = jev-stencil endif + ! Store the previous velocities for time-filtered transports and OBCs. + do j=jsv,jev ; do I=isv-2,iev+1 + ubt_prev(I,j) = ubt(I,j) + enddo ; enddo + do J=jsv-2,jev+1 ; do i=isv,iev + vbt_prev(i,J) = vbt(i,J) + enddo ; enddo + if (integral_BT_cont) then !$OMP parallel do default(shared) do j=jsv-1,jev+1 ; do I=isv-2,iev+1 @@ -2563,22 +2577,22 @@ subroutine btstep_timeloop(eta, ubt, vbt, uhbt0, Datu, BTCL_u, vhbt0, Datv, BTCL if (v_first) then ! On odd-steps, update v first. call btloop_update_v(dtbt, ubt, vbt, v_accel_bt, Cor_v, PFv, isv-1, iev+1, jsv-1, jev, & - f_4_v, bt_rem_v, BT_force_v, vbt_prev, Cor_ref_v, Rayleigh_v, & + f_4_v, bt_rem_v, BT_force_v, Cor_ref_v, Rayleigh_v, & wt_accel(n), G, US, CS) ! Now update the zonal velocity. call btloop_update_u(dtbt, ubt, vbt, u_accel_bt, Cor_u, PFu, isv-1, iev, jsv, jev, & - f_4_u, bt_rem_u, BT_force_u, ubt_prev, Cor_ref_u, Rayleigh_u, & + f_4_u, bt_rem_u, BT_force_u, Cor_ref_u, Rayleigh_u, & wt_accel(n), G, US, CS) else ! On even steps, update u first. call btloop_update_u(dtbt, ubt, vbt, u_accel_bt, Cor_u, PFu, isv-1, iev, jsv-1, jev+1, & - f_4_u, bt_rem_u, BT_force_u, ubt_prev, Cor_ref_u, Rayleigh_u, & + f_4_u, bt_rem_u, BT_force_u, Cor_ref_u, Rayleigh_u, & wt_accel(n), G, US, CS) ! Now update the meridional velocity. call btloop_update_v(dtbt, ubt, vbt, v_accel_bt, Cor_v, PFv, isv, iev, jsv-1, jev, & - f_4_v, bt_rem_v, BT_force_v, vbt_prev, Cor_ref_v, Rayleigh_v, & + f_4_v, bt_rem_v, BT_force_v, Cor_ref_v, Rayleigh_v, & wt_accel(n), G, US, CS, Cor_bracket_bug=CS%use_old_coriolis_bracket_bug) endif @@ -3190,7 +3204,7 @@ end subroutine btloop_add_dyn_PF !> Update meridional velocity. subroutine btloop_update_v(dtbt, ubt, vbt, v_accel_bt, & Cor_v, PFv, is_v, ie_v, Js_v, Je_v, f_4_v, & - bt_rem_v, BT_force_v, vbt_prev, Cor_ref_v, Rayleigh_v, & + bt_rem_v, BT_force_v, Cor_ref_v, Rayleigh_v, & wt_accel_n, G, US, CS, Cor_bracket_bug) type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure. type(barotropic_CS), intent(inout) :: CS !< Barotropic control structure @@ -3216,8 +3230,6 @@ subroutine btloop_update_v(dtbt, ubt, vbt, v_accel_bt, & integer, intent(in) :: ie_v !< The ending i-index of the range of v-point values to calculate integer, intent(in) :: Js_v !< The starting j-index of the range of v-point values to calculate integer, intent(in) :: Je_v !< The ending j-index of the range of v-point values to calculate - real, dimension(SZIW_(CS),SZJBW_(CS)), intent(inout) :: & - vbt_prev !< The previous velocity, stored for time-filtered transports and OBCs [L T-1 ~> m s-1] real, dimension(SZIW_(CS),SZJBW_(CS)), intent(in) :: & bt_rem_v !< The fraction of the barotropic meridional velocity that !! remains after a time step, the rest being lost to bottom @@ -3265,7 +3277,6 @@ subroutine btloop_update_v(dtbt, ubt, vbt, v_accel_bt, & !$OMP do schedule(static) ! This updates the v-velocity, except at OBC points. do J=Js_v,Je_v ; do i=is_v,ie_v - vbt_prev(i,J) = vbt(i,J) vbt(i,J) = bt_rem_v(i,J) * (vbt(i,J) + & dtbt * ((BT_force_v(i,J) + Cor_v(i,J)) + PFv(i,J))) if (abs(vbt(i,J)) < CS%vel_underflow) vbt(i,J) = 0.0 @@ -3290,7 +3301,7 @@ end subroutine btloop_update_v !> Update zonal velocity. subroutine btloop_update_u(dtbt, ubt, vbt, u_accel_bt, & Cor_u, PFu, Is_u, Ie_u, js_u, je_u, f_4_u, & - bt_rem_u, BT_force_u, ubt_prev, Cor_ref_u, Rayleigh_u, & + bt_rem_u, BT_force_u, Cor_ref_u, Rayleigh_u, & wt_accel_n, G, US, CS) type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure. type(barotropic_CS), intent(inout) :: CS !< Barotropic control structure @@ -3324,8 +3335,6 @@ subroutine btloop_update_u(dtbt, ubt, vbt, u_accel_bt, & real, dimension(SZIBW_(CS),SZJW_(CS)), intent(in) :: & BT_force_u !< The vertical average of all of the v-accelerations that are !! not explicitly included in the barotropic equation [L T-2 ~> m s-2]. - real, dimension(SZIBW_(CS),SZJW_(CS)), intent(inout) :: & - ubt_prev !< The previous velocity, stored for time-filtered transports and OBCs [L T-1 ~> m s-1] real, dimension(SZIBW_(CS),SZJW_(CS)), intent(in) :: & Cor_ref_u !< The meridional barotropic Coriolis acceleration due !! to the reference velocities [L T-2 ~> m s-2]. @@ -3347,7 +3356,6 @@ subroutine btloop_update_u(dtbt, ubt, vbt, u_accel_bt, & ((f_4_u(3,I,j) * vbt(i,J)) + (f_4_u(2,I,j) * vbt(i+1,J-1)))) - & Cor_ref_u(I,j) - ubt_prev(I,j) = ubt(I,j) ubt(I,j) = bt_rem_u(I,j) * (ubt(I,j) + & dtbt * ((BT_force_u(I,j) + Cor_u(I,j)) + PFu(I,j))) if (abs(ubt(I,j)) < CS%vel_underflow) ubt(I,j) = 0.0 @@ -4022,6 +4030,8 @@ subroutine initialize_BT_OBC(OBC, BT_OBC, G, CS) ! converted to real numbers to work with the MOM6 halo update code [nondim] real :: OBC_sign ! A sign encoding the direction of the OBC being used at a point [nondim] real :: OBC_type ! A real copy of the integer encoding the type of OBC being used at a point [nondim] + logical :: reversed_OBCs ! True of there any OBCs in the opposite halo on this PE, e.g. points + ! with a southern OBC in a northern halo. integer :: i, j, isdw, iedw, jsdw, jedw integer :: l_seg, Flather_OBC_in_halo @@ -4121,6 +4131,14 @@ subroutine initialize_BT_OBC(OBC, BT_OBC, G, CS) BT_OBC%u_OBCs_on_PE = ((BT_OBC%Is_u_E_obc <= iedw) .or. (BT_OBC%Is_u_W_obc <= iedw)) BT_OBC%v_OBCs_on_PE = ((BT_OBC%is_v_N_obc <= iedw) .or. (BT_OBC%is_v_S_obc <= iedw)) + + ! Determine whether there are any OBCs in the opposite halo on any processors in the domain, e.g., + ! points with OBC_DIRECTION_S in a northern halo. + reversed_OBCs = (BT_OBC%u_OBCs_on_PE .and. ((BT_OBC%Is_u_E_obc <= G%isc-1) .or. (BT_OBC%Ie_u_W_obc >= G%iec))) .or. & + (BT_OBC%v_OBCs_on_PE .and. ((BT_OBC%Js_v_N_obc <= G%jsc-1) .or. (BT_OBC%Je_v_S_obc >= G%jec))) + BT_OBC%wide_stencil = any_across_PEs(reversed_OBCs) + if (BT_OBC%wide_stencil) call MOM_mesg("OBCs in an opposite halo require the use of a wider stencil.", 3) + ! Allocate time-varying arrays that will be used for open boundary conditions. ! This pair is used with either Flather or specified OBCs. @@ -5804,6 +5822,7 @@ subroutine barotropic_init(u, v, h, Time, G, GV, US, param_file, diag, CS, & else ! There are no OBC points anywhere. CS%BT_OBC%u_OBCs_on_PE = .false. CS%BT_OBC%v_OBCs_on_PE = .false. + CS%BT_OBC%wide_stencil = .false. CS%integral_OBCs = .false. endif From 61a852e62641f3fb008af94d8a8074aa118844b1 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Fri, 13 Jun 2025 12:14:12 -0400 Subject: [PATCH 2/2] +Add BT_WIDE_HALO_MIN_STENCIL & DEBUG_BT_WIDE_HALOS This commit adds the new runtime parameters BT_WIDE_HALO_MIN_STENCIL and DEBUG_BT_WIDE_HALOS to control the stencil width used with wide halos and to specify whether the checksums for the entire active halo of values are written when DEBUG_BT is enabled, or just the symmetric computational domain. Setting `BT_WIDE_HALO_MIN_STENCIL = 2` and `DEBUG_BT_WIDE_HALOS = False` provides for cleaner checksums in some cases because they are invariant to the use of wide halos, the width of those wide halos, or the relative positions of interior- domain OBCs relative to PE boundaries. This commit also changes a number of velocity point checksums in btstep_timeloop to do symmetric output. All model solutions are bitwise identical regardless of how these new parameters are set, but they can alter some debugging output. --- src/core/MOM_barotropic.F90 | 84 ++++++++++++++++++++++--------------- 1 file changed, 50 insertions(+), 34 deletions(-) diff --git a/src/core/MOM_barotropic.F90 b/src/core/MOM_barotropic.F90 index 46ffdf871b..3042449314 100644 --- a/src/core/MOM_barotropic.F90 +++ b/src/core/MOM_barotropic.F90 @@ -93,9 +93,6 @@ module MOM_barotropic integer, allocatable :: v_OBC_type(:,:) !< An integer encoding the type and direction of v-point OBCs logical :: u_OBCs_on_PE !< True if this PE has an open boundary at any u-points. logical :: v_OBCs_on_PE !< True if this PE has an open boundary at any v-points. - logical :: wide_stencil !< If true, a wider stencil must be used in the barotropic iteration because - !! of the presence of OBCs in the opposite halo (e.g., southern OBCs in a - !! northern halo) on any processors anywhere in the domain. This is uncommon. !>@{ Index ranges on the local PE for the open boundary conditions in various directions integer :: Is_u_W_obc, Ie_u_W_obc, js_u_W_obc, je_u_W_obc integer :: Is_u_E_obc, Ie_u_E_obc, js_u_E_obc, je_u_E_obc @@ -269,12 +266,18 @@ module MOM_barotropic !! velocities. The streaming band-pass filter must be turned on. logical :: use_wide_halos !< If true, use wide halos and march in during the !! barotropic time stepping for efficiency. + integer :: min_stencil !< The minimum stencil width to use with the wide halo iterations. + !! A nonzero value may reflect the distribution of OBC faces or it + !! may be useful for debugging purposes. logical :: clip_velocity !< If true, limit any velocity components that are !! are large enough for a CFL number to exceed !! CFL_trunc. This should only be used as a !! desperate debugging measure. logical :: debug !< If true, write verbose checksums for debugging purposes. - logical :: debug_bt !< If true, write verbose checksums for debugging purposes. + logical :: debug_bt !< If true, write verbose checksums from within the barotropic + !! time-stepping loop for debugging purposes. + logical :: debug_wide_halos !< If true, write the checksums on the full wide halos. Otherwise + !! only the output for the final computational domain is written. real :: vel_underflow !< Velocity components smaller than vel_underflow !! are set to 0 [L T-1 ~> m s-1]. real :: maxvel !< Velocity components greater than maxvel are @@ -760,9 +763,9 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce, interp_eta_PF = associated(eta_PF_start) ! Figure out the fullest arrays that could be updated. - stencil = 1 + stencil = max(1, CS%min_stencil) if ((.not.use_BT_cont) .and. CS%Nonlinear_continuity .and. & - (CS%Nonlin_cont_update_period > 0)) stencil = 2 + (CS%Nonlin_cont_update_period > 0)) stencil = max(2, CS%min_stencil) find_etaav = present(etaav) @@ -2399,6 +2402,7 @@ subroutine btstep_timeloop(eta, ubt, vbt, uhbt0, Datu, BTCL_u, vhbt0, Datv, BTCL integer :: stencil ! The stencil size of the algorithm, often 1 or 2. integer :: err_count ! A counter to limit the volume of error messages written to stdout. integer :: i, j, n, is, ie, js, je + integer :: debug_halo ! The halo size to use for debugging checksums integer :: isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec @@ -2407,10 +2411,9 @@ subroutine btstep_timeloop(eta, ubt, vbt, uhbt0, Datu, BTCL_u, vhbt0, Datv, BTCL err_count = 0 ! Figure out the fullest arrays that could be updated. - stencil = 1 + stencil = max(1, CS%min_stencil) if ((.not.use_BT_cont) .and. CS%Nonlinear_continuity .and. (CS%Nonlin_cont_update_period > 0)) & - stencil = 2 - if (CS%BT_OBC%wide_stencil) stencil = 2 + stencil = max(2, CS%min_stencil) num_cycles = 1 if (CS%use_wide_halos) & @@ -2646,23 +2649,24 @@ subroutine btstep_timeloop(eta, ubt, vbt, uhbt0, Datu, BTCL_u, vhbt0, Datv, BTCL ! This might need to be moved outside of the OMP do loop directives. if (CS%debug_bt) then write(mesg,'("BT vel update ",I4)') n - call uvchksum(trim(mesg)//" PF[uv]", PFu, PFv, CS%debug_BT_HI, haloshift=iev-ie, & - unscale=US%L_T_to_m_s*US%s_to_T) - call uvchksum(trim(mesg)//" Cor_[uv]", Cor_u, Cor_v, CS%debug_BT_HI, haloshift=iev-ie, & - unscale=US%L_T_to_m_s*US%s_to_T) - call uvchksum(trim(mesg)//" BT_force_[uv]", BT_force_u, BT_force_v, CS%debug_BT_HI, haloshift=iev-ie, & - unscale=US%L_T_to_m_s*US%s_to_T) - call uvchksum(trim(mesg)//" BT_rem_[uv]", BT_rem_u, BT_rem_v, CS%debug_BT_HI, & - haloshift=iev-ie, scalar_pair=.true.) - call uvchksum(trim(mesg)//" [uv]bt", ubt, vbt, CS%debug_BT_HI, haloshift=iev-ie, & - unscale=US%L_T_to_m_s) - call uvchksum(trim(mesg)//" [uv]bt_trans", ubt_trans, vbt_trans, CS%debug_BT_HI, haloshift=iev-ie, & - unscale=US%L_T_to_m_s) - call uvchksum(trim(mesg)//" [uv]hbt", uhbt, vhbt, CS%debug_BT_HI, haloshift=iev-ie, & - unscale=US%s_to_T*US%L_to_m**2*GV%H_to_m) + debug_halo = 0 ; if (CS%debug_wide_halos) debug_halo = iev - ie + call uvchksum(trim(mesg)//" PF[uv]", PFu, PFv, CS%debug_BT_HI, haloshift=debug_halo, & + symmetric=.true., unscale=US%L_T_to_m_s*US%s_to_T) + call uvchksum(trim(mesg)//" Cor_[uv]", Cor_u, Cor_v, CS%debug_BT_HI, haloshift=debug_halo, & + symmetric=.true., unscale=US%L_T_to_m_s*US%s_to_T) + call uvchksum(trim(mesg)//" BT_force_[uv]", BT_force_u, BT_force_v, CS%debug_BT_HI, haloshift=debug_halo, & + symmetric=.true., unscale=US%L_T_to_m_s*US%s_to_T) + call uvchksum(trim(mesg)//" BT_rem_[uv]", BT_rem_u, BT_rem_v, CS%debug_BT_HI, haloshift=debug_halo, & + symmetric=.true., scalar_pair=.true.) + call uvchksum(trim(mesg)//" [uv]bt", ubt, vbt, CS%debug_BT_HI, haloshift=debug_halo, & + symmetric=.true., unscale=US%L_T_to_m_s) + call uvchksum(trim(mesg)//" [uv]bt_trans", ubt_trans, vbt_trans, CS%debug_BT_HI, haloshift=debug_halo, & + symmetric=.true., unscale=US%L_T_to_m_s) + call uvchksum(trim(mesg)//" [uv]hbt", uhbt, vhbt, CS%debug_BT_HI, haloshift=debug_halo, & + symmetric=.true., unscale=US%s_to_T*US%L_to_m**2*GV%H_to_m) if (integral_BT_cont) & - call uvchksum(trim(mesg)//" [uv]hbt_int", uhbt_int, vhbt_int, CS%debug_BT_HI, haloshift=iev-ie, & - unscale=US%L_to_m**2*GV%H_to_m) + call uvchksum(trim(mesg)//" [uv]hbt_int", uhbt_int, vhbt_int, CS%debug_BT_HI, haloshift=debug_halo, & + symmetric=.true., unscale=US%L_to_m**2*GV%H_to_m) endif ! Apply open boundary condition considerations to revise the updated velocities and transports. @@ -2699,11 +2703,11 @@ subroutine btstep_timeloop(eta, ubt, vbt, uhbt0, Datu, BTCL_u, vhbt0, Datv, BTCL !$OMP end do nowait if (CS%debug_bt) then - call uvchksum("BT [uv]hbt just after OBC", uhbt, vhbt, CS%debug_BT_HI, haloshift=iev-ie, & - unscale=US%s_to_T*US%L_to_m**2*GV%H_to_m) + call uvchksum("BT [uv]hbt just after OBC", uhbt, vhbt, CS%debug_BT_HI, haloshift=debug_halo, & + symmetric=.true., unscale=US%s_to_T*US%L_to_m**2*GV%H_to_m) if (integral_BT_cont) & call uvchksum("BT [uv]hbt_int just after OBC", uhbt_int, vhbt_int, CS%debug_BT_HI, & - haloshift=iev-ie, unscale=US%L_to_m**2*GV%H_to_m) + haloshift=debug_halo, symmetric=.true., unscale=US%L_to_m**2*GV%H_to_m) endif ! Update eta in a corrector step using the barotropic continuity equation. @@ -2725,9 +2729,9 @@ subroutine btstep_timeloop(eta, ubt, vbt, uhbt0, Datu, BTCL_u, vhbt0, Datv, BTCL if (CS%debug_bt) then write(mesg,'("BT step ",I4)') n - call uvchksum(trim(mesg)//" [uv]bt", ubt, vbt, CS%debug_BT_HI, haloshift=iev-ie, & - unscale=US%L_T_to_m_s) - call hchksum(eta, trim(mesg)//" eta", CS%debug_BT_HI, haloshift=iev-ie, unscale=GV%H_to_MKS) + call uvchksum(trim(mesg)//" [uv]bt", ubt, vbt, CS%debug_BT_HI, haloshift=debug_halo, & + symmetric=.true., unscale=US%L_T_to_m_s) + call hchksum(eta, trim(mesg)//" eta", CS%debug_BT_HI, haloshift=debug_halo, unscale=GV%H_to_MKS) endif ! Issue warnings if there are unphysical values of the sea surface height or total water column mass. @@ -4032,6 +4036,7 @@ subroutine initialize_BT_OBC(OBC, BT_OBC, G, CS) real :: OBC_type ! A real copy of the integer encoding the type of OBC being used at a point [nondim] logical :: reversed_OBCs ! True of there any OBCs in the opposite halo on this PE, e.g. points ! with a southern OBC in a northern halo. + logical :: any_reversed_OBCs integer :: i, j, isdw, iedw, jsdw, jedw integer :: l_seg, Flather_OBC_in_halo @@ -4136,8 +4141,9 @@ subroutine initialize_BT_OBC(OBC, BT_OBC, G, CS) ! points with OBC_DIRECTION_S in a northern halo. reversed_OBCs = (BT_OBC%u_OBCs_on_PE .and. ((BT_OBC%Is_u_E_obc <= G%isc-1) .or. (BT_OBC%Ie_u_W_obc >= G%iec))) .or. & (BT_OBC%v_OBCs_on_PE .and. ((BT_OBC%Js_v_N_obc <= G%jsc-1) .or. (BT_OBC%Je_v_S_obc >= G%jec))) - BT_OBC%wide_stencil = any_across_PEs(reversed_OBCs) - if (BT_OBC%wide_stencil) call MOM_mesg("OBCs in an opposite halo require the use of a wider stencil.", 3) + any_reversed_OBCs = any_across_PEs(reversed_OBCs) + if (any_reversed_OBCs) call MOM_mesg("OBCs in an opposite halo require the use of a wider stencil.", 5) + if (any_reversed_OBCs) CS%min_stencil = max(CS%min_stencil, 2) ! Allocate time-varying arrays that will be used for open boundary conditions. @@ -5446,6 +5452,11 @@ subroutine barotropic_init(u, v, h, Time, G, GV, US, param_file, diag, CS, & call get_param(param_file, mdl, "BTHALO", bt_halo_sz, & "The minimum halo size for the barotropic solver.", default=0, & layoutParam=.true.) + call get_param(param_file, mdl, "BT_WIDE_HALO_MIN_STENCIL", CS%min_stencil, & + "The minimum stencil width to use with the wide halo iterations. "//& + "A nonzero value may be useful for debugging purposes, but at the "//& + "cost of reducing the efficiency gain from BT_USE_WIDE_HALOS.", & + default=0, layoutParam=.true., do_not_log=.not.CS%use_wide_halos) #ifdef STATIC_MEMORY_ if ((bt_halo_sz > 0) .and. (bt_halo_sz /= BTHALO_)) call MOM_error(FATAL, & "barotropic_init: Run-time values of BTHALO must agree with the "//& @@ -5664,6 +5675,12 @@ subroutine barotropic_init(u, v, h, Time, G, GV, US, param_file, diag, CS, & "barotropic time-stepping loop. The data volume can be "//& "quite large if this is true.", default=CS%debug, & debuggingParam=.true.) + call get_param(param_file, mdl, "DEBUG_BT_WIDE_HALOS", CS%debug_wide_halos, & + "If true, write the checksums on the full wide halos. Otherwise only the "//& + "output for the final computational domain is written. This can be valuable "//& + "for debugging certain cases where the stencil used in the wide halo "//& + "iterations depends on which opoen boundary conditions are in the halos.", & + default=.true., do_not_log=.not.(CS%debug_bt.and.CS%use_wide_halos), debuggingParam=.true.) call get_param(param_file, mdl, "LINEARIZED_BT_CORIOLIS", CS%linearized_BT_PV, & "If true use the bottom depth instead of the total water column thickness "//& @@ -5822,7 +5839,6 @@ subroutine barotropic_init(u, v, h, Time, G, GV, US, param_file, diag, CS, & else ! There are no OBC points anywhere. CS%BT_OBC%u_OBCs_on_PE = .false. CS%BT_OBC%v_OBCs_on_PE = .false. - CS%BT_OBC%wide_stencil = .false. CS%integral_OBCs = .false. endif