diff --git a/doc/RTD/source/ParameterFiles/output_selection.rst b/doc/RTD/source/ParameterFiles/output_selection.rst
index b9f4ceb6e9a31d376e85d35cff6511f62b0e29af..44fe7bcb362fb55cdd0f8c4c12b5e3e78a33523e 100644
--- a/doc/RTD/source/ParameterFiles/output_selection.rst
+++ b/doc/RTD/source/ParameterFiles/output_selection.rst
@@ -145,7 +145,7 @@ selection:
      basename: bh
      subdir: snip
 
-This will put the outputs corresponding to the ``BlackHolesOnly`` selection into
+This will put the outputs corresponding to this "BlackHolesOnly" selection into
 a sub-directory called ``snip`` and have the files themselves called
 ``bh_0000.hdf5`` where the number corresponds to the global number of
 snapshots. The counter is global and not reset for each type of selection.
@@ -155,6 +155,19 @@ The sub-directories are created when writing the first snapshot of a given
 category; the onus is hence on the user to ensure correct writing permissions
 ahead of that time.
 
+Finally, it is possible to specify individual sub-sampling ratios for each
+output selection:
+
+.. code:: YAML
+
+   Default:
+     subsample: [0, 1, 0, 0, 0, 0, 1]  # Sub-sample the DM and neutrinos
+     subsmple_fractions: [0, 0.01, 0, 0, 0, 0, 0.1]
+
+If these keywords are omitted then the code will use the default values
+specified in the ``Snapshots`` section of the main parameter file.
+
+
 Combining Output Lists and Output Selection
 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
@@ -187,6 +200,8 @@ cousins. To do this, we will define two top-level sections in our
     Element_Mass_Fractions_Gas: off
     Densities_Gas: FMantissa9
     basename: snip
+    subsample: [0, 1, 0, 0, 0, 0, 1]  # Sub-sample the DM and neutrinos
+    subsmple_fractions: [0, 0.01, 0, 0, 0, 0, 0.1]
     ...
 
 To then select which outputs are 'snapshots' and which are 'snipshots', you
diff --git a/doc/RTD/source/ParameterFiles/parameter_description.rst b/doc/RTD/source/ParameterFiles/parameter_description.rst
index 2d80536d59b5a4d6a9aa362ff4ab2523dc41dd07..b97d399108cc841561fdc5a61c72727270a19f3a 100644
--- a/doc/RTD/source/ParameterFiles/parameter_description.rst
+++ b/doc/RTD/source/ParameterFiles/parameter_description.rst
@@ -792,6 +792,21 @@ of MPI ranks used in a given run. The individual files of snapshot 1234 will
 have the name ``base_name_1234.x.hdf5`` where when running on N MPI ranks, ``x``
 runs from 0 to N-1.
 
+Users can optionally ask to randomly sub-sample the particles in the snapshots.
+This is specified for each particle type individually:
+
+* Whether to switch on sub-sampling: ``subsample``   
+* Whether to switch on sub-sampling: ``subsample_fraction`` 
+
+These are arrays of 7 elements defaulting to seven 0s if left unspecified. Each
+entry corresponds to the particle type used in the initial conditions and
+snapshots [#f3]_.  The ``subsample`` array is made of ``0`` and ``1`` to indicate which
+particle types to subsample. The other array is a float between ``0`` and ``1``
+indicating the fraction of particles to keep in the outputs.  Note that the
+selection of particles is selected randomly for each individual
+snapshot. Particles can hence not be traced back from output to output when this
+is switched on.
+  
 Users can optionally specify the level of compression used by the HDF5 library
 using the parameter:
 
@@ -855,6 +870,7 @@ would have:
      time_first:          0.01
      delta_time:          0.005
      invoke_stf:          0
+     invoke_fof:          1
      compression:         3
      distributed:         1
      UnitLength_in_cgs:   1.  # Use cm in outputs
@@ -862,6 +878,10 @@ would have:
      UnitVelocity_in_cgs: 1.  # Use cm/s in outputs
      UnitCurrent_in_cgs:  1.  # Use Ampere in outputs
      UnitTemp_in_cgs:     1.  # Use Kelvin in outputs
+     subsample:           [0, 1, 0, 0, 0, 0, 1]   # Sub-sample the DM and neutrinos
+     subsample_fraction:  [0, 0.01, 0, 0, 0, 0, 0.1]  # Write 1% of the DM parts and 10% of the neutrinos
+     run_on_dump:         1
+     dump_command:        ./submit_analysis.sh
 
 Some additional specific options for the snapshot outputs are described in the
 following pages:
@@ -1486,14 +1506,6 @@ and all the gparts are not active during the timestep of the snapshot dump, the
 exact forces computation is performed on the first timestep at which all the
 gparts are active after that snapshot output timestep.
 
-
-------------------------
-
-.. [#f1] The thorough reader (or overly keen SWIFT tester) would find  that the speed of light is :math:`c=1.8026\times10^{12}\,\rm{fur}\,\rm{ftn}^{-1}`, Newton's constant becomes :math:`G_N=4.896735\times10^{-4}~\rm{fur}^3\,\rm{fir}^{-1}\,\rm{ftn}^{-2}` and Planck's constant turns into :math:`h=4.851453\times 10^{-34}~\rm{fur}^2\,\rm{fir}\,\rm{ftn}^{-1}`.
-
-
-.. [#f2] which would translate into a constant :math:`G_N=1.5517771\times10^{-9}~cm^{3}\,g^{-1}\,s^{-2}` if expressed in the CGS system.
-
 Neutrinos
 ---------
 
@@ -1514,3 +1526,13 @@ A complete specification of the model looks like
     generate_ics:  1    # Replace neutrino particle velocities with random Fermi-Dirac momenta at the start
     use_delta_f:   1    # Use the delta-f method for shot noise reduction
     neutrino_seed: 1234 # A random seed used for the Fermi-Dirac momenta
+
+
+------------------------
+    
+.. [#f1] The thorough reader (or overly keen SWIFT tester) would find  that the speed of light is :math:`c=1.8026\times10^{12}\,\rm{fur}\,\rm{ftn}^{-1}`, Newton's constant becomes :math:`G_N=4.896735\times10^{-4}~\rm{fur}^3\,\rm{fir}^{-1}\,\rm{ftn}^{-2}` and Planck's constant turns into :math:`h=4.851453\times 10^{-34}~\rm{fur}^2\,\rm{fir}\,\rm{ftn}^{-1}`.
+
+
+.. [#f2] which would translate into a constant :math:`G_N=1.5517771\times10^{-9}~cm^{3}\,g^{-1}\,s^{-2}` if expressed in the CGS system.
+
+.. [#f3] The mapping is 0 --> gas, 1 --> dark matter, 2 --> background dark matter, 3 --> sinks, 4 --> stars, 5 --> black holes, 6 --> neutrinos.
diff --git a/doc/RTD/source/Snapshots/index.rst b/doc/RTD/source/Snapshots/index.rst
index 64cb94fae991870f33c1999011315c9cb4ec7af7..880be8ad24cc5f10024f026c343b14f9a69e2aab 100644
--- a/doc/RTD/source/Snapshots/index.rst
+++ b/doc/RTD/source/Snapshots/index.rst
@@ -35,16 +35,38 @@ this sub-snapshot file when the user asked for distributed snapshots (see
 ``NumPart_Total``. Note however, that there is no high-word for this field. We
 store it as a 64-bits integer [#f1]_. The field ``NumFilesPerSnapshot`` specifies the
 number of sub-snapshot files (always 1 unless a distributed snapshot was asked
-for). 
+for) and ``ThisFile`` the id of that specific file (always 0 unless a distributed
+snapshot was asked for). 
 
 The field ``InitialMassTable`` contains the *mean* initial mass of each of the
 particle types present in the initial conditions. This can be used as estimator
 of the mass resolution of the run. The masses are expressed in internal units.
 
+The field ``OutputType`` contains information about the kind of output this
+snapshot is. The possible values are:
+
++---------------------+-----------------------------------------------------+
+| OutputType          | Definition                                          |
++=====================+=====================================================+
+| ``FullVolume``      | Regular vanilla snapshot                            |
++---------------------+-----------------------------------------------------+
+| ``SubSampled``      | Snapshot where some particle types were sub-sampled |
++---------------------+-----------------------------------------------------+
+| ``LineOfSight``     | Line-of-sight snapshot                              |
++---------------------+-----------------------------------------------------+
+| ``FOF``             | Friends-Of-Friends Halo Catalogue                   |
++---------------------+-----------------------------------------------------+
+
+
 The ``RunName`` field contains the name of the simulation that was specified as
 the ``run_name`` in the :ref:`Parameters_meta_data` section of the YAML
 parameter file.
 
+The ``System`` field contains the name of the machine where the MPI rank 0 was
+placed. This name is whatever UNIX's ``gethostname()`` function returns on that
+system. Similarly, the ``SnapshotDate`` field contains the date and time when
+the file was written.
+
 The ``TimeBase_dloga`` field contains the change in logarithm of the
 scale-factor corresponding to a time-step of length 1 on the integer
 time-line. This is the smallest time-step size that the code can use. This field
@@ -54,6 +76,17 @@ would be the increase in time a particle in the time-bin one would have. Note
 that in cosmological runs this quantity evolves with redhsift as the (logarithm
 of the) scale-factor is used on the integer time-line.
 
+The field ``SelectOutput`` will contain the name of the
+:ref:`Output_selection_label` used for this specific output and will take the value
+``Default`` if no such selection (or the default one) was used.
+
+If a sub-sampling of the particle fields was used, then the header additionally
+contains a field describing the fraction of the particles of each type that were
+written to the snapshot. Note, however, that when sub-sampling the fields 
+``NumPart_Total``, ``NumPart_HighWord``, and ``NumPart_ThisFile`` contain the number
+of particles actually written (i.e. after sub-sampling), not the total number of
+particles in the run.
+
 Meta-data about the code and run
 --------------------------------
 
@@ -365,6 +398,10 @@ In the case of a single-file snapshot, the ``Files`` array is just an array of
 zeroes since all the particles will be in the 0-th file. Note also that in the
 case of a multi-files snapshot, a cell is always contained in a single file.
 
+If a snapshot used a sub-sampled output, then the counts and offsets are
+adjusted accordingly and correspond to the actual content of the file
+(i.e. after the sub-sampling was applied).
+
 As an example, if one is interested in retriving all the densities of the gas
 particles in the cell around the position `[1, 1, 1]` in a single-file
 snapstshot one could use a piece of code similar to:
diff --git a/src/common_io.c b/src/common_io.c
index d61fe52e07dce14e6d3b3c45bf6ed4738cfee7dc..9fc9fccd4f3985fda7c9f0258b30c89dfd4e5edc 100644
--- a/src/common_io.c
+++ b/src/common_io.c
@@ -907,20 +907,6 @@ void io_prepare_dm_background_gparts(struct threadpool* tp,
                  sizeof(struct gpart), threadpool_auto_chunk_size, NULL);
 }
 
-size_t io_count_dm_background_gparts(const struct gpart* const gparts,
-                                     const size_t Ndm) {
-
-  swift_declare_aligned_ptr(const struct gpart, gparts_array, gparts,
-                            SWIFT_STRUCT_ALIGNMENT);
-
-  size_t count = 0;
-  for (size_t i = 0; i < Ndm; ++i) {
-    if (gparts_array[i].type == swift_type_dark_matter_background) ++count;
-  }
-
-  return count;
-}
-
 void io_prepare_dm_neutrino_gparts_mapper(void* restrict data, int Ndm,
                                           void* dummy) {
 
@@ -957,20 +943,6 @@ void io_prepare_dm_neutrino_gparts(struct threadpool* tp,
                  sizeof(struct gpart), threadpool_auto_chunk_size, NULL);
 }
 
-size_t io_count_dm_neutrino_gparts(const struct gpart* const gparts,
-                                   const size_t Ndm) {
-
-  swift_declare_aligned_ptr(const struct gpart, gparts_array, gparts,
-                            SWIFT_STRUCT_ALIGNMENT);
-
-  size_t count = 0;
-  for (size_t i = 0; i < Ndm; ++i) {
-    if (gparts_array[i].type == swift_type_neutrino) ++count;
-  }
-
-  return count;
-}
-
 struct duplication_data {
 
   struct part* parts;
@@ -1221,12 +1193,17 @@ void io_duplicate_black_holes_gparts(struct threadpool* tp,
 /**
  * @brief Copy every non-inhibited #part into the parts_written array.
  *
+ * Also takes into account possible downsampling.
+ *
  * @param parts The array of #part containing all particles.
  * @param xparts The array of #xpart containing all particles.
  * @param parts_written The array of #part to fill with particles we want to
  * write.
  * @param xparts_written The array of #xpart  to fill with particles we want to
  * write.
+ * @param subsample Are we subsampling the particles?
+ * @param subsample_ratio The fraction of particles to write if subsampling.
+ * @param snap_num The snapshot ID (used to seed the RNG when sub-sampling).
  * @param Nparts The total number of #part.
  * @param Nparts_written The total number of #part to write.
  */
@@ -1234,7 +1211,8 @@ void io_collect_parts_to_write(const struct part* restrict parts,
                                const struct xpart* restrict xparts,
                                struct part* restrict parts_written,
                                struct xpart* restrict xparts_written,
-                               const size_t Nparts,
+                               const int subsample, const float subsample_ratio,
+                               const int snap_num, const size_t Nparts,
                                const size_t Nparts_written) {
 
   size_t count = 0;
@@ -1246,6 +1224,14 @@ void io_collect_parts_to_write(const struct part* restrict parts,
     if (parts[i].time_bin != time_bin_inhibited &&
         parts[i].time_bin != time_bin_not_created) {
 
+      /* When subsampling, select particles at random */
+      if (subsample) {
+        const float r = random_unit_interval(parts[i].id, snap_num,
+                                             random_number_snapshot_sampling);
+
+        if (r > subsample_ratio) continue;
+      }
+
       parts_written[count] = parts[i];
       xparts_written[count] = xparts[i];
       count++;
@@ -1261,14 +1247,21 @@ void io_collect_parts_to_write(const struct part* restrict parts,
 /**
  * @brief Copy every non-inhibited #spart into the sparts_written array.
  *
+ * Also takes into account possible downsampling.
+ *
  * @param sparts The array of #spart containing all particles.
  * @param sparts_written The array of #spart to fill with particles we want to
  * write.
+ * @param subsample Are we subsampling the particles?
+ * @param subsample_ratio The fraction of particles to write if subsampling.
+ * @param snap_num The snapshot ID (used to seed the RNG when sub-sampling).
  * @param Nsparts The total number of #part.
  * @param Nsparts_written The total number of #part to write.
  */
 void io_collect_sparts_to_write(const struct spart* restrict sparts,
                                 struct spart* restrict sparts_written,
+                                const int subsample,
+                                const float subsample_ratio, const int snap_num,
                                 const size_t Nsparts,
                                 const size_t Nsparts_written) {
 
@@ -1281,6 +1274,14 @@ void io_collect_sparts_to_write(const struct spart* restrict sparts,
     if (sparts[i].time_bin != time_bin_inhibited &&
         sparts[i].time_bin != time_bin_not_created) {
 
+      /* When subsampling, select particles at random */
+      if (subsample) {
+        const float r = random_unit_interval(sparts[i].id, snap_num,
+                                             random_number_snapshot_sampling);
+
+        if (r > subsample_ratio) continue;
+      }
+
       sparts_written[count] = sparts[i];
       count++;
     }
@@ -1295,15 +1296,21 @@ void io_collect_sparts_to_write(const struct spart* restrict sparts,
 /**
  * @brief Copy every non-inhibited #sink into the sinks_written array.
  *
+ * Also takes into account possible downsampling.
+ *
  * @param sinks The array of #sink containing all particles.
  * @param sinks_written The array of #sink to fill with particles we want to
  * write.
+ * @param subsample Are we subsampling the particles?
+ * @param subsample_ratio The fraction of particles to write if subsampling.
+ * @param snap_num The snapshot ID (used to seed the RNG when sub-sampling).
  * @param Nsinks The total number of #sink.
  * @param Nsinks_written The total number of #sink to write.
  */
 void io_collect_sinks_to_write(const struct sink* restrict sinks,
                                struct sink* restrict sinks_written,
-                               const size_t Nsinks,
+                               const int subsample, const float subsample_ratio,
+                               const int snap_num, const size_t Nsinks,
                                const size_t Nsinks_written) {
 
   size_t count = 0;
@@ -1315,6 +1322,14 @@ void io_collect_sinks_to_write(const struct sink* restrict sinks,
     if (sinks[i].time_bin != time_bin_inhibited &&
         sinks[i].time_bin != time_bin_not_created) {
 
+      /* When subsampling, select particles at random */
+      if (subsample) {
+        const float r = random_unit_interval(sinks[i].id, snap_num,
+                                             random_number_snapshot_sampling);
+
+        if (r > subsample_ratio) continue;
+      }
+
       sinks_written[count] = sinks[i];
       count++;
     }
@@ -1329,14 +1344,21 @@ void io_collect_sinks_to_write(const struct sink* restrict sinks,
 /**
  * @brief Copy every non-inhibited #bpart into the bparts_written array.
  *
+ * Also takes into account possible downsampling.
+ *
  * @param bparts The array of #bpart containing all particles.
  * @param bparts_written The array of #bpart to fill with particles we want to
  * write.
+ * @param subsample Are we subsampling the particles?
+ * @param subsample_ratio The fraction of particles to write if subsampling.
+ * @param snap_num The snapshot ID (used to seed the RNG when sub-sampling).
  * @param Nbparts The total number of #part.
  * @param Nbparts_written The total number of #part to write.
  */
 void io_collect_bparts_to_write(const struct bpart* restrict bparts,
                                 struct bpart* restrict bparts_written,
+                                const int subsample,
+                                const float subsample_ratio, const int snap_num,
                                 const size_t Nbparts,
                                 const size_t Nbparts_written) {
 
@@ -1349,6 +1371,14 @@ void io_collect_bparts_to_write(const struct bpart* restrict bparts,
     if (bparts[i].time_bin != time_bin_inhibited &&
         bparts[i].time_bin != time_bin_not_created) {
 
+      /* When subsampling, select particles at random */
+      if (subsample) {
+        const float r = random_unit_interval(bparts[i].id, snap_num,
+                                             random_number_snapshot_sampling);
+
+        if (r > subsample_ratio) continue;
+      }
+
       bparts_written[count] = bparts[i];
       count++;
     }
@@ -1364,12 +1394,17 @@ void io_collect_bparts_to_write(const struct bpart* restrict bparts,
  * @brief Copy every non-inhibited regulat DM #gpart into the gparts_written
  * array.
  *
+ * Also takes into account possible downsampling.
+ *
  * @param gparts The array of #gpart containing all particles.
  * @param vr_data The array of gpart-related VELOCIraptor output.
  * @param gparts_written The array of #gpart to fill with particles we want to
  * write.
  * @param vr_data_written The array of gpart-related VELOCIraptor with particles
  * we want to write.
+ * @param subsample Are we subsampling the particles?
+ * @param subsample_ratio The fraction of particles to write if subsampling.
+ * @param snap_num The snapshot ID (used to seed the RNG when sub-sampling).
  * @param Ngparts The total number of #part.
  * @param Ngparts_written The total number of #part to write.
  * @param with_stf Are we running with STF? i.e. do we want to collect vr data?
@@ -1379,6 +1414,7 @@ void io_collect_gparts_to_write(
     const struct velociraptor_gpart_data* restrict vr_data,
     struct gpart* restrict gparts_written,
     struct velociraptor_gpart_data* restrict vr_data_written,
+    const int subsample, const float subsample_ratio, const int snap_num,
     const size_t Ngparts, const size_t Ngparts_written, const int with_stf) {
 
   size_t count = 0;
@@ -1391,6 +1427,15 @@ void io_collect_gparts_to_write(
         (gparts[i].time_bin != time_bin_not_created) &&
         (gparts[i].type == swift_type_dark_matter)) {
 
+      /* When subsampling, select particles at random */
+      if (subsample) {
+        const float r =
+            random_unit_interval(gparts[i].id_or_neg_offset, snap_num,
+                                 random_number_snapshot_sampling);
+
+        if (r > subsample_ratio) continue;
+      }
+
       if (with_stf) vr_data_written[count] = vr_data[i];
 
       gparts_written[count] = gparts[i];
@@ -1408,12 +1453,17 @@ void io_collect_gparts_to_write(
  * @brief Copy every non-inhibited background DM #gpart into the gparts_written
  * array.
  *
+ * Also takes into account possible downsampling.
+ *
  * @param gparts The array of #gpart containing all particles.
  * @param vr_data The array of gpart-related VELOCIraptor output.
  * @param gparts_written The array of #gpart to fill with particles we want to
  * write.
  * @param vr_data_written The array of gpart-related VELOCIraptor with particles
  * we want to write.
+ * @param subsample Are we subsampling the particles?
+ * @param subsample_ratio The fraction of particles to write if subsampling.
+ * @param snap_num The snapshot ID (used to seed the RNG when sub-sampling).
  * @param Ngparts The total number of #part.
  * @param Ngparts_written The total number of #part to write.
  * @param with_stf Are we running with STF? i.e. do we want to collect vr data?
@@ -1423,6 +1473,7 @@ void io_collect_gparts_background_to_write(
     const struct velociraptor_gpart_data* restrict vr_data,
     struct gpart* restrict gparts_written,
     struct velociraptor_gpart_data* restrict vr_data_written,
+    const int subsample, const float subsample_ratio, const int snap_num,
     const size_t Ngparts, const size_t Ngparts_written, const int with_stf) {
 
   size_t count = 0;
@@ -1435,6 +1486,15 @@ void io_collect_gparts_background_to_write(
         (gparts[i].time_bin != time_bin_not_created) &&
         (gparts[i].type == swift_type_dark_matter_background)) {
 
+      /* When subsampling, select particles at random */
+      if (subsample) {
+        const float r =
+            random_unit_interval(gparts[i].id_or_neg_offset, snap_num,
+                                 random_number_snapshot_sampling);
+
+        if (r > subsample_ratio) continue;
+      }
+
       if (with_stf) vr_data_written[count] = vr_data[i];
 
       gparts_written[count] = gparts[i];
@@ -1452,12 +1512,17 @@ void io_collect_gparts_background_to_write(
  * @brief Copy every non-inhibited neutrino DM #gpart into the gparts_written
  * array.
  *
+ * Also takes into account possible downsampling.
+ *
  * @param gparts The array of #gpart containing all particles.
  * @param vr_data The array of gpart-related VELOCIraptor output.
  * @param gparts_written The array of #gpart to fill with particles we want to
  * write.
  * @param vr_data_written The array of gpart-related VELOCIraptor with particles
  * we want to write.
+ * @param subsample Are we subsampling the particles?
+ * @param subsample_ratio The fraction of particles to write if subsampling.
+ * @param snap_num The snapshot ID (used to seed the RNG when sub-sampling).
  * @param Ngparts The total number of #part.
  * @param Ngparts_written The total number of #part to write.
  * @param with_stf Are we running with STF? i.e. do we want to collect vr data?
@@ -1467,6 +1532,7 @@ void io_collect_gparts_neutrino_to_write(
     const struct velociraptor_gpart_data* restrict vr_data,
     struct gpart* restrict gparts_written,
     struct velociraptor_gpart_data* restrict vr_data_written,
+    const int subsample, const float subsample_ratio, const int snap_num,
     const size_t Ngparts, const size_t Ngparts_written, const int with_stf) {
 
   size_t count = 0;
@@ -1479,6 +1545,15 @@ void io_collect_gparts_neutrino_to_write(
         (gparts[i].time_bin != time_bin_not_created) &&
         (gparts[i].type == swift_type_neutrino)) {
 
+      /* When subsampling, select particles at random */
+      if (subsample) {
+        const float r =
+            random_unit_interval(gparts[i].id_or_neg_offset, snap_num,
+                                 random_number_snapshot_sampling);
+
+        if (r > subsample_ratio) continue;
+      }
+
       if (with_stf) vr_data_written[count] = vr_data[i];
 
       gparts_written[count] = gparts[i];
diff --git a/src/common_io.h b/src/common_io.h
index fadfd8e83b2f2f12f707cbefdbeff3e193add838..73719f35c920203960afc094ea3060e04456e968 100644
--- a/src/common_io.h
+++ b/src/common_io.h
@@ -34,6 +34,7 @@
 
 /* Avoid cyclic inclusion problems */
 struct cell;
+struct space;
 struct part;
 struct gpart;
 struct velociraptor_gpart_data;
@@ -110,6 +111,9 @@ void io_write_cell_offsets(hid_t h_grp, const int cdim[3], const double dim[3],
                            const struct cell* cells_top, const int nr_cells,
                            const double width[3], const int nodeID,
                            const int distributed,
+                           const int subsample[swift_type_count],
+                           const float subsample_fraction[swift_type_count],
+                           const int snap_num,
                            const long long global_counts[swift_type_count],
                            const long long global_offsets[swift_type_count],
                            const int num_fields[swift_type_count],
@@ -132,42 +136,85 @@ void io_copy_temp_buffer(void* temp, const struct engine* e,
 size_t io_sizeof_type(enum IO_DATA_TYPE type);
 int io_is_double_precision(enum IO_DATA_TYPE type);
 
+long long io_count_gas_to_write(const struct space* s, const int subsample,
+                                const float subsample_ratio,
+                                const int snap_num);
+
+long long io_count_dark_matter_to_write(const struct space* s,
+                                        const int subsample,
+                                        const float subsample_ratio,
+                                        const int snap_num);
+
+long long io_count_background_dark_matter_to_write(const struct space* s,
+                                                   const int subsample,
+                                                   const float subsample_ratio,
+                                                   const int snap_num);
+
+long long io_count_stars_to_write(const struct space* s, const int subsample,
+                                  const float subsample_ratio,
+                                  const int snap_num);
+
+long long io_count_sinks_to_write(const struct space* s, const int subsample,
+                                  const float subsample_ratio,
+                                  const int snap_num);
+
+long long io_count_black_holes_to_write(const struct space* s,
+                                        const int subsample,
+                                        const float subsample_ratio,
+                                        const int snap_num);
+
+long long io_count_neutrinos_to_write(const struct space* s,
+                                      const int subsample,
+                                      const float subsample_ratio,
+                                      const int snap_num);
+
 void io_collect_parts_to_write(const struct part* restrict parts,
                                const struct xpart* restrict xparts,
                                struct part* restrict parts_written,
                                struct xpart* restrict xparts_written,
-                               const size_t Nparts,
+                               const int subsample, const float subsample_ratio,
+                               const int snap_num, const size_t Nparts,
                                const size_t Nparts_written);
 void io_collect_sinks_to_write(const struct sink* restrict sinks,
                                struct sink* restrict sinks_written,
-                               const size_t Nsinks,
+                               const int subsample, const float subsample_ratio,
+                               const int snap_num, const size_t Nsinks,
                                const size_t Nsinks_written);
 void io_collect_sparts_to_write(const struct spart* restrict sparts,
                                 struct spart* restrict sparts_written,
+                                const int subsample,
+                                const float subsample_ratio, const int snap_num,
                                 const size_t Nsparts,
                                 const size_t Nsparts_written);
 void io_collect_bparts_to_write(const struct bpart* restrict bparts,
                                 struct bpart* restrict bparts_written,
+                                const int subsample,
+                                const float subsample_ratio, const int snap_num,
                                 const size_t Nbparts,
                                 const size_t Nbparts_written);
 void io_collect_gparts_to_write(const struct gpart* restrict gparts,
                                 const struct velociraptor_gpart_data* vr_data,
                                 struct gpart* restrict gparts_written,
                                 struct velociraptor_gpart_data* vr_data_written,
+                                const int subsample,
+                                const float subsample_ratio, const int snap_num,
                                 const size_t Ngparts,
                                 const size_t Ngparts_written, int with_stf);
 void io_collect_gparts_background_to_write(
     const struct gpart* restrict gparts,
     const struct velociraptor_gpart_data* vr_data,
     struct gpart* restrict gparts_written,
-    struct velociraptor_gpart_data* vr_data_written, const size_t Ngparts,
+    struct velociraptor_gpart_data* vr_data_written, const int subsample,
+    const float subsample_ratio, const int snap_num, const size_t Ngparts,
     const size_t Ngparts_written, int with_stf);
 void io_collect_gparts_neutrino_to_write(
     const struct gpart* restrict gparts,
     const struct velociraptor_gpart_data* vr_data,
     struct gpart* restrict gparts_written,
-    struct velociraptor_gpart_data* vr_data_written, const size_t Ngparts,
+    struct velociraptor_gpart_data* vr_data_written, const int subsample,
+    const float subsample_ratio, const int snap_num, const size_t Ngparts,
     const size_t Ngparts_written, int with_stf);
+
 void io_prepare_dm_gparts(struct threadpool* tp, struct gpart* const gparts,
                           size_t Ndm);
 void io_prepare_dm_background_gparts(struct threadpool* tp,
diff --git a/src/common_io_cells.c b/src/common_io_cells.c
index 6a3c945efca41e0df4eda964cec14bb3d37f85e1..1d76aea1708b529786e3b7b629c5620bf2337004 100644
--- a/src/common_io_cells.c
+++ b/src/common_io_cells.c
@@ -25,32 +25,38 @@
 
 /* Local includes. */
 #include "cell.h"
+#include "random.h"
 #include "timeline.h"
 #include "units.h"
 
-#if defined(HAVE_HDF5)
-
-#include <hdf5.h>
-
-/* MPI headers. */
-#ifdef WITH_MPI
-#include <mpi.h>
-#endif
-
-static long long cell_count_non_inhibited_gas(const struct cell* c) {
+static long long cell_count_non_inhibited_gas(const struct cell* c,
+                                              const int subsample,
+                                              const float subsample_ratio,
+                                              const int snap_num) {
   const int total_count = c->hydro.count;
   struct part* parts = c->hydro.parts;
   long long count = 0;
   for (int i = 0; i < total_count; ++i) {
     if ((parts[i].time_bin != time_bin_inhibited) &&
         (parts[i].time_bin != time_bin_not_created)) {
+
+      /* When subsampling, select particles at random */
+      if (subsample) {
+        const float r = random_unit_interval(parts[i].id, snap_num,
+                                             random_number_snapshot_sampling);
+        if (r > subsample_ratio) continue;
+      }
+
       ++count;
     }
   }
   return count;
 }
 
-static long long cell_count_non_inhibited_dark_matter(const struct cell* c) {
+static long long cell_count_non_inhibited_dark_matter(
+    const struct cell* c, const int subsample, const float subsample_ratio,
+    const int snap_num) {
+
   const int total_count = c->grav.count;
   struct gpart* gparts = c->grav.parts;
   long long count = 0;
@@ -58,6 +64,14 @@ static long long cell_count_non_inhibited_dark_matter(const struct cell* c) {
     if ((gparts[i].time_bin != time_bin_inhibited) &&
         (gparts[i].time_bin != time_bin_not_created) &&
         (gparts[i].type == swift_type_dark_matter)) {
+
+      /* When subsampling, select particles at random */
+      if (subsample) {
+        const float r =
+            random_unit_interval(gparts[i].id_or_neg_offset, snap_num,
+                                 random_number_snapshot_sampling);
+        if (r > subsample_ratio) continue;
+      }
       ++count;
     }
   }
@@ -65,7 +79,9 @@ static long long cell_count_non_inhibited_dark_matter(const struct cell* c) {
 }
 
 static long long cell_count_non_inhibited_background_dark_matter(
-    const struct cell* c) {
+    const struct cell* c, const int subsample, const float subsample_ratio,
+    const int snap_num) {
+
   const int total_count = c->grav.count;
   struct gpart* gparts = c->grav.parts;
   long long count = 0;
@@ -73,52 +89,97 @@ static long long cell_count_non_inhibited_background_dark_matter(
     if ((gparts[i].time_bin != time_bin_inhibited) &&
         (gparts[i].time_bin != time_bin_not_created) &&
         (gparts[i].type == swift_type_dark_matter_background)) {
+
+      /* When subsampling, select particles at random */
+      if (subsample) {
+        const float r =
+            random_unit_interval(gparts[i].id_or_neg_offset, snap_num,
+                                 random_number_snapshot_sampling);
+        if (r > subsample_ratio) continue;
+      }
+
       ++count;
     }
   }
   return count;
 }
 
-static long long cell_count_non_inhibited_stars(const struct cell* c) {
+static long long cell_count_non_inhibited_stars(const struct cell* c,
+                                                const int subsample,
+                                                const float subsample_ratio,
+                                                const int snap_num) {
   const int total_count = c->stars.count;
   struct spart* sparts = c->stars.parts;
   long long count = 0;
   for (int i = 0; i < total_count; ++i) {
     if ((sparts[i].time_bin != time_bin_inhibited) &&
         (sparts[i].time_bin != time_bin_not_created)) {
+
+      /* When subsampling, select particles at random */
+      if (subsample) {
+        const float r = random_unit_interval(sparts[i].id, snap_num,
+                                             random_number_snapshot_sampling);
+        if (r > subsample_ratio) continue;
+      }
+
       ++count;
     }
   }
   return count;
 }
 
-static long long cell_count_non_inhibited_black_holes(const struct cell* c) {
+static long long cell_count_non_inhibited_black_holes(
+    const struct cell* c, const int subsample, const float subsample_ratio,
+    const int snap_num) {
   const int total_count = c->black_holes.count;
   struct bpart* bparts = c->black_holes.parts;
   long long count = 0;
   for (int i = 0; i < total_count; ++i) {
     if ((bparts[i].time_bin != time_bin_inhibited) &&
         (bparts[i].time_bin != time_bin_not_created)) {
+
+      /* When subsampling, select particles at random */
+      if (subsample) {
+        const float r = random_unit_interval(bparts[i].id, snap_num,
+                                             random_number_snapshot_sampling);
+        if (r > subsample_ratio) continue;
+      }
+
       ++count;
     }
   }
   return count;
 }
 
-static long long cell_count_non_inhibited_sinks(const struct cell* c) {
+static long long cell_count_non_inhibited_sinks(const struct cell* c,
+                                                const int subsample,
+                                                const float subsample_ratio,
+                                                const int snap_num) {
   const int total_count = c->sinks.count;
   struct sink* sinks = c->sinks.parts;
   long long count = 0;
   for (int i = 0; i < total_count; ++i) {
     if ((sinks[i].time_bin != time_bin_inhibited) &&
         (sinks[i].time_bin != time_bin_not_created)) {
+
+      /* When subsampling, select particles at random */
+      if (subsample) {
+        const float r = random_unit_interval(sinks[i].id, snap_num,
+                                             random_number_snapshot_sampling);
+        if (r > subsample_ratio) continue;
+      }
+
       ++count;
     }
   }
   return count;
 }
 
-static long long cell_count_non_inhibited_neutrinos(const struct cell* c) {
+static long long cell_count_non_inhibited_neutrinos(const struct cell* c,
+                                                    const int subsample,
+                                                    const float subsample_ratio,
+                                                    const int snap_num) {
+
   const int total_count = c->grav.count;
   struct gpart* gparts = c->grav.parts;
   long long count = 0;
@@ -126,12 +187,205 @@ static long long cell_count_non_inhibited_neutrinos(const struct cell* c) {
     if ((gparts[i].time_bin != time_bin_inhibited) &&
         (gparts[i].time_bin != time_bin_not_created) &&
         (gparts[i].type == swift_type_neutrino)) {
+
+      /* When subsampling, select particles at random */
+      if (subsample) {
+        const float r =
+            random_unit_interval(gparts[i].id_or_neg_offset, snap_num,
+                                 random_number_snapshot_sampling);
+        if (r > subsample_ratio) continue;
+      }
+
       ++count;
     }
   }
   return count;
 }
 
+/**
+ * @brief Count the number of local non-inhibited gas particles to write.
+ *
+ * Takes into account downsampling.
+ *
+ * @param s The #space.
+ * @param subsample Are we subsampling?
+ * @param subsample_ratio The fraction of particle to keep when subsampling.
+ * @param snap_num The snapshot number to use as random seed.
+ */
+long long io_count_gas_to_write(const struct space* s, const int subsample,
+                                const float subsample_ratio,
+                                const int snap_num) {
+
+  long long count = 0;
+  for (int i = 0; i < s->nr_local_cells; ++i) {
+
+    const struct cell* c = &s->cells_top[s->local_cells_top[i]];
+    count +=
+        cell_count_non_inhibited_gas(c, subsample, subsample_ratio, snap_num);
+  }
+  return count;
+}
+
+/**
+ * @brief Count the number of local non-inhibited dark matter particles to
+ * write.
+ *
+ * Takes into account downsampling.
+ *
+ * @param s The #space.
+ * @param subsample Are we subsampling?
+ * @param subsample_ratio The fraction of particle to keep when subsampling.
+ * @param snap_num The snapshot number to use as random seed.
+ */
+long long io_count_dark_matter_to_write(const struct space* s,
+                                        const int subsample,
+                                        const float subsample_ratio,
+                                        const int snap_num) {
+
+  long long count = 0;
+  for (int i = 0; i < s->nr_local_cells; ++i) {
+
+    const struct cell* c = &s->cells_top[s->local_cells_top[i]];
+    count += cell_count_non_inhibited_dark_matter(c, subsample, subsample_ratio,
+                                                  snap_num);
+  }
+  return count;
+}
+
+/**
+ * @brief Count the number of local non-inhibited background dark matter
+ * particles to write.
+ *
+ * Takes into account downsampling.
+ *
+ * @param s The #space.
+ * @param subsample Are we subsampling?
+ * @param subsample_ratio The fraction of particle to keep when subsampling.
+ * @param snap_num The snapshot number to use as random seed.
+ */
+long long io_count_background_dark_matter_to_write(const struct space* s,
+                                                   const int subsample,
+                                                   const float subsample_ratio,
+                                                   const int snap_num) {
+
+  long long count = 0;
+  for (int i = 0; i < s->nr_local_cells; ++i) {
+
+    const struct cell* c = &s->cells_top[s->local_cells_top[i]];
+    count += cell_count_non_inhibited_background_dark_matter(
+        c, subsample, subsample_ratio, snap_num);
+  }
+  return count;
+}
+
+/**
+ * @brief Count the number of local non-inhibited stars particles to write.
+ *
+ * Takes into account downsampling.
+ *
+ * @param s The #space.
+ * @param subsample Are we subsampling?
+ * @param subsample_ratio The fraction of particle to keep when subsampling.
+ * @param snap_num The snapshot number to use as random seed.
+ */
+long long io_count_stars_to_write(const struct space* s, const int subsample,
+                                  const float subsample_ratio,
+                                  const int snap_num) {
+
+  long long count = 0;
+  for (int i = 0; i < s->nr_local_cells; ++i) {
+
+    const struct cell* c = &s->cells_top[s->local_cells_top[i]];
+    count +=
+        cell_count_non_inhibited_stars(c, subsample, subsample_ratio, snap_num);
+  }
+  return count;
+}
+
+/**
+ * @brief Count the number of local non-inhibited sinks particles to write.
+ *
+ * Takes into account downsampling.
+ *
+ * @param s The #space.
+ * @param subsample Are we subsampling?
+ * @param subsample_ratio The fraction of particle to keep when subsampling.
+ * @param snap_num The snapshot number to use as random seed.
+ */
+long long io_count_sinks_to_write(const struct space* s, const int subsample,
+                                  const float subsample_ratio,
+                                  const int snap_num) {
+
+  long long count = 0;
+  for (int i = 0; i < s->nr_local_cells; ++i) {
+
+    const struct cell* c = &s->cells_top[s->local_cells_top[i]];
+    count +=
+        cell_count_non_inhibited_sinks(c, subsample, subsample_ratio, snap_num);
+  }
+  return count;
+}
+
+/**
+ * @brief Count the number of local non-inhibited black holes particles to
+ * write.
+ *
+ * Takes into account downsampling.
+ *
+ * @param s The #space.
+ * @param subsample Are we subsampling?
+ * @param subsample_ratio The fraction of particle to keep when subsampling.
+ * @param snap_num The snapshot number to use as random seed.
+ */
+long long io_count_black_holes_to_write(const struct space* s,
+                                        const int subsample,
+                                        const float subsample_ratio,
+                                        const int snap_num) {
+
+  long long count = 0;
+  for (int i = 0; i < s->nr_local_cells; ++i) {
+
+    const struct cell* c = &s->cells_top[s->local_cells_top[i]];
+    count += cell_count_non_inhibited_black_holes(c, subsample, subsample_ratio,
+                                                  snap_num);
+  }
+  return count;
+}
+
+/**
+ * @brief Count the number of local non-inhibited neutrinos particles to write.
+ *
+ * Takes into account downsampling.
+ *
+ * @param s The #space.
+ * @param subsample Are we subsampling?
+ * @param subsample_ratio The fraction of particle to keep when subsampling.
+ * @param snap_num The snapshot number to use as random seed.
+ */
+long long io_count_neutrinos_to_write(const struct space* s,
+                                      const int subsample,
+                                      const float subsample_ratio,
+                                      const int snap_num) {
+
+  long long count = 0;
+  for (int i = 0; i < s->nr_local_cells; ++i) {
+
+    const struct cell* c = &s->cells_top[s->local_cells_top[i]];
+    count += cell_count_non_inhibited_neutrinos(c, subsample, subsample_ratio,
+                                                snap_num);
+  }
+  return count;
+}
+
+#if defined(HAVE_HDF5)
+
+#include <hdf5.h>
+
+/* MPI headers. */
+#ifdef WITH_MPI
+#include <mpi.h>
+#endif
+
 /**
  * @brief Write a single 1D array to a hdf5 group.
  *
@@ -190,10 +444,32 @@ void io_write_array(hid_t h_grp, const int n, const void* array,
   H5Sclose(h_space);
 }
 
+/**
+ * @brief Compute and write the top-level cell counts and offsets meta-data.
+ *
+ * @param h_grp the hdf5 group to write to.
+ * @param cdim The number of top-level cells along each axis.
+ * @param dim The box size.
+ * @param cells_top The top-level cells.
+ * @param nr_cells The number of top-level cells.
+ * @param distributed Is this a distributed snapshot?
+ * @param subsample Are we subsampling the different particle types?
+ * @param subsample_fraction The fraction of particles to keep when subsampling.
+ * @param snap_num The snapshot number used as subsampling random seed.
+ * @param global_counts The total number of particles across all nodes.
+ * @param global_offsets The offsets of this node into the global list of
+ * particles.
+ * @param numFields The number of fields to write for each particle type.
+ * @param internal_units The internal unit system.
+ * @param snapshot_units The snapshot unit system.
+ */
 void io_write_cell_offsets(hid_t h_grp, const int cdim[3], const double dim[3],
                            const struct cell* cells_top, const int nr_cells,
                            const double width[3], const int nodeID,
                            const int distributed,
+                           const int subsample[swift_type_count],
+                           const float subsample_fraction[swift_type_count],
+                           const int snap_num,
                            const long long global_counts[swift_type_count],
                            const long long global_offsets[swift_type_count],
                            const int num_fields[swift_type_count],
@@ -288,14 +564,34 @@ void io_write_cell_offsets(hid_t h_grp, const int cdim[3], const double dim[3],
       centres[i * 3 + 2] = box_wrap(centres[i * 3 + 2], 0.0, dim[2]);
 
       /* Count real particles that will be written */
-      count_part[i] = cell_count_non_inhibited_gas(&cells_top[i]);
-      count_gpart[i] = cell_count_non_inhibited_dark_matter(&cells_top[i]);
+      count_part[i] = cell_count_non_inhibited_gas(
+          &cells_top[i], subsample[swift_type_gas],
+          subsample_fraction[swift_type_gas], snap_num);
+
+      count_gpart[i] = cell_count_non_inhibited_dark_matter(
+          &cells_top[i], subsample[swift_type_dark_matter],
+          subsample_fraction[swift_type_dark_matter], snap_num);
+
       count_background_gpart[i] =
-          cell_count_non_inhibited_background_dark_matter(&cells_top[i]);
-      count_spart[i] = cell_count_non_inhibited_stars(&cells_top[i]);
-      count_bpart[i] = cell_count_non_inhibited_black_holes(&cells_top[i]);
-      count_sink[i] = cell_count_non_inhibited_sinks(&cells_top[i]);
-      count_nupart[i] = cell_count_non_inhibited_neutrinos(&cells_top[i]);
+          cell_count_non_inhibited_background_dark_matter(
+              &cells_top[i], subsample[swift_type_dark_matter_background],
+              subsample_fraction[swift_type_dark_matter_background], snap_num);
+
+      count_spart[i] = cell_count_non_inhibited_stars(
+          &cells_top[i], subsample[swift_type_stars],
+          subsample_fraction[swift_type_stars], snap_num);
+
+      count_bpart[i] = cell_count_non_inhibited_black_holes(
+          &cells_top[i], subsample[swift_type_black_hole],
+          subsample_fraction[swift_type_black_hole], snap_num);
+
+      count_sink[i] = cell_count_non_inhibited_sinks(
+          &cells_top[i], subsample[swift_type_sink],
+          subsample_fraction[swift_type_sink], snap_num);
+
+      count_nupart[i] = cell_count_non_inhibited_neutrinos(
+          &cells_top[i], subsample[swift_type_neutrino],
+          subsample_fraction[swift_type_neutrino], snap_num);
 
       /* Offsets including the global offset of all particles on this MPI rank
        * Note that in the distributed case, the global offsets are 0 such that
diff --git a/src/common_io_fields.c b/src/common_io_fields.c
index 00fa41aaf6c0f6fc824cbd7feea2c1a9778133e7..fe29638824d6e693bc4ca4e0d3bd1a10b8ce848a 100644
--- a/src/common_io_fields.c
+++ b/src/common_io_fields.c
@@ -168,6 +168,14 @@ void io_prepare_output_fields(struct output_options* output_options,
     /* Default snapshot subdir name for this output selection */
     char subdir_name[FILENAME_BUFFER_SIZE] = select_output_default_subdir_name;
 
+    int subsample[swift_type_count];
+    for (int i = 0; i < swift_type_count; ++i)
+      subsample[i] = select_output_default_subsample;
+
+    float subsample_fraction[swift_type_count];
+    for (int i = 0; i < swift_type_count; ++i)
+      subsample_fraction[i] = select_output_default_subsample_fraction;
+
     /* Initialise section-specific writing counters for each particle type.
      * If default is 'write', then we start from the total to deduct any fields
      * that are switched off. If the default is 'off', we have to start from
@@ -220,6 +228,20 @@ void io_prepare_output_fields(struct output_options* output_options,
         continue;
       }
 
+      /* Deal with a possible non-standard subsampling option */
+      if (strstr(param_name, ":subsample") != NULL) {
+        parser_get_param_int_array(params, param_name, swift_type_count,
+                                   subsample);
+        continue;
+      }
+
+      /* Deal with a possible non-standard subsampling fraction */
+      if (strstr(param_name, ":subsample_fraction") != NULL) {
+        parser_get_param_float_array(params, param_name, swift_type_count,
+                                     subsample_fraction);
+        continue;
+      }
+
       /* Get the particle type for current parameter
        * (raises an error if it could not determine it) */
       const int param_ptype = io_get_param_ptype(param_name);
@@ -300,6 +322,10 @@ void io_prepare_output_fields(struct output_options* output_options,
     /* Also copy the output names */
     strcpy(output_options->basenames[section_id], basename);
     strcpy(output_options->subdir_names[section_id], subdir_name);
+    memcpy(output_options->subsample[section_id], subsample,
+           swift_type_count * sizeof(int));
+    memcpy(output_options->subsample_fractions[section_id], subsample_fraction,
+           swift_type_count * sizeof(float));
 
   } /* Ends loop over sections, for different output classes */
 
@@ -314,6 +340,12 @@ void io_prepare_output_fields(struct output_options* output_options,
             select_output_default_basename);
     sprintf(output_options->subdir_names[default_id], "%s",
             select_output_default_subdir_name);
+    for (int i = 0; i < swift_type_count; ++i)
+      output_options->subsample[default_id][i] =
+          select_output_default_subsample;
+    for (int i = 0; i < swift_type_count; ++i)
+      output_options->subsample_fractions[default_id][i] =
+          select_output_default_subsample_fraction;
   }
 }
 
diff --git a/src/distributed_io.c b/src/distributed_io.c
index ff7944beee922b78ce014da5d6033a54897b7453..2f5bc125e5e9e2405c94a1d80f67352df7fced68 100644
--- a/src/distributed_io.c
+++ b/src/distributed_io.c
@@ -321,34 +321,6 @@ void write_output_distributed(struct engine* e,
   const size_t Nstars = e->s->nr_sparts;
   const size_t Nblackholes = e->s->nr_bparts;
 
-  size_t Ndm_background = 0;
-  if (with_DM_background) {
-    Ndm_background = io_count_dm_background_gparts(gparts, Ntot);
-  }
-  size_t Ndm_neutrino = 0;
-  if (with_neutrinos) {
-    Ndm_neutrino = io_count_dm_neutrino_gparts(gparts, Ntot);
-  }
-
-  /* Number of particles that we will write in this file.
-   * Recall that background particles are never inhibited and have no extras */
-  const size_t Ntot_written =
-      e->s->nr_gparts - e->s->nr_inhibited_gparts - e->s->nr_extra_gparts;
-  const size_t Ngas_written =
-      e->s->nr_parts - e->s->nr_inhibited_parts - e->s->nr_extra_parts;
-  const size_t Nsinks_written =
-      e->s->nr_sinks - e->s->nr_inhibited_sinks - e->s->nr_extra_sinks;
-  const size_t Nstars_written =
-      e->s->nr_sparts - e->s->nr_inhibited_sparts - e->s->nr_extra_sparts;
-  const size_t Nblackholes_written =
-      e->s->nr_bparts - e->s->nr_inhibited_bparts - e->s->nr_extra_bparts;
-  const size_t Nbaryons_written =
-      Ngas_written + Nstars_written + Nblackholes_written + Nsinks_written;
-  const size_t Ndm_written =
-      Ntot_written > 0
-          ? Ntot_written - Nbaryons_written - Ndm_background - Ndm_neutrino
-          : 0;
-
   /* Determine if we are writing a reduced snapshot, and if so which
    * output selection type to use */
   char current_selection_name[FIELD_BUFFER_SIZE] =
@@ -404,6 +376,86 @@ void write_output_distributed(struct engine* e,
   if (mpi_rank == 0) safe_checkdir(dirName, /*create=*/1);
   MPI_Barrier(comm);
 
+  /* Do we want to sub-sample any of the arrays */
+  int subsample[swift_type_count];
+  float subsample_fraction[swift_type_count];
+  for (int i = 0; i < swift_type_count; ++i) {
+    subsample[i] = 0;
+    subsample_fraction[i] = 1.f;
+  }
+
+  output_options_get_subsampling(
+      output_options, current_selection_name, e->snapshot_subsample,
+      e->snapshot_subsample_fraction, subsample, subsample_fraction);
+
+  /* Is any particle type being subsampled? */
+  int subsample_any = 0;
+  for (int i = 0; i < swift_type_count; ++i) {
+    subsample_any += subsample[i];
+    if (!subsample[i]) subsample_fraction[i] = 1.f;
+  }
+
+  /* Number of particles that we will write */
+  size_t Ngas_written, Ndm_written, Ndm_background, Ndm_neutrino,
+      Nsinks_written, Nstars_written, Nblackholes_written;
+
+  if (subsample[swift_type_gas]) {
+    Ngas_written = io_count_gas_to_write(e->s, /*subsample=*/1,
+                                         subsample_fraction[swift_type_gas],
+                                         e->snapshot_output_count);
+  } else {
+    Ngas_written =
+        e->s->nr_parts - e->s->nr_inhibited_parts - e->s->nr_extra_parts;
+  }
+
+  if (subsample[swift_type_stars]) {
+    Nstars_written = io_count_stars_to_write(
+        e->s, /*subsample=*/1, subsample_fraction[swift_type_stars],
+        e->snapshot_output_count);
+  } else {
+    Nstars_written =
+        e->s->nr_sparts - e->s->nr_inhibited_sparts - e->s->nr_extra_sparts;
+  }
+
+  if (subsample[swift_type_black_hole]) {
+    Nblackholes_written = io_count_black_holes_to_write(
+        e->s, /*subsample=*/1, subsample_fraction[swift_type_black_hole],
+        e->snapshot_output_count);
+  } else {
+    Nblackholes_written =
+        e->s->nr_bparts - e->s->nr_inhibited_bparts - e->s->nr_extra_bparts;
+  }
+
+  if (subsample[swift_type_sink]) {
+    Nsinks_written = io_count_sinks_to_write(
+        e->s, /*subsample=*/1, subsample_fraction[swift_type_sink],
+        e->snapshot_output_count);
+  } else {
+    Nsinks_written =
+        e->s->nr_sinks - e->s->nr_inhibited_sinks - e->s->nr_extra_sinks;
+  }
+
+  Ndm_written = io_count_dark_matter_to_write(
+      e->s, subsample[swift_type_dark_matter],
+      subsample_fraction[swift_type_dark_matter], e->snapshot_output_count);
+
+  if (with_DM_background) {
+    Ndm_background = io_count_dark_matter_to_write(
+        e->s, subsample[swift_type_dark_matter_background],
+        subsample_fraction[swift_type_dark_matter_background],
+        e->snapshot_output_count);
+  } else {
+    Ndm_background = 0;
+  }
+
+  if (with_neutrinos) {
+    Ndm_neutrino = io_count_neutrinos_to_write(
+        e->s, subsample[swift_type_neutrino],
+        subsample_fraction[swift_type_neutrino], e->snapshot_output_count);
+  } else {
+    Ndm_neutrino = 0;
+  }
+
   /* Compute offset in the file and total number of particles */
   const long long N[swift_type_count] = {
       Ngas_written,   Ndm_written,         Ndm_background, Nsinks_written,
@@ -467,7 +519,7 @@ void write_output_distributed(struct engine* e,
   struct tm* timeinfo = localtime(&tm);
   char snapshot_date[64];
   strftime(snapshot_date, 64, "%T %F %Z", timeinfo);
-  io_write_attribute_s(h_grp, "Snapshot date", snapshot_date);
+  io_write_attribute_s(h_grp, "SnapshotDate", snapshot_date);
 
   /* GADGET-2 legacy values:  Number of particles of each type */
   long long numParticlesThisFile[swift_type_count] = {0};
@@ -507,9 +559,16 @@ void write_output_distributed(struct engine* e,
                      swift_type_count);
   io_write_attribute_i(h_grp, "NumFilesPerSnapshot", numFiles);
   io_write_attribute_i(h_grp, "ThisFile", mpi_rank);
-  io_write_attribute_s(h_grp, "OutputType", "FullVolume");
   io_write_attribute_s(h_grp, "SelectOutput", current_selection_name);
 
+  if (subsample_any) {
+    io_write_attribute_s(h_grp, "OutputType", "SubSampled");
+    io_write_attribute(h_grp, "SubSampleFractions", FLOAT, subsample_fraction,
+                       swift_type_count);
+  } else {
+    io_write_attribute_s(h_grp, "OutputType", "FullVolume");
+  }
+
   /* Close header */
   H5Gclose(h_grp);
 
@@ -527,8 +586,9 @@ void write_output_distributed(struct engine* e,
   /* Write the location of the particles in the arrays */
   io_write_cell_offsets(h_grp, e->s->cdim, e->s->dim, e->s->cells_top,
                         e->s->nr_cells, e->s->width, mpi_rank,
-                        /*distributed=*/1, N_total, global_offsets, numFields,
-                        internal_units, snapshot_units);
+                        /*distributed=*/1, subsample, subsample_fraction,
+                        e->snapshot_output_count, N_total, global_offsets,
+                        numFields, internal_units, snapshot_units);
   H5Gclose(h_grp);
 
   /* Loop over all particle types */
@@ -599,8 +659,10 @@ void write_output_distributed(struct engine* e,
             error("Error while allocating temporary memory for xparts");
 
           /* Collect the particles we want to write */
-          io_collect_parts_to_write(parts, xparts, parts_written,
-                                    xparts_written, Ngas, Ngas_written);
+          io_collect_parts_to_write(
+              parts, xparts, parts_written, xparts_written,
+              subsample[swift_type_gas], subsample_fraction[swift_type_gas],
+              e->snapshot_output_count, Ngas, Ngas_written);
 
           /* Select the fields to write */
           io_select_hydro_fields(parts_written, xparts_written, with_cosmology,
@@ -642,9 +704,11 @@ void write_output_distributed(struct engine* e,
           }
 
           /* Collect the non-inhibited DM particles from gpart */
-          io_collect_gparts_to_write(gparts, e->s->gpart_group_data,
-                                     gparts_written, gpart_group_data_written,
-                                     Ntot, Ndm_written, with_stf);
+          io_collect_gparts_to_write(
+              gparts, e->s->gpart_group_data, gparts_written,
+              gpart_group_data_written, subsample[swift_type_dark_matter],
+              subsample_fraction[swift_type_dark_matter],
+              e->snapshot_output_count, Ntot, Ndm_written, with_stf);
 
           /* Select the fields to write */
           io_select_dm_fields(gparts_written, gpart_group_data_written,
@@ -676,7 +740,10 @@ void write_output_distributed(struct engine* e,
         /* Collect the non-inhibited DM particles from gpart */
         io_collect_gparts_background_to_write(
             gparts, e->s->gpart_group_data, gparts_written,
-            gpart_group_data_written, Ntot, Ndm_background, with_stf);
+            gpart_group_data_written,
+            subsample[swift_type_dark_matter_background],
+            subsample_fraction[swift_type_dark_matter_background],
+            e->snapshot_output_count, Ntot, Ndm_background, with_stf);
 
         /* Select the fields to write */
         io_select_dm_fields(gparts_written, gpart_group_data_written, with_fof,
@@ -707,7 +774,9 @@ void write_output_distributed(struct engine* e,
         /* Collect the non-inhibited DM particles from gpart */
         io_collect_gparts_neutrino_to_write(
             gparts, e->s->gpart_group_data, gparts_written,
-            gpart_group_data_written, Ntot, Ndm_neutrino, with_stf);
+            gpart_group_data_written, subsample[swift_type_neutrino],
+            subsample_fraction[swift_type_neutrino], e->snapshot_output_count,
+            Ntot, Ndm_neutrino, with_stf);
 
         /* Select the fields to write */
         io_select_neutrino_fields(gparts_written, gpart_group_data_written,
@@ -736,8 +805,10 @@ void write_output_distributed(struct engine* e,
             error("Error while allocating temporary memory for sinks");
 
           /* Collect the particles we want to write */
-          io_collect_sinks_to_write(sinks, sinks_written, Nsinks,
-                                    Nsinks_written);
+          io_collect_sinks_to_write(
+              sinks, sinks_written, subsample[swift_type_sink],
+              subsample_fraction[swift_type_sink], e->snapshot_output_count,
+              Nsinks, Nsinks_written);
 
           /* Select the fields to write */
           io_select_sink_fields(sinks_written, with_cosmology, with_fof,
@@ -767,8 +838,10 @@ void write_output_distributed(struct engine* e,
             error("Error while allocating temporary memory for sparts");
 
           /* Collect the particles we want to write */
-          io_collect_sparts_to_write(sparts, sparts_written, Nstars,
-                                     Nstars_written);
+          io_collect_sparts_to_write(
+              sparts, sparts_written, subsample[swift_type_stars],
+              subsample_fraction[swift_type_stars], e->snapshot_output_count,
+              Nstars, Nstars_written);
 
           /* Select the fields to write */
           io_select_star_fields(sparts_written, with_cosmology, with_fof,
@@ -798,8 +871,10 @@ void write_output_distributed(struct engine* e,
             error("Error while allocating temporary memory for bparts");
 
           /* Collect the particles we want to write */
-          io_collect_bparts_to_write(bparts, bparts_written, Nblackholes,
-                                     Nblackholes_written);
+          io_collect_bparts_to_write(
+              bparts, bparts_written, subsample[swift_type_black_hole],
+              subsample_fraction[swift_type_black_hole],
+              e->snapshot_output_count, Nblackholes, Nblackholes_written);
 
           /* Select the fields to write */
           io_select_bh_fields(bparts_written, with_cosmology, with_fof,
diff --git a/src/engine.c b/src/engine.c
index 896b0d0a856cfdf8418ed6e5203ba169c8cf9342..74c3a84ea16f3973b92c6f985fd0ce852b2b8446 100644
--- a/src/engine.c
+++ b/src/engine.c
@@ -127,6 +127,8 @@ const char *engine_policy_names[] = {"none",
                                      "sink",
                                      "rt"};
 
+const int engine_default_snapshot_subsample[swift_type_count] = {0};
+
 /** The rank of the engine as a global variable (for messages). */
 int engine_rank;
 
@@ -2825,6 +2827,11 @@ void engine_init(
   parser_get_param_string(params, "Snapshots:basename", e->snapshot_base_name);
   parser_get_opt_param_string(params, "Snapshots:subdir", e->snapshot_subdir,
                               engine_default_snapshot_subdir);
+  parser_get_opt_param_int_array(params, "Snapshots:subsample",
+                                 swift_type_count, e->snapshot_subsample);
+  parser_get_opt_param_float_array(params, "Snapshots:subsample_fraction",
+                                   swift_type_count,
+                                   e->snapshot_subsample_fraction);
   e->snapshot_run_on_dump =
       parser_get_opt_param_int(params, "Snapshots:run_on_dump", 0);
   if (e->snapshot_run_on_dump) {
diff --git a/src/engine.h b/src/engine.h
index c69e3788ffe3b4b0973b1c33731f2c1acd808197..1ef13fa954678875f69d4987e41cdd0366caa213 100644
--- a/src/engine.h
+++ b/src/engine.h
@@ -315,6 +315,8 @@ struct engine {
   char snapshot_base_name[PARSER_MAX_LINE_SIZE];
   char snapshot_subdir[PARSER_MAX_LINE_SIZE];
   char snapshot_dump_command[PARSER_MAX_LINE_SIZE];
+  int snapshot_subsample[swift_type_count];
+  float snapshot_subsample_fraction[swift_type_count];
   int snapshot_run_on_dump;
   int snapshot_distributed;
   int snapshot_compression;
diff --git a/src/fof_catalogue_io.c b/src/fof_catalogue_io.c
index b0ebeb8ae494376f04c63fa29fbba35d7b6c1be0..151345970c8d8c8b776ea70e6d88125d13327153 100644
--- a/src/fof_catalogue_io.c
+++ b/src/fof_catalogue_io.c
@@ -82,7 +82,7 @@ void write_fof_hdf5_header(hid_t h_file, const struct engine* e,
   struct tm* timeinfo = localtime(&tm);
   char snapshot_date[64];
   strftime(snapshot_date, 64, "%T %F %Z", timeinfo);
-  io_write_attribute_s(h_grp, "Snapshot date", snapshot_date);
+  io_write_attribute_s(h_grp, "SnapshotDate", snapshot_date);
 
   /* GADGET-2 legacy values */
   /* Number of particles of each type */
diff --git a/src/line_of_sight.c b/src/line_of_sight.c
index ca0045d1099a90be1ae17f90e034fa8b738d68c0..49af137901f45bf1130e50b22a409322a542c543 100644
--- a/src/line_of_sight.c
+++ b/src/line_of_sight.c
@@ -504,7 +504,7 @@ void write_hdf5_header(hid_t h_file, const struct engine *e,
   struct tm *timeinfo = localtime(&tm);
   char snapshot_date[64];
   strftime(snapshot_date, 64, "%T %F %Z", timeinfo);
-  io_write_attribute_s(h_grp, "Snapshot date", snapshot_date);
+  io_write_attribute_s(h_grp, "SnapshotDate", snapshot_date);
 
   /* GADGET-2 legacy values */
   /* Number of particles of each type */
diff --git a/src/output_options.c b/src/output_options.c
index 28c8ce32302223a31122405e1c7afca4809f0abb..214dc9d101f551202350ede6b21a623914353b27 100644
--- a/src/output_options.c
+++ b/src/output_options.c
@@ -114,6 +114,14 @@ void output_options_struct_dump(struct output_options* output_options,
   restart_write_blocks(output_options->subdir_names, count_chars * sizeof(char),
                        1, stream, "output_options_subdir_names",
                        "output options");
+
+  const size_t count_sub =
+      (OUTPUT_LIST_MAX_NUM_OF_SELECT_OUTPUT_STYLES + 1) * swift_type_count;
+  restart_write_blocks(output_options->subsample, count_sub * sizeof(int), 1,
+                       stream, "output_options_subsample", "output options");
+  restart_write_blocks(output_options->subsample_fractions,
+                       count_sub * sizeof(float), 1, stream,
+                       "output_options_subsample_fractions", "output options");
 }
 
 /**
@@ -140,9 +148,17 @@ void output_options_struct_restore(struct output_options* output_options,
   const size_t count_chars =
       (OUTPUT_LIST_MAX_NUM_OF_SELECT_OUTPUT_STYLES + 1) * FILENAME_BUFFER_SIZE;
   restart_read_blocks(output_options->basenames, count_chars * sizeof(char), 1,
-                      stream, NULL, "output options_basenames");
+                      stream, NULL, "output_options_basenames");
   restart_read_blocks(output_options->subdir_names, count_chars * sizeof(char),
-                      1, stream, NULL, "output options_subdir_names");
+                      1, stream, NULL, "output_options_subdir_names");
+
+  const size_t count_sub =
+      (OUTPUT_LIST_MAX_NUM_OF_SELECT_OUTPUT_STYLES + 1) * swift_type_count;
+  restart_read_blocks(output_options->subsample, count_sub * sizeof(int), 1,
+                      stream, NULL, "output_options_subsample");
+  restart_read_blocks(output_options->subsample_fractions,
+                      count_sub * sizeof(float), 1, stream, NULL,
+                      "output_options_subsample_fractions");
 }
 
 /**
@@ -309,7 +325,7 @@ int output_options_get_num_fields_to_write(
  * @param default_subdirname The default general sub-directory name.
  * @param default_basename The default general snapshot base name.
  * @param subdir_name (return) The sub-directory name to use for this dump.
- * @param basename (return) The snapshot base name to use for this dump,
+ * @param basename (return) The snapshot base name to use for this dump.
  */
 void output_options_get_basename(const struct output_options* output_options,
                                  const char* selection_name,
@@ -347,3 +363,55 @@ void output_options_get_basename(const struct output_options* output_options,
     sprintf(subdir_name, "%s", output_options->subdir_names[selection_id]);
   }
 }
+
+/**
+ * @brief Return the subsampling fraction for the current output selection.
+ *
+ * @param output_options The #output_options structure.
+ * @param selection_name The current output selection name.
+ * @param default_subsample The default general subsampling status.
+ * @param default_subsample_fraction The default general subsampling fraction.
+ * @param subsample (return) The subsampling status to use for this dump.
+ * @param subsample_fraction (return) The subsampling fraction to use for this
+ * dump.
+ */
+void output_options_get_subsampling(
+    const struct output_options* output_options, const char* selection_name,
+    const int default_subsample[swift_type_count],
+    const float default_subsample_fraction[swift_type_count],
+    int subsample[swift_type_count],
+    float subsample_fraction[swift_type_count]) {
+
+  /* Get the ID of the output selection in the structure */
+  int selection_id =
+      parser_get_section_id(output_options->select_output, selection_name);
+
+  /* Special treatment for absent `Default` section */
+  if (selection_id < 0) {
+    selection_id = output_options->select_output->sectionCount;
+  }
+
+  /* If the default keyword is found, we use the subsample flag provided
+   * in the param file (not the select output!), aka. the argument
+   * of the function. */
+  if (output_options->subsample[selection_id][0] ==
+      select_output_default_subsample) {
+    memcpy(subsample, default_subsample, sizeof(int) * swift_type_count);
+  } else {
+    memcpy(subsample, output_options->subsample[selection_id],
+           sizeof(int) * swift_type_count);
+  }
+
+  /* If the default keyword is found, we use the subsample fraction provided
+   * in the param file (not the select output!), aka. the argument
+   * of the function. */
+  if (output_options->subsample_fractions[selection_id][0] ==
+      select_output_default_subsample_fraction) {
+    memcpy(subsample_fraction, default_subsample_fraction,
+           sizeof(float) * swift_type_count);
+  } else {
+    memcpy(subsample_fraction,
+           output_options->subsample_fractions[selection_id],
+           sizeof(float) * swift_type_count);
+  }
+}
diff --git a/src/output_options.h b/src/output_options.h
index d4ab0627d92b7b789d7398a0a3717614a4c86654..60cc9631b0b886427eeed530620506deec6819a0 100644
--- a/src/output_options.h
+++ b/src/output_options.h
@@ -32,6 +32,8 @@
 #define select_output_header_default_name "Default"
 #define select_output_default_basename "Standard"
 #define select_output_default_subdir_name "Standard"
+#define select_output_default_subsample -1
+#define select_output_default_subsample_fraction -1.f
 
 /**
  * @brief Output selection properties, including the parsed files.
@@ -58,6 +60,18 @@ struct output_options {
    * snapshot section of the param file. */
   char subdir_names[OUTPUT_LIST_MAX_NUM_OF_SELECT_OUTPUT_STYLES + 1]
                    [FILENAME_BUFFER_SIZE];
+
+  /*! Subsample choices to use for the snapshots of the given selection
+   * If set to the default, the code will use the name set in the
+   * snapshot section of the param file. */
+  int subsample[OUTPUT_LIST_MAX_NUM_OF_SELECT_OUTPUT_STYLES + 1]
+               [swift_type_count];
+
+  /*! Subsample fractions to use for the snapshots of the given selection
+   * If set to the default, the code will use the name set in the
+   * snapshot section of the param file. */
+  float subsample_fractions[OUTPUT_LIST_MAX_NUM_OF_SELECT_OUTPUT_STYLES + 1]
+                           [swift_type_count];
 };
 
 /* Create and destroy */
@@ -93,4 +107,10 @@ void output_options_get_basename(const struct output_options* output_options,
                                  char subdir_name[FILENAME_BUFFER_SIZE],
                                  char snap_basename[FILENAME_BUFFER_SIZE]);
 
+void output_options_get_subsampling(
+    const struct output_options* output_options, const char* selection_name,
+    const int default_subsample[swift_type_count],
+    const float default_subsample_fraction[swift_type_count],
+    int subsample[swift_type_count],
+    float subsample_fraction[swift_type_count]);
 #endif
diff --git a/src/parallel_io.c b/src/parallel_io.c
index 4dda323a4e5688b83cb10d8453b33d304c7cd366..1bdb23b3f7ad23bb63e490f6bacafc5078d8f1b6 100644
--- a/src/parallel_io.c
+++ b/src/parallel_io.c
@@ -1214,7 +1214,7 @@ void prepare_file(struct engine* e, const char* fileName,
   struct tm* timeinfo = localtime(&tm);
   char snapshot_date[64];
   strftime(snapshot_date, 64, "%T %F %Z", timeinfo);
-  io_write_attribute_s(h_grp, "Snapshot date", snapshot_date);
+  io_write_attribute_s(h_grp, "SnapshotDate", snapshot_date);
 
   /* GADGET-2 legacy values */
   /* Number of particles of each type */
@@ -1429,52 +1429,6 @@ void write_output_parallel(struct engine* e,
   const size_t Nstars = e->s->nr_sparts;
   const size_t Nsinks = e->s->nr_sinks;
   const size_t Nblackholes = e->s->nr_bparts;
-  // const size_t Nbaryons = Ngas + Nstars;
-  // const size_t Ndm = Ntot > 0 ? Ntot - Nbaryons : 0;
-
-  size_t Ndm_background = 0;
-  if (with_DM_background) {
-    Ndm_background = io_count_dm_background_gparts(gparts, Ntot);
-  }
-  size_t Ndm_neutrino = 0;
-  if (with_neutrinos) {
-    Ndm_neutrino = io_count_dm_neutrino_gparts(gparts, Ntot);
-  }
-
-  /* Number of particles that we will write */
-  const size_t Ntot_written =
-      e->s->nr_gparts - e->s->nr_inhibited_gparts - e->s->nr_extra_gparts;
-  const size_t Ngas_written =
-      e->s->nr_parts - e->s->nr_inhibited_parts - e->s->nr_extra_parts;
-  const size_t Nsinks_written =
-      e->s->nr_sinks - e->s->nr_inhibited_sinks - e->s->nr_extra_sinks;
-  const size_t Nstars_written =
-      e->s->nr_sparts - e->s->nr_inhibited_sparts - e->s->nr_extra_sparts;
-  const size_t Nblackholes_written =
-      e->s->nr_bparts - e->s->nr_inhibited_bparts - e->s->nr_extra_bparts;
-  const size_t Nbaryons_written =
-      Ngas_written + Nstars_written + Nblackholes_written + Nsinks_written;
-  const size_t Ndm_written =
-      Ntot_written > 0
-          ? Ntot_written - Nbaryons_written - Ndm_background - Ndm_neutrino
-          : 0;
-
-  /* Compute offset in the file and total number of particles */
-  size_t N[swift_type_count] = {
-      Ngas_written,   Ndm_written,         Ndm_background, Nsinks_written,
-      Nstars_written, Nblackholes_written, Ndm_neutrino};
-  long long N_total[swift_type_count] = {0};
-  long long offset[swift_type_count] = {0};
-  MPI_Exscan(N, offset, swift_type_count, MPI_LONG_LONG_INT, MPI_SUM, comm);
-  for (int ptype = 0; ptype < swift_type_count; ++ptype)
-    N_total[ptype] = offset[ptype] + N[ptype];
-
-  /* The last rank now has the correct N_total. Let's
-   * broadcast from there */
-  MPI_Bcast(N_total, swift_type_count, MPI_LONG_LONG_INT, mpi_size - 1, comm);
-
-  /* Now everybody konws its offset and the total number of
-   * particles of each type */
 
 #ifdef IO_SPEED_MEASUREMENT
   ticks tic = getticks();
@@ -1508,6 +1462,25 @@ void write_output_parallel(struct engine* e,
   /* Create the directory */
   if (mpi_rank == 0) safe_checkdir(snapshot_subdir_name, /*create=*/1);
 
+  /* Do we want to sub-sample any of the arrays */
+  int subsample[swift_type_count];
+  float subsample_fraction[swift_type_count];
+  for (int i = 0; i < swift_type_count; ++i) {
+    subsample[i] = 0;
+    subsample_fraction[i] = 1.f;
+  }
+
+  output_options_get_subsampling(
+      output_options, current_selection_name, e->snapshot_subsample,
+      e->snapshot_subsample_fraction, subsample, subsample_fraction);
+
+  /* Is any particle type being subsampled? */
+  int subsample_any = 0;
+  for (int i = 0; i < swift_type_count; ++i) {
+    subsample_any += subsample[i];
+    if (!subsample[i]) subsample_fraction[i] = 1.f;
+  }
+
   /* Total number of fields to write per ptype */
   int numFields[swift_type_count] = {0};
   for (int ptype = 0; ptype < swift_type_count; ++ptype) {
@@ -1515,6 +1488,84 @@ void write_output_parallel(struct engine* e,
         output_options, current_selection_name, ptype);
   }
 
+  /* Number of particles that we will write */
+  size_t Ngas_written, Ndm_written, Ndm_background, Ndm_neutrino,
+      Nsinks_written, Nstars_written, Nblackholes_written;
+
+  if (subsample[swift_type_gas]) {
+    Ngas_written = io_count_gas_to_write(e->s, /*subsample=*/1,
+                                         subsample_fraction[swift_type_gas],
+                                         e->snapshot_output_count);
+  } else {
+    Ngas_written =
+        e->s->nr_parts - e->s->nr_inhibited_parts - e->s->nr_extra_parts;
+  }
+
+  if (subsample[swift_type_stars]) {
+    Nstars_written = io_count_stars_to_write(
+        e->s, /*subsample=*/1, subsample_fraction[swift_type_stars],
+        e->snapshot_output_count);
+  } else {
+    Nstars_written =
+        e->s->nr_sparts - e->s->nr_inhibited_sparts - e->s->nr_extra_sparts;
+  }
+
+  if (subsample[swift_type_black_hole]) {
+    Nblackholes_written = io_count_black_holes_to_write(
+        e->s, /*subsample=*/1, subsample_fraction[swift_type_black_hole],
+        e->snapshot_output_count);
+  } else {
+    Nblackholes_written =
+        e->s->nr_bparts - e->s->nr_inhibited_bparts - e->s->nr_extra_bparts;
+  }
+
+  if (subsample[swift_type_sink]) {
+    Nsinks_written = io_count_sinks_to_write(
+        e->s, /*subsample=*/1, subsample_fraction[swift_type_sink],
+        e->snapshot_output_count);
+  } else {
+    Nsinks_written =
+        e->s->nr_sinks - e->s->nr_inhibited_sinks - e->s->nr_extra_sinks;
+  }
+
+  Ndm_written = io_count_dark_matter_to_write(
+      e->s, subsample[swift_type_dark_matter],
+      subsample_fraction[swift_type_dark_matter], e->snapshot_output_count);
+
+  if (with_DM_background) {
+    Ndm_background = io_count_dark_matter_to_write(
+        e->s, subsample[swift_type_dark_matter_background],
+        subsample_fraction[swift_type_dark_matter_background],
+        e->snapshot_output_count);
+  } else {
+    Ndm_background = 0;
+  }
+
+  if (with_neutrinos) {
+    Ndm_neutrino = io_count_neutrinos_to_write(
+        e->s, subsample[swift_type_neutrino],
+        subsample_fraction[swift_type_neutrino], e->snapshot_output_count);
+  } else {
+    Ndm_neutrino = 0;
+  }
+
+  /* Compute offset in the file and total number of particles */
+  size_t N[swift_type_count] = {
+      Ngas_written,   Ndm_written,         Ndm_background, Nsinks_written,
+      Nstars_written, Nblackholes_written, Ndm_neutrino};
+  long long N_total[swift_type_count] = {0};
+  long long offset[swift_type_count] = {0};
+  MPI_Exscan(N, offset, swift_type_count, MPI_LONG_LONG_INT, MPI_SUM, comm);
+  for (int ptype = 0; ptype < swift_type_count; ++ptype)
+    N_total[ptype] = offset[ptype] + N[ptype];
+
+  /* The last rank now has the correct N_total. Let's
+   * broadcast from there */
+  MPI_Bcast(N_total, swift_type_count, MPI_LONG_LONG_INT, mpi_size - 1, comm);
+
+  /* Now everybody konws its offset and the total number of
+   * particles of each type */
+
   /* Rank 0 prepares the file */
   if (mpi_rank == 0)
     prepare_file(e, fileName, xmfFileName, N_total, numFields,
@@ -1548,7 +1599,8 @@ void write_output_parallel(struct engine* e,
   /* Write the location of the particles in the arrays */
   io_write_cell_offsets(h_grp_cells, e->s->cdim, e->s->dim, e->s->cells_top,
                         e->s->nr_cells, e->s->width, mpi_rank,
-                        /*distributed=*/0, N_total, offset, numFields,
+                        /*distributed=*/0, subsample, subsample_fraction,
+                        e->snapshot_output_count, N_total, offset, numFields,
                         internal_units, snapshot_units);
 
   /* Close everything */
@@ -1672,8 +1724,10 @@ void write_output_parallel(struct engine* e,
             error("Error while allocating temporary memory for xparts");
 
           /* Collect the particles we want to write */
-          io_collect_parts_to_write(parts, xparts, parts_written,
-                                    xparts_written, Ngas, Ngas_written);
+          io_collect_parts_to_write(
+              parts, xparts, parts_written, xparts_written,
+              subsample[swift_type_gas], subsample_fraction[swift_type_gas],
+              e->snapshot_output_count, Ngas, Ngas_written);
 
           /* Select the fields to write */
           io_select_hydro_fields(parts_written, xparts_written, with_cosmology,
@@ -1714,9 +1768,11 @@ void write_output_parallel(struct engine* e,
           }
 
           /* Collect the non-inhibited DM particles from gpart */
-          io_collect_gparts_to_write(gparts, e->s->gpart_group_data,
-                                     gparts_written, gpart_group_data_written,
-                                     Ntot, Ndm_written, with_stf);
+          io_collect_gparts_to_write(
+              gparts, e->s->gpart_group_data, gparts_written,
+              gpart_group_data_written, subsample[swift_type_dark_matter],
+              subsample_fraction[swift_type_dark_matter],
+              e->snapshot_output_count, Ntot, Ndm_written, with_stf);
 
           /* Select the fields to write */
           io_select_dm_fields(gparts_written, gpart_group_data_written,
@@ -1748,7 +1804,10 @@ void write_output_parallel(struct engine* e,
         /* Collect the non-inhibited DM particles from gpart */
         io_collect_gparts_background_to_write(
             gparts, e->s->gpart_group_data, gparts_written,
-            gpart_group_data_written, Ntot, Ndm_background, with_stf);
+            gpart_group_data_written,
+            subsample[swift_type_dark_matter_background],
+            subsample_fraction[swift_type_dark_matter_background],
+            e->snapshot_output_count, Ntot, Ndm_background, with_stf);
 
         /* Select the fields to write */
         io_select_dm_fields(gparts_written, gpart_group_data_written, with_fof,
@@ -1779,7 +1838,9 @@ void write_output_parallel(struct engine* e,
         /* Collect the non-inhibited DM particles from gpart */
         io_collect_gparts_neutrino_to_write(
             gparts, e->s->gpart_group_data, gparts_written,
-            gpart_group_data_written, Ntot, Ndm_neutrino, with_stf);
+            gpart_group_data_written, subsample[swift_type_neutrino],
+            subsample_fraction[swift_type_neutrino], e->snapshot_output_count,
+            Ntot, Ndm_neutrino, with_stf);
 
         /* Select the fields to write */
         io_select_neutrino_fields(gparts_written, gpart_group_data_written,
@@ -1809,8 +1870,10 @@ void write_output_parallel(struct engine* e,
             error("Error while allocating temporary memory for sink");
 
           /* Collect the particles we want to write */
-          io_collect_sinks_to_write(sinks, sinks_written, Nsinks,
-                                    Nsinks_written);
+          io_collect_sinks_to_write(
+              sinks, sinks_written, subsample[swift_type_sink],
+              subsample_fraction[swift_type_sink], e->snapshot_output_count,
+              Nsinks, Nsinks_written);
 
           /* Select the fields to write */
           io_select_sink_fields(sinks_written, with_cosmology, with_fof,
@@ -1840,8 +1903,10 @@ void write_output_parallel(struct engine* e,
             error("Error while allocating temporary memory for sparts");
 
           /* Collect the particles we want to write */
-          io_collect_sparts_to_write(sparts, sparts_written, Nstars,
-                                     Nstars_written);
+          io_collect_sparts_to_write(
+              sparts, sparts_written, subsample[swift_type_stars],
+              subsample_fraction[swift_type_stars], e->snapshot_output_count,
+              Nstars, Nstars_written);
 
           /* Select the fields to write */
           io_select_star_fields(sparts_written, with_cosmology, with_fof,
@@ -1871,8 +1936,10 @@ void write_output_parallel(struct engine* e,
             error("Error while allocating temporary memory for bparts");
 
           /* Collect the particles we want to write */
-          io_collect_bparts_to_write(bparts, bparts_written, Nblackholes,
-                                     Nblackholes_written);
+          io_collect_bparts_to_write(
+              bparts, bparts_written, subsample[swift_type_black_hole],
+              subsample_fraction[swift_type_black_hole],
+              e->snapshot_output_count, Nblackholes, Nblackholes_written);
 
           /* Select the fields to write */
           io_select_bh_fields(bparts_written, with_cosmology, with_fof,
diff --git a/src/random.h b/src/random.h
index 9c4e35f641a508df1c94549f9635017b8b620eb0..87f1d05ede2a8322f0963907a9708467c72eb92f 100644
--- a/src/random.h
+++ b/src/random.h
@@ -42,7 +42,6 @@
  * generator.
  * In case new numbers need to be added other possible
  * numbers could be:
- * 65610001
  * 126247697
  * 193877777
  * 303595777
@@ -64,6 +63,7 @@ enum random_number_type {
   random_number_BH_feedback = 1640531371LL,
   random_number_BH_swallow = 4947009007LL,
   random_number_BH_reposition = 59969537LL,
+  random_number_snapshot_sampling = 6561001LL,
 };
 
 #ifndef __APPLE__
diff --git a/src/serial_io.c b/src/serial_io.c
index ef5567066781f06dc812b7633fa8c911359a1986..b6f5e28b294ec4a14b96397a83343b06d432bed7 100644
--- a/src/serial_io.c
+++ b/src/serial_io.c
@@ -980,34 +980,6 @@ void write_output_serial(struct engine* e,
   // const size_t Nbaryons = Ngas + Nstars;
   // const size_t Ndm = Ntot > 0 ? Ntot - Nbaryons : 0;
 
-  size_t Ndm_background = 0;
-  if (with_DM_background) {
-    Ndm_background = io_count_dm_background_gparts(gparts, Ntot);
-  }
-  size_t Ndm_neutrino = 0;
-  if (with_neutrinos) {
-    Ndm_neutrino = io_count_dm_neutrino_gparts(gparts, Ntot);
-  }
-
-  /* Number of particles that we will write
-   * Recall that background particles are never inhibited and have no extras */
-  const size_t Ntot_written =
-      e->s->nr_gparts - e->s->nr_inhibited_gparts - e->s->nr_extra_gparts;
-  const size_t Ngas_written =
-      e->s->nr_parts - e->s->nr_inhibited_parts - e->s->nr_extra_parts;
-  const size_t Nsinks_written =
-      e->s->nr_sinks - e->s->nr_inhibited_sinks - e->s->nr_extra_sinks;
-  const size_t Nstars_written =
-      e->s->nr_sparts - e->s->nr_inhibited_sparts - e->s->nr_extra_sparts;
-  const size_t Nblackholes_written =
-      e->s->nr_bparts - e->s->nr_inhibited_bparts - e->s->nr_extra_bparts;
-  const size_t Nbaryons_written =
-      Ngas_written + Nstars_written + Nblackholes_written + Nsinks_written;
-  const size_t Ndm_written =
-      Ntot_written > 0
-          ? Ntot_written - Nbaryons_written - Ndm_background - Ndm_neutrino
-          : 0;
-
   /* 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. */
@@ -1037,6 +1009,86 @@ void write_output_serial(struct engine* e,
   /* Create the directory */
   if (mpi_rank == 0) safe_checkdir(snapshot_subdir_name, /*create=*/1);
 
+  /* Do we want to sub-sample any of the arrays */
+  int subsample[swift_type_count];
+  float subsample_fraction[swift_type_count];
+  for (int i = 0; i < swift_type_count; ++i) {
+    subsample[i] = 0;
+    subsample_fraction[i] = 1.f;
+  }
+
+  output_options_get_subsampling(
+      output_options, current_selection_name, e->snapshot_subsample,
+      e->snapshot_subsample_fraction, subsample, subsample_fraction);
+
+  /* Is any particle type being subsampled? */
+  int subsample_any = 0;
+  for (int i = 0; i < swift_type_count; ++i) {
+    subsample_any += subsample[i];
+    if (!subsample[i]) subsample_fraction[i] = 1.f;
+  }
+
+  /* Number of particles that we will write */
+  size_t Ngas_written, Ndm_written, Ndm_background, Ndm_neutrino,
+      Nsinks_written, Nstars_written, Nblackholes_written;
+
+  if (subsample[swift_type_gas]) {
+    Ngas_written = io_count_gas_to_write(e->s, /*subsample=*/1,
+                                         subsample_fraction[swift_type_gas],
+                                         e->snapshot_output_count);
+  } else {
+    Ngas_written =
+        e->s->nr_parts - e->s->nr_inhibited_parts - e->s->nr_extra_parts;
+  }
+
+  if (subsample[swift_type_stars]) {
+    Nstars_written = io_count_stars_to_write(
+        e->s, /*subsample=*/1, subsample_fraction[swift_type_stars],
+        e->snapshot_output_count);
+  } else {
+    Nstars_written =
+        e->s->nr_sparts - e->s->nr_inhibited_sparts - e->s->nr_extra_sparts;
+  }
+
+  if (subsample[swift_type_black_hole]) {
+    Nblackholes_written = io_count_black_holes_to_write(
+        e->s, /*subsample=*/1, subsample_fraction[swift_type_black_hole],
+        e->snapshot_output_count);
+  } else {
+    Nblackholes_written =
+        e->s->nr_bparts - e->s->nr_inhibited_bparts - e->s->nr_extra_bparts;
+  }
+
+  if (subsample[swift_type_sink]) {
+    Nsinks_written = io_count_sinks_to_write(
+        e->s, /*subsample=*/1, subsample_fraction[swift_type_sink],
+        e->snapshot_output_count);
+  } else {
+    Nsinks_written =
+        e->s->nr_sinks - e->s->nr_inhibited_sinks - e->s->nr_extra_sinks;
+  }
+
+  Ndm_written = io_count_dark_matter_to_write(
+      e->s, subsample[swift_type_dark_matter],
+      subsample_fraction[swift_type_dark_matter], e->snapshot_output_count);
+
+  if (with_DM_background) {
+    Ndm_background = io_count_dark_matter_to_write(
+        e->s, subsample[swift_type_dark_matter_background],
+        subsample_fraction[swift_type_dark_matter_background],
+        e->snapshot_output_count);
+  } else {
+    Ndm_background = 0;
+  }
+
+  if (with_neutrinos) {
+    Ndm_neutrino = io_count_neutrinos_to_write(
+        e->s, subsample[swift_type_neutrino],
+        subsample_fraction[swift_type_neutrino], e->snapshot_output_count);
+  } else {
+    Ndm_neutrino = 0;
+  }
+
   /* Total number of fields to write per ptype */
   int numFields[swift_type_count] = {0};
   for (int ptype = 0; ptype < swift_type_count; ++ptype) {
@@ -1122,7 +1174,7 @@ void write_output_serial(struct engine* e,
     struct tm* timeinfo = localtime(&tm);
     char snapshot_date[64];
     strftime(snapshot_date, 64, "%T %F %Z", timeinfo);
-    io_write_attribute_s(h_grp, "Snapshot date", snapshot_date);
+    io_write_attribute_s(h_grp, "SnapshotDate", snapshot_date);
 
     /* GADGET-2 legacy values: Number of particles of each type */
     long long numParticlesThisFile[swift_type_count] = {0};
@@ -1156,8 +1208,14 @@ void write_output_serial(struct engine* e,
                        swift_type_count);
     io_write_attribute(h_grp, "NumFilesPerSnapshot", INT, &numFiles, 1);
     io_write_attribute_i(h_grp, "ThisFile", 0);
-    io_write_attribute_s(h_grp, "OutputType", "FullVolume");
     io_write_attribute_s(h_grp, "SelectOutput", current_selection_name);
+    if (subsample_any) {
+      io_write_attribute_s(h_grp, "OutputType", "SubSampled");
+      io_write_attribute(h_grp, "SubSampleFractions", FLOAT, subsample_fraction,
+                         swift_type_count);
+    } else {
+      io_write_attribute_s(h_grp, "OutputType", "FullVolume");
+    }
 
     /* Close header */
     H5Gclose(h_grp);
@@ -1217,7 +1275,8 @@ void write_output_serial(struct engine* e,
   /* Write the location of the particles in the arrays */
   io_write_cell_offsets(h_grp_cells, e->s->cdim, e->s->dim, e->s->cells_top,
                         e->s->nr_cells, e->s->width, mpi_rank,
-                        /*distributed=*/0, N_total, offset, numFields,
+                        /*distributed=*/0, subsample, subsample_fraction,
+                        e->snapshot_output_count, N_total, offset, numFields,
                         internal_units, snapshot_units);
 
   /* Close everything */
@@ -1299,8 +1358,10 @@ void write_output_serial(struct engine* e,
                 error("Error while allocating temporary memory for xparts");
 
               /* Collect the particles we want to write */
-              io_collect_parts_to_write(parts, xparts, parts_written,
-                                        xparts_written, Ngas, Ngas_written);
+              io_collect_parts_to_write(
+                  parts, xparts, parts_written, xparts_written,
+                  subsample[swift_type_gas], subsample_fraction[swift_type_gas],
+                  e->snapshot_output_count, Ngas, Ngas_written);
 
               /* Select the fields to write */
               io_select_hydro_fields(parts_written, xparts_written,
@@ -1346,7 +1407,9 @@ void write_output_serial(struct engine* e,
               /* Collect the non-inhibited DM particles from gpart */
               io_collect_gparts_to_write(
                   gparts, e->s->gpart_group_data, gparts_written,
-                  gpart_group_data_written, Ntot, Ndm_written, with_stf);
+                  gpart_group_data_written, subsample[swift_type_dark_matter],
+                  subsample_fraction[swift_type_dark_matter],
+                  e->snapshot_output_count, Ntot, Ndm_written, with_stf);
 
               /* Select the fields to write */
               io_select_dm_fields(gparts_written, gpart_group_data_written,
@@ -1379,7 +1442,10 @@ void write_output_serial(struct engine* e,
             /* Collect the non-inhibited DM particles from gpart */
             io_collect_gparts_background_to_write(
                 gparts, e->s->gpart_group_data, gparts_written,
-                gpart_group_data_written, Ntot, Ndm_background, with_stf);
+                gpart_group_data_written,
+                subsample[swift_type_dark_matter_background],
+                subsample_fraction[swift_type_dark_matter_background],
+                e->snapshot_output_count, Ntot, Ndm_background, with_stf);
 
             /* Select the fields to write */
             io_select_dm_fields(gparts_written, gpart_group_data_written,
@@ -1412,7 +1478,9 @@ void write_output_serial(struct engine* e,
             /* Collect the non-inhibited DM particles from gpart */
             io_collect_gparts_neutrino_to_write(
                 gparts, e->s->gpart_group_data, gparts_written,
-                gpart_group_data_written, Ntot, Ndm_neutrino, with_stf);
+                gpart_group_data_written, subsample[swift_type_neutrino],
+                subsample_fraction[swift_type_neutrino],
+                e->snapshot_output_count, Ntot, Ndm_neutrino, with_stf);
 
             /* Select the fields to write */
             io_select_neutrino_fields(gparts_written, gpart_group_data_written,
@@ -1441,8 +1509,10 @@ void write_output_serial(struct engine* e,
                 error("Error while allocating temporary memory for sinks");
 
               /* Collect the particles we want to write */
-              io_collect_sinks_to_write(sinks, sinks_written, Nsinks,
-                                        Nsinks_written);
+              io_collect_sinks_to_write(
+                  sinks, sinks_written, subsample[swift_type_sink],
+                  subsample_fraction[swift_type_sink], e->snapshot_output_count,
+                  Nsinks, Nsinks_written);
 
               /* Select the fields to write */
               io_select_sink_fields(sinks_written, with_cosmology, with_fof,
@@ -1472,8 +1542,10 @@ void write_output_serial(struct engine* e,
                 error("Error while allocating temporary memory for sparts");
 
               /* Collect the particles we want to write */
-              io_collect_sparts_to_write(sparts, sparts_written, Nstars,
-                                         Nstars_written);
+              io_collect_sparts_to_write(
+                  sparts, sparts_written, subsample[swift_type_stars],
+                  subsample_fraction[swift_type_stars],
+                  e->snapshot_output_count, Nstars, Nstars_written);
 
               /* Select the fields to write */
               io_select_star_fields(sparts_written, with_cosmology, with_fof,
@@ -1503,8 +1575,10 @@ void write_output_serial(struct engine* e,
                 error("Error while allocating temporary memory for bparts");
 
               /* Collect the particles we want to write */
-              io_collect_bparts_to_write(bparts, bparts_written, Nblackholes,
-                                         Nblackholes_written);
+              io_collect_bparts_to_write(
+                  bparts, bparts_written, subsample[swift_type_black_hole],
+                  subsample_fraction[swift_type_black_hole],
+                  e->snapshot_output_count, Nblackholes, Nblackholes_written);
 
               /* Select the fields to write */
               io_select_bh_fields(bparts_written, with_cosmology, with_fof,
diff --git a/src/single_io.c b/src/single_io.c
index adba7c6e1f29dacecb8b959190f336e6fabbe4b5..3e2666602dd72262de285ea9777ee8827c308534 100644
--- a/src/single_io.c
+++ b/src/single_io.c
@@ -824,43 +824,6 @@ void write_output_single(struct engine* e,
   const size_t Nstars = e->s->nr_sparts;
   const size_t Nsinks = e->s->nr_sinks;
   const size_t Nblackholes = e->s->nr_bparts;
-  // const size_t Nbaryons = Ngas + Nstars;
-  // const size_t Ndm = Ntot > 0 ? Ntot - Nbaryons : 0;
-
-  size_t Ndm_background = 0;
-  if (with_DM_background) {
-    Ndm_background = io_count_dm_background_gparts(gparts, Ntot);
-  }
-  size_t Ndm_neutrino = 0;
-  if (with_neutrinos) {
-    Ndm_neutrino = io_count_dm_neutrino_gparts(gparts, Ntot);
-  }
-
-  /* Number of particles that we will write
-   * Recall that background particles are never inhibited and have no extras */
-  const size_t Ntot_written =
-      e->s->nr_gparts - e->s->nr_inhibited_gparts - e->s->nr_extra_gparts;
-  const size_t Ngas_written =
-      e->s->nr_parts - e->s->nr_inhibited_parts - e->s->nr_extra_parts;
-  const size_t Nstars_written =
-      e->s->nr_sparts - e->s->nr_inhibited_sparts - e->s->nr_extra_sparts;
-  const size_t Nsinks_written =
-      e->s->nr_sinks - e->s->nr_inhibited_sinks - e->s->nr_extra_sinks;
-  const size_t Nblackholes_written =
-      e->s->nr_bparts - e->s->nr_inhibited_bparts - e->s->nr_extra_bparts;
-  const size_t Nbaryons_written =
-      Ngas_written + Nstars_written + Nblackholes_written + Nsinks_written;
-  const size_t Ndm_written =
-      Ntot_written > 0
-          ? Ntot_written - Nbaryons_written - Ndm_background - Ndm_neutrino
-          : 0;
-
-  /* Format things in a Gadget-friendly array */
-  long long N_total[swift_type_count] = {
-      (long long)Ngas_written,   (long long)Ndm_written,
-      (long long)Ndm_background, (long long)Nsinks_written,
-      (long long)Nstars_written, (long long)Nblackholes_written,
-      (long long)Ndm_neutrino};
 
   /* Determine if we are writing a reduced snapshot, and if so which
    * output selection type to use */
@@ -890,6 +853,25 @@ void write_output_single(struct engine* e,
   /* Create the directory */
   safe_checkdir(snapshot_subdir_name, /*create=*/1);
 
+  /* Do we want to sub-sample any of the arrays */
+  int subsample[swift_type_count];
+  float subsample_fraction[swift_type_count];
+  for (int i = 0; i < swift_type_count; ++i) {
+    subsample[i] = 0;
+    subsample_fraction[i] = 1.f;
+  }
+
+  output_options_get_subsampling(
+      output_options, current_selection_name, e->snapshot_subsample,
+      e->snapshot_subsample_fraction, subsample, subsample_fraction);
+
+  /* Is any particle type being subsampled? */
+  int subsample_any = 0;
+  for (int i = 0; i < swift_type_count; ++i) {
+    subsample_any += subsample[i];
+    if (!subsample[i]) subsample_fraction[i] = 1.f;
+  }
+
   /* First time, we need to create the XMF file */
   if (e->snapshot_output_count == 0) xmf_create_file(xmfFileName);
 
@@ -900,6 +882,74 @@ void write_output_single(struct engine* e,
   /* Write the part corresponding to this specific output */
   xmf_write_outputheader(xmfFile, fileName, e->time);
 
+  /* Number of particles that we will write */
+  size_t Ngas_written, Ndm_written, Ndm_background, Ndm_neutrino,
+      Nsinks_written, Nstars_written, Nblackholes_written;
+
+  if (subsample[swift_type_gas]) {
+    Ngas_written = io_count_gas_to_write(e->s, /*subsample=*/1,
+                                         subsample_fraction[swift_type_gas],
+                                         e->snapshot_output_count);
+  } else {
+    Ngas_written =
+        e->s->nr_parts - e->s->nr_inhibited_parts - e->s->nr_extra_parts;
+  }
+
+  if (subsample[swift_type_stars]) {
+    Nstars_written = io_count_stars_to_write(
+        e->s, /*subsample=*/1, subsample_fraction[swift_type_stars],
+        e->snapshot_output_count);
+  } else {
+    Nstars_written =
+        e->s->nr_sparts - e->s->nr_inhibited_sparts - e->s->nr_extra_sparts;
+  }
+
+  if (subsample[swift_type_black_hole]) {
+    Nblackholes_written = io_count_black_holes_to_write(
+        e->s, /*subsample=*/1, subsample_fraction[swift_type_black_hole],
+        e->snapshot_output_count);
+  } else {
+    Nblackholes_written =
+        e->s->nr_bparts - e->s->nr_inhibited_bparts - e->s->nr_extra_bparts;
+  }
+
+  if (subsample[swift_type_sink]) {
+    Nsinks_written = io_count_sinks_to_write(
+        e->s, /*subsample=*/1, subsample_fraction[swift_type_sink],
+        e->snapshot_output_count);
+  } else {
+    Nsinks_written =
+        e->s->nr_sinks - e->s->nr_inhibited_sinks - e->s->nr_extra_sinks;
+  }
+
+  Ndm_written = io_count_dark_matter_to_write(
+      e->s, subsample[swift_type_dark_matter],
+      subsample_fraction[swift_type_dark_matter], e->snapshot_output_count);
+
+  if (with_DM_background) {
+    Ndm_background = io_count_dark_matter_to_write(
+        e->s, subsample[swift_type_dark_matter_background],
+        subsample_fraction[swift_type_dark_matter_background],
+        e->snapshot_output_count);
+  } else {
+    Ndm_background = 0;
+  }
+
+  if (with_neutrinos) {
+    Ndm_neutrino = io_count_neutrinos_to_write(
+        e->s, subsample[swift_type_neutrino],
+        subsample_fraction[swift_type_neutrino], e->snapshot_output_count);
+  } else {
+    Ndm_neutrino = 0;
+  }
+
+  /* Format things in a Gadget-friendly array */
+  long long N_total[swift_type_count] = {
+      (long long)Ngas_written,   (long long)Ndm_written,
+      (long long)Ndm_background, (long long)Nsinks_written,
+      (long long)Nstars_written, (long long)Nblackholes_written,
+      (long long)Ndm_neutrino};
+
   /* Open file */
   /* message("Opening file '%s'.", fileName); */
   h_file = H5Fcreate(fileName, H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT);
@@ -949,7 +999,7 @@ void write_output_single(struct engine* e,
   struct tm* timeinfo = localtime(&tm);
   char snapshot_date[64];
   strftime(snapshot_date, 64, "%T %F %Z", timeinfo);
-  io_write_attribute_s(h_grp, "Snapshot date", snapshot_date);
+  io_write_attribute_s(h_grp, "SnapshotDate", snapshot_date);
 
   /* GADGET-2 legacy values: number of particles of each type */
   long long numParticlesThisFile[swift_type_count] = {0};
@@ -989,9 +1039,16 @@ void write_output_single(struct engine* e,
                      swift_type_count);
   io_write_attribute(h_grp, "NumFilesPerSnapshot", INT, &numFiles, 1);
   io_write_attribute_i(h_grp, "ThisFile", 0);
-  io_write_attribute_s(h_grp, "OutputType", "FullVolume");
   io_write_attribute_s(h_grp, "SelectOutput", current_selection_name);
 
+  if (subsample_any) {
+    io_write_attribute_s(h_grp, "OutputType", "SubSampled");
+    io_write_attribute(h_grp, "SubSampleFractions", FLOAT, subsample_fraction,
+                       swift_type_count);
+  } else {
+    io_write_attribute_s(h_grp, "OutputType", "FullVolume");
+  }
+
   /* Close header */
   H5Gclose(h_grp);
 
@@ -1006,8 +1063,9 @@ void write_output_single(struct engine* e,
   /* Write the location of the particles in the arrays */
   io_write_cell_offsets(h_grp, e->s->cdim, e->s->dim, e->s->cells_top,
                         e->s->nr_cells, e->s->width, e->nodeID,
-                        /*distributed=*/0, N_total, global_offsets, numFields,
-                        internal_units, snapshot_units);
+                        /*distributed=*/0, subsample, subsample_fraction,
+                        e->snapshot_output_count, N_total, global_offsets,
+                        numFields, internal_units, snapshot_units);
   H5Gclose(h_grp);
 
   /* Loop over all particle types */
@@ -1082,8 +1140,10 @@ void write_output_single(struct engine* e,
             error("Error while allocating temporary memory for xparts");
 
           /* Collect the particles we want to write */
-          io_collect_parts_to_write(parts, xparts, parts_written,
-                                    xparts_written, Ngas, Ngas_written);
+          io_collect_parts_to_write(
+              parts, xparts, parts_written, xparts_written,
+              subsample[swift_type_gas], subsample_fraction[swift_type_gas],
+              e->snapshot_output_count, Ngas, Ngas_written);
 
           /* Select the fields to write */
           io_select_hydro_fields(parts_written, xparts_written, with_cosmology,
@@ -1125,9 +1185,11 @@ void write_output_single(struct engine* e,
           }
 
           /* Collect the non-inhibited DM particles from gpart */
-          io_collect_gparts_to_write(gparts, e->s->gpart_group_data,
-                                     gparts_written, gpart_group_data_written,
-                                     Ntot, Ndm_written, with_stf);
+          io_collect_gparts_to_write(
+              gparts, e->s->gpart_group_data, gparts_written,
+              gpart_group_data_written, subsample[swift_type_dark_matter],
+              subsample_fraction[swift_type_dark_matter],
+              e->snapshot_output_count, Ntot, Ndm_written, with_stf);
 
           /* Select the fields to write */
           io_select_dm_fields(gparts_written, gpart_group_data_written,
@@ -1159,7 +1221,10 @@ void write_output_single(struct engine* e,
         /* Collect the non-inhibited DM particles from gpart */
         io_collect_gparts_background_to_write(
             gparts, e->s->gpart_group_data, gparts_written,
-            gpart_group_data_written, Ntot, Ndm_background, with_stf);
+            gpart_group_data_written,
+            subsample[swift_type_dark_matter_background],
+            subsample_fraction[swift_type_dark_matter_background],
+            e->snapshot_output_count, Ntot, Ndm_background, with_stf);
 
         /* Select the fields to write */
         io_select_dm_fields(gparts_written, gpart_group_data_written, with_fof,
@@ -1191,7 +1256,9 @@ void write_output_single(struct engine* e,
         /* Collect the non-inhibited DM particles from gpart */
         io_collect_gparts_neutrino_to_write(
             gparts, e->s->gpart_group_data, gparts_written,
-            gpart_group_data_written, Ntot, Ndm_neutrino, with_stf);
+            gpart_group_data_written, subsample[swift_type_neutrino],
+            subsample_fraction[swift_type_neutrino], e->snapshot_output_count,
+            Ntot, Ndm_neutrino, with_stf);
 
         /* Select the fields to write */
         io_select_neutrino_fields(gparts_written, gpart_group_data_written,
@@ -1220,8 +1287,10 @@ void write_output_single(struct engine* e,
             error("Error while allocating temporary memory for sinks");
 
           /* Collect the particles we want to write */
-          io_collect_sinks_to_write(sinks, sinks_written, Nsinks,
-                                    Nsinks_written);
+          io_collect_sinks_to_write(
+              sinks, sinks_written, subsample[swift_type_sink],
+              subsample_fraction[swift_type_sink], e->snapshot_output_count,
+              Nsinks, Nsinks_written);
 
           /* Select the fields to write */
           io_select_sink_fields(sinks_written, with_cosmology, with_fof,
@@ -1251,8 +1320,10 @@ void write_output_single(struct engine* e,
             error("Error while allocating temporary memory for sparts");
 
           /* Collect the particles we want to write */
-          io_collect_sparts_to_write(sparts, sparts_written, Nstars,
-                                     Nstars_written);
+          io_collect_sparts_to_write(
+              sparts, sparts_written, subsample[swift_type_stars],
+              subsample_fraction[swift_type_stars], e->snapshot_output_count,
+              Nstars, Nstars_written);
 
           /* Select the fields to write */
           io_select_star_fields(sparts_written, with_cosmology, with_fof,
@@ -1282,8 +1353,10 @@ void write_output_single(struct engine* e,
             error("Error while allocating temporary memory for bparts");
 
           /* Collect the particles we want to write */
-          io_collect_bparts_to_write(bparts, bparts_written, Nblackholes,
-                                     Nblackholes_written);
+          io_collect_bparts_to_write(
+              bparts, bparts_written, subsample[swift_type_black_hole],
+              subsample_fraction[swift_type_black_hole],
+              e->snapshot_output_count, Nblackholes, Nblackholes_written);
 
           /* Select the fields to write */
           io_select_bh_fields(bparts_written, with_cosmology, with_fof,