From 5134f9929a87cbc479965c11cdb3acb9259ffa50 Mon Sep 17 00:00:00 2001
From: Erik van Sebille <e.vansebille@uu.nl>
Date: Thu, 22 Oct 2026 17:17:10 +0200
Subject: [PATCH 01/15] First attempt at cleaning up unit tests

---
 src/parcels/_core/particlefile.py | 10 +++++---
 src/parcels/_repr_utils.py        |  2 +-
 tests/test_particlefile.py        | 41 +++++++++++++++++++------------
 3 files changed, 32 insertions(+), 21 deletions(-)

diff --git a/src/parcels/_core/particlefile.py b/src/parcels/_core/particlefile.py
index fb02885c6..ff397f17f 100644
--- a/src/parcels/_core/particlefile.py
+++ b/src/parcels/_core/particlefile.py
@@ -166,7 +166,6 @@ def write(self, pset: ParticleSet | ParticleSetView, t, indices=None):
             )
 
         if isinstance(t, (np.timedelta64, np.datetime64)):
-            print(t)
             if self._time_interval is None:
                 raise ValueError(
                     "The time interval is unknown, so datetime/timedelta times cannot be converted. "
@@ -174,10 +173,13 @@ def write(self, pset: ParticleSet | ParticleSetView, t, indices=None):
                 )
             t = timedelta_to_float(t - self._time_interval.left)
         vars_to_write = _get_vars_to_write(pclass)
-        if indices is None:
-            indices_to_write = _to_write_particles(particle_data, t)
-        else:
+        if indices is not None:
             indices_to_write = indices
+        elif isinstance(pset, ParticleSetView):
+            # A view already selects the particles to write, so `t` is not used to filter
+            indices_to_write = np.flatnonzero(pset._index)
+        else:
+            indices_to_write = _to_write_particles(particle_data, t)
 
         self._writer.write_table(
             pa.table({v.name: pa.array(particle_data[v.name][indices_to_write]) for v in vars_to_write}),
diff --git a/src/parcels/_repr_utils.py b/src/parcels/_repr_utils.py
index 159494d31..21197334d 100644
--- a/src/parcels/_repr_utils.py
+++ b/src/parcels/_repr_utils.py
@@ -104,7 +104,7 @@ def particleset_repr(pset: ParticleSet) -> str:
 
 def particlesetview_repr(pview: Any) -> str:
     """Return a pretty repr for ParticleSetView"""
-    time_string = "not_yet_set" if pview.time is None or np.isnan(pview.time) else f"{pview.time:f}"
+    time_string = "not_yet_set" if pview.t is None or np.isnan(pview.t) else f"{pview.t:f}"
     out = f"P[{pview.particle_id}]: time={time_string}, z={pview.z:f}, y={pview.y:f}, x={pview.x:f}"
     vars = [v.name for v in pview._pclass.variables if v.to_write is True and v.name not in ["z", "y", "x", "time"]]
     for var in vars:
diff --git a/tests/test_particlefile.py b/tests/test_particlefile.py
index 1bec5a89d..2d1b83f73 100755
--- a/tests/test_particlefile.py
+++ b/tests/test_particlefile.py
@@ -21,6 +21,7 @@
     StatusCode,
     Variable,
     convert,
+    particlefile_to_v3_zarr,
 )
 from parcels._core.particle import get_default_particle
 from parcels._core.particlefile import get_schema
@@ -156,44 +157,52 @@ def test_write_dtypes_pfile(fieldset, tmp_parquet):
         assert tab[f"v_{d.__name__}"].type == pa.from_numpy_dtype(d)
 
 
-@pytest.mark.skip(reason="Pending ParticleFile refactor; see issue #2386")
 @pytest.mark.parametrize("dt", [-np.timedelta64(1, "s"), np.timedelta64(1, "s")])
 @pytest.mark.parametrize("maxvar", [2, 4, 10])
-def test_pset_repeated_release_delayed_adding_deleting(fieldset, tmp_parquet, dt, maxvar):
+def test_pset_repeated_release_delayed_adding_deleting(fieldset, tmp_parquet, tmp_zarr, dt, maxvar):
     """Tests that if particles are released and deleted based on age that resulting output file is correct."""
     npart = 10
     fieldset.add_context("maxvar", maxvar)
 
-    MyParticle = Particle.add_variable(
-        [Variable("sample_var", initial=0.0), Variable("v_once", dtype=np.float64, initial=0.0)]
-    )
+    MyParticle = Particle.add_variable([Variable("sample_var", initial=0.0)])
 
     pset = ParticleSet(
         fieldset,
+        pclass=MyParticle,
         x=np.zeros(npart),
         y=np.zeros(npart),
-        pclass=MyParticle,
         t=fieldset.time_interval.left + [np.timedelta64(i + 1, "s") for i in range(npart)],
     )
     pfile = ParticleFile(tmp_parquet, outputdt=abs(dt))
 
     def IncrLon(particles, fieldset):  # pragma: no cover
+        ptcls_to_write = particles[(particles.sample_var == 0) & (particles.t != 0)]
+        if len(ptcls_to_write.t) > 0:
+            pfile.write(ptcls_to_write, ptcls_to_write.t)
         particles.sample_var += 1.0
-        particles.state = np.where(
-            particles.sample_var > fieldset.maxvar,
-            StatusCode.Delete,
-            particles.state,
-        )
+        particles.state[particles.sample_var > fieldset.maxvar] = StatusCode.Delete
 
-    for _ in range(npart):
-        pset.execute(IncrLon, dt=dt, runtime=np.timedelta64(1, "s"), output_file=pfile)
+    pset.execute(IncrLon, dt=dt, runtime=np.timedelta64(npart, "s"), output_file=pfile)
+
+    df = parcels.read_particlefile(tmp_parquet)
+    for traj in df.partition_by("particle_id", maintain_order=True):
+        print(traj["sample_var"][:])
+
+    samplevar = df["sample_var"][:]
+    np.testing.assert_array_equal(samplevar.shape, (npart, min(maxvar, npart) + 1))
+
+    particlefile_to_v3_zarr(tmp_parquet, tmp_zarr)
+    ds = xr.open_zarr(tmp_zarr)
 
-    ds = xr.open_zarr(tmp_parquet)
     samplevar = ds["sample_var"][:]
-    assert samplevar.shape == (npart, min(maxvar, npart + 1))
+    print(samplevar.values)
+    print(ds["time"].values)
+    print("")
+    np.testing.assert_array_equal(samplevar.shape, (npart, min(maxvar, npart) + 1))
     # test whether samplevar[:, k] = k
     for k in range(samplevar.shape[1]):
-        assert np.allclose([p for p in samplevar[:, k] if np.isfinite(p)], k + 1)
+        print([p.values for p in samplevar[1:, k] if np.isfinite(p)])
+        assert np.allclose([p for p in samplevar[1:, k] if np.isfinite(p)], k + 1)
 
 
 def test_file_warnings(fieldset, tmp_parquet):

From 9197c57e6631ee995ef8a87799f6fd15a8c1abca Mon Sep 17 00:00:00 2001
From: Erik van Sebille <e.vansebille@uu.nl>
Date: Fri, 23 Oct 2026 16:24:27 +0200
Subject: [PATCH 02/15] Removing xfail on passing test in
 test_particleset_execute

---
 tests/test_particleset_execute.py | 13 ++++++-------
 1 file changed, 6 insertions(+), 7 deletions(-)

diff --git a/tests/test_particleset_execute.py b/tests/test_particleset_execute.py
index 9343307a8..a8c5fd3e6 100644
--- a/tests/test_particleset_execute.py
+++ b/tests/test_particleset_execute.py
@@ -571,24 +571,23 @@ def test_uxstommelgyre_multiparticle_pset_execute():
     )
 
 
-@pytest.mark.xfail(reason="Output file not implemented yet")
 def test_uxstommelgyre_pset_execute_output():
     ds = datasets_unstructured["stommel_gyre_delaunay"]
     fieldset = FieldSet.from_ugrid_conventions(ds)
     pset = ParticleSet(
         fieldset,
-        x=[30.0],
-        y=[5.0],
-        z=[50.0],
-        t=[0.0],
+        x=30.0,
+        y=5.0,
+        z=50.0,
+        t=np.timedelta64(0, "s"),
         pclass=Particle,
     )
     output_file = ParticleFile(
-        name="stommel_uxarray_particles.zarr",  # the file name
+        path="stommel_uxarray_particles.parquet",  # the file name
         outputdt=np.timedelta64(5, "m"),  # the time step of the outputs
     )
     pset.execute(
-        runtime=np.timedelta64(10, "m"), dt=np.timedelta64(60, "s"), kernels=AdvectionEE, output_file=output_file
+        runtime=np.timedelta64(10, "m"), dt=np.timedelta64(150, "s"), kernels=AdvectionEE, output_file=output_file
     )
 
 

From fc051470fe1945530c4f9d9b8f49aa92c00a90af Mon Sep 17 00:00:00 2001
From: Erik van Sebille <e.vansebille@uu.nl>
Date: Fri, 23 Oct 2026 16:24:59 +0200
Subject: [PATCH 03/15] Fixing localize function to use new position dictionary

---
 src/parcels/_core/xgrid.py | 45 +++++++++++++++++++++-----------------
 tests/test_xgrid.py        | 10 ++++++---
 2 files changed, 32 insertions(+), 23 deletions(-)

diff --git a/src/parcels/_core/xgrid.py b/src/parcels/_core/xgrid.py
index 80bada25d..edcb2599e 100644
--- a/src/parcels/_core/xgrid.py
+++ b/src/parcels/_core/xgrid.py
@@ -257,8 +257,8 @@ def get_axis_dim(self, axis: ptyping.XgridAxis) -> int:
         return get_cell_count_along_dim(self._ds, self.sgrid_metadata.vertical_dimensions[0])
 
     def localize(
-        self, position: dict[ptyping.XgridAxis, tuple[int, float]], dims: list[str]
-    ) -> dict[str, tuple[int, float]]:
+        self, position: dict[ptyping.XgridAxis, dict[str, np.ndarray]], dims: list[str]
+    ) -> dict[str, dict[str, np.ndarray]]:
         """
         Uses the grid context (i.e., the staggering of the grid) to convert a position relative
         to the f-points in the grid to a position relative to the staggered grid the array
@@ -271,24 +271,28 @@ def localize(
         Parameters
         ----------
         position : dict
-            A mapping of the axis to a tuple of (index, barycentric coordinate) for the
-            f-points in the grid.
+            A mapping of the axis to a dict with the arrays ``"index"`` and ``"bcoord"`` (barycentric
+            coordinate) for the f-points in the grid, as returned by ``search``.
         dims : list[str]
             A list of dimension names that the DataArray is defined on. This is used to determine
             the staggering of the grid and which axis each dimension corresponds to.
 
         Returns
         -------
-        dict[str, tuple[int, float]]
-            A mapping of the dimension names to a tuple of (index, barycentric coordinate) for
+        dict[str, dict[str, np.ndarray]]
+            A mapping of the dimension names to a dict with the arrays ``"index"`` and ``"bcoord"`` for
             the staggered grid the DataArray is defined on.
 
         Example
         -------
-        >>> position = {"X": (5, 0.51), "Y": (10, 0.25), "Z": (3, 0.75)}
+        >>> position = {
+        ...     "X": {"index": np.array([5]), "bcoord": np.array([0.51])},
+        ...     "Y": {"index": np.array([10]), "bcoord": np.array([0.25])},
+        ...     "Z": {"index": np.array([3]), "bcoord": np.array([0.75])},
+        ... }
         >>> dims = ["time", "depth", "YC", "XC"]
         >>> grid.localize(position, dims)
-        {'depth': (3, 0.75), 'YC': (9, 0.75), 'XC': (5, 0.01)}
+        {'depth': {'index': array([3]), 'bcoord': array([0.75])}, 'YC': {'index': array([9]), 'bcoord': array([0.75])}, 'XC': {'index': array([5]), 'bcoord': array([0.01])}}
         """
         dim_to_axis = _get_dim_to_axis_mapping(self.sgrid_metadata)
         axis_to_var = {dim_to_axis[dim]: dim for dim in dims if dim in dim_to_axis}
@@ -297,12 +301,13 @@ def localize(
         }
         return {
             axis_to_var[axis]: _convert_center_pos_to_fpoint(
-                index=index,
-                bcoord=bcoord,
+                index=pos["index"],
+                bcoord=pos["bcoord"],
                 position=var_positions[axis],
                 f_point_position=self._fpoint_info[axis],
             )
-            for axis, (index, bcoord) in position.items()
+            for axis, pos in position.items()
+            if axis in axis_to_var
         }
 
     @property
@@ -591,11 +596,11 @@ def assert_valid_depth(da_depth):
 
 def _convert_center_pos_to_fpoint(
     *,
-    index: int,
-    bcoord: float,
+    index: np.ndarray,
+    bcoord: np.ndarray,
     position: ptyping.GridPosition,
     f_point_position: sgrid.Padding,
-) -> tuple[int, float]:
+) -> dict[str, np.ndarray]:
     """Converts a physical position relative to the cell edges defined in the grid to be relative to the face center.
 
     This is used to "localize" a position to be relative to the staggered grid at which the field is defined, so that
@@ -604,16 +609,16 @@ def _convert_center_pos_to_fpoint(
     This also handles different model input cell edges and centers are staggered in different directions (e.g., with NEMO and MITgcm).
     """
     if position != "face":  # Data is already defined on the F nodes
-        return index, bcoord
+        return {"index": index, "bcoord": bcoord}
 
     bcoord = bcoord - 0.5
-    if bcoord < 0:
-        bcoord += 1.0
-        index -= 1
+    wrapped = bcoord < 0
+    bcoord = np.where(wrapped, bcoord + 1.0, bcoord)
+    index = np.where(wrapped, index - 1, index)
 
     # Correct relative to the f-point position
     # Padding.BOTH was "inner", Padding.LOW was "right" in xgcm vocabulary
     if f_point_position in (sgrid.Padding.BOTH, sgrid.Padding.LOW):
-        index += 1
+        index = index + 1
 
-    return index, bcoord
+    return {"index": index, "bcoord": bcoord}
diff --git a/tests/test_xgrid.py b/tests/test_xgrid.py
index 465284215..6553a6d91 100644
--- a/tests/test_xgrid.py
+++ b/tests/test_xgrid.py
@@ -304,7 +304,6 @@ def test_search_1d_array_some_out_of_bounds(search_1d, array, x, expected_xi):
     np.testing.assert_array_equal(xi, expected_xi)
 
 
-@pytest.mark.xfail(reason="grid.localize not yet adapted to xi, xsi, yi, eta, zi, zeta position keys")
 @pytest.mark.parametrize(
     "ds, da_name, expected",
     [
@@ -352,9 +351,14 @@ def test_search_1d_array_some_out_of_bounds(search_1d, array, x, expected_xi):
 )
 def test_xgrid_localize_zero_position(ds, da_name, expected):
     """Test localize function using left and right datasets."""
-    grid = XGrid.from_dataset(ds)
+    fieldset = FieldSet.from_sgrid_conventions(ds)
+    grid = fieldset.models[0].grid
+
     da = ds[da_name]
     position = grid.search(0, 0, 0)
 
     local_position = grid.localize(position, da.dims)
-    assert local_position == expected, f"Expected {expected}, got {local_position}"
+    assert local_position.keys() == expected.keys()
+    for dim, (index, bcoord) in expected.items():
+        np.testing.assert_array_equal(local_position[dim]["index"], [index])
+        np.testing.assert_allclose(local_position[dim]["bcoord"], [bcoord])

From c9e430fd7bcfb8a0b23dfb09aaeca24c1d8788be Mon Sep 17 00:00:00 2001
From: Erik van Sebille <e.vansebille@uu.nl>
Date: Fri, 23 Oct 2026 16:26:19 +0200
Subject: [PATCH 04/15] Renaming localize to UK spelling localise

---
 docs/user_guide/examples/tutorial_schism.ipynb | 4 ++--
 src/parcels/_core/xgrid.py                     | 6 +++---
 tests/test_xgrid.py                            | 6 +++---
 3 files changed, 8 insertions(+), 8 deletions(-)

diff --git a/docs/user_guide/examples/tutorial_schism.ipynb b/docs/user_guide/examples/tutorial_schism.ipynb
index 406938c96..a81b4cb97 100644
--- a/docs/user_guide/examples/tutorial_schism.ipynb
+++ b/docs/user_guide/examples/tutorial_schism.ipynb
@@ -26,7 +26,7 @@
     "1. The mesh coordinates are in a **projected coordinate system (meters)**, not longitude/latitude, so\n",
     "   we build the grid as a *flat* (Cartesian) mesh.\n",
     "2. Velocities are stored at **mesh nodes** over a vertical column of layers, ordered **bottom → surface**.\n",
-    "3. SCHISM uses a **localized vertical grid (LSC2)**: the number of valid vertical levels varies from\n",
+    "3. SCHISM uses a **localised vertical grid (LSC2)**: the number of valid vertical levels varies from\n",
     "   node to node, and levels below the seabed are stored as `NaN`. We use this to flag particles that\n",
     "   end up below the local bathymetry and stop advecting them.\n",
     "\n",
@@ -147,7 +147,7 @@
     "\n",
     "* **Layer ordering.** SCHISM stores levels **bottom → surface**, while Parcels expects depth increasing\n",
     "  downward from the surface, so we reverse the layer axis.\n",
-    "* **Vertical coordinate.** SCHISM's LSC2 vertical grid (*Localized Sigma Coordinates with Shaved\n",
+    "* **Vertical coordinate.** SCHISM's LSC2 vertical grid (*Localised Sigma Coordinates with Shaved\n",
     "  cells*) varies with horizontal position: each node has its own layer depths, and even its own\n",
     "  *number* of levels (the true depths are written to a separate `zCoordinates` output, not used here).\n",
     "  **Parcels does not currently support a vertical grid that varies with lateral position.** `UxGrid`\n",
diff --git a/src/parcels/_core/xgrid.py b/src/parcels/_core/xgrid.py
index edcb2599e..b58783086 100644
--- a/src/parcels/_core/xgrid.py
+++ b/src/parcels/_core/xgrid.py
@@ -256,7 +256,7 @@ def get_axis_dim(self, axis: ptyping.XgridAxis) -> int:
         assert self.sgrid_metadata.vertical_dimensions is not None
         return get_cell_count_along_dim(self._ds, self.sgrid_metadata.vertical_dimensions[0])
 
-    def localize(
+    def localise(
         self, position: dict[ptyping.XgridAxis, dict[str, np.ndarray]], dims: list[str]
     ) -> dict[str, dict[str, np.ndarray]]:
         """
@@ -291,7 +291,7 @@ def localize(
         ...     "Z": {"index": np.array([3]), "bcoord": np.array([0.75])},
         ... }
         >>> dims = ["time", "depth", "YC", "XC"]
-        >>> grid.localize(position, dims)
+        >>> grid.localise(position, dims)
         {'depth': {'index': array([3]), 'bcoord': array([0.75])}, 'YC': {'index': array([9]), 'bcoord': array([0.75])}, 'XC': {'index': array([5]), 'bcoord': array([0.01])}}
         """
         dim_to_axis = _get_dim_to_axis_mapping(self.sgrid_metadata)
@@ -603,7 +603,7 @@ def _convert_center_pos_to_fpoint(
 ) -> dict[str, np.ndarray]:
     """Converts a physical position relative to the cell edges defined in the grid to be relative to the face center.
 
-    This is used to "localize" a position to be relative to the staggered grid at which the field is defined, so that
+    This is used to "localise" a position to be relative to the staggered grid at which the field is defined, so that
     it can be easily interpolated.
 
     This also handles different model input cell edges and centers are staggered in different directions (e.g., with NEMO and MITgcm).
diff --git a/tests/test_xgrid.py b/tests/test_xgrid.py
index 6553a6d91..2a365e0eb 100644
--- a/tests/test_xgrid.py
+++ b/tests/test_xgrid.py
@@ -349,15 +349,15 @@ def test_search_1d_array_some_out_of_bounds(search_1d, array, x, expected_xi):
         ),
     ],
 )
-def test_xgrid_localize_zero_position(ds, da_name, expected):
-    """Test localize function using left and right datasets."""
+def test_xgrid_localise_zero_position(ds, da_name, expected):
+    """Test localise function using left and right datasets."""
     fieldset = FieldSet.from_sgrid_conventions(ds)
     grid = fieldset.models[0].grid
 
     da = ds[da_name]
     position = grid.search(0, 0, 0)
 
-    local_position = grid.localize(position, da.dims)
+    local_position = grid.localise(position, da.dims)
     assert local_position.keys() == expected.keys()
     for dim, (index, bcoord) in expected.items():
         np.testing.assert_array_equal(local_position[dim]["index"], [index])

From 81e5d9a69b9cb59f89c35450c832d72336d6e2c6 Mon Sep 17 00:00:00 2001
From: Erik van Sebille <e.vansebille@uu.nl>
Date: Fri, 23 Oct 2026 16:26:34 +0200
Subject: [PATCH 05/15] Remove irrelevant test

---
 tests/test_particlefile.py | 50 --------------------------------------
 1 file changed, 50 deletions(-)

diff --git a/tests/test_particlefile.py b/tests/test_particlefile.py
index 5de1c465c..f5b13a865 100755
--- a/tests/test_particlefile.py
+++ b/tests/test_particlefile.py
@@ -8,7 +8,6 @@
 import pyarrow as pa
 import pyarrow.parquet as pq
 import pytest
-import xarray as xr
 from re_assert import Matches
 
 import parcels.tutorial
@@ -21,7 +20,6 @@
     StatusCode,
     Variable,
     convert,
-    particlefile_to_v3_zarr,
 )
 from parcels._core.particle import get_default_particle
 from parcels._core.particlefile import get_schema
@@ -157,54 +155,6 @@ def test_write_dtypes_pfile(fieldset, tmp_parquet):
         assert tab[f"v_{d.__name__}"].type == pa.from_numpy_dtype(d)
 
 
-@pytest.mark.parametrize("dt", [-np.timedelta64(1, "s"), np.timedelta64(1, "s")])
-@pytest.mark.parametrize("maxvar", [2, 4, 10])
-def test_pset_repeated_release_delayed_adding_deleting(fieldset, tmp_parquet, tmp_zarr, dt, maxvar):
-    """Tests that if particles are released and deleted based on age that resulting output file is correct."""
-    npart = 10
-    fieldset.add_context("maxvar", maxvar)
-
-    MyParticle = Particle.add_variable([Variable("sample_var", initial=0.0)])
-
-    pset = ParticleSet(
-        fieldset,
-        pclass=MyParticle,
-        x=np.zeros(npart),
-        y=np.zeros(npart),
-        t=fieldset.time_interval.left + [np.timedelta64(i + 1, "s") for i in range(npart)],
-    )
-    pfile = ParticleFile(tmp_parquet, outputdt=abs(dt))
-
-    def IncrLon(particles, fieldset):  # pragma: no cover
-        ptcls_to_write = particles[(particles.sample_var == 0) & (particles.t != 0)]
-        if len(ptcls_to_write.t) > 0:
-            pfile.write(ptcls_to_write, ptcls_to_write.t)
-        particles.sample_var += 1.0
-        particles.state[particles.sample_var > fieldset.maxvar] = StatusCode.Delete
-
-    pset.execute(IncrLon, dt=dt, runtime=np.timedelta64(npart, "s"), output_file=pfile)
-
-    df = parcels.read_particlefile(tmp_parquet)
-    for traj in df.partition_by("particle_id", maintain_order=True):
-        print(traj["sample_var"][:])
-
-    samplevar = df["sample_var"][:]
-    np.testing.assert_array_equal(samplevar.shape, (npart, min(maxvar, npart) + 1))
-
-    particlefile_to_v3_zarr(tmp_parquet, tmp_zarr)
-    ds = xr.open_zarr(tmp_zarr)
-
-    samplevar = ds["sample_var"][:]
-    print(samplevar.values)
-    print(ds["time"].values)
-    print("")
-    np.testing.assert_array_equal(samplevar.shape, (npart, min(maxvar, npart) + 1))
-    # test whether samplevar[:, k] = k
-    for k in range(samplevar.shape[1]):
-        print([p.values for p in samplevar[1:, k] if np.isfinite(p)])
-        assert np.allclose([p for p in samplevar[1:, k] if np.isfinite(p)], k + 1)
-
-
 def test_file_warnings(fieldset, tmp_parquet):
     pset = ParticleSet(fieldset, x=[0, 0], y=[0, 0], t=[np.timedelta64(0, "s"), np.timedelta64(1, "s")])
     pfile = ParticleFile(tmp_parquet, outputdt=np.timedelta64(2, "s"))

From 12016fe990c0bc938189737e389c49c1b4f87f8e Mon Sep 17 00:00:00 2001
From: Erik van Sebille <e.vansebille@uu.nl>
Date: Mon, 26 Oct 2026 15:44:54 +0100
Subject: [PATCH 06/15] Using temporary file for output

---
 tests/test_particleset_execute.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/tests/test_particleset_execute.py b/tests/test_particleset_execute.py
index a8c5fd3e6..93c7893e3 100644
--- a/tests/test_particleset_execute.py
+++ b/tests/test_particleset_execute.py
@@ -571,7 +571,7 @@ def test_uxstommelgyre_multiparticle_pset_execute():
     )
 
 
-def test_uxstommelgyre_pset_execute_output():
+def test_uxstommelgyre_pset_execute_output(tmp_parquet):
     ds = datasets_unstructured["stommel_gyre_delaunay"]
     fieldset = FieldSet.from_ugrid_conventions(ds)
     pset = ParticleSet(
@@ -583,7 +583,7 @@ def test_uxstommelgyre_pset_execute_output():
         pclass=Particle,
     )
     output_file = ParticleFile(
-        path="stommel_uxarray_particles.parquet",  # the file name
+        path=tmp_parquet,  # the file name
         outputdt=np.timedelta64(5, "m"),  # the time step of the outputs
     )
     pset.execute(

From f71ac911a274728014788f55ba8a7a6ecc15c916 Mon Sep 17 00:00:00 2001
From: Erik van Sebille <e.vansebille@uu.nl>
Date: Mon, 26 Oct 2026 16:42:31 +0100
Subject: [PATCH 07/15] Updating unit tests in xgrid

---
 tests/test_xgrid.py | 36 +++++++++++++-----------------------
 1 file changed, 13 insertions(+), 23 deletions(-)

diff --git a/tests/test_xgrid.py b/tests/test_xgrid.py
index 2a365e0eb..483ca25aa 100644
--- a/tests/test_xgrid.py
+++ b/tests/test_xgrid.py
@@ -50,13 +50,10 @@ def assert_equal(actual, expected):
         assert_allclose(actual, expected)
 
 
-@pytest.mark.skip(
-    "Needs updating after refactoring from https://github.com/Parcels-code/Parcels/pull/2646"
-)  # TODO: Remove or replace
 @pytest.mark.parametrize("ds", [datasets["ds_2d_left"]])
 def test_grid_init_param_types(ds):
-    with pytest.raises(ValueError, match="Invalid value 'invalid'. Valid options are.*"):
-        XGrid.from_dataset(ds, mesh="invalid")
+    with pytest.raises(TypeError, match="mesh must be a FlatMesh or SphericalMesh object. Got mesh.*"):
+        FieldSet.from_sgrid_conventions(ds, mesh="invalid")
 
 
 def test_xgrid_repr(fieldset):
@@ -85,20 +82,17 @@ def test_xgrid_mesh_inferred_from_units(ds, mesh):
     assert isinstance(grid.mesh, mesh)
 
 
-@pytest.mark.skip(
-    "Needs updating after refactoring from https://github.com/Parcels-code/Parcels/pull/2646"
-)  # TODO: Remove or replace
 @pytest.mark.parametrize("ds", [datasets["ds_2d_left"]])
 def test_transpose_xfield_data_to_tzyx(ds):
     da = ds["data_g"]
-    grid = XGrid.from_dataset(ds)
+    grid = FieldSet.from_sgrid_conventions(ds).data_g.grid
 
     all_combinations = (itertools.combinations(da.dims, n) for n in range(len(da.dims)))
     all_combinations = itertools.chain(*all_combinations)
     for subset_dims in all_combinations:
         isel = {dim: 0 for dim in subset_dims}
         da_subset = da.isel(isel, drop=True)
-        da_test = _transpose_xfield_data_to_tzyx(da_subset, grid.xgcm_grid)
+        da_test = _transpose_xfield_data_to_tzyx(da_subset, grid.sgrid_metadata)
         utils.assert_valid_field_data(da_test, grid)
 
 
@@ -114,42 +108,41 @@ def test_invalid_xgrid_field_array():
     ...
 
 
-@pytest.mark.skip(
-    "Needs updating after refactoring from https://github.com/Parcels-code/Parcels/pull/2646"
-)  # TODO: Remove or replace
 def test_invalid_lon_lat():
     """Stress test the grid initialiser by creating incompatible datasets that test the edge cases"""
     ds = datasets["ds_2d_left"].copy()
     ds["lon"], ds["lat"] = xr.broadcast(ds["YC"], ds["XC"])
+    ds["lon"].attrs["units"] = "meters"
+    ds["lat"].attrs["units"] = "meters"
 
     with pytest.raises(
         ValueError,
-        match=r".*is defined on the center of the grid, but must be defined on the F points\.",
+        match=r".*is defined on the faces of the grid, but must be defined on the F nodes\.",
     ):
-        XGrid.from_dataset(ds)
+        FieldSet.from_sgrid_conventions(ds)
 
     ds = datasets["ds_2d_left"].copy()
     ds["lon"], _ = xr.broadcast(ds["YG"], ds["XG"])
+    ds["lon"].attrs["units"] = "meters"
     with pytest.raises(
         ValueError,
         match=r".*have different dimensionalities\.",
     ):
-        XGrid.from_dataset(ds)
+        FieldSet.from_sgrid_conventions(ds)
 
     ds = datasets["ds_2d_left"].copy()
     ds["lon"], ds["lat"] = xr.broadcast(ds["YG"], ds["XG"])
     ds["lon"], ds["lat"] = ds["lon"].transpose(), ds["lat"].transpose()
+    ds["lon"].attrs["units"] = "meters"
+    ds["lat"].attrs["units"] = "meters"
 
     with pytest.raises(
         ValueError,
         match=r".*must be defined on the X and Y axes and transposed to have dimensions in order of Y, X\.",
     ):
-        XGrid.from_dataset(ds)
+        FieldSet.from_sgrid_conventions(ds)
 
 
-@pytest.mark.skip(
-    "Needs updating after refactoring from https://github.com/Parcels-code/Parcels/pull/2646"
-)  # TODO: Remove or replace
 def test_invalid_depth():
     ds = datasets["ds_2d_left"].copy()
     ds = ds.reindex({"ZG": ds.ZG[::-1]})
@@ -184,9 +177,6 @@ def test_dim_with_duplicate_axis():
         FieldSet.from_sgrid_conventions(ds)
 
 
-# TODO restructure: Look into the test below
-# remove: eval with timedelta64 time fails; TimeInterval.is_all_time_in_interval expects float (seconds) not timedelta64; time arg type requirement changed
-@pytest.mark.skip("remove: see comment above")
 @pytest.mark.parametrize("ds", [datasets["ds_2d_left"]])
 def test_vertical1D_field(ds):
     ds = ds.drop_vars(set(ds.data_vars) - {"grid"})

From 8dce2a8cdc217c6d2a085c674f8b81e776453e2a Mon Sep 17 00:00:00 2001
From: Erik van Sebille <e.vansebille@uu.nl>
Date: Mon, 26 Oct 2026 16:44:21 +0100
Subject: [PATCH 08/15] Removing unit test that I think is not relevant anymore

As I think that duplicate axes are now supported by model.py
---
 tests/test_xgrid.py | 19 +------------------
 1 file changed, 1 insertion(+), 18 deletions(-)

diff --git a/tests/test_xgrid.py b/tests/test_xgrid.py
index 483ca25aa..0e52b505c 100644
--- a/tests/test_xgrid.py
+++ b/tests/test_xgrid.py
@@ -20,7 +20,7 @@
     XGrid,
     _transpose_xfield_data_to_tzyx,
 )
-from parcels._datasets.structured.generic import X, Y, Z, datasets, datasets_sgrid
+from parcels._datasets.structured.generic import X, Y, Z, datasets
 from tests import utils
 
 GridTestCase = namedtuple("GridTestCase", ["ds", "attr", "expected"])
@@ -148,23 +148,6 @@ def test_invalid_depth():
     ds = ds.reindex({"ZG": ds.ZG[::-1]})
 
     with pytest.raises(ValueError, match="Depth DataArray .* must be strictly increasing*"):
-        XGrid.from_dataset(ds)
-
-
-@pytest.mark.skip(
-    "Needs updating after refactoring from https://github.com/Parcels-code/Parcels/pull/2646"
-)  # TODO: I think we can just rely on the SGRID metadata for this (which already has robust error messaging). How should discrepencies between SGRID and axis attr be handled?
-def test_dim_with_duplicate_axis():
-    ds = datasets_sgrid["ds_2d_padded_low"].copy()
-
-    # Add an extra Z axis
-    ds = ds[["data_g", "grid"]]
-    ds["data_g"] = ds["data_g"].expand_dims("vertical_dimensions_dim2", 1)
-
-    z = ds["vertical_dimensions_dim2"]
-    z.attrs.update({"axis": "Z"})
-    ds["vertical_dimensions_dim2"] = z
-
     # TODO: Clean up this attribute setting (really, this should be on the source datasets)
     lon = ds["lon"]
     lat = ds["lat"]

From c3bf3be52701a8b6750d67b5c327b9d36cb9c1d2 Mon Sep 17 00:00:00 2001
From: Erik van Sebille <e.vansebille@uu.nl>
Date: Mon, 26 Oct 2026 16:44:37 +0100
Subject: [PATCH 09/15] Update test_xgrid.py

---
 tests/test_xgrid.py | 9 ---------
 1 file changed, 9 deletions(-)

diff --git a/tests/test_xgrid.py b/tests/test_xgrid.py
index 0e52b505c..2e224d967 100644
--- a/tests/test_xgrid.py
+++ b/tests/test_xgrid.py
@@ -148,15 +148,6 @@ def test_invalid_depth():
     ds = ds.reindex({"ZG": ds.ZG[::-1]})
 
     with pytest.raises(ValueError, match="Depth DataArray .* must be strictly increasing*"):
-    # TODO: Clean up this attribute setting (really, this should be on the source datasets)
-    lon = ds["lon"]
-    lat = ds["lat"]
-    lon.attrs.update({"units": "metres"})
-    lat.attrs.update({"units": "metres"})
-    ds["lon"] = lon
-    ds["lat"] = lat
-
-    with pytest.raises(ValueError, match="Two dimensions .*provide values in the axis direction 'Z'."):
         FieldSet.from_sgrid_conventions(ds)
 
 

From 023c33904f14b5571106a97953b6636029b31579 Mon Sep 17 00:00:00 2001
From: Erik van Sebille <e.vansebille@uu.nl>
Date: Mon, 26 Oct 2026 16:57:49 +0100
Subject: [PATCH 10/15] removing unused/unnecessary tests

---
 tests/test_particleset_execute.py | 6 ------
 tests/test_xgrid.py               | 5 -----
 2 files changed, 11 deletions(-)

diff --git a/tests/test_particleset_execute.py b/tests/test_particleset_execute.py
index 93c7893e3..2b37e325a 100644
--- a/tests/test_particleset_execute.py
+++ b/tests/test_particleset_execute.py
@@ -390,12 +390,6 @@ def FieldAccessOutsideTime(particles, fieldset):  # pragma: no cover
         pset.execute(FieldAccessOutsideTime, runtime=np.timedelta64(1, "D"), dt=np.timedelta64(10, "D"))
 
 
-def test_raise_grid_searching_error(): ...
-
-
-def test_raise_general_error(): ...
-
-
 def test_errorinterpolation(fieldset):
     class NaNInterpolator(VectorInterpolator):  # pragma: no cover
         def interp(self, particle_positions, grid_positions, field):
diff --git a/tests/test_xgrid.py b/tests/test_xgrid.py
index 2e224d967..17b2efccd 100644
--- a/tests/test_xgrid.py
+++ b/tests/test_xgrid.py
@@ -103,11 +103,6 @@ def test_xgrid_get_axis_dim(fieldset):
     assert grid.get_axis_dim("X") == X - 1
 
 
-def test_invalid_xgrid_field_array():
-    """Stress test initialiser by creating incompatible datasets that test the edge cases"""
-    ...
-
-
 def test_invalid_lon_lat():
     """Stress test the grid initialiser by creating incompatible datasets that test the edge cases"""
     ds = datasets["ds_2d_left"].copy()

From 2b4f6cf9a6efb0433f77dd21de5802c99231475c Mon Sep 17 00:00:00 2001
From: Erik van Sebille <e.vansebille@uu.nl>
Date: Mon, 26 Oct 2026 16:58:58 +0100
Subject: [PATCH 11/15] Expanding grid combing tests

---
 tests/test_fieldset.py | 18 +++++++++++++++---
 1 file changed, 15 insertions(+), 3 deletions(-)

diff --git a/tests/test_fieldset.py b/tests/test_fieldset.py
index 10f7c7045..482c19c33 100644
--- a/tests/test_fieldset.py
+++ b/tests/test_fieldset.py
@@ -17,7 +17,7 @@
 from parcels._core.mesh import SphericalMesh
 from parcels._core.model import _default_vector_field_components, _infer_mesh_from_sgrid_data
 from parcels._core.utils.time import TimeInterval
-from parcels._datasets.structured.generated import simple_UV_sigma_dataset
+from parcels._datasets.structured.generated import simple_UV_dataset, simple_UV_sigma_dataset
 from parcels._datasets.structured.generic import datasets as datasets_structured
 from parcels._datasets.structured.generic import datasets_sgrid
 from parcels._datasets.structured.strategies import reduced_dims_dataset, sgrid_datasets_varying_dims
@@ -295,7 +295,11 @@ def test_datetime_to_msg(input_, expected):
 
 def test_fieldset_samegrids_UV():
     """Test that if a simple fieldset with U and V is created, that only one grid object is defined."""
-    ...
+    ds = simple_UV_dataset(mesh_str="flat")
+    fieldset = FieldSet.from_sgrid_conventions(ds)
+
+    assert fieldset.U.grid is fieldset.V.grid
+    assert len(fieldset.gridset) == 1
 
 
 def test_fieldset_grid_deduplication():
@@ -305,7 +309,15 @@ def test_fieldset_grid_deduplication():
     When grid deduplication is actually implemented, this might need to be refactored
     into multiple tests (/more might be needed).
     """
-    ...
+    ds = simple_UV_dataset(mesh_str="flat")
+    fieldset = FieldSet.from_sgrid_conventions(ds)
+
+    ds2 = simple_UV_dataset(mesh_str="flat").rename({"U": "U2", "V": "V2"})
+    fieldset += FieldSet.from_sgrid_conventions(ds2)
+    assert fieldset.U.grid is fieldset.V.grid
+    assert fieldset.U2.grid is fieldset.V2.grid
+    assert fieldset.U.grid is not fieldset.U2.grid
+    assert fieldset.U.grid is not fieldset.V2.grid
 
 
 def test_fieldset_from_icon():

From a7b8af18ef9fca0e359a77ece83cd8d8316c706a Mon Sep 17 00:00:00 2001
From: Erik van Sebille <e.vansebille@uu.nl>
Date: Mon, 26 Oct 2026 16:59:25 +0100
Subject: [PATCH 12/15] Adding test for mismatched calendar errors

---
 src/parcels/_core/fieldset.py |  2 +-
 tests/test_fieldset.py        | 13 +++++++++++--
 2 files changed, 12 insertions(+), 3 deletions(-)

diff --git a/src/parcels/_core/fieldset.py b/src/parcels/_core/fieldset.py
index 9355c145e..ca6f3aea7 100644
--- a/src/parcels/_core/fieldset.py
+++ b/src/parcels/_core/fieldset.py
@@ -70,12 +70,12 @@ def __init__(self, models: list[ModelData]):
         for model in models:
             if not isinstance(model, ModelData):
                 raise ValueError(f"Expected `model` to be a ModelData object. Got {model}")
-        # assert_compatible_calendars(fields)
 
         self.models = models
         self.constant_model: StructuredModelData | None = None
         self._fields: dict[str, Field | VectorField] | None = None
         self.reconstruct_fields()
+        assert_compatible_calendars(self.fields.values())
         self.context: dict[str, float] = {}
         assert_models_have_same_mesh(self.models)
 
diff --git a/tests/test_fieldset.py b/tests/test_fieldset.py
index 482c19c33..ea126e664 100644
--- a/tests/test_fieldset.py
+++ b/tests/test_fieldset.py
@@ -13,7 +13,7 @@
 import parcels.tutorial
 import tests
 from parcels import Particle, ParticleFile, ParticleSet, Variable, convert
-from parcels._core.fieldset import FieldSet, IncompatibleMeshesException, _datetime_to_msg
+from parcels._core.fieldset import CalendarError, FieldSet, IncompatibleMeshesException, _datetime_to_msg
 from parcels._core.mesh import SphericalMesh
 from parcels._core.model import _default_vector_field_components, _infer_mesh_from_sgrid_data
 from parcels._core.utils.time import TimeInterval
@@ -272,7 +272,16 @@ def test_multi_model_nonoverlapping_time_interval():
 
 def test_fieldset_add_incompatible_calendars():
     # tests the adding of fieldsets that have incompatible calendars
-    ...
+    ds1 = datasets_structured["ds_2d_left"][["U_A_grid", "V_A_grid", "grid"]]
+    fieldset = FieldSet.from_sgrid_conventions(ds1)
+
+    ds2 = ds1.copy().rename({"U_A_grid": "U2", "V_A_grid": "V2"})
+    # Replace the numpy datetime64 time axis by an array of cftime objects with a different calendar
+    times = [cftime.DatetimeNoLeap(2000, 1, 1 + i) for i in range(ds2["time"].size)]
+    ds2["time"] = (ds2["time"].dims, np.array(times, dtype=object), ds2["time"].attrs)
+
+    with pytest.raises(CalendarError):
+        fieldset + FieldSet.from_sgrid_conventions(ds2)
 
 
 @pytest.mark.parametrize(

From 0118e1c569f470df2787ca086d1b38eb7a41db4d Mon Sep 17 00:00:00 2001
From: Erik van Sebille <e.vansebille@uu.nl>
Date: Mon, 26 Oct 2026 17:07:43 +0100
Subject: [PATCH 13/15] Adding test for incompatible calendars in VectorField

---
 tests/test_field.py | 11 ++++++++++-
 1 file changed, 10 insertions(+), 1 deletion(-)

diff --git a/tests/test_field.py b/tests/test_field.py
index 2272a5cb3..ede3e865e 100644
--- a/tests/test_field.py
+++ b/tests/test_field.py
@@ -2,6 +2,7 @@
 
 import numpy as np
 import pytest
+import xarray as xr
 from re_assert import Matches
 
 from parcels import Field, VectorField
@@ -82,7 +83,15 @@ def test_field_time_interval():
 
 def test_vectorfield_init_different_time_intervals():
     # Tests that a VectorField raises a ValueError if the component fields have different time domains.
-    ...
+    dsU = simple_UV_dataset(mesh_str="flat")[["U", "grid"]]
+    dsV = simple_UV_dataset(mesh_str="flat")[["V"]]
+
+    dsV = dsV.rename({"time": "timeV"})
+    dsV["timeV"] = dsV["timeV"] + np.timedelta64(1, "D")
+
+    ds = xr.merge([dsU, dsV])
+    with pytest.raises(ValueError, match="Multiple time coordinates"):
+        FieldSet.from_sgrid_conventions(ds)
 
 
 def test_field_invalid_interpolator():

From 24eb9177fa6e53531764c7305b334ef5f77d4c28 Mon Sep 17 00:00:00 2001
From: Erik van Sebille <e.vansebille@uu.nl>
Date: Mon, 26 Oct 2026 17:07:53 +0100
Subject: [PATCH 14/15] Removing irrelevant tests

---
 tests/test_field.py | 12 ------------
 1 file changed, 12 deletions(-)

diff --git a/tests/test_field.py b/tests/test_field.py
index ede3e865e..9edb6042c 100644
--- a/tests/test_field.py
+++ b/tests/test_field.py
@@ -251,15 +251,3 @@ def test_field_eval_out_of_bounds_unstructured(field_name, location, expected):
         match="Some interpolated values are out-of-bounds. These values are set to 0. Treat carefully.",
     ):
         np.testing.assert_allclose(field.eval(*location), expected)
-
-
-def test_field_unstructured_grid_creation(): ...
-
-
-def test_field_interpolation(): ...
-
-
-def test_field_interpolation_out_of_spatial_bounds(): ...
-
-
-def test_field_interpolation_out_of_time_bounds(): ...

From aec5edaf746c9ce7d8b55df8e69198b83443eaaf Mon Sep 17 00:00:00 2001
From: Erik van Sebille <e.vansebille@uu.nl>
Date: Mon, 26 Oct 2026 17:31:18 +0100
Subject: [PATCH 15/15] Adding explicit axis to timeV

Since it failed on min-versions test
---
 tests/test_field.py | 1 +
 1 file changed, 1 insertion(+)

diff --git a/tests/test_field.py b/tests/test_field.py
index 9edb6042c..7def3ca8b 100644
--- a/tests/test_field.py
+++ b/tests/test_field.py
@@ -88,6 +88,7 @@ def test_vectorfield_init_different_time_intervals():
 
     dsV = dsV.rename({"time": "timeV"})
     dsV["timeV"] = dsV["timeV"] + np.timedelta64(1, "D")
+    dsV["timeV"].attrs["axis"] = "T"
 
     ds = xr.merge([dsU, dsV])
     with pytest.raises(ValueError, match="Multiple time coordinates"):