g=sigma fix. Plotting cleanup

docs/src/IceNucleationParcel0D.md

@@ -296,7 +296,7 @@ include("../../parcel/Example_Liquid_only.jl")
 The plots below are the results of the adiabatic parcel model
   with immersion freezing, condensation growth, and deposition growth for
   both a monodisperse and gamma size distribution. Note that this has not
-  yet been validated against literature. Results are very sensitive to 
+  yet been validated against literature. Results are very sensitive to
   initial conditions.
 
 ```@example
@@ -339,15 +339,23 @@ include("../../parcel/Example_P3_ice_nuc.jl")
 ## Immersion Freezing Parametrization based on Frostenberg et al. 2023
 
 The parcel model also includes the parametrization of immersion freezing based on [Frostenberg2023](@cite).
-The concentration of ice nucleating particles (INPs) is randomly selected based on a lognormal relative frequency distribution, depending only on the temperature. If the random INP concentration exceeds the existing concentration of ice crystals, new ice crystals are created, such that their total concentration is equal to the new INP concentration.
+The concentration of ice nucleating particles (INPs) is randomly selected
+  based on a lognormal relative frequency distribution,
+  and depends only on the temperature.
+If the INP concentration exceeds the existing concentration of ice crystals,
+  new ice crystals are created, such that their total concentration is equal to the new INP concentration.
 
 Three different implementations of this parametrization are used in the parcel model:
 - **Frostenberg_mean**, in which the concentration of INPs is equal to its mean value defined in [Frostenberg2023](@cite) for a given temperature,
-- **Frostenberg_random**, in which the concentration of INPs is drawn randomly from the distribution defined in [Frostenberg2023](@cite). The time between draws is set by the *drawing_interval* parameter,
-- **Frostenberg_stochastic**, in which the freezing is treated as a stochastic process, with mean defined in [Frostenberg2023](@cite) and Wiener noise. The timescale of the process is set by the `\gamma` parameter. 
+- **Frostenberg_random**, in which the concentration of INPs is drawn randomly
+  from the distribution defined in [Frostenberg2023](@cite).
+  The time between draws is set by the *drawing_interval* parameter,
+- **Frostenberg_stochastic**, in which the freezing is treated as a stochastic process,
+  with mean defined in [Frostenberg2023](@cite) and Wiener noise.
+  The timescale of the process is set by the `\gamma` parameter.
 
-
-The stochastic implementation is based on the equation for a stochastic process $x$ with Gaussian random noise:
+The stochastic implementation is based on the equation for a stochastic process $x$
+  with Gaussian random noise:
 
 ```math
 \begin{equation}
@@ -399,16 +407,13 @@ We map this notation onto our problem:
   x(t) = {\bf{N}}(x_0 e^{-\gamma t} + (1-e^{-\gamma t}) \mu(T), \frac{g^2}{2\gamma} (1-e^{-2\gamma t}))
 \end{equation}
 ```
-If ``\gamma >> 1``, then
-```math
-\begin{equation}
-  x(t) = {\bf{N}}(\mu(T), \frac{g^2}{2\gamma})
-\end{equation}
-```
-Hence ``g = \sigma * \sqrt{2 \gamma}``.
-
+In the limit of ``\gamma >> 1``, we want ``x(t) = \mu(T)``.
+Hence ``g = \sigma``.
 
-The following plot shows resuls of the parcel model with Frostenberg_mean parametrization, Frostenberg_random parametrization with different drawing intervals `t_d` and Frostenberg_stochastic parametrization with different time scales `\gamma`. The timestep of the model is 1 s.
+The following plot shows resuls of the parcel model with the `mean` (black line)
+ `random` parametrization with different drawing intervals `t_d` (dotted lines)
+ `stochastic` parametrization with different time scales `\gamma` (solid lines).
+The timestep of the model is 1 s.
 
 ```@example
 include("../../parcel/Example_Frostenberg_Immersion_Freezing.jl")

parcel/Example_Frostenberg_Immersion_Freezing.jl

@@ -3,8 +3,10 @@ import CairoMakie as MK
 import Thermodynamics as TD
 import CloudMicrophysics as CM
 import ClimaParams as CP
+
 import Random as RAND
-random_seeds = [0, 1234, 5678, 3443]
+RAND.seed!(44)
+N_ensemble = 32
 
 # definition of the ODE problem for parcel model
 include(joinpath(pkgdir(CM), "parcel", "Parcel.jl"))
@@ -22,7 +24,7 @@ r₀ = FT(1e-6)
 p₀ = FT(800 * 1e2)
 T₀ = FT(251)
 qᵥ = FT(8.1e-4)
-qₗ = Nₗ * 4 / 3 * FT(π) * r₀^3 * ρₗ / FT(1.2) # 1.2 should be ρₐ
+qₗ = Nₗ * 4 / 3 * FT(π) * r₀^3 * ρₗ / FT(1.2) # ρₐ ~ 1.2
 qᵢ = FT(0)
 x_sulph = FT(0.01)
 
@@ -39,13 +41,14 @@ IC = [Sₗ, p₀, T₀, qᵥ, qₗ, qᵢ, Nₐ, Nₗ, Nᵢ, x_sulph]
 w = FT(0.7)                                # updraft speed
 const_dt = FT(1)                           # model timestep
 t_max = FT(1200)
-aerosol = CMP.Illite(FT)
 condensation_growth = "Condensation"
 deposition_growth = "Deposition"
 DSD = "Monodisperse"
+drawing_interval_range = [1, 10, 500]
+γ_range = [1, 10, 500] .* const_dt
 
 # Plotting
-fig = MK.Figure(resolution = (900, 700))
+fig = MK.Figure(size = (900, 700))
 ax1 = MK.Axis(fig[1, 1], ylabel = "Ice Supersaturation [-]")
 ax2 = MK.Axis(fig[1, 2], ylabel = "Temperature [K]")
 ax3 = MK.Axis(fig[2, 1], ylabel = "q_ice [g/kg]")
@@ -121,7 +124,6 @@ end
 params = parcel_params{FT}(
     const_dt = const_dt,
     w = w,
-    aerosol = aerosol,
     heterogeneous = "Frostenberg_mean",
     condensation_growth = condensation_growth,
     deposition_growth = deposition_growth,
@@ -131,64 +133,48 @@ params = parcel_params{FT}(
 sol = run_parcel(IC, FT(0), t_max, params)
 plot_results(sol.t, sol, "mean")
 
-
 # Frostenberg_random with different drawing frequencies
-drawing_interval_range = range(FT(1), stop = FT(5), length = 3)
-
 for (drawing_interval, color) in zip(drawing_interval_range, colors)
-
     # creating an ensamble of solutions
     solutions = []
-    for random_seed in random_seeds
-
-        RAND.seed!(trunc(Int, random_seed)) #set the random seed
-
-        params = parcel_params{FT}(
+    for ensemble_member in range(1, length=N_ensemble)
+        local params = parcel_params{FT}(
             const_dt = const_dt,
             w = w,
-            aerosol = aerosol,
             heterogeneous = "Frostenberg_random",
             condensation_growth = condensation_growth,
             deposition_growth = deposition_growth,
             size_distribution = DSD,
             drawing_interval = drawing_interval,
         )
         # solve ODE
-        sol = run_parcel(IC, FT(0), t_max, params)
+        local sol = run_parcel(IC, FT(0), t_max, params)
         push!(solutions, sol)
     end
     mean_sol = sum(solutions) / length(solutions)
-    plot_results(sol.t, mean_sol, drawing_interval, "t_d=", " s", :dot, color)
+    plot_results(sol.t, mean_sol, drawing_interval, "f=", " dt", :dot, color)
 end
 
-
 # Frostenberg_stochastic with different timescales γ
-γ_range = range(FT(1), stop = FT(5), length = 3)
-
 for (γ, color) in zip(γ_range, colors)
-
     # creating an ensamble of solutions
     solutions = []
-    for random_seed in random_seeds
-
-        RAND.seed!(trunc(Int, random_seed)) #set the random seed
-
-        params = parcel_params{FT}(
+    for ensemble_member in range(1, length=N_ensemble)
+        local params = parcel_params{FT}(
             const_dt = const_dt,
             w = w,
-            aerosol = aerosol,
             heterogeneous = "Frostenberg_stochastic",
             condensation_growth = condensation_growth,
             deposition_growth = deposition_growth,
             size_distribution = DSD,
             γ = γ,
         )
         # solve ODE
-        sol = run_parcel(IC, FT(0), t_max, params)
+        local sol = run_parcel(IC, FT(0), t_max, params)
         push!(solutions, sol)
     end
     mean_sol = sum(solutions) / length(solutions)
-    plot_results(sol.t, mean_sol, γ, "γ=", " s", :solid, color)
+    plot_results(sol.t, mean_sol, γ, "γ=", " 1/s", :solid, color)
 end
 
 MK.save("frostenberg_immersion_freezing.svg", fig)

parcel/ParcelTendencies.jl

@@ -133,7 +133,7 @@ function immersion_freezing(params::Frostenberg_stochastic, PSD, state)
     (; ip, γ) = params
     (; T, Nₗ, Nᵢ, t) = state
     μ = CMI_het.INP_concentration_mean(T)
-    g = ip.σ * sqrt(2 * γ)
+    g = ip.σ
     mean = (1 - exp(-γ * t)) * μ
     st_dev = sqrt(g^2 / (2 * γ) * (1 - exp(-2 * γ * t)))
     INPC = exp(rand(DS.Normal(mean, st_dev)))