diff --git a/src/Grids/level.jl b/src/Grids/level.jl
index 6f84e4a6b2..ba187bed15 100644
--- a/src/Grids/level.jl
+++ b/src/Grids/level.jl
@@ -6,8 +6,7 @@ struct LevelGrid{
     level::L
 end
 
-quadrature_style(levelgrid::LevelGrid) =
-    quadrature_style(levelgrid.full_grid.horizontal_grid)
+quadrature_style(levelgrid::LevelGrid) = quadrature_style(levelgrid.full_grid)
 
 level(
     grid::AbstractExtrudedFiniteDifferenceGrid,
@@ -23,6 +22,7 @@ dss_weights(levelgrid::LevelGrid, _) = dss_weights(levelgrid.full_grid, nothing)
 
 local_geometry_type(::Type{LevelGrid{G, L}}) where {G, L} =
     local_geometry_type(G)
+
 local_geometry_data(levelgrid::LevelGrid{<:Any, Int}, ::Nothing) = level(
     local_geometry_data(levelgrid.full_grid, CellCenter()),
     levelgrid.level,
@@ -32,7 +32,7 @@ local_geometry_data(levelgrid::LevelGrid{<:Any, PlusHalf{Int}}, ::Nothing) =
         local_geometry_data(levelgrid.full_grid, CellFace()),
         levelgrid.level + half,
     )
-global_geometry(levlgrid::LevelGrid) = global_geometry(levlgrid.full_grid)
+global_geometry(levelgrid::LevelGrid) = global_geometry(levelgrid.full_grid)
 
 ## GPU compatibility
 Adapt.adapt_structure(to, grid::LevelGrid) =
diff --git a/src/InputOutput/readers.jl b/src/InputOutput/readers.jl
index 12376f545a..4aaec83427 100644
--- a/src/InputOutput/readers.jl
+++ b/src/InputOutput/readers.jl
@@ -522,6 +522,7 @@ function read_grid_new(reader, name)
         else
             level = attrs(group)["level_half"] + half
         end
+
         return Grids.LevelGrid(full_grid, level)
     else
         error("Unsupported grid type $type")
diff --git a/test/Fields/unit_field.jl b/test/Fields/unit_field.jl
index 3401bff879..ae47438f35 100644
--- a/test/Fields/unit_field.jl
+++ b/test/Fields/unit_field.jl
@@ -689,6 +689,35 @@ end
     end
 end
 
+@testset "Levels of nonlocal Fields and nonlocal Field broadcasts" begin
+    FT = Float64
+    gradh = Operators.Gradient()
+    # Todo: Make this work over all spaces; currently broken for everything else.
+    for space in (
+        TU.CenterExtrudedFiniteDifferenceSpace(FT),
+        TU.FaceExtrudedFiniteDifferenceSpace(FT),
+    )
+        TU.levelable(space) || continue
+        field = fill((; x = FT(1)), space)
+
+        op_on_level_of_field =
+            gradh.(
+                Fields.Field(
+                    Spaces.level(Fields.field_values(field.x), 1),
+                    Spaces.level(space, TU.fc_index(1, space)),
+                )
+            )
+
+        @test op_on_level_of_field ==
+              (Spaces.level(gradh.(field.x), TU.fc_index(1, space)))
+
+        @test_broken op_on_level_of_field == Base.materialize((Spaces.level(
+            lazy.(gradh.(field.x)),
+            TU.fc_index(1, space),
+        )),)
+    end
+end
+
 @testset "Columns of Fields and Field broadcasts" begin
     FT = Float64
     for space in TU.all_spaces(FT)
diff --git a/test/InputOutput/unit_hybrid3dcubedsphere_topography.jl b/test/InputOutput/unit_hybrid3dcubedsphere_topography.jl
index 0260d34cf1..f11d3d8384 100644
--- a/test/InputOutput/unit_hybrid3dcubedsphere_topography.jl
+++ b/test/InputOutput/unit_hybrid3dcubedsphere_topography.jl
@@ -13,6 +13,11 @@ using ClimaCore:
     InputOutput,
     Grids
 
+@isdefined(TU) || include(
+    joinpath(pkgdir(ClimaCore), "test", "TestUtilities", "TestUtilities.jl"),
+);
+import .TestUtilities as TU;
+
 using ClimaComms
 const comms_ctx = ClimaComms.context(ClimaComms.CPUSingleThreaded())
 pid, nprocs = ClimaComms.init(comms_ctx)
@@ -54,7 +59,7 @@ end
                 z_max / 8 .* (
                     cosd.(Fields.coordinate_field(h_space).lat) .+
                     cosd.(Fields.coordinate_field(h_space).long) .+ 1
-                )
+                ),
             )
 
         z_mesh = Meshes.IntervalMesh(z_domain, nelems = z_elem)
@@ -90,3 +95,50 @@ end
         end
     end
 end
+
+
+@testset "HDF5 restart test for a Named Tuple of Levels of a 3D hybrid cubed sphere for deep" begin
+    # This I/O is used for the computation of the topographic drag
+    FT = Float32
+
+    for space in (
+        TU.CenterExtrudedFiniteDifferenceSpace(FT, context = comms_ctx),
+        TU.FaceExtrudedFiniteDifferenceSpace(FT, context = comms_ctx),
+    )
+        TU.levelable(space) || continue
+        field = fill((; x = FT(1)), space)
+
+        level_of_field = Fields.Field(
+            Spaces.level(Fields.field_values(field.x), 1),
+            Spaces.level(space, TU.fc_index(1, space)),
+        )
+
+        fake_drag = fill(
+            (;
+                t11 = FT(0.0),
+                t12 = FT(0.0),
+                t21 = FT(0.0),
+                t22 = FT(0.0),
+                hmin = FT(0.0),
+                hmax = FT(0.0),
+            ),
+            axes(level_of_field),
+        )
+
+        # write field vector to hdf5 file
+        InputOutput.HDF5Writer(filename, comms_ctx) do writer
+            InputOutput.write!(writer, fake_drag, "fake_drag")
+        end
+
+        InputOutput.HDF5Reader(filename, comms_ctx) do reader
+            restart_fake_drag = InputOutput.read_field(reader, "fake_drag") # read fieldvector from hdf5 file
+
+            # The underlying space is of a different instance, so we cannot use == to check for equivalence.
+            # Instead, we make sure that the values and types are the same.
+            @test typeof(restart_fake_drag) == typeof(fake_drag)
+            @test maximum(
+                abs.(parent(fake_drag.t21) .- parent(restart_fake_drag.t21)),
+            ) == 0.0f0
+        end
+    end
+end