diff --git a/src/cache/prognostic_edmf_precomputed_quantities.jl b/src/cache/prognostic_edmf_precomputed_quantities.jl index df85dbcb16..a4ed36a42c 100644 --- a/src/cache/prognostic_edmf_precomputed_quantities.jl +++ b/src/cache/prognostic_edmf_precomputed_quantities.jl @@ -256,6 +256,7 @@ function set_prognostic_edmf_precomputed_quantities_closures!(Y, p, t) ) end + # TODO add the 1-moment microphysics option here @. ᶜS_q_tot⁰ = q_tot_precipitation_sources( precip_model, thermo_params, diff --git a/src/parameterized_tendencies/microphysics/microphysics_wrappers.jl b/src/parameterized_tendencies/microphysics/microphysics_wrappers.jl index c2d85a9590..4e9c3807cc 100644 --- a/src/parameterized_tendencies/microphysics/microphysics_wrappers.jl +++ b/src/parameterized_tendencies/microphysics/microphysics_wrappers.jl @@ -3,64 +3,280 @@ import Thermodynamics as TD import CloudMicrophysics.Microphysics0M as CM0 +# define some aliases and functions to make the code more readable +const Iₗ = TD.internal_energy_liquid +const Iᵢ = TD.internal_energy_ice +const Lf = TD.latent_heat_fusion +const Tₐ = TD.air_temperature +const PP = TD.PhasePartition +const qᵥ = TD.vapor_specific_humidity +qₗ(thp, ts) = TD.PhasePartition(thp, ts).liq +qᵢ(thp, ts) = TD.PhasePartition(thp, ts).ice + +# helper function to safely get precipitation from state +function qₚ(ρqₚ::FT, ρ::FT) where {FT} + return max(FT(0), ρqₚ / ρ) +end + +# helper function to limit the tendency +function limit(q::FT, dt::FT) where {FT} + return q / dt / 5 +end + """ - q_tot_precipitation_sources(precip_model, t_prs, m_prs, dt, q_tot, ts) + q_tot_precipitation_sources(precip_model, thp, cmp, dt, qₜ, ts) - precip_model - a type for precipitation scheme choice - - t_prs, m_prs - stricts with thermodynamic and microphysics parameters + - thp, cmp - structs with thermodynamic and microphysics parameters - dt - model time step - - q_tot - total water specific humidity - - ts - thermodynamic state (see td package for details) + - qₜ - total water specific humidity + - ts - thermodynamic state (see Thermodynamics.jl package for details) -Returns the q_tot source term due to precipitation formation +Returns the qₜ source term due to precipitation formation defined as Δm_tot / (m_dry + m_tot) """ function q_tot_precipitation_sources( ::NoPrecipitation, - t_prs, - m_prs, + thp, + cmp, dt, - q_tot::FT, + qₜ::FT, ts, ) where {FT <: Real} return FT(0) end function q_tot_precipitation_sources( ::Microphysics0Moment, - t_prs, - m_prs, + thp, + cmp, dt, - q_tot::FT, + qₜ::FT, ts, ) where {FT <: Real} return -min( - max(q_tot, 0) / dt, + max(qₜ, 0) / dt, -CM0.remove_precipitation( - m_prs, - TD.PhasePartition(t_prs, ts), - TD.q_vap_saturation(t_prs, ts), + cmp, + PP(thp, ts), + TD.q_vap_saturation(thp, ts), ), ) end """ - e_tot_0M_precipitation_sources_helper(t_prs, ts, Φ) + e_tot_0M_precipitation_sources_helper(thp, ts, Φ) - - t_prs - set with thermodynamics parameters + - thp - set with thermodynamics parameters - ts - thermodynamic state (see td package for details) - Φ - geopotential Returns the total energy source term multiplier from precipitation formation +for the 0-moment scheme """ function e_tot_0M_precipitation_sources_helper( - t_prs, + thp, ts, Φ::FT, ) where {FT <: Real} - λ::FT = TD.liquid_fraction(t_prs, ts) - I_l::FT = TD.internal_energy_liquid(t_prs, ts) - I_i::FT = TD.internal_energy_ice(t_prs, ts) + λ::FT = TD.liquid_fraction(thp, ts) + + return λ * Iₗ(thp, ts) + (1 - λ) * Iᵢ(thp, ts) + Φ +end + +""" + compute_precipitation_sources!(Sᵖ, Sᵖ_snow, Sqₜᵖ, Sqᵣᵖ, Sqₛᵖ, Seₜᵖ, ρ, ρqᵣ, ρqₛ, ts, Φ, dt, mp, thp) + + - Sᵖ, Sᵖ_snow - temporary containters to help compute precipitation source terms + - Sqₜᵖ, Sqᵣᵖ, Sqₛᵖ, Seₜᵖ - cached storage for precipitation source terms + - ρ - air density + - ρqᵣ, ρqₛ - precipitation (rain and snow) densities + - ts - thermodynamic state (see td package for details) + - Φ - geopotential + - dt - model time step + - thp, cmp - structs with thermodynamic and microphysics parameters + +Returns the q source terms due to precipitation formation from the 1-moment scheme. +The specific humidity source terms are defined as defined as Δmᵢ / (m_dry + m_tot) +where i stands for total, rain or snow. +Also returns the total energy source term due to the microphysics processes. +""" +function compute_precipitation_sources!( + Sᵖ, + Sᵖ_snow, + Sqₜᵖ, + Sqᵣᵖ, + Sqₛᵖ, + Seₜᵖ, + ρ, + ρqᵣ, + ρqₛ, + ts, + Φ, + dt, + mp, + thp, +) + FT = eltype(Sqₜᵖ) + #! format: off + # rain autoconversion: q_liq -> q_rain + @. Sᵖ = min( + limit(qₗ(thp, ts), dt), + CM1.conv_q_liq_to_q_rai(mp.pr.acnv1M, qₗ(thp, ts), true), + ) + @. Sqₜᵖ -= Sᵖ + @. Sqᵣᵖ += Sᵖ + @. Seₜᵖ -= Sᵖ * (Iₗ(thp, ts) + Φ) + + # snow autoconversion assuming no supersaturation: q_ice -> q_snow + @. Sᵖ = min( + limit(qᵢ(thp, ts), dt), + CM1.conv_q_ice_to_q_sno_no_supersat(mp.ps.acnv1M, qᵢ(thp, ts), true), + ) + @. Sqₜᵖ -= Sᵖ + @. Sqₛᵖ += Sᵖ + @. Seₜᵖ -= Sᵖ * (Iᵢ(thp, ts) + Φ) + + # accretion: q_liq + q_rain -> q_rain + @. Sᵖ = min( + limit(qₗ(thp, ts), dt), + CM1.accretion(mp.cl, mp.pr, mp.tv.rain, mp.ce, qₗ(thp, ts), qₚ(ρqᵣ, ρ), ρ), + ) + @. Sqₜᵖ -= Sᵖ + @. Sqᵣᵖ += Sᵖ + @. Seₜᵖ -= Sᵖ * (Iₗ(thp, ts) + Φ) + + # accretion: q_ice + q_snow -> q_snow + @. Sᵖ = min( + limit(qᵢ(thp, ts), dt), + CM1.accretion(mp.ci, mp.ps, mp.tv.snow, mp.ce, qᵢ(thp, ts), qₚ(ρqₛ, ρ), ρ), + ) + @. Sqₜᵖ -= Sᵖ + @. Sqₛᵖ += Sᵖ + @. Seₜᵖ -= Sᵖ * (Iᵢ(thp, ts) + Φ) + + # accretion: q_liq + q_sno -> q_sno or q_rai + # sink of cloud water via accretion cloud water + snow + @. Sᵖ = min( + limit(qₗ(thp, ts), dt), + CM1.accretion(mp.cl, mp.ps, mp.tv.snow, mp.ce, qₗ(thp, ts), qₚ(ρqₛ, ρ), ρ), + ) + # 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) + @. Sᵖ_snow = ifelse( + Tₐ(thp, ts) < mp.ps.T_freeze, + Sᵖ, + FT(-1) * min(Sᵖ * α(thp, ts), limit(qₚ(ρqₛ, ρ), dt)), + ) + @. Sqₛᵖ += Sᵖ_snow + @. Sqₜᵖ -= Sᵖ + @. Sqᵣᵖ += ifelse(Tₐ(thp, ts) < mp.ps.T_freeze, FT(0), Sᵖ - Sᵖ_snow) + @. Seₜᵖ -= ifelse( + Tₐ(thp, ts) < mp.ps.T_freeze, + Sᵖ * (Iᵢ(thp, ts) + Φ), + Sᵖ * (Iₗ(thp, ts) + Φ) - Sᵖ_snow * (Iₗ(thp, ts) - Iᵢ(thp, ts)), + ) + + # accretion: q_ice + q_rai -> q_sno + @. Sᵖ = min( + limit(qᵢ(thp, ts), dt), + CM1.accretion(mp.ci, mp.pr, mp.tv.rain, mp.ce, qᵢ(thp, ts), qₚ(ρqᵣ, ρ), ρ), + ) + @. Sqₜᵖ -= Sᵖ + @. Sqₛᵖ += Sᵖ + @. Seₜᵖ -= Sᵖ * (Iᵢ(thp, ts) + Φ) + # sink of rain via accretion cloud ice - rain + @. Sᵖ = min( + limit(qₚ(ρqᵣ, ρ), dt), + CM1.accretion_rain_sink(mp.pr, mp.ci, mp.tv.rain, mp.ce, qᵢ(thp, ts), qₚ(ρqᵣ, ρ), ρ), + ) + @. Sqᵣᵖ -= Sᵖ + @. Sqₛᵖ += Sᵖ + @. Seₜᵖ += Sᵖ * Lf(thp, ts) + + # accretion: q_rai + q_sno -> q_rai or q_sno + @. Sᵖ = ifelse( + Tₐ(thp, ts) < mp.ps.T_freeze, + min( + limit(qₚ(ρqᵣ, ρ), dt), + CM1.accretion_snow_rain(mp.ps, mp.pr, mp.tv.rain, mp.tv.snow, mp.ce, qₚ(ρqₛ, ρ), qₚ(ρqᵣ, ρ), ρ), + ), + -min( + limit(qₚ(ρqₛ, ρ), dt), + CM1.accretion_snow_rain(mp.pr, mp.ps, mp.tv.snow, mp.tv.rain, mp.ce, qₚ(ρqᵣ, ρ), qₚ(ρqₛ, ρ), ρ), + ), + ) + @. Sqₛᵖ += Sᵖ + @. Sqᵣᵖ -= Sᵖ + @. Seₜᵖ += Sᵖ * Lf(thp, ts) + #! format: on +end + +""" + compute_precipitation_sinks!(Sᵖ, Sqₜᵖ, Sqᵣᵖ, Sqₛᵖ, Seₜᵖ, ρ, ρqᵣ, ρqₛ, ts, Φ, dt, mp, thp) - return λ * I_l + (1 - λ) * I_i + Φ + - Sᵖ - a temporary containter to help compute precipitation source terms + - Sqₜᵖ, Sqᵣᵖ, Sqₛᵖ, Seₜᵖ - cached storage for precipitation source terms + - ρ - air density + - ρqᵣ, ρqₛ - precipitation (rain and snow) densities + - ts - thermodynamic state (see td package for details) + - Φ - geopotential + - dt - model time step + - thp, cmp - structs with thermodynamic and microphysics parameters + +Returns the q source terms due to precipitation sinks from the 1-moment scheme. +The specific humidity source terms are defined as defined as Δmᵢ / (m_dry + m_tot) +where i stands for total, rain or snow. +Also returns the total energy source term due to the microphysics processes. +""" +function compute_precipitation_sinks!( + Sᵖ, + Sqₜᵖ, + Sqᵣᵖ, + Sqₛᵖ, + Seₜᵖ, + ρ, + ρqᵣ, + ρqₛ, + ts, + Φ, + dt, + mp, + thp, +) + FT = eltype(Sqₜᵖ) + sps = (mp.ps, mp.tv.snow, mp.aps, thp) + rps = (mp.pr, mp.tv.rain, mp.aps, thp) + + #! format: off + # evaporation: q_rai -> q_vap + @. Sᵖ = -min( + limit(qₚ(ρqᵣ, ρ), dt), + -CM1.evaporation_sublimation(rps..., PP(thp, ts), qₚ(ρqᵣ, ρ), ρ, Tₐ(thp, ts)), + ) + @. Sqₜᵖ -= Sᵖ + @. Sqᵣᵖ += Sᵖ + @. Seₜᵖ -= Sᵖ * (Iₗ(thp, ts) + Φ) + + # melting: q_sno -> q_rai + @. Sᵖ = min( + limit(qₚ(ρqₛ, ρ), dt), + CM1.snow_melt(sps..., qₚ(ρqₛ, ρ), ρ, Tₐ(thp, ts)), + ) + @. Sqᵣᵖ += Sᵖ + @. Sqₛᵖ -= Sᵖ + @. Seₜᵖ -= Sᵖ * Lf(thp, ts) + + # deposition/sublimation: q_vap <-> q_sno + @. Sᵖ = CM1.evaporation_sublimation(sps..., PP(thp, ts), qₚ(ρqₛ, ρ), ρ, Tₐ(thp, ts)) + @. Sᵖ = ifelse( + Sᵖ > FT(0), + min(limit(qᵥ(thp, ts), dt), Sᵖ), + -min(limit(qₚ(ρqₛ, ρ), dt), FT(-1) * Sᵖ), + ) + @. Sqₜᵖ -= Sᵖ + @. Sqₛᵖ += Sᵖ + @. Seₜᵖ -= Sᵖ * (Iᵢ(thp, ts) + Φ) + #! format: on end diff --git a/src/parameterized_tendencies/microphysics/precipitation.jl b/src/parameterized_tendencies/microphysics/precipitation.jl index 0e84ad11e9..e820c16fe6 100644 --- a/src/parameterized_tendencies/microphysics/precipitation.jl +++ b/src/parameterized_tendencies/microphysics/precipitation.jl @@ -21,7 +21,7 @@ precipitation_cache(Y, precip_model::NoPrecipitation) = (;) precipitation_tendency!(Yₜ, Y, p, t, colidx, ::NoPrecipitation) = nothing ##### -##### 0-Moment without sgs scheme or with diagnostic edmf +##### 0-Moment without sgs scheme or with diagnostic/prognostic edmf ##### function precipitation_cache(Y, precip_model::Microphysics0Moment) @@ -40,16 +40,19 @@ function precipitation_cache(Y, precip_model::Microphysics0Moment) end function compute_precipitation_cache!(Y, p, colidx, ::Microphysics0Moment, _) - (; params) = p + (; params, dt) = p (; ᶜts) = p.precomputed (; ᶜS_ρq_tot) = p.precipitation cm_params = CAP.microphysics_params(params) thermo_params = CAP.thermodynamics_params(params) - #TODO missing limiting by q_tot/Δt @. ᶜS_ρq_tot[colidx] = - Y.c.ρ[colidx] * CM.Microphysics0M.remove_precipitation( + Y.c.ρ[colidx] * q_tot_precipitation_sources( + Microphysics0Moment(), + thermo_params, cm_params, - TD.PhasePartition(thermo_params, ᶜts[colidx]), + dt, + Y.c.ρq_tot[colidx] / Y.c.ρ[colidx], + ᶜts[colidx], ) end @@ -185,268 +188,62 @@ function precipitation_cache(Y, precip_model::Microphysics1Moment) ) end -function qₚ(ρqₚ::FT, ρ::FT) where {FT} - return max(FT(0), ρqₚ / ρ) -end - -function limit(q::FT, dt::FT) where {FT} - return q / dt / 5 -end - function compute_precipitation_cache!(Y, p, colidx, ::Microphysics1Moment, _) FT = Spaces.undertype(axes(Y.c)) - (; params) = p (; dt) = p (; ᶜts) = p.precomputed (; ᶜΦ) = p.core (; ᶜSqₜᵖ, ᶜSqᵣᵖ, ᶜSqₛᵖ, ᶜSeₜᵖ) = p.precipitation - ᶜz = Fields.coordinate_field(Y.c).z - ᶜSᵖ = p.scratch.ᶜtemp_scalar + ᶜSᵖ_snow = p.scratch.ᶜtemp_scalar_2 # get thermodynamics and 1-moment microphysics params + (; params) = p cmp = CAP.microphysics_params(params) thp = CAP.thermodynamics_params(params) - # some helper functions to make the code more readable - Iₗ(ts) = TD.internal_energy_liquid(thp, ts) - Iᵢ(ts) = TD.internal_energy_ice(thp, ts) - qₗ(ts) = TD.PhasePartition(thp, ts).liq - qᵢ(ts) = TD.PhasePartition(thp, ts).ice - qᵥ(ts) = TD.vapor_specific_humidity(thp, ts) - Lf(ts) = TD.latent_heat_fusion(thp, ts) - Tₐ(ts) = TD.air_temperature(thp, ts) - α(ts) = TD.Parameters.cv_l(thp) / Lf(ts) * (Tₐ(ts) - cmp.ps.T_freeze) - - # zero out the source terms + # zero out the helper source terms @. ᶜSqₜᵖ[colidx] = FT(0) - @. ᶜSeₜᵖ[colidx] = FT(0) @. ᶜSqᵣᵖ[colidx] = FT(0) @. ᶜSqₛᵖ[colidx] = FT(0) + @. ᶜSeₜᵖ[colidx] = FT(0) - # All the tendencies are individually limited - # by the available condensate (q_ / dt). - - # rain autoconversion: q_liq -> q_rain - @. ᶜSᵖ[colidx] = min( - limit(qₗ(ᶜts[colidx]), dt), - CM1.conv_q_liq_to_q_rai(cmp.pr.acnv1M, qₗ(ᶜts[colidx]), true), - ) - @. ᶜSqₜᵖ[colidx] -= ᶜSᵖ[colidx] - @. ᶜSqᵣᵖ[colidx] += ᶜSᵖ[colidx] - @. ᶜSeₜᵖ[colidx] -= ᶜSᵖ[colidx] * (Iₗ(ᶜts[colidx]) + ᶜΦ[colidx]) - - # snow autoconversion assuming no supersaturation: q_ice -> q_snow - @. ᶜSᵖ[colidx] = min( - limit(qᵢ(ᶜts[colidx]), dt), - CM1.conv_q_ice_to_q_sno_no_supersat( - cmp.ps.acnv1M, - qᵢ(ᶜts[colidx]), - true, - ), - ) - @. ᶜSqₜᵖ[colidx] -= ᶜSᵖ[colidx] - @. ᶜSqₛᵖ[colidx] += ᶜSᵖ[colidx] - @. ᶜSeₜᵖ[colidx] -= ᶜSᵖ[colidx] * (Iᵢ(ᶜts[colidx]) + ᶜΦ[colidx]) - - # accretion: q_liq + q_rain -> q_rain - @. ᶜSᵖ[colidx] = min( - limit(qₗ(ᶜts[colidx]), dt), - CM1.accretion( - cmp.cl, - cmp.pr, - cmp.tv.rain, - cmp.ce, - qₗ(ᶜts[colidx]), - qₚ(Y.c.ρq_rai[colidx], Y.c.ρ[colidx]), - Y.c.ρ[colidx], - ), - ) - @. ᶜSqₜᵖ[colidx] -= ᶜSᵖ[colidx] - @. ᶜSqᵣᵖ[colidx] += ᶜSᵖ[colidx] - @. ᶜSeₜᵖ[colidx] -= ᶜSᵖ[colidx] * (Iₗ(ᶜts[colidx]) + ᶜΦ[colidx]) - - # accretion: q_ice + q_snow -> q_snow - @. ᶜSᵖ[colidx] = min( - limit(qᵢ(ᶜts[colidx]), dt), - CM1.accretion( - cmp.ci, - cmp.ps, - cmp.tv.snow, - cmp.ce, - qᵢ(ᶜts[colidx]), - qₚ(Y.c.ρq_sno[colidx], Y.c.ρ[colidx]), - Y.c.ρ[colidx], - ), - ) - @. ᶜSqₜᵖ[colidx] -= ᶜSᵖ[colidx] - @. ᶜSqₛᵖ[colidx] += ᶜSᵖ[colidx] - @. ᶜSeₜᵖ[colidx] -= ᶜSᵖ[colidx] * (Iᵢ(ᶜts[colidx]) + ᶜΦ[colidx]) - - # accretion: q_liq + q_sno -> q_sno or q_rai - # sink of cloud water via accretion cloud water + snow - @. ᶜSᵖ[colidx] = min( - limit(qₗ(ᶜts[colidx]), dt), - CM1.accretion( - cmp.cl, - cmp.ps, - cmp.tv.snow, - cmp.ce, - qₗ(ᶜts[colidx]), - qₚ(Y.c.ρq_sno[colidx], Y.c.ρ[colidx]), - Y.c.ρ[colidx], - ), - ) - # if T < T_freeze cloud droplets freeze to become snow - # else the snow melts and both cloud water and snow become rain - ᶜSᵖ_snow = p.scratch.ᶜtemp_scalar_2 - @. ᶜSᵖ_snow[colidx] = ifelse( - Tₐ(ᶜts[colidx]) < cmp.ps.T_freeze, + # compute precipitation source terms + # TODO - need version of this for EDMF SGS + compute_precipitation_sources!( ᶜSᵖ[colidx], - FT(-1) * min( - ᶜSᵖ[colidx] * α(ᶜts[colidx]), - limit(qₚ(Y.c.ρq_sno[colidx], Y.c.ρ[colidx]), dt), - ), - ) - @. ᶜSqₛᵖ[colidx] += ᶜSᵖ_snow[colidx] - @. ᶜSqₜᵖ[colidx] -= ᶜSᵖ[colidx] - @. ᶜSqᵣᵖ[colidx] += ifelse( - Tₐ(ᶜts[colidx]) < cmp.ps.T_freeze, - FT(0), - ᶜSᵖ[colidx] - ᶜSᵖ_snow[colidx], - ) - @. ᶜSeₜᵖ[colidx] -= ifelse( - Tₐ(ᶜts[colidx]) < cmp.ps.T_freeze, - ᶜSᵖ[colidx] * (Iᵢ(ᶜts[colidx]) + ᶜΦ[colidx]), - ᶜSᵖ[colidx] * (Iₗ(ᶜts[colidx]) + ᶜΦ[colidx]) - - ᶜSᵖ_snow[colidx] * (Iₗ(ᶜts[colidx]) - Iᵢ(ᶜts[colidx])), - ) - - # accretion: q_ice + q_rai -> q_sno - @. ᶜSᵖ[colidx] = min( - limit(qᵢ(ᶜts[colidx]), dt), - CM1.accretion( - cmp.ci, - cmp.pr, - cmp.tv.rain, - cmp.ce, - qᵢ(ᶜts[colidx]), - qₚ(Y.c.ρq_rai[colidx], Y.c.ρ[colidx]), - Y.c.ρ[colidx], - ), - ) - @. ᶜSqₜᵖ[colidx] -= ᶜSᵖ[colidx] - @. ᶜSqₛᵖ[colidx] += ᶜSᵖ[colidx] - @. ᶜSeₜᵖ[colidx] -= ᶜSᵖ[colidx] * (Iᵢ(ᶜts[colidx]) + ᶜΦ[colidx]) - # sink of rain via accretion cloud ice - rain - @. ᶜSᵖ[colidx] = min( - limit(qₚ(Y.c.ρq_rai[colidx], Y.c.ρ[colidx]), dt), - CM1.accretion_rain_sink( - cmp.pr, - cmp.ci, - cmp.tv.rain, - cmp.ce, - qᵢ(ᶜts[colidx]), - qₚ(Y.c.ρq_rai[colidx], Y.c.ρ[colidx]), - Y.c.ρ[colidx], - ), - ) - @. ᶜSqᵣᵖ[colidx] -= ᶜSᵖ[colidx] - @. ᶜSqₛᵖ[colidx] += ᶜSᵖ[colidx] - @. ᶜSeₜᵖ[colidx] += ᶜSᵖ[colidx] * Lf(ᶜts[colidx]) - - # accretion: q_rai + q_sno -> q_rai or q_sno - @. ᶜSᵖ[colidx] = ifelse( - Tₐ(ᶜts[colidx]) < cmp.ps.T_freeze, - min( - limit(qₚ(Y.c.ρq_rai[colidx], Y.c.ρ[colidx]), dt), - CM1.accretion_snow_rain( - cmp.ps, - cmp.pr, - cmp.tv.rain, - cmp.tv.snow, - cmp.ce, - qₚ(Y.c.ρq_sno[colidx], Y.c.ρ[colidx]), - qₚ(Y.c.ρq_rai[colidx], Y.c.ρ[colidx]), - Y.c.ρ[colidx], - ), - ), - -min( - limit(qₚ(Y.c.ρq_sno[colidx], Y.c.ρ[colidx]), dt), - CM1.accretion_snow_rain( - cmp.pr, - cmp.ps, - cmp.tv.snow, - cmp.tv.rain, - cmp.ce, - qₚ(Y.c.ρq_rai[colidx], Y.c.ρ[colidx]), - qₚ(Y.c.ρq_sno[colidx], Y.c.ρ[colidx]), - Y.c.ρ[colidx], - ), - ), - ) - @. ᶜSqₛᵖ[colidx] += ᶜSᵖ[colidx] - @. ᶜSqᵣᵖ[colidx] -= ᶜSᵖ[colidx] - @. ᶜSeₜᵖ[colidx] += ᶜSᵖ[colidx] * Lf(ᶜts[colidx]) - - # evaporation: q_rai -> q_vap - @. ᶜSᵖ[colidx] = - -min( - limit(qₚ(Y.c.ρq_rai[colidx], Y.c.ρ[colidx]), dt), - -CM1.evaporation_sublimation( - cmp.pr, - cmp.tv.rain, - cmp.aps, - thp, - TD.PhasePartition(thp, ᶜts[colidx]), - qₚ(Y.c.ρq_rai[colidx], Y.c.ρ[colidx]), - Y.c.ρ[colidx], - Tₐ(ᶜts[colidx]), - ), - ) - @. ᶜSqₜᵖ[colidx] -= ᶜSᵖ[colidx] - @. ᶜSqᵣᵖ[colidx] += ᶜSᵖ[colidx] - @. ᶜSeₜᵖ[colidx] -= ᶜSᵖ[colidx] * (Iₗ(ᶜts[colidx]) + ᶜΦ[colidx]) - - # melting: q_sno -> q_rai - @. ᶜSᵖ[colidx] = min( - limit(qₚ(Y.c.ρq_sno[colidx], Y.c.ρ[colidx]), dt), - CM1.snow_melt( - cmp.ps, - cmp.tv.snow, - cmp.aps, - thp, - qₚ(Y.c.ρq_sno[colidx], Y.c.ρ[colidx]), - Y.c.ρ[colidx], - Tₐ(ᶜts[colidx]), - ), + ᶜSᵖ_snow[colidx], + ᶜSqₜᵖ[colidx], + ᶜSqᵣᵖ[colidx], + ᶜSqₛᵖ[colidx], + ᶜSeₜᵖ[colidx], + Y.c.ρ[colidx], + Y.c.ρq_rai[colidx], + Y.c.ρq_sno[colidx], + ᶜts[colidx], + ᶜΦ[colidx], + dt, + cmp, + thp, ) - @. ᶜSqᵣᵖ[colidx] += ᶜSᵖ[colidx] - @. ᶜSqₛᵖ[colidx] -= ᶜSᵖ[colidx] - @. ᶜSeₜᵖ[colidx] -= ᶜSᵖ[colidx] * Lf(ᶜts[colidx]) - # deposition/sublimation: q_vap <-> q_sno - @. ᶜSᵖ[colidx] = CM1.evaporation_sublimation( - cmp.ps, - cmp.tv.snow, - cmp.aps, - thp, - TD.PhasePartition(thp, ᶜts[colidx]), - qₚ(Y.c.ρq_sno[colidx], Y.c.ρ[colidx]), + # compute precipitation sinks + # (For now only done on the grid mean) + compute_precipitation_sinks!( + ᶜSᵖ[colidx], + ᶜSqₜᵖ[colidx], + ᶜSqᵣᵖ[colidx], + ᶜSqₛᵖ[colidx], + ᶜSeₜᵖ[colidx], Y.c.ρ[colidx], - Tₐ(ᶜts[colidx]), - ) - @. ᶜSᵖ[colidx] = ifelse( - ᶜSᵖ[colidx] > FT(0), - min(limit(qᵥ(ᶜts[colidx]), dt), ᶜSᵖ[colidx]), - -min( - limit(qₚ(Y.c.ρq_sno[colidx], Y.c.ρ[colidx]), dt), - FT(-1) * ᶜSᵖ[colidx], - ), + Y.c.ρq_rai[colidx], + Y.c.ρq_sno[colidx], + ᶜts[colidx], + ᶜΦ[colidx], + dt, + cmp, + thp, ) - @. ᶜSqₜᵖ[colidx] -= ᶜSᵖ[colidx] - @. ᶜSqₛᵖ[colidx] += ᶜSᵖ[colidx] - @. ᶜSeₜᵖ[colidx] -= ᶜSᵖ[colidx] * (Iᵢ(ᶜts[colidx]) + ᶜΦ[colidx]) end function precipitation_tendency!( @@ -457,17 +254,12 @@ function precipitation_tendency!( colidx, precip_model::Microphysics1Moment, ) - FT = Spaces.undertype(axes(Y.c)) - (; ᶜSqₜᵖ, ᶜSqᵣᵖ, ᶜSqₛᵖ, ᶜSeₜᵖ) = p.precipitation + (; turbconv_model) = p.atmos + compute_precipitation_cache!(Y, p, colidx, precip_model, turbconv_model) - compute_precipitation_cache!( - Y, - p, - colidx, - precip_model, - p.atmos.turbconv_model, - ) + (; ᶜSqₜᵖ, ᶜSqᵣᵖ, ᶜSqₛᵖ, ᶜSeₜᵖ) = p.precipitation + # Update grid mean tendencies @. Yₜ.c.ρ[colidx] += Y.c.ρ[colidx] * ᶜSqₜᵖ[colidx] @. Yₜ.c.ρq_tot[colidx] += Y.c.ρ[colidx] * ᶜSqₜᵖ[colidx] @. Yₜ.c.ρe_tot[colidx] += Y.c.ρ[colidx] * ᶜSeₜᵖ[colidx] diff --git a/src/prognostic_equations/remaining_tendency.jl b/src/prognostic_equations/remaining_tendency.jl index dd87abe296..92dbb89942 100644 --- a/src/prognostic_equations/remaining_tendency.jl +++ b/src/prognostic_equations/remaining_tendency.jl @@ -81,6 +81,7 @@ NVTX.@annotate function additional_tendency!(Yₜ, Y, p, t) p.atmos.turbconv_model, ) edmfx_tke_tendency!(Yₜ, Y, p, t, colidx, p.atmos.turbconv_model) + # TODO - add the 1-moment precipitation microphysics here edmfx_precipitation_tendency!( Yₜ, Y,