diff --git a/examples/CosmoVolume/cosmoVolume.yml b/examples/CosmoVolume/cosmoVolume.yml
index 2da45221043187402bb42028f4d03723cbb26688..1bc3b3d3a6522710a2f41745a81cac7bae88c575 100644
--- a/examples/CosmoVolume/cosmoVolume.yml
+++ b/examples/CosmoVolume/cosmoVolume.yml
@@ -1,6 +1,6 @@
 
 # Define the system of units to use internally. 
-UnitSystem:
+InternalUnitSystem:
   UnitMass_in_cgs:     1   # Grams
   UnitLength_in_cgs:   1   # Centimeters
   UnitVelocity_in_cgs: 1   # Centimeters per second
@@ -21,6 +21,17 @@ TimeIntegration:
   dt_min:     1e-6  # The minimal time-step size of the simulation (in internal units).
   dt_max:     1e-2  # The maximal time-step size of the simulation (in internal units).
 
+# Parameters governing the snapshots
+Snapshots:
+  basename:            cosmo # Common part of the name of output files
+  time_first:          0.    # Time of the first output (in internal units)
+  delta_time:          0.05  # Time difference between consecutive outputs (in internal units)
+  UnitMass_in_cgs:     1   # Grams
+  UnitLength_in_cgs:   1   # Centimeters
+  UnitVelocity_in_cgs: 1   # Centimeters per second
+  UnitCurrent_in_cgs:  1   # Amperes
+  UnitTemp_in_cgs:     1   # Kelvin
+
 # Parameters for the hydrodynamics scheme
 SPH:
   resolution_eta:        1.2349   # Target smoothing length in units of the mean inter-particle separation (1.2349 == 48Ngbs with the cubic spline kernel).
diff --git a/examples/ExternalPointMass/externalPointMass.yml b/examples/ExternalPointMass/externalPointMass.yml
index cf96535b64c078d063b6f773823ca085876c2ecc..f834cb7015ff8f94a5595ed754f9d473a4e2d4f8 100644
--- a/examples/ExternalPointMass/externalPointMass.yml
+++ b/examples/ExternalPointMass/externalPointMass.yml
@@ -1,6 +1,6 @@
 
 # Define the system of units to use internally. 
-UnitSystem:
+InternalUnitSystem:
   UnitMass_in_cgs:     1.9885e33     # Grams
   UnitLength_in_cgs:   3.0856776e21  # Centimeters
   UnitVelocity_in_cgs: 1e5           # Centimeters per second
@@ -21,6 +21,17 @@ TimeIntegration:
   dt_min:     1e-6  # The minimal time-step size of the simulation (in internal units).
   dt_max:     1e-3  # The maximal time-step size of the simulation (in internal units).
 
+# Parameters governing the snapshots
+Snapshots:
+  basename:            pointMass # Common part of the name of output files
+  time_first:          0.        # Time of the first output (in internal units)
+  delta_time:          0.02      # Time difference between consecutive outputs (in internal units)
+  UnitMass_in_cgs:     1.9885e33     # Grams
+  UnitLength_in_cgs:   3.0856776e21  # Centimeters
+  UnitVelocity_in_cgs: 1e5           # Centimeters per second
+  UnitCurrent_in_cgs:  1   # Amperes
+  UnitTemp_in_cgs:     1   # Kelvin
+
 # Parameters for the hydrodynamics scheme
 SPH:
   resolution_eta:        1.2349   # Target smoothing length in units of the mean inter-particle separation (1.2349 == 48Ngbs with the cubic spline kernel).
diff --git a/examples/SedovBlast/sedov.yml b/examples/SedovBlast/sedov.yml
index 55974b03b823befde8365cddab187f5a18c5bbb7..feeee40a1433d5a2f57f9a80c00b18a550a4ebee 100644
--- a/examples/SedovBlast/sedov.yml
+++ b/examples/SedovBlast/sedov.yml
@@ -1,6 +1,6 @@
 
 # Define the system of units to use internally. 
-UnitSystem:
+InternalUnitSystem:
   UnitMass_in_cgs:     1   # Grams
   UnitLength_in_cgs:   1   # Centimeters
   UnitVelocity_in_cgs: 1   # Centimeters per second
@@ -21,6 +21,17 @@ TimeIntegration:
   dt_min:     1e-7  # The minimal time-step size of the simulation (in internal units).
   dt_max:     1e-2  # The maximal time-step size of the simulation (in internal units).
 
+# Parameters governing the snapshots
+Snapshots:
+  basename:            sedov # Common part of the name of output files
+  time_first:          0.    # Time of the first output (in internal units)
+  delta_time:          0.1   # Time difference between consecutive outputs (in internal units)
+  UnitMass_in_cgs:     1   # Grams
+  UnitLength_in_cgs:   1   # Centimeters
+  UnitVelocity_in_cgs: 1   # Centimeters per second
+  UnitCurrent_in_cgs:  1   # Amperes
+  UnitTemp_in_cgs:     1   # Kelvin
+
 # Parameters for the hydrodynamics scheme
 SPH:
   resolution_eta:        1.2349   # Target smoothing length in units of the mean inter-particle separation (1.2349 == 48Ngbs with the cubic spline kernel).
diff --git a/examples/SodShock/sodShock.yml b/examples/SodShock/sodShock.yml
index ab43d6682b2a16304d364784efee530ad5289cef..6a253a29610768c72dbbf1b9cb0ff35955da477e 100644
--- a/examples/SodShock/sodShock.yml
+++ b/examples/SodShock/sodShock.yml
@@ -1,6 +1,6 @@
 
 # Define the system of units to use internally. 
-UnitSystem:
+InternalUnitSystem:
   UnitMass_in_cgs:     1   # Grams
   UnitLength_in_cgs:   1   # Centimeters
   UnitVelocity_in_cgs: 1   # Centimeters per second
@@ -21,6 +21,17 @@ TimeIntegration:
   dt_min:     1e-7  # The minimal time-step size of the simulation (in internal units).
   dt_max:     1e-2  # The maximal time-step size of the simulation (in internal units).
 
+# Parameters governing the snapshots
+Snapshots:
+  basename:            sod # Common part of the name of output files
+  time_first:          0.  # Time of the first output (in internal units)
+  delta_time:          0.1 # Time difference between consecutive outputs (in internal units)
+  UnitMass_in_cgs:     1   # Grams
+  UnitLength_in_cgs:   1   # Centimeters
+  UnitVelocity_in_cgs: 1   # Centimeters per second
+  UnitCurrent_in_cgs:  1   # Amperes
+  UnitTemp_in_cgs:     1   # Kelvin
+
 # Parameters for the hydrodynamics scheme
 SPH:
   resolution_eta:        1.2349   # Target smoothing length in units of the mean inter-particle separation (1.2349 == 48Ngbs with the cubic spline kernel).
diff --git a/examples/UniformBox/uniformBox.yml b/examples/UniformBox/uniformBox.yml
index 0474b0f8202effa73210ee2b459806f2376a37f2..5ff07586f3fe035ec5ccf6474edc732e4483bf10 100644
--- a/examples/UniformBox/uniformBox.yml
+++ b/examples/UniformBox/uniformBox.yml
@@ -1,6 +1,6 @@
 
 # Define the system of units to use internally. 
-UnitSystem:
+InternalUnitSystem:
   UnitMass_in_cgs:     1   # Grams
   UnitLength_in_cgs:   1   # Centimeters
   UnitVelocity_in_cgs: 1   # Centimeters per second
@@ -21,6 +21,17 @@ TimeIntegration:
   dt_min:     1e-6  # The minimal time-step size of the simulation (in internal units).
   dt_max:     1e-2  # The maximal time-step size of the simulation (in internal units).
 
+# Parameters governing the snapshots
+Snapshots:
+  basename:            uniformBox # Common part of the name of output files
+  time_first:          0.         # Time of the first output (in internal units)
+  delta_time:          0.01       # Time difference between consecutive outputs (in internal units)
+  UnitMass_in_cgs:     1   # Grams
+  UnitLength_in_cgs:   1   # Centimeters
+  UnitVelocity_in_cgs: 1   # Centimeters per second
+  UnitCurrent_in_cgs:  1   # Amperes
+  UnitTemp_in_cgs:     1   # Kelvin
+
 # Parameters for the hydrodynamics scheme
 SPH:
   resolution_eta:        1.2349   # Target smoothing length in units of the mean inter-particle separation (1.2349 == 48Ngbs with the cubic spline kernel).
diff --git a/examples/main.c b/examples/main.c
index e8341bfe787b3051b712f5cf906a0931a6589369..f4dd3d93e30e58bbefa6d243d4e76d32b1194709 100644
--- a/examples/main.c
+++ b/examples/main.c
@@ -291,7 +291,7 @@ int main(int argc, char *argv[]) {
   /* Initialize unit system and constants */
   struct UnitSystem us;
   struct phys_const prog_const;
-  units_init(&us, params);
+  units_init(&us, params, "InternalUnitSystem");
   phys_const_init(&us, &prog_const);
   if (myrank == 0 && verbose > 0) {
     message("Unit system: U_M = %e g.", us.UnitMass_in_cgs);
@@ -347,7 +347,8 @@ int main(int argc, char *argv[]) {
   if (!with_hydro) {
     free(parts);
     parts = NULL;
-    for (size_t k = 0; k < Ngpart; ++k) if(gparts[k].id > 0) error("Linking problem");
+    for (size_t k = 0; k < Ngpart; ++k)
+      if (gparts[k].id > 0) error("Linking problem");
     Ngas = 0;
   }
 
@@ -421,33 +422,13 @@ int main(int argc, char *argv[]) {
     fflush(stdout);
   }
 
+  /* Write the state of the system before starting time integration. */
+  if (!dry_run) engine_dump_snapshot(&e);
+
   /* Now that everything is ready, no need for the parameters any more */
   free(params);
   params = NULL;
 
-  int with_outputs = 1;
-  if (with_outputs && !dry_run) {
-    /* Write the state of the system before starting time integration. */
-    if (myrank == 0) clocks_gettime(&tic);
-#if defined(WITH_MPI)
-#if defined(HAVE_PARALLEL_HDF5)
-    write_output_parallel(&e, &us, myrank, nr_nodes, MPI_COMM_WORLD,
-                          MPI_INFO_NULL);
-#else
-    write_output_serial(&e, &us, myrank, nr_nodes, MPI_COMM_WORLD,
-                        MPI_INFO_NULL);
-#endif
-#else
-    write_output_single(&e, &us);
-#endif
-    if (myrank == 0 && verbose) {
-      clocks_gettime(&toc);
-      message("writing particle properties took %.3f %s.",
-              clocks_diff(&tic, &toc), clocks_getunit());
-      fflush(stdout);
-    }
-  }
-
 /* Init the runner history. */
 #ifdef HIST
   for (k = 0; k < runner_hist_N; k++) runner_hist_bins[k] = 0;
@@ -456,10 +437,10 @@ int main(int argc, char *argv[]) {
   /* Get some info to the user. */
   if (myrank == 0) {
     message(
-        "Running on %lld gas particles and %lld DM particles until t=%.3e with "
-        "%i threads and %i queues (dt_min=%.3e, dt_max=%.3e)...",
-        N_total[0], N_total[1], e.timeEnd, e.nr_threads, e.sched.nr_queues,
-        e.dt_min, e.dt_max);
+        "Running on %lld gas particles and %lld DM particles from t=%.3e until "
+        "t=%.3e with %d threads and %d queues (dt_min=%.3e, dt_max=%.3e)...",
+        N_total[0], N_total[1], e.timeBegin, e.timeEnd, e.nr_threads,
+        e.sched.nr_queues, e.dt_min, e.dt_max);
     fflush(stdout);
   }
 
@@ -504,28 +485,6 @@ int main(int argc, char *argv[]) {
     /* Take a step. */
     engine_step(&e);
 
-    if (with_outputs && j % 100 == 0) {
-
-      if (myrank == 0) clocks_gettime(&tic);
-#if defined(WITH_MPI)
-#if defined(HAVE_PARALLEL_HDF5)
-      write_output_parallel(&e, &us, myrank, nr_nodes, MPI_COMM_WORLD,
-                            MPI_INFO_NULL);
-#else
-      write_output_serial(&e, &us, myrank, nr_nodes, MPI_COMM_WORLD,
-                          MPI_INFO_NULL);
-#endif
-#else
-      write_output_single(&e, &us);
-#endif
-      if (myrank == 0 && verbose) {
-        clocks_gettime(&toc);
-        message("writing particle properties took %.3f %s.",
-                clocks_diff(&tic, &toc), clocks_getunit());
-        fflush(stdout);
-      }
-    }
-
     /* Dump the task data using the given frequency. */
     if (dump_tasks && (dump_tasks == 1 || j % dump_tasks == 1)) {
 #ifdef WITH_MPI
@@ -607,28 +566,8 @@ int main(int argc, char *argv[]) {
            (double)runner_hist_bins[k]);
 #endif
 
-  if (with_outputs) {
-
-    if (myrank == 0) clocks_gettime(&tic);
-/* Write final output. */
-#if defined(WITH_MPI)
-#if defined(HAVE_PARALLEL_HDF5)
-    write_output_parallel(&e, &us, myrank, nr_nodes, MPI_COMM_WORLD,
-                          MPI_INFO_NULL);
-#else
-    write_output_serial(&e, &us, myrank, nr_nodes, MPI_COMM_WORLD,
-                        MPI_INFO_NULL);
-#endif
-#else
-    write_output_single(&e, &us);
-#endif
-    if (myrank == 0 && verbose) {
-      clocks_gettime(&toc);
-      message("writing particle properties took %.3f %s.",
-              clocks_diff(&tic, &toc), clocks_getunit());
-      fflush(stdout);
-    }
-  }
+  /* Write final output. */
+  engine_dump_snapshot(&e);
 
 #ifdef WITH_MPI
   if ((res = MPI_Finalize()) != MPI_SUCCESS)
diff --git a/examples/parameter_example.yml b/examples/parameter_example.yml
index 228e748f3b2899639303597c92ea29ad09987313..e629ba78f48c0975ebd4405c071fc70ad3623dce 100644
--- a/examples/parameter_example.yml
+++ b/examples/parameter_example.yml
@@ -1,6 +1,6 @@
 
 # Define the system of units to use internally. 
-UnitSystem:
+InternalUnitSystem:
   UnitMass_in_cgs:     1   # Grams
   UnitLength_in_cgs:   1   # Centimeters
   UnitVelocity_in_cgs: 1   # Centimeters per second
@@ -21,6 +21,17 @@ TimeIntegration:
   dt_min:     1e-6  # The minimal time-step size of the simulation (in internal units).
   dt_max:     1e-2  # The maximal time-step size of the simulation (in internal units).
 
+# Parameters governing the snapshots
+Snapshots:
+  basename:   output      # Common part of the name of output files
+  time_first: 0.          # Time of the first output (in internal units)
+  delta_time: 0.01        # Time difference between consecutive outputs (in internal units)
+  UnitMass_in_cgs:     1  # Unit system for the outputs (Grams)
+  UnitLength_in_cgs:   1  # Unit system for the outputs (Centimeters)
+  UnitVelocity_in_cgs: 1  # Unit system for the outputs (Centimeters per second)
+  UnitCurrent_in_cgs:  1  # Unit system for the outputs (Amperes)
+  UnitTemp_in_cgs:     1  # Unit system for the outputs (Kelvin)
+
 # Parameters for the hydrodynamics scheme
 SPH:
   resolution_eta:        1.2349   # Target smoothing length in units of the mean inter-particle separation (1.2349 == 48Ngbs with the cubic spline kernel).
diff --git a/src/common_io.c b/src/common_io.c
index 2a635723d5bd4db7bce0a0172e8c083bf479ac32..143c1ca94160f2a43c40aa44384803ce721e4dbf 100644
--- a/src/common_io.c
+++ b/src/common_io.c
@@ -335,11 +335,15 @@ void writeCodeDescription(hid_t h_file) {
  *
  * @todo Use a proper XML library to avoid stupid copies.
  */
-FILE* prepareXMFfile() {
+FILE* prepareXMFfile(const char* baseName) {
   char buffer[1024];
 
-  FILE* xmfFile = fopen("output.xmf", "r");
-  FILE* tempFile = fopen("output_temp.xmf", "w");
+  char fileName[FILENAME_BUFFER_SIZE];
+  char tempFileName[FILENAME_BUFFER_SIZE];
+  snprintf(fileName, FILENAME_BUFFER_SIZE, "%s.xmf", baseName);
+  snprintf(tempFileName, FILENAME_BUFFER_SIZE, "%s_temp.xmf", baseName);
+  FILE* xmfFile = fopen(fileName, "r");
+  FILE* tempFile = fopen(tempFileName, "w");
 
   if (xmfFile == NULL) error("Unable to open current XMF file.");
 
@@ -355,8 +359,8 @@ FILE* prepareXMFfile() {
   fclose(xmfFile);
 
   /* We then copy the XMF file back up to the closing lines */
-  xmfFile = fopen("output.xmf", "w");
-  tempFile = fopen("output_temp.xmf", "r");
+  xmfFile = fopen(fileName, "w");
+  tempFile = fopen(tempFileName, "r");
 
   if (xmfFile == NULL) error("Unable to open current XMF file.");
 
@@ -369,7 +373,7 @@ FILE* prepareXMFfile() {
   }
   fprintf(xmfFile, "\n");
   fclose(tempFile);
-  remove("output_temp.xmf");
+  remove(tempFileName);
 
   return xmfFile;
 }
@@ -380,8 +384,11 @@ FILE* prepareXMFfile() {
  * @todo Exploit the XML nature of the XMF format to write a proper XML writer
  *and simplify all the XMF-related stuff.
  */
-void createXMFfile() {
-  FILE* xmfFile = fopen("output.xmf", "w");
+void createXMFfile(const char* baseName) {
+
+  char fileName[FILENAME_BUFFER_SIZE];
+  snprintf(fileName, FILENAME_BUFFER_SIZE, "%s.xmf", baseName);
+  FILE* xmfFile = fopen(fileName, "w");
 
   fprintf(xmfFile, "<?xml version=\"1.0\" ?> \n");
   fprintf(xmfFile, "<!DOCTYPE Xdmf SYSTEM \"Xdmf.dtd\" []> \n");
diff --git a/src/common_io.h b/src/common_io.h
index b7f3a1a317d69937dde8692eead8f00c75649477..70ed25993cb110a1faf42a41c08d316c81b54271 100644
--- a/src/common_io.h
+++ b/src/common_io.h
@@ -97,8 +97,8 @@ void writeAttribute_i(hid_t grp, char* name, int data);
 void writeAttribute_l(hid_t grp, char* name, long data);
 void writeAttribute_s(hid_t grp, char* name, const char* str);
 
-void createXMFfile();
-FILE* prepareXMFfile();
+void createXMFfile(const char* baseName);
+FILE* prepareXMFfile(const char* baseName);
 void writeXMFoutputheader(FILE* xmfFile, char* hdfFileName, float time);
 void writeXMFoutputfooter(FILE* xmfFile, int outputCount, float time);
 void writeXMFgroupheader(FILE* xmfFile, char* hdfFileName, size_t N,
diff --git a/src/engine.c b/src/engine.c
index 96c0285bfc4a7b833f2382889c2206a408fb339f..eabe688f293577c0fac1c1fb61d3a5f0780e25d4 100644
--- a/src/engine.c
+++ b/src/engine.c
@@ -58,6 +58,9 @@
 #include "minmax.h"
 #include "part.h"
 #include "partition.h"
+#include "parallel_io.h"
+#include "serial_io.h"
+#include "single_io.h"
 #include "timers.h"
 
 const char *engine_policy_names[13] = {
@@ -2046,6 +2049,25 @@ void engine_step(struct engine *e) {
   }
 #endif
 
+  /* Check for output */
+  while (ti_end_min >= e->ti_nextSnapshot && e->ti_nextSnapshot > 0) {
+
+    e->ti_old = e->ti_current;
+    e->ti_current = e->ti_nextSnapshot;
+    e->time = e->ti_current * e->timeBase + e->timeBegin;
+    e->timeOld = e->ti_old * e->timeBase + e->timeBegin;
+    e->timeStep = (e->ti_current - e->ti_old) * e->timeBase;
+
+    /* Drift everybody to the snapshot position */
+    engine_launch(e, e->nr_threads, 1 << task_type_drift, 0);
+
+    /* Dump... */
+    engine_dump_snapshot(e);
+
+    /* ... and find the next output time */
+    engine_compute_next_snapshot_time(e);
+  }
+
   /* Move forward in time */
   e->ti_old = e->ti_current;
   e->ti_current = ti_end_min;
@@ -2324,6 +2346,38 @@ void engine_split(struct engine *e, struct partition *initial_partition) {
 #endif
 }
 
+/**
+ * @brief Writes a snapshot with the current state of the engine
+ *
+ * @param e The #engine.
+ */
+void engine_dump_snapshot(struct engine *e) {
+
+  struct clocks_time time1, time2;
+  clocks_gettime(&time1);
+
+  if (e->verbose)
+    message("writing snapshot at t=%f", e->time);
+  
+/* Dump... */
+#if defined(WITH_MPI)
+#if defined(HAVE_PARALLEL_HDF5)
+  write_output_parallel(e, e->snapshotBaseName, e->snapshotUnits, e->nodeID,
+                        e->nr_nodes, MPI_COMM_WORLD, MPI_INFO_NULL);
+#else
+  write_output_serial(e, e->snapshotBaseName, e->snapshotUnits, e->nodeID,
+                      e->nr_nodes, MPI_COMM_WORLD, MPI_INFO_NULL);
+#endif
+#else
+  write_output_single(e, e->snapshotBaseName, e->snapshotUnits);
+#endif
+
+  clocks_gettime(&time2);
+  if (e->verbose)
+    message("writing particle properties took %.3f %s.",
+            (float)clocks_diff(&time1, &time2), clocks_getunit());
+}
+
 #if defined(HAVE_LIBNUMA) && defined(_GNU_SOURCE)
 static bool hyperthreads_present(void) {
 #ifdef __linux__
@@ -2387,6 +2441,15 @@ void engine_init(struct engine *e, struct space *s,
   e->ti_old = 0;
   e->ti_current = 0;
   e->timeStep = 0.;
+  e->timeBase = 0.;
+  e->timeFirstSnapshot =
+      parser_get_param_double(params, "Snapshots:time_first");
+  e->deltaTimeSnapshot =
+      parser_get_param_double(params, "Snapshots:delta_time");
+  e->ti_nextSnapshot = 0;
+  parser_get_param_string(params, "Snapshots:basename", e->snapshotBaseName);
+  e->snapshotUnits = malloc(sizeof(struct UnitSystem));
+  units_init(e->snapshotUnits, params, "Snapshots");
   e->dt_min = parser_get_param_double(params, "TimeIntegration:dt_min");
   e->dt_max = parser_get_param_double(params, "TimeIntegration:dt_max");
   e->file_stats = NULL;
@@ -2563,6 +2626,19 @@ void engine_init(struct engine *e, struct space *s,
     error("Maximal time-step size larger than the simulation run time t=%e",
           e->timeEnd - e->timeBegin);
 
+  /* Deal with outputs */
+  if (e->deltaTimeSnapshot < 0.)
+    error("Time between snapshots (%e) must be positive.",
+          e->deltaTimeSnapshot);
+
+  if (e->timeFirstSnapshot < e->timeBegin)
+    error(
+        "Time of first snapshot (%e) must be after the simulation start t=%e.",
+        e->timeFirstSnapshot, e->timeBegin);
+
+  /* Find the time of the first output */
+  engine_compute_next_snapshot_time(e);
+
 /* Construct types for MPI communications */
 #ifdef WITH_MPI
   part_create_mpi_types();
@@ -2667,3 +2743,25 @@ void engine_print_policy(struct engine *e) {
   fflush(stdout);
 #endif
 }
+
+/**
+ * @brief Computes the next time (on the time line) for a dump
+ *
+ * @param e The #engine.
+ */
+void engine_compute_next_snapshot_time(struct engine *e) {
+
+  for (double time = e->timeFirstSnapshot; time < e->timeEnd;
+       time += e->deltaTimeSnapshot) {
+
+    /* Output time on the integer timeline */
+    e->ti_nextSnapshot = (time - e->timeBegin) / e->timeBase;
+
+    if (e->ti_nextSnapshot > e->ti_current) break;
+  }
+
+  /* Be nice, talk... */
+  const float next_snapshot_time =
+      e->ti_nextSnapshot * e->timeBase + e->timeBegin;
+  if (e->verbose) message("Next output time set to t=%f", next_snapshot_time);
+}
diff --git a/src/engine.h b/src/engine.h
index 82e76481f8d1b17c6f23ee70b15b2013bb0fce64..f4bc3be00cacf128cbaf58ac5735e17f4fe000e4 100644
--- a/src/engine.h
+++ b/src/engine.h
@@ -47,6 +47,7 @@
 #include "partition.h"
 #include "physical_constants.h"
 #include "potentials.h"
+#include "units.h"
 
 /* Some constants. */
 enum engine_policy {
@@ -131,6 +132,13 @@ struct engine {
   /* Time base */
   double timeBase;
 
+  /* Snapshot information */
+  double timeFirstSnapshot;
+  double deltaTimeSnapshot;
+  int ti_nextSnapshot;
+  char snapshotBaseName[200];
+  struct UnitSystem *snapshotUnits;
+
   /* File for statistics */
   FILE *file_stats;
 
@@ -181,6 +189,8 @@ struct engine {
 
 /* Function prototypes. */
 void engine_barrier(struct engine *e, int tid);
+void engine_compute_next_snapshot_time(struct engine *e);
+void engine_dump_snapshot(struct engine *e);
 void engine_init(struct engine *e, struct space *s,
                  const struct swift_params *params, int nr_nodes, int nodeID,
                  int nr_threads, int policy, int verbose,
diff --git a/src/parallel_io.c b/src/parallel_io.c
index d1c739b59021f38b2259f82dd06c547e0e7c147d..4579be8f04ae687140413383e061988c9783b570 100644
--- a/src/parallel_io.c
+++ b/src/parallel_io.c
@@ -509,8 +509,12 @@ void read_ic_parallel(char* fileName, double dim[3], struct part** parts,
  *its XMF descriptor
  *
  * @param e The engine containing all the system.
- * @param us The UnitSystem used for the conversion of units
- *in the output
+ * @param baseName The common part of the snapshot file name.
+ * @param us The UnitSystem used for the conversion of units in the output.
+ * @param mpi_rank The MPI rank of this node.
+ * @param mpi_size The number of MPI ranks.
+ * @param comm The MPI communicator.
+ * @param info The MPI information object
  *
  * Creates an HDF5 output file and writes the particles
  *contained
@@ -522,9 +526,9 @@ void read_ic_parallel(char* fileName, double dim[3], struct part** parts,
  * Calls #error() if an error occurs.
  *
  */
-void write_output_parallel(struct engine* e, struct UnitSystem* us,
-                           int mpi_rank, int mpi_size, MPI_Comm comm,
-                           MPI_Info info) {
+void write_output_parallel(struct engine* e, const char* baseName,
+                           struct UnitSystem* us, int mpi_rank, int mpi_size,
+                           MPI_Comm comm, MPI_Info info) {
   hid_t h_file = 0, h_grp = 0, h_grpsph = 0;
   const size_t Ngas = e->s->nr_parts;
   const size_t Ntot = e->s->nr_gparts;
@@ -536,22 +540,19 @@ void write_output_parallel(struct engine* e, struct UnitSystem* us,
   static int outputCount = 0;
   FILE* xmfFile = 0;
 
-  /* Number of particles of each type */
-  // const size_t Ndm = Ntot - Ngas;
-
-  /* MATTHIEU: Temporary fix to preserve master */
+  /* Number of unassociated gparts */
   const size_t Ndm = Ntot > 0 ? Ntot - Ngas : 0;
-  /* MATTHIEU: End temporary fix */
 
   /* File name */
   char fileName[FILENAME_BUFFER_SIZE];
-  snprintf(fileName, FILENAME_BUFFER_SIZE, "output_%03i.hdf5", outputCount);
+  snprintf(fileName, FILENAME_BUFFER_SIZE, "%s_%03i.hdf5", baseName,
+           outputCount);
 
   /* First time, we need to create the XMF file */
-  if (outputCount == 0 && mpi_rank == 0) createXMFfile();
+  if (outputCount == 0 && mpi_rank == 0) createXMFfile(baseName);
 
   /* Prepare the XMF file for the new entry */
-  if (mpi_rank == 0) xmfFile = prepareXMFfile();
+  if (mpi_rank == 0) xmfFile = prepareXMFfile(baseName);
 
   /* Open HDF5 file */
   hid_t plist_id = H5Pcreate(H5P_FILE_ACCESS);
diff --git a/src/parallel_io.h b/src/parallel_io.h
index f3691cb29b8d5e7f17382f1f81ba230c3898a929..970ad8c41dcbc2df3a85b178f836e16926147788 100644
--- a/src/parallel_io.h
+++ b/src/parallel_io.h
@@ -36,9 +36,9 @@ void read_ic_parallel(char* fileName, double dim[3], struct part** parts,
                       int* periodic, int mpi_rank, int mpi_size, MPI_Comm comm,
                       MPI_Info info, int dry_run);
 
-void write_output_parallel(struct engine* e, struct UnitSystem* us,
-                           int mpi_rank, int mpi_size, MPI_Comm comm,
-                           MPI_Info info);
+void write_output_parallel(struct engine* e, const char* baseName,
+                           struct UnitSystem* us, int mpi_rank, int mpi_size,
+                           MPI_Comm comm, MPI_Info info);
 
 #endif
 
diff --git a/src/serial_io.c b/src/serial_io.c
index 10eab97f1bf118a842e274b521056d0d81b32db1..f2c0dcce1cc258e7c9a43027fbdd45e299b47ae6 100644
--- a/src/serial_io.c
+++ b/src/serial_io.c
@@ -590,6 +590,7 @@ void read_ic_serial(char* fileName, double dim[3], struct part** parts,
  * @brief Writes an HDF5 output file (GADGET-3 type) with its XMF descriptor
  *
  * @param e The engine containing all the system.
+ * @param baseName The common part of the snapshot file name.
  * @param us The UnitSystem used for the conversion of units in the output.
  * @param mpi_rank The MPI rank of this node.
  * @param mpi_size The number of MPI ranks.
@@ -604,8 +605,9 @@ void read_ic_serial(char* fileName, double dim[3], struct part** parts,
  * Calls #error() if an error occurs.
  *
  */
-void write_output_serial(struct engine* e, struct UnitSystem* us, int mpi_rank,
-                         int mpi_size, MPI_Comm comm, MPI_Info info) {
+void write_output_serial(struct engine* e, const char* baseName,
+                         struct UnitSystem* us, int mpi_rank, int mpi_size,
+                         MPI_Comm comm, MPI_Info info) {
   hid_t h_file = 0, h_grp = 0, h_grpsph = 0;
   const size_t Ngas = e->s->nr_parts;
   const size_t Ntot = e->s->nr_gparts;
@@ -617,16 +619,13 @@ void write_output_serial(struct engine* e, struct UnitSystem* us, int mpi_rank,
   static int outputCount = 0;
   FILE* xmfFile = 0;
 
-  /* Number of particles of each type */
-  // const size_t Ndm = Ntot - Ngas;
-
-  /* MATTHIEU: Temporary fix to preserve master */
+  /* Number of unassociated gparts */
   const size_t Ndm = Ntot > 0 ? Ntot - Ngas : 0;
-  /* MATTHIEU: End temporary fix */
 
   /* File name */
   char fileName[FILENAME_BUFFER_SIZE];
-  snprintf(fileName, FILENAME_BUFFER_SIZE, "output_%03i.hdf5", outputCount);
+  snprintf(fileName, FILENAME_BUFFER_SIZE, "%s_%03i.hdf5", baseName,
+           outputCount);
 
   /* Compute offset in the file and total number of particles */
   size_t N[NUM_PARTICLE_TYPES] = {Ngas, Ndm, 0};
@@ -646,10 +645,10 @@ void write_output_serial(struct engine* e, struct UnitSystem* us, int mpi_rank,
   if (mpi_rank == 0) {
 
     /* First time, we need to create the XMF file */
-    if (outputCount == 0) createXMFfile();
+    if (outputCount == 0) createXMFfile(baseName);
 
     /* Prepare the XMF file for the new entry */
-    xmfFile = prepareXMFfile();
+    xmfFile = prepareXMFfile(baseName);
 
     /* Write the part corresponding to this specific output */
     writeXMFoutputheader(xmfFile, fileName, e->time);
diff --git a/src/serial_io.h b/src/serial_io.h
index 74ab8326dbeeb955e354687059cdd595657285f0..b7ed6eb62d823829a473f828696c291e552effa3 100644
--- a/src/serial_io.h
+++ b/src/serial_io.h
@@ -36,8 +36,9 @@ void read_ic_serial(char* fileName, double dim[3], struct part** parts,
                     int* periodic, int mpi_rank, int mpi_size, MPI_Comm comm,
                     MPI_Info info, int dry_run);
 
-void write_output_serial(struct engine* e, struct UnitSystem* us, int mpi_rank,
-                         int mpi_size, MPI_Comm comm, MPI_Info info);
+void write_output_serial(struct engine* e, const char* baseName,
+                         struct UnitSystem* us, int mpi_rank, int mpi_size,
+                         MPI_Comm comm, MPI_Info info);
 
 #endif
 
diff --git a/src/single_io.c b/src/single_io.c
index 1dc71087e102ff884dba7b7d4b6dcd6339335cac..d545fb086fa4bc63fe58dee2bfe85d9586997850 100644
--- a/src/single_io.c
+++ b/src/single_io.c
@@ -456,7 +456,8 @@ void read_ic_single(char* fileName, double dim[3], struct part** parts,
  * @brief Writes an HDF5 output file (GADGET-3 type) with its XMF descriptor
  *
  * @param e The engine containing all the system.
- * @param us The UnitSystem used for the conversion of units in the output
+ * @param baseName The common part of the snapshot file name.
+ * @param us The UnitSystem used for the conversion of units in the output.
  *
  * Creates an HDF5 output file and writes the particles contained
  * in the engine. If such a file already exists, it is erased and replaced
@@ -466,7 +467,8 @@ void read_ic_single(char* fileName, double dim[3], struct part** parts,
  * Calls #error() if an error occurs.
  *
  */
-void write_output_single(struct engine* e, struct UnitSystem* us) {
+void write_output_single(struct engine* e, const char* baseName,
+                         struct UnitSystem* us) {
 
   hid_t h_file = 0, h_grp = 0;
   const size_t Ngas = e->s->nr_parts;
@@ -478,25 +480,22 @@ void write_output_single(struct engine* e, struct UnitSystem* us) {
   struct gpart* dmparts = NULL;
   static int outputCount = 0;
 
-  /* Number of particles of each type */
-  // const size_t Ndm = Ntot - Ngas;
-
-  /* MATTHIEU: Temporary fix to preserve master */
+  /* Number of unassociated gparts */
   const size_t Ndm = Ntot > 0 ? Ntot - Ngas : 0;
-  /* MATTHIEU: End temporary fix */
 
   long long N_total[NUM_PARTICLE_TYPES] = {Ngas, Ndm, 0};
 
   /* File name */
   char fileName[FILENAME_BUFFER_SIZE];
-  snprintf(fileName, FILENAME_BUFFER_SIZE, "output_%03i.hdf5", outputCount);
+  snprintf(fileName, FILENAME_BUFFER_SIZE, "%s_%03i.hdf5", baseName,
+           outputCount);
 
   /* First time, we need to create the XMF file */
-  if (outputCount == 0) createXMFfile();
+  if (outputCount == 0) createXMFfile(baseName);
 
   /* Prepare the XMF file for the new entry */
   FILE* xmfFile = 0;
-  xmfFile = prepareXMFfile();
+  xmfFile = prepareXMFfile(baseName);
 
   /* Write the part corresponding to this specific output */
   writeXMFoutputheader(xmfFile, fileName, e->time);
diff --git a/src/single_io.h b/src/single_io.h
index 587ebe07b6fa2b984b964baf282e7ceb1003ad29..adfc5b43941b2c0d69d9ce0924164aff56864a23 100644
--- a/src/single_io.h
+++ b/src/single_io.h
@@ -30,7 +30,8 @@ void read_ic_single(char* fileName, double dim[3], struct part** parts,
                     struct gpart** gparts, size_t* Ngas, size_t* Ndm,
                     int* periodic, int dry_run);
 
-void write_output_single(struct engine* e, struct UnitSystem* us);
+void write_output_single(struct engine* e, const char* baseName,
+                         struct UnitSystem* us);
 
 #endif
 
diff --git a/src/units.c b/src/units.c
index c0276fd4fa74136f46f2f4b7d719c124c92fa48c..8e0ff08f1d92f47afbf44366acef7038fc8675c5 100644
--- a/src/units.c
+++ b/src/units.c
@@ -46,20 +46,23 @@
  *
  * @param us The UnitSystem to initialize.
  * @param params The parsed parameter file.
+ * @param category The section of the parameter file to read from.
  */
-void units_init(struct UnitSystem* us, const struct swift_params* params) {
-
-  us->UnitMass_in_cgs =
-      parser_get_param_double(params, "UnitSystem:UnitMass_in_cgs");
-  us->UnitLength_in_cgs =
-      parser_get_param_double(params, "UnitSystem:UnitLength_in_cgs");
-  const double unitVelocity =
-      parser_get_param_double(params, "UnitSystem:UnitVelocity_in_cgs");
+void units_init(struct UnitSystem* us, const struct swift_params* params,
+                const char* category) {
+
+  char buffer[200];
+  sprintf(buffer, "%s:UnitMass_in_cgs", category);
+  us->UnitMass_in_cgs = parser_get_param_double(params, buffer);
+  sprintf(buffer, "%s:UnitLength_in_cgs", category);
+  us->UnitLength_in_cgs = parser_get_param_double(params, buffer);
+  sprintf(buffer, "%s:UnitVelocity_in_cgs", category);
+  const double unitVelocity = parser_get_param_double(params, buffer);
   us->UnitTime_in_cgs = us->UnitLength_in_cgs / unitVelocity;
-  us->UnitCurrent_in_cgs =
-      parser_get_param_double(params, "UnitSystem:UnitCurrent_in_cgs");
-  us->UnitTemperature_in_cgs =
-      parser_get_param_double(params, "UnitSystem:UnitTemp_in_cgs");
+  sprintf(buffer, "%s:UnitCurrent_in_cgs", category);
+  us->UnitCurrent_in_cgs = parser_get_param_double(params, buffer);
+  sprintf(buffer, "%s:UnitTemp_in_cgs", category);
+  us->UnitTemperature_in_cgs = parser_get_param_double(params, buffer);
 }
 
 /**
diff --git a/src/units.h b/src/units.h
index 197d3e1f1dac2e01015d44758dd079ac21c6a0b7..24e37e177480d7f84e41df1b73e2036aa00b7220 100644
--- a/src/units.h
+++ b/src/units.h
@@ -92,7 +92,8 @@ enum UnitConversionFactor {
   UNIT_CONV_TEMPERATURE
 };
 
-void units_init(struct UnitSystem*, const struct swift_params*);
+void units_init(struct UnitSystem*, const struct swift_params*,
+                const char* category);
 double units_get_base_unit(const struct UnitSystem*, enum BaseUnits);
 const char* units_get_base_unit_symbol(enum BaseUnits);
 const char* units_get_base_unit_CGS_symbol(enum BaseUnits);