diff --git a/.dev/Project.toml b/.dev/Project.toml index d03b3145..71708c80 100644 --- a/.dev/Project.toml +++ b/.dev/Project.toml @@ -2,4 +2,4 @@ JuliaFormatter = "98e50ef6-434e-11e9-1051-2b60c6c9e899" [compat] -JuliaFormatter = "0.22" +JuliaFormatter = "1" diff --git a/Project.toml b/Project.toml index d45284ce..18aa0e3b 100644 --- a/Project.toml +++ b/Project.toml @@ -4,10 +4,12 @@ authors = ["Climate Modeling Alliance"] version = "0.12.0" [deps] +CLIMAParameters = "6eacf6c3-8458-43b9-ae03-caf5306d3d53" +ClimaParams = "5c42b081-d73a-476f-9059-fd94b934656c" DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae" +QuadGK = "1fd47b50-473d-5c70-9696-f719f8f3bcdc" RootSolvers = "7181ea78-2dcb-4de3-ab41-2b8ab5a31e74" Thermodynamics = "b60c26fb-14c3-4610-9d3e-2d17fe7ff00c" -ClimaParams = "5c42b081-d73a-476f-9059-fd94b934656c" [weakdeps] ClimaParams = "5c42b081-d73a-476f-9059-fd94b934656c" @@ -18,6 +20,7 @@ CreateParametersExt = "ClimaParams" [compat] ClimaParams = "0.10" DocStringExtensions = "0.8, 0.9" +QuadGK = "2.11.2" RootSolvers = "0.4" Thermodynamics = "0.12.8" julia = "1.9" diff --git a/docs/bibliography.bib b/docs/bibliography.bib index 0af76ef9..f744e9ef 100644 --- a/docs/bibliography.bib +++ b/docs/bibliography.bib @@ -2,6 +2,29 @@ # Last author name (titlecase), followed by # (no characters in-between) the year. + +@article{Bonan2019, + title = {Climate Change and Terrestrial Ecosystem Modeling}, + author = {Bonan, G}, + journal = {Cambridge University Press}, + doi = {10.1017/9781107339217}, + year = {2019}, + publisher = {Cambridge University Press} +} + + +@article{Physick1995, + title = {Incorporation of a high-roughness lower boundary into a mesoscale model for studies of dry deposition over complex terrain}, + author = {Physick, W and Garratt, R}, + journal = {Boundary-Layer Meteorology}, + volume = {74}, + number = {1}, + doi = {10.1007/BF00715710}, + pages = {55-71}, + year = {1995}, + publisher = {Springer} +} + @article{Cheng2005, title = {Flux-profile Relationships for Wind Speed and Temperature in the Stable Atmospheric Boundary Layer}, author = {Chenge, Y and Brutsaert, W}, diff --git a/docs/make.jl b/docs/make.jl index e21d98f6..48b44162 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -11,7 +11,8 @@ pages = Any[ "Home" => "index.md", "References" => "References.md", "Equations" => "SurfaceFluxes.md", - "Universal Functions" => "UniversalFunctions.md" + "Universal Functions" => "UniversalFunctions.md", + "Sample Profiles" => "SampleProfiles.md" ] mathengine = MathJax(Dict( diff --git a/docs/src/SampleProfiles.md b/docs/src/SampleProfiles.md new file mode 100644 index 00000000..22e8fa9a --- /dev/null +++ b/docs/src/SampleProfiles.md @@ -0,0 +1,29 @@ +# Sample Vertical Profiles of Wind Speed and Temperature + +SurfaceFluxes.jl provides profile recovery functions for the roughness sublayer using the Physick and Garratt (1995) formulation. Here, we use these functions to reproduce vertical profiles of wind speed and temperature for different conditions of atmospheric stability, reproducing Bonan 2019 Figure 6.4 with and without canopy correction. + +```@example +include("plot_profiles.jl") +``` + +# Fig 6.4 (a) + +![](Fig6.4a_profile.svg) +![](Fig6.4a_canopy_profile.svg) + +# Fig 6.4 (b) + +![](Fig6.4b_profile.svg) +![](Fig6.4b_canopy_profile.svg) + +# Fig 6.4 (c) + +![](Fig6.4c_profile.svg) +![](Fig6.4c_canopy_profile.svg) + +# Fig 6.4 (d) + +![](Fig6.4d_profile.svg) +![](Fig6.4d_canopy_profile.svg) + + diff --git a/docs/src/plot_profiles.jl b/docs/src/plot_profiles.jl new file mode 100644 index 00000000..f6f7fd33 --- /dev/null +++ b/docs/src/plot_profiles.jl @@ -0,0 +1,292 @@ +using Plots + +using SurfaceFluxes +const SF = SurfaceFluxes +SurfaceFluxes.error_on_non_convergence() = true + +import SurfaceFluxes.UniversalFunctions as UF +import Thermodynamics +Thermodynamics.print_warning() = false + +include(joinpath(pkgdir(SurfaceFluxes), "parameters", "create_parameters.jl")) +const SFP = SurfaceFluxes.Parameters +const FT = Float64; + +### Generate parameter lists +toml_dict = CP.create_toml_dict(FT; dict_type = "alias") +param_set = create_parameters(toml_dict, UF.Gryanik()) +thermo_params = SFP.thermodynamics_params(param_set) +uft = SFP.universal_func_type(param_set) +### + +""" +Here, we reproduce Figure 6.4 from Bonan (2019) Chapter 6. + +# References +- [Bonan2019](@cite) (Chapter 6 Figure 6.4) + +# Original Research +- [Physick2019](@cite) + +""" +T₀ = FT(273.15) +z_star = FT(49) +d = FT(19) +h_c = FT(22) +z0m = FT(0.6) # Figure 6.4 +z0b = FT(0.135) * z0m + +u_star_stable = FT(0.13) +u_star_unstable = FT(0.4) +θ_star_stable = FT(0.06) +θ_star_unstable = FT(-0.5) + +κ = SFP.von_karman_const(param_set) + +# For these measured values, see Chapter 6.6 page 89 of Bonan (2019) +z_measured = FT[21, 29, 21]; +u_measured = FT[1.0, 2.1, 1.0]; +θ_measured = FT[29.0, 28.1, 29.5] .+ T₀ + +p_sfc = FT(99340) # Pa +q_sfc = FT(0.0107) # kg/kg + +θ_sfc = FT(289.7) # K +ts_sfc_test = Thermodynamics.PhaseEquil_pθq(thermo_params, 99340.0, 289.7, 0.0107) +ts_int_test = Thermodynamics.PhaseEquil_pθq(thermo_params, 95342.0, 298.0, 0.0085) +state_in = SurfaceFluxes.InteriorValues(FT(100), (FT(1.0), FT(0)), ts_int_test) +state_sfc = SurfaceFluxes.SurfaceValues(FT(0), (FT(0), FT(0)), ts_sfc_test) +sc = SurfaceFluxes.ValuesOnly{FT}(; state_in, state_sfc, z0m, z0b) +ts_sfc_test = Thermodynamics.PhaseEquil_pθq(thermo_params, 100000.0, 289.7, 0.0) +ts_int_test = Thermodynamics.PhaseEquil_pθq(thermo_params, 99990.0, 298.0, 0.0) +Z = collect(range(FT(d + 10^-10), stop = FT(100), length = 250)) + +""" +save_profile(param_set, sc, ca, L_MOs, Z, X_sfc, transport, uft, scheme, x_star, d) + +Saves profiles of variable X given values of Z coordinates. Follows Nishizawa equation (21,22) + +## Arguments + - param_set: Abstract Parameter Set containing physical, thermodynamic parameters. + - sc: Container for surface conditions based on known combination + of the state vector, and {fluxes, friction velocity, exchange coefficients} for a given experiment + - L_MOs: Monin-Obukhov length(s) + - Z: Z coordinate(s) (within surface layer) for which variable values are required + - X_sfc: For variable X, values at interior and surface nodes + - transport: Transport type, (e.g. Momentum or Heat, used to determine physical scale coefficients) + - uft: A Universal Function type, (returned by, e.g., Businger()) + - scheme: Discretization scheme (currently supports FD and FV) + - rsl : Roughness Sublayer Formulation (e.g. NoRSL, PhysickRSL, DeRidderRSL) + - x_star: characteristic scale for variable x + - d: Displacement height (measure of the spatial lengthscale of the effect of the canopy roughness on near-wall turbulence) +""" +function save_profile( + param_set::SurfaceFluxes.APS, + sc::SurfaceFluxes.AbstractSurfaceConditions, + L_MOs::Array{FT, 1}, + Z::Array{FT, 1}, + X_sfc, + transport, + uft::UF.AUFT, + scheme::Union{SurfaceFluxes.FVScheme, SurfaceFluxes.FDScheme}, + rsl::SurfaceFluxes.AbstractRoughnessSublayerType, + x_star, + d; + title = nothing, + xlims = nothing, + ylims = nothing, + xlabel = "u(z)", + ylabel = "z", + fig_prefix = "", + xaxis = :identity, + yaxis = :identity, +) + Plots.plot() + for L_MO in L_MOs + x_i = map(Z) do z + Zi = typeof(rsl) == SurfaceFluxes.NoRSL ? FT(z - d) : FT(z) + state_in = SurfaceFluxes.InteriorValues(FT(Zi), (FT(1.0), FT(0)), ts_int_test) + state_sfc = SurfaceFluxes.SurfaceValues(FT(0), (FT(0), FT(0)), ts_sfc_test) + sc = SurfaceFluxes.ValuesOnly{FT}(; state_in, state_sfc, z0m, z0b) + rsc = SurfaceFluxes.surface_conditions(param_set, sc, SurfaceFluxes.FDScheme()) + dx = SurfaceFluxes.recover_profile(param_set, sc, L_MO, Zi, X_sfc, x_star, transport, uft, scheme, rsl) + end + + uf = UF.universal_func(uft, L_MO, SFP.uf_params(param_set)) + _π_group = FT(UF.π_group(uf, transport)) + + Δx = @. (x_i - X_sfc) + + Plots.plot!(Δx, Z, label = "L_MO = $L_MO") + Plots.plot!(; title, xlabel, ylabel, ylims, xlims, grid = :off, legend = :outerright, titlefontalign = :center) + + Plots.savefig("$(fig_prefix)_profile.png") + end +end + + + +# Save profile for unstable and stable L_MO set +Stable_L_MOs = FT[30, 50, 1000] +Unstable_L_MOs = FT[-10, -50, -1000] + +testcanopy = SurfaceFluxes.SparseCanopy{FT}(d, z_star) +PhysickRSL = SurfaceFluxes.PhysickRSL(testcanopy) +NoRSL = SurfaceFluxes.NoRSL() + + +# Bonan2019 Fig. 6.4a (With PG95 RSL) +save_profile( + param_set, + sc, + Unstable_L_MOs, + Z, + FT(0), + UF.MomentumTransport(), + uft, + SurfaceFluxes.FDScheme(), + PhysickRSL, + u_star_unstable, + d; + xlims = (0, 4), + ylims = (15, 50), + fig_prefix = "Fig6.4a_canopy", + title = "Vertical Profile of Wind Velocity (PG95 RSL)", +) + +# Bonan2019 Fig. 6.4a (No RSL Model) +save_profile( + param_set, + sc, + Unstable_L_MOs, + Z, + FT(0), + UF.MomentumTransport(), + uft, + SurfaceFluxes.FDScheme(), + NoRSL, + u_star_unstable, + d; + xlims = (0, 4), + ylims = (15, 50), + fig_prefix = "Fig6.4a", + title = "Vertical Profile of Wind Velocity (No RSL Model)", +) + +# Bonan2019 Fig. 6.4b (With PG95 RSL) +save_profile( + param_set, + sc, + Unstable_L_MOs, + Z, + θ_sfc, + UF.HeatTransport(), + uft, + SurfaceFluxes.FDScheme(), + PhysickRSL, + θ_star_unstable, + d; + xlims = (-8, 0), + ylims = (15, 50), + xlabel = "θ - θ_sfc", + fig_prefix = "Fig6.4b_canopy", + title = "Vertical Profile of Temperature (PG95 RSL)", +) + +# Bonan2019 Fig. 6.4b (No RSL Model) +save_profile( + param_set, + sc, + Unstable_L_MOs, + Z, + θ_sfc, + UF.HeatTransport(), + uft, + SurfaceFluxes.FDScheme(), + NoRSL, + θ_star_unstable, + d; + xlims = (-8, 0), + ylims = (15, 50), + xlabel = "θ- θ_sfc", + fig_prefix = "Fig6.4b", + title = "Vertical Profile of Temperature (No RSL Model)", +) + +# Bonan2019 Fig. 6.4c (With PG95 RSL) +save_profile( + param_set, + sc, + Stable_L_MOs, + Z, + FT(0), + UF.MomentumTransport(), + uft, + SurfaceFluxes.FDScheme(), + PhysickRSL, + u_star_stable, + d; + xlims = (0, 4), + ylims = (15, 50), + fig_prefix = "Fig6.4c_canopy", + title = "Vertical Profile of Wind Velocity (PG95 RSL)", +) + +# Bonan2019 Fig. 6.4c (No RSL Model) +save_profile( + param_set, + sc, + Stable_L_MOs, + Z, + FT(0), + UF.MomentumTransport(), + uft, + SurfaceFluxes.FDScheme(), + NoRSL, + u_star_stable, + d; + xlims = (0, 4), + ylims = (15, 50), + fig_prefix = "Fig6.4c", + title = "Vertical Profile of Wind Velocity (No RSL Model)", +) + +# Bonan2019 Fig. 6.4d (With PG95 RSL) +save_profile( + param_set, + sc, + Stable_L_MOs, + Z, + θ_sfc, + UF.HeatTransport(), + uft, + SurfaceFluxes.FDScheme(), + PhysickRSL, + θ_star_stable, + d; + xlims = (0, 2), + ylims = (15, 50), + xlabel = "θ-θ_sfc", + fig_prefix = "Fig6.4d_canopy", + title = "Vertical Profile of Temperature (PG95 RSL)", +) + +# Bonan2019 Fig. 6.4d (No RSL Model) +save_profile( + param_set, + sc, + Stable_L_MOs, + Z, + θ_sfc, + UF.HeatTransport(), + uft, + SurfaceFluxes.FDScheme(), + NoRSL, + θ_star_stable, + d; + xlims = (0, 2), + ylims = (15, 50), + xlabel = "θ-θ_sfc", + fig_prefix = "Fig6.4d", + title = "Vertical Profile of Temperature (No RSL Model)", +) diff --git a/docs/src/plot_universal_functions.jl b/docs/src/plot_universal_functions.jl index 7adcf436..3345b708 100644 --- a/docs/src/plot_universal_functions.jl +++ b/docs/src/plot_universal_functions.jl @@ -27,7 +27,7 @@ function save_ϕ_figs(ufps, ζ; ylims = nothing, fig_prefix = "", xaxis = :ident uf = UF.universal_func(uft, L, ufp) ϕ_m = UF.phi.(uf, ζ, UF.MomentumTransport()) label = "$(typeof(uf).name)" - Plots.plot!(ζ, ϕ_m; xlabel = "ζ", ylabel = "ϕ_m", label, ylims, xaxis, yaxis) + Plots.plot!(ζ, ϕ_m; xlabel = "ζ = z/L_MO", ylabel = "ϕ_m", label, ylims, xaxis, yaxis) end Plots.savefig("$(fig_prefix)_phi_m.svg") Plots.plot() @@ -36,7 +36,7 @@ function save_ϕ_figs(ufps, ζ; ylims = nothing, fig_prefix = "", xaxis = :ident uf = UF.universal_func(uft, L, ufp) ϕ_h = UF.phi.(uf, ζ, UF.HeatTransport()) label = "$(typeof(uf).name)" - Plots.plot!(ζ, ϕ_h; xlabel = "ζ", ylabel = "ϕ_h", label, ylims, xaxis, yaxis) + Plots.plot!(ζ, ϕ_h; xlabel = "ζ = z/L_MO", ylabel = "ϕ_h", label, ylims, xaxis, yaxis) end Plots.savefig("$(fig_prefix)_phi_h.svg") end @@ -47,7 +47,7 @@ function save_ψ_figs(ufps, ζ; ylims = nothing, fig_prefix = "", xaxis = :ident uf = UF.universal_func(uft, L, ufp) ψ_m = UF.psi.(uf, ζ, UF.MomentumTransport()) label = "$(typeof(uf).name)" - Plots.plot!(ζ, ψ_m; xlabel = "ζ", ylabel = "ψ_m", label, ylims, xaxis, yaxis) + Plots.plot!(ζ, ψ_m; xlabel = "ζ = z/L_MO", ylabel = "ψ_m", label, ylims, xaxis, yaxis) end Plots.savefig("$(fig_prefix)_psi_m.svg") Plots.plot() @@ -56,7 +56,7 @@ function save_ψ_figs(ufps, ζ; ylims = nothing, fig_prefix = "", xaxis = :ident uf = UF.universal_func(uft, L, ufp) ψ_h = UF.psi.(uf, ζ, UF.HeatTransport()) label = "$(typeof(uf).name)" - Plots.plot!(ζ, ψ_h; xlabel = "ζ", ylabel = "ψ_h", label, ylims, xaxis, yaxis) + Plots.plot!(ζ, ψ_h; xlabel = "ζ = z/L_MO", ylabel = "ψ_h", label, ylims, xaxis, yaxis) end Plots.savefig("$(fig_prefix)_psi_h.svg") end diff --git a/src/SurfaceFluxes.jl b/src/SurfaceFluxes.jl index 915aa808..4b001f0c 100644 --- a/src/SurfaceFluxes.jl +++ b/src/SurfaceFluxes.jl @@ -37,6 +37,43 @@ const SFP = Parameters const APS = SFP.AbstractSurfaceFluxesParameters abstract type SolverScheme end +struct FVScheme <: SolverScheme end +struct FDScheme <: SolverScheme end + +using QuadGK + +abstract type CanopyType end +struct SparseCanopy{FT} <: CanopyType + d::FT + z_star::FT +end +struct DenseCanopy{FT} <: CanopyType + d::FT + z_star::FT +end + +abstract type AbstractRoughnessSublayerType end +struct NoRSL <: AbstractRoughnessSublayerType end +struct PhysickRSL <: AbstractRoughnessSublayerType + canopy::CanopyType +end +struct DeRidderRSL <: AbstractRoughnessSublayerType + canopy::CanopyType +end +struct HarmonRSL <: AbstractRoughnessSublayerType + canopy::CanopyType +end + +# Allow users to skip error on non-convergence +# by importing: +# ```julia +# import SurfaceFluxes +# SurfaceFluxes.error_on_non_convergence() = false +# ``` +# Error on convergence must be the default +# behavior because this can result in printing +# very large logs resulting in CI to seemingly hang. +error_on_non_convergence() = true struct LayerAverageScheme <: SolverScheme end struct PointValueScheme <: SolverScheme end @@ -557,8 +594,9 @@ compute_ustar(param_set, L_MO, sc::Fluxes, uft, scheme) = compute_ustar(param_set, L_MO, sc::Coefficients, uft, scheme) = sqrt(sc.Cd) * (windspeed(sc)) -compute_ustar(param_set, L_MO, sc::ValuesOnly, uft, scheme) = - windspeed(sc) * compute_physical_scale_coeff(param_set, sc, L_MO, UF.MomentumTransport(), uft, scheme) +function compute_ustar(param_set, L_MO, sc::ValuesOnly, uft, scheme) + return windspeed(sc) * compute_physical_scale_coeff(param_set, sc, L_MO, UF.MomentumTransport(), uft, scheme) +end """ momentum_exchange_coefficient(param_set, L_MO, sc, uft, scheme) @@ -843,16 +881,17 @@ Recover profiles of variable X given values of Z coordinates. Follows Nishizawa - X_sfc: For variable X, values at surface nodes - transport: Transport type, (e.g. Momentum or Heat, used to determine physical scale coefficients) - scheme: Discretization scheme (currently supports FD and FV) - -# TODO: add tests + - rsl : Roughness Sublayer Formulation (e.g. NoRSL, PhysickRSL, DeRidderRSL) + # TODO: add tests + # TODO: Verify that all current RSL models fall into this general code pattern, then 𝜙 can be abstracted """ function recover_profile( param_set::APS, sc::AbstractSurfaceConditions, L_MO, Z, - X_star, X_sfc, + X_star, transport, scheme::Union{LayerAverageScheme, PointValueScheme}, ) @@ -867,6 +906,145 @@ function recover_profile( return Σnum * X_star / von_karman_const + X_sfc end + +""" + recover_profile(param_set, sc, L_MO, Z, X_sfc, transport, uft, scheme, rsl) + +Recover profiles of variable X given values of Z coordinates, as long as Z > d. +Follows Nishizawa equation (21, 22) +Canopy modification follows Physick and Garratt (1995) equation (8, 9) +## Arguments + - param_set: Abstract Parameter Set containing physical, thermodynamic parameters. + - sc: Container for surface conditions based on known combination + of the state vector, and {fluxes, friction velocity, exchange coefficients} for a given experiment + - L_MO: Monin-Obukhov length + - Z: Z coordinate(s) (within surface layer) for which variable values are required + - X_sfc: For variable X, values at surface nodes + - transport: Transport type, (e.g. Momentum or Heat, used to determine physical scale coefficients) + - uft: A Universal Function type, (returned by, e.g., Businger()) + - scheme: Discretization scheme (currently supports FD and FV) + - rsl : Roughness Sublayer Formulation (e.g. NoRSL, PhysickRSL, DeRidderRSL) + +# TODO: add tests +# TODO: Verify that all current RSL models fall into this general code pattern, then 𝜙 can be abstracted +""" +function recover_profile( + param_set::APS, + sc::AbstractSurfaceConditions, + L_MO::FT, + Z, + X_sfc, + X_star, + transport, + uft::UF.AUFT, + scheme::Union{FVScheme, FDScheme}, + rsl::PhysickRSL, +) where {FT} + z_star = rsl.canopy.z_star + d = rsl.canopy.d + uf = UF.universal_func(uft, L_MO, SFP.uf_params(param_set)) + von_karman_const::FT = SFP.von_karman_const(param_set) + _π_group = FT(UF.π_group(uf, transport)) + _π_group⁻¹ = (1 / _π_group) + num1 = log((Z - d) / z0(sc, transport)) + num2 = -UF.psi(uf, (Z - d) / L_MO, transport) + num3 = UF.psi(uf, z0(sc, transport) / L_MO, transport) + Σnum = num1 + num2 + num3 + ### Protoype : QuadGK integration to evaluate the canopy RSL effect given by Physick and Garratt (1995) + ### PG95 assume the same 𝜙 function for momentum and scalar (velocity, heat, moisture etc.) variables + ### For a model level `Z`, we have the offset coordinate `z = Z-d` over which the functions in PG95 are defined + + # Physick 1995 RSL Term + function ψ_RSL_term(z) + if (z < d || z > z_star) + return 0 + end + + function integrand(x) + ϕ = UF.phi(uf, (x - d) / L_MO, transport) + + # Physick1995 Eq. 18 + ϕ_italic = FT(0.5) * exp(log(2) * (x - d) / (z_star - d)) + return ϕ * (1 - ϕ_italic) / (x - d) + end + + integral, error = quadgk(integrand, z, z_star) + return integral + end + rsl_physick = ψ_RSL_term(Z) + Σnum += rsl_physick + + return Σnum * X_star / von_karman_const + X_sfc +end + +""" + recover_profile(param_set, sc, L_MO, Z, X_sfc, transport, uft, scheme, rsl) + +Recover profiles of variable X given values of Z coordinates, as long as Z > d. +Follows Nishizawa equation (21, 22) +Canopy modification follows De Ridder (2010) equation (6, 7, 13) + +## Arguments + - param_set: Abstract Parameter Set containing physical, thermodynamic parameters. + - sc: Container for surface conditions based on known combination + of the state vector, and {fluxes, friction velocity, exchange coefficients} for a given experiment + - L_MO: Monin-Obukhov length + - Z: Z coordinate(s) (within surface layer) for which variable values are required + - X_sfc: For variable X, values at surface nodes + - transport: Transport type, (e.g. Momentum or Heat, used to determine physical scale coefficients) + - uft: A Universal Function type, (returned by, e.g., Businger()) + - scheme: Discretization scheme (currently supports FD and FV) + - rsl : Roughness Sublayer Formulation (e.g. NoRSL, PhysickRSL, DeRidderRSL) + +# TODO: add tests +# TODO: Verify that all current RSL models fall into this general code pattern, then 𝜙 can be abstracted +""" +function recover_profile( + param_set::APS, + sc::AbstractSurfaceConditions, + L_MO::FT, + Z, + X_sfc, + X_star, + transport, + uft::UF.AUFT, + scheme::Union{FVScheme, FDScheme}, + rsl::DeRidderRSL, +) where {FT} + z_star = rsl.canopy.z_star + d = rsl.canopy.d + uf = UF.universal_func(uft, L_MO, SFP.uf_params(param_set)) + von_karman_const::FT = SFP.von_karman_const(param_set) + _π_group = FT(UF.π_group(uf, transport)) + _π_group⁻¹ = (1 / _π_group) + num1 = log((Z - d) / z0(sc, transport)) + num2 = -UF.psi(uf, (Z - d) / L_MO, transport) + num3 = UF.psi(uf, z0(sc, transport) / L_MO, transport) + Σnum = num1 + num2 + num3 + + ## De Ridder (2010) RSL Term + function ψ_analytic(z) + # TODO Move to ClimaParameters + # De Ridder Eq. 13 + + ν = FT(0.5) + if transport isa UF.MomentumTransport + μ = FT(2.59) # Momentum, μₕ = 0.95 + elseif transport isa UF.HeatTransport + μ = FT(0.95) # Momentum, μₕ = 0.95 + end + λ = FT(1.5) + return UF.phi(uf, z / L_MO * (1 + ν / (μ * z / z_star)), transport) * + (1 / λ) * + log(1 + (λ / (μ * z / z_star))) * + exp(-μ * z / z_star) + end + rsl_deridder = ψ_analytic(Z - d) + Σnum += rsl_deridder + + return Σnum * X_star / von_karman_const + X_sfc +end + # For backwards compatibility with package extensions if !isdefined(Base, :get_extension) include(joinpath("..", "ext", "CreateParametersExt.jl")) diff --git a/test/runtests.jl b/test/runtests.jl index 59d6106c..26f955d7 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -141,6 +141,7 @@ end @test all(x -> x <= FloatType(0.005), abs.(rdiff_sol)) end + @testset "Test profiles" begin include("test_profiles.jl") end @@ -153,10 +154,15 @@ end @testset "Test generated thermodynamic states" begin include("test_convergence.jl") end +<<<<<<< HEAD +@testset "Test roughness sublayer profiles" begin + include("test_rsl.jl") +======= @testset "Test generated thermodynamic states" begin include("test_convergence.jl") end @testset "Quality assurance" begin include("aqua.jl") +>>>>>>> origin/main end diff --git a/test/test_profiles.jl b/test/test_profiles.jl index 38a5d3f7..77776732 100644 --- a/test/test_profiles.jl +++ b/test/test_profiles.jl @@ -7,10 +7,19 @@ import ArtifactWrappers const AW = ArtifactWrappers const TD = Thermodynamics +<<<<<<< HEAD +include(joinpath(pkgdir(SurfaceFluxes), "parameters", "create_parameters.jl")) +FT = Float32; +toml_dict = CP.create_toml_dict(FT; dict_type = "alias") +param_set = create_parameters(toml_dict, UF.BusingerType()) +thermo_params = SFP.thermodynamics_params(param_set) + +======= FT = Float32 param_set = SFP.SurfaceFluxesParameters(FT, BusingerParams) thermo_params = param_set.thermo_params uft = UF.BusingerType() +>>>>>>> origin/main #! format: off PyCLES_output_dataset = AW.ArtifactWrapper( diff --git a/test/test_rsl.jl b/test/test_rsl.jl new file mode 100644 index 00000000..a55c2afe --- /dev/null +++ b/test/test_rsl.jl @@ -0,0 +1,125 @@ +using Test + +import SurfaceFluxes +const SF = SurfaceFluxes +SurfaceFluxes.error_on_non_convergence() = true + +import SurfaceFluxes.UniversalFunctions as UF +using StaticArrays +using Thermodynamics +Thermodynamics.print_warning() = false + +include(joinpath(pkgdir(SurfaceFluxes), "parameters", "create_parameters.jl")) +const SFP = SF.Parameters +FT = Float64; + +### Generate parameter lists +toml_dict = CP.create_toml_dict(FT; dict_type = "alias") +param_set = create_parameters(toml_dict, UF.Gryanik()) +thermo_params = SFP.thermodynamics_params(param_set) +uft = SFP.universal_func_type(param_set) +### + +""" +Parameters from Bonan (2019) Chapter 6 + +# References +- [Bonan2019](@cite) (Chapter 6) + +# Original Research +- [Physick1995](@cite) + +""" +T₀ = FT(273.15) +z_star = FT(49) +d = FT(19) +h_c = FT(22) +z0m = FT(0.6) # Figure 6.4 +z0b = FT(0.135) * z0m +κ = SFP.von_karman_const(param_set) + +u_sfc = FT(3.7) +u_in = FT(0) +θ_sfc = FT(289.7) # K +θ_in = FT(298.0) # K + +ts_sfc_test = Thermodynamics.PhaseEquil_pθq(thermo_params, 99340.0, 289.7, 0.0107) +ts_int_test = Thermodynamics.PhaseEquil_pθq(thermo_params, 95342.0, 298.0, 0.0085) +state_in = SF.InteriorValues(FT(350), (FT(0.0), FT(0)), ts_int_test) +state_sfc = SF.SurfaceValues(FT(0), (FT(0), FT(0)), ts_sfc_test) +sc = SF.ValuesOnly{FT}(; state_in, state_sfc, z0m, z0b) +Z = collect(range(FT(d + 0.1), stop = FT(100), length = 250)) + + +# TODO: expand test to include multiple canopies +@testset "De Ridder and Physick RSL Profiles Agree" begin + testcanopy = SurfaceFluxes.SparseCanopy{FT}(d, z_star) + PhysickRSL = SurfaceFluxes.PhysickRSL(testcanopy) + DeRidderRSL = SurfaceFluxes.DeRidderRSL(testcanopy) + + @testset "Wind Profiles" begin + for (iz, z) in enumerate(Z) + for L_MO in FT[-1000, -50, -10, 30, 50, 1000] + du_ridder = SF.recover_profile( + param_set, + sc, + L_MO, + FT(z), + u_in, + u_sfc, + UF.MomentumTransport(), + uft, + SF.FDScheme(), + DeRidderRSL, + ) + du_physick = SF.recover_profile( + param_set, + sc, + L_MO, + FT(z), + u_in, + u_sfc, + UF.MomentumTransport(), + uft, + SF.FDScheme(), + PhysickRSL, + ) + + @test isapprox(du_ridder, du_physick, rtol = 0.25) + end + end + end + + @testset "Temperature Profiles" begin + for (iz, z) in enumerate(Z) + for L_MO in FT[-1000, -50, -10, 30, 50, 1000] + dθ_ridder = SF.recover_profile( + param_set, + sc, + L_MO, + FT(z), + θ_in, + θ_sfc, + UF.HeatTransport(), + uft, + SF.FDScheme(), + DeRidderRSL, + ) + dθ_physick = SF.recover_profile( + param_set, + sc, + L_MO, + FT(z), + θ_in, + θ_sfc, + UF.HeatTransport(), + uft, + SF.FDScheme(), + PhysickRSL, + ) + + @test isapprox(dθ_ridder, dθ_physick, rtol = 0.25) + end + end + end +end diff --git a/test/test_universal_functions.jl b/test/test_universal_functions.jl index f8ec9830..92f15006 100644 --- a/test/test_universal_functions.jl +++ b/test/test_universal_functions.jl @@ -6,9 +6,21 @@ import SurfaceFluxes as SF import SurfaceFluxes.UniversalFunctions as UF import ClimaParams as CP +<<<<<<< HEAD +import CLIMAParameters +const CP = CLIMAParameters + +include(joinpath(pkgdir(SurfaceFluxes), "parameters", "create_parameters.jl")) +FT = Float32; +toml_dict = CP.create_toml_dict(FT; dict_type = "alias") +param_set = create_parameters(toml_dict, UF.BusingerType()) + +universal_functions(uft, L) = UF.universal_func(uft, L, create_uf_parameters(toml_dict, uft)) +======= FT = Float32 param_set = SFP.SurfaceFluxesParameters(FT, BusingerParams) thermo_params = param_set.thermo_params +>>>>>>> origin/main # TODO: Right now, we test these functions for # type stability and correctness in the asymptotic