diff --git a/src/parameterizations/lateral/MOM_self_attr_load.F90 b/src/parameterizations/lateral/MOM_self_attr_load.F90
index 7f6e35008c..b9b17df997 100644
--- a/src/parameterizations/lateral/MOM_self_attr_load.F90
+++ b/src/parameterizations/lateral/MOM_self_attr_load.F90
@@ -199,6 +199,8 @@ subroutine SAL_init(h, tv, G, GV, US, param_file, CS, restart_CS)
   character(len=40)  :: mdl = "MOM_self_attr_load" ! This module's name.
   integer :: lmax ! Total modes of the real spherical harmonics [nondim]
   real :: rhoE    ! The average density of Earth [R ~> kg m-3].
+  real :: rho_ocean_SAL ! The mean ocean density used for calculating the
+                        ! self-attraction and loading [R ~> kg m-3]
   character(len=20)  :: bpa_config ! String for reference bottom pressure config option
   real :: tmp(G%isd:G%ied, G%jsd:G%jed) ! Temporary field storing mass returned by find_col_mass
                                         ! [R Z ~> kg m-2]
@@ -305,8 +307,14 @@ subroutine SAL_init(h, tv, G, GV, US, param_file, CS, restart_CS)
                  default=0, do_not_log=(.not. CS%use_sal_sht))
   call get_param(param_file, mdl, "RHO_SOLID_EARTH", rhoE, &
                  "The mean solid earth density.  This is used for calculating the "// &
-                 "self-attraction and loading term.", units="kg m-3", &
-                 default=5517.0, scale=US%kg_m3_to_R, do_not_log=(.not. CS%use_sal_sht))
+                 "self-attraction and loading term.", &
+                 units="kg m-3", default=5517.0, scale=US%kg_m3_to_R, &
+                 do_not_log=(.not. CS%use_sal_sht))
+  call get_param(param_file, mdl, "RHO_OCEAN_SAL", rho_ocean_SAL, &
+                 "The mean ocean density used for calculating the self-attraction and loading.", &
+                 units="kg m-3", default=US%R_to_kg_m3*GV%Rho0, scale=US%kg_m3_to_R, &
+                 do_not_log=.not.(CS%use_sal_sht.or.CS%use_bpa))
+
 
   ! Set scaling coefficients for scalar approximation
   if (CS%use_sal_scalar .or. CS%use_tidal_sal_prev) then
@@ -316,7 +324,7 @@ subroutine SAL_init(h, tv, G, GV, US, param_file, CS, restart_CS)
       CS%eta_prop = sal_scalar_value
     endif
     if (CS%use_bpa) then
-      CS%linear_scaling = CS%eta_prop / (GV%Rho0 * GV%g_Earth)
+      CS%linear_scaling = CS%eta_prop / (rho_ocean_SAL * GV%g_Earth)
     else
       CS%linear_scaling = CS%eta_prop
     endif
@@ -331,7 +339,7 @@ subroutine SAL_init(h, tv, G, GV, US, param_file, CS, restart_CS)
     allocate(CS%Snm_Im(lmax), source=0.0)
 
     allocate(CS%Love_scaling(lmax), source=0.0)
-    call calc_love_scaling(GV%Rho0, rhoE, GV%g_Earth, CS)
+    call calc_love_scaling(rho_ocean_SAL, rhoE, GV%g_Earth, CS)
 
     allocate(CS%sht)
     call spherical_harmonics_init(G, param_file, CS%sht)