diff --git a/NEWS.md b/NEWS.md
index 629987be79..4a863740f5 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -4,6 +4,24 @@ ClimaCore.jl Release Notes
 main
 -------
 
+### ![][badge-✨feature/enhancement] Various improvements to `Remapper` [2060](https://github.com/CliMA/ClimaCore.jl/pull/2060)
+
+The `ClimaCore.Remapping` module received two improvements. First, `Remapper` is
+now compatible with purely vertical `Space`s (performing a linear
+interpolation), making it compatible with column setups. Second, a new set of
+simplified interpolation functions are provided.
+
+Now, interpolating a `Field` `field` is as easy as
+```julia
+import ClimaCore.Remapping: interpolate
+output_array = interpolate(field)
+```
+The target coordinates are automatically determined, but can also be customized.
+Refer to the [documentation](https://clima.github.io/ClimaCore.jl/dev/remapping/)
+for more information.
+
+v0.14.20
+---------
  - We've added new convenience constructors for spaces PR [2082](https://github.com/CliMA/ClimaCore.jl/pull/2082). Here are links to the new constructors:
    - [ExtrudedCubedSphereSpace]()
    - [CubedSphereSpace]()
@@ -12,8 +30,6 @@ main
    - [SliceXZSpace]()
    - [RectangleXYSpace]()
 
-v0.14.20
-
  - We've added new convenience constructors for grids PR [1848](https://github.com/CliMA/ClimaCore.jl/pull/1848). Here are links to the new constructors:
    - [ExtrudedCubedSphereGrid](https://github.com/CliMA/ClimaCore.jl/blob/cbb193042fac3b4bef33251fbc0f232427bfe506/src/CommonGrids/CommonGrids.jl#L85-L144)
    - [CubedSphereGrid](https://github.com/CliMA/ClimaCore.jl/blob/cbb193042fac3b4bef33251fbc0f232427bfe506/src/CommonGrids/CommonGrids.jl#L200-L235)
diff --git a/docs/make.jl b/docs/make.jl
index a2496aed68..66bed38583 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -77,6 +77,7 @@ withenv("GKSwstype" => "nul") do
             "Installation and How-to Guides" => "installation_instructions.md",
             "Geometry" => "geometry.md",
             "Operators" => "operators.md",
+            "Remapping" => "remapping.md",
             "MatrixFields" => "matrix_fields.md",
             "API" => "api.md",
             "Developer docs" => ["Performance tips" => "performance_tips.md"],
diff --git a/docs/src/remapping.md b/docs/src/remapping.md
new file mode 100644
index 0000000000..16bcccc8cf
--- /dev/null
+++ b/docs/src/remapping.md
@@ -0,0 +1,159 @@
+# Remapping to regular grids
+
+`ClimaCore` horizontal domains are spectral elements. Points are not distributed
+uniformly within an element, and elements are also not necessarily organized in
+a simple way. For these reasons, remapping to regular grids becomes a
+fundamental operations when inspecting the simulation output. In this section,
+we describe the remappers currently available in `ClimaCore`.
+
+Broadly speaking, we can classify remappers in two categories: conservative, and
+non-conservative. Conservative remappers preserve areas (and masses) when going
+from the spectral grid to Cartesian ones. Conservative remappers are non-local
+operations (meaning that they require communication between different elements)
+and are more expensive, so they are typically reserved to operations where
+physical conservation is important (e.g., exchange between component models in a
+coupled simulation). On the other hand, non-conservative remappers are local to
+an element and faster to evaluate, which makes them suitable to operations like
+diagnostics and plotting, where having perfect physical conservation is not as
+important.
+
+## Non-conservative remapping
+
+Non-conservative remappers are fast and do not require communication, but they
+are not as accurate as conservative remappers, especially with large elements
+with sharp gradients. These remappers are optimally suited for diagnostics and
+plots.
+
+The main non-conservative remapper currently implemented utilizes a Lagrange
+interpolation with the barycentric formula in [Berrut2004], equation (3.2), for
+the horizontal interpolation. Vertical interpolation is linear except in the
+boundary elements where it is 0th order.
+
+### Quick start
+
+Assuming you have a `ClimaCore` `Field` with name `field`, the simplest way to
+interpolate onto a uniform grid is with
+```julia
+julia> import ClimaCore.Remapping
+julia> Remapping.interpolate(field)
+```
+
+This will return an `Array` (or a `CuArray`) with the `field` interpolated on
+some uniform grid that is automatically determined based on the `Space` of
+definition of `field`. To obtain such coordinates, you can call the
+[`Remapping.default_target_hcoords`](@ref) and
+[`Remapping.default_target_zcoords`](@ref) functions. These functions return an
+`Array` with the coordinates over which interpolation will occur. These arrays are
+of `Geometry.Point`s.
+
+[`Remapping.interpolate`](@ref) allocates new output arrays. As such, it is not
+suitable for performance-critical applications. [`Remapping.interpolate!`](@ref)
+performs interpolation in-place. When using the in-place version`, the
+`dest`ination has to be the same array type as the device in use (e.g.,
+`CuArray` for CUDA runs). For performance-critical applications, it is preferable to a
+[`Remapping.Remapper`](@ref) and use it directly (see next Section).
+
+#### Example
+
+Given `field`, a `Field` defined on a cubed sphere.
+
+By default, a target uniform grid is chosen (with resolution `hresolution` and
+`vresolution`), so remapping is
+```julia
+interpolated_array = interpolate(field, hcoords, zcoords)
+```
+Coordinates can be specified:
+
+```julia
+longpts = range(-180.0, 180.0, 21)
+latpts = range(-80.0, 80.0, 21)
+zpts = range(0.0, 1000.0, 21)
+
+hcoords = [Geometry.LatLongPoint(lat, long) for long in longpts, lat in latpts]
+zcoords = [Geometry.ZPoint(z) for z in zpts]
+
+interpolated_array = interpolate(field, hcoords, zcoords)
+```
+The output is defined on the Cartesian product of `hcoords` with `zcoords`.
+
+If the default target coordinates are being used, it is possible to broadcast
+`ClimaCore.Geometry.components` to extract them as a vector of tuples (and then
+broadcast `getindex` to extract the respective coordinates as vectors).
+
+### The `Remapper` object
+
+A [`Remapping.Remapper`](@ref) is an object that is tied to a specified `Space`
+and can interpolate scalar `Field`s defined on that space onto a predefined
+target grid. The grid does not have to be regular, but it has to be defined as a
+Cartesian product between some horizontal and vertical coordinates (meaning, for
+each horizontal point, there is a fixed column of vertical coordinates).
+
+Let us create our first remapper, assuming we have `space` defined on the surface of the sphere
+```julia
+import ClimaCore.Geometry: LatLongPoint, ZPoint 
+import ClimaCore.Remapping: Remapper
+
+hcoords = [Geometry.LatLongPoint(lat, long) for long in -180.:180., lat in -90.:90.]
+remapper = Remapper(space, target_hcoords)
+```
+This `remapper` object knows can interpolate `Field`s defined on `space` with the same `interpolate` and `interpolate!` functions. 
+```julia
+import ClimaCore.Fields: coordinate_field
+import ClimaCore.Remapping: interpolate, interpolate!
+
+example_field = coordinate_field(space)
+interpolated_array = interpolate(remapper, example_field) 
+
+# Interpolate in place
+interpolate!(interpolated_array, remapper, example_field)
+```
+
+!!! note
+
+    `interpolate` allocates new arrays and has some internal
+     type-instability, `interpolate!` is non-allocating and type-stable. 
+     When using `interpolate!`, the destination has to be defined on the 
+     root process and has to be `nothing` for the other MPI processes.
+
+
+Multiple fields defined on the same space can be interpolate at the same time
+```julia
+example_field2 = cosd.(example_field)
+interpolated_arrays = interpolate(remapper, [example_field, example_field2]) 
+```
+
+When interpolating multiple fields, greater performance can be achieved by
+creating the `Remapper` with a larger internal buffer to store intermediate
+values for interpolation. Effectively, this controls how many fields can be
+remapped simultaneously in `interpolate`. When more fields than `buffer_length`
+are passed, the remapper will batch the work in sizes of `buffer_length`. The optimal number of fields passed is the `buffer_length` of the `remapper`. If
+more fields are passed, the `remapper` will batch work with size up to its
+`buffer_length`.
+ 
+#### Example
+
+Given `field1`,`field2`, two `Field` defined on a cubed sphere.
+
+```julia
+longpts = range(-180.0, 180.0, 21)
+latpts = range(-80.0, 80.0, 21)
+zpts = range(0.0, 1000.0, 21)
+
+hcoords = [Geometry.LatLongPoint(lat, long) for long in longpts, lat in latpts]
+zcoords = [Geometry.ZPoint(z) for z in zpts]
+
+space = axes(field1)
+
+remapper = Remapper(space, hcoords, zcoords)
+
+int1 = interpolate(remapper, field1)
+int2 = interpolate(remapper, field2)
+
+# Or
+int12 = interpolate(remapper, [field1, field2])
+# With int1 = int12[1, :, :, :]
+```
+
+## Conservative remapping with `TempestRemap`
+
+This section hasn't been written yet. You can help by writing it.
diff --git a/ext/cuda/remapping_distributed.jl b/ext/cuda/remapping_distributed.jl
index 70246c47d4..eb88896fa2 100644
--- a/ext/cuda/remapping_distributed.jl
+++ b/ext/cuda/remapping_distributed.jl
@@ -17,8 +17,14 @@ function _set_interpolated_values_device!(
 )
     # FIXME: Avoid allocation of tuple
     field_values = tuple(map(f -> Fields.field_values(f), fields)...)
-    nblocks, _ = size(interpolation_matrix[1])
-    nthreads = length(vert_interpolation_weights)
+
+    purely_vertical_space = isnothing(interpolation_matrix)
+    num_horizontal_points =
+        purely_vertical_space ? 1 : prod(size(local_horiz_indices))
+    num_points = num_horizontal_points * length(vert_interpolation_weights)
+    max_threads = 256
+    nthreads = min(num_points, max_threads)
+    nblocks = cld(num_points, nthreads)
     args = (
         out,
         interpolation_matrix,
@@ -131,6 +137,41 @@ function set_interpolated_values_kernel!(
     return nothing
 end
 
+# GPU, vertical case
+function set_interpolated_values_kernel!(
+    out::AbstractArray,
+    ::Nothing,
+    ::Nothing,
+    vert_interpolation_weights,
+    vert_bounding_indices,
+    field_values,
+)
+    # TODO: Check the memory access pattern. This was not optimized and likely inefficient!
+    num_fields = length(field_values)
+    num_vert = length(vert_bounding_indices)
+
+    vindex = (blockIdx().y - Int32(1)) * blockDim().y + threadIdx().y
+    findex = (blockIdx().z - Int32(1)) * blockDim().z + threadIdx().z
+
+    totalThreadsY = gridDim().y * blockDim().y
+    totalThreadsZ = gridDim().z * blockDim().z
+
+    CI = CartesianIndex
+    for j in vindex:totalThreadsY:num_vert
+        v_lo, v_hi = vert_bounding_indices[j]
+        A, B = vert_interpolation_weights[j]
+        for k in findex:totalThreadsZ:num_fields
+            if j ≤ num_vert && k ≤ num_fields
+                out[j, k] = (
+                    A * field_values[k][CI(1, 1, 1, v_lo, 1)] +
+                    B * field_values[k][CI(1, 1, 1, v_hi, 1)]
+                )
+            end
+        end
+    end
+    return nothing
+end
+
 function _set_interpolated_values_device!(
     out::AbstractArray,
     fields::AbstractArray{<:Fields.Field},
diff --git a/src/Remapping/Remapping.jl b/src/Remapping/Remapping.jl
index 9639c35932..e3e01a8bec 100644
--- a/src/Remapping/Remapping.jl
+++ b/src/Remapping/Remapping.jl
@@ -5,6 +5,7 @@ using LinearAlgebra, StaticArrays
 import ClimaComms
 import ..DataLayouts,
     ..Geometry,
+    ..Domains,
     ..Spaces,
     ..Grids,
     ..Topologies,
diff --git a/src/Remapping/distributed_remapping.jl b/src/Remapping/distributed_remapping.jl
index 5f22c9636d..b5eeba7f37 100644
--- a/src/Remapping/distributed_remapping.jl
+++ b/src/Remapping/distributed_remapping.jl
@@ -78,6 +78,9 @@
 #   process-local interpolate points in the correct shape with respect to the global
 #   interpolation (and where to collect results)
 #
+# Horizontal and vertical interpolation can be switch off, so that we interpolate purely
+# horizontal/vertical Fields.
+#
 # To process multiple Fields at the same time, some of the scratch spaces gain an extra
 # dimension (buffer_length). With this extra dimension, we can batch the work and process up
 # to buffer_length fields at the same time. This reduces the number of kernel launches and
@@ -116,62 +119,72 @@ function target_hcoords_pid_bitmask(target_hcoords, topology, pid)
     return pid_hcoord.(target_hcoords) .== pid
 end
 
+
+# TODO: Define an inner construct and restrict types, as was done in
+#       https://github.com/CliMA/RRTMGP.jl/pull/352
+#       to avoid potential compilation issues.
 struct Remapper{
     CC <: ClimaComms.AbstractCommsContext,
     SPACE <: Spaces.AbstractSpace,
-    T1 <: AbstractArray,
+    T1, # <: Union{AbstractArray, Nothing},
     TARG_Z <: Union{Nothing, AA1} where {AA1 <: AbstractArray},
-    T3 <: AbstractArray,
-    T4 <: Tuple,
-    T5 <: AbstractArray,
+    T3, # <: Union{AbstractArray, Nothing},
+    T4, # <: Union{Tuple, Nothing},
+    T5, # <: Union{AbstractArray, Nothing},
     VERT_W <: Union{Nothing, AA2} where {AA2 <: AbstractArray},
     VERT_IND <: Union{Nothing, AA3} where {AA3 <: AbstractArray},
-    T8 <: AbstractArray,
-    T9 <: AbstractArray,
+    T8, # <: AbstractArray,
+    T9, # <: AbstractArray,
     T10 <: AbstractArray,
     T11 <: Union{Tuple{Colon}, Tuple{Colon, Colon}, Tuple{Colon, Colon, Colon}},
 }
-
+    # The ClimaComms context
     comms_ctx::CC
 
     # Space over which the remapper is defined
     space::SPACE
 
-    # Target points that are on the process where this object is defined
+    # Target points that are on the process where this object is defined.
     # local_target_hcoords is stored as a 1D array (we store it as 1D array because in
     # general there is no structure to this object, especially for cubed sphere, which have
-    # points spread all over the place)
+    # points spread all over the place). This is nothing when remapping purely vertical
+    # spaces.
     local_target_hcoords::T1
 
-    # Target coordinates in the vertical direction. zcoords are the same for all the processes.
-    # target_zcoords are always assumed to be "reference z"s.
+    # Target coordinates in the vertical direction. zcoords are the same for all the
+    # processes. target_zcoords are always assumed to be "reference z"s. is nothing when
+    # remapping purely horizontal spaces.
     target_zcoords::TARG_Z
 
     # bitmask that identifies which target points are on process where the object is
     # defined. In general, local_target_points_bitmask is a 2D matrix and it is used to fill
     # the correct values in the final output. Every time we find a 1, we are going to stick
-    # the vertical column of the interpolated data.
+    # the vertical column of the interpolated data. This is nothing when remapping purely
+    # vertical spaces.
     local_target_hcoords_bitmask::T3
 
     # Tuple of arrays of weights that performs horizontal interpolation (fixed the
     # horizontal and target spaces). It contains 1 element for grids with one horizontal
-    # dimension, and 2 elements for grids with two horizontal dimensions.
+    # dimension, and 2 elements for grids with two horizontal dimensions. This is nothing
+    # when remapping purely vertical spaces.
     local_horiz_interpolation_weights::T4
 
     # Local indices the element of each local_target_hcoords in the given topology. This is
-    # a linear array with the same length as local_target_hcoords.
+    # a linear array with the same length as local_target_hcoords. This is nothing when
+    # remapping purely vertical spaces.
     local_horiz_indices::T5
 
     # Given the target_zcoords, vert_reference_coordinates contains the reference coordinate
     # in the element. For center spaces, the reference coordinate is shifted to be in (0, 1)
     # when the point is in the lower half of the cell, and in (-1, 0) when it is in the
     # upper half. This shift is needed to directly use the reference coordinate in linear
-    # vertical interpolation. Array of tuples or Nothing.
+    # vertical interpolation. Array of tuples or Nothing. This is nothing when remapping
+    # purely horizontal spaces.
     vert_interpolation_weights::VERT_W
 
     # Given the target_zcoords, vert_bounding_indices contain the vertical indices of the
-    # neighboring elements that are required for vertical interpolation.
-    # Array of tuples or Nothing.
+    # neighboring elements that are required for vertical interpolation. Array of tuples or
+    # Nothing. This is nothing when remapping purely horizontal spaces.
     vert_bounding_indices::VERT_IND
 
     # Scratch space where we save the process-local interpolated values. We keep overwriting
@@ -181,8 +194,8 @@ struct Remapper{
 
     # Scratch space where we save the process-local field value. We keep overwriting this to
     # avoid extra allocations. Ideally, we wouldn't need this and we would use views for
-    # everything. This has dimensions (Nq, ) or (Nq, Nq, )
-    # depending if the horizontal space is 1D or 2D.
+    # everything. This has dimensions (Nq, ) or (Nq, Nq, ) depending if the horizontal space
+    # is 1D or 2D. This is nothing when remapping purely vertical spaces.
     _field_values::T9
 
     # Storage area where the interpolated values are saved. This is meaningful only for the
@@ -200,19 +213,24 @@ struct Remapper{
 end
 
 """
-   Remapper(space, target_hcoords, target_zcoords; buffer_length = 1)
+   Remapper(space, target_hcoords, target_zcoords, buffer_length = 1)
+   Remapper(space; target_hcoords, target_zcoords, buffer_length = 1)
    Remapper(space, target_hcoords; buffer_length = 1)
+   Remapper(space, target_zcoords; buffer_length = 1)
 
 Return a `Remapper` responsible for interpolating any `Field` defined on the given `space`
 to the Cartesian product of `target_hcoords` with `target_zcoords`.
 
-`target_zcoords` can be `nothing` for interpolation on horizontal spaces.
+`target_zcoords` can be `nothing` for interpolation on horizontal spaces. Similarly,
+`target_hcoords` can be `nothing` for interpolation on vertical spaces.
 
 The `Remapper` is designed to not be tied to any particular `Field`. You can use the same
 `Remapper` for any `Field` as long as they are all defined on the same `topology`.
 
 `Remapper` is the main argument to the `interpolate` function.
 
+If you want to quickly remap something, you can call directly [`interpolate`](@ref).
+
 Keyword arguments
 =================
 
@@ -220,15 +238,42 @@ Keyword arguments
 for interpolation. Effectively, this controls how many fields can be remapped simultaneously
 in `interpolate`. When more fields than `buffer_length` are passed, the remapper will batch
 the work in sizes of `buffer_length`.
-
 """
-function Remapper(
+function Remapper end
+
+Remapper(
+    space::Spaces.AbstractSpace;
+    target_hcoords::Union{AbstractArray, Nothing} = nothing,
+    target_zcoords::Union{AbstractArray, Nothing} = nothing,
+    buffer_length::Int = 1,
+) = _Remapper(space; target_zcoords, target_hcoords, buffer_length)
+
+Remapper(
     space::Spaces.AbstractSpace,
-    target_hcoords::AbstractArray,
+    target_hcoords::Union{AbstractArray, Nothing},
     target_zcoords::Union{AbstractArray, Nothing};
     buffer_length::Int = 1,
-)
+) = _Remapper(space; target_hcoords, target_zcoords, buffer_length)
 
+Remapper(
+    space::Spaces.FiniteDifferenceSpace;
+    target_zcoords::AbstractArray,
+    buffer_length::Int = 1,
+) = _Remapper(space; target_zcoords, target_hcoords = nothing, buffer_length)
+
+Remapper(
+    space::Spaces.AbstractSpectralElementSpace,
+    target_hcoords::AbstractArray;
+    buffer_length::Int = 1,
+) = _Remapper(space; target_zcoords = nothing, target_hcoords, buffer_length)
+
+# Constructor for the case with horizontal spaces
+function _Remapper(
+    space::Spaces.AbstractSpace;
+    target_zcoords::Union{AbstractArray, Nothing},
+    target_hcoords::AbstractArray,
+    buffer_length::Int = 1,
+)
     comms_ctx = ClimaComms.context(space)
     pid = ClimaComms.mypid(comms_ctx)
     FT = Spaces.undertype(space)
@@ -367,11 +412,45 @@ function Remapper(
     )
 end
 
-Remapper(
-    space::Spaces.AbstractSpace,
-    target_hcoords::AbstractArray;
+# Constructor for the case with vertical spaces
+function _Remapper(
+    space::Spaces.FiniteDifferenceSpace;
+    target_zcoords::AbstractArray,
+    target_hcoords::Nothing,
     buffer_length::Int = 1,
-) = Remapper(space, target_hcoords, nothing; buffer_length)
+)
+    comms_ctx = ClimaComms.context(space)
+    FT = Spaces.undertype(space)
+    ArrayType = ClimaComms.array_type(space)
+
+    vert_interpolation_weights =
+        ArrayType(vertical_interpolation_weights(space, target_zcoords))
+    vert_bounding_indices =
+        ArrayType(vertical_bounding_indices(space, target_zcoords))
+
+    local_interpolated_values =
+        ArrayType(zeros(FT, (length(target_zcoords), buffer_length)))
+    interpolated_values =
+        ArrayType(zeros(FT, (length(target_zcoords), buffer_length)))
+    colons = (:,)
+
+    return Remapper(
+        comms_ctx,
+        space,
+        nothing, # local_target_hcoords,
+        target_zcoords,
+        nothing, # local_target_hcoords_bitmask,
+        nothing, # local_horiz_interpolation_weights,
+        nothing, # local_horiz_indices,
+        vert_interpolation_weights,
+        vert_bounding_indices,
+        local_interpolated_values,
+        nothing, # field_values,
+        interpolated_values,
+        buffer_length,
+        colons,
+    )
+end
 
 """
     _set_interpolated_values!(remapper, field)
@@ -439,6 +518,38 @@ function set_interpolated_values_cpu_kernel!(
     end
 end
 
+# CPU, vertical case
+function set_interpolated_values_cpu_kernel!(
+    out::AbstractArray,
+    fields::AbstractArray{<:Fields.Field},
+    ::Nothing,
+    ::Nothing,
+    vert_interpolation_weights,
+    vert_bounding_indices,
+    ::Nothing,
+)
+    space = axes(first(fields))
+    FT = Spaces.undertype(space)
+    for (field_index, field) in enumerate(fields)
+        field_values = Fields.field_values(field)
+
+        # Reading values from field_values is expensive, so we try to limit the number of reads. We can do
+        # this because multiple target points might be all contained in the same element.
+        prev_vindex = -1
+        @inbounds for (vindex, (A, B)) in enumerate(vert_interpolation_weights)
+            (v_lo, v_hi) = vert_bounding_indices[vindex]
+            # If we are no longer in the same element, read the field values again
+            if prev_vindex != vindex
+                out[vindex, field_index] = (
+                    A * field_values[CartesianIndex(1, 1, 1, v_lo, 1)] +
+                    B * field_values[CartesianIndex(1, 1, 1, v_hi, 1)]
+                )
+                prev_vindex = vindex
+            end
+        end
+    end
+end
+
 # CPU, 2D case
 function set_interpolated_values_cpu_kernel!(
     out::AbstractArray,
@@ -559,7 +670,6 @@ function _set_interpolated_values_device!(
     ::Nothing,
     ::ClimaComms.AbstractDevice,
 )
-
     space = axes(first(fields))
     FT = Spaces.undertype(space)
     quad = Spaces.quadrature_style(space)
@@ -659,7 +769,6 @@ function _collect_interpolated_values!(
     index_field_end::Int;
     only_one_field,
 )
-
     if only_one_field
         ClimaComms.reduce!(
             remapper.comms_ctx,
@@ -757,6 +866,8 @@ function interpolate(remapper::Remapper, fields)
             error("Field is defined on a different space than remapper")
     end
 
+    isa_vertical_space = remapper.space isa Spaces.FiniteDifferenceSpace
+
     index_field_begin, index_field_end =
         1, min(length(fields), remapper.buffer_length)
 
@@ -777,13 +888,18 @@ function interpolate(remapper::Remapper, fields)
             remapper,
             view(fields, index_field_begin:index_field_end),
         )
-        # Reshape the output so that it is a nice grid.
-        _apply_mpi_bitmask!(remapper, num_fields)
+
+        if !isa_vertical_space
+            # For spaces with an horizontal component, reshape the output so that it is a nice grid.
+            _apply_mpi_bitmask!(remapper, num_fields)
+        else
+            # For purely vertical spaces, just move to _interpolated_values
+            remapper._interpolated_values .= remapper._local_interpolated_values
+        end
+
         # Finally, we have to send all the _interpolated_values to root and sum them up to
-        # obtain the final answer. Only the root will contain something useful. This also
-        # moves the data off the GPU
-        ret = _collect_and_return_interpolated_values!(remapper, num_fields)
-        return ret
+        # obtain the final answer. Only the root will contain something useful.
+        return _collect_and_return_interpolated_values!(remapper, num_fields)
     end
 
     # Non-root processes
@@ -793,39 +909,6 @@ function interpolate(remapper::Remapper, fields)
            interpolated_values
 end
 
-"""
-       interpolate(field::ClimaCore.Fields, target_hcoords, target_zcoords)
-
-Interpolate the given fields on the Cartesian product of `target_hcoords` with
-`target_zcoords` (if not empty).
-
-Coordinates have to be `ClimaCore.Geometry.Points`.
-
-Note: do not use this method when performance is important. Instead, define a `Remapper` and
-call `interpolate(remapper, fields)`. Different `Field`s defined on the same `Space` can
-share a `Remapper`, so that interpolation can be optimized.
-
-Example
-========
-
-Given `field`, a `Field` defined on a cubed sphere.
-
-```julia
-longpts = range(-180.0, 180.0, 21)
-latpts = range(-80.0, 80.0, 21)
-zpts = range(0.0, 1000.0, 21)
-
-hcoords = [Geometry.LatLongPoint(lat, long) for long in longpts, lat in latpts]
-zcoords = [Geometry.ZPoint(z) for z in zpts]
-
-interpolate(field, hcoords, zcoords)
-```
-"""
-function interpolate(field::Fields.Field, target_hcoords, target_zcoords)
-    remapper = Remapper(axes(field), target_hcoords, target_zcoords)
-    return interpolate(remapper, field)
-end
-
 # dest has to be allowed to be nothing because interpolation happens only on the root
 # process
 function interpolate!(
@@ -837,6 +920,7 @@ function interpolate!(
     if only_one_field
         fields = [fields]
     end
+    isa_vertical_space = remapper.space isa Spaces.FiniteDifferenceSpace
 
     if !isnothing(dest)
         # !isnothing(dest) means that this is the root process, in this case, the size have
@@ -869,11 +953,17 @@ function interpolate!(
             remapper,
             view(fields, index_field_begin:index_field_end),
         )
-        # Reshape the output so that it is a nice grid.
-        _apply_mpi_bitmask!(remapper, num_fields)
+
+        if !isa_vertical_space
+            # For spaces with an horizontal component, reshape the output so that it is a nice grid.
+            _apply_mpi_bitmask!(remapper, num_fields)
+        else
+            # For purely vertical spaces, just move to _interpolated_values
+            remapper._interpolated_values .= remapper._local_interpolated_values
+        end
+
         # Finally, we have to send all the _interpolated_values to root and sum them up to
-        # obtain the final answer. Only the root will contain something useful. This also
-        # moves the data off the GPU
+        # obtain the final answer.
         _collect_interpolated_values!(
             dest,
             remapper,
@@ -889,3 +979,229 @@ function interpolate!(
     end
     return nothing
 end
+
+"""
+       interpolate(field::ClimaCore.Fields;
+                   hresolution = 180,
+                   zresolution = 50,
+                   target_hcoords = default_target_hcoords(space; hresolution),
+                   target_zcoords = default_target_vcoords(space; zresolution)
+                   )
+
+Interpolate the given fields on the Cartesian product of `target_hcoords` with
+`target_zcoords` (if not empty).
+
+Coordinates have to be `ClimaCore.Geometry.Points`.
+
+Note: do not use this method when performance is important. Instead, define a `Remapper` and
+call `interpolate(remapper, fields)`. Different `Field`s defined on the same `Space` can
+share a `Remapper`, so that interpolation can be optimized.
+
+Example
+========
+
+Given `field`, a `Field` defined on a cubed sphere.
+
+By default, a target uniform grid is chosen (with resolution `hresolution` and
+`zresolution`), so remapping is simply
+```julia
+julia> interpolate(field)
+```
+This will return an array of interpolated values.
+
+Resolution can be specified
+```julia
+julia> interpolate(field; hresolution = 100, zresolution = 50)
+```
+Coordinates can be also specified directly:
+```julia
+julia> longpts = range(-180.0, 180.0, 21)
+julia> latpts = range(-80.0, 80.0, 21)
+julia> zpts = range(0.0, 1000.0, 21)
+
+julia> hcoords = [Geometry.LatLongPoint(lat, long) for long in longpts, lat in latpts]
+julia> zcoords = [Geometry.ZPoint(z) for z in zpts]
+
+julia> interpolate(field, target_hcoords, target_zcoords)
+```
+
+If you need the array of coordinates, you can call
+[`default_target_hcoords`](@ref) (or [`default_target_vcoords`](@ref)) passing
+`axes(field)`. This will return an array of `Geometry.Point`s. The functions
+[`Geometry.Components`](@ref) and [`Geometry.Component`](@ref)  can be used to extract the components as numeric
+values. For example,
+```julia
+julia> Geometry.components.(Geometry.components.([
+           Geometry.LatLongPoint(x, y) for x in range(-180.0, 180.0, length = 180),
+           y in range(-90.0, 90.0, length = 180)
+       ]))
+180×180 Matrix{StaticArraysCore.SVector{2, Float64}}:
+ [-180.0, -90.0]    [-180.0, -88.9944]    …  [-180.0, 88.9944]    [-180.0, 90.0]
+ [-177.989, -90.0]  [-177.989, -88.9944]     [-177.989, 88.9944]  [-177.989, 90.0]
+ [-175.978, -90.0]  [-175.978, -88.9944]     [-175.978, 88.9944]  [-175.978, 90.0]
+ [-173.966, -90.0]  [-173.966, -88.9944]     [-173.966, 88.9944]  [-173.966, 90.0]
+ [-171.955, -90.0]  [-171.955, -88.9944]     [-171.955, 88.9944]  [-171.955, 90.0]
+ [-169.944, -90.0]  [-169.944, -88.9944]  …  [-169.944, 88.9944]  [-169.944, 90.0]
+ [-167.933, -90.0]  [-167.933, -88.9944]     [-167.933, 88.9944]  [-167.933, 90.0]
+ [-165.922, -90.0]  [-165.922, -88.9944]     [-165.922, 88.9944]  [-165.922, 90.0]
+ [-163.911, -90.0]  [-163.911, -88.9944]     [-163.911, 88.9944]  [-163.911, 90.0]
+ [-161.899, -90.0]  [-161.899, -88.9944]     [-161.899, 88.9944]  [-161.899, 90.0]
+ [-159.888, -90.0]  [-159.888, -88.9944]  …  [-159.888, 88.9944]  [-159.888, 90.0]
+ [-157.877, -90.0]  [-157.877, -88.9944]     [-157.877, 88.9944]  [-157.877, 90.0]
+ [-155.866, -90.0]  [-155.866, -88.9944]     [-155.866, 88.9944]  [-155.866, 90.0]
+ [-153.855, -90.0]  [-153.855, -88.9944]     [-153.855, 88.9944]  [-153.855, 90.0]
+ [-151.844, -90.0]  [-151.844, -88.9944]     [-151.844, 88.9944]  [-151.844, 90.0]
+ [-149.832, -90.0]  [-149.832, -88.9944]  …  [-149.832, 88.9944]  [-149.832, 90.0]
+ [-147.821, -90.0]  [-147.821, -88.9944]     [-147.821, 88.9944]  [-147.821, 90.0]
+ [-145.81, -90.0]   [-145.81, -88.9944]      [-145.81, 88.9944]   [-145.81, 90.0]
+ [-143.799, -90.0]  [-143.799, -88.9944]     [-143.799, 88.9944]  [-143.799, 90.0]
+ [-141.788, -90.0]  [-141.788, -88.9944]     [-141.788, 88.9944]  [-141.788, 90.0]
+ [-139.777, -90.0]  [-139.777, -88.9944]  …  [-139.777, 88.9944]  [-139.777, 90.0]
+ [-137.765, -90.0]  [-137.765, -88.9944]     [-137.765, 88.9944]  [-137.765, 90.0]
+ ⋮                                        ⋱
+ [137.765, -90.0]   [137.765, -88.9944]      [137.765, 88.9944]   [137.765, 90.0]
+ [139.777, -90.0]   [139.777, -88.9944]      [139.777, 88.9944]   [139.777, 90.0]
+ [141.788, -90.0]   [141.788, -88.9944]   …  [141.788, 88.9944]   [141.788, 90.0]
+ [143.799, -90.0]   [143.799, -88.9944]      [143.799, 88.9944]   [143.799, 90.0]
+ [145.81, -90.0]    [145.81, -88.9944]       [145.81, 88.9944]    [145.81, 90.0]
+ [147.821, -90.0]   [147.821, -88.9944]      [147.821, 88.9944]   [147.821, 90.0]
+ [149.832, -90.0]   [149.832, -88.9944]      [149.832, 88.9944]   [149.832, 90.0]
+ [151.844, -90.0]   [151.844, -88.9944]   …  [151.844, 88.9944]   [151.844, 90.0]
+ [153.855, -90.0]   [153.855, -88.9944]      [153.855, 88.9944]   [153.855, 90.0]
+ [155.866, -90.0]   [155.866, -88.9944]      [155.866, 88.9944]   [155.866, 90.0]
+ [157.877, -90.0]   [157.877, -88.9944]      [157.877, 88.9944]   [157.877, 90.0]
+ [159.888, -90.0]   [159.888, -88.9944]      [159.888, 88.9944]   [159.888, 90.0]
+ [161.899, -90.0]   [161.899, -88.9944]   …  [161.899, 88.9944]   [161.899, 90.0]
+ [163.911, -90.0]   [163.911, -88.9944]      [163.911, 88.9944]   [163.911, 90.0]
+ [165.922, -90.0]   [165.922, -88.9944]      [165.922, 88.9944]   [165.922, 90.0]
+ [167.933, -90.0]   [167.933, -88.9944]      [167.933, 88.9944]   [167.933, 90.0]
+ [169.944, -90.0]   [169.944, -88.9944]      [169.944, 88.9944]   [169.944, 90.0]
+ [171.955, -90.0]   [171.955, -88.9944]   …  [171.955, 88.9944]   [171.955, 90.0]
+ [173.966, -90.0]   [173.966, -88.9944]      [173.966, 88.9944]   [173.966, 90.0]
+ [175.978, -90.0]   [175.978, -88.9944]      [175.978, 88.9944]   [175.978, 90.0]
+ [177.989, -90.0]   [177.989, -88.9944]      [177.989, 88.9944]   [177.989, 90.0]
+ [180.0, -90.0]     [180.0, -88.9944]        [180.0, 88.9944]     [180.0, 90.0]
+```
+To extract only long or lat, one can broadcast `getindex`
+```julia
+julia> lats = getindex.(Geometry.components.([Geometry.LatLongPoint(x, y)
+                                              for x in range(-180.0, 180.0, length = 180),
+                                                  y in range(-90.0, 90.0, length = 180)
+                                             ]),
+                        1)
+```
+This can be used directly for plotting.
+"""
+function interpolate(
+    field::Fields.Field;
+    zresolution = 50,
+    hresolution = 180,
+    target_hcoords = default_target_hcoords(axes(field); hresolution),
+    target_zcoords = default_target_zcoords(axes(field); zresolution),
+)
+    return interpolate(field, axes(field); hresolution, zresolution)
+end
+
+function interpolate(field::Fields.Field, target_hcoords, target_zcoords)
+    remapper = Remapper(axes(field), target_hcoords, target_zcoords)
+    return interpolate(remapper, field)
+end
+
+"""
+    default_target_hcoords(space::Spaces.AbstractSpace; hresolution)
+
+Return an Array with the Geometry.Points to interpolate uniformly the horizontal
+component of the given `space`.
+"""
+function default_target_hcoords(space::Spaces.AbstractSpace; hresolution = 180)
+    return default_target_hcoords(Spaces.horizontal_space(space); hresolution)
+end
+
+"""
+    default_target_hcoords_as_vectors(space::Spaces.AbstractSpace; hresolution)
+
+Return an Vectors with the coordinate to interpolate uniformly the horizontal
+component of the given `space`.
+"""
+function default_target_hcoords_as_vectors(
+    space::Spaces.AbstractSpace;
+    hresolution = 180,
+)
+    return default_target_hcoords_as_vectors(
+        Spaces.horizontal_space(space);
+        hresolution,
+    )
+end
+
+function default_target_hcoords(
+    space::Spaces.SpectralElementSpace2D;
+    hresolution = 180,
+)
+    topology = Spaces.topology(space)
+    domain = Meshes.domain(topology.mesh)
+    xrange, yrange = default_target_hcoords_as_vectors(space; hresolution)
+    PointType =
+        domain isa Domains.SphereDomain ? Geometry.LatLongPoint :
+        Topologies.coordinate_type(topology)
+    return [PointType(x, y) for x in xrange, y in yrange]
+end
+
+function default_target_hcoords_as_vectors(
+    space::Spaces.SpectralElementSpace2D;
+    hresolution = 180,
+)
+    FT = Spaces.undertype(space)
+    topology = Spaces.topology(space)
+    domain = Meshes.domain(topology.mesh)
+    if domain isa Domains.SphereDomain
+        return FT.(range(-180.0, 180.0, hresolution)),
+        FT.(range(-90.0, 90.0, hresolution))
+    else
+        x1min = Geometry.component(domain.interval1.coord_min, 1)
+        x2min = Geometry.component(domain.interval2.coord_min, 1)
+        x1max = Geometry.component(domain.interval1.coord_max, 1)
+        x2max = Geometry.component(domain.interval2.coord_max, 1)
+        return collect(range(x1min, x1max, hresolution)),
+        collect(range(x2min, x2max, hresolution))
+    end
+end
+
+function default_target_hcoords(
+    space::Spaces.SpectralElementSpace1D;
+    hresolution = 180,
+)
+    topology = Spaces.topology(space)
+    PointType = Topologies.coordinate_type(topology)
+    return PointType.(default_target_hcoords_as_vectors(space; hresolution))
+end
+
+function default_target_hcoords_as_vectors(
+    space::Spaces.SpectralElementSpace1D;
+    hresolution = 180,
+)
+    topology = Spaces.topology(space)
+    domain = Meshes.domain(topology.mesh)
+    xmin = Geometry.component(domain.coord_min, 1)
+    xmax = Geometry.component(domain.coord_max, 1)
+    return collect(range(xmin, xmax, hresolution))
+end
+
+
+"""
+    default_target_zcoords(space::Spaces.AbstractSpace; zresolution)
+
+Return an Array with the Geometry.Points to interpolate uniformly the vertical component of
+the given `space`.
+"""
+function default_target_zcoords(space; zresolution = 50)
+    return Geometry.ZPoint.(
+        default_target_zcoords_as_vectors(space; zresolution)
+    )
+
+end
+
+function default_target_zcoords_as_vectors(space; zresolution = 50)
+    return collect(
+        range(Domains.z_min(space), Domains.z_max(space), zresolution),
+    )
+end
diff --git a/src/Remapping/interpolate_array.jl b/src/Remapping/interpolate_array.jl
index 67d41df3d8..83935e0bfd 100644
--- a/src/Remapping/interpolate_array.jl
+++ b/src/Remapping/interpolate_array.jl
@@ -120,6 +120,7 @@ function vertical_indices_ref_coordinate(
         end
         ξ3 = ξ3 - 1
     end
+    @show v_lo, v_hi, ξ3
     return v_lo, v_hi, ξ3
 end
 
diff --git a/src/Remapping/remapping_utils.jl b/src/Remapping/remapping_utils.jl
index 84122836d5..23f461ffca 100644
--- a/src/Remapping/remapping_utils.jl
+++ b/src/Remapping/remapping_utils.jl
@@ -27,7 +27,9 @@ function vertical_reference_coordinates(space, zcoords)
     vert_topology = Spaces.vertical_topology(space)
     vert_mesh = vert_topology.mesh
 
-    is_cell_center = space isa Spaces.CenterExtrudedFiniteDifferenceSpace
+    is_cell_center =
+        space isa Spaces.CenterExtrudedFiniteDifferenceSpace ||
+        space isa Spaces.CenterFiniteDifferenceSpace
 
     ξ3s = map(zcoords) do zcoord
         velem = Meshes.containing_element(vert_mesh, zcoord)
@@ -51,10 +53,13 @@ element in a column, no interpolation is performed and the value at the cell cen
 returned. Effectively, this means that the interpolation is first-order accurate across the
 column, but zeroth-order accurate close to the boundaries.
 """
-function vertical_bounding_indices(space, zcoords) end
+function vertical_bounding_indices end
 
 function vertical_bounding_indices(
-    space::Spaces.FaceExtrudedFiniteDifferenceSpace,
+    space::Union{
+        Spaces.FaceExtrudedFiniteDifferenceSpace,
+        Spaces.FaceFiniteDifferenceSpace,
+    },
     zcoords,
 )
     vert_topology = Spaces.vertical_topology(space)
@@ -64,7 +69,10 @@ function vertical_bounding_indices(
 end
 
 function vertical_bounding_indices(
-    space::Spaces.CenterExtrudedFiniteDifferenceSpace,
+    space::Union{
+        Spaces.CenterExtrudedFiniteDifferenceSpace,
+        Spaces.CenterFiniteDifferenceSpace,
+    },
     zcoords,
 )
     vert_topology = Spaces.vertical_topology(space)
diff --git a/test/MatrixFields/field_names.jl b/test/MatrixFields/field_names.jl
index 66dde985f6..305cf18a7b 100644
--- a/test/MatrixFields/field_names.jl
+++ b/test/MatrixFields/field_names.jl
@@ -45,7 +45,7 @@ const x = (; foo = Foo(0), a = (; b = 1, c = ((; d = 2), (;), ((), nothing))))
     @test_all MatrixFields.broadcasted_has_field(typeof(x), @name(a.c.:(1).d))
     @test_all !MatrixFields.broadcasted_has_field(
         typeof(x),
-        @name(foo.invalid_name),
+        @name(foo.invalid_name)
     )
 
     @test_all MatrixFields.broadcasted_get_field(x, @name()) == x
@@ -65,7 +65,7 @@ const x = (; foo = Foo(0), a = (; b = 1, c = ((; d = 2), (;), ((), nothing))))
               @name(c.:(1).d)
     @test_throws "is not a child name" MatrixFields.extract_internal_name(
         @name(a.c.:(1).d),
-        @name(foo),
+        @name(foo)
     )
 
     @test_all MatrixFields.append_internal_name(@name(a), @name(c.:(1).d)) ==
@@ -125,21 +125,21 @@ end
 
     @testset "FieldNameSet Constructors" begin
         @test_throws "Invalid FieldNameSet value" vector_keys(
-            @name(invalid_name),
+            @name(invalid_name)
         )
         @test_throws "Invalid FieldNameSet value" matrix_keys((
             @name(invalid_name),
-            @name(a.c),
+            @name(a.c)
         ),)
 
         for constructor in (vector_keys, vector_keys_no_tree)
             @test_throws "Duplicate FieldNameSet values" constructor(
                 @name(foo),
-                @name(foo),
+                @name(foo)
             )
             @test_throws "Overlapping FieldNameSet values" constructor(
                 @name(foo),
-                @name(foo.value),
+                @name(foo.value)
             )
         end
         for constructor in (matrix_keys, matrix_keys_no_tree)
@@ -166,7 +166,7 @@ end
         @name(foo.value),
         @name(a.c.:(1)),
         @name(a.c.:(2)),
-        @name(a.c.:(3)),
+        @name(a.c.:(3))
     )
     m_set3 = matrix_keys(
         (@name(foo.value), @name(a.c.:(1))),
@@ -376,7 +376,7 @@ end
                               @name(a.b),
                               @name(a.c.:(1)),
                               @name(a.c.:(2)),
-                              @name(a.c.:(3)),
+                              @name(a.c.:(3))
                           )
                 @test_all setdiff(set1, set4) ==
                           vector_keys_no_tree(@name(foo), @name(a.c.:(3)))
@@ -511,7 +511,7 @@ end
         ) == vector_keys_no_tree(
             @name(a.c.:(1)),
             @name(a.c.:(2)),
-            @name(a.c.:(3)),
+            @name(a.c.:(3))
         )
         @test_all MatrixFields.matrix_product_keys(
             matrix_keys_no_tree(
diff --git a/test/Remapping/distributed_remapping.jl b/test/Remapping/distributed_remapping.jl
index 0748934dc4..e3c5fdc6b5 100644
--- a/test/Remapping/distributed_remapping.jl
+++ b/test/Remapping/distributed_remapping.jl
@@ -1,3 +1,7 @@
+#=
+julia --project
+using Revise; include("test/Remapping/distributed_remapping.jl")
+=#
 using Logging
 using Test
 using IntervalSets
@@ -36,7 +40,7 @@ end
 end
 
 on_gpu = device isa ClimaComms.CUDADevice
-broken = false
+with_mpi = context isa ClimaComms.MPICommsContext
 
 if !on_gpu
     @testset "2D extruded" begin
@@ -657,3 +661,238 @@ end
         @test interp_sin_long ≈ dest[:, :, 3]
     end
 end
+
+if !with_mpi
+    @testset "Purely vertical space" begin
+        vertdomain = Domains.IntervalDomain(
+            Geometry.ZPoint(0.0),
+            Geometry.ZPoint(1000.0);
+            boundary_names = (:bottom, :top),
+        )
+
+        vertmesh = Meshes.IntervalMesh(vertdomain, nelems = 30)
+        verttopo = Topologies.IntervalTopology(
+            ClimaComms.SingletonCommsContext(ClimaComms.device()),
+            vertmesh,
+        )
+        cspace = Spaces.CenterFiniteDifferenceSpace(verttopo)
+        fspace = Spaces.FaceFiniteDifferenceSpace(cspace)
+
+        zpts = range(0.0, 1000.0, 21)
+        zcoords = [Geometry.ZPoint(z) for z in zpts]
+
+        # Center space
+        cremapper = Remapping.Remapper(
+            cspace;
+            target_zcoords = zcoords,
+            buffer_length = 2,
+        )
+        ccoords = Fields.coordinate_field(cspace)
+        cinterp_z = Remapping.interpolate(cremapper, ccoords.z)
+        cexpected_z = ArrayType(zpts)
+        ClimaComms.allowscalar(device) do
+            cexpected_z[1] = 1000.0 * (0 / 30 + 1 / 30) / 2
+            cexpected_z[end] = 1000.0 * (29 / 30 + 30 / 30) / 2
+        end
+        @test cinterp_z ≈ cexpected_z
+
+        # Face space
+        fremapper = Remapping.Remapper(
+            fspace;
+            target_zcoords = zcoords,
+            buffer_length = 2,
+        )
+        fcoords = Fields.coordinate_field(fspace)
+        finterp_z = Remapping.interpolate(fremapper, fcoords.z)
+        fexpected_z = ArrayType(zpts)
+        finterp_z ≈ fexpected_z
+
+        # Remapping two fields
+        cinterp = Remapping.interpolate(cremapper, [ccoords.z, ccoords.z])
+        @test cexpected_z ≈ cinterp[:, 1]
+        @test cexpected_z ≈ cinterp[:, 2]
+
+        finterp = Remapping.interpolate(fremapper, [fcoords.z, fcoords.z])
+        @test fexpected_z ≈ finterp[:, 1]
+        @test fexpected_z ≈ finterp[:, 2]
+
+        # Remapping three fields (more than the buffer length)
+        cinterp =
+            Remapping.interpolate(cremapper, [ccoords.z, ccoords.z, ccoords.z])
+        @test cexpected_z ≈ cinterp[:, 1]
+        @test cexpected_z ≈ cinterp[:, 2]
+        @test cexpected_z ≈ cinterp[:, 3]
+
+        finterp =
+            Remapping.interpolate(fremapper, [fcoords.z, fcoords.z, fcoords.z])
+        @test fexpected_z ≈ finterp[:, 1]
+        @test fexpected_z ≈ finterp[:, 2]
+        @test fexpected_z ≈ finterp[:, 3]
+
+        # Remapping in-place one field
+        dest = ArrayType(zeros(21))
+        Remapping.interpolate!(dest, cremapper, ccoords.z)
+        @test cexpected_z ≈ dest
+
+        Remapping.interpolate!(dest, fremapper, fcoords.z)
+        @test fexpected_z ≈ dest
+
+        # Three fields (more than buffer length)
+        dest = ArrayType(zeros(21, 3))
+        Remapping.interpolate!(
+            dest,
+            cremapper,
+            [ccoords.z, ccoords.z, ccoords.z],
+        )
+        @test cexpected_z ≈ dest[:, 1]
+        @test cexpected_z ≈ dest[:, 2]
+        @test cexpected_z ≈ dest[:, 3]
+
+        Remapping.interpolate!(
+            dest,
+            fremapper,
+            [fcoords.z, fcoords.z, fcoords.z],
+        )
+        @test fexpected_z ≈ dest[:, 1]
+        @test fexpected_z ≈ dest[:, 2]
+        @test fexpected_z ≈ dest[:, 3]
+    end
+
+end
+
+@testset "Convenience interpolate" begin
+
+    # 3D Sphere space
+    vertdomain = Domains.IntervalDomain(
+        Geometry.ZPoint(0.0),
+        Geometry.ZPoint(1000.0);
+        boundary_names = (:bottom, :top),
+    )
+
+    vertmesh = Meshes.IntervalMesh(vertdomain, nelems = 30)
+    verttopo = Topologies.IntervalTopology(
+        ClimaComms.SingletonCommsContext(ClimaComms.device()),
+        vertmesh,
+    )
+    vert_center_space = Spaces.CenterFiniteDifferenceSpace(verttopo)
+
+    horzdomain = Domains.SphereDomain(1e6)
+
+    quad = Quadratures.GLL{4}()
+    horzmesh = Meshes.EquiangularCubedSphere(horzdomain, 6)
+    horztopology = Topologies.Topology2D(context, horzmesh)
+    horzspace = Spaces.SpectralElementSpace2D(horztopology, quad)
+
+    hv_center_space =
+        Spaces.ExtrudedFiniteDifferenceSpace(horzspace, vert_center_space)
+
+    @test all(
+        Remapping.default_target_zcoords(hv_center_space) .≈
+        Geometry.ZPoint.(range(0.0, 1000; length = 50)),
+    )
+
+    @test Remapping.default_target_hcoords(hv_center_space) == [
+        Geometry.LatLongPoint(x, y) for x in range(-180.0, 180.0, length = 180),
+        y in range(-90.0, 90.0, length = 180)
+    ]
+
+    # Purely horizontal 2D space sphere
+    @test Remapping.default_target_hcoords(horzspace) == [
+        Geometry.LatLongPoint(x, y) for x in range(-180.0, 180.0, length = 180),
+        y in range(-90.0, 90.0, length = 180)
+    ]
+
+    # Purely vertical spaces
+
+    @test all(
+        Remapping.default_target_zcoords(vert_center_space) .≈
+        Geometry.ZPoint.(range(0.0, 1000; length = 50)),
+    )
+    @test all(
+        Remapping.default_target_zcoords(
+            Spaces.FaceFiniteDifferenceSpace(vert_center_space),
+        ) .≈ Geometry.ZPoint.(range(0.0, 1000; length = 50)),
+    )
+
+    # 3D box space
+
+    vertdomain = Domains.IntervalDomain(
+        Geometry.ZPoint(0.0),
+        Geometry.ZPoint(1000.0);
+        boundary_names = (:bottom, :top),
+    )
+
+    vertmesh = Meshes.IntervalMesh(vertdomain, nelems = 30)
+    verttopo = Topologies.IntervalTopology(
+        ClimaComms.SingletonCommsContext(device),
+        vertmesh,
+    )
+    vert_center_space = Spaces.CenterFiniteDifferenceSpace(verttopo)
+
+    horzdomain = Domains.RectangleDomain(
+        Geometry.XPoint(-500.0) .. Geometry.XPoint(500.0),
+        Geometry.YPoint(-300.0) .. Geometry.YPoint(300.0),
+        x1periodic = true,
+        x2periodic = true,
+    )
+
+    quad = Quadratures.GLL{4}()
+    horzmesh = Meshes.RectilinearMesh(horzdomain, 10, 10)
+    horztopology = Topologies.Topology2D(
+        ClimaComms.SingletonCommsContext(device),
+        horzmesh,
+    )
+    horzspace = Spaces.SpectralElementSpace2D(horztopology, quad)
+
+    hv_center_space =
+        Spaces.ExtrudedFiniteDifferenceSpace(horzspace, vert_center_space)
+
+    @test all(
+        Remapping.default_target_hcoords(hv_center_space) .≈ [
+            Geometry.XYPoint(x, y) for x in range(-500.0, 500.0, length = 180),
+            y in range(-300.0, 300.0, length = 180)
+        ],
+    )
+
+    # Purely horizontal 2D space box
+    @test all(
+        Remapping.default_target_hcoords(horzspace) .≈ [
+            Geometry.XYPoint(x, y) for x in range(-500.0, 500.0, length = 180),
+            y in range(-300.0, 300.0, length = 180)
+        ],
+    )
+
+    # 2D slice space
+    if !on_gpu
+        horzdomain = Domains.IntervalDomain(
+            Geometry.XPoint(-500.0) .. Geometry.XPoint(500.0),
+            periodic = true,
+        )
+        quad = Quadratures.GLL{4}()
+        horzmesh = Meshes.IntervalMesh(horzdomain, nelems = 10)
+        horztopology = Topologies.IntervalTopology(
+            ClimaComms.SingletonCommsContext(device),
+            horzmesh,
+        )
+        horzspace = Spaces.SpectralElementSpace1D(horztopology, quad)
+
+        hv_center_space =
+            Spaces.ExtrudedFiniteDifferenceSpace(horzspace, vert_center_space)
+
+        @test all(
+            Remapping.default_target_hcoords(hv_center_space) .≈
+            [Geometry.XPoint(x) for x in range(-500.0, 500.0, length = 180)],
+        )
+
+        @test all(
+            Remapping.default_target_zcoords(hv_center_space) .≈
+            Geometry.ZPoint.(range(0.0, 1000; length = 50)),
+        )
+
+        # Purely horizontal 1D space
+        @test all(
+            Remapping.default_target_hcoords(horzspace) .≈
+            [Geometry.XPoint(x) for x in range(-500.0, 500.0, length = 180)],
+        )
+    end
+end
diff --git a/test/Remapping/interpolate_array.jl b/test/Remapping/interpolate_array.jl
index 948343e047..eaabd28e61 100644
--- a/test/Remapping/interpolate_array.jl
+++ b/test/Remapping/interpolate_array.jl
@@ -47,9 +47,8 @@ device = ClimaComms.CPUSingleThreaded()
 
     coords = Fields.coordinate_field(hv_center_space)
 
-    xpts = range(Geometry.XPoint(-500.0), Geometry.XPoint(500.0), length = 21)
-    zpts = range(Geometry.ZPoint(0.0), Geometry.ZPoint(1000.0), length = 21)
-
+    xpts = range(Geometry.XPoint(-500.0), Geometry.XPoint(500.0), length = 1)
+    zpts = range(Geometry.ZPoint(0.0), Geometry.ZPoint(1000.0), length = 3)
 
     interp_x = Remapping.interpolate_array(coords.x, xpts, zpts)
     @test interp_x ≈ [x.x for x in xpts, z in zpts]