Merge pull request #312 from CliMA/aj/plots_docs_cleanup

dont show plot examples code in docs
CliMA · Feb 5, 2024 · 0d680b0 · 0d680b0
2 parents 9e38a97 + 0266d43
commit 0d680b0
Show file tree

Hide file tree

Showing 9 changed files with 1,022 additions and 479 deletions.
diff --git a/docs/bibliography.bib b/docs/bibliography.bib
@@ -190,7 +190,7 @@ @article{Jensen2022
 number = {17},
 pages = {e2022JD036535},
 keywords = {cirrus, nucleation, freezing},
-doi = {https://doi.org/10.1029/2022JD036535},
+doi = {10.1029/2022JD036535},
 url = {https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/2022JD036535},
 note = {e2022JD036535 2022JD036535},
 year = {2022}
@@ -616,7 +616,7 @@ @article{Vehkamaki2002
   volume = {107},
   number = {D22},
   year = {2002},
-  doi = {110.1029/2002JD002184}
+  doi = {10.1029/2002JD002184}
 }
 
 @article{Wang2006,

diff --git a/docs/src/AerosolActivation.md b/docs/src/AerosolActivation.md
@@ -257,104 +257,16 @@ where:
   - ``S_{ci}`` is the mode critical supersaturation,
   - ``S_{max}`` is the maximum supersaturation.
 
-## Example figures
-
-```@example
-import Plots
-
-import CloudMicrophysics
-import CLIMAParameters
-import Thermodynamics
-
-const PL = Plots
-const AM = CloudMicrophysics.AerosolModel
-const AA = CloudMicrophysics.AerosolActivation
-const CMP = CloudMicrophysics.Parameters
-const TD = Thermodynamics
-
-FT = Float64
-tps = Thermodynamics.Parameters.ThermodynamicsParameters(FT)
-aip = CMP.AirProperties(FT)
-ap = CMP.AerosolActivationParameters(FT)
-
-# Atmospheric conditions
-T = 294.0         # air temperature
-p = 1000.0 *1e2   # air pressure
-w = 0.5           # vertical velocity
-
-# We need the phase partition here only so that we can compute the
-# moist R_m and cp_m in aerosol activation module.
-# We are assuming here saturated conditions and no liquid water or ice.
-# This is consistent with the assumptions of the aerosol activation scheme.
-p_vs = TD.saturation_vapor_pressure(tps, T, TD.Liquid())
-q_vs = 1 / (1 - TD.Parameters.molmass_ratio(tps) * (p_vs - p) / p_vs)
-q = TD.PhasePartition(q_vs, 0.0, 0.0)
-
-# Abdul-Razzak and Ghan 2000 Figure 1 mode 1
-# https://doi.org/10.1029/1999JD901161
-r_dry = 0.05 * 1e-6 # um
-stdev = 2.0         # -
-N_1 = 100.0 * 1e6   # 1/m3
-
-# Sulfate - universal parameters
-sulfate = CMP.Sulfate(FT)
-
-n_components_1 = 1
-mass_fractions_1 = (1.0,)
-paper_mode_1_B = AM.Mode_B(
-    r_dry,
-    stdev,
-    N_1,
-    mass_fractions_1,
-    (sulfate.ϵ,),
-    (sulfate.ϕ,),
-    (sulfate.M,),
-    (sulfate.ν,),
-    (sulfate.ρ,),
-    n_components_1,
-)
-
-N_2_range = range(0, stop=5000 * 1e6, length=100)
-N_act_frac_B = Vector{Float64}(undef, 100)
-
-it = 1
-for N_2 in N_2_range
-        n_components_2 = 1
-        mass_fractions_2 = (1.0,)
-        paper_mode_2_B = AM.Mode_B(
-          r_dry,
-          stdev,
-          N_2,
-          mass_fractions_2,
-          (sulfate.ϵ,),
-          (sulfate.ϕ,),
-          (sulfate.M,),
-          (sulfate.ν,),
-          (sulfate.ρ,),
-          n_components_2,
-        )
-        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
-        global it += 1
-end
-
-# data read from Fig 1 in Abdul-Razzak and Ghan 2000
-# using https://automeris.io/WebPlotDigitizer/
-include("plots/ARGdata.jl")
-
-PL.plot(N_2_range * 1e-6, N_act_frac_B, label="CliMA-B", 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.savefig("Abdul-Razzak_and_Ghan_fig_1.svg")
-```
-![](Abdul-Razzak_and_Ghan_fig_1.svg)
-
 ## Example of Aerosol Activation from ARG2000
 
 The figures below have been reproduced from [Abdul-RazzakandGhan2000](@cite)
 using the aerosol activation parameterization implemented in this project.
 
+```@example
+include("plots/ARGplots_fig1.jl")
+```
+![](Abdul-Razzak_and_Ghan_fig_1.svg)
+
 ```@example
 include("plots/ARGplots.jl")
 make_ARG_figX(2)

diff --git a/docs/src/Microphysics1M.md b/docs/src/Microphysics1M.md
@@ -751,173 +751,8 @@ If ``T > T_{freeze}``:
 ```
 ## Example figures
 
-```@example example_figures
-import Plots
-
-import Thermodynamics
-import CloudMicrophysics
-import CLIMAParameters
-
-const PL = Plots
-const CM1 = CloudMicrophysics.Microphysics1M
-const TD = Thermodynamics
-const CMP = CloudMicrophysics.Parameters
-
-FT = Float64
-
-const tps = Thermodynamics.Parameters.ThermodynamicsParameters(FT)
-const aps = CMP.AirProperties(FT)
-const liquid = CMP.CloudLiquid(FT)
-const ice = CMP.CloudIce(FT)
-const rain = CMP.Rain(FT)
-const snow = CMP.Snow(FT)
-const Chen2022 = CMP.Chen2022VelType(FT)
-const Blk1MVel = CMP.Blk1MVelType(FT)
-const ce = CMP.CollisionEff(FT)
-
-# eq. 5d in [Grabowski1996](@cite)
-function terminal_velocity_empirical(q_rai::DT, q_tot::DT, ρ::DT, ρ_air_ground::DT) where {DT<:Real}
-    rr  = q_rai / (DT(1) - q_tot)
-    vel = DT(14.34) * ρ_air_ground^DT(0.5) * ρ^-DT(0.3654) * rr^DT(0.1346)
-    return vel
-end
-
-# eq. 5b in [Grabowski1996](@cite)
-function accretion_empirical(q_rai::DT, q_liq::DT, q_tot::DT) where {DT<:Real}
-    rr  = q_rai / (DT(1) - q_tot)
-    rl  = q_liq / (DT(1) - q_tot)
-    return DT(2.2) * rl * rr^DT(7/8)
-end
-
-# eq. 5c in [Grabowski1996](@cite)
-function rain_evap_empirical(q_rai::DT, q::TD.PhasePartition, T::DT, p::DT, ρ::DT) where {DT<:Real}
-
-    ts_neq = TD.PhaseNonEquil_ρTq(tps, ρ, T, q)
-    q_sat  = TD.q_vap_saturation(tps, ts_neq)
-    q_vap  = q.tot - q.liq
-    rr     = q_rai / (DT(1) - q.tot)
-    rv_sat = q_sat / (DT(1) - q.tot)
-    S      = q_vap/q_sat - DT(1)
-
-    ag, bg = 5.4 * 1e2, 2.55 * 1e5
-    G = DT(1) / (ag + bg / p / rv_sat) / ρ
-
-    av, bv = 1.6, 124.9
-    F = av * (ρ/DT(1e3))^DT(0.525)  * rr^DT(0.525) + bv * (ρ/DT(1e3))^DT(0.7296) * rr^DT(0.7296)
-
-    return DT(1) / (DT(1) - q.tot) * S * F * G
-end
-
-# example values
-q_liq_range  = range(1e-8, stop=5e-3, length=100)
-q_ice_range  = range(1e-8, stop=5e-3, length=100)
-q_rain_range = range(1e-8, stop=5e-3, length=100)
-q_snow_range = range(1e-8, stop=5e-3, length=100)
-ρ_air, ρ_air_ground = 1.2, 1.22
-q_liq, q_ice, q_tot = 5e-4, 5e-4, 20e-3
-
-PL.plot( q_rain_range * 1e3, [CM1.terminal_velocity(rain, Blk1MVel.rain, ρ_air, q_rai) for q_rai in q_rain_range], linewidth=3, label="Rain-1M-default", color=:blue, xlabel="q_rain/ice/snow [g/kg]", ylabel="terminal velocity [m/s]")
-PL.plot!(q_rain_range * 1e3, [CM1.terminal_velocity(rain, Chen2022.rain, ρ_air, q_rai) for q_rai in q_rain_range], linewidth=3, label="Rain-Chen",       color=:blue, linestyle=:dot)
-PL.plot!(q_rain_range * 1e3, [terminal_velocity_empirical(q_rai, q_tot, ρ_air, ρ_air_ground) for q_rai in q_rain_range],        linewidth=3, label="Rain-Empirical",  color=:green)
-PL.plot!(q_ice_range  * 1e3, [CM1.terminal_velocity(ice,  Chen2022.snow_ice, ρ_air, q_ice) for q_ice in q_ice_range],  linewidth=3, label="Ice-Chen",        color=:pink, linestyle=:dot)
-PL.plot!(q_snow_range * 1e3, [CM1.terminal_velocity(snow, Blk1MVel.snow, ρ_air, q_sno) for q_sno in q_snow_range], linewidth=3, label="Snow-1M-default", color=:red)
-PL.plot!(q_snow_range * 1e3, [CM1.terminal_velocity(snow, Chen2022.snow_ice, ρ_air, q_sno) for q_sno in q_snow_range], linewidth=3, label="Snow-Chen",       color=:red, linestyle=:dot)
-PL.savefig("terminal_velocity.svg") # hide
-
-T = 273.15
-PL.plot( q_liq_range * 1e3, [CM1.conv_q_liq_to_q_rai(rain.acnv1M, q_liq) for q_liq in q_liq_range], linewidth=3, xlabel="q_liq or q_ice [g/kg]", ylabel="autoconversion rate [1/s]", label="Rain")
-PL.plot!(q_ice_range * 1e3, [CM1.conv_q_ice_to_q_sno(ice, aps, tps, TD.PhasePartition(q_tot, 0., q_ice), ρ_air, T-5)  for q_ice in q_ice_range], linewidth=3, label="Snow T=-5C")
-PL.plot!(q_ice_range * 1e3, [CM1.conv_q_ice_to_q_sno(ice, aps, tps, TD.PhasePartition(q_tot, 0., q_ice), ρ_air, T-10) for q_ice in q_ice_range], linewidth=3, label="Snow T=-10C")
-PL.plot!(q_ice_range * 1e3, [CM1.conv_q_ice_to_q_sno(ice, aps, tps, TD.PhasePartition(q_tot, 0., q_ice), ρ_air, T-15) for q_ice in q_ice_range], linewidth=3, label="Snow T=-15C")
-PL.savefig("autoconversion_rate.svg") # hide
-
-PL.plot( q_rain_range * 1e3, [CM1.accretion(liquid, rain, Blk1MVel.rain, ce, q_liq, q_rai, ρ_air) for q_rai in q_rain_range], linewidth=3, xlabel="q_rain or q_snow [g/kg]", ylabel="accretion rate [1/s]", label="Liq+Rain-CLIMA")
-PL.plot!(q_rain_range * 1e3, [CM1.accretion(ice,    rain, Blk1MVel.rain, ce, q_ice, q_rai, ρ_air) for q_rai in q_rain_range], linewidth=3, label="Ice+Rain-CLIMA")
-PL.plot!(q_snow_range * 1e3, [CM1.accretion(liquid, snow, Blk1MVel.snow, ce, q_liq, q_sno, ρ_air) for q_sno in q_snow_range], linewidth=3, label="Liq+Snow-CLIMA")
-PL.plot!(q_snow_range * 1e3, [CM1.accretion(ice,    snow, Blk1MVel.snow, ce, q_ice, q_sno, ρ_air) for q_sno in q_snow_range], linewidth=4, linestyle=:dash, label="Ice+Snow-CLIMA")
-PL.plot!(q_rain_range * 1e3, [accretion_empirical(q_rai, q_liq, q_tot) for q_rai in q_rain_range], linewidth=3, label="Liq+Rain-Empirical")
-PL.savefig("accretion_rate.svg") # hide
-
-q_ice = 1e-6
-PL.plot(q_rain_range * 1e3, [CM1.accretion_rain_sink(rain, ice, Blk1MVel.rain, ce, q_ice, q_rai, ρ_air) for q_rai in q_rain_range], linewidth=3, xlabel="q_rain or q_snow [g/kg]", ylabel="accretion rain sink rate [1/s]", label="q_ice=1e-6")
-q_ice = 1e-5
-PL.plot!(q_rain_range * 1e3, [CM1.accretion_rain_sink(rain, ice, Blk1MVel.rain, ce, q_ice, q_rai, ρ_air) for q_rai in q_rain_range], linewidth=3, xlabel="q_rain or q_snow [g/kg]", ylabel="accretion rain sink rate [1/s]", label="q_ice=1e-5")
-q_ice = 1e-4
-PL.plot!(q_rain_range * 1e3, [CM1.accretion_rain_sink(rain, ice, Blk1MVel.rain, ce, q_ice, q_rai, ρ_air) for q_rai in q_rain_range], linewidth=3, xlabel="q_rain or q_snow [g/kg]", ylabel="accretion rain sink rate [1/s]", label="q_ice=1e-4")
-PL.savefig("accretion_rain_sink_rate.svg") # hide
-
-q_sno = 1e-6
-PL.plot(q_rain_range * 1e3, [CM1.accretion_snow_rain(rain, snow, Blk1MVel.rain, Blk1MVel.snow, ce, q_rai, q_sno, ρ_air) for q_rai in q_rain_range], linewidth=3, xlabel="q_rain [g/kg]", ylabel="snow-rain accretion rate [1/s] T>0", label="q_snow=1e-6")
-q_sno = 1e-5
-PL.plot!(q_rain_range * 1e3, [CM1.accretion_snow_rain(rain, snow, Blk1MVel.rain, Blk1MVel.snow, ce, q_rai, q_sno, ρ_air) for q_rai in q_rain_range], linewidth=3, label="q_snow=1e-5")
-q_sno = 1e-4
-PL.plot!(q_rain_range * 1e3, [CM1.accretion_snow_rain(rain, snow, Blk1MVel.rain, Blk1MVel.snow, ce, q_rai, q_sno, ρ_air) for q_rai in q_rain_range], linewidth=3, label="q_snow=1e-4")
-PL.savefig("accretion_snow_rain_above_freeze.svg") # hide
-
-q_rai = 1e-6
-PL.plot(q_snow_range * 1e3, [CM1.accretion_snow_rain(snow, rain, Blk1MVel.snow, Blk1MVel.rain, ce, q_sno, q_rai, ρ_air) for q_sno in q_snow_range], linewidth=3, xlabel="q_snow [g/kg]", ylabel="snow-rain accretion rate [1/s] T<0", label="q_rain=1e-6")
-q_rai = 1e-5
-PL.plot!(q_snow_range * 1e3, [CM1.accretion_snow_rain(snow, rain, Blk1MVel.snow, Blk1MVel.rain, ce, q_sno, q_rai, ρ_air) for q_sno in q_snow_range], linewidth=3, label="q_rain=1e-5")
-q_rai = 1e-4
-PL.plot!(q_snow_range * 1e3, [CM1.accretion_snow_rain(snow, rain, Blk1MVel.snow, Blk1MVel.rain, ce, q_sno, q_rai, ρ_air) for q_sno in q_snow_range], linewidth=3, label="q_rain=1e-4")
-PL.savefig("accretion_snow_rain_below_freeze.svg") # hide
-
-# example values
-T, p = 273.15 + 15, 90000.
-ϵ = 1. / TD.Parameters.molmass_ratio(tps)
-p_sat = TD.saturation_vapor_pressure(tps, T, TD.Liquid())
-q_sat = ϵ * p_sat / (p + p_sat * (ϵ - 1.))
-q_rain_range = range(1e-8, stop=5e-3, length=100)
-q_tot = 15e-3
-q_vap = 0.15 * q_sat
-q_ice = 0.
-q_liq = q_tot - q_vap - q_ice
-q = TD.PhasePartition(q_tot, q_liq, q_ice)
-R = TD.gas_constant_air(tps, q)
-ρ = p / R / T
-
-PL.plot(q_rain_range * 1e3,  [CM1.evaporation_sublimation(rain, Blk1MVel.rain, aps, tps, q, q_rai, ρ, T) for q_rai in q_rain_range], xlabel="q_rain [g/kg]", linewidth=3, ylabel="rain evaporation rate [1/s]", label="ClimateMachine")
-PL.plot!(q_rain_range * 1e3, [rain_evap_empirical(q_rai, q, T, p, ρ) for q_rai in q_rain_range], linewidth=3, label="empirical")
-PL.savefig("rain_evaporation_rate.svg") # hide
-
-# example values
-T, p = 273.15 - 15, 90000.
-ϵ = 1. / TD.Parameters.molmass_ratio(tps)
-p_sat = TD.saturation_vapor_pressure(tps, T, TD.Ice())
-q_sat = ϵ * p_sat / (p + p_sat * (ϵ - 1.))
-q_snow_range = range(1e-8, stop=5e-3, length=100)
-q_tot = 15e-3
-q_vap = 0.15 * q_sat
-q_liq = 0.
-q_ice = q_tot - q_vap - q_liq
-q = TD.PhasePartition(q_tot, q_liq, q_ice)
-R = TD.gas_constant_air(tps, q)
-ρ = p / R / T
-
-PL.plot(q_snow_range * 1e3,  [CM1.evaporation_sublimation(snow, Blk1MVel.snow, aps, tps, q, q_sno, ρ, T) for q_sno in q_snow_range], xlabel="q_snow [g/kg]", linewidth=3, ylabel="snow deposition sublimation rate [1/s]", label="T<0")
-
-T, p = 273.15 + 15, 90000.
-ϵ = 1. / TD.Parameters.molmass_ratio(tps)
-p_sat = TD.saturation_vapor_pressure(tps, T, TD.Ice())
-q_sat = ϵ * p_sat / (p + p_sat * (ϵ - 1.))
-q_snow_range = range(1e-8, stop=5e-3, length=100)
-q_tot = 15e-3
-q_vap = 0.15 * q_sat
-q_liq = 0.
-q_ice = q_tot - q_vap - q_liq
-q = TD.PhasePartition(q_tot, q_liq, q_ice)
-R = TD.gas_constant_air(tps, q)
-ρ = p / R / T
-
-PL.plot!(q_snow_range * 1e3,  [CM1.evaporation_sublimation(snow, Blk1MVel.snow, aps, tps, q, q_sno, ρ, T) for q_sno in q_snow_range], xlabel="q_snow [g/kg]", linewidth=3, ylabel="snow deposition sublimation rate [1/s]", label="T>0")
-PL.savefig("snow_sublimation_deposition_rate.svg") # hide
-
-T=273.15
-PL.plot( q_snow_range * 1e3,  [CM1.snow_melt(snow, Blk1MVel.snow, aps, tps, q_sno, ρ, T+2) for q_sno in q_snow_range], xlabel="q_snow [g/kg]", linewidth=3, ylabel="snow melt rate [1/s]", label="T=2C")
-PL.plot!(q_snow_range * 1e3,  [CM1.snow_melt(snow, Blk1MVel.snow, aps, tps, q_sno, ρ, T+4) for q_sno in q_snow_range], xlabel="q_snow [g/kg]", linewidth=3, label="T=4C")
-PL.plot!(q_snow_range * 1e3,  [CM1.snow_melt(snow, Blk1MVel.snow, aps, tps, q_sno, ρ, T+6) for q_sno in q_snow_range], xlabel="q_snow [g/kg]", linewidth=3, label="T=6C")
-PL.savefig("snow_melt_rate.svg") # hide
-
+```@example
+include("plots/Microphysics1M_plots.jl")
 ```
 ![](terminal_velocity.svg)
 ![](autoconversion_rate.svg)