Tully works

parcel/Tully_et_al_2023.jl

@@ -45,9 +45,6 @@ function Tully_et_al_2023(FT)
 
     # get free parameters
     tps = TD.Parameters.ThermodynamicsParameters(FT)
-    aps = CMP.AirProperties(FT)
-    wps = CMP.WaterProperties(FT)
-    ip = CMP.IceNucleationParameters(FT)
     # Initial conditions for 1st period
     N_aerosol = FT(2000 * 1e3)
     N_droplets = FT(0)
@@ -69,15 +66,13 @@ function Tully_et_al_2023(FT)
     t_max = 30 * 60
 
     # Simulation parameters passed into ODE solver
-    r_nuc = FT(0.5 * 1.e-4 * 1e-6)             # assumed size of nucleated particles
-    w = FT(3.5 * 1e-2)                         # updraft speed
-    α_m = FT(0.5)                              # accomodation coefficient
-    const_dt = 0.1                             # model timestep
-    aerosol = CMP.DesertDust(FT)               # aerosol type
-    ice_nucleation_modes_list =                # ice nucleation modes
-        [["MohlerAF_Deposition"], ["MohlerRate_Deposition"]]
-    growth_modes = ["Deposition"]              # switch on deposition growth
-    droplet_size_distribution = ["Monodisperse"]
+    r_nuc = FT(0.5 * 1.e-4 * 1e-6)                    # assumed size of nucleated particles
+    w = FT(3.5 * 1e-2)                                # updraft speed
+    const_dt = 0.1                                    # model timestep
+    aerosol = CMP.DesertDust(FT)                      # aerosol type
+    ice_nucleation_modes = ["MohlerAF", "MohlerRate"] # ice nucleation modes
+    growth_modes = ["Deposition"]                     # switch on deposition growth
+    size_distribution = "Monodisperse"
 
     # Plots
     fig = MK.Figure(resolution = (800, 600))
@@ -96,22 +91,16 @@ function Tully_et_al_2023(FT)
         TD.saturation_vapor_pressure(tps, T, TD.Ice())
     S_i(T, S_liq) = ξ(T) * S_liq - 1
 
-    for ice_nucleation_modes in ice_nucleation_modes_list
-        p = (;
-            wps,
-            aps,
-            tps,
-            ip,
-            const_dt,
-            r_nuc,
-            w,
-            α_m,
-            aerosol,
-            ice_nucleation_modes,
-            growth_modes,
-            droplet_size_distribution,
+    for mode in ice_nucleation_modes
+        params = parcel_params{FT}(
+            const_dt = const_dt,
+            r_nuc = r_nuc,
+            w = w,
+            aerosol = aerosol,
+            deposition = mode,
+            growth_modes = growth_modes,
+            size_distribution = size_distribution,
         )
-
         # Simulation 1
         IC1 = get_initial_condition(
             tps,
@@ -125,8 +114,8 @@ function Tully_et_al_2023(FT)
             N_0,
             x_sulph,
         )
-        sol1 = run_parcel(IC1, 0, t_max, p)
-        if ice_nucleation_modes == ["MohlerAF_Deposition"]
+        sol1 = run_parcel(IC1, 0, t_max, params)
+        if mode == "MohlerAF"
             # Simulation 2
             # (alternatively set T and take q_vap from the previous simulation)
             #IC2 = get_initial_condition(sol1[2, end], sol1[3, end], T2, sol1[5, end], 0.0, sol1[6, end], sol1[7, end])
@@ -143,7 +132,7 @@ function Tully_et_al_2023(FT)
                 sol1[9, end],
                 x_sulph,
             )
-            sol2 = run_parcel(IC2, sol1.t[end], sol1.t[end] + t_max, p)
+            sol2 = run_parcel(IC2, sol1.t[end], sol1.t[end] + t_max, params)
 
             # Simulation 3
             # (alternatively set T and take q_vap from the previous simulation)
@@ -161,7 +150,7 @@ function Tully_et_al_2023(FT)
                 sol2[9, end],
                 x_sulph,
             )
-            sol3 = run_parcel(IC3, sol2.t[end], sol2.t[end] + t_max, p)
+            sol3 = run_parcel(IC3, sol2.t[end], sol2.t[end] + t_max, params)
 
             # Plot results
             sol = [sol1, sol2, sol3]
@@ -184,7 +173,7 @@ function Tully_et_al_2023(FT)
             end
             #! format: on
 
-        elseif ice_nucleation_modes == ["MohlerRate_Deposition"]
+        elseif mode == "MohlerRate"
             MK.lines!(ax1, sol1.t * w, sol1[1, :] .- 1, color = :lightblue)
             MK.lines!(ax2, sol1.t * w, sol1[3, :], color = :lightblue)
             MK.lines!(ax3, sol1.t * w, sol1[9, :] * 1e-3, color = :lightblue)

parcel/parcel.jl

@@ -47,6 +47,7 @@ struct Lognormal{FT} <: CMP.ParametersType{FT}
 end
 # Size distributiom moments
 function distribution_moments(::Monodisperse, qₗ, Nₗ, ρₗ, ρ_air, qᵢ, Nᵢ, ρᵢ)
+    FT = typeof(qₗ)
     # liquid droplet
     if Nₗ == FT(0)
         rₗ = FT(0)
@@ -141,16 +142,16 @@ function deposition_nucleation(params::MohlerAF, state, dY)
     (; ips, aerosol, const_dt, tps) = params
     (; Sₗ, T, Nₐ, Nᵢ) = state
     Sᵢ = ξ(tps, T) * Sₗ
-
-    if S_i >= ips.deposition.Sᵢ_max
+    FT = eltype(state)
+    if Sᵢ >= ips.deposition.Sᵢ_max
         @warn("Supersaturation exceeds Sᵢ_max. No dust will be activated.")
     end
 
     AF = Sᵢ >= ips.deposition.Sᵢ_max ?
         FT(0) :
         AF = CMI_het.dust_activated_number_fraction(
             aerosol,
-            ips.deposition_params,
+            ips.deposition,
             Sᵢ,
             T,
         )
@@ -309,6 +310,7 @@ end
 function deposition(params::DepParams, PSD, state, ρ_air)
     (; aps, tps) = params
     (; T, Sₗ, Nᵢ) = state
+    FT = eltype(state)
     Sᵢ = ξ(tps, T) * Sₗ
     Gᵢ = CMO.G_func(aps, tps, T, TD.Ice())
     return 4 * FT(π) / ρ_air * (Sᵢ - 1) * Gᵢ * PSD.rᵢ * Nᵢ
@@ -360,11 +362,11 @@ function parcel_model(dY, Y, p, t)
     ρ_air = TD.air_density(tps, ts)
 
     # Adiabatic parcel coefficients
-    a1 = L_vap * grav / cp_air / T^2 / R_v - grav / R_air / T
+    a1 = L_vap * grav / cp_air / T^2 / Rᵥ - grav / R_air / T
     a2 = 1 / qᵥ
-    a3 = L_vap^2 / R_v / T^2 / cp_air
-    a4 = L_vap * L_subl / R_v / T^2 / cp_air
-    a5 = L_vap * L_fus / R_v / cp_air / (T^2)
+    a3 = L_vap^2 / Rᵥ / T^2 / cp_air
+    a4 = L_vap * L_subl / Rᵥ / T^2 / cp_air
+    a5 = L_vap * L_fus / Rᵥ / cp_air / (T^2)
 
     # Mean radius, area and volume of liquid droplets and ice crystals
     PSD = distribution_moments(distr, qₗ, Nₗ, ρₗ, ρ_air, qᵢ, Nᵢ, ρᵢ)