update with main

CliMA · Oct 24, 2023 · e5202b5 · e5202b5
2 parents 4514699 + b5ad6da
commit e5202b5
Showing 1 changed file with 30 additions and 7 deletions.
diff --git a/parcel/parcel.jl b/parcel/parcel.jl
@@ -1,8 +1,8 @@
 import Thermodynamics as TD
-import CloudMicrophysics as CM
 
 import CloudMicrophysics.Common as CMO
 import CloudMicrophysics.HetIceNucleation as CMI_het
+import CloudMicrophysics.HomIceNucleation as CMI_hom
 import CloudMicrophysics.Parameters as CMP
 
 #! format: off
@@ -87,13 +87,13 @@ function parcel_model(dY, Y, p, t)
             CMO.a_w_xT(H2SO4_prs, tps, x_sulph, T) - CMO.a_w_ice(tps, 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 && "ImmersionFreezing" in ice_nucleation_modes
+        if "Monodisperse" in droplet_size_distribution
             r_l = cbrt(q_liq / N_liq / (4 / 3 * π) / ρ_liq * ρ_air)
             A_aer = 4 * π * 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
         end
-        if "Gamma" in droplet_size_distribution && "ImmersionFreezing" in ice_nucleation_modes
+        if "Gamma" in droplet_size_distribution
             λ = cbrt(32 * π * N_liq / q_liq * ρ_liq / ρ_air)
             #A = N_liq* λ^2
             r_l = 2 / λ
@@ -103,9 +103,32 @@ function parcel_model(dY, Y, p, t)
         end
     end
 
-    dN_ice_dt = dN_act_dt_depo + dN_act_dt_immersion
+    dN_act_dt_hom = FT(0)
+    dqi_dt_new_hom = FT(0)
+    if "HomogeneousFreezing" in ice_nucleation_modes
+        Δa_w = T > FT(185) && T < FT(235) ?
+            CMO.a_w_xT(H2SO4_prs, tps, x_sulph, T) - CMO.a_w_ice(tps, T) :
+            CMO.a_w_eT(tps, e, T) - CMO.a_w_ice(tps, T)
+        J_hom = CMI_hom.homogeneous_J(ip, Δa_w)
+        if "Monodisperse" in droplet_size_distribution
+            r_l = cbrt(q_liq / N_liq / (4 / 3 * π) / ρ_liq * ρ_air)
+            V_l = 4 / 3 * π * 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
+        end
+        if "Gamma" in droplet_size_distribution
+            λ = cbrt(32 * π * N_liq / q_liq * ρ_liq / ρ_air)
+            #A = N_liq* λ^2
+            r_l = 2 / λ
+            V_l = 4 /3 * π * 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
+        end
+    end
+
+    dN_ice_dt = dN_act_dt_depo + dN_act_dt_immersion + dN_act_dt_hom
     dN_aer_dt = -dN_act_dt_depo
-    dN_liq_dt = -dN_act_dt_immersion
+    dN_liq_dt = -dN_act_dt_immersion - dN_act_dt_hom
 
     # Growth
     dqi_dt_depo = FT(0)
@@ -136,8 +159,8 @@ function parcel_model(dY, Y, p, t)
     end
 
     # Update the tendecies
-    dq_ice_dt = dqi_dt_new_depo + dqi_dt_depo + dqi_dt_new_immers
-    dq_liq_dt = dql_dt_cond - dqi_dt_new_immers
+    dq_ice_dt = dqi_dt_new_depo + dqi_dt_depo + dqi_dt_new_immers + dqi_dt_new_hom
+    dq_liq_dt = dql_dt_cond - dqi_dt_new_immers - dqi_dt_new_hom
     dS_l_dt = a1 * w * S_liq - (a2 + a3) * S_liq * dq_liq_dt - (a2 + a4) * S_liq * dq_ice_dt
     dp_a_dt = -p_a * grav / R_a / T * w
     dT_dt = -grav / cp_a * w + L_vap / cp_a * dq_liq_dt + L_subl / cp_a * dq_ice_dt