FINALLY reverse the k index (#462)

FINALLY reverse the k index

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

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/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

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

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,