From 9a075951976f2cc82998ace00f09a735d6a786e7 Mon Sep 17 00:00:00 2001 From: Charles Kawczynski Date: Thu, 9 Jan 2020 09:44:43 -0800 Subject: [PATCH] Improve SurfaceFluxes --- Manifest.toml | 91 +++-- Project.toml | 1 + docs/Project.toml | 1 + docs/src/Atmos/SurfaceFluxes.md | 91 +---- .../SurfaceFluxes/SurfaceFluxes.jl | 360 +++++++----------- .../SurfaceFluxes/runtests.jl | 93 ++--- 6 files changed, 239 insertions(+), 398 deletions(-) diff --git a/Manifest.toml b/Manifest.toml index fb535d75f9a..c65cf2838b5 100644 --- a/Manifest.toml +++ b/Manifest.toml @@ -26,9 +26,9 @@ version = "0.3.2" [[ArrayInterface]] deps = ["LinearAlgebra", "Requires", "SparseArrays"] -git-tree-sha1 = "487fc79b7b77cef8391675c0000ccd8e933e3533" +git-tree-sha1 = "656fd4bcdf204ea96945d9bc1068c0056013438a" uuid = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9" -version = "2.0.0" +version = "2.3.1" [[ArrayLayouts]] deps = ["FillArrays", "LinearAlgebra"] @@ -93,9 +93,14 @@ uuid = "49dc2e85-a5d0-5ad3-a950-438e2897f1b9" version = "0.5.1" [[Cassette]] -git-tree-sha1 = "da85d135b6048d3e78603e277cf9a4609f7e0673" +git-tree-sha1 = "36bd4e0088652b0b2d25a03e531f0d04258feb78" uuid = "7057c7e9-c182-5462-911a-8362d720325c" -version = "0.2.6" +version = "0.3.0" + +[[ChainRulesCore]] +git-tree-sha1 = "8f08bb658bdf59f311a416ca3b3c765a221aee8e" +uuid = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4" +version = "0.4.0" [[CodecZlib]] deps = ["BinaryProvider", "Libdl", "TranscodingStreams"] @@ -146,9 +151,9 @@ version = "1.1.0" [[DataStructures]] deps = ["InteractiveUtils", "OrderedCollections"] -git-tree-sha1 = "a1b652fb77ae8ca7ea328fa7ba5aa151036e5c10" +git-tree-sha1 = "f784254f428fb8fd7ac15982e5862a38a44523d3" uuid = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8" -version = "0.17.6" +version = "0.17.7" [[Dates]] deps = ["Printf"] @@ -156,9 +161,9 @@ uuid = "ade2ca70-3891-5945-98fb-dc099432e06a" [[DelayDiffEq]] deps = ["DataStructures", "DiffEqBase", "LinearAlgebra", "Logging", "OrdinaryDiffEq", "Parameters", "Printf", "RecursiveArrayTools", "Reexport", "Roots"] -git-tree-sha1 = "afa98d9b36df419dbf9586e8f3aa0bf1f4f7cf1d" +git-tree-sha1 = "c88083ef22a570fc155802fc30e21e0d0112d28c" uuid = "bcd4f6db-9728-5f36-b5f7-82caef46ccdb" -version = "5.18.0" +version = "5.19.0" [[DelimitedFiles]] deps = ["Mmap"] @@ -171,10 +176,10 @@ uuid = "39dd38d3-220a-591b-8e3c-4c3a8c710a94" version = "0.4.1" [[DiffEqBase]] -deps = ["ArrayInterface", "Compat", "DataStructures", "DiffEqDiffTools", "Distributed", "DocStringExtensions", "FunctionWrappers", "IterativeSolvers", "IteratorInterfaceExtensions", "LinearAlgebra", "MuladdMacro", "Parameters", "Printf", "RecipesBase", "RecursiveArrayTools", "RecursiveFactorization", "Requires", "Roots", "SparseArrays", "StaticArrays", "Statistics", "SuiteSparse", "TableTraits", "TreeViews"] -git-tree-sha1 = "a484bf139eb4e15df47b2ef82f6e2af0fbb4ad01" +deps = ["ArrayInterface", "ChainRulesCore", "Compat", "DataStructures", "DiffEqDiffTools", "Distributed", "DocStringExtensions", "FunctionWrappers", "IterativeSolvers", "IteratorInterfaceExtensions", "LinearAlgebra", "MuladdMacro", "Parameters", "Printf", "RecipesBase", "RecursiveArrayTools", "RecursiveFactorization", "Requires", "Roots", "SparseArrays", "StaticArrays", "Statistics", "SuiteSparse", "TableTraits", "TreeViews", "ZygoteRules"] +git-tree-sha1 = "2bb222b5dd4b138a3ce285693b1360c7b04e3fb8" uuid = "2b5f629d-d688-5b77-993f-72d75c75574e" -version = "6.7.0" +version = "6.10.2" [[DiffEqCallbacks]] deps = ["DataStructures", "DiffEqBase", "ForwardDiff", "LinearAlgebra", "NLsolve", "OrdinaryDiffEq", "RecipesBase", "RecursiveArrayTools", "StaticArrays"] @@ -184,9 +189,9 @@ version = "2.10.0" [[DiffEqDiffTools]] deps = ["ArrayInterface", "LinearAlgebra", "Requires", "SparseArrays", "StaticArrays"] -git-tree-sha1 = "6d83f9b2c2a552bf5ce29c20a526c0cbfd9e5270" +git-tree-sha1 = "a4ed8a740484627ea41b47f7e1a25dd909a28353" uuid = "01453d9d-ee7c-5054-8395-0335cb756afa" -version = "1.5.0" +version = "1.7.0" [[DiffEqFinancial]] deps = ["DiffEqBase", "DiffEqNoiseProcess", "LinearAlgebra", "Markdown", "RandomNumbers"] @@ -195,22 +200,22 @@ uuid = "5a0ffddc-d203-54b0-88ba-2c03c0fc2e67" version = "2.2.1" [[DiffEqJump]] -deps = ["Compat", "DataStructures", "DiffEqBase", "FunctionWrappers", "LinearAlgebra", "Parameters", "PoissonRandom", "Random", "RandomNumbers", "RecursiveArrayTools", "TreeViews"] -git-tree-sha1 = "d556c231a60f331a2468ee454f3b05afbb5a454b" +deps = ["Compat", "DataStructures", "DiffEqBase", "FunctionWrappers", "LinearAlgebra", "Parameters", "PoissonRandom", "Random", "RandomNumbers", "RecursiveArrayTools", "Statistics", "TreeViews"] +git-tree-sha1 = "b76e5511c2b6675445abb3b35892c18ade9abf74" uuid = "c894b116-72e5-5b58-be3c-e6d8d4ac2b12" -version = "6.3.0" +version = "6.4.0" [[DiffEqNoiseProcess]] deps = ["DataStructures", "DiffEqBase", "LinearAlgebra", "Random", "RandomNumbers", "RecipesBase", "RecursiveArrayTools", "Requires", "ResettableStacks", "StaticArrays", "Statistics"] -git-tree-sha1 = "02eb7e2e202f23a5834727e05eeac40422f0165c" +git-tree-sha1 = "0641d5d6a1dda2a871698155a50c567526974040" uuid = "77a26b50-5914-5dd7-bc55-306e6241c503" -version = "3.6.0" +version = "3.7.0" [[DiffEqPhysics]] deps = ["DiffEqBase", "DiffEqCallbacks", "ForwardDiff", "LinearAlgebra", "Printf", "Random", "RecipesBase", "RecursiveArrayTools", "Reexport", "StaticArrays"] -git-tree-sha1 = "3d7f15f2805bc2f95cfb06d10a920746d11c76d7" +git-tree-sha1 = "abc7d7767994ced55cd2b0caf3057548a8702021" uuid = "055956cb-9e8b-5191-98cc-73ae4a59e68a" -version = "3.3.0" +version = "3.4.0" [[DiffResults]] deps = ["Compat", "StaticArrays"] @@ -371,9 +376,9 @@ version = "0.21.0" [[LLVM]] deps = ["CEnum", "Libdl", "Printf", "Unicode"] -git-tree-sha1 = "74fe444b8b6d1ac01d639b2f9eaf395bcc2e24fc" +git-tree-sha1 = "1d08d7e4250f452f6cb20e4574daaebfdbee0ff7" uuid = "929cbde3-209d-540e-8aea-75f648917ca0" -version = "1.3.2" +version = "1.3.3" [[LazyArrays]] deps = ["ArrayLayouts", "FillArrays", "LinearAlgebra", "MacroTools", "StaticArrays"] @@ -457,9 +462,9 @@ version = "0.2.1" [[MultiScaleArrays]] deps = ["DiffEqBase", "RecursiveArrayTools", "Statistics", "StochasticDiffEq", "TreeViews"] -git-tree-sha1 = "cf6bc444a74ddcb781b04d0be8988bbff83c91df" +git-tree-sha1 = "c69e4a0db8d482a825c7ac76a433759ade2a5af2" uuid = "f9640e96-87f6-5992-9c3b-0743c6a49ffa" -version = "1.5.0" +version = "1.6.0" [[NLSolversBase]] deps = ["Calculus", "DiffEqDiffTools", "DiffResults", "Distributed", "ForwardDiff"] @@ -474,10 +479,10 @@ uuid = "2774e3e8-f4cf-5e23-947b-6d7e65073b56" version = "4.2.0" [[NNlib]] -deps = ["Libdl", "LinearAlgebra", "Requires", "Statistics", "TimerOutputs"] -git-tree-sha1 = "0c667371391fc6bb31f7f12f96a56a17098b3de8" +deps = ["BinaryProvider", "Libdl", "LinearAlgebra", "Requires", "Statistics"] +git-tree-sha1 = "135c0de4794d5e214b06f1fb4787af4a72896e61" uuid = "872c559c-99b0-510c-b3b7-b6c96a88d5cd" -version = "0.6.0" +version = "0.6.2" [[NaNMath]] git-tree-sha1 = "928b8ca9b2791081dc71a51c55347c27c618760f" @@ -492,9 +497,9 @@ version = "1.1.0" [[OrdinaryDiffEq]] deps = ["ArrayInterface", "DataStructures", "DiffEqBase", "DiffEqDiffTools", "ExponentialUtilities", "ForwardDiff", "GenericSVD", "LinearAlgebra", "Logging", "MacroTools", "MuladdMacro", "Parameters", "RecursiveArrayTools", "Reexport", "SparseArrays", "SparseDiffTools", "StaticArrays"] -git-tree-sha1 = "d75c2b12afed168d5ebd60170ebbf46746beac77" +git-tree-sha1 = "143b4d228d3a2e84c45a4acf937f1b0d6788da65" uuid = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" -version = "5.26.4" +version = "5.26.8" [[PDMats]] deps = ["Arpack", "LinearAlgebra", "SparseArrays", "SuiteSparse", "Test"] @@ -556,9 +561,9 @@ uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" [[RandomNumbers]] deps = ["Random", "Requires"] -git-tree-sha1 = "1417be19c15706c1584d01e32662eb640a4cc908" +git-tree-sha1 = "441e6fc35597524ada7f85e13df1f4e10137d16f" uuid = "e6cf234a-135c-5ec9-84dd-332b85af5143" -version = "1.3.0" +version = "1.4.0" [[RecipesBase]] git-tree-sha1 = "7bdce29bc9b2f5660a6e5e64d64d91ec941f6aa2" @@ -567,9 +572,9 @@ version = "0.7.0" [[RecursiveArrayTools]] deps = ["ArrayInterface", "RecipesBase", "Requires", "StaticArrays", "Statistics"] -git-tree-sha1 = "cb0fcf68e2e19e76b84c892b22887f02f10f3d9a" +git-tree-sha1 = "35a01ee8529e84d1ca4c7131fee30e7f61047fc8" uuid = "731186ca-8d62-57ce-b412-fbd966d074cd" -version = "1.2.0" +version = "1.2.1" [[RecursiveFactorization]] deps = ["LinearAlgebra"] @@ -637,10 +642,10 @@ deps = ["LinearAlgebra", "Random"] uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" [[SparseDiffTools]] -deps = ["Adapt", "ArrayInterface", "Compat", "DataStructures", "DiffEqDiffTools", "ForwardDiff", "LightGraphs", "LinearAlgebra", "Requires", "SparseArrays", "VertexSafeGraphs"] -git-tree-sha1 = "15d981e62c77b117887c30c38976c3f5c8e168e5" +deps = ["Adapt", "ArrayInterface", "Compat", "DataStructures", "DiffEqDiffTools", "ForwardDiff", "InteractiveUtils", "LightGraphs", "LinearAlgebra", "Requires", "SparseArrays", "VertexSafeGraphs"] +git-tree-sha1 = "2efd0a78e58302caedf3aed0b5a9daaa7bca7439" uuid = "47a9eef4-7e08-11e9-0b38-333d64bd3804" -version = "1.0.0" +version = "1.2.0" [[SpecialFunctions]] deps = ["BinDeps", "BinaryProvider", "Libdl"] @@ -666,9 +671,9 @@ version = "0.32.0" [[StatsFuns]] deps = ["Rmath", "SpecialFunctions"] -git-tree-sha1 = "938f2520e6c093aaf5a7e7c38489e03c27a8c136" +git-tree-sha1 = "79982835d2ff3970685cb704500909c94189bde9" uuid = "4c63d2b9-4356-54db-8cca-17b64c39e42c" -version = "0.9.2" +version = "0.9.3" [[SteadyStateDiffEq]] deps = ["DiffEqBase", "DiffEqCallbacks", "LinearAlgebra", "NLsolve", "Reexport"] @@ -678,9 +683,9 @@ version = "1.5.0" [[StochasticDiffEq]] deps = ["ArrayInterface", "DataStructures", "DiffEqBase", "DiffEqDiffTools", "DiffEqNoiseProcess", "FillArrays", "ForwardDiff", "LinearAlgebra", "Logging", "MuladdMacro", "NLsolve", "Parameters", "Random", "RandomNumbers", "RecursiveArrayTools", "Reexport", "SparseArrays", "SparseDiffTools", "StaticArrays"] -git-tree-sha1 = "1fda79b9919ae2310e969c4000784858abf32f45" +git-tree-sha1 = "238db5de41192bf9695539c6eea528cd8b10f6df" uuid = "789caeaf-c7a9-5a7d-9973-96adeb23e2a0" -version = "6.15.0" +version = "6.16.1" [[SuiteSparse]] deps = ["Libdl", "LinearAlgebra", "Serialization", "SparseArrays"] @@ -744,3 +749,9 @@ deps = ["Base64", "CodecZlib", "LightXML", "Random", "TranscodingStreams"] git-tree-sha1 = "3ecc73f1c3175218b7ab37fe487dbe4cb1148e40" uuid = "64499a7a-5c06-52f2-abe2-ccb03c286192" version = "1.2.2" + +[[ZygoteRules]] +deps = ["MacroTools"] +git-tree-sha1 = "b3b4882cc9accf6731a08cc39543fbc6b669dca8" +uuid = "700de1a5-db45-46bc-99cf-38207098b444" +version = "0.2.0" diff --git a/Project.toml b/Project.toml index 545a0e74a47..1cda7b960d7 100644 --- a/Project.toml +++ b/Project.toml @@ -27,6 +27,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" Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f" Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" diff --git a/docs/Project.toml b/docs/Project.toml index 685644993da..fa466c7f9bd 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -16,6 +16,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" StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" diff --git a/docs/src/Atmos/SurfaceFluxes.md b/docs/src/Atmos/SurfaceFluxes.md index b8e440cb818..c5448a01a5a 100644 --- a/docs/src/Atmos/SurfaceFluxes.md +++ b/docs/src/Atmos/SurfaceFluxes.md @@ -4,98 +4,33 @@ CurrentModule = CLIMA.SurfaceFluxes ``` -Surface flux functions, e.g. for buoyancy flux, friction velocity, and exchange coefficients. +This module provides a means to compute surface fluxes given several variables, described in [`surface_conditions`](@ref). -## `Byun1990` - -Compute surface fluxes using the approach in Byun (1990). - -### Plots - -```@example byun1990 -using CLIMA.SurfaceFluxes.Byun1990 -using Plots, LaTeXStrings - -Ri_range = range(-1.2, stop=0.24, length=100) -scales = [50,200,600,1000,10_000] - -z_0 = 1.0 -γ_m, γ_h = 15.0, 9.0 -β_m, β_h = 4.8, 7.8 -Pr_0 = 0.74 - -plot(Ri_range, - [Byun1990.compute_exchange_coefficients(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(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 -``` -![](exchange_byun1990_fig4b.svg) - -Recreation of Figure 4(b) from Byun (1990) - -## `Nishizawa2018` - -### Plots -```@example -using CLIMA.SurfaceFluxes.Nishizawa2018 -using Plots, LaTeXStrings - -a = 4.7 -θ = 350 -z_0 = 10 -u_ave = 10 -flux = 1 -Δz = range(10.0, stop=100.0, length=100) -Ψ_m_tol, tol_abs, iter_max = 1e-3, 1e-3, 10 -u_star = Nishizawa2018.compute_friction_velocity.( - 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) 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 ## API ```@docs -compute_buoyancy_flux -Byun1990.compute_MO_len -Byun1990.compute_friction_velocity -Byun1990.compute_exchange_coefficients -Nishizawa2018.compute_MO_len -Nishizawa2018.compute_friction_velocity -Nishizawa2018.compute_exchange_coefficients +surface_conditions ``` ## 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](http://dx.doi.org/10.1175/1520-0450(1990)029<0652:OTASOF>2.0.CO;2) diff --git a/src/Atmos/Parameterizations/SurfaceFluxes/SurfaceFluxes.jl b/src/Atmos/Parameterizations/SurfaceFluxes/SurfaceFluxes.jl index 4479c683a05..8654a05cce4 100644 --- a/src/Atmos/Parameterizations/SurfaceFluxes/SurfaceFluxes.jl +++ b/src/Atmos/Parameterizations/SurfaceFluxes/SurfaceFluxes.jl @@ -1,39 +1,27 @@ """ 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: - - [`compute_MO_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 - - Ref. Byun1990: - - 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. - http://dx.doi.org/10.1175/1520-0450(1990)029<0652:OTASOF>2.0.CO;2 - - - Ref. Wyngaard1975: - - Wyngaard, John C. "Modeling the planetary boundary layer-Extension - to the stable case." Boundary-Layer Meteorology 9.4 (1975): 441-460. - https://link.springer.com/article/10.1007/BF00223393 - - Ref. Nishizawa2018: - 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. https://doi.org/10.1029/2018MS001534 + - Ref. Wyngaard1975: + - Wyngaard, John C. "Modeling the planetary boundary layer-Extension + to the stable case." Boundary-Layer Meteorology 9.4 (1975): 441-460. + https://link.springer.com/article/10.1007/BF00223393 + - Ref. Businger1971: - Businger, Joost A., et al. "Flux-profile relationships in the atmospheric surface layer." Journal of the atmospheric Sciences @@ -45,147 +33,130 @@ module SurfaceFluxes using ..RootSolvers using ..MoistThermodynamics using ..PlanetParameters +using DocStringExtensions +using NLsolve -# export compute_buoyancy_flux - -""" - compute_buoyancy_flux(shf, lhf, T_b, qt_b, ql_b, qi_b, alpha0_0) +export surface_conditions -Computes buoyancy flux given sensible heat flux `shf`, -latent heat flux `lhf`, surface boundary temperature `T_b`, -total specific humidity `qt_b`, liquid specific humidity `ql_b`, -ice specific humidity `qi_b` and specific `alpha0_0`. """ -function compute_buoyancy_flux(shf, lhf, T_b, qt_b, ql_b, qi_b, alpha0_0) - cp_ = cp_m(PhasePartition(qt_b, ql_b, qi_b)) - lv = latent_heat_vapor(T_b) - temp1 = (molmass_ratio-1) - temp2 = (shf + temp1 * cp_ * T_b * lhf /lv) - return (grav * alpha0_0 / cp_ / T_b * temp2) -end - -module Byun1990 - -using ...RootSolvers -using ...MoistThermodynamics -using ...PlanetParameters + SurfaceFluxConditions{T} -""" Computes ψ_m for stable case. See Eq. 12 Ref. Byun1990 """ -ψ_m_stable(ζ, ζ_0, β_m) = -β_m * (ζ - ζ_0) +Surface flux conditions, returned from `surface_conditions`. -""" Computes ψ_h for stable case. See Eq. 13 Ref. Byun1990 """ -ψ_h_stable(ζ, ζ_0, β_h) = -β_h * (ζ - ζ_0) +# Fields -""" Computes ψ_m for unstable case. See Eq. 14 Ref. Byun1990 """ -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 +$(DocStringExtensions.FIELDS) +""" +struct SurfaceFluxConditions{T} + L_MO::T + pottemp_flux_star::T + flux::Vector{T} + x_star::Vector{T} + K_exchange::Vector{T} end -""" Computes ψ_h for unstable case. See Eq. 15 Ref. Byun1990 """ -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 +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_MO_len(u, flux) - -Computes the Monin-Obukhov length (Eq. 3 Ref. Byun1990) -""" -compute_MO_len(u, flux) = - u^3 / (flux * k_Karman) - -""" - compute_friction_velocity(u_ave, flux, z_0, z_1, β_m, γ_m, tol_abs, iter_max) - -Computes roots of friction velocity equation (Eq. 10 in Ref. Byun1990) - -`u_ave = u_* ( ln(z/z_0) - ψ_m(z/L, z_0/L) ) /κ Eq. 10 in Ref. Byun1990` -""" -function compute_friction_velocity(u_ave, flux, z_0, z_1, β_m, γ_m, tol_abs, iter_max) - - ustar_0 = u_ave * k_Karman / 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 - - # use neutral condition as first guess - stable = z_1/compute_MO_len(ustar_0, flux) >= 0 - function compute_ψ_m(u) - L_MO = compute_MO_len(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) - return (log(z_1 / z_0) - compute_ψ_m(u)) / k_Karman # Eq. 10 in Ref. Byun1990 - 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), - ustar_0, ustar_1, SecantMethod(), CompactSolution(), - tol_abs, iter_max) - ustar = sol.root - end - +struct Momentum end +struct Heat end + +""" + surface_conditions(x_initial, x_ave, x_s, z_0, F_exchange, dimensionless_number, θ_bar, Δz, z, a) + +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. Otherwise +""" +function surface_conditions(x_initial::Vector{T}, + x_ave::Vector{T}, + x_s::Vector{T}, + z_0::Vector{T}, + F_exchange::Vector{T}, + dimensionless_number::Vector{T}, + θ_bar::T, + Δz::T, + z::T, + a::T, + pottemp_flux_given::Union{Nothing,T}=nothing + ) where T<:AbstractFloat + + 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 + 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 - compute_MO_len(u, θ_bar, pottemp_flux) + for i in 1:n_vars + ϕ = x_vec[i] + transport = i==1 ? Momentum() : Heat() + F[i+1] = ϕ - compute_physical_scale(ϕ, u, Δz, a, θ_bar, pottemp_flux, z_0[i], x_ave[i], x_s[i], dimensionless_number[i], transport) + end end - - return ustar -end - -""" - compute_exchange_coefficients(Ri, z_b, z_0, γ_m, γ_h, β_m, β_h, Pr_0) - -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) -""" -function compute_exchange_coefficients(Ri, z_b, z_0, γ_m, γ_h, β_m, β_h, Pr_0) - logz = log(z_b/z_0) - zfactor = z_b/(z_b-z_0)*logz - s_b = Ri/Pr_0 - if Ri > 0 - temp = ((1-2*β_h*Ri)-sqrt(1+4*(β_h - β_m)*s_b)) - ζ = zfactor/(2*β_h*(β_m*Ri - 1))*temp # Eq. 19 in Ref. Byun1990 - L_mo = z_b/ζ # LHS of Eq. 3 in Byun1990 - ζ_0 = z_0/L_mo - ψ_m = ψ_m_stable(ζ, ζ_0, β_m) - ψ_h = ψ_h_stable(ζ, ζ_0, β_h) + sol = nlsolve(f!, x_initial, autodiff = :forward) + if converged(sol) + L_MO, x_star = sol.zero[1], sol.zero[2:end] + u_star, θ_star = x_star[1], x_star[2] else - Q_b = 1/9 * (1 /(γ_m^2) + 3 * γ_h/γ_m * s_b^2) # Eq. 31 in Ref. Byun1990 - P_b = 1/54 * (-2/(γ_m^3) + 9/γ_m * (-γ_h/γ_m + 3)*s_b^2) # Eq. 32 in Ref. Byun1990 - crit = Q_b^3 - P_b^2 - if crit < 0 - T_b = cbrt(sqrt(-crit) + abs(P_b)) # Eq. 34 in Ref. Byun1990 - ζ = zfactor * (1/(3*γ_m)-(T_b + Q_b/T_b)) # Eq. 29 in Ref. Byun1990 - else - θ_b = acos(P_b/sqrt(Q_b^3)) # Eq. 33 in Ref. Byun1990 - ζ = zfactor * (-2 * sqrt(Q_b) * cos(θ_b/3)+1/(3*γ_m)) # Eq. 28 in Ref. Byun1990 - end - L_mo = z_b/ζ # LHS of Eq. 3 in Byun1990 - ζ_0 = z_0/L_mo - ψ_m = ψ_m_unstable(ζ, ζ_0, γ_m) - ψ_h = ψ_h_unstable(ζ, ζ_0, γ_h) + @show sol + error("Unconverged surface fluxes") end - cu = k_Karman/(logz-ψ_m) # Eq. 10 in Ref. Byun1990, solved for u^* - cth = k_Karman/(logz-ψ_h)/Pr_0 # Eq. 11 in Ref. Byun1990, solved for h^* - C_D = cu^2 # Eq. 36 in Byun1990 - C_H = cu*cth # Eq. 37 in Byun1990 - return C_D, C_H, L_mo + + pottemp_flux_star = -u_star^3*θ_bar/(k_Karman*grav*L_MO) + flux = -u_star*x_star + K_exchange = compute_exchange_coefficients(z, F_exchange, a, x_star, θ_bar, dimensionless_number, L_MO) + + return SurfaceFluxConditions(L_MO, + pottemp_flux_star, + flux, + x_star, + K_exchange) end -end # Byun1990 module -module Nishizawa2018 -using ...RootSolvers -using ...MoistThermodynamics -using ...PlanetParameters +function compute_physical_scale(x, u, Δz, a, θ_bar, pottemp_flux, z_0, x_ave, x_s, dimensionless_number, transport) + L = compute_MO_len(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)*k_Karman/(temp1+temp2+temp3+temp4)*(x_ave-x_s) +end """ Computes R_z0 expression, defined after Eq. 15 Ref. Nishizawa2018 """ compute_R_z0(z_0, Δz) = 1 - z_0/Δz @@ -196,19 +167,19 @@ compute_f_m(ζ) = sqrt(sqrt(1-15*ζ)) """ Computes f_h in Eq. A8 Ref. Nishizawa2018 """ compute_f_h(ζ) = sqrt(1-9*ζ) -""" Computes ψ_m in Eq. A3 Ref. Nishizawa2018 """ -function compute_ψ_m(ζ, L, a) +""" Computes psi_m in Eq. A3 Ref. Nishizawa2018 """ +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)+π/2 end -""" Computes ψ_h in Eq. A4 Ref. Nishizawa2018 """ -compute_ψ_h(ζ, L, a, Pr) = L>=0 ? -a*ζ/Pr : 2*log((1+compute_f_h(ζ))/2) +""" Computes psi_h in Eq. A4 Ref. Nishizawa2018 """ +compute_psi(ζ, L, a, dimensionless_number, ::Heat) = L>=0 ? -a*ζ/dimensionless_number : 2*log((1+compute_f_h(ζ))/2) -""" Computes Ψ_m in Eq. A5 Ref. Nishizawa2018 """ -function compute_Ψ_m(ζ, L, a, tol) - if ζ < tol - return ζ>=0 ? -a*ζ/2 : -15*ζ/8 # Computes Ψ_m in Eq. A13 Ref. Nishizawa2018 +""" Computes Psi_m in Eq. A5 Ref. Nishizawa2018 """ +function compute_Psi(ζ, L, a, dimensionless_number, ::Momentum) + if abs(ζ) < eps(typeof(L)) + return ζ>=0 ? -a*ζ/2 : -15*ζ/8 # Computes Psi_m in Eq. A13 Ref. Nishizawa2018 else f_m = compute_f_m(ζ) # Note that "1-f^3" in Ref. Nishizawa2018 is a typo, it is supposed to be "1-f_m^3". @@ -217,13 +188,13 @@ function compute_Ψ_m(ζ, L, a, tol) end end -""" Computes Ψ_h in Eq. A6 Ref. Nishizawa2018 """ -function compute_Ψ_h(ζ, L, a, Pr, tol) - if ζ < tol - return ζ>=0 ? -a*ζ/(2*Pr) : -9*ζ/4 # Computes Ψ_h in Eq. A14 Ref. Nishizawa2018 +""" Computes Psi_h in Eq. A6 Ref. Nishizawa2018 """ +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 Ref. Nishizawa2018 else f_h = compute_f_h(ζ) - return L>=0 ? -a*ζ/(2*Pr) : 2*log((1+f_h)/2) + 2*(1-f_h)/(9*ζ) + return L>=0 ? -a*ζ/(2*dimensionless_number) : 2*log((1+f_h)/2) + 2*(1-f_h)/(9*ζ) end end @@ -232,67 +203,26 @@ end Computes Monin-Obukhov length. Eq. 3 Ref. Nishizawa2018 """ -compute_MO_len(u, θ, flux) = - u^3* θ / (k_Karman * grav * flux) - -""" - compute_friction_velocity(u_ave, θ, flux, Δz, z_0, a, Ψ_m_tol, tol_abs, iter_max) - -Computes friction velocity, in Eq. 12 in -Ref. Nishizawa2018, 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 `compute_MO_len`). -""" -function compute_friction_velocity(u_ave, θ, flux, Δz, z_0, a, Ψ_m_tol, tol_abs, iter_max) - ustar_0 = u_ave * k_Karman / log(Δz / z_0) - ustar = ustar_0 - let u_ave=u_ave, θ=θ, 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 = compute_MO_len(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) / k_Karman - 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), - ustar_0, ustar_1, SecantMethod(), CompactSolution(), - tol_abs, iter_max) - ustar = sol.root - end - return ustar -end +compute_MO_len(u, θ_bar, flux) = - u^3* θ_bar / (k_Karman * grav * flux) """ - compute_exchange_coefficients(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 + - `K_exchange` exchange coefficients """ -function compute_exchange_coefficients(z, F_m, F_h, a, u_star, θ, flux, Pr) +function compute_exchange_coefficients(z, F_exchange, a, x_star, θ_bar, dimensionless_number, L_MO) - L_mo = compute_MO_len(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*k_Karman*z/(u_star * ψ_m) # Eq. 19 in Ref. Nishizawa2018 - K_h = -F_h*k_Karman*z/(Pr * θ * ψ_h) # Eq. 20 in Ref. Nishizawa2018 + N = length(F_exchange) + K_exchange = Vector{typeof(z)}(undef, N) + 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]*k_Karman*z/(x_star[i] * psi) # Eq. 19 in Ref. Nishizawa2018 + end - return K_m, K_h, L_mo + return K_exchange end -end # Nishizawa2018 module - end # SurfaceFluxes module diff --git a/test/Atmos/Parameterizations/SurfaceFluxes/runtests.jl b/test/Atmos/Parameterizations/SurfaceFluxes/runtests.jl index 915fef5942c..e4e6477b14c 100644 --- a/test/Atmos/Parameterizations/SurfaceFluxes/runtests.jl +++ b/test/Atmos/Parameterizations/SurfaceFluxes/runtests.jl @@ -1,8 +1,6 @@ using Test using CLIMA.SurfaceFluxes -using CLIMA.SurfaceFluxes.Nishizawa2018 -using CLIMA.SurfaceFluxes.Byun1990 using CLIMA.MoistThermodynamics using CLIMA.RootSolvers @@ -12,69 +10,34 @@ using CLIMA.RootSolvers @testset "SurfaceFluxes" begin - shf, lhf, T_b, qt_b, ql_b, qi_b, alpha0_0 = 60.0, 50.0, 350.0, 0.01, 0.002, 0.0001, 1.0 - buoyancy_flux = SurfaceFluxes.compute_buoyancy_flux(shf, - lhf, - T_b, - qt_b, - ql_b, - qi_b, - alpha0_0 - ) - @test buoyancy_flux ≈ 0.0017808608107074118 -end - -@testset "SurfaceFluxes.Byun1990" begin - u, flux = 0.1, 0.2 - MO_len = Byun1990.compute_MO_len(u, flux) - @test MO_len ≈ -0.0125 - - u_ave, buoyancy_flux, z_0, z_1 = 0.1, 0.2, 2.0, 5.0 - γ_m, β_m = 15.0, 4.8 - tol_abs, iter_max = 1e-3, 10 - u_star = Byun1990.compute_friction_velocity(u_ave, - buoyancy_flux, - z_0, - z_1, - β_m, - γ_m, - tol_abs, - iter_max - ) - @test u_star ≈ 0.201347256193615 atol=tol_abs - - - Ri, z_b, z_0, Pr_0 = 10, 2.0, 5.0, 0.74 - γ_m, γ_h, β_m, β_h = 15.0, 9.0, 4.8, 7.8 - cm, ch, L_mo = Byun1990.compute_exchange_coefficients(Ri, z_b, z_0, γ_m, γ_h, β_m, β_h, Pr_0) - @test cm ≈ 19.700348427787368 - @test ch ≈ 3.3362564728997803 - @test L_mo ≈ -14.308268023583906 - - Ri, z_b, z_0, Pr_0 = -10, 10.0, 1.0, 0.74 - γ_m, γ_h, β_m, β_h = 15.0, 9.0, 4.8, 7.8 - cm, ch, L_mo = Byun1990.compute_exchange_coefficients(Ri, z_b, z_0, γ_m, γ_h, β_m, β_h, Pr_0) - @test cm ≈ 0.33300280321092746 - @test ch ≈ 1.131830939627489 - @test L_mo ≈ -0.3726237964444814 + x_initial = [100, 15.0, 350.0] + F_exchange = [2.0, 3.0] + z_0 = [1.0, 1.0] + dimensionless_number = [1.0, 0.74] + x_ave = [5.0, 350.0] + x_s = [0.0, 300.0] + + Δz = 2.0 + z = 0.5 + θ_bar = 300.0 + a = 4.7 + pottemp_flux_given = -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(args[1:end-1]...) + + @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(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 -@testset "SurfaceFluxes.Nishizawa2018" begin - u, θ, flux = 2, 350, 20 - MO_len = Nishizawa2018.compute_MO_len(u, θ, flux) - @test MO_len ≈ -35.67787971457696 - - u_ave, θ, flux, Δz, z_0, a = 110, 350, 20, 100, 0.01, 5 - Ψ_m_tol, tol_abs, iter_max = 1e-3, 1e-3, 10 - u_star = Nishizawa2018.compute_friction_velocity(u_ave, θ, flux, Δz, z_0, a, Ψ_m_tol, tol_abs, iter_max) - @test u_star ≈ 5.526644550864822 atol=tol_abs - - z, F_m, F_h, a, u_star, θ, flux, Pr = 1.0, 2.0, 3.0, 5, 110, 350, 20, 0.74 - - K_m, K_h, L_mo = Nishizawa2018.compute_exchange_coefficients(z,F_m,F_h,a,u_star,θ,flux,Pr) - @test K_m ≈ -11512.071612359368 - @test K_h ≈ -6111.6196776263805 - @test L_mo ≈ -5.935907237512742e6 - -end