diff --git a/src/ALE/MOM_ALE.F90 b/src/ALE/MOM_ALE.F90
index 2f763c7827..da93c37c10 100644
--- a/src/ALE/MOM_ALE.F90
+++ b/src/ALE/MOM_ALE.F90
@@ -1835,7 +1835,7 @@ subroutine ALE_initThicknessToCoord( CS, G, GV, h, height_units )
   scale = GV%Z_to_H
   if (present(height_units)) then ; if (height_units) scale = 1.0 ; endif
   do j = G%jsd,G%jed ; do i = G%isd,G%ied
-    h(i,j,:) = scale * getStaticThickness( CS%regridCS, 0., G%bathyT(i,j)+G%Z_ref )
+    h(i,j,:) = scale * getStaticThickness( CS%regridCS, 0., max(G%meanSL(i,j)+G%bathyT(i,j), 0.0) )
   enddo ; enddo
 
 end subroutine ALE_initThicknessToCoord
diff --git a/src/ALE/MOM_hybgen_unmix.F90 b/src/ALE/MOM_hybgen_unmix.F90
index d87762b721..5d8b8c9c10 100644
--- a/src/ALE/MOM_hybgen_unmix.F90
+++ b/src/ALE/MOM_hybgen_unmix.F90
@@ -219,18 +219,18 @@ subroutine hybgen_unmix(G, GV, US, CS, tv, Reg, ntr, h)
       do k=1,nk
         dz_tot = dz_tot + GV%H_to_RZ * tv%SpV_avg(i,j,k) * h_col(k)
       enddo
-      if (dz_tot <= CS%min_dilate*(G%bathyT(i,j)+G%Z_ref)) then
+      if (dz_tot <= CS%min_dilate * (G%meanSL(i,j) + G%bathyT(i,j))) then
         dilate = CS%min_dilate
-      elseif (dz_tot >= CS%max_dilate*(G%bathyT(i,j)+G%Z_ref)) then
+      elseif (dz_tot >= CS%max_dilate * (G%meanSL(i,j) + G%bathyT(i,j))) then
         dilate = CS%max_dilate
       else
-        dilate = dz_tot / (G%bathyT(i,j)+G%Z_ref)
+        dilate = dz_tot / (G%meanSL(i,j) + G%bathyT(i,j))
       endif
     else
-      nominalDepth = (G%bathyT(i,j)+G%Z_ref)*GV%Z_to_H
-      if (h_tot <= CS%min_dilate*nominalDepth) then
+      nominalDepth = (G%meanSL(i,j) + G%bathyT(i,j)) * GV%Z_to_H
+      if (h_tot <= CS%min_dilate * nominalDepth) then
         dilate = CS%min_dilate
-      elseif (h_tot >= CS%max_dilate*nominalDepth) then
+      elseif (h_tot >= CS%max_dilate * nominalDepth) then
         dilate = CS%max_dilate
       else
         dilate = h_tot / nominalDepth
diff --git a/src/ALE/MOM_regridding.F90 b/src/ALE/MOM_regridding.F90
index 872766bbab..2564b9a3bf 100644
--- a/src/ALE/MOM_regridding.F90
+++ b/src/ALE/MOM_regridding.F90
@@ -871,7 +871,7 @@ subroutine initialize_regridding(CS, G, GV, US, max_depth, param_file, mdl, &
       endif
       do i=G%isc-1,G%iec+1; do j=G%jsc-1,G%jec+1
         if (G%mask2dT(i,j)>0.) then
-          nominalDepth = (G%bathyT(i,j)+G%Z_ref)*US%Z_to_m
+          nominalDepth = max(G%meanSL(i,j) + G%bathyT(i,j), 0.0) * US%Z_to_m
           if (nominalDepth <= depth_s) then
             do k= 1,n_sigma
               dz_3d(i,j,k) = dz_shallow(k)
@@ -1251,15 +1251,15 @@ subroutine regridding_main( remapCS, CS, G, GV, US, h, tv, h_new, dzInterface, &
       tot_dz(i,j) = tot_dz(i,j) + GV%H_to_RZ * tv%SpV_avg(i,j,k) * h(i,j,k)
     enddo ; enddo ; enddo
     do j=G%jsc-1,G%jec+1 ; do i=G%isc-1,G%iec+1
-      if ((tot_dz(i,j) > 0.0) .and. (G%bathyT(i,j)+G%Z_ref > 0.0)) then
-        nom_depth_H(i,j) = (G%bathyT(i,j)+G%Z_ref) * (tot_h(i,j) / tot_dz(i,j))
+      if (tot_dz(i,j) > 0.0) then
+        nom_depth_H(i,j) = max(G%meanSL(i,j) + G%bathyT(i,j), 0.0) * (tot_h(i,j) / tot_dz(i,j))
       else
         nom_depth_H(i,j) = 0.0
       endif
     enddo ; enddo
   else
     do j=G%jsc-1,G%jec+1 ; do i=G%isc-1,G%iec+1
-      nom_depth_H(i,j) = max((G%bathyT(i,j)+G%Z_ref) * Z_to_H, 0.0)
+      nom_depth_H(i,j) = max(G%meanSL(i,j) + G%bathyT(i,j), 0.0) * Z_to_H
     enddo ; enddo
   endif
 
diff --git a/src/core/MOM.F90 b/src/core/MOM.F90
index c2b8104827..2350791f4d 100644
--- a/src/core/MOM.F90
+++ b/src/core/MOM.F90
@@ -4285,8 +4285,8 @@ subroutine extract_surface_state(CS, sfc_state_in)
     do j=js,je ; do i=is,ie
       if (G%mask2dT(i,j)>0.) then
         localError = sfc_state%sea_lev(i,j) < -G%bathyT(i,j) - G%Z_ref &
-                .or. sfc_state%sea_lev(i,j) >=  CS%bad_val_ssh_max  &
-                .or. sfc_state%sea_lev(i,j) <= -CS%bad_val_ssh_max  &
+                .or. sfc_state%sea_lev(i,j) >=  CS%bad_val_ssh_max + (G%meanSL(i,j) - G%Z_ref) &
+                .or. sfc_state%sea_lev(i,j) <= -CS%bad_val_ssh_max + (G%meanSL(i,j) - G%Z_ref) &
                 .or. sfc_state%sea_lev(i,j) + G%bathyT(i,j) + G%Z_ref < CS%bad_val_col_thick
         if (use_temperature) localError = localError &
                 .or. sfc_state%SSS(i,j)<0.                        &
diff --git a/src/core/MOM_PressureForce_FV.F90 b/src/core/MOM_PressureForce_FV.F90
index 0f5e536c38..6a16392c3a 100644
--- a/src/core/MOM_PressureForce_FV.F90
+++ b/src/core/MOM_PressureForce_FV.F90
@@ -438,8 +438,9 @@ subroutine PressureForce_FV_nonBouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, AD
     else
       !$OMP parallel do default(shared)
       do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
-        SSH(i,j) = (za(i,j,1) - alpha_ref*p(i,j,1)) * I_gEarth - G%Z_ref &
-                  - max(-G%bathyT(i,j)-G%Z_ref, 0.0)
+        SSH(i,j) = (za(i,j,1) - alpha_ref*p(i,j,1)) * I_gEarth - G%Z_ref
+        ! Remove above sea level topography at floodable cells
+        SSH(i,j) = SSH(i,j) - max(-G%bathyT(i,j)-G%meanSL(i,j), 0.0)
       enddo ; enddo
       call calc_SAL(SSH, e_sal, G, CS%SAL_CSp, tmp_scale=US%Z_to_m)
     endif
@@ -1163,7 +1164,7 @@ subroutine PressureForce_FV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, ADp,
     !$OMP parallel do default(shared)
     do j=Jsq,Jeq+1
       do i=Isq,Ieq+1
-        SSH(i,j) = min(-G%bathyT(i,j) - G%Z_ref, 0.0)
+        SSH(i,j) = min(-G%bathyT(i,j) - G%meanSL(i,j), 0.0)
       enddo
       do k=1,nz ; do i=Isq,Ieq+1
         SSH(i,j) = SSH(i,j) + h(i,j,k)*GV%H_to_Z
@@ -1275,7 +1276,7 @@ subroutine PressureForce_FV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, ADp,
     enddo ; enddo
   else
     do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
-      Z_0p(i,j) = G%Z_ref
+      Z_0p(i,j) = G%meanSL(i,j)
     enddo ; enddo
   endif
 
@@ -1359,7 +1360,9 @@ subroutine PressureForce_FV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, ADp,
     else
       !$OMP parallel do default(shared)
       do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
-        SSH(i,j) = e(i,j,1) - max(-G%bathyT(i,j) - G%Z_ref, 0.0) ! Remove topography above sea level
+        SSH(i,j) = e(i,j,1) - G%Z_ref
+        ! Remove above sea level topography at floodable cells
+        SSH(i,j) = SSH(i,j) - max(-G%bathyT(i,j)-G%meanSL(i,j), 0.0)
       enddo ; enddo
       call calc_SAL(SSH, e_sal, G, CS%SAL_CSp, tmp_scale=US%Z_to_m)
     endif
diff --git a/src/core/MOM_PressureForce_Montgomery.F90 b/src/core/MOM_PressureForce_Montgomery.F90
index 1529af9d83..0098470502 100644
--- a/src/core/MOM_PressureForce_Montgomery.F90
+++ b/src/core/MOM_PressureForce_Montgomery.F90
@@ -197,7 +197,7 @@ subroutine PressureForce_Mont_nonBouss(h, tv, PFu, PFv, G, GV, US, CS, p_atm, pb
     ! of self-attraction and loading.
     !$OMP parallel do default(shared)
     do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
-      SSH(i,j) = min(-G%bathyT(i,j) - G%Z_ref, 0.0)
+      SSH(i,j) = min(-G%bathyT(i,j) - G%meanSL(i,j), 0.0)
     enddo ; enddo
     if (use_EOS) then
       !$OMP parallel do default(shared)
@@ -476,7 +476,7 @@ subroutine PressureForce_Mont_Bouss(h, tv, PFu, PFv, G, GV, US, CS, p_atm, pbce,
     ! barotropic tides.
     !$OMP parallel do default(shared)
     do j=Jsq,Jeq+1
-      do i=Isq,Ieq+1 ; SSH(i,j) = min(-G%bathyT(i,j) - G%Z_ref, 0.0) ; enddo
+      do i=Isq,Ieq+1 ; SSH(i,j) = min(-G%bathyT(i,j) - G%meanSL(i,j), 0.0) ; enddo
       do k=1,nz ; do i=Isq,Ieq+1
         SSH(i,j) = SSH(i,j) + h(i,j,k)*GV%H_to_Z
       enddo ; enddo
@@ -707,7 +707,7 @@ subroutine Set_pbce_Bouss(e, tv, G, GV, US, Rho0, GFS_scale, pbce, rho_star)
       do j=Jsq,Jeq+1
         do i=Isq,Ieq+1
           Ihtot(i) = GV%H_to_Z / ((e(i,j,1)-e(i,j,nz+1)) + dz_neglect)
-          press(i) = -Rho0xG*(e(i,j,1) - G%Z_ref)
+          press(i) = -Rho0xG*(e(i,j,1) - G%meanSL(i,j))
         enddo
         call calculate_density(tv%T(:,j,1), tv%S(:,j,1), press, rho_in_situ, &
                                tv%eqn_of_state, EOSdom)
@@ -716,7 +716,7 @@ subroutine Set_pbce_Bouss(e, tv, G, GV, US, Rho0, GFS_scale, pbce, rho_star)
         enddo
         do k=2,nz
           do i=Isq,Ieq+1
-            press(i) = -Rho0xG*(e(i,j,K) - G%Z_ref)
+            press(i) = -Rho0xG*(e(i,j,K) - G%meanSL(i,j))
             T_int(i) = 0.5*(tv%T(i,j,k-1)+tv%T(i,j,k))
             S_int(i) = 0.5*(tv%S(i,j,k-1)+tv%S(i,j,k))
           enddo
diff --git a/src/core/MOM_barotropic.F90 b/src/core/MOM_barotropic.F90
index 2c7e19be58..156ab6234c 100644
--- a/src/core/MOM_barotropic.F90
+++ b/src/core/MOM_barotropic.F90
@@ -5373,14 +5373,14 @@ subroutine find_face_areas(Datu, Datv, G, GV, US, CS, MS, halo, eta, add_max)
         H1 = CS%bathyT(i,j)*GV%Z_to_H + eta(i,j) ; H2 = CS%bathyT(i+1,j)*GV%Z_to_H + eta(i+1,j)
         Datu(I,j) = 0.0 ; if ((H1 > 0.0) .and. (H2 > 0.0)) &
         Datu(I,j) = CS%dy_Cu(I,j) * (2.0 * H1 * H2) / (H1 + H2)
-!       Datu(I,j) = CS%dy_Cu(I,j) * 0.5 * (H1 + H2)
+        ! Datu(I,j) = CS%dy_Cu(I,j) * 0.5 * (H1 + H2)
       enddo ; enddo
       !$OMP do
       do J=js-1-hs,je+hs ; do i=is-hs,ie+hs
         H1 = CS%bathyT(i,j)*GV%Z_to_H + eta(i,j) ; H2 = CS%bathyT(i,j+1)*GV%Z_to_H + eta(i,j+1)
         Datv(i,J) = 0.0 ; if ((H1 > 0.0) .and. (H2 > 0.0)) &
         Datv(i,J) = CS%dx_Cv(i,J) * (2.0 * H1 * H2) / (H1 + H2)
-!       Datv(i,J) = CS%dy_v(i,J) * 0.5 * (H1 + H2)
+        ! Datv(i,J) = CS%dy_v(i,J) * 0.5 * (H1 + H2)
       enddo ; enddo
     else
       !$OMP do
@@ -5403,27 +5403,31 @@ subroutine find_face_areas(Datu, Datv, G, GV, US, CS, MS, halo, eta, add_max)
 
     !$OMP do
     do j=js-hs,je+hs ; do I=is-1-hs,ie+hs
-      Datu(I,j) = CS%dy_Cu(I,j) * Z_to_H * &
-                 max(max(CS%bathyT(i+1,j), CS%bathyT(i,j)) + (G%Z_ref + add_max), 0.0)
+      H1 = max((G%meanSL(i+1,j) + add_max) + G%bathyT(i+1,j), 0.0)
+      H2 = max((G%meanSL(i,j) + add_max) + G%bathyT(i,j), 0.0)
+      Datu(I,j) = CS%dy_Cu(I,j) * Z_to_H * max(H1, H2)
     enddo ; enddo
     !$OMP do
     do J=js-1-hs,je+hs ; do i=is-hs,ie+hs
-      Datv(i,J) = CS%dx_Cv(i,J) * Z_to_H * &
-                 max(max(CS%bathyT(i,j+1), CS%bathyT(i,j)) + (G%Z_ref + add_max), 0.0)
+      H1 = max((G%meanSL(i,j+1) + add_max) + G%bathyT(i,j+1), 0.0)
+      H2 = max((G%meanSL(i,j) + add_max) + G%bathyT(i,j), 0.0)
+      Datv(i,J) = CS%dx_Cv(i,J) * Z_to_H * max(H1, H2)
     enddo ; enddo
   else
     Z_to_H = GV%Z_to_H ; if (.not.GV%Boussinesq) Z_to_H = GV%RZ_to_H * CS%Rho_BT_lin
 
     !$OMP do
     do j=js-hs,je+hs ; do I=is-1-hs,ie+hs
-      H1 = (CS%bathyT(i,j) + G%Z_ref) * Z_to_H ; H2 = (CS%bathyT(i+1,j) + G%Z_ref) * Z_to_H
+      H1 = max(G%meanSL(i,j) + G%bathyT(i,j), 0.0) * Z_to_H
+      H2 = max(G%meanSL(i+1,j) + G%bathyT(i+1,j), 0.0) * Z_to_H
       Datu(I,j) = 0.0
       if ((H1 > 0.0) .and. (H2 > 0.0)) &
         Datu(I,j) = CS%dy_Cu(I,j) * (2.0 * H1 * H2) / (H1 + H2)
     enddo ; enddo
     !$OMP do
     do J=js-1-hs,je+hs ; do i=is-hs,ie+hs
-      H1 = (CS%bathyT(i,j) + G%Z_ref) * Z_to_H ; H2 = (CS%bathyT(i,j+1) + G%Z_ref) * Z_to_H
+      H1 = max(G%meanSL(i,j) + G%bathyT(i,j), 0.0) * Z_to_H
+      H2 = max(G%meanSL(i,j+1) + G%bathyT(i,j+1), 0.0) * Z_to_H
       Datv(i,J) = 0.0
       if ((H1 > 0.0) .and. (H2 > 0.0)) &
         Datv(i,J) = CS%dx_Cv(i,J) * (2.0 * H1 * H2) / (H1 + H2)
@@ -5546,8 +5550,6 @@ subroutine barotropic_init(u, v, h, Time, G, GV, US, param_file, diag, CS, &
                                        ! name in wave_drag_file.
   character(len=80)  :: wave_drag_v    ! The wave drag piston velocity variable
                                        ! name in wave_drag_file.
-  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 :: htot        ! Total column thickness used when BT_NONLIN_STRESS is false [Z ~> m].
   real :: Z_to_H      ! A local unit conversion factor [H Z-1 ~> nondim or kg m-3]
   real :: H_to_Z      ! A local unit conversion factor [Z H-1 ~> nondim or m3 kg-1]
@@ -6150,25 +6152,26 @@ subroutine barotropic_init(u, v, h, Time, G, GV, US, param_file, diag, CS, &
 
     Z_to_H = GV%Z_to_H ; if (.not.GV%Boussinesq) Z_to_H = GV%RZ_to_H * CS%Rho_BT_lin
 
-    Mean_SL = G%Z_ref
     do j=js,je ; do I=is-1,ie
-      CS%D_u_Cor(I,j) = 0.5 * (max(Mean_SL+G%bathyT(i+1,j),0.0) + max(Mean_SL+G%bathyT(i,j),0.0)) * Z_to_H
+      CS%D_u_Cor(I,j) = 0.5 * ( max(G%meanSL(i+1,j) + G%bathyT(i+1,j), 0.0) &
+                              + max(G%meanSL(i,j) + G%bathyT(i,j), 0.0) ) * Z_to_H
     enddo ; enddo
     if (CS%interior_OBC_PV .and. CS%BT_OBC%u_OBCs_on_PE) then ; do j=js,je ; do I=is-1,ie
       if (CS%BT_OBC%u_OBC_type(I,j) < 0) & ! Western boundary condition
-        CS%D_u_Cor(I,j) = max(Mean_SL+G%bathyT(i+1,j),0.0) * Z_to_H
+        CS%D_u_Cor(I,j) = max(G%meanSL(i+1,j) + G%bathyT(i+1,j), 0.0) * Z_to_H
       if (CS%BT_OBC%u_OBC_type(I,j) > 0) & ! Eastern boundary condition
-        CS%D_u_Cor(I,j) = max(Mean_SL+G%bathyT(i,j),0.0) * Z_to_H
+        CS%D_u_Cor(I,j) = max(G%meanSL(i,j) + G%bathyT(i,j), 0.0) * Z_to_H
     enddo ; enddo ; endif
 
     do J=js-1,je ; do i=is,ie
-      CS%D_v_Cor(i,J) = 0.5 * (max(Mean_SL+G%bathyT(i,j+1),0.0) + max(Mean_SL+G%bathyT(i,j),0.0)) * Z_to_H
+      CS%D_v_Cor(i,J) = 0.5 * ( max(G%meanSL(i,j+1) + G%bathyT(i,j+1), 0.0) &
+                              + max(G%meanSL(i,j) + G%bathyT(i,j), 0.0) )  * Z_to_H
     enddo ; enddo
     if (CS%interior_OBC_PV .and. CS%BT_OBC%v_OBCs_on_PE) then ; do J=js-1,je ; do i=is,ie
       if (CS%BT_OBC%v_OBC_type(i,J) < 0) & ! Southern boundary condition
-        CS%D_v_Cor(i,J) = max(Mean_SL+G%bathyT(i,j+1),0.0) * Z_to_H
+        CS%D_v_Cor(i,J) = max(G%meanSL(i,j+1) + G%bathyT(i,j+1), 0.0) * Z_to_H
       if (CS%BT_OBC%v_OBC_type(i,J) > 0) & ! Northern boundary condition
-        CS%D_v_Cor(i,J) = max(Mean_SL+G%bathyT(i,j),0.0) * Z_to_H
+        CS%D_v_Cor(i,J) = max(G%meanSL(i,j) + G%bathyT(i,j), 0.0) * Z_to_H
     enddo ; enddo ; endif
 
     h_a_neglect = GV%H_subroundoff * 1.0 * US%m_to_L**2
@@ -6176,10 +6179,11 @@ subroutine barotropic_init(u, v, h, Time, G, GV, US, param_file, diag, CS, &
       if ((CS%q_wt(1,I,J) + CS%q_wt(4,I,J)) + (CS%q_wt(2,I,J) + CS%q_wt(3,I,J)) > 0.) then
         CS%q_D(I,J) = 0.25 * (CS%BT_Coriolis_scale * G%CoriolisBu(I,J)) * &
            ((CS%q_wt(1,I,J) + CS%q_wt(4,I,J)) + (CS%q_wt(2,I,J) + CS%q_wt(3,I,J))) / &
-           max(Z_to_H * (((CS%q_wt(1,I,J) * max(Mean_SL+G%bathyT(i,j),0.0)) + &
-                          (CS%q_wt(4,I,J) * max(Mean_SL+G%bathyT(i+1,j+1),0.0))) + &
-                         ((CS%q_wt(2,I,J) * max(Mean_SL+G%bathyT(i+1,j),0.0)) + &
-                          (CS%q_wt(3,I,J) * max(Mean_SL+G%bathyT(i,j+1),0.0)))), h_a_neglect)
+           max(Z_to_H * (((CS%q_wt(1,I,J) * max(G%meanSL(i,j) + G%bathyT(i,j), 0.0)) + &
+                          (CS%q_wt(4,I,J) * max(G%meanSL(i+1,j+1) + G%bathyT(i+1,j+1), 0.0))) + &
+                         ((CS%q_wt(2,I,J) * max(G%meanSL(i+1,j) + G%bathyT(i+1,j), 0.0)) + &
+                          (CS%q_wt(3,I,J) * max(G%meanSL(i,j+1) + G%bathyT(i,j+1), 0.0)))), &
+               h_a_neglect)
       else ! All four h points are masked out so q_D(I,J) is meaningless
         CS%q_D(I,J) = 0.
       endif
@@ -6450,10 +6454,9 @@ subroutine barotropic_init(u, v, h, Time, G, GV, US, param_file, diag, CS, &
   ! Calculate other constants which are used for btstep.
 
   if (.not.CS%nonlin_stress) then
-    Mean_SL = G%Z_ref
     Z_to_H = GV%Z_to_H ; if (.not.GV%Boussinesq) Z_to_H = GV%RZ_to_H * CS%Rho_BT_lin
     do j=js,je ; do I=is-1,ie
-      htot = max(G%bathyT(i+1,j) + G%Z_ref, 0.0) + max(G%bathyT(i,j) + G%Z_ref, 0.0)
+      htot = max(G%meanSL(i+1,j) + G%bathyT(i+1,j), 0.0) + max(G%meanSL(i,j) + G%bathyT(i,j), 0.0)
       if (G%OBCmaskCu(I,j) * htot > 0.) then
         CS%IDatu(I,j) = G%OBCmaskCu(I,j) * 2.0 / (Z_to_H * htot)
       else ! Both neighboring H points are masked out or this is an OBC face so IDatu(I,j) is unused
@@ -6461,7 +6464,7 @@ subroutine barotropic_init(u, v, h, Time, G, GV, US, param_file, diag, CS, &
       endif
     enddo ; enddo
     do J=js-1,je ; do i=is,ie
-      htot = max(G%bathyT(i,j+1) + G%Z_ref, 0.0) + max(G%bathyT(i,j) + G%Z_ref, 0.0)
+      htot = max(G%meanSL(i,j+1) + G%bathyT(i,j+1), 0.0) + max(G%meanSL(i,j) + G%bathyT(i,j), 0.0)
       if (G%OBCmaskCv(i,J) * htot > 0.) then
         CS%IDatv(i,J) = G%OBCmaskCv(i,J) * 2.0 / (Z_to_H * htot)
       else ! Both neighboring H points are masked out or this is an OBC face so IDatv(i,J) is unused
diff --git a/src/core/MOM_grid.F90 b/src/core/MOM_grid.F90
index 94583673c2..4e8c1a9cc2 100644
--- a/src/core/MOM_grid.F90
+++ b/src/core/MOM_grid.F90
@@ -158,7 +158,16 @@ module MOM_grid
     y_ax_unit_short     !< A short description of the y-axis units for documenting parameter units
 
   real ALLOCABLE_, dimension(NIMEM_,NJMEM_) :: &
-    bathyT           !< Ocean bottom depth at tracer points, in depth units [Z ~> m].
+    bathyT           !< Ocean bottom depth, referenced to Z_ref at tracer points. bathyT is in
+                     !! depth units and positive *below* Z_ref [Z ~> m].
+  real ALLOCABLE_, dimension(NIMEM_,NJMEM_) :: &
+    meanSL           !< Spatially varying time mean sea level, referenced to Z_ref at tracer points.
+                     !! meanSL is in height units and positive *above* Z_ref. It is used
+                     !! a) as the height where p = p_atm or zero;
+                     !! b) to calculate time mean thickness of the water column, where
+                     !!    mean thickness = max(meanSL + bathyT, 0.0).
+                     !! meanSL is 2D for the consideration of a domain with spatically varying mean
+                     !! height, e.g. the Great Lakes system [Z ~> m].
   real    :: Z_ref   !< A reference value for all geometric height fields, such as bathyT [Z ~> m].
 
   logical :: bathymetry_at_vel  !< If true, there are separate values for the
@@ -584,6 +593,7 @@ subroutine allocate_metrics(G)
   ALLOC_(G%IareaCv(isd:ied,JsdB:JedB)) ; G%IareaCv(:,:) = 0.0
 
   ALLOC_(G%bathyT(isd:ied, jsd:jed)) ; G%bathyT(:,:) = -G%Z_ref
+  ALLOC_(G%meanSL(isd:ied, jsd:jed)) ; G%meanSL(:,:) = G%Z_ref
   ALLOC_(G%CoriolisBu(IsdB:IedB, JsdB:JedB)) ; G%CoriolisBu(:,:) = 0.0
   ALLOC_(G%Coriolis2Bu(IsdB:IedB, JsdB:JedB)) ; G%Coriolis2Bu(:,:) = 0.0
   ALLOC_(G%dF_dx(isd:ied, jsd:jed)) ; G%dF_dx(:,:) = 0.0
@@ -631,9 +641,10 @@ subroutine MOM_grid_end(G)
 
   DEALLOC_(G%dx_Cv) ; DEALLOC_(G%dy_Cu)
 
-  DEALLOC_(G%bathyT)  ; DEALLOC_(G%CoriolisBu) ; DEALLOC_(G%Coriolis2Bu)
-  DEALLOC_(G%dF_dx)   ; DEALLOC_(G%dF_dy)
-  DEALLOC_(G%sin_rot) ; DEALLOC_(G%cos_rot)
+  DEALLOC_(G%bathyT)     ; DEALLOC_(G%meanSL)
+  DEALLOC_(G%CoriolisBu) ; DEALLOC_(G%Coriolis2Bu)
+  DEALLOC_(G%dF_dx)      ; DEALLOC_(G%dF_dy)
+  DEALLOC_(G%sin_rot)    ; DEALLOC_(G%cos_rot)
 
   DEALLOC_(G%porous_DminU) ; DEALLOC_(G%porous_DmaxU) ; DEALLOC_(G%porous_DavgU)
   DEALLOC_(G%porous_DminV) ; DEALLOC_(G%porous_DmaxV) ; DEALLOC_(G%porous_DavgV)
diff --git a/src/core/MOM_transcribe_grid.F90 b/src/core/MOM_transcribe_grid.F90
index d9ca19985f..5fd28164a2 100644
--- a/src/core/MOM_transcribe_grid.F90
+++ b/src/core/MOM_transcribe_grid.F90
@@ -56,6 +56,7 @@ subroutine copy_dyngrid_to_MOM_grid(dG, oG, US)
     oG%dyT(i,j) = dG%dyT(i+ido,j+jdo)
     oG%areaT(i,j) = dG%areaT(i+ido,j+jdo)
     oG%bathyT(i,j) = dG%bathyT(i+ido,j+jdo) - oG%Z_ref
+    oG%meanSL(i,j) = dG%meanSL(i+ido,j+jdo) + oG%Z_ref
 
     oG%dF_dx(i,j) = dG%dF_dx(i+ido,j+jdo)
     oG%dF_dy(i,j) = dG%dF_dy(i+ido,j+jdo)
@@ -145,6 +146,7 @@ subroutine copy_dyngrid_to_MOM_grid(dG, oG, US)
 ! Update the halos in case the dynamic grid has smaller halos than the ocean grid.
   call pass_var(oG%areaT, oG%Domain)
   call pass_var(oG%bathyT, oG%Domain)
+  call pass_var(oG%meanSL, oG%Domain)
   call pass_var(oG%geoLonT, oG%Domain)
   call pass_var(oG%geoLatT, oG%Domain)
   call pass_vector(oG%dxT, oG%dyT, oG%Domain, To_All+Scalar_Pair, AGRID)
@@ -217,6 +219,7 @@ subroutine copy_MOM_grid_to_dyngrid(oG, dG, US)
     dG%dyT(i,j) = oG%dyT(i+ido,j+jdo)
     dG%areaT(i,j) = oG%areaT(i+ido,j+jdo)
     dG%bathyT(i,j) = oG%bathyT(i+ido,j+jdo) + oG%Z_ref
+    dG%meanSL(i,j) = oG%meanSL(i+ido,j+jdo) - oG%Z_ref
 
     dG%dF_dx(i,j) = oG%dF_dx(i+ido,j+jdo)
     dG%dF_dy(i,j) = oG%dF_dy(i+ido,j+jdo)
@@ -307,6 +310,7 @@ subroutine copy_MOM_grid_to_dyngrid(oG, dG, US)
 ! Update the halos in case the dynamic grid has smaller halos than the ocean grid.
   call pass_var(dG%areaT, dG%Domain)
   call pass_var(dG%bathyT, dG%Domain)
+  call pass_var(dG%meanSL, dG%Domain)
   call pass_var(dG%geoLonT, dG%Domain)
   call pass_var(dG%geoLatT, dG%Domain)
   call pass_vector(dG%dxT, dG%dyT, dG%Domain, To_All+Scalar_Pair, AGRID)
diff --git a/src/diagnostics/MOM_wave_speed.F90 b/src/diagnostics/MOM_wave_speed.F90
index be5b29ac0e..09bb083c72 100644
--- a/src/diagnostics/MOM_wave_speed.F90
+++ b/src/diagnostics/MOM_wave_speed.F90
@@ -550,8 +550,9 @@ subroutine wave_speed(h, tv, G, GV, US, cg1, CS, halo_size, use_ebt_mode, mono_N
                 ! Determine whether N2 estimates should not be allowed to increase with depth.
                 if (l_mono_N2_column_fraction>0.) then
                   if (GV%Boussinesq .or. GV%semi_Boussinesq) then
-                    below_mono_N2_frac = (max(G%bathyT(i,j)+G%Z_ref, 0.0) - GV%H_to_Z*sum_hc < &
-                                          l_mono_N2_column_fraction*max(G%bathyT(i,j)+G%Z_ref, 0.0))
+                    below_mono_N2_frac = &
+                        (max(G%meanSL(i,j) + G%bathyT(i,j), 0.0) - GV%H_to_Z * sum_hc < &
+                         l_mono_N2_column_fraction * max(G%meanSL(i,j) + G%bathyT(i,j), 0.0))
                   else
                     below_mono_N2_frac = (htot(i) - sum_hc < l_mono_N2_column_fraction*htot(i))
                   endif
diff --git a/src/framework/MOM_dyn_horgrid.F90 b/src/framework/MOM_dyn_horgrid.F90
index 59edd72f9c..e5c6c1dab5 100644
--- a/src/framework/MOM_dyn_horgrid.F90
+++ b/src/framework/MOM_dyn_horgrid.F90
@@ -157,7 +157,16 @@ module MOM_dyn_horgrid
     y_ax_unit_short     !< A short description of the y-axis units for documenting parameter units
 
   real, allocatable, dimension(:,:) :: &
-    bathyT        !< Ocean bottom depth at tracer points, in depth units [Z ~> m].
+    bathyT        !< Ocean bottom depth, referenced to a zero reference height at tracer points.
+                  !! bathyT is in depth units and positive *below* the reference height [Z ~> m].
+  real, allocatable, dimension(:,:) :: &
+    meanSL        !< Spatially varying time mean sea level, referenced to a zero reference height
+                  !! at tracer points. meanSL is in height units and positive *above* zero. It is used
+                  !! a) as the height where p = p_atm or zero;
+                  !! b) to calculate time mean thickness of the water column, where
+                  !!    mean thickness = max(meanSL + bathyT, 0.0).
+                  !! meanSL is 2D for the consideration of a domain with spatically varying mean
+                  !! height, e.g. the Great Lakes system [Z ~> m].
 
   logical :: bathymetry_at_vel  !< If true, there are separate values for the
                   !! basin depths at velocity points.  Otherwise the effects of
@@ -290,8 +299,8 @@ subroutine create_dyn_horgrid(G, HI, bathymetry_at_vel)
   allocate(G%porous_DmaxV(isd:ied,JsdB:JedB), source=0.0)
   allocate(G%porous_DavgV(isd:ied,JsdB:JedB), source=0.0)
 
-
   allocate(G%bathyT(isd:ied, jsd:jed), source=0.0)
+  allocate(G%meanSL(isd:ied, jsd:jed), source=0.0)
   allocate(G%CoriolisBu(IsdB:IedB, JsdB:JedB), source=0.0)
   allocate(G%Coriolis2Bu(IsdB:IedB, JsdB:JedB), source=0.0)
   allocate(G%dF_dx(isd:ied, jsd:jed), source=0.0)
@@ -333,6 +342,7 @@ subroutine rotate_dyn_horgrid(G_in, G, US, turns)
   call rotate_array_pair(G_in%dxT, G_in%dyT, turns, G%dxT, G%dyT)
   call rotate_array(G_in%areaT, turns, G%areaT)
   call rotate_array(G_in%bathyT, turns, G%bathyT)
+  call rotate_array(G_in%meanSL, turns, G%meanSL)
 
   call rotate_array_pair(G_in%df_dx, G_in%df_dy, turns, G%df_dx, G%df_dy)
   call rotate_array(G_in%sin_rot, turns, G%sin_rot)
@@ -435,6 +445,7 @@ subroutine rescale_dyn_horgrid_bathymetry(G, m_in_new_units)
   rescale = 1.0 / m_in_new_units
   do j=jsd,jed ; do i=isd,ied
     G%bathyT(i,j) = rescale*G%bathyT(i,j)
+    G%meanSL(i,j) = rescale*G%meanSL(i,j)
   enddo ; enddo
   if (G%bathymetry_at_vel) then ; do j=jsd,jed ; do I=IsdB,IedB
     G%Dblock_u(I,j) = rescale*G%Dblock_u(I,j) ; G%Dopen_u(I,j) = rescale*G%Dopen_u(I,j)
@@ -519,7 +530,7 @@ subroutine destroy_dyn_horgrid(G)
   deallocate(G%areaT)  ; deallocate(G%IareaT)
   deallocate(G%areaBu) ; deallocate(G%IareaBu)
   deallocate(G%areaCu) ; deallocate(G%IareaCu)
-  deallocate(G%areaCv)  ; deallocate(G%IareaCv)
+  deallocate(G%areaCv) ; deallocate(G%IareaCv)
 
   deallocate(G%mask2dT)  ; deallocate(G%mask2dCu) ; deallocate(G%OBCmaskCu)
   deallocate(G%mask2dCv) ; deallocate(G%OBCmaskCv) ; deallocate(G%mask2dBu)
@@ -534,9 +545,10 @@ subroutine destroy_dyn_horgrid(G)
   deallocate(G%porous_DminU) ; deallocate(G%porous_DmaxU) ; deallocate(G%porous_DavgU)
   deallocate(G%porous_DminV) ; deallocate(G%porous_DmaxV) ; deallocate(G%porous_DavgV)
 
-  deallocate(G%bathyT)  ; deallocate(G%CoriolisBu) ; deallocate(G%Coriolis2Bu)
-  deallocate(G%dF_dx)   ; deallocate(G%dF_dy)
-  deallocate(G%sin_rot) ; deallocate(G%cos_rot)
+  deallocate(G%bathyT)     ; deallocate(G%meanSL)
+  deallocate(G%CoriolisBu) ; deallocate(G%Coriolis2Bu)
+  deallocate(G%dF_dx)      ; deallocate(G%dF_dy)
+  deallocate(G%sin_rot)    ; deallocate(G%cos_rot)
 
   if (allocated(G%Dblock_u)) deallocate(G%Dblock_u)
   if (allocated(G%Dopen_u)) deallocate(G%Dopen_u)
diff --git a/src/initialization/MOM_fixed_initialization.F90 b/src/initialization/MOM_fixed_initialization.F90
index 78559c72f2..78266f2749 100644
--- a/src/initialization/MOM_fixed_initialization.F90
+++ b/src/initialization/MOM_fixed_initialization.F90
@@ -25,6 +25,7 @@ module MOM_fixed_initialization
 use MOM_shared_initialization, only : read_face_length_list, set_velocity_depth_max, set_velocity_depth_min
 use MOM_shared_initialization, only : set_subgrid_topo_at_vel_from_file
 use MOM_shared_initialization, only : compute_global_grid_integrals
+use MOM_shared_initialization, only : set_meanSL_from_file
 use MOM_unit_scaling, only : unit_scale_type
 
 use user_initialization, only : user_initialize_topography
@@ -61,7 +62,8 @@ subroutine MOM_initialize_fixed(G, US, OBC, PF)
   ! Local variables
   character(len=200) :: inputdir   ! The directory where NetCDF input files are.
   character(len=200) :: config
-  logical            :: read_porous_file, OBC_projection_bug, open_corners, enable_bugs
+  logical            :: OBC_projection_bug, open_corners, enable_bugs
+  logical            :: read_porous_file, read_meanSL_file
   character(len=40)  :: mdl = "MOM_fixed_initialization" ! This module's name.
   integer :: I, J
   logical :: debug
@@ -79,10 +81,17 @@ subroutine MOM_initialize_fixed(G, US, OBC, PF)
   ! Set up the parameters of the physical domain (i.e. the grid), G
   call set_grid_metrics(G, PF, US)
 
+  ! Calculate time mean ocean total thickness
+  call get_param(PF, mdl, "READ_MEAN_SEA_LEVEL", read_meanSL_file, &
+                "If true, use a 2D map for time mean sea level, which is used to calculate "// &
+                "time mean ocean total thickness.", default=.False.)
+  if (read_meanSL_file) &
+    call set_meanSL_from_file(G%meanSL, G, PF, US)
+
   ! Set up the bottom depth, G%bathyT either analytically or from file
   ! This also sets G%max_depth based on the input parameter MAXIMUM_DEPTH,
   ! or, if absent, is diagnosed as G%max_depth = max( G%D(:,:) )
-  call MOM_initialize_topography(G%bathyT, G%max_depth, G, PF, US)
+  call MOM_initialize_topography(G%bathyT, G%max_depth, G, PF, US, meanSL=G%meanSL)
 
   ! To initialize masks, the bathymetry in halo regions must be filled in
   call pass_var(G%bathyT, G%Domain)
@@ -103,6 +112,12 @@ subroutine MOM_initialize_fixed(G, US, OBC, PF)
                  default=enable_bugs, do_not_log=.not.associated(OBC))
   open_corners = .not.OBC_projection_bug
 
+  if (associated(OBC) .and. OBC_projection_bug .and. read_meanSL_file) &
+    ! OBC_projection_bug modifies bathyT outside of the open boundaries, so meanSL would have to be
+    ! modified as well.
+    call MOM_error(FATAL, "MOM_initialize_fixed: To read mean sea level file, "//&
+                   "OBC_PROJECTION_BUG needs to be False.")
+
   ! This call sets masks that prohibit flow over any point interpreted as land
   if (associated(OBC)) then
     if (OBC_projection_bug) &
@@ -194,21 +209,29 @@ subroutine MOM_initialize_fixed(G, US, OBC, PF)
 end subroutine MOM_initialize_fixed
 
 !> MOM_initialize_topography makes the appropriate call to set up the bathymetry in units of [Z ~> m].
-subroutine MOM_initialize_topography(D, max_depth, G, PF, US)
+subroutine MOM_initialize_topography(D, max_depth, G, PF, US, meanSL)
   type(dyn_horgrid_type),           intent(in)  :: G  !< The dynamic horizontal grid type
   real, dimension(G%isd:G%ied,G%jsd:G%jed), &
                                     intent(out) :: D  !< Ocean bottom depth [Z ~> m]
   type(param_file_type),            intent(in)  :: PF !< Parameter file structure
-  real,                             intent(out) :: max_depth !< Maximum depth of model [Z ~> m]
+  real,                             intent(out) :: max_depth !< Maximum depth or geometric thickness,
+                                                             !! with meanSL present, of model [Z ~> m]
   type(unit_scale_type),            intent(in)  :: US !< A dimensional unit scaling type
+  real, dimension(G%isd:G%ied,G%jsd:G%jed), &
+                          optional, intent(in)  :: meanSL !< Mean sea level [Z ~> m]
 
   ! This subroutine makes the appropriate call to set up the bottom depth.
   ! This is a separate subroutine so that it can be made public and shared with
   ! the ice-sheet code or other components.
 
   ! Local variables
+  real :: max_depth_default = -1.e9 ! Default value of MAXIMUM_DEPTH parameter [m]
   character(len=40)  :: mdl = "MOM_initialize_topography" ! This subroutine's name.
   character(len=200) :: config
+  real, dimension(G%isd:G%ied, G%jsd:G%jed) :: D_meanSL ! depth (positive below meanSL) referenced
+                                        ! to meanSL. A temporary field used to diagnose maximum
+                                        ! static column thickness. D_meanSL = D + meanSL [Z ~> m].
+  integer :: i, j
 
   call get_param(PF, mdl, "TOPO_CONFIG", config, &
                  "This specifies how bathymetry is specified: \n"//&
@@ -238,7 +261,8 @@ subroutine MOM_initialize_topography(D, max_depth, G, PF, US)
                  " \t dense - Denmark Strait-like dense water formation and overflow.\n"//&
                  " \t USER - call a user modified routine.", &
                  fail_if_missing=.true.)
-  call get_param(PF, mdl, "MAXIMUM_DEPTH", max_depth, units="m", default=-1.e9, scale=US%m_to_Z, do_not_log=.true.)
+  call get_param(PF, mdl, "MAXIMUM_DEPTH", max_depth, units="m", default=max_depth_default, &
+                 scale=US%m_to_Z, do_not_log=.true.)
   select case ( trim(config) )
     case ("file");      call initialize_topography_from_file(D, G, PF, US)
     case ("flat");      call initialize_topography_named(D, G, PF, config, max_depth, US)
@@ -262,17 +286,27 @@ subroutine MOM_initialize_topography(D, max_depth, G, PF, US)
     case default ;      call MOM_error(FATAL,"MOM_initialize_topography: "// &
       "Unrecognized topography setup '"//trim(config)//"'")
   end select
-  if (max_depth>0.) then
+  if (max_depth /= max_depth_default * US%m_to_Z) then
     call log_param(PF, mdl, "MAXIMUM_DEPTH", max_depth, &
                    "The maximum depth of the ocean.", units="m", unscale=US%Z_to_m)
+    if (trim(config) /= "DOME") then
+      call limit_topography(D, G, PF, max_depth, US)
+    endif
   else
-    max_depth = diagnoseMaximumDepth(D,G)
+    if (present(meanSL)) then
+      D_meanSL(:,:) = 0.0
+      do j=G%jsc,G%jec ; do i=G%isc,G%iec ; D_meanSL(i,j) = D(i,j) + meanSL(i,j) ; enddo ; enddo
+      max_depth = diagnoseMaximumDepth(D_meanSL, G)
+    else
+      max_depth = diagnoseMaximumDepth(D, G)
+    endif
     call log_param(PF, mdl, "!MAXIMUM_DEPTH", max_depth, &
                    "The (diagnosed) maximum depth of the ocean.", &
                    units="m", unscale=US%Z_to_m, like_default=.true.)
-  endif
-  if (trim(config) /= "DOME") then
-    call limit_topography(D, G, PF, max_depth, US)
+    if (trim(config) /= "DOME") then
+      ! MAXIMUM_DEPTH is not set and topography does not need to be trimmed by its maximum depth.
+      call limit_topography(D, G, PF, -max_depth_default * US%m_to_Z, US)
+    endif
   endif
 
 end subroutine MOM_initialize_topography
diff --git a/src/initialization/MOM_shared_initialization.F90 b/src/initialization/MOM_shared_initialization.F90
index 6f8ed1ed8e..132dc069e7 100644
--- a/src/initialization/MOM_shared_initialization.F90
+++ b/src/initialization/MOM_shared_initialization.F90
@@ -30,6 +30,7 @@ module MOM_shared_initialization
 public read_face_length_list, set_velocity_depth_max, set_velocity_depth_min
 public set_subgrid_topo_at_vel_from_file
 public compute_global_grid_integrals, write_ocean_geometry_file
+public set_meanSL_from_file
 
 ! A note on unit descriptions in comments: MOM6 uses units that can be rescaled for dimensional
 ! consistency testing. These are noted in comments with units like Z, H, L, and T, along with
@@ -136,6 +137,42 @@ function diagnoseMaximumDepth(D, G)
   call max_across_PEs(diagnoseMaximumDepth)
 end function diagnoseMaximumDepth
 
+!> Read time mean ocean sea level from a file
+subroutine set_meanSL_from_file(meanSL, G, param_file, US)
+  type(dyn_horgrid_type),           intent(in)  :: G !< The dynamic horizontal grid type
+  real, dimension(G%isd:G%ied,G%jsd:G%jed), &
+                                    intent(out) :: meanSL !< Mean sea level referenced to a zero
+                                                          !! reference height at tracer points [Z ~> m].
+  type(param_file_type),            intent(in)  :: param_file !< Parameter file structure
+  type(unit_scale_type),            intent(in)  :: US !< A dimensional unit scaling type
+  ! Local variables
+  logical :: read_meanSL_file
+  character(len=200) :: filename, file, inputdir ! Strings for file/path
+  character(len=200) :: varname                  ! Variable name in file
+  character(len=40)  :: mdl = "set_meanSL_from_file" ! This subroutine's name.
+  integer :: i, j
+
+  call callTree_enter(trim(mdl)//"(), MOM_shared_initialization.F90")
+
+  call get_param(param_file, mdl, "INPUTDIR", inputdir, default=".")
+  inputdir = slasher(inputdir)
+  call get_param(param_file, mdl, "MEAN_SEA_LEVEL_FILE", file, &
+                 "The file from which the mean sea level is read.", &
+                 default="mean_sea_level.nc")
+  call get_param(param_file, mdl, "MEAN_SEA_LEVEL_VARNAME", varname, &
+                 "The name of the mean sea level variable in MEAN_SEA_LEVEL_FILE.", &
+                 default="meanSL")
+  filename = trim(inputdir)//trim(file)
+  call log_param(param_file, mdl, "INPUTDIR/TOPO_FILE", filename)
+
+  if (.not.file_exists(filename, G%Domain)) &
+    call MOM_error(FATAL, " "//mdl//": Unable to open "//trim(filename))
+
+  call MOM_read_data(filename, trim(varname), meanSL, G%Domain, scale=US%m_to_Z)
+  call pass_var(meanSL, G%Domain)
+
+  call callTree_leave(trim(mdl)//'()')
+end subroutine set_meanSL_from_file
 
 !> Read gridded depths from file
 subroutine initialize_topography_from_file(D, G, param_file, US)
diff --git a/src/parameterizations/lateral/MOM_MEKE.F90 b/src/parameterizations/lateral/MOM_MEKE.F90
index 7abc24fe2d..5f3ad1e19d 100644
--- a/src/parameterizations/lateral/MOM_MEKE.F90
+++ b/src/parameterizations/lateral/MOM_MEKE.F90
@@ -372,12 +372,12 @@ subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, US, CS, hu, h
       if (GV%Boussinesq) then
         !$OMP parallel do default(shared)
         do j=js-1,je+1 ; do i=is-1,ie+1
-          depth_tot(i,j) = max(G%bathyT(i,j) + G%Z_ref, 0.0) * GV%Z_to_H
+          depth_tot(i,j) = max(G%meanSL(i,j) + G%bathyT(i,j), 0.0) * GV%Z_to_H
         enddo ; enddo
       else
         !$OMP parallel do default(shared)
         do j=js-1,je+1 ; do i=is-1,ie+1
-          depth_tot(i,j) = max(G%bathyT(i,j) + G%Z_ref, 0.0) * CS%rho_fixed_total_depth * GV%RZ_to_H
+          depth_tot(i,j) = max(G%meanSL(i,j) + G%bathyT(i,j), 0.0) * CS%rho_fixed_total_depth * GV%RZ_to_H
         enddo ; enddo
       endif
     else
diff --git a/src/parameterizations/lateral/MOM_internal_tides.F90 b/src/parameterizations/lateral/MOM_internal_tides.F90
index 05c471cc17..73ff57adbc 100644
--- a/src/parameterizations/lateral/MOM_internal_tides.F90
+++ b/src/parameterizations/lateral/MOM_internal_tides.F90
@@ -3773,7 +3773,7 @@ subroutine internal_tides_init(Time, G, GV, US, param_file, diag, CS)
   do j=G%jsc,G%jec ; do i=G%isc,G%iec
     ! Restrict RMS topographic roughness to a fraction (10 percent by default) of the column depth.
     if (RMS_roughness_frac >= 0.0) then
-      h2(i,j) = max(min((RMS_roughness_frac * max(G%bathyT(i,j)+G%Z_ref, 0.0))**2, h2(i,j)), 0.0)
+      h2(i,j) = max(min((RMS_roughness_frac * max(G%meanSL(i,j) + G%bathyT(i,j), 0.0))**2, h2(i,j)), 0.0)
     else
       h2(i,j) = max(h2(i,j), 0.0)
     endif
diff --git a/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 b/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90
index baa4eb3d0c..2a1bea4b0b 100644
--- a/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90
+++ b/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90
@@ -206,6 +206,7 @@ subroutine calc_depth_function(G, CS)
   integer :: is, ie, js, je, Isq, Ieq, Jsq, Jeq
   integer :: i, j
   real    :: H0   ! The depth above which KHTH is linearly scaled away [Z ~> m]
+  real    :: h1, h2  ! Temporary total thicknesses [Z ~> m]
   real    :: expo ! exponent used in the depth dependent scaling [nondim]
   is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec
   Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB
@@ -224,13 +225,15 @@ subroutine calc_depth_function(G, CS)
   expo = CS%depth_scaled_khth_exp
 !$OMP do
   do j=js,je ; do I=is-1,Ieq
-    CS%Depth_fn_u(I,j) = (MIN(1.0, &
-      (0.5 * (max(G%bathyT(i,j) + G%Z_ref, 0.0) + max(G%bathyT(i+1,j) + G%Z_ref, 0.0))) / H0))**expo
+    h1 = max(G%meanSL(i,j) + G%bathyT(i,j), 0.0)
+    h2 = max(G%meanSL(i+1,j) + G%bathyT(i+1,j), 0.0)
+    CS%Depth_fn_u(I,j) = (MIN(1.0, (0.5 * (h1 + h2)) / H0))**expo
   enddo ; enddo
 !$OMP do
   do J=js-1,Jeq ; do i=is,ie
-    CS%Depth_fn_v(i,J) = (MIN(1.0, &
-      (0.5 * (max(G%bathyT(i,j) + G%Z_ref, 0.0) + max(G%bathyT(i,j+1) + G%Z_ref, 0.0))) / H0))**expo
+    h1 = max(G%meanSL(i,j) + G%bathyT(i,j), 0.0)
+    h2 = max(G%meanSL(i,j+1) + G%bathyT(i,j+1), 0.0)
+    CS%Depth_fn_v(i,J) = (MIN(1.0, (0.5 * (h1 + h2)) / H0))**expo
   enddo ; enddo
 
 end subroutine calc_depth_function
@@ -1194,6 +1197,7 @@ subroutine calc_slope_functions_using_just_e(h, G, GV, US, CS, e)
   ! real :: dz(SZI_(G),SZJ_(G),SZK_(GV)) ! The vertical distance across each layer [Z ~> m]
   real :: H_cutoff      ! Local estimate of a minimum thickness for masking [H ~> m or kg m-2]
   real :: dZ_cutoff     ! A minimum water column depth for masking [H ~> m or kg m-2]
+  real :: h1, h2        ! Temporary total thicknesses [Z ~> m]
   real :: h_neglect     ! A thickness that is so small it is usually lost
                         ! in roundoff and can be neglected [H ~> m or kg m-2].
   real :: S2            ! Interface slope squared [Z2 L-2 ~> nondim]
@@ -1303,9 +1307,10 @@ subroutine calc_slope_functions_using_just_e(h, G, GV, US, CS, e)
       enddo
     else
       do I=is-1,ie
-        if ( min(G%bathyT(i,j), G%bathyT(i+1,j)) + G%Z_ref > dZ_cutoff ) then
-          CS%SN_u(I,j) = G%OBCmaskCu(I,j) * sqrt( CS%SN_u(I,j) / &
-                                                  (max(G%bathyT(i,j), G%bathyT(i+1,j)) + G%Z_ref) )
+        h1 = max(G%meanSL(i,j) + G%bathyT(i,j), 0.0)
+        h2 = max(G%meanSL(i+1,j) + G%bathyT(i+1,j), 0.0)
+        if ( min(h1, h2) > dZ_cutoff ) then
+          CS%SN_u(I,j) = G%OBCmaskCu(I,j) * sqrt( CS%SN_u(I,j) / max(h1, h2) )
         else
           CS%SN_u(I,j) = 0.0
         endif
@@ -1328,9 +1333,10 @@ subroutine calc_slope_functions_using_just_e(h, G, GV, US, CS, e)
         ! There is a primordial horizontal indexing bug on the following line from the previous
         ! versions of the code.  This comment should be deleted by the end of 2024.
         ! if ( min(G%bathyT(i,j), G%bathyT(i+1,j)) + G%Z_ref > dZ_cutoff ) then
-        if ( min(G%bathyT(i,j), G%bathyT(i,j+1)) + G%Z_ref > dZ_cutoff ) then
-          CS%SN_v(i,J) = G%OBCmaskCv(i,J) * sqrt( CS%SN_v(i,J) / &
-                                                  (max(G%bathyT(i,j), G%bathyT(i,j+1)) + G%Z_ref) )
+        h1 = max(G%meanSL(i,j) + G%bathyT(i,j), 0.0)
+        h2 = max(G%meanSL(i,j+1) + G%bathyT(i,j+1), 0.0)
+        if ( min(h1, h2) > dZ_cutoff ) then
+          CS%SN_v(i,J) = G%OBCmaskCv(i,J) * sqrt( CS%SN_v(i,J) / max(h1, h2) )
         else
           CS%SN_v(i,J) = 0.0
         endif
diff --git a/src/parameterizations/vertical/MOM_internal_tide_input.F90 b/src/parameterizations/vertical/MOM_internal_tide_input.F90
index 2a7a3a2a7c..4bf0351039 100644
--- a/src/parameterizations/vertical/MOM_internal_tide_input.F90
+++ b/src/parameterizations/vertical/MOM_internal_tide_input.F90
@@ -551,13 +551,13 @@ subroutine int_tide_input_init(Time, G, GV, US, param_file, diag, CS, itide)
 
   do fr=1,num_freq ; do j=js,je ; do i=is,ie
     mask_itidal = 1.0
-    if (G%bathyT(i,j) + G%Z_ref < min_zbot_itides) mask_itidal = 0.0
+    if (G%meanSL(i,j) + G%bathyT(i,j) < min_zbot_itides) mask_itidal = 0.0
 
     CS%tideamp(i,j,fr) = CS%tideamp(i,j,fr) * mask_itidal * G%mask2dT(i,j)
 
     ! Restrict rms topo to a fraction (often 10 percent) of the column depth.
     if (max_frac_rough >= 0.0) &
-      itide%h2(i,j) = min((max_frac_rough * max(G%bathyT(i,j)+G%Z_ref, 0.0))**2, itide%h2(i,j))
+      itide%h2(i,j) = min((max_frac_rough * max(G%meanSL(i,j) + G%bathyT(i,j), 0.0))**2, itide%h2(i,j))
 
     ! Compute the fixed part of internal tidal forcing; units are [R Z4 H-1 T-2 ~> J m-2 or J m kg-1] here.
     CS%TKE_itidal_coef(i,j,fr) = 0.5*US%L_to_Z*kappa_h2_factor * GV%H_to_RZ * &
diff --git a/src/parameterizations/vertical/MOM_tidal_mixing.F90 b/src/parameterizations/vertical/MOM_tidal_mixing.F90
index 13b76a77a1..714a4efd50 100644
--- a/src/parameterizations/vertical/MOM_tidal_mixing.F90
+++ b/src/parameterizations/vertical/MOM_tidal_mixing.F90
@@ -507,16 +507,16 @@ logical function tidal_mixing_init(Time, G, GV, US, param_file, int_tide_CSp, di
                  units="nondim", default=0.1)
 
     do j=js,je ; do i=is,ie
-      if (G%bathyT(i,j)+G%Z_ref < CS%min_zbot_itides) CS%mask_itidal(i,j) = 0.0
+      if (max(G%meanSL(i,j) + G%bathyT(i,j), 0.0) < CS%min_zbot_itides) CS%mask_itidal(i,j) = 0.0
       CS%tideamp(i,j) = CS%tideamp(i,j) * CS%mask_itidal(i,j) * G%mask2dT(i,j)
 
       ! Restrict rms topo to a fraction (often 10 percent) of the column depth.
       if ((CS%tidal_answer_date < 20190101) .and. (max_frac_rough >= 0.0)) then
-        hamp = min(max_frac_rough*(G%bathyT(i,j)+G%Z_ref), sqrt(CS%h2(i,j)))
+        hamp = min(max_frac_rough * max(G%meanSL(i,j) + G%bathyT(i,j), 0.0), sqrt(CS%h2(i,j)))
         CS%h2(i,j) = hamp*hamp
       else
         if (max_frac_rough >= 0.0) &
-          CS%h2(i,j) = min((max_frac_rough * max(G%bathyT(i,j)+G%Z_ref, 0.0))**2, CS%h2(i,j))
+          CS%h2(i,j) = min((max_frac_rough * max(G%meanSL(i,j) + G%bathyT(i,j), 0.0))**2, CS%h2(i,j))
       endif
 
       utide = CS%tideamp(i,j)