diff --git a/examples/tutorial/blob2d.ipynb b/examples/tutorial/blob2d.ipynb
index 563aa202..ca4db220 100644
--- a/examples/tutorial/blob2d.ipynb
+++ b/examples/tutorial/blob2d.ipynb
@@ -198,7 +198,7 @@
    "outputs": [],
    "source": [
     "# dz is a scalar, so gets stored in 'metadata'\n",
-    "dz = xr.DataArray(ds.metadata['dz']).expand_dims({'z': ds.dims['z']})\n",
+    "dz = xr.DataArray(ds.metadata['dz']).expand_dims({'z': ds.sizes['z']})\n",
     "z = dz.cumsum(dim='z')\n",
     "ds = ds.assign_coords({'z': z})\n",
     "\n",
diff --git a/xbout/geometries.py b/xbout/geometries.py
index 2299f54b..12ef5d20 100644
--- a/xbout/geometries.py
+++ b/xbout/geometries.py
@@ -144,7 +144,7 @@ def apply_geometry(ds, geometry_name, *, coordinates=None, grid=None):
         # 'dx' may not be consistent between different regions (e.g. core and PFR).
         # For some geometries xcoord may have already been created by
         # add_geometry_coords, in which case we do not need this.
-        nx = updated_ds.dims[xcoord]
+        nx = updated_ds.sizes[xcoord]
 
         # can't use commented out version, uncommented one works around xarray bug
         # removing attrs
@@ -181,7 +181,7 @@ def apply_geometry(ds, geometry_name, *, coordinates=None, grid=None):
     if zcoord in updated_ds.dims and zcoord not in updated_ds.coords:
         # Generates a coordinate whose value is 0 on the first grid point, not dz/2, to
         # match how BOUT++ generates fields from input file expressions.
-        nz = updated_ds.dims[zcoord]
+        nz = updated_ds.sizes[zcoord]
 
         # In BOUT++ v5, dz is either a Field2D or Field3D.
         # We can use it as a 1D coordinate if it's a Field3D, _or_ if nz == 1
diff --git a/xbout/load.py b/xbout/load.py
index 3a82a617..851b29cb 100644
--- a/xbout/load.py
+++ b/xbout/load.py
@@ -810,10 +810,10 @@ def get_nonnegative_scalar(ds, key, default=1, info=True):
     mxg = get_nonnegative_scalar(ds, "MXG", default=2, info=info)
     myg = get_nonnegative_scalar(ds, "MYG", default=0, info=info)
     mxsub = get_nonnegative_scalar(
-        ds, "MXSUB", default=ds.dims["x"] - 2 * mxg, info=info
+        ds, "MXSUB", default=ds.sizes["x"] - 2 * mxg, info=info
     )
     mysub = get_nonnegative_scalar(
-        ds, "MYSUB", default=ds.dims["y"] - 2 * myg, info=info
+        ds, "MYSUB", default=ds.sizes["y"] - 2 * myg, info=info
     )
 
     # Check whether this is a single file squashed from the multiple output files of a
@@ -828,8 +828,8 @@ def get_nonnegative_scalar(ds, key, default=1, info=True):
     else:
         # Workaround for older data files
         ny = ds["MYSUB"].values * ds["NYPE"].values
-    nx_file = ds.dims["x"]
-    ny_file = ds.dims["y"]
+    nx_file = ds.sizes["x"]
+    ny_file = ds.sizes["y"]
     is_squashed_doublenull = False
     if nxpe > 1 or nype > 1:
         # if nxpe = nype = 1, was only one process anyway, so no need to check for