From 6057b0d4643ae9922dade202d7b520343ba04b55 Mon Sep 17 00:00:00 2001 From: Caterina Croci Date: Fri, 24 May 2024 09:16:57 -0700 Subject: [PATCH] Make SB2006 limiters optional, add SB2006 cloud size distribution parameters --- docs/src/API.md | 3 +- docs/src/Microphysics2M.md | 102 ++-- docs/src/plots/Microphysics2M_plots.jl | 20 +- docs/src/plots/RainEvapoartionSB2006.jl | 2 +- src/Microphysics2M.jl | 267 ++++++--- src/PrecipitationSusceptibility.jl | 4 +- src/parameters/Microphysics2M.jl | 61 +- src/parameters/toml/SB2006_limiters.toml | 4 +- test/gpu_tests.jl | 106 ++-- test/microphysics2M_tests.jl | 694 +++++++++++++---------- test/performance_tests.jl | 76 +-- 11 files changed, 835 insertions(+), 504 deletions(-) diff --git a/docs/src/API.md b/docs/src/API.md index d04352f6e..1cbf1f381 100644 --- a/docs/src/API.md +++ b/docs/src/API.md @@ -41,7 +41,8 @@ Microphysics1M.snow_melt ```@docs Microphysics2M Microphysics2M.LiqRaiRates -Microphysics2M.raindrops_limited_vars +Microphysics2M.pdf_cloud +Microphysics2M.pdf_rain Microphysics2M.autoconversion Microphysics2M.accretion Microphysics2M.liquid_self_collection diff --git a/docs/src/Microphysics2M.md b/docs/src/Microphysics2M.md index 5ff159de5..8feb56713 100644 --- a/docs/src/Microphysics2M.md +++ b/docs/src/Microphysics2M.md @@ -40,40 +40,72 @@ where ``x`` and ``y`` are drop masses and ``x^*`` is the mass threshold chosen t The default value of ``x^*=6.54\times 10^{-11} kg`` corresponds to the drop radius ``r^* \approx 25 \mu m``. -The cloud droplets Gamma distribution function is described by +### Assumed size distributions + +The cloud droplet mass follows a Gamma distribution ```math \begin{align} - f_c(x)=A_cx^\nu_c e^{-B_cx},\quad \nu_c=\text{2.0}. + f_c(x) = A_c x^{\nu_c} \exp\left(-B_c x^{\mu_c} \right), \quad \nu_c = 2, \;\; \mu_c=1. \end{align} ``` -Here, - - ``A_c=\frac{N_{liq} \, B_c^{\nu_c+1}}{\Gamma(\nu_c+1)}``, - - ``B_c=\left(\frac{\Gamma(\nu_c+1)}{\Gamma(\nu_c+2)}x_c\right)^{-1}``, - - ``x_c=\frac{L_c}{N_{liq}}`` - averaged mass of cloud particles, - - ``L_c = \rho q_{liq}`` - cloud liquid water content. -The raindrops exponential distribution is expressed as +The rain drop diameter ``D`` follows an exponential distribution +```math +\begin{align} + f_r(D) = N_0 \exp\left(- \lambda D \right), +\end{align} +``` +which is a special case of a Gamma distribution ```math \begin{align} - f_r(D)=\alpha e^{-\beta D}, + f_r(x) = A_r x^{\nu_r} \exp\left(-B_r x^{\mu_r} \right), \quad \nu_r = -\frac{2}{3}, \;\; \mu_r=\frac{1}{3}. \end{align} ``` -where ``D`` is the drop diameter which is proportional to ``x^{1/3}``. -Through this proportionality relation, we can express the raindrops exponential distribution as a gamma distribution function +The free parameters can be found analytically by + integrating over the assumed mass distributions to find the prognostic variables ```math \begin{align} - f_r(x)=A_rx^{-\nu_r}_r e^{-B_rx^{\mu_r}}, + N = \int_0^\infty f(x) dx, \quad q = \frac{1}{\rho_a} \int_0^\infty x \; f(x) dx \end{align} ``` -where ``\nu_r=-\frac{2}{3}``, and ``\mu_r=\frac{1}{3}``. -In this case, - - ``A_r=\frac{\mu_r \, N_{rai}}{\Gamma(\frac{\nu_r+1}{\mu_r})}B_r^{\frac{\nu_r+1}{\mu_r}}``, - - ``B_r=\left(\frac{\Gamma(\frac{\nu_r+1}{\mu_r})}{\Gamma(\frac{\nu_r+2}{\mu_r})}x_r\right)^{(-\mu_r)}``, - - ``x_r=\frac{L_r}{N_{rai}}`` - averaged mass of rain particles, - - ``L_r = \rho q_{rai}`` - rain water content. +where + - ``\rho_a`` is the air density, + - ``\overline{x} = \frac{\rho_a q}{N}`` is the mean droplet mass, + - ``B = \left(\overline{x} \; \frac{\Gamma\left(\frac{\nu+1}{\mu}\right)}{\Gamma\left(\frac{\nu+2}{\mu}\right)}\right)^{-\mu}``, + - ``A = \frac{\mu \, N \, }{\Gamma(\frac{\nu+1}{\mu})} \; B^{\frac{\nu+1}{\mu}}``. + +Similarily, for the exponential rain diameter distribution and assuming spherical rain drops +```math +\begin{align} + N = \int_0^\infty f(D) dD, \quad q = \frac{\pi \rho_w}{6 \rho_a} \int_0^\infty D^3 \; f(D) dD +\end{align} +``` +where + - ``\rho_w`` is the density of water, + - ``\lambda = \left( \frac{\pi \; N_r \; \rho_w}{q_r \rho_a} \right)^{\frac{1}{3}}``, + - ``N_0 = \lambda \; N``. + +Lastly, when computing collision rates between the P3 scheme representation of snow and ice, + it is useful to convert the cloud droplet mass distribution to a diameter distribution. +Again, assuming spherical droplets results in +```math +\begin{align} + f_c(D) = C_c D^{\phi_c} \exp\left(-E_c D^{\psi_c} \right), +\end{align} +``` +where + - ``C_c = \frac{A_c}{3^{\nu_c}} \; \left(\frac{\pi \rho_w}{2}\right)^{\nu_c + 1}`` + - ``E_c = B_c \; \left(\frac{\pi \rho_w}{6}\right)^{\mu_c}`` + - ``\phi_c = 3\nu_c + 2`` + - ``\psi_c = 3\mu_c``. !!! note - In the derivation of the parametrization, it is assumed that the cloud droplet distribution ``f_c(x)`` does not contain a significant number of droplets with masses almost equal or larger than ``x^*``. This is reffered to as the undeveloped cloud droplet spectrum assumption. Similarly the raindrop distribution does not contain a significant number of rain drops with masses almost equal or smaller than ``x^*``. These assumptions allow us to simplify the calculation of moments of the distributions by integrating from zero to infinity. + In the derivation of the parametrization, it is assumed that the cloud droplet distribution + ``f_c(x)`` does not contain a significant number of droplets with masses almost equal + or larger than ``x^*``. This is reffered to as the undeveloped cloud droplet spectrum assumption. + Similarly the raindrop distribution does not contain a significant number of rain drops + with masses almost equal or smaller than ``x^*``. These assumptions allow us + to simplify the calculation of moments of the distributions by integrating from zero to infinity. ### Autoconversion @@ -311,7 +343,7 @@ The default free parameter values are: |``N_{0,\, min}`` | ``3.5 \times 10^{5} \, m^{-4}`` | |``N_{0,\, max}`` | ``2 \times 10^{10} \, m^{-4}`` | |``\lambda_{min}`` | ``1 \times 10^{3} \, m^{-1}`` | -|``\lambda_{max}`` | ``3.5 \times 10^{4} \, m^{-1}`` | +|``\lambda_{max}`` | ``4 \times 10^{4} \, m^{-1}`` | ### Rain evaporation @@ -383,7 +415,7 @@ include("plots/RainEvapoartionSB2006.jl") ``` ![](SB2006_rain_evaporation.svg) -### Radar reflectivity +### Radar reflectivity The radar reflectivity factor (``Z``) is used to measure the power returned by a radar signal when it encounters atmospheric particles (cloud and rain droplets), and it is defined as the sixth moment of the particles distributions (``n(r)``) : ```math @@ -392,23 +424,24 @@ Z = {\int_0^\infty r^{6} \, n(r) \, dr}. \label{eq:Z} \end{equation} ``` -To take into consideration the effect of both cloud and rain droplets, we integrate separately over the two distributions defined in equations (2) and (3). +To take into consideration the effect of both cloud and rain droplets, we integrate separately over the two distributions defined in equations (2) and (4). -For cloud droplets, integrating over the assumed Gamma distribution (eq. 6) leads to +After defining ``C = \frac{4}{3} \pi \rho_w``, integrating over the assumed cloud droplets Gamma distribution (eq. 6) leads to ```math \begin{equation} -Z_c = \frac{24 \, A_c}{(B_c)^{5} \, (\frac{4}{3} \, \pi \, \rho_w)^{2}}, +Z_c = A_c C^{\nu_c+1} \frac{ (B_c C^{\mu_c})^{-\frac{3+\nu_c}{\mu_c}} \, \Gamma \left(\frac{3+\nu_c}{\mu_c}\right)}{\mu_c}, \end{equation} ``` -where ``\rho_w`` is the liquid water density. +where ``\Gamma \,(x) = \int_{0}^{\infty} \! t^{x - 1} e^{-t} \mathrm{d}t`` is the gamma function. + By performing an analogous integration for the rain droplets exponential distribution, we obtain ```math \begin{equation} -Z_r = 3 \, \frac{6! \, A_r}{(B_r)^{7} \, (\frac{4}{3} \, \pi \, \rho_w)^{2}}. +Z_r = A_r C^{\nu_r+1} \frac{ (B_r C^{\mu_r})^{-\frac{3+\nu_r}{\mu_r}} \, \Gamma \left(\frac{3+\nu_r}{\mu_r}\right)}{\mu_r}, \end{equation} ``` -To obtain the logarithmic radar reflectivity ``L_z``, which is commonly used to refer to the radar reflectivity values, we sum ``Z_c`` and ``Z_r``, and we normalize the result with the equivalent return of a ``1 mm`` drop in a volume of a meter cube (``Z_0``). +To obtain the logarithmic radar reflectivity ``L_z``, which is commonly used to refer to the radar reflectivity values, we sum ``Z_c`` and ``Z_r``, and we normalize the result with the equivalent return of a ``1 mm`` drop in a volume of a meter cube (``Z_0``). Then, we apply the decimal logarithm, and multiply the result by ``10``: ```math \begin{equation} @@ -427,23 +460,22 @@ r_{eff} = \frac{(M_{3}^c + M_{3}^r)}{M_{2}^c + M_{2}^r} = \frac{{\int_0^\infty r \label{eq:reff} \end{equation} ``` -By computing each integral, and calling ``C = \frac{4}{3} \pi \rho``, the numerator becomes +By computing each integral, the numerator becomes ```math \begin{equation} -M_{3}^c + M_{3}^r = \frac{6}{\frac{4}{3} \, \pi \, \rho_w} \, \left(\frac{A_c}{B_c^4} + \frac{3 \, A_r}{B_r^4}\right). +M_{3}^c + M_{3}^r = A_c C^{\nu_c+1} \frac{ (B_c C^{\mu_c})^{-\frac{2+\nu_c}{\mu_c}} \, \Gamma \left(\frac{2+\nu_c}{\mu_c}\right)}{\mu_c} + A_r C^{\nu_r+1} \frac{ (B_r C^{\mu_r})^{-\frac{2+\nu_r}{\mu_r}} \, \Gamma \left(\frac{2+\nu_r}{\mu_r}\right)}{\mu_r}. \end{equation} ``` -Analogously, at the denominator we have +Analogously, at the denominator we have ```math \begin{equation} -M_{2}^c + M_{2}^r = \frac{6}{(\frac{4}{3} \, \pi \, \rho_w)^{\frac{2}{3}}} \, \left(\frac{A_c \, \Gamma(\frac{11}{3})}{(B_c)^{\frac{11}{3}}} + \frac{6 \, A_r}{B_r^3}\right), +M_{2}^c + M_{2}^r = A_c C^{\nu_c+1} \frac{ (B_c C^{\mu_c})^{-\frac{5+3\nu_c}{3\mu_c}} \, \Gamma \left(\frac{5+3\nu_c}{3\mu_c}\right)}{\mu_c} + A_r C^{\nu_r+1} \frac{ (B_r C^{\mu_r})^{-\frac{5+3\nu_r}{3\mu_r}} \, \Gamma \left(\frac{5+3\nu_r}{3\mu_r}\right)}{\mu_r}. \end{equation} ``` -where ``\Gamma \,(x) = \int_{0}^{\infty} \! t^{x - 1} e^{-t} \mathrm{d}t`` is the gamma function. ### Effective radius - Liu and Hallett, 1997 An additional option for the parametrization of the effective radius is obtained following [Liu1997](@cite). -In this case +In this case ```math \begin{equation} r_{eff} = \frac{r_{vol}}{k^{\frac{1}{3}}}, @@ -451,7 +483,7 @@ In this case ``` where: - ``r_{vol}`` represents the volume-averaged radius, - - ``k = 0.8``. + - ``k = 0.8``. Where the volume-averaged radius is computed using ```math @@ -461,7 +493,7 @@ Where the volume-averaged radius is computed using ``` where: - ``L = \rho q``, is the liquid water content, - - ``N = N_{liq} + N_{rai}``. + - ``N = N_{liq} + N_{rai}``. ## Additional 2-moment microphysics options diff --git a/docs/src/plots/Microphysics2M_plots.jl b/docs/src/plots/Microphysics2M_plots.jl index 6fe4dcd5c..1cace8aa7 100644 --- a/docs/src/plots/Microphysics2M_plots.jl +++ b/docs/src/plots/Microphysics2M_plots.jl @@ -54,8 +54,14 @@ q_liq_VarTimeScaleAcnv = [ CM2.conv_q_liq_to_q_rai(VarTSc, q_liq, ρ_air, N_d) for q_liq in q_liq_range ] q_liq_SB2006 = [ - CM2.autoconversion(SB2006.acnv, q_liq, q_rai, ρ_air, N_d).dq_rai_dt for - q_liq in q_liq_range + CM2.autoconversion( + SB2006.acnv, + SB2006.pdf_c, + q_liq, + q_rai, + ρ_air, + N_d, + ).dq_rai_dt for q_liq in q_liq_range ] q_liq_K1969 = [CM1.conv_q_liq_to_q_rai(rain.acnv1M, q_liq) for q_liq in q_liq_range] @@ -71,8 +77,14 @@ N_d_LD2004 = N_d_VarTimeScaleAcnv = [CM2.conv_q_liq_to_q_rai(VarTSc, q_liq, ρ_air, N_d) for N_d in N_d_range] N_d_SB2006 = [ - CM2.autoconversion(SB2006.acnv, q_liq, q_rai, ρ_air, N_d).dq_rai_dt for - N_d in N_d_range + CM2.autoconversion( + SB2006.acnv, + SB2006.pdf_c, + q_liq, + q_rai, + ρ_air, + N_d, + ).dq_rai_dt for N_d in N_d_range ] accKK2000_q_liq = diff --git a/docs/src/plots/RainEvapoartionSB2006.jl b/docs/src/plots/RainEvapoartionSB2006.jl index 3ff2e49d1..e517ea9a4 100644 --- a/docs/src/plots/RainEvapoartionSB2006.jl +++ b/docs/src/plots/RainEvapoartionSB2006.jl @@ -30,7 +30,7 @@ function rain_evaporation_CPU(SB2006, aps, tps, q, q_rai, ρ, N_rai, T) x_star = SB2006.pdf_r.xr_min G = CO.G_func(aps, tps, T, TD.Liquid()) - xr = CM2.raindrops_limited_vars(SB2006.pdf_r, q_rai, ρ, N_rai).xr + xr = CM2.pdf_rain(SB2006.pdf_r, q_rai, ρ, N_rai).xr Dr = (FT(6) / FT(π) / ρw)^FT(1 / 3) * xr^FT(1 / 3) t_star = (FT(6) * x_star / xr)^FT(1 / 3) diff --git a/src/Microphysics2M.jl b/src/Microphysics2M.jl index ec59ef43d..e8f7a6ec4 100644 --- a/src/Microphysics2M.jl +++ b/src/Microphysics2M.jl @@ -46,19 +46,92 @@ end # mean terminal velocity of raindrops, and rain evaporation rates from Seifert and Beheng 2001 """ - raindrops_limited_vars(pdf_r, q_rai, ρ, N_rai) + pdf_cloud(pdf_c, qₗ, ρₐ, Nₗ) + + - `pdf_c` - a struct with SB2006 cloud droplets size distribution parameters + - `qₗ` - cloud water specific humidity + - `ρₐ` - air density + - `Nₗ` cloud droplet number concentration + + Returns the mean mass of cloud droplets xc [μg], + a multiplier needed to convert xc to base SI units χ [-], + and the two cloud droplet distribution parameters + Ac [m^-3, μg^-3] and Bc [μg^-1]. +""" +function pdf_cloud( + pdf_c::CMP.CloudParticlePDF_SB2006{FT}, + qₗ::FT, + ρₐ::FT, + Nₗ::FT, +) where {FT} + + (; νc, μc, ρw) = pdf_c + ϵ_N = FT(10) # TODO - should it be a parameter + + χ = 9 # convert mean droplet mass from kg to μg + + if qₗ < eps(FT) || Nₗ < ϵ_N + return ( + xc = FT(0), + Ac = FT(0), + Bc = FT(0), + Cc = FT(0), + Ec = FT(0), + ϕc = FT(0), + ψc = FT(0), + ) + else + xc = qₗ * ρₐ / Nₗ * FT(10)^χ # μg + Bc = (xc * SF.gamma((νc + 1) / μc) / SF.gamma((νc + 2) / μc))^(-μc) # 1/μg + Ac = μc * Nₗ * Bc^((νc + 1) / μc) / SF.gamma((νc + 1) / μc) # 1/m3 1/μg3 + + _Ac = Ac * FT(1e-9) # 1/mm3 1/μg3 + Cc = _Ac * ((FT(π) * ρw * FT(10)^(χ - 9)) / FT(2))^(νc + 1) / FT(3)^νc # (1/mm3)^4 + Ec = Bc * (FT(π / 6) * ρw * FT(10)^(χ - 9))^μc # 1/mm3 + ϕc = 3 * νc + 2 + ψc = 3 * μc + + return (; xc, χ, Ac, Bc, Cc, Ec, ϕc, ψc) + end +end + +""" + pdf_rain(pdf_r, q_rai, ρ, N_rai) - `pdf_r` - a struct with SB2006 raindrops size distribution parameters - `q_rai` - rain water specific humidity - `ρ` - air density - - `N_rai` raindrops number density + - `N_rai` rain drop number concentration - Returns the mean mass of raindrops, xr, the rate parameter of the assumed size distribution - of raindrops (based on drops diameter), λr, limited within prescribed ranges, and the two - rain distribution parameters, A and B. + Returns the mean mass of rain drops xr [kg], + and the two rain drop mass distribution parameters A and B. + Also returns the rate parameter of the rain drop diameter distribution λr + (optionally limited within prescribed ranges). """ -function raindrops_limited_vars( +function pdf_rain( pdf_r::CMP.RainParticlePDF_SB2006{FT}, + qᵣ::FT, + ρₐ::FT, + Nᵣ::FT, +) where {FT} + (; νr, μr, ρw) = pdf_r + ϵ_N = FT(10) # TODO - should it be a parameter? + + if qᵣ < eps(FT) || Nᵣ < ϵ_N + return (λr = FT(0), xr = FT(0), Ar = FT(0), Br = FT(0)) + else + λr = (Nᵣ * FT(π) * ρw / ρₐ / qᵣ)^FT(1 / 3) + αr = Nᵣ * λr + + xr = qᵣ * ρₐ / Nᵣ + Br = (xr * SF.gamma((νr + 1) / μr) / SF.gamma((νr + 2) / μr))^(-μr) + Ar = μr * Nᵣ * Br^((νr + 1) / μr) / SF.gamma(FT(νr + 1) / μr) + + return (; λr, αr, xr, Ar, Br) + end +end +function pdf_rain( + pdf_r::CMP.RainParticlePDF_SB2006_limited{FT}, q_rai::FT, ρ::FT, N_rai::FT, @@ -71,19 +144,20 @@ function raindrops_limited_vars( N0 = max(N0_min, min(N0_max, N_rai * (FT(π) * ρw / xr_hat)^FT(1 / 3))) λr = max(λ_min, min(λ_max, (FT(π) * ρw * N0 / L_rai)^FT(1 / 4))) xr = max(xr_min, min(xr_max, L_rai * λr / N0)) + αr = N_rai * λr Br = (xr == 0) ? FT(0) : (SF.gamma(FT(νr + 1) / μr) * xr / SF.gamma(FT(νr + 2) / μr))^(-μr) Ar = μr * N_rai * Br^(FT(νr + 1) / μr) / SF.gamma(FT(νr + 1) / μr) - return (; λr, xr, Ar, Br) + return (; λr, αr, xr, Ar, Br) end """ autoconversion(scheme, q_liq, q_rai, ρ, N_liq) - - `scheme` - type for 2-moment rain autoconversion parameterization + - `acnv`, `pdf_c` - structs with autoconversion and cloud size distribution parameters - `q_liq` - cloud water specific humidity - `q_rai` - rain water specific humidity - `ρ` - air density @@ -94,6 +168,7 @@ collisions between cloud droplets (autoconversion) for `scheme == SB2006Type` """ function autoconversion( acnv::CMP.AcnvSB2006{FT}, + pdf_c::CMP.CloudParticlePDF_SB2006{FT}, q_liq, q_rai, ρ, @@ -104,7 +179,8 @@ function autoconversion( return LiqRaiRates{FT}() end - (; kcc, νc, x_star, ρ0, A, a, b) = acnv + (; kcc, x_star, ρ0, A, a, b) = acnv + (; νc) = pdf_c L_liq = ρ * q_liq x_liq = min(x_star, L_liq / N_liq) @@ -181,6 +257,7 @@ that produce larger cloud droplets (self-collection) for `scheme == SB2006Type` """ function liquid_self_collection( acnv::CMP.AcnvSB2006{FT}, + pdf_c::CMP.CloudParticlePDF_SB2006{FT}, q_liq::FT, ρ::FT, dN_liq_dt_au::FT, @@ -189,7 +266,8 @@ function liquid_self_collection( if q_liq < eps(FT) return FT(0) end - (; kcc, ρ0, νc) = acnv + (; kcc, ρ0) = acnv + (; νc) = pdf_c L_liq = ρ * q_liq @@ -212,15 +290,15 @@ Returns a named tupple containing a LiqRaiRates object for the autoconversion ra the liquid self-collection rate for `scheme == SB2006Type` """ function autoconversion_and_liquid_self_collection( - (; acnv)::CMP.SB2006{FT}, + (; acnv, pdf_c)::CMP.SB2006{FT}, q_liq::FT, q_rai::FT, ρ::FT, N_liq::FT, ) where {FT} - au = autoconversion(acnv, q_liq, q_rai, ρ, N_liq) - sc = liquid_self_collection(acnv, q_liq, ρ, au.dN_liq_dt) + au = autoconversion(acnv, pdf_c, q_liq, q_rai, ρ, N_liq) + sc = liquid_self_collection(acnv, pdf_c, q_liq, ρ, au.dN_liq_dt) return (; au, sc) end @@ -237,7 +315,10 @@ Returns the raindrops number density tendency due to collisions of raindrops that produce larger raindrops (self-collection) for `scheme == SB2006Type` """ function rain_self_collection( - pdf::CMP.RainParticlePDF_SB2006{FT}, + pdf::Union{ + CMP.RainParticlePDF_SB2006{FT}, + CMP.RainParticlePDF_SB2006_limited{FT}, + }, self::CMP.SelfColSB2006{FT}, q_rai::FT, ρ::FT, @@ -252,9 +333,10 @@ function rain_self_collection( (; ρ0, ρw) = pdf L_rai = ρ * q_rai - λr = - raindrops_limited_vars(pdf, q_rai, ρ, N_rai).λr * - (SF.gamma(FT(4)) / FT(π) / ρw)^FT(1 / 3) + #λr = + # pdf_rain(pdf, q_rai, ρ, N_rai).λr * + # (SF.gamma(FT(4)) / FT(π) / ρw)^FT(1 / 3) + λr = pdf_rain(pdf, q_rai, ρ, N_rai).Br dN_rai_dt_sc = -krr * N_rai * L_rai * sqrt(ρ0 / ρ) * (1 + κrr / λr)^d @@ -274,7 +356,10 @@ Returns the raindrops number density tendency due to breakup of raindrops that produce smaller raindrops for `scheme == SB2006Type` """ function rain_breakup( - pdf::CMP.RainParticlePDF_SB2006{FT}, + pdf::Union{ + CMP.RainParticlePDF_SB2006{FT}, + CMP.RainParticlePDF_SB2006_limited{FT}, + }, brek::CMP.BreakupSB2006{FT}, q_rai::FT, ρ::FT, @@ -287,7 +372,7 @@ function rain_breakup( end (; Deq, Dr_th, kbr, κbr) = brek ρw = pdf.ρw - xr = raindrops_limited_vars(pdf, q_rai, ρ, N_rai).xr + xr = pdf_rain(pdf, q_rai, ρ, N_rai).xr Dr = (xr * 6 / FT(π) / ρw)^FT(1 / 3) ΔD = Dr - Deq phi_br = @@ -300,7 +385,7 @@ end """ rain_self_collection_and_breakup(SB2006, q_rai, ρ, N_rai) - - `SB2006` - a struct with SB2006 parameters for raindrops size + - `SB2006` - a struct with SB2006 parameters for raindrops size distribution, self collection, and breakup - `q_rai` - rain water specific humidity - `ρ` - air density @@ -348,7 +433,7 @@ function rain_terminal_velocity( if q_rai < eps(FT) return (FT(0), FT(0)) end - λr = raindrops_limited_vars(pdf_r, q_rai, ρ, N_rai).λr + λr = pdf_rain(pdf_r, q_rai, ρ, N_rai).λr vt0 = max(FT(0), sqrt(ρ0 / ρ) * (aR - bR / (1 + cR / λr))) vt1 = max(FT(0), sqrt(ρ0 / ρ) * (aR - bR / (1 + cR / λr)^FT(4))) return (vt0, vt1) @@ -366,7 +451,7 @@ function rain_terminal_velocity( # coefficients from Table B1 from Chen et. al. 2022 aiu, bi, ciu = CO.Chen2022_vel_coeffs_small(vel, ρ) # size distribution parameter - λ = raindrops_limited_vars(pdf_r, q_rai, ρ, N_rai).λr + λ = pdf_rain(pdf_r, q_rai, ρ, N_rai).λr # eq 20 from Chen et al 2022 vt0 = sum(CO.Chen2022_vel_add.(aiu, bi, ciu, λ, 0)) @@ -429,7 +514,7 @@ function rain_evaporation( ρw = pdf_r.ρw G = CO.G_func(aps, tps, T, TD.Liquid()) - xr = raindrops_limited_vars(pdf_r, q_rai, ρ, N_rai).xr + xr = pdf_rain(pdf_r, q_rai, ρ, N_rai).xr Dr = (FT(6) / FT(π) / ρw)^FT(1 / 3) * xr^FT(1 / 3) t_star = (FT(6) * x_star / xr)^FT(1 / 3) @@ -453,18 +538,16 @@ function rain_evaporation( return (; evap_rate_0, evap_rate_1) end -""" +""" radar_reflectivity(structs, q_liq, q_rai, N_liq, N_rai, ρ_air, τ_q, τ_N) - - `structs` - structs with SB2006 cloud droplets and raindrops + - `structs` - structs with SB2006 cloud droplets and raindrops size distributions parameters - - `q_liq` - cloud water specific humidity + - `q_liq` - cloud water specific humidity - `q_rai` - rain water specific humidity - `N_liq` - cloud droplet number density - `N_rai` - rain droplet number density - `ρ_air` - air density - - `τ_q` - threshold for minimum specific humidity value - - `τ_N` - threshold for minimum number density value Returns logarithmic radar reflectivity from the assumed cloud and rain particle size distribuions normalized by the reflectivty of 1 millimiter drop in a volume of @@ -477,50 +560,48 @@ function radar_reflectivity( N_liq::FT, N_rai::FT, ρ_air::FT, - τ_q::FT, - τ_N::FT, ) where {FT} (; νc, μc) = pdf_c - (; ρw) = pdf_r + (; νr, μr, ρw) = pdf_r C = FT(4 / 3 * π * ρw) - Z₀ = FT(1e-18) - - # change of units for N_liq - # from m^-3 to mm^-3 for accuracy - N_liq *= FT(1e-9) - q_liq = (q_liq < τ_q) ? FT(0) : q_liq - - # computation rain parameters and change units - # to have the same units that we have for clouds - Br = raindrops_limited_vars(pdf_r, q_rai, ρ_air, N_rai).Br * FT(1e-3) - Ar = raindrops_limited_vars(pdf_r, q_rai, ρ_air, N_rai).Ar * FT(1e-12) - - xc = (N_liq < τ_N) ? FT(0) : ((q_liq * ρ_air) / N_liq) - Bc = - (xc == 0) ? FT(0) : - (SF.gamma(FT(νc + 1) / μc) * xc / SF.gamma(FT(νc + 2) / μc))^(-μc) - Ac = μc * N_liq * Bc^(FT(νc + 1) / μc) / SF.gamma(FT(νc + 1) / μc) - - Zc = (Bc == 0) ? FT(0) : (FT(24) * Ac / (Bc^FT(5) * C^FT(2))) - Zr = (Br == 0) ? FT(0) : (FT(2160) * Ar / (Br^FT(7) * C^FT(2))) - - return ((Zc + Zr) == 0) ? FT(-135) : - FT(10) * (log10((Zc + Zr) / Z₀) + log10(FT(1e-9))) + log_10_Z₀ = -18 + + # rain size distribution parameters + Ar = pdf_rain(pdf_r, q_rai, ρ_air, N_rai).Ar + Br = pdf_rain(pdf_r, q_rai, ρ_air, N_rai).Br + + # cloud size distribution parameters (χ converts from μg to base SI) + Ac = pdf_cloud(pdf_c, q_liq, ρ_air, N_liq).Ac # 1/m3 1/μg3 + Bc = pdf_cloud(pdf_c, q_liq, ρ_air, N_liq).Bc # 1/μg + χ = pdf_cloud(pdf_c, q_liq, ρ_air, N_liq).χ + + Zc = + Bc == 0 ? FT(0) : + Ac * + C^(νc + 1) * + (Bc * C^μc)^(-(3 + νc) / μc) * + SF.gamma((3 + νc) / μc) / μc * FT(10)^(-2 * χ) + Zr = + Br == 0 ? FT(0) : + Ar * + C^(νr + 1) * + (Br * C^μr)^(-(3 + νr) / μr) * + SF.gamma((3 + νr) / μr) / μr + + return Zc + Zr == 0 ? FT(-135) : 10 * (log10(Zc + Zr) - log_10_Z₀) end -""" - effective_radius(structs, q_liq, q_rai, N_liq, N_rai, ρ_air, τ_q, τ_N) +""" + effective_radius(structs, q_liq, q_rai, N_liq, N_rai, ρ_air) - - `structs` - structs with SB2006 cloud droplets and raindrops + - `structs` - structs with SB2006 cloud droplets and raindrops size distribution parameters - - `q_liq` - cloud water specific humidity + - `q_liq` - cloud water specific humidity - `q_rai` - rain water specific humidity - `N_liq` - cloud droplet number density - `N_rai` - rain droplet number density - `ρ_air` - air density - - `τ_q` - threshold for minimum specific humidity value - - `τ_N` - threshold for minimum number density value Returns effective radius using the 2-moment scheme cloud and rain particle size distributions @@ -532,40 +613,48 @@ function effective_radius( N_liq::FT, N_rai::FT, ρ_air::FT, - τ_q::FT, - τ_N::FT, ) where {FT} (; νc, μc) = pdf_c - (; ρw) = pdf_r + (; νr, μr, ρw) = pdf_r C = FT(4 / 3 * π * ρw) - # change of units for N_liq - # from m^-3 to mm^-3 for accuracy - N_liq *= FT(1e-9) - q_liq = (q_liq < τ_q) ? FT(0) : q_liq - - # computation rain parameters and change units - # to have the same units that we have for clouds - Br = raindrops_limited_vars(pdf_r, q_rai, ρ_air, N_rai).Br * FT(1e-3) - Ar = raindrops_limited_vars(pdf_r, q_rai, ρ_air, N_rai).Ar * FT(1e-12) - - xc = (N_liq < τ_N) ? FT(0) : ((q_liq * ρ_air) / N_liq) - Bc = - (xc == 0) ? FT(0) : - (SF.gamma(FT(νc + 1) / μc) * xc / SF.gamma(FT(νc + 2) / μc))^(-μc) - Ac = μc * N_liq * Bc^(FT(νc + 1) / μc) / SF.gamma(FT(νc + 1) / μc) - - cloud_3moment = (Bc == 0) ? FT(0) : (FT(6) * Ac * C^3 / (Bc * C)^FT(4)) - cloud_2moment = - (Bc == 0) ? FT(0) : - (Ac * SF.gamma(FT(11 / 3)) / (Bc^FT(11 / 3) * C^FT(2 / 3))) - rain_3moment = (Br == 0) ? FT(0) : (FT(24) * Ar / (Br^FT(4) * C)) - rain_2moment = (Br == 0) ? FT(0) : (FT(6) * Ar / (Br^3 * C^FT(2 / 3))) - - return ((cloud_2moment + rain_2moment) == 0) ? FT(0) : - ((cloud_3moment + rain_3moment) / (cloud_2moment + rain_2moment)) * - FT(1e-3) + # rain size distribution parameters + Ar = pdf_rain(pdf_r, q_rai, ρ_air, N_rai).Ar + Br = pdf_rain(pdf_r, q_rai, ρ_air, N_rai).Br + + # cloud size distribution parameters (χ converts from μg to base SI) + Ac = pdf_cloud(pdf_c, q_liq, ρ_air, N_liq).Ac # 1/m3 1/μg3 + Bc = pdf_cloud(pdf_c, q_liq, ρ_air, N_liq).Bc # 1/μg + χ = pdf_cloud(pdf_c, q_liq, ρ_air, N_liq).χ + + M3_c = + Bc == 0 ? FT(0) : + Ac * + C^(νc + 1) * + (Bc * C^μc)^(-(2 + νc) / μc) * + SF.gamma((2 + νc) / μc) / μc / FT(10)^χ + M3_r = + Br == 0 ? FT(0) : + Ar * + C^(νr + 1) * + (Br * C^μr)^(-(2 + νr) / μr) * + SF.gamma((2 + νr) / μr) / μr + + M2_c = + Bc == 0 ? FT(0) : + Ac * + C^(νc + 1) * + (Bc * C^μc)^(-(5 + 3 * νc) / (3 * μc)) * + SF.gamma((5 + 3 * νc) / (3 * μc)) / μc / FT(10)^(χ * 2 / 3) + M2_r = + Br == 0 ? FT(0) : + Ar * + C^(νr + 1) * + (Br * C^μr)^(-(5 + 3 * νr) / (3 * μr)) * + SF.gamma((5 + 3 * νr) / (3 * μr)) / μr + + return M2_c + M2_r == 0 ? FT(0) : (M3_c + M3_r) / (M2_c + M2_r) end """ diff --git a/src/PrecipitationSusceptibility.jl b/src/PrecipitationSusceptibility.jl index 7056c80f5..fea1cebea 100644 --- a/src/PrecipitationSusceptibility.jl +++ b/src/PrecipitationSusceptibility.jl @@ -41,7 +41,9 @@ function precipitation_susceptibility_autoconversion( N_liq::FT, ) where {FT} grad = ForwardDiff.gradient( - x -> log(CM2.autoconversion(scheme.acnv, exp.(x)...).dq_rai_dt), + x -> log( + CM2.autoconversion(scheme.acnv, scheme.pdf_c, exp.(x)...).dq_rai_dt, + ), log.(abs.([q_liq, q_rai, ρ, N_liq])), ) return precip_susceptibility_rates( diff --git a/src/parameters/Microphysics2M.jl b/src/parameters/Microphysics2M.jl index 2caa6a1e0..61d007836 100644 --- a/src/parameters/Microphysics2M.jl +++ b/src/parameters/Microphysics2M.jl @@ -336,14 +336,15 @@ function VarTimescaleAcnv(td::CP.AbstractTOMLDict) end """ - RainParticlePDF_SB2006 + RainParticlePDF_SB2006_limited -Rain size distribution parameters from SB2006 +Rain size distribution parameters from SB2006 including the limiters +on drop maximum mass and the size distribution coefficinets N0 and lambda # Fields $(DocStringExtensions.FIELDS) """ -Base.@kwdef struct RainParticlePDF_SB2006{FT} <: ParametersType{FT} +Base.@kwdef struct RainParticlePDF_SB2006_limited{FT} <: ParametersType{FT} "Raindrop size distribution coefficient νr" νr::FT "Raindrop size distribution coefficient μr" @@ -366,7 +367,7 @@ Base.@kwdef struct RainParticlePDF_SB2006{FT} <: ParametersType{FT} ρ0::FT end -function RainParticlePDF_SB2006(td::CP.AbstractTOMLDict) +function RainParticlePDF_SB2006_limited(td::CP.AbstractTOMLDict) name_map = (; :SB2006_rain_distribution_coeff_nu => :νr, :SB2006_rain_distribution_coeff_mu => :μr, @@ -381,6 +382,40 @@ function RainParticlePDF_SB2006(td::CP.AbstractTOMLDict) ) parameters = CP.get_parameter_values(td, name_map, "CloudMicrophysics") FT = CP.float_type(td) + return RainParticlePDF_SB2006_limited{FT}(; parameters...) +end + +""" + RainParticlePDF_SB2006 + +Rain size distribution parameters from SB2006 but without the limiters + +# Fields +$(DocStringExtensions.FIELDS) +""" +Base.@kwdef struct RainParticlePDF_SB2006{FT} <: ParametersType{FT} + "Raindrop size distribution coefficient νr" + νr::FT + "Raindrop size distribution coefficient μr" + μr::FT + "Cloud liquid water density [kg/m3]" + ρw::FT + "Reference air density [kg/m3]" + ρ0::FT + "Raindrop minimal mass" + xr_min::FT +end + +function RainParticlePDF_SB2006(td::CP.AbstractTOMLDict) + name_map = (; + :SB2006_rain_distribution_coeff_nu => :νr, + :SB2006_rain_distribution_coeff_mu => :μr, + :density_liquid_water => :ρw, + :SB2006_reference_air_density => :ρ0, + :SB2006_raindrops_min_mass => :xr_min, + ) + parameters = CP.get_parameter_values(td, name_map, "CloudMicrophysics") + FT = CP.float_type(td) return RainParticlePDF_SB2006{FT}(; parameters...) end @@ -397,12 +432,15 @@ Base.@kwdef struct CloudParticlePDF_SB2006{FT} <: ParametersType{FT} νc::FT "Cloud droplet size distribution coefficient μc" μc::FT + "Cloud liquid water density [kg/m3]" + ρw::FT end function CloudParticlePDF_SB2006(td::CP.AbstractTOMLDict) name_map = (; :SB2006_cloud_gamma_distribution_parameter => :νc, :SB2006_cloud_gamma_distribution_coeff_mu => :μc, + :density_liquid_water => :ρw, ) parameters = CP.get_parameter_values(td, name_map, "CloudMicrophysics") FT = CP.float_type(td) @@ -420,8 +458,6 @@ $(DocStringExtensions.FIELDS) Base.@kwdef struct AcnvSB2006{FT} <: ParametersType{FT} "Collection kernel coefficient" kcc::FT - "Cloud gamma distribution coefficient" - νc::FT "Minimum mass of rain droplets" x_star::FT "Reference air density [kg/m3]" @@ -437,7 +473,6 @@ end function AcnvSB2006(td::CP.AbstractTOMLDict) name_map = (; :SB2006_collection_kernel_coeff_kcc => :kcc, - :SB2006_cloud_gamma_distribution_parameter => :νc, :SB2006_raindrops_min_mass => :x_star, :SB2006_reference_air_density => :ρ0, :SB2006_autoconversion_correcting_function_coeff_A => :A, @@ -578,7 +613,8 @@ end SB2006 The type and parameters for 2-moment precipitation formation by -Seifert and Beheng (2006) +Seifert and Beheng (2006). The pdf_r type choses between running with or without +limiters on raindrop size distribution parameters DOI: 10.1007/s00703-005-0112-4 @@ -602,11 +638,14 @@ struct SB2006{FT, PDc, PDr, AV, AR, SC, BR, EV} <: Precipitation2MType{FT} evap::EV end -SB2006(::Type{FT}) where {FT <: AbstractFloat} = SB2006(CP.create_toml_dict(FT)) +SB2006(::Type{FT}, is_limited = true) where {FT <: AbstractFloat} = + SB2006(CP.create_toml_dict(FT), is_limited) -function SB2006(toml_dict::CP.AbstractTOMLDict) +function SB2006(toml_dict::CP.AbstractTOMLDict, is_limited = true) pdf_c = CloudParticlePDF_SB2006(toml_dict) - pdf_r = RainParticlePDF_SB2006(toml_dict) + pdf_r = + is_limited ? RainParticlePDF_SB2006_limited(toml_dict) : + RainParticlePDF_SB2006(toml_dict) acnv = AcnvSB2006(toml_dict) accr = AccrSB2006(toml_dict) self = SelfColSB2006(toml_dict) diff --git a/src/parameters/toml/SB2006_limiters.toml b/src/parameters/toml/SB2006_limiters.toml index e034fb2a9..fbeab7139 100644 --- a/src/parameters/toml/SB2006_limiters.toml +++ b/src/parameters/toml/SB2006_limiters.toml @@ -5,7 +5,7 @@ value = 6.54e-11 value = 3.5e5 [SB2006_raindrops_size_distribution_coeff_N0_max] -value = 2e10 +value = 2e11 [SB2006_raindrops_size_distribution_coeff_lambda_max] -value = 3.5e4 \ No newline at end of file +value = 4e4 \ No newline at end of file diff --git a/test/gpu_tests.jl b/test/gpu_tests.jl index 6b4d51be3..56413ffa0 100644 --- a/test/gpu_tests.jl +++ b/test/gpu_tests.jl @@ -722,6 +722,7 @@ function test_gpu(FT) # 2-moment microphysics SB2006 = CMP.SB2006(FT) + SB2006_no_limiters = CMP.SB2006(FT, false) KK2000 = CMP.KK2000(FT) B1994 = CMP.B1994(FT) TC1980 = CMP.TC1980(FT) @@ -905,49 +906,74 @@ function test_gpu(FT) @test Array(output)[2] ≈ FT(7.2e-6) @test Array(output)[3] ≈ FT(4.7e-6) - dims = (15, 1) - (; output, ndrange) = setup_output(dims, FT) + for SB in [SB2006, SB2006_no_limiters] + dims = (15, 1) + (; output, ndrange) = setup_output(dims, FT) - T = ArrayType([FT(290)]) - qt = ArrayType([FT(7e-3)]) - ql = ArrayType([FT(2e-3)]) - qr = ArrayType([FT(5e-4)]) - ρ = ArrayType([FT(1.2)]) - Nl = ArrayType([FT(1e8)]) - Nr = ArrayType([FT(1e7)]) + T = ArrayType([FT(290)]) + qt = ArrayType([FT(7e-3)]) + ql = ArrayType([FT(2e-3)]) + qr = ArrayType([FT(5e-4)]) + ρ = ArrayType([FT(1.2)]) + Nl = ArrayType([FT(1e8)]) + Nr = ArrayType([FT(1e7)]) - kernel! = test_2_moment_SB2006_kernel!(backend, work_groups) - kernel!( - aps, - tps, - SB2006, - SB2006Vel, - output, - qt, - ql, - qr, - Nl, - Nr, - ρ, - T, - ndrange = ndrange, - ) - @test isapprox(Array(output)[1], FT(-4.083606e-7), rtol = 1e-6) - @test isapprox(Array(output)[2], FT(-3769.4827), rtol = 1e-6) - @test isapprox(Array(output)[3], FT(4.083606e-7), rtol = 1e-6) - @test isapprox(Array(output)[4], FT(1884.7413), rtol = 1e-6) - @test isapprox(Array(output)[5], FT(-31040.115), rtol = 1e-6) - @test isapprox(Array(output)[6], FT(-6.358926e-6), rtol = 1e-6) - @test isapprox(Array(output)[7], FT(-317946.28), rtol = 1e-6) - @test isapprox(Array(output)[8], FT(6.358926e-6), rtol = 1e-6) - @test isapprox(Array(output)[9], FT(0.0), rtol = 1e-6) - @test isapprox(Array(output)[10], FT(-21187.494), rtol = 1e-6) - @test isapprox(Array(output)[11], FT(14154.027), rtol = 1e-6) - @test isapprox(Array(output)[12], FT(0.9868878), rtol = 1e-6) - @test isapprox(Array(output)[13], FT(4.517734), rtol = 1e-6) - @test isapprox(Array(output)[14], FT(-243259.75126), rtol = 1e-6) - @test isapprox(Array(output)[15], FT(-0.0034601581), rtol = 1e-6) + kernel! = test_2_moment_SB2006_kernel!(backend, work_groups) + kernel!( + aps, + tps, + SB, + SB2006Vel, + output, + qt, + ql, + qr, + Nl, + Nr, + ρ, + T, + ndrange = ndrange, + ) + + @test isapprox(Array(output)[1], FT(-4.083606e-7), rtol = 1e-6) + @test isapprox(Array(output)[2], FT(-3769.4827), rtol = 1e-6) + @test isapprox(Array(output)[3], FT(4.083606e-7), rtol = 1e-6) + @test isapprox(Array(output)[4], FT(1884.7413), rtol = 1e-6) + @test isapprox(Array(output)[5], FT(-31040.115), rtol = 1e-6) + @test isapprox(Array(output)[6], FT(-6.358926e-6), rtol = 1e-6) + @test isapprox(Array(output)[7], FT(-317946.28), rtol = 1e-6) + @test isapprox(Array(output)[8], FT(6.358926e-6), rtol = 1e-6) + @test isapprox(Array(output)[9], FT(0.0), rtol = 1e-6) + if SB == SB2006 + @test isapprox(Array(output)[10], FT(-21187.494), rtol = 1e-6) + @test isapprox(Array(output)[11], FT(14154.027), rtol = 1e-6) + @test isapprox(Array(output)[12], FT(0.9868878), rtol = 1e-6) + @test isapprox(Array(output)[13], FT(4.517734), rtol = 1e-6) + @test isapprox( + Array(output)[14], + FT(-243259.75126), + rtol = 1e-6, + ) + @test isapprox( + Array(output)[15], + FT(-0.0034601581), + rtol = 1e-6, + ) + end + if SB == SB2006_no_limiters + @test isapprox(Array(output)[10], FT(-40447.855), rtol = 1e-6) + @test isapprox(Array(output)[11], FT(0), rtol = 1e-6) + @test isapprox(Array(output)[12], FT(0), rtol = 1e-6) + @test isapprox(Array(output)[13], FT(0), rtol = 1e-6) + @test isapprox(Array(output)[14], FT(-52903.817), rtol = 1e-6) + @test isapprox( + Array(output)[15], + FT(-9.3601206e-5), + rtol = 1e-6, + ) + end + end end @testset "Common Kernels" begin diff --git a/test/microphysics2M_tests.jl b/test/microphysics2M_tests.jl index a04cc6158..710f8a9f7 100644 --- a/test/microphysics2M_tests.jl +++ b/test/microphysics2M_tests.jl @@ -8,6 +8,9 @@ import CloudMicrophysics.Parameters as CMP import CloudMicrophysics.Common as CMC import CloudMicrophysics.Microphysics2M as CM2 +import QuadGK as QGK +import SpecialFunctions as SF + @info "Microphysics 2M Tests" function test_microphysics2M(FT) @@ -29,6 +32,7 @@ function test_microphysics2M(FT) ) toml_dict = CP.create_toml_dict(FT; override_file) SB2006 = CMP.SB2006(toml_dict) + SB2006_no_limiters = CMP.SB2006(toml_dict, false) # Thermodynamics and air properties parameters aps = CMP.AirProperties(FT) @@ -155,8 +159,8 @@ function test_microphysics2M(FT) for Nr in N_rai for qr in q_rai #action - λ = CM2.raindrops_limited_vars(SB2006.pdf_r, qr, ρ, Nr).λr - xr = CM2.raindrops_limited_vars(SB2006.pdf_r, qr, ρ, Nr).xr + λ = CM2.pdf_rain(SB2006.pdf_r, qr, ρ, Nr).λr + xr = CM2.pdf_rain(SB2006.pdf_r, qr, ρ, Nr).xr #test TT.@test λ_min <= λ <= λ_max @@ -173,71 +177,87 @@ function test_microphysics2M(FT) N_liq = FT(1e8) q_rai = FT(1e-6) - (; kcc, x_star, νc, ρ0) = SB2006.acnv - - #action - au = CM2.autoconversion(SB2006.acnv, q_liq, q_rai, ρ, N_liq) - sc = CM2.liquid_self_collection(SB2006.acnv, q_liq, ρ, au.dN_liq_dt) - au_sc = CM2.autoconversion_and_liquid_self_collection( - SB2006, - q_liq, - q_rai, - ρ, - N_liq, - ) - - Lc = ρ * q_liq - Lr = ρ * q_rai - xc = min(x_star, Lc / N_liq) - τ = 1 - Lc / (Lc + Lr) - ϕ_au = 400 * τ^0.7 * (1 - τ^0.7)^3 - dqrdt_au = - kcc / 20 / x_star * (νc + 2) * (νc + 4) / (νc + 1)^2 * - Lc^2 * - xc^2 * - (1 + ϕ_au / (1 - τ)^2) * - (ρ0 / ρ) / ρ - dqcdt_au = -dqrdt_au - dNcdt_au = 2 / x_star * ρ * dqcdt_au - dNrdt_au = -0.5 * dNcdt_au - dNcdt_sc = -kcc * (νc + 2) / (νc + 1) * (ρ0 / ρ) * Lc^2 - au.dN_liq_dt - - #test - TT.@test au isa CM2.LiqRaiRates - TT.@test au.dq_liq_dt ≈ dqcdt_au rtol = 1e-6 - TT.@test au.dq_rai_dt ≈ dqrdt_au rtol = 1e-6 - TT.@test au.dN_liq_dt ≈ dNcdt_au rtol = 1e-6 - TT.@test au.dN_rai_dt ≈ dNrdt_au rtol = 1e-6 - TT.@test sc ≈ dNcdt_sc rtol = 1e-6 - TT.@test au_sc isa NamedTuple - TT.@test au_sc.au.dq_liq_dt ≈ dqcdt_au rtol = 1e-6 - TT.@test au_sc.au.dq_rai_dt ≈ dqrdt_au rtol = 1e-6 - TT.@test au_sc.au.dN_liq_dt ≈ dNcdt_au rtol = 1e-6 - TT.@test au_sc.au.dN_rai_dt ≈ dNrdt_au rtol = 1e-6 - TT.@test au_sc.sc ≈ dNcdt_sc rtol = 1e-6 - - #action - au = CM2.autoconversion(SB2006.acnv, FT(0), FT(0), ρ, N_liq) - sc = CM2.liquid_self_collection(SB2006.acnv, FT(0), ρ, au.dN_liq_dt) - au_sc = CM2.autoconversion_and_liquid_self_collection( - SB2006, - FT(0), - FT(0), - ρ, - N_liq, - ) + for SB in [SB2006, SB2006_no_limiters] + (; kcc, x_star, ρ0) = SB.acnv + (; νc) = SB.pdf_c - #test - TT.@test au.dq_liq_dt ≈ FT(0) atol = eps(FT) - TT.@test au.dq_rai_dt ≈ FT(0) atol = eps(FT) - TT.@test au.dN_liq_dt ≈ FT(0) atol = eps(FT) - TT.@test au.dN_rai_dt ≈ FT(0) atol = eps(FT) - TT.@test sc ≈ FT(0) atol = eps(FT) - TT.@test au_sc.au.dq_liq_dt ≈ FT(0) atol = eps(FT) - TT.@test au_sc.au.dq_rai_dt ≈ FT(0) atol = eps(FT) - TT.@test au_sc.au.dN_liq_dt ≈ FT(0) atol = eps(FT) - TT.@test au_sc.au.dN_rai_dt ≈ FT(0) atol = eps(FT) - TT.@test au_sc.sc ≈ FT(0) atol = eps(FT) + #action + au = CM2.autoconversion(SB.acnv, SB.pdf_c, q_liq, q_rai, ρ, N_liq) + sc = CM2.liquid_self_collection( + SB.acnv, + SB.pdf_c, + q_liq, + ρ, + au.dN_liq_dt, + ) + au_sc = CM2.autoconversion_and_liquid_self_collection( + SB, + q_liq, + q_rai, + ρ, + N_liq, + ) + + Lc = ρ * q_liq + Lr = ρ * q_rai + xc = min(x_star, Lc / N_liq) + τ = 1 - Lc / (Lc + Lr) + ϕ_au = 400 * τ^0.7 * (1 - τ^0.7)^3 + dqrdt_au = + kcc / 20 / x_star * (νc + 2) * (νc + 4) / (νc + 1)^2 * + Lc^2 * + xc^2 * + (1 + ϕ_au / (1 - τ)^2) * + (ρ0 / ρ) / ρ + dqcdt_au = -dqrdt_au + dNcdt_au = 2 / x_star * ρ * dqcdt_au + dNrdt_au = -0.5 * dNcdt_au + dNcdt_sc = + -kcc * (νc + 2) / (νc + 1) * (ρ0 / ρ) * Lc^2 - au.dN_liq_dt + + #test + TT.@test au isa CM2.LiqRaiRates + TT.@test au.dq_liq_dt ≈ dqcdt_au rtol = 1e-6 + TT.@test au.dq_rai_dt ≈ dqrdt_au rtol = 1e-6 + TT.@test au.dN_liq_dt ≈ dNcdt_au rtol = 1e-6 + TT.@test au.dN_rai_dt ≈ dNrdt_au rtol = 1e-6 + TT.@test sc ≈ dNcdt_sc rtol = 1e-6 + TT.@test au_sc isa NamedTuple + TT.@test au_sc.au.dq_liq_dt ≈ dqcdt_au rtol = 1e-6 + TT.@test au_sc.au.dq_rai_dt ≈ dqrdt_au rtol = 1e-6 + TT.@test au_sc.au.dN_liq_dt ≈ dNcdt_au rtol = 1e-6 + TT.@test au_sc.au.dN_rai_dt ≈ dNrdt_au rtol = 1e-6 + TT.@test au_sc.sc ≈ dNcdt_sc rtol = 1e-6 + + #action + au = CM2.autoconversion(SB.acnv, SB.pdf_c, FT(0), FT(0), ρ, N_liq) + sc = CM2.liquid_self_collection( + SB.acnv, + SB.pdf_c, + FT(0), + ρ, + au.dN_liq_dt, + ) + au_sc = CM2.autoconversion_and_liquid_self_collection( + SB, + FT(0), + FT(0), + ρ, + N_liq, + ) + + #test + TT.@test au.dq_liq_dt ≈ FT(0) atol = eps(FT) + TT.@test au.dq_rai_dt ≈ FT(0) atol = eps(FT) + TT.@test au.dN_liq_dt ≈ FT(0) atol = eps(FT) + TT.@test au.dN_rai_dt ≈ FT(0) atol = eps(FT) + TT.@test sc ≈ FT(0) atol = eps(FT) + TT.@test au_sc.au.dq_liq_dt ≈ FT(0) atol = eps(FT) + TT.@test au_sc.au.dq_rai_dt ≈ FT(0) atol = eps(FT) + TT.@test au_sc.au.dN_liq_dt ≈ FT(0) atol = eps(FT) + TT.@test au_sc.au.dN_rai_dt ≈ FT(0) atol = eps(FT) + TT.@test au_sc.sc ≈ FT(0) atol = eps(FT) + end end TT.@testset "2M_microphysics - Seifert and Beheng 2006 accretion" begin @@ -248,37 +268,39 @@ function test_microphysics2M(FT) q_rai = FT(1e-6) N_rai = FT(1e4) - (; kcr, ρ0) = SB2006.accr - - #action - ac = CM2.accretion(SB2006, q_liq, q_rai, ρ, N_liq) - - Lc = ρ * q_liq - Lr = ρ * q_rai - xc = Lc / N_liq - τ = 1 - Lc / (Lc + Lr) - ϕ_ac = (τ / (τ + 5e-5))^4 - - dqrdt_ac = kcr * Lc * Lr * ϕ_ac * sqrt(ρ0 / ρ) / ρ - dqcdt_ac = -dqrdt_ac - dNcdt_ac = 1 / xc * ρ * dqcdt_ac - dNrdt_ac = FT(0) - - #test - TT.@test ac isa CM2.LiqRaiRates - TT.@test ac.dq_liq_dt ≈ dqcdt_ac rtol = FT(1e-6) - TT.@test ac.dq_rai_dt ≈ dqrdt_ac rtol = FT(1e-6) - TT.@test ac.dN_liq_dt ≈ dNcdt_ac rtol = FT(1e-6) - TT.@test ac.dN_rai_dt ≈ dNrdt_ac rtol = FT(1e-6) - - #action - ac = CM2.accretion(SB2006, FT(0), FT(0), ρ, N_liq) - - #test - TT.@test ac.dq_liq_dt ≈ FT(0) atol = eps(FT) - TT.@test ac.dq_rai_dt ≈ FT(0) atol = eps(FT) - TT.@test ac.dN_liq_dt ≈ FT(0) atol = eps(FT) - TT.@test ac.dN_rai_dt ≈ FT(0) atol = eps(FT) + for SB in [SB2006, SB2006_no_limiters] + (; kcr, ρ0) = SB.accr + + #action + ac = CM2.accretion(SB, q_liq, q_rai, ρ, N_liq) + + Lc = ρ * q_liq + Lr = ρ * q_rai + xc = Lc / N_liq + τ = 1 - Lc / (Lc + Lr) + ϕ_ac = (τ / (τ + 5e-5))^4 + + dqrdt_ac = kcr * Lc * Lr * ϕ_ac * sqrt(ρ0 / ρ) / ρ + dqcdt_ac = -dqrdt_ac + dNcdt_ac = 1 / xc * ρ * dqcdt_ac + dNrdt_ac = FT(0) + + #test + TT.@test ac isa CM2.LiqRaiRates + TT.@test ac.dq_liq_dt ≈ dqcdt_ac rtol = FT(1e-6) + TT.@test ac.dq_rai_dt ≈ dqrdt_ac rtol = FT(1e-6) + TT.@test ac.dN_liq_dt ≈ dNcdt_ac rtol = FT(1e-6) + TT.@test ac.dN_rai_dt ≈ dNrdt_ac rtol = FT(1e-6) + + #action + ac = CM2.accretion(SB, FT(0), FT(0), ρ, N_liq) + + #test + TT.@test ac.dq_liq_dt ≈ FT(0) atol = eps(FT) + TT.@test ac.dq_rai_dt ≈ FT(0) atol = eps(FT) + TT.@test ac.dN_liq_dt ≈ FT(0) atol = eps(FT) + TT.@test ac.dN_rai_dt ≈ FT(0) atol = eps(FT) + end end TT.@testset "2M_microphysics - Seifert and Beheng 2006 rain self-collection and breakup" begin @@ -287,65 +309,68 @@ function test_microphysics2M(FT) q_rai = FT(1e-6) N_rai = FT(1e4) - (; krr, κrr) = SB2006.self - (; Deq, kbr, κbr) = SB2006.brek - ρ0 = SB2006.pdf_r.ρ0 - - #action - sc_rai = - CM2.rain_self_collection(SB2006.pdf_r, SB2006.self, q_rai, ρ, N_rai) - br_rai = - CM2.rain_breakup(SB2006.pdf_r, SB2006.brek, q_rai, ρ, N_rai, sc_rai) - sc_br_rai = - CM2.rain_self_collection_and_breakup(SB2006, q_rai, ρ, N_rai) - - λr = - CM2.raindrops_limited_vars(SB2006.pdf_r, q_rai, ρ, N_rai).λr * - FT(6 / π / 1000)^FT(1 / 3) - dNrdt_sc = -krr * N_rai * ρ * q_rai * (1 + κrr / λr)^-5 * sqrt(ρ0 / ρ) - - Dr = - ( - CM2.raindrops_limited_vars(SB2006.pdf_r, q_rai, ρ, N_rai).xr / - 1000 / FT(π) * 6 - )^FT(1 / 3) - ΔDr = Dr - Deq - ϕ_br = - Dr < 0.35e-3 ? FT(-1) : - ((Dr < 0.9e-3) ? kbr * ΔDr : 2 * (exp(κbr * ΔDr) - 1)) - - dNrdt_br = -(ϕ_br + 1) * sc_rai - - #test - TT.@test sc_rai ≈ dNrdt_sc rtol = 1e-6 - TT.@test CM2.rain_self_collection( - SB2006.pdf_r, - SB2006.self, - FT(0), - ρ, - N_rai, - ) ≈ FT(0) atol = eps(FT) - TT.@test br_rai ≈ dNrdt_br rtol = 1e-6 - TT.@test sc_br_rai isa NamedTuple - TT.@test sc_br_rai.sc ≈ dNrdt_sc rtol = 1e-6 - TT.@test sc_br_rai.br ≈ dNrdt_br rtol = 1e-6 - - #setup - q_rai = FT(0) - - #action - sc_rai = - CM2.rain_self_collection(SB2006.pdf_r, SB2006.self, q_rai, ρ, N_rai) - br_rai = - CM2.rain_breakup(SB2006.pdf_r, SB2006.brek, q_rai, ρ, N_rai, sc_rai) - sc_br_rai = - CM2.rain_self_collection_and_breakup(SB2006, q_rai, ρ, N_rai) - - #test - TT.@test sc_rai ≈ FT(0) atol = eps(FT) - TT.@test br_rai ≈ FT(0) atol = eps(FT) - TT.@test sc_br_rai.sc ≈ FT(0) atol = eps(FT) - TT.@test sc_br_rai.br ≈ FT(0) atol = eps(FT) + for SB in [SB2006, SB2006_no_limiters] + (; krr, κrr) = SB.self + (; Deq, kbr, κbr) = SB.brek + ρ0 = SB.pdf_r.ρ0 + + #action + sc_rai = + CM2.rain_self_collection(SB.pdf_r, SB.self, q_rai, ρ, N_rai) + br_rai = + CM2.rain_breakup(SB.pdf_r, SB.brek, q_rai, ρ, N_rai, sc_rai) + sc_br_rai = + CM2.rain_self_collection_and_breakup(SB, q_rai, ρ, N_rai) + + λr = + CM2.pdf_rain(SB.pdf_r, q_rai, ρ, N_rai).λr * + FT(6 / π / 1000)^FT(1 / 3) + dNrdt_sc = + -krr * N_rai * ρ * q_rai * (1 + κrr / λr)^-5 * sqrt(ρ0 / ρ) + + Dr = + ( + CM2.pdf_rain(SB.pdf_r, q_rai, ρ, N_rai).xr / 1000 / FT(π) * + 6 + )^FT(1 / 3) + ΔDr = Dr - Deq + ϕ_br = + Dr < 0.35e-3 ? FT(-1) : + ((Dr < 0.9e-3) ? kbr * ΔDr : 2 * (exp(κbr * ΔDr) - 1)) + + dNrdt_br = -(ϕ_br + 1) * sc_rai + + #test + TT.@test sc_rai ≈ dNrdt_sc rtol = 1e-6 + TT.@test CM2.rain_self_collection( + SB.pdf_r, + SB.self, + FT(0), + ρ, + N_rai, + ) ≈ FT(0) atol = eps(FT) + TT.@test br_rai ≈ dNrdt_br rtol = 1e-6 + TT.@test sc_br_rai isa NamedTuple + TT.@test sc_br_rai.sc ≈ dNrdt_sc rtol = 1e-6 + TT.@test sc_br_rai.br ≈ dNrdt_br rtol = 1e-6 + + #setup + q_rai = FT(0) + + #action + sc_rai = + CM2.rain_self_collection(SB.pdf_r, SB.self, q_rai, ρ, N_rai) + br_rai = + CM2.rain_breakup(SB.pdf_r, SB.brek, q_rai, ρ, N_rai, sc_rai) + sc_br_rai = + CM2.rain_self_collection_and_breakup(SB, q_rai, ρ, N_rai) + + #test + TT.@test sc_rai ≈ FT(0) atol = eps(FT) + TT.@test br_rai ≈ FT(0) atol = eps(FT) + TT.@test sc_br_rai.sc ≈ FT(0) atol = eps(FT) + TT.@test sc_br_rai.br ≈ FT(0) atol = eps(FT) + end end TT.@testset "2M_microphysics - Seifert and Beheng 2006 rain terminal velocity" begin @@ -356,31 +381,34 @@ function test_microphysics2M(FT) (; ρ0, aR, bR, cR) = SB2006Vel - #action - vt_rai = CM2.rain_terminal_velocity(SB2006, SB2006Vel, q_rai, ρ, N_rai) + for SB in [SB2006, SB2006_no_limiters] - λr = CM2.raindrops_limited_vars(SB2006.pdf_r, q_rai, ρ, N_rai).λr - vt0 = max(0, sqrt(ρ0 / ρ) * (aR - bR / (1 + cR / λr))) - vt1 = max(0, sqrt(ρ0 / ρ) * (aR - bR / (1 + cR / λr)^4)) + #action + vt_rai = CM2.rain_terminal_velocity(SB, SB2006Vel, q_rai, ρ, N_rai) - #test - TT.@test vt_rai isa Tuple - TT.@test vt_rai[1] ≈ vt0 rtol = 1e-6 - TT.@test vt_rai[2] ≈ vt1 rtol = 1e-6 - TT.@test CM2.rain_terminal_velocity( - SB2006, - SB2006Vel, - FT(0), - ρ, - N_rai, - )[1] ≈ 0 atol = eps(FT) - TT.@test CM2.rain_terminal_velocity( - SB2006, - SB2006Vel, - FT(0), - ρ, - N_rai, - )[2] ≈ 0 atol = eps(FT) + λr = CM2.pdf_rain(SB.pdf_r, q_rai, ρ, N_rai).λr + vt0 = max(0, sqrt(ρ0 / ρ) * (aR - bR / (1 + cR / λr))) + vt1 = max(0, sqrt(ρ0 / ρ) * (aR - bR / (1 + cR / λr)^4)) + + #test + TT.@test vt_rai isa Tuple + TT.@test vt_rai[1] ≈ vt0 rtol = 1e-6 + TT.@test vt_rai[2] ≈ vt1 rtol = 1e-6 + TT.@test CM2.rain_terminal_velocity( + SB, + SB2006Vel, + FT(0), + ρ, + N_rai, + )[1] ≈ 0 atol = eps(FT) + TT.@test CM2.rain_terminal_velocity( + SB, + SB2006Vel, + FT(0), + ρ, + N_rai, + )[2] ≈ 0 atol = eps(FT) + end end TT.@testset "2M_microphysics - Chen 2022 rain terminal velocity" begin @@ -389,34 +417,36 @@ function test_microphysics2M(FT) q_rai = FT(5e-4) N_rai = FT(1e4) - #action - vt_rai = - CM2.rain_terminal_velocity(SB2006, Chen2022Vel, q_rai, ρ, N_rai) - v_bigger = - CM2.rain_terminal_velocity(SB2006, Chen2022Vel, q_rai * 2, ρ, N_rai) - - #test - TT.@test vt_rai isa Tuple - TT.@test vt_rai[1] ≈ 1.0738503635546666 rtol = 1e-6 - TT.@test vt_rai[2] ≈ 4.00592218028957 rtol = 1e-6 - - TT.@test CM2.rain_terminal_velocity( - SB2006, - Chen2022Vel, - FT(0), - ρ, - N_rai, - )[1] ≈ 0 atol = eps(FT) - TT.@test CM2.rain_terminal_velocity( - SB2006, - Chen2022Vel, - FT(0), - ρ, - N_rai, - )[2] ≈ 0 atol = eps(FT) + for SB in [SB2006, SB2006_no_limiters] + #action + vt_rai = + CM2.rain_terminal_velocity(SB, Chen2022Vel, q_rai, ρ, N_rai) + v_bigger = + CM2.rain_terminal_velocity(SB, Chen2022Vel, q_rai * 2, ρ, N_rai) + + #test + TT.@test vt_rai isa Tuple + TT.@test vt_rai[1] ≈ 1.0738503635546666 rtol = 1e-6 + TT.@test vt_rai[2] ≈ 4.00592218028957 rtol = 1e-6 + + TT.@test CM2.rain_terminal_velocity( + SB, + Chen2022Vel, + FT(0), + ρ, + N_rai, + )[1] ≈ 0 atol = eps(FT) + TT.@test CM2.rain_terminal_velocity( + SB, + Chen2022Vel, + FT(0), + ρ, + N_rai, + )[2] ≈ 0 atol = eps(FT) - TT.@test v_bigger[1] > vt_rai[1] - TT.@test v_bigger[2] > vt_rai[2] + TT.@test v_bigger[1] > vt_rai[1] + TT.@test v_bigger[2] > vt_rai[2] + end end TT.@testset "2M_microphysics - Seifert and Beheng 2006 rain evaporation" begin @@ -428,53 +458,56 @@ function test_microphysics2M(FT) q_tot = FT(1e-3) q = TD.PhasePartition(q_tot) - (; av, bv, α, β, ρ0) = SB2006.evap - (; ν_air, D_vapor) = aps - - #action - evap = CM2.rain_evaporation(SB2006, aps, tps, q, q_rai, ρ, N_rai, T) - - G = CMC.G_func(aps, tps, T, TD.Liquid()) - S = TD.supersaturation(tps, q, ρ, T, TD.Liquid()) - - xr = CM2.raindrops_limited_vars(SB2006.pdf_r, q_rai, ρ, N_rai).xr - Dr = FT(6 / π / 1000.0)^FT(1 / 3) * xr^FT(1 / 3) - N_Re = α * xr^β * sqrt(ρ0 / ρ) * Dr / ν_air - - a_vent_0 = av * FT(0.15344374450453543) - b_vent_0 = bv * FT(0.17380986321413017) - Fv0 = a_vent_0 + b_vent_0 * (ν_air / D_vapor)^FT(1 / 3) * sqrt(N_Re) - a_vent_1 = av * FT(0.5503212081491045) - b_vent_1 = bv * FT(0.5873135598802672) - Fv1 = a_vent_1 + b_vent_1 * (ν_air / D_vapor)^FT(1 / 3) * sqrt(N_Re) - - evap0 = 2 * FT(π) * G * S * N_rai * Dr * Fv0 / xr - evap1 = 2 * FT(π) * G * S * N_rai * Dr * Fv1 / ρ - - #test - TT.@test evap isa NamedTuple - TT.@test evap.evap_rate_0 ≈ evap0 rtol = 1e-4 - TT.@test evap.evap_rate_1 ≈ evap1 rtol = 1e-5 - TT.@test CM2.rain_evaporation( - SB2006, - aps, - tps, - q, - FT(0), - ρ, - N_rai, - T, - ).evap_rate_0 ≈ 0 atol = eps(FT) - TT.@test CM2.rain_evaporation( - SB2006, - aps, - tps, - q, - FT(0), - ρ, - N_rai, - T, - ).evap_rate_1 ≈ 0 atol = eps(FT) + for SB in [SB2006, SB2006_no_limiters] + + (; av, bv, α, β, ρ0) = SB.evap + (; ν_air, D_vapor) = aps + + #action + evap = CM2.rain_evaporation(SB, aps, tps, q, q_rai, ρ, N_rai, T) + + G = CMC.G_func(aps, tps, T, TD.Liquid()) + S = TD.supersaturation(tps, q, ρ, T, TD.Liquid()) + + xr = CM2.pdf_rain(SB.pdf_r, q_rai, ρ, N_rai).xr + Dr = FT(6 / π / 1000.0)^FT(1 / 3) * xr^FT(1 / 3) + N_Re = α * xr^β * sqrt(ρ0 / ρ) * Dr / ν_air + + a_vent_0 = av * FT(0.15344374450453543) + b_vent_0 = bv * FT(0.17380986321413017) + Fv0 = a_vent_0 + b_vent_0 * (ν_air / D_vapor)^FT(1 / 3) * sqrt(N_Re) + a_vent_1 = av * FT(0.5503212081491045) + b_vent_1 = bv * FT(0.5873135598802672) + Fv1 = a_vent_1 + b_vent_1 * (ν_air / D_vapor)^FT(1 / 3) * sqrt(N_Re) + + evap0 = 2 * FT(π) * G * S * N_rai * Dr * Fv0 / xr + evap1 = 2 * FT(π) * G * S * N_rai * Dr * Fv1 / ρ + + #test + TT.@test evap isa NamedTuple + TT.@test evap.evap_rate_0 ≈ evap0 rtol = 1e-4 + TT.@test evap.evap_rate_1 ≈ evap1 rtol = 1e-5 + TT.@test CM2.rain_evaporation( + SB, + aps, + tps, + q, + FT(0), + ρ, + N_rai, + T, + ).evap_rate_0 ≈ 0 atol = eps(FT) + TT.@test CM2.rain_evaporation( + SB, + aps, + tps, + q, + FT(0), + ρ, + N_rai, + T, + ).evap_rate_1 ≈ 0 atol = eps(FT) + end end TT.@testset "2M_microphysics - Seifert and Beheng 2006 radar reflectivity" begin @@ -484,23 +517,14 @@ function test_microphysics2M(FT) N_liq = FT(15053529) q_rai = FT(1.573e-4) N_rai = FT(510859) - τ_q = FT(1e-12) - τ_N = FT(1e-18) - - #action - rr = CM2.radar_reflectivity( - SB2006, - q_liq, - q_rai, - N_liq, - N_rai, - ρ_air, - τ_q, - τ_N, - ) - TT.@test rr ≈ FT(-13) atol = FT(0.5) + for SB in [SB2006, SB2006_no_limiters] + #action + rr = CM2.radar_reflectivity(SB, q_liq, q_rai, N_liq, N_rai, ρ_air) + # test + TT.@test rr ≈ FT(-13) atol = FT(0.5) + end end TT.@testset "2M_microphysics - Seifert and Beheng 2006 effective radius" begin @@ -511,24 +535,14 @@ function test_microphysics2M(FT) N_liq = FT(15053529) q_rai = FT(1.573e-4) N_rai = FT(510859) - τ_q = FT(1e-12) - τ_N = FT(1e-18) - - #action - reff = CM2.effective_radius( - SB2006, - q_liq, - q_rai, - N_liq, - N_rai, - ρ_air, - τ_q, - τ_N, - ) - #test - TT.@test reff ≈ FT(2.66e-05) atol = FT(3e-7) + for SB in [SB2006, SB2006_no_limiters] + #action + reff = CM2.effective_radius(SB, q_liq, q_rai, N_liq, N_rai, ρ_air) + #test + TT.@test reff ≈ FT(2.66e-05) atol = FT(8e-6) + end end TT.@testset "2M_microphysics - '1/3' power law from Liu and Hallett (1997)" begin @@ -554,6 +568,118 @@ function test_microphysics2M(FT) TT.@test reff ≈ FT(2.66e-05) atol = FT(8e-6) end + + TT.@testset "2M_microphysics - Seifert and Beheng 2006 rain distribution sanity checks" begin + + # air and liquid water densities + ρₐ = FT(1.2) + ρₗ = SB2006.pdf_r.ρw + + # example number concentration and specific humidity + Nᵣ = FT(5e5) + qᵣ = FT(5e-4) + + # limited distribution parameters for rain + Ar_l = CM2.pdf_rain(SB2006.pdf_r, qᵣ, ρₐ, Nᵣ).Ar + Br_l = CM2.pdf_rain(SB2006.pdf_r, qᵣ, ρₐ, Nᵣ).Br + λr_l = CM2.pdf_rain(SB2006.pdf_r, qᵣ, ρₐ, Nᵣ).λr + αr_l = CM2.pdf_rain(SB2006.pdf_r, qᵣ, ρₐ, Nᵣ).αr + # not limited distribution parameters for rain + (; λr, αr, Ar, Br) = CM2.pdf_rain(SB2006_no_limiters.pdf_r, qᵣ, ρₐ, Nᵣ) + + # mass of liquid droplet as a function of its diameter + m(D) = FT(π / 6) * ρₗ * D^3 + + # rain drop diameter distribution (eq.(3) from 2M docs) + f_D(D) = αr * exp(-λr * D) + # rain drop diameter distribution, but using SB2006 limiters + f_D_limited(D) = αr_l * exp(-λr_l * D) + + # rain drop mass distribution (eq.(4) from 2M docs) + f_x(x) = Ar * x^(-2 / 3) * exp(-Br * x^(1 / 3)) + # rain drop mass distribution, but using the SB2006 limiters + f_x_limited(x) = Ar_l * x^(-2 / 3) * exp(-Br_l * x^(1 / 3)) + + # integral bounds + D₀ = 1e-7 + D∞ = 1e-2 + m₀ = m(D₀) + m∞ = m(D∞) + + # Sanity checks for number concentrations for rain + ND = QGK.quadgk(x -> f_D(x), D₀, D∞)[1] + Nx = QGK.quadgk(x -> f_x(x), m₀, m∞)[1] + ND_lim = QGK.quadgk(x -> f_D_limited(x), D₀, D∞)[1] + Nx_lim = QGK.quadgk(x -> f_x_limited(x), m₀, m∞)[1] + TT.@test ND ≈ Nᵣ rtol = FT(1e-2) + TT.@test Nx ≈ Nᵣ rtol = FT(1e-2) + TT.@test ND_lim ≈ Nᵣ rtol = FT(1e-2) + TT.@test Nx_lim ≈ Nᵣ rtol = FT(1e-2) + + # Sanity checks for specific humidities for rain + qD = QGK.quadgk(x -> m(x) * f_D(x), D₀, D∞)[1] / ρₐ + qx = QGK.quadgk(x -> x * f_x(x), m₀, m∞)[1] / ρₐ + qD_lim = QGK.quadgk(x -> m(x) * f_D_limited(x), D₀, D∞)[1] / ρₐ + qx_lim = QGK.quadgk(x -> x * f_x_limited(x), m₀, m∞)[1] / ρₐ + TT.@test qD ≈ qᵣ atol = FT(1e-6) + TT.@test qx ≈ qᵣ atol = FT(1e-6) + + # The mass integrals don't work with limiters + TT.@test qx_lim ≈ qᵣ atol = FT(1e-6) + TT.@test qD_lim ≈ qx_lim atol = FT(1e-6) + end + + TT.@testset "2M_microphysics - Seifert and Beheng 2006 cloud distribution sanity checks" begin + + # example number concentration and specific humidity + Nₗ = FT(1e9) + qₗ = FT(1e-3) + + # air and liquid water densities in μg/m3 + ρₐ = FT(1.2) + ρₗ = SB2006.pdf_r.ρw + + # distribution parameters for cloud, units: [Bc] = 1/μg, [Ac] = 1/m3 1/μg3 + (; χ, Ac, Bc, Cc, Ec, ϕc, ψc) = + CM2.pdf_cloud(SB2006_no_limiters.pdf_c, qₗ, ρₐ, Nₗ) + + # mass of liquid droplet as a function of its diameter in μg + m(D) = FT(π / 6) * ρₗ * FT(10)^χ * D^3 + + # cloud droplet mass distribution (eq.(2) from 2M docs) + f_x(x) = Ac * x^(2) * exp(-Bc * x^(1)) + + # cloud droplet diameter distribution (eq.(7) from 2M docs) + f_D(D) = Cc * D^ϕc * exp(-Ec * D^ψc) + + # integral bounds + D₀ = 1e-8 + D∞ = 1e-4 + m₀ = m(D₀) + m∞ = m(D∞) + + # Sanity checks specific humidity and number concentration with mass distribution + # Sanity checks for number concentrations for cloud + qx = QGK.quadgk(x -> x * f_x(x), m₀, m∞)[1] / (ρₐ * FT(10)^χ) + Nx = QGK.quadgk(x -> f_x(x), m₀, m∞)[1] + TT.@test qx ≈ qₗ rtol = FT(1e-6) + TT.@test Nx ≈ Nₗ rtol = FT(1e-6) + + # mass of liquid droplets as a function of its diameter in mm and μg + _m(D) = FT(π / 6) * ρₗ * FT(10)^(χ - 9) * D^3 + + # integral bounds in millimiters + _D₀ = 1e-8 * FT(1e3) + _D∞ = 1e-4 * FT(1e3) + + # Sanity checks specific humidity and number concentration with diameter distribution + qD = + QGK.quadgk(x -> _m(x) * f_D(x), _D₀, _D∞)[1] / (ρₐ * FT(10)^(χ - 9)) + ND = QGK.quadgk(x -> f_D(x), _D₀, _D∞)[1] * FT(1e9) + TT.@test ND ≈ Nₗ rtol = FT(1e-6) + TT.@test qD ≈ qₗ rtol = FT(1e-6) + + end end println("Testing Float64") diff --git a/test/performance_tests.jl b/test/performance_tests.jl index b4f8c2291..1e0f9a9b2 100644 --- a/test/performance_tests.jl +++ b/test/performance_tests.jl @@ -64,6 +64,8 @@ function benchmark_test(FT) ) toml_dict = CP.create_toml_dict(FT; override_file) sb2006 = CMP.SB2006(toml_dict) + sb2006_no_limiters = CMP.SB2006(toml_dict, false) + # P3 scheme p3 = CMP.ParametersP3(FT) # terminal velocity @@ -206,44 +208,46 @@ function benchmark_test(FT) bench_press(CM1.radar_reflectivity, (rain, q_rai, ρ_air), 250) # 2-moment - bench_press( - CM2.autoconversion_and_liquid_self_collection, - (sb2006, q_liq, q_rai, ρ_air, N_liq), - 250, - ) - bench_press( - CM2.rain_self_collection_and_breakup, - (sb2006, q_rai, ρ_air, N_rai), - 1200, - ) - bench_press( - CM2.rain_evaporation, - (sb2006, aps, tps, q, q_rai, ρ_air, N_rai, T_air), - 2000, - ) - bench_press( - CM2.rain_terminal_velocity, - (sb2006, sb2006vel, q_rai, ρ_air, N_rai), - 700, - ) - bench_press( - CM2.rain_terminal_velocity, - (sb2006, ch2022.rain, q_rai, ρ_air, N_rai), - 2000, - ) - bench_press( - CM2.radar_reflectivity, - (sb2006, q_liq, q_rai, N_liq, N_rai, FT(1), FT(1e-12), FT(1e-18)), - 1200, - ) - bench_press( - CM2.effective_radius, - (sb2006, q_liq, q_rai, N_liq, N_rai, FT(1), FT(1e-12), FT(1e-18)), - 1200, - ) + for sb in [sb2006, sb2006_no_limiters] + bench_press( + CM2.autoconversion_and_liquid_self_collection, + (sb, q_liq, q_rai, ρ_air, N_liq), + 250, + ) + bench_press( + CM2.rain_self_collection_and_breakup, + (sb, q_rai, ρ_air, N_rai), + 1200, + ) + bench_press( + CM2.rain_evaporation, + (sb, aps, tps, q, q_rai, ρ_air, N_rai, T_air), + 2000, + ) + bench_press( + CM2.rain_terminal_velocity, + (sb, sb2006vel, q_rai, ρ_air, N_rai), + 700, + ) + bench_press( + CM2.rain_terminal_velocity, + (sb, ch2022.rain, q_rai, ρ_air, N_rai), + 2000, + ) + bench_press( + CM2.radar_reflectivity, + (sb, q_liq, q_rai, N_liq, N_rai, ρ_air), + 2000, + ) + bench_press( + CM2.effective_radius, + (sb, q_liq, q_rai, N_liq, N_rai, ρ_air), + 2000, + ) + end bench_press( CM2.effective_radius_Liu_Hallet_97, - (q_liq, q_rai, N_liq, N_rai, FT(1), FT(1000)), + (q_liq, q_rai, N_liq, N_rai, ρ_air, FT(1000)), 300, ) # Homogeneous Nucleation