diff --git a/docs/src/API.md b/docs/src/API.md index 6b7e209f92..3b4798b00c 100644 --- a/docs/src/API.md +++ b/docs/src/API.md @@ -93,6 +93,7 @@ HetIceNucleation.ABIFM_J HetIceNucleation.P3_deposition_N_i HetIceNucleation.P3_het_N_i HetIceNucleation.INP_concentration_frequency +HetIceNucleation.INP_concentration_mean ``` # Homogeneous ice nucleation diff --git a/parcel/Frostenberg_different_dt.jl b/parcel/Frostenberg_different_dt.jl new file mode 100644 index 0000000000..b43ec75452 --- /dev/null +++ b/parcel/Frostenberg_different_dt.jl @@ -0,0 +1,112 @@ +import OrdinaryDiffEq as ODE +import CairoMakie as MK +import Thermodynamics as TD +import CloudMicrophysics as CM +import CLIMAParameters as CP + +# definition of the ODE problem for parcel model +include(joinpath(pkgdir(CM), "parcel", "Parcel.jl")) +FT = Float32 +# get free parameters +tps = TD.Parameters.ThermodynamicsParameters(FT) +wps = CMP.WaterProperties(FT) + +# Initial conditions +ρₗ = wps.ρw +Nₐ = FT(0) +Nₗ = FT(500 * 1e3) +Nᵢ = FT(0) +r₀ = FT(1e-6) +p₀ = FT(800 * 1e2) +T₀ = FT(251) +qᵥ = FT(8.1e-4) +qₗ = Nₗ * 4 / 3 * FT(π) * r₀^3 * ρₗ / FT(1.2) # 1.2 should be ρₐ +qᵢ = FT(0) +x_sulph = FT(0.01) + +# Moisture dependent initial conditions +q = TD.PhasePartition.(qᵥ + qₗ + qᵢ, qₗ, qᵢ) +R_v = TD.Parameters.R_v(tps) +Rₐ = TD.gas_constant_air(tps, q) +eₛ = TD.saturation_vapor_pressure(tps, T₀, TD.Liquid()) +e = eᵥ(qᵥ, p₀, Rₐ, R_v) +Sₗ = FT(e / eₛ) +IC = [Sₗ, p₀, T₀, qᵥ, qₗ, qᵢ, Nₐ, Nₗ, Nᵢ, x_sulph] + +# Simulation parameters passed into ODE solver +w = FT(0.7) # updraft speed +t_max = FT(1200) +aerosol = CMP.Illite(FT) +heterogeneous = "Frostenberg_random" +condensation_growth = "Condensation" +deposition_growth = "Deposition" +DSD = "Monodisperse" +const_dt_range = range(FT(0.1), stop = FT(1), length = 4) #changing the time step and drawing in every time step + +# Plotting +fig = MK.Figure(resolution = (800, 600)) +ax1 = MK.Axis(fig[1, 1], ylabel = "Ice Supersaturation [-]") +ax2 = MK.Axis(fig[1, 2], ylabel = "Temperature [K]") +ax3 = MK.Axis(fig[2, 1], ylabel = "q_ice [g/kg]") +ax4 = MK.Axis(fig[2, 2], ylabel = "q_liq [g/kg]") +ax5 = MK.Axis(fig[3, 1], xlabel = "Time [min]", ylabel = "N_liq") +ax6 = MK.Axis(fig[3, 2], xlabel = "Time [min]", ylabel = "N_ice") + +for const_dt in const_dt_range + + params = parcel_params{FT}( + const_dt = const_dt, + w = w, + aerosol = aerosol, + heterogeneous = heterogeneous, + condensation_growth = condensation_growth, + deposition_growth = deposition_growth, + size_distribution = DSD, + drawing_interval = const_dt, #drawing in every time step + ) + # solve ODE + sol = run_parcel(IC, FT(0), t_max, params) + + # Plot results + MK.lines!( + ax1, + sol.t / 60, + S_i.(tps, sol[3, :], sol[1, :]) .- 1, + label = "dt=" * string(const_dt) * " s", + ) + MK.lines!( + ax2, + sol.t / 60, + sol[3, :], + label = "dt=" * string(const_dt) * " s", + ) + MK.lines!( + ax3, + sol.t / 60, + sol[6, :] * 1e3, + label = "dt=" * string(const_dt) * " s", + ) + MK.lines!( + ax4, + sol.t / 60, + sol[5, :] * 1e3, + label = "dt=" * string(const_dt) * " s", + ) + sol_Nₗ = sol[8, :] + sol_Nᵢ = sol[9, :] + MK.lines!(ax5, sol.t / 60, sol_Nₗ, label = "dt=" * string(const_dt) * " s") + MK.lines!(ax6, sol.t / 60, sol_Nᵢ, label = "dt=" * string(const_dt) * " s") +end + +for ax in [ax1, ax2, ax3, ax4, ax5, ax6] + MK.axislegend( + ax, + framevisible = false, + labelsize = 12, + orientation = :horizontal, + nbanks = 2, + position = :lc, + ) +end + +MK.save("dt.svg", fig) diff --git a/parcel/Frostenberg_different_frequency.jl b/parcel/Frostenberg_different_frequency.jl new file mode 100644 index 0000000000..29a74a9789 --- /dev/null +++ b/parcel/Frostenberg_different_frequency.jl @@ -0,0 +1,135 @@ +import OrdinaryDiffEq as ODE +import CairoMakie as MK +import Thermodynamics as TD +import CloudMicrophysics as CM +import CLIMAParameters as CP + +# definition of the ODE problem for parcel model +include(joinpath(pkgdir(CM), "parcel", "Parcel.jl")) +FT = Float32 +# get free parameters +tps = TD.Parameters.ThermodynamicsParameters(FT) +wps = CMP.WaterProperties(FT) + +# Initial conditions +ρₗ = wps.ρw +Nₐ = FT(0) +Nₗ = FT(500 * 1e3) +Nᵢ = FT(0) +r₀ = FT(1e-6) +p₀ = FT(800 * 1e2) +T₀ = FT(251) +qᵥ = FT(8.1e-4) +qₗ = Nₗ * 4 / 3 * FT(π) * r₀^3 * ρₗ / FT(1.2) # 1.2 should be ρₐ +qᵢ = FT(0) +x_sulph = FT(0.01) + +# Moisture dependent initial conditions +q = TD.PhasePartition.(qᵥ + qₗ + qᵢ, qₗ, qᵢ) +R_v = TD.Parameters.R_v(tps) +Rₐ = TD.gas_constant_air(tps, q) +eₛ = TD.saturation_vapor_pressure(tps, T₀, TD.Liquid()) +e = eᵥ(qᵥ, p₀, Rₐ, R_v) +Sₗ = FT(e / eₛ) +IC = [Sₗ, p₀, T₀, qᵥ, qₗ, qᵢ, Nₐ, Nₗ, Nᵢ, x_sulph] + +# Simulation parameters passed into ODE solver +w = FT(0.7) # updraft speed +t_max = FT(1200) +aerosol = CMP.Illite(FT) +heterogeneous = "Frostenberg_random" +condensation_growth = "Condensation" +deposition_growth = "Deposition" +DSD = "Monodisperse" +const_dt = FT(1) #constent time step +drawing_interval_range = range(FT(1), stop = FT(5), length = 4) # changing the time between drawing + +# Plotting +fig = MK.Figure(resolution = (800, 600)) +ax1 = MK.Axis(fig[1, 1], ylabel = "Ice Supersaturation [-]") +ax2 = MK.Axis(fig[1, 2], ylabel = "Temperature [K]") +ax3 = MK.Axis(fig[2, 1], ylabel = "q_ice [g/kg]") +ax4 = MK.Axis(fig[2, 2], ylabel = "q_liq [g/kg]") +ax5 = MK.Axis(fig[3, 1], xlabel = "Time [min]", ylabel = "N_liq") +ax6 = MK.Axis(fig[3, 2], xlabel = "Time [min]", ylabel = "N_ice") + +for drawing_interval in drawing_interval_range + + params = parcel_params{FT}( + const_dt = const_dt, + w = w, + aerosol = aerosol, + heterogeneous = heterogeneous, + condensation_growth = condensation_growth, + deposition_growth = deposition_growth, + size_distribution = DSD, + drawing_interval = drawing_interval, + ) + # solve ODE + sol = run_parcel(IC, FT(0), t_max, params) + + # Plot results + MK.lines!( + ax1, + sol.t / 60, + S_i.(tps, sol[3, :], sol[1, :]) .- 1, + label = "f =" * + string(round(const_dt / drawing_interval, digits = 1)) * + "/s", + ) + MK.lines!( + ax2, + sol.t / 60, + sol[3, :], + label = "f =" * + string(round(const_dt / drawing_interval, digits = 1)) * + "/s", + ) + MK.lines!( + ax3, + sol.t / 60, + sol[6, :] * 1e3, + label = "f =" * + string(round(const_dt / drawing_interval, digits = 1)) * + "/s", + ) + MK.lines!( + ax4, + sol.t / 60, + sol[5, :] * 1e3, + label = "f =" * + string(round(const_dt / drawing_interval, digits = 1)) * + "/s", + ) + sol_Nₗ = sol[8, :] + sol_Nᵢ = sol[9, :] + MK.lines!( + ax5, + sol.t / 60, + sol_Nₗ, + label = "f =" * + string(round(const_dt / drawing_interval, digits = 1)) * + "/s", + ) + MK.lines!( + ax6, + sol.t / 60, + sol_Nᵢ, + label = "f =" * + string(round(const_dt / drawing_interval, digits = 1)) * + "/s", + ) +end + +for ax in [ax1, ax2, ax3, ax4, ax5, ax6] + MK.axislegend( + ax, + framevisible = false, + labelsize = 12, + orientation = :horizontal, + nbanks = 2, + position = :lc, + ) +end + +MK.save("frequency.svg", fig)