diff --git a/docs/bibliography.bib b/docs/bibliography.bib index 22b52ccee..e74090b0e 100644 --- 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 index 50dd7a144..a753f0da4 100644 --- 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 index 42f7de27e..8e7e18daf 100644 --- 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) diff --git a/docs/src/Microphysics2M.md b/docs/src/Microphysics2M.md index e6ddb38fc..0c96be287 100644 --- a/docs/src/Microphysics2M.md +++ b/docs/src/Microphysics2M.md @@ -633,136 +633,7 @@ and the default free parameter values are: ## Example figures -```@example example_figures - -using CairoMakie -CairoMakie.activate!(type = "svg") - -import CLIMAParameters -import Thermodynamics as TD -import CloudMicrophysics -import CloudMicrophysics.Microphysics1M as CM1 -import CloudMicrophysics.Microphysics2M as CM2 -import CloudMicrophysics.Parameters as CMP - -FT = Float64 - -const tps = TD.Parameters.ThermodynamicsParameters(FT) -const aps = CMP.AirProperties(FT) - -const KK2000 = CMP.KK2000(FT) -const B1994 = CMP.B1994(FT) -const TC1980 = CMP.TC1980(FT) -const LD2004 = CMP.LD2004(FT) -const VarTSc = CMP.VarTimescaleAcnv(FT) -const SB2006 = CMP.SB2006(FT) - -const ce = CMP.CollisionEff(FT) - -const liquid = CMP.CloudLiquid(FT) -const rain = CMP.Rain(FT) -const blk1mvel = CMP.Blk1MVelType(FT) - -include(joinpath(pkgdir(CloudMicrophysics), "docs", "src", "plots", "Wooddata.jl")) - -# Example values -q_liq_range = range(1e-8, stop=1e-3, length=1000) -q_rai_range = range(1e-8, stop=1e-3, length=1000) -N_d_range = range(1e7, stop=1e9, length=1000) -q_liq = 5e-4 -q_rai = 5e-4 -ρ_air = 1.0 # kg m^-3 - -q_liq_KK2000 = [CM2.conv_q_liq_to_q_rai(KK2000, q_liq, ρ_air, N_d = 1e8) for q_liq in q_liq_range] -q_liq_B1994 = [CM2.conv_q_liq_to_q_rai(B1994, q_liq, ρ_air, N_d = 1e8) for q_liq in q_liq_range] -q_liq_TC1980 = [CM2.conv_q_liq_to_q_rai(TC1980, q_liq, ρ_air, N_d = 1e8) for q_liq in q_liq_range] -q_liq_LD2004 = [CM2.conv_q_liq_to_q_rai(LD2004, q_liq, ρ_air, N_d = 1e8) for q_liq in q_liq_range] -q_liq_VarTimeScaleAcnv = [CM2.conv_q_liq_to_q_rai(VarTSc, q_liq, ρ_air, N_d = 1e8) for q_liq in q_liq_range] -q_liq_SB2006 = [CM2.autoconversion(SB2006.acnv, q_liq, q_rai, ρ_air, 1e8).dq_rai_dt for q_liq in q_liq_range] -q_liq_K1969 = [CM1.conv_q_liq_to_q_rai(rain.acnv1M, q_liq) for q_liq in q_liq_range] - -N_d_KK2000 = [CM2.conv_q_liq_to_q_rai(KK2000, 5e-4, ρ_air, N_d = N_d) for N_d in N_d_range] -N_d_B1994 = [CM2.conv_q_liq_to_q_rai(B1994, 5e-4, ρ_air, N_d = N_d) for N_d in N_d_range] -N_d_TC1980 = [CM2.conv_q_liq_to_q_rai(TC1980, 5e-4, ρ_air, N_d = N_d) for N_d in N_d_range] -N_d_LD2004 = [CM2.conv_q_liq_to_q_rai(LD2004, 5e-4, ρ_air, N_d = N_d) for N_d in N_d_range] -N_d_VarTimeScaleAcnv = [CM2.conv_q_liq_to_q_rai(VarTSc, 5e-4, ρ_air, N_d = N_d) for N_d in N_d_range] -N_d_SB2006 = [CM2.autoconversion(SB2006.acnv, q_liq, q_rai, ρ_air, N_d).dq_rai_dt for N_d in N_d_range] - -accKK2000_q_liq = [CM2.accretion(KK2000, q_liq, q_rai, ρ_air) for q_liq in q_liq_range] -accB1994_q_liq = [CM2.accretion(B1994, q_liq, q_rai, ρ_air) for q_liq in q_liq_range] -accTC1980_q_liq = [CM2.accretion(TC1980, q_liq, q_rai) for q_liq in q_liq_range] -accSB2006_q_liq = [CM2.accretion(SB2006, q_liq, q_rai, ρ_air, 1e8).dq_rai_dt for q_liq in q_liq_range] -accK1969_q_liq = [CM1.accretion(liquid, rain, blk1mvel.rain, ce, q_liq, q_rai, ρ_air) for q_liq in q_liq_range] - -accKK2000_q_rai = [CM2.accretion(KK2000, q_liq, q_rai, ρ_air) for q_rai in q_rai_range] -accB1994_q_rai = [CM2.accretion(B1994, q_liq, q_rai, ρ_air) for q_rai in q_rai_range] -accTC1980_q_rai = [CM2.accretion(TC1980, q_liq, q_rai) for q_rai in q_rai_range] -accSB2006_q_rai = [CM2.accretion(SB2006, q_liq, q_rai, ρ_air, 1e8).dq_rai_dt for q_rai in q_rai_range] -accK1969_q_rai = [CM1.accretion(liquid, rain, blk1mvel.rain, ce, q_liq, q_rai, ρ_air) for q_rai in q_rai_range] - -fig = Figure() - -ax1 = Axis(fig[1, 1]; yscale = log10) -ax2 = Axis(fig[1, 2]; xscale = log10, yscale = log10) -ax3 = Axis(fig[2, 1]; yscale = log10) -ax4 = Axis(fig[2, 2]; yscale = log10) - -ylims!(ax1, [1e-13, 1e-5]) -ylims!(ax2, [1e-13, 1e-5]) -ylims!(ax3, [1e-7, 5e-6]) -ylims!(ax4, [1e-7, 5e-6]) - -l1 = lines!(ax1, q_liq_range * 1e3, q_liq_KK2000, color = :red) -l2 = lines!(ax1, q_liq_range * 1e3, q_liq_B1994, color = :green) -l3 = lines!(ax1, q_liq_range * 1e3, q_liq_TC1980, color = :blue) -l4 = lines!(ax1, q_liq_range * 1e3, q_liq_LD2004, color = :purple) -l5 = lines!(ax1, q_liq_range * 1e3, q_liq_K1969, color = :black) -l6 = lines!(ax1, KK2000_x_q_liq, KK2000_y_q_liq, color = :red, linestyle = :dash) -l7 = lines!(ax1, B1994_x_q_liq, B1994_y_q_liq, color = :green, linestyle = :dash) -l8 = lines!(ax1, TC1980_x_q_liq, TC1980_y_q_liq, color = :blue, linestyle = :dash) -l9 = lines!(ax1, LD2004_x_q_liq, LD2004_y_q_liq, color = :purple, linestyle = :dash) - -l10 = lines!(ax2, N_d_range * 1e-6, N_d_KK2000, color = :red) -l11 = lines!(ax2, N_d_range * 1e-6, N_d_B1994, color = :green) -l12 = lines!(ax2, N_d_range * 1e-6, N_d_TC1980, color = :blue) -l13 = lines!(ax2, N_d_range * 1e-6, N_d_LD2004, color = :purple) -l14 = lines!(ax2, KK2000_x_N_d, KK2000_y_N_d, color = :red, linestyle = :dash) -l15 = lines!(ax2, B1994_x_N_d, B1994_y_N_d, color = :green, linestyle = :dash) -l16 = lines!(ax2, TC1980_x_N_d, TC1980_y_N_d, color = :blue, linestyle = :dash) -l17 = lines!(ax2, LD2004_x_N_d, LD2004_y_N_d, color = :purple, linestyle = :dash) - -l18 = lines!(ax3, q_liq_range * 1e3, accKK2000_q_liq, color = :red) -l19 = lines!(ax3, q_liq_range * 1e3, accB1994_q_liq, color = :green) -l20 = lines!(ax3, q_liq_range * 1e3, accTC1980_q_liq, color = :blue) -l21 = lines!(ax3, q_liq_range * 1e3, accK1969_q_liq, color = :black) - -l22 = lines!(ax4, q_rai_range * 1e3, accKK2000_q_rai, color = :red) -l23 = lines!(ax4, q_rai_range * 1e3, accB1994_q_rai, color = :green) -l24 = lines!(ax4, q_rai_range * 1e3, accTC1980_q_rai, color = :blue) -l25 = lines!(ax4, q_rai_range * 1e3, accK1969_q_rai, color = :black) - -l26 = lines!(ax1, q_liq_range * 1e3, q_liq_SB2006, color = :cyan) -l27 = lines!(ax2, N_d_range * 1e-6, N_d_SB2006, color = :cyan) -l28 = lines!(ax3, q_liq_range * 1e3, accSB2006_q_liq, color = :cyan) -l28 = lines!(ax4, q_rai_range * 1e3, accSB2006_q_rai, color = :cyan) - -l29 = lines!(ax1, q_liq_range * 1e3, q_liq_VarTimeScaleAcnv, color = :orange) -l30 = lines!(ax2, N_d_range * 1e-6, N_d_VarTimeScaleAcnv, color = :orange) - -ax1.xlabel = "q_liq [g/kg]" -ax1.ylabel = "autoconversion rate [1/s]" -ax2.xlabel = "N_d [1/cm3]" -ax2.ylabel = "autoconversion rate [1/s]" -ax3.xlabel = "q_liq [g/kg] (q_rai = 0.5 g/kg)" -ax3.ylabel = "accretion rate [1/s]" -ax4.xlabel = "q_rai [g/kg] (q_liq = 0.5 g/kg)" -ax4.ylabel = "accretion rate [1/s]" - -Legend( - fig[1, 3], - [l1, l2, l3, l4, l26, l5, l6, l7, l8, l9, l29], - ["KK2000", "B1994", "TC1980", "LD2004", "SB2006", "K1969", "Wood_KK2000", "Wood_B1994", "Wood_TC1980", "Wood_LD2004", "SA2023"] -) -save("Autoconversion_accretion.svg", fig) +```@example +include("plots/Microphysics2M_plots.jl") ``` ![](Autoconversion_accretion.svg) diff --git a/docs/src/ThresholdsTransition.md b/docs/src/ThresholdsTransition.md index a3f74d602..b7419a444 100644 --- a/docs/src/ThresholdsTransition.md +++ b/docs/src/ThresholdsTransition.md @@ -29,92 +29,8 @@ The curve defined by the above equation smoothly transition from ``f(x) = 0``, f ## Example figures -```@example example_figures -import Plots -import CloudMicrophysics -import CLIMAParameters - -const PL = Plots -const CM1 = CloudMicrophysics.Microphysics1M -const CM2 = CloudMicrophysics.Microphysics2M -const CP = CLIMAParameters -const CMP = CloudMicrophysics.Parameters - -FT = Float64 - -rain = [] -B1994 = [] -TC1980 = [] -LD2004 = [] - -k_thrshld_stpnss_values = [5.0, 2.0, 12.0] -for i in 1:3 - override_file = joinpath("override_dict.toml") - open(override_file, "w") do io - println(io, "[threshold_smooth_transition_steepness]") - println(io, "alias = \"k_thrshld_stpnss\"") - println(io, "value = " * string(k_thrshld_stpnss_values[i])) - println(io, "type = \"float\"") - end - toml_dict = CP.create_toml_dict(FT; override_file) - isfile(override_file) && rm(override_file; force=true) - - push!(rain, CMP.Rain(FT, toml_dict)) - push!(B1994, CMP.B1994(FT, toml_dict)) - push!(TC1980, CMP.TC1980(FT, toml_dict)) - push!(LD2004, CMP.LD2004(FT, toml_dict)) -end - -# example values -q_liq_range = range(1e-8, stop=1.5e-3, length=1000) -N_d_range = range(1e7, stop=1e9, length=1000) -ρ_air = 1.0 # kg m^-3 - -q_liq_K1969 = [CM1.conv_q_liq_to_q_rai(rain[1].acnv1M, q_liq) for q_liq in q_liq_range] -q_liq_K1969_s = [CM1.conv_q_liq_to_q_rai(rain[1].acnv1M, q_liq, smooth_transition = true) for q_liq in q_liq_range] - -q_liq_TC1980 = [CM2.conv_q_liq_to_q_rai(TC1980[2], q_liq, ρ_air) for q_liq in q_liq_range] -q_liq_TC1980_s = [CM2.conv_q_liq_to_q_rai(TC1980[2], q_liq, ρ_air, smooth_transition = true) for q_liq in q_liq_range] -q_liq_LD2004 = [CM2.conv_q_liq_to_q_rai(LD2004[2], q_liq, ρ_air) for q_liq in q_liq_range] -q_liq_LD2004_s = [CM2.conv_q_liq_to_q_rai(LD2004[2], q_liq, ρ_air, smooth_transition = true) for q_liq in q_liq_range] - -N_d_B1994 = [CM2.conv_q_liq_to_q_rai(B1994[3], 5e-4, ρ_air, N_d = N_d) for N_d in N_d_range] -N_d_B1994_s = [CM2.conv_q_liq_to_q_rai(B1994[3], 5e-4, ρ_air, N_d = N_d, smooth_transition = true) for N_d in N_d_range] - - -PL.plot(q_liq_range * 1e3, q_liq_K1969, - linewidth=2, - xlabel="q_liq [g/kg]", - ylabel="autoconversion rate [1/s]", - label="K1969 without smoothing" -) -PL.plot!(q_liq_range * 1e3, q_liq_K1969_s, linewidth=2, label="K1969 with smoothing") -PL.savefig("q_liq_K1969.svg") # hide - -PL.plot(q_liq_range * 1e3, q_liq_TC1980, - linewidth=2, - xlabel="q_liq [g/kg]", - ylabel="autoconversion rate [1/s]", - label="TC1980 without smoothing", - yaxis=:log, - ylim=(1e-10, 1e-5) -) -PL.plot!(q_liq_range * 1e3, q_liq_TC1980_s, linewidth=2, label="TC1980 with smoothing") -PL.plot!(q_liq_range * 1e3, q_liq_LD2004, linewidth=2, label="LD2004 without smoothing") -PL.plot!(q_liq_range * 1e3, q_liq_LD2004_s, linewidth=2, label="LD2004 with smoothing") -PL.savefig("q_liq_TC1980_LD2004.svg") # hide - -PL.plot(N_d_range * 1e-6, N_d_B1994, - linewidth=2, - xlabel="N_d [1/cm3]", - ylabel="autoconversion rate [1/s]", - label="B1994 without smoothing", - xaxis=:log, - yaxis=:log, - ylim=(1e-13, 1e-5) -) -PL.plot!(N_d_range * 1e-6, N_d_B1994_s, linewidth=2, label="B1994 with smoothing") -PL.savefig("N_d_B1994.svg") # hide +```@example +include("plots/Thersholds_transitions.jl") ``` ![](q_liq_K1969.svg) ![](q_liq_TC1980_LD2004.svg) diff --git a/docs/src/plots/ARGplots_fig1.jl b/docs/src/plots/ARGplots_fig1.jl new file mode 100644 index 000000000..15d2c2a36 --- /dev/null +++ b/docs/src/plots/ARGplots_fig1.jl @@ -0,0 +1,104 @@ +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") diff --git a/docs/src/plots/Microphysics1M_plots.jl b/docs/src/plots/Microphysics1M_plots.jl new file mode 100644 index 000000000..d0cab59ef --- /dev/null +++ b/docs/src/plots/Microphysics1M_plots.jl @@ -0,0 +1,564 @@ +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.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.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.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.0 +ϵ = 1.0 / TD.Parameters.molmass_ratio(tps) +p_sat = TD.saturation_vapor_pressure(tps, T, TD.Liquid()) +q_sat = ϵ * p_sat / (p + p_sat * (ϵ - 1.0)) +q_rain_range = range(1e-8, stop = 5e-3, length = 100) +q_tot = 15e-3 +q_vap = 0.15 * q_sat +q_ice = 0.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.0 +ϵ = 1.0 / TD.Parameters.molmass_ratio(tps) +p_sat = TD.saturation_vapor_pressure(tps, T, TD.Ice()) +q_sat = ϵ * p_sat / (p + p_sat * (ϵ - 1.0)) +q_snow_range = range(1e-8, stop = 5e-3, length = 100) +q_tot = 15e-3 +q_vap = 0.15 * q_sat +q_liq = 0.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.0 +ϵ = 1.0 / TD.Parameters.molmass_ratio(tps) +p_sat = TD.saturation_vapor_pressure(tps, T, TD.Ice()) +q_sat = ϵ * p_sat / (p + p_sat * (ϵ - 1.0)) +q_snow_range = range(1e-8, stop = 5e-3, length = 100) +q_tot = 15e-3 +q_vap = 0.15 * q_sat +q_liq = 0.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 diff --git a/docs/src/plots/Microphysics2M_plots.jl b/docs/src/plots/Microphysics2M_plots.jl new file mode 100644 index 000000000..d8a223628 --- /dev/null +++ b/docs/src/plots/Microphysics2M_plots.jl @@ -0,0 +1,206 @@ +using CairoMakie +CairoMakie.activate!(type = "svg") + +import CLIMAParameters +import Thermodynamics as TD +import CloudMicrophysics +import CloudMicrophysics.Microphysics1M as CM1 +import CloudMicrophysics.Microphysics2M as CM2 +import CloudMicrophysics.Parameters as CMP + +FT = Float64 + +const tps = TD.Parameters.ThermodynamicsParameters(FT) +const aps = CMP.AirProperties(FT) + +const KK2000 = CMP.KK2000(FT) +const B1994 = CMP.B1994(FT) +const TC1980 = CMP.TC1980(FT) +const LD2004 = CMP.LD2004(FT) +const VarTSc = CMP.VarTimescaleAcnv(FT) +const SB2006 = CMP.SB2006(FT) + +const ce = CMP.CollisionEff(FT) + +const liquid = CMP.CloudLiquid(FT) +const rain = CMP.Rain(FT) +const blk1mvel = CMP.Blk1MVelType(FT) + +include( + joinpath(pkgdir(CloudMicrophysics), "docs", "src", "plots", "Wooddata.jl"), +) + +# Example values +q_liq_range = range(1e-8, stop = 1e-3, length = 1000) +q_rai_range = range(1e-8, stop = 1e-3, length = 1000) +N_d_range = range(1e7, stop = 1e9, length = 1000) +q_liq = 5e-4 +q_rai = 5e-4 +ρ_air = 1.0 # kg m^-3 + +q_liq_KK2000 = [ + CM2.conv_q_liq_to_q_rai(KK2000, q_liq, ρ_air, N_d = 1e8) for + q_liq in q_liq_range +] +q_liq_B1994 = [ + CM2.conv_q_liq_to_q_rai(B1994, q_liq, ρ_air, N_d = 1e8) for + q_liq in q_liq_range +] +q_liq_TC1980 = [ + CM2.conv_q_liq_to_q_rai(TC1980, q_liq, ρ_air, N_d = 1e8) for + q_liq in q_liq_range +] +q_liq_LD2004 = [ + CM2.conv_q_liq_to_q_rai(LD2004, q_liq, ρ_air, N_d = 1e8) for + q_liq in q_liq_range +] +q_liq_VarTimeScaleAcnv = [ + CM2.conv_q_liq_to_q_rai(VarTSc, q_liq, ρ_air, N_d = 1e8) for + q_liq in q_liq_range +] +q_liq_SB2006 = [ + CM2.autoconversion(SB2006.acnv, q_liq, q_rai, ρ_air, 1e8).dq_rai_dt for + q_liq in q_liq_range +] +q_liq_K1969 = + [CM1.conv_q_liq_to_q_rai(rain.acnv1M, q_liq) for q_liq in q_liq_range] + +N_d_KK2000 = [ + CM2.conv_q_liq_to_q_rai(KK2000, 5e-4, ρ_air, N_d = N_d) for N_d in N_d_range +] +N_d_B1994 = [ + CM2.conv_q_liq_to_q_rai(B1994, 5e-4, ρ_air, N_d = N_d) for N_d in N_d_range +] +N_d_TC1980 = [ + CM2.conv_q_liq_to_q_rai(TC1980, 5e-4, ρ_air, N_d = N_d) for N_d in N_d_range +] +N_d_LD2004 = [ + CM2.conv_q_liq_to_q_rai(LD2004, 5e-4, ρ_air, N_d = N_d) for N_d in N_d_range +] +N_d_VarTimeScaleAcnv = [ + CM2.conv_q_liq_to_q_rai(VarTSc, 5e-4, ρ_air, N_d = N_d) for N_d in N_d_range +] +N_d_SB2006 = [ + CM2.autoconversion(SB2006.acnv, q_liq, q_rai, ρ_air, N_d).dq_rai_dt for + N_d in N_d_range +] + +accKK2000_q_liq = + [CM2.accretion(KK2000, q_liq, q_rai, ρ_air) for q_liq in q_liq_range] +accB1994_q_liq = + [CM2.accretion(B1994, q_liq, q_rai, ρ_air) for q_liq in q_liq_range] +accTC1980_q_liq = [CM2.accretion(TC1980, q_liq, q_rai) for q_liq in q_liq_range] +accSB2006_q_liq = [ + CM2.accretion(SB2006, q_liq, q_rai, ρ_air, 1e8).dq_rai_dt for + q_liq in q_liq_range +] +accK1969_q_liq = [ + CM1.accretion(liquid, rain, blk1mvel.rain, ce, q_liq, q_rai, ρ_air) for + q_liq in q_liq_range +] + +accKK2000_q_rai = + [CM2.accretion(KK2000, q_liq, q_rai, ρ_air) for q_rai in q_rai_range] +accB1994_q_rai = + [CM2.accretion(B1994, q_liq, q_rai, ρ_air) for q_rai in q_rai_range] +accTC1980_q_rai = [CM2.accretion(TC1980, q_liq, q_rai) for q_rai in q_rai_range] +accSB2006_q_rai = [ + CM2.accretion(SB2006, q_liq, q_rai, ρ_air, 1e8).dq_rai_dt for + q_rai in q_rai_range +] +accK1969_q_rai = [ + CM1.accretion(liquid, rain, blk1mvel.rain, ce, q_liq, q_rai, ρ_air) for + q_rai in q_rai_range +] + +fig = Figure(resolution = (900, 600)) + +ax1 = Axis(fig[1, 1]; yscale = log10) +ax2 = Axis(fig[1, 2]; xscale = log10, yscale = log10) +ax3 = Axis(fig[2, 1]; yscale = log10) +ax4 = Axis(fig[2, 2]; yscale = log10) + +ylims!(ax1, [1e-13, 1e-5]) +ylims!(ax2, [1e-13, 1e-5]) +ylims!(ax3, [1e-7, 5e-6]) +ylims!(ax4, [1e-7, 5e-6]) + +l1 = lines!(ax1, q_liq_range * 1e3, q_liq_KK2000, color = :red) +l2 = lines!(ax1, q_liq_range * 1e3, q_liq_B1994, color = :green) +l3 = lines!(ax1, q_liq_range * 1e3, q_liq_TC1980, color = :blue) +l4 = lines!(ax1, q_liq_range * 1e3, q_liq_LD2004, color = :purple) +l5 = lines!(ax1, q_liq_range * 1e3, q_liq_K1969, color = :black) +l6 = + lines!(ax1, KK2000_x_q_liq, KK2000_y_q_liq, color = :red, linestyle = :dash) +l7 = + lines!(ax1, B1994_x_q_liq, B1994_y_q_liq, color = :green, linestyle = :dash) +l8 = lines!( + ax1, + TC1980_x_q_liq, + TC1980_y_q_liq, + color = :blue, + linestyle = :dash, +) +l9 = lines!( + ax1, + LD2004_x_q_liq, + LD2004_y_q_liq, + color = :purple, + linestyle = :dash, +) + +l10 = lines!(ax2, N_d_range * 1e-6, N_d_KK2000, color = :red) +l11 = lines!(ax2, N_d_range * 1e-6, N_d_B1994, color = :green) +l12 = lines!(ax2, N_d_range * 1e-6, N_d_TC1980, color = :blue) +l13 = lines!(ax2, N_d_range * 1e-6, N_d_LD2004, color = :purple) +l14 = lines!(ax2, KK2000_x_N_d, KK2000_y_N_d, color = :red, linestyle = :dash) +l15 = lines!(ax2, B1994_x_N_d, B1994_y_N_d, color = :green, linestyle = :dash) +l16 = lines!(ax2, TC1980_x_N_d, TC1980_y_N_d, color = :blue, linestyle = :dash) +l17 = + lines!(ax2, LD2004_x_N_d, LD2004_y_N_d, color = :purple, linestyle = :dash) + +l18 = lines!(ax3, q_liq_range * 1e3, accKK2000_q_liq, color = :red) +l19 = lines!(ax3, q_liq_range * 1e3, accB1994_q_liq, color = :green) +l20 = lines!(ax3, q_liq_range * 1e3, accTC1980_q_liq, color = :blue) +l21 = lines!(ax3, q_liq_range * 1e3, accK1969_q_liq, color = :black) + +l22 = lines!(ax4, q_rai_range * 1e3, accKK2000_q_rai, color = :red) +l23 = lines!(ax4, q_rai_range * 1e3, accB1994_q_rai, color = :green) +l24 = lines!(ax4, q_rai_range * 1e3, accTC1980_q_rai, color = :blue) +l25 = lines!(ax4, q_rai_range * 1e3, accK1969_q_rai, color = :black) + +l26 = lines!(ax1, q_liq_range * 1e3, q_liq_SB2006, color = :cyan) +l27 = lines!(ax2, N_d_range * 1e-6, N_d_SB2006, color = :cyan) +l28 = lines!(ax3, q_liq_range * 1e3, accSB2006_q_liq, color = :cyan) +l28 = lines!(ax4, q_rai_range * 1e3, accSB2006_q_rai, color = :cyan) + +l29 = lines!(ax1, q_liq_range * 1e3, q_liq_VarTimeScaleAcnv, color = :orange) +l30 = lines!(ax2, N_d_range * 1e-6, N_d_VarTimeScaleAcnv, color = :orange) + +ax1.xlabel = "q_liq [g/kg]" +ax1.ylabel = "autoconversion rate [1/s]" +ax2.xlabel = "N_d [1/cm3]" +ax2.ylabel = "autoconversion rate [1/s]" +ax3.xlabel = "q_liq [g/kg] (q_rai = 0.5 g/kg)" +ax3.ylabel = "accretion rate [1/s]" +ax4.xlabel = "q_rai [g/kg] (q_liq = 0.5 g/kg)" +ax4.ylabel = "accretion rate [1/s]" + +Legend( + fig[1, 3], + [l1, l2, l3, l4, l26, l5, l6, l7, l8, l9, l29], + [ + "KK2000", + "B1994", + "TC1980", + "LD2004", + "SB2006", + "K1969", + "Wood_KK2000", + "Wood_B1994", + "Wood_TC1980", + "Wood_LD2004", + "SA2023", + ], +) +save("Autoconversion_accretion.svg", fig) diff --git a/docs/src/plots/Thersholds_transitions.jl b/docs/src/plots/Thersholds_transitions.jl new file mode 100644 index 000000000..99a4331c5 --- /dev/null +++ b/docs/src/plots/Thersholds_transitions.jl @@ -0,0 +1,135 @@ +import Plots +import CloudMicrophysics +import CLIMAParameters + +const PL = Plots +const CM1 = CloudMicrophysics.Microphysics1M +const CM2 = CloudMicrophysics.Microphysics2M +const CP = CLIMAParameters +const CMP = CloudMicrophysics.Parameters + +FT = Float64 + +rain = [] +B1994 = [] +TC1980 = [] +LD2004 = [] + +k_thrshld_stpnss_values = [5.0, 2.0, 12.0] +for i in 1:3 + override_file = joinpath("override_dict.toml") + open(override_file, "w") do io + println(io, "[threshold_smooth_transition_steepness]") + println(io, "alias = \"k_thrshld_stpnss\"") + println(io, "value = " * string(k_thrshld_stpnss_values[i])) + println(io, "type = \"float\"") + end + toml_dict = CP.create_toml_dict(FT; override_file) + isfile(override_file) && rm(override_file; force = true) + + push!(rain, CMP.Rain(FT, toml_dict)) + push!(B1994, CMP.B1994(FT, toml_dict)) + push!(TC1980, CMP.TC1980(FT, toml_dict)) + push!(LD2004, CMP.LD2004(FT, toml_dict)) +end + +# example values +q_liq_range = range(1e-8, stop = 1.5e-3, length = 1000) +N_d_range = range(1e7, stop = 1e9, length = 1000) +ρ_air = 1.0 # kg m^-3 + +q_liq_K1969 = + [CM1.conv_q_liq_to_q_rai(rain[1].acnv1M, q_liq) for q_liq in q_liq_range] +q_liq_K1969_s = [ + CM1.conv_q_liq_to_q_rai(rain[1].acnv1M, q_liq, smooth_transition = true) for q_liq in q_liq_range +] + +q_liq_TC1980 = + [CM2.conv_q_liq_to_q_rai(TC1980[2], q_liq, ρ_air) for q_liq in q_liq_range] +q_liq_TC1980_s = [ + CM2.conv_q_liq_to_q_rai(TC1980[2], q_liq, ρ_air, smooth_transition = true) for q_liq in q_liq_range +] +q_liq_LD2004 = + [CM2.conv_q_liq_to_q_rai(LD2004[2], q_liq, ρ_air) for q_liq in q_liq_range] +q_liq_LD2004_s = [ + CM2.conv_q_liq_to_q_rai(LD2004[2], q_liq, ρ_air, smooth_transition = true) for q_liq in q_liq_range +] + +N_d_B1994 = [ + CM2.conv_q_liq_to_q_rai(B1994[3], 5e-4, ρ_air, N_d = N_d) for + N_d in N_d_range +] +N_d_B1994_s = [ + CM2.conv_q_liq_to_q_rai( + B1994[3], + 5e-4, + ρ_air, + N_d = N_d, + smooth_transition = true, + ) for N_d in N_d_range +] + +PL.plot( + q_liq_range * 1e3, + q_liq_K1969, + linewidth = 2, + xlabel = "q_liq [g/kg]", + ylabel = "autoconversion rate [1/s]", + label = "K1969 without smoothing", +) +PL.plot!( + q_liq_range * 1e3, + q_liq_K1969_s, + linewidth = 2, + label = "K1969 with smoothing", +) +PL.savefig("q_liq_K1969.svg") # hide + +PL.plot( + q_liq_range * 1e3, + q_liq_TC1980, + linewidth = 2, + xlabel = "q_liq [g/kg]", + ylabel = "autoconversion rate [1/s]", + label = "TC1980 without smoothing", + yaxis = :log, + ylim = (1e-10, 1e-5), +) +PL.plot!( + q_liq_range * 1e3, + q_liq_TC1980_s, + linewidth = 2, + label = "TC1980 with smoothing", +) +PL.plot!( + q_liq_range * 1e3, + q_liq_LD2004, + linewidth = 2, + label = "LD2004 without smoothing", +) +PL.plot!( + q_liq_range * 1e3, + q_liq_LD2004_s, + linewidth = 2, + label = "LD2004 with smoothing", +) +PL.savefig("q_liq_TC1980_LD2004.svg") # hide + +PL.plot( + N_d_range * 1e-6, + N_d_B1994, + linewidth = 2, + xlabel = "N_d [1/cm3]", + ylabel = "autoconversion rate [1/s]", + label = "B1994 without smoothing", + xaxis = :log, + yaxis = :log, + ylim = (1e-13, 1e-5), +) +PL.plot!( + N_d_range * 1e-6, + N_d_B1994_s, + linewidth = 2, + label = "B1994 with smoothing", +) +PL.savefig("N_d_B1994.svg") # hide