Refactor parcel

docs/src/IceNucleationParcel0D.md

@@ -205,7 +205,7 @@ There are multiple ways of running deposition nucleation in the parcel.
   above this value will result in nucleation in a different mode.
   `"ActivityBasedDeposition"` will trigger a water activity based approach
   from [Alpert2022](@cite). In this approach, ice production rate ``P_{ice, depo}``
-  is calculated from 
+  is calculated from
 ```math
 \begin{equation}
   P_{ice, depo} = \left[ \frac{dN_i}{dt} \right]_{depo} = J_{depo}\;A_{aero}\;N_{aero}
@@ -252,14 +252,14 @@ Here we show various example simulation results from the adiabatic parcel
   and homogeneous freezing with deposition growth.
 
 We start with deposition freezing on dust.
-The model is run three times using the `"MohlerAF_Deposition"` approach 
+The model is run three times using the `"MohlerAF_Deposition"` approach
   for 30 minutes simulation time, (shown by three different colors on the plot).
 Between each run the water vapor specific humidity is changed,
   while keeping all other state variables the same as at the last time step
   of the previous run.
 The prescribed vertical velocity is equal to 3.5 cm/s.
 Supersaturation is plotted for both liquid (solid lines) and ice (dashed lines).
-The pale blue line uses the `"MohlerRate_Deposition"` approach. 
+The pale blue line uses the `"MohlerRate_Deposition"` approach.
   We only run it for the first GCM timestep because the rate approach requires
   the change in ice saturation over time. With the discontinuous jump in saturation,
   the parameterization is unable to determine a proper nucleation rate. When we force
@@ -268,7 +268,7 @@ The pale blue line uses the `"MohlerRate_Deposition"` approach.
   in the `"MohlerAF_Deposition"` approach for the first GCM timestep.
 
 ```@example
-include("../../parcel/Tully_et_al_2023.jl")
+include("../../parcel/Example_Tully_et_al_2023.jl")
 ```
 ![](cirrus_box.svg)
 
@@ -280,7 +280,7 @@ The water activity based parameterization for deposition nucleation shows
   is no common aerosol type between the two parameterizations.
 
 ```@example
-include("../../parcel/Deposition_Nucleation.jl")
+include("../../parcel/Example_Deposition_Nucleation.jl")
 ```
 ![](deposition_nucleation.svg)
 
@@ -289,7 +289,7 @@ In the plots below, the parcel model is ran with only condensation (no ice or fr
 It is compared to [Rogers1975](@cite).
 
 ```@example
-include("../../parcel/Liquid_only.jl")
+include("../../parcel/Example_Liquid_only.jl")
 ```
 ![](liquid_only_parcel.svg)
 
@@ -299,7 +299,7 @@ The plots below are the results of the adiabatic parcel model
   yet been validated against literature.
 
 ```@example
-include("../../parcel/Immersion_Freezing.jl")
+include("../../parcel/Example_Immersion_Freezing.jl")
 ```
 ![](immersion_freezing.svg)
 
@@ -312,13 +312,13 @@ The following plots show the parcel model running with homogeneous freezing and
   for demonstrative purposes.
 
 ```@example
-include("../../parcel/Jensen_et_al_2022.jl")
+include("../../parcel/Example_Jensen_et_al_2022.jl")
 ```
 ![](Jensen_et_al_2022.svg)
 
 ## P3 Ice Nucleation Parameterizations
 The parcel also includes ice nucleation parameterizations used in
-  the P3 scheme as described in [MorrisonMilbrandt2015](@cite). 
+  the P3 scheme as described in [MorrisonMilbrandt2015](@cite).
   Deposition nucleation is based on the ice crystal number parameterization
   from Cooper (1986). The heterogeneous freezing parameterization, which
   follows Bigg(1953) with parameters from Barklie aand Gokhale (1959), is
@@ -330,6 +330,6 @@ Shown below are three separate parcel simulations for deposition nucleation,
   INP number, while the other two freezing modes do. Updraft velocity is
   set to 0.5 m/s.
 ```@example
-include("../../parcel/P3_ice_nuc.jl")
+include("../../parcel/Example_P3_ice_nuc.jl")
 ```
 ![](P3_ice_nuc.svg)

parcel/Deposition_Nucleation.jl → parcel/Example_Deposition_Nucleation.jl

@@ -5,18 +5,11 @@ import CloudMicrophysics as CM
 import CLIMAParameters as CP
 
 # definition of the ODE problem for parcel model
-include(joinpath(pkgdir(CM), "parcel", "parcel.jl"))
+include(joinpath(pkgdir(CM), "parcel", "Parcel.jl"))
 FT = Float32
+
 # get free parameters
 tps = TD.Parameters.ThermodynamicsParameters(FT)
-aps = CMP.AirProperties(FT)
-wps = CMP.WaterProperties(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(2000 * 1e3)
@@ -34,60 +27,45 @@ 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)
+Rᵥ = TD.Parameters.R_v(tps)
 eₛ = TD.saturation_vapor_pressure(tps, T₀, TD.Liquid())
-e = qᵥ * p₀ * R_v / Rₐ
+e = eᵥ(qᵥ, p₀, Rₐ, Rᵥ)
 
 Sₗ = FT(e / eₛ)
 IC = [Sₗ, p₀, T₀, qᵥ, qₗ, qᵢ, Nₐ, Nₗ, Nᵢ, x_sulph]
 
-ξ(T) =
-    TD.saturation_vapor_pressure(tps, T, TD.Liquid()) /
-    TD.saturation_vapor_pressure(tps, T, TD.Ice())
-S_i(T, S_liq) = ξ(T) * S_liq
-
 # Simulation parameters passed into ODE solver
 r_nuc = FT(1.25e-6)                     # assumed size of nucleated particles
 w = FT(3.5 * 1e-2)                      # updraft speed
-α_m = FT(0.5)                           # accomodation coefficient
 const_dt = FT(0.1)                      # model timestep
 t_max = FT(100)
-aerosol_1 = [CMP.DesertDust(FT), CMP.ArizonaTestDust(FT)]                     # aersol type for DustDeposition
-aerosol_2 = [CMP.Feldspar(FT), CMP.Ferrihydrite(FT), CMP.Kaolinite(FT)]   # aersol type for DepositionNucleation
-ice_nucleation_modes_list =
-    [["MohlerRate_Deposition"], ["ActivityBasedDeposition"]]
-growth_modes = ["Deposition"]
-droplet_size_distribution_list = [["Monodisperse"]]
+aerosol_1 = [CMP.DesertDust(FT), CMP.ArizonaTestDust(FT)]               # aersol type for DustDeposition
+aerosol_2 = [CMP.Feldspar(FT), CMP.Ferrihydrite(FT), CMP.Kaolinite(FT)] # aersol type for DepositionNucleation
+deposition_modes = ["MohlerRate", "ABDINM"]
+deposition_growth = "Deposition"
+size_distribution = "Monodisperse"
 
 # Plotting
-fig = MK.Figure(resolution = (800, 600))
+fig = MK.Figure(size = (800, 600))
 ax1 = MK.Axis(fig[1, 1], ylabel = "Ice Saturation [-]")
 ax2 = MK.Axis(fig[1, 2], ylabel = "Temperature [K]")
 ax3 = MK.Axis(fig[2, 1], ylabel = "q_ice [g/kg]", xlabel = "time")
 ax4 = MK.Axis(fig[2, 2], ylabel = "N_ice [m^-3]", xlabel = "time")
 
-for ice_nucleation_modes in ice_nucleation_modes_list
-    nuc_mode = ice_nucleation_modes[1]
-    droplet_size_distribution = droplet_size_distribution_list[1]
-
-    if nuc_mode == "MohlerRate_Deposition"
+for deposition in deposition_modes
+    if deposition == "MohlerRate"
         for aerosol in aerosol_1
-            p = (;
-                wps,
-                aps,
-                tps,
-                ip,
-                const_dt,
-                r_nuc,
-                w,
-                α_m,
-                aerosol,
-                ice_nucleation_modes,
-                growth_modes,
-                droplet_size_distribution,
+            local params = parcel_params{FT}(
+                const_dt = const_dt,
+                r_nuc = r_nuc,
+                w = w,
+                aerosol = aerosol,
+                deposition = deposition,
+                deposition_growth = deposition_growth,
+                size_distribution = size_distribution,
             )
-
             # solve ODE
-            sol = run_parcel(IC, FT(0), t_max, p)
+            local sol = run_parcel(IC, FT(0), t_max, params)
 
             # Plot results
             if aerosol == CMP.DesertDust(FT)
@@ -98,33 +76,27 @@ for ice_nucleation_modes in ice_nucleation_modes_list
             MK.lines!(                               # saturation
                 ax1,
                 sol.t,
-                S_i.(sol[3, :], sol[1, :]),
+                S_i.(tps, sol[3, :], sol[1, :]),
                 label = aero_label,
             )
             MK.lines!(ax2, sol.t, sol[3, :])         # temperature
             MK.lines!(ax3, sol.t, sol[6, :] * 1e3)   # q_ice
             MK.lines!(ax4, sol.t, sol[9, :])         # N_ice
         end
 
-    elseif nuc_mode == "ActivityBasedDeposition"
+    elseif deposition == "ABDINM"
         for aerosol in aerosol_2
-            p = (;
-                wps,
-                aps,
-                tps,
-                ip,
-                const_dt,
-                r_nuc,
-                w,
-                α_m,
-                aerosol,
-                ice_nucleation_modes,
-                growth_modes,
-                droplet_size_distribution,
+            local params = parcel_params{FT}(
+                const_dt = const_dt,
+                r_nuc = r_nuc,
+                w = w,
+                aerosol = aerosol,
+                deposition = deposition,
+                deposition_growth = deposition_growth,
+                size_distribution = size_distribution,
             )
-
             # solve ODE
-            sol = run_parcel(IC, FT(0), t_max, p)
+            local sol = run_parcel(IC, FT(0), t_max, params)
 
             # Plot results
             if aerosol == CMP.Feldspar(FT)
@@ -137,7 +109,7 @@ for ice_nucleation_modes in ice_nucleation_modes_list
             MK.lines!(                                                      # saturation
                 ax1,
                 sol.t,
-                S_i.(sol[3, :], sol[1, :]),
+                S_i.(tps, sol[3, :], sol[1, :]),
                 label = aero_label,
                 linestyle = :dash,
             )

parcel/Immersion_Freezing.jl → parcel/Example_Immersion_Freezing.jl

@@ -5,20 +5,14 @@ import CloudMicrophysics as CM
 import CLIMAParameters as CP
 
 # definition of the ODE problem for parcel model
-include(joinpath(pkgdir(CM), "parcel", "parcel.jl"))
+include(joinpath(pkgdir(CM), "parcel", "Parcel.jl"))
 FT = Float32
 # get free parameters
 tps = TD.Parameters.ThermodynamicsParameters(FT)
-aps = CMP.AirProperties(FT)
 wps = CMP.WaterProperties(FT)
-ip = CMP.IceNucleationParameters(FT)
-
-# Constants
-ρₗ = wps.ρw
-R_v = TD.Parameters.R_v(tps)
-R_d = TD.Parameters.R_d(tps)
 
 # Initial conditions
+ρₗ = wps.ρw
 Nₐ = FT(0)
 Nₗ = FT(500 * 1e3)
 Nᵢ = FT(0)
@@ -31,29 +25,26 @@ qᵢ = FT(0)
 x_sulph = FT(0.01)
 
 # 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)
+q = TD.PhasePartition.(qᵥ + qₗ + qᵢ, qₗ, qᵢ)
+R_v = TD.Parameters.R_v(tps)
 Rₐ = TD.gas_constant_air(tps, q)
 eₛ = TD.saturation_vapor_pressure(tps, T₀, TD.Liquid())
-e = qᵥ * p₀ * R_v / Rₐ
-
+e = eᵥ(qᵥ, p₀, Rₐ, R_v)
 Sₗ = FT(e / eₛ)
 IC = [Sₗ, p₀, T₀, qᵥ, qₗ, qᵢ, Nₐ, Nₗ, Nᵢ, x_sulph]
 
 # Simulation parameters passed into ODE solver
-r_nuc = FT(0.5 * 1.e-4 * 1e-6)             # assumed size of nucleated particles
 w = FT(0.7)                                # updraft speed
-α_m = FT(0.5)                              # accomodation coefficient
 const_dt = FT(1)                           # model timestep
 t_max = FT(1200)
 aerosol = CMP.Illite(FT)
-ice_nucleation_modes = ["ImmersionFreezing"]
-growth_modes = ["Condensation", "Deposition"]
-droplet_size_distribution_list = [["Monodisperse"], ["Gamma"]]
+heterogeneous = "ABIFM"
+condensation_growth = "Condensation"
+deposition_growth = "Deposition"
+size_distribution_list = ["Monodisperse", "Gamma"]
 
 # Plotting
-fig = MK.Figure(resolution = (800, 600))
+fig = MK.Figure(size = (800, 600))
 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]")
@@ -67,47 +58,31 @@ ax6 = MK.Axis(
 )
 MK.ylims!(ax6, 1e-6, 1)
 
-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,
+for DSD in size_distribution_list
+    local params = parcel_params{FT}(
+        const_dt = const_dt,
+        w = w,
+        aerosol = aerosol,
+        heterogeneous = heterogeneous,
+        condensation_growth = condensation_growth,
+        deposition_growth = deposition_growth,
+        size_distribution = DSD,
     )
     # solve ODE
-    sol = run_parcel(IC, FT(0), t_max, p)
-
-    DSD = droplet_size_distribution[1]
+    local sol = run_parcel(IC, FT(0), t_max, params)
 
     # Plot results
-    ξ(T) =
-        TD.saturation_vapor_pressure(tps, T, TD.Liquid()) /
-        TD.saturation_vapor_pressure(tps, T, TD.Ice())
-    S_i(T, S_liq) = ξ(T) * S_liq - 1
-
-    MK.lines!(ax1, sol.t / 60, S_i.(sol[3, :], sol[1, :]), label = DSD)
+    MK.lines!(
+        ax1,
+        sol.t / 60,
+        S_i.(tps, sol[3, :], sol[1, :]) .- 1,
+        label = DSD,
+    )
     MK.lines!(ax2, sol.t / 60, sol[3, :])
     MK.lines!(ax3, sol.t / 60, sol[6, :] * 1e3)
     MK.lines!(ax4, sol.t / 60, sol[5, :] * 1e3)
-
     sol_Nₗ = sol[8, :]
     sol_Nᵢ = sol[9, :]
-    sol_qₗ = sol[5, :]
-    if DSD == "Monodisperse"
-        rₗ = cbrt.(sol_qₗ ./ sol_Nₗ ./ (4 / 3 * FT(π)) / ρₗ * ρₐ)
-    end
-    if DSD == "Gamma"
-        λ = cbrt.(32 .* FT(π) .* sol_Nₗ ./ sol_qₗ * ρₗ / ρₐ)
-        rₗ = 2 ./ λ
-    end
     MK.lines!(ax5, sol.t / 60, sol_Nₗ)
     MK.lines!(ax6, sol.t / 60, sol_Nᵢ ./ (sol_Nₗ .+ sol_Nᵢ))
 end
@@ -118,7 +93,7 @@ MK.axislegend(
     labelsize = 12,
     orientation = :horizontal,
     nbanks = 2,
-    position = :rb,
+    position = :rt,
 )
 
 MK.save("immersion_freezing.svg", fig)