diff --git a/src/core/MOM_PressureForce_FV.F90 b/src/core/MOM_PressureForce_FV.F90
index 64df200f31..5fb3ade634 100644
--- a/src/core/MOM_PressureForce_FV.F90
+++ b/src/core/MOM_PressureForce_FV.F90
@@ -319,7 +319,7 @@ subroutine PressureForce_FV_nonBouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_
       SSH(i,j) = (za(i,j) - alpha_ref*p(i,j,1)) * I_gEarth - G%Z_ref &
                  - max(-G%bathyT(i,j)-G%Z_ref, 0.0)
     enddo ; enddo
-    call calc_SAL(SSH, e_sal, G, CS%SAL_CSp)
+    call calc_SAL(SSH, e_sal, G, CS%SAL_CSp, tmp_scale=US%Z_to_m)
 
     if ((CS%tides_answer_date>20230630) .or. (.not.GV%semi_Boussinesq) .or. (.not.CS%tides)) then
       !$OMP parallel do default(shared)
@@ -587,7 +587,7 @@ subroutine PressureForce_FV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_atm
           SSH(i,j) = SSH(i,j) + h(i,j,k)*GV%H_to_Z
         enddo ; enddo
       enddo
-      call calc_SAL(SSH, e_sal, G, CS%SAL_CSp)
+      call calc_SAL(SSH, e_sal, G, CS%SAL_CSp, tmp_scale=US%Z_to_m)
       !$OMP parallel do default(shared)
       do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
         e(i,j,nz+1) = e(i,j,nz+1) - e_sal(i,j)
@@ -618,7 +618,7 @@ subroutine PressureForce_FV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_atm
           SSH(i,j) = SSH(i,j) + h(i,j,k)*GV%H_to_Z
         enddo ; enddo
       enddo
-      call calc_SAL(SSH, e_sal, G, CS%SAL_CSp)
+      call calc_SAL(SSH, e_sal, G, CS%SAL_CSp, tmp_scale=US%Z_to_m)
     else
       !$OMP parallel do default(shared)
       do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
diff --git a/src/core/MOM_PressureForce_Montgomery.F90 b/src/core/MOM_PressureForce_Montgomery.F90
index 3de713c801..6d982bc7e3 100644
--- a/src/core/MOM_PressureForce_Montgomery.F90
+++ b/src/core/MOM_PressureForce_Montgomery.F90
@@ -216,7 +216,7 @@ subroutine PressureForce_Mont_nonBouss(h, tv, PFu, PFv, G, GV, US, CS, p_atm, pb
       enddo ; enddo ; enddo
     endif
 
-    call calc_SAL(SSH, e_sal, G, CS%SAL_CSp)
+    call calc_SAL(SSH, e_sal, G, CS%SAL_CSp, tmp_scale=US%Z_to_m)
     !$OMP parallel do default(shared)
     do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
       geopot_bot(i,j) = geopot_bot(i,j) - GV%g_Earth*e_sal(i,j)
@@ -481,7 +481,7 @@ subroutine PressureForce_Mont_Bouss(h, tv, PFu, PFv, G, GV, US, CS, p_atm, pbce,
         SSH(i,j) = SSH(i,j) + h(i,j,k)*GV%H_to_Z
       enddo ; enddo
     enddo
-    call calc_SAL(SSH, e_sal, G, CS%SAL_CSp)
+    call calc_SAL(SSH, e_sal, G, CS%SAL_CSp, tmp_scale=US%Z_to_m)
     !$OMP parallel do default(shared)
     do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
       e(i,j,nz+1) = e(i,j,nz+1) - e_sal(i,j)
diff --git a/src/parameterizations/lateral/MOM_self_attr_load.F90 b/src/parameterizations/lateral/MOM_self_attr_load.F90
index 20d239eb53..7f7215c9d8 100644
--- a/src/parameterizations/lateral/MOM_self_attr_load.F90
+++ b/src/parameterizations/lateral/MOM_self_attr_load.F90
@@ -42,19 +42,21 @@ module MOM_self_attr_load
 !! be changed into bottom pressure anomaly in the future. Note that the SAL calculation applies to all motions
 !! across the spectrum. Tidal-specific methods that assume periodicity, i.e. iterative and read-in SAL, are
 !! stored in MOM_tidal_forcing module.
-subroutine calc_SAL(eta, eta_sal, G, CS)
+subroutine calc_SAL(eta, eta_sal, G, CS, tmp_scale)
   type(ocean_grid_type), intent(in)  :: G  !< The ocean's grid structure.
   real, dimension(SZI_(G),SZJ_(G)), intent(in)  :: eta     !< The sea surface height anomaly from
                                                            !! a time-mean geoid [Z ~> m].
   real, dimension(SZI_(G),SZJ_(G)), intent(out) :: eta_sal !< The sea surface height anomaly from
                                                            !! self-attraction and loading [Z ~> m].
   type(SAL_CS), intent(inout) :: CS !< The control structure returned by a previous call to SAL_init.
+  real,                   optional, intent(in)  :: tmp_scale !< A rescaling factor to temporarily convert eta
+                                                            !! to MKS units in reproducing sumes [m Z-1 ~> 1]
 
   ! Local variables
   integer :: n, m, l
   integer :: Isq, Ieq, Jsq, Jeq
   integer :: i, j
-  real :: eta_prop
+  real :: eta_prop ! The scalar constant of proportionality between eta and eta_sal [nondim]
 
   call cpu_clock_begin(id_clock_SAL)
 
@@ -69,7 +71,7 @@ subroutine calc_SAL(eta, eta_sal, G, CS)
 
   ! use the spherical harmonics method
   elseif (CS%use_sal_sht) then
-    call spherical_harmonics_forward(G, CS%sht, eta, CS%Snm_Re, CS%Snm_Im, CS%sal_sht_Nd)
+    call spherical_harmonics_forward(G, CS%sht, eta, CS%Snm_Re, CS%Snm_Im, CS%sal_sht_Nd, tmp_scale=tmp_scale)
 
     ! Multiply scaling factors to each mode
     do m = 0,CS%sal_sht_Nd
@@ -119,8 +121,8 @@ subroutine calc_love_scaling(nlm, rhoW, rhoE, Love_Scaling)
   real, dimension(:), intent(out) :: Love_Scaling !< Scaling factors for inverse SHT [nondim]
 
   ! Local variables
-  real, dimension(:), allocatable :: HDat, LDat, KDat ! Love numbers converted in CF reference frames
-  real :: H1, L1, K1 ! Temporary variables to store degree 1 Love numbers
+  real, dimension(:), allocatable :: HDat, LDat, KDat ! Love numbers converted in CF reference frames [nondim]
+  real :: H1, L1, K1 ! Temporary variables to store degree 1 Love numbers [nondim]
   integer :: n_tot ! Size of the stored Love numbers
   integer :: n, m, l
 
@@ -163,7 +165,7 @@ subroutine SAL_init(G, US, param_file, CS)
 
   logical :: calculate_sal
   logical :: tides, use_tidal_sal_file
-  real :: tide_sal_scalar_value
+  real :: tide_sal_scalar_value ! Scaling SAL factor [nondim]
 
   ! Read all relevant parameters and write them to the model log.
   call log_version(param_file, mdl, version, "")
diff --git a/src/parameterizations/lateral/MOM_spherical_harmonics.F90 b/src/parameterizations/lateral/MOM_spherical_harmonics.F90
index 2a72d26a20..26258e6b8e 100644
--- a/src/parameterizations/lateral/MOM_spherical_harmonics.F90
+++ b/src/parameterizations/lateral/MOM_spherical_harmonics.F90
@@ -42,7 +42,7 @@ module MOM_spherical_harmonics
 contains
 
 !> Calculates forward spherical harmonics transforms
-subroutine spherical_harmonics_forward(G, CS, var, Snm_Re, Snm_Im, Nd)
+subroutine spherical_harmonics_forward(G, CS, var, Snm_Re, Snm_Im, Nd, tmp_scale)
   type(ocean_grid_type), intent(in)    :: G            !< The ocean's grid structure.
   type(sht_CS),          intent(inout) :: CS           !< Control structure for SHT
   real, dimension(SZI_(G),SZJ_(G)), &
@@ -51,13 +51,20 @@ subroutine spherical_harmonics_forward(G, CS, var, Snm_Re, Snm_Im, Nd)
   real,                  intent(out)   :: Snm_Im(:)    !< SHT coefficients for the imaginary modes (sine) [A]
   integer,     optional, intent(in)    :: Nd           !< Maximum degree of the spherical harmonics
                                                        !! overriding ndegree in the CS [nondim]
+  real,        optional, intent(in)    :: tmp_scale    !< A temporary rescaling factor to convert
+                                                       !! var to MKS units during the reproducing
+                                                       !! sums [a A-1 ~> 1]
   ! local variables
-  integer :: Nmax ! Local copy of the maximum degree of the spherical harmonics [nondim]
-  integer :: Ltot ! Local copy of the number of spherical harmonics [nondim]
+  integer :: Nmax ! Local copy of the maximum degree of the spherical harmonics
+  integer :: Ltot ! Local copy of the number of spherical harmonics
   real, dimension(SZI_(G),SZJ_(G)) :: &
     pmn,   & ! Current associated Legendre polynomials of degree n and order m [nondim]
     pmnm1, & ! Associated Legendre polynomials of degree n-1 and order m [nondim]
     pmnm2    ! Associated Legendre polynomials of degree n-2 and order m [nondim]
+  real :: scale ! A rescaling factor to temporarily convert var to MKS units during the
+                ! reproducing sums [a A-1 ~> 1]
+  real :: I_scale ! The inverse of scale [A a-1 ~> 1]
+  real :: sum_tot ! The total of all components output by the reproducing sum in arbitrary units [a]
   integer :: i, j, k
   integer :: is, ie, js, je, isd, ied, jsd, jed
   integer :: m, n, l
@@ -81,12 +88,13 @@ subroutine spherical_harmonics_forward(G, CS, var, Snm_Re, Snm_Im, Nd)
   do l=1,Ltot ; Snm_Re(l) = 0.0; Snm_Im(l) = 0.0 ; enddo
 
   if (CS%reprod_sum) then
+    scale = 1.0 ; if (present(tmp_scale)) scale = tmp_scale
     do m=0,Nmax
       l = order2index(m, Nmax)
 
       do j=js,je ; do i=is,ie
-        CS%Snm_Re_raw(i,j,l) = var(i,j) * CS%Pmm(i,j,m+1) * CS%cos_lonT_wtd(i,j,m+1)
-        CS%Snm_Im_raw(i,j,l) = var(i,j) * CS%Pmm(i,j,m+1) * CS%sin_lonT_wtd(i,j,m+1)
+        CS%Snm_Re_raw(i,j,l) = (scale*var(i,j)) * CS%Pmm(i,j,m+1) * CS%cos_lonT_wtd(i,j,m+1)
+        CS%Snm_Im_raw(i,j,l) = (scale*var(i,j)) * CS%Pmm(i,j,m+1) * CS%sin_lonT_wtd(i,j,m+1)
         pmnm2(i,j) = 0.0
         pmnm1(i,j) = CS%Pmm(i,j,m+1)
       enddo ; enddo
@@ -94,8 +102,8 @@ subroutine spherical_harmonics_forward(G, CS, var, Snm_Re, Snm_Im, Nd)
       do n = m+1, Nmax ; do j=js,je ; do i=is,ie
         pmn(i,j) = &
           CS%a_recur(n+1,m+1) * CS%cos_clatT(i,j) * pmnm1(i,j) - CS%b_recur(n+1,m+1) * pmnm2(i,j)
-        CS%Snm_Re_raw(i,j,l+n-m) = var(i,j) * pmn(i,j) * CS%cos_lonT_wtd(i,j,m+1)
-        CS%Snm_Im_raw(i,j,l+n-m) = var(i,j) * pmn(i,j) * CS%sin_lonT_wtd(i,j,m+1)
+        CS%Snm_Re_raw(i,j,l+n-m) = (scale*var(i,j)) * pmn(i,j) * CS%cos_lonT_wtd(i,j,m+1)
+        CS%Snm_Im_raw(i,j,l+n-m) = (scale*var(i,j)) * pmn(i,j) * CS%sin_lonT_wtd(i,j,m+1)
         pmnm2(i,j) = pmnm1(i,j)
         pmnm1(i,j) = pmn(i,j)
       enddo ; enddo ; enddo
@@ -125,10 +133,15 @@ subroutine spherical_harmonics_forward(G, CS, var, Snm_Re, Snm_Im, Nd)
   if (id_clock_sht_global_sum>0) call cpu_clock_begin(id_clock_sht_global_sum)
 
   if (CS%reprod_sum) then
-    do l=1,Ltot
-      Snm_Re(l) = reproducing_sum(CS%Snm_Re_raw(:,:,l))
-      Snm_Im(l) = reproducing_sum(CS%Snm_Im_raw(:,:,l))
-    enddo
+    sum_tot = reproducing_sum(CS%Snm_Re_raw(:,:,1:Ltot), sums=Snm_Re(1:Ltot))
+    sum_tot = reproducing_sum(CS%Snm_Im_raw(:,:,1:Ltot), sums=Snm_Im(1:Ltot))
+    if (scale /= 1.0) then
+      I_scale = 1.0 / scale
+      do l=1,Ltot
+        Snm_Re(l) = I_scale * Snm_Re(l)
+        Snm_Im(l) = I_scale * Snm_Im(l)
+      enddo
+    endif
   else
     call sum_across_PEs(Snm_Re, Ltot)
     call sum_across_PEs(Snm_Im, Ltot)
@@ -240,8 +253,9 @@ subroutine spherical_harmonics_init(G, param_file, CS)
   allocate(CS%a_recur(CS%ndegree+1, CS%ndegree+1)); CS%a_recur(:,:) = 0.0
   allocate(CS%b_recur(CS%ndegree+1, CS%ndegree+1)); CS%b_recur(:,:) = 0.0
   do m=0,CS%ndegree ; do n=m+1,CS%ndegree
+    ! These expressione will give NaNs with 32-bit integers for n > 23170, but this is trapped elsewhere.
     CS%a_recur(n+1,m+1) = sqrt(real((2*n-1) * (2*n+1)) / real((n-m) * (n+m)))
-    CS%b_recur(n+1,m+1) = sqrt(real((2*n+1) * (n+m-1) * (n-m-1)) / real((n-m) * (n+m) * (2*n-3)))
+    CS%b_recur(n+1,m+1) = sqrt((real(2*n+1) * real((n+m-1) * (n-m-1))) / (real((n-m) * (n+m)) * real(2*n-3)))
   enddo ; enddo
 
   ! Calculate complex exponential factors
@@ -253,8 +267,8 @@ subroutine spherical_harmonics_init(G, param_file, CS)
     do j=js,je ; do i=is,ie
       CS%cos_lonT(i,j,m+1)     = cos(real(m) * (G%geolonT(i,j)*RADIAN))
       CS%sin_lonT(i,j,m+1)     = sin(real(m) * (G%geolonT(i,j)*RADIAN))
-      CS%cos_lonT_wtd(i,j,m+1) = CS%cos_lonT(i,j,m+1) * G%areaT(i,j) / G%Rad_Earth**2
-      CS%sin_lonT_wtd(i,j,m+1) = CS%sin_lonT(i,j,m+1) * G%areaT(i,j) / G%Rad_Earth**2
+      CS%cos_lonT_wtd(i,j,m+1) = CS%cos_lonT(i,j,m+1) * G%areaT(i,j) / G%Rad_Earth_L**2
+      CS%sin_lonT_wtd(i,j,m+1) = CS%sin_lonT(i,j,m+1) * G%areaT(i,j) / G%Rad_Earth_L**2
     enddo ; enddo
   enddo