feat: add face density method

Author: haakon-e

src/prognostic_equations/advection.jl

@@ -222,7 +222,10 @@ NVTX.@annotate function explicit_vertical_advection_tendency!(Yₜ, Y, p, t)
     advect_tke = use_prognostic_tke(turbconv_model)
     point_type = eltype(Fields.coordinate_field(Y.c))
     (; dt) = p
-    ᶜJ = Fields.local_geometry_field(Y.c).J
+    ᶜJ, ᶠJ = Fields.local_geometry_field(Y.c).J, Fields.local_geometry_field(Y.f).J
+    ᶜρ, ᶠρ = Y.c.ρ, ᶠface_density(Y.c.ρ)
+    ᶠρJ = @. lazy(ᶠρ * ᶠJ)
+    ᶜρJ = @. lazy(ᶜρ * ᶜJ)
     (; ᶜf³, ᶠf¹², ᶜΦ) = p.core
     (; ᶜu, ᶠu³, ᶜK, ᶜts) = p.precomputed
     (; edmfx_mse_q_tot_upwinding) = n > 0 || advect_tke ? p.atmos.numerics : all_nothing
@@ -241,14 +244,14 @@ NVTX.@annotate function explicit_vertical_advection_tendency!(Yₜ, Y, p, t)
         advect_tke ?
         (
             turbconv_model isa PrognosticEDMFX ?
-            (@. lazy(ρa⁰(Y.c.ρ, Y.c.sgsʲs, turbconv_model))) : Y.c.ρ
+            (@. lazy(ρa⁰(ᶜρ, Y.c.sgsʲs, turbconv_model))) : ᶜρ
         ) : nothing
     ᶜρ⁰ = if advect_tke
         if n > 0
             (; ᶜts⁰) = p.precomputed
             @. lazy(TD.air_density(thermo_params, ᶜts⁰))
         else
-            Y.c.ρ
+            ᶜρ
         end
     else
         nothing
@@ -278,8 +281,6 @@ NVTX.@annotate function explicit_vertical_advection_tendency!(Yₜ, Y, p, t)
     end
     # Without the CT12(), the right-hand side would be a CT1 or CT2 in 2D space.
 
-    ᶜρ = Y.c.ρ
-
     # Full vertical advection of passive tracers (like liq, rai, etc) ...
     # If sgs_mass_flux is true, the advection term is computed from the sum of SGS fluxes
     if !(
@@ -288,7 +289,7 @@ NVTX.@annotate function explicit_vertical_advection_tendency!(Yₜ, Y, p, t)
     )
         foreach_gs_tracer(Yₜ, Y) do ᶜρχₜ, ᶜρχ, ρχ_name
             if !(ρχ_name in (@name(ρe_tot), @name(ρq_tot)))
-                ᶜχ = @. lazy(specific(ᶜρχ, Y.c.ρ))
+                ᶜχ = @. lazy(specific(ᶜρχ, ᶜρ))
                 vtt = vertical_transport(ᶜρ, ᶠu³, ᶜχ, dt, tracer_upwinding)
                 @. ᶜρχₜ += vtt
             end
@@ -310,7 +311,7 @@ NVTX.@annotate function explicit_vertical_advection_tendency!(Yₜ, Y, p, t)
     end
 
     if !(p.atmos.moisture_model isa DryModel) && energy_q_tot_upwinding != Val(:none)
-        ᶜq_tot = @. lazy(specific(Y.c.ρq_tot, Y.c.ρ))
+        ᶜq_tot = @. lazy(specific(Y.c.ρq_tot, ᶜρ))
         vtt = vertical_transport(ᶜρ, ᶠu³, ᶜq_tot, dt, energy_q_tot_upwinding)
         vtt_central = vertical_transport(ᶜρ, ᶠu³, ᶜq_tot, dt, Val(:none))
         @. Yₜ.c.ρq_tot += vtt - vtt_central
@@ -322,9 +323,7 @@ NVTX.@annotate function explicit_vertical_advection_tendency!(Yₜ, Y, p, t)
 
     if isnothing(ᶠf¹²)
         # shallow atmosphere
-        @. Yₜ.c.uₕ -=
-            ᶜinterp(ᶠω¹² × (ᶠinterp(Y.c.ρ * ᶜJ) * ᶠu³)) / (Y.c.ρ * ᶜJ) +
-            (ᶜf³ + ᶜω³) × CT12(ᶜu)
+        @. Yₜ.c.uₕ -= ᶜinterp(ᶠω¹² × (ᶠρJ * ᶠu³)) / ᶜρJ + (ᶜf³ + ᶜω³) × CT12(ᶜu)
         @. Yₜ.f.u₃ -= ᶠω¹² × ᶠinterp(CT12(ᶜu)) + ᶠgradᵥ(ᶜK)
         for j in 1:n
             @. Yₜ.f.sgsʲs.:($$j).u₃ -=
@@ -333,9 +332,7 @@ NVTX.@annotate function explicit_vertical_advection_tendency!(Yₜ, Y, p, t)
         end
     else
         # deep atmosphere
-        @. Yₜ.c.uₕ -=
-            ᶜinterp((ᶠf¹² + ᶠω¹²) × (ᶠinterp(Y.c.ρ * ᶜJ) * ᶠu³)) /
-            (Y.c.ρ * ᶜJ) + (ᶜf³ + ᶜω³) × CT12(ᶜu)
+        @. Yₜ.c.uₕ -= ᶜinterp((ᶠf¹² + ᶠω¹²) × (ᶠρJ * ᶠu³)) / ᶜρJ + (ᶜf³ + ᶜω³) × CT12(ᶜu)
         @. Yₜ.f.u₃ -= (ᶠf¹² + ᶠω¹²) × ᶠinterp(CT12(ᶜu)) + ᶠgradᵥ(ᶜK)
         for j in 1:n
             @. Yₜ.f.sgsʲs.:($$j).u₃ -=
@@ -380,11 +377,7 @@ Modifies EDMF updraft fields in `Yₜ.c.sgsʲs` and `Yₜ.f.sgsʲs`.
 """
 edmfx_sgs_vertical_advection_tendency!(Yₜ, Y, p, t, turbconv_model) = nothing
 
-function edmfx_sgs_vertical_advection_tendency!(
-    Yₜ,
-    Y,
-    p,
-    t,
+function edmfx_sgs_vertical_advection_tendency!(Yₜ, Y, p, t,
     turbconv_model::PrognosticEDMFX,
 )
     (; params) = p
@@ -426,27 +419,21 @@ function edmfx_sgs_vertical_advection_tendency!(
     end
 
     for j in 1:n
-        ᶜa = (@. lazy(draft_area(Y.c.sgsʲs.:($$j).ρa, ᶜρʲs.:($$j))))
+        ᶜρʲ = ᶜρʲs.:($j)
+        ᶜa = (@. lazy(draft_area(Y.c.sgsʲs.:($$j).ρa, ᶜρʲ)))
+        ᶠu³ʲ = ᶠu³ʲs.:($j)
+        ᶜsgsʲ = Y.c.sgsʲs.:($j)
 
         # Flux form vertical advection of area farction with the grid mean velocity
-        vtt =
-            vertical_transport(ᶜρʲs.:($j), ᶠu³ʲs.:($j), ᶜa, dt, edmfx_mse_q_tot_upwinding)
+        vtt = vertical_transport(ᶜρʲ, ᶠu³ʲ, ᶜa, dt, edmfx_mse_q_tot_upwinding)
         @. Yₜ.c.sgsʲs.:($$j).ρa += vtt
 
         # Advective form advection of mse and q_tot with the grid mean velocity
         # Note: This allocates because the function is too long
-        va = vertical_advection(
-            ᶠu³ʲs.:($j),
-            Y.c.sgsʲs.:($j).mse,
-            edmfx_mse_q_tot_upwinding,
-        )
+        va = vertical_advection(ᶠu³ʲ, ᶜsgsʲ.mse, edmfx_mse_q_tot_upwinding)
         @. Yₜ.c.sgsʲs.:($$j).mse += va
 
-        va = vertical_advection(
-            ᶠu³ʲs.:($j),
-            Y.c.sgsʲs.:($j).q_tot,
-            edmfx_mse_q_tot_upwinding,
-        )
+        va = vertical_advection(ᶠu³ʲ, ᶜsgsʲ.q_tot, edmfx_mse_q_tot_upwinding)
         @. Yₜ.c.sgsʲs.:($$j).q_tot += va
 
         if p.atmos.moisture_model isa NonEquilMoistModel && (
@@ -464,15 +451,8 @@ function edmfx_sgs_vertical_advection_tendency!(
             ᶜ∂ρ∂t_sed = p.scratch.ᶜtemp_scalar_3
             @. ᶜ∂ρ∂t_sed = 0
 
-            ᶜinv_ρ̂ = (@. lazy(
-                specific(
-                    FT(1),
-                    Y.c.sgsʲs.:($$j).ρa,
-                    FT(0),
-                    ᶜρʲs.:($$j),
-                    turbconv_model,
-                ),
-            ))
+            ᶜinv_ρ̂ =
+                @. lazy(specific(FT(1), Y.c.sgsʲs.:($$j).ρa, FT(0), ᶜρʲ, turbconv_model))
 
             # Sedimentation
             # TODO - lazify ᶜwₗʲs computation. No need to cache it.
@@ -493,21 +473,11 @@ function edmfx_sgs_vertical_advection_tendency!(
                 ᶜwʲ = MatrixFields.get_field(p.precomputed, wʲ_name)
 
                 # Advective form advection of tracers with updraft velocity
-                va = vertical_advection(
-                    ᶠu³ʲs.:($j),
-                    ᶜqʲ,
-                    edmfx_tracer_upwinding,
-                )
+                va = vertical_advection(ᶠu³ʲs.:($j), ᶜqʲ, edmfx_tracer_upwinding)
                 @. ᶜqʲₜ += va
 
                 # Flux form sedimentation of tracers
-                vtt = updraft_sedimentation(
-                    ᶜρʲs.:($j),
-                    ᶜwʲ,
-                    ᶜa,
-                    ᶜqʲ,
-                    ᶠJ,
-                )
+                vtt = updraft_sedimentation(ᶜρʲs.:($j), ᶜwʲ, ᶜa, ᶜqʲ, ᶠJ)
                 @. ᶜqʲₜ += ᶜinv_ρ̂ * vtt
                 @. Yₜ.c.sgsʲs.:($$j).q_tot += ᶜinv_ρ̂ * vtt
                 @. ᶜ∂ρ∂t_sed += vtt
@@ -524,15 +494,8 @@ function edmfx_sgs_vertical_advection_tendency!(
                 else
                     error("Unsupported moisture tracer variable")
                 end
-                vtt = updraft_sedimentation(
-                    ᶜρʲs.:($j),
-                    ᶜwʲ,
-                    ᶜa,
-                    ᶜqʲ .* ᶜmse_li,
-                    ᶠJ,
-                )
-                @. Yₜ.c.sgsʲs.:($$j).mse +=
-                    ᶜinv_ρ̂ * vtt
+                vtt = updraft_sedimentation(ᶜρʲ, ᶜwʲ, ᶜa, ᶜqʲ .* ᶜmse_li, ᶠJ)
+                @. Yₜ.c.sgsʲs.:($$j).mse += ᶜinv_ρ̂ * vtt
             end
 
             # Contribution of density variation due to sedimentation
@@ -573,21 +536,11 @@ function edmfx_sgs_vertical_advection_tendency!(
                 ᶜwʲ = MatrixFields.get_field(p.precomputed, wʲ_name)
 
                 # Advective form advection of tracers with updraft velocity
-                va = vertical_transport(
-                    ᶠu³ʲs.:($j),
-                    ᶜχʲ,
-                    edmfx_tracer_upwinding,
-                )
+                va = vertical_transport(ᶠu³ʲs.:($j), ᶜχʲ, edmfx_tracer_upwinding)
                 @. ᶜχʲₜ += va
 
                 # Flux form sedimentation of tracers
-                vtt = updraft_sedimentation(
-                    ᶜρʲs.:($j),
-                    ᶜwʲ,
-                    ᶜa,
-                    ᶜχʲ,
-                    ᶠJ,
-                )
+                vtt = updraft_sedimentation(ᶜρʲ, ᶜwʲ, ᶜa, ᶜχʲ, ᶠJ)
                 @. ᶜχʲₜ += ᶜinv_ρ̂ * vtt
 
                 # Contribution of density variation due to sedimentation

src/prognostic_equations/implicit/implicit_tendency.jl

@@ -73,28 +73,22 @@ end
 # the implicit tendency function. Since dt >= dtγ, we can safely use dt for now.
 
 function vertical_transport(ᶜρ, ᶠu³, ᶜχ, dt, ::Val{:none})
-    ᶜJ = Fields.local_geometry_field(axes(ᶜρ)).J
-    ᶠJ = Fields.local_geometry_field(axes(ᶠu³)).J
-    return @. lazy(-(ᶜadvdivᵥ(ᶠinterp(ᶜρ * ᶜJ) / ᶠJ * ᶠu³ * ᶠinterp(ᶜχ))))
+    ᶠρ = ᶠface_density(ᶜρ)
+    return @. lazy(-(ᶜadvdivᵥ(ᶠρ * ᶠu³ * ᶠinterp(ᶜχ))))
 end
 function vertical_transport(ᶜρ, ᶠu³, ᶜχ, dt, ::Val{:first_order})
-    ᶜJ = Fields.local_geometry_field(axes(ᶜρ)).J
-    ᶠJ = Fields.local_geometry_field(axes(ᶠu³)).J
-    return @. lazy(-(ᶜadvdivᵥ(ᶠinterp(ᶜρ * ᶜJ) / ᶠJ * ᶠupwind1(ᶠu³, ᶜχ))))
+    ᶠρ = ᶠface_density(ᶜρ)
+    return @. lazy(-(ᶜadvdivᵥ(ᶠρ * ᶠupwind1(ᶠu³, ᶜχ))))
 end
 @static if pkgversion(ClimaCore) ≥ v"0.14.22"
     function vertical_transport(ᶜρ, ᶠu³, ᶜχ, dt, ::Val{:vanleer_limiter})
-        ᶜJ = Fields.local_geometry_field(axes(ᶜρ)).J
-        ᶠJ = Fields.local_geometry_field(axes(ᶠu³)).J
-        return @. lazy(
-            -(ᶜadvdivᵥ(ᶠinterp(ᶜρ * ᶜJ) / ᶠJ * ᶠlin_vanleer(ᶠu³, ᶜχ, dt))),
-        )
+        ᶠρ = ᶠface_density(ᶜρ)
+        return @. lazy(-(ᶜadvdivᵥ(ᶠρ * ᶠlin_vanleer(ᶠu³, ᶜχ, dt))))
     end
 end
 function vertical_transport(ᶜρ, ᶠu³, ᶜχ, dt, ::Val{:third_order})
-    ᶜJ = Fields.local_geometry_field(axes(ᶜρ)).J
-    ᶠJ = Fields.local_geometry_field(axes(ᶠu³)).J
-    return @. lazy(-(ᶜadvdivᵥ(ᶠinterp(ᶜρ * ᶜJ) / ᶠJ * ᶠupwind3(ᶠu³, ᶜχ))))
+    ᶠρ = ᶠface_density(ᶜρ)
+    return @. lazy(-(ᶜadvdivᵥ(ᶠρ * ᶠupwind3(ᶠu³, ᶜχ))))
 end
 
 vertical_advection(ᶠu³, ᶜχ, ::Val{:none}) =
@@ -105,32 +99,25 @@ vertical_advection(ᶠu³, ᶜχ, ::Val{:third_order}) =
     @. lazy(-(ᶜadvdivᵥ(ᶠupwind3(ᶠu³, ᶜχ)) - ᶜχ * ᶜadvdivᵥ(ᶠu³)))
 
 function implicit_vertical_advection_tendency!(Yₜ, Y, p, t)
-    (; moisture_model, turbconv_model, rayleigh_sponge, microphysics_model) =
-        p.atmos
+    (; moisture_model, turbconv_model, rayleigh_sponge, microphysics_model) = p.atmos
     (; params, dt) = p
     n = n_mass_flux_subdomains(turbconv_model)
-    ᶜJ = Fields.local_geometry_field(axes(Y.c)).J
-    ᶠJ = Fields.local_geometry_field(axes(Y.f)).J
+    ᶜρ, ᶠρ = Y.c.ρ, ᶠface_density(Y.c.ρ)
     (; ᶠgradᵥ_ᶜΦ) = p.core
-    (; ᶠu³, ᶜp, ᶜts) = p.precomputed
+    (; ᶠu³, ᶜts) = p.precomputed
     thermo_params = CAP.thermodynamics_params(params)
     cp_d = CAP.cp_d(params)
-    ᶜh_tot = @. lazy(
-        TD.total_specific_enthalpy(
-            thermo_params,
-            ᶜts,
-            specific(Y.c.ρe_tot, Y.c.ρ),
-        ),
-    )
+    ᶜe_tot = @. lazy(specific(Y.c.ρe_tot, ᶜρ))
+    ᶜh_tot = @. lazy(TD.total_specific_enthalpy(thermo_params, ᶜts, ᶜe_tot))
 
-    @. Yₜ.c.ρ -= ᶜdivᵥ(ᶠinterp(Y.c.ρ * ᶜJ) / ᶠJ * ᶠu³)
+    @. Yₜ.c.ρ -= ᶜdivᵥ(ᶠρ * ᶠu³)
 
     # Central vertical advection of active tracers (e_tot and q_tot)
-    vtt = vertical_transport(Y.c.ρ, ᶠu³, ᶜh_tot, dt, Val(:none))
+    vtt = vertical_transport(ᶜρ, ᶠu³, ᶜh_tot, dt, Val(:none))
     @. Yₜ.c.ρe_tot += vtt
     if !(moisture_model isa DryModel)
-        ᶜq_tot = @. lazy(specific(Y.c.ρq_tot, Y.c.ρ))
-        vtt = vertical_transport(Y.c.ρ, ᶠu³, ᶜq_tot, dt, Val(:none))
+        ᶜq_tot = @. lazy(specific(Y.c.ρq_tot, ᶜρ))
+        vtt = vertical_transport(ᶜρ, ᶠu³, ᶜq_tot, dt, Val(:none))
         @. Yₜ.c.ρq_tot += vtt
     end
 
@@ -140,63 +127,39 @@ function implicit_vertical_advection_tendency!(Yₜ, Y, p, t)
     # using downward biasing and free outflow bottom boundary condition
     if moisture_model isa NonEquilMoistModel
         (; ᶜwₗ, ᶜwᵢ) = p.precomputed
-        @. Yₜ.c.ρq_liq -= ᶜprecipdivᵥ(
-            ᶠinterp(Y.c.ρ * ᶜJ) / ᶠJ * ᶠright_bias(
-                Geometry.WVector(-(ᶜwₗ)) * specific(Y.c.ρq_liq, Y.c.ρ),
-            ),
-        )
-        @. Yₜ.c.ρq_ice -= ᶜprecipdivᵥ(
-            ᶠinterp(Y.c.ρ * ᶜJ) / ᶠJ * ᶠright_bias(
-                Geometry.WVector(-(ᶜwᵢ)) * specific(Y.c.ρq_ice, Y.c.ρ),
-            ),
-        )
+        q_liq = @. lazy(specific(Y.c.ρq_liq, ᶜρ))
+        q_ice = @. lazy(specific(Y.c.ρq_ice, ᶜρ))
+        @. Yₜ.c.ρq_liq -= ᶜprecipdivᵥ(ᶠρ * ᶠright_bias(-WVec(ᶜwₗ) * q_liq))
+        @. Yₜ.c.ρq_ice -= ᶜprecipdivᵥ(ᶠρ * ᶠright_bias(-WVec(ᶜwᵢ) * q_ice))
     end
     if microphysics_model isa Microphysics1Moment
         (; ᶜwᵣ, ᶜwₛ) = p.precomputed
-        @. Yₜ.c.ρq_rai -= ᶜprecipdivᵥ(
-            ᶠinterp(Y.c.ρ * ᶜJ) / ᶠJ * ᶠright_bias(
-                Geometry.WVector(-(ᶜwᵣ)) * specific(Y.c.ρq_rai, Y.c.ρ),
-            ),
-        )
-        @. Yₜ.c.ρq_sno -= ᶜprecipdivᵥ(
-            ᶠinterp(Y.c.ρ * ᶜJ) / ᶠJ * ᶠright_bias(
-                Geometry.WVector(-(ᶜwₛ)) * specific(Y.c.ρq_sno, Y.c.ρ),
-            ),
-        )
+        q_rai = @. lazy(specific(Y.c.ρq_rai, ᶜρ))
+        q_sno = @. lazy(specific(Y.c.ρq_sno, ᶜρ))
+        @. Yₜ.c.ρq_rai -= ᶜprecipdivᵥ(ᶠρ * ᶠright_bias(-WVec(ᶜwᵣ) * q_rai))
+        @. Yₜ.c.ρq_sno -= ᶜprecipdivᵥ(ᶠρ * ᶠright_bias(-WVec(ᶜwₛ) * q_sno))
     end
     if microphysics_model isa Microphysics2Moment
         (; ᶜwₙₗ, ᶜwₙᵣ, ᶜwᵣ, ᶜwₛ) = p.precomputed
-        @. Yₜ.c.ρn_liq -= ᶜprecipdivᵥ(
-            ᶠinterp(Y.c.ρ * ᶜJ) / ᶠJ * ᶠright_bias(
-                Geometry.WVector(-(ᶜwₙₗ)) * specific(Y.c.ρn_liq, Y.c.ρ),
-            ),
-        )
-        @. Yₜ.c.ρn_rai -= ᶜprecipdivᵥ(
-            ᶠinterp(Y.c.ρ * ᶜJ) / ᶠJ * ᶠright_bias(
-                Geometry.WVector(-(ᶜwₙᵣ)) * specific(Y.c.ρn_rai, Y.c.ρ),
-            ),
-        )
-        @. Yₜ.c.ρq_rai -= ᶜprecipdivᵥ(
-            ᶠinterp(Y.c.ρ * ᶜJ) / ᶠJ * ᶠright_bias(
-                Geometry.WVector(-(ᶜwᵣ)) * specific(Y.c.ρq_rai, Y.c.ρ),
-            ),
-        )
-        @. Yₜ.c.ρq_sno -= ᶜprecipdivᵥ(
-            ᶠinterp(Y.c.ρ * ᶜJ) / ᶠJ * ᶠright_bias(
-                Geometry.WVector(-(ᶜwₛ)) * specific(Y.c.ρq_sno, Y.c.ρ),
-            ),
-        )
+        n_liq = @. lazy(specific(Y.c.ρn_liq, ᶜρ))
+        n_rai = @. lazy(specific(Y.c.ρn_rai, ᶜρ))
+        q_rai = @. lazy(specific(Y.c.ρq_rai, ᶜρ))
+        q_sno = @. lazy(specific(Y.c.ρq_sno, ᶜρ))
+        @. Yₜ.c.ρn_liq -= ᶜprecipdivᵥ(ᶠρ * ᶠright_bias(-WVec(ᶜwₙₗ) * n_liq))
+        @. Yₜ.c.ρn_rai -= ᶜprecipdivᵥ(ᶠρ * ᶠright_bias(-WVec(ᶜwₙᵣ) * n_rai))
+        @. Yₜ.c.ρq_rai -= ᶜprecipdivᵥ(ᶠρ * ᶠright_bias(-WVec(ᶜwᵣ) * q_rai))
+        @. Yₜ.c.ρq_sno -= ᶜprecipdivᵥ(ᶠρ * ᶠright_bias(-WVec(ᶜwₛ) * q_sno))
     end
     if microphysics_model isa Microphysics2MomentP3
-        (; ρ, ρn_ice, ρq_rim, ρb_rim) = Y.c
-        ᶜwnᵢ = @. lazy(Geometry.WVector(p.precomputed.ᶜwnᵢ))
-        ᶜwᵢ = @. lazy(Geometry.WVector(p.precomputed.ᶜwᵢ))
-        ᶠρ = @. lazy(ᶠinterp(ρ * ᶜJ) / ᶠJ)
+        (; ᶜwnᵢ, ᶜwᵢ) = p.precomputed
 
         # Note: `ρq_ice` is handled above, in `moisture_model isa NonEquilMoistModel`
-        @. Yₜ.c.ρn_ice -= ᶜprecipdivᵥ(ᶠρ * ᶠright_bias(- ᶜwnᵢ * specific(ρn_ice, ρ)))
-        @. Yₜ.c.ρq_rim -= ᶜprecipdivᵥ(ᶠρ * ᶠright_bias(- ᶜwᵢ * specific(ρq_rim, ρ)))
-        @. Yₜ.c.ρb_rim -= ᶜprecipdivᵥ(ᶠρ * ᶠright_bias(- ᶜwᵢ * specific(ρb_rim, ρ)))
+        n_ice = @. lazy(specific(Y.c.ρn_ice, ᶜρ))
+        q_rim = @. lazy(specific(Y.c.ρq_rim, ᶜρ))
+        b_rim = @. lazy(specific(Y.c.ρb_rim, ᶜρ))
+        @. Yₜ.c.ρn_ice -= ᶜprecipdivᵥ(ᶠρ * ᶠright_bias(-WVec(ᶜwnᵢ) * n_ice))
+        @. Yₜ.c.ρq_rim -= ᶜprecipdivᵥ(ᶠρ * ᶠright_bias(-WVec(ᶜwᵢ) * q_rim))
+        @. Yₜ.c.ρb_rim -= ᶜprecipdivᵥ(ᶠρ * ᶠright_bias(-WVec(ᶜwᵢ) * b_rim))
     end
 
     # TODO - decide if this needs to be explicit or implicit
@@ -207,8 +170,7 @@ function implicit_vertical_advection_tendency!(Yₜ, Y, p, t)
     ᶜθ_v = @. lazy(theta_v(thermo_params, ᶜts))
     ᶜθ_vr = @. lazy(theta_vr(thermo_params, ᶜts))
     ᶜΠ = @. lazy(dry_exner_function(thermo_params, ᶜts))
-    @. Yₜ.f.u₃ -= ᶠgradᵥ_ᶜΦ - ᶠgradᵥ(ᶜΦ_r) +
-                  cp_d * (ᶠinterp(ᶜθ_v - ᶜθ_vr)) * ᶠgradᵥ(ᶜΠ)
+    @. Yₜ.f.u₃ -= ᶠgradᵥ_ᶜΦ - ᶠgradᵥ(ᶜΦ_r) + cp_d * (ᶠinterp(ᶜθ_v - ᶜθ_vr)) * ᶠgradᵥ(ᶜΠ)
 
     if rayleigh_sponge isa RayleighSponge
         ᶠz = Fields.coordinate_field(Y.f).z
@@ -217,8 +179,7 @@ function implicit_vertical_advection_tendency!(Yₜ, Y, p, t)
         @. Yₜ.f.u₃ -= β_rayleigh_w(rs, ᶠz, zmax) * Y.f.u₃
         if turbconv_model isa PrognosticEDMFX
             for j in 1:n
-                @. Yₜ.f.sgsʲs.:($$j).u₃ -=
-                    β_rayleigh_w(rs, ᶠz, zmax) * Y.f.sgsʲs.:($$j).u₃
+                @. Yₜ.f.sgsʲs.:($$j).u₃ -= β_rayleigh_w(rs, ᶠz, zmax) * Y.f.sgsʲs.:($$j).u₃
             end
         end
     end