diff --git a/src/Oceananigans.jl b/src/Oceananigans.jl
index ae11c5176e..d957c9bac8 100644
--- a/src/Oceananigans.jl
+++ b/src/Oceananigans.jl
@@ -49,7 +49,7 @@ export
 
     # Model output writers
     Checkpointer, restore_from_checkpoint, read_output,
-    JLD2OutputWriter, FieldOutput, FieldOutputs,
+    JLD2OutputWriter, NetCDFOutputWriter, FieldOutput, FieldOutputs,
     write_grid, NetCDFOutputWriter,
 
     # Model diagnostics
diff --git a/src/Operators/ops_regular_cartesian_grid.jl b/src/Operators/ops_regular_cartesian_grid.jl
index cac8c33371..6edb90e9f5 100644
--- a/src/Operators/ops_regular_cartesian_grid.jl
+++ b/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
diff --git a/src/TurbulenceClosures/closure_operators.jl b/src/TurbulenceClosures/closure_operators.jl
index f9869b5e0b..7e417cc7dc 100644
--- a/src/TurbulenceClosures/closure_operators.jl
+++ b/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
diff --git a/src/boundary_conditions.jl b/src/boundary_conditions.jl
index 941e5d3591..8239290129 100644
--- a/src/boundary_conditions.jl
+++ b/src/boundary_conditions.jl
@@ -139,20 +139,17 @@ const CBC = CoordinateBoundaryConditions
 PeriodicBCs() = CBC(PeriodicBC(), PeriodicBC())
 
 # Here we overload setproperty! and getproperty to permit users to call
-# the 'right' and 'left' bcs in the z-direction 'bottom' and 'top'.
-#
-# Note that 'right' technically corresponds to face point N+1. Thus
-# the fact that right == bottom is associated with the reverse z-indexing
-# convention. With ordinary indexing, right == top.
+# the 'left' and 'right' bcs in the z-direction 'bottom' and 'top'.
+# Note that 'right' technically corresponds to face point N+1.
 Base.setproperty!(cbc::CBC, side::Symbol, bc) = setbc!(cbc, Val(side), bc)
 setbc!(cbc::CBC, ::Val{S}, bc) where S = setfield!(cbc, S, bc)
-setbc!(cbc::CBC, ::Val{:bottom}, bc) = setfield!(cbc, :right, bc)
-setbc!(cbc::CBC, ::Val{:top}, bc) = setfield!(cbc, :left, bc)
+setbc!(cbc::CBC, ::Val{:bottom}, bc) = setfield!(cbc, :left, bc)
+setbc!(cbc::CBC, ::Val{:top}, bc) = setfield!(cbc, :right, bc)
 
 Base.getproperty(cbc::CBC, side::Symbol) = getbc(cbc, Val(side))
 getbc(cbc::CBC, ::Val{S}) where S = getfield(cbc, S)
-getbc(cbc::CBC, ::Val{:bottom}) = getfield(cbc, :right)
-getbc(cbc::CBC, ::Val{:top}) = getfield(cbc, :left)
+getbc(cbc::CBC, ::Val{:bottom}) = getfield(cbc, :left)
+getbc(cbc::CBC, ::Val{:top}) = getfield(cbc, :right)
 
 #####
 ##### Boundary conditions for Fields
@@ -189,7 +186,7 @@ function HorizontallyPeriodicBCs(;    top = BoundaryCondition(Flux, nothing),
 
     x = PeriodicBCs()
     y = PeriodicBCs()
-    z = CoordinateBoundaryConditions(top, bottom)
+    z = CoordinateBoundaryConditions(bottom, top)
 
     return FieldBoundaryConditions(x, y, z)
 end
@@ -214,7 +211,7 @@ function ChannelBCs(;  north = BoundaryCondition(Flux, nothing),
 
     x = PeriodicBCs()
     y = CoordinateBoundaryConditions(south, north)
-    z = CoordinateBoundaryConditions(top, bottom)
+    z = CoordinateBoundaryConditions(bottom, top)
 
     return FieldBoundaryConditions(x, y, z)
 end
@@ -228,12 +225,12 @@ default_tracer_bcs(tracers, solution_bcs) = DefaultTracerBoundaryConditions(solu
 """
     SolutionBoundaryConditions(tracers, proposal_bcs)
 
-Construct a `NamedTuple` of `FieldBoundaryConditions` for a model with 
-fields `u`, `v`, `w`, and `tracers` from the proposal boundary conditions 
-`proposal_bcs`, which must contain the boundary conditions on `u`, `v`, and `w` 
+Construct a `NamedTuple` of `FieldBoundaryConditions` for a model with
+fields `u`, `v`, `w`, and `tracers` from the proposal boundary conditions
+`proposal_bcs`, which must contain the boundary conditions on `u`, `v`, and `w`
 and may contain some or all of the boundary conditions on `tracers`.
 """
-SolutionBoundaryConditions(tracers, proposal_bcs) = 
+SolutionBoundaryConditions(tracers, proposal_bcs) =
     with_tracers(tracers, proposal_bcs, default_tracer_bcs, with_velocities=true)
 
 """
@@ -282,7 +279,7 @@ const BoundaryConditions = HorizontallyPeriodicSolutionBCs
 #####
 ##### Tracer, tendency and pressure boundary condition "translators":
 #####
-#####   * Default boundary conditions on tracers are periodic or no flux and 
+#####   * Default boundary conditions on tracers are periodic or no flux and
 #####     can be derived from boundary conditions on any field
 #####
 #####   * Boundary conditions on tendency terms are
@@ -296,7 +293,7 @@ const BoundaryConditions = HorizontallyPeriodicSolutionBCs
 DefaultTracerBC(::BC) = BoundaryCondition(Flux, nothing)
 DefaultTracerBC(::PBC) = PeriodicBC()
 
-DefaultTracerCoordinateBCs(bcs) = 
+DefaultTracerCoordinateBCs(bcs) =
     CoordinateBoundaryConditions(DefaultTracerBC(bcs.left), DefaultTracerBC(bcs.right))
 
 DefaultTracerBoundaryConditions(field_bcs) =
@@ -310,7 +307,7 @@ TendencyBC(::NPBC) = NoPenetrationBC()
 TendencyCoordinateBCs(bcs) =
     CoordinateBoundaryConditions(TendencyBC(bcs.left), TendencyBC(bcs.right))
 
-TendencyFieldBCs(field_bcs) = 
+TendencyFieldBCs(field_bcs) =
     FieldBoundaryConditions(Tuple(TendencyCoordinateBCs(bcs) for bcs in field_bcs))
 
 TendenciesBoundaryConditions(solution_bcs) =
@@ -365,7 +362,7 @@ the source term `Gc` at the top and bottom.
 """
 function apply_z_bcs!(Gc, arch, grid, top_bc, bottom_bc, args...)
     @launch device(arch) config=launch_config(grid, 2) _apply_z_bcs!(Gc, grid, top_bc, bottom_bc, args...)
-    return
+    return nothing
 end
 
 # Fall back functions for boundary conditions that are not of type Flux.
@@ -390,7 +387,7 @@ If `top_bc.condition` is a function, the function must have the signature
     `top_bc.condition(i, j, grid, boundary_condition_args...)`
 """
 @inline apply_z_top_bc!(Gc, top_flux::BC{<:Flux}, i, j, grid, args...) =
-    @inbounds Gc[i, j, 1] -= getbc(top_flux, i, j, grid, args...) / grid.Δz
+    @inbounds Gc[i, j, grid.Nz] -= getbc(top_flux, i, j, grid, args...) / grid.Δz
 
 """
     apply_z_bottom_bc!(Gc, bottom_flux::BC{<:Flux}, i, j, grid, args...)
@@ -406,18 +403,18 @@ If `bottom_bc.condition` is a function, the function must have the signature
     `bottom_bc.condition(i, j, grid, boundary_condition_args...)`
 """
 @inline apply_z_bottom_bc!(Gc, bottom_flux::BC{<:Flux}, i, j, grid, args...) =
-    @inbounds Gc[i, j, grid.Nz] += getbc(bottom_flux, i, j, grid, args...) / grid.Δz
+    @inbounds Gc[i, j, 1] += getbc(bottom_flux, i, j, grid, args...) / grid.Δz
 
 """
     _apply_z_bcs!(Gc, grid, top_bc, bottom_bc, args...)
 
 Apply a top and/or bottom boundary condition to variable `c`.
 """
-function _apply_z_bcs!(Gc, grid, top_bc, bottom_bc, args...)
+function _apply_z_bcs!(Gc, grid, bottom_bc, top_bc, args...)
     @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)
-               apply_z_top_bc!(Gc, top_bc,    i, j, grid, args...)
             apply_z_bottom_bc!(Gc, bottom_bc, i, j, grid, args...)
+               apply_z_top_bc!(Gc, top_bc,    i, j, grid, args...)
         end
     end
 end
diff --git a/src/diagnostics.jl b/src/diagnostics.jl
index e4f2f2ba40..b76e95532a 100644
--- a/src/diagnostics.jl
+++ b/src/diagnostics.jl
@@ -79,7 +79,7 @@ function run_diagnostic(model, havg::HorizontalAverage{NTuple{1}})
     zero_halo_regions!(havg.fields[1], model.grid)
     sum!(havg.profile, havg.fields[1])
     normalize_horizontal_sum!(havg, model.grid)
-    return
+    return nothing
 end
 
 function run_diagnostic(model::Model, havg::HorizontalAverage)
@@ -93,7 +93,7 @@ function run_diagnostic(model::Model, havg::HorizontalAverage)
     sum!(havg.profile, tmp)
     normalize_horizontal_sum!(havg, model.grid)
 
-    return
+    return nothing
 end
 
 function (havg::HorizontalAverage{F, Nothing})(model) where F
diff --git a/src/grids.jl b/src/grids.jl
index 4b265e75fc..291762689f 100644
--- a/src/grids.jl
+++ b/src/grids.jl
@@ -106,13 +106,13 @@ function RegularCartesianGrid(T, N, L)
 
     V = Δx*Δy*Δz
 
-    xC = Δx/2:Δx:Lx
-    yC = Δy/2:Δy:Ly
-    zC = -Δz/2:-Δz:-Lz
+    xC = range(Δx/2, Lx-Δx/2; length=Nx)
+    yC = range(Δy/2, Ly-Δy/2; length=Ny)
+    zC = range(-Lz+Δz/2, -Δz/2; length=Nz)
 
-    xF = 0:Δx:Lx
-    yF = 0:Δy:Ly
-    zF = 0:-Δz:-Lz
+    xF = range(0, Lx; length=Nx+1)
+    yF = range(0, Ly; length=Ny+1)
+    zF = range(-Lz, 0; length=Nz+1)
 
     RegularCartesianGrid{T, typeof(xC)}(Nx, Ny, Nz, Hx, Hy, Hz, Tx, Ty, Tz,
                                         Lx, Ly, Lz, Δx, Δy, Δz, Ax, Ay, Az, V,
diff --git a/src/halo_regions.jl b/src/halo_regions.jl
index 51a6f45e5d..4271bd7711 100644
--- a/src/halo_regions.jl
+++ b/src/halo_regions.jl
@@ -7,33 +7,33 @@
 # ranges are used to reference the data copied into halos, as this produces views of the correct
 # dimension (eg size = (1, Ny, Nz) for the west halos).
 
- _fill_west_halo!(c, ::FBC, H, N) = @views @. c.parent[1:H, :, :] = c.parent[1+H:1+H,  :, :]
-_fill_south_halo!(c, ::FBC, H, N) = @views @. c.parent[:, 1:H, :] = c.parent[:, 1+H:1+H,  :]
-  _fill_top_halo!(c, ::FBC, H, N) = @views @. c.parent[:, :, 1:H] = c.parent[:, :,  1+H:1+H]
+  _fill_west_halo!(c, ::FBC, H, N) = @views @. c.parent[1:H, :, :] = c.parent[1+H:1+H,  :, :]
+ _fill_south_halo!(c, ::FBC, H, N) = @views @. c.parent[:, 1:H, :] = c.parent[:, 1+H:1+H,  :]
+_fill_bottom_halo!(c, ::FBC, H, N) = @views @. c.parent[:, :, 1:H] = c.parent[:, :,  1+H:1+H]
 
-  _fill_east_halo!(c, ::FBC, H, N) = @views @. c.parent[N+H+1:N+2H, :, :] = c.parent[N+H:N+H, :, :]
- _fill_north_halo!(c, ::FBC, H, N) = @views @. c.parent[:, N+H+1:N+2H, :] = c.parent[:, N+H:N+H, :]
-_fill_bottom_halo!(c, ::FBC, H, N) = @views @. c.parent[:, :, N+H+1:N+2H] = c.parent[:, :, N+H:N+H]
+ _fill_east_halo!(c, ::FBC, H, N) = @views @. c.parent[N+H+1:N+2H, :, :] = c.parent[N+H:N+H, :, :]
+_fill_north_halo!(c, ::FBC, H, N) = @views @. c.parent[:, N+H+1:N+2H, :] = c.parent[:, N+H:N+H, :]
+  _fill_top_halo!(c, ::FBC, H, N) = @views @. c.parent[:, :, N+H+1:N+2H] = c.parent[:, :, N+H:N+H]
 
 # Periodic boundary conditions
- _fill_west_halo!(c, ::PBC, H, N) = @views @. c.parent[1:H, :, :] = c.parent[N+1:N+H, :, :]
-_fill_south_halo!(c, ::PBC, H, N) = @views @. c.parent[:, 1:H, :] = c.parent[:, N+1:N+H, :]
-  _fill_top_halo!(c, ::PBC, H, N) = @views @. c.parent[:, :, 1:H] = c.parent[:, :, N+1:N+H]
+  _fill_west_halo!(c, ::PBC, H, N) = @views @. c.parent[1:H, :, :] = c.parent[N+1:N+H, :, :]
+ _fill_south_halo!(c, ::PBC, H, N) = @views @. c.parent[:, 1:H, :] = c.parent[:, N+1:N+H, :]
+_fill_bottom_halo!(c, ::PBC, H, N) = @views @. c.parent[:, :, 1:H] = c.parent[:, :, N+1:N+H]
 
-  _fill_east_halo!(c, ::PBC, H, N) = @views @. c.parent[N+H+1:N+2H, :, :] = c.parent[1+H:2H, :, :]
- _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]
+ _fill_east_halo!(c, ::PBC, H, N) = @views @. c.parent[N+H+1:N+2H, :, :] = c.parent[1+H:2H, :, :]
+_fill_north_halo!(c, ::PBC, H, N) = @views @. c.parent[:, N+H+1:N+2H, :] = c.parent[:, 1+H:2H, :]
+  _fill_top_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 first grid point (k=1) in an array with a no-penetration boundary
+# condition lies on the boundary, where as the last grid point (k=Nz) 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:1+H, :, :] = 0
+ _fill_south_halo!(c, ::NPBC, H, N) = @views @. c.parent[:, 1:1+H, :] = 0
+_fill_bottom_halo!(c, ::NPBC, H, N) = @views @. c.parent[:, :, 1: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+1:N+2H, :, :] = 0
+_fill_north_halo!(c, ::NPBC, H, N) = @views @. c.parent[:, N+H+1:N+2H, :] = 0
+  _fill_top_halo!(c, ::NPBC, H, N) = @views @. c.parent[:, :, N+H+1:N+2H] = 0
 
 # Generate functions that implement flux, periodic, and no-penetration boundary conditions
 sides = (:west, :east, :south, :north, :top, :bottom)
@@ -50,71 +50,69 @@ 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 * Δ
 
-@inline top_gradient(bc::GBC, c¹, Δ, i, j, args...) = getbc(bc, i, j, args...)
-@inline bottom_gradient(bc::GBC, cᴺ, Δ, i, j, args...) = getbc(bc, i, j, args...)
+@inline top_gradient(bc::GBC, cᴺ, Δ, i, j, args...) = getbc(bc, i, j, args...)
+@inline bottom_gradient(bc::GBC, c¹, Δ, i, j, args...) = getbc(bc, i, j, args...)
 
-@inline top_gradient(bc::VBC, c¹, Δ, i, j, args...) =    ( getbc(bc, i, j, args...) - c¹ ) / (Δ/2)
-@inline bottom_gradient(bc::VBC, cᴺ, Δ, i, j, args...) = ( cᴺ - getbc(bc, i, j, args...) ) / (Δ/2)
+@inline top_gradient(bc::VBC, cᴺ, Δ, i, j, args...) =    ( getbc(bc, i, j, args...) - cᴺ ) / (Δ/2)
+@inline bottom_gradient(bc::VBC, c¹, Δ, i, j, args...) = ( c¹ - getbc(bc, i, j, args...) ) / (Δ/2)
 
 function _fill_top_halo!(c, bc::Union{VBC, GBC}, grid, args...)
     @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 ∇c = top_gradient(bc, c[i, j, 1], grid.Δz, i, j, grid, args...)
-            @unroll for k in (1-grid.Hz):0
-                Δ = (1-k) * grid.Δz
-                @inbounds c[i, j, k] = linearly_extrapolate(c[i, j, 1], ∇c, Δ)
+            @inbounds ∇c = top_gradient(bc, c[i, j, grid.Nz], grid.Δz, i, j, grid, args...)
+            @unroll for k in (grid.Nz + 1) : (grid.Nz + grid.Hz)
+                Δ = (k - grid.Nz) * grid.Δz
+                @inbounds c[i, j, k] = linearly_extrapolate(c[i, j, grid.Nz], ∇c, Δ)
             end
         end
     end
-    return
+    return nothing
 end
 
 function _fill_bottom_halo!(c, bc::Union{VBC, GBC}, grid, args...)
     @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 ∇c = bottom_gradient(bc, c[i, j, grid.Nz], grid.Δz, i, j, grid, args...)
-            @unroll for k in grid.Nz+1:grid.Nz+grid.Hz
-                Δ = (grid.Nz-k) * grid.Δz # separation between bottom grid cell and halo is negative
-                @inbounds c[i, j, k] = linearly_extrapolate(c[i, j, grid.Nz], ∇c, Δ)
+            @inbounds ∇c = bottom_gradient(bc, c[i, j, 1], grid.Δz, i, j, grid, args...)
+            @unroll for k in (1 - grid.Hz):0
+                Δ = (k - 1) * grid.Δz  # separation between bottom grid cell and halo is negative
+                @inbounds c[i, j, k] = linearly_extrapolate(c[i, j, 1], ∇c, Δ)
             end
         end
     end
-    return
+    return nothing
 end
 
 function fill_top_halo!(c, bc::Union{VBC, GBC}, arch, grid, args...)
     @launch device(arch) config=launch_config(grid, 2) _fill_top_halo!(c, bc, grid, args...)
-    return
+    return nothing
 end
 
 function fill_bottom_halo!(c, bc::Union{VBC, GBC}, arch::AbstractArchitecture, grid::AbstractGrid, args...)
     @launch device(arch) config=launch_config(grid, 2) _fill_bottom_halo!(c, bc, grid, args...)
-    return
+    return nothing
 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...)
-
       fill_west_halo!(c, fieldbcs.x.left,  arch, grid, args...)
       fill_east_halo!(c, fieldbcs.x.right, arch, grid, args...)
 
      fill_south_halo!(c, fieldbcs.y.left,  arch, grid, args...)
      fill_north_halo!(c, fieldbcs.y.right, arch, grid, args...)
 
-       fill_top_halo!(c, fieldbcs.z.left,  arch, grid, args...)
-    fill_bottom_halo!(c, fieldbcs.z.right, arch, grid, args...)
-
-    return
+     fill_bottom_halo!(c, fieldbcs.z.bottom, arch, grid, args...)
+        fill_top_halo!(c, fieldbcs.z.top,    arch, grid, args...)
+    return nothing
 end
 
 """
@@ -127,7 +125,7 @@ function fill_halo_regions!(fields::NamedTuple, bcs, arch, grid, args...)
     for (field, fieldbcs) in zip(fields, bcs)
         fill_halo_regions!(field, fieldbcs, arch, grid, args...)
     end
-    return
+    return nothing
 end
 
 """
@@ -140,6 +138,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
@@ -148,13 +147,13 @@ fill_halo_regions!(::Nothing, args...) = nothing
 #### Halo zeroing functions
 ####
 
- zero_west_halo!(c, H, N) = @views @. c[1:H, :, :] = 0
-zero_south_halo!(c, H, N) = @views @. c[:, 1:H, :] = 0
-  zero_top_halo!(c, H, N) = @views @. c[:, :, 1:H] = 0
+  zero_west_halo!(c, H, N) = @views @. c[1:H, :, :] = 0
+ zero_south_halo!(c, H, N) = @views @. c[:, 1:H, :] = 0
+zero_bottom_halo!(c, H, N) = @views @. c[:, :, 1:H] = 0
 
-  zero_east_halo!(c, H, N) = @views @. c[N+H+1:N+2H, :, :] = 0
- zero_north_halo!(c, H, N) = @views @. c[:, N+H+1:N+2H, :] = 0
-zero_bottom_halo!(c, H, N) = @views @. c[:, :, N+H+1:N+2H] = 0
+ zero_east_halo!(c, H, N) = @views @. c[N+H+1:N+2H, :, :] = 0
+zero_north_halo!(c, H, N) = @views @. c[:, N+H+1:N+2H, :] = 0
+  zero_top_halo!(c, H, N) = @views @. c[:, :, N+H+1:N+2H] = 0
 
 function zero_halo_regions!(c::AbstractArray, grid)
       zero_west_halo!(c, grid.Hx, grid.Nx)
@@ -163,12 +162,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
diff --git a/src/output_writers.jl b/src/output_writers.jl
index 7faed74563..d75299bfab 100644
--- a/src/output_writers.jl
+++ b/src/output_writers.jl
@@ -621,4 +621,3 @@ function restore_from_checkpoint(filepath; kwargs=Dict())
 
     return model
 end
-
diff --git a/src/time_steppers.jl b/src/time_steppers.jl
index 64749ea691..690cbbc162 100644
--- a/src/time_steppers.jl
+++ b/src/time_steppers.jl
@@ -59,7 +59,7 @@ function adams_bashforth_time_step!(model, U, C, P, K, Gⁿ, G⁻, Δt, χ)
     fill_halo_regions!(P.pHY′, model.boundary_conditions.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, model.coriolis, model.closure, U, C, P.pHY′, K, model.forcing, 
+    calculate_interior_source_terms!(Gⁿ, arch, grid, model.coriolis, model.closure, U, C, P.pHY′, K, model.forcing,
                                      model.parameters, model.clock.time)
     calculate_boundary_source_terms!(Gⁿ, model.boundary_conditions.solution, arch, grid, bc_args...)
 
@@ -68,7 +68,7 @@ function adams_bashforth_time_step!(model, U, C, P, K, Gⁿ, G⁻, Δt, χ)
 
     # Start pressure correction substep with a pressure solve:
     fill_halo_regions!(Gⁿ[1:3], model.boundary_conditions.tendency[1:3], arch, grid)
-    @launch device(arch) config=launch_config(grid, 3) calculate_poisson_right_hand_side!(model.poisson_solver.storage, arch, grid, 
+    @launch device(arch) config=launch_config(grid, 3) calculate_poisson_right_hand_side!(model.poisson_solver.storage, arch, grid,
                                                                                           model.poisson_solver.bcs, Δt, U, Gⁿ)
     solve_for_pressure!(P.pNHS, arch, grid, model.poisson_solver, model.poisson_solver.storage)
     fill_halo_regions!(P.pNHS, model.boundary_conditions.pressure, arch, grid)
@@ -156,22 +156,22 @@ function calculate_interior_source_terms!(G, arch, grid, coriolis, closure, U, C
 
     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.u, grid, coriolis, closure, U, C, K, F, 
+    @launch device(arch) threads=(Tx, Ty) blocks=(Bx, By, Bz) calculate_Gu!(G.u, grid, coriolis, closure, U, C, K, F,
                                                                             pHY′, parameters, time)
 
-    @launch device(arch) threads=(Tx, Ty) blocks=(Bx, By, Bz) calculate_Gv!(G.v, grid, coriolis, closure, U, C, K, F, 
+    @launch device(arch) threads=(Tx, Ty) blocks=(Bx, By, Bz) calculate_Gv!(G.v, grid, coriolis, closure, U, C, K, F,
                                                                             pHY′, parameters, time)
 
-    @launch device(arch) threads=(Tx, Ty) blocks=(Bx, By, Bz) calculate_Gw!(G.w, grid, coriolis, closure, U, C, K, F, 
+    @launch device(arch) threads=(Tx, Ty) blocks=(Bx, By, Bz) calculate_Gw!(G.w, grid, coriolis, closure, U, C, K, F,
                                                                             parameters, time)
 
     for tracer_index in 1:length(C)
         @inbounds Gc = G[tracer_index+3]
         @inbounds Fc = F[tracer_index+3]
         @inbounds  c = C[tracer_index]
-        @launch device(arch) threads=(Tx, Ty) blocks=(Bx, By, Bz) calculate_Gc!(Gc, grid, closure, c, Val(tracer_index), 
+        @launch device(arch) threads=(Tx, Ty) blocks=(Bx, By, Bz) calculate_Gc!(Gc, grid, closure, c, Val(tracer_index),
                                                                                 U, C, K, Fc, parameters, time)
-                                                                                
+
     end
 
     return nothing
@@ -182,12 +182,12 @@ function calculate_boundary_source_terms!(Gⁿ, bcs, arch, grid, args...)
 
     # Velocity fields
     for (i, ubcs) in enumerate(bcs[1:3])
-        apply_z_bcs!(Gⁿ[i], arch, grid, ubcs.z.left, ubcs.z.right, args...)
+        apply_z_bcs!(Gⁿ[i], arch, grid, ubcs.z.bottom, ubcs.z.top, args...)
     end
 
     # Tracer fields
     for (i, cbcs) in enumerate(bcs[4:end])
-        apply_z_bcs!(Gⁿ[i+3], arch, grid, cbcs.z.left, cbcs.z.right, args...)
+        apply_z_bcs!(Gⁿ[i+3], arch, grid, cbcs.z.bottom, cbcs.z.top, args...)
     end
 
     return nothing
@@ -201,14 +201,14 @@ end
 function solve_for_pressure!(pressure, ::CPU, grid, poisson_solver, ϕ)
     solve_poisson_3d!(poisson_solver, grid)
     view(pressure, 1:grid.Nx, 1:grid.Ny, 1:grid.Nz) .= real.(ϕ)
-    return nothing 
+    return nothing
 end
 
 "Solve the Poisson equation for non-hydrostatic pressure on the GPU."
 function solve_for_pressure!(pressure, ::GPU, grid, poisson_solver, ϕ)
     solve_poisson_3d!(poisson_solver, grid)
     @launch device(GPU()) config=launch_config(grid, 3) idct_permute!(pressure, grid, poisson_solver.bcs, ϕ)
-    return nothing 
+    return nothing
 end
 
 """
@@ -220,10 +220,10 @@ the `buoyancy_perturbation` downwards:
 function update_hydrostatic_pressure!(pHY′, grid, buoyancy, C)
     @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, C) * 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, C) * grid.Δz
+            @inbounds pHY′[i, j, grid.Nz] = - ▶z_aaf(i, j, grid.Nz, grid, buoyancy_perturbation, buoyancy, C) * 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, C) * grid.Δz
             end
         end
     end
@@ -454,7 +454,7 @@ function adams_bashforth_update_tracer_source_term!(Gcⁿ, grid::AbstractGrid{FT
 end
 
 """
-Evaluate the right-hand-side terms for velocity fields and tracer fields 
+Evaluate the right-hand-side terms for velocity fields and tracer fields
 at time step n+½ using a weighted 2nd-order Adams-Bashforth method.
 """
 function adams_bashforth_update_source_terms!(Gⁿ, arch, grid, χ, G⁻)
@@ -522,20 +522,19 @@ function update_solution!(U, C, arch, grid, G, pNHS, Δt)
 end
 
 """
-Compute the vertical velocity w by integrating the continuity equation downwards
+Compute the vertical velocity w by integrating the continuity equation from the bottom upwards
 
     `w^{n+1} = -∫ [∂/∂x (u^{n+1}) + ∂/∂y (v^{n+1})] dz`
 """
 function compute_w_from_continuity!(U, grid)
     @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
+            # U.w[i, j, 0] = 0 is enforced via halo regions.
             @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, k] = U.w[i, j, k-1] - grid.Δz * ∇h_u(i, j, k-1, grid, U.u, U.v)
             end
         end
     end
 
     return nothing
 end
-
diff --git a/src/utils.jl b/src/utils.jl
index b641943d8b..f92423ef78 100644
--- a/src/utils.jl
+++ b/src/utils.jl
@@ -258,7 +258,7 @@ Ensure that frequency and interval are not both `nothing`.
 """
 function validate_interval(frequency, interval)
     isnothing(frequency) && isnothing(interval) && @error "Must specify a frequency or interval!"
-    return
+    return nothing
 end
 
 has_interval(obj) = :interval in propertynames(obj) && obj.interval != nothing
diff --git a/test/regression_tests/ocean_large_eddy_simulation_regression_test.jl b/test/regression_tests/ocean_large_eddy_simulation_regression_test.jl
index b95fac4d0f..27209fc739 100644
--- a/test/regression_tests/ocean_large_eddy_simulation_regression_test.jl
+++ b/test/regression_tests/ocean_large_eddy_simulation_regression_test.jl
@@ -1,36 +1,40 @@
 using Oceananigans.TurbulenceClosures: VerstappenAnisotropicMinimumDissipation
 
-interiordata(a, grid) = view(a, grid.Hx+1:grid.Nx+grid.Hx, 
-                                grid.Hy+1:grid.Ny+grid.Hy, 
-                                grid.Hz+1:grid.Nz+grid.Hz, 
+interiordata(a, grid) = view(a, grid.Hx+1:grid.Nx+grid.Hx,
+                                grid.Hy+1:grid.Ny+grid.Hy,
+                                grid.Hz+1:grid.Nz+grid.Hz,
                             )
 
+# Temporary until k index is reversed
+const Nx, Ny, Nz = 16, 16, 16
+const Hx, Hy, Hz = 1, 1, 1
+
 function get_fields_from_checkpoint(filename)
 
     file = jldopen(filename)
 
     solution = (
-        u = file["velocities/u"],
-        v = file["velocities/v"],
-        w = file["velocities/w"],
-        T = file["tracers/T"],
-        S = file["tracers/S"]
+        u = reverse(file["velocities/u"]; dims=3),
+        v = reverse(file["velocities/v"]; dims=3),
+        w = cat(zeros(Nx+2Hx, Ny+2Hy), reverse(file["velocities/w"][:, :, 2:Nz+2Hz]; dims=3); dims=3),
+        T = reverse(file["tracers/T"]; dims=3),
+        S = reverse(file["tracers/S"]; dims=3)
     )
 
     Gⁿ = (
-        u = file["timestepper/Gⁿ/Gu"],
-        v = file["timestepper/Gⁿ/Gv"],
-        w = file["timestepper/Gⁿ/Gw"],
-        T = file["timestepper/Gⁿ/GT"],
-        S = file["timestepper/Gⁿ/GS"]
+        u = reverse(file["timestepper/Gⁿ/Gu"]; dims=3),
+        v = reverse(file["timestepper/Gⁿ/Gv"]; dims=3),
+        w = cat(zeros(Nx+2Hx, Ny+2Hy), reverse(file["timestepper/Gⁿ/Gw"][:, :, 2:Nz+2Hz]; dims=3); dims=3),
+        T = reverse(file["timestepper/Gⁿ/GT"]; dims=3),
+        S = reverse(file["timestepper/Gⁿ/GS"]; dims=3)
     )
-    
+
     G⁻ = (
-        u = file["timestepper/G⁻/Gu"],
-        v = file["timestepper/G⁻/Gv"],
-        w = file["timestepper/G⁻/Gw"],
-        T = file["timestepper/G⁻/GT"],
-        S = file["timestepper/G⁻/GS"]
+        u = reverse(file["timestepper/G⁻/Gu"]; dims=3),
+        v = reverse(file["timestepper/G⁻/Gv"]; dims=3),
+        w = cat(zeros(Nx+2Hx, Ny+2Hy), reverse(file["timestepper/G⁻/Gw"][:, :, 2:Nz+2Hz]; dims=3); dims=3),
+        T = reverse(file["timestepper/G⁻/GT"]; dims=3),
+        S = reverse(file["timestepper/G⁻/GS"]; dims=3)
     )
 
     close(file)
@@ -54,12 +58,12 @@ function run_ocean_large_eddy_simulation_regression_test(arch, closure)
     u_bcs = HorizontallyPeriodicBCs(    top = BoundaryCondition(Flux, Qᵘ)       )
     T_bcs = HorizontallyPeriodicBCs(    top = BoundaryCondition(Flux, Qᵀ),
                                      bottom = BoundaryCondition(Gradient, ∂T∂z) )
-    S_bcs = HorizontallyPeriodicBCs(    top = BoundaryCondition(Flux, 5e-8)   )
+    S_bcs = HorizontallyPeriodicBCs(    top = BoundaryCondition(Flux, 5e-8)     )
 
     # Model instantiation
     model = Model(
-             architecture = arch,  
-                     grid = RegularCartesianGrid(N=(16, 16, 16), L=(16, 16, 16)), 
+             architecture = arch,
+                     grid = RegularCartesianGrid(N=(16, 16, 16), L=(16, 16, 16)),
                  coriolis = FPlane(f=1e-4),
                  buoyancy = SeawaterBuoyancy(equation_of_state=LinearEquationOfState(α=2e-4, β=8e-4)),
                   closure = closure,
@@ -82,7 +86,7 @@ function run_ocean_large_eddy_simulation_regression_test(arch, closure)
                                                           prefix = name * "_",
                                                              dir = joinpath(dirname(@__FILE__), "data")
                                                       )
-                                                       
+
     time_step!(model, 2test_steps, Δt)
     pop!(model.output_writers, :checkpointer)
     =#
@@ -90,7 +94,7 @@ function run_ocean_large_eddy_simulation_regression_test(arch, closure)
     initial_filename = joinpath(dirname(@__FILE__), "data", name * "_$spinup_steps.jld2")
 
     solution₀, Gⁿ₀, G⁻₀ = get_fields_from_checkpoint(initial_filename)
-    
+
     model.velocities.u.data.parent .= ArrayType(solution₀.u)
     model.velocities.v.data.parent .= ArrayType(solution₀.v)
     model.velocities.w.data.parent .= ArrayType(solution₀.w)
diff --git a/test/regression_tests/rayleigh_benard_regression_test.jl b/test/regression_tests/rayleigh_benard_regression_test.jl
index e34a2495f2..d913c605b3 100644
--- a/test/regression_tests/rayleigh_benard_regression_test.jl
+++ b/test/regression_tests/rayleigh_benard_regression_test.jl
@@ -1,3 +1,4 @@
+
 function run_rayleigh_benard_regression_test(arch)
 
     #####
@@ -26,14 +27,16 @@ function run_rayleigh_benard_regression_test(arch)
     c★(x, z) = exp(4z) * sin(2π/Lx * x)
     Fc(i, j, k, grid, time, U, C, params) = 1/10 * (c★(grid.xC[i], grid.zC[k]) - C.c[i, j, k])
 
+    bbcs = HorizontallyPeriodicBCs(   top = BoundaryCondition(Value, 0.0),
+                                   bottom = BoundaryCondition(Value, Δb))
+
     model = Model(
                architecture = arch,
                        grid = RegularCartesianGrid(N=(Nx, Ny, Nz), L=(Lx, Ly, Lz)),
                     closure = ConstantIsotropicDiffusivity(ν=ν, κ=κ),
                     tracers = (:b, :c),
-                   buoyancy = BuoyancyTracer(), 
-        boundary_conditions = BoundaryConditions(b=HorizontallyPeriodicBCs(top=BoundaryCondition(Value, 0.0), 
-                                                                           bottom=BoundaryCondition(Value, Δb))),
+                   buoyancy = BuoyancyTracer(),
+        boundary_conditions = BoundaryConditions(b=bbcs),
                     forcing = ModelForcing(c=Fc)
     )
 
@@ -49,9 +52,8 @@ function run_rayleigh_benard_regression_test(arch)
     outputfields = Dict(:U=>output_U, :Φ=>output_Φ, :G=>output_G)
 
     prefix = "data_rayleigh_benard_regression"
-    outputwriter = JLD2OutputWriter(model, outputfields; dir=joinpath(dirname(@__FILE__), "data"), 
+    outputwriter = JLD2OutputWriter(model, outputfields; dir=joinpath(dirname(@__FILE__), "data"),
                                     prefix=prefix, frequency=test_steps, including=[])
-                                   
 
     #####
     ##### Initial condition and spinup steps for creating regression test data
@@ -80,6 +82,22 @@ function run_rayleigh_benard_regression_test(arch)
     T₀, S₀ = get_output_tuple(outputwriter, spinup_steps, :Φ)
     Gu, Gv, Gw, GT, GS = get_output_tuple(outputwriter, spinup_steps, :G)
 
+    u₀ = reverse(u₀; dims=3)
+    v₀ = reverse(v₀; dims=3)
+    T₀ = reverse(T₀; dims=3)
+    S₀ = reverse(S₀; dims=3)
+
+    w₀ = reverse(w₀[:, :, 2:Nz]; dims=3)
+    w₀ = cat(zeros(Nx, Ny), w₀; dims=3)
+
+    Gu = reverse(Gu; dims=3)
+    Gv = reverse(Gv; dims=3)
+    GT = reverse(GT; dims=3)
+    GS = reverse(GS; dims=3)
+
+    Gw = reverse(Gw[:, :, 2:Nz]; dims=3)
+    Gw = cat(zeros(Nx, Ny), Gw; dims=3)
+
     data(model.velocities.u) .= ArrayType(u₀)
     data(model.velocities.v) .= ArrayType(v₀)
     data(model.velocities.w) .= ArrayType(w₀)
@@ -102,6 +120,14 @@ function run_rayleigh_benard_regression_test(arch)
     u₁, v₁, w₁ = get_output_tuple(outputwriter, spinup_steps+test_steps, :U)
     T₁, S₁ = get_output_tuple(outputwriter, spinup_steps+test_steps, :Φ)
 
+    u₁ = reverse(u₁; dims=3)
+    v₁ = reverse(v₁; dims=3)
+    T₁ = reverse(T₁; dims=3)
+    S₁ = reverse(S₁; dims=3)
+
+    w₁ = reverse(w₁[:, :, 2:Nz]; dims=3)
+    w₁ = cat(zeros(Nx, Ny), w₁; dims=3)
+
     field_names = ["u", "v", "w", "T", "S"]
     fields = [model.velocities.u, model.velocities.v, model.velocities.w, model.tracers.b, model.tracers.c]
     fields_gm = [u₁, v₁, w₁, T₁, S₁]
@@ -118,9 +144,8 @@ function run_rayleigh_benard_regression_test(arch)
     # Now test that the model state matches the regression output.
     @test all(Array(data(model.velocities.u)) .≈ u₁)
     @test all(Array(data(model.velocities.v)) .≈ v₁)
-    @test all(Array(data(model.velocities.w)) .≈ w₁)
+    @test all(Array(data(model.velocities.w))[:, :, 2:Nz] .≈ w₁[:, :, 2:Nz])
     @test all(Array(data(model.tracers.b))    .≈ T₁)
     @test all(Array(data(model.tracers.c))    .≈ S₁)
     return nothing
 end
-
diff --git a/test/regression_tests/thermal_bubble_regression_test.jl b/test/regression_tests/thermal_bubble_regression_test.jl
index 76bc728c04..51080b64d1 100644
--- a/test/regression_tests/thermal_bubble_regression_test.jl
+++ b/test/regression_tests/thermal_bubble_regression_test.jl
@@ -1,11 +1,10 @@
-
 function run_thermal_bubble_regression_test(arch)
     Nx, Ny, Nz = 16, 16, 16
     Lx, Ly, Lz = 100, 100, 100
     Δt = 6
 
-    model = BasicModel(N=(Nx, Ny, Nz), L=(Lx, Ly, Lz), architecture=arch, ν=4e-2, κ=4e-2,
-                       coriolis=FPlane(f=1e-4))
+    model = BasicModel(architecture = arch, N = (Nx, Ny, Nz), L = (Lx, Ly, Lz),
+                       ν = 4e-2, κ = 4e-2, coriolis = FPlane(f = 1e-4))
 
     model.tracers.T.data.parent .= 9.85
     model.tracers.S.data.parent .= 35.0
@@ -15,27 +14,38 @@ function run_thermal_bubble_regression_test(arch)
     i1, i2 = round(Int, Nx/4), round(Int, 3Nx/4)
     j1, j2 = round(Int, Ny/4), round(Int, 3Ny/4)
     k1, k2 = round(Int, Nz/4), round(Int, 3Nz/4)
-    model.tracers.T.data[i1:i2, j1:j2, k1:k2] .+= 0.01
-
-    outputs = Dict("v"=>model.velocities.v,
-                   "u"=>model.velocities.u,
-                   "w"=>model.velocities.w,
-                   "T"=>model.tracers.T,
-                   "S"=>model.tracers.S)
-    nc_writer = NetCDFOutputWriter(model, outputs,
-                                   filename="thermal_bubble_regression_test.nc",
-                                   frequency=10)
-    push!(model.output_writers, nc_writer)
+    model.tracers.T.data[i1:i2, j1:j2, k1+1:k2+1] .+= 0.01
+
+    # Uncomment to generate regression data.
+    # outputs = Dict("v" => model.velocities.v,
+    #                "u" => model.velocities.u,
+    #                "w" => model.velocities.w,
+    #                "T" => model.tracers.T,
+    #                "S" => model.tracers.S)
+    #
+    # nc_writer = NetCDFOutputWriter(model, outputs,
+    #                                filename="data/thermal_bubble_regression_10.nc",
+    #                                frequency=10)
+    # push!(model.output_writers, nc_writer)
 
     time_step!(model, 10, Δt)
 
-    close(nc_writer)
+    regression_data_filepath = joinpath(dirname(@__FILE__), "data", "thermal_bubble_regression_10.nc")
+    ds = Dataset(regression_data_filepath, "r")
+
+    u = ds["u"][:]
+    v = ds["v"][:]
+    w = ds["w"][:]
+    T = ds["T"][:]
+    S = ds["S"][:]
 
-    u = read_output(nc_writer, "u")
-    v = read_output(nc_writer, "v")
-    w = read_output(nc_writer, "w")
-    T = read_output(nc_writer, "T")
-    S = read_output(nc_writer, "S")
+    u = reverse(u; dims=3)
+    v = reverse(v; dims=3)
+    T = reverse(T; dims=3)
+    S = reverse(S; dims=3)
+
+    w = reverse(w[:, :, 2:Nz]; dims=3)
+    w = cat(zeros(Nx, Ny), w; dims=3)
 
     field_names = ["u", "v", "w", "T", "S"]
     fields = [model.velocities.u, model.velocities.v, model.velocities.w, model.tracers.T, model.tracers.S]
@@ -53,9 +63,7 @@ function run_thermal_bubble_regression_test(arch)
     # Now test that the model state matches the regression output.
     @test all(Array(data(model.velocities.u)) .≈ u)
     @test all(Array(data(model.velocities.v)) .≈ v)
-    @test all(Array(data(model.velocities.w)) .≈ w)
+    @test all(Array(data(model.velocities.w))[:, :, 2:Nz] .≈ w[:, :, 2:Nz])
     @test all(Array(data(model.tracers.T))    .≈ T)
     @test all(Array(data(model.tracers.S))    .≈ S)
 end
-
-
diff --git a/test/runtests.jl b/test/runtests.jl
index 49e47e64f3..1526d1c233 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -15,6 +15,7 @@ using
 using Statistics: mean
 using LinearAlgebra: norm
 using GPUifyLoops: @launch, @loop
+using NCDatasets: Dataset
 
 using Oceananigans: PoissonSolver, PPN, PNN, solve_poisson_3d!,
                     velocity_div!, compute_w_from_continuity!,
diff --git a/test/test_dynamics.jl b/test/test_dynamics.jl
index 61e8750ab0..97d4603c53 100644
--- a/test/test_dynamics.jl
+++ b/test/test_dynamics.jl
@@ -147,18 +147,20 @@ function passive_tracer_advection_test(; N=128, κ=1e-12, Nt=100)
 end
 
 """
-Pearson vortex test
-See p. 310 of "Nodal Discontinuous Galerkin Methods: Algorithms, Analysis, and Application" by Hesthaven & Warburton.
+Taylor-Green vortex test
+See: https://en.wikipedia.org/wiki/Taylor%E2%80%93Green_vortex#Taylor%E2%80%93Green_vortex_solution
+     and p. 310 of "Nodal Discontinuous Galerkin Methods: Algorithms, Analysis, and Application" by Hesthaven & Warburton.
 """
-function pearson_vortex_test(arch; FT=Float64, N=64, Nt=10)
+function taylor_green_vortex_test(arch; FT=Float64, N=64, Nt=10)
     Nx, Ny, Nz = N, N, 2
     Lx, Ly, Lz = 1, 1, 1
     ν = 1
 
-    # Choose a very small time step ~O(1/Δx²) as we are diffusion-limited in this test.
-    Δt = (Lx/Nx)^2 / (10π*ν)
+    # Choose a very small time step as we are diffusion-limited in this test: Δt ≤ Δx² / 2ν
+    Δx = Lx / Nx
+    Δt = (1/10π) * Δx^2 / ν
 
-    # Pearson vortex analytic solution.
+    # Taylor-Green vortex analytic solution.
     @inline u(x, y, z, t) = -sin(2π*y) * exp(-4π^2 * ν * t)
     @inline v(x, y, z, t) =  sin(2π*x) * exp(-4π^2 * ν * t)
 
@@ -189,7 +191,7 @@ function pearson_vortex_test(arch; FT=Float64, N=64, Nt=10)
     v_rel_err_avg = mean(v_rel_err)
     v_rel_err_max = maximum(v_rel_err)
 
-    @info "Pearson vortex test ($arch, $FT) with Nx=Ny=$N @ Nt=$Nt: " *
+    @info "Taylor-Green vortex test ($arch, $FT) with Nx=Ny=$N @ Nt=$Nt: " *
           @sprintf("Δu: (avg=%6.3g, max=%6.3g), Δv: (avg=%6.3g, max=%6.3g)\n",
                    u_rel_err_avg, u_rel_err_max, v_rel_err_avg, v_rel_err_max)
 
@@ -234,8 +236,8 @@ end
         @test internal_wave_test()
     end
 
-    @testset "Pearson vortex" begin
-        println("  Testing Pearson vortex...")
-        @test pearson_vortex_test(CPU())
+    @testset "Taylor-Green vortex" begin
+        println("  Testing Taylor-Green vortex...")
+        @test taylor_green_vortex_test(CPU())
     end
 end
diff --git a/test/test_grids.jl b/test/test_grids.jl
index 94ca4c178f..e3f135ae8b 100644
--- a/test/test_grids.jl
+++ b/test/test_grids.jl
@@ -1,31 +1,46 @@
-function correct_grid_size(ft::DataType)
-    g = RegularCartesianGrid(ft, (4, 6, 8), (2π, 4π, 9π))
+function correct_grid_size(FT)
+    grid = RegularCartesianGrid(FT, (4, 6, 8), (2π, 4π, 9π))
 
     # Checking ≈ as the grid could be storing Float32 values.
-    (g.Nx ≈ 4  && g.Ny ≈ 6  && g.Nz ≈ 8 &&
-     g.Lx ≈ 2π && g.Ly ≈ 4π && g.Lz ≈ 9π)
+    return (grid.Nx ≈ 4  && grid.Ny ≈ 6  && grid.Nz ≈ 8 &&
+            grid.Lx ≈ 2π && grid.Ly ≈ 4π && grid.Lz ≈ 9π)
 end
 
-function correct_cell_volume(ft::DataType)
+function correct_cell_volume(FT)
     Nx, Ny, Nz = 19, 13, 7
     Δx, Δy, Δz = 0.1, 0.2, 0.3
     Lx, Ly, Lz = Nx*Δx, Ny*Δy, Nz*Δz
-    g = RegularCartesianGrid(ft, (Nx, Ny, Nz), (Lx, Ly, Lz))
+    grid = RegularCartesianGrid(FT, (Nx, Ny, Nz), (Lx, Ly, Lz))
 
     # Checking ≈ as the grid could be storing Float32 values.
-    g.V ≈ Δx*Δy*Δz
+    return grid.V ≈ Δx*Δy*Δz
 end
 
-function faces_start_at_zero(ft::DataType)
-    g = RegularCartesianGrid(ft, (10, 10, 10), (2π, 2π, 2π))
-    g.xF[1] == 0 && g.yF[1] == 0 && g.zF[1] == 0
+function faces_start_at_zero(FT)
+    grid = RegularCartesianGrid(FT, (10, 10, 10), (2π, 2π, 2π))
+    return grid.xF[1] == 0 && grid.yF[1] == 0 && grid.zF[end] == 0
 end
 
-function end_faces_match_grid_length(ft::DataType)
-    g = RegularCartesianGrid(ft, (12, 13, 14), (π, π^2, π^3))
-    (g.xF[end] - g.xF[1] ≈ π   &&
-     g.yF[end] - g.yF[1] ≈ π^2 &&
-     g.zF[1] - g.zF[end] ≈ π^3)
+function end_faces_match_grid_length(FT)
+    grid = RegularCartesianGrid(FT, (12, 13, 14), (π, π^2, π^3))
+    return (grid.xF[end] - grid.xF[1] ≈ π   &&
+            grid.yF[end] - grid.yF[1] ≈ π^2 &&
+            grid.zF[end] - grid.zF[1] ≈ π^3)
+end
+
+function ranges_have_correct_length(FT)
+    Nx, Ny, Nz = 8, 9, 10
+    grid = RegularCartesianGrid(FT, (Nx, Ny, Nz), (1, 1, 1))
+    return (length(grid.xC) == Nx && length(grid.xF) == Nx+1 &&
+            length(grid.yC) == Ny && length(grid.yF) == Ny+1 &&
+            length(grid.zC) == Nz && length(grid.zF) == Nz+1)
+end
+
+# See: https://github.com/climate-machine/Oceananigans.jl/issues/480
+function no_roundoff_error_in_ranges(FT)
+    Nx, Ny, Nz = 1, 1, 64
+    grid = RegularCartesianGrid(FT, (Nx, Ny, Nz), (1, 1, π/2))
+    return length(grid.zC) == Nz
 end
 
 @testset "Grids" begin
@@ -38,6 +53,8 @@ end
             @test correct_cell_volume(FT)
             @test faces_start_at_zero(FT)
             @test end_faces_match_grid_length(FT)
+            @test ranges_have_correct_length(FT)
+            @test no_roundoff_error_in_ranges(FT)
         end
     end
 
diff --git a/test/test_time_stepping.jl b/test/test_time_stepping.jl
index 05c284a1e2..4da6789ab9 100644
--- a/test/test_time_stepping.jl
+++ b/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"
diff --git a/test/test_turbulence_closures.jl b/test/test_turbulence_closures.jl
index 2376f233eb..9aa39dee54 100644
--- a/test/test_turbulence_closures.jl
+++ b/test/test_turbulence_closures.jl
@@ -93,8 +93,7 @@ function test_anisotropic_diffusivity_fluxdiv(FT=Float64; νh=FT(0.3), κh=FT(0.
     return (   ∇_κ_∇c(2, 1, 3, grid, C.T, Val(1), closure, nothing) == 8κh + 10κv &&
             ∂ⱼ_2ν_Σ₁ⱼ(2, 1, 3, grid, closure, U..., nothing) == 2νh + 4νv &&
             ∂ⱼ_2ν_Σ₂ⱼ(2, 1, 3, grid, closure, U..., nothing) == 4νh + 6νv &&
-            ∂ⱼ_2ν_Σ₃ⱼ(2, 1, 3, grid, closure, U..., nothing) == 6νh + 8νv
-            )
+            ∂ⱼ_2ν_Σ₃ⱼ(2, 1, 3, grid, closure, U..., nothing) == 6νh + 8νv)
 end
 
 function test_function_interpolation(T=Float64)
@@ -136,9 +135,8 @@ function test_function_differentiation(T=Float64)
     ∂y_ϕ_f = ϕ²[2, 2, 2] - ϕ²[2, 1, 2]
     ∂y_ϕ_c = ϕ²[2, 3, 2] - ϕ²[2, 2, 2]
 
-    # Note reverse indexing here!
-    ∂z_ϕ_f = ϕ²[2, 2, 1] - ϕ²[2, 2, 2]
-    ∂z_ϕ_c = ϕ²[2, 2, 2] - ϕ²[2, 2, 3]
+    ∂z_ϕ_f = ϕ²[2, 2, 2] - ϕ²[2, 2, 1]
+    ∂z_ϕ_c = ϕ²[2, 2, 3] - ϕ²[2, 2, 2]
 
     f(i, j, k, grid, ϕ) = ϕ[i, j, k]^2