diff --git a/Manifest.toml b/Manifest.toml index 7a74b9f7652..db06ce21be7 100644 --- a/Manifest.toml +++ b/Manifest.toml @@ -44,9 +44,9 @@ version = "3.5.0+3" [[ArrayInterface]] deps = ["LinearAlgebra", "Requires", "SparseArrays"] -git-tree-sha1 = "649c08a5a3a513f4662673d3777fe6ccb4df9f5d" +git-tree-sha1 = "066d1e7a9eb4873660791db7f0d8c7902600b81c" uuid = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9" -version = "2.8.7" +version = "2.11.0" [[ArrayLayouts]] deps = ["FillArrays", "LinearAlgebra"] @@ -99,15 +99,9 @@ version = "0.3.3" [[ChainRulesCore]] deps = ["MuladdMacro"] -git-tree-sha1 = "32e2c6e44d4fdd985b5688b5e85c1f6892cf3d15" +git-tree-sha1 = "c384e0e4fe6bfeb6bec0d41f71cc5e391cd110ba" uuid = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4" -version = "0.8.0" - -[[CodeTracking]] -deps = ["InteractiveUtils", "UUIDs"] -git-tree-sha1 = "9c173f62af93cce8af2bd3527d160b6ddd6eaf81" -uuid = "da1fd8a2-8d9e-5ec2-8556-3022fb5608a2" -version = "1.0.0" +version = "0.8.1" [[CodecZlib]] deps = ["TranscodingStreams", "Zlib_jll"] @@ -128,9 +122,9 @@ version = "0.3.0" [[Compat]] deps = ["Base64", "Dates", "DelimitedFiles", "Distributed", "InteractiveUtils", "LibGit2", "Libdl", "LinearAlgebra", "Markdown", "Mmap", "Pkg", "Printf", "REPL", "Random", "SHA", "Serialization", "SharedArrays", "Sockets", "SparseArrays", "Statistics", "Test", "UUIDs", "Unicode"] -git-tree-sha1 = "054993b6611376ddb40203e973e954fd9d1d1902" +git-tree-sha1 = "a6a8197ae253f2c1a22b2ae17c2dfaf5812c03aa" uuid = "34da2185-b29b-5c13-b0c7-acf172513d20" -version = "3.12.0" +version = "3.13.0" [[CompilerSupportLibraries_jll]] deps = ["Libdl", "Pkg"] @@ -174,12 +168,6 @@ git-tree-sha1 = "f0464e499ab9973b43c20f8216d088b61fda80c6" uuid = "adafc99b-e345-5852-983c-f28acb93d879" version = "0.2.2" -[[Cthulhu]] -deps = ["CodeTracking", "InteractiveUtils", "REPL", "UUIDs", "Unicode"] -git-tree-sha1 = "3f4601d6ec0e967a24a9926ea7a981c7fe63e9a9" -uuid = "f68482b8-f384-11e8-15f7-abe071a5a75f" -version = "1.1.2" - [[DataAPI]] git-tree-sha1 = "176e23402d80e7743fc26c19c681bfb11246af32" uuid = "9a962f9c-6df0-11e9-0e5d-c546b8b5ee8a" @@ -252,10 +240,10 @@ uuid = "497a8b3b-efae-58df-a0af-a86822472b78" version = "1.1.12" [[ExponentialUtilities]] -deps = ["LinearAlgebra", "Printf", "SparseArrays"] -git-tree-sha1 = "1672dedeacaab85345fd359ad56dde8fb5d48a45" +deps = ["LinearAlgebra", "Printf", "Requires", "SparseArrays"] +git-tree-sha1 = "91f7498b66205431fe3e35833cda97a22b1ab6a5" uuid = "d4d017d3-3776-5f7e-afef-a10c40355c18" -version = "1.6.0" +version = "1.7.0" [[ExprTools]] git-tree-sha1 = "6f0517056812fd6aa3af23d4b70d5325a2ae4e95" @@ -282,9 +270,9 @@ version = "0.8.13" [[FiniteDiff]] deps = ["ArrayInterface", "LinearAlgebra", "Requires", "SparseArrays", "StaticArrays"] -git-tree-sha1 = "fec7c2cb45c27071ef487fa7cae4fcac7509aa10" +git-tree-sha1 = "43b397bf66d07d113cc2cc012dd2029d336c3894" uuid = "6a86dc24-6348-571c-b903-95158fe2bd41" -version = "2.3.2" +version = "2.5.1" [[Formatting]] deps = ["Printf"] @@ -335,9 +323,9 @@ version = "0.4.0" [[HTTP]] deps = ["Base64", "Dates", "IniFile", "MbedTLS", "Sockets"] -git-tree-sha1 = "eca61b35cdd8cd2fcc5eec1eda766424a995b02f" +git-tree-sha1 = "2ac03263ce44be4222342bca1c51c36ce7566161" uuid = "cd3eb016-35fb-5094-929b-558a96fad6f3" -version = "0.8.16" +version = "0.8.17" [[Inflate]] git-tree-sha1 = "f5fc07d4e706b84f72d54eedcc1c13d92fb0871c" @@ -356,9 +344,9 @@ uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240" [[Intervals]] deps = ["Dates", "Printf", "RecipesBase", "Serialization", "TimeZones"] -git-tree-sha1 = "8ffd41e6ddb0b569744aac36f8f5a4c92d28fc9e" +git-tree-sha1 = "f8fd8065b6bccfbc2f9e7dd1fa5c0dde81d09185" uuid = "d8418881-c3e1-53bb-8760-2df7ec849ed5" -version = "1.4.0" +version = "1.4.1" [[IterativeSolvers]] deps = ["LinearAlgebra", "Printf", "Random", "RecipesBase", "SparseArrays"] @@ -395,17 +383,17 @@ git-tree-sha1 = "d9c6e1efcaa6c2fcd043da812a62b3e489a109a3" uuid = "929cbde3-209d-540e-8aea-75f648917ca0" version = "1.7.0" +[[LabelledArrays]] +deps = ["ArrayInterface", "LinearAlgebra", "MacroTools", "StaticArrays"] +git-tree-sha1 = "5e04374019448f8509349948ab504f117e3b575a" +uuid = "2ee39098-c373-598a-b85f-a56591580800" +version = "1.3.0" + [[LambertW]] git-tree-sha1 = "2d9f4009c486ef676646bca06419ac02061c088e" uuid = "984bce1d-4616-540c-a9ee-88d1112d94c9" version = "0.4.5" -[[LabelledArrays]] -deps = ["ArrayInterface", "LinearAlgebra", "MacroTools", "StaticArrays"] -git-tree-sha1 = "f6def2c9c88908fdde81b9a39c236cf69de94450" -uuid = "2ee39098-c373-598a-b85f-a56591580800" -version = "1.2.2" - [[LazyArrays]] deps = ["ArrayLayouts", "FillArrays", "LinearAlgebra", "MacroTools", "MatrixFactorizations", "StaticArrays"] git-tree-sha1 = "eb02bd7606db9516205b400da8ba01e829ed5f7b" @@ -463,15 +451,16 @@ version = "2.5.0" uuid = "56ddb016-857b-54e1-b83d-db4d58db5568" [[LoggingExtras]] -git-tree-sha1 = "b60616c70eff0cc2c0831b6aace75940aeb0939d" +deps = ["Dates"] +git-tree-sha1 = "03289aba73c0abc25ff0229bed60f2a4129cd15c" uuid = "e6f89c97-d47a-5376-807f-9c37f3926c36" -version = "0.4.1" +version = "0.4.2" [[LoopVectorization]] deps = ["DocStringExtensions", "LinearAlgebra", "OffsetArrays", "SIMDPirates", "SLEEFPirates", "UnPack", "VectorizationBase"] -git-tree-sha1 = "c130c0ec960fa489dc87041d60834c459b307d0c" +git-tree-sha1 = "25b1363081880dcd6a9e27736f4b1951ea4b8933" uuid = "bdcacae8-1622-11e9-2a5c-532679323890" -version = "0.8.18" +version = "0.8.20" [[MPI]] deps = ["Distributed", "DocStringExtensions", "Libdl", "MPICH_jll", "MicrosoftMPI_jll", "OpenMPI_jll", "Pkg", "Random", "Requires", "Serialization", "Sockets"] @@ -545,12 +534,6 @@ git-tree-sha1 = "271e4670a6b30ec1e5daf6ba8b30240b7d85f535" uuid = "85f8d34a-cbdd-5861-8df4-14fed0d494ab" version = "0.10.2" -[[NNPACK_jll]] -deps = ["Libdl", "Pkg"] -git-tree-sha1 = "c3d1a616362645754b18e12dbba96ec311b0867f" -uuid = "a6bfbf70-4841-5cb9-aa18-3a8ad3c413ee" -version = "2018.6.22+0" - [[NLSolversBase]] deps = ["DiffResults", "Distributed", "FiniteDiff", "ForwardDiff"] git-tree-sha1 = "7c4e66c47848562003250f28b579c584e55becc0" @@ -563,6 +546,12 @@ git-tree-sha1 = "ea172c86745810136d744fc67650d2e2de669c4f" uuid = "2774e3e8-f4cf-5e23-947b-6d7e65073b56" version = "4.4.0" +[[NNPACK_jll]] +deps = ["Libdl", "Pkg"] +git-tree-sha1 = "c3d1a616362645754b18e12dbba96ec311b0867f" +uuid = "a6bfbf70-4841-5cb9-aa18-3a8ad3c413ee" +version = "2018.6.22+0" + [[NNlib]] deps = ["Libdl", "LinearAlgebra", "LoopVectorization", "NNPACK_jll", "Pkg", "Requires", "Statistics"] git-tree-sha1 = "1d8128735fdf3ab1643dd8bc9499e4b34ccb718d" @@ -575,9 +564,9 @@ uuid = "77ba4419-2d1f-58cd-9bb1-8ffee604a2e3" version = "0.3.4" [[OffsetArrays]] -git-tree-sha1 = "4ba4cd84c88df8340da1c3e2d8dcb9d18dd1b53b" +git-tree-sha1 = "2066e16af994955287f2e03ba1d9e890eb43b0dd" uuid = "6fe1bfb0-de20-5000-8ca7-80f57d26f881" -version = "1.1.1" +version = "1.1.2" [[OpenBLAS_jll]] deps = ["CompilerSupportLibraries_jll", "Libdl", "Pkg"] @@ -648,15 +637,15 @@ uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7" [[ProgressLogging]] deps = ["Logging", "SHA", "UUIDs"] -git-tree-sha1 = "e47914361d124d8760f5e356403daeb7f3b81633" +git-tree-sha1 = "59398022b661b6fd569f25de6b18fde39843196a" uuid = "33c8b6b6-d38a-422a-b730-caa89a2f386c" -version = "0.1.2" +version = "0.1.3" [[ProgressMeter]] deps = ["Distributed", "Printf"] -git-tree-sha1 = "3e1784c27847bba115815d4d4e668b99873985e5" +git-tree-sha1 = "2de4cddc0ceeddafb6b143b5b6cd9c659b64507c" uuid = "92933f4c-e287-5a05-a399-4b506db050ca" -version = "1.3.1" +version = "1.3.2" [[QuadGK]] deps = ["DataStructures", "LinearAlgebra"] @@ -679,27 +668,27 @@ deps = ["Serialization"] uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" [[RecipesBase]] -git-tree-sha1 = "54f8ceb165a0f6d083f0d12cb4996f5367c6edbc" +git-tree-sha1 = "58de8f7e33b7fda6ee39eff65169cd1e19d0c107" 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" +version = "1.0.2" [[RecursiveArrayTools]] deps = ["ArrayInterface", "LinearAlgebra", "RecipesBase", "Requires", "StaticArrays", "Statistics", "ZygoteRules"] -git-tree-sha1 = "96e71928efa701fa5a6df0f88b51f05ceed70f2c" +git-tree-sha1 = "0ffe36b65f0fc4967a42a673c1a9ffa65724dee6" uuid = "731186ca-8d62-57ce-b412-fbd966d074cd" -version = "2.4.4" +version = "2.5.0" [[RecursiveFactorization]] -deps = ["LinearAlgebra", "LoopVectorization"] -git-tree-sha1 = "09217cb106dd826de9960986207175b52e3035f2" +deps = ["LinearAlgebra", "LoopVectorization", "VectorizationBase"] +git-tree-sha1 = "4ca0bdad1d69abbd59c35af89a9a2ab6cd5ef0f1" uuid = "f2c3362d-daeb-58d1-803e-2bc74f2840b4" -version = "0.1.2" +version = "0.1.4" + +[[Reexport]] +deps = ["Pkg"] +git-tree-sha1 = "7b1d07f411bc8ddb7977ec7f377b97b158514fe0" +uuid = "189a3867-3050-52da-a836-e630ba90ab69" +version = "0.2.0" [[Requires]] deps = ["UUIDs"] @@ -727,24 +716,24 @@ version = "0.2.0" [[Roots]] deps = ["Printf"] -git-tree-sha1 = "c2f7348c55d1433d1cab0159b4d2c6d27af36fc4" +git-tree-sha1 = "069e68c2173b4e4d0c37ffb3268d37f168ad719c" uuid = "f2b01f46-fcfa-551c-844a-d8ac1e96c665" -version = "1.0.2" +version = "1.0.4" [[SHA]] uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce" [[SIMDPirates]] deps = ["VectorizationBase"] -git-tree-sha1 = "cbe8797ac354d0b1dfe75d983429938db834b706" +git-tree-sha1 = "19880eef12759d7d049516fcb11f1ed8440963d6" uuid = "21efa798-c60a-11e8-04d3-e1a92915a26a" -version = "0.8.16" +version = "0.8.21" [[SLEEFPirates]] deps = ["Libdl", "SIMDPirates", "VectorizationBase"] -git-tree-sha1 = "c750d618b7c8268a97e55c70e8c88e56080d30fa" +git-tree-sha1 = "67ae90a18aa8c22bf159318300e765fbd89ddf6e" uuid = "476501e8-09a2-5ece-8869-fb82de89a1fa" -version = "0.5.4" +version = "0.5.5" [[Serialization]] uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b" @@ -774,9 +763,9 @@ uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" [[SparseDiffTools]] deps = ["Adapt", "ArrayInterface", "Compat", "DataStructures", "FiniteDiff", "ForwardDiff", "LightGraphs", "LinearAlgebra", "Requires", "SparseArrays", "VertexSafeGraphs"] -git-tree-sha1 = "bfe68e0d914952932594b3c838f08463b0841037" +git-tree-sha1 = "93666e93899d995ec37abddde4811f533e49c074" uuid = "47a9eef4-7e08-11e9-0b38-333d64bd3804" -version = "1.8.0" +version = "1.9.1" [[SpecialFunctions]] deps = ["OpenSpecFun_jll"] @@ -824,9 +813,9 @@ version = "1.0.0" [[TerminalLoggers]] deps = ["LeftChildRightSiblingTrees", "Logging", "Markdown", "Printf", "ProgressLogging", "UUIDs"] -git-tree-sha1 = "8c05be75dfe73d90e5dfb6293e0c852013f7282d" +git-tree-sha1 = "cbea752b5eef52a3e1188fb31580c3e4fa0cbc35" uuid = "5d786b92-1e48-4d6f-9151-6b4477ca9bed" -version = "0.1.1" +version = "0.1.2" [[Test]] deps = ["Distributed", "InteractiveUtils", "Logging", "Random"] @@ -881,9 +870,9 @@ uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5" [[VectorizationBase]] deps = ["CpuId", "Libdl", "LinearAlgebra"] -git-tree-sha1 = "95c0c737c307dfd4f65ad50a79856b343fdb7959" +git-tree-sha1 = "46afc4dc02f534e149149173f38f8f59cc23c165" uuid = "3d5dd08c-fd9d-11e8-17fa-ed2836048c2f" -version = "0.12.24" +version = "0.12.29" [[VersionParsing]] git-tree-sha1 = "80229be1f670524750d905f8fc8148e5a8c4537f" diff --git a/Project.toml b/Project.toml index 5a156e95c36..3d0b6e3fb41 100644 --- a/Project.toml +++ b/Project.toml @@ -30,6 +30,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" OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" PackageCompiler = "9b87118b-4619-50d2-8e1e-99f35a4d4d9d" diff --git a/docs/src/APIs/Common/SurfaceFluxes.md b/docs/src/APIs/Common/SurfaceFluxes.md index 7ea5cb9406a..290641bac10 100644 --- a/docs/src/APIs/Common/SurfaceFluxes.md +++ b/docs/src/APIs/Common/SurfaceFluxes.md @@ -1,23 +1,11 @@ # Surface Fluxes -```@meta -CurrentModule = ClimateMachine.SurfaceFluxes -``` - -## Surface Fluxes - ```@docs -SurfaceFluxes +ClimateMachine.SurfaceFluxes ``` -## Methods +## API ```@docs -compute_buoyancy_flux -Byun1990.monin_obukhov_len -Byun1990.compute_friction_velocity -Byun1990.compute_exchange_coefficients -Nishizawa2018.monin_obukhov_len -Nishizawa2018.compute_friction_velocity -Nishizawa2018.compute_exchange_coefficients +ClimateMachine.SurfaceFluxes.surface_conditions ``` diff --git a/docs/src/Theory/Common/SurfaceFluxes.md b/docs/src/Theory/Common/SurfaceFluxes.md index 6075ec0955a..04998ec0dc5 100644 --- a/docs/src/Theory/Common/SurfaceFluxes.md +++ b/docs/src/Theory/Common/SurfaceFluxes.md @@ -1,95 +1,18 @@ # 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 SurfaceFluxes.surface_conditions). -## `Nishizawa2018` - -### Plots - -```@example nishizawa2018 -using ClimateMachine.SurfaceFluxes.Nishizawa2018 -using Plots, LaTeXStrings -using CLIMAParameters -struct EarthParameterSet <: AbstractEarthParameterSet end -param_set = EarthParameterSet() - -FT = Float64 - -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 -``` -![](friction_velocity.svg) +## Interface + - [`surface_conditions`](@ref SurfaceFluxes.surface_conditions) 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 @@ -100,6 +23,10 @@ nothing # hide Similarity for Finite Volume Models." Journal of Advances in Modeling Earth Systems 10.12 (2018): 3159-3175. +- 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. diff --git a/src/Common/SurfaceFluxes/SurfaceFluxes.jl b/src/Common/SurfaceFluxes/SurfaceFluxes.jl index af9b22b881d..cbe4c888d75 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` computes the Monin-Obukhov length - - `compute_friction_velocity` computes the friction velocity - - `compute_exchange_coefficients` 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,223 @@ pages={652--657}, year={1990} } - """ module SurfaceFluxes -using RootSolvers using ..Thermodynamics +using KernelAbstractions: @print +using DocStringExtensions +using NLsolve +using CLIMAParameters: AbstractParameterSet using CLIMAParameters using CLIMAParameters.Planet: molmass_ratio, grav - -""" - 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 - -using RootSolvers -using ...Thermodynamics -using CLIMAParameters using CLIMAParameters.SubgridScale: von_karman_const -""" Computes ψ_m for stable case. See Eq. 12 """ -ψ_m_stable(ζ, ζ_0, β_m) = -β_m * (ζ - ζ_0) - -""" Computes ψ_h for stable case. See Eq. 13 """ -ψ_h_stable(ζ, ζ_0, β_h) = -β_h * (ζ - ζ_0) +const APS = AbstractParameterSet +abstract type SurfaceFluxesModel end -""" 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 +struct Momentum end +struct Heat end -""" 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 +export surface_conditions """ - monin_obukhov_len(param_set, u, flux) + SurfaceFluxConditions{FT, VFT} + +Surface flux conditions, returned from `surface_conditions`. -Computes the Monin-Obukhov length (Eq. 3) +# Fields + +$(DocStringExtensions.FIELDS) """ -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) +struct SurfaceFluxConditions{FT, VFT} + L_MO::FT + pottemp_flux_star::FT + flux::VFT + x_star::VFT + K_exchange::VFT 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( - 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 +function surface_conditions( + param_set::APS, + 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, APS <: 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] + @print("Non-converged 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( - 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 +function compute_physical_scale( + param_set::APS, + 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) + term1 = log(Δz / z_0) + term2 = -compute_Psi(Δz / L, L, a, dimensionless_number, transport) + term3 = + z_0 / Δz * compute_Psi(z_0 / L, L, a, dimensionless_number, transport) + term4 = + R_z0 * (compute_psi(z_0 / L, L, a, dimensionless_number, transport) - 1) + return (1 / dimensionless_number) * _von_karman_const / + (term1 + term2 + term3 + term4) * (x_ave - x_s) end -end # Byun1990 module - -module Nishizawa2018 -using RootSolvers -using ...Thermodynamics -using CLIMAParameters -using CLIMAParameters.Planet: grav -using CLIMAParameters.SubgridScale: von_karman_const - """ 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) - + 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 +255,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::APS, 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) - -Computes exchange transfer coefficients: + compute_exchange_coefficients(z, F_exchange, a, x_star, θ_bar, dimensionless_number, L_MO) - - `K_D` momentum exchange coefficient - - `K_H` thermodynamic exchange coefficient - - `L_mo` Monin-Obukhov length +Computes exchange transfer coefficients -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 e82dad04f3f..2f4199e55fd 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.Thermodynamics -using RootSolvers -using CLIMAParameters + +using CLIMAParameters: AbstractEarthParameterSet struct EarthParameterSet <: AbstractEarthParameterSet end const param_set = EarthParameterSet() @@ -13,147 +10,48 @@ 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.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