From ece33a12b092f5900cf00212adec84e8bf361b67 Mon Sep 17 00:00:00 2001 From: Anna Jaruga Date: Mon, 20 Nov 2023 09:33:11 -0800 Subject: [PATCH] Add 1-moment microphysics --- docs/src/equations.md | 118 +++++++- examples/hybrid/driver.jl | 8 +- src/ClimaAtmos.jl | 1 + .../precipitation_precomputed_quantities.jl | 25 ++ src/cache/precomputed_quantities.jl | 10 +- .../microphysics/precipitation.jl | 253 +++++++++++++++++- .../implicit/implicit_tendency.jl | 4 +- 7 files changed, 401 insertions(+), 18 deletions(-) create mode 100644 src/cache/precipitation_precomputed_quantities.jl diff --git a/docs/src/equations.md b/docs/src/equations.md index 62f073a34f..d306947149 100644 --- a/docs/src/equations.md +++ b/docs/src/equations.md @@ -100,13 +100,13 @@ We make use of the following operators Follows the continuity equation ```math -\frac{\partial}{\partial t} \rho = - \nabla \cdot(\rho \boldsymbol{u}) . +\frac{\partial}{\partial t} \rho = - \nabla \cdot(\rho \boldsymbol{u}) + \rho \mathcal{S}_{qt}. ``` This is discretized using the following ```math \frac{\partial}{\partial t} \rho -= - \hat{\mathcal{D}}_h[ \rho \bar{\boldsymbol{u}}] - \mathcal{D}^c_v \left[WI^f( J, \rho) \tilde{\boldsymbol{u}} \right] += - \hat{\mathcal{D}}_h[ \rho \bar{\boldsymbol{u}}] - \mathcal{D}^c_v \left[WI^f( J, \rho) \tilde{\boldsymbol{u}} \right] + \rho \mathcal{S}_{qt} ``` with the @@ -139,7 +139,7 @@ This is stabilized with the addition of 4th-order vector hyperviscosity ```math -\nu_u \, \nabla_h^2 (\nabla_h^2(\boldsymbol{\overbar{u}})), ``` -projected onto the first two contravariant directions, where ``\nabla_{h}^2(\boldsymbol{v})`` is the horizontal vector Laplacian. For grid scale hyperdiffusion, ``\boldsymbol{v}`` is identical to ``\boldsymbol{\overbar{u}}``, the cell-center valued velocity vector. +projected onto the first two contravariant directions, where ``\nabla_{h}^2(\boldsymbol{v})`` is the horizontal vector Laplacian. For grid scale hyperdiffusion, ``\boldsymbol{v}`` is identical to ``\boldsymbol{\overbar{u}}``, the cell-center valued velocity vector. ```math \nabla_h^2(\boldsymbol{v}) = \nabla_h(\nabla_{h} \cdot \boldsymbol{v}) - \nabla_{h} \times (\nabla_{h} \times \boldsymbol{v}). ``` @@ -199,7 +199,7 @@ projected onto the third contravariant direction. ### Total energy ```math -\frac{\partial}{\partial t} \rho e = - \nabla \cdot((\rho e + p) \boldsymbol{u} + \boldsymbol{F}_R), +\frac{\partial}{\partial t} \rho e = - \nabla \cdot((\rho e + p) \boldsymbol{u} + \boldsymbol{F}_R) + \rho \mathcal{S}_{e}, ``` which is stabilized with the addition of a 4th-order hyperdiffusion term on total enthalpy: ```math @@ -236,7 +236,7 @@ is treated implicitly. For an arbitrary scalar ``\chi``, the density-weighted scalar ``\rho\chi`` follows the continuity equation ```math -\frac{\partial}{\partial t} \rho \chi = - \nabla \cdot(\rho \chi \boldsymbol{u}). +\frac{\partial}{\partial t} \rho \chi = - \nabla \cdot(\rho \chi \boldsymbol{u}) + \rho \mathcal{S}_{\chi}. ``` This is stabilized with the addition of a 4th-order hyperdiffusion term ```math @@ -267,3 +267,111 @@ is treated implicitly. - \mathcal{D}^c_v\left[WI^f(J, \rho) U^f\left(I^f(\boldsymbol{u}_h) + \boldsymbol{u}_v, \frac{\rho \chi}{\rho} \right) \right] ``` is treated implicitly. + +## Microphysics source terms + +Sources from cloud microphysics ``\mathcal{S}`` represent the transfer of mass + between the working fluid (dry air, water vapor cloud liquid and cloud ice) + and precipitation (rain and snow), + as well as the latent heat release due to phase changes. + +The scalars ``\rho q_{rai}`` and ``\rho q_{sno}`` are part of the state vector + when running simulations with 1-moment microphysics scheme, + and represent the specific humidity of liquid and solid precipitation + (i.e. rain and snow). + +```math +q_{rai} := \frac{m_{rai}}{m_{dry} + m_{vap} + m_{liq} + m_{ice}}\; , +\;\;\;\; +q_{sno} := \frac{m_{sno}}{m_{dry} + m_{vap} + m_{liq} + m_{ice}} +``` + +The different source terms are provided by + [CloudMicrophysics.jl](https://github.com/CliMA/CloudMicrophysics.jl) library + and are defined as the change of mass of one of the cloud condensate or + precipitation species normalised by the mass of the working fluid. +See the [CloudMicrophysics.jl docs](https://clima.github.io/CloudMicrophysics.jl/dev/) + for more details. + +!!! todo + Throughout the rest of the derivations we are assuming that the volume + of the working fluid is constant (not the pressure). + This is strange for phase changes and needs more thinking. + +### Case 1: Mass of the working fluid is changed + +When the phase change is happening within the working fluid + (for example condensation from water vapor to liquid water), + there is no change to any of the state variables. +Considering the transition from + ``x \rightarrow y`` where ``x`` is either + water vapor, cloud liquid water or cloud ice and + ``y`` is either rain or snow +```math +\mathcal{S}_{x \rightarrow y} := \frac{\frac{dm_x}{dt}}{m_{dry} + m_{vap} + m_{liq} + m_{ice}} +``` +```math +\frac{d}{dt} \rho = +\frac{d}{dt} \rho q_{tot} = +\rho \mathcal{S}_{x \rightarrow y} = +- \frac{d}{dt} \rho q_y +``` +```math +\frac{d}{dt} \rho e = \rho \mathcal{S}_{x \rightarrow y} (I_{y} + \Phi) +``` +where ``I_{y}`` is the internal energy of the ``y`` phase. +This formula applies to the majority of microphysics processes. +Namely, it is valid for processes where ``T=const`` such as + autoconversion and accretion between species of the same phase. +It is also valid for rain evaporation, deposition/sublimation, and + accretion of cloud water and snow in temperatures below freezing + (which result in snow). + +### Case 2: Phase change outside of the working fluid + +For cases where both ``x`` and ``y`` are not part of the working fluid + (melting of snow, freezing of rain) +```math +\mathcal{S}_{x \rightarrow y} := \frac{\frac{dm_{x}}{dt}}{m_{dry} + m_{vap} + m_{liq} + m_{ice}} +``` +```math +\frac{d}{dt} \rho q_{x} = +- \frac{d}{dt} \rho q_{y} = +\rho \mathcal{S}_{x \rightarrow y} +``` +```math +\frac{d}{dt} \rho = \frac{d}{dt} \rho q_{tot} = 0 +``` +```math +\frac{d}{dt} \rho e = - \rho \mathcal{S}_{x \rightarrow y} L_{f} +``` +where ``L_f`` is the latent heat of fusion. +The sign in the last equation assumes ``x`` stands for rain and ``y`` for snow. + +### Additional cases + +Accretion of cloud ice by rain results in snow. +This process combines the effects from the loss of working fluid ``q_{ice}`` + (described by case 1) + and the phase change from rain to snow + (described by case 2). + +Accretion of cloud liquid by snow in temperatures above freezing results in rain. +It is assummed that some fraction ``\alpha`` of snow is melted during the process + and both cloud liquid and melted snow are turned into rain. + +```math +\mathcal{S}_{acc} := \frac{\frac{dm_{liq}}{dt}}{m_{dry} + m_{vap} + m_{liq} + m_{ice}} +``` +```math +\frac{d}{dt} \rho = \frac{d}{dt} \rho q_{tot} = \rho S_{acc} +``` +```math +\frac{d}{dt} \rho q_{sno} = \rho \alpha S_{acc} +``` +```math +\frac{d}{dt} \rho q_{rai} = - \rho (1 + \alpha) S_{acc} +``` +```math +\frac{d}{dt} \rho e = \rho \mathcal{S}_{acc} ((1+\alpha) I_{liq} - \alpha I_{ice} + \Phi) +``` diff --git a/examples/hybrid/driver.jl b/examples/hybrid/driver.jl index 1c1f5e52ef..7e289d0041 100644 --- a/examples/hybrid/driver.jl +++ b/examples/hybrid/driver.jl @@ -152,12 +152,12 @@ if config.parsed_args["check_precipitation"] @assert !any(isnan, Yₜ.c.ρe_tot[colidx]) @assert !any(isnan, Yₜ.c.ρq_rai[colidx]) @assert !any(isnan, Yₜ.c.ρq_sno[colidx]) - @assert !any(isnan, sol.prob.p.precipitation.ᶜwᵣ[colidx]) - @assert !any(isnan, sol.prob.p.precipitation.ᶜwₛ[colidx]) + @assert !any(isnan, sol.prob.p.precomputed.ᶜwᵣ[colidx]) + @assert !any(isnan, sol.prob.p.precomputed.ᶜwₛ[colidx]) # treminal velocity is positive - @test minimum(sol.prob.p.precipitation.ᶜwᵣ[colidx]) >= FT(0) - @test minimum(sol.prob.p.precipitation.ᶜwₛ[colidx]) >= FT(0) + @test minimum(sol.prob.p.precomputed.ᶜwᵣ[colidx]) >= FT(0) + @test minimum(sol.prob.p.precomputed.ᶜwₛ[colidx]) >= FT(0) # checking for water budget conservation # in the presence of precipitation sinks diff --git a/src/ClimaAtmos.jl b/src/ClimaAtmos.jl index 0999585a9d..2e62dfefce 100644 --- a/src/ClimaAtmos.jl +++ b/src/ClimaAtmos.jl @@ -23,6 +23,7 @@ include(joinpath("parameterized_tendencies", "radiation", "radiation.jl")) include(joinpath("cache", "prognostic_edmf_precomputed_quantities.jl")) include(joinpath("cache", "diagnostic_edmf_precomputed_quantities.jl")) +include(joinpath("cache", "precipitation_precomputed_quantities.jl")) include(joinpath("cache", "precomputed_quantities.jl")) include(joinpath("initial_conditions", "InitialConditions.jl")) diff --git a/src/cache/precipitation_precomputed_quantities.jl b/src/cache/precipitation_precomputed_quantities.jl new file mode 100644 index 0000000000..24b74e5a45 --- /dev/null +++ b/src/cache/precipitation_precomputed_quantities.jl @@ -0,0 +1,25 @@ +##### +##### Precomputed quantities +##### +import CloudMicrophysics.Microphysics1M as CM1 + +""" + set_precipitation_precomputed_quantities!(Y, p, t) + +Updates the precipitation terminal velocity stored in `p` +for the 1-moment microphysics scheme +""" +function set_precipitation_precomputed_quantities!(Y, p, t) + @assert (p.atmos.precip_model isa Microphysics1Moment) + + (; ᶜwᵣ, ᶜwₛ) = p.precomputed + + cmp = CAP.microphysics_params(p.params) + + # compute the precipitation terminal velocity [m/s] + @. ᶜwᵣ = + CM1.terminal_velocity(cmp.pr, cmp.tv.rain, Y.c.ρ, Y.c.ρq_rai / Y.c.ρ) + @. ᶜwₛ = + CM1.terminal_velocity(cmp.ps, cmp.tv.snow, Y.c.ρ, Y.c.ρq_sno / Y.c.ρ) + return nothing +end diff --git a/src/cache/precomputed_quantities.jl b/src/cache/precomputed_quantities.jl index ddb73bde86..876bac8bfc 100644 --- a/src/cache/precomputed_quantities.jl +++ b/src/cache/precomputed_quantities.jl @@ -109,10 +109,14 @@ function precomputed_quantities(Y, atmos) ᶜK_h = similar(Y.c, FT), ρatke_flux = similar(Fields.level(Y.f, half), C3{FT}), ) : (;) + precipitation_quantities = + atmos.precip_model isa Microphysics1Moment ? + (; ᶜwᵣ = similar(Y.c, FT), ᶜwₛ = similar(Y.c, FT)) : (;) return (; gs_quantities..., advective_sgs_quantities..., diagnostic_sgs_quantities..., + precipitation_quantities..., ) end @@ -285,7 +289,7 @@ function instead of recomputing the value yourself. Otherwise, it will be difficult to ensure that the duplicated computations are consistent. """ NVTX.@annotate function set_precomputed_quantities!(Y, p, t) - (; moisture_model, turbconv_model) = p.atmos + (; moisture_model, turbconv_model, precip_model) = p.atmos thermo_params = CAP.thermodynamics_params(p.params) n = n_mass_flux_subdomains(turbconv_model) thermo_args = (thermo_params, moisture_model) @@ -340,6 +344,10 @@ NVTX.@annotate function set_precomputed_quantities!(Y, p, t) set_diagnostic_edmf_precomputed_quantities_env_closures!(Y, p, t) end + if precip_model isa Microphysics1Moment + set_precipitation_precomputed_quantities!(Y, p, t) + end + return nothing end diff --git a/src/parameterized_tendencies/microphysics/precipitation.jl b/src/parameterized_tendencies/microphysics/precipitation.jl index d64271f665..cf1753927b 100644 --- a/src/parameterized_tendencies/microphysics/precipitation.jl +++ b/src/parameterized_tendencies/microphysics/precipitation.jl @@ -156,23 +156,264 @@ function precipitation_cache(Y, precip_model::Microphysics1Moment) ᶜSqᵣᵖ = similar(Y.c, FT), ᶜSqₛᵖ = similar(Y.c, FT), ᶜSeₜᵖ = similar(Y.c, FT), - ᶜwᵣ = similar(Y.c, FT), - ᶜwₛ = similar(Y.c, FT), ) end function compute_precipitation_cache!(Y, p, colidx, ::Microphysics1Moment, _) FT = Spaces.undertype(axes(Y.c)) + (; params) = p + (; dt) = p.simulation + (; ᶜts) = p.precomputed + (; ᶜΦ) = p.core + (; ᶜSqₜᵖ, ᶜSqᵣᵖ, ᶜSqₛᵖ, ᶜSeₜᵖ) = p.precipitation + ᶜz = Fields.coordinate_field(Y.c).z + + ᶜSᵖ = p.scratch.ᶜtemp_scalar - (; ᶜSqₜᵖ, ᶜSqᵣᵖ, ᶜSqₛᵖ, ᶜSeₜᵖ, ᶜwᵣ, ᶜwₛ) = p.precipitation + # get thermodynamics and 1-moment microphysics params + cmp = CAP.microphysics_params(params) + thp = CAP.thermodynamics_params(params) - # TODO in the next PR - @. ᶜwᵣ[colidx] = FT(0) - @. ᶜwₛ[colidx] = FT(0) + # 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 @. ᶜSqₜᵖ[colidx] = FT(0) @. ᶜSeₜᵖ[colidx] = FT(0) @. ᶜSqᵣᵖ[colidx] = FT(0) @. ᶜSqₛᵖ[colidx] = FT(0) + + # All the tendencies are individually limited + # by the available condensate (q_ / dt). + + # rain autoconversion: q_liq -> q_rain + @. ᶜSᵖ[colidx] = min( + qₗ(ᶜts[colidx]) / dt, + CM1.conv_q_liq_to_q_rai( + cmp.pr.acnv1M, + qₗ(ᶜts[colidx]), + smooth_transition = 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( + qᵢ(ᶜts[colidx]) / dt, + CM1.conv_q_ice_to_q_sno_no_supersat( + cmp.ps.acnv1M, + qᵢ(ᶜts[colidx]), + smooth_transition = 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( + qₗ(ᶜts[colidx]) / dt, + CM1.accretion( + cmp.cl, + cmp.pr, + cmp.tv.rain, + cmp.ce, + qₗ(ᶜts[colidx]), + 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( + qᵢ(ᶜts[colidx]) / dt, + CM1.accretion( + cmp.ci, + cmp.ps, + cmp.tv.snow, + cmp.ce, + qᵢ(ᶜts[colidx]), + 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( + qₗ(ᶜts[colidx]) / dt, + CM1.accretion( + cmp.cl, + cmp.ps, + cmp.tv.snow, + cmp.ce, + qₗ(ᶜts[colidx]), + 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, + ᶜSᵖ[colidx], + FT(-1) * min( + ᶜSᵖ[colidx] * α(ᶜts[colidx]), + 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( + qᵢ(ᶜts[colidx]) / dt, + CM1.accretion( + cmp.ci, + cmp.pr, + cmp.tv.rain, + cmp.ce, + qᵢ(ᶜts[colidx]), + Y.c.ρq_rai[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( + 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]), + 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( + 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, + Y.c.ρq_sno[colidx] / Y.c.ρ[colidx], + Y.c.ρq_rai[colidx] / Y.c.ρ[colidx], + Y.c.ρ[colidx], + ), + ), + -min( + 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, + Y.c.ρq_rai[colidx] / Y.c.ρ[colidx], + 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( + 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]), + 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( + Y.c.ρq_sno[colidx] / Y.c.ρ[colidx] / dt, + CM1.snow_melt( + cmp.ps, + cmp.tv.snow, + cmp.aps, + thp, + Y.c.ρq_sno[colidx] / Y.c.ρ[colidx], + Y.c.ρ[colidx], + Tₐ(ᶜts[colidx]), + ), + ) + @. ᶜ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]), + Y.c.ρq_sno[colidx] / Y.c.ρ[colidx], + Y.c.ρ[colidx], + Tₐ(ᶜts[colidx]), + ) + @. ᶜSᵖ[colidx] = ifelse( + ᶜSᵖ[colidx] > FT(0), + min(qᵥ(ᶜts[colidx]) / dt, ᶜSᵖ[colidx]), + -min(Y.c.ρq_sno[colidx] / Y.c.ρ[colidx] / dt, FT(-1) * ᶜSᵖ[colidx]), + ) + @. ᶜSqₜᵖ[colidx] -= ᶜSᵖ[colidx] + @. ᶜSqₛᵖ[colidx] += ᶜSᵖ[colidx] + @. ᶜSeₜᵖ[colidx] -= ᶜSᵖ[colidx] * (Iᵢ(ᶜts[colidx]) + ᶜΦ[colidx]) end function precipitation_tendency!( diff --git a/src/prognostic_equations/implicit/implicit_tendency.jl b/src/prognostic_equations/implicit/implicit_tendency.jl index 5a4db2cef3..086992c508 100644 --- a/src/prognostic_equations/implicit/implicit_tendency.jl +++ b/src/prognostic_equations/implicit/implicit_tendency.jl @@ -118,7 +118,7 @@ function implicit_vertical_advection_tendency!(Yₜ, Y, p, t, colidx) @. ᶠu³ₚ[colidx] = FT(-1) * - ᶠinterp(p.precipitation.ᶜwᵣ[colidx]) * + ᶠinterp(p.precomputed.ᶜwᵣ[colidx]) * CT3(unit_basis_vector_data(CT3, lgf[colidx])) @. ᶜqₚ[colidx] = Y.c.ρq_rai[colidx] / Y.c.ρ[colidx] @@ -140,7 +140,7 @@ function implicit_vertical_advection_tendency!(Yₜ, Y, p, t, colidx) @. ᶠu³ₚ[colidx] = FT(-1) * - ᶠinterp(p.precipitation.ᶜwₛ[colidx]) * + ᶠinterp(p.precomputed.ᶜwₛ[colidx]) * CT3(unit_basis_vector_data(CT3, lgf[colidx])) @. ᶜqₚ[colidx] = Y.c.ρq_sno[colidx] / Y.c.ρ[colidx] vertical_transport!(