Add precipitation heating terms

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

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(

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)
 

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!(