diff --git a/components/mpas-ocean/bld/build-namelist b/components/mpas-ocean/bld/build-namelist
index 19418582863e..9354c7afb733 100755
--- a/components/mpas-ocean/bld/build-namelist
+++ b/components/mpas-ocean/bld/build-namelist
@@ -718,7 +718,7 @@ add_default($nl, 'config_enable_shortwave_energy_fixer');
 ###########################################
 
 add_default($nl, 'config_use_self_attraction_loading');
-add_default($nl, 'config_self_attraction_loading_smoothing_width');
+add_default($nl, 'config_self_attraction_loading_depth_cutoff');
 add_default($nl, 'config_mpas_to_grid_weights_file');
 add_default($nl, 'config_grid_to_mpas_weights_file');
 add_default($nl, 'config_self_attraction_loading_compute_interval');
@@ -818,8 +818,11 @@ add_default($nl, 'config_loglaw_bottom_drag_min');
 add_default($nl, 'config_loglaw_bottom_drag_max');
 add_default($nl, 'config_explicit_bottom_drag_coeff');
 add_default($nl, 'config_use_topographic_wave_drag');
+add_default($nl, 'config_topographic_wave_drag_scheme');
 add_default($nl, 'config_topographic_wave_drag_coeff');
 add_default($nl, 'config_thickness_drag_type');
+add_default($nl, 'config_topographic_wave_drag_cutoff_depth');
+add_default($nl, 'config_topographic_wave_drag_cutoff_width');
 
 ####################################
 # Namelist group: Rayleigh_damping #
diff --git a/components/mpas-ocean/bld/build-namelist-section b/components/mpas-ocean/bld/build-namelist-section
index 0c834f05cec1..7f28313ca01f 100644
--- a/components/mpas-ocean/bld/build-namelist-section
+++ b/components/mpas-ocean/bld/build-namelist-section
@@ -243,7 +243,7 @@ add_default($nl, 'config_enable_shortwave_energy_fixer');
 ###########################################
 
 add_default($nl, 'config_use_self_attraction_loading');
-add_default($nl, 'config_self_attraction_loading_smoothing_width');
+add_default($nl, 'config_self_attraction_loading_depth_cutoff');
 add_default($nl, 'config_mpas_to_grid_weights_file');
 add_default($nl, 'config_grid_to_mpas_weights_file');
 add_default($nl, 'config_self_attraction_loading_compute_interval');
@@ -335,8 +335,11 @@ add_default($nl, 'config_loglaw_bottom_drag_min');
 add_default($nl, 'config_loglaw_bottom_drag_max');
 add_default($nl, 'config_explicit_bottom_drag_coeff');
 add_default($nl, 'config_use_topographic_wave_drag');
+add_default($nl, 'config_topographic_wave_drag_scheme');
 add_default($nl, 'config_topographic_wave_drag_coeff');
 add_default($nl, 'config_thickness_drag_type');
+add_default($nl, 'config_topographic_wave_drag_cutoff_depth');
+add_default($nl, 'config_topographic_wave_drag_cutoff_width');
 
 ####################################
 # Namelist group: Rayleigh_damping #
diff --git a/components/mpas-ocean/bld/namelist_files/namelist_defaults_mpaso.xml b/components/mpas-ocean/bld/namelist_files/namelist_defaults_mpaso.xml
index 7ff1b6744d87..2932f91c287f 100644
--- a/components/mpas-ocean/bld/namelist_files/namelist_defaults_mpaso.xml
+++ b/components/mpas-ocean/bld/namelist_files/namelist_defaults_mpaso.xml
@@ -295,7 +295,7 @@
 
 <!-- self_attraction_loading -->
 <config_use_self_attraction_loading>.false.</config_use_self_attraction_loading>
-<config_self_attraction_loading_smoothing_width>1.0</config_self_attraction_loading_smoothing_width>
+<config_self_attraction_loading_depth_cutoff>10.0</config_self_attraction_loading_depth_cutoff>
 <config_mpas_to_grid_weights_file>'mpas_to_grid.nc'</config_mpas_to_grid_weights_file>
 <config_grid_to_mpas_weights_file>'grid_to_mpas.nc'</config_grid_to_mpas_weights_file>
 <config_self_attraction_loading_compute_interval>'0000-00-00_00:30:00'</config_self_attraction_loading_compute_interval>
@@ -399,8 +399,11 @@
 <config_loglaw_bottom_drag_max>1.0e-1</config_loglaw_bottom_drag_max>
 <config_explicit_bottom_drag_coeff>1.0e-3</config_explicit_bottom_drag_coeff>
 <config_use_topographic_wave_drag>.false.</config_use_topographic_wave_drag>
+<config_topographic_wave_drag_scheme>'ZAE'</config_topographic_wave_drag_scheme>
 <config_topographic_wave_drag_coeff>5.0e-4</config_topographic_wave_drag_coeff>
 <config_thickness_drag_type>'centered'</config_thickness_drag_type>
+<config_topographic_wave_drag_cutoff_depth>500</config_topographic_wave_drag_cutoff_depth>
+<config_topographic_wave_drag_cutoff_width>10</config_topographic_wave_drag_cutoff_width>
 
 <!-- Rayleigh_damping -->
 <config_Rayleigh_damping_coeff>0.0</config_Rayleigh_damping_coeff>
diff --git a/components/mpas-ocean/bld/namelist_files/namelist_definition_mpaso.xml b/components/mpas-ocean/bld/namelist_files/namelist_definition_mpaso.xml
index 196d4abf730c..ae872caf67d9 100644
--- a/components/mpas-ocean/bld/namelist_files/namelist_definition_mpaso.xml
+++ b/components/mpas-ocean/bld/namelist_files/namelist_definition_mpaso.xml
@@ -1229,9 +1229,9 @@ Valid values: .true. or .false.
 Default: Defined in namelist_defaults.xml
 </entry>
 
-<entry id="config_self_attraction_loading_smoothing_width" type="real"
+<entry id="config_self_attraction_loading_depth_cutoff" type="real"
 	category="self_attraction_loading" group="self_attraction_loading">
-Defines region over which ssh is smoothed to zero at coasts for SAL calculation.
+Defines depths over which ssh is smoothed to zero for SAL calculation.
 
 Valid values: Any positive real number.
 Default: Defined in namelist_defaults.xml
@@ -1780,11 +1780,19 @@ Valid values: .true. or .false.
 Default: Defined in namelist_defaults.xml
 </entry>
 
+<entry id="config_topographic_wave_drag_scheme" type="char*1024"
+	category="bottom_drag" group="bottom_drag">
+Which of the following three wave drag schemes to use: JSL: Jayne and St. Laurent; ZAE: Zaron and Egbert; or LGF: Local generation formula
+
+Valid values: JSL or ZAE or LGF
+Default: Defined in namelist_defaults.xml
+</entry>
+
 <entry id="config_topographic_wave_drag_coeff" type="real"
 	category="bottom_drag" group="bottom_drag">
 Dimensionless topographic wave drag coefficient, $c_{topo}$.
 
-Valid values: any positive real, typically 5.0e-4
+Valid values: any positive real
 Default: Defined in namelist_defaults.xml
 </entry>
 
@@ -1796,6 +1804,22 @@ Valid values: 'harmonic-mean', 'centered'
 Default: Defined in namelist_defaults.xml
 </entry>
 
+<entry id="config_topographic_wave_drag_cutoff_depth" type="real"
+	category="bottom_drag" group="bottom_drag">
+Specifies minimum depth at which topographic wave drag is applied
+
+Valid values: any positive real
+Default: Defined in namelist_defaults.xml
+</entry>
+
+<entry id="config_topographic_wave_drag_cutoff_width" type="real"
+	category="bottom_drag" group="bottom_drag">
+Width over which the topographic wave drag is progressively turned off
+
+Valid values: any positive real
+Default: Defined in namelist_defaults.xml
+</entry>
+
 
 <!-- Rayleigh_damping -->
 
diff --git a/components/mpas-ocean/src/Registry.xml b/components/mpas-ocean/src/Registry.xml
index bb399577eee3..af10e1d1548c 100644
--- a/components/mpas-ocean/src/Registry.xml
+++ b/components/mpas-ocean/src/Registry.xml
@@ -829,8 +829,8 @@
 					description="Controls if self-attraction and loading is applied to ssh"
 					possible_values=".true. or .false."
 		/>
-		<nml_option name="config_self_attraction_loading_smoothing_width" type="real" default_value="1.0" units="km"
-					description="Defines region over which ssh is smoothed to zero at coasts for SAL calculation."
+		<nml_option name="config_self_attraction_loading_depth_cutoff" type="real" default_value="10.0" units="m"
+					description="Defines depths over which ssh is smoothed to zero for SAL calculation."
 					possible_values="Any positive real number."
 		/>
 		<nml_option name="config_mpas_to_grid_weights_file" type="character" default_value="mpas_to_grid.nc"
@@ -1108,14 +1108,26 @@
 					description="If true, topographic wave drag is used on the momentum equation."
 					possible_values=".true. or .false."
 		/>
-		<nml_option name="config_topographic_wave_drag_coeff" type="real" default_value="5.0e-4"
+		<nml_option name="config_topographic_wave_drag_scheme" type="character" default_value="ZAE" units="unitless"
+					description="Which of the following three wave drag schemes to use: JSL: Jayne and St. Laurent; ZAE: Zaron and Egbert; or LGF: Local generation formula"
+					possible_values="JSL or ZAE or LGF"
+		/>
+		<nml_option name="config_topographic_wave_drag_coeff" type="real" default_value="0.3" units="unitless"
 					description="Dimensionless topographic wave drag coefficient, $c_{topo}$."
-					possible_values="any positive real, typically 5.0e-4"
+					possible_values="any positive real"
 		/>
 		<nml_option name="config_thickness_drag_type" type="character" default_value="centered"
 					description="The type of layerThickness averaging to use on the drag term. The standard MPAS-O approach is 'centered'."
 					possible_values="'harmonic-mean', 'centered'"
 		/>
+		<nml_option name="config_topographic_wave_drag_cutoff_depth" type="real" default_value="500" units="m"
+					description="Specifies minimum depth at which topographic wave drag is applied"
+					possible_values="any positive real"
+		/>
+		<nml_option name="config_topographic_wave_drag_cutoff_width" type="real" default_value="10" units="m"
+					description="Width over which the topographic wave drag is progressively turned off"
+					possible_values="any positive real"
+		/>
 	</nml_record>
 	<nml_record name="Rayleigh_damping" mode="forward">
 		<nml_option name="config_Rayleigh_damping_coeff" type="real" default_value="1.0e-4" units="s^-1"
@@ -1962,7 +1974,7 @@
             <var name="gotmEpsbTopOfCell"/>
             <var name="gotmDissTopOfCell"/>
             <var name="gotmLengthTopOfCell"/>
-            <var name="ssh_sal"/>
+            <var name="pgf_sal"/>
 		</stream>
 
 		<stream name="output"
@@ -2262,7 +2274,13 @@
 				filename_template="topographic_wave_drag.nc"
 				runtime_format="single_file"
 				mode="forward;analysis">
-			<var name="topographic_wave_drag"/>
+			<var name="topo_rinv"/>
+			<var name="topo_buoyancy_N1V"/>
+			<var name="topo_buoyancy_N1B"/>
+			<var name="bathy_stddev"/>
+			<var name="bed_slope_edges"/>
+			<var name="lonGradEdge"/>
+			<var name="latGradEdge"/>
 		</stream>
 
 	</streams>
@@ -3233,10 +3251,34 @@
 		<var name="surfaceFluxAttenuationCoefficientRunoff" type="real" dimensions="nCells Time" units="m"
 			description="The spatially-dependent length scale of exponential decay of river runoff. Fluxes are multiplied by $e^{z/\gamma}$, where this coefficient is $\gamma$."
 		/>
-		<var name="topographic_wave_drag" type="real" dimensions="nEdges" units="1"
+		<var name="topo_rinv" type="real" dimensions="nEdges" units="unitless"
 			description="wave drag coefficient or 1/(rinv) where rinv is the e-folding time used in HyCOM"
 			packages="topographicWaveDragPKG"
 		/>
+		<var name="topo_buoyancy_N1V" type="real" dimensions="nEdges" units="unitless"
+			description="Vertically averaged buoyancy frequencies from WOA 2013 data"
+			packages="topographicWaveDragPKG"
+		/>
+		<var name="topo_buoyancy_N1B" type="real" dimensions="nEdges" units="unitless"
+			description="Bottom buoyancy frequencies from WOA 2013 data"
+			packages="topographicWaveDragPKG"
+		/>
+		<var name="bathy_stddev" type="real" dimensions="nEdges" units="unitless"
+			description="Standard deviation of the bathymetry within a given cell"
+			packages="topographicWaveDragPKG"
+		/>
+		<var name="bed_slope_edges" type="real" dimensions="nEdges" units="unitless"
+			description="The integral gradient average per cell"
+			packages="topographicWaveDragPKG"
+		/>
+		<var name="lonGradEdge" type="real" dimensions="nEdges" units="unitless"
+			description="Meridional gradient of bathymetry"
+			packages="topographicWaveDragPKG"
+		/>
+		<var name="latGradEdge" type="real" dimensions="nEdges" units="unitless"
+			description="Zonal gradient of bathymetry"
+			packages="topographicWaveDragPKG"
+		/>
 		<!-- diagnostic fields for land-ice fluxes -->
 		<var name="landIceFrictionVelocity" type="real" dimensions="nCells Time" units="m s^-1"
 			description="The friction velocity $u_*$ under land ice"
@@ -3419,12 +3461,12 @@
 			description="A vector used in the split-implicit barotropic mode solver"
 			packages="semiImplicitTimePKG"
 		/>
-		<var name="ssh_sal" type="real" dimensions="nCells Time" units="m"
-			description="sea surface height perturbation from self-attraction and loading"
+		<var name="pgf_sal" type="real" dimensions="nCells Time" units="m"
+			description="Barotropic pressure perturbation due to self-attraction and loading"
 			packages="tidalPotentialForcingPKG"
 		/>
-		<var name="ssh_sal_grad" type="real" dimensions="nEdges Time" units=""
-			description="gradient of sea surface height perturbation from self-attraction and loading"
+		<var name="pgf_sal_grad" type="real" dimensions="nEdges Time" units=""
+			description="Gradient of barotropic pressure perturbation due to self-attraction and loading"
 			packages="tidalPotentialForcingPKG"
 		/>
 		<!-- FIELD split_explicit AB2 time stepping -->
@@ -3436,6 +3478,22 @@
                         description="Coriolis term at the previous time step used in the split-explicit AB2 time stepping"
 			packages="splitAB2TimeIntegrator"
 		/>
+		<var name="temp_twd" type="real" dimensions="nEdges Time" units=""
+			description="tendency perturbation due to topographic wave drag"
+			packages="tidalPotentialForcingPKG"
+		/>
+		<var name="normalCoeffTWD" type="real" dimensions="nEdges" units=""
+			description="Coefficient multiplied by the normal velocity in the topographic wave drag"
+			packages="tidalPotentialForcingPKG"
+		/>
+		<var name="tangentialCoeffTWD" type="real" dimensions="nEdges" units=""
+			description="Coefficient multiplied by the tangential velocity in the topographic wave drag"
+			packages="tidalPotentialForcingPKG"
+		/>
+		<var name="topographicWaveDrag" type="real" dimensions="nEdges" units=""
+			description="topographic wave drag multiplicative factor"
+			packages="tidalPotentialForcingPKG"
+		/>
 	</var_struct>
 	<var_struct name="shortwave" time_levs="1">
 		<!-- **********************************************************************
diff --git a/components/mpas-ocean/src/analysis_members/Registry_harmonic_analysis.xml b/components/mpas-ocean/src/analysis_members/Registry_harmonic_analysis.xml
index 9337cc1e6e23..a43f3c0536d3 100644
--- a/components/mpas-ocean/src/analysis_members/Registry_harmonic_analysis.xml
+++ b/components/mpas-ocean/src/analysis_members/Registry_harmonic_analysis.xml
@@ -80,6 +80,9 @@
 		<var name="analysisConstituentNodalPhase" type="real" dimensions="maxTidalConstituents Time" units="rad/s"
 			description="Frequency of each constituent"
 		/>
+		<var name="sshInit" type="real" dimensions="nCells" units="m"
+			description="Copy of the initial ssh distribution. Used to compute analysis with respect to ssh anomalies, ensuring cells with initial ssh offsets are handled correctly"
+		/>
 		<var name="leastSquaresLHSMatrix" type="real" dimensions="maxTidalConstituentsX2 maxTidalConstituentsX2 Time" units="rad/s"
 			description="Frequency of each constituent"
 		/>
@@ -154,6 +157,7 @@
 		                 output_interval="stream:restart:output_interval"
 		                 immutable="false"
 		                 mode="forward;analysis">
+		    <var name="sshInit"/>
 			<var name="leastSquaresLHSMatrix"/>
 			<var name="leastSquaresRHSVector"/>
 		</stream>
diff --git a/components/mpas-ocean/src/analysis_members/mpas_ocn_harmonic_analysis.F b/components/mpas-ocean/src/analysis_members/mpas_ocn_harmonic_analysis.F
index 6f53e885f784..ee1495a3e59a 100644
--- a/components/mpas-ocean/src/analysis_members/mpas_ocn_harmonic_analysis.F
+++ b/components/mpas-ocean/src/analysis_members/mpas_ocn_harmonic_analysis.F
@@ -117,11 +117,12 @@ subroutine ocn_init_harmonic_analysis(domain, err)!{{{
 
       type (dm_info) :: dminfo
       type (block_type), pointer :: block
-      type (mpas_pool_type), pointer :: harmonicAnalysisAMPool 
+      type (mpas_pool_type), pointer :: statePool, harmonicAnalysisAMPool 
 
-      integer, pointer :: nAnalysisConstituents 
+      integer, pointer :: nAnalysisConstituents, nCellsSolve
       real (kind=RKIND), pointer :: harmonicAnalysisStart
       real (kind=RKIND), pointer :: harmonicAnalysisEnd
+      real (kind=RKIND), dimension(:),   pointer :: ssh, sshInit
       real (kind=RKIND), dimension(:,:), pointer :: leastSquaresLHSMatrix
       real (kind=RKIND), dimension(:,:), pointer :: leastSquaresRHSVector
       real (kind=RKIND), dimension(:,:), pointer :: decomposedConstituentAmplitude 
@@ -135,7 +136,7 @@ subroutine ocn_init_harmonic_analysis(domain, err)!{{{
       integer, dimension(:), allocatable :: tidalConstituentType
 
       type (MPAS_Time_Type) :: refTime
-      integer :: iCon
+      integer :: iCon, iCell
 
       err = 0
 
@@ -145,8 +146,13 @@ subroutine ocn_init_harmonic_analysis(domain, err)!{{{
 
       block => domain % blocklist
       do while (associated(block))
+         call mpas_pool_get_subpool(block % structs, 'state', statePool)
          call mpas_pool_get_subpool(block % structs, 'harmonicAnalysisAM', harmonicAnalysisAMPool)
 
+         call mpas_pool_get_dimension(block % dimensions, 'nCellsSolve', nCellsSolve)
+         call mpas_pool_get_array(statePool, 'ssh', ssh, 1)
+
+         call mpas_pool_get_array(harmonicAnalysisAMPool, 'sshInit', sshInit)
          call mpas_pool_get_array(harmonicAnalysisAMPool, 'nAnalysisConstituents', nAnalysisConstituents)
          call mpas_pool_get_array(harmonicAnalysisAMPool, 'leastSquaresLHSMatrix', leastSquaresLHSMatrix)
          call mpas_pool_get_array(harmonicAnalysisAMPool, 'leastSquaresRHSVector', leastSquaresRHSVector)
@@ -214,11 +220,14 @@ subroutine ocn_init_harmonic_analysis(domain, err)!{{{
                                         tidalConstituentType, &
                                         err)
 
-          if (.not. config_do_restart) then
-            leastSquaresRHSVector = 0.0_RKIND
-            leastSquaresLHSMatrix = 0.0_RKIND
-          end if
+         if (.not. config_do_restart) then
+           leastSquaresRHSVector = 0.0_RKIND
+           leastSquaresLHSMatrix = 0.0_RKIND
+         end if
            
+         do iCell = 1,nCellSSolve
+           sshInit(iCell) = ssh(iCell)
+         end do
 
          do iCon = 1,nAnalysisConstituents  
            call mpas_log_write('Constituent '//constituentList(iCon)%constituent)
@@ -291,7 +300,8 @@ subroutine ocn_compute_harmonic_analysis(domain, timeLevel, err)!{{{
 
       integer, pointer :: nCellsSolve      
       integer, pointer :: nAnalysisConstituents
-      real (kind=RKIND), dimension(:),   pointer :: ssh
+      real (kind=RKIND), dimension(:),   pointer :: ssh, sshInit
+      real (kind=RKIND), dimension(:), allocatable :: sshDiff
       real (kind=RKIND), dimension(:,:), pointer :: leastSquaresLHSMatrix
       real (kind=RKIND), dimension(:,:), pointer :: leastSquaresRHSVector
       real (kind=RKIND), dimension(:),   pointer :: analysisConstituentFrequency 
@@ -334,6 +344,7 @@ subroutine ocn_compute_harmonic_analysis(domain, timeLevel, err)!{{{
 
          call mpas_pool_get_dimension(block % dimensions, 'nCellsSolve', nCellsSolve)
          call mpas_pool_get_array(statePool, 'ssh', ssh, timeLevel)
+         call mpas_pool_get_array(harmonicAnalysisAMPool, 'sshInit', sshInit)
          call mpas_pool_get_array(harmonicAnalysisAMPool, 'nAnalysisConstituents', nAnalysisConstituents)
          call mpas_pool_get_array(harmonicAnalysisAMPool, 'leastSquaresLHSMatrix', leastSquaresLHSMatrix)
          call mpas_pool_get_array(harmonicAnalysisAMPool, 'leastSquaresRHSVector', leastSquaresRHSVector)
@@ -355,6 +366,12 @@ subroutine ocn_compute_harmonic_analysis(domain, timeLevel, err)!{{{
 
          ! update if within harmonic analysis period 
          if (daysSinceStartOfSim .le. harmonicAnalysisEnd) then
+           allocate(sshDiff(nCellsSolve))
+
+           do iCell = 1,nCellsSolve
+             sshDiff(iCell) = ssh(iCell) - sshInit(iCell)
+           end do
+
            CALL update_least_squares_LHS_matrix(nAnalysisConstituents, &
                                                 time, &
                                                 analysisConstituentFrequency, &
@@ -364,8 +381,9 @@ subroutine ocn_compute_harmonic_analysis(domain, timeLevel, err)!{{{
                                                 time, &
                                                 nCellsSolve, &
                                                 analysisConstituentFrequency, &
-                                                ssh, &
+                                                sshDiff, &
                                                 leastSquaresRHSVector)
+           deallocate(sshDiff)
            call mpas_log_write('harmonicAnalysisAM update')
          end if
 
diff --git a/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_si.F b/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_si.F
index 743c250ed306..9529fbdc6ea1 100644
--- a/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_si.F
+++ b/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_si.F
@@ -155,7 +155,7 @@ module ocn_time_integration_si
       SIvec_zh1_field,SIvec_wh0_field,SIvec_rh1_field   
                          ! field pointer for halo update
 
-   integer :: ssh_sal_on
+   real(kind=RKIND) :: pgf_sal_on
 
 #ifdef USE_GPU_AWARE_MPI
    ! SI - halo exch related variables ----------------------------------
@@ -715,7 +715,7 @@ subroutine ocn_time_integrator_si(domain, dt)!{{{
       endif
 
       if (config_use_tidal_potential_forcing) then
-         !$acc enter data copyin(ssh_sal)
+         !$acc enter data copyin(pgf_sal)
          !$acc enter data create(sshSubcycleCurWithTides,tidalPotentialEta)
       endif
 #endif
@@ -1160,7 +1160,7 @@ subroutine ocn_time_integrator_si(domain, dt)!{{{
  
                !$acc parallel loop gang  &
                !$acc    present(sshSubcycleCurWithTides,sshSubcycleCur, &
-               !$acc       tidalPotentialEta,ssh_sal,sshSubcycleCur)
+               !$acc       tidalPotentialEta,pgf_sal,sshSubcycleCur)
 #else
                !$omp parallel
                !$omp do schedule(runtime)
@@ -1169,8 +1169,8 @@ subroutine ocn_time_integrator_si(domain, dt)!{{{
                   sshSubcycleCurWithTides(iCell) = &
                            sshSubcycleCur(iCell) - &
                         tidalPotentialEta(iCell) - &
-                     ssh_sal_on * ssh_sal(iCell) - &
-                     (1 - ssh_sal_on) * self_attraction_and_loading_beta* &
+                     pgf_sal_on * pgf_sal(iCell) - &
+                     (1.0_RKIND - pgf_sal_on) * self_attraction_and_loading_beta* &
                            sshSubcycleCur(iCell)
                end do
 #ifndef MPAS_OPENACC
@@ -2801,7 +2801,7 @@ subroutine ocn_time_integrator_si(domain, dt)!{{{
 
 #ifdef MPAS_OPENACC
       if (config_use_tidal_potential_forcing) then
-         !$acc exit data delete(ssh_sal,tidalPotentialEta)
+         !$acc exit data delete(pgf_sal,tidalPotentialEta)
          !$acc exit data copyout(sshSubcycleCurWithTides)
       endif
 #endif
@@ -3501,10 +3501,10 @@ subroutine ocn_time_integration_si_init(domain, dt)!{{{
       !--------------------------------------------------------------------
 
       if (config_use_self_attraction_loading) then
-         ! set ssh_sal_on to 1
-         ssh_sal_on = 1
+         ! set pgf_sal_on to 1
+         pgf_sal_on = 1.0_RKIND
       else
-         ssh_sal_on = 0
+         pgf_sal_on = 0.0_RKIND
       endif
 
       self_attraction_and_loading_beta = config_self_attraction_and_loading_beta
diff --git a/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_split.F b/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_split.F
index 4f24b544033c..653d88198131 100644
--- a/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_split.F
+++ b/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_split.F
@@ -91,12 +91,12 @@ module ocn_time_integration_split
       numTSIterations, &! number of outer timestep iterations
       nBtrSubcycles     ! number of barotropic subcycles
 
-   integer :: ssh_sal_on
    integer :: &
       thickEdgeFluxChoice ! choice of thickness flux type
    integer, parameter :: &
       thickEdgeFluxCenter   = 1, &! use mean thickness of cell neighbors
       thickEdgeFluxUpwind   = 2   ! use upwind cell thickness at edge
+   real (kind=RKIND) :: pgf_sal_on
 
    real (kind=RKIND) :: &
       useVelocityCorrection,&! mask for velocity correction
@@ -1020,8 +1020,8 @@ subroutine ocn_time_integrator_split(domain, dt)!{{{
                         sshSubcycleCurWithTides(iCell) = &
                                  sshSubcycleCur(iCell) - &
                               tidalPotentialEta(iCell) - &
-                           ssh_sal_on * ssh_sal(iCell) - &
-                           (1 - ssh_sal_on) * self_attraction_and_loading_beta* &
+                           pgf_sal_on * pgf_sal(iCell) - &
+                           (1.0_RKIND - pgf_sal_on) * self_attraction_and_loading_beta* &
                                  sshSubcycleCur(iCell)
                      end do
                      !$omp end do
@@ -1268,15 +1268,15 @@ subroutine ocn_time_integrator_split(domain, dt)!{{{
                         sshSubcycleCurWithTides(iCell) = &
                                  sshSubcycleCur(iCell) - &
                               tidalPotentialEta(iCell) - &
-                           ssh_sal_on * ssh_sal(iCell) - &
-                           (1 - ssh_sal_on) * self_attraction_and_loading_beta* &
+                           pgf_sal_on * pgf_sal(iCell) - &
+                           (1.0_RKIND - pgf_sal_on) * self_attraction_and_loading_beta* &
                                  sshSubcycleCur(iCell)
 
                         sshSubcycleNewWithTides(iCell) = &
                                  sshSubcycleNew(iCell) - &
                               tidalPotentialEta(iCell) - &
-                           ssh_sal_on * ssh_sal(iCell) - &
-                           (1 - ssh_sal_on) * self_attraction_and_loading_beta* &
+                           pgf_sal_on * pgf_sal(iCell) - &
+                           (1.0_RKIND - pgf_sal_on) * self_attraction_and_loading_beta* &
                                  sshSubcycleNew(iCell)
                      end do
                      !$omp end do
@@ -3071,10 +3071,10 @@ subroutine ocn_time_integration_split_init(domain)!{{{
 
       !-----------------------------------------------------------------
       if (config_use_self_attraction_loading) then
-         ! set ssh_sal_on to 1
-         ssh_sal_on = 1
+         ! set pgf_sal_on to 1
+         pgf_sal_on = 1.0_RKIND
       else
-         ssh_sal_on = 0
+         pgf_sal_on = 0.0_RKIND
       endif
 
       self_attraction_and_loading_beta = config_self_attraction_and_loading_beta
diff --git a/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_split_ab2.F b/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_split_ab2.F
index a719bd24e135..c6040aeba1d6 100644
--- a/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_split_ab2.F
+++ b/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_split_ab2.F
@@ -92,7 +92,7 @@ module ocn_time_integration_split_ab2
       numTSIterations, &! number of outer timestep iterations
       nBtrSubcycles     ! number of barotropic subcycles
 
-   integer :: ssh_sal_on
+   real (kind=RKIND) :: pgf_sal_on
    integer :: &
       thickEdgeFluxChoice ! choice of thickness flux type
    integer, parameter :: &
@@ -1076,8 +1076,8 @@ subroutine ocn_time_integrator_split_ab2(domain, dt)!{{{
                      sshSubcycleCurWithTides(iCell) = &
                               sshSubcycleCur(iCell) - &
                            tidalPotentialEta(iCell) - &
-                        ssh_sal_on * ssh_sal(iCell) - &
-                        (1 - ssh_sal_on) * self_attraction_and_loading_beta* &
+                        pgf_sal_on * pgf_sal(iCell) - &
+                        (1.0_RKIND - pgf_sal_on) * self_attraction_and_loading_beta* &
                               sshSubcycleCur(iCell)
                   end do
                   !$omp end do
@@ -1294,15 +1294,15 @@ subroutine ocn_time_integrator_split_ab2(domain, dt)!{{{
                      sshSubcycleCurWithTides(iCell) = &
                               sshSubcycleCur(iCell) - &
                            tidalPotentialEta(iCell) - &
-                        ssh_sal_on * ssh_sal(iCell) - &
-                        (1 - ssh_sal_on) * self_attraction_and_loading_beta* &
+                        pgf_sal_on * pgf_sal(iCell) - &
+                        (1.0_RKIND - pgf_sal_on) * self_attraction_and_loading_beta* &
                               sshSubcycleCur(iCell)
 
                      sshSubcycleNewWithTides(iCell) = &
                               sshSubcycleNew(iCell) - &
                            tidalPotentialEta(iCell) - &
-                        ssh_sal_on * ssh_sal(iCell) - &
-                        (1 - ssh_sal_on) * self_attraction_and_loading_beta* &
+                        pgf_sal_on * pgf_sal(iCell) - &
+                        (1.0_RKIND - pgf_sal_on) * self_attraction_and_loading_beta* &
                               sshSubcycleNew(iCell)
                   end do
                   !$omp end do
@@ -3215,10 +3215,10 @@ subroutine ocn_time_integration_split_ab2_init(domain)!{{{
 
       !-----------------------------------------------------------------
       if (config_use_self_attraction_loading) then
-         ! set ssh_sal_on to 1
-         ssh_sal_on = 1
+         ! set pgf_sal_on to 1
+         pgf_sal_on = 1.0_RKIND
       else
-         ssh_sal_on = 0
+         pgf_sal_on = 0.0_RKIND
       endif
 
       self_attraction_and_loading_beta = config_self_attraction_and_loading_beta
diff --git a/components/mpas-ocean/src/mode_init/mpas_ocn_init_global_ocean.F b/components/mpas-ocean/src/mode_init/mpas_ocn_init_global_ocean.F
index 926400345e65..1f75bd1c515d 100644
--- a/components/mpas-ocean/src/mode_init/mpas_ocn_init_global_ocean.F
+++ b/components/mpas-ocean/src/mode_init/mpas_ocn_init_global_ocean.F
@@ -893,7 +893,8 @@ subroutine ocn_init_setup_global_ocean_create_model_topo(domain, iErr)!{{{
              end if
              if (isOcean) then
                 ! Enforce minimum depth
-                bottomDepth(iCell) = max(bottomDepthObserved(iCell), refBottomDepth(minimum_levels))
+                bottomDepth(iCell) = max(bottomDepthObserved(iCell), config_global_ocean_minimum_depth)
+                bottomDepth(iCell) = min(bottomDepth(iCell), refBottomDepth(minimum_levels))
 
                 maxLevelCell(iCell) = -1
                 do k = 1, nVertLevels
diff --git a/components/mpas-ocean/src/mode_init/mpas_ocn_init_test_sht.F b/components/mpas-ocean/src/mode_init/mpas_ocn_init_test_sht.F
index 33eaf9af0127..ae3450695435 100644
--- a/components/mpas-ocean/src/mode_init/mpas_ocn_init_test_sht.F
+++ b/components/mpas-ocean/src/mode_init/mpas_ocn_init_test_sht.F
@@ -89,7 +89,7 @@ subroutine ocn_init_setup_test_sht(domain, iErr)!{{{
       integer, pointer :: nCellsSolve
       real (kind=RKIND), dimension(:), pointer :: lonCell, latCell
       real (kind=RKIND), dimension(:), pointer :: xCell, yCell, zCell
-      real (kind=RKIND), dimension(:), pointer :: ssh, ssh_sal, ssh_sal_grad
+      real (kind=RKIND), dimension(:), pointer :: ssh, surfacePressure, pgf_sal, pgf_sal_grad
       real (kind=RKIND), dimension(:), pointer :: areaCell, dcEdge
       real (kind=RKIND), pointer :: config_test_sht_cosine_bell_lat_center, &
                                     config_test_sht_cosine_bell_lon_center, &
@@ -188,9 +188,11 @@ subroutine ocn_init_setup_test_sht(domain, iErr)!{{{
         call mpas_pool_get_array(meshPool, 'cellsOnEdge', cellsOnEdge)
 
         call mpas_pool_get_array(statePool, 'ssh', ssh)
-        call mpas_pool_get_array(diagnosticsPool, 'ssh_sal', ssh_sal)
-        call mpas_pool_get_array(diagnosticsPool, 'ssh_sal_grad', ssh_sal_grad)
+        call mpas_pool_get_array(diagnosticsPool, 'surfacePressure', surfacePressure)
+        call mpas_pool_get_array(diagnosticsPool, 'pgf_sal', pgf_sal)
+        call mpas_pool_get_array(diagnosticsPool, 'pgf_sal_grad', pgf_sal_grad)
 
+        surfacePressure(:) = 0.0_RKIND
         
         ! Setup function to approximate
         if (config_test_sht_function_option == 1) then
@@ -283,22 +285,23 @@ subroutine ocn_init_setup_test_sht(domain, iErr)!{{{
         call ocn_vel_self_attraction_loading_init(domain,err)
         do i = 1,config_test_sht_n_iterations
           call mpas_log_write('Iteration = $i',intArgs=(/i/))
-          call ocn_compute_self_attraction_loading(domain, forcingPool, domain % dminfo, ssh, err)
+          call ocn_compute_self_attraction_loading(domain, forcingPool, domain % dminfo, &
+                                                   ssh, surfacePressure, err)
         enddo
 
         ! Compute gradient
         do iEdge = 1,nEdgesArray(1)
           cell1 = cellsOnEdge(1,iEdge)
           cell2 = cellsOnEdge(2,iEdge)
-          ssh_sal_grad(iEdge) = (ssh_sal(cell2) - ssh_sal(cell1))/dcEdge(iEdge) 
+          pgf_sal_grad(iEdge) = (pgf_sal(cell2) - pgf_sal(cell1))/dcEdge(iEdge) 
         enddo
 
         ! Compute RMS Error
         error_local_area = 0.0_RKIND
         error_local = 0.0_RKIND
         do iCell = 1,nCellsSolve
-          error_local_area = error_local_area + areaCell(iCell)*(ssh(iCell)-ssh_sal(iCell))**2
-          error_local = error_local + (ssh(iCell)-ssh_sal(iCell))**2
+          error_local_area = error_local_area + areaCell(iCell)*(ssh(iCell)-pgf_sal(iCell))**2
+          error_local = error_local + (ssh(iCell)-pgf_sal(iCell))**2
         enddo 
         call mpas_dmpar_sum_real(domain % dminfo,error_local_area,error_global_area)
         call mpas_dmpar_sum_real(domain % dminfo,error_local,error_global)
diff --git a/components/mpas-ocean/src/shared/mpas_ocn_diagnostics_variables.F b/components/mpas-ocean/src/shared/mpas_ocn_diagnostics_variables.F
index 567c3efc183c..e4392ffa7737 100644
--- a/components/mpas-ocean/src/shared/mpas_ocn_diagnostics_variables.F
+++ b/components/mpas-ocean/src/shared/mpas_ocn_diagnostics_variables.F
@@ -72,7 +72,11 @@ module ocn_diagnostics_variables
 
    real (kind=RKIND), dimension(:), pointer :: gradSSH
    real (kind=RKIND), dimension(:), pointer :: pressureAdjustedSSH
-   real (kind=RKIND), dimension(:), pointer :: ssh_sal
+   real (kind=RKIND), dimension(:), pointer :: pgf_sal
+   real (kind=RKIND), dimension(:), pointer :: temp_twd
+   real (kind=RKIND), dimension(:), pointer :: normalCoeffTWD 
+   real (kind=RKIND), dimension(:), pointer :: tangentialCoeffTWD 
+   real (kind=RKIND), dimension(:), pointer :: topographicWaveDrag 
 
    real (kind=RKIND), dimension(:,:), pointer :: tracersSurfaceValue
    real (kind=RKIND), dimension(:,:), pointer :: tracersSurfaceLayerValue
@@ -85,7 +89,13 @@ module ocn_diagnostics_variables
    real (kind=RKIND), dimension(:), pointer :: landIceFrictionVelocity
    real (kind=RKIND), dimension(:), pointer :: topDrag
    real (kind=RKIND), dimension(:), pointer :: topDragMag
-   real (kind=RKIND), dimension(:), pointer :: topographic_wave_drag
+   real (kind=RKIND), dimension(:), pointer :: topo_rinv
+   real (kind=RKIND), dimension(:), pointer :: topo_buoyancy_N1V
+   real (kind=RKIND), dimension(:), pointer :: topo_buoyancy_N1B
+   real (kind=RKIND), dimension(:), pointer :: bathy_stddev
+   real (kind=RKIND), dimension(:), pointer :: bed_slope_edges
+   real (kind=RKIND), dimension(:), pointer :: lonGradEdge
+   real (kind=RKIND), dimension(:), pointer :: latGradEdge
    real (kind=RKIND), dimension(:,:), pointer :: landIceBoundaryLayerTracers
    real (kind=RKIND), dimension(:,:), pointer :: landIceTracerTransferVelocities
 
@@ -254,7 +264,13 @@ subroutine ocn_diagnostics_variables_init(domain, jenkinsOn, hollandJenkinsOn, e
       call mpas_pool_get_subpool(domain % blocklist % structs, 'diagnostics', diagnosticsPool)
 
       if (config_use_topographic_wave_drag) then
-         call mpas_pool_get_array(diagnosticsPool, 'topographic_wave_drag', topographic_wave_drag)
+         call mpas_pool_get_array(diagnosticsPool, 'topo_rinv', topo_rinv)
+         call mpas_pool_get_array(diagnosticsPool, 'topo_buoyancy_N1V', topo_buoyancy_N1V)
+         call mpas_pool_get_array(diagnosticsPool, 'topo_buoyancy_N1B', topo_buoyancy_N1B)
+         call mpas_pool_get_array(diagnosticsPool, 'bathy_stddev', bathy_stddev)
+         call mpas_pool_get_array(diagnosticsPool, 'bed_slope_edges', bed_slope_edges)
+         call mpas_pool_get_array(diagnosticsPool, 'lonGradEdge', lonGradEdge)
+         call mpas_pool_get_array(diagnosticsPool, 'latGradEdge', latGradEdge)
       end if
 
       if ( trim(config_ocean_run_mode) == 'forward' ) then
@@ -368,17 +384,48 @@ subroutine ocn_diagnostics_variables_init(domain, jenkinsOn, hollandJenkinsOn, e
       end if
 
       if (config_use_topographic_wave_drag) then
+         call mpas_pool_get_array(diagnosticsPool, 'topo_rinv', topo_rinv)
          call mpas_pool_get_array(diagnosticsPool, &
-                  'topographic_wave_drag', &
-                   topographic_wave_drag)
+                    'topo_buoyancy_N1V', &
+                     topo_buoyancy_N1V)
+         call mpas_pool_get_array(diagnosticsPool, &
+                    'topo_buoyancy_N1B', &
+                     topo_buoyancy_N1B)
+         call mpas_pool_get_array(diagnosticsPool, &
+                    'bathy_stddev', &
+                     bathy_stddev)
+         call mpas_pool_get_array(diagnosticsPool, &
+                    'bed_slope_edges', &
+                     bed_slope_edges)
+         call mpas_pool_get_array(diagnosticsPool, &
+                  'lonGradEdge', &
+                   lonGradEdge)
+         call mpas_pool_get_array(diagnosticsPool, &
+                  'latGradEdge', &
+                   latGradEdge)
       end if
 
       if (config_use_tidal_potential_forcing) then
          call mpas_pool_get_array(diagnosticsPool, &
-                  'ssh_sal', &
-                   ssh_sal)
+                  'pgf_sal', &
+                   pgf_sal)
+      end if
+      if (config_use_topographic_wave_drag) then
+         call mpas_pool_get_array(diagnosticsPool, &
+                  'temp_twd', &
+                   temp_twd)
+         call mpas_pool_get_array(diagnosticsPool, &
+                  'normalCoeffTWD', &
+                   normalCoeffTWD)
+         call mpas_pool_get_array(diagnosticsPool, &
+                  'tangentialCoeffTWD', &
+                   tangentialCoeffTWD)
+      end if
+      if (config_use_topographic_wave_drag) then
+         call mpas_pool_get_array(diagnosticsPool, &
+                  'topographicWaveDrag', &
+                   topographicWaveDrag)
       end if
-
       if ( trim(config_ocean_run_mode) == 'init' .or. &
                 config_AM_debugDiagnostics_enable ) then
          call mpas_pool_get_array(diagnosticsPool, 'rx1MaxCell', &
@@ -688,10 +735,20 @@ subroutine ocn_diagnostics_variables_init(domain, jenkinsOn, hollandJenkinsOn, e
          !$acc                   landIceFrictionVelocity)
       end if
       if (config_use_topographic_wave_drag) then
-         !$acc enter data create(topographic_wave_drag)
+         !$acc enter data create(topo_rinv)
+         !$acc enter data create(topo_buoyancy_N1V)
+         !$acc enter data create(topo_buoyancy_N1B)
+         !$acc enter data create(bathy_stddev)
+         !$acc enter data create(bed_slope_edges)
+         !$acc enter data create(lonGradEdge)
+         !$acc enter data create(latGradEdge)
+         !$acc enter data create(temp_twd)
+         !$acc enter data create(normalCoeffTWD)
+         !$acc enter data create(tangentialCoeffTWD)
+         !$acc enter data create(topographicWaveDrag)
       end if
       if (config_use_self_attraction_loading) then
-         !$acc enter data copyin(ssh_sal)
+         !$acc enter data copyin(pgf_sal)
       end if
       if ( trim(config_ocean_run_mode) == 'init' ) then
          !$acc enter data create(                          &
@@ -936,9 +993,13 @@ subroutine ocn_diagnostics_variables_destroy(err) !{{{
       end if
       if (config_use_topographic_wave_drag) then
          !$acc exit data delete(topographic_wave_drag)
+         !$acc exit data delete(temp_twd)
+         !$acc exit data delete(normalCoeffTWD)
+         !$acc exit data delete(tangentialCoeffTWD)
+         !$acc exit data delete(topographicWaveDrag)
       end if
       if (config_use_self_attraction_loading) then
-         !$acc exit data delete(ssh_sal)
+         !$acc exit data delete(pgf_sal)
       end if
       if ( trim(config_ocean_run_mode) == 'init' ) then
          !$acc exit data delete(                           &
@@ -1275,7 +1336,13 @@ subroutine ocn_diagnostics_variables_destroy(err) !{{{
          nullify(salinitySurfaceRestoringTendency)
       end if
       if ( config_use_self_attraction_loading ) then
-         nullify(ssh_sal)
+         nullify(pgf_sal)
+      end if
+      if ( config_use_topographic_wave_drag) then
+         nullify(temp_twd)
+         nullify(normalCoeffTWD)
+         nullify(tangentialCoeffTWD)
+         nullify(topographicWaveDrag)
       end if
       if ( trim(config_time_integrator) == 'split_implicit' ) then
          nullify(SIvec_r0,              &
diff --git a/components/mpas-ocean/src/shared/mpas_ocn_tendency.F b/components/mpas-ocean/src/shared/mpas_ocn_tendency.F
index 91f182db9993..a16e52f14499 100644
--- a/components/mpas-ocean/src/shared/mpas_ocn_tendency.F
+++ b/components/mpas-ocean/src/shared/mpas_ocn_tendency.F
@@ -426,12 +426,13 @@ subroutine ocn_tend_vel(domain, tendPool, statePool, forcingPool, &
                                thermExpCoeff, salineContractCoeff, &
                                tendVel, err)
 
-      ! Compute ssh additions due to self-attraction and loading
+      ! Compute pgf additions due to self-attraction and loading
       if(mpas_is_alarm_ringing(domain % clock, 'salComputeAlarm', ierr=err)) then 
 #ifdef MPAS_DEBUG
           call mpas_log_write('       Computing SAL')
 #endif
-          call ocn_compute_self_attraction_loading(domain, forcingPool, dminfo, ssh, err)
+          call ocn_compute_self_attraction_loading(domain, forcingPool, dminfo, &
+                                                   ssh, surfacePressure, err)
           call mpas_reset_clock_alarm(domain % clock, 'salComputeAlarm', ierr=err)
       endif
 
diff --git a/components/mpas-ocean/src/shared/mpas_ocn_vel_forcing_topographic_wave_drag.F b/components/mpas-ocean/src/shared/mpas_ocn_vel_forcing_topographic_wave_drag.F
index 6cb7362af143..8f641ca4f31c 100644
--- a/components/mpas-ocean/src/shared/mpas_ocn_vel_forcing_topographic_wave_drag.F
+++ b/components/mpas-ocean/src/shared/mpas_ocn_vel_forcing_topographic_wave_drag.F
@@ -10,7 +10,7 @@
 !  ocn_vel_forcing_topographic_wave_drag
 !
 !> \brief MPAS ocean topographic wave drag
-!> \author Nairita Pal
+!> \author Nairita Pal, Kristin Barton
 !> \date   Oct 2020
 !> \details
 !>  This module contains the routine for computing
@@ -55,8 +55,9 @@ module ocn_vel_forcing_topographic_wave_drag
    !
    !--------------------------------------------------------------------
 
-   logical :: topographicWaveDragOn
+   logical :: topographicWaveDragOn, tensorScheme 
    real (kind=RKIND) :: topographicWaveDragCoeff
+!   real (kind=RKIND), dimension(:), allocatable :: topographicWaveDrag
 
 !***********************************************************************
 
@@ -67,7 +68,7 @@ module ocn_vel_forcing_topographic_wave_drag
 !  routine ocn_vel_forcing_topographic_wave_drag_tend
 !
 !> \brief   Computes tendency term from topographic wave drag
-!> \author  Nairita Pal
+!> \author  Nairita Pal, Kristin Barton
 !> \date    15 Oct 2020
 !> \details
 !>  This routine computes the topographic wave drag tendency for momentum
@@ -117,19 +118,32 @@ subroutine ocn_vel_forcing_topographic_wave_drag_tend(normalVelocity, & !{{{
 
       call mpas_timer_start('vel topographic wave drag')
 
-      !$omp do schedule(runtime) private(k)
-      do iEdge = 1, nEdgesOwned
-        k =  maxLevelEdgeTop(iEdge)
-        ! topographic wave drag term:
-        ! du/dt = ... 1/rinv * u
-        ! rinv is the e-folding time use in HyCOM. Here
-        ! topographic_wave_drag = 1/rinv
-        if (k>0) then
-          tend(k,iEdge) = tend(k,iEdge) - topographicWaveDragCoeff &
-             * topographic_wave_drag(iEdge) * normalVelocity(k,iEdge)
-        endif
-      enddo
-      !$omp end do
+      ! JSL: C = pi std(H)^2 N_B / L
+      ! ZAE: C = gamma H (grad(H))^2 N_B bar(N) / (8 * pi^2 * omega)
+      !          N(z) = N_0 exp(z/1300)
+      ! Drag: rho_0 C u
+      
+      if ( .not. tensorScheme) then
+          !$omp do schedule(runtime) private(k)
+          do iEdge = 1, nEdgesOwned
+            k =  maxLevelEdgeTop(iEdge)
+            if (k>0) then
+              temp_twd(iEdge) = topographicWaveDrag(iEdge) * normalVelocity(k,iEdge)
+              tend(k,iEdge) = tend(k,iEdge) - temp_twd(iEdge)
+            endif
+          enddo
+          !$omp end do
+      else
+          !$omp do schedule(runtime) private(k)
+          do iEdge = 1, nEdgesOwned
+            k =  maxLevelEdgeTop(iEdge)
+            if (k>0) then
+              temp_twd(iEdge) = topographicWaveDrag(iEdge) * (normalVelocity(k,iEdge) * normalCoeffTWD(iEdge) + tangentialVelocity(k,iEdge) * tangentialCoeffTWD(iEdge))
+              tend(k,iEdge) = tend(k,iEdge) - temp_twd(iEdge)
+            endif
+          enddo
+          !$omp end do
+      endif
 
       call mpas_timer_stop('vel topographic wave drag')
 
@@ -142,7 +156,7 @@ end subroutine ocn_vel_forcing_topographic_wave_drag_tend!}}}
 !  routine ocn_vel_forcing_topographic_wave_drag_init
 !
 !> \brief   Initializes ocean topographic wave drag forcing
-!> \author  Nairita Pal
+!> \author  Nairita Pal, Kristin Barton
 !> \date    Oct 2020
 !> \details
 !>  This routine initializes quantities related to topographic wave  drag
@@ -162,13 +176,113 @@ subroutine ocn_vel_forcing_topographic_wave_drag_init(err)!{{{
 
       integer, intent(out) :: err !< Output: error flag
 
+      real (kind=RKIND) :: topoDragCutoffDepth, topoDragCutoffWidth
+      real (kind=RKIND) :: omegaM2 = 1.405189e-4_RKIND
+      real (kind=RKIND) :: kappa
+      real (kind=RKIND) :: N0, L, Nb, Nbar, gam
+      real (kind=RKIND) :: topoDragTurnOff
+
+      integer :: iEdge, k, kmax
+
       err = 0
 
+
       topographicWaveDragCoeff = 0.0_RKIND
 
       if (config_use_topographic_wave_drag) then
           topographicWaveDragOn = .true.
           topographicWaveDragCoeff = config_topographic_wave_drag_coeff
+          topoDragCutoffDepth = config_topographic_wave_drag_cutoff_depth
+          topoDragCutoffWidth = config_topographic_wave_drag_cutoff_width
+          if (config_topographic_wave_drag_scheme.EQ."JSL") then
+              call mpas_log_write("")
+              call mpas_log_write(" Topographic wave drag scheme is: Jayne and St. Laurent")
+              call mpas_log_write("")
+              tensorScheme = .false.
+              kappa = pii/10000.0_RKIND
+
+              do iEdge = 1, nEdgesAll
+                  kmax =  maxLevelEdgeTop(iEdge)
+                  topographicWaveDrag(iEdge) = 0.0_RKIND
+                  if (kmax > 0) then
+                      k =  maxLevelEdgeTop(iEdge)
+                      topographicWaveDrag(iEdge) = topographicWaveDragCoeff*kappa &
+                                                 * bathy_stddev(iEdge)**2 &
+                                                 * topo_buoyancy_N1B(iEdge) &
+                                                 / layerThickEdgeMean(k,iEdge)
+                  endif
+              end do
+
+          else if (config_topographic_wave_drag_scheme.EQ."ZAE") then
+              tensorScheme = .false.
+              call mpas_log_write("")
+              call mpas_log_write(" Topographic wave drag scheme is: Zaron and Egbert")
+              call mpas_log_write("")
+
+              ! Cze = Gamma H grad(H)^2 Nb Nbar / (8 pi^2 omega)
+              ! Gamma = 50
+              ! N(z) = N0 exp(z/1300)
+              ! N0 = 5.24 x 10^-3
+              ! z = 0 at surface
+              ! Nb is N(z) evaluated at z = -H
+              ! Nbar = 1300 N0 (1 - exp(-H/1300))/H
+
+              N0 = 5.24e-3_RKIND
+              L = 1300.0_RKIND
+              gam = 50.0_RKIND
+              do iEdge = 1, nEdgesAll
+                  kmax =  maxLevelEdgeTop(iEdge)
+                  topographicWaveDrag(iEdge) = 0.0_RKIND
+                  if (kmax > 0) then
+                      k = 1 ! Note: Only works for single layer case
+
+                      Nb = N0*exp(-layerThickEdgeMean(k,iEdge)/L)
+                      Nbar = N0*L*(1.0_RKIND - exp(-layerThickEdgeMean(k,iEdge)/L)) / layerThickEdgeMean(k,iEdge)
+
+                      topographicWaveDrag(iEdge) = topographicWaveDragCoeff * gam  &
+                                                 * bed_slope_edges(iEdge)**2 * Nbar * Nb &
+                                                 / (8.0_RKIND * omegaM2 * pii**2)
+                  endif
+              enddo
+          else if (config_topographic_wave_drag_scheme.EQ."LGF") then
+              tensorScheme = .true.
+              call mpas_log_write("")
+              call mpas_log_write(" Topographic wave drag scheme is: Local Generation Formula")
+              call mpas_log_write("")
+              do iEdge = 1, nEdgesAll
+                  kmax =  maxLevelEdgeTop(iEdge)
+                  topographicWaveDrag(iEdge) = 0.0_RKIND
+                  if (kmax > 0) then
+                      topographicWaveDrag(iEdge) = topographicWaveDragCoeff  &
+                                                 * (sqrt((topo_buoyancy_N1B(iEdge)**2 - omegaM2**2)   &
+                                                       * (topo_buoyancy_N1V(iEdge)**2 - omegaM2**2))) &
+                                                 / (4.0_RKIND*pii*omegaM2)
+                      normalCoeffTWD(iEdge) = (lonGradEdge(iEdge)**2) &
+                                         * (cos(angleEdge(iEdge))**2) + (latGradEdge(iEdge)**2)*(sin(angleEdge(iEdge))**2) &
+                                         - 2.0_RKIND*latGradEdge(iEdge)*lonGradEdge(iEdge)*sin(angleEdge(iEdge))*(cos(angleEdge(iEdge)))
+                      tangentialCoeffTWD(iEdge) = sin(angleEdge(iEdge)) &
+                                             * cos(angleEdge(iEdge)) &
+                                             * (lonGradEdge(iEdge)**2 - latGradEdge(iEdge)**2) &
+                                             + latGradEdge(iEdge)*lonGradEdge(iEdge)*(cos(angleEdge(iEdge))**2 - sin(angleEdge(iEdge))**2)
+                  endif
+              enddo
+          else 
+              call mpas_log_write("")
+
+              call mpas_log_write("Invalid parameter for config_topopgraphic_wave_drag_scheme. It must be one of: 'JSL', 'ZAE', or 'LGF'" , MPAS_LOG_CRIT)
+              call mpas_log_write("")
+          end if
+
+          ! Apply ramping factor to wave drag
+          do iEdge = 1, nEdgesAll
+              kmax =  maxLevelEdgeTop(iEdge)
+              if (kmax > 0) then
+                  k = 1
+                  topoDragTurnOff = (TANH((layerThickEdgeMean(k,iEdge)-topoDragCutoffDepth)/topoDragCutoffWidth)+1.0_RKIND)/2.0_RKIND
+                  topographicWaveDrag(iEdge) = topoDragTurnOff * topographicWaveDrag(iEdge)
+              endif
+          enddo
+          
       endif
 
       if (config_disable_vel_topographic_wave_drag) topographicWaveDragOn = .false.
diff --git a/components/mpas-ocean/src/shared/mpas_ocn_vel_self_attraction_loading.F b/components/mpas-ocean/src/shared/mpas_ocn_vel_self_attraction_loading.F
index 6b91daac4d37..90d1a4e4630d 100644
--- a/components/mpas-ocean/src/shared/mpas_ocn_vel_self_attraction_loading.F
+++ b/components/mpas-ocean/src/shared/mpas_ocn_vel_self_attraction_loading.F
@@ -103,7 +103,7 @@ module ocn_vel_self_attraction_loading
     real(kind=RKIND), dimension(:),   allocatable :: Snm, SnmRe, SnmIm
     real(kind=RKIND), dimension(:),   allocatable :: Snm_local, SnmRe_local, SnmIm_local
     real(kind=RKIND), dimension(:,:), allocatable :: Snm_local_reproSum, SnmRe_local_reproSum, SnmIm_local_reproSum
-    real(kind=RKIND), dimension(:),   allocatable :: sshSmoothed
+    real(kind=RKIND), dimension(:),   allocatable :: btrPressure
     real(kind=RKIND), dimension(:), pointer :: coastalSmoothingFactor
     integer, dimension(:,:), allocatable :: blockIdxForward
     integer, dimension(:,:), allocatable :: blockIdxInverse
@@ -128,7 +128,8 @@ module ocn_vel_self_attraction_loading
 !
 !-----------------------------------------------------------------------
 
-    subroutine ocn_compute_self_attraction_loading(domain, forcingPool, dminfo, ssh, err)!{{{
+    subroutine ocn_compute_self_attraction_loading(domain, forcingPool, &
+        dminfo, ssh, surfacePressure, err)!{{{
 
         !-----------------------------------------------------------------
         !
@@ -137,6 +138,7 @@ subroutine ocn_compute_self_attraction_loading(domain, forcingPool, dminfo, ssh,
         !-----------------------------------------------------------------
         type (dm_info), intent(in) :: dminfo
         real (kind=RKIND), dimension(:), intent(in) :: ssh
+        real (kind=RKIND), dimension(:), intent(in) :: surfacePressure
 
         !-----------------------------------------------------------------
         !
@@ -161,7 +163,8 @@ subroutine ocn_compute_self_attraction_loading(domain, forcingPool, dminfo, ssh,
         !-----------------------------------------------------------------
 
         integer :: iCell
-
+        real (kind=RKIND) :: rho0gInv
+        
         call mpas_timer_start('SAL Calculation')
 
         err = 0
@@ -170,8 +173,12 @@ subroutine ocn_compute_self_attraction_loading(domain, forcingPool, dminfo, ssh,
             return
         endif
 
+        rho0gInv = 1.0_RKIND / rho_sw / gravity
+
         do iCell = 1,nCellsAll
-          sshSmoothed(iCell) = coastalSmoothingFactor(iCell)*ssh(iCell)
+          ! compute SAL wrt. the full barotropic bottom pressure
+          btrPressure(iCell) = coastalSmoothingFactor(iCell) * ( &
+              ssh(iCell) + rho0gInv * surfacePressure(iCell) )
         enddo
 
         if (.not. config_use_parallel_self_attraction_loading) then
@@ -184,9 +191,9 @@ subroutine ocn_compute_self_attraction_loading(domain, forcingPool, dminfo, ssh,
 
         endif
 
-        ! ssh_sal currently computed on host but must be transferred
+        ! pgf_sal currently computed on host but must be transferred
         ! to device for later tendency calculations
-        !$acc update device(ssh_sal)
+        !$acc update device(pgf_sal)
 
         call mpas_timer_stop('SAL Calculation')
 
@@ -236,7 +243,7 @@ subroutine self_attraction_loading_serial(domain, dminfo, err)!{{{
         !-----------------------------------------------------------------
 
         ! SAL variables
-        real (kind=RKIND), dimension(:), pointer :: ssh_gg
+        real (kind=RKIND), dimension(:), pointer :: pgf_gg
 
         ! MPI Variables
         integer :: iCell, ilm, curProc
@@ -253,8 +260,8 @@ subroutine self_attraction_loading_serial(domain, dminfo, err)!{{{
             allocate(globalArray(nCellsGlobal), gatheredArray(nCellsGlobal))
         endif
 
-        ! Gather only the nCellsOwned from ssh (does not include Halos)
-        call MPI_GATHERV(sshSmoothed, nCellsOwned, MPI_DOUBLE, gatheredArray, nCellsPerProc, &
+        ! Gather only the nCellsOwned from pgf (does not include Halos)
+        call MPI_GATHERV(btrPressure, nCellsOwned, MPI_DOUBLE, gatheredArray, nCellsPerProc, &
                         nCellsDisplacement, MPI_DOUBLE, 0, dminfo % comm, err_tmp)
         err = ior(err, err_tmp)
         call mpas_timer_stop('Serial SAL: Gather')
@@ -262,20 +269,20 @@ subroutine self_attraction_loading_serial(domain, dminfo, err)!{{{
         ! Perform SAL calculation only on process 0
         if (curProc.eq.0) then
             call mpas_timer_start('Serial SAL: Forward')
-            allocate(ssh_gg(nGrid))
-            ssh_gg(:) = 0.0_dp
+            allocate(pgf_gg(nGrid))
+            pgf_gg(:) = 0.0_dp
             Slm(:)=(0.0_dp, 0.0_dp)
 
-            ! Rearrange ssh into CellID order
+            ! Rearrange into CellID order
             do iCell = 1,nCellsGlobal
                globalArray(indexToCellIDGathered(iCell)) = gatheredArray(iCell)
             enddo
 
-            ! Interpolate ssh onto Gaussian Grid
-            call interpolate( toColValues, toRowValues, toSValues, globalArray, ssh_gg)
+            ! Interpolate pgf onto Gaussian Grid
+            call interpolate( toColValues, toRowValues, toSValues, globalArray, pgf_gg)
             ! Perform spherical harmonic transform
 #ifdef USE_SHTNS
-            call Spat_to_SH(shtns_c, ssh_gg, Slm)
+            call Spat_to_SH(shtns_c, pgf_gg, Slm)
             call mpas_timer_stop('Serial SAL: Forward')
             call mpas_timer_start('Serial SAL: Inverse')
             ! Multiply each harmonic coefficient by the scaling factor
@@ -283,10 +290,10 @@ subroutine self_attraction_loading_serial(domain, dminfo, err)!{{{
                 Slm(ilm) = Slm(ilm) * LoveScaling(ilm)
             enddo
             ! Perform inverse spherical harmonic transform
-            call SH_to_spat(shtns_c, Slm, ssh_gg)
+            call SH_to_spat(shtns_c, Slm, pgf_gg)
 #endif
             ! Interpolate back to MPAS mesh
-            call interpolate( fromColValues, fromRowValues, fromSValues, ssh_gg, globalArray)
+            call interpolate( fromColValues, fromRowValues, fromSValues, pgf_gg, globalArray)
 
             ! Rearrange back to index order
             do iCell = 1,nCellsGlobal
@@ -300,18 +307,18 @@ subroutine self_attraction_loading_serial(domain, dminfo, err)!{{{
             call mpas_timer_stop('Serial SAL: Inverse')
         endif
 
-        ! Scatter back to ssh_sal
+        ! Scatter back to pgf_sal
         call mpas_timer_start('Serial SAL: Scatter')
         call MPI_SCATTERV(gatheredArray, nCellsPerProc, nCellsDisplacement, MPI_DOUBLE, &
-                    ssh_sal, nCellsAll, MPI_DOUBLE, 0, dminfo % comm, err_tmp)
+                    pgf_sal, nCellsAll, MPI_DOUBLE, 0, dminfo % comm, err_tmp)
         err = ior(err, err_tmp)
 
         if (curProc.eq.0) then
-          deallocate(globalArray, gatheredArray, ssh_gg)
+          deallocate(globalArray, gatheredArray, pgf_gg)
         endif
 
         ! Perform Halo exchange update
-        call mpas_dmpar_field_halo_exch(domain,'ssh_sal')
+        call mpas_dmpar_field_halo_exch(domain,'pgf_sal')
         call mpas_timer_stop('Serial SAL: Scatter')
 
     end subroutine self_attraction_loading_serial!}}}
@@ -394,13 +401,13 @@ subroutine self_attraction_loading_parallel(domain, dminfo)!{{{
                 ! Compute local integral contribution
                 if (config_parallel_self_attraction_loading_bfb) then
                   do iCell = endIdx, startIdx, -1
-                      SnmRe_local_reproSum(iCell,l) = SnmRe_local_reproSum(iCell,l) + sshSmoothed(iCell)*pmnm2(iCell)*complexFactorRe(iCell,m+1)
-                      SnmIm_local_reproSum(iCell,l) = SnmIm_local_reproSum(iCell,l) + sshSmoothed(iCell)*pmnm2(iCell)*complexFactorIm(iCell,m+1)
+                      SnmRe_local_reproSum(iCell,l) = SnmRe_local_reproSum(iCell,l) + btrPressure(iCell)*pmnm2(iCell)*complexFactorRe(iCell,m+1)
+                      SnmIm_local_reproSum(iCell,l) = SnmIm_local_reproSum(iCell,l) + btrPressure(iCell)*pmnm2(iCell)*complexFactorIm(iCell,m+1)
                   enddo
                 else
                   do iCell = endIdx, startIdx, -1
-                      SnmRe_local(l) = SnmRe_local(l) + sshSmoothed(iCell)*pmnm2(iCell)*complexFactorRe(iCell,m+1)
-                      SnmIm_local(l) = SnmIm_local(l) + sshSmoothed(iCell)*pmnm2(iCell)*complexFactorIm(iCell,m+1)
+                      SnmRe_local(l) = SnmRe_local(l) + btrPressure(iCell)*pmnm2(iCell)*complexFactorRe(iCell,m+1)
+                      SnmIm_local(l) = SnmIm_local(l) + btrPressure(iCell)*pmnm2(iCell)*complexFactorIm(iCell,m+1)
                   enddo
                 endif
 
@@ -416,13 +423,13 @@ subroutine self_attraction_loading_parallel(domain, dminfo)!{{{
                     ! Compute local integral contribution
                     if (config_parallel_self_attraction_loading_bfb) then
                       do iCell = endIdx, startIdx, -1
-                          SnmRe_local_reproSum(iCell,l) = SnmRe_local_reproSum(iCell,l) + sshSmoothed(iCell)*pmnm1(iCell)*complexFactorRe(iCell,m+1)
-                          SnmIm_local_reproSum(iCell,l) = SnmIm_local_reproSum(iCell,l) + sshSmoothed(iCell)*pmnm1(iCell)*complexFactorIm(iCell,m+1)
+                          SnmRe_local_reproSum(iCell,l) = SnmRe_local_reproSum(iCell,l) + btrPressure(iCell)*pmnm1(iCell)*complexFactorRe(iCell,m+1)
+                          SnmIm_local_reproSum(iCell,l) = SnmIm_local_reproSum(iCell,l) + btrPressure(iCell)*pmnm1(iCell)*complexFactorIm(iCell,m+1)
                       enddo
                     else
                       do iCell = endIdx, startIdx, -1
-                          SnmRe_local(l) = SnmRe_local(l) + sshSmoothed(iCell)*pmnm1(iCell)*complexFactorRe(iCell,m+1)
-                          SnmIm_local(l) = SnmIm_local(l) + sshSmoothed(iCell)*pmnm1(iCell)*complexFactorIm(iCell,m+1)
+                          SnmRe_local(l) = SnmRe_local(l) + btrPressure(iCell)*pmnm1(iCell)*complexFactorRe(iCell,m+1)
+                          SnmIm_local(l) = SnmIm_local(l) + btrPressure(iCell)*pmnm1(iCell)*complexFactorIm(iCell,m+1)
                       enddo
                     endif
 
@@ -447,13 +454,13 @@ subroutine self_attraction_loading_parallel(domain, dminfo)!{{{
                     ! Compute local integral contribution
                     if (config_parallel_self_attraction_loading_bfb) then
                       do iCell = endIdx, startIdx, -1
-                          SnmRe_local_reproSum(iCell,l) = SnmRe_local_reproSum(iCell,l) + sshSmoothed(iCell)*pmn(iCell)*complexFactorRe(iCell,m+1)
-                          SnmIm_local_reproSum(iCell,l) = SnmIm_local_reproSum(iCell,l) + sshSmoothed(iCell)*pmn(iCell)*complexFactorIm(iCell,m+1)
+                          SnmRe_local_reproSum(iCell,l) = SnmRe_local_reproSum(iCell,l) + btrPressure(iCell)*pmn(iCell)*complexFactorRe(iCell,m+1)
+                          SnmIm_local_reproSum(iCell,l) = SnmIm_local_reproSum(iCell,l) + btrPressure(iCell)*pmn(iCell)*complexFactorIm(iCell,m+1)
                       enddo
                     else
                       do iCell = endIdx, startIdx, -1
-                          SnmRe_local(l) = SnmRe_local(l) + sshSmoothed(iCell)*pmn(iCell)*complexFactorRe(iCell,m+1)
-                          SnmIm_local(l) = SnmIm_local(l) + sshSmoothed(iCell)*pmn(iCell)*complexFactorIm(iCell,m+1)
+                          SnmRe_local(l) = SnmRe_local(l) + btrPressure(iCell)*pmn(iCell)*complexFactorRe(iCell,m+1)
+                          SnmIm_local(l) = SnmIm_local(l) + btrPressure(iCell)*pmn(iCell)*complexFactorIm(iCell,m+1)
                       enddo
                     endif
 
@@ -514,7 +521,7 @@ subroutine self_attraction_loading_parallel(domain, dminfo)!{{{
         !!!!!!!!!!!!!!!!!!!!
 
         do iCell = 1,nCellsAll
-           ssh_sal(iCell) = 0.0_RKIND
+           pgf_sal(iCell) = 0.0_RKIND
         enddo
 
         do m = 0,nOrder
@@ -539,7 +546,7 @@ subroutine self_attraction_loading_parallel(domain, dminfo)!{{{
 
                 ! Sum together product of spherical harmonic functions and coefficients
                 do iCell = endIdx, startIdx, -1
-                    ssh_sal(iCell) = ssh_sal(iCell) + mFac*pmnm2(iCell)*(SnmRe(l)*complexExpRe(iCell,m+1)+SnmIm(l)*complexExpIm(iCell,m+1))
+                    pgf_sal(iCell) = pgf_sal(iCell) + mFac*pmnm2(iCell)*(SnmRe(l)*complexExpRe(iCell,m+1)+SnmIm(l)*complexExpIm(iCell,m+1))
                 enddo
 
                 !------------
@@ -553,7 +560,7 @@ subroutine self_attraction_loading_parallel(domain, dminfo)!{{{
 
                     ! Sum together product of spherical harmonic functions and coefficients
                     do iCell = endIdx, startIdx, -1
-                        ssh_sal(iCell) = ssh_sal(iCell) + mFac*pmnm1(iCell)*(SnmRe(l)*complexExpRe(iCell,m+1)+SnmIm(l)*complexExpIm(iCell,m+1))
+                        pgf_sal(iCell) = pgf_sal(iCell) + mFac*pmnm1(iCell)*(SnmRe(l)*complexExpRe(iCell,m+1)+SnmIm(l)*complexExpIm(iCell,m+1))
                     enddo
 
                 endif
@@ -576,7 +583,7 @@ subroutine self_attraction_loading_parallel(domain, dminfo)!{{{
 
                     ! Sum together product of spherical harmonic functions and coefficients
                     do iCell = endIdx, startIdx, -1
-                        ssh_sal(iCell) = ssh_sal(iCell) + mFac*pmn(iCell)*(SnmRe(l)*complexExpRe(iCell,m+1)+SnmIm(l)*complexExpIm(iCell,m+1))
+                        pgf_sal(iCell) = pgf_sal(iCell) + mFac*pmn(iCell)*(SnmRe(l)*complexExpRe(iCell,m+1)+SnmIm(l)*complexExpIm(iCell,m+1))
                     enddo
 
                 enddo ! n loop
@@ -609,6 +616,7 @@ subroutine ocn_vel_self_attraction_loading_init(domain,err)!{{{
         ! Config variables
         type (block_type), pointer :: block_ptr
         type (mpas_pool_type), pointer :: forcingPool
+        type (mpas_pool_type), pointer :: statePool
 
         ! NetCDF and weights file variables
         integer :: toNcId, toNsDimId, toRowId, toColId, toSId
@@ -641,8 +649,8 @@ subroutine ocn_vel_self_attraction_loading_init(domain,err)!{{{
         integer :: nCellBlock
         integer :: err_tmp
 
-        real(kind=RKIND) :: d
-        real(kind=RKIND) :: trans_start, trans_width
+        real(kind=RKIND), dimension(:), pointer :: ssh
+        real(kind=RKIND) :: depth, depth_tol
 
         call mpas_timer_start('SAL Init')
 
@@ -665,24 +673,24 @@ subroutine ocn_vel_self_attraction_loading_init(domain,err)!{{{
 
         block_ptr => domain % blocklist
         call mpas_pool_get_subpool(block_ptr % structs, 'forcing', forcingPool)
+        call mpas_pool_get_subpool(block_ptr % structs, 'state', statePool)
 
         call mpas_pool_get_array(forcingPool, 'coastalSmoothingFactor', &
                                                coastalSmoothingFactor)
+        call mpas_pool_get_array(statePool, 'ssh', ssh, 1)
 
-        ! Set up coastal ssh smoothing
-        allocate(sshSmoothed(nCellsAll))
-
-        if (config_self_attraction_loading_smoothing_width > 2.0_RKIND) then
-           trans_width = config_self_attraction_loading_smoothing_width*1000.0_RKIND
-           trans_start = 0.0_RKIND
-           do iCell = 1,nCellsAll
-               d = distanceToCoast(iCell)
-               coastalSmoothingFactor(iCell) = 0.5_RKIND*(tanh((d - trans_start - 0.5_RKIND * trans_width) / (0.2_RKIND *trans_width)) + 1.0_RKIND)
-           enddo
-        else
-           coastalSmoothingFactor(:) = 1.0_RKIND
-        endif
+        ! Set up coastal smoothing
+        allocate(btrPressure(nCellsAll))
 
+        ! Set-up smoothing mask for shallow regions where SAL should not
+        ! be computed across e.g. wet-dry cells
+        depth_tol = config_self_attraction_loading_depth_cutoff
+        do iCell = 1,nCellsAll
+           depth = ssh(iCell) + bottomDepth(iCell)
+           coastalSmoothingFactor(iCell) = 0.5_RKIND*(tanh( &
+              (depth - 0.5_RKIND*depth_tol) / (0.2_RKIND*depth_tol)) + 1.0_RKIND)
+        enddo
+        
         ! Setup Alarm for SAL
         alarmTime = mpas_get_clock_time(domain % clock, MPAS_START_TIME, ierr=err)
         call mpas_set_timeInterval(alarmTimeStep, timeString = &
diff --git a/components/mpas-ocean/src/shared/mpas_ocn_vel_tidal_potential.F b/components/mpas-ocean/src/shared/mpas_ocn_vel_tidal_potential.F
index 092c52c54d09..211538ac0a73 100644
--- a/components/mpas-ocean/src/shared/mpas_ocn_vel_tidal_potential.F
+++ b/components/mpas-ocean/src/shared/mpas_ocn_vel_tidal_potential.F
@@ -95,8 +95,8 @@ module ocn_vel_tidal_potential
    type(char_array), dimension(37) :: constituentList
    public :: char_array
 
-   ! Flag for ssh_sal
-   real(kind=RKIND) :: ssh_sal_on
+   ! Flag for pgf_sal
+   real(kind=RKIND) :: pgf_sal_on
 
 !***********************************************************************
 
@@ -162,7 +162,7 @@ subroutine ocn_vel_tidal_potential_tend(ssh, surfacePressure, tend, err)!{{{
 #ifdef MPAS_OPENACC
          !$acc parallel loop &
          !$acc    present(cellsOnEdge, maxLevelEdgeTop, dcEdge, edgeMask, &
-         !$acc            tidalPotEta, ssh, surfacePressure, tend, ssh_sal) &
+         !$acc            tidalPotEta, ssh, surfacePressure, tend, pgf_sal) &
          !$acc    private(cell1, cell2, invdcEdge, potentialGrad, k, kMax)
 #else
          !$omp parallel do schedule(runtime) &
@@ -176,8 +176,8 @@ subroutine ocn_vel_tidal_potential_tend(ssh, surfacePressure, tend, err)!{{{
 
             potentialGrad = - gravity*invdcEdge* &
                ( tidalPotEta(cell2) - tidalPotEta(cell1) &
-                + (1.0_RKIND - ssh_sal_on) * betaSelfAttrLoad * ( rho0gInv*(surfacePressure(cell2) - surfacePressure(cell1))) &
-                + ssh_sal_on * (ssh_sal(cell2) - ssh_sal(cell1)) )
+                + (1.0_RKIND - pgf_sal_on) * betaSelfAttrLoad * ( rho0gInv*(surfacePressure(cell2) - surfacePressure(cell1))) &
+                + pgf_sal_on * (pgf_sal(cell2) - pgf_sal(cell1)) )
 
             do k=1,kMax
                tend(k,iEdge) = tend(k,iEdge) - &
@@ -193,7 +193,7 @@ subroutine ocn_vel_tidal_potential_tend(ssh, surfacePressure, tend, err)!{{{
 #ifdef MPAS_OPENACC
          !$acc parallel loop &
          !$acc    present(cellsOnEdge, maxLevelEdgeTop, dcEdge, edgeMask, &
-         !$acc            tidalPotEta, ssh, surfacePressure, tend, ssh_sal) &
+         !$acc            tidalPotEta, ssh, surfacePressure, tend, pgf_sal) &
          !$acc    private(cell1, cell2, invdcEdge, potentialGrad, k, kMax)
 #else
          !$omp parallel do schedule(runtime) &
@@ -207,9 +207,9 @@ subroutine ocn_vel_tidal_potential_tend(ssh, surfacePressure, tend, err)!{{{
 
             potentialGrad = - gravity*invdcEdge* &
                (  tidalPotEta(cell2) - tidalPotEta(cell1) &
-                + (1.0_RKIND - ssh_sal_on) * betaSelfAttrLoad *( (ssh(cell2) - ssh(cell1)) &
-                                                         + rho0gInv*(surfacePressure(cell2) - surfacePressure(cell1))) &
-                + ssh_sal_on * (ssh_sal(cell2) - ssh_sal(cell1)) )
+               + (1.0_RKIND - pgf_sal_on) * betaSelfAttrLoad *( (ssh(cell2) - ssh(cell1)) & 
+                                                        + rho0gInv*(surfacePressure(cell2) - surfacePressure(cell1))) &
+               + pgf_sal_on * (pgf_sal(cell2) - pgf_sal(cell1)) )
 
             do k=1,kMax
                tend(k,iEdge) = tend(k,iEdge) - &
@@ -366,9 +366,9 @@ subroutine ocn_vel_tidal_potential_init(domain,err)!{{{
 
          tidalPotentialOff = .false.
          if (config_use_self_attraction_loading) then
-            ssh_sal_on = 1.0_RKIND
+            pgf_sal_on = 1.0_RKIND
          else
-            ssh_sal_on = 0.0_RKIND
+            pgf_sal_on = 0.0_RKIND
          endif
 
          betaSelfAttrLoad = config_self_attraction_and_loading_beta