diff --git a/Manifest.toml b/Manifest.toml index 2a098064cf8..de5ecf418fe 100644 --- a/Manifest.toml +++ b/Manifest.toml @@ -30,11 +30,17 @@ git-tree-sha1 = "e214a9b9bd1b4e1b4f15b22c0994862b66af7ff7" uuid = "68821587-b530-5797-8361-c406ea357684" version = "3.5.0+3" +[[ArrayInterface]] +deps = ["LinearAlgebra", "Requires", "SparseArrays"] +git-tree-sha1 = "649c08a5a3a513f4662673d3777fe6ccb4df9f5d" +uuid = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9" +version = "2.8.7" + [[ArrayLayouts]] deps = ["FillArrays", "LinearAlgebra"] -git-tree-sha1 = "cee952c726065d21a6b9a0abe04e201debdffc20" +git-tree-sha1 = "0517f50df07250903194a94c2bc793c2bc239a76" uuid = "4c555306-a7a7-4459-81d9-ec55ddd5c99a" -version = "0.3.0" +version = "0.3.1" [[Base64]] uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f" @@ -176,6 +182,12 @@ git-tree-sha1 = "eb0c34204c8410888844ada5359ac8b96292cfd1" uuid = "b552c78f-8df3-52c6-915a-8e097449b14b" version = "1.0.1" +[[Distances]] +deps = ["LinearAlgebra", "Statistics"] +git-tree-sha1 = "23717536c81b63e250f682b0e0933769eecd1411" +uuid = "b4f34e82-e78d-54a5-968a-f98e89d6e8f7" +version = "0.8.2" + [[Distributed]] deps = ["Random", "Serialization", "Sockets"] uuid = "8ba89e20-285c-5b6f-9357-94700520ee1b" @@ -217,9 +229,15 @@ version = "1.3.0" [[FillArrays]] deps = ["LinearAlgebra", "Random", "SparseArrays"] -git-tree-sha1 = "5322d34d7600d3429665b37bcf7628dc602a28cc" +git-tree-sha1 = "6c89d5b673e59b8173c546c84127e5f623d865f6" uuid = "1a297f60-69ca-5386-bcde-b61e274b549b" -version = "0.8.8" +version = "0.8.9" + +[[FiniteDiff]] +deps = ["ArrayInterface", "LinearAlgebra", "Requires", "SparseArrays", "StaticArrays"] +git-tree-sha1 = "e65805de69d457029940acff64dd92e57b93c8a5" +uuid = "6a86dc24-6348-571c-b903-95158fe2bd41" +version = "2.3.1" [[ForwardDiff]] deps = ["CommonSubexpressions", "DiffResults", "DiffRules", "NaNMath", "Random", "SpecialFunctions", "StaticArrays"] @@ -327,6 +345,12 @@ git-tree-sha1 = "be855e3c975b89746b09952407c156b5e4a33a1d" uuid = "9c8b4983-aa76-5018-a973-4c85ecc9e179" version = "0.8.1" +[[LineSearches]] +deps = ["LinearAlgebra", "NLSolversBase", "NaNMath", "Parameters", "Printf", "Test"] +git-tree-sha1 = "54eb90e8dbe745d617c78dee1d6ae95c7f6f5779" +uuid = "d3d80556-e9d4-5f37-9878-2ab0fcc64255" +version = "7.0.1" + [[LinearAlgebra]] deps = ["Libdl"] uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" @@ -395,6 +419,18 @@ git-tree-sha1 = "271e4670a6b30ec1e5daf6ba8b30240b7d85f535" uuid = "85f8d34a-cbdd-5861-8df4-14fed0d494ab" version = "0.10.2" +[[NLSolversBase]] +deps = ["DiffResults", "Distributed", "FiniteDiff", "ForwardDiff"] +git-tree-sha1 = "7c4e66c47848562003250f28b579c584e55becc0" +uuid = "d41bc354-129a-5804-8e4c-c37616107c6c" +version = "7.6.1" + +[[NLsolve]] +deps = ["Distances", "LineSearches", "LinearAlgebra", "NLSolversBase", "Printf", "Reexport"] +git-tree-sha1 = "cce0463af83a0f36c7bfa5820e373ac090cc46ad" +uuid = "2774e3e8-f4cf-5e23-947b-6d7e65073b56" +version = "4.3.0" + [[NNlib]] deps = ["BinaryProvider", "Libdl", "LinearAlgebra", "Requires", "Statistics"] git-tree-sha1 = "d9f196d911f55aeaff11b11f681b135980783824" @@ -429,6 +465,12 @@ git-tree-sha1 = "2fc6f50ddd959e462f0a2dbc802ddf2a539c6e35" uuid = "90014a1f-27ba-587c-ab20-58faa44d9150" version = "0.9.12" +[[Parameters]] +deps = ["OrderedCollections", "UnPack"] +git-tree-sha1 = "38b2e970043613c187bd56a995fe2e551821eb4a" +uuid = "d96e819e-fc66-5662-9728-84c9c7592b0a" +version = "0.12.1" + [[Parsers]] deps = ["Dates", "Test"] git-tree-sha1 = "72c3451932513427caffbd8bab15643ad693804b" @@ -474,6 +516,12 @@ git-tree-sha1 = "54f8ceb165a0f6d083f0d12cb4996f5367c6edbc" uuid = "3cdcf5f2-1ef4-517c-9805-6587b60abb01" version = "1.0.1" +[[Reexport]] +deps = ["Pkg"] +git-tree-sha1 = "7b1d07f411bc8ddb7977ec7f377b97b158514fe0" +uuid = "189a3867-3050-52da-a836-e630ba90ab69" +version = "0.2.0" + [[Requires]] deps = ["UUIDs"] git-tree-sha1 = "d37400976e98018ee840e0ca4f9d20baa231dc6b" @@ -586,6 +634,11 @@ version = "0.4.1" deps = ["Random", "SHA"] uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4" +[[UnPack]] +git-tree-sha1 = "bc9ef72a4a826740895bf2772b48c21f9a1c13a7" +uuid = "3a884ed6-31ef-47d7-9d2a-63182c4928ed" +version = "1.0.0" + [[Unicode]] uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5" diff --git a/Project.toml b/Project.toml index cf608282678..94d5b7a440c 100644 --- a/Project.toml +++ b/Project.toml @@ -29,6 +29,7 @@ Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306" Logging = "56ddb016-857b-54e1-b83d-db4d58db5568" MPI = "da04e1cc-30fd-572f-bb4f-1f8673147195" NCDatasets = "85f8d34a-cbdd-5861-8df4-14fed0d494ab" +NLsolve = "2774e3e8-f4cf-5e23-947b-6d7e65073b56" OrderedCollections = "bac558e1-5e72-5ebc-8fee-abe8a469f55d" Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f" Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7" diff --git a/docs/Project.toml b/docs/Project.toml index 6f70a2e992c..a96c937f473 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -18,6 +18,7 @@ LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306" Logging = "56ddb016-857b-54e1-b83d-db4d58db5568" MPI = "da04e1cc-30fd-572f-bb4f-1f8673147195" +NLsolve = "2774e3e8-f4cf-5e23-947b-6d7e65073b56" Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7" SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b" diff --git a/docs/src/HowToGuides/Common/SurfaceFluxes.md b/docs/src/HowToGuides/Common/SurfaceFluxes.md index 7632cd52916..98462fc421a 100644 --- a/docs/src/HowToGuides/Common/SurfaceFluxes.md +++ b/docs/src/HowToGuides/Common/SurfaceFluxes.md @@ -1,94 +1,36 @@ # Surface Fluxes -Surface flux functions, e.g. for buoyancy flux, friction velocity, and exchange coefficients. - -## `Byun1990` - -Compute surface fluxes using the approach in Byun (1990). - -### Plots - -```@example byun1990 -using ClimateMachine.SurfaceFluxes.Byun1990 -using Plots, LaTeXStrings -using CLIMAParameters -struct EarthParameterSet <: AbstractEarthParameterSet end -param_set = EarthParameterSet() - -FT = Float64 -Ri_range = range(FT(-1.2), stop=FT(0.24), length=100) -scales = FT[50,200,600,1000,10_000] - -z_0 = FT(1.0) -γ_m, γ_h = FT(15.0), FT(9.0) -β_m, β_h = FT(4.8), FT(7.8) -Pr_0 = FT(0.74) - -plot(Ri_range, - [Byun1990.compute_exchange_coefficients(param_set, Ri,scale*z_0,z_0,γ_m,γ_h,β_m,β_h,Pr_0)[1] - for Ri in Ri_range, scale in scales], - xlabel="Bulk Richardson number (Ri_b)", ylabel="Drag coefficient", title="Momentum exchange coefficient", - labels=scales, legendtitle=L"z/z_0") - -savefig("exchange_byun1990_fig4a.svg") # hide -nothing # hide -``` -![](exchange_byun1990_fig4a.svg) - -Recreation of Figure 4(a) from Byun (1990) - -```@example byun1990 -plot(Ri_range, - [Byun1990.compute_exchange_coefficients(param_set, Ri,scale*z_0,z_0,γ_m,γ_h,β_m,β_h,Pr_0)[2] - for Ri in Ri_range, scale in scales], - xlabel="Bulk Richardson number (Ri_b)", ylabel="Exchange coefficient", title="Heat exchange coefficient", - labels=scales, legendtitle=L"z/z_0") - -savefig("exchange_byun1990_fig4b.svg") # hide -nothing # hide +```@meta +CurrentModule = CLIMA.SurfaceFluxes ``` -![](exchange_byun1990_fig4b.svg) -Recreation of Figure 4(b) from Byun (1990) +This module provides a means to compute surface fluxes given several variables, described in [`surface_conditions`](@ref). -## `Nishizawa2018` +## Interface + - [`surface_conditions`](@ref) computes + - Monin-Obukhov length + - Potential temperature flux (if not given) using Monin-Obukhov theory + - transport fluxes using Monin-Obukhov theory + - friction velocity/temperature scale/tracer scales + - exchange coefficients -### Plots -```@example -using ClimateMachine.SurfaceFluxes.Nishizawa2018 -using Plots, LaTeXStrings -using CLIMAParameters -struct EarthParameterSet <: AbstractEarthParameterSet end -param_set = EarthParameterSet() -FT = Float64 +## API -a = FT(4.7) -θ = FT(350) -z_0 = FT(10) -u_ave = FT(10) -flux = FT(1) -Δz = range(FT(10.0), stop=FT(100.0), length=100) -Ψ_m_tol, tol_abs, iter_max = FT(1e-3), FT(1e-3), 10 -u_star = Nishizawa2018.compute_friction_velocity.( - Ref(param_set), - u_ave, θ, flux, Δz, z_0, a, Ψ_m_tol, tol_abs, iter_max) -plot(u_star, Δz, title = "Friction velocity vs dz", xlabel = "Friction velocity", ylabel = "dz") -savefig("friction_velocity.svg") # hide -nothing # hide +```@docs +surface_conditions ``` -![](friction_velocity.svg) ## References -- Businger, Joost A., et al. "Flux-profile relationships in the atmospheric surface - layer." Journal of the atmospheric Sciences 28.2 (1971): 181-189. - doi: [10.1175/1520-0469(1971)028<0181:FPRITA>2.0.CO;2](https://doi.org/10.1175/1520-0469(1971)028<0181:FPRITA>2.0.CO;2) - - Nishizawa, S., and Y. Kitamura. "A Surface Flux Scheme Based on the Monin-Obukhov Similarity for Finite Volume Models." Journal of Advances in Modeling Earth Systems 10.12 (2018): 3159-3175. doi: [10.1029/2018MS001534](https://doi.org/10.1029/2018MS001534) +- Businger, Joost A., et al. "Flux-profile relationships in the atmospheric surface + layer." Journal of the atmospheric Sciences 28.2 (1971): 181-189. + doi: [10.1175/1520-0469(1971)028<0181:FPRITA>2.0.CO;2](https://doi.org/10.1175/1520-0469(1971)028<0181:FPRITA>2.0.CO;2) + - Byun, Daewon W. "On the analytical solutions of flux-profile relationships for the atmospheric surface layer." Journal of Applied Meteorology 29.7 (1990): 652-657. doi: [10.1175/1520-0450(1990)029<0652:OTASOF>2.0.CO;2](https://doi.org/10.1175/1520-0450(1990)029<0652:OTASOF>2.0.CO;2) diff --git a/src/Common/SurfaceFluxes/SurfaceFluxes.jl b/src/Common/SurfaceFluxes/SurfaceFluxes.jl index f54e37ede59..a955b25212a 100644 --- a/src/Common/SurfaceFluxes/SurfaceFluxes.jl +++ b/src/Common/SurfaceFluxes/SurfaceFluxes.jl @@ -1,19 +1,13 @@ """ SurfaceFluxes - Surface flux functions, e.g., for buoyancy flux, - friction velocity, and exchange coefficients. - -## Sub-modules - - module Byun1990 - - module Nishizawa2018 - ## Interface - - [`compute_buoyancy_flux`](@ref) computes the buoyancy flux - - In addition, each sub-module has the following functions: - - [`monin_obukhov_len`](@ref) computes the Monin-Obukhov length - - [`compute_friction_velocity`](@ref) computes the friction velocity - - [`compute_exchange_coefficients`](@ref) computes the exchange coefficients + - [`surface_conditions`](@ref) computes + - Monin-Obukhov length + - Potential temperature flux (if not given) using Monin-Obukhov theory + - transport fluxes using Monin-Obukhov theory + - friction velocity/temperature scale/tracer scales + - exchange coefficients ## References @@ -37,255 +31,231 @@ pages={652--657}, year={1990} } - """ module SurfaceFluxes using RootSolvers using ..MoistThermodynamics -using CLIMAParameters +using DocStringExtensions +using NLsolve +using CLIMAParameters: AbstractParameterSet using CLIMAParameters.Planet: molmass_ratio, grav +using CLIMAParameters.SubgridScale: von_karman_const -""" - compute_buoyancy_flux(param_set, shf, lhf, T_b, q, α_0) - -Computes buoyancy flux given - - `shf` sensible heat flux - - `lhf` latent heat flux - - `T_b` surface boundary temperature - - `q` phase partition (see [`PhasePartition`](@ref)) - - `α_0` specific volume -""" -function compute_buoyancy_flux( - param_set::AbstractParameterSet, - shf, - lhf, - T_b, - q, - α_0, -) - FT = typeof(shf) - _molmass_ratio::FT = molmass_ratio(param_set) - _grav::FT = grav(param_set) - cp_ = cp_m(param_set, q) - lv = latent_heat_vapor(param_set, T_b) - temp1 = (_molmass_ratio - 1) - temp2 = (shf + temp1 * cp_ * T_b * lhf / lv) - return (_grav * α_0 / cp_ / T_b * temp2) -end - -module Byun1990 +abstract type SurfaceFluxesModel end -using RootSolvers -using ...MoistThermodynamics -using CLIMAParameters -using CLIMAParameters.SubgridScale: von_karman_const +struct Momentum end +struct Heat end -""" Computes ψ_m for stable case. See Eq. 12 """ -ψ_m_stable(ζ, ζ_0, β_m) = -β_m * (ζ - ζ_0) +export surface_conditions -""" Computes ψ_h for stable case. See Eq. 13 """ -ψ_h_stable(ζ, ζ_0, β_h) = -β_h * (ζ - ζ_0) +""" + SurfaceFluxConditions{FT} -""" Computes ψ_m for unstable case. See Eq. 14 """ -function ψ_m_unstable(ζ, ζ_0, γ_m) - x(ζ) = sqrt(sqrt(1 - γ_m * ζ)) - temp(ζ, ζ_0) = log1p((ζ - ζ_0) / (1 + ζ_0)) # log((1 + ζ)/(1 + ζ_0)) - ψ_m = ( - 2 * temp(x(ζ), x(ζ_0)) + temp(x(ζ)^2, x(ζ_0)^2) - 2 * atan(x(ζ)) + - 2 * atan(x(ζ_0)) - ) - return ψ_m -end +Surface flux conditions, returned from `surface_conditions`. -""" Computes ψ_h for unstable case. See Eq. 15 """ -function ψ_h_unstable(ζ, ζ_0, γ_h) - y(ζ) = sqrt(1 - γ_h * ζ) - ψ_h = 2 * log1p((y(ζ) - y(ζ_0)) / (1 + y(ζ_0))) # log((1 + y(ζ))/(1 + y(ζ_0))) - return ψ_h -end +# Fields +$(DocStringExtensions.FIELDS) """ - monin_obukhov_len(param_set, u, flux) +struct SurfaceFluxConditions{FT} + L_MO::FT + pottemp_flux_star::FT + flux::Vector{FT} + x_star::Vector{FT} + K_exchange::Vector{FT} +end -Computes the Monin-Obukhov length (Eq. 3) -""" -function monin_obukhov_len(param_set::AbstractParameterSet, u, flux) - FT = typeof(u) - _von_karman_const::FT = von_karman_const(param_set) - return -u^3 / (flux * _von_karman_const) +function Base.show(io::IO, sfc::SurfaceFluxConditions) + println(io, "----------------------- SurfaceFluxConditions") + println(io, "L_MO = ", sfc.L_MO) + println(io, "pottemp_flux_star = ", sfc.pottemp_flux_star) + println(io, "flux = ", sfc.flux) + println(io, "x_star = ", sfc.x_star) + println(io, "K_exchange = ", sfc.K_exchange) + println(io, "-----------------------") end """ - compute_friction_velocity(param_set, flux, z_0, z_1, β_m, γ_m, tol_abs, iter_max) - -Computes roots of friction velocity equation (Eq. 10) - -`u_ave = u_* ( ln(z/z_0) - ψ_m(z/L, z_0/L) ) /κ` + surface_conditions + +Surface conditions given + - `x_initial` initial guess for solution (`L_MO, u_star, θ_star, ϕ_star, ...`) + - `x_ave` volume-averaged value for variable `x` + - `x_s` surface value for variable `x` + - `z_0` roughness length for variable `x` + - `F_exchange` flux at the top for variable `x` + - `dimensionless_number` dimensionless turbulent transport coefficient: + - Momentum: 1 + - Heat: Turbulent Prantl number at neutral stratification + - Mass: Turbulent Schmidt number + - ... + - `θ_bar` basic potential temperature + - `Δz` layer thickness (not spatial discretization) + - `z` coordinate axis + - `a` free model parameter with prescribed value of 4.7 + - `pottemp_flux_given` potential temperature flux (optional) + +If `pottemp_flux` is not given, then it is computed by iteration +of equations 3, 17, and 18 in Nishizawa2018. """ -function compute_friction_velocity( +function surface_conditions( param_set::AbstractParameterSet, - u_ave::FT, - flux::FT, - z_0::FT, - z_1::FT, - β_m::FT, - γ_m::FT, - tol_abs::FT, - iter_max::IT, -) where {FT, IT} - - _von_karman_const = FT(von_karman_const(param_set)) - - ustar_0 = u_ave * _von_karman_const / log(z_1 / z_0) - ustar = ustar_0 - let u_ave = u_ave, - flux = flux, - z_0 = z_0, - z_1 = z_1, - β_m = β_m, - γ_m = γ_m, - param_set = param_set - - # use neutral condition as first guess - stable = z_1 / monin_obukhov_len(param_set, ustar_0, flux) >= 0 - function compute_ψ_m(u) - L_MO = monin_obukhov_len(param_set, u, flux) - ζ = z_1 / L_MO - ζ_0 = z_0 / L_MO - return stable ? ψ_m_stable(ζ, ζ_0, β_m) : ψ_m_unstable(ζ, ζ_0, γ_m) - end - function compute_u_ave_over_ustar(u) - _von_karman_const = FT(von_karman_const(param_set)) - return (log(z_1 / z_0) - compute_ψ_m(u)) / _von_karman_const # Eq. 10 - end - compute_ustar(u) = u_ave / compute_u_ave_over_ustar(u) - - if (abs(flux) > 0) - ustar_1 = compute_ustar(ustar_0) - sol = RootSolvers.find_zero( - u -> u_ave - u * compute_u_ave_over_ustar(u), - SecantMethod(ustar_0, ustar_1), - CompactSolution(), - SolutionTolerance(tol_abs), - iter_max, - ) - ustar = sol.root + x_initial::Vector{FT}, + x_ave::Vector{FT}, + x_s::Vector{FT}, + z_0::Vector{FT}, + F_exchange::Vector{FT}, + dimensionless_number::Vector{FT}, + θ_bar::FT, + Δz::FT, + z::FT, + a::FT, + pottemp_flux_given::Union{Nothing, FT} = nothing, +) where {FT <: AbstractFloat, AbstractParameterSet} + + n_vars = length(x_initial) - 1 + @assert length(x_initial) == n_vars + 1 + @assert length(x_ave) == n_vars + @assert length(x_s) == n_vars + @assert length(z_0) == n_vars + @assert length(F_exchange) == n_vars + @assert length(dimensionless_number) == n_vars + local sol + let param_set = param_set + function f!(F, x_all) + L_MO, x_vec = x_all[1], x_all[2:end] + u, θ = x_vec[1], x_vec[2] + pottemp_flux = + pottemp_flux_given == nothing ? -u * θ : pottemp_flux_given + F[1] = L_MO - monin_obukhov_len(param_set, u, θ_bar, pottemp_flux) + for i in 1:n_vars + ϕ = x_vec[i] + transport = i == 1 ? Momentum() : Heat() + F[i + 1] = + ϕ - compute_physical_scale( + param_set, + ϕ, + u, + Δz, + a, + θ_bar, + pottemp_flux, + z_0[i], + x_ave[i], + x_s[i], + dimensionless_number[i], + transport, + ) + end end - + sol = nlsolve(f!, x_initial, autodiff = :forward) + end + if converged(sol) + L_MO, x_star = sol.zero[1], sol.zero[2:end] + u_star, θ_star = x_star[1], x_star[2] + else + L_MO, x_star = sol.zero[1], sol.zero[2:end] + u_star, θ_star = x_star[1], x_star[2] + println("Unconverged surface fluxes") end - return ustar -end + _grav::FT = grav(param_set) + _von_karman_const::FT = von_karman_const(param_set) + pottemp_flux_star = -u_star^3 * θ_bar / (_von_karman_const * _grav * L_MO) + flux = -u_star * x_star + K_exchange = compute_exchange_coefficients( + param_set, + z, + F_exchange, + a, + x_star, + θ_bar, + dimensionless_number, + L_MO, + ) -""" - compute_exchange_coefficients(param_set, Ri, z_b, z_0, γ_m, γ_h, β_m, β_h, Pr_0) + return SurfaceFluxConditions( + L_MO, + pottemp_flux_star, + flux, + x_star, + K_exchange, + ) +end -Computes exchange transfer coefficients: - - `C_D` momentum exchange coefficient (Eq. 36) - - `C_H` thermodynamic exchange coefficient (Eq. 37) - - `L_mo` Monin-Obukhov length (re-arranged Eq. 3) -TODO: `Pr_0` should come from CLIMAParameters -""" -function compute_exchange_coefficients( +function compute_physical_scale( param_set::AbstractParameterSet, - Ri::FT, - z_b::FT, - z_0::FT, - γ_m::FT, - γ_h::FT, - β_m::FT, - β_h::FT, - Pr_0::FT, -) where {FT} - logz = log(z_b / z_0) - zfactor = z_b / (z_b - z_0) * logz - s_b = Ri / Pr_0 + x, + u, + Δz, + a, + θ_bar, + pottemp_flux, + z_0, + x_ave, + x_s, + dimensionless_number, + transport, +) + FT = typeof(u) _von_karman_const::FT = von_karman_const(param_set) - if Ri > 0 - temp = ((1 - 2 * β_h * Ri) - sqrt(1 + 4 * (β_h - β_m) * s_b)) - ζ = zfactor / (2 * β_h * (β_m * Ri - 1)) * temp # Eq. 19 - L_mo = z_b / ζ # LHS of Eq. 3 - ζ_0 = z_0 / L_mo - ψ_m = ψ_m_stable(ζ, ζ_0, β_m) - ψ_h = ψ_h_stable(ζ, ζ_0, β_h) - else - Q_b = FT(1 / 9) * (1 / (γ_m^2) + 3 * γ_h / γ_m * s_b^2) # Eq. 31 - P_b = FT(1 / 54) * (-2 / (γ_m^3) + 9 / γ_m * (-γ_h / γ_m + 3) * s_b^2) # Eq. 32 - crit = Q_b^3 - P_b^2 - if crit < 0 - T_b = cbrt(sqrt(-crit) + abs(P_b)) # Eq. 34 - ζ = zfactor * (1 / (3 * γ_m) - (T_b + Q_b / T_b)) # Eq. 29 - else - θ_b = acos(P_b / sqrt(Q_b^3)) # Eq. 33 - ζ = zfactor * (-2 * sqrt(Q_b) * cos(θ_b / 3) + 1 / (3 * γ_m)) # Eq. 28 - end - L_mo = z_b / ζ # LHS of Eq. 3 - ζ_0 = z_0 / L_mo - ψ_m = ψ_m_unstable(ζ, ζ_0, γ_m) - ψ_h = ψ_h_unstable(ζ, ζ_0, γ_h) - end - cu = _von_karman_const / (logz - ψ_m) # Eq. 10, solved for u^* - cth = _von_karman_const / (logz - ψ_h) / Pr_0 # Eq. 11, solved for h^* - C_D = cu^2 # Eq. 36 - C_H = cu * cth # Eq. 37 - return C_D, C_H, L_mo + L = monin_obukhov_len(param_set, u, θ_bar, pottemp_flux) + R_z0 = compute_R_z0(z_0, Δz) + temp1 = log(Δz / z_0) + temp2 = -compute_Psi(Δz / L, L, a, dimensionless_number, transport) + temp3 = + z_0 / Δz * compute_Psi(z_0 / L, L, a, dimensionless_number, transport) + temp4 = + R_z0 * (compute_psi(z_0 / L, L, a, dimensionless_number, transport) - 1) + return (1 / dimensionless_number) * _von_karman_const / + (temp1 + temp2 + temp3 + temp4) * (x_ave - x_s) end -end # Byun1990 module - -module Nishizawa2018 -using RootSolvers -using ...MoistThermodynamics -using CLIMAParameters -using CLIMAParameters.Planet: grav -using CLIMAParameters.SubgridScale: von_karman_const - -""" Computes `R_z0` expression, defined after Eq. 15 """ +""" Computes R_z0 expression, defined after Eq. 15 """ compute_R_z0(z_0, Δz) = 1 - z_0 / Δz -""" Computes `f_m` in Eq. A7 """ +""" Computes f_m in Eq. A7 """ compute_f_m(ζ) = sqrt(sqrt(1 - 15 * ζ)) -""" Computes `f_h` in Eq. A8 """ +""" Computes f_h in Eq. A8 """ compute_f_h(ζ) = sqrt(1 - 9 * ζ) -""" Computes `ψ_m` in Eq. A3 """ -function compute_ψ_m(ζ, L, a) - FT = typeof(ζ) +""" Computes psi_m in Eq. A3 """ +function compute_psi(ζ, L, a, dimensionless_number, ::Momentum) f_m = compute_f_m(ζ) return L >= 0 ? -a * ζ : - log((1 + f_m)^2 * (1 + f_m^2) / 8) - 2 * atan(f_m) + FT(π / 2) + log((1 + f_m)^2 * (1 + f_m^2) / 8) - 2 * atan(f_m) + FT(π) / 2 end -""" Computes `ψ_h` in Eq. A4 """ -compute_ψ_h(ζ, L, a, Pr) = - L >= 0 ? -a * ζ / Pr : 2 * log((1 + compute_f_h(ζ)) / 2) +""" Computes psi_h in Eq. A4 """ +function compute_psi(ζ, L, a, dimensionless_number, ::Heat) + L >= 0 ? -a * ζ / dimensionless_number : 2 * log((1 + compute_f_h(ζ)) / 2) +end -""" Computes `Ψ_m` in Eq. A5 """ -function compute_Ψ_m(ζ, L, a, tol) - FT = typeof(ζ) - if ζ < tol - return ζ >= 0 ? -a * ζ / 2 : -FT(15) * ζ / FT(8) # Computes Ψ_m in Eq. A13 +""" Computes Psi_m in Eq. A5 """ +function compute_Psi(ζ, L, a, dimensionless_number, ::Momentum) + FT = typeof(L) + if abs(ζ) < eps(FT) + return ζ >= 0 ? -a * ζ / 2 : -FT(15) * ζ / FT(8) # Computes Psi_m in Eq. A13 else f_m = compute_f_m(ζ) # Note that "1-f^3" in is a typo, it is supposed to be "1-f_m^3". # This was confirmed by communication with the author. return L >= 0 ? -a * ζ / 2 : - log((1 + f_m)^2 * (1 + f_m^2) / 8) - 2 * atan(f_m) + FT(π / 2) - - 1 + (1 - f_m^3) / (12 * ζ) + log((1 + f_m)^2 * (1 + f_m^2) / 8) - 2 * atan(f_m) + FT(π) / 2 - 1 + + (1 - f_m^3) / (12 * ζ) end end -""" Computes `Ψ_h` in Eq. A6 """ -function compute_Ψ_h(ζ, L, a, Pr, tol) - FT = typeof(ζ) - if ζ < tol - return ζ >= 0 ? -a * ζ / (2 * Pr) : -9 * ζ / 4 # Computes Ψ_h in Eq. A14 +""" Computes Psi_h in Eq. A6 """ +function compute_Psi(ζ, L, a, dimensionless_number, ::Heat) + if abs(ζ) < eps(typeof(L)) + return ζ >= 0 ? -a * ζ / (2 * dimensionless_number) : -9 * ζ / 4 # Computes Psi_h in Eq. A14 else f_h = compute_f_h(ζ) - return L >= 0 ? -a * ζ / (2 * Pr) : + return L >= 0 ? -a * ζ / (2 * dimensionless_number) : 2 * log((1 + f_h) / 2) + 2 * (1 - f_h) / (9 * ζ) end end @@ -293,110 +263,45 @@ end """ monin_obukhov_len(param_set, u, θ, flux) -Computes Monin-Obukhov length. Eq. 3 +Monin-Obukhov length. Eq. 3 """ -function monin_obukhov_len(param_set::AbstractParameterSet, u, θ, flux) +function monin_obukhov_len(param_set::AbstractParameterSet, u, θ_bar, flux) FT = typeof(u) - _von_karman_const::FT = von_karman_const(param_set) _grav::FT = grav(param_set) - return -u^3 * θ / (_von_karman_const * _grav * flux) -end - -""" - compute_friction_velocity(param_set, u_ave, θ, flux, Δz, z_0, a, Ψ_m_tol, tol_abs, iter_max) - -Computes friction velocity, in Eq. 12 in, by solving the -non-linear equation: - -`u_ave = ustar/κ * ( ln(Δz/z_0) - Ψ_m(Δz/L) + z_0/Δz * Ψ_m(z_0/L) + R_z0 [ψ_m(z_0/L) - 1] )` - -where `L` is a non-linear function of `ustar` (see [`monin_obukhov_len`](@ref)). -""" -function compute_friction_velocity( - param_set::AbstractParameterSet, - u_ave::FT, - θ::FT, - flux::FT, - Δz::FT, - z_0::FT, - a::FT, - Ψ_m_tol::FT, - tol_abs::FT, - iter_max::IT, -) where {FT, IT} _von_karman_const::FT = von_karman_const(param_set) - ustar_0 = u_ave * _von_karman_const / log(Δz / z_0) - ustar = ustar_0 - let u_ave = u_ave, - _von_karman_const = _von_karman_const, - θ = θ, - flux = flux, - Δz = Δz, - z_0 = z_0, - a = a, - Ψ_m_tol = Ψ_m_tol, - tol_abs = tol_abs, - iter_max = iter_max - # Note the lowercase psi (ψ) and uppercase psi (Ψ): - Ψ_m_closure(ζ, L) = compute_Ψ_m(ζ, L, a, Ψ_m_tol) - ψ_m_closure(ζ, L) = compute_ψ_m(ζ, L, a) - function compute_u_ave_over_ustar(u) - L = monin_obukhov_len(param_set, u, θ, flux) - R_z0 = compute_R_z0(z_0, Δz) - temp1 = log(Δz / z_0) - temp2 = -Ψ_m_closure(Δz / L, L) - temp3 = z_0 / Δz * Ψ_m_closure(z_0 / L, L) - temp4 = R_z0 * (ψ_m_closure(z_0 / L, L) - 1) - return (temp1 + temp2 + temp3 + temp4) / _von_karman_const - end - compute_ustar(u) = u_ave / compute_u_ave_over_ustar(u) - ustar_1 = compute_ustar(ustar_0) - sol = RootSolvers.find_zero( - u -> u_ave - u * compute_u_ave_over_ustar(u), - SecantMethod(ustar_0, ustar_1), - CompactSolution(), - SolutionTolerance(tol_abs), - iter_max, - ) - ustar = sol.root - end - return ustar + return -u^3 * θ_bar / (_von_karman_const * _grav * flux) end """ - compute_exchange_coefficients(param_set, z, F_m, F_h, a, u_star, θ, flux, Pr) + compute_exchange_coefficients(z, F_exchange, a, x_star, θ_bar, dimensionless_number, L_MO) -Computes exchange transfer coefficients: +Computes exchange transfer coefficients - - `K_D` momentum exchange coefficient - - `K_H` thermodynamic exchange coefficient - - `L_mo` Monin-Obukhov length - -TODO: `Pr` should come from CLIMAParameters + - `K_exchange` exchange coefficients """ function compute_exchange_coefficients( - param_set::AbstractParameterSet, - z::FT, - F_m::FT, - F_h::FT, - a::FT, - u_star::FT, - θ::FT, - flux::FT, - Pr::FT, -) where {FT} + param_set, + z, + F_exchange, + a, + x_star, + θ_bar, + dimensionless_number, + L_MO, +) + N = length(F_exchange) + FT = typeof(z) + K_exchange = Vector{FT}(undef, N) _von_karman_const::FT = von_karman_const(param_set) - L_mo = monin_obukhov_len(param_set, u_star, θ, flux) - ψ_m = compute_ψ_m(z / L_mo, L_mo, a) - ψ_h = compute_ψ_h(z / L_mo, L_mo, a, Pr) - - K_m = -F_m * _von_karman_const * z / (u_star * ψ_m) # Eq. 19 - K_h = -F_h * _von_karman_const * z / (Pr * θ * ψ_h) # Eq. 20 + for i in 1:N + transport = i == 1 ? Momentum() : Heat() + psi = compute_psi(z / L_MO, L_MO, a, dimensionless_number[i], transport) + K_exchange[i] = + -F_exchange[i] * _von_karman_const * z / (x_star[i] * psi) # Eq. 19 in + end - return K_m, K_h, L_mo + return K_exchange end -end # Nishizawa2018 module - end # SurfaceFluxes module diff --git a/test/Common/SurfaceFluxes/runtests.jl b/test/Common/SurfaceFluxes/runtests.jl index 95459b50dd6..afe1f317466 100644 --- a/test/Common/SurfaceFluxes/runtests.jl +++ b/test/Common/SurfaceFluxes/runtests.jl @@ -1,11 +1,8 @@ using Test using ClimateMachine.SurfaceFluxes -using ClimateMachine.SurfaceFluxes.Nishizawa2018 -using ClimateMachine.SurfaceFluxes.Byun1990 -using ClimateMachine.MoistThermodynamics -using RootSolvers -using CLIMAParameters + +using CLIMAParameters: AbstractEarthParameterSet struct EarthParameterSet <: AbstractEarthParameterSet end const param_set = EarthParameterSet() @@ -13,147 +10,50 @@ const param_set = EarthParameterSet() # These tests have been run to ensure they do not fail, # but they need further testing for correctness. -FT = Float32 -rtol = 10 * eps(FT) - @testset "SurfaceFluxes" begin - shf, lhf, T_b, q_pt, α_0 = - FT(60), FT(50), FT(350), PhasePartition{FT}(0.01, 0.002, 0.0001), FT(1) - buoyancy_flux = - SurfaceFluxes.compute_buoyancy_flux(param_set, shf, lhf, T_b, q_pt, α_0) - @test buoyancy_flux ≈ 0.0017808608107074118 - @test buoyancy_flux isa FT -end - -@testset "SurfaceFluxes.Byun1990" begin - u, flux = FT(0.1), FT(0.2) - MO_len = Byun1990.monin_obukhov_len(param_set, u, flux) - @test MO_len ≈ -0.0125 - - u_ave, buoyancy_flux, z_0, z_1 = FT(0.1), FT(0.2), FT(2), FT(5) - γ_m, β_m = FT(15), FT(4.8) - tol_abs, iter_max = FT(1e-3), 10 - u_star = Byun1990.compute_friction_velocity( - param_set, - u_ave, - buoyancy_flux, - z_0, - z_1, - β_m, - γ_m, - tol_abs, - iter_max, - ) - @test u_star ≈ 0.201347256193615 atol = tol_abs - @test u_star isa FT - - - Ri, z_b, z_0, Pr_0 = FT(10), FT(2), FT(5.0), FT(0.74) - γ_m, γ_h, β_m, β_h = FT(15), FT(9), FT(4.8), FT(7.8) - cm, ch, L_mo = Byun1990.compute_exchange_coefficients( - param_set, - Ri, - z_b, - z_0, - γ_m, - γ_h, - β_m, - β_h, - Pr_0, - ) - @test cm ≈ 19.700348427787368 - @test ch ≈ 3.3362564728997803 - @test L_mo ≈ -14.308268023583906 - @test cm isa FT - @test ch isa FT - @test L_mo isa FT - - Ri, z_b, z_0, Pr_0 = FT(-10), FT(10.0), FT(1.0), FT(0.74) - γ_m, γ_h, β_m, β_h = FT(15.0), FT(9.0), FT(4.8), FT(7.8) - cm, ch, L_mo = Byun1990.compute_exchange_coefficients( - param_set, - Ri, - z_b, - z_0, - γ_m, - γ_h, - β_m, - β_h, - Pr_0, - ) - @test cm ≈ 0.33300280321092746 rtol = rtol - @test ch ≈ 1.131830939627489 rtol = rtol - @test L_mo ≈ -0.3726237964444814 rtol = rtol - @test cm isa FT - @test ch isa FT - @test L_mo isa FT -end - -@testset "SurfaceFluxes.Nishizawa2018" begin - FT = Float32 - rtol = 10 * eps(FT) - u, θ, flux = FT(2), FT(350), FT(20) - MO_len = Nishizawa2018.monin_obukhov_len(param_set, u, θ, flux) - @test MO_len ≈ -35.67787971457696 - - u_ave, θ, flux, Δz, z_0, a = - FT(110.0), FT(350.0), FT(20.0), FT(100.0), FT(0.01), FT(5.0) - Ψ_m_tol, tol_abs, iter_max = FT(1e-3), FT(1e-3), 10 - u_star = Nishizawa2018.compute_friction_velocity( - param_set, - u_ave, - θ, - flux, - Δz, - z_0, - a, - Ψ_m_tol, - tol_abs, - iter_max, - ) - @test u_star ≈ 5.526644550864822 atol = tol_abs - @test u_star isa FT - - FT = Float64 - rtol = 10 * eps(FT) - z, F_m, F_h, a, u_star, θ, flux, Pr = - FT(1), FT(2), FT(3), FT(5), FT(110), FT(350), FT(20), FT(0.74) - - K_m, K_h, L_mo = Nishizawa2018.compute_exchange_coefficients( - param_set, - z, - F_m, - F_h, - a, - u_star, - θ, - flux, - Pr, - ) - @test K_m ≈ -11512.071612359368 rtol = rtol - @test K_h ≈ -6111.6196776263805 rtol = rtol - @test L_mo ≈ -5.935907237512742e6 rtol = rtol - - # Test type-stability: FT = Float32 - rtol = 10 * eps(FT) - z, F_m, F_h, a, u_star, θ, flux, Pr = - FT(1), FT(2), FT(3), FT(5), FT(110), FT(350), FT(20), FT(0.74) - - K_m, K_h, L_mo = Nishizawa2018.compute_exchange_coefficients( - param_set, - z, - F_m, - F_h, - a, - u_star, - θ, - flux, - Pr, - ) - @test K_m isa FT - @test K_h isa FT - @test L_mo isa FT + x_initial = FT[100, 15.0, 350.0] + F_exchange = FT[2.0, 3.0] + z_0 = FT[1.0, 1.0] + dimensionless_number = FT[1.0, 0.74] + x_ave = FT[5.0, 350.0] + x_s = FT[0.0, 300.0] + + Δz = FT(2.0) + z = FT(0.5) + θ_bar = FT(300.0) + a = FT(4.7) + pottemp_flux_given = FT(-200.0) + args = x_initial, + x_ave, + x_s, + z_0, + F_exchange, + dimensionless_number, + θ_bar, + Δz, + z, + a, + pottemp_flux_given + sfc = surface_conditions(param_set, args[1:(end - 1)]...) + @test sfc.L_MO isa FT + @test sfc.pottemp_flux_star isa FT + + # @test sfc.parameters == FT + + @test sfc.L_MO ≈ 54.563405359719404 + @test sfc.pottemp_flux_star ≈ -1132.9097989525164 + @test all(sfc.flux .≈ [-86.79000158448329, -1132.9097989138224]) + @test all(sfc.x_star .≈ [9.316115155175106, 121.60753490520031]) + @test all(sfc.K_exchange .≈ [0.9969164175880834, 0.08477271454000379]) + + sfc = surface_conditions(param_set, args...) + + @test sfc.L_MO ≈ 405.8862509767147 + @test sfc.pottemp_flux_star ≈ -199.99999999999997 + @test all(sfc.flux .≈ [-104.07858678549196, -1399.2078659154145]) + @test all(sfc.x_star .≈ [10.20189133374258, 137.15181039887756]) + @test all(sfc.K_exchange .≈ [6.771981702484999, 0.5591365770421184]) end