diff --git a/budget/README.md b/budget/README.md
index 89453840..08a9c56d 100644
--- a/budget/README.md
+++ b/budget/README.md
@@ -1,7 +1,9 @@
 # To close budget on particles
 
-The pdf in this directory describes, in theory, how to close budget on Lagrangian particles using seaduck. A paper will (hopefully) be published soon on this, and there will be a proper way to cite then. 
+> We are activatively tidying up the code for budget closure. The process is quite slow because the boy need to prioritize on getting his PhD :P.
 
-The corresponding code exists, but is a lot less tidy than this package, and is therefore in a private repo. If you would like to use this functionality, please email me (Wenrui Jiang) at wjiang33@jhu.edu. I am happy to share the code, explain things, and potentially corroborate on your project. 
+The pdf in this directory describes, in theory, how to close budget on Lagrangian particles using seaduck. A paper will (hopefully) be published soon on this, and there will be a proper way to cite then.
 
-If you want to contribute to the package by helping me make the budget code public, also email me. I don't know how to thank you enough.  
+The corresponding code exists, but is a lot less tidy than this package, and is therefore in a private repo. If you would like to use this functionality, please email me (Wenrui Jiang) at wjiang33@jhu.edu. I am happy to share the code, explain things, and potentially corroborate on your project.
+
+If you want to contribute to the package by helping me make the budget code public, also email me. I don't know how to thank you enough.
diff --git a/docs/get_started.md b/docs/get_started.md
index beabddee..f5bafe8d 100644
--- a/docs/get_started.md
+++ b/docs/get_started.md
@@ -1,6 +1,6 @@
 # Get started
 
-For a extremely quick start, read this [one minute guide](./one_min_guide.md). Or, learn about how the objects and functions are related to each other [here](./network_of_object.md).
+For a extremely quick start, read this [one minute guide](./one_min_guide.ipynb). Or, learn about how the objects and functions are related to each other [here](./network_of_object.md).
 
 If you have questions about any function, you can always go to [API references](./public_api_reference.md). Or you could checkout the following examples. More example can be found [here](ideal_test.md) and [here](ocean_example.md)
 
diff --git a/pyproject.toml b/pyproject.toml
index 51be0116..6b1cacab 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -40,7 +40,9 @@ module = [
   "numba",
   "pooch",
   "scipy",
-  "scipy.interpolate"
+  "scipy.interpolate",
+  "xgcm",
+  "zarr"
 ]
 
 [tool.ruff]
diff --git a/seaduck/eulerian_budget.py b/seaduck/eulerian_budget.py
new file mode 100644
index 00000000..0a52ced0
--- /dev/null
+++ b/seaduck/eulerian_budget.py
@@ -0,0 +1,349 @@
+import numpy as np
+
+try:  # pragma: no cover
+    import xgcm
+except ImportError:  # pragma: no cover
+    pass
+
+
+def _raise_if_no_xgcm():
+    try:
+        import xgcm
+
+        xgcm
+    except ImportError:
+        raise ImportError(
+            "The python package xgcm is needed."
+            "You can install it with:"
+            "conda install -c xgcm"
+        )
+
+
+def create_ecco_grid(ds):
+    face_connections = {
+        "face": {
+            0: {"X": ((12, "Y", False), (3, "X", False)), "Y": (None, (1, "Y", False))},
+            1: {
+                "X": ((11, "Y", False), (4, "X", False)),
+                "Y": ((0, "Y", False), (2, "Y", False)),
+            },
+            2: {
+                "X": ((10, "Y", False), (5, "X", False)),
+                "Y": ((1, "Y", False), (6, "X", False)),
+            },
+            3: {"X": ((0, "X", False), (9, "Y", False)), "Y": (None, (4, "Y", False))},
+            4: {
+                "X": ((1, "X", False), (8, "Y", False)),
+                "Y": ((3, "Y", False), (5, "Y", False)),
+            },
+            5: {
+                "X": ((2, "X", False), (7, "Y", False)),
+                "Y": ((4, "Y", False), (6, "Y", False)),
+            },
+            6: {
+                "X": ((2, "Y", False), (7, "X", False)),
+                "Y": ((5, "Y", False), (10, "X", False)),
+            },
+            7: {
+                "X": ((6, "X", False), (8, "X", False)),
+                "Y": ((5, "X", False), (10, "Y", False)),
+            },
+            8: {
+                "X": ((7, "X", False), (9, "X", False)),
+                "Y": ((4, "X", False), (11, "Y", False)),
+            },
+            9: {"X": ((8, "X", False), None), "Y": ((3, "X", False), (12, "Y", False))},
+            10: {
+                "X": ((6, "Y", False), (11, "X", False)),
+                "Y": ((7, "Y", False), (2, "X", False)),
+            },
+            11: {
+                "X": ((10, "X", False), (12, "X", False)),
+                "Y": ((8, "Y", False), (1, "X", False)),
+            },
+            12: {
+                "X": ((11, "X", False), None),
+                "Y": ((9, "Y", False), (0, "X", False)),
+            },
+        }
+    }
+
+    grid = xgcm.Grid(
+        ds,
+        periodic=False,
+        face_connections=face_connections,
+        coords={
+            "X": {"center": "X", "left": "Xp1"},
+            "Y": {"center": "Y", "left": "Yp1"},
+            "Z": {"center": "Z", "left": "Zl"},
+            "time": {"center": "time", "inner": "time_midp"},
+        },
+    )
+    return grid
+
+
+def hor_div(tub, grid, xfluxname, yfluxname):
+    """Calculate horizontal divergence using xgcm.
+
+    Parameters
+    ----------
+    tub: sd.OceData or xr.Dataset
+        The dataset to calculate data from
+    grid: xgcm.Grid
+        The Grid of the dataset
+    xfluxname, yfluxname: string
+        The name of the variables corresponding to the horizontal fluxes
+        in concentration m^3/s
+    """
+    try:
+        tub["Vol"]
+    except KeyError:
+        tub._add_missing_vol()
+    xy_diff = grid.diff_2d_vector(
+        {"X": tub[xfluxname].fillna(0), "Y": tub[yfluxname].fillna(0)},
+        boundary="fill",
+        fill_value=0.0,
+    )
+    x_diff = xy_diff["X"]
+    y_diff = xy_diff["Y"]
+    hConv = -(x_diff + y_diff) / tub["Vol"]
+    return hConv
+
+
+def ver_div(tub, grid, zfluxname):
+    """Calculate horizontal divergence using xgcm.
+
+    Parameters
+    ----------
+    tub: sd.OceData or xr.Dataset
+        The dataset to calculate data from
+    grid: xgcm.Grid
+        The Grid of the dataset
+    xfluxname, yfluxname, zfluxname: string
+        The name of the variables corresponding to the fluxes
+        in concentration m^3/s
+    """
+    try:
+        tub["Vol"]
+    except KeyError:
+        tub._add_missing_vol()
+    vConv = (
+        grid.diff(tub[zfluxname].fillna(0), "Z", boundary="fill", fill_value=0.0)
+        / tub["Vol"]
+    )
+    return vConv
+
+
+def total_div(tub, grid, xfluxname, yfluxname, zfluxname):
+    """Calculate 3D divergence using xgcm.
+
+    Parameters
+    ----------
+    tub: sd.OceData or xr.Dataset
+        The dataset to calculate data from
+    grid: xgcm.Grid
+        The Grid of the dataset
+    zfluxname: string
+        The name of the variables corresponding to the vertical flux
+        in concentration m^3/s
+    """
+    hConv = hor_div(tub, grid, xfluxname, yfluxname)
+    vConv = ver_div(tub, grid, zfluxname)
+    return hConv + vConv
+
+
+def _slice_corner(array, fc, iy1, iy2, ix1, ix2):
+    left = np.minimum(ix1, ix2)
+    righ = np.maximum(ix1, ix2)
+    down = np.minimum(iy1, iy2)
+    uppp = np.maximum(iy1, iy2)
+    return array[..., fc, down : uppp + 1, left : righ + 1]
+
+
+def _right90(array):
+    return array[..., ::-1].transpose([0, 2, 1])
+
+
+def _left90(array):
+    return array[..., ::-1, :].transpose([0, 2, 1])
+
+
+def buffer_x_withface(s, face, lm, rm, tp):
+    """Create buffer zone in x direction for one face of an array.
+
+    Parameters
+    ----------
+    s: numpy.ndarray
+        the center field, the last dimension being X,
+        the second dimension being face.
+    face: int
+        which face to create buffer for.
+    lm: int
+        the size of the margin to the left.
+    rm: int
+        the size of the margin to the right.
+    tp: seaduck.Topology
+        the topology object of the
+    """
+    shape = list(s.shape)
+    shape.pop(1)
+    shape[-1] += lm + rm
+    xbuffer = np.zeros(shape)
+    xbuffer[..., lm:-rm] = s[:, face]
+    try:
+        fc1, iy1, ix1 = tp.ind_moves((face, tp.iymax, 0), [2 for i in range(lm)])
+        fc2, iy2, ix2 = tp.ind_moves((face, 0, 0), [2])
+        left = _slice_corner(s, fc1, iy1, iy2, ix1, ix2)
+    except IndexError:
+        left = np.zeros_like(xbuffer[..., :lm])
+
+    try:
+        fc3, iy3, ix3 = tp.ind_moves((face, tp.iymax, tp.ixmax), [3 for i in range(rm)])
+        fc4, iy4, ix4 = tp.ind_moves((face, 0, tp.ixmax), [3])
+        righ = _slice_corner(s, fc3, iy3, iy4, ix3, ix4)
+    except IndexError:
+        righ = np.zeros_like(xbuffer[..., -rm:])
+    try:
+        xbuffer[..., :lm] = left
+    except ValueError:
+        xbuffer[..., :lm] = _right90(left)
+
+    try:
+        xbuffer[..., -rm:] = righ
+    except ValueError:
+        xbuffer[..., -rm:] = _right90(righ)
+
+    return xbuffer
+
+
+def buffer_y_withface(s, face, lm, rm, tp):
+    """Create buffer zone in x direction for one face of an array.
+
+    Parameters
+    ----------
+    s: numpy.ndarray
+        the center field, the last dimension being X,
+        the second dimension being face.
+    face: int
+        which face to create buffer for.
+    lm: int
+        the size of the margin to the bottom.
+    rm: int
+        the size of the margin to the top.
+    tp: seaduck.Topology
+        the topology object of the
+    """
+    shape = list(s.shape)
+    shape.pop(1)
+    shape[-2] += lm + rm
+    ybuffer = np.zeros(shape)
+    ybuffer[..., lm:-rm, :] = s[:, face]
+    try:
+        fc1, iy1, ix1 = tp.ind_moves((face, 0, tp.ixmax), [1 for i in range(lm)])
+        fc2, iy2, ix2 = tp.ind_moves((face, 0, 0), [1])
+        left = _slice_corner(s, fc1, iy1, iy2, ix1, ix2)
+    except IndexError:
+        left = np.zeros_like(ybuffer[..., :lm, :])
+
+    try:
+        fc3, iy3, ix3 = tp.ind_moves((face, tp.iymax, tp.ixmax), [0 for i in range(rm)])
+        fc4, iy4, ix4 = tp.ind_moves((face, tp.iymax, 0), [0])
+        righ = _slice_corner(s, fc3, iy3, iy4, ix3, ix4)
+    except IndexError:
+        righ = np.zeros_like(ybuffer[..., -rm:, :])
+    try:
+        ybuffer[..., :lm, :] = left
+    except ValueError:
+        ybuffer[..., :lm, :] = _left90(left)
+
+    try:
+        ybuffer[..., -rm:, :] = righ
+    except ValueError:
+        ybuffer[..., -rm:, :] = _left90(righ)
+    return ybuffer
+
+
+def third_order_upwind_z(s, w):
+    """Get interpolated salinity in the vertical.
+
+    for more info, see
+    https://mitgcm.readthedocs.io/en/latest/algorithm/adv-schemes.html#third-order-upwind-bias-advection
+
+    Parameters
+    ----------
+    s: numpy.ndarray
+        the tracer field, the first dimension being Z
+    w: numpy.ndarray
+        the vertical velocity field, the first dimension being Zl
+    """
+    w[0] = 0
+    deltas = np.nan_to_num(s[1:] - s[:-1], 0)
+    Rj = np.zeros_like(s)
+    Rj[1:] = deltas
+    Rjp = np.roll(Rj, -1, axis=0)
+    Rjm = np.roll(Rj, 1, axis=0)
+    Rjjp = Rjp - Rj
+    Rjjm = Rj - Rjm
+
+    ssum = np.zeros_like(s)
+    ssum[1:] = s[1:] + s[:-1]
+    sz = 0.5 * (ssum - 1 / 6 * (Rjjp + Rjjm) - np.sign(w) * 1 / 6 * (Rjjp - Rjjm))
+    return sz
+
+
+def third_order_DST_x(xbuffer, u_cfl):
+    """Get interpolated salinity in X.
+
+    for more info, see
+    https://mitgcm.readthedocs.io/en/latest/algorithm/adv-schemes.html#third-order-direct-space-time
+
+    Parameters
+    ----------
+    xbuffer: numpy.ndarray
+        the tracer field with buffer zone added (lm=2,rm=1), the last dimension being X
+    u_cfl: numpy.ndarray
+        the u velocity normalized by grid size in x, in s^-1.
+    """
+    lm = 2
+    deltas = np.nan_to_num(np.diff(xbuffer, axis=-1), 0)
+    Rjp = deltas[..., lm:]
+    Rj = deltas[..., lm - 1 : -1]
+    Rjm = deltas[..., lm - 2 : -2]
+
+    d0 = (2.0 - u_cfl) * (1.0 - u_cfl) / 6
+    d1 = (1.0 - u_cfl * u_cfl) / 6
+
+    sx = 0.5 * (
+        (1 + np.sign(u_cfl)) * (xbuffer[..., lm - 1 : -2] + (d0 * Rj + d1 * Rjm))
+        + (1 - np.sign(u_cfl)) * (xbuffer[..., lm:-1] - (d0 * Rj + d1 * Rjp))
+    )
+    return sx
+
+
+def third_order_DST_y(ybuffer, u_cfl):
+    """Get interpolated salinity in Y.
+
+    for more info, see
+    https://mitgcm.readthedocs.io/en/latest/algorithm/adv-schemes.html#third-order-direct-space-time
+
+    Parameters
+    ----------
+    ybuffer: numpy.ndarray
+        the tracer field with buffer zone added (lm=2,rm=1), the second last dimension being Y
+    u_cfl: numpy.ndarray
+        the v velocity normalized by grid size in y, in s^-1.
+    """
+    lm = 2
+    deltas = np.nan_to_num(np.diff(ybuffer, axis=-2), 0)
+    Rjp = deltas[..., lm:, :]
+    Rj = deltas[..., lm - 1 : -1, :]
+    Rjm = deltas[..., lm - 2 : -2, :]
+
+    d0 = (2.0 - u_cfl) * (1.0 - u_cfl) / 6
+    d1 = (1.0 - u_cfl * u_cfl) / 6
+
+    sy = 0.5 * (
+        (1 + np.sign(u_cfl)) * (ybuffer[..., lm - 1 : -2, :] + (d0 * Rj + d1 * Rjm))
+        + (1 - np.sign(u_cfl)) * (ybuffer[..., lm:-1, :] - (d0 * Rj + d1 * Rjp))
+    )
+    return sy
diff --git a/seaduck/lagrangian.py b/seaduck/lagrangian.py
index eafc9166..353cfe93 100644
--- a/seaduck/lagrangian.py
+++ b/seaduck/lagrangian.py
@@ -6,136 +6,22 @@
 
 from seaduck.eulerian import Position
 from seaduck.kernel_weight import KnW
+from seaduck.lagrangian_budget import store_lists
 from seaduck.ocedata import RelCoord
-from seaduck.runtime_conf import compileable
-from seaduck.utils import find_rel, find_rx_ry_oceanparcel, rel2latlon, to_180
+from seaduck.utils import (
+    _stationary,
+    _stationary_time,
+    _time2wall,
+    _which_early,
+    find_rel,
+    find_rx_ry_oceanparcel,
+    rel2latlon,
+    to_180,
+)
 
 DEG2M = 6271e3 * np.pi / 180
 
 
-@compileable
-def _increment(t, u, du):
-    """Find how far it will travel in duration t.
-
-    For a one dimensional particle with speed u and speed derivative du,
-    find how far it will travel in duration t.
-
-    Parameters
-    ----------
-    t: float, numpy.ndarray
-        The time duration
-    u: float, numpy.ndarray
-        The velocity defined at the starting point.
-    du: float, numpy.ndarray
-        The velocity gradient. Assumed to be constant.
-    """
-    incr = u / du * (np.exp(du * t) - 1)
-    no_gradient = np.abs(du) < 1e-12
-    incr[no_gradient] = (u * t)[no_gradient]
-    return incr
-
-
-def _stationary(t, u, du, x0):
-    """Find the final position after time t.
-
-    For a one dimensional Particle with speed u and speed derivative du
-    starting at x0, find the final position after time t.
-    "Stationary" means that we are assuming there is no time dependency.
-
-    Parameters
-    ----------
-    t: float, numpy.ndarray
-        The time duration
-    u: float, numpy.ndarray
-        The velocity defined at the starting point.
-    du: float, numpy.ndarray
-        The velocity gradient. Assumed to be constant.
-    x0: float, numpy.ndarray
-        The starting position.
-    """
-    incr = _increment(t, u, du)
-    return incr + x0
-
-
-@compileable
-def _stationary_time(u, du, x0):
-    """Find the amount of time to leave the cell.
-
-    Find the amount of time it needs for a Particle to hit x = -0.5 and 0.5.
-    The time could be negative.
-
-    Parameters
-    ----------
-    u: numpy.ndarray
-        The velocity defined at the starting point.
-    du: numpy.ndarray
-        The velocity gradient. Assumed to be constant.
-    x0: numpy.ndarray
-        The starting position.
-
-    Returns
-    -------
-    tl: numpy.ndarray
-        The time it takes to hit -0.5.
-    tr: numpy.ndarray
-        The time it takes to hit 0.5
-    """
-    tl = np.log(1 - du / u * (0.5 + x0)) / du
-    tr = np.log(1 + du / u * (0.5 - x0)) / du
-    no_gradient = np.abs(du) < 1e-12
-    tl[no_gradient] = (-x0[no_gradient] - 0.5) / u[no_gradient]
-    tr[no_gradient] = (0.5 - x0[no_gradient]) / u[no_gradient]
-    return tl, tr
-
-
-@compileable
-def _uleftright_from_udu(u, du, x0):
-    """Calculate the velocity at -0.5 and 0.5."""
-    u_left = u - (x0 + 0.5) * du
-    u_right = u + (0.5 - x0) * du
-    return u_left, u_right
-
-
-def _time2wall(pos_list, u_list, du_list, tf):
-    """Apply stationary_time three times for all three dimensions."""
-    ts = []
-    for i in range(3):
-        tl, tr = _stationary_time(u_list[i], du_list[i], pos_list[i])
-        ul, ur = _uleftright_from_udu(u_list[i], du_list[i], pos_list[i])
-        sign = np.sign(tf)
-        cannot_left = -ul * sign <= 1e-12  # aroung 30000 years
-        tl[cannot_left] = -sign[cannot_left]
-        cannot_right = ur * sign <= 1e-12
-        tr[cannot_right] = -sign[cannot_right]
-        ts.append(tl)
-        ts.append(tr)
-    return ts
-
-
-def _which_early(tf, ts):
-    """Find out which event happens first.
-
-    We are trying to integrate the Particle to time tf.
-    The first event is either:
-    1. tf is reached before reaching a wall
-    2. ts[i] is reached, and a Particle hit a wall. ts[i]*tf>0.
-
-    Parameters
-    ----------
-    tf: float, numpy.ndarray
-        The final time
-    ts: list
-        The list of events calculated using _time2wall
-    """
-    ts.append(np.ones(len(ts[0])) * tf)  # float or array both ok
-    t_directed = np.array(ts) * np.sign(tf)
-    t_directed[np.isnan(t_directed)] = np.inf
-    t_directed[t_directed < 0] = np.inf
-    tend = t_directed.argmin(axis=0)
-    t_event = np.array([ts[te][i] for i, te in enumerate(tend)])
-    return tend, t_event
-
-
 uvkernel = np.array([[0, 0], [1, 0], [0, 1]])
 ukernel = np.array([[0, 0], [1, 0]])
 vkernel = np.array([[0, 0], [0, 1]])
@@ -255,12 +141,7 @@ def __init__(
             try:
                 self.ocedata["Vol"]
             except KeyError:
-                if self.ocedata.readiness["Zl"]:
-                    self.ocedata["Vol"] = np.array(
-                        self.ocedata._ds["drF"] * self.ocedata._ds["rA"]
-                    )
-                else:
-                    self.ocedata["Vol"] = np.array(self.ocedata._ds["rA"])
+                self.ocedata._add_missing_vol(as_numpy=True)
 
         # whether or not setting the w at the surface
         # just to prevent particles taking off
@@ -912,7 +793,12 @@ def to_next_stop(self, t_stop):
             self.it[~before_first] += 1
 
     def to_list_of_time(
-        self, normal_stops, update_stops="default", return_in_between=True
+        self,
+        normal_stops,
+        update_stops="default",
+        return_in_between=True,
+        dump_filename=False,
+        store_kwarg={},
     ):
         """Integrate the particles to a list of time.
 
@@ -944,6 +830,8 @@ def to_list_of_time(
         """
         t_min = np.minimum(np.min(normal_stops), self.t[0])
         t_max = np.maximum(np.max(normal_stops), self.t[0])
+        if not self.save_raw and dump_filename:
+            raise ValueError("saving to file only available if save_raw = True")
 
         if "time" not in self.ocedata[self.uname].dims:
             pass
@@ -978,7 +866,9 @@ def to_list_of_time(
         self.get_u_du()
         to_return = []
         for i, tl in enumerate(stops):
-            logging.info(np.datetime64(round(tl), "s"))
+            timestr = str(np.datetime64(round(tl), "s"))
+            # logging.info(timestr)
+            print(timestr)
             if self.save_raw:
                 # save the very start of everything.
                 self.note_taking(stamp=15)
@@ -992,9 +882,16 @@ def to_list_of_time(
                     self.update_uvw_array()
                     self.get_u_du()
                 if return_in_between:
-                    to_return.append(self.deepcopy())
+                    if dump_filename:
+                        store_lists(self, dump_filename + timestr, **store_kwarg)
+                    else:
+                        to_return.append(self.deepcopy())
             else:
-                to_return.append(self.deepcopy())
+                if dump_filename:
+                    store_lists(self, dump_filename + timestr, **store_kwarg)
+                else:
+                    to_return.append(self.deepcopy())
             if self.save_raw:
                 self.empty_lists()
-        return stops, to_return
+        if not dump_filename:
+            return stops, to_return
diff --git a/seaduck/lagrangian_budget.py b/seaduck/lagrangian_budget.py
new file mode 100644
index 00000000..d2a92eab
--- /dev/null
+++ b/seaduck/lagrangian_budget.py
@@ -0,0 +1,435 @@
+import copy
+
+import numpy as np
+import xarray as xr
+
+from seaduck.eulerian import Position
+from seaduck.runtime_conf import compileable
+from seaduck.utils import (
+    _time2wall,
+    _uleftright_from_udu,
+    _which_early,
+    parallelpointinpolygon,
+)
+
+try:  # pragma: no cover
+    import zarr
+except ImportError:  # pragma: no cover
+    pass
+
+UV_DIC = {"U": 0, "V": 1}
+MOVE_DIC = {
+    0: np.array([0, 0]),  # delta_y,delta_x
+    1: np.array([0, 1]),
+    2: np.array([0, 0]),
+    3: np.array([1, 0]),
+    4: np.array([0, 0]),
+    5: np.array([0, 0]),
+}
+
+
+def read_from_ds(particle_ds, oce):
+    temp = Position.__new__(Position)
+    temp.ocedata = oce
+    temp.tp = temp.ocedata.tp
+
+    # it = np.array(particle_ds.it)
+    izl = np.array(particle_ds.iz)
+    fc = np.array(particle_ds.fc)
+    iy = np.array(particle_ds.iy)
+    ix = np.array(particle_ds.ix)
+    rzl = np.array(particle_ds.rz)
+    ry = np.array(particle_ds.ry)
+    rx = np.array(particle_ds.rx)
+
+    # temp.it  = it .astype(int)
+    temp.izl_lin = izl.astype(int)
+    temp.iz = (izl - 1).astype(int)
+    temp.face = fc.astype(int)
+    temp.iy = iy.astype(int)
+    temp.ix = ix.astype(int)
+    temp.rzl_lin = rzl
+    temp.ry = ry
+    temp.rx = rx
+
+    temp.N = len(temp.ix)
+
+    uu = np.array(particle_ds.uu)
+    vv = np.array(particle_ds.vv)
+    ww = np.array(particle_ds.ww)
+    du = np.array(particle_ds.du)
+    dv = np.array(particle_ds.dv)
+    dw = np.array(particle_ds.dw)
+
+    temp.u = uu
+    temp.v = vv
+    temp.w = ww
+    temp.du = du
+    temp.dv = dv
+    temp.dw = dw
+
+    temp.lon = np.array(particle_ds.xx)
+    temp.lat = np.array(particle_ds.yy)
+    temp.dep = np.array(particle_ds.zz)
+    temp.t = np.array(particle_ds.tt)
+
+    temp.vs = np.array(particle_ds.vs)
+
+    temp.shapes = list(particle_ds.shapes)
+
+    return temp
+
+
+def which_wall(pt):
+    distl = np.abs(0.5 + pt.rx)
+    distr = np.abs(0.5 - pt.rx)
+    distd = np.abs(0.5 + pt.ry)
+    distt = np.abs(0.5 - pt.ry)
+    distdeep = np.abs(pt.rzl_lin)
+    distshal = np.abs(1 - pt.rzl_lin)
+
+    distance = np.vstack([distl, distr, distd, distt, distdeep, distshal])
+    return np.argmin(distance, axis=0)
+
+
+def pseudo_motion(pt):
+    us = [pt.u, pt.v, pt.w]
+    dus = [pt.du, pt.dv, pt.dw]
+    xs = [pt.rx, pt.ry, pt.rzl_lin - 1 / 2]
+    ts_forward = _time2wall(xs, us, dus, 1e80 * np.ones(pt.N))
+    ts_backward = _time2wall(xs, us, dus, -1e80 * np.ones(pt.N))
+    tendf, tf = _which_early(1e80, ts_forward)
+    tendb, tb = _which_early(-1e80, ts_backward)
+
+    return tendf, tf, tendb, tb
+
+
+@compileable
+def fast_cumsum(shapes):
+    return np.cumsum(shapes)
+
+
+def first_last_neither(shapes, return_neither=True):
+    acc = fast_cumsum(shapes)
+    last = acc - 1
+    first = np.roll(acc, 1)
+    first[0] = 0
+    if not return_neither:
+        return first, last
+    else:
+        neither = np.array(
+            [
+                first[i] + j
+                for i, length in enumerate(shapes)
+                for j in range(1, length - 1)
+            ]
+        )
+        return first, last, neither
+
+
+def pt_ulist(pt):
+    ul, ur = _uleftright_from_udu(pt.u, pt.du, pt.rx)
+    vl, vr = _uleftright_from_udu(pt.v, pt.dv, pt.ry)
+    wl, wr = _uleftright_from_udu(pt.w, pt.dw, pt.rzl_lin - 0.5)
+    return np.array([ul, ur, vl, vr, wl, wr])
+
+
+def residence_time(pt):
+    u_list = pt_ulist(pt)
+    return 2 / np.sum(np.abs(u_list), axis=0)
+
+
+def tres_update(tres0, temp, first, last, fraction_first, fraction_last):
+    fracs = np.ones(temp.N)
+    # fracs_a = np.ones(temp.N)
+    fracs[first + 1] = fraction_first
+    # fracs_a[first+1] = fraction_first
+    fracs[last] -= fraction_last
+    tres = tres0 * fracs
+
+    tres[temp.vs > 6] = 0.0
+    # mask = np.logical_and(tres==0, temp.vs<7)
+    # assert (tres[temp.vs<7]>0).all(), (tres0[mask],fracs[mask], fracs_a[mask], np.where(mask))
+
+    return tres
+
+
+def tres_fraction(temp, first, last, fraction_first, fraction_last):
+    tres0 = residence_time(temp)
+    tres = tres_update(tres0, temp, first, last, fraction_first, fraction_last)
+    return tres
+
+
+def ind_tend_uv(ind, tp):
+    """Return the index of the velocity node.
+
+    The node is not of the same index as the center point.
+    iw determines where to return iw
+    """
+    iw, face, iy, ix = ind
+    if iw == 1:
+        new_wall, new_ind = tp._ind_tend_V((face, iy, ix), 0)
+    elif iw == 0:
+        new_wall, new_ind = tp._ind_tend_U((face, iy, ix), 3)
+    else:
+        raise ValueError("illegal iw")
+    niw = UV_DIC[new_wall]
+    nfc, niy, nix = new_ind
+    return niw, nfc, niy, nix
+
+
+def deepcopy_inds(temp):
+    iz = copy.deepcopy(temp.izl_lin)
+    iy = copy.deepcopy(temp.iy)
+    ix = copy.deepcopy(temp.ix)
+    face = copy.deepcopy(temp.face)
+    # assert (iz>=1).all(),iz
+    return iz, face, iy, ix
+
+
+def wall_index(inds, iwall, tp):
+    iw = iwall // 2
+    iz, face, iy, ix = copy.deepcopy(inds)
+    # assert (iz>=1).all(),iz
+
+    ind = copy.deepcopy(np.array([face, iy, ix]))
+    old_ind = copy.deepcopy(ind)
+    naive_move = np.array([MOVE_DIC[i] for i in iwall], dtype=int).T
+    iy += naive_move[0]
+    ix += naive_move[1]
+    ind = np.array([face, iy, ix])
+    illegal = tp.check_illegal(ind, cuvwg="C")
+    redo = np.array(np.where(illegal)).T
+    for num, loc in enumerate(redo):
+        j = loc[0]
+        ind = (iw[j],) + tuple(old_ind[:, j])
+        new_ind = ind_tend_uv(ind, tp)
+        iw[j], face[j], iy[j], ix[j] = new_ind
+
+    iz[iwall == 4] += 1
+    iz -= 1
+    return np.array([iw, iz, face, iy, ix]).astype(int)
+
+
+def redo_index(pt):
+    # assert (pt.izl_lin>=1).all()
+    inds = deepcopy_inds(pt)
+    iz, face, iy, ix = inds
+    # assert (iz>=1).all(),iz
+    tendf, tf, tendb, tb = pseudo_motion(pt)
+
+    funderflow = np.where(tendf == 6)
+    bunderflow = np.where(tendb == 6)
+    tendf[funderflow] = 0
+    tendb[bunderflow] = 0
+    vf = wall_index(inds, tendf, pt.ocedata.tp)
+    vb = wall_index(inds, tendb, pt.ocedata.tp)
+    # vf[:,funderflow] = vb[:,funderflow]
+    # vb[:,bunderflow] = vf[:,bunderflow]
+    tim = tf - tb
+    frac = -tb / tim
+    assert (~np.isnan(tim)).any(), [
+        i[np.isnan(tim)] for i in [pt.rx, pt.ry, pt.rzl_lin - 1 / 2]
+    ]
+    assert (tim != 0).all(), [i[tim == 0] for i in [pt.rx, pt.ry, pt.rzl_lin - 1 / 2]]
+    at_corner = np.where(tim == 0)
+    frac[at_corner] = 1
+    return vf, vb, frac
+
+
+def find_ind_frac_tres(neo, oce, region_names=False, region_polys=None):
+    temp = read_from_ds(neo, oce)
+    temp.shapes = list(temp.shapes)
+    if region_names:
+        masks = []
+        for reg in region_polys:
+            mask = parallelpointinpolygon(temp.lon, temp.lat, reg)
+            # mask = np.where(mask)[0]
+            masks.append(mask)
+    first, last, neither = first_last_neither(np.array(temp.shapes))
+
+    ind1 = np.zeros((5, temp.N), "int16")
+    ind2 = np.ones((5, temp.N), "int16")
+    frac = np.ones(temp.N)
+
+    # ind1[:, wrong_ind] = lookup[:, lookup_ind]
+
+    neithers = temp.subset(neither)
+    neither_inds = deepcopy_inds(neithers)
+    iwalls = which_wall(neithers)
+    ind1[:, neither] = wall_index(neither_inds, iwalls, temp.ocedata.tp)
+
+    firsts = temp.subset(first)
+    lasts = temp.subset(last)
+    ind1[:, first], ind2[:, first], frac[first] = redo_index(firsts)
+    ind1[:, last], ind2[:, last], frac[last] = redo_index(lasts)
+
+    tres = tres_fraction(temp, first, last, frac[first], frac[last])
+    if region_names:
+        return ind1, ind2, frac, masks, tres, last, first
+    else:
+        return ind1, ind2, frac, tres, last, first
+
+
+def flatten(lstoflst, shapes=None):
+    if shapes is None:
+        shapes = [len(i) for i in lstoflst]
+    suffix = np.cumsum(shapes)
+    thething = np.zeros(suffix[-1])
+    thething[: suffix[0]] = lstoflst[0]
+    for i in range(1, len(lstoflst)):
+        thething[suffix[i - 1] : suffix[i]] = lstoflst[i]
+    return thething
+
+
+def particle2xarray(p):
+    shapes = [len(i) for i in p.xxlist]
+    # it = flatten(p.itlist,shapes = shapes)
+    fc = flatten(p.fclist, shapes=shapes)
+    iy = flatten(p.iylist, shapes=shapes)
+    iz = flatten(p.izlist, shapes=shapes)
+    ix = flatten(p.ixlist, shapes=shapes)
+    rx = flatten(p.rxlist, shapes=shapes)
+    ry = flatten(p.rylist, shapes=shapes)
+    rz = flatten(p.rzlist, shapes=shapes)
+    tt = flatten(p.ttlist, shapes=shapes)
+    uu = flatten(p.uulist, shapes=shapes)
+    vv = flatten(p.vvlist, shapes=shapes)
+    ww = flatten(p.wwlist, shapes=shapes)
+    du = flatten(p.dulist, shapes=shapes)
+    dv = flatten(p.dvlist, shapes=shapes)
+    dw = flatten(p.dwlist, shapes=shapes)
+    xx = flatten(p.xxlist, shapes=shapes)
+    yy = flatten(p.yylist, shapes=shapes)
+    zz = flatten(p.zzlist, shapes=shapes)
+    vs = flatten(p.vslist, shapes=shapes)
+
+    ds = xr.Dataset(
+        coords=dict(shapes=(["shapes"], shapes), nprof=(["nprof"], np.arange(len(xx)))),
+        data_vars=dict(
+            # it = (['nprof'],it),
+            fc=(["nprof"], fc),
+            iy=(["nprof"], iy),
+            iz=(["nprof"], iz),
+            ix=(["nprof"], ix),
+            rx=(["nprof"], rx),
+            ry=(["nprof"], ry),
+            rz=(["nprof"], rz),
+            tt=(["nprof"], tt),
+            uu=(["nprof"], uu),
+            vv=(["nprof"], vv),
+            ww=(["nprof"], ww),
+            du=(["nprof"], du),
+            dv=(["nprof"], dv),
+            dw=(["nprof"], dw),
+            xx=(["nprof"], xx),
+            yy=(["nprof"], yy),
+            zz=(["nprof"], zz),
+            vs=(["nprof"], vs),
+        ),
+    )
+    return ds
+
+
+def dump_to_zarr(neo, oce, filename, region_names=False, region_polys=None):
+    if region_names:
+        (ind1, ind2, frac, masks, tres, last, first) = find_ind_frac_tres(
+            neo, oce, region_names=region_names, region_polys=region_polys
+        )
+    else:
+        ind1, ind2, frac, tres, last, first = find_ind_frac_tres(neo, oce)
+
+    neo["five"] = xr.DataArray(["iw", "iz", "face", "iy", "ix"], dims="five")
+    if region_names:
+        for ir, reg in enumerate(region_names):
+            neo[reg] = xr.DataArray(masks[ir].astype(bool), dims="nprof")
+        # neo['gulf'] = xr.DataArray(gulf_ind.astype(bool),dims = 'nprof')
+        # neo['labr'] = xr.DataArray(labr_ind.astype(bool),dims = 'nprof')
+        # neo['gdbk'] = xr.DataArray(gdbk_ind.astype(bool),dims = 'nprof')
+        # neo['nace'] = xr.DataArray(nace_ind.astype(bool),dims = 'nprof')
+        # neo['egrl'] = xr.DataArray(egrl_ind.astype(bool),dims = 'nprof')
+        # region_ind = np.concatenate([gulf_ind,labr_ind,gdbk_ind,nace_ind,egrl_ind]).astype('int32')
+        # neo.attrs['region_shape'] = [len(i) for i in [gulf_ind,labr_ind,gdbk_ind,nace_ind,egrl_ind]]
+        # if len(region_ind)>0:
+        #     neo = neo.assign_coords(region_ind = xr.DataArray(region_ind,dims = 'region_ind'))
+
+    neo["ind1"] = xr.DataArray(ind1.astype("int16"), dims=["five", "nprof"])
+    neo["ind2"] = xr.DataArray(ind2.astype("int16"), dims=["five", "nprof"])
+    neo["frac"] = xr.DataArray(frac, dims="nprof")
+    neo["tres"] = xr.DataArray(tres, dims="nprof")
+    # neo['last'] = xr.DataArray(last.astype('int64'), dims = 'shapes')
+    # neo['first'] = xr.DataArray(first.astype('int64'), dims = 'shapes')
+
+    neo["face"] = neo["fc"].astype("int16")
+    neo["ix"] = neo["ix"].astype("int16")
+    neo["iy"] = neo["iy"].astype("int16")
+    neo["iz"] = neo["iz"].astype("int16")
+    neo["vs"] = neo["vs"].astype("int16")
+
+    neo = neo.drop_vars(["rx", "ry", "rz", "uu", "vv", "ww", "du", "dv", "dw", "fc"])
+
+    neo.to_zarr(filename, mode="w")
+    zarr.consolidate_metadata(filename)
+
+
+def store_lists(pt, name, region_names=False, region_polys=None):
+    neo = particle2xarray(pt)
+    dump_to_zarr(
+        neo, pt.ocedata, name, region_names=region_names, region_polys=region_polys
+    )
+
+
+def prefetch_scalar(ds_slc, scalar_names):
+    prefetch = {}
+    for var in scalar_names:
+        # print(var, end = ' ')
+        prefetch[var] = np.array(ds_slc[var])
+    return prefetch
+
+
+def prefetch_vector(ds_slc, xname="sxprime", yname="syprime", zname="szprime"):
+    return np.array(tuple(np.array(ds_slc[i]) for i in [xname, yname, zname]))
+
+
+def read_prefetched_scalar(ind, scalar_names, prefetch):
+    res = {}
+    for var in scalar_names:
+        res[var] = prefetch[var][ind]
+    return res
+
+
+def lhs_contribution(t, scalar_dic, last, lhs_name="lhs"):
+    deltat = np.nan_to_num(np.diff(t))
+    deltat[last[:-1]] = 0
+    lhs_scalar = scalar_dic[lhs_name][:-1]
+    correction = deltat * lhs_scalar
+    return correction
+
+
+def contr_p_relaxed(deltas, tres, step_dic, termlist, p=1):
+    nds = len(deltas)
+    # if len(wrong_ind)>0:
+    #     if wrong_ind[-1] == len(deltas):
+    #         wrong_ind = wrong_ind[:-1]
+    dic = {"error": np.zeros(nds)}
+    # dic['error'][wrong_ind] = deltas[wrong_ind]
+    # deltas[wrong_ind] = 0
+    # tres[wrong_ind] = 0
+
+    deno = np.zeros(nds)
+    sums = np.zeros(nds)
+    for var in termlist:
+        deno += step_dic[var][:-1] ** (p + 1)
+        sums += step_dic[var][:-1]
+    disparity = deltas - sums * tres
+    total = np.zeros(nds)
+    # dic['quality'] = np.nan_to_num(np.abs((disparity/tres)**(p+1)/deno))
+    # mask = (dic['quality']<=1).astype(int)
+    for var in termlist:
+        ratio = step_dic[var][:-1] ** (p + 1) / deno
+        dic[var] = step_dic[var][:-1] * tres + ratio * disparity
+        total += dic[var]
+    final_correction = deltas - total
+    assert np.allclose(final_correction, 0)
+    dic["error"] += final_correction
+    return dic
diff --git a/seaduck/ocedata.py b/seaduck/ocedata.py
index e0797bf7..c8752001 100644
--- a/seaduck/ocedata.py
+++ b/seaduck/ocedata.py
@@ -279,6 +279,17 @@ def _add_missing_cs_sn(self):
             self["CS"] = cs
             self["SN"] = sn
 
+    def _add_missing_vol(self, as_numpy=False):
+        if self.readiness["Zl"]:
+            vol = self._ds["drF"] * self._ds["rA"]
+        else:
+            vol = self._ds["rA"]
+
+        if as_numpy:
+            self["Vol"] = np.array(vol)
+        else:
+            self["Vol"] = vol
+
     def _hgrid2array(self):
         """Extract the horizontal grid data into numpy arrays.
 
diff --git a/seaduck/runtime_conf.py b/seaduck/runtime_conf.py
index e6cc706d..05cd9770 100644
--- a/seaduck/runtime_conf.py
+++ b/seaduck/runtime_conf.py
@@ -17,3 +17,11 @@ def compileable(func):  # pragma: no cover
         return njit(func)
     else:
         return func
+
+
+def compileable_parallel(func):  # pragma: no cover
+    """Decorate function to compile them (parallel) using numba when available."""
+    if rcParam["compilable"]:
+        return njit(func, parallel=True)
+    else:
+        return func
diff --git a/seaduck/utils.py b/seaduck/utils.py
index 61e70e28..682de80a 100644
--- a/seaduck/utils.py
+++ b/seaduck/utils.py
@@ -7,7 +7,12 @@
 import xarray as xr
 from scipy import spatial
 
-from seaduck.runtime_conf import compileable
+from seaduck.runtime_conf import compileable, compileable_parallel
+
+try:
+    import numba
+except ImportError:
+    pass
 
 try:
     import pooch
@@ -90,6 +95,9 @@ def process_ecco(ds):
         np.stack([etan_3d * 1.5, etan_3d * 2.5], axis=0),
         dims=("time_midp", "face", "Y", "X"),
     )
+    ds["utrans"] = ds["UVELMASS"] * ds["drF"] * ds["dyG"]
+    ds["vtrans"] = ds["VVELMASS"] * ds["drF"] * ds["dxG"]
+    ds["wtrans"] = ds["WVELMASS"] * ds["rA"]
     return ds
 
 
@@ -713,3 +721,157 @@ def easy_3d_cube(lon, lat, dep, tim, print_total_number=False):
     if print_total_number:
         print(f"A total of {len(x)} positions are defined.")
     return x, y, z, t
+
+
+@compileable
+def pointinpolygon(x, y, poly):
+    n = len(poly)
+    inside = False
+    p2x = 0.0
+    p2y = 0.0
+    xints = 0.0
+    p1x, p1y = poly[0]
+    for i in numba.prange(n + 1):
+        p2x, p2y = poly[i % n]
+        if y > min(p1y, p2y):
+            if y <= max(p1y, p2y):
+                if x <= max(p1x, p2x):
+                    if p1y != p2y:
+                        xints = (y - p1y) * (p2x - p1x) / (p2y - p1y) + p1x
+                    if p1x == p2x or x <= xints:
+                        inside = not inside
+        p1x, p1y = p2x, p2y
+
+    return inside
+
+
+@compileable_parallel
+def parallelpointinpolygon(xs, ys, poly):
+    D = np.empty(len(xs), dtype=numba.boolean)
+    for i in numba.prange(0, len(D)):
+        D[i] = pointinpolygon(xs[i], ys[i], poly)
+    return D
+
+
+# functions used in Lagrangian models to avoid cyclic import
+@compileable
+def _increment(t, u, du):
+    """Find how far it will travel in duration t.
+
+    For a one dimensional particle with speed u and speed derivative du,
+    find how far it will travel in duration t.
+
+    Parameters
+    ----------
+    t: float, numpy.ndarray
+        The time duration
+    u: float, numpy.ndarray
+        The velocity defined at the starting point.
+    du: float, numpy.ndarray
+        The velocity gradient. Assumed to be constant.
+    """
+    incr = u / du * (np.exp(du * t) - 1)
+    no_gradient = np.abs(du) < 1e-12
+    incr[no_gradient] = (u * t)[no_gradient]
+    return incr
+
+
+def _stationary(t, u, du, x0):
+    """Find the final position after time t.
+
+    For a one dimensional Particle with speed u and speed derivative du
+    starting at x0, find the final position after time t.
+    "Stationary" means that we are assuming there is no time dependency.
+
+    Parameters
+    ----------
+    t: float, numpy.ndarray
+        The time duration
+    u: float, numpy.ndarray
+        The velocity defined at the starting point.
+    du: float, numpy.ndarray
+        The velocity gradient. Assumed to be constant.
+    x0: float, numpy.ndarray
+        The starting position.
+    """
+    incr = _increment(t, u, du)
+    return incr + x0
+
+
+@compileable
+def _stationary_time(u, du, x0):
+    """Find the amount of time to leave the cell.
+
+    Find the amount of time it needs for a Particle to hit x = -0.5 and 0.5.
+    The time could be negative.
+
+    Parameters
+    ----------
+    u: numpy.ndarray
+        The velocity defined at the starting point.
+    du: numpy.ndarray
+        The velocity gradient. Assumed to be constant.
+    x0: numpy.ndarray
+        The starting position.
+
+    Returns
+    -------
+    tl: numpy.ndarray
+        The time it takes to hit -0.5.
+    tr: numpy.ndarray
+        The time it takes to hit 0.5
+    """
+    tl = np.log(1 - du / u * (0.5 + x0)) / du
+    tr = np.log(1 + du / u * (0.5 - x0)) / du
+    no_gradient = np.abs(du) < 1e-12
+    tl[no_gradient] = (-x0[no_gradient] - 0.5) / u[no_gradient]
+    tr[no_gradient] = (0.5 - x0[no_gradient]) / u[no_gradient]
+    return tl, tr
+
+
+@compileable
+def _uleftright_from_udu(u, du, x0):
+    """Calculate the velocity at -0.5 and 0.5."""
+    u_left = u - (x0 + 0.5) * du
+    u_right = u + (0.5 - x0) * du
+    return u_left, u_right
+
+
+def _time2wall(pos_list, u_list, du_list, tf):
+    """Apply stationary_time three times for all three dimensions."""
+    ts = []
+    for i in range(3):
+        tl, tr = _stationary_time(u_list[i], du_list[i], pos_list[i])
+        ul, ur = _uleftright_from_udu(u_list[i], du_list[i], pos_list[i])
+        sign = np.sign(tf)
+        cannot_left = -ul * sign <= 1e-12  # aroung 30000 years
+        tl[cannot_left] = -sign[cannot_left]
+        cannot_right = ur * sign <= 1e-12
+        tr[cannot_right] = -sign[cannot_right]
+        ts.append(tl)
+        ts.append(tr)
+    return ts
+
+
+def _which_early(tf, ts):
+    """Find out which event happens first.
+
+    We are trying to integrate the Particle to time tf.
+    The first event is either:
+    1. tf is reached before reaching a wall
+    2. ts[i] is reached, and a Particle hit a wall. ts[i]*tf>0.
+
+    Parameters
+    ----------
+    tf: float, numpy.ndarray
+        The final time
+    ts: list
+        The list of events calculated using _time2wall
+    """
+    ts.append(np.ones(len(ts[0])) * tf)  # float or array both ok
+    t_directed = np.array(ts) * np.sign(tf)
+    t_directed[np.isnan(t_directed)] = np.inf
+    t_directed[t_directed < 0] = np.inf
+    tend = t_directed.argmin(axis=0)
+    t_event = np.array([ts[te][i] for i, te in enumerate(tend)])
+    return tend, t_event
diff --git a/tests/test_lagrangian_budget.py b/tests/test_lagrangian_budget.py
new file mode 100644
index 00000000..ed1b30f2
--- /dev/null
+++ b/tests/test_lagrangian_budget.py
@@ -0,0 +1,133 @@
+import numpy as np
+import pytest
+
+import seaduck as sd
+from seaduck import utils
+from seaduck.lagrangian_budget import (
+    find_ind_frac_tres,
+    first_last_neither,
+    flatten,
+    ind_tend_uv,
+    particle2xarray,
+    redo_index,
+    store_lists,
+)
+
+
+@pytest.fixture
+def custom_pt():
+    x = np.linspace(-50, -15, 200)
+    y = np.ones_like(x) * 52.0
+    z = np.ones_like(x) * (-9)
+    t = np.ones_like(x)
+    od = sd.OceData(utils.get_dataset("ecco"))
+    pt = sd.Particle(
+        x=x,
+        y=y,
+        z=z,
+        t=t,
+        data=od,
+        uname="utrans",
+        vname="vtrans",
+        wname="wtrans",
+        transport=True,
+        save_raw=True,
+    )
+    pt.to_next_stop(t[0] + 1e7)
+    return pt
+
+
+@pytest.fixture
+def region_info():
+    GULF = np.array(
+        [
+            [-92.5, 25],
+            [-75, 25],
+            [-55, 40],
+            [-72.5, 40],
+        ]
+    )
+    LABR = np.array([[-72, 63.5], [-60, 50], [-49, 53], [-61, 66.5]])
+    GDBK = np.array(
+        [
+            # GULF[2],
+            GULF[3],
+            LABR[1],
+            LABR[2],
+            [-42, 50],
+            [-42, 40],
+        ]
+    )
+    NACE = np.array([GDBK[-2], GDBK[-1], [-21, 47], [-21, 57]])
+    EGRL = np.array([[-22.5, 72], [-44, 63], [-44, 57], [-22.5, 66]])
+
+    return ["gulf", "labr", "gdbk", "nace", "egrl"], [GULF, LABR, GDBK, NACE, EGRL]
+
+
+@pytest.mark.parametrize(
+    "ind,exp",
+    [
+        ((0, 2, 45, 45), (0, 2, 45, 46)),
+        ((1, 2, 45, 45), (1, 2, 46, 45)),
+        ((1, 2, 89, 45), (0, 6, 44, 0)),
+        ((0, 5, 45, 89), (1, 7, 0, 44)),
+    ],
+)
+@pytest.mark.parametrize("od", ["ecco"], indirect=True)
+def test_ind_tend_uv(ind, exp, od):
+    tub = od
+    tp = tub.tp
+    ans = ind_tend_uv(ind, tp)
+    assert exp == ans
+
+
+def test_ind_tend_uv_error():
+    ind = (4, 3, 2, 1)
+    with pytest.raises(ValueError):
+        ind_tend_uv(ind, None)
+
+
+def test_redo_index(custom_pt):
+    vf, vb, frac = redo_index(custom_pt)
+    assert (frac <= 1).all()
+    assert (frac >= 0).all()
+
+
+def test_flatten_list():
+    lst = [[0, 1], [1, 1, 1, 1]]
+    new = flatten(lst)
+    assert (new == np.array([0, 1, 1, 1, 1, 1])).all()
+    assert isinstance(new, np.ndarray)
+
+
+@pytest.mark.parametrize("od", ["ecco"], indirect=True)
+def test_ind_frac_find(custom_pt, od):
+    particle_datasets = particle2xarray(custom_pt)
+    tub = od
+    ind1, ind2, frac, tres, last, first = find_ind_frac_tres(particle_datasets, tub)
+    assert ind1.shape[0] == 5
+    assert (frac != 1).any()
+    assert (tres >= 0).all()
+
+
+def test_store_lists(custom_pt):
+    store_lists(custom_pt, "PleaseIgnore_dump.zarr")
+
+
+def test_store_lists_with_region(custom_pt, region_info):
+    pytest.importorskip("numba")
+    region_names, region_polys = region_info
+    store_lists(
+        custom_pt,
+        "PleaseIgnore_dump.zarr",
+        region_names=region_names,
+        region_polys=region_polys,
+    )
+
+
+def test_first_last_neither():
+    shapes = np.random.randint(2, 5, 10)
+    first, last, neither = first_last_neither(shapes)
+    merged = np.sort(np.concatenate([first, last, neither]))
+    sums = np.sum(shapes)
+    assert np.allclose(np.arange(sums), merged)
diff --git a/tests/test_particle.py b/tests/test_particle.py
index 053b3a0d..e5345774 100644
--- a/tests/test_particle.py
+++ b/tests/test_particle.py
@@ -138,7 +138,7 @@ def test_underflow_u(one_p, seed):
 @pytest.mark.parametrize("seed", [0])
 def test_u_du_uwall_conversion(one_p, seed):
     one_p = random_p(one_p, seed)
-    ul, ur = sd.lagrangian._uleftright_from_udu(one_p.u, one_p.du, one_p.rx)
+    ul, ur = sd.utils._uleftright_from_udu(one_p.u, one_p.du, one_p.rx)
     u, du = u_du_from_uwall(ul, ur, one_p.rx)
     assert np.allclose(u, one_p.u)
     assert np.allclose(du, one_p.du)