@testset "Δz_field" begin
    FT = Float64
    context = ClimaComms.SingletonCommsContext()
    x = FT(1)
    y = FT(2)
    z = FT(3)
    lat, long = FT(4), FT(5)
    x1 = FT(1)
    x2 = FT(2)
    x3 = FT(3)
    coords = [
        Geometry.ZPoint(z),
        Geometry.XZPoint(x, z),
        Geometry.XYZPoint(x, y, z),
        Geometry.LatLongZPoint(lat, long, z),
        Geometry.Cartesian3Point(x3),
        Geometry.Cartesian13Point(x1, x3),
        Geometry.Cartesian123Point(x1, x2, x3),
    ]
    all_components = [
        SMatrix{1,1}(FT[1]),
        SMatrix{2,2}(FT[1 2; 3 4]),
        SMatrix{3,3}(FT[1 2 10; 4 5 6; 7 8 9]),
        SMatrix{3,3}(FT[1 2 10; 4 5 6; 7 8 9]),
        SMatrix{2,2}(FT[1 2; 3 4]),
        SMatrix{3,3}(FT[1 2 10; 4 5 6; 7 8 9])
    ]

    expected_dzs = [1.0, 4.0, 9.0, 9.0, 1.0, 2.0, 9.0]

    for (components, coord, expected_dz) in
        zip(all_components, coords, expected_dzs)
        CoordType = typeof(coord)
        AIdx = Geometry.coordinate_axis(CoordType)
        at = Geometry.AxisTensor(
            (Geometry.LocalAxis{AIdx}(), Geometry.CovariantAxis{AIdx}()),
            components,
        )
        local_geometry = Geometry.LocalGeometry(coord, FT(1.0), FT(1.0), at)
        space = Spaces.PointSpace(context, local_geometry)
        dz_computed = Array(parent(Fields.Δz_field(space)))
        @test length(dz_computed) == 1
        @test dz_computed[1] == expected_dz
    end
end