diff --git a/gusto/timestepping/timestepper.py b/gusto/timestepping/timestepper.py
index d7c80391b..3bc241a90 100644
--- a/gusto/timestepping/timestepper.py
+++ b/gusto/timestepping/timestepper.py
@@ -25,10 +25,12 @@ def __init__(self, equation, io):
         Args:
             equation (:class:`PrognosticEquation`): the prognostic equation.
             io (:class:`IO`): the model's object for controlling input/output.
+            init_io (:bool): whether or not to set up the IO
         """
 
         self.equation = equation
         self.io = io
+        self.init_io = True    # flag so that IO is only set up once
         self.dt = self.equation.domain.dt
         self.t = self.equation.domain.t
         self.reference_profiles_initialised = False
@@ -189,16 +191,23 @@ def run(self, t, tmax, pick_up=False):
             tmax (float): the end time of the run
             pick_up: (bool): specify whether to pick_up from a previous run
         """
+        if self.init_io:
+            # Set up diagnostics, which may set up some fields necessary to pick up
+            self.io.setup_diagnostics(self.fields)
+            self.io.setup_log_courant(self.fields)
+            if self.equation.domain.mesh.extruded:
+                self.io.setup_log_courant(self.fields, component='horizontal')
+                self.io.setup_log_courant(self.fields, component='vertical')
+            if self.transporting_velocity != "prognostic":
+                self.io.setup_log_courant(self.fields, name='transporting_velocity',
+                                          expression=self.transporting_velocity)
+
+            # Set up dump, which may also include an initial dump
+            with timed_stage("Dump output"):
+                logger.debug('Dumping output to disk')
+                self.io.setup_dump(self.fields, t, pick_up)
 
-        # Set up diagnostics, which may set up some fields necessary to pick up
-        self.io.setup_diagnostics(self.fields)
-        self.io.setup_log_courant(self.fields)
-        if self.equation.domain.mesh.extruded:
-            self.io.setup_log_courant(self.fields, component='horizontal')
-            self.io.setup_log_courant(self.fields, component='vertical')
-        if self.transporting_velocity != "prognostic":
-            self.io.setup_log_courant(self.fields, name='transporting_velocity',
-                                      expression=self.transporting_velocity)
+            self.init_io = False
 
         if pick_up:
             # Pick up fields, and return other info to be picked up
@@ -213,11 +222,6 @@ def run(self, t, tmax, pick_up=False):
         else:
             self.step = 1
 
-        # Set up dump, which may also include an initial dump
-        with timed_stage("Dump output"):
-            logger.debug('Dumping output to disk')
-            self.io.setup_dump(self.fields, t, pick_up)
-
         self.log_field_stats()
 
         self.t.assign(t)