Skip to content

Commit

Permalink
moving code around
Browse files Browse the repository at this point in the history
  • Loading branch information
@amylu00
amylu00 committed Aug 31, 2023
1 parent 8659459 commit 47effef
Show file tree
Hide file tree
Showing 3 changed files with 94 additions and 61 deletions.
54 changes: 2 additions & 52 deletions docs/src/IceNucleation.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -117,57 +117,7 @@ This parameterization is valid only when ``0.26 < \Delta a_w < 0.36`` and ``185K

## ABIFM Example Figures
```@example
import Plots
import CloudMicrophysics
import CLIMAParameters
import Thermodynamics
const PL = Plots
const IN = CloudMicrophysics.HetIceNucleation
const CMP = CloudMicrophysics.Parameters
const CT = CloudMicrophysics.CommonTypes
const CO = CloudMicrophysics.Common
const CP = CLIMAParameters
const TD = Thermodynamics
include(joinpath(pkgdir(CloudMicrophysics), "test", "create_parameters.jl"))
FT = Float64
toml_dict = CP.create_toml_dict(FT; dict_type = "alias")
const prs = cloud_microphysics_parameters(toml_dict)
thermo_params = CMP.thermodynamics_params(prs)
# Initializing
temp = collect(210.0:2:232.0) # air temperature
x = 0.1 # wt% sulphuric acid in droplets
Delta_a = Vector{Float64}(undef, length(temp))
J = Vector{Float64}(undef, length(temp))
# Knopf and Alpert 2013 Figure 4A
# https://doi.org/10.1039/C3FD00035D
dust_type = CT.KaoliniteType()
it = 1
for T in temp
Delta_a[it] = CO.Delta_a_w(prs, x, T)
J[it] = IN.ABIFM_J(prs, dust_type, Delta_a[it])
global it += 1
end
log10J_converted = @. log10(J*1e-4)
# data read from Fig 4 in Knopf & Alpert 2013
# using https://automeris.io/WebPlotDigitizer/
KA13_Delta_a_obs = [0.13641, 0.16205, 0.21538, 0.23897, 0.24513, 0.24718, 0.25026, 0.25128, 0.25231, 0.25333, 0.25538, 0.25744, 0.25846, 0.25949, 0.26051, 0.26051, 0.26462, 0.26462, 0.26872, 0.26974, 0.27077, 0.27077, 0.27179, 0.27385, 0.27692, 0.27795, 0.27795, 0.27795, 0.28308, 0.28410, 0.28410, 0.28615, 0.28718, 0.28718, 0.29128, 0.29128, 0.29231, 0.29333, 0.29744, 0.29744, 0.29744, 0.29949, 0.30359, 0.30462, 0.30564, 0.30667, 0.31077, 0.31077, 0.31077]
KA13_log10J_obs = [-3.51880, -3.20301, 2.21053, 2.57143, 2.25564, 3.56391, 3.20301, 2.25564, 3.78947, 4.42105, 3.51880, 2.84211, 4.15038, 3.24812, 3.78947, 4.37594, 3.38346, 4.46617, 4.06015, 4.73684, 4.06015, 3.60902, 6.13534, 4.51128, 4.37594, 4.82707, 4.96241, 5.23308, 3.92481, 5.36842, 5.63910, 5.81955, 4.60150, 4.96241, 5.50376, 6.00000, 5.14286, 5.77444, 5.41353, 6.09023, 5.77444, 5.14286, 6.18045, 5.86466, 5.54887, 5.27820, 6.09023, 5.77444, 5.54887]
KA13_Delta_a_param = [0.10256, 0.35692, 0.21949]
KA13_log10J_param = [-4.91729, 8.97744, 1.44361]
PL.plot(Delta_a, log10J_converted, label="CliMA", xlabel="Delta a_w [unitless]", ylabel="log10(J) [cm^-2 s^-1]")
PL.scatter!(KA13_Delta_a_obs, KA13_log10J_obs, markercolor = :black, label="paper observations")
PL.plot!(KA13_Delta_a_param, KA13_log10J_param, linecolor = :red, label="paper parameterization")
PL.savefig("Knopf_Alpert_fig_1.svg")
include("ice_nucleation_plots/KnopfAlpert2013_fig1.jl")
```
![](Knopf_Alpert_fig_1.svg)

Expand All @@ -176,6 +126,7 @@ The following plot shows J as a function of temperature as compared to figure 5a
```@example
include("ice_nucleation_plots/KnopfAlpert2013_fig5.jl")
```
![](KnopfAlpert2013_fig5.svg)
Note that water activity of the droplet was assumed equal to relative humidity so that:
```math
\begin{equation}
Expand All @@ -187,4 +138,3 @@ where `T_dew` is the dewpoint (in this example it is -45C).
It is also important to note that this plot is reflective of cirrus clouds
and shows only a very small temperature range. The two curves are slightly
off because of small differences in parameterizations for vapor pressures.

51 changes: 51 additions & 0 deletions docs/src/ice_nucleation_plots/KnopfAlpert2013_fig1.jl
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,51 @@
import Plots

import CloudMicrophysics
import CLIMAParameters
import Thermodynamics

const PL = Plots
const IN = CloudMicrophysics.HetIceNucleation
const CMP = CloudMicrophysics.Parameters
const CT = CloudMicrophysics.CommonTypes
const CO = CloudMicrophysics.Common
const CP = CLIMAParameters
const TD = Thermodynamics

include(joinpath(pkgdir(CloudMicrophysics), "test", "create_parameters.jl"))
FT = Float64
toml_dict = CP.create_toml_dict(FT; dict_type = "alias")
const prs = cloud_microphysics_parameters(toml_dict)
thermo_params = CMP.thermodynamics_params(prs)

# Initializing
temp = collect(210.0:2:232.0) # air temperature
x = 0.1 # wt% sulphuric acid in droplets
Delta_a = Vector{Float64}(undef, length(temp))
J = Vector{Float64}(undef, length(temp))

# Knopf and Alpert 2013 Figure 4A
# https://doi.org/10.1039/C3FD00035D
dust_type = CT.KaoliniteType()

it = 1
for T in temp
Delta_a[it] = CO.Delta_a_w(prs, x, T)
J[it] = IN.ABIFM_J(prs, dust_type, Delta_a[it])
global it += 1
end
log10J_converted = @. log10(J*1e-4)

# data read from Fig 4 in Knopf & Alpert 2013
# using https://automeris.io/WebPlotDigitizer/
KA13_Delta_a_obs = [0.13641, 0.16205, 0.21538, 0.23897, 0.24513, 0.24718, 0.25026, 0.25128, 0.25231, 0.25333, 0.25538, 0.25744, 0.25846, 0.25949, 0.26051, 0.26051, 0.26462, 0.26462, 0.26872, 0.26974, 0.27077, 0.27077, 0.27179, 0.27385, 0.27692, 0.27795, 0.27795, 0.27795, 0.28308, 0.28410, 0.28410, 0.28615, 0.28718, 0.28718, 0.29128, 0.29128, 0.29231, 0.29333, 0.29744, 0.29744, 0.29744, 0.29949, 0.30359, 0.30462, 0.30564, 0.30667, 0.31077, 0.31077, 0.31077]
KA13_log10J_obs = [-3.51880, -3.20301, 2.21053, 2.57143, 2.25564, 3.56391, 3.20301, 2.25564, 3.78947, 4.42105, 3.51880, 2.84211, 4.15038, 3.24812, 3.78947, 4.37594, 3.38346, 4.46617, 4.06015, 4.73684, 4.06015, 3.60902, 6.13534, 4.51128, 4.37594, 4.82707, 4.96241, 5.23308, 3.92481, 5.36842, 5.63910, 5.81955, 4.60150, 4.96241, 5.50376, 6.00000, 5.14286, 5.77444, 5.41353, 6.09023, 5.77444, 5.14286, 6.18045, 5.86466, 5.54887, 5.27820, 6.09023, 5.77444, 5.54887]

KA13_Delta_a_param = [0.10256, 0.35692, 0.21949]
KA13_log10J_param = [-4.91729, 8.97744, 1.44361]

PL.plot(Delta_a, log10J_converted, label="CliMA", xlabel="Delta a_w [unitless]", ylabel="log10(J) [cm^-2 s^-1]")
PL.scatter!(KA13_Delta_a_obs, KA13_log10J_obs, markercolor = :black, label="paper observations")
PL.plot!(KA13_Delta_a_param, KA13_log10J_param, linecolor = :red, label="paper parameterization")

PL.savefig("Knopf_Alpert_fig_1.svg")
50 changes: 41 additions & 9 deletions docs/src/ice_nucleation_plots/KnopfAlpert2013_fig5.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -18,8 +18,25 @@ thermo_params = CMP.thermodynamics_params(prs)

# 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]
KA13_fig5A_J = [70170382.86704, 12101528.74384, 1277935.32665, 52710.05764, 6040.39151, 293.39697, 18.97171, 4.18121]
temp = [228.20357,
228.33571,
228.50357,
228.75000,
228.92143,
229.16786,
229.39643,
229.52143,
]
KA13_fig5A_J = [
70170382.86704,
12101528.74384,
1277935.32665,
52710.05764,
6040.39151,
293.39697,
18.97171,
4.18121,
]

# Our parameterization
dust_type = CMT.IlliteType()
Expand All @@ -30,21 +47,36 @@ J_ABIFM = Vector{Float64}(undef, length(temp))
it = 1
for T in temp

a_sol = TD.saturation_vapor_pressure(thermo_params, T_dew, TD.Liquid()) / TD.saturation_vapor_pressure(thermo_params, T, TD.Liquid())
a_ice = TD.saturation_vapor_pressure(thermo_params, T, TD.Ice()) / TD.saturation_vapor_pressure(thermo_params, T, TD.Liquid())
a_sol =
TD.saturation_vapor_pressure(thermo_params, T_dew, TD.Liquid()) /
TD.saturation_vapor_pressure(thermo_params, T, TD.Liquid())
a_ice =
TD.saturation_vapor_pressure(thermo_params, T, TD.Ice()) /
TD.saturation_vapor_pressure(thermo_params, T, TD.Liquid())
Δa_w = max(abs(a_sol - a_ice), FT(0.))
J_ABIFM[it] = CMI.ABIFM_J(prs, dust_type, Δa_w) * 1e-4
J_ABIFM[it] = CMI.ABIFM_J(prs, dust_type, Δa_w) * 1e-4 # converted from SI units to cm^-2 s^-1

global it += 1
end

# Plot results
fig = MK.Figure(resolution = (1000, 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!(ax1, temp, KA13_fig5A_J, label = "KA13 Figure 5A", linestyle = :dash, color = :green)
MK.lines!(ax1,
temp,
KA13_fig5A_J,
label = "KA13 Figure 5A",
linestyle = :dash,
color = :green,
)
MK.lines!(ax1, temp, J_ABIFM, label = "CliMA", color = :green)
fig[1,2] = MK.Legend(fig, ax1)
fig[1, 2] = MK.Legend(fig, ax1)

MK.save("ABIFM_cirrus.svg", fig)
MK.save("KnopfAlpert2013_fig5.svg", fig)

0 comments on commit 47effef

Please sign in to comment.