Merge #148

148: Separate parcel from Tully example, update docs to final paper r=amylu00 a=trontrytel This PR separates the ODE definition of the parcel model from the example that uses it. This should help with using the parcel model to reproduce 0-dimensional simulations from different modelling studies. It also updates the documentation to the final published version of Tully et. al. 2023 Co-authored-by: Anna Jaruga <[email protected]>
CliMA · Jul 21, 2023 · 2fdd694 · 2fdd694
2 parents 91932bd + 9713f43
commit 2fdd694
Show file tree

Hide file tree

Showing 4 changed files with 203 additions and 194 deletions.
diff --git 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
@@ -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)
+![](cirrus_box.svg)
diff --git a/parcel/Tully_et_al_2023.jl 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)