diff --git a/components/mpas-ocean/bld/build-namelist b/components/mpas-ocean/bld/build-namelist
index 8138bbb0a733..6f09fa34ba4a 100755
--- a/components/mpas-ocean/bld/build-namelist
+++ b/components/mpas-ocean/bld/build-namelist
@@ -580,7 +580,8 @@ add_default($nl, 'config_Redi_min_layers_diag_terms');
 add_default($nl, 'config_use_GM');
 add_default($nl, 'config_GM_closure');
 add_default($nl, 'config_GM_constant_kappa');
-add_default($nl, 'config_GM_constant_gravWaveSpeed');
+add_default($nl, 'config_GM_constant_bclModeSpeed');
+add_default($nl, 'config_GM_minBclModeSpeed_method');
 add_default($nl, 'config_GM_spatially_variable_min_kappa');
 add_default($nl, 'config_GM_spatially_variable_max_kappa');
 add_default($nl, 'config_GM_spatially_variable_baroclinic_mode');
diff --git a/components/mpas-ocean/bld/build-namelist-section b/components/mpas-ocean/bld/build-namelist-section
index c612ebb0de10..96e924bf0207 100644
--- a/components/mpas-ocean/bld/build-namelist-section
+++ b/components/mpas-ocean/bld/build-namelist-section
@@ -123,7 +123,8 @@ add_default($nl, 'config_Redi_min_layers_diag_terms');
 add_default($nl, 'config_use_GM');
 add_default($nl, 'config_GM_closure');
 add_default($nl, 'config_GM_constant_kappa');
-add_default($nl, 'config_GM_constant_gravWaveSpeed');
+add_default($nl, 'config_GM_constant_bclModeSpeed');
+add_default($nl, 'config_GM_minBclModeSpeed_method');
 add_default($nl, 'config_GM_spatially_variable_min_kappa');
 add_default($nl, 'config_GM_spatially_variable_max_kappa');
 add_default($nl, 'config_GM_spatially_variable_baroclinic_mode');
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 1aa360ab0839..fe16b02f9efe 100644
--- a/components/mpas-ocean/bld/namelist_files/namelist_defaults_mpaso.xml
+++ b/components/mpas-ocean/bld/namelist_files/namelist_defaults_mpaso.xml
@@ -174,7 +174,8 @@
 <config_GM_constant_kappa ocn_forcing="datm_forced_restoring" ocn_grid="WCAtl12to45E2r4">600.0</config_GM_constant_kappa>
 <config_GM_constant_kappa ocn_forcing="datm_forced_restoring" ocn_grid="SOwISC12to60E2r4">600.0</config_GM_constant_kappa>
 <config_GM_constant_kappa ocn_forcing="datm_forced_restoring" ocn_grid="ECwISC30to60E2r1">600.0</config_GM_constant_kappa>
-<config_GM_constant_gravWaveSpeed>0.3</config_GM_constant_gravWaveSpeed>
+<config_GM_constant_bclModeSpeed>0.3</config_GM_constant_bclModeSpeed>
+<config_GM_minBclModeSpeed_method>'constant'</config_GM_minBclModeSpeed_method>
 <config_GM_spatially_variable_min_kappa>300.0</config_GM_spatially_variable_min_kappa>
 <config_GM_spatially_variable_max_kappa>1800.0</config_GM_spatially_variable_max_kappa>
 <config_GM_spatially_variable_baroclinic_mode>3.0</config_GM_spatially_variable_baroclinic_mode>
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 f91f4bb4a984..d9b4674b79ac 100644
--- a/components/mpas-ocean/bld/namelist_files/namelist_definition_mpaso.xml
+++ b/components/mpas-ocean/bld/namelist_files/namelist_definition_mpaso.xml
@@ -487,7 +487,7 @@ Default: Defined in namelist_defaults.xml
 
 <entry id="config_Redi_N2_based_taper_limit_term1" type="logical"
 	category="Redi_isopycnal_mixing" group="Redi_isopycnal_mixing">
-If true, the N2 lim iting is applied to the horizontal diffusion term
+If true, the N2 limiting is applied to the horizontal diffusion term
 
 Valid values: .true. or .false.
 Default: Defined in namelist_defaults.xml
@@ -528,14 +528,22 @@ Valid values: MISSING POSSIBLE VALUES
 Default: Defined in namelist_defaults.xml
 </entry>
 
-<entry id="config_GM_constant_gravWaveSpeed" type="real"
+<entry id="config_GM_constant_bclModeSpeed" type="real"
 	category="GM_eddy_parameterization" group="GM_eddy_parameterization">
-Gravity wave speed for the vertical stream function boundary value problem. This appears as $c$ in eqn 16a of Ferrari et al. 2010 (https://doi.org/10.1016/j.ocemod.2010.01.004).
+The parameter setting for the first baroclinic mode speed for the vertical stream function boundary value problem. This appears as $c$ in eqn 16a of Ferrari et al. 2010 (https://doi.org/10.1016/j.ocemod.2010.01.004).
 
 Valid values: Positive real numbers
 Default: Defined in namelist_defaults.xml
 </entry>
 
+<entry id="config_GM_minBclModeSpeed_method" type="char*1024"
+	category="GM_eddy_parameterization" group="GM_eddy_parameterization">
+Determines how the GM setting for the minimum of the first baroclinic mode speed is computed. If 'constant' then use config_GM_constant_bclModeSpeed. If 'computed' then compute at every edge at every time step using the Brunt-Vaisala frequency
+
+Valid values: 'constant' and 'computed'
+Default: Defined in namelist_defaults.xml
+</entry>
+
 <entry id="config_GM_spatially_variable_min_kappa" type="real"
 	category="GM_eddy_parameterization" group="GM_eddy_parameterization">
 minimum value of bolus diffusivity for spatially variable GM schemes. Used for all choices of config_GM_closure other than 'constant'.
diff --git a/components/mpas-ocean/src/Registry.xml b/components/mpas-ocean/src/Registry.xml
index b5b77183e342..25c8cb34db75 100644
--- a/components/mpas-ocean/src/Registry.xml
+++ b/components/mpas-ocean/src/Registry.xml
@@ -357,10 +357,14 @@
 		<nml_option name="config_GM_constant_kappa" type="real" default_value="600.0" units="m^2 s^{-1}"
 					description="Coefficient of standard GM parametrization of eddy transport (Bolus component), $\kappa$. Only used when config_GM_closure is set to constant."
 		/>
-		<nml_option name="config_GM_constant_gravWaveSpeed" type="real" default_value="0.3" units="m/s"
-					description="Gravity wave speed for the vertical stream function boundary value problem. This appears as $c$ in eqn 16a of Ferrari et al. 2010 (https://doi.org/10.1016/j.ocemod.2010.01.004)."
+		<nml_option name="config_GM_constant_bclModeSpeed" type="real" default_value="0.3" units="m/s"
+					description="The parameter setting for the first baroclinic mode speed for the vertical stream function boundary value problem. This appears as $c$ in eqn 16a of Ferrari et al. 2010 (https://doi.org/10.1016/j.ocemod.2010.01.004)."
 					possible_values="Positive real numbers"
 		/>
+		<nml_option name="config_GM_minBclModeSpeed_method" type="character" default_value="constant" units="NA"
+					description="Determines how the GM setting for the minimum of the first baroclinic mode speed is computed. If 'constant' then use config_GM_constant_bclModeSpeed. If 'computed' then compute at every edge at every time step using the Brunt-Vaisala frequency"
+					possible_values="'constant' and 'computed'"
+		/>
 		<nml_option name="config_GM_spatially_variable_min_kappa" type="real" default_value="300.0" units="m^2 s^{-1}"
 					description="minimum value of bolus diffusivity for spatially variable GM schemes. Used for all choices of config_GM_closure other than 'constant'."
 					possible_values="values around 100s"
diff --git a/components/mpas-ocean/src/shared/mpas_ocn_gm.F b/components/mpas-ocean/src/shared/mpas_ocn_gm.F
index c971a4737adb..5c6bbf500c28 100644
--- a/components/mpas-ocean/src/shared/mpas_ocn_gm.F
+++ b/components/mpas-ocean/src/shared/mpas_ocn_gm.F
@@ -53,6 +53,8 @@ module ocn_gm
    logical :: local_config_GM_kappa_lat_depth_variable
    logical :: local_config_GM_compute_EdenGreatbatch
    real(kind=RKIND) :: slopeTaperFactor, sfcTaperFactor, rediGMinitValue
+   real(kind=RKIND) :: local_config_GM_constant_bclModeSpeed
+   real(kind=RKIND) :: gm_minBclModeSpeed_constant, gm_minBclModeSpeed_compute_on
 
 !***********************************************************************
 
@@ -132,7 +134,7 @@ subroutine ocn_GM_compute_Bolus_velocity(statePool, &
       real(kind=RKIND)                 :: sumRi, RiTopOfEdge, zEdge, zMLD, sfcTaper
       real(kind=RKIND) :: dcEdgeInv, drhoDx, drhoDT, drhoDS, dTdx, dSdx
       real(kind=RKIND) :: slopeTaperUp, slopeTaperDown, sfcTaperUp, sfcTaperDown, invAreaCell
-      real(kind=RKIND) :: lt1, lt2
+      real(kind=RKIND) :: lt1, lt2, c_min
       real(kind=RKIND) :: sigma, Lr, Length, L_rhines, shearEdgeInv
       real(kind=RKIND), dimension(:), allocatable :: dzTop, dTdzTop, dSdzTop, k33Norm
       real(kind=RKIND) :: c_Visbeck   ! baroclinic wave speed from Visbeck parameterization
@@ -543,12 +545,14 @@ subroutine ocn_GM_compute_Bolus_velocity(statePool, &
          ! config_GM_spatially_variable_baroclinic_mode
          if (local_config_GM_lat_variable_c2) then
             !$omp parallel
-            !$omp do schedule(runtime) private(k, cell1, cell2, sumN2, lt1, lt2)
+            !$omp do schedule(runtime) private(k, cell1, cell2, sumN2, lt1, lt2, ltSum, c_min)
             do iEdge = 1, nEdges
                cell1 = cellsOnEdge(1, iEdge)
                cell2 = cellsOnEdge(2, iEdge)
-               sumN2 = 0.0
-
+               sumN2 = 0.0_RKIND
+               lt1 = 0.0_RKIND
+               lt2 = 0.0_RKIND
+               ltSum = epsGM
                do k = minLevelEdgeBot(iEdge)+1, maxLevelEdgeTop(iEdge)-1
 
                   lt1 = 0.5_RKIND*(layerThickness(k,cell1) + layerThickness(k-1,cell1))
@@ -556,9 +560,22 @@ subroutine ocn_GM_compute_Bolus_velocity(statePool, &
                   sumN2 = sumN2 + 0.5_RKIND*(lt1*sqrt(max(BruntVaisalaFreqTop(k,cell1),0.0_RKIND)) + &
                                 lt2*sqrt(max(BruntVaisalaFreqTop(k,cell2),0.0_RKIND)))
 
+                  ltSum = ltSum + 0.5*(lt1+lt2)
                end do
 
-               cGMphaseSpeed(iEdge) = max(config_GM_constant_gravWaveSpeed,                          &
+               ! Compute the minimum allowed speed of the first baroclinic mode. 
+               ! If config_GM_minBclModeSpeed_method='constant' use local_config_GM_constant_bclModeSpeed.
+               ! If config_GM_minBclModeSpeed_method='computed' use Brunt-Vaisala frequency on this edge.
+               ! See initialization of these variables in this modules init routine.
+               ! The value 0.5*(lt1+lt2) is meant to be an estimate of the largest thickness
+               ! in the column, and c_min is constructed such that we satisfy
+               ! the condition c/N < dz.  The condition is an estimate, but
+               ! tries to insure that the streamfunction is well resolved even
+               ! for the thickest layers.
+               c_min = gm_minBclModeSpeed_constant + gm_minBclModeSpeed_compute_on*max(0.01_RKIND,sumN2/ltSum*(0.5*(lt1+lt2)))
+
+               ! Compute the speed of the first baroclinic mode from the Brunt-Vaisala frequency.
+               cGMphaseSpeed(iEdge) = max(c_min,                          &
                                           sumN2/(config_GM_spatially_variable_baroclinic_mode*3.141592_RKIND))
 
             end do
@@ -569,7 +586,7 @@ subroutine ocn_GM_compute_Bolus_velocity(statePool, &
             !$omp parallel
             !$omp do schedule(runtime)
             do iEdge = 1, nEdges
-               cGMphaseSpeed(iEdge) = config_GM_constant_gravWaveSpeed
+               cGMphaseSpeed(iEdge) = local_config_GM_constant_bclModeSpeed
             end do
             !$omp end do
             !$omp end parallel
@@ -910,6 +927,19 @@ subroutine ocn_GM_init(domain, err)!{{{
             sfcTaperFactor = 0.0_RKIND
          end if
 
+         local_config_GM_constant_bclModeSpeed = config_GM_constant_bclModeSpeed
+
+         if (config_GM_minBclModeSpeed_method=='constant') then
+            gm_minBclModeSpeed_compute_on = 0.0_RKIND
+            gm_minBclModeSpeed_constant = config_GM_constant_bclModeSpeed
+         elseif (config_GM_minBclModeSpeed_method=='computed') then
+            gm_minBclModeSpeed_compute_on = 1.0_RKIND
+            gm_minBclModeSpeed_constant = 0.0_RKIND
+         else
+            call mpas_log_write('config_GM_minBclModeSpeed_method must be constant or computed', &
+                                 MPAS_LOG_CRIT)
+         end if
+
          RediGMinitValue = 1.0_RKIND
          if (config_GM_closure == 'constant') then
             local_config_GM_lat_variable_c2 = .false.
diff --git a/components/mpas-ocean/src/shared/mpas_ocn_test.F b/components/mpas-ocean/src/shared/mpas_ocn_test.F
index 48ac21644934..3b6859da6fab 100644
--- a/components/mpas-ocean/src/shared/mpas_ocn_test.F
+++ b/components/mpas-ocean/src/shared/mpas_ocn_test.F
@@ -316,7 +316,7 @@ subroutine ocn_init_gm_test_functions(meshPool, scratchPool)!{{{
       ! zBot is location we apply boundary conditions on the ODE for stream function.
       zBot = zMax
 
-      L = config_GM_constant_gravWaveSpeed * sqrt(rho_sw * zMax / gravity / config_eos_linear_alpha / config_gm_analytic_temperature3)
+      L = config_GM_constant_bclModeSpeed * sqrt(rho_sw * zMax / gravity / config_eos_linear_alpha / config_gm_analytic_temperature3)
       R = - config_GM_constant_kappa * config_gm_analytic_temperature2 * zMax / config_gm_analytic_temperature3 &
         / config_gm_analytic_ymax
       c1 = R*(1-exp(-zBot/L))/(exp(zBot/L) - exp(-zBot/L))