Merge pull request #1 from navidcy/saturation-adjustment

Implement saturation adjustment for Boussinesq LES with condensation

Author: glwagner

Project.toml

@@ -4,7 +4,11 @@ authors = ["Navid C. Constantinou <[email protected]>, Madi, Greg
 version = "0.1.0"
 
 [deps]
+JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819"
 Oceananigans = "9e8cae18-63c1-5223-a75c-80ca9d6e9a09"
+RootSolvers = "7181ea78-2dcb-4de3-ab41-2b8ab5a31e74"
 
 [compat]
+JLD2 = "0.5.13"
 Oceananigans = "0.96"
+RootSolvers = "0.4.4"

examples/JRA55_saturation_specific_humidity.jl

@@ -0,0 +1,39 @@
+using JLD2
+using AquaSkyLES
+using AquaSkyLES: saturation_specific_humidity, AtmosphereThermodynamics
+using GLMakie
+
+@load "JRA55_atmospheric_state_Jan_1_1991.jld2" q T p
+
+thermo = AtmosphereThermodynamics()
+
+ρ = 1.2
+qᵛ★ = saturation_specific_humidity.(T, 1.2, Ref(thermo))
+qˡ = @. max(0, q - q★)
+
+fig = Figure(size=(1200, 600))
+ax1 = Axis(fig[1, 2], title="Temperature")
+ax2 = Axis(fig[2, 1], title="Saturation specific humidity")
+ax3 = Axis(fig[1, 1], title="Total specific humidity")
+ax4 = Axis(fig[2, 2], title="Liquid specific humidity")
+# heatmap!(ax1, T, colormap=:magma)
+# heatmap!(ax2, qᵛ★)
+# heatmap!(ax3, q, colormap=:grays)
+# heatmap!(ax4, qˡ, colormap=:grays, colorrange=(0, 0.001))
+
+# Compute cloudiness for instantaneous drop
+θ⁻ = T .- 10
+Ψ = AquaSkyLES.ThermodynamicState{Float64}.(θ⁻, q, 0)
+ℛ = AquaSkyLES.ReferenceState{Float64}(101325, 20)
+T⁻ = AquaSkyLES.temperature.(Ψ, Ref(ℛ), Ref(thermo))
+qᵛ★ = saturation_specific_humidity.(T⁻, 1.2, Ref(thermo))
+qˡ = @. max(0, q - qᵛ★)
+Π = AquaSkyLES.exner_function.(Ψ, Ref(ℛ), Ref(thermo))
+Tu = @. Π * θ⁻
+ΔT = T⁻ - Tu
+
+heatmap!(ax1, ΔT, colormap=:magma)
+heatmap!(ax2, qᵛ★)
+heatmap!(ax3, q, colormap=:grays)
+heatmap!(ax4, qˡ, colormap=:grays, colorrange=(0, 0.01))
+display(fig)

examples/free_convection.jl

@@ -1,46 +1,95 @@
 using Oceananigans
 using Oceananigans.Units
+using Printf
 using AquaSkyLES
 
-Nx = Nz = 256
-Lz = 3 * 1024
+Nx = Nz = 128
+Lz = 4 * 1024
 grid = RectilinearGrid(size=(Nx, Nz), x=(0, 2Lz), z=(0, Lz), topology=(Periodic, Flat, Bounded))
 
-ρ₀ = 1   # air density
-cₚ = 1e3 # air specific heat
-Q = 1000 # heat flux in W / m²
-Jθ = Q / (ρ₀ * cₚ)
-heating = FluxBoundaryCondition(Jθ)
-θ_bcs = FieldBoundaryConditions(bottom=heating)
+p₀ = 101325 # Pa
+θ₀ = 288 # K
+reference_state = AquaSkyLES.ReferenceState(base_pressure=p₀, potential_temperature=θ₀)
+buoyancy = AquaSkyLES.MoistAirBuoyancy(; reference_state)
+
+ρ₀ = AquaSkyLES.base_density(buoyancy) # 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-3)
+q_bcs = FieldBoundaryConditions(bottom=vapor_flux)
 
 advection = WENO() #(momentum=WENO(), θ=WENO(), q=WENO(bounds=(0, 1)))
 tracers = (:θ, :q)
-buoyancy = AquaSkyLES.MoistAirBuoyancy()
-model = NonhydrostaticModel(; grid, advection, tracers, buoyancy, boundary_conditions=(; θ=θ_bcs))
+model = NonhydrostaticModel(; grid, advection, buoyancy,
+                            tracers = (:θ, :q),
+                            boundary_conditions = (θ=θ_bcs, q=q_bcs))
 
 Lz = grid.Lz
-Δθ = 10 # ᵒC
-Tₛ = 20 # ᵒC
+Δθ = 5 # K
+Tₛ = reference_state.θ # K
 θᵢ(x, z) = Tₛ + Δθ * z / Lz + 1e-2 * Δθ * randn()
-qᵢ(x, z) = 1e-4 * rand()
+qᵢ(x, z) = 1e-2 + 1e-5 * rand()
 set!(model, θ=θᵢ, q=qᵢ)
 
-simulation = Simulation(model, Δt=10, stop_time=1hour)
+simulation = Simulation(model, Δt=10, stop_time=2hours)
 conjure_time_step_wizard!(simulation, cfl=0.7)
-progress(sim) = @info string(iteration(sim), ": ", prettytime(sim))
+
+T = AquaSkyLES.TemperatureField(model)
+qˡ = AquaSkyLES.CondensateField(model, T)
+qᵛ★ = AquaSkyLES.SaturationField(model, T)
+δ = Field(model.tracers.q - qᵛ★)
+
+function progress(sim)
+    compute!(T)
+    compute!(qˡ)
+    compute!(δ)
+    q = sim.model.tracers.q
+    θ = sim.model.tracers.θ
+    u, v, w = sim.model.velocities
+
+    umax = maximum(u)
+    vmax = maximum(v)
+    wmax = maximum(w)
+
+    qmin = minimum(q)
+    qmax = maximum(q)
+    qˡmax = maximum(qˡ)
+    δmax = maximum(δ)
+
+    θmin = minimum(θ)
+    θmax = maximum(θ)
+
+    msg = @sprintf("Iter: %d, t = %s, max|u|: (%.2e, %.2e, %.2e)",
+                    iteration(sim), prettytime(sim), umax, vmax, wmax)
+
+    msg *= @sprintf(", extrema(q): (%.2e, %.2e), max(qˡ): %.2e, min(δ): %.2e, extrema(θ): (%.2e, %.2e)",
+                     qmin, qmax, qˡmax, δmax, θmin, θmax)
+
+    @info msg
+
+    return nothing
+end
+
 add_callback!(simulation, progress, IterationInterval(10))
 
-ow = JLD2Writer(model, merge(model.velocities, model.tracers),
+outputs = merge(model.velocities, model.tracers, (; T, qˡ, qᵛ★))
+
+ow = JLD2Writer(model, outputs,
                 filename = "free_convection.jld2",
-                schedule = TimeInterval(1minutes),
+                schedule = IterationInterval(10),
                 overwrite_existing = true)
 
 simulation.output_writers[:jld2] = ow
 
 run!(simulation)
 
 θt = FieldTimeSeries("free_convection.jld2", "θ")
+Tt = FieldTimeSeries("free_convection.jld2", "T")
 qt = FieldTimeSeries("free_convection.jld2", "q")
+qˡt = FieldTimeSeries("free_convection.jld2", "qˡ")
 times = qt.times
 Nt = length(θt)
 
@@ -50,22 +99,35 @@ n = Observable(length(θt))
 
 θn = @lift θt[$n]
 qn = @lift qt[$n]
+Tn = @lift Tt[$n]
+qˡn = @lift qˡt[$n]
 title = @lift "t = $(prettytime(times[$n]))"
 
-fig = Figure(size=(1600, 400), fontsize=22)
+fig = Figure(size=(800, 400), fontsize=12)
 axθ = Axis(fig[1, 1], xlabel="x (m)", ylabel="z (m)")
 axq = Axis(fig[1, 2], xlabel="x (m)", ylabel="z (m)")
+axT = Axis(fig[2, 1], xlabel="x (m)", ylabel="z (m)")
+axqˡ = Axis(fig[2, 2], xlabel="x (m)", ylabel="z (m)")
 
 fig[0, :] = Label(fig, title, fontsize=22, tellwidth=false)
 
-hmθ = heatmap!(axθ, θn, colorrange=(Tₛ, Tₛ+Δθ))
-hmq = heatmap!(axq, qn, colorrange=(0, 8e-5), colormap=:magma)
-
-Label(fig[0, 1], "θ", tellwidth=false)
-Label(fig[0, 2], "q", tellwidth=false)
+Tmin = minimum(Tt)
+Tmax = maximum(Tt)
 
-Colorbar(fig[2, 1], hmθ, label = "[ᵒC]", vertical=false)
-Colorbar(fig[2, 2], hmq, label = "", vertical=false)
+hmθ = heatmap!(axθ, θn, colorrange=(Tₛ, Tₛ+Δθ))
+hmq = heatmap!(axq, qn, colorrange=(0, 2e-2), colormap=:magma)
+hmT = heatmap!(axT, Tn, colorrange=(Tmin, Tmax))
+hmqˡ = heatmap!(axqˡ, qˡn, colorrange=(0, 2e-4), colormap=:magma)
+
+# Label(fig[0, 1], "θ", tellwidth=false)
+# Label(fig[0, 2], "q", tellwidth=false)
+# Label(fig[0, 1], "θ", tellwidth=false)
+# Label(fig[0, 2], "q", tellwidth=false)
+
+Colorbar(fig[1, 0], hmθ, label = "θ [K]", vertical=true)
+Colorbar(fig[1, 3], hmq, label = "q", vertical=true)
+Colorbar(fig[2, 0], hmT, label = "T [K]", vertical=true)
+Colorbar(fig[2, 3], hmqˡ, label = "qˡ", vertical=true)
 
 fig
 

examples/large_yeager_saturation_specific_humidity.jl

@@ -0,0 +1,19 @@
+using GLMakie
+using AquaSkyLES
+
+thermo = AquaSkyLES.AtmosphereThermodynamics()
+
+saturation_specific_humidity_large_yeager(T, ρ) = 640380 * exp(-5107.4 / T) / ρ
+
+ρ = 1.2
+T₀ = 273.15
+T = collect(T₀:0.01:T₀+50)
+q★_large_yeager = saturation_specific_humidity_large_yeager.(T, ρ)
+q★_aqua_sky = AquaSkyLES.saturation_specific_humidity.(T, ρ, Ref(thermo))
+
+fig = Figure()
+ax = Axis(fig[1, 1], xlabel = "Temperature (K)", ylabel = "Saturation Specific Humidity (g/kg)")
+lines!(ax, T, q★_large_yeager, label = "Large and Yeager (2009)")
+lines!(ax, T, q★_aqua_sky, label = "AquaSkyLES")
+Legend(fig[1, 2], ax)
+display(fig)