Repository: OceanParcels/Parcels
Pull request #1895 (2025-02-27 17:11:39)

Moving the time interpolation down to within _Interpolator functions
Parcels previously would run an interpolation on time index `ti` and one on `ti+1`, and then at the end linearly interpolate the results.

In this PR, we change that to do the time-interpolation closer to the `InterpolationContext`. We need to discuss whether this is the right place now, or whether we want to push it even into the interpolators themselves. The latter would open possible support also non-linear time-interpolation

Hence, this is more a proof of concept than a full-fledged PR to be merged...


- [x] Chose the correct base branch (`main` for v3 changes, `v4-dev` for v4 changes)



Patch:
diff --git a/parcels/_interpolation.py b/parcels/_interpolation.py
index a533fd9d..88a2864d 100644
--- a/parcels/_interpolation.py
+++ b/parcels/_interpolation.py
@@ -4,6 +4,7 @@ from dataclasses import dataclass
 import numpy as np
 
 from parcels._typing import GridIndexingType
+from parcels.tools._helpers import should_calculate_next_ti
 
 
 @dataclass
@@ -14,6 +15,8 @@ class InterpolationContext2D:
     ----------
     data: np.ndarray
         field data of shape (time, y, x)
+    tau: float
+        time interpolation coordinate in unit length
     eta: float
         y-direction interpolation coordinate in unit cube (between 0 and 1)
     xsi: float
@@ -28,6 +31,7 @@ class InterpolationContext2D:
     """
 
     data: np.ndarray
+    tau: float
     eta: float
     xsi: float
     ti: int
@@ -45,6 +49,8 @@ class InterpolationContext3D:
         field data of shape (time, z, y, x). This needs to be complete in the vertical
         direction as some interpolation methods need to know whether they are at the
         surface or bottom.
+    tau: float
+        time interpolation coordinate in unit length
     zeta: float
         vertical interpolation coordinate in unit cube
     eta: float
@@ -65,6 +71,7 @@ class InterpolationContext3D:
     """
 
     data: np.ndarray
+    tau: float
     zeta: float
     eta: float
     xsi: float
@@ -110,7 +117,11 @@ def register_3d_interpolator(name: str):
 def _nearest_2d(ctx: InterpolationContext2D) -> float:
     xii = ctx.xi if ctx.xsi <= 0.5 else ctx.xi + 1
     yii = ctx.yi if ctx.eta <= 0.5 else ctx.yi + 1
-    return ctx.data[ctx.ti, yii, xii]
+    ft0 = ctx.data[ctx.ti, yii, xii]
+    if not should_calculate_next_ti(ctx.ti, ctx.tau, ctx.data.shape[0]):
+        return ft0
+    ft1 = ctx.data[ctx.ti + 1, yii, xii]
+    return (1 - ctx.tau) * ft0 + ctx.tau * ft1
 
 
 def _interp_on_unit_square(*, eta: float, xsi: float, data: np.ndarray, yi: int, xi: int) -> float:
@@ -128,7 +139,11 @@ def _interp_on_unit_square(*, eta: float, xsi: float, data: np.ndarray, yi: int,
 @register_2d_interpolator("partialslip")
 @register_2d_interpolator("freeslip")
 def _linear_2d(ctx: InterpolationContext2D) -> float:
-    return _interp_on_unit_square(eta=ctx.eta, xsi=ctx.xsi, data=ctx.data[ctx.ti, :, :], yi=ctx.yi, xi=ctx.xi)
+    ft0 = _interp_on_unit_square(eta=ctx.eta, xsi=ctx.xsi, data=ctx.data[ctx.ti, :, :], yi=ctx.yi, xi=ctx.xi)
+    if not should_calculate_next_ti(ctx.ti, ctx.tau, ctx.data.shape[0]):
+        return ft0
+    ft1 = _interp_on_unit_square(eta=ctx.eta, xsi=ctx.xsi, data=ctx.data[ctx.ti + 1, :, :], yi=ctx.yi, xi=ctx.xi)
+    return (1 - ctx.tau) * ft0 + ctx.tau * ft1
 
 
 @register_2d_interpolator("linear_invdist_land_tracer")
@@ -142,6 +157,13 @@ def _linear_invdist_land_tracer_2d(ctx: InterpolationContext2D) -> float:
     land = np.isclose(data[ti, yi : yi + 2, xi : xi + 2], 0.0)
     nb_land = np.sum(land)
 
+    def _get_data_temporalinterp(*, ti, yi, xi):
+        dt0 = data[ti, yi, xi]
+        if not should_calculate_next_ti(ctx.ti, ctx.tau, ctx.data.shape[0]):
+            return dt0
+        dt1 = data[ti + 1, yi, xi]
+        return (1 - ctx.tau) * dt0 + ctx.tau * dt1
+
     if nb_land == 4:
         return 0
     elif nb_land > 0:
@@ -154,9 +176,9 @@ def _linear_invdist_land_tracer_2d(ctx: InterpolationContext2D) -> float:
                     if land[j][i] == 1:  # index search led us directly onto land
                         return 0
                     else:
-                        return data[ti, yi + j, xi + i]
+                        return _get_data_temporalinterp(ti=ti, yi=yi + j, xi=xi + i)
                 elif land[j][i] == 0:
-                    val += data[ti, yi + j, xi + i] / distance
+                    val += _get_data_temporalinterp(ti=ti, yi=yi + j, xi=xi + i) / distance
                     w_sum += 1 / distance
         return val / w_sum
     else:
@@ -166,7 +188,11 @@ def _linear_invdist_land_tracer_2d(ctx: InterpolationContext2D) -> float:
 @register_2d_interpolator("cgrid_tracer")
 @register_2d_interpolator("bgrid_tracer")
 def _tracer_2d(ctx: InterpolationContext2D) -> float:
-    return ctx.data[ctx.ti, ctx.yi + 1, ctx.xi + 1]
+    ft0 = ctx.data[ctx.ti, ctx.yi + 1, ctx.xi + 1]
+    if not should_calculate_next_ti(ctx.ti, ctx.tau, ctx.data.shape[0]):
+        return ft0
+    ft1 = ctx.data[ctx.ti + 1, ctx.yi + 1, ctx.xi + 1]
+    return (1 - ctx.tau) * ft0 + ctx.tau * ft1
 
 
 @register_3d_interpolator("nearest")
@@ -174,25 +200,52 @@ def _nearest_3d(ctx: InterpolationContext3D) -> float:
     xii = ctx.xi if ctx.xsi <= 0.5 else ctx.xi + 1
     yii = ctx.yi if ctx.eta <= 0.5 else ctx.yi + 1
     zii = ctx.zi if ctx.zeta <= 0.5 else ctx.zi + 1
-    return ctx.data[ctx.ti, zii, yii, xii]
+    ft0 = ctx.data[ctx.ti, zii, yii, xii]
+    if not should_calculate_next_ti(ctx.ti, ctx.tau, ctx.data.shape[0]):
+        return ft0
+    ft1 = ctx.data[ctx.ti + 1, zii, yii, xii]
+    return (1 - ctx.tau) * ft0 + ctx.tau * ft1
+
+
+def _get_cgrid_depth_point(*, zeta: float, data: np.ndarray, zi: int, yi: int, xi: int) -> float:
+    f0 = data[zi, yi, xi]
+    f1 = data[zi + 1, yi, xi]
+    return (1 - zeta) * f0 + zeta * f1
 
 
 @register_3d_interpolator("cgrid_velocity")
-def _cgrid_velocity_3d(ctx: InterpolationContext3D) -> float:
+def _cgrid_W_velocity_3d(ctx: InterpolationContext3D) -> float:
     # evaluating W velocity in c_grid
     if ctx.gridindexingtype == "nemo":
-        f0 = ctx.data[ctx.ti, ctx.zi, ctx.yi + 1, ctx.xi + 1]
-        f1 = ctx.data[ctx.ti, ctx.zi + 1, ctx.yi + 1, ctx.xi + 1]
+        ft0 = _get_cgrid_depth_point(
+            zeta=ctx.zeta, data=ctx.data[ctx.ti, :, :, :], zi=ctx.zi, yi=ctx.yi + 1, xi=ctx.xi + 1
+        )
     elif ctx.gridindexingtype in ["mitgcm", "croco"]:
-        f0 = ctx.data[ctx.ti, ctx.zi, ctx.yi, ctx.xi]
-        f1 = ctx.data[ctx.ti, ctx.zi + 1, ctx.yi, ctx.xi]
-    return (1 - ctx.zeta) * f0 + ctx.zeta * f1
+        ft0 = _get_cgrid_depth_point(zeta=ctx.zeta, data=ctx.data[ctx.ti, :, :, :], zi=ctx.zi, yi=ctx.yi, xi=ctx.xi)
+    if not should_calculate_next_ti(ctx.ti, ctx.tau, ctx.data.shape[0]):
+        return ft0
+
+    if ctx.gridindexingtype == "nemo":
+        ft1 = _get_cgrid_depth_point(
+            zeta=ctx.zeta, data=ctx.data[ctx.ti + 1, :, :, :], zi=ctx.zi, yi=ctx.yi + 1, xi=ctx.xi + 1
+        )
+    elif ctx.gridindexingtype in ["mitgcm", "croco"]:
+        ft1 = _get_cgrid_depth_point(zeta=ctx.zeta, data=ctx.data[ctx.ti + 1, :, :, :], zi=ctx.zi, yi=ctx.yi, xi=ctx.xi)
+    return (1 - ctx.tau) * ft0 + ctx.tau * ft1
 
 
 @register_3d_interpolator("linear_invdist_land_tracer")
 def _linear_invdist_land_tracer_3d(ctx: InterpolationContext3D) -> float:
     land = np.isclose(ctx.data[ctx.ti, ctx.zi : ctx.zi + 2, ctx.yi : ctx.yi + 2, ctx.xi : ctx.xi + 2], 0.0)
     nb_land = np.sum(land)
+
+    def _get_data_temporalinterp(*, ti, zi, yi, xi):
+        dt0 = ctx.data[ti, zi, yi, xi]
+        if not should_calculate_next_ti(ctx.ti, ctx.tau, ctx.data.shape[0]):
+            return dt0
+        dt1 = data[ti + 1, zi, yi, xi]
+        return (1 - ctx.tau) * dt0 + ctx.tau * dt1
+
     if nb_land == 8:
         return 0
     elif nb_land > 0:
@@ -206,9 +259,11 @@ def _linear_invdist_land_tracer_3d(ctx: InterpolationContext3D) -> float:
                         if land[k][j][i] == 1:  # index search led us directly onto land
                             return 0
                         else:
-                            return ctx.data[ctx.ti, ctx.zi + k, ctx.yi + j, ctx.xi + i]
+                            return _get_data_temporalinterp(ti=ctx.ti, zi=ctx.zi + k, yi=ctx.yi + j, xi=ctx.xi + i)
                     elif land[k][j][i] == 0:
-                        val += ctx.data[ctx.ti, ctx.zi + k, ctx.yi + j, ctx.xi + i] / distance
+                        val += (
+                            _get_data_temporalinterp(ti=ctx.ti, zi=ctx.zi + k, yi=ctx.yi + j, xi=ctx.xi + i) / distance
+                        )
                         w_sum += 1 / distance
         return val / w_sum
     else:
@@ -253,9 +308,15 @@ def _z_layer_interp(
 def _linear_3d(ctx: InterpolationContext3D) -> float:
     zdim = ctx.data.shape[1]
     data_3d = ctx.data[ctx.ti, :, :, :]
-    f0, f1 = _get_3d_f0_f1(eta=ctx.eta, xsi=ctx.xsi, data=data_3d, zi=ctx.zi, yi=ctx.yi, xi=ctx.xi)
+    fz0, fz1 = _get_3d_f0_f1(eta=ctx.eta, xsi=ctx.xsi, data=data_3d, zi=ctx.zi, yi=ctx.yi, xi=ctx.xi)
+    if should_calculate_next_ti(ctx.ti, ctx.tau, ctx.data.shape[0]):
+        data_3d = ctx.data[ctx.ti + 1, :, :, :]
+        fz0_t1, fz1_t1 = _get_3d_f0_f1(eta=ctx.eta, xsi=ctx.xsi, data=data_3d, zi=ctx.zi, yi=ctx.yi, xi=ctx.xi)
+        fz0 = (1 - ctx.tau) * fz0 + ctx.tau * fz0_t1
+        if fz1_t1 is not None and fz1 is not None:
+            fz1 = (1 - ctx.tau) * fz1 + ctx.tau * fz1_t1
 
-    return _z_layer_interp(zeta=ctx.zeta, f0=f0, f1=f1, zi=ctx.zi, zdim=zdim, gridindexingtype=ctx.gridindexingtype)
+    return _z_layer_interp(zeta=ctx.zeta, f0=fz0, f1=fz1, zi=ctx.zi, zdim=zdim, gridindexingtype=ctx.gridindexingtype)
 
 
 @register_3d_interpolator("bgrid_velocity")
@@ -277,4 +338,8 @@ def _linear_3d_bgrid_w_velocity(ctx: InterpolationContext3D) -> float:
 @register_3d_interpolator("bgrid_tracer")
 @register_3d_interpolator("cgrid_tracer")
 def _tracer_3d(ctx: InterpolationContext3D) -> float:
-    return ctx.data[ctx.ti, ctx.zi, ctx.yi + 1, ctx.xi + 1]
+    ft0 = ctx.data[ctx.ti, ctx.zi, ctx.yi + 1, ctx.xi + 1]
+    if not should_calculate_next_ti(ctx.ti, ctx.tau, ctx.data.shape[0]):
+        return ft0
+    ft1 = ctx.data[ctx.ti + 1, ctx.zi, ctx.yi + 1, ctx.xi + 1]
+    return (1 - ctx.tau) * ft0 + ctx.tau * ft1
diff --git a/parcels/application_kernels/advection.py b/parcels/application_kernels/advection.py
index d35d3063..faa50be7 100644
--- a/parcels/application_kernels/advection.py
+++ b/parcels/application_kernels/advection.py
@@ -177,16 +177,14 @@ def AdvectionAnalytical(particle, fieldset, time):  # pragma: no cover
     direction = 1.0 if particle.dt > 0 else -1.0
     withW = True if "W" in [f.name for f in fieldset.get_fields()] else False
     withTime = True if len(fieldset.U.grid.time_full) > 1 else False
-    ti = fieldset.U._time_index(time)
+    tau, zeta, eta, xsi, ti, zi, yi, xi = fieldset.U._search_indices(
+        time, particle.depth, particle.lat, particle.lon, particle=particle
+    )
     ds_t = particle.dt
     if withTime:
-        tau = (time - fieldset.U.grid.time[ti]) / (fieldset.U.grid.time[ti + 1] - fieldset.U.grid.time[ti])
         time_i = np.linspace(0, fieldset.U.grid.time[ti + 1] - fieldset.U.grid.time[ti], I_s)
         ds_t = min(ds_t, time_i[np.where(time - fieldset.U.grid.time[ti] < time_i)[0][0]])
 
-    zeta, eta, xsi, zi, yi, xi = fieldset.U._search_indices(
-        time, particle.depth, particle.lat, particle.lon, ti, particle=particle
-    )
     if withW:
         if abs(xsi - 1) < tol:
             if fieldset.U.data[0, zi + 1, yi + 1, xi + 1] > 0:
diff --git a/parcels/field.py b/parcels/field.py
index e360fae2..09521431 100644
--- a/parcels/field.py
+++ b/parcels/field.py
@@ -25,7 +25,7 @@ from parcels._typing import (
     assert_valid_gridindexingtype,
     assert_valid_interp_method,
 )
-from parcels.tools._helpers import default_repr, field_repr
+from parcels.tools._helpers import default_repr, field_repr, should_calculate_next_ti
 from parcels.tools.converters import (
     TimeConverter,
     UnitConverter,
@@ -542,17 +542,21 @@ class Field:
 
         return data
 
-    def _search_indices(self, time, z, y, x, ti, particle=None, search2D=False):
+    def _search_indices(self, time, z, y, x, particle=None, search2D=False):
+        tau, ti = self._search_time_index(time)
+
         if self.grid._gtype in [GridType.RectilinearSGrid, GridType.RectilinearZGrid]:
-            return _search_indices_rectilinear(self, time, z, y, x, ti, particle=particle, search2D=search2D)
+            (zeta, eta, xsi, zi, yi, xi) = _search_indices_rectilinear(
+                self, time, z, y, x, ti, particle=particle, search2D=search2D
+            )
         else:
-            return _search_indices_curvilinear(self, time, z, y, x, ti, particle=particle, search2D=search2D)
+            (zeta, eta, xsi, zi, yi, xi) = _search_indices_curvilinear(
+                self, time, z, y, x, ti, particle=particle, search2D=search2D
+            )
+        return (tau, zeta, eta, xsi, ti, zi, yi, xi)
 
-    def _interpolator2D(self, time, z, y, x, ti, particle=None):
+    def _interpolator2D(self, time, z, y, x, particle=None):
         """Impelement 2D interpolation with coordinate transformations as seen in Delandmeter and Van Sebille (2019), 10.5194/gmd-12-3571-2019.."""
-        (_, eta, xsi, _, yi, xi) = self._search_indices(time, z, y, x, ti, particle=particle)
-        ctx = InterpolationContext2D(self.data, eta, xsi, ti, yi, xi)
-
         try:
             f = get_2d_interpolator_registry()[self.interp_method]
         except KeyError:
@@ -562,26 +566,31 @@ class Field:
                 )
             else:
                 raise RuntimeError(self.interp_method + " is not implemented for 2D grids")
+
+        (tau, _, eta, xsi, ti, _, yi, xi) = self._search_indices(time, z, y, x, particle=particle)
+
+        ctx = InterpolationContext2D(self.data, tau, eta, xsi, ti, yi, xi)
         return f(ctx)
 
-    def _interpolator3D(self, time, z, y, x, ti, particle=None):
+    def _interpolator3D(self, time, z, y, x, particle=None):
         """Impelement 3D interpolation with coordinate transformations as seen in Delandmeter and Van Sebille (2019), 10.5194/gmd-12-3571-2019.."""
-        (zeta, eta, xsi, zi, yi, xi) = self._search_indices(time, z, y, x, ti, particle=particle)
-        ctx = InterpolationContext3D(self.data, zeta, eta, xsi, ti, zi, yi, xi, self.gridindexingtype)
-
         try:
             f = get_3d_interpolator_registry()[self.interp_method]
         except KeyError:
             raise RuntimeError(self.interp_method + " is not implemented for 3D grids")
+
+        (tau, zeta, eta, xsi, ti, zi, yi, xi) = self._search_indices(time, z, y, x, particle=particle)
+
+        ctx = InterpolationContext3D(self.data, tau, zeta, eta, xsi, ti, zi, yi, xi, self.gridindexingtype)
         return f(ctx)
 
-    def _spatial_interpolation(self, time, z, y, x, ti, particle=None):
+    def _interpolate(self, time, z, y, x, particle=None):
         """Interpolate spatial field values."""
         try:
             if self.grid.zdim == 1:
-                val = self._interpolator2D(time, z, y, x, ti, particle=particle)
+                val = self._interpolator2D(time, z, y, x, particle=particle)
             else:
-                val = self._interpolator3D(time, z, y, x, ti, particle=particle)
+                val = self._interpolator3D(time, z, y, x, particle=particle)
 
             if np.isnan(val):
                 # Detect Out-of-bounds sampling and raise exception
@@ -593,8 +602,8 @@ class Field:
             e = add_note(e, f"Error interpolating field '{self.name}'.", before=True)
             raise e
 
-    def _time_index(self, time):
-        """Find the index in the time array associated with a given time.
+    def _search_time_index(self, time):
+        """Find and return the index and relative coordinate in the time array associated with a given time.
 
         Note that we normalize to either the first or the last index
         if the sampled value is outside the time value range.
@@ -606,13 +615,24 @@ class Field:
         if time_index.all():
             # If given time > last known field time, use
             # the last field frame without interpolation
-            return len(self.grid.time) - 1
+            ti = len(self.grid.time) - 1
         elif np.logical_not(time_index).all():
             # If given time < any time in the field, use
             # the first field frame without interpolation
-            return 0
+            ti = 0
+        else:
+            ti = time_index.argmin() - 1 if time_index.any() else 0
+        if self.grid.tdim == 1:
+            tau = 0
+        elif ti == len(self.grid.time) - 1:
+            tau = 1
         else:
-            return time_index.argmin() - 1 if time_index.any() else 0
+            tau = (
+                (time - self.grid.time[ti]) / (self.grid.time[ti + 1] - self.grid.time[ti])
+                if self.grid.time[ti] != self.grid.time[ti + 1]
+                else 0
+            )
+        return tau, ti
 
     def _check_velocitysampling(self):
         if self.name in ["U", "V", "W"]:
@@ -639,20 +659,10 @@ class Field:
         conversion to the result. Note that we defer to
         scipy.interpolate to perform spatial interpolation.
         """
-        ti = self._time_index(time)
         if self.gridindexingtype == "croco" and self not in [self.fieldset.H, self.fieldset.Zeta]:
             z = _croco_from_z_to_sigma_scipy(self.fieldset, time, z, y, x, particle=particle)
-        if ti < self.grid.tdim - 1 and time > self.grid.time[ti]:
-            f0 = self._spatial_interpolation(time, z, y, x, ti, particle=particle)
-            f1 = self._spatial_interpolation(time, z, y, x, ti + 1, particle=particle)
-            t0 = self.grid.time[ti]
-            t1 = self.grid.time[ti + 1]
-            value = f0 + (f1 - f0) * ((time - t0) / (t1 - t0))
-        else:
-            # Skip temporal interpolation if time is outside
-            # of the defined time range or if we have hit an
-            # exact value in the time array.
-            value = self._spatial_interpolation(self.grid.time[ti], z, y, x, ti, particle=particle)
+
+        value = self._interpolate(time, z, y, x, particle=particle)
 
         if applyConversion:
             return self.units.to_target(value, z, y, x)
@@ -805,9 +815,9 @@ class VectorField:
             and np.allclose(grid1.time_full, grid2.time_full)
         )
 
-    def spatial_c_grid_interpolation2D(self, ti, z, y, x, time, particle=None, applyConversion=True):
+    def c_grid_interpolation2D(self, time, z, y, x, particle=None, applyConversion=True):
         grid = self.U.grid
-        (_, eta, xsi, zi, yi, xi) = self.U._search_indices(time, z, y, x, ti, particle=particle)
+        (tau, _, eta, xsi, ti, zi, yi, xi) = self.U._search_indices(time, z, y, x, particle=particle)
 
         if grid._gtype in [GridType.RectilinearSGrid, GridType.RectilinearZGrid]:
             px = np.array([grid.lon[xi], grid.lon[xi + 1], grid.lon[xi + 1], grid.lon[xi]])
@@ -827,59 +837,68 @@ class VectorField:
         c2 = i_u._geodetic_distance(py[1], py[2], px[1], px[2], grid.mesh, np.dot(i_u.phi2D_lin(eta, 1.0), py))
         c3 = i_u._geodetic_distance(py[2], py[3], px[2], px[3], grid.mesh, np.dot(i_u.phi2D_lin(1.0, xsi), py))
         c4 = i_u._geodetic_distance(py[3], py[0], px[3], px[0], grid.mesh, np.dot(i_u.phi2D_lin(eta, 0.0), py))
-        if grid.zdim == 1:
-            if self.gridindexingtype == "nemo":
-                U0 = self.U.data[ti, yi + 1, xi] * c4
-                U1 = self.U.data[ti, yi + 1, xi + 1] * c2
-                V0 = self.V.data[ti, yi, xi + 1] * c1
-                V1 = self.V.data[ti, yi + 1, xi + 1] * c3
-            elif self.gridindexingtype in ["mitgcm", "croco"]:
-                U0 = self.U.data[ti, yi, xi] * c4
-                U1 = self.U.data[ti, yi, xi + 1] * c2
-                V0 = self.V.data[ti, yi, xi] * c1
-                V1 = self.V.data[ti, yi + 1, xi] * c3
-        else:
-            if self.gridindexingtype == "nemo":
-                U0 = self.U.data[ti, zi, yi + 1, xi] * c4
-                U1 = self.U.data[ti, zi, yi + 1, xi + 1] * c2
-                V0 = self.V.data[ti, zi, yi, xi + 1] * c1
-                V1 = self.V.data[ti, zi, yi + 1, xi + 1] * c3
-            elif self.gridindexingtype in ["mitgcm", "croco"]:
-                U0 = self.U.data[ti, zi, yi, xi] * c4
-                U1 = self.U.data[ti, zi, yi, xi + 1] * c2
-                V0 = self.V.data[ti, zi, yi, xi] * c1
-                V1 = self.V.data[ti, zi, yi + 1, xi] * c3
-        U = (1 - xsi) * U0 + xsi * U1
-        V = (1 - eta) * V0 + eta * V1
-        rad = np.pi / 180.0
-        deg2m = 1852 * 60.0
-        if applyConversion:
-            meshJac = (deg2m * deg2m * math.cos(rad * y)) if grid.mesh == "spherical" else 1
-        else:
-            meshJac = deg2m if grid.mesh == "spherical" else 1
-
-        jac = i_u._compute_jacobian_determinant(py, px, eta, xsi) * meshJac
-
-        u = (
-            (-(1 - eta) * U - (1 - xsi) * V) * px[0]
-            + ((1 - eta) * U - xsi * V) * px[1]
-            + (eta * U + xsi * V) * px[2]
-            + (-eta * U + (1 - xsi) * V) * px[3]
-        ) / jac
-        v = (
-            (-(1 - eta) * U - (1 - xsi) * V) * py[0]
-            + ((1 - eta) * U - xsi * V) * py[1]
-            + (eta * U + xsi * V) * py[2]
-            + (-eta * U + (1 - xsi) * V) * py[3]
-        ) / jac
-        if isinstance(u, da.core.Array):
-            u = u.compute()
-            v = v.compute()
+
+        def _calc_UV(ti, yi, xi):
+            if grid.zdim == 1:
+                if self.gridindexingtype == "nemo":
+                    U0 = self.U.data[ti, yi + 1, xi] * c4
+                    U1 = self.U.data[ti, yi + 1, xi + 1] * c2
+                    V0 = self.V.data[ti, yi, xi + 1] * c1
+                    V1 = self.V.data[ti, yi + 1, xi + 1] * c3
+                elif self.gridindexingtype in ["mitgcm", "croco"]:
+                    U0 = self.U.data[ti, yi, xi] * c4
+                    U1 = self.U.data[ti, yi, xi + 1] * c2
+                    V0 = self.V.data[ti, yi, xi] * c1
+                    V1 = self.V.data[ti, yi + 1, xi] * c3
+            else:
+                if self.gridindexingtype == "nemo":
+                    U0 = self.U.data[ti, zi, yi + 1, xi] * c4
+                    U1 = self.U.data[ti, zi, yi + 1, xi + 1] * c2
+                    V0 = self.V.data[ti, zi, yi, xi + 1] * c1
+                    V1 = self.V.data[ti, zi, yi + 1, xi + 1] * c3
+                elif self.gridindexingtype in ["mitgcm", "croco"]:
+                    U0 = self.U.data[ti, zi, yi, xi] * c4
+                    U1 = self.U.data[ti, zi, yi, xi + 1] * c2
+                    V0 = self.V.data[ti, zi, yi, xi] * c1
+                    V1 = self.V.data[ti, zi, yi + 1, xi] * c3
+            U = (1 - xsi) * U0 + xsi * U1
+            V = (1 - eta) * V0 + eta * V1
+            rad = np.pi / 180.0
+            deg2m = 1852 * 60.0
+            if applyConversion:
+                meshJac = (deg2m * deg2m * math.cos(rad * y)) if grid.mesh == "spherical" else 1
+            else:
+                meshJac = deg2m if grid.mesh == "spherical" else 1
+
+            jac = i_u._compute_jacobian_determinant(py, px, eta, xsi) * meshJac
+
+            u = (
+                (-(1 - eta) * U - (1 - xsi) * V) * px[0]
+                + ((1 - eta) * U - xsi * V) * px[1]
+                + (eta * U + xsi * V) * px[2]
+                + (-eta * U + (1 - xsi) * V) * px[3]
+            ) / jac
+            v = (
+                (-(1 - eta) * U - (1 - xsi) * V) * py[0]
+                + ((1 - eta) * U - xsi * V) * py[1]
+                + (eta * U + xsi * V) * py[2]
+                + (-eta * U + (1 - xsi) * V) * py[3]
+            ) / jac
+            if isinstance(u, da.core.Array):
+                u = u.compute()
+                v = v.compute()
+            return (u, v)
+
+        u, v = _calc_UV(ti, yi, xi)
+        if should_calculate_next_ti(ti, tau, self.U.grid.tdim):
+            ut1, vt1 = _calc_UV(ti + 1, yi, xi)
+            u = (1 - tau) * u + tau * ut1
+            v = (1 - tau) * v + tau * vt1
         return (u, v)
 
-    def spatial_c_grid_interpolation3D_full(self, ti, z, y, x, time, particle=None):
+    def c_grid_interpolation3D_full(self, time, z, y, x, particle=None):
         grid = self.U.grid
-        (zeta, eta, xsi, zi, yi, xi) = self.U._search_indices(time, z, y, x, ti, particle=particle)
+        (tau, zeta, eta, xsi, ti, zi, yi, xi) = self.U._search_indices(time, z, y, x, particle=particle)
 
         if grid._gtype in [GridType.RectilinearSGrid, GridType.RectilinearZGrid]:
             px = np.array([grid.lon[xi], grid.lon[xi + 1], grid.lon[xi + 1], grid.lon[xi]])
@@ -932,6 +951,14 @@ class VectorField:
         w0 = self.W.data[ti, zi, yi + 1, xi + 1]
         w1 = self.W.data[ti, zi + 1, yi + 1, xi + 1]
 
+        if should_calculate_next_ti(ti, tau, self.U.grid.tdim):
+            u0 = (1 - tau) * u0 + tau * self.U.data[ti + 1, zi, yi + 1, xi]
+            u1 = (1 - tau) * u1 + tau * self.U.data[ti + 1, zi, yi + 1, xi + 1]
+            v0 = (1 - tau) * v0 + tau * self.V.data[ti + 1, zi, yi, xi + 1]
+            v1 = (1 - tau) * v1 + tau * self.V.data[ti + 1, zi, yi + 1, xi + 1]
+            w0 = (1 - tau) * w0 + tau * self.W.data[ti + 1, zi, yi + 1, xi + 1]
+            w1 = (1 - tau) * w1 + tau * self.W.data[ti + 1, zi + 1, yi + 1, xi + 1]
+
         U0 = u0 * i_u.jacobian3D_lin_face(pz, py, px, zeta, eta, 0, "zonal", grid.mesh)
         U1 = u1 * i_u.jacobian3D_lin_face(pz, py, px, zeta, eta, 1, "zonal", grid.mesh)
         V0 = v0 * i_u.jacobian3D_lin_face(pz, py, px, zeta, 0, xsi, "meridional", grid.mesh)
@@ -1082,7 +1109,7 @@ class VectorField:
             w = w.compute()
         return (u, v, w)
 
-    def spatial_c_grid_interpolation3D(self, ti, z, y, x, time, particle=None, applyConversion=True):
+    def c_grid_interpolation3D(self, ti, z, y, x, time, particle=None, applyConversion=True):
         """Perform C grid interpolation in 3D. ::
 
             +---+---+---+
@@ -1099,11 +1126,11 @@ class VectorField:
         Curvilinear grids are treated properly, since the element is projected to a rectilinear parent element.
         """
         if self.U.grid._gtype in [GridType.RectilinearSGrid, GridType.CurvilinearSGrid]:
-            (u, v, w) = self.spatial_c_grid_interpolation3D_full(ti, z, y, x, time, particle=particle)
+            (u, v, w) = self.c_grid_interpolation3D_full(time, z, y, x, particle=particle)
         else:
             if self.gridindexingtype == "croco":
                 z = _croco_from_z_to_sigma_scipy(self.fieldset, time, z, y, x, particle=particle)
-            (u, v) = self.spatial_c_grid_interpolation2D(ti, z, y, x, time, particle=particle)
+            (u, v) = self.c_grid_interpolation2D(time, z, y, x, particle=particle)
             w = self.W.eval(time, z, y, x, particle=particle, applyConversion=False)
             if applyConversion:
                 w = self.W.units.to_target(w, z, y, x)
@@ -1121,8 +1148,8 @@ class VectorField:
             else:
                 return True
 
-    def spatial_slip_interpolation(self, ti, z, y, x, time, particle=None, applyConversion=True):
-        (zeta, eta, xsi, zi, yi, xi) = self.U._search_indices(time, z, y, x, ti, particle=particle)
+    def slip_interpolation(self, time, z, y, x, particle=None, applyConversion=True):
+        (_, zeta, eta, xsi, ti, zi, yi, xi) = self.U._search_indices(time, z, y, x, particle=particle)
         di = ti if self.U.grid.zdim == 1 else zi  # general third dimension
 
         f_u, f_v, f_w = 1, 1, 1
@@ -1223,66 +1250,22 @@ class VectorField:
             return u, v
 
     def eval(self, time, z, y, x, particle=None, applyConversion=True):
-        if self.U.interp_method not in ["cgrid_velocity", "partialslip", "freeslip"]:
+        if self.U.interp_method in ["partialslip", "freeslip"]:
+            return self.slip_interpolation(time, z, y, x, particle=particle, applyConversion=applyConversion)
+
+        if self.U.interp_method not in ["cgrid_velocity"]:
             u = self.U.eval(time, z, y, x, particle=particle, applyConversion=False)
             v = self.V.eval(time, z, y, x, particle=particle, applyConversion=False)
             if applyConversion:
                 u = self.U.units.to_target(u, z, y, x)
                 v = self.V.units.to_target(v, z, y, x)
-            if "3D" in self.vector_type:
-                w = self.W.eval(time, z, y, x, particle=particle, applyConversion=False)
-                if applyConversion:
-                    w = self.W.units.to_target(w, z, y, x)
-                return (u, v, w)
-            else:
-                return (u, v)
+        elif self.U.interp_method == "cgrid_velocity":
+            (u, v) = self.c_grid_interpolation2D(time, z, y, x, particle=particle, applyConversion=applyConversion)
+        if "3D" in self.vector_type:
+            w = self.W.eval(time, z, y, x, particle=particle, applyConversion=applyConversion)
+            return (u, v, w)
         else:
-            interp = {
-                "cgrid_velocity": {
-                    "2D": self.spatial_c_grid_interpolation2D,
-                    "3D": self.spatial_c_grid_interpolation3D,
-                },
-                "partialslip": {"2D": self.spatial_slip_interpolation, "3D": self.spatial_slip_interpolation},
-                "freeslip": {"2D": self.spatial_slip_interpolation, "3D": self.spatial_slip_interpolation},
-            }
-            grid = self.U.grid
-            ti = self.U._time_index(time)
-            if ti < grid.tdim - 1 and time > grid.time[ti]:
-                t0 = grid.time[ti]
-                t1 = grid.time[ti + 1]
-                if "3D" in self.vector_type:
-                    (u0, v0, w0) = interp[self.U.interp_method]["3D"](
-                        ti, z, y, x, time, particle=particle, applyConversion=applyConversion
-                    )
-                    (u1, v1, w1) = interp[self.U.interp_method]["3D"](
-                        ti + 1, z, y, x, time, particle=particle, applyConversion=applyConversion
-                    )
-                    w = w0 + (w1 - w0) * ((time - t0) / (t1 - t0))
-                else:
-                    (u0, v0) = interp[self.U.interp_method]["2D"](
-                        ti, z, y, x, time, particle=particle, applyConversion=applyConversion
-                    )
-                    (u1, v1) = interp[self.U.interp_method]["2D"](
-                        ti + 1, z, y, x, time, particle=particle, applyConversion=applyConversion
-                    )
-                u = u0 + (u1 - u0) * ((time - t0) / (t1 - t0))
-                v = v0 + (v1 - v0) * ((time - t0) / (t1 - t0))
-                if "3D" in self.vector_type:
-                    return (u, v, w)
-                else:
-                    return (u, v)
-            else:
-                # Skip temporal interpolation if time is outside
-                # of the defined time range or if we have hit an
-                # exact value in the time array.
-                if "3D" in self.vector_type:
-                    return interp[self.U.interp_method]["3D"](
-                        ti, z, y, x, grid.time[ti], particle=particle, applyConversion=applyConversion
-                    )
-                else:
-                    return interp[self.U.interp_method]["2D"](
-                        ti, z, y, x, grid.time[ti], particle=particle, applyConversion=applyConversion
-                    )
+            return (u, v)
 
     def __getitem__(self, key):
         try:
diff --git a/parcels/tools/_helpers.py b/parcels/tools/_helpers.py
index c49bbef9..d3dd7fb7 100644
--- a/parcels/tools/_helpers.py
+++ b/parcels/tools/_helpers.py
@@ -152,3 +152,8 @@ def timedelta_to_float(dt: float | timedelta | np.timedelta64) -> float:
     if isinstance(dt, np.timedelta64):
         return float(dt / np.timedelta64(1, "s"))
     return float(dt)
+
+
+def should_calculate_next_ti(ti: int, tau: float, tdim: int):
+    """Check if the time is beyond the last time in the field"""
+    return np.greater(tau, 0) and ti < tdim - 1


Test patch:
diff --git a/tests/test_advection.py b/tests/test_advection.py
index b926c25e..07ae63ad 100644
--- a/tests/test_advection.py
+++ b/tests/test_advection.py
@@ -231,7 +231,7 @@ def test_conversion_3DCROCO():
 
 
 @pytest.mark.v4alpha
-@pytest.mark.xfail(reason="When refactoring fieldfilebuffer croco support was dropped. This will be fixed in v4.")
+@pytest.mark.xfail(reason="CROCO 3D interpolation is not yet implemented correctly in v4. ")
 def test_advection_3DCROCO():
     fieldset = FieldSet.from_modulefile(TEST_DATA / "fieldset_CROCO3D.py")
 
diff --git a/tests/test_fieldset.py b/tests/test_fieldset.py
index f7174513..26618605 100644
--- a/tests/test_fieldset.py
+++ b/tests/test_fieldset.py
@@ -452,7 +452,7 @@ def test_fieldset_write(tmp_zarrfile):
     def UpdateU(particle, fieldset, time):  # pragma: no cover
         tmp1, tmp2 = fieldset.UV[particle]
         _, yi, xi = fieldset.U.unravel_index(particle.ei)
-        ti = fieldset.U._time_index(time)
+        _, ti = fieldset.U._search_time_index(time)
         fieldset.U.data[ti, yi, xi] += 1
         fieldset.U.grid.time[0] = time
 
diff --git a/tests/test_grids.py b/tests/test_grids.py
index 21a3d2ea..7a18010b 100644
--- a/tests/test_grids.py
+++ b/tests/test_grids.py
@@ -457,7 +457,8 @@ def test_cgrid_uniform_2dvel(time):
     assert (pset[0].meridional - 1) < 1e-6
 
 
-@pytest.mark.parametrize("vert_mode", ["zlev", "slev1", "slev2"])
+@pytest.mark.v4alpha
+@pytest.mark.parametrize("vert_mode", ["zlev"])  # , "slev1", "slev2"])  # TODO v4: time-varying depth not supported yet
 @pytest.mark.parametrize("time", [True, False])
 def test_cgrid_uniform_3dvel(vert_mode, time):
     lon = np.array([[0, 2], [0.4, 1.5]])
diff --git a/tests/test_interpolation.py b/tests/test_interpolation.py
index 07c290ed..9b764815 100644
--- a/tests/test_interpolation.py
+++ b/tests/test_interpolation.py
@@ -68,9 +68,9 @@ def data_2d():
 )
 def test_raw_2d_interpolation(data_2d, func, eta, xsi, expected):
     """Test the 2D interpolation functions on the raw arrays."""
-    ti = 0
+    tau, ti = 0, 0
     yi, xi = 1, 1
-    ctx = interpolation.InterpolationContext2D(data_2d, eta, xsi, ti, yi, xi)
+    ctx = interpolation.InterpolationContext2D(data_2d, tau, eta, xsi, ti, yi, xi)
     assert func(ctx) == expected