diff --git a/.buildkite/Manifest-v1.11.toml b/.buildkite/Manifest-v1.11.toml
index 891d84bf96..f0aabe571c 100644
--- a/.buildkite/Manifest-v1.11.toml
+++ b/.buildkite/Manifest-v1.11.toml
@@ -381,9 +381,9 @@ weakdeps = ["CUDA", "MPI"]
 
 [[deps.ClimaCore]]
 deps = ["Adapt", "BandedMatrices", "BlockArrays", "ClimaComms", "CubedSphere", "DataStructures", "ForwardDiff", "GaussQuadrature", "GilbertCurves", "HDF5", "InteractiveUtils", "IntervalSets", "KrylovKit", "LazyBroadcast", "LinearAlgebra", "MultiBroadcastFusion", "NVTX", "PkgVersion", "RecursiveArrayTools", "RootSolvers", "SparseArrays", "StaticArrays", "Statistics", "UnrolledUtilities"]
-git-tree-sha1 = "14d3d5810ce1e3c990450a2ce7abc6a1e162855f"
+git-tree-sha1 = "a273452127dbb052f2963e3c1095730a996d49a6"
 uuid = "d414da3d-4745-48bb-8d80-42e94e092884"
-version = "0.14.35"
+version = "0.14.36"
 weakdeps = ["CUDA", "Krylov"]
 
     [deps.ClimaCore.extensions]
diff --git a/.buildkite/pipeline.yml b/.buildkite/pipeline.yml
index 274e3767b0..d8eb0f63e0 100644
--- a/.buildkite/pipeline.yml
+++ b/.buildkite/pipeline.yml
@@ -83,7 +83,7 @@ steps:
           --config_file $CONFIG_PATH/single_column_precipitation_test.yml
           --job_id single_column_precipitation_test
         artifact_paths: "single_column_precipitation_test/output_active/*"
-      
+
       - label: ":umbrella: 2-moment precipitation sanity test single column"
         command: >
           julia --color=yes --project=.buildkite .buildkite/ci_driver.jl
@@ -854,6 +854,15 @@ steps:
         agents:
           slurm_mem: 20GB
 
+      - label: ":genie: Prognostic EDMFX Dycoms RF02 in a column"
+        command: >
+          julia --color=yes --project=.buildkite .buildkite/ci_driver.jl
+          --config_file $CONFIG_PATH/prognostic_edmfx_dycoms_rf02_column.yml
+          --job_id prognostic_edmfx_dycoms_rf02_column
+        artifact_paths: "prognostic_edmfx_dycoms_rf02_column/output_active/*"
+        agents:
+          slurm_mem: 20GB
+
       - label: ":umbrella: Prognostic EDMFX Rico in a column"
         command: >
           julia --color=yes --project=.buildkite .buildkite/ci_driver.jl
@@ -1035,9 +1044,9 @@ steps:
         agents:
           slurm_gpus: 1
           slurm_mem: 16GB
-      
+
       - label: ":earth_americas: GPU: Diagnostic EDMFX ERA5 Weather Model Initial Condition"
-        command: > 
+        command: >
           julia --color=yes --project=.buildkite .buildkite/ci_driver.jl
           --config_file $CONFIG_PATH/diagnostic_edmfx_era5_initial_condition.yml
           --job_id diagnostic_edmfx_era5_initial_condition
diff --git a/Project.toml b/Project.toml
index bdefca8f3f..c2bd1f065d 100644
--- a/Project.toml
+++ b/Project.toml
@@ -43,7 +43,7 @@ ArgParse = "1"
 Artifacts = "1"
 AtmosphericProfilesLibrary = "0.1.7"
 ClimaComms = "0.6.8"
-ClimaCore = "0.14.35"
+ClimaCore = "0.14.36"
 ClimaDiagnostics = "0.2.12"
 ClimaInterpolations = "0.1.0"
 ClimaParams = "0.10.35"
diff --git a/config/model_configs/prognostic_edmfx_dycoms_rf02_column.yml b/config/model_configs/prognostic_edmfx_dycoms_rf02_column.yml
new file mode 100644
index 0000000000..e14c8fdcfa
--- /dev/null
+++ b/config/model_configs/prognostic_edmfx_dycoms_rf02_column.yml
@@ -0,0 +1,43 @@
+initial_condition: DYCOMS_RF02
+subsidence: DYCOMS
+scm_coriolis: DYCOMS_RF02
+rad: DYCOMS
+surface_setup: DYCOMS_RF02
+turbconv: "prognostic_edmfx"
+implicit_diffusion: false
+approximate_linear_solve_iters: 2
+edmfx_upwinding: first_order
+edmfx_entr_model: "Generalized"
+edmfx_detr_model: "Generalized"
+edmfx_sgs_mass_flux: true
+edmfx_sgs_diffusive_flux: true
+edmfx_nh_pressure: true
+edmfx_filter: true
+prognostic_tke: true
+moist: "nonequil"
+cloud_model: "quadrature_sgs"
+precip_model: "1M"
+call_cloud_diagnostics_per_stage: true
+config: "column"
+x_elem: 2
+y_elem: 2
+z_elem: 30
+z_max: 1500
+z_stretch: false
+perturb_initstate: false
+dt: 10secs
+t_end: 6hours
+dt_save_state_to_disk: 10mins
+toml: [toml/prognostic_edmfx_1M.toml]
+netcdf_interpolation_num_points: [8, 8, 30]
+diagnostics:
+  - short_name: [ts, ta, thetaa, ha, pfull, rhoa, ua, va, wa, hur, hus, cl, clw, cli, hussfc, evspsbl, pr]
+    period: 10mins
+  - short_name: [arup, waup, taup, thetaaup, haup, husup, hurup, clwup, cliup, husraup, hussnup]
+    period: 10mins
+  - short_name: [waen, taen, thetaaen, haen, husen, huren, clwen, clien, husraen, hussnen, tke]
+    period: 10mins
+  - short_name: [entr, detr, lmix, bgrad, strain, edt, evu]
+    period: 10mins
+  - short_name: [husra, hussn]
+    period: 10mins
diff --git a/post_processing/ci_plots.jl b/post_processing/ci_plots.jl
index d5c346e0d5..69ce28b990 100644
--- a/post_processing/ci_plots.jl
+++ b/post_processing/ci_plots.jl
@@ -1400,6 +1400,7 @@ EDMFBoxPlots = Union{
 EDMFBoxPlotsWithPrecip = Union{
     Val{:prognostic_edmfx_rico_column},
     Val{:prognostic_edmfx_trmm_column},
+    Val{:prognostic_edmfx_dycoms_rf02_column},
 }
 """
     plot_edmf_vert_profile!(grid_loc, var_group)
diff --git a/src/cache/precomputed_quantities.jl b/src/cache/precomputed_quantities.jl
index 04d6e3bca6..33d77ccf32 100644
--- a/src/cache/precomputed_quantities.jl
+++ b/src/cache/precomputed_quantities.jl
@@ -165,6 +165,12 @@ function precomputed_quantities(Y, atmos)
             ᶜSqᵢᵖʲs = similar(Y.c, NTuple{n, FT}),
             ᶜSqᵣᵖʲs = similar(Y.c, NTuple{n, FT}),
             ᶜSqₛᵖʲs = similar(Y.c, NTuple{n, FT}),
+            ᶜwₗʲs = similar(Y.c, NTuple{n, FT}),
+            ᶜwᵢʲs = similar(Y.c, NTuple{n, FT}),
+            ᶜwᵣʲs = similar(Y.c, NTuple{n, FT}),
+            ᶜwₛʲs = similar(Y.c, NTuple{n, FT}),
+            ᶜwₜʲs = similar(Y.c, NTuple{n, FT}),
+            ᶜwₕʲs = similar(Y.c, NTuple{n, FT}),
             ᶜSqₗᵖ⁰ = similar(Y.c, FT),
             ᶜSqᵢᵖ⁰ = similar(Y.c, FT),
             ᶜSqᵣᵖ⁰ = similar(Y.c, FT),
diff --git a/src/cache/prognostic_edmf_precomputed_quantities.jl b/src/cache/prognostic_edmf_precomputed_quantities.jl
index 131a780ccd..bc9262c52d 100644
--- a/src/cache/prognostic_edmf_precomputed_quantities.jl
+++ b/src/cache/prognostic_edmf_precomputed_quantities.jl
@@ -608,13 +608,73 @@ NVTX.@annotate function set_prognostic_edmf_precomputed_quantities_precipitation
     (; ᶜSqₗᵖ⁰, ᶜSqᵢᵖ⁰, ᶜSqᵣᵖ⁰, ᶜSqₛᵖ⁰, ᶜρ⁰, ᶜts⁰) = p.precomputed
     (; ᶜq_tot⁰, ᶜq_liq⁰, ᶜq_ice⁰, ᶜq_rai⁰, ᶜq_sno⁰) = p.precomputed
 
+    (; ᶜwₗʲs, ᶜwᵢʲs, ᶜwᵣʲs, ᶜwₛʲs, ᶜwₜʲs, ᶜwₕʲs) = p.precomputed
+
     # TODO - can I re-use them between js and env?
     ᶜSᵖ = p.scratch.ᶜtemp_scalar
     ᶜSᵖ_snow = p.scratch.ᶜtemp_scalar_2
 
     n = n_mass_flux_subdomains(p.atmos.turbconv_model)
+    FT = eltype(params)
 
     for j in 1:n
+
+        # compute terminal velocity for precipitation
+        @. ᶜwᵣʲs.:($$j) = CM1.terminal_velocity(
+            cmp.pr,
+            cmp.tv.rain,
+            ᶜρʲs.:($$j),
+            max(zero(Y.c.ρ), Y.c.sgsʲs.:($$j).q_rai),
+        )
+        @. ᶜwₛʲs.:($$j) = CM1.terminal_velocity(
+            cmp.ps,
+            cmp.tv.snow,
+            ᶜρʲs.:($$j),
+            max(zero(Y.c.ρ), Y.c.sgsʲs.:($$j).q_sno),
+        )
+        # compute sedimentation velocity for cloud condensate [m/s]
+        @. ᶜwₗʲs.:($$j) = CMNe.terminal_velocity(
+            cmc.liquid,
+            cmc.Ch2022.rain,
+            ᶜρʲs.:($$j),
+            max(zero(Y.c.ρ), Y.c.sgsʲs.:($$j).q_liq),
+        )
+        @. ᶜwᵢʲs.:($$j) = CMNe.terminal_velocity(
+            cmc.ice,
+            cmc.Ch2022.small_ice,
+            ᶜρʲs.:($$j),
+            max(zero(Y.c.ρ), Y.c.sgsʲs.:($$j).q_ice),
+        )
+        # compute their contirbutions to energy and total water advection
+        @. ᶜwₜʲs.:($$j) = ifelse(
+            Y.c.sgsʲs.:($$j).ρa * Y.c.sgsʲs.:($$j).q_tot > FT(0),
+            (
+                ᶜwₗʲs.:($$j) * Y.c.sgsʲs.:($$j).q_liq +
+                ᶜwᵢʲs.:($$j) * Y.c.sgsʲs.:($$j).q_ice +
+                ᶜwᵣʲs.:($$j) * Y.c.sgsʲs.:($$j).q_rai +
+                ᶜwₛʲs.:($$j) * Y.c.sgsʲs.:($$j).q_sno
+            ) / Y.c.sgsʲs.:($$j).ρa / Y.c.sgsʲs.:($$j).q_tot,
+            FT(0),
+        )
+        @. ᶜwₕʲs.:($$j) = ifelse(
+            Y.c.sgsʲs.:($$j).ρa * abs(Y.c.sgsʲs.:($$j).mse) > FT(0),
+            (
+                ᶜwₗʲs.:($$j) *
+                Y.c.sgsʲs.:($$j).q_liq *
+                (Iₗ(thp, ᶜtsʲs.:($$j)) + ᶜΦ) +
+                ᶜwᵢʲs.:($$j) *
+                Y.c.sgsʲs.:($$j).q_ice *
+                (Iᵢ(thp, ᶜtsʲs.:($$j)) + ᶜΦ) +
+                ᶜwᵣʲs.:($$j) *
+                Y.c.sgsʲs.:($$j).q_rai *
+                (Iₗ(thp, ᶜtsʲs.:($$j)) + ᶜΦ) +
+                ᶜwₛʲs.:($$j) *
+                Y.c.sgsʲs.:($$j).q_sno *
+                (Iᵢ(thp, ᶜtsʲs.:($$j)) + ᶜΦ)
+            ) / Y.c.sgsʲs.:($$j).ρa / abs(Y.c.sgsʲs.:($$j).mse),
+            FT(0),
+        )
+
         # Precipitation sources and sinks from the updrafts
         compute_precipitation_sources!(
             ᶜSᵖ,
diff --git a/src/prognostic_equations/advection.jl b/src/prognostic_equations/advection.jl
index fb31a95910..f64d3a7837 100644
--- a/src/prognostic_equations/advection.jl
+++ b/src/prognostic_equations/advection.jl
@@ -10,7 +10,7 @@ import ClimaCore.Geometry as Geometry
     horizontal_dynamics_tendency!(Yₜ, Y, p, t)
 
 Computes tendencies due to horizontal advection for prognostic variables of the
-grid mean and EDMFX subdomains, and also applies horizontal pressure gradient and 
+grid mean and EDMFX subdomains, and also applies horizontal pressure gradient and
 gravitational acceleration terms for horizontal momentum.
 
 Specifically, this function calculates:
@@ -339,6 +339,7 @@ function edmfx_sgs_vertical_advection_tendency!(
     end
 
     for j in 1:n
+        # Flux form vertical advection of area farction with the grid mean velocity
         @. ᶜa_scalar = draft_area(Y.c.sgsʲs.:($$j).ρa, ᶜρʲs.:($$j))
         vtt = vertical_transport(
             ᶜρʲs.:($j),
@@ -349,24 +350,28 @@ function edmfx_sgs_vertical_advection_tendency!(
         )
         @. Yₜ.c.sgsʲs.:($$j).ρa += vtt
 
+        # Advective form advection of mse and q_tot with the grid mean velocity
         va = vertical_advection(
             ᶠu³ʲs.:($j),
             Y.c.sgsʲs.:($j).mse,
             edmfx_upwinding,
         )
         @. Yₜ.c.sgsʲs.:($$j).mse += va
-
         va = vertical_advection(
             ᶠu³ʲs.:($j),
             Y.c.sgsʲs.:($j).q_tot,
             edmfx_upwinding,
         )
         @. Yₜ.c.sgsʲs.:($$j).q_tot += va
+
         if p.atmos.moisture_model isa NonEquilMoistModel &&
            p.atmos.microphysics_model isa Microphysics1Moment
-            # TODO - add precipitation terminal velocity
-            # TODO - add cloud sedimentation velocity
-            # TODO - add their contributions to mean energy and mass
+            # TODO - add precipitation terminal velocity in implicit solver/tendency with if/else
+            # TODO - add cloud sedimentation velocity in implicit solver/tendency with if/else
+
+            (; ᶜwₗʲs, ᶜwᵢʲs, ᶜwᵣʲs, ᶜwₛʲs, ᶜwₜʲs, ᶜwₕʲs) = p.precomputed
+
+            # Advective form advection of moisture tracers with the grid mean velocity
             va = vertical_advection(
                 ᶠu³ʲs.:($j),
                 Y.c.sgsʲs.:($j).q_liq,
@@ -391,6 +396,284 @@ function edmfx_sgs_vertical_advection_tendency!(
                 edmfx_upwinding,
             )
             @. Yₜ.c.sgsʲs.:($$j).q_sno += va
+
+            FT = eltype(params)
+            (; ᶜwₗʲs, ᶜwᵢʲs, ᶜwᵣʲs, ᶜwₛʲs, ᶜtsʲs) = p.precomputed
+
+            ᶠwₗ³ʲs =
+                (@. lazy(CT3(ᶠinterp(Geometry.WVector(-1 * ᶜwₗʲs.:($$j))))))
+            ᶠwᵢ³ʲs =
+                (@. lazy(CT3(ᶠinterp(Geometry.WVector(-1 * ᶜwᵢʲs.:($$j))))))
+            ᶠwᵣ³ʲs =
+                (@. lazy(CT3(ᶠinterp(Geometry.WVector(-1 * ᶜwᵣʲs.:($$j))))))
+            ᶠwₛ³ʲs =
+                (@. lazy(CT3(ᶠinterp(Geometry.WVector(-1 * ᶜwₛʲs.:($$j))))))
+
+            ᶜa = (@. lazy(draft_area(Y.c.sgsʲs.:($$j).ρa, ᶜρʲs.:($$j))))
+
+            # Flux form vertical advection of rho * area with sedimentation contributions
+            vtt = vertical_transport(
+                ᶜρʲs.:($j),
+                ᶠwₗ³ʲs,
+                (@. lazy(ᶜa * Y.c.sgsʲs.:($$j).q_liq)),
+                dt,
+                edmfx_upwinding,
+            )
+            @. Yₜ.c.sgsʲs.:($$j).ρa += vtt
+            vtt = vertical_transport(
+                ᶜρʲs.:($j),
+                ᶠwᵢ³ʲs,
+                (@. lazy(ᶜa * Y.c.sgsʲs.:($$j).q_ice)),
+                dt,
+                edmfx_upwinding,
+            )
+            @. Yₜ.c.sgsʲs.:($$j).ρa += vtt
+            vtt = vertical_transport(
+                ᶜρʲs.:($j),
+                ᶠwᵣ³ʲs,
+                (@. lazy(ᶜa * Y.c.sgsʲs.:($$j).q_rai)),
+                dt,
+                edmfx_upwinding,
+            )
+            @. Yₜ.c.sgsʲs.:($$j).ρa += vtt
+            vtt = vertical_transport(
+                ᶜρʲs.:($j),
+                ᶠwₛ³ʲs,
+                (@. lazy(ᶜa * Y.c.sgsʲs.:($$j).q_sno)),
+                dt,
+                edmfx_upwinding,
+            )
+            @. Yₜ.c.sgsʲs.:($$j).ρa += vtt
+
+            # q_tot and moisture tracers advective form advection with sedimentation velocity
+            va = vertical_advection(
+                ᶠwₗ³ʲs,
+                Y.c.sgsʲs.:($j).q_liq,
+                edmfx_upwinding,
+            )
+            @. Yₜ.c.sgsʲs.:($$j).q_tot += (1 - Y.c.sgsʲs.:($$j).q_tot) * va
+            @. Yₜ.c.sgsʲs.:($$j).q_liq += va
+            va = vertical_advection(
+                ᶠwᵢ³ʲs,
+                Y.c.sgsʲs.:($j).q_ice,
+                edmfx_upwinding,
+            )
+            @. Yₜ.c.sgsʲs.:($$j).q_tot += (1 - Y.c.sgsʲs.:($$j).q_tot) * va
+            @. Yₜ.c.sgsʲs.:($$j).q_ice += va
+            va = vertical_advection(
+                ᶠwᵣ³ʲs,
+                Y.c.sgsʲs.:($j).q_rai,
+                edmfx_upwinding,
+            )
+            @. Yₜ.c.sgsʲs.:($$j).q_tot += (1 - Y.c.sgsʲs.:($$j).q_tot) * va
+            @. Yₜ.c.sgsʲs.:($$j).q_rai += va
+            va = vertical_advection(
+                ᶠwₛ³ʲs,
+                Y.c.sgsʲs.:($j).q_sno,
+                edmfx_upwinding,
+            )
+            @. Yₜ.c.sgsʲs.:($$j).q_tot += (1 - Y.c.sgsʲs.:($$j).q_tot) * va
+            @. Yₜ.c.sgsʲs.:($$j).q_sno += va
+
+            # mse advective form advection with sedimentation velocity
+            (; ᶜΦ) = p.core
+            thp = CAP.thermodynamics_params(params)
+            ᶜmseₗ = (@. lazy(
+                Y.c.sgsʲs.:($$j).q_liq * (
+                    TD.internal_energy_liquid(thp, ᶜtsʲs.:($$j)) +
+                    TD.gas_constant_air(thp, ᶜtsʲs.:($$j)) *
+                    TD.air_temperature(thp, ᶜtsʲs.:($$j)) +
+                    ᶜΦ
+                ),
+            ))
+            ᶜmseᵢ = (@. lazy(
+                Y.c.sgsʲs.:($$j).q_ice * (
+                    TD.internal_energy_ice(thp, ᶜtsʲs.:($$j)) +
+                    TD.gas_constant_air(thp, ᶜtsʲs.:($$j)) *
+                    TD.air_temperature(thp, ᶜtsʲs.:($$j)) +
+                    ᶜΦ
+                ),
+            ))
+            ᶜmseᵣ = (@. lazy(
+                Y.c.sgsʲs.:($$j).q_rai * (
+                    TD.internal_energy_liquid(thp, ᶜtsʲs.:($$j)) +
+                    TD.gas_constant_air(thp, ᶜtsʲs.:($$j)) *
+                    TD.air_temperature(thp, ᶜtsʲs.:($$j)) +
+                    ᶜΦ
+                ),
+            ))
+            ᶜmseₛ = (@. lazy(
+                Y.c.sgsʲs.:($$j).q_sno * (
+                    TD.internal_energy_ice(thp, ᶜtsʲs.:($$j)) +
+                    TD.gas_constant_air(thp, ᶜtsʲs.:($$j)) *
+                    TD.air_temperature(thp, ᶜtsʲs.:($$j)) +
+                    ᶜΦ
+                ),
+            ))
+
+            # RICO works well but we are missing mse * q_ advection term
+            va = vertical_advection(ᶠwₗ³ʲs, ᶜmseₗ, edmfx_upwinding)
+            @. Yₜ.c.sgsʲs.:($$j).mse += va
+            va = vertical_advection(ᶠwᵢ³ʲs, ᶜmseᵢ, edmfx_upwinding)
+            @. Yₜ.c.sgsʲs.:($$j).mse += va
+            va = vertical_advection(ᶠwᵣ³ʲs, ᶜmseᵣ, edmfx_upwinding)
+            @. Yₜ.c.sgsʲs.:($$j).mse += va
+            va = vertical_advection(ᶠwₛ³ʲs, ᶜmseₛ, edmfx_upwinding)
+            @. Yₜ.c.sgsʲs.:($$j).mse += va
+            #=
+                        va = vertical_advection(ᶠwₗ³ʲs, Y.c.sgsʲs.:($j).q_liq, edmfx_upwinding)
+                        @. Yₜ.c.sgsʲs.:($$j).mse -= Y.c.sgsʲs.:($$j).mse * va
+                        va = vertical_advection(ᶠwᵢ³ʲs, Y.c.sgsʲs.:($j).q_ice, edmfx_upwinding)
+                        @. Yₜ.c.sgsʲs.:($$j).mse -= Y.c.sgsʲs.:($$j).mse * va
+                        va = vertical_advection(ᶠwᵣ³ʲs, Y.c.sgsʲs.:($j).q_rai, edmfx_upwinding)
+                        @. Yₜ.c.sgsʲs.:($$j).mse -= Y.c.sgsʲs.:($$j).mse * va
+                        va = vertical_advection(ᶠwₛ³ʲs, Y.c.sgsʲs.:($j).q_sno, edmfx_upwinding)
+                        @. Yₜ.c.sgsʲs.:($$j).mse -= Y.c.sgsʲs.:($$j).mse * va
+            =#
+
+            # mse, q_tot and moisture tracers terms proportional to 1/ρ̂ ∂zρ̂
+            ᶜinv_ρ̂ = (@. lazy(
+                divide_by_ρa(
+                    FT(1),
+                    Y.c.sgsʲs.:($$j).ρa,
+                    FT(0),
+                    Y.c.ρ,
+                    turbconv_model,
+                ),
+            ))
+            ᶜ∂ρ̂∂zₗ = (@. lazy(
+                upwind_biased_grad(
+                    -1 * Geometry.WVector(ᶜwₗʲs.:($$j)),
+                    Y.c.sgsʲs.:($$j).ρa,
+                ),
+            ))
+            ᶜ∂ρ̂∂zᵢ = (@. lazy(
+                upwind_biased_grad(
+                    -1 * Geometry.WVector(ᶜwᵢʲs.:($$j)),
+                    Y.c.sgsʲs.:($$j).ρa,
+                ),
+            ))
+            ᶜ∂ρ̂∂zᵣ = (@. lazy(
+                upwind_biased_grad(
+                    -1 * Geometry.WVector(ᶜwᵣʲs.:($$j)),
+                    Y.c.sgsʲs.:($$j).ρa,
+                ),
+            ))
+            ᶜ∂ρ̂∂zₛ = (@. lazy(
+                upwind_biased_grad(
+                    -1 * Geometry.WVector(ᶜwₛʲs.:($$j)),
+                    Y.c.sgsʲs.:($$j).ρa,
+                ),
+            ))
+
+            @. Yₜ.c.sgsʲs.:($$j).q_tot -=
+                dot(
+                    ᶜinv_ρ̂ * ᶜ∂ρ̂∂zₗ,
+                    CT3(Geometry.WVector(-1 * ᶜwₗʲs.:($$j))),
+                ) *
+                Y.c.sgsʲs.:($$j).q_liq *
+                (1 - Y.c.sgsʲs.:($$j).q_tot)
+            @. Yₜ.c.sgsʲs.:($$j).q_liq -=
+                dot(
+                    ᶜinv_ρ̂ * ᶜ∂ρ̂∂zₗ,
+                    CT3(Geometry.WVector(-1 * ᶜwₗʲs.:($$j))),
+                ) * Y.c.sgsʲs.:($$j).q_liq
+            @. Yₜ.c.sgsʲs.:($$j).mse -=
+                dot(
+                    ᶜinv_ρ̂ * ᶜ∂ρ̂∂zₗ,
+                    CT3(Geometry.WVector(-1 * ᶜwₗʲs.:($$j))),
+                ) * (Y.c.sgsʲs.:($$j).q_liq * Y.c.sgsʲs.:($$j).mse - ᶜmseₗ)
+
+            @. Yₜ.c.sgsʲs.:($$j).q_tot -=
+                dot(
+                    ᶜinv_ρ̂ * ᶜ∂ρ̂∂zᵢ,
+                    CT3(Geometry.WVector(-1 * ᶜwᵢʲs.:($$j))),
+                ) *
+                Y.c.sgsʲs.:($$j).q_ice *
+                (1 - Y.c.sgsʲs.:($$j).q_tot)
+            @. Yₜ.c.sgsʲs.:($$j).q_ice -=
+                dot(
+                    ᶜinv_ρ̂ * ᶜ∂ρ̂∂zᵢ,
+                    CT3(Geometry.WVector(-1 * ᶜwᵢʲs.:($$j))),
+                ) * Y.c.sgsʲs.:($$j).q_ice
+            @. Yₜ.c.sgsʲs.:($$j).mse -=
+                dot(
+                    ᶜinv_ρ̂ * ᶜ∂ρ̂∂zᵢ,
+                    CT3(Geometry.WVector(-1 * ᶜwᵢʲs.:($$j))),
+                ) * (Y.c.sgsʲs.:($$j).q_ice * Y.c.sgsʲs.:($$j).mse - ᶜmseᵢ)
+
+            @. Yₜ.c.sgsʲs.:($$j).q_tot -=
+                dot(
+                    ᶜinv_ρ̂ * ᶜ∂ρ̂∂zᵣ,
+                    CT3(Geometry.WVector(-1 * ᶜwᵣʲs.:($$j))),
+                ) *
+                Y.c.sgsʲs.:($$j).q_rai *
+                (1 - Y.c.sgsʲs.:($$j).q_tot)
+            @. Yₜ.c.sgsʲs.:($$j).q_rai -=
+                dot(
+                    ᶜinv_ρ̂ * ᶜ∂ρ̂∂zᵣ,
+                    CT3(Geometry.WVector(-1 * ᶜwᵣʲs.:($$j))),
+                ) * Y.c.sgsʲs.:($$j).q_rai
+            @. Yₜ.c.sgsʲs.:($$j).mse -=
+                dot(
+                    ᶜinv_ρ̂ * ᶜ∂ρ̂∂zᵣ,
+                    CT3(Geometry.WVector(-1 * ᶜwᵣʲs.:($$j))),
+                ) * (Y.c.sgsʲs.:($$j).q_rai * Y.c.sgsʲs.:($$j).mse - ᶜmseᵣ)
+
+            @. Yₜ.c.sgsʲs.:($$j).q_tot -=
+                dot(
+                    ᶜinv_ρ̂ * ᶜ∂ρ̂∂zₛ,
+                    CT3(Geometry.WVector(-1 * ᶜwₛʲs.:($$j))),
+                ) *
+                Y.c.sgsʲs.:($$j).q_sno *
+                (1 - Y.c.sgsʲs.:($$j).q_tot)
+            @. Yₜ.c.sgsʲs.:($$j).q_sno -=
+                dot(
+                    ᶜinv_ρ̂ * ᶜ∂ρ̂∂zₛ,
+                    CT3(Geometry.WVector(-1 * ᶜwₛʲs.:($$j))),
+                ) * Y.c.sgsʲs.:($$j).q_sno
+            @. Yₜ.c.sgsʲs.:($$j).mse -=
+                dot(
+                    ᶜinv_ρ̂ * ᶜ∂ρ̂∂zₛ,
+                    CT3(Geometry.WVector(-1 * ᶜwₛʲs.:($$j))),
+                ) * (Y.c.sgsʲs.:($$j).q_sno * Y.c.sgsʲs.:($$j).mse - ᶜmseₛ)
+
+            # mse, q_tot and moisture tracer terms proportional to velocity gradients
+            @. Yₜ.c.sgsʲs.:($$j).q_tot -=
+                ᶜdivᵥ(ᶠwₗ³ʲs) *
+                Y.c.sgsʲs.:($$j).q_liq *
+                (1 - Y.c.sgsʲs.:($$j).q_tot)
+            @. Yₜ.c.sgsʲs.:($$j).q_liq -= ᶜdivᵥ(ᶠwₗ³ʲs) * Y.c.sgsʲs.:($$j).q_liq
+            @. Yₜ.c.sgsʲs.:($$j).mse -=
+                ᶜdivᵥ(ᶠwₗ³ʲs) *
+                (Y.c.sgsʲs.:($$j).q_liq * Y.c.sgsʲs.:($$j).mse - ᶜmseₗ)
+
+            @. Yₜ.c.sgsʲs.:($$j).q_tot -=
+                ᶜdivᵥ(ᶠwᵢ³ʲs) *
+                Y.c.sgsʲs.:($$j).q_ice *
+                (1 - Y.c.sgsʲs.:($$j).q_tot)
+            @. Yₜ.c.sgsʲs.:($$j).q_ice -= ᶜdivᵥ(ᶠwᵢ³ʲs) * Y.c.sgsʲs.:($$j).q_ice
+            @. Yₜ.c.sgsʲs.:($$j).mse -=
+                ᶜdivᵥ(ᶠwᵢ³ʲs) *
+                (Y.c.sgsʲs.:($$j).q_ice * Y.c.sgsʲs.:($$j).mse - ᶜmseᵢ)
+
+            @. Yₜ.c.sgsʲs.:($$j).q_tot -=
+                ᶜdivᵥ(ᶠwᵣ³ʲs) *
+                Y.c.sgsʲs.:($$j).q_rai *
+                (1 - Y.c.sgsʲs.:($$j).q_tot)
+            @. Yₜ.c.sgsʲs.:($$j).q_rai -= ᶜdivᵥ(ᶠwᵣ³ʲs) * Y.c.sgsʲs.:($$j).q_rai
+            @. Yₜ.c.sgsʲs.:($$j).mse -=
+                ᶜdivᵥ(ᶠwᵣ³ʲs) *
+                (Y.c.sgsʲs.:($$j).q_rai * Y.c.sgsʲs.:($$j).mse - ᶜmseᵣ)
+
+            @. Yₜ.c.sgsʲs.:($$j).q_tot -=
+                ᶜdivᵥ(ᶠwₛ³ʲs) *
+                Y.c.sgsʲs.:($$j).q_sno *
+                (1 - Y.c.sgsʲs.:($$j).q_tot)
+            @. Yₜ.c.sgsʲs.:($$j).q_sno -= ᶜdivᵥ(ᶠwₛ³ʲs) * Y.c.sgsʲs.:($$j).q_sno
+            @. Yₜ.c.sgsʲs.:($$j).mse -=
+                ᶜdivᵥ(ᶠwₛ³ʲs) *
+                (Y.c.sgsʲs.:($$j).q_sno * Y.c.sgsʲs.:($$j).mse - ᶜmseₛ)
         end
     end
 end
diff --git a/src/prognostic_equations/implicit/implicit_tendency.jl b/src/prognostic_equations/implicit/implicit_tendency.jl
index 1468746c67..f4fd137053 100644
--- a/src/prognostic_equations/implicit/implicit_tendency.jl
+++ b/src/prognostic_equations/implicit/implicit_tendency.jl
@@ -70,32 +70,32 @@ end
 # expressions are less convoluted?
 
 function vertical_transport(ᶜρ, ᶠu³, ᶜχ, dt, ::Val{:none})
-    ᶜJ = Fields.local_geometry_field(ᶜρ).J
-    ᶠJ = Fields.local_geometry_field(ᶠu³).J
+    ᶜJ = Fields.local_geometry_field(axes(ᶜρ)).J
+    ᶠJ = Fields.local_geometry_field(axes(ᶠu³)).J
     return @. lazy(-(ᶜadvdivᵥ(ᶠinterp(ᶜρ * ᶜJ) / ᶠJ * ᶠu³ * ᶠinterp(ᶜχ))))
 end
 function vertical_transport(ᶜρ, ᶠu³, ᶜχ, dt, ::Val{:first_order})
-    ᶜJ = Fields.local_geometry_field(ᶜρ).J
-    ᶠJ = Fields.local_geometry_field(ᶠu³).J
+    ᶜJ = Fields.local_geometry_field(axes(ᶜρ)).J
+    ᶠJ = Fields.local_geometry_field(axes(ᶠu³)).J
     return @. lazy(-(ᶜadvdivᵥ(ᶠinterp(ᶜρ * ᶜJ) / ᶠJ * ᶠupwind1(ᶠu³, ᶜχ))))
 end
 @static if pkgversion(ClimaCore) ≥ v"0.14.22"
     function vertical_transport(ᶜρ, ᶠu³, ᶜχ, dt, ::Val{:vanleer_limiter})
-        ᶜJ = Fields.local_geometry_field(ᶜρ).J
-        ᶠJ = Fields.local_geometry_field(ᶠu³).J
+        ᶜ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))),
         )
     end
 end
 function vertical_transport(ᶜρ, ᶠu³, ᶜχ, dt, ::Val{:third_order})
-    ᶜJ = Fields.local_geometry_field(ᶜρ).J
-    ᶠJ = Fields.local_geometry_field(ᶠu³).J
+    ᶜJ = Fields.local_geometry_field(axes(ᶜρ)).J
+    ᶠJ = Fields.local_geometry_field(axes(ᶠu³)).J
     return @. lazy(-(ᶜadvdivᵥ(ᶠinterp(ᶜρ * ᶜJ) / ᶠJ * ᶠupwind3(ᶠu³, ᶜχ))))
 end
 function vertical_transport(ᶜρ, ᶠu³, ᶜχ, dt, ::Val{:boris_book})
-    ᶜJ = Fields.local_geometry_field(ᶜρ).J
-    ᶠJ = Fields.local_geometry_field(ᶠu³).J
+    ᶜJ = Fields.local_geometry_field(axes(ᶜρ)).J
+    ᶠJ = Fields.local_geometry_field(axes(ᶠu³)).J
     return @. lazy(
         -(ᶜadvdivᵥ(
             ᶠinterp(ᶜρ * ᶜJ) / ᶠJ * (
@@ -109,8 +109,8 @@ function vertical_transport(ᶜρ, ᶠu³, ᶜχ, dt, ::Val{:boris_book})
     )
 end
 function vertical_transport(ᶜρ, ᶠu³, ᶜχ, dt, ::Val{:zalesak})
-    ᶜJ = Fields.local_geometry_field(ᶜρ).J
-    ᶠJ = Fields.local_geometry_field(ᶠu³).J
+    ᶜJ = Fields.local_geometry_field(axes(ᶜρ)).J
+    ᶠJ = Fields.local_geometry_field(axes(ᶠu³)).J
     return @. lazy(
         -(ᶜadvdivᵥ(
             ᶠinterp(ᶜρ * ᶜJ) / ᶠJ * (
@@ -137,8 +137,8 @@ function implicit_vertical_advection_tendency!(Yₜ, Y, p, t)
         p.atmos
     (; dt) = p
     n = n_mass_flux_subdomains(turbconv_model)
-    ᶜJ = Fields.local_geometry_field(Y.c).J
-    ᶠJ = Fields.local_geometry_field(Y.f).J
+    ᶜJ = Fields.local_geometry_field(axes(Y.c)).J
+    ᶠJ = Fields.local_geometry_field(axes(Y.f)).J
     (; ᶠgradᵥ_ᶜΦ) = p.core
     (; ᶜh_tot, ᶠu³, ᶜp) = p.precomputed
 
diff --git a/src/utils/abbreviations.jl b/src/utils/abbreviations.jl
index 186bcedf39..a52ecef711 100644
--- a/src/utils/abbreviations.jl
+++ b/src/utils/abbreviations.jl
@@ -72,7 +72,7 @@ const ᶠcurlᵥ = Operators.CurlC2F(
     bottom = Operators.SetCurl(CT12(0, 0)),
     top = Operators.SetCurl(CT12(0, 0)),
 )
-
+const upwind_biased_grad = Operators.UpwindBiasedGradient()
 const ᶠupwind1 = Operators.UpwindBiasedProductC2F()
 const ᶠupwind3 = Operators.Upwind3rdOrderBiasedProductC2F(
     bottom = Operators.ThirdOrderOneSided(),