diff --git a/docs/src/plots/ARGdata.jl b/docs/src/plots/ARGdata.jl index 7cee34886..a71eb25e7 100644 --- a/docs/src/plots/ARGdata.jl +++ b/docs/src/plots/ARGdata.jl @@ -531,3 +531,19 @@ Fig5b_y_param = [ 0.7030508474576271, 0.7505084745762711, ] + +function get_ARG_data(fig_num) + if fig_num == 1 + return Fig1_x_obs, Fig1_y_obs, Fig1_x_param, Fig1_y_param + elseif fig_num == 2 + return Fig2a_x_obs, Fig2a_y_obs, Fig2a_x_param, Fig2a_y_param + elseif fig_num == 3 + return Fig3a_x_obs, Fig3a_y_obs, Fig3a_x_param, Fig3a_y_param + elseif fig_num == 4 + return Fig4a_x_obs, Fig4a_y_obs, Fig4a_x_param, Fig4a_y_param + elseif fig_num == 5 + return Fig5a_x_obs, Fig5a_y_obs, Fig5a_x_param, Fig5a_y_param + else + error("Figure $fig_num does not exist") + end +end diff --git a/docs/src/plots/ARGplots.jl b/docs/src/plots/ARGplots.jl index e4aab2820..b73af5da8 100644 --- a/docs/src/plots/ARGplots.jl +++ b/docs/src/plots/ARGplots.jl @@ -1,4 +1,4 @@ -import Plots as PL +import CairoMakie as MK import CloudMicrophysics as CM import CloudMicrophysics.AerosolModel as AM @@ -36,23 +36,112 @@ M_insol = 0.044 # molar mass of insol ρ_insol = 1770.0 # density of insol κ_insol = 0.0 # hygroscopicity of insol +all_components_B = ( + (sulfate.ϵ, ϵ_insol), + (sulfate.ϕ, ϕ_insol), + (sulfate.M, M_insol), + (sulfate.ν, ν_insol), + (sulfate.ρ, ρ_insol), +) + +all_components_κ = ( + (sulfate.M, M_insol), + (sulfate.κ, κ_insol), +) + +first_n_values(tup, n) = getindex(tup, 1:n) + +function make_mode_B(r_dry, stdev, N, mass_mixing_ratios) + n_components = length(mass_mixing_ratios) + return AM.Mode_B( + r_dry, stdev, N, mass_mixing_ratios, + first_n_values.(all_components_B, n_components)..., + ) +end + +function make_mode_κ(r_dry, stdev, N, vol_mixing_ratios, mass_mixing_ratios) + n_components = length(mass_mixing_ratios) + return AM.Mode_κ( + r_dry, stdev, N, vol_mixing_ratios, mass_mixing_ratios, + first_n_values.(all_components_κ, n_components)..., + ) +end + function mass2vol(mass_mixing_ratios) - if length(mass_mixing_ratios) == 2 - densities = (sulfate.ρ, ρ_insol) - else - densities = (sulfate.ρ,) - end - volfractions = - (mass_mixing_ratios ./ densities) ./ - sum(mass_mixing_ratios ./ densities) + n_components = length(mass_mixing_ratios) + densities = first_n_values((sulfate.ρ, ρ_insol), n_components) + volfractions = @. mass_mixing_ratios / densities / $sum(mass_mixing_ratios / densities) return volfractions end +function compute_activation_fractions(mode_1, mode_2, w) + AD = AM.AerosolDistribution((mode_1, mode_2)) + act_per_mode = AA.N_activated_per_mode(ap, AD, aip, tps, T, p, w, q_vs, FT(0), FT(0)) + return act_per_mode +end + +# function get_fig_params(X, ) +# @assert X in (1, 2, 3, 4, 5) +# xvar = Dict( +# 1 => :N_2, +# 2 => :N_2, +# 3 => :sulfate_fraction, +# 4 => :r_dry_2, +# 5 => :w, +# )[X] +# mass_mixing_ratios_1 = X in (1, 4) ? (1.0,) : (1.0, 0.0) +# default_params = (; +# mass_mixing_ratios_1, +# r_dry_1 = 0.05 * 1e-6, # um +# stdev_1 = 2.0, # - +# N_1 = 100.0 * 1e6, # 1/m3 +# stdev_2 = 2.0, +# r_dry_2 = 0.05 * 1e-6, # um +# N_2 = 100.0 * 1e6, # 1/m3 +# w = 0.5, # vertical velocity, m/s +# ) + +# if X == 1 +# return (; +# default_params..., +# mass_mixing_ratios_2 = (1.0,), +# ) +# elseif X == 2 +# return (; +# default_params..., +# mass_mixing_ratios_2 = (0.1, 0.9), +# ) +# elseif X == 3 +# return (; +# default_params..., +# mass_mixing_ratios_2 = [(i, 1 - i) for i in range(0.1, stop = 1, length = 100)], +# ) +# elseif X == 4 +# if X in (1, 4) +# return (; +# mass_mixing_ratios_1 = (1.0,), +# w = 0.5, +# r_dry_2 = 0.05 * 1e-6, +# N_2 = range(100, stop = 5000, length = 100) * 1e6, +# mass_mixing_ratios_2 = X == 1 ? (1.0,) : (0.1, 0.9), +# ) +# elseif X in (2, 3, 5) +# return (; +# mass_mixing_ratios_1 = (1.0, 0.0), +# ) +# end +# end + # Abdul-Razzak and Ghan 2000 # https://doi.org/10.1029/1999JD901161 function make_ARG_figX(X) - p1 = PL.plot() - p2 = PL.plot() + # Create figure with subplots + fig = MK.Figure(; size = (800, 600)) + ax1 = MK.Axis(fig[1, 1]; ylabel = "Mode 1", title = "ARG2000 Fig $X") + ax2 = MK.Axis(fig[2, 1]; ylabel = "Mode 2", xlabel = "Mode 2 aerosol number concentration [1/cm3]") + MK.Label(fig[:, 0], "Activation fraction"; rotation = π / 2) + MK.linkaxes!(ax1, ax2) + MK.limits!(ax2, (0, nothing), (0, 1)) for v_B in (true, false) # mode 1 definitions @@ -65,25 +154,13 @@ function make_ARG_figX(X) mass_mixing_ratios_1 = (1.0,) if v_B paper_mode_1 = AM.Mode_B( - r_dry_1, - stdev_1, - N_1, - mass_mixing_ratios_1, - (sulfate.ϵ,), - (sulfate.ϕ,), - (sulfate.M,), - (sulfate.ν,), - (sulfate.ρ,), + r_dry_1, stdev_1, N_1, mass_mixing_ratios_1, + (sulfate.ϵ,), (sulfate.ϕ,), (sulfate.M,), (sulfate.ν,), (sulfate.ρ,), ) else paper_mode_1 = AM.Mode_κ( - r_dry_1, - stdev_1, - N_1, - vol_mixing_ratios_1, - mass_mixing_ratios_1, - (sulfate.M,), - (sulfate.κ,), + r_dry_1, stdev_1, N_1, vol_mixing_ratios_1, mass_mixing_ratios_1, + (sulfate.M,), (sulfate.κ,), ) end end @@ -92,33 +169,14 @@ function make_ARG_figX(X) vol_mixing_ratios_1 = (1.0, 0.0) mass_mixing_ratios_1 = (1.0, 0.0) if v_B - paper_mode_1 = AM.Mode_B( - r_dry_1, - stdev_1, - N_1, - mass_mixing_ratios_1, - (sulfate.ϵ, ϵ_insol), - (sulfate.ϕ, ϕ_insol), - (sulfate.M, M_insol), - (sulfate.ν, ν_insol), - (sulfate.ρ, ρ_insol), - ) + paper_mode_1 = make_mode_B(r_dry_1, stdev_1, N_1, mass_mixing_ratios_1) else - paper_mode_1 = AM.Mode_κ( - r_dry_1, - stdev_1, - N_1, - vol_mixing_ratios_1, - mass_mixing_ratios_1, - (sulfate.M, M_insol), - (sulfate.κ, κ_insol), - ) + paper_mode_1 = make_mode_κ(r_dry_1, stdev_1, N_1, vol_mixing_ratios_1, mass_mixing_ratios_1) end end # parcel and mode 2 definitions len = 100 - global it = 1 act_frac1 = Vector{Float64}(undef, len) act_frac2 = Vector{Float64}(undef, len) @@ -130,49 +188,19 @@ function make_ARG_figX(X) mass_mixing_ratios_2 = (1.0,) # all sulfate vol_mixing_ratios_2 = (1.0,) # all sulfate - for N2i in N_2 + for (it, N2i) in enumerate(N_2) if v_B - paper_mode_2 = AM.Mode_B( - r_dry_2, - stdev_2, - N2i, - mass_mixing_ratios_2, - (sulfate.ϵ,), - (sulfate.ϕ,), - (sulfate.M,), - (sulfate.ν,), - (sulfate.ρ,), - ) + paper_mode_2 = make_mode_B(r_dry_2, stdev_2, N2i, mass_mixing_ratios_2) else - paper_mode_2 = AM.Mode_κ( - r_dry_2, - stdev_2, - N2i, - vol_mixing_ratios_2, - mass_mixing_ratios_2, - (sulfate.M,), - (sulfate.κ,), - ) + paper_mode_2 = make_mode_κ(r_dry_2, stdev_2, N2i, vol_mixing_ratios_2, mass_mixing_ratios_2) 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_vs, FT(0), FT(0))[1] / - N_1 - act_frac2[it] = - AA.N_activated_per_mode(ap, AD, aip, tps, T, p, w, q_vs, FT(0), FT(0))[2] / - N2i - global it += 1 + act_per_mode = AA.N_activated_per_mode(ap, AD, aip, tps, T, p, w, q_vs, FT(0), FT(0)) + act_frac1[it] = act_per_mode[1] / N_1 + act_frac2[it] = act_per_mode[2] / N2i end - x1_obs = Fig1_x_obs - y1_obs = Fig1_y_obs - x1_param = Fig1_x_param - y1_param = Fig1_y_param - - x2_obs = Fig1_x_obs - y2_obs = Fig1_y_obs - x2_param = Fig1_x_param - y2_param = Fig1_y_param + xvar = N_2 * 1e-6 xlabel = "Mode 2 aerosol number concentration [1/cm3]" @@ -183,49 +211,18 @@ function make_ARG_figX(X) mass_mixing_ratios_2 = (0.1, 0.9) # 10% sulfate, 90% insoluble vol_mixing_ratios_2 = mass2vol(mass_mixing_ratios_2) - for N2i in N_2 + for (it, N2i) in enumerate(N_2) if v_B - paper_mode_2 = AM.Mode_B( - r_dry_2, - stdev_2, - N2i, - mass_mixing_ratios_2, - (sulfate.ϵ, ϵ_insol), - (sulfate.ϕ, ϕ_insol), - (sulfate.M, M_insol), - (sulfate.ν, ν_insol), - (sulfate.ρ, ρ_insol), - ) + paper_mode_2 = make_mode_B(r_dry_2, stdev_2, N2i, mass_mixing_ratios_2) else - paper_mode_2 = AM.Mode_κ( - r_dry_2, - stdev_2, - N2i, - vol_mixing_ratios_2, - mass_mixing_ratios_2, - (sulfate.M, M_insol), - (sulfate.κ, κ_insol), - ) + paper_mode_2 = make_mode_κ(r_dry_2, stdev_2, N2i, vol_mixing_ratios_2, mass_mixing_ratios_2) 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_vs, FT(0), FT(0))[1] / - N_1 - act_frac2[it] = - AA.N_activated_per_mode(ap, AD, aip, tps, T, p, w, q_vs, FT(0), FT(0))[2] / - N2i - global it += 1 + act_per_mode = AA.N_activated_per_mode(ap, AD, aip, tps, T, p, w, q_vs, FT(0), FT(0)) + act_frac1[it] = act_per_mode[1] / N_1 + act_frac2[it] = act_per_mode[2] / N2i end - x1_obs = Fig2a_x_obs - y1_obs = Fig2a_y_obs - x1_param = Fig2a_x_param - y1_param = Fig2a_y_param - - x2_obs = Fig2b_x_obs - y2_obs = Fig2b_y_obs - x2_param = Fig2b_x_param - y2_param = Fig2b_y_param xvar = N_2 * 1e-6 xlabel = "Mode 2 aerosol number concentration [1/cm3]" @@ -237,50 +234,19 @@ function make_ARG_figX(X) xvar = range(0.1, stop = 1, length = len) mass_mixing_ratios_2 = [(i, 1 - i) for i in xvar] - for mmr2i in mass_mixing_ratios_2 + for (it, mmr2i) in enumerate(mass_mixing_ratios_2) vmr2i = mass2vol(mmr2i) if v_B - paper_mode_2 = AM.Mode_B( - r_dry_2, - stdev_2, - N_2, - mmr2i, - (sulfate.ϵ, ϵ_insol), - (sulfate.ϕ, ϕ_insol), - (sulfate.M, M_insol), - (sulfate.ν, ν_insol), - (sulfate.ρ, ρ_insol), - ) + paper_mode_2 = make_mode_B(r_dry_2, stdev_2, N_2, mmr2i) else - paper_mode_2 = AM.Mode_κ( - r_dry_2, - stdev_2, - N_2, - vmr2i, - mmr2i, - (sulfate.M, M_insol), - (sulfate.κ, κ_insol), - ) + paper_mode_2 = make_mode_κ(r_dry_2, stdev_2, N_2, vmr2i, mmr2i) 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_vs, FT(0), FT(0))[1] / - N_1 - act_frac2[it] = - AA.N_activated_per_mode(ap, AD, aip, tps, T, p, w, q_vs, FT(0), FT(0))[2] / - N_2 - global it += 1 + act_per_mode = AA.N_activated_per_mode(ap, AD, aip, tps, T, p, w, q_vs, FT(0), FT(0)) + act_frac1[it] = act_per_mode[1] / N_1 + act_frac2[it] = act_per_mode[2] / N_2 end - x1_obs = Fig3a_x_obs - y1_obs = Fig3a_y_obs - x1_param = Fig3a_x_param - y1_param = Fig3a_y_param - - x2_obs = Fig3b_x_obs - y2_obs = Fig3b_y_obs - x2_param = Fig3b_x_param - y2_param = Fig3b_y_param xlabel = "Mode 2 soluble fraction (fraction of sulfate)" elseif X == 4 @@ -290,49 +256,18 @@ function make_ARG_figX(X) mass_mixing_ratios_2 = (1.0,) # all sulfate vol_mixing_ratios_2 = mass2vol(mass_mixing_ratios_2) - for rd2i in r_dry_2 + for (it, rd2i) in enumerate(r_dry_2) if v_B - paper_mode_2 = AM.Mode_B( - rd2i, - stdev_2, - N_2, - mass_mixing_ratios_2, - (sulfate.ϵ,), - (sulfate.ϕ,), - (sulfate.M,), - (sulfate.ν,), - (sulfate.ρ,), - ) + paper_mode_2 = make_mode_B(rd2i, stdev_2, N_2, mass_mixing_ratios_2) else - paper_mode_2 = AM.Mode_κ( - rd2i, - stdev_2, - N_2, - vol_mixing_ratios_2, - mass_mixing_ratios_2, - (sulfate.M,), - (sulfate.κ,), - ) + paper_mode_2 = make_mode_κ(rd2i, stdev_2, N_2, vol_mixing_ratios_2, mass_mixing_ratios_2) 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_vs, FT(0), FT(0))[1] / - N_1 - act_frac2[it] = - AA.N_activated_per_mode(ap, AD, aip, tps, T, p, w, q_vs, FT(0), FT(0))[2] / - N_2 - global it += 1 + act_per_mode = AA.N_activated_per_mode(ap, AD, aip, tps, T, p, w, q_vs, FT(0), FT(0)) + act_frac1[it] = act_per_mode[1] / N_1 + act_frac2[it] = act_per_mode[2] / N_2 end - x1_obs = Fig4a_x_obs - y1_obs = Fig4a_y_obs - x1_param = Fig4a_x_param - y1_param = Fig4a_y_param - - x2_obs = Fig4b_x_obs - y2_obs = Fig4b_y_obs - x2_param = Fig4b_x_param - y2_param = Fig4b_y_param xvar = r_dry_2 * 1e6 xlabel = "Mode 2 mean radius [um]" @@ -343,93 +278,42 @@ function make_ARG_figX(X) mass_mixing_ratios_2 = (0.1, 0.9) # 10% sulfate, 90% insoluble vol_mixing_ratios_2 = mass2vol(mass_mixing_ratios_2) - if v_B - paper_mode_2 = AM.Mode_B( - r_dry_2, - stdev_2, - N_2, - mass_mixing_ratios_2, - (sulfate.ϵ, ϵ_insol), - (sulfate.ϕ, ϕ_insol), - (sulfate.M, M_insol), - (sulfate.ν, ν_insol), - (sulfate.ρ, ρ_insol), - ) - else - paper_mode_2 = AM.Mode_κ( - r_dry_2, - stdev_2, - N_2, - vol_mixing_ratios_2, - mass_mixing_ratios_2, - (sulfate.M, M_insol), - (sulfate.κ, κ_insol), - ) - end + paper_mode_2_B = make_mode_B(r_dry_2, stdev_2, N_2, mass_mixing_ratios_2) + paper_mode_2_κ = make_mode_κ(r_dry_2, stdev_2, N_2, vol_mixing_ratios_2, mass_mixing_ratios_2) AD = AM.AerosolDistribution((paper_mode_1, paper_mode_2)) - for wi in w - act_frac1[it] = - AA.N_activated_per_mode(ap, AD, aip, tps, T, p, wi, q_vs, FT(0), FT(0))[1] / - N_1 - act_frac2[it] = - AA.N_activated_per_mode(ap, AD, aip, tps, T, p, wi, q_vs, FT(0), FT(0))[2] / - N_2 - global it += 1 + for (it, wi) in enumerate(w) + act_per_mode = AA.N_activated_per_mode(ap, AD, aip, tps, T, p, wi, q_vs, FT(0), FT(0)) + act_frac1[it] = act_per_mode[1] / N_1 + act_frac2[it] = act_per_mode[2] / N_2 end - x1_obs = Fig5a_x_obs - y1_obs = Fig5a_y_obs - x1_param = Fig5a_x_param - y1_param = Fig5a_y_param - - x2_obs = Fig5b_x_obs - y2_obs = Fig5b_y_obs - x2_param = Fig5b_x_param - y2_param = Fig5b_y_param xvar = w xlabel = "Vertical velocity, w [m/s]" end + + # Update xlabel for the second subplot + ax2.xlabel[] = xlabel + v_B ? label = "CliMA-B" : label = "CliMA-κ" - PL.plot!( - p1, - xvar, - act_frac1, - label = label, - ylim = [0, 1], - ylabel = "Mode 1 act frac", - title = "ARG2000 Fig " * string(X), - ) - PL.plot!( - p2, - xvar, - act_frac2, - legend = false, - ylim = [0, 1], - xlabel = xlabel, - ylabel = "Mode 2 act frac", - ) - if v_B == false - PL.scatter!( - p1, - x1_obs, - y1_obs, - markercolor = :black, - label = "ARG2000 observations", - ) - PL.plot!( - p1, - x1_param, - y1_param, - linecolor = :black, - label = "ARG2000 parameterization", - ) - PL.scatter!(p2, x2_obs, y2_obs, markercolor = :black) - PL.plot!(p2, x2_param, y2_param, linecolor = :black) - end + + # Plot lines + MK.lines!(ax1, xvar, act_frac1, label = label) + MK.lines!(ax2, xvar, act_frac2) end - PL.plot(p1, p2, layout = (2, 1)) - PL.savefig("Abdul-Razzak_and_Ghan_fig_" * string(X) * ".svg") + x_obs, y_obs, x_param, y_param = get_ARG_data(X) + MK.scatter!(ax1, x_obs, y_obs, color = :black, label = "ARG2000 observations") + MK.lines!(ax1, x_param, y_param, color = :black, label = "ARG2000 parameterization") + MK.scatter!(ax2, x_obs, y_obs, color = :black) + MK.lines!(ax2, x_param, y_param, color = :black) + + # Add legend to the first subplot + MK.axislegend(ax1, position = :rt) + + # Save the figure + # MK.save("Abdul-Razzak_and_Ghan_fig_$X.svg", fig) + + return fig end