diff --git a/src/parcels/_core/index_search.py b/src/parcels/_core/index_search.py
index 5e2c223d0..5e84eee2e 100644
--- a/src/parcels/_core/index_search.py
+++ b/src/parcels/_core/index_search.py
@@ -92,56 +92,151 @@ def _search_time_index(field: Field, time: np.ndarray):
 
 
 def curvilinear_point_in_cell(grid, y: np.ndarray, x: np.ndarray, yi: np.ndarray, xi: np.ndarray):
-    xsi = eta = -1.0 * np.ones(len(x), dtype=float)
-    invA = np.array(
-        [
-            [1, 0, 0, 0],
-            [-1, 1, 0, 0],
-            [-1, 0, 0, 1],
-            [1, -1, 1, -1],
-        ]
-    )
+    """Locate points within 2D curvilinear grid cells.
 
-    px = np.array([grid.lon[yi, xi], grid.lon[yi, xi + 1], grid.lon[yi + 1, xi + 1], grid.lon[yi + 1, xi]])
+    For spherical grids, the bilinear inverse runs in a tangent plane on a sphere
+    through each cell's centroid — robust across the antimeridian and the
+    polar cap. For flat grids it runs directly in (lon, lat). In both
+    cases the returned (xsi, eta) are bilinear weights in the space the inverse
+    was solved in; downstream interpolators must agree on that same space.
+    """
+    clon = np.asarray(
+        [grid.lon[yi, xi], grid.lon[yi, xi + 1], grid.lon[yi + 1, xi + 1], grid.lon[yi + 1, xi]],
+        dtype=float,
+    )
+    clat = np.asarray(
+        [grid.lat[yi, xi], grid.lat[yi, xi + 1], grid.lat[yi + 1, xi + 1], grid.lat[yi + 1, xi]],
+        dtype=float,
+    )
 
     if grid._mesh == "spherical":
-        # Map grid and particle longitudes to [-180,180)
-        px = ((px + 180.0) % 360.0) - 180.0
-        x = ((x + 180.0) % 360.0) - 180.0
+        xsi, eta = _bilinear_inverse_tangent_plane(clon, clat, x, y)
+        is_in_cell = np.where((xsi >= 0) & (xsi <= 1) & (eta >= 0) & (eta <= 1), 1, 0)
+    else:
+        xsi, eta = _bilinear_inverse_latlon(clon, clat, x, y, apply_antimeridian_shift=False)
+        is_in_cell = np.where((xsi >= 0) & (xsi <= 1) & (eta >= 0) & (eta <= 1), 1, 0)
 
-        # Create a mask for antimeridian cells
-        lon_span = px.max(axis=0) - px.min(axis=0)
-        antimeridian_cell = lon_span > 180.0
+    return is_in_cell, np.column_stack((xsi, eta))
 
-        if np.any(antimeridian_cell):
-            # For any antimeridian cell ...
-            #  If particle longitude is closer to 180.0, then negative cell longitudes need to be bumped by +360
-            mask = (px < 0.0) & antimeridian_cell[np.newaxis, :] & (x[np.newaxis, :] > 0.0)
-            px[mask] += 360.0
-            #  If particle longitude is closer to -180.0, then positive cell longitudes need to be bumped by -360
-            mask = (px > 0.0) & antimeridian_cell[np.newaxis, :] & (x[np.newaxis, :] < 0.0)
-            px[mask] -= 360.0
 
-    py = np.array([grid.lat[yi, xi], grid.lat[yi, xi + 1], grid.lat[yi + 1, xi + 1], grid.lat[yi + 1, xi]])
+_invA = np.array(
+    [
+        [1, 0, 0, 0],
+        [-1, 1, 0, 0],
+        [-1, 0, 0, 1],
+        [1, -1, 1, -1],
+    ]
+)
+
 
-    a, b = np.dot(invA, px), np.dot(invA, py)
+def _bilinear_inverse(px, py, xq, yq):
+    """Solve for (xsi, eta) s.t. bilinear blend of (px, py) corners equals (xq, yq)."""
+    eta_init = -1.0 * np.ones(len(xq), dtype=float)
+    a, b = np.dot(_invA, px), np.dot(_invA, py)
     aa = a[3] * b[2] - a[2] * b[3]
-    bb = a[3] * b[0] - a[0] * b[3] + a[1] * b[2] - a[2] * b[1] + x * b[3] - y * a[3]
-    cc = a[1] * b[0] - a[0] * b[1] + x * b[1] - y * a[1]
+    bb = a[3] * b[0] - a[0] * b[3] + a[1] * b[2] - a[2] * b[1] + xq * b[3] - yq * a[3]
+    cc = a[1] * b[0] - a[0] * b[1] + xq * b[1] - yq * a[1]
     det2 = bb * bb - 4 * aa * cc
 
     with np.errstate(divide="ignore", invalid="ignore"):
-        det = np.where(det2 > 0, np.sqrt(det2), eta)
-        eta = np.where(abs(aa) < 1e-12, -cc / bb, np.where(det2 > 0, (-bb + det) / (2 * aa), eta))
-
+        det = np.where(det2 > 0, np.sqrt(det2), eta_init)
+        eta = np.where(abs(aa) < 1e-12, -cc / bb, np.where(det2 > 0, (-bb + det) / (2 * aa), eta_init))
         xsi = np.where(
             abs(a[1] + a[3] * eta) < 1e-12,
-            ((y - py[0]) / (py[1] - py[0]) + (y - py[3]) / (py[2] - py[3])) * 0.5,
-            (x - a[0] - a[2] * eta) / (a[1] + a[3] * eta),
+            ((yq - py[0]) / (py[1] - py[0]) + (yq - py[3]) / (py[2] - py[3])) * 0.5,
+            (xq - a[0] - a[2] * eta) / (a[1] + a[3] * eta),
         )
-    is_in_cell = np.where((xsi >= 0) & (xsi <= 1) & (eta >= 0) & (eta <= 1), 1, 0)
+    return xsi, eta
 
-    return is_in_cell, np.column_stack((xsi, eta))
+
+def _bilinear_inverse_latlon(clon, clat, x, y, apply_antimeridian_shift=True):
+    """Bilinear inverse in raw lat/lon, with antimeridian corner shifting."""
+    px = clon.copy()
+    xq = np.asarray(x, dtype=float).copy()
+    if apply_antimeridian_shift:
+        px = ((px + 180.0) % 360.0) - 180.0
+        xq = ((xq + 180.0) % 360.0) - 180.0
+        lon_span = px.max(axis=0) - px.min(axis=0)
+        antimeridian_cell = lon_span > 180.0
+        if np.any(antimeridian_cell):
+            mask = (px < 0.0) & antimeridian_cell[np.newaxis, :] & (xq[np.newaxis, :] > 0.0)
+            px[mask] += 360.0
+            mask = (px > 0.0) & antimeridian_cell[np.newaxis, :] & (xq[np.newaxis, :] < 0.0)
+            px[mask] -= 360.0
+    py = clat
+    yq = np.asarray(y, dtype=float)
+    return _bilinear_inverse(px, py, xq, yq)
+
+
+def _bilinear_inverse_tangent_plane(clon, clat, x, y):
+    """Bilinear inverse after tangent-plane projection through the cell centroid.
+
+    Avoids lat/lon singularities at poles and on the antimeridian.
+    """
+    px, py, xq, yq = _spherical_project_cell_and_query(clon, clat, x, y)
+    return _bilinear_inverse(px, py, xq, yq)
+
+
+def _spherical_project_cell_and_query(clon, clat, x, y):
+    """Project 4 cell corners and the query onto a plane defined by the cell.
+
+    The plane's basis (e_u, e_v) is built from the cell's edge-midpoint
+    difference vectors, then orthonormalized using Gramm-Schmidt
+
+    Inputs
+    ------
+    clon, clat : (4, N) arrays of corner lon/lat in degrees
+    x, y       : (N,) query lon/lat in degrees
+
+    Returns
+    -------
+    px, py : (4, N) projected corner coordinates (u, v)
+    xq, yq : (N,) projected query coordinates in the same (u, v)
+
+    Notes
+    -----
+    Orthogonal projection is 2-to-1 on the full sphere: a query and its
+    antipode map to the same (u, v). In practice this is harmless because the
+    downstream bilinear gate ``(xsi, eta) in [0, 1]`` only admits queries whose
+    projected magnitude is ≲ the cell's angular size.
+    """
+    cX, cY, cZ = _latlon_rad_to_xyz(np.deg2rad(clat), np.deg2rad(clon))  # (4, N) each
+    qX, qY, qZ = _latlon_rad_to_xyz(np.deg2rad(np.asarray(y, dtype=float)), np.deg2rad(np.asarray(x, dtype=float)))
+
+    # Cell-intrinsic tangent basis from edge-midpoint difference vectors.
+    # Corner order is CCW: c0=(yi,xi), c1=(yi,xi+1), c2=(yi+1,xi+1), c3=(yi+1,xi).
+    # u-axis (along xsi): right-edge midpoint - left-edge midpoint = (c1+c2) - (c0+c3)
+    # v-axis (along eta): top-edge   midpoint - bottom-edge midpoint = (c2+c3) - (c0+c1)
+    ux = (cX[1] + cX[2]) - (cX[0] + cX[3])
+    uy = (cY[1] + cY[2]) - (cY[0] + cY[3])
+    uz = (cZ[1] + cZ[2]) - (cZ[0] + cZ[3])
+    u_norm = np.sqrt(ux * ux + uy * uy + uz * uz)
+    u_norm = np.where(u_norm == 0.0, 1.0, u_norm)
+    e_ux, e_uy, e_uz = ux / u_norm, uy / u_norm, uz / u_norm
+
+    vx = (cX[2] + cX[3]) - (cX[0] + cX[1])
+    vy = (cY[2] + cY[3]) - (cY[0] + cY[1])
+    vz = (cZ[2] + cZ[3]) - (cZ[0] + cZ[1])
+    # Gram-Schmidt: make e_v perpendicular to e_u so the bilinear inverse
+    # in (u, v) is well-conditioned even on skewed cells.
+    v_dot_eu = vx * e_ux + vy * e_uy + vz * e_uz
+    vx = vx - v_dot_eu * e_ux
+    vy = vy - v_dot_eu * e_uy
+    vz = vz - v_dot_eu * e_uz
+    v_norm = np.sqrt(vx * vx + vy * vy + vz * vz)
+    v_norm = np.where(v_norm == 0.0, 1.0, v_norm)
+    e_vx, e_vy, e_vz = vx / v_norm, vy / v_norm, vz / v_norm
+
+    # Orthogonal projection onto span(e_u, e_v): two dot products per point.
+    # Same pattern as uxgrid_point_in_cell (drop the component normal to the plane).
+    def _project(vx, vy, vz):
+        u = vx * e_ux + vy * e_uy + vz * e_uz
+        v = vx * e_vx + vy * e_vy + vz * e_vz
+        return u, v
+
+    px_u, px_v = _project(cX, cY, cZ)  # (4, N) each
+    xq_u, xq_v = _project(qX, qY, qZ)  # (N,) each
+    return px_u, px_v, xq_u, xq_v
 
 
 def _search_indices_curvilinear_2d(
diff --git a/src/parcels/interpolators/_xinterpolators.py b/src/parcels/interpolators/_xinterpolators.py
index 7c0eebdcb..725422d02 100644
--- a/src/parcels/interpolators/_xinterpolators.py
+++ b/src/parcels/interpolators/_xinterpolators.py
@@ -177,7 +177,6 @@ def CGrid_Velocity(
     yi, eta = grid_positions["Y"]["index"], grid_positions["Y"]["bcoord"]
     zi, zeta = grid_positions["Z"]["index"], grid_positions["Z"]["bcoord"]
     ti, tau = grid_positions["T"]["index"], grid_positions["T"]["bcoord"]
-    lon = particle_positions["lon"]
 
     U = vectorfield.U.data
     V = vectorfield.V.data
@@ -305,13 +304,6 @@ def _compute_corner_data(data, selection_dict) -> np.ndarray:
         u /= conversion
         v /= conversion
 
-    # check whether the grid conversion has been applied correctly
-    xx = (1 - xsi) * (1 - eta) * px[0] + xsi * (1 - eta) * px[1] + xsi * eta * px[2] + (1 - xsi) * eta * px[3]
-    dlon = xx - lon
-    if grid._mesh == "spherical":
-        dlon = ((dlon + 180.0) % 360.0) - 180.0
-    u = np.where(np.abs(dlon / lon) > 1e-4, np.nan, u)
-
     if vectorfield.W:
         W = vectorfield.W.data
 
diff --git a/tests/test_index_search.py b/tests/test_index_search.py
index bb7ec3f3b..b85e6e3b6 100644
--- a/tests/test_index_search.py
+++ b/tests/test_index_search.py
@@ -4,7 +4,7 @@
 
 import parcels.tutorial
 from parcels import Field, XGrid
-from parcels._core.index_search import _search_indices_curvilinear_2d
+from parcels._core.index_search import _latlon_rad_to_xyz, _search_indices_curvilinear_2d
 from parcels._datasets.structured.generic import datasets
 from parcels.interpolators import XLinear
 
@@ -73,6 +73,25 @@ def test_indexing_nemo_curvilinear():
     for i in range(lons.shape[0]):
         np.testing.assert_allclose(px[1, i], px[:, i], atol=5)
 
-    # Reconstruct lons values from cornerpoints
-    xx = (1 - xsi) * (1 - eta) * px[0] + xsi * (1 - eta) * px[1] + xsi * eta * px[2] + (1 - xsi) * eta * px[3]
-    np.testing.assert_allclose(xx, lons, atol=1e-6)
+    # Each query should have been located inside some cell: xsi, eta in [0, 1]
+    assert np.all((xsi >= 0) & (xsi <= 1)), f"xsi out of [0,1]: {xsi}"
+    assert np.all((eta >= 0) & (eta <= 1)), f"eta out of [0,1]: {eta}"
+
+    # Reconstruct query lat/lon by bilinear-blending the 4 corner 3D unit-sphere
+    # vectors with (xsi, eta) and renormalizing onto the unit sphere. Tolerance
+    # reflects the residual spherical curvature for a 1/4° NEMO cell.
+    clat = np.array([grid.lat[yi, xi], grid.lat[yi, xi + 1], grid.lat[yi + 1, xi + 1], grid.lat[yi + 1, xi]])
+    clon = np.array([grid.lon[yi, xi], grid.lon[yi, xi + 1], grid.lon[yi + 1, xi + 1], grid.lon[yi + 1, xi]])
+    cX, cY, cZ = _latlon_rad_to_xyz(np.deg2rad(clat), np.deg2rad(clon))
+    w = np.array([(1 - xsi) * (1 - eta), xsi * (1 - eta), xsi * eta, (1 - xsi) * eta])
+    x = np.sum(w * cX, axis=0)
+    y = np.sum(w * cY, axis=0)
+    z = np.sum(w * cZ, axis=0)
+    n = np.sqrt(x * x + y * y + z * z)
+    x, y, z = x / n, y / n, z / n
+    lat_recon = np.rad2deg(np.arcsin(z))
+    lon_recon = np.rad2deg(np.arctan2(y, x))
+    lons_wrapped = ((lons + 180) % 360) - 180
+    dlon = ((lon_recon - lons_wrapped + 180) % 360) - 180
+    np.testing.assert_allclose(lat_recon, lats, atol=1e-5)
+    np.testing.assert_allclose(dlon, 0.0, atol=1e-5)