Skip to content

use the correct a_m in mixing length #4135

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
szy21 merged 1 commit into from
Dec 4, 2025
Merged

use the correct a_m in mixing length #4135

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
6 changes: 4 additions & 2 deletions src/parameters/Parameters.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand Down Expand Up @@ -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
Expand Down
4 changes: 3 additions & 1 deletion src/parameters/create_parameters.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down Expand Up @@ -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
Expand Down
15 changes: 7 additions & 8 deletions src/prognostic_equations/eddy_diffusion_closures.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down
Loading