diff --git a/docs/src/breeze.bib b/docs/src/breeze.bib index 071e2ca2..5933f52b 100644 --- a/docs/src/breeze.bib +++ b/docs/src/breeze.bib @@ -28,3 +28,14 @@ @book{Straka2009 year = {2009}, doi = {10.1017/CBO9780511581168} } + +@article{Pressel2015, + title={Large-eddy simulation in an anelastic framework with closed water and entropy balances}, + author={Pressel, Kyle G and Kaul, Colleen M and Schneider, Tapio and Tan, Zhihong and Mishra, Siddhartha}, + journal={Journal of Advances in Modeling Earth Systems}, + volume={7}, + number={3}, + pages={1425--1456}, + year={2015}, + publisher={Wiley Online Library} +} \ No newline at end of file diff --git a/docs/src/index.md b/docs/src/index.md index 7ab44e6a..7e80d990 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -39,28 +39,27 @@ seed!(42) Nx = Nz = 64 Lz = 4 * 1024 -grid = RectilinearGrid(CPU(), size=(Nx, Nz), x=(0, 2Lz), z=(0, Lz), topology=(Periodic, Flat, Bounded)) +grid = RectilinearGrid(size=(Nx, Nz), x=(0, 2Lz), z=(0, Lz), topology=(Periodic, Flat, Bounded)) -reference_constants = Breeze.Thermodynamics.ReferenceStateConstants(base_pressure=1e5, potential_temperature=288) -buoyancy = Breeze.MoistAirBuoyancy(; reference_constants) +p₀, θᵣ = 1e5, 288 # reference state parameters +buoyancy = Breeze.MoistAirBuoyancy(grid, base_pressure=p₀, reference_potential_temperature=θᵣ) -Q₀ = 1000 # heat flux in W / m² -ρ₀ = Breeze.MoistAirBuoyancies.base_density(buoyancy) # air density at z=0 +thermo = buoyancy.thermodynamics +ρ₀ = Breeze.Thermodynamics.base_density(p₀, θᵣ, thermo) cₚ = buoyancy.thermodynamics.dry_air.heat_capacity +Q₀ = 1000 # heat flux in W / m² θ_bcs = FieldBoundaryConditions(bottom=FluxBoundaryCondition(Q₀ / (ρ₀ * cₚ))) -q_bcs = FieldBoundaryConditions(bottom=FluxBoundaryCondition(1e-2)) +qᵗ_bcs = FieldBoundaryConditions(bottom=FluxBoundaryCondition(1e-2)) advection = WENO() -tracers = (:θ, :q) model = NonhydrostaticModel(; grid, advection, buoyancy, - tracers = (:θ, :q), - boundary_conditions = (θ=θ_bcs, q=q_bcs)) + tracers = (:θ, :qᵗ), + boundary_conditions = (θ=θ_bcs, qᵗ=qᵗ_bcs)) Δθ = 2 # ᵒK -Tₛ = reference_constants.reference_potential_temperature # K +Tₛ = buoyancy.reference_state.potential_temperature # K θᵢ(x, z) = Tₛ + Δθ * z / grid.Lz + 2e-2 * Δθ * (rand() - 0.5) -qᵢ(x, z) = 0 # 1e-2 + 1e-5 * rand() -set!(model, θ=θᵢ, q=qᵢ) +set!(model, θ=θᵢ) simulation = Simulation(model, Δt=10, stop_time=2hours) conjure_time_step_wizard!(simulation, cfl=0.7) diff --git a/docs/src/microphysics/saturation_adjustment.md b/docs/src/microphysics/saturation_adjustment.md index e88fb7a1..1fc4fa8d 100644 --- a/docs/src/microphysics/saturation_adjustment.md +++ b/docs/src/microphysics/saturation_adjustment.md @@ -1,11 +1,11 @@ # Warm-phase saturation adjustment Warm-phase saturation adjustment is a model for water droplet nucleation that assumes that water vapor in excess of the saturation specific humidity is instantaneously converted to liquid water. -Mixed-phase saturation adjustment is described by [Chammas2023](@citet). +Mixed-phase saturation adjustment is described by [Pressel2015](@citet). Saturation adjustment may be formulated as a nonlinear algebraic equation that relates temperature, potential temperature, and total specific humidity, derived from the definition of liquid potential temperature, ```math -θ = \frac{T}{Π} \left \{1 - \frac{ℒᵥ₀}{cᵖᵐ T} \max \left [0, qᵗ - qᵛ⁺(T) \right ] \right \} , +θ = \frac{1}{Π} \left \{T - \frac{ℒᵥ₀}{cᵖᵐ} \max \left [0, qᵗ - qᵛ⁺(T) \right ] \right \} , ``` where ``Π`` is the Exner function, ``θ`` is potential temperature, ``T`` is temperature, @@ -32,14 +32,15 @@ The saturation specific humidity is then ```@example microphysics using Breeze -using Breeze.MoistAirBuoyancies: saturation_specific_humidity, HeightReferenceThermodynamicState +using Breeze.Thermodynamics: saturation_specific_humidity thermo = ThermodynamicConstants() -ref = ReferenceStateConstants(base_pressure=101325, potential_temperature=288) -z = 0.0 # [m] height -θ = 290.0 # [ᵒK] potential temperature -qᵛ⁺₀ = saturation_specific_humidity(θ, z, ref, thermo, thermo.liquid) +p₀ = 101325.0 +θ₀ = 288.0 +Rᵈ = Breeze.Thermodynamics.dry_air_gas_constant(thermo) +ρ₀ = p₀ / (Rᵈ * θ₀) +qᵛ⁺₀ = saturation_specific_humidity(θ₀, ρ₀, thermo, thermo.liquid) ``` Recall that the specific humidity is unitless, or has units "``kg / kg``": kg of water vapor @@ -49,17 +50,20 @@ given a total specific humidity slightly above saturation: ```@example microphysics using Breeze.MoistAirBuoyancies: temperature +using Breeze.Thermodynamics: PotentialTemperatureState, MoistureMassFractions +z = 0.0 qᵗ = 0.012 # [kg kg⁻¹] total specific humidity -U = HeightReferenceThermodynamicState(θ, qᵗ, z) -T = temperature(U, ref, thermo) +q = MoistureMassFractions(qᵗ, zero(qᵗ), zero(qᵗ)) +U = PotentialTemperatureState(θ₀, q, z, p₀, p₀, ρ₀) +T = temperature(U, thermo) ``` Finally, we recover the amount of liquid condensate by subtracting the saturation specific humidity from the total: ```@example microphysics -qᵛ⁺ = saturation_specific_humidity(T, z, ref, thermo, thermo.liquid) +qᵛ⁺ = saturation_specific_humidity(T, ρ₀, thermo, thermo.liquid) qˡ = qᵗ - qᵛ⁺ ``` @@ -69,14 +73,15 @@ As a second example, we examine the dependence of temperature on total specific when the potential temperature is constant: ```@example microphysics -qᵗ = 0:1e-4:0.04 # [kg kg⁻¹] total specific humidity -U = [HeightReferenceThermodynamicState(θ, qᵗⁱ, z) for qᵗⁱ in qᵗ] -T = [temperature(Uⁱ, ref, thermo) for Uⁱ in U] +qᵗ = 0:1e-4:0.035 # [kg kg⁻¹] total specific humidity +q = [MoistureMassFractions(qᵗⁱ, 0.0, 0.0) for qᵗⁱ in qᵗ] +U = [PotentialTemperatureState(θ₀, qⁱ, z, p₀, p₀, ρ₀) for qⁱ in q] +T = [temperature(Uⁱ, thermo) for Uⁱ in U] ## Compare with a simple piecewise linear model -ℒᵥ₀ = thermo.liquid.latent_heat +ℒᵥ₀ = thermo.liquid.reference_latent_heat cᵖᵈ = thermo.dry_air.heat_capacity -T̃ = [290 + ℒᵥ₀ / cᵖᵈ * max(0, qᵗⁱ - qᵛ⁺₀) for qᵗⁱ in qᵗ] +T̃ = [288 + ℒᵥ₀ / cᵖᵈ * max(0, qᵗⁱ - qᵛ⁺₀) for qᵗⁱ in qᵗ] using CairoMakie @@ -94,13 +99,32 @@ For a third example, we consider a state with constant potential temperature and but at varying heights: ```@example microphysics +grid = RectilinearGrid(size=100, z=(0, 1e4), topology=(Flat, Flat, Bounded)) +thermo = ThermodynamicConstants() +reference_state = ReferenceState(grid, thermo) + +θᵣ = reference_state.potential_temperature +p₀ = reference_state.base_pressure qᵗ = 0.005 -z = 0:100:10e3 +q = MoistureMassFractions(qᵗ, 0.0, 0.0) + +z = znodes(grid, Center()) +T = zeros(grid.Nz) +qᵛ⁺ = zeros(grid.Nz) +qˡ = zeros(grid.Nz) +rh = zeros(grid.Nz) + +for k = 1:grid.Nz + ρᵣ = reference_state.density[1, 1, k] + pᵣ = reference_state.pressure[1, 1, k] + + U = PotentialTemperatureState(θᵣ, q, z[k], p₀, pᵣ, ρᵣ) + T[k] = temperature(U, thermo) -T = [temperature(HeightReferenceThermodynamicState(θ, qᵗ, zᵏ), ref, thermo) for zᵏ in z] -qᵛ⁺ = [saturation_specific_humidity(T[k], z[k], ref, thermo, thermo.liquid) for k = 1:length(z)] -qˡ = [max(0, qᵗ - qᵛ⁺ᵏ) for qᵛ⁺ᵏ in qᵛ⁺] -rh = [100 * min(qᵗ, qᵛ⁺ᵏ) / qᵛ⁺ᵏ for qᵛ⁺ᵏ in qᵛ⁺] + qᵛ⁺[k] = saturation_specific_humidity(T[k], ρᵣ, thermo, thermo.liquid) + qˡ[k] = max(0, qᵗ - qᵛ⁺[k]) + rh[k] = 100 * min(qᵗ, qᵛ⁺[k]) / qᵛ⁺[k] +end cᵖᵈ = thermo.dry_air.heat_capacity g = thermo.gravitational_acceleration diff --git a/docs/src/thermodynamics.md b/docs/src/thermodynamics.md index 046a331d..0de463a9 100644 --- a/docs/src/thermodynamics.md +++ b/docs/src/thermodynamics.md @@ -187,23 +187,19 @@ with ``Rᵈ = 286.71 \; \mathrm{J} \, \mathrm{K}^{-1}``: ```@example reference_state using Breeze -using Breeze.Thermodynamics: reference_pressure, reference_density using CairoMakie -thermo = ThermodynamicConstants() -constants = ReferenceStateConstants(base_pressure=101325, potential_temperature=288) grid = RectilinearGrid(size=160, z=(0, 12_000), topology=(Flat, Flat, Bounded)) +thermo = ThermodynamicConstants() +reference_state = ReferenceState(grid, thermo, base_pressure=101325, potential_temperature=288) -pᵣ = CenterField(grid) -ρᵣ = CenterField(grid) - -set!(pᵣ, z -> reference_pressure(z, constants, thermo)) -set!(ρᵣ, z -> reference_density(z, constants, thermo)) +pᵣ = reference_state.pressure +ρᵣ = reference_state.density Rᵈ = Breeze.Thermodynamics.dry_air_gas_constant(thermo) cᵖᵈ = thermo.dry_air.heat_capacity -p₀ = constants.base_pressure -θ₀ = constants.reference_potential_temperature +p₀ = reference_state.base_pressure +θ₀ = reference_state.potential_temperature g = thermo.gravitational_acceleration # Verify that Tᵣ = θ₀ (1 - g z / (cᵖᵈ θ₀)) @@ -280,17 +276,18 @@ More generally, ``qᵈ = 1 - qᵛ - qᶜ``, where ``qᶜ`` is the total mass ratio of condensed species. In most situations on Earth, ``qᶜ ≪ qᵛ``. ```@example thermo +using Breeze.Thermodynamics: MoistureMassFractions + # Compute mixture properties for air with 0.01 specific humidity -qᵛ = 0.01 # 1% water vapor by mass -Rᵐ = mixture_gas_constant(qᵛ, thermo) +qᵗ = 0.01 # 1% water vapor by mass +q = MoistureMassFractions(qᵗ, zero(qᵗ), zero(qᵗ)) +Rᵐ = mixture_gas_constant(q, thermo) ``` We likewise define a mixture heat capacity via ``cᵖᵐ = qᵈ cᵖᵈ + qᵛ cᵖᵛ``, - ```@example thermo -q = 0.01 # 1% water vapor by mass -cᵖᵐ = mixture_heat_capacity(qᵛ, thermo) +cᵖᵐ = mixture_heat_capacity(q, thermo) ``` ## Liquid-ice potential temperature @@ -313,7 +310,7 @@ heats, the latent heat of a phase transition is linear in temperature. For example, for phase change from vapor to liquid, ```math -ℒˡ(T) = ℒˡ(T=0) + \big ( \underbrace{cᵖᵛ - cᵖˡ}_{≡Δcˡ} \big ) T , +ℒˡ(T) = ℒˡ(T=0) + \big ( \underbrace{cᵖᵛ - cˡ}_{≡Δcˡ} \big ) T , ``` where ``ℒˡ(T=0)`` is the latent heat at absolute zero, ``T = 0 \; \mathrm{K}``. @@ -327,13 +324,13 @@ Consider parameters for liquid water, ```@example thermo using Breeze.Thermodynamics: CondensedPhase -liquid_water = CondensedPhase(latent_heat=2500800, heat_capacity=4181) +liquid_water = CondensedPhase(reference_latent_heat=2500800, heat_capacity=4181) ``` or water ice, ```@example thermo -water_ice = CondensedPhase(latent_heat=2834000, heat_capacity=2108) +water_ice = CondensedPhase(reference_latent_heat=2834000, heat_capacity=2108) ``` The saturation vapor pressure is @@ -369,14 +366,19 @@ and this is what it looks like: ```@example using Breeze -using Breeze.MoistAirBuoyancies: saturation_specific_humidity +using Breeze.Thermodynamics: saturation_specific_humidity thermo = ThermodynamicConstants() -ref = ReferenceStateConstants(base_pressure=101325, potential_temperature=288) -z = 0 +p₀ = 101325 +Rᵈ = Breeze.Thermodynamics.dry_air_gas_constant(thermo) T = collect(273.2:0.1:313.2) -qᵛ⁺ = [saturation_specific_humidity(Tⁱ, z, ref, thermo, thermo.liquid) for Tⁱ in T] +qᵛ⁺ = zeros(length(T)) + +for i = 1:length(T) + ρ = p₀ / (Rᵈ * T[i]) + qᵛ⁺[i] = saturation_specific_humidity(T[i], ρ, thermo, thermo.liquid) +end using CairoMakie diff --git a/examples/JRA55_saturation_specific_humidity.jl b/examples/JRA55_saturation_specific_humidity.jl index a584f449..d04641f8 100644 --- a/examples/JRA55_saturation_specific_humidity.jl +++ b/examples/JRA55_saturation_specific_humidity.jl @@ -24,7 +24,7 @@ ax4 = Axis(fig[2, 2], title="Liquid specific humidity") # Compute cloudiness for instantaneous drop θ⁻ = T .- 10 Ψ = Breeze.ThermodynamicState{Float64}.(θ⁻, q, 0) -ℛ = Breeze.ReferenceStateConstants{Float64}(101325, 20) +ℛ = Breeze.ReferenceState{Float64}(101325, 20) T⁻ = Breeze.temperature.(Ψ, Ref(ℛ), Ref(thermo)) qᵛ★ = saturation_specific_humidity.(T⁻, 1.2, Ref(thermo)) qˡ = @. max(0, q - qᵛ★) diff --git a/examples/atmos_convection.jl b/examples/atmos_convection.jl index 2bafd187..b5940b5d 100644 --- a/examples/atmos_convection.jl +++ b/examples/atmos_convection.jl @@ -11,7 +11,7 @@ grid = RectilinearGrid(size=(Nx, Nz), x=(0, Lx), z=(0, Lz), topology=(Periodic, ρe_bcs = FieldBoundaryConditions(bottom=FluxBoundaryCondition(1000)) model = AtmosphereModel(grid, advection=WENO(), boundary_conditions=(; ρe=ρe_bcs)) -Tₛ = model.formulation.constants.reference_potential_temperature +Tₛ = model.formulation.reference_state.potential_temperature g = model.thermodynamics.gravitational_acceleration N² = 1e-6 dθdz = N² * Tₛ / g # Background potential temperature gradient diff --git a/examples/bomex.jl b/examples/bomex.jl index 23ffed31..4667ad77 100644 --- a/examples/bomex.jl +++ b/examples/bomex.jl @@ -42,14 +42,13 @@ u_bomex = AtmosphericProfilesLibrary.Bomex_u(FT) p₀ = 101500 # Pa θ₀ = 299.1 # K -reference_constants = Breeze.Thermodynamics.ReferenceStateConstants(base_pressure=p₀, potential_temperature=θ₀) -buoyancy = Breeze.MoistAirBuoyancy(; reference_constants) #, microphysics) +buoyancy = Breeze.MoistAirBuoyancy(grid, base_pressure=p₀, reference_potential_temperature=θ₀) # Simple precipitation scheme from CloudMicrophysics FT = eltype(grid) microphysics = CloudMicrophysics.Parameters.Parameters0M{FT}(τ_precip=600, S_0=0, qc_0=0.02) -@inline precipitation(x, y, z, t, q, params) = remove_precipitation(params, q, 0) -q_precip_forcing = Forcing(precipitation, field_dependencies=:q, parameters=microphysics) +@inline precipitation(x, y, z, t, qᵗ, params) = remove_precipitation(params, qᵗ, 0) +q_precip_forcing = Forcing(precipitation, field_dependencies=:qᵗ, parameters=microphysics) θ_bcs = FieldBoundaryConditions(bottom=FluxBoundaryCondition(8e-3)) q_bcs = FieldBoundaryConditions(bottom=FluxBoundaryCondition(5.2e-5)) @@ -143,17 +142,11 @@ set!(Fθ_field, z -> dTdt_bomex(1, z)) θ_radiation_forcing = Forcing(Fθ_field) θ_forcing = (θ_radiation_forcing, θ_subsidence_forcing) -# Note that most of the models that participated in Siebesma et al 2003 -# use 2nd order advection together with either TKE or Smag-Lilly closure. -advection = WENO(order=9) -closure = nothing -# advection = Centered(order=2) -# closure = SmagorinskyLilly() - -model = NonhydrostaticModel(; grid, advection, buoyancy, coriolis, closure, - tracers = (:θ, :q), +model = NonhydrostaticModel(; grid, buoyancy, coriolis, + advection = WENO(order=5), + tracers = (:θ, :qᵗ), forcing = (; q=q_forcing, u=u_forcing, v=v_forcing, θ=θ_forcing), - boundary_conditions = (θ=θ_bcs, q=q_bcs, u=u_bcs, v=v_bcs)) + boundary_conditions = (θ=θ_bcs, qᵗ=q_bcs, u=u_bcs, v=v_bcs)) # Values for the initial perturbations can be found in Appendix B # of Siebesma et al 2003, 3rd paragraph @@ -162,7 +155,7 @@ qϵ = 2.5e-5 θᵢ(x, y, z) = θ_bomex(z) + θϵ * randn() qᵢ(x, y, z) = q_bomex(z) + qϵ * randn() uᵢ(x, y, z) = u_bomex(z) -set!(model, θ=θᵢ, q=qᵢ, u=uᵢ) +set!(model, θ=θᵢ, qᵗ=qᵢ, u=uᵢ) simulation = Simulation(model; Δt=1, stop_time) conjure_time_step_wizard!(simulation, cfl=0.7) @@ -171,7 +164,7 @@ conjure_time_step_wizard!(simulation, cfl=0.7) u_avg = Field(Average(model.velocities.u, dims=(1, 2))) v_avg = Field(Average(model.velocities.v, dims=(1, 2))) θ_avg = Field(Average(model.tracers.θ, dims=(1, 2))) -q_avg = Field(Average(model.tracers.q, dims=(1, 2))) +q_avg = Field(Average(model.tracers.qᵗ, dims=(1, 2))) function compute_averages!(sim) compute!(u_avg) @@ -190,7 +183,7 @@ add_callback!(simulation, compute_averages!) T = Breeze.TemperatureField(model) qˡ = Breeze.CondensateField(model, T) qᵛ⁺ = Breeze.SaturationField(model, T) -qᵛ = model.tracers.q - qˡ +qᵛ = model.tracers.qᵗ - qˡ rh = Field(qᵛ / qᵛ⁺) # relative humidity T_avg = Field(Average(T, dims=(1, 2))) @@ -246,8 +239,8 @@ function progress(sim) umax = maximum(abs, u_avg) vmax = maximum(abs, v_avg) - q = sim.model.tracers.q - qmax = maximum(q) + qᵗ = sim.model.tracers.qᵗ + qᵗmax = maximum(qᵗ) θ = sim.model.tracers.θ θmin = minimum(θ) @@ -256,8 +249,8 @@ function progress(sim) msg = @sprintf("Iter: %d, t: %s, Δt: %s, max|ū|: (%.2e, %.2e), max(rh): %.2f", iteration(sim), prettytime(sim), prettytime(sim.Δt), umax, vmax, rhmax) - msg *= @sprintf(", max(q): %.2e, max(qˡ): %.2e, extrema(θ): (%.3e, %.3e)", - qmax, qˡmax, θmin, θmax) + msg *= @sprintf(", max(qᵗ): %.2e, max(qˡ): %.2e, extrema(θ): (%.3e, %.3e)", + qᵗmax, qˡmax, θmin, θmax) @info msg @@ -296,14 +289,14 @@ using CairoMakie if get(ENV, "CI", "false") == "false" θt = FieldTimeSeries(averages_filename, "θ") Tt = FieldTimeSeries(averages_filename, "T") - qt = FieldTimeSeries(averages_filename, "q") + qᵗt = FieldTimeSeries(averages_filename, "qᵗ") qˡt = FieldTimeSeries(averages_filename, "qˡ") - times = qt.times + times = qᵗt.times Nt = length(θt) fig = Figure(size=(1200, 800), fontsize=12) axθ = Axis(fig[1, 1], xlabel="θ (ᵒK)", ylabel="z (m)") - axq = Axis(fig[1, 2], xlabel="q (kg/kg)", ylabel="z (m)") + axqᵗ = Axis(fig[1, 2], xlabel="qᵗ (kg/kg)", ylabel="z (m)") axT = Axis(fig[2, 1], xlabel="T (ᵒK)", ylabel="z (m)") axqˡ = Axis(fig[2, 2], xlabel="qˡ (kg/kg)", ylabel="z (m)") @@ -312,8 +305,8 @@ if get(ENV, "CI", "false") == "false" n = slider.value #Observable(length(θt)) θn = @lift interior(θt[$n], 1, 1, :) - qn = @lift interior(qt[$n], 1, 1, :) Tn = @lift interior(Tt[$n], 1, 1, :) + qᵗn = @lift interior(qᵗt[$n], 1, 1, :) qˡn = @lift interior(qˡt[$n], 1, 1, :) z = znodes(θt) title = @lift "t = $(prettytime(times[$n]))" @@ -321,8 +314,8 @@ if get(ENV, "CI", "false") == "false" fig[0, :] = Label(fig, title, fontsize=22, tellwidth=false) hmθ = lines!(axθ, θn, z) - hmq = lines!(axq, qn, z) hmT = lines!(axT, Tn, z) + hmqᵗ = lines!(axqᵗ, qᵗn, z) hmqˡ = lines!(axqˡ, qˡn, z) xlims!(axqˡ, -1e-4, 1.5e-3) diff --git a/examples/free_convection.jl b/examples/free_convection.jl index b64e4835..01c1fb21 100644 --- a/examples/free_convection.jl +++ b/examples/free_convection.jl @@ -5,60 +5,52 @@ using Breeze arch = CPU() -Nx = Nz = 128 +Nx = Nz = 64 Lz = 4 * 1024 grid = RectilinearGrid(arch, size=(Nx, Nz), x=(0, 2Lz), z=(0, Lz), topology=(Periodic, Flat, Bounded)) p₀ = 101325 # Pa θ₀ = 288 # K -reference_constants = Breeze.Thermodynamics.ReferenceStateConstants(base_pressure=p₀, potential_temperature=θ₀) -buoyancy = Breeze.MoistAirBuoyancy(; reference_constants) +buoyancy = Breeze.MoistAirBuoyancy(grid, base_pressure=p₀, reference_potential_temperature=θ₀) -# Simple precipitation scheme from CloudMicrophysics -using CloudMicrophysics -using CloudMicrophysics.Microphysics0M: remove_precipitation - -FT = eltype(grid) -microphysics = CloudMicrophysics.Parameters.Parameters0M{FT}(τ_precip=600, S_0=0, qc_0=0.02) -@inline precipitation(x, z, t, q, params) = remove_precipitation(params, q, 0) -q_forcing = Forcing(precipitation, field_dependencies=:q, parameters=microphysics) - -ρ₀ = Breeze.MoistAirBuoyancies.base_density(buoyancy) # air density at z=0 +θᵣ = buoyancy.reference_state.potential_temperature +p₀ = buoyancy.reference_state.base_pressure +Rᵈ = Breeze.Thermodynamics.dry_air_gas_constant(buoyancy.thermodynamics) +ρ₀ = p₀ / (Rᵈ * θᵣ) # air density at z=0 cₚ = buoyancy.thermodynamics.dry_air.heat_capacity Q₀ = 1000 # heat flux in W / m² Jθ = Q₀ / (ρ₀ * cₚ) # temperature flux θ_bcs = FieldBoundaryConditions(bottom=FluxBoundaryCondition(Jθ)) vapor_flux = FluxBoundaryCondition(1e-2) -q_bcs = FieldBoundaryConditions(bottom=vapor_flux) +qᵗ_bcs = FieldBoundaryConditions(bottom=vapor_flux) advection = WENO() #(momentum=WENO(), θ=WENO(), q=WENO(bounds=(0, 1))) tracers = (:θ, :q) model = NonhydrostaticModel(; grid, advection, buoyancy, - tracers = (:θ, :q), - forcing = (; q=q_forcing), - boundary_conditions = (θ=θ_bcs, q=q_bcs)) + tracers = (:θ, :qᵗ), + boundary_conditions = (θ=θ_bcs, qᵗ=qᵗ_bcs)) Lz = grid.Lz Δθ = 5 # K -Tₛ = reference_constants.reference_potential_temperature # K +Tₛ = buoyancy.reference_state.potential_temperature # K θᵢ(x, z) = Tₛ + Δθ * z / Lz + 1e-2 * Δθ * randn() -qᵢ(x, z) = 0 # 1e-2 + 1e-5 * rand() -set!(model, θ=θᵢ, q=qᵢ) +qᵗᵢ(x, z) = 0 # 1e-2 + 1e-5 * rand() +set!(model, θ=θᵢ, qᵗ=qᵗᵢ) -simulation = Simulation(model, Δt=10, stop_time=4hours) +simulation = Simulation(model, Δt=10, stop_iteration=1000) #stop_time=4hours) conjure_time_step_wizard!(simulation, cfl=0.7) T = Breeze.TemperatureField(model) qˡ = Breeze.CondensateField(model, T) -qᵛ★ = Breeze.SaturationField(model, T) -δ = Field(model.tracers.q - qᵛ★) +qᵛ⁺ = Breeze.SaturationField(model, T) +δ = Field(model.tracers.qᵗ - qᵛ⁺) function progress(sim) compute!(T) compute!(qˡ) compute!(δ) - q = sim.model.tracers.q + qᵗ = sim.model.tracers.qᵗ θ = sim.model.tracers.θ u, v, w = sim.model.velocities @@ -66,8 +58,8 @@ function progress(sim) vmax = maximum(v) wmax = maximum(w) - qmin = minimum(q) - qmax = maximum(q) + qᵗmin = minimum(qᵗ) + qᵗmax = maximum(qᵗ) qˡmax = maximum(qˡ) δmax = maximum(δ) @@ -77,8 +69,8 @@ function progress(sim) msg = @sprintf("Iter: %d, t: %s, Δt: %s, max|u|: (%.2e, %.2e, %.2e)", iteration(sim), prettytime(sim), prettytime(sim.Δt), umax, vmax, wmax) - msg *= @sprintf(", extrema(q): (%.2e, %.2e), max(qˡ): %.2e, min(δ): %.2e, extrema(θ): (%.2e, %.2e)", - qmin, qmax, qˡmax, δmax, θmin, θmax) + msg *= @sprintf(", max(qˡ): %.2e, min(δ): %.2e, extrema(θ): (%.2e, %.2e)", + qᵗmin, qᵗmax, qˡmax, δmax, θmin, θmax) @info msg @@ -88,8 +80,8 @@ end add_callback!(simulation, progress, IterationInterval(10)) using Oceananigans.Models: ForcingOperation -Sʳ = ForcingOperation(:q, model) -outputs = merge(model.velocities, model.tracers, (; T, qˡ, qᵛ★, Sʳ)) +Sʳ = ForcingOperation(:qᵗ, model) +outputs = merge(model.velocities, model.tracers, (; T, qˡ, qᵛ⁺, Sʳ)) ow = JLD2Writer(model, outputs, filename = "free_convection.jld2", diff --git a/examples/mountain_wave.jl b/examples/mountain_wave.jl index 20394d56..7f4aece4 100644 --- a/examples/mountain_wave.jl +++ b/examples/mountain_wave.jl @@ -19,7 +19,7 @@ grid = ImmersedBoundaryGrid(underlying_grid, PartialCellBottom(hill)) model = AtmosphereModel(grid, advection = WENO()) # Initial conditions -θ₀ = model.formulation.constants.reference_potential_temperature +θ₀ = model.formulation.constants.potential_temperature g = model.thermodynamics.gravitational_acceleration N² = 1e-4 # Brunt-Väisälä frequency squared (s⁻²) dθdz = N² * θ₀ / g # Background potential temperature gradient diff --git a/examples/thermal_bubble.jl b/examples/thermal_bubble.jl index 4b94c176..7e221465 100644 --- a/examples/thermal_bubble.jl +++ b/examples/thermal_bubble.jl @@ -1,5 +1,5 @@ # Thermal Bubble Simulation following Ahmad & Lindeman (2007) and Sridhar et al (2022) -# This simulates a circular potential temperature perturbation in 2D +# This simulates a circular moist static energy perturbation in 2D using Breeze using Oceananigans.Units @@ -25,33 +25,29 @@ grid = RectilinearGrid(arch, halo = (5, 5), topology = (Periodic, Flat, Bounded)) -# Thermodynamic setup -p₀ = 101325 # Pa - standard atmospheric pressure -θ₀ = 300.0 # K - reference potential temperature -reference_constants = ReferenceStateConstants(base_pressure=p₀, potential_temperature=θ₀) -buoyancy = MoistAirBuoyancy(; reference_constants) - -# Advection scheme - WENO for high-order accuracy advection = WENO(order=5) - -# Create the model model = AtmosphereModel(grid; advection) -# Thermal bubble parameters +# Thermal bubble parameters (moist static energy) x₀ = Lx / 2 # Center of bubble in x z₀ = 4e3 # Center of bubble in z (2 km height) r₀ = 2e3 # Initial radius of bubble (2 km) -Δθ = 2 # Potential temperature perturbation (ᵒK) +Δθ = 10 # K + +# Background stratification (used to construct a gently increasing MSE with height via θ) +N² = 1e-6 # Brunt-Väisälä frequency squared (s⁻²) +θ₀ = model.formulation.reference_state.potential_temperature +dθdz = N² * θ₀ / 9.81 # Background potential temperature gradient used to build MSE -# Background stratification -N² = 1e-8 # Brunt-Väisälä frequency squared (s⁻²) -dθdz = N² * θ₀ / 9.81 # Background potential temperature gradient +# Initial conditions (set moist static energy directly) +# We form MSE ρe ≈ ρᵣ cᵖᵈ θ with a localized perturbation. +ρʳ = model.formulation.reference_state.density +cᵖᵈ = model.thermodynamics.dry_air.heat_capacity -# Initial conditions function θᵢ(x, z) - θ̄ = θ₀ + dθdz * z # background stratification + θ̄ = θ₀ + dθdz * z # background potential temperature used to build MSE r = sqrt((x - x₀)^2 + (z - z₀)^2) # distance from bubble center - θ′ = Δθ * max(0, 1 - r / r₀) # bubble + θ′ = Δθ * max(0, 1 - r / r₀) return θ̄ + θ′ end @@ -67,7 +63,7 @@ conjure_time_step_wizard!(simulation, cfl=0.7) # Progress monitoring function progress(sim) ρe = sim.model.energy - u, w = sim.model.velocities + u, v, w = sim.model.velocities ρe_max = maximum(ρe) ρe_min = minimum(ρe) @@ -92,7 +88,7 @@ u, v, w = model.velocities # Temperature perturbation from background state ρE = Field{Nothing, Nothing, Center}(grid) -set!(E, Field(Average(model.energy, dims=(1, 2)))) +set!(ρE, Field(Average(model.energy, dims=(1, 2)))) ρe′ = model.energy - ρE ρe = model.energy T = model.temperature @@ -107,7 +103,7 @@ writer = JLD2Writer(model, outputs; filename, simulation.output_writers[:jld2] = writer @info "Running thermal bubble simulation on grid: $grid" -@info "Bubble parameters: center=($x₀, $z₀), radius=$r₀, Δθ=$Δθ K" +@info "Bubble parameters: center=($x₀, $z₀), radius=$r₀, Δe=$Δe J/m³" @info "Domain: $(Lx/1000) km × $(Lz/1000) km, resolution: $Nx × $Nz" run!(simulation) @@ -127,7 +123,7 @@ if get(ENV, "CI", "false") == "false" times = ρet.times Nt = length(ρet) - # Create visualization - expanded to 6 panels + # Create visualization - 4 panels fig = Figure(size=(1500, 1000), fontsize=12) # Subplot layout - 3 rows, 2 columns @@ -139,7 +135,7 @@ if get(ENV, "CI", "false") == "false" axw = Axis(fig[3, 2], xlabel="x (km)", ylabel="z (km)", title="Vertical Velocity w (m/s)") # Time slider - slider = Slider(fig[4, 1:2], range=1:Nt, startvalue=1) + slider = Slider(fig[3, 1:2], range=1:Nt, startvalue=1) n = slider.value # Observable fields @@ -163,14 +159,13 @@ if get(ENV, "CI", "false") == "false" ρe′_range = (minimum(ρe′t), maximum(ρe′t)) T_range = (minimum(Tt), maximum(Tt)) ζ_range = maximum(abs, ζt) - u_range = maximum(abs, ut) w_range = maximum(abs, wt) hme = heatmap!(axe, x, z, ρen, colorrange=ρe_range, colormap=:thermal) hme′ = heatmap!(axe′, x, z, ρe′n, colorrange=ρe′_range, colormap=:balance) hmT = heatmap!(axT, x, z, Tn, colorrange=T_range, colormap=:thermal) hmζ = heatmap!(axζ, x, z, ζn, colorrange=(-ζ_range, ζ_range), colormap=:balance) - hmu = heatmap!(axu, x, z, un, colorrange=(-u_range, u_range), colormap=:balance) + hmu = heatmap!(axu, x, z, un, colorrange=(-w_range, w_range), colormap=:balance) hmw = heatmap!(axw, x, z, wn, colorrange=(-w_range, w_range), colormap=:balance) # Add colorbars diff --git a/examples/time_step_atmos_model.jl b/examples/time_step_atmos_model.jl index daf098d4..c94512b3 100644 --- a/examples/time_step_atmos_model.jl +++ b/examples/time_step_atmos_model.jl @@ -15,7 +15,7 @@ model = AtmosphereModel(grid, advection=WENO(), boundary_conditions=(; e=e_bcs)) Lz = grid.Lz Δθ = 5 # K -Tₛ = model.formulation.constants.reference_potential_temperature +Tₛ = model.formulation.constants.potential_temperature # θᵢ(x, y, z) = Tₛ + Δθ * z / Lz qᵢ(x, y, z) = 0 Ξᵢ(x, y, z) = 1e-2 * randn() diff --git a/src/AtmosphereModels/anelastic_formulation.jl b/src/AtmosphereModels/anelastic_formulation.jl index 2045d97c..4f16f495 100644 --- a/src/AtmosphereModels/anelastic_formulation.jl +++ b/src/AtmosphereModels/anelastic_formulation.jl @@ -1,8 +1,8 @@ using ..Thermodynamics: + MoistureMassFractions, ThermodynamicConstants, - ReferenceStateConstants, - reference_pressure, - reference_density, + ReferenceState, + AnelasticThermodynamicState, mixture_gas_constant, mixture_heat_capacity, dry_air_gas_constant @@ -22,17 +22,15 @@ import Oceananigans.TimeSteppers: compute_pressure_correction!, make_pressure_co ##### Formulation definition ##### -struct AnelasticFormulation{FT, F} - constants :: ReferenceStateConstants{FT} - reference_pressure :: F - reference_density :: F +struct AnelasticFormulation{R} + reference_state :: R end const AnelasticModel = AtmosphereModel{<:AnelasticFormulation} function Base.summary(formulation::AnelasticFormulation) p₀ = formulation.constants.base_pressure - θᵣ = formulation.constants.reference_potential_temperature + θᵣ = formulation.constants.potential_temperature return string("AnelasticFormulation(p₀=", prettysummary(p₀), ", θᵣ=", prettysummary(θᵣ), ")") end @@ -41,14 +39,6 @@ Base.show(io::IO, formulation::AnelasticFormulation) = print(io, "AnelasticFormu field_names(::AnelasticFormulation, tracer_names) = (:ρu, :ρv, :ρw, :ρe, :ρqᵗ, tracer_names...) -struct AnelasticThermodynamicState{FT} - potential_temperature :: FT - specific_humidity :: FT - reference_density :: FT - reference_pressure :: FT - exner_function :: FT -end - function AnelasticFormulation(grid, state_constants, thermo) pᵣ = Field{Nothing, Nothing, Center}(grid) ρᵣ = Field{Nothing, Nothing, Center}(grid) @@ -66,19 +56,20 @@ end function thermodynamic_state(i, j, k, grid, formulation::AnelasticFormulation, thermo, energy, absolute_humidity) @inbounds begin e = energy[i, j, k] - pᵣ = formulation.reference_pressure[1, 1, k] - ρᵣ = formulation.reference_density[1, 1, k] - ρq = absolute_humidity[i, j, k] + pᵣ = formulation.reference_state.pressure[i, j, k] + ρᵣ = formulation.reference_state.density[i, j, k] + ρqᵗ = absolute_humidity[i, j, k] end cᵖᵈ = thermo.dry_air.heat_capacity θ = e / (cᵖᵈ * ρᵣ) - q = ρq / ρᵣ + qᵗ = ρqᵗ / ρᵣ + q = MoistureMassFractions(qᵗ, zero(qᵗ), zero(qᵗ)) # assuming non-condensed state Rᵐ = mixture_gas_constant(q, thermo) cᵖᵐ = mixture_heat_capacity(q, thermo) - p₀ = formulation.constants.base_pressure + p₀ = formulation.reference_state.base_pressure Π = (pᵣ / p₀)^(Rᵐ / cᵖᵐ) return AnelasticThermodynamicState(θ, q, ρᵣ, pᵣ, Π) @@ -86,20 +77,29 @@ end @inline function specific_volume(i, j, k, grid, formulation, temperature, specific_humidity, thermo) @inbounds begin - q = specific_humidity[i, j, k] - pᵣ = formulation.reference_pressure[1, 1, k] + qᵗ = specific_humidity[i, j, k] + pᵣ = formulation.reference_state.pressure[i, j, k] T = temperature[i, j, k] end + # TODO: fix this assumption of non-condensed state + q = MoistureMassFractions(qᵗ, zero(qᵗ), zero(qᵗ)) Rᵐ = mixture_gas_constant(q, thermo) return Rᵐ * T / pᵣ end +@inline function buoyancy(i, j, k, grid, formulation, T, qᵗ, thermo) + α = specific_volume(i, j, k, grid, formulation, T, qᵗ, thermo) + αᵣ = reference_specific_volume(i, j, k, grid, formulation, thermo) + g = thermo.gravitational_acceleration + return g * (α - αᵣ) / αᵣ +end + @inline function reference_specific_volume(i, j, k, grid, formulation, thermo) Rᵈ = dry_air_gas_constant(thermo) - pᵣ = @inbounds formulation.reference_pressure[1, 1, k] - θᵣ = formulation.constants.reference_potential_temperature + pᵣ = @inbounds formulation.reference_state.pressure[i, j, k] + θᵣ = formulation.reference_state.potential_temperature return Rᵈ * θᵣ / pᵣ end @@ -143,7 +143,7 @@ end tridiagonal_direction(formulation::AnelasticTridiagonalSolverFormulation) = ZDirection() function formulation_pressure_solver(anelastic_formulation::AnelasticFormulation, grid) - reference_density = anelastic_formulation.reference_density + reference_density = anelastic_formulation.reference_state.density tridiagonal_formulation = AnelasticTridiagonalSolverFormulation(reference_density) solver = if grid isa Oceananigans.ImmersedBoundaries.ImmersedBoundaryGrid @@ -212,7 +212,7 @@ function compute_pressure_correction!(model::AnelasticModel, Δt) foreach(mask_immersed_field!, model.momentum) fill_halo_regions!(model.momentum, model.clock, fields(model)) - ρᵣ = model.formulation.reference_density + ρᵣ = model.formulation.reference_state.density ρŨ = model.momentum solver = model.pressure_solver pₙ = model.nonhydrostatic_pressure @@ -284,7 +284,7 @@ function make_pressure_correction!(model::AnelasticModel, Δt) model.grid, Δt, model.nonhydrostatic_pressure, - model.formulation.reference_density) + model.formulation.reference_state.density) return nothing end diff --git a/src/AtmosphereModels/atmosphere_model.jl b/src/AtmosphereModels/atmosphere_model.jl index f6ffdfde..1edce098 100644 --- a/src/AtmosphereModels/atmosphere_model.jl +++ b/src/AtmosphereModels/atmosphere_model.jl @@ -1,11 +1,4 @@ -using ..Thermodynamics: - ThermodynamicConstants, - ReferenceStateConstants, - reference_pressure, - reference_density, - mixture_gas_constant, - mixture_heat_capacity, - dry_air_gas_constant +using ..Thermodynamics: ThermodynamicConstants, ReferenceState using Oceananigans: AbstractModel, Center, CenterField, Clock, Field using Oceananigans: WENO, XFaceField, YFaceField, ZFaceField @@ -60,11 +53,8 @@ mutable struct AtmosphereModel{Frm, Arc, Tst, Grd, Clk, Thm, Den, Mom, Eng, Wat, end function default_formulation(grid, thermo) - FT = eltype(grid) - base_pressure = convert(FT, 101325) - potential_temperature = convert(FT, 288) - constants = ReferenceStateConstants(base_pressure, potential_temperature) - return AnelasticFormulation(grid, constants, thermo) + reference_state = ReferenceState(grid, thermo) + return AnelasticFormulation(reference_state) end """ diff --git a/src/AtmosphereModels/dynamics_kernel_functions.jl b/src/AtmosphereModels/dynamics_kernel_functions.jl index 8a454732..0891b55a 100644 --- a/src/AtmosphereModels/dynamics_kernel_functions.jl +++ b/src/AtmosphereModels/dynamics_kernel_functions.jl @@ -6,13 +6,6 @@ using Oceananigans.Operators: ∂xᶠᶜᶜ, ∂yᶜᶠᶜ, ∂zᶜᶜᶠ, ℑz ##### Some key functions ##### -@inline function buoyancy(i, j, k, grid, formulation, temperature, specific_humidity, thermo) - α = specific_volume(i, j, k, grid, formulation, temperature, specific_humidity, thermo) - αᵣ = reference_specific_volume(i, j, k, grid, formulation, thermo) - g = thermo.gravitational_acceleration - return g * (α - αᵣ) / αᵣ -end - @inline function ρ_bᶜᶜᶠ(i, j, k, grid, ρ, T, q, formulation, thermo) ρᶜᶜᶠ = ℑzᵃᵃᶠ(i, j, k, grid, ρ) bᶜᶜᶠ = ℑzᵃᵃᶠ(i, j, k, grid, buoyancy, formulation, T, q, thermo) diff --git a/src/AtmosphereModels/saturation_adjustment.jl b/src/AtmosphereModels/saturation_adjustment.jl index 05cfc214..cc193e16 100644 --- a/src/AtmosphereModels/saturation_adjustment.jl +++ b/src/AtmosphereModels/saturation_adjustment.jl @@ -22,7 +22,7 @@ end r₁ = saturation_adjustment_residual(T₁, qˡ₁, state, thermo) q = state.specific_humidity - ℒ = thermo.liquid.latent_heat + ℒ = thermo.liquid.reference_latent_heat cᵖᵐ = mixture_heat_capacity(q, thermo) T₂ = (T₁ + sqrt(T₁^2 + 4 * ℒ * qˡ₁ / cᵖᵐ)) / 2 qˡ₂ = condensate_specific_humidity(T₂, state, thermo) @@ -53,7 +53,7 @@ end end @inline function saturation_adjustment_residual(T, qˡ, state, thermo) - ℒᵛ = thermo.liquid.latent_heat + ℒᵛ = thermo.liquid.reference_latent_heat q = state.specific_humidity θ = state.potential_temperature Π = state.exner_function diff --git a/src/AtmosphereModels/set_atmosphere_model.jl b/src/AtmosphereModels/set_atmosphere_model.jl index 3d5eef9a..f9fb75cf 100644 --- a/src/AtmosphereModels/set_atmosphere_model.jl +++ b/src/AtmosphereModels/set_atmosphere_model.jl @@ -1,21 +1,47 @@ +using Oceananigans.Grids: znode, Center +using Oceananigans.TimeSteppers: update_state! using Oceananigans.BoundaryConditions: fill_halo_regions! using Oceananigans.TimeSteppers: compute_pressure_correction!, make_pressure_correction!, update_state! +using ..Thermodynamics: + PotentialTemperatureState, + exner_function, + mixture_heat_capacity + import Oceananigans.Fields: set! +const c = Center() + +move_to_front(names, name) = tuple(name, filter(n -> n != name, names)...) + +function prioritize_names(names) + for n in (:w, :ρw, :v, :ρv, :u, :ρu, :qᵗ, :ρqᵗ) + if n ∈ names + names = move_to_front(names, n) + end + end + + return names +end + function set!(model::AtmosphereModel; enforce_mass_conservation=true, kw...) - for (name, value) in kw + names = collect(keys(kw)) + prioritized = prioritize_names(names) + + for name in prioritized + value = kw[name] # Prognostic variables if name ∈ propertynames(model.momentum) - ϕ = getproperty(model.momentum, name) - set!(ϕ, value) + ρu = getproperty(model.momentum, name) + set!(ρu, value) elseif name ∈ propertynames(model.tracers) - ϕ = getproperty(model.tracers, name) + c = getproperty(model.tracers, name) + set!(c, value) elseif name == :ρe - ϕ = model.energy - elseif name == :ρq - ϕ = model.absolute_humidity + set!(model.energy, value) + elseif name == :ρqᵗ + set!(model.absolute_humidity, value) end # Setting diagnostic variables @@ -23,34 +49,37 @@ function set!(model::AtmosphereModel; enforce_mass_conservation=true, kw...) θ = model.temperature # use scratch set!(θ, value) - ρᵣ = model.formulation.reference_density - cᵖᵈ = model.thermodynamics.dry_air.heat_capacity - ϕ = model.energy - value = ρᵣ * cᵖᵈ * θ - elseif name == :q - q = model.specific_humidity - set!(q, value) - - ρᵣ = model.formulation.reference_density - ϕ = model.absolute_humidity - value = ρᵣ * q + grid = model.grid + arch = grid.architecture + thermo = model.thermodynamics + formulation = model.formulation + energy = model.energy + specific_humidity = model.specific_humidity + launch!(arch, grid, :xyz, _energy_from_potential_temperature!, energy, grid, + θ, specific_humidity, formulation, thermo) + + elseif name == :qᵗ + qᵗ = model.specific_humidity + set!(qᵗ, value) + ρᵣ = model.formulation.reference_state.density + ρqᵗ = model.absolute_humidity + set!(ρqᵗ, ρᵣ * qᵗ) + elseif name ∈ (:u, :v, :w) u = model.velocities[name] set!(u, value) - ρᵣ = model.formulation.reference_density + ρᵣ = model.formulation.reference_state.density ϕ = model.momentum[Symbol(:ρ, name)] value = ρᵣ * u + set!(ϕ, value) end - - set!(ϕ, value) - fill_halo_regions!(ϕ, model.clock, fields(model)) end # Apply a mask - # foreach(mask_immersed_field!, prognostic_fields(model)) + foreach(mask_immersed_field!, prognostic_fields(model)) update_state!(model, compute_tendencies=false) - + if enforce_mass_conservation FT = eltype(model.grid) Δt = one(FT) @@ -59,5 +88,38 @@ function set!(model::AtmosphereModel; enforce_mass_conservation=true, kw...) update_state!(model, compute_tendencies=false) end + fill_halo_regions!(model.energy) + return nothing end + +@kernel function _energy_from_potential_temperature!(moist_static_energy, grid, + potential_temperature, + specific_humidity, + formulation, + thermo) + i, j, k = @index(Global, NTuple) + + @inbounds begin + ρᵣ = formulation.reference_state.density[i, j, k] + qᵗ = specific_humidity[i, j, k] + pᵣ = formulation.reference_state.pressure[i, j, k] + θ = potential_temperature[i, j, k] + end + + p₀ = formulation.reference_state.base_pressure + z = znode(i, j, k, grid, c, c, c) + + # Assuming a state with no condensate? + q = MoistureMassFractions(qᵗ, zero(qᵗ), zero(qᵗ)) + + 𝒰 = PotentialTemperatureState(θ, q, z, p₀, pᵣ, ρᵣ) + Π = exner_function(𝒰, thermo) + T = Π * θ + + ℒ₀ = thermo.liquid.reference_latent_heat + g = thermo.gravitational_acceleration + cᵖᵐ = mixture_heat_capacity(q, thermo) + + @inbounds moist_static_energy[i, j, k] = ρᵣ * (cᵖᵐ * T + g * z + qᵗ * ℒ₀) +end \ No newline at end of file diff --git a/src/AtmosphereModels/update_atmosphere_model_state.jl b/src/AtmosphereModels/update_atmosphere_model_state.jl index 8a96711d..554121b1 100644 --- a/src/AtmosphereModels/update_atmosphere_model_state.jl +++ b/src/AtmosphereModels/update_atmosphere_model_state.jl @@ -1,5 +1,6 @@ using ..Thermodynamics: saturation_specific_humidity, + total_specific_humidity, mixture_heat_capacity, mixture_gas_constant @@ -64,8 +65,8 @@ end ρv = momentum.ρv[i, j, k] ρw = momentum.ρw[i, j, k] - ρᵣᵃᵃᶜ = formulation.reference_density[i, j, k] - ρᵣᵃᵃᶠ = ℑzᵃᵃᶠ(i, j, k, grid, formulation.reference_density) + ρᵣᵃᵃᶜ = formulation.reference_state.density[i, j, k] + ρᵣᵃᵃᶠ = ℑzᵃᵃᶠ(i, j, k, grid, formulation.reference_state.density) velocities.u[i, j, k] = ρu / ρᵣᵃᵃᶜ velocities.v[i, j, k] = ρv / ρᵣᵃᵃᶜ velocities.w[i, j, k] = ρw / ρᵣᵃᵃᶠ @@ -83,7 +84,7 @@ end i, j, k = @index(Global, NTuple) 𝒰 = thermodynamic_state(i, j, k, grid, formulation, thermo, energy, absolute_humidity) - @inbounds specific_humidity[i, j, k] = 𝒰.specific_humidity + @inbounds specific_humidity[i, j, k] = total_specific_humidity(𝒰) # Possibly perform saturation adjustment # Note, we will make this much prettier in the future @@ -113,7 +114,7 @@ function compute_tendencies!(model::AnelasticModel) fields(model)) pₕ′ = model.hydrostatic_pressure_anomaly - ρᵣ = model.formulation.reference_density + ρᵣ = model.formulation.reference_state.density u_args = tuple(common_args..., model.forcing.ρu, pₕ′, ρᵣ) v_args = tuple(common_args..., model.forcing.ρv, pₕ′, ρᵣ) w_args = tuple(common_args..., model.forcing.ρw, ρᵣ, @@ -133,10 +134,10 @@ function compute_tendencies!(model::AnelasticModel) model.specific_humidity, model.thermodynamics, model.condensates, model.microphysics) launch!(arch, grid, :xyz, compute_moist_static_energy_tendency!, Gρe, grid, ρe_args) - ρq = model.absolute_humidity + ρqᵗ = model.absolute_humidity Gρqᵗ = model.timestepper.Gⁿ.ρqᵗ Fρqᵗ = model.forcing.ρqᵗ - ρq_args = tuple(ρq, Fρqᵗ, scalar_args...) + ρq_args = tuple(ρqᵗ, Fρqᵗ, scalar_args...) launch!(arch, grid, :xyz, compute_scalar_tendency!, Gρqᵗ, grid, ρq_args) return nothing diff --git a/src/Breeze.jl b/src/Breeze.jl index 1eab2927..ddd8d13c 100644 --- a/src/Breeze.jl +++ b/src/Breeze.jl @@ -8,7 +8,7 @@ module Breeze export MoistAirBuoyancy, ThermodynamicConstants, - ReferenceStateConstants, + ReferenceState, AnelasticFormulation, AtmosphereModel, TemperatureField, diff --git a/src/MoistAirBuoyancies.jl b/src/MoistAirBuoyancies.jl index eb2c499f..a282f7cd 100644 --- a/src/MoistAirBuoyancies.jl +++ b/src/MoistAirBuoyancies.jl @@ -1,10 +1,21 @@ module MoistAirBuoyancies -export MoistAirBuoyancy -export UnsaturatedMoistAirBuoyancy -export TemperatureField -export CondensateField -export SaturationField +export + MoistAirBuoyancy, + TemperatureField, + CondensateField, + SaturationField + +using ..Thermodynamics: + PotentialTemperatureState, + MoistureMassFractions, + total_specific_humidity, + dry_air_gas_constant, + vapor_gas_constant, + with_moisture, + saturation_vapor_pressure, + density, + exner_function using Oceananigans: Oceananigans, Center, Field, KernelFunctionOperation using Oceananigans.Grids: AbstractGrid @@ -12,352 +23,373 @@ using Oceananigans.Operators: ∂zᶜᶜᶠ using Adapt: Adapt, adapt -import Oceananigans.BuoyancyFormulations: AbstractBuoyancyFormulation, - buoyancy_perturbationᶜᶜᶜ, - ∂z_b, - required_tracers +import Oceananigans.BuoyancyFormulations: + AbstractBuoyancyFormulation, + buoyancy_perturbationᶜᶜᶜ, + ∂z_b, + required_tracers + +import ..Thermodynamics: saturation_specific_humidity using ..Thermodynamics: ThermodynamicConstants, - ReferenceStateConstants, + ReferenceState, mixture_heat_capacity, - mixture_gas_constant, - reference_specific_volume, - reference_pressure + mixture_gas_constant -import ..Thermodynamics: - base_density, - saturation_specific_humidity, - condensate_specific_humidity - -struct MoistAirBuoyancy{FT, AT} <: AbstractBuoyancyFormulation{Nothing} - reference_constants :: ReferenceStateConstants{FT} +struct MoistAirBuoyancy{RS, AT} <: AbstractBuoyancyFormulation{Nothing} + reference_state :: RS thermodynamics :: AT end """ - MoistAirBuoyancy(FT=Oceananigans.defaults.FloatType; - thermodynamics = ThermodynamicConstants(FT), - reference_constants = ReferenceStateConstants{FT}(101325, 290)) + MoistAirBuoyancy(grid; + base_pressure = 101325, + reference_potential_temperature = 288, + thermodynamics = ThermodynamicConstants(FT)) -Return a MoistAirBuoyancy formulation that can be provided as input to an `AtmosphereModel` -or an `Oceananigans.NonhydrostaticModel`. +Return a MoistAirBuoyancy formulation that can be provided as input to an +`Oceananigans.NonhydrostaticModel`. !!! note "Required tracers" - `MoistAirBuoyancy` requires tracers `q` and `θ` to be included in the model. + `MoistAirBuoyancy` requires tracers `θ` and `qᵗ`. Example ======= -```jldoctest -julia> using Breeze, Oceananigans +```jldoctest mab +using Breeze, Oceananigans -julia> buoyancy = MoistAirBuoyancy() -MoistAirBuoyancy -├── reference_constants: Breeze.Thermodynamics.ReferenceStateConstants{Float64} -└── thermodynamics: ThermodynamicConstants +grid = RectilinearGrid(size=(1, 1, 8), extent=(1, 1, 3e3)) +buoyancy = MoistAirBuoyancy(grid) -julia> model = NonhydrostaticModel(; grid = RectilinearGrid(size=(8, 8, 8), extent=(1, 2, 3)), - buoyancy, tracers = (:θ, :q)) +# output +MoistAirBuoyancy: +├── reference_state: ReferenceState{Float64}(p₀=101325.0, θᵣ=288.0) +└── thermodynamics: ThermodynamicConstants{Float64} +``` + +To build a model with MoistAirBuoyancy, we include potential temperature and total specific humidity +tracers `θ` and `qᵗ` to the model. + +```jldoctest mab +model = NonhydrostaticModel(; grid, buoyancy, tracers = (:θ, :qᵗ)) + +# output NonhydrostaticModel{CPU, RectilinearGrid}(time = 0 seconds, iteration = 0) -├── grid: 8×8×8 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo +├── grid: 1×1×8 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 1×1×3 halo ├── timestepper: RungeKutta3TimeStepper ├── advection scheme: Centered(order=2) -├── tracers: (θ, q) +├── tracers: (θ, qᵗ) ├── closure: Nothing ├── buoyancy: MoistAirBuoyancy with ĝ = NegativeZDirection() └── coriolis: Nothing ``` """ -function MoistAirBuoyancy(FT=Oceananigans.defaults.FloatType; - thermodynamics = ThermodynamicConstants(FT), - reference_constants = ReferenceStateConstants{FT}(101325, 290)) - - AT = typeof(thermodynamics) - return MoistAirBuoyancy{FT, AT}(reference_constants, thermodynamics) +function MoistAirBuoyancy(grid; + base_pressure = 101325, + reference_potential_temperature = 288, + thermodynamics = ThermodynamicConstants(eltype(grid))) + + reference_state = ReferenceState(grid, thermodynamics; + base_pressure, + potential_temperature = reference_potential_temperature) + + return MoistAirBuoyancy(reference_state, thermodynamics) end Base.summary(b::MoistAirBuoyancy) = "MoistAirBuoyancy" function Base.show(io::IO, b::MoistAirBuoyancy) - print(io, summary(b), "\n", - "├── reference_constants: ", summary(b.reference_constants), "\n", + print(io, summary(b), ":\n", + "├── reference_state: ", summary(b.reference_state), "\n", "└── thermodynamics: ", summary(b.thermodynamics)) end -required_tracers(::MoistAirBuoyancy) = (:θ, :q) -reference_density(z, mb::MoistAirBuoyancy) = reference_density(z, mb.reference_constants, mb.thermodynamics) -base_density(mb::MoistAirBuoyancy) = base_density(mb.reference_constants, mb.thermodynamics) - -##### -##### -##### +required_tracers(::MoistAirBuoyancy) = (:θ, :qᵗ) const c = Center() + @inline function buoyancy_perturbationᶜᶜᶜ(i, j, k, grid, mb::MoistAirBuoyancy, tracers) + @inbounds begin + pᵣ = mb.reference_state.pressure[i, j, k] + ρᵣ = mb.reference_state.density[i, j, k] + θ = tracers.θ[i, j, k] + qᵗ = tracers.qᵗ[i, j, k] + end + z = Oceananigans.Grids.znode(i, j, k, grid, c, c, c) - θ = @inbounds tracers.θ[i, j, k] - q = @inbounds tracers.q[i, j, k] - 𝒰 = HeightReferenceThermodynamicState(θ, q, z) + p₀ = mb.reference_state.base_pressure + q = MoistureMassFractions(qᵗ, zero(qᵗ), zero(qᵗ)) + 𝒰 = PotentialTemperatureState(θ, q, z, p₀, pᵣ, ρᵣ) - # Perform Saturation adjustment - α = specific_volume(𝒰, mb.reference_constants, mb.thermodynamics) + # Perform saturation adjustment + T = temperature(𝒰, mb.thermodynamics) - # Compute reference specific volume - αᵣ = reference_specific_volume(z, mb.reference_constants, mb.thermodynamics) - g = mb.thermodynamics.gravitational_acceleration + # Compute specific volume + Rᵐ = mixture_gas_constant(q, mb.thermodynamics) + α = Rᵐ * T / pᵣ - # Formulation in terms of base density: - # ρ₀ = base_density(mb.reference_constants, mb.thermodynamics) - # return ρ₀ * g * (α - αᵣ) + g = mb.thermodynamics.gravitational_acceleration - return g * (α - αᵣ) / αᵣ + # b = g * (α - αᵣ) / αᵣ + return g * (ρᵣ * α - 1) end @inline ∂z_b(i, j, k, grid, mb::MoistAirBuoyancy, tracers) = ∂zᶜᶜᶠ(i, j, k, grid, buoyancy_perturbationᶜᶜᶜ, mb, tracers) -const c = Center() - -##### -##### Temperature -##### - -function temperature(i, j, k, grid::AbstractGrid, mb::MoistAirBuoyancy, θ, q) - z = Oceananigans.Grids.znode(i, j, k, grid, c, c, c) - θi = @inbounds θ[i, j, k] - qi = @inbounds q[i, j, k] - 𝒰 = HeightReferenceThermodynamicState(θi, qi, z) - return temperature(𝒰, mb.reference_constants, mb.thermodynamics) -end - -struct TemperatureKernelFunction end - -@inline (::TemperatureKernelFunction)(i, j, k, grid, buoyancy, θ, q) = - temperature(i, j, k, grid, buoyancy, θ, q) - -function TemperatureField(model) - func = TemperatureKernelFunction() - grid = model.grid - buoyancy = model.buoyancy.formulation - θ = model.tracers.θ - q = model.tracers.q - op = KernelFunctionOperation{Center, Center, Center}(func, grid, buoyancy, θ, q) - return Field(op) -end - -##### -##### Saturation specific humidity -##### - -@inline function saturation_specific_humidity(i, j, k, grid, mb::MoistAirBuoyancy, T, condensed_phase) - z = Oceananigans.Grids.znode(i, j, k, grid, c, c, c) - Ti = @inbounds T[i, j, k] - return saturation_specific_humidity(Ti, z, mb.reference_constants, mb.thermodynamics, condensed_phase) -end - -struct PhaseTransitionConstantsKernel{T, P} - condensed_phase :: P - temperature :: T -end - -Adapt.adapt_structure(to, sk::PhaseTransitionConstantsKernel) = - PhaseTransitionConstantsKernel(adapt(to, sk.condensed_phase), - adapt(to, sk.temperature)) - -@inline function (kernel::PhaseTransitionConstantsKernel)(i, j, k, grid, buoyancy) - T = kernel.temperature - return saturation_specific_humidity(i, j, k, grid, buoyancy, T, kernel.condensed_phase) -end - -function SaturationField(model, - T = TemperatureField(model); - condensed_phase = model.buoyancy.formulation.thermodynamics.liquid) - func = PhaseTransitionConstantsKernel(condensed_phase, T) - grid = model.grid - buoyancy = model.buoyancy.formulation - op = KernelFunctionOperation{Center, Center, Center}(func, grid, buoyancy) - return Field(op) -end - -##### -##### Condensate -##### - -struct CondensateKernel{T} - temperature :: T -end - -Adapt.adapt_structure(to, ck::CondensateKernel) = CondensateKernel(adapt(to, ck.temperature)) - -@inline function condensate_specific_humidity(i, j, k, grid, mb::MoistAirBuoyancy, T, q) - z = Oceananigans.Grids.znode(i, j, k, grid, c, c, c) - Ti = @inbounds T[i, j, k] - qi = @inbounds q[i, j, k] - qˡ = condensate_specific_humidity(Ti, qi, z, mb.reference_constants, mb.thermodynamics) - return qˡ -end - -@inline function (kernel::CondensateKernel)(i, j, k, grid, buoyancy, q) - T = kernel.temperature - return condensate_specific_humidity(i, j, k, grid, buoyancy, T, q) -end - -function CondensateField(model, T=TemperatureField(model)) - func = CondensateKernel(T) - grid = model.grid - buoyancy = model.buoyancy.formulation - q = model.tracers.q - op = KernelFunctionOperation{Center, Center, Center}(func, grid, buoyancy, q) - return Field(op) -end - ##### ##### Saturation adjustment ##### -# Organizing information about the state is a WIP -struct HeightReferenceThermodynamicState{FT} - θ :: FT - q :: FT - z :: FT -end - -Base.summary(::HeightReferenceThermodynamicState{FT}) where FT = "HeightReferenceThermodynamicState{$FT}" - -function Base.show(io::IO, hrts::HeightReferenceThermodynamicState) - print(io, summary(hrts), ":", '\n', - "├── θ: ", hrts.θ, "\n", - "├── q: ", hrts.q, "\n", - "└── z: ", hrts.z) -end - -condensate_specific_humidity(T, state::HeightReferenceThermodynamicState, ref, thermo) = - condensate_specific_humidity(T, state.q, state.z, ref, thermo) - - # Solve # θ = T/Π ( 1 - ℒ qˡ / (cᵖᵐ T)) # for temperature T with qˡ = max(0, q - qᵛ⁺). # root of: f(T) = T - Π θ - ℒ qˡ / cᵖᵐ """ - temperature(state::HeightReferenceThermodynamicState, ref, thermo) + temperature(state::PotentialTemperatureState, ref, thermo) Return the temperature ``T`` that satisfies saturation adjustment, that is, the temperature for which ```math -θ = [1 - ℒ qˡ / (cᵖᵐ T)] T / Π , +θ = [1 - ℒˡᵣ qˡ / (cᵖᵐ T)] T / Π , ``` -with ``qˡ = \\max(0, qᵗ - qᵛ⁺)`` the condensate specific humidity, where ``qᵗ`` is the +with ``ℒˡᵣ`` the latent heat at the reference temperature ``Tᵣ``, ``cᵖᵐ`` the mixture +specific heat, ``Π`` the Exner function, ``qˡ = \\max(0, qᵗ - qᵛ⁺)`` +the condensate specific humidity, ``qᵗ`` is the total specific humidity, ``qᵛ⁺`` is the saturation specific humidity. The saturation adjustment temperature is obtained by solving ``r(T)``, where ```math -r(T) ≡ T - θ Π - ℒ qˡ / (cᵖᵐ T) . +r(T) ≡ T - θ Π - ℒˡᵣ qˡ / cᵖᵐ . ``` Solution of ``r(T) = 0`` is found via the [secant method](https://en.wikipedia.org/wiki/Secant_method). """ -@inline function temperature(state::HeightReferenceThermodynamicState{FT}, ref, thermo) where FT - state.θ == 0 && return zero(FT) - - # Generate guess for unsaturated conditions - Π = exner_function(state, ref, thermo) - T₁ = Π * state.θ - qˡ₁ = condensate_specific_humidity(T₁, state, ref, thermo) - qˡ₁ <= 0 && return T₁ - - # If we made it this far, we have condensation - r₁ = saturation_adjustment_residual(T₁, Π, qˡ₁, state, thermo) - - ℒᵛ = thermo.liquid.latent_heat - cᵖᵐ = mixture_heat_capacity(state.q, thermo) - T₂ = T₁ + ℒᵛ * qˡ₁ / cᵖᵐ - qˡ₂ = condensate_specific_humidity(T₂, state, ref, thermo) - r₂ = saturation_adjustment_residual(T₂, Π, qˡ₂, state, thermo) - - # Saturation adjustment - R = sqrt(max(T₂, T₁)) - ϵ = convert(FT, 1e-9) - δ = ϵ * R +@inline function temperature(𝒰₀::PotentialTemperatureState{FT}, thermo) where FT + θ = 𝒰₀.potential_temperature + θ == 0 && return zero(FT) + + # Generate guess for unsaturated conditions; if dry, return T₁ + qᵗ = total_specific_humidity(𝒰₀) + q₁ = MoistureMassFractions(qᵗ, zero(qᵗ), zero(qᵗ)) + 𝒰₁ = with_moisture(𝒰₀, q₁) + Π₁ = exner_function(𝒰₀, thermo) + T₁ = Π₁ * θ + + pᵣ = 𝒰₀.reference_pressure + ρ₁ = density(pᵣ, T₁, q₁, thermo) + qᵛ⁺₁ = saturation_specific_humidity(T₁, ρ₁, thermo, thermo.liquid) + qᵗ <= qᵛ⁺₁ && return T₁ + + # If we made it this far, the state is saturated. + # T₁ then provides a lower bound, and our state 𝒰₁ + # has to be modified to consistently include the liquid mass fraction. + # Subsequent computations will assume that the specific humidity + # is given by the saturation specific humidity, eg ``qᵛ = qᵛ⁺``. + qᵛ⁺₁ = adjustment_saturation_specific_humidity(T₁, 𝒰₁, thermo) + qˡ₁ = qᵗ - qᵛ⁺₁ + q₁ = MoistureMassFractions(qᵛ⁺₁, qˡ₁, zero(qˡ₁)) + 𝒰₁ = with_moisture(𝒰₀, q₁) + + # We generate a second guess simply by adding 1 K to T₁... + + # NOTE: We could also generate a second guess to start a secant iteration + # by applying the potential temperature assuming a liquid fraction + # associated with T₁. This should represent an _overestimate_, + # since ``qᵛ⁺₁(T₁)`` underestimates the saturation specific humidity, + # and therefore qˡ₁ is overestimated. This is similar to an approach + # used in Pressel et al 2015. However, it doesn't work for large liquid fractions. + T₂ = T₁ + 1 + + #= + ℒˡᵣ = thermo.liquid.reference_latent_heat + cᵖᵐ = mixture_heat_capacity(q₁, thermo) + T₂ = T₁ + ℒˡᵣ * qˡ₁ / cᵖᵐ + =# + + 𝒰₂ = adjust_state(𝒰₁, T₂, thermo) + + # Initialize saturation adjustment + r₁ = saturation_adjustment_residual(T₁, 𝒰₁, thermo) + r₂ = saturation_adjustment_residual(T₂, 𝒰₂, thermo) + δ = convert(FT, 1e-3) iter = 0 - while abs(r₂ - r₁) > δ + while abs(T₂ - T₁) > δ # Compute slope ΔTΔr = (T₂ - T₁) / (r₂ - r₁) # Store previous values r₁ = r₂ T₁ = T₂ + 𝒰₁ = 𝒰₂ - # Update T₂ -= r₂ * ΔTΔr - qˡ₂ = condensate_specific_humidity(T₂, state, ref, thermo) - r₂ = saturation_adjustment_residual(T₂, Π, qˡ₂, state, thermo) + 𝒰₂ = adjust_state(𝒰₂, T₂, thermo) + r₂ = saturation_adjustment_residual(T₂, 𝒰₂, thermo) + iter += 1 end return T₂ end -@inline function saturation_adjustment_residual(T, Π, qˡ, state::HeightReferenceThermodynamicState, thermo) - ℒᵛ = thermo.liquid.latent_heat - cᵖᵐ = mixture_heat_capacity(state.q, thermo) - return T - ℒᵛ * qˡ / cᵖᵐ - Π * state.θ +# This estimate assumes that the specific humidity is itself the saturation +# specific humidity, eg ``qᵛ = qᵛ⁺``. Knowledge of the specific humidity +# is needed to compute the mixture gas constant, and thus density, +# which in turn is needed to compute the _saturation_ specific humidity. +# This consideration culminates in a new expression for saturation specific humidity +# used below, and also written in Pressel et al 2015, equation 37. +# (There is an error in the description below it, but the equation 37 is correct.) +@inline function adjustment_saturation_specific_humidity(T, 𝒰, thermo) + pᵛ⁺ = saturation_vapor_pressure(T, thermo, thermo.liquid) + pᵣ = 𝒰.reference_pressure + qᵗ = total_specific_humidity(𝒰) + Rᵈ = dry_air_gas_constant(thermo) + Rᵛ = vapor_gas_constant(thermo) + ϵᵈᵛ = Rᵈ / Rᵛ + return ϵᵈᵛ * (1 - qᵗ) * pᵛ⁺ / (pᵣ - pᵛ⁺) end -@inline function specific_volume(state::HeightReferenceThermodynamicState, ref, thermo) - T = temperature(state, ref, thermo) - Rᵐ = mixture_gas_constant(state.q, thermo) - pᵣ = reference_pressure(state.z, ref, thermo) - return Rᵐ * T / pᵣ +@inline function adjust_state(𝒰₀, T, thermo) + qᵛ⁺ = adjustment_saturation_specific_humidity(T, 𝒰₀, thermo) + qᵗ = total_specific_humidity(𝒰₀) + qˡ = max(0, qᵗ - qᵛ⁺) + q₁ = MoistureMassFractions(qᵛ⁺, qˡ, zero(qˡ)) + return with_moisture(𝒰₀, q₁) end -@inline function exner_function(state::HeightReferenceThermodynamicState, ref, thermo) - Rᵐ = mixture_gas_constant(state.q, thermo) - cᵖᵐ = mixture_heat_capacity(state.q, thermo) - inv_ϰᵐ = Rᵐ / cᵖᵐ - pᵣ = reference_pressure(state.z, ref, thermo) - p₀ = ref.base_pressure - return (pᵣ / p₀)^inv_ϰᵐ +@inline function saturation_adjustment_residual(T, 𝒰, thermo) + Π = exner_function(𝒰, thermo) + q = 𝒰.moisture_fractions + θ = 𝒰.potential_temperature + ℒˡᵣ = thermo.liquid.reference_latent_heat + cᵖᵐ = mixture_heat_capacity(q, thermo) + qˡ = q.liquid + θ = 𝒰.potential_temperature + return T - Π * θ - ℒˡᵣ * qˡ / cᵖᵐ end ##### -##### Reference implementation of an "unsaturated" moist air buoyancy model, -##### which assumes unsaturated air +##### Diagnostics ##### -struct UnsaturatedMoistAirBuoyancy{FT} <: AbstractBuoyancyFormulation{Nothing} - expansion_coefficient :: FT - reference_potential_temperature :: FT - gas_constant_ratio :: FT +const c = Center() + +# Temperature +@inline function temperature(i, j, k, grid::AbstractGrid, mb::MoistAirBuoyancy, θ, qᵗ) + @inbounds begin + θᵢ = θ[i, j, k] + qᵗᵢ = qᵗ[i, j, k] + pᵣ = mb.reference_state.pressure[i, j, k] + ρᵣ = mb.reference_state.density[i, j, k] + end + z = Oceananigans.Grids.znode(i, j, k, grid, c, c, c) + p₀ = mb.reference_state.base_pressure + q = MoistureMassFractions(qᵗᵢ, zero(qᵗᵢ), zero(qᵗᵢ)) + 𝒰 = PotentialTemperatureState(θᵢ, q, z, p₀, pᵣ, ρᵣ) + return temperature(𝒰, mb.thermodynamics) end -function UnsaturatedMoistAirBuoyancy(FT=Oceananigans.defaults.FloatType; - expansion_coefficient = 3.27e-2, - reference_potential_temperature = 0, - gas_constant_ratio = 1.61) +struct TemperatureKernelFunction end + +@inline (::TemperatureKernelFunction)(i, j, k, grid, buoyancy, θ, qᵗ) = + temperature(i, j, k, grid, buoyancy, θ, qᵗ) + +function TemperatureField(model) + func = TemperatureKernelFunction() + grid = model.grid + buoyancy = model.buoyancy.formulation + θ = model.tracers.θ + qᵗ = model.tracers.qᵗ + op = KernelFunctionOperation{Center, Center, Center}(func, grid, buoyancy, θ, qᵗ) + return Field(op) +end - return UnsaturatedMoistAirBuoyancy{FT}(expansion_coefficient, - reference_potential_temperature, - gas_constant_ratio) +# Saturation specific humidity +@inline function saturation_specific_humidity(i, j, k, grid, mb::MoistAirBuoyancy, T, qᵗ, phase) + @inbounds begin + Tᵢ = T[i, j, k] + qᵗᵢ = qᵗ[i, j, k] + pᵣ = mb.reference_state.pressure[i, j, k] + end + q = MoistureMassFractions(qᵗᵢ, zero(qᵗᵢ), zero(qᵗᵢ)) + ρ = density(pᵣ, Tᵢ, q, mb.thermodynamics) + return saturation_specific_humidity(Tᵢ, ρ, mb.thermodynamics, phase) end -required_tracers(::UnsaturatedMoistAirBuoyancy) = (:θ, :q) - -@inline function buoyancy_perturbationᶜᶜᶜ(i, j, k, grid, mb::UnsaturatedMoistAirBuoyancy, tracers) - β = mb.expansion_coefficient - θ₀ = mb.reference_potential_temperature - ϵᵥ = mb.gas_constant_ratio - δ = ϵᵥ - 1 - θ = @inbounds tracers.θ[i, j, k] - q = @inbounds tracers.q[i, j, k] - θᵥ = θ * (1 + δ * q) - return β * (θᵥ - θ₀) +struct PhaseTransitionConstantsKernel{T, P} + condensed_phase :: P + temperature :: T +end + +Adapt.adapt_structure(to, sk::PhaseTransitionConstantsKernel) = + PhaseTransitionConstantsKernel(adapt(to, sk.condensed_phase), + adapt(to, sk.temperature)) + +@inline function (kernel::PhaseTransitionConstantsKernel)(i, j, k, grid, buoyancy, qᵗ) + T = kernel.temperature + return saturation_specific_humidity(i, j, k, grid, buoyancy, T, qᵗ, kernel.condensed_phase) +end + +function SaturationField(model, T = TemperatureField(model); + condensed_phase = model.buoyancy.formulation.thermodynamics.liquid) + func = PhaseTransitionConstantsKernel(condensed_phase, T) + grid = model.grid + buoyancy = model.buoyancy.formulation + qᵗ = model.tracers.qᵗ + op = KernelFunctionOperation{Center, Center, Center}(func, grid, buoyancy, qᵗ) + return Field(op) +end + +# Condensate +struct CondensateKernel{T} + temperature :: T +end + +Adapt.adapt_structure(to, ck::CondensateKernel) = CondensateKernel(adapt(to, ck.temperature)) + +@inline function liquid_mass_fraction(i, j, k, grid, mb::MoistAirBuoyancy, T, θ, qᵗ) + @inbounds begin + Tᵢ = T[i, j, k] + θᵢ = θ[i, j, k] + qᵗᵢ = qᵗ[i, j, k] + pᵣ = mb.reference_state.pressure[i, j, k] + ρᵣ = mb.reference_state.density[i, j, k] + end + + # First assume non-saturation. + z = Oceananigans.Grids.znode(i, j, k, grid, c, c, c) + p₀ = mb.reference_state.base_pressure + q = MoistureMassFractions(qᵗᵢ, zero(qᵗᵢ), zero(qᵗᵢ)) + 𝒰 = PotentialTemperatureState(Tᵢ, q, z, p₀, pᵣ, ρᵣ) + Π = exner_function(𝒰, mb.thermodynamics) + Tᵢ <= Π * θᵢ + 10 * eps(Tᵢ) && return zero(qᵗᵢ) + + # Next assume a saturation value + qᵛ⁺ = adjustment_saturation_specific_humidity(Tᵢ, 𝒰, mb.thermodynamics) + return max(0, qᵗᵢ - qᵛ⁺) +end + +@inline function (kernel::CondensateKernel)(i, j, k, grid, buoyancy, θ, qᵗ) + T = kernel.temperature + return liquid_mass_fraction(i, j, k, grid, buoyancy, T, θ, qᵗ) +end + +function CondensateField(model, T=TemperatureField(model)) + func = CondensateKernel(T) + grid = model.grid + buoyancy = model.buoyancy.formulation + qᵗ = model.tracers.qᵗ + θ = model.tracers.θ + op = KernelFunctionOperation{Center, Center, Center}(func, grid, buoyancy, θ, qᵗ) + return Field(op) end end # module diff --git a/src/Thermodynamics/Thermodynamics.jl b/src/Thermodynamics/Thermodynamics.jl index 9950eaf6..783f98ce 100644 --- a/src/Thermodynamics/Thermodynamics.jl +++ b/src/Thermodynamics/Thermodynamics.jl @@ -1,11 +1,12 @@ module Thermodynamics -export ThermodynamicConstants, ReferenceStateConstants, IdealGas, +export ThermodynamicConstants, ReferenceState, IdealGas, CondensedPhase, mixture_gas_constant, mixture_heat_capacity -include("atmosphere_thermodynamics.jl") +include("thermodynamics_constants.jl") include("vapor_saturation.jl") include("reference_states.jl") +include("dynamic_states.jl") end # module diff --git a/src/Thermodynamics/dynamic_states.jl b/src/Thermodynamics/dynamic_states.jl new file mode 100644 index 00000000..b1456564 --- /dev/null +++ b/src/Thermodynamics/dynamic_states.jl @@ -0,0 +1,41 @@ +struct PotentialTemperatureState{FT} + potential_temperature :: FT + moisture_fractions :: MoistureMassFractions{FT} + height :: FT + base_pressure :: FT + reference_pressure :: FT + reference_density :: FT +end + +@inline function exner_function(𝒰::PotentialTemperatureState, thermo::ThermodynamicConstants) + q = 𝒰.moisture_fractions + z = 𝒰.height + Rᵐ = mixture_gas_constant(q, thermo) + cᵖᵐ = mixture_heat_capacity(q, thermo) + pᵣ = 𝒰.reference_pressure + p₀ = 𝒰.base_pressure + return (pᵣ / p₀)^(Rᵐ / cᵖᵐ) +end + +@inline total_specific_humidity(state::PotentialTemperatureState) = + total_specific_humidity(state.moisture_fractions) + +@inline function with_moisture(𝒰::PotentialTemperatureState, q::MoistureMassFractions) + return PotentialTemperatureState(𝒰.potential_temperature, + q, + 𝒰.height, + 𝒰.base_pressure, + 𝒰.reference_pressure, + 𝒰.reference_density) +end + +# TODO: deprecate this +struct AnelasticThermodynamicState{FT} + potential_temperature :: FT + moisture_fractions :: MoistureMassFractions{FT} + reference_density :: FT + reference_pressure :: FT + exner_function :: FT +end + +@inline total_specific_humidity(state::AnelasticThermodynamicState) = total_specific_humidity(state.moisture_fractions) diff --git a/src/Thermodynamics/reference_states.jl b/src/Thermodynamics/reference_states.jl index 5856740b..cd7eb30a 100644 --- a/src/Thermodynamics/reference_states.jl +++ b/src/Thermodynamics/reference_states.jl @@ -1,72 +1,78 @@ -using Oceananigans: Oceananigans - +using Oceananigans: Oceananigans, Center, Field, set!, fill_halo_regions! + ##### ##### Reference state computations for Boussinesq and Anelastic models ##### -struct ReferenceStateConstants{FT} +struct ReferenceState{FT, F} base_pressure :: FT # base pressure: reference pressure at z=0 - reference_potential_temperature :: FT # constant reference potential temperature + potential_temperature :: FT # constant reference potential temperature + pressure :: F + density :: F end -function ReferenceStateConstants(FT = Oceananigans.defaults.FloatType; - base_pressure = 101325, - potential_temperature = 288) +Base.eltype(::ReferenceState{FT}) where FT = FT - return ReferenceStateConstants{FT}(convert(FT, base_pressure), - convert(FT, potential_temperature)) +function Base.summary(ref::ReferenceState) + FT = eltype(ref) + return string("ReferenceState{$FT}(p₀=", prettysummary(ref.base_pressure), + ", θᵣ=", prettysummary(ref.potential_temperature), ")") end -""" - reference_density(z, ref::ReferenceStateConstants, thermo) +Base.show(io::IO, ref::ReferenceState) = print(io, summary(ref)) -Compute the reference density at height `z` that associated with the reference pressure and -potential temperature. The reference density is defined as the density of dry air at the -reference pressure and temperature. -""" -@inline function reference_density(z, ref::ReferenceStateConstants, thermo) - Rᵈ = dry_air_gas_constant(thermo) - cᵖᵈ = thermo.dry_air.heat_capacity - pᵣ = reference_pressure(z, ref, thermo) - ρ₀ = base_density(ref, thermo) - p₀ = ref.base_pressure - return ρ₀ * (pᵣ / p₀)^(1 - Rᵈ / cᵖᵈ) -end +##### +##### How to compute the reference state +##### -@inline function base_density(ref::ReferenceStateConstants, thermo) +@inline function base_density(p₀, θᵣ, thermo) Rᵈ = dry_air_gas_constant(thermo) - p₀ = ref.base_pressure - θᵣ = ref.reference_potential_temperature return p₀ / (Rᵈ * θᵣ) end -@inline function reference_specific_volume(z, ref::ReferenceStateConstants, thermo) - Rᵈ = dry_air_gas_constant(thermo) - pᵣ = reference_pressure(z, ref, thermo) - θᵣ = ref.reference_potential_temperature - return Rᵈ * θᵣ / pᵣ -end +""" + adiabatic_hydrostatic_pressure(z, p₀, θᵣ, thermo) -@inline function reference_pressure(z, ref::ReferenceStateConstants, thermo) +Compute the reference pressure at height `z` that associated with the reference pressure and +potential temperature. The reference pressure is defined as the pressure of dry air at the +reference pressure and temperature. +""" +@inline function adiabatic_hydrostatic_pressure(z, p₀, θᵣ, thermo) cᵖᵈ = thermo.dry_air.heat_capacity Rᵈ = dry_air_gas_constant(thermo) g = thermo.gravitational_acceleration - θᵣ = ref.reference_potential_temperature - p₀ = ref.base_pressure return p₀ * (1 - g * z / (cᵖᵈ * θᵣ))^(cᵖᵈ / Rᵈ) end -@inline function saturation_specific_humidity(T, z, ref::ReferenceStateConstants, thermo, condensed_phase) - ρ = reference_density(z, ref, thermo) - return saturation_specific_humidity(T, ρ, thermo, condensed_phase) -end +""" + adiabatic_hydrostatic_density(z, p₀, θᵣ, thermo) -function condensate_specific_humidity(T, q, z, ref::ReferenceStateConstants, thermo) - qᵛ⁺ = saturation_specific_humidity(T, z, ref, thermo, thermo.liquid) - return max(0, q - qᵛ⁺) +Compute the reference density at height `z` that associated with the reference pressure and +potential temperature. The reference density is defined as the density of dry air at the +reference pressure and temperature. +""" +@inline function adiabatic_hydrostatic_density(z, p₀, θᵣ, thermo) + Rᵈ = dry_air_gas_constant(thermo) + cᵖᵈ = thermo.dry_air.heat_capacity + pᵣ = adiabatic_hydrostatic_pressure(z, p₀, θᵣ, thermo) + ρ₀ = base_density(p₀, θᵣ, thermo) + return ρ₀ * (pᵣ / p₀)^(1 - Rᵈ / cᵖᵈ) end -function ice_specific_humidity(T, q, z, ref::ReferenceStateConstants, thermo) - qi⁺ = saturation_specific_humidity(T, z, ref, thermo, thermo.solid) - return max(0, q - qi⁺) +function ReferenceState(grid, thermo; + base_pressure = 101325, + potential_temperature = 288) + + FT = eltype(grid) + p₀ = convert(FT, base_pressure) + θᵣ = convert(FT, potential_temperature) + + pᵣ = Field{Nothing, Nothing, Center}(grid) + ρᵣ = Field{Nothing, Nothing, Center}(grid) + set!(pᵣ, z -> adiabatic_hydrostatic_pressure(z, p₀, θᵣ, thermo)) + set!(ρᵣ, z -> adiabatic_hydrostatic_density(z, p₀, θᵣ, thermo)) + fill_halo_regions!(pᵣ) + fill_halo_regions!(ρᵣ) + + return ReferenceState(p₀, θᵣ, pᵣ, ρᵣ) end diff --git a/src/Thermodynamics/atmosphere_thermodynamics.jl b/src/Thermodynamics/thermodynamics_constants.jl similarity index 70% rename from src/Thermodynamics/atmosphere_thermodynamics.jl rename to src/Thermodynamics/thermodynamics_constants.jl index 8d9c8f91..050172d8 100644 --- a/src/Thermodynamics/atmosphere_thermodynamics.jl +++ b/src/Thermodynamics/thermodynamics_constants.jl @@ -48,48 +48,48 @@ function IdealGas(FT = Oceananigans.defaults.FloatType; end struct CondensedPhase{FT} - latent_heat :: FT + reference_latent_heat :: FT heat_capacity :: FT end function Base.summary(ph::CondensedPhase{FT}) where FT return string("CondensedPhase{", FT, "}(", - "latent_heat=", prettysummary(ph.latent_heat), ", ", + "reference_latent_heat=", prettysummary(ph.reference_latent_heat), ", ", "heat_capacity=", prettysummary(ph.heat_capacity), ")") end Base.show(io::IO, ph::CondensedPhase) = print(io, summary(ph)) Adapt.adapt_structure(to, pt::CondensedPhase) = - CondensedPhase(adapt(to, pt.latent_heat), + CondensedPhase(adapt(to, pt.reference_latent_heat), adapt(to, pt.heat_capacity)) """ - CondensedPhase(FT = Oceananigans.defaults.FloatType; latent_heat, heat_capacity) + CondensedPhase(FT = Oceananigans.defaults.FloatType; reference_latent_heat, heat_capacity) Returns `CondensedPhase` with specified parameters converted to `FT`. Two examples of `CondensedPhase` are liquid and solid. When matter is converted from vapor to liquid, water molecules in the gas phase cluster together and slow down to form liquid with `heat_capacity`, -The lost of molecular kinetic energy is called the `latent_heat`. +The lost of molecular kinetic energy is called the `reference_latent_heat`. Likewise, during deposition, water molecules in the gas phase cluster into ice crystals. Arguments ========= - `FT`: Float type to use (defaults to Oceananigans.defaults.FloatType) -- `latent_heat`: Difference between the internal energy of the gaseous phase at +- `reference_latent_heat`: Difference between the internal energy of the gaseous phase at the `energy_reference_temperature`. - `heat_capacity`: Heat capacity of the phase of matter. """ -function CondensedPhase(FT = Oceananigans.defaults.FloatType; latent_heat, heat_capacity) - return CondensedPhase{FT}(convert(FT, latent_heat), +function CondensedPhase(FT = Oceananigans.defaults.FloatType; reference_latent_heat, heat_capacity) + return CondensedPhase{FT}(convert(FT, reference_latent_heat), convert(FT, heat_capacity)) end -liquid_water(FT) = CondensedPhase(FT; latent_heat=2500800, heat_capacity=4181) -water_ice(FT) = CondensedPhase(FT; latent_heat=2834000, heat_capacity=2108) +liquid_water(FT) = CondensedPhase(FT; reference_latent_heat=2500800, heat_capacity=4181) +water_ice(FT) = CondensedPhase(FT; reference_latent_heat=2834000, heat_capacity=2108) struct ThermodynamicConstants{FT, C, S} molar_gas_constant :: FT @@ -148,7 +148,7 @@ end ThermodynamicConstants(FT = Oceananigans.defaults.FloatType; molar_gas_constant = 8.314462618, gravitational_acceleration = 9.81, - energy_reference_temperature = 273.16, + energy_reference_temperature = 273.15, triple_point_temperature = 273.16, triple_point_pressure = 611.657, dry_air_molar_mass = 0.02897, @@ -200,7 +200,7 @@ can then be used for both condensation (vapor → liquid) and deposition (vapor function ThermodynamicConstants(FT = Oceananigans.defaults.FloatType; molar_gas_constant = 8.314462618, gravitational_acceleration = 9.81, - energy_reference_temperature = 273.16, + energy_reference_temperature = 273.15, triple_point_temperature = 273.16, triple_point_pressure = 611.657, dry_air_molar_mass = 0.02897, @@ -228,134 +228,85 @@ function ThermodynamicConstants(FT = Oceananigans.defaults.FloatType; end const TC = ThermodynamicConstants -const IG = IdealGas @inline vapor_gas_constant(thermo::TC) = thermo.molar_gas_constant / thermo.vapor.molar_mass @inline dry_air_gas_constant(thermo::TC) = thermo.molar_gas_constant / thermo.dry_air.molar_mass -const NonCondensingThermodynamicConstants{FT} = ThermodynamicConstants{FT, Nothing, Nothing} +##### +##### Mixtures of dry air with vapor, liquid, and ice +##### + +struct MoistureMassFractions{FT} + vapor :: FT + liquid :: FT + ice :: FT +end + +const MMF = MoistureMassFractions + +@inline total_specific_humidity(q::MMF) = q.vapor + q.liquid + q.ice +@inline dry_air_mass_fraction(q::MMF) = 1 - total_specific_humidity(q) """ - mixture_gas_constant(q, thermo) + mixture_gas_constant(q::MoistureMassFractions, thermo) -Compute the gas constant of moist air given the specific humidity `q` and +Compute the gas constant of moist air given the specific humidity `q` and thermodynamic parameters `thermo`. The mixture gas constant is calculated as a weighted average of the dry air -and water vapor gas constants: +and water vapor gas thermo: ```math -R_m = R_d (1 - q) + R_v q +Rᵐ = qᵈ * Rᵈ + qᵛ * Rᵛ ``` where: -- `R_d` is the dry air gas constant -- `R_v` is the water vapor gas constant -- `q` is the specific humidity (mass fraction of water vapor) +- `Rᵈ` is the dry air gas constant +- `Rᵛ` is the water vapor gas constant +- `qᵈ` is the mass fraction of dry air +- `qᵛ` is the mass fraction of water vapor # Arguments -- `q`: Specific humidity (dimensionless) -- `thermo`: `ThermodynamicConstants` instance containing gas constants +- `qᵈ`: Mass fraction of dry air (dimensionless) +- `qᵛ`: Mass fraction of water vapor (dimensionless) +- `thermo`: `ThermodynamicConstants` instance containing gas thermo # Returns - Gas constant of the moist air mixture in J/(kg·K) """ -@inline function mixture_gas_constant(q, thermo::TC) +@inline function mixture_gas_constant(q::MMF, thermo::TC) + qᵈ = dry_air_mass_fraction(q) + qᵛ = q.vapor Rᵈ = dry_air_gas_constant(thermo) Rᵛ = vapor_gas_constant(thermo) - return Rᵈ * (1 - q) + Rᵛ * q + return qᵈ * Rᵈ + qᵛ * Rᵛ end """ - mixture_heat_capacity(q, thermo) + mixture_heat_capacity(qᵈ, qᵛ, thermo) Compute the heat capacity of state air given the total specific humidity q and assuming that condensate mass ratio qᶜ ≪ q, where qℓ is the mass ratio of liquid condensate. """ -@inline function mixture_heat_capacity(q, thermo::TC) +@inline function mixture_heat_capacity(q::MMF, thermo::TC) + qᵈ = dry_air_mass_fraction(q) + qᵛ = q.vapor cᵖᵈ = thermo.dry_air.heat_capacity cᵖᵛ = thermo.vapor.heat_capacity - return cᵖᵈ * (1 - q) + cᵖᵛ * q + return qᵈ * cᵖᵈ + qᵛ * cᵖᵛ end ##### -##### state thermodynamics for a Boussinesq model +##### Equation of state ##### -# Organizing information about the state is a WIP -struct ThermodynamicState{FT} - θ :: FT - q :: FT - z :: FT -end - -struct ReferenceState{FT} - p₀ :: FT # base pressure: reference pressure at z=0 - θ :: FT # constant reference potential temperature +@inline function density(p, T, q::MMF, thermo::TC) + Rᵐ = mixture_gas_constant(q, thermo) + return p / (Rᵐ * T) end -Adapt.adapt_structure(to, ref::ReferenceState) = - ReferenceState(adapt(to, ref.p₀), - adapt(to, ref.θ)) - -function ReferenceState(FT = Oceananigans.defaults.FloatType; - base_pressure = 101325, - potential_temperature = 288) - - return ReferenceState{FT}(convert(FT, base_pressure), - convert(FT, potential_temperature)) -end - -""" - reference_density(z, ref, thermo) - -Compute the reference density associated with the reference pressure and potential temperature. -The reference density is defined as the density of dry air at the reference pressure and temperature. -""" -@inline function reference_density(z, ref, thermo) - Rᵈ = dry_air_gas_constant(thermo) - p = reference_pressure(z, ref, thermo) - return p / (Rᵈ * ref.θ) -end - -@inline function base_density(ref, thermo) - Rᵈ = dry_air_gas_constant(thermo) - return ref.p₀ / (Rᵈ * ref.θ) -end - -@inline function reference_specific_volume(z, ref, thermo) - Rᵈ = dry_air_gas_constant(thermo) - p = reference_pressure(z, ref, thermo) - return Rᵈ * ref.θ / p -end - -@inline function reference_pressure(z, ref, thermo) - cᵖᵈ = thermo.dry_air.heat_capacity - Rᵈ = dry_air_gas_constant(thermo) - inv_ϰᵈ = Rᵈ / cᵖᵈ - g = thermo.gravitational_acceleration - return ref.p₀ * (1 - g * z / (cᵖᵈ * ref.θ))^inv_ϰᵈ -end - -@inline function saturation_specific_humidity(T, z, ref::ReferenceState, thermo, condensed_phase) - ρ = reference_density(z, ref, thermo) - return saturation_specific_humidity(T, ρ, thermo, condensed_phase) -end - -@inline function exner_function(state, ref, thermo) - Rᵐ = mixture_gas_constant(state.q, thermo) - cᵖᵐ = mixture_heat_capacity(state.q, thermo) - inv_ϰᵐ = Rᵐ / cᵖᵐ - pᵣ = reference_pressure(state.z, ref, thermo) - p₀ = ref.base_pressure - return (pᵣ / p₀)^inv_ϰᵐ -end - -condensate_specific_humidity(T, state, ref, thermo) = - condensate_specific_humidity(T, state.q, state.z, ref, thermo) - -function condensate_specific_humidity(T, q, z, ref, thermo) - qᵛ⁺ = saturation_specific_humidity(T, z, ref, thermo, thermo.liquid) - return max(0, q - qᵛ⁺) +@inline function specific_volume(p, T, q::MMF, thermo::TC) + Rᵐ = mixture_gas_constant(q, mb.thermodynamics) + return Rᵐ * T / p end diff --git a/src/Thermodynamics/vapor_saturation.jl b/src/Thermodynamics/vapor_saturation.jl index 6eb735f2..874fae78 100644 --- a/src/Thermodynamics/vapor_saturation.jl +++ b/src/Thermodynamics/vapor_saturation.jl @@ -1,54 +1,62 @@ """ saturation_vapor_pressure(T, thermo, phase::CondensedPhase) -Compute the saturation vapor pressure ``pᵛ⁺`` over a fluid surface at -`phase` (i.e., liquid or solid) from the Clausius-Clapeyron relation, +Compute the saturation vapor pressure ``pᵛ⁺`` over a planar surface +composed of the "``β``-phase" (liquid, or ice) +using the Clausius-Clapeyron relation, ```math -dpᵛ⁺ / dT = pᵛ⁺ ℒᵛ(T) / (Rᵛ T^2) , +dpᵛ⁺ / dT = pᵛ⁺ ℒᵝ(T) / (Rᵛ T^2) , ``` -where ``ℒˡ(T) = ℒˡ(T=0) + Δcˡ T``, with ``Δcˡ ≡ (cᵖᵛ - cᵖˡ)`` . +where the latent heat is ``ℒᵝ(T) = ℒᵝ(T=0) + Δcᵝ T``, with ``Δcᵝ ≡ (cᵖᵛ - cᵝ)`` . The saturation vapor pressure ``pᵛ⁺`` is obtained after integrating the above from the triple point, i.e., ``p(Tᵗʳ) = pᵗʳ`` to get ```math -pᵛ⁺(T) = pᵗʳ \\left ( \\frac{T}{Tᵗʳ} \\right )^{Δcˡ / Rᵛ} \\exp \\left [ (1/Tᵗʳ - 1/T) ℒˡ(T=0) / Rᵛ \\right ] . +pᵛ⁺(T) = pᵗʳ \\left ( \\frac{T}{Tᵗʳ} \\right )^{Δcᵝ / Rᵛ} \\exp \\left [ (1/Tᵗʳ - 1/T) ℒᵝ(T=0) / Rᵛ \\right ] . ``` -Note that latent heat ``ℒ₀`` is at the reference temperature ``T₀`` -and that ``ℒ(T=0) = ℒ₀ - Δcˡ T₀``. +Note that latent heat ``ℒᵝ(T=0)`` is the difference between the enthalpy of water vapor +and the phase ``β`` when the temperature is ``T = 0``K, or absolute zero. +We define the latent heat using its value ``ℒᵝᵣ = ℒᵝ(T=Tᵣ)`` at the "energy reference temperature" +``T = Tᵣ`` (usually ``Tᵣ ≡ 273.15``K ``= 0^∘``C), such that + +```math +ℒᵝ(T) = ℒᵝᵣ + Δcᵝ (T - Tᵣ), \\quad \text{and} \\quad ℒᵝ(T=0) = ℒᵝᵣ - Δcᵝ Tᵣ`` . +``` """ @inline function saturation_vapor_pressure(T, thermo, phase::CondensedPhase) - ℒ₀ = phase.latent_heat # at thermo.energy_reference_temperature - cᵖˡ = phase.heat_capacity - T₀ = thermo.energy_reference_temperature + ℒᵣ = phase.reference_latent_heat # at thermo.energy_reference_temperature + Tᵣ = thermo.energy_reference_temperature + Tᵗʳ = thermo.triple_point_temperature pᵗʳ = thermo.triple_point_pressure - cᵖᵛ = thermo.vapor.heat_capacity Rᵛ = vapor_gas_constant(thermo) - Δcˡ = cᵖᵛ - cᵖˡ - return pᵗʳ * (T / Tᵗʳ)^(Δcˡ / Rᵛ) * exp((1/Tᵗʳ - 1/T) * (ℒ₀ - Δcˡ * T₀) / Rᵛ) + cᵖᵛ = thermo.vapor.heat_capacity + cᵝ = phase.heat_capacity + Δcᵝ = cᵖᵛ - cᵝ + + return pᵗʳ * (T / Tᵗʳ)^(Δcᵝ / Rᵛ) * exp((1/Tᵗʳ - 1/T) * (ℒᵣ - Δcᵝ * Tᵣ) / Rᵛ) end -# Over a liquid surface """ - saturation_specific_humidity(T, ρ, thermo, condensed_phase::CondensedPhase) + saturation_specific_humidity(T, ρ, thermo, phase::CondensedPhase) Compute the saturation specific humidity for a gas at temperature `T`, total -density `ρ`, `thermo`dynamics, and `condensed_phase` via: +density `ρ`, `thermo`dynamics, and `phase` via: ```math qᵛ⁺ = pᵛ⁺ / (ρ Rᵛ T) , ``` -where ``pᵛ⁺`` is the [`saturation_vapor_pressure`](@ref), and ``Rᵛ`` is the specific gas -constant for water vapor. +where ``pᵛ⁺`` is the [`saturation_vapor_pressure`](@ref), ``ρ`` is total density, +and ``Rᵛ`` is the specific gas constant for water vapor. """ -@inline function saturation_specific_humidity(T, ρ, thermo, condensed_phase::CondensedPhase) - p★ = saturation_vapor_pressure(T, thermo, condensed_phase) +@inline function saturation_specific_humidity(T, ρ, thermo::ThermodynamicConstants, phase::CondensedPhase) + pᵛ⁺ = saturation_vapor_pressure(T, thermo, phase) Rᵛ = vapor_gas_constant(thermo) - return p★ / (ρ * Rᵛ * T) + return pᵛ⁺ / (ρ * Rᵛ * T) end diff --git a/src/microphysics.jl b/src/microphysics.jl index f1f1aa0c..e69de29b 100644 --- a/src/microphysics.jl +++ b/src/microphysics.jl @@ -1,62 +0,0 @@ -##### -##### Microphysics schemes -##### - -# Solve -# θ = T/Π ( 1 - ℒ qˡ / (cᵖᵐ T)) -# for temperature T with qˡ = max(0, q - qᵛ⁺). -# root of: f(T) = T² - Π θ T - ℒ qˡ / cᵖᵐ -@inline function temperature(state::ThermodynamicState{FT}, ref, thermo) where FT - state.θ == 0 && return zero(FT) - - # Generate guess for unsaturated conditions - Π = exner_function(state, ref, thermo) - T₁ = Π * state.θ - qˡ₁ = condensate_specific_humidity(T₁, state, ref, thermo) - qˡ₁ <= 0 && return T₁ - - # If we made it this far, we have condensation - r₁ = saturation_adjustment_residual(T₁, Π, qˡ₁, state, thermo) - - ℒ = thermo.liquid.latent_heat - cᵖᵐ = mixture_heat_capacity(state.q, thermo) - T₂ = (T₁ + sqrt(T₁^2 + 4 * ℒ * qˡ₁ / cᵖᵐ)) / 2 - qˡ₂ = condensate_specific_humidity(T₂, state, ref, thermo) - r₂ = saturation_adjustment_residual(T₂, Π, qˡ₂, state, thermo) - - # Saturation adjustment - R = sqrt(max(T₂, T₁)) - ϵ = convert(FT, 1e-4) - δ = ϵ * R - iter = 0 - - while abs(r₂ - r₁) > δ - # Compute slope - ΔTΔr = (T₂ - T₁) / (r₂ - r₁) - - # Store previous values - r₁ = r₂ - T₁ = T₂ - - # Update - T₂ -= r₂ * ΔTΔr - qˡ₂ = condensate_specific_humidity(T₂, state, ref, thermo) - r₂ = saturation_adjustment_residual(T₂, Π, qˡ₂, state, thermo) - iter += 1 - end - - return T₂ -end - -@inline function saturation_adjustment_residual(T, Π, qˡ, state, thermo) - ℒᵛ = thermo.liquid.latent_heat - cᵖᵐ = mixture_heat_capacity(state.q, thermo) - return T^2 - ℒᵛ * qˡ / cᵖᵐ - Π * state.θ * T -end - -@inline function specific_volume(state, ref, thermo) - T = temperature(state, ref, thermo) - Rᵐ = mixture_gas_constant(state.q, thermo) - pᵣ = reference_pressure(state.z, ref, thermo) - return Rᵐ * T / pᵣ -end diff --git a/test/anelastic_pressure_solver_analytic.jl b/test/anelastic_pressure_solver_analytic.jl index 558a323f..05376f5b 100644 --- a/test/anelastic_pressure_solver_analytic.jl +++ b/test/anelastic_pressure_solver_analytic.jl @@ -6,8 +6,8 @@ using Oceananigans.Fields: fill_halo_regions! @testset "Anelastic pressure solver recovers analytic solution [$FT]" for FT in (Float32, Float64) grid = RectilinearGrid(FT; size=48, z=(0, 1), topology=(Flat, Flat, Bounded)) thermodynamics = ThermodynamicConstants(FT) - reference_constants = ReferenceStateConstants(FT; base_pressure=101325.0, potential_temperature=288.0) - formulation = AnelasticFormulation(grid, reference_constants, thermodynamics) + reference_state = ReferenceState(grid, thermodynamics, base_pressure=101325, potential_temperature=288) + formulation = AnelasticFormulation(reference_state) #= ρᵣ = 2 + cos(π z / 2) @@ -28,8 +28,8 @@ using Oceananigans.Fields: fill_halo_regions! ⟹ ρw = z² - z³ =# - set!(formulation.reference_density, z -> z) - fill_halo_regions!(formulation.reference_density) + set!(formulation.reference_state.density, z -> z) + fill_halo_regions!(formulation.reference_state.density) model = AtmosphereModel(grid; thermodynamics, formulation) set!(model, ρw = z -> z^2 - z^3) diff --git a/test/anelastic_pressure_solver_nonhydrostatic.jl b/test/anelastic_pressure_solver_nonhydrostatic.jl index aae55a58..2ef829e5 100644 --- a/test/anelastic_pressure_solver_nonhydrostatic.jl +++ b/test/anelastic_pressure_solver_nonhydrostatic.jl @@ -7,12 +7,12 @@ using Oceananigans z = 0:(1/Nz):1 grid = RectilinearGrid(FT; size=(Nx, Ny, Nz), x=(0, 1), y=(0, 1), z) thermodynamics = ThermodynamicConstants(FT) - reference_constants = ReferenceStateConstants(FT; base_pressure=101325, potential_temperature=288) + reference_state = ReferenceState(grid, thermodynamics) - formulation = AnelasticFormulation(grid, reference_constants, thermodynamics) - parent(formulation.reference_density) .= 1 + formulation = AnelasticFormulation(reference_state) + parent(formulation.reference_state.density) .= 1 - anelastic = AtmosphereModel(grid; thermodynamics=thermodynamics, formulation) + anelastic = AtmosphereModel(grid; thermodynamics, formulation) boussinesq = NonhydrostaticModel(; grid) uᵢ = rand(size(grid)...) diff --git a/test/atmosphere_model_unit.jl b/test/atmosphere_model_unit.jl index 3c4d25d7..3ce7a5f7 100644 --- a/test/atmosphere_model_unit.jl +++ b/test/atmosphere_model_unit.jl @@ -9,12 +9,12 @@ using Test for p₀ in (101325, 100000) for θ₀ in (288, 300) - constants = Breeze.Thermodynamics.ReferenceStateConstants(p₀, θ₀) - formulation = AnelasticFormulation(grid, constants, thermo) - model = AtmosphereModel(grid; formulation) + reference_state = ReferenceState(grid, thermo, base_pressure=p₀, potential_temperature=θ₀) + formulation = AnelasticFormulation(reference_state) + model = AtmosphereModel(grid; thermodynamics=thermo, formulation) # test set! - ρᵣ = model.formulation.reference_density + ρᵣ = model.formulation.reference_state.density cᵖᵈ = model.thermodynamics.dry_air.heat_capacity ρeᵢ = ρᵣ * cᵖᵈ * θ₀ @@ -22,7 +22,7 @@ using Test ρe₁ = deepcopy(model.energy) set!(model; ρe = ρeᵢ) - @test model.energy == ρe₁ + @test model.energy ≈ ρe₁ end end end diff --git a/test/moist_air_buoyancy.jl b/test/moist_air_buoyancy.jl index 5128753b..f6250798 100644 --- a/test/moist_air_buoyancy.jl +++ b/test/moist_air_buoyancy.jl @@ -3,18 +3,16 @@ using Oceananigans using Test @testset "NonhydrostaticModel with MoistAirBuoyancy" begin - reference_constants = ReferenceStateConstants(potential_temperature=300) - buoyancy = MoistAirBuoyancy(; reference_constants) - grid = RectilinearGrid(size=(8, 8, 8), x=(0, 400), y=(0, 400), z=(0, 400)) - model = NonhydrostaticModel(; grid, buoyancy, tracers = (:θ, :q)) + buoyancy = MoistAirBuoyancy(grid; reference_potential_temperature=300) + model = NonhydrostaticModel(; grid, buoyancy, tracers = (:θ, :qᵗ)) - θ₀ = reference_constants.reference_potential_temperature + θ₀ = buoyancy.reference_state.potential_temperature Δθ = 2 Lz = grid.Lz θᵢ(x, y, z) = θ₀ + Δθ * z / Lz - set!(model; θ = θᵢ, q = 0) + set!(model; θ = θᵢ, qᵗ = 0) # Can time-step success = try diff --git a/test/saturation_adjustment.jl b/test/saturation_adjustment.jl new file mode 100644 index 00000000..b00bf2b0 --- /dev/null +++ b/test/saturation_adjustment.jl @@ -0,0 +1,89 @@ +using Breeze +using Oceananigans +using Test + +using Breeze.Thermodynamics: + MoistureMassFractions, + PotentialTemperatureState, + exner_function, + density, + with_moisture, + saturation_specific_humidity + +using Breeze.MoistAirBuoyancies: temperature + +@testset "Saturation adjustment (MoistAirBuoyancies)" begin + for FT in (Float32, Float64) + # Minimal grid and reference state + # grid = RectilinearGrid(FT, size=(), topology=(Flat, Flat, Flat)) + grid = RectilinearGrid(FT, size=(1, 1, 1), x=(0, 1), y=(0, 1), z=(0, 1)) + thermo = ThermodynamicConstants(FT) + reference_state = ReferenceState(grid, thermo; base_pressure=101325, potential_temperature=288) + + # Sample a single cell + pᵣ = reference_state.pressure[1, 1, 1] + ρᵣ = reference_state.density[1, 1, 1] + p₀ = reference_state.base_pressure + z = FT(0.5) + + # Case 0: Absolute zero potential temperature returns zero temperature + θ₀ = zero(FT) + q₀ = MoistureMassFractions(zero(FT), zero(FT), zero(FT)) + 𝒰₀ = PotentialTemperatureState(θ₀, q₀, z, p₀, pᵣ, ρᵣ) + T₀ = temperature(𝒰₀, thermo) + @test T₀ == 0 + + # Helper for tolerances + atol_T = FT === Float64 ? 1e-6 : FT(1e-3) + + # Case 1: Unsaturated, dry (qᵗ = 0) + θ₁ = FT(300) + qᵗ₁ = zero(FT) + q₁ = MoistureMassFractions(qᵗ₁, zero(FT), zero(FT)) + 𝒰₁ = PotentialTemperatureState(θ₁, q₁, z, p₀, pᵣ, ρᵣ) + Π₁ = exner_function(𝒰₁, thermo) + T_dry₁ = Π₁ * θ₁ + + T₁ = temperature(𝒰₁, thermo) + @test isapprox(T₁, T_dry₁; atol=atol_T) + + # Case 2: Unsaturated, humid but below saturation at dry temperature + θ₂ = FT(300) + q₂ = MoistureMassFractions(zero(FT), zero(FT), zero(FT)) + 𝒰₂ = PotentialTemperatureState(θ₂, q₂, z, p₀, pᵣ, ρᵣ) + Π₂ = exner_function(𝒰₂, thermo) + T_dry₂ = Π₂ * θ₂ + + # Choose qᵗ well below saturation at T_dry₂ + ρ₂ = density(pᵣ, T_dry₂, q₂, thermo) + qᵛ⁺₂ = saturation_specific_humidity(T_dry₂, ρ₂, thermo, thermo.liquid) + qᵗ₂ = qᵛ⁺₂ / 2 + q₂ = MoistureMassFractions(qᵗ₂, zero(FT), zero(FT)) + 𝒰₂ = with_moisture(𝒰₂, q₂) + + T₂ = temperature(𝒰₂, thermo) + Π₂ = exner_function(𝒰₂, thermo) + T_dry₂ = Π₂ * θ₂ + @test isapprox(T₂, T_dry₂; atol=atol_T) + + # Case 3: Saturated, humid (qᵗ = qᵛ⁺) + T₃ = θ̃ = FT(300) + qᵗ = FT(0.025) + q̃ = MoistureMassFractions(qᵗ, zero(FT), zero(FT)) + 𝒰 = PotentialTemperatureState(θ̃, q̃, z, p₀, pᵣ, ρᵣ) + qᵛ⁺ = Breeze.MoistAirBuoyancies.adjustment_saturation_specific_humidity(T₃, 𝒰, thermo) + @test qᵗ > qᵛ⁺ # otherwise the test is wrong + + qˡ = qᵗ - qᵛ⁺ + q₃ = MoistureMassFractions(qᵛ⁺, qˡ, zero(FT)) + 𝒰₃ = with_moisture(𝒰, q₃) + Π₃ = exner_function(𝒰₃, thermo) + cᵖᵐ = mixture_heat_capacity(q₃, thermo) + ℒˡᵣ = thermo.liquid.reference_latent_heat + θ₃ = (T₃ - ℒˡᵣ / cᵖᵐ * qˡ) / Π₃ + 𝒰₃ = PotentialTemperatureState(θ₃, q₃, z, p₀, pᵣ, ρᵣ) + + T₃_solve = temperature(𝒰₃, thermo) + @test isapprox(T₃_solve, T₃; atol=atol_T) + end +end