diff --git a/components/mpas-ocean/bld/build-namelist b/components/mpas-ocean/bld/build-namelist
index 1bee62303852..271f46c6cc54 100755
--- a/components/mpas-ocean/bld/build-namelist
+++ b/components/mpas-ocean/bld/build-namelist
@@ -579,6 +579,7 @@ add_default($nl, 'config_Redi_N2_based_taper_limit_term1');
 add_default($nl, 'config_use_GM');
 add_default($nl, 'config_GM_closure');
 add_default($nl, 'config_GM_constant_kappa');
+add_default($nl, 'config_GM_gravWaveSpeed');
 add_default($nl, 'config_GM_constant_gravWaveSpeed');
 add_default($nl, 'config_GM_spatially_variable_min_kappa');
 add_default($nl, 'config_GM_spatially_variable_max_kappa');
diff --git a/components/mpas-ocean/bld/build-namelist-section b/components/mpas-ocean/bld/build-namelist-section
index 7812c81af9f0..26e2410dac9b 100644
--- a/components/mpas-ocean/bld/build-namelist-section
+++ b/components/mpas-ocean/bld/build-namelist-section
@@ -122,6 +122,7 @@ add_default($nl, 'config_Redi_N2_based_taper_limit_term1');
 add_default($nl, 'config_use_GM');
 add_default($nl, 'config_GM_closure');
 add_default($nl, 'config_GM_constant_kappa');
+add_default($nl, 'config_GM_gravWaveSpeed');
 add_default($nl, 'config_GM_constant_gravWaveSpeed');
 add_default($nl, 'config_GM_spatially_variable_min_kappa');
 add_default($nl, 'config_GM_spatially_variable_max_kappa');
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 1740b01a7e2d..7820c2fe081e 100644
--- a/components/mpas-ocean/bld/namelist_files/namelist_defaults_mpaso.xml
+++ b/components/mpas-ocean/bld/namelist_files/namelist_defaults_mpaso.xml
@@ -173,6 +173,7 @@
 <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_gravWaveSpeed>'constant'</config_GM_gravWaveSpeed>
 <config_GM_constant_gravWaveSpeed>0.3</config_GM_constant_gravWaveSpeed>
 <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>
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 931515e45375..4d7d2c74c891 100644
--- a/components/mpas-ocean/bld/namelist_files/namelist_definition_mpaso.xml
+++ b/components/mpas-ocean/bld/namelist_files/namelist_definition_mpaso.xml
@@ -520,6 +520,14 @@ Valid values: MISSING POSSIBLE VALUES
 Default: Defined in namelist_defaults.xml
 </entry>
 
+<entry id="+config_GM_gravWaveSpeed" type="char*1024"
+	category="GM_eddy_parameterization" group="GM_eddy_parameterization">
+Determines how GM gravity wave speed is computed. If 'constant' then use config_GM_constant_gravWaveSpeed. If 'computed' then compute at every edge at every time step using the Brunt-Vaisala frequency
+
+Valid values: 'constant', 'computed'
+Default: Defined in namelist_defaults.xml
+</entry>
+
 <entry id="config_GM_constant_gravWaveSpeed" 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).
diff --git a/components/mpas-ocean/src/Registry.xml b/components/mpas-ocean/src/Registry.xml
index 957b9f29909b..5ac3449b49d8 100644
--- a/components/mpas-ocean/src/Registry.xml
+++ b/components/mpas-ocean/src/Registry.xml
@@ -353,6 +353,10 @@
 		<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_gravWaveSpeed" type="character" default_value="constant" units="NA"
+					description="Determines how GM gravity wave speed is computed. If 'constant' then use config_GM_constant_gravWaveSpeed. 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_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)."
 					possible_values="Positive real numbers"
diff --git a/components/mpas-ocean/src/shared/mpas_ocn_gm.F b/components/mpas-ocean/src/shared/mpas_ocn_gm.F
index c971a4737adb..a429ba3c8c87 100644
--- a/components/mpas-ocean/src/shared/mpas_ocn_gm.F
+++ b/components/mpas-ocean/src/shared/mpas_ocn_gm.F
@@ -53,6 +53,7 @@ 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) :: gm_constant_gravWaveSpeed, gm_computed_gravWaveSpeed_on
 
 !***********************************************************************
 
@@ -132,7 +133,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 +544,12 @@ 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
+               ltSum = epsGM
                do k = minLevelEdgeBot(iEdge)+1, maxLevelEdgeTop(iEdge)-1
 
                   lt1 = 0.5_RKIND*(layerThickness(k,cell1) + layerThickness(k-1,cell1))
@@ -556,9 +557,13 @@ 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 gravity wave speed, either a specified constant, or
+               ! computed from the Brunt-Vaisala frequency.
+               c_min = gm_constant_gravWaveSpeed + gm_computed_gravWaveSpeed_on*max(0.01_RKIND,sumN2/ltSum*(0.5*(lt1+lt2)))
+               cGMphaseSpeed(iEdge) = max(c_min,                          &
                                           sumN2/(config_GM_spatially_variable_baroclinic_mode*3.141592_RKIND))
 
             end do
@@ -910,6 +915,17 @@ subroutine ocn_GM_init(domain, err)!{{{
             sfcTaperFactor = 0.0_RKIND
          end if
 
+         if (config_gm_gravWaveSpeed=='constant') then
+            gm_constant_gravWaveSpeed = config_gm_constant_gravWaveSpeed
+            gm_computed_gravWaveSpeed_on = 0.0_RKIND
+         elseif (config_gm_gravWaveSpeed=='computed') then
+            gm_constant_gravWaveSpeed = 0.0_RKIND
+            gm_computed_gravWaveSpeed_on = 1.0_RKIND
+         else
+            call mpas_log_write('config_gm_gravWaveSpeed 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.