diff --git a/doc/RTD/source/ParameterFiles/output_selection.rst b/doc/RTD/source/ParameterFiles/output_selection.rst
index 8100c97abdc8f2839cffa29c9683bcf68bea00f0..41ae9cbb4f366b99e36da1dcc043cbd01485f593 100644
--- a/doc/RTD/source/ParameterFiles/output_selection.rst
+++ b/doc/RTD/source/ParameterFiles/output_selection.rst
@@ -60,34 +60,89 @@ CGS. Entries in the file look like:
 
 .. code:: YAML
 
-   SelectOutput:
-     # Particle Type Gas
-     Coordinates_Gas:      1  # Co-moving positions of the particles. ::: Conversion to physical CGS: a U_L  [ cm ]
-     Velocities_Gas:       1  # Peculiar velocities of the stars. This is (a * dx/dt) where x is the co-moving positions of the particles. ::: Conversion to physical CGS: U_L U_t^-1  [ cm s^-1 ]
-     Masses_Gas:           1  # Masses of the particles. ::: Conversion to physical CGS: U_M  [ g ]
-     SmoothingLengths_Gas: 1  # Co-moving smoothing lengths (FWHM of the kernel) of the particles. ::: Conversion to physical CGS: a U_L  [ cm ]
-
-Users can select the particle fields to output in snapshot using the
-YAML parameter file.  In section ``SelectOutput``, users can demand to
-remove a field by adding a parameter formatted in the following way
-``field_parttype`` where ``field`` is the name of the field that is to
-be removed (e.g. ``Masses``) and ``parttype`` is the type of particles
-that contains this field (``Gas``, ``DM``, ``Stars`` or ``BH``).  For
-a parameter, the only values accepted are ``0`` (skip this field when
-writing) or ``1`` (default, do not skip this field when writing). By
-default all fields are written.
-
-This field is mostly used to remove unnecessary output by listing them
-with 0's. A classic use-case for this feature is a DM-only simulation
-(pure n-body) where all particles have the same mass. Outputting the
-mass field in the snapshots results in extra i/o time and unnecessary
-waste of disk space. The corresponding section of the ``yaml``
-parameter file would look like:
+  Default:
+    # Particle Type Gas
+    Coordinates_Gas: off  # Co-moving positions of the particles : a U_L  [ cm ]
+    Velocities_Gas: on  # Peculiar velocities of the stars. This is (a * dx/dt) where x is the co-moving positions of the particles : U_L U_t^-1  [ cm s^-1 ]
+    Masses_Gas: on  # Masses of the particles : U_M  [ g ]
+    SmoothingLengths_Gas: on  # Co-moving smoothing lengths (FWHM of the kernel) of the particles : a U_L  [ cm ]
+    ...
+
+Users can select the particle fields to output in snapshot using a (separate)
+YAML parameter file. By default, you can define a section `Default` at the
+top level of this file (in the exact same way as the file dumped by using the
+`-o` option in SWIFT). By default, all fields are written, but by using the
+"off" string, you can force the code to skip over unwanted outputs.
+
+You must then, in the regular SWIFT parameter file, select the following
+options:
 
 .. code:: YAML
-	  
-  SelectOutput:
-    Masses_DM:   0
+
+  Snapshots:
+    select_output_on: 1
+    select_output: your_select_output_yaml.yml
+
+This field is mostly used to remove unnecessary output by listing them with
+0's. A classic use-case for this feature is a DM-only simulation (pure
+n-body) where all particles have the same mass. Outputting the mass field in
+the snapshots results in extra i/o time and unnecessary waste of disk space.
+The corresponding section of the YAML file would look like:
+
+.. code:: YAML
+
+  Default:
+    Masses_DM: off
 
 Entries can simply be copied from the ``output.yml`` generated by the
 ``-o`` runtime flag. 
+
+
+Combining Output Lists and Output Selection
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+It is possible to combine the behaviour of the output list and the select
+output file. To do so, you will need to enable both the ``select_output`` and
+``output_list`` options in your main ``parameter_file.yml`` as follows:
+
+.. code:: YAML
+
+  Snapshots:
+    output_list_on: 1
+    output_list: "output_list.txt"
+    select_output_on: 1
+    select_output: "select_output.yml"
+
+A typical use case for such a scenario is the dumping of 'snapshots' and
+so-called 'snipshots', containing less information than their full snapshot
+cousins. To do this, we will define two top-level sections in our
+``select_output.yml`` file as follows:
+
+.. code:: YAML
+
+  # Only turn off DM masses in snapshots, everything else is turned on
+  Snapshot:
+    Masses_DM: off
+
+  # Turn off lots of stuff in snipshots!
+  Snipshot:
+    Metal_Mass_Fractions_Gas: off
+    Element_Mass_Fractions_Gas: off
+    ...
+
+To then select which outputs are 'snapshots' and which are 'snipshots', you
+will need to add the ``Select Output`` column to the ``output_list.txt`` as
+follows::
+
+  # Redshift, Select Output
+  100.0, Snapshot
+  90.0, Snipshot
+  80.0, Snipshot
+  70.0, Snipshot
+  60.0, Snapshot
+  ...
+
+This will enable your simulation to perform partial dumps only at the outputs
+labelled as ``Snipshot``. The name of the output selection that corresponds
+to your choice in the output list will be written to the snapshot header as
+``Header/SelectOutput``.
diff --git a/examples/SmallCosmoVolume/SmallCosmoVolume_Snipshots/README b/examples/SmallCosmoVolume/SmallCosmoVolume_Snipshots/README
new file mode 100644
index 0000000000000000000000000000000000000000..9be4c76f304eff3962d9ba19c2159ac72dde7e0a
--- /dev/null
+++ b/examples/SmallCosmoVolume/SmallCosmoVolume_Snipshots/README
@@ -0,0 +1,17 @@
+Small LCDM cosmological simulation generated by C. Power. Cosmology
+is WMAP9 and the box is 100Mpc/h in size with 64^3 particles.
+We use a softening length of 1/25th of the mean inter-particle separation.
+
+The ICs have been generated to run with Gadget-2 so we need to switch
+on the options to cancel the h-factors and a-factors at reading time.
+We generate gas from the ICs using SWIFT's internal mechanism and set the
+temperature to the expected gas temperature at this redshift.
+
+The 'plotTempEvolution.py' plots the temperature evolution of the gas
+in the simulated volume.
+
+This version uses an output list that has 'snapshots' and 'snipshots'
+as a useful example for this functionality.
+
+MD5 checksum of the ICs:
+08736c3101fd738e22f5159f78e6022b  small_cosmo_volume.hdf5
diff --git a/examples/SmallCosmoVolume/SmallCosmoVolume_Snipshots/getIC.sh b/examples/SmallCosmoVolume/SmallCosmoVolume_Snipshots/getIC.sh
new file mode 100755
index 0000000000000000000000000000000000000000..3b8136cc5aca00a25792655c6c505cfeeb0f2bc9
--- /dev/null
+++ b/examples/SmallCosmoVolume/SmallCosmoVolume_Snipshots/getIC.sh
@@ -0,0 +1,3 @@
+#!/bin/bash
+wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/ICs/small_cosmo_volume.hdf5
+
diff --git a/examples/SmallCosmoVolume/SmallCosmoVolume_Snipshots/output_list.txt b/examples/SmallCosmoVolume/SmallCosmoVolume_Snipshots/output_list.txt
new file mode 100644
index 0000000000000000000000000000000000000000..89c07c7b908a388c70982f513b217854884abb6c
--- /dev/null
+++ b/examples/SmallCosmoVolume/SmallCosmoVolume_Snipshots/output_list.txt
@@ -0,0 +1,252 @@
+# Scale Factor, Select Output
+0.019607843137254943, Snapshot
+0.02352941176470591, Snipshot
+0.027450980392156876, Snipshot
+0.03137254901960784, Snipshot
+0.03529411764705881, Snipshot
+0.039215686274509776, Snapshot
+0.04313725490196085, Snipshot
+0.04705882352941182, Snipshot
+0.050980392156862786, Snipshot
+0.05490196078431375, Snipshot
+0.05882352941176472, Snapshot
+0.06274509803921569, Snipshot
+0.06666666666666665, Snipshot
+0.07058823529411762, Snipshot
+0.07450980392156858, Snipshot
+0.07843137254901966, Snapshot
+0.08235294117647063, Snipshot
+0.0862745098039216, Snipshot
+0.09019607843137256, Snipshot
+0.09411764705882353, Snipshot
+0.0980392156862745, Snapshot
+0.10196078431372546, Snipshot
+0.10588235294117654, Snipshot
+0.1098039215686275, Snipshot
+0.11372549019607847, Snipshot
+0.11764705882352944, Snapshot
+0.1215686274509804, Snipshot
+0.12549019607843137, Snipshot
+0.12941176470588234, Snipshot
+0.1333333333333333, Snipshot
+0.13725490196078427, Snapshot
+0.14117647058823535, Snipshot
+0.14509803921568631, Snipshot
+0.14901960784313728, Snipshot
+0.15294117647058825, Snipshot
+0.1568627450980392, Snapshot
+0.16078431372549018, Snipshot
+0.16470588235294115, Snipshot
+0.16862745098039222, Snipshot
+0.1725490196078432, Snipshot
+0.17647058823529416, Snapshot
+0.18039215686274512, Snipshot
+0.1843137254901961, Snipshot
+0.18823529411764706, Snipshot
+0.19215686274509802, Snipshot
+0.196078431372549, Snapshot
+0.19999999999999996, Snipshot
+0.20392156862745103, Snipshot
+0.207843137254902, Snipshot
+0.21176470588235297, Snipshot
+0.21568627450980393, Snapshot
+0.2196078431372549, Snipshot
+0.22352941176470587, Snipshot
+0.22745098039215683, Snipshot
+0.2313725490196079, Snipshot
+0.23529411764705888, Snapshot
+0.23921568627450984, Snipshot
+0.2431372549019608, Snipshot
+0.24705882352941178, Snipshot
+0.25098039215686274, Snipshot
+0.2549019607843137, Snapshot
+0.2588235294117647, Snipshot
+0.26274509803921564, Snipshot
+0.2666666666666667, Snipshot
+0.2705882352941177, Snipshot
+0.27450980392156865, Snapshot
+0.2784313725490196, Snipshot
+0.2823529411764706, Snipshot
+0.28627450980392155, Snipshot
+0.2901960784313725, Snipshot
+0.2941176470588236, Snapshot
+0.29803921568627456, Snipshot
+0.3019607843137255, Snipshot
+0.3058823529411765, Snipshot
+0.30980392156862746, Snipshot
+0.3137254901960784, Snapshot
+0.3176470588235294, Snipshot
+0.32156862745098036, Snipshot
+0.3254901960784313, Snipshot
+0.3294117647058824, Snipshot
+0.33333333333333337, Snapshot
+0.33725490196078434, Snipshot
+0.3411764705882353, Snipshot
+0.34509803921568627, Snipshot
+0.34901960784313724, Snipshot
+0.3529411764705882, Snapshot
+0.3568627450980393, Snipshot
+0.36078431372549025, Snipshot
+0.3647058823529412, Snipshot
+0.3686274509803922, Snipshot
+0.37254901960784315, Snapshot
+0.3764705882352941, Snipshot
+0.3803921568627451, Snipshot
+0.38431372549019605, Snipshot
+0.388235294117647, Snipshot
+0.3921568627450981, Snapshot
+0.39607843137254906, Snipshot
+0.4, Snipshot
+0.403921568627451, Snipshot
+0.40784313725490196, Snipshot
+0.4117647058823529, Snapshot
+0.4156862745098039, Snipshot
+0.41960784313725497, Snipshot
+0.42352941176470593, Snipshot
+0.4274509803921569, Snipshot
+0.43137254901960786, Snapshot
+0.43529411764705883, Snipshot
+0.4392156862745098, Snipshot
+0.44313725490196076, Snipshot
+0.44705882352941173, Snipshot
+0.4509803921568627, Snapshot
+0.4549019607843138, Snipshot
+0.45882352941176474, Snipshot
+0.4627450980392157, Snipshot
+0.4666666666666667, Snipshot
+0.47058823529411764, Snapshot
+0.4745098039215686, Snipshot
+0.4784313725490196, Snipshot
+0.48235294117647065, Snipshot
+0.4862745098039216, Snipshot
+0.4901960784313726, Snapshot
+0.49411764705882355, Snipshot
+0.4980392156862745, Snipshot
+0.5019607843137255, Snipshot
+0.5058823529411764, Snipshot
+0.5098039215686274, Snapshot
+0.5137254901960784, Snipshot
+0.5176470588235293, Snipshot
+0.5215686274509804, Snipshot
+0.5254901960784314, Snipshot
+0.5294117647058824, Snapshot
+0.5333333333333333, Snipshot
+0.5372549019607843, Snipshot
+0.5411764705882354, Snipshot
+0.5450980392156863, Snipshot
+0.5490196078431373, Snapshot
+0.5529411764705883, Snipshot
+0.5568627450980392, Snipshot
+0.5607843137254902, Snipshot
+0.5647058823529412, Snipshot
+0.5686274509803921, Snapshot
+0.5725490196078431, Snipshot
+0.5764705882352941, Snipshot
+0.580392156862745, Snipshot
+0.5843137254901961, Snipshot
+0.5882352941176471, Snapshot
+0.592156862745098, Snipshot
+0.596078431372549, Snipshot
+0.6, Snipshot
+0.603921568627451, Snipshot
+0.607843137254902, Snapshot
+0.611764705882353, Snipshot
+0.615686274509804, Snipshot
+0.6196078431372549, Snipshot
+0.6235294117647059, Snipshot
+0.6274509803921569, Snapshot
+0.6313725490196078, Snipshot
+0.6352941176470588, Snipshot
+0.6392156862745098, Snipshot
+0.6431372549019607, Snipshot
+0.6470588235294118, Snapshot
+0.6509803921568628, Snipshot
+0.6549019607843137, Snipshot
+0.6588235294117647, Snipshot
+0.6627450980392157, Snipshot
+0.6666666666666667, Snapshot
+0.6705882352941177, Snipshot
+0.6745098039215687, Snipshot
+0.6784313725490196, Snipshot
+0.6823529411764706, Snipshot
+0.6862745098039216, Snapshot
+0.6901960784313725, Snipshot
+0.6941176470588235, Snipshot
+0.6980392156862745, Snipshot
+0.7019607843137254, Snipshot
+0.7058823529411764, Snapshot
+0.7098039215686275, Snipshot
+0.7137254901960784, Snipshot
+0.7176470588235294, Snipshot
+0.7215686274509804, Snipshot
+0.7254901960784313, Snapshot
+0.7294117647058824, Snipshot
+0.7333333333333334, Snipshot
+0.7372549019607844, Snipshot
+0.7411764705882353, Snipshot
+0.7450980392156863, Snapshot
+0.7490196078431373, Snipshot
+0.7529411764705882, Snipshot
+0.7568627450980392, Snipshot
+0.7607843137254902, Snipshot
+0.7647058823529411, Snapshot
+0.7686274509803921, Snipshot
+0.7725490196078432, Snipshot
+0.7764705882352941, Snipshot
+0.7803921568627451, Snipshot
+0.7843137254901961, Snapshot
+0.788235294117647, Snipshot
+0.7921568627450981, Snipshot
+0.7960784313725491, Snipshot
+0.8, Snipshot
+0.803921568627451, Snapshot
+0.807843137254902, Snipshot
+0.8117647058823529, Snipshot
+0.8156862745098039, Snipshot
+0.8196078431372549, Snipshot
+0.8235294117647058, Snapshot
+0.8274509803921568, Snipshot
+0.8313725490196078, Snipshot
+0.8352941176470589, Snipshot
+0.8392156862745098, Snipshot
+0.8431372549019608, Snapshot
+0.8470588235294118, Snipshot
+0.8509803921568627, Snipshot
+0.8549019607843138, Snipshot
+0.8588235294117648, Snipshot
+0.8627450980392157, Snapshot
+0.8666666666666667, Snipshot
+0.8705882352941177, Snipshot
+0.8745098039215686, Snipshot
+0.8784313725490196, Snipshot
+0.8823529411764706, Snapshot
+0.8862745098039215, Snipshot
+0.8901960784313725, Snipshot
+0.8941176470588236, Snipshot
+0.8980392156862745, Snipshot
+0.9019607843137255, Snapshot
+0.9058823529411765, Snipshot
+0.9098039215686274, Snipshot
+0.9137254901960784, Snipshot
+0.9176470588235294, Snipshot
+0.9215686274509804, Snapshot
+0.9254901960784314, Snipshot
+0.9294117647058824, Snipshot
+0.9333333333333333, Snipshot
+0.9372549019607843, Snipshot
+0.9411764705882353, Snapshot
+0.9450980392156862, Snipshot
+0.9490196078431372, Snipshot
+0.9529411764705882, Snipshot
+0.9568627450980393, Snipshot
+0.9607843137254902, Snapshot
+0.9647058823529412, Snipshot
+0.9686274509803922, Snipshot
+0.9725490196078431, Snipshot
+0.9764705882352941, Snipshot
+0.9803921568627451, Snapshot
+0.9843137254901961, Snipshot
+0.9882352941176471, Snipshot
+0.9921568627450981, Snipshot
+0.996078431372549, Snipshot
+1.0, Snapshot
diff --git a/examples/SmallCosmoVolume/SmallCosmoVolume_Snipshots/plotTempEvolution.py b/examples/SmallCosmoVolume/SmallCosmoVolume_Snipshots/plotTempEvolution.py
new file mode 100644
index 0000000000000000000000000000000000000000..d707f70450471f2d2fc589dbc382366280e0e7f3
--- /dev/null
+++ b/examples/SmallCosmoVolume/SmallCosmoVolume_Snipshots/plotTempEvolution.py
@@ -0,0 +1,182 @@
+################################################################################
+# This file is part of SWIFT.
+# Copyright (c) 2018 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU Lesser General Public License as published
+# by the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public License
+# along with this program.  If not, see <http://www.gnu.org/licenses/>.
+#
+################################################################################
+
+# Computes the temperature evolution of the gas in a cosmological box
+
+# Physical constants needed for internal energy to temperature conversion
+k_in_J_K = 1.38064852e-23
+mH_in_kg = 1.6737236e-27
+
+# Number of snapshots generated
+n_snapshots = 200
+
+import matplotlib
+matplotlib.use("Agg")
+from pylab import *
+import h5py
+import os.path
+
+# Plot parameters
+params = {'axes.labelsize': 10,
+'axes.titlesize': 10,
+'font.size': 9,
+'legend.fontsize': 9,
+'xtick.labelsize': 10,
+'ytick.labelsize': 10,
+'text.usetex': True,
+ 'figure.figsize' : (3.15,3.15),
+'figure.subplot.left'    : 0.14,
+'figure.subplot.right'   : 0.99,
+'figure.subplot.bottom'  : 0.12,
+'figure.subplot.top'     : 0.99,
+'figure.subplot.wspace'  : 0.15,
+'figure.subplot.hspace'  : 0.12,
+'lines.markersize' : 6,
+'lines.linewidth' : 2.,
+'text.latex.unicode': True
+}
+rcParams.update(params)
+rc('font',**{'family':'sans-serif','sans-serif':['Times']})
+
+# Read the simulation data
+sim = h5py.File("snap_0000.hdf5", "r")
+boxSize = sim["/Header"].attrs["BoxSize"][0]
+time = sim["/Header"].attrs["Time"][0]
+scheme = sim["/HydroScheme"].attrs["Scheme"][0]
+kernel = sim["/HydroScheme"].attrs["Kernel function"][0]
+neighbours = sim["/HydroScheme"].attrs["Kernel target N_ngb"][0]
+eta = sim["/HydroScheme"].attrs["Kernel eta"][0]
+alpha = sim["/HydroScheme"].attrs["Alpha viscosity"][0]
+H_mass_fraction = sim["/HydroScheme"].attrs["Hydrogen mass fraction"][0]
+H_transition_temp = sim["/HydroScheme"].attrs["Hydrogen ionization transition temperature"][0]
+T_initial = sim["/HydroScheme"].attrs["Initial temperature"][0]
+T_minimal = sim["/HydroScheme"].attrs["Minimal temperature"][0]
+git = sim["Code"].attrs["Git Revision"]
+
+# Cosmological parameters
+H_0 = sim["/Cosmology"].attrs["H0 [internal units]"][0]
+gas_gamma = sim["/HydroScheme"].attrs["Adiabatic index"][0]
+
+unit_length_in_cgs = sim["/Units"].attrs["Unit length in cgs (U_L)"]
+unit_mass_in_cgs = sim["/Units"].attrs["Unit mass in cgs (U_M)"]
+unit_time_in_cgs = sim["/Units"].attrs["Unit time in cgs (U_t)"]
+
+unit_length_in_si = 0.01 * unit_length_in_cgs
+unit_mass_in_si = 0.001 * unit_mass_in_cgs
+unit_time_in_si = unit_time_in_cgs
+
+# Primoridal ean molecular weight as a function of temperature
+def mu(T, H_frac=H_mass_fraction, T_trans=H_transition_temp):
+    if T > T_trans:
+        return 4. / (8. - 5. * (1. - H_frac))
+    else:
+        return 4. / (1. + 3. * H_frac)
+    
+# Temperature of some primoridal gas with a given internal energy
+def T(u, H_frac=H_mass_fraction, T_trans=H_transition_temp):
+    T_over_mu = (gas_gamma - 1.) * u * mH_in_kg / k_in_J_K
+    ret = np.ones(np.size(u)) * T_trans
+
+    # Enough energy to be ionized?
+    mask_ionized = (T_over_mu > (T_trans+1) / mu(T_trans+1, H_frac, T_trans))
+    if np.sum(mask_ionized)  > 0:
+        ret[mask_ionized] = T_over_mu[mask_ionized] * mu(T_trans*10, H_frac, T_trans)
+
+    # Neutral gas?
+    mask_neutral = (T_over_mu < (T_trans-1) / mu((T_trans-1), H_frac, T_trans))
+    if np.sum(mask_neutral)  > 0:
+        ret[mask_neutral] = T_over_mu[mask_neutral] * mu(0, H_frac, T_trans)
+        
+    return ret
+
+
+z = np.zeros(n_snapshots)
+a = np.zeros(n_snapshots)
+T_mean = np.zeros(n_snapshots)
+T_std = np.zeros(n_snapshots)
+T_log_mean = np.zeros(n_snapshots)
+T_log_std = np.zeros(n_snapshots)
+T_median = np.zeros(n_snapshots)
+T_min = np.zeros(n_snapshots)
+T_max = np.zeros(n_snapshots)
+
+# Loop over all the snapshots
+for i in range(n_snapshots):
+    sim = h5py.File("snap_%04d.hdf5"%i, "r")
+
+    z[i] = sim["/Cosmology"].attrs["Redshift"][0]
+    a[i] = sim["/Cosmology"].attrs["Scale-factor"][0]
+
+    u = sim["/PartType0/InternalEnergies"][:]
+
+    # Compute the temperature
+    u *= (unit_length_in_si**2 / unit_time_in_si**2)
+    u /= a[i]**(3 * (gas_gamma - 1.))
+    Temp = T(u)
+
+    # Gather statistics
+    T_median[i] = np.median(Temp)
+    T_mean[i] = Temp.mean()
+    T_std[i] = Temp.std()
+    T_log_mean[i] = np.log10(Temp).mean()
+    T_log_std[i] = np.log10(Temp).std()
+    T_min[i] = Temp.min()
+    T_max[i] = Temp.max()
+
+# CMB evolution
+a_evol = np.logspace(-3, 0, 60)
+T_cmb = (1. / a_evol)**2 * 2.72
+
+# Plot the interesting quantities
+figure()
+subplot(111, xscale="log", yscale="log")
+
+fill_between(a, T_mean-T_std, T_mean+T_std, color='C0', alpha=0.1)
+plot(a, T_max, ls='-.', color='C0', lw=1., label="${\\rm max}~T$")
+plot(a, T_min, ls=':', color='C0', lw=1., label="${\\rm min}~T$")
+plot(a, T_mean, color='C0', label="${\\rm mean}~T$", lw=1.5)
+fill_between(a, 10**(T_log_mean-T_log_std), 10**(T_log_mean+T_log_std), color='C1', alpha=0.1)
+plot(a, 10**T_log_mean, color='C1', label="${\\rm mean}~{\\rm log} T$", lw=1.5)
+plot(a, T_median, color='C2', label="${\\rm median}~T$", lw=1.5)
+
+legend(loc="upper left", frameon=False, handlelength=1.5)
+
+# Expected lines
+plot([1e-10, 1e10], [H_transition_temp, H_transition_temp], 'k--', lw=0.5, alpha=0.7)
+text(2.5e-2, H_transition_temp*1.07, "$T_{\\rm HII\\rightarrow HI}$", va="bottom", alpha=0.7, fontsize=8)
+plot([1e-10, 1e10], [T_minimal, T_minimal], 'k--', lw=0.5, alpha=0.7)
+text(1e-2, T_minimal*0.8, "$T_{\\rm min}$", va="top", alpha=0.7, fontsize=8)
+plot(a_evol, T_cmb, 'k--', lw=0.5, alpha=0.7)
+text(a_evol[20], T_cmb[20]*0.55, "$(1+z)^2\\times T_{\\rm CMB,0}$", rotation=-34, alpha=0.7, fontsize=8, va="top", bbox=dict(facecolor='w', edgecolor='none', pad=1.0, alpha=0.9))
+
+
+redshift_ticks = np.array([0., 1., 2., 5., 10., 20., 50., 100.])
+redshift_labels = ["$0$", "$1$", "$2$", "$5$", "$10$", "$20$", "$50$", "$100$"]
+a_ticks = 1. / (redshift_ticks + 1.)
+
+xticks(a_ticks, redshift_labels)
+minorticks_off()
+
+xlabel("${\\rm Redshift}~z$", labelpad=0)
+ylabel("${\\rm Temperature}~T~[{\\rm K}]$", labelpad=0)
+xlim(9e-3, 1.1)
+ylim(20, 2.5e7)
+
+savefig("Temperature_evolution.png", dpi=200)
+
diff --git a/examples/SmallCosmoVolume/SmallCosmoVolume_Snipshots/run.sh b/examples/SmallCosmoVolume/SmallCosmoVolume_Snipshots/run.sh
new file mode 100755
index 0000000000000000000000000000000000000000..b2585d70b7cd2b717af02f005d690d0e8a9f932e
--- /dev/null
+++ b/examples/SmallCosmoVolume/SmallCosmoVolume_Snipshots/run.sh
@@ -0,0 +1,14 @@
+#!/bin/bash
+
+ # Generate the initial conditions if they are not present.
+if [ ! -e small_cosmo_volume.hdf5 ]
+then
+    echo "Fetching initial conditions for the small cosmological volume example..."
+    ./getIC.sh
+fi
+
+# Run SWIFT
+../../swift --cosmology --hydro --self-gravity --threads=8 small_cosmo_volume.yml 2>&1 | tee output.log
+
+# Plot the temperature evolution
+python plotTempEvolution.py
diff --git a/examples/SmallCosmoVolume/SmallCosmoVolume_Snipshots/select_output.yml b/examples/SmallCosmoVolume/SmallCosmoVolume_Snipshots/select_output.yml
new file mode 100644
index 0000000000000000000000000000000000000000..5bfd8cf3169ec8533acffdbe00c14a1d17bb0079
--- /dev/null
+++ b/examples/SmallCosmoVolume/SmallCosmoVolume_Snipshots/select_output.yml
@@ -0,0 +1,49 @@
+# All fields are written by default so in some sense the repeated 'on's are
+# unnecessary, but they are included for completeness.
+Snapshot:
+  # Particle Type Gas
+  Coordinates_Gas: on  # Co-moving positions of the particles : a U_L  [ cm ]
+  Velocities_Gas: on  # Peculiar velocities of the stars. This is (a * dx/dt) where x is the co-moving positions of the particles : U_L U_t^-1  [ cm s^-1 ]
+  Masses_Gas: on  # Masses of the particles : U_M  [ g ]
+  SmoothingLengths_Gas: on  # Co-moving smoothing lengths (FWHM of the kernel) of the particles : a U_L  [ cm ]
+  Entropies_Gas: on  # Co-moving entropies per unit mass of the particles : U_M^-1.6667 U_L^4 U_t^-2  [ g^-1.6667 cm^4 s^-2 ]
+  ParticleIDs_Gas: on  # Unique IDs of the particles : [ - ] 
+  Densities_Gas: on  # Co-moving mass densities of the particles : a^-3 U_M U_L^-3  [ g cm^-3 ]
+  InternalEnergies_Gas: on  # Co-moving thermal energies per unit mass of the particles : a^-2 U_L^2 U_t^-2  [ cm^2 s^-2 ]
+  Pressures_Gas: on  # Co-moving pressures of the particles : a^-5 U_M U_L^-1 U_t^-2  [ g cm^-1 s^-2 ]
+  Potentials_Gas: on  # Co-moving gravitational potential at position of the particles : a^-1 U_L^2 U_t^-2  [ cm^2 s^-2 ]
+  Temperatures_Gas: on  # Temperature of the particles : U_T  [ K ]
+  VELOCIraptorGroupIDs_Gas: on  # Group IDs of the particles in the VELOCIraptor catalogue : [ - ] 
+
+  # Particle Type DM
+  Coordinates_DM: on  # Co-moving position of the particles : a U_L  [ cm ]
+  Velocities_DM: on  # Peculiar velocities of the stars. This is a * dx/dt where x is the co-moving position of the particles. : U_L U_t^-1  [ cm s^-1 ]
+  Masses_DM: on  # Masses of the particles : U_M  [ g ]
+  ParticleIDs_DM: on  # Unique ID of the particles : [ - ] 
+  Softenings_DM: on  # Co-moving Plummer-equivalent softening lengths of the particles. : a U_L  [ cm ]
+  VELOCIraptorGroupIDs_DM: on  # Group IDs of the particles in the VELOCIraptor catalogue : [ - ] 
+
+
+Snipshot:
+  # Particle Type Gas
+  Coordinates_Gas: on  # Co-moving positions of the particles : a U_L  [ cm ]
+  Velocities_Gas: on  # Peculiar velocities of the stars. This is (a * dx/dt) where x is the co-moving positions of the particles : U_L U_t^-1  [ cm s^-1 ]
+  Masses_Gas: on  # Masses of the particles : U_M  [ g ]
+  SmoothingLengths_Gas: on  # Co-moving smoothing lengths (FWHM of the kernel) of the particles : a U_L  [ cm ]
+  Entropies_Gas: off  # Co-moving entropies per unit mass of the particles : U_M^-1.6667 U_L^4 U_t^-2  [ g^-1.6667 cm^4 s^-2 ]
+  ParticleIDs_Gas: on  # Unique IDs of the particles : [ - ] 
+  Densities_Gas: off  # Co-moving mass densities of the particles : a^-3 U_M U_L^-3  [ g cm^-3 ]
+  InternalEnergies_Gas: off  # Co-moving thermal energies per unit mass of the particles : a^-2 U_L^2 U_t^-2  [ cm^2 s^-2 ]
+  Pressures_Gas: off  # Co-moving pressures of the particles : a^-5 U_M U_L^-1 U_t^-2  [ g cm^-1 s^-2 ]
+  Potentials_Gas: off  # Co-moving gravitational potential at position of the particles : a^-1 U_L^2 U_t^-2  [ cm^2 s^-2 ]
+  Temperatures_Gas: off  # Temperature of the particles : U_T  [ K ]
+  VELOCIraptorGroupIDs_Gas: off  # Group IDs of the particles in the VELOCIraptor catalogue : [ - ] 
+
+
+  # Particle Type DM
+  Coordinates_DM: on  # Co-moving position of the particles : a U_L  [ cm ]
+  Velocities_DM: on  # Peculiar velocities of the stars. This is a * dx/dt where x is the co-moving position of the particles. : U_L U_t^-1  [ cm s^-1 ]
+  Masses_DM: off  # Masses of the particles : U_M  [ g ]
+  ParticleIDs_DM: on  # Unique ID of the particles : [ - ] 
+  Softenings_DM: off  # Co-moving Plummer-equivalent softening lengths of the particles. : a U_L  [ cm ]
+  VELOCIraptorGroupIDs_DM: off # Group IDs of the particles in the VELOCIraptor catalogue : [ - ] 
diff --git a/examples/SmallCosmoVolume/SmallCosmoVolume_Snipshots/small_cosmo_volume.yml b/examples/SmallCosmoVolume/SmallCosmoVolume_Snipshots/small_cosmo_volume.yml
new file mode 100644
index 0000000000000000000000000000000000000000..dc554e3dd8182b717803a902d1fb8b9f698d2f8e
--- /dev/null
+++ b/examples/SmallCosmoVolume/SmallCosmoVolume_Snipshots/small_cosmo_volume.yml
@@ -0,0 +1,79 @@
+# Define the system of units to use internally. 
+InternalUnitSystem:
+  UnitMass_in_cgs:     1.98841e43    # 10^10 M_sun
+  UnitLength_in_cgs:   3.08567758e24 # 1 Mpc
+  UnitVelocity_in_cgs: 1e5           # 1 km/s
+  UnitCurrent_in_cgs:  1             # Amperes
+  UnitTemp_in_cgs:     1             # Kelvin
+
+Cosmology:                      # WMAP9 cosmology
+  Omega_m:        0.276
+  Omega_lambda:   0.724
+  Omega_b:        0.0455
+  h:              0.703
+  a_begin:        0.019607843	# z_ini = 50.
+  a_end:          1.0		# z_end = 0.
+
+# Parameters governing the time integration
+TimeIntegration:
+  dt_min:     1e-6 
+  dt_max:     1e-2 
+
+# Parameters for the self-gravity scheme
+Gravity:
+  eta:          0.025         
+  theta:        0.5           
+  comoving_DM_softening:         0.0889     # 1/25th of the mean inter-particle separation: 88.9 kpc
+  max_physical_DM_softening:     0.0889     # 1/25th of the mean inter-particle separation: 88.9 kpc
+  comoving_baryon_softening:     0.0889     # 1/25th of the mean inter-particle separation: 88.9 kpc
+  max_physical_baryon_softening: 0.0889     # 1/25th of the mean inter-particle separation: 88.9 kpc
+  mesh_side_length:       64
+
+# Parameters of the hydro scheme
+SPH:
+  resolution_eta:      1.2348   # "48 Ngb" with the cubic spline kernel
+  h_min_ratio:         0.1
+  CFL_condition:       0.1
+  initial_temperature: 7075.    # (1 + z_ini)^2 * 2.72K
+  minimal_temperature: 100.
+
+# Parameters governing the snapshots
+Snapshots:
+  basename:            snap
+  delta_time:          1.05
+  scale_factor_first:  0.05
+  compression:         4
+  select_output_on:    1
+  select_output:       select_output.yml
+  output_list_on:      1
+  output_list:         output_list.txt
+  
+# Parameters governing the conserved quantities statistics
+Statistics:
+  delta_time:          1.02
+  scale_factor_first:  0.02
+  
+Scheduler:
+  max_top_level_cells: 8
+  cell_split_size:     50
+  
+# Parameters related to the initial conditions
+InitialConditions:
+  file_name:  small_cosmo_volume.hdf5
+  periodic:                    1
+  cleanup_h_factors:           1    
+  cleanup_velocity_factors:    1  
+  generate_gas_in_ics:         1      # Generate gas particles from the DM-only ICs
+  cleanup_smoothing_lengths:   1      # Since we generate gas, make use of the (expensive) cleaning-up procedure.
+
+# Parameters for the EAGLE "equation of state"
+EAGLEEntropyFloor:
+  Jeans_density_threshold_H_p_cm3: 0.1       # Physical density above which the EAGLE Jeans limiter entropy floor kicks in expressed in Hydrogen atoms per cm^3.
+  Jeans_over_density_threshold:    10.       # Overdensity above which the EAGLE Jeans limiter entropy floor can kick in.
+  Jeans_temperature_norm_K:        8000      # Temperature of the EAGLE Jeans limiter entropy floor at the density threshold expressed in Kelvin.
+  Jeans_gamma_effective:           1.3333333 # Slope the of the EAGLE Jeans limiter entropy floor
+  Cool_density_threshold_H_p_cm3: 1e-5       # Physical density above which the EAGLE Cool limiter entropy floor kicks in expressed in Hydrogen atoms per cm^3.
+  Cool_over_density_threshold:    10.        # Overdensity above which the EAGLE Cool limiter entropy floor can kick in.
+  Cool_temperature_norm_K:        8000       # Temperature of the EAGLE Cool limiter entropy floor at the density threshold expressed in Kelvin.
+  Cool_gamma_effective:           1.         # Slope the of the EAGLE Cool limiter entropy floor
+
diff --git a/examples/main.c b/examples/main.c
index 9a898f91206ad224dad83396418bfcf773cba5dc..90ec834dce338537196657b2e61f943104bf0c67 100644
--- a/examples/main.c
+++ b/examples/main.c
@@ -346,7 +346,7 @@ int main(int argc, char *argv[]) {
 
   /* Write output parameter file */
   if (myrank == 0 && output_parameters_filename != NULL) {
-    io_write_output_field_parameter(output_parameters_filename);
+    io_write_output_field_parameter(output_parameters_filename, with_cosmology);
     printf("End of run.\n");
     return 0;
   }
@@ -630,11 +630,20 @@ int main(int argc, char *argv[]) {
         parser_set_param(params, cmdps.param[k]);
     }
   }
+
 #ifdef WITH_MPI
   /* Broadcast the parameter file */
   MPI_Bcast(params, sizeof(struct swift_params), MPI_BYTE, 0, MPI_COMM_WORLD);
 #endif
 
+  /* Read the provided output selection file, if available. Best to
+   * do this after broadcasting the parameters as there may be code in this
+   * function that is repeated on each node based on the parameter file. */
+
+  struct output_options *output_options =
+      (struct output_options *)malloc(sizeof(struct output_options));
+  output_options_init(params, myrank, output_options);
+
   /* Temporary early aborts for modes not supported over MPI. */
 #ifdef WITH_MPI
   if (with_mpole_reconstruction && nr_nodes > 1)
@@ -1079,9 +1088,6 @@ int main(int argc, char *argv[]) {
     const int with_DM_background_particles =
         N_total[swift_type_dark_matter_background] > 0;
 
-    /* Verify that the fields to dump actually exist */
-    if (myrank == 0) io_check_output_fields(params, N_total);
-
     /* Initialize the space with these data. */
     if (myrank == 0) clocks_gettime(&tic);
     space_init(&s, params, &cosmo, dim, parts, gparts, sparts, bparts, Ngas,
@@ -1151,6 +1157,12 @@ int main(int argc, char *argv[]) {
     N_total[swift_type_black_hole] = s.nr_bparts;
 #endif
 
+    /* Verify that the fields to dump actually exist - this must be done after
+     * space_init so we know whether or not we have gas particles. */
+    if (myrank == 0)
+      io_check_output_fields(output_options->select_output, N_total,
+                             with_cosmology);
+
     /* Say a few nice things about the space we just created. */
     if (myrank == 0) {
       message("space dimensions are [ %.3f %.3f %.3f ].", s.dim[0], s.dim[1],
@@ -1219,14 +1231,15 @@ int main(int argc, char *argv[]) {
 
     /* Initialize the engine with the space and policies. */
     if (myrank == 0) clocks_gettime(&tic);
-    engine_init(
-        &e, &s, params, N_total[swift_type_gas], N_total[swift_type_count],
-        N_total[swift_type_stars], N_total[swift_type_black_hole],
-        N_total[swift_type_dark_matter_background], engine_policies, talking,
-        &reparttype, &us, &prog_const, &cosmo, &hydro_properties,
-        &entropy_floor, &gravity_properties, &stars_properties,
-        &black_holes_properties, &feedback_properties, &mesh, &potential,
-        &cooling_func, &starform, &chemistry, &fof_properties, &los_properties);
+    engine_init(&e, &s, params, output_options, N_total[swift_type_gas],
+                N_total[swift_type_count], N_total[swift_type_stars],
+                N_total[swift_type_black_hole],
+                N_total[swift_type_dark_matter_background], engine_policies,
+                talking, &reparttype, &us, &prog_const, &cosmo,
+                &hydro_properties, &entropy_floor, &gravity_properties,
+                &stars_properties, &black_holes_properties,
+                &feedback_properties, &mesh, &potential, &cooling_func,
+                &starform, &chemistry, &fof_properties, &los_properties);
     engine_config(/*restart=*/0, /*fof=*/0, &e, params, nr_nodes, myrank,
                   nr_threads, with_aff, talking, restart_file);
 
@@ -1559,6 +1572,7 @@ int main(int argc, char *argv[]) {
   if (with_feedback) feedback_clean(e.feedback_props);
   engine_clean(&e, /*fof=*/0, restart);
   free(params);
+  free(output_options);
 
 #ifdef WITH_MPI
   if ((res = MPI_Finalize()) != MPI_SUCCESS)
diff --git a/examples/main_fof.c b/examples/main_fof.c
index ac48a3bac84ca5d2ca6e043b737ccd1e2a6f4056..1ca3ab8402674d19745c5e646407919d025c36fc 100644
--- a/examples/main_fof.c
+++ b/examples/main_fof.c
@@ -257,7 +257,8 @@ int main(int argc, char *argv[]) {
 
   /* Write output parameter file */
   if (myrank == 0 && output_parameters_filename != NULL) {
-    io_write_output_field_parameter(output_parameters_filename);
+    io_write_output_field_parameter(output_parameters_filename,
+                                    /*with_cosmology=*/1);
     printf("End of run.\n");
     return 0;
   }
@@ -345,6 +346,13 @@ int main(int argc, char *argv[]) {
   MPI_Bcast(params, sizeof(struct swift_params), MPI_BYTE, 0, MPI_COMM_WORLD);
 #endif
 
+  /* "Read" the Select Output file - this should actually do nothing, but
+   * we need to mock the struct up for passing to `engine_init` */
+
+  struct output_options *output_options =
+      (struct output_options *)malloc(sizeof(struct output_options));
+  output_options_init(params, myrank, output_options);
+
   /* Check that we can write the snapshots by testing if the output
    * directory exists and is searchable and writable. */
   char basename[PARSER_MAX_LINE_SIZE];
@@ -580,16 +588,17 @@ int main(int argc, char *argv[]) {
 
   /* Initialize the engine with the space and policies. */
   if (myrank == 0) clocks_gettime(&tic);
-  engine_init(
-      &e, &s, params, N_total[swift_type_gas], N_total[swift_type_count],
-      N_total[swift_type_stars], N_total[swift_type_black_hole],
-      N_total[swift_type_dark_matter_background], engine_policies, talking,
-      &reparttype, &us, &prog_const, &cosmo,
-      /*hydro_properties=*/NULL, /*entropy_floor=*/NULL, &gravity_properties,
-      /*stars_properties=*/NULL, /*black_holes_properties=*/NULL,
-      /*feedback_properties=*/NULL, &mesh, /*potential=*/NULL,
-      /*cooling_func=*/NULL, /*starform=*/NULL, /*chemistry=*/NULL,
-      &fof_properties, /*los_properties=*/NULL);
+  engine_init(&e, &s, params, output_options, N_total[swift_type_gas],
+              N_total[swift_type_count], N_total[swift_type_stars],
+              N_total[swift_type_black_hole],
+              N_total[swift_type_dark_matter_background], engine_policies,
+              talking, &reparttype, &us, &prog_const, &cosmo,
+              /*hydro_properties=*/NULL, /*entropy_floor=*/NULL,
+              &gravity_properties,
+              /*stars_properties=*/NULL, /*black_holes_properties=*/NULL,
+              /*feedback_properties=*/NULL, &mesh, /*potential=*/NULL,
+              /*cooling_func=*/NULL, /*starform=*/NULL, /*chemistry=*/NULL,
+              &fof_properties, /*los_properties=*/NULL);
   engine_config(/*restart=*/0, /*fof=*/1, &e, params, nr_nodes, myrank,
                 nr_threads, with_aff, talking, NULL);
 
@@ -700,6 +709,7 @@ int main(int argc, char *argv[]) {
   pm_mesh_clean(&mesh);
   engine_clean(&e, /*fof=*/1, /*restart=*/0);
   free(params);
+  free(output_options);
 
   /* Say goodbye. */
   if (myrank == 0) message("done. Bye.");
diff --git a/examples/parameter_example.yml b/examples/parameter_example.yml
index fbd44cc56172a9cc1fb629afab9c7529e19a397e..1bc907c17f36c62c822ceee7b6163be012991f03 100644
--- a/examples/parameter_example.yml
+++ b/examples/parameter_example.yml
@@ -142,6 +142,8 @@ Snapshots:
   UnitTemp_in_cgs:     1  # (Optional) Unit system for the outputs (Kelvin)
   output_list_on:      0  # (Optional) Enable the output list
   output_list:         snaplist.txt # (Optional) File containing the output times (see documentation in "Parameter File" section)
+  select_output_on:    0  # (Optional) Enable the output selection behaviour
+  select_output:       selectoutput.yml # (Optional) File containing information to select outputs with (see documentation in the "Output Selection" section)
 
 # Parameters governing the logger snapshot system
 Logger:
diff --git a/src/Makefile.am b/src/Makefile.am
index 09b3ae7a073a003d49265ac0573c3e452af72716..17794b3d1e204bfdf1ec296918ae08e2e33681e6 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -50,7 +50,7 @@ include_HEADERS = space.h runner.h queue.h task.h lock.h cell.h part.h const.h \
     dump.h logger.h active.h timeline.h xmf.h gravity_properties.h gravity_derivatives.h \
     gravity_softened_derivatives.h vector_power.h collectgroup.h hydro_space.h sort_part.h \
     chemistry.h chemistry_io.h chemistry_struct.h cosmology.h restart.h space_getsid.h utilities.h \
-    mesh_gravity.h cbrt.h exp10.h velociraptor_interface.h swift_velociraptor_part.h outputlist.h \
+    mesh_gravity.h cbrt.h exp10.h velociraptor_interface.h swift_velociraptor_part.h output_list.h \
     logger_io.h tracers_io.h tracers.h tracers_struct.h star_formation_io.h fof.h fof_struct.h fof_io.h \
     multipole.h multipole_struct.h sincos.h \
     star_formation_struct.h star_formation.h star_formation_iact.h \
@@ -108,8 +108,8 @@ AM_SOURCES = space.c runner_main.c runner_doiact_hydro.c runner_doiact_limiter.c
     part_type.c xmf.c gravity_properties.c gravity.c \
     collectgroup.c hydro_space.c equation_of_state.c \
     chemistry.c cosmology.c restart.c mesh_gravity.c velociraptor_interface.c \
-    outputlist.c velociraptor_dummy.c logger_io.c memuse.c mpiuse.c memuse_rnodes.c fof.c \
-    hashmap.c pressure_floor.c space_unique_id.c line_of_sight.c \
+    output_list.c velociraptor_dummy.c logger_io.c memuse.c mpiuse.c memuse_rnodes.c fof.c \
+    hashmap.c pressure_floor.c space_unique_id.c output_options.c line_of_sight.c \
     $(QLA_COOLING_SOURCES) \
     $(EAGLE_COOLING_SOURCES) $(EAGLE_FEEDBACK_SOURCES) \
     $(GRACKLE_COOLING_SOURCES) $(GEAR_FEEDBACK_SOURCES)
@@ -122,7 +122,7 @@ nobase_noinst_HEADERS = align.h approx_math.h atomic.h barrier.h cycle.h error.h
 		 runner_doiact_functions_limiter.h runner_doiact_limiter.h units.h intrinsics.h minmax.h \
                  kick.h timestep.h drift.h adiabatic_index.h io_properties.h dimension.h part_type.h periodic.h memswap.h \
                  timestep_limiter.h timestep_limiter_iact.h timestep_sync.h timestep_sync_part.h timestep_limiter_struct.h \
-                 dump.h logger.h sign.h logger_io.h hashmap.h gravity.h gravity_io.h gravity_cache.h \
+                 dump.h logger.h sign.h logger_io.h hashmap.h gravity.h gravity_io.h gravity_cache.h output_options.h \
 		 gravity/Default/gravity.h gravity/Default/gravity_iact.h gravity/Default/gravity_io.h \
 		 gravity/Default/gravity_debug.h gravity/Default/gravity_part.h  \
 		 gravity/Potential/gravity.h gravity/Potential/gravity_iact.h gravity/Potential/gravity_io.h \
diff --git a/src/common_io.c b/src/common_io.c
index 8846c8847b609cd84c5594cae38aea3edfcbf07d..c7d78b0086675fbd7bbf5c242efe6a98fe0f13ab 100644
--- a/src/common_io.c
+++ b/src/common_io.c
@@ -2271,115 +2271,148 @@ void io_collect_gparts_background_to_write(
 /**
  * @brief Verify the io parameter file
  *
- * @param params The #swift_params
+ * @param params The #swift_params instance corresponding to the select_output
+ *               file.
  * @param N_total The total number of each particle type.
+ * @param with_cosmolgy Ran with cosmology?
  */
-void io_check_output_fields(const struct swift_params* params,
-                            const long long N_total[swift_type_count]) {
+void io_check_output_fields(struct swift_params* params,
+                            const long long N_total[swift_type_count],
+                            int with_cosmology) {
 
-  /* Loop over all particle types to check the fields */
-  for (int ptype = 0; ptype < swift_type_count; ptype++) {
-
-    int num_fields = 0;
-    struct io_props list[100];
-
-    /* Don't do anything if no particle of this kind */
-    if (N_total[ptype] == 0) continue;
-
-    /* Gather particle fields from the particle structures */
-    switch (ptype) {
-
-      case swift_type_gas:
-        hydro_write_particles(NULL, NULL, list, &num_fields);
-        num_fields += chemistry_write_particles(NULL, list + num_fields);
-        num_fields +=
-            cooling_write_particles(NULL, NULL, list + num_fields, NULL);
-        num_fields += tracers_write_particles(NULL, NULL, list + num_fields,
-                                              /*with_cosmology=*/1);
-        num_fields +=
-            star_formation_write_particles(NULL, NULL, list + num_fields);
-        num_fields += fof_write_parts(NULL, NULL, list + num_fields);
-        num_fields += velociraptor_write_parts(NULL, NULL, list + num_fields);
-        break;
-
-      case swift_type_dark_matter:
-        darkmatter_write_particles(NULL, list, &num_fields);
-        num_fields += fof_write_gparts(NULL, list + num_fields);
-        num_fields += velociraptor_write_gparts(NULL, list + num_fields);
-        break;
+  /* Loop over each section */
+  for (int section_id = 0; section_id < params->sectionCount; section_id++) {
+    char section_name[FIELD_BUFFER_SIZE];
+    sprintf(section_name, "%s", params->section[section_id].name);
 
-      case swift_type_dark_matter_background:
-        darkmatter_write_particles(NULL, list, &num_fields);
-        num_fields += fof_write_gparts(NULL, list + num_fields);
-        num_fields += velociraptor_write_gparts(NULL, list + num_fields);
-        break;
-
-      case swift_type_stars:
-        stars_write_particles(NULL, list, &num_fields, /*with_cosmology=*/1);
-        num_fields += chemistry_write_sparticles(NULL, list + num_fields);
-        num_fields += tracers_write_sparticles(NULL, list + num_fields,
-                                               /*with_cosmology=*/1);
-        num_fields += star_formation_write_sparticles(NULL, list + num_fields);
-        num_fields += fof_write_sparts(NULL, list + num_fields);
-        num_fields += velociraptor_write_sparts(NULL, list + num_fields);
-        break;
-
-      case swift_type_black_hole:
-        black_holes_write_particles(NULL, list, &num_fields,
-                                    /*with_cosmology=*/1);
-        num_fields += chemistry_write_bparticles(NULL, list + num_fields);
-        num_fields += fof_write_bparts(NULL, list + num_fields);
-        num_fields += velociraptor_write_bparts(NULL, list + num_fields);
-        break;
-
-      default:
-        error("Particle Type %d not yet supported. Aborting", ptype);
-    }
-
-    /* loop over each parameter */
+    /* Loop over each parameter */
     for (int param_id = 0; param_id < params->paramCount; param_id++) {
+
       const char* param_name = params->data[param_id].name;
 
-      char section_name[PARSER_MAX_LINE_SIZE];
+      char comparison_section_name[FIELD_BUFFER_SIZE];
 
       /* Skip if wrong section */
-      sprintf(section_name, "SelectOutput:");
-      if (strstr(param_name, section_name) == NULL) continue;
+      sprintf(comparison_section_name, "%s", "SelectOutput:");
+      if (strstr(param_name, comparison_section_name) != NULL) {
+        error(
+            "Output selection files no longer require the use of top level "
+            "SelectOutput; see the documentation for changes.");
+        continue;
+      }
 
-      /* Skip if wrong particle type */
-      sprintf(section_name, "_%s", part_type_names[ptype]);
-      if (strstr(param_name, section_name) == NULL) continue;
+      /* Skip if top-level section */
+      sprintf(comparison_section_name, "%s", section_name);
+      if (strstr(param_name, comparison_section_name) == NULL) continue;
 
+      /* Loop over all particle types to check the fields */
       int found = 0;
+      for (int ptype = 0; ptype < swift_type_count; ptype++) {
+
+        /* Skip if wrong particle type */
+        sprintf(comparison_section_name, "_%s", part_type_names[ptype]);
+        if (strstr(param_name, section_name) == NULL) continue;
+
+        int num_fields = 0;
+        struct io_props list[100];
+
+        /* Don't do anything if no particle of this kind */
+        if (N_total[ptype] == 0) continue;
+
+        /* Gather particle fields from the particle structures */
+        switch (ptype) {
+
+          case swift_type_gas:
+            hydro_write_particles(NULL, NULL, list, &num_fields);
+            num_fields += chemistry_write_particles(NULL, list + num_fields);
+            num_fields +=
+                cooling_write_particles(NULL, NULL, list + num_fields, NULL);
+            num_fields += tracers_write_particles(NULL, NULL, list + num_fields,
+                                                  with_cosmology);
+            num_fields +=
+                star_formation_write_particles(NULL, NULL, list + num_fields);
+            num_fields += fof_write_parts(NULL, NULL, list + num_fields);
+            num_fields +=
+                velociraptor_write_parts(NULL, NULL, list + num_fields);
+            break;
+
+          case swift_type_dark_matter:
+            darkmatter_write_particles(NULL, list, &num_fields);
+            num_fields += fof_write_gparts(NULL, list + num_fields);
+            num_fields += velociraptor_write_gparts(NULL, list + num_fields);
+            break;
+
+          case swift_type_dark_matter_background:
+            darkmatter_write_particles(NULL, list, &num_fields);
+            num_fields += fof_write_gparts(NULL, list + num_fields);
+            num_fields += velociraptor_write_gparts(NULL, list + num_fields);
+            break;
+
+          case swift_type_stars:
+            stars_write_particles(NULL, list, &num_fields, with_cosmology);
+            num_fields += chemistry_write_sparticles(NULL, list + num_fields);
+            num_fields += tracers_write_sparticles(NULL, list + num_fields,
+                                                   with_cosmology);
+            num_fields +=
+                star_formation_write_sparticles(NULL, list + num_fields);
+            num_fields += fof_write_sparts(NULL, list + num_fields);
+            num_fields += velociraptor_write_sparts(NULL, list + num_fields);
+            break;
+
+          case swift_type_black_hole:
+            black_holes_write_particles(NULL, list, &num_fields,
+                                        with_cosmology);
+            num_fields += chemistry_write_bparticles(NULL, list + num_fields);
+            num_fields += fof_write_bparts(NULL, list + num_fields);
+            num_fields += velociraptor_write_bparts(NULL, list + num_fields);
+            break;
+
+          default:
+            error("Particle Type %d not yet supported. Aborting", ptype);
+        }
 
-      /* loop over each possible output field */
-      for (int field_id = 0; field_id < num_fields; field_id++) {
-        char field_name[PARSER_MAX_LINE_SIZE];
-        sprintf(field_name, "SelectOutput:%.*s_%s", FIELD_BUFFER_SIZE,
-                list[field_id].name, part_type_names[ptype]);
-
-        if (strcmp(param_name, field_name) == 0) {
-          found = 1;
-          /* check if correct input */
-          int retParam = 0;
-          char str[PARSER_MAX_LINE_SIZE];
-          sscanf(params->data[param_id].value, "%d%s", &retParam, str);
-
-          /* Check that we have a 0 or 1 */
-          if (retParam != 0 && retParam != 1)
-            message(
-                "WARNING: Unexpected input for %s. Received %i but expect 0 or "
-                "1. ",
-                field_name, retParam);
-
-          /* Found it, so move to the next one. */
-          break;
+        /* For this particle type, loop over each possible output field */
+        for (int field_id = 0; field_id < num_fields; field_id++) {
+          char field_name[PARSER_MAX_LINE_SIZE];
+          /* Note that section_name includes a : */
+          sprintf(field_name, "%s%.*s_%s", section_name, FIELD_BUFFER_SIZE,
+                  list[field_id].name, part_type_names[ptype]);
+
+          if (strcmp(param_name, field_name) == 0) {
+            found = 1;
+
+            /* Perform a correctness check on the _value_ of that
+             * parameter */
+            char field_value[FIELD_BUFFER_SIZE];
+            parser_get_param_string(params, field_name, &field_value[0]);
+
+            int value_is_valid = 0;
+
+            for (int allowed_value_index = 0;
+                 allowed_value_index < compression_level_count;
+                 allowed_value_index++) {
+              if (strcmp(field_value,
+                         compression_level_names[allowed_value_index]) == 0) {
+                value_is_valid = 1;
+                break;
+              }
+            }
+
+            if (value_is_valid) {
+              /* Found value and it is correct, so move to the next one. */
+              break;
+            } else {
+              error("Choice of output selection parameter %s:%s is invalid.",
+                    field_name, field_value);
+            }
+          }
         }
       }
       if (!found)
         message(
-            "WARNING: Trying to dump particle field '%s' (read from '%s') that "
-            "does not exist.",
+            "WARNING: Trying to change behaviour of field '%s' (read from "
+            "'%s') that does not exist. This may because you are not running "
+            "with all of the physics that you compiled the code with.",
             param_name, params->fileName);
     }
   }
@@ -2389,8 +2422,9 @@ void io_check_output_fields(const struct swift_params* params,
  * @brief Write the output field parameters file
  *
  * @param filename The file to write.
+ * @param with_cosmology Use cosmological name variant?
  */
-void io_write_output_field_parameter(const char* filename) {
+void io_write_output_field_parameter(const char* filename, int with_cosmology) {
 
   FILE* file = fopen(filename, "w");
   if (file == NULL) error("Error opening file '%s'", filename);
@@ -2400,7 +2434,7 @@ void io_write_output_field_parameter(const char* filename) {
   units_init_cgs(&snapshot_units);
 
   /* Loop over all particle types */
-  fprintf(file, "SelectOutput:\n");
+  fprintf(file, "Default:\n");
   for (int ptype = 0; ptype < swift_type_count; ptype++) {
 
     int num_fields = 0;
@@ -2415,7 +2449,7 @@ void io_write_output_field_parameter(const char* filename) {
         num_fields +=
             cooling_write_particles(NULL, NULL, list + num_fields, NULL);
         num_fields += tracers_write_particles(NULL, NULL, list + num_fields,
-                                              /*with_cosmology=*/1);
+                                              with_cosmology);
         num_fields +=
             star_formation_write_particles(NULL, NULL, list + num_fields);
         num_fields += fof_write_parts(NULL, NULL, list + num_fields);
@@ -2435,18 +2469,17 @@ void io_write_output_field_parameter(const char* filename) {
         break;
 
       case swift_type_stars:
-        stars_write_particles(NULL, list, &num_fields, /*with_cosmology=*/1);
+        stars_write_particles(NULL, list, &num_fields, with_cosmology);
         num_fields += chemistry_write_sparticles(NULL, list + num_fields);
-        num_fields += tracers_write_sparticles(NULL, list + num_fields,
-                                               /*with_cosmology=*/1);
+        num_fields +=
+            tracers_write_sparticles(NULL, list + num_fields, with_cosmology);
         num_fields += star_formation_write_sparticles(NULL, list + num_fields);
         num_fields += fof_write_sparts(NULL, list + num_fields);
         num_fields += velociraptor_write_sparts(NULL, list + num_fields);
         break;
 
       case swift_type_black_hole:
-        black_holes_write_particles(NULL, list, &num_fields,
-                                    /*with_cosmology=*/1);
+        black_holes_write_particles(NULL, list, &num_fields, with_cosmology);
         num_fields += chemistry_write_bparticles(NULL, list + num_fields);
         num_fields += fof_write_bparts(NULL, list + num_fields);
         num_fields += velociraptor_write_bparts(NULL, list + num_fields);
@@ -2464,14 +2497,25 @@ void io_write_output_field_parameter(const char* filename) {
     /* Write all the fields of this particle type */
     for (int i = 0; i < num_fields; ++i) {
 
-      char buffer[FIELD_BUFFER_SIZE] = {0};
-      units_cgs_conversion_string(buffer, &snapshot_units, list[i].units,
+      char unit_buffer[FIELD_BUFFER_SIZE] = {0};
+      units_cgs_conversion_string(unit_buffer, &snapshot_units, list[i].units,
                                   list[i].scale_factor_exponent);
 
-      fprintf(file,
-              "  %s_%s: %*d \t # %s. ::: Conversion to physical CGS: %s\n",
-              list[i].name, part_type_names[ptype],
-              (int)(28 - strlen(list[i].name)), 1, list[i].description, buffer);
+      /* Need to buffer with a maximal size - otherwise we can't read in again
+       * because comments are too long */
+      char comment_write_buffer[PARSER_MAX_LINE_SIZE / 2];
+
+      sprintf(comment_write_buffer, "%.*s", PARSER_MAX_LINE_SIZE / 2 - 1,
+              list[i].description);
+
+      /* If our string is too long, replace the last few characters (before
+       * \0) with ... for 'fancy printing' */
+      if (strlen(comment_write_buffer) > PARSER_MAX_LINE_SIZE / 2 - 3) {
+        strcpy(&comment_write_buffer[PARSER_MAX_LINE_SIZE / 2 - 4], "...");
+      }
+
+      fprintf(file, "  %s_%s: %s  # %s : %s\n", list[i].name,
+              part_type_names[ptype], "on", comment_write_buffer, unit_buffer);
     }
 
     fprintf(file, "\n");
diff --git a/src/common_io.h b/src/common_io.h
index 2a6392dc69fe03f0fba5511ca6936e5d1e0ce8e9..9d9512b18688ba12f0b2a19dd94e8b21ea59a6a8 100644
--- a/src/common_io.h
+++ b/src/common_io.h
@@ -167,10 +167,10 @@ void io_duplicate_black_holes_gparts(struct threadpool* tp,
                                      struct gpart* const gparts, size_t Nstars,
                                      size_t Ndm);
 
-void io_check_output_fields(const struct swift_params* params,
-                            const long long N_total[3]);
+void io_check_output_fields(struct swift_params* params,
+                            const long long N_total[3], int with_cosmology);
 
-void io_write_output_field_parameter(const char* filename);
+void io_write_output_field_parameter(const char* filename, int with_cosmology);
 
 void io_make_snapshot_subdir(const char* dirname);
 
diff --git a/src/distributed_io.c b/src/distributed_io.c
index 355fdb05ac33b40ffdb0c2b9630b37a2e4b136c5..a2e6e5f0f3929703e2fe83146b2fa87d2c224d22 100644
--- a/src/distributed_io.c
+++ b/src/distributed_io.c
@@ -50,6 +50,8 @@
 #include "hydro_properties.h"
 #include "io_properties.h"
 #include "memuse.h"
+#include "output_list.h"
+#include "output_options.h"
 #include "part.h"
 #include "part_type.h"
 #include "star_formation_io.h"
@@ -246,7 +248,8 @@ void write_output_distributed(struct engine* e,
   const struct gpart* gparts = e->s->gparts;
   const struct spart* sparts = e->s->sparts;
   const struct bpart* bparts = e->s->bparts;
-  struct swift_params* params = e->parameter_file;
+  struct output_options* output_options = e->output_options;
+  struct output_list* output_list = e->output_list_snapshots;
   const int with_cosmology = e->policy & engine_policy_cosmology;
   const int with_cooling = e->policy & engine_policy_cooling;
   const int with_temperature = e->policy & engine_policy_temperature;
@@ -353,6 +356,16 @@ void write_output_distributed(struct engine* e,
                          e->s->dim[1] * factor_length,
                          e->s->dim[2] * factor_length};
 
+  /* Determine if we are writing a reduced snapshot, and if so which
+   * output selection type to use */
+  char current_selection_name[FIELD_BUFFER_SIZE] =
+      select_output_header_default_name;
+  if (output_list) {
+    /* Users could have specified a different Select Output scheme for each
+     * snapshot. */
+    output_list_get_current_select_output(output_list, current_selection_name);
+  }
+
   /* Print the relevant information and print status */
   io_write_attribute(h_grp, "BoxSize", DOUBLE, dim, 3);
   io_write_attribute(h_grp, "Time", DOUBLE, &dblTime, 1);
@@ -392,6 +405,7 @@ void write_output_distributed(struct engine* e,
   io_write_attribute_i(h_grp, "NumFilesPerSnapshot", numFiles);
   io_write_attribute_i(h_grp, "ThisFile", mpi_rank);
   io_write_attribute_s(h_grp, "OutputType", "Snapshot");
+  io_write_attribute_s(h_grp, "SelectOutput", current_selection_name);
 
   /* Close header */
   H5Gclose(h_grp);
@@ -704,13 +718,12 @@ void write_output_distributed(struct engine* e,
     }
 
     /* Write everything that is not cancelled */
+
     for (int i = 0; i < num_fields; ++i) {
 
       /* Did the user cancel this field? */
-      char field[PARSER_MAX_LINE_SIZE];
-      sprintf(field, "SelectOutput:%.*s_%s", FIELD_BUFFER_SIZE, list[i].name,
-              part_type_names[ptype]);
-      int should_write = parser_get_opt_param_int(params, field, 1);
+      const int should_write = output_options_should_write_field(
+          output_options, current_selection_name, list[i].name, ptype);
 
       if (should_write)
         write_distributed_array(e, h_grp, fileName, partTypeGroupName, list[i],
diff --git a/src/engine.c b/src/engine.c
index 48aa501d165eaff1514ab2dcdfc2f346b35cc131..83e1dff8ee2f7e6e5ffc2d05a9440ef95ecb8ec4 100644
--- a/src/engine.c
+++ b/src/engine.c
@@ -77,7 +77,8 @@
 #include "memuse.h"
 #include "minmax.h"
 #include "mpiuse.h"
-#include "outputlist.h"
+#include "output_list.h"
+#include "output_options.h"
 #include "parallel_io.h"
 #include "part.h"
 #include "partition.h"
@@ -3818,6 +3819,7 @@ static void engine_dumper_init(struct engine *e) {
  * @param e The #engine.
  * @param s The #space in which this #runner will run.
  * @param params The parsed parameter file.
+ * @param output_options Output options for snapshots.
  * @param Ngas total number of gas particles in the simulation.
  * @param Ngparts total number of gravity particles in the simulation.
  * @param Nstars total number of star particles in the simulation.
@@ -3843,9 +3845,10 @@ static void engine_dumper_init(struct engine *e) {
  * @param fof_properties The #fof_props.
  */
 void engine_init(struct engine *e, struct space *s, struct swift_params *params,
-                 long long Ngas, long long Ngparts, long long Nstars,
-                 long long Nblackholes, long long Nbackground_gparts,
-                 int policy, int verbose, struct repartition *reparttype,
+                 struct output_options *output_options, long long Ngas,
+                 long long Ngparts, long long Nstars, long long Nblackholes,
+                 long long Nbackground_gparts, int policy, int verbose,
+                 struct repartition *reparttype,
                  const struct unit_system *internal_units,
                  const struct phys_const *physical_constants,
                  struct cosmology *cosmo, struct hydro_props *hydro,
@@ -3939,6 +3942,7 @@ void engine_init(struct engine *e, struct space *s, struct swift_params *params,
   e->chemistry = chemistry;
   e->fof_properties = fof_properties;
   e->parameter_file = params;
+  e->output_options = output_options;
   e->stf_this_timestep = 0;
   e->los_properties = los_properties;
 #ifdef WITH_MPI
@@ -4790,7 +4794,7 @@ void engine_print_policy(struct engine *e) {
  */
 void engine_compute_next_snapshot_time(struct engine *e) {
 
-  /* Do outputlist file case */
+  /* Do output_list file case */
   if (e->output_list_snapshots) {
     output_list_read_next_time(e->output_list_snapshots, e, "snapshots",
                                &e->ti_next_snapshot);
@@ -5339,6 +5343,8 @@ void engine_clean(struct engine *e, const int fof, const int restart) {
   output_list_clean(&e->output_list_stats);
   output_list_clean(&e->output_list_stf);
 
+  output_options_clean(e->output_options);
+
   swift_free("links", e->links);
 #if defined(WITH_LOGGER)
   if (e->policy & engine_policy_logger) {
@@ -5379,6 +5385,7 @@ void engine_clean(struct engine *e, const int fof, const int restart) {
      in engine_struct_restore() */
   if (restart) {
     free((void *)e->parameter_file);
+    free((void *)e->output_options);
     free((void *)e->external_potential);
     free((void *)e->black_holes_properties);
     free((void *)e->stars_properties);
@@ -5454,6 +5461,7 @@ void engine_struct_dump(struct engine *e, FILE *stream) {
 #endif
   los_struct_dump(e->los_properties, stream);
   parser_struct_dump(e->parameter_file, stream);
+  output_options_struct_dump(e->output_options, stream);
   if (e->output_list_snapshots)
     output_list_struct_dump(e->output_list_snapshots, stream);
   if (e->output_list_stats)
@@ -5595,6 +5603,11 @@ void engine_struct_restore(struct engine *e, FILE *stream) {
   parser_struct_restore(parameter_file, stream);
   e->parameter_file = parameter_file;
 
+  struct output_options *output_options =
+      (struct output_options *)malloc(sizeof(struct output_options));
+  output_options_struct_restore(output_options, stream);
+  e->output_options = output_options;
+
   if (e->output_list_snapshots) {
     struct output_list *output_list_snapshots =
         (struct output_list *)malloc(sizeof(struct output_list));
diff --git a/src/engine.h b/src/engine.h
index e1ace452ecec54ff28d1fa6b474ff064183c2897..957bf86e0748c601f4c84d7ca8ab818c56a9358b 100644
--- a/src/engine.h
+++ b/src/engine.h
@@ -38,6 +38,7 @@
 #include "collectgroup.h"
 #include "dump.h"
 #include "mesh_gravity.h"
+#include "output_options.h"
 #include "parser.h"
 #include "partition.h"
 #include "potential.h"
@@ -455,6 +456,9 @@ struct engine {
   /* The (parsed) parameter file */
   struct swift_params *parameter_file;
 
+  /* The output selection options */
+  struct output_options *output_options;
+
   /* Temporary struct to hold a group of deferable properties (in MPI mode
    * these are reduced together, but may not be required just yet). */
   struct collectgroup1 collect_group1;
@@ -537,9 +541,10 @@ void engine_collect_end_of_step(struct engine *e, int apply);
 void engine_dump_snapshot(struct engine *e);
 void engine_init_output_lists(struct engine *e, struct swift_params *params);
 void engine_init(struct engine *e, struct space *s, struct swift_params *params,
-                 long long Ngas, long long Ngparts, long long Nstars,
-                 long long Nblackholes, long long Nbackground_gparts,
-                 int policy, int verbose, struct repartition *reparttype,
+                 struct output_options *output_options, long long Ngas,
+                 long long Ngparts, long long Nstars, long long Nblackholes,
+                 long long Nbackground_gparts, int policy, int verbose,
+                 struct repartition *reparttype,
                  const struct unit_system *internal_units,
                  const struct phys_const *physical_constants,
                  struct cosmology *cosmo, struct hydro_props *hydro,
diff --git a/src/outputlist.c b/src/output_list.c
similarity index 55%
rename from src/outputlist.c
rename to src/output_list.c
index da841103d2ce35b849112d62ffd523f632cc457b..fdc4fa7835ee906adab2f3540985073bda9c53dd 100644
--- a/src/outputlist.c
+++ b/src/output_list.c
@@ -21,7 +21,7 @@
 #include "../config.h"
 
 /* This object's header. */
-#include "outputlist.h"
+#include "output_list.h"
 
 /* Local includes. */
 #include "cosmology.h"
@@ -36,12 +36,12 @@
 /**
  * @brief Read a file containing a list of time
  *
- * @param outputlist The #output_list to fill.
+ * @param output_list The #output_list to fill.
  * @param filename The file to read.
  * @param cosmo The #cosmology model.
  */
-void output_list_read_file(struct output_list *outputlist, const char *filename,
-                           struct cosmology *cosmo) {
+void output_list_read_file(struct output_list *output_list,
+                           const char *filename, struct cosmology *cosmo) {
 
   /* Open file */
   FILE *file = fopen(filename, "r");
@@ -53,16 +53,20 @@ void output_list_read_file(struct output_list *outputlist, const char *filename,
   size_t nber_line = 0;
   while (getline(&line, &len, file) != -1) nber_line++;
 
-  outputlist->size = nber_line - 1; /* Do not count header */
+  output_list->size = nber_line - 1; /* Do not count header */
 
   /* Return to start of file and initialize time array */
   fseek(file, 0, SEEK_SET);
-  outputlist->times = (double *)malloc(sizeof(double) * outputlist->size);
-  if (!outputlist->times)
+  output_list->times = (double *)malloc(sizeof(double) * output_list->size);
+  output_list->select_output_indices =
+      (int *)malloc(sizeof(int) * output_list->size);
+
+  if ((!output_list->times) || (!output_list->select_output_indices)) {
     error(
         "Unable to malloc output_list. "
         "Try reducing the number of lines in %s",
         filename);
+  }
 
   /* Read header */
   if (getline(&line, &len, file) == -1)
@@ -73,16 +77,29 @@ void output_list_read_file(struct output_list *outputlist, const char *filename,
 
   /* Find type of data in file */
   int type = -1;
+  output_list->select_output_on = 0;
+  output_list->select_output_number_of_names = 0;
 
   trim_trailing(line);
-  if (strcasecmp(line, "# Redshift") == 0)
+
+  if (strcasecmp(line, "# Redshift") == 0) {
     type = OUTPUT_LIST_REDSHIFT;
-  else if (strcasecmp(line, "# Time") == 0)
+  } else if (strcasecmp(line, "# Time") == 0) {
     type = OUTPUT_LIST_AGE;
-  else if (strcasecmp(line, "# Scale Factor") == 0)
+  } else if (strcasecmp(line, "# Scale Factor") == 0) {
     type = OUTPUT_LIST_SCALE_FACTOR;
-  else
+  } else if (strcasecmp(line, "# Redshift, Select Output") == 0) {
+    type = OUTPUT_LIST_REDSHIFT;
+    output_list->select_output_on = 1;
+  } else if (strcasecmp(line, "# Time, Select Output") == 0) {
+    type = OUTPUT_LIST_AGE;
+    output_list->select_output_on = 1;
+  } else if (strcasecmp(line, "# Scale Factor, Select Output") == 0) {
+    type = OUTPUT_LIST_SCALE_FACTOR;
+    output_list->select_output_on = 1;
+  } else {
     error("Unable to interpret the header (%s) in file '%s'", line, filename);
+  }
 
   if (!cosmo &&
       (type == OUTPUT_LIST_SCALE_FACTOR || type == OUTPUT_LIST_REDSHIFT))
@@ -93,14 +110,26 @@ void output_list_read_file(struct output_list *outputlist, const char *filename,
 
   /* Read file */
   size_t ind = 0;
+  int read_successfully = 0;
+  int found_select_output = 0;
+  char select_output_buffer[FIELD_BUFFER_SIZE] =
+      select_output_header_default_name;
   while (getline(&line, &len, file) != -1) {
-    double *time = &outputlist->times[ind];
-    /* Write data to outputlist */
-    if (sscanf(line, "%lf", time) != 1)
+    double *time = &output_list->times[ind];
+    /* Write data to output_list */
+    if (output_list->select_output_on) {
+      read_successfully =
+          sscanf(line, "%lf, %s", time, select_output_buffer) == 2;
+    } else {
+      read_successfully = sscanf(line, "%lf", time) == 1;
+    }
+
+    if (!read_successfully) {
       error(
-          "Tried parsing double but found '%s' with illegal double "
+          "Tried parsing output_list but found '%s' with illegal "
           "characters in file '%s'.",
           line, filename);
+    }
 
     /* Transform input into correct time (e.g. ages or scale factor) */
     if (type == OUTPUT_LIST_REDSHIFT) *time = 1. / (1. + *time);
@@ -108,6 +137,32 @@ void output_list_read_file(struct output_list *outputlist, const char *filename,
     if (cosmo && type == OUTPUT_LIST_AGE)
       *time = cosmology_get_scale_factor(cosmo, *time);
 
+    /* Search to find index for select output - select_output_index is the index
+     * in the select_output_names array that corresponds to this select output
+     * name. */
+    found_select_output = 0;
+    for (int select_output_index = 0;
+         select_output_index < output_list->select_output_number_of_names;
+         select_output_index++) {
+      if (!strcmp(select_output_buffer,
+                  output_list->select_output_names[select_output_index])) {
+        /* We already have this select output list string in the buffer! */
+        output_list->select_output_indices[ind] = select_output_index;
+        found_select_output = 1;
+      }
+    }
+    /* If we did not assign it above, we haven't encountered this name before
+     * and we need to create this name in the array */
+    if (!found_select_output) {
+      strcpy(
+          output_list
+              ->select_output_names[output_list->select_output_number_of_names],
+          select_output_buffer);
+      output_list->select_output_indices[ind] =
+          output_list->select_output_number_of_names;
+      output_list->select_output_number_of_names += 1;
+    }
+
     /* Update size */
     ind += 1;
   }
@@ -115,30 +170,30 @@ void output_list_read_file(struct output_list *outputlist, const char *filename,
   /* Cleanup */
   free(line);
 
-  if (ind != outputlist->size)
+  if (ind != output_list->size)
     error("Did not read the correct number of output times.");
 
   /* Check that the list is in monotonic order */
-  for (size_t i = 1; i < outputlist->size; ++i) {
+  for (size_t i = 1; i < output_list->size; ++i) {
 
     if ((type == OUTPUT_LIST_REDSHIFT) &&
-        (outputlist->times[i] <= outputlist->times[i - 1]))
+        (output_list->times[i] <= output_list->times[i - 1]))
       error("Output list not having monotonically decreasing redshifts.");
 
     if ((type == OUTPUT_LIST_AGE) &&
-        (outputlist->times[i] <= outputlist->times[i - 1]))
+        (output_list->times[i] <= output_list->times[i - 1]))
       error("Output list not having monotonically increasing ages.");
 
     if ((type == OUTPUT_LIST_SCALE_FACTOR) &&
-        (outputlist->times[i] <= outputlist->times[i - 1]))
+        (output_list->times[i] <= output_list->times[i - 1]))
       error("Output list not having monotonically increasing scale-factors.");
   }
 
   /* set current indice to 0 */
-  outputlist->cur_ind = 0;
+  output_list->cur_ind = 0;
 
   /* We check if this is true later */
-  outputlist->final_step_dump = 0;
+  output_list->final_step_dump = 0;
 
   fclose(file);
 }
@@ -162,7 +217,7 @@ void output_list_read_next_time(struct output_list *t, const struct engine *e,
   else
     time_end = e->time_end;
 
-  /* Find next snasphot above current time */
+  /* Find next snapshot above current time */
   double time = t->times[t->cur_ind];
   size_t ind = t->cur_ind;
   while (time <= time_end) {
@@ -183,7 +238,11 @@ void output_list_read_next_time(struct output_list *t, const struct engine *e,
     t->cur_ind = ind;
   }
 
-  /* Do we need to do a dump at the end of the last timestep? */
+  /* Do we need to do a dump at the end of the last timestep?
+   * Note that what this really does is given that a t=t_max, z=0,
+   * or a=1 is found in output_list.txt set the flag `final_step_dump`
+   * to 1 - this is not special behaviour that is controlled by a
+   * parameter file flag. */
   if (time == time_end) {
     t->final_step_dump = 1;
     if (e->verbose) {
@@ -198,8 +257,12 @@ void output_list_read_next_time(struct output_list *t, const struct engine *e,
   /* Deal with last statistics */
   if (*ti_next >= max_nr_timesteps || ind == t->size || time >= time_end) {
     *ti_next = -1;
-    if (e->verbose && t->final_step_dump != 1)
+    if (e->verbose && t->final_step_dump != 1) {
       message("No further output time for %s.", name);
+
+      /* Do not print anything about the output style (below) */
+      return;
+    }
   } else {
 
     /* Be nice, talk... */
@@ -214,6 +277,28 @@ void output_list_read_next_time(struct output_list *t, const struct engine *e,
         message("Next output time for %s set to t=%e.", name, next_time);
     }
   }
+
+  /* Finally, talk if we are a snapshot and we are using SelectOutput */
+  if (e->verbose) {
+    if (t->select_output_number_of_names > 1) {
+      char select_output_style[FIELD_BUFFER_SIZE];
+      output_list_get_current_select_output(t, select_output_style);
+      message("Next output style for %s set to %s.", name, select_output_style);
+    }
+  }
+}
+
+/**
+ * @brief Copys the string describing the current output name into the
+ *        buffer described by select_output_name.
+ *
+ * @param t The #output_list
+ * @param select_output_name updated to the current Select Output choice.
+ **/
+void output_list_get_current_select_output(struct output_list *t,
+                                           char *select_output_name) {
+  strcpy(select_output_name,
+         t->select_output_names[t->select_output_indices[t->cur_ind]]);
 }
 
 /**
@@ -223,8 +308,8 @@ void output_list_read_next_time(struct output_list *t, const struct engine *e,
  * @param e The #engine
  * @param name The name of the section in params
  * @param delta_time updated to the initial delta time
- * @param time_first updated to the time of first output (scale factor or cosmic
- * time)
+ * @param time_first updated to the time of first output (scale factor or
+ * cosmic time)
  */
 void output_list_init(struct output_list **list, const struct engine *e,
                       const char *name, double *delta_time,
@@ -238,13 +323,14 @@ void output_list_init(struct output_list **list, const struct engine *e,
   /* Read output on/off */
   char param_name[PARSER_MAX_LINE_SIZE];
   sprintf(param_name, "%s:output_list_on", name);
-  int outputlist_on = parser_get_opt_param_int(params, param_name, 0);
+  int output_list_on = parser_get_opt_param_int(params, param_name, 0);
 
-  /* Check if read outputlist */
-  if (!outputlist_on) return;
+  /* Check if read output_list */
+  if (!output_list_on) return;
 
-  /* Read outputlist for snapshots */
+  /* Read output_list for snapshots */
   *list = (struct output_list *)malloc(sizeof(struct output_list));
+  (*list)->output_list_on = output_list_on;
 
   /* Read filename */
   char filename[PARSER_MAX_LINE_SIZE];
@@ -270,26 +356,39 @@ void output_list_init(struct output_list **list, const struct engine *e,
 /**
  * @brief Print an #output_list
  */
-void output_list_print(const struct output_list *outputlist) {
+void output_list_print(const struct output_list *output_list) {
 
-  printf("/*\t Time Array\t */\n");
-  printf("Number of Line: %zu\n", outputlist->size);
-  for (size_t ind = 0; ind < outputlist->size; ind++) {
-    if (ind == outputlist->cur_ind)
-      printf("\t%lf <-- Current\n", outputlist->times[ind]);
+  printf("/*\t Total number of Select Output options: %d \t */\n",
+         output_list->select_output_number_of_names);
+  printf("/*\t Select Output options found in output list: \t */\n");
+  for (int i = 0; i < output_list->select_output_number_of_names; i++) {
+    printf("\t Index: %d, Name: %s\n", i, output_list->select_output_names[i]);
+  }
+  printf("\n/*\t Time Array, Select Output \t */\n");
+  printf("Number of Lines: %zu\n", output_list->size);
+  for (size_t ind = 0; ind < output_list->size; ind++) {
+    if (ind == output_list->cur_ind)
+      printf(
+          "\t %lf, %s <-- Current\n", output_list->times[ind],
+          output_list
+              ->select_output_names[output_list->select_output_indices[ind]]);
     else
-      printf("\t%lf\n", outputlist->times[ind]);
+      printf(
+          "\t %lf, %s\n", output_list->times[ind],
+          output_list
+              ->select_output_names[output_list->select_output_indices[ind]]);
   }
 }
 
 /**
  * @brief Clean an #output_list
  */
-void output_list_clean(struct output_list **outputlist) {
-  if (*outputlist) {
-    free((*outputlist)->times);
-    free(*outputlist);
-    *outputlist = NULL;
+void output_list_clean(struct output_list **output_list) {
+  if (*output_list) {
+    free((*output_list)->times);
+    free((*output_list)->select_output_indices);
+    free(*output_list);
+    *output_list = NULL;
   }
 }
 
@@ -302,6 +401,9 @@ void output_list_struct_dump(struct output_list *list, FILE *stream) {
 
   restart_write_blocks(list->times, list->size, sizeof(double), stream,
                        "output_list", "times");
+
+  restart_write_blocks(list->select_output_indices, list->size, sizeof(int),
+                       stream, "output_list", "select_output_indices");
 }
 
 /**
@@ -314,4 +416,8 @@ void output_list_struct_restore(struct output_list *list, FILE *stream) {
   list->times = (double *)malloc(sizeof(double) * list->size);
   restart_read_blocks(list->times, list->size, sizeof(double), stream, NULL,
                       "times");
+
+  list->select_output_indices = (int *)malloc(sizeof(int) * list->size);
+  restart_read_blocks(list->select_output_indices, list->size, sizeof(int),
+                      stream, NULL, "select_output_indices");
 }
diff --git a/src/outputlist.h b/src/output_list.h
similarity index 65%
rename from src/outputlist.h
rename to src/output_list.h
index 8ca6bc6fecba772d3bae41bbb63358e6f33eb0a5..7f8cfc14a5a651896a25bbde80ac4a0ecd1a46cb 100644
--- a/src/outputlist.h
+++ b/src/output_list.h
@@ -23,8 +23,11 @@
 #include "../config.h"
 
 /* Local includes */
+#include "common_io.h"
 #include "cosmology.h"
 
+#define OUTPUT_LIST_MAX_NUM_OF_SELECT_OUTPUT_STYLES 8
+
 struct engine;
 
 /**
@@ -41,27 +44,47 @@ enum output_list_type {
  */
 struct output_list {
 
+  /* Select output names. */
+  char select_output_names[OUTPUT_LIST_MAX_NUM_OF_SELECT_OUTPUT_STYLES]
+                          [FIELD_BUFFER_SIZE];
+
   /* Time array */
   double *times;
 
-  /* Size of the time array */
+  /* Select output indices - each index corresponds to a string
+   * in select_output. Chosen to be this instead of an array of
+   * pointers because of restarts. */
+  int *select_output_indices;
+
+  /* Total number of currently used select output names */
+  int select_output_number_of_names;
+
+  /* Size of the time array (i.e. number of outputs) */
   size_t size;
 
   /* Current index */
   size_t cur_ind;
 
+  /* Was the Select Output option used? */
+  int select_output_on;
+
+  /* Is this output list activated? */
+  int output_list_on;
+
   /* Dump on final timestep? */
   int final_step_dump;
 };
 
-void output_list_read_file(struct output_list *outputlist, const char *filename,
-                           struct cosmology *cosmo);
+void output_list_read_file(struct output_list *output_list,
+                           const char *filename, struct cosmology *cosmo);
 void output_list_read_next_time(struct output_list *t, const struct engine *e,
                                 const char *name, integertime_t *ti_next);
+void output_list_get_current_select_output(struct output_list *t,
+                                           char *select_output_name);
 void output_list_init(struct output_list **list, const struct engine *e,
                       const char *name, double *delta_time, double *time_first);
-void output_list_print(const struct output_list *outputlist);
-void output_list_clean(struct output_list **outputlist);
+void output_list_print(const struct output_list *output_list);
+void output_list_clean(struct output_list **output_list);
 void output_list_struct_dump(struct output_list *list, FILE *stream);
 void output_list_struct_restore(struct output_list *list, FILE *stream);
 int output_list_check_duplicates(const struct output_list *list_a,
diff --git a/src/output_options.c b/src/output_options.c
new file mode 100644
index 0000000000000000000000000000000000000000..23938223767213c3e3ff37c92c991cf91072c2dc
--- /dev/null
+++ b/src/output_options.c
@@ -0,0 +1,156 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2020 Josh Borrow (josh.borrow@durham.ac.uk)
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published
+ * by the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+#include "../config.h"
+
+/* Some standard headers. */
+#include <stdlib.h>
+
+/* MPI headers. */
+#ifdef WITH_MPI
+#include <mpi.h>
+#endif
+
+/* This object's header. */
+#include "output_options.h"
+
+/* Local headers. */
+#include "parser.h"
+#include "part_type.h"
+#include "swift.h"
+
+/* Compression level names. */
+const char* compression_level_names[compression_level_count] = {
+    "off", "on", "low", "med", "high"};
+
+/**
+ * @brief Initialise the output options struct with the information read
+ *        from file. Only rank 0 reads from file; this data is then broadcast
+ *        to all nodes.
+ *
+ * @param parameter_file the pre-parsed parameter file.
+ * @param mpi_rank the MPI rank of this node.
+ * @param output_options the empty output_options struct (return).
+ **/
+void output_options_init(struct swift_params* parameter_file, int mpi_rank,
+                         struct output_options* output_options) {
+
+  /* Load select_output */
+  struct swift_params* select_output =
+      (struct swift_params*)malloc(sizeof(struct swift_params));
+  if (select_output == NULL)
+    error("Error allocating memory for select output options.");
+
+  if (mpi_rank == 0) {
+    const int select_output_on = parser_get_opt_param_int(
+        parameter_file, "Snapshots:select_output_on", 0);
+
+    if (select_output_on) {
+      char select_param_filename[PARSER_MAX_LINE_SIZE];
+      parser_get_param_string(parameter_file, "Snapshots:select_output",
+                              select_param_filename);
+      message("Reading select output parameters from file '%s'",
+              select_param_filename);
+      parser_read_file(select_param_filename, select_output);
+    }
+  }
+
+#ifdef WITH_MPI
+  MPI_Bcast(select_output, sizeof(struct swift_params), MPI_BYTE, 0,
+            MPI_COMM_WORLD);
+#endif
+
+  output_options->select_output = select_output;
+}
+
+/**
+ * @breif Destroys an output_options instance.
+ *
+ * @param output_options the output_options struct to free the contents of.
+ **/
+void output_options_clean(struct output_options* output_options) {
+  if (output_options) {
+    free(output_options->select_output);
+  }
+}
+
+/**
+ * @brief Dumps the output_options struct to restart file
+ *
+ * @param output_options pointer to output options struct
+ * @param stream bytestream to write to
+ **/
+void output_options_struct_dump(struct output_options* output_options,
+                                FILE* stream) {
+  parser_struct_dump(output_options->select_output, stream);
+}
+
+/**
+ * @brief Loads the output_options struct from the restart file
+ *
+ * @param output_options pointer to the output options struct
+ * @param stream bytestream to read from
+ **/
+void output_options_struct_restore(struct output_options* output_options,
+                                   FILE* stream) {
+  struct swift_params* select_output =
+      (struct swift_params*)malloc(sizeof(struct swift_params));
+  parser_struct_restore(select_output, stream);
+
+  output_options->select_output = select_output;
+}
+
+/**
+ * @brief Decides whether or not a given field should be written. Returns
+ *        a truthy value if yes, falsey if not.
+ *
+ * @param output_options pointer to the output options struct
+ * @param snapshot_type pointer to a char array containing the type of
+ *        output (i.e. top level section in the yaml).
+ * @param field_name pointer to a char array containing the name of the
+ *        relevant field.
+ * @param part_type integer particle type
+ *
+ * @return should_write integer determining whether this field should be
+ *         written
+ **/
+int output_options_should_write_field(struct output_options* output_options,
+                                      char* snapshot_type, char* field_name,
+                                      enum part_type part_type) {
+  /* Full name for the field path */
+  char field[PARSER_MAX_LINE_SIZE];
+  sprintf(field, "%.*s:%.*s_%s", FIELD_BUFFER_SIZE, snapshot_type,
+          FIELD_BUFFER_SIZE, field_name, part_type_names[part_type]);
+
+  char compression_level[FIELD_BUFFER_SIZE];
+  parser_get_opt_param_string(
+      output_options->select_output, field, compression_level,
+      compression_level_names[compression_level_default]);
+
+  int should_write = strcmp(compression_level_names[compression_do_not_write],
+                            compression_level);
+
+#ifdef SWIFT_DEBUG_CHECKS
+  message(
+      "Check for whether %s should be written returned %s from a provided "
+      "value of \"%s\"",
+      field, should_write ? "True" : "False", compression_level);
+#endif
+
+  return should_write;
+}
diff --git a/src/output_options.h b/src/output_options.h
new file mode 100644
index 0000000000000000000000000000000000000000..02ba151b1103fe09a50740c99a693ded0c4eec1c
--- /dev/null
+++ b/src/output_options.h
@@ -0,0 +1,76 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2020 Josh Borrow (josh.borrow@durham.ac.uk)
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published
+ * by the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+#ifndef SWIFT_OUTPUT_OPTIONS_H
+#define SWIFT_OUTPUT_OPTIONS_H
+
+#include "config.h"
+#include "parser.h"
+#include "part_type.h"
+
+/* Compression level names */
+enum compression_levels {
+  compression_do_not_write = 0,
+  compression_write_lossless,
+  compression_write_low_lossy,
+  compression_write_med_lossy,
+  compression_write_high_lossy,
+  /* Counter, always leave last */
+  compression_level_count,
+};
+
+/* Default value for SelectOutput */
+#define compression_level_default compression_write_lossless
+
+/* Default name for the SelectOutput header */
+#define select_output_header_default_name "Default"
+
+/**
+ * @brief Names of the compression levels, used in the select_output.yml
+ *        parameter file.
+ **/
+extern const char* compression_level_names[];
+
+/**
+ * @brief Output selection properties, including the parsed files.
+ **/
+struct output_options {
+
+  /*! Select output file, parsed */
+  struct swift_params* select_output;
+
+  /* Pass-through struct for now but may need more later. */
+};
+
+/* Create and destroy */
+void output_options_init(struct swift_params* parameter_file, int mpi_rank,
+                         struct output_options* output_options);
+void output_options_clean(struct output_options* output_options);
+
+/* Restart files */
+void output_options_struct_dump(struct output_options* output_options,
+                                FILE* stream);
+void output_options_struct_restore(struct output_options* output_options,
+                                   FILE* stream);
+
+/* Logic functions */
+int output_options_should_write_field(struct output_options* output_options,
+                                      char* snapshot_type, char* field_name,
+                                      enum part_type part_type);
+
+#endif
\ No newline at end of file
diff --git a/src/parallel_io.c b/src/parallel_io.c
index be820e2ecd3434f7d41635f8b0397ff5a1ef368c..994fe100949d1bdf6badaa367307013a42d0708e 100644
--- a/src/parallel_io.c
+++ b/src/parallel_io.c
@@ -51,6 +51,8 @@
 #include "hydro_properties.h"
 #include "io_properties.h"
 #include "memuse.h"
+#include "output_list.h"
+#include "output_options.h"
 #include "part.h"
 #include "part_type.h"
 #include "star_formation_io.h"
@@ -1055,7 +1057,8 @@ void prepare_file(struct engine* e, const char* fileName,
   const struct gpart* gparts = e->s->gparts;
   const struct spart* sparts = e->s->sparts;
   const struct bpart* bparts = e->s->bparts;
-  struct swift_params* params = e->parameter_file;
+  struct output_options* output_options = e->output_options;
+  struct output_list* output_list = e->output_list_snapshots;
   const int with_cosmology = e->policy & engine_policy_cosmology;
   const int with_cooling = e->policy & engine_policy_cooling;
   const int with_temperature = e->policy & engine_policy_temperature;
@@ -1100,6 +1103,16 @@ void prepare_file(struct engine* e, const char* fileName,
                          e->s->dim[1] * factor_length,
                          e->s->dim[2] * factor_length};
 
+  /* Determine if we are writing a reduced snapshot, and if so which
+   * output selection type to use */
+  char current_selection_name[FIELD_BUFFER_SIZE] =
+      select_output_header_default_name;
+  if (output_list) {
+    /* Users could have specified a different Select Output scheme for each
+     * snapshot. */
+    output_list_get_current_select_output(output_list, current_selection_name);
+  }
+
   /* Print the relevant information and print status */
   io_write_attribute(h_grp, "BoxSize", DOUBLE, dim, 3);
   io_write_attribute(h_grp, "Time", DOUBLE, &dblTime, 1);
@@ -1140,6 +1153,7 @@ void prepare_file(struct engine* e, const char* fileName,
   io_write_attribute(h_grp, "NumFilesPerSnapshot", INT, &numFiles, 1);
   io_write_attribute_i(h_grp, "ThisFile", 0);
   io_write_attribute_s(h_grp, "OutputType", "Snapshot");
+  io_write_attribute_s(h_grp, "SelectOutput", current_selection_name);
 
   /* Close header */
   H5Gclose(h_grp);
@@ -1257,13 +1271,12 @@ void prepare_file(struct engine* e, const char* fileName,
     }
 
     /* Prepare everything that is not cancelled */
+
     for (int i = 0; i < num_fields; ++i) {
 
       /* Did the user cancel this field? */
-      char field[PARSER_MAX_LINE_SIZE];
-      sprintf(field, "SelectOutput:%.*s_%s", FIELD_BUFFER_SIZE, list[i].name,
-              part_type_names[ptype]);
-      int should_write = parser_get_opt_param_int(params, field, 1);
+      const int should_write = output_options_should_write_field(
+          output_options, current_selection_name, list[i].name, ptype);
 
       if (should_write)
         prepare_array_parallel(e, h_grp, fileName, xmfFile, partTypeGroupName,
@@ -1315,7 +1328,8 @@ void write_output_parallel(struct engine* e,
   const struct gpart* gparts = e->s->gparts;
   const struct spart* sparts = e->s->sparts;
   const struct bpart* bparts = e->s->bparts;
-  struct swift_params* params = e->parameter_file;
+  struct output_options* output_options = e->output_options;
+  struct output_list* output_list = e->output_list_snapshots;
   const int with_cosmology = e->policy & engine_policy_cosmology;
   const int with_cooling = e->policy & engine_policy_cooling;
   const int with_temperature = e->policy & engine_policy_temperature;
@@ -1767,13 +1781,19 @@ void write_output_parallel(struct engine* e,
     }
 
     /* Write everything that is not cancelled */
+    char current_selection_name[FIELD_BUFFER_SIZE] =
+        select_output_header_default_name;
+    if (output_list) {
+      /* Users could have specified a different Select Output scheme for each
+       * snapshot. */
+      output_list_get_current_select_output(output_list,
+                                            current_selection_name);
+    }
     for (int i = 0; i < num_fields; ++i) {
 
       /* Did the user cancel this field? */
-      char field[PARSER_MAX_LINE_SIZE];
-      sprintf(field, "SelectOutput:%.*s_%s", FIELD_BUFFER_SIZE, list[i].name,
-              part_type_names[ptype]);
-      int should_write = parser_get_opt_param_int(params, field, 1);
+      const int should_write = output_options_should_write_field(
+          output_options, current_selection_name, list[i].name, ptype);
 
       if (should_write)
         write_array_parallel(e, h_grp, fileName, partTypeGroupName, list[i],
diff --git a/src/serial_io.c b/src/serial_io.c
index 478058299a7b777b1351adcf64ddbf96b7bfb2bc..00d2f3c84dc923c7c359b6ab8549119c5836e201 100644
--- a/src/serial_io.c
+++ b/src/serial_io.c
@@ -51,6 +51,8 @@
 #include "hydro_properties.h"
 #include "io_properties.h"
 #include "memuse.h"
+#include "output_list.h"
+#include "output_options.h"
 #include "part.h"
 #include "part_type.h"
 #include "star_formation_io.h"
@@ -868,7 +870,8 @@ void write_output_serial(struct engine* e,
   const struct gpart* gparts = e->s->gparts;
   const struct spart* sparts = e->s->sparts;
   const struct bpart* bparts = e->s->bparts;
-  struct swift_params* params = e->parameter_file;
+  struct output_options* output_options = e->output_options;
+  struct output_list* output_list = e->output_list_snapshots;
   const int with_cosmology = e->policy & engine_policy_cosmology;
   const int with_cooling = e->policy & engine_policy_cooling;
   const int with_temperature = e->policy & engine_policy_temperature;
@@ -919,6 +922,17 @@ void write_output_serial(struct engine* e,
                            e->snapshot_output_count, e->snapshot_subdir,
                            e->snapshot_base_name);
 
+  /* Determine if we are writing a reduced snapshot, and if so which
+   * output selection type to use. Can just create a copy of this on
+   * each rank. */
+  char current_selection_name[FIELD_BUFFER_SIZE] =
+      select_output_header_default_name;
+  if (output_list) {
+    /* Users could have specified a different Select Output scheme for each
+     * snapshot. */
+    output_list_get_current_select_output(output_list, current_selection_name);
+  }
+
   /* Compute offset in the file and total number of particles */
   size_t N[swift_type_count] = {Ngas_written,   Ndm_written,
                                 Ndm_background, 0,
@@ -1007,6 +1021,7 @@ void write_output_serial(struct engine* e,
     io_write_attribute(h_grp, "NumFilesPerSnapshot", INT, &numFiles, 1);
     io_write_attribute_i(h_grp, "ThisFile", 0);
     io_write_attribute_s(h_grp, "OutputType", "Snapshot");
+    io_write_attribute_s(h_grp, "SelectOutput", current_selection_name);
 
     /* Close header */
     H5Gclose(h_grp);
@@ -1386,13 +1401,12 @@ void write_output_serial(struct engine* e,
         }
 
         /* Write everything that is not cancelled */
+
         for (int i = 0; i < num_fields; ++i) {
 
           /* Did the user cancel this field? */
-          char field[PARSER_MAX_LINE_SIZE];
-          sprintf(field, "SelectOutput:%.*s_%s", FIELD_BUFFER_SIZE,
-                  list[i].name, part_type_names[ptype]);
-          int should_write = parser_get_opt_param_int(params, field, 1);
+          const int should_write = output_options_should_write_field(
+              output_options, current_selection_name, list[i].name, ptype);
 
           if (should_write)
             write_array_serial(e, h_grp, fileName, xmfFile, partTypeGroupName,
diff --git a/src/single_io.c b/src/single_io.c
index d25da91279e5f89ed456c9519d913a39ff05596a..5560b90823c174d6bec63911cf12d5bc16e1bc83 100644
--- a/src/single_io.c
+++ b/src/single_io.c
@@ -50,6 +50,8 @@
 #include "hydro_properties.h"
 #include "io_properties.h"
 #include "memuse.h"
+#include "output_list.h"
+#include "output_options.h"
 #include "part.h"
 #include "part_type.h"
 #include "star_formation_io.h"
@@ -727,7 +729,8 @@ void write_output_single(struct engine* e,
   const struct gpart* gparts = e->s->gparts;
   const struct spart* sparts = e->s->sparts;
   const struct bpart* bparts = e->s->bparts;
-  struct swift_params* params = e->parameter_file;
+  struct output_options* output_options = e->output_options;
+  struct output_list* output_list = e->output_list_snapshots;
   const int with_cosmology = e->policy & engine_policy_cosmology;
   const int with_cooling = e->policy & engine_policy_cooling;
   const int with_temperature = e->policy & engine_policy_temperature;
@@ -812,6 +815,16 @@ void write_output_single(struct engine* e,
                          e->s->dim[1] * factor_length,
                          e->s->dim[2] * factor_length};
 
+  /* Determine if we are writing a reduced snapshot, and if so which
+   * output selection type to use */
+  char current_selection_name[FIELD_BUFFER_SIZE] =
+      select_output_header_default_name;
+  if (output_list) {
+    /* Users could have specified a different Select Output scheme for each
+     * snapshot. */
+    output_list_get_current_select_output(output_list, current_selection_name);
+  }
+
   /* Print the relevant information and print status */
   io_write_attribute(h_grp, "BoxSize", DOUBLE, dim, 3);
   io_write_attribute(h_grp, "Time", DOUBLE, &dblTime, 1);
@@ -852,6 +865,7 @@ void write_output_single(struct engine* e,
   io_write_attribute(h_grp, "NumFilesPerSnapshot", INT, &numFiles, 1);
   io_write_attribute_i(h_grp, "ThisFile", 0);
   io_write_attribute_s(h_grp, "OutputType", "Snapshot");
+  io_write_attribute_s(h_grp, "SelectOutput", current_selection_name);
 
   /* Close header */
   H5Gclose(h_grp);
@@ -1041,7 +1055,7 @@ void write_output_single(struct engine* e,
         if (swift_memalign("gparts_written", (void**)&gparts_written,
                            gpart_align,
                            Ndm_background * sizeof(struct gpart)) != 0)
-          error("Error while allocating temporart memory for gparts");
+          error("Error while allocating temporary memory for gparts");
 
         if (with_stf) {
           if (swift_memalign(
@@ -1049,7 +1063,7 @@ void write_output_single(struct engine* e,
                   gpart_align,
                   Ndm_background * sizeof(struct velociraptor_gpart_data)) != 0)
             error(
-                "Error while allocating temporart memory for gparts STF "
+                "Error while allocating temporary memory for gparts STF "
                 "data");
         }
 
@@ -1072,7 +1086,7 @@ void write_output_single(struct engine* e,
       case swift_type_stars: {
         if (Nstars == Nstars_written) {
 
-          /* No inhibted particles: easy case */
+          /* No inhibited particles: easy case */
           N = Nstars;
           stars_write_particles(sparts, list, &num_fields, with_cosmology);
           num_fields += chemistry_write_sparticles(sparts, list + num_fields);
@@ -1123,7 +1137,7 @@ void write_output_single(struct engine* e,
       case swift_type_black_hole: {
         if (Nblackholes == Nblackholes_written) {
 
-          /* No inhibted particles: easy case */
+          /* No inhibited particles: easy case */
           N = Nblackholes;
           black_holes_write_particles(bparts, list, &num_fields,
                                       with_cosmology);
@@ -1169,17 +1183,17 @@ void write_output_single(struct engine* e,
     }
 
     /* Write everything that is not cancelled */
+
     for (int i = 0; i < num_fields; ++i) {
 
       /* Did the user cancel this field? */
-      char field[PARSER_MAX_LINE_SIZE];
-      sprintf(field, "SelectOutput:%.*s_%s", FIELD_BUFFER_SIZE, list[i].name,
-              part_type_names[ptype]);
-      int should_write = parser_get_opt_param_int(params, field, 1);
+      const int should_write = output_options_should_write_field(
+          output_options, current_selection_name, list[i].name, ptype);
 
-      if (should_write)
+      if (should_write) {
         write_array_single(e, h_grp, fileName, xmfFile, partTypeGroupName,
                            list[i], N, internal_units, snapshot_units);
+      }
     }
 
     /* Free temporary arrays */
diff --git a/src/swift.h b/src/swift.h
index 3c659aa7dce1d78d9887187da31e4ee500ad6848..ad16acf3758abfe783a7bf986e15691f1c9096c6 100644
--- a/src/swift.h
+++ b/src/swift.h
@@ -60,7 +60,8 @@
 #include "minmax.h"
 #include "mpiuse.h"
 #include "multipole.h"
-#include "outputlist.h"
+#include "output_list.h"
+#include "output_options.h"
 #include "parallel_io.h"
 #include "parser.h"
 #include "part.h"
diff --git a/tests/selectOutput.yml b/tests/selectOutput.yml
index 1778935146b19992e25efcb320d8cc523c6472a5..bbcbadff8f66d4b0c6259c8725f430869bb31ac7 100644
--- a/tests/selectOutput.yml
+++ b/tests/selectOutput.yml
@@ -1,12 +1,7 @@
-SelectOutput:
+Default:
   # Particle Type 0
-  Coordinates_Gas: 1   # check if written when specified
-  Velocities_Gas:  0   # check if not written when specified
-  Masses_Gas:     -5   # check warning if not 0 or 1 and if written
-  Pot_Gas:         1   # check warning if wrong name
+  Coordinates_Gas: on   # check if written when specified
+  Velocities_Gas:  off   # check if not written when specified
+  Pot_Gas:         on   # check warning if wrong name
   # Density_Gas:   1   # check if written when not specified
 
-# Parameters for the hydrodynamics scheme
-SPH:
-  resolution_eta:        1.2348   # Target smoothing length in units of the mean inter-particle separation (1.2348 == 48Ngbs with the cubic spline kernel).
-  CFL_condition:         0.1      # Courant-Friedrich-Levy condition for time integration.
diff --git a/tests/selectOutputParameters.yml b/tests/selectOutputParameters.yml
new file mode 100644
index 0000000000000000000000000000000000000000..984b105b02259205415f68b704bd97c0b24de424
--- /dev/null
+++ b/tests/selectOutputParameters.yml
@@ -0,0 +1,8 @@
+# Parameters for the hydrodynamics scheme
+SPH:
+  resolution_eta:        1.2348   # Target smoothing length in units of the mean inter-particle separation (1.2348 == 48Ngbs with the cubic spline kernel).
+  CFL_condition:         0.1      # Courant-Friedrich-Levy condition for time integration.
+
+Snapshots:
+  select_output_on:      1
+  select_output:         selectOutput.yml
diff --git a/tests/testSelectOutput.c b/tests/testSelectOutput.c
index 4e16df60d882f2bde04d91ace429f5526ba5d7f2..6ff3f59e1813528d3139f1a14a2db22b1a6e9934 100644
--- a/tests/testSelectOutput.c
+++ b/tests/testSelectOutput.c
@@ -26,12 +26,14 @@
 void select_output_engine_init(struct engine *e, struct space *s,
                                struct cosmology *cosmo,
                                struct swift_params *params,
+                               struct output_options *output,
                                struct cooling_function_data *cooling,
                                struct hydro_props *hydro_properties) {
   /* set structures */
   e->s = s;
   e->cooling_func = cooling;
   e->parameter_file = params;
+  e->output_options = output;
   e->cosmology = cosmo;
   e->policy = engine_policy_hydro;
   e->hydro_properties = hydro_properties;
@@ -96,9 +98,12 @@ int main(int argc, char *argv[]) {
   /* parse parameters */
   message("Reading parameters.");
   struct swift_params param_file;
-  const char *input_file = "selectOutput.yml";
+  const char *input_file = "selectOutputParameters.yml";
   parser_read_file(input_file, &param_file);
 
+  struct output_options output_options;
+  output_options_init(&param_file, 0, &output_options);
+
   /* Default unit system */
   message("Initialization of the unit system.");
   struct unit_system us;
@@ -149,14 +154,14 @@ int main(int argc, char *argv[]) {
   e.physical_constants = &prog_const;
   sprintf(e.snapshot_base_name, "testSelectOutput");
   sprintf(e.run_name, "Select Output Test");
-  select_output_engine_init(&e, &s, &cosmo, &param_file, &cooling,
-                            &hydro_properties);
+  select_output_engine_init(&e, &s, &cosmo, &param_file, &output_options,
+                            &cooling, &hydro_properties);
 
   /* check output selection */
   message("Checking output parameters.");
   long long N_total[swift_type_count] = {
       (long long)Ngas, (long long)Ngpart, 0, 0, (long long)Nspart, 0};
-  io_check_output_fields(&param_file, N_total);
+  io_check_output_fields(&param_file, N_total, /*with_cosmology=*/0);
 
   /* write output file */
   message("Writing output.");
diff --git a/tests/testSelectOutput.py b/tests/testSelectOutput.py
index 97e1c865d133b161c5661c7ac63f728065461bd6..79316d867a35c10214e7138618bfe61943fa2470 100644
--- a/tests/testSelectOutput.py
+++ b/tests/testSelectOutput.py
@@ -1,22 +1,22 @@
 ###############################################################################
- # This file is part of SWIFT.
- # Copyright (c) 2015 Bert Vandenbroucke (bert.vandenbroucke@ugent.be)
- #                    Matthieu Schaller (matthieu.schaller@durham.ac.uk)
- # 
- # This program is free software: you can redistribute it and/or modify
- # it under the terms of the GNU Lesser General Public License as published
- # by the Free Software Foundation, either version 3 of the License, or
- # (at your option) any later version.
- # 
- # This program is distributed in the hope that it will be useful,
- # but WITHOUT ANY WARRANTY; without even the implied warranty of
- # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
- # GNU General Public License for more details.
- # 
- # You should have received a copy of the GNU Lesser General Public License
- # along with this program.  If not, see <http://www.gnu.org/licenses/>.
- # 
- ##############################################################################
+# This file is part of SWIFT.
+# Copyright (c) 2015 Bert Vandenbroucke (bert.vandenbroucke@ugent.be)
+#                    Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU Lesser General Public License as published
+# by the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public License
+# along with this program.  If not, see <http://www.gnu.org/licenses/>.
+#
+##############################################################################
 
 # Check the output done with swift
 
@@ -47,8 +47,8 @@ if "Densities" not in part0:
 with open(log_filename, "r") as f:
     data = f.read()
 
-if "SelectOutput:Masses_Gas" not in data:
-    raise Exception("Input error in `SelectOuput:Masses_Gas` not detected")
+if "Default:Masses_Gas" not in data:
+    raise Exception("Input error in `Default:Masses_Gas` not detected")
 
-if "SelectOutput:Pot_Gas" not in data:
+if "Default:Pot_Gas" not in data:
     raise Exception("Parameter name error not detected for `SelectOutput:Pot_Gas`")