diff --git a/.clang-tidy b/.clang-tidy
index a142355fc4..23167e93d0 100644
--- a/.clang-tidy
+++ b/.clang-tidy
@@ -18,6 +18,9 @@ CheckOptions:
 
   - key: readability-qualified-auto.AllowedTypes
     value: 'MPI_Comm'
+
+  - key: misc-include-cleaner.IgnoreHeaders
+    value: 'adios2/.*;bits/.*'
 ---
 
 Disabled checks and reasons:
diff --git a/include/bout/adios_object.hxx b/include/bout/adios_object.hxx
index 07482a0dcb..fcade4f55c 100755
--- a/include/bout/adios_object.hxx
+++ b/include/bout/adios_object.hxx
@@ -16,6 +16,8 @@
 
 #if BOUT_HAS_ADIOS2
 
+#include "bout/boutexception.hxx"
+
 #include <adios2.h>
 #include <memory>
 #include <mpi.h>
@@ -36,14 +38,12 @@ IOPtr GetIOPtr(const std::string IOName);
 class ADIOSStream {
 public:
   adios2::IO io;
-  adios2::Engine engine;
   adios2::Variable<double> vTime;
   adios2::Variable<int> vStep;
   int adiosStep = 0;
-  bool isInStep = false; // true if BeginStep was called and EndStep was not yet called
 
   /** create or return the ADIOSStream based on the target file name */
-  static ADIOSStream& ADIOSGetStream(const std::string& fname);
+  static ADIOSStream& ADIOSGetStream(const std::string& fname, adios2::Mode mode);
 
   ~ADIOSStream();
 
@@ -74,14 +74,63 @@ public:
     return v;
   }
 
+  auto engine() -> adios2::Engine& {
+    if (not engine_) {
+      engine_ = io.Open(fname, file_mode);
+      if (not engine_) {
+        throw BoutException("Could not open ADIOS file '{:s}' for writing", fname);
+      }
+    }
+    return engine_;
+  }
+
+  void beginStep() {
+    if (not isInStep) {
+      engine().BeginStep();
+      isInStep = true;
+      adiosStep = static_cast<int>(engine().CurrentStep());
+    }
+  }
+
+  void endStep() {
+    if (isInStep) {
+      engine().EndStep();
+      isInStep = false;
+    }
+  }
+
+  void finish() {
+    if (engine_) {
+      engine().EndStep();
+      engine().Close();
+    }
+  }
+
 private:
-  ADIOSStream(const std::string fname) : fname(fname){};
+  ADIOSStream(const std::string& fname, adios2::Mode mode)
+      : fname(fname), file_mode(mode) {
+
+    ADIOSPtr adiosp = GetADIOSPtr();
+    std::string ioname = "write_" + fname;
+    try {
+      io = adiosp->AtIO(ioname);
+    } catch (const std::invalid_argument& e) {
+      io = adiosp->DeclareIO(ioname);
+      io.SetEngine("BP5");
+    }
+  };
+
   std::string fname;
+  adios2::Mode file_mode;
+  adios2::Engine engine_;
+
+  /// true if BeginStep was called and EndStep was not yet called
+  bool isInStep = false;
 };
 
 /** Set user parameters for an IO group */
-void ADIOSSetParameters(const std::string& input, const char delimKeyValue,
-                        const char delimItem, adios2::IO& io);
+void ADIOSSetParameters(const std::string& input, char delimKeyValue, char delimItem,
+                        adios2::IO& io);
 
 } // namespace bout
 
diff --git a/src/sys/adios_object.cxx b/src/sys/adios_object.cxx
index 53d8ccaa84..1b3e0a12e2 100644
--- a/src/sys/adios_object.cxx
+++ b/src/sys/adios_object.cxx
@@ -5,8 +5,8 @@
 #include "bout/adios_object.hxx"
 #include "bout/boutexception.hxx"
 
-#include <exception>
-#include <iostream>
+#include <adios2.h>
+
 #include <unordered_map>
 
 namespace bout {
@@ -48,25 +48,25 @@ IOPtr GetIOPtr(const std::string IOName) {
 }
 
 ADIOSStream::~ADIOSStream() {
-  if (engine) {
+  if (engine_) {
     if (isInStep) {
-      engine.EndStep();
+      engine_.EndStep();
       isInStep = false;
     }
-    engine.Close();
+    engine_.Close();
   }
 }
 
-ADIOSStream& ADIOSStream::ADIOSGetStream(const std::string& fname) {
+ADIOSStream& ADIOSStream::ADIOSGetStream(const std::string& fname, adios2::Mode mode) {
   auto it = adiosStreams.find(fname);
   if (it == adiosStreams.end()) {
-    it = adiosStreams.emplace(fname, ADIOSStream(fname)).first;
+    it = adiosStreams.emplace(fname, ADIOSStream(fname, mode)).first;
   }
   return it->second;
 }
 
-void ADIOSSetParameters(const std::string& input, const char delimKeyValue,
-                        const char delimItem, adios2::IO& io) {
+void ADIOSSetParameters(const std::string& input, char delimKeyValue, char delimItem,
+                        adios2::IO& io) {
   auto lf_Trim = [](std::string& input) {
     input.erase(0, input.find_first_not_of(" \n\r\t")); // prefixing spaces
     input.erase(input.find_last_not_of(" \n\r\t") + 1); // suffixing spaces
@@ -76,7 +76,7 @@ void ADIOSSetParameters(const std::string& input, const char delimKeyValue,
   std::string parameter;
   while (std::getline(inputSS, parameter, delimItem)) {
     const size_t position = parameter.find(delimKeyValue);
-    if (position == parameter.npos) {
+    if (position == std::string::npos) {
       throw BoutException("ADIOSSetParameters(): wrong format for IO parameter "
                           + parameter + ", format must be key" + delimKeyValue
                           + "value for each entry");
diff --git a/src/sys/options/options_adios.cxx b/src/sys/options/options_adios.cxx
index 4a41a014eb..2e64423824 100644
--- a/src/sys/options/options_adios.cxx
+++ b/src/sys/options/options_adios.cxx
@@ -1,50 +1,172 @@
 #include "bout/build_defines.hxx"
-#include "bout/traits.hxx"
 
 #if BOUT_HAS_ADIOS2
 
 #include "options_adios.hxx"
-#include "bout/adios_object.hxx"
 
-#include "bout/bout.hxx"
+#include "bout/adios_object.hxx"
+#include "bout/array.hxx"
+#include "bout/assert.hxx"
+#include "bout/bout_types.hxx"
 #include "bout/boutexception.hxx"
+#include "bout/field2d.hxx"
 #include "bout/globals.hxx"
 #include "bout/mesh.hxx"
+#include "bout/options.hxx"
+#include "bout/options_io.hxx"
+#include "bout/output.hxx"
 #include "bout/sys/timer.hxx"
+#include "bout/sys/variant.hxx"
+#include "bout/traits.hxx"
+#include "bout/utils.hxx"
 
-#include "adios2.h"
+#include <adios2.h> // IWYU pragma: keep
+#include <fmt/format.h>
+#include <fmt/ranges.h>
 
+#include <algorithm>
 #include <cstddef>
+#include <cstdint>
 #include <exception>
-#include <tuple>
+#include <iterator>
+#include <stdexcept>
+#include <string>
 #include <utility>
 #include <vector>
 
-namespace bout {
-/// Name of the attribute used to track individual variable's time indices
-constexpr auto current_time_index_name = "current_time_index";
+namespace {
+auto to_int_dims(const adios2::Dims& dims) {
+  std::vector<int> int_dims;
+  int_dims.reserve(dims.size());
+  std::transform(dims.begin(), dims.end(), std::back_inserter(int_dims),
+                 [](auto dim) { return static_cast<int>(dim); });
+  return int_dims;
+}
 
-OptionsADIOS::OptionsADIOS(Options& options) : OptionsIO(options) {
-  if (options["file"].doc("File name. Defaults to <path>/<prefix>.pb").isSet()) {
-    filename = options["file"].as<std::string>();
-  } else {
-    // Both path and prefix must be set
-    filename = fmt::format("{}/{}.bp", options["path"].as<std::string>(),
-                           options["prefix"].as<std::string>());
+// Helper class to construct Adios hyperslices
+struct Selection {
+  // Offset of this processor's data into the global array
+  adios2::Dims start;
+  // The size of the mapped region
+  adios2::Dims count;
+  // Where the actual data starts in data pointer (to exclude ghost cells)
+  adios2::Dims mem_start;
+  // The actual size of data pointer in memory (including ghost cells)
+  adios2::Dims mem_count;
+  // Global shape, including boundaries but not guard cells
+  adios2::Dims shape;
+  // Shape of the local variable to read into
+  std::vector<int> dims;
+
+  // Distributed Field/Array/Matrix/Tensor
+  bool should_set_selection{false};
+
+  Selection(const std::vector<std::string>& dim_names, const std::vector<int>& dim_sizes,
+            const Mesh& mesh) {
+    const auto ndims = dim_names.size();
+    const bool dim0_is_rank = ndims > 0 ? (dim_names[0] == "rank") : false;
+    const bool dim0_is_x = ndims > 0 ? (dim_names[0] == "x") : false;
+    const bool dim1_is_y = ndims > 1 ? (dim_names[1] == "y") : false;
+    const bool dim1_is_z = ndims > 1 ? (dim_names[1] == "z") : false;
+    const bool dim2_is_z = ndims > 2 ? (dim_names[2] == "z") : false;
+
+    should_set_selection = dim0_is_rank or (ndims == 2 and dim0_is_x and dim1_is_y)
+                           or (ndims == 2 and dim0_is_x and dim1_is_z)
+                           or (ndims == 3 and dim0_is_x and dim1_is_y and dim2_is_z);
+
+    if (dim0_is_rank) {
+      const auto ndim_sizes = dim_sizes.size();
+      ASSERT3(ndim_sizes > 1);
+
+      // This is a distributed array, so the local variable is going
+      // to be shape (dim_sizes[1]...) (that is, drop the rank)
+      dims.push_back(dim_sizes[1]);
+      // but we tell adios to read our rank's bit with the full ndims
+      start = {static_cast<std::size_t>(BoutComm::rank()), 0};
+      count = {std::size_t{1}, static_cast<std::size_t>(dim_sizes[0])};
+
+      if (ndim_sizes > 2) {
+        dims.push_back(dim_sizes[2]);
+        start.push_back(0);
+        count.push_back(dim_sizes[1]);
+      }
+      if (ndim_sizes > 3) {
+        dims.push_back(dim_sizes[3]);
+        start.push_back(0);
+        count.push_back(dim_sizes[2]);
+      }
+      mem_count = count;
+      mem_start = start;
+      return;
+    }
+
+    if (ndims > 0) {
+      dims.push_back(dim0_is_x ? mesh.LocalNx : dim_sizes[0]);
+    }
+    if (ndims > 1) {
+      if (dim1_is_y) {
+        dims.push_back(mesh.LocalNy);
+      } else if (dim1_is_z) {
+        dims.push_back(mesh.LocalNz);
+      } else {
+        dims.push_back(dim_sizes[1]);
+      }
+    }
+    if (ndims > 2) {
+      dims.push_back(dim2_is_z ? mesh.LocalNz : dim_sizes[2]);
+    }
+
+    shape.push_back(static_cast<std::size_t>(mesh.GlobalNx));
+    start.push_back(static_cast<std::size_t>(mesh.MapGlobalX));
+    count.push_back(static_cast<std::size_t>(mesh.MapCountX));
+    mem_start.push_back(static_cast<std::size_t>(mesh.MapLocalX));
+    mem_count.push_back(static_cast<std::size_t>(mesh.LocalNx));
+
+    if (dim1_is_y) {
+      shape.push_back(static_cast<std::size_t>(mesh.GlobalNy));
+      start.push_back(static_cast<std::size_t>(mesh.MapGlobalY));
+      count.push_back(static_cast<std::size_t>(mesh.MapCountY));
+      mem_start.push_back(static_cast<std::size_t>(mesh.MapLocalY));
+      mem_count.push_back(static_cast<std::size_t>(mesh.LocalNy));
+    } else if (dim1_is_z) {
+      shape.push_back(static_cast<std::size_t>(mesh.GlobalNz));
+      start.push_back(static_cast<std::size_t>(mesh.MapGlobalZ));
+      count.push_back(static_cast<std::size_t>(mesh.MapCountZ));
+      mem_start.push_back(static_cast<std::size_t>(mesh.MapLocalZ));
+      mem_count.push_back(static_cast<std::size_t>(mesh.LocalNz));
+    }
+
+    if (dim2_is_z) {
+      shape.push_back(static_cast<std::size_t>(mesh.GlobalNz));
+      start.push_back(static_cast<std::size_t>(mesh.MapGlobalZ));
+      count.push_back(static_cast<std::size_t>(mesh.MapCountZ));
+      mem_start.push_back(static_cast<std::size_t>(mesh.MapLocalZ));
+      mem_count.push_back(static_cast<std::size_t>(mesh.LocalNz));
+    }
   }
 
-  file_mode = (options["append"].doc("Append to existing file?").withDefault<bool>(false))
-                  ? FileMode::append
-                  : FileMode::replace;
+  auto selection() const { return adios2::Box<adios2::Dims>{start, count}; }
+  auto memorySelection() const { return adios2::Box<adios2::Dims>{mem_start, mem_count}; }
+};
 
-  singleWriteFile = options["singleWriteFile"].withDefault<bool>(false);
+template <template <class> class T, class U>
+auto read_variable(adios2::IO& io, adios2::Engine& reader, const std::string& name,
+                   const Selection& selection, T<U>& value) {
+  auto variable = io.InquireVariable<U>(name);
+
+  if (selection.should_set_selection) {
+    variable.SetSelection(selection.selection());
+    variable.SetMemorySelection(selection.memorySelection());
+  }
+
+  reader.Get<U>(variable, value.begin(), adios2::Mode::Sync);
+  return Options(value);
 }
 
 template <class T>
 Options readVariable(adios2::Engine& reader, adios2::IO& io, const std::string& name,
                      const std::string& type) {
-  std::vector<T> data;
-  adios2::Variable<T> variable = io.InquireVariable<T>(name);
+  auto variable = io.InquireVariable<T>(name);
 
   using bout::globals::mesh;
 
@@ -55,207 +177,183 @@ Options readVariable(adios2::Engine& reader, adios2::IO& io, const std::string&
   }
 
   if (variable.ShapeID() == adios2::ShapeID::LocalArray) {
-    throw std::invalid_argument(
-        "ADIOS reader did not implement reading local arrays like " + type + " " + name
-        + " in file " + reader.Name());
+    throw BoutException(
+        "ADIOS reader did not implement reading local arrays like `{}` '{}' in file '{}'",
+        type, name, reader.Name());
   }
 
-  if (type != "double" && type != "float") {
-    throw std::invalid_argument(
-        "ADIOS reader did not implement reading arrays that are not double/float type. "
-        "Found "
-        + type + " " + name + " in file " + reader.Name());
+  if (type != "double" && type != "float" && type != "int32_t") {
+    throw BoutException("ADIOS reader did not implement reading arrays that are not "
+                        "`double`/`float`/`int32_t` type. "
+                        "Found `{}` '{}' in file '{}'",
+                        type, name, reader.Name());
   }
 
   if (type == "double" && sizeof(BoutReal) != sizeof(double)) {
-    throw std::invalid_argument(
+    throw BoutException(
         "ADIOS does not allow for implicit type conversions. BoutReal type is "
-        "float but found "
-        + type + " " + name + " in file " + reader.Name());
+        "float but found `{}` '{}' in file '{}'",
+        type, name, reader.Name());
   }
 
   if (type == "float" && sizeof(BoutReal) != sizeof(float)) {
-    throw std::invalid_argument(
+    throw BoutException(
         "ADIOS reader does not allow for implicit type conversions. BoutReal type is "
-        "double but found "
-        + type + " " + name + " in file " + reader.Name());
+        "double but found `{}` '{}' in file '{}'",
+        type, name, reader.Name());
   }
 
-  auto dims = variable.Shape();
-  auto ndims = dims.size();
-  adios2::Variable<BoutReal> variableD = io.InquireVariable<BoutReal>(name);
+  const auto dims = to_int_dims(variable.Shape());
+  auto dims_attr =
+      io.InquireAttribute<std::string>(fmt::format("{}/__xarray_dimensions__", name))
+          .Data();
+
+  const auto selection = Selection(dims_attr, dims, *mesh);
+  const auto ndims = selection.dims.size();
 
   switch (ndims) {
   case 1: {
-    Array<BoutReal> value(static_cast<int>(dims[0]));
-    BoutReal* data = value.begin();
-    reader.Get<BoutReal>(variableD, data, adios2::Mode::Sync);
-    return Options(value);
+    if (type == "float" or type == "double") {
+      Array<BoutReal> value(dims[0]);
+      return read_variable(io, reader, name, selection, value);
+    }
+    if (type == "int32_t") {
+      Array<int> value(dims[0]);
+      return read_variable(io, reader, name, selection, value);
+    }
+    break;
   }
   case 2: {
-    // This could be a Field2D (XY) or FieldPerp (XZ)
-    // Here we look for array sizes, but if Y and Z dimensions are the same size
-    // then it's ambiguous. This method also depends on the global Mesh object.
-    // Some possibilities:
-    // - Add an attribute to specify field type or dimension labels
-    // - Load all the data, and select a region when converting to a Field in Options
-    // - Add a lazy loading type to Options, and load data when needed
-    if ((static_cast<int>(dims[0]) == mesh->GlobalNx)
-        and (static_cast<int>(dims[1]) == mesh->GlobalNy)) {
-      // Probably a Field2D
-
-      // Read just the local piece of the array
-      Matrix<BoutReal> value(mesh->LocalNx, mesh->LocalNy);
-
-      // Offset of this processor's data into the global array
-      adios2::Dims start = {static_cast<size_t>(mesh->MapGlobalX),
-                            static_cast<size_t>(mesh->MapGlobalY)};
-
-      // The size of the mapped region
-      adios2::Dims count = {static_cast<size_t>(mesh->MapCountX),
-                            static_cast<size_t>(mesh->MapCountY)};
-
-      // Where the actual data starts in data pointer (to exclude ghost cells)
-      adios2::Dims memStart = {static_cast<size_t>(mesh->MapLocalX),
-                               static_cast<size_t>(mesh->MapLocalY)};
-
-      // The actual size of data pointer in memory (including ghost cells)
-      adios2::Dims memCount = {static_cast<size_t>(mesh->LocalNx),
-                               static_cast<size_t>(mesh->LocalNy)};
-      variableD.SetSelection({start, count});
-      variableD.SetMemorySelection({memStart, memCount});
-      BoutReal* data = value.begin();
-      reader.Get<BoutReal>(variableD, data, adios2::Mode::Sync);
-      return Options(value);
+    if (type == "float" or type == "double") {
+      Matrix<BoutReal> value(selection.dims[0], selection.dims[1]);
+      return read_variable(io, reader, name, selection, value);
     }
-    if ((static_cast<int>(dims[0]) == mesh->GlobalNx)
-        and (static_cast<int>(dims[2]) == mesh->GlobalNz)) {
-      // Probably a FieldPerp
-
-      // Read just the local piece of the array
-      Matrix<BoutReal> value(mesh->LocalNx, mesh->LocalNz);
-
-      // Offset of this processor's data into the global array
-      adios2::Dims start = {static_cast<size_t>(mesh->MapGlobalX),
-                            static_cast<size_t>(mesh->MapGlobalZ)};
-
-      // The size of the mapped region
-      adios2::Dims count = {static_cast<size_t>(mesh->MapCountX),
-                            static_cast<size_t>(mesh->MapCountZ)};
-
-      // Where the actual data starts in data pointer (to exclude ghost cells)
-      adios2::Dims memStart = {static_cast<size_t>(mesh->MapLocalX),
-                               static_cast<size_t>(mesh->MapLocalZ)};
-
-      // The actual size of data pointer in memory (including ghost cells)
-      adios2::Dims memCount = {static_cast<size_t>(mesh->LocalNx),
-                               static_cast<size_t>(mesh->LocalNz)};
-      variableD.SetSelection({start, count});
-      variableD.SetMemorySelection({memStart, memCount});
-      BoutReal* data = value.begin();
-      reader.Get<BoutReal>(variableD, data, adios2::Mode::Sync);
-      return Options(value);
+    if (type == "int32_t") {
+      Matrix<int> value(selection.dims[0], selection.dims[1]);
+      return read_variable(io, reader, name, selection, value);
     }
-    Matrix<BoutReal> value(static_cast<int>(dims[0]), static_cast<int>(dims[1]));
-    BoutReal* data = value.begin();
-    reader.Get<BoutReal>(variableD, data, adios2::Mode::Sync);
-    return Options(value);
+    break;
   }
   case 3: {
-    if ((static_cast<int>(dims[0]) == mesh->GlobalNx)
-        and (static_cast<int>(dims[1]) == mesh->GlobalNy)
-        and (static_cast<int>(dims[2]) == mesh->GlobalNz)) {
-      // Global array. Read just this processor's part of it
-
-      Tensor<BoutReal> value(mesh->LocalNx, mesh->LocalNy, mesh->LocalNz);
-
-      // Offset of this processor's data into the global array
-      adios2::Dims start = {static_cast<size_t>(mesh->MapGlobalX),
-                            static_cast<size_t>(mesh->MapGlobalY),
-                            static_cast<size_t>(mesh->MapGlobalZ)};
-
-      // The size of the mapped region
-      adios2::Dims count = {static_cast<size_t>(mesh->MapCountX),
-                            static_cast<size_t>(mesh->MapCountY),
-                            static_cast<size_t>(mesh->MapCountZ)};
-
-      // Where the actual data starts in data pointer (to exclude ghost cells)
-      adios2::Dims memStart = {static_cast<size_t>(mesh->MapLocalX),
-                               static_cast<size_t>(mesh->MapLocalY),
-                               static_cast<size_t>(mesh->MapLocalZ)};
-
-      // The actual size of data pointer in memory (including ghost cells)
-      adios2::Dims memCount = {static_cast<size_t>(mesh->LocalNx),
-                               static_cast<size_t>(mesh->LocalNy),
-                               static_cast<size_t>(mesh->LocalNz)};
-
-      variableD.SetSelection({start, count});
-      variableD.SetMemorySelection({memStart, memCount});
-      BoutReal* data = value.begin();
-      reader.Get<BoutReal>(variableD, data, adios2::Mode::Sync);
-      return Options(value);
+    if (type == "float" or type == "double") {
+      Tensor<BoutReal> value(selection.dims[0], selection.dims[1], selection.dims[2]);
+      return read_variable(io, reader, name, selection, value);
     }
-    // Doesn't match global array size.
-    // Read the entire array, in case it can be handled later
-    Tensor<BoutReal> value(static_cast<int>(dims[0]), static_cast<int>(dims[1]),
-                           static_cast<int>(dims[2]));
-    BoutReal* data = value.begin();
-    reader.Get<BoutReal>(variableD, data, adios2::Mode::Sync);
-    return Options(value);
+    if (type == "int32_t") {
+      Tensor<int> value(selection.dims[0], selection.dims[1], selection.dims[2]);
+      return read_variable(io, reader, name, selection, value);
+    }
+    break;
   }
+  // Throw below
+  default:
+    break;
   }
-  throw BoutException("ADIOS reader failed to read '{}' of dimension {} in file '{}'",
-                      name, ndims, reader.Name());
+  auto dims_str = fmt::format("[{}]", fmt::join(dims, ", "));
+  throw BoutException(
+      "ADIOS reader failed to read '{}' (shape: {}, type: '{}') in file '{}'", name,
+      dims_str, type, reader.Name());
 }
 
 Options readVariable(adios2::Engine& reader, adios2::IO& io, const std::string& name,
                      const std::string& type) {
-#define declare_template_instantiation(T)           \
-  if (type == adios2::GetType<T>()) {               \
-    return readVariable<T>(reader, io, name, type); \
-  }
-  ADIOS2_FOREACH_ATTRIBUTE_PRIMITIVE_STDTYPE_1ARG(declare_template_instantiation)
-  declare_template_instantiation(std::string)
-#undef declare_template_instantiation
-      output_warn.write("ADIOS readVariable can't read type '{}' (variable '{}')", type,
-                        name);
+  if (type == adios2::GetType<std::string>()) {
+    return readVariable<std::string>(reader, io, name, type);
+  }
+  if (type == adios2::GetType<char>()) {
+    return readVariable<char>(reader, io, name, type);
+  }
+  if (type == adios2::GetType<int8_t>()) {
+    return readVariable<int8_t>(reader, io, name, type);
+  }
+  if (type == adios2::GetType<int16_t>()) {
+    return readVariable<int16_t>(reader, io, name, type);
+  }
+  if (type == adios2::GetType<int32_t>()) {
+    return readVariable<int32_t>(reader, io, name, type);
+  }
+  if (type == adios2::GetType<int64_t>()) {
+    return readVariable<int64_t>(reader, io, name, type);
+  }
+  if (type == adios2::GetType<uint8_t>()) {
+    return readVariable<uint8_t>(reader, io, name, type);
+  }
+  if (type == adios2::GetType<uint16_t>()) {
+    return readVariable<uint16_t>(reader, io, name, type);
+  }
+  if (type == adios2::GetType<uint32_t>()) {
+    return readVariable<uint32_t>(reader, io, name, type);
+  }
+  if (type == adios2::GetType<uint64_t>()) {
+    return readVariable<uint64_t>(reader, io, name, type);
+  }
+  if (type == adios2::GetType<float>()) {
+    return readVariable<float>(reader, io, name, type);
+  }
+  if (type == adios2::GetType<double>()) {
+    return readVariable<double>(reader, io, name, type);
+  }
+  if (type == adios2::GetType<long double>()) {
+    return readVariable<long double>(reader, io, name, type);
+  }
+
+  output_warn.write("ADIOS readVariable can't read type '{}' (variable '{}')", type,
+                    name);
   return Options{};
 }
 
 bool readAttribute(adios2::IO& io, const std::string& name, const std::string& type,
                    Options& result) {
   // Attribute is the part of 'name' after the last '/' separator
-  std::string attrname;
+  std::string attrname = name;
   auto pos = name.find_last_of('/');
-  if (pos == std::string::npos) {
-    attrname = name;
-  } else {
+  if (pos != std::string::npos) {
     attrname = name.substr(pos + 1);
   }
 
-#define declare_template_instantiation(T)                  \
-  if (type == adios2::GetType<T>()) {                      \
-    adios2::Attribute<T> a = io.InquireAttribute<T>(name); \
-    result.attributes[attrname] = *a.Data().data();        \
-    return true;                                           \
+  if (type == adios2::GetType<int>()) {
+    result.attributes[attrname] = *io.InquireAttribute<int>(name).Data().data();
+    return true;
+  }
+  if (type == adios2::GetType<BoutReal>()) {
+    result.attributes[attrname] = *io.InquireAttribute<BoutReal>(name).Data().data();
+    return true;
+  }
+  if (type == adios2::GetType<std::string>()) {
+    result.attributes[attrname] = *io.InquireAttribute<std::string>(name).Data().data();
+    return true;
   }
-  // Only some types of attributes are supported
-  //declare_template_instantiation(bool)
-  declare_template_instantiation(int) declare_template_instantiation(BoutReal)
-      declare_template_instantiation(std::string)
-#undef declare_template_instantiation
-          output_warn.write("ADIOS readAttribute can't read type '{}' (variable '{}')",
-                            type, name);
+
+  output_warn.write("ADIOS readAttribute can't read type '{}' (variable '{}')", type,
+                    name);
   return false;
 }
+} // namespace
+
+namespace bout {
+OptionsADIOS::OptionsADIOS(Options& options) : OptionsIO(options) {
+  if (options["file"].doc("File name. Defaults to <path>/<prefix>.pb").isSet()) {
+    filename = options["file"].as<std::string>();
+  } else {
+    // Both path and prefix must be set
+    filename = fmt::format("{}/{}.bp", options["path"].as<std::string>(),
+                           options["prefix"].as<std::string>());
+  }
+
+  file_mode = (options["append"].doc("Append to existing file?").withDefault<bool>(false))
+                  ? adios2::Mode::Append
+                  : adios2::Mode::Write;
+
+  singleWriteFile = options["singleWriteFile"].withDefault<bool>(false);
+}
 
 Options OptionsADIOS::read([[maybe_unused]] bool lazy) {
-  Timer timer("io");
+  const Timer timer("io");
 
   // Open file
-  ADIOSPtr adiosp = GetADIOSPtr();
+  ADIOSPtr const adiosp = GetADIOSPtr();
   adios2::IO io;
-  std::string ioname = "read_" + filename;
+  const std::string ioname = "read_" + filename;
   try {
     io = adiosp->AtIO(ioname);
   } catch (const std::invalid_argument& e) {
@@ -305,13 +403,19 @@ Options OptionsADIOS::read([[maybe_unused]] bool lazy) {
 }
 
 void OptionsADIOS::verifyTimesteps() const {
-  ADIOSStream& stream = ADIOSStream::ADIOSGetStream(filename);
-  stream.engine.EndStep();
-  stream.isInStep = false;
-  return;
+  // This doesn't _verify_ the timesteps, but does flush to disk.
+  // Maybe this is fine, because ADIOS2 doesn't require every variable
+  // to be in sync?
+  if (singleWriteFile) {
+    return;
+  }
+
+  ADIOSStream& stream = ADIOSStream::ADIOSGetStream(filename, file_mode);
+
+  stream.endStep();
 }
+} // namespace bout
 
-const std::vector<std::string> DIMS_NONE;
 const std::vector<std::string> DIMS_X = {"x"};
 const std::vector<std::string> DIMS_XY = {"x", "y"};
 const std::vector<std::string> DIMS_XZ = {"x", "z"};
@@ -321,10 +425,9 @@ namespace {
 using bout::utils::tuple_index_sequence;
 
 template <class Tuple, std::size_t... I>
-auto make_shape_impl(std::size_t first, Tuple&& t,
+auto make_shape_impl(std::size_t first, const Tuple& t,
                      std::index_sequence<I...> /* index */) {
-  return adios2::Dims{first,
-                      static_cast<std::size_t>(std::get<I>(std::forward<Tuple>(t)))...};
+  return adios2::Dims{first, static_cast<std::size_t>(std::get<I>(t))...};
 }
 // Return an `adios2::Dims` with value ``{first, value.shape()[0]...}``
 template <class T>
@@ -334,11 +437,11 @@ auto make_shape(std::size_t first, const T& value) {
 }
 
 template <class Tuple, std::size_t... I>
-auto make_start_impl(Tuple&& t, std::index_sequence<I...> /* index */) {
+auto make_start_impl(const Tuple& t, std::index_sequence<I...> /* index */) {
   // Hey look, a legitimate use of the comma operator to get a bunch
   // of zeros the length of the index_sequence!
   return adios2::Dims{static_cast<std::size_t>(BoutComm::rank()),
-                      (std::get<I>(std::forward<Tuple>(t)), std::size_t{0})...};
+                      (std::get<I>(t), std::size_t{0})...};
 }
 // Return an `adios2::Dims` with value ``{rank, 0...}``, with as many zeros as the dimension of ``T``
 template <class T>
@@ -346,206 +449,135 @@ auto make_start(const T& value) {
   const auto shape = value.shape();
   return make_start_impl(shape, tuple_index_sequence<decltype(shape)>{});
 }
-} // namespace
+
+template <std::size_t... I>
+auto make_dims_impl(std::index_sequence<I...> /*index*/) {
+  using namespace std::string_literals;
+  return std::vector{"rank"s, (fmt::format("dim_{}", I))...};
+}
+// Return vector of dimension names: `{"rank", "dim_0", ...}`
+template <class T>
+auto make_dims(const T& value) {
+  return make_dims_impl(tuple_index_sequence<decltype(value.shape())>{});
+}
 
 /// Visit a variant type, and put the data into a NcVar
 struct ADIOSPutVarVisitor {
-  ADIOSPutVarVisitor(const std::string& name, ADIOSStream& stream)
+  ADIOSPutVarVisitor(const std::string& name, bout::ADIOSStream& stream)
       : varname(name), stream(stream) {}
   template <typename T>
   void operator()(const T& value) {
     adios2::Variable<T> var = stream.GetValueVariable<T>(varname);
-    stream.engine.Put(var, value);
+    stream.engine().Put(var, value);
   }
 
-  void operator()(const Array<int>& value) { put_helper(value); }
-  void operator()(const Array<BoutReal>& value) { put_helper(value); }
-  void operator()(const Matrix<int>& value) { put_helper(value); }
-  void operator()(const Matrix<BoutReal>& value) { put_helper(value); }
-  void operator()(const Tensor<int>& value) { put_helper(value); }
-  void operator()(const Tensor<BoutReal>& value) { put_helper(value); }
-
-private:
-  const std::string& varname;
-  ADIOSStream& stream;
-
-  // helper for `Array`, `Matrix`, `Tensor`
-  template <template <class> class T, class U>
-  void put_helper(const T<U>& value) {
-    const auto shape = make_shape(static_cast<std::size_t>(BoutComm::size()), value);
-    auto var = stream.GetArrayVariable<U>(varname, shape, DIMS_NONE, BoutComm::rank());
-    var.SetSelection({make_start(value), make_shape(1, value)});
-    stream.engine.Put<U>(var, value.begin());
+  void operator()(const Array<int>& value) { amt_put_helper(value); }
+  void operator()(const Array<BoutReal>& value) { amt_put_helper(value); }
+  void operator()(const Matrix<int>& value) { amt_put_helper(value); }
+  void operator()(const Matrix<BoutReal>& value) { amt_put_helper(value); }
+  void operator()(const Tensor<int>& value) { amt_put_helper(value); }
+  void operator()(const Tensor<BoutReal>& value) { amt_put_helper(value); }
+  void operator()(bool value) {
+    // Scalars are only written from processor 0
+    if (BoutComm::rank() != 0) {
+      return;
+    }
+    stream.engine().Put(stream.GetValueVariable<int>(varname), static_cast<int>(value));
   }
-};
-
-template <>
-void ADIOSPutVarVisitor::operator()<bool>(const bool& value) {
-  // Scalars are only written from processor 0
-  if (BoutComm::rank() != 0) {
-    return;
+  void operator()(int value) {
+    // Scalars are only written from processor 0
+    if (BoutComm::rank() != 0) {
+      return;
+    }
+    stream.engine().Put(stream.GetValueVariable<int>(varname), value);
   }
-  adios2::Variable<int> var = stream.GetValueVariable<int>(varname);
-  stream.engine.Put<int>(var, static_cast<int>(value));
-}
-
-template <>
-void ADIOSPutVarVisitor::operator()<int>(const int& value) {
-  // Scalars are only written from processor 0
-  if (BoutComm::rank() != 0) {
-    return;
+  void operator()(BoutReal value) {
+    // Scalars are only written from processor 0
+    if (BoutComm::rank() != 0) {
+      return;
+    }
+    stream.engine().Put(stream.GetValueVariable<BoutReal>(varname), value);
   }
-  adios2::Variable<int> var = stream.GetValueVariable<int>(varname);
-  stream.engine.Put<int>(var, value);
-}
 
-template <>
-void ADIOSPutVarVisitor::operator()<BoutReal>(const BoutReal& value) {
-  // Scalars are only written from processor 0
-  if (BoutComm::rank() != 0) {
-    return;
+  void operator()(const std::string& value) {
+    // Scalars are only written from processor 0
+    if (BoutComm::rank() != 0) {
+      return;
+    }
+    stream.engine().Put<std::string>(stream.GetValueVariable<std::string>(varname), value,
+                                     adios2::Mode::Sync);
   }
-  adios2::Variable<BoutReal> var = stream.GetValueVariable<BoutReal>(varname);
-  stream.engine.Put<BoutReal>(var, value);
-}
 
-template <>
-void ADIOSPutVarVisitor::operator()<std::string>(const std::string& value) {
-  // Scalars are only written from processor 0
-  if (BoutComm::rank() != 0) {
-    return;
+  void operator()(const Field2D& value) {
+    // Empty dim_sizes is fine because we provide full list of dim_names
+    auto selection = Selection(DIMS_XY, {}, *value.getMesh());
+    auto var = stream.GetArrayVariable<BoutReal>(varname, selection.shape, DIMS_XY,
+                                                 BoutComm::rank());
+    var.SetSelection(selection.selection());
+    var.SetMemorySelection(selection.memorySelection());
+    stream.engine().Put(var, &value(0, 0));
   }
-  adios2::Variable<std::string> var = stream.GetValueVariable<std::string>(varname);
-  stream.engine.Put<std::string>(var, value, adios2::Mode::Sync);
-}
-
-template <>
-void ADIOSPutVarVisitor::operator()<Field2D>(const Field2D& value) {
-  // Get the mesh that describes how local data relates to global arrays
-  auto mesh = value.getMesh();
-
-  // The global size of this array includes boundary cells but not communication guard cells.
-  // In general this array will be sparse because it may have gaps.
-  adios2::Dims shape = {static_cast<size_t>(mesh->GlobalNx),
-                        static_cast<size_t>(mesh->GlobalNy)};
-
-  // Offset of this processor's data into the global array
-  adios2::Dims start = {static_cast<size_t>(mesh->MapGlobalX),
-                        static_cast<size_t>(mesh->MapGlobalY)};
-
-  // The size of the mapped region
-  adios2::Dims count = {static_cast<size_t>(mesh->MapCountX),
-                        static_cast<size_t>(mesh->MapCountY)};
-
-  // Where the actual data starts in data pointer (to exclude ghost cells)
-  adios2::Dims memStart = {static_cast<size_t>(mesh->MapLocalX),
-                           static_cast<size_t>(mesh->MapLocalY)};
-
-  // The actual size of data pointer in memory (including ghost cells)
-  adios2::Dims memCount = {static_cast<size_t>(value.getNx()),
-                           static_cast<size_t>(value.getNy())};
-
-  adios2::Variable<BoutReal> var =
-      stream.GetArrayVariable<BoutReal>(varname, shape, DIMS_XY, BoutComm::rank());
-  var.SetSelection({start, count});
-  var.SetMemorySelection({memStart, memCount});
-  stream.engine.Put<BoutReal>(var, &value(0, 0));
-}
-
-template <>
-void ADIOSPutVarVisitor::operator()<Field3D>(const Field3D& value) {
-  // Get the mesh that describes how local data relates to global arrays
-  auto mesh = value.getMesh();
-
-  // The global size of this array includes boundary cells but not communication guard cells.
-  // In general this array will be sparse because it may have gaps.
-  adios2::Dims shape = {static_cast<size_t>(mesh->GlobalNx),
-                        static_cast<size_t>(mesh->GlobalNy),
-                        static_cast<size_t>(mesh->GlobalNz)};
-
-  // Offset of this processor's data into the global array
-  adios2::Dims start = {static_cast<size_t>(mesh->MapGlobalX),
-                        static_cast<size_t>(mesh->MapGlobalY),
-                        static_cast<size_t>(mesh->MapGlobalZ)};
-
-  // The size of the mapped region
-  adios2::Dims count = {static_cast<size_t>(mesh->MapCountX),
-                        static_cast<size_t>(mesh->MapCountY),
-                        static_cast<size_t>(mesh->MapCountZ)};
-
-  // Where the actual data starts in data pointer (to exclude ghost cells)
-  adios2::Dims memStart = {static_cast<size_t>(mesh->MapLocalX),
-                           static_cast<size_t>(mesh->MapLocalY),
-                           static_cast<size_t>(mesh->MapLocalZ)};
-
-  // The actual size of data pointer in memory (including ghost cells)
-  adios2::Dims memCount = {static_cast<size_t>(value.getNx()),
-                           static_cast<size_t>(value.getNy()),
-                           static_cast<size_t>(value.getNz())};
-
-  adios2::Variable<BoutReal> var =
-      stream.GetArrayVariable<BoutReal>(varname, shape, DIMS_XYZ, BoutComm::rank());
-  var.SetSelection({start, count});
-  var.SetMemorySelection({memStart, memCount});
-  stream.engine.Put<BoutReal>(var, &value(0, 0, 0));
-}
-
-template <>
-void ADIOSPutVarVisitor::operator()<FieldPerp>(const FieldPerp& value) {
-  // Get the mesh that describes how local data relates to global arrays
-  auto mesh = value.getMesh();
 
-  // The global size of this array includes boundary cells but not communication guard cells.
-  // In general this array will be sparse because it may have gaps.
-  adios2::Dims shape = {static_cast<size_t>(mesh->GlobalNx),
-                        static_cast<size_t>(mesh->GlobalNz)};
-
-  // Offset of this processor's data into the global array
-  adios2::Dims start = {static_cast<size_t>(mesh->MapGlobalX),
-                        static_cast<size_t>(mesh->MapGlobalZ)};
+  void operator()(const Field3D& value) {
+    // Empty dim_sizes is fine because we provide full list of dim_names
+    auto selection = Selection(DIMS_XYZ, {}, *value.getMesh());
+    auto var = stream.GetArrayVariable<BoutReal>(varname, selection.shape, DIMS_XYZ,
+                                                 BoutComm::rank());
+    var.SetSelection(selection.selection());
+    var.SetMemorySelection(selection.memorySelection());
+    stream.engine().Put(var, &value(0, 0, 0));
+  }
 
-  // The size of the mapped region
-  adios2::Dims count = {static_cast<size_t>(mesh->MapCountX),
-                        static_cast<size_t>(mesh->MapCountZ)};
+  void operator()(const FieldPerp& value) {
+    // Empty dim_sizes is fine because we provide full list of dim_names
+    auto selection = Selection(DIMS_XZ, {}, *value.getMesh());
+    auto var = stream.GetArrayVariable<BoutReal>(varname, selection.shape, DIMS_XZ,
+                                                 BoutComm::rank());
+    var.SetSelection(selection.selection());
+    var.SetMemorySelection(selection.memorySelection());
+    stream.engine().Put<BoutReal>(var, &value(0, 0));
+  }
 
-  // Where the actual data starts in data pointer (to exclude ghost cells)
-  adios2::Dims memStart = {static_cast<size_t>(mesh->MapLocalX),
-                           static_cast<size_t>(mesh->MapLocalZ)};
+private:
+  // Ok to keep const/refs here as visitor is only a temporary
+  const std::string& varname; // NOLINT(*-avoid-const-or-ref-data-members)
+  bout::ADIOSStream& stream;  // NOLINT(*-avoid-const-or-ref-data-members)
 
-  // The actual size of data pointer in memory (including ghost cells)
-  adios2::Dims memCount = {static_cast<size_t>(value.getNx()),
-                           static_cast<size_t>(value.getNz())};
-
-  adios2::Variable<BoutReal> var =
-      stream.GetArrayVariable<BoutReal>(varname, shape, DIMS_XZ, BoutComm::rank());
-  var.SetSelection({start, count});
-  var.SetMemorySelection({memStart, memCount});
-  stream.engine.Put<BoutReal>(var, &value(0, 0));
-}
+  // helper for `Array`, `Matrix`, `Tensor`
+  template <template <class> class T, class U>
+  void amt_put_helper(const T<U>& value) {
+    const auto shape = make_shape(static_cast<std::size_t>(BoutComm::size()), value);
+    auto var =
+        stream.GetArrayVariable<U>(varname, shape, make_dims(value), BoutComm::rank());
+    var.SetSelection({make_start(value), make_shape(1, value)});
+    stream.engine().Put<U>(var, value.begin());
+  }
+};
 
 /// Visit a variant type, and put the data into a NcVar
 struct ADIOSPutAttVisitor {
   ADIOSPutAttVisitor(const std::string& varname, const std::string& attrname,
-                     ADIOSStream& stream)
+                     bout::ADIOSStream& stream)
       : varname(varname), attrname(attrname), stream(stream) {}
   template <typename T>
   void operator()(const T& value) {
     stream.io.DefineAttribute<T>(attrname, value, varname, "/", false);
   }
 
+  void operator()(bool value) {
+    stream.io.DefineAttribute<int>(attrname, static_cast<int>(value), varname, "/",
+                                   false);
+  }
+
 private:
-  const std::string& varname;
-  const std::string& attrname;
-  ADIOSStream& stream;
+  // Ok to keep const/refs here as visitor is only a temporary
+  const std::string& varname;  // NOLINT(*-avoid-const-or-ref-data-members)
+  const std::string& attrname; // NOLINT(*-avoid-const-or-ref-data-members)
+  bout::ADIOSStream& stream;   // NOLINT(*-avoid-const-or-ref-data-members)
 };
 
-template <>
-void ADIOSPutAttVisitor::operator()<bool>(const bool& value) {
-  stream.io.DefineAttribute<int>(attrname, (int)value, varname, "/", false);
-}
-
-void writeGroup(const Options& options, ADIOSStream& stream, const std::string& groupname,
-                const std::string& time_dimension) {
+void writeGroup(const Options& options, bout::ADIOSStream& stream,
+                const std::string& groupname, const std::string& time_dimension) {
 
   for (const auto& childpair : options.getChildren()) {
     const auto& name = childpair.first;
@@ -579,11 +611,12 @@ void writeGroup(const Options& options, ADIOSStream& stream, const std::string&
 
         // Write the variable
         // Note: ADIOS2 uses '/' to as a group separator; BOUT++ uses ':'
-        std::string varname = groupname.empty() ? name : groupname + "/" + name;
+        const std::string varname =
+            groupname.empty() ? name : fmt::format("{}/{}", groupname, name);
         bout::utils::visit(ADIOSPutVarVisitor(varname, stream), child.value);
 
         // Write attributes
-        if (!BoutComm::rank()) {
+        if (BoutComm::rank() == 0) {
           for (const auto& attribute : child.attributes) {
             const std::string& att_name = attribute.first;
             const auto& att = attribute.second;
@@ -598,54 +631,26 @@ void writeGroup(const Options& options, ADIOSStream& stream, const std::string&
     }
   }
 }
+} // namespace
 
+namespace bout {
 /// Write options to file
 void OptionsADIOS::write(const Options& options, const std::string& time_dim) {
-  Timer timer("io");
+  const Timer timer("io");
 
   // ADIOSStream is just a BOUT++ object, it does not create anything inside ADIOS
-  ADIOSStream& stream = ADIOSStream::ADIOSGetStream(filename);
-
-  // Need to have an adios2::IO object first, which can only be created once.
-  if (!stream.io) {
-    ADIOSPtr adiosp = GetADIOSPtr();
-    std::string ioname = "write_" + filename;
-    try {
-      stream.io = adiosp->AtIO(ioname);
-    } catch (const std::invalid_argument& e) {
-      stream.io = adiosp->DeclareIO(ioname);
-      stream.io.SetEngine("BP5");
-    }
-  }
+  ADIOSStream& stream = ADIOSStream::ADIOSGetStream(filename, file_mode);
 
-  /* Open file once and keep it open, close in stream desctructor
-     or close after writing if singleWriteFile == true
-    */
-  if (!stream.engine) {
-    adios2::Mode amode =
-        (file_mode == FileMode::append ? adios2::Mode::Append : adios2::Mode::Write);
-    stream.engine = stream.io.Open(filename, amode);
-    if (!stream.engine) {
-      throw BoutException("Could not open ADIOS file '{:s}' for writing", filename);
-    }
-  }
-
-  /* Multiple write() calls allowed in a single adios step to output multiple
-     Options objects in the same step. verifyTimesteps() will indicate the 
-     completion of the step (and adios will publish the step).
-  */
-  if (!stream.isInStep) {
-    stream.engine.BeginStep();
-    stream.isInStep = true;
-    stream.adiosStep = stream.engine.CurrentStep();
-  }
+  // Multiple write() calls allowed in a single adios step to output multiple
+  // Options objects in the same step. verifyTimesteps() will indicate the
+  // completion of the step (and adios will publish the step).
+  stream.beginStep();
 
   writeGroup(options, stream, "", time_dim);
 
-  /* In singleWriteFile mode, we complete the step and close the file */
+  // In singleWriteFile mode, we complete the step and close the file
   if (singleWriteFile) {
-    stream.engine.EndStep();
-    stream.engine.Close();
+    stream.finish();
   }
 }
 
diff --git a/src/sys/options/options_adios.hxx b/src/sys/options/options_adios.hxx
index 6e81508ca3..d3066378ce 100644
--- a/src/sys/options/options_adios.hxx
+++ b/src/sys/options/options_adios.hxx
@@ -5,19 +5,21 @@
 #define OPTIONS_ADIOS_H
 
 #include "bout/build_defines.hxx"
-#include "bout/options.hxx"
 #include "bout/options_io.hxx"
 
 #if !BOUT_HAS_ADIOS2
 
 namespace {
-bout::RegisterUnavailableOptionsIO
+const bout::RegisterUnavailableOptionsIO
     registerunavailableoptionsadios("adios", "BOUT++ was not configured with ADIOS2");
 }
 
 #else
 
-#include <memory>
+#include "bout/options.hxx"
+
+#include <adios2.h> // IWYU pragma: keep
+
 #include <string>
 
 namespace bout {
@@ -44,7 +46,7 @@ public:
 
   OptionsADIOS(const OptionsADIOS&) = delete;
   OptionsADIOS(OptionsADIOS&&) noexcept = default;
-  ~OptionsADIOS() = default;
+  ~OptionsADIOS() override = default;
 
   OptionsADIOS& operator=(const OptionsADIOS&) = delete;
   OptionsADIOS& operator=(OptionsADIOS&&) noexcept = default;
@@ -61,20 +63,15 @@ public:
   void verifyTimesteps() const override;
 
 private:
-  enum class FileMode {
-    replace, ///< Overwrite file when writing
-    append   ///< Append to file when writing
-  };
-
   /// Name of the file on disk
   std::string filename;
   /// How to open the file for writing
-  FileMode file_mode{FileMode::replace};
+  adios2::Mode file_mode{adios2::Mode::Write};
   bool singleWriteFile = false;
 };
 
 namespace {
-RegisterOptionsIO<OptionsADIOS> registeroptionsadios("adios");
+const RegisterOptionsIO<OptionsADIOS> registeroptionsadios("adios");
 }
 
 } // namespace bout
diff --git a/tests/integrated/CMakeLists.txt b/tests/integrated/CMakeLists.txt
index e11403efb9..3696ef0372 100644
--- a/tests/integrated/CMakeLists.txt
+++ b/tests/integrated/CMakeLists.txt
@@ -32,7 +32,6 @@ add_subdirectory(test-laplacexz)
 add_subdirectory(test-multigrid_laplace)
 add_subdirectory(test-naulin-laplace)
 add_subdirectory(test-options-netcdf)
-add_subdirectory(test-options-adios)
 add_subdirectory(test-petsc_laplace)
 add_subdirectory(test-petsc_laplace_MAST-grid)
 add_subdirectory(test-restart-io)
diff --git a/tests/integrated/test-options-adios/CMakeLists.txt b/tests/integrated/test-options-adios/CMakeLists.txt
deleted file mode 100644
index cc61fabe57..0000000000
--- a/tests/integrated/test-options-adios/CMakeLists.txt
+++ /dev/null
@@ -1,6 +0,0 @@
-bout_add_integrated_test(test-options-adios
-  SOURCES test-options-adios.cxx
-  USE_RUNTEST
-  USE_DATA_BOUT_INP
-  REQUIRES BOUT_HAS_ADIOS2
-  )
diff --git a/tests/integrated/test-options-adios/data/BOUT.inp b/tests/integrated/test-options-adios/data/BOUT.inp
deleted file mode 100644
index fa0f6d3681..0000000000
--- a/tests/integrated/test-options-adios/data/BOUT.inp
+++ /dev/null
@@ -1,6 +0,0 @@
-
-
-[mesh]
-nx = 5
-ny = 2
-nz = 2
diff --git a/tests/integrated/test-options-adios/runtest b/tests/integrated/test-options-adios/runtest
deleted file mode 100755
index 03a83fc0ba..0000000000
--- a/tests/integrated/test-options-adios/runtest
+++ /dev/null
@@ -1,74 +0,0 @@
-#!/usr/bin/env python3
-
-# Note: This test requires NCDF4, whereas on Travis NCDF is used
-# requires: netcdf
-# requires: adios
-# requires: not legacy_netcdf
-
-from boututils.datafile import DataFile
-from boututils.run_wrapper import build_and_log, shell, launch
-from boutdata.data import BoutOptionsFile
-
-import math
-import numpy as np
-
-build_and_log("options-netcdf test")
-shell("rm -f test-out.ini")
-shell("rm -f test-out.nc")
-
-# Create a NetCDF input file
-with DataFile("test.nc", create=True, format="NETCDF4") as f:
-    f.write("int", 42)
-    f.write("real", 3.1415)
-    f.write("string", "hello")
-
-# run BOUT++
-launch("./test-options-adios", nproc=1, mthread=1)
-
-# Check the output INI file
-result = BoutOptionsFile("test-out.ini")
-
-print(result)
-
-assert result["int"] == 42
-assert math.isclose(result["real"], 3.1415)
-assert result["string"] == "hello"
-
-print("Checking saved ADIOS2 test-out file -- Not implemented")
-
-# Check the output NetCDF file
-# with DataFile("test-out.nc") as f:
-#    assert f["int"] == 42
-#    assert math.isclose(f["real"], 3.1415)
-#    assert result["string"] == "hello"
-
-print("Checking saved settings.ini")
-
-# Check the settings.ini file, coming from BOUT.inp
-# which is converted to NetCDF, read in, then written again
-settings = BoutOptionsFile("settings.ini")
-
-assert settings["mesh"]["nx"] == 5
-assert settings["mesh"]["ny"] == 2
-
-print("Checking saved fields.bp -- Not implemented")
-
-# with DataFile("fields.nc") as f:
-#    assert f["f2d"].shape == (5, 6)  # Field2D
-#    assert f["f3d"].shape == (5, 6, 2)  # Field3D
-#    assert f["fperp"].shape == (5, 2)  # FieldPerp
-#    assert np.allclose(f["f2d"], 1.0)
-#    assert np.allclose(f["f3d"], 2.0)
-#    assert np.allclose(f["fperp"], 3.0)
-
-print("Checking saved fields2.bp -- Not implemented")
-
-# with DataFile("fields2.nc") as f:
-#    assert f["f2d"].shape == (5, 6)  # Field2D
-#    assert f["f3d"].shape == (5, 6, 2)  # Field3D
-#    assert f["fperp"].shape == (5, 2)  # FieldPerp
-#    assert np.allclose(f["f2d"], 1.0)
-#    assert np.allclose(f["f3d"], 2.0)
-#    assert np.allclose(f["fperp"], 3.0)
-
-print(" => Passed")
diff --git a/tests/integrated/test-options-adios/test-options-adios.cxx b/tests/integrated/test-options-adios/test-options-adios.cxx
deleted file mode 100644
index ba0bbe3898..0000000000
--- a/tests/integrated/test-options-adios/test-options-adios.cxx
+++ /dev/null
@@ -1,125 +0,0 @@
-
-#include "bout/bout.hxx"
-
-#include "bout/options_io.hxx"
-#include "bout/optionsreader.hxx"
-#include "bout/utils.hxx"
-
-using bout::OptionsIO;
-
-int main(int argc, char** argv) {
-  BoutInitialise(argc, argv);
-
-  // Read values from a NetCDF file
-  auto file = bout::OptionsIO::create("test.nc");
-
-  auto values = file->read();
-
-  values.printUnused();
-
-  // Write to an INI text file
-  OptionsReader* reader = OptionsReader::getInstance();
-  reader->write(&values, "test-out.ini");
-
-  // Write to ADIOS file, by setting file type "adios"
-  OptionsIO::create({{"file", "test-out.bp"},
-                     {"type", "adios"},
-                     {"append", false},
-                     {"singleWriteFile", true}})
-      ->write(values);
-
-  ///////////////////////////
-
-  // Write the BOUT.inp settings to ADIOS file
-  OptionsIO::create({{"file", "settings.bp"},
-                     {"type", "adios"},
-                     {"append", false},
-                     {"singleWriteFile", true}})
-      ->write(Options::root());
-
-  // Read back in
-  auto settings = OptionsIO::create({{"file", "settings.bp"}, {"type", "adios"}})->read();
-
-  // Write to INI file
-  reader->write(&settings, "settings.ini");
-
-  ///////////////////////////
-  // Write fields
-
-  {
-    constexpr int nx = 2;
-    constexpr int ny = 3;
-    Matrix<int> matrix_in(nx, ny);
-    int count = 0;
-
-    for (int i = 0; i < nx; ++i) {
-      for (int j = 0; j < ny; ++j) {
-        matrix_in(i, j) = ++count;
-      }
-    }
-
-    Options fields{{"f2d", Field2D(1.0)},
-                   {"f3d", Field3D(2.0)},
-                   {"fperp", FieldPerp(3.0)},
-                   {"matrix", matrix_in}};
-    auto f = OptionsIO::create({{"file", "fields.bp"}, {"type", "adios"}});
-    /*
-     write() for adios only buffers data but does not guarantee writing to disk
-     unless singleWriteFile is set to true 
-     */
-    f->write(fields);
-    // indicate completion of step, required to get data on disk
-    f->verifyTimesteps();
-  }
-
-  ///////////////////////////
-  // Read fields
-
-  Options fields_in =
-      OptionsIO::create({{"file", "fields.bp"}, {"type", "adios"}})->read();
-
-  auto f2d = fields_in["f2d"].as<Field2D>(bout::globals::mesh);
-  auto f3d = fields_in["f3d"].as<Field3D>(bout::globals::mesh);
-  auto fperp = fields_in["fperp"].as<FieldPerp>(bout::globals::mesh);
-
-  Options fields2;
-  fields2["f2d"] = f2d;
-  fields2["f3d"] = f3d;
-  fields2["fperp"] = fperp;
-
-  // Write out again
-  auto f2 = OptionsIO::create({{"file", "fields2.bp"},
-                               {"type", "adios"},
-                               {"append", false},
-                               {"singleWriteFile", true}});
-  f2->write(fields2);
-
-  ///////////////////////////
-  // Time dependent values
-
-  Options data;
-  data["scalar"] = 1.0;
-  data["scalar"].attributes["time_dimension"] = "t";
-
-  data["field"] = Field3D(2.0);
-  data["field"].attributes["time_dimension"] = "t";
-
-  OptionsIO::create({{"file", "time.bp"},
-                     {"type", "adios"},
-                     {"append", false},
-                     {"singleWriteFile", true}})
-      ->write(data);
-
-  // Update time-dependent values
-  data["scalar"] = 2.0;
-  data["field"] = Field3D(3.0);
-
-  // Append data to file
-  OptionsIO::create({{"file", "time.bp"},
-                     {"type", "adios"},
-                     {"append", true},
-                     {"singleWriteFile", true}})
-      ->write(data);
-
-  BoutFinalise();
-};
diff --git a/tests/unit/CMakeLists.txt b/tests/unit/CMakeLists.txt
index 47253c508f..40aa207dea 100644
--- a/tests/unit/CMakeLists.txt
+++ b/tests/unit/CMakeLists.txt
@@ -80,6 +80,7 @@ set(serial_tests_source
   ./sys/test_expressionparser.cxx
   ./sys/test_msg_stack.cxx
   ./sys/test_options.cxx
+  ./sys/test_options_adios2.cxx
   ./sys/test_options_fields.cxx
   ./sys/test_options_netcdf.cxx
   ./sys/test_optionsreader.cxx
diff --git a/tests/unit/fake_mesh.hxx b/tests/unit/fake_mesh.hxx
index e6f78f8767..b945b46524 100644
--- a/tests/unit/fake_mesh.hxx
+++ b/tests/unit/fake_mesh.hxx
@@ -50,6 +50,17 @@ public:
     OffsetY = 0;
     OffsetZ = 0;
 
+    // These bits only for ADIOS2, also boring due to single process
+    MapCountX = nx - 2;
+    MapCountY = ny - 2;
+    MapCountZ = nz;
+    MapGlobalX = nx;
+    MapGlobalY = ny;
+    MapGlobalZ = nz;
+    MapLocalX = nx - 2;
+    MapLocalY = ny - 2;
+    MapLocalZ = nz;
+
     // Small "inner" region
     xstart = 1;
     xend = nx - 2;
diff --git a/tests/unit/sys/test_options_adios2.cxx b/tests/unit/sys/test_options_adios2.cxx
new file mode 100644
index 0000000000..13bbe3db54
--- /dev/null
+++ b/tests/unit/sys/test_options_adios2.cxx
@@ -0,0 +1,345 @@
+// Test reading and writing to ADIOS2
+
+#include "bout/build_defines.hxx"
+
+#if BOUT_HAS_ADIOS2
+
+#include "gtest/gtest.h"
+
+#include "test_extras.hxx"
+#include "test_tmpfiles.hxx"
+#include "bout/adios_object.hxx"
+#include "bout/field3d.hxx"
+#include "bout/mesh.hxx"
+#include "bout/options_io.hxx"
+#include "bout/utils.hxx"
+
+using bout::OptionsIO;
+
+#include <string>
+
+#include "fake_mesh_fixture.hxx"
+
+// Reuse the "standard" fixture for FakeMesh
+class OptionsAdios2Test : public FakeMeshFixture {
+public:
+  static void SetUpTestSuite() { bout::ADIOSInit(BoutComm::get()); }
+  static void TearDownTestSuite() { bout::ADIOSFinalize(); }
+
+  ~OptionsAdios2Test() override {
+    bout::ADIOSStream::ADIOSGetStream(filename, adios2::Mode::Deferred).engine().Close();
+  }
+
+  // A temporary filename
+  bout::testing::TempFile filename;
+  WithQuietOutput quiet{output_info};
+  Options file_options{
+      {"file", filename.string()},
+      {"type", "adios"},
+      {"singleWriteFile", true},
+  };
+  Options no_write_options{
+      {"file", filename.string()},
+      {"type", "adios"},
+      {"singleWriteFile", false},
+  };
+};
+
+TEST_F(OptionsAdios2Test, ReadWriteInt) {
+  {
+    Options options;
+    options["test"] = 42;
+
+    // Write the file
+    OptionsIO::create(file_options)->write(options);
+  }
+
+  // Read again
+  Options data = OptionsIO::create(file_options)->read();
+
+  EXPECT_EQ(data["test"], 42);
+}
+
+TEST_F(OptionsAdios2Test, ReadWriteString) {
+  {
+    Options options;
+    options["test"] = std::string{"hello"};
+
+    // Write file
+    OptionsIO::create(file_options)->write(options);
+  }
+
+  // Read file
+  Options data = OptionsIO::create(file_options)->read();
+
+  EXPECT_EQ(data["test"], std::string("hello"));
+}
+
+TEST_F(OptionsAdios2Test, ReadWriteField2D) {
+  {
+    Options options;
+    options["test"] = Field2D(1.0);
+
+    // Write file
+    OptionsIO::create(file_options)->write(options);
+  }
+
+  // Read file
+  Options data = OptionsIO::create(file_options)->read();
+
+  Field2D value = data["test"].as<Field2D>(bout::globals::mesh);
+
+  EXPECT_DOUBLE_EQ(value(0, 1), 1.0);
+  EXPECT_DOUBLE_EQ(value(1, 0), 1.0);
+}
+
+TEST_F(OptionsAdios2Test, ReadWriteField3D) {
+  {
+    Options options;
+    options["test"] = Field3D(2.4);
+
+    // Write file
+    OptionsIO::create(file_options)->write(options);
+  }
+
+  // Read file
+  Options data = OptionsIO::create(file_options)->read();
+
+  Field3D value = data["test"].as<Field3D>(bout::globals::mesh);
+
+  EXPECT_DOUBLE_EQ(value(0, 1, 0), 2.4);
+  EXPECT_DOUBLE_EQ(value(1, 0, 1), 2.4);
+  EXPECT_DOUBLE_EQ(value(1, 1, 1), 2.4);
+}
+
+TEST_F(OptionsAdios2Test, ReadWriteMatrixInt) {
+  constexpr int nx = 2;
+  constexpr int ny = 3;
+  Matrix<int> matrix_in(nx, ny);
+  int count = 0;
+
+  for (int i = 0; i < nx; ++i) {
+    for (int j = 0; j < ny; ++j) {
+      matrix_in(i, j) = ++count;
+    }
+  }
+
+  {
+    Options options{{"int_matrix", matrix_in}};
+    OptionsIO::create(file_options)->write(options);
+  }
+
+  Options data = OptionsIO::create(file_options)->read();
+
+  const auto matrix_out = data["int_matrix"].as<Matrix<int>>();
+  ASSERT_EQ(matrix_out(nx - 1, ny - 1), matrix_in(nx - 1, ny - 1));
+}
+
+TEST_F(OptionsAdios2Test, Groups) {
+  {
+    Options options;
+    options["test"]["key"] = 42;
+
+    // Write file
+    OptionsIO::create(file_options)->write(options);
+  }
+
+  // Read file
+  Options data = OptionsIO::create(file_options)->read();
+  EXPECT_EQ(data["test"]["key"], 42);
+}
+
+TEST_F(OptionsAdios2Test, AttributeInt) {
+  {
+    Options options;
+    options["test"] = 3;
+    options["test"].attributes["thing"] = 4;
+
+    // Write file
+    OptionsIO::create(file_options)->write(options);
+  }
+
+  // Read file
+  Options data = OptionsIO::create(file_options)->read();
+  EXPECT_EQ(data["test"].attributes["thing"].as<int>(), 4);
+}
+
+TEST_F(OptionsAdios2Test, AttributeBoutReal) {
+  {
+    Options options;
+    options["test"] = 3;
+    options["test"].attributes["thing"] = 3.14;
+
+    // Write file
+    OptionsIO::create(file_options)->write(options);
+  }
+
+  // Read file
+  Options data = OptionsIO::create(file_options)->read();
+  EXPECT_DOUBLE_EQ(data["test"].attributes["thing"].as<BoutReal>(), 3.14);
+}
+
+TEST_F(OptionsAdios2Test, AttributeString) {
+  {
+    Options options;
+    options["test"] = 3;
+    options["test"].attributes["thing"] = "hello";
+
+    // Write file
+    OptionsIO::create(file_options)->write(options);
+  }
+
+  // Read file
+  Options data = OptionsIO::create(file_options)->read();
+  EXPECT_EQ(data["test"].attributes["thing"].as<std::string>(), "hello");
+}
+
+TEST_F(OptionsAdios2Test, Field2DWriteCellCentre) {
+  {
+    Options options;
+    options["f2d"] = Field2D(2.0);
+
+    // Write file
+    OptionsIO::create(file_options)->write(options);
+  }
+
+  // Read file
+  Options data = OptionsIO::create(file_options)->read();
+
+  EXPECT_EQ(data["f2d"].attributes["cell_location"].as<std::string>(),
+            toString(CELL_CENTRE));
+}
+
+TEST_F(OptionsAdios2Test, Field2DWriteCellYLow) {
+  {
+    Options options;
+    options["f2d"] = Field2D(2.0, mesh_staggered).setLocation(CELL_YLOW);
+
+    // Write file
+    OptionsIO::create(file_options)->write(options);
+  }
+
+  // Read file
+  Options data = OptionsIO::create(file_options)->read();
+
+  EXPECT_EQ(data["f2d"].attributes["cell_location"].as<std::string>(),
+            toString(CELL_YLOW));
+}
+
+TEST_F(OptionsAdios2Test, Field3DWriteCellCentre) {
+  {
+    Options options;
+    options["f3d"] = Field3D(2.0);
+
+    // Write file
+    OptionsIO::create(file_options)->write(options);
+  }
+
+  // Read file
+  Options data = OptionsIO::create(file_options)->read();
+
+  EXPECT_EQ(data["f3d"].attributes["cell_location"].as<std::string>(),
+            toString(CELL_CENTRE));
+}
+
+TEST_F(OptionsAdios2Test, Field3DWriteCellYLow) {
+  {
+    Options options;
+    options["f3d"] = Field3D(2.0, mesh_staggered).setLocation(CELL_YLOW);
+
+    // Write file
+    OptionsIO::create(file_options)->write(options);
+  }
+
+  // Read file
+  Options data = OptionsIO::create(file_options)->read();
+
+  EXPECT_EQ(data["f3d"].attributes["cell_location"].as<std::string>(),
+            toString(CELL_YLOW));
+}
+
+TEST_F(OptionsAdios2Test, FieldPerpWriteCellCentre) {
+  {
+    Options options;
+    FieldPerp fperp(3.0);
+    fperp.setIndex(2);
+    options["fperp"] = fperp;
+
+    // Ensure MPI is initialised, otherwise we end up creating threads while
+    // the file is open, and the lock is not removed on closing.
+    fperp.getMesh()->getXcomm();
+
+    // Write file
+    OptionsIO::create(file_options)->write(options);
+  }
+
+  // Read file
+  Options data = OptionsIO::create(file_options)->read();
+
+  EXPECT_EQ(data["fperp"].attributes["cell_location"].as<std::string>(),
+            toString(CELL_CENTRE));
+  EXPECT_EQ(data["fperp"].attributes["yindex_global"].as<int>(), 2);
+}
+
+TEST_F(OptionsAdios2Test, VerifyTimesteps) {
+  {
+    Options options;
+    options["thing1"] = 1.0;
+    options["thing1"].attributes["time_dimension"] = "t";
+
+    auto f = OptionsIO::create(no_write_options);
+    f->write(options);
+    EXPECT_NO_THROW(f->verifyTimesteps());
+  }
+
+  {
+    Options options;
+    options["thing1"] = 2.0;
+    options["thing1"].attributes["time_dimension"] = "t";
+
+    options["thing2"] = 3.0;
+    options["thing2"].attributes["time_dimension"] = "t";
+
+    OptionsIO::create({{"type", "adios"}, {"file", filename.string()}, {"append", true}})
+        ->write(options);
+  }
+
+  // Doesn't throw, but it's less of a problem for ADIOS2 than netCDF?
+  // EXPECT_THROW(OptionsIO::create(no_write_options)->verifyTimesteps(), BoutException);
+}
+
+TEST_F(OptionsAdios2Test, WriteTimeDimension) {
+  {
+    Options options;
+    options["thing1"].assignRepeat(1.0);       // default time dim
+    options["thing2"].assignRepeat(2.0, "t2"); // non-default
+
+    // Only write non-default time dim
+    OptionsIO::create(file_options)->write(options, "t2");
+  }
+
+  Options data = OptionsIO::create(file_options)->read();
+
+  EXPECT_FALSE(data.isSet("thing1"));
+  EXPECT_TRUE(data.isSet("thing2"));
+}
+
+TEST_F(OptionsAdios2Test, WriteMultipleTimeDimensions) {
+  {
+    Options options;
+    options["thing1_t1"].assignRepeat(1.0); // default time dim
+    options["thing2_t1"].assignRepeat(1.0); // default time dim
+
+    options["thing3_t2"].assignRepeat(2.0, "t2"); // non-default
+    options["thing4_t2"].assignRepeat(2.0, "t2"); // non-default
+
+    // Write the non-default time dim twice
+    OptionsIO::create(no_write_options)->write(options, "t2");
+    OptionsIO::create(no_write_options)->write(options, "t2");
+    OptionsIO::create(no_write_options)->write(options, "t");
+  }
+
+  EXPECT_NO_THROW(OptionsIO::create(no_write_options)->verifyTimesteps());
+}
+
+#endif // Bout_Has_ADIOS2
diff --git a/tests/unit/test_tmpfiles.hxx b/tests/unit/test_tmpfiles.hxx
index 35d7579df2..6c37e7c792 100644
--- a/tests/unit/test_tmpfiles.hxx
+++ b/tests/unit/test_tmpfiles.hxx
@@ -30,7 +30,7 @@ public:
   TempFile& operator=(const TempFile&) = delete;
   TempFile& operator=(TempFile&&) = delete;
 
-  ~TempFile() { std::filesystem::remove(filename); }
+  ~TempFile() { std::filesystem::remove_all(filename); }
 
   // Enable conversions to std::string / const char*
   operator std::string() const { return filename.string(); }