diff --git a/amr-wind/equation_systems/icns/source_terms/DragForcing.cpp b/amr-wind/equation_systems/icns/source_terms/DragForcing.cpp
index 36655d04ec..2a6c6654c1 100644
--- a/amr-wind/equation_systems/icns/source_terms/DragForcing.cpp
+++ b/amr-wind/equation_systems/icns/source_terms/DragForcing.cpp
@@ -4,6 +4,7 @@
 #include "AMReX_Gpu.H"
 #include "AMReX_Random.H"
 #include "amr-wind/wind_energy/ABL.H"
+#include "amr-wind/utilities/linear_interpolation.H"
 
 namespace amr_wind::pde::icns {
 
@@ -115,23 +116,21 @@ void DragForcing::operator()(
         yi_end = sponge_north * std::max(yi_end, 0.0);
         ystart_damping = sponge_strength * yi_start * yi_start;
         yend_damping = sponge_strength * yi_end * yi_end;
-        amrex::Real spongeVelX = 0.0;
-        amrex::Real spongeVelY = 0.0;
-        amrex::Real spongeVelZ = 0.0;
-        amrex::Real residual = 1000;
-        amrex::Real height_error = 0.0;
-        for (unsigned ii = 0; ii < vsize; ++ii) {
-            height_error = std::abs(z - device_vel_ht[ii]);
-            if (height_error < residual) {
-                residual = height_error;
-                const unsigned ix = 3 * ii;
-                const unsigned iy = 3 * ii + 1;
-                const unsigned iz = 3 * ii + 2;
-                spongeVelX = device_vel_vals[ix];
-                spongeVelY = device_vel_vals[iy];
-                spongeVelZ = device_vel_vals[iz];
-            }
-        }
+
+        const auto idx =
+            interp::bisection_search(device_vel_ht, device_vel_ht + vsize, z);
+        const amrex::Real spongeVelX =
+            (vsize > 0) ? interp::linear_impl(
+                              device_vel_ht, device_vel_vals, z, idx, 3, 0)
+                        : 0.0;
+        const amrex::Real spongeVelY =
+            (vsize > 0) ? interp::linear_impl(
+                              device_vel_ht, device_vel_vals, z, idx, 3, 1)
+                        : 0.0;
+        const amrex::Real spongeVelZ =
+            (vsize > 0) ? interp::linear_impl(
+                              device_vel_ht, device_vel_vals, z, idx, 3, 2)
+                        : 0.0;
         amrex::Real Dxz = 0.0;
         amrex::Real Dyz = 0.0;
         amrex::Real bc_forcing_x = 0;
diff --git a/amr-wind/utilities/linear_interpolation.H b/amr-wind/utilities/linear_interpolation.H
index 40c5994a83..926dcdf5e0 100644
--- a/amr-wind/utilities/linear_interpolation.H
+++ b/amr-wind/utilities/linear_interpolation.H
@@ -61,6 +61,23 @@ bisection_search(const It begin, const It end, const T& x)
     return idx;
 }
 
+template <typename It, typename T>
+AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Index
+nearest_search(const It begin, const It end, const T& x)
+{
+
+    auto idx = bisection_search(begin, end, x);
+    if ((idx.lim == Limits::LOWLIM) || (idx.lim == Limits::UPLIM)) {
+        return idx;
+    }
+
+    if ((x - begin[idx.idx]) > (begin[idx.idx + 1] - x)) {
+        idx.idx += 1;
+    }
+
+    return idx;
+}
+
 template <typename It, typename T>
 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Index
 find_index(const It begin, const It end, const T& x, const int hint = 1)
@@ -82,36 +99,63 @@ find_index(const It begin, const It end, const T& x, const int hint = 1)
 template <typename C1, typename C2>
 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
     typename std::iterator_traits<C2>::value_type
-    linear(
+    linear_impl(
         const C1 xbegin,
-        const C1 xend,
         const C2 yinp,
-        const typename std::iterator_traits<C1>::value_type& xout)
+        const typename std::iterator_traits<C1>::value_type& xout,
+        const Index& idx,
+        const int ncomp = 1,
+        const int comp = 0)
 {
     using DType1 = typename std::iterator_traits<C1>::value_type;
-    const auto idx = bisection_search(xbegin, xend, xout);
 
     if ((idx.lim == Limits::LOWLIM) || (idx.lim == Limits::UPLIM)) {
-        return yinp[idx.idx];
+        return yinp[ncomp * idx.idx + comp];
     }
     static constexpr DType1 eps = 1.0e-8;
     const int j = idx.idx;
     const auto denom = (xbegin[j + 1] - xbegin[j]);
     const auto facR = (denom > eps) ? ((xout - xbegin[j]) / denom) : 1.0;
     const auto facL = static_cast<DType1>(1.0) - facR;
-    return facL * yinp[j] + facR * yinp[j + 1];
+    return facL * yinp[ncomp * j + comp] + facR * yinp[ncomp * (j + 1) + comp];
 }
 
 template <typename C1, typename C2>
-inline typename C2::value_type
-linear(const C1& xinp, const C2& yinp, const typename C1::value_type& xout)
+AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
+    typename std::iterator_traits<C2>::value_type
+    linear(
+        const C1 xbegin,
+        const C1 xend,
+        const C2 yinp,
+        const typename std::iterator_traits<C1>::value_type& xout,
+        const int ncomp = 1,
+        const int comp = 0)
+{
+    const auto idx = bisection_search(xbegin, xend, xout);
+    return linear_impl(xbegin, yinp, xout, idx, ncomp, comp);
+}
+
+template <typename C1, typename C2>
+inline typename C2::value_type linear(
+    const C1& xinp,
+    const C2& yinp,
+    const typename C1::value_type& xout,
+    const int ncomp = 1,
+    const int comp = 0)
 {
-    return linear(xinp.data(), (xinp.data() + xinp.size()), yinp.data(), xout);
+    return linear(
+        xinp.data(), (xinp.data() + xinp.size()), yinp.data(), xout, ncomp,
+        comp);
 }
 
 template <typename C1, typename C2>
-inline void
-linear_monotonic(const C1& xinp, const C2& yinp, const C1& xout, C2& yout)
+inline void linear_monotonic(
+    const C1& xinp,
+    const C2& yinp,
+    const C1& xout,
+    C2& yout,
+    const int ncomp = 1,
+    const int comp = 0)
 {
     static constexpr typename C1::value_type eps = 1.0e-8;
     AMREX_ASSERT(xinp.size() == yinp.size());
@@ -125,20 +169,28 @@ linear_monotonic(const C1& xinp, const C2& yinp, const C1& xout, C2& yout)
             find_index(xinp.data(), (xinp.data() + xinp.size()), x, hint);
 
         if ((idx.lim == Limits::LOWLIM) || (idx.lim == Limits::UPLIM)) {
-            yout[i] = yinp[idx.idx];
+            yout[i] = yinp[ncomp * idx.idx + comp];
         } else if (idx.lim == Limits::VALID) {
             int j = idx.idx;
             const auto denom = (xinp[j + 1] - xinp[j]);
             const auto facR = (denom > eps) ? ((x - xinp[j]) / denom) : 1.0;
             const auto facL = static_cast<typename C1::value_type>(1.0) - facR;
-            yout[i] = facL * yinp[j] + facR * yinp[j + 1];
+            yout[i] = facL * yinp[ncomp * j + comp] +
+                      facR * yinp[ncomp * (j + 1) + comp];
+            ;
         }
         hint = idx.idx + 1;
     }
 }
 
 template <typename C1, typename C2>
-inline void linear(const C1& xinp, const C2& yinp, const C1& xout, C2& yout)
+inline void linear(
+    const C1& xinp,
+    const C2& yinp,
+    const C1& xout,
+    C2& yout,
+    const int ncomp = 1,
+    const int comp = 0)
 {
     AMREX_ASSERT(
         static_cast<amrex::Long>(xinp.size()) ==
@@ -149,7 +201,7 @@ inline void linear(const C1& xinp, const C2& yinp, const C1& xout, C2& yout)
 
     int npts = xout.size();
     for (int i = 0; i < npts; ++i) {
-        yout[i] = linear(xinp, yinp, xout[i]);
+        yout[i] = linear(xinp, yinp, xout[i], ncomp, comp);
     }
 }
 
diff --git a/unit_tests/utilities/test_linear_interpolation.cpp b/unit_tests/utilities/test_linear_interpolation.cpp
index a18cde5eb1..3ef880e1be 100644
--- a/unit_tests/utilities/test_linear_interpolation.cpp
+++ b/unit_tests/utilities/test_linear_interpolation.cpp
@@ -60,6 +60,46 @@ TEST(LinearInterpolation, bisection_search)
     }
 }
 
+TEST(LinearInterpolation, nearest_search)
+{
+    namespace interp = amr_wind::interp;
+    std::vector<amrex::Real> xvec(10);
+    std::iota(xvec.begin(), xvec.end(), 0.0);
+
+    const auto* start = xvec.data();
+    const auto* end = xvec.data() + xvec.size();
+    {
+        const auto idx = interp::nearest_search(start, end, 5.0);
+        EXPECT_EQ(idx.idx, 5);
+        EXPECT_EQ(idx.lim, interp::Limits::VALID);
+    }
+    {
+        const auto idx = interp::nearest_search(start, end, 1.1);
+        EXPECT_EQ(idx.idx, 1);
+        EXPECT_EQ(idx.lim, interp::Limits::VALID);
+    }
+    {
+        const auto idx = interp::nearest_search(start, end, 1.6);
+        EXPECT_EQ(idx.idx, 2);
+        EXPECT_EQ(idx.lim, interp::Limits::VALID);
+    }
+    {
+        const auto idx = interp::nearest_search(start, end, 3.5);
+        EXPECT_EQ(idx.idx, 3);
+        EXPECT_EQ(idx.lim, interp::Limits::VALID);
+    }
+    {
+        const auto idx = interp::nearest_search(start, end, -1.0);
+        EXPECT_EQ(idx.idx, 0);
+        EXPECT_EQ(idx.lim, interp::Limits::LOWLIM);
+    }
+    {
+        const auto idx = interp::nearest_search(start, end, 9.1);
+        EXPECT_EQ(idx.idx, 9);
+        EXPECT_EQ(idx.lim, interp::Limits::UPLIM);
+    }
+}
+
 TEST(LinearInterpolation, find_index)
 {
     namespace interp = amr_wind::interp;
@@ -108,6 +148,31 @@ TEST(LinearInterpolation, lin_interp_single)
     }
 }
 
+TEST(LinearInterpolation, lin_interp_single_multicomponent)
+{
+    namespace interp = amr_wind::interp;
+
+    const int ncomp = 3;
+    std::vector<amrex::Real> xvec(10), yvec(ncomp * 10);
+    std::iota(xvec.begin(), xvec.end(), 0.0);
+    const amrex::Vector<amrex::Real> mult_facs = {
+        2.0 + 10.0 * amrex::Random(), 2.0 + 10.0 * amrex::Random(),
+        2.0 + 10.0 * amrex::Random()};
+    for (int i = 0; i < static_cast<int>(xvec.size()); i++) {
+        for (int n = 0; n < ncomp; n++) {
+            yvec[ncomp * i + n] = mult_facs[n] * xvec[i];
+        }
+    }
+
+    std::vector<amrex::Real> xtest{2.5, 4.5, 6.3, 8.8};
+    for (const auto& x : xtest) {
+        for (int n = 0; n < ncomp; n++) {
+            const auto y = interp::linear(xvec, yvec, x, 3, n);
+            EXPECT_NEAR(y, mult_facs[n] * x, 1.0e-12);
+        }
+    }
+}
+
 TEST(LinearInterpolation, lin_interp)
 {
     namespace interp = amr_wind::interp;