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 cd15a29e84a6..0e410e7a1f74 100644
--- a/components/mpas-ocean/bld/namelist_files/namelist_defaults_mpaso.xml
+++ b/components/mpas-ocean/bld/namelist_files/namelist_defaults_mpaso.xml
@@ -147,9 +147,9 @@
 <config_Redi_horizontal_ramp_max ocn_grid="WCAtl12to45E2r4">40e3</config_Redi_horizontal_ramp_max>
 
 <!-- submesoscale_eddy_parameterization -->
-<config_submesoscale_enable>.false.</config_submesoscale_enable>
+<config_submesoscale_enable>.true.</config_submesoscale_enable>
 <config_submesoscale_tau>172800</config_submesoscale_tau>
-<config_submesoscale_Ce>0.06</config_submesoscale_Ce>
+<config_submesoscale_Ce>0.08</config_submesoscale_Ce>
 <config_submesoscale_Lfmin>1000.0</config_submesoscale_Lfmin>
 <config_submesoscale_ds_max>100000.0</config_submesoscale_ds_max>
 
diff --git a/components/mpas-ocean/cime_config/buildnml b/components/mpas-ocean/cime_config/buildnml
index c6b328ff730b..33553b239939 100755
--- a/components/mpas-ocean/cime_config/buildnml
+++ b/components/mpas-ocean/cime_config/buildnml
@@ -1118,6 +1118,8 @@ def buildnml(case, caseroot, compname):
         lines.append('    <var name="vertVelocityTop"/>')
         lines.append('    <var name="zMid"/>')
         lines.append('    <var name="atmosphericPressure"/>')
+        lines.append('    <var name="normalMLEvelocity"/>')
+        lines.append('    <var name="vertMLEBolusVelocityTop"/>')
         if not ocn_grid.startswith("oRRS1"):
             lines.append('    <var name="normalGMBolusVelocity"/>')
             lines.append('    <var name="vertGMBolusVelocityTop"/>')
diff --git a/components/mpas-ocean/src/Registry.xml b/components/mpas-ocean/src/Registry.xml
index a8ae9695b45e..3b5e0568795e 100644
--- a/components/mpas-ocean/src/Registry.xml
+++ b/components/mpas-ocean/src/Registry.xml
@@ -2875,8 +2875,8 @@
 			description="count of times redi limiter is invoked on a timestep"
 			packages="gm"
 		/>
-		<var name="gradBuoyEddy" type="real" dimensions="nVertLevelsP1 nEdges Time" units="s^-2"
-			description="horizontal gradient of buoyancy at cell edge and interfaces in vertical.  this is used for the GM parameterization and the submesoscale eddy parameterization"
+		<var name="gradDensityEddy" type="real" dimensions="nVertLevelsP1 nEdges Time" units="kg m^-4"
+			description="horizontal gradient of density at cell edge and interfaces in vertical.  this is used for the GM parameterization and the submesoscale eddy parameterization"
 			 packages="gm;submeso"
 		/>
 		<var name="mleVelocityZonal" type="real" dimensions="nVertLevels nCells Time" units="m s^-1"
diff --git a/components/mpas-ocean/src/analysis_members/Registry_time_series_stats_monthly_mean.xml b/components/mpas-ocean/src/analysis_members/Registry_time_series_stats_monthly_mean.xml
index 11942da077af..247dac423e6d 100644
--- a/components/mpas-ocean/src/analysis_members/Registry_time_series_stats_monthly_mean.xml
+++ b/components/mpas-ocean/src/analysis_members/Registry_time_series_stats_monthly_mean.xml
@@ -134,6 +134,8 @@
 			<var name="atmosphericPressure"/>
 			<var name="normalGMBolusVelocity"/>
 			<var name="vertGMBolusVelocityTop"/>
+			<var name="normalMLEvelocity"/>
+			<var name="vertMLEBolusVelocityTop"/>
 			<var name="GMBolusVelocityZonal"/>
 			<var name="GMBolusVelocityMeridional"/>
 			<var name="cGMphaseSpeed"/>
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 3850155fb4a3..2ecae166ff34 100644
--- a/components/mpas-ocean/src/shared/mpas_ocn_diagnostics_variables.F
+++ b/components/mpas-ocean/src/shared/mpas_ocn_diagnostics_variables.F
@@ -122,7 +122,7 @@ module ocn_diagnostics_variables
    real(kind=RKIND), dimension(:,:), pointer :: RediKappaScaling
    real(kind=RKIND), dimension(:,:), pointer :: RediKappaSfcTaper
    real(kind=RKIND), dimension(:,:), pointer :: k33
-   real(kind=RKIND), dimension(:,:), pointer :: gradBuoyEddy
+   real(kind=RKIND), dimension(:,:), pointer :: gradDensityEddy
    real(kind=RKIND), dimension(:), pointer   :: gmBolusKappa
    real(kind=RKIND), dimension(:), pointer   :: cGMphaseSpeed
    real(kind=RKIND), dimension(:,:,:), pointer :: slopeTriadUp, slopeTriadDown
@@ -271,8 +271,8 @@ subroutine ocn_diagnostics_variables_init(domain, jenkinsOn, hollandJenkinsOn, e
       end if
 
       if (config_use_GM .or. config_submesoscale_enable) then
-         call mpas_pool_get_array(diagnosticsPool, 'gradBuoyEddy', &
-                   gradBuoyEddy)
+         call mpas_pool_get_array(diagnosticsPool, 'gradDensityEddy', &
+                   gradDensityEddy)
       end if
 
       if ( config_compute_active_tracer_budgets ) then
@@ -781,7 +781,7 @@ subroutine ocn_diagnostics_variables_init(domain, jenkinsOn, hollandJenkinsOn, e
          !$acc                   )
       end if
       if (config_use_GM.or.config_submesoscale_enable) then
-         !$acc enter data create(gradBuoyEddy)
+         !$acc enter data create(gradDensityEddy)
       end if
       if ( config_compute_active_tracer_budgets ) then
          !$acc enter data create(                                         &
@@ -1029,7 +1029,7 @@ subroutine ocn_diagnostics_variables_destroy(err) !{{{
          !$acc                   )
       end if
       if (config_use_GM.or.config_submesoscale_enable) then
-         !$acc exit data delete(gradBuoyEddy)
+         !$acc exit data delete(gradDensityEddy)
       end if
       if ( config_compute_active_tracer_budgets ) then
          !$acc exit data delete(                                          &
@@ -1230,7 +1230,7 @@ subroutine ocn_diagnostics_variables_destroy(err) !{{{
                  RossbyRadius)
       end if
       if (config_use_GM.or.config_submesoscale_enable) then
-         nullify(gradBuoyEddy)
+         nullify(gradDensityEddy)
       end if
       if ( config_compute_active_tracer_budgets ) then
          nullify(activeTracerHorMixTendency,              &
diff --git a/components/mpas-ocean/src/shared/mpas_ocn_eddy_parameterization_helpers.F b/components/mpas-ocean/src/shared/mpas_ocn_eddy_parameterization_helpers.F
index 236da8d38849..8ac494147374 100644
--- a/components/mpas-ocean/src/shared/mpas_ocn_eddy_parameterization_helpers.F
+++ b/components/mpas-ocean/src/shared/mpas_ocn_eddy_parameterization_helpers.F
@@ -166,7 +166,7 @@ subroutine ocn_eddy_compute_buoyancy_gradient()
      do iEdge = 1, nEdges
         if (maxLevelEdgeTop(iEdge) .GE. minLevelEdgeBot(iEdge)) then
            do k = minLevelEdgeBot(iEdge), maxLevelEdgeTop(iEdge)+1
-              gradBuoyEddy(k,iEdge) = (gradDensityTopOfEdge(k,iEdge) - dDensityDzTopOfEdge(k,iEdge) &
+              gradDensityEddy(k,iEdge) = (gradDensityTopOfEdge(k,iEdge) - dDensityDzTopOfEdge(k,iEdge) &
                                              * gradZMidTopOfEdge(k,iEdge))
            end do
         end if
diff --git a/components/mpas-ocean/src/shared/mpas_ocn_gm.F b/components/mpas-ocean/src/shared/mpas_ocn_gm.F
index 336a160e3b9f..13bb0279de84 100644
--- a/components/mpas-ocean/src/shared/mpas_ocn_gm.F
+++ b/components/mpas-ocean/src/shared/mpas_ocn_gm.F
@@ -716,7 +716,7 @@ subroutine ocn_GM_compute_Bolus_velocity(statePool, &
                tridiagC(k - 1) = 2.0_RKIND*cGMphaseSpeed(iEdge)**2/layerThickEdgeMean(k, iEdge) &
                                  /(layerThickEdgeMean(k - 1, iEdge) + layerThickEdgeMean(k, iEdge))
                rightHandSide(k - 1) = gmBolusKappa(iEdge)*gmKappaScaling(k,iEdge)*gravity/rho_sw* &
-                                            gradBuoyEddy(k,iEdge) 
+                                            gradDensityEddy(k,iEdge) 
 
                ! Second to next to the last rows
                do k = minLevelEdgeBot(iEdge) + 2, maxLevelEdgeTop(iEdge) - 1
@@ -729,7 +729,7 @@ subroutine ocn_GM_compute_Bolus_velocity(statePool, &
                   tridiagC(k - 1) = 2.0_RKIND*cGMphaseSpeed(iEdge)**2/layerThickEdgeMean(k, iEdge) &
                                     /(layerThickEdgeMean(k - 1, iEdge) + layerThickEdgeMean(k, iEdge))
                   rightHandSide(k - 1) = gmBolusKappa(iEdge)*gmKappaScaling(k,iEdge)*gravity/rho_sw* &
-                                             gradBuoyEddy(k,iEdge)
+                                             gradDensityEddy(k,iEdge)
                end do
 
                ! Last row
@@ -741,7 +741,7 @@ subroutine ocn_GM_compute_Bolus_velocity(statePool, &
                tridiagB(k - 1) = -2.0_RKIND*cGMphaseSpeed(iEdge)**2/(layerThickEdgeMean(k - 1, iEdge) &
                                                                      *layerThickEdgeMean(k, iEdge)) - BruntVaisalaFreqTopEdge
                rightHandSide(k - 1) = gmBolusKappa(iEdge)*gmKappaScaling(k,iEdge)*gravity/rho_sw* &
-                                            gradBuoyEddy(k,iEdge)
+                                            gradDensityEddy(k,iEdge)
                ! Total number of rows
                N = maxLevelEdgeTop(iEdge) - minLevelEdgeBot(iEdge)
 
diff --git a/components/mpas-ocean/src/shared/mpas_ocn_submesoscale_eddies.F b/components/mpas-ocean/src/shared/mpas_ocn_submesoscale_eddies.F
index 99bab62abb5e..1236336d8b44 100644
--- a/components/mpas-ocean/src/shared/mpas_ocn_submesoscale_eddies.F
+++ b/components/mpas-ocean/src/shared/mpas_ocn_submesoscale_eddies.F
@@ -96,7 +96,7 @@ subroutine ocn_submesoscale_compute_velocity()
      !This is thickness weighted to allow for non uniform grid spacing, 
      !and centered on layer interfaces. The mixed layer goes from the 
      !ocean top layer to mid-depth of layer(indMLDedge)
-     !NOTE, gradBuoyEddy and BFV are defined at the top cell interface 
+     !NOTE, gradDensityEddy and BFV are defined at the top cell interface 
      !and not valid for minLevelCell, so properties for the top layer 
      !rely on the values at (minLevelCell+1).
      !$omp parallel
@@ -118,7 +118,7 @@ subroutine ocn_submesoscale_compute_velocity()
         hML = 0.5_RKIND*layerThickEdgeMean(minLevelEdgeBot(iEdge),iEdge)
         bvfML = bvfML + hML*bvfAv
         if (minLevelEdgeTop(iEdge) .ge. 1) then
-           gradBuoyML = hML*gradBuoyEddy(minLevelEdgeTop(iEdge)+1,iEdge)
+           gradBuoyML = hML*gradDensityEddy(minLevelEdgeTop(iEdge)+1,iEdge)
            bvfML = hML*bvfAv
         else
            cycle
@@ -131,10 +131,11 @@ subroutine ocn_submesoscale_compute_velocity()
                               max(BruntVaisalaFreqTop(k,cell2),1.0E-20_RKIND)))
            bvfML = bvfML + hAv*bvfAv
            hML = hML + hAv
-           gradBuoyML = gradBuoyML + hAv*gradBuoyEddy(k,iEdge)
+           gradBuoyML = gradBuoyML + hAv*gradDensityEddy(k,iEdge)
         end do
         bvfML = bvfML / (1.0E-20_RKIND + hML)
         gradBuoyML = gradBuoyML / (1.0E-20_RKIND + hML)
+        gradBuoyML = gravity*gradBuoyML/rho_sw !convert from density unit to buoyancy
 
         !compute depths and shape function