adding suggestions

docs/src/IceNucleationParcel0D.md

@@ -197,7 +197,7 @@ where:
 - ``N_{act}`` is the number of activated ice particles.
 
 ``N_{act}`` can be computed for example from
-  [activated fraction](https://clima.github.io/CloudMicrophysics.jl/previews/PR103/IceNucleation/#Activated-fraction-for-deposition-freezing-on-dust) ``f_i``
+  [activated fraction](https://clima.github.io/CloudMicrophysics.jl/dev/IceNucleation/#Activated-fraction-for-deposition-freezing-on-dust) ``f_i``
 ```math
 \begin{equation}
   N_{act} = N_{aer} f_i
@@ -212,7 +212,7 @@ Following the water activity based immersion freezing model (ABIFM), the ABIFM d
   per second via immersion freezing can then be calculating using
 ```math
 \begin{equation}
-  P_{ice} = [\frac{dN_i}{dt}]_{immer} = J_{immer}A(N_{liq})
+  P_{ice} = \left[ \frac{dN_i}{dt} \right]_{immer} = J_{immer}A(N_{liq})
   \label{eq:ABIFM_P_ice}
 \end{equation}
 ```
@@ -221,12 +221,12 @@ where ``N_{liq}`` is total number of ice nuclei containing droplets and
 
 ## Homogeneous Freezing
 Homogeneous freezing follows the water-activity based model described in the
-  ``Ice Nucleation`` section which gives a nucleation rate coefficient of
+  [Ice Nucleation](https://clima.github.io/CloudMicrophysics.jl/dev/IceNucleation/) section which gives a nucleation rate coefficient of
   units ``cm^{-3}s^{-1}``.
 The ice production rate from homogeneous freezing can then be determined:
 ```math
 \begin{equation}
-  P_{ice} = [\frac{dN_i}{dt}]_{hom} = J_{hom}V(N_{liq})
+  P_{ice} = \left[ \frac{dN_i}{dt} \right]_{hom} = J_{hom}V(N_{liq})
   \label{eq:hom_P_ice}
 \end{equation}
 ```

parcel/Immersion_Freezing.jl

@@ -26,7 +26,7 @@ r₀ = FT(1e-6)
 p₀ = FT(800 * 1e2)
 T₀ = FT(251)
 qᵥ = FT(8.1e-4)
-qₗ = Nₗ * 4 / 3 * π * r₀^3 * ρₗ / 1.2 # 1.2 should be ρₐ
+qₗ = Nₗ * 4 / 3 * FT(π) * r₀^3 * ρₗ / 1.2 # 1.2 should be ρₐ
 qᵢ = FT(0)
 x_sulph = FT(0.01)
 
@@ -102,10 +102,10 @@ for droplet_size_distribution in droplet_size_distribution_list
     sol_Nᵢ = sol[9, :]
     sol_qₗ = sol[5, :]
     if DSD == "Monodisperse"
-        rₗ = cbrt.(sol_qₗ ./ sol_Nₗ ./ (4 / 3 * π) / ρₗ * ρₐ)
+        rₗ = cbrt.(sol_qₗ ./ sol_Nₗ ./ (4 / 3 * FT(π)) / ρₗ * ρₐ)
     end
     if DSD == "Gamma"
-        λ = cbrt.(32 .* π .* sol_Nₗ ./ sol_qₗ * ρₗ / ρₐ)
+        λ = cbrt.(32 .* FT(π) .* sol_Nₗ ./ sol_qₗ * ρₗ / ρₐ)
         rₗ = 2 ./ λ
     end
     MK.lines!(ax5, sol.t / 60, sol_Nₗ)

parcel/Jensen_et_al_2022.jl

@@ -27,7 +27,7 @@ R_d = TD.Parameters.R_d(tps)
 ξ(T) =
     TD.saturation_vapor_pressure(tps, T, TD.Liquid()) /
     TD.saturation_vapor_pressure(tps, T, TD.Ice())
-s_i(T, s_liq) = @. s_liq * ξ(T) # saturation ratio
+S_i(T, S_liq) = @. S_liq * ξ(T) # saturation ratio
 
 # Initial conditions
 Nₐ = FT(300 * 1e6)
@@ -46,12 +46,12 @@ qᵢ = FT(0)
 q = TD.PhasePartition(qᵥ + qₗ + qᵢ, qₗ, qᵢ)
 R_a = TD.gas_constant_air(tps, q)
 e = qᵥ * p₀ * R_v / R_a
-sₗ = e / eₛ
+Sₗ = e / eₛ
 
 ## Moisture dependent initial conditions
 ts = TD.PhaseNonEquil_pTq(tps, p₀, T₀, q)
 
-IC = [sₗ, p₀, T₀, qᵥ, qₗ, qᵢ, Nₐ, Nₗ, Nᵢ, x_sulph]
+IC = [Sₗ, p₀, T₀, qᵥ, qₗ, qᵢ, Nₐ, Nₗ, Nᵢ, x_sulph]
 
 # Simulation parameters passed into ODE solver
 r_nuc = r₀          # assumed size of nucleated particles
@@ -118,7 +118,7 @@ for droplet_size_distribution in droplet_size_distribution_list
     DSD = droplet_size_distribution[1]
 
     # Plot results
-    MK.lines!(ax1, sol.t, s_i(sol[3, :], (sol[1, :])), label = "ice", color = :blue)
+    MK.lines!(ax1, sol.t, S_i(sol[3, :], (sol[1, :])), label = "ice", color = :blue)
     MK.lines!(ax1, sol.t, (sol[1, :]), label = "liquid", color = :red) # liq saturation
     MK.lines!(ax2, sol.t, sol[3, :], label = "CM.jl")                  # temperature
     MK.lines!(ax3, sol.t, sol[4, :] * 1e3) # q_vap

parcel/Liquid_only.jl

@@ -36,7 +36,7 @@ e = eₛ
 # 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
+ml_v = Nₗ * 4 / 3 * FT(π) * ρₗ * r₀^3
 
 qᵥ = mv_v / (md_v + mv_v + ml_v)
 qₗ = ml_v / (md_v + mv_v + ml_v)

parcel/Project.toml

@@ -10,3 +10,6 @@ OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
 RootSolvers = "7181ea78-2dcb-4de3-ab41-2b8ab5a31e74"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 Thermodynamics = "b60c26fb-14c3-4610-9d3e-2d17fe7ff00c"
+
+[compat]
+CLIMAParameters = "0.7.12"

parcel/parcel.jl

@@ -5,7 +5,6 @@ import CloudMicrophysics.Common as CMO
 import CloudMicrophysics.HetIceNucleation as CMI_het
 import CloudMicrophysics.HomIceNucleation as CMI_hom
 import CloudMicrophysics.Parameters as CMP
-import Random as RD
 
 #! format: off
 """
@@ -28,7 +27,7 @@ function parcel_model(dY, Y, p, t)
     # TODO - We are solving for both q_ and supersaturation.
     # We should decide if we want to solve for S (as in the cirrus box model)
     # or for q_ which is closer to what ClimaAtmos is doing.
-    s_liq = Y[1]      # saturation ratio over liquid water
+    S_liq = Y[1]      # saturation ratio over liquid water
     p_a = Y[2]        # pressure
     T = Y[3]          # temperature
     q_vap = Y[4]      # vapor specific humidity
@@ -54,7 +53,7 @@ function parcel_model(dY, Y, p, t)
     e_si = TD.saturation_vapor_pressure(tps, T, TD.Ice())
     e_sl = TD.saturation_vapor_pressure(tps, T, TD.Liquid())
     ξ = e_sl / e_si
-    s_i = ξ * s_liq
+    S_i = ξ * S_liq
 
     # Constants and variables that depend on the moisture content
     R_a = TD.gas_constant_air(tps, q)
@@ -80,11 +79,11 @@ function parcel_model(dY, Y, p, t)
         AF = CMI_het.dust_activated_number_fraction(
             aerosol,
             ip.deposition,
-            s_i,
+            S_i,
             T,
         )
         dN_act_dt_depo = max(FT(0), AF * N_aer - N_ice) / const_dt
-        dqi_dt_new_depo = dN_act_dt_depo * 4 / 3 * π * r_nuc^3 * ρ_ice / ρ_air
+        dqi_dt_new_depo = dN_act_dt_depo * 4 / 3 * FT(π) * r_nuc^3 * ρ_ice / ρ_air
     end
 
     dN_act_dt_immersion = FT(0)
@@ -95,15 +94,15 @@ function parcel_model(dY, Y, p, t)
             CMO.a_w_eT(tps, e, T) - CMO.a_w_ice(tps, T)
         J_immersion = CMI_het.ABIFM_J(aerosol, Δa_w)
         if "Monodisperse" in droplet_size_distribution
-            r_l = cbrt(q_liq / N_liq / (4 / 3 * π) / ρ_liq * ρ_air)
+            r_l = cbrt(q_liq / N_liq / (4 / 3 * FT(π)) / ρ_liq * ρ_air)
         elseif "Gamma" in droplet_size_distribution
-            λ = cbrt(32 * π * N_liq / q_liq * ρ_liq / ρ_air)
+            λ = cbrt(32 * FT(π) * N_liq / q_liq * ρ_liq / ρ_air)
             #A = N_liq* λ^2
             r_l = 2 / λ
         end
-        A_aer = 4 * π * r_l^2
+        A_aer = 4 * FT(π) * r_l^2
         dN_act_dt_immersion = max(FT(0), J_immersion * N_liq * A_aer)
-        dqi_dt_new_immers = dN_act_dt_immersion * 4 / 3 * π * r_l^3 * ρ_ice / ρ_air
+        dqi_dt_new_immers = dN_act_dt_immersion * 4 / 3 * FT(π) * r_l^3 * ρ_ice / ρ_air
     end
 
     dN_act_dt_hom = FT(0)
@@ -113,33 +112,29 @@ function parcel_model(dY, Y, p, t)
         Δa_w = CMO.a_w_eT(tps, e, T) - a_w_ice
         if "Monodisperse" in droplet_size_distribution
             J_hom = CMI_hom.homogeneous_J(ip.homogeneous, Δa_w)
-            r_l = cbrt(q_liq / N_liq / (4 / 3 * π) / ρ_liq * ρ_air)
-            V_l = 4 /3 * π * r_l^3
+            r_l = cbrt(q_liq / N_liq / (4 / 3 * FT(π)) / ρ_liq * ρ_air)
+            V_l = 4 /3 * FT(π) * r_l^3
             dN_act_dt_hom = max(FT(0), J_hom * N_liq * V_l)
-            dqi_dt_new_hom = dN_act_dt_hom * 4 / 3 * π * r_l^3 * ρ_ice / ρ_air
+            dqi_dt_new_hom = dN_act_dt_hom * 4 / 3 * FT(π) * r_l^3 * ρ_ice / ρ_air
         elseif "Gamma" in droplet_size_distribution
             J_hom = CMI_hom.homogeneous_J(ip.homogeneous, Δa_w)
-            λ = cbrt(32 * π * N_liq / q_liq * ρ_liq / ρ_air)
+            λ = cbrt(32 * FT(π) * N_liq / q_liq * ρ_liq / ρ_air)
             #A = N_liq* λ^2
             r_l = 2 / λ
-            V_l = 4 / 3 * π * r_l^3
+            V_l = 4 / 3 * FT(π) * r_l^3
             dN_act_dt_hom = max(FT(0), J_hom * N_aer * V_l)
-            dqi_dt_new_hom = dN_act_dt_hom * 4 / 3 * π * r_l^3 * ρ_ice / ρ_air
+            dqi_dt_new_hom = dN_act_dt_hom * 4 / 3 * FT(π) * r_l^3 * ρ_ice / ρ_air
         elseif "Lognormal" in droplet_size_distribution
             J_hom = CMI_hom.homogeneous_J(ip.homogeneous, Δa_w)
-            dN_act_dt_hom = max(FT(0), J_hom * N_liq * 4 / 3 * π * exp((6*log(R_mode_liq) + 9 * σ^2)/(2)))
+            dN_act_dt_hom = max(FT(0), J_hom * N_liq * 4 / 3 * FT(π) * exp((6*log(R_mode_liq) + 9 * σ^2)/(2)))
             r_l = R_mode_liq * exp(σ)
-            dqi_dt_new_hom = dN_act_dt_hom * 4 / 3 * π * r_l^3 * ρ_ice / ρ_air
+            dqi_dt_new_hom = dN_act_dt_hom * 4 / 3 * FT(π) * r_l^3 * ρ_ice / ρ_air
         elseif "Jensen_Example" in droplet_size_distribution
-            if Δa_w < ip.homogeneous.Δa_w_min
-                Δa_w = ip.homogeneous.Δa_w_min
-            elseif Δa_w > ip.homogeneous.Δa_w_max
-                Δa_w = ip.homogeneous.Δa_w_max
-            end
+            Δa_w = ifelse(Δa_w < ip.homogeneous.Δa_w_min, ip.homogeneous.Δa_w_min, ifelse(Δa_w > ip.homogeneous.Δa_w_max, ip.homogeneous.Δa_w_max, Δa_w))
             J_hom = CMI_hom.homogeneous_J(ip.homogeneous, Δa_w)
-            dN_act_dt_hom = max(FT(0), J_hom * N_aer * 4 / 3 * π * exp((6*log(r_nuc) + 9 * σ^2)/(2)))
+            dN_act_dt_hom = max(FT(0), J_hom * N_aer * 4 / 3 * FT(π) * exp((6*log(r_nuc) + 9 * σ^2)/(2)))
             r_aero = r_nuc * exp(σ)
-            dqi_dt_new_hom = dN_act_dt_hom * 4 / 3 * π * r_aero^3 * ρ_ice / ρ_air
+            dqi_dt_new_hom = dN_act_dt_hom * 4 / 3 * FT(π) * r_aero^3 * ρ_ice / ρ_air
         end
     end
 
@@ -158,19 +153,19 @@ function parcel_model(dY, Y, p, t)
         G_i = CMO.G_func(aps, tps, T, TD.Ice())
         if "Monodisperse" in droplet_size_distribution
             # Deposition on existing crystals (assuming all are the same...)
-            r_i = cbrt(q_ice / N_ice / (4 / 3 * π) / ρ_ice * ρ_air)
+            r_i = cbrt(q_ice / N_ice / (4 / 3 * FT(π)) / ρ_ice * ρ_air)
             C_i = r_i
-            dqi_dt_depo = 4 * π / ρ_air * (s_i - 1) * G_i * r_i * N_ice
+            dqi_dt_depo = 4 * FT(π) / ρ_air * (S_i - 1) * G_i * r_i * N_ice
         elseif "Gamma" in droplet_size_distribution
         # Condensation on existing droplets assuming n(r) = A r exp(-λr)
-            λ = cbrt(32 * π * N_ice / q_ice * ρ_ice / ρ_air)
+            λ = cbrt(32 * FT(π) * N_ice / q_ice * ρ_ice / ρ_air)
             #A = N_ice* λ^2
             r_i = 2 / λ
             C_i = r_i
-            dqi_dt_depo = 4 * π / ρ_air * (s_i - 1) * G_i * r_i * N_ice
+            dqi_dt_depo = 4 * FT(π) / ρ_air * (S_i - 1) * G_i * r_i * N_ice
         elseif "Lognormal" in droplet_size_distribution || "Jensen_Example" in droplet_size_distribution
-            r_i = cbrt(q_ice / N_ice / (4 / 3 * π) / ρ_ice * ρ_air) * exp(σ)
-            dqi_dt_depo = 4 * π / ρ_air * (s_i - 1) * G_i * r_i * N_ice
+            r_i = cbrt(q_ice / N_ice / (4 / 3 * FT(π)) / ρ_ice * ρ_air) * exp(σ)
+            dqi_dt_depo = 4 * FT(π) / ρ_air * (S_i - 1) * G_i * r_i * N_ice
         end
         # TODO - check if r_i is non-zero if initial conditions say q_ice = 0
     end
@@ -180,14 +175,14 @@ function parcel_model(dY, Y, p, t)
         G_l = CMO.G_func(aps, tps, T, TD.Liquid())
         if "Monodisperse" in droplet_size_distribution
         # Condensation on existing droplets assuming all are the same
-            r_l = cbrt(q_liq / N_liq / (4 / 3 * π) / ρ_liq * ρ_air)
+            r_l = cbrt(q_liq / N_liq / (4 / 3 * FT(π)) / ρ_liq * ρ_air)
         elseif "Gamma" in droplet_size_distribution
         # Condensation on existing droplets assuming n(r) = A r exp(-λr)
-            λ = cbrt(32 * π * N_liq / q_liq * ρ_liq / ρ_air)
+            λ = cbrt(32 * FT(π) * N_liq / q_liq * ρ_liq / ρ_air)
             #A = N_liq* λ^2
             r_l = 2 / λ
         end
-        dql_dt_cond = 4 * π / ρ_air * (s_liq - 1) * G_l * r_l * N_liq
+        dql_dt_cond = 4 * FT(π) / ρ_air * (S_liq - 1) * G_l * r_l * N_liq
     end
 
     dq_liq_dt_vap_to_liq = dql_dt_cond   # change in q_liq from vap-liq transitions
@@ -207,12 +202,12 @@ function parcel_model(dY, Y, p, t)
     dq_ice_dt = dq_ice_dt_vap_to_ice + dq_ice_dt_liq_to_ice
     dq_liq_dt = dq_liq_dt_vap_to_liq - dq_ice_dt_liq_to_ice
     dq_vap_dt = -dq_ice_dt - dq_liq_dt
-    ds_l_dt = a1 * w * s_liq - (a2 + a3) * s_liq * dq_liq_dt_vap_to_liq - (a2 + a4) * s_liq * dq_ice_dt_vap_to_ice - a5 * s_liq * dq_ice_dt_liq_to_ice
+    dS_l_dt = a1 * w * S_liq - (a2 + a3) * S_liq * dq_liq_dt_vap_to_liq - (a2 + a4) * S_liq * dq_ice_dt_vap_to_ice - a5 * S_liq * dq_ice_dt_liq_to_ice
     dp_a_dt = -p_a * grav / R_a / T * w
     dT_dt = -grav / cp_a * w + L_vap / cp_a * dq_liq_dt_vap_to_liq + L_fus / cp_a * dq_ice_dt_liq_to_ice + L_subl / cp_a * dq_ice_dt_vap_to_ice
 
     # Set tendencies
-    dY[1] = ds_l_dt        # saturation ratio over liquid water
+    dY[1] = dS_l_dt        # saturation ratio over liquid water
     dY[2] = dp_a_dt        # pressure
     dY[3] = dT_dt          # temperature
     dY[4] = dq_vap_dt      # vapor specific humidity
@@ -234,13 +229,13 @@ Returns the solution of an ODE probelm defined by the parcel model.
 
 Inputs:
  - IC - A vector with the initial conditions for
-   [s_l, p_a, T, q_vap, q_liq, q_ice, N_aer, N_liq, N_ice, x_sulph]
+   [S_l, p_a, T, q_vap, q_liq, q_ice, N_aer, N_liq, N_ice, x_sulph]
  - t_0 - simulation start time
  - t_end - simulation end time
  - p - a named tuple with simulation parameters.
 
 Initial condition contains (all in base SI units):
- - s_l - saturation ratio over liquid water,
+ - S_l - saturation ratio over liquid water,
  - p_a - atmospheric pressure
  - T - temperature
  - q_vap - water vapor specific humidity