diff --git a/docs/src/P3Scheme.md b/docs/src/P3Scheme.md index 9647668bcd..b2864b5019 100644 --- a/docs/src/P3Scheme.md +++ b/docs/src/P3Scheme.md @@ -11,28 +11,38 @@ Traditionally, cloud ice microphysics schemes use various predefined a single ice category and evolving its properties. This simplification aids in attempts to constrain the scheme's free parameters. +The prognostic variables are: + - ``N_{ice}`` - number concentration 1/m3 + - ``q_{ice}`` - ice mass density kg/m3 + - ``q_{rim}`` - rime mass density kg/m3 + - ``B_{rim}`` - rime volume - TODO - units? + +!!! note + TODO - At some point we should switch to specific humidities... + ## Assumed particle size distribution (PSD) Following [MorrisonMilbrandt2015](@cite), the scheme assumes a gamma distribution for the concentration of particles per unit volume based on particle size measurements obtained by [Heymsfield2003](@cite) - in tropical and midlatitude ice clouds and implemented by + in tropical and midlatitude ice clouds and implemented by [MorrisonGrabowski2008](@cite): ```math N'(D) = N_{0} D^\mu \, e^{-\lambda \, D} ``` - where: - ``N'`` is the number concentration (``m^{-4}``), - ``D`` is the maximum particle dimension (``m``), - ``N_0`` is the intercept parameter (``m^{-4}``), - - ``\mu`` is the shape parameter (unitless I believe), + - ``\mu`` is the shape parameter (dimensionless), - ``\lambda`` is the slope parameter (``m^{-1}``). -We assume ``\mu \ = 0.00191 \lambda \ ^{0.8} - 2`` for ``\mu \ \in (0,6)``, which occurs only for ``\frac{1}{\gamma} < ~0.17 mm``. - -``N_0`` and ``\lambda`` are found using different moments of the PSD, total, cumulative number concentration ``N`` and mass mixing ratio ``q``, where +We assume ``\mu \ = 0.00191 \; \lambda \ ^{0.8} - 2``. +Following [MorrisonGrabowski2008](@cite) we limit ``\mu \ \in (0,6)``. +A non-zero ``\mu`` can occur only for very small mean particle sizes``\frac{1}{\lambda} < ~0.17 mm``. +``N_0`` and ``\lambda`` can be found using different moments of the PSD, +namely the total number concentration ``N`` and mass mixing ratio ``q``, where ```math N = \int_{0}^{\infty} \! N'(D) \mathrm{d}D @@ -42,44 +52,47 @@ N = \int_{0}^{\infty} \! N'(D) \mathrm{d}D q = \int_{0}^{\infty} \! m(D) N'(D) \mathrm{d}D ``` +!!! note + TODO - write the formulae for ``N_0`` and ``\lambda`` + ## Assumed particle mass relationships The mass ``m`` of particles as a function of maximum particle dimension ``D`` - is a piecewise function with variable breakpoints described + is a piecewise function with variable thresholds described by the following table. -| particle properties | condition(s) | m(D) relation | -|:--------------------|:----------------------|:----------------------| -|small, spherical ice | ``D < D_{th}`` | ``\frac{\pi}{6} \rho_i \ D^3`` | -|large, unrimed ice | ``q_{rim} = 0`` and ``D > D_{th}`` | ``\alpha_{va} \ D^{\beta_{va}}`` | -|dense nonspherical ice | ``q_{rim} > 0`` and ``D_{th} < D < D_{gr}`` | ``\alpha_{va} \ D^{\beta_{va}}`` | -|partially rimed ice | ``q_{rim} > 0`` and ``D > D_{cr}`` | ``\frac{\alpha_{va}}{1-F_r} D^{\beta_{va}}`` | -|graupel (completely rimed, spherical)| ``q_{rim} > 0``and ``D_{gr} < D < D_{cr}`` | ``\frac{\pi}{6} \rho_g \ D^3`` | - -with physical constant ``\rho_i \ = 916.7 kg m^{-3}``: see [CliMA Parameters](https://github.com/CliMA/CLIMAParameters.jl)'s `ρ_cloud_ice` ; -empirical coefficients: +| particle properties | condition(s) | m(D) relation | +|:-------------------------------------|:---------------------------------------------|:---------------------------------------------| +|small, spherical ice | ``D < D_{th}`` | ``\frac{\pi}{6} \rho_i \ D^3`` | +|large, unrimed ice | ``q_{rim} = 0`` and ``D > D_{th}`` | ``\alpha_{va} \ D^{\beta_{va}}`` | +|dense nonspherical ice | ``q_{rim} > 0`` and ``D_{th} < D < D_{gr}`` | ``\alpha_{va} \ D^{\beta_{va}}`` | +|graupel (completely rimed, spherical) | ``q_{rim} > 0`` and ``D_{gr} < D < D_{cr}`` | ``\frac{\pi}{6} \rho_g \ D^3`` | +|partially rimed ice | ``q_{rim} > 0`` and ``D > D_{cr}`` | ``\frac{\alpha_{va}}{1-F_r} D^{\beta_{va}}`` | - - ``\alpha_{va} = 7.38 \times 10^{-11} \times 10^{6 \beta_{va} - 3}`` (``kg m^{-β_va}``), a parameter from [BrownFrancis1995](@cite), derived from measurements of mass grown by vapor diffusion and aggregation in midlatitude cirrus, and modified to agree with CliMA's units of ``kg`` and ``m`` rather than [BrownFrancis1995](@cite)'s ``g`` and ``\mu m``, - - ``\beta_{va} = 1.9``, another parameter from [BrownFrancis1995](@cite); - -which together determine ``D_{th} = (\frac{\pi \rho_i}{6\alpha_{va}})^{\frac{1}{\beta_{va} - 3}}`` (``m``), the threshold particle dimension between small spherical and large, nonspherical unrimed ice. -Relevant model state variables here are: - - - ``q_{rim}``, rime mass concentration; - - ``q_{ice}``, total ice mass concentration; - - ``B_{rim}``, bulk rime volume; - -from which rime mass fraction ``F_r = \frac{q_rim}{q_ice}`` and predicted rime density ``\rho_{r} = \frac{q_rim}{B_rim}`` are derived. -The following four thresholds (``D_{i}``) and densities (``\rho_{i}``) which form a nonlinear system solved by `thresholds(ρ_r, F_r, u0)` , which employs [`NonlinearSolve.jl`](https://docs.sciml.ai/NonlinearSolve/stable/): - - - ``D_{gr} = (\frac{6\alpha_{va}}{\pi \rho_g})^{\frac{1}{3 - \beta_{va}}}`` (``m``), a threshold defined to ensure continuity and reasonable masses for smaller D; - - ``D_{cr} = [ (\frac{1}{1-F_r}) \frac{6 \alpha_{va}}{\pi \rho_g} ]^{\frac{1}{3 - \beta_{va}}}`` (``m``), the threshold separating partially rimed ice from graupel; - - ``\rho_g = \rho_r F_r + (1 - F_r) \rho_d`` (``kg m^{-3}``), bulk density of graupel, calculated with the rime mass fraction weighted average of the predicted rime density and the density of the unrimed part, where - - ``\rho_d = \frac{6\alpha_{va}(D_{cr}^{\beta{va} \ - 2} - D_{gr}^{\beta{va} \ - 2})}{\pi \ (\beta_{va} \ - 2)(D_{cr}-D_{gr})}`` (``kg m^{-3}``) is the density of the unrimed part. - -Below are graphical representations of the m(D) regime, replicating - Figures 1 (a) and (b) of [MorrisonMilbrandt2015](@cite): +where: + - ``D_{th}``, ``D_{gr}``, ``D_{cr}`` are particle size thresholds in ``m``, + - ``\rho_i`` is cloud ice density in ``kg m^{-3}``, + - ``\beta_{va} = 1.9`` is a dimensionless parameter from [BrownFrancis1995](@cite) (based on measurements of vapor diffusion and aggregation in midlatitude cirrus), + - ``\alpha_{va} = 7.38 \; 10^{-11} \; 10^{6 \beta_{va} - 3}`` in ``kg m^{-β_{va}}`` is a parameter from [BrownFrancis1995](@cite) in base SI units, + - ``\rho_g`` is the bulk density of graupel in ``kg m^{-3}``. + +The first threshold is solely determined by the free parameters: + ``D_{th} = (\frac{\pi \rho_i}{6\alpha_{va}})^{\frac{1}{\beta_{va} - 3}}``. +The remaining thresholds: ``D_{gr}``, ``D_{cr}``, as well as the + bulk density of graupel ``\rho_{g}``, + and the bulk density of the unrimed part ``\rho_d`` + form a nonlinear system: + - ``D_{gr} = (\frac{6\alpha_{va}}{\pi \rho_g})^{\frac{1}{3 - \beta_{va}}}`` + - ``D_{cr} = [ (\frac{1}{1-F_r}) \frac{6 \alpha_{va}}{\pi \rho_g} ]^{\frac{1}{3 - \beta_{va}}}`` + - ``\rho_g = \rho_r F_r + (1 - F_r) \rho_d`` + - ``\rho_d = \frac{6\alpha_{va}(D_{cr}^{\beta{va} \ - 2} - D_{gr}^{\beta{va} \ - 2})}{\pi \ (\beta_{va} \ - 2)(D_{cr}-D_{gr})}`` +where + - ``F_r = \frac{q_{rim}}{q_{ice}}`` is the rime mass fraction, + - ``\rho_{r} = \frac{q_{rim}}{B_{rim}}`` is the predicted rime density. +The system is solved using [`NonlinearSolve.jl`](https://docs.sciml.ai/NonlinearSolve/stable/). + +Below we show the m(D) regime, replicating Figures 1 (a) and (b) from [MorrisonMilbrandt2015](@cite). ```@example include("P3SchemePlots.jl") @@ -89,4 +102,4 @@ p3_m_plot2(["cyan2", "cyan4", "midnightblue"], ["hotpink", "magenta3", "purple4" ![](MorrisonandMilbrandtFig1a.svg) -![](MorrisonandMilbrandtFig1b.svg) \ No newline at end of file +![](MorrisonandMilbrandtFig1b.svg) diff --git a/docs/src/P3SchemePlots.jl b/docs/src/P3SchemePlots.jl index 9413933449..7d1a70ac66 100644 --- a/docs/src/P3SchemePlots.jl +++ b/docs/src/P3SchemePlots.jl @@ -13,6 +13,7 @@ thermo_params = CMP.thermodynamics_params(param_set) # bulk density of ice const ρ_i::FT = CMP.ρ_cloud_ice(param_set) +# TODO - use from CLIMAParameters const β_va::FT = 1.9 const α_va::FT = (7.38e-11) * 10^((6 * β_va) - 3) const D_th::FT = ((FT(π) * ρ_i) / (6 * α_va))^(1 / (β_va - 3)) @@ -106,6 +107,8 @@ function p3_m_plot1( ) D_range = range(1e-5, stop = 1e-2, length = len_D_range) + lw = 4 + fig1_a = Plt.Figure() ax1_a = Plt.Axis( @@ -128,6 +131,7 @@ function p3_m_plot1( D_range * 1e3, [mass(D, 0.0) for D in D_range], color = colors[1], + linewidth = lw, ) fig1_a_5 = Plt.lines!( @@ -135,6 +139,7 @@ function p3_m_plot1( D_range * 1e3, [mass(D, 0.5, P3.thresholds(400.0, 0.5)) for D in D_range], color = colors[2], + linewidth = lw, ) fig1_a_8 = Plt.lines!( @@ -142,6 +147,7 @@ function p3_m_plot1( D_range * 1e3, [mass(D, 0.8, P3.thresholds(400.0, 0.8)) for D in D_range], color = colors[3], + linewidth = lw, ) d_cr_5 = Plt.lines!( @@ -150,6 +156,7 @@ function p3_m_plot1( range(1e-14, stop = 1e-4, length = len_D_range), linestyle = "---", color = threshold_colors[2], + linewidth = lw, ) d_cr_8 = Plt.lines!( @@ -158,6 +165,7 @@ function p3_m_plot1( range(1e-14, stop = 1e-4, length = len_D_range), linestyle = "---", color = threshold_colors[3], + linewidth = lw, ) d_gr_5 = Plt.lines!( @@ -166,7 +174,7 @@ function p3_m_plot1( range(1e-14, stop = 1e-4, length = len_D_range), linestyle = "...", color = threshold_colors[2], - linewidth = 1.5, + linewidth = lw, ) d_gr_8 = Plt.lines!( @@ -175,7 +183,7 @@ function p3_m_plot1( range(1e-14, stop = 1e-4, length = len_D_range), linestyle = "...", color = threshold_colors[3], - linewidth = 1.5, + linewidth = lw, ) d_tha = Plt.lines!( @@ -184,6 +192,7 @@ function p3_m_plot1( range(1e-14, stop = 1e-4, length = len_D_range), linestyle = "---", color = "red", + linewidth = lw, ) leg1_a = Plt.Legend( @@ -227,6 +236,8 @@ function p3_m_plot2( fig1_b = Plt.Figure() + lw = 4 + ax1_b = Plt.Axis( fig1_b[1:10, 1:11], title = Plt.L"m(D) regime for $F_r = 0.95$", @@ -247,6 +258,7 @@ function p3_m_plot2( D_range * 1e3, [mass(D, 0.95, P3.thresholds(200.0, 0.95)) for D in D_range], color = colors[1], + linewidth = lw, ) fig1_b400 = Plt.lines!( @@ -254,6 +266,7 @@ function p3_m_plot2( D_range * 1e3, [mass(D, 0.95, P3.thresholds(400.0, 0.95)) for D in D_range], color = colors[2], + linewidth = lw, ) fig1_b800 = Plt.lines!( @@ -261,6 +274,7 @@ function p3_m_plot2( D_range * 1e3, [mass(D, 0.95, P3.thresholds(800.0, 0.95)) for D in D_range], color = colors[3], + linewidth = lw, ) d_thb = Plt.lines!( @@ -269,6 +283,7 @@ function p3_m_plot2( range(1e-14, stop = 1e-4, length = len_D_range), linestyle = "---", color = "red", + linewidth = lw, ) d_cr_200 = Plt.lines!( @@ -277,6 +292,7 @@ function p3_m_plot2( range(1e-14, stop = 1e-4, length = len_D_range), linestyle = "---", color = threshold_colors[1], + linewidth = lw, ) d_cr_400 = Plt.lines!( @@ -285,6 +301,7 @@ function p3_m_plot2( range(1e-14, stop = 1e-4, length = len_D_range), linestyle = "---", color = threshold_colors[2], + linewidth = lw, ) d_cr_800 = Plt.lines!( @@ -293,6 +310,7 @@ function p3_m_plot2( range(1e-14, stop = 1e-4, length = len_D_range), linestyle = "---", color = threshold_colors[3], + linewidth = lw, ) d_gr_200 = Plt.lines!( @@ -301,7 +319,7 @@ function p3_m_plot2( range(1e-14, stop = 1e-4, length = len_D_range), linestyle = "...", color = threshold_colors[1], - linewidth = 1.5, + linewidth = lw, ) d_gr_400 = Plt.lines!( @@ -310,6 +328,7 @@ function p3_m_plot2( range(1e-14, stop = 1e-4, length = len_D_range), linestyle = "...", color = threshold_colors[2], + linewidth = lw, ) d_gr_800 = Plt.lines!( @@ -318,7 +337,7 @@ function p3_m_plot2( range(1e-14, stop = 1e-4, length = len_D_range), linestyle = "...", color = threshold_colors[3], - linewidth = 1.5, + linewidth = lw, ) leg1_b = Plt.Legend( diff --git a/test/performance_tests.jl b/test/performance_tests.jl index 54c70eb70e..c9429b89ac 100644 --- a/test/performance_tests.jl +++ b/test/performance_tests.jl @@ -22,7 +22,7 @@ include(joinpath(pkgdir(CM), "test", "create_parameters.jl")) @info "Performance Tests" -function bench_press(foo, args, min_run_time) +function bench_press(foo, args, min_run_time, min_memory=0.0, min_allocs=0.0) println("Testing ", "$foo") # Calling foo once before benchmarking @@ -34,8 +34,8 @@ function bench_press(foo, args, min_run_time) println("\n") TT.@test BT.minimum(trail).time < min_run_time - TT.@test trail.memory <= 4e6 - TT.@test trail.allocs <= 4e4 + TT.@test trail.memory <= min_memory + TT.@test trail.allocs <= min_allocs end function benchmark_test(FT) @@ -90,7 +90,7 @@ function benchmark_test(FT) Delta_a_w = FT(0.27) # P3 scheme - bench_press(P3.thresholds, (ρ_r, F_r), 12e6) + bench_press(P3.thresholds, (ρ_r, F_r), 12e6, 4e6, 4e4) # aerosol activation bench_press(