diff --git a/docs/src/API.md b/docs/src/API.md index 5d39ca4804..ead6d061d0 100644 --- a/docs/src/API.md +++ b/docs/src/API.md @@ -96,6 +96,9 @@ Common.logistic_function Common.logistic_function_integral Common.H2SO4_soln_saturation_vapor_pressure Common.Delta_a_w +Common.Chen2022_snow_ice_coeffs +Common.Chen2022_vel_add +Common.Chen2022_vel_coeffs ``` # Common utility types @@ -115,6 +118,9 @@ CommonTypes.B1994Type CommonTypes.TC1980Type CommonTypes.LD2004Type CommonTypes.SB2006Type +CommonTypes.AbstractTerminalVelocityType +CommonTypes.Blk1MVelType +CommonTypes.SB2006VelType CommonTypes.Chen2022Type CommonTypes.AbstractAerosolType CommonTypes.ArizonaTestDustType diff --git a/docs/src/Microphysics1M.md b/docs/src/Microphysics1M.md index 13ec10a308..4de0bcdc0c 100644 --- a/docs/src/Microphysics1M.md +++ b/docs/src/Microphysics1M.md @@ -109,13 +109,99 @@ where: - ``\rho`` is the density of air. !!! note + Assuming a constant drag coefficient is an approximation and it should + be size and flow dependent, see for example + [here](https://www1.grc.nasa.gov/beginners-guide-to-aeronautics/drag-of-a-sphere/). + We are in the process of updating the 1-moment microphysics scheme to formulae from [Chen2022](@cite). + Other possibilities: [Khvorostyanov2002](@cite) or [Karrer2020](@cite) + +[Chen2022](@cite) provides a terminal velocity parameterisation + based on an empirical fit to a high accuracy model. +The terminal velocity depends on particle shape, size and density, + consideres the deformation effects of large rain drops, + as well as size-specific air density dependence. +The fall speed of individual particles $v(D)$ is parameterized as: +```math +\begin{equation} + v_{term}(D) = \phi^{\kappa} \sum_{i=1}^{j} \; a_i D^{b_i} e^{-c_i \; D} +\end{equation} +``` +where: + - ``D`` is the particle diameter, + - ``a_i``, ``b_i``, ``c_i`` are the free parameters, + - ``\phi`` is the aspect ratio, and + - ``\kappa`` is a parameter that depends on the particle shape (``\kappa=0`` for spheres, ``\kappa=-1/3`` for oblate and ``\kappa=1/6`` for prolate spheroids). + +For ice and snow ``j=2`` and for rain ``j=3``, to account for deformation at larger sizes. +For rain and ice we assume ``\phi=1`` (spherical). +For snow we assume ``\kappa = -1/3`` and + find the aspect ratio that is consistent with the assumed ``m(r)`` and ``a(r)`` relationships. +The aspect ratio is defined as: +```math +\begin{equation} + \phi \equiv c/a +\end{equation} +``` +where: + - ``a`` is the basal plane axial half-length, and + - ``c`` is perpendicular to the basal plane. + +The volume of a spheroid can be represented as ``V_p = 4\pi/3 \; a^2 c`` + and the area can be represented as ``A_p = \pi a c``. +It follows that + ``c = (4A_p^2) / (3 \pi V_p)``, + ``a = (3V_p) / (4A_p)``, and + ``\phi = (16 A_p^3) / (9 \pi V_p^2)``. +The volume and area are defined by the assumed power-law size relations + ``V_p = m(r) / (\rho_{ice})``, ``A_p = a(r)``. +As a result the terminal velocity of individual snow particles as: +```math +\begin{equation} + v_{term}(r) = \left(\frac{16 \; \rho_{ice}^2 \; a_0^3 \; (r/r_0)^{3a_e}}{9 \pi \; m_0^2 \; (r/r_0)^{2 m_e}} \right)^{\kappa} \sum_{i=1}^{2} \; a_i (2r)^{b_i} e^{-2 c_i r} +\end{equation} +``` +where $r$ is the radius of the particle. + +Here, we plot the aspect ratio as a function of $r$, the radius of the particle. +```@example +include("TerminalVelocityComparisons.jl") +``` +![](aspect_ratio_1M_snow.svg) + +The rain ``a_i``, ``b_i``, and ``c_i`` are listed in the table below. +The formula is applicable when ``D > 0.1 mm``, +$q$ refers to ``q = e^{0.115231 \; \rho_a}``, where ``\rho_a`` is air density [kg/m3]. +The units are: [v] = m/s, [D]=mm, [``a_i``] = ``mm^{-b_i} m/s``, [``b_i``] is dimensionless, [``c_i``] = 1/mm. + +| ``i`` | ``a_i`` | ``b_i`` | ``c_i`` | +|---------|-------------------------------------------|--------------------------------------|---------------------| +| 1 | `` 0.044612 \; q`` | ``2.2955 \; -0.038465 \; \rho_a`` | ``0`` | +| 2 | ``-0.263166 \; q`` | ``2.2955 \; -0.038465 \; \rho_a`` | ``0.184325`` | +| 3 | ``4.7178 \; q \; (\rho_a)^{-0.47335}`` | ``1.1451 \; -0.038465 \; \rho_a`` | ``0.184325`` | + +The ice and snow ``a_i``, ``b_i``, and ``c_i`` are listed in the table below. +The formula is applicable when ``D < 0.625 mm``. + +| ``i`` | ``a_i`` | ``b_i`` | ``c_i`` | +|-------|-----------------------------|-------------------------|-----------| +| 1 | ``E_s (\rho_a)^{A_s}`` | ``B_s + C_s \rho_a`` | ``0`` | +| 2 | ``F_s (\rho_a)^{A_s}`` | ``B_s + C_s \rho_a`` | ``G_s`` | + +| Coefficient | Formula | +|--------------|------------------------------------------------------------------------------------------------| +| ``A_s`` | ``0.00174079 \log{(\rho_{ice})}^2 − 0.0378769 \log{(\rho_{ice})} - 0.263503`` | +| ``B_s`` | ``(0.575231 + 0.0909307 \log{(\rho_{ice})} + 0.515579 / \sqrt{\rho_{ice}})^{-1}`` | +| ``C_s`` | ``-0.345387 + 0.177362 \, \exp{(-0.000427794 \rho_{ice})} + 0.00419647 \sqrt{\rho_{ice}}`` | +| ``E_s`` | ``-0.156593 - 0.0189334 \log{(\rho_{ice})}^2 + 0.1377817 \sqrt{\rho_{ice}}`` | +| ``F_s`` | ``- \exp{[-3.35641 - 0.0156199 \log{\rho_{ice}}^2 + 0.765337 \log{\rho_{ice}}]}`` | +| ``G_s`` | ``(-0.0309715 + 1.55054 / \log{(\rho_{ice})} - 0.518349 log{(\rho_{ice})} / \rho_{ice})^{-1}`` | + +The terminal velocity formulae from the current default 1-moment scheme and [Chen2022](@cite) are plotted below. +```@example +include("TerminalVelocityComparisons.jl") +``` +![](1M_individual_terminal_velocity_comparisons.svg) - It would be great to replace the above simple power laws - with more accurate relationships. - For example: - [Khvorostyanov2002](@cite) - or - [Karrer2020](@cite) ## Assumed particle size distributions @@ -259,8 +345,7 @@ F(r) = a_{vent} + ## Terminal velocity -The mass weighted terminal velocity ``v_t`` is defined following -[Ogura1971](@cite): +The mass weighted terminal velocity ``v_t`` (following [Ogura1971](@cite)) is: ```math \begin{equation} v_t = \frac{\int_0^\infty n(r) \, m(r) \, v_{term}(r) \, dr} @@ -268,8 +353,8 @@ The mass weighted terminal velocity ``v_t`` is defined following \label{eq:vt} \end{equation} ``` -Integrating over the assumed Marshall-Palmer distribution and using the - ``m(r)`` and ``v_{term}(r)`` relationships results in +Integrating the default 1-moment ``m(r)`` and ``v_{term}(r)`` relationships + over the assumed Marshall-Palmer distribution results in group terminal velocity: ```math \begin{equation} v_t = \chi_v \, v_0 \, \left(\frac{1}{r_0 \, \lambda}\right)^{v_e + \Delta_v} @@ -278,15 +363,30 @@ Integrating over the assumed Marshall-Palmer distribution and using the \end{equation} ``` -!!! note - - Assuming a constant drag coefficient is an approximation and it should - be size and flow dependent, see for example - [here](https://www1.grc.nasa.gov/beginners-guide-to-aeronautics/drag-of-a-sphere/). - In general we should implement these terminal velocity parameterizations: - [Khvorostyanov2002](@cite) - or - [Karrer2020](@cite) +Integrating [Chen2022](@cite) formulae for rain and ice + over the assumed Marshall-Palmer size distribution, + results in group terminal velocity: +```math +\begin{equation} + v_t = \sum_{i=1}^{j} \frac{a_i \lambda^{\delta} \Gamma(b_i + \delta)}{(\lambda + c_i)^{b_i + 4} \; \Gamma(4)} +\end{equation} +``` +Finally, integrating [Chen2022](@cite) formulae for snow + over the assumed Marshall-Palmer distribution, + results in group terminal velocity: +```math +\begin{equation} + v_t = \sum_{i=1}^{2} t_i \frac{\Gamma(3a_e \kappa - 2 m_e \kappa + b_i + k + 1)}{\Gamma(k+1)} +\end{equation} +``` +where: +```math +\begin{equation} +t_i = \frac{[16 a_0^3 \rho_{ice}^2]^{\kappa} \; a_i \; 2^{b_i} [2 c_i \lambda]^{-(3 a_e \kappa - 2 m_e \kappa + b_i + k)-1}} +{[9 \pi m_0^2]^{\kappa} \; r_0^{3 a_e \kappa - 2 m_e \kappa} \lambda^{-k-1}} +\end{equation} +``` +and $k = 3$. ## Rain autoconversion @@ -679,6 +779,8 @@ const liquid = CMT.LiquidType() const ice = CMT.IceType() const rain = CMT.RainType() const snow = CMT.SnowType() +const Chen2022 = CMT.Chen2022Type() +const Blk1MVel = CMT.Blk1MVelType() # eq. 5d in [Grabowski1996](@cite) function terminal_velocity_empirical(q_rai::DT, q_tot::DT, ρ::DT, ρ_air_ground::DT) where {DT<:Real} @@ -721,9 +823,12 @@ q_snow_range = range(1e-8, stop=5e-3, length=100) ρ_air, ρ_air_ground = 1.2, 1.22 q_liq, q_ice, q_tot = 5e-4, 5e-4, 20e-3 -PL.plot( q_rain_range * 1e3, [CM1.terminal_velocity(param_set, rain, ρ_air, q_rai) for q_rai in q_rain_range], linewidth=3, xlabel="q_rain or q_snow [g/kg]", ylabel="terminal velocity [m/s]", label="Rain-CLIMA") -PL.plot!(q_snow_range * 1e3, [CM1.terminal_velocity(param_set, snow, ρ_air, q_sno) for q_sno in q_snow_range], linewidth=3, label="Snow-CLIMA") -PL.plot!(q_rain_range * 1e3, [terminal_velocity_empirical(q_rai, q_tot, ρ_air, ρ_air_ground) for q_rai in q_rain_range], linewidth=3, label="Rain-Empirical") +PL.plot( q_rain_range * 1e3, [CM1.terminal_velocity(param_set, rain, Blk1MVel, ρ_air, q_rai) for q_rai in q_rain_range], linewidth=3, label="Rain-1M-default", color=:blue, xlabel="q_rain/ice/snow [g/kg]", ylabel="terminal velocity [m/s]") +PL.plot!(q_rain_range * 1e3, [CM1.terminal_velocity(param_set, rain, Chen2022, ρ_air, q_rai) for q_rai in q_rain_range], linewidth=3, label="Rain-Chen", color=:blue, linestyle=:dot) +PL.plot!(q_rain_range * 1e3, [terminal_velocity_empirical(q_rai, q_tot, ρ_air, ρ_air_ground) for q_rai in q_rain_range], linewidth=3, label="Rain-Empirical", color=:green) +PL.plot!(q_ice_range * 1e3, [CM1.terminal_velocity(param_set, ice, Chen2022, ρ_air, q_ice) for q_ice in q_ice_range], linewidth=3, label="Ice-Chen", color=:pink, linestyle=:dot) +PL.plot!(q_snow_range * 1e3, [CM1.terminal_velocity(param_set, snow, Blk1MVel, ρ_air, q_sno) for q_sno in q_snow_range], linewidth=3, label="Snow-1M-default", color=:red) +PL.plot!(q_snow_range * 1e3, [CM1.terminal_velocity(param_set, snow, Chen2022, ρ_air, q_sno) for q_sno in q_snow_range], linewidth=3, label="Snow-Chen", color=:red, linestyle=:dot) PL.savefig("terminal_velocity.svg") # hide T = 273.15 @@ -830,4 +935,3 @@ PL.savefig("snow_melt_rate.svg") # hide ![](rain_evaporation_rate.svg) ![](snow_sublimation_deposition_rate.svg) ![](snow_melt_rate.svg) - diff --git a/docs/src/Microphysics2M.md b/docs/src/Microphysics2M.md index 2dfce1531a..3a60a33bdd 100644 --- a/docs/src/Microphysics2M.md +++ b/docs/src/Microphysics2M.md @@ -403,12 +403,10 @@ Below, we compare the individual terminal velocity formulas for [Chen2022](@cite We also compare bulk number weighted [ND] and mass weighted [MD] terminal velocities for both formulas integrated over the size distribution from [SeifertBeheng2006](@cite). We also show the mass weighted terminal velocity from the 1-moment scheme. - ```@example include("TerminalVelocityComparisons.jl") ``` -![](terminal_velocity_bulk_comparisons.svg) -![](terminal_velocity_individual_raindrop.svg) +![](2M_terminal_velocity_comparisons.svg) ### Accretion and Autoconversion diff --git a/docs/src/TerminalVelocityComparisons.jl b/docs/src/TerminalVelocityComparisons.jl index d8c9a15af3..d653c9d264 100644 --- a/docs/src/TerminalVelocityComparisons.jl +++ b/docs/src/TerminalVelocityComparisons.jl @@ -1,37 +1,41 @@ import Plots as PL -FT = Float64 - import CloudMicrophysics as CM import CLIMAParameters as CP +FT = Float64 + const CMT = CM.CommonTypes const CM1 = CM.Microphysics1M const CM2 = CM.Microphysics2M const CMP = CM.Parameters -include(joinpath(pkgdir(CM), "test", "create_parameters.jl")) -toml_dict = CP.create_toml_dict(FT; dict_type = "alias") -const param_set = cloud_microphysics_parameters(toml_dict) - const rain = CMT.RainType() +const liquid = CMT.LiquidType() +const ice = CMT.IceType() +const snow = CMT.SnowType() const SB2006 = CMT.SB2006Type() const Chen2022 = CMT.Chen2022Type() +const SB2006Vel = CMT.SB2006VelType() +const Blk1MVel = CMT.Blk1MVelType() +const APS = CMP.AbstractCloudMicrophysicsParameters + +include(joinpath(pkgdir(CM), "test", "create_parameters.jl")) +toml_dict = CP.create_toml_dict(FT; dict_type = "alias") +const param_set = cloud_microphysics_parameters(toml_dict) """ - rain_terminal_velocity_individual_C(ρ, scheme, D_r) + rain_terminal_velocity_individual_Chen(param_set, ρ, D_r) - - `ρ`` - air density - - `scheme` - type for parameterization + - `param_set` - set with free parameters + - `ρ` - air density - `D_r` - diameter of the raindrops - Returns the fall velocity of a raindrop using multiple gamma-function terms - to increase accuracy for `scheme == Chen2022Type` +Returns the fall velocity of a raindrop from Chen et al 2022 """ function rain_terminal_velocity_individual_Chen( param_set, ρ::FT, - scheme::CMT.Chen2022Type, D_r::FT, #in m ) where {FT <: Real} @@ -63,6 +67,15 @@ function rain_terminal_velocity_individual_Chen( return v end +""" + rain_terminal_velocity_individual_SB(param_set, ρ, D_r) + + - `param_set` - set with free parameters + - `ρ` - air density + - `D_r` - diameter of the raindrops + +Returns the fall velocity of a raindrop from Seifert and Beheng 2006 +""" function rain_terminal_velocity_individual_SB( param_set, ρ::FT, @@ -78,95 +91,184 @@ function rain_terminal_velocity_individual_SB( return v end -q_liq_range = range(1e-8, stop = 1e-3, length = 1000) -q_rai_range = range(1e-8, stop = 1e-3, length = 1000) -D_r_range = range(50e-6, stop = 9e-3, length = 1000) -N_d_range = range(1e7, stop = 1e9, length = 1000) -q_liq = 5e-4 -# q_rai = 5e-4 -ρ_air = 1.0 # kg m^-3 -N_rai = 1e4 #this value matters for the individual terminal velocity - -PL.plot( - q_rai_range * 1e3, - [ - CM2.rain_terminal_velocity(param_set, SB2006, q_rai, ρ_air, N_rai)[1] - for q_rai in q_rai_range - ], - linewidth = 3, - xlabel = "q_rain [g/kg]", - ylabel = "terminal velocity [m/s]", - label = "Rain-SB2006 [ND]", -) -PL.plot!( - q_rai_range * 1e3, - [ - CM2.rain_terminal_velocity(param_set, Chen2022, q_rai, ρ_air, N_rai)[1] - for q_rai in q_rai_range - ], - linewidth = 3, - xlabel = "q_rain [g/kg]", - ylabel = "terminal velocity [m/s]", - label = "Rain-Chen2022 [ND]", -) -PL.plot!( - q_rai_range * 1e3, - [ - CM2.rain_terminal_velocity(param_set, SB2006, q_rai, ρ_air, N_rai)[2] - for q_rai in q_rai_range - ], - linewidth = 3, - xlabel = "q_rain [g/kg]", - ylabel = "terminal velocity [m/s]", - label = "Rain-SB2006 [M]", -) -PL.plot!( - q_rai_range * 1e3, - [ - CM2.rain_terminal_velocity(param_set, Chen2022, q_rai, ρ_air, N_rai)[2] - for q_rai in q_rai_range - ], - linewidth = 3, - xlabel = "q_rain [g/kg]", - ylabel = "terminal velocity [m/s]", - label = "Rain-Chen2022 [M]", -) -PL.plot!( - q_rai_range * 1e3, - [ - CM1.terminal_velocity(param_set, rain, ρ_air, q_rai) for - q_rai in q_rai_range - ], - linewidth = 3, - xlabel = "q_rain [g/kg]", - ylabel = "terminal velocity [m/s]", - label = "Rain-Ogura-1-moment", -) - -PL.title!("Terminal Velocity vs. Specific Humidity of Rain", titlefontsize = 9) -PL.savefig("terminal_velocity_bulk_comparisons.svg") - -PL.plot( - D_r_range, - [ - rain_terminal_velocity_individual_SB(param_set, ρ_air, D_r) for - D_r in D_r_range - ], - linewidth = 3, - xlabel = "D [m]", - ylabel = "terminal velocity [m/s]", - label = "Rain-SB2006", -) -PL.plot!( - D_r_range, - [ - rain_terminal_velocity_individual_Chen(param_set, ρ_air, Chen2022, D_r) - for D_r in D_r_range - ], - linewidth = 3, - xlabel = "D [m]", - ylabel = "terminal velocity [m/s]", - label = "Rain-Chen2022", -) -PL.title!("Terminal Velocity vs. Diameter", titlefontsize = 13) -PL.savefig("terminal_velocity_individual_raindrop.svg") +""" + terminal_velocity_individual_1M(prs, precip, ρ, D_r) + + - `prs` - set with free parameters + - `precip` - precipitation type (rain or snow) + - `ρ` - air density + - `D_r` - particle diameter + +Returns the fall velocity of a raindrop or snow from 1-moment scheme +""" +function terminal_velocity_individual_1M( + prs::APS, + precip::CMT.AbstractPrecipType, + ρ::FT, + D_r::FT, +) where {FT <: Real} + _χv::FT = CM1.χv(prs, precip) + _v0::FT = CM1.v0(prs, ρ, precip) + _ve::FT = CM1.ve(prs, precip) + _r0::FT = CM1.r0(prs, precip) + _Δv::FT = CM1.Δv(prs, precip) + vt = _χv * _v0 * (D_r / (2 * _r0))^(_Δv + _ve) + return vt +end + +""" + ice_terminal_velocity_individual_Chen(param_set, ρ, D_r) + + - `param_set` - set with free parameters + - `ρ` - air density + - `D_r` - diameter of the raindrops + +Returns the fall velocity of an ice particle from Chen et al 2022 +""" +function ice_terminal_velocity_individual_Chen( + prs::APS, + ρ::FT, + D_r::FT, #in m +) where {FT <: Real} + + ρ_i::FT = CMP.ρ_cloud_ice(prs) + + _As, _Bs, _Cs, _Es, _Fs, _Gs = CMO.Chen2022_snow_ice_coeffs(prs, ρ_i) + + ai = [_Es * ρ^_As, _Fs * ρ^_As] + bi = [_Bs + ρ * _Cs, _Bs + ρ * _Cs] + ci = [0, _Gs] + + D_r = D_r * 1000 #D_r is in mm in the paper --> multiply D_r by 1000 + v = 0 + for i in 1:2 + v += (ai[i] * (D_r)^(bi[i]) * exp(-D_r * ci[i])) + end + return v +end + +""" + snow_terminal_velocity_individual_Chen(prs, ρ, D_r) + + - `prs` - set with free parameters + - `ρ` - air density + - `D_r` - particle diameter + +Returns the fall velocity of snow from Chen et al 2022 +""" +function snow_terminal_velocity_individual_Chen( + prs::APS, + ρ::FT, + D_r::FT, +) where {FT <: Real} + D_r = D_r * 1000 + _m0::FT = CM1.m0(prs, snow) + _me::FT = CM1.me(prs, snow) + _Δm::FT = CM1.Δm(prs, snow) + _a0::FT = CM1.a0(prs, snow) + _ae::FT = CM1.ae(prs, snow) + _Δa::FT = CM1.Δa(prs, snow) + _χm::FT = CM1.χm(prs, snow) + _χa::FT = CM1.χa(prs, snow) + _r0::FT = CM1.r0(prs, snow) + ρ_i::FT = CMP.ρ_cloud_ice(prs) + m0_comb::FT = _m0 * _χm + a0_comb::FT = _a0 * _χa + me_comb::FT = _me + _Δm + ae_comb::FT = _ae + _Δa + α = FT(-1 / 3) + + _As, _Bs, _Cs, _Es, _Fs, _Gs = CMO.Chen2022_snow_ice_coeffs(prs, ρ_i) + + ai = [_Es * ρ^_As, _Fs * ρ^_As] + bi = [_Bs + ρ * _Cs, _Bs + ρ * _Cs] + ci = [0, _Gs] + + aspect_ratio = + ((16 * (a0_comb)^3 * (ρ_i)^2) / (9 * π * (m0_comb)^2)) * + (D_r / (2 * _r0 * 1000))^(3 * ae_comb - 2 * me_comb) + aspect_ratio = aspect_ratio^(α) + vt = 0 + for i in 1:2 + vt = vt + aspect_ratio * ai[i] * ((D_r))^(bi[i]) * exp(-1 * ci[i] * D_r) + end + return vt +end + +function aspect_ratio_snow_1M(prs::APS, D_r::FT) where {FT <: Real} + _m0::FT = CM1.m0(prs, snow) + _me::FT = CM1.me(prs, snow) + _Δm::FT = CM1.Δm(prs, snow) + _a0::FT = CM1.a0(prs, snow) + _ae::FT = CM1.ae(prs, snow) + _Δa::FT = CM1.Δa(prs, snow) + _χm::FT = CM1.χm(prs, snow) + _χa::FT = CM1.χa(prs, snow) + _r0::FT = CM1.r0(prs, snow) + ρ_i::FT = CMP.ρ_cloud_ice(prs) + + m0_comb::FT = _m0 * _χm + a0_comb::FT = _a0 * _χa + me_comb::FT = _me + _Δm + ae_comb::FT = _ae + _Δa + + aspect_ratio = + ((16 * (a0_comb)^3 * (ρ_i)^2) / (9 * π * (m0_comb)^2)) * + (D_r / (2 * _r0))^(3 * ae_comb - 2 * me_comb) + return aspect_ratio +end + +ρ_air, ρ_air_ground = 1.2, 1.22 +q_liq, q_ice, q_tot = 5e-4, 5e-4, 20e-3 +N_rai = 1e7 #this value matters for the individual terminal velocity +q_rain_range = range(1e-8, stop = 5e-3, length = 100) +D_r_range = range(1e-6, stop = 5e-3, length = 1000) +N_d_range = range(1e6, stop = 1e9, length = 1000) +D_r_ar_range = range(1e-6, stop = 6.25e-4, length = 1000) + +#! format: off + +SB_rain_bN = [CM2.rain_terminal_velocity(param_set, SB2006, SB2006Vel, q_rai, ρ_air, N_rai)[1] for q_rai in q_rain_range] +Ch_rain_bN = [CM2.rain_terminal_velocity(param_set, SB2006, Chen2022, q_rai, ρ_air, N_rai)[1] for q_rai in q_rain_range] +SB_rain_bM = [CM2.rain_terminal_velocity(param_set, SB2006, SB2006Vel, q_rai, ρ_air, N_rai)[2] for q_rai in q_rain_range] +Ch_rain_bM = [CM2.rain_terminal_velocity(param_set, SB2006, Chen2022, q_rai, ρ_air, N_rai)[2] for q_rai in q_rain_range] +M1_rain_bM = [CM1.terminal_velocity(param_set, rain, Blk1MVel, ρ_air, q_rai) for q_rai in q_rain_range] + +SB_rain = [rain_terminal_velocity_individual_SB(param_set, ρ_air, D_r) for D_r in D_r_range] +Ch_rain = [rain_terminal_velocity_individual_Chen(param_set, ρ_air, D_r) for D_r in D_r_range] +M1_rain = [terminal_velocity_individual_1M(param_set, rain, ρ_air, D_r) for D_r in D_r_range] +Ch_snow = [snow_terminal_velocity_individual_Chen(param_set, ρ_air, D_r) for D_r in D_r_range] +M1_snow = [terminal_velocity_individual_1M(param_set, snow, ρ_air, D_r) for D_r in D_r_range] +Ch_ice = [ice_terminal_velocity_individual_Chen(param_set, ρ_air, D_r) for D_r in D_r_range] + +Aspect_Ratio = [aspect_ratio_snow_1M(param_set, D_r) for D_r in D_r_ar_range] + +# 2 Moment Scheme Figures + +# single drop comparison +p1 = PL.plot(D_r_range *1e3, SB_rain, linewidth = 3, xlabel = "D [mm]", ylabel = "terminal velocity [m/s]", label = "Rain-SB2006", color = :blue) +p1 = PL.plot!(D_r_range * 1e3, Ch_rain, linewidth = 3, xlabel = "D [mm]", ylabel = "terminal velocity [m/s]", label = "Rain-Chen2022", color = :red) +p1 = PL.plot!(D_r_range * 1e3, M1_rain, linewidth = 3, xlabel = "D [mm]", ylabel = "terminal velocity [m/s]", label = "Rain-1M", color = :green) + +# group velocity comparison +p2 = PL.plot(q_rain_range * 1e3, SB_rain_bN, linewidth = 3, xlabel = "q_rain [g/kg]", ylabel = "terminal velocity [m/s]", label = "Rain-SB2006 [ND]", linestyle = :dot, color = :blue) +p2 = PL.plot!(q_rain_range * 1e3, Ch_rain_bN, linewidth = 3, xlabel = "q_rain [g/kg]", ylabel = "terminal velocity [m/s]", label = "Rain-Chen2022 [ND]", linestyle = :dot, color = :red) +p2 = PL.plot!(q_rain_range * 1e3, SB_rain_bM, linewidth = 3, xlabel = "q_rain [g/kg]", ylabel = "terminal velocity [m/s]", label = "Rain-SB2006 [M]", color = :blue) +p2 = PL.plot!(q_rain_range * 1e3, Ch_rain_bM, linewidth = 3, xlabel = "q_rain [g/kg]", ylabel = "terminal velocity [m/s]", label = "Rain-Chen2022 [M]", color = :red) +p2 = PL.plot!(q_rain_range * 1e3, M1_rain_bM, linewidth = 3, xlabel = "q_rain [g/kg]", ylabel = "terminal velocity [m/s]", label = "Rain-default-1M [M]", color = :green) +PL.plot(p1, p2, layout = (1, 2), size = (800, 500), dpi = 300) +PL.savefig("2M_terminal_velocity_comparisons.svg") + +# 1 Moment Scheme Figures + +# individual particle comparison +PL.plot(D_r_range * 1e3, Ch_rain, linewidth = 3, xlabel = "D [mm]", ylabel = "terminal velocity [m/s]", label = "Rain-Chen", color = :red) +PL.plot!(D_r_range * 1e3, M1_rain, linewidth = 3, xlabel = "D [mm]", ylabel = "terminal velocity [m/s]", label = "Rain-default-1M", color = :green) +PL.plot!(D_r_range * 1e3, Ch_snow, linewidth = 3, xlabel = "D [mm]", ylabel = "terminal velocity [m/s]", label = "Snow-Chen", color = :red, linestyle =:dot) +PL.plot!(D_r_range * 1e3, M1_snow, linewidth = 3, xlabel = "D [mm]", ylabel = "terminal velocity [m/s]", label = "Snow-default-1M", color = :green, linestyle = :dot) +PL.plot!(D_r_range * 1e3, Ch_ice, linewidth = 3, xlabel = "D [mm]", ylabel = "terminal velocity [m/s]", label = "Ice-Chen", color = :blue, linestyle =:dot) +PL.savefig("1M_individual_terminal_velocity_comparisons.svg") + +PL.plot(D_r_ar_range * 1e3, Aspect_Ratio, linewidth = 3, xlabel = "D [mm]", ylabel = "aspect ratio", label = "Aspect Ratio", color = :red) +PL.savefig("aspect_ratio_1M_snow.svg") +#! format: on diff --git a/src/Common.jl b/src/Common.jl index 01f3523be5..b1a81f1617 100644 --- a/src/Common.jl +++ b/src/Common.jl @@ -3,16 +3,21 @@ """ module Common -import Thermodynamics -const TD = Thermodynamics +import SpecialFunctions as SF -import ..Parameters -const CMP = Parameters -const APS = Parameters.AbstractCloudMicrophysicsParameters +import Thermodynamics as TD + +import ..Parameters as CMP +const APS = CMP.AbstractCloudMicrophysicsParameters + +import ..CommonTypes as CT export G_func export H2SO4_soln_saturation_vapor_pressure export ABIFM_Delta_a_w +export Chen_snow_ice_coeffs +export Chen2022_vel_add +export Chen2022_vel_coeffs """ G_func(param_set, T, Liquid()) @@ -121,7 +126,7 @@ end """ H2SO4_soln_saturation_vapor_pressure(x, T) - - `x` - wt percent sulphuric acid [unitless] + - `x` - wt percent sulphuric acid [unitless] - `T` - air temperature [K]. Returns the saturation vapor pressure above a sulphuric acid solution droplet in Pa. @@ -166,4 +171,118 @@ function Delta_a_w(prs::APS, x::FT, T::FT) where {FT <: Real} return min(Δa_w, FT(1)) end +""" + Chen2022_snow_ice_coeffs(prs, ρ_i) + + - prs - set with model parameters + - ρ_i - cloud ice density + +Returns the coefficients from Appendix B, Table B3 in Chen et al 2022 +DOI: 10.1016/j.atmosres.2022.106171 +needed for snow and ice terminal velocity +""" +function Chen2022_snow_ice_coeffs(prs::APS, ρ_i::FT) where {FT <: Real} + + As_1::FT = CMP.As_coeff_1_Ch2022(prs) + As_2::FT = CMP.As_coeff_2_Ch2022(prs) + As_3::FT = CMP.As_coeff_3_Ch2022(prs) + Bs_1::FT = CMP.Bs_coeff_1_Ch2022(prs) + Bs_2::FT = CMP.Bs_coeff_2_Ch2022(prs) + Bs_3::FT = CMP.Bs_coeff_3_Ch2022(prs) + Cs_1::FT = CMP.Cs_coeff_1_Ch2022(prs) + Cs_2::FT = CMP.Cs_coeff_2_Ch2022(prs) + Cs_3::FT = CMP.Cs_coeff_3_Ch2022(prs) + Cs_4::FT = CMP.Cs_coeff_4_Ch2022(prs) + Es_1::FT = CMP.Es_coeff_1_Ch2022(prs) + Es_2::FT = CMP.Es_coeff_2_Ch2022(prs) + Es_3::FT = CMP.Es_coeff_3_Ch2022(prs) + Fs_1::FT = CMP.Fs_coeff_1_Ch2022(prs) + Fs_2::FT = CMP.Fs_coeff_2_Ch2022(prs) + Fs_3::FT = CMP.Fs_coeff_3_Ch2022(prs) + Gs_1::FT = CMP.Gs_coeff_1_Ch2022(prs) + Gs_2::FT = CMP.Gs_coeff_2_Ch2022(prs) + Gs_3::FT = CMP.Gs_coeff_3_Ch2022(prs) + + As = As_1 * (log(ρ_i))^2 − As_2 * log(ρ_i) - As_3 + Bs = FT(1) / (Bs_1 + Bs_2 * log(ρ_i) + Bs_3 / sqrt(ρ_i)) + Cs = Cs_1 + Cs_2 * exp(Cs_3 * ρ_i) + Cs_4 * sqrt(ρ_i) + Es = Es_1 - Es_2 * (log(ρ_i))^2 + Es_3 * sqrt(ρ_i) + Fs = -exp(Fs_1 - Fs_2 * (log(ρ_i))^2 + Fs_3 * log(ρ_i)) + Gs = FT(1) / (Gs_1 + Gs_2 / (log(ρ_i)) - Gs_3 * log(ρ_i) / ρ_i) + + return (As, Bs, Cs, Es, Fs, Gs) +end + +""" + Chen2022_vel_coeffs(prs, precip_type, ρ) + + - prs - set with free parameters + - precip_type - type for ice, rain or snow + - ρ - air density + +Returns the coefficients from Appendix B in Chen et al 2022 +DOI: 10.1016/j.atmosres.2022.106171 +""" +function Chen2022_vel_coeffs(prs::APS, ::CT.RainType, ρ::FT) where {FT <: Real} + + ρ0::FT = CMP.q_coeff_rain_Ch2022(prs) + a1::FT = CMP.a1_coeff_rain_Ch2022(prs) + a2::FT = CMP.a2_coeff_rain_Ch2022(prs) + a3::FT = CMP.a3_coeff_rain_Ch2022(prs) + a3_pow::FT = CMP.a3_pow_coeff_rain_Ch2022(prs) + b1::FT = CMP.b1_coeff_rain_Ch2022(prs) + b2::FT = CMP.b2_coeff_rain_Ch2022(prs) + b3::FT = CMP.b3_coeff_rain_Ch2022(prs) + b_ρ::FT = CMP.b_rho_coeff_rain_Ch2022(prs) + c1::FT = CMP.c1_coeff_rain_Ch2022(prs) + c2::FT = CMP.c2_coeff_rain_Ch2022(prs) + c3::FT = CMP.c3_coeff_rain_Ch2022(prs) + + q = exp(ρ0 * ρ) + ai = (a1 * q, a2 * q, a3 * q * ρ^a3_pow) + bi = (b1 - b_ρ * ρ, b2 - b_ρ * ρ, b3 - b_ρ * ρ) + ci = (c1, c2, c3) + + # unit conversions + aiu = ai .* 1000 .^ bi + ciu = ci .* 1000 + + return (aiu, bi, ciu) +end +function Chen2022_vel_coeffs( + prs::APS, + ::Union{CT.IceType, CT.SnowType}, + ρ::FT, +) where {FT <: Real} + + ρ_i::FT = CMP.ρ_cloud_ice(prs) + + _As, _Bs, _Cs, _Es, _Fs, _Gs = Chen2022_snow_ice_coeffs(prs, ρ_i) + + ai = (_Es * ρ^_As, _Fs * ρ^_As) + bi = (_Bs + ρ * _Cs, _Bs + ρ * _Cs) + ci = (FT(0), _Gs) + # unit conversions + aiu = ai .* 1000 .^ bi + ciu = ci .* 1000 + + return (aiu, bi, ciu) +end + +""" + Chen2022_vel_add(a, b, c, λ, k) + + - a, b, c, - free parameters defined in Chen etl al 2022 + - λ - size distribution parameter + - k - size distribution moment for which we compute the bulk fall speed + +Returns the addends of the bulk fall speed of rain or ice particles +following Chen et al 2022 DOI: 10.1016/j.atmosres.2022.106171 in [m/s]. +We are assuming exponential size distribution and hence μ=0. +""" +function Chen2022_vel_add(a::FT, b::FT, c::FT, λ::FT, k::Int) where {FT <: Real} + μ = 0 # Exponantial instaed of gamma distribution + δ = μ + k + 1 + return a * λ^δ * SF.gamma(b + δ) / (λ + c)^(b + δ) / SF.gamma(δ) +end end diff --git a/src/CommonTypes.jl b/src/CommonTypes.jl index c1782f43ca..342f027c1f 100644 --- a/src/CommonTypes.jl +++ b/src/CommonTypes.jl @@ -114,12 +114,33 @@ The type for 2-moment precipitation formation by Seifert and Beheng (2006) """ struct SB2006Type <: Abstract2MPrecipType end +""" + AbstractTerminalVelocityType + +The top-level super-type for terminal velocity parameterizations +""" +abstract type AbstractTerminalVelocityType end + +""" + Blk1MVelType + +The type for precipitation terminal velocity from the simple 1-moment scheme +""" +struct Blk1MVelType <: AbstractTerminalVelocityType end + +""" + SB2006VelType + +The type for precipitation terminal velocity from Seifert and Beheng (2006) +""" +struct SB2006VelType <: AbstractTerminalVelocityType end + """ Chen2022Type -The type for 2-moment precipitation terminal velocity by Chen et al 2022 +The type for precipitation terminal velocity from Chen et. al. 2022 """ -struct Chen2022Type <: Abstract2MPrecipType end +struct Chen2022Type <: AbstractTerminalVelocityType end """ AbstractAerosolType diff --git a/src/Microphysics1M.jl b/src/Microphysics1M.jl index ea795c9db8..eeac9c6046 100644 --- a/src/Microphysics1M.jl +++ b/src/Microphysics1M.jl @@ -132,59 +132,66 @@ ve(prs::APS, ::CT.SnowType) = CMP.ve_sno(prs) """ lambda(q, ρ, n0, m0, me, r0, χm, Δm) + - `prs` - set with free parameters + - `precip` - a type for cloud ice, rain or snow - `q` - specific humidity of rain, ice or snow - `ρ` - air density - - `n0` - size distribution parameter - - `m0`, `me`, `χm`, `Δm`, `r0` - mass(radius) parameters Returns the rate parameter of the assumed size distribution of particles (rain drops, ice crystals, snow crystals). """ function lambda( + prs::APS, + precip::Union{CT.IceType, CT.RainType, CT.SnowType}, q::FT, ρ::FT, - n0::FT, - m0::FT, - me::FT, - r0::FT, - χm::FT, - Δm::FT, ) where {FT <: Real} + _n0::FT = n0(prs, q, ρ, precip) + _r0::FT = r0(prs, precip) + _m0::FT = m0(prs, precip) + _me::FT = me(prs, precip) + _Δm::FT = Δm(prs, precip) + _χm::FT = χm(prs, precip) + λ::FT = FT(0) if q > FT(0) λ = ( - χm * m0 * n0 * SF.gamma(me + Δm + FT(1)) / ρ / q / r0^(me + Δm) - )^FT(1 / (me + Δm + 1)) + _χm * _m0 * _n0 * SF.gamma(_me + _Δm + FT(1)) / ρ / q / + _r0^(_me + _Δm) + )^FT(1 / (_me + _Δm + 1)) end return λ end """ - terminal_velocity(prs, precip, ρ, q_) + terminal_velocity(prs, precip, velo_scheme, ρ, q_) - `prs` - abstract set with Earth parameters - - `precip` - a type for rain or snow + - `precip` - a type for ice, rain or snow + - `velo_scheme` - type for terminal velocity parameterization - `ρ` - air density - `q_` - rain or snow specific humidity Returns the mass weighted average terminal velocity assuming -a Marshall-Palmer (1948) distribution of rain drops and snow crystals. +a Marshall-Palmer (1948) distribution of particles. +Fall velocity of individual rain drops is parameterized: + - assuming an empirical power-law relations for `velo_scheme == Blk1MVelType` + - following Chen et. al 2022, DOI: 10.1016/j.atmosres.2022.106171, for `velo_scheme == Chen2022Type` """ function terminal_velocity( prs::APS, precip::CT.AbstractPrecipType, + velo_scheme::CT.Blk1MVelType, ρ::FT, q_::FT, ) where {FT <: Real} fall_w = FT(0) if q_ > FT(0) - _n0::FT = n0(prs, q_, ρ, precip) _r0::FT = r0(prs, precip) - _m0::FT = m0(prs, precip) _me::FT = me(prs, precip) _Δm::FT = Δm(prs, precip) _χm::FT = χm(prs, precip) @@ -192,7 +199,7 @@ function terminal_velocity( _v0::FT = v0(prs, ρ, precip) _ve::FT = ve(prs, precip) _Δv::FT = Δv(prs, precip) - _λ::FT = lambda(q_, ρ, _n0, _m0, _me, _r0, _χm, _Δm) + _λ::FT = lambda(prs, precip, q_, ρ) fall_w = _χv * @@ -204,6 +211,95 @@ function terminal_velocity( return fall_w end +function terminal_velocity( + prs::APS, + precip::CT.RainType, + velo_scheme::CT.Chen2022Type, + ρ::FT, + q_::FT, +) where {FT <: Real} + fall_w = FT(0) + if q_ > FT(0) + + # coefficients from Table B1 from Chen et. al. 2022 + aiu, bi, ciu = CO.Chen2022_vel_coeffs(prs, precip, ρ) + # size distribution parameter + _λ::FT = lambda(prs, precip, q_, ρ) + + # eq 20 from Chen et al 2022 + fall_w = sum(CO.Chen2022_vel_add.(aiu, bi, ciu, _λ, 3)) + # It should be ϕ^κ * fall_w, but for rain drops ϕ = 1 and κ = 0 + fall_w = max(FT(0), fall_w) + end + return fall_w +end +function terminal_velocity( + prs::APS, + precip::CT.IceType, + velo_scheme::CT.Chen2022Type, + ρ::FT, + q_::FT, +) where {FT <: Real} + fall_w = FT(0) + if q_ > FT(0) + + ρ_i::FT = CMP.ρ_cloud_ice(prs) + _λ::FT = lambda(prs, precip, q_, ρ) + + # coefficients from Appendix B from Chen et. al. 2022 + aiu, bi, ciu = CO.Chen2022_vel_coeffs(prs, precip, ρ) + + # eq 20 from Chen et al 2022 + fall_w = sum(CO.Chen2022_vel_add.(aiu, bi, ciu, _λ, 3)) + # It should be ϕ^κ * fall_w, but for rain drops ϕ = 1 and κ = 0 + fall_w = max(FT(0), fall_w) + end + return fall_w +end +function terminal_velocity( + prs::APS, + precip::CT.SnowType, + velo_scheme::CT.Chen2022Type, + ρ::FT, + q_::FT, +) where {FT <: Real} + fall_w = FT(0) + if q_ > FT(0) + + _r0::FT = r0(prs, precip) + _λ::FT = lambda(prs, precip, q_, ρ) + + m0c::FT = m0(prs, precip) * χm(prs, precip) + a0c::FT = a0(prs, precip) * χa(prs, precip) + mec::FT = me(prs, precip) + Δm(prs, precip) + aec::FT = ae(prs, precip) + Δa(prs, precip) + + ρ_i::FT = CMP.ρ_cloud_ice(prs) + + # coefficients from Appendix B from Chen et. al. 2022 + aiu, bi, ciu = CO.Chen2022_vel_coeffs(prs, precip, ρ) + + κ = FT(-1 / 3) #oblate + k = 3 # mass weighted + + tmp = + _λ^(k + 1) * + ((16 * a0c^3 * ρ_i^2) / (9 * π * m0c^2 * _r0^(3 * aec - 2 * mec)))^κ + ci_pow = + (2 .* ciu .+ _λ) .^ + (.-(3 .* aec .* κ .- 2 .* mec .* κ .+ bi .+ k .+ 1)) + + ti = tmp .* aiu .* FT(2) .^ bi .* ci_pow + + Chen2022_vel_add_sno(t, b, aec, mec, κ, k) = + t * SF.gamma(3 * κ * aec - 2 * κ * mec + b + k + 1) / + SF.gamma(k + 1) + + fall_w = sum(Chen2022_vel_add_sno.(ti, bi, aec, mec, κ, k)) + fall_w = max(FT(0), fall_w) + end + return fall_w +end """ conv_q_liq_to_q_rai(prs, q_liq) @@ -292,12 +388,9 @@ function conv_q_ice_to_q_sno( _r_ice_snow::FT = CMP.r_ice_snow(prs) _n0::FT = n0(prs, FT(0), ρ, CT.IceType()) - _r0::FT = r0(prs, CT.IceType()) - _m0::FT = m0(prs, CT.IceType()) _me::FT = me(prs, CT.IceType()) _Δm::FT = Δm(prs, CT.IceType()) - _χm::FT = χm(prs, CT.IceType()) - _λ::FT = lambda(q.ice, ρ, _n0, _m0, _me, _r0, _χm, _Δm) + _λ::FT = lambda(prs, CT.IceType(), q.ice, ρ) acnv_rate = 4 * FT(π) * _S * _G * _n0 / ρ * @@ -337,10 +430,6 @@ function accretion( _n0::FT = n0(prs, q_pre, ρ, precip) _r0::FT = r0(prs, precip) - _m0::FT = m0(prs, precip) - _me::FT = me(prs, precip) - _Δm::FT = Δm(prs, precip) - _χm::FT = χm(prs, precip) _χv::FT = χv(prs, precip) _v0::FT = v0(prs, ρ, precip) _ve::FT = ve(prs, precip) @@ -349,7 +438,7 @@ function accretion( _ae::FT = ae(prs, precip) _χa::FT = χa(prs, precip) _Δa::FT = Δa(prs, precip) - _λ::FT = lambda(q_pre, ρ, _n0, _m0, _me, _r0, _χm, _Δm) + _λ::FT = lambda(prs, precip, q_pre, ρ) _E::FT = E(prs, cloud, precip) accr_rate = @@ -381,13 +470,7 @@ function accretion_rain_sink( accr_rate = FT(0) if (q_ice > FT(0) && q_rai > FT(0)) - _r_ice_snow::FT = CMP.r_ice_snow(prs) _n0_ice::FT = n0(prs, FT(0), ρ, CT.IceType()) - _r0_ice::FT = r0(prs, CT.IceType()) - _m0_ice::FT = m0(prs, CT.IceType()) - _me_ice::FT = me(prs, CT.IceType()) - _Δm_ice::FT = Δm(prs, CT.IceType()) - _χm_ice::FT = χm(prs, CT.IceType()) _n0_rai::FT = n0(prs, q_rai, ρ, CT.RainType()) _r0_rai::FT = r0(prs, CT.RainType()) _m0_rai::FT = m0(prs, CT.RainType()) @@ -404,26 +487,8 @@ function accretion_rain_sink( _Δa_rai::FT = Δa(prs, CT.RainType()) _E::FT = E(prs, CT.IceType(), CT.RainType()) - _λ_rai::FT = lambda( - q_rai, - ρ, - _n0_rai, - _m0_rai, - _me_rai, - _r0_rai, - _χm_rai, - _Δm_rai, - ) - _λ_ice::FT = lambda( - q_ice, - ρ, - _n0_ice, - _m0_ice, - _me_ice, - _r0_ice, - _χm_ice, - _Δm_ice, - ) + _λ_rai::FT = lambda(prs, CT.RainType(), q_rai, ρ) + _λ_ice::FT = lambda(prs, CT.IceType(), q_ice, ρ) accr_rate = _E / ρ * @@ -479,42 +544,21 @@ function accretion_snow_rain( if (q_i > FT(0) && q_j > FT(0)) _n0_i::FT = n0(prs, q_i, ρ, type_i) - _r0_i::FT = r0(prs, type_i) - _m0_i::FT = m0(prs, type_i) - _me_i::FT = me(prs, type_i) - _Δm_i::FT = Δm(prs, type_i) - _χm_i::FT = χm(prs, type_i) - _χv_i::FT = χv(prs, type_i) - _v0_i::FT = v0(prs, ρ, type_i) - _ve_i::FT = ve(prs, type_i) - _Δv_i::FT = Δv(prs, type_i) - _a0_i::FT = a0(prs, type_i) - _ae_i::FT = ae(prs, type_i) - _χa_i::FT = χa(prs, type_i) - _Δa_i::FT = Δa(prs, type_i) - _n0_j::FT = n0(prs, q_j, ρ, type_j) + _r0_j::FT = r0(prs, type_j) _m0_j::FT = m0(prs, type_j) _me_j::FT = me(prs, type_j) _Δm_j::FT = Δm(prs, type_j) _χm_j::FT = χm(prs, type_j) - _χv_j::FT = χv(prs, type_j) - _v0_j::FT = v0(prs, ρ, type_j) - _ve_j::FT = ve(prs, type_j) - _Δv_j::FT = Δv(prs, type_j) - _a0_j::FT = a0(prs, type_j) - _ae_j::FT = ae(prs, type_j) - _χa_j::FT = χa(prs, type_j) - _Δa_j::FT = Δa(prs, type_j) _E_ij::FT = E(prs, type_i, type_j) - _λ_i::FT = lambda(q_i, ρ, _n0_i, _m0_i, _me_i, _r0_i, _χm_i, _Δm_i) - _λ_j::FT = lambda(q_j, ρ, _n0_j, _m0_j, _me_j, _r0_j, _χm_j, _Δm_j) + _λ_i::FT = lambda(prs, type_i, q_i, ρ) + _λ_j::FT = lambda(prs, type_j, q_j, ρ) - _v_ti = terminal_velocity(prs, type_i, ρ, q_i) - _v_tj = terminal_velocity(prs, type_j, ρ, q_j) + _v_ti = terminal_velocity(prs, type_i, CT.Blk1MVelType(), ρ, q_i) + _v_tj = terminal_velocity(prs, type_j, CT.Blk1MVelType(), ρ, q_j) accr_rate = FT(π) / ρ * @@ -571,23 +615,15 @@ function evaporation_sublimation( _n0::FT = n0(prs, q_rai, ρ, rain) _r0::FT = r0(prs, rain) - _m0::FT = m0(prs, rain) - _me::FT = me(prs, rain) - _Δm::FT = Δm(prs, rain) - _χm::FT = χm(prs, rain) _χv::FT = χv(prs, rain) _v0::FT = v0(prs, ρ, rain) _ve::FT = ve(prs, rain) _Δv::FT = Δv(prs, rain) - _a0::FT = a0(prs, rain) - _ae::FT = ae(prs, rain) - _χa::FT = χa(prs, rain) - _Δa::FT = Δa(prs, rain) _a_vent::FT = a_vent(prs, rain) _b_vent::FT = b_vent(prs, rain) - _λ::FT = lambda(q_rai, ρ, _n0, _m0, _me, _r0, _χm, _Δm) + _λ::FT = lambda(prs, rain, q_rai, ρ) evap_subl_rate = 4 * FT(π) * _n0 / ρ * _S * _G / _λ^FT(2) * ( @@ -620,23 +656,15 @@ function evaporation_sublimation( _n0::FT = n0(prs, q_sno, ρ, snow) _r0::FT = r0(prs, snow) - _m0::FT = m0(prs, snow) - _me::FT = me(prs, snow) - _Δm::FT = Δm(prs, snow) - _χm::FT = χm(prs, snow) _χv::FT = χv(prs, snow) _v0::FT = v0(prs, ρ, snow) _ve::FT = ve(prs, snow) _Δv::FT = Δv(prs, snow) - _a0::FT = a0(prs, snow) - _ae::FT = ae(prs, snow) - _χa::FT = χa(prs, snow) - _Δa::FT = Δa(prs, snow) _a_vent::FT = a_vent(prs, snow) _b_vent::FT = b_vent(prs, snow) - _λ::FT = lambda(q_sno, ρ, _n0, _m0, _me, _r0, _χm, _Δm) + _λ::FT = lambda(prs, snow, q_sno, ρ) evap_subl_rate = 4 * FT(π) * _n0 / ρ * _S * _G / _λ^FT(2) * ( @@ -678,23 +706,15 @@ function snow_melt(prs::APS, q_sno::FT, ρ::FT, T::FT) where {FT <: Real} _n0::FT = n0(prs, q_sno, ρ, snow) _r0::FT = r0(prs, snow) - _m0::FT = m0(prs, snow) - _me::FT = me(prs, snow) - _Δm::FT = Δm(prs, snow) - _χm::FT = χm(prs, snow) _χv::FT = χv(prs, snow) _v0::FT = v0(prs, ρ, snow) _ve::FT = ve(prs, snow) _Δv::FT = Δv(prs, snow) - _a0::FT = a0(prs, snow) - _ae::FT = ae(prs, snow) - _χa::FT = χa(prs, snow) - _Δa::FT = Δa(prs, snow) _a_vent::FT = a_vent(prs, snow) _b_vent::FT = b_vent(prs, snow) - _λ::FT = lambda(q_sno, ρ, _n0, _m0, _me, _r0, _χm, _Δm) + _λ::FT = lambda(prs, snow, q_sno, ρ) snow_melt_rate = 4 * FT(π) * _n0 / ρ * _K_therm / L * (T - _T_freeze) / _λ^FT(2) * ( diff --git a/src/Microphysics2M.jl b/src/Microphysics2M.jl index 1ebf15f9ce..94303d53a1 100644 --- a/src/Microphysics2M.jl +++ b/src/Microphysics2M.jl @@ -27,7 +27,6 @@ export rain_self_collection export rain_breakup export rain_self_collection_and_breakup export rain_terminal_velocity -export rain_terminal_velocity_reduce export rain_evaporation export conv_q_liq_to_q_rai @@ -376,26 +375,29 @@ function rain_self_collection_and_breakup( end """ - rain_terminal_velocity(param_set, scheme, q_rai, ρ, N_rai) + rain_terminal_velocity(param_set, scheme, velo_scheme, q_rai, ρ, N_rai) - `param_set` - abstract set with Earth parameters - - `scheme` - type for 2-moment liquid self-collection parameterization - - `q_rai` - rain water specific humidity - - `ρ` - air density - - `N_rai` - raindrops number density + - `scheme` - type for 2-moment parameterization + - `velo_scheme` - type for terminal velocity parameterization + - `q_rai` - rain water specific humidity [kg/kg] + - `ρ` - air density [kg/m^3] + - `N_rai` - raindrops number density [1/m^3] -Returns a tuple containing the number and mass weigthed mean fall velocities of raindrops, -assuming an empirical relation similar to Rogers (1993) for fall velocity of individual drops -and an exponential size distribution, for `scheme == SB2006Type` +Returns a tuple containing the number and mass weigthed mean fall velocities of raindrops in [m/s]. +Assuming an exponential size distribution from Seifert and Beheng 2006 for `scheme == SB2006Type` +Fall velocity of individual rain drops is parameterized: + - assuming an empirical relation similar to Rogers (1993) for `velo_scheme == SB2006VelType` + - following Chen et. al 2022, DOI: 10.1016/j.atmosres.2022.106171 for `velo_scheme == Chen2022Type` """ function rain_terminal_velocity( param_set::APS, scheme::CT.SB2006Type, + velo_scheme::CT.SB2006VelType, q_rai::FT, ρ::FT, N_rai::FT, ) where {FT <: Real} - if q_rai < eps(FT) return (FT(0), FT(0)) end @@ -411,24 +413,10 @@ function rain_terminal_velocity( return (vt0, vt1) end - -""" - rain_terminal_velocity(param_set, scheme, q_rai, ρ, N_rai) - - - `param_set` - abstract set with Earth parameters - - `scheme` - type for parameterization - - `q_rai` - rain water specific humidity [kg/kg] - - `ρ` - air density [kg/m^3] - - `N_rai` - raindrops number density [1/m^3] - -Returns a tuple containing the number and mass weigthed men fall velocities of raindrops. -Individual rain drop terminal velocity is computed based on multiple gamma-funciton terms -from Chen et. al 2022, see https://doi.org/10.1016/j.atmosres.2022.106171 -Assuming an exponential size distribution from Seifert and Beheng 2006). -""" function rain_terminal_velocity( param_set::APS, - scheme::CT.Chen2022Type, + scheme::CT.SB2006Type, + velo_scheme::CT.Chen2022Type, q_rai::FT, ρ::FT, N_rai::FT, @@ -436,39 +424,14 @@ function rain_terminal_velocity( if q_rai < eps(FT) return (FT(0), FT(0)) end - # coefficients from Table B1 from Chen et. al. 2022 - ρ0::FT = CMP.q_coeff_rain_Ch2022(param_set) - a1::FT = CMP.a1_coeff_rain_Ch2022(param_set) - a2::FT = CMP.a2_coeff_rain_Ch2022(param_set) - a3::FT = CMP.a3_coeff_rain_Ch2022(param_set) - a3_pow::FT = CMP.a3_pow_coeff_rain_Ch2022(param_set) - b1::FT = CMP.b1_coeff_rain_Ch2022(param_set) - b2::FT = CMP.b2_coeff_rain_Ch2022(param_set) - b3::FT = CMP.b3_coeff_rain_Ch2022(param_set) - b_ρ::FT = CMP.b_rho_coeff_rain_Ch2022(param_set) - c1::FT = CMP.c1_coeff_rain_Ch2022(param_set) - c2::FT = CMP.c2_coeff_rain_Ch2022(param_set) - c3::FT = CMP.c3_coeff_rain_Ch2022(param_set) - q = exp(ρ0 * ρ) - ai = (a1 * q, a2 * q, a3 * q * ρ^a3_pow) - bi = (b1 - b_ρ * ρ, b2 - b_ρ * ρ, b3 - b_ρ * ρ) - ci = (c1, c2, c3) - - #unit conversions - aiu = ai .* 1000 .^ bi - ciu = ci .* 1000 - - # size distribution coefficients + aiu, bi, ciu = CO.Chen2022_vel_coeffs(param_set, CT.RainType(), ρ) + # size distribution parameter λ = raindrops_limited_vars(param_set, q_rai, ρ, N_rai).λr - μ = FT(0) # eq 20 from Chen et al 2022 - # k = 0 for the number density, k = 3 for the mass - v_i(a, b, c, λ, δ) = - a .* λ^δ .* SF.gamma.(b .+ δ) ./ (λ .+ c) .^ (b .+ δ) ./ SF.gamma(δ) - vt0 = sum(v_i(aiu, bi, ciu, λ, μ + 1)) - vt3 = sum(v_i(aiu, bi, ciu, λ, μ + 4)) + vt0 = sum(CO.Chen2022_vel_add.(aiu, bi, ciu, λ, 0)) + vt3 = sum(CO.Chen2022_vel_add.(aiu, bi, ciu, λ, 3)) vt0 = max(FT(0), vt0) vt3 = max(FT(0), vt3) diff --git a/src/Parameters.jl b/src/Parameters.jl index 778905f308..6de925c9b1 100644 --- a/src/Parameters.jl +++ b/src/Parameters.jl @@ -135,6 +135,25 @@ Base.@kwdef struct CloudMicrophysicsParameters{FT, TP} <: c1_coeff_rain_Ch2022::FT c2_coeff_rain_Ch2022::FT c3_coeff_rain_Ch2022::FT + As_coeff_1_Ch2022::FT + As_coeff_2_Ch2022::FT + As_coeff_3_Ch2022::FT + Bs_coeff_1_Ch2022::FT + Bs_coeff_2_Ch2022::FT + Bs_coeff_3_Ch2022::FT + Cs_coeff_1_Ch2022::FT + Cs_coeff_2_Ch2022::FT + Cs_coeff_3_Ch2022::FT + Cs_coeff_4_Ch2022::FT + Es_coeff_1_Ch2022::FT + Es_coeff_2_Ch2022::FT + Es_coeff_3_Ch2022::FT + Fs_coeff_1_Ch2022::FT + Fs_coeff_2_Ch2022::FT + Fs_coeff_3_Ch2022::FT + Gs_coeff_1_Ch2022::FT + Gs_coeff_2_Ch2022::FT + Gs_coeff_3_Ch2022::FT Si_max_Mohler2006::FT T_thr_Mohler2006::FT S0_warm_ATD_Mohler2006::FT diff --git a/test/Project.toml b/test/Project.toml index 4ccfb70aa9..01b5531af6 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -6,6 +6,5 @@ CUDAKernels = "72cfdca4-0801-4ab0-bf6a-d52aa10adc57" KernelAbstractions = "63c18a36-062a-441e-b654-da1e3ab1ce7c" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" Thermodynamics = "b60c26fb-14c3-4610-9d3e-2d17fe7ff00c" - [compat] -CLIMAParameters = ">=0.7.5" +CLIMAParameters = ">=0.7.9" diff --git a/test/microphysics_tests.jl b/test/microphysics_tests.jl index 7332185628..fe676659bc 100644 --- a/test/microphysics_tests.jl +++ b/test/microphysics_tests.jl @@ -18,10 +18,13 @@ const liquid = CMT.LiquidType() const ice = CMT.IceType() const rain = CMT.RainType() const snow = CMT.SnowType() +const SB2006 = CMT.SB2006Type() const KK2000 = CMT.KK2000Type() const B1994 = CMT.B1994Type() const TC1980 = CMT.TC1980Type() const LD2004 = CMT.LD2004Type() +const Blk1MVel = CMT.Blk1MVelType() +const SB2006Vel = CMT.SB2006VelType() const Ch2022 = CMT.Chen2022Type() @info "Microphysics Tests" @@ -96,13 +99,65 @@ function test_microphysics(FT) ρ_air, q_tot, ρ_air_ground = FT(1.2), FT(20 * 1e-3), FT(1.22) for q_rai in q_rain_range - TT.@test CM1.terminal_velocity(prs, rain, ρ_air, q_rai) ≈ + TT.@test CM1.terminal_velocity(prs, rain, Blk1MVel, ρ_air, q_rai) ≈ terminal_velocity_empir(q_rai, q_tot, ρ_air, ρ_air_ground) atol = 0.2 * terminal_velocity_empir(q_rai, q_tot, ρ_air, ρ_air_ground) end end + TT.@testset "1M_microphysics - Chen 2022 rain terminal velocity" begin + #setup + ρ = FT(1.2) + q_rai = FT(5e-4) + N_rai = FT(1e4) + + #action + vt_rai = CM1.terminal_velocity(prs, rain, Ch2022, ρ, q_rai) + v_bigger = CM1.terminal_velocity(prs, rain, Ch2022, ρ, q_rai * 2) + + #test + TT.@test vt_rai ≈ 3.0721397260869705 rtol = 2e-6 + TT.@test CM1.terminal_velocity(prs, rain, Ch2022, ρ, FT(0))[1] ≈ 0 atol = + eps(FT) + TT.@test v_bigger > vt_rai + end + + TT.@testset "1M_microphysics - Chen 2022 ice terminal velocity" begin + #setup + ρ = FT(1.2) + q_ice = FT(5e-4) + N_rai = FT(1e4) + + #action + vt_ice = CM1.terminal_velocity(prs, ice, Ch2022, ρ, q_ice) + v_bigger = CM1.terminal_velocity(prs, ice, Ch2022, ρ, q_ice * 2) + + #test + TT.@test vt_ice ≈ 2.2309476335899388 rtol = 2e-6 + TT.@test CM1.terminal_velocity(prs, ice, Ch2022, ρ, FT(0)) ≈ 0 atol = + eps(FT) + TT.@test v_bigger > vt_ice + end + + TT.@testset "1M_microphysics - Chen 2022 snow terminal velocity" begin + #setup + ρ = FT(1.1) + q_sno = FT(5e-4) + N_rai = FT(1e4) + + #action + vt_sno = CM1.terminal_velocity(prs, snow, Ch2022, ρ, q_sno) + v_bigger = CM1.terminal_velocity(prs, snow, Ch2022, ρ, q_sno * 2) + + #test + TT.@test vt_sno ≈ 1.4722373984934243 rtol = 2e-6 + TT.@test CM1.terminal_velocity(prs, snow, Ch2022, ρ, FT(0)) ≈ 0 atol = + eps(FT) + TT.@test v_bigger > vt_sno + end + + TT.@testset "CloudLiquidCondEvap" begin q_liq_sat = FT(5e-3) @@ -587,17 +642,11 @@ function test_microphysics(FT) ρ0 = FT(1.225) #action - au = CM2.autoconversion(prs, CMT.SB2006Type(), q_liq, q_rai, ρ, N_liq) - sc = CM2.liquid_self_collection( - prs, - CMT.SB2006Type(), - q_liq, - ρ, - au.dN_liq_dt, - ) + au = CM2.autoconversion(prs, SB2006, q_liq, q_rai, ρ, N_liq) + sc = CM2.liquid_self_collection(prs, SB2006, q_liq, ρ, au.dN_liq_dt) au_sc = CM2.autoconversion_and_liquid_self_collection( prs, - CMT.SB2006Type(), + SB2006, q_liq, q_rai, ρ, @@ -635,17 +684,11 @@ function test_microphysics(FT) TT.@test au_sc.sc ≈ dNcdt_sc rtol = 1e-6 #action - au = CM2.autoconversion(prs, CMT.SB2006Type(), FT(0), FT(0), ρ, N_liq) - sc = CM2.liquid_self_collection( - prs, - CMT.SB2006Type(), - FT(0), - ρ, - au.dN_liq_dt, - ) + au = CM2.autoconversion(prs, SB2006, FT(0), FT(0), ρ, N_liq) + sc = CM2.liquid_self_collection(prs, SB2006, FT(0), ρ, au.dN_liq_dt) au_sc = CM2.autoconversion_and_liquid_self_collection( prs, - CMT.SB2006Type(), + SB2006, FT(0), FT(0), ρ, @@ -677,7 +720,7 @@ function test_microphysics(FT) ρ0 = FT(1.225) #action - ac = CM2.accretion(prs, CMT.SB2006Type(), q_liq, q_rai, ρ, N_liq) + ac = CM2.accretion(prs, SB2006, q_liq, q_rai, ρ, N_liq) Lc = ρ * q_liq Lr = ρ * q_rai @@ -698,7 +741,7 @@ function test_microphysics(FT) TT.@test ac.dN_rai_dt ≈ dNrdt_ac rtol = FT(1e-6) #action - ac = CM2.accretion(prs, CMT.SB2006Type(), FT(0), FT(0), ρ, N_liq) + ac = CM2.accretion(prs, SB2006, FT(0), FT(0), ρ, N_liq) #test TT.@test ac.dq_liq_dt ≈ FT(0) atol = eps(FT) @@ -722,17 +765,10 @@ function test_microphysics(FT) ρ0 = FT(1.225) #action - sc_rai = - CM2.rain_self_collection(prs, CMT.SB2006Type(), q_rai, ρ, N_rai) - br_rai = - CM2.rain_breakup(prs, CMT.SB2006Type(), q_rai, ρ, N_rai, sc_rai) - sc_br_rai = CM2.rain_self_collection_and_breakup( - prs, - CMT.SB2006Type(), - q_rai, - ρ, - N_rai, - ) + sc_rai = CM2.rain_self_collection(prs, SB2006, q_rai, ρ, N_rai) + br_rai = CM2.rain_breakup(prs, SB2006, q_rai, ρ, N_rai, sc_rai) + sc_br_rai = + CM2.rain_self_collection_and_breakup(prs, SB2006, q_rai, ρ, N_rai) λr = CM2.raindrops_limited_vars(prs, q_rai, ρ, N_rai).λr * @@ -753,13 +789,8 @@ function test_microphysics(FT) #test TT.@test sc_rai ≈ dNrdt_sc rtol = 1e-6 - TT.@test CM2.rain_self_collection( - prs, - CMT.SB2006Type(), - FT(0), - ρ, - N_rai, - ) ≈ FT(0) atol = eps(FT) + TT.@test CM2.rain_self_collection(prs, SB2006, 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 @@ -769,17 +800,10 @@ function test_microphysics(FT) q_rai = FT(0) #action - sc_rai = - CM2.rain_self_collection(prs, CMT.SB2006Type(), q_rai, ρ, N_rai) - br_rai = - CM2.rain_breakup(prs, CMT.SB2006Type(), q_rai, ρ, N_rai, sc_rai) - sc_br_rai = CM2.rain_self_collection_and_breakup( - prs, - CMT.SB2006Type(), - q_rai, - ρ, - N_rai, - ) + sc_rai = CM2.rain_self_collection(prs, SB2006, q_rai, ρ, N_rai) + br_rai = CM2.rain_breakup(prs, SB2006, q_rai, ρ, N_rai, sc_rai) + sc_br_rai = + CM2.rain_self_collection_and_breakup(prs, SB2006, q_rai, ρ, N_rai) #test TT.@test sc_rai ≈ FT(0) atol = eps(FT) @@ -801,7 +825,7 @@ function test_microphysics(FT) #action vt_rai = - CM2.rain_terminal_velocity(prs, CMT.SB2006Type(), q_rai, ρ, N_rai) + CM2.rain_terminal_velocity(prs, SB2006, SB2006Vel, q_rai, ρ, N_rai) λr = CM2.raindrops_limited_vars(prs, q_rai, ρ, N_rai).λr vt0 = max(0, sqrt(ρ0 / ρ) * (aR - bR / (1 + cR / λr))) @@ -813,14 +837,16 @@ function test_microphysics(FT) TT.@test vt_rai[2] ≈ vt1 rtol = 1e-6 TT.@test CM2.rain_terminal_velocity( prs, - CMT.SB2006Type(), + SB2006, + SB2006Vel, FT(0), ρ, N_rai, )[1] ≈ 0 atol = eps(FT) TT.@test CM2.rain_terminal_velocity( prs, - CMT.SB2006Type(), + SB2006, + SB2006Vel, FT(0), ρ, N_rai, @@ -834,18 +860,32 @@ function test_microphysics(FT) N_rai = FT(1e4) #action - vt_rai = CM2.rain_terminal_velocity(prs, Ch2022, q_rai, ρ, N_rai) - v_bigger = CM2.rain_terminal_velocity(prs, Ch2022, q_rai * 2, ρ, N_rai) + vt_rai = + CM2.rain_terminal_velocity(prs, SB2006, Ch2022, q_rai, ρ, N_rai) + v_bigger = + CM2.rain_terminal_velocity(prs, SB2006, Ch2022, q_rai * 2, ρ, N_rai) #test TT.@test vt_rai isa Tuple TT.@test vt_rai[1] ≈ 1.2591475834547752 rtol = 1e-6 TT.@test vt_rai[2] ≈ 4.552478635185714 rtol = 1e-6 - TT.@test CM2.rain_terminal_velocity(prs, Ch2022, FT(0), ρ, N_rai)[1] ≈ 0 atol = - eps(FT) - TT.@test CM2.rain_terminal_velocity(prs, Ch2022, FT(0), ρ, N_rai)[2] ≈ 0 atol = - eps(FT) + TT.@test CM2.rain_terminal_velocity( + prs, + SB2006, + Ch2022, + FT(0), + ρ, + N_rai, + )[1] ≈ 0 atol = eps(FT) + TT.@test CM2.rain_terminal_velocity( + prs, + SB2006, + Ch2022, + 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] @@ -869,8 +909,7 @@ function test_microphysics(FT) ρ0 = FT(1.225) #action - evap = - CM2.rain_evaporation(prs, CMT.SB2006Type(), q, q_rai, ρ, N_rai, T) + evap = CM2.rain_evaporation(prs, SB2006, q, q_rai, ρ, N_rai, T) G = CMC.G_func(prs, T, TD.Liquid()) thermo_params = CMP.thermodynamics_params(prs) @@ -894,24 +933,10 @@ function test_microphysics(FT) TT.@test evap isa Tuple TT.@test evap[1] ≈ evap0 rtol = 1e-5 TT.@test evap[2] ≈ evap1 rtol = 1e-5 - TT.@test CM2.rain_evaporation( - prs, - CMT.SB2006Type(), - q, - FT(0), - ρ, - N_rai, - T, - )[1] ≈ 0 atol = eps(FT) - TT.@test CM2.rain_evaporation( - prs, - CMT.SB2006Type(), - q, - FT(0), - ρ, - N_rai, - T, - )[2] ≈ 0 atol = eps(FT) + TT.@test CM2.rain_evaporation(prs, SB2006, q, FT(0), ρ, N_rai, T)[1] ≈ 0 atol = + eps(FT) + TT.@test CM2.rain_evaporation(prs, SB2006, q, FT(0), ρ, N_rai, T)[2] ≈ 0 atol = + eps(FT) end end diff --git a/test/performance_tests.jl b/test/performance_tests.jl index 2e5d877c10..12741e01cc 100644 --- a/test/performance_tests.jl +++ b/test/performance_tests.jl @@ -45,6 +45,7 @@ function benchmark_test(FT) liquid = CMT.LiquidType() rain = CMT.RainType() sb2006 = CMT.SB2006Type() + sb2006vel = CMT.SB2006VelType() ch2022 = CMT.Chen2022Type() dust = CMT.DesertDustType() @@ -136,12 +137,12 @@ function benchmark_test(FT) ) bench_press( CM2.rain_terminal_velocity, - (prs, sb2006, q_rai, ρ_air, N_rai), + (prs, sb2006, sb2006vel, q_rai, ρ_air, N_rai), 300, ) bench_press( CM2.rain_terminal_velocity, - (prs, ch2022, q_rai, ρ_air, N_rai), + (prs, sb2006, ch2022, q_rai, ρ_air, N_rai), 1700, ) # Homogeneous Nucleation