Merge branch 'main' of https://github.com/CliMA/CloudMicrophysics.jl …

…into al/p3_ice_nuc

parcel/Immersion_Freezing.jl

@@ -6,7 +6,7 @@ import CLIMAParameters as CP
 
 # definition of the ODE problem for parcel model
 include(joinpath(pkgdir(CM), "parcel", "parcel.jl"))
-FT = Float64
+FT = Float32
 # get free parameters
 tps = TD.Parameters.ThermodynamicsParameters(FT)
 aps = CMP.AirProperties(FT)
@@ -26,7 +26,7 @@ r₀ = FT(1e-6)
 p₀ = FT(800 * 1e2)
 T₀ = FT(251)
 qᵥ = FT(8.1e-4)
-qₗ = Nₗ * 4 / 3 * FT(π) * r₀^3 * ρₗ / 1.2 # 1.2 should be ρₐ
+qₗ = Nₗ * 4 / 3 * FT(π) * r₀^3 * ρₗ / FT(1.2) # 1.2 should be ρₐ
 qᵢ = FT(0)
 x_sulph = FT(0.01)
 

parcel/Jensen_et_al_2022.jl

@@ -9,7 +9,7 @@ import CloudMicrophysics.Parameters as CMP
 # definition of the ODE problem for parcel model
 include(joinpath(pkgdir(CM), "parcel", "parcel.jl"))
 
-FT = Float64
+FT = Float32
 
 # Get free parameters
 tps = TD.Parameters.ThermodynamicsParameters(FT)
@@ -42,7 +42,7 @@ x_sulph = FT(0)
 # saturation and partial pressure
 eₛ = TD.saturation_vapor_pressure(tps, T₀, TD.Liquid())
 qᵥ = ϵₘ / (ϵₘ - 1 + 1 / cᵥ₀)
-qₗ = Nₗ * 4 / 3 * FT(π) * r₀^3 * ρₗ / 1.2
+qₗ = Nₗ * 4 / 3 * FT(π) * r₀^3 * ρₗ / FT(1.2)
 qᵢ = FT(0)
 q = TD.PhasePartition(qᵥ + qₗ + qᵢ, qₗ, qᵢ)
 R_a = TD.gas_constant_air(tps, q)

parcel/Liquid_only.jl

@@ -7,7 +7,7 @@ import CLIMAParameters as CP
 # definition of the ODE problem for parcel model
 include(joinpath(pkgdir(CM), "parcel", "parcel.jl"))
 
-FT = Float64
+FT = Float32
 
 # Get free parameters
 tps = TD.Parameters.ThermodynamicsParameters(FT)

parcel/Tully_et_al_2023.jl

@@ -203,4 +203,4 @@ function Tully_et_al_2023(FT)
     end
     MK.save("cirrus_box.svg", fig)
 end
-Tully_et_al_2023(Float64)
+Tully_et_al_2023(Float32)

parcel/parcel.jl

@@ -36,7 +36,7 @@ function parcel_model(dY, Y, p, t)
     N_aer = Y[7]      # number concentration of interstitial aerosol
     N_liq = Y[8]      # number concentration of existing water droplets
     N_ice = Y[9]      # number concentration of activated ice crystals
-    x_sulph = Y[10]   # percent mass sulphuric acid    
+    x_sulph = Y[10]   # percent mass sulphuric acid
 
     # Constants
     R_v = TD.Parameters.R_v(tps)
@@ -404,8 +404,8 @@ function run_parcel(IC, t_0, t_end, p)
         problem,
         timestepper,
         dt = const_dt,
-        reltol = 10 * eps(FT),
-        abstol = 10 * eps(FT),
+        reltol = 100 * eps(FT),
+        abstol = 100 * eps(FT),
     )
     return sol
 end

src/AerosolActivation.jl

@@ -165,13 +165,13 @@ function max_supersaturation(
 
         mode_i = ad.Modes[i]
 
-        f::FT = 0.5 * exp(2.5 * (log(mode_i.stdev))^2)
-        g::FT = 1 + 0.25 * log(mode_i.stdev)
+        f::FT = ap.f1 * exp(ap.f2 * (log(mode_i.stdev))^2)
+        g::FT = ap.g1 + ap.g2 * log(mode_i.stdev)
         η::FT = (α * w / G)^FT(3 / 2) / (FT(2 * pi) * ap.ρ_w * γ * mode_i.N)
 
         tmp +=
             1 / (Sm[i])^2 *
-            (f * (ζ / η)^FT(3 / 2) + g * (Sm[i]^2 / (η + 3 * ζ))^FT(3 / 4))
+            (f * (ζ / η)^ap.p1 + g * (Sm[i]^2 / (η + 3 * ζ))^ap.p2)
     end
 
     return FT(1) / sqrt(tmp)

src/parameters/AerosolActivation.jl

@@ -20,6 +20,18 @@ Base.@kwdef struct AerosolActivationParameters{FT} <: ParametersType{FT}
     σ::FT
     "gravitational_acceleration [m/s2]"
     g::FT
+    "scaling coefficient in Abdul-Razzak and Ghan 2000 [-]"
+    f1::FT
+    "scaling coefficient in Abdul-Razzak and Ghan 2000 [-]"
+    f2::FT
+    "scaling coefficient in Abdul-Razzak and Ghan 2000 [-]"
+    g1::FT
+    "scaling coefficient in Abdul-Razzak and Ghan 2000 [-]"
+    g2::FT
+    "power of (zeta / eta) in Abdul-Razzak and Ghan 2000 [-]"
+    p1::FT
+    "power of (S_m^2 / (zeta + 3 * eta)) in Abdul-Razzak and Ghan 2000 [-]"
+    p2::FT
 end
 
 AerosolActivationParameters(::Type{FT}) where {FT <: AbstractFloat} =
@@ -32,6 +44,12 @@ function AerosolActivationParameters(td::CP.AbstractTOMLDict)
         :density_liquid_water => :ρ_w,
         :surface_tension_water => :σ,
         :gravitational_acceleration => :g,
+        :ARG2000_f_coeff_1 => :f1,
+        :ARG2000_f_coeff_2 => :f2,
+        :ARG2000_g_coeff_1 => :g1,
+        :ARG2000_g_coeff_2 => :g2,
+        :ARG2000_pow_1 => :p1,
+        :ARG2000_pow_2 => :p2,
     )
     parameters = CP.get_parameter_values(td, name_map, "CloudMicrophysics")
     FT = CP.float_type(td)

src/parameters/toml/ARG2000.toml

@@ -0,0 +1,17 @@
+[ARG2000_f_coeff_1]
+value = 0.26583888195264627
+
+[ARG2000_f_coeff_2]
+value = 2.3851515425961853
+
+[ARG2000_g_coeff_1]
+value = 0.779519468021862
+
+[ARG2000_g_coeff_2]
+value = 0.10571967167118024
+
+[ARG2000_pow_1]
+value = 1.6523365679298359
+
+[ARG2000_pow_2]
+value = 0.7578626397779737