Merge pull request #4143 from CliMA/zs/mixing_length_uf

use universal functions in mixing length

Author: szy21

Project.toml

@@ -68,8 +68,8 @@ SciMLBase = "~2.108, ~2.109, ~2.110, ~2.111, ~2.112"
 SparseMatrixColorings = "0.4.20"
 StaticArrays = "1.9"
 Statistics = "1"
-SurfaceFluxes = "0.12.3, 0.13, 0.14.1"
-Thermodynamics = "0.14.1, 0.15"
+SurfaceFluxes = "0.14.1"
+Thermodynamics = "0.14.2, 0.15"
 UnrolledUtilities = "0.1.9"
 YAML = "0.4"
 julia = "1.9"

reproducibility_tests/ref_counter.jl

@@ -1,4 +1,4 @@
-283
+284
 
 # **README**
 #
@@ -20,6 +20,10 @@
 
 
 #=
+284
+- Use universal functions in SurfaceFluxes in the mixing length closure,
+which changes the closure in the stable case (from Gryanik to Businger)
+
 283
 - Change the Jacobian terms related to dp_drhoq_tot
 
@@ -39,7 +43,8 @@
 - Add ∂/∂q elements to Jacobian
 
 277
-- Update to use [email protected]. Supports Charnock-parameterization for aerodynamic roughness (default is still user-prescribed ScalarRoughness). SurfaceFluxes catch for neutrally-stable boundary layers (ζ ≈ 0) removed. 
+- Update to use [email protected]. Supports Charnock-parameterization for aerodynamic roughness 
+(default is still user-prescribed ScalarRoughness). SurfaceFluxes catch for neutrally-stable boundary layers (ζ ≈ 0) removed. 
 
 276
 - Update prognostic EDMF boundary conditions: apply equal surface fluxes to the

src/parameters/Parameters.jl

@@ -120,9 +120,6 @@ Base.@kwdef struct ClimaAtmosParameters{
     α_hyperdiff_tracer::FT
     # Vertical diffusion
     α_vert_diff_tracer::FT
-    # Gryanik coefficient
-    coeff_a_m_gryanik::FT
-    coeff_b_m_gryanik::FT
 end
 
 Base.eltype(::ClimaAtmosParameters{FT}) where {FT} = FT
@@ -151,11 +148,6 @@ von_karman_const(ps::ACAP) =
 
 # ------ MOST (Monin–Obukhov) stability-function coefficients ------
 
-# 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
 day(ps::ACAP) = IP.day(insolation_params(ps))
 tot_solar_irrad(ps::ACAP) = IP.tot_solar_irrad(insolation_params(ps))

src/parameters/create_parameters.jl

@@ -42,10 +42,6 @@ function ClimaAtmosParameters(
         SurfaceFluxesParameters(toml_dict, UF.BusingerParams)
     SFP = typeof(surface_fluxes_params)
 
-    # 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)
     STP = typeof(surface_temp_params)
 
@@ -111,8 +107,6 @@ 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
 

src/prognostic_equations/eddy_diffusion_closures.jl

@@ -6,6 +6,7 @@ import StaticArrays as SA
 import Thermodynamics.Parameters as TDP
 import ClimaCore.Geometry as Geometry
 import ClimaCore.Fields as Fields
+import SurfaceFluxes.UniversalFunctions as UF
 
 """
     buoyancy_gradients(
@@ -294,17 +295,9 @@ function mixing_length_lopez_gomez_2020(
 
     c_m = CAP.tke_ed_coeff(turbconv_params)
     c_d = CAP.tke_diss_coeff(turbconv_params)
-    smin_ub = CAP.smin_ub(turbconv_params)
-    smin_rm = CAP.smin_rm(turbconv_params)
     c_b = CAP.static_stab_coeff(turbconv_params)
     vkc = CAP.von_karman_const(params)
 
-    # 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
-    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
     # Ensure l_z is non-negative when ᶜz is numerically smaller than z_sfc.
@@ -330,45 +323,12 @@ function mixing_length_lopez_gomez_2020(
     # and approaches 0 when l_z → 0.
     l_W_base = vkc * l_z / l_W_denom
 
-    if obukhov_length < FT(0) # Unstable case
-        obukhov_len_safe = min(obukhov_length, -eps_FT) # Ensure L < 0
-        zeta = l_z / obukhov_len_safe # Stability parameter zeta = z/L (<0)
-
-        # Calculate MOST term (1 - b_m * zeta)
-        # Since zeta is negative, this term is > 1
-        inner_term = 1 - most_b_m * zeta
-
-        # Numerical safety check – by theory the value is ≥ 1.
-        inner_term_safe = max(inner_term, eps_FT)
-
-        # Unstable-regime correction factor:
-        #     (1 − b_m ζ)^(1/4) = φ_m⁻¹,
-        # where φ_m is the Businger stability function φ_m = (1 − b_m ζ)^(-1/4).
-        stability_correction = sqrt(sqrt(inner_term_safe))
-        l_W = l_W_base * stability_correction
+    obukhov_len_safe =
+        obukhov_length < FT(0) ? min(obukhov_length, -eps_FT) : max(obukhov_length, eps_FT)
+    zeta = l_z / obukhov_len_safe # Stability parameter zeta
+    phi_m = UF.phi(sf_params.ufp, zeta, UF.MomentumTransport())
+    l_W = l_W_base / max(phi_m, eps_FT)
 
-    else # Neutral or stable case
-        # Ensure L > 0 for Monin-Obukhov length
-        obukhov_len_safe_stable = max(obukhov_length, eps_FT)
-        zeta = l_z / obukhov_len_safe_stable # zeta >= 0
-
-        # Stable/neutral-regime correction after Gryanik (2020):
-        #     φ_m = 1 + a_m ζ / (1 + b_m ζ)^(2/3),
-        # a nonlinear refinement to the Businger formulation.
-        # 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 b_m > 0, retained for robustness).
-        phi_m_denom = max(cbrt((1 + gryanik_b_m * zeta)^2), eps_FT)
-
-        phi_m = 1 + (gryanik_a_m * zeta) / phi_m_denom
-
-        # Stable-regime correction factor: 1 / φ_m.
-        # phi_m should be >= 1 for stable/neutral
-        stability_correction = 1 / max(phi_m, eps_FT)
-
-        # Apply the correction factor
-        l_W = l_W_base * stability_correction
-    end
     l_W = max(l_W, FT(0)) # Ensure non-negative
 
     # --- l_TKE: TKE production-dissipation balance scale ---