diff --git a/xbout/load.py b/xbout/load.py
index 6a8effd2..4f565f8b 100644
--- a/xbout/load.py
+++ b/xbout/load.py
@@ -18,6 +18,11 @@
     _is_dir,
 )
 
+# Override file reading engine.
+# Use: xbout.load.file_engine = "netcdf4" or "h5netcdf"
+file_engine = None
+
+
 _BOUT_GEOMETRY_VARS = [
     "ixseps1",
     "ixseps2",
@@ -584,8 +589,30 @@ def _check_dataset_type(datapath):
 
     filepaths, filetype = _expand_filepaths(datapath)
 
-    ds = xr.open_dataset(filepaths[0], engine=filetype)
-    ds.close()
+    try:
+        ds = xr.open_dataset(filepaths[0], engine=file_engine or filetype)
+        ds.close()
+    except RuntimeError as e:
+        if "H5DSget_num_scales" in str(e):
+            msg = (
+                "\n\nFailed to open dataset due to an HDF5 compatibility error between\n"
+                "h5py and h5netcdf, likely because both were installed with pip.\n"
+                "See: https://github.com/boutproject/xBOUT/issues/329\n"
+                "Also see: https://github.com/HDFGroup/hdf5/issues/6268\n\n"
+                "There are three possible fixes:\n"
+                "  1. Install both from source against a single shared HDF5:\n"
+                "       sudo apt install libhdf5-dev libnetcdf-dev\n"
+                "       pip install --no-binary netCDF4,h5py netCDF4 h5py\n\n"
+                "  2. Install both from your distribution package manager,\n"
+                "     like apt or dnf or install with conda or Spack\n\n"
+                "  3. Switch to the netcdf4 engine:\n"
+                "       import xbout\n"
+                "       xbout.load.file_engine = 'netcdf4'\n"
+                "     netcdf4 may however cause segfaults on file closure\n\n"
+                f"Original error:\n\t{e}"
+            )
+            raise RuntimeError(msg) from e
+        raise
     if "metadata:keep_yboundaries" in ds.attrs:
         # (i)
         return "reload"
@@ -653,7 +680,7 @@ def _auto_open_mfboutdataset(
             concat_dim=concat_dims,
             combine="nested",
             preprocess=_preprocess,
-            engine=filetype,
+            engine=file_engine or filetype,
             chunks=chunks,
             # Only data variables in which the dimension already
             # appears are concatenated.
@@ -790,7 +817,7 @@ def get_nonnegative_scalar(ds, key, default=1, info=True):
                 print(f"{key} not found, setting to {default}")
             if default < 0:
                 raise ValueError(
-                    f"Default for {key} is {val}," f" but negative values are not valid"
+                    f"Default for {key} is {val}, but negative values are not valid"
                 )
             return default
 
@@ -1117,7 +1144,9 @@ def _open_grid(datapath, chunks, keep_xboundaries, keep_yboundaries, mxg=2, **kw
     if _is_path(datapath):
         gridfilepath = Path(datapath)
         grid = xr.open_dataset(
-            gridfilepath, engine=_check_filetype(gridfilepath), **kwargs
+            gridfilepath,
+            engine=(file_engine or _check_filetype(gridfilepath)),
+            **kwargs,
         )
     else:
         grid = datapath