Variance dissipation computation

src/Advection/adapt_advection_order.jl

@@ -39,7 +39,7 @@ function adapt_advection_order(advection, grid::AbstractGrid)
         @info "Using the advection scheme $(summary(new_advection.y)) in the y-direction because size(grid, 2) = $(size(grid, 2))"
     end
     if changed_z
-        @info "Using the advection scheme $(summary(new_advection.z)) in the x-direction because size(grid, 3) = $(size(grid, 3))"
+        @info "Using the advection scheme $(summary(new_advection.z)) in the z-direction because size(grid, 3) = $(size(grid, 3))"
     end
 
     return ifelse(changed_advection, new_advection, advection)

src/Models/Models.jl

@@ -188,4 +188,9 @@ default_nan_checker(::OnlyParticleTrackingModel) = nothing
 # applicable to both `NonhydrostaticModel` and  `HydrostaticFreeSurfaceModel`
 include("seawater_density.jl")
 
+# Implementation of the diagnostic for computing the dissipation rate
+include("VarianceDissipationComputation/VarianceDissipationComputation.jl")
+
+using .VarianceDissipationComputation
+
 end # module

src/Models/VarianceDissipationComputation/VarianceDissipationComputation.jl

@@ -0,0 +1,118 @@
+module VarianceDissipationComputation
+
+export VarianceDissipation, get_dissipation_fields
+
+using Oceananigans.Grids: architecture
+using Oceananigans.Utils
+using Oceananigans.TimeSteppers
+using Oceananigans.Fields
+using Oceananigans.Fields: Field, VelocityFields
+using Oceananigans.Operators
+using Oceananigans.BoundaryConditions
+using Oceananigans.TurbulenceClosures: viscosity,
+                                       diffusivity, 
+                                       ScalarDiffusivity, 
+                                       ScalarBiharmonicDiffusivity,
+                                       AbstractTurbulenceClosure,
+                                       HorizontalFormulation
+
+using Oceananigans.Advection: _advective_tracer_flux_x, 
+                              _advective_tracer_flux_y, 
+                              _advective_tracer_flux_z,
+                              horizontal_advection_U, 
+                              horizontal_advection_V 
+
+using Oceananigans.Operators: volume
+using KernelAbstractions: @kernel, @index
+
+struct VarianceDissipation{P, K, A, D, S, G}
+    advective_production :: P
+    diffusive_production :: K
+    advective_fluxes :: A
+    diffusive_fluxes :: D
+    previous_state :: S
+    gradient_squared :: G
+end
+
+include("dissipation_utils.jl")
+
+function VarianceDissipation(model; 
+                             tracers = propertynames(model.tracers), 
+                             include_vorticity = true)
+
+    if !(model.timestepper isa QuasiAdamsBashforth2TimeStepper)
+        throw(ArgumentError("DissipationComputation requires a QuasiAdamsBashforth2TimeStepper"))
+    end
+
+    tracers = tupleit(tracers)
+    diffusivities = model.diffusivity_fields
+    closure       = model.closure
+    grid          = model.grid
+
+    P    = NamedTuple{tracers}(tracer_fluxes(grid) for tracer in tracers)
+
+    K    = NamedTuple{tracers}(tracer_closure_dissipation(grid, diffusivities, closure, id) for id in eachindex(tracers))
+    Vⁿ   = NamedTuple{tracers}(tracer_closure_dissipation(grid, diffusivities, closure, id) for id in eachindex(tracers))
+    Vⁿ⁻¹ = NamedTuple{tracers}(tracer_closure_dissipation(grid, diffusivities, closure, id) for id in eachindex(tracers))    
+
+    K    = NamedTuple{tracers}(tracer_fluxes(grid) for tracer in tracers)
+    Fⁿ   = NamedTuple{tracers}(tracer_fluxes(grid) for tracer in tracers)
+    Fⁿ⁻¹ = NamedTuple{tracers}(tracer_fluxes(grid) for tracer in tracers)
+
+    Uⁿ⁻¹ = VelocityFields(grid)
+    Uⁿ   = VelocityFields(grid)
+
+    cⁿ⁻¹ =  NamedTuple{tracers}(CenterField(grid) for tracer in tracers)
+
+    if include_vorticity
+        Fζⁿ   = vorticity_fluxes(grid)
+        Fζⁿ⁻¹ = vorticity_fluxes(grid)
+        Pζ    = vorticity_fluxes(grid)
+        ζⁿ⁻¹  = Field{Face, Face, Center}(grid)
+
+        P    = merge(P,    (; ζ = Pζ))
+        Fⁿ   = merge(Fⁿ,   (; ζ = Fζⁿ))
+        Fⁿ⁻¹ = merge(Fⁿ⁻¹, (; ζ = Fζⁿ⁻¹))
+        cⁿ⁻¹ = merge(cⁿ⁻¹, (; ζ = ζⁿ⁻¹))
+
+        Kζ    = enstrophy_closure_dissipation(grid, diffusivities, closure)
+        Vζⁿ   = enstrophy_closure_dissipation(grid, diffusivities, closure)
+        Vζⁿ⁻¹ = enstrophy_closure_dissipation(grid, diffusivities, closure)
+
+        K    = merge(K,    (; ζ = Kζ))
+        Vⁿ   = merge(Vⁿ,   (; ζ = Vζⁿ))
+        Vⁿ⁻¹ = merge(Vⁿ⁻¹, (; ζ = Vζⁿ⁻¹))
+    end
+
+    previous_state   = merge(cⁿ⁻¹, (; Uⁿ⁻¹, Uⁿ))
+    advective_fluxes = (; Fⁿ, Fⁿ⁻¹)
+    diffusive_fluxes = (; Vⁿ, Vⁿ⁻¹)
+
+    gradients = deepcopy(P)
+
+    return VarianceDissipation(P, K, advective_fluxes, diffusive_fluxes, previous_state, gradients)
+end
+
+# Function to call in a callback
+# Note: This works only if the callback is called with an IterationInterval(1), if not the
+# previous fluxes and velocities will not be correct
+# TODO: make sure that the correct velocities and fluxes are used even if 
+# the callback is not called with an IterationInterval(1)
+function (ϵ::VarianceDissipation)(simulation)
+
+    # We first assemble values for Pⁿ⁻¹
+    assemble_dissipation!(simulation, ϵ)
+
+    # Then we update the fluxes to be used in the next time step
+    update_fluxes!(simulation, ϵ)
+
+    return nothing
+end
+
+include("get_dissipation_fields.jl")
+include("update_fluxes.jl")
+include("advective_fluxes.jl")
+include("diffusive_fluxes.jl")
+include("assemble_dissipation.jl")
+
+end

src/Models/VarianceDissipationComputation/advective_dissipation.jl

@@ -0,0 +1,61 @@
+# TODO: This is only for AB2, figure out how to generalize this for other timesteppers for example RK3
+@kernel function _assemble_advective_tracer_dissipation!(P, grid, χ, Fⁿ, Fⁿ⁻¹, Uⁿ⁺¹, Uⁿ, Uⁿ⁻¹, cⁿ⁺¹, cⁿ)
+    i, j, k = @index(Global, NTuple)
+
+    δˣc★ = δxᶠᶜᶜ(i, j, k, grid, c★, cⁿ⁺¹, cⁿ)
+    δˣc² = δxᶠᶜᶜ(i, j, k, grid, c², cⁿ⁺¹, cⁿ)
+
+    δʸc★ = δyᶜᶠᶜ(i, j, k, grid, c★, cⁿ⁺¹, cⁿ)
+    δʸc² = δyᶜᶠᶜ(i, j, k, grid, c², cⁿ⁺¹, cⁿ)
+
+    δᶻc★ = δzᶜᶜᶠ(i, j, k, grid, c★, cⁿ⁺¹, cⁿ)
+    δᶻc² = δzᶜᶜᶠ(i, j, k, grid, c², cⁿ⁺¹, cⁿ)
+
+    C₁  = convert(eltype(grid), 1.5 + χ)
+    C₂  = convert(eltype(grid), 0.5 + χ)
+
+    @inbounds begin
+        u₁ = C₁ *   Uⁿ.u[i, j, k] 
+        u₂ = C₂ * Uⁿ⁻¹.u[i, j, k] 
+        v₁ = C₁ *   Uⁿ.v[i, j, k] 
+        v₂ = C₂ * Uⁿ⁻¹.v[i, j, k] 
+        w₁ = C₁ *   Uⁿ.w[i, j, k] 
+        w₂ = C₂ * Uⁿ⁻¹.w[i, j, k] 
+
+        fx₁ = C₁ * Fⁿ.x[i, j, k]
+        fx₂ = C₂ * Fⁿ⁻¹.x[i, j, k]
+        fy₁ = C₁ * Fⁿ.y[i, j, k]
+        fy₂ = C₂ * Fⁿ⁻¹.y[i, j, k]
+        fz₁ = C₁ * Fⁿ.z[i, j, k]
+        fz₂ = C₂ * Fⁿ⁻¹.z[i, j, k]
+
+        P.x[i, j, k] = 2 * δˣc★ * (fx₁ - fx₂) - δˣc² * (u₁ - u₂)
+        P.y[i, j, k] = 2 * δʸc★ * (fy₁ - fy₂) - δʸc² * (v₁ - v₂)
+        P.z[i, j, k] = 2 * δᶻc★ * (fz₁ - fz₂) - δᶻc² * (w₁ - w₂)
+    end
+end
+
+@kernel function _assemble_advective_vorticity_dissipation!(P, grid, χ, Fⁿ, Fⁿ⁻¹, Uⁿ⁺¹, Uⁿ, Uⁿ⁻¹, c, ζⁿ)
+    i, j, k = @index(Global, NTuple)
+
+    δˣζ★ = δxᶠᶜᶜ(i, j, k, grid, ζ★, Uⁿ⁺¹.u, Uⁿ⁺¹.v, ζⁿ)
+    δˣζ² = δxᶠᶜᶜ(i, j, k, grid, ζ², Uⁿ⁺¹.u, Uⁿ⁺¹.v, ζⁿ)
+
+    δʸζ★ = δyᶜᶠᶜ(i, j, k, grid, ζ★, Uⁿ⁺¹.u, Uⁿ⁺¹.v, ζⁿ)
+    δʸζ² = δyᶜᶠᶜ(i, j, k, grid, ζ², Uⁿ⁺¹.u, Uⁿ⁺¹.v, ζⁿ)
+
+    @inbounds begin
+        u₁ = C₁ * ℑxyᶜᶠᵃ(i, j, k, grid, Δy_qᶠᶜᶜ, Uⁿ.u)   / Δyᶠᶜᶜ(i, j, k, grid)
+        u₂ = C₂ * ℑxyᶜᶠᵃ(i, j, k, grid, Δy_qᶠᶜᶜ, Uⁿ⁻¹.u) / Δyᶠᶜᶜ(i, j, k, grid)
+        v₁ = C₁ * ℑxyᶠᶜᵃ(i, j, k, grid, Δx_qᶜᶠᶜ, Uⁿ.v)   / Δxᶜᶠᶜ(i, j, k, grid)
+        v₂ = C₂ * ℑxyᶠᶜᵃ(i, j, k, grid, Δx_qᶜᶠᶜ, Uⁿ⁻¹.v) / Δxᶜᶠᶜ(i, j, k, grid)
+
+        fx₁ = C₁ * Fⁿ.x[i, j, k]
+        fx₂ = C₂ * Fⁿ⁻¹.x[i, j, k]
+        fy₁ = C₁ * Fⁿ.y[i, j, k]
+        fy₂ = C₂ * Fⁿ⁻¹.y[i, j, k]
+
+        P.x[i, j, k] = 2 * δˣζ★ * (fx₁ - fx₂) - δˣζ² * (u₁ - u₂)
+        P.y[i, j, k] = 2 * δʸζ★ * (fy₁ - fy₂) - δʸζ² * (v₁ - v₂)
+    end
+end

src/Models/VarianceDissipationComputation/advective_fluxes.jl

@@ -0,0 +1,42 @@
+@kernel function _update_advective_tracer_fluxes!(Gⁿ, Fⁿ, Fⁿ⁻¹, cⁿ⁻¹, grid, advection, U, c)
+    i, j, k = @index(Global, NTuple)
+    u, v, w = U
+
+    @inbounds begin
+        # Save previous advective fluxes
+        Fⁿ⁻¹.x[i, j, k] = Fⁿ.x[i, j, k]
+        Fⁿ⁻¹.y[i, j, k] = Fⁿ.y[i, j, k]
+        Fⁿ⁻¹.z[i, j, k] = Fⁿ.z[i, j, k]
+
+        cⁿ⁻¹[i, j, k] = c[i, j, k]
+
+        # Calculate new advective fluxes
+        Fⁿ.x[i, j, k] = _advective_tracer_flux_x(i, j, k, grid, advection, u, c) 
+        Fⁿ.y[i, j, k] = _advective_tracer_flux_y(i, j, k, grid, advection, v, c) 
+        Fⁿ.z[i, j, k] = _advective_tracer_flux_z(i, j, k, grid, advection, w, c) 
+
+        Gⁿ.x[i, j, k] = Axᶠᶜᶜ(i, j, k, grid) * δxᶠᶜᶜ(i, j, k, grid, c)^2 / Δxᶠᶜᶜ(i, j, k, grid)
+        Gⁿ.y[i, j, k] = Ayᶜᶠᶜ(i, j, k, grid) * δyᶜᶠᶜ(i, j, k, grid, c)^2 / Δyᶜᶠᶜ(i, j, k, grid)
+        Gⁿ.z[i, j, k] = Azᶜᶜᶠ(i, j, k, grid) * δzᶜᶜᶠ(i, j, k, grid, c)^2 / Δzᶜᶜᶠ(i, j, k, grid)
+    end
+end
+
+@kernel function _update_advective_vorticity_fluxes!(Gⁿ, Fⁿ, Fⁿ⁻¹, ζⁿ⁻¹, grid, advection, U, c)
+    i, j, k = @index(Global, NTuple)
+    u, v, w = U
+
+    @inbounds begin
+        # Save previous advective fluxes
+        Fⁿ⁻¹.x[i, j, k] = Fⁿ.x[i, j, k]
+        Fⁿ⁻¹.y[i, j, k] = Fⁿ.y[i, j, k]
+
+        ζⁿ⁻¹[i, j, k] = ζ₃ᶠᶠᶜ(i, j, k, grid, U.u, U.v)
+
+        # Calculate new advective fluxes
+        Fⁿ.x[i, j, k] =   horizontal_advection_V(i, j, k, grid, advection, u, ζ) * Axᶜᶠᶜ(i, j, k, grid) 
+        Fⁿ.y[i, j, k] = - horizontal_advection_U(i, j, k, grid, advection, v, ζ) * Ayᶠᶜᶜ(i, j, k, grid)
+
+        Gⁿ.x[i, j, k] = Axᶜᶠᶜ(i, j, k, grid) * δxᶜᶠᶜ(i, j, k, grid, ζ₃ᶠᶠᶜ, U.u)^2 / Δxᶜᶠᶜ(i, j, k, grid)
+        Gⁿ.y[i, j, k] = Ayᶠᶜᶜ(i, j, k, grid) * δyᶠᶜᶜ(i, j, k, grid, ζ₃ᶠᶠᶜ, U.v)^2 / Δyᶠᶜᶜ(i, j, k, grid)
+    end
+end

src/Models/VarianceDissipationComputation/assemble_dissipation.jl

@@ -0,0 +1,61 @@
+function assemble_dissipation!(simulation, dissipation)
+    model = simulation.model
+
+    for tracer_name in keys(dissipation.advective_production)
+        assemble_dissipation!(dissipation, model, tracer_name)
+    end
+
+    return nothing
+end
+
+@inline c★(i, j, k, grid, cⁿ⁺¹, cⁿ) = @inbounds (cⁿ⁺¹[i, j, k] + cⁿ[i, j, k]) / 2
+@inline c²(i, j, k, grid, cⁿ⁺¹, cⁿ) = @inbounds (cⁿ⁺¹[i, j, k] * cⁿ[i, j, k])
+
+@inline ζ★(i, j, k, grid, uⁿ⁺¹, vⁿ⁺¹, ζⁿ) = @inbounds (ζ₃ᶠᶠᶜ(i, j, k, grid, uⁿ⁺¹, vⁿ⁺¹) + ζⁿ[i, j, k]) / 2
+@inline ζ²(i, j, k, grid, uⁿ⁺¹, vⁿ⁺¹, ζⁿ) = @inbounds (ζ₃ᶠᶠᶜ(i, j, k, grid, uⁿ⁺¹, vⁿ⁺¹) * ζⁿ[i, j, k])
+
+function assemble_dissipation!(dissipation, model, tracer_name::Symbol)
+
+
+    arch = architecture(grid)
+    χ = simulation.model.timestepper.χ
+
+    # General velocities
+    Uⁿ⁺¹ = model.velocities
+    Uⁿ   = dissipation.previous_state.Uⁿ
+    Uⁿ⁻¹ = dissipation.previous_state.Uⁿ⁻¹
+
+    cⁿ⁺¹ = tracer_symbol == :ζ ? nothing : model.tracers[tracer_name]
+    cⁿ   = dissipation.previous_state[tracer_name]
+
+
+    _assemble_advective_dissipation! = assemble_advective_dissipation_kernel(Val(tracer_name))
+    _assemble_diffusive_dissipation! = assemble_diffusive_dissipation_kernel(Val(tracer_name))
+
+    ####
+    #### Assemble the advective dissipation
+    ####
+
+    P    = dissipation.advective_production[tracer_name]
+    Fⁿ   = dissipation.advective_fluxes.Fⁿ[tracer_name]
+    Fⁿ⁻¹ = dissipation.advective_fluxes.Fⁿ⁻¹[tracer_name]
+
+    launch!(arch, grid, :xyz, _assemble_advective_dissipation!, P, grid, χ, Fⁿ, Fⁿ⁻¹, Uⁿ⁺¹, Uⁿ, Uⁿ⁻¹, cⁿ⁺¹, cⁿ)
+
+    ####
+    #### Assemble the diffusive dissipation
+    #### 
+
+    K    = dissipation.diffusive_production[tracer_name]
+    Vⁿ   = dissipation.advective_fluxes.Vⁿ[tracer_name]
+    Vⁿ⁻¹ = dissipation.advective_fluxes.Vⁿ⁻¹[tracer_name]
+
+    launch!(arch, grid, params, _assemble_diffusive_dissipation!, K, grid, χ, Vⁿ, Vⁿ⁻¹, Uⁿ⁺¹, cⁿ⁺¹, cⁿ)
+
+    return nothing
+end
+
+assemble_advective_dissipation_kernel(val_tracer_name) = _assemble_advective_tracer_dissipation!
+assemble_advective_dissipation_kernel(::Val{:ζ})       = _assemble_advective_vorticity_dissipation!
+assemble_diffusive_dissipation_kernel(val_tracer_name) = _assemble_diffusive_tracer_dissipation!
+assemble_diffusive_dissipation_kernel(::Val{:ζ})       = _assemble_diffusive_vorticity_dissipation!

src/Models/VarianceDissipationComputation/diffusive_dissipation.jl

@@ -0,0 +1,49 @@
+@kernel function _assemble_diffusive_tracer_dissipation!(K, grid, χ, Vⁿ, Vⁿ⁻¹, Uⁿ⁺¹, cⁿ⁺¹, cⁿ)
+    i, j, k = @index(Global, NTuple)
+    compute_diffusive_dissipation!(K, i, j, k, grid, Vⁿ, Vⁿ⁻¹, χ, cⁿ⁺¹, cⁿ)
+end
+
+@inline function compute_diffusive_tracer_dissipation!(K::Tuple, i, j, k, grid, Vⁿ, Vⁿ⁻¹, χ, cⁿ⁺¹, cⁿ)
+    for n in eachindex(K)
+        compute_diffusive_dissipation!(K[n], i, j, k, grid, Vⁿ[n], Vⁿ⁻¹[n], χ, cⁿ⁺¹, cⁿ)
+    end
+end
+
+@inline function compute_diffusive_tracer_dissipation!(K, i, j, k, grid, Vⁿ, Vⁿ⁻¹, χ, cⁿ⁺¹, cⁿ)
+    C₁  = convert(eltype(grid), 1.5 + χ)
+    C₂  = convert(eltype(grid), 0.5 + χ)
+
+    δˣc★ = δxᶠᶜᶜ(i, j, k, grid, c★, cⁿ⁺¹, cⁿ)    
+    δʸc★ = δyᶜᶠᶜ(i, j, k, grid, c★, cⁿ⁺¹, cⁿ)
+    δᶻc★ = δzᶜᶜᶠ(i, j, k, grid, c★, cⁿ⁺¹, cⁿ)
+
+    @inbounds begin
+        K.x[i, j, k] = 2 * δˣc★ * (C₁ * Vⁿ.x[i, j, k] - C₂ * Vⁿ⁻¹.x[i, j, k])
+        K.y[i, j, k] = 2 * δʸc★ * (C₁ * Vⁿ.y[i, j, k] - C₂ * Vⁿ⁻¹.y[i, j, k])
+        K.z[i, j, k] = 2 * δᶻc★ * (C₁ * Vⁿ.z[i, j, k] - C₂ * Vⁿ⁻¹.z[i, j, k])
+    end
+end
+
+@kernel function _assemble_diffusive_vorticity_dissipation!(K, grid, χ, Vⁿ, Vⁿ⁻¹, Uⁿ⁺¹, cⁿ⁺¹, ζⁿ)
+    i, j, k = @index(Global, NTuple)
+    compute_diffusive_vorticity_dissipation!(K, i, j, k, grid, Vⁿ, Vⁿ⁻¹, χ, Uⁿ⁺¹, ζⁿ)
+end
+
+@inline function compute_diffusive_vorticity_dissipation!(K::Tuple, i, j, k, grid, Vⁿ, Vⁿ⁻¹, χ, Uⁿ⁺¹, ζⁿ)
+    for n in eachindex(K)
+        compute_diffusive_vorticity_dissipation!(K[n], i, j, k, grid, Vⁿ[n], Vⁿ⁻¹[n], χ, Uⁿ⁺¹, ζⁿ)
+    end
+end
+
+@inline function compute_diffusive_vorticity_dissipation!(K, i, j, k, grid, Vⁿ, Vⁿ⁻¹, χ, Uⁿ⁺¹, ζⁿ)
+    C₁  = convert(eltype(grid), 1.5 + χ)
+    C₂  = convert(eltype(grid), 0.5 + χ)
+
+    δˣζ★ = δxᶠᶜᶜ(i, j, k, grid, ζ★, Uⁿ⁺¹.u, Uⁿ⁺¹.v, ζⁿ)
+    δʸζ★ = δyᶜᶠᶜ(i, j, k, grid, ζ★, Uⁿ⁺¹.u, Uⁿ⁺¹.v, ζⁿ)
+
+    @inbounds begin
+        K.x[i, j, k] = 2 * δˣζ★ * (C₁ * Vⁿ.x[i, j, k] - C₂ * Vⁿ⁻¹.x[i, j, k])
+        K.y[i, j, k] = 2 * δʸζ★ * (C₁ * Vⁿ.y[i, j, k] - C₂ * Vⁿ⁻¹.y[i, j, k])
+    end
+end

src/Models/VarianceDissipationComputation/diffusive_fluxes.jl

@@ -0,0 +1,40 @@
+
+@kenrel function _update_diffusive_tracer_fluxes!(Vⁿ, Vⁿ⁻¹, grid, closure, diffusivity, bouyancy, c, tracer_id, clk, model_fields) 
+    i, j, k = @index(Global, NTuple)
+    compute_diffusive_tracer_fluxes!(Vⁿ, Vⁿ⁻¹, i, j, k, grid, closure, diffusivity, bouyancy, c, tracer_id, clk, model_fields)
+end
+
+@inline function compute_diffusive_tracer_fluxes!(Vⁿ, Vⁿ⁻¹, i, j, k, grid, closure::Tuple, K, args...) 
+    for n in eachindex(closure)
+        compute_diffusive_tracer_fluxes!(Vⁿ[n], Vⁿ⁻¹[n], i, j, k, grid, closure[n], K[n], args...)
+    end
+end
+
+@inline function compute_diffusive_tracer_fluxes!(Vⁿ, Vⁿ⁻¹, i, j, k, grid, clo, K, b, c, c_id, clk, fields)
+    Vⁿ⁻¹.x[i, j, k] = Vⁿ.x[i, j, k]
+    Vⁿ⁻¹.y[i, j, k] = Vⁿ.y[i, j, k]
+    Vⁿ⁻¹.z[i, j, k] = Vⁿ.z[i, j, k]
+
+    Vⁿ.x[i, j, k] = _diffusive_tracer_flux_x(i, j, k, grid, clo, K, Val(c_id), c, clk, fields, b) * Axᶠᶜᶜ(i, j, k, grid)
+    Vⁿ.y[i, j, k] = _diffusive_tracer_flux_y(i, j, k, grid, clo, K, Val(c_id), c, clk, fields, b) * Ayᶜᶠᶜ(i, j, k, grid)
+    Vⁿ.z[i, j, k] = _diffusive_tracer_flux_z(i, j, k, grid, clo, K, Val(c_id), c, clk, fields, b) * Azᶜᶜᶠ(i, j, k, grid)
+end
+
+@kenrel function _update_diffusive_vorticity_fluxes!(Vⁿ, Vⁿ⁻¹, grid, closure, diffusivity, bouyancy, c, tracer_id, clk, model_fields)
+    i, j, k = @index(Global, NTuple)
+    compute_diffusive_vorticity_fluxes!(Vⁿ, Vⁿ⁻¹, i, j, k, grid, closure, diffusivity, bouyancy, clk, model_fields)
+end
+
+@inline function compute_diffusive_vorticity_fluxes!(Vⁿ, Vⁿ⁻¹, i, j, k, grid, closure::Tuple, K, args...) 
+    for n in eachindex(closure)
+        compute_diffusive_vorticity_fluxes!(Vⁿ[n], Vⁿ⁻¹[n], i, j, k, grid, closure[n], K[n], args...)
+    end
+end
+
+@inline function compute_diffusive_vorticity_fluxes!(Vⁿ, Vⁿ⁻¹, i, j, k, grid, clo, K, b, clk, fields) 
+    Vⁿ⁻¹.x[i, j, k] = Vⁿ.x[i, j, k]
+    Vⁿ⁻¹.y[i, j, k] = Vⁿ.y[i, j, k]
+
+    Vⁿ.x[i, j, k] =   ∂ⱼ_τ₂ⱼ(i, j, k, grid, clo, K, clk, fields, b) * Axᶜᶠᶜ(i, j, k, grid)
+    Vⁿ.y[i, j, k] = - ∂ⱼ_τ₁ⱼ(i, j, k, grid, clo, K, clk, fields, b) * Ayᶠᶜᶜ(i, j, k, grid)
+end