Skip to content

Commit

Permalink
same opts as climaatmos
Browse files Browse the repository at this point in the history
  • Loading branch information
haakon-e committed Feb 28, 2025
1 parent 40659dc commit 1229da4
Show file tree
Hide file tree
Showing 79 changed files with 1,741 additions and 504 deletions.
4 changes: 3 additions & 1 deletion .JuliaFormatter.toml
Original file line number Diff line number Diff line change
@@ -1,3 +1,5 @@
margin = 80
indent = 4
always_for_in = true
whitespace_typedefs = true
whitespace_ops_in_indices = true
whitespace_ops_in_indices = true
7 changes: 6 additions & 1 deletion .dev/up_deps.jl
Original file line number Diff line number Diff line change
Expand Up @@ -4,7 +4,12 @@ files in all of our environments.
=#

root = dirname(@__DIR__)
dirs = (root, joinpath(root, "test"), joinpath(root, "docs"), joinpath(root, "parcel"))
dirs = (
root,
joinpath(root, "test"),
joinpath(root, "docs"),
joinpath(root, "parcel"),
)

cd(root) do
for dir in dirs
Expand Down
8 changes: 6 additions & 2 deletions box/Alpert_Knopf_2016_forward.jl
Original file line number Diff line number Diff line change
Expand Up @@ -49,8 +49,12 @@ p_def = (; tps, A_aero, aerosol, cooling_rate, N₀)

# run the box model without and with distribution sampling
sol_cst = run_box(IC, t_0, t_end, (p_def..., const_dt = dt, flag = false))
sol_var =
run_box(IC, t_0, t_end, (p_def..., const_dt = dt, flag = true, Aj_sorted = Aj_sorted))
sol_var = run_box(
IC,
t_0,
t_end,
(p_def..., const_dt = dt, flag = true, Aj_sorted = Aj_sorted),
)

# Compute the immersion freezing rate for the no-sampling case
Δa = @. FT(1) - CMO.a_w_ice(tps, sol_cst[1, :])
Expand Down
7 changes: 6 additions & 1 deletion docs/make.jl
Original file line number Diff line number Diff line change
Expand Up @@ -72,7 +72,12 @@ pages = Any[
]

mathengine = MathJax(
Dict(:TeX => Dict(:equationNumbers => Dict(:autoNumber => "AMS"), :Macros => Dict())),
Dict(
:TeX => Dict(
:equationNumbers => Dict(:autoNumber => "AMS"),
:Macros => Dict(),
),
),
)

format = Documenter.HTML(
Expand Down
3 changes: 2 additions & 1 deletion docs/src/guides/GettingStarted.jl
Original file line number Diff line number Diff line change
Expand Up @@ -35,5 +35,6 @@ nothing #hide

# Finally, we call `accretion`, which will return the accretion rates for
# cloud and rain water specific humidities, as well as cloud and rain water number concentrations.
(; dq_rai_dt, dq_liq_dt, dN_rai_dt, dN_liq_dt) = CM2.accretion(SB2006, qₗ, qᵣ, ρₐ, Nₗ)
(; dq_rai_dt, dq_liq_dt, dN_rai_dt, dN_liq_dt) =
CM2.accretion(SB2006, qₗ, qᵣ, ρₐ, Nₗ)
@info("Accretion rates: ", dq_rai_dt, dq_liq_dt, dN_rai_dt, dN_liq_dt)
6 changes: 4 additions & 2 deletions docs/src/guides/Parameters.jl
Original file line number Diff line number Diff line change
Expand Up @@ -46,8 +46,10 @@ nothing #hide
# Overwriting the parameters can also be done using dictionaries instead of toml files.
# As an example we create a dictionary were we define a different value of the
# rain autoconversion time scale and create another parameter struct based on it.
override_file =
Dict("rain_autoconversion_timescale" => Dict("value" => 13, "type" => "float"))
override_file = Dict(
"rain_autoconversion_timescale" =>
Dict("value" => 13, "type" => "float"),
)
toml_dict2 = CP.create_toml_dict(FT; override_file)
const overwrite2 = CMP.Rain(toml_dict2)

Expand Down
34 changes: 23 additions & 11 deletions docs/src/plots/ARGplots.jl
Original file line number Diff line number Diff line change
Expand Up @@ -48,7 +48,9 @@ function mass2vol(mass_mixing_ratios)
else
densities = (sulfate.ρ,)
end
volfractions = (mass_mixing_ratios ./ densities) ./ sum(mass_mixing_ratios ./ densities)
volfractions =
(mass_mixing_ratios ./ densities) ./
sum(mass_mixing_ratios ./ densities)
return volfractions
end

Expand Down Expand Up @@ -160,9 +162,11 @@ function make_ARG_figX(X)
end
AD = AM.AerosolDistribution((paper_mode_1, paper_mode_2))
act_frac1[it] =
AA.N_activated_per_mode(ap, AD, aip, tps, T, p, w, q)[1] / N_1
AA.N_activated_per_mode(ap, AD, aip, tps, T, p, w, q)[1] /
N_1
act_frac2[it] =
AA.N_activated_per_mode(ap, AD, aip, tps, T, p, w, q)[2] / N2i
AA.N_activated_per_mode(ap, AD, aip, tps, T, p, w, q)[2] /
N2i
global it += 1
end

Expand Down Expand Up @@ -211,9 +215,11 @@ function make_ARG_figX(X)
end
AD = AM.AerosolDistribution((paper_mode_1, paper_mode_2))
act_frac1[it] =
AA.N_activated_per_mode(ap, AD, aip, tps, T, p, w, q)[1] / N_1
AA.N_activated_per_mode(ap, AD, aip, tps, T, p, w, q)[1] /
N_1
act_frac2[it] =
AA.N_activated_per_mode(ap, AD, aip, tps, T, p, w, q)[2] / N2i
AA.N_activated_per_mode(ap, AD, aip, tps, T, p, w, q)[2] /
N2i
global it += 1
end

Expand Down Expand Up @@ -264,9 +270,11 @@ function make_ARG_figX(X)
end
AD = AM.AerosolDistribution((paper_mode_1, paper_mode_2))
act_frac1[it] =
AA.N_activated_per_mode(ap, AD, aip, tps, T, p, w, q)[1] / N_1
AA.N_activated_per_mode(ap, AD, aip, tps, T, p, w, q)[1] /
N_1
act_frac2[it] =
AA.N_activated_per_mode(ap, AD, aip, tps, T, p, w, q)[2] / N_2
AA.N_activated_per_mode(ap, AD, aip, tps, T, p, w, q)[2] /
N_2
global it += 1
end

Expand Down Expand Up @@ -314,9 +322,11 @@ function make_ARG_figX(X)
end
AD = AM.AerosolDistribution((paper_mode_1, paper_mode_2))
act_frac1[it] =
AA.N_activated_per_mode(ap, AD, aip, tps, T, p, w, q)[1] / N_1
AA.N_activated_per_mode(ap, AD, aip, tps, T, p, w, q)[1] /
N_1
act_frac2[it] =
AA.N_activated_per_mode(ap, AD, aip, tps, T, p, w, q)[2] / N_2
AA.N_activated_per_mode(ap, AD, aip, tps, T, p, w, q)[2] /
N_2
global it += 1
end

Expand Down Expand Up @@ -366,9 +376,11 @@ function make_ARG_figX(X)

for wi in w
act_frac1[it] =
AA.N_activated_per_mode(ap, AD, aip, tps, T, p, wi, q)[1] / N_1
AA.N_activated_per_mode(ap, AD, aip, tps, T, p, wi, q)[1] /
N_1
act_frac2[it] =
AA.N_activated_per_mode(ap, AD, aip, tps, T, p, wi, q)[2] / N_2
AA.N_activated_per_mode(ap, AD, aip, tps, T, p, wi, q)[2] /
N_2
global it += 1
end

Expand Down
17 changes: 14 additions & 3 deletions docs/src/plots/ARGplots_fig1.jl
Original file line number Diff line number Diff line change
Expand Up @@ -68,7 +68,8 @@ for N_2 in N_2_range
(sulfate.ρ,),
)
AD_B = AM.AerosolDistribution((paper_mode_1_B, paper_mode_2_B))
N_act_frac_B[it] = AA.N_activated_per_mode(ap, AD_B, aip, tps, T, p, w, q)[1] / N_1
N_act_frac_B[it] =
AA.N_activated_per_mode(ap, AD_B, aip, tps, T, p, w, q)[1] / N_1
global it += 1
end

Expand All @@ -83,7 +84,17 @@ PL.plot(
xlabel = "Mode 2 aerosol number concentration [1/cm3]",
ylabel = "Mode 1 number fraction activated",
)
PL.scatter!(Fig1_x_obs, Fig1_y_obs, markercolor = :black, label = "paper observations")
PL.plot!(Fig1_x_param, Fig1_y_param, linecolor = :black, label = "paper parameterization")
PL.scatter!(
Fig1_x_obs,
Fig1_y_obs,
markercolor = :black,
label = "paper observations",
)
PL.plot!(
Fig1_x_param,
Fig1_y_param,
linecolor = :black,
label = "paper parameterization",
)

PL.savefig("Abdul-Razzak_and_Ghan_fig_1.svg")
15 changes: 10 additions & 5 deletions docs/src/plots/Baumgartner2022_fig5.jl
Original file line number Diff line number Diff line change
Expand Up @@ -35,22 +35,27 @@ sol = RS.find_zero(fun, RS.SecantMethod(FT(0), FT(0.5)), RS.CompactSolution())
a_w = Δa_crit .+ a_w_ice

# various ways Baumgartner calculates vapor pressures
Baumgartner_p_ice(T) = exp(9.550426 - 5723.265 / T + 3.53068 * log(T) - 0.00728332 * T)
Baumgartner_p_ice(T) =
exp(9.550426 - 5723.265 / T + 3.53068 * log(T) - 0.00728332 * T)
MK_p_liq(T) = exp(
54.842763 - 6763.22 / T - 4.21 * log(T) +
0.000367 * T +
tanh(0.0415 * (T - 218.8)) * (53.878 - 1331.22 / T - 9.44523 * log(T) + 0.014025 * T),
tanh(0.0415 * (T - 218.8)) *
(53.878 - 1331.22 / T - 9.44523 * log(T) + 0.014025 * T),
)
Tab_p_liq(T) = exp(
100 * (18.4524 - 3505.15788 / T - 330918.55 / (T^2) + 12725068.26 / (T^3)),
)
Tab_p_liq(T) =
exp(100 * (18.4524 - 3505.15788 / T - 330918.55 / (T^2) + 12725068.26 / (T^3)))
Nach_p_liq(T) = exp(74.8727 - 7167.405 / T - 7.77107 * log(T) + 0.00505 * T)

a_w_MK = [Δa_crit + (Baumgartner_p_ice(T)) / (MK_p_liq(T)) for T in T_range]
a_w_Tab = [Δa_crit + (Baumgartner_p_ice(T)) / (Tab_p_liq(T)) for T in T_range]
a_w_Nach = [Δa_crit + (Baumgartner_p_ice(T)) / (Nach_p_liq(T)) for T in T_range]
a_w_Luo = [
Δa_crit +
(Baumgartner_p_ice(T)) / (CMO.H2SO4_soln_saturation_vapor_pressure(H2SO4_prs, 0.0, T)) for T in T_range
(Baumgartner_p_ice(T)) /
(CMO.H2SO4_soln_saturation_vapor_pressure(H2SO4_prs, 0.0, T)) for
T in T_range
]

#! format: off
Expand Down
36 changes: 30 additions & 6 deletions docs/src/plots/Cholette2019_fig1.jl
Original file line number Diff line number Diff line change
Expand Up @@ -67,14 +67,20 @@ function fig1()
F_liqs = (FT(0), FT(0.33), FT(0.67), FT(1))
F_rs = (FT(0), FT(1 - eps(FT)))

labels = ["F_liq = 0", "F_liq = 0.33", "F_liq = 0.67", "F_liq = 1", "ϕᵢ = 1"]
labels =
["F_liq = 0", "F_liq = 0.33", "F_liq = 0.67", "F_liq = 1", "ϕᵢ = 1"]
colors = [:black, :blue, :red, :green]
res = 50

fig = Plt.Figure(size = (1200, 400))

#F_r = 0
ax1 = Plt.Axis(fig[1:7, 1:9], title = "Fᵣ = 0", xlabel = "Dₚ (m)", ylabel = "V (m s⁻¹)")
ax1 = Plt.Axis(
fig[1:7, 1:9],
title = "Fᵣ = 0",
xlabel = "Dₚ (m)",
ylabel = "V (m s⁻¹)",
)

for i in 1:4
D_ps, V_0, V_0_ϕ = get_values(
Expand All @@ -87,12 +93,23 @@ function fig1()
res,
)
Plt.lines!(ax1, D_ps, V_0_ϕ, color = colors[i], label = labels[i])
Plt.lines!(ax1, D_ps, V_0, color = colors[i], label = labels[i], linestyle = :dash)
Plt.lines!(
ax1,
D_ps,
V_0,
color = colors[i],
label = labels[i],
linestyle = :dash,
)
end

# F_r = 1
ax2 =
Plt.Axis(fig[1:7, 10:18], title = "Fᵣ = 1", xlabel = "Dₚ (m)", ylabel = "V (m s⁻¹)")
ax2 = Plt.Axis(
fig[1:7, 10:18],
title = "Fᵣ = 1",
xlabel = "Dₚ (m)",
ylabel = "V (m s⁻¹)",
)
for i in 1:4
D_ps, V_1, V_1_ϕ = get_values(
p3,
Expand All @@ -104,7 +121,14 @@ function fig1()
res,
)
Plt.lines!(ax2, D_ps, V_1_ϕ, color = colors[i], label = labels[i])
Plt.lines!(ax2, D_ps, V_1, color = colors[i], label = labels[i], linestyle = :dash)
Plt.lines!(
ax2,
D_ps,
V_1,
color = colors[i],
label = labels[i],
linestyle = :dash,
)
end

Plt.Legend(
Expand Down
9 changes: 7 additions & 2 deletions docs/src/plots/CloudDiagnostics.jl
Original file line number Diff line number Diff line change
Expand Up @@ -15,8 +15,13 @@ rain = CMP.Rain(FT)
cloud_liquid = CMP.CloudLiquid(FT)
cloud_ice = CMP.CloudIce(FT)

override_file =
joinpath(pkgdir(CloudMicrophysics), "src", "parameters", "toml", "SB2006_limiters.toml")
override_file = joinpath(
pkgdir(CloudMicrophysics),
"src",
"parameters",
"toml",
"SB2006_limiters.toml",
)
toml_dict = CP.create_toml_dict(FT; override_file)
SB = CMP.SB2006(toml_dict)
SB_no_limiters = CMP.SB2006(toml_dict, false)
Expand Down
12 changes: 9 additions & 3 deletions docs/src/plots/Frostenberg_fig1.jl
Original file line number Diff line number Diff line change
Expand Up @@ -13,8 +13,8 @@ INPC_range = 10.0 .^ (range(-5, stop = 7, length = 500)) #ice nucleating particl

frequency = [
IN.INP_concentration_frequency(ip, INPC, T) > 0.015 ?
IN.INP_concentration_frequency(ip, INPC, T) : missing for INPC in INPC_range,
T in T_range
IN.INP_concentration_frequency(ip, INPC, T) : missing for
INPC in INPC_range, T in T_range
]
mu = [exp(IN.INP_concentration_mean(T)) for T in T_range]

Expand All @@ -33,7 +33,13 @@ PL.contourf(
lw = 0,
)

PL.plot!(T_range, mu, label = "median INPC", legend = :topright, color = :darkred)
PL.plot!(
T_range,
mu,
label = "median INPC",
legend = :topright,
color = :darkred,
)

PL.plot!(
repeat([257], 50),
Expand Down
15 changes: 12 additions & 3 deletions docs/src/plots/HomFreezingPlots.jl
Original file line number Diff line number Diff line change
Expand Up @@ -18,9 +18,18 @@ T_range = range(229.0, stop = 234.5, length = 50) # air temperature
x_sulph = Vector{FT}([0.03, 0.04, 0.06]) # wt% sulphuric acid in droplets

# Solving for Δa and J values
Δa1 = [CMO.a_w_xT(H2SO4_prs, tps, x_sulph[1], T) - CMO.a_w_ice(tps, T) for T in T_range]
Δa2 = [CMO.a_w_xT(H2SO4_prs, tps, x_sulph[2], T) - CMO.a_w_ice(tps, T) for T in T_range]
Δa3 = [CMO.a_w_xT(H2SO4_prs, tps, x_sulph[3], T) - CMO.a_w_ice(tps, T) for T in T_range]
Δa1 = [
CMO.a_w_xT(H2SO4_prs, tps, x_sulph[1], T) - CMO.a_w_ice(tps, T) for
T in T_range
]
Δa2 = [
CMO.a_w_xT(H2SO4_prs, tps, x_sulph[2], T) - CMO.a_w_ice(tps, T) for
T in T_range
]
Δa3 = [
CMO.a_w_xT(H2SO4_prs, tps, x_sulph[3], T) - CMO.a_w_ice(tps, T) for
T in T_range
]

J1 = @. CMI.homogeneous_J_cubic(ip.homogeneous, Δa1)
J2 = @. CMI.homogeneous_J_cubic(ip.homogeneous, Δa2)
Expand Down
19 changes: 16 additions & 3 deletions docs/src/plots/KnopfAlpert2013_fig5.jl
Original file line number Diff line number Diff line change
Expand Up @@ -14,8 +14,16 @@ illite = CMP.Illite(FT) # dust type

# Knopf and Alpert 2013 Figure 5A: Cirrus
# https://doi.org/10.1039/C3FD00035D
temp =
[228.20357, 228.33571, 228.50357, 228.75000, 228.92143, 229.16786, 229.39643, 229.52143]
temp = [
228.20357,
228.33571,
228.50357,
228.75000,
228.92143,
229.16786,
229.39643,
229.52143,
]
KA13_fig5A_J = [
70170382.86704,
12101528.74384,
Expand Down Expand Up @@ -43,7 +51,12 @@ J_ABIFM = @. CMI.ABIFM_J(illite, Δa_w) * 1e-4 # converted from SI units to cm^-

# Plot results
fig = MK.Figure(resolution = (800, 600))
ax1 = MK.Axis(fig[1, 1], ylabel = "J_het [cm^-2 s^-1]", xlabel = "Temp [K]", yscale = log10)
ax1 = MK.Axis(
fig[1, 1],
ylabel = "J_het [cm^-2 s^-1]",
xlabel = "Temp [K]",
yscale = log10,
)

MK.ylims!(ax1, FT(1), FT(1e8))
MK.lines!(
Expand Down
12 changes: 8 additions & 4 deletions docs/src/plots/MarshallPalmer_distribution.jl
Original file line number Diff line number Diff line change
Expand Up @@ -36,10 +36,14 @@ qᵣ_3 = FT(1.0e-3)

r_range = range(25e-6, stop = 1e-2, length = 1000)

n_r_0 = [Marshall_Palmer_distribution(CMP.Rain(FT), qᵣ_0, ρ, r) for r in r_range]
n_r_1 = [Marshall_Palmer_distribution(CMP.Rain(FT), qᵣ_1, ρ, r) for r in r_range]
n_r_2 = [Marshall_Palmer_distribution(CMP.Rain(FT), qᵣ_2, ρ, r) for r in r_range]
n_r_3 = [Marshall_Palmer_distribution(CMP.Rain(FT), qᵣ_3, ρ, r) for r in r_range]
n_r_0 =
[Marshall_Palmer_distribution(CMP.Rain(FT), qᵣ_0, ρ, r) for r in r_range]
n_r_1 =
[Marshall_Palmer_distribution(CMP.Rain(FT), qᵣ_1, ρ, r) for r in r_range]
n_r_2 =
[Marshall_Palmer_distribution(CMP.Rain(FT), qᵣ_2, ρ, r) for r in r_range]
n_r_3 =
[Marshall_Palmer_distribution(CMP.Rain(FT), qᵣ_3, ρ, r) for r in r_range]

fig = MK.Figure(resolution = (1100, 600))
ax1 = MK.Axis(
Expand Down
Loading

0 comments on commit 1229da4

Please sign in to comment.