diff --git a/.gitignore b/.gitignore
index 4c704f0665ed0b73d3abb0ca576bfd4709e080e2..5a986acbd59a818b151540fb9303eadb4f926f77 100644
--- a/.gitignore
+++ b/.gitignore
@@ -43,6 +43,7 @@ tests/testTimeIntegration
 tests/testSPHStep
 tests/testKernel
 tests/testParser
+tests/parser_output.yml
 
 theory/latex/swift.pdf
 theory/kernel/kernels.pdf
diff --git a/README b/README
index 320df3f8ca6880776d338408c2c71ea82b1414c8..0c57e3f5656268c71bb7732af933302cbde9547b 100644
--- a/README
+++ b/README
@@ -1,4 +1,4 @@
- Welcome to the cosmological code
+ Welcome to the cosmological hydrodynamical code
     ______       _________________
    / ___/ |     / /  _/ ___/_  __/
    \__ \| | /| / // // /_   / /   
@@ -6,8 +6,26 @@
  /____/ |__/|__/___/_/    /_/     
  SPH With Inter-dependent Fine-grained Tasking
 
-Website: www.swiftsim.com
-Twitter: @SwiftSimulation
+ Website: www.swiftsim.com
+ Twitter: @SwiftSimulation
 
-See INSTALL.swift for instructions.
+See INSTALL.swift for install instructions.
 
+Usage: swift [OPTION] PARAMFILE
+
+Valid options are:
+  -c          Run with cosmological time integration
+  -d          Dry run. Read the parameter file, allocate memory but does not read 
+              the particles from ICs and exit before the start of time integration.
+              Allows user to check validy of parameter and IC files as well as memory limits.
+  -e          Enable floating-point exceptions (debugging mode)
+  -f    {int} Overwrite the CPU frequency (Hz) to be used for time measurements
+  -g          Run with an external gravitational potential
+  -G          Run with self-gravity
+  -s          Run with SPH
+  -v     [12] Increase the level of verbosity 1: MPI-rank 0 writes 
+              2: All MPI-ranks write
+  -y    {int} Time-step frequency at which task graphs are dumped
+  -h          Print this help message and exit
+
+See the file examples/parameter_example.yml for an example of parameter file.
diff --git a/configure.ac b/configure.ac
index e5d44fda300f15088c282b93b25499ecb242e24f..11ad6550d87f6764570f48449719292bcec3704d 100644
--- a/configure.ac
+++ b/configure.ac
@@ -287,11 +287,11 @@ AC_SUBST([METIS_LIBS])
 AC_SUBST([METIS_INCS])
 AM_CONDITIONAL([HAVEMETIS],[test -n "$METIS_LIBS"])
 
-# Check for zlib.
-AC_CHECK_LIB([z],[gzopen],[
-    AC_DEFINE([HAVE_LIBZ],[1],[Set to 1 if zlib is installed.])
-    LDFLAGS="$LDFLAGS -lz"
-    ],[])
+# # Check for zlib.
+# AC_CHECK_LIB([z],[gzopen],[
+#     AC_DEFINE([HAVE_LIBZ],[1],[Set to 1 if zlib is installed.])
+#     LDFLAGS="$LDFLAGS -lz"
+#     ],[])
 
 
 # Check for HDF5. This is required.
diff --git a/examples/CosmoVolume/cosmoVolume.yml b/examples/CosmoVolume/cosmoVolume.yml
new file mode 100644
index 0000000000000000000000000000000000000000..20d5febb280748a208633f75351d523b79286035
--- /dev/null
+++ b/examples/CosmoVolume/cosmoVolume.yml
@@ -0,0 +1,48 @@
+
+# Define the system of units to use internally. 
+UnitSystem:
+  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 task scheduling
+Scheduler:
+  nr_threads:       16        # The number of threads per MPI rank to use.
+  nr_queues:        0        # The number of task queues to use. Use 0  to let the system decide.
+  cell_max_size:    8000000  # Maximal number of interactions per task (this is the default value).
+  cell_sub_size:    5000     # Maximal number of interactions per sub-task  (this is the default value).
+  cell_split_size:  400      # Maximal number of particles per cell (this is the default value).
+
+# Parameters governing the time integration
+TimeIntegration:
+  time_begin: 0.    # The starting time of the simulation (in internal units).
+  time_end:   1.    # The end time of the simulation (in internal units).
+  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 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).
+  delta_neighbours:      1.       # The tolerance for the targetted number of neighbours.
+  CFL_condition:         0.1      # Courant-Friedrich-Levy condition for time integration.
+  max_ghost_iterations:  30       # Maximal number of iterations allowed to converge towards the smoothing length.
+  max_smoothing_length:  0.6      # Maximal smoothing length allowed (in internal units).
+
+# Parameters related to the initial conditions
+InitialConditions:
+  file_name:  ./cosmoVolume.hdf5     # The file to read
+  h_scaling:  1.                    # A scaling factor to apply to all smoothing lengths in the ICs.
+  shift_x:    0.                    # A shift to apply to all particles read from the ICs (in internal units).
+  shift_y:    0.
+  shift_z:    0.
+
+# Parameters govering domain decomposition
+DomainDecomposition:
+  initial_type:       m     # The initial strategy ("g", "m", "w", or "v"). See documentation for details.
+  initial_grid_x:    10     # Grid size if the 'g' strategy is chosen.
+  initial_grid_y:    10
+  initial_grid_z:    10
+  repartition_type:   b     # The re-decomposition strategy ("n", "b", "v", "e" or "x"). See documentation for details.
+ 
diff --git a/examples/CosmoVolume/run.sh b/examples/CosmoVolume/run.sh
new file mode 100755
index 0000000000000000000000000000000000000000..a788a35c76a7c0b205297a7de922a9a7e833243a
--- /dev/null
+++ b/examples/CosmoVolume/run.sh
@@ -0,0 +1,10 @@
+#!/bin/bash
+
+ # Generate the initial conditions if they are not present.
+if [ ! -e cosmoVolume.hdf5 ]
+then
+    echo "Fetching initial conditions for the cosmo volume example..."
+    ./getIC.sh
+fi
+
+../swift -s cosmoVolume.yml
diff --git a/examples/SedovBlast/run.sh b/examples/SedovBlast/run.sh
new file mode 100755
index 0000000000000000000000000000000000000000..58646cf42eecc3f31fdb8a63ca2108c02d9580ba
--- /dev/null
+++ b/examples/SedovBlast/run.sh
@@ -0,0 +1,10 @@
+#!/bin/bash
+
+# Generate the initial conditions if they are not present.
+if [ ! -e sedov.hdf5 ]
+then
+    echo "Generating initial conditions for the SedovBlast example..."
+    python makeIC_fcc.py
+fi
+
+../swift -s sedov.yml
diff --git a/examples/SedovBlast/sedov.yml b/examples/SedovBlast/sedov.yml
new file mode 100644
index 0000000000000000000000000000000000000000..f354ef5679eb5b6176ab90298bb307c6c2b27f0e
--- /dev/null
+++ b/examples/SedovBlast/sedov.yml
@@ -0,0 +1,48 @@
+
+# Define the system of units to use internally. 
+UnitSystem:
+  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 task scheduling
+Scheduler:
+  nr_threads:       16       # The number of threads per MPI rank to use.
+  nr_queues:        0        # The number of task queues to use. Use 0  to let the system decide.
+  cell_max_size:    8000000  # Maximal number of interactions per task (this is the default value).
+  cell_sub_size:    5000     # Maximal number of interactions per sub-task  (this is the default value).
+  cell_split_size:  400      # Maximal number of particles per cell (this is the default value).
+
+# Parameters governing the time integration
+TimeIntegration:
+  time_begin: 0.    # The starting time of the simulation (in internal units).
+  time_end:   1.    # The end time of the simulation (in internal units).
+  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 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).
+  delta_neighbours:      1.       # The tolerance for the targetted number of neighbours.
+  CFL_condition:         0.1      # Courant-Friedrich-Levy condition for time integration.
+  max_ghost_iterations:  30       # Maximal number of iterations allowed to converge towards the smoothing length.
+  max_smoothing_length:  1.       # Maximal smoothing length allowed (in internal units).
+
+# Parameters related to the initial conditions
+InitialConditions:
+  file_name:  ./sedov.hdf5          # The file to read
+  h_scaling:  1.                    # A scaling factor to apply to all smoothing lengths in the ICs.
+  shift_x:    0.                    # A shift to apply to all particles read from the ICs (in internal units).
+  shift_y:    0.
+  shift_z:    0.
+
+# Parameters govering domain decomposition
+DomainDecomposition:
+  initial_type:       m     # The initial strategy ("g", "m", "w", or "v"). See documentation for details.
+  initial_grid_x:    10     # Grid size if the 'g' strategy is chosen.
+  initial_grid_y:    10
+  initial_grid_z:    10
+  repartition_type:   b     # The re-decomposition strategy ("n", "b", "v", "e" or "x"). See documentation for details.
+ 
diff --git a/examples/SodShock/run.sh b/examples/SodShock/run.sh
new file mode 100755
index 0000000000000000000000000000000000000000..646f1e3a337170e2e406c24e7505e42b81de364b
--- /dev/null
+++ b/examples/SodShock/run.sh
@@ -0,0 +1,10 @@
+#!/bin/bash
+
+# Generate the initial conditions if they are not present.
+if [ ! -e sodShock.hdf5 ]
+then
+    echo "Generating initial conditions for the SodShock example..."
+    python makeIC.py
+fi
+
+../swift -s sodShock.yml
diff --git a/examples/SodShock/sodShock.yml b/examples/SodShock/sodShock.yml
new file mode 100644
index 0000000000000000000000000000000000000000..5fe7be7b9fc13bb5bc67556d79d8ff9d9eff81d9
--- /dev/null
+++ b/examples/SodShock/sodShock.yml
@@ -0,0 +1,48 @@
+
+# Define the system of units to use internally. 
+UnitSystem:
+  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 task scheduling
+Scheduler:
+  nr_threads:       16        # The number of threads per MPI rank to use.
+  nr_queues:        0        # The number of task queues to use. Use 0  to let the system decide.
+  cell_max_size:    8000000  # Maximal number of interactions per task (this is the default value).
+  cell_sub_size:    5000     # Maximal number of interactions per sub-task.
+  cell_split_size:  400      # Maximal number of particles per cell (this is the default value).
+
+# Parameters governing the time integration
+TimeIntegration:
+  time_begin: 0.    # The starting time of the simulation (in internal units).
+  time_end:   1.    # The end time of the simulation (in internal units).
+  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 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).
+  delta_neighbours:      1.       # The tolerance for the targetted number of neighbours.
+  CFL_condition:         0.1      # Courant-Friedrich-Levy condition for time integration.
+  max_ghost_iterations:  30       # Maximal number of iterations allowed to converge towards the smoothing length.
+  max_smoothing_length:  0.01     # Maximal smoothing length allowed (in internal units).
+
+# Parameters related to the initial conditions
+InitialConditions:
+  file_name:  ./sodShock.hdf5       # The file to read
+  h_scaling:  1.                    # A scaling factor to apply to all smoothing lengths in the ICs.
+  shift_x:    0.                    # A shift to apply to all particles read from the ICs (in internal units).
+  shift_y:    0.
+  shift_z:    0.
+
+# Parameters govering domain decomposition
+DomainDecomposition:
+  initial_type:       m     # The initial strategy ("g", "m", "w", or "v"). See documentation for details.
+  initial_grid_x:    10     # Grid size if the 'g' strategy is chosen.
+  initial_grid_y:    10
+  initial_grid_z:    10
+  repartition_type:   b     # The re-decomposition strategy ("n", "b", "v", "e" or "x"). See documentation for details.
+ 
diff --git a/examples/UniformBox/run.sh b/examples/UniformBox/run.sh
new file mode 100755
index 0000000000000000000000000000000000000000..ca78b0ac0425bf1b3f6dd9d30bfc95d35083739f
--- /dev/null
+++ b/examples/UniformBox/run.sh
@@ -0,0 +1,10 @@
+#!/bin/bash
+
+# Generate the initial conditions if they are not present.
+if [ ! -e uniformBox.hdf5 ]
+then
+    echo "Generating initial conditions for the uniform box example..."
+    python makeIC.py 100
+fi
+
+../swift -s uniformBox.yml
diff --git a/examples/UniformBox/uniformBox.yml b/examples/UniformBox/uniformBox.yml
new file mode 100644
index 0000000000000000000000000000000000000000..2d5512815b60511b5dbc373df43fae4658272093
--- /dev/null
+++ b/examples/UniformBox/uniformBox.yml
@@ -0,0 +1,48 @@
+
+# Define the system of units to use internally. 
+UnitSystem:
+  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 task scheduling
+Scheduler:
+  nr_threads:       16        # The number of threads per MPI rank to use.
+  nr_queues:        0        # The number of task queues to use. Use 0  to let the system decide.
+  cell_max_size:    8000000  # Maximal number of interactions per task (this is the default value).
+  cell_sub_size:    5000     # Maximal number of interactions per sub-task  (this is the default value).
+  cell_split_size:  400      # Maximal number of particles per cell (this is the default value).
+
+# Parameters governing the time integration
+TimeIntegration:
+  time_begin: 0.    # The starting time of the simulation (in internal units).
+  time_end:   1.    # The end time of the simulation (in internal units).
+  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 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).
+  delta_neighbours:      1.       # The tolerance for the targetted number of neighbours.
+  CFL_condition:         0.1      # Courant-Friedrich-Levy condition for time integration.
+  max_ghost_iterations:  30       # Maximal number of iterations allowed to converge towards the smoothing length.
+  max_smoothing_length:  0.1      # Maximal smoothing length allowed (in internal units).
+
+# Parameters related to the initial conditions
+InitialConditions:
+  file_name:  ./uniformBox.hdf5     # The file to read
+  h_scaling:  1.                    # A scaling factor to apply to all smoothing lengths in the ICs.
+  shift_x:    0.                    # A shift to apply to all particles read from the ICs (in internal units).
+  shift_y:    0.
+  shift_z:    0.
+
+# Parameters govering domain decomposition
+DomainDecomposition:
+  initial_type:       m     # The initial strategy ("g", "m", "w", or "v"). See documentation for details.
+  initial_grid_x:    10     # Grid size if the 'g' strategy is chosen.
+  initial_grid_y:    10
+  initial_grid_z:    10
+  repartition_type:   b     # The re-decomposition strategy ("n", "b", "v", "e" or "x"). See documentation for details.
+ 
diff --git a/examples/main.c b/examples/main.c
index 9523af49ed30c54d256d287ea2846a854650bc05..5cfae5efba9157ba7b727115b03ac467287edc3d 100644
--- a/examples/main.c
+++ b/examples/main.c
@@ -23,20 +23,11 @@
 #include "../config.h"
 
 /* Some standard headers. */
+#include <fenv.h>
+#include <unistd.h>
 #include <stdio.h>
 #include <stdlib.h>
-#include <unistd.h>
 #include <string.h>
-#include <pthread.h>
-#include <math.h>
-#include <float.h>
-#include <limits.h>
-#include <fenv.h>
-
-/* Conditional headers. */
-#ifdef HAVE_LIBZ
-#include <zlib.h>
-#endif
 
 /* MPI headers. */
 #ifdef WITH_MPI
@@ -51,59 +42,51 @@
 #define ENGINE_POLICY engine_policy_none
 #endif
 
+/**
+ * @brief Help messages for the command line parameters.
+ */
+void print_help_message() {
+
+  printf("\nUsage: swift [OPTION] PARAMFILE\n\n");
+  printf("Valid options are:\n");
+  printf("  %2s %8s %s\n", "-c", "", "Run with cosmological time integration");
+  printf(
+      "  %2s %8s %s\n", "-d", "",
+      "Dry run. Read the parameter file, allocate memory but does not read ");
+  printf(
+      "  %2s %8s %s\n", "", "",
+      "the particles from ICs and exit before the start of time integration.");
+  printf("  %2s %8s %s\n", "", "",
+         "Allows user to check validy of parameter and IC files as well as "
+         "memory limits.");
+  printf("  %2s %8s %s\n", "-e", "",
+         "Enable floating-point exceptions (debugging mode)");
+  printf("  %2s %8s %s\n", "-f", "{int}",
+         "Overwrite the CPU frequency (Hz) to be used for time measurements");
+  printf("  %2s %8s %s\n", "-g", "",
+         "Run with an external gravitational potential");
+  printf("  %2s %8s %s\n", "-G", "", "Run with self-gravity");
+  printf("  %2s %8s %s\n", "-s", "", "Run with SPH");
+  printf("  %2s %8s %s\n", "-v", "[12]",
+         "Increase the level of verbosity 1: MPI-rank 0 writes ");
+  printf("  %2s %8s %s\n", "", "", "2: All MPI-ranks write");
+  printf("  %2s %8s %s\n", "-y", "{int}",
+         "Time-step frequency at which task graphs are dumped");
+  printf("  %2s %8s %s\n", "-h", "", "Print this help message and exit");
+  printf(
+      "\nSee the file parameter_example.yml for an example of "
+      "parameter file.\n");
+}
+
 /**
  * @brief Main routine that loads a few particles and generates some output.
  *
  */
 int main(int argc, char *argv[]) {
 
-  int c, icount, periodic = 1;
-  size_t Ngas = 0, Ngpart = 0;
-  long long N_total[2] = {0, 0};
-  int nr_threads = 1, nr_queues = -1;
-  int dump_tasks = 0;
-  int data[2];
-  double dim[3] = {1.0, 1.0, 1.0}, shift[3] = {0.0, 0.0, 0.0};
-  double h_max = -1.0, scaling = 1.0;
-  double time_end = DBL_MAX;
-  struct part *parts = NULL;
-  struct gpart *gparts = NULL;
-  struct space s;
-  struct engine e;
-  struct UnitSystem us;
   struct clocks_time tic, toc;
-  char ICfileName[200] = "";
-  char dumpfile[30];
-  float dt_max = 0.0f, dt_min = 0.0f;
-  int nr_nodes = 1, myrank = 0;
-  FILE *file_thread;
-  int with_outputs = 1;
-  int with_external_gravity = 0;
-  int with_self_gravity = 0;
-  int engine_policies = 0;
-  int verbose = 0, talking = 0;
-  unsigned long long cpufreq = 0;
 
-#ifdef WITH_MPI
-  struct partition initial_partition;
-  enum repartition_type reparttype = REPART_NONE;
-
-  initial_partition.type = INITPART_GRID;
-  initial_partition.grid[0] = 1;
-  initial_partition.grid[1] = 1;
-  initial_partition.grid[2] = 1;
-#ifdef HAVE_METIS
-  /* Defaults make use of METIS. */
-  reparttype = REPART_METIS_BOTH;
-  initial_partition.type = INITPART_METIS_NOWEIGHT;
-#endif
-#endif
-
-  /* Choke on FP-exceptions. */
-  // feenableexcept( FE_DIVBYZERO | FE_INVALID | FE_OVERFLOW );
-
-  /* Initialize CPU frequency, this also starts time. */
-  clocks_set_cpufreq(cpufreq);
+  int nr_nodes = 1, myrank = 0;
 
 #ifdef WITH_MPI
   /* Start by initializing MPI. */
@@ -122,70 +105,63 @@ int main(int argc, char *argv[]) {
   if ((res = MPI_Comm_set_errhandler(MPI_COMM_WORLD, MPI_ERRORS_RETURN)) !=
       MPI_SUCCESS)
     error("Call to MPI_Comm_set_errhandler failed with error %i.", res);
-  if (myrank == 0) message("MPI is up and running with %i node(s).", nr_nodes);
+  if (myrank == 0)
+    printf("[0000][00000.0] MPI is up and running with %i node(s).\n",
+           nr_nodes);
+  if (nr_nodes == 1) {
+    message("WARNING: you are running with one MPI rank.");
+    message("WARNING: you should use the non-MPI version of this program.");
+  }
   fflush(stdout);
-
-  /* Set a default grid so that grid[0]*grid[1]*grid[2] == nr_nodes. */
-  factor(nr_nodes, &initial_partition.grid[0], &initial_partition.grid[1]);
-  factor(nr_nodes / initial_partition.grid[1], &initial_partition.grid[0],
-         &initial_partition.grid[2]);
-  factor(initial_partition.grid[0] * initial_partition.grid[1],
-         &initial_partition.grid[1], &initial_partition.grid[0]);
 #endif
 
-  /* Greeting message */
-  if (myrank == 0) greetings();
-
 #if defined(HAVE_SETAFFINITY) && defined(HAVE_LIBNUMA) && defined(_GNU_SOURCE)
   if ((ENGINE_POLICY) & engine_policy_setaffinity) {
     /* Ensure the NUMA node on which we initialise (first touch) everything
      * doesn't change before engine_init allocates NUMA-local workers.
-     * Otherwise,
-     * we may be scheduled elsewhere between the two times.
+     * Otherwise, we may be scheduled elsewhere between the two times.
      */
     cpu_set_t affinity;
     CPU_ZERO(&affinity);
     CPU_SET(sched_getcpu(), &affinity);
     if (sched_setaffinity(0, sizeof(cpu_set_t), &affinity) != 0) {
-      message("failed to set entry thread's affinity");
-    } else {
-      message("set entry thread's affinity");
+      error("failed to set entry thread's affinity");
     }
   }
 #endif
 
-  /* Init the space. */
-  bzero(&s, sizeof(struct space));
+  /* Welcome to SWIFT, you made the right choice */
+  if (myrank == 0) greetings();
 
-  /* Parse the options */
-  while ((c = getopt(argc, argv, "a:c:d:e:f:gGh:m:oP:q:R:s:t:v:w:y:z:")) != -1)
-    switch (c) {
-      case 'a':
-        if (sscanf(optarg, "%lf", &scaling) != 1)
-          error("Error parsing cutoff scaling.");
-        if (myrank == 0) message("scaling cutoff by %.3f.", scaling);
-        fflush(stdout);
-        break;
+  int dry_run = 0;
+  int dump_tasks = 0;
+  int with_cosmology = 0;
+  int with_external_gravity = 0;
+  int with_self_gravity = 0;
+  int with_hydro = 0;
+  int with_fp_exceptions = 0;
+  int verbose = 0;
+  char paramFileName[200] = "";
+  unsigned long long cpufreq = 0;
+
+  /* Parse the parameters */
+  int c;
+  while ((c = getopt(argc, argv, "cdef:gGhsv:y")) != -1) switch (c) {
       case 'c':
-        if (sscanf(optarg, "%lf", &time_end) != 1)
-          error("Error parsing final time.");
-        if (myrank == 0) message("time_end set to %.3e.", time_end);
-        fflush(stdout);
+        with_cosmology = 1;
         break;
       case 'd':
-        if (sscanf(optarg, "%f", &dt_min) != 1)
-          error("Error parsing minimal timestep.");
-        if (myrank == 0) message("dt_min set to %e.", dt_min);
-        fflush(stdout);
+        dry_run = 1;
         break;
       case 'e':
-        if (sscanf(optarg, "%f", &dt_max) != 1)
-          error("Error parsing maximal timestep.");
-        if (myrank == 0) message("dt_max set to %e.", dt_max);
-        fflush(stdout);
+        with_fp_exceptions = 1;
         break;
       case 'f':
-        if (!strcpy(ICfileName, optarg)) error("Error parsing IC file name.");
+        if (sscanf(optarg, "%llu", &cpufreq) != 1) {
+          if (myrank == 0) printf("Error parsing CPU frequency (-f).\n");
+          if (myrank == 0) print_help_message();
+          return 1;
+        }
         break;
       case 'g':
         with_external_gravity = 1;
@@ -194,257 +170,177 @@ int main(int argc, char *argv[]) {
         with_self_gravity = 1;
         break;
       case 'h':
-        if (sscanf(optarg, "%llu", &cpufreq) != 1)
-          error("Error parsing CPU frequency.");
-        if (myrank == 0) message("CPU frequency set to %llu.", cpufreq);
-        fflush(stdout);
-        break;
-      case 'm':
-        if (sscanf(optarg, "%lf", &h_max) != 1) error("Error parsing h_max.");
-        if (myrank == 0) message("maximum h set to %e.", h_max);
-        fflush(stdout);
-        break;
-      case 'o':
-        with_outputs = 0;
-        break;
-      case 'P':
-/* Partition type is one of "g", "m", "w", or "v"; "g" can be
- * followed by three numbers defining the grid. */
-#ifdef WITH_MPI
-        switch (optarg[0]) {
-          case 'g':
-            initial_partition.type = INITPART_GRID;
-            if (strlen(optarg) > 2) {
-              if (sscanf(optarg, "g %i %i %i", &initial_partition.grid[0],
-                         &initial_partition.grid[1],
-                         &initial_partition.grid[2]) != 3)
-                error("Error parsing grid.");
-            }
-            break;
-#ifdef HAVE_METIS
-          case 'm':
-            initial_partition.type = INITPART_METIS_NOWEIGHT;
-            break;
-          case 'w':
-            initial_partition.type = INITPART_METIS_WEIGHT;
-            break;
-#endif
-          case 'v':
-            initial_partition.type = INITPART_VECTORIZE;
-            break;
-        }
-#endif
-        break;
-      case 'q':
-        if (sscanf(optarg, "%d", &nr_queues) != 1)
-          error("Error parsing number of queues.");
-        break;
-      case 'R':
-/* Repartition type "n", "b", "v", "e" or "x".
- * Note only none is available without METIS. */
-#ifdef WITH_MPI
-        switch (optarg[0]) {
-          case 'n':
-            reparttype = REPART_NONE;
-            break;
-#ifdef HAVE_METIS
-          case 'b':
-            reparttype = REPART_METIS_BOTH;
-            break;
-          case 'v':
-            reparttype = REPART_METIS_VERTEX;
-            break;
-          case 'e':
-            reparttype = REPART_METIS_EDGE;
-            break;
-          case 'x':
-            reparttype = REPART_METIS_VERTEX_EDGE;
-            break;
-#endif
-        }
-#endif
-        break;
+        if (myrank == 0) print_help_message();
+        return 0;
       case 's':
-        if (sscanf(optarg, "%lf %lf %lf", &shift[0], &shift[1], &shift[2]) != 3)
-          error("Error parsing shift.");
-        if (myrank == 0)
-          message("will shift parts by [ %.3f %.3f %.3f ].", shift[0], shift[1],
-                  shift[2]);
-        break;
-      case 't':
-        if (sscanf(optarg, "%d", &nr_threads) != 1)
-          error("Error parsing number of threads.");
+        with_hydro = 1;
         break;
       case 'v':
-        /* verbose = 1: MPI rank 0 writes
-           verbose = 2: all MPI ranks write */
-        if (sscanf(optarg, "%d", &verbose) != 1)
-          error("Error parsing verbosity level.");
-        break;
-      case 'w':
-        if (sscanf(optarg, "%d", &space_subsize) != 1)
-          error("Error parsing sub size.");
-        if (myrank == 0) message("sub size set to %i.", space_subsize);
+        if (sscanf(optarg, "%d", &verbose) != 1) {
+          if (myrank == 0) printf("Error parsing verbosity level (-v).\n");
+          if (myrank == 0) print_help_message();
+          return 1;
+        }
         break;
       case 'y':
-        if (sscanf(optarg, "%d", &dump_tasks) != 1)
-          error("Error parsing dump_tasks (-y)");
-        break;
-      case 'z':
-        if (sscanf(optarg, "%d", &space_splitsize) != 1)
-          error("Error parsing split size.");
-        if (myrank == 0) message("split size set to %i.", space_splitsize);
+        if (sscanf(optarg, "%d", &dump_tasks) != 1) {
+          if (myrank == 0) printf("Error parsing dump_tasks (-y). \n");
+          if (myrank == 0) print_help_message();
+          return 1;
+        }
         break;
       case '?':
-        error("Unknown option.");
+        if (myrank == 0) print_help_message();
+        return 1;
         break;
     }
+  if (optind == argc - 1) {
+    if (!strcpy(paramFileName, argv[optind++]))
+      error("Error reading parameter file name.");
+  } else if (optind > argc - 1) {
+    if (myrank == 0) printf("Error: A parameter file name must be provided\n");
+    if (myrank == 0) print_help_message();
+    return 1;
+  } else {
+    if (myrank == 0) printf("Error: Too many parameters given\n");
+    if (myrank == 0) print_help_message();
+    return 1;
+  }
+  if (!with_self_gravity && !with_hydro && !with_external_gravity) {
+    if (myrank == 0)
+      printf("Error: At least one of -s, -g or -G must be chosen.\n");
+    if (myrank == 0) print_help_message();
+    return 1;
+  }
 
-#ifdef WITH_MPI
+  /* Genesis 1.1: And then, there was time ! */
+  clocks_set_cpufreq(cpufreq);
+
+  if (myrank == 0 && dry_run)
+    message(
+        "Executing a dry run. No i/o or time integration will be performed.");
+
+  /* Report CPU frequency. */
   if (myrank == 0) {
-    message("Running with %i thread(s) per node.", nr_threads);
-    message("Using initial partition %s",
-            initial_partition_name[initial_partition.type]);
-    if (initial_partition.type == INITPART_GRID)
-      message("grid set to [ %i %i %i ].", initial_partition.grid[0],
-              initial_partition.grid[1], initial_partition.grid[2]);
-    message("Using %s repartitioning", repartition_name[reparttype]);
+    cpufreq = clocks_get_cpufreq();
+    message("CPU frequency used for tick conversion: %llu Hz", cpufreq);
+  }
 
-    if (nr_nodes == 1) {
-      message("WARNING: you are running with one MPI rank.");
-      message("WARNING: you should use the non-MPI version of this program.");
-    }
-    fflush(stdout);
+  /* Do we choke on FP-exceptions ? */
+  if (with_fp_exceptions) {
+    feenableexcept(FE_DIVBYZERO | FE_INVALID | FE_OVERFLOW);
+    if (myrank == 0) message("Floating point exceptions will be reported.");
   }
-#else
-  if (myrank == 0) message("Running with %i thread(s).", nr_threads);
-#endif
 
   /* How large are the parts? */
   if (myrank == 0) {
-    message("sizeof(struct part) is %li bytes.", (long int)sizeof(struct part));
-    message("sizeof(struct xpart) is %li bytes.",
-            (long int)sizeof(struct xpart));
-    message("sizeof(struct gpart) is %li bytes.",
-            (long int)sizeof(struct gpart));
+    message("sizeof(struct part)  is %4zi bytes.", sizeof(struct part));
+    message("sizeof(struct xpart) is %4zi bytes.", sizeof(struct xpart));
+    message("sizeof(struct gpart) is %4zi bytes.", sizeof(struct gpart));
   }
 
+  /* How vocal are we ? */
+  const int talking = (verbose == 1 && myrank == 0) || (verbose == 2);
+
+  /* Read the parameter file */
+  struct swift_params *params = malloc(sizeof(struct swift_params));
+  if (params == NULL) error("Error allocating memory for the parameter file.");
+  if (myrank == 0) {
+    message("Reading parameters from file '%s'", paramFileName);
+    parser_read_file(paramFileName, params);
+    // parser_print_params(&params);
+    parser_write_params_to_file(params, "used_parameters.yml");
+  }
+#ifdef WITH_MPI
+  /* Broadcast the parameter file */
+  MPI_Bcast(params, sizeof(struct swift_params), MPI_BYTE, 0, MPI_COMM_WORLD);
+#endif
+
   /* Initialize unit system */
-  initUnitSystem(&us);
+  struct UnitSystem us;
+  units_init(&us, params);
   if (myrank == 0) {
     message("Unit system: U_M = %e g.", us.UnitMass_in_cgs);
     message("Unit system: U_L = %e cm.", us.UnitLength_in_cgs);
     message("Unit system: U_t = %e s.", us.UnitTime_in_cgs);
     message("Unit system: U_I = %e A.", us.UnitCurrent_in_cgs);
     message("Unit system: U_T = %e K.", us.UnitTemperature_in_cgs);
-    message("Density units: %e a^%f h^%f.",
-            conversionFactor(&us, UNIT_CONV_DENSITY),
-            aFactor(&us, UNIT_CONV_DENSITY), hFactor(&us, UNIT_CONV_DENSITY));
-    message("Entropy units: %e a^%f h^%f.",
-            conversionFactor(&us, UNIT_CONV_ENTROPY),
-            aFactor(&us, UNIT_CONV_ENTROPY), hFactor(&us, UNIT_CONV_ENTROPY));
   }
 
-  /* Report CPU frequency. */
+/* Prepare the domain decomposition scheme */
+#ifdef WITH_MPI
+  struct partition initial_partition;
+  enum repartition_type reparttype;
+  partition_init(&initial_partition, &reparttype, params, nr_nodes);
+
+  /* Let's report what we did */
   if (myrank == 0) {
-    cpufreq = clocks_get_cpufreq();
-    message("CPU frequency used for tick conversion: %llu Hz", cpufreq);
+    message("Using initial partition %s",
+            initial_partition_name[initial_partition.type]);
+    if (initial_partition.type == INITPART_GRID)
+      message("grid set to [ %i %i %i ].", initial_partition.grid[0],
+              initial_partition.grid[1], initial_partition.grid[2]);
+    message("Using %s repartitioning", repartition_name[reparttype]);
   }
+#endif
 
-  /* Check whether an IC file has been provided */
-  if (strcmp(ICfileName, "") == 0)
-    error("An IC file name must be provided via the option -f");
-
-  /* Read particles and space information from (GADGET) IC */
-
+  /* Read particles and space information from (GADGET) ICs */
+  char ICfileName[200] = "";
+  parser_get_param_string(params, "InitialConditions:file_name", ICfileName);
+  if (myrank == 0) message("Reading ICs from file '%s'", ICfileName);
+  struct part *parts = NULL;
+  struct gpart *gparts = NULL;
+  size_t Ngas = 0, Ngpart = 0;
+  double dim[3] = {0., 0., 0.};
+  int periodic = 0;
   if (myrank == 0) clocks_gettime(&tic);
 #if defined(WITH_MPI)
 #if defined(HAVE_PARALLEL_HDF5)
   read_ic_parallel(ICfileName, dim, &parts, &gparts, &Ngas, &Ngpart, &periodic,
-                   myrank, nr_nodes, MPI_COMM_WORLD, MPI_INFO_NULL);
+                   myrank, nr_nodes, MPI_COMM_WORLD, MPI_INFO_NULL, dry_run);
 #else
   read_ic_serial(ICfileName, dim, &parts, &gparts, &Ngas, &Ngpart, &periodic,
-                 myrank, nr_nodes, MPI_COMM_WORLD, MPI_INFO_NULL);
+                 myrank, nr_nodes, MPI_COMM_WORLD, MPI_INFO_NULL, dry_run);
 #endif
 #else
-  read_ic_single(ICfileName, dim, &parts, &gparts, &Ngas, &Ngpart, &periodic);
+  read_ic_single(ICfileName, dim, &parts, &gparts, &Ngas, &Ngpart, &periodic,
+                 dry_run);
 #endif
-
   if (myrank == 0) {
     clocks_gettime(&toc);
-    message("reading particle properties took %.3f %s.",
-            clocks_diff(&tic, &toc), clocks_getunit());
+    message("Reading initial conditions took %.3f %s.", clocks_diff(&tic, &toc),
+            clocks_getunit());
     fflush(stdout);
   }
 
-#if defined(WITH_MPI)
-  long long N_long[2] = {Ngas, Ngpart};
-  MPI_Reduce(&N_long, &N_total, 2, MPI_LONG_LONG, MPI_SUM, 0, MPI_COMM_WORLD);
-  N_total[1] -= N_total[0];
-  if (myrank == 0)
-    message("Read %lld gas particles and %lld DM particles from the ICs",
-            N_total[0], N_total[1]);
-#else
-  N_total[0] = Ngas;
-  N_total[1] = Ngpart - Ngas;
-  message("Read %lld gas particles and %lld DM particles from the ICs",
-          N_total[0], N_total[1]);
-#endif
-
-  /* MATTHIEU: Temporary fix to preserve master */
+  /* Discard gparts if we don't have gravity
+   * (Better implementation of i/o will come)*/
   if (!with_external_gravity && !with_self_gravity) {
     free(gparts);
     gparts = NULL;
     for (size_t k = 0; k < Ngas; ++k) parts[k].gpart = NULL;
     Ngpart = 0;
+  }
+
+  /* Get the total number of particles across all nodes. */
+  long long N_total[2] = {0, 0};
 #if defined(WITH_MPI)
-    N_long[0] = Ngas;
-    N_long[1] = Ngpart;
-    MPI_Reduce(&N_long, &N_total, 2, MPI_LONG_LONG, MPI_SUM, 0, MPI_COMM_WORLD);
-    if (myrank == 0)
-      message(
-          "AFTER FIX: Read %lld gas particles and %lld DM particles from the "
-          "ICs",
-          N_total[0], N_total[1]);
+  long long N_long[2] = {Ngas, Ngpart};
+  MPI_Reduce(&N_long, &N_total, 2, MPI_LONG_LONG, MPI_SUM, 0, MPI_COMM_WORLD);
 #else
-    N_total[0] = Ngas;
-    N_total[1] = Ngpart;
-    message(
-        "AFTER FIX: Read %lld gas particles and %lld DM particles from the ICs",
-        N_total[0], N_total[1]);
+  N_total[0] = Ngas;
+  N_total[1] = Ngpart;
 #endif
-  }
-  /* MATTHIEU: End temporary fix */
-
-  /* Apply h scaling */
-  if (scaling != 1.0)
-    for (size_t k = 0; k < Ngas; k++) parts[k].h *= scaling;
-
-  /* Apply shift */
-  if (shift[0] != 0 || shift[1] != 0 || shift[2] != 0) {
-    for (size_t k = 0; k < Ngas; k++) {
-      parts[k].x[0] += shift[0];
-      parts[k].x[1] += shift[1];
-      parts[k].x[2] += shift[2];
-    }
-    for (size_t k = 0; k < Ngpart; k++) {
-      gparts[k].x[0] += shift[0];
-      gparts[k].x[1] += shift[1];
-      gparts[k].x[2] += shift[2];
-    }
-  }
-
-  /* Set default number of queues. */
-  if (nr_queues < 0) nr_queues = nr_threads;
-
-  /* How vocal are we ? */
-  talking = (verbose == 1 && myrank == 0) || (verbose == 2);
+  if (myrank == 0)
+    message("Read %lld gas particles and %lld gparts from the ICs.", N_total[0],
+            N_total[1]);
 
-  /* Initialize the space with this data. */
+  /* Initialize the space with these data. */
   if (myrank == 0) clocks_gettime(&tic);
-  space_init(&s, dim, parts, gparts, Ngas, Ngpart, periodic, h_max,
-             myrank == 0);
-  if (myrank == 0 && verbose) {
+  struct space s;
+  space_init(&s, params, dim, parts, gparts, Ngas, Ngpart, periodic, talking,
+             dry_run);
+  if (myrank == 0) {
     clocks_gettime(&toc);
     message("space_init took %.3f %s.", clocks_diff(&tic, &toc),
             clocks_getunit());
@@ -461,49 +357,47 @@ int main(int argc, char *argv[]) {
     message("%zi parts in %i cells.", s.nr_parts, s.tot_cells);
     message("%zi gparts in %i cells.", s.nr_gparts, s.tot_cells);
     message("maximum depth is %d.", s.maxdepth);
-    // message( "cutoffs in [ %g %g ]." , s.h_min , s.h_max ); fflush(stdout);
   }
 
   /* Verify that each particle is in it's proper cell. */
-  if (myrank == 0) {
-    icount = 0;
+  if (talking && !dry_run) {
+    int icount = 0;
     space_map_cells_pre(&s, 0, &map_cellcheck, &icount);
     message("map_cellcheck picked up %i parts.", icount);
   }
 
-  if (myrank == 0) {
-    data[0] = s.maxdepth;
-    data[1] = 0;
+  /* Verify the maximal depth of cells. */
+  if (talking && !dry_run) {
+    int data[2] = {s.maxdepth, 0};
     space_map_cells_pre(&s, 0, &map_maxdepth, data);
     message("nr of cells at depth %i is %i.", data[0], data[1]);
   }
 
   /* Construct the engine policy */
-  engine_policies = ENGINE_POLICY | engine_policy_steal | engine_policy_hydro;
-  if (with_external_gravity) engine_policies |= engine_policy_external_gravity;
+  int engine_policies = ENGINE_POLICY | engine_policy_steal;
+  if (with_hydro) engine_policies |= engine_policy_hydro;
   if (with_self_gravity) engine_policies |= engine_policy_self_gravity;
+  if (with_external_gravity) engine_policies |= engine_policy_external_gravity;
+  if (with_cosmology) engine_policies |= engine_policy_cosmology;
 
-  /* Initialize the engine with this space. */
+  /* Initialize the engine with the space and policies. */
   if (myrank == 0) clocks_gettime(&tic);
-  if (myrank == 0) message("nr_nodes is %i.", nr_nodes);
-  engine_init(&e, &s, dt_max, nr_threads, nr_queues, nr_nodes, myrank,
-              engine_policies, 0, time_end, dt_min, dt_max, talking);
-  if (myrank == 0 && verbose) {
+  struct engine e;
+  engine_init(&e, &s, params, nr_nodes, myrank, engine_policies, talking);
+  if (myrank == 0) {
     clocks_gettime(&toc);
     message("engine_init took %.3f %s.", clocks_diff(&tic, &toc),
             clocks_getunit());
     fflush(stdout);
   }
 
-#ifdef WITH_MPI
-  /* Split the space. */
-  engine_split(&e, &initial_partition);
-  engine_redistribute(&e);
-#endif
+  /* Now that everything is ready, no need for the parameters any more */
+  free(params);
+  params = NULL;
 
-  if (with_outputs) {
-    /* Write the state of the system as it is before starting time integration.
-     */
+  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)
@@ -529,26 +423,42 @@ int main(int argc, char *argv[]) {
   for (k = 0; k < runner_hist_N; k++) runner_hist_bins[k] = 0;
 #endif
 
+  /* 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], time_end, e.nr_threads, e.sched.nr_queues,
+        N_total[0], N_total[1], e.timeEnd, e.nr_threads, e.sched.nr_queues,
         e.dt_min, e.dt_max);
     fflush(stdout);
   }
 
+  /* Time to say good-bye if this was not a serious run. */
+  if (dry_run) {
+#ifdef WITH_MPI
+    if ((res = MPI_Finalize()) != MPI_SUCCESS)
+      error("call to MPI_Finalize failed with error %i.", res);
+#endif
+    if (myrank == 0)
+      message("Time integration ready to start. End of dry-run.");
+    return 0;
+  }
+
+#ifdef WITH_MPI
+  /* Split the space. */
+  engine_split(&e, &initial_partition);
+  engine_redistribute(&e);
+#endif
+
   /* Initialise the particles */
   engine_init_particles(&e);
 
   /* Legend */
   if (myrank == 0)
-    printf(
-        "# Step  Time  time-step  Number of updates  Number of updates "
-        "CPU Wall-clock time [%s]\n",
-        clocks_getunit());
+    printf("# %6s %14s %14s %10s %10s %16s [%s]\n", "Step", "Time", "Time-step",
+           "Updates", "g-Updates", "Wall-clock time", clocks_getunit());
 
-  /* Let loose a runner on the space. */
+  /* Main simulation loop */
   for (int j = 0; !engine_is_done(&e); j++) {
 
 /* Repartition the space amongst the nodes? */
@@ -591,7 +501,9 @@ int main(int argc, char *argv[]) {
 #ifdef WITH_MPI
 
       /* Make sure output file is empty, only on one rank. */
-      sprintf(dumpfile, "thread_info_MPI-step%d.dat", j);
+      char dumpfile[30];
+      snprintf(dumpfile, 30, "thread_info_MPI-step%d.dat", j);
+      FILE *file_thread;
       if (myrank == 0) {
         file_thread = fopen(dumpfile, "w");
         fclose(file_thread);
@@ -636,7 +548,9 @@ int main(int argc, char *argv[]) {
       }
 
 #else
-      sprintf(dumpfile, "thread_info-step%d.dat", j);
+      char dumpfile[30];
+      snprintf(dumpfile, 30, "thread_info-step%d.dat", j);
+      FILE *file_thread;
       file_thread = fopen(dumpfile, "w");
       for (int l = 0; l < e.sched.nr_tasks; l++)
         if (!e.sched.tasks[l].skip && !e.sched.tasks[l].implicit)
@@ -650,26 +564,6 @@ int main(int argc, char *argv[]) {
       fclose(file_thread);
 #endif
     }
-
-    /* Dump a line of aggregate output. */
-    /*     if (myrank == 0) { */
-    /*       printf("%i %e %.16e %.16e %.16e %.3e %.3e %i %.3e %.3e", j, e.time,
-     */
-    /*              e.ekin + e.epot, e.ekin, e.epot, e.dt, e.dt_step,
-     * e.count_step, */
-    /*              e.dt_min, e.dt_max); */
-    /*       for (k = 0; k < timer_count; k++) */
-    /*         printf(" %.3f", clocks_from_ticks(timers[k]); */
-    /*       printf("\n"); */
-    /*       fflush(stdout); */
-    /*     } */
-
-    /* if (myrank == 0) { */
-    /*   printf("%i %e", j, e.time); */
-    /*   printf(" %.3f", clocks_from_ticks(timers[timer_count - 1]); */
-    /*   printf("\n"); */
-    /*   fflush(stdout); */
-    /* } */
   }
 
 /* Print the values of the runner histogram. */
@@ -707,7 +601,7 @@ int main(int argc, char *argv[]) {
   }
 
 #ifdef WITH_MPI
-  if (MPI_Finalize() != MPI_SUCCESS)
+  if ((res = MPI_Finalize()) != MPI_SUCCESS)
     error("call to MPI_Finalize failed with error %i.", res);
 #endif
 
diff --git a/examples/parameter_example.yml b/examples/parameter_example.yml
new file mode 100644
index 0000000000000000000000000000000000000000..b91e99baf383a399b72bfb73f1791ab7ac6f3d91
--- /dev/null
+++ b/examples/parameter_example.yml
@@ -0,0 +1,48 @@
+
+# Define the system of units to use internally. 
+UnitSystem:
+  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 task scheduling
+Scheduler:
+  nr_threads:       2        # The number of threads per MPI rank to use.
+  nr_queues:        0        # The number of task queues to use. Use 0  to let the system decide.
+  cell_max_size:    8000000  # Maximal number of interactions per task (this is the default value).
+  cell_sub_size:    8000000  # Maximal number of interactions per sub-task  (this is the default value).
+  cell_split_size:  400      # Maximal number of particles per cell (this is the default value).
+
+# Parameters governing the time integration
+TimeIntegration:
+  time_begin: 0.    # The starting time of the simulation (in internal units).
+  time_end:   1.    # The end time of the simulation (in internal units).
+  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 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).
+  delta_neighbours:      1.       # The tolerance for the targetted number of neighbours.
+  CFL_condition:         0.1      # Courant-Friedrich-Levy condition for time integration.
+  max_ghost_iterations:  30       # Maximal number of iterations allowed to converge towards the smoothing length.
+  max_smoothing_length:  3.       # Maximal smoothing length allowed (in internal units).
+
+# Parameters related to the initial conditions
+InitialConditions:
+  file_name:  SedovBlast/sedov.hdf5 # The file to read
+  h_scaling:  1.                    # A scaling factor to apply to all smoothing lengths in the ICs.
+  shift_x:    0.                    # A shift to apply to all particles read from the ICs (in internal units).
+  shift_y:    0.
+  shift_z:    0.
+
+# Parameters govering domain decomposition
+DomainDecomposition:
+  initial_type:       m     # The initial strategy ("g", "m", "w", or "v"). See documentation for details.
+  initial_grid_x:    10     # Grid size if the 'g' strategy is chosen.
+  initial_grid_y:    10
+  initial_grid_z:    10
+  repartition_type:   b     # The re-decomposition strategy ("n", "b", "v", "e" or "x"). See documentation for details.
+ 
diff --git a/examples/runs.sh b/examples/runs.sh
deleted file mode 100755
index 339d8659675843f2491068ed8d30b528cb147c34..0000000000000000000000000000000000000000
--- a/examples/runs.sh
+++ /dev/null
@@ -1,53 +0,0 @@
-#!/bin/bash
-
-# Set some global stuff
-export OMP_WAIT_POLICY=PASSIVE
-
-# Generate the initial conditions if they are not present.
-if [ ! -e SodShock/sodShock.hdf5 ]
-then
-    echo "Generating initial conditions for the SodShock example..."
-    cd SodShock
-    python makeIC.py
-    cd ..
-fi
-if [ ! -e SedovBlast/sedov.hdf5 ]
-then
-    echo "Generating initial conditions for the SedovBlast example..."
-    cd SedovBlast/
-    python makeIC_fcc.py
-    cd ..
-fi
-if [ ! -e CosmoVolume/cosmoVolume.hdf5 ]
-then
-    echo "Downloading initial conditions for the CosmoVolume example..."
-    cd CosmoVolume
-    ./getIC.sh
-    cd ..
-fi
-
-
-# Loop over number of cores
-for cpu in {1..32}
-do
-
-    # Sod-Shock runs
-    if [ ! -e SodShock_${cpu}.dump ]
-    then
-        ./swift -t $cpu -f SodShock/sodShock.hdf5 -m 0.01 -w 5000 -c 1. -d 1e-7 -e 0.01 > SodShock_fixed_${cpu}.dump
-    fi
-    
-    # Sedov blast
-    if [ ! -e SedovBlast_${cpu}.dump ]
-    then
-        ./swift -t $cpu -f SedovBlast/sedov.hdf5 -m 0.02 -w 5000 -c 1. -d 1e-7 -e 0.01 > SedovBlast_fixed_${cpu}.dump
-    fi
-    
-    # Cosmological volume
-    if [ ! -e CosmoVolume_${cpu}.dump ]
-    then
-        ./swift -t $cpu -f CosmoVolume/cosmoVolume.hdf5 -m 0.6 -w 5000 -c 1. -d 1e-7 -e 0.01 > CosmoVolume_fixed_${cpu}.dump
-    fi
-
-done
-
diff --git a/src/common_io.c b/src/common_io.c
index bc981ecab31c56925ce08bfafb1a3a16aeee104b..2a635723d5bd4db7bce0a0172e8c083bf479ac32 100644
--- a/src/common_io.c
+++ b/src/common_io.c
@@ -282,15 +282,15 @@ void writeUnitSystem(hid_t h_file, struct UnitSystem* us) {
   if (h_grpunit < 0) error("Error while creating Unit System group");
 
   writeAttribute_d(h_grpunit, "Unit mass in cgs (U_M)",
-                   getBaseUnit(us, UNIT_MASS));
+                   units_get_base_unit(us, UNIT_MASS));
   writeAttribute_d(h_grpunit, "Unit length in cgs (U_L)",
-                   getBaseUnit(us, UNIT_LENGTH));
+                   units_get_base_unit(us, UNIT_LENGTH));
   writeAttribute_d(h_grpunit, "Unit time in cgs (U_t)",
-                   getBaseUnit(us, UNIT_TIME));
+                   units_get_base_unit(us, UNIT_TIME));
   writeAttribute_d(h_grpunit, "Unit current in cgs (U_I)",
-                   getBaseUnit(us, UNIT_CURRENT));
+                   units_get_base_unit(us, UNIT_CURRENT));
   writeAttribute_d(h_grpunit, "Unit temperature in cgs (U_T)",
-                   getBaseUnit(us, UNIT_TEMPERATURE));
+                   units_get_base_unit(us, UNIT_TEMPERATURE));
 
   H5Gclose(h_grpunit);
 }
diff --git a/src/const.h b/src/const.h
index 3bd9edff8227a87d040ec7309998364c946307af..6a52ec4796a4904629a57ffa8b32a3107bde263e 100644
--- a/src/const.h
+++ b/src/const.h
@@ -70,9 +70,4 @@
 #define GADGET2_SPH
 //#define DEFAULT_SPH
 
-/* System of units */
-#define const_unit_length_in_cgs 1   /* 3.08567810e16  /\* 1Mpc *\/ */
-#define const_unit_mass_in_cgs 1     /* 1.9891e33      /\* 1 M_sun *\/ */
-#define const_unit_velocity_in_cgs 1 /* 1e5            /\* km s^-1 *\/ */
-
 #endif /* SWIFT_CONST_H */
diff --git a/src/engine.c b/src/engine.c
index aa97dc42aee79e53f533ccbc7f061b3f87582638..e49d6da779d4333a00a60da920144d92a9241305 100644
--- a/src/engine.c
+++ b/src/engine.c
@@ -56,10 +56,11 @@
 #include "partition.h"
 #include "timers.h"
 
-const char *engine_policy_names[12] = {
-    "none",          "rand",   "steal",        "keep",
-    "block",         "fix_dt", "cpu_tight",    "mpi",
-    "numa_affinity", "hydro",  "self_gravity", "external_gravity"};
+const char *engine_policy_names[13] = {
+    "none",                 "rand",   "steal",        "keep",
+    "block",                "fix_dt", "cpu_tight",    "mpi",
+    "numa_affinity",        "hydro",  "self_gravity", "external_gravity",
+    "cosmology_integration"};
 
 /** The rank of the engine as a global variable (for messages). */
 int engine_rank;
@@ -1650,7 +1651,7 @@ void engine_rebuild(struct engine *e) {
     error("engine_marktasks failed after space_rebuild.");
 
   /* Print the status of the system */
-  engine_print_task_counts(e);
+  if (e->verbose) engine_print_task_counts(e);
 
   if (e->verbose)
     message("took %.3f %s.", clocks_from_ticks(getticks() - tic),
@@ -2026,8 +2027,8 @@ void engine_step(struct engine *e) {
   if (e->nodeID == 0) {
 
     /* Print some information to the screen */
-    printf("%d %e %e %d %d %.3f\n", e->step, e->time, e->timeStep, updates,
-           g_updates, e->wallclock_time);
+    printf("  %6d %14e %14e %10d %10d %21.3f\n", e->step, e->time, e->timeStep,
+           updates, g_updates, e->wallclock_time);
     fflush(stdout);
 
     /* Write some energy statistics */
@@ -2320,30 +2321,25 @@ static bool hyperthreads_present(void) {
  *
  * @param e The #engine.
  * @param s The #space in which this #runner will run.
- * @param dt The initial time step to use.
- * @param nr_threads The number of threads to spawn.
- * @param nr_queues The number of task queues to create.
+ * @param params The parsed parameter file.
  * @param nr_nodes The number of MPI ranks.
  * @param nodeID The MPI rank of this node.
  * @param policy The queuing policy to use.
- * @param timeBegin Time at the begininning of the simulation.
- * @param timeEnd Time at the end of the simulation.
- * @param dt_min Minimal allowed timestep (unsed with fixdt policy)
- * @param dt_max Maximal allowed timestep
  * @param verbose Is this #engine talkative ?
  */
 
-void engine_init(struct engine *e, struct space *s, float dt, int nr_threads,
-                 int nr_queues, int nr_nodes, int nodeID, int policy,
-                 float timeBegin, float timeEnd, float dt_min, float dt_max,
-                 int verbose) {
+void engine_init(struct engine *e, struct space *s,
+                 const struct swift_params *params, int nr_nodes, int nodeID,
+                 int policy, int verbose) {
+
+  /* Clean-up everything */
+  bzero(e, sizeof(struct engine));
 
   /* Store the values. */
   e->s = s;
-  e->nr_threads = nr_threads;
+  e->nr_threads = parser_get_param_int(params, "Scheduler:nr_threads");
   e->policy = policy;
   e->step = 0;
-  e->nullstep = 0;
   e->nr_nodes = nr_nodes;
   e->nodeID = nodeID;
   e->proxy_ind = NULL;
@@ -2352,15 +2348,15 @@ void engine_init(struct engine *e, struct space *s, float dt, int nr_threads,
   e->forcerepart = REPART_NONE;
   e->links = NULL;
   e->nr_links = 0;
-  e->timeBegin = timeBegin;
-  e->timeEnd = timeEnd;
-  e->timeOld = timeBegin;
-  e->time = timeBegin;
+  e->timeBegin = parser_get_param_double(params, "TimeIntegration:time_begin");
+  e->timeEnd = parser_get_param_double(params, "TimeIntegration:time_end");
+  e->timeOld = e->timeBegin;
+  e->time = e->timeBegin;
   e->ti_old = 0;
   e->ti_current = 0;
   e->timeStep = 0.;
-  e->dt_min = dt_min;
-  e->dt_max = dt_max;
+  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;
   e->verbose = verbose;
   e->count_step = 0;
@@ -2370,6 +2366,11 @@ void engine_init(struct engine *e, struct space *s, float dt, int nr_threads,
   /* Make the space link back to the engine. */
   s->e = e;
 
+  /* Get the number of queues */
+  int nr_queues = parser_get_param_int(params, "Scheduler:nr_queues");
+  if (nr_queues <= 0) nr_queues = e->nr_threads;
+  s->nr_queues = nr_queues;
+
 #if defined(HAVE_SETAFFINITY)
   const int nr_cores = sysconf(_SC_NPROCESSORS_ONLN);
   int cpuid[nr_cores];
@@ -2465,15 +2466,19 @@ void engine_init(struct engine *e, struct space *s, float dt, int nr_threads,
   engine_print_policy(e);
 
   /* Print information about the hydro scheme */
-  if (e->nodeID == 0) message("Hydrodynamic scheme: %s", SPH_IMPLEMENTATION);
-  if (e->nodeID == 0) message("Hydrodynamic kernel: %s", kernel_name);
+  if ((e->policy & engine_policy_hydro) == engine_policy_hydro) {
+    if (e->nodeID == 0) message("Hydrodynamic scheme: %s.", SPH_IMPLEMENTATION);
+    if (e->nodeID == 0)
+      message("Hydrodynamic kernel: %s with %.2f +/- %.2f neighbours.",
+              kernel_name, kernel_nwneigh, const_delta_nwneigh);
+  }
 
   /* Check we have sensible time bounds */
-  if (timeBegin >= timeEnd)
+  if (e->timeBegin >= e->timeEnd)
     error(
         "Final simulation time (t_end = %e) must be larger than the start time "
         "(t_beg = %e)",
-        timeEnd, timeBegin);
+        e->timeEnd, e->timeBegin);
 
   /* Check we have sensible time-step values */
   if (e->dt_min > e->dt_max)
@@ -2483,7 +2488,7 @@ void engine_init(struct engine *e, struct space *s, float dt, int nr_threads,
         e->dt_min, e->dt_max);
 
   /* Deal with timestep */
-  e->timeBase = (timeEnd - timeBegin) / max_nr_timesteps;
+  e->timeBase = (e->timeEnd - e->timeBegin) / max_nr_timesteps;
   e->ti_current = 0;
 
   /* Fixed time-step case */
@@ -2502,12 +2507,12 @@ void engine_init(struct engine *e, struct space *s, float dt, int nr_threads,
     if (e->nodeID == 0) {
       message("Absolute minimal timestep size: %e", e->timeBase);
 
-      float dt_min = timeEnd - timeBegin;
+      float dt_min = e->timeEnd - e->timeBegin;
       while (dt_min > e->dt_min) dt_min /= 2.f;
 
       message("Minimal timestep size (on time-line): %e", dt_min);
 
-      float dt_max = timeEnd - timeBegin;
+      float dt_max = e->timeEnd - e->timeBegin;
       while (dt_max > e->dt_max) dt_max /= 2.f;
 
       message("Maximal timestep size (on time-line): %e", dt_max);
@@ -2541,10 +2546,9 @@ void engine_init(struct engine *e, struct space *s, float dt, int nr_threads,
   e->barrier_launchcount = 0;
 
   /* Init the scheduler with enough tasks for the initial sorting tasks. */
-  int nr_tasks = 2 * s->tot_cells + e->nr_threads;
+  const int nr_tasks = 2 * s->tot_cells + 2 * e->nr_threads;
   scheduler_init(&e->sched, e->s, nr_tasks, nr_queues, scheduler_flag_steal,
                  e->nodeID);
-  s->nr_queues = nr_queues;
 
   /* Create the sorting tasks. */
   for (int i = 0; i < e->nr_threads; i++) {
@@ -2558,10 +2562,10 @@ void engine_init(struct engine *e, struct space *s, float dt, int nr_threads,
   scheduler_ranktasks(&e->sched);
 
   /* Allocate and init the threads. */
-  if ((e->runners =
-           (struct runner *)malloc(sizeof(struct runner) * nr_threads)) == NULL)
+  if ((e->runners = (struct runner *)malloc(sizeof(struct runner) *
+                                            e->nr_threads)) == NULL)
     error("Failed to allocate threads array.");
-  for (int k = 0; k < nr_threads; k++) {
+  for (int k = 0; k < e->nr_threads; k++) {
     e->runners[k].id = k;
     e->runners[k].e = e;
     e->barrier_running += 1;
@@ -2573,7 +2577,7 @@ void engine_init(struct engine *e, struct space *s, float dt, int nr_threads,
 
       /* Set a reasonable queue ID. */
       e->runners[k].cpuid = cpuid[k % nr_cores];
-      if (nr_queues < nr_threads)
+      if (nr_queues < e->nr_threads)
         e->runners[k].qid = cpuid[k % nr_cores] * nr_queues / nr_cores;
       else
         e->runners[k].qid = k;
@@ -2592,7 +2596,7 @@ void engine_init(struct engine *e, struct space *s, float dt, int nr_threads,
 #endif
     } else {
       e->runners[k].cpuid = k;
-      e->runners[k].qid = k * nr_queues / nr_threads;
+      e->runners[k].qid = k * nr_queues / e->nr_threads;
     }
     // message( "runner %i on cpuid=%i with qid=%i." , e->runners[k].id ,
     // e->runners[k].cpuid , e->runners[k].qid );
diff --git a/src/engine.h b/src/engine.h
index 4d1860b9eed0203bf9bf75711ec6e6549d837fe7..e1c3f61d1293fc01e24b9bcb0673d75fa3ce4648 100644
--- a/src/engine.h
+++ b/src/engine.h
@@ -38,6 +38,7 @@
 #include "scheduler.h"
 #include "space.h"
 #include "task.h"
+#include "parser.h"
 #include "partition.h"
 
 /* Some constants. */
@@ -53,7 +54,8 @@ enum engine_policy {
   engine_policy_setaffinity = (1 << 7),
   engine_policy_hydro = (1 << 8),
   engine_policy_self_gravity = (1 << 9),
-  engine_policy_external_gravity = (1 << 10)
+  engine_policy_external_gravity = (1 << 10),
+  engine_policy_cosmology = (1 << 11)
 };
 
 extern const char *engine_policy_names[];
@@ -126,7 +128,7 @@ struct engine {
   FILE *file_stats;
 
   /* The current step number. */
-  int step, nullstep;
+  int step;
 
   /* The number of particles updated in the previous step. */
   int count_step;
@@ -166,10 +168,9 @@ struct engine {
 
 /* Function prototypes. */
 void engine_barrier(struct engine *e, int tid);
-void engine_init(struct engine *e, struct space *s, float dt, int nr_threads,
-                 int nr_queues, int nr_nodes, int nodeID, int policy,
-                 float timeBegin, float timeEnd, float dt_min, float dt_max,
-                 int verbose);
+void engine_init(struct engine *e, struct space *s,
+                 const struct swift_params *params, int nr_nodes, int nodeID,
+                 int policy, int verbose);
 void engine_launch(struct engine *e, int nr_runners, unsigned int mask,
                    unsigned int submask);
 void engine_prepare(struct engine *e);
diff --git a/src/parallel_io.c b/src/parallel_io.c
index 0076c225e1c5361287280f8a567c8062aefd914e..d1c739b59021f38b2259f82dd06c547e0e7c147d 100644
--- a/src/parallel_io.c
+++ b/src/parallel_io.c
@@ -274,11 +274,11 @@ void writeArrayBackEnd(hid_t grp, char* fileName, FILE* xmfFile,
                  type);
 
   /* Write unit conversion factors for this data set */
-  conversionString(buffer, us, convFactor);
+  units_conversion_string(buffer, us, convFactor);
   writeAttribute_d(h_data, "CGS conversion factor",
-                   conversionFactor(us, convFactor));
-  writeAttribute_f(h_data, "h-scale exponent", hFactor(us, convFactor));
-  writeAttribute_f(h_data, "a-scale exponent", aFactor(us, convFactor));
+                   units_conversion_factor(us, convFactor));
+  writeAttribute_f(h_data, "h-scale exponent", units_h_factor(us, convFactor));
+  writeAttribute_f(h_data, "a-scale exponent", units_a_factor(us, convFactor));
   writeAttribute_s(h_data, "Conversion factor", buffer);
 
   /* Free and close everything */
@@ -349,6 +349,7 @@ void writeArrayBackEnd(hid_t grp, char* fileName, FILE* xmfFile,
  * @param parts (output) The array of #part read from the file.
  * @param N (output) The number of particles read from the file.
  * @param periodic (output) 1 if the volume is periodic, 0 if not.
+ * @param dry_run If 1, don't read the particle. Only allocates the arrays.
  *
  * Opens the HDF5 file fileName and reads the particles contained
  * in the parts array. N is the returned number of particles found
@@ -363,7 +364,7 @@ void writeArrayBackEnd(hid_t grp, char* fileName, FILE* xmfFile,
 void read_ic_parallel(char* fileName, double dim[3], struct part** parts,
                       struct gpart** gparts, size_t* Ngas, size_t* Ngparts,
                       int* periodic, int mpi_rank, int mpi_size, MPI_Comm comm,
-                      MPI_Info info) {
+                      MPI_Info info, int dry_run) {
   hid_t h_file = 0, h_grp = 0;
   /* GADGET has only cubic boxes (in cosmological mode) */
   double boxSize[3] = {0.0, -1.0, -1.0};
@@ -469,13 +470,15 @@ void read_ic_parallel(char* fileName, double dim[3], struct part** parts,
     switch (ptype) {
 
       case GAS:
-        hydro_read_particles(h_grp, N[ptype], N_total[ptype], offset[ptype],
-                             *parts);
+        if (!dry_run)
+          hydro_read_particles(h_grp, N[ptype], N_total[ptype], offset[ptype],
+                               *parts);
         break;
 
       case DM:
-        darkmatter_read_particles(h_grp, N[ptype], N_total[ptype],
-                                  offset[ptype], *gparts);
+        if (!dry_run)
+          darkmatter_read_particles(h_grp, N[ptype], N_total[ptype],
+                                    offset[ptype], *gparts);
         break;
 
       default:
@@ -487,10 +490,10 @@ void read_ic_parallel(char* fileName, double dim[3], struct part** parts,
   }
 
   /* Prepare the DM particles */
-  prepare_dm_gparts(*gparts, Ndm);
+  if (!dry_run) prepare_dm_gparts(*gparts, Ndm);
 
   /* Now duplicate the hydro particle into gparts */
-  duplicate_hydro_gparts(*parts, *gparts, *Ngas, Ndm);
+  if (!dry_run) duplicate_hydro_gparts(*parts, *gparts, *Ngas, Ndm);
 
   /* message("Done Reading particles..."); */
 
diff --git a/src/parallel_io.h b/src/parallel_io.h
index 663f0aabac44c08682b964512839b925673ea5c5..f3691cb29b8d5e7f17382f1f81ba230c3898a929 100644
--- a/src/parallel_io.h
+++ b/src/parallel_io.h
@@ -34,7 +34,7 @@
 void read_ic_parallel(char* fileName, double dim[3], struct part** parts,
                       struct gpart** gparts, size_t* Ngas, size_t* Ngparts,
                       int* periodic, int mpi_rank, int mpi_size, MPI_Comm comm,
-                      MPI_Info info);
+                      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,
diff --git a/src/parser.c b/src/parser.c
index 06dc819842d54d952704e4e0c40ebec5b561f691..0f767bc434ef596df403fb12d3ae0f77ea546df3 100644
--- a/src/parser.c
+++ b/src/parser.c
@@ -24,6 +24,7 @@
 /* Needs to be included so that strtok returns char * instead of a int *. */
 #include <string.h>
 #include <stdlib.h>
+#include <ctype.h>
 
 /* This object's header. */
 #include "parser.h"
@@ -31,9 +32,23 @@
 /* Local headers. */
 #include "error.h"
 
+#define PARSER_COMMENT_STRING "#"
+#define PARSER_COMMENT_CHAR '#'
+#define PARSER_VALUE_CHAR ':'
+#define PARSER_VALUE_STRING ":"
+#define PARSER_START_OF_FILE "---"
+#define PARSER_END_OF_FILE "..."
+
 /* Private functions. */
-static int count_char(char *str, char val);
-static void parse_line(FILE *fp, struct swift_params *params);
+static int count_char(const char *str, char val);
+static int is_empty(const char *str);
+static int count_indentation(const char *str);
+static void parse_line(char *line, struct swift_params *params);
+static void parse_value(char *line, struct swift_params *params);
+static void parse_section_param(char *line, int *isFirstParam,
+                                char *sectionName, struct swift_params *params);
+
+static int lineNumber = 0;
 
 /**
  * @brief Reads an input file and stores each parameter in a structure.
@@ -43,24 +58,29 @@ static void parse_line(FILE *fp, struct swift_params *params);
  */
 
 void parser_read_file(const char *file_name, struct swift_params *params) {
+  /* Open file for reading */
+  FILE *file = fopen(file_name, "r");
 
-  FILE *fp;
+  /* Line to parsed. */
+  char line[PARSER_MAX_LINE_SIZE];
 
+  /* Initialise parameter count. */
   params->count = 0;
 
-  /* Open file for reading */
-  fp = fopen(file_name, "r");
-
-  if (fp == NULL) {
+  /* Check if parameter file exits. */
+  if (file == NULL) {
     error("Error opening parameter file: %s", file_name);
   }
 
   /* Read until the end of the file is reached.*/
-  while (!feof(fp)) {
-    parse_line(fp, params);
+  while (!feof(file)) {
+    if (fgets(line, PARSER_MAX_LINE_SIZE, file) != NULL) {
+      lineNumber++;
+      parse_line(line, params);
+    }
   }
 
-  fclose(fp);
+  fclose(file);
 }
 
 /**
@@ -68,10 +88,11 @@ void parser_read_file(const char *file_name, struct swift_params *params) {
  *
  * @param str String to be checked
  * @param val Character to be counted
+ *
+ * @return Number of occurrences of val inside str
  */
 
-static int count_char(char *str, char val) {
-
+static int count_char(const char *str, char val) {
   int count = 0;
 
   /* Check if the line contains the character */
@@ -83,110 +104,277 @@ static int count_char(char *str, char val) {
 }
 
 /**
- * @brief Parses a line from a file and stores any parameters in a structure.
+ * @brief Counts the number of white spaces that prefix a string.
  *
- * @param fp File pointer to file to be read
- * @param params Structure to be populated from file
+ * @param str String to be checked
  *
+ * @return Number of white spaces prefixing str
  */
 
-static void parse_line(FILE *fp, struct swift_params *params) {
+static int count_indentation(const char *str) {
+  int count = 0;
 
-  char line[PARSER_MAX_LINE_SIZE];
-  char trim_line[PARSER_MAX_LINE_SIZE];
+  /* Check if the line contains the character */
+  while (*(++str) == ' ') {
+    count++;
+  }
+  return count;
+}
 
-  /* Read a line of the file */
-  if (fgets(line, PARSER_MAX_LINE_SIZE, fp) != NULL) {
+/**
+ * @brief Checks if a string is empty.
+ *
+ * @param str String to be checked
+ *
+ * @return Returns 1 if str is empty, 0 otherwise
+ */
 
+static int is_empty(const char *str) {
+  int retParam = 1;
+  while (*str != '\0') {
+    if (!isspace(*str)) {
+      retParam = 0;
+      break;
+    }
+    str++;
+  }
+
+  return retParam;
+}
+
+/**
+ * @brief Parses a line from a file and stores any parameters in a structure.
+ *
+ * @param line Line to be parsed.
+ * @param params Structure to be populated from file.
+ */
+static void parse_line(char *line, struct swift_params *params) {
+  /* Parse line if it doesn't begin with a comment. */
+  if (*line != PARSER_COMMENT_CHAR) {
+    char trim_line[PARSER_MAX_LINE_SIZE];
+    char tmp_str[PARSER_MAX_LINE_SIZE];
     char *token;
-    /* Remove comments */
-    token = strtok(line, PARSER_COMMENT_CHAR);
-    strcpy(trim_line, token);
-
-    /* Check if the line contains a value */
-    if (strchr(trim_line, PARSER_VALUE_CHAR)) {
-      /* Check for more than one parameter on the same line. */
-      if (count_char(trim_line, PARSER_VALUE_CHAR) > 1) {
-        error("Found more than one parameter in '%s', only one allowed.", line);
-      } else {
-        /* Take first token as the parameter name. */
-        token = strtok(trim_line, PARSER_VALUE_STRING);
-        strcpy(params->data[params->count].name, token);
-
-        /* Take second token as the parameter value. */
-        token = strtok(NULL, " #\n");
-        strcpy(params->data[params->count++].value, token);
+
+    /* Remove comments at the end of a line. */
+    token = strtok(line, PARSER_COMMENT_STRING);
+    strcpy(tmp_str, token);
+
+    /* Check if the line is just white space. */
+    if (!is_empty(tmp_str)) {
+      /* Trim '\n' characters from string. */
+      token = strtok(tmp_str, "\n");
+      strcpy(trim_line, token);
+
+      /* Check if the line contains a value and parse it. */
+      if (strchr(trim_line, PARSER_VALUE_CHAR)) {
+        parse_value(trim_line, params);
+      }
+      /* Check for invalid lines,not including the start and end of file. */
+      /* Note: strcmp returns 0 if both strings are the same.*/
+      else if (strcmp(trim_line, PARSER_START_OF_FILE) &&
+               strcmp(trim_line, PARSER_END_OF_FILE)) {
+        error("Invalid line:%d '%s'.", lineNumber, trim_line);
       }
     }
   }
 }
 
+/**
+ * @brief Performs error checking and stores a parameter in a structure.
+ *
+ * @param line Line containing the parameter
+ * @param params Structure to be written to
+ *
+ */
+
+static void parse_value(char *line, struct swift_params *params) {
+  static int inSection = 0;
+  static char section[PARSER_MAX_LINE_SIZE]; /* Keeps track of current section
+                                                name. */
+  static int isFirstParam = 1;
+  char tmpStr[PARSER_MAX_LINE_SIZE];
+
+  char *token;
+
+  /* Check for more than one value on the same line. */
+  if (count_char(line, PARSER_VALUE_CHAR) > 1) {
+    error("Inavlid line:%d '%s', only one value allowed per line.", lineNumber,
+          line);
+  }
+
+  /* Check that standalone parameters have correct indentation. */
+  if (!inSection && *line == ' ') {
+    error(
+        "Invalid line:%d '%s', standalone parameter defined with incorrect "
+        "indentation.",
+        lineNumber, line);
+  }
+
+  /* Check that it is a parameter inside a section.*/
+  if (*line == ' ' || *line == '\t') {
+    parse_section_param(line, &isFirstParam, section, params);
+  } else {/*Else it is the start of a new section or standalone parameter. */
+    /* Take first token as the parameter name. */
+    token = strtok(line, " :\t");
+    strcpy(tmpStr, token);
+
+    /* Take second token as the parameter value. */
+    token = strtok(NULL, " #\n");
+
+    /* If second token is NULL then the line must be a section heading. */
+    if (token == NULL) {
+      strcat(tmpStr, PARSER_VALUE_STRING);
+      strcpy(section, tmpStr);
+      inSection = 1;
+      isFirstParam = 1;
+    } else {
+      /* Must be a standalone parameter so no need to prefix name with a
+       * section. */
+      strcpy(params->data[params->count].name, tmpStr);
+      strcpy(params->data[params->count++].value, token);
+      inSection = 0;
+      isFirstParam = 1;
+    }
+  }
+}
+
+/**
+ * @brief Parses a parameter that appears in a section and stores it in a
+ *structure.
+ *
+ * @param line Line containing the parameter
+ * @param isFirstParam Shows if the first parameter of a section has been found
+ * @param sectionName String containing the current section name
+ * @param params Structure to be written to
+ *
+ */
+
+static void parse_section_param(char *line, int *isFirstParam,
+                                char *sectionName,
+                                struct swift_params *params) {
+  static int sectionIndent = 0;
+  char tmpStr[PARSER_MAX_LINE_SIZE];
+  char paramName[PARSER_MAX_LINE_SIZE];
+  char *token;
+
+  /* Count indentation of each parameter and check that it
+   * is consistent with the first parameter in the section. */
+  if (*isFirstParam) {
+    sectionIndent = count_indentation(line);
+    *isFirstParam = 0;
+  } else if (count_indentation(line) != sectionIndent) {
+    error("Invalid line:%d '%s', parameter has incorrect indentation.",
+          lineNumber, line);
+  }
+
+  /* Take first token as the parameter name and trim leading white space. */
+  token = strtok(line, " :\t");
+  strcpy(tmpStr, token);
+
+  /* Take second token as the parameter value. */
+  token = strtok(NULL, " #\n");
+
+  /* Prefix the parameter name with its section name and
+   * copy it into the parameter structure. */
+  strcpy(paramName, sectionName);
+  strcat(paramName, tmpStr);
+  strcpy(params->data[params->count].name, paramName);
+  strcpy(params->data[params->count++].value, token);
+}
+
 /**
  * @brief Retrieve integer parameter from structure.
  *
  * @param params Structure that holds the parameters
  * @param name Name of the parameter to be found
- * @param retParam Value of the parameter found
- *
+ * @return Value of the parameter found
  */
+int parser_get_param_int(const struct swift_params *params, const char *name) {
 
-void parser_get_param_int(struct swift_params *params, char *name,
-                          int *retParam) {
-
-  char str[128];
+  char str[PARSER_MAX_LINE_SIZE];
+  int retParam = 0;
 
   for (int i = 0; i < params->count; i++) {
-
     /*strcmp returns 0 if both strings are the same.*/
     if (!strcmp(name, params->data[i].name)) {
-
       /* Check that exactly one number is parsed. */
-      if (sscanf(params->data[i].value, "%d%s", retParam, str) != 1) {
+      if (sscanf(params->data[i].value, "%d%s", &retParam, str) != 1) {
         error(
             "Tried parsing int '%s' but found '%s' with illegal integer "
             "characters '%s'.",
             params->data[i].name, params->data[i].value, str);
       }
 
-      return;
+      return retParam;
     }
   }
 
-  message("Cannot find '%s' in the structure.", name);
+  error("Cannot find '%s' in the structure.", name);
+  return 0;
 }
 
 /**
- * @brief Retrieve float parameter from structure.
+ * @brief Retrieve char parameter from structure.
  *
  * @param params Structure that holds the parameters
  * @param name Name of the parameter to be found
- * @param retParam Value of the parameter found
- *
+ * @return Value of the parameter found
  */
+char parser_get_param_char(const struct swift_params *params,
+                           const char *name) {
 
-void parser_get_param_float(struct swift_params *params, char *name,
-                            float *retParam) {
-
-  char str[128];
+  char str[PARSER_MAX_LINE_SIZE];
+  char retParam = 0;
 
   for (int i = 0; i < params->count; i++) {
-
     /*strcmp returns 0 if both strings are the same.*/
     if (!strcmp(name, params->data[i].name)) {
+      /* Check that exactly one number is parsed. */
+      if (sscanf(params->data[i].value, "%c%s", &retParam, str) != 1) {
+        error(
+            "Tried parsing char '%s' but found '%s' with illegal char "
+            "characters '%s'.",
+            params->data[i].name, params->data[i].value, str);
+      }
 
+      return retParam;
+    }
+  }
+
+  error("Cannot find '%s' in the structure.", name);
+  return 0;
+}
+
+/**
+ * @brief Retrieve float parameter from structure.
+ *
+ * @param params Structure that holds the parameters
+ * @param name Name of the parameter to be found
+ * @return Value of the parameter found
+ */
+float parser_get_param_float(const struct swift_params *params,
+                             const char *name) {
+
+  char str[PARSER_MAX_LINE_SIZE];
+  float retParam = 0.f;
+
+  for (int i = 0; i < params->count; i++) {
+    /*strcmp returns 0 if both strings are the same.*/
+    if (!strcmp(name, params->data[i].name)) {
       /* Check that exactly one number is parsed. */
-      if (sscanf(params->data[i].value, "%f%s", retParam, str) != 1) {
+      if (sscanf(params->data[i].value, "%f%s", &retParam, str) != 1) {
         error(
             "Tried parsing float '%s' but found '%s' with illegal float "
             "characters '%s'.",
             params->data[i].name, params->data[i].value, str);
       }
 
-      return;
+      return retParam;
     }
   }
 
-  message("Cannot find '%s' in the structure.", name);
+  error("Cannot find '%s' in the structure.", name);
+  return 0.f;
 }
 
 /**
@@ -194,33 +382,30 @@ void parser_get_param_float(struct swift_params *params, char *name,
  *
  * @param params Structure that holds the parameters
  * @param name Name of the parameter to be found
- * @param retParam Value of the parameter found
- *
+ * @return Value of the parameter found
  */
+double parser_get_param_double(const struct swift_params *params,
+                               const char *name) {
 
-void parser_get_param_double(struct swift_params *params, char *name,
-                             double *retParam) {
-
-  char str[128];
+  char str[PARSER_MAX_LINE_SIZE];
+  double retParam = 0.;
 
   for (int i = 0; i < params->count; i++) {
-
     /*strcmp returns 0 if both strings are the same.*/
     if (!strcmp(name, params->data[i].name)) {
-
       /* Check that exactly one number is parsed. */
-      if (sscanf(params->data[i].value, "%lf", retParam) != 1) {
+      if (sscanf(params->data[i].value, "%lf%s", &retParam, str) != 1) {
         error(
             "Tried parsing double '%s' but found '%s' with illegal double "
             "characters '%s'.",
             params->data[i].name, params->data[i].value, str);
       }
-
-      return;
+      return retParam;
     }
   }
 
-  message("Cannot find '%s' in the structure.", name);
+  error("Cannot find '%s' in the structure.", name);
+  return 0.;
 }
 
 /**
@@ -228,32 +413,27 @@ void parser_get_param_double(struct swift_params *params, char *name,
  *
  * @param params Structure that holds the parameters
  * @param name Name of the parameter to be found
- * @param retParam Value of the parameter found
- *
+ * @param retParam (return) Value of the parameter found
  */
-
-void parser_get_param_string(struct swift_params *params, char *name,
-                             char *retParam) {
-
+void parser_get_param_string(const struct swift_params *params,
+                             const char *name, char *retParam) {
   for (int i = 0; i < params->count; i++) {
-
     /*strcmp returns 0 if both strings are the same.*/
     if (!strcmp(name, params->data[i].name)) {
       strcpy(retParam, params->data[i].value);
       return;
     }
   }
+
+  error("Cannot find '%s' in the structure.", name);
 }
 
 /**
  * @brief Prints the contents of the parameter structure.
  *
  * @param params Structure that holds the parameters
- *
  */
-
-void parser_print_params(struct swift_params *params) {
-
+void parser_print_params(const struct swift_params *params) {
   printf("\n--------------------------\n");
   printf("|  SWIFT Parameter File  |\n");
   printf("--------------------------\n");
@@ -263,3 +443,51 @@ void parser_print_params(struct swift_params *params) {
     printf("Parameter value: %s\n", params->data[i].value);
   }
 }
+
+/**
+ * @brief Write the contents of the parameter structure to a file in YAML
+ *format.
+ *
+ * @param params Structure that holds the parameters
+ * @param file_name Name of file to be written
+ */
+void parser_write_params_to_file(const struct swift_params *params,
+                                 const char *file_name) {
+  FILE *file = fopen(file_name, "w");
+  char section[PARSER_MAX_LINE_SIZE];
+  char param_name[PARSER_MAX_LINE_SIZE];
+  char *token;
+
+  /* Start of file identifier in YAML. */
+  fprintf(file, "%s\n", PARSER_START_OF_FILE);
+
+  for (int i = 0; i < params->count; i++) {
+    /* Check that the parameter name contains a section name. */
+    if (strchr(params->data[i].name, PARSER_VALUE_CHAR)) {
+      /* Copy the parameter name into a temporary string and find the section
+       * name. */
+      strcpy(param_name, params->data[i].name);
+      token = strtok(param_name, PARSER_VALUE_STRING);
+
+      /* If a new section name is found print it to the file. */
+      if (strcmp(token, section)) {
+        strcpy(section, token);
+        fprintf(file, "\n%s%c\n", section, PARSER_VALUE_CHAR);
+      }
+
+      /* Remove white space from parameter name and write it to the file. */
+      token = strtok(NULL, " #\n");
+
+      fprintf(file, "\t%s%c %s\n", token, PARSER_VALUE_CHAR,
+              params->data[i].value);
+    } else {
+      fprintf(file, "\n%s%c %s\n", params->data[i].name, PARSER_VALUE_CHAR,
+              params->data[i].value);
+    }
+  }
+
+  /* End of file identifier in YAML. */
+  fprintf(file, PARSER_END_OF_FILE);
+
+  fclose(file);
+}
diff --git a/src/parser.h b/src/parser.h
index 2fb4148944cd423da016341744cb6d58e222182e..7b2088ae12cdd5136a96baeabd01dd80255c8a3b 100644
--- a/src/parser.h
+++ b/src/parser.h
@@ -19,21 +19,20 @@
 #ifndef SWIFT_PARSER_H
 #define SWIFT_PARSER_H
 
-#include <stdio.h>
+/* Config parameters. */
+#include "../config.h"
 
-#define PARSER_MAX_LINE_SIZE 128
+/* Some constants. */
+#define PARSER_MAX_LINE_SIZE 256
 #define PARSER_MAX_NO_OF_PARAMS 512
 
-#define PARSER_COMMENT_CHAR "#"
-#define PARSER_VALUE_CHAR ':'
-#define PARSER_VALUE_STRING ":"
-#define PARSER_END_OF_FILE "..."
-
+/* A parameter in the input file */
 struct parameter {
   char name[PARSER_MAX_LINE_SIZE];
   char value[PARSER_MAX_LINE_SIZE];
 };
 
+/* The array of parameters read from a file */
 struct swift_params {
   struct parameter data[PARSER_MAX_NO_OF_PARAMS];
   int count;
@@ -41,14 +40,16 @@ struct swift_params {
 
 /* Public API. */
 void parser_read_file(const char *file_name, struct swift_params *params);
-void parser_print_params(struct swift_params *params);
-void parser_get_param_int(struct swift_params *params, char *name,
-                          int *retParam);
-void parser_get_param_float(struct swift_params *params, char *name,
-                            float *retParam);
-void parser_get_param_double(struct swift_params *params, char *name,
-                             double *retParam);
-void parser_get_param_string(struct swift_params *params, char *name,
-                             char *retParam);
+void parser_print_params(const struct swift_params *params);
+void parser_write_params_to_file(const struct swift_params *params,
+                                 const char *file_name);
+char parser_get_param_char(const struct swift_params *params, const char *name);
+int parser_get_param_int(const struct swift_params *params, const char *name);
+float parser_get_param_float(const struct swift_params *params,
+                             const char *name);
+double parser_get_param_double(const struct swift_params *params,
+                               const char *name);
+void parser_get_param_string(const struct swift_params *params,
+                             const char *name, char *retParam);
 
 #endif /* SWIFT_PARSER_H */
diff --git a/src/partition.c b/src/partition.c
index e9ecae60ad3945c918e4783602c521de5d1ae12d..e559273fac96c9ecb860126d2a0a05601861c343 100644
--- a/src/partition.c
+++ b/src/partition.c
@@ -52,6 +52,7 @@
 #include "error.h"
 #include "partition.h"
 #include "space.h"
+#include "tools.h"
 
 /* Maximum weight used for METIS. */
 #define metis_maxweight 10000.0f
@@ -908,6 +909,111 @@ void partition_initial_partition(struct partition *initial_partition,
   }
 }
 
+/**
+ * @brief Initialises the partition and re-partition scheme from the parameter
+ *file
+ *
+ * @param partition The #partition scheme to initialise.
+ * @param reparttype The repartition scheme to initialise.
+ * @param params The parsed parameter file.
+ * @param nr_nodes The number of MPI nodes we are running on.
+ */
+void partition_init(struct partition *partition,
+                    enum repartition_type *reparttype,
+                    const struct swift_params *params, int nr_nodes) {
+
+#ifdef WITH_MPI
+
+/* Defaults make use of METIS if available */
+#ifdef HAVE_METIS
+  *reparttype = REPART_METIS_BOTH;
+  partition->type = INITPART_METIS_NOWEIGHT;
+#else
+  *reparttype = REPART_NONE;
+  partition->type = INITPART_GRID;
+#endif
+
+  /* Set a default grid so that grid[0]*grid[1]*grid[2] == nr_nodes. */
+  factor(nr_nodes, &partition->grid[0], &partition->grid[1]);
+  factor(nr_nodes / partition->grid[1], &partition->grid[0],
+         &partition->grid[2]);
+  factor(partition->grid[0] * partition->grid[1], &partition->grid[1],
+         &partition->grid[0]);
+
+  /* Now let's check what the user wants as an initial domain*/
+  const char part_type =
+      parser_get_param_char(params, "DomainDecomposition:initial_type");
+
+  switch (part_type) {
+    case 'g':
+      partition->type = INITPART_GRID;
+      break;
+    case 'v':
+      partition->type = INITPART_VECTORIZE;
+      break;
+#ifdef HAVE_METIS
+    case 'm':
+      partition->type = INITPART_METIS_NOWEIGHT;
+      break;
+    case 'w':
+      partition->type = INITPART_METIS_WEIGHT;
+      break;
+    default:
+      message("Invalid choice of initial partition type '%c'.", part_type);
+      error("Permitted values are: 'g','m','v' or 'w'.");
+#else
+    default:
+      message("Invalid choice of initial partition type '%c'.", part_type);
+      error("Permitted values are: 'g' or 'v' when compiled without metis.");
+#endif
+  }
+
+  /* In case of grid, read more parameters */
+  if (part_type == 'g') {
+    partition->grid[0] =
+        parser_get_param_int(params, "DomainDecomposition:initial_grid_x");
+    partition->grid[1] =
+        parser_get_param_int(params, "DomainDecomposition:initial_grid_y");
+    partition->grid[2] =
+        parser_get_param_int(params, "DomainDecomposition:initial_grid_z");
+  }
+
+  /* Now let's check what the user wants as a repartition strategy */
+  const char repart_type =
+      parser_get_param_char(params, "DomainDecomposition:repartition_type");
+
+  switch (repart_type) {
+    case 'n':
+      *reparttype = REPART_NONE;
+      break;
+#ifdef HAVE_METIS
+    case 'b':
+      *reparttype = REPART_METIS_BOTH;
+      break;
+    case 'e':
+      *reparttype = REPART_METIS_EDGE;
+      break;
+    case 'v':
+      *reparttype = REPART_METIS_VERTEX;
+      break;
+    case 'x':
+      *reparttype = REPART_METIS_VERTEX_EDGE;
+      break;
+    default:
+      message("Invalid choice of re-partition type '%c'.", repart_type);
+      error("Permitted values are: 'b','e','n', 'v' or 'x'.");
+#else
+    default:
+      message("Invalid choice of re-partition type '%c'.", repart_type);
+      error("Permitted values are: 'n' when compiled without metis.");
+#endif
+  }
+
+#else
+  error("SWIFT was not compiled with MPI support");
+#endif
+}
+
 /*  General support */
 /*  =============== */
 
diff --git a/src/partition.h b/src/partition.h
index 3ab5f5a817bf5b77d45c7fb3313158b83d98e251..af3d0be9242d75ce2c319de034428dc934e3a925 100644
--- a/src/partition.h
+++ b/src/partition.h
@@ -19,6 +19,7 @@
 #ifndef SWIFT_PARTITION_H
 #define SWIFT_PARTITION_H
 
+#include "parser.h"
 #include "space.h"
 #include "task.h"
 
@@ -56,4 +57,8 @@ void partition_repartition(enum repartition_type reparttype, int nodeID,
 void partition_initial_partition(struct partition *initial_partition,
                                  int nodeID, int nr_nodes, struct space *s);
 
+void partition_init(struct partition *partition,
+                    enum repartition_type *reparttypestruct,
+                    const struct swift_params *params, int nr_nodes);
+
 #endif /* SWIFT_PARTITION_H */
diff --git a/src/serial_io.c b/src/serial_io.c
index 40bd2b1c8921f4acbfa0950984d6915ebd3d241e..10eab97f1bf118a842e274b521056d0d81b32db1 100644
--- a/src/serial_io.c
+++ b/src/serial_io.c
@@ -237,11 +237,11 @@ void prepareArray(hid_t grp, char* fileName, FILE* xmfFile,
   writeXMFline(xmfFile, fileName, partTypeGroupName, name, N_total, dim, type);
 
   /* Write unit conversion factors for this data set */
-  conversionString(buffer, us, convFactor);
+  units_conversion_string(buffer, us, convFactor);
   writeAttribute_d(h_data, "CGS conversion factor",
-                   conversionFactor(us, convFactor));
-  writeAttribute_f(h_data, "h-scale exponent", hFactor(us, convFactor));
-  writeAttribute_f(h_data, "a-scale exponent", aFactor(us, convFactor));
+                   units_conversion_factor(us, convFactor));
+  writeAttribute_f(h_data, "h-scale exponent", units_h_factor(us, convFactor));
+  writeAttribute_f(h_data, "a-scale exponent", units_a_factor(us, convFactor));
   writeAttribute_s(h_data, "Conversion factor", buffer);
 
   H5Pclose(h_prop);
@@ -412,6 +412,7 @@ void writeArrayBackEnd(hid_t grp, char* fileName, FILE* xmfFile,
  * @param mpi_size The number of MPI ranks
  * @param comm The MPI communicator
  * @param info The MPI information object
+ * @param dry_run If 1, don't read the particle. Only allocates the arrays.
  *
  * Opens the HDF5 file fileName and reads the particles contained
  * in the parts array. N is the returned number of particles found
@@ -420,13 +421,11 @@ void writeArrayBackEnd(hid_t grp, char* fileName, FILE* xmfFile,
  * @warning Can not read snapshot distributed over more than 1 file !!!
  * @todo Read snapshots distributed in more than one file.
  *
- * Calls #error() if an error occurs.
- *
  */
 void read_ic_serial(char* fileName, double dim[3], struct part** parts,
                     struct gpart** gparts, size_t* Ngas, size_t* Ngparts,
                     int* periodic, int mpi_rank, int mpi_size, MPI_Comm comm,
-                    MPI_Info info) {
+                    MPI_Info info, int dry_run) {
   hid_t h_file = 0, h_grp = 0;
   /* GADGET has only cubic boxes (in cosmological mode) */
   double boxSize[3] = {0.0, -1.0, -1.0};
@@ -516,10 +515,12 @@ void read_ic_serial(char* fileName, double dim[3], struct part** parts,
 
   /* message("Allocated %8.2f MB for particles.", *N * sizeof(struct part) / */
   /* 	  (1024.*1024.)); */
-
   /* message("BoxSize = %lf", dim[0]); */
   /* message("NumPart = [%zd, %zd] Total = %zd", *Ngas, Ndm, *Ngparts); */
 
+  /* For dry runs, only need to do this on rank 0 */
+  if (dry_run) mpi_size = 1;
+
   /* Now loop over ranks and read the data */
   for (int rank = 0; rank < mpi_size; ++rank) {
 
@@ -549,13 +550,15 @@ void read_ic_serial(char* fileName, double dim[3], struct part** parts,
         switch (ptype) {
 
           case GAS:
-            hydro_read_particles(h_grp, N[ptype], N_total[ptype], offset[ptype],
-                                 *parts);
+            if (!dry_run)
+              hydro_read_particles(h_grp, N[ptype], N_total[ptype],
+                                   offset[ptype], *parts);
             break;
 
           case DM:
-            darkmatter_read_particles(h_grp, N[ptype], N_total[ptype],
-                                      offset[ptype], *gparts);
+            if (!dry_run)
+              darkmatter_read_particles(h_grp, N[ptype], N_total[ptype],
+                                        offset[ptype], *gparts);
             break;
 
           default:
@@ -575,10 +578,10 @@ void read_ic_serial(char* fileName, double dim[3], struct part** parts,
   }
 
   /* Prepare the DM particles */
-  prepare_dm_gparts(*gparts, Ndm);
+  if (!dry_run) prepare_dm_gparts(*gparts, Ndm);
 
   /* Now duplicate the hydro particle into gparts */
-  duplicate_hydro_gparts(*parts, *gparts, *Ngas, Ndm);
+  if (!dry_run) duplicate_hydro_gparts(*parts, *gparts, *Ngas, Ndm);
 
   /* message("Done Reading particles..."); */
 }
diff --git a/src/serial_io.h b/src/serial_io.h
index 5a34d420cfabd88d4147e3f3630e0efe89951c41..74ab8326dbeeb955e354687059cdd595657285f0 100644
--- a/src/serial_io.h
+++ b/src/serial_io.h
@@ -34,7 +34,7 @@
 void read_ic_serial(char* fileName, double dim[3], struct part** parts,
                     struct gpart** gparts, size_t* Ngas, size_t* Ngparts,
                     int* periodic, int mpi_rank, int mpi_size, MPI_Comm comm,
-                    MPI_Info info);
+                    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);
diff --git a/src/single_io.c b/src/single_io.c
index 801428433ef5170082b68dec425e52f845bb41ae..1dc71087e102ff884dba7b7d4b6dcd6339335cac 100644
--- a/src/single_io.c
+++ b/src/single_io.c
@@ -241,11 +241,11 @@ void writeArrayBackEnd(hid_t grp, char* fileName, FILE* xmfFile,
   writeXMFline(xmfFile, fileName, partTypeGroupName, name, N, dim, type);
 
   /* Write unit conversion factors for this data set */
-  conversionString(buffer, us, convFactor);
+  units_conversion_string(buffer, us, convFactor);
   writeAttribute_d(h_data, "CGS conversion factor",
-                   conversionFactor(us, convFactor));
-  writeAttribute_f(h_data, "h-scale exponent", hFactor(us, convFactor));
-  writeAttribute_f(h_data, "a-scale exponent", aFactor(us, convFactor));
+                   units_conversion_factor(us, convFactor));
+  writeAttribute_f(h_data, "h-scale exponent", units_h_factor(us, convFactor));
+  writeAttribute_f(h_data, "a-scale exponent", units_a_factor(us, convFactor));
   writeAttribute_s(h_data, "Conversion factor", buffer);
 
   /* Free and close everything */
@@ -319,6 +319,7 @@ void writeArrayBackEnd(hid_t grp, char* fileName, FILE* xmfFile,
  * @param Ngas (output) number of Gas particles read.
  * @param Ngparts (output) The number of #gpart read.
  * @param periodic (output) 1 if the volume is periodic, 0 if not.
+ * @param dry_run If 1, don't read the particle. Only allocates the arrays.
  *
  * Opens the HDF5 file fileName and reads the particles contained
  * in the parts array. N is the returned number of particles found
@@ -327,12 +328,10 @@ void writeArrayBackEnd(hid_t grp, char* fileName, FILE* xmfFile,
  * @warning Can not read snapshot distributed over more than 1 file !!!
  * @todo Read snapshots distributed in more than one file.
  *
- * Calls #error() if an error occurs.
- *
  */
 void read_ic_single(char* fileName, double dim[3], struct part** parts,
                     struct gpart** gparts, size_t* Ngas, size_t* Ngparts,
-                    int* periodic) {
+                    int* periodic, int dry_run) {
   hid_t h_file = 0, h_grp = 0;
   /* GADGET has only cubic boxes (in cosmological mode) */
   double boxSize[3] = {0.0, -1.0, -1.0};
@@ -426,11 +425,11 @@ void read_ic_single(char* fileName, double dim[3], struct part** parts,
     switch (ptype) {
 
       case GAS:
-        hydro_read_particles(h_grp, *Ngas, *Ngas, 0, *parts);
+        if (!dry_run) hydro_read_particles(h_grp, *Ngas, *Ngas, 0, *parts);
         break;
 
       case DM:
-        darkmatter_read_particles(h_grp, Ndm, Ndm, 0, *gparts);
+        if (!dry_run) darkmatter_read_particles(h_grp, Ndm, Ndm, 0, *gparts);
         break;
 
       default:
@@ -442,10 +441,10 @@ void read_ic_single(char* fileName, double dim[3], struct part** parts,
   }
 
   /* Prepare the DM particles */
-  prepare_dm_gparts(*gparts, Ndm);
+  if (!dry_run) prepare_dm_gparts(*gparts, Ndm);
 
   /* Now duplicate the hydro particle into gparts */
-  duplicate_hydro_gparts(*parts, *gparts, *Ngas, Ndm);
+  if (!dry_run) duplicate_hydro_gparts(*parts, *gparts, *Ngas, Ndm);
 
   /* message("Done Reading particles..."); */
 
diff --git a/src/single_io.h b/src/single_io.h
index c5250280e82e1801b2a4a6136d404d09093dd0ec..587ebe07b6fa2b984b964baf282e7ceb1003ad29 100644
--- a/src/single_io.h
+++ b/src/single_io.h
@@ -28,7 +28,7 @@
 
 void read_ic_single(char* fileName, double dim[3], struct part** parts,
                     struct gpart** gparts, size_t* Ngas, size_t* Ndm,
-                    int* periodic);
+                    int* periodic, int dry_run);
 
 void write_output_single(struct engine* e, struct UnitSystem* us);
 
diff --git a/src/space.c b/src/space.c
index 62db872cc0f436f8025e4c03af008a4b78d2c637..ca93490cb45b50ad0a74d7170813613f972cdef2 100644
--- a/src/space.c
+++ b/src/space.c
@@ -1295,14 +1295,15 @@ struct cell *space_getcell(struct space *s) {
  * @brief Split the space into cells given the array of particles.
  *
  * @param s The #space to initialize.
+ * @param params The parsed parameter file.
  * @param dim Spatial dimensions of the domain.
  * @param parts Array of Gas particles.
  * @param gparts Array of Gravity particles.
- * @param Ngas The number of Gas particles in the space.
+ * @param Npart The number of Gas particles in the space.
  * @param Ngpart The number of Gravity particles in the space.
  * @param periodic flag whether the domain is periodic or not.
- * @param h_max The maximal interaction radius.
  * @param verbose Print messages to stdout or not
+ * @param dry_run If 1, just initialise stuff, don't do anything with the parts.
  *
  * Makes a grid of edge length > r_max and fills the particles
  * into the respective cells. Cells containing more than #space_splitsize
@@ -1310,65 +1311,114 @@ struct cell *space_getcell(struct space *s) {
  * recursively.
  */
 
-void space_init(struct space *s, double dim[3], struct part *parts,
-                struct gpart *gparts, size_t Ngas, size_t Ngpart, int periodic,
-                double h_max, int verbose) {
+void space_init(struct space *s, const struct swift_params *params,
+                double dim[3], struct part *parts, struct gpart *gparts,
+                size_t Npart, size_t Ngpart, int periodic, int verbose,
+                int dry_run) {
+
+  /* Clean-up everything */
+  bzero(s, sizeof(struct space));
 
   /* Store everything in the space. */
   s->dim[0] = dim[0];
   s->dim[1] = dim[1];
   s->dim[2] = dim[2];
   s->periodic = periodic;
-  s->nr_parts = Ngas;
-  s->size_parts = Ngas;
+  s->nr_parts = Npart;
+  s->size_parts = Npart;
   s->parts = parts;
   s->nr_gparts = Ngpart;
   s->size_gparts = Ngpart;
   s->gparts = gparts;
-  s->cell_min = h_max;
-  s->nr_queues = 1;
+  s->cell_min = parser_get_param_double(params, "SPH:max_smoothing_length");
+  s->nr_queues = 1; /* Temporary value until engine construction */
   s->size_parts_foreign = 0;
 
-  /* Check that all the gas particle positions are reasonable, wrap if periodic.
-   */
-  if (periodic) {
-    for (int k = 0; k < Ngas; k++)
-      for (int j = 0; j < 3; j++) {
-        while (parts[k].x[j] < 0) parts[k].x[j] += dim[j];
-        while (parts[k].x[j] >= dim[j]) parts[k].x[j] -= dim[j];
-      }
-  } else {
-    for (int k = 0; k < Ngas; k++)
-      for (int j = 0; j < 3; j++)
-        if (parts[k].x[j] < 0 || parts[k].x[j] >= dim[j])
-          error("Not all particles are within the specified domain.");
+  /* Get the constants for the scheduler */
+  space_maxsize = parser_get_param_int(params, "Scheduler:cell_max_size");
+  space_subsize = parser_get_param_int(params, "Scheduler:cell_sub_size");
+  space_splitsize = parser_get_param_int(params, "Scheduler:cell_split_size");
+  if(verbose)
+    message("max_size set to %d, sub_size set to %d, split_size set to %d",
+	    space_maxsize, space_subsize, space_splitsize);
+
+  /* Check that we have enough cells */
+  if (s->cell_min * 3 > dim[0] || s->cell_min * 3 > dim[1] ||
+      s->cell_min * 3 > dim[2])
+    error(
+        "Maximal smoothing length (%e) too large. Needs to be "
+        "smaller than 1/3 the simulation box size [%e %e %e]",
+        s->cell_min, dim[0], dim[1], dim[2]);
+
+  /* Apply h scaling */
+  const double scaling =
+      parser_get_param_double(params, "InitialConditions:h_scaling");
+  if (scaling != 1.0 && !dry_run) {
+    message("Re-scaling smoothing lengths by a factor %e", scaling);
+    for (size_t k = 0; k < Npart; k++) parts[k].h *= scaling;
   }
 
-  /* Same for the gparts */
-  if (periodic) {
-    for (int k = 0; k < Ngpart; k++)
-      for (int j = 0; j < 3; j++) {
-        while (gparts[k].x[j] < 0) gparts[k].x[j] += dim[j];
-        while (gparts[k].x[j] >= dim[j]) gparts[k].x[j] -= dim[j];
-      }
-  } else {
-    for (int k = 0; k < Ngpart; k++)
-      for (int j = 0; j < 3; j++)
-        if (gparts[k].x[j] < 0 || gparts[k].x[j] >= dim[j])
-          error("Not all particles are within the specified domain.");
+  /* Apply shift */
+  double shift[3] = {0.0, 0.0, 0.0};
+  shift[0] = parser_get_param_double(params, "InitialConditions:shift_x");
+  shift[1] = parser_get_param_double(params, "InitialConditions:shift_y");
+  shift[2] = parser_get_param_double(params, "InitialConditions:shift_z");
+  if ((shift[0] != 0 || shift[1] != 0 || shift[2] != 0) && !dry_run) {
+    message("Shifting particles by [%e %e %e]", shift[0], shift[1], shift[2]);
+    for (size_t k = 0; k < Npart; k++) {
+      parts[k].x[0] += shift[0];
+      parts[k].x[1] += shift[1];
+      parts[k].x[2] += shift[2];
+    }
+    for (size_t k = 0; k < Ngpart; k++) {
+      gparts[k].x[0] += shift[0];
+      gparts[k].x[1] += shift[1];
+      gparts[k].x[2] += shift[2];
+    }
+  }
+
+  if (!dry_run) {
+
+    /* Check that all the part positions are reasonable, wrap if periodic. */
+    if (periodic) {
+      for (int k = 0; k < Npart; k++)
+        for (int j = 0; j < 3; j++) {
+          while (parts[k].x[j] < 0) parts[k].x[j] += dim[j];
+          while (parts[k].x[j] >= dim[j]) parts[k].x[j] -= dim[j];
+        }
+    } else {
+      for (int k = 0; k < Npart; k++)
+        for (int j = 0; j < 3; j++)
+          if (parts[k].x[j] < 0 || parts[k].x[j] >= dim[j])
+            error("Not all particles are within the specified domain.");
+    }
+
+    /* Same for the gparts */
+    if (periodic) {
+      for (int k = 0; k < Ngpart; k++)
+        for (int j = 0; j < 3; j++) {
+          while (gparts[k].x[j] < 0) gparts[k].x[j] += dim[j];
+          while (gparts[k].x[j] >= dim[j]) gparts[k].x[j] -= dim[j];
+        }
+    } else {
+      for (int k = 0; k < Ngpart; k++)
+        for (int j = 0; j < 3; j++)
+          if (gparts[k].x[j] < 0 || gparts[k].x[j] >= dim[j])
+            error("Not all g-particles are within the specified domain.");
+    }
   }
 
   /* Allocate the extra parts array. */
   if (posix_memalign((void *)&s->xparts, xpart_align,
-                     Ngas * sizeof(struct xpart)) != 0)
+                     Npart * sizeof(struct xpart)) != 0)
     error("Failed to allocate xparts.");
-  bzero(s->xparts, Ngas * sizeof(struct xpart));
+  bzero(s->xparts, Npart * sizeof(struct xpart));
 
   /* Init the space lock. */
   if (lock_init(&s->lock) != 0) error("Failed to create space spin-lock.");
 
   /* Build the cells and the tasks. */
-  space_regrid(s, h_max, verbose);
+  if (!dry_run) space_regrid(s, s->cell_min, verbose);
 }
 
 /**
diff --git a/src/space.h b/src/space.h
index 8eb166d82fac2e1a551665680af95c2cf4a188ac..88e2f6f52774651217c4ff24e25f549d8ae1e347 100644
--- a/src/space.h
+++ b/src/space.h
@@ -24,6 +24,7 @@
 
 /* Local includes. */
 #include "cell.h"
+#include "parser.h"
 #include "part.h"
 
 /* Forward-declare the engine to avoid cyclic includes. */
@@ -132,9 +133,10 @@ void space_gparts_sort(struct space *s, int *ind, size_t N, int min, int max,
 struct cell *space_getcell(struct space *s);
 int space_getsid(struct space *s, struct cell **ci, struct cell **cj,
                  double *shift);
-void space_init(struct space *s, double dim[3], struct part *parts,
-                struct gpart *gparts, size_t N, size_t Ngpart, int periodic,
-                double h_max, int verbose);
+void space_init(struct space *s, const struct swift_params *params,
+                double dim[3], struct part *parts, struct gpart *gparts,
+                size_t Npart, size_t Ngpart, int periodic, int verbose,
+                int dry_run);
 void space_map_cells_pre(struct space *s, int full,
                          void (*fun)(struct cell *c, void *data), void *data);
 void space_map_parts(struct space *s,
diff --git a/src/units.c b/src/units.c
index 8c9fd14452e9e1fdfe029ac89d22d7cd43aa0ef7..184dbe8a0df000008dba1d7003558d83b1f08cad 100644
--- a/src/units.c
+++ b/src/units.c
@@ -43,17 +43,24 @@
 
 /**
  * @brief Initialises the UnitSystem structure with the constants given in
- * const.h
- * @param us The UnitSystem to initialize
+ * rhe parameter file.
+ *
+ * @param us The UnitSystem to initialize.
+ * @param params The parsed parameter file.
  */
-
-void initUnitSystem(struct UnitSystem* us) {
-  us->UnitMass_in_cgs = const_unit_mass_in_cgs;
-  us->UnitLength_in_cgs = const_unit_length_in_cgs;
-  us->UnitTime_in_cgs = 1. / ((double)const_unit_velocity_in_cgs /
-                              ((double)const_unit_length_in_cgs));
-  us->UnitCurrent_in_cgs = 1.;
-  us->UnitTemperature_in_cgs = 1.;
+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");
+  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");
 }
 
 /**
@@ -61,7 +68,8 @@ void initUnitSystem(struct UnitSystem* us) {
  * @param us The UnitSystem used
  * @param baseUnit The base unit
  */
-double getBaseUnit(struct UnitSystem* us, enum BaseUnits baseUnit) {
+double units_get_base_unit(const struct UnitSystem* us,
+                           enum BaseUnits baseUnit) {
   switch (baseUnit) {
     case UNIT_MASS:
       return us->UnitMass_in_cgs;
@@ -83,7 +91,7 @@ double getBaseUnit(struct UnitSystem* us, enum BaseUnits baseUnit) {
  * @brief Returns the base unit symbol
  * @param baseUnit The base unit
  */
-const char* getBaseUnitSymbol(enum BaseUnits baseUnit) {
+const char* units_get_base_unit_symbol(enum BaseUnits baseUnit) {
   switch (baseUnit) {
     case UNIT_MASS:
       return "U_M";
@@ -105,7 +113,7 @@ const char* getBaseUnitSymbol(enum BaseUnits baseUnit) {
  * @brief Returns the base unit symbol in the cgs system
  * @param baseUnit The base unit
  */
-const char* getBaseUnitCGSSymbol(enum BaseUnits baseUnit) {
+const char* units_get_base_unit_CGS_symbol(enum BaseUnits baseUnit) {
   switch (baseUnit) {
     case UNIT_MASS:
       return "g";
@@ -123,8 +131,8 @@ const char* getBaseUnitCGSSymbol(enum BaseUnits baseUnit) {
   return "";
 }
 
-void getBaseUnitExponantsArray(float baseUnitsExp[5],
-                               enum UnitConversionFactor unit) {
+void units_get_base_unit_exponants_array(float baseUnitsExp[5],
+                                         enum UnitConversionFactor unit) {
   switch (unit) {
     case UNIT_CONV_NO_UNITS:
       break;
@@ -265,12 +273,13 @@ void getBaseUnitExponantsArray(float baseUnitsExp[5],
  * @param us The system of units in use
  * @param unit The unit to convert
  */
-double conversionFactor(struct UnitSystem* us, enum UnitConversionFactor unit) {
+double units_conversion_factor(const struct UnitSystem* us,
+                               enum UnitConversionFactor unit) {
   float baseUnitsExp[5] = {0.f};
 
-  getBaseUnitExponantsArray(baseUnitsExp, unit);
+  units_get_base_unit_exponants_array(baseUnitsExp, unit);
 
-  return generalConversionFactor(us, baseUnitsExp);
+  return units_general_conversion_factor(us, baseUnitsExp);
 }
 
 /**
@@ -278,12 +287,13 @@ double conversionFactor(struct UnitSystem* us, enum UnitConversionFactor unit) {
  * @param us The system of units in use
  * @param unit The unit to convert
  */
-float hFactor(struct UnitSystem* us, enum UnitConversionFactor unit) {
+float units_h_factor(const struct UnitSystem* us,
+                     enum UnitConversionFactor unit) {
   float baseUnitsExp[5] = {0.f};
 
-  getBaseUnitExponantsArray(baseUnitsExp, unit);
+  units_get_base_unit_exponants_array(baseUnitsExp, unit);
 
-  return generalhFactor(us, baseUnitsExp);
+  return units_general_h_factor(us, baseUnitsExp);
 }
 
 /**
@@ -291,25 +301,26 @@ float hFactor(struct UnitSystem* us, enum UnitConversionFactor unit) {
  * @param us The system of units in use
  * @param unit The unit to convert
  */
-float aFactor(struct UnitSystem* us, enum UnitConversionFactor unit) {
+float units_a_factor(const struct UnitSystem* us,
+                     enum UnitConversionFactor unit) {
   float baseUnitsExp[5] = {0.f};
 
-  getBaseUnitExponantsArray(baseUnitsExp, unit);
+  units_get_base_unit_exponants_array(baseUnitsExp, unit);
 
-  return generalaFactor(us, baseUnitsExp);
+  return units_general_a_factor(us, baseUnitsExp);
 }
 
 /**
  * @brief Returns a string containing the exponents of the base units making up
  * the conversion factors
  */
-void conversionString(char* buffer, struct UnitSystem* us,
-                      enum UnitConversionFactor unit) {
+void units_conversion_string(char* buffer, const struct UnitSystem* us,
+                             enum UnitConversionFactor unit) {
   float baseUnitsExp[5] = {0.f};
 
-  getBaseUnitExponantsArray(baseUnitsExp, unit);
+  units_get_base_unit_exponants_array(baseUnitsExp, unit);
 
-  generalConversionString(buffer, us, baseUnitsExp);
+  units_general_conversion_string(buffer, us, baseUnitsExp);
 }
 
 /**
@@ -319,14 +330,14 @@ void conversionString(char* buffer, struct UnitSystem* us,
  * @param baseUnitsExponants The exponent of each base units required to form
  * the desired quantity. See conversionFactor() for a working example
  */
-double generalConversionFactor(struct UnitSystem* us,
-                               float baseUnitsExponants[5]) {
+double units_general_conversion_factor(const struct UnitSystem* us,
+                                       float baseUnitsExponants[5]) {
   double factor = 1.;
   int i;
 
   for (i = 0; i < 5; ++i)
     if (baseUnitsExponants[i] != 0)
-      factor *= pow(getBaseUnit(us, i), baseUnitsExponants[i]);
+      factor *= pow(units_get_base_unit(us, i), baseUnitsExponants[i]);
   return factor;
 }
 
@@ -337,7 +348,8 @@ double generalConversionFactor(struct UnitSystem* us,
  * @param baseUnitsExponants The exponent of each base units required to form
  * the desired quantity. See conversionFactor() for a working example
  */
-float generalhFactor(struct UnitSystem* us, float baseUnitsExponants[5]) {
+float units_general_h_factor(const struct UnitSystem* us,
+                             float baseUnitsExponants[5]) {
   float factor_exp = 0.f;
 
   factor_exp += -baseUnitsExponants[UNIT_MASS];
@@ -354,7 +366,8 @@ float generalhFactor(struct UnitSystem* us, float baseUnitsExponants[5]) {
  * @param baseUnitsExponants The exponent of each base units required to form
  * the desired quantity. See conversionFactor() for a working example
  */
-float generalaFactor(struct UnitSystem* us, float baseUnitsExponants[5]) {
+float units_general_a_factor(const struct UnitSystem* us,
+                             float baseUnitsExponants[5]) {
   float factor_exp = 0.f;
 
   factor_exp += baseUnitsExponants[UNIT_LENGTH];
@@ -371,11 +384,11 @@ float generalaFactor(struct UnitSystem* us, float baseUnitsExponants[5]) {
  * @param baseUnitsExponants The exponent of each base units required to form
  * the desired quantity. See conversionFactor() for a working example
  */
-void generalConversionString(char* buffer, struct UnitSystem* us,
-                             float baseUnitsExponants[5]) {
+void units_general_conversion_string(char* buffer, const struct UnitSystem* us,
+                                     float baseUnitsExponants[5]) {
   char temp[14];
-  double a_exp = generalaFactor(us, baseUnitsExponants);
-  double h_exp = generalhFactor(us, baseUnitsExponants);
+  double a_exp = units_general_a_factor(us, baseUnitsExponants);
+  double h_exp = units_general_h_factor(us, baseUnitsExponants);
   int i;
 
   /* Check whether we are unitless or not */
@@ -415,12 +428,13 @@ void generalConversionString(char* buffer, struct UnitSystem* us,
       if (baseUnitsExponants[i] == 0.)
         sprintf(temp, " ");
       else if (baseUnitsExponants[i] == 1.)
-        sprintf(temp, "%s ", getBaseUnitSymbol(i));
+        sprintf(temp, "%s ", units_get_base_unit_symbol(i));
       else if (remainder(baseUnitsExponants[i], 1.) == 0)
-        sprintf(temp, "%s^%d ", getBaseUnitSymbol(i),
+        sprintf(temp, "%s^%d ", units_get_base_unit_symbol(i),
                 (int)baseUnitsExponants[i]);
       else
-        sprintf(temp, "%s^%7.4f ", getBaseUnitSymbol(i), baseUnitsExponants[i]);
+        sprintf(temp, "%s^%7.4f ", units_get_base_unit_symbol(i),
+                baseUnitsExponants[i]);
       strncat(buffer, temp, 12);
     }
 
@@ -432,12 +446,12 @@ void generalConversionString(char* buffer, struct UnitSystem* us,
       if (baseUnitsExponants[i] == 0.)
         continue;
       else if (baseUnitsExponants[i] == 1.)
-        sprintf(temp, "%s ", getBaseUnitCGSSymbol(i));
+        sprintf(temp, "%s ", units_get_base_unit_CGS_symbol(i));
       else if (remainder(baseUnitsExponants[i], 1.) == 0)
-        sprintf(temp, "%s^%d ", getBaseUnitCGSSymbol(i),
+        sprintf(temp, "%s^%d ", units_get_base_unit_CGS_symbol(i),
                 (int)baseUnitsExponants[i]);
       else
-        sprintf(temp, "%s^%7.4f ", getBaseUnitCGSSymbol(i),
+        sprintf(temp, "%s^%7.4f ", units_get_base_unit_CGS_symbol(i),
                 baseUnitsExponants[i]);
       strncat(buffer, temp, 12);
     }
diff --git a/src/units.h b/src/units.h
index 1b977529784c1ef3069e1e932b16fd0b87073786..3e349dc16787cd4052a3e9205b21dce3c3732448 100644
--- a/src/units.h
+++ b/src/units.h
@@ -19,6 +19,12 @@
 #ifndef SWIFT_UNITS_H
 #define SWIFT_UNITS_H
 
+/* Config parameters. */
+#include "../config.h"
+
+/* Local includes. */
+#include "parser.h"
+
 /**
  * @brief The unit system used internally.
  *
@@ -86,74 +92,25 @@ enum UnitConversionFactor {
   UNIT_CONV_TEMPERATURE
 };
 
-/**
- * @brief Initialises the UnitSystem structure with the constants given in
- * const.h
- */
-void initUnitSystem(struct UnitSystem*);
-
-/**
- * @brief Returns the base unit conversion factor for a given unit system
- */
-double getBaseUnit(struct UnitSystem*, enum BaseUnits);
-
-/**
- * @brief Returns the base unit symbol in the cgs system
- */
-const char* getBaseUnitSymbol(enum BaseUnits);
-
-/**
- * @brief Returns the base unit symbol in the cgs system
- */
-const char* getBaseUnitCGSSymbol(enum BaseUnits);
-
-/**
- * @brief Returns the conversion factor for a given unit (expressed in terms of
- * the 5 fundamental units) in the chosen unit system
- */
-double generalConversionFactor(struct UnitSystem* us,
-                               float baseUnitsExponants[5]);
-
-/**
- * @brief Returns the conversion factor for a given unit in the chosen unit
- * system
- */
-double conversionFactor(struct UnitSystem* us, enum UnitConversionFactor unit);
-
-/**
- * @brief Returns the h factor for a given unit (expressed in terms of the 5
- * fundamental units) in the chosen unit system
- */
-float generalhFactor(struct UnitSystem* us, float baseUnitsExponants[5]);
-
-/**
- * @brief Returns the h factor for a given unit in the chosen unit system
- */
-float hFactor(struct UnitSystem* us, enum UnitConversionFactor unit);
-
-/**
- * @brief Returns the scaling factor for a given unit (expressed in terms of the
- * 5 fundamental units) in the chosen unit system
- */
-float generalaFactor(struct UnitSystem* us, float baseUnitsExponants[5]);
-
-/**
- * @brief Returns the scaling factor for a given unit in the chosen unit system
- */
-float aFactor(struct UnitSystem* us, enum UnitConversionFactor unit);
-
-/**
- * @brief Returns a string containing the exponents of the base units making up
- * the conversion factors (expressed in terms of the 5 fundamental units)
- */
-void generalConversionString(char* buffer, struct UnitSystem* us,
+void units_init(struct UnitSystem*, const struct swift_params*);
+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);
+double units_general_conversion_factor(const struct UnitSystem* us,
+                                       float baseUnitsExponants[5]);
+double units_conversion_factor(const struct UnitSystem* us,
+                               enum UnitConversionFactor unit);
+float units_general_h_factor(const struct UnitSystem* us,
                              float baseUnitsExponants[5]);
-
-/**
- * @brief Returns a string containing the exponents of the base units making up
- * the conversion factors
- */
-void conversionString(char* buffer, struct UnitSystem* us,
-                      enum UnitConversionFactor unit);
+float units_h_factor(const struct UnitSystem* us,
+                     enum UnitConversionFactor unit);
+float units_general_a_factor(const struct UnitSystem* us,
+                             float baseUnitsExponants[5]);
+float units_a_factor(const struct UnitSystem* us,
+                     enum UnitConversionFactor unit);
+void units_general_conversion_string(char* buffer, const struct UnitSystem* us,
+                                     float baseUnitsExponants[5]);
+void units_conversion_string(char* buffer, const struct UnitSystem* us,
+                             enum UnitConversionFactor unit);
 
 #endif /* SWIFT_UNITS_H */
diff --git a/src/version.c b/src/version.c
index 6aeee2d8bcbc4652f679bbb786e9e512ebc4caa6..27841a16019a69442e66b21c327f4241e440fb12 100644
--- a/src/version.c
+++ b/src/version.c
@@ -241,7 +241,7 @@ const char *metis_version(void) {
  */
 void greetings(void) {
 
-  printf(" Welcome to the cosmological code\n");
+  printf(" Welcome to the cosmological hydrodynamical code\n");
   printf("    ______       _________________\n");
   printf("   / ___/ |     / /  _/ ___/_  __/\n");
   printf("   \\__ \\| | /| / // // /_   / /   \n");
diff --git a/tests/testParser.c b/tests/testParser.c
index a4b8789fca056fef659bca78eae9d0effb2ceb66..0b08d20c9e2d48de1858877cf186eaa9d0ac84c0 100644
--- a/tests/testParser.c
+++ b/tests/testParser.c
@@ -20,48 +20,51 @@
 #include "parser.h"
 #include <assert.h>
 #include <string.h>
+#include <stdio.h>
 #include <math.h>
 
 int main(int argc, char *argv[]) {
-
   const char *input_file = argv[1];
 
   /* Create a structure to read file into. */
   struct swift_params param_file;
 
-  /* Create variables that will be set from the parameter file. */
-  int no_of_threads = 0;
-  int no_of_time_steps = 0;
-  float max_h = 0.0f;
-  double start_time = 0.0;
-  char ic_file[PARSER_MAX_LINE_SIZE];
-
   /* Read the parameter file. */
   parser_read_file(input_file, &param_file);
 
-  /* Print the contents of the structure. */
+  /* Print the contents of the structure to stdout. */
   parser_print_params(&param_file);
 
+  /* Print the contents of the structure to a file in YAML format. */
+  parser_write_params_to_file(&param_file, "parser_output.yml");
+
   /* Retrieve parameters and store them in variables defined above.
    * Have to specify the name of the parameter as it appears in the
    * input file: testParserInput.yaml.*/
-  parser_get_param_int(&param_file, "no_of_threads", &no_of_threads);
-  parser_get_param_int(&param_file, "no_of_time_steps", &no_of_time_steps);
-  parser_get_param_float(&param_file, "max_h", &max_h);
-  parser_get_param_double(&param_file, "start_time", &start_time);
-  parser_get_param_string(&param_file, "ic_file", ic_file);
+  const int no_of_threads =
+      parser_get_param_int(&param_file, "Scheduler:no_of_threads");
+  const int no_of_time_steps =
+      parser_get_param_int(&param_file, "Simulation:no_of_time_steps");
+  const float max_h = parser_get_param_float(&param_file, "Simulation:max_h");
+  const double start_time =
+      parser_get_param_double(&param_file, "Simulation:start_time");
+  const int kernel = parser_get_param_int(&param_file, "kernel");
+
+  char ic_file[PARSER_MAX_LINE_SIZE];
+  parser_get_param_string(&param_file, "IO:ic_file", ic_file);
 
   /* Print the variables to check their values are correct. */
   printf(
       "no_of_threads: %d, no_of_time_steps: %d, max_h: %f, start_time: %lf, "
-      "ic_file: %s\n",
-      no_of_threads, no_of_time_steps, max_h, start_time, ic_file);
+      "ic_file: %s, kernel: %d\n",
+      no_of_threads, no_of_time_steps, max_h, start_time, ic_file, kernel);
 
   assert(no_of_threads == 16);
   assert(no_of_time_steps == 10);
   assert(fabs(max_h - 1.1255) < 0.00001);
   assert(fabs(start_time - 1.23456789) < 0.00001);
   assert(strcmp(ic_file, "ic_file.ini") == 0); /*strcmp returns 0 if correct.*/
+  assert(kernel == 4);
 
   return 0;
 }
diff --git a/tests/testParser.sh b/tests/testParser.sh
index 3dad7f386f792ff2beb6e94eb093bad4085023a4..53d2bbe4e0230032666ace228449f913f03e0464 100755
--- a/tests/testParser.sh
+++ b/tests/testParser.sh
@@ -1,3 +1,4 @@
 #!/bin/bash
 
+rm parser_output.yml
 ./testParser testParserInput.yaml
diff --git a/tests/testParserInput.yaml b/tests/testParserInput.yaml
index d695e6a8ddd327e31224f36a6e34767ea8d36408..c7fefb3242ab4e140756789aad9979d024f83906 100644
--- a/tests/testParserInput.yaml
+++ b/tests/testParserInput.yaml
@@ -1,9 +1,19 @@
 ---
-no_of_threads:      16 # The number of threads that will be used. 
-no_of_time_steps:   10
-max_h:              1.1255
-start_time:         1.23456789
-#Input file
-ic_file:            ic_file.ini
+#section_1:
+#  var_a: 4.5e10
+#  var_b: Hello World!
 
+Scheduler:
+  no_of_threads:      16 # The number of threads that will be used. 
+
+kernel: 4
+
+Simulation:    
+  no_of_time_steps:   10
+  max_h:              1.1255
+  start_time:         1.23456789
+
+IO:
+  #Input file
+  ic_file:            ic_file.ini
 ...
diff --git a/tests/testReading.c b/tests/testReading.c
index 9dda4c7bad75d35a8a93e0c2acb0619409a91afd..33aeb5095ba499bc0fd18ba15b513e351692432e 100644
--- a/tests/testReading.c
+++ b/tests/testReading.c
@@ -35,7 +35,8 @@ int main() {
   const double rho = 2.;
 
   /* Read data */
-  read_ic_single("input.hdf5", dim, &parts, &gparts, &Ngas, &Ngpart, &periodic);
+  read_ic_single("input.hdf5", dim, &parts, &gparts, &Ngas, &Ngpart, &periodic,
+                 0);
 
   /* Check global properties read are correct */
   assert(dim[0] == boxSize);
diff --git a/tests/testSingle.c b/tests/testSingle.c
index 8771fba0c1912905d3936562fa9dad0223d89220..eb49a570b93b14734c9e6af37d3d8a2b90d04078 100644
--- a/tests/testSingle.c
+++ b/tests/testSingle.c
@@ -91,8 +91,8 @@ int main(int argc, char *argv[]) {
   p2.force.POrho2 = p2.u * (const_hydro_gamma - 1.0f) / p2.rho;
 
   /* Dump a header. */
-  //printParticle_single(&p1, NULL);
-  //printParticle_single(&p2, NULL);
+  // printParticle_single(&p1, NULL);
+  // printParticle_single(&p2, NULL);
   printf("# r a_1 udt_1 a_2 udt_2\n");
 
   /* Loop over the different radii. */