diff --git a/docs/src/equations.md b/docs/src/equations.md index 68e14a00fc..2358177ce0 100644 --- a/docs/src/equations.md +++ b/docs/src/equations.md @@ -382,6 +382,23 @@ It is assummed that some fraction ``\alpha`` of snow is melted during the proces \frac{d}{dt} \rho e = \rho \mathcal{S}_{acc} ((1+\alpha) I_{liq} - \alpha I_{ice} + \Phi) ``` +### Precipitation temperature + +Precipitation is assumed to have the same temperature as the surrounding air ``T_a``. +The corresponding energy sink associated with heat transfer + between air and precipitating species can be written as + +```math +\frac{d}{dt} \rho e = - \rho q_p (\boldsymbol{u} - v_p) c_p \nabla T_a +``` +where ``q_p``, ``\boldsymbol{u}``, ``v_p``, ``c_p `` are the precipitation specific humidity, +air velocity, precipitation terminal velocity assumed to be along the gravity axis, +specific heat of precipitating species. + +!!! todo + We should consider replacing ``T_a`` with + some approximation of wet bulb temperature. + ### Stability and positivity All source terms are individually limited such that they don't exceed the diff --git a/src/cache/temporary_quantities.jl b/src/cache/temporary_quantities.jl index 49fafaa23f..04356c41c2 100644 --- a/src/cache/temporary_quantities.jl +++ b/src/cache/temporary_quantities.jl @@ -35,6 +35,7 @@ function temporary_quantities(Y, atmos) ᶜtemp_C12 = Fields.Field(C12{FT}, center_space), # ᶜuₕ_mean ᶜtemp_C3 = Fields.Field(C3{FT}, center_space), # ᶜ∇Φ₃ ᶜtemp_CT3 = Fields.Field(CT3{FT}, center_space), # ᶜω³, ᶜ∇Φ³ + ᶜtemp_CT123 = Fields.Field(CT123{FT}, center_space), ᶠtemp_CT3 = Fields.Field(CT3{FT}, face_space), # ᶠuₕ³ ᶠtemp_CT12 = Fields.Field(CT12{FT}, face_space), # ᶠω¹² ᶠtemp_CT12ʲs = Fields.Field( diff --git a/src/parameterized_tendencies/microphysics/microphysics_wrappers.jl b/src/parameterized_tendencies/microphysics/microphysics_wrappers.jl index e543bdaa0b..5c73a0735f 100644 --- a/src/parameterized_tendencies/microphysics/microphysics_wrappers.jl +++ b/src/parameterized_tendencies/microphysics/microphysics_wrappers.jl @@ -16,6 +16,8 @@ const qᵥ = TD.vapor_specific_humidity qₜ(thp, ts) = TD.PhasePartition(thp, ts).tot qₗ(thp, ts) = TD.PhasePartition(thp, ts).liq qᵢ(thp, ts) = TD.PhasePartition(thp, ts).ice +cᵥₗ(thp) = TD.Parameters.cv_l(thp) +cᵥᵢ(thp) = TD.Parameters.cv_i(thp) # helper function to limit the tendency function limit(q::FT, dt::FT, n::Int) where {FT} @@ -211,7 +213,7 @@ function compute_precipitation_sources!( ) # if T < T_freeze cloud droplets freeze to become snow # else the snow melts and both cloud water and snow become rain - α(thp, ts) = TD.Parameters.cv_l(thp) / Lf(thp, ts) * (Tₐ(thp, ts) - mp.ps.T_freeze) + α(thp, ts) = cᵥₗ(thp) / Lf(thp, ts) * (Tₐ(thp, ts) - mp.ps.T_freeze) @. Sᵖ_snow = ifelse( Tₐ(thp, ts) < mp.ps.T_freeze, Sᵖ, @@ -261,6 +263,44 @@ function compute_precipitation_sources!( #! format: on end +""" + compute_precipitation_heating(Seₜᵖ, ᶜwᵣ, ᶜwₛ, ᶜu, qᵣ, qₛ, ᶜts, thp) + + - Seₜᵖ - cached storage for precipitation energy source terms + - ᶜwᵣ, ᶜwₛ - rain and snow terminal velocities + - ᶜu - air velocity + - qᵣ, qₛ - precipitation (rain and snow) specific humidities + - ᶜts - thermodynamic state (see td package for details) + - ᶜ∇T - cached temporary variable to store temperature gradient + - thp - structs with thermodynamic and microphysics parameters + + Augments the energy source terms with heat exchange between air + and precipitating species, following eq. 36 from Raymond 2013 + doi:10.1002/jame.20050 and assuming that precipitation has the same + temperature as the surrounding air. +""" +function compute_precipitation_heating!( + ᶜSeₜᵖ, + ᶜwᵣ, + ᶜwₛ, + ᶜu, + ᶜqᵣ, + ᶜqₛ, + ᶜts, + ᶜ∇T, + thp, +) + # TODO - at some point we want to switch to assuming that precipitation + # is at wet bulb temperature + + # compute full temperature gradient + @. ᶜ∇T = CT123(ᶜgradᵥ(ᶠinterp(Tₐ(thp, ᶜts)))) + @. ᶜ∇T += CT123(gradₕ(Tₐ(thp, ᶜts))) + # dot product with effective velocity of precipitation + # (times q and specific heat) + @. ᶜSeₜᵖ -= dot(ᶜ∇T, (ᶜu - C123(Geometry.WVector(ᶜwᵣ)))) * cᵥₗ(thp) * ᶜqᵣ + @. ᶜSeₜᵖ -= dot(ᶜ∇T, (ᶜu - C123(Geometry.WVector(ᶜwₛ)))) * cᵥᵢ(thp) * ᶜqₛ +end """ compute_precipitation_sinks!(Sᵖ, Sqₜᵖ, Sqᵣᵖ, Sqₛᵖ, Seₜᵖ, ρ, qᵣ, qₛ, ts, Φ, dt, mp, thp) diff --git a/src/parameterized_tendencies/microphysics/precipitation.jl b/src/parameterized_tendencies/microphysics/precipitation.jl index 10ad9ffac7..a2b9db6273 100644 --- a/src/parameterized_tendencies/microphysics/precipitation.jl +++ b/src/parameterized_tendencies/microphysics/precipitation.jl @@ -184,11 +184,13 @@ end function compute_precipitation_cache!(Y, p, ::Microphysics1Moment, _) FT = Spaces.undertype(axes(Y.c)) (; dt) = p - (; ᶜts, ᶜqᵣ, ᶜqₛ) = p.precomputed + (; ᶜts, ᶜqᵣ, ᶜqₛ, ᶜwᵣ, ᶜwₛ, ᶜu) = p.precomputed (; ᶜΦ) = p.core (; ᶜSqₜᵖ, ᶜSqᵣᵖ, ᶜSqₛᵖ, ᶜSeₜᵖ) = p.precipitation + ᶜSᵖ = p.scratch.ᶜtemp_scalar ᶜSᵖ_snow = p.scratch.ᶜtemp_scalar_2 + ᶜ∇T = p.scratch.ᶜtemp_CT123 # get thermodynamics and 1-moment microphysics params (; params) = p @@ -230,7 +232,10 @@ function compute_precipitation_cache!(Y, p, ::Microphysics1Moment, _) cmp, thp, ) + # first term of eq 36 from Raymond 2013 + compute_precipitation_heating!(ᶜSeₜᵖ, ᶜwᵣ, ᶜwₛ, ᶜu, ᶜqᵣ, ᶜqₛ, ᶜts, ᶜ∇T, thp) end + function compute_precipitation_cache!( Y, p, @@ -239,12 +244,13 @@ function compute_precipitation_cache!( ) FT = Spaces.undertype(axes(Y.c)) (; dt) = p - (; ᶜts, ᶜqᵣ, ᶜqₛ) = p.precomputed + (; ᶜts, ᶜqᵣ, ᶜqₛ, ᶜwᵣ, ᶜwₛ, ᶜu) = p.precomputed (; ᶜΦ) = p.core # Grid mean precipitation sinks (; ᶜSqₜᵖ, ᶜSqᵣᵖ, ᶜSqₛᵖ, ᶜSeₜᵖ) = p.precipitation # additional scratch storage ᶜSᵖ = p.scratch.ᶜtemp_scalar + ᶜ∇T = p.scratch.ᶜtemp_CT123 # get thermodynamics and 1-moment microphysics params (; params) = p @@ -273,6 +279,8 @@ function compute_precipitation_cache!( cmp, thp, ) + # first term of eq 36 from Raymond 2013 + compute_precipitation_heating!(ᶜSeₜᵖ, ᶜwᵣ, ᶜwₛ, ᶜu, ᶜqᵣ, ᶜqₛ, ᶜts, ᶜ∇T, thp) end function precipitation_tendency!(