From 9713f43330be11cd55656b7fb9a56918a4503eec Mon Sep 17 00:00:00 2001 From: Anna Jaruga Date: Thu, 20 Jul 2023 04:07:48 -0700 Subject: [PATCH] Separate parcel from Tully example, update docs to final paper --- docs/bibliography.bib | 15 +-- docs/src/IceNucleationBox.md | 6 +- parcel/Tully_et_al_2023.jl | 191 +++++++++++++++++++++++++++++++++++ parcel/parcel.jl | 185 +-------------------------------- 4 files changed, 203 insertions(+), 194 deletions(-) create mode 100644 parcel/Tully_et_al_2023.jl diff --git a/docs/bibliography.bib b/docs/bibliography.bib index 57f96124c..629e62a73 100644 --- a/docs/bibliography.bib +++ b/docs/bibliography.bib @@ -439,14 +439,15 @@ @article{TripoliCotton1980 doi = {10.1175/1520-0450(1980)019<1037:ANIOSF>2.0.CO;2} } -@article{Tully2022, +@article{Tully2023, + title = {Assessing predicted cirrus ice properties between two deterministic ice formation parameterizations}, author = {Tully, C. and Neubauer, D. and Lohmann, U.}, - title = {Technical Note: assessing predicted cirrus ice properties between two deterministic ice formation parameterizations}, - journal = {EGUsphere}, - volume = {2022}, - year = {2022}, - pages = {1--25}, - doi = {10.5194/egusphere-2022-1057} + journal = {Geoscientific Model Development}, + volume = {16}, + year = {2023}, + number = {10}, + pages = {2957--2973}, + doi = {10.5194/gmd-16-2957-2023} } @article{Vehkamaki2002, diff --git a/docs/src/IceNucleationBox.md b/docs/src/IceNucleationBox.md index ddad8a73b..206de11c0 100644 --- a/docs/src/IceNucleationBox.md +++ b/docs/src/IceNucleationBox.md @@ -4,7 +4,7 @@ This is a 0-dimensional adiabatic parcel model for testing nucleation schemes. It is based on [KorolevMazin2003](@cite), as well as the cirrus box model - [Karcher2006](@cite), [Tully2022](@cite). + [Karcher2006](@cite), [Tully2023](@cite). The model solves for saturation in a 0-dimensional adiabatic parcel raising with constant velocity. @@ -172,6 +172,6 @@ Between each run the water vapor specific humidity is changed, The prescribed vertical velocity is equal to 3.5 cm/s. ```@example -include("../../parcel/parcel.jl") +include("../../parcel/Tully_et_al_2023.jl") ``` -![](cirrus_box.svg) \ No newline at end of file +![](cirrus_box.svg) diff --git a/parcel/Tully_et_al_2023.jl b/parcel/Tully_et_al_2023.jl new file mode 100644 index 000000000..eb7751a2d --- /dev/null +++ b/parcel/Tully_et_al_2023.jl @@ -0,0 +1,191 @@ +import OrdinaryDiffEq as ODE +import CairoMakie as MK + +import Thermodynamics as TD +import CloudMicrophysics as CM +import CLIMAParameters as CP + +# boilerplate code to get free parameter values +include(joinpath(pkgdir(CM), "test", "create_parameters.jl")) + +# definition of the ODE problem for parcel model +include(joinpath(pkgdir(CM), "parcel", "parcel.jl")) + +""" + Wrapper for initial condition +""" +function get_initial_condition( + prs, + N_act, + p_a, + T, + q_vap, + q_liq, + q_ice, + N_aerosol, +) + thermo_params = CMP.thermodynamics_params(prs) + q = TD.PhasePartition(q_vap + q_liq + q_ice, q_liq, q_ice) + R_a = TD.gas_constant_air(thermo_params, q) + R_v = CMP.R_v(prs) + e_si = TD.saturation_vapor_pressure(thermo_params, T, TD.Ice()) + e = q_vap * p_a * R_v / R_a + S_i = e / e_si + + return [S_i, N_act, p_a, T, q_vap, q_ice, N_aerosol] +end + +""" + Wrapper for running the simulation following the same framework as in + Tully et al 2023 (https://doi.org/10.5194/gmd-16-2957-2023) + - the simulation consists of 3 periods mimicking 3 large scale model steps + - each period is 30 minutes long + - each period is run with user specified constant timestep + - the large scale initial conditions between each period are different +""" +function run_parcel(FT) + + # Boiler plate code to have access to model parameters and constants + toml_dict = CP.create_toml_dict(FT; dict_type = "alias") + prs = cloud_microphysics_parameters(toml_dict) + thermo_params = CMP.thermodynamics_params(prs) + + # Initial conditions for 1st period + N_aerosol = FT(2000 * 1e3) + N_0 = FT(0) + p_0 = FT(20000) + T_0 = FT(230) + q_vap_0 = FT(0.0003345) + q_liq_0 = FT(0) + q_ice_0 = FT(0) + # Initial conditions for the 2nd period + T2 = FT(229.25) + q_vap2 = FT(0.00034) + # Initial conditions for the 3rd period + T3 = FT(228.55) + q_vap3 = FT(0.000345) + + # Simulation time + t_max = 30 * 60 + + # Simulation parameters passed into ODE solver + r_nuc = FT(0.5 * 1.e-4 * 1e-6) # assumed size of nucleated particles + w = FT(3.5 * 1e-2) # updraft speed + α_m = FT(0.5) # accomodation coefficient + const_dt = 0.1 # model timestep + p = (; prs, const_dt, r_nuc, w, α_m) + + # Simulation 1 + IC1 = get_initial_condition( + prs, + N_0, + p_0, + T_0, + q_vap_0, + q_liq_0, + q_ice_0, + N_aerosol, + ) + prob1 = ODE.ODEProblem(parcel_model, IC1, (FT(0), t_max), p) + sol1 = ODE.solve( + prob1, + ODE.Euler(), + dt = const_dt, + reltol = 10 * eps(FT), + abstol = 10 * eps(FT), + ) + + # Simulation 2 + # (alternatively set T and take q_vap from the previous simulation) + #IC2 = get_initial_condition(sol1[2, end], sol1[3, end], T2, sol1[5, end], 0.0, sol1[6, end], sol1[7, end]) + IC2 = get_initial_condition( + prs, + sol1[2, end], + sol1[3, end], + sol1[4, end], + q_vap2, + q_liq_0, + sol1[6, end], + sol1[7, end], + ) + prob2 = ODE.ODEProblem( + parcel_model, + IC2, + (FT(sol1.t[end]), sol1.t[end] + t_max), + p, + ) + sol2 = ODE.solve( + prob2, + ODE.Euler(), + dt = const_dt, + reltol = 10 * eps(FT), + abstol = 10 * eps(FT), + ) + + # Simulation 3 + # (alternatively set T and take q_vap from the previous simulation) + #IC3 = get_initial_condition(sol2[2, end], sol2[3, end], T3, sol2[5, end], 0.0, sol2[6, end], sol2[7, end]) + IC3 = get_initial_condition( + prs, + sol2[2, end], + sol2[3, end], + sol2[4, end], + q_vap3, + q_liq_0, + sol2[6, end], + sol2[7, end], + ) + prob3 = ODE.ODEProblem( + parcel_model, + IC3, + (FT(sol2.t[end]), sol2.t[end] + t_max), + p, + ) + sol3 = ODE.solve( + prob3, + ODE.Euler(), + dt = const_dt, + reltol = 10 * eps(FT), + abstol = 10 * eps(FT), + ) + + # Plot results + fig = MK.Figure(resolution = (800, 600)) + ax1 = MK.Axis(fig[1, 1], ylabel = "Supersaturation [-]") + ax2 = MK.Axis(fig[1, 2], ylabel = "Temperature [K]") + ax3 = MK.Axis(fig[2, 1], ylabel = "N act [1/dm3]", yscale = log10) + ax4 = MK.Axis(fig[2, 2], ylabel = "N areo [1/dm3]") + ax5 = MK.Axis(fig[3, 1], ylabel = "q_vap [g/kg]", xlabel = "Height [m]") + ax6 = MK.Axis(fig[3, 2], ylabel = "q_ice [g/kg]", xlabel = "Height [m]") + + MK.ylims!(ax1, 1.0, 1.5) + MK.ylims!(ax3, 3, 2e3) + + MK.lines!(ax1, sol1.t * w, sol1[1, :]) + MK.lines!(ax1, sol2.t * w, sol2[1, :]) + MK.lines!(ax1, sol3.t * w, sol3[1, :]) + + MK.lines!(ax2, sol1.t * w, sol1[4, :]) + MK.lines!(ax2, sol2.t * w, sol2[4, :]) + MK.lines!(ax2, sol3.t * w, sol3[4, :]) + + MK.lines!(ax3, sol1.t * w, sol1[2, :] * 1e-3) + MK.lines!(ax3, sol2.t * w, sol2[2, :] * 1e-3) + MK.lines!(ax3, sol3.t * w, sol3[2, :] * 1e-3) + + MK.lines!(ax4, sol1.t * w, sol1[7, :] * 1e-3) + MK.lines!(ax4, sol2.t * w, sol2[7, :] * 1e-3) + MK.lines!(ax4, sol3.t * w, sol3[7, :] * 1e-3) + + MK.lines!(ax5, sol1.t * w, sol1[5, :] * 1e3) + MK.lines!(ax5, sol2.t * w, sol2[5, :] * 1e3) + MK.lines!(ax5, sol3.t * w, sol3[5, :] * 1e3) + + MK.lines!(ax6, sol1.t * w, sol1[6, :] * 1e3) + MK.lines!(ax6, sol2.t * w, sol2[6, :] * 1e3) + MK.lines!(ax6, sol3.t * w, sol3[6, :] * 1e3) + + MK.save("cirrus_box.svg", fig) +end + +run_parcel(Float64) diff --git a/parcel/parcel.jl b/parcel/parcel.jl index c4995a99d..6c91bbeca 100644 --- a/parcel/parcel.jl +++ b/parcel/parcel.jl @@ -1,9 +1,5 @@ -import OrdinaryDiffEq as ODE -import CairoMakie as MK - import Thermodynamics as TD import CloudMicrophysics as CM -import CLIMAParameters as CP const CMT = CM.CommonTypes const CMO = CM.Common @@ -15,7 +11,7 @@ include(joinpath(pkgdir(CM), "test", "create_parameters.jl")) """ ODE problem definition """ -function cirrus_box(dY, Y, p, t) +function parcel_model(dY, Y, p, t) # Get simulation parameters (; prs, const_dt, r_nuc, w, α_m) = p @@ -99,182 +95,3 @@ function cirrus_box(dY, Y, p, t) # TODO - add diagnostics output (radius, S, etc) end - -""" - Wrapper for initial condition -""" -function get_initial_condition( - prs, - N_act, - p_a, - T, - q_vap, - q_liq, - q_ice, - N_aerosol, -) - thermo_params = CMP.thermodynamics_params(prs) - q = TD.PhasePartition(q_vap + q_liq + q_ice, q_liq, q_ice) - R_a = TD.gas_constant_air(thermo_params, q) - R_v = CMP.R_v(prs) - e_si = TD.saturation_vapor_pressure(thermo_params, T, TD.Ice()) - e = q_vap * p_a * R_v / R_a - S_i = e / e_si - - return [S_i, N_act, p_a, T, q_vap, q_ice, N_aerosol] -end - -""" - Wrapper for running the simulation following the same framework as in - Tully et al 2022 (10.5194/egusphere-2022-1057) - - the simulation consists of 3 periods mimicking 3 large scale model steps - - each period is 30 minutes long - - each period is run with user specified constant timestep - - the large scale initial conditions between each period are different -""" -function run_parcel(FT) - - # Boiler plate code to have access to model parameters and constants - toml_dict = CP.create_toml_dict(FT; dict_type = "alias") - prs = cloud_microphysics_parameters(toml_dict) - thermo_params = CMP.thermodynamics_params(prs) - - # Initial conditions for 1st period - N_aerosol = FT(2000 * 1e3) - N_0 = FT(0) - p_0 = FT(20000) - T_0 = FT(230) - q_vap_0 = FT(0.0003345) - q_liq_0 = FT(0) - q_ice_0 = FT(0) - # Initial conditions for the 2nd period - T2 = FT(229.25) - q_vap2 = FT(0.00034) - # Initial conditions for the 3rd period - T3 = FT(228.55) - q_vap3 = FT(0.000345) - - # Simulation time - t_max = 30 * 60 - - # Simulation parameters passed into ODE solver - r_nuc = FT(0.5 * 1.e-4 * 1e-6) # assumed size of nucleated particles - w = FT(3.5 * 1e-2) # updraft speed - α_m = FT(0.5) # accomodation coefficient - const_dt = 0.1 # model timestep - p = (; prs, const_dt, r_nuc, w, α_m) - - # Simulation 1 - IC1 = get_initial_condition( - prs, - N_0, - p_0, - T_0, - q_vap_0, - q_liq_0, - q_ice_0, - N_aerosol, - ) - prob1 = ODE.ODEProblem(cirrus_box, IC1, (FT(0), t_max), p) - sol1 = ODE.solve( - prob1, - ODE.Euler(), - dt = const_dt, - reltol = 10 * eps(FT), - abstol = 10 * eps(FT), - ) - - # Simulation 2 - # (alternatively set T and take q_vap from the previous simulation) - #IC2 = get_initial_condition(sol1[2, end], sol1[3, end], T2, sol1[5, end], 0.0, sol1[6, end], sol1[7, end]) - IC2 = get_initial_condition( - prs, - sol1[2, end], - sol1[3, end], - sol1[4, end], - q_vap2, - q_liq_0, - sol1[6, end], - sol1[7, end], - ) - prob2 = ODE.ODEProblem( - cirrus_box, - IC2, - (FT(sol1.t[end]), sol1.t[end] + t_max), - p, - ) - sol2 = ODE.solve( - prob2, - ODE.Euler(), - dt = const_dt, - reltol = 10 * eps(FT), - abstol = 10 * eps(FT), - ) - - # Simulation 3 - # (alternatively set T and take q_vap from the previous simulation) - #IC3 = get_initial_condition(sol2[2, end], sol2[3, end], T3, sol2[5, end], 0.0, sol2[6, end], sol2[7, end]) - IC3 = get_initial_condition( - prs, - sol2[2, end], - sol2[3, end], - sol2[4, end], - q_vap3, - q_liq_0, - sol2[6, end], - sol2[7, end], - ) - prob3 = ODE.ODEProblem( - cirrus_box, - IC3, - (FT(sol2.t[end]), sol2.t[end] + t_max), - p, - ) - sol3 = ODE.solve( - prob3, - ODE.Euler(), - dt = const_dt, - reltol = 10 * eps(FT), - abstol = 10 * eps(FT), - ) - - # Plot results - fig = MK.Figure(resolution = (800, 600)) - ax1 = MK.Axis(fig[1, 1], ylabel = "Supersaturation [-]") - ax2 = MK.Axis(fig[1, 2], ylabel = "Temperature [K]") - ax3 = MK.Axis(fig[2, 1], ylabel = "N act [1/dm3]", yscale = log10) - ax4 = MK.Axis(fig[2, 2], ylabel = "N areo [1/dm3]") - ax5 = MK.Axis(fig[3, 1], ylabel = "q_vap [g/kg]", xlabel = "Height [m]") - ax6 = MK.Axis(fig[3, 2], ylabel = "q_ice [g/kg]", xlabel = "Height [m]") - - MK.ylims!(ax1, 1.0, 1.5) - MK.ylims!(ax3, 3, 2e3) - - MK.lines!(ax1, sol1.t * w, sol1[1, :]) - MK.lines!(ax1, sol2.t * w, sol2[1, :]) - MK.lines!(ax1, sol3.t * w, sol3[1, :]) - - MK.lines!(ax2, sol1.t * w, sol1[4, :]) - MK.lines!(ax2, sol2.t * w, sol2[4, :]) - MK.lines!(ax2, sol3.t * w, sol3[4, :]) - - MK.lines!(ax3, sol1.t * w, sol1[2, :] * 1e-3) - MK.lines!(ax3, sol2.t * w, sol2[2, :] * 1e-3) - MK.lines!(ax3, sol3.t * w, sol3[2, :] * 1e-3) - - MK.lines!(ax4, sol1.t * w, sol1[7, :] * 1e-3) - MK.lines!(ax4, sol2.t * w, sol2[7, :] * 1e-3) - MK.lines!(ax4, sol3.t * w, sol3[7, :] * 1e-3) - - MK.lines!(ax5, sol1.t * w, sol1[5, :] * 1e3) - MK.lines!(ax5, sol2.t * w, sol2[5, :] * 1e3) - MK.lines!(ax5, sol3.t * w, sol3[5, :] * 1e3) - - MK.lines!(ax6, sol1.t * w, sol1[6, :] * 1e3) - MK.lines!(ax6, sol2.t * w, sol2[6, :] * 1e3) - MK.lines!(ax6, sol3.t * w, sol3[6, :] * 1e3) - - MK.save("cirrus_box.svg", fig) -end - -run_parcel(Float64)