diff --git a/src/core/MOM_barotropic.F90 b/src/core/MOM_barotropic.F90
index abc2e228f6..f912ff3275 100644
--- a/src/core/MOM_barotropic.F90
+++ b/src/core/MOM_barotropic.F90
@@ -225,6 +225,9 @@ module MOM_barotropic
                              !! pressure [nondim].  Stable values are < ~1.0.
                              !! The default is 0.9.
   logical :: tides           !< If true, apply tidal momentum forcing.
+  logical :: tidal_sal_bug   !< If true, the tidal self-attraction and loading anomaly in the
+                             !! barotropic solver has the wrong sign, replicating a long-standing
+                             !! bug.
   real    :: G_extra         !< A nondimensional factor by which gtot is enhanced.
   integer :: hvel_scheme     !< An integer indicating how the thicknesses at
                              !! velocity points are calculated. Valid values are
@@ -1048,7 +1051,11 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce,
 
   if (CS%tides) then
     call tidal_forcing_sensitivity(G, CS%tides_CSp, det_de)
-    dgeo_de = 1.0 + det_de + CS%G_extra
+    if (CS%tidal_sal_bug) then
+      dgeo_de = 1.0 + det_de + CS%G_extra
+    else
+      dgeo_de = (1.0 - det_de) + CS%G_extra
+    endif
   else
     dgeo_de = 1.0 + CS%G_extra
   endif
@@ -2792,7 +2799,11 @@ subroutine set_dtbt(G, GV, US, CS, eta, pbce, BT_cont, gtot_est, SSH_add)
 
   det_de = 0.0
   if (CS%tides) call tidal_forcing_sensitivity(G, CS%tides_CSp, det_de)
-  dgeo_de = 1.0 + max(0.0, det_de + CS%G_extra)
+  if (CS%tidal_sal_bug) then
+    dgeo_de = 1.0 + max(0.0, det_de + CS%G_extra)
+  else
+    dgeo_de = 1.0 + max(0.0, CS%G_extra - det_de)
+  endif
   if (present(pbce)) then
     do j=js,je ; do i=is,ie
       gtot_E(i,j) = 0.0 ; gtot_W(i,j) = 0.0
@@ -4329,6 +4340,9 @@ subroutine barotropic_init(u, v, h, eta, Time, G, GV, US, param_file, diag, CS,
                       ! a restart file to the internal representation in this run.
   real :: mean_SL     ! The mean sea level that is used along with the bathymetry to estimate the
                       ! geometry when LINEARIZED_BT_CORIOLIS is true or BT_NONLIN_STRESS is false [Z ~> m].
+  real :: det_de      ! The partial derivative due to self-attraction and loading of the reference
+                      ! geopotential with the sea surface height when tides are enabled.
+                      ! This is typically ~0.09 or less.
   real, allocatable, dimension(:,:) :: lin_drag_h
   type(memory_size_type) :: MS
   type(group_pass_type) :: pass_static_data, pass_q_D_Cor
@@ -4483,6 +4497,14 @@ subroutine barotropic_init(u, v, h, eta, Time, G, GV, US, param_file, diag, CS,
 
   call get_param(param_file, mdl, "TIDES", CS%tides, &
                  "If true, apply tidal momentum forcing.", default=.false.)
+  det_de = 0.0
+  if (CS%tides .and. associated(CS%tides_CSp)) &
+    call tidal_forcing_sensitivity(G, CS%tides_CSp, det_de)
+  call get_param(param_file, mdl, "BAROTROPIC_TIDAL_SAL_BUG", CS%tidal_sal_bug, &
+                 "If true, the tidal self-attraction and loading anomaly in the barotropic "//&
+                 "solver has the wrong sign, replicating a long-standing bug with a scalar "//&
+                 "self-attraction and loading term or the SAL term from a previous simulation.", &
+                 default=.true., do_not_log=(det_de==0.0))
   call get_param(param_file, mdl, "SADOURNY", CS%Sadourny, &
                  "If true, the Coriolis terms are discretized with the "//&
                  "Sadourny (1975) energy conserving scheme, otherwise "//&