-
Notifications
You must be signed in to change notification settings - Fork 8
Commit
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
testing stuff out more cleanup add Base back Bump version for a new release No kwargs in parameters constructors Update ci to borsless life Goodbye bors and thank you fix docs builds Add a box model to test immersion freezing
- Loading branch information
Showing
2 changed files
with
179 additions
and
8 deletions.
There are no files selected for viewing
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,135 @@ | ||
| import OrdinaryDiffEq as ODE | ||
| import CairoMakie as MK | ||
| import Thermodynamics as TD | ||
| import CloudMicrophysics as CM | ||
| import CLIMAParameters as CP | ||
|
|
||
| # definition of the ODE problem for parcel model | ||
| include(joinpath(pkgdir(CM), "parcel", "parcel.jl")) | ||
|
|
||
| FT = Float64 | ||
|
|
||
| # Get free parameters | ||
| tps = CMP.ThermodynamicsParameters(FT) | ||
| wps = CMP.WaterProperties(FT) | ||
| aps = CMP.AirProperties(FT) | ||
| ip = CMP.IceNucleationParameters(FT) | ||
|
|
||
| # Constants | ||
| ρₗ = wps.ρw | ||
| R_v = TD.Parameters.R_v(tps) | ||
| R_d = TD.Parameters.R_d(tps) | ||
|
|
||
| # Initial conditions | ||
| Nₐ = FT(0) | ||
| Nₗ = FT(300 * 1e6) # Jensen2022 starts with Nₐ = 300 * 1e6 and activates them within parcel | ||
| Nᵢ = FT(0) | ||
| r₀ = FT(25 * 1e-9) # Value of dry aerosol r from Jensen2022 | ||
| p₀ = FT(800 * 1e2) | ||
| T₀ = FT(190) | ||
| x_sulph = FT(0) | ||
|
|
||
|
|
||
| eₛ = TD.saturation_vapor_pressure(tps, T₀, TD.Liquid()) | ||
| ξ(T) = | ||
| TD.saturation_vapor_pressure(tps, T, TD.Liquid()) / | ||
| TD.saturation_vapor_pressure(tps, T, TD.Ice()) | ||
| S_l(T, S_i) = @. (S_i + 1) / ξ(T) | ||
| Sᵢ = FT(1.55) | ||
| Sₗ = S_l(T₀, Sᵢ) | ||
| e = eₛ * Sₗ | ||
| ϵ = R_d / R_v | ||
|
|
||
| # mass per volume for dry air, vapor and liquid | ||
| md_v = (p₀ - e) / R_d / T₀ | ||
| mv_v = e / R_v / T₀ | ||
| ml_v = Nₗ * 4 / 3 * π * ρₗ * r₀^3 | ||
|
|
||
| qᵥ = mv_v / (md_v + mv_v + ml_v) | ||
| qₗ = ml_v / (md_v + mv_v + ml_v) | ||
| qᵢ = FT(0) | ||
|
|
||
| # Moisture dependent initial conditions | ||
| q = TD.PhasePartition(qᵥ + qₗ + qᵢ, qₗ, qᵢ) | ||
| ts = TD.PhaseNonEquil_pTq(tps, p₀, T₀, q) | ||
| ρₐ = TD.air_density(tps, ts) | ||
| Rₐ = TD.gas_constant_air(tps, q) | ||
|
|
||
| #e = qᵥ * p₀ * R_v / Rₐ | ||
|
|
||
| println("S_l initial = ", Sₗ ) | ||
| IC = [Sₗ, p₀, T₀, qᵥ, qₗ, qᵢ, Nₐ, Nₗ, Nᵢ, x_sulph] | ||
|
|
||
| # Simulation parameters passed into ODE solver | ||
| r_nuc = FT(25 * 1e-9) # assumed size of nucleated particles | ||
| w = FT(1) # updraft speed | ||
| α_m = FT(0.5) # accomodation coefficient | ||
| const_dt = FT(0.5) # model timestep | ||
| #t_max = FT(120) | ||
| t_max = FT(30) | ||
| aerosol = [] | ||
| ice_nucleation_modes = ["HomogeneousFreezing"] # homogeneous freezing only | ||
| growth_modes = [] # no growth | ||
| droplet_size_distribution_list = [["Monodisperse"]] | ||
|
|
||
| # Data from Jensen(2022) Figure 1 | ||
| # https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2022JD036535 | ||
| #! format: off | ||
| Jensen_t_sat = [0, 62.71, 70.52, 76.87, 82.4, 84.84, 88.1, 92, 96.07, 100.63, 105.35, 112.51, 119.83] | ||
| Jensen_sat = [1.55, 1.694, 1.7107, 1.7208, 1.725, 1.726, 1.7259, 1.722, 1.715, 1.702, 1.686, 1.653, 1.6126] | ||
| Jensen_t_T = [0, 120] | ||
| Jensen_T = [190, 189] | ||
| Jensen_t_ICNC = [0.217, 42.69, 50.02, 54.41, 58.97, 65.316, 72.477, 82.08, 92.658, 94.123, 95.5877, 119.84] | ||
| Jensen_ICNC = [0, 0, 0.282, 0.789, 1.804, 4.1165, 7.218, 12.12, 16.35, 16.8, 16.97, 17.086] | ||
| #! format: on | ||
|
|
||
| fig = MK.Figure(resolution = (800, 600)) | ||
| ax1 = MK.Axis(fig[1, 1], ylabel = "Ice Saturation") | ||
| ax2 = MK.Axis(fig[3, 1], xlabel = "Time [s]", ylabel = "Temperature [K]") | ||
| ax3 = MK.Axis(fig[2, 1], ylabel = "q_vap [g/kg]") | ||
| ax4 = MK.Axis(fig[2, 2], xlabel = "Time [s]", ylabel = "q_liq [g/kg]") | ||
| ax5 = MK.Axis(fig[1, 2], ylabel = "ICNC [cm^-3]") | ||
|
|
||
| MK.lines!(ax1, Jensen_t_sat, Jensen_sat, label = "Jensen 2022", color = :green) | ||
| MK.lines!(ax2, Jensen_t_T, Jensen_T, color = :green) | ||
| MK.lines!(ax5, Jensen_t_ICNC, Jensen_ICNC, color = :green) | ||
|
|
||
| for droplet_size_distribution in droplet_size_distribution_list | ||
| p = (; | ||
| wps, | ||
| aps, | ||
| tps, | ||
| ip, | ||
| const_dt, | ||
| r_nuc, | ||
| w, | ||
| α_m, | ||
| aerosol, | ||
| ice_nucleation_modes, | ||
| growth_modes, | ||
| droplet_size_distribution, | ||
| ) | ||
| # solve ODE | ||
| sol = run_parcel(IC, FT(0), t_max, p) | ||
|
|
||
| DSD = droplet_size_distribution[1] | ||
|
|
||
| # Plot results | ||
| MK.lines!(ax1, sol.t, S_i(sol[3, :], (sol[1, :])), label = DSD) | ||
| MK.lines!(ax2, sol.t, sol[3, :]) | ||
| MK.lines!(ax3, sol.t, sol[4, :] * 1e3) | ||
| MK.lines!(ax4, sol.t, sol[5, :] * 1e3) | ||
| MK.lines!(ax5, sol.t, sol[9, :] * 1e6) | ||
| end | ||
|
|
||
| #fig[3, 2] = MK.Legend(fig, ax1, framevisible = false) | ||
| MK.axislegend( | ||
| ax1, | ||
| framevisible = false, | ||
| labelsize = 12, | ||
| orientation = :horizontal, | ||
| nbanks = 2, | ||
| position = :rb, | ||
| ) | ||
|
|
||
| MK.save("Jensen_et_al_2022.svg", fig) |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters