FINALLY reverse the k index

src/Operators/ops_regular_cartesian_grid.jl

@@ -10,10 +10,10 @@ using Oceananigans: AbstractGrid, RegularCartesianGrid
 @inline δy_e2f(g::RegularCartesianGrid, f, i, j, k) = @inbounds f[i, j+1, k] - f[i, j,   k]
 @inline δy_f2e(g::RegularCartesianGrid, f, i, j, k) = @inbounds f[i, j,   k] - f[i, j-1, k]
 
-@inline δz_c2f(g::RegularCartesianGrid, f, i, j, k) = @inbounds f[i, j, k-1] - f[i, j,   k]
-@inline δz_f2c(g::RegularCartesianGrid, f, i, j, k) = @inbounds f[i, j,   k] - f[i, j, k+1]
-@inline δz_e2f(g::RegularCartesianGrid, f, i, j, k) = @inbounds f[i, j,   k] - f[i, j, k+1]
-@inline δz_f2e(g::RegularCartesianGrid, f, i, j, k) = @inbounds f[i, j, k-1] - f[i, j,   k]
+@inline δz_c2f(g::RegularCartesianGrid, f, i, j, k) = @inbounds f[i, j,   k] - f[i, j, k-1]
+@inline δz_f2c(g::RegularCartesianGrid, f, i, j, k) = @inbounds f[i, j, k+1] - f[i, j,   k]
+@inline δz_e2f(g::RegularCartesianGrid, f, i, j, k) = @inbounds f[i, j, k+1] - f[i, j,   k]
+@inline δz_f2e(g::RegularCartesianGrid, f, i, j, k) = @inbounds f[i, j,   k] - f[i, j, k-1]
 
 @inline avgx_c2f(g::RegularCartesianGrid{T}, f, i, j, k) where T = @inbounds T(0.5) * (f[i,   j, k] + f[i-1,  j, k])
 @inline avgx_f2c(g::RegularCartesianGrid{T}, f, i, j, k) where T = @inbounds T(0.5) * (f[i+1, j, k] + f[i,    j, k])
@@ -23,9 +23,9 @@ using Oceananigans: AbstractGrid, RegularCartesianGrid
 @inline avgy_f2c(g::RegularCartesianGrid{T}, f, i, j, k) where T = @inbounds T(0.5) * (f[i, j+1, k] + f[i,    j, k])
 @inline avgy_f2e(g::RegularCartesianGrid{T}, f, i, j, k) where T = @inbounds T(0.5) * (f[i,   j, k] + f[i,  j-1, k])
 
-@inline avgz_c2f(g::RegularCartesianGrid{T}, f, i, j, k) where T = @inbounds T(0.5) * (f[i, j, k] + f[i, j, k-1])
-@inline avgz_f2c(g::RegularCartesianGrid{T}, f, i, j, k) where T = @inbounds T(0.5) * (f[i, j, k+1] + f[i, j, k])
-@inline avgz_f2e(g::RegularCartesianGrid{T}, f, i, j, k) where T = @inbounds T(0.5) * (f[i, j, k] + f[i, j, k-1])
+@inline avgz_c2f(g::RegularCartesianGrid{T}, f, i, j, k) where T = @inbounds T(0.5) * (f[i, j,   k] + f[i, j, k-1])
+@inline avgz_f2c(g::RegularCartesianGrid{T}, f, i, j, k) where T = @inbounds T(0.5) * (f[i, j, k+1] + f[i, j,   k])
+@inline avgz_f2e(g::RegularCartesianGrid{T}, f, i, j, k) where T = @inbounds T(0.5) * (f[i, j,   k] + f[i, j, k-1])
 
 @inline ∂x_p(i, j, k, grid::RegularCartesianGrid, p) = δx_c2f(grid, p, i, j, k) / grid.Δx
 @inline ∂y_p(i, j, k, grid::RegularCartesianGrid, p) = δy_c2f(grid, p, i, j, k) / grid.Δy
@@ -48,12 +48,12 @@ end
 end
 
 @inline function δz_f2c_ab̄ᶻ(g::RegularCartesianGrid, a, b, i, j, k)
-    @inbounds (a[i, j,   k] * avgz_c2f(g, b, i, j,   k) -
-               a[i, j, k+1] * avgz_c2f(g, b, i, j, k+1))
+    @inbounds (a[i, j, k+1] * avgz_c2f(g, b, i, j, k+1) -
+               a[i, j,   k] * avgz_c2f(g, b, i, j,   k))
 end
 
 @inline div_flux(g::RegularCartesianGrid, u, v, w, Q, i, j, k) =
-	(δx_f2c_ab̄ˣ(g, u, Q, i, j, k) / g.Δx) + (δy_f2c_ab̄ʸ(g, v, Q, i, j, k) / g.Δy) + (δz_f2c_ab̄ᶻ(g, w, Q, i, j, k) / g.Δz)
+    (δx_f2c_ab̄ˣ(g, u, Q, i, j, k) / g.Δx) + (δy_f2c_ab̄ʸ(g, v, Q, i, j, k) / g.Δy) + (δz_f2c_ab̄ᶻ(g, w, Q, i, j, k) / g.Δz)
 
 @inline δx_c2f_ūˣūˣ(g::RegularCartesianGrid, u, i, j, k) =
     avgx_f2c(g, u, i, j, k)^2 - avgx_f2c(g, u, i-1, j, k)^2
@@ -64,8 +64,8 @@ end
 end
 
 @inline function δz_e2f_w̄ˣūᶻ(g::RegularCartesianGrid, u, w, i, j, k)
-    @inbounds avgx_f2e(g, w, i, j,   k) * avgz_f2e(g, u, i, j,   k) -
-              avgx_f2e(g, w, i, j, k+1) * avgz_f2e(g, u, i, j, k+1)
+    @inbounds avgx_f2e(g, w, i, j, k+1) * avgz_f2e(g, u, i, j, k+1) -
+              avgx_f2e(g, w, i, j,   k) * avgz_f2e(g, u, i, j,   k)
 end
 
 @inline u∇u(g::RegularCartesianGrid, u, v, w, i, j, k) =
@@ -80,8 +80,8 @@ end
     avgy_f2c(g, v, i, j, k)^2 - avgy_f2c(g, v, i, j-1, k)^2
 
 @inline function δz_e2f_w̄ʸv̄ᶻ(g::RegularCartesianGrid, v, w, i, j, k)
-    @inbounds avgy_f2e(g, w, i, j,   k) * avgz_f2e(g, v, i, j,   k) -
-              avgy_f2e(g, w, i, j, k+1) * avgz_f2e(g, v, i, j, k+1)
+    @inbounds avgy_f2e(g, w, i, j, k+1) * avgz_f2e(g, v, i, j, k+1) -
+              avgy_f2e(g, w, i, j,   k) * avgz_f2e(g, v, i, j,   k)
 end
 
 @inline u∇v(g::RegularCartesianGrid, u, v, w, i, j, k) =
@@ -98,26 +98,25 @@ end
 end
 
 @inline δz_c2f_w̄ᶻw̄ᶻ(g::RegularCartesianGrid, w, i, j, k) =
-	avgz_f2c(g, w, i, j, k-1)^2 - avgz_f2c(g, w, i, j, k)^2
+    avgz_f2c(g, w, i, j, k)^2 - avgz_f2c(g, w, i, j, k-1)^2
 
 @inline u∇w(g::RegularCartesianGrid, u, v, w, i, j, k) =
     (δx_e2f_ūᶻw̄ˣ(g, u, w, i, j, k) / g.Δx) + (δy_e2f_v̄ᶻw̄ʸ(g, v, w, i, j, k) / g.Δy) + (δz_c2f_w̄ᶻw̄ᶻ(g, w, i, j, k) / g.Δz)
 
 @inline δx²_c2f2c(g::RegularCartesianGrid, f, i, j, k) = δx_c2f(g, f, i+1, j, k) - δx_c2f(g, f, i, j, k)
 @inline δy²_c2f2c(g::RegularCartesianGrid, f, i, j, k) = δy_c2f(g, f, i, j+1, k) - δy_c2f(g, f, i, j, k)
-@inline δz²_c2f2c(g::RegularCartesianGrid, f, i, j, k) = δz_c2f(g, f, i, j,   k) - δz_c2f(g, f, i, j, k+1)
+@inline δz²_c2f2c(g::RegularCartesianGrid, f, i, j, k) = δz_c2f(g, f, i, j, k+1) - δz_c2f(g, f, i, j, k)
 
 @inline κ∇²(g::RegularCartesianGrid, Q, κh, κv, i, j, k) =
     ((κh/g.Δx^2) * δx²_c2f2c(g, Q, i, j, k)) + ((κh/g.Δy^2) * δy²_c2f2c(g, Q, i, j, k)) + ((κv/g.Δz^2) * δz²_c2f2c(g, Q, i, j, k))
 
 @inline δx²_f2c2f(g::RegularCartesianGrid, f, i, j, k) = δx_f2c(g, f, i, j, k) - δx_f2c(g, f, i-1, j, k)
 @inline δy²_f2c2f(g::RegularCartesianGrid, f, i, j, k) = δy_f2c(g, f, i, j, k) - δy_f2c(g, f, i, j-1, k)
+@inline δz²_f2c2f(g::RegularCartesianGrid, f, i, j, k) = δz_f2c(g, f, i, j, k) - δz_f2c(g, f, i, j, k-1)
 
 @inline δx²_f2e2f(g::RegularCartesianGrid, f, i, j, k) = δx_f2e(g, f, i+1, j, k) - δx_f2e(g, f, i, j, k)
 @inline δy²_f2e2f(g::RegularCartesianGrid, f, i, j, k) = δy_f2e(g, f, i, j+1, k) - δy_f2e(g, f, i, j, k)
-
-@inline δz²_f2e2f(g::RegularCartesianGrid, f, i, j, k) = δz_f2e(g, f, i, j,   k) - δz_f2e(g, f, i, j, k+1)
-@inline δz²_f2c2f(g::RegularCartesianGrid, f, i, j, k) = δz_f2c(g, f, i, j, k-1) - δz_f2c(g, f, i, j, k)
+@inline δz²_f2e2f(g::RegularCartesianGrid, f, i, j, k) = δz_f2e(g, f, i, j, k+1) - δz_f2e(g, f, i, j, k)
 
 @inline 𝜈∇²u(g::RegularCartesianGrid, u, 𝜈h, 𝜈v, i, j, k) =
     ((𝜈h/g.Δx^2) * δx²_f2c2f(g, u, i, j, k)) + ((𝜈h/g.Δy^2) * δy²_f2e2f(g, u, i, j, k)) + ((𝜈v/g.Δz^2) * δz²_f2e2f(g, u, i, j, k))
@@ -129,6 +128,6 @@ end
     ((𝜈h/g.Δx^2) * δx²_f2e2f(g, w, i, j, k)) + ((𝜈h/g.Δy^2) * δy²_f2e2f(g, w, i, j, k)) + ((𝜈v/g.Δz^2) * δz²_f2c2f(g, w, i, j, k))
 
 @inline ∇²(g::RegularCartesianGrid, f, i, j, k) =
-	(δx²_c2f2c(g, f, i, j, k) / g.Δx^2) + (δy²_c2f2c(g, f, i, j, k) / g.Δy^2) + (δz²_c2f2c(g, f, i, j, k) / g.Δz^2)
+    (δx²_c2f2c(g, f, i, j, k) / g.Δx^2) + (δy²_c2f2c(g, f, i, j, k) / g.Δy^2) + (δz²_c2f2c(g, f, i, j, k) / g.Δz^2)
 
 @inline ∇h_u(i, j, k, grid, u, v) = δx_f2c(grid, u, i, j, k) / grid.Δx + δy_f2c(grid, v, i, j, k) / grid.Δy

src/TurbulenceClosures/closure_operators.jl

@@ -8,8 +8,8 @@
 @inline ∂y_aca(i, j, k, grid, v::AbstractArray) = @inbounds (v[i, j+1, k] - v[i, j, k]) / grid.Δy
 @inline ∂y_afa(i, j, k, grid, c::AbstractArray) = @inbounds (c[i, j, k] - c[i, j-1, k]) / grid.Δy
 
-@inline ∂z_aac(i, j, k, grid, w::AbstractArray) = @inbounds (w[i, j, k] - w[i, j, k+1]) / grid.Δz
-@inline ∂z_aaf(i, j, k, grid, c::AbstractArray) = @inbounds (c[i, j, k-1] - c[i, j, k]) / grid.Δz
+@inline ∂z_aac(i, j, k, grid, w::AbstractArray) = @inbounds (w[i, j, k+1] - w[i, j, k]) / grid.Δz
+@inline ∂z_aaf(i, j, k, grid, c::AbstractArray) = @inbounds (c[i, j, k] - c[i, j, k-1]) / grid.Δz
 
 #
 # Differentiation and interpolation operators for functions
@@ -108,7 +108,7 @@ Differentiate the function or callable object
 located at `aac` in `z`, across `aaf`.
 """
 @inline ∂z_aaf(i, j, k, grid::AbstractGrid, F::TF, args...) where TF<:Function =
-    (F(i, j, k-1, grid, args...) - F(i, j, k, grid, args...)) / grid.Δz
+    (F(i, j, k, grid, args...) - F(i, j, k-1, grid, args...)) / grid.Δz
 
 """
     ∂z_aac(i, j, k, grid, F, args...)
@@ -120,7 +120,7 @@ Differentiate the function or callable object
 located at `aaf` in `z`, across `aac`.
 """
 @inline ∂z_aac(i, j, k, grid, F::TF, args...) where TF<:Function =
-    (F(i, j, k, grid, args...) - F(i, j, k+1, grid, args...)) / grid.Δz
+    (F(i, j, k+1, grid, args...) - F(i, j, k, grid, args...)) / grid.Δz
 
 #####
 ##### Double differentiation

src/grids.jl

@@ -105,11 +105,11 @@ function RegularCartesianGrid(T, N, L)
 
     xC = Δx/2:Δx:Lx
     yC = Δy/2:Δy:Ly
-    zC = -Δz/2:-Δz:-Lz
+    zC = -Lz+Δz/2:Δz:-Δz/2
 
     xF = 0:Δx:Lx
     yF = 0:Δy:Ly
-    zF = 0:-Δz:-Lz
+    zF = -Lz:Δz:0
 
     RegularCartesianGrid{T, typeof(xC)}(Nx, Ny, Nz, Hx, Hy, Hz, Tx, Ty, Tz,
                                         Lx, Ly, Lz, Δx, Δy, Δz, Ax, Ay, Az, V,

src/halo_regions.jl

@@ -24,16 +24,16 @@ _fill_south_halo!(c, ::PBC, H, N) = @views @. c.parent[:, 1:H, :] = c.parent[:,
  _fill_north_halo!(c, ::PBC, H, N) = @views @. c.parent[:, N+H+1:N+2H, :] = c.parent[:, 1+H:2H, :]
 _fill_bottom_halo!(c, ::PBC, H, N) = @views @. c.parent[:, :, N+H+1:N+2H] = c.parent[:, :, 1+H:2H]
 
-# Recall that, by convention, the first grid point in an array with no penetration boundary
-# condition lies on the boundary, where as the final grid point lies in the domain.
+# Recall that, by convention, the last grid point (k=Nz) in an array with a no-penetration boundary
+# condition lies on the boundary, where as the first grid point (k=1) lies in the domain.
 
- _fill_west_halo!(c, ::NPBC, H, N) = @views @. c.parent[1:H+1, :, :] = 0
-_fill_south_halo!(c, ::NPBC, H, N) = @views @. c.parent[:, 1:H+1, :] = 0
-  _fill_top_halo!(c, ::NPBC, H, N) = @views @. c.parent[:, :, 1:H+1] = 0
+ _fill_west_halo!(c, ::NPBC, H, N) = @views @. c.parent[1:H, :, :] = 0
+_fill_south_halo!(c, ::NPBC, H, N) = @views @. c.parent[:, 1:H, :] = 0
+  _fill_top_halo!(c, ::NPBC, H, N) = @views @. c.parent[:, :, 1:H] = 0
 
-  _fill_east_halo!(c, ::NPBC, H, N) = @views @. c.parent[N+H+1:N+2H, :, :] = 0
- _fill_north_halo!(c, ::NPBC, H, N) = @views @. c.parent[:, N+H+1:N+2H, :] = 0
-_fill_bottom_halo!(c, ::NPBC, H, N) = @views @. c.parent[:, :, N+H+1:N+2H] = 0
+  _fill_east_halo!(c, ::NPBC, H, N) = @views @. c.parent[N+H:N+2H, :, :] = 0
+ _fill_north_halo!(c, ::NPBC, H, N) = @views @. c.parent[:, N+H:N+2H, :] = 0
+_fill_bottom_halo!(c, ::NPBC, H, N) = @views @. c.parent[:, :, N+H:N+2H] = 0
 
 # Generate functions that implement flux, periodic, and no-penetration boundary conditions
 sides = (:west, :east, :south, :north, :top, :bottom)
@@ -50,9 +50,9 @@ for (x, side) in zip(coords, sides)
     end
 end
 
-#####
-##### Halo filling for value and gradient boundary conditions
-#####
+####
+#### Halo filling for value and gradient boundary conditions
+####
 
 @inline linearly_extrapolate(c₀, ∇c, Δ) = c₀ + ∇c * Δ
 
@@ -98,9 +98,9 @@ function fill_bottom_halo!(c, bc::Union{VBC, GBC}, arch::AbstractArchitecture, g
     return
 end
 
-#####
-##### General halo filling functions
-#####
+####
+#### General halo filling functions
+####
 
 "Fill halo regions in x, y, and z for a given field."
 function fill_halo_regions!(c::AbstractArray, fieldbcs, arch, grid, args...)
@@ -140,6 +140,7 @@ function fill_halo_regions!(fields::NamedTuple, bcs::FieldBoundaryConditions, ar
     for field in fields
         fill_halo_regions!(field, bcs, arch, grid, args...)
     end
+    return nothing
 end
 
 fill_halo_regions!(::Nothing, args...) = nothing
@@ -163,12 +164,12 @@ function zero_halo_regions!(c::AbstractArray, grid)
      zero_north_halo!(c, grid.Hy, grid.Ny)
        zero_top_halo!(c, grid.Hz, grid.Nz)
     zero_bottom_halo!(c, grid.Hz, grid.Nz)
-    return
+    return nothing
 end
 
 function zero_halo_regions!(fields::Tuple, grid)
     for field in fields
         zero_halo_regions!(field, grid)
     end
-    return
+    return nothing
 end

src/time_steppers.jl

@@ -28,8 +28,8 @@ function time_step!(model, Nt, Δt; init_with_euler=true)
     for n in 1:Nt
         χ = ifelse(init_with_euler && n==1, FT(-0.5), model.timestepper.χ)
 
-        adams_bashforth_time_step!(model, model.architecture, model.grid, model.buoyancy, model.coriolis, 
-                                   model.closure, model.forcing, model.boundary_conditions, 
+        adams_bashforth_time_step!(model, model.architecture, model.grid, model.buoyancy, model.coriolis,
+                                   model.closure, model.forcing, model.boundary_conditions,
                                    U, Φ, p, K, RHS, Gⁿ,  G⁻, Δt, χ)
 
         [ time_to_run(model.clock, diag) && run_diagnostic(model, diag) for diag in values(model.diagnostics) ]
@@ -61,7 +61,7 @@ function adams_bashforth_time_step!(model, arch, grid, buoyancy, coriolis, closu
     fill_halo_regions!(p.pHY′, bcs.pressure, arch, grid)
 
     # Calculate tendency terms (minus non-hydrostatic pressure, which is updated in a pressure correction step):
-    calculate_interior_source_terms!(Gⁿ, arch, grid, coriolis, closure, U, Φ, p.pHY′, K, forcing, 
+    calculate_interior_source_terms!(Gⁿ, arch, grid, coriolis, closure, U, Φ, p.pHY′, K, forcing,
                                      model.parameters, model.clock.time)
     calculate_boundary_source_terms!(Gⁿ, arch, grid, bcs.solution, boundary_condition_args...)
 
@@ -128,10 +128,10 @@ the `buoyancy_perturbation` downwards:
 function update_hydrostatic_pressure!(pHY′, grid, buoyancy, Φ)
     @loop for j in (1:grid.Ny; (blockIdx().y - 1) * blockDim().y + threadIdx().y)
         @loop for i in (1:grid.Nx; (blockIdx().x - 1) * blockDim().x + threadIdx().x)
-            @inbounds pHY′[i, j, 1] = - ▶z_aaf(i, j, 1, grid, buoyancy_perturbation, buoyancy, Φ) * grid.Δz
-            @unroll for k in 2:grid.Nz
-                @inbounds pHY′[i, j, k] = 
-                    pHY′[i, j, k-1] - ▶z_aaf(i, j, k, grid, buoyancy_perturbation, buoyancy, Φ) * grid.Δz
+            @inbounds pHY′[i, j, grid.Nz] = - ▶z_aaf(i, j, grid.Nz, grid, buoyancy_perturbation, buoyancy, Φ) * grid.Δz
+            @unroll for k in grid.Nz-1 : -1 : 1
+                @inbounds pHY′[i, j, k] =
+                    pHY′[i, j, k+1] - ▶z_aaf(i, j, k+1, grid, buoyancy_perturbation, buoyancy, Φ) * grid.Δz
             end
         end
     end
@@ -212,19 +212,19 @@ end
 function calculate_interior_source_terms!(G, arch, grid, coriolis, closure, U, Φ, pHY′, K, F, parameters, time)
 
     Bx, By, Bz = floor(Int, grid.Nx/Tx), floor(Int, grid.Ny/Ty), grid.Nz  # Blocks in grid
-    @launch device(arch) threads=(Tx, Ty) blocks=(Bx, By, Bz) calculate_Gu!(G.Gu, grid, coriolis, closure, U, Φ, pHY′, 
+    @launch device(arch) threads=(Tx, Ty) blocks=(Bx, By, Bz) calculate_Gu!(G.Gu, grid, coriolis, closure, U, Φ, pHY′,
                                                                             K, F, parameters, time)
 
-    @launch device(arch) threads=(Tx, Ty) blocks=(Bx, By, Bz) calculate_Gv!(G.Gv, grid, coriolis, closure, U, Φ, pHY′, 
+    @launch device(arch) threads=(Tx, Ty) blocks=(Bx, By, Bz) calculate_Gv!(G.Gv, grid, coriolis, closure, U, Φ, pHY′,
                                                                             K, F, parameters, time)
 
-    @launch device(arch) threads=(Tx, Ty) blocks=(Bx, By, Bz) calculate_Gw!(G.Gw, grid, coriolis, closure, U, Φ, pHY′, 
+    @launch device(arch) threads=(Tx, Ty) blocks=(Bx, By, Bz) calculate_Gw!(G.Gw, grid, coriolis, closure, U, Φ, pHY′,
                                                                             K, F, parameters, time)
 
-    @launch device(arch) threads=(Tx, Ty) blocks=(Bx, By, Bz) calculate_GT!(G.GT, grid, closure, U, Φ, pHY′, 
+    @launch device(arch) threads=(Tx, Ty) blocks=(Bx, By, Bz) calculate_GT!(G.GT, grid, closure, U, Φ, pHY′,
                                                                             K, F, parameters, time)
 
-    @launch device(arch) threads=(Tx, Ty) blocks=(Bx, By, Bz) calculate_GS!(G.GS, grid, closure, U, Φ, pHY′, 
+    @launch device(arch) threads=(Tx, Ty) blocks=(Bx, By, Bz) calculate_GS!(G.GS, grid, closure, U, Φ, pHY′,
                                                                             K, F, parameters, time)
 end
 
@@ -392,7 +392,7 @@ function update_velocities_and_tracers!(grid::AbstractGrid, U, Φ, pNHS, Gⁿ, 
     @loop for k in (1:grid.Nz; (blockIdx().z - 1) * blockDim().z + threadIdx().z)
         @loop for j in (1:grid.Ny; (blockIdx().y - 1) * blockDim().y + threadIdx().y)
             @loop for i in (1:grid.Nx; (blockIdx().x - 1) * blockDim().x + threadIdx().x)
-                                            
+
                 @inbounds U.u[i, j, k] = U.u[i, j, k] + (Gⁿ.Gu[i, j, k] - ∂x_p(i, j, k, grid, pNHS)) * Δt
                 @inbounds U.v[i, j, k] = U.v[i, j, k] + (Gⁿ.Gv[i, j, k] - ∂y_p(i, j, k, grid, pNHS)) * Δt
                 @inbounds Φ.T[i, j, k] = Φ.T[i, j, k] + (Gⁿ.GT[i, j, k] * Δt)
@@ -410,9 +410,9 @@ Compute the vertical velocity w by integrating the continuity equation downwards
 function compute_w_from_continuity!(grid::AbstractGrid, U)
     @loop for j in (1:grid.Ny; (blockIdx().y - 1) * blockDim().y + threadIdx().y)
         @loop for i in (1:grid.Nx; (blockIdx().x - 1) * blockDim().x + threadIdx().x)
-            @inbounds U.w[i, j, 1] = 0
-            @unroll for k in 2:grid.Nz
-                @inbounds U.w[i, j, k] = U.w[i, j, k-1] + grid.Δz * ∇h_u(i, j, k-1, grid, U.u, U.v)
+            @inbounds U.w[i, j, grid.Nz] = 0
+            @unroll for k in grid.Nz-1 : -1 : 1
+                @inbounds U.w[i, j, k] = U.w[i, j, k+1] + grid.Δz * ∇h_u(i, j, k, grid, U.u, U.v)
             end
         end
     end

test/test_time_stepping.jl

@@ -47,12 +47,12 @@ function compute_w_from_continuity(arch, FT)
     fill_halo_regions!(w.data, bcs.w, arch, grid)
     velocity_div!(grid, u.data, v.data, w.data, div_u.data)
 
-    # Set div_u to zero at the bottom because the initial velocity field is not divergence-free
-    # so we end up some divergence at the bottom if we don't do this.
-    data(div_u)[:, :, end] .= zero(FT)
+    # Set div_u to zero at the top because the initial velocity field is not divergence-free
+    # so we end up some divergence at the top if we don't do this.
+    data(div_u)[:, :, Nz] .= zero(FT)
 
     min_div = minimum(data(div_u))
-    max_div = minimum(data(div_u))
+    max_div = maximum(data(div_u))
     sum_div = sum(data(div_u))
     abs_sum_div = sum(abs.(data(div_u)))
     @info "Velocity divergence after recomputing w ($arch, $FT): min=$min_div, max=$max_div, sum=$sum_div, abs_sum=$abs_sum_div"