diff --git a/p3_sandbox/Project.toml b/p3_sandbox/Project.toml new file mode 100644 index 0000000000..5c6f2144f5 --- /dev/null +++ b/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" diff --git a/p3_sandbox/p3_sandbox.jl b/p3_sandbox/p3_sandbox.jl new file mode 100644 index 0000000000..b1e4a3a30d --- /dev/null +++ b/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, :]) diff --git a/parcel/Immersion_Freezing.jl b/parcel/Immersion_Freezing.jl index 3c2630a5e4..2f279a812f 100644 --- a/parcel/Immersion_Freezing.jl +++ b/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) diff --git a/parcel/Jensen_et_al_2022.jl b/parcel/Jensen_et_al_2022.jl index e6c728dca7..e09ee88c19 100644 --- a/parcel/Jensen_et_al_2022.jl +++ b/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) diff --git a/parcel/Liquid_only.jl b/parcel/Liquid_only.jl index 7f6480d42b..1377752f8a 100644 --- a/parcel/Liquid_only.jl +++ b/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) diff --git a/parcel/Tully_et_al_2023.jl b/parcel/Tully_et_al_2023.jl index eba2d70491..73df491f02 100644 --- a/parcel/Tully_et_al_2023.jl +++ b/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) diff --git a/parcel/parcel.jl b/parcel/parcel.jl index afd017719e..b85a8bb17c 100644 --- a/parcel/parcel.jl +++ b/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 diff --git a/src/AerosolActivation.jl b/src/AerosolActivation.jl index 482a7bb783..71ac1ae374 100644 --- a/src/AerosolActivation.jl +++ b/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) diff --git a/src/parameters/AerosolActivation.jl b/src/parameters/AerosolActivation.jl index 889891f865..f460650af5 100644 --- a/src/parameters/AerosolActivation.jl +++ b/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) diff --git a/src/parameters/toml/ARG2000.toml b/src/parameters/toml/ARG2000.toml new file mode 100644 index 0000000000..011884bc9a --- /dev/null +++ b/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 diff --git a/test/aerosol_activation_tests.jl b/test/aerosol_activation_tests.jl index 0e2534afe7..cbb5b3fb56 100644 --- a/test/aerosol_activation_tests.jl +++ b/test/aerosol_activation_tests.jl @@ -1,6 +1,6 @@ import Test as TT -import CLIMAParameters +import CLIMAParameters as CP import CloudMicrophysics as CM import Thermodynamics as TD @@ -14,7 +14,13 @@ function test_aerosol_activation(FT) tps = TD.Parameters.ThermodynamicsParameters(FT) aip = CMP.AirProperties(FT) - ap = CMP.AerosolActivationParameters(FT) + #default ARG2000 parameter values + ap_default = CMP.AerosolActivationParameters(FT) + # calibrated values based on PySDM + override_file = + joinpath(pkgdir(CM), "src", "parameters", "toml", "ARG2000.toml") + toml_dict = CP.create_toml_dict(FT; override_file) + ap_calibrated = CMP.AerosolActivationParameters(toml_dict) # Atmospheric conditions T = FT(294) # air temperature K @@ -132,124 +138,142 @@ function test_aerosol_activation(FT) TT.@testset "callable and positive" begin for AM_t in (AM_3_B, AM_3_κ) - TT.@test all(AA.mean_hygroscopicity_parameter(ap, AM_t) .> 0.0) - TT.@test AA.max_supersaturation(ap, AM_t, aip, tps, T, p, w, q) > - 0.0 - TT.@test all( - AA.N_activated_per_mode(ap, AM_t, aip, tps, T, p, w, q) .> 0.0, - ) - TT.@test all( - AA.M_activated_per_mode(ap, AM_t, aip, tps, T, p, w, q) .> 0.0, - ) - TT.@test AA.total_N_activated(ap, AM_t, aip, tps, T, p, w, q) > 0.0 - TT.@test AA.total_M_activated(ap, AM_t, aip, tps, T, p, w, q) > 0.0 + for ap in (ap_default, ap_calibrated) + TT.@test all(AA.mean_hygroscopicity_parameter(ap, AM_t) .> 0.0) + TT.@test AA.max_supersaturation( + ap, + AM_t, + aip, + tps, + T, + p, + w, + q, + ) > 0.0 + TT.@test all( + AA.N_activated_per_mode(ap, AM_t, aip, tps, T, p, w, q) .> + 0.0, + ) + TT.@test all( + AA.M_activated_per_mode(ap, AM_t, aip, tps, T, p, w, q) .> + 0.0, + ) + TT.@test AA.total_N_activated(ap, AM_t, aip, tps, T, p, w, q) > + 0.0 + TT.@test AA.total_M_activated(ap, AM_t, aip, tps, T, p, w, q) > + 0.0 + end end end TT.@testset "same mean hygroscopicity for the same aerosol" begin + for ap in (ap_default, ap_calibrated) - TT.@test AA.mean_hygroscopicity_parameter(ap, AM_3_B)[1] == - AA.mean_hygroscopicity_parameter(ap, AM_1_B)[1] - - TT.@test AA.mean_hygroscopicity_parameter(ap, AM_3_B)[2] == - AA.mean_hygroscopicity_parameter(ap, AM_1_B)[1] + TT.@test AA.mean_hygroscopicity_parameter(ap, AM_3_B)[1] == + AA.mean_hygroscopicity_parameter(ap, AM_1_B)[1] - TT.@test AA.mean_hygroscopicity_parameter(ap, AM_3_κ)[1] == - AA.mean_hygroscopicity_parameter(ap, AM_1_κ)[1] + TT.@test AA.mean_hygroscopicity_parameter(ap, AM_3_B)[2] == + AA.mean_hygroscopicity_parameter(ap, AM_1_B)[1] - TT.@test AA.mean_hygroscopicity_parameter(ap, AM_3_κ)[2] == - AA.mean_hygroscopicity_parameter(ap, AM_1_κ)[1] + TT.@test AA.mean_hygroscopicity_parameter(ap, AM_3_κ)[1] == + AA.mean_hygroscopicity_parameter(ap, AM_1_κ)[1] + TT.@test AA.mean_hygroscopicity_parameter(ap, AM_3_κ)[2] == + AA.mean_hygroscopicity_parameter(ap, AM_1_κ)[1] + end end TT.@testset "B and kappa hygroscopicities are equivalent" begin - - TT.@test all( - isapprox( - AA.mean_hygroscopicity_parameter(ap, AM_3_κ)[2], - AA.mean_hygroscopicity_parameter(ap, AM_3_B)[2], - rtol = 0.1, - ), - ) + for ap in (ap_default, ap_calibrated) + TT.@test all( + isapprox( + AA.mean_hygroscopicity_parameter(ap, AM_3_κ)[2], + AA.mean_hygroscopicity_parameter(ap, AM_3_B)[2], + rtol = 0.1, + ), + ) + end end TT.@testset "order of modes does not matter" begin - - TT.@test AA.total_N_activated(ap, AM_3_B, aip, tps, T, p, w, q) == - AA.total_N_activated(ap, AM_2_B, aip, tps, T, p, w, q) - TT.@test AA.total_M_activated(ap, AM_3_B, aip, tps, T, p, w, q) == - AA.total_M_activated(ap, AM_2_B, aip, tps, T, p, w, q) - - TT.@test AA.total_N_activated(ap, AM_3_κ, aip, tps, T, p, w, q) == - AA.total_N_activated(ap, AM_2_κ, aip, tps, T, p, w, q) - TT.@test AA.total_M_activated(ap, AM_3_κ, aip, tps, T, p, w, q) == - AA.total_M_activated(ap, AM_2_κ, aip, tps, T, p, w, q) + for ap in (ap_default, ap_calibrated) + TT.@test AA.total_N_activated(ap, AM_3_B, aip, tps, T, p, w, q) == + AA.total_N_activated(ap, AM_2_B, aip, tps, T, p, w, q) + TT.@test AA.total_M_activated(ap, AM_3_B, aip, tps, T, p, w, q) == + AA.total_M_activated(ap, AM_2_B, aip, tps, T, p, w, q) + + TT.@test AA.total_N_activated(ap, AM_3_κ, aip, tps, T, p, w, q) == + AA.total_N_activated(ap, AM_2_κ, aip, tps, T, p, w, q) + TT.@test AA.total_M_activated(ap, AM_3_κ, aip, tps, T, p, w, q) == + AA.total_M_activated(ap, AM_2_κ, aip, tps, T, p, w, q) + end end TT.@testset "Abdul-Razzak and Ghan 2000 Fig 1" begin - - # data read from Fig 1 in Abdul-Razzak and Ghan 2000 - # using https://automeris.io/WebPlotDigitizer/ - N_2_obs = [ - FT(18.74716810149539), - FT(110.41572270049846), - FT(416.00589034889026), - FT(918.1014952424102), - FT(1914.816492976891), - FT(4919.913910285455), - ] - N_act_obs = [ - FT(0.7926937018577255), - FT(0.7161078386950611), - FT(0.5953670140462167), - FT(0.4850589034888989), - FT(0.34446080652469424), - FT(0.162630267331219), - ] - - N_act_frac_B = Vector{FT}(undef, 6) - N_act_frac_κ = Vector{FT}(undef, 6) - - it = 1 - for N_2_paper in N_2_obs - paper_mode_2_B = AM.Mode_B( - r_dry_paper, - stdev_paper, - N_2_paper * FT(1e6), - (FT(1.0),), - (sulfate.ϵ,), - (sulfate.ϕ,), - (sulfate.M,), - (sulfate.ν,), - (sulfate.ρ,), - 1, - ) - paper_mode_2_κ = AM.Mode_κ( - r_dry_paper, - stdev_paper, - N_2_paper * FT(1e6), - (FT(1.0),), - (FT(1.0),), - (sulfate.M,), - (sulfate.κ,), - 1, - ) - - AD_B = AM.AerosolDistribution((paper_mode_1_B, paper_mode_2_B)) - AD_κ = AM.AerosolDistribution((paper_mode_1_κ, paper_mode_2_κ)) - - N_act_frac_B[it] = - AA.N_activated_per_mode(ap, AD_B, aip, tps, T, p, w, q)[1] / - N_1_paper - N_act_frac_κ[it] = - AA.N_activated_per_mode(ap, AD_κ, aip, tps, T, p, w, q)[1] / - N_1_paper - it += 1 + for ap in (ap_default,) + + # data read from Fig 1 in Abdul-Razzak and Ghan 2000 + # using https://automeris.io/WebPlotDigitizer/ + N_2_obs = [ + FT(18.74716810149539), + FT(110.41572270049846), + FT(416.00589034889026), + FT(918.1014952424102), + FT(1914.816492976891), + FT(4919.913910285455), + ] + N_act_obs = [ + FT(0.7926937018577255), + FT(0.7161078386950611), + FT(0.5953670140462167), + FT(0.4850589034888989), + FT(0.34446080652469424), + FT(0.162630267331219), + ] + + N_act_frac_B = Vector{FT}(undef, 6) + N_act_frac_κ = Vector{FT}(undef, 6) + + it = 1 + for N_2_paper in N_2_obs + paper_mode_2_B = AM.Mode_B( + r_dry_paper, + stdev_paper, + N_2_paper * FT(1e6), + (FT(1.0),), + (sulfate.ϵ,), + (sulfate.ϕ,), + (sulfate.M,), + (sulfate.ν,), + (sulfate.ρ,), + 1, + ) + paper_mode_2_κ = AM.Mode_κ( + r_dry_paper, + stdev_paper, + N_2_paper * FT(1e6), + (FT(1.0),), + (FT(1.0),), + (sulfate.M,), + (sulfate.κ,), + 1, + ) + + AD_B = AM.AerosolDistribution((paper_mode_1_B, paper_mode_2_B)) + AD_κ = AM.AerosolDistribution((paper_mode_1_κ, paper_mode_2_κ)) + + N_act_frac_B[it] = + AA.N_activated_per_mode(ap, AD_B, aip, tps, T, p, w, q)[1] / + N_1_paper + N_act_frac_κ[it] = + AA.N_activated_per_mode(ap, AD_κ, aip, tps, T, p, w, q)[1] / + N_1_paper + it += 1 + end + + TT.@test all(isapprox(N_act_frac_B, N_act_obs, rtol = 0.05)) + TT.@test all(isapprox(N_act_frac_κ, N_act_obs, rtol = 0.1)) end - - TT.@test all(isapprox(N_act_frac_B, N_act_obs, rtol = 0.05)) - TT.@test all(isapprox(N_act_frac_κ, N_act_obs, rtol = 0.1)) - end end