From 4f7255e558a2141a3aa3767c3b0f0809b79630e1 Mon Sep 17 00:00:00 2001 From: Charles Kawczynski Date: Mon, 4 May 2020 08:29:23 -0700 Subject: [PATCH] Synchronize with clima, add docs and formatter --- .dev/Project.toml | 5 + .dev/clima_formatter_options.jl | 8 + .dev/format.jl | 18 + .github/workflows/Docs.yml | 26 ++ .gitignore | 18 +- Manifest.toml | 82 +---- Project.toml | 6 - README.md | 10 +- azure-pipelines.yml | 40 --- bors.toml | 3 +- docs/Project.toml | 2 + docs/make.jl | 46 +++ docs/src/API.md | 53 +++ docs/src/HowToGuide.md | 89 +++++ docs/src/Installation.md | 7 + docs/src/TestedProfiles.md | 85 +++++ docs/src/assets/logo.svg | 13 + docs/src/index.md | 37 ++ env/Plots/Project.toml | 2 + env/test/Project.toml | 5 + src/MoistThermodynamics.jl | 22 +- src/isentropic.jl | 34 +- src/profiles.jl | 81 +++++ src/relations.jl | 610 ++++++++++++++++++-------------- src/states.jl | 160 +++------ test/data_tests.jl | 20 +- test/runtests.jl | 398 ++++++++++++++------- 27 files changed, 1227 insertions(+), 653 deletions(-) create mode 100644 .dev/Project.toml create mode 100644 .dev/clima_formatter_options.jl create mode 100644 .dev/format.jl create mode 100644 .github/workflows/Docs.yml create mode 100644 docs/Project.toml create mode 100644 docs/make.jl create mode 100644 docs/src/API.md create mode 100644 docs/src/HowToGuide.md create mode 100644 docs/src/Installation.md create mode 100644 docs/src/TestedProfiles.md create mode 100644 docs/src/assets/logo.svg create mode 100644 docs/src/index.md create mode 100644 env/Plots/Project.toml create mode 100644 env/test/Project.toml create mode 100644 src/profiles.jl diff --git a/.dev/Project.toml b/.dev/Project.toml new file mode 100644 index 00000000..d86279dc --- /dev/null +++ b/.dev/Project.toml @@ -0,0 +1,5 @@ +[deps] +JuliaFormatter = "98e50ef6-434e-11e9-1051-2b60c6c9e899" + +[compat] +JuliaFormatter = "0.3" diff --git a/.dev/clima_formatter_options.jl b/.dev/clima_formatter_options.jl new file mode 100644 index 00000000..1db65e8c --- /dev/null +++ b/.dev/clima_formatter_options.jl @@ -0,0 +1,8 @@ +clima_formatter_options = ( + indent = 4, + margin = 80, + always_for_in = true, + whitespace_typedefs = true, + whitespace_ops_in_indices = true, + remove_extra_newlines = false, +) diff --git a/.dev/format.jl b/.dev/format.jl new file mode 100644 index 00000000..a1da12a6 --- /dev/null +++ b/.dev/format.jl @@ -0,0 +1,18 @@ +#!/usr/bin/env julia + +using Pkg +Pkg.activate(@__DIR__) +Pkg.instantiate() + +using JuliaFormatter + +include("clima_formatter_options.jl") + +headbranch = get(ARGS, 1, "master") + +for filename in + readlines(`git diff --name-only --diff-filter=AM $headbranch...`) + endswith(filename, ".jl") || continue + + format(filename; clima_formatter_options...) +end diff --git a/.github/workflows/Docs.yml b/.github/workflows/Docs.yml new file mode 100644 index 00000000..2a47c639 --- /dev/null +++ b/.github/workflows/Docs.yml @@ -0,0 +1,26 @@ +name: Documentation + +on: + push: + branches: + - master + - trying + - staging + tags: '*' + pull_request: + +jobs: + docs-build: + runs-on: ubuntu-latest + steps: + - uses: actions/checkout@v2 + - uses: julia-actions/setup-julia@latest + with: + version: 1.3 + - name: Install dependencies + run: julia --project=docs/ -e 'using Pkg; Pkg.develop(PackageSpec(path=pwd())); Pkg.instantiate()' + - name: Build and deploy + env: + GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }} # For authentication with GitHub Actions token + DOCUMENTER_KEY: ${{ secrets.DOCUMENTER_KEY }} # For authentication with SSH deploy key + run: julia --project=docs/ docs/make.jl \ No newline at end of file diff --git a/.gitignore b/.gitignore index 421cbf2a..d37fb255 100644 --- a/.gitignore +++ b/.gitignore @@ -1,4 +1,11 @@ +# Temporary *.jl.cov +*.DS_Store +*.tar.gz +*.jl.*.cov +*.jl.mem + +# Data *.vtk *.dat *.csv @@ -6,11 +13,12 @@ *.pvtu *.swp *.png +*.nc + +# Docs !docs/src/assets/*.png -*.DS_Store docs/build/ -data/ docs/site/ -*.jl.*.cov -*.jl.mem -deps/deps.jl + +# Deps +Manifest.toml diff --git a/Manifest.toml b/Manifest.toml index 92194a31..bc3b2537 100644 --- a/Manifest.toml +++ b/Manifest.toml @@ -3,23 +3,11 @@ [[Base64]] uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f" -[[BinDeps]] -deps = ["Libdl", "Pkg", "SHA", "URIParser", "Unicode"] -git-tree-sha1 = "66158ad56b4bf6cc8413b37d0b7bc52402682764" -uuid = "9e28174c-4ba2-5203-b857-d8d62c4213ee" -version = "1.0.0" - -[[CFTime]] -deps = ["Dates", "Printf"] -git-tree-sha1 = "143ef231a14c2ab1ac7e838c559ee1e45c0d1b57" -uuid = "179af706-886a-5703-950a-314cd64e0468" -version = "0.1.0" - [[CLIMAParameters]] deps = ["Test"] -git-tree-sha1 = "e07ffe98ec361798673445962f211659d14205ea" +git-tree-sha1 = "7b80821bffd8702ccc2a712aef804e30b3bd7012" uuid = "6eacf6c3-8458-43b9-ae03-caf5306d3d53" -version = "0.1.1" +version = "0.1.3" [[CommonSubexpressions]] deps = ["Test"] @@ -33,24 +21,6 @@ git-tree-sha1 = "7c4f882c41faa72118841185afc58a2eb00ef612" uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae" version = "0.3.3+0" -[[Conda]] -deps = ["JSON", "VersionParsing"] -git-tree-sha1 = "7a58bb32ce5d85f8bf7559aa7c2842f9aecf52fc" -uuid = "8f4d0f93-b110-5947-807f-2305c1781a2d" -version = "1.4.1" - -[[CondaBinDeps]] -deps = ["BinDeps", "Conda"] -git-tree-sha1 = "25f750df2893991f2c9b18425bfac6f2ce855154" -uuid = "a9693cdc-2bc8-5703-a9cd-1da358117377" -version = "0.2.0" - -[[DataStructures]] -deps = ["InteractiveUtils", "OrderedCollections"] -git-tree-sha1 = "5a431d46abf2ef2a4d5d00bd0ae61f651cf854c8" -uuid = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8" -version = "0.17.10" - [[Dates]] deps = ["Printf"] uuid = "ade2ca70-3891-5945-98fb-dc099432e06a" @@ -87,12 +57,6 @@ version = "0.10.10" deps = ["Markdown"] uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240" -[[JSON]] -deps = ["Dates", "Mmap", "Parsers", "Unicode"] -git-tree-sha1 = "b34d7cef7b337321e97d22242c3c2b91f476748e" -uuid = "682c06a0-de6a-54ab-a142-c8b1cf79cde6" -version = "0.21.0" - [[LibGit2]] uuid = "76f85450-5226-5b5a-8eaa-529ad045b433" @@ -110,15 +74,6 @@ uuid = "56ddb016-857b-54e1-b83d-db4d58db5568" deps = ["Base64"] uuid = "d6f4376e-aef5-505a-96c1-9c027394607a" -[[Mmap]] -uuid = "a63ad114-7e13-5084-954f-fe012c677804" - -[[NCDatasets]] -deps = ["BinDeps", "CFTime", "CondaBinDeps", "DataStructures", "Dates", "Libdl", "Printf"] -git-tree-sha1 = "89d5c2259694a459194983105be10b43f6fbefc8" -uuid = "85f8d34a-cbdd-5861-8df4-14fed0d494ab" -version = "0.10.1" - [[NaNMath]] git-tree-sha1 = "928b8ca9b2791081dc71a51c55347c27c618760f" uuid = "77ba4419-2d1f-58cd-9bb1-8ffee604a2e3" @@ -130,18 +85,6 @@ git-tree-sha1 = "d51c416559217d974a1113522d5919235ae67a87" uuid = "efe28fd5-8261-553b-a9e1-b2916fc3738e" version = "0.5.3+3" -[[OrderedCollections]] -deps = ["Random", "Serialization", "Test"] -git-tree-sha1 = "c4c13474d23c60d20a67b217f1d7f22a40edf8f1" -uuid = "bac558e1-5e72-5ebc-8fee-abe8a469f55d" -version = "1.1.0" - -[[Parsers]] -deps = ["Dates", "Test"] -git-tree-sha1 = "d6d82d5bdbb75048e574cd2d2c89dfbf2c74250c" -uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0" -version = "1.0.0" - [[Pkg]] deps = ["Dates", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "Test", "UUIDs"] uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f" @@ -159,10 +102,10 @@ deps = ["Serialization"] uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" [[RootSolvers]] -deps = ["ForwardDiff", "Test"] -git-tree-sha1 = "5e047b4bea5cb80adc6a965082595b5a16bd318d" +deps = ["DocStringExtensions", "ForwardDiff", "Test"] +git-tree-sha1 = "0e5b394adc5c6fb39b3964bce2a259a44cc312d3" uuid = "7181ea78-2dcb-4de3-ab41-2b8ab5a31e74" -version = "0.1.0" +version = "0.2.0" [[SHA]] uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce" @@ -185,9 +128,9 @@ version = "0.10.0" [[StaticArrays]] deps = ["LinearAlgebra", "Random", "Statistics"] -git-tree-sha1 = "5a3bcb6233adabde68ebc97be66e95dcb787424c" +git-tree-sha1 = "5c06c0aeb81bef54aed4b3f446847905eb6cbda0" uuid = "90137ffa-7385-5640-81b9-e52037218182" -version = "0.12.1" +version = "0.12.3" [[Statistics]] deps = ["LinearAlgebra", "SparseArrays"] @@ -197,20 +140,9 @@ uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" deps = ["Distributed", "InteractiveUtils", "Logging", "Random"] uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40" -[[URIParser]] -deps = ["Test", "Unicode"] -git-tree-sha1 = "6ddf8244220dfda2f17539fa8c9de20d6c575b69" -uuid = "30578b45-9adc-5946-b283-645ec420af67" -version = "0.4.0" - [[UUIDs]] deps = ["Random", "SHA"] uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4" [[Unicode]] uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5" - -[[VersionParsing]] -git-tree-sha1 = "80229be1f670524750d905f8fc8148e5a8c4537f" -uuid = "81def892-9a0e-5fdd-b105-ffc91e053289" -version = "1.2.0" diff --git a/Project.toml b/Project.toml index 5aa9ffdd..fe137190 100644 --- a/Project.toml +++ b/Project.toml @@ -6,15 +6,9 @@ version = "0.1.0" [deps] CLIMAParameters = "6eacf6c3-8458-43b9-ae03-caf5306d3d53" DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae" -NCDatasets = "85f8d34a-cbdd-5861-8df4-14fed0d494ab" -Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f" -Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" RootSolvers = "7181ea78-2dcb-4de3-ab41-2b8ab5a31e74" -Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [compat] CLIMAParameters = "0.1.0" DocStringExtensions = "0.8.1" -NCDatasets = "0.10.1" -RootSolvers = "0.1.0" julia = "1.3" diff --git a/README.md b/README.md index 10475c63..a29e9f40 100644 --- a/README.md +++ b/README.md @@ -3,13 +3,17 @@ A package containing a library of moist thermodynamic relations ||| |---------------------:|:----------------------------------------------| -| **Documentation** | [![latest][docs-latest-img]][docs-latest-url] | +| **Docs Build** | [![docs build][docs-bld-img]][docs-bld-url] | +| **Documentation** | [![dev][docs-dev-img]][docs-dev-url] | | **Azure Build** | [![azure][azure-img]][azure-url] | | **Code Coverage** | [![codecov][codecov-img]][codecov-url] | | **Bors** | [![Bors enabled][bors-img]][bors-url] | -[docs-latest-img]: https://img.shields.io/badge/docs-latest-blue.svg -[docs-latest-url]: https://climate-machine.github.io/MoistThermodynamics.jl/latest/ +[docs-dev-img]: https://img.shields.io/badge/docs-dev-blue.svg +[docs-dev-url]: https://climate-machine.github.io/MoistThermodynamics.jl/dev/ + +[docs-dev-img]: https://img.shields.io/badge/docs-dev-blue.svg +[docs-dev-url]: https://climate-machine.github.io/MoistThermodynamics.jl/dev/ [azure-img]: https://dev.azure.com/climate-machine/MoistThermodynamics.jl/_apis/build/status/climate-machine.MoistThermodynamics.jl?branchName=master [azure-url]: https://dev.azure.com/climate-machine/MoistThermodynamics.jl/_build/latest?definitionId=1&branchName=master diff --git a/azure-pipelines.yml b/azure-pipelines.yml index bf07dd6d..9f0a2a3e 100644 --- a/azure-pipelines.yml +++ b/azure-pipelines.yml @@ -111,43 +111,3 @@ jobs: C:\julia-$(JULIA_VERSION)\bin\julia.exe --project=@. -e 'using Pkg; Pkg.instantiate()' C:\julia-$(JULIA_VERSION)\bin\julia.exe --project=@. -e 'using Pkg; Pkg.test()' displayName: 'Run the tests' - -# - job: Documentation - -# pool: -# vmImage: 'ubuntu-16.04' - -# strategy: -# matrix: -# Julia 1.3: -# JULIA_VERSION: '1.3' - -# steps: -# - bash: | -# set -o xtrace -# wget -nv https://julialang-s3.julialang.org/bin/linux/x64/$(JULIA_VERSION)/julia-$(JULIA_VERSION)-latest-linux-x86_64.tar.gz -# mkdir julia-$(JULIA_VERSION) -# tar zxf julia-$(JULIA_VERSION)-latest-linux-x86_64.tar.gz -C julia-$(JULIA_VERSION) --strip-components 1 -# displayName: 'Download and extract Julia' -# - bash: | -# set -o xtrace -# sudo apt-get update -# displayName: 'Install dependencies' -# - bash: | -# set -o xtrace -# export TRAVIS_REPO_SLUG="$BUILD_REPOSITORY_NAME" -# export TRAVIS_PULL_REQUEST="${SYSTEM_PULLREQUEST_PULLREQUESTNUMBER:-false}" -# if [[ $BUILD_SOURCEBRANCH == refs/tags/* ]]; then -# export TRAVIS_TAG="${BUILD_SOURCEBRANCH:10}" -# fi -# if [[ $BUILD_SOURCEBRANCH == refs/heads/* ]]; then -# export TRAVIS_BRANCH="${BUILD_SOURCEBRANCH:11}" -# fi -# ./julia-$(JULIA_VERSION)/bin/julia -e 'using InteractiveUtils; versioninfo()' -# ./julia-$(JULIA_VERSION)/bin/julia --project=docs/ -e 'using Pkg; Pkg.instantiate(); -# Pkg.develop(PackageSpec(path=pwd())); -# Pkg.build()' -# ./julia-$(JULIA_VERSION)/bin/julia --project=docs/ docs/make.jl -# env: -# DOCUMENTER_KEY: $(documenter_key) -# displayName: 'Build and deploy docs' diff --git a/bors.toml b/bors.toml index 7919a0dc..aaced1c9 100644 --- a/bors.toml +++ b/bors.toml @@ -1,5 +1,6 @@ status = [ - "climate-machine.MoistThermodynamics.jl" + "climate-machine.MoistThermodynamics.jl", + "docs-build" ] delete_merged_branches = true timeout_sec = 86400 diff --git a/docs/Project.toml b/docs/Project.toml new file mode 100644 index 00000000..dfa65cd1 --- /dev/null +++ b/docs/Project.toml @@ -0,0 +1,2 @@ +[deps] +Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" diff --git a/docs/make.jl b/docs/make.jl new file mode 100644 index 00000000..3d2eef5e --- /dev/null +++ b/docs/make.jl @@ -0,0 +1,46 @@ + +rm(joinpath(@__DIR__, "Manifest.toml"), force = true) # Remove local Manifest + +push!(LOAD_PATH, joinpath(@__DIR__, "..", "env", "Plots")) # add Plots env +push!(LOAD_PATH, joinpath(@__DIR__, "..")) # add MoistThermodynamics env + +using Pkg +Pkg.activate(@__DIR__) +Pkg.instantiate(; verbose = true) + +using MoistThermodynamics, Documenter + +pages = Any[ + "Home" => "index.md", + "Installation" => "Installation.md", + "API" => "API.md", + "How-to-guide" => "HowToGuide.md", + "Tested profiles" => "TestedProfiles.md", +] + +mathengine = MathJax(Dict( + :TeX => Dict( + :equationNumbers => Dict(:autoNumber => "AMS"), + :Macros => Dict(), + ), +)) + +format = Documenter.HTML( + prettyurls = get(ENV, "CI", nothing) == "true", + mathengine = mathengine, + collapselevel = 1, +) + +makedocs( + sitename = "MoistThermodynamics.jl", + format = format, + clean = true, + modules = [Documenter, MoistThermodynamics], + pages = pages, +) + +deploydocs( + repo = "github.com/climate-machine/MoistThermodynamics.jl.git", + target = "build", + push_preview = true, +) diff --git a/docs/src/API.md b/docs/src/API.md new file mode 100644 index 00000000..80321090 --- /dev/null +++ b/docs/src/API.md @@ -0,0 +1,53 @@ +# API + +```@meta +CurrentModule = MoistThermodynamics +``` + +## Thermodynamic State Constructors + +```@docs +PhasePartition +PhasePartition_equil +ThermodynamicState +PhaseDry +PhaseEquil +PhaseNonEquil +TemperatureSHumEquil +LiquidIcePotTempSHumEquil +LiquidIcePotTempSHumNonEquil +LiquidIcePotTempSHumNonEquil_given_pressure +``` + +## Thermodynamic state methods +```@docs +air_density +air_pressure +air_temperature +air_temperature_from_liquid_ice_pottemp +cp_m +cv_m +dry_pottemp +exner +gas_constant_air +Ice +internal_energy +internal_energy_sat +latent_heat_fusion +latent_heat_sublim +latent_heat_vapor +Liquid +liquid_fraction +liquid_ice_pottemp +liquid_ice_pottemp_sat +gas_constants +saturation_adjustment +saturation_excess +q_vap_saturation +q_vap_saturation_generic +saturation_vapor_pressure +soundspeed_air +specific_volume +total_energy +virtual_pottemp +``` diff --git a/docs/src/HowToGuide.md b/docs/src/HowToGuide.md new file mode 100644 index 00000000..ac17cce7 --- /dev/null +++ b/docs/src/HowToGuide.md @@ -0,0 +1,89 @@ +# How-to-guide + +## Usage + +Users are encouraged to first establish a thermodynamic state with one of our [Thermodynamic State Constructors](@ref). For example, we would construct a moist thermodynamic state using + +```julia +ts = PhaseEquil(param_set, e_int, ρ, q_tot); +``` + +here, `ρ` is the density of the moist air, and the internal energy `e_int = e_tot - e_kin - geopotential` is the total energy `e_tot` minus kinetic energy `e_kin` and potential energy `geopotential` (all energies per unit mass). Once we've established a thermodynamic state, we can call [Thermodynamic state methods](@ref) that support thermodynamic states: + +```julia +T = air_temperature(ts); +q = PhasePartition(ts); +``` + +No changes to the "right-hand sides" of the dynamical equations are needed for a moist dynamical core that supports clouds, as long as they do not precipitate. Additional source-sink terms arise from precipitation. + +Schematically, the workflow in such a core would look as follows: +```julia +# initialize +geopotential = grav * z +q_tot = ... +ρ = ... + +(u, v, w) = ... +e_kin = 0.5 * (u^2 + v^2 + w^2) + +e_tot = total_energy(e_kin, geopotential, T, q_tot) + +do timestep # timestepping loop + + # advance dynamical variables by a timestep (temperature typically + # appears in terms on the rhs, such as radiative transfer) + advance(u, v, w, ρ, e_tot, q_tot) + + # compute internal energy from dynamic variables + e_int = e_tot - 0.5 * (u^2 + v^2 + w^2) - geopotential + + # compute temperature, pressure and condensate specific humidities, + ts = PhaseEquil(param_set, e_int, ρ, q_tot); + T = air_temperature(ts); + q = PhasePartition(ts); + p = air_pressure(ts); + +end +``` + +For a dynamical core that additionally uses the liquid and ice specific humidities `q.liq` and `q.ice` as prognostic variables, and thus explicitly allows the presence of non-equilibrium phases such as supercooled water, the saturation adjustment in the above workflow is replaced calling a non-equilibrium moist thermodynamic state: +```julia +q_tot, q_liq, q_ice = ... +ts = PhaseNonEquil(param_set, e_int, ρ, PhasePartition(q_tot, q_liq, q_ice)); +T = air_temperature(ts); +p = air_pressure(ts); +``` + +## Extending + +If MoistThermodynamics.jl does not have a particular thermodynamic constructor that is needed, one can implement a new one in `MoistThermodynamics/src/states.jl`. In this constructor, one must add whichever arguments they wish to offer as inputs, then translate this thermodynamic state into one of: + + - `PhaseDry` a dry thermodynamic state, uniquely determined by two independent thermodynamic properties + - `PhaseEquil` a moist thermodynamic state in thermodynamic equilibrium, uniquely determined by three independent thermodynamic properties + - `PhaseNonEquil` a moist thermodynamic state in thermodynamic non-equilibrium, uniquely determined by four independent thermodynamic properties + +For example, to add a thermodynamic state constructor that accepts temperature, density and total specific humidity, we could add the following code to states: + +``` +""" + TemperatureSHumEquil_given_density(param_set, T, ρ, q_tot) + +Constructs a [`PhaseEquil`](@ref) thermodynamic state from temperature. + + - `param_set` parameter set, used to dispatch planet parameter function calls + - `T` temperature + - `ρ` density + - `q_tot` total specific humidity +""" +function TemperatureSHumEquil( + param_set::APS, + T::FT, + ρ::FT, + q_tot::FT, +) where {FT <: Real} + q = PhasePartition_equil(param_set, T, ρ, q_tot) + e_int = internal_energy(param_set, T, q) + return PhaseEquil{FT, typeof(param_set)}(param_set, e_int, ρ, q_tot, T) +end +``` diff --git a/docs/src/Installation.md b/docs/src/Installation.md new file mode 100644 index 00000000..0fe5280c --- /dev/null +++ b/docs/src/Installation.md @@ -0,0 +1,7 @@ +# Installation + +MoistThermodynamics.jl is a Julia registered package, and can be added from the Julia Pkg manager: + +```julia +(v1.x) pkg> add MoistThermodynamics +``` diff --git a/docs/src/TestedProfiles.md b/docs/src/TestedProfiles.md new file mode 100644 index 00000000..d95654c9 --- /dev/null +++ b/docs/src/TestedProfiles.md @@ -0,0 +1,85 @@ +# Tested Profiles + +MoistThermodynamics.jl is tested using a set of profiles, or thermodynamic state regimes, specified in [`tested_profiles`](@ref MoistThermodynamics.tested_profiles). + +## Dry Phase + +```@example +using MoistThermodynamics +MT = MoistThermodynamics +using CLIMAParameters +using CLIMAParameters.Planet +using Plots + +struct EarthParameterSet <: AbstractEarthParameterSet end; +const param_set = EarthParameterSet(); + +FT = Float64; +e_int, ρ, q_tot, q_pt, T, p, θ_liq_ice = MT.tested_profiles(param_set, 50, FT); + +mask_dry = q_tot .≈ 0; +ρ_dry = ρ[mask_dry]; +T_dry = T[mask_dry]; + +scatter(ρ_dry, T_dry, xlabel="density [kg/m^3]", ylabel="T [K]", title="Tested states for dry thermodynamic phase", legend=false); +savefig("tested_profiles_dry.svg"); +``` +![](tested_profiles_dry.svg) + +## Moist Phase in thermodynamic equilibrium + +Here is a 2D representation: +```@example +using MoistThermodynamics +MT = MoistThermodynamics +using CLIMAParameters +using CLIMAParameters.Planet +using Plots + +struct EarthParameterSet <: AbstractEarthParameterSet end; +const param_set = EarthParameterSet(); + +FT = Float64; +e_int, ρ, q_tot, q_pt, T, p, θ_liq_ice = MT.tested_profiles(param_set, 50, FT); +scatter(ρ, T, xlabel="density [kg/m^3]", ylabel="T [K]", marker_z=q_tot, title="Tested states for moist thermodynamic phase", label="q_tot"); +savefig("tested_profiles.svg") +``` +![](tested_profiles.svg) + +And here is a 3D representation: +```@example +using MoistThermodynamics +MT = MoistThermodynamics +using CLIMAParameters +using CLIMAParameters.Planet +using Plots + +struct EarthParameterSet <: AbstractEarthParameterSet end; +const param_set = EarthParameterSet(); + +FT = Float64; +e_int, ρ, q_tot, q_pt, T, p, θ_liq_ice = MT.tested_profiles(param_set, 50, FT); + +# initialize a 3D plot with 1 empty series +plt = plot3d( + 1, + xlim = (min(ρ...), max(ρ...)), + ylim = (min(T...), max(T...)), + zlim = (min(q_tot...), max(q_tot...)), + xlabel="density [kg/m^3]", + ylabel="T [K]", + zlabel="total specific humidity []", + legend=false, + title = "Tested states for moist thermodynamic phase", + marker = 2, +) + +# build an animated gif by pushing new points to the plot, saving every nth frame +@gif for i=1:length(ρ) + push!(plt, ρ[i], T[i], q_tot[i]) +end every 5 +``` + +## Moist Phase in thermodynamic non-equilibrium + +In progress... diff --git a/docs/src/assets/logo.svg b/docs/src/assets/logo.svg new file mode 100644 index 00000000..ad72caa6 --- /dev/null +++ b/docs/src/assets/logo.svg @@ -0,0 +1,13 @@ + + + + + + + + + + + + + \ No newline at end of file diff --git a/docs/src/index.md b/docs/src/index.md new file mode 100644 index 00000000..a8fa0183 --- /dev/null +++ b/docs/src/index.md @@ -0,0 +1,37 @@ +# MoistThermodynamics + +MoistThermodynamics.jl provides all thermodynamic functions needed for the atmosphere and functions shared across model components. The functions are general for a moist atmosphere that includes suspended cloud condensate in the working fluid; the special case of a dry atmosphere is obtained for zero specific humidities (or simply by omitting the optional specific humidity arguments in the functions that are needed for a dry atmosphere). The general formulation assumes that there are tracers for specific humidity `q`, partitioned into + + - `q.tot` total water specific humidity + - `q.liq` liquid specific humidity + - `q.ice` ice specific humidity + +to characterize the thermodynamic state and composition of moist air. + +There are several types of functions: + +1. Equation of state (ideal gas law): + * `air_pressure` +2. Specific gas constant and isobaric and isochoric specific heats of moist air: + * `gas_constant_air` + * `cp_m` + * `cv_m` +3. Specific latent heats of vaporization, fusion, and sublimation: + * `latent_heat_vapor` + * `latent_heat_fusion` + * `latent_heat_sublim` +4. Saturation vapor pressure and specific humidity over liquid and ice: + * `sat_vapor_press_liquid` + * `sat_vapor_press_ice` + * `sat_shum` +5. Functions computing energies and inverting them to obtain temperatures + * `total_energy` + * `internal_energy` + * `air_temperature` +6. Functions to compute temperatures and partitioning of water into phases in thermodynamic equilibrium (when Gibbs' phase rule implies that the entire thermodynamic state of moist air, including the liquid and ice specific humidities, can be calculated from the 3 thermodynamic state variables, such as energy, pressure, and total specific humidity) + * `liquid_fraction` (fraction of condensate that is liquid) + * `saturation_adjustment` (compute temperature from energy, density, and total specific humidity) +7. Auxiliary functions for diagnostic purposes, e.g., other thermodynamic quantities + * `liquid_ice_pottemp` (liquid-ice potential temperature) + +A moist dynamical core that assumes equilibrium thermodynamics can be obtained from a dry dynamical core with total energy as a prognostic variable by including a tracer for the total specific humidity `q.tot`, using the functions, e.g., for the energies in the module, and computing the temperature `T` and the liquid and ice specific humidities (`q.liq` and `q.ice`) from the internal energy `e_int` by saturation adjustment. diff --git a/env/Plots/Project.toml b/env/Plots/Project.toml new file mode 100644 index 00000000..43aec5b7 --- /dev/null +++ b/env/Plots/Project.toml @@ -0,0 +1,2 @@ +[deps] +Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" diff --git a/env/test/Project.toml b/env/test/Project.toml new file mode 100644 index 00000000..83b1b9cb --- /dev/null +++ b/env/test/Project.toml @@ -0,0 +1,5 @@ +[deps] +NCDatasets = "85f8d34a-cbdd-5861-8df4-14fed0d494ab" +Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f" +Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" +Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" diff --git a/src/MoistThermodynamics.jl b/src/MoistThermodynamics.jl index 81b8893d..b0bb365c 100644 --- a/src/MoistThermodynamics.jl +++ b/src/MoistThermodynamics.jl @@ -4,17 +4,37 @@ Moist thermodynamic functions, e.g., for air pressure (atmosphere equation of state), latent heats of phase transitions, saturation vapor pressures, and saturation specific humidities. + + +## AbstractParameterSet's + +Many functions defined in this module rely on CLIMAParameters.jl. +CLIMAParameters.jl defines several functions (e.g., many planet +parameters). For example, to compute the mole-mass ratio: + +```julia +using CLIMAParameters.Planet: molmass_ratio +using CLIMAParameters: AbstractEarthParameterSet +struct EarthParameterSet <: AbstractEarthParameterSet end +param_set = EarthParameterSet() +_molmass_ratio = molmass_ratio(param_set) +``` + +Because these parameters are widely used throughout this module, +`param_set` is an argument for many MoistThermodynamics functions. """ module MoistThermodynamics using DocStringExtensions using RootSolvers -using CLIMAParameters +using CLIMAParameters: AbstractParameterSet using CLIMAParameters.Planet +const APS = AbstractParameterSet @inline q_pt_0(::Type{FT}) where {FT} = PhasePartition{FT}(FT(0), FT(0), FT(0)) include("states.jl") +include("profiles.jl") include("isentropic.jl") include("relations.jl") diff --git a/src/isentropic.jl b/src/isentropic.jl index 2435b095..e6ea46c7 100644 --- a/src/isentropic.jl +++ b/src/isentropic.jl @@ -14,54 +14,64 @@ For dispatching to isentropic formulas struct DryAdiabaticProcess end """ - air_pressure_given_θ(θ::FT, Φ::FT, ::DryAdiabaticProcess) + air_pressure_given_θ(param_set, θ::FT, Φ::FT, ::DryAdiabaticProcess) The air pressure for an isentropic process, where + - `param_set` an `AbstractParameterSet`, see the [`MoistThermodynamics`](@ref) for more details - `θ` potential temperature - `Φ` gravitational potential """ function air_pressure_given_θ( - param_set::PS, + param_set::APS, θ::FT, Φ::FT, ::DryAdiabaticProcess, -) where {FT, PS} - return FT(MSLP(param_set)) * (1 - Φ / (θ * FT(cp_d(param_set))))^(FT(cp_d(param_set)) / FT(R_d(param_set))) +) where {FT <: AbstractFloat} + _MSLP::FT = MSLP(param_set) + _R_d::FT = R_d(param_set) + _cp_d::FT = cp_d(param_set) + return _MSLP * (1 - Φ / (θ * _cp_d))^(_cp_d / _R_d) end """ - air_pressure(T::FT, T∞::FT, p∞::FT, ::DryAdiabaticProcess) + air_pressure(param_set, T::FT, T∞::FT, p∞::FT, ::DryAdiabaticProcess) The air pressure for an isentropic process, where + - `param_set` an `AbstractParameterSet`, see the [`MoistThermodynamics`](@ref) for more details - `T` temperature - `T∞` ambient temperature - `p∞` ambient pressure """ function air_pressure( - param_set::PS, + param_set::APS, T::FT, T∞::FT, p∞::FT, ::DryAdiabaticProcess, -) where {FT, PS} - return p∞ * (T / T∞)^(FT(1) / FT(kappa_d(param_set))) +) where {FT <: AbstractFloat} + _kappa_d::FT = kappa_d(param_set) + return p∞ * (T / T∞)^(FT(1) / _kappa_d) end """ - air_temperature(p::FT, θ::FT, Φ::FT, ::DryAdiabaticProcess) + air_temperature(param_set, p::FT, θ::FT, Φ::FT, ::DryAdiabaticProcess) The air temperature for an isentropic process, where + - `param_set` an `AbstractParameterSet`, see the [`MoistThermodynamics`](@ref) for more details - `p` pressure - `θ` potential temperature """ function air_temperature( - param_set::PS, + param_set::APS, p::FT, θ::FT, ::DryAdiabaticProcess, -) where {FT, PS} - return (p / FT(MSLP(param_set)))^(FT(R_d(param_set)) / FT(cp_d(param_set))) * θ +) where {FT <: AbstractFloat} + _R_d::FT = R_d(param_set) + _cp_d::FT = cp_d(param_set) + _MSLP::FT = MSLP(param_set) + return (p / _MSLP)^(_R_d / _cp_d) * θ end diff --git a/src/profiles.jl b/src/profiles.jl new file mode 100644 index 00000000..1e090ad7 --- /dev/null +++ b/src/profiles.jl @@ -0,0 +1,81 @@ +#### Reference state profiles + +""" + fixed_lapse_rate_ref_state( + param_set::APS, + z::FT, + T_surface::FT, + T_min::FT, + ) where {FT <: AbstractFloat} + +Fixed lapse rate hydrostatic reference state +""" +function fixed_lapse_rate_ref_state( + param_set::APS, + z::FT, + T_surface::FT, + T_min::FT, +) where {FT <: AbstractFloat} + _grav::FT = grav(param_set) + _cp_d::FT = cp_d(param_set) + _R_d::FT = R_d(param_set) + _MSLP::FT = MSLP(param_set) + Γ = _grav / _cp_d + z_tropopause = (T_surface - T_min) / Γ + H_min = _R_d * T_min / _grav + T = max(T_surface - Γ * z, T_min) + p = _MSLP * (T / T_surface)^(_grav / (_R_d * Γ)) + T == T_min && (p = p * exp(-(z - z_tropopause) / H_min)) + ρ = p / (_R_d * T) + return T, p, ρ +end + +""" + tested_profiles(param_set::APS, n::Int, ::Type{FT}) where {FT} + +A range of input arguments to thermodynamic state constructors + - `param_set` an `AbstractParameterSet`, see the [`MoistThermodynamics`](@ref) for more details + - `e_int` internal energy + - `ρ` (moist-)air density + - `q_tot` total specific humidity + - `q_pt` phase partition + - `T` air temperature + - `θ_liq_ice` liquid-ice potential temperature +that are tested for convergence in saturation adjustment. + +Note that the output vectors are of size ``n*n_RH``, and they +should span the input arguments to all of the constructors. +""" +function tested_profiles(param_set::APS, n::Int, ::Type{FT}) where {FT} + n_RS1 = 10 + n_RS2 = 20 + n_RS = n_RS1 + n_RS2 + z_range = range(FT(0), stop = FT(2.5e4), length = n) + relative_sat1 = range(FT(0), stop = FT(1), length = n_RS1) + relative_sat2 = range(FT(1), stop = FT(1.02), length = n_RS2) + relative_sat = [relative_sat1..., relative_sat2...] + T_min = FT(150) + T_surface = FT(350) + + args = + fixed_lapse_rate_ref_state.( + Ref(param_set), + z_range, + Ref(T_surface), + Ref(T_min), + ) + T, p, ρ = getindex.(args, 1), getindex.(args, 2), getindex.(args, 3) + + p = collect(Iterators.flatten([p for RS in 1:n_RS])) + ρ = collect(Iterators.flatten([ρ for RS in 1:n_RS])) + T = collect(Iterators.flatten([T for RS in 1:n_RS])) + relative_sat = collect(Iterators.flatten([relative_sat for RS in 1:n])) + + # Additional variables + q_sat = q_vap_saturation.(Ref(param_set), T, ρ) + q_tot = min.(relative_sat .* q_sat, FT(1)) + q_pt = PhasePartition_equil.(Ref(param_set), T, ρ, q_tot) + e_int = internal_energy.(Ref(param_set), T, q_pt) + θ_liq_ice = liquid_ice_pottemp.(Ref(param_set), T, ρ, q_pt) + return e_int, ρ, q_tot, q_pt, T, p, θ_liq_ice +end diff --git a/src/relations.jl b/src/relations.jl index 583280a8..8ec10700 100644 --- a/src/relations.jl +++ b/src/relations.jl @@ -1,4 +1,3 @@ - # Atmospheric equation of state export air_pressure, air_temperature, air_density, specific_volume, soundspeed_air @@ -34,27 +33,24 @@ export air_temperature_from_liquid_ice_pottemp, export air_temperature_from_liquid_ice_pottemp_non_linear export vapor_specific_humidity + """ - gas_constant_air([q::PhasePartition]) + gas_constant_air(param_set, [q::PhasePartition]) The specific gas constant of moist air given + - `param_set` an `AbstractParameterSet`, see the [`MoistThermodynamics`](@ref) for more details - `q` [`PhasePartition`](@ref). Without this argument, the results are for dry air. """ -function gas_constant_air( - param_set::PS, - q::PhasePartition{FT}=q_pt_0(FT), -) where {FT, PS} - return FT(R_d(param_set)) * - (1 + (FT(molmass_ratio(param_set)) - 1) * q.tot - - FT(molmass_ratio(param_set)) * (q.liq + q.ice)) -end -function gas_constant_air( - param_set::PS, - ::Type{FT}, -) where {FT, PS} - return gas_constant_air(param_set, q_pt_0(FT)) +function gas_constant_air(param_set::APS, q::PhasePartition{FT}) where {FT} + _R_d::FT = R_d(param_set) + _molmass_ratio::FT = molmass_ratio(param_set) + return _R_d * + (1 + (_molmass_ratio - 1) * q.tot - _molmass_ratio * (q.liq + q.ice)) end +gas_constant_air(param_set::APS, ::Type{FT}) where {FT} = + gas_constant_air(param_set, q_pt_0(FT)) + """ gas_constant_air(ts::ThermodynamicState) @@ -75,23 +71,24 @@ vapor_specific_humidity(ts::ThermodynamicState) = vapor_specific_humidity(PhasePartition(ts)) """ - air_pressure(T, ρ[, q::PhasePartition]) + air_pressure(param_set, T, ρ[, q::PhasePartition]) The air pressure from the equation of state (ideal gas law) where + - `param_set` an `AbstractParameterSet`, see the [`MoistThermodynamics`](@ref) for more details - `T` air temperature - `ρ` (moist-)air density and, optionally, - `q` [`PhasePartition`](@ref). Without this argument, the results are for dry air. """ function air_pressure( - param_set::PS, + param_set::APS, T::FT, ρ::FT, q::PhasePartition{FT} = q_pt_0(FT), -) where {FT <: Real, PS} - return gas_constant_air(param_set, q) * ρ * T +) where {FT <: Real} + return gas_constant_air(param_set, q) * ρ * T end """ @@ -107,24 +104,26 @@ air_pressure(ts::ThermodynamicState) = air_pressure( PhasePartition(ts), ) + """ - air_density(T, p[, q::PhasePartition]) + air_density(param_set, T, p[, q::PhasePartition]) The (moist-)air density from the equation of state (ideal gas law) where + - `param_set` an `AbstractParameterSet`, see the [`MoistThermodynamics`](@ref) for more details - `T` air temperature - `p` pressure and, optionally, - `q` [`PhasePartition`](@ref). Without this argument, the results are for dry air. """ function air_density( - param_set::PS, + param_set::APS, T::FT, p::FT, q::PhasePartition{FT} = q_pt_0(FT), -) where {FT <: Real, PS} - return p / (gas_constant_air(param_set, q) * T) +) where {FT <: Real} + return p / (gas_constant_air(param_set, q) * T) end """ @@ -151,28 +150,26 @@ total_specific_humidity(ts::PhaseDry{FT}) where {FT} = FT(0) total_specific_humidity(ts::PhaseNonEquil) = ts.q.tot """ - cp_m([q::PhasePartition]) + cp_m(param_set, [q::PhasePartition]) The isobaric specific heat capacity of moist air where, optionally, + - `param_set` an `AbstractParameterSet`, see the [`MoistThermodynamics`](@ref) for more details - `q` [`PhasePartition`](@ref). Without this argument, the results are for dry air. """ -function cp_m( - param_set::PS, - q::PhasePartition{FT}=q_pt_0(FT), -) where {FT <: Real, PS} - return FT(cp_d(param_set)) + - (FT(cp_v(param_set)) - FT(cp_d(param_set))) * q.tot + - (FT(cp_l(param_set)) - FT(cp_v(param_set))) * q.liq + - (FT(cp_i(param_set)) - FT(cp_v(param_set))) * q.ice +function cp_m(param_set::APS, q::PhasePartition{FT}) where {FT <: Real} + _cp_d::FT = cp_d(param_set) + _cp_v::FT = cp_v(param_set) + _cp_l::FT = cp_l(param_set) + _cp_i::FT = cp_i(param_set) + return _cp_d + + (_cp_v - _cp_d) * q.tot + + (_cp_l - _cp_v) * q.liq + + (_cp_i - _cp_v) * q.ice end -function cp_m( - param_set::PS, - ::Type{FT}, -) where {FT <: Real, PS} - return cp_m(param_set, q_pt_0(FT)) -end +cp_m(param_set::APS, ::Type{FT}) where {FT <: Real} = + cp_m(param_set, q_pt_0(FT)) """ cp_m(ts::ThermodynamicState) @@ -184,28 +181,26 @@ cp_m(ts::ThermodynamicState) = cp_m(ts.param_set, PhasePartition(ts)) cp_m(ts::PhaseDry{FT}) where {FT <: Real} = FT(cp_d(ts.param_set)) """ - cv_m([q::PhasePartition]) + cv_m(param_set, [q::PhasePartition]) The isochoric specific heat capacity of moist air where optionally, + - `param_set` an `AbstractParameterSet`, see the [`MoistThermodynamics`](@ref) for more details - `q` [`PhasePartition`](@ref). Without this argument, the results are for dry air. """ -function cv_m( - param_set::PS, - q::PhasePartition{FT}=q_pt_0(FT), -) where {FT <: Real, PS} - return FT(cv_d(param_set)) + - (FT(cv_v(param_set)) - FT(cv_d(param_set))) * q.tot + - (FT(cv_l(param_set)) - FT(cv_v(param_set))) * q.liq + - (FT(cv_i(param_set)) - FT(cv_v(param_set))) * q.ice +function cv_m(param_set::APS, q::PhasePartition{FT}) where {FT <: Real} + _cv_d::FT = cv_d(param_set) + _cv_v::FT = cv_v(param_set) + _cv_l::FT = cv_l(param_set) + _cv_i::FT = cv_i(param_set) + return _cv_d + + (_cv_v - _cv_d) * q.tot + + (_cv_l - _cv_v) * q.liq + + (_cv_i - _cv_v) * q.ice end -function cv_m( - param_set::PS, - ::Type{FT}, -) where {FT <: Real, PS} - return cv_m(param_set, q_pt_0(FT)) -end +cv_m(param_set::APS, ::Type{FT}) where {FT <: Real} = + cv_m(param_set, q_pt_0(FT)) """ cv_m(ts::ThermodynamicState) @@ -218,9 +213,10 @@ cv_m(ts::PhaseDry{FT}) where {FT <: Real} = FT(cv_d(ts.param_set)) """ - (R_m, cp_m, cv_m, γ_m) = gas_constants([q::PhasePartition]) + (R_m, cp_m, cv_m, γ_m) = gas_constants(param_set, [q::PhasePartition]) Wrapper to compute all gas constants at once, where optionally, + - `param_set` an `AbstractParameterSet`, see the [`MoistThermodynamics`](@ref) for more details - `q` [`PhasePartition`](@ref). Without this argument, the results are for dry air. The function returns a tuple of @@ -230,21 +226,15 @@ The function returns a tuple of - `γ_m = cp_m/cv_m` """ function gas_constants( - param_set::PS, + param_set::APS, q::PhasePartition{FT} = q_pt_0(FT), -) where {FT <: Real, PS} +) where {FT <: Real} R_gas = gas_constant_air(param_set, q) cp = cp_m(param_set, q) cv = cv_m(param_set, q) γ = cp / cv return (R_gas, cp, cv, γ) end -function gas_constants( - param_set::PS, - ::Type{FT}, -) where {FT <: Real, PS} - return gas_constants(param_set, q_pt_0(FT)) -end """ (R_m, cp_m, cv_m, γ_m) = gas_constants(ts::ThermodynamicState) @@ -262,23 +252,26 @@ gas_constants(ts::ThermodynamicState) = gas_constants(ts.param_set, PhasePartition(ts)) """ - air_temperature(e_int, q::PhasePartition) + air_temperature(param_set, e_int, q::PhasePartition) The air temperature, where + - `param_set` an `AbstractParameterSet`, see the [`MoistThermodynamics`](@ref) for more details - `e_int` internal energy per unit mass and, optionally, - `q` [`PhasePartition`](@ref). Without this argument, the results are for dry air. """ function air_temperature( - param_set::PS, + param_set::APS, e_int::FT, q::PhasePartition{FT} = q_pt_0(FT), -) where {FT <: Real, PS} - return FT(T_0(param_set)) + - ( - e_int - (q.tot - q.liq) * FT(e_int_v0(param_set)) + - q.ice * (FT(e_int_v0(param_set)) + FT(e_int_i0(param_set))) +) where {FT <: Real} + _T_0::FT = T_0(param_set) + _e_int_v0::FT = e_int_v0(param_set) + _e_int_i0::FT = e_int_i0(param_set) + return _T_0 + + ( + e_int - (q.tot - q.liq) * _e_int_v0 + q.ice * (_e_int_v0 + _e_int_i0) ) / cv_m(param_set, q) end @@ -293,21 +286,25 @@ air_temperature(ts::PhaseEquil) = ts.T """ - internal_energy(T[, q::PhasePartition]) + internal_energy(param_set, T[, q::PhasePartition]) The internal energy per unit mass, given a thermodynamic state `ts` or + - `param_set` an `AbstractParameterSet`, see the [`MoistThermodynamics`](@ref) for more details - `T` temperature and, optionally, - `q` [`PhasePartition`](@ref). Without this argument, the results are for dry air. """ function internal_energy( - param_set::PS, + param_set::APS, T::FT, q::PhasePartition{FT} = q_pt_0(FT), -) where {FT <: Real, PS} - return cv_m(param_set, q) * (T - FT(T_0(param_set))) + (q.tot - q.liq) * FT(e_int_v0(param_set)) - - q.ice * (FT(e_int_v0(param_set)) + FT(e_int_i0(param_set))) +) where {FT <: Real} + _T_0::FT = T_0(param_set) + _e_int_v0::FT = e_int_v0(param_set) + _e_int_i0::FT = e_int_i0(param_set) + return cv_m(param_set, q) * (T - _T_0) + (q.tot - q.liq) * _e_int_v0 - + q.ice * (_e_int_v0 + _e_int_i0) end """ @@ -341,21 +338,26 @@ The internal energy per unit mass, given end """ - internal_energy_sat(T, ρ, q_tot) + internal_energy_sat(param_set, T, ρ, q_tot) The internal energy per unit mass in thermodynamic equilibrium at saturation where + - `param_set` an `AbstractParameterSet`, see the [`MoistThermodynamics`](@ref) for more details - `T` temperature - `ρ` (moist-)air density - `q_tot` total specific humidity """ function internal_energy_sat( - param_set::PS, + param_set::APS, T::FT, ρ::FT, q_tot::FT, -) where {FT <: Real, PS} - return internal_energy(param_set, T, PhasePartition_equil(param_set, T, ρ, q_tot)) +) where {FT <: Real} + return internal_energy( + param_set, + T, + PhasePartition_equil(param_set, T, ρ, q_tot), + ) end """ @@ -374,10 +376,11 @@ internal_energy_sat(ts::ThermodynamicState) = internal_energy_sat( """ - total_energy(e_kin, e_pot, T[, q::PhasePartition]) + total_energy(param_set, e_kin, e_pot, T[, q::PhasePartition]) The total energy per unit mass, given + - `param_set` an `AbstractParameterSet`, see the [`MoistThermodynamics`](@ref) for more details - `e_kin` kinetic energy per unit mass - `e_pot` potential energy per unit mass - `T` temperature @@ -386,13 +389,13 @@ and, optionally, """ function total_energy( - param_set::PS, + param_set::APS, e_kin::FT, e_pot::FT, T::FT, q::PhasePartition{FT} = q_pt_0(FT), -) where {FT <: Real, PS} - return e_kin + e_pot + internal_energy(param_set, T, q) +) where {FT <: Real} + return e_kin + e_pot + internal_energy(param_set, T, q) end """ @@ -410,18 +413,19 @@ function total_energy( end """ - soundspeed_air(T[, q::PhasePartition]) + soundspeed_air(param_set, T[, q::PhasePartition]) The speed of sound in unstratified air, where + - `param_set` an `AbstractParameterSet`, see the [`MoistThermodynamics`](@ref) for more details - `T` temperature and, optionally, - `q` [`PhasePartition`](@ref). Without this argument, the results are for dry air. """ function soundspeed_air( - param_set::PS, + param_set::APS, T::FT, q::PhasePartition{FT} = q_pt_0(FT), -) where {FT <: Real, PS} +) where {FT <: Real} γ = cp_m(param_set, q) / cv_m(param_set, q) R_m = gas_constant_air(param_set, q) return sqrt(γ * R_m * T) @@ -437,13 +441,18 @@ soundspeed_air(ts::ThermodynamicState) = """ - latent_heat_vapor(T::FT) where {FT<:Real} + latent_heat_vapor(param_set, T::FT) where {FT<:Real} The specific latent heat of vaporization where + - `param_set` an `AbstractParameterSet`, see the [`MoistThermodynamics`](@ref) for more details - `T` temperature """ -latent_heat_vapor(param_set::PS, T::FT) where {FT <: Real, PS} = - latent_heat_generic(param_set, T, FT(LH_v0(param_set)), FT(cp_v(param_set)) - FT(cp_l(param_set))) +function latent_heat_vapor(param_set::APS, T::FT) where {FT <: Real} + _cp_l::FT = cp_l(param_set) + _cp_v::FT = cp_v(param_set) + _LH_v0::FT = LH_v0(param_set) + return latent_heat_generic(param_set, T, _LH_v0, _cp_v - _cp_l) +end """ latent_heat_vapor(ts::ThermodynamicState) @@ -455,13 +464,18 @@ latent_heat_vapor(ts::ThermodynamicState) = latent_heat_vapor(ts.param_set, air_temperature(ts)) """ - latent_heat_sublim(T::FT) where {FT<:Real} + latent_heat_sublim(param_set, T::FT) where {FT<:Real} The specific latent heat of sublimation where + - `param_set` an `AbstractParameterSet`, see the [`MoistThermodynamics`](@ref) for more details - `T` temperature """ -latent_heat_sublim(param_set::PS, T::FT) where {FT <: Real, PS} = - latent_heat_generic(param_set, T, FT(LH_s0(param_set)), FT(cp_v(param_set)) - FT(cp_i(param_set))) +function latent_heat_sublim(param_set::APS, T::FT) where {FT <: Real} + _LH_s0::FT = LH_s0(param_set) + _cp_v::FT = cp_v(param_set) + _cp_i::FT = cp_i(param_set) + return latent_heat_generic(param_set, T, _LH_s0, _cp_v - _cp_i) +end """ latent_heat_sublim(ts::ThermodynamicState) @@ -473,13 +487,18 @@ latent_heat_sublim(ts::ThermodynamicState) = latent_heat_sublim(ts.param_set, air_temperature(ts)) """ - latent_heat_fusion(T::FT) where {FT<:Real} + latent_heat_fusion(param_set, T::FT) where {FT<:Real} The specific latent heat of fusion where + - `param_set` an `AbstractParameterSet`, see the [`MoistThermodynamics`](@ref) for more details - `T` temperature """ -latent_heat_fusion(param_set::PS, T::FT) where {FT <: Real, PS} = - latent_heat_generic(param_set, T, FT(LH_f0(param_set)), FT(cp_l(param_set)) - FT(cp_i(param_set))) +function latent_heat_fusion(param_set::APS, T::FT) where {FT <: Real} + _LH_f0::FT = LH_f0(param_set) + _cp_l::FT = cp_l(param_set) + _cp_i::FT = cp_i(param_set) + return latent_heat_generic(param_set, T, _LH_f0, _cp_l - _cp_i) +end """ latent_heat_fusion(ts::ThermodynamicState) @@ -487,16 +506,17 @@ latent_heat_fusion(param_set::PS, T::FT) where {FT <: Real, PS} = The specific latent heat of fusion given a thermodynamic state `ts`. """ -latent_heat_fusion(ts::ThermodynamicState{FT}) where {FT <: Real} = +latent_heat_fusion(ts::ThermodynamicState) = latent_heat_fusion(ts.param_set, air_temperature(ts)) """ - latent_heat_generic(T::FT, LH_0::FT, Δcp::FT) where {FT<:Real} + latent_heat_generic(param_set, T::FT, LH_0::FT, Δcp::FT) where {FT<:Real} The specific latent heat of a generic phase transition between two phases, computed using Kirchhoff's relation with constant isobaric specific heat capacities of the two phases, given + - `param_set` an `AbstractParameterSet`, see the [`MoistThermodynamics`](@ref) for more details - `T` temperature - `LH_0` latent heat of the phase transition at `T_0` - `Δcp` difference between the isobaric specific heat capacities @@ -504,12 +524,13 @@ isobaric specific heat capacities of the two phases, given in the lower-temperature phase). """ function latent_heat_generic( - param_set::PS, + param_set::APS, T::FT, LH_0::FT, Δcp::FT, -) where {FT <: Real, PS} - return LH_0 + Δcp * (T - FT(T_0(param_set))) +) where {FT <: Real} + _T_0::FT = T_0(param_set) + return LH_0 + Δcp * (T - _T_0) end @@ -541,17 +562,19 @@ An ice phase, to dispatch over struct Ice <: Phase end """ - saturation_vapor_pressure(T, Liquid()) + saturation_vapor_pressure(param_set, T, Liquid()) -Return the saturation vapor pressure over a plane liquid surface at -temperature `T`. +Return the saturation vapor pressure over a plane liquid surface given + - `T` temperature + - `param_set` an `AbstractParameterSet`, see the [`MoistThermodynamics`](@ref) for more details - saturation_vapor_pressure(T, Ice()) + saturation_vapor_pressure(param_set, T, Ice()) -Return the saturation vapor pressure over a plane ice surface at -temperature `T`. +Return the saturation vapor pressure over a plane ice surface given + - `T` temperature + - `param_set` an `AbstractParameterSet`, see the [`MoistThermodynamics`](@ref) for more details - saturation_vapor_pressure(T, LH_0, Δcp) + saturation_vapor_pressure(param_set, T, LH_0, Δcp) Compute the saturation vapor pressure over a plane surface by integration of the Clausius-Clapeyron relation. @@ -573,65 +596,81 @@ the triple point pressure `press_triple`. """ function saturation_vapor_pressure( - param_set::PS, + param_set::APS, T::FT, ::Liquid, -) where {FT <: Real, PS} - return saturation_vapor_pressure(param_set, T, FT(LH_v0(param_set)), FT(cp_v(param_set)) - FT(cp_l(param_set))) +) where {FT <: Real} + _LH_v0::FT = LH_v0(param_set) + _cp_v::FT = cp_v(param_set) + _cp_l::FT = cp_l(param_set) + return saturation_vapor_pressure(param_set, T, _LH_v0, _cp_v - _cp_l) end function saturation_vapor_pressure( ts::ThermodynamicState{FT}, ::Liquid, ) where {FT <: Real} - + _LH_v0::FT = LH_v0(ts.param_set) + _cp_v::FT = cp_v(ts.param_set) + _cp_l::FT = cp_l(ts.param_set) return saturation_vapor_pressure( ts.param_set, air_temperature(ts), - FT(LH_v0(ts.param_set)), - FT(cp_v(ts.param_set)) - FT(cp_l(ts.param_set)), + _LH_v0, + _cp_v - _cp_l, ) end function saturation_vapor_pressure( - param_set::PS, + param_set::APS, T::FT, ::Ice, -) where {FT <: Real, PS} - return saturation_vapor_pressure(param_set, T, FT(LH_s0(param_set)), FT(cp_v(param_set)) - FT(cp_i(param_set))) +) where {FT <: Real} + _LH_s0::FT = LH_s0(param_set) + _cp_v::FT = cp_v(param_set) + _cp_i::FT = cp_i(param_set) + return saturation_vapor_pressure(param_set, T, _LH_s0, _cp_v - _cp_i) end function saturation_vapor_pressure( ts::ThermodynamicState{FT}, ::Ice, -) where {FT <: Real, PS} +) where {FT <: Real} + _LH_s0::FT = LH_s0(ts.param_set) + _cp_v::FT = cp_v(ts.param_set) + _cp_i::FT = cp_i(ts.param_set) return saturation_vapor_pressure( - ts.param_set, - air_temperature(ts), - FT(LH_s0(ts.param_set)), - FT(cp_v(ts.param_set)) - FT(cp_i(ts.param_set)), -) + ts.param_set, + air_temperature(ts), + _LH_s0, + _cp_v - _cp_i, + ) end function saturation_vapor_pressure( - param_set::PS, + param_set::APS, T::FT, LH_0::FT, Δcp::FT, -) where {FT <: Real, PS} +) where {FT <: Real} + _press_triple::FT = press_triple(param_set) + _R_v::FT = R_v(param_set) + _T_triple::FT = T_triple(param_set) + _T_0::FT = T_0(param_set) - return FT(press_triple(param_set)) * - (T / FT(T_triple(param_set)))^(Δcp / FT(R_v(param_set))) * - exp((FT(LH_0) - Δcp * FT(T_0(param_set))) / FT(R_v(param_set)) * (1 / FT(T_triple(param_set)) - 1 / T)) + return _press_triple * + (T / _T_triple)^(Δcp / _R_v) * + exp((LH_0 - Δcp * _T_0) / _R_v * (1 / _T_triple - 1 / T)) end """ - q_vap_saturation_generic(T, ρ[; phase=Liquid()]) + q_vap_saturation_generic(param_set, T, ρ[; phase=Liquid()]) Compute the saturation specific humidity over a plane surface of condensate, given + - `param_set` an `AbstractParameterSet`, see the [`MoistThermodynamics`](@ref) for more details - `T` temperature - `ρ` (moist-)air density and, optionally, @@ -639,20 +678,21 @@ and, optionally, - `Ice()` indicating condensate is ice """ function q_vap_saturation_generic( - param_set::PS, + param_set::APS, T::FT, ρ::FT; phase::Phase = Liquid(), -) where {FT <: Real, PS} +) where {FT <: Real} p_v_sat = saturation_vapor_pressure(param_set, T, phase) return q_vap_saturation_from_pressure(param_set, T, ρ, p_v_sat) end """ - q_vap_saturation(T, ρ[, q::PhasePartition]) + q_vap_saturation(param_set, T, ρ[, q::PhasePartition]) Compute the saturation specific humidity, given + - `param_set` an `AbstractParameterSet`, see the [`MoistThermodynamics`](@ref) for more details - `T` temperature - `ρ` (moist-)air density and, optionally, @@ -672,21 +712,25 @@ fraction of liquid given by temperature dependent `liquid_fraction(T)` and the fraction of ice by the complement `1 - liquid_fraction(T)`. """ function q_vap_saturation( - param_set::PS, + param_set::APS, T::FT, ρ::FT, q::PhasePartition{FT} = q_pt_0(FT), -) where {FT <: Real, PS} +) where {FT <: Real} + _LH_v0::FT = LH_v0(param_set) + _LH_s0::FT = LH_s0(param_set) + _cp_v::FT = cp_v(param_set) + _cp_l::FT = cp_l(param_set) + _cp_i::FT = cp_i(param_set) # get phase partitioning _liquid_frac = liquid_fraction(param_set, T, q) _ice_frac = 1 - _liquid_frac # effective latent heat at T_0 and effective difference in isobaric specific # heats of the mixture - LH_0 = _liquid_frac * FT(LH_v0(param_set)) + _ice_frac * FT(LH_s0(param_set)) - Δcp = - _liquid_frac * (FT(cp_v(param_set)) - FT(cp_l(param_set))) + _ice_frac * (FT(cp_v(param_set)) - FT(cp_i(param_set))) + LH_0 = _liquid_frac * _LH_v0 + _ice_frac * _LH_s0 + Δcp = _liquid_frac * (_cp_v - _cp_l) + _ice_frac * (_cp_v - _cp_i) # saturation vapor pressure over possible mixture of liquid and ice p_v_sat = saturation_vapor_pressure(param_set, T, LH_0, Δcp) @@ -708,28 +752,31 @@ q_vap_saturation(ts::ThermodynamicState) = q_vap_saturation( ) """ - q_vap_saturation_from_pressure(T, ρ, p_v_sat) + q_vap_saturation_from_pressure(param_set, T, ρ, p_v_sat) Compute the saturation specific humidity, given + - `param_set` an `AbstractParameterSet`, see the [`MoistThermodynamics`](@ref) for more details - `T` temperature, - `ρ` (moist-)air density - `p_v_sat` saturation vapor pressure """ function q_vap_saturation_from_pressure( - param_set::PS, + param_set::APS, T::FT, ρ::FT, p_v_sat::FT, -) where {FT <: Real, PS} - return p_v_sat / (ρ * FT(R_v(param_set)) * T) +) where {FT <: Real} + _R_v::FT = R_v(param_set) + return p_v_sat / (ρ * _R_v * T) end """ - saturation_excess(T, ρ, q::PhasePartition) + saturation_excess(param_set, T, ρ, q::PhasePartition) The saturation excess in equilibrium where + - `param_set` an `AbstractParameterSet`, see the [`MoistThermodynamics`](@ref) for more details - `T` temperature - `ρ` (moist-)air density - `q` [`PhasePartition`](@ref) @@ -739,11 +786,11 @@ and the saturation specific humidity in equilibrium, and it is defined to be nonzero only if this difference is positive. """ function saturation_excess( - param_set::PS, + param_set::APS, T::FT, ρ::FT, q::PhasePartition{FT}, -) where {FT <: Real, PS} +) where {FT <: Real} return max(0, q.tot - q_vap_saturation(param_set, T, ρ, q)) end @@ -761,10 +808,11 @@ saturation_excess(ts::ThermodynamicState) = saturation_excess( ) """ - liquid_fraction(T[, q::PhasePartition]) + liquid_fraction(param_set, T[, q::PhasePartition]) The fraction of condensate that is liquid where + - `param_set` an `AbstractParameterSet`, see the [`MoistThermodynamics`](@ref) for more details - `T` temperature - `q` [`PhasePartition`](@ref) @@ -775,17 +823,18 @@ Otherwise, phase equilibrium is assumed so that the fraction of liquid is a function that is 1 above `T_freeze` and goes to zero below `T_freeze`. """ function liquid_fraction( - param_set::PS, + param_set::APS, T::FT, q::PhasePartition{FT} = q_pt_0(FT), -) where {FT <: Real, PS} +) where {FT <: Real} + _T_freeze::FT = T_freeze(param_set) q_c = q.liq + q.ice # condensate specific humidity if q_c > 0 return q.liq / q_c else # For now: Heaviside function for partitioning into liquid and ice: all liquid # for T > T_freeze; all ice for T <= T_freeze - return FT(T > FT(T_freeze(param_set))) + return FT(T > _T_freeze) end end @@ -798,11 +847,12 @@ liquid_fraction(ts::ThermodynamicState) = liquid_fraction(ts.param_set, air_temperature(ts), PhasePartition(ts)) """ - PhasePartition_equil(T, ρ, q_tot) + PhasePartition_equil(param_set, T, ρ, q_tot) Partition the phases in equilibrium, returning a [`PhasePartition`](@ref) object using the [`liquid_fraction`](@ref) function where + - `param_set` an `AbstractParameterSet`, see the [`MoistThermodynamics`](@ref) for more details - `T` temperature - `ρ` (moist-)air density - `q_tot` total specific humidity @@ -810,11 +860,11 @@ Partition the phases in equilibrium, returning a [`PhasePartition`](@ref) object The residual `q.tot - q.liq - q.ice` is the vapor specific humidity. """ function PhasePartition_equil( - param_set::PS, + param_set::APS, T::FT, ρ::FT, q_tot::FT, -) where {FT <: Real, PS} +) where {FT <: Real} _liquid_frac = liquid_fraction(param_set, T) # fraction of condensate that is liquid q_c = saturation_excess(param_set, T, ρ, PhasePartition(q_tot)) # condensate specific humidity q_liq = _liquid_frac * q_c # liquid specific humidity @@ -840,29 +890,39 @@ PhasePartition(ts::PhaseEquil) = PhasePartition_equil( PhasePartition(ts::PhaseNonEquil) = ts.q function ∂e_int_∂T( - param_set::PS, + param_set::APS, T::FT, e_int::FT, ρ::FT, q_tot::FT, -) where {FT <: Real, PS} +) where {FT <: Real} + _LH_v0::FT = LH_v0(param_set) + _LH_s0::FT = LH_s0(param_set) + _R_v::FT = R_v(param_set) + _T_0::FT = T_0(param_set) + _cv_v::FT = cv_v(param_set) + _cv_l::FT = cv_l(param_set) + _cv_i::FT = cv_i(param_set) + _e_int_v0::FT = e_int_v0(param_set) + _e_int_i0::FT = e_int_i0(param_set) + cvm = cv_m(param_set, PhasePartition_equil(param_set, T, ρ, q_tot)) q_vap_sat = q_vap_saturation(param_set, T, ρ) λ = liquid_fraction(param_set, T) - L = λ * FT(LH_v0(param_set)) + (1 - λ) * FT(LH_s0(param_set)) - ∂q_vap_sat_∂T = q_vap_sat * L / (FT(R_v(param_set)) * T^2) - T0 = FT(T_0(param_set)) - dcvm_dq_vap = FT(cv_v(param_set)) - λ * FT(cv_l(param_set)) - (1 - λ) * FT(cv_i(param_set)) + L = λ * _LH_v0 + (1 - λ) * _LH_s0 + ∂q_vap_sat_∂T = q_vap_sat * L / (_R_v * T^2) + dcvm_dq_vap = _cv_v - λ * _cv_l - (1 - λ) * _cv_i return cvm + - (FT(e_int_v0(param_set)) + (1 - λ) * FT(e_int_i0(param_set)) + (T - T0) * dcvm_dq_vap) * + (_e_int_v0 + (1 - λ) * _e_int_i0 + (T - _T_0) * dcvm_dq_vap) * ∂q_vap_sat_∂T end """ - saturation_adjustment(e_int, ρ, q_tot) + saturation_adjustment(param_set, e_int, ρ, q_tot) Compute the temperature that is consistent with + - `param_set` an `AbstractParameterSet`, see the [`MoistThermodynamics`](@ref) for more details - `e_int` internal energy - `ρ` (moist-)air density - `q_tot` total specific humidity @@ -871,34 +931,33 @@ Compute the temperature that is consistent with by finding the root of -``e_int - internal_energy_sat(T,ρ,q_tot) = 0`` +``e_int - internal_energy_sat(param_set, T, ρ, q_tot) = 0`` using Newtons method with analytic gradients. See also [`saturation_adjustment`](@ref). """ function saturation_adjustment( - param_set::PS, + param_set::APS, e_int::FT, ρ::FT, q_tot::FT, maxiter::Int, tol::FT, -) where {FT <: Real, PS} - T_1 = - max(FT(T_min(param_set)), air_temperature(param_set, e_int, PhasePartition(q_tot))) # Assume all vapor +) where {FT <: Real} + _T_min::FT = T_min(param_set) + + T_1 = max(_T_min, air_temperature(param_set, e_int, PhasePartition(q_tot))) # Assume all vapor q_v_sat = q_vap_saturation(param_set, T_1, ρ) unsaturated = q_tot <= q_v_sat - if unsaturated && T_1 > FT(T_min(param_set)) + if unsaturated && T_1 > _T_min return T_1 else sol = find_zero( T -> internal_energy_sat(param_set, T, ρ, q_tot) - e_int, - T_ -> ∂e_int_∂T(param_set, T_, e_int, ρ, q_tot), - T_1, - NewtonsMethod(), + NewtonsMethod(T_1, T_ -> ∂e_int_∂T(param_set, T_, e_int, ρ, q_tot)), CompactSolution(), - tol, + SolutionTolerance(tol), maxiter, ) if !sol.converged @@ -934,10 +993,11 @@ saturation adjustment using Secant method end """ - saturation_adjustment_SecantMethod(e_int, ρ, q_tot) + saturation_adjustment_SecantMethod(param_set, e_int, ρ, q_tot) Compute the temperature `T` that is consistent with + - `param_set` an `AbstractParameterSet`, see the [`MoistThermodynamics`](@ref) for more details - `e_int` internal energy - `ρ` (moist-)air density - `q_tot` total specific humidity @@ -946,23 +1006,23 @@ Compute the temperature `T` that is consistent with by finding the root of -``e_int - internal_energy_sat(T,ρ,q_tot) = 0`` +``e_int - internal_energy_sat(param_set, T, ρ, q_tot) = 0`` See also [`saturation_adjustment_q_tot_θ_liq_ice`](@ref). """ function saturation_adjustment_SecantMethod( - param_set::PS, + param_set::APS, e_int::FT, ρ::FT, q_tot::FT, maxiter::Int, tol::FT, -) where {FT <: Real, PS} - T_1 = - max(FT(T_min(param_set)), air_temperature(param_set, e_int, PhasePartition(q_tot))) # Assume all vapor +) where {FT <: Real} + _T_min::FT = T_min(param_set) + T_1 = max(_T_min, air_temperature(param_set, e_int, PhasePartition(q_tot))) # Assume all vapor q_v_sat = q_vap_saturation(param_set, T_1, ρ) unsaturated = q_tot <= q_v_sat - if unsaturated && T_1 > FT(T_min(param_set)) + if unsaturated && T_1 > _T_min return T_1 else # FIXME here: need to revisit bounds for saturation adjustment to guarantee bracketing of zero. @@ -974,11 +1034,9 @@ function saturation_adjustment_SecantMethod( T_2 = bound_upper_temperature(T_1, T_2) sol = find_zero( T -> internal_energy_sat(param_set, T, ρ, q_tot) - e_int, - T_1, - T_2, - SecantMethod(), + SecantMethod(T_1, T_2), CompactSolution(), - tol, + SolutionTolerance(tol), maxiter, ) if !sol.converged @@ -989,10 +1047,11 @@ function saturation_adjustment_SecantMethod( end """ - saturation_adjustment_q_tot_θ_liq_ice(θ_liq_ice, ρ, q_tot) + saturation_adjustment_q_tot_θ_liq_ice(param_set, θ_liq_ice, ρ, q_tot) Compute the temperature `T` that is consistent with + - `param_set` an `AbstractParameterSet`, see the [`MoistThermodynamics`](@ref) for more details - `θ_liq_ice` liquid-ice potential temperature - `q_tot` total specific humidity - `ρ` (moist-)air density @@ -1002,21 +1061,22 @@ Compute the temperature `T` that is consistent with by finding the root of `` - θ_{liq_ice} - liquid_ice_pottemp_sat(T, ρ, q_tot) = 0 + θ_{liq_ice} - liquid_ice_pottemp_sat(param_set, T, ρ, q_tot) = 0 `` See also [`saturation_adjustment`](@ref). """ function saturation_adjustment_q_tot_θ_liq_ice( - param_set::PS, + param_set::APS, θ_liq_ice::FT, ρ::FT, q_tot::FT, maxiter::Int, tol::FT, -) where {FT <: Real, PS} +) where {FT <: Real} + _T_min::FT = T_min(param_set) T_1 = max( - FT(T_min(param_set)), + _T_min, air_temperature_from_liquid_ice_pottemp( param_set, θ_liq_ice, @@ -1026,7 +1086,7 @@ function saturation_adjustment_q_tot_θ_liq_ice( ) # Assume all vapor q_v_sat = q_vap_saturation(param_set, T_1, ρ) unsaturated = q_tot <= q_v_sat - if unsaturated && T_1 > FT(T_min(param_set)) + if unsaturated && T_1 > _T_min return T_1 else T_2 = air_temperature_from_liquid_ice_pottemp( @@ -1038,11 +1098,9 @@ function saturation_adjustment_q_tot_θ_liq_ice( T_2 = bound_upper_temperature(T_1, T_2) sol = find_zero( T -> liquid_ice_pottemp_sat(param_set, T, ρ, q_tot) - θ_liq_ice, - T_1, - T_2, - SecantMethod(), + SecantMethod(T_1, T_2), CompactSolution(), - tol, + SolutionTolerance(tol), maxiter, ) if !sol.converged @@ -1053,10 +1111,11 @@ function saturation_adjustment_q_tot_θ_liq_ice( end """ - saturation_adjustment_q_tot_θ_liq_ice_given_pressure(θ_liq_ice, p, q_tot) + saturation_adjustment_q_tot_θ_liq_ice_given_pressure(param_set, θ_liq_ice, p, q_tot) Compute the temperature `T` that is consistent with + - `param_set` an `AbstractParameterSet`, see the [`MoistThermodynamics`](@ref) for more details - `θ_liq_ice` liquid-ice potential temperature - `q_tot` total specific humidity - `p` pressure @@ -1066,19 +1125,20 @@ Compute the temperature `T` that is consistent with by finding the root of `` - θ_{liq_ice} - liquid_ice_pottemp_sat(T, air_density(T, p, PhasePartition(q_tot)), q_tot) = 0 + θ_{liq_ice} - liquid_ice_pottemp_sat(param_set, T, air_density(param_set, T, p, PhasePartition(q_tot)), q_tot) = 0 `` See also [`saturation_adjustment`](@ref). """ function saturation_adjustment_q_tot_θ_liq_ice_given_pressure( - param_set::PS, + param_set::APS, θ_liq_ice::FT, p::FT, q_tot::FT, maxiter::Int, tol::FT, -) where {FT <: Real, PS} +) where {FT <: Real} + _T_min::FT = T_min(param_set) T_1 = air_temperature_from_liquid_ice_pottemp_given_pressure( param_set, θ_liq_ice, @@ -1088,7 +1148,7 @@ function saturation_adjustment_q_tot_θ_liq_ice_given_pressure( ρ = air_density(param_set, T_1, p, PhasePartition(q_tot)) q_v_sat = q_vap_saturation(param_set, T_1, ρ) unsaturated = q_tot <= q_v_sat - if unsaturated && T_1 > FT(T_min(param_set)) + if unsaturated && T_1 > _T_min return T_1 else T_2 = air_temperature_from_liquid_ice_pottemp( @@ -1106,11 +1166,9 @@ function saturation_adjustment_q_tot_θ_liq_ice_given_pressure( air_density(param_set, T, p, PhasePartition(q_tot)), q_tot, ) - θ_liq_ice, - T_1, - T_2, - SecantMethod(), + SecantMethod(T_1, T_2), CompactSolution(), - tol, + SolutionTolerance(tol), maxiter, ) if !sol.converged @@ -1121,34 +1179,39 @@ function saturation_adjustment_q_tot_θ_liq_ice_given_pressure( end """ - latent_heat_liq_ice(q::PhasePartition{FT}) + latent_heat_liq_ice(param_set, q::PhasePartition{FT}) Effective latent heat of condensate (weighted sum of liquid and ice), -with specific latent heat evaluated at reference temperature `T_0`. +with specific latent heat evaluated at reference temperature `T_0` given + - `param_set` an `AbstractParameterSet`, see the [`MoistThermodynamics`](@ref) for more details + - `q` [`PhasePartition`](@ref). Without this argument, the results are for dry air. """ function latent_heat_liq_ice( - param_set::PS, + param_set::APS, q::PhasePartition{FT} = q_pt_0(FT), -) where {FT <: Real, PS} - return FT(LH_v0(param_set)) * q.liq + FT(LH_s0(param_set)) * q.ice +) where {FT <: Real} + _LH_v0::FT = LH_v0(param_set) + _LH_s0::FT = LH_s0(param_set) + return _LH_v0 * q.liq + _LH_s0 * q.ice end """ - liquid_ice_pottemp_given_pressure(T, p, q::PhasePartition) + liquid_ice_pottemp_given_pressure(param_set, T, p, q::PhasePartition) The liquid-ice potential temperature where + - `param_set` an `AbstractParameterSet`, see the [`MoistThermodynamics`](@ref) for more details - `T` temperature - `p` pressure and, optionally, - `q` [`PhasePartition`](@ref). Without this argument, the results are for dry air. """ function liquid_ice_pottemp_given_pressure( - param_set::PS, + param_set::APS, T::FT, p::FT, q::PhasePartition{FT} = q_pt_0(FT), -) where {FT <: Real, PS} +) where {FT <: Real} # liquid-ice potential temperature, approximating latent heats # of phase transitions as constants return dry_pottemp_given_pressure(param_set, T, p, q) * @@ -1157,26 +1220,27 @@ end """ - liquid_ice_pottemp(T, ρ, q::PhasePartition) + liquid_ice_pottemp(param_set, T, ρ, q::PhasePartition) The liquid-ice potential temperature where + - `param_set` an `AbstractParameterSet`, see the [`MoistThermodynamics`](@ref) for more details - `T` temperature - `ρ` (moist-)air density and, optionally, - `q` [`PhasePartition`](@ref). Without this argument, the results are for dry air. """ function liquid_ice_pottemp( - param_set::PS, + param_set::APS, T::FT, ρ::FT, q::PhasePartition{FT} = q_pt_0(FT), -) where {FT <: Real, PS} - return liquid_ice_pottemp_given_pressure( - param_set, - T, - air_pressure(param_set, T, ρ, q), - q, -) +) where {FT <: Real} + return liquid_ice_pottemp_given_pressure( + param_set, + T, + air_pressure(param_set, T, ρ, q), + q, + ) end """ @@ -1193,41 +1257,43 @@ liquid_ice_pottemp(ts::ThermodynamicState) = liquid_ice_pottemp( ) """ - dry_pottemp(T, ρ[, q::PhasePartition]) + dry_pottemp(param_set, T, ρ[, q::PhasePartition]) The dry potential temperature where + - `param_set` an `AbstractParameterSet`, see the [`MoistThermodynamics`](@ref) for more details - `T` temperature - `ρ` (moist-)air density and, optionally, - `q` [`PhasePartition`](@ref). Without this argument, the results are for dry air. """ function dry_pottemp( - param_set::PS, + param_set::APS, T::FT, ρ::FT, q::PhasePartition{FT} = q_pt_0(FT), -) where {FT <: Real, PS} - return T / exner(param_set, T, ρ, q) +) where {FT <: Real} + return T / exner(param_set, T, ρ, q) end """ - dry_pottemp_given_pressure(T, p[, q::PhasePartition]) + dry_pottemp_given_pressure(param_set, T, p[, q::PhasePartition]) The dry potential temperature where + - `param_set` an `AbstractParameterSet`, see the [`MoistThermodynamics`](@ref) for more details - `T` temperature - `p` pressure and, optionally, - `q` [`PhasePartition`](@ref). Without this argument, the results are for dry air. """ function dry_pottemp_given_pressure( - param_set::PS, + param_set::APS, T::FT, p::FT, q::PhasePartition{FT} = q_pt_0(FT), -) where {FT <: Real, PS} - return T / exner_given_pressure(param_set, p, q) +) where {FT <: Real} + return T / exner_given_pressure(param_set, p, q) end """ @@ -1243,36 +1309,39 @@ dry_pottemp(ts::ThermodynamicState) = dry_pottemp( ) """ - air_temperature_from_liquid_ice_pottemp(θ_liq_ice, ρ, q::PhasePartition) + air_temperature_from_liquid_ice_pottemp(param_set, θ_liq_ice, ρ, q::PhasePartition) The temperature given + - `param_set` an `AbstractParameterSet`, see the [`MoistThermodynamics`](@ref) for more details - `θ_liq_ice` liquid-ice potential temperature - `ρ` (moist-)air density and, optionally, - `q` [`PhasePartition`](@ref). Without this argument, the results are for dry air. """ function air_temperature_from_liquid_ice_pottemp( - param_set::PS, + param_set::APS, θ_liq_ice::FT, ρ::FT, q::PhasePartition{FT} = q_pt_0(FT), -) where {FT <: Real, PS} +) where {FT <: Real} + _MSLP::FT = MSLP(param_set) cvm = cv_m(param_set, q) cpm = cp_m(param_set, q) R_m = gas_constant_air(param_set, q) κ = 1 - cvm / cpm - T_u = (ρ * R_m * θ_liq_ice / FT(MSLP(param_set)))^(R_m / cvm) * θ_liq_ice + T_u = (ρ * R_m * θ_liq_ice / _MSLP)^(R_m / cvm) * θ_liq_ice T_1 = latent_heat_liq_ice(param_set, q) / cvm T_2 = -κ / (2 * T_u) * (latent_heat_liq_ice(param_set, q) / cvm)^2 return T_u + T_1 + T_2 end """ - air_temperature_from_liquid_ice_pottemp_non_linear(θ_liq_ice, ρ, q::PhasePartition) + air_temperature_from_liquid_ice_pottemp_non_linear(param_set, θ_liq_ice, ρ, q::PhasePartition) Computes temperature `T` given + - `param_set` an `AbstractParameterSet`, see the [`MoistThermodynamics`](@ref) for more details - `θ_liq_ice` liquid-ice potential temperature - `ρ` (moist-)air density - `tol` tolerance for non-linear equation solve @@ -1282,17 +1351,19 @@ and, optionally, by finding the root of `` - T - air_temperature_from_liquid_ice_pottemp_given_pressure(θ_liq_ice, air_pressure(T, ρ, q), q) = 0 + T - air_temperature_from_liquid_ice_pottemp_given_pressure(param_set, θ_liq_ice, air_pressure(param_set, T, ρ, q), q) = 0 `` """ function air_temperature_from_liquid_ice_pottemp_non_linear( - param_set::PS, + param_set::APS, θ_liq_ice::FT, ρ::FT, maxiter::Int, tol::FT, q::PhasePartition{FT} = q_pt_0(FT), -) where {FT <: Real, PS} +) where {FT <: Real} + _T_min::FT = T_min(param_set) + _T_max::FT = T_max(param_set) sol = find_zero( T -> T - air_temperature_from_liquid_ice_pottemp_given_pressure( @@ -1301,11 +1372,9 @@ function air_temperature_from_liquid_ice_pottemp_non_linear( air_pressure(param_set, T, ρ, q), q, ), - FT(T_min(param_set)), - FT(T_max(param_set)), - SecantMethod(), + SecantMethod(_T_min, _T_max), CompactSolution(), - tol, + SolutionTolerance(tol), maxiter, ) if !sol.converged @@ -1315,42 +1384,46 @@ function air_temperature_from_liquid_ice_pottemp_non_linear( end """ - air_temperature_from_liquid_ice_pottemp_given_pressure(θ_liq_ice, p[, q::PhasePartition]) + air_temperature_from_liquid_ice_pottemp_given_pressure(param_set, θ_liq_ice, p[, q::PhasePartition]) The air temperature where + - `param_set` an `AbstractParameterSet`, see the [`MoistThermodynamics`](@ref) for more details - `θ_liq_ice` liquid-ice potential temperature - `p` pressure and, optionally, - `q` [`PhasePartition`](@ref). Without this argument, the results are for dry air. """ function air_temperature_from_liquid_ice_pottemp_given_pressure( - param_set::PS, + param_set::APS, θ_liq_ice::FT, p::FT, q::PhasePartition{FT} = q_pt_0(FT), -) where {FT <: Real, PS} +) where {FT <: Real} return θ_liq_ice * exner_given_pressure(param_set, p, q) + latent_heat_liq_ice(param_set, q) / cp_m(param_set, q) end """ - virtual_pottemp(T, ρ[, q::PhasePartition]) + virtual_pottemp(param_set, T, ρ[, q::PhasePartition]) The virtual temperature where + - `param_set` an `AbstractParameterSet`, see the [`MoistThermodynamics`](@ref) for more details - `T` temperature - `ρ` (moist-)air density and, optionally, - `q` [`PhasePartition`](@ref). Without this argument, the results are for dry air. """ function virtual_pottemp( - param_set::PS, + param_set::APS, T::FT, ρ::FT, q::PhasePartition{FT} = q_pt_0(FT), -) where {FT <: Real, PS} - return gas_constant_air(param_set, q) / FT(R_d(param_set)) * dry_pottemp(param_set, T, ρ, q) +) where {FT <: Real} + _R_d::FT = R_d(param_set) + return gas_constant_air(param_set, q) / _R_d * + dry_pottemp(param_set, T, ρ, q) end """ @@ -1367,46 +1440,48 @@ virtual_pottemp(ts::ThermodynamicState) = virtual_pottemp( ) """ - liquid_ice_pottemp_sat(T, ρ[, q::PhasePartition]) + liquid_ice_pottemp_sat(param_set, T, ρ[, q::PhasePartition]) The saturated liquid ice potential temperature where + - `param_set` an `AbstractParameterSet`, see the [`MoistThermodynamics`](@ref) for more details - `T` temperature - `ρ` (moist-)air density and, optionally, - `q` [`PhasePartition`](@ref). Without this argument, the results are for dry air. """ function liquid_ice_pottemp_sat( - param_set::PS, + param_set::APS, T::FT, ρ::FT, q::PhasePartition{FT} = q_pt_0(FT), -) where {FT <: Real, PS} +) where {FT <: Real} q_v_sat = q_vap_saturation(param_set, T, ρ, q) return liquid_ice_pottemp(param_set, T, ρ, PhasePartition(q_v_sat)) end """ - liquid_ice_pottemp_sat(T, ρ, q_tot) + liquid_ice_pottemp_sat(param_set, T, ρ, q_tot) The saturated liquid ice potential temperature where + - `param_set` an `AbstractParameterSet`, see the [`MoistThermodynamics`](@ref) for more details - `T` temperature - `ρ` (moist-)air density - `q_tot` total specific humidity """ function liquid_ice_pottemp_sat( - param_set::PS, + param_set::APS, T::FT, ρ::FT, q_tot::FT, -) where {FT <: Real, PS} - return liquid_ice_pottemp( - param_set, - T, - ρ, - PhasePartition_equil(param_set, T, ρ, q_tot), -) +) where {FT <: Real} + return liquid_ice_pottemp( + param_set, + T, + ρ, + PhasePartition_equil(param_set, T, ρ, q_tot), + ) end """ @@ -1422,40 +1497,43 @@ liquid_ice_pottemp_sat(ts::ThermodynamicState) = liquid_ice_pottemp_sat( ) """ - exner_given_pressure(p[, q::PhasePartition]) + exner_given_pressure(param_set, p[, q::PhasePartition]) The Exner function where + - `param_set` an `AbstractParameterSet`, see the [`MoistThermodynamics`](@ref) for more details - `p` pressure and, optionally, - `q` [`PhasePartition`](@ref). Without this argument, the results are for dry air. """ function exner_given_pressure( - param_set::PS, + param_set::APS, p::FT, q::PhasePartition{FT} = q_pt_0(FT), -) where {FT <: Real, PS} +) where {FT <: Real} + _MSLP::FT = MSLP(param_set) # gas constant and isobaric specific heat of moist air _R_m = gas_constant_air(param_set, q) _cp_m = cp_m(param_set, q) - return (p / FT(MSLP(param_set)))^(_R_m / _cp_m) + return (p / _MSLP)^(_R_m / _cp_m) end """ - exner(T, ρ[, q::PhasePartition)]) + exner(param_set, T, ρ[, q::PhasePartition)]) The Exner function where + - `param_set` an `AbstractParameterSet`, see the [`MoistThermodynamics`](@ref) for more details - `T` temperature - `ρ` (moist-)air density and, optionally, - `q` [`PhasePartition`](@ref). Without this argument, the results are for dry air. """ function exner( - param_set::PS, + param_set::APS, T::FT, ρ::FT, q::PhasePartition{FT} = q_pt_0(FT), -) where {FT <: Real, PS} +) where {FT <: Real} return exner_given_pressure(param_set, air_pressure(param_set, T, ρ, q), q) end @@ -1472,27 +1550,28 @@ exner(ts::ThermodynamicState) = exner( ) """ - relative_humidity(T, p, e_int, q::PhasePartition) + relative_humidity(param_set, T, p, e_int, q::PhasePartition) The relative humidity, given - - `T` temperature + - `param_set` an `AbstractParameterSet`, see the [`MoistThermodynamics`](@ref) for more details - `p` pressure - `e_int` internal energy per unit mass and, optionally, - `q` [`PhasePartition`](@ref). Without this argument, the results are for dry air. """ function relative_humidity( - param_set::PS, + param_set::APS, T::FT, p::FT, e_int::FT, q::PhasePartition{FT} = q_pt_0(FT), -) where {FT <: Real, PS} +) where {FT <: Real} + _R_v::FT = R_v(param_set) q_vap = q.tot - q.liq - q.ice p_vap = q_vap * air_density(param_set, T, p, q) * - FT(R_v(param_set)) * + _R_v * air_temperature(param_set, e_int, q) liq_frac = liquid_fraction(param_set, T, q) p_vap_sat = @@ -1514,4 +1593,3 @@ relative_humidity(ts::ThermodynamicState{FT}) where {FT <: Real} = internal_energy(ts), PhasePartition(ts), ) - diff --git a/src/states.jl b/src/states.jl index 90adc40e..f54b49cb 100644 --- a/src/states.jl +++ b/src/states.jl @@ -60,14 +60,14 @@ may be needed). # Constructors - PhaseEquil(e_int, ρ, q_tot) + PhaseEquil(param_set, e_int, ρ, q_tot) # Fields $(DocStringExtensions.FIELDS) """ struct PhaseEquil{FT, PS} <: ThermodynamicState{FT} - "parameter set (e.g., planet parameters)" + "parameter set, used to dispatch planet parameter function calls" param_set::PS "internal energy" e_int::FT @@ -79,19 +79,19 @@ struct PhaseEquil{FT, PS} <: ThermodynamicState{FT} T::FT end function PhaseEquil( - param_set::PS, + param_set::APS, e_int::FT, ρ::FT, q_tot::FT, maxiter::Int = 3, tol::FT = FT(1e-1), sat_adjust::Function = saturation_adjustment, -) where {FT <: Real, PS} +) where {FT <: Real} # TODO: Remove these safety nets, or at least add warnings # waiting on fix: github.com/vchuravy/GPUifyLoops.jl/issues/104 q_tot_safe = clamp(q_tot, FT(0), FT(1)) T = sat_adjust(param_set, e_int, ρ, q_tot_safe, maxiter, tol) - return PhaseEquil{FT, PS}(param_set, e_int, ρ, q_tot_safe, T) + return PhaseEquil{FT, typeof(param_set)}(param_set, e_int, ρ, q_tot_safe, T) end """ @@ -101,45 +101,45 @@ A dry thermodynamic state (`q_tot = 0`). # Constructors - PhaseDry(e_int, ρ) + PhaseDry(param_set, e_int, ρ) # Fields $(DocStringExtensions.FIELDS) """ struct PhaseDry{FT, PS} <: ThermodynamicState{FT} - "parameter set (e.g., planet parameters)" + "parameter set, used to dispatch planet parameter function calls" param_set::PS "internal energy" e_int::FT "density of dry air" ρ::FT end +PhaseDry(param_set::APS, e_int::FT, ρ::FT) where {FT} = + PhaseDry{FT, typeof(param_set)}(param_set, e_int, ρ) """ - PhaseDry_given_pT(p, T) + PhaseDry_given_pT(param_set, p, T) Constructs a [`PhaseDry`](@ref) thermodynamic state from: + - `param_set` an `AbstractParameterSet`, see the [`MoistThermodynamics`](@ref) for more details - `p` pressure - `T` temperature """ -function PhaseDry_given_pT( - param_set::PS, - p::FT, - T::FT, -) where {FT <: Real, PS} +function PhaseDry_given_pT(param_set::APS, p::FT, T::FT) where {FT <: Real} e_int = internal_energy(param_set, T) ρ = air_density(param_set, T, p) - return PhaseDry{FT, PS}(param_set, e_int, ρ) + return PhaseDry{FT, typeof(param_set)}(param_set, e_int, ρ) end """ - LiquidIcePotTempSHumEquil(θ_liq_ice, ρ, q_tot) + LiquidIcePotTempSHumEquil(param_set, θ_liq_ice, ρ, q_tot) Constructs a [`PhaseEquil`](@ref) thermodynamic state from: + - `param_set` an `AbstractParameterSet`, see the [`MoistThermodynamics`](@ref) for more details - `θ_liq_ice` liquid-ice potential temperature - `ρ` (moist-)air density - `q_tot` total specific humidity @@ -147,13 +147,13 @@ Constructs a [`PhaseEquil`](@ref) thermodynamic state from: - `maxiter` maximum iterations for saturation adjustment """ function LiquidIcePotTempSHumEquil( - param_set::PS, + param_set::APS, θ_liq_ice::FT, ρ::FT, q_tot::FT, maxiter::Int = 30, tol::FT = FT(1e-1), -) where {FT <: Real, PS} +) where {FT <: Real} T = saturation_adjustment_q_tot_θ_liq_ice( param_set, θ_liq_ice, @@ -164,15 +164,15 @@ function LiquidIcePotTempSHumEquil( ) q_pt = PhasePartition_equil(param_set, T, ρ, q_tot) e_int = internal_energy(param_set, T, q_pt) - return PhaseEquil{FT, PS}(param_set, e_int, ρ, q_tot, T) + return PhaseEquil{FT, typeof(param_set)}(param_set, e_int, ρ, q_tot, T) end - """ - LiquidIcePotTempSHumEquil_given_pressure(θ_liq_ice, p, q_tot) + LiquidIcePotTempSHumEquil_given_pressure(param_set, θ_liq_ice, p, q_tot) Constructs a [`PhaseEquil`](@ref) thermodynamic state from: + - `param_set` an `AbstractParameterSet`, see the [`MoistThermodynamics`](@ref) for more details - `θ_liq_ice` liquid-ice potential temperature - `p` pressure - `q_tot` total specific humidity @@ -180,13 +180,13 @@ Constructs a [`PhaseEquil`](@ref) thermodynamic state from: - `maxiter` maximum iterations for saturation adjustment """ function LiquidIcePotTempSHumEquil_given_pressure( - param_set::PS, + param_set::APS, θ_liq_ice::FT, p::FT, q_tot::FT, maxiter::Int = 30, tol::FT = FT(1e-1), -) where {FT <: Real, PS} +) where {FT <: Real} T = saturation_adjustment_q_tot_θ_liq_ice_given_pressure( param_set, θ_liq_ice, @@ -198,28 +198,29 @@ function LiquidIcePotTempSHumEquil_given_pressure( ρ = air_density(param_set, T, p, PhasePartition(q_tot)) q = PhasePartition_equil(param_set, T, ρ, q_tot) e_int = internal_energy(param_set, T, q) - return PhaseEquil{FT, PS}(param_set, e_int, ρ, q.tot, T) + return PhaseEquil{FT, typeof(param_set)}(param_set, e_int, ρ, q.tot, T) end """ - TemperatureSHumEquil(T, p, q_tot) + TemperatureSHumEquil(param_set, T, p, q_tot) Constructs a [`PhaseEquil`](@ref) thermodynamic state from temperature. + - `param_set` an `AbstractParameterSet`, see the [`MoistThermodynamics`](@ref) for more details - `T` temperature - `p` pressure - `q_tot` total specific humidity """ function TemperatureSHumEquil( - param_set::PS, + param_set::APS, T::FT, p::FT, q_tot::FT, -) where {FT <: Real, PS} +) where {FT <: Real} ρ = air_density(param_set, T, p, PhasePartition(q_tot)) q = PhasePartition_equil(param_set, T, ρ, q_tot) e_int = internal_energy(param_set, T, q) - return PhaseEquil{FT, PS}(param_set, e_int, ρ, q_tot, T) + return PhaseEquil{FT, typeof(param_set)}(param_set, e_int, ρ, q_tot, T) end """ @@ -230,7 +231,7 @@ be computed directly). # Constructors - PhaseNonEquil(e_int, q::PhasePartition, ρ) + PhaseNonEquil(param_set, e_int, q::PhasePartition, ρ) # Fields @@ -238,7 +239,7 @@ $(DocStringExtensions.FIELDS) """ struct PhaseNonEquil{FT, PS} <: ThermodynamicState{FT} - "parameter set (e.g., planet parameters)" + "parameter set, used to dispatch planet parameter function calls" param_set::PS "internal energy" e_int::FT @@ -248,18 +249,20 @@ struct PhaseNonEquil{FT, PS} <: ThermodynamicState{FT} q::PhasePartition{FT} end function PhaseNonEquil( - param_set::PS, + param_set::APS, e_int::FT, ρ::FT, -) where {FT, PS} - return PhaseNonEquil{FT, PS}(param_set, e_int, ρ, q_pt_0(FT)) + q::PhasePartition{FT} = q_pt_0(FT), +) where {FT} + return PhaseNonEquil{FT, typeof(param_set)}(param_set, e_int, ρ, q) end """ - LiquidIcePotTempSHumNonEquil(θ_liq_ice, ρ, q_pt) + LiquidIcePotTempSHumNonEquil(param_set, θ_liq_ice, ρ, q_pt) Constructs a [`PhaseNonEquil`](@ref) thermodynamic state from: + - `param_set` an `AbstractParameterSet`, see the [`MoistThermodynamics`](@ref) for more details - `θ_liq_ice` liquid-ice potential temperature - `ρ` (moist-)air density - `q_pt` phase partition @@ -268,13 +271,13 @@ and, optionally - `maxiter` maximum iterations for non-linear equation solve """ function LiquidIcePotTempSHumNonEquil( - param_set::PS, + param_set::APS, θ_liq_ice::FT, ρ::FT, q_pt::PhasePartition{FT}, maxiter::Int = 5, tol::FT = FT(1e-1), -) where {FT <: Real, PS} +) where {FT <: Real} T = air_temperature_from_liquid_ice_pottemp_non_linear( param_set, θ_liq_ice, @@ -284,24 +287,25 @@ function LiquidIcePotTempSHumNonEquil( q_pt, ) e_int = internal_energy(param_set, T, q_pt) - return PhaseNonEquil{FT, PS}(param_set, e_int, ρ, q_pt) + return PhaseNonEquil{FT, typeof(param_set)}(param_set, e_int, ρ, q_pt) end """ - LiquidIcePotTempSHumNonEquil_given_pressure(θ_liq_ice, p, q_pt) + LiquidIcePotTempSHumNonEquil_given_pressure(param_set, θ_liq_ice, p, q_pt) Constructs a [`PhaseNonEquil`](@ref) thermodynamic state from: + - `param_set` an `AbstractParameterSet`, see the [`MoistThermodynamics`](@ref) for more details - `θ_liq_ice` liquid-ice potential temperature - `p` pressure - `q_pt` phase partition """ function LiquidIcePotTempSHumNonEquil_given_pressure( - param_set::PS, + param_set::APS, θ_liq_ice::FT, p::FT, q_pt::PhasePartition{FT}, -) where {FT <: Real, PS} +) where {FT <: Real} T = air_temperature_from_liquid_ice_pottemp_given_pressure( param_set, θ_liq_ice, @@ -310,81 +314,5 @@ function LiquidIcePotTempSHumNonEquil_given_pressure( ) ρ = air_density(param_set, T, p, q_pt) e_int = internal_energy(param_set, T, q_pt) - return PhaseNonEquil{FT, PS}(param_set, e_int, ρ, q_pt) -end - -""" - fixed_lapse_rate_ref_state(z::FT, - T_surface::FT, - T_min::FT) where {FT<:AbstractFloat} - -Fixed lapse rate hydrostatic reference state -""" -function fixed_lapse_rate_ref_state( - param_set::PS, - z::FT, - T_surface::FT, - T_min::FT, -) where {FT <: AbstractFloat, PS} - Γ = FT(grav(param_set)) / FT(cp_d(param_set)) - z_tropopause = (T_surface - T_min) / Γ - H_min = FT(R_d(param_set)) * T_min / FT(grav(param_set)) - T = max(T_surface - Γ * z, T_min) - p = FT(MSLP(param_set)) * (T / T_surface)^(FT(grav(param_set)) / (FT(R_d(param_set)) * Γ)) - T == T_min && (p = p * exp(-(z - z_tropopause) / FT(H_min))) - ρ = p / (FT(R_d(param_set)) * T) - return T, p, ρ -end - -""" - tested_convergence_range(FT, n) - -A range of input arguments to thermodynamic state constructors - - `e_int` internal energy - - `ρ` (moist-)air density - - `q_tot` total specific humidity - - `q_pt` phase partition - - `T` air temperature - - `θ_liq_ice` liquid-ice potential temperature -that are tested for convergence in saturation adjustment. - -Note that the output vectors are of size ``n*n_RH``, and they -should span the input arguments to all of the constructors. -""" -function tested_convergence_range( - param_set::PS, - n::Int, - ::Type{FT}, -) where {FT, PS} - n_RS1 = 10 - n_RS2 = 20 - n_RS = n_RS1 + n_RS2 - z_range = range(FT(0), stop = FT(2.5e4), length = n) - relative_sat1 = range(FT(0), stop = FT(1), length = n_RS1) - relative_sat2 = range(FT(1), stop = FT(1.02), length = n_RS2) - relative_sat = [relative_sat1..., relative_sat2...] - T_min = FT(150) - T_surface = FT(350) - - args = - fixed_lapse_rate_ref_state.( - Ref(param_set), - z_range, - Ref(T_surface), - Ref(T_min), - ) - T, p, ρ = getindex.(args, 1), getindex.(args, 2), getindex.(args, 3) - - p = collect(Iterators.flatten([p for RS in 1:n_RS])) - ρ = collect(Iterators.flatten([ρ for RS in 1:n_RS])) - T = collect(Iterators.flatten([T for RS in 1:n_RS])) - relative_sat = collect(Iterators.flatten([relative_sat for RS in 1:n])) - - # Additional variables - q_sat = q_vap_saturation.(Ref(param_set), T, ρ) - q_tot = min.(relative_sat .* q_sat, FT(1)) - q_pt = PhasePartition_equil.(Ref(param_set), T, ρ, q_tot) - e_int = internal_energy.(Ref(param_set), T, q_pt) - θ_liq_ice = liquid_ice_pottemp.(Ref(param_set), T, ρ, q_pt) - return e_int, ρ, q_tot, q_pt, T, p, θ_liq_ice + return PhaseNonEquil{FT, typeof(param_set)}(param_set, e_int, ρ, q_pt) end diff --git a/test/data_tests.jl b/test/data_tests.jl index b992a660..f92f4d67 100644 --- a/test/data_tests.jl +++ b/test/data_tests.jl @@ -4,24 +4,26 @@ include("ArtifactWrappers.jl") using .ArtifactWrappers # Get dycoms dataset folder: -dycoms_dataset = - ArtifactWrapper(joinpath(@__DIR__, "Artifacts.toml"), "dycoms", - ArtifactFile[ - ArtifactFile(url="https://caltech.box.com/shared/static/bxau6i46y6ikxn2sy9krgz0sw5vuptfo.nc", - filename="test_data_PhaseEquil.nc"), - ] - ) +dycoms_dataset = ArtifactWrapper( + joinpath(@__DIR__, "Artifacts.toml"), + "dycoms", + ArtifactFile[ArtifactFile( + url = "https://caltech.box.com/shared/static/bxau6i46y6ikxn2sy9krgz0sw5vuptfo.nc", + filename = "test_data_PhaseEquil.nc", + ),], +) dycoms_dataset_path = get_data_folder(dycoms_dataset) @testset "Data tests" begin FT = Float64 - e_int, ρ, q_tot, q_pt, T, p, θ_liq_ice = MT.tested_convergence_range(param_set, 50, FT) + data = joinpath(dycoms_dataset_path, "test_data_PhaseEquil.nc") ds_PhaseEquil = Dataset(data, "r") e_int = Array{FT}(ds_PhaseEquil["e_int"][:]) ρ = Array{FT}(ds_PhaseEquil["ρ"][:]) q_tot = Array{FT}(ds_PhaseEquil["q_tot"][:]) - # ts = PhaseEquil.(Ref(param_set), e_int, ρ, q_tot) # Fails + ts = PhaseEquil.(Ref(param_set), e_int, ρ, q_tot, 4) + # ts = PhaseEquil.(Ref(param_set), e_int, ρ, q_tot, 3) # Fails end diff --git a/test/runtests.jl b/test/runtests.jl index 1361b110..97a9eb8f 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,146 +1,245 @@ +push!(LOAD_PATH, joinpath(@__DIR__, "..", "env", "test")) + +rm(joinpath(@__DIR__, "..", "Manifest.toml"), force=true) +using Pkg +Pkg.activate(joinpath(@__DIR__, "..")) +Pkg.instantiate(; verbose = true) + using Test using MoistThermodynamics +using NCDatasets +using Random MT = MoistThermodynamics + using CLIMAParameters using CLIMAParameters.Planet -using NCDatasets -using Random struct EarthParameterSet <: AbstractEarthParameterSet end -param_set = EarthParameterSet() +const param_set = EarthParameterSet() float_types = [Float32, Float64] include("data_tests.jl") -@testset "Isentropic processes" begin +@testset "moist thermodynamics - isentropic processes" begin for FT in [Float64] + _R_d = FT(R_d(param_set)) + _molmass_ratio = FT(molmass_ratio(param_set)) + _cp_d = FT(cp_d(param_set)) + _cp_v = FT(cp_v(param_set)) + _cp_l = FT(cp_l(param_set)) + _cv_d = FT(cv_d(param_set)) + _cv_v = FT(cv_v(param_set)) + _cv_l = FT(cv_l(param_set)) + _cv_i = FT(cv_i(param_set)) + _T_0 = FT(T_0(param_set)) + _e_int_v0 = FT(e_int_v0(param_set)) + _e_int_i0 = FT(e_int_i0(param_set)) + _LH_v0 = FT(LH_v0(param_set)) + _LH_s0 = FT(LH_s0(param_set)) + _cp_i = FT(cp_i(param_set)) + _LH_f0 = FT(LH_f0(param_set)) + _press_triple = FT(press_triple(param_set)) + _R_v = FT(R_v(param_set)) + _T_triple = FT(T_triple(param_set)) + _T_freeze = FT(T_freeze(param_set)) + _T_min = FT(T_min(param_set)) + _MSLP = FT(MSLP(param_set)) + _T_max = FT(T_max(param_set)) + _kappa_d = FT(kappa_d(param_set)) + # for FT in float_types e_int, ρ, q_tot, q_pt, T, p, θ_liq_ice = - MT.tested_convergence_range(param_set, 50, FT) + MT.tested_profiles(param_set, 50, FT) Φ = FT(1) Random.seed!(15) perturbation = FT(0.1) * rand(length(T)) # TODO: Use reasonable values for ambient temperature/pressure T∞, p∞ = T .* perturbation, p .* perturbation - - T_air = (p ./ FT(MSLP(param_set))) .^ (FT(R_d(param_set)) / FT(cp_d(param_set))) .* θ_liq_ice - @test air_temperature.(Ref(param_set), p, θ_liq_ice, Ref(DryAdiabaticProcess())) ≈ T_air - - p_air = FT(MSLP(param_set)) .* (1 .- Φ ./ (θ_liq_ice .* FT(cp_d(param_set)))) .^ (FT(cp_d(param_set)) / FT(R_d(param_set))) - @test air_pressure_given_θ.(Ref(param_set), θ_liq_ice, Φ, Ref(DryAdiabaticProcess())) ≈ p_air - - p_air = p∞ .* (T ./ T∞) .^ (FT(1) / FT(kappa_d(param_set))) - @test air_pressure.(Ref(param_set), T, T∞, p∞, Ref(DryAdiabaticProcess())) ≈ p_air + @test air_temperature.( + Ref(param_set), + p, + θ_liq_ice, + Ref(DryAdiabaticProcess()), + ) ≈ (p ./ _MSLP) .^ (_R_d / _cp_d) .* θ_liq_ice + @test air_pressure_given_θ.( + Ref(param_set), + θ_liq_ice, + Φ, + Ref(DryAdiabaticProcess()), + ) ≈ _MSLP .* (1 .- Φ ./ (θ_liq_ice .* _cp_d)) .^ (_cp_d / _R_d) + @test air_pressure.( + Ref(param_set), + T, + T∞, + p∞, + Ref(DryAdiabaticProcess()), + ) ≈ p∞ .* (T ./ T∞) .^ (FT(1) / _kappa_d) end end -@testset "Correctness" begin +@testset "moist thermodynamics - correctness" begin FT = Float64 + _R_d = FT(R_d(param_set)) + _molmass_ratio = FT(molmass_ratio(param_set)) + _cp_d = FT(cp_d(param_set)) + _cp_v = FT(cp_v(param_set)) + _cp_l = FT(cp_l(param_set)) + _cv_d = FT(cv_d(param_set)) + _cv_v = FT(cv_v(param_set)) + _cv_l = FT(cv_l(param_set)) + _cv_i = FT(cv_i(param_set)) + _T_0 = FT(T_0(param_set)) + _e_int_v0 = FT(e_int_v0(param_set)) + _e_int_i0 = FT(e_int_i0(param_set)) + _LH_v0 = FT(LH_v0(param_set)) + _LH_s0 = FT(LH_s0(param_set)) + _cp_i = FT(cp_i(param_set)) + _LH_f0 = FT(LH_f0(param_set)) + _press_triple = FT(press_triple(param_set)) + _R_v = FT(R_v(param_set)) + _T_triple = FT(T_triple(param_set)) + _T_freeze = FT(T_freeze(param_set)) + _T_min = FT(T_min(param_set)) + _MSLP = FT(MSLP(param_set)) + _T_max = FT(T_max(param_set)) + _kappa_d = FT(kappa_d(param_set)) + _T_icenuc = FT(T_icenuc(param_set)) + # ideal gas law - @test air_pressure(param_set, FT(1), FT(1), PhasePartition(FT(1))) === FT(R_v(param_set)) - @test air_pressure(param_set, FT(1), FT(2), PhasePartition(FT(1), FT(0.5), FT(0))) === - FT(R_v(param_set)) - @test air_pressure(param_set, FT(1), FT(1)) === FT(R_d(param_set)) - @test air_pressure(param_set, FT(1), FT(2)) === 2 * FT(R_d(param_set)) - @test air_density(param_set, FT(1), FT(1)) === 1 / FT(R_d(param_set)) - @test air_density(param_set, FT(1), FT(2)) === 2 / FT(R_d(param_set)) + @test air_pressure(param_set, FT(1), FT(1), PhasePartition(FT(1))) === _R_v + @test air_pressure( + param_set, + FT(1), + FT(2), + PhasePartition(FT(1), FT(0.5), FT(0)), + ) === _R_v + @test air_pressure(param_set, FT(1), FT(1)) === _R_d + @test air_pressure(param_set, FT(1), FT(2)) === 2 * _R_d + @test air_density(param_set, FT(1), FT(1)) === 1 / _R_d + @test air_density(param_set, FT(1), FT(2)) === 2 / _R_d # gas constants and heat capacities - @test gas_constant_air(param_set, PhasePartition(FT(0))) === FT(R_d(param_set)) - @test gas_constant_air(param_set, PhasePartition(FT(1))) === FT(R_v(param_set)) - @test gas_constant_air(param_set, PhasePartition(FT(0.5), FT(0.5))) ≈ FT(R_d(param_set)) / 2 - @test gas_constant_air(param_set, FT) == FT(R_d(param_set)) - - @test cp_m(param_set, PhasePartition(FT(0))) === FT(cp_d(param_set)) - @test cp_m(param_set, PhasePartition(FT(1))) === FT(cp_v(param_set)) - @test cp_m(param_set, PhasePartition(FT(1), FT(1))) === FT(cp_l(param_set)) - @test cp_m(param_set, PhasePartition(FT(1), FT(0), FT(1))) === FT(cp_i(param_set)) - @test cp_m(param_set, FT) == FT(cp_d(param_set)) - - @test cv_m(param_set, PhasePartition(FT(0))) === FT(cp_d(param_set) - R_d(param_set)) - @test cv_m(param_set, PhasePartition(FT(1))) === FT(cp_v(param_set) - R_v(param_set)) - @test cv_m(param_set, PhasePartition(FT(1), FT(1))) === FT(cv_l(param_set)) - @test cv_m(param_set, PhasePartition(FT(1), FT(0), FT(1))) === FT(cv_i(param_set)) - @test cv_m(param_set, FT) == FT(cv_d(param_set)) + @test gas_constant_air(param_set, PhasePartition(FT(0))) === _R_d + @test gas_constant_air(param_set, PhasePartition(FT(1))) === _R_v + @test gas_constant_air(param_set, PhasePartition(FT(0.5), FT(0.5))) ≈ + _R_d / 2 + @test gas_constant_air(param_set, FT) == _R_d + + @test cp_m(param_set, PhasePartition(FT(0))) === _cp_d + @test cp_m(param_set, PhasePartition(FT(1))) === _cp_v + @test cp_m(param_set, PhasePartition(FT(1), FT(1))) === _cp_l + @test cp_m(param_set, PhasePartition(FT(1), FT(0), FT(1))) === _cp_i + @test cp_m(param_set, FT) == _cp_d + + @test cv_m(param_set, PhasePartition(FT(0))) === _cp_d - _R_d + @test cv_m(param_set, PhasePartition(FT(1))) === _cp_v - _R_v + @test cv_m(param_set, PhasePartition(FT(1), FT(1))) === _cv_l + @test cv_m(param_set, PhasePartition(FT(1), FT(0), FT(1))) === _cv_i + @test cv_m(param_set, FT) == _cv_d # speed of sound - soundspeedair = sqrt(cp_d(param_set) / cv_d(param_set) * R_d(param_set) * (T_0(param_set) + 20)) - @test soundspeed_air(param_set, T_0(param_set) + 20, PhasePartition(FT(0))) == soundspeedair - soundspeedair = sqrt(cp_v(param_set) / cv_v(param_set) * R_v(param_set) * (T_0(param_set) + 100)) - @test soundspeed_air(param_set, T_0(param_set) + 100, PhasePartition(FT(1))) == soundspeedair + @test soundspeed_air(param_set, _T_0 + 20, PhasePartition(FT(0))) == + sqrt(_cp_d / _cv_d * _R_d * (_T_0 + 20)) + @test soundspeed_air(param_set, _T_0 + 100, PhasePartition(FT(1))) == + sqrt(_cp_v / _cv_v * _R_v * (_T_0 + 100)) # specific latent heats - @test latent_heat_vapor(param_set, FT(T_0(param_set))) ≈ LH_v0(param_set) - @test latent_heat_fusion(param_set, FT(T_0(param_set))) ≈ LH_f0(param_set) - @test latent_heat_sublim(param_set, FT(T_0(param_set))) ≈ LH_s0(param_set) + @test latent_heat_vapor(param_set, _T_0) ≈ _LH_v0 + @test latent_heat_fusion(param_set, _T_0) ≈ _LH_f0 + @test latent_heat_sublim(param_set, _T_0) ≈ _LH_s0 # saturation vapor pressure and specific humidity p = FT(1.e5) q_tot = FT(0.23) ρ = FT(1.0) - ρ_v_triple = press_triple(param_set) / R_v(param_set) / T_triple(param_set) - @test saturation_vapor_pressure(param_set, FT(T_triple(param_set)), Liquid()) ≈ press_triple(param_set) - @test saturation_vapor_pressure(param_set, FT(T_triple(param_set)), Ice()) ≈ press_triple(param_set) + ρ_v_triple = _press_triple / _R_v / _T_triple + @test saturation_vapor_pressure(param_set, _T_triple, Liquid()) ≈ + _press_triple + @test saturation_vapor_pressure(param_set, _T_triple, Ice()) ≈ _press_triple - @test q_vap_saturation(param_set, FT(T_triple(param_set)), ρ, PhasePartition(FT(0))) == - ρ_v_triple / ρ - @test q_vap_saturation(param_set, FT(T_triple(param_set)), ρ, PhasePartition(q_tot, q_tot)) == + @test q_vap_saturation(param_set, _T_triple, ρ, PhasePartition(FT(0))) == ρ_v_triple / ρ + @test q_vap_saturation( + param_set, + _T_triple, + ρ, + PhasePartition(q_tot, q_tot), + ) == ρ_v_triple / ρ - @test q_vap_saturation_generic(param_set, FT(T_triple(param_set)), ρ; phase = Liquid()) == + @test q_vap_saturation_generic(param_set, _T_triple, ρ; phase = Liquid()) == ρ_v_triple / ρ - @test q_vap_saturation_generic(param_set, FT(T_triple(param_set)), ρ; phase = Ice()) == + @test q_vap_saturation_generic(param_set, _T_triple, ρ; phase = Ice()) == ρ_v_triple / ρ - @test q_vap_saturation_generic(param_set, FT(T_triple(param_set) - 20), ρ; phase = Liquid()) >= - q_vap_saturation_generic(param_set, FT(T_triple(param_set) - 20), ρ; phase = Ice()) - - @test saturation_excess(param_set, FT(T_triple(param_set)), ρ, PhasePartition(q_tot)) == + @test q_vap_saturation_generic( + param_set, + _T_triple - 20, + ρ; + phase = Liquid(), + ) >= q_vap_saturation_generic(param_set, _T_triple - 20, ρ; phase = Ice()) + + @test saturation_excess(param_set, _T_triple, ρ, PhasePartition(q_tot)) == q_tot - ρ_v_triple / ρ - @test saturation_excess(param_set, FT(T_triple(param_set)), ρ, PhasePartition(q_tot / 1000)) == - 0.0 + @test saturation_excess( + param_set, + _T_triple, + ρ, + PhasePartition(q_tot / 1000), + ) == 0.0 # energy functions and inverse (temperature) T = FT(300) e_kin = FT(11) e_pot = FT(13) - @test air_temperature(param_set, FT(cv_d(param_set) * (T - T_0(param_set)))) === FT(T) - @test air_temperature(param_set, FT(cv_d(param_set) * (T - T_0(param_set))), PhasePartition(FT(0))) === FT(T) + @test air_temperature(param_set, _cv_d * (T - _T_0)) === FT(T) + @test air_temperature( + param_set, + _cv_d * (T - _T_0), + PhasePartition(FT(0)), + ) === FT(T) - @test air_temperature(param_set, - cv_m(param_set, PhasePartition(FT(0))) * (T - T_0(param_set)), + @test air_temperature( + param_set, + cv_m(param_set, PhasePartition(FT(0))) * (T - _T_0), PhasePartition(FT(0)), ) === FT(T) - @test air_temperature(param_set, - cv_m(param_set, PhasePartition(FT(q_tot))) * (T - T_0(param_set)) + q_tot * e_int_v0(param_set), + @test air_temperature( + param_set, + cv_m(param_set, PhasePartition(FT(q_tot))) * (T - _T_0) + + q_tot * _e_int_v0, PhasePartition(q_tot), ) ≈ FT(T) - @test total_energy(param_set, FT(e_kin), FT(e_pot), FT(T_0(param_set))) === FT(e_kin) + FT(e_pot) + @test total_energy(param_set, FT(e_kin), FT(e_pot), _T_0) === + FT(e_kin) + FT(e_pot) @test total_energy(param_set, FT(e_kin), FT(e_pot), FT(T)) ≈ - FT(e_kin) + FT(e_pot) + cv_d(param_set) * (T - T_0(param_set)) - @test total_energy(param_set, FT(0), FT(0), FT(T_0(param_set)), PhasePartition(q_tot)) ≈ - q_tot * e_int_v0(param_set) + FT(e_kin) + FT(e_pot) + _cv_d * (T - _T_0) + @test total_energy(param_set, FT(0), FT(0), _T_0, PhasePartition(q_tot)) ≈ + q_tot * _e_int_v0 # phase partitioning in equilibrium q_liq = FT(0.1) - T = FT(T_icenuc(param_set) - 10) + T = FT(_T_icenuc - 10) ρ = FT(1.0) q_tot = FT(0.21) @test liquid_fraction(param_set, T) === FT(0) - @test liquid_fraction(param_set, T, PhasePartition(q_tot, q_liq, q_liq)) === FT(0.5) + @test liquid_fraction(param_set, T, PhasePartition(q_tot, q_liq, q_liq)) === + FT(0.5) q = PhasePartition_equil(param_set, T, ρ, q_tot) @test q.liq ≈ FT(0) @test 0 < q.ice <= q_tot - T = FT(T_freeze(param_set) + 10) + T = FT(_T_freeze + 10) ρ = FT(0.1) q_tot = FT(0.60) @test liquid_fraction(param_set, T) === FT(1) - @test liquid_fraction(param_set, T, PhasePartition(q_tot, q_liq, q_liq / 2)) === - FT(2 / 3) + @test liquid_fraction( + param_set, + T, + PhasePartition(q_tot, q_liq, q_liq / 2), + ) === FT(2 / 3) q = PhasePartition_equil(param_set, T, ρ, q_tot) @test 0 < q.liq <= q_tot @test q.ice ≈ 0 @@ -150,7 +249,8 @@ end tol_T = 1e-1 q_tot = FT(0) ρ = FT(1) - @test MT.saturation_adjustment_SecantMethod(param_set, + @test MT.saturation_adjustment_SecantMethod( + param_set, internal_energy_sat(param_set, 300.0, ρ, q_tot), ρ, q_tot, @@ -158,7 +258,8 @@ end 1e-2, ) ≈ 300.0 @test abs( - MT.saturation_adjustment(param_set, + MT.saturation_adjustment( + param_set, internal_energy_sat(param_set, 300.0, ρ, q_tot), ρ, q_tot, @@ -169,7 +270,8 @@ end q_tot = FT(0.21) ρ = FT(0.1) - @test MT.saturation_adjustment_SecantMethod(param_set, + @test MT.saturation_adjustment_SecantMethod( + param_set, internal_energy_sat(param_set, 200.0, ρ, q_tot), ρ, q_tot, @@ -177,7 +279,8 @@ end 1e-2, ) ≈ 200.0 @test abs( - MT.saturation_adjustment(param_set, + MT.saturation_adjustment( + param_set, internal_energy_sat(param_set, 200.0, ρ, q_tot), ρ, q_tot, @@ -187,8 +290,8 @@ end ) < tol_T q = PhasePartition_equil(param_set, T, ρ, q_tot) @test q.tot - q.liq - q.ice ≈ - vapor_specific_humidity(q) ≈ - q_vap_saturation(param_set, T, ρ) + vapor_specific_humidity(q) ≈ + q_vap_saturation(param_set, T, ρ) ρ = FT(1) ρu = FT[1, 2, 3] @@ -198,21 +301,24 @@ end # potential temperatures T = FT(300) - @test liquid_ice_pottemp_given_pressure(param_set, T, FT(MSLP(param_set))) === T - @test liquid_ice_pottemp_given_pressure(param_set, T, FT(MSLP(param_set)) / 10) ≈ - T * 10^(R_d(param_set) / cp_d(param_set)) - @test liquid_ice_pottemp_given_pressure(param_set, + @test liquid_ice_pottemp_given_pressure(param_set, T, _MSLP) === T + @test liquid_ice_pottemp_given_pressure(param_set, T, _MSLP / 10) ≈ + T * 10^(_R_d / _cp_d) + @test liquid_ice_pottemp_given_pressure( + param_set, T, - FT(MSLP(param_set)) / 10, + _MSLP / 10, PhasePartition(FT(1)), - ) ≈ T * 10^(R_v(param_set) / cp_v(param_set)) + ) ≈ T * 10^(_R_v / _cp_v) # dry potential temperatures. FIXME: add correctness tests T = FT(300) p = FT(1.e5) q_tot = FT(0.23) - @test dry_pottemp_given_pressure(param_set, T, p, PhasePartition(q_tot)) isa typeof(p) - @test air_temperature_from_liquid_ice_pottemp_given_pressure(param_set, + @test dry_pottemp_given_pressure(param_set, T, p, PhasePartition(q_tot)) isa + typeof(p) + @test air_temperature_from_liquid_ice_pottemp_given_pressure( + param_set, dry_pottemp_given_pressure(param_set, T, p, PhasePartition(q_tot)), p, PhasePartition(q_tot), @@ -221,18 +327,19 @@ end # Exner function. FIXME: add correctness tests p = FT(1.e5) q_tot = FT(0.23) - @test exner_given_pressure(param_set, p, PhasePartition(q_tot)) isa typeof(p) + @test exner_given_pressure(param_set, p, PhasePartition(q_tot)) isa + typeof(p) end -@testset "Default behavior accuracy" begin +@testset "moist thermodynamics - default behavior accuracy" begin # Input arguments should be accurate within machine precision # Temperature is approximated via saturation adjustment, and should be within a physical tolerance for FT in float_types rtol = FT(1e-2) e_int, ρ, q_tot, q_pt, T, p, θ_liq_ice = - MT.tested_convergence_range(param_set, 50, FT) + MT.tested_profiles(param_set, 50, FT) # PhaseEquil ts_exact = PhaseEquil.(Ref(param_set), e_int, ρ, q_tot, 100, FT(1e-4)) @@ -297,7 +404,15 @@ end )) # LiquidIcePotTempSHumEquil - ts_exact = LiquidIcePotTempSHumEquil.(Ref(param_set), θ_liq_ice, ρ, q_tot, 40, FT(1e-3)) + ts_exact = + LiquidIcePotTempSHumEquil.( + Ref(param_set), + θ_liq_ice, + ρ, + q_tot, + 40, + FT(1e-3), + ) ts = LiquidIcePotTempSHumEquil.(Ref(param_set), θ_liq_ice, ρ, q_tot) # Should be machine accurate: @test all(air_density.(ts) .≈ air_density.(ts_exact)) @@ -332,11 +447,12 @@ end 40, FT(1e-3), ) - ts = LiquidIcePotTempSHumEquil_given_pressure.( - Ref(param_set), - θ_liq_ice, - p, - q_tot, + ts = + LiquidIcePotTempSHumEquil_given_pressure.( + Ref(param_set), + θ_liq_ice, + p, + q_tot, ) # Should be machine accurate: @test all( @@ -374,13 +490,8 @@ end q_pt, 40, FT(1e-3), - ) - ts = LiquidIcePotTempSHumNonEquil.( - Ref(param_set), - θ_liq_ice, - ρ, - q_pt, ) + ts = LiquidIcePotTempSHumNonEquil.(Ref(param_set), θ_liq_ice, ρ, q_pt) # Should be machine accurate: @test all( getproperty.(PhasePartition.(ts), :tot) .≈ @@ -416,14 +527,17 @@ end end -@testset "Constructor consistency" begin +@testset "moist thermodynamics - constructor consistency" begin # Make sure `ThermodynamicState` arguments are returned unchanged for FT in float_types rtol = FT(1e-2) + + _MSLP = FT(MSLP(param_set)) + e_int, ρ, q_tot, q_pt, T, p, θ_liq_ice = - MT.tested_convergence_range(param_set, 50, FT) + MT.tested_profiles(param_set, 50, FT) # PhaseDry ts = PhaseDry.(Ref(param_set), e_int, ρ) @@ -469,7 +583,8 @@ end @test all(air_density.(ts) .≈ ρ) # air_temperature_from_liquid_ice_pottemp_given_pressure-liquid_ice_pottemp inverse - θ_liq_ice_ = liquid_ice_pottemp_given_pressure.(Ref(param_set), T, p, q_pt) + θ_liq_ice_ = + liquid_ice_pottemp_given_pressure.(Ref(param_set), T, p, q_pt) @test all( air_temperature_from_liquid_ice_pottemp_given_pressure.( Ref(param_set), @@ -487,7 +602,10 @@ end p, q_pt, ) - @test all(liquid_ice_pottemp_given_pressure.(Ref(param_set), T, p, q_pt) .≈ θ_liq_ice) + @test all( + liquid_ice_pottemp_given_pressure.(Ref(param_set), T, p, q_pt) .≈ + θ_liq_ice, + ) # Accurate but expensive `LiquidIcePotTempSHumNonEquil` constructor (Non-linear temperature from θ_liq_ice) T_non_linear = @@ -500,7 +618,12 @@ end q_pt, ) T_expansion = - air_temperature_from_liquid_ice_pottemp.(Ref(param_set), θ_liq_ice, ρ, q_pt) + air_temperature_from_liquid_ice_pottemp.( + Ref(param_set), + θ_liq_ice, + ρ, + q_pt, + ) @test all(isapprox.(T_non_linear, T_expansion, rtol = rtol)) e_int_ = internal_energy.(Ref(param_set), T_non_linear, q_pt) ts = PhaseNonEquil.(Ref(param_set), e_int_, ρ, q_pt) @@ -508,7 +631,15 @@ end @test all(θ_liq_ice .≈ liquid_ice_pottemp.(ts)) # LiquidIcePotTempSHumEquil - ts = LiquidIcePotTempSHumEquil.(Ref(param_set), θ_liq_ice, ρ, q_tot, 40, FT(1e-3)) + ts = + LiquidIcePotTempSHumEquil.( + Ref(param_set), + θ_liq_ice, + ρ, + q_tot, + 40, + FT(1e-3), + ) @test all(isapprox.(liquid_ice_pottemp.(ts), θ_liq_ice, atol = 1e-1)) @test all(isapprox.(air_density.(ts), ρ, rtol = rtol)) @test all(getproperty.(PhasePartition.(ts), :tot) .≈ q_tot) @@ -532,10 +663,16 @@ end @test all( getproperty.(PhasePartition.(ts), :tot) .≈ getproperty.(q_pt, :tot), ) - @test all(isapprox.(air_pressure.(ts), p, atol = FT(MSLP(param_set)) * 2e-2)) + @test all(isapprox.(air_pressure.(ts), p, atol = _MSLP * 2e-2)) # LiquidIcePotTempSHumNonEquil_given_pressure - ts = LiquidIcePotTempSHumNonEquil_given_pressure.(Ref(param_set), θ_liq_ice, p, q_pt) + ts = + LiquidIcePotTempSHumNonEquil_given_pressure.( + Ref(param_set), + θ_liq_ice, + p, + q_pt, + ) @test all(liquid_ice_pottemp.(ts) .≈ θ_liq_ice) @test all(air_pressure.(ts) .≈ p) @test all( @@ -549,7 +686,15 @@ end ) # LiquidIcePotTempSHumNonEquil - ts = LiquidIcePotTempSHumNonEquil.(Ref(param_set), θ_liq_ice, ρ, q_pt, 5, FT(1e-3)) + ts = + LiquidIcePotTempSHumNonEquil.( + Ref(param_set), + θ_liq_ice, + ρ, + q_pt, + 5, + FT(1e-3), + ) @test all(θ_liq_ice .≈ liquid_ice_pottemp.(ts)) @test all(air_density.(ts) .≈ ρ) @test all( @@ -566,12 +711,13 @@ end end -@testset "Type-stability" begin +@testset "moist thermodynamics - type-stability" begin # NOTE: `Float32` saturation adjustment tends to have more difficulty # with converging to the same tolerances as `Float64`, so they're relaxed here. FT = Float32 - e_int, ρ, q_tot, q_pt, T, p, θ_liq_ice = MT.tested_convergence_range(param_set, 50, FT) + e_int, ρ, q_tot, q_pt, T, p, θ_liq_ice = + MT.tested_profiles(param_set, 50, FT) ρu = FT[1.0, 2.0, 3.0] e_pot = FT(100.0) @@ -590,7 +736,14 @@ end ) ts_neq = PhaseNonEquil.(Ref(param_set), e_int, ρ, q_pt) ts_θ_liq_ice_eq = - LiquidIcePotTempSHumEquil.(Ref(param_set), θ_liq_ice, ρ, q_tot, 40, FT(1e-3)) + LiquidIcePotTempSHumEquil.( + Ref(param_set), + θ_liq_ice, + ρ, + q_tot, + 40, + FT(1e-3), + ) ts_θ_liq_ice_eq_p = LiquidIcePotTempSHumEquil_given_pressure.( Ref(param_set), @@ -600,9 +753,15 @@ end 40, FT(1e-3), ) - ts_θ_liq_ice_neq = LiquidIcePotTempSHumNonEquil.(Ref(param_set), θ_liq_ice, ρ, q_pt) + ts_θ_liq_ice_neq = + LiquidIcePotTempSHumNonEquil.(Ref(param_set), θ_liq_ice, ρ, q_pt) ts_θ_liq_ice_neq_p = - LiquidIcePotTempSHumNonEquil_given_pressure.(Ref(param_set), θ_liq_ice, p, q_pt) + LiquidIcePotTempSHumNonEquil_given_pressure.( + Ref(param_set), + θ_liq_ice, + p, + q_pt, + ) for ts in ( ts_dry, @@ -645,10 +804,11 @@ end end -@testset "Dry limit" begin +@testset "moist thermodynamics - dry limit" begin FT = Float64 - e_int, ρ, q_tot, q_pt, T, p, θ_liq_ice = MT.tested_convergence_range(param_set, 50, FT) + e_int, ρ, q_tot, q_pt, T, p, θ_liq_ice = + MT.tested_profiles(param_set, 50, FT) # PhasePartition test is noisy, so do this only once: ts_dry = PhaseDry(param_set, first(e_int), first(ρ))