diff --git a/doc/RTD/source/ParameterFiles/parameter_description.rst b/doc/RTD/source/ParameterFiles/parameter_description.rst
index 387556e2b8adc0b65c2ec4164b5c163b7fcaa373..eab1e670cd8e3dcdef153b2ce5b3e3c58565df14 100644
--- a/doc/RTD/source/ParameterFiles/parameter_description.rst
+++ b/doc/RTD/source/ParameterFiles/parameter_description.rst
@@ -631,6 +631,15 @@ the case of non-cosmological runs, the time of the first snapshot is expressed
 in the internal units of time. Users also have to provide the difference in time
 (or scale-factor) between consecutive outputs:
 
+* Directory in which to write snapshots: ``subdir``.
+  (default: empty string).
+
+If this is set then the full path to the snapshot files will be generated by
+taking this value and appending a slash and then the snapshot file name
+described above - e.g. ``subdir/base_name_1234.hdf5``. The directory is
+created if necessary. Any VELOCIraptor output produced by the run is also written
+to this directory.
+
 * Time difference between consecutive outputs: ``delta_time``.
 
 In non-cosmological runs this is also expressed in internal units. For
@@ -1193,13 +1202,28 @@ scale-factor of the next call. This implies that the outputs are
 equally spaced in :math:`\log(a)` (See :ref:`Output_list_label` to have
 calls not regularly spaced in time).
 
+Since VELOCIraptor produces many small output files when running with MPI,
+it can be useful to make a separate directory for each output time:
+
+* Base name of directory created for each VELOCIraptor output: ``subdir_per_output``
+  (default: empty string).
+
+If this is set then a new directory is created each time VELOCIraptor is run.
+The directory name will be subdir_per_output followed by the same output number
+used in the filenames. Note that this directory is relative to the ``subdir`` parameter
+from the Snapshots section if that is set.
+
+By default this is an empty string, which means that all VELOCIraptor outputs will
+be written to a single directory.
+
 Showing all the parameters for a basic cosmologica test-case, one would have:
 
 .. code:: YAML
 
    StructureFinding:
     config_file_name:     my_stf_configuration_file.cfg  # See the VELOCIraptor manual for the content of this file.
-    basename:             ./haloes/                      # Write the catalogs in this sub-directory
+    basename:             haloes                         # Base name for VELOCIraptor output files
+    subdir_per_output:    stf                            # Make a stf_XXXX subdirectory for each output
     scale_factor_first:   0.1                            # Scale-factor of the first output
     delta_time:           1.1                            # Delta log-a between outputs
 
diff --git a/examples/parameter_example.yml b/examples/parameter_example.yml
index 7e753c4a1cf06f4dc28f9e80534df63be1565cb4..fcba1f33795edcd1de2b12b123016f10f218ad7f 100644
--- a/examples/parameter_example.yml
+++ b/examples/parameter_example.yml
@@ -120,7 +120,8 @@ TimeIntegration:
 
 # Parameters governing the snapshots
 Snapshots:
-  basename:   output      # Common part of the name of output files
+  basename:   output      # Common part of the name of output files.
+  subdir:     dir         # (Optional) Sub-directory in which to write the snapshots. Defaults to "" (i.e. the directory where SWIFT is run).
   scale_factor_first: 0.1 # (Optional) Scale-factor of the first snapshot if cosmological time-integration.
   time_first: 0.          # (Optional) Time of the first output if non-cosmological time-integration (in internal units)
   delta_time: 0.01        # Time difference between consecutive outputs (in internal units)
@@ -202,7 +203,8 @@ DomainDecomposition:
 # Structure finding options (requires velociraptor)
 StructureFinding:
   config_file_name:     stf_input.cfg # Name of the STF config file.
-  basename:             ./stf         # Common part of the name of output files.
+  basename:             haloes        # Common part of the name of output files.
+  subdir_per_output:    stf           # (Optional) Sub-directory (within Snapshots:subdir) to use for each invocation of STF. Defaults to "" (i.e. no sub-directory) 
   scale_factor_first:   0.92          # (Optional) Scale-factor of the first snaphot (cosmological run)
   time_first:           0.01          # (Optional) Time of the first structure finding output (in internal units).
   delta_time:           1.10          # (Optional) Time difference between consecutive structure finding outputs (in internal units) in simulation time intervals.
diff --git a/src/engine.c b/src/engine.c
index e93855c9d26c2a1e7be9ba30d7a77017576d90ab..737a84daecffe097a20a3f5df7606375b82a4a6b 100644
--- a/src/engine.c
+++ b/src/engine.c
@@ -33,6 +33,8 @@
 #include <stdio.h>
 #include <stdlib.h>
 #include <string.h>
+#include <sys/stat.h>
+#include <sys/types.h>
 #include <unistd.h>
 
 /* MPI headers. */
@@ -3211,21 +3213,41 @@ void engine_dump_snapshot(struct engine *e) {
   engine_collect_stars_counter(e);
 #endif
 
+  /* Determine snapshot location */
+  char snapshotBase[FILENAME_BUFFER_SIZE];
+  if (strnlen(e->snapshot_subdir, PARSER_MAX_LINE_SIZE) > 0) {
+    if (snprintf(snapshotBase, FILENAME_BUFFER_SIZE, "%s/%s",
+                 e->snapshot_subdir,
+                 e->snapshot_base_name) >= FILENAME_BUFFER_SIZE) {
+      error(
+          "FILENAME_BUFFER_SIZE is too small for snapshot path and file name");
+    }
+      /* Try to ensure the directory exists */
+#ifdef WITH_MPI
+    if (engine_rank == 0) mkdir(e->snapshot_subdir, 0777);
+    MPI_Barrier(MPI_COMM_WORLD);
+#else
+    mkdir(e->snapshot_subdir, 0777);
+#endif
+  } else {
+    if (snprintf(snapshotBase, FILENAME_BUFFER_SIZE, "%s",
+                 e->snapshot_base_name) >= FILENAME_BUFFER_SIZE) {
+      error("FILENAME_BUFFER_SIZE is too small for snapshot file name");
+    }
+  }
+
 /* Dump... */
 #if defined(HAVE_HDF5)
 #if defined(WITH_MPI)
 #if defined(HAVE_PARALLEL_HDF5)
-  write_output_parallel(e, e->snapshot_base_name, e->internal_units,
-                        e->snapshot_units, e->nodeID, e->nr_nodes,
-                        MPI_COMM_WORLD, MPI_INFO_NULL);
+  write_output_parallel(e, snapshotBase, e->internal_units, e->snapshot_units,
+                        e->nodeID, e->nr_nodes, MPI_COMM_WORLD, MPI_INFO_NULL);
 #else
-  write_output_serial(e, e->snapshot_base_name, e->internal_units,
-                      e->snapshot_units, e->nodeID, e->nr_nodes, MPI_COMM_WORLD,
-                      MPI_INFO_NULL);
+  write_output_serial(e, snapshotBase, e->internal_units, e->snapshot_units,
+                      e->nodeID, e->nr_nodes, MPI_COMM_WORLD, MPI_INFO_NULL);
 #endif
 #else
-  write_output_single(e, e->snapshot_base_name, e->internal_units,
-                      e->snapshot_units);
+  write_output_single(e, snapshotBase, e->internal_units, e->snapshot_units);
 #endif
 #endif
 
@@ -3412,6 +3434,8 @@ void engine_init(struct engine *e, struct space *s, struct swift_params *params,
       parser_get_opt_param_double(params, "Snapshots:delta_time", -1.);
   e->ti_next_snapshot = 0;
   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);
   e->snapshot_compression =
       parser_get_opt_param_int(params, "Snapshots:compression", 0);
   e->snapshot_int_time_label_on =
@@ -3505,6 +3529,9 @@ void engine_init(struct engine *e, struct space *s, struct swift_params *params,
   if (e->policy & engine_policy_structure_finding) {
     parser_get_param_string(params, "StructureFinding:basename",
                             e->stf_base_name);
+    parser_get_opt_param_string(params, "StructureFinding:subdir_per_output",
+                                e->stf_subdir_per_output,
+                                engine_default_stf_subdir_per_output);
     parser_get_param_string(params, "StructureFinding:config_file_name",
                             e->stf_config_file_name);
 
diff --git a/src/engine.h b/src/engine.h
index 8213212ee7fda802d595d1d009fb43fd803e948e..ca10d34467f07fa86fbb1544cc61606c655c6dfe 100644
--- a/src/engine.h
+++ b/src/engine.h
@@ -113,6 +113,8 @@ enum engine_step_properties {
 #define engine_max_sparts_per_ghost_default 1000
 #define engine_star_resort_task_depth_default 2
 #define engine_tasks_per_cell_margin 1.2
+#define engine_default_stf_subdir_per_output ""
+#define engine_default_snapshot_subdir ""
 
 /**
  * @brief The rank of the engine as a global variable (for messages).
@@ -286,6 +288,7 @@ struct engine {
   integertime_t ti_next_snapshot;
 
   char snapshot_base_name[PARSER_MAX_LINE_SIZE];
+  char snapshot_subdir[PARSER_MAX_LINE_SIZE];
   int snapshot_compression;
   int snapshot_int_time_label_on;
   int snapshot_invoke_stf;
@@ -305,6 +308,7 @@ struct engine {
 
   char stf_config_file_name[PARSER_MAX_LINE_SIZE];
   char stf_base_name[PARSER_MAX_LINE_SIZE];
+  char stf_subdir_per_output[PARSER_MAX_LINE_SIZE];
   int stf_output_count;
 
   /* FoF black holes seeding information */
diff --git a/src/velociraptor_interface.c b/src/velociraptor_interface.c
index 6e2f4e298f9ab963fcf9871d27388a1ef0a2eeb6..6f6f76bc12ea9a78fd661fe8664c054e500c2576 100644
--- a/src/velociraptor_interface.c
+++ b/src/velociraptor_interface.c
@@ -21,6 +21,10 @@
 #include "../config.h"
 
 /* Some standard headers. */
+#include <stdio.h>
+#include <string.h>
+#include <sys/stat.h>
+#include <sys/types.h>
 #include <unistd.h>
 
 /* This object's header. */
@@ -566,15 +570,52 @@ void velociraptor_invoke(struct engine *e, const int linked_with_snap) {
         "VELOCIraptor conf: MPI rank %d sending %zu gparts to VELOCIraptor.",
         engine_rank, nr_gparts);
 
-  /* Append base name with the current output number */
-  char outputFileName[PARSER_MAX_LINE_SIZE + 128];
+  /* Generate directory name for this output - start with snapshot directory, if
+   * specified */
+  char outputDirName[FILENAME_BUFFER_SIZE] = "";
+  if (strnlen(e->snapshot_subdir, PARSER_MAX_LINE_SIZE) > 0) {
+    if (snprintf(outputDirName, FILENAME_BUFFER_SIZE, "%s/",
+                 e->snapshot_subdir) >= FILENAME_BUFFER_SIZE) {
+      error("FILENAME_BUFFER_SIZE is to small for snapshot directory name!");
+    }
+#ifdef WITH_MPI
+    if (engine_rank == 0) mkdir(outputDirName, 0777);
+    MPI_Barrier(MPI_COMM_WORLD);
+#else
+    mkdir(outputDirName, 0777);
+#endif
+  }
+
+  /* Then create output-specific subdirectory if necessary */
+  char subDirName[FILENAME_BUFFER_SIZE] = "";
+  if (strnlen(e->stf_subdir_per_output, PARSER_MAX_LINE_SIZE) > 0) {
+    if (snprintf(subDirName, FILENAME_BUFFER_SIZE, "%s%s_%04i/", outputDirName,
+                 e->stf_subdir_per_output,
+                 e->snapshot_output_count) >= FILENAME_BUFFER_SIZE) {
+      error(
+          "FILENAME_BUFFER_SIZE is to small for Velociraptor directory name!");
+    }
+#ifdef WITH_MPI
+    if (engine_rank == 0) mkdir(subDirName, 0777);
+    MPI_Barrier(MPI_COMM_WORLD);
+#else
+    mkdir(subDirName, 0777);
+#endif
+  } else {
+    /* Not making separate directories so subDirName=outputDirName */
+    strncpy(subDirName, outputDirName, FILENAME_BUFFER_SIZE);
+  }
 
   /* What is the snapshot number? */
   int snapnum = e->stf_output_count;
 
   /* What should the filename be? */
-  snprintf(outputFileName, PARSER_MAX_LINE_SIZE + 128, "%s_%04i.VELOCIraptor",
-           e->stf_base_name, snapnum);
+  char outputFileName[FILENAME_BUFFER_SIZE];
+  if (snprintf(outputFileName, FILENAME_BUFFER_SIZE, "%s%s_%04i.VELOCIraptor",
+               subDirName, e->stf_base_name,
+               e->stf_output_count) >= FILENAME_BUFFER_SIZE) {
+    error("FILENAME_BUFFER_SIZE is too small for Velociraptor file name!");
+  }
 
   tic = getticks();
 
diff --git a/src/xmf.c b/src/xmf.c
index 3c8bb257cca224a06e5b516be46e92dd0006183a..d3d88627b8fe0c0e40df53c6eb79960ef82f1519 100644
--- a/src/xmf.c
+++ b/src/xmf.c
@@ -107,6 +107,7 @@ void xmf_create_file(const char* baseName) {
   char fileName[FILENAME_BUFFER_SIZE];
   snprintf(fileName, FILENAME_BUFFER_SIZE, "%s.xmf", baseName);
   FILE* xmfFile = fopen(fileName, "w");
+  if (xmfFile == NULL) error("Unable to create XMF file.");
 
   fprintf(xmfFile, "<?xml version=\"1.0\" ?> \n");
   fprintf(xmfFile, "<!DOCTYPE Xdmf SYSTEM \"Xdmf.dtd\" []> \n");