diff --git a/ext/OceananigansNCDatasetsExt.jl b/ext/OceananigansNCDatasetsExt.jl
index cdbfec9e125..4d915d5219b 100644
--- a/ext/OceananigansNCDatasetsExt.jl
+++ b/ext/OceananigansNCDatasetsExt.jl
@@ -1,21 +1,22 @@
 module OceananigansNCDatasetsExt
 
 using NCDatasets
+import NCDatasets: defVar
 
 using Dates: AbstractTime, UTC, now
 using Printf: @sprintf
 using OrderedCollections: OrderedDict
 
 using Oceananigans: initialize!, prettytime, pretty_filesize, AbstractModel
-using Oceananigans.Architectures: CPU, GPU
-using Oceananigans.AbstractOperations: KernelFunctionOperation
+using Oceananigans.Architectures: CPU, GPU, on_architecture
+using Oceananigans.AbstractOperations: KernelFunctionOperation, AbstractOperation
 using Oceananigans.BuoyancyFormulations: BuoyancyForce, BuoyancyTracer, SeawaterBuoyancy, LinearEquationOfState
 using Oceananigans.Fields
-using Oceananigans.Fields: reduced_dimensions, reduced_location, location
+using Oceananigans.Fields: Reduction, reduced_dimensions, reduced_location, location
 using Oceananigans.Grids: Center, Face, Flat, Periodic, Bounded,
                           AbstractGrid, RectilinearGrid, LatitudeLongitudeGrid, StaticVerticalDiscretization,
                           topology, halo_size, xspacings, yspacings, zspacings, λspacings, φspacings,
-                          parent_index_range, ξnodes, ηnodes, rnodes, validate_index, peripheral_node,
+                          parent_index_range, nodes, ξnodes, ηnodes, rnodes, validate_index, peripheral_node,
                           constructor_arguments
 using Oceananigans.ImmersedBoundaries: ImmersedBoundaryGrid, GridFittedBottom, GFBIBG, GridFittedBoundary, PartialCellBottom, PCBIBG
 using Oceananigans.Models: ShallowWaterModel, LagrangianParticles
@@ -41,13 +42,71 @@ using Oceananigans.OutputWriters:
     show_array_type
 
 import Oceananigans: write_output!
-import Oceananigans.OutputWriters: NetCDFWriter, write_grid_reconstruction_data!, convert_for_netcdf, materialize_from_netcdf, reconstruct_grid
+import Oceananigans.OutputWriters:
+    NetCDFWriter,
+    write_grid_reconstruction_data!,
+    convert_for_netcdf,
+    materialize_from_netcdf,
+    reconstruct_grid,
+    trilocation_dim_name
 
 const c = Center()
 const f = Face()
 const BoussinesqSeawaterBuoyancy = SeawaterBuoyancy{FT, <:BoussinesqEquationOfState, T, S} where {FT, T, S}
 const BuoyancyBoussinesqEOSModel = BuoyancyForce{<:BoussinesqSeawaterBuoyancy, g} where {g}
 
+#####
+##### Extend defVar to be able to write fields to NetCDF directly
+#####
+
+defVar(ds, name, op::AbstractOperation; kwargs...) = defVar(ds, name, Field(op); kwargs...)
+defVar(ds, name, op::Reduction; kwargs...) = defVar(ds, name, Field(op); kwargs...)
+
+function defVar(ds, name, field::AbstractField;
+                time_dependent=false,
+                dimension_name_generator = trilocation_dim_name,
+                kwargs...)
+    field_cpu = on_architecture(CPU(), field) # Need to bring field to CPU in order to write it to NetCDF
+    field_data = Array{eltype(field)}(field_cpu)
+    dims = field_dimensions(field, dimension_name_generator)
+    all_dims = time_dependent ? (dims..., "time") : dims
+
+    # Validate that all dimensions exist and match the field
+    create_field_dimensions!(ds, field, all_dims, dimension_name_generator)
+    defVar(ds, name, field_data, all_dims; kwargs...)
+end
+
+#####
+##### Dimension validation
+#####
+
+"""
+    create_field_dimensions!(ds, field::AbstractField, all_dims, dimension_name_generator)
+
+Creates all dimensions for the given `field` in the NetCDF dataset `ds`. If the dimensions
+already exist, they are validated to match the expected dimensions for the given `field`.
+
+Arguments:
+- `ds`: NetCDF dataset
+- `field`: AbstractField being written  
+- `all_dims`: Tuple of dimension names to create/validate
+- `dimension_name_generator`: Function to generate dimension names
+"""
+function create_field_dimensions!(ds, field::AbstractField, all_dims, dimension_name_generator)
+    dimension_attributes = default_dimension_attributes(field.grid, dimension_name_generator)
+    spatial_dims = all_dims[1:end-(("time" in all_dims) ? 1 : 0)]
+
+    spatial_dims_dict = Dict(dim_name => dim_data for (dim_name, dim_data) in zip(spatial_dims, nodes(field)))
+    create_spatial_dimensions!(ds, spatial_dims_dict, dimension_attributes; array_type=Array{eltype(field)})
+
+    # Create time dimension if needed
+    if "time" in all_dims && "time" ∉ keys(ds.dim)
+        create_time_dimension!(ds)
+    end
+    
+    return nothing
+end
+
 #####
 ##### Utils
 #####
@@ -67,73 +126,42 @@ function collect_dim(ξ, ℓ, T, N, H, inds, with_halos)
     end
 end
 
-#####
-##### Dimension name generators
-#####
-
-function suffixed_dim_name_generator(var_name, grid::AbstractGrid{FT, TX}, LX, LY, LZ, dim::Val{:x}; connector="_", location_letters) where {FT, TX}
-    if TX == Flat || isnothing(LX)
-        return ""
-    else
-        return "$(var_name)" * connector * location_letters
+function create_time_dimension!(dataset; attrib=nothing)
+    if "time" ∉ keys(dataset.dim)
+        # Create an unlimited dimension "time"
+        # Time should always be Float64 to be extra safe from rounding errors.
+        # See: https://github.com/CliMA/Oceananigans.jl/issues/3056
+        defDim(dataset, "time", Inf)
+        defVar(dataset, "time", Float64, ("time",), attrib=attrib)
     end
 end
 
-function suffixed_dim_name_generator(var_name, grid::AbstractGrid{FT, TX, TY}, LX, LY, LZ, dim::Val{:y}; connector="_", location_letters) where {FT, TX, TY}
-    if TY == Flat || isnothing(LY)
-        return ""
-    else
-        return "$(var_name)" * connector * location_letters
-    end
-end
 
-function suffixed_dim_name_generator(var_name, grid::AbstractGrid{FT, TX, TY, TZ}, LX, LY, LZ, dim::Val{:z}; connector="_", location_letters) where {FT, TX, TY, TZ}
-    if TZ == Flat || isnothing(LZ)
-        return ""
-    else
-        return "$(var_name)" * connector * location_letters
+"""
+    create_spatial_dimensions!(dataset, dims, attributes_dict; array_type=Array{Float32}, kwargs...)
+
+Create spatial dimensions in the NetCDF dataset and define corresponding variables to store
+their coordinate values. Each dimension variable has itself as its sole dimension (e.g., the
+`x` variable has dimension `x`). The dimensions are created if they don't exist, and validated
+against provided arrays if they do exist. An error is thrown if the dimension already exists
+but is different from the provided array.
+"""
+function create_spatial_dimensions!(dataset, dims, attributes_dict; array_type=Array{Float32}, kwargs...)
+    for (i, (dim_name, dim_array)) in enumerate(dims)
+        if dim_name ∉ keys(dataset.dim)
+            # Create missing dimension
+            defVar(dataset, dim_name, array_type(dim_array), (dim_name,), attrib=attributes_dict[dim_name]; kwargs...)
+        else
+            # Validate existing dimension
+            if dataset[dim_name] != dim_array
+                throw(ArgumentError("Dimension '$dim_name' already exists in dataset but is different from expected.\n" *
+                                    "  Actual:   $(dataset[dim_name]) (length=$(length(dataset[dim_name])))\n" *
+                                    "  Expected: $(dim_array) (length=$(length(dim_array)))"))
+            end
+        end
     end
 end
 
-suffixed_dim_name_generator(var_name, ::StaticVerticalDiscretization, LX, LY, LZ, dim::Val{:z}; connector="_", location_letters) = var_name * connector * location_letters
-
-loc2letter(::Face, full=true) = "f"
-loc2letter(::Center, full=true) = "c"
-loc2letter(::Nothing, full=true) = full ? "a" : ""
-
-minimal_location_string(::RectilinearGrid, LX, LY, LZ, ::Val{:x}) = loc2letter(LX, false)
-minimal_location_string(::RectilinearGrid, LX, LY, LZ, ::Val{:y}) = loc2letter(LY, false)
-
-minimal_location_string(::LatitudeLongitudeGrid, LX, LY, LZ, ::Val{:x}) = loc2letter(LX, false) * loc2letter(LY, false)
-minimal_location_string(::LatitudeLongitudeGrid, LX, LY, LZ, ::Val{:y}) = loc2letter(LX, false) * loc2letter(LY, false)
-
-minimal_location_string(grid::AbstractGrid,             LX, LY, LZ, dim::Val{:z}) = minimal_location_string(grid.z, LX, LY, LZ, dim)
-minimal_location_string(::StaticVerticalDiscretization, LX, LY, LZ, dim::Val{:z}) = loc2letter(LZ, false)
-minimal_location_string(grid,                           LX, LY, LZ, dim)          = loc2letter(LX, false) * loc2letter(LY, false) * loc2letter(LZ, false)
-
-minimal_dim_name(var_name, grid, LX, LY, LZ, dim) =
-    suffixed_dim_name_generator(var_name, grid, LX, LY, LZ, dim, connector="_", location_letters=minimal_location_string(grid, LX, LY, LZ, dim))
-
-minimal_dim_name(var_name, grid::ImmersedBoundaryGrid, args...) = minimal_dim_name(var_name, grid.underlying_grid, args...)
-
-
-
-trilocation_location_string(::RectilinearGrid, LX, LY, LZ, ::Val{:x}) = loc2letter(LX) * "aa"
-trilocation_location_string(::RectilinearGrid, LX, LY, LZ, ::Val{:y}) = "a" * loc2letter(LY) * "a"
-
-trilocation_location_string(::LatitudeLongitudeGrid, LX, LY, LZ, ::Val{:x}) = loc2letter(LX) * loc2letter(LY) * "a"
-trilocation_location_string(::LatitudeLongitudeGrid, LX, LY, LZ, ::Val{:y}) = loc2letter(LX) * loc2letter(LY) * "a"
-
-trilocation_location_string(grid::AbstractGrid,             LX, LY, LZ, dim::Val{:z}) = trilocation_location_string(grid.z, LX, LY, LZ, dim)
-trilocation_location_string(::StaticVerticalDiscretization, LX, LY, LZ, dim::Val{:z}) = "aa" * loc2letter(LZ)
-trilocation_location_string(grid,                           LX, LY, LZ, dim)          = loc2letter(LX) * loc2letter(LY) * loc2letter(LZ)
-
-trilocation_dim_name(var_name, grid, LX, LY, LZ, dim) =
-    suffixed_dim_name_generator(var_name, grid, LX, LY, LZ, dim, connector="_", location_letters=trilocation_location_string(grid, LX, LY, LZ, dim))
-
-trilocation_dim_name(var_name, grid::ImmersedBoundaryGrid, args...) = trilocation_dim_name(var_name, grid.underlying_grid, args...)
-
-
 #####
 ##### Gathering of grid dimensions
 #####
@@ -417,8 +445,8 @@ function field_dimensions(field::AbstractField, grid::LatitudeLongitudeGrid, dim
     LΛ, LΦ, LZ = location(field)
 
     λ_dim_name = dim_name_generator("λ", grid, LΛ(), nothing, nothing, Val(:x))
-    φ_dim_name = dim_name_generator("φ",  grid, nothing, LΦ(), nothing, Val(:y))
-    z_dim_name = dim_name_generator("z",         grid, nothing, nothing, LZ(), Val(:z))
+    φ_dim_name = dim_name_generator("φ", grid, nothing, LΦ(), nothing, Val(:y))
+    z_dim_name = dim_name_generator("z", grid, nothing, nothing, LZ(), Val(:z))
 
     λ_dim_name = isempty(λ_dim_name) ? tuple() : tuple(λ_dim_name)
     φ_dim_name = isempty(φ_dim_name) ? tuple() : tuple(φ_dim_name)
@@ -1146,18 +1174,8 @@ function initialize_nc_file(model,
             Dict("long_name" => "Time", "units" => "seconds since 2000-01-01 00:00:00") :
             Dict("long_name" => "Time", "units" => "seconds")
 
-        # Create an unlimited dimension "time"
-        # Time should always be Float64 to be extra safe from rounding errors.
-        # See: https://github.com/CliMA/Oceananigans.jl/issues/3056
-        defDim(dataset, "time", Inf)
-        defVar(dataset, "time", Float64, ("time",), attrib=time_attrib)
-
-        # Create spatial dimensions as variables whose dimensions are themselves.
-        # Each should already have a default attribute.
-        for (dim_name, dim_array) in dims
-            defVar(dataset, dim_name, array_type(dim_array), (dim_name,),
-                   deflatelevel=deflatelevel, attrib=output_attributes[dim_name])
-        end
+        create_time_dimension!(dataset, attrib=time_attrib)
+        create_spatial_dimensions!(dataset, dims, output_attributes; deflatelevel=1, array_type)
 
         time_independent_vars = Dict()
 
diff --git a/src/OutputWriters/netcdf_writer.jl b/src/OutputWriters/netcdf_writer.jl
index fb323e8335a..e1a8a90bcba 100644
--- a/src/OutputWriters/netcdf_writer.jl
+++ b/src/OutputWriters/netcdf_writer.jl
@@ -4,6 +4,66 @@
 ##### NetCDFWriter functionality is implemented in ext/OceananigansNCDatasetsExt
 #####
 
+using Oceananigans.Grids: topology, Flat, StaticVerticalDiscretization
+using Oceananigans.ImmersedBoundaries: ImmersedBoundaryGrid
+
+#####
+##### Dimension name generators
+#####
+
+function suffixed_dim_name_generator(var_name, grid::AbstractGrid{FT, TX}, LX, LY, LZ, dim::Val{:x}; connector="_", location_letters) where {FT, TX}
+    if TX == Flat || isnothing(LX)
+        return ""
+    else
+        return "$(var_name)" * connector * location_letters
+    end
+end
+
+function suffixed_dim_name_generator(var_name, grid::AbstractGrid{FT, TX, TY}, LX, LY, LZ, dim::Val{:y}; connector="_", location_letters) where {FT, TX, TY}
+    if TY == Flat || isnothing(LY)
+        return ""
+    else
+        return "$(var_name)" * connector * location_letters
+    end
+end
+
+function suffixed_dim_name_generator(var_name, grid::AbstractGrid{FT, TX, TY, TZ}, LX, LY, LZ, dim::Val{:z}; connector="_", location_letters) where {FT, TX, TY, TZ}
+    if TZ == Flat || isnothing(LZ)
+        return ""
+    else
+        return "$(var_name)" * connector * location_letters
+    end
+end
+
+suffixed_dim_name_generator(var_name, ::StaticVerticalDiscretization, LX, LY, LZ, dim::Val{:z}; connector="_", location_letters) = var_name * connector * location_letters
+
+loc2letter(::Face, full=true) = "f"
+loc2letter(::Center, full=true) = "c"
+loc2letter(::Nothing, full=true) = full ? "a" : ""
+
+minimal_location_string(::RectilinearGrid, LX, LY, LZ, ::Val{:x}) = loc2letter(LX, false)
+minimal_location_string(::RectilinearGrid, LX, LY, LZ, ::Val{:y}) = loc2letter(LY, false)
+
+minimal_dim_name(var_name, grid, LX, LY, LZ, dim) =
+minimal_dim_name(var_name, grid::ImmersedBoundaryGrid, args...) = minimal_dim_name(var_name, grid.underlying_grid, args...)
+
+trilocation_location_string(::RectilinearGrid, LX, LY, LZ, ::Val{:x}) = loc2letter(LX) * "aa"
+trilocation_location_string(::RectilinearGrid, LX, LY, LZ, ::Val{:y}) = "a" * loc2letter(LY) * "a"
+
+trilocation_location_string(::LatitudeLongitudeGrid, LX, LY, LZ, ::Val{:x}) = loc2letter(LX) * loc2letter(LY) * "a"
+trilocation_location_string(::LatitudeLongitudeGrid, LX, LY, LZ, ::Val{:y}) = loc2letter(LX) * loc2letter(LY) * "a"
+
+trilocation_location_string(grid::AbstractGrid,             LX, LY, LZ, dim::Val{:z}) = trilocation_location_string(grid.z, LX, LY, LZ, dim)
+trilocation_location_string(::StaticVerticalDiscretization, LX, LY, LZ, dim::Val{:z}) = "aa" * loc2letter(LZ)
+trilocation_location_string(grid,                           LX, LY, LZ, dim)          = loc2letter(LX) * loc2letter(LY) * loc2letter(LZ)
+
+trilocation_dim_name(var_name, grid, LX, LY, LZ, dim) =
+    suffixed_dim_name_generator(var_name, grid, LX, LY, LZ, dim, connector="_", location_letters=trilocation_location_string(grid, LX, LY, LZ, dim))
+
+trilocation_dim_name(var_name, grid::ImmersedBoundaryGrid, args...) = trilocation_dim_name(var_name, grid.underlying_grid, args...)
+
+
+
 mutable struct NetCDFWriter{G, D, O, T, A, FS, DN} <: AbstractOutputWriter
     grid :: G
     filepath :: String
diff --git a/test/test_netcdf_writer.jl b/test/test_netcdf_writer.jl
index bdbe04be7a3..8dfe0431275 100644
--- a/test/test_netcdf_writer.jl
+++ b/test/test_netcdf_writer.jl
@@ -11,6 +11,8 @@ using SeawaterPolynomials.SecondOrderSeawaterPolynomials: RoquetEquationOfState
 
 using Oceananigans: Clock
 using Oceananigans.Models.HydrostaticFreeSurfaceModels: VectorInvariant
+using Oceananigans.OutputWriters: trilocation_dim_name
+using Oceananigans.Grids: ξname, ηname, rname, ξnodes, ηnodes
 
 function test_datetime_netcdf_output(arch)
     grid = RectilinearGrid(arch, size=(1, 1, 1), extent=(1, 1, 1))
@@ -33,7 +35,6 @@ function test_datetime_netcdf_output(arch)
                                                       filename = filepath,
                                                       schedule = IterationInterval(1),
                                                       include_grid_metrics = false)
-
     run!(simulation)
 
     ds = NCDataset(filepath)
@@ -2676,7 +2677,6 @@ function test_netcdf_free_surface_mixed_output(arch)
 end
 
 function test_netcdf_buoyancy_force(arch)
-
     Nx, Nz = 8, 8
     Hx, Hz = 2, 3
     Lx, H  = 2, 1
@@ -2697,7 +2697,6 @@ function test_netcdf_buoyancy_force(arch)
                       RoquetEquationOfState(:SimplestRealistic))
 
     for eos in Boussinesq_eos
-
         model = NonhydrostaticModel(; grid,
                                     closure = ScalarDiffusivity(ν=4e-2, κ=4e-2),
                                     buoyancy = SeawaterBuoyancy(equation_of_state=eos),
@@ -2726,10 +2725,176 @@ function test_netcdf_buoyancy_force(arch)
     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
+
 for arch in archs
     @testset "NetCDF output writer [$(typeof(arch))]" begin
         @info "  Testing NetCDF output writer [$(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))
+
         test_datetime_netcdf_output(arch)
         test_timedate_netcdf_output(arch)
 
@@ -2775,5 +2940,16 @@ for arch in archs
         test_netcdf_free_surface_mixed_output(arch)
 
         test_netcdf_buoyancy_force(arch)
+
+        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
 end