diff --git a/include/bout/solver.hxx b/include/bout/solver.hxx
index 446acefecb..14df410928 100644
--- a/include/bout/solver.hxx
+++ b/include/bout/solver.hxx
@@ -451,19 +451,19 @@ protected:
   ///
   /// There are two important things to note about how \p iter is
   /// passed along to each monitor:
-  /// - The solvers all start their iteration numbering from zero, so the
-  ///   initial state is calculated at \p iter = -1
+  /// - The initial state is written at \p iter = 0, and solver output
+  ///   steps are numbered from 1 to NOUT
   /// - Secondly, \p iter is passed along to each monitor *relative to
   ///   that monitor's period*
   ///
   /// In practice, this means that each monitor is called like:
   ///
   ///     monitor->call(solver, simulation_time,
-  ///                   ((iter + 1) / monitor->period) - 1,
+  ///                   iter / monitor->period,
   ///                   NOUT / monitor->period);
   ///
-  /// e.g. for a monitor with period 10, passing \p iter = 9 will
-  /// result in it being called with a value of `(9 + 1)/10 - 1 == 0`
+  /// e.g. for a monitor with period 10, passing \p iter = 10 will
+  /// result in it being called with a value of `10/10 == 1`
   int call_monitors(BoutReal simtime, int iter, int NOUT);
 
   /// Should timesteps be monitored?
diff --git a/manual/sphinx/user_docs/time_integration.rst b/manual/sphinx/user_docs/time_integration.rst
index 37bf9a6abf..2e0d0f1b32 100644
--- a/manual/sphinx/user_docs/time_integration.rst
+++ b/manual/sphinx/user_docs/time_integration.rst
@@ -1278,7 +1278,9 @@ implement the outputMonitor method of PhysicsModel::
     int outputMonitor(BoutReal simtime, int iter, int nout)
 
 The first input is the current simulation time, the second is the output
-number, and the last is the total number of outputs requested.
+number, and the last is the total number of outputs requested. If an initial
+dump is written, it is output number ``0``. Solver output steps are numbered
+from ``1`` to ``nout``, so ``iter == nout`` indicates the final output.
 This method is called by a monitor object PhysicsModel::modelMonitor, which
 writes the restart files at the same time. You can change the frequency at which
 the monitor is called by calling, in PhysicsModel::init::
@@ -1303,7 +1305,9 @@ returns an int::
 
 The first input is the solver object, the second is the current
 simulation time, the third is the output number, and the last is the
-total number of outputs requested. To get the solver to call this
+total number of outputs requested. As for ``outputMonitor()``, output number
+``0`` is reserved for the initial dump when it is written, and solver output
+steps are numbered from ``1`` to ``NOUT``. To get the solver to call this
 function every output time, define a `MyOutputMonitor` object as a member of your
 PhysicsModel::
 
diff --git a/src/solver/impls/adams_bashforth/adams_bashforth.cxx b/src/solver/impls/adams_bashforth/adams_bashforth.cxx
index 00bd8f4ff2..0cd39a5b6d 100644
--- a/src/solver/impls/adams_bashforth/adams_bashforth.cxx
+++ b/src/solver/impls/adams_bashforth/adams_bashforth.cxx
@@ -355,7 +355,7 @@ int AdamsBashforthSolver::run() {
   [[maybe_unused]] int nwasted = 0;
   [[maybe_unused]] int nwasted_following_fail = 0;
 
-  for (int s = 0; s < getNumberOutputSteps(); s++) {
+  for (int s = 1; s <= getNumberOutputSteps(); s++) {
     BoutReal target = simtime + getOutputTimestep();
 
     bool running = true;
diff --git a/src/solver/impls/arkode/arkode.cxx b/src/solver/impls/arkode/arkode.cxx
index 804e67a803..11e3561d0c 100644
--- a/src/solver/impls/arkode/arkode.cxx
+++ b/src/solver/impls/arkode/arkode.cxx
@@ -510,7 +510,7 @@ int ArkodeSolver::run() {
     throw BoutException("ArkodeSolver not initialised\n");
   }
 
-  for (int i = 0; i < getNumberOutputSteps(); i++) {
+  for (int i = 1; i <= getNumberOutputSteps(); i++) {
 
     /// Run the solver for one output timestep
     simtime = run(simtime + getOutputTimestep());
diff --git a/src/solver/impls/cvode/cvode.cxx b/src/solver/impls/cvode/cvode.cxx
index 4f72e65933..2710534a62 100644
--- a/src/solver/impls/cvode/cvode.cxx
+++ b/src/solver/impls/cvode/cvode.cxx
@@ -489,7 +489,7 @@ int CvodeSolver::run() {
     throw BoutException("CvodeSolver not initialised\n");
   }
 
-  for (int i = 0; i < getNumberOutputSteps(); i++) {
+  for (int i = 1; i <= getNumberOutputSteps(); i++) {
 
     /// Run the solver for one output timestep
     simtime = run(simtime + getOutputTimestep());
diff --git a/src/solver/impls/euler/euler.cxx b/src/solver/impls/euler/euler.cxx
index 9754a68826..599d556c00 100644
--- a/src/solver/impls/euler/euler.cxx
+++ b/src/solver/impls/euler/euler.cxx
@@ -66,7 +66,7 @@ int EulerSolver::init() {
 
 int EulerSolver::run() {
 
-  for (int s = 0; s < getNumberOutputSteps(); s++) {
+  for (int s = 1; s <= getNumberOutputSteps(); s++) {
     BoutReal target = simtime + getOutputTimestep();
 
     bool running = true;
diff --git a/src/solver/impls/ida/ida.cxx b/src/solver/impls/ida/ida.cxx
index 84a0bd4abd..68ddf7f90d 100644
--- a/src/solver/impls/ida/ida.cxx
+++ b/src/solver/impls/ida/ida.cxx
@@ -236,7 +236,7 @@ int IdaSolver::run() {
     throw BoutException("IdaSolver not initialised\n");
   }
 
-  for (int i = 0; i < getNumberOutputSteps(); i++) {
+  for (int i = 1; i <= getNumberOutputSteps(); i++) {
 
     /// Run the solver for one output timestep
     simtime = run(simtime + getOutputTimestep());
diff --git a/src/solver/impls/imex-bdf2/imex-bdf2.cxx b/src/solver/impls/imex-bdf2/imex-bdf2.cxx
index 4ccf27981b..e30ae4c743 100644
--- a/src/solver/impls/imex-bdf2/imex-bdf2.cxx
+++ b/src/solver/impls/imex-bdf2/imex-bdf2.cxx
@@ -751,7 +751,7 @@ int IMEXBDF2::run() {
 
   int internalCounter = 0; // Cumulative number of successful internal iterations
 
-  for (int s = 0; s < getNumberOutputSteps(); s++) {
+  for (int s = 1; s <= getNumberOutputSteps(); s++) {
     BoutReal cumulativeTime = 0.;
     int counter = 0; // How many iterations in this output step
 
diff --git a/src/solver/impls/petsc/petsc.cxx b/src/solver/impls/petsc/petsc.cxx
index 1e1e05596b..0d5c16aebd 100644
--- a/src/solver/impls/petsc/petsc.cxx
+++ b/src/solver/impls/petsc/petsc.cxx
@@ -234,7 +234,7 @@ PetscErrorCode PetscMonitor(TS ts, PetscInt UNUSED(step), PetscReal t, Vec X, vo
     s->load_vars(const_cast<BoutReal*>(x));
     PetscCall(VecRestoreArrayRead(interpolatedX, &x));
 
-    if (s->call_monitors(output_time, i++, s->getNumberOutputSteps()) != 0) {
+    if (s->call_monitors(output_time, ++i, s->getNumberOutputSteps()) != 0) {
       PetscFunctionReturn(1);
     }
 
diff --git a/src/solver/impls/power/power.cxx b/src/solver/impls/power/power.cxx
index 5dfccfedda..744bdfe46e 100644
--- a/src/solver/impls/power/power.cxx
+++ b/src/solver/impls/power/power.cxx
@@ -45,7 +45,7 @@ int PowerSolver::run() {
   // Make sure that f0 has a norm of 1
   divide(f0, norm(f0));
 
-  for (int s = 0; s < getNumberOutputSteps(); s++) {
+  for (int s = 1; s <= getNumberOutputSteps(); s++) {
 
     load_vars(std::begin(f0));
     run_rhs(curtime);
diff --git a/src/solver/impls/pvode/pvode.cxx b/src/solver/impls/pvode/pvode.cxx
index 38901a2a57..0723faa99c 100644
--- a/src/solver/impls/pvode/pvode.cxx
+++ b/src/solver/impls/pvode/pvode.cxx
@@ -294,7 +294,7 @@ int PvodeSolver::run() {
     throw BoutException("PvodeSolver not initialised\n");
   }
 
-  for (int i = 0; i < getNumberOutputSteps(); i++) {
+  for (int i = 1; i <= getNumberOutputSteps(); i++) {
 
     /// Run the solver for one output timestep
     simtime = run(simtime + getOutputTimestep());
diff --git a/src/solver/impls/rk3-ssp/rk3-ssp.cxx b/src/solver/impls/rk3-ssp/rk3-ssp.cxx
index 7c00d44b34..26319eaa9a 100644
--- a/src/solver/impls/rk3-ssp/rk3-ssp.cxx
+++ b/src/solver/impls/rk3-ssp/rk3-ssp.cxx
@@ -61,7 +61,7 @@ int RK3SSP::init() {
 
 int RK3SSP::run() {
 
-  for (int s = 0; s < getNumberOutputSteps(); s++) {
+  for (int s = 1; s <= getNumberOutputSteps(); s++) {
     BoutReal target = simtime + getOutputTimestep();
 
     BoutReal dt;
diff --git a/src/solver/impls/rk4/rk4.cxx b/src/solver/impls/rk4/rk4.cxx
index 35b3cc2268..bb7c14fab2 100644
--- a/src/solver/impls/rk4/rk4.cxx
+++ b/src/solver/impls/rk4/rk4.cxx
@@ -73,7 +73,7 @@ int RK4Solver::init() {
 }
 
 int RK4Solver::run() {
-  for (int s = 0; s < getNumberOutputSteps(); s++) {
+  for (int s = 1; s <= getNumberOutputSteps(); s++) {
     BoutReal target = simtime + getOutputTimestep();
 
     BoutReal dt;
diff --git a/src/solver/impls/rkgeneric/rkgeneric.cxx b/src/solver/impls/rkgeneric/rkgeneric.cxx
index edf8ed50e3..20f4cc39d6 100644
--- a/src/solver/impls/rkgeneric/rkgeneric.cxx
+++ b/src/solver/impls/rkgeneric/rkgeneric.cxx
@@ -80,7 +80,7 @@ void RKGenericSolver::resetInternalFields() {
 }
 
 int RKGenericSolver::run() {
-  for (int s = 0; s < getNumberOutputSteps(); s++) {
+  for (int s = 1; s <= getNumberOutputSteps(); s++) {
     BoutReal target = simtime + getOutputTimestep();
 
     BoutReal dt;
diff --git a/src/solver/impls/slepc/slepc.cxx b/src/solver/impls/slepc/slepc.cxx
index 2c52cb3b38..9887f184eb 100644
--- a/src/solver/impls/slepc/slepc.cxx
+++ b/src/solver/impls/slepc/slepc.cxx
@@ -579,7 +579,7 @@ void SlepcSolver::monitor(PetscInt its, PetscInt nconv, PetscScalar eigr[],
   static bool first = true;
   if (eigenValOnly && first) {
     first = false;
-    resetIterationCounter();
+    resetIterationCounter(1);
   }
 
   // Temporary eigenvalues, converted from the SLEPc eigenvalues
@@ -701,7 +701,7 @@ void SlepcSolver::analyseResults() {
   output << "Converged eigenvalues :\n"
             "\tIndex\tSlepc eig (mag.)\t\t\tBOUT eig (mag.)\n";
 
-  resetIterationCounter();
+  resetIterationCounter(1);
 
   // Declare and create vectors to store eigenfunctions
   Vec vecReal, vecImag;
diff --git a/src/solver/impls/snes/snes.cxx b/src/solver/impls/snes/snes.cxx
index 5b28ff413c..2667ad55ff 100644
--- a/src/solver/impls/snes/snes.cxx
+++ b/src/solver/impls/snes/snes.cxx
@@ -886,7 +886,7 @@ int SNESSolver::run() {
 
   BoutReal target = simtime;
   recent_failure_rate = 0.0;
-  for (int s = 0; s < getNumberOutputSteps(); s++) {
+  for (int s = 1; s <= getNumberOutputSteps(); s++) {
     target += getOutputTimestep();
 
     bool looping = true;
diff --git a/src/solver/impls/split-rk/split-rk.cxx b/src/solver/impls/split-rk/split-rk.cxx
index ac4d62ebe2..e5fb0f6a8e 100644
--- a/src/solver/impls/split-rk/split-rk.cxx
+++ b/src/solver/impls/split-rk/split-rk.cxx
@@ -95,7 +95,7 @@ int SplitRK::init() {
 
 int SplitRK::run() {
 
-  for (int step = 0; step < getNumberOutputSteps(); step++) {
+  for (int step = 1; step <= getNumberOutputSteps(); step++) {
     // Take an output step
 
     BoutReal target = simtime + getOutputTimestep();
diff --git a/tests/unit/solver/test_fakesolver.hxx b/tests/unit/solver/test_fakesolver.hxx
index 3d96b717e3..c8600e4649 100644
--- a/tests/unit/solver/test_fakesolver.hxx
+++ b/tests/unit/solver/test_fakesolver.hxx
@@ -19,6 +19,13 @@ public:
     if ((*options)["throw_run"].withDefault(false)) {
       throw BoutException("Deliberate exception in FakeSolver::run");
     }
+    if ((*options)["call_final_monitor"].withDefault(false)) {
+      const int ret =
+          call_monitors(simtime, getNumberOutputSteps(), getNumberOutputSteps());
+      if (ret != 0) {
+        return ret;
+      }
+    }
     return (*options)["fail_run"].withDefault(0);
   }
   bool run_called{false};
diff --git a/tests/unit/solver/test_solver.cxx b/tests/unit/solver/test_solver.cxx
index e9a612103c..525d52767f 100644
--- a/tests/unit/solver/test_solver.cxx
+++ b/tests/unit/solver/test_solver.cxx
@@ -971,6 +971,26 @@ TEST_F(SolverTest, BasicSolve) {
   EXPECT_TRUE(solver.run_called);
 }
 
+TEST_F(SolverTest, SolveCleansUpMonitorsAtFinalIteration) {
+  Options options;
+  options["call_final_monitor"] = true;
+  FakeSolver solver{&options};
+
+  StrictMock<MockMonitor> monitor;
+  solver.addMonitor(&monitor);
+
+  NiceMock<MockPhysicsModel> model{};
+  solver.setModel(&model);
+
+  Options::cleanup();
+
+  EXPECT_CALL(monitor, call(_, _, 0, nout));
+  EXPECT_CALL(monitor, call(_, _, nout, nout));
+  EXPECT_CALL(monitor, cleanup());
+
+  solver.solve(nout, timestep);
+}
+
 TEST_F(SolverTest, GetRunID) {
   Options options;
   FakeSolver solver{&options};