From 92e3a7c1865e3a1754c485ef2c32ad089f4a7a7d Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Thu, 16 Apr 2026 13:34:51 +0200 Subject: [PATCH 01/11] Add coefficient-based bulk flux formulation (Large & Yeager 2004) Introduce `CoefficientBasedFluxes` with a unified `transfer_coefficients` interface that accepts either `SimilarityScales` (constant or callable coefficients) or `LargeYeagerTransferCoefficients` (full L&Y 2004/2009 stability-corrected algorithm). Add `PolynomialNeutralDragCoefficient`, `LinearStableStabilityFunction`, and `large_yeager_stability_functions`. Co-Authored-By: Claude Opus 4.6 (1M context) --- Project.toml | 4 + docs/src/earth_system_model.md | 2 +- docs/src/interface_fluxes.md | 64 +++- .../InterfaceComputations.jl | 5 + .../coefficient_based_turbulent_fluxes.jl | 297 +++++++++++++++--- .../similarity_theory_turbulent_fluxes.jl | 46 +++ src/NumericalEarth.jl | 1 + test/test_coefficient_based_fluxes.jl | 179 +++++++++++ 8 files changed, 535 insertions(+), 63 deletions(-) create mode 100644 test/test_coefficient_based_fluxes.jl diff --git a/Project.toml b/Project.toml index 0f01f8982..cc4a9fd63 100644 --- a/Project.toml +++ b/Project.toml @@ -11,6 +11,8 @@ ClimaSeaIce = "6ba0ff68-24e6-4315-936c-2e99227c95a4" DataDeps = "124859b0-ceae-595e-8997-d05f6a7a8dfe" Dates = "ade2ca70-3891-5945-98fb-dc099432e06a" DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae" +Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" +DocumenterCitations = "daee34ce-89f3-4625-b898-19384cb65244" Downloads = "f43a241f-c20a-4ad4-852c-f6b1247861c6" GPUArraysCore = "46192b85-c4d5-4398-a991-12ede77f4527" Glob = "c27321d9-0574-5035-807b-f59d2c89b15c" @@ -61,6 +63,8 @@ CopernicusMarine = "0.1.1" DataDeps = "0.7" Dates = "<0.0.1, 1" DocStringExtensions = "0.9" +Documenter = "1.17.0" +DocumenterCitations = "1.4.1" Downloads = "<0.0.1, 1" GPUArraysCore = "0.2.0" Glob = "1" diff --git a/docs/src/earth_system_model.md b/docs/src/earth_system_model.md index e405437b1..513aeabce 100644 --- a/docs/src/earth_system_model.md +++ b/docs/src/earth_system_model.md @@ -120,7 +120,7 @@ To change the defaults, construct `ComponentInterfaces` yourself and pass it in. For example, to use constant transfer coefficients instead of similarity theory: ```jldoctest esm -atmosphere_ocean_fluxes = CoefficientBasedFluxes(drag_coefficient=2e-3) +atmosphere_ocean_fluxes = CoefficientBasedFluxes(transfer_coefficients = SimilarityScales(2e-3, 2e-3, 2e-3)) interfaces = NumericalEarth.EarthSystemModels.ComponentInterfaces(nothing, ocean; atmosphere_ocean_fluxes) model = OceanOnlyModel(ocean; interfaces) diff --git a/docs/src/interface_fluxes.md b/docs/src/interface_fluxes.md index b776bc475..3613b1999 100644 --- a/docs/src/interface_fluxes.md +++ b/docs/src/interface_fluxes.md @@ -137,9 +137,29 @@ A comprehensive example is given below, but we note briefly here that ```@example interface_fluxes using NumericalEarth -coefficient_fluxes = CoefficientBasedFluxes(drag_coefficient=2e-3, - heat_transfer_coefficient=2e-3, - vapor_flux_coefficient=1e-3) +coefficient_fluxes = CoefficientBasedFluxes(transfer_coefficients = SimilarityScales(2e-3, 2e-3, 1e-3)) +``` + +Alternatively, the drag coefficient can be specified as a wind-speed-dependent polynomial +following Large & Yeager (2004). In this case `CoefficientBasedFluxes` evaluates the polynomial +at each iteration rather than using a constant: + +```@example interface_fluxes +poly_drag = PolynomialNeutralDragCoefficient() +poly_fluxes = CoefficientBasedFluxes(transfer_coefficients = SimilarityScales(poly_drag, 1e-3, 1e-3)) +``` + +For the full Large & Yeager (2004) bulk algorithm with stability corrections, +use `LargeYeagerTransferCoefficients`. This computes all three transfer coefficients +(drag, Stanton, Dalton) from the neutral drag polynomial with Monin-Obukhov +stability corrections (L&Y eqs. 6c-6d, 10a-10c): + +```@example interface_fluxes +using NumericalEarth.EarthSystemModels.InterfaceComputations: FixedIterations + +ly = LargeYeagerTransferCoefficients() +ly_fluxes = CoefficientBasedFluxes(transfer_coefficients = ly, + solver_stop_criteria = FixedIterations(5)) ``` ### Similarity theory for neutral boundary layers @@ -335,13 +355,18 @@ neutral_similarity_fluxes = SimilarityTheoryFluxes(stability_functions=nothing; interfaces = ComponentInterfaces(atmosphere, ocean; atmosphere_ocean_fluxes=neutral_similarity_fluxes) increased_roughness_model = OceanOnlyModel(ocean; atmosphere, interfaces) -coefficient_fluxes = CoefficientBasedFluxes(drag_coefficient=2e-3) +coefficient_fluxes = CoefficientBasedFluxes(transfer_coefficients = SimilarityScales(2e-3, 2e-3, 2e-3)) interfaces = ComponentInterfaces(atmosphere, ocean; atmosphere_ocean_fluxes=coefficient_fluxes) coefficient_model = OceanOnlyModel(ocean; atmosphere, interfaces) + +ly_fluxes = CoefficientBasedFluxes(transfer_coefficients = LargeYeagerTransferCoefficients(), + solver_stop_criteria = FixedIterations(5)) +interfaces = ComponentInterfaces(atmosphere, ocean; atmosphere_ocean_fluxes=ly_fluxes) +ly_model = OceanOnlyModel(ocean; atmosphere, interfaces) ``` Note that `EarthSystemModel` computes fluxes upon instantiation, so after constructing -the two models we are ready to analyze the results. +the models we are ready to analyze the results. We first verify that the similarity model friction velocity has been computed successfully, ```@example interface_fluxes @@ -370,7 +395,16 @@ Cᴰ_coeff = @. (u★_coeff / uᵃᵗ)^2 extrema(Cᴰ_coeff) ``` -We'll compare the computed fluxes and drag coefficients from our two models with +We also extract the drag coefficient from the Large & Yeager transfer coefficient model: + +```@example interface_fluxes +u★_ly = ly_model.interfaces.atmosphere_ocean_interface.fluxes.friction_velocity +u★_ly = interior(u★_ly, :, 1, 1) +Cᴰ_ly = @. (u★_ly / uᵃᵗ)^2 +extrema(Cᴰ_ly) +``` + +We'll compare the computed fluxes and drag coefficients from our models with a polynomial expression due to [large2009global](@citet), and an expression reported by [edson2013exchange](@citet) that was developed at ECMWF, @@ -406,16 +440,18 @@ Cᴰ_rough = @. (u★_rough / uᵃᵗ)^2 fig = Figure(size=(800, 400)) axu = Axis(fig[1:2, 1], xlabel="uᵃᵗ (m s⁻¹) at 10 m", ylabel="u★ (m s⁻¹)") -lines!(axu, uᵃᵗ, u★, label="NumericalEarth default") -lines!(axu, uᵃᵗ, u★_rough, label="Increased roughness model") -lines!(axu, uᵃᵗ, u★_LY, label="Large and Yeager (2009) polynomial fit") -lines!(axu, uᵃᵗ, u★_EC, label="ECMWF polynomial fit (Edson et al. 2013)") +lines!(axu, uᵃᵗ, u★, label="SimilarityTheoryFluxes (default)") +lines!(axu, uᵃᵗ, u★_rough, label="SimilarityTheoryFluxes (increased roughness)") +lines!(axu, uᵃᵗ, u★_ly, label="LargeYeagerTransferCoefficients (L&Y 2004)") +lines!(axu, uᵃᵗ, u★_LY, label="Large and Yeager (2009) polynomial fit", linestyle=:dash) +lines!(axu, uᵃᵗ, u★_EC, label="ECMWF polynomial fit (Edson et al. 2013)", linestyle=:dash) axd = Axis(fig[1:2, 2], xlabel="uᵃᵗ (m s⁻¹) at 10 m", ylabel="1000 × Cᴰ") -lines!(axd, uᵃᵗ, 1000 .* Cᴰ_default, label="NumericalEarth default") -lines!(axd, uᵃᵗ, 1000 .* Cᴰ_rough, label="Increased roughness model") -lines!(axd, uᵃᵗ, 1000 .* Cᴰ_LY, label="Large and Yeager (2009) polynomial fit") -lines!(axd, uᵃᵗ, 1000 .* Cᴰ_EC, label="ECMWF polynomial fit (Edson et al. 2013)") +lines!(axd, uᵃᵗ, 1000 .* Cᴰ_default, label="SimilarityTheoryFluxes (default)") +lines!(axd, uᵃᵗ, 1000 .* Cᴰ_rough, label="SimilarityTheoryFluxes (increased roughness)") +lines!(axd, uᵃᵗ, 1000 .* Cᴰ_ly, label="LargeYeagerTransferCoefficients (L&Y 2004)") +lines!(axd, uᵃᵗ, 1000 .* Cᴰ_LY, label="Large and Yeager (2009) polynomial fit", linestyle=:dash) +lines!(axd, uᵃᵗ, 1000 .* Cᴰ_EC, label="ECMWF polynomial fit (Edson et al. 2013)", linestyle=:dash) Legend(fig[3, 1:2], axd, nbanks = 2) diff --git a/src/EarthSystemModels/InterfaceComputations/InterfaceComputations.jl b/src/EarthSystemModels/InterfaceComputations/InterfaceComputations.jl index ed93f8c9c..5c28ec9b5 100644 --- a/src/EarthSystemModels/InterfaceComputations/InterfaceComputations.jl +++ b/src/EarthSystemModels/InterfaceComputations/InterfaceComputations.jl @@ -13,10 +13,15 @@ export MomentumRoughnessLength, ScalarRoughnessLength, CoefficientBasedFluxes, + SimilarityScales, + PolynomialNeutralDragCoefficient, + LargeYeagerTransferCoefficients, + LinearStableStabilityFunction, SkinTemperature, BulkTemperature, atmosphere_ocean_stability_functions, atmosphere_sea_ice_stability_functions, + large_yeager_stability_functions, compute_atmosphere_ocean_fluxes!, compute_atmosphere_sea_ice_fluxes!, compute_sea_ice_ocean_fluxes!, diff --git a/src/EarthSystemModels/InterfaceComputations/coefficient_based_turbulent_fluxes.jl b/src/EarthSystemModels/InterfaceComputations/coefficient_based_turbulent_fluxes.jl index b192e9baa..9dbef8976 100644 --- a/src/EarthSystemModels/InterfaceComputations/coefficient_based_turbulent_fluxes.jl +++ b/src/EarthSystemModels/InterfaceComputations/coefficient_based_turbulent_fluxes.jl @@ -1,61 +1,197 @@ using DocStringExtensions +##### +##### Polynomial neutral drag coefficient (Large & Yeager 2004) +##### + +""" + PolynomialNeutralDragCoefficient{FT} + +Neutral 10-m drag coefficient from the Large & Yeager (2004) polynomial: + + C_dN = (a/U + b + c*U - d*U^6) × 10⁻³ for U < high_wind_speed_threshold + C_dN = high_wind_drag_coefficient otherwise + +Called as a functor: `C_dN = drag(U)`. + +References: +- Large, W.G. & Yeager, S.G. (2004): NCAR/TN-460+STR +""" +struct PolynomialNeutralDragCoefficient{FT} + a :: FT + b :: FT + c :: FT + d :: FT + high_wind_speed_threshold :: FT + high_wind_drag_coefficient :: FT + minimum_wind_speed :: FT +end + +function PolynomialNeutralDragCoefficient(FT = Oceananigans.defaults.FloatType; + a = 2.7, + b = 0.142, + c = 1 / 13.09, + d = 3.14807e-10, + high_wind_speed_threshold = 33, + high_wind_drag_coefficient = 2.34e-3, + minimum_wind_speed = 0.5) + + return PolynomialNeutralDragCoefficient(convert(FT, a), + convert(FT, b), + convert(FT, c), + convert(FT, d), + convert(FT, high_wind_speed_threshold), + convert(FT, high_wind_drag_coefficient), + convert(FT, minimum_wind_speed)) +end + +@inline function (p::PolynomialNeutralDragCoefficient)(U) + U = max(U, p.minimum_wind_speed) + Cd = ifelse(U < p.high_wind_speed_threshold, + (p.a / U + p.b + p.c * U - p.d * U^6) / 1000, + p.high_wind_drag_coefficient) + return Cd +end + +Base.summary(::PolynomialNeutralDragCoefficient{FT}) where FT = + "PolynomialNeutralDragCoefficient{$FT}" + +##### +##### Large & Yeager transfer coefficients (L&Y 2004 with stability corrections) +##### + +""" + LargeYeagerTransferCoefficients{FT, D, SF} + +Wind-dependent bulk transfer coefficients following the NCAR/Large & Yeager (2004, 2009) +algorithm as used in OMIP-2 simulations. + +Computes all three transfer coefficients (drag, Stanton, Dalton) from a neutral drag +polynomial with stability corrections via Monin-Obukhov similarity functions. + +Pass as `transfer_coefficients` to `CoefficientBasedFluxes`: + + ly = LargeYeagerTransferCoefficients() + fluxes = CoefficientBasedFluxes(transfer_coefficients = ly, + solver_stop_criteria = FixedIterations(5)) + +References: +- Large, W.G. & Yeager, S.G. (2004): NCAR/TN-460+STR +- Large, W.G. & Yeager, S.G. (2009): Climate Dynamics, 33, 341-364 +""" +struct LargeYeagerTransferCoefficients{FT, D, SF} + von_karman_constant :: FT + neutral_drag_coefficient :: D + stability_functions :: SF + reference_height :: FT + stanton_stable :: FT + stanton_unstable :: FT + dalton :: FT +end + +function LargeYeagerTransferCoefficients(FT = Oceananigans.defaults.FloatType; + von_karman_constant = 0.4, + neutral_drag_coefficient = PolynomialNeutralDragCoefficient(FT), + stability_functions = large_yeager_stability_functions(FT), + reference_height = 10, + stanton_stable = 18, + stanton_unstable = 32.7, + dalton = 34.6) + + return LargeYeagerTransferCoefficients(convert(FT, von_karman_constant), + neutral_drag_coefficient, + stability_functions, + convert(FT, reference_height), + convert(FT, stanton_stable), + convert(FT, stanton_unstable), + convert(FT, dalton)) +end + +Adapt.adapt_structure(to, n::LargeYeagerTransferCoefficients) = + LargeYeagerTransferCoefficients(Adapt.adapt(to, n.von_karman_constant), + Adapt.adapt(to, n.neutral_drag_coefficient), + Adapt.adapt(to, n.stability_functions), + Adapt.adapt(to, n.reference_height), + Adapt.adapt(to, n.stanton_stable), + Adapt.adapt(to, n.stanton_unstable), + Adapt.adapt(to, n.dalton)) + +Base.summary(::LargeYeagerTransferCoefficients{FT}) where FT = "LargeYeagerTransferCoefficients{$FT}" + +function Base.show(io::IO, n::LargeYeagerTransferCoefficients{FT}) where FT + print(io, "LargeYeagerTransferCoefficients{$FT}", '\n') + print(io, "├── von_karman_constant: ", n.von_karman_constant, '\n') + print(io, "├── neutral_drag_coefficient: ", summary(n.neutral_drag_coefficient), '\n') + print(io, "├── stability_functions: ", summary(n.stability_functions), '\n') + print(io, "├── reference_height: ", n.reference_height, '\n') + print(io, "├── stanton (stable/unstable): ", n.stanton_stable, " / ", n.stanton_unstable, '\n') + print(io, "└── dalton: ", n.dalton) +end + +##### +##### CoefficientBasedFluxes +##### + """ - struct CoefficientBasedFluxes{CD, CH, CQ, S} + struct CoefficientBasedFluxes{C, S} -A structure for computing turbulent fluxes using constant bulk transfer coefficients. +A structure for computing turbulent fluxes using bulk transfer coefficients. $(TYPEDFIELDS) """ -struct CoefficientBasedFluxes{CD, CH, CQ, S} - "Coefficient for momentum transfer" - drag_coefficient :: CD - "Coefficient for sensible heat transfer" - heat_transfer_coefficient :: CH - "Coefficient for latent heat transfer" - vapor_flux_coefficient :: CQ - "Criteria for iterative solver convergence" - solver_stop_criteria :: S +struct CoefficientBasedFluxes{C, S} + transfer_coefficients :: C # `SimilarityScales` with constant or callable coefficients, or an `LargeYeagerTransferCoefficients`." + solver_stop_criteria :: S # "Criteria for iterative solver convergence." end +Adapt.adapt_structure(to, f::CoefficientBasedFluxes) = CoefficientBasedFluxes(Adapt.adapt(to, f.transfer_coefficients), f.solver_stop_criteria) + Base.summary(flux_formulation::CoefficientBasedFluxes) = "CoefficientBasedFluxes" function Base.show(io::IO, flux_formulation::CoefficientBasedFluxes) print(io, summary(flux_formulation), '\n') - print(io, "├── drag_coefficient: ", prettysummary(flux_formulation.drag_coefficient), '\n') - print(io, "├── heat_transfer_coefficient: ", prettysummary(flux_formulation.heat_transfer_coefficient), '\n') - print(io, "└── vapor_flux_coefficient: ", prettysummary(flux_formulation.vapor_flux_coefficient), '\n') + print(io, "├── transfer_coefficients: ", summary(flux_formulation.transfer_coefficients), '\n') + print(io, "└── solver_stop_criteria: ", summary(flux_formulation.solver_stop_criteria)) end convert_if_number(FT, a::Number) = convert(FT, a) convert_if_number(FT, a) = a +convert_transfer_coefficients(FT, c) = c +convert_transfer_coefficients(FT, c::SimilarityScales) = SimilarityScales(convert_if_number(FT, c.momentum), + convert_if_number(FT, c.temperature), + convert_if_number(FT, c.water_vapor)) + """ CoefficientBasedFluxes(FT = Oceananigans.defaults.FloatType; - drag_coefficient = 1e-3, - heat_transfer_coefficient = drag_coefficient, - vapor_flux_coefficient = drag_coefficient, - solver_stop_criteria = nothing, - solver_tolerance = 1e-8, - solver_maxiter = 20) - -Return the structure for computing turbulent fluxes using constant bulk transfer coefficients. + transfer_coefficients = SimilarityScales(1e-3, 1e-3, 1e-3), + solver_stop_criteria = nothing, + solver_tolerance = 1e-8, + solver_maxiter = 20) + +Return the structure for computing turbulent fluxes using bulk transfer coefficients. Used in bulk flux calculations to determine the exchange of momentum, heat, and moisture -between the ocean/ice surface and the atmosphere using constant transfer coefficients. +between the ocean/ice surface and the atmosphere. + +Arguments +========= -# Arguments - `FT`: (optional) Float type for the coefficients, defaults to Oceananigans.defaults.FloatType -# Keyword Arguments -- `drag_coefficient`: Coefficient for momentum transfer (`Cᵈ`), defaults to 1e-3. -- `heat_transfer_coefficient`: Coefficient for heat transfer (`Cʰ`), defaults to `drag_coefficient` -- `vapor_flux_coefficient`: Coefficient for moisture transfer (`Cᵍ`), defaults to `drag_coefficient` +Keyword Arguments +================= + +- `transfer_coefficients`: Transfer coefficients for momentum, heat, and moisture. + Can be a `SimilarityScales` with constant or callable entries, or an `LargeYeagerTransferCoefficients`. + Defaults to `SimilarityScales(1e-3, 1e-3, 1e-3)`. - `solver_stop_criteria`: Criteria for iterative solver convergence. If `nothing`, creates new criteria using `solver_tolerance` and `solver_maxiter`. - `solver_tolerance`: Tolerance for solver convergence when creating new stop criteria, defaults to 1e-8. - `solver_maxiter`: Maximum iterations for solver when creating new stop criteria, defaults to 20 -# Example +Example +======== + ```jldoctest using Oceananigans using NumericalEarth @@ -63,9 +199,7 @@ using NumericalEarth grid = RectilinearGrid(size=3, z=(-1, 0), topology=(Flat, Flat, Bounded)) ocean = ocean_simulation(grid; timestepper = :QuasiAdamsBashforth2) -ao_fluxes = CoefficientBasedFluxes(drag_coefficient = 1e-2, - heat_transfer_coefficient = 1e-3, - vapor_flux_coefficient = 1e-3) +ao_fluxes = CoefficientBasedFluxes(transfer_coefficients = SimilarityScales(1e-2, 1e-3, 1e-3)) interfaces = ComponentInterfaces(nothing, ocean; atmosphere_ocean_fluxes=ao_fluxes) @@ -74,9 +208,7 @@ ComponentInterfaces ``` """ function CoefficientBasedFluxes(FT = Oceananigans.defaults.FloatType; - drag_coefficient = 1e-3, - heat_transfer_coefficient = drag_coefficient, - vapor_flux_coefficient = drag_coefficient, + transfer_coefficients = SimilarityScales(1e-3, 1e-3, 1e-3), solver_stop_criteria = nothing, solver_tolerance = 1e-8, solver_maxiter = 20) @@ -86,16 +218,84 @@ function CoefficientBasedFluxes(FT = Oceananigans.defaults.FloatType; solver_stop_criteria = ConvergenceStopCriteria(solver_tolerance, solver_maxiter) end - drag_coefficient = convert_if_number(FT, drag_coefficient) - heat_transfer_coefficient = convert_if_number(FT, heat_transfer_coefficient) - vapor_flux_coefficient = convert_if_number(FT, vapor_flux_coefficient) + transfer_coefficients = convert_transfer_coefficients(FT, transfer_coefficients) - return CoefficientBasedFluxes(drag_coefficient, - heat_transfer_coefficient, - vapor_flux_coefficient, - solver_stop_criteria) + return CoefficientBasedFluxes(transfer_coefficients, solver_stop_criteria) end +##### +##### Evaluate transfer coefficients (dispatch on coefficient type) +##### + +# Constant coefficients: no stability computation needed +@inline evaluate_coefficients(coeffs::SimilarityScales{<:Number, <:Number, <:Number}, args...) = (coeffs.momentum, coeffs.temperature, coeffs.water_vapor) + +# Polynomial drag coefficient (no stability correction): wind-dependent Cd, constant Ch and Cq +@inline function evaluate_coefficients(coeffs::SimilarityScales{<:PolynomialNeutralDragCoefficient, <:Number, <:Number}, ΔU, args...) + Cd = coeffs.momentum(ΔU) + return Cd, coeffs.temperature, coeffs.water_vapor +end + +# NCAR transfer coefficients: full L&Y stability correction (eqs. 6c-6d, 10a-10c) +@inline function evaluate_coefficients(ly::LargeYeagerTransferCoefficients, + ΔU, Tₛ, qₛ, Δh, + approximate_interface_state, + atmosphere_properties) + + FT = eltype(approximate_interface_state) + ℂᵃᵗ = atmosphere_properties.thermodynamics_parameters + g = atmosphere_properties.gravitational_acceleration + + κ = ly.von_karman_constant + u★ = approximate_interface_state.u★ + θ★ = approximate_interface_state.θ★ + q★ = approximate_interface_state.q★ + + neutral_drag = ly.neutral_drag_coefficient + h₀ = ly.reference_height + Umin = neutral_drag.minimum_wind_speed + ΔU = max(ΔU, Umin) + + # Monin-Obukhov length from previous iteration scales + b★ = buoyancy_scale(θ★, q★, ℂᵃᵗ, Tₛ, qₛ, g) + L★ = ifelse(b★ == 0, convert(FT, Inf), u★^2 / (κ * b★)) + ζ = Δh / L★ + + # Stability functions + ψₘ = stability_profile(ly.stability_functions.momentum, ζ) + ψₕ = stability_profile(ly.stability_functions.temperature, ζ) + + # Neutral 10-m wind speed + Cdp = ifelse(u★ == 0, neutral_drag(ΔU), u★^2 / ΔU^2) + UN10 = ΔU / (1 + sqrt(Cdp) / κ * (log(Δh / h₀) - ψₘ)) + UN10 = max(UN10, Umin) + + # Neutral transfer coefficients + CdN = neutral_drag(UN10) + sqrtCdN = sqrt(CdN) + + # Stability-dependent neutral Stanton/Dalton numbers (L&Y eq. 6c-6d) + stable = ζ > 0 + ChN = sqrtCdN / 1000 * ifelse(stable, ly.stanton_stable, ly.stanton_unstable) + CqN = sqrtCdN / 1000 * ly.dalton + + # Stability corrections (L&Y eq. 10a-10c) + ξₘ = sqrtCdN / κ * (log(Δh / h₀) - ψₘ) + Cd = CdN / (1 + ξₘ)^2 + + ξₕ = sqrtCdN / κ * (log(Δh / h₀) - ψₕ) + ratio = sqrt(Cd) / sqrtCdN + + Ch = ChN * ratio / (1 + ChN * ξₕ) + Cq = CqN * ratio / (1 + CqN * ξₕ) + + return Cd, Ch, Cq +end + +##### +##### Iteration +##### + @inline function iterate_interface_fluxes(flux_formulation::CoefficientBasedFluxes, Tₛ, qₛ, Δθ, Δq, Δh, approximate_interface_state, @@ -103,14 +303,15 @@ end interface_properties, atmosphere_properties) - Ψₐ = atmosphere_state - Ψ̃ᵢ = approximate_interface_state - Δu, Δv = velocity_difference(interface_properties.velocity_formulation, Ψₐ, Ψ̃ᵢ) + Δu, Δv = velocity_difference(interface_properties.velocity_formulation, + atmosphere_state, + approximate_interface_state) ΔU = sqrt(Δu^2 + Δv^2) - Cd = flux_formulation.drag_coefficient - Ch = flux_formulation.heat_transfer_coefficient - Cq = flux_formulation.vapor_flux_coefficient + Cd, Ch, Cq = evaluate_coefficients(flux_formulation.transfer_coefficients, + ΔU, Tₛ, qₛ, Δh, + approximate_interface_state, + atmosphere_properties) u★ = sqrt(Cd) * ΔU θ★ = Ch / sqrt(Cd) * Δθ diff --git a/src/EarthSystemModels/InterfaceComputations/similarity_theory_turbulent_fluxes.jl b/src/EarthSystemModels/InterfaceComputations/similarity_theory_turbulent_fluxes.jl index dc02ebfd9..cdc27fe89 100644 --- a/src/EarthSystemModels/InterfaceComputations/similarity_theory_turbulent_fluxes.jl +++ b/src/EarthSystemModels/InterfaceComputations/similarity_theory_turbulent_fluxes.jl @@ -568,6 +568,52 @@ Base.show(io::IO, ss::SplitStabilityFunction) = print(io, "SplitStabilityFunctio return ifelse(stable, Ψ_stable, Ψ_unstable) end +##### +##### Linear stable stability function (ψ = -c ζ, bounded) +##### + +""" + LinearStableStabilityFunction{FT} + +A simple linear stability function for stable conditions: ``ψ = -c ζ``, +bounded at ``|ζ| ≤ ζ_{max}``. + +Used by the NCAR/Large-Yeager (2004) bulk formulae with ``c = 5`` and ``ζ_{max} = 10``. + +References: +- Large, W.G. & Yeager, S.G. (2004): NCAR/TN-460+STR +""" +@kwdef struct LinearStableStabilityFunction{FT} <: AbstractStabilityFunction + coefficient :: FT = 5.0 + maximum_stability_parameter :: FT = 10.0 +end + +@inline function stability_profile(ψ::LinearStableStabilityFunction, ζ) + c = ψ.coefficient + ζmax = ψ.maximum_stability_parameter + ζ⁺ = max(zero(ζ), ζ) + return -c * min(ζ⁺, ζmax) +end + +Base.summary(::LinearStableStabilityFunction{FT}) where FT = "LinearStableStabilityFunction{$FT}" +Base.show(io::IO, ::LinearStableStabilityFunction{FT}) where FT = print(io, "LinearStableStabilityFunction{$FT}") + +""" + large_yeager_stability_functions(FT = Float64) + +NCAR/Large-Yeager (2004) stability functions combining: +- Unstable: Paulson (1970) with γ = 16 +- Stable: linear ψ = -5ζ, bounded at |ζ| ≤ 10 + +Used for OMIP-2 protocol compliance. +""" +function large_yeager_stability_functions(FT=Oceananigans.defaults.FloatType) + stable = LinearStableStabilityFunction{FT}() + momentum = SplitStabilityFunction(stable, PaulsonMomentumStabilityFunction{FT}()) + scalar = SplitStabilityFunction(stable, PaulsonScalarStabilityFunction{FT}()) + return SimilarityScales(momentum, scalar, scalar) +end + function atmosphere_sea_ice_stability_functions(FT=Oceananigans.defaults.FloatType) unstable_momentum = PaulsonMomentumStabilityFunction{FT}() stable_momentum = ShebaMomentumStabilityFunction{FT}() diff --git a/src/NumericalEarth.jl b/src/NumericalEarth.jl index 1421eadc5..3cdea607c 100644 --- a/src/NumericalEarth.jl +++ b/src/NumericalEarth.jl @@ -19,6 +19,7 @@ export LatitudeDependentAlbedo, SimilarityTheoryFluxes, CoefficientBasedFluxes, + SimilarityScales, MomentumRoughnessLength, ScalarRoughnessLength, ComponentInterfaces, diff --git a/test/test_coefficient_based_fluxes.jl b/test/test_coefficient_based_fluxes.jl new file mode 100644 index 000000000..00b6a2e2f --- /dev/null +++ b/test/test_coefficient_based_fluxes.jl @@ -0,0 +1,179 @@ +include("runtests_setup.jl") + +using NumericalEarth.EarthSystemModels.InterfaceComputations: + CoefficientBasedFluxes, + PolynomialNeutralDragCoefficient, + LargeYeagerTransferCoefficients, + LinearStableStabilityFunction, + large_yeager_stability_functions, + SimilarityScales, + ComponentInterfaces, + stability_profile, + FixedIterations + +using NumericalEarth.DataWrangling: all_dates + +@testset "PolynomialNeutralDragCoefficient" begin + p = PolynomialNeutralDragCoefficient() + @test p isa PolynomialNeutralDragCoefficient{Float64} + + # Basic evaluation + @test p(3.0) > 0 + @test p(3.0) < 5e-3 + + # High wind cap + @test p(40.0) ≈ 2.34e-3 + + # Wind floor + @test p(0.0) == p(0.5) + + # Monotonicity in moderate winds + @test p(20.0) > p(5.0) + + # Float32 + p32 = PolynomialNeutralDragCoefficient(Float32) + @test p32 isa PolynomialNeutralDragCoefficient{Float32} + @test p32(10f0) isa Float32 +end + +@testset "LinearStableStabilityFunction and large_yeager_stability_functions" begin + ψ = LinearStableStabilityFunction{Float64}() + @test stability_profile(ψ, 0.0) ≈ 0.0 + @test stability_profile(ψ, 1.0) ≈ -5.0 + @test stability_profile(ψ, -1.0) ≈ 0.0 + @test stability_profile(ψ, 20.0) ≈ -50.0 # bounded at ζ_max = 10 + + sf = large_yeager_stability_functions() + @test sf isa SimilarityScales + + # Momentum: Paulson unstable + linear stable + @test stability_profile(sf.momentum, -1.0) > 0 + @test stability_profile(sf.momentum, 1.0) ≈ -5.0 + @test stability_profile(sf.momentum, 0.0) ≈ 0.0 atol=1e-10 + + # Scalar: same structure + @test stability_profile(sf.temperature, -1.0) > 0 + @test stability_profile(sf.temperature, 1.0) ≈ -5.0 +end + +@testset "LargeYeagerTransferCoefficients constructor" begin + ly = LargeYeagerTransferCoefficients() + @test ly isa LargeYeagerTransferCoefficients{Float64} + @test ly.reference_height ≈ 10.0 + @test ly.stanton_stable ≈ 18.0 + @test ly.stanton_unstable ≈ 32.7 + @test ly.dalton ≈ 34.6 + @test ly.neutral_drag_coefficient isa PolynomialNeutralDragCoefficient{Float64} + + ly32 = LargeYeagerTransferCoefficients(Float32) + @test ly32 isa LargeYeagerTransferCoefficients{Float32} +end + +@testset "CoefficientBasedFluxes with constant coefficients" begin + for arch in test_architectures + grid = LatitudeLongitudeGrid(arch; + size = 1, + latitude = 10, + longitude = 10, + z = (-1, 0), + topology = (Flat, Flat, Bounded)) + + ocean = ocean_simulation(grid; + momentum_advection = nothing, + tracer_advection = nothing, + closure = nothing, + bottom_drag_coefficient = 0) + + dates = all_dates(RepeatYearJRA55(), :temperature) + atmosphere = JRA55PrescribedAtmosphere(arch, Float64; end_date=dates[2], backend=InMemory()) + + constant_fluxes = CoefficientBasedFluxes(transfer_coefficients = SimilarityScales(2e-3, 2e-3, 2e-3)) + interfaces = ComponentInterfaces(atmosphere, ocean; + atmosphere_ocean_fluxes=constant_fluxes) + + set!(ocean.model, T=15, S=35) + coupled_model = OceanOnlyModel(ocean; atmosphere, interfaces) + fluxes = coupled_model.interfaces.atmosphere_ocean_interface.fluxes + + CUDA.@allowscalar begin + @test isfinite(fluxes.sensible_heat[1, 1, 1]) + @test isfinite(fluxes.latent_heat[1, 1, 1]) + @test isfinite(fluxes.water_vapor[1, 1, 1]) + end + end +end + +@testset "CoefficientBasedFluxes with PolynomialNeutralDragCoefficient" begin + for arch in test_architectures + grid = LatitudeLongitudeGrid(arch; + size = 1, + latitude = 10, + longitude = 10, + z = (-1, 0), + topology = (Flat, Flat, Bounded)) + + ocean = ocean_simulation(grid; + momentum_advection = nothing, + tracer_advection = nothing, + closure = nothing, + bottom_drag_coefficient = 0) + + dates = all_dates(RepeatYearJRA55(), :temperature) + atmosphere = JRA55PrescribedAtmosphere(arch, Float64; end_date=dates[2], backend=InMemory()) + + poly_drag = PolynomialNeutralDragCoefficient() + poly_fluxes = CoefficientBasedFluxes(transfer_coefficients = SimilarityScales(poly_drag, 1e-3, 1e-3)) + + interfaces = ComponentInterfaces(atmosphere, ocean; + atmosphere_ocean_fluxes=poly_fluxes) + + set!(ocean.model, T=15, S=35) + coupled_model = OceanOnlyModel(ocean; atmosphere, interfaces) + fluxes = coupled_model.interfaces.atmosphere_ocean_interface.fluxes + + CUDA.@allowscalar begin + @test isfinite(fluxes.sensible_heat[1, 1, 1]) + @test isfinite(fluxes.latent_heat[1, 1, 1]) + @test isfinite(fluxes.water_vapor[1, 1, 1]) + @test fluxes.friction_velocity[1, 1, 1] > 0 + end + end +end + +@testset "CoefficientBasedFluxes with LargeYeagerTransferCoefficients" begin + for arch in test_architectures + grid = LatitudeLongitudeGrid(arch; + size = 1, + latitude = 10, + longitude = 10, + z = (-1, 0), + topology = (Flat, Flat, Bounded)) + + ocean = ocean_simulation(grid; + momentum_advection = nothing, + tracer_advection = nothing, + closure = nothing, + bottom_drag_coefficient = 0) + + dates = all_dates(RepeatYearJRA55(), :temperature) + atmosphere = JRA55PrescribedAtmosphere(arch, Float64; end_date=dates[2], backend=InMemory()) + + ly = LargeYeagerTransferCoefficients() + ly_fluxes = CoefficientBasedFluxes(transfer_coefficients = ly, + solver_stop_criteria = FixedIterations(5)) + + interfaces = ComponentInterfaces(atmosphere, ocean; + atmosphere_ocean_fluxes=ly_fluxes) + + set!(ocean.model, T=15, S=35) + coupled_model = OceanOnlyModel(ocean; atmosphere, interfaces) + fluxes = coupled_model.interfaces.atmosphere_ocean_interface.fluxes + + CUDA.@allowscalar begin + @test isfinite(fluxes.sensible_heat[1, 1, 1]) + @test isfinite(fluxes.latent_heat[1, 1, 1]) + @test isfinite(fluxes.water_vapor[1, 1, 1]) + @test fluxes.friction_velocity[1, 1, 1] > 0 + end + end +end From d5b495dfe9720036fa96a0202052fccf41d757d5 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Thu, 16 Apr 2026 13:43:18 +0200 Subject: [PATCH 02/11] a tad more generic --- .../coefficient_based_turbulent_fluxes.jl | 52 ++++++++++--------- 1 file changed, 27 insertions(+), 25 deletions(-) diff --git a/src/EarthSystemModels/InterfaceComputations/coefficient_based_turbulent_fluxes.jl b/src/EarthSystemModels/InterfaceComputations/coefficient_based_turbulent_fluxes.jl index 9dbef8976..ca5239157 100644 --- a/src/EarthSystemModels/InterfaceComputations/coefficient_based_turbulent_fluxes.jl +++ b/src/EarthSystemModels/InterfaceComputations/coefficient_based_turbulent_fluxes.jl @@ -90,31 +90,31 @@ struct LargeYeagerTransferCoefficients{FT, D, SF} end function LargeYeagerTransferCoefficients(FT = Oceananigans.defaults.FloatType; - von_karman_constant = 0.4, - neutral_drag_coefficient = PolynomialNeutralDragCoefficient(FT), - stability_functions = large_yeager_stability_functions(FT), - reference_height = 10, - stanton_stable = 18, - stanton_unstable = 32.7, - dalton = 34.6) + von_karman_constant = 0.4, + neutral_drag_coefficient = PolynomialNeutralDragCoefficient(FT), + stability_functions = large_yeager_stability_functions(FT), + reference_height = 10, + stanton_stable = 18, + stanton_unstable = 32.7, + dalton = 34.6) return LargeYeagerTransferCoefficients(convert(FT, von_karman_constant), - neutral_drag_coefficient, - stability_functions, - convert(FT, reference_height), - convert(FT, stanton_stable), - convert(FT, stanton_unstable), - convert(FT, dalton)) + neutral_drag_coefficient, + stability_functions, + convert(FT, reference_height), + convert(FT, stanton_stable), + convert(FT, stanton_unstable), + convert(FT, dalton)) end Adapt.adapt_structure(to, n::LargeYeagerTransferCoefficients) = LargeYeagerTransferCoefficients(Adapt.adapt(to, n.von_karman_constant), - Adapt.adapt(to, n.neutral_drag_coefficient), - Adapt.adapt(to, n.stability_functions), - Adapt.adapt(to, n.reference_height), - Adapt.adapt(to, n.stanton_stable), - Adapt.adapt(to, n.stanton_unstable), - Adapt.adapt(to, n.dalton)) + Adapt.adapt(to, n.neutral_drag_coefficient), + Adapt.adapt(to, n.stability_functions), + Adapt.adapt(to, n.reference_height), + Adapt.adapt(to, n.stanton_stable), + Adapt.adapt(to, n.stanton_unstable), + Adapt.adapt(to, n.dalton)) Base.summary(::LargeYeagerTransferCoefficients{FT}) where FT = "LargeYeagerTransferCoefficients{$FT}" @@ -227,13 +227,15 @@ end ##### Evaluate transfer coefficients (dispatch on coefficient type) ##### -# Constant coefficients: no stability computation needed -@inline evaluate_coefficients(coeffs::SimilarityScales{<:Number, <:Number, <:Number}, args...) = (coeffs.momentum, coeffs.temperature, coeffs.water_vapor) +@inline evaluate_coefficient(C::Number, args...) = C +@inline evaluate_coefficient(C::Function, args...) = C(args...) +@inline evaluate_coefficient(C::PolynomialNeutralDragCoefficient, ΔU, args...) = C(ΔU) -# Polynomial drag coefficient (no stability correction): wind-dependent Cd, constant Ch and Cq -@inline function evaluate_coefficients(coeffs::SimilarityScales{<:PolynomialNeutralDragCoefficient, <:Number, <:Number}, ΔU, args...) - Cd = coeffs.momentum(ΔU) - return Cd, coeffs.temperature, coeffs.water_vapor +@inline function evaluate_coefficients(coeffs::SimilarityScales, args...) + Cd = evaluate_coefficient(coeffs.momentum, args...) + Ch = evaluate_coefficient(coeffs.temperature, args...) + Cq = evaluate_coefficient(coeffs.water_vapor, arhs...) + return Cd, Ch, Cq end # NCAR transfer coefficients: full L&Y stability correction (eqs. 6c-6d, 10a-10c) From bd3cdd10902e24fca402e6836459c4a33cb5089b Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Thu, 16 Apr 2026 16:39:49 +0200 Subject: [PATCH 03/11] ensure a minimum speed --- .../coefficient_based_turbulent_fluxes.jl | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/src/EarthSystemModels/InterfaceComputations/coefficient_based_turbulent_fluxes.jl b/src/EarthSystemModels/InterfaceComputations/coefficient_based_turbulent_fluxes.jl index ca5239157..e72068b4d 100644 --- a/src/EarthSystemModels/InterfaceComputations/coefficient_based_turbulent_fluxes.jl +++ b/src/EarthSystemModels/InterfaceComputations/coefficient_based_turbulent_fluxes.jl @@ -227,6 +227,10 @@ end ##### Evaluate transfer coefficients (dispatch on coefficient type) ##### +# Minimum wind speed floor: only LargeYeager needs it; constant coefficients don't. +@inline minimum_wind_speed(::SimilarityScales) = 0 +@inline minimum_wind_speed(ly::LargeYeagerTransferCoefficients) = ly.neutral_drag_coefficient.minimum_wind_speed + @inline evaluate_coefficient(C::Number, args...) = C @inline evaluate_coefficient(C::Function, args...) = C(args...) @inline evaluate_coefficient(C::PolynomialNeutralDragCoefficient, ΔU, args...) = C(ΔU) @@ -308,7 +312,9 @@ end Δu, Δv = velocity_difference(interface_properties.velocity_formulation, atmosphere_state, approximate_interface_state) - ΔU = sqrt(Δu^2 + Δv^2) + + Umin = minimum_wind_speed(flux_formulation.transfer_coefficients) + ΔU = max(sqrt(Δu^2 + Δv^2), Umin) Cd, Ch, Cq = evaluate_coefficients(flux_formulation.transfer_coefficients, ΔU, Tₛ, qₛ, Δh, From a4461d836e81c9bc6c8328e7804d8a2826f0236c Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Thu, 16 Apr 2026 17:00:44 +0200 Subject: [PATCH 04/11] apply a guard based on the drag coefficient --- .../coefficient_based_turbulent_fluxes.jl | 15 +++++++-------- 1 file changed, 7 insertions(+), 8 deletions(-) diff --git a/src/EarthSystemModels/InterfaceComputations/coefficient_based_turbulent_fluxes.jl b/src/EarthSystemModels/InterfaceComputations/coefficient_based_turbulent_fluxes.jl index e72068b4d..b28857474 100644 --- a/src/EarthSystemModels/InterfaceComputations/coefficient_based_turbulent_fluxes.jl +++ b/src/EarthSystemModels/InterfaceComputations/coefficient_based_turbulent_fluxes.jl @@ -278,19 +278,18 @@ end # Neutral transfer coefficients CdN = neutral_drag(UN10) - sqrtCdN = sqrt(CdN) # Stability-dependent neutral Stanton/Dalton numbers (L&Y eq. 6c-6d) stable = ζ > 0 - ChN = sqrtCdN / 1000 * ifelse(stable, ly.stanton_stable, ly.stanton_unstable) - CqN = sqrtCdN / 1000 * ly.dalton + ChN = sqrt(CdN) / 1000 * ifelse(stable, ly.stanton_stable, ly.stanton_unstable) + CqN = sqrt(CdN) / 1000 * ly.dalton # Stability corrections (L&Y eq. 10a-10c) - ξₘ = sqrtCdN / κ * (log(Δh / h₀) - ψₘ) + ξₘ = sqrt(CdN) / κ * (log(Δh / h₀) - ψₘ) Cd = CdN / (1 + ξₘ)^2 - ξₕ = sqrtCdN / κ * (log(Δh / h₀) - ψₕ) - ratio = sqrt(Cd) / sqrtCdN + ξₕ = sqrt(CdN) / κ * (log(Δh / h₀) - ψₕ) + ratio = sqrt(Cd) / sqrt(CdN) Ch = ChN * ratio / (1 + ChN * ξₕ) Cq = CqN * ratio / (1 + CqN * ξₕ) @@ -322,8 +321,8 @@ end atmosphere_properties) u★ = sqrt(Cd) * ΔU - θ★ = Ch / sqrt(Cd) * Δθ - q★ = Cq / sqrt(Cd) * Δq + θ★ = ifelse(Cd == 0, zero(Δθ), Ch / sqrt(Cd) * Δθ) + q★ = ifelse(Cd == 0, zero(Δq), Cq / sqrt(Cd) * Δq) return u★, θ★, q★ end From c07de50f6b52b41233b85824c97bc803217be8cf Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Mon, 20 Apr 2026 19:10:51 +0200 Subject: [PATCH 05/11] typo arhs -> args --- Project.toml | 2 +- .../coefficient_based_turbulent_fluxes.jl | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/Project.toml b/Project.toml index 6aaea8e24..ed3ef1bd3 100644 --- a/Project.toml +++ b/Project.toml @@ -11,7 +11,7 @@ ClimaSeaIce = "6ba0ff68-24e6-4315-936c-2e99227c95a4" DataDeps = "124859b0-ceae-595e-8997-d05f6a7a8dfe" Dates = "ade2ca70-3891-5945-98fb-dc099432e06a" DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae" -Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" +arhsDocumenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" DocumenterCitations = "daee34ce-89f3-4625-b898-19384cb65244" Downloads = "f43a241f-c20a-4ad4-852c-f6b1247861c6" GPUArraysCore = "46192b85-c4d5-4398-a991-12ede77f4527" diff --git a/src/EarthSystemModels/InterfaceComputations/coefficient_based_turbulent_fluxes.jl b/src/EarthSystemModels/InterfaceComputations/coefficient_based_turbulent_fluxes.jl index b28857474..fbb8c5e29 100644 --- a/src/EarthSystemModels/InterfaceComputations/coefficient_based_turbulent_fluxes.jl +++ b/src/EarthSystemModels/InterfaceComputations/coefficient_based_turbulent_fluxes.jl @@ -236,9 +236,9 @@ end @inline evaluate_coefficient(C::PolynomialNeutralDragCoefficient, ΔU, args...) = C(ΔU) @inline function evaluate_coefficients(coeffs::SimilarityScales, args...) - Cd = evaluate_coefficient(coeffs.momentum, args...) + Cd = evaluate_coefficient(coeffs.momentum, args...) Ch = evaluate_coefficient(coeffs.temperature, args...) - Cq = evaluate_coefficient(coeffs.water_vapor, arhs...) + Cq = evaluate_coefficient(coeffs.water_vapor, args...) return Cd, Ch, Cq end From 321cdc5c5f4ee2cb21f1486b6844668372c34083 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Mon, 20 Apr 2026 19:11:21 +0200 Subject: [PATCH 06/11] Rename arhsDocumenter to Documenter in Project.toml --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index ed3ef1bd3..6aaea8e24 100644 --- a/Project.toml +++ b/Project.toml @@ -11,7 +11,7 @@ ClimaSeaIce = "6ba0ff68-24e6-4315-936c-2e99227c95a4" DataDeps = "124859b0-ceae-595e-8997-d05f6a7a8dfe" Dates = "ade2ca70-3891-5945-98fb-dc099432e06a" DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae" -arhsDocumenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" +Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" DocumenterCitations = "daee34ce-89f3-4625-b898-19384cb65244" Downloads = "f43a241f-c20a-4ad4-852c-f6b1247861c6" GPUArraysCore = "46192b85-c4d5-4398-a991-12ede77f4527" From 68abcf80ed61279426c8cd447503fc6bcb7916a4 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Mon, 20 Apr 2026 20:40:46 +0200 Subject: [PATCH 07/11] better with validation --- docs/src/earth_system_model.md | 2 +- .../coefficient_based_turbulent_fluxes.jl | 57 ++++++++++++++++--- 2 files changed, 50 insertions(+), 9 deletions(-) diff --git a/docs/src/earth_system_model.md b/docs/src/earth_system_model.md index 513aeabce..a88c4d839 100644 --- a/docs/src/earth_system_model.md +++ b/docs/src/earth_system_model.md @@ -120,7 +120,7 @@ To change the defaults, construct `ComponentInterfaces` yourself and pass it in. For example, to use constant transfer coefficients instead of similarity theory: ```jldoctest esm -atmosphere_ocean_fluxes = CoefficientBasedFluxes(transfer_coefficients = SimilarityScales(2e-3, 2e-3, 2e-3)) +atmosphere_ocean_fluxes = CoefficientBasedFluxes(transfer_coefficients = (2e-3, 2e-3, 2e-3)) interfaces = NumericalEarth.EarthSystemModels.ComponentInterfaces(nothing, ocean; atmosphere_ocean_fluxes) model = OceanOnlyModel(ocean; interfaces) diff --git a/src/EarthSystemModels/InterfaceComputations/coefficient_based_turbulent_fluxes.jl b/src/EarthSystemModels/InterfaceComputations/coefficient_based_turbulent_fluxes.jl index fbb8c5e29..37bed31f9 100644 --- a/src/EarthSystemModels/InterfaceComputations/coefficient_based_turbulent_fluxes.jl +++ b/src/EarthSystemModels/InterfaceComputations/coefficient_based_turbulent_fluxes.jl @@ -158,13 +158,18 @@ convert_if_number(FT, a::Number) = convert(FT, a) convert_if_number(FT, a) = a convert_transfer_coefficients(FT, c) = c +convert_transfer_coefficients(FT, c::Tuple) = convert_if_number.(FT, c) +convert_transfer_coefficients(FT, c::NamedTuple) = (; momentum = convert_if_number(FT, c.momentum) + temperature = convert_if_number(FT, c.temperature), + water_vapor = convert_if_number(FT, c.water_vapor)) + convert_transfer_coefficients(FT, c::SimilarityScales) = SimilarityScales(convert_if_number(FT, c.momentum), convert_if_number(FT, c.temperature), convert_if_number(FT, c.water_vapor)) """ CoefficientBasedFluxes(FT = Oceananigans.defaults.FloatType; - transfer_coefficients = SimilarityScales(1e-3, 1e-3, 1e-3), + transfer_coefficients = (1e-3, 1e-3, 1e-3), solver_stop_criteria = nothing, solver_tolerance = 1e-8, solver_maxiter = 20) @@ -182,8 +187,8 @@ Keyword Arguments ================= - `transfer_coefficients`: Transfer coefficients for momentum, heat, and moisture. - Can be a `SimilarityScales` with constant or callable entries, or an `LargeYeagerTransferCoefficients`. - Defaults to `SimilarityScales(1e-3, 1e-3, 1e-3)`. + Can be a `SimilarityScales`, a `Tuple`, or a `NamedTuple` with constant or callable entries, + or an `LargeYeagerTransferCoefficients`. Defaults to `(1e-3, 1e-3, 1e-3)`. - `solver_stop_criteria`: Criteria for iterative solver convergence. If `nothing`, creates new criteria using `solver_tolerance` and `solver_maxiter`. - `solver_tolerance`: Tolerance for solver convergence when creating new stop criteria, defaults to 1e-8. @@ -199,7 +204,7 @@ using NumericalEarth grid = RectilinearGrid(size=3, z=(-1, 0), topology=(Flat, Flat, Bounded)) ocean = ocean_simulation(grid; timestepper = :QuasiAdamsBashforth2) -ao_fluxes = CoefficientBasedFluxes(transfer_coefficients = SimilarityScales(1e-2, 1e-3, 1e-3)) +ao_fluxes = CoefficientBasedFluxes(transfer_coefficients = (1e-2, 1e-3, 1e-3)) interfaces = ComponentInterfaces(nothing, ocean; atmosphere_ocean_fluxes=ao_fluxes) @@ -208,21 +213,50 @@ ComponentInterfaces ``` """ function CoefficientBasedFluxes(FT = Oceananigans.defaults.FloatType; - transfer_coefficients = SimilarityScales(1e-3, 1e-3, 1e-3), + transfer_coefficients = (1e-3, 1e-3, 1e-3), solver_stop_criteria = nothing, solver_tolerance = 1e-8, solver_maxiter = 20) + transfer_coefficients = validate_coefficients(FT, transfer_coefficients) + if isnothing(solver_stop_criteria) solver_tolerance = convert(FT, solver_tolerance) solver_stop_criteria = ConvergenceStopCriteria(solver_tolerance, solver_maxiter) end - transfer_coefficients = convert_transfer_coefficients(FT, transfer_coefficients) - return CoefficientBasedFluxes(transfer_coefficients, solver_stop_criteria) end +validate_coefficients(FT, ly::LargeYeagerTransferCoefficients) = ly +validate_coefficients(FT, sc::SimilarityScales) = convert_transfer_coefficients(FT, sc) + +function validate_coefficients(FT, nt::NamedTuple) + required = (:momentum, :temperature, :water_vapor) + missing = filter(k -> !haskey(nt, k), required) + + if !isempty(missing) + throw(ArgumentError( + "Transfer coefficients NamedTuple must contain keys $(required). " * + "Missing keys: $(missing). Received: $(keys(nt))" + )) + end + + return convert_transfer_coefficients(FT, nt) +end + +function validate_coefficients(FT, tc::Tuple) + if length(tc) != 3 + throw(ArgumentError( + "Transfer coefficients must be a tuple of length 3: " * + "(momentum, temperature, water_vapor). " * + "Got length $(length(tc)) with value $(tc)." + )) + end + + return convert_transfer_coefficients(FT, tc) +end + ##### ##### Evaluate transfer coefficients (dispatch on coefficient type) ##### @@ -235,13 +269,20 @@ end @inline evaluate_coefficient(C::Function, args...) = C(args...) @inline evaluate_coefficient(C::PolynomialNeutralDragCoefficient, ΔU, args...) = C(ΔU) -@inline function evaluate_coefficients(coeffs::SimilarityScales, args...) +@inline function evaluate_coefficients(coeffs::Union{SimilarityScales, NamedTuple}, args...) Cd = evaluate_coefficient(coeffs.momentum, args...) Ch = evaluate_coefficient(coeffs.temperature, args...) Cq = evaluate_coefficient(coeffs.water_vapor, args...) return Cd, Ch, Cq end +@inline function evaluate_coefficients(coeffs::Tuple, args...) + Cd = evaluate_coefficient(coeffs[1], args...) + Ch = evaluate_coefficient(coeffs[2], args...) + Cq = evaluate_coefficient(coeffs[3], args...) + return Cd, Ch, Cq +end + # NCAR transfer coefficients: full L&Y stability correction (eqs. 6c-6d, 10a-10c) @inline function evaluate_coefficients(ly::LargeYeagerTransferCoefficients, ΔU, Tₛ, qₛ, Δh, From 08b5dca9420859be3a4135adc4ad3825b54da139 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Mon, 20 Apr 2026 21:09:21 +0200 Subject: [PATCH 08/11] fix typo --- .../coefficient_based_turbulent_fluxes.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/EarthSystemModels/InterfaceComputations/coefficient_based_turbulent_fluxes.jl b/src/EarthSystemModels/InterfaceComputations/coefficient_based_turbulent_fluxes.jl index 37bed31f9..3191db6c0 100644 --- a/src/EarthSystemModels/InterfaceComputations/coefficient_based_turbulent_fluxes.jl +++ b/src/EarthSystemModels/InterfaceComputations/coefficient_based_turbulent_fluxes.jl @@ -159,7 +159,7 @@ convert_if_number(FT, a) = a convert_transfer_coefficients(FT, c) = c convert_transfer_coefficients(FT, c::Tuple) = convert_if_number.(FT, c) -convert_transfer_coefficients(FT, c::NamedTuple) = (; momentum = convert_if_number(FT, c.momentum) +convert_transfer_coefficients(FT, c::NamedTuple) = (; momentum = convert_if_number(FT, c.momentum), temperature = convert_if_number(FT, c.temperature), water_vapor = convert_if_number(FT, c.water_vapor)) @@ -277,7 +277,7 @@ end end @inline function evaluate_coefficients(coeffs::Tuple, args...) - Cd = evaluate_coefficient(coeffs[1], args...) + Cd = evaluate_coefficient(coeffs[1], args...) Ch = evaluate_coefficient(coeffs[2], args...) Cq = evaluate_coefficient(coeffs[3], args...) return Cd, Ch, Cq From c75c70fbebee5dda6f2395ea307feecb486445ee Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Mon, 20 Apr 2026 21:50:45 +0200 Subject: [PATCH 09/11] remove SimilarityScale --- docs/src/interface_fluxes.md | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/docs/src/interface_fluxes.md b/docs/src/interface_fluxes.md index 3613b1999..0b1081331 100644 --- a/docs/src/interface_fluxes.md +++ b/docs/src/interface_fluxes.md @@ -137,7 +137,7 @@ A comprehensive example is given below, but we note briefly here that ```@example interface_fluxes using NumericalEarth -coefficient_fluxes = CoefficientBasedFluxes(transfer_coefficients = SimilarityScales(2e-3, 2e-3, 1e-3)) +coefficient_fluxes = CoefficientBasedFluxes(transfer_coefficients = (2e-3, 2e-3, 1e-3)) ``` Alternatively, the drag coefficient can be specified as a wind-speed-dependent polynomial @@ -146,7 +146,7 @@ at each iteration rather than using a constant: ```@example interface_fluxes poly_drag = PolynomialNeutralDragCoefficient() -poly_fluxes = CoefficientBasedFluxes(transfer_coefficients = SimilarityScales(poly_drag, 1e-3, 1e-3)) +poly_fluxes = CoefficientBasedFluxes(transfer_coefficients = (poly_drag, 1e-3, 1e-3)) ``` For the full Large & Yeager (2004) bulk algorithm with stability corrections, @@ -355,7 +355,7 @@ neutral_similarity_fluxes = SimilarityTheoryFluxes(stability_functions=nothing; interfaces = ComponentInterfaces(atmosphere, ocean; atmosphere_ocean_fluxes=neutral_similarity_fluxes) increased_roughness_model = OceanOnlyModel(ocean; atmosphere, interfaces) -coefficient_fluxes = CoefficientBasedFluxes(transfer_coefficients = SimilarityScales(2e-3, 2e-3, 2e-3)) +coefficient_fluxes = CoefficientBasedFluxes(transfer_coefficients = (2e-3, 2e-3, 2e-3)) interfaces = ComponentInterfaces(atmosphere, ocean; atmosphere_ocean_fluxes=coefficient_fluxes) coefficient_model = OceanOnlyModel(ocean; atmosphere, interfaces) From 7578fa8f72ffb6db7fd66b5b9af8021acf75bec6 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Mon, 20 Apr 2026 23:11:08 +0200 Subject: [PATCH 10/11] import PolynomialNeutralDragCoefficient --- docs/src/interface_fluxes.md | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/docs/src/interface_fluxes.md b/docs/src/interface_fluxes.md index 0b1081331..7854442d7 100644 --- a/docs/src/interface_fluxes.md +++ b/docs/src/interface_fluxes.md @@ -145,6 +145,8 @@ following Large & Yeager (2004). In this case `CoefficientBasedFluxes` evaluates at each iteration rather than using a constant: ```@example interface_fluxes +using NumericalEarth.EarthSystemModels.InterfaceComputations: PolynomialNeutralDragCoefficient + poly_drag = PolynomialNeutralDragCoefficient() poly_fluxes = CoefficientBasedFluxes(transfer_coefficients = (poly_drag, 1e-3, 1e-3)) ``` @@ -155,7 +157,7 @@ use `LargeYeagerTransferCoefficients`. This computes all three transfer coeffici stability corrections (L&Y eqs. 6c-6d, 10a-10c): ```@example interface_fluxes -using NumericalEarth.EarthSystemModels.InterfaceComputations: FixedIterations +using NumericalEarth.EarthSystemModels.InterfaceComputations: FixedIterations, LargeYeagerTransferCoefficients ly = LargeYeagerTransferCoefficients() ly_fluxes = CoefficientBasedFluxes(transfer_coefficients = ly, From c571c574a340df59c8ec32ac7210c52342ceb49d Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Tue, 21 Apr 2026 07:44:22 +0200 Subject: [PATCH 11/11] fix method --- .../InterfaceComputations/coefficient_based_turbulent_fluxes.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/src/EarthSystemModels/InterfaceComputations/coefficient_based_turbulent_fluxes.jl b/src/EarthSystemModels/InterfaceComputations/coefficient_based_turbulent_fluxes.jl index 3191db6c0..245c06149 100644 --- a/src/EarthSystemModels/InterfaceComputations/coefficient_based_turbulent_fluxes.jl +++ b/src/EarthSystemModels/InterfaceComputations/coefficient_based_turbulent_fluxes.jl @@ -262,6 +262,7 @@ end ##### # Minimum wind speed floor: only LargeYeager needs it; constant coefficients don't. +@inline minimum_wind_speed(::Union{Tuple, NTuple}) = 0 @inline minimum_wind_speed(::SimilarityScales) = 0 @inline minimum_wind_speed(ly::LargeYeagerTransferCoefficients) = ly.neutral_drag_coefficient.minimum_wind_speed