diff --git a/src/parameters/Parameters.jl b/src/parameters/Parameters.jl
index 652e48f49a..6a89854556 100644
--- a/src/parameters/Parameters.jl
+++ b/src/parameters/Parameters.jl
@@ -120,7 +120,8 @@ Base.@kwdef struct ClimaAtmosParameters{
     α_hyperdiff_tracer::FT
     # Vertical diffusion
     α_vert_diff_tracer::FT
-    # Gryanik b_m coefficient
+    # Gryanik coefficient
+    coeff_a_m_gryanik::FT
     coeff_b_m_gryanik::FT
 end
 
@@ -150,8 +151,9 @@ von_karman_const(ps::ACAP) =
 
 # ------ MOST (Monin–Obukhov) stability-function coefficients ------
 
-# Gryanik b_m
+# Gryanik
 # needed because surface_fluxes_params defaults to BusingerParams
+coefficient_a_m_gryanik(ps::ACAP) = ps.coeff_a_m_gryanik
 coefficient_b_m_gryanik(ps::ACAP) = ps.coeff_b_m_gryanik
 
 # Insolation parameters
diff --git a/src/parameters/create_parameters.jl b/src/parameters/create_parameters.jl
index fad267ddb4..3a4bb76afa 100644
--- a/src/parameters/create_parameters.jl
+++ b/src/parameters/create_parameters.jl
@@ -42,7 +42,8 @@ function ClimaAtmosParameters(
         SurfaceFluxesParameters(toml_dict, UF.BusingerParams)
     SFP = typeof(surface_fluxes_params)
 
-    # Fetch Gryanik b_m coefficient (since surface_fluxes_params defaults to BusingerParams)
+    # Fetch Gryanik coefficients (since surface_fluxes_params defaults to BusingerParams)
+    coeff_a_m_gryanik_val = UF.GryanikParams(FT).a_m
     coeff_b_m_gryanik_val = UF.GryanikParams(FT).b_m
 
     surface_temp_params = SurfaceTemperatureParameters(toml_dict)
@@ -110,6 +111,7 @@ function ClimaAtmosParameters(
         surface_temp_params,
         vert_diff_params,
         external_forcing_params,
+        coeff_a_m_gryanik = coeff_a_m_gryanik_val,
         coeff_b_m_gryanik = coeff_b_m_gryanik_val,
     )
 end
diff --git a/src/prognostic_equations/eddy_diffusion_closures.jl b/src/prognostic_equations/eddy_diffusion_closures.jl
index af65232463..4b0642168e 100644
--- a/src/prognostic_equations/eddy_diffusion_closures.jl
+++ b/src/prognostic_equations/eddy_diffusion_closures.jl
@@ -309,7 +309,8 @@ function mixing_length_lopez_gomez_2020(
     # MOST stability function coefficients
     most_a_m = sf_params.ufp.a_m # Businger a_m
     most_b_m = sf_params.ufp.b_m # Businger b_m
-    most_g_m = CAP.coefficient_b_m_gryanik(params)  # Gryanik b_m
+    gryanik_a_m = CAP.coefficient_a_m_gryanik(params)  # Gryanik a_m
+    gryanik_b_m = CAP.coefficient_b_m_gryanik(params)  # Gryanik b_m
 
     # l_z: Geometric distance from the surface
     l_z = ᶜz - z_sfc
@@ -359,16 +360,14 @@ function mixing_length_lopez_gomez_2020(
         zeta = l_z / obukhov_len_safe_stable # zeta >= 0
 
         # Stable/neutral-regime correction after Gryanik (2020):
-        #     φ_m = 1 + a_m ζ / (1 + g_m ζ)^(2/3),
+        #     φ_m = 1 + a_m ζ / (1 + b_m ζ)^(2/3),
         # a nonlinear refinement to the Businger formulation.
-        phi_m_denom_term = (1 + most_g_m * zeta)
+        # Optimization: (1+b_m*ζ)^(2/3) -> cbrt((1+b_m*ζ)^2)
         # Guard against a negative base in the fractional power
-        # (theoretically impossible for ζ ≥ 0 and g_m > 0, retained for robustness).
-        phi_m_denom_cubed_sqrt = cbrt(phi_m_denom_term)
-        phi_m_denom =
-            max(phi_m_denom_cubed_sqrt * phi_m_denom_cubed_sqrt, eps_FT) # (val)^(2/3)
+        # (theoretically impossible for ζ ≥ 0 and b_m > 0, retained for robustness).
+        phi_m_denom = max(cbrt((1 + gryanik_b_m * zeta)^2), eps_FT)
 
-        phi_m = 1 + (most_a_m * zeta) / phi_m_denom
+        phi_m = 1 + (gryanik_a_m * zeta) / phi_m_denom
 
         # Stable-regime correction factor: 1 / φ_m.
         # phi_m should be >= 1 for stable/neutral