diff --git a/src/Hypsography/Hypsography.jl b/src/Hypsography/Hypsography.jl
index 78603a00da..833e049edc 100644
--- a/src/Hypsography/Hypsography.jl
+++ b/src/Hypsography/Hypsography.jl
@@ -23,10 +23,10 @@ end
 
 Locate vertical levels using an exponential function between the surface field and the top
 of the domain, using the method of [Schar2002](@cite). This method is modified
-such no warping is applied above some user defined parameter 0 ≤ ηₕ < 1.0, where the lower and upper 
-bounds represent the domain bottom and top respectively. `s` governs the decay rate. 
-If the decay-scale is poorly specified (i.e., `s * zₜ` is lower than the maximum 
-surface elevation), a warning is thrown and `s` is adjusted such that it `szₜ > maximum(z_surface)`. 
+such no warping is applied above some user defined parameter 0 ≤ ηₕ < 1.0, where the lower and upper
+bounds represent the domain bottom and top respectively. `s` governs the decay rate.
+If the decay-scale is poorly specified (i.e., `s * zₜ` is lower than the maximum
+surface elevation), a warning is thrown and `s` is adjusted such that it `szₜ > maximum(z_surface)`.
 """
 struct SLEVEAdaption{F <: Fields.Field, FT <: Real} <: HypsographyAdaption
     # Union can be removed once deprecation removed.
@@ -201,15 +201,15 @@ end
 
 Option for 2nd order diffusive smoothing of generated terrain.
 Mutate (smooth) a given elevation profile `f` before assigning the surface
-elevation to the `HypsographyAdaption` type. A spectral second-order diffusion 
+elevation to the `HypsographyAdaption` type. A spectral second-order diffusion
 operator is applied with forward-Euler updates to generate
 profiles for each new iteration. Steps to generate smoothed terrain (
-represented as a ClimaCore Field) are as follows: 
+represented as a ClimaCore Field) are as follows:
 - Compute discrete elevation profile f
 - Compute diffuse_surface_elevation!(f, κ, iter). f is mutated.
 - Define `Hypsography.LinearAdaption(f)`
 - Define `ExtrudedFiniteDifferenceSpace` with new surface elevation.
-Default diffusion parameters are appropriate for spherical arrangements. 
+Default diffusion parameters are appropriate for spherical arrangements.
 For `zmax-zsfc` == 𝒪(10^4), κ == 𝒪(10^8), dt == 𝒪(10⁻¹).
 """
 function diffuse_surface_elevation!(
diff --git a/src/Remapping/Remapping.jl b/src/Remapping/Remapping.jl
index 08d8765de8..77ff1060e7 100644
--- a/src/Remapping/Remapping.jl
+++ b/src/Remapping/Remapping.jl
@@ -4,7 +4,13 @@ using LinearAlgebra, StaticArrays
 
 import ClimaComms
 import ..DataLayouts,
-    ..Geometry, ..Spaces, ..Topologies, ..Meshes, ..Operators, ..Fields
+    ..Geometry,
+    ..Spaces,
+    ..Topologies,
+    ..Meshes,
+    ..Operators,
+    ..Fields,
+    ..Hypsography
 
 using ..RecursiveApply
 
diff --git a/src/Remapping/distributed_remapping.jl b/src/Remapping/distributed_remapping.jl
index 3099169b1a..5818453d41 100644
--- a/src/Remapping/distributed_remapping.jl
+++ b/src/Remapping/distributed_remapping.jl
@@ -124,7 +124,6 @@ The `Remapper` is designed to not be tied to any particular `Field`. You can use
 `Remapper` for any `Field` as long as they are all defined on the same `topology`.
 
 `Remapper` is the main argument to the `interpolate` function.
-
 """
 function Remapper(target_hcoords, target_zcoords, space)
 
@@ -192,10 +191,19 @@ end
 
 
 """
-   interpolate(remapper, field)
+   interpolate(remapper, field; physical_z = false)
 
 Interpolate the given `field` as prescribed by `remapper`.
 
+
+Keyword arguments
+==================
+
+- `physical_z`: When `true`, interpolate to the physical z coordinates, taking into account
+                hypsography. If `false` (the default), interpolation will be based on the
+                reference z coordinate, i.e. will correspond to constant model levels.
+                `NaN`s are returned for values that are below the surface.
+
 Example
 ========
 
@@ -218,7 +226,11 @@ int2 = interpolate(remapper, field2)
 
 ```
 """
-function interpolate(remapper::Remapper, field::T) where {T <: Fields.Field}
+function interpolate(
+    remapper::Remapper,
+    field::T;
+    physical_z = false,
+) where {T <: Fields.Field}
 
     axes(field) == remapper.space ||
         error("Field is defined on a different space than remapper")
@@ -257,7 +269,8 @@ function interpolate(remapper::Remapper, field::T) where {T <: Fields.Field}
                     field,
                     remapper.target_zcoords,
                     weights,
-                    gidx,
+                    gidx;
+                    physical_z,
                 ) for (gidx, weights) in
                 zip(remapper.local_indices, remapper.interpolation_coeffs)
             )'
diff --git a/src/Remapping/interpolate_array.jl b/src/Remapping/interpolate_array.jl
index 06fde5df90..6a2467b30f 100644
--- a/src/Remapping/interpolate_array.jl
+++ b/src/Remapping/interpolate_array.jl
@@ -32,6 +32,24 @@ function interpolate_slab(
     return x
 end
 
+"""
+    interpolate_slab_level(
+                           field::Fields.Field,
+                           h::Integer,
+                           Is::Tuple,
+                           zcoord,
+                           )
+
+Vertically interpolate the given `field` on `zcoord`.
+
+The field is linearly interpolated across two neighboring vertical elements.
+
+For centered-valued fields, if `zcoord` is in the top (bottom) half of a top (bottom)
+element in a column, no interpolation is performed and the value at the cell center is
+returned. Effectively, this means that the interpolation is first-order accurate across the
+column, but zeroth-order accurate close to the boundaries.
+
+"""
 function interpolate_slab_level(
     field::Fields.Field,
     h::Integer,
@@ -197,12 +215,56 @@ and global index in the topology.
 
 The coefficients `weights` are computed with `Spaces.Quadratures.interpolation_matrix`.
 See also `interpolate_array`.
+
+Keyword arguments
+==================
+
+- `physical_z`: When `true`, the given `zpts` are interpreted as "from mean sea
+                            level" (instead of "from surface") and hypsography is
+                            interpolated. `NaN`s are returned for values that are below the
+                            surface.
+
 """
 function interpolate_column(
     field::Fields.ExtrudedFiniteDifferenceField,
     zpts,
     weights,
-    gidx,
+    gidx;
+    physical_z = false,
 )
-    return [interpolate_slab_level(field, gidx, weights, z) for z in zpts]
+    space = axes(field)
+    FT = Spaces.undertype(space)
+
+    # When we don't have hypsography, there is no notion of "interpolating hypsography". In
+    # this case, the reference vertical points coincide with the physical ones. Setting
+    # physical_z = false ensures that zpts_ref = zpts
+    if space.hypsography isa Spaces.Flat
+        physical_z = false
+    end
+
+    # If we physical_z, we have to move the z coordinates from physical to
+    # reference ones.
+    if physical_z
+        # We are hardcoding the transformation from Hypsography.LinearAdaption
+        space.hypsography isa Hypsography.LinearAdaption ||
+            error("Cannot interpolate $(space.hypsography) hypsography")
+
+        z_surface = interpolate_slab(
+            space.hypsography.surface,
+            Fields.SlabIndex(nothing, gidx),
+            weights,
+        )
+        z_top = Spaces.vertical_topology(space).mesh.domain.coord_max.z
+        zpts_ref = [
+            Geometry.ZPoint((z.z - z_surface) / (1 - z_surface / z_top)) for
+            z in zpts
+        ]
+    else
+        zpts_ref = zpts
+    end
+
+    return [
+        z.z >= 0 ? interpolate_slab_level(field, gidx, weights, z) : FT(NaN) for
+        z in zpts_ref
+    ]
 end
diff --git a/test/Remapping/distributed_remapping.jl b/test/Remapping/distributed_remapping.jl
index f666a3e62b..451bf756dc 100644
--- a/test/Remapping/distributed_remapping.jl
+++ b/test/Remapping/distributed_remapping.jl
@@ -3,7 +3,15 @@ using Test
 using IntervalSets
 
 import ClimaCore:
-    Domains, Fields, Geometry, Meshes, Operators, Spaces, Topologies, Remapping
+    Domains,
+    Fields,
+    Geometry,
+    Meshes,
+    Operators,
+    Spaces,
+    Topologies,
+    Remapping,
+    Hypsography
 using ClimaComms
 const context = ClimaComms.context()
 const pid, nprocs = ClimaComms.init(context)
@@ -68,6 +76,42 @@ end
         @test interp_z[:, end] ≈
               [1000.0 * (29 / 30 + 30 / 30) / 2 for x in xpts]
     end
+
+    # With topography
+    z_surface = cosd.(Fields.coordinate_field(horzspace).x) .+ 20
+    hv_center_space = Spaces.ExtrudedFiniteDifferenceSpace(
+        horzspace,
+        vert_center_space,
+        Hypsography.LinearAdaption(z_surface),
+    )
+
+    coords = Fields.coordinate_field(hv_center_space)
+
+    remapper = Remapping.Remapper(hcoords, zcoords, hv_center_space;)
+
+    function ref_z(x, z)
+        z_surface = cosd(x) + 20
+        z_top = vertdomain.coord_max.z
+        return (z - z_surface) / (1 - z_surface / z_top)
+    end
+
+    val_or_NaN(x, z, val) = ref_z(x, z) >= 0 ? val : NaN
+
+    interp_x = Remapping.interpolate(remapper, coords.x; physical_z = true)
+    if ClimaComms.iamroot(context)
+        @test interp_x ≈ [val_or_NaN(x, z, x) for x in xpts, z in zpts] nans = true
+    end
+    interp_z = Remapping.interpolate(remapper, coords.z; physical_z = true)
+    expected_z = [val_or_NaN(x, z, z) for x in xpts, z in zpts]
+    if ClimaComms.iamroot(context)
+        @test interp_z[:, 2:(end - 1)] ≈ expected_z[:, 2:(end - 1)] nans = true
+
+        z_bottom = 1000.0 * (0 / 30 + 1 / 30) / 2
+        @test interp_z[:, 1] ≈ [val_or_NaN(x, z_bottom, z_bottom) for x in xpts] nans = true
+
+        z_top = 1000.0 * (29 / 30 + 30 / 30) / 2
+        @test interp_z[:, end] ≈ [val_or_NaN(x, z_top, z_top) for x in xpts] rtol = 0.02 nans = true
+    end
 end
 
 @testset "3D box" begin
@@ -147,6 +191,60 @@ end
     if ClimaComms.iamroot(context)
         @test interp_y ≈ [y for x in xpts, y in ypts]
     end
+
+
+    # With topography
+    z_surface =
+        cosd.(Fields.coordinate_field(horzspace).x) .+
+        cosd.(Fields.coordinate_field(horzspace).y) .+ 20
+    hv_center_space = Spaces.ExtrudedFiniteDifferenceSpace(
+        horzspace,
+        vert_center_space,
+        Hypsography.LinearAdaption(z_surface),
+    )
+
+    coords = Fields.coordinate_field(hv_center_space)
+
+    remapper = Remapping.Remapper(hcoords, zcoords, hv_center_space;)
+
+    function ref_z(x, y, z)
+        z_surface = cosd(x) + cosd(y) + 20
+        z_top = vertdomain.coord_max.z
+        return (z - z_surface) / (1 - z_surface / z_top)
+    end
+
+    val_or_NaN(x, y, z, val) = ref_z(x, y, z) >= 0 ? val : NaN
+
+    interp_x = Remapping.interpolate(remapper, coords.x; physical_z = true)
+    if ClimaComms.iamroot(context)
+        @test interp_x ≈
+              [val_or_NaN(x, y, z, x) for x in xpts, y in ypts, z in zpts] nans =
+            true
+    end
+
+    interp_y = Remapping.interpolate(remapper, coords.y; physical_z = true)
+    if ClimaComms.iamroot(context)
+        @test interp_y ≈
+              [val_or_NaN(x, y, z, y) for x in xpts, y in ypts, z in zpts] nans =
+            true
+    end
+
+    interp_z = Remapping.interpolate(remapper, coords.z; physical_z = true)
+    expected_z = [val_or_NaN(x, y, z, z) for x in xpts, y in ypts, z in zpts]
+    if ClimaComms.iamroot(context)
+        @test interp_z[:, :, 2:(end - 1)] ≈ expected_z[:, :, 2:(end - 1)] nans =
+            true
+
+        z_bottom = 1000.0 * (0 / 30 + 1 / 30) / 2
+        @test interp_z[:, :, 1] ≈
+              [val_or_NaN(x, y, z_bottom, z_bottom) for x in xpts, y in ypts] nans =
+            true
+
+        z_top = 1000.0 * (29 / 30 + 30 / 30) / 2
+        @test interp_z[:, :, end] ≈
+              [val_or_NaN(x, y, z_top, z_top) for x in xpts, y in ypts] rtol =
+            0.02 nans = true
+    end
 end
 
 
@@ -229,6 +327,59 @@ end
     if ClimaComms.iamroot(context)
         @test interp_y ≈ [y for x in xpts, y in ypts]
     end
+
+    # With topography
+    z_surface =
+        cosd.(Fields.coordinate_field(horzspace).x) .+
+        cosd.(Fields.coordinate_field(horzspace).y) .+ 20
+    hv_center_space = Spaces.ExtrudedFiniteDifferenceSpace(
+        horzspace,
+        vert_center_space,
+        Hypsography.LinearAdaption(z_surface),
+    )
+
+    coords = Fields.coordinate_field(hv_center_space)
+
+    remapper = Remapping.Remapper(hcoords, zcoords, hv_center_space;)
+
+    function ref_z(x, y, z)
+        z_surface = cosd(x) + cosd(y) + 20
+        z_top = vertdomain.coord_max.z
+        return (z - z_surface) / (1 - z_surface / z_top)
+    end
+
+    val_or_NaN(x, y, z, val) = ref_z(x, y, z) >= 0 ? val : NaN
+
+    interp_x = Remapping.interpolate(remapper, coords.x; physical_z = true)
+    if ClimaComms.iamroot(context)
+        @test interp_x ≈
+              [val_or_NaN(x, y, z, x) for x in xpts, y in ypts, z in zpts] nans =
+            true
+    end
+
+    interp_y = Remapping.interpolate(remapper, coords.y; physical_z = true)
+    if ClimaComms.iamroot(context)
+        @test interp_y ≈
+              [val_or_NaN(x, y, z, y) for x in xpts, y in ypts, z in zpts] nans =
+            true
+    end
+
+    interp_z = Remapping.interpolate(remapper, coords.z; physical_z = true)
+    expected_z = [val_or_NaN(x, y, z, z) for x in xpts, y in ypts, z in zpts]
+    if ClimaComms.iamroot(context)
+        @test interp_z[:, :, 2:(end - 1)] ≈ expected_z[:, :, 2:(end - 1)] nans =
+            true
+
+        z_bottom = 1000.0 * (0 / 30 + 1 / 30) / 2
+        @test interp_z[:, :, 1] ≈
+              [val_or_NaN(x, y, z_bottom, z_bottom) for x in xpts, y in ypts] nans =
+            true
+
+        z_top = 1000.0 * (29 / 30 + 30 / 30) / 2
+        @test interp_z[:, :, end] ≈
+              [val_or_NaN(x, y, z_top, z_top) for x in xpts, y in ypts] rtol =
+            0.02 nans = true
+    end
 end
 
 @testset "3D sphere" begin
@@ -305,4 +456,69 @@ end
     if ClimaComms.iamroot(context)
         @test interp_sin_lat ≈ [sind(y) for x in longpts, y in latpts] rtol = 0.01
     end
+
+    # With topography
+    #
+    # When we have topography, we should allow for larger error because elements close to
+    # the surface are not reconstructed well
+    z_surface =
+        cosd.(Fields.coordinate_field(horzspace).lat) .+
+        cosd.(Fields.coordinate_field(horzspace).long) .+ 20
+    hv_center_space = Spaces.ExtrudedFiniteDifferenceSpace(
+        horzspace,
+        vert_center_space,
+        Hypsography.LinearAdaption(z_surface),
+    )
+
+    coords = Fields.coordinate_field(hv_center_space)
+
+    remapper = Remapping.Remapper(hcoords, zcoords, hv_center_space)
+
+    function ref_z(x, y, z)
+        z_surface = cosd(x) + cosd(y) + 20
+        z_top = vertdomain.coord_max.z
+        return (z - z_surface) / (1 - z_surface / z_top)
+    end
+
+    val_or_NaN(x, y, z, val) = ref_z(x, y, z) >= 0 ? val : NaN
+
+    interp_sin_long =
+        Remapping.interpolate(remapper, sind.(coords.long); physical_z = true)
+
+    if ClimaComms.iamroot(context)
+        @test interp_sin_long ≈ [
+            val_or_NaN(x, y, z, sind(x)) for x in longpts, y in latpts,
+            z in zpts
+        ] rtol = 0.06 nans = true
+    end
+
+    interp_sin_lat =
+        Remapping.interpolate(remapper, sind.(coords.lat); physical_z = true)
+    if ClimaComms.iamroot(context)
+        # We have to add atol because we have an element that is identically zero
+        @test interp_sin_lat ≈ [
+            val_or_NaN(x, y, z, sind(y)) for x in longpts, y in latpts,
+            z in zpts
+        ] rtol = 0.01 nans = true atol = 1e-15
+    end
+
+    interp_z = Remapping.interpolate(remapper, coords.z, physical_z = true)
+
+    FT = Spaces.undertype(hv_center_space)
+
+    expected_z = [z for x in longpts, y in latpts, z in zpts]
+
+    if ClimaComms.iamroot(context)
+        @test interp_z[:, :, 2:(end - 1)] ≈ expected_z[:, :, 2:(end - 1)]
+
+        z_bottom = 1000.0 * (0 / 30 + 1 / 30) / 2
+        @test interp_z[:, :, 1] ≈ [
+            val_or_NaN(x, y, z_bottom, z_bottom) for x in longpts, y in latpts
+        ] nans = true
+
+        z_top = 1000.0 * (29 / 30 + 30 / 30) / 2
+        @test interp_z[:, :, end] ≈
+              [val_or_NaN(x, y, z_top, z_top) for x in longpts, y in latpts] nans =
+            true rtol = 0.02
+    end
 end