make sure implicit operator works

src/Models/HydrostaticFreeSurfaceModels/hydrostatic_free_surface_rk3_step.jl

@@ -107,7 +107,7 @@ function rk3_substep_tracers!(tracers, model, Δt, γⁿ, ζⁿ)
         closure = model.closure
 
         launch!(architecture(grid), grid, :xyz,
-                _split_rk3_substep_field!, tracer_field, Δt, γⁿ, ζⁿ, Gⁿ, Ψ⁻)
+                _split_rk3_substep_tracer_field!, tracer_field, Δt, γⁿ, ζⁿ, Gⁿ, Ψ⁻)
 
         implicit_step!(tracer_field,
                        model.timestepper.implicit_solver,

src/TurbulenceClosures/abstract_scalar_diffusivity_closure.jl

@@ -259,7 +259,6 @@ end
 ##### Products of viscosity and stress, divergence, vorticity
 #####
 
-
 @inline κ_σᶠᶜᶜ(i, j, k, grid, closure, K, id, clock, fields, σᶠᶜᶜ, args...) = κᶠᶜᶜ(i, j, k, grid, closure, K, id, clock, fields) * σᶠᶜᶜ(i, j, k, grid, args...)
 @inline κ_σᶜᶠᶜ(i, j, k, grid, closure, K, id, clock, fields, σᶜᶠᶜ, args...) = κᶜᶠᶜ(i, j, k, grid, closure, K, id, clock, fields) * σᶜᶠᶜ(i, j, k, grid, args...)
 @inline κ_σᶜᶜᶠ(i, j, k, grid, closure, K, id, clock, fields, σᶜᶜᶠ, args...) = κᶜᶜᶠ(i, j, k, grid, closure, K, id, clock, fields) * σᶜᶜᶠ(i, j, k, grid, args...)

src/TurbulenceClosures/vertically_implicit_diffusion_solver.jl

@@ -45,13 +45,19 @@ implicit_diffusion_solver(::ExplicitTimeDiscretization, args...; kwargs...) = no
 const c = Center()
 const f = Face()
 
+# The vertical spacing used here is Δz for velocities and Δr for tracers, since the 
+# implicit solver operator is applied to the scaled tracer e₃θ instead of just θ
+
+@inline vertical_spacing(i, j, k, grid, ℓx, ℓy, ℓz) = Δz(i, j, k, grid, ℓx, ℓy, ℓz)
+@inline vertical_spacing(i, j, k, grid, ::Center, ::Center, ℓz) = Δr(i, j, k, grid, c, c, ℓz)
+
 # Tracers and horizontal velocities at cell centers in z
 @inline function ivd_upper_diagonal(i, j, k, grid, closure, K, id, ℓx, ℓy, ::Center, clock, Δt, κz)
     closure_ij = getclosure(i, j, closure)
     κᵏ⁺¹   = κz(i, j, k+1, grid, closure_ij, K, id, clock)
-    Δzᶜₖ   = Δz(i, j, k,   grid, ℓx, ℓy, c)
-    Δzᶠₖ₊₁ = Δz(i, j, k+1, grid, ℓx, ℓy, f)
-    du     = - Δt * κᵏ⁺¹ / (Δzᶜₖ * Δzᶠₖ₊₁)
+    Δzᶜₖ   = vertical_spacing(i, j, k,   grid, ℓx, ℓy, c)
+    Δzᶠₖ₊₁ = vertical_spacing(i, j, k+1, grid, ℓx, ℓy, f)
+    du     = - Δt * κᵏ⁺¹ / (Δrᶜₖ * Δrᶠₖ₊₁)
 
     # This conditional ensures the diagonal is correct
     return ifelse(k > grid.Nz-1, zero(grid), du)
@@ -61,9 +67,9 @@ end
     k = k′ + 1 # Shift index to match LinearAlgebra.Tridiagonal indexing convenction
     closure_ij = getclosure(i, j, closure)  
     κᵏ   = κz(i, j, k, grid, closure_ij, K, id, clock)
-    Δzᶜₖ = Δz(i, j, k, grid, ℓx, ℓy, c)
-    Δzᶠₖ = Δz(i, j, k, grid, ℓx, ℓy, f)
-    dl   = - Δt * κᵏ / (Δzᶜₖ * Δzᶠₖ)
+    Δrᶜₖ = vertical_spacing(i, j, k, grid, ℓx, ℓy, c)
+    Δrᶠₖ = vertical_spacing(i, j, k, grid, ℓx, ℓy, f)
+    dl   = - Δt * κᵏ / (Δrᶜₖ * Δrᶠₖ)
 
     # This conditional ensures the diagonal is correct: the lower diagonal does not
     # exist for k′ = 0. (Note we use LinearAlgebra.Tridiagonal indexing convention,
@@ -80,19 +86,19 @@ end
 @inline function ivd_upper_diagonal(i, j, k, grid, closure, K, id, ℓx, ℓy, ::Face, clock, Δt, νzᶜᶜᶜ) 
     closure_ij = getclosure(i, j, closure)  
     νᵏ = νzᶜᶜᶜ(i, j, k, grid, closure_ij, K, clock)
-    Δzᶜₖ = Δz(i, j, k, grid, ℓx, ℓy, c)
-    Δzᶠₖ = Δz(i, j, k, grid, ℓx, ℓy, f)
-    du   = - Δt * νᵏ / (Δzᶜₖ * Δzᶠₖ)
+    Δrᶜₖ = vertical_spacing(i, j, k, grid, ℓx, ℓy, c)
+    Δrᶠₖ = vertical_spacing(i, j, k, grid, ℓx, ℓy, f)
+    du   = - Δt * νᵏ / (Δrᶜₖ * Δrᶠₖ)
     return ifelse(k < 1, zero(grid), du)
 end
 
 @inline function ivd_lower_diagonal(i, j, k, grid, closure, K, id, ℓx, ℓy, ::Face, clock, Δt, νzᶜᶜᶜ)
     k′ = k + 2 # Shift to adjust for Tridiagonal indexing convention
     closure_ij = getclosure(i, j, closure)  
     νᵏ⁻¹   = νzᶜᶜᶜ(i, j, k′-1, grid, closure_ij, K, clock)
-    Δzᶜₖ   = Δz(i, j, k′,   grid, ℓx, ℓy, c)
-    Δzᶠₖ₋₁ = Δz(i, j, k′-1, grid, ℓx, ℓy, f)
-    dl     = - Δt * νᵏ⁻¹ / (Δzᶜₖ * Δzᶠₖ₋₁)
+    Δrᶜₖ   = vertical_spacing(i, j, k′,   grid, ℓx, ℓy, c)
+    Δrᶠₖ₋₁ = vertical_spacing(i, j, k′-1, grid, ℓx, ℓy, f)
+    dl     = - Δt * νᵏ⁻¹ / (Δrᶜₖ * Δrᶠₖ₋₁)
     return ifelse(k < 1, zero(grid), dl)
 end