diff --git a/NEWS.md b/NEWS.md deleted file mode 100644 index 4f2a9d10e..000000000 --- a/NEWS.md +++ /dev/null @@ -1,31 +0,0 @@ -CloudMicrophysics.jl Release Notes -======================== - -main ------- - - -- Add aerosol activation to parcel. ([#429](https://github.com/CliMA/CloudMicrophysics.jl/pull/429)) - -- Bugfix; fixed the evaporation scheme of SB2006 to prevent returning NaN values when limiters are not applied. ([#420](https://github.com/CliMA/CloudMicrophysics.jl/pull/415)) - -v0.20.0 ------- -- Added modified terminal velocity without limiters for SB2006 ([#415](https://github.com/CliMA/CloudMicrophysics.jl/pull/415)) - -- Add lognormal droplet distribution in parcel model ([#414](https://github.com/CliMA/CloudMicrophysics.jl/pull/414)) - -- Added calling function to grab artifacts ([#413](https://github.com/CliMA/CloudMicrophysics.jl/pull/413)) - -- API change: - - return a named tuple in 2-moment microphysics rain evaporation - - changed aerosol size distribution constructor - - [#392](https://github.com/CliMA/CloudMicrophysics.jl/pull/392) - -- Added AIDA homogeneous ice nucleation data as artifacts ([#388](https://github.com/CliMA/CloudMicrophysics.jl/pull/388)) - -- Generalize calibration functions in ice_nucleation_2024 ([#380](https://github.com/CliMA/CloudMicrophysics.jl/pull/380)) - -v0.19.0 ------- -- Started changelog ([#383](https://github.com/CliMA/CloudMicrophysics.jl/pull/383)) diff --git a/Project.toml b/Project.toml index c2b0aa48d..2dd6d733e 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "CloudMicrophysics" uuid = "6a9e3e04-43cd-43ba-94b9-e8782df3c71b" authors = ["Climate Modeling Alliance"] -version = "0.22.1" +version = "0.22.3" [deps] ClimaParams = "5c42b081-d73a-476f-9059-fd94b934656c" diff --git a/docs/bibliography.bib b/docs/bibliography.bib index c17a19c6f..1350e5299 100644 --- a/docs/bibliography.bib +++ b/docs/bibliography.bib @@ -87,6 +87,17 @@ @article{Bigg1953 doi = {10.1088/0370-1301/66/8/309} } +@article{Blahak2010, + title = {Efficient approximation of the incomplete gamma function for use in cloud model applications}, + author = {Blahak, U.}, + journal = {Geoscientific Model Development}, + volume = {3}, + year = {2010}, + number = {2}, + pages = {329--336}, + doi = {10.5194/gmd-3-329-2010} +} + @article{BrownFrancis1995, title = {Improved Measurements of the Ice Water Content in Cirrus Using a Total-Water Probe}, author = {Philip R. A. Brown and Peter N. Francis}, @@ -120,6 +131,16 @@ @article{China2017 year = {2017} } +@article{Choletteetal2019, + author = {Mélissa Cholette et al}, + title = {Parameterization of the Bulk Liquid Fraction on Mixed-Phase Particles in the Predicted Particle Properties (P3) Scheme: Description and Idealized Simulations}, + journal = {Journal of the Atmospheric Sciences}, + year = {2019}, + volume = {76}, + doi = {10.1175/JAS-D-18-0278.1}, + pages= {561--582} +} + @article{Desai2019, title = {Aerosol-Mediated Glaciation of Mixed-Phase Clouds: Steady-State Laboratory Measurements}, author = {Desai, Neel and Chandrakar, KK and Kinney, G and Cantrell, W and Shaw, RA}, @@ -616,6 +637,17 @@ @article{SeifertBeheng2006 publisher = {Springer} } +@article{Simmeletal2002, + title = {Numerical solution of the stochastic collection equation—comparison of the Linear Discrete Method with other methods}, + author = {Simmel, Martin and Trautmann, Thomas and Tetzlaff, Gerd}, + journal = {Atmoshperic Research}, + volume = {61}, + number = {2}, + year = {2002}, + doi = {10.1016/S0169-8095(01)00131-4}, + publisher = {Elsevier} +} + @article{Spichtinger2023, AUTHOR = {Spichtinger, P. and Marschalik, P. and Baumgartner, M.}, TITLE = {Impact of formulations of the homogeneous nucleation rate on ice nucleation events in cirrus}, diff --git a/docs/src/API.md b/docs/src/API.md index a47abd72f..8447fc1b2 100644 --- a/docs/src/API.md +++ b/docs/src/API.md @@ -10,6 +10,7 @@ MicrophysicsNonEq MicrophysicsNonEq.τ_relax MicrophysicsNonEq.conv_q_vap_to_q_liq_ice MicrophysicsNonEq.conv_q_vap_to_q_liq_ice_MM2015 +MicrophysicsNonEq.conv_q_vap_to_q_liq_ice_MM2015_timeintegrator ``` # 0-moment precipitation microphysics @@ -69,6 +70,7 @@ P3Scheme.thresholds P3Scheme.distribution_parameter_solver P3Scheme.ice_terminal_velocity P3Scheme.het_ice_nucleation +P3Scheme.ice_melt ``` # Aerosol model diff --git a/docs/src/MicrophysicsNonEq.md b/docs/src/MicrophysicsNonEq.md index 22603d973..c96fdd5c2 100644 --- a/docs/src/MicrophysicsNonEq.md +++ b/docs/src/MicrophysicsNonEq.md @@ -93,25 +93,145 @@ where: - ``R_v`` is the gas constant of water vapor, - ``L_v`` and ``L_s`` is the latent heat of vaporization and sublimation. - Note that these forms of condensation/sublimation and deposition/sublimation are equivalent to those described in the adiabatic parcel model with some rearrangements and assumptions. It is just necessary to use the definitions of ```\tau```, ```q_{sl}```, and the thermal diffusivity ```D_v```: - - ```math - \begin{equation} - \tau = 4 \pi N_{tot} \bar{r}, \;\;\;\;\;\;\;\; - q_{sl} = \frac{e_{sl}}{\rho R_v T}, \;\;\;\;\;\;\;\; - D_v = \frac{K}{\rho c_p} - \end{equation} - ``` - and if we assume that the supersaturation S can be approximated by specific humidities (this is only exactly true for mass mixing ratios): - ```math - \begin{equation} - S_l = \frac{q_{vap}}{q_{sl}} - \end{equation} - ``` - then we can write: - ```math - \begin{equation} - q_{vap} - q_{sl} = q_{sl}(S_l - 1) - \end{equation} - ``` - and ```Gamma``` is equivalent to the ```G``` function used in our parcel and parameterizations. \ No newline at end of file +Note that these forms of condensation/sublimation and deposition/sublimation are equivalent to those described in the adiabatic parcel model with some rearrangements and assumptions. It is just necessary to use the definitions of ```\tau```, ```q_{sat}```, and the thermal diffusivity ```D_v```: + +```math +\begin{equation} + \tau = 4 \pi N_{tot} \bar{r}, \;\;\;\;\;\;\;\; + q_{sat} = \frac{e_{sat}}{\rho R_v T}, \;\;\;\;\;\;\;\; + D_v = \frac{K}{\rho c_p} +\end{equation} +``` +and if we assume that the supersaturation S can be approximated by specific humidities (this is only exactly true for mass mixing ratios): +```math +\begin{equation} + S = \frac{q_{vap}}{q_{sat}} +\end{equation} +``` +then we can write: +```math +\begin{equation} + q_{vap} - q_{sat} = q_{sat}(S - 1) +\end{equation} +``` +and ```Gamma``` is equivalent to the ```G``` function used in our parcel and parameterizations. + + +### Time integrated version + +Finally, we can also use the time integrated formulation of condensation/evaporation and deposition/sublimation described by [MorrisonMilbrandt2015](@cite). We assume that the condensation/evaporation and deposition/sublimation over the course of the timestep can be calculated using the average difference between the $q_v$ value and those of saturation. Ie, they write: + +```math +\begin{equation} + \left( \frac{d q_l}{dt} \right) _{cond} = \frac{\bar{\delta_l}}{\tau_l \Gamma_l} +\end{equation} +``` +where $\delta$ is defined as +```math +\begin{equation} + \delta = q_v - q_{sl} +\end{equation} +``` + +The same can be done for deposition/sublimation, with an addition of the correction for the difference between the specific humidities of supersaturation for liquid and for ice: + +```math +\begin{equation} + \left( \frac{d q_i}{dt} \right) _{cond} = \frac{\bar{\delta_l}}{\tau_i \Gamma_i} + \frac{q_{sl} - q_{si}}{\tau_i \Gamma_i} +\end{equation} +``` + +In order to calculate $\bar{\delta}$, we consider the derivative of delta: + +```math +\begin{equation} + \frac{d \delta}{dt} = \frac{dq_v}{dt} - \frac{dq_{sl}}{dt} +\end{equation} +``` + +We can write the derivation for the saturated specific humidity: + +```math +\begin{equation} + \frac{d q_{s}}{dt} = \frac{dq_{s}}{dT} \frac{dT}{dt} + \frac{dq_{s}}{dp} \frac{dp}{dt} = \frac{dq_{s}}{dT} \frac{dT}{dt} + \frac{q_{s} \rho_a g w}{p - e} +\end{equation} +``` +where +- $\rho_a$ is the air density +- $p$ is the air pressure +- $e$ is the saturation vapor pressure +- This makes sense given the assumptions that $\frac{dp}{dt} = -\rho_a gw$ from hydrostatic balance, and $q_s = \frac{e_s}{p-e_s}$ for water vapor, where $e_s$ is the saturation water vapor pressure. Then $\frac{d q_s}{dp} = - \frac{e_s}{(p-e_s)^2} = - \frac{q_s}{p-e_s}$ . + +and they write the change in temperature with time as: + +```math +\begin{equation} + \frac{dT}{dt} = - \frac{g w}{c_p} + \frac{L_v}{c_p} \frac{\delta}{\tau_l \Gamma_l} + \frac{L_s}{c_p} \left(\frac{\delta}{\tau_i \Gamma_i} + \frac{q_{sl} - q_{si}}{\tau_i \Gamma_i} \right) +\end{equation} +``` + +Then + +```math +\begin{equation} + \frac{d q_{sl}}{dt} = \frac{dq_{sl}}{dT} \left[- \frac{g w}{c_p} + \frac{L_v}{c_p} \frac{\delta}{\tau_l \Gamma_l} + \frac{L_s}{c_p} \left(\frac{\delta}{\tau_i \Gamma_i} + \frac{q_{sl} - q_{si}}{\tau_i \Gamma_i} \right) \right] + \frac{q_{sl} \rho_a g w}{p - e} +\end{equation} +``` + +!!! note + we neglect terms due to radiation and mixing (included in the [MorrisonGrabowski2008_supersat](@cite) and [MorrisonMilbrandt2015](@cite) versions) in the parcel case. + +Putting these back into the $\delta$ equations and rearranging based on the definitions of $\Gamma_l$ and $\Gamma_i$, we get: + +```math +\begin{equation} + \frac{d \delta}{dt} = \frac{- q_{sl} \rho_a g w}{p - e_{sl}} + \frac{gw}{c_p} \frac{dq_{sl}}{dT} - \frac{\delta}{\tau} +\end{equation} +``` + +where $\tau$ is + +```math +\begin{equation} + \tau = \left( \tau_{l}^{-1} + \left( 1 + \frac{L_{s}}{c_p} \frac{dq_{sl}}{dT} \right) \frac{\tau_i^{-1}}{\Gamma_i} \right)^{-1} +\end{equation} +``` + +where + +- $\tau_l$ is the relaxation timescale for liquid droplets +- $\tau_i$ is the relaxation timescale for ice droplets +- $L_s$ is the latent heat of sublimation +- $c_p$ is the specific heat of air at constant pressure +- $T$ is temperature + +If we assume that changes in $\frac{dq_s}{dT}$ and $\tau$ are small in comparison to their magnitudes, then we can approximate them as constant and both $\delta$ equations are linear differential equations with a solution of the form: + +```math +\begin{equation} + \delta (t) = A_c \tau + (\delta_{t=0} - A_c \tau) e^{-t/\tau} +\end{equation} +``` + +Then the $A_c$ values can be calculated based on the previous derivatives. + +Condensation Ac: +```math +\begin{equation} + A_c = - \frac{q_{sl} \rho g w}{p - e_s} + \frac{dq_{sl}}{dT} \frac{wg}{c_p} - \frac{(q_{sl} - q_{si})}{\tau_i \Gamma_i} \left( 1 + \frac{L_s}{c_p} \frac{dq_{sl}}{dT} \right) +\end{equation} +``` + +Finally, to get the values for condensation/evaporation and deposition/sublimation over the timestep, we calculate the time average of $\delta$ by dividing by $\Delta t$ (the timestep) and then by the relaxation time as described in the above equations. + +```math +\begin{equation} + \left( \frac{d q_l}{dt} \right)_{cond} = \frac{A_c \tau}{\tau_i \Gamma_l} + (\delta_{t=0} - A_c \tau) \frac{\tau}{\Delta t \tau_i \Gamma_l} (1-e^{-\Delta t / \tau} ) +\end{equation} +``` + +```math +\begin{equation} + \left( \frac{d q_{si}}{dt} \right)_{dep} = A_c \frac{\tau}{\tau_i \Gamma_i} + (\delta_{t=0} - A_c \tau) \frac{\tau}{\Delta t \tau_i \Gamma_i} (1-e^{-\Delta t / \tau} ) + \frac{(q_{sl} - q_{si})}{\tau_i \Gamma_i} +\end{equation} +``` \ No newline at end of file diff --git a/docs/src/P3Scheme.md b/docs/src/P3Scheme.md index 76f36896b..1f03e4b37 100644 --- a/docs/src/P3Scheme.md +++ b/docs/src/P3Scheme.md @@ -4,7 +4,7 @@ The `P3Scheme.jl` module implements the predicted particle properties (P3) scheme for ice-phase microphysics developed by [MorrisonMilbrandt2015](@cite). The P3 scheme is a 2-moment, bulk scheme involving a single ice-phase category with 4 degrees of freedom: total ice content, - rime content, rime volume, and number mixing ratios. + rime content, rime volume, and number concentration. Traditionally, cloud ice microphysics schemes use various predefined categories (such as ice, graupel, or hail) to represent ice modes, but the P3 scheme sidesteps the problem of prescribing transitions between ice categories by adopting @@ -13,7 +13,7 @@ This simplification aids in attempts to constrain the scheme's free parameters. The prognostic variables are: - ``N_{ice}`` - number concentration 1/m3 - - ``L_{ice}`` - ice mass content kg/m3 + - ``L_{p3, tot}`` - total ice particle mass content kg/m3 - ``L_{rim}`` - rime mass content kg/m3 - ``B_{rim}`` - rime volume - (volume of rime per total air volume: dimensionless) @@ -39,7 +39,7 @@ The mass ``m`` and projected area ``A`` of particles |large, unrimed ice | ``L_{rim} = 0`` and ``D > D_{th}`` | ``\alpha_{va} \ D^{\beta_{va}}`` | ``\gamma \ D^{\sigma}`` | |dense nonspherical ice | ``L_{rim} > 0`` and ``D_{gr} > D > D_{th}`` | ``\alpha_{va} \ D^{\beta_{va}}`` | ``\gamma \ D^{\sigma}`` | |graupel (completely rimed, spherical) | ``L_{rim} > 0`` and ``D_{cr} > D > D_{gr}`` | ``\frac{\pi}{6} \rho_g \ D^3`` | ``\frac{\pi}{4} D^2`` | -|partially rimed ice | ``L_{rim} > 0`` and ``D > D_{cr}`` | ``\frac{\alpha_{va}}{1-F_r} D^{\beta_{va}}`` | ``F_{r} \frac{\pi}{4} D^2 + (1-F_{r})\gamma \ D^{\sigma}`` | +|partially rimed ice | ``L_{rim} > 0`` and ``D > D_{cr}`` | ``\frac{\alpha_{va}}{1-F_{rim}} D^{\beta_{va}}`` | ``F_{rim} \frac{\pi}{4} D^2 + (1-F_{rim})\gamma \ D^{\sigma}`` | where: - ``D_{th}``, ``D_{gr}``, ``D_{cr}`` are particle size thresholds in ``m``, @@ -57,11 +57,11 @@ The remaining thresholds: ``D_{gr}``, ``D_{cr}``, as well as the 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`` + - ``D_{cr} = [ (\frac{1}{1-F_{rim}}) \frac{6 \alpha_{va}}{\pi \rho_g} ]^{\frac{1}{3 - \beta_{va}}}`` + - ``\rho_g = \rho_r F_{rim} + (1 - F_{rim}) \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{L_{rim}}{L_{ice}}`` is the rime mass fraction, + - ``F_{rim} = \frac{L_{rim}}{L_{ice}}`` is the rime mass fraction and ``L_{ice}`` the solid ice mass content, - ``\rho_{r} = \frac{L_{rim}}{B_{rim}}`` is the predicted rime density. !!! note @@ -106,7 +106,7 @@ The model predicted ice number concentration and ice content are defined as N_{ice} = \int_{0}^{\infty} \! N'(D) \mathrm{d}D ``` ```math -L_{ice} = \int_{0}^{\infty} \! m(D) N'(D) \mathrm{d}D +L_{p3, tot} = \int_{0}^{\infty} \! m(D) N'(D) \mathrm{d}D ``` ## Calculating shape parameters @@ -119,21 +119,29 @@ Solving for ``N_{ice}`` is relatively straightforward: N_{ice} = \int_{0}^{\infty} \! N'(D) \mathrm{d}D = \int_{0}^{\infty} \! N_{0} D^\mu \, e^{-\lambda \, D} \mathrm{d}D = N_{0} (\lambda \,^{-(\mu \, + 1)} \Gamma \,(1 + \mu \,)) ``` -``L_{ice}`` depends on the variable mass-size relation ``m(D)`` defined above. -We solve for ``L_{ice}`` in a piece-wise fashion defined by the same thresholds as ``m(D)``. -As a result ``L_{ice}`` can be expressed as a sum of incomplete gamma functions, +``L_{p3, tot}`` depends on the variable mass-size relation ``m(D)`` defined above. +We solve for ``L_{p3, tot}`` in a piece-wise fashion defined by the same thresholds as ``m(D)``. +As a result ``L_{p3, tot}`` can be expressed as a sum of incomplete gamma functions, and the shape parameters are found using iterative solver. -| condition(s) | ``L_{ice} = \int \! m(D) N'(D) \mathrm{d}D`` | gamma representation | +| condition(s) | ``L_{p3, tot} = \int \! m(D) N'(D) \mathrm{d}D`` | gamma representation | |:---------------------------------------------|:-----------------------------------------------------------------------------------------|:---------------------------------------------| | ``D < D_{th}`` | ``\int_{0}^{D_{th}} \! \frac{\pi}{6} \rho_i \ D^3 N'(D) \mathrm{d}D`` | ``\frac{\pi}{6} \rho_i N_0 \lambda \,^{-(\mu \, + 4)} (\Gamma \,(\mu \, + 4) - \Gamma \,(\mu \, + 4, \lambda \,D_{th}))``| | ``L_{rim} = 0`` and ``D > D_{th}`` | ``\int_{D_{th}}^{\infty} \! \alpha_{va} \ D^{\beta_{va}} N'(D) \mathrm{d}D`` | ``\alpha_{va} \ N_0 \lambda \,^{-(\mu \, + \beta_{va} \, + 1)} (\Gamma \,(\mu \, + \beta_{va} \, + 1, \lambda \,D_{th}))`` | | ``L_{rim} > 0`` and ``D_{gr} > D > D_{th}`` | ``\int_{D_{th}}^{D_{gr}} \! \alpha_{va} \ D^{\beta_{va}} N'(D) \mathrm{d}D`` | ``\alpha_{va} \ N_0 \lambda \,^{-(\mu \, + \beta_{va} \, + 1)} (\Gamma \,(\mu \, + \beta_{va} \, + 1, \lambda \,D_{th}) - \Gamma \,(\mu \, + \beta_{va} \, + 1, \lambda \,D_{gr}))`` | | ``L_{rim} > 0`` and ``D_{cr} > D > D_{gr}`` | ``\int_{D_{gr}}^{D_{cr}} \! \frac{\pi}{6} \rho_g \ D^3 N'(D) \mathrm{d}D`` | ``\frac{\pi}{6} \rho_g N_0 \lambda \,^{-(\mu \, + 4)} (\Gamma \,(\mu \, + 4, \lambda \,D_{gr}) - \Gamma \,(\mu \, + 4, \lambda \,D_{cr}))`` | -| ``L_{rim} > 0`` and ``D > D_{cr}`` | ``\int_{D_{cr}}^{\infty} \! \frac{\alpha_{va}}{1-F_r} D^{\beta_{va}} N'(D) \mathrm{d}D`` | ``\frac{\alpha_{va}}{1-F_r} N_0 \lambda \,^{-(\mu \, + \beta_{va} \, + 1)} (\Gamma \,(\mu \, + \beta_{va} \, + 1, \lambda \,D_{cr}))`` | +| ``L_{rim} > 0`` and ``D > D_{cr}`` | ``\int_{D_{cr}}^{\infty} \! \frac{\alpha_{va}}{1-F_{rim}} D^{\beta_{va}} N'(D) \mathrm{d}D`` | ``\frac{\alpha_{va}}{1-F_{rim}} N_0 \lambda \,^{-(\mu \, + \beta_{va} \, + 1)} (\Gamma \,(\mu \, + \beta_{va} \, + 1, \lambda \,D_{cr}))`` | where ``\Gamma \,(a, z) = \int_{z}^{\infty} \! t^{a - 1} e^{-t} \mathrm{d}D`` and ``\Gamma \,(a) = \Gamma \,(a, 0)`` for simplicity. +We rely on [Blahak2010](@cite) parameterization to compute the incomplete gamma functions. +Below plot roughly reproduce Figure 4 from [Blahak2010](@cite). +The relative error differences are mostly due to too small values returned by + our reference incomplete gamma function from Julia SpecialFunctions package. +```@example +include("plots/P3GammaAprox.jl") +``` +![](P3_GammaInc_error.svg) In our solver, we approximate ``\mu`` from ``L/N`` and keep it constant throughout the solving step. We approximate ``\mu`` by an exponential function given by the ``L/N`` points @@ -174,14 +182,20 @@ include("plots/P3LambdaErrorPlots.jl") ## Terminal Velocity -We use the [Chen2022](@cite) velocity parametrization: +We use the [Chen2022](@cite) velocity parameterization: ```math V(D) = \phi^{\kappa} \sum_{i=1}^{j} \; a_i D^{b_i} e^{-c_i \; D} ``` -where ``\phi = (16 \rho_{ice}^2 A(D)^3) / (9 \pi m(D)^2)`` is the aspect ratio, -and ``\kappa``, ``a_i``, ``b_i`` and ``c_i`` are the free parameters. +where +```math +\phi = (16 \rho_{ice}(D)^2 A(D)^3) / (9 \pi m(D)^2) +``` +is the aspect ratio, and ``\kappa``, ``a_i``, ``b_i`` and ``c_i`` are the free parameters. -The aspect ratio of a spheroid is defined ``\phi = \frac{a}{c}``, where ``a`` is the equatorial radius and ``c`` is the distance from the pole to the center. In terms of ``a`` and ``c``, a spheroid's volume can be represented as ``V = \frac{4}{3} \pi a^2 c``, and its cross-sectional area can be assumed ``A(a, c) = \pi a c``. We use ``m(D)`` and ``A(D)`` from P3, so by substituting ``m(a, c) = \rho_{ice} * V(a, c)``, ``A(a, c)`` for ``m(D)``, ``A(D)`` into the formulation of aspect ratio above, we demonstrate agreement with the definition ``\phi = \frac{a}{c}``. +The aspect ratio of a spheroid is defined ``\phi = \frac{a}{c}``, where ``a`` is the equatorial radius and ``c`` is the distance from the pole to the center. + In terms of ``a`` and ``c``, a spheroid's volume can be represented as ``V(a, c) = \frac{4}{3} \pi a^2 c``, and its cross-sectional area can be assumed ``A(a, c) = \pi a c``. + We use ``m(D)`` and ``A(D)`` from P3, so by substituting ``m(a, c) = \rho_{ice} V(a, c)``, ``A(a, c)`` for ``m(D)``, ``A(D)`` into the formulation of aspect ratio above, + we demonstrate agreement with the definition ``\phi = \frac{a}{c}``. Note that ``\phi = 1`` corresponds to spherical particles (small spherical ice (``D < D_{th}``) and graupel (``D_{gr} < D < D_{cr}``)). @@ -208,13 +222,100 @@ Below we provide plots of these relationships for small, medium, and large ``D_m the first row highlights the particle size regime, the second displays ``D_m`` of the particles, the third shows the aspect ratio ``\phi (D_m)``, - and the final row exhibits ``V_m``. + the fourth shows ``V_m`` without using aspect ratio in the computation (i.e. ``\phi = 1``), + and the final row exhibits ``V_m``, computed using ``\phi(D)``. They can be compared with Figure 2 from [MorrisonMilbrandt2015](@cite). ```@example include("plots/P3TerminalVelocityPlots.jl") ``` ![](MorrisonandMilbrandtFig2.svg) +## Liquid Fraction + +To allow for the modeling of mixed-phase particles with P3, a new prognostic variable can be introduced: + ``L_{liq}``, the content of liquid on mixed-phase particles. As described in [Choletteetal2019](@cite), + this addition to the framework of P3 opens the door to tracking the gradual melting (and refreezing) of a + particle population. Here, we describe the characteristics of the scheme with the addition of ``L_{liq}``. + +Liquid fraction, analogous to ``F_{rim}`` for rime, is defined ``F_{liq} = \frac{L_{liq}}{L_{p3, tot}}`` + where ``L_{p3, tot} = L_{ice} + L_{liq}`` and where ``L_{ice}`` is the solid ice mass content which includes + rime mass content and mass grown by vapor deposition. This is important notably for + ``F_{rim} = \frac{L_{rim}}{L_{ice}}`` — which differs now from ``F_{rim} = \frac{L_{rim}}{L_{p3, tot}}`` + because we want to normalize only by solid ice mass content for ``F_{rim}``. + +Based on Fig. 1 from [Choletteetal2019](@cite), we can expect the accumulation of liquid on an ice core to increase + velocities of small particles for all ``F_{rim}`` values. Below, we reproduce this figure with ``\rho_{r} = 900 kg m^{-3}``, + and, notably, we continue to use the terminal velocity parameterizations from [Chen2022](@cite) + (described [here](https://clima.github.io/CloudMicrophysics.jl/dev/Microphysics2M/#Terminal-Velocity) for rain), whereas [Choletteetal2019](@cite) uses + [MitchellHeymsfield2005](@cite) for snow and ice and [Simmeletal2002](@cite) for rain. Despite these different choices, we + reproduce similar behavior with, and we include a dashed line for the velocity computed without aspect ratio. This dashed line mimics velocity + with ``F_{rim} = 1``, since both ``\phi = 1`` and ``F_{rim} = 1`` shift us into a spherical particle regime. + +```@example +include("plots/Cholette2019_fig1.jl") +``` +![](Choletteetal2019_fig1.svg) + +The addition of the liquid fraction does not change the thresholds ``D_{th}``, ``D_{gr}``, ``D_{cr}``, + since the threshold regime depends only on ice core properties. + +However, the assumed particle properties become ``F_{liq}``-weighted averages of particles' solid and liquid + components: + +```math +m(D, F_{liq}) = (1 - F_{liq}) m(D, F_{liq} = 0) + F_{liq} m_{liq}(D) +``` +```math +A(D, F_{liq}) = (1 - F_liq) A(D, F_{liq} = 0) + F_{liq} A_{liq}(D) +``` + +where ``m_{liq}(D) = \frac{\pi}{6} \rho_{liq} D^3`` and ``A_{liq}(D) = \frac{\pi}{4} D^2``. + +When calculating shape parameters and integrating over the particle size distribution (PSD), it is important to + keep in mind whether the desired moment of the PSD is tied only to ice (in which case we concern ourselves with the + ice core diameter ``D_{core}``) or to the whole mixed-phase particle (in which case we need ``D_p``). If calculating + the ice core parameters, ``F_{liq} = 0`` is passed into the solver framework as indicated in the code. + +For the above particle properties and for terminal velocity, we use the PSD corresponding to the whole mixed-phase particle, + so our terminal velocity of a mixed-phase particle is: + +```math +V(D_{p}, F_{liq}) = (1 - F_{liq}) V_{ice}(D_{p}) + F_{liq} V_{rain}(D_{p}) +``` + +To continue with the same plotting format as we see above for terminal velocity, below, we show terminal velocity with + ``F_{liq} = 0.0, 0.5, 0.9`` using the mass-weighted terminal velocity with aspect ratio. Clearly, the velocity increases + with ``F_{liq}`` and becomes decreasingly dependent on ``F_{rim}`` and ``\rho_{r}`` as we shrink the size of the ice core. + Of note relating to our calculation of terminal velocity with liquid fraction, we let ``F_{liq} = 0`` in our calculation of + ``\phi``. This is because we want ``V_{ice}(D_{p})`` (velocity of a mixed-phase particle treating it as an ice particle + with the same D). + +```@example +include("plots/P3TerminalVelocityPlots_WithFliq.jl") +``` +![](MorrisonandMilbrandtFig2_with_F_liq.svg) + +Visualizing mass-weighted terminal velocity as a function of ``F_{liq}``, ``F_{rim}`` with ``\rho_{r} = 900 kg m^{-3}`` for small, + medium, and large particles, we have mostly graupel (``D_{gr} < D < D_{cr}``) for small ``D_m`` and mostly partially rimed ice + (``D > D_{cr}``) for medium and large ``D_m``. Thus, we can attribute the non-monotonic behavior of velocity with ``F_liq`` in the + medium and large ``D_m`` plots to the variations in ``\phi`` caused by nonspherical particle shape, whereas the small ``D_m`` plot + confirms a more monotonic change in ``V_m`` for spherical ice. The ``L`` and ``N`` values used to generate small, medium, and large ``D_{m}`` + are the same as in the plot above. + +```@example +include("plots/P3TerminalVelocity_F_liq_rim.jl") +``` +![](P3TerminalVelocity_F_liq_rim.svg) + +When modifying process rates, we now need to consider whether they are concerned with + the ice core or the whole particle, in addition to whether they become sources + and sinks of different prognostic variables with the inclusion of ``F_{liq}``. + With the addition of liquid fraction, too, + come new process rates. + +!!! note + TODO - Implement process rates, complete docs. + ## Heterogeneous Freezing Immersion freezing is parameterized based on water activity and follows the ABIFM @@ -243,6 +344,36 @@ include("plots/P3ImmersionFreezing.jl") ``` ![](P3_het_ice_nucleation.svg) +## Melting + +Melting rate is derived in the same way as in the + [1-moment scheme](https://clima.github.io/CloudMicrophysics.jl/dev/Microphysics1M/#Snow-melt). +We assume the same ventilation factor parameterization as in [SeifertBeheng2006](@cite), + and use the terminal velocity parameterization from [Chen2022](@cite). +The ``dm/dD`` derivative is computed for each P3 size regime. +The bulk melting rate is computed by numerically integrating over the particle + size distribution: +```math +\begin{equation} + \left. \frac{dL}{dt} \right|_{melt} = \frac{4 \, K_{thermo}}{L_f} \left(T - T_{freeze}\right) + \int_{0}^{\infty} \frac{dm(D)}{dD} \frac{F_v(D) N(D)}{D} +\end{equation} +``` +The melting rate for number concentration is assumed to be proportional to the + ice content melting rate. +```math +\begin{equation} + \left. \frac{dN}{dt} \right|_{melt} = \frac{N}{L} \left. \frac{dL}{dt} \right|_{melt} +\end{equation} +``` +Both rates are limited by the total available ice content and number concentration + divided by model time step length. + +```@example +include("plots/P3Melting.jl") +``` +![](P3_ice_melt.svg) + ## Acknowledgments Click on the P3 mascot duck to be taken to the repository diff --git a/docs/src/plots/Cholette2019_fig1.jl b/docs/src/plots/Cholette2019_fig1.jl new file mode 100644 index 000000000..16cce4136 --- /dev/null +++ b/docs/src/plots/Cholette2019_fig1.jl @@ -0,0 +1,162 @@ +import ClimaParams +import CloudMicrophysics as CM +import CloudMicrophysics.P3Scheme as P3 +import CloudMicrophysics.Parameters as CMP +import CloudMicrophysics.Common as CO +import CairoMakie as Plt + +const PSP3 = CMP.ParametersP3 + +FT = Float64 + +p3 = CMP.ParametersP3(FT) + +# Testing terminal velocity with liquid fraction + +function get_values( + p3::PSP3, + Chen2022::CMP.Chen2022VelType, + F_liq::FT, + F_r::FT, + ρ_a::FT, + th, + res::Int, +) where {FT <: Real} + + # particle dimension from 0 to 1 cm + D_ps = range(FT(0), stop = FT(0.01), length = res) + + # ρ_r = 900 + ρ_r = FT(900) + + V = zeros(res) + Vϕ = zeros(res) + do_not_use_aspect = false + use_aspect_ratio = true + + for i in 1:res + Vϕ[i] = P3.p3_particle_terminal_velocity( + p3, + D_ps[i], + Chen2022, + ρ_a, + F_r, + F_liq, + th, + use_aspect_ratio, + ) + V[i] = P3.p3_particle_terminal_velocity( + p3, + D_ps[i], + Chen2022, + ρ_a, + F_r, + F_liq, + th, + do_not_use_aspect, + ) + end + return (D_ps, V, Vϕ) +end + +function fig1() + Chen2022 = CMP.Chen2022VelType(FT) + ρ_a = FT(1.2) + ρ_r = FT(900) + + F_liqs = (FT(0), FT(0.33), FT(0.67), FT(1)) + F_rs = (FT(0), FT(1 - eps(FT))) + + labels = + ["F_liq = 0", "F_liq = 0.33", "F_liq = 0.67", "F_liq = 1", "ϕᵢ = 1"] + colors = [:black, :blue, :red, :green] + res = 50 + + fig = Plt.Figure(size = (1200, 400)) + + #F_r = 0 + ax1 = Plt.Axis( + fig[1:7, 1:9], + title = "Fᵣ = 0", + xlabel = "Dₚ (m)", + ylabel = "V (m s⁻¹)", + ) + + for i in 1:4 + D_ps, V_0, V_0_ϕ = get_values( + p3, + Chen2022, + F_liqs[i], + F_rs[1], + ρ_a, + P3.thresholds(p3, ρ_r, F_rs[1]), + res, + ) + Plt.lines!(ax1, D_ps, V_0_ϕ, color = colors[i], label = labels[i]) + Plt.lines!( + ax1, + D_ps, + V_0, + color = colors[i], + label = labels[i], + linestyle = :dash, + ) + end + + # F_r = 1 + ax2 = Plt.Axis( + fig[1:7, 10:18], + title = "Fᵣ = 1", + xlabel = "Dₚ (m)", + ylabel = "V (m s⁻¹)", + ) + for i in 1:4 + D_ps, V_1, V_1_ϕ = get_values( + p3, + Chen2022, + F_liqs[i], + F_rs[2], + ρ_a, + P3.thresholds(p3, ρ_r, F_rs[2]), + res, + ) + Plt.lines!(ax2, D_ps, V_1_ϕ, color = colors[i], label = labels[i]) + Plt.lines!( + ax2, + D_ps, + V_1, + color = colors[i], + label = labels[i], + linestyle = :dash, + ) + end + + Plt.Legend( + fig[4, 19:22], + [ + Plt.LineElement(color = :black), + Plt.LineElement(color = :blue), + Plt.LineElement(color = :red), + Plt.LineElement(color = :green), + Plt.LineElement(color = :black, linestyle = :dash), + ], + labels, + framevisible = false, + ) + + Plt.linkaxes!(ax1, ax2) + ax1.xticks = ( + [0, 0.002, 0.004, 0.006, 0.008, 0.01], + string.([0, 0.002, 0.004, 0.006, 0.008, 0.01]), + ) + ax2.xticks = ( + [0, 0.002, 0.004, 0.006, 0.008, 0.01], + string.([0, 0.002, 0.004, 0.006, 0.008, 0.01]), + ) + + + Plt.resize_to_layout!(fig) + Plt.save("Choletteetal2019_fig1.svg", fig) +end + +fig1() diff --git a/docs/src/plots/P3AspectRatioPlot.jl b/docs/src/plots/P3AspectRatioPlot.jl index 9ae6090e9..816c99a7f 100644 --- a/docs/src/plots/P3AspectRatioPlot.jl +++ b/docs/src/plots/P3AspectRatioPlot.jl @@ -50,7 +50,7 @@ function p3_aspect() # define plot axis #[row, column] ax1 = define_axis(fig, 1:7, 1:9, CMK.L"ϕᵢ(D) regime for $ρ_r = 400 kg m^{-3}$", CMK.L"$ϕᵢ$ (-)", [0.0, 0.5, 1.0], 1.9, logscale = false) - ax3 = define_axis(fig, 1:7, 10:18, CMK.L"ϕᵢ(D) regime for $F_r = 0.4$", CMK.L"$ϕᵢ$ (-)", [0.0, 0.5, 1.0], 1.8, logscale = false) + ax3 = define_axis(fig, 1:7, 10:18, CMK.L"ϕᵢ(D) regime for $F_rim = 0.4$", CMK.L"$ϕᵢ$ (-)", [0.0, 0.5, 1.0], 1.8, logscale = false) # Get thresholds sol4_0 = P3.thresholds(p3, 400.0, 0.0) @@ -87,10 +87,10 @@ function p3_aspect() global d_gr_800 = CMK.vlines!(ax, sol_8[2] * 1e3, linestyle = :dash, color = cl[3], linewidth = lw) end # add legend - CMK.Legend(fig[8:9, 1], [fig1_0, fig1_5, fig1_8], [CMK.L"$F_{r} = 0.0$", CMK.L"$F_{r} = 0.5$", CMK.L"$F_{r} = 0.8$"], framevisible = false) + CMK.Legend(fig[8:9, 1], [fig1_0, fig1_5, fig1_8], [CMK.L"$F_{rim} = 0.0$", CMK.L"$F_{rim} = 0.5$", CMK.L"$F_{rim} = 0.8$"], framevisible = false) CMK.Legend(fig[8:9, 3], [d_tha], [CMK.L"$D_{th}$"], framevisible = false) - CMK.Legend(fig[8:9, 7], [d_cr_5, d_cr_8], [CMK.L"$D_{cr}$ for $F_{r} = 0.5$", CMK.L"$D_{cr}$ for $F_{r} = 0.8$"], framevisible = false) - CMK.Legend(fig[8:9, 5], [d_gr_5, d_gr_8], [CMK.L"$D_{gr}$ for $F_{r} = 0.5$", CMK.L"$D_{gr}$ for $F_{r} = 0.8$"], framevisible = false) + CMK.Legend(fig[8:9, 7], [d_cr_5, d_cr_8], [CMK.L"$D_{cr}$ for $F_{rim} = 0.5$", CMK.L"$D_{cr}$ for $F_{rim} = 0.8$"], framevisible = false) + CMK.Legend(fig[8:9, 5], [d_gr_5, d_gr_8], [CMK.L"$D_{gr}$ for $F_{rim} = 0.5$", CMK.L"$D_{gr}$ for $F_{rim} = 0.8$"], framevisible = false) CMK.Legend(fig[8:9, 13], [fig2_200, fig2_400, fig2_800], [CMK.L"$\rho_{r} = 200.0 kg m^{-3}$", CMK.L"$\rho_{r} = 400.0 kg m^{-3}$", CMK.L"$\rho_{r} = 800.0 kg m^{-3}$",], framevisible = false) CMK.Legend(fig[8:9, 14], [d_thb], [CMK.L"$D_{th}$"], framevisible = false) diff --git a/docs/src/plots/P3GammaAprox.jl b/docs/src/plots/P3GammaAprox.jl new file mode 100644 index 000000000..3779776e0 --- /dev/null +++ b/docs/src/plots/P3GammaAprox.jl @@ -0,0 +1,145 @@ +import CairoMakie as PL +PL.activate!(type = "svg") + +import SpecialFunctions as SF + +import Thermodynamics as TD +import CloudMicrophysics.Parameters as CMP +import CloudMicrophysics.P3Scheme as P3 + +#! format: off + +FT = Float64 +res = 100 # plot resolution + +# Plot Fig 3 from Blahak 2010 +# --------------------------- + +a_range = range(1e-1, stop = 30, length = res) +c1 = [P3.c₁(a) for a in a_range] +c2 = [P3.c₂(a) for a in a_range] +c3 = [P3.c₃(a) for a in a_range] +c4 = [P3.c₄(a) for a in a_range] + +fig = PL.Figure(size = (1500, 1000), fontsize=22, linewidth=3) + +ax1 = PL.Axis(fig[1, 1]) +ax2 = PL.Axis(fig[1, 2]) +ax3 = PL.Axis(fig[2, 1]) +ax4 = PL.Axis(fig[2, 2]) + +ax1.ylabel = "c₁" +ax2.ylabel = "c₂" +ax3.ylabel = "c₃" +ax4.ylabel = "c₄" + +ax3.xlabel = "a" +ax4.xlabel = "a" + +PL.ylims!(ax2, 0, 2) +PL.ylims!(ax4, 0, 10) + +c1l = PL.lines!(ax1, a_range, c1) +c2l = PL.lines!(ax2, a_range, c2) +c3l = PL.lines!(ax3, a_range, c3) +c4l = PL.lines!(ax4, a_range, c4) + +PL.save("P3_GammaInc_c_coeffs.svg", fig) + + +# Plot figure 1 from Blahak 2010 +# ------------------------------ + +fig2 = PL.Figure(size = (1000, 500), fontsize=22) +ax = PL.Axis(fig2[1,1]) +ax.xlabel = "z" +ax.ylabel = "P" + +z_range = range(1e-1, stop = 25, length=res) + +P_a1 = [P3.Γ_lower(FT(1), z) / SF.gamma(FT(1)) for z in z_range] +P_a3 = [P3.Γ_lower(FT(3), z) / SF.gamma(FT(3)) for z in z_range] +P_a5 = [P3.Γ_lower(FT(5), z) / SF.gamma(FT(5)) for z in z_range] +P_a10 = [P3.Γ_lower(FT(10), z) / SF.gamma(FT(10)) for z in z_range] +P_a15 = [P3.Γ_lower(FT(15), z) / SF.gamma(FT(15)) for z in z_range] +P_a30 = [P3.Γ_lower(FT(30), z) / SF.gamma(FT(30)) for z in z_range] + +G_a1 = [SF.gamma_inc(FT(1), z)[1] for z in z_range] +G_a3 = [SF.gamma_inc(FT(3), z)[1] for z in z_range] +G_a5 = [SF.gamma_inc(FT(5), z)[1] for z in z_range] +G_a10 = [SF.gamma_inc(FT(10), z)[1] for z in z_range] +G_a15 = [SF.gamma_inc(FT(15), z)[1] for z in z_range] +G_a30 = [SF.gamma_inc(FT(30), z)[1] for z in z_range] + +l1 = PL.lines!(ax, z_range, P_a1, label="a=1", color="midnightblue", linewidth=3) +l3 = PL.lines!(ax, z_range, P_a3, label="a=3", color="slateblue", linewidth=3) +l5 = PL.lines!(ax, z_range, P_a5, label="a=5", color="steelblue1", linewidth=3) +l10 = PL.lines!(ax, z_range, P_a10, label="a=10", color="orange", linewidth=3) +l15 = PL.lines!(ax, z_range, P_a15, label="a=15", color="orangered2", linewidth=3) +l30 = PL.lines!(ax, z_range, P_a30, label="a=30", color="brown", linewidth=3) + +g1 = PL.lines!(ax, z_range, G_a1, label="a=1", color="midnightblue", linewidth=3, linestyle=:dash) +g3 = PL.lines!(ax, z_range, G_a3, label="a=1", color="slateblue", linewidth=3, linestyle=:dash) +g5 = PL.lines!(ax, z_range, G_a5, label="a=1", color="steelblue1", linewidth=3, linestyle=:dash) +g10 = PL.lines!(ax, z_range, G_a10, label="a=1", color="orange", linewidth=3, linestyle=:dash) +g15 = PL.lines!(ax, z_range, G_a15, label="a=1", color="orangered2", linewidth=3, linestyle=:dash) +g30 = PL.lines!(ax, z_range, G_a30, label="a=1", color="brown", linewidth=3, linestyle=:dash) + +PL.Legend( + fig2[1, 2], + [l1, l3, l5, l10, l15, l30], + ["a=1", "a=3", "a=5", "a=10", "a=15", "a=30"], +) +PL.save("P3_GammaInc_P_check.svg", fig2) + +# Plot figure 4 from Blahak 2010 +# ------------------------------ +# The relative error blows up for cases where the Γ_lower is very small. +# The reference value we obtain from Julia Special Functions is off +# by several orders of magnitude when compared with Wolfram Mathematica + +y_range = range(1e-1, stop = 2.5, length = res) +z_range = y_range .* (a_range .+ 1) + +val = zeros(res, res) +err_rel = zeros(res, res) +err_abs = zeros(res, res) +for it in range(1, stop=res, step=1) # a + for jt in range(1, stop=res, step=1) # z + a = FT(a_range[it]) + z = FT(z_range[jt]) + val[it, jt] = P3.Γ_lower(a, z) / SF.gamma(a) + err_abs[it, jt] = (P3.Γ_lower(a, z) / SF.gamma(a)) - SF.gamma_inc(a, z)[1] + err_rel[it, jt] = (P3.Γ_lower(a, z) / SF.gamma(a)) / SF.gamma_inc(a, z)[1] - 1 + end +end + +fig3 = PL.Figure(size = (1800, 600), fontsize=22) +ax1 = PL.Axis(fig3[1, 1], xlabel = "a", ylabel = "z / (a + 1)", title = "Relative error") +ax2 = PL.Axis(fig3[1, 3], xlabel = "a", title = "Absolute error") +ax3 = PL.Axis(fig3[1, 5], xlabel = "a", title = "Γ_lower(a, z) / Γ(a)") + +PL.ylims!(ax1, 0, 2.5) +PL.ylims!(ax2, 0, 2.5) +PL.ylims!(ax3, 0, 2.5) + +cplot1 = PL.heatmap!( + ax1, a_range, y_range, err_rel, + colormap = PL.cgrad(:viridis, 20, categorical=true), + colorrange = (-0.1, 10), highclip=:red, lowclip=:orange, +) +PL.Colorbar(fig3[1,2], cplot1) + +cplot2 = PL.contourf!( + ax2, a_range, y_range, err_abs, + levels=-0.03:0.005:0.03, + colormap="RdBu" +) +PL.Colorbar(fig3[1,4], cplot2) + +cplot3 = PL.contourf!( + ax3, a_range, y_range, val, +) +PL.Colorbar(fig3[1,6], cplot3) + +PL.save("P3_GammaInc_error.svg", fig3) diff --git a/docs/src/plots/P3LambdaErrorPlots.jl b/docs/src/plots/P3LambdaErrorPlots.jl index 5aaeba16c..a692b3821 100644 --- a/docs/src/plots/P3LambdaErrorPlots.jl +++ b/docs/src/plots/P3LambdaErrorPlots.jl @@ -10,19 +10,22 @@ FT = Float64 const PSP3 = CMP.ParametersP3 p3 = CMP.ParametersP3(FT) +F_liq = FT(0) # preserving original P3 +# TODO: investigate for F_liq != 0 -function λ_diff(F_r::FT, ρ_r::FT, N::FT, λ_ex::FT, p3::PSP3) where {FT} +function λ_diff(F_rim::FT, ρ_r::FT, N::FT, λ_ex::FT, p3::PSP3) where {FT} # Find the P3 scheme thresholds - th = P3.thresholds(p3, ρ_r, F_r) + th = P3.thresholds(p3, ρ_r, F_rim) # Get μ corresponding to λ μ = P3.DSD_μ(p3, λ_ex) # Convert λ to ensure it remains positive x = log(λ_ex) # Compute mass density based on input shape parameters - L_calc = N * P3.L_over_N_gamma(p3, F_r, x, μ, th) + L_calc = N * P3.L_over_N_gamma(p3, F_rim, F_liq, x, μ, th) - (λ_calculated,) = P3.distribution_parameter_solver(p3, L_calc, N, ρ_r, F_r) + (λ_calculated,) = + P3.distribution_parameter_solver(p3, L_calc, N, ρ_r, F_rim, F_liq) return abs(λ_ex - λ_calculated) end @@ -30,31 +33,31 @@ function get_errors( p3::PSP3, log10_λ_min::FT, log10_λ_max::FT, - F_r_min::FT, - F_r_max::FT, + F_rim_min::FT, + F_rim_max::FT, ρ_r::FT, N::FT, λSteps::Int, - F_rSteps::Int, + F_rimSteps::Int, ) where {FT} logλs = range(FT(log10_λ_min), stop = log10_λ_max, length = λSteps) λs = [10^logλ for logλ in logλs] - F_rs = range(F_r_min, stop = F_r_max, length = F_rSteps) - E = zeros(λSteps, F_rSteps) + F_rims = range(F_rim_min, stop = F_rim_max, length = F_rimSteps) + E = zeros(λSteps, F_rimSteps) for i in 1:λSteps - for j in 1:F_rSteps + for j in 1:F_rimSteps λ = λs[i] - F_r = F_rs[j] + F_rim = F_rims[j] - er = log(λ_diff(F_r, ρ_r, N, λ, p3) / λ) + er = log(λ_diff(F_rim, ρ_r, N, λ, p3) / λ) er = er == Inf ? 9999.99 : er er = er == -Inf ? -9999.99 : er E[i, j] = er end end - return (; λs, F_rs, E) + return (; λs, F_rims, E) end #! format: off @@ -62,12 +65,12 @@ function plot_relerrors( N::FT, log10_λ_min::FT, log10_λ_max::FT, - F_r_min::FT, - F_r_max::FT, + F_rim_min::FT, + F_rim_max::FT, ρ_r_min::FT, ρ_r_max::FT, λSteps::Int, - F_rSteps::Int, + F_rimSteps::Int, numPlots::Int, p3::PSP3, ) where {FT} @@ -81,10 +84,10 @@ function plot_relerrors( i = 1 ρ = ρ_rs[i] - (λs, F_rs, E) = get_errors(p3, log10_λ_min, log10_λ_max, F_r_min, F_r_max, ρ, N, λSteps, F_rSteps) + (λs, F_rims, E) = get_errors(p3, log10_λ_min, log10_λ_max, F_rim_min, F_rim_max, ρ, N, λSteps, F_rimSteps) - CMK.Axis(f[x, y], xlabel = "λ", ylabel = "F_r", title = string("log(relative error calculated λ) for ρ_r = ", string(ρ)), width = 400, height = 300, xscale = log10) - hm = CMK.heatmap!(λs, F_rs, E, colormap = CMK.cgrad(:viridis, 20, categorical=true), colorrange = (-10, 0), highclip = :red, lowclip = :indigo) + CMK.Axis(f[x, y], xlabel = "λ", ylabel = "F_rim", title = string("log(relative error calculated λ) for ρ_r = ", string(ρ)), width = 400, height = 300, xscale = log10) + hm = CMK.heatmap!(λs, F_rims, E, colormap = CMK.cgrad(:viridis, 20, categorical=true), colorrange = (-10, 0), highclip = :red, lowclip = :indigo) CMK.Colorbar(f[x, y + 1], hm) y = y + 2 @@ -98,7 +101,7 @@ function plot_relerrors( end #! format: on -function μ_approximation_effects(F_r::FT, ρ_r::FT) where {FT} +function μ_approximation_effects(F_rim::FT, ρ_r::FT) where {FT} f = CMK.Figure(size = (800, 500)) @@ -106,7 +109,7 @@ function μ_approximation_effects(F_r::FT, ρ_r::FT) where {FT} f[1, 1], xlabel = "L/N", ylabel = "μ", - title = string("μ vs L/N for F_r = ", F_r, " ρ_r = ", ρ_r), + title = string("μ vs L/N for F_rim = ", F_rim, " ρ_r = ", ρ_r), width = 400, height = 300, xscale = log10, @@ -116,7 +119,7 @@ function μ_approximation_effects(F_r::FT, ρ_r::FT) where {FT} f[1, 2], xlabel = "λ", ylabel = "μ", - title = string("μ vs λ for F_r = ", F_r, " ρ_r = ", ρ_r), + title = string("μ vs λ for F_rim = ", F_rim, " ρ_r = ", ρ_r), width = 400, height = 300, xscale = log10, @@ -126,7 +129,7 @@ function μ_approximation_effects(F_r::FT, ρ_r::FT) where {FT} f[1, 3], xlabel = "λ", ylabel = "L/N", - title = string("L/N vs λ for F_r = ", F_r, " ρ_r = ", ρ_r), + title = string("L/N vs λ for F_rim = ", F_rim, " ρ_r = ", ρ_r), width = 400, height = 300, xscale = log10, @@ -139,7 +142,7 @@ function μ_approximation_effects(F_r::FT, ρ_r::FT) where {FT} numpts = 100 # Set up vectors - th = P3.thresholds(p3, ρ_r, F_r) + th = P3.thresholds(p3, ρ_r, F_rim) log_λs = range(FT(3.6), stop = FT(4.6), length = numpts) λs = [10^log_λ for log_λ in log_λs] μs = [P3.DSD_μ(p3, λ) for λ in λs] @@ -149,12 +152,19 @@ function μ_approximation_effects(F_r::FT, ρ_r::FT) where {FT} λ_solved = [FT(0) for λ in λs] for i in 1:numpts - L_over_N = P3.L_over_N_gamma(p3, F_r, log(λs[i]), μs[i], th) + L_over_N = P3.L_over_N_gamma(p3, F_rim, F_liq, log(λs[i]), μs[i], th) Ls[i] = L_over_N N = FT(1e6) - (L, N) = P3.distribution_parameter_solver(p3, L_over_N * N, N, ρ_r, F_r) + (L, N) = P3.distribution_parameter_solver( + p3, + L_over_N * N, + N, + ρ_r, + F_rim, + F_liq, + ) λ_solved[i] = L - μs_approx[i] = P3.DSD_μ_approx(p3, N * L_over_N, N, ρ_r, F_r) + μs_approx[i] = P3.DSD_μ_approx(p3, N * L_over_N, N, ρ_r, F_rim, F_liq) end # Plot @@ -179,31 +189,31 @@ end log10_λ_min = FT(2) log10_λ_max = FT(6) -F_r_min = FT(0) -F_r_max = FT(0.9) +F_rim_min = FT(0) +F_rim_max = FT(0.9) ρ_r_min = FT(100) ρ_r_max = FT(900) N = FT(1e8) λ_Steps = 40 -F_r_Steps = 40 +F_rim_Steps = 40 NumPlots = 9 plot_relerrors( N, log10_λ_min, log10_λ_max, - F_r_min, - F_r_max, + F_rim_min, + F_rim_max, ρ_r_min, ρ_r_max, λ_Steps, - F_r_Steps, + F_rim_Steps, NumPlots, p3, ) -F_r = FT(0.5) +F_rim = FT(0.5) ρ_r = FT(500) -μ_approximation_effects(F_r, ρ_r) +μ_approximation_effects(F_rim, ρ_r) diff --git a/docs/src/plots/P3Melting.jl b/docs/src/plots/P3Melting.jl new file mode 100644 index 000000000..ad8a8ba75 --- /dev/null +++ b/docs/src/plots/P3Melting.jl @@ -0,0 +1,89 @@ +import CairoMakie as PL +PL.activate!(type = "svg") + +import Thermodynamics as TD +import CloudMicrophysics.Parameters as CMP +import CloudMicrophysics.P3Scheme as P3 + +FT = Float64 + +# parameters +tps = TD.Parameters.ThermodynamicsParameters(FT) +pps = CMP.ParametersP3(FT) +vel = CMP.Chen2022VelType(FT) +aps = CMP.AirProperties(FT) +tps = TD.Parameters.ThermodynamicsParameters(FT) + +# model time step (for limiting) +dt = FT(1) + +# initial ice content and number concentration +Lᵢ = FT(1e-4) +Nᵢ = FT(2e5) + +# tested temperature range +ΔT_range = range(1e-4, stop = 0.025, length = 1000) + +# limiters to not melt more mass and number than we have +max_dLdt = [Lᵢ / dt for ΔT in ΔT_range] +max_dNdt = [Nᵢ / dt for ΔT in ΔT_range] + +#! format: off + +ρₐ1 = FT(1.2) +Fᵣ1 = FT(0.8) +ρᵣ1 = FT(800) +dLdt1 = [P3.ice_melt(pps, vel.snow_ice, aps, tps, Lᵢ, Nᵢ, pps.T_freeze .+ ΔT, ρₐ1, Fᵣ1, ρᵣ1, dt).dLdt for ΔT in ΔT_range] +dNdt1 = [P3.ice_melt(pps, vel.snow_ice, aps, tps, Lᵢ, Nᵢ, pps.T_freeze .+ ΔT, ρₐ1, Fᵣ1, ρᵣ1, dt).dNdt for ΔT in ΔT_range] + +Fᵣ2 = FT(0.2) +dLdt2 = [P3.ice_melt(pps, vel.snow_ice, aps, tps, Lᵢ, Nᵢ, pps.T_freeze .+ ΔT, ρₐ1, Fᵣ2, ρᵣ1, dt).dLdt for ΔT in ΔT_range] +dNdt2 = [P3.ice_melt(pps, vel.snow_ice, aps, tps, Lᵢ, Nᵢ, pps.T_freeze .+ ΔT, ρₐ1, Fᵣ2, ρᵣ1, dt).dNdt for ΔT in ΔT_range] + +ρᵣ2 = FT(200) +dLdt3 = [P3.ice_melt(pps, vel.snow_ice, aps, tps, Lᵢ, Nᵢ, pps.T_freeze .+ ΔT, ρₐ1, Fᵣ2, ρᵣ2, dt).dLdt for ΔT in ΔT_range] +dNdt3 = [P3.ice_melt(pps, vel.snow_ice, aps, tps, Lᵢ, Nᵢ, pps.T_freeze .+ ΔT, ρₐ1, Fᵣ2, ρᵣ2, dt).dNdt for ΔT in ΔT_range] + +ρₐ2 = FT(0.5) +dLdt4 = [P3.ice_melt(pps, vel.snow_ice, aps, tps, Lᵢ, Nᵢ, pps.T_freeze .+ ΔT, ρₐ2, Fᵣ2, ρᵣ2, dt).dLdt for ΔT in ΔT_range] +dNdt4 = [P3.ice_melt(pps, vel.snow_ice, aps, tps, Lᵢ, Nᵢ, pps.T_freeze .+ ΔT, ρₐ2, Fᵣ2, ρᵣ2, dt).dNdt for ΔT in ΔT_range] + +# plotting +fig = PL.Figure(size = (1500, 500), fontsize=22, linewidth=3) + +ax1 = PL.Axis(fig[1, 1]; yscale = log10) +ax2 = PL.Axis(fig[1, 2]; yscale = log10) + +ax1.xlabel = "T [C]" +ax1.ylabel = "ice mass melting rate [g/m3/s]" +ax2.xlabel = "T [C]" +ax2.ylabel = "ice number melting rate [1/cm3/s]" + +l_max_dLdt = PL.lines!(ax1, ΔT_range, max_dLdt * 1e3, color = :thistle) +l_max_dNdt = PL.lines!(ax2, ΔT_range, max_dNdt * 1e-6, color = :thistle) + +l_dLdt1 = PL.lines!(ax1, ΔT_range, dLdt1 * 1e3, color = :skyblue) +l_dNdt1 = PL.lines!(ax2, ΔT_range, dNdt1 * 1e-6, color = :skyblue) + +l_dLdt2 = PL.lines!(ax1, ΔT_range, dLdt2 * 1e3, color = :blue3) +l_dNdt2 = PL.lines!(ax2, ΔT_range, dNdt2 * 1e-6, color = :blue3) + +l_dLdt3 = PL.lines!(ax1, ΔT_range, dLdt3 * 1e3, color = :orchid) +l_dNdt3 = PL.lines!(ax2, ΔT_range, dNdt3 * 1e-6, color = :orchid) + +l_dLdt4 = PL.lines!(ax1, ΔT_range, dLdt4 * 1e3, color = :purple) +l_dNdt4 = PL.lines!(ax2, ΔT_range, dNdt4 * 1e-6, color = :purple) + +PL.Legend( + fig[1, 3], + [l_max_dNdt, l_dNdt1, l_dNdt2, l_dNdt3, l_dNdt4], + [ + "limit", + "ρₐ=1.2 kg/m3, Fᵣ=0.8, ρᵣ=800kg/m3", + "ρₐ=1.2 kg/m3, Fᵣ=0.2, ρᵣ=800kg/m3", + "ρₐ=1.2 kg/m3, Fᵣ=0.2, ρᵣ=200kg/m3", + "ρₐ=0.5 kg/m3, Fᵣ=0.2, ρᵣ=200kg/m3", + ], + framevisible = false, +) +PL.save("P3_ice_melt.svg", fig) diff --git a/docs/src/plots/P3SchemePlots.jl b/docs/src/plots/P3SchemePlots.jl index 77eb66a05..a06c42b25 100644 --- a/docs/src/plots/P3SchemePlots.jl +++ b/docs/src/plots/P3SchemePlots.jl @@ -8,6 +8,7 @@ FT = Float64 const PSP3 = CMP.ParametersP3 p3 = CMP.ParametersP3(FT) +F_liq = FT(0) function define_axis( fig, @@ -50,24 +51,24 @@ function p3_relations_plot() # define plot axis #[row, column] ax1 = define_axis(fig, 1:7, 1:9, CMK.L"m(D) regime for $ρ_r = 400 kg m^{-3}$", CMK.L"$m$ (kg)", [1e-10, 1e-8, 1e-6], 1.9) - ax3 = define_axis(fig, 1:7, 10:18, CMK.L"m(D) regime for $F_r = 0.95$", CMK.L"$m$ (kg)", [1e-10, 1e-8, 1e-6], 1.8) + ax3 = define_axis(fig, 1:7, 10:18, CMK.L"m(D) regime for $F_rim = 0.95$", CMK.L"$m$ (kg)", [1e-10, 1e-8, 1e-6], 1.8) ax2 = define_axis(fig, 8:15, 1:9, CMK.L"A(D) regime for $ρ_r = 400 kg m^{-3}$", CMK.L"$A$ ($m^2$)", [1e-8, 1e-6, 1e-4], 1.7) - ax4 = define_axis(fig, 8:15, 10:18, CMK.L"A(D) regime for $F_r = 0.95$", CMK.L"$A$ ($m^2$)", [1e-8, 1e-6, 1e-4], 1.6) + ax4 = define_axis(fig, 8:15, 10:18, CMK.L"A(D) regime for $F_rim = 0.95$", CMK.L"$A$ ($m^2$)", [1e-8, 1e-6, 1e-4], 1.6) ax5 = define_axis(fig, 16:22, 1:9, CMK.L"ρ(D) regime for $ρ_r = 400 kg m^{-3}$", CMK.L"$ρ$ ($kg m^{-3}$)", [100, 500, 900], 2, logscale = false) - ax6 = define_axis(fig, 16:22, 10:18, CMK.L"ρ(D) regime for $F_r = 0.95$", CMK.L"$ρ$ ($kg m^{-3}$)", [100, 500, 900], 1.9, logscale = false) + ax6 = define_axis(fig, 16:22, 10:18, CMK.L"ρ(D) regime for $F_rim = 0.95$", CMK.L"$ρ$ ($kg m^{-3}$)", [100, 500, 900], 1.9, logscale = false) # Get thresholds sol4_0 = P3.thresholds(p3, 400.0, 0.0) sol4_5 = P3.thresholds(p3, 400.0, 0.5) sol4_8 = P3.thresholds(p3, 400.0, 0.8) # m(D) - fig1_0 = CMK.lines!(ax1, D_range * 1e3, [P3.p3_mass(p3, D, 0.0 ) for D in D_range], color = cl[1], linewidth = lw) - fig1_5 = CMK.lines!(ax1, D_range * 1e3, [P3.p3_mass(p3, D, 0.5, sol4_5) for D in D_range], color = cl[2], linewidth = lw) - fig1_8 = CMK.lines!(ax1, D_range * 1e3, [P3.p3_mass(p3, D, 0.8, sol4_8) for D in D_range], color = cl[3], linewidth = lw) + fig1_0 = CMK.lines!(ax1, D_range * 1e3, [P3.p3_mass(p3, D, 0.0, F_liq ) for D in D_range], color = cl[1], linewidth = lw) + fig1_5 = CMK.lines!(ax1, D_range * 1e3, [P3.p3_mass(p3, D, 0.5, F_liq, sol4_5) for D in D_range], color = cl[2], linewidth = lw) + fig1_8 = CMK.lines!(ax1, D_range * 1e3, [P3.p3_mass(p3, D, 0.8, F_liq, sol4_8) for D in D_range], color = cl[3], linewidth = lw) # a(D) - fig2_0 = CMK.lines!(ax2, D_range * 1e3, [P3.p3_area(p3, D, 0.0 ) for D in D_range], color = cl[1], linewidth = lw) - fig2_5 = CMK.lines!(ax2, D_range * 1e3, [P3.p3_area(p3, D, 0.5, sol4_5) for D in D_range], color = cl[2], linewidth = lw) - fig2_8 = CMK.lines!(ax2, D_range * 1e3, [P3.p3_area(p3, D, 0.8, sol4_8) for D in D_range], color = cl[3], linewidth = lw) + fig2_0 = CMK.lines!(ax2, D_range * 1e3, [P3.p3_area(p3, D, 0.0, F_liq ) for D in D_range], color = cl[1], linewidth = lw) + fig2_5 = CMK.lines!(ax2, D_range * 1e3, [P3.p3_area(p3, D, 0.5, F_liq, sol4_5) for D in D_range], color = cl[2], linewidth = lw) + fig2_8 = CMK.lines!(ax2, D_range * 1e3, [P3.p3_area(p3, D, 0.8, F_liq, sol4_8) for D in D_range], color = cl[3], linewidth = lw) # ρ(D) density_0 = CMK.lines!(ax5, D_range * 1e3, [P3.p3_density(p3, D, 0.0, sol4_0) for D in D_range], color = cl[1], linewidth = lw) density_5 = CMK.lines!(ax5, D_range * 1e3, [P3.p3_density(p3, D, 0.5, sol4_5) for D in D_range], color = cl[2], linewidth = lw) @@ -86,13 +87,13 @@ function p3_relations_plot() sol_4 = P3.thresholds(p3, 400.0, 0.95) sol_8 = P3.thresholds(p3, 800.0, 0.95) # m(D) - fig3_200 = CMK.lines!(ax3, D_range * 1e3, [P3.p3_mass(p3, D, 0.95, sol_2) for D in D_range], color = cl[1], linewidth = lw) - fig3_400 = CMK.lines!(ax3, D_range * 1e3, [P3.p3_mass(p3, D, 0.95, sol_4) for D in D_range], color = cl[2], linewidth = lw) - fig3_800 = CMK.lines!(ax3, D_range * 1e3, [P3.p3_mass(p3, D, 0.95, sol_8) for D in D_range], color = cl[3], linewidth = lw) + fig3_200 = CMK.lines!(ax3, D_range * 1e3, [P3.p3_mass(p3, D, 0.95, F_liq, sol_2) for D in D_range], color = cl[1], linewidth = lw) + fig3_400 = CMK.lines!(ax3, D_range * 1e3, [P3.p3_mass(p3, D, 0.95, F_liq, sol_4) for D in D_range], color = cl[2], linewidth = lw) + fig3_800 = CMK.lines!(ax3, D_range * 1e3, [P3.p3_mass(p3, D, 0.95, F_liq, sol_8) for D in D_range], color = cl[3], linewidth = lw) # a(D) - fig3_200 = CMK.lines!(ax4, D_range * 1e3, [P3.p3_area(p3, D, 0.5, sol_2) for D in D_range], color = cl[1], linewidth = lw) - fig3_400 = CMK.lines!(ax4, D_range * 1e3, [P3.p3_area(p3, D, 0.5, sol_4) for D in D_range], color = cl[2], linewidth = lw) - fig3_800 = CMK.lines!(ax4, D_range * 1e3, [P3.p3_area(p3, D, 0.5, sol_8) for D in D_range], color = cl[3], linewidth = lw) + fig3_200 = CMK.lines!(ax4, D_range * 1e3, [P3.p3_area(p3, D, 0.5, F_liq, sol_2) for D in D_range], color = cl[1], linewidth = lw) + fig3_400 = CMK.lines!(ax4, D_range * 1e3, [P3.p3_area(p3, D, 0.5, F_liq, sol_4) for D in D_range], color = cl[2], linewidth = lw) + fig3_800 = CMK.lines!(ax4, D_range * 1e3, [P3.p3_area(p3, D, 0.5, F_liq, sol_8) for D in D_range], color = cl[3], linewidth = lw) # ρ(D) density_200 = CMK.lines!(ax6, D_range * 1e3, [P3.p3_density(p3, D, 0.5, sol_2) for D in D_range], color = cl[1], linewidth = lw) density_400 = CMK.lines!(ax6, D_range * 1e3, [P3.p3_density(p3, D, 0.5, sol_4) for D in D_range], color = cl[2], linewidth = lw) @@ -108,10 +109,10 @@ function p3_relations_plot() global d_gr_800 = CMK.vlines!(ax, sol_8[2] * 1e3, linestyle = :dash, color = cl[3], linewidth = lw) end # add legend - CMK.Legend(fig[24:25, 1], [fig1_0, fig1_5, fig1_8], [CMK.L"$F_{r} = 0.0$", CMK.L"$F_{r} = 0.5$", CMK.L"$F_{r} = 0.8$"], framevisible = false) + CMK.Legend(fig[24:25, 1], [fig1_0, fig1_5, fig1_8], [CMK.L"$F_{rim} = 0.0$", CMK.L"$F_{rim} = 0.5$", CMK.L"$F_{rim} = 0.8$"], framevisible = false) CMK.Legend(fig[24:25, 3], [d_tha], [CMK.L"$D_{th}$"], framevisible = false) - CMK.Legend(fig[24:25, 7], [d_cr_5, d_cr_8], [CMK.L"$D_{cr}$ for $F_{r} = 0.5$", CMK.L"$D_{cr}$ for $F_{r} = 0.8$"], framevisible = false) - CMK.Legend(fig[24:25, 5], [d_gr_5, d_gr_8], [CMK.L"$D_{gr}$ for $F_{r} = 0.5$", CMK.L"$D_{gr}$ for $F_{r} = 0.8$"], framevisible = false) + CMK.Legend(fig[24:25, 7], [d_cr_5, d_cr_8], [CMK.L"$D_{cr}$ for $F_{rim} = 0.5$", CMK.L"$D_{cr}$ for $F_{rim} = 0.8$"], framevisible = false) + CMK.Legend(fig[24:25, 5], [d_gr_5, d_gr_8], [CMK.L"$D_{gr}$ for $F_{rim} = 0.5$", CMK.L"$D_{gr}$ for $F_{rim} = 0.8$"], framevisible = false) CMK.Legend(fig[24:25, 13], [fig3_200, fig3_400, fig3_800], [CMK.L"$\rho_{r} = 200.0 kg m^{-3}$", CMK.L"$\rho_{r} = 400.0 kg m^{-3}$", CMK.L"$\rho_{r} = 800.0 kg m^{-3}$",], framevisible = false) CMK.Legend(fig[24:25, 14], [d_thb], [CMK.L"$D_{th}$"], framevisible = false) diff --git a/docs/src/plots/P3ShapeSolverPlots.jl b/docs/src/plots/P3ShapeSolverPlots.jl index 5106ea2cc..ab81c87ab 100644 --- a/docs/src/plots/P3ShapeSolverPlots.jl +++ b/docs/src/plots/P3ShapeSolverPlots.jl @@ -8,6 +8,8 @@ FT = Float64 const PSP3 = CMP.ParametersP3 p3 = CMP.ParametersP3(FT) +F_liq = FT(0) # preserving original P3 +# TODO: investigate for F_liq != 0 function guess_value(λ::FT, p1::FT, p2::FT, L1::FT, L2::FT) return L1 * (λ / p1)^((log(L1) - log(L2)) / (log(p1) - log(p2))) @@ -16,31 +18,42 @@ end function lambda_guess_plot() N = FT(1e6) - F_r_s = [FT(0.0), FT(0.5), FT(0.8)] + F_rim_s = [FT(0.0), FT(0.5), FT(0.8)] ρ_r_s = [FT(200), FT(400), FT(800)] f = Plt.Figure() - for i in 1:length(F_r_s) + for i in 1:length(F_rim_s) for j in 1:length(ρ_r_s) - F_r = F_r_s[i] + F_rim = F_rim_s[i] ρ_r = ρ_r_s[j] Plt.Axis( f[i, j], xlabel = "log(L/N)", ylabel = "log(λ)", - title = string("λ vs L/N for F_r = ", F_r, " and ρ_r = ", ρ_r), + title = string( + "λ vs L/N for F_rim = ", + F_rim, + " and ρ_r = ", + ρ_r, + ), height = 300, width = 400, ) logλs = FT(1):FT(0.01):FT(6) λs = [10^logλ for logλ in logλs] - th = P3.thresholds(p3, ρ_r, F_r) + th = P3.thresholds(p3, ρ_r, F_rim) Ls_over_N = [ - P3.L_over_N_gamma(p3, F_r, log(λ), P3.DSD_μ(p3, λ), th) for - λ in λs + P3.L_over_N_gamma( + p3, + F_rim, + F_liq, + log(λ), + P3.DSD_μ(p3, λ), + th, + ) for λ in λs ] guesses = [FT(0) for λ in λs] @@ -48,8 +61,9 @@ function lambda_guess_plot() (min,) = P3.get_bounds( N, Ls_over_N[i] * N, - P3.DSD_μ_approx(p3, Ls_over_N[i] * N, N, ρ_r, F_r), - F_r, + P3.DSD_μ_approx(p3, Ls_over_N[i] * N, N, ρ_r, F_rim, F_liq), + F_rim, + F_liq, p3, th, ) diff --git a/docs/src/plots/P3TerminalVelocityPlots.jl b/docs/src/plots/P3TerminalVelocityPlots.jl index db7d52cfe..cad714dd4 100644 --- a/docs/src/plots/P3TerminalVelocityPlots.jl +++ b/docs/src/plots/P3TerminalVelocityPlots.jl @@ -10,30 +10,33 @@ const PSP3 = CMP.ParametersP3 FT = Float64 p3 = CMP.ParametersP3(FT) +F_liq = FT(0) # set to 0 for now function get_values( p3::PSP3, - Chen2022::CMP.Chen2022VelTypeSnowIce, + Chen2022::CMP.Chen2022VelType, L::FT, N::FT, ρ_a::FT, x_resolution::Int, y_resolution::Int, ) where {FT} - F_rs = range(FT(0), stop = FT(1 - eps(FT)), length = x_resolution) + F_rims = range(FT(0), stop = FT(1 - eps(FT)), length = x_resolution) ρ_rs = range(FT(25), stop = FT(975), length = y_resolution) V_m = zeros(x_resolution, y_resolution) + V_m_ϕ = zeros(x_resolution, y_resolution) D_m = zeros(x_resolution, y_resolution) D_m_regimes = zeros(x_resolution, y_resolution) ϕᵢ = zeros(x_resolution, y_resolution) for i in 1:x_resolution for j in 1:y_resolution - F_r = F_rs[i] + F_rim = F_rims[i] ρ_r = ρ_rs[j] - aspect_ratio = true - th = P3.thresholds(p3, ρ_r, F_r) + use_aspect_ratio = true + do_not_use_aspect = false + th = P3.thresholds(p3, ρ_r, F_rim) V_m[i, j] = P3.ice_terminal_velocity( p3, @@ -41,21 +44,34 @@ function get_values( L, N, ρ_r, - F_r, + F_rim, + F_liq, ρ_a, - aspect_ratio, + do_not_use_aspect, )[2] - D_m[i, j] = P3.D_m(p3, L, N, ρ_r, F_r) + V_m_ϕ[i, j] = P3.ice_terminal_velocity( + p3, + Chen2022, + L, + N, + ρ_r, + F_rim, + F_liq, + ρ_a, + use_aspect_ratio, + )[2] + + D_m[i, j] = P3.D_m(p3, L, N, ρ_r, F_rim, F_liq) D_m_regimes[i, j] = D_m[i, j] - ϕᵢ[i, j] = P3.ϕᵢ(p3, D_m[i, j], F_r, th) + ϕᵢ[i, j] = P3.ϕᵢ(p3, D_m[i, j], F_rim, th) # plot the regimes D_th = P3.D_th_helper(p3) if D_th > D_m[i, j] # small spherical ice D_m_regimes[i, j] = 0 - elseif F_r == 0 + elseif F_rim == 0 # large nonspherical unrimed ice D_m_regimes[i, j] = 0 elseif th.D_gr > D_m[i, j] >= D_th @@ -71,13 +87,13 @@ function get_values( end end D_m *= 1e3 - return (; F_rs, ρ_rs, D_m_regimes, D_m, ϕᵢ, V_m) + return (; F_rims, ρ_rs, D_m_regimes, D_m, ϕᵢ, V_m, V_m_ϕ) end function make_axis(fig, row, col, title) return Plt.Axis( fig[row, col], - xlabel = "F_r", + xlabel = "F_rim", ylabel = "ρ_r", title = title, width = 350, @@ -109,12 +125,12 @@ function figure_2() xres = 100 yres = 100 - (F_rs, ρ_rs, D_m_regimes_s, D_m_s, ϕᵢ_s, V_m_s) = - get_values(p3, Chen2022.snow_ice, L_s, N_s, ρ_a, xres, yres) - (F_rm, ρ_rm, D_m_regimes_m, D_m_m, ϕᵢ_m, V_m_m) = - get_values(p3, Chen2022.snow_ice, L_m, N_m, ρ_a, xres, yres) - (F_rl, ρ_rl, D_m_regimes_l, D_m_l, ϕᵢ_l, V_m_l) = - get_values(p3, Chen2022.snow_ice, L_l, N_l, ρ_a, xres, yres) + (F_rims, ρ_rs, D_m_regimes_s, D_m_s, ϕᵢ_s, V_m_s, V_m_ϕ_s) = + get_values(p3, Chen2022, L_s, N_s, ρ_a, xres, yres) + (F_rimm, ρ_rm, D_m_regimes_m, D_m_m, ϕᵢ_m, V_m_m, V_m_ϕ_m) = + get_values(p3, Chen2022, L_m, N_m, ρ_a, xres, yres) + (F_riml, ρ_rl, D_m_regimes_l, D_m_l, ϕᵢ_l, V_m_l, V_m_ϕ_l) = + get_values(p3, Chen2022, L_l, N_l, ρ_a, xres, yres) fig = Plt.Figure() @@ -131,7 +147,7 @@ function figure_2() ax1 = make_axis(fig, 1, 1, "Particle regimes with small Dₘ") hm = Plt.contourf!( ax1, - F_rs, + F_rims, ρ_rs, D_m_regimes_s, levels = 3, @@ -149,7 +165,7 @@ function figure_2() ax2 = make_axis(fig, 1, 2, "Particle regimes with medium Dₘ") hm = Plt.contourf!( ax2, - F_rm, + F_rimm, ρ_rm, D_m_regimes_m, levels = 3, @@ -167,7 +183,7 @@ function figure_2() ax3 = make_axis(fig, 1, 3, "Particle regimes with large Dₘ") hm = Plt.contourf!( ax3, - F_rl, + F_riml, ρ_rl, D_m_regimes_l, levels = 3, @@ -183,47 +199,62 @@ function figure_2() ) ax4 = make_axis(fig, 3, 1, "Dₘ (mm) (L = 8e-4 kgm⁻³, N = 1e6 m⁻³)") - hm = Plt.contourf!(ax4, F_rs, ρ_rs, D_m_s) + hm = Plt.contourf!(ax4, F_rims, ρ_rs, D_m_s) Plt.Colorbar(fig[4, 1], hm, vertical = false) ax5 = make_axis(fig, 3, 2, "Dₘ (mm) (L = 0.22 kgm⁻³, N = 1e6 m⁻³)") - hm = Plt.contourf!(ax5, F_rm, ρ_rm, D_m_m) + hm = Plt.contourf!(ax5, F_rimm, ρ_rm, D_m_m) Plt.Colorbar(fig[4, 2], hm, vertical = false) ax6 = make_axis(fig, 3, 3, "Dₘ (mm) (L = 0.7 kgm⁻³, N = 1e6 m⁻³)") - hm = Plt.contourf!(ax6, F_rl, ρ_rl, D_m_l) + hm = Plt.contourf!(ax6, F_riml, ρ_rl, D_m_l) Plt.Colorbar(fig[4, 3], hm, vertical = false) ax7 = make_axis(fig, 5, 1, "ϕᵢ with small Dₘ") - hm = Plt.contourf!(ax7, F_rs, ρ_rs, ϕᵢ_s, levels = 20) + hm = Plt.contourf!(ax7, F_rims, ρ_rs, ϕᵢ_s, levels = 20) Plt.Colorbar(fig[6, 1], hm, vertical = false) - Plt.contour!(ax7, F_rs, ρ_rs, D_m_s, label = "D_m"; contour_args...) + Plt.contour!(ax7, F_rims, ρ_rs, D_m_s, label = "D_m"; contour_args...) ax8 = make_axis(fig, 5, 2, "ϕᵢ with medium Dₘ") - hm = Plt.contourf!(ax8, F_rm, ρ_rm, ϕᵢ_m) + hm = Plt.contourf!(ax8, F_rimm, ρ_rm, ϕᵢ_m) Plt.Colorbar(fig[6, 2], hm, vertical = false) - Plt.contour!(ax8, F_rm, ρ_rm, D_m_m, label = "D_m"; contour_args...) + Plt.contour!(ax8, F_rimm, ρ_rm, D_m_m, label = "D_m"; contour_args...) ax9 = make_axis(fig, 5, 3, "ϕᵢ with large Dₘ") - hm = Plt.contourf!(ax9, F_rl, ρ_rl, ϕᵢ_l) + hm = Plt.contourf!(ax9, F_riml, ρ_rl, ϕᵢ_l) Plt.Colorbar(fig[6, 3], hm, vertical = false) - Plt.contour!(ax9, F_rl, ρ_rl, D_m_l, label = "D_m"; contour_args...) + Plt.contour!(ax9, F_riml, ρ_rl, D_m_l, label = "D_m"; contour_args...) - ax10 = make_axis(fig, 7, 1, "Vₘ with small Dₘ") - hm = Plt.contourf!(ax10, F_rs, ρ_rs, V_m_s) - Plt.contour!(ax10, F_rs, ρ_rs, D_m_s; contour_args...) + ax10 = make_axis(fig, 7, 1, "Vₘ (ϕᵢ = 1) with small Dₘ") + hm = Plt.contourf!(ax10, F_rims, ρ_rs, V_m_s) + Plt.contour!(ax10, F_rims, ρ_rs, D_m_s; contour_args...) Plt.Colorbar(fig[8, 1], hm, vertical = false) - ax11 = make_axis(fig, 7, 2, "Vₘ with medium Dₘ") - hm = Plt.contourf!(ax11, F_rm, ρ_rm, V_m_m) - Plt.contour!(ax11, F_rm, ρ_rm, D_m_m; contour_args...) + ax11 = make_axis(fig, 7, 2, "Vₘ (ϕᵢ = 1) with medium Dₘ") + hm = Plt.contourf!(ax11, F_rimm, ρ_rm, V_m_m) + Plt.contour!(ax11, F_rimm, ρ_rm, D_m_m; contour_args...) Plt.Colorbar(fig[8, 2], hm, vertical = false) - ax12 = make_axis(fig, 7, 3, "Vₘ with large Dₘ") - hm = Plt.contourf!(ax12, F_rl, ρ_rl, V_m_l) - Plt.contour!(ax12, F_rl, ρ_rl, D_m_l; contour_args...) + ax12 = make_axis(fig, 7, 3, "Vₘ (ϕᵢ = 1) with large Dₘ") + hm = Plt.contourf!(ax12, F_riml, ρ_rl, V_m_l) + Plt.contour!(ax12, F_riml, ρ_rl, D_m_l; contour_args...) Plt.Colorbar(fig[8, 3], hm, vertical = false) + ax13 = make_axis(fig, 9, 1, "Vₘ (using ϕᵢ) with small Dₘ") + hm = Plt.contourf!(ax13, F_rims, ρ_rs, V_m_ϕ_s) + Plt.contour!(ax13, F_rims, ρ_rs, D_m_s; contour_args...) + Plt.Colorbar(fig[10, 1], hm, vertical = false) + + ax14 = make_axis(fig, 9, 2, "Vₘ (using ϕᵢ) with medium Dₘ") + hm = Plt.contourf!(ax14, F_rimm, ρ_rm, V_m_ϕ_m) + Plt.contour!(ax14, F_rimm, ρ_rm, D_m_m; contour_args...) + Plt.Colorbar(fig[10, 2], hm, vertical = false) + + ax15 = make_axis(fig, 9, 3, "Vₘ (using ϕᵢ) with large Dₘ") + hm = Plt.contourf!(ax15, F_riml, ρ_rl, V_m_ϕ_l) + Plt.contour!(ax15, F_riml, ρ_rl, D_m_l; contour_args...) + Plt.Colorbar(fig[10, 3], hm, vertical = false) + Plt.linkxaxes!(ax1, ax2) Plt.linkxaxes!(ax2, ax3) Plt.linkyaxes!(ax1, ax2) @@ -246,6 +277,12 @@ function figure_2() Plt.linkyaxes!(ax1, ax10) Plt.linkyaxes!(ax10, ax11) Plt.linkyaxes!(ax11, ax12) + Plt.linkxaxes!(ax1, ax13) + Plt.linkxaxes!(ax13, ax14) + Plt.linkxaxes!(ax14, ax15) + Plt.linkyaxes!(ax1, ax13) + Plt.linkyaxes!(ax13, ax14) + Plt.linkyaxes!(ax14, ax15) Plt.resize_to_layout!(fig) Plt.save("MorrisonandMilbrandtFig2.svg", fig) diff --git a/docs/src/plots/P3TerminalVelocityPlots_WithFliq.jl b/docs/src/plots/P3TerminalVelocityPlots_WithFliq.jl new file mode 100644 index 000000000..2a7ebc619 --- /dev/null +++ b/docs/src/plots/P3TerminalVelocityPlots_WithFliq.jl @@ -0,0 +1,190 @@ +import ClimaParams +import CloudMicrophysics as CM +import CloudMicrophysics.P3Scheme as P3 +import CloudMicrophysics.Parameters as CMP +import CairoMakie as Plt + +const PSP3 = CMP.ParametersP3 + +FT = Float64 + +p3 = CMP.ParametersP3(FT) + +# Testing terminal velocity with liquid fraction + +function get_values( + p3::PSP3, + Chen2022::CMP.Chen2022VelType, + L::FT, + N::FT, + F_liq::FT, + ρ_a::FT, + x_resolution::Int, + y_resolution::Int, +) where {FT <: Real} + F_rims = range(FT(0), stop = FT(1 - eps(FT)), length = x_resolution) + ρ_rs = range(FT(25), stop = FT(975), length = y_resolution) + + V_m = zeros(x_resolution, y_resolution) + D_m = zeros(x_resolution, y_resolution) + use_aspect_ratio = true + + for i in 1:x_resolution + for j in 1:y_resolution + F_rim = F_rims[i] + ρ_r = ρ_rs[j] + + V_m[i, j] = P3.ice_terminal_velocity( + p3, + Chen2022, + L, + N, + ρ_r, + F_rim, + F_liq, + ρ_a, + use_aspect_ratio, + )[2] + # get D_m in mm for plots + D_m[i, j] = 1e3 * P3.D_m(p3, L, N, ρ_r, F_rim, F_liq) + end + end + return (; F_rims, ρ_rs, V_m, D_m) +end + +function make_axis(fig, row, col, title) + return Plt.Axis( + fig[row, col], + xlabel = "F_rim", + ylabel = "ρ_r", + title = title, + width = 350, + height = 350, + limits = (0, 1.0, 25, 975), + xticks = [0.0, 0.2, 0.4, 0.6, 0.8, 1.0], + yticks = [200, 400, 600, 800], + ) +end +#! format: off +function figure_2() + + Chen2022 = CMP.Chen2022VelType(FT) + # density of air in kg/m^3 + ρ_a = FT(1.2) #FT(1.293) + # small D_m + L_s_0 = FT(0.0008) + L_s_5 = FT(0.00091) + L_s_9 = FT(0.00096) + N_s = FT(1e6) + # medium D_m + L_m_0 = FT(0.22) + L_m_5 = FT(0.555) + L_m_9 = FT(0.635) + N_m = FT(1e6) + # large D_m + L_l_0 = FT(0.7) + L_l_5 = FT(2.6) + L_l_9 = FT(3) + N_l = FT(1e6) + # get V_m and D_m + xres = 100 + yres = 100 + + (F_rims_0, ρ_rs_0, V_ms_0, D_ms_0) = get_values(p3, Chen2022, L_s_0, N_s, FT(0), ρ_a, xres, yres) + (F_rimm_0, ρ_rm_0, V_mm_0, D_mm_0) = get_values(p3, Chen2022, L_m_0, N_m, FT(0), ρ_a, xres, yres) + (F_riml_0, ρ_rl_0, V_ml_0, D_ml_0) = get_values(p3, Chen2022, L_l_0, N_l, FT(0), ρ_a, xres, yres) + + (F_rims_5, ρ_rs_5, V_ms_5, D_ms_5) = get_values(p3, Chen2022, L_s_5, N_s, FT(0.5), ρ_a, xres, yres) + (F_rimm_5, ρ_rm_5, V_mm_5, D_mm_5) = get_values(p3, Chen2022, L_m_5, N_m, FT(0.5), ρ_a, xres, yres) + (F_riml_5, ρ_rl_5, V_ml_5, D_ml_5) = get_values(p3, Chen2022, L_l_5, N_l, FT(0.5), ρ_a, xres, yres) + + (F_rims_9, ρ_rs_9, V_ms_9, D_ms_9) = get_values(p3, Chen2022, L_s_9, N_s, FT(0.9), ρ_a, xres, yres) + (F_rimm_9, ρ_rm_9, V_mm_9, D_mm_9) = get_values(p3, Chen2022, L_m_9, N_m, FT(0.9), ρ_a, xres, yres) + (F_riml_9, ρ_rl_9, V_ml_9, D_ml_9) = get_values(p3, Chen2022, L_l_9, N_l, FT(0.9), ρ_a, xres, yres) + + fig = Plt.Figure() + + # Plot velocities as in Fig 2 in Morrison and Milbrandt 2015 + + args = (color = :black, labels = true, levels = 4, linewidth = 1.5, labelsize = 18) + + # set the colorscale to be the same for all figures + crange_small(len) = range(0.39, 0.62, length = len) + crange_med(len) = range(3, 7.5, length = len) + crange_large(len) = range(5, 10, length = len) + + ax1 = make_axis(fig, 1, 1, "Vₘ with small Dₘ, F_liq = 0") + hm = Plt.contourf!(ax1, F_rims_0, ρ_rs_0, V_ms_0, levels = crange_small(10), extendlow = :auto, extendhigh = :auto) + Plt.contour!(ax1, F_rims_0, ρ_rs_0, D_ms_0; args...) + # Plt.Colorbar(fig[2, 1], hm, vertical = false) + + ax2 = make_axis(fig, 1, 2, "Vₘ with medium Dₘ, F_liq = 0") + hm = Plt.contourf!(ax2, F_rimm_0, ρ_rm_0, V_mm_0, levels = crange_med(10), extendlow = :auto, extendhigh = :auto) + Plt.contour!(ax2, F_rimm_0, ρ_rm_0, D_mm_0; args...) + # Plt.Colorbar(fig[2, 2], hm, vertical = false) + + ax3 = make_axis(fig, 1, 3, "Vₘ with large Dₘ, F_liq = 0") + hm = Plt.contourf!(ax3, F_riml_0, ρ_rl_0, V_ml_0, levels = crange_large(10), extendlow = :auto, extendhigh = :auto) + Plt.contour!(ax3, F_riml_0, ρ_rl_0, D_ml_0; args...) + # Plt.Colorbar(fig[2, 3], hm, vertical = false) + + Plt.linkxaxes!(ax1, ax2) + Plt.linkxaxes!(ax2, ax3) + Plt.linkyaxes!(ax1, ax2) + Plt.linkyaxes!(ax2, ax3) + + ## Plot F_liq = 0.5 as second row of comparisons + + ax4 = make_axis(fig, 3, 1, "Vₘ with small Dₘ, F_liq = 0.5") + hm = Plt.contourf!(ax4, F_rims_5, ρ_rs_5, V_ms_5, levels = crange_small(25), extendlow = :auto, extendhigh = :auto) + Plt.contour!(ax4, F_rims_5, ρ_rs_5, D_ms_5; args...) + # Plt.Colorbar(fig[4, 1], hm, vertical = false) + + ax5 = make_axis(fig, 3, 2, "Vₘ with medium Dₘ, F_liq = 0.5") + hm = Plt.contourf!(ax5, F_rimm_5, ρ_rm_5, V_mm_5, levels = crange_med(25), extendlow = :auto, extendhigh = :auto) + Plt.contour!(ax5, F_rimm_5, ρ_rm_5, D_mm_5; args...) + # Plt.Colorbar(fig[4, 2], hm, vertical = false) + + ax6 = make_axis(fig, 3, 3, "Vₘ with large Dₘ, F_liq = 0.5") + hm = Plt.contourf!(ax6, F_riml_5, ρ_rl_5, V_ml_5, levels = crange_large(30), extendlow = :auto, extendhigh = :auto) + Plt.contour!(ax6, F_riml_5, ρ_rl_5, D_ml_5; args...) + # Plt.Colorbar(fig[4, 3], hm, vertical = false) + + Plt.linkxaxes!(ax1, ax4) + Plt.linkxaxes!(ax4, ax5) + Plt.linkxaxes!(ax5, ax6) + Plt.linkyaxes!(ax1, ax4) + Plt.linkyaxes!(ax4, ax5) + Plt.linkyaxes!(ax5, ax6) + + ## Plot F_liq = 0.9 as second row of comparisons + + ax7 = make_axis(fig, 5, 1, "Vₘ with small Dₘ, F_liq = 0.9") + hm = Plt.contourf!(ax7, F_rims_9, ρ_rs_9, V_ms_9, levels = crange_small(45), extendlow = :auto, extendhigh = :auto) + Plt.contour!(ax7, F_rims_9, ρ_rs_9, D_ms_9; args...) + Plt.Colorbar(fig[6, 1], hm, vertical = false) + + ax8 = make_axis(fig, 5, 2, "Vₘ with medium Dₘ, F_liq = 0.9") + hm = Plt.contourf!(ax8, F_rimm_9, ρ_rm_9, V_mm_9, levels = crange_med(45), extendlow = :auto, extendhigh = :auto) + Plt.contour!(ax8, F_rimm_9, ρ_rm_9, D_mm_9; args...) + Plt.Colorbar(fig[6, 2], hm, vertical = false) + + ax9 = make_axis(fig, 5, 3, "Vₘ with large Dₘ, F_liq = 0.9") + hm = Plt.contourf!(ax9, F_riml_9, ρ_rl_9, V_ml_9, levels = crange_large(45), extendlow = :auto, extendhigh = :auto) + Plt.contour!(ax9, F_riml_9, ρ_rl_9, D_ml_9; args...) + Plt.Colorbar(fig[6, 3], hm, vertical = false) + + Plt.linkxaxes!(ax1, ax7) + Plt.linkxaxes!(ax7, ax8) + Plt.linkxaxes!(ax8, ax9) + Plt.linkyaxes!(ax1, ax7) + Plt.linkyaxes!(ax7, ax8) + Plt.linkyaxes!(ax8, ax9) + + + Plt.resize_to_layout!(fig) + # Plt.display(fig) + Plt.save("MorrisonandMilbrandtFig2_with_F_liq.svg", fig) +end + +figure_2() diff --git a/docs/src/plots/P3TerminalVelocity_F_liq_rim.jl b/docs/src/plots/P3TerminalVelocity_F_liq_rim.jl new file mode 100644 index 000000000..55f2d62d8 --- /dev/null +++ b/docs/src/plots/P3TerminalVelocity_F_liq_rim.jl @@ -0,0 +1,234 @@ +import ClimaParams +import CloudMicrophysics as CM +import CloudMicrophysics.P3Scheme as P3 +import CloudMicrophysics.Parameters as CMP +import CairoMakie as Plt + +const PSP3 = CMP.ParametersP3 + +FT = Float64 + +p3 = CMP.ParametersP3(FT) + +# Testing terminal velocity with liquid fraction + +function get_values( + p3::PSP3, + Chen2022::CMP.Chen2022VelType, + L::FT, + N::FT, + ρ_a::FT, + x_resolution::Int, + y_resolution::Int, +) where {FT <: Real} + F_rims = range(FT(0), stop = FT(1 - eps(FT)), length = x_resolution) + F_liqs = range(FT(0), stop = FT(1), length = y_resolution) + ρ_r = FT(900) + use_aspect_ratio = true + + V_m = zeros(x_resolution, y_resolution) + D_m_regimes = zeros(x_resolution, y_resolution) + D_m = zeros(x_resolution, y_resolution) + + for i in 1:x_resolution + for j in 1:y_resolution + F_rim = F_rims[i] + F_liq = F_liqs[j] + th = P3.thresholds(p3, ρ_r, F_rim) + + V_m[i, j] = P3.ice_terminal_velocity( + p3, + Chen2022, + L, + N, + ρ_r, + F_rim, + F_liq, + ρ_a, + use_aspect_ratio, + )[2] + D_m[i, j] = P3.D_m(p3, L, N, ρ_r, F_rim, F_liq) + D_m_regimes[i, j] = D_m[i, j] + + # plot the regimes + D_th = P3.D_th_helper(p3) + if D_th > D_m[i, j] + # small spherical ice + D_m_regimes[i, j] = 0 + elseif F_rim == 0 + # large nonspherical unrimed ice + D_m_regimes[i, j] = 0 + elseif th.D_gr > D_m[i, j] >= D_th + # dense nonspherical ice + D_m_regimes[i, j] = 1 + elseif th.D_cr > D_m[i, j] >= th.D_gr + # graupel + D_m_regimes[i, j] = 2 + else #elseif D >= th.D_cr + # partially rimed ice + D_m_regimes[i, j] = 3 + end + end + end + D_m *= 1e3 + return (; F_rims, F_liqs, V_m, D_m, D_m_regimes) +end + +function make_axis(fig, row, col, title) + return Plt.Axis( + fig[row, col], + xlabel = "F_rim", + ylabel = "F_liq", + title = title, + width = 350, + height = 350, + limits = (0, 1.0, 0, 1.0), + xticks = [0.0, 0.2, 0.4, 0.6, 0.8, 1.0], + yticks = [0.0, 0.2, 0.4, 0.6, 0.8, 1.0], + ) +end + +function fig1() + Chen2022 = CMP.Chen2022VelType(FT) + ρ_a = FT(1.2) + + # small D_m + L_s = FT(0.0008) + N_s = FT(1e6) + # medium D_m + L_m = FT(0.22) + N_m = FT(1e6) + # large D_m + L_l = FT(0.7) + N_l = FT(1e6) + + crange_small = range(0.39, 0.6, 20) + crange_med = range(3, 6, 20) + crange_large = range(4, 7, 20) + + # get V_m and D_m + xres = 100 + yres = 100 + + (F_rims, F_liqs, V_ms, D_ms, D_ms_regimes) = + get_values(p3, Chen2022, L_s, N_s, ρ_a, xres, yres) + (F_rimm, F_liqm, V_mm, D_mm, D_mm_regimes) = + get_values(p3, Chen2022, L_m, N_m, ρ_a, xres, yres) + (F_riml, F_liql, V_ml, D_ml, D_ml_regimes) = + get_values(p3, Chen2022, L_l, N_l, ρ_a, xres, yres) + + args = ( + color = :black, + labels = true, + levels = 3, + linewidth = 1.5, + labelsize = 18, + ) + + + fig = Plt.Figure() + + ax1_ = make_axis(fig, 1, 1, "Particle regimes with small Dₘ") + hm = Plt.contourf!( + ax1_, + F_rims, + F_liqs, + D_ms_regimes, + levels = 3, + colormap = Plt.cgrad(:PuBuGn_3, 3, categorical = true), + ) + Plt.Colorbar( + fig[2, 1], + hm, + vertical = false, + label = "graupel (2), partially rimed ice (3)", + ticks = [1, 2, 3], + labelsize = 14, + ) + + ax2_ = make_axis(fig, 1, 2, "Particle regimes with medium Dₘ") + hm = Plt.contourf!( + ax2_, + F_rimm, + F_liqm, + D_mm_regimes, + levels = 3, + colormap = Plt.cgrad(:PuBuGn_3, 3, categorical = true), + ) + Plt.Colorbar( + fig[2, 2], + hm, + vertical = false, + label = "graupel (2), partially rimed ice (3)", + ticks = [1, 2, 3], + labelsize = 14, + ) + + ax3_ = make_axis(fig, 1, 3, "Particle regimes with large Dₘ") + hm = Plt.contourf!( + ax3_, + F_riml, + F_liql, + D_ml_regimes, + levels = 3, + colormap = Plt.cgrad(:PuBuGn_3, 3, categorical = true), + ) + Plt.Colorbar( + fig[2, 3], + hm, + vertical = false, + label = "graupel (2), partially rimed ice (3)", + ticks = [1, 2, 3], + labelsize = 14, + ) + + ax1 = make_axis(fig, 3, 1, "Vₘ with small Dₘ") + hm = Plt.contourf!( + ax1, + F_rims, + F_liqs, + V_ms, + levels = crange_small, + extendlow = :auto, + extendhigh = :auto, + ) + Plt.Colorbar(fig[4, 1], hm, vertical = false) + Plt.contour!(ax1, F_rims, F_liqs, D_ms; args...) + + + ax2 = make_axis(fig, 3, 2, "Vₘ with medium Dₘ") + hm = Plt.contourf!( + ax2, + F_rimm, + F_liqm, + V_mm, + levels = crange_med, + extendlow = :auto, + extendhigh = :auto, + ) + Plt.Colorbar(fig[4, 2], hm, vertical = false) + Plt.contour!(ax2, F_rimm, F_liqm, D_mm; args...) + + + ax3 = make_axis(fig, 3, 3, "Vₘ with large Dₘ") + hm = Plt.contourf!( + ax3, + F_riml, + F_liql, + V_ml, + levels = crange_large, + extendlow = :auto, + extendhigh = :auto, + ) + Plt.Colorbar(fig[4, 3], hm, vertical = false) + Plt.contour!(ax3, F_riml, F_liql, D_ml; args...) + + + # Plt.linkaxes!(ax1_, ax2_, ax3_) + Plt.linkaxes!(ax1, ax2, ax3) + + Plt.resize_to_layout!(fig) + Plt.save("P3TerminalVelocity_F_liq_rim.svg", fig) +end + +fig1() diff --git a/p3_sandbox/p3_sandbox.jl b/p3_sandbox/p3_sandbox.jl index 51f29a797..4898cd365 100644 --- a/p3_sandbox/p3_sandbox.jl +++ b/p3_sandbox/p3_sandbox.jl @@ -23,7 +23,7 @@ function p3_sandbox(dY, Y, p, t) qᵣ = Y[3] # Rime mass mixing ratio Bᵣ = Y[4] # Rime volume - F_r = qᵣ / qᵢ + F_rim = qᵣ / qᵢ ρ_r = ifelse(Bᵣ == 0, FT(0), qᵣ / Bᵣ) q = TD.PhasePartition(qᵥ + qₗ + qᵢ, qₗ, qᵢ) @@ -38,7 +38,7 @@ function p3_sandbox(dY, Y, p, t) dNᵢ_dt = J_immersion * Nₗ * 4 * π * rₗ^2 println("J = ", J_immersion) - sol = P3.thresholds(p3, ρ_r, F_r) + sol = P3.thresholds(p3, ρ_r, F_rim) println(" ") println("D_cr = ", sol.D_cr) println("D_gr = ", sol.D_gr) diff --git a/src/MicrophysicsNonEq.jl b/src/MicrophysicsNonEq.jl index 24fe82f92..4e95ebfab 100644 --- a/src/MicrophysicsNonEq.jl +++ b/src/MicrophysicsNonEq.jl @@ -14,6 +14,7 @@ import ..Parameters as CMP export τ_relax export conv_q_vap_to_q_liq_ice export conv_q_vap_to_q_liq_ice_MM2015 +export conv_q_vap_to_q_liq_ice_MM2015_timeintegrator """ τ_relax(liquid) @@ -91,7 +92,7 @@ function conv_q_vap_to_q_liq_ice_MM2015( dqsldT = qᵥ_sat_liq * (Lᵥ / (Rᵥ * T^2) - 1 / T) Γₗ = FT(1) + (Lᵥ / cₚ_air) * dqsldT - return (qᵥ - qᵥ_sat_liq) / τ_relax * Γₗ + return (qᵥ - qᵥ_sat_liq) / (τ_relax * Γₗ) end function conv_q_vap_to_q_liq_ice_MM2015( (; τ_relax)::CMP.CloudIce{FT}, @@ -111,7 +112,95 @@ function conv_q_vap_to_q_liq_ice_MM2015( dqsidT = qᵥ_sat_ice * (Lₛ / (Rᵥ * T^2) - 1 / T) Γᵢ = FT(1) + (Lₛ / cₚ_air) * dqsidT - return (qᵥ - qᵥ_sat_ice) / τ_relax * Γᵢ + return (qᵥ - qᵥ_sat_ice) / (τ_relax * Γᵢ) end -end #module MicrophysicsNonEq.jl +""" + conv_q_vap_to_q_liq_ice_MM2015_timeintegrator(liquid, ice, tps, q, ρ, T, w, p_air, const_dt, type) + +- `liquid` - a struct with cloud water free parameters +- `ice` - a struct with cloud ice free parameters +- `tps` - thermodynamics parameters struct +- `q` - current PhasePartition +- `ρ` - air density [kg/m3] +- `T` - air temperature [K] +- `w` - vertical velocity [m/s] +- `p_air` - air pressure [Pa] +- `const_dt` - length of time step [s] +- `"condensation"` or `"deposition"` - type of process to calculate and output [str] + +Returns the cloud water tendency due to condensation and evaporation +or cloud ice tendency due to sublimation and vapor deposition. +The formulation is based on Morrison and Grabowski 2008 and +Morrison and Milbrandt 2015 +""" +function conv_q_vap_to_q_liq_ice_MM2015_timeintegrator( + liquid::CMP.CloudLiquid{FT}, + ice::CMP.CloudIce{FT}, + tps::TDP.ThermodynamicsParameters{FT}, + q::TD.PhasePartition{FT}, + ρ::FT, + T::FT, + w::FT, + p_air::FT, + const_dt::FT, + ::Val{type}, +) where {FT, type} + + cₚ_air = TD.cp_m(tps, q) + Lₛ = TD.latent_heat_sublim(tps, T) + Lᵥ = TD.latent_heat_vapor(tps, T) + Rᵥ = TD.Parameters.R_v(tps) + g = TD.Parameters.grav(tps) + qᵥ = TD.vapor_specific_humidity(q) + + pᵥ_sat_liq = TD.saturation_vapor_pressure(tps, T, TD.Liquid()) + qᵥ_sat_liq = TD.q_vap_saturation_from_density(tps, T, ρ, pᵥ_sat_liq) + + pᵥ_sat_ice = TD.saturation_vapor_pressure(tps, T, TD.Ice()) + qᵥ_sat_ice = TD.q_vap_saturation_from_density(tps, T, ρ, pᵥ_sat_ice) + + dqsldT = qᵥ_sat_liq * (Lᵥ / (Rᵥ * T^2) - 1 / T) + dqsidT = qᵥ_sat_ice * (Lₛ / (Rᵥ * T^2) - 1 / T) + + Γₗ = FT(1) + (Lᵥ / cₚ_air) * dqsldT + Γᵢ = FT(1) + (Lₛ / cₚ_air) * dqsidT + + A_c_WBF = + (qᵥ_sat_liq - qᵥ_sat_ice) / (ice.τ_relax * Γᵢ) * + (1 + (Lₛ / cₚ_air) * dqsldT) + + A_c_uplift = + -(qᵥ_sat_liq * g * w * ρ) / (p_air - pᵥ_sat_liq) + + dqsldT * w * g / cₚ_air + + A_c = A_c_uplift - A_c_WBF + + τ = + ( + liquid.τ_relax^(-1) + + (1 + (Lₛ / cₚ_air) * dqsldT) * ice.τ_relax^(-1) / Γᵢ + )^(-1) + + δ_0 = qᵥ - qᵥ_sat_liq + + if type == :condensation + cond_rate = + A_c * τ / (liquid.τ_relax * Γₗ) + + (δ_0 - A_c * τ) * τ / (const_dt * liquid.τ_relax * Γₗ) * + (FT(1) - exp(-const_dt / τ)) + return cond_rate + elseif type == :deposition + dep_rate = + A_c * τ / (ice.τ_relax * Γᵢ) + + (δ_0 - A_c * τ) * τ / (const_dt * ice.τ_relax * Γᵢ) * + (FT(1) - exp(-const_dt / τ)) + + (qᵥ_sat_liq - qᵥ_sat_ice) / (ice.τ_relax * Γᵢ) + return dep_rate + else + error("please specify condensation or deposition") + end + +end + +end diff --git a/src/P3_helpers.jl b/src/P3_helpers.jl index fc91733b8..f5a560d38 100644 --- a/src/P3_helpers.jl +++ b/src/P3_helpers.jl @@ -1,8 +1,99 @@ -# Some wrappers to cast types from SF.gamma +# Wrapper to cast types from SF.gamma # (which returns Float64 even when the input is Float32) -Γ(a::FT, z::FT) where {FT <: Real} = FT(SF.gamma(a, z)) +# TODO - replace with parameterization of our own Γ(a::FT) where {FT <: Real} = FT(SF.gamma(a)) +# helper functions for Γ_approx +function c₁(a::FT) where {FT <: Real} + p1 = FT(9.4368392235e-03) + p2 = -FT(1.0782666481e-04) + p3 = -FT(5.8969657295e-06) + p4 = FT(2.8939523781e-07) + p5 = FT(1.0043326298e-01) + p6 = FT(5.5637848465e-01) + p = [p1, p2, p3, p4, p5, p6] + + ret = 1 + p[5] * (exp(-p[6] * a) - 1) + for it in range(1, 4, step = 1) + ret += p[it] * a^it + end + return ret +end +function c₂(a::FT) where {FT <: Real} + q1 = FT(1.1464706419e-01) + q2 = FT(2.6963429121) + q3 = -FT(2.9647038257) + q4 = FT(2.1080724954) + q = [q1, q2, q3, q4] + + ret = q[1] + for it in range(2, 4, step = 1) + ret += q[it] / a^(it - 1) + end + return ret +end +function c₃(a::FT) where {FT <: Real} + r1 = FT(0) + r2 = FT(1.1428716184) + r3 = -FT(6.6981186438e-03) + r4 = FT(1.0480765092e-04) + r = [r1, r2, r3, r4] + + ret = 0 + for it in range(1, 4, step = 1) + ret += r[it] * a^(it - 1) + end + return ret +end +function c₄(a::FT) where {FT <: Real} + s1 = FT(1.0356711153) + s2 = FT(2.3423452308) + s3 = -FT(3.6174503174e-01) + s4 = -FT(3.1376557650) + s5 = FT(2.9092306039) + s = [s1, s2, s3, s4, s5] + + ret = 0 + for it in range(1, 5, step = 1) + ret += s[it] / a^(it - 1) + end + return ret +end + +""" + Γ_lower(a, z) + +An approximated lower incomplete Gamma function based on Blahak 2010: +doi:10.5194/gmd-3-329-2010 +https://gmd.copernicus.org/articles/3/329/2010/gmd-3-329-2010.pdf +""" +function Γ_lower(a::FT, z::FT) where {FT <: Real} + + if isnan(z) || z <= FT(0) + return FT(0) + else + W(a, z) = FT(0.5) + FT(0.5) * tanh(c₂(a) * (z - c₃(a))) + + return exp(-z) * + z^a * + ( + 1 / a + + c₁(a) * z / a / (a + 1) + + (c₁(a) * z)^2 / a / (a + 1) / (a + 2) + ) * + (1 - W(a, z)) + Γ(a) * W(a, z) * (1 - c₄(a)^(-z)) + end +end + +""" + Γ_upper(a, z) + +An approximated upper incomplete Gamma function computed as Γ(a) - Γ_lower(a, z) +""" +function Γ_upper(a::FT, z::FT) where {FT <: Real} + return Γ(a) - Γ_lower(a, z) +end + """ ∫_Γ(x₀, x_end, c1, c2, c3) @@ -15,12 +106,14 @@ f(D, c1, c2, c3) = c1 * D ^ (c2) * exp(-c3 * D) Integrates f(D, c1, c2, c3) dD from x₀ to x_end """ function ∫_Γ(x₀::FT, x_end::FT, c1::FT, c2::FT, c3::FT) where {FT} - if x_end == Inf - return c1 * c3^(-c2 - 1) * (Γ(1 + c2, x₀ * c3)) + if x_end == FT(Inf) + return c1 * c3^(-c2 - 1) * (Γ_upper(1 + c2, x₀ * c3)) elseif x₀ == 0 - return c1 * c3^(-c2 - 1) * (Γ(1 + c2) - Γ(1 + c2, x_end * c3)) + return c1 * c3^(-c2 - 1) * (Γ_lower(1 + c2, x_end * c3)) else - return c1 * c3^(-c2 - 1) * (Γ(1 + c2, x₀ * c3) - Γ(1 + c2, x_end * c3)) + return c1 * + c3^(-c2 - 1) * + (Γ_upper(1 + c2, x₀ * c3) - Γ_upper(1 + c2, x_end * c3)) end end diff --git a/src/P3_particle_properties.jl b/src/P3_particle_properties.jl index 9401301a6..2cf61001c 100644 --- a/src/P3_particle_properties.jl +++ b/src/P3_particle_properties.jl @@ -24,18 +24,18 @@ D_th_helper(p3::PSP3{FT}) where {FT} = (FT(π) * p3.ρ_i / 6 / α_va_si(p3))^(1 / (p3.β_va - 3)) """ - D_cr_helper(p3, F_r, ρ_g) + D_cr_helper(p3, F_rim, ρ_g) - p3 - a struct with P3 scheme parameters - - F_r - rime mass fraction (L_rim/L_ice) [-] + - F_rim - rime mass fraction L_rim / L_ice) [-] - ρ_g - is the effective density of a spherical graupel particle [kg/m^3] Returns the size of equal mass for graupel and partially rimed ice, in meters. Eq. 14 in Morrison and Milbrandt (2015). """ -function D_cr_helper(p3::PSP3{FT}, F_r::FT, ρ_g::FT) where {FT} +function D_cr_helper(p3::PSP3{FT}, F_rim::FT, ρ_g::FT) where {FT} α_va = α_va_si(p3) - return (1 / (1 - F_r) * 6 * α_va / FT(π) / ρ_g)^(1 / (3 - p3.β_va)) + return (1 / (1 - F_rim) * 6 * α_va / FT(π) / ρ_g)^(1 / (3 - p3.β_va)) end """ @@ -53,16 +53,17 @@ function D_gr_helper(p3::PSP3{FT}, ρ_g::FT) where {FT} end """ - ρ_g_helper(ρ_r, F_r, ρ_d) + ρ_g_helper(ρ_r, F_rim, ρ_d) - ρ_r - rime density (L_rim/B_rim) [kg/m^3] - - F_r - rime mass fraction (L_rim/L_ice) [-] + - F_rim - rime mass fraction (L_rim / L_ice) [-] - ρ_g - is the effective density of a spherical graupel particle [kg/m^3] Returns the density of total (deposition + rime) ice mass for graupel, in kg/m3 Eq. 16 in Morrison and Milbrandt (2015). """ -ρ_g_helper(ρ_r::FT, F_r::FT, ρ_d::FT) where {FT} = F_r * ρ_r + (1 - F_r) * ρ_d +ρ_g_helper(ρ_r::FT, F_rim::FT, ρ_d::FT) where {FT} = + F_rim * ρ_r + (1 - F_rim) * ρ_d """ ρ_d_helper(p3, D_cr, D_gr) @@ -82,11 +83,11 @@ function ρ_d_helper(p3::PSP3{FT}, D_cr::FT, D_gr::FT) where {FT} end """ - thresholds(p3, ρ_r, F_r) + thresholds(p3, ρ_r, F_rim) - p3 - a struct with P3 scheme parameters - ρ_r - rime density (L_rim/B_rim) [kg/m^3] - - F_r - rime mass fraction (L_rim/L_ice) [-] + - F_rim - rime mass fraction (L_rim / L_ice) [-] Solves the nonlinear system consisting of D_cr, D_gr, ρ_g, ρ_d for a given rime density and rime mass fraction. @@ -96,12 +97,12 @@ Returns a named tuple containing: - ρ_g - is the effective density of a spherical graupel particle [kg/m3], - ρ_d - is the density of the unrimed portion of the particle [kg/m3], """ -function thresholds(p3::PSP3{FT}, ρ_r::FT, F_r::FT) where {FT} +function thresholds(p3::PSP3{FT}, ρ_r::FT, F_rim::FT) where {FT} - @assert F_r >= FT(0) # rime mass fraction must be positive ... - @assert F_r < FT(1) # ... and there must always be some unrimed part + @assert F_rim >= FT(0) # rime mass fraction must be positive ... + @assert F_rim < FT(1) # ... and there must always be some unrimed part - if F_r == FT(0) + if F_rim == FT(0) return (; D_cr = FT(0), D_gr = FT(0), ρ_g = FT(0), ρ_d = FT(0)) else @assert ρ_r > FT(0) # rime density must be positive ... @@ -110,8 +111,8 @@ function thresholds(p3::PSP3{FT}, ρ_r::FT, F_r::FT) where {FT} P3_problem(ρ_d) = ρ_d - ρ_d_helper( p3, - D_cr_helper(p3, F_r, ρ_g_helper(ρ_r, F_r, ρ_d)), - D_gr_helper(p3, ρ_g_helper(ρ_r, F_r, ρ_d)), + D_cr_helper(p3, F_rim, ρ_g_helper(ρ_r, F_rim, ρ_d)), + D_gr_helper(p3, ρ_g_helper(ρ_r, F_rim, ρ_d)), ) ρ_d = @@ -120,10 +121,10 @@ function thresholds(p3::PSP3{FT}, ρ_r::FT, F_r::FT) where {FT} RS.SecantMethod(FT(0), FT(1000)), RS.CompactSolution(), ).root - ρ_g = ρ_g_helper(ρ_r, F_r, ρ_d) + ρ_g = ρ_g_helper(ρ_r, F_rim, ρ_d) return (; - D_cr = D_cr_helper(p3, F_r, ρ_g), + D_cr = D_cr_helper(p3, F_rim, ρ_g), D_gr = D_gr_helper(p3, ρ_g), ρ_g, ρ_d, @@ -132,26 +133,39 @@ function thresholds(p3::PSP3{FT}, ρ_r::FT, F_r::FT) where {FT} end """ -p3_density(p3, D, F_r, th) + p3_F_liq_average(F_liq, X_ice, X_liq) + + - F_liq - liquid fraction (L_liq / L_p3_tot) + - X_ice - ice core parameterization (i.e. mass, etc) + - X_liq - liquid part parameterization + +Returns the liquid fraction weighted average of X_ice and X_liq. +""" +function p3_F_liq_average(F_liq::FT, X_ice::FT, X_liq::FT) where {FT} + return (1 - F_liq) * X_ice + F_liq * X_liq +end + +""" + p3_density(p3, D, F_rim, th) - p3 - a struct with P3 parameters - D - maximum particle dimension [m] -- F_r - rime mass fraction (L_rim / L_ice) [-] +- F_rim - rime mass fraction (L_rim / L_ice) [-] - th - P3 scheme thresholds() output tuple (D_cr, D_gr, ρ_g, ρ_d) Returns the density of a particle based on where it falls in the particle-size-based properties regime. Following Morrison and Milbrandt (2015), the density of nonspherical particles is assumed to be the particle mass divided by the volume of a sphere with the same D. +Needed for aspect ratio calculation, so we assume zero liquid fraction. """ -function p3_density(p3::PSP3, D::FT, F_r::FT, th) where {FT} - D_th = D_th_helper(p3) - if D_th > D +function p3_density(p3::PSP3, D::FT, F_rim::FT, th) where {FT} + if D_th_helper(p3) > D # small spherical ice return p3.ρ_i - elseif F_r == 0 + elseif F_rim == 0 # large nonspherical unrimed ice return (6 * α_va_si(p3)) / FT(π) * D^(p3.β_va - 3) - elseif th.D_gr > D >= D_th + elseif th.D_gr > D >= D_th_helper(p3) # dense nonspherical ice return (6 * α_va_si(p3)) / FT(π) * D^(p3.β_va - 3) elseif th.D_cr > D >= th.D_gr @@ -159,39 +173,70 @@ function p3_density(p3::PSP3, D::FT, F_r::FT, th) where {FT} return th.ρ_g else #elseif D >= th.D_cr # partially rimed ice - return (6 * α_va_si(p3)) / (FT(π) * (1 - F_r)) * D^(p3.β_va - 3) + return (6 * α_va_si(p3)) / (FT(π) * (1 - F_rim)) * D^(p3.β_va - 3) end end """ - mass_(p3, D, ρ, F_r) + mass_(p3, D, ρ, F_rim) - p3 - a struct with P3 scheme parameters - D - maximum particle dimension [m] - ρ - bulk ice density (ρ_i for small ice, ρ_g for graupel) [kg/m3] - - F_r - rime mass fraction [L_rim/L_ice] + - F_rim - rime mass fraction [L_rim/L_ice] Returns mass as a function of size for differen particle regimes [kg] """ -# for spherical ice (small ice or completely rimed ice) +# for spherical particles (small ice, completely rimed ice, or liquid on ice) mass_s(D::FT, ρ::FT) where {FT <: Real} = FT(π) / 6 * ρ * D^3 # for large nonspherical ice (used for unrimed and dense types) mass_nl(p3::PSP3, D::FT) where {FT <: Real} = α_va_si(p3) * D^p3.β_va # for partially rimed ice -mass_r(p3::PSP3, D::FT, F_r::FT) where {FT <: Real} = - α_va_si(p3) / (1 - F_r) * D^p3.β_va +mass_r(p3::PSP3, D::FT, F_rim::FT) where {FT <: Real} = + α_va_si(p3) / (1 - F_rim) * D^p3.β_va """ - p3_mass(p3, D, F_r, th) + p3_mass(p3, D, F_rim, F_liq, th) - p3 - a struct with P3 scheme parameters - D - maximum particle dimension - - F_r - rime mass fraction (L_rim/L_ice) + - F_rim - rime mass fraction (L_rim / L_ice) + - F_liq - liquid fraction (L_liq / L_p3_tot) - th - P3 scheme thresholds() output tuple (D_cr, D_gr, ρ_g, ρ_d) Returns mass(D) regime, used to create figures for the docs page. """ function p3_mass( + p3::PSP3, + D::FT, + F_rim::FT, + F_liq::FT, + th = (; D_cr = FT(0), D_gr = FT(0), ρ_g = FT(0), ρ_d = FT(0)), +) where {FT} + if D_th_helper(p3) > D + return p3_F_liq_average(F_liq, mass_s(D, p3.ρ_i), mass_s(D, p3.ρ_l)) # small spherical ice + elseif F_rim == 0 + return p3_F_liq_average(F_liq, mass_nl(p3, D), mass_s(D, p3.ρ_l)) # large nonspherical unrimed ice + elseif th.D_gr > D >= D_th_helper(p3) + return p3_F_liq_average(F_liq, mass_nl(p3, D), mass_s(D, p3.ρ_l)) # dense nonspherical ice + elseif th.D_cr > D >= th.D_gr + return p3_F_liq_average(F_liq, mass_s(D, th.ρ_g), mass_s(D, p3.ρ_l)) # graupel + else #elseif D >= th.D_cr + return p3_F_liq_average(F_liq, mass_r(p3, D, F_rim), mass_s(D, p3.ρ_l)) # partially rimed ice + end +end + +""" + p3_dmdD(p3, D, F_r, th) + + - p3 - a struct containing p3 parameters + - D - maximum dimension of the particle + - F_r - rime mass fraction (q_rim/ q_i) + - th - thresholds as calculated by thresholds() + +Returns dm(D)/dD for each particle regime, used in snow melting +""" +function p3_dmdD( p3::PSP3, D::FT, F_r::FT, @@ -199,15 +244,15 @@ function p3_mass( ) where {FT <: Real} D_th = D_th_helper(p3) if D_th > D - return mass_s(D, p3.ρ_i) # small spherical ice + return FT(π) / 2 * p3.ρ_i * D^2 elseif F_r == 0 - return mass_nl(p3, D) # large nonspherical unrimed ice + return α_va_si(p3) * p3.β_va * D^(p3.β_va - 1) elseif th.D_gr > D >= D_th - return mass_nl(p3, D) # dense nonspherical ice + return α_va_si(p3) * p3.β_va * D^(p3.β_va - 1) elseif th.D_cr > D >= th.D_gr - return mass_s(D, th.ρ_g) # graupel - else #elseif D >= th.D_cr - return mass_r(p3, D, F_r) # partially rimed ice + return FT(π) / 2 * th.ρ_g * D^2 + else # D >= th.D_cr + return α_va_si(p3) / (1 - F_r) * p3.β_va * D^(p3.β_va - 1) end end @@ -220,19 +265,20 @@ end Returns particle projected area as a function of size for different particle regimes """ # for spherical particles -A_s(D::FT) where {FT <: Real} = FT(π) / 4 * D^2 +A_s(D::FT) where {FT} = FT(π) / 4 * D^2 # for nonspherical particles -A_ns(p3::PSP3, D::FT) where {FT <: Real} = p3.γ * D^p3.σ +A_ns(p3::PSP3, D::FT) where {FT} = p3.γ * D^p3.σ # partially rimed ice -A_r(p3::PSP3, F_r::FT, D::FT) where {FT <: Real} = - F_r * A_s(D) + (1 - F_r) * A_ns(p3, D) +A_r(p3::PSP3, F_rim::FT, D::FT) where {FT} = + F_rim * A_s(D) + (1 - F_rim) * A_ns(p3, D) """ - p3_area(p3, D, F_r, th) + p3_area(p3, D, F_rim, F_liq, th) - p3 - a struct with P3 scheme parameters - D - maximum particle dimension - - F_r - rime mass fraction (L_rim/L_ice) + - F_rim - rime mass fraction (L_rim / L_ice) + - F_liq - liquid fraction (L_liq / L_p3_tot) - th - P3 scheme thresholds() output tuple (D_cr, D_gr, ρ_g, ρ_d) Returns area(D), used to create figures for the documentation. @@ -240,21 +286,44 @@ Returns area(D), used to create figures for the documentation. function p3_area( p3::PSP3, D::FT, - F_r::FT, + F_rim::FT, + F_liq::FT, th = (; D_cr = FT(0), D_gr = FT(0), ρ_g = FT(0), ρ_d = FT(0)), -) where {FT <: Real} +) where {FT} # Area regime: if D_th_helper(p3) > D - return A_s(D) # small spherical ice - elseif F_r == 0 - return A_ns(p3, D) # large nonspherical unrimed ice + return A_s(D) # small spherical ice + elseif F_rim == 0 + return p3_F_liq_average(F_liq, A_ns(p3, D), A_s(D)) # large nonspherical unrimed ice elseif th.D_gr > D >= D_th_helper(p3) - return A_ns(p3, D) # dense nonspherical ice + return p3_F_liq_average(F_liq, A_ns(p3, D), A_s(D)) # dense nonspherical ice elseif th.D_cr > D >= th.D_gr - return A_s(D) # graupel - elseif D >= th.D_cr - return A_r(p3, F_r, D) # partially rimed ice - else - throw("D not in range") + return A_s(D) # graupel + else # D >= th.D_cr + return p3_F_liq_average(F_liq, A_r(p3, F_rim, D), A_s(D)) # partially rimed ice end end + +""" + ϕᵢ(p3, D, F_rim, th) + + - p3 - a struct containing P3 parameters + - D - maximum dimension of ice particle [m] + - F_rim - rime mass fraction (L_rim/ L_ice) [-] + - th - P3 particle properties thresholds + +Returns the aspect ratio (ϕ) for an ice particle with mass, cross-sectional area, +and ice density determined using the size-dependent particle property regimes +following Morrison and Milbrandt (2015). The density of nonspherical +particles is assumed to be equal to the particle mass divided by the volume of a +spherical particle with the same D_max. +Assuming zero liquid fraction. +""" +function ϕᵢ(p3::PSP3, D::FT, F_rim::FT, th) where {FT} + F_liq = FT(0) + mᵢ = p3_mass(p3, D, F_rim, F_liq, th) + aᵢ = p3_area(p3, D, F_rim, F_liq, th) + ρᵢ = p3_density(p3, D, F_rim, th) + + return ifelse(D == 0, FT(0), 16 * ρᵢ^2 * aᵢ^3 / (9 * FT(π) * mᵢ^2)) +end diff --git a/src/P3_processes.jl b/src/P3_processes.jl index b2798c33d..07e528770 100644 --- a/src/P3_processes.jl +++ b/src/P3_processes.jl @@ -45,3 +45,72 @@ function het_ice_nucleation( return (; dNdt, dLdt) end + +""" + ice_melt(p3, Chen2022, aps, tps, L_ice, N_ice, Tₐ, ρₐ, F_rim, ρ_rim, dt) + + - p3 - a struct containing p3 parameters + - Chen2022 - struct containing Chen 2022 velocity parameters + - aps - air properties + - tps - thermodynamics parameters + - L_ice - ice content + - N_ice - ice number concentration + - T - temperature (K) + - ρ_a - air density + - F_r - rime mass fraction (q_rim/ q_i) + - ρ_r - rime density (q_rim/B_rim) + - dt - model time step (for limiting the tendnecy) + +Returns the melting rate of ice (QIMLT in Morrison and Mildbrandt (2015)). +""" +function ice_melt( + p3::PSP3, + Chen2022::CMP.Chen2022VelTypeSnowIce, + aps::CMP.AirProperties{FT}, + tps::TDP.ThermodynamicsParameters{FT}, + L_ice::FT, + N_ice::FT, + Tₐ::FT, + ρₐ::FT, + F_rim::FT, + ρ_rim::FT, + dt::FT, +) where {FT} + # process not dependent on F_liq + # (we want ice core shape params) + F_liq_ = FT(0) + # Get constants + (; ν_air, D_vapor, K_therm) = aps + L_f = TD.latent_heat_fusion(tps, Tₐ) + N_sc = ν_air / D_vapor + + # Get the P3 diameter distribution... + th = thresholds(p3, ρ_rim, F_rim) + (λ, N_0) = + distribution_parameter_solver(p3, L_ice, N_ice, ρ_rim, F_rim, F_liq_) + N(D) = N′ice(p3, D, λ, N_0) + # ... and D_max for the integral + bound = get_ice_bound(p3, λ, N_ice, FT(1e-6)) + + # Ice particle terminal velocity + v(D) = ice_particle_terminal_velocity(p3, D, Chen2022, ρₐ, F_rim, th) + # Reynolds number + N_Re(D) = D * v(D) / ν_air + # Ventillation factor + F_v(D) = p3.vent_a + p3.vent_b * N_sc^(1 / 3) * N_Re(D)^(1 / 2) + dmdD(D) = p3_dmdD(p3, D, F_rim, th) + + f(D) = 4 * K_therm / L_f * (Tₐ - p3.T_freeze) * dmdD(D) / D * F_v(D) * N(D) + # Integrate + (dLdt, error) = QGK.quadgk(d -> f(d), 0, bound, rtol = FT(1e-6)) + + # only consider melting (not fusion) + dLdt = max(0, dLdt) + # compute change of N_ice proportional to change in L + dNdt = N_ice / L_ice * dLdt + + # ... and dont exceed the available number and mass of water droplets + dNdt = min(dNdt, N_ice / dt) + dLdt = min(dLdt, L_ice / dt) + return (; dNdt, dLdt) +end diff --git a/src/P3_size_distribution.jl b/src/P3_size_distribution.jl index a20417ab0..848a47c5a 100644 --- a/src/P3_size_distribution.jl +++ b/src/P3_size_distribution.jl @@ -50,51 +50,64 @@ end N'(D) = N0 * D ^ μ * exp(-λD)) at given D """ function N′ice(p3::PSP3, D::FT, λ::FT, N_0::FT) where {FT} - return N_0 * D^DSD_μ(p3, λ) * exp(-λ * D) + #@assert D > FT(0) + return ifelse(D > 0, N_0 * D^DSD_μ(p3, λ) * exp(-λ * D), FT(0)) end """ - get_ice_bound(p3, λ, tolerance) + get_ice_bound(p3, λ, N, rtol) - p3 - a struct containing p3 parameters - λ - shape parameters of ice distribution - - tolerance - tolerance for how much of the distribution we want to integrate over - - Returns the bound on the distribution that would guarantee that 1-tolerance - of the ice distribution is integrated over. This is calculated by setting - N_0(1 - tolerance) = ∫ N'(D) dD from 0 to bound and solving for bound. - This was further simplified to cancel out the N_0 from both sides. - The guess was calculated through a linear approximation extrapolated from - numerical solutions. + - N - ice number concentration + - rtol - relative tolerance for root solver + + Ice size distribution upper bound for numerical integrals, such that + total number concentration is equal to the integral over the size distribution + with relative tolerance rtol. """ -function get_ice_bound(p3::PSP3, λ::FT, tolerance::FT) where {FT} - ice_problem(x) = - tolerance - Γ(1 + DSD_μ(p3, λ), FT(exp(x)) * λ) / Γ(1 + DSD_μ(p3, λ)) - guess = log(19 / 6 * (DSD_μ(p3, λ) - 1) + 39) - log(λ) - log_ice_x = - RS.find_zero( - ice_problem, - RS.SecantMethod(guess - 1, guess), - RS.CompactSolution(), - RS.RelativeSolutionTolerance(eps(FT)), - 5, - ).root +function get_ice_bound(p3::PSP3, λ::FT, N::FT, rtol::FT) where {FT} + #return FT(1e-2) - From what I'm seeing so far, this is not such a bad guess + #We might want to reconsider using it in the future (after more testing). + μ = DSD_μ(p3, λ) + N₀ = DSD_N₀(p3, N, λ) + + ice_bound_problem(D_max_hat) = + N - N₀ / λ^(μ + 1) * Γ_lower(μ + 1, D_max_hat) + + D_guess_low = FT(1e-6) + D_guess_high = FT(1e-1) - return exp(log_ice_x) + if ice_bound_problem(D_guess_low * λ) <= rtol + return D_guess_low + elseif ice_bound_problem(D_guess_high * λ) <= rtol + return D_guess_high + else + bound = + RS.find_zero( + ice_bound_problem, + RS.SecantMethod(D_guess_low * λ, D_guess_high * λ), + RS.CompactSolution(), + RS.RelativeSolutionTolerance(rtol), + 5, + ).root + return bound / λ + end end """ - L_(p3, ρ, F_r, λ, μ, D_min, D_max) + L_(p3, ρ, F_rim, F_liq, λ, μ, D_min, D_max) - p3 - a struct with P3 scheme parameters - ρ - bulk ice density (ρ_i for small ice, ρ_g for graupel) [kg/m^3] - - F_r - rime mass fraction [L_rim/L_ice] + - F_rim - rime mass fraction (L_rim / L_ice) [-] - μ - shape parameter of N′ gamma distribution - λ - slope parameter of N′ gamma distribution - D_min - minimum bound for regime - D_max - maximum bound for regime (if not specified, then infinity) - Returns ice content for a given m(D) regime +Returns ice content for a given m(D) regime. +Liquid fraction weighted averaging is not computed here but in L_over_N_gamma(). """ # small, spherical ice or graupel (completely rimed, spherical) # D_min = 0, D_max = D_th, ρ = ρᵢ @@ -113,15 +126,21 @@ function L_n(p3::PSP3, μ::FT, λ::FT, D_min::FT, D_max::FT) where {FT} end # partially rimed ice or large unrimed ice (upper bound on D is infinity) # L_rim > 0 and D_min = D_cr, D_max = inf -function L_r(p3::PSP3, F_r::FT, μ::FT, λ::FT, D_min::FT) where {FT} - return ∫_Γ(D_min, FT(Inf), α_va_si(p3) / (1 - F_r), μ + p3.β_va, λ) +function L_r(p3::PSP3, F_rim::FT, μ::FT, λ::FT, D_min::FT) where {FT} + return ∫_Γ(D_min, FT(Inf), α_va_si(p3) / (1 - F_rim), μ + p3.β_va, λ) +end +# F_liq != 0 (liquid mass on mixed-phase particles for D in [D_min, D_max]) +function L_liq(p3::PSP3, μ::FT, λ::FT) where {FT} + return ∫_Γ(FT(0), FT(Inf), (FT(π) / 6) * p3.ρ_l, μ + 3, λ) end """ - L_over_N_gamma(p3, F_r, λ, th) + L_over_N_gamma(p3, F_rim, F_liq, log_λ, th) - p3 - a struct with P3 scheme parameters - - F_r - rime mass fraction [L_rim/L_ice] + - F_rim - rime mass fraction (L_rim / L_ice) [-] + - F_liq - liquid fraction (L_liq / L_p3_tot) [-]: + - zero if solving for ice core shape parameters - log_λ - logarithm of the slope parameter of N′ gamma distribution - μ - shape parameter of N′ gamma distribution - th - P3 scheme thresholds() output tuple (D_cr, D_gr, ρ_g, ρ_d) @@ -131,7 +150,8 @@ Eq. 5 in Morrison and Milbrandt (2015). """ function L_over_N_gamma( p3::PSP3, - F_r::FT, + F_rim::FT, + F_liq::FT, log_λ::FT, μ::FT, th = (; D_cr = FT(0), D_gr = FT(0), ρ_g = FT(0), ρ_d = FT(0)), @@ -140,41 +160,54 @@ function L_over_N_gamma( D_th = D_th_helper(p3) λ = exp(log_λ) N = Γ(1 + μ) / (λ^(1 + μ)) - return ifelse( - F_r == FT(0), - (L_s(p3, p3.ρ_i, μ, λ, FT(0), D_th) + L_rz(p3, μ, λ, D_th)) / N, - ( + F_rim == FT(0), + (p3_F_liq_average( + F_liq, + L_s(p3, p3.ρ_i, μ, λ, FT(0), D_th) + L_rz(p3, μ, λ, D_th), + L_liq(p3, μ, λ), + )) / N, + (p3_F_liq_average( + F_liq, L_s(p3, p3.ρ_i, μ, λ, FT(0), D_th) + L_n(p3, μ, λ, D_th, th.D_gr) + L_s(p3, th.ρ_g, μ, λ, th.D_gr, th.D_cr) + - L_r(p3, F_r, μ, λ, th.D_cr) - ) / N, + L_r(p3, F_rim, μ, λ, th.D_cr), + L_liq(p3, μ, λ), + )) / N, ) end """ - DSD_μ_approx(p3, L_ice, N, ρ_r, F_r) + DSD_μ_approx(p3, L, N, ρ_r, F_rim, F_liq) - p3 - a struct with P3 scheme parameters - - L_ice - ice content [kg/m3] + - L - ice content [kg/m3] - N - total ice number concentration [1/m3] - ρ_r - rime density (L_rim/B_rim) [kg/m^3] - - F_r - rime mass fraction (L_rim/L_ice) - + - F_rim - rime mass fraction (L_rim / L_ice) [-] + - F_liq - liquid fraction (L_liq / L_p3_tot) [-]: + - zero if solving for ice core shape parameters Returns the approximated shape parameter μ for a given q and N value """ -function DSD_μ_approx(p3::PSP3, L::FT, N::FT, ρ_r::FT, F_r::FT) where {FT} - # Get thresholds for given F_r, ρ_r - th = thresholds(p3, ρ_r, F_r) +function DSD_μ_approx( + p3::PSP3, + L::FT, + N::FT, + ρ_r::FT, + F_rim::FT, + F_liq::FT, +) where {FT} + # Get thresholds for given F_rim, ρ_r + th = thresholds(p3, ρ_r, F_rim) # Get min and max lambda values λ_0 = μ_to_λ(p3, FT(0)) λ_6 = μ_to_λ(p3, p3.μ_max) - # Get corresponding L/N values at given F_r - L_over_N_min = log(L_over_N_gamma(p3, F_r, log(λ_0), FT(0), th)) - L_over_N_max = log(L_over_N_gamma(p3, F_r, log(λ_6), p3.μ_max, th)) + # Get corresponding L/N values at given F_rim + L_over_N_min = log(L_over_N_gamma(p3, F_rim, F_liq, log(λ_0), FT(0), th)) + L_over_N_max = log(L_over_N_gamma(p3, F_rim, F_liq, log(λ_6), p3.μ_max, th)) # Return approximation between them μ = (p3.μ_max / (L_over_N_max - L_over_N_min)) * (log(L / N) - L_over_N_min) @@ -184,13 +217,14 @@ function DSD_μ_approx(p3::PSP3, L::FT, N::FT, ρ_r::FT, F_r::FT) where {FT} end """ - get_bounds(N, L, F_r, p3, th) + get_bounds(N, L, F_rim, F_liq, p3, th) - N - ice number concentration [1/m3] - L - ice content [kg/m3] - μ - shape parameter of N′ gamma distribution - - F_r -rime mass fraction [L_rim/L_ice] - - p3 - a struct with P3 scheme parameters + - F_rim - rime mass fraction (L_rim / L_ice) [-] + - F_liq - liquid fraction (L_liq / L_p3_tot) [-]: + - zero if solving for ice core shape parameters - p3 - a struct with P3 scheme parameters - th - P3 scheme thresholds() output tuple (D_cr, D_gr, ρ_g, ρ_d) Returns estimated guess for λ from L to be used in distribution_parameter_solver() @@ -199,7 +233,8 @@ function get_bounds( N::FT, L::FT, μ::FT, - F_r::FT, + F_rim::FT, + F_liq::FT, p3::PSP3, th = (; D_cr = FT(0), D_gr = FT(0), ρ_g = FT(0), ρ_d = FT(0)), ) where {FT} @@ -219,8 +254,8 @@ function get_bounds( radius = FT(0.2) end - Ll = L_over_N_gamma(p3, F_r, log(left), μ, th) - Lr = L_over_N_gamma(p3, F_r, log(right), μ, th) + Ll = L_over_N_gamma(p3, F_rim, F_liq, log(left), μ, th) + Lr = L_over_N_gamma(p3, F_rim, F_liq, log(right), μ, th) guess = left * (goal / (Ll))^((log(right) - log(left)) / (log(Lr) - log(Ll))) @@ -235,10 +270,12 @@ end distrbution_parameter_solver() - p3 - a struct with P3 scheme parameters - - L - mass mixing ratio - - N - number mixing ratio + - L - ice content [kg/m3] + - N - number concentration [1/m3] - ρ_r - rime density (L_rim/B_rim) [kg/m^3] - - F_r - rime mass fraction (L_rim/L_ice) + - F_rim - rime mass fraction (L_rim / L_ice) [-] + - F_liq - liquid fraction (L_liq / L_p3_tot) [-]: + - zero if solving for ice core shape parameters - p3 - a struct with P3 scheme parameters Solves the nonlinear system consisting of N_0 and λ for P3 prognostic variables Returns a named tuple containing: @@ -250,19 +287,20 @@ function distribution_parameter_solver( L::FT, N::FT, ρ_r::FT, - F_r::FT, + F_rim::FT, + F_liq::FT, ) where {FT} # Get the thresholds for different particles regimes - th = thresholds(p3, ρ_r, F_r) + th = thresholds(p3, ρ_r, F_rim) # Get μ given L and N - μ = DSD_μ_approx(p3, L, N, ρ_r, F_r) + μ = DSD_μ_approx(p3, L, N, ρ_r, F_rim, F_liq) # To ensure that λ is positive solve for x such that λ = exp(x) - shape_problem(x) = L / N - L_over_N_gamma(p3, F_r, x, μ, th) + shape_problem(x) = L / N - L_over_N_gamma(p3, F_rim, F_liq, x, μ, th) # Get intial guess for solver - (; min, max) = get_bounds(N, L, μ, F_r, p3, th) + (; min, max) = get_bounds(N, L, μ, F_rim, F_liq, p3, th) # Find slope parameter x = @@ -278,39 +316,46 @@ function distribution_parameter_solver( end """ - D_m (p3, L, N, ρ_r, F_r) + D_m (p3, L, N, ρ_r, F_rim) - p3 - a struct with P3 scheme parameters - - L - mass mixing ratio - - N - number mixing ratio + - L - ice mass content [kg/m3] + - N - number concentration [1/m3] - ρ_r - rime density (L_rim/B_rim) [kg/m^3] - - F_r - rime mass fraction (L_rim/L_ice) + - F_rim - rime mass fraction (L_rim / L_ice) [-] + - F_liq - liquid fraction (L_liq / L_p3_tot) [-]: + - zero if solving for ice core D_m + - p3 - a struct with P3 scheme parameters Return the mass weighted mean particle size [m] """ -function D_m(p3::PSP3, L::FT, N::FT, ρ_r::FT, F_r::FT) where {FT} +function D_m(p3::PSP3, L::FT, N::FT, ρ_r::FT, F_rim::FT, F_liq::FT) where {FT} # Get the thresholds for different particles regimes - th = thresholds(p3, ρ_r, F_r) + th = thresholds(p3, ρ_r, F_rim) D_th = D_th_helper(p3) # Get the shape parameters - (λ, N_0) = distribution_parameter_solver(p3, L, N, ρ_r, F_r) + (λ, N_0) = distribution_parameter_solver(p3, L, N, ρ_r, F_rim, F_liq) μ = DSD_μ(p3, λ) # Redefine α_va to be in si units α_va = α_va_si(p3) # Calculate numerator - n = 0 - if F_r == 0 - n += ∫_Γ(FT(0), D_th, π / 6 * p3.ρ_i * N_0, μ + 4, λ) - n += ∫_Γ(D_th, Inf, α_va * N_0, μ + p3.β_va + 1, λ) + n_nl = 0 + if F_rim == 0 + n_nl += ∫_Γ(FT(0), D_th, π / 6 * p3.ρ_i * N_0, μ + 4, λ) + n_nl += ∫_Γ(D_th, FT(Inf), α_va * N_0, μ + p3.β_va + 1, λ) + n_l = ∫_Γ(FT(0), FT(Inf), N_0 * p3.ρ_l * (FT(π) / 6), μ + 4, λ) else - n += ∫_Γ(FT(0), D_th, π / 6 * p3.ρ_i * N_0, μ + 4, λ) - n += ∫_Γ(D_th, th.D_gr, α_va * N_0, μ + p3.β_va + 1, λ) - n += ∫_Γ(th.D_gr, th.D_cr, π / 6 * th.ρ_g * N_0, μ + 4, λ) - n += ∫_Γ(th.D_cr, Inf, α_va / (1 - F_r) * N_0, μ + p3.β_va + 1, λ) + n_nl += ∫_Γ(FT(0), D_th, π / 6 * p3.ρ_i * N_0, μ + 4, λ) + n_nl += ∫_Γ(D_th, th.D_gr, α_va * N_0, μ + p3.β_va + 1, λ) + n_nl += ∫_Γ(th.D_gr, th.D_cr, π / 6 * th.ρ_g * N_0, μ + 4, λ) + n_nl += + ∫_Γ(th.D_cr, FT(Inf), α_va / (1 - F_rim) * N_0, μ + p3.β_va + 1, λ) + n_l = ∫_Γ(FT(0), FT(Inf), N_0 * p3.ρ_l * (FT(π) / 6), μ + 4, λ) end - # Normalize by L - return n / L + + # compute F_liq-weighted average and normalize by L + return ifelse(L < eps(FT), FT(0), p3_F_liq_average(F_liq, n_nl, n_l) / L) end diff --git a/src/P3_terminal_velocity.jl b/src/P3_terminal_velocity.jl index 0621e1121..c20569a22 100644 --- a/src/P3_terminal_velocity.jl +++ b/src/P3_terminal_velocity.jl @@ -1,33 +1,11 @@ """ - ϕᵢ(p3, D, F_r, th) - - - p3 - a struct containing P3 parameters - - D - maximum dimension of ice particle [m] - - F_r - rime mass fraction (L_rim/ L_ice) [-] - - th - P3 particle properties thresholds - -Returns the aspect ratio (ϕ) for an ice particle with mass, cross-section area, -and ice density determined using the size-dependent particle property regimes -following Morrison and Milbrandt (2015). The density of nonspherical -particles is assumed to be equal to the particle mass divided by the volume of a -spherical particle with the same D_max. -""" -function ϕᵢ(p3::PSP3, D::FT, F_r::FT, th) where {FT} - - mᵢ = p3_mass(p3, D, F_r, th) - aᵢ = p3_area(p3, D, F_r, th) - ρᵢ = p3_density(p3, D, F_r, th) - - return ifelse(D == 0, FT(0), 16 * ρᵢ^2 * aᵢ^3 / (9 * FT(π) * mᵢ^2)) -end - -""" - ice_particle_terminal_velocity(p3, D, Chen2022, ρₐ, F_r, th, use_aspect_ratio) + ice_particle_terminal_velocity(p3, D, Chen2022, ρₐ, F_rim, th, use_aspect_ratio) + - p3 - p3 parameters - D - maximum particle dimension - Chen2022 - a struct with terminal velocity parameters from Chen 2022 - ρₐ - air density - - F_r - rime mass fraction (L_rim/ L_ice) + - F_rim - rime mass fraction (L_rim/L_ice) [-] - th - P3 particle properties thresholds - use_aspect_ratio - Bool flag set to true if we want to consider the effects of particle aspect ratio on its terminal velocity (default: true) @@ -40,7 +18,7 @@ function ice_particle_terminal_velocity( D::FT, Chen2022::CMP.Chen2022VelTypeSnowIce, ρₐ::FT, - F_r::FT, + F_rim::FT, th, use_aspect_ratio = true, ) where {FT} @@ -51,11 +29,51 @@ function ice_particle_terminal_velocity( end v = sum(@. sum(ai * D^bi * exp(-ci * D))) - return ifelse(use_aspect_ratio, ϕᵢ(p3, D, F_r, th)^FT(1 / 3) * v, v) + return ifelse(use_aspect_ratio, ϕᵢ(p3, D, F_rim, th)^FT(1 / 3) * v, v) end """ - velocity_difference(type, Dₗ, Dᵢ, p3, Chen2022, ρ_a, F_r, th, use_aspect_ratio) + p3_particle_terminal_velocity(p3, D, Chen2022, ρₐ, F_rim, F_liq, th, use_aspect_ratio) + + - p3 - p3 parameters + - D - maximum particle dimension + - Chen2022 - struct with terminal velocity parameters from Chen 2022 + - ρₐ - air density + - F_rim - rime mass fraction (L_rim/L_ice) [-] + - F_liq - liquid fraction (L_liq / L_p3_tot) [-] + - th - thresholds as calculated by thresholds() + - use_aspect_ratio - Bool flag set to true if we want to consider the effects + of particle aspect ratio on its terminal velocity (default: true) + +Returns the terminal velocity of a single mixed-phase particle using the Chen 2022 +parametrizations by computing an F_liq-weighted average of solid and liquid +phase terminal velocities, using the maximum dimension of the whole particle for both. +""" +function p3_particle_terminal_velocity( + p3::PSP3, + D::FT, + Chen2022::CMP.Chen2022VelType, + ρₐ::FT, + F_rim::FT, + F_liq::FT, + th, + use_aspect_ratio = true, +) where {FT} + v_i = ice_particle_terminal_velocity( + p3, + D, + Chen2022.snow_ice, + ρₐ, + F_rim, + th, + use_aspect_ratio, + ) + v_r = CM2.rain_particle_terminal_velocity(D, Chen2022.rain, ρₐ) + return p3_F_liq_average(F_liq, v_i, v_r) +end + +""" + velocity_difference(type, Dₗ, Dᵢ, p3, Chen2022, ρₐ, F_rim, F_liq, th, aspect_ratio) - type - a struct containing the size distribution parameters of the particle colliding with ice - Dₗ - maximum dimension of the particle colliding with ice @@ -63,12 +81,13 @@ end - p3 - a struct containing P3 parameters - Chen2022 - a struct containing Chen 2022 velocity parameters - ρ_a - density of air - - F_r - rime mass fraction (L_rim/ L_ice) + - F_rim - rime mass fraction (L_rim/L_ice) [-] + - F_liq - liquid fraction (L_liq / L_p3_tot) [-] - th - P3 particle properties thresholds - use_aspect_ratio - Bool flag set to true if we want to consider the effects of particle aspect ratio on its terminal velocity (default: true) -Returns the absolute value of the velocity difference between an ice particle and +Returns the absolute value of the velocity difference between a mixed-phase particle and cloud or rain drop as a function of their sizes. It uses Chen 2022 velocity parameterization for ice and rain and assumes no sedimentation of cloud droplets. """ @@ -82,18 +101,19 @@ function velocity_difference( p3::PSP3, Chen2022::CMP.Chen2022VelType, ρₐ::FT, - F_r::FT, + F_rim::FT, th, use_aspect_ratio = true, ) where {FT} # velocity difference for rain-ice collisions return abs( - ice_particle_terminal_velocity( + p3_particle_terminal_velocity( p3, Dᵢ, - Chen2022.snow_ice, + Chen2022, ρₐ, - F_r, + F_rim, + F_liq, th, use_aspect_ratio, ) - CM2.rain_particle_terminal_velocity(Dₗ, Chen2022.rain, ρₐ), @@ -106,18 +126,20 @@ function velocity_difference( p3::PSP3, Chen2022::CMP.Chen2022VelType, ρₐ::FT, - F_r::FT, + F_rim::FT, + F_liq::FT, th, use_aspect_ratio = true, ) where {FT} # velocity difference for cloud-ice collisions return abs( - ice_particle_terminal_velocity( + p3_particle_terminal_velocity( p3, Dᵢ, - Chen2022.snow_ice, + Chen2022, ρₐ, - F_r, + F_rim, + F_liq, th, use_aspect_ratio, ), @@ -125,15 +147,15 @@ function velocity_difference( end """ - ice_terminal_velocity(p3, Chen2022, L, N, ρ_r, F_r, ρ_a, use_aspect_ratio) + ice_terminal_velocity(p3, Chen2022, L, N, ρ_r, F_rim, ρₐ, use_aspect_ratio) - p3 - a struct with P3 scheme parameters - Chen2022 - a struch with terminal velocity parameters as in Chen(2022) - - L - mass mixing ratio - - N - number mixing ratio + - L - ice mass content [kg/m3] + - N - number concentration [1/m3] - ρ_r - rime density (L_rim/B_rim) [kg/m^3] - - F_r - rime mass fraction (L_rim/q_ice) - - ρ_a - density of air + - F_rim - rime mass fraction (L_rim/L_ice) + - ρₐ - density of air - use_aspect_ratio - Bool flag set to true if we want to consider the effects of particle aspect ratio on its terminal velocity (default: true) @@ -142,11 +164,12 @@ eq C10 of Morrison and Milbrandt (2015) and using Chen 2022 terminal velocity sc """ function ice_terminal_velocity( p3::PSP3, - Chen2022::CMP.Chen2022VelTypeSnowIce, + Chen2022::CMP.Chen2022VelType, L::FT, N::FT, ρ_r::FT, - F_r::FT, + F_rim::FT, + F_liq::FT, ρₐ::FT, use_aspect_ratio = true, ) where {FT} @@ -154,23 +177,24 @@ function ice_terminal_velocity( return FT(0), FT(0) else # get the particle properties thresholds - th = thresholds(p3, ρ_r, F_r) + th = thresholds(p3, ρ_r, F_rim) # get the size distribution parameters - (λ, N₀) = distribution_parameter_solver(p3, L, N, ρ_r, F_r) + (λ, N₀) = distribution_parameter_solver(p3, L, N, ρ_r, F_rim, F_liq) # get the integral limit - D_max = get_ice_bound(p3, λ, eps(FT)) + D_max = get_ice_bound(p3, λ, N, 1e-8) # ∫N(D) m(D) v(D) dD v_m = QGK.quadgk( D -> N′ice(p3, D, λ, N₀) * - p3_mass(p3, D, F_r, th) * - ice_particle_terminal_velocity( + p3_mass(p3, D, F_rim, F_liq, th) * + p3_particle_terminal_velocity( p3, D, Chen2022, ρₐ, - F_r, + F_rim, + F_liq, th, use_aspect_ratio, ), @@ -182,12 +206,13 @@ function ice_terminal_velocity( # ∫N(D) v(D) dD v_n = QGK.quadgk( D -> - N′ice(p3, D, λ, N₀) * ice_particle_terminal_velocity( + N′ice(p3, D, λ, N₀) * p3_particle_terminal_velocity( p3, D, Chen2022, ρₐ, - F_r, + F_rim, + F_liq, th, use_aspect_ratio, ), diff --git a/test/gpu_tests.jl b/test/gpu_tests.jl index ba0717f5c..99d2b8574 100644 --- a/test/gpu_tests.jl +++ b/test/gpu_tests.jl @@ -111,6 +111,9 @@ end qᵥ_sl, qᵢ, qᵢ_s, + w, + p, + const_dt, ) where {FT} i = @index(Group, Linear) @@ -123,7 +126,19 @@ end ρ[i], T[i], ) - output[2, i] = CMN.conv_q_vap_to_q_liq_ice( + output[2, i] = CMN.conv_q_vap_to_q_liq_ice_MM2015_timeintegrator( + liquid, + ice, + tps, + TD.PhasePartition(qᵥ_sl[i]), + ρ[i], + T[i], + w[i], + p[i], + const_dt[i], + Val(:condensation), + ) + output[3, i] = CMN.conv_q_vap_to_q_liq_ice( ice, TD.PhasePartition(FT(0), FT(0), qᵢ_s[i]), TD.PhasePartition(FT(0), FT(0), qᵢ[i]), @@ -706,15 +721,15 @@ end @kernel function P3_scheme_kernel!( p3, output::AbstractArray{FT}, - F_r, + F_rim, ρ_r, ) where {FT} i = @index(Group, Linear) @inbounds begin - output[1, i] = P3.thresholds(p3, ρ_r[i], F_r[i])[1] - output[2, i] = P3.thresholds(p3, ρ_r[i], F_r[i])[2] + output[1, i] = P3.thresholds(p3, ρ_r[i], F_rim[i])[1] + output[2, i] = P3.thresholds(p3, ρ_r[i], F_rim[i])[2] end end @@ -812,7 +827,8 @@ function test_gpu(FT) end @testset "non-equilibrium microphysics kernels" begin - dims = (2, 1) + + dims = (3, 1) (; output, ndrange) = setup_output(dims, FT) ρ = ArrayType([FT(0.8)]) @@ -820,13 +836,34 @@ function test_gpu(FT) qᵥ_sl = ArrayType([FT(0.0035)]) qᵢ = ArrayType([FT(0.003)]) qᵢ_s = ArrayType([FT(0.002)]) + w = ArrayType([FT(1)]) + p = ArrayType([FT(800 * 1e2)]) + const_dt = ArrayType([FT(0.1)]) + w = ArrayType([FT(1)]) + p = ArrayType([FT(800 * 1e2)]) + const_dt = ArrayType([FT(0.1)]) kernel! = test_noneq_micro_kernel!(backend, work_groups) - kernel!(liquid, ice, tps, output, ρ, T, qᵥ_sl, qᵢ, qᵢ_s, ; ndrange) + kernel!( + liquid, + ice, + tps, + output, + ρ, + T, + qᵥ_sl, + qᵢ, + qᵢ_s, + w, + p, + const_dt; + ndrange, + ) # test that nonequilibrium cloud formation is callable and returns a reasonable value - @test Array(output)[1] ≈ FT(9.043587231238157e-5) + @test Array(output)[1] ≈ FT(3.763783850665844e-5) @test Array(output)[2] ≈ FT(-1e-4) + @test Array(output)[3] ≈ FT(-1e-4) end @testset "0-moment microphysics kernels" begin @@ -1306,11 +1343,11 @@ function test_gpu(FT) dims = (2, 2) (; output, ndrange) = setup_output(dims, FT) - F_r = ArrayType([FT(0.5), FT(0.95)]) + F_rim = ArrayType([FT(0.5), FT(0.95)]) ρ_r = ArrayType([FT(400), FT(800)]) kernel! = P3_scheme_kernel!(backend, work_groups) - kernel!(p3, output, F_r, ρ_r; ndrange) + kernel!(p3, output, F_rim, ρ_r; ndrange) # test if all output is positive... @test all(Array(output) .> FT(0)) diff --git a/test/microphysics_noneq_tests.jl b/test/microphysics_noneq_tests.jl index 712e130d1..dc91b5d41 100644 --- a/test/microphysics_noneq_tests.jl +++ b/test/microphysics_noneq_tests.jl @@ -54,7 +54,7 @@ function test_microphysics_noneq(FT) end end - TT.@testset "CondEvap_DepSub_MG2008" begin + TT.@testset "CondEvap_DepSub_MM2015" begin ρ = FT(0.8) T = FT(273 - 10) @@ -78,8 +78,8 @@ function test_microphysics_noneq(FT) TT.@test CMNe.conv_q_vap_to_q_liq_ice_MM2015(ice, tps, qₚ( qᵥ_si), ρ, T) ≈ FT(0) # smoke test for values - TT.@test CMNe.conv_q_vap_to_q_liq_ice_MM2015(liquid, tps, qₚ(FT(1.2 * qᵥ_sl)), ρ, T) ≈ 9.0419475e-5 rtol = 1e-6 - TT.@test CMNe.conv_q_vap_to_q_liq_ice_MM2015(ice, tps, qₚ(FT(1.2 * qᵥ_si)), ρ, T) ≈ 8.627814e-5 rtol = 1e-6 + TT.@test CMNe.conv_q_vap_to_q_liq_ice_MM2015(liquid, tps, qₚ(FT(1.2 * qᵥ_sl)), ρ, T) ≈ 3.7631e-5 rtol = 1e-6 + TT.@test CMNe.conv_q_vap_to_q_liq_ice_MM2015(ice, tps, qₚ(FT(1.2 * qᵥ_si)), ρ, T) ≈ 3.2356777e-5 rtol = 1e-6 # ice grows faster than liquid TT.@test CMNe.conv_q_vap_to_q_liq_ice_MM2015(liquid, tps, qₚ(FT(1.2 * qᵥ_sl)), ρ, T) < @@ -87,6 +87,44 @@ function test_microphysics_noneq(FT) #! format: on end + + TT.@testset "CondEvap_DepSub_MM2015_timeintegrator" begin + + ρ = FT(0.8) + T = FT(273 - 10) + w = FT(1) + p_air = FT(800 * 1e2) + const_dt = FT(0.1) + + pᵥ_sl = TD.saturation_vapor_pressure(tps, T, TD.Liquid()) + qᵥ_sl = TD.q_vap_saturation_from_density(tps, T, ρ, pᵥ_sl) + + pᵥ_si = TD.saturation_vapor_pressure(tps, T, TD.Ice()) + qᵥ_si = TD.q_vap_saturation_from_density(tps, T, ρ, pᵥ_si) + + qₚ(qᵥ) = TD.PhasePartition(FT(qᵥ)) + + #! format: off + # test sign + TT.@test CMNe.conv_q_vap_to_q_liq_ice_MM2015_timeintegrator(liquid, ice, tps, qₚ(FT(0.5 * qᵥ_sl)), ρ, T, w, p_air, const_dt, Val(:condensation)) < FT(0) + TT.@test CMNe.conv_q_vap_to_q_liq_ice_MM2015_timeintegrator(liquid, ice, tps, qₚ(FT(1.5 * qᵥ_sl)), ρ, T, w, p_air, const_dt, Val(:condensation)) > FT(0) + # wouldnt be zero as a result of WBF + #TT.@test CMNe.conv_q_vap_to_q_liq_ice_MM2015_timeintegrator(liquid, ice, tps, qₚ( qᵥ_sl), ρ, T, FT(0), p_air, const_dt, "condensation") ≈ FT(0) rtol = 1e-6 + + TT.@test CMNe.conv_q_vap_to_q_liq_ice_MM2015_timeintegrator(liquid, ice, tps, qₚ(FT(0.5 * qᵥ_si)), ρ, T, w, p_air, const_dt, Val(:deposition)) < FT(0) + TT.@test CMNe.conv_q_vap_to_q_liq_ice_MM2015_timeintegrator(liquid, ice, tps, qₚ(FT(1.5 * qᵥ_si)), ρ, T, w, p_air, const_dt, Val(:deposition)) > FT(0) + #TT.@test CMNe.conv_q_vap_to_q_liq_ice_MM2015_timeintegrator(liquid, ice, tps, qₚ( qᵥ_si), ρ, T, FT(0), p_air, const_dt, "deposition") ≈ FT(0) rtol = 1e-6 + + # smoke test for values + TT.@test CMNe.conv_q_vap_to_q_liq_ice_MM2015_timeintegrator(liquid, ice, tps, qₚ(FT(1.2 * qᵥ_sl)), ρ, T, w, p_air, const_dt, Val(:condensation)) ≈ 3.7177455e-5 rtol = 5e-6 + TT.@test CMNe.conv_q_vap_to_q_liq_ice_MM2015_timeintegrator(liquid, ice, tps, qₚ(FT(1.2 * qᵥ_si)), ρ, T, w, p_air, const_dt, Val(:deposition)) ≈ 3.2125972e-5 rtol = 1e-6 + + # ice grows faster than liquid + TT.@test CMNe.conv_q_vap_to_q_liq_ice_MM2015_timeintegrator(liquid, ice, tps, qₚ(FT(1.2 * qᵥ_sl)), ρ, T, w, p_air, const_dt, Val(:condensation)) < + CMNe.conv_q_vap_to_q_liq_ice_MM2015_timeintegrator(liquid, ice, tps, qₚ(FT(1.2 * qᵥ_sl)), ρ, T, w, p_air, const_dt, Val(:deposition)) + + #! format: on + end end test_microphysics_noneq(Float32) diff --git a/test/p3_tests.jl b/test/p3_tests.jl index d52188ed2..8d9b02538 100644 --- a/test/p3_tests.jl +++ b/test/p3_tests.jl @@ -8,74 +8,97 @@ import Thermodynamics as TD import QuadGK as QGK -@info "P3 Scheme Tests" +@info("P3 Scheme Tests") -function test_p3_thresholds(FT) +function test_incomplete_gamma_function_approximation(FT) + # A smoke tests to see if the reference values are not changed + TT.@testset "Incomplete Γ function approximation" begin + + Γ_lower_reference = [ + [0.6171173028812539, 0.9999647375921944, 0.9999999999999561], + [66.43238525128415, 192195.93434392574, 362553.66f0], + [9.801186821098075e11, 8.309791f14, 1.1845235f17], + ] + + a = [FT(1), FT(10), FT(20)] + z = [FT(1), FT(10), FT(30)] + + for it in [1, 2, 3] + for jt in [1, 2, 3] + TT.@test P3.Γ_lower(a[it], z[jt]) ≈ Γ_lower_reference[it][jt] rtol = + 1e-3 + end + end + end +end + +function test_thresholds_solver(FT) p3 = CMP.ParametersP3(FT) - TT.@testset "thresholds (nonlinear solver function)" begin + TT.@testset "Thresholds - nonlinear solver" begin # initialize test values: ρ_r = FT(400) - F_r = FT(0.8) + F_rim = FT(0.8) ρ_r_good = (FT(200), FT(400), FT(800)) # representative ρ_r values - F_r_good = (FT(0.5), FT(0.8), FT(0.95)) # representative F_r values + F_rim_good = (FT(0.5), FT(0.8), FT(0.95)) # representative F_rim values # test asserts for _ρ_r in (FT(0), FT(-1)) TT.@test_throws AssertionError("ρ_r > FT(0)") P3.thresholds( p3, _ρ_r, - F_r, + F_rim, ) end for _ρ_r in (FT(0), FT(-1)) TT.@test_throws AssertionError("ρ_r <= p3.ρ_l") P3.thresholds( p3, FT(1200), - F_r, + F_rim, ) end - for _F_r in (FT(-1 * eps(FT)), FT(-1)) - TT.@test_throws AssertionError("F_r >= FT(0)") P3.thresholds( + for _F_rim in (FT(-1 * eps(FT)), FT(-1)) + TT.@test_throws AssertionError("F_rim >= FT(0)") P3.thresholds( p3, ρ_r, - _F_r, + _F_rim, ) end - for _F_r in (FT(1), FT(1.5)) - TT.@test_throws AssertionError("F_r < FT(1)") P3.thresholds( + for _F_rim in (FT(1), FT(1.5)) + TT.@test_throws AssertionError("F_rim < FT(1)") P3.thresholds( p3, ρ_r, - _F_r, + _F_rim, ) end # Test if the P3 scheme solution satisifies the conditions # from eqs. 14-17 in Morrison and Milbrandt 2015 - for F_r in F_r_good + for F_rim in F_rim_good for ρ_r in ρ_r_good - sol = P3.thresholds(p3, ρ_r, F_r) + sol = P3.thresholds(p3, ρ_r, F_rim) atol = 5e3 * eps(FT) - TT.@test sol.D_cr ≈ P3.D_cr_helper(p3, F_r, sol[3]) atol = atol + TT.@test sol.D_cr ≈ P3.D_cr_helper(p3, F_rim, sol[3]) atol = + atol TT.@test sol.D_gr ≈ P3.D_gr_helper(p3, sol[3]) atol = atol - TT.@test sol.ρ_g ≈ P3.ρ_g_helper(ρ_r, F_r, sol[4]) atol = atol + TT.@test sol.ρ_g ≈ P3.ρ_g_helper(ρ_r, F_rim, sol[4]) atol = atol TT.@test sol.ρ_d ≈ P3.ρ_d_helper(p3, sol[1], sol[2]) atol = atol end end # Check that the P3 scheme solution matches the published values - function diff(ρ_r, F_r, el, gold, rtol = 2e-2) - TT.@test P3.thresholds(p3, ρ_r, F_r)[el] * 1e3 ≈ gold rtol = rtol + function diff(ρ_r, F_rim, el, gold, rtol = 2e-2) + TT.@test P3.thresholds(p3, ρ_r, F_rim)[el] * 1e3 ≈ gold rtol = rtol end # D_cr and D_gr vs Fig. 1a Morrison and Milbrandt 2015 D_cr_fig_1a_ref = [FT(0.4946323381999426), FT(1.0170979628696817)] D_gr_fig_1a_ref = [FT(0.26151186272014415), FT(0.23392868352755775)] for val in [1, 2] - diff(ρ_r_good[2], F_r_good[val], 1, D_cr_fig_1a_ref[val]) - diff(ρ_r_good[2], F_r_good[val], 2, D_gr_fig_1a_ref[val]) + diff(ρ_r_good[2], F_rim_good[val], 1, D_cr_fig_1a_ref[val]) + diff(ρ_r_good[2], F_rim_good[val], 2, D_gr_fig_1a_ref[val]) end # D_cr and D_gr vs Fig. 1b Morrison and Milbrandt 2015 #! format: off @@ -83,19 +106,20 @@ function test_p3_thresholds(FT) D_gr_fig_1b_ref = [FT(0.39875043123651077), FT(0.2147085163169669), FT(0.11516682512848)] #! format: on for val in [1, 2, 3] - diff(ρ_r_good[val], F_r_good[3], 1, D_cr_fig_1b_ref[val]) - diff(ρ_r_good[val], F_r_good[3], 2, D_gr_fig_1b_ref[val]) + diff(ρ_r_good[val], F_rim_good[3], 1, D_cr_fig_1b_ref[val]) + diff(ρ_r_good[val], F_rim_good[3], 2, D_gr_fig_1b_ref[val]) end end - TT.@testset "mass, area, density and aspect ratio tests" begin + TT.@testset "Thresholds - mass, area, density, aspect ratio" begin # values ρ_r = FT(500) - F_r = FT(0.5) + F_rim = FT(0.5) + F_liq = FT(0) # testing F_liq = 0 # get thresholds D_th = P3.D_th_helper(p3) - th = P3.thresholds(p3, ρ_r, F_r) + th = P3.thresholds(p3, ρ_r, F_rim) (; D_gr, D_cr) = th # define in between values @@ -104,78 +128,134 @@ function test_p3_thresholds(FT) D_3 = (D_gr + D_cr) / 2 # test area - TT.@test P3.p3_area(p3, D_1, F_r, th) == P3.A_s(D_1) - TT.@test P3.p3_area(p3, D_2, F_r, th) == P3.A_ns(p3, D_2) - TT.@test P3.p3_area(p3, D_3, F_r, th) == P3.A_s(D_3) - TT.@test P3.p3_area(p3, D_cr, F_r, th) == P3.A_r(p3, F_r, D_cr) + TT.@test P3.p3_area(p3, D_1, F_rim, F_liq, th) == P3.A_s(D_1) + TT.@test P3.p3_area(p3, D_2, F_rim, F_liq, th) == P3.A_ns(p3, D_2) + TT.@test P3.p3_area(p3, D_3, F_rim, F_liq, th) == P3.A_s(D_3) + TT.@test P3.p3_area(p3, D_cr, F_rim, F_liq, th) == + P3.A_r(p3, F_rim, D_cr) # test mass - TT.@test P3.p3_mass(p3, D_1, F_r, th) == P3.mass_s(D_1, p3.ρ_i) - TT.@test P3.p3_mass(p3, D_2, F_r, th) == P3.mass_nl(p3, D_2) - TT.@test P3.p3_mass(p3, D_3, F_r, th) == P3.mass_s(D_3, th.ρ_g) - TT.@test P3.p3_mass(p3, D_cr, F_r, th) == P3.mass_r(p3, D_cr, F_r) + TT.@test P3.p3_mass(p3, D_1, F_rim, F_liq, th) == P3.mass_s(D_1, p3.ρ_i) + TT.@test P3.p3_mass(p3, D_2, F_rim, F_liq, th) == P3.mass_nl(p3, D_2) + TT.@test P3.p3_mass(p3, D_3, F_rim, F_liq, th) == P3.mass_s(D_3, th.ρ_g) + TT.@test P3.p3_mass(p3, D_cr, F_rim, F_liq, th) == + P3.mass_r(p3, D_cr, F_rim) # test density - TT.@test P3.p3_density(p3, D_1, F_r, th) ≈ p3.ρ_i - TT.@test P3.p3_density(p3, D_2, F_r, th) ≈ 544.916989830 - TT.@test P3.p3_density(p3, D_3, F_r, th) ≈ th.ρ_g - TT.@test P3.p3_density(p3, D_cr, F_r, th) ≈ 383.33480937 + TT.@test P3.p3_density(p3, D_1, F_rim, th) ≈ p3.ρ_i + TT.@test P3.p3_density(p3, D_2, F_rim, th) ≈ 544.916989830 + TT.@test P3.p3_density(p3, D_3, F_rim, th) ≈ th.ρ_g + TT.@test P3.p3_density(p3, D_cr, F_rim, th) ≈ 383.33480937 # test aspect ratio - TT.@test P3.ϕᵢ(p3, D_1, F_r, th) ≈ 1 - TT.@test P3.ϕᵢ(p3, D_2, F_r, th) ≈ 0.5777887690 - TT.@test P3.ϕᵢ(p3, D_3, F_r, th) ≈ 1 - TT.@test P3.ϕᵢ(p3, D_cr, F_r, th) ≈ 0.662104776 - - # test F_r = 0 and D > D_th - F_r = FT(0) - TT.@test P3.p3_area(p3, D_2, F_r, th) == P3.A_ns(p3, D_2) - TT.@test P3.p3_mass(p3, D_2, F_r, th) == P3.mass_nl(p3, D_2) + TT.@test P3.ϕᵢ(p3, D_1, F_rim, th) ≈ 1 + TT.@test P3.ϕᵢ(p3, D_2, F_rim, th) ≈ 0.5777887690 + TT.@test P3.ϕᵢ(p3, D_3, F_rim, th) ≈ 1 + TT.@test P3.ϕᵢ(p3, D_cr, F_rim, th) ≈ 0.662104776 + + # test F_rim = 0 and D > D_th + F_rim = FT(0) + TT.@test P3.p3_area(p3, D_2, F_rim, F_liq, th) == P3.A_ns(p3, D_2) + TT.@test P3.p3_mass(p3, D_2, F_rim, F_liq, th) == P3.mass_nl(p3, D_2) + + # test F_liq != 0 + F_liq = FT(0.5) + + # get D_i values again + # values + ρ_r = FT(500) + F_rim = FT(0.5) + + # get thresholds + D_th = P3.D_th_helper(p3) + th = P3.thresholds(p3, ρ_r, F_rim) + (; D_gr, D_cr) = th + + # define in between values + D_1 = D_th / 2 + D_2 = (D_th + D_gr) / 2 + D_3 = (D_gr + D_cr) / 2 + + # test area + TT.@test P3.p3_area(p3, D_1, F_rim, F_liq, th) == P3.A_s(D_1) + TT.@test P3.p3_area(p3, D_2, F_rim, F_liq, th) == + (1 - F_liq) * P3.A_ns(p3, D_2) + F_liq * P3.A_s(D_2) + TT.@test P3.p3_area(p3, D_3, F_rim, F_liq, th) == P3.A_s(D_3) + TT.@test P3.p3_area(p3, D_cr, F_rim, F_liq, th) == + (1 - F_liq) * P3.A_r(p3, F_rim, D_cr) + F_liq * P3.A_s(D_cr) + + # test mass + TT.@test P3.p3_mass(p3, D_1, F_rim, F_liq, th) == + (1 - F_liq) * P3.mass_s(D_1, p3.ρ_i) + + F_liq * P3.mass_s(D_1, p3.ρ_l) + TT.@test P3.p3_mass(p3, D_2, F_rim, F_liq, th) == + (1 - F_liq) * P3.mass_nl(p3, D_2) + + F_liq * P3.mass_s(D_2, p3.ρ_l) + TT.@test P3.p3_mass(p3, D_3, F_rim, F_liq, th) == + (1 - F_liq) * P3.mass_s(D_3, th.ρ_g) + + F_liq * P3.mass_s(D_3, p3.ρ_l) + TT.@test P3.p3_mass(p3, D_cr, F_rim, F_liq, th) == + (1 - F_liq) * P3.mass_r(p3, D_cr, F_rim) + + F_liq * P3.mass_s(D_cr, p3.ρ_l) + + # test F_rim = 0 and D > D_th + F_rim = FT(0) + TT.@test P3.p3_area(p3, D_2, F_rim, F_liq, th) == + (1 - F_liq) * P3.A_ns(p3, D_2) + F_liq * P3.A_s(D_2) + TT.@test P3.p3_mass(p3, D_2, F_rim, F_liq, th) == + (1 - F_liq) * P3.mass_nl(p3, D_2) + + F_liq * P3.mass_s(D_2, p3.ρ_l) end end -function test_p3_shape_solver(FT) +function test_shape_solver(FT) p3 = CMP.ParametersP3(FT) - TT.@testset "shape parameters (nonlinear solver function)" begin + TT.@testset "Shape parameters - nonlinear solver" begin # initialize test values: ep = 1 #1e4 * eps(FT) - N_test = (FT(1e7), FT(1e8), FT(1e9), FT(1e10)) # N values - λ_test = (FT(1e1), FT(1e2), FT(1e3), FT(1e4), FT(1e5), FT(1e6)) # test λ values in range also do 15000, 20000 - ρ_r_test = (FT(200), FT(400), FT(600), FT(800)) # representative ρ_r values - F_r_test = (FT(0), FT(0.5), FT(0.8), FT(0.95)) # representative F_r values + N_test = (FT(1e7), FT(1e8), FT(1e9), FT(1e10)) # N values + λ_test = (FT(1e1), FT(1e2), FT(1e3), FT(1e4), FT(1e5), FT(1e6)) # test λ values in range also do 15000, 20000 + ρ_r_test = (FT(200), FT(400), FT(600), FT(800)) # representative ρ_r values + F_rim_test = (FT(0), FT(0.5), FT(0.8), FT(0.95)) # representative F_rim values + F_liq_test = (FT(0), FT(0.33), FT(0.67), FT(1)) # representative F_rim values # check that the shape solution solves to give correct values for N in N_test for λ_ex in λ_test for ρ_r in ρ_r_test - for F_r in F_r_test - # Compute the shape parameters that correspond to the - # input test values - μ_ex = P3.DSD_μ(p3, λ_ex) - N₀_ex = P3.DSD_N₀(p3, N, λ_ex) - # Find the P3 scheme thresholds - th = P3.thresholds(p3, ρ_r, F_r) - # Convert λ to ensure it remains positive - x = log(λ_ex) - # Compute mass density based on input shape parameters - L_calc = N * P3.L_over_N_gamma(p3, F_r, x, μ_ex, th) - - if L_calc < FT(1) - # Solve for shape parameters - (λ, N₀) = P3.distribution_parameter_solver( - p3, - L_calc, - N, - ρ_r, - F_r, - ) - - # Compare solved values with the input expected values - TT.@test λ ≈ λ_ex rtol = ep - TT.@test N₀ ≈ N₀_ex rtol = ep + for F_rim in F_rim_test + for F_liq in F_liq_test + # Compute the shape parameters that correspond to the + # input test values + μ_ex = P3.DSD_μ(p3, λ_ex) + N₀_ex = P3.DSD_N₀(p3, N, λ_ex) + # Find the P3 scheme thresholds + th = P3.thresholds(p3, ρ_r, F_rim) + # Convert λ to ensure it remains positive + x = log(λ_ex) + # Compute mass density based on input shape parameters + L_calc = + N * + P3.L_over_N_gamma(p3, F_liq, F_rim, x, μ_ex, th) + + if L_calc < FT(1) + # Solve for shape parameters + (λ, N₀) = P3.distribution_parameter_solver( + p3, + L_calc, + N, + ρ_r, + F_liq, + F_rim, + ) + + # Compare solved values with the input expected values + TT.@test λ ≈ λ_ex rtol = ep + TT.@test N₀ ≈ N₀_ex rtol = ep + end end end end @@ -190,7 +270,7 @@ function test_particle_terminal_velocities(FT) Chen2022 = CMP.Chen2022VelType(FT) ρ_a = FT(1.2) - TT.@testset "Chen 2022 - Rain" begin + TT.@testset "Smoke tests for rain particle terminal vel from Chen 2022" begin Ds = range(FT(1e-6), stop = FT(1e-5), length = 5) expected = [0.002508, 0.009156, 0.01632, 0.02377, 0.03144] for i in axes(Ds, 1) @@ -200,10 +280,11 @@ function test_particle_terminal_velocities(FT) end end - TT.@testset "Chen 2022 - Ice" begin - F_r = FT(0.5) + TT.@testset "Smoke tests for ice particle terminal vel from Chen 2022" begin + F_rim = FT(0.5) ρ_r = FT(500) - th = P3.thresholds(p3, ρ_r, F_r) + F_liq = FT(0) + th = P3.thresholds(p3, ρ_r, F_rim) use_aspect_ratio = false # Allow for a D falling into every regime of the P3 Scheme Ds = range(FT(0.5e-4), stop = FT(4.5e-4), length = 5) @@ -215,7 +296,71 @@ function test_particle_terminal_velocities(FT) D, Chen2022.snow_ice, ρ_a, - F_r, + F_rim, + th, + use_aspect_ratio, + ) + TT.@test vel >= 0 + TT.@test vel ≈ expected[i] rtol = 1e-3 + end + + use_aspect_ratio = true + # Allow for a D falling into every regime of the P3 Scheme + Ds = range(FT(0.5e-4), stop = FT(4.5e-4), length = 5) + expected = [0.08109, 0.34441, 0.79121, 1.155, 1.289] + for i in axes(Ds, 1) + D = Ds[i] + vel = P3.ice_particle_terminal_velocity( + p3, + D, + Chen2022.snow_ice, + ρ_a, + F_rim, + th, + use_aspect_ratio, + ) + TT.@test vel >= 0 + TT.@test vel ≈ expected[i] rtol = 1e-3 + end + end + + TT.@testset "Smoke tests for mixed phase particle terminal velocity" begin + F_rim = FT(0.5) + F_liq = FT(0.5) + ρ_r = FT(500) + th = P3.thresholds(p3, ρ_r, F_rim) + use_aspect_ratio = true + # Allow for a D falling into every regime of the P3 Scheme + Ds = range(FT(0.5e-4), stop = FT(4.5e-4), length = 5) + expected = [0.13192, 0.471, 0.90753, 1.3015, 1.5766] + for i in axes(Ds, 1) + D = Ds[i] + vel = P3.p3_particle_terminal_velocity( + p3, + D, + Chen2022, + ρ_a, + F_rim, + F_liq, + th, + use_aspect_ratio, + ) + TT.@test vel >= 0 + TT.@test vel ≈ expected[i] rtol = 1e-3 + end + use_aspect_ratio = false + # Allow for a D falling into every regime of the P3 Scheme + Ds = range(FT(0.5e-4), stop = FT(4.5e-4), length = 5) + expected = [0.13191, 0.50457, 0.90753, 1.301499, 1.67569] + for i in axes(Ds, 1) + D = Ds[i] + vel = P3.p3_particle_terminal_velocity( + p3, + D, + Chen2022, + ρ_a, + F_rim, + F_liq, th, use_aspect_ratio, ) @@ -231,98 +376,206 @@ function test_bulk_terminal_velocities(FT) L = FT(0.22) N = FT(1e6) ρ_a = FT(1.2) - ρ_rs = [FT(200), FT(400), FT(600), FT(800)] - F_rs = [FT(0), FT(0.2), FT(0.4), FT(0.6), FT(0.8)] + ρ_r = FT(800) + F_rims = [FT(0), FT(0.6)] + F_liqs = [FT(0.5), FT(1)] TT.@testset "Mass and number weighted terminal velocities" begin - reference_vals_m = [ - [7.79, 7.27, 6.66, 5.94, 5.25], - [7.79, 7.26, 6.62, 5.83, 4.82], - [7.79, 7.25, 6.62, 5.81, 4.7], - [7.79, 7.25, 6.62, 5.81, 4.65], - ] - reference_vals_n = [ - [3.65, 3.37, 3.05, 2.64, 2.14], - [3.64, 3.37, 3.04, 2.62, 2.04], - [3.65, 3.37, 3.04, 2.62, 2.02], - [3.64, 3.37, 3.04, 2.61, 2.01], - ] - reference_vals_m_ϕ = [ - [4.23, 4.65, 4.90, 4.94, 4.89], - [4.23, 4.65, 4.88, 4.84, 4.44], - [4.23, 4.65, 4.87, 4.82, 4.33], - [4.23, 4.64, 4.87, 4.81, 4.28], - ] - reference_vals_n_ϕ = [ - [2.21, 2.34, 2.38, 2.31, 2.04], - [2.21, 2.33, 2.37, 2.27, 2.04], - [2.21, 2.33, 2.35, 2.25, 1.91], - [2.21, 2.33, 2.36, 2.24, 1.9], - ] - for i in 1:length(ρ_rs) - for j in 1:length(F_rs) - ρ_r = ρ_rs[i] - F_r = F_rs[j] + # Liquid fraction = 0 + + ref_v_n = [3.637320375996216, 2.6130243256593313] + ref_v_n_ϕ = [2.2058692424028923, 2.24224820764301] + ref_v_m = [7.7208191246659705, 5.751352472270078] + ref_v_m_ϕ = [4.193799087495265, 4.773592179767702] + + for k in 1:length(F_rims) + F_rim = F_rims[k] + F_liq = FT(0) + vel = P3.ice_terminal_velocity( + p3, + Chen2022, + L, + N, + ρ_r, + F_rim, + F_liq, + ρ_a, + false, + ) + vel_ϕ = P3.ice_terminal_velocity( + p3, + Chen2022, + L, + N, + ρ_r, + F_rim, + F_liq, + ρ_a, + true, + ) + + # number weighted + TT.@test vel[1] > 0 + TT.@test vel_ϕ[1] > 0 + TT.@test vel[1] ≈ ref_v_n[k] rtol = 1e-6 + TT.@test vel_ϕ[1] ≈ ref_v_n_ϕ[k] rtol = 1e-6 + + # mass weighted + TT.@test vel[2] > 0 + TT.@test vel_ϕ[2] > 0 + TT.@test vel[2] ≈ ref_v_m[k] rtol = 1e-6 + TT.@test vel_ϕ[2] ≈ ref_v_m_ϕ[k] rtol = 1e-6 + + # slower with aspect ratio + TT.@test vel_ϕ[1] <= vel[1] + TT.@test vel_ϕ[2] <= vel[2] + end + + # Liquid fraction != 0 + ref_v_n = [1.674591925057434, 1.6180970319460353] + ref_v_n_ϕ = [1.549777478756061, 1.6180970319460353] + ref_v_m = [5.126648558302173, 5.416679316254198] + ref_v_m_ϕ = [4.6358422594886495, 5.416679316254198] + + for k in 1:length(F_liqs) + F_liq = F_liqs[k] + F_rim = FT(0.4) + vel = P3.ice_terminal_velocity( + p3, + Chen2022, + L, + N, + ρ_r, + F_rim, + F_liq, + ρ_a, + false, + ) + vel_ϕ = P3.ice_terminal_velocity( + p3, + Chen2022, + L, + N, + ρ_r, + F_rim, + F_liq, + ρ_a, + true, + ) + # number weighted + TT.@test vel[1] > 0 + TT.@test vel_ϕ[1] > 0 + TT.@test vel[1] ≈ ref_v_n[k] rtol = 1e-6 + TT.@test vel_ϕ[1] ≈ ref_v_n_ϕ[k] rtol = 1e-6 + + # mass weighted + TT.@test vel[2] > 0 + TT.@test vel_ϕ[2] > 0 + TT.@test vel[2] ≈ ref_v_m[k] rtol = 1e-6 + TT.@test vel_ϕ[2] ≈ ref_v_m_ϕ[k] rtol = 1e-6 + + # slower with aspect ratio + TT.@test vel_ϕ[1] <= vel[1] + TT.@test vel_ϕ[2] <= vel[2] + end + end + TT.@testset "Mass-weighted mean diameters" begin + F_liq = FT(0) + ref_vals = [0.00536934536778844, 0.003321590762128599] + for j in 1:length(F_rims) + Dₘ = P3.D_m(p3, L, N, ρ_r, F_rims[j], F_liq) + TT.@test Dₘ > 0 + TT.@test Dₘ ≈ ref_vals[j] + end + + # nonzero F_liq + F_rim = F_rims[2] + F_liqs = [FT(0.33), FT(1)] + ref_vals = [FT(0.0021371920600012184), FT(0.0016487352655895715)] + for i in eachindex(F_liqs) + Dₘ = P3.D_m(p3, L, N, ρ_r, F_rim, F_liqs[i]) + TT.@test Dₘ ≈ ref_vals[i] + end + end +end + +function test_numerical_integrals(FT) + p3 = CMP.ParametersP3(FT) + Chen2022 = CMP.Chen2022VelType(FT) + + N = FT(1e8) + Ls = range(FT(0.001), stop = FT(0.005), length = 5) + ρ_r = FT(500) + F_rims = [FT(0), FT(0.5)] + ρ_a = FT(1.2) + use_aspect_ratio = false + tolerance = eps(FT) + F_liq = FT(0) # test F_liq = 0 + + TT.@testset "Numerical integrals sanity checks for N, velocity and diameter" begin + for F_rim in F_rims + for i in axes(Ls, 1) + L = Ls[i] - calculated_vel = P3.ice_terminal_velocity( + # Get shape parameters, thresholds and intergal bounds + λ, N_0 = P3.distribution_parameter_solver( p3, - Chen2022.snow_ice, L, N, ρ_r, - F_r, - ρ_a, - false, + F_rim, + F_liq, ) - calculated_vel_ϕ = P3.ice_terminal_velocity( + th = P3.thresholds(p3, ρ_r, F_rim) + ice_bound = P3.get_ice_bound(p3, λ, N, tolerance) + + # Number concentration comparison + f(d) = P3.N′ice(p3, d, λ, N_0) + qgk_N, = QGK.quadgk(d -> f(d), FT(0), ice_bound) + TT.@test N ≈ qgk_N rtol = sqrt(eps(FT)) + + # Bulk velocity comparison + vel_N, vel_m = P3.ice_terminal_velocity( p3, - Chen2022.snow_ice, + Chen2022, L, N, ρ_r, - F_r, + F_rim, + F_liq, ρ_a, - true, + use_aspect_ratio, ) - # number weighted - TT.@test calculated_vel[1] > 0 - TT.@test reference_vals_n[i][j] ≈ calculated_vel[1] atol = 0.1 - TT.@test calculated_vel_ϕ[1] > 0 - TT.@test reference_vals_n_ϕ[i][j] ≈ calculated_vel_ϕ[1] atol = - 0.1 - - # mass weighted - TT.@test calculated_vel[2] > 0 - TT.@test reference_vals_m[i][j] ≈ calculated_vel[2] atol = 0.1 - TT.@test calculated_vel_ϕ[2] > 0 - TT.@test reference_vals_m_ϕ[i][j] ≈ calculated_vel_ϕ[2] atol = - 0.1 - - TT.@test calculated_vel_ϕ[1] <= calculated_vel[1] - TT.@test calculated_vel_ϕ[2] <= calculated_vel[2] - end - end - end - - TT.@testset "Mass-weighted mean diameters" begin - paper_vals = [ - [5, 5, 5, 5, 5], - [4.5, 4.5, 4.5, 4.5, 4.5], - [3.5, 3.5, 3.5, 3.5, 3.5], - [3.5, 3.5, 2.5, 2.5, 2.5], - ] - for i in 1:length(ρ_rs) - for j in 1:length(F_rs) - ρ_r = ρ_rs[i] - F_r = F_rs[j] - - calculated_dm = P3.D_m(p3, L, N, ρ_r, F_r) * 1e3 + particle_vel(d) = P3.ice_particle_terminal_velocity( + p3, + d, + Chen2022.snow_ice, + ρ_a, + F_rim, + th, + use_aspect_ratio, + ) + g(d) = particle_vel(d) * P3.N′ice(p3, d, λ, N_0) - TT.@test calculated_dm > 0 - TT.@test paper_vals[i][j] ≈ calculated_dm atol = 3.14 + qgk_vel_N, = QGK.quadgk(d -> g(d) / N, FT(0), ice_bound) + qgk_vel_m, = QGK.quadgk( + d -> g(d) * P3.p3_mass(p3, d, F_rim, F_liq, th) / L, + FT(0), + ice_bound, + ) + TT.@test vel_N ≈ qgk_vel_N rtol = 1e-6 + TT.@test vel_m ≈ qgk_vel_m rtol = 1e-6 + # Dₘ comparisons + D_m = P3.D_m(p3, L, N, ρ_r, F_rim, F_liq) + h(d) = + d * + P3.p3_mass(p3, d, F_rim, F_liq, th) * + P3.N′ice(p3, d, λ, N_0) + qgk_D_m, = QGK.quadgk(d -> h(d) / L, FT(0), ice_bound) + TT.@test D_m ≈ qgk_D_m rtol = 1e-2 end end end @@ -461,71 +714,6 @@ end # #end -function test_integrals(FT) - p3 = CMP.ParametersP3(FT) - Chen2022 = CMP.Chen2022VelType(FT) - - N = FT(1e8) - Ls = range(0.001, stop = 0.005, length = 5) - ρ_r = FT(500) - F_rs = [FT(0), FT(0.5)] - ρ_a = FT(1.2) - use_aspect_ratio = false - tolerance = eps(FT) - - TT.@testset "Gamma vs Integral Comparison" begin - for F_r in F_rs - for i in axes(Ls, 1) - L = Ls[i] - - # Velocity comparisons - vel_N, vel_m = P3.ice_terminal_velocity( - p3, - Chen2022.snow_ice, - L, - N, - ρ_r, - F_r, - ρ_a, - use_aspect_ratio, - ) - - λ, N_0 = P3.distribution_parameter_solver(p3, L, N, ρ_r, F_r) - th = P3.thresholds(p3, ρ_r, F_r) - ice_bound = P3.get_ice_bound(p3, λ, tolerance) - vel(d) = P3.ice_particle_terminal_velocity( - p3, - d, - Chen2022.snow_ice, - ρ_a, - F_r, - th, - use_aspect_ratio, - ) - f(d) = vel(d) * P3.N′ice(p3, d, λ, N_0) - - qgk_vel_N, = QGK.quadgk(d -> f(d) / N, FT(0), ice_bound) - qgk_vel_m, = QGK.quadgk( - d -> f(d) * P3.p3_mass(p3, d, F_r, th) / L, - FT(0), - ice_bound, - ) - - TT.@test vel_N ≈ qgk_vel_N rtol = 1e-5 - TT.@test vel_m ≈ qgk_vel_m rtol = 1e-5 - - # Dₘ comparisons - D_m = P3.D_m(p3, L, N, ρ_r, F_r) - f_d(d) = - d * P3.p3_mass(p3, d, F_r, th) * P3.N′ice(p3, d, λ, N_0) - qgk_D_m, = QGK.quadgk(d -> f_d(d) / L, FT(0), ice_bound) - - TT.@test D_m ≈ qgk_D_m rtol = 1e-8 - end - end - end -end - function test_p3_het_freezing(FT) TT.@testset "Heterogeneous Freezing Smoke Test" begin @@ -544,7 +732,7 @@ function test_p3_het_freezing(FT) expected_freeze_L = [1.544e-22, 1.068e-6, 0.0001428, 0.0001428, 0.0001428, 0.0001428] expected_freeze_N = [1.082e-10, 747647.5, N_liq, N_liq, N_liq, N_liq] - qᵥ_range = range(0.5e-3, stop = 1.5e-3, length = 6) + qᵥ_range = range(FT(0.5e-3), stop = FT(1.5e-3), length = 6) for it in range(1, 6) qₚ = TD.PhasePartition(FT(qᵥ_range[it]), FT(2e-4), FT(0)) @@ -564,19 +752,108 @@ function test_p3_het_freezing(FT) end end -println("Testing Float32") -test_p3_thresholds(Float32) -test_particle_terminal_velocities(Float32) -#TODO - only works for Float64 now. We should switch the units inside the solver -# from SI base to something more managable -#test_p3_shape_solver(Float32) -# test_bulk_terminal_velocities(Float32) - -println("Testing Float64") -test_p3_thresholds(Float64) -test_p3_shape_solver(Float64) -test_particle_terminal_velocities(Float64) -test_bulk_terminal_velocities(Float64) -test_integrals(Float64) -test_p3_het_freezing(Float64) -#test_tendencies(Float64) +function test_p3_melting(FT) + + TT.@testset "Melting Smoke Test" begin + + p3 = CMP.ParametersP3(FT) + vel = CMP.Chen2022VelType(FT) + aps = CMP.AirProperties(FT) + tps = TD.Parameters.ThermodynamicsParameters(FT) + + ρₐ = FT(1.2) + qᵢ = FT(1e-4) + Lᵢ = qᵢ * ρₐ + Nᵢ = FT(2e5) * ρₐ + F_rim = FT(0.8) + ρ_rim = FT(800) + dt = FT(1) + F_liq = FT(0) # process not dependent on F_liq + + T_cold = FT(273.15 - 0.01) + rate = P3.ice_melt( + p3, + vel.snow_ice, + aps, + tps, + Lᵢ, + Nᵢ, + T_cold, + ρₐ, + F_rim, + ρ_rim, + dt, + ) + + TT.@test rate.dNdt == 0 + TT.@test rate.dLdt == 0 + + T_warm = FT(273.15 + 0.01) + rate = P3.ice_melt( + p3, + vel.snow_ice, + aps, + tps, + Lᵢ, + Nᵢ, + T_warm, + ρₐ, + F_rim, + ρ_rim, + dt, + ) + + TT.@test rate.dNdt >= 0 + TT.@test rate.dLdt >= 0 + + TT.@test rate.dNdt ≈ 171326.17992979713 rtol = 1e-6 + TT.@test rate.dLdt ≈ 8.566308996489856e-5 rtol = 1e-6 + + T_vwarm = FT(273.15 + 0.1) + rate = P3.ice_melt( + p3, + vel.snow_ice, + aps, + tps, + Lᵢ, + Nᵢ, + T_vwarm, + ρₐ, + F_rim, + ρ_rim, + dt, + ) + + TT.@test rate.dNdt == Nᵢ + TT.@test rate.dLdt == Lᵢ + end +end + +for FT in [Float32, Float64] + @info("Testing " * string(FT)) + + # numerics + test_incomplete_gamma_function_approximation(FT) + test_thresholds_solver(FT) + + # velocity + test_particle_terminal_velocities(FT) + + # processes + test_p3_het_freezing(FT) +end + +# TODO - make work for Float32 +for FT in [Float64] + @info("Tests that only work for Float64") + + # numerics + test_shape_solver(FT) + test_numerical_integrals(FT) + + # velocity + test_bulk_terminal_velocities(FT) + + # processes + test_p3_melting(FT) +end diff --git a/test/performance_tests.jl b/test/performance_tests.jl index aae5e31ff..b19fbcd1f 100644 --- a/test/performance_tests.jl +++ b/test/performance_tests.jl @@ -114,8 +114,10 @@ function benchmark_test(FT) e = FT(600) ρ_r = FT(400.0) - F_r = FT(0.95) + F_rim = FT(0.95) + F_liq = FT(0.33) N = FT(1e8) + L = ρ_air * q_ice T_air_2 = FT(250) RH_2 = FT(1.5) @@ -149,45 +151,56 @@ function benchmark_test(FT) INPC = FT(1e5) - # P3 scheme - bench_press(P3.thresholds, (p3, ρ_r, F_r), 12e6, 2048, 80) - if FT == Float64 - bench_press( - P3.distribution_parameter_solver, - (p3, q_ice, N, ρ_r, F_r), - 1e5, - ) - bench_press( - P3.ice_terminal_velocity, - (p3, ch2022.snow_ice, q_ice, N, ρ_r, F_r, ρ_air, false), - 2.5e5, - 800, - 4, - ) - bench_press( - P3.ice_terminal_velocity, - (p3, ch2022.snow_ice, q_ice, N, ρ_r, F_r, ρ_air, true), - 2.5e5, - 800, - 4, - ) - bench_press(P3.D_m, (p3, q_ice, N, ρ_r, F_r), 1e5) - end + # P3 scheme - TODO - bring them back once we optimize P3 scheme + #bench_press(P3.thresholds, (p3, ρ_r, F_rim), 12e6, 2048, 80) + #if FT == Float64 + # bench_press( + # P3.distribution_parameter_solver, + # (p3, q_ice, N, ρ_r, F_r), + # 1e5, + # ) + # bench_press( + # P3.ice_terminal_velocity, + # (p3, ch2022.snow_ice, q_ice, N, ρ_r, F_r, ρ_air, false), + # 2.5e5, + # 800, + # 4, + # ) + # bench_press( + # P3.ice_terminal_velocity, + # (p3, ch2022.snow_ice, q_ice, N, ρ_r, F_r, ρ_air, true), + # 2.5e5, + # 800, + # 4, + # ) + # bench_press(P3.D_m, (p3, q_ice, N, ρ_r, F_r), 1e5) + #end # P3 ice nucleation - bench_press( - P3.het_ice_nucleation, - ( - kaolinite, - tps, - TD.PhasePartition(q_tot, q_liq, q_ice), - N_liq, - RH_2, - T_air_2, - ρ_air, - Δt, - ), - 150, - ) + #bench_press( + # P3.het_ice_nucleation, + # ( + # kaolinite, + # tps, + # TD.PhasePartition(q_tot, q_liq, q_ice), + # N_liq, + # RH_2, + # T_air_2, + # ρ_air, + # Δt, + # ), + # 200, + #) + #if FT == Float64 + # bench_press( + # P3.ice_melt, + # (p3, ch2022.snow_ice, aps, tps, L, N, T_air, ρ_air, F_rim, ρ_r, Δt), + # 3.7e5, + # 2e3, + # 3, + # ) + #end + #bench_press(CMI_het.P3_deposition_N_i, (ip.p3, T_air_cold), 230) + #bench_press(CMI_het.P3_het_N_i, (ip.p3, T_air_cold, N_liq, V_liq, Δt), 230) # aerosol activation bench_press( @@ -218,9 +231,7 @@ function benchmark_test(FT) 80, ) bench_press(CMI_het.deposition_J, (kaolinite, Delta_a_w), 230) - bench_press(CMI_het.P3_deposition_N_i, (ip.p3, T_air_cold), 230) bench_press(CMI_het.ABIFM_J, (desert_dust, Delta_a_w), 230) - bench_press(CMI_het.P3_het_N_i, (ip.p3, T_air_cold, N_liq, V_liq, Δt), 230) bench_press( CMI_het.INP_concentration_frequency, (ip_frostenberg, INPC, T_air_cold),