Skip to content

Commit

Permalink
Refactor parcel
Browse files Browse the repository at this point in the history
  • Loading branch information
@trontrytel
trontrytel committed Feb 23, 2024
1 parent 7fbe344 commit 5ea2106
Show file tree
Hide file tree
Showing 16 changed files with 960 additions and 906 deletions.
92 changes: 32 additions & 60 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,18 +5,11 @@ 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)
aps = CMP.AirProperties(FT)
wps = CMP.WaterProperties(FT)
ip = CMP.IceNucleationParameters(FT)

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

# Initial conditions
Nₐ = FT(2000 * 1e3)
Expand All @@ -34,60 +27,45 @@ 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
α_m = FT(0.5) # accomodation coefficient
const_dt = FT(0.1) # model timestep
t_max = FT(100)
aerosol_1 = [CMP.DesertDust(FT), CMP.ArizonaTestDust(FT)] # aersol type for DustDeposition
aerosol_2 = [CMP.Feldspar(FT), CMP.Ferrihydrite(FT), CMP.Kaolinite(FT)] # aersol type for DepositionNucleation
ice_nucleation_modes_list =
[["MohlerRate_Deposition"], ["ActivityBasedDeposition"]]
aerosol_1 = [CMP.DesertDust(FT), CMP.ArizonaTestDust(FT)] # aersol type for DustDeposition
aerosol_2 = [CMP.Feldspar(FT), CMP.Ferrihydrite(FT), CMP.Kaolinite(FT)] # aersol type for DepositionNucleation
deposition_modes = ["MohlerRate", "ActivityBased"]
growth_modes = ["Deposition"]
droplet_size_distribution_list = [["Monodisperse"]]
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")
ax4 = MK.Axis(fig[2, 2], ylabel = "N_ice [m^-3]", xlabel = "time")

for ice_nucleation_modes in ice_nucleation_modes_list
nuc_mode = ice_nucleation_modes[1]
droplet_size_distribution = droplet_size_distribution_list[1]

if nuc_mode == "MohlerRate_Deposition"
for deposition in deposition_modes
if deposition == "MohlerRate"
for aerosol in aerosol_1
p = (;
wps,
aps,
tps,
ip,
const_dt,
r_nuc,
w,
α_m,
aerosol,
ice_nucleation_modes,
growth_modes,
droplet_size_distribution,
params = parcel_params{FT}(
const_dt = const_dt,
r_nuc = r_nuc,
w = w,
aerosol = aerosol,
deposition = deposition,
growth_modes = growth_modes,
size_distribution = size_distribution,
)

# solve ODE
sol = run_parcel(IC, FT(0), t_max, p)
sol = run_parcel(IC, FT(0), t_max, params)

# Plot results
if aerosol == CMP.DesertDust(FT)
Expand All @@ -98,33 +76,27 @@ for ice_nucleation_modes in ice_nucleation_modes_list
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
MK.lines!(ax3, sol.t, sol[6, :] * 1e3) # q_ice
MK.lines!(ax4, sol.t, sol[9, :]) # N_ice
end

elseif nuc_mode == "ActivityBasedDeposition"
elseif deposition == "ActivityBased"
for aerosol in aerosol_2
p = (;
wps,
aps,
tps,
ip,
const_dt,
r_nuc,
w,
α_m,
aerosol,
ice_nucleation_modes,
growth_modes,
droplet_size_distribution,
params = parcel_params{FT}(
const_dt = const_dt,
r_nuc = r_nuc,
w = w,
aerosol = aerosol,
deposition = deposition,
growth_modes = growth_modes,
size_distribution = size_distribution,
)

# solve ODE
sol = run_parcel(IC, FT(0), t_max, p)
sol = run_parcel(IC, FT(0), t_max, params)

# Plot results
if aerosol == CMP.Feldspar(FT)
Expand All @@ -137,7 +109,7 @@ for ice_nucleation_modes in ice_nucleation_modes_list
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
66 changes: 18 additions & 48 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,20 +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)
aps = CMP.AirProperties(FT)
wps = CMP.WaterProperties(FT)
ip = CMP.IceNucleationParameters(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 @@ -31,26 +25,23 @@ 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]

# Simulation parameters passed into ODE solver
r_nuc = FT(0.5 * 1.e-4 * 1e-6) # assumed size of nucleated particles
w = FT(0.7) # updraft speed
α_m = FT(0.5) # accomodation coefficient
const_dt = FT(1) # model timestep
t_max = FT(1200)
aerosol = CMP.Illite(FT)
ice_nucleation_modes = ["ImmersionFreezing"]
heterogeneous = "ImmersionFreezing"
growth_modes = ["Condensation", "Deposition"]
droplet_size_distribution_list = [["Monodisperse"], ["Gamma"]]
size_distribution_list = ["Monodisperse", "Gamma"]

# Plotting
fig = MK.Figure(resolution = (800, 600))
Expand All @@ -67,47 +58,26 @@ ax6 = MK.Axis(
)
MK.ylims!(ax6, 1e-6, 1)

for droplet_size_distribution in droplet_size_distribution_list
p = (;
wps,
aps,
tps,
ip,
const_dt,
r_nuc,
w,
α_m,
aerosol,
ice_nucleation_modes,
growth_modes,
droplet_size_distribution,
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 = DSD,
)
# solve ODE
sol = run_parcel(IC, FT(0), t_max, p)

DSD = droplet_size_distribution[1]
sol = run_parcel(IC, FT(0), t_max, params)

# 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
Loading

0 comments on commit 5ea2106

Please sign in to comment.