include("dependencies_for_runtests.jl")

using TimesDates: TimeDate
using Dates: DateTime, Nanosecond, Millisecond
using TimesDates: TimeDate

using CUDA
using NCDatasets
using SeawaterPolynomials.TEOS10: TEOS10EquationOfState
using SeawaterPolynomials.SecondOrderSeawaterPolynomials: RoquetEquationOfState

using Oceananigans: Clock
using Oceananigans.Models.HydrostaticFreeSurfaceModels: VectorInvariant, ImplicitFreeSurface
using Oceananigans.OutputWriters: trilocation_dim_name, vertical_coordinate_name
using Oceananigans.Grids: ξname, ηname, rname, ξnodes, ηnodes
using Oceananigans.Fields: interpolate!

function test_datetime_netcdf_output(arch)
    grid = RectilinearGrid(arch, size=(1, 1, 1), extent=(1, 1, 1))

    clock = Clock(time=DateTime(2021, 1, 1))

    model = NonhydrostaticModel(grid; clock,
                                timestepper = :QuasiAdamsBashforth2,
                                buoyancy = SeawaterBuoyancy(),
                                tracers=(:T, :S))

    Δt = 5days + 3hours + 44.123seconds
    simulation = Simulation(model; Δt, stop_time=DateTime(2021, 2, 1))

    Arch = typeof(arch)
    filepath = "test_datetime_$Arch.nc"
    isfile(filepath) && rm(filepath)

    simulation.output_writers[:netcdf] = NetCDFWriter(model, fields(model);
                                                      filename = filepath,
                                                      schedule = IterationInterval(1),
                                                      include_grid_metrics = false)
    run!(simulation)

    ds = NCDataset(filepath)
    @test ds["time"].attrib["units"] == "seconds since 2000-01-01 00:00:00"

    Nt = length(ds["time"])
    @test Nt == 8 # There should be 8 outputs total.

    for n in 1:Nt-1
        @test ds["time"][n] == DateTime(2021, 1, 1) + (n-1) * Millisecond(1000Δt)
    end

    @test ds["time"][Nt] == DateTime(2021, 2, 1)

    close(ds)
    rm(filepath)

    return nothing
end

function test_timedate_netcdf_output(arch)
    grid = RectilinearGrid(arch, size=(1, 1, 1), extent=(1, 1, 1))

    clock = Clock(time=TimeDate(2021, 1, 1))

    model = NonhydrostaticModel(grid; clock,
                                timestepper = :QuasiAdamsBashforth2,
                                buoyancy = SeawaterBuoyancy(),
                                tracers = (:T, :S))

    Δt = 5days + 3hours + 44.123seconds
    simulation = Simulation(model, Δt=Δt, stop_time=TimeDate(2021, 2, 1))

    Arch = typeof(arch)
    filepath = "test_timedate_$Arch.nc"
    isfile(filepath) && rm(filepath)

    simulation.output_writers[:netcdf] = NetCDFWriter(model, fields(model);
                                                      filename = filepath,
                                                      schedule = IterationInterval(1),
                                                      include_grid_metrics = false)

    run!(simulation)

    ds = NCDataset(filepath)
    @test ds["time"].attrib["units"] == "seconds since 2000-01-01 00:00:00"

    Nt = length(ds["time"])
    @test Nt == 8 # There should be 8 outputs total.

    for n in 1:Nt-1
        @test ds["time"][n] == DateTime(2021, 1, 1) + (n-1) * Millisecond(1000Δt)
    end

    @test ds["time"][Nt] == DateTime(2021, 2, 1)

    close(ds)
    rm(filepath)

    return nothing
end

function test_netcdf_grid_metrics_rectilinear(arch, FT)
    Nx, Ny, Nz = 8, 6, 4
    Hx, Hy, Hz = 1, 2, 3

    grid = RectilinearGrid(arch,
                           topology = (Periodic, Bounded, Bounded),
                           size = (Nx, Ny, Nz),
                           halo = (Hx, Hy, Hz),
                           x = (0, 1), y = (0, 2), z = LinRange(0, 3, Nz + 1))

    model = NonhydrostaticModel(grid;
                                closure = ScalarDiffusivity(ν=4e-2, κ=4e-2),
                                buoyancy = SeawaterBuoyancy(),
                                tracers = (:T, :S))

    Nt = 10
    simulation = Simulation(model, Δt=0.1, stop_iteration=Nt)

    Arch = typeof(arch)
    filepath_metrics_halos = "test_grid_metrics_rectilinear_halos_$(Arch)_$FT.nc"
    isfile(filepath_metrics_halos) && rm(filepath_metrics_halos)

    simulation.output_writers[:with_metrics_and_halos] = NetCDFWriter(model, fields(model),
                                                                      filename = filepath_metrics_halos,
                                                                      schedule = IterationInterval(1),
                                                                      array_type = Array{FT},
                                                                      with_halos = true,
                                                                      include_grid_metrics = true,
                                                                      verbose = true)

    filepath_metrics_nohalos = "test_grid_metrics_rectilinear_nohalos_$(Arch)_$FT.nc"
    isfile(filepath_metrics_nohalos) && rm(filepath_metrics_nohalos)

    simulation.output_writers[:with_metrics_no_halos] = NetCDFWriter(model, fields(model),
                                                                     filename = filepath_metrics_nohalos,
                                                                     schedule = IterationInterval(1),
                                                                     array_type = Array{FT},
                                                                     with_halos = false,
                                                                     include_grid_metrics = true,
                                                                     verbose = true)

    filepath_nometrics = "test_grid_metrics_rectilinear_nometrics_$(Arch)_$FT.nc"
    isfile(filepath_nometrics) && rm(filepath_nometrics)

    simulation.output_writers[:no_metrics] = NetCDFWriter(model, fields(model),
                                                          filename = filepath_nometrics,
                                                          schedule = IterationInterval(1),
                                                          array_type = Array{FT},
                                                          with_halos = true,
                                                          include_grid_metrics = false,
                                                          verbose = true)

    i_slice = Colon()
    j_slice = 2:4
    k_slice = Nz

    nx = Nx
    ny = length(j_slice)
    nz = 1

    filepath_sliced = "test_grid_metrics_rectilinear_sliced_$(Arch)_$FT.nc"
    isfile(filepath_sliced) && rm(filepath_sliced)

    simulation.output_writers[:sliced] = NetCDFWriter(model, fields(model),
                                                      filename = filepath_sliced,
                                                      indices = (i_slice, j_slice, k_slice),
                                                      schedule = IterationInterval(1),
                                                      array_type = Array{FT},
                                                      with_halos = false,
                                                      include_grid_metrics = true,
                                                      verbose = true)

    run!(simulation)

    # Test NetCDF output with metrics and halos
    ds_mh = NCDataset(filepath_metrics_halos)

    @test haskey(ds_mh, "time")
    @test eltype(ds_mh["time"]) == Float64 # All dimensions should be Float64 by default

    dims = ("x_faa", "x_caa", "y_afa", "y_aca", "z_aaf", "z_aac")
    metrics = ("Δx_faa", "Δx_caa", "Δy_afa", "Δy_aca", "Δy_afa", "Δy_aca")
    vars = ("u", "v", "w", "T", "S")

    for var in (metrics..., vars...)
        @test haskey(ds_mh, var)
        @test haskey(ds_mh[var].attrib, "long_name")
        @test haskey(ds_mh[var].attrib, "units")
        @test eltype(ds_mh[var]) == FT
    end

    for var in dims
        @test haskey(ds_mh, var)
        @test haskey(ds_mh[var].attrib, "long_name")
        @test haskey(ds_mh[var].attrib, "units")
        @test eltype(ds_mh[var]) == Float64 # All dimensions should be Float64 by default
    end


    @test dimsize(ds_mh["time"]) == (time=Nt + 1,)

    @test dimsize(ds_mh[:x_faa]) == (x_faa=Nx + 2Hx,)
    @test dimsize(ds_mh[:x_caa]) == (x_caa=Nx + 2Hx,)
    @test dimsize(ds_mh[:y_afa]) == (y_afa=Ny + 2Hy + 1,)
    @test dimsize(ds_mh[:y_aca]) == (y_aca=Ny + 2Hy,)
    @test dimsize(ds_mh[:z_aaf]) == (z_aaf=Nz + 2Hz + 1,)
    @test dimsize(ds_mh[:z_aac]) == (z_aac=Nz + 2Hz,)

    @test dimsize(ds_mh[:Δx_faa]) == (x_faa=Nx + 2Hx,)
    @test dimsize(ds_mh[:Δx_caa]) == (x_caa=Nx + 2Hx,)
    @test dimsize(ds_mh[:Δy_afa]) == (y_afa=Ny + 2Hy + 1,)
    @test dimsize(ds_mh[:Δy_aca]) == (y_aca=Ny + 2Hy,)
    @test dimsize(ds_mh[:Δz_aaf]) == (z_aaf=Nz + 2Hz + 1,)
    @test dimsize(ds_mh[:Δz_aac]) == (z_aac=Nz + 2Hz,)

    @test dimsize(ds_mh[:u]) == (x_faa=Nx + 2Hx, y_aca=Ny + 2Hy,     z_aac=Nz + 2Hz,     time=Nt + 1)
    @test dimsize(ds_mh[:v]) == (x_caa=Nx + 2Hx, y_afa=Ny + 2Hy + 1, z_aac=Nz + 2Hz,     time=Nt + 1)
    @test dimsize(ds_mh[:w]) == (x_caa=Nx + 2Hx, y_aca=Ny + 2Hy,     z_aaf=Nz + 2Hz + 1, time=Nt + 1)
    @test dimsize(ds_mh[:T]) == (x_caa=Nx + 2Hx, y_aca=Ny + 2Hy,     z_aac=Nz + 2Hz,     time=Nt + 1)
    @test dimsize(ds_mh[:S]) == (x_caa=Nx + 2Hx, y_aca=Ny + 2Hy,     z_aac=Nz + 2Hz,     time=Nt + 1)

    close(ds_mh)
    rm(filepath_metrics_halos)

    # Test NetCDF output with metrics but not halos
    ds_m = NCDataset(filepath_metrics_nohalos)

    @test haskey(ds_m, "time")
    @test eltype(ds_m["time"]) == Float64 # All dimensions should be Float64 by default

    for var in (metrics..., vars...)
        @test haskey(ds_m, var)
        @test haskey(ds_m[var].attrib, "long_name")
        @test haskey(ds_m[var].attrib, "units")
        @test eltype(ds_m[var]) == FT
    end

    for var in dims
        @test haskey(ds_m, var)
        @test haskey(ds_m[var].attrib, "long_name")
        @test haskey(ds_m[var].attrib, "units")
        @test eltype(ds_m[var]) == Float64 # All dimensions should be Float64 by default
    end

    @test dimsize(ds_m[:x_faa]) == (x_faa=Nx,)
    @test dimsize(ds_m[:x_caa]) == (x_caa=Nx,)
    @test dimsize(ds_m[:y_afa]) == (y_afa=Ny + 1,)
    @test dimsize(ds_m[:y_aca]) == (y_aca=Ny,)
    @test dimsize(ds_m[:z_aaf]) == (z_aaf=Nz + 1,)
    @test dimsize(ds_m[:z_aac]) == (z_aac=Nz,)

    @test dimsize(ds_m[:Δx_faa]) == (x_faa=Nx,)
    @test dimsize(ds_m[:Δx_caa]) == (x_caa=Nx,)
    @test dimsize(ds_m[:Δy_afa]) == (y_afa=Ny + 1,)
    @test dimsize(ds_m[:Δy_aca]) == (y_aca=Ny,)
    @test dimsize(ds_m[:Δz_aaf]) == (z_aaf=Nz + 1,)
    @test dimsize(ds_m[:Δz_aac]) == (z_aac=Nz,)

    @test dimsize(ds_m[:u]) == (x_faa=Nx, y_aca=Ny,     z_aac=Nz,     time=Nt + 1)
    @test dimsize(ds_m[:v]) == (x_caa=Nx, y_afa=Ny + 1, z_aac=Nz,     time=Nt + 1)
    @test dimsize(ds_m[:w]) == (x_caa=Nx, y_aca=Ny,     z_aaf=Nz + 1, time=Nt + 1)
    @test dimsize(ds_m[:T]) == (x_caa=Nx, y_aca=Ny,     z_aac=Nz,     time=Nt + 1)
    @test dimsize(ds_m[:S]) == (x_caa=Nx, y_aca=Ny,     z_aac=Nz,     time=Nt + 1)

    close(ds_m)
    rm(filepath_metrics_nohalos)

    # Test NetCDF output with no metrics (with halos)
    ds_h = NCDataset(filepath_nometrics)

    @test haskey(ds_h, "time")
    @test eltype(ds_h["time"]) == Float64 # All dimensions should be Float64 by default

    for var in vars
        @test haskey(ds_h, var)
        @test haskey(ds_h[var].attrib, "long_name")
        @test haskey(ds_h[var].attrib, "units")
        @test eltype(ds_h[var]) == FT
    end

    for var in dims
        @test haskey(ds_h, var)
        @test haskey(ds_h[var].attrib, "long_name")
        @test haskey(ds_h[var].attrib, "units")
        @test eltype(ds_h[var]) == Float64 # All dimensions should be Float64 by default
    end

    for metric in metrics
        @test !haskey(ds_h, metric)
    end

    @test dimsize(ds_h["time"]) == (time=Nt + 1,)

    @test dimsize(ds_h[:x_faa]) == (x_faa=Nx + 2Hx,)
    @test dimsize(ds_h[:x_caa]) == (x_caa=Nx + 2Hx,)
    @test dimsize(ds_h[:y_afa]) == (y_afa=Ny + 2Hy + 1,)
    @test dimsize(ds_h[:y_aca]) == (y_aca=Ny + 2Hy,)
    @test dimsize(ds_h[:z_aaf]) == (z_aaf=Nz + 2Hz + 1,)
    @test dimsize(ds_h[:z_aac]) == (z_aac=Nz + 2Hz,)

    @test dimsize(ds_h[:u]) == (x_faa=Nx + 2Hx, y_aca=Ny + 2Hy,     z_aac=Nz + 2Hz,     time=Nt + 1)
    @test dimsize(ds_h[:v]) == (x_caa=Nx + 2Hx, y_afa=Ny + 2Hy + 1, z_aac=Nz + 2Hz,     time=Nt + 1)
    @test dimsize(ds_h[:w]) == (x_caa=Nx + 2Hx, y_aca=Ny + 2Hy,     z_aaf=Nz + 2Hz + 1, time=Nt + 1)
    @test dimsize(ds_h[:T]) == (x_caa=Nx + 2Hx, y_aca=Ny + 2Hy,     z_aac=Nz + 2Hz,     time=Nt + 1)
    @test dimsize(ds_h[:S]) == (x_caa=Nx + 2Hx, y_aca=Ny + 2Hy,     z_aac=Nz + 2Hz,     time=Nt + 1)

    close(ds_h)
    rm(filepath_nometrics)

    # Test NetCDF sliced output with metrics
    ds_s = NCDataset(filepath_sliced)

    @test haskey(ds_s, "time")
    @test eltype(ds_s["time"]) == Float64 # All dimensions should be Float64 by default

    for var in (metrics..., vars...)
        @test haskey(ds_s, var)
        @test haskey(ds_s[var].attrib, "long_name")
        @test haskey(ds_s[var].attrib, "units")
        @test eltype(ds_s[var]) == FT
    end

    for var in dims
        @test haskey(ds_s, var)
        @test haskey(ds_s[var].attrib, "long_name")
        @test haskey(ds_s[var].attrib, "units")
        @test eltype(ds_s[var]) == Float64 # All dimensions should be Float64 by default
    end

    @test dimsize(ds_s[:x_faa]) == (x_faa=nx,)
    @test dimsize(ds_s[:x_caa]) == (x_caa=nx,)
    @test dimsize(ds_s[:y_afa]) == (y_afa=ny,)
    @test dimsize(ds_s[:y_aca]) == (y_aca=ny,)
    @test dimsize(ds_s[:z_aaf]) == (z_aaf=nz,)
    @test dimsize(ds_s[:z_aac]) == (z_aac=nz,)

    @test dimsize(ds_s[:Δx_faa]) == (x_faa=nx,)
    @test dimsize(ds_s[:Δx_caa]) == (x_caa=nx,)
    @test dimsize(ds_s[:Δy_afa]) == (y_afa=ny,)
    @test dimsize(ds_s[:Δy_aca]) == (y_aca=ny,)
    @test dimsize(ds_s[:Δz_aaf]) == (z_aaf=nz,)
    @test dimsize(ds_s[:Δz_aac]) == (z_aac=nz,)

    @test dimsize(ds_s[:u]) == (x_faa=nx, y_aca=ny, z_aac=nz, time=Nt + 1)
    @test dimsize(ds_s[:v]) == (x_caa=nx, y_afa=ny, z_aac=nz, time=Nt + 1)
    @test dimsize(ds_s[:w]) == (x_caa=nx, y_aca=ny, z_aaf=nz, time=Nt + 1)
    @test dimsize(ds_s[:T]) == (x_caa=nx, y_aca=ny, z_aac=nz, time=Nt + 1)
    @test dimsize(ds_s[:S]) == (x_caa=nx, y_aca=ny, z_aac=nz, time=Nt + 1)

    close(ds_s)
    rm(filepath_sliced)

    return nothing
end

function test_netcdf_grid_metrics_latlon(arch, FT)
    Nλ, Nφ, Nz = 10, 9, 8
    Hλ, Hφ, Hz = 4, 3, 2

    grid = LatitudeLongitudeGrid(arch,
                                 topology = (Bounded, Bounded, Bounded),
                                 size = (Nλ, Nφ, Nz),
                                 halo = (Hλ, Hφ, Hz),
                                 longitude = (-15, 15),
                                 latitude = (-10, 10),
                                 z = LinRange(-1000, 0, Nz + 1))

    model = HydrostaticFreeSurfaceModel(grid;
                                        momentum_advection = VectorInvariant(),
                                        buoyancy = SeawaterBuoyancy(),
                                        tracers = (:T, :S))

    Nt = 5
    simulation = Simulation(model, Δt=0.1, stop_iteration=Nt)

    outputs = merge(model.velocities, model.tracers)

    Arch = typeof(arch)
    filepath_metrics_halos = "test_grid_metrics_latlon_halos_$(Arch)_$FT.nc"
    isfile(filepath_metrics_halos) && rm(filepath_metrics_halos)

    # Test with halos and metrics
    simulation.output_writers[:with_metrics_and_halos] = NetCDFWriter(model, outputs,
                                                                      filename = filepath_metrics_halos,
                                                                      schedule = IterationInterval(1),
                                                                      array_type = Array{FT},
                                                                      with_halos = true,
                                                                      include_grid_metrics = true,
                                                                      verbose = true)

    # Test with halos but no metrics
    filepath_nometrics = "test_grid_metrics_latlon_nometrics_$(Arch)_$FT.nc"
    isfile(filepath_nometrics) && rm(filepath_nometrics)

    simulation.output_writers[:no_metrics] = NetCDFWriter(model, outputs,
                                                          filename = filepath_nometrics,
                                                          schedule = IterationInterval(1),
                                                          array_type = Array{FT},
                                                          with_halos = true,
                                                          include_grid_metrics = false,
                                                          verbose = true)

    # Test without halos but with metrics
    filepath_metrics_nohalos = "test_grid_metrics_latlon_nohalos_$(Arch)_$FT.nc"
    isfile(filepath_metrics_nohalos) && rm(filepath_metrics_nohalos)

    simulation.output_writers[:with_metrics_no_halos] = NetCDFWriter(model, outputs,
                                                                     filename = filepath_metrics_nohalos,
                                                                     schedule = IterationInterval(1),
                                                                     array_type = Array{FT},
                                                                     with_halos = false,
                                                                     include_grid_metrics = true,
                                                                     verbose = true)

    # Test a slice of the domain
    i_slice = Colon()
    j_slice = 3:7
    k_slice = Nz

    nx = Nλ
    ny = length(j_slice)
    nz = 1

    filepath_sliced = "test_grid_metrics_latlon_sliced_$(Arch)_$FT.nc"
    isfile(filepath_sliced) && rm(filepath_sliced)

    simulation.output_writers[:sliced] = NetCDFWriter(model, outputs,
                                                      filename = filepath_sliced,
                                                      indices = (i_slice, j_slice, k_slice),
                                                      schedule = IterationInterval(1),
                                                      array_type = Array{FT},
                                                      with_halos = false,
                                                      include_grid_metrics = true,
                                                      verbose = true)

    run!(simulation)

    # Test NetCDF output with metrics and halos
    ds_mh = NCDataset(filepath_metrics_halos)

    @test haskey(ds_mh, "time")
    @test eltype(ds_mh["time"]) == Float64 # All dimensions should be Float64 by default

    dims = ("λ_faa", "λ_caa", "φ_afa", "φ_aca", "z_aaf", "z_aac")
    metrics = ("Δλ_faa", "Δλ_caa", "Δλ_afa", "Δλ_aca", "Δz_aaf", "Δz_aac",
               "Δx_ffa", "Δx_fca", "Δx_cfa", "Δx_cca",
               "Δy_ffa", "Δy_fca", "Δy_cfa", "Δy_cca")
    vars = ("u", "v", "w", "T", "S")

    for var in (metrics..., vars...)
        @test haskey(ds_mh, var)
        @test haskey(ds_mh[var].attrib, "long_name")
        @test haskey(ds_mh[var].attrib, "units")
        @test eltype(ds_mh[var]) == FT
    end

    for var in dims
        @test haskey(ds_mh, var)
        @test haskey(ds_mh[var].attrib, "long_name")
        @test haskey(ds_mh[var].attrib, "units")
        @test eltype(ds_mh[var]) == Float64 # All dimensions should be Float64 by default
    end

    @test dimsize(ds_mh["time"]) == (time=Nt + 1,)

    @test dimsize(ds_mh[:λ_faa]) == (λ_faa=Nλ + 2Hλ + 1,)
    @test dimsize(ds_mh[:λ_caa]) == (λ_caa=Nλ + 2Hλ,)
    @test dimsize(ds_mh[:φ_afa]) == (φ_afa=Nφ + 2Hφ + 1,)
    @test dimsize(ds_mh[:φ_aca]) == (φ_aca=Nφ + 2Hφ,)
    @test dimsize(ds_mh[:z_aaf]) == (z_aaf=Nz + 2Hz + 1,)
    @test dimsize(ds_mh[:z_aac]) == (z_aac=Nz + 2Hz,)

    @test dimsize(ds_mh[:Δλ_faa]) == (λ_faa=Nλ + 2Hλ + 1,)
    @test dimsize(ds_mh[:Δλ_caa]) == (λ_caa=Nλ + 2Hλ,)
    @test dimsize(ds_mh[:Δλ_afa]) == (φ_afa=Nφ + 2Hφ + 1,)
    @test dimsize(ds_mh[:Δλ_aca]) == (φ_aca=Nφ + 2Hφ,)
    @test dimsize(ds_mh[:Δz_aaf]) == (z_aaf=Nz + 2Hz + 1,)
    @test dimsize(ds_mh[:Δz_aac]) == (z_aac=Nz + 2Hz,)

    @test dimsize(ds_mh[:Δx_ffa]) == (λ_faa=Nλ + 2Hλ + 1, φ_afa=Nφ + 2Hφ + 1)
    @test dimsize(ds_mh[:Δx_fca]) == (λ_faa=Nλ + 2Hλ + 1, φ_aca=Nφ + 2Hφ)
    @test dimsize(ds_mh[:Δx_cfa]) == (λ_caa=Nλ + 2Hλ,     φ_afa=Nφ + 2Hφ + 1)
    @test dimsize(ds_mh[:Δx_cca]) == (λ_caa=Nλ + 2Hλ,     φ_aca=Nφ + 2Hφ)

    @test dimsize(ds_mh[:Δy_ffa]) == (λ_faa=Nλ + 2Hλ + 1, φ_afa=Nφ + 2Hφ + 1)
    @test dimsize(ds_mh[:Δy_fca]) == (λ_faa=Nλ + 2Hλ + 1, φ_aca=Nφ + 2Hφ)
    @test dimsize(ds_mh[:Δy_cfa]) == (λ_caa=Nλ + 2Hλ,     φ_afa=Nφ + 2Hφ + 1)
    @test dimsize(ds_mh[:Δy_cca]) == (λ_caa=Nλ + 2Hλ,     φ_aca=Nφ + 2Hφ)

    @test dimsize(ds_mh[:u]) == (λ_faa=Nλ + 2Hλ + 1, φ_aca=Nφ + 2Hφ,     z_aac=Nz + 2Hz,     time=Nt + 1)
    @test dimsize(ds_mh[:v]) == (λ_caa=Nλ + 2Hλ,     φ_afa=Nφ + 2Hφ + 1, z_aac=Nz + 2Hz,     time=Nt + 1)
    @test dimsize(ds_mh[:w]) == (λ_caa=Nλ + 2Hλ,     φ_aca=Nφ + 2Hφ,     z_aaf=Nz + 2Hz + 1, time=Nt + 1)
    @test dimsize(ds_mh[:T]) == (λ_caa=Nλ + 2Hλ,     φ_aca=Nφ + 2Hφ,     z_aac=Nz + 2Hz,     time=Nt + 1)
    @test dimsize(ds_mh[:S]) == (λ_caa=Nλ + 2Hλ,     φ_aca=Nφ + 2Hφ,     z_aac=Nz + 2Hz,     time=Nt + 1)

    close(ds_mh)
    rm(filepath_metrics_halos)

    # Test NetCDF output with halos but no metrics
    ds_h = NCDataset(filepath_nometrics)

    @test haskey(ds_h, "time")
    @test eltype(ds_h["time"]) == Float64 # All dimensions should be Float64 by default

    for var in vars
        @test haskey(ds_h, var)
        @test haskey(ds_h[var].attrib, "long_name")
        @test haskey(ds_h[var].attrib, "units")
        @test eltype(ds_h[var]) == FT
    end

    for var in dims
        @test haskey(ds_h, var)
        @test haskey(ds_h[var].attrib, "long_name")
        @test haskey(ds_h[var].attrib, "units")
        @test eltype(ds_h[var]) == Float64 # All dimensions should be Float64 by default
    end

    # Verify that metrics are not present
    for metric in metrics
        @test !haskey(ds_h, metric)
    end

    @test dimsize(ds_h["time"]) == (time=Nt + 1,)

    @test dimsize(ds_h[:λ_faa]) == (λ_faa=Nλ + 2Hλ + 1,)
    @test dimsize(ds_h[:λ_caa]) == (λ_caa=Nλ + 2Hλ,)
    @test dimsize(ds_h[:φ_afa]) == (φ_afa=Nφ + 2Hφ + 1,)
    @test dimsize(ds_h[:φ_aca]) == (φ_aca=Nφ + 2Hφ,)
    @test dimsize(ds_h[:z_aaf]) == (z_aaf=Nz + 2Hz + 1,)
    @test dimsize(ds_h[:z_aac]) == (z_aac=Nz + 2Hz,)

    @test dimsize(ds_h[:u]) == (λ_faa=Nλ + 2Hλ + 1, φ_aca=Nφ + 2Hφ,     z_aac=Nz + 2Hz,     time=Nt + 1)
    @test dimsize(ds_h[:v]) == (λ_caa=Nλ + 2Hλ,     φ_afa=Nφ + 2Hφ + 1, z_aac=Nz + 2Hz,     time=Nt + 1)
    @test dimsize(ds_h[:w]) == (λ_caa=Nλ + 2Hλ,     φ_aca=Nφ + 2Hφ,     z_aaf=Nz + 2Hz + 1, time=Nt + 1)
    @test dimsize(ds_h[:T]) == (λ_caa=Nλ + 2Hλ,     φ_aca=Nφ + 2Hφ,     z_aac=Nz + 2Hz,     time=Nt + 1)
    @test dimsize(ds_h[:S]) == (λ_caa=Nλ + 2Hλ,     φ_aca=Nφ + 2Hφ,     z_aac=Nz + 2Hz,     time=Nt + 1)

    close(ds_h)
    rm(filepath_nometrics)

    # Test NetCDF output with metrics but no halos
    ds_m = NCDataset(filepath_metrics_nohalos)

    for var in (metrics..., vars...)
        @test haskey(ds_m, var)
        @test haskey(ds_m[var].attrib, "long_name")
        @test haskey(ds_m[var].attrib, "units")
        @test eltype(ds_m[var]) == FT
    end

    for var in dims
        @test haskey(ds_m, var)
        @test haskey(ds_m[var].attrib, "long_name")
        @test haskey(ds_m[var].attrib, "units")
        @test eltype(ds_m[var]) == Float64 # All dimensions should be Float64 by default
    end

    @test dimsize(ds_m[:λ_faa]) == (λ_faa=Nλ + 1,)
    @test dimsize(ds_m[:λ_caa]) == (λ_caa=Nλ,)
    @test dimsize(ds_m[:φ_afa]) == (φ_afa=Nφ + 1,)
    @test dimsize(ds_m[:φ_aca]) == (φ_aca=Nφ,)
    @test dimsize(ds_m[:z_aaf]) == (z_aaf=Nz + 1,)
    @test dimsize(ds_m[:z_aac]) == (z_aac=Nz,)

    @test dimsize(ds_m[:Δλ_faa]) == (λ_faa=Nλ + 1,)
    @test dimsize(ds_m[:Δλ_caa]) == (λ_caa=Nλ,)
    @test dimsize(ds_m[:Δλ_afa]) == (φ_afa=Nφ + 1,)
    @test dimsize(ds_m[:Δλ_aca]) == (φ_aca=Nφ,)
    @test dimsize(ds_m[:Δz_aaf]) == (z_aaf=Nz + 1,)
    @test dimsize(ds_m[:Δz_aac]) == (z_aac=Nz,)

    @test dimsize(ds_m[:Δx_ffa]) == (λ_faa=Nλ + 1, φ_afa=Nφ + 1)
    @test dimsize(ds_m[:Δx_fca]) == (λ_faa=Nλ + 1, φ_aca=Nφ)
    @test dimsize(ds_m[:Δx_cfa]) == (λ_caa=Nλ,     φ_afa=Nφ + 1)
    @test dimsize(ds_m[:Δx_cca]) == (λ_caa=Nλ,     φ_aca=Nφ)

    @test dimsize(ds_m[:Δy_ffa]) == (λ_faa=Nλ + 1, φ_afa=Nφ + 1)
    @test dimsize(ds_m[:Δy_fca]) == (λ_faa=Nλ + 1, φ_aca=Nφ)
    @test dimsize(ds_m[:Δy_cfa]) == (λ_caa=Nλ,     φ_afa=Nφ + 1)
    @test dimsize(ds_m[:Δy_cca]) == (λ_caa=Nλ,     φ_aca=Nφ)

    @test dimsize(ds_m[:u]) == (λ_faa=Nλ + 1, φ_aca=Nφ,     z_aac=Nz,     time=Nt + 1)
    @test dimsize(ds_m[:v]) == (λ_caa=Nλ,     φ_afa=Nφ + 1, z_aac=Nz,     time=Nt + 1)
    @test dimsize(ds_m[:w]) == (λ_caa=Nλ,     φ_aca=Nφ,     z_aaf=Nz + 1, time=Nt + 1)
    @test dimsize(ds_m[:T]) == (λ_caa=Nλ,     φ_aca=Nφ,     z_aac=Nz,     time=Nt + 1)
    @test dimsize(ds_m[:S]) == (λ_caa=Nλ,     φ_aca=Nφ,     z_aac=Nz,     time=Nt + 1)

    close(ds_m)
    rm(filepath_metrics_nohalos)

    # Test NetCDF sliced output with metrics
    ds_s = NCDataset(filepath_sliced)

    for var in (metrics..., vars...)
        @test haskey(ds_s, var)
        @test haskey(ds_s[var].attrib, "long_name")
        @test haskey(ds_s[var].attrib, "units")
        @test eltype(ds_s[var]) == FT
    end

    for var in dims
        @test haskey(ds_s, var)
        @test haskey(ds_s[var].attrib, "long_name")
        @test haskey(ds_s[var].attrib, "units")
        @test eltype(ds_s[var]) == Float64 # All dimensions should be Float64 by default
    end

    @test dimsize(ds_s[:λ_faa]) == (λ_faa=nx + 1,)
    @test dimsize(ds_s[:λ_caa]) == (λ_caa=nx,)
    @test dimsize(ds_s[:φ_afa]) == (φ_afa=ny,)
    @test dimsize(ds_s[:φ_aca]) == (φ_aca=ny,)
    @test dimsize(ds_s[:z_aaf]) == (z_aaf=nz,)
    @test dimsize(ds_s[:z_aac]) == (z_aac=nz,)

    @test dimsize(ds_s[:Δλ_faa]) == (λ_faa=nx + 1,)
    @test dimsize(ds_s[:Δλ_caa]) == (λ_caa=nx,)
    @test dimsize(ds_s[:Δλ_afa]) == (φ_afa=ny,)
    @test dimsize(ds_s[:Δλ_aca]) == (φ_aca=ny,)
    @test dimsize(ds_s[:Δz_aaf]) == (z_aaf=nz,)
    @test dimsize(ds_s[:Δz_aac]) == (z_aac=nz,)

    @test dimsize(ds_s[:Δx_ffa]) == (λ_faa=nx + 1, φ_afa=ny)
    @test dimsize(ds_s[:Δx_fca]) == (λ_faa=nx + 1, φ_aca=ny)
    @test dimsize(ds_s[:Δx_cfa]) == (λ_caa=nx,     φ_afa=ny)
    @test dimsize(ds_s[:Δx_cca]) == (λ_caa=nx,     φ_aca=ny)

    @test dimsize(ds_s[:Δy_ffa]) == (λ_faa=nx + 1, φ_afa=ny)
    @test dimsize(ds_s[:Δy_fca]) == (λ_faa=nx + 1, φ_aca=ny)
    @test dimsize(ds_s[:Δy_cfa]) == (λ_caa=nx,     φ_afa=ny)
    @test dimsize(ds_s[:Δy_cca]) == (λ_caa=nx,     φ_aca=ny)

    @test dimsize(ds_s[:u]) == (λ_faa=nx + 1, φ_aca=ny, z_aac=nz, time=Nt + 1)
    @test dimsize(ds_s[:v]) == (λ_caa=nx,     φ_afa=ny, z_aac=nz, time=Nt + 1)
    @test dimsize(ds_s[:w]) == (λ_caa=nx,     φ_aca=ny, z_aaf=nz, time=Nt + 1)
    @test dimsize(ds_s[:T]) == (λ_caa=nx,     φ_aca=ny, z_aac=nz, time=Nt + 1)
    @test dimsize(ds_s[:S]) == (λ_caa=nx,     φ_aca=ny, z_aac=nz, time=Nt + 1)

    close(ds_s)
    rm(filepath_sliced)

    return nothing
end

function test_netcdf_rectilinear_grid_fitted_bottom(arch, bottom_boundary_type)
    Nx, Ny, Nz = 16, 16, 16
    Hx, Hy, Hz = 2, 3, 4

    Lx, Ly, H = 1, 1, 1

    underlying_grid = RectilinearGrid(arch;
                                      topology = (Bounded, Bounded, Bounded),
                                      size = (Nx, Ny, Nz),
                                      halo = (Hx, Hy, Hz),
                                      x = (-Lx, Lx),
                                      y = (-Ly, Ly),
                                      z = (-H, 0))

    height = H / 2
    width = Lx / 3
    mount(x, y) = height * exp(-x^2 / 2width^2) * exp(-y^2 / 2width^2)
    bottom(x, y) = -H + mount(x, y)

    grid = ImmersedBoundaryGrid(underlying_grid, bottom_boundary_type(bottom))

    model = NonhydrostaticModel(grid;
                                closure = ScalarDiffusivity(ν=4e-2, κ=4e-2),
                                buoyancy = SeawaterBuoyancy(),
                                tracers = (:T, :S))

    Nt = 10
    simulation = Simulation(model, Δt=0.1, stop_iteration=Nt)

    Arch = typeof(arch)
    filepath_with_halos = "test_immersed_grid_rectilinear_with_halos_$Arch.nc"
    isfile(filepath_with_halos) && rm(filepath_with_halos)

    simulation.output_writers[:with_halos] = NetCDFWriter(model, fields(model),
                                                          filename = filepath_with_halos,
                                                          schedule = IterationInterval(1),
                                                          array_type = Array{Float64},
                                                          with_halos = true,
                                                          include_grid_metrics = true,
                                                          verbose = true)

    filepath_no_halos = "test_immersed_grid_rectilinear_no_halos_$Arch.nc"
    isfile(filepath_no_halos) && rm(filepath_no_halos)

    simulation.output_writers[:no_halos] = NetCDFWriter(model, fields(model),
                                                        filename = filepath_no_halos,
                                                        schedule = IterationInterval(1),
                                                        array_type = Array{Float32},
                                                        with_halos = false,
                                                        include_grid_metrics = true,
                                                        verbose = true)

    filepath_sliced = "test_immersed_grid_rectilinear_sliced_$Arch.nc"
    isfile(filepath_sliced) && rm(filepath_sliced)

    i_slice = 4:12
    j_slice = 6:8
    k_slice = Nz

    nx = length(i_slice)
    ny = length(j_slice)
    nz = 1

    simulation.output_writers[:sliced] = NetCDFWriter(model, fields(model),
                                                      filename = filepath_sliced,
                                                      schedule = IterationInterval(1),
                                                      array_type = Array{Float32},
                                                      indices = (i_slice, j_slice, k_slice),
                                                      with_halos = false,
                                                      include_grid_metrics = true,
                                                      verbose = true)

    run!(simulation)

    # Test NetCDF output with halos
    ds_h = NCDataset(filepath_with_halos)

    @test haskey(ds_h, "bottom_height")
    @test eltype(ds_h[:bottom_height]) == Float64
    @test dimsize(ds_h[:bottom_height]) == (x_caa=Nx + 2Hx, y_aca=Ny + 2Hy)

    for loc in ("ccc", "fcc", "cfc", "ccf")
        @test haskey(ds_h, "peripheral_nodes_$loc")
        @test eltype(ds_h["peripheral_nodes_$loc"]) == Float64
        @test haskey(ds_h, "inactive_nodes_$loc")
        @test eltype(ds_h["inactive_nodes_$loc"]) == Float64
    end

    @test dimsize(ds_h[:peripheral_nodes_ccc]) == (x_caa=Nx + 2Hx,     y_aca=Ny + 2Hy,     z_aac=Nz + 2Hz)
    @test dimsize(ds_h[:peripheral_nodes_fcc]) == (x_faa=Nx + 2Hx + 1, y_aca=Ny + 2Hy,     z_aac=Nz + 2Hz)
    @test dimsize(ds_h[:peripheral_nodes_cfc]) == (x_caa=Nx + 2Hx,     y_afa=Ny + 2Hy + 1, z_aac=Nz + 2Hz)
    @test dimsize(ds_h[:peripheral_nodes_ccf]) == (x_caa=Nx + 2Hx,     y_aca=Ny + 2Hy,     z_aaf=Nz + 2Hz + 1)

    @test dimsize(ds_h[:inactive_nodes_ccc]) == (x_caa=Nx + 2Hx,     y_aca=Ny + 2Hy,     z_aac=Nz + 2Hz)
    @test dimsize(ds_h[:inactive_nodes_fcc]) == (x_faa=Nx + 2Hx + 1, y_aca=Ny + 2Hy,     z_aac=Nz + 2Hz)
    @test dimsize(ds_h[:inactive_nodes_cfc]) == (x_caa=Nx + 2Hx,     y_afa=Ny + 2Hy + 1, z_aac=Nz + 2Hz)
    @test dimsize(ds_h[:inactive_nodes_ccf]) == (x_caa=Nx + 2Hx,     y_aca=Ny + 2Hy,     z_aaf=Nz + 2Hz + 1)

    @test all(ds_h[:bottom_height][:, :] .≈ Array(parent(grid.immersed_boundary.bottom_height)))

    close(ds_h)
    rm(filepath_with_halos)

    # Test NetCDF output without halos
    ds_n = NCDataset(filepath_no_halos)

    @test haskey(ds_n, "bottom_height")
    @test eltype(ds_n[:bottom_height]) == Float32
    @test dimsize(ds_n[:bottom_height]) == (x_caa=Nx, y_aca=Ny)

    for loc in ("ccc", "fcc", "cfc", "ccf")
        @test haskey(ds_n, "peripheral_nodes_$loc")
        @test eltype(ds_n["peripheral_nodes_$loc"]) == Float32
        @test haskey(ds_n, "inactive_nodes_$loc")
        @test eltype(ds_n["inactive_nodes_$loc"]) == Float32
    end

    @test dimsize(ds_n[:peripheral_nodes_ccc]) == (x_caa=Nx,     y_aca=Ny,     z_aac=Nz)
    @test dimsize(ds_n[:peripheral_nodes_fcc]) == (x_faa=Nx + 1, y_aca=Ny,     z_aac=Nz)
    @test dimsize(ds_n[:peripheral_nodes_cfc]) == (x_caa=Nx,     y_afa=Ny + 1, z_aac=Nz)
    @test dimsize(ds_n[:peripheral_nodes_ccf]) == (x_caa=Nx,     y_aca=Ny,     z_aaf=Nz + 1)

    @test dimsize(ds_n[:inactive_nodes_ccc]) == (x_caa=Nx,     y_aca=Ny,     z_aac=Nz)
    @test dimsize(ds_n[:inactive_nodes_fcc]) == (x_faa=Nx + 1, y_aca=Ny,     z_aac=Nz)
    @test dimsize(ds_n[:inactive_nodes_cfc]) == (x_caa=Nx,     y_afa=Ny + 1, z_aac=Nz)
    @test dimsize(ds_n[:inactive_nodes_ccf]) == (x_caa=Nx,     y_aca=Ny,     z_aaf=Nz + 1)

    @test all(ds_n[:bottom_height][:, :] .≈ Array(bottom_height_interior(grid.immersed_boundary.bottom_height)))

    close(ds_n)
    rm(filepath_no_halos)

    # Test NetCDF sliced output
    ds_s = NCDataset(filepath_sliced)

    @test haskey(ds_s, "bottom_height")
    @test eltype(ds_s[:bottom_height]) == Float32
    @test dimsize(ds_s[:bottom_height]) == (x_caa=nx, y_aca=ny)

    for loc in ("ccc", "fcc", "cfc", "ccf")
        @test haskey(ds_s, "peripheral_nodes_$loc")
        @test eltype(ds_s["peripheral_nodes_$loc"]) == Float32
        @test haskey(ds_s, "inactive_nodes_$loc")
        @test eltype(ds_s["inactive_nodes_$loc"]) == Float32
    end

    @test dimsize(ds_s[:peripheral_nodes_ccc]) == (x_caa=nx, y_aca=ny, z_aac=nz)
    @test dimsize(ds_s[:peripheral_nodes_fcc]) == (x_faa=nx, y_aca=ny, z_aac=nz)
    @test dimsize(ds_s[:peripheral_nodes_cfc]) == (x_caa=nx, y_afa=ny, z_aac=nz)
    @test dimsize(ds_s[:peripheral_nodes_ccf]) == (x_caa=nx, y_aca=ny, z_aaf=nz)

    @test dimsize(ds_s[:inactive_nodes_ccc]) == (x_caa=nx, y_aca=ny, z_aac=nz)
    @test dimsize(ds_s[:inactive_nodes_fcc]) == (x_faa=nx, y_aca=ny, z_aac=nz)
    @test dimsize(ds_s[:inactive_nodes_cfc]) == (x_caa=nx, y_afa=ny, z_aac=nz)
    @test dimsize(ds_s[:inactive_nodes_ccf]) == (x_caa=nx, y_aca=ny, z_aaf=nz)

    @test all(ds_s[:bottom_height][:, :] .≈ Array(bottom_height_interior(grid.immersed_boundary.bottom_height)[i_slice, j_slice, :]))

    close(ds_s)
    rm(filepath_sliced)

    return nothing
end

function test_netcdf_latlon_grid_fitted_bottom(arch, bottom_boundary_type)
    Nλ, Nφ, Nz = 16, 16, 16
    Hλ, Hφ, Hz = 2, 3, 4
    Lλ, Lφ, H = 20, 10, 1000

    underlying_grid = LatitudeLongitudeGrid(arch;
                                            topology = (Bounded, Bounded, Bounded),
                                            size = (Nλ, Nφ, Nz),
                                            halo = (Hλ, Hφ, Hz),
                                            longitude = (-Lλ, Lλ),
                                            latitude = (-Lφ, Lφ),
                                            z = (-H, 0))

    # Create a Gaussian seamount
    height = H / 2
    λ_width = Lλ / 3
    φ_width = Lφ / 3
    seamount(λ, φ) = height * exp(-λ^2 / 2λ_width^2) * exp(-φ^2 / 2φ_width^2)
    bottom(λ, φ) = -H + seamount(λ, φ)

    grid = ImmersedBoundaryGrid(underlying_grid, bottom_boundary_type(bottom))

    model = HydrostaticFreeSurfaceModel(grid;
                                        momentum_advection = VectorInvariant(),
                                        buoyancy = SeawaterBuoyancy(),
                                        tracers = (:T, :S))

    Nt = 5
    simulation = Simulation(model, Δt=0.1, stop_iteration=Nt)

    outputs = merge(model.velocities, model.tracers)

    Arch = typeof(arch)
    filepath_with_halos = "test_immersed_grid_latlon_with_halos_$Arch.nc"
    isfile(filepath_with_halos) && rm(filepath_with_halos)

    # Test with halos with Float64 output
    simulation.output_writers[:with_halos] = NetCDFWriter(model, outputs,
                                                          filename = filepath_with_halos,
                                                          schedule = IterationInterval(1),
                                                          array_type = Array{Float64},
                                                          with_halos = true,
                                                          include_grid_metrics = true,
                                                          verbose = true)

    # Test without halos with Float32 output
    filepath_no_halos = "test_immersed_grid_latlon_no_halos_$Arch.nc"
    isfile(filepath_no_halos) && rm(filepath_no_halos)

    simulation.output_writers[:no_halos] = NetCDFWriter(model, outputs,
                                                        filename = filepath_no_halos,
                                                        schedule = IterationInterval(1),
                                                        array_type = Array{Float32},
                                                        with_halos = false,
                                                        include_grid_metrics = true,
                                                        verbose = true)

    # Test with slice
    filepath_sliced = "test_immersed_grid_latlon_sliced_$Arch.nc"
    isfile(filepath_sliced) && rm(filepath_sliced)

    i_slice = 4:12
    j_slice = 6:8
    k_slice = Nz

    nλ = length(i_slice)
    nφ = length(j_slice)
    nz = 1

    simulation.output_writers[:sliced] = NetCDFWriter(model, outputs,
                                                      filename = filepath_sliced,
                                                      schedule = IterationInterval(1),
                                                      array_type = Array{Float32},
                                                      indices = (i_slice, j_slice, k_slice),
                                                      with_halos = false,
                                                      include_grid_metrics = true,
                                                      verbose = true)

    run!(simulation)

    # Test NetCDF output with halos
    ds_h = NCDataset(filepath_with_halos)

    @test haskey(ds_h, "bottom_height")
    @test eltype(ds_h[:bottom_height]) == Float64
    @test dimsize(ds_h[:bottom_height]) == (λ_caa=Nλ + 2Hλ, φ_aca=Nφ + 2Hφ)

    for loc in ("ccc", "fcc", "cfc", "ccf")
        @test haskey(ds_h, "peripheral_nodes_$loc")
        @test eltype(ds_h["peripheral_nodes_$loc"]) == Float64
        @test haskey(ds_h, "inactive_nodes_$loc")
        @test eltype(ds_h["inactive_nodes_$loc"]) == Float64
    end

    @test dimsize(ds_h[:peripheral_nodes_ccc]) == (λ_caa=Nλ + 2Hλ,     φ_aca=Nφ + 2Hφ,     z_aac=Nz + 2Hz)
    @test dimsize(ds_h[:peripheral_nodes_fcc]) == (λ_faa=Nλ + 2Hλ + 1, φ_aca=Nφ + 2Hφ,     z_aac=Nz + 2Hz)
    @test dimsize(ds_h[:peripheral_nodes_cfc]) == (λ_caa=Nλ + 2Hλ,     φ_afa=Nφ + 2Hφ + 1, z_aac=Nz + 2Hz)
    @test dimsize(ds_h[:peripheral_nodes_ccf]) == (λ_caa=Nλ + 2Hλ,     φ_aca=Nφ + 2Hφ,     z_aaf=Nz + 2Hz + 1)

    @test dimsize(ds_h[:inactive_nodes_ccc]) == (λ_caa=Nλ + 2Hλ,     φ_aca=Nφ + 2Hφ,     z_aac=Nz + 2Hz)
    @test dimsize(ds_h[:inactive_nodes_fcc]) == (λ_faa=Nλ + 2Hλ + 1, φ_aca=Nφ + 2Hφ,     z_aac=Nz + 2Hz)
    @test dimsize(ds_h[:inactive_nodes_cfc]) == (λ_caa=Nλ + 2Hλ,     φ_afa=Nφ + 2Hφ + 1, z_aac=Nz + 2Hz)
    @test dimsize(ds_h[:inactive_nodes_ccf]) == (λ_caa=Nλ + 2Hλ,     φ_aca=Nφ + 2Hφ,     z_aaf=Nz + 2Hz + 1)

    @test all(ds_h[:bottom_height][:, :] .≈ Array(parent(grid.immersed_boundary.bottom_height)))

    close(ds_h)
    rm(filepath_with_halos)

    # Test NetCDF output without halos
    ds_n = NCDataset(filepath_no_halos)

    @test haskey(ds_n, "bottom_height")
    @test eltype(ds_n[:bottom_height]) == Float32
    @test dimsize(ds_n[:bottom_height]) == (λ_caa=Nλ, φ_aca=Nφ)

    for loc in ("ccc", "fcc", "cfc", "ccf")
        @test haskey(ds_n, "peripheral_nodes_$loc")
        @test eltype(ds_n["peripheral_nodes_$loc"]) == Float32
        @test haskey(ds_n, "inactive_nodes_$loc")
        @test eltype(ds_n["inactive_nodes_$loc"]) == Float32
    end

    @test dimsize(ds_n[:peripheral_nodes_ccc]) == (λ_caa=Nλ,     φ_aca=Nφ,     z_aac=Nz)
    @test dimsize(ds_n[:peripheral_nodes_fcc]) == (λ_faa=Nλ + 1, φ_aca=Nφ,     z_aac=Nz)
    @test dimsize(ds_n[:peripheral_nodes_cfc]) == (λ_caa=Nλ,     φ_afa=Nφ + 1, z_aac=Nz)
    @test dimsize(ds_n[:peripheral_nodes_ccf]) == (λ_caa=Nλ,     φ_aca=Nφ,     z_aaf=Nz + 1)

    @test dimsize(ds_n[:inactive_nodes_ccc]) == (λ_caa=Nλ,     φ_aca=Nφ,     z_aac=Nz)
    @test dimsize(ds_n[:inactive_nodes_fcc]) == (λ_faa=Nλ + 1, φ_aca=Nφ,     z_aac=Nz)
    @test dimsize(ds_n[:inactive_nodes_cfc]) == (λ_caa=Nλ,     φ_afa=Nφ + 1, z_aac=Nz)
    @test dimsize(ds_n[:inactive_nodes_ccf]) == (λ_caa=Nλ,     φ_aca=Nφ,     z_aaf=Nz + 1)

    @test all(ds_n[:bottom_height][:, :] .≈ Array(bottom_height_interior(grid.immersed_boundary.bottom_height)))

    close(ds_n)
    rm(filepath_no_halos)

    # Test NetCDF sliced output
    ds_s = NCDataset(filepath_sliced)

    @test haskey(ds_s, "bottom_height")
    @test eltype(ds_s[:bottom_height]) == Float32
    @test dimsize(ds_s[:bottom_height]) == (λ_caa=nλ, φ_aca=nφ)

    for loc in ("ccc", "fcc", "cfc", "ccf")
        @test haskey(ds_s, "peripheral_nodes_$loc")
        @test eltype(ds_s["peripheral_nodes_$loc"]) == Float32
        @test haskey(ds_s, "inactive_nodes_$loc")
        @test eltype(ds_s["inactive_nodes_$loc"]) == Float32
    end

    @test dimsize(ds_s[:peripheral_nodes_ccc]) == (λ_caa=nλ, φ_aca=nφ, z_aac=nz)
    @test dimsize(ds_s[:peripheral_nodes_fcc]) == (λ_faa=nλ, φ_aca=nφ, z_aac=nz)
    @test dimsize(ds_s[:peripheral_nodes_cfc]) == (λ_caa=nλ, φ_afa=nφ, z_aac=nz)
    @test dimsize(ds_s[:peripheral_nodes_ccf]) == (λ_caa=nλ, φ_aca=nφ, z_aaf=nz)

    @test dimsize(ds_s[:inactive_nodes_ccc]) == (λ_caa=nλ, φ_aca=nφ, z_aac=nz)
    @test dimsize(ds_s[:inactive_nodes_fcc]) == (λ_faa=nλ, φ_aca=nφ, z_aac=nz)
    @test dimsize(ds_s[:inactive_nodes_cfc]) == (λ_caa=nλ, φ_afa=nφ, z_aac=nz)
    @test dimsize(ds_s[:inactive_nodes_ccf]) == (λ_caa=nλ, φ_aca=nφ, z_aaf=nz)

    @test all(ds_s[:bottom_height][:, :] .≈ Array(bottom_height_interior(grid.immersed_boundary.bottom_height)[i_slice, j_slice, :]))

    close(ds_s)
    rm(filepath_sliced)

    return nothing
end

function test_netcdf_rectilinear_flat_xy(arch)
    Nx, Ny = 8, 8
    Hx, Hy = 2, 3

    grid = RectilinearGrid(arch,
                           topology = (Periodic, Bounded, Flat),
                           size = (Nx, Ny),
                           halo = (Hx, Hy),
                           extent = (π, 7))

    model = NonhydrostaticModel(grid;
                                closure = ScalarDiffusivity(ν=4e-2, κ=4e-2),
                                buoyancy = SeawaterBuoyancy(),
                                tracers = (:T, :S))

    Nt = 7
    simulation = Simulation(model, Δt=0.1, stop_iteration=Nt)

    Arch = typeof(arch)
    filepath_with_halos = "test_netcdf_rectilinear_flat_xy_$Arch.nc"
    isfile(filepath_with_halos) && rm(filepath_with_halos)

    simulation.output_writers[:with_halos] = NetCDFWriter(model, fields(model),
                                                          filename = filepath_with_halos,
                                                          schedule = IterationInterval(1),
                                                          array_type = Array{Float64},
                                                          with_halos = true,
                                                          include_grid_metrics = true,
                                                          verbose = true)

    i_slice = 3:6
    j_slice = Ny

    nx = length(i_slice)
    ny = 1

    filepath_sliced = "test_netcdf_rectilinear_flat_xy_sliced_$(Arch).nc"
    isfile(filepath_sliced) && rm(filepath_sliced)

    simulation.output_writers[:sliced] = NetCDFWriter(model, fields(model),
                                                      filename = filepath_sliced,
                                                      indices = (i_slice, j_slice, :),
                                                      schedule = IterationInterval(1),
                                                      array_type = Array{Float64},
                                                      with_halos = false,
                                                      include_grid_metrics = true,
                                                      verbose = true)

    run!(simulation)

    # Test NetCDF output with halos
    ds_h = NCDataset(filepath_with_halos)

    dims = ("x_faa", "x_caa", "y_afa", "y_aca")
    not_dims = ("z_aaf", "z_aac")

    metrics = ("Δx_faa", "Δx_caa", "Δy_afa", "Δy_aca")
    not_metrics = ("Δz_aaf", "Δz_aac")

    vars = ("u", "v", "w", "T", "S")

    for var in (dims..., metrics..., vars...)
        @test haskey(ds_h, var)
        @test haskey(ds_h[var].attrib, "long_name")
        @test haskey(ds_h[var].attrib, "units")
        @test eltype(ds_h[var]) == Float64
    end

    for var in (not_dims..., not_metrics...)
        @test !haskey(ds_h, var)
    end

    @test dimsize(ds_h[:x_faa]) == (x_faa=Nx + 2Hx,)
    @test dimsize(ds_h[:x_caa]) == (x_caa=Nx + 2Hx,)
    @test dimsize(ds_h[:y_afa]) == (y_afa=Ny + 2Hy + 1,)
    @test dimsize(ds_h[:y_aca]) == (y_aca=Ny + 2Hy,)

    @test dimsize(ds_h[:Δx_faa]) == (x_faa=Nx + 2Hx,)
    @test dimsize(ds_h[:Δx_caa]) == (x_caa=Nx + 2Hx,)
    @test dimsize(ds_h[:Δy_afa]) == (y_afa=Ny + 2Hy + 1,)
    @test dimsize(ds_h[:Δy_aca]) == (y_aca=Ny + 2Hy,)

    @test dimsize(ds_h[:u]) == (x_faa=Nx + 2Hx, y_aca=Ny + 2Hy,     time=Nt + 1)
    @test dimsize(ds_h[:v]) == (x_caa=Nx + 2Hx, y_afa=Ny + 2Hy + 1, time=Nt + 1)
    @test dimsize(ds_h[:w]) == (x_caa=Nx + 2Hx, y_aca=Ny + 2Hy,     time=Nt + 1)
    @test dimsize(ds_h[:T]) == (x_caa=Nx + 2Hx, y_aca=Ny + 2Hy,     time=Nt + 1)
    @test dimsize(ds_h[:S]) == (x_caa=Nx + 2Hx, y_aca=Ny + 2Hy,     time=Nt + 1)

    close(ds_h)
    rm(filepath_with_halos)

    # Test NetCDF sliced output
    ds_s = NCDataset(filepath_sliced)

    for var in (dims..., metrics..., vars...)
        @test haskey(ds_s, var)
        @test haskey(ds_s[var].attrib, "long_name")
        @test haskey(ds_s[var].attrib, "units")
        @test eltype(ds_s[var]) == Float64
    end

    for var in (not_dims..., not_metrics...)
        @test !haskey(ds_s, var)
    end

    @test dimsize(ds_s[:x_faa]) == (x_faa=nx,)
    @test dimsize(ds_s[:x_caa]) == (x_caa=nx,)
    @test dimsize(ds_s[:y_afa]) == (y_afa=ny,)
    @test dimsize(ds_s[:y_aca]) == (y_aca=ny,)

    @test dimsize(ds_s[:Δx_faa]) == (x_faa=nx,)
    @test dimsize(ds_s[:Δx_caa]) == (x_caa=nx,)
    @test dimsize(ds_s[:Δy_afa]) == (y_afa=ny,)
    @test dimsize(ds_s[:Δy_aca]) == (y_aca=ny,)

    @test dimsize(ds_s[:u]) == (x_faa=nx, y_aca=ny, time=Nt + 1)
    @test dimsize(ds_s[:v]) == (x_caa=nx, y_afa=ny, time=Nt + 1)
    @test dimsize(ds_s[:w]) == (x_caa=nx, y_aca=ny, time=Nt + 1)
    @test dimsize(ds_s[:T]) == (x_caa=nx, y_aca=ny, time=Nt + 1)
    @test dimsize(ds_s[:S]) == (x_caa=nx, y_aca=ny, time=Nt + 1)

    close(ds_s)
    rm(filepath_sliced)

    return nothing
end

function test_netcdf_rectilinear_flat_xz(arch; immersed)
    Nx, Nz = 8, 8
    Hx, Hz = 2, 3
    Lx, H  = 2, 1

    grid = RectilinearGrid(arch,
                           topology = (Periodic, Flat, Bounded),
                           size = (Nx, Nz),
                           halo = (Hx, Hz),
                           x = (-Lx, Lx),
                           z = (-H, 0))

    if immersed
        height = H / 2
        width = Lx / 3
        mount(x) = height * exp(-x^2 / 2width^2)
        bottom(x) = -H + mount(x)

        grid = ImmersedBoundaryGrid(grid, GridFittedBottom(bottom))
    end

    model = NonhydrostaticModel(grid;
                                closure = ScalarDiffusivity(ν=4e-2, κ=4e-2),
                                buoyancy = SeawaterBuoyancy(),
                                tracers = (:T, :S))

    Nt = 7
    simulation = Simulation(model, Δt=0.1, stop_iteration=Nt)

    Arch = typeof(arch)
    filepath_with_halos = "test_netcdf_rectilinear_flat_xz_$Arch.nc"
    isfile(filepath_with_halos) && rm(filepath_with_halos)

    simulation.output_writers[:with_halos] = NetCDFWriter(model, fields(model),
                                                          filename = filepath_with_halos,
                                                          schedule = IterationInterval(1),
                                                          array_type = Array{Float64},
                                                          with_halos = true,
                                                          include_grid_metrics = true,
                                                          verbose = true)

    i_slice = 3:6
    k_slice = Nz

    nx = length(i_slice)
    nz = 1

    filepath_sliced = "test_netcdf_rectilinear_flat_xz_sliced_$(Arch).nc"
    isfile(filepath_sliced) && rm(filepath_sliced)

    simulation.output_writers[:sliced] = NetCDFWriter(model, fields(model),
                                                      filename = filepath_sliced,
                                                      indices = (i_slice, :, k_slice),
                                                      schedule = IterationInterval(1),
                                                      array_type = Array{Float64},
                                                      with_halos = false,
                                                      include_grid_metrics = true,
                                                      verbose = true)

    run!(simulation)

    # Test NetCDF output with halos
    ds_h = NCDataset(filepath_with_halos)

    dims = ("x_faa", "x_caa", "z_aaf", "z_aac")
    not_dims = ("y_afa", "y_aca")

    metrics = ("Δx_faa", "Δx_caa", "Δz_aaf", "Δz_aac")
    not_metrics = ("Δy_afa", "Δy_aca")

    vars = ("u", "v", "w", "T", "S")

    for var in (dims..., metrics..., vars...)
        @test haskey(ds_h, var)
        @test haskey(ds_h[var].attrib, "long_name")
        @test haskey(ds_h[var].attrib, "units")
        @test eltype(ds_h[var]) == Float64
    end

    for var in (not_dims..., not_metrics...)
        @test !haskey(ds_h, var)
    end

    @test dimsize(ds_h[:x_faa]) == (x_faa=Nx + 2Hx,)
    @test dimsize(ds_h[:x_caa]) == (x_caa=Nx + 2Hx,)
    @test dimsize(ds_h[:z_aaf]) == (z_aaf=Nz + 2Hz + 1,)
    @test dimsize(ds_h[:z_aac]) == (z_aac=Nz + 2Hz,)

    @test dimsize(ds_h[:Δx_faa]) == (x_faa=Nx + 2Hx,)
    @test dimsize(ds_h[:Δx_caa]) == (x_caa=Nx + 2Hx,)
    @test dimsize(ds_h[:Δz_aaf]) == (z_aaf=Nz + 2Hz + 1,)
    @test dimsize(ds_h[:Δz_aac]) == (z_aac=Nz + 2Hz,)

    @test dimsize(ds_h[:u]) == (x_faa=Nx + 2Hx, z_aac=Nz + 2Hz,     time=Nt + 1)
    @test dimsize(ds_h[:v]) == (x_caa=Nx + 2Hx, z_aac=Nz + 2Hz,     time=Nt + 1)
    @test dimsize(ds_h[:w]) == (x_caa=Nx + 2Hx, z_aaf=Nz + 2Hz + 1, time=Nt + 1)
    @test dimsize(ds_h[:T]) == (x_caa=Nx + 2Hx, z_aac=Nz + 2Hz,     time=Nt + 1)
    @test dimsize(ds_h[:S]) == (x_caa=Nx + 2Hx, z_aac=Nz + 2Hz,     time=Nt + 1)

    close(ds_h)
    rm(filepath_with_halos)

    # Test NetCDF sliced output
    ds_s = NCDataset(filepath_sliced)

    for var in (dims..., metrics..., vars...)
        @test haskey(ds_s, var)
        @test haskey(ds_s[var].attrib, "long_name")
        @test haskey(ds_s[var].attrib, "units")
        @test eltype(ds_s[var]) == Float64
    end

    for var in (not_dims..., not_metrics...)
        @test !haskey(ds_s, var)
    end

    @test dimsize(ds_s[:x_faa]) == (x_faa=nx,)
    @test dimsize(ds_s[:x_caa]) == (x_caa=nx,)
    @test dimsize(ds_s[:z_aaf]) == (z_aaf=nz,)
    @test dimsize(ds_s[:z_aac]) == (z_aac=nz,)

    @test dimsize(ds_s[:Δx_faa]) == (x_faa=nx,)
    @test dimsize(ds_s[:Δx_caa]) == (x_caa=nx,)
    @test dimsize(ds_s[:Δz_aaf]) == (z_aaf=nz,)
    @test dimsize(ds_s[:Δz_aac]) == (z_aac=nz,)

    @test dimsize(ds_s[:u]) == (x_faa=nx, z_aac=nz, time=Nt + 1)
    @test dimsize(ds_s[:v]) == (x_caa=nx, z_aac=nz, time=Nt + 1)
    @test dimsize(ds_s[:w]) == (x_caa=nx, z_aaf=nz, time=Nt + 1)
    @test dimsize(ds_s[:T]) == (x_caa=nx, z_aac=nz, time=Nt + 1)
    @test dimsize(ds_s[:S]) == (x_caa=nx, z_aac=nz, time=Nt + 1)

    close(ds_s)
    rm(filepath_sliced)

    return nothing
end

function test_netcdf_rectilinear_flat_yz(arch; immersed)
    Ny, Nz = 8, 8
    Hy, Hz = 2, 3
    Ly, H  = 2, 1

    grid = RectilinearGrid(arch,
                           topology = (Flat, Periodic, Bounded),
                           size = (Ny, Nz),
                           halo = (Hy, Hz),
                           y = (-Ly, Ly),
                           z = (-H, 0))

    if immersed
        height = H / 2
        width = Ly / 3
        mount(x) = height * exp(-x^2 / 2width^2)
        bottom(x) = -H + mount(x)

        grid = ImmersedBoundaryGrid(grid, GridFittedBottom(bottom))
    end

    model = NonhydrostaticModel(grid;
                                closure = ScalarDiffusivity(ν=4e-2, κ=4e-2),
                                buoyancy = SeawaterBuoyancy(),
                                tracers = (:T, :S))

    Nt = 7
    simulation = Simulation(model, Δt=0.1, stop_iteration=Nt)

    Arch = typeof(arch)
    filepath_with_halos = "test_netcdf_rectilinear_flat_yz_$Arch.nc"
    isfile(filepath_with_halos) && rm(filepath_with_halos)

    simulation.output_writers[:with_halos] =
        NetCDFWriter(model, fields(model),
            filename = filepath_with_halos,
            schedule = IterationInterval(1),
            array_type = Array{Float64},
            with_halos = true,
            include_grid_metrics = true,
            verbose = true)

    j_slice = 3:6
    k_slice = Nz

    ny = length(j_slice)
    nz = 1

    filepath_sliced = "test_netcdf_rectilinear_flat_yz_sliced_$(Arch).nc"
    isfile(filepath_sliced) && rm(filepath_sliced)

    simulation.output_writers[:sliced] =
        NetCDFWriter(model, fields(model),
            filename = filepath_sliced,
            indices = (:, j_slice, k_slice),
            schedule = IterationInterval(1),
            array_type = Array{Float64},
            with_halos = false,
            include_grid_metrics = true,
            verbose = true)

    run!(simulation)

    # Test NetCDF output with halos
    ds_h = NCDataset(filepath_with_halos)

    dims = ("y_afa", "y_aca", "z_aaf", "z_aac")
    not_dims = ("x_faa", "x_caa")

    metrics = ("Δy_afa", "Δy_aca", "Δz_aaf", "Δz_aac")
    not_metrics = ("Δx_faa", "Δx_caa")

    vars = ("u", "v", "w", "T", "S")

    for var in (dims..., metrics..., vars...)
        @test haskey(ds_h, var)
        @test haskey(ds_h[var].attrib, "long_name")
        @test haskey(ds_h[var].attrib, "units")
        @test eltype(ds_h[var]) == Float64
    end

    for var in (not_dims..., not_metrics...)
        @test !haskey(ds_h, var)
    end

    @test dimsize(ds_h[:y_afa]) == (y_afa=Ny + 2Hy,)
    @test dimsize(ds_h[:y_aca]) == (y_aca=Ny + 2Hy,)
    @test dimsize(ds_h[:z_aaf]) == (z_aaf=Nz + 2Hz + 1,)
    @test dimsize(ds_h[:z_aac]) == (z_aac=Nz + 2Hz,)

    @test dimsize(ds_h[:Δy_afa]) == (y_afa=Ny + 2Hy,)
    @test dimsize(ds_h[:Δy_aca]) == (y_aca=Ny + 2Hy,)
    @test dimsize(ds_h[:Δz_aaf]) == (z_aaf=Nz + 2Hz + 1,)
    @test dimsize(ds_h[:Δz_aac]) == (z_aac=Nz + 2Hz,)

    @test dimsize(ds_h[:u]) == (y_aca=Ny + 2Hy, z_aac=Nz + 2Hz,     time=Nt + 1)
    @test dimsize(ds_h[:v]) == (y_afa=Ny + 2Hy, z_aac=Nz + 2Hz,     time=Nt + 1)
    @test dimsize(ds_h[:w]) == (y_aca=Ny + 2Hy, z_aaf=Nz + 2Hz + 1, time=Nt + 1)
    @test dimsize(ds_h[:T]) == (y_aca=Ny + 2Hy, z_aac=Nz + 2Hz,     time=Nt + 1)
    @test dimsize(ds_h[:S]) == (y_aca=Ny + 2Hy, z_aac=Nz + 2Hz,     time=Nt + 1)

    close(ds_h)
    rm(filepath_with_halos)

    # Test NetCDF sliced output
    ds_s = NCDataset(filepath_sliced)

    for var in (dims..., metrics..., vars...)
        @test haskey(ds_s, var)
        @test haskey(ds_s[var].attrib, "long_name")
        @test haskey(ds_s[var].attrib, "units")
        @test eltype(ds_s[var]) == Float64
    end

    for var in (not_dims..., not_metrics...)
        @test !haskey(ds_s, var)
    end

    @test dimsize(ds_s[:y_afa]) == (y_afa=ny,)
    @test dimsize(ds_s[:y_aca]) == (y_aca=ny,)
    @test dimsize(ds_s[:z_aaf]) == (z_aaf=nz,)
    @test dimsize(ds_s[:z_aac]) == (z_aac=nz,)

    @test dimsize(ds_s[:Δy_afa]) == (y_afa=ny,)
    @test dimsize(ds_s[:Δy_aca]) == (y_aca=ny,)
    @test dimsize(ds_s[:Δz_aaf]) == (z_aaf=nz,)
    @test dimsize(ds_s[:Δz_aac]) == (z_aac=nz,)

    @test dimsize(ds_s[:u]) == (y_aca=ny, z_aac=nz, time=Nt + 1)
    @test dimsize(ds_s[:v]) == (y_afa=ny, z_aac=nz, time=Nt + 1)
    @test dimsize(ds_s[:w]) == (y_aca=ny, z_aaf=nz, time=Nt + 1)
    @test dimsize(ds_s[:T]) == (y_aca=ny, z_aac=nz, time=Nt + 1)
    @test dimsize(ds_s[:S]) == (y_aca=ny, z_aac=nz, time=Nt + 1)

    close(ds_s)
    rm(filepath_sliced)

    return nothing
end

function test_netcdf_rectilinear_column(arch)
    N = 17
    H = 2

    grid = RectilinearGrid(arch,
                           topology = (Flat, Flat, Bounded),
                           size = N,
                           halo = H,
                           z = (-55, 0))

    model = NonhydrostaticModel(grid;
                                closure = ScalarDiffusivity(ν=4e-2, κ=4e-2),
                                buoyancy = SeawaterBuoyancy(),
                                tracers = (:T, :S))

    Nt = 5
    simulation = Simulation(model, Δt=0.1, stop_iteration=Nt)

    Arch = typeof(arch)
    filepath_with_halos = "test_netcdf_rectilinear_column_$Arch.nc"
    isfile(filepath_with_halos) && rm(filepath_with_halos)

    simulation.output_writers[:with_halos] =
        NetCDFWriter(model, fields(model),
            filename = filepath_with_halos,
            schedule = IterationInterval(1),
            array_type = Array{Float64},
            with_halos = true,
            include_grid_metrics = true,
            verbose = true)

    k_slice = 5:10
    n = length(k_slice)

    filepath_sliced = "test_netcdf_rectilinear_column_sliced_$(Arch).nc"
    isfile(filepath_sliced) && rm(filepath_sliced)

    simulation.output_writers[:sliced] =
        NetCDFWriter(model, fields(model),
            filename = filepath_sliced,
            indices = (:, :, k_slice),
            schedule = IterationInterval(1),
            array_type = Array{Float64},
            with_halos = false,
            include_grid_metrics = true,
            verbose = true)

    run!(simulation)

    # Test NetCDF output with halos
    ds_h = NCDataset(filepath_with_halos)

    dims = ("z_aaf", "z_aac")
    not_dims = ("x_faa", "x_caa", "y_afa", "y_aca")

    metrics = ("Δz_aaf", "Δz_aac")
    not_metrics = ("Δx_faa", "Δx_caa", "Δy_afa", "Δy_aca")

    vars = ("u", "v", "w", "T", "S")

    for var in (dims..., metrics..., vars...)
        @test haskey(ds_h, var)
        @test haskey(ds_h[var].attrib, "long_name")
        @test haskey(ds_h[var].attrib, "units")
        @test eltype(ds_h[var]) == Float64
    end

    for var in (not_dims..., not_metrics...)
        @test !haskey(ds_h, var)
    end

    @test dimsize(ds_h[:z_aaf]) == (z_aaf=N + 2H + 1,)
    @test dimsize(ds_h[:z_aac]) == (z_aac=N + 2H,)

    @test dimsize(ds_h[:Δz_aaf]) == (z_aaf=N + 2H + 1,)
    @test dimsize(ds_h[:Δz_aac]) == (z_aac=N + 2H,)

    @test dimsize(ds_h[:u]) == (z_aac=N + 2H,     time=Nt + 1)
    @test dimsize(ds_h[:v]) == (z_aac=N + 2H,     time=Nt + 1)
    @test dimsize(ds_h[:w]) == (z_aaf=N + 2H + 1, time=Nt + 1)
    @test dimsize(ds_h[:T]) == (z_aac=N + 2H,     time=Nt + 1)
    @test dimsize(ds_h[:S]) == (z_aac=N + 2H,     time=Nt + 1)

    close(ds_h)
    rm(filepath_with_halos)

    # Test NetCDF sliced output
    ds_s = NCDataset(filepath_sliced)

    for var in (dims..., metrics..., vars...)
        @test haskey(ds_s, var)
        @test haskey(ds_s[var].attrib, "long_name")
        @test haskey(ds_s[var].attrib, "units")
        @test eltype(ds_s[var]) == Float64
    end

    for var in (not_dims..., not_metrics...)
        @test !haskey(ds_s, var)
    end

    @test dimsize(ds_s[:z_aaf]) == (z_aaf=n,)
    @test dimsize(ds_s[:z_aac]) == (z_aac=n,)

    @test dimsize(ds_s[:Δz_aaf]) == (z_aaf=n,)
    @test dimsize(ds_s[:Δz_aac]) == (z_aac=n,)

    @test dimsize(ds_s[:u]) == (z_aac=n, time=Nt + 1)
    @test dimsize(ds_s[:v]) == (z_aac=n, time=Nt + 1)
    @test dimsize(ds_s[:w]) == (z_aaf=n, time=Nt + 1)
    @test dimsize(ds_s[:T]) == (z_aac=n, time=Nt + 1)
    @test dimsize(ds_s[:S]) == (z_aac=n, time=Nt + 1)

    close(ds_s)
    rm(filepath_sliced)

    return nothing
end

function test_thermal_bubble_netcdf_output(arch, FT; with_halos=false)
    Nx, Ny, Nz = 16, 16, 16
    Lx, Ly, Lz = 100, 100, 100
    Hx, Hy, Hz = 4, 3, 2

    grid = RectilinearGrid(arch,
                           topology = (Periodic, Periodic, Bounded),
                           size = (Nx, Ny, Nz),
                           halo = (Hx, Hy, Hz),
                           extent = (Lx, Ly, Lz))

    model = NonhydrostaticModel(grid;
                                closure = ScalarDiffusivity(ν=4e-2, κ=4e-2),
                                buoyancy = SeawaterBuoyancy(),
                                tracers = (:T, :S))

    simulation = Simulation(model, Δt=6, stop_iteration=10)

    # Add a cube-shaped warm temperature anomaly that takes up the middle 50%
    # of the domain volume.
    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)
    view(model.tracers.T, 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
    )

    Arch = typeof(arch)
    halo_suffix = with_halos ? "_with_halos" : ""
    nc_filepath = "test_thermal_bubble$(halo_suffix)_$(Arch)_$FT.nc"
    isfile(nc_filepath) && rm(nc_filepath)

    simulation.output_writers[:nc_writer] = NetCDFWriter(model, outputs,
        filename = nc_filepath,
        schedule = IterationInterval(10),
        array_type = Array{FT},
        with_halos = with_halos,
        include_grid_metrics = false,
        verbose = true)

    i_slice = 1:10
    j_slice = 13
    k_slice = 9:11
    indices = (i_slice, j_slice, k_slice)
    j_slice = j_slice:j_slice  # So we can correctly index with it for later tests.

    nc_sliced_filepath = "test_thermal_bubble_sliced_$(Arch)_$FT.nc"
    isfile(nc_sliced_filepath) && rm(nc_sliced_filepath)

    simulation.output_writers[:nc_sliced_writer] = NetCDFWriter(model, outputs,
        filename = nc_sliced_filepath,
        schedule = IterationInterval(10),
        array_type = Array{FT},
        indices = indices,
        include_grid_metrics = false,
        verbose = true)

    run!(simulation)

    ds = NCDataset(nc_filepath)

    @test haskey(ds.attrib, "date")
    @test haskey(ds.attrib, "Julia")
    @test haskey(ds.attrib, "Oceananigans")
    @test haskey(ds.attrib, "schedule")
    @test haskey(ds.attrib, "interval")
    @test haskey(ds.attrib, "output iteration interval")

    @test !isnothing(ds.attrib["date"])
    @test !isnothing(ds.attrib["Julia"])
    @test !isnothing(ds.attrib["Oceananigans"])
    @test ds.attrib["schedule"] == "IterationInterval"
    @test ds.attrib["interval"] == 10
    @test !isnothing(ds.attrib["output iteration interval"])

    # Dimensions should be Float64 by default
    @test eltype(ds["time"])  == Float64
    @test eltype(ds["x_caa"]) == Float64
    @test eltype(ds["x_faa"]) == Float64
    @test eltype(ds["y_aca"]) == Float64
    @test eltype(ds["y_afa"]) == Float64
    @test eltype(ds["z_aac"]) == Float64
    @test eltype(ds["z_aaf"]) == Float64

    if with_halos
        @test length(ds["x_caa"]) == Nx+2Hx
        @test length(ds["y_aca"]) == Ny+2Hy
        @test length(ds["z_aac"]) == Nz+2Hz
        @test length(ds["x_faa"]) == Nx+2Hx
        @test length(ds["y_afa"]) == Ny+2Hy
        @test length(ds["z_aaf"]) == Nz+2Hz+1  # z is Bounded

        @test ds["x_caa"][1] == grid.xᶜᵃᵃ[1-Hx]
        @test ds["x_faa"][1] == grid.xᶠᵃᵃ[1-Hx]
        @test ds["y_aca"][1] == grid.yᵃᶜᵃ[1-Hy]
        @test ds["y_afa"][1] == grid.yᵃᶠᵃ[1-Hy]
        @test ds["z_aac"][1] == grid.z.cᵃᵃᶜ[1-Hz]
        @test ds["z_aaf"][1] == grid.z.cᵃᵃᶠ[1-Hz]

        @test ds["x_caa"][end] == grid.xᶜᵃᵃ[Nx+Hx]
        @test ds["x_faa"][end] == grid.xᶠᵃᵃ[Nx+Hx]
        @test ds["y_aca"][end] == grid.yᵃᶜᵃ[Ny+Hy]
        @test ds["y_afa"][end] == grid.yᵃᶠᵃ[Ny+Hy]
        @test ds["z_aac"][end] == grid.z.cᵃᵃᶜ[Nz+Hz]
        @test ds["z_aaf"][end] == grid.z.cᵃᵃᶠ[Nz+Hz+1]  # z is Bounded
    else
        @test length(ds["x_caa"]) == Nx
        @test length(ds["y_aca"]) == Ny
        @test length(ds["z_aac"]) == Nz
        @test length(ds["x_faa"]) == Nx
        @test length(ds["y_afa"]) == Ny
        @test length(ds["z_aaf"]) == Nz+1  # z is Bounded

        @test ds["x_caa"][1] == grid.xᶜᵃᵃ[1]
        @test ds["x_faa"][1] == grid.xᶠᵃᵃ[1]
        @test ds["y_aca"][1] == grid.yᵃᶜᵃ[1]
        @test ds["y_afa"][1] == grid.yᵃᶠᵃ[1]
        @test ds["z_aac"][1] == grid.z.cᵃᵃᶜ[1]
        @test ds["z_aaf"][1] == grid.z.cᵃᵃᶠ[1]

        @test ds["x_caa"][end] == grid.xᶜᵃᵃ[Nx]
        @test ds["x_faa"][end] == grid.xᶠᵃᵃ[Nx]
        @test ds["y_aca"][end] == grid.yᵃᶜᵃ[Ny]
        @test ds["y_afa"][end] == grid.yᵃᶠᵃ[Ny]
        @test ds["z_aac"][end] == grid.z.cᵃᵃᶜ[Nz]
        @test ds["z_aaf"][end] == grid.z.cᵃᵃᶠ[Nz+1]  # z is Bounded
    end

    @test eltype(ds["u"]) == FT
    @test eltype(ds["v"]) == FT
    @test eltype(ds["w"]) == FT
    @test eltype(ds["T"]) == FT
    @test eltype(ds["S"]) == FT

    u = ds["u"][:, :, :, end]
    v = ds["v"][:, :, :, end]
    w = ds["w"][:, :, :, end]
    T = ds["T"][:, :, :, end]
    S = ds["S"][:, :, :, end]

    close(ds)

    if with_halos
        @test all(u .≈ Array(model.velocities.u.data.parent))
        @test all(v .≈ Array(model.velocities.v.data.parent))
        @test all(w .≈ Array(model.velocities.w.data.parent))
        @test all(T .≈ Array(model.tracers.T.data.parent))
        @test all(S .≈ Array(model.tracers.S.data.parent))
    else
        @test all(u .≈ Array(interior(model.velocities.u)))
        @test all(v .≈ Array(interior(model.velocities.v)))
        @test all(w .≈ Array(interior(model.velocities.w)))
        @test all(T .≈ Array(interior(model.tracers.T)))
        @test all(S .≈ Array(interior(model.tracers.S)))
    end

    rm(nc_filepath)

    # Test sliced output only when not using halos
    if !with_halos
        ds2 = NCDataset(nc_sliced_filepath)

        @test haskey(ds2.attrib, "date")
        @test haskey(ds2.attrib, "Julia")
        @test haskey(ds2.attrib, "Oceananigans")
        @test haskey(ds2.attrib, "schedule")
        @test haskey(ds2.attrib, "interval")
        @test haskey(ds2.attrib, "output iteration interval")

        @test !isnothing(ds2.attrib["date"])
        @test !isnothing(ds2.attrib["Julia"])
        @test !isnothing(ds2.attrib["Oceananigans"])
        @test ds2.attrib["schedule"] == "IterationInterval"
        @test ds2.attrib["interval"] == 10
        @test !isnothing(ds2.attrib["output iteration interval"])

        # Dimensions should be Float64 by default
        @test eltype(ds2["time"])  == Float64
        @test eltype(ds2["x_caa"]) == Float64
        @test eltype(ds2["x_faa"]) == Float64
        @test eltype(ds2["y_aca"]) == Float64
        @test eltype(ds2["y_afa"]) == Float64
        @test eltype(ds2["z_aac"]) == Float64
        @test eltype(ds2["z_aaf"]) == Float64

        @test length(ds2["x_caa"]) == length(i_slice)
        @test length(ds2["x_faa"]) == length(i_slice)
        @test length(ds2["y_aca"]) == length(j_slice)
        @test length(ds2["y_afa"]) == length(j_slice)
        @test length(ds2["z_aac"]) == length(k_slice)
        @test length(ds2["z_aaf"]) == length(k_slice)

        @test ds2["x_caa"][1] == grid.xᶜᵃᵃ[i_slice[1]]
        @test ds2["x_faa"][1] == grid.xᶠᵃᵃ[i_slice[1]]
        @test ds2["y_aca"][1] == grid.yᵃᶜᵃ[j_slice[1]]
        @test ds2["y_afa"][1] == grid.yᵃᶠᵃ[j_slice[1]]
        @test ds2["z_aac"][1] == grid.z.cᵃᵃᶜ[k_slice[1]]
        @test ds2["z_aaf"][1] == grid.z.cᵃᵃᶠ[k_slice[1]]

        @test ds2["x_caa"][end] == grid.xᶜᵃᵃ[i_slice[end]]
        @test ds2["x_faa"][end] == grid.xᶠᵃᵃ[i_slice[end]]
        @test ds2["y_aca"][end] == grid.yᵃᶜᵃ[j_slice[end]]
        @test ds2["y_afa"][end] == grid.yᵃᶠᵃ[j_slice[end]]
        @test ds2["z_aac"][end] == grid.z.cᵃᵃᶜ[k_slice[end]]
        @test ds2["z_aaf"][end] == grid.z.cᵃᵃᶠ[k_slice[end]]

        @test eltype(ds2["u"]) == FT
        @test eltype(ds2["v"]) == FT
        @test eltype(ds2["w"]) == FT
        @test eltype(ds2["T"]) == FT
        @test eltype(ds2["S"]) == FT

        u_sliced = ds2["u"][:, :, :, end]
        v_sliced = ds2["v"][:, :, :, end]
        w_sliced = ds2["w"][:, :, :, end]
        T_sliced = ds2["T"][:, :, :, end]
        S_sliced = ds2["S"][:, :, :, end]

        close(ds2)

        @test all(u_sliced .≈ Array(interior(model.velocities.u))[i_slice, j_slice, k_slice])
        @test all(v_sliced .≈ Array(interior(model.velocities.v))[i_slice, j_slice, k_slice])
        @test all(w_sliced .≈ Array(interior(model.velocities.w))[i_slice, j_slice, k_slice])
        @test all(T_sliced .≈ Array(interior(model.tracers.T))[i_slice, j_slice, k_slice])
        @test all(S_sliced .≈ Array(interior(model.tracers.S))[i_slice, j_slice, k_slice])

        rm(nc_sliced_filepath)
    end

    return nothing
end

function test_netcdf_size_file_splitting(arch)
    grid = RectilinearGrid(arch,
                           size = (16, 16, 16),
                           extent = (1, 1, 1),
                           halo = (1, 1, 1))

    model = NonhydrostaticModel(grid;
                                buoyancy = SeawaterBuoyancy(),
                                tracers = (:T, :S))

    simulation = Simulation(model, Δt=1, stop_iteration=10)

    fake_attributes = Dict("fake_attribute" => "fake_attribute")

    Arch = typeof(arch)

    simulation.output_writers[:nc_writer] = NetCDFWriter(model, (; u=model.velocities.u);
        dir = ".",
        filename = "test_size_file_splitting_$Arch",
        schedule = IterationInterval(1),
        array_type = Array{Float64},
        with_halos = true,
        global_attributes = fake_attributes,
        file_splitting = FileSizeLimit(200KiB),
        overwrite_existing = true)

    # 531 KiB of output will be written which should get split into 3 files.
    run!(simulation)

    # Test that files has been split according to size as expected.
    @test filesize("test_size_file_splitting_$(Arch)_part1.nc") > 200KiB
    @test filesize("test_size_file_splitting_$(Arch)_part2.nc") > 200KiB
    @test filesize("test_size_file_splitting_$(Arch)_part3.nc") < 200KiB
    @test !isfile("test_size_file_splitting_$(Arch)_part4.nc")

    for n in string.(1:3)
        filename = "test_size_file_splitting_$(Arch)_part$n.nc"
        ds = NCDataset(filename, "r")

        # Test that all files contain the same dimensions.
        num_dims = length(keys(ds.dim))
        @test num_dims == 7

        # Test that all files contain the user defined attributes.
        @test ds.attrib["fake_attribute"] == "fake_attribute"

        close(ds)

        # Leave test directory clean.
        rm(filename)
    end

    return nothing
end

function test_netcdf_time_file_splitting(arch)
    grid = RectilinearGrid(arch,
                           size = (16, 16, 16),
                           extent = (1, 1, 1),
                           halo = (1, 1, 1))

    model = NonhydrostaticModel(grid;
                                buoyancy = SeawaterBuoyancy(),
                                tracers=(:T, :S))

    simulation = Simulation(model, Δt=1, stop_time=12seconds)

    fake_attributes = Dict("fake_attribute" => "fake_attribute")

    Arch = typeof(arch)
    simulation.output_writers[:nc_writer] = NetCDFWriter(model, (; u=model.velocities.u);
        dir = ".",
        filename = "test_time_file_splitting_$Arch",
        schedule = IterationInterval(2),
        array_type = Array{Float64},
        with_halos = true,
        global_attributes = fake_attributes,
        file_splitting = TimeInterval(4seconds),
        overwrite_existing = true)

    run!(simulation)

    for n in string.(1:3)
        filename = "test_time_file_splitting_$(Arch)_part$n.nc"
        ds = NCDataset(filename, "r")

        # Test that all files contain the same dimensions.
        num_dims = length(ds["time"])
        @test num_dims == 2

        # Test that all files contain the user defined attributes.
        @test ds.attrib["fake_attribute"] == "fake_attribute"

        close(ds)

        # Leave test directory clean.
        rm(filename)
    end

    # TODO: Why is there a part 4? We should add more tests.
    rm("test_time_file_splitting_$(Arch)_part4.nc")

    return nothing
end

function test_netcdf_function_output(arch)
    Nx = Ny = Nz = N = 16
    L = 1
    Δt = 1.25
    iters = 3

    grid = RectilinearGrid(arch,
                           size = (Nx, Ny, Nz),
                           extent = (L, 2L, 3L))

    model = NonhydrostaticModel(grid;
                                timestepper = :QuasiAdamsBashforth2,
                                buoyancy = SeawaterBuoyancy(),
                                tracers = (:T, :S))

    simulation = Simulation(model, Δt=Δt, stop_iteration=iters)

    # Define scalar, vector, and 2D slice outputs
    f(model) = model.clock.time^2

    g(model) = model.clock.time .* exp.(znodes(grid, Center()))

    xC = xnodes(grid, Center())
    yF = ynodes(grid, Face())

    XC = [xC[i] for i in 1:Nx, j in 1:Ny]
    YF = [yF[j] for i in 1:Nx, j in 1:Ny]

    h(model) = @. model.clock.time * sin(XC) * cos(YF) # xy slice output

    outputs = (scalar=f, profile=g, slice=h)
    dims = (scalar=(), profile=("z_aac",), slice=("x_caa", "y_aca"))

    output_attributes = (
        scalar = (long_name="Some scalar", units="bananas"),
        profile = (long_name="Some vertical profile", units="watermelons"),
        slice = (long_name="Some slice", units="mushrooms")
    )

    global_attributes = (
        location = "Bay of Fundy",
        onions = 7
    )

    Arch = typeof(arch)
    nc_filepath = "test_function_outputs_$Arch.nc"

    simulation.output_writers[:food] =
        NetCDFWriter(model, outputs;
            global_attributes,
            output_attributes,
            filename = nc_filepath,
            schedule = TimeInterval(Δt),
            dimensions = dims,
            array_type = Array{Float64},
            include_grid_metrics = false,
            verbose = true,
            overwrite_existing = true)

    run!(simulation)

    ds = NCDataset(nc_filepath, "r")

    @test haskey(ds.attrib, "date")
    @test haskey(ds.attrib, "Julia")
    @test haskey(ds.attrib, "Oceananigans")
    @test haskey(ds.attrib, "schedule")
    @test haskey(ds.attrib, "interval")
    @test haskey(ds.attrib, "output time interval")

    @test !isnothing(ds.attrib["date"])
    @test !isnothing(ds.attrib["Julia"])
    @test !isnothing(ds.attrib["Oceananigans"])
    @test !isnothing(ds.attrib["schedule"])
    @test !isnothing(ds.attrib["interval"])
    @test !isnothing(ds.attrib["output time interval"])

    @test eltype(ds["time"]) == Float64

    @test eltype(ds["x_caa"]) == Float64
    @test eltype(ds["x_faa"]) == Float64
    @test eltype(ds["y_aca"]) == Float64
    @test eltype(ds["y_afa"]) == Float64
    @test eltype(ds["z_aac"]) == Float64
    @test eltype(ds["z_aaf"]) == Float64

    @test length(ds["x_caa"]) == N
    @test length(ds["y_aca"]) == N
    @test length(ds["z_aac"]) == N
    @test length(ds["x_faa"]) == N
    @test length(ds["y_afa"]) == N
    @test length(ds["z_aaf"]) == N+1  # z is Bounded

    @test ds["x_caa"][1] == grid.xᶜᵃᵃ[1]
    @test ds["x_faa"][1] == grid.xᶠᵃᵃ[1]
    @test ds["y_aca"][1] == grid.yᵃᶜᵃ[1]
    @test ds["y_afa"][1] == grid.yᵃᶠᵃ[1]
    @test ds["z_aac"][1] == grid.z.cᵃᵃᶜ[1]
    @test ds["z_aaf"][1] == grid.z.cᵃᵃᶠ[1]

    @test ds["x_caa"][end] == grid.xᶜᵃᵃ[N]
    @test ds["y_aca"][end] == grid.yᵃᶜᵃ[N]
    @test ds["x_faa"][end] == grid.xᶠᵃᵃ[N]
    @test ds["y_afa"][end] == grid.yᵃᶠᵃ[N]
    @test ds["z_aac"][end] == grid.z.cᵃᵃᶜ[N]
    @test ds["z_aaf"][end] == grid.z.cᵃᵃᶠ[N+1]  # z is Bounded

    @test ds.attrib["location"] == "Bay of Fundy"
    @test ds.attrib["onions"] == 7

    @test eltype(ds["scalar"]) == Float64
    @test eltype(ds["profile"]) == Float64
    @test eltype(ds["slice"]) == Float64

    @test length(ds["time"]) == iters+1
    @test ds["time"][:] == [n*Δt for n in 0:iters]

    @test length(ds["scalar"]) == iters+1
    @test ds["scalar"].attrib["long_name"] == "Some scalar"
    @test ds["scalar"].attrib["units"] == "bananas"
    @test ds["scalar"][:] == [(n*Δt)^2 for n in 0:iters]
    @test dimnames(ds["scalar"]) == ("time",)

    @test ds["profile"].attrib["long_name"] == "Some vertical profile"
    @test ds["profile"].attrib["units"] == "watermelons"
    @test size(ds["profile"]) == (N, iters+1)
    @test dimnames(ds["profile"]) == ("z_aac", "time")

    for n in 0:iters
        @test ds["profile"][:, n+1] == n*Δt .* exp.(znodes(grid, Center()))
    end

    @test ds["slice"].attrib["long_name"] == "Some slice"
    @test ds["slice"].attrib["units"] == "mushrooms"
    @test size(ds["slice"]) == (N, N, iters+1)
    @test dimnames(ds["slice"]) == ("x_caa", "y_aca", "time")

    for n in 0:iters
        @test ds["slice"][:, :, n+1] == n*Δt .* sin.(XC) .* cos.(YF)
    end

    close(ds)

    #####
    ##### Take 1 more time step and test that appending to a NetCDF file works
    #####

    iters += 1
    simulation = Simulation(model, Δt=Δt, stop_iteration=iters)

    simulation.output_writers[:food] =
        NetCDFWriter(model, outputs;
            global_attributes,
            output_attributes,
            filename = nc_filepath,
            overwrite_existing = false,
            schedule = IterationInterval(1),
            array_type = Array{Float64},
            dimensions = dims,
            verbose = true)

    run!(simulation)

    ds = NCDataset(nc_filepath, "r")

    @test length(ds["time"]) == iters+1
    @test length(ds["scalar"]) == iters+1
    @test size(ds["profile"]) == (N, iters+1)
    @test size(ds["slice"]) == (N, N, iters+1)

    @test ds["time"][:] == [n*Δt for n in 0:iters]
    @test ds["scalar"][:] == [(n*Δt)^2 for n in 0:iters]

    for n in 0:iters
        @test ds["profile"][:, n+1] ≈ n*Δt .* exp.(znodes(grid, Center()))
        @test ds["slice"][:, :, n+1] ≈ n*Δt .* (sin.(XC) .* cos.(YF))
    end

    close(ds)

    rm(nc_filepath)

    return nothing
end

function test_netcdf_spatial_average(arch)
    topo = (Periodic, Periodic, Periodic)
    domain = (x=(0, 1), y=(0, 1), z=(0, 1))
    grid = RectilinearGrid(arch, topology=topo, size=(4, 4, 4); domain...)

    model = NonhydrostaticModel(grid; tracers = :c)

    set!(model, c=1)

    Δt = 0.01 # Floating point number chosen conservatively to flag rounding errors

    simulation = Simulation(model, Δt=Δt, stop_iteration=10)

    ∫c_dx = Field(Average(model.tracers.c, dims=(1)))
    ∫∫c_dxdy = Field(Average(model.tracers.c, dims=(1, 2)))
    ∫∫∫c_dxdydz = Field(Average(model.tracers.c, dims=(1, 2, 3)))

    Arch = typeof(arch)
    nc_filepath = "volume_averaged_field_test_$Arch.nc"

    simulation.output_writers[:averages] =
        NetCDFWriter(model,
            (; ∫c_dx, ∫∫c_dxdy, ∫∫∫c_dxdydz),
            array_type = Array{Float64},
            verbose = true,
            filename = nc_filepath,
            schedule = IterationInterval(2),
            include_grid_metrics = false)

    run!(simulation)

    ds = NCDataset(nc_filepath)

    for (n, t) in enumerate(ds["time"])
        @test all(ds["∫c_dx"][:,:, n] .≈ 1)
        @test all(ds["∫∫c_dxdy"][:, n] .≈ 1)
        @test all(ds["∫∫∫c_dxdydz"][n] .≈ 1)
    end

    close(ds)
    rm(nc_filepath)

    return nothing
end

function test_netcdf_time_averaging(arch)
    # Test for both "nice" floating point number and one that is more susceptible
    # to rounding errors
    for Δt in (1/64, 0.01)
        # Results should be very close (rtol < 1e-5) for stride = 1.
        # stride > 2 is currently not robust and can give inconsistent
        # results due to floating number errors that can result in vanishingly
        # small timesteps, which essentially decouples the clock time from
        # the iteration number.
        # Can add stride > 1 cases to the following line to test them.
        # for (stride, rtol) in zip((1, 2), (1e-5, 1e-3))
        for (stride, rtol) in zip((1), (1e-5))
            @info "  Testing time-averaging of NetCDF outputs [$(typeof(arch))] with " *
                  "timestep of $(Δt), stride of $(stride), and relative tolerance of $(rtol)."

            topo = (Periodic, Periodic, Periodic)
            domain = (x=(0, 1), y=(0, 1), z=(0, 1))
            grid = RectilinearGrid(arch, topology=topo, size=(4, 4, 4); domain...)

            λ1(x, y, z) = x + (1 - y)^2 + tanh(z)
            λ2(x, y, z) = x + (1 - y)^2 + tanh(4z)

            Fc1(x, y, z, t, c1) = - λ1(x, y, z) * c1
            Fc2(x, y, z, t, c2) = - λ2(x, y, z) * c2

            c1_forcing = Forcing(Fc1, field_dependencies=:c1)
            c2_forcing = Forcing(Fc2, field_dependencies=:c2)

            model = NonhydrostaticModel(grid;
                                        tracers = (:c1, :c2),
                                        forcing = (c1=c1_forcing, c2=c2_forcing))

            set!(model, c1=1, c2=1)

            # Floating point number chosen conservatively to flag rounding errors
            simulation = Simulation(model, Δt=Δt, stop_time=50Δt)

            ∫c1_dxdy = Field(Average(model.tracers.c1, dims=(1, 2)))
            ∫c2_dxdy = Field(Average(model.tracers.c2, dims=(1, 2)))

            nc_outputs = Dict(
                "c1" => ∫c1_dxdy,
                "c2" => ∫c2_dxdy
            )

            Arch = typeof(arch)
            horizontal_average_nc_filepath = "decay_averaged_field_test_$Arch.nc"

            simulation.output_writers[:horizontal_average] =
                NetCDFWriter(model,
                             nc_outputs,
                             array_type = Array{Float64},
                             verbose = true,
                             filename = horizontal_average_nc_filepath,
                             schedule = TimeInterval(10Δt),
                             include_grid_metrics = false,
                             overwrite_existing = true)

            multiple_time_average_nc_filepath = "decay_windowed_time_average_test_$Arch.nc"
            single_time_average_nc_filepath = "single_decay_windowed_time_average_test_$Arch.nc"
            window = 6Δt

            single_nc_output = Dict("c1" => ∫c1_dxdy)

            simulation.output_writers[:single_output_time_average] =
                NetCDFWriter(model,
                             single_nc_output,
                             array_type = Array{Float64},
                             verbose = true,
                             filename = single_time_average_nc_filepath,
                             schedule = AveragedTimeInterval(10Δt; window, stride),
                             include_grid_metrics = false,
                             overwrite_existing = true)

            simulation.output_writers[:multiple_output_time_average] =
                NetCDFWriter(model,
                             nc_outputs,
                             array_type = Array{Float64},
                             verbose = true,
                             filename = multiple_time_average_nc_filepath,
                             schedule = AveragedTimeInterval(10Δt; window, stride),
                             include_grid_metrics = false,
                             overwrite_existing = true)

            run!(simulation)

            ##### For each λ, the horizontal average should evaluate to
            #####
            #####     c̄(z, t) = ∫₀¹ ∫₀¹ exp{- λ(x, y, z) * t} dx dy
            #####             = 1 / (Nx*Ny) * Σᵢ₌₁ᴺˣ Σⱼ₌₁ᴺʸ exp{- λ(i, j, k) * t}
            #####
            ##### which we can compute analytically.

            ds = NCDataset(horizontal_average_nc_filepath)

            Nx, Ny, Nz = size(grid)
            xs, ys, zs = nodes(model.tracers.c1)

            c̄1(z, t) = 1 / (Nx * Ny) * sum(exp(-λ1(x, y, z) * t) for x in xs for y in ys)
            c̄2(z, t) = 1 / (Nx * Ny) * sum(exp(-λ2(x, y, z) * t) for x in xs for y in ys)

            for (n, t) in enumerate(ds["time"])
                @test all(isapprox.(ds["c1"][:, n], c̄1.(zs, t), rtol=rtol))
                @test all(isapprox.(ds["c2"][:, n], c̄2.(zs, t), rtol=rtol))
            end

            close(ds)

            # Compute time averages...
            c̄1(ts) = 1/length(ts) * sum(c̄1.(zs, t) for t in ts)
            c̄2(ts) = 1/length(ts) * sum(c̄2.(zs, t) for t in ts)

            #####
            ##### Test strided windowed time average against analytic solution
            ##### for *single* NetCDF output
            #####

            single_ds = NCDataset(single_time_average_nc_filepath)

            attribute_names = (
                "schedule", "interval", "output time interval",
                "time_averaging_window", "time averaging window",
                "time_averaging_stride", "time averaging stride"
            )

            for name in attribute_names
                @test haskey(single_ds.attrib, name)
                @test !isnothing(single_ds.attrib[name])
            end

            window_size = Int(window/Δt)

            @info "    Testing time-averaging of a single NetCDF output [$(typeof(arch))]..."

            for (n, t) in enumerate(single_ds["time"][2:end])
                averaging_times = [t - n*Δt for n in 0:stride:window_size-1 if t - n*Δt >= 0]
                @test all(isapprox.(single_ds["c1"][:, n+1], c̄1(averaging_times), rtol=rtol, atol=rtol))
            end

            close(single_ds)

            #####
            ##### Test strided windowed time average against analytic solution
            ##### for *multiple* NetCDF outputs
            #####

            ds = NCDataset(multiple_time_average_nc_filepath)

            @info "    Testing time-averaging of multiple NetCDF outputs [$(typeof(arch))]..."

            for (n, t) in enumerate(ds["time"][2:end])
                averaging_times = [t - n*Δt for n in 0:stride:window_size-1 if t - n*Δt >= 0]
                @test all(isapprox.(ds["c2"][:, n+1], c̄2(averaging_times), rtol=rtol))
            end

            close(ds)

            rm(horizontal_average_nc_filepath)
            rm(single_time_average_nc_filepath)
            rm(multiple_time_average_nc_filepath)
        end
    end

    return nothing
end

function test_netcdf_output_alignment(arch)
    grid = RectilinearGrid(size=(1, 1, 1), extent=(1, 1, 1))

    model = NonhydrostaticModel(grid;
                                timestepper = :QuasiAdamsBashforth2,
                                buoyancy = SeawaterBuoyancy(),
                                tracers = (:T, :S))

    simulation = Simulation(model, Δt=0.2, stop_time=40)

    Arch = typeof(arch)
    test_filename1 = "test_output_alignment1_$Arch.nc"
    simulation.output_writers[:stuff] =
        NetCDFWriter(model,
            model.velocities,
            filename = test_filename1,
            schedule = TimeInterval(7.3),
            include_grid_metrics = false)

    test_filename2 = "test_output_alignment2_$Arch.nc"
    simulation.output_writers[:something] =
        NetCDFWriter(model,
            model.tracers,
            filename = test_filename2,
            schedule = TimeInterval(3.0),
            include_grid_metrics = false)

    run!(simulation)

    NCDataset(test_filename1, "r") do ds
        @test all(ds["time"] .== 0:7.3:40)
    end

    NCDataset(test_filename2, "r") do ds
        @test all(ds["time"] .== 0:3.0:40)
    end

    rm(test_filename1)
    rm(test_filename2)

    return nothing
end

function test_netcdf_output_just_particles(arch)
    grid = RectilinearGrid(arch;
                           topology = (Periodic, Periodic, Bounded),
                           size = (5, 5, 5),
                           x = (-1, 1),
                           y = (-1, 1),
                           z = (-1, 0))

    Np = 10
    xs = on_architecture(arch, 0.25 * ones(Np))
    ys = on_architecture(arch, -0.4 * ones(Np))
    zs = on_architecture(arch, -0.12 * ones(Np))

    particles = LagrangianParticles(x=xs, y=ys, z=zs)

    model = NonhydrostaticModel(grid; particles)

    set!(model, u=1, v=1)

    Nt = 10
    simulation = Simulation(model, Δt=1e-2, stop_iteration=Nt)

    Arch = typeof(arch)
    filepath = "test_just_particles_$Arch.nc"
    simulation.output_writers[:particles_nc] =
        NetCDFWriter(model,
            (; particles = model.particles),
            filename = filepath,
            schedule = IterationInterval(1),
            include_grid_metrics = false)

    run!(simulation)

    ds = NCDataset(filepath)

    @test haskey(ds, "time")
    @test length(ds[:time]) == Nt + 1

    @test haskey(ds, "particle_id")
    @test length(ds[:particle_id]) == Np

    @test haskey(ds, "x")
    @test haskey(ds, "y")
    @test haskey(ds, "z")

    @test size(ds[:x]) == (Np, Nt+1)
    @test size(ds[:y]) == (Np, Nt+1)
    @test size(ds[:z]) == (Np, Nt+1)

    close(ds)
    rm(filepath)

    return nothing
end

function test_netcdf_output_particles_and_fields(arch)
    N = 5

    grid = RectilinearGrid(arch;
                           topology = (Periodic, Periodic, Bounded),
                           size = (N, N, N),
                           x = (-1, 1),
                           y = (-1, 1),
                           z = (-1, 0))

    Np = 10
    xs = on_architecture(arch, 0.25 * ones(Np))
    ys = on_architecture(arch, -0.4 * ones(Np))
    zs = on_architecture(arch, -0.12 * ones(Np))

    particles = LagrangianParticles(x=xs, y=ys, z=zs)

    model = NonhydrostaticModel(grid; particles)

    set!(model, u=1, v=1)

    Nt = 10
    simulation = Simulation(model, Δt=1e-2, stop_iteration=Nt)

    outputs = (particles = model.particles,
               u = model.velocities.u,
               v = model.velocities.v,
               w = model.velocities.w)

    Arch = typeof(arch)
    filepath = "test_particles_and_fields_$Arch.nc"
    simulation.output_writers[:particles_nc] =
        NetCDFWriter(model, outputs,
            filename = filepath,
            schedule = IterationInterval(1),
            include_grid_metrics = false)

    run!(simulation)

    ds = NCDataset(filepath)

    @test haskey(ds, "time")
    @test length(ds[:time]) == Nt + 1

    @test haskey(ds, "particle_id")
    @test length(ds[:particle_id]) == Np

    @test haskey(ds, "x")
    @test haskey(ds, "y")
    @test haskey(ds, "z")

    @test size(ds[:x]) == (Np, Nt+1)
    @test size(ds[:y]) == (Np, Nt+1)
    @test size(ds[:z]) == (Np, Nt+1)

    @test haskey(ds, "u")
    @test haskey(ds, "v")
    @test haskey(ds, "w")

    @test size(ds[:u]) == (N, N, N,   Nt+1)
    @test size(ds[:v]) == (N, N, N,   Nt+1)
    @test size(ds[:w]) == (N, N, N+1, Nt+1)

    close(ds)
    rm(filepath)

    return nothing
end

function test_netcdf_vertically_stretched_grid_output(arch)
    Nx = Ny = 8
    Nz = 16

    z_faces = [k^2 for k in 0:Nz]

    grid = RectilinearGrid(arch;
                           size = (Nx, Ny, Nz),
                           x = (0, 1),
                           y = (-π, π),
                           z = z_faces)

    model = NonhydrostaticModel(grid;
                                buoyancy = SeawaterBuoyancy(),
                                tracers = (:T, :S))

    simulation = Simulation(model, Δt=1.25, stop_iteration=3)

    Arch = typeof(arch)
    nc_filepath = "test_netcdf_vertically_stretched_grid_output_$Arch.nc"

    simulation.output_writers[:fields] =
        NetCDFWriter(model,
            merge(model.velocities, model.tracers),
            filename = nc_filepath,
            schedule = IterationInterval(1),
            array_type = Array{Float64},
            include_grid_metrics = false,
            verbose = true)

    run!(simulation)

    ds = NCDataset(nc_filepath)

    @test length(ds["x_caa"]) == Nx
    @test length(ds["y_aca"]) == Ny
    @test length(ds["z_aac"]) == Nz
    @test length(ds["x_faa"]) == Nx
    @test length(ds["y_afa"]) == Ny
    @test length(ds["z_aaf"]) == Nz+1  # z is Bounded

    @test ds["x_caa"][1] == grid.xᶜᵃᵃ[1]
    @test ds["x_faa"][1] == grid.xᶠᵃᵃ[1]
    @test ds["y_aca"][1] == grid.yᵃᶜᵃ[1]
    @test ds["y_afa"][1] == grid.yᵃᶠᵃ[1]

    @test @allowscalar ds["z_aac"][1] == grid.z.cᵃᵃᶜ[1]
    @test @allowscalar ds["z_aaf"][1] == grid.z.cᵃᵃᶠ[1]

    @test ds["x_caa"][end] == grid.xᶜᵃᵃ[Nx]
    @test ds["x_faa"][end] == grid.xᶠᵃᵃ[Nx]
    @test ds["y_aca"][end] == grid.yᵃᶜᵃ[Ny]
    @test ds["y_afa"][end] == grid.yᵃᶠᵃ[Ny]

    @test @allowscalar ds["z_aac"][end] == grid.z.cᵃᵃᶜ[Nz]
    @test @allowscalar ds["z_aaf"][end] == grid.z.cᵃᵃᶠ[Nz+1]  # z is Bounded

    close(ds)
    rm(nc_filepath)

    return nothing
end

function test_netcdf_overriding_attributes(arch)
    arch = CPU()

    grid = LatitudeLongitudeGrid(arch;
                                 topology = (Bounded, Bounded, Bounded),
                                 size = (4, 4, 4),
                                 longitude = (-1, 1),
                                 latitude = (-1, 1),
                                 z = (-100, 0))

    model = HydrostaticFreeSurfaceModel(grid;
                                        closure = ScalarDiffusivity(ν=4e-2, κ=4e-2),
                                        buoyancy = SeawaterBuoyancy(),
                                        tracers = (:T, :S))

    outputs = merge(model.velocities, model.tracers)

    Arch = typeof(arch)
    nc_filepath = "test_attributes_$Arch.nc"
    isfile(nc_filepath) && rm(nc_filepath)

    global_attributes = Dict(
        "date" => "yesterday", # This should override the default date attribute.
        "fruit" => "papaya"    # This should show up in the NetCDF file.
    )

    output_attributes = Dict(
        "φ_afa" => Dict("units" => "No units for you!"),
        "u" => Dict("long_name" => "zonal velocity", "units" => "miles/fortnight")
    )

    nc_writer = NetCDFWriter(model, outputs;
        filename = nc_filepath,
        schedule = IterationInterval(1),
        global_attributes,
        output_attributes)

    ds = NCDataset(nc_filepath)

    @test ds.attrib["date"] == "yesterday"
    @test ds.attrib["fruit"] == "papaya"

    @test ds["φ_afa"].attrib["units"] == "No units for you!"
    @test !haskey(ds["φ_afa"].attrib, "long_name")

    @test ds["u"].attrib["long_name"] == "zonal velocity"
    @test ds["u"].attrib["units"] == "miles/fortnight"

    return nothing
end

function test_netcdf_hydrostatic_free_surface_only_output(arch; immersed=false, vertically_stretched=false)
    Nλ, Nφ, Nz = 8, 8, 4
    Hλ, Hφ, Hz = 3, 4, 2

    z = vertically_stretched ? [k^2 - 100 for k in 0:Nz] : (-100, 0)

    underlying_grid = LatitudeLongitudeGrid(arch;
                                            topology = (Bounded, Bounded, Bounded),
                                            size = (Nλ, Nφ, Nz),
                                            halo = (Hλ, Hφ, Hz),
                                            longitude = (-1, 1),
                                            latitude = (-1, 1),
                                            z)

    grid = immersed ? ImmersedBoundaryGrid(underlying_grid, GridFittedBottom(-50)) : underlying_grid

    model = HydrostaticFreeSurfaceModel(grid;
                                        closure = ScalarDiffusivity(ν=4e-2, κ=4e-2),
                                        buoyancy = SeawaterBuoyancy(),
                                        tracers = (:T, :S))

    Nt = 5
    simulation = Simulation(model, Δt=0.1, stop_iteration=Nt)

    outputs = (; model.free_surface.displacement,)

    Arch = typeof(arch)
    immersed_str = immersed ? "_immersed" : ""
    stretched_str = vertically_stretched ? "_stretched" : ""
    filepath_with_halos = "test_free_surface_with_halos_$(Arch)$(immersed_str)$(stretched_str).nc"
    isfile(filepath_with_halos) && rm(filepath_with_halos)

    simulation.output_writers[:with_halos] =
        NetCDFWriter(model, outputs;
            filename = filepath_with_halos,
            schedule = IterationInterval(1),
            with_halos = true)

    filepath_no_halos = "test_free_surface_no_halos_$(Arch)$(immersed_str)$(stretched_str).nc"
    isfile(filepath_no_halos) && rm(filepath_no_halos)

    simulation.output_writers[:no_halos] =
        NetCDFWriter(model, outputs;
            filename = filepath_no_halos,
            schedule = IterationInterval(1),
            with_halos = false)

    run!(simulation)

    ds_h = NCDataset(filepath_with_halos)

    @test haskey(ds_h, "displacement")
    @test dimsize(ds_h["displacement"]) == (λ_caa=Nλ + 2Hλ, φ_aca=Nφ + 2Hφ, z_aaf_displacement=1, time=Nt + 1)

    close(ds_h)
    rm(filepath_with_halos)

    ds_n = NCDataset(filepath_no_halos)

    @test haskey(ds_n, "displacement")
    @test dimsize(ds_n["displacement"]) == (λ_caa=Nλ, φ_aca=Nφ, z_aaf_displacement=1, time=Nt + 1)

    close(ds_n)
    rm(filepath_no_halos)

    return nothing
end

function test_netcdf_hydrostatic_free_surface_mixed_output(arch; immersed=false, vertically_stretched=false)
    Nλ, Nφ, Nz = 8, 8, 4
    Hλ, Hφ, Hz = 3, 4, 2

    z = vertically_stretched ? [k^2 - 100 for k in 0:Nz] : (-100, 0)

    underlying_grid = LatitudeLongitudeGrid(arch;
                                            topology = (Bounded, Bounded, Bounded),
                                            size = (Nλ, Nφ, Nz),
                                            halo = (Hλ, Hφ, Hz),
                                            longitude = (-1, 1),
                                            latitude = (-1, 1),
                                            z)

    grid = immersed ? ImmersedBoundaryGrid(underlying_grid, GridFittedBottom(-50)) : underlying_grid

    model = HydrostaticFreeSurfaceModel(grid;
                                        closure = ScalarDiffusivity(ν=4e-2, κ=4e-2),
                                        buoyancy = SeawaterBuoyancy(),
                                        tracers = (:T, :S))

    Nt = 5
    simulation = Simulation(model, Δt=0.1, stop_iteration=Nt)

    free_surface_outputs = (; model.free_surface.displacement,)

    outputs = merge(model.velocities, model.tracers, free_surface_outputs)

    Arch = typeof(arch)
    immersed_str = immersed ? "_immersed" : ""
    stretched_str = vertically_stretched ? "_stretched" : ""
    filepath_with_halos = "test_mixed_free_surface_with_halos_$(Arch)$(immersed_str)$(stretched_str).nc"
    isfile(filepath_with_halos) && rm(filepath_with_halos)

    simulation.output_writers[:with_halos] =
        NetCDFWriter(model, outputs;
            filename = filepath_with_halos,
            schedule = IterationInterval(1),
            with_halos = true,
            overwrite_existing = true)

    filepath_no_halos = "test_mixed_free_surface_no_halos_$(Arch)$(immersed_str)$(stretched_str).nc"
    isfile(filepath_no_halos) && rm(filepath_no_halos)

    simulation.output_writers[:no_halos] =
        NetCDFWriter(model, outputs;
            filename = filepath_no_halos,
            schedule = IterationInterval(1),
            with_halos = false,
            overwrite_existing = true)

    run!(simulation)

    ds_h = NCDataset(filepath_with_halos)

    @test haskey(ds_h, "displacement")
    @test dimsize(ds_h["displacement"]) == (λ_caa=Nλ + 2Hλ, φ_aca=Nφ + 2Hφ, z_aaf_displacement=1, time=Nt + 1)

    @test dimsize(ds_h[:u]) == (λ_faa=Nλ + 2Hλ + 1, φ_aca=Nφ + 2Hφ,     z_aac=Nz + 2Hz,     time=Nt + 1)
    @test dimsize(ds_h[:v]) == (λ_caa=Nλ + 2Hλ,     φ_afa=Nφ + 2Hφ + 1, z_aac=Nz + 2Hz,     time=Nt + 1)
    @test dimsize(ds_h[:w]) == (λ_caa=Nλ + 2Hλ,     φ_aca=Nφ + 2Hφ,     z_aaf=Nz + 2Hz + 1, time=Nt + 1)
    @test dimsize(ds_h[:T]) == (λ_caa=Nλ + 2Hλ,     φ_aca=Nφ + 2Hφ,     z_aac=Nz + 2Hz,     time=Nt + 1)
    @test dimsize(ds_h[:S]) == (λ_caa=Nλ + 2Hλ,     φ_aca=Nφ + 2Hφ,     z_aac=Nz + 2Hz,     time=Nt + 1)

    close(ds_h)
    rm(filepath_with_halos)

    ds_n = NCDataset(filepath_no_halos)

    @test haskey(ds_n, "displacement")
    @test dimsize(ds_n["displacement"]) == (λ_caa=Nλ, φ_aca=Nφ, z_aaf_displacement=1, time=Nt + 1)

    @test dimsize(ds_n[:u]) == (λ_faa=Nλ + 1, φ_aca=Nφ,     z_aac=Nz,     time=Nt + 1)
    @test dimsize(ds_n[:v]) == (λ_caa=Nλ,     φ_afa=Nφ + 1, z_aac=Nz,     time=Nt + 1)
    @test dimsize(ds_n[:w]) == (λ_caa=Nλ,     φ_aca=Nφ,     z_aaf=Nz + 1, time=Nt + 1)
    @test dimsize(ds_n[:T]) == (λ_caa=Nλ,     φ_aca=Nφ,     z_aac=Nz,     time=Nt + 1)
    @test dimsize(ds_n[:S]) == (λ_caa=Nλ,     φ_aca=Nφ,     z_aac=Nz,     time=Nt + 1)

    close(ds_n)
    rm(filepath_no_halos)

    return nothing
end

function test_netcdf_nonhydrostatic_free_surface_only_output(arch; immersed=false, vertically_stretched=false)
    Nx, Ny, Nz = 8, 8, 4
    Hx, Hy, Hz = 3, 4, 2

    z = vertically_stretched ? [k^2 - 100 for k in 0:Nz] : (-100, 0)

    underlying_grid = RectilinearGrid(arch;
                                      topology = (Bounded, Bounded, Bounded),
                                      size = (Nx, Ny, Nz),
                                      halo = (Hx, Hy, Hz),
                                      x = (-1, 1),
                                      y = (-1, 1),
                                      z)

    grid = immersed ? ImmersedBoundaryGrid(underlying_grid, GridFittedBottom(-50)) : underlying_grid

    model = NonhydrostaticModel(grid;
                                free_surface = ImplicitFreeSurface(),
                                closure = ScalarDiffusivity(ν=4e-2, κ=4e-2),
                                buoyancy = SeawaterBuoyancy(),
                                tracers = (:T, :S))

    Nt = 5
    simulation = Simulation(model, Δt=0.1, stop_iteration=Nt)

    outputs = (; model.free_surface.displacement)

    Arch = typeof(arch)
    immersed_str = immersed ? "_immersed" : ""
    stretched_str = vertically_stretched ? "_stretched" : ""
    filepath_with_halos = "test_nonhydrostatic_free_surface_with_halos_$(Arch)$(immersed_str)$(stretched_str).nc"
    isfile(filepath_with_halos) && rm(filepath_with_halos)

    simulation.output_writers[:with_halos] =
        NetCDFWriter(model, outputs;
            filename = filepath_with_halos,
            schedule = IterationInterval(1),
            with_halos = true)

    filepath_no_halos = "test_nonhydrostatic_free_surface_no_halos_$(Arch)$(immersed_str)$(stretched_str).nc"
    isfile(filepath_no_halos) && rm(filepath_no_halos)

    simulation.output_writers[:no_halos] =
        NetCDFWriter(model, outputs;
            filename = filepath_no_halos,
            schedule = IterationInterval(1),
            with_halos = false)

    run!(simulation)

    ds_h = NCDataset(filepath_with_halos)
    @test haskey(ds_h, "displacement")
    @test dimsize(ds_h["displacement"]) == (x_caa=Nx + 2Hx, y_aca=Ny + 2Hy, z_aaf_displacement=1, time=Nt + 1)
    close(ds_h)
    rm(filepath_with_halos)

    ds_n = NCDataset(filepath_no_halos)
    @test haskey(ds_n, "displacement")
    @test dimsize(ds_n["displacement"]) == (x_caa=Nx, y_aca=Ny, z_aaf_displacement=1, time=Nt + 1)
    close(ds_n)
    rm(filepath_no_halos)

    return nothing
end

function test_netcdf_nonhydrostatic_free_surface_mixed_output(arch; immersed=false, vertically_stretched=false)
    Nx, Ny, Nz = 8, 8, 4
    Hx, Hy, Hz = 3, 4, 2

    z = vertically_stretched ? [k^2 - 100 for k in 0:Nz] : (-100, 0)

    underlying_grid = RectilinearGrid(arch;
                                      topology = (Bounded, Bounded, Bounded),
                                      size = (Nx, Ny, Nz),
                                      halo = (Hx, Hy, Hz),
                                      x = (-1, 1),
                                      y = (-1, 1),
                                      z)

    grid = immersed ? ImmersedBoundaryGrid(underlying_grid, GridFittedBottom(-50)) : underlying_grid

    model = NonhydrostaticModel(grid;
                                free_surface = ImplicitFreeSurface(),
                                closure = ScalarDiffusivity(ν=4e-2, κ=4e-2),
                                buoyancy = SeawaterBuoyancy(),
                                tracers = (:T, :S))

    Nt = 5
    simulation = Simulation(model, Δt=0.1, stop_iteration=Nt)

    free_surface_outputs = (; model.free_surface.displacement)
    outputs = merge(model.velocities, model.tracers, free_surface_outputs)

    Arch = typeof(arch)
    immersed_str = immersed ? "_immersed" : ""
    stretched_str = vertically_stretched ? "_stretched" : ""
    filepath_with_halos = "test_nonhydrostatic_mixed_free_surface_with_halos_$(Arch)$(immersed_str)$(stretched_str).nc"
    isfile(filepath_with_halos) && rm(filepath_with_halos)

    simulation.output_writers[:with_halos] =
        NetCDFWriter(model, outputs;
            filename = filepath_with_halos,
            schedule = IterationInterval(1),
            with_halos = true,
            overwrite_existing = true)

    filepath_no_halos = "test_nonhydrostatic_mixed_free_surface_no_halos_$(Arch)$(immersed_str)$(stretched_str).nc"
    isfile(filepath_no_halos) && rm(filepath_no_halos)

    simulation.output_writers[:no_halos] =
        NetCDFWriter(model, outputs;
            filename = filepath_no_halos,
            schedule = IterationInterval(1),
            with_halos = false,
            overwrite_existing = true)

    run!(simulation)

    ds_h = NCDataset(filepath_with_halos)

    @test haskey(ds_h, "displacement")
    @test dimsize(ds_h["displacement"]) == (x_caa=Nx + 2Hx, y_aca=Ny + 2Hy, z_aaf_displacement=1, time=Nt + 1)

    @test dimsize(ds_h[:u]) == (x_faa=Nx + 2Hx + 1, y_aca=Ny + 2Hy,     z_aac=Nz + 2Hz,     time=Nt + 1)
    @test dimsize(ds_h[:v]) == (x_caa=Nx + 2Hx,     y_afa=Ny + 2Hy + 1, z_aac=Nz + 2Hz,     time=Nt + 1)
    @test dimsize(ds_h[:w]) == (x_caa=Nx + 2Hx,     y_aca=Ny + 2Hy,     z_aaf=Nz + 2Hz + 1, time=Nt + 1)
    @test dimsize(ds_h[:T]) == (x_caa=Nx + 2Hx,     y_aca=Ny + 2Hy,     z_aac=Nz + 2Hz,     time=Nt + 1)
    @test dimsize(ds_h[:S]) == (x_caa=Nx + 2Hx,     y_aca=Ny + 2Hy,     z_aac=Nz + 2Hz,     time=Nt + 1)

    close(ds_h)
    rm(filepath_with_halos)

    ds_n = NCDataset(filepath_no_halos)

    @test haskey(ds_n, "displacement")
    @test dimsize(ds_n["displacement"]) == (x_caa=Nx, y_aca=Ny, z_aaf_displacement=1, time=Nt + 1)

    @test dimsize(ds_n[:u]) == (x_faa=Nx + 1, y_aca=Ny,     z_aac=Nz,     time=Nt + 1)
    @test dimsize(ds_n[:v]) == (x_caa=Nx,     y_afa=Ny + 1, z_aac=Nz,     time=Nt + 1)
    @test dimsize(ds_n[:w]) == (x_caa=Nx,     y_aca=Ny,     z_aaf=Nz + 1, time=Nt + 1)
    @test dimsize(ds_n[:T]) == (x_caa=Nx,     y_aca=Ny,     z_aac=Nz,     time=Nt + 1)
    @test dimsize(ds_n[:S]) == (x_caa=Nx,     y_aca=Ny,     z_aac=Nz,     time=Nt + 1)

    close(ds_n)
    rm(filepath_no_halos)

    return nothing
end

function test_netcdf_buoyancy_force(arch)
    Nx, Nz = 8, 8
    Hx, Hz = 2, 3
    Lx, H  = 2, 1

    grid = RectilinearGrid(arch,
                           topology = (Periodic, Flat, Bounded),
                           size = (Nx, Nz),
                           halo = (Hx, Hz),
                           x = (-Lx, Lx),
                           z = (-H, 0))

    Boussinesq_eos = (TEOS10EquationOfState(),
                      RoquetEquationOfState(:Linear),
                      RoquetEquationOfState(:Cabbeling),
                      RoquetEquationOfState(:CabbelingThermobaricity),
                      RoquetEquationOfState(:Freezing),
                      RoquetEquationOfState(:SecondOrder),
                      RoquetEquationOfState(:SimplestRealistic))

    for eos in Boussinesq_eos
        model = NonhydrostaticModel(grid;
                                    closure = ScalarDiffusivity(ν=4e-2, κ=4e-2),
                                    buoyancy = SeawaterBuoyancy(equation_of_state=eos),
                                    tracers = (:T, :S))

        Nt = 7
        simulation = Simulation(model, Δt=0.1, stop_iteration=Nt)

        simulation.output_writers[:b_eos] = NetCDFWriter(model, fields(model),
                                                         filename = string(eos)*"_.nc",
                                                         schedule = IterationInterval(1),
                                                         array_type = Array{Float64},
                                                         include_grid_metrics = true,
                                                         verbose = true)
        # only tests that the writer builds, produces a file at filepath and sets attributes
        @test simulation.output_writers[:b_eos] isa NetCDFWriter
        @test isfile(simulation.output_writers[:b_eos].filepath)
        ds = NCDataset(simulation.output_writers[:b_eos].filepath)
        @test ds["T"].attrib["long_name"] == "Conservative temperature"
        @test ds["T"].attrib["units"] == "°C"
        @test ds["S"].attrib["long_name"] == "Absolute salinity"
        @test ds["S"].attrib["units"] == "g/kg"
        close(ds)
        rm(simulation.output_writers[:b_eos].filepath)
    end
    return nothing
end

function test_netcdf_single_field_defvar(grid; immersed=false)
    grid = immersed ? ImmersedBoundaryGrid(grid, GridFittedBottom(-1/2)) : grid

    c = CenterField(grid)
    set!(c, (x, y, z) -> x + y + z) # Make sure field is non-trivial

    filepath = "test_single_field_defvar.nc"
    isfile(filepath) && rm(filepath)

    ds = NCDataset(filepath, "c")
    c_long_name = "Center field"
    c_units = "kg/m³"
    defVar(ds, "c", c, attrib=Dict("long_name" => c_long_name, "units" => c_units))

    # Write an abstract operation
    c² = c^2
    c²_long_name = "Center field squared"
    defVar(ds, "c²", c², attrib=Dict("long_name" => c²_long_name))

    # Write a reduction
    c̄ = Average(c, dims=3)
    c̄_long_name = "Average center field"
    c̄_units = "kg/m³"
    defVar(ds, "c̄", c̄, attrib=Dict("long_name" => c̄_long_name, "units" => c̄_units))

    close(ds)
    ds = NCDataset(filepath, "r")

    @test "c" ∈ keys(ds)
    @test ds["c"] == Array(interior(c))
    @test ds["c"].attrib["long_name"] == c_long_name
    @test ds["c"].attrib["units"] == c_units

    @test "c²" ∈ keys(ds)
    @test ds["c²"] == Array(interior(Field(c²)))
    @test ds["c²"].attrib["long_name"] == c²_long_name

    @test "c̄" ∈ keys(ds)
    @test all(ds["c̄"] .== Array(interior(Field(c̄))))
    @test ds["c̄"].attrib["long_name"] == c̄_long_name
    @test ds["c̄"].attrib["units"] == c̄_units

    close(ds)
    rm(filepath)

    return nothing
end

function test_netcdf_field_dimension_validation(grid; immersed=false)
    grid = immersed ? ImmersedBoundaryGrid(grid, GridFittedBottom(-1/2)) : grid
    c = CenterField(grid)

    filepath = "test_dimension_validation_success.nc"
    isfile(filepath) && rm(filepath)
    ds = NCDataset(filepath, "c")

    # Get dimension names in a way that works for multiple grid types
    ξ_center_name = string(ξname(grid)) * "_caa"
    η_center_name = string(ηname(grid)) * "_aca"
    r_center_name = string(rname(grid)) * "_aac"

    # Define dimensions first
    defVar(ds, ξ_center_name, ξnodes(grid, Center()), (ξ_center_name,))
    # We "forget" to define the η dimension. Should work regardless
    defVar(ds, r_center_name, rnodes(grid, Center()), (r_center_name,))

    # Write variable to disk after dimensions are created and match the
    defVar(ds, "c", c, time_dependent=false)
    @test ds[ξ_center_name][:] == ξnodes(grid, Center())
    @test ds[η_center_name][:] == ηnodes(grid, Center())
    @test ds[r_center_name][:] == rnodes(grid, Center())

    close(ds)
    rm(filepath)

    # Test that wrong dimension should throw error
    filepath_wrong_size = "test_dimension_validation_wrong_size.nc"
    isfile(filepath_wrong_size) && rm(filepath_wrong_size)

    ds_wrong = NCDataset(filepath_wrong_size, "c")

    # Define dimensions with wrong sizes
    defVar(ds_wrong, ξ_center_name, ξnodes(grid, Center()), (ξ_center_name,))
    defVar(ds_wrong, η_center_name, ηnodes(grid, Center()), (η_center_name,))
    defVar(ds_wrong, r_center_name, rnodes(grid, Face()), (r_center_name,)) # ⚠ wrongly write faces instead of center

    @test_throws ArgumentError defVar(ds_wrong, "c", c)

    close(ds_wrong)
    rm(filepath_wrong_size)

    return nothing
end

function test_netcdf_multiple_grids_defvar(grid1, grid2; immersed=false)
    if immersed
        grid1 = ImmersedBoundaryGrid(grid1, GridFittedBottom(-1/2))
        grid2 = ImmersedBoundaryGrid(grid2, GridFittedBottom(-0.8))
    end

    c1 = CenterField(grid1)
    c2 = CenterField(grid2)

    # Set different values for each field
    set!(c1, (x, y, z) -> x + y + z)
    set!(c2, (x, y, z) -> 2 * (x^2 + y^2 + z^2))

    filepath = "test_multiple_grids_defvar.nc"
    isfile(filepath) && rm(filepath)
    ds = NCDataset(filepath, "c")

    # Define variables from grid1
    suffixed_dim_name_grid1(args...) = trilocation_dim_name(args...) * "_grid1"
    defVar(ds, "c1", c1, attrib=Dict("long_name" => "Field from grid 1"), dimension_name_generator=suffixed_dim_name_grid1)

    # Define variables from grid2 - this should create new dimension names
    suffixed_dim_name_grid2(args...) = trilocation_dim_name(args...) * "_grid2"
    defVar(ds, "c2", c2, attrib=Dict("long_name" => "Field from grid 2"), dimension_name_generator=suffixed_dim_name_grid2)

    close(ds)

    # Verify the file was created correctly
    ds = NCDataset(filepath, "r")
    @test "c1" ∈ keys(ds)
    @test "c2" ∈ keys(ds)

    @test size(ds["c1"]) == size(grid1)
    @test size(ds["c2"]) == size(grid2)

    # Check that the coordinate variables were created with appropriate names
    # Get dimension names in a way that works for multiple grid types
    ξ_center_name1 = string(ξname(grid1)) * "_caa_grid1"
    η_center_name1 = string(ηname(grid1)) * "_aca_grid1"
    r_center_name1 = string(rname(grid1)) * "_aac_grid1"

    @test ξ_center_name1 ∈ keys(ds)
    @test η_center_name1 ∈ keys(ds)
    @test r_center_name1 ∈ keys(ds)

    ξ_center_name2 = string(ξname(grid2)) * "_caa_grid2"
    η_center_name2 = string(ηname(grid2)) * "_aca_grid2"
    r_center_name2 = string(rname(grid2)) * "_aac_grid2"

    @test ξ_center_name2 ∈ keys(ds)
    @test η_center_name2 ∈ keys(ds)
    @test r_center_name2 ∈ keys(ds)

    # Check that data values match
    @test ds["c1"] == Array(interior(c1))
    @test ds["c2"] == Array(interior(c2))

    close(ds)
    rm(filepath)

    return nothing
end

function test_netcdf_writer_different_grid(arch)

    grid = RectilinearGrid(arch, size=(1, 1, 8), extent=(1, 1, 1))
    model = NonhydrostaticModel(grid)

    coarse_grid = RectilinearGrid(arch, size=(grid.Nx, grid.Ny, grid.Nz÷2), extent=(grid.Lx, grid.Ly, grid.Lz))
    coarse_u = Field{Face, Center, Center}(coarse_grid)

    outputs = (; u = coarse_u)

    Arch = typeof(arch)
    filepath = "test_coarse_u_$Arch.nc"
    isfile(filepath) && rm(filepath)

    # NetCDFWriter should automatically use coarse_grid for dimensions
    # since coarse_u is a field on coarse_grid
    output_writer = NetCDFWriter(model, outputs;
                                 filename = filepath,
                                 schedule = IterationInterval(1),
                                 overwrite_existing = true)

    # Run simulation to write output
    simulation = Simulation(model, Δt=0.1, stop_iteration=2)
    simulation.output_writers[:coarse] = output_writer
    run!(simulation)

    # Verify that dimensions match the coarse grid, not the model grid
    ds = NCDataset(filepath)
    @test length(ds["z_aac"]) == grid.Nz÷2  # Should be 4, not 8
    @test length(ds["z_aaf"]) == grid.Nz÷2 + 1  # Should be 5, not 9

    # Verify that the z coordinates match the coarse grid
    expected_z_centers = znodes(coarse_grid, Center())
    expected_z_faces = znodes(coarse_grid, Face())
    @test ds["z_aac"][:] ≈ Array(expected_z_centers)
    @test ds["z_aaf"][:] ≈ Array(expected_z_faces)

    close(ds)
    rm(filepath)

    return nothing
end

function test_singleton_dimension_behavior(arch)
    Nx, Nz = 8, 8
    Hx, Hz = 2, 3
    Lx, H = 2, 1

    # Create grid with Flat in y dimension
    grid = RectilinearGrid(arch,
                           topology = (Periodic, Flat, Bounded),
                           size = (Nx, Nz),
                           halo = (Hx, Hz),
                           x = (-Lx, Lx),
                           z = (-H, 0))

    model = NonhydrostaticModel(grid;
                                  closure = ScalarDiffusivity(ν=4e-2, κ=4e-2),
                                  buoyancy = SeawaterBuoyancy(),
                                  tracers = (:T, :S))
    u_xavg = Field(Average(model.velocities.u, dims=(1)))

    Nt = 5
    simulation = Simulation(model, Δt=0.1, stop_iteration=Nt)

    Arch = typeof(arch)
    filepath_full = "test_singleton_full_u_$Arch.nc"
    simulation.output_writers[:full] = NetCDFWriter(model, (; u = model.velocities.u, u_xavg),
                                                    filename = filepath_full,
                                                    schedule = IterationInterval(1),
                                                    array_type = Array{Float64},
                                                    with_halos = false,
                                                    include_grid_metrics = false,
                                                    overwrite_existing = true)

    filepath_slice = "test_singleton_slice_u_$Arch.nc"
    simulation.output_writers[:slice] = NetCDFWriter(model, (; u = model.velocities.u, u_xavg),
                                                     filename = filepath_slice,
                                                     indices = (:, :, 1),
                                                     schedule = IterationInterval(1),
                                                     array_type = Array{Float64},
                                                     with_halos = false,
                                                     include_grid_metrics = false,
                                                     overwrite_existing = true)

    run!(simulation)

    # Verify full u output
    ds_full = NCDataset(filepath_full)
    @test haskey(ds_full, "u")
    @test haskey(ds_full, "u_xavg")
    @test dimsize(ds_full[:u]) == (x_faa=Nx, z_aac=Nz, time=Nt + 1)
    @test dimsize(ds_full[:u_xavg]) == (z_aac=Nz, time=Nt + 1)
    @test !haskey(ds_full, "y_aca")  # y dimension should not exist (Flat)
    @test !haskey(ds_full, "y_afa")  # y dimension should not exist (Flat)
    close(ds_full)
    rm(filepath_full)

    # Verify sliced u output (:, :, 1)
    ds_slice = NCDataset(filepath_slice)
    @test haskey(ds_slice, "u")
    @test haskey(ds_slice, "u_xavg")
    @test dimsize(ds_slice[:u]) == (x_faa=Nx, z_aac=1, time=Nt + 1)
    @test dimsize(ds_slice[:u_xavg]) == (z_aac=1, time=Nt + 1)
    @test !haskey(ds_slice, "y_aca")  # y dimension should not exist (Flat)
    @test !haskey(ds_slice, "y_afa")  # y dimension should not exist (Flat)
    close(ds_slice)
    rm(filepath_slice)

    return nothing
end

function test_netcdf_reduced_field_time_series(arch)
    # Regression test for reading reduced (Nothing-location) fields back as a
    # FieldTimeSeries from NetCDF. The on-disk variable is squeezed to fewer than 3
    # spatial dimensions, and the reader must re-inflate those dimensions (via
    # `inflate_nothing_dimensions`) before building the OffsetArray — otherwise
    # `offset_data` is handed a lower-dimensional array than the grid expects.
    Nx, Ny, Nz = 4, 5, 6
    grid = RectilinearGrid(arch, size=(Nx, Ny, Nz), extent=(1, 1, 1))
    model = NonhydrostaticModel(grid; tracers=:c)
    set!(model, c=(x, y, z) -> x + 2y + 3z)

    cx  = Field(Average(model.tracers.c, dims=1))        # (Nothing, Center, Center) → 2D on disk
    cxy = Field(Average(model.tracers.c, dims=(1, 2)))   # (Nothing, Nothing, Center) → 1D on disk
    cyz = Field(Average(model.tracers.c, dims=(2, 3)))   # (Center, Nothing, Nothing) → 1D on disk

    Nt = 3
    sim = Simulation(model; Δt=1e-3, stop_iteration=Nt)

    Arch = typeof(arch)
    fp = "test_reduced_fts_$Arch.nc"
    isfile(fp) && rm(fp)
    sim.output_writers[:nc] = NetCDFWriter(model, (; cx, cxy, cyz);
                                           filename=fp, schedule=IterationInterval(1),
                                           array_type=Array{Float64}, overwrite_existing=true)

    # Snapshot the reduced fields (and the save times) in memory on the same schedule as
    # the writer so the round-trip can be compared value-for-value.
    save_times = Float64[]
    snaps = (cx=Array{Float64,3}[], cxy=Array{Float64,3}[], cyz=Array{Float64,3}[])
    save = s -> begin
        push!(save_times, s.model.clock.time)
        push!(snaps.cx,  Array(interior(compute!(cx))))
        push!(snaps.cxy, Array(interior(compute!(cxy))))
        push!(snaps.cyz, Array(interior(compute!(cyz))))
    end
    sim.callbacks[:save] = Callback(save, IterationInterval(1))

    run!(sim)

    for (name, ref_loc, ref) in (("cx",  (Nothing, Center, Center), snaps.cx),
                                 ("cxy", (Nothing, Nothing, Center), snaps.cxy),
                                 ("cyz", (Center, Nothing, Nothing), snaps.cyz))
        fts = FieldTimeSeries(fp, name; architecture=arch)
        @test fts isa FieldTimeSeries
        @test location(fts) == ref_loc
        @test fts.times ≈ save_times
        @test size(fts) == (size(ref[1])..., length(ref))
        for k in 1:length(ref)
            rec = Array(interior(fts[k]))
            @test size(rec) == size(ref[k])
            @test rec ≈ ref[k]
        end
    end

    rm(fp)
    return nothing
end

function test_netcdf_dimension_type(arch)
    grid = RectilinearGrid(arch, size=(4, 4, 4), extent=(1, 1, 1))
    model = NonhydrostaticModel(grid)

    Nt = 3
    simulation = Simulation(model, Δt=0.1, stop_iteration=Nt)

    Arch = typeof(arch)

    # Test with default dimension_type (Float64)
    filepath_float64 = "test_dimension_type_float64_$Arch.nc"
    simulation.output_writers[:float64] = NetCDFWriter(model, (; u = model.velocities.u),
                                                       filename = filepath_float64,
                                                       schedule = IterationInterval(1),
                                                       include_grid_metrics = false,
                                                       overwrite_existing = true)

    # Test with Float32 dimension_type
    filepath_float32 = "test_dimension_type_float32_$Arch.nc"
    simulation.output_writers[:float32] = NetCDFWriter(model, (; u = model.velocities.u),
                                                        filename = filepath_float32,
                                                        schedule = IterationInterval(1),
                                                        dimension_type = Float32,
                                                        include_grid_metrics = false,
                                                        overwrite_existing = true)

    run!(simulation)

    # Verify Float64 dimensions (default)
    ds_float64 = NCDataset(filepath_float64)
    @test eltype(ds_float64["time"]) == Float64
    @test eltype(ds_float64["x_faa"]) == Float64
    @test eltype(ds_float64["x_caa"]) == Float64
    @test eltype(ds_float64["y_aca"]) == Float64
    @test eltype(ds_float64["y_afa"]) == Float64
    @test eltype(ds_float64["z_aac"]) == Float64
    @test eltype(ds_float64["z_aaf"]) == Float64
    close(ds_float64)
    rm(filepath_float64)

    # Verify Float32 dimensions
    ds_float32 = NCDataset(filepath_float32)
    @test eltype(ds_float32["time"]) == Float32
    @test eltype(ds_float32["x_faa"]) == Float32
    @test eltype(ds_float32["x_caa"]) == Float32
    @test eltype(ds_float32["y_aca"]) == Float32
    @test eltype(ds_float32["y_afa"]) == Float32
    @test eltype(ds_float32["z_aac"]) == Float32
    @test eltype(ds_float32["z_aaf"]) == Float32
    close(ds_float32)
    rm(filepath_float32)

    return nothing
end

#####
##### OrthogonalSphericalShellGrid output tests
#####
#
# These exercise the CF §5.2 auxiliary-coordinate path for OSSG (TripolarGrid,
# RotatedLatitudeLongitudeGrid, ConformalCubedSpherePanelGrid). The key cross-
# validation test uses a non-rotated RotatedLatitudeLongitudeGrid, whose data
# must match the equivalent LatitudeLongitudeGrid to roundoff.
#

using Oceananigans.OrthogonalSphericalShellGrids: TripolarGrid, RotatedLatitudeLongitudeGrid,
                                                  ConformalCubedSpherePanelGrid
using Oceananigans.OutputWriters: reconstruct_grid, materialize_from_netcdf
using Oceananigans.ImmersedBoundaries: ImmersedBoundaryGrid, GridFittedBottom, bottom_height_interior,
                                       CenterImmersedCondition, InterfaceImmersedCondition

function test_netcdf_tripolar_grid_output(arch)
    grid = TripolarGrid(arch, size=(20, 16, 4), z=(-100, 0))
    fs = SplitExplicitFreeSurface(grid; substeps=10)
    model = HydrostaticFreeSurfaceModel(grid; free_surface=fs, tracers=(:T,))

    simulation = Simulation(model; Δt=1, stop_iteration=2)

    Arch = typeof(arch)
    filepath = "test_tripolar_netcdf_$Arch.nc"
    isfile(filepath) && rm(filepath)

    simulation.output_writers[:nc] = NetCDFWriter(model, (; T=model.tracers.T, u=model.velocities.u);
                                                  filename = filepath,
                                                  schedule = IterationInterval(1),
                                                  overwrite_existing = true,
                                                  include_grid_metrics = true)
    run!(simulation)

    ds = NCDataset(filepath)

    # Bare horizontal dimensions: i_*, j_* should exist as dims but have no coord variable.
    for dname in ("i_caa", "i_faa", "j_aca", "j_afa")
        @test dname ∈ keys(ds.dim)
        @test dname ∉ keys(ds)
    end

    # All eight 2D λ/φ auxiliary coordinate variables.
    for v in ("λ_cca", "λ_fca", "λ_cfa", "λ_ffa", "φ_cca", "φ_fca", "φ_cfa", "φ_ffa")
        @test v ∈ keys(ds)
        @test ndims(ds[v]) == 2
    end

    # Center-cell λ/φ have the (Nx, Ny) shape.
    @test size(ds["λ_cca"]) == (size(grid, 1), size(grid, 2))
    @test size(ds["φ_cca"]) == (size(grid, 1), size(grid, 2))

    # `coordinates` attribute on data fields points at the right (λ, φ, z) trio.
    @test ds["T"].attrib["coordinates"] == "λ_cca φ_cca z_aac"
    @test ds["u"].attrib["coordinates"] == "λ_fca φ_fca z_aac"

    # Grid metrics: 4 × Δx, 4 × Δy, 4 × Az, plus Δz.
    for v in ("Δx_cca", "Δx_fca", "Δx_cfa", "Δx_ffa",
              "Δy_cca", "Δy_fca", "Δy_cfa", "Δy_ffa",
              "Az_cca", "Az_fca", "Az_cfa", "Az_ffa",
              "Δz_aac", "Δz_aaf")
        @test v ∈ keys(ds)
    end

    # Field shape: (Nx, Ny, Nz, Nt).
    @test size(ds["T"]) == (size(grid, 1), size(grid, 2), size(grid, 3), length(ds["time"]))

    close(ds)
    rm(filepath)
    return nothing
end

function test_netcdf_rotated_llg_matches_llg(arch)
    # The linchpin: with `north_pole = (0, 90)` the rotated grid is numerically
    # equivalent to a LatitudeLongitudeGrid with the same lat/lon/z. The two
    # NetCDF representations differ (1D coord vars vs. 2D aux coords) but the
    # data and coordinate values should agree to roundoff.

    N = (16, 12, 3)
    common = (size=N, longitude=(0, 360), latitude=(-60, 60), z=(-100, 0))

    llg  = LatitudeLongitudeGrid(arch; common...)
    rllg = RotatedLatitudeLongitudeGrid(arch; common..., north_pole=(0, 90))

    fs_llg  = SplitExplicitFreeSurface(llg;  substeps=5)
    fs_rllg = SplitExplicitFreeSurface(rllg; substeps=5)

    m_llg  = HydrostaticFreeSurfaceModel(llg;  free_surface=fs_llg,  tracers=(:T,))
    m_rllg = HydrostaticFreeSurfaceModel(rllg; free_surface=fs_rllg, tracers=(:T,))

    # Identical analytic initial condition. Must be 360°-periodic in λ so that the LLG
    # convention (λ ∈ [0, 360]) and the rotated-grid convention (λ ∈ [-180, 180]) yield
    # the same value at every grid point.
    T0(λ, φ, z) = cos(deg2rad(λ)) * exp(-φ^2 / 400)
    set!(m_llg.tracers.T,  T0)
    set!(m_rllg.tracers.T, T0)

    s_llg  = Simulation(m_llg;  Δt=1, stop_iteration=2)
    s_rllg = Simulation(m_rllg; Δt=1, stop_iteration=2)

    Arch = typeof(arch)
    p_llg  = "test_llg_$(Arch).nc"
    p_rllg = "test_rllg_$(Arch).nc"
    isfile(p_llg)  && rm(p_llg)
    isfile(p_rllg) && rm(p_rllg)

    s_llg.output_writers[:nc]  = NetCDFWriter(m_llg,  (; T=m_llg.tracers.T);
                                              filename=p_llg,  schedule=IterationInterval(1),
                                              overwrite_existing=true, include_grid_metrics=false)
    s_rllg.output_writers[:nc] = NetCDFWriter(m_rllg, (; T=m_rllg.tracers.T);
                                              filename=p_rllg, schedule=IterationInterval(1),
                                              overwrite_existing=true, include_grid_metrics=false)

    run!(s_llg)
    run!(s_rllg)

    ds_llg  = NCDataset(p_llg)
    ds_rllg = NCDataset(p_rllg)

    # Coordinates: 2D RotatedLLG aux coords match 1D LLG coords (broadcasted).
    # RotatedLLG normalizes λ to (-180, 180], so compare modulo 360.
    λ_llg  = collect(ds_llg["λ_caa"])
    φ_llg  = collect(ds_llg["φ_aca"])
    λ_rllg = collect(ds_rllg["λ_cca"])  # 2D (Nx, Ny)
    φ_rllg = collect(ds_rllg["φ_cca"])  # 2D (Nx, Ny)

    mod360(x) = mod(x + 360, 360)  # normalize to [0, 360)
    λ_llg_norm  = mod360.(λ_llg)
    λ_rllg_norm = mod360.(λ_rllg)
    @test all(j -> all(isapprox.(λ_rllg_norm[:, j], λ_llg_norm; atol=1e-8)), 1:N[2])
    @test all(i -> all(isapprox.(φ_rllg[i, :], φ_llg; atol=1e-10)), 1:N[1])

    # Field values agree at every (i, j, k, t).
    T_llg  = collect(ds_llg["T"])
    T_rllg = collect(ds_rllg["T"])
    @test size(T_llg) == size(T_rllg)
    @test isapprox(T_llg, T_rllg; atol=1e-6)

    close(ds_llg)
    close(ds_rllg)
    rm(p_llg)
    rm(p_rllg)
    return nothing
end

function test_netcdf_tripolar_field_time_series(arch)
    # Use a generous initial condition that stays well within the dynamic range and
    # `stop_iteration=1` to avoid any chance of numerical instability propagating into
    # the round-trip comparison — we're testing the I/O round-trip, not the physics.
    Nx, Ny, Nz = 12, 10, 3
    grid = TripolarGrid(arch, size=(Nx, Ny, Nz), z=(-100, 0))
    fs = SplitExplicitFreeSurface(grid; substeps=10)
    model = HydrostaticFreeSurfaceModel(grid; free_surface=fs, tracers=(:T,))

    # Smooth, bounded initial condition.
    set!(model.tracers.T, (λ, φ, z) -> sin(deg2rad(λ)) * cos(deg2rad(φ)))

    sim = Simulation(model; Δt=1, stop_iteration=1)

    Arch = typeof(arch)
    fp = "test_tripolar_fts_$Arch.nc"
    isfile(fp) && rm(fp)
    # Metrics-based OSSG reconstruction (used by `FieldTimeSeries`) requires
    # `include_grid_metrics=true` on the writer (the default) so the file carries
    # the Δx/Δy/Az arrays.
    sim.output_writers[:nc] = NetCDFWriter(model, (; T=model.tracers.T);
                                            filename=fp, schedule=IterationInterval(1),
                                            overwrite_existing=true, include_grid_metrics=true)

    # Capture in-memory snapshots that align with the file's time series — the
    # callback fires on the same IterationInterval(1) as the writer, so both record
    # the state after iterations 1, 2, …, stop_iteration.
    T_snapshots = Array{eltype(model.tracers.T), 3}[]
    saved = (s) -> push!(T_snapshots, Array(interior(s.model.tracers.T)))
    sim.callbacks[:save] = Callback(saved, IterationInterval(1))

    run!(sim)

    # Reconstruct as FieldTimeSeries. OSSG variants are reconstructed by reading the
    # saved λ/φ/Δx/Δy/Az/z arrays back directly (bypassing the user-facing constructor)
    # plus the serialized `conformal_mapping`, which preserves the TripolarGrid type
    # alias so default boundary conditions still pick `Zipper` for the north fold.
    fts = FieldTimeSeries(fp, "T"; architecture=arch)
    @test fts.grid isa TripolarGrid
    @test size(fts.grid) == size(grid)
    @test (fts.grid.Hx, fts.grid.Hy, fts.grid.Hz) == (grid.Hx, grid.Hy, grid.Hz)
    @test topology(fts.grid) == topology(grid)
    @test Array(fts.grid.λᶜᶜᵃ[1:Nx, 1:Ny]) ≈ Array(grid.λᶜᶜᵃ[1:Nx, 1:Ny]) atol=1e-10
    @test Array(fts.grid.φᶜᶜᵃ[1:Nx, 1:Ny]) ≈ Array(grid.φᶜᶜᵃ[1:Nx, 1:Ny]) atol=1e-10

    # Time values match.
    @test length(fts.times) == length(T_snapshots)

    # Field values agree at each time index. We use a NaN-tolerant comparison so the
    # test isn't sensitive to any halo-region NaNs that might leak into the interior
    # under unusual conditions — what matters here is that the on-disk values match
    # what was in memory, regardless of whether some cells happen to be NaN in both.
    nan_or_close(a, b; atol) = (isnan(a) && isnan(b)) || isapprox(a, b; atol=atol)
    for k in 1:min(length(fts.times), length(T_snapshots))
        rec = Array(interior(fts[k]))
        ref = T_snapshots[k]
        @test size(rec) == size(ref)
        @test all(nan_or_close.(rec, ref; atol=1e-4))
    end

    rm(fp)
    return nothing
end

function test_netcdf_tripolar_variable_z_output(arch)
    # Variable vertical spacing on a TripolarGrid. `StaticVerticalDiscretization` here
    # holds `OffsetVector{Float64, Vector{Float64}}` for the z-faces rather than
    # `StepRangeLen{Float64, TwicePrecision{…}}` (the uniform case), so the grid type
    # signature is meaningfully different — we want to make sure metrics + reconstruction
    # both work on this path too.
    Nx, Ny, Nz = 16, 12, 3
    z_faces = collect(range(-100.0, 0.0; length=Nz + 1))
    grid = TripolarGrid(arch, size=(Nx, Ny, Nz), z=z_faces)
    fs = SplitExplicitFreeSurface(grid; substeps=10)
    model = HydrostaticFreeSurfaceModel(grid; free_surface=fs, tracers=(:T,))
    sim = Simulation(model; Δt=1, stop_iteration=1)

    Arch = typeof(arch)
    fp = "test_tripolar_variable_z_$Arch.nc"
    isfile(fp) && rm(fp)
    sim.output_writers[:nc] = NetCDFWriter(model, (; T=model.tracers.T);
                                           filename=fp, schedule=IterationInterval(1),
                                           overwrite_existing=true, include_grid_metrics=true)
    run!(sim)

    ds = NCDataset(fp)
    @test "λ_cca" ∈ keys(ds)
    @test "Δz_aac" ∈ keys(ds)
    # The 1D z faces in the file should match what we put in.
    @test collect(ds["z_aaf"]) ≈ z_faces
    close(ds)

    # Metrics-based reconstruction should round-trip the variable z faithfully.
    recon = reconstruct_grid(fp; architecture=arch)
    @test recon isa TripolarGrid
    @test Array(recon.z.cᵃᵃᶠ[1:Nz+1]) ≈ z_faces atol=1e-10

    rm(fp)
    return nothing
end

function test_netcdf_tripolar_mvd_output(arch)
    # `MutableVerticalDiscretization` (z-star) on a `TripolarGrid` — combines two
    # mechanisms in this PR (OSSG aux-coord plumbing + `r`-named reference vertical
    # for MVD) plus the `Zipper` BC dispatch on the tripolar fold. Smoke-test that
    # all of the dim names land where we expect.
    Nx, Ny, Nz = 16, 12, 3
    zcoord = Oceananigans.Grids.MutableVerticalDiscretization(collect(range(-100.0, 0.0; length=Nz + 1)))
    grid = TripolarGrid(arch, size=(Nx, Ny, Nz), z=zcoord)
    fs = SplitExplicitFreeSurface(grid; substeps=10)
    ZStar = Oceananigans.Models.HydrostaticFreeSurfaceModels.ZStarCoordinate
    model = HydrostaticFreeSurfaceModel(grid; free_surface=fs, tracers=(:T,),
                                        vertical_coordinate=ZStar())
    sim = Simulation(model; Δt=1, stop_iteration=1)

    Arch = typeof(arch)
    fp = "test_tripolar_mvd_$Arch.nc"
    isfile(fp) && rm(fp)
    sim.output_writers[:nc] = NetCDFWriter(model, (; T=model.tracers.T);
                                           filename=fp, schedule=IterationInterval(1),
                                           overwrite_existing=true, include_grid_metrics=true)
    run!(sim)

    ds = NCDataset(fp)
    # Reference vertical coord is `r_*` (not `z_*`) on MVD grids, regardless of
    # horizontal grid type.
    @test "r_aac" ∈ keys(ds.dim) && "r_aaf" ∈ keys(ds.dim)
    @test "z_aac" ∉ keys(ds.dim) && "z_aaf" ∉ keys(ds.dim)
    @test "Δr_aac" ∈ keys(ds) && "Δr_aaf" ∈ keys(ds)
    # OSSG horizontal aux coords still present.
    for v in ("λ_cca", "φ_cca", "Δx_cca", "Δy_cca", "Az_cca")
        @test v ∈ keys(ds)
    end
    # Field dim signature picks up `r_aac` as the vertical dimension.
    @test dimnames(ds["T"]) == ("i_caa", "j_aca", "r_aac", "time")
    @test ds["T"].attrib["coordinates"] == "λ_cca φ_cca r_aac"
    close(ds)
    rm(fp)
    return nothing
end

function test_netcdf_tripolar_immersed_output(arch)
    # ImmersedBoundaryGrid wrapping a TripolarGrid: exercises the OSSG path AND the
    # immersed-boundary reconstruction-data writer (which `defVar`s `bottom_height`
    # in a subgroup whose dimension scope is local — `create_field_coord_variables!`
    # for OSSG fields has to defDim the bare `i_*`/`j_*` dims in the subgroup).
    Nx, Ny, Nz = 20, 16, 3
    ug = TripolarGrid(arch, size=(Nx, Ny, Nz), z=(-100, 0))
    grid = ImmersedBoundaryGrid(ug, GridFittedBottom((λ, φ) -> -50))

    fs = SplitExplicitFreeSurface(grid; substeps=10)
    model = HydrostaticFreeSurfaceModel(grid; free_surface=fs, tracers=(:T,))
    sim = Simulation(model; Δt=1, stop_iteration=1)

    Arch = typeof(arch)
    fp = "test_tripolar_ibg_$Arch.nc"
    isfile(fp) && rm(fp)
    sim.output_writers[:nc] = NetCDFWriter(model, (; T=model.tracers.T);
                                           filename=fp, schedule=IterationInterval(1),
                                           overwrite_existing=true, include_grid_metrics=true)
    run!(sim)

    ds = NCDataset(fp)
    @test "T" ∈ keys(ds)
    @test "λ_cca" ∈ keys(ds)
    @test "immersed_grid_reconstruction_args" ∈ keys(ds.group)
    @test "bottom_height" ∈ keys(ds.group["immersed_grid_reconstruction_args"])

    # Reconstruct via FieldTimeSeries: underlying TripolarGrid is rebuilt + wrapped
    # in an `ImmersedBoundaryGrid`, with the saved `bottom_height` populated.
    fts = FieldTimeSeries(fp, "T"; architecture=arch)
    @test fts.grid isa ImmersedBoundaryGrid
    @test fts.grid.underlying_grid isa TripolarGrid
    @test size(fts.grid) == size(grid)

    close(ds)
    rm(fp)
    return nothing
end

function test_netcdf_cubed_sphere_panel_output(arch)
    # `ConformalCubedSpherePanelGrid` is another `OrthogonalSphericalShellGrid` alias.
    # It carries a different `conformal_mapping` (CubedSphereConformalMapping) and a
    # very deep type signature; verify that basic write + the OSSG dim/aux-coord
    # plumbing handles it correctly.
    Nx, Ny, Nz = 8, 8, 3
    grid = ConformalCubedSpherePanelGrid(arch, size=(Nx, Ny, Nz), z=(-100, 0), halo=(2, 2, 2))
    model = HydrostaticFreeSurfaceModel(grid; tracers=(:T,))
    sim = Simulation(model; Δt=1, stop_iteration=1)

    Arch = typeof(arch)
    fp = "test_ccspg_$Arch.nc"
    isfile(fp) && rm(fp)
    sim.output_writers[:nc] = NetCDFWriter(model, (; T=model.tracers.T);
                                           filename=fp, schedule=IterationInterval(1),
                                           overwrite_existing=true, include_grid_metrics=true)
    run!(sim)

    ds = NCDataset(fp)
    # All eight 2D λ/φ auxiliary coords are written.
    for v in ("λ_cca", "λ_fca", "λ_cfa", "λ_ffa", "φ_cca", "φ_fca", "φ_cfa", "φ_ffa")
        @test v ∈ keys(ds)
    end
    @test ds["T"].attrib["coordinates"] == "λ_cca φ_cca z_aac"
    # All metrics, including vertical Δz, are written. (The vertical Δz Field is
    # constructed as a plain Field rather than via `Field(zspacings(grid, …))` to
    # avoid a Julia 1.12 specialization bug on the deep CCSPG conformal-mapping type;
    # see `ossg_vertical_spacing_field`.)
    for v in ("Δx_cca", "Δy_cca", "Az_cca", "Δz_aac", "Δz_aaf")
        @test v ∈ keys(ds)
    end

    # The reconstructed grid is a generic `OrthogonalSphericalShellGrid` — the
    # `CubedSphereConformalMapping` (with its ξ/η ranges + rotation) isn't yet
    # round-trippable through the conformal_mapping serialization, but the metric
    # arrays themselves are faithful.
    recon = reconstruct_grid(fp; architecture=arch)
    @test recon isa Oceananigans.OrthogonalSphericalShellGrids.OrthogonalSphericalShellGrid
    @test size(recon) == size(grid)
    @test Array(recon.λᶜᶜᵃ[1:Nx, 1:Ny]) ≈ Array(grid.λᶜᶜᵃ[1:Nx, 1:Ny]) atol=1e-10

    close(ds)
    rm(fp)
    return nothing
end

function test_netcdf_cubed_sphere_panel_immersed_output(arch)
    # ImmersedBoundaryGrid wrapping a ConformalCubedSpherePanelGrid. Exercises the
    # same OSSG-fields-in-subgroup path as `test_netcdf_tripolar_immersed_output`,
    # but with the deeper CCSPG type signature.
    Nx, Ny, Nz = 8, 8, 3
    ug = ConformalCubedSpherePanelGrid(arch, size=(Nx, Ny, Nz), z=(-100, 0), halo=(2, 2, 2))
    grid = ImmersedBoundaryGrid(ug, GridFittedBottom((λ, φ) -> -50))
    model = HydrostaticFreeSurfaceModel(grid; tracers=(:T,))
    sim = Simulation(model; Δt=1, stop_iteration=1)

    Arch = typeof(arch)
    fp = "test_ccspg_ibg_$Arch.nc"
    isfile(fp) && rm(fp)
    sim.output_writers[:nc] = NetCDFWriter(model, (; T=model.tracers.T);
                                           filename=fp, schedule=IterationInterval(1),
                                           overwrite_existing=true, include_grid_metrics=true)
    run!(sim)

    ds = NCDataset(fp)
    @test "T" ∈ keys(ds)
    @test "λ_cca" ∈ keys(ds)
    @test "immersed_grid_reconstruction_args" ∈ keys(ds.group)
    @test "bottom_height" ∈ keys(ds.group["immersed_grid_reconstruction_args"])
    close(ds)
    rm(fp)
    return nothing
end

function test_netcdf_rectilinear_mvd_output(arch)
    # Smoke test: RectilinearGrid + MutableVerticalDiscretization should write
    # `r_aac`/`r_aaf` (not `z_*`) and use them as the vertical field dimension.
    zcoord = Oceananigans.Grids.MutableVerticalDiscretization(collect(range(-100.0, 0.0; length=6)))
    grid = RectilinearGrid(arch; size=(8, 6, 5), x=(0, 1), y=(0, 1), z=zcoord,
                           topology=(Periodic, Periodic, Bounded))
    ZStar = Oceananigans.Models.HydrostaticFreeSurfaceModels.ZStarCoordinate
    model = HydrostaticFreeSurfaceModel(grid;
                                        free_surface=SplitExplicitFreeSurface(grid; substeps=10),
                                        tracers=(:T,),
                                        vertical_coordinate=ZStar())
    sim = Simulation(model; Δt=0.5, stop_iteration=2)

    Arch = typeof(arch)
    fp = "test_rect_mvd_$Arch.nc"
    isfile(fp) && rm(fp)
    sim.output_writers[:nc] = NetCDFWriter(model, (; T=model.tracers.T);
                                            filename=fp, schedule=IterationInterval(1),
                                            overwrite_existing=true, include_grid_metrics=true)
    run!(sim)

    ds = NCDataset(fp)
    # Vertical dimension is `r_*` (not `z_*`), since the discretization is mutable.
    @test "r_aac" ∈ keys(ds.dim) && "r_aaf" ∈ keys(ds.dim)
    @test "z_aac" ∉ keys(ds.dim) && "z_aaf" ∉ keys(ds.dim)
    @test "r_aac" ∈ keys(ds) && "r_aaf" ∈ keys(ds)

    # Vertical metric variables also use the `Δr_*` name.
    @test "Δr_aac" ∈ keys(ds) && "Δr_aaf" ∈ keys(ds)
    @test "Δz_aac" ∉ keys(ds) && "Δz_aaf" ∉ keys(ds)

    # Field dim signature ends in `r_aac` for a Center field.
    @test dimnames(ds["T"]) == ("x_caa", "y_aca", "r_aac", "time")

    close(ds)
    rm(fp)
    return nothing
end

# Minimal stand-in for a non-static vertical coordinate defined outside Oceananigans
# (e.g. in a downstream package); carries the reference-coordinate fields the NetCDF
# writer reads, shared with `MutableVerticalDiscretization`.
struct ReferenceTestVerticalCoordinate{C, D, E, F} <: Oceananigans.Grids.AbstractVerticalCoordinate
    cᵃᵃᶠ :: C
    cᵃᵃᶜ :: D
    Δᵃᵃᶠ :: E
    Δᵃᵃᶜ :: F
end

function test_netcdf_abstract_vertical_coordinate_name(arch)
    # The NetCDF writer labels the vertical axis from the coordinate's *type*: a
    # `StaticVerticalDiscretization` writes the physical `z`, while every other
    # `AbstractVerticalCoordinate` writes the reference `r`.
    custom = ReferenceTestVerticalCoordinate(0:3, 0.5:2.5, 1, 1)
    @test vertical_coordinate_name(custom) == "r"        # generic fallback (under test)

    mvd = Oceananigans.Grids.MutableVerticalDiscretization(collect(0.0:3.0))
    @test vertical_coordinate_name(mvd) == "r"           # specific MutableVerticalDiscretization method still wins

    static_grid = RectilinearGrid(arch; size=(1, 1, 2), x=(0, 1), y=(0, 1), z=(0, 1),
                                  topology=(Periodic, Periodic, Bounded))
    @test vertical_coordinate_name(static_grid.z) == "z" # static coordinate still writes z

    return nothing
end

function test_netcdf_tripolar_grid_reconstruction(arch)
    Nx, Ny, Nz = 20, 16, 3
    grid = TripolarGrid(arch, size=(Nx, Ny, Nz), z=(-100, 0))
    fs = SplitExplicitFreeSurface(grid; substeps=10)
    model = HydrostaticFreeSurfaceModel(grid; free_surface=fs, tracers=(:T,))
    simulation = Simulation(model; Δt=1, stop_iteration=1)

    Arch = typeof(arch)
    filepath = "test_tripolar_reconstruct_$Arch.nc"
    isfile(filepath) && rm(filepath)

    # Metrics-based OSSG reconstruction requires `include_grid_metrics=true` on the
    # writer (the default) so that the Δx/Δy/Az arrays are present in the file.
    simulation.output_writers[:nc] = NetCDFWriter(model, (; T=model.tracers.T);
                                                  filename=filepath,
                                                  schedule=IterationInterval(1),
                                                  overwrite_existing=true,
                                                  include_grid_metrics=true)
    run!(simulation)

    reconstructed = reconstruct_grid(filepath; architecture=arch)

    # OSSG reconstruction is metrics-based: the grid is rebuilt from the saved
    # λ/φ/Δx/Δy/Az/z arrays directly (bypassing the user-facing constructor) plus the
    # serialized conformal mapping, which preserves the `TripolarGrid` type alias.
    @test reconstructed isa TripolarGrid
    @test size(reconstructed) == size(grid)
    @test (reconstructed.Hx, reconstructed.Hy, reconstructed.Hz) == (grid.Hx, grid.Hy, grid.Hz)
    @test topology(reconstructed) == topology(grid)

    # Coordinate arrays agree (interior portion).
    @test Array(reconstructed.λᶜᶜᵃ[1:Nx, 1:Ny]) ≈ Array(grid.λᶜᶜᵃ[1:Nx, 1:Ny]) atol=1e-10
    @test Array(reconstructed.φᶜᶜᵃ[1:Nx, 1:Ny]) ≈ Array(grid.φᶜᶜᵃ[1:Nx, 1:Ny]) atol=1e-10

    rm(filepath)
    return nothing
end

function test_materialize_from_netcdf_strings()
    @test materialize_from_netcdf("(Oceananigans.Grids.Periodic, Oceananigans.Grids.Bounded)") === (Periodic, Bounded)
    @test materialize_from_netcdf("Oceananigans.Architectures.CPU()") isa CPU
    @test materialize_from_netcdf("Float64") === Float64

    @test materialize_from_netcdf("CenterImmersedCondition()") isa CenterImmersedCondition
    @test materialize_from_netcdf("InterfaceImmersedCondition()") isa InterfaceImmersedCondition
    return nothing
end

@testset "NetCDF output writer" begin
    @info "Testing NetCDF output writer..."

    for arch in archs
        A = typeof(arch)

        # Pre-build RectilinearGrids with different sizes
        N = 4
        rectilinear_grid1 = RectilinearGrid(arch, size=(N, N, N), extent=(1, 1, 1))
        rectilinear_grid2 = RectilinearGrid(arch, size=(6, 8, 5), extent=(2, 3, 1.5))

        # Pre-build LatitudeLongitudeGrids with different sizes
        latlon_grid1 = LatitudeLongitudeGrid(arch, size=(N, N, N), longitude=(-180, 180), latitude=(-90, 90), z=(-1, 0))
        latlon_grid2 = LatitudeLongitudeGrid(arch, size=(8, 6, 4), longitude=(-120, 60), latitude=(-60, 60), z=(-2, 0))

        @testset "DateTime and TimeDate output [$A]" begin
            @info "  Testing DateTime and TimeDate output [$A]..."
            test_datetime_netcdf_output(arch)
            test_timedate_netcdf_output(arch)
            test_netcdf_dimension_type(arch)
        end

        @testset "Grid metrics [$A]" begin
            @info "  Testing grid metrics [$A]..."

            @testset "Rectilinear grid metrics [$A]" begin
                test_netcdf_grid_metrics_rectilinear(arch, Float64)
                test_netcdf_grid_metrics_rectilinear(arch, Float32)
            end

            @testset "LatitudeLongitude grid metrics [$A]" begin
                test_netcdf_grid_metrics_latlon(arch, Float64)
                test_netcdf_grid_metrics_latlon(arch, Float32)
            end
        end

        @testset "Immersed boundary grids [$A]" begin
            @info "  Testing immersed boundary grids [$A]..."

            @testset "Rectilinear grid fitted bottom [$A]" begin
                test_netcdf_rectilinear_grid_fitted_bottom(arch, GridFittedBottom)
                test_netcdf_rectilinear_grid_fitted_bottom(arch, PartialCellBottom)
            end

            @testset "LatitudeLongitude grid fitted bottom [$A]" begin
                test_netcdf_latlon_grid_fitted_bottom(arch, GridFittedBottom)
                test_netcdf_latlon_grid_fitted_bottom(arch, PartialCellBottom)
            end
        end

        @testset "Flat dimensions [$A]" begin
            @info "  Testing flat dimensions [$A]..."
            test_netcdf_rectilinear_flat_xy(arch)
            test_netcdf_rectilinear_flat_xz(arch, immersed=false)
            test_netcdf_rectilinear_flat_xz(arch, immersed=true)
            test_netcdf_rectilinear_flat_yz(arch, immersed=false)
            test_netcdf_rectilinear_flat_yz(arch, immersed=true)
            test_netcdf_rectilinear_column(arch)
        end

        @testset "Thermal bubble output [$A]" begin
            @info "  Testing thermal bubble output [$A]..."
            test_thermal_bubble_netcdf_output(arch, Float64)
            test_thermal_bubble_netcdf_output(arch, Float32)
            test_thermal_bubble_netcdf_output(arch, Float64, with_halos=true)
            test_thermal_bubble_netcdf_output(arch, Float32, with_halos=true)
        end

        @testset "File splitting [$A]" begin
            @info "  Testing file splitting [$A]..."
            test_netcdf_size_file_splitting(arch)
            test_netcdf_time_file_splitting(arch)
        end

        @testset "Function and alignment output [$A]" begin
            @info "  Testing function and alignment output [$A]..."
            test_netcdf_function_output(arch)
            test_netcdf_output_alignment(arch)
        end

        @testset "Averaging [$A]" begin
            @info "  Testing averaging [$A]..."
            test_netcdf_spatial_average(arch)
            test_netcdf_time_averaging(arch)
        end

        @testset "Lagrangian particles [$A]" begin
            @info "  Testing Lagrangian particles [$A]..."
            test_netcdf_output_just_particles(arch)
            test_netcdf_output_particles_and_fields(arch)
        end

        @testset "Vertically stretched grid [$A]" begin
            @info "  Testing vertically stretched grid [$A]..."
            test_netcdf_vertically_stretched_grid_output(arch)
        end

        @testset "Overriding attributes [$A]" begin
            @info "  Testing overriding attributes [$A]..."
            test_netcdf_overriding_attributes(arch)
        end

        @testset "Free surface output [$A]" begin
            @info "  Testing free surface output [$A]..."
            for immersed in (false, true), vertically_stretched in (false, true)
                test_netcdf_hydrostatic_free_surface_only_output(arch; immersed, vertically_stretched)
                test_netcdf_hydrostatic_free_surface_mixed_output(arch; immersed, vertically_stretched)
                test_netcdf_nonhydrostatic_free_surface_only_output(arch; immersed, vertically_stretched)
                test_netcdf_nonhydrostatic_free_surface_mixed_output(arch; immersed, vertically_stretched)
            end
        end

        @testset "Buoyancy force [$A]" begin
            @info "  Testing buoyancy force [$A]..."
            test_netcdf_buoyancy_force(arch)
        end

        @testset "Different grid output [$A]" begin
            @info "  Testing different grid output [$A]..."
            test_netcdf_writer_different_grid(arch)
        end

        @testset "Singleton dimension behavior [$A]" begin
            @info "  Testing singleton dimension behavior [$A]..."
            test_singleton_dimension_behavior(arch)
        end

        @testset "Reduced FieldTimeSeries round-trip [$A]" begin
            @info "  Testing reduced FieldTimeSeries round-trip [$A]..."
            test_netcdf_reduced_field_time_series(arch)
        end

        @testset "Field defvar and dimension validation [$A]" begin
            @info "  Testing field defvar and dimension validation [$A]..."
            for grids in ((rectilinear_grid1, rectilinear_grid2),
                          (latlon_grid1, latlon_grid2))
                grid1, grid2 = grids
                test_netcdf_single_field_defvar(grid1, immersed=false)
                test_netcdf_single_field_defvar(grid1, immersed=true)
                test_netcdf_field_dimension_validation(grid1, immersed=false)
                test_netcdf_field_dimension_validation(grid1, immersed=true)
                test_netcdf_multiple_grids_defvar(grid1, grid2, immersed=false)
                test_netcdf_multiple_grids_defvar(grid1, grid2, immersed=true)
            end
        end

        @testset "OrthogonalSphericalShellGrid output [$A]" begin
            @info "  Testing OrthogonalSphericalShellGrid output [$A]..."
            test_netcdf_tripolar_grid_output(arch)
            test_netcdf_rotated_llg_matches_llg(arch)
            test_netcdf_tripolar_grid_reconstruction(arch)
            test_netcdf_tripolar_field_time_series(arch)
            test_netcdf_tripolar_variable_z_output(arch)
            test_netcdf_tripolar_mvd_output(arch)
            test_netcdf_tripolar_immersed_output(arch)
            test_netcdf_cubed_sphere_panel_output(arch)
            test_netcdf_cubed_sphere_panel_immersed_output(arch)
        end

        @testset "Non-static vertical coordinate output [$A]" begin
            @info "  Testing non-static vertical coordinate output [$A]..."
            test_netcdf_rectilinear_mvd_output(arch)
            test_netcdf_abstract_vertical_coordinate_name(arch)
        end
    end

    @testset "Materialize grid-reconstruction strings from NetCDF" begin
        @info "  Testing materialize_from_netcdf with qualified and bare type names..."
        test_materialize_from_netcdf_strings()
    end
end
