Skip to content

*Fix PE-layout reproducibility with interior OBCs #916

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
Hallberg-NOAA merged 2 commits into from
Jun 17, 2025
Merged

*Fix PE-layout reproducibility with interior OBCs #916

Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
5d3b339
*Fix PE-layout reproducibility with interior OBCs
Hallberg-NOAA Jun 10, 2025
61a852e
+Add BT_WIDE_HALO_MIN_STENCIL & DEBUG_BT_WIDE_HALOS
Hallberg-NOAA Jun 13, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
113 changes: 74 additions & 39 deletions src/core/MOM_barotropic.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -265,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
Expand Down Expand Up @@ -756,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)

Expand Down Expand Up @@ -2395,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
Expand All @@ -2403,9 +2411,10 @@ 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
stencil = max(2, CS%min_stencil)

num_cycles = 1
if (CS%use_wide_halos) &
num_cycles = min((is-CS%isdw) / stencil, (js-CS%jsdw) / stencil)
Expand Down Expand Up @@ -2508,6 +2517,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
Expand Down Expand Up @@ -2563,22 +2580,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

Expand Down Expand Up @@ -2632,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.
Expand Down Expand Up @@ -2685,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.
Expand All @@ -2711,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.
Expand Down Expand Up @@ -3190,7 +3208,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
Expand All @@ -3216,8 +3234,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
Expand Down Expand Up @@ -3265,7 +3281,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
Expand All @@ -3290,7 +3305,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
Expand Down Expand Up @@ -3324,8 +3339,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].
Expand All @@ -3347,7 +3360,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
Expand Down Expand Up @@ -4022,6 +4034,9 @@ 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.
logical :: any_reversed_OBCs
integer :: i, j, isdw, iedw, jsdw, jedw
integer :: l_seg, Flather_OBC_in_halo

Expand Down Expand Up @@ -4121,6 +4136,15 @@ 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)))
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.

! This pair is used with either Flather or specified OBCs.
Expand Down Expand Up @@ -5428,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 "//&
Expand Down Expand Up @@ -5646,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 "//&
Expand Down