Fix kinetic energy conservation with topography

src/prognostic_equations/advection.jl

@@ -73,6 +73,7 @@ NVTX.@annotate function explicit_vertical_advection_tendency!(Yₜ, Y, p, t)
     point_type = eltype(Fields.coordinate_field(Y.c))
     (; dt) = p
     ᶜJ = Fields.local_geometry_field(Y.c).J
+    ᶠJ = Fields.local_geometry_field(Y.f).J
     (; ᶜf³, ᶠf¹², ᶜΦ) = p.core
     (; ᶜu, ᶠu³, ᶜK) = p.precomputed
     (; edmfx_upwinding) = n > 0 || advect_tke ? p.atmos.numerics : all_nothing
@@ -116,6 +117,7 @@ NVTX.@annotate function explicit_vertical_advection_tendency!(Yₜ, Y, p, t)
                     coeff,
                     Yₜ.c.ρe_tot[colidx],
                     ᶜJ[colidx],
+                    ᶠJ[colidx],
                     Y.c.ρ[colidx],
                     ᶠu³[colidx],
                     ᶜh_tot[colidx],
@@ -132,6 +134,7 @@ NVTX.@annotate function explicit_vertical_advection_tendency!(Yₜ, Y, p, t)
                     coeff,
                     ᶜρχₜ[colidx],
                     ᶜJ[colidx],
+                    ᶠJ[colidx],
                     Y.c.ρ[colidx],
                     ᶠu³[colidx],
                     ᶜχ[colidx],
@@ -184,6 +187,7 @@ NVTX.@annotate function explicit_vertical_advection_tendency!(Yₜ, Y, p, t)
             vertical_transport!(
                 Yₜ.c.sgs⁰.ρatke[colidx],
                 ᶜJ[colidx],
+                ᶠJ[colidx],
                 ᶜρ⁰[colidx],
                 ᶠu³⁰[colidx],
                 ᶜa_scalar[colidx],
@@ -209,6 +213,7 @@ function edmfx_sgs_vertical_advection_tendency!(
     n = n_prognostic_mass_flux_subdomains(turbconv_model)
     (; dt) = p
     ᶜJ = Fields.local_geometry_field(Y.c).J
+    ᶠJ = Fields.local_geometry_field(Y.f).J
     (; edmfx_upwinding) = p.atmos.numerics
     (; ᶠu³ʲs, ᶠKᵥʲs, ᶜρʲs) = p.precomputed
     (; ᶠgradᵥ_ᶜΦ) = p.core
@@ -244,6 +249,7 @@ function edmfx_sgs_vertical_advection_tendency!(
         vertical_transport!(
             Yₜ.c.sgsʲs.:($j).ρa[colidx],
             ᶜJ[colidx],
+            ᶠJ[colidx],
             ᶜρʲs.:($j)[colidx],
             ᶠu³ʲs.:($j)[colidx],
             ᶜa_scalar[colidx],

src/prognostic_equations/edmfx_sgs_flux.jl

@@ -20,6 +20,7 @@ function edmfx_sgs_mass_flux_tendency!(
     (; ᶜρa⁰, ᶜρ⁰, ᶠu³⁰, ᶜK⁰, ᶜmse⁰, ᶜq_tot⁰) = p.precomputed
     (; dt) = p
     ᶜJ = Fields.local_geometry_field(Y.c).J
+    ᶠJ = Fields.local_geometry_field(Y.f).J
 
     if p.atmos.edmfx_sgs_mass_flux
         # energy
@@ -35,6 +36,7 @@ function edmfx_sgs_mass_flux_tendency!(
             vertical_transport!(
                 Yₜ.c.ρe_tot[colidx],
                 ᶜJ[colidx],
+                ᶠJ[colidx],
                 ᶜρʲs.:($j)[colidx],
                 ᶠu³_diff_colidx,
                 ᶜa_scalar_colidx,
@@ -49,6 +51,7 @@ function edmfx_sgs_mass_flux_tendency!(
         vertical_transport!(
             Yₜ.c.ρe_tot[colidx],
             ᶜJ[colidx],
+            ᶠJ[colidx],
             ᶜρ⁰[colidx],
             ᶠu³_diff_colidx,
             ᶜa_scalar_colidx,
@@ -66,6 +69,7 @@ function edmfx_sgs_mass_flux_tendency!(
                 vertical_transport!(
                     Yₜ.c.ρq_tot[colidx],
                     ᶜJ[colidx],
+                    ᶠJ[colidx],
                     ᶜρʲs.:($j)[colidx],
                     ᶠu³_diff_colidx,
                     ᶜa_scalar_colidx,
@@ -80,6 +84,7 @@ function edmfx_sgs_mass_flux_tendency!(
             vertical_transport!(
                 Yₜ.c.ρq_tot[colidx],
                 ᶜJ[colidx],
+                ᶠJ[colidx],
                 ᶜρ⁰[colidx],
                 ᶠu³_diff_colidx,
                 ᶜa_scalar_colidx,
@@ -111,6 +116,7 @@ function edmfx_sgs_mass_flux_tendency!(
     (; ᶜρaʲs, ᶜρʲs, ᶠu³ʲs, ᶜKʲs, ᶜmseʲs, ᶜq_totʲs) = p.precomputed
     (; dt) = p
     ᶜJ = Fields.local_geometry_field(Y.c).J
+    ᶠJ = Fields.local_geometry_field(Y.f).J
     FT = eltype(Y)
 
     if p.atmos.edmfx_sgs_mass_flux
@@ -140,6 +146,7 @@ function edmfx_sgs_mass_flux_tendency!(
             vertical_transport!(
                 Yₜ.c.ρe_tot[colidx],
                 ᶜJ[colidx],
+                ᶠJ[colidx],
                 ᶜρʲs.:($j)[colidx],
                 ᶠu³_diff_colidx,
                 ᶜa_scalar_colidx,
@@ -175,6 +182,7 @@ function edmfx_sgs_mass_flux_tendency!(
                 vertical_transport!(
                     Yₜ.c.ρq_tot[colidx],
                     ᶜJ[colidx],
+                    ᶠJ[colidx],
                     ᶜρʲs.:($j)[colidx],
                     ᶠu³_diff_colidx,
                     ᶜa_scalar_colidx,

src/prognostic_equations/implicit/implicit_solver.jl

@@ -499,6 +499,7 @@ function update_implicit_equation_jacobian!(A, Y, p, dtγ, colidx)
     ᶜuₕ = Y.c.uₕ
     ᶠu₃ = Y.f.u₃
     ᶜJ = Fields.local_geometry_field(Y.c).J
+    ᶠJ = Fields.local_geometry_field(Y.f).J
     ᶜgⁱʲ = Fields.local_geometry_field(Y.c).gⁱʲ
     ᶠgⁱʲ = Fields.local_geometry_field(Y.f).gⁱʲ
 
@@ -507,14 +508,8 @@ function update_implicit_equation_jacobian!(A, Y, p, dtγ, colidx)
         TD.gas_constant_air(thermo_params, ᶜts[colidx]) /
         TD.cv_m(thermo_params, ᶜts[colidx])
 
-    if use_derivative(topography_flag)
-        @. ∂ᶜK_∂ᶜuₕ[colidx] = DiagonalMatrixRow(
-            adjoint(CTh(ᶜuₕ[colidx])) +
-            adjoint(ᶜinterp(ᶠu₃[colidx])) * g³ʰ(ᶜgⁱʲ[colidx]),
-        )
-    else
-        @. ∂ᶜK_∂ᶜuₕ[colidx] = DiagonalMatrixRow(adjoint(CTh(ᶜuₕ[colidx])))
-    end
+    @. ∂ᶜK_∂ᶜuₕ[colidx] =
+        DiagonalMatrixRow(adjoint(CTh(ᶜuₕ[colidx]) + CTh(ᶜinterp(ᶠu₃[colidx]))))
     @. ∂ᶜK_∂ᶠu₃[colidx] =
         ᶜinterp_matrix() ⋅ DiagonalMatrixRow(adjoint(CT3(ᶠu₃[colidx]))) +
         DiagonalMatrixRow(adjoint(CT3(ᶜuₕ[colidx]))) ⋅ ᶜinterp_matrix()
@@ -524,7 +519,7 @@ function update_implicit_equation_jacobian!(A, Y, p, dtγ, colidx)
 
     @. ᶜadvection_matrix[colidx] =
         -(ᶜadvdivᵥ_matrix()) ⋅
-        DiagonalMatrixRow(ᶠwinterp(ᶜJ[colidx], ᶜρ[colidx]))
+        DiagonalMatrixRow(ᶠinterp(ᶜJ[colidx] * ᶜρ[colidx]) / ᶠJ[colidx])
 
     if use_derivative(topography_flag)
         ∂ᶜρ_err_∂ᶜuₕ = matrix[@name(c.ρ), @name(c.uₕ)]
@@ -716,7 +711,7 @@ function update_implicit_equation_jacobian!(A, Y, p, dtγ, colidx)
             ᶠtmp = p.ᶠtemp_CT3
             @. ᶠtmp[colidx] =
                 CT3(unit_basis_vector_data(CT3, ᶠlg[colidx])) *
-                ᶠwinterp(ᶜJ[colidx], ᶜρ[colidx])
+                ᶠinterp(ᶜJ[colidx] * ᶜρ[colidx]) / ᶠJ[colidx]
             @. ∂ᶜρqₚ_err_∂ᶜρqₚ[colidx] +=
                 dtγ * -(ᶜprecipdivᵥ_matrix()) ⋅ DiagonalMatrixRow(ᶠtmp) ⋅
                 ᶠright_bias_matrix() ⋅
@@ -802,14 +797,15 @@ function update_implicit_equation_jacobian!(A, Y, p, dtγ, colidx)
                                 ᶜρʲs.:(1)[colidx],
                             ),
                         ),
-                    ),
-                ) ⋅ ᶠwinterp_matrix(ᶜJ[colidx]) ⋅ DiagonalMatrixRow(
-                    ᶜkappa_mʲ[colidx] * (ᶜρʲs.:(1)[colidx])^2 /
+                    ) / ᶠJ[colidx],
+                ) ⋅ ᶠinterp_matrix() ⋅ DiagonalMatrixRow(
+                    ᶜJ[colidx] * ᶜkappa_mʲ[colidx] * (ᶜρʲs.:(1)[colidx])^2 /
                     ((ᶜkappa_mʲ[colidx] + 1) * ᶜp[colidx]) * e_int_v0,
                 )
             @. ᶠbidiagonal_matrix_ct3_2[colidx] =
-                DiagonalMatrixRow(ᶠwinterp(ᶜJ[colidx], ᶜρʲs.:(1)[colidx])) ⋅
-                ᶠset_upwind_matrix_bcs(ᶠupwind_matrix(ᶠu³ʲs.:(1)[colidx])) ⋅
+                DiagonalMatrixRow(
+                    ᶠinterp(ᶜJ[colidx] * ᶜρʲs.:(1)[colidx]) / ᶠJ[colidx],
+                ) ⋅ ᶠset_upwind_matrix_bcs(ᶠupwind_matrix(ᶠu³ʲs.:(1)[colidx])) ⋅
                 DiagonalMatrixRow(
                     Y.c.sgsʲs.:(1).ρa[colidx] * ᶜkappa_mʲ[colidx] /
                     ((ᶜkappa_mʲ[colidx] + 1) * ᶜp[colidx]) * e_int_v0,
@@ -832,14 +828,15 @@ function update_implicit_equation_jacobian!(A, Y, p, dtγ, colidx)
                                 ᶜρʲs.:(1)[colidx],
                             ),
                         ),
-                    ),
-                ) ⋅ ᶠwinterp_matrix(ᶜJ[colidx]) ⋅ DiagonalMatrixRow(
-                    ᶜkappa_mʲ[colidx] * (ᶜρʲs.:(1)[colidx])^2 /
+                    ) / ᶠJ[colidx],
+                ) ⋅ ᶠinterp_matrix() ⋅ DiagonalMatrixRow(
+                    ᶜJ[colidx] * ᶜkappa_mʲ[colidx] * (ᶜρʲs.:(1)[colidx])^2 /
                     ((ᶜkappa_mʲ[colidx] + 1) * ᶜp[colidx]),
                 )
             @. ᶠbidiagonal_matrix_ct3_2[colidx] =
-                DiagonalMatrixRow(ᶠwinterp(ᶜJ[colidx], ᶜρʲs.:(1)[colidx])) ⋅
-                ᶠset_upwind_matrix_bcs(ᶠupwind_matrix(ᶠu³ʲs.:(1)[colidx])) ⋅
+                DiagonalMatrixRow(
+                    ᶠinterp(ᶜJ[colidx] * ᶜρʲs.:(1)[colidx]) / ᶠJ[colidx],
+                ) ⋅ ᶠset_upwind_matrix_bcs(ᶠupwind_matrix(ᶠu³ʲs.:(1)[colidx])) ⋅
                 DiagonalMatrixRow(
                     Y.c.sgsʲs.:(1).ρa[colidx] * ᶜkappa_mʲ[colidx] /
                     ((ᶜkappa_mʲ[colidx] + 1) * ᶜp[colidx]),
@@ -852,8 +849,9 @@ function update_implicit_equation_jacobian!(A, Y, p, dtγ, colidx)
             ∂ᶜρaʲ_err_∂ᶜρaʲ =
                 matrix[@name(c.sgsʲs.:(1).ρa), @name(c.sgsʲs.:(1).ρa)]
             @. ᶜadvection_matrix[colidx] =
-                -(ᶜadvdivᵥ_matrix()) ⋅
-                DiagonalMatrixRow(ᶠwinterp(ᶜJ[colidx], ᶜρʲs.:(1)[colidx]))
+                -(ᶜadvdivᵥ_matrix()) ⋅ DiagonalMatrixRow(
+                    ᶠinterp(ᶜJ[colidx] * ᶜρʲs.:(1)[colidx]) / ᶠJ[colidx],
+                )
             @. ∂ᶜρaʲ_err_∂ᶜρaʲ[colidx] =
                 dtγ * ᶜadvection_matrix[colidx] ⋅
                 ᶠset_upwind_matrix_bcs(ᶠupwind_matrix(ᶠu³ʲs.:(1)[colidx])) ⋅

src/prognostic_equations/implicit/implicit_tendency.jl

@@ -56,49 +56,64 @@ end
 # the implicit tendency function. Since dt >= dtγ, we can safely use dt for now.
 # TODO: Can we rewrite ᶠfct_boris_book and ᶠfct_zalesak so that their broadcast
 # expressions are less convoluted?
-vertical_transport!(ᶜρχₜ, ᶜJ, ᶜρ, ᶠu³, ᶜχ, dt, upwinding::Val, ᶜdivᵥ) =
-    vertical_transport!(1, ᶜρχₜ, ᶜJ, ᶜρ, ᶠu³, ᶜχ, dt, upwinding, ᶜdivᵥ)
-vertical_transport!(ᶜρχₜ, ᶜJ, ᶜρ, ᶠu³, ᶜχ, dt, upwinding::Val) =
-    vertical_transport!(1, ᶜρχₜ, ᶜJ, ᶜρ, ᶠu³, ᶜχ, dt, upwinding, ᶜadvdivᵥ)
+vertical_transport!(ᶜρχₜ, ᶜJ, ᶠJ, ᶜρ, ᶠu³, ᶜχ, dt, upwinding::Val, ᶜdivᵥ) =
+    vertical_transport!(1, ᶜρχₜ, ᶜJ, ᶠJ, ᶜρ, ᶠu³, ᶜχ, dt, upwinding, ᶜdivᵥ)
+vertical_transport!(ᶜρχₜ, ᶜJ, ᶠJ, ᶜρ, ᶠu³, ᶜχ, dt, upwinding::Val) =
+    vertical_transport!(1, ᶜρχₜ, ᶜJ, ᶠJ, ᶜρ, ᶠu³, ᶜχ, dt, upwinding, ᶜadvdivᵥ)
 vertical_transport!(
     coeff::Int,
     ᶜρχₜ,
     ᶜJ,
+    ᶠJ,
     ᶜρ,
     ᶠu³,
     ᶜχ,
     dt::Real,
     upwinding::Val,
-) = vertical_transport!(coeff, ᶜρχₜ, ᶜJ, ᶜρ, ᶠu³, ᶜχ, dt, upwinding, ᶜadvdivᵥ)
+) = vertical_transport!(
+    coeff,
+    ᶜρχₜ,
+    ᶜJ,
+    ᶠJ,
+    ᶜρ,
+    ᶠu³,
+    ᶜχ,
+    dt,
+    upwinding,
+    ᶜadvdivᵥ,
+)
 
-vertical_transport!(coeff, ᶜρχₜ, ᶜJ, ᶜρ, ᶠu³, ᶜχ, dt, ::Val{:none}, ᶜdivᵥ) =
-    @. ᶜρχₜ += -coeff * (ᶜdivᵥ(ᶠwinterp(ᶜJ, ᶜρ) * ᶠu³ * ᶠinterp(ᶜχ)))
+vertical_transport!(coeff, ᶜρχₜ, ᶜJ, ᶠJ, ᶜρ, ᶠu³, ᶜχ, dt, ::Val{:none}, ᶜdivᵥ) =
+    @. ᶜρχₜ += -coeff * (ᶜdivᵥ(ᶠinterp(ᶜJ * ᶜρ) * ᶠu³ * ᶠinterp(ᶜχ) / ᶠJ))
 vertical_transport!(
     coeff,
     ᶜρχₜ,
     ᶜJ,
+    ᶠJ,
     ᶜρ,
     ᶠu³,
     ᶜχ,
     dt,
     ::Val{:first_order},
     ᶜdivᵥ,
-) = @. ᶜρχₜ += -coeff * (ᶜdivᵥ(ᶠwinterp(ᶜJ, ᶜρ) * ᶠupwind1(ᶠu³, ᶜχ)))
+) = @. ᶜρχₜ += -coeff * (ᶜdivᵥ(ᶠinterp(ᶜJ * ᶜρ) * ᶠupwind1(ᶠu³, ᶜχ) / ᶠJ))
 vertical_transport!(
     coeff,
     ᶜρχₜ,
     ᶜJ,
+    ᶠJ,
     ᶜρ,
     ᶠu³,
     ᶜχ,
     dt,
     ::Val{:third_order},
     ᶜdivᵥ,
-) = @. ᶜρχₜ += -coeff * (ᶜdivᵥ(ᶠwinterp(ᶜJ, ᶜρ) * ᶠupwind3(ᶠu³, ᶜχ)))
+) = @. ᶜρχₜ += -coeff * (ᶜdivᵥ(ᶠinterp(ᶜJ * ᶜρ) * ᶠupwind3(ᶠu³, ᶜχ) / ᶠJ))
 vertical_transport!(
     coeff,
     ᶜρχₜ,
     ᶜJ,
+    ᶠJ,
     ᶜρ,
     ᶠu³,
     ᶜχ,
@@ -107,24 +122,34 @@ vertical_transport!(
     ᶜdivᵥ,
 ) = @. ᶜρχₜ +=
     -coeff * (ᶜdivᵥ(
-        ᶠwinterp(ᶜJ, ᶜρ) * (
+        ᶠinterp(ᶜJ * ᶜρ) * (
             ᶠupwind1(ᶠu³, ᶜχ) + ᶠfct_boris_book(
                 ᶠupwind3(ᶠu³, ᶜχ) - ᶠupwind1(ᶠu³, ᶜχ),
-                ᶜχ / dt - ᶜdivᵥ(ᶠwinterp(ᶜJ, ᶜρ) * ᶠupwind1(ᶠu³, ᶜχ)) / ᶜρ,
+                ᶜχ / dt - ᶜdivᵥ(ᶠinterp(ᶜJ * ᶜρ) * ᶠupwind1(ᶠu³, ᶜχ) / ᶠJ) / ᶜρ,
+            )
+        ) / ᶠJ,
+    ))
+vertical_transport!(
+    coeff,
+    ᶜρχₜ,
+    ᶜJ,
+    ᶠJ,
+    ᶜρ,
+    ᶠu³,
+    ᶜχ,
+    dt,
+    ::Val{:zalesak},
+    ᶜdivᵥ,
+) = @. ᶜρχₜ +=
+    -coeff * (ᶜdivᵥ(
+        ᶠinterp(ᶜJ * ᶜρ) * (
+            ᶠupwind1(ᶠu³, ᶜχ) + ᶠfct_zalesak(
+                ᶠupwind3(ᶠu³, ᶜχ) - ᶠupwind1(ᶠu³, ᶜχ),
+                ᶜχ / dt,
+                ᶜχ / dt - ᶜdivᵥ(ᶠinterp(ᶜJ * ᶜρ) * ᶠupwind1(ᶠu³, ᶜχ) / ᶠJ) / ᶜρ,
             )
-        ),
+        ) / ᶠJ,
     ))
-vertical_transport!(coeff, ᶜρχₜ, ᶜJ, ᶜρ, ᶠu³, ᶜχ, dt, ::Val{:zalesak}, ᶜdivᵥ) =
-    @. ᶜρχₜ +=
-        -coeff * (ᶜdivᵥ(
-            ᶠwinterp(ᶜJ, ᶜρ) * (
-                ᶠupwind1(ᶠu³, ᶜχ) + ᶠfct_zalesak(
-                    ᶠupwind3(ᶠu³, ᶜχ) - ᶠupwind1(ᶠu³, ᶜχ),
-                    ᶜχ / dt,
-                    ᶜχ / dt - ᶜdivᵥ(ᶠwinterp(ᶜJ, ᶜρ) * ᶠupwind1(ᶠu³, ᶜχ)) / ᶜρ,
-                )
-            ),
-        ))
 
 vertical_advection!(ᶜρχₜ, ᶠu³, ᶜχ, ::Val{:none}) =
     @. ᶜρχₜ -= ᶜadvdivᵥ(ᶠu³ * ᶠinterp(ᶜχ)) - ᶜχ * ᶜadvdivᵥ(ᶠu³)
@@ -138,16 +163,18 @@ function implicit_vertical_advection_tendency!(Yₜ, Y, p, t, colidx)
     (; 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
     (; ᶠgradᵥ_ᶜΦ, ᶜρ_ref, ᶜp_ref) = p.core
     (; ᶜh_tot, ᶜspecific, ᶠu³, ᶜp) = p.precomputed
 
     @. Yₜ.c.ρ[colidx] -=
-        ᶜdivᵥ(ᶠwinterp(ᶜJ[colidx], Y.c.ρ[colidx]) * ᶠu³[colidx])
+        ᶜdivᵥ(ᶠinterp(ᶜJ[colidx] * Y.c.ρ[colidx]) * ᶠu³[colidx] / ᶠJ[colidx])
 
     # Central advection of active tracers (e_tot and q_tot)
     vertical_transport!(
         Yₜ.c.ρe_tot[colidx],
         ᶜJ[colidx],
+        ᶠJ[colidx],
         Y.c.ρ[colidx],
         ᶠu³[colidx],
         ᶜh_tot[colidx],
@@ -158,6 +185,7 @@ function implicit_vertical_advection_tendency!(Yₜ, Y, p, t, colidx)
         vertical_transport!(
             Yₜ.c.ρq_tot[colidx],
             ᶜJ[colidx],
+            ᶠJ[colidx],
             Y.c.ρ[colidx],
             ᶠu³[colidx],
             ᶜspecific.q_tot[colidx],
@@ -175,13 +203,13 @@ function implicit_vertical_advection_tendency!(Yₜ, Y, p, t, colidx)
         ᶠlg = Fields.local_geometry_field(Y.f)
         @. Yₜ.c.ρq_rai[colidx] -= ᶜprecipdivᵥ(
             CT3(unit_basis_vector_data(CT3, ᶠlg[colidx])) *
-            ᶠwinterp(ᶜJ[colidx], Y.c.ρ[colidx]) *
-            ᶠright_bias(-p.precomputed.ᶜwᵣ[colidx] * ᶜspecific.q_rai[colidx]),
+            ᶠinterp(ᶜJ[colidx] * Y.c.ρ[colidx]) *
+            ᶠright_bias(-p.precomputed.ᶜwᵣ[colidx] * ᶜspecific.q_rai[colidx]) / ᶠJ[colidx],
         )
         @. Yₜ.c.ρq_sno[colidx] -= ᶜprecipdivᵥ(
             CT3(unit_basis_vector_data(CT3, ᶠlg[colidx])) *
-            ᶠwinterp(ᶜJ[colidx], Y.c.ρ[colidx]) *
-            ᶠright_bias(-p.precomputed.ᶜwₛ[colidx] * ᶜspecific.q_sno[colidx]),
+            ᶠinterp(ᶜJ[colidx] * Y.c.ρ[colidx]) *
+            ᶠright_bias(-p.precomputed.ᶜwₛ[colidx] * ᶜspecific.q_sno[colidx]) / ᶠJ[colidx],
         )
     end
 

src/utils/utilities.jl

@@ -61,11 +61,11 @@ function compute_kinetic!(κ::Fields.Field, uₕ::Fields.Field, uᵥ::Fields.Fie
     @assert eltype(uₕ) <: Union{C1, C2, C12}
     @assert eltype(uᵥ) <: C3
     @. κ =
-        1 / 2 * (
-            dot(C123(uₕ), CT123(uₕ)) +
-            ᶜinterp(dot(C123(uᵥ), CT123(uᵥ))) +
-            2 * dot(CT123(uₕ), ᶜinterp(C123(uᵥ)))
-        )
+        (
+            dot(uₕ, CT12(uₕ) + CT12(ᶜinterp(uᵥ))) +
+            dot(ᶜinterp(uᵥ), CT3(uₕ)) +
+            ᶜinterp(dot(uᵥ, CT3(uᵥ)))
+        ) / 2
 end
 
 """