Skip to content

Commit

Permalink
final version, if only Jensen was the same
Browse files Browse the repository at this point in the history
  • Loading branch information
@trontrytel
trontrytel committed Feb 23, 2024
1 parent 764b31d commit d043a9e
Show file tree
Hide file tree
Showing 14 changed files with 733 additions and 788 deletions.
20 changes: 6 additions & 14 deletions parcel/Deposition_Nucleation.jl → parcel/Example_Deposition_Nucleation.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5,16 +5,12 @@ import CloudMicrophysics as CM
import CLIMAParameters as CP

# definition of the ODE problem for parcel model
include(joinpath(pkgdir(CM), "parcel", "parcel.jl"))
include(joinpath(pkgdir(CM), "parcel", "Parcel.jl"))
FT = Float32

# get free parameters
tps = TD.Parameters.ThermodynamicsParameters(FT)

# Constants
R_v = TD.Parameters.R_v(tps)
R_d = TD.Parameters.R_d(tps)

# Initial conditions
Nₐ = FT(2000 * 1e3)
Nₗ = FT(0)
Expand All @@ -31,17 +27,13 @@ q = TD.PhasePartition(qᵥ + qₗ + qᵢ, qₗ, qᵢ)
ts = TD.PhaseNonEquil_pTq(tps, p₀, T₀, q)
ρₐ = TD.air_density(tps, ts)
Rₐ = TD.gas_constant_air(tps, q)
Rᵥ = TD.Parameters.R_v(tps)
eₛ = TD.saturation_vapor_pressure(tps, T₀, TD.Liquid())
e = qᵥ * p₀ * R_v / Rₐ
e = eᵥ(qᵥ, p₀, Rₐ, Rᵥ)

Sₗ = FT(e / eₛ)
IC = [Sₗ, p₀, T₀, qᵥ, qₗ, qᵢ, Nₐ, Nₗ, Nᵢ, x_sulph]

ξ(T) =
TD.saturation_vapor_pressure(tps, T, TD.Liquid()) /
TD.saturation_vapor_pressure(tps, T, TD.Ice())
S_i(T, S_liq) = ξ(T) * S_liq

# Simulation parameters passed into ODE solver
r_nuc = FT(1.25e-6) # assumed size of nucleated particles
w = FT(3.5 * 1e-2) # updraft speed
Expand All @@ -54,7 +46,7 @@ growth_modes = ["Deposition"]
size_distribution = "Monodisperse"

# Plotting
fig = MK.Figure(resolution = (800, 600))
fig = MK.Figure(size = (800, 600))
ax1 = MK.Axis(fig[1, 1], ylabel = "Ice Saturation [-]")
ax2 = MK.Axis(fig[1, 2], ylabel = "Temperature [K]")
ax3 = MK.Axis(fig[2, 1], ylabel = "q_ice [g/kg]", xlabel = "time")
Expand Down Expand Up @@ -84,7 +76,7 @@ for deposition in deposition_modes
MK.lines!( # saturation
ax1,
sol.t,
S_i.(sol[3, :], sol[1, :]),
S_i.(tps, sol[3, :], sol[1, :]),
label = aero_label,
)
MK.lines!(ax2, sol.t, sol[3, :]) # temperature
Expand Down Expand Up @@ -117,7 +109,7 @@ for deposition in deposition_modes
MK.lines!( # saturation
ax1,
sol.t,
S_i.(sol[3, :], sol[1, :]),
S_i.(tps, sol[3, :], sol[1, :]),
label = aero_label,
linestyle = :dash,
)
Expand Down
38 changes: 8 additions & 30 deletions parcel/Immersion_Freezing.jl → parcel/Example_Immersion_Freezing.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5,18 +5,14 @@ import CloudMicrophysics as CM
import CLIMAParameters as CP

# definition of the ODE problem for parcel model
include(joinpath(pkgdir(CM), "parcel", "parcel.jl"))
include(joinpath(pkgdir(CM), "parcel", "Parcel.jl"))
FT = Float32
# get free parameters
tps = TD.Parameters.ThermodynamicsParameters(FT)
wps = CMP.WaterProperties(FT)

# Constants
ρₗ = wps.ρw
R_v = TD.Parameters.R_v(tps)
R_d = TD.Parameters.R_d(tps)

# Initial conditions
ρₗ = wps.ρw
Nₐ = FT(0)
Nₗ = FT(500 * 1e3)
Nᵢ = FT(0)
Expand All @@ -29,13 +25,11 @@ qᵢ = FT(0)
x_sulph = FT(0.01)

# Moisture dependent initial conditions
q = TD.PhasePartition(qᵥ + qₗ + qᵢ, qₗ, qᵢ)
ts = TD.PhaseNonEquil_pTq(tps, p₀, T₀, q)
ρₐ = TD.air_density(tps, ts)
q = TD.PhasePartition.(qᵥ + qₗ + qᵢ, qₗ, qᵢ)
R_v = TD.Parameters.R_v(tps)
Rₐ = TD.gas_constant_air(tps, q)
eₛ = TD.saturation_vapor_pressure(tps, T₀, TD.Liquid())
e = qᵥ * p₀ * R_v / Rₐ

e = eᵥ(qᵥ, p₀, Rₐ, R_v)
Sₗ = FT(e / eₛ)
IC = [Sₗ, p₀, T₀, qᵥ, qₗ, qᵢ, Nₐ, Nₗ, Nᵢ, x_sulph]

Expand Down Expand Up @@ -64,42 +58,26 @@ ax6 = MK.Axis(
)
MK.ylims!(ax6, 1e-6, 1)

for size_distribution in size_distribution_list
for DSD in size_distribution_list
params = parcel_params{FT}(
const_dt = const_dt,
r_nuc = r_nuc,
w = w,
aerosol = aerosol,
heterogeneous = heterogeneous,
growth_modes = growth_modes,
size_distribution = size_distribution,
size_distribution = DSD,
)
# solve ODE
sol = run_parcel(IC, FT(0), t_max, params)

DSD = size_distribution

# Plot results
ξ(T) =
TD.saturation_vapor_pressure(tps, T, TD.Liquid()) /
TD.saturation_vapor_pressure(tps, T, TD.Ice())
S_i(T, S_liq) = ξ(T) * S_liq - 1

MK.lines!(ax1, sol.t / 60, S_i.(sol[3, :], sol[1, :]), label = DSD)
MK.lines!(ax1, sol.t / 60, S_i.(sol[3, :], sol[1, :]) .- 1, label = DSD)
MK.lines!(ax2, sol.t / 60, sol[3, :])
MK.lines!(ax3, sol.t / 60, sol[6, :] * 1e3)
MK.lines!(ax4, sol.t / 60, sol[5, :] * 1e3)

sol_Nₗ = sol[8, :]
sol_Nᵢ = sol[9, :]
sol_qₗ = sol[5, :]
if DSD == "Monodisperse"
rₗ = cbrt.(sol_qₗ ./ sol_Nₗ ./ (4 / 3 * FT(π)) / ρₗ * ρₐ)
end
if DSD == "Gamma"
λ = cbrt.(32 .* FT(π) .* sol_Nₗ ./ sol_qₗ * ρₗ / ρₐ)
rₗ = 2 ./ λ
end
MK.lines!(ax5, sol.t / 60, sol_Nₗ)
MK.lines!(ax6, sol.t / 60, sol_Nᵢ ./ (sol_Nₗ .+ sol_Nᵢ))
end
Expand Down
30 changes: 9 additions & 21 deletions parcel/Jensen_et_al_2022.jl → parcel/Example_Jensen_et_al_2022.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -7,26 +7,13 @@ import Distributions as DS
import CloudMicrophysics.Parameters as CMP

# definition of the ODE problem for parcel model
include(joinpath(pkgdir(CM), "parcel", "parcel.jl"))

include(joinpath(pkgdir(CM), "parcel", "Parcel.jl"))
FT = Float32

# Get free parameters
tps = TD.Parameters.ThermodynamicsParameters(FT)
wps = CMP.WaterProperties(FT)

# Constants
ρₗ = wps.ρw
R_v = TD.Parameters.R_v(tps)
R_d = TD.Parameters.R_d(tps)
ϵₘ = R_d / R_v

# Saturation conversion equations
ξ(T) =
TD.saturation_vapor_pressure(tps, T, TD.Liquid()) /
TD.saturation_vapor_pressure(tps, T, TD.Ice())
S_i(T, S_liq) = @. S_liq * ξ(T) # saturation ratio

# Initial conditions
Nₐ = FT(300 * 1e6)
Nₗ = FT(0)
Expand All @@ -36,18 +23,19 @@ T₀ = FT(190)
cᵥ₀ = FT(5 * 1e-6)
x_sulph = FT(0)

# saturation and partial pressure
# Constants
ρₗ = wps.ρw
R_v = TD.Parameters.R_v(tps)
R_d = TD.Parameters.R_d(tps)
ϵₘ = R_d / R_v
eₛ = TD.saturation_vapor_pressure(tps, T₀, TD.Liquid())
qᵥ = ϵₘ / (ϵₘ - 1 + 1 / cᵥ₀)
qₗ = Nₗ * 4 / 3 * FT(π) * r₀^3 * ρₗ / FT(1.2)
qᵢ = FT(0)
q = TD.PhasePartition(qᵥ + qₗ + qᵢ, qₗ, qᵢ)
Sᵢ = FT(1.55)
Sₗ = Sᵢ / ξ(T₀)
Sₗ = Sᵢ / ξ(tps, T₀)
e = Sₗ * eₛ
p₀ = e / cᵥ₀ #TODO - ask

## Moisture dependent initial conditions
IC = [Sₗ, p₀, T₀, qᵥ, qₗ, qᵢ, Nₐ, Nₗ, Nᵢ, x_sulph]

# Simulation parameters passed into ODE solver
Expand All @@ -71,7 +59,7 @@ Jensen_t_ICNC = [0.217, 42.69, 50.02, 54.41, 58.97, 65.316, 72.477, 82.08, 92.65
Jensen_ICNC = [0, 0, 0.282, 0.789, 1.804, 4.1165, 7.218, 12.12, 16.35, 16.8, 16.97, 17.086]
#! format: on

fig = MK.Figure(resolution = (1000, 1000))
fig = MK.Figure(size = (1000, 1000))
ax1 = MK.Axis(fig[1, 1], ylabel = "Saturation")
ax2 = MK.Axis(fig[3, 1], xlabel = "Time [s]", ylabel = "Temperature [K]")
ax3 = MK.Axis(fig[2, 1], ylabel = "q_vap [g/kg]")
Expand Down Expand Up @@ -102,7 +90,7 @@ params = parcel_params{FT}(
sol = run_parcel(IC, FT(0), t_max, params)

# Plot results
MK.lines!(ax1, sol.t, S_i(sol[3, :], (sol[1, :])), label = "ice", color = :blue)
MK.lines!(ax1, sol.t, S_i.(tps, sol[3, :], (sol[1, :])), label = "ice", color = :blue)
MK.lines!(ax1, sol.t, (sol[1, :]), label = "liquid", color = :red) # liq saturation
MK.lines!(ax2, sol.t, sol[3, :], label = "CM.jl") # temperature
MK.lines!(ax3, sol.t, sol[4, :] * 1e3) # q_vap
Expand Down
53 changes: 25 additions & 28 deletions parcel/Liquid_only.jl → parcel/Example_Liquid_only.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,19 +2,19 @@ import OrdinaryDiffEq as ODE
import CairoMakie as MK
import Thermodynamics as TD
import CloudMicrophysics as CM
import CloudMicrophysics.Parameters as CMP
import CLIMAParameters as CP

# definition of the ODE problem for parcel model
include(joinpath(pkgdir(CM), "parcel", "parcel.jl"))
include(joinpath(pkgdir(CM), "parcel", "Parcel.jl"))

FT = Float32

# Get free parameters
tps = TD.Parameters.ThermodynamicsParameters(FT)
wps = CMP.WaterProperties(FT)

# Constants
ρₗ = wps.ρw
ρᵢ = wps.ρi
R_v = TD.Parameters.R_v(tps)
R_d = TD.Parameters.R_d(tps)

Expand All @@ -26,25 +26,14 @@ r₀ = FT(8e-6)
p₀ = FT(800 * 1e2)
T₀ = FT(273.15 + 7.0)
x_sulph = FT(0)

# mass per volume for dry air, vapor and liquid
eₛ = TD.saturation_vapor_pressure(tps, T₀, TD.Liquid())
e = eₛ
e = TD.saturation_vapor_pressure(tps, T₀, TD.Liquid())
Sₗ = FT(1)
md_v = (p₀ - e) / R_d / T₀
mv_v = e / R_v / T₀
ml_v = Nₗ * 4 / 3 * FT(π) * ρₗ * r₀^3

qᵥ = mv_v / (md_v + mv_v + ml_v)
qₗ = ml_v / (md_v + mv_v + ml_v)
qᵢ = FT(0)

# Moisture dependent initial conditions
q = TD.PhasePartition(qᵥ + qₗ + qᵢ, qₗ, qᵢ)
ts = TD.PhaseNonEquil_pTq(tps, p₀, T₀, q)
ρₐ = TD.air_density(tps, ts)

Rₐ = TD.gas_constant_air(tps, q)
Sₗ = FT(e / eₛ)
IC = [Sₗ, p₀, T₀, qᵥ, qₗ, qᵢ, Nₐ, Nₗ, Nᵢ, x_sulph]

# Simulation parameters passed into ODE solver
Expand All @@ -62,44 +51,52 @@ Rogers_time_radius = [0.561, 2, 3.99, 10.7, 14.9, 19.9]
Rogers_radius = [8.0, 8.08, 8.26, 8.91, 9.26, 9.68]
#! format: on

fig = MK.Figure(resolution = (800, 600))
# Setup the plots
fig = MK.Figure(size = (800, 600))
ax1 = MK.Axis(fig[1, 1], ylabel = "Supersaturation [%]")
ax2 = MK.Axis(fig[3, 1], xlabel = "Time [s]", ylabel = "Temperature [K]")
ax3 = MK.Axis(fig[2, 1], ylabel = "q_vap [g/kg]")
ax4 = MK.Axis(fig[2, 2], xlabel = "Time [s]", ylabel = "q_liq [g/kg]")
ax5 = MK.Axis(fig[1, 2], ylabel = "radius [μm]")

MK.lines!(ax1, Rogers_time_supersat, Rogers_supersat, label = "Rogers_1975")
MK.lines!(ax5, Rogers_time_radius, Rogers_radius)

for size_distribution in size_distribution_list
for DSD in size_distribution_list
params = parcel_params{FT}(
size_distribution = size_distribution,
size_distribution = DSD,
const_dt = const_dt,
w = w,
)
# solve ODE
sol = run_parcel(IC, FT(0), t_max, params)

DSD = size_distribution

# Plot results
MK.lines!(ax1, sol.t, (sol[1, :] .- 1) * 100.0, label = DSD)
MK.lines!(ax2, sol.t, sol[3, :])
MK.lines!(ax3, sol.t, sol[4, :] * 1e3)
MK.lines!(ax4, sol.t, sol[5, :] * 1e3)

sol_Nₗ = sol[8, :]
sol_Nᵢ = sol[9, :]
# Compute the current air density
sol_T = sol[3, :]
sol_p = sol[2, :]
sol_qᵥ = sol[4, :]
sol_qₗ = sol[5, :]
rₗ = DSD == "Monodisperse" ?
cbrt.(sol_qₗ ./ sol_Nₗ ./ (4 / 3 * FT(π)) / ρₗ * ρₐ) :
DSD == "Gamma" ?
2 ./ (cbrt.(32 .* FT(π) .* sol_Nₗ ./ sol_qₗ * ρₗ / ρₐ)) :
FT(0)
sol_qᵢ = sol[6, :]
q = TD.PhasePartition.(sol_qᵥ + sol_qₗ + sol_qᵢ, sol_qₗ, sol_qᵢ)
ts = TD.PhaseNonEquil_pTq.(tps, sol_p, sol_T, q)
ρₐ = TD.air_density.(tps, ts)
# Compute the mean particle size based on the distribution
distr = sol.prob.p.distr
moms = distribution_moments.(distr, sol_qₗ, sol_Nₗ, ρₗ, ρₐ, sol_qᵢ, Nᵢ, ρᵢ)
rₗ = similar(sol_T)
for it in range(1, length(sol_T))
rₗ[it] = moms[it].rₗ
end
MK.lines!(ax5, sol.t, rₗ * 1e6)
end

#fig[3, 2] = MK.Legend(fig, ax1, framevisible = false)
MK.axislegend(
ax1,
framevisible = false,
Expand Down
Loading

0 comments on commit d043a9e

Please sign in to comment.