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

…into ap/lambdas

p3_sandbox/Project.toml

@@ -0,0 +1,12 @@
+[deps]
+CLIMAParameters = "6eacf6c3-8458-43b9-ae03-caf5306d3d53"
+CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0"
+CloudMicrophysics = "6a9e3e04-43cd-43ba-94b9-e8782df3c71b"
+KernelAbstractions = "63c18a36-062a-441e-b654-da1e3ab1ce7c"
+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.9"

p3_sandbox/p3_sandbox.jl

@@ -0,0 +1,95 @@
+import OrdinaryDiffEq as ODE
+
+import CLIMAParameters
+import CloudMicrophysics.Parameters as CMP
+import CloudMicrophysics.P3Scheme as P3
+import CloudMicrophysics.HetIceNucleation as CMI_het
+import CloudMicrophysics.Common as CMO
+import Thermodynamics as TD
+
+"""
+    ODE problem definitions
+"""
+function p3_sandbox(dY, Y, p, t)
+
+    # Numerical precision used in the simulation
+    FT = eltype(Y)
+    (; const_dt, tps, T, pₐ, qᵥ, qₗ, Nₗ, rₗ, aero_type) = p
+    p3 = CMP.ParametersP3(FT)
+
+    # Our state vector
+    Nᵢ = Y[1]      # Ice number concentration
+    qᵢ = Y[2]      # Ice mixing ratio
+    qᵣ = Y[3]      # Rime mass mixing ratio
+    Bᵣ = Y[4]      # Rime volume
+
+    F_r = qᵣ / qᵢ
+    ρ_r = ifelse(Bᵣ == 0, FT(0), qᵣ / Bᵣ)
+    q = TD.PhasePartition(qᵥ + qₗ + qᵢ, qₗ, qᵢ)
+
+    Rₐ = TD.gas_constant_air(tps, TD.PhasePartition(qᵥ))
+    R_v = TD.Parameters.R_v(tps)
+    e = qᵥ * pₐ * R_v / Rₐ
+
+    a_w = CMO.a_w_eT(tps, e, T)
+    a_w_ice = CMO.a_w_ice(tps, T)
+    Δa_w = a_w - a_w_ice
+    J_immersion = CMI_het.ABIFM_J(aero_type, Δa_w)
+    dNᵢ_dt = J_immersion * Nₗ * 4 * π * rₗ^2
+    println("J = ", J_immersion)
+
+    sol = P3.thresholds(p3, ρ_r, F_r)
+    println(" ")
+    println("D_cr = ", sol.D_cr)
+    println("D_gr = ", sol.D_gr)
+    println("ρ_g =  ", sol.ρ_g)
+    println("ρ_dr = ", sol.ρ_d)
+
+    # TODO - compute the corresponding N0 and λ
+    # TODO - add ice nucleation source terms
+    # TODO - allow for zero rimed mass and volume
+
+    # Set tendencies
+    dY[1] = dNᵢ_dt       # ice number concentration
+    dY[2] = FT(0)        # ice mixing ratio
+    dY[3] = FT(0)        # rime mixing ratio
+    dY[4] = FT(0)        # rime volume
+end
+
+FT = Float64
+
+const_dt = 1.0
+t_0 = 0
+t_end = 2
+
+tps = TD.Parameters.ThermodynamicsParameters(FT) # thermodynamics free parameters
+wps = CMP.WaterProperties(FT)
+
+T = FT(251)
+pₐ = FT(800 * 1e2)   # air pressure
+ρₗ = wps.ρw
+Nₗ = FT(500 * 1e3)   # number of cloud droplets
+rₗ = FT(1e-6)        # radius of dropletes
+qᵥ = FT(8.1e-4)      # mixing ratio of water vapor
+qₗ = Nₗ * 4 / 3 * π * rₗ^3 * ρₗ / 1.2 # 1.2 should be ρₐ
+aero_type = CMP.Illite(FT)
+
+p = (; const_dt, tps, T, pₐ, qᵥ, qₗ, Nₗ, rₗ, aero_type)
+
+Nᵢ_0 = 100 * 1e6
+qᵢ_0 = 1e-3
+qᵣ_0 = 1e-4
+Bᵣ_0 = 1 / 200 * 1e-4
+
+IC = [FT(Nᵢ_0), FT(qᵢ_0), FT(qᵣ_0), FT(Bᵣ_0)]
+
+problem = ODE.ODEProblem(p3_sandbox, IC, (FT(t_0), FT(t_end)), p)
+sol = ODE.solve(
+    problem,
+    ODE.Euler(),
+    dt = const_dt,
+    reltol = 10 * eps(FT),
+    abstol = 10 * eps(FT),
+)
+
+println("number of ice crystals = ", sol[1, :])

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)
@@ -350,8 +350,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