diff --git a/src/black_holes/EAGLE/black_holes_io.h b/src/black_holes/EAGLE/black_holes_io.h
index cc6cb962db61f9e45ded990b9dca7f00db075733..7f67f340e5e76ff62b09f9f804709cdd9dd5dfc6 100644
--- a/src/black_holes/EAGLE/black_holes_io.h
+++ b/src/black_holes/EAGLE/black_holes_io.h
@@ -29,9 +29,9 @@
  * @param list The list of i/o properties to read.
  * @param num_fields The number of i/o fields to read.
  */
-INLINE static void black_holes_read_particles(struct bpart *bparts,
-                                              struct io_props *list,
-                                              int *num_fields) {
+INLINE static void black_holes_read_particles(struct bpart* bparts,
+                                              struct io_props* list,
+                                              int* num_fields) {
 
   /* Say how much we want to read */
   *num_fields = 6;
@@ -51,6 +51,53 @@ INLINE static void black_holes_read_particles(struct bpart *bparts,
                                 UNIT_CONV_ENERGY, bparts, energy_reservoir);
 }
 
+INLINE static void convert_bpart_pos(const struct engine* e,
+                                     const struct bpart* bp, double* ret) {
+
+  if (e->s->periodic) {
+    ret[0] = box_wrap(bp->x[0], 0.0, e->s->dim[0]);
+    ret[1] = box_wrap(bp->x[1], 0.0, e->s->dim[1]);
+    ret[2] = box_wrap(bp->x[2], 0.0, e->s->dim[2]);
+  } else {
+    ret[0] = bp->x[0];
+    ret[1] = bp->x[1];
+    ret[2] = bp->x[2];
+  }
+}
+
+INLINE static void convert_bpart_vel(const struct engine* e,
+                                     const struct bpart* bp, float* ret) {
+
+  const int with_cosmology = (e->policy & engine_policy_cosmology);
+  const struct cosmology* cosmo = e->cosmology;
+  const integertime_t ti_current = e->ti_current;
+  const double time_base = e->time_base;
+
+  const integertime_t ti_beg = get_integer_time_begin(ti_current, bp->time_bin);
+  const integertime_t ti_end = get_integer_time_end(ti_current, bp->time_bin);
+
+  /* Get time-step since the last kick */
+  float dt_kick_grav;
+  if (with_cosmology) {
+    dt_kick_grav = cosmology_get_grav_kick_factor(cosmo, ti_beg, ti_current);
+    dt_kick_grav -=
+        cosmology_get_grav_kick_factor(cosmo, ti_beg, (ti_beg + ti_end) / 2);
+  } else {
+    dt_kick_grav = (ti_current - ((ti_beg + ti_end) / 2)) * time_base;
+  }
+
+  /* Extrapolate the velocites to the current time */
+  const struct gpart* gp = bp->gpart;
+  ret[0] = gp->v_full[0] + gp->a_grav[0] * dt_kick_grav;
+  ret[1] = gp->v_full[1] + gp->a_grav[1] * dt_kick_grav;
+  ret[2] = gp->v_full[2] + gp->a_grav[2] * dt_kick_grav;
+
+  /* Conversion from internal units to peculiar velocities */
+  ret[0] *= cosmo->a_inv;
+  ret[1] *= cosmo->a_inv;
+  ret[2] *= cosmo->a_inv;
+}
+
 /**
  * @brief Specifies which b-particle fields to write to a dataset
  *
@@ -58,40 +105,72 @@ INLINE static void black_holes_read_particles(struct bpart *bparts,
  * @param list The list of i/o properties to write.
  * @param num_fields The number of i/o fields to write.
  */
-INLINE static void black_holes_write_particles(const struct bpart *bparts,
-                                               struct io_props *list,
-                                               int *num_fields) {
+INLINE static void black_holes_write_particles(const struct bpart* bparts,
+                                               struct io_props* list,
+                                               int* num_fields,
+                                               int with_cosmology) {
 
   /* Say how much we want to write */
   *num_fields = 12;
 
   /* List what we want to write */
-  list[0] = io_make_output_field("Coordinates", DOUBLE, 3, UNIT_CONV_LENGTH,
-                                 bparts, x);
-  list[1] =
-      io_make_output_field("Velocities", FLOAT, 3, UNIT_CONV_SPEED, bparts, v);
+  list[0] = io_make_output_field_convert_bpart(
+      "Coordinates", DOUBLE, 3, UNIT_CONV_LENGTH, 1.f, bparts,
+      convert_bpart_pos, "Co-moving position of the particles");
+
+  list[1] = io_make_output_field_convert_bpart(
+      "Velocities", FLOAT, 3, UNIT_CONV_SPEED, 0.f, bparts, convert_bpart_vel,
+      "Peculiar velocities of the particles. This is a * dx/dt where x is the "
+      "co-moving position of the particles.");
+
   list[2] =
-      io_make_output_field("Masses", FLOAT, 1, UNIT_CONV_MASS, bparts, mass);
+      io_make_output_field("DynamicalMasses", FLOAT, 1, UNIT_CONV_MASS, 0.f,
+                           bparts, mass, "Dynamical masses of the particles");
+
   list[3] = io_make_output_field("ParticleIDs", LONGLONG, 1, UNIT_CONV_NO_UNITS,
-                                 bparts, id);
-  list[4] = io_make_output_field("SmoothingLength", FLOAT, 1, UNIT_CONV_LENGTH,
-                                 bparts, h);
-  list[5] = io_make_output_field("SubgridMasses", FLOAT, 1, UNIT_CONV_MASS,
-                                 bparts, subgrid_mass);
-  list[6] = io_make_output_field("FormationTime", FLOAT, 1, UNIT_CONV_TIME,
-                                 bparts, formation_time);
-  list[7] = io_make_output_field("GasDensity", FLOAT, 1, UNIT_CONV_DENSITY,
-                                 bparts, rho_gas);
-  list[8] = io_make_output_field("GasSoundSpeed", FLOAT, 1, UNIT_CONV_SPEED,
-                                 bparts, sound_speed_gas);
-  list[9] = io_make_output_field("EnergyReservoir", FLOAT, 1, UNIT_CONV_ENERGY,
-                                 bparts, energy_reservoir);
-  list[10] = io_make_output_field("AccretionRate", FLOAT, 1,
-                                  UNIT_CONV_MASS_PER_UNIT_TIME, bparts,
-                                  accretion_rate);
-  list[11] = io_make_output_field("TotalAccretedMass", FLOAT, 1,
-                                  UNIT_CONV_MASS_PER_UNIT_TIME, bparts,
-                                  total_accreted_mass);
+                                 0.f, bparts, id, "Unique ID of the particles");
+
+  list[4] = io_make_output_field(
+      "SmoothingLengths", FLOAT, 1, UNIT_CONV_LENGTH, 1.f, bparts, h,
+      "Co-moving smoothing lengths (FWHM of the kernel) of the particles");
+
+  list[5] = io_make_output_field("SubgridMasses", FLOAT, 1, UNIT_CONV_MASS, 0.f,
+                                 bparts, subgrid_mass,
+                                 "Subgrid masses of the particles");
+
+  if (with_cosmology) {
+    list[6] = io_make_output_field(
+        "FormationScaleFactors", FLOAT, 1, UNIT_CONV_NO_UNITS, 0.f, bparts,
+        formation_scale_factor, "Scale-factors at which the BHs were formed");
+  } else {
+    list[6] = io_make_output_field("FormationTimes", FLOAT, 1, UNIT_CONV_TIME,
+                                   0.f, bparts, formation_time,
+                                   "Times at which the BHs were formed");
+  }
+
+  list[7] = io_make_output_field(
+      "GasDensities", FLOAT, 1, UNIT_CONV_DENSITY, 0.f, bparts, rho_gas,
+      "Co-moving densities of the gas around the particles");
+
+  list[8] = io_make_output_field(
+      "GasSoundSpeeds", FLOAT, 1, UNIT_CONV_SPEED, 1.5f * hydro_gamma_minus_one,
+      bparts, sound_speed_gas,
+      "Co-moving sound-speeds of the gas around the particles");
+
+  list[9] = io_make_output_field(
+      "EnergyReservoirs", FLOAT, 1, UNIT_CONV_ENERGY, 0.f, bparts,
+      energy_reservoir,
+      "Physcial energy contained in the feedback reservoir of the particles");
+
+  list[10] = io_make_output_field(
+      "AccretionRates", FLOAT, 1, UNIT_CONV_MASS_PER_UNIT_TIME, 0.f, bparts,
+      accretion_rate,
+      "Physical instantaneous accretion rates of the particles");
+
+  list[11] = io_make_output_field(
+      "TotalAccretedMasses", FLOAT, 1, UNIT_CONV_MASS_PER_UNIT_TIME, 0.f,
+      bparts, total_accreted_mass,
+      "Total masses accreted onto the particles since its birth");
 
 #ifdef DEBUG_INTERACTIONS_BLACK_HOLES
 
diff --git a/src/black_holes/EAGLE/black_holes_part.h b/src/black_holes/EAGLE/black_holes_part.h
index b6d15dfa06be1c4572845d1b9656b053439f70d1..e79acb722aec1857d4d8f5fc26b0b0cf445c086c 100644
--- a/src/black_holes/EAGLE/black_holes_part.h
+++ b/src/black_holes/EAGLE/black_holes_part.h
@@ -19,6 +19,8 @@
 #ifndef SWIFT_EAGLE_BLACK_HOLE_PART_H
 #define SWIFT_EAGLE_BLACK_HOLE_PART_H
 
+#include "timeline.h"
+
 /**
  * @brief Particle fields for the black hole particles.
  *
diff --git a/src/chemistry/EAGLE/chemistry_io.h b/src/chemistry/EAGLE/chemistry_io.h
index 62de85d47329a588d658cc61db72510337988638..810eb2dd857addcb3bc6707c9f90be9ee4fc5c4f 100644
--- a/src/chemistry/EAGLE/chemistry_io.h
+++ b/src/chemistry/EAGLE/chemistry_io.h
@@ -58,50 +58,72 @@ INLINE static int chemistry_write_particles(const struct part* parts,
                                             struct io_props* list) {
 
   /* List what we want to write */
-  list[0] = io_make_output_field("ElementAbundance", FLOAT,
-                                 chemistry_element_count, UNIT_CONV_NO_UNITS,
-                                 parts, chemistry_data.metal_mass_fraction);
+  list[0] = io_make_output_field(
+      "ElementMassFractions", FLOAT, chemistry_element_count,
+      UNIT_CONV_NO_UNITS, 0.f, parts, chemistry_data.metal_mass_fraction,
+      "Fractions of the particles' masses that are in the given element");
 
   list[1] = io_make_output_field(
       "SmoothedElementAbundance", FLOAT, chemistry_element_count,
-      UNIT_CONV_NO_UNITS, parts, chemistry_data.smoothed_metal_mass_fraction);
+      UNIT_CONV_NO_UNITS, 0.f, parts,
+      chemistry_data.smoothed_metal_mass_fraction,
+      "Smoothed fractions of the particles' masses that are "
+      "in the given element");
 
-  list[2] =
-      io_make_output_field("Metallicity", FLOAT, 1, UNIT_CONV_NO_UNITS, parts,
-                           chemistry_data.metal_mass_fraction_total);
+  list[2] = io_make_output_field(
+      "MetalMassFractions", FLOAT, 1, UNIT_CONV_NO_UNITS, 0.f, parts,
+      chemistry_data.metal_mass_fraction_total,
+      "Fractions of the particles' masses that are in metals");
 
   list[3] = io_make_output_field(
-      "SmoothedMetallicity", FLOAT, 1, UNIT_CONV_NO_UNITS, parts,
-      chemistry_data.smoothed_metal_mass_fraction_total);
-
-  list[4] = io_make_output_field("TotalMassFromSNIa", FLOAT, 1, UNIT_CONV_MASS,
-                                 parts, chemistry_data.mass_from_SNIa);
-
-  list[5] = io_make_output_field("MetalMassFracFromSNIa", FLOAT, 1,
-                                 UNIT_CONV_NO_UNITS, parts,
-                                 chemistry_data.metal_mass_fraction_from_SNIa);
-
-  list[6] = io_make_output_field("TotalMassFromAGB", FLOAT, 1, UNIT_CONV_MASS,
-                                 parts, chemistry_data.mass_from_AGB);
-
-  list[7] =
-      io_make_output_field("MetalMassFracFromAGB", FLOAT, 1, UNIT_CONV_NO_UNITS,
-                           parts, chemistry_data.metal_mass_fraction_from_AGB);
-
-  list[8] = io_make_output_field("TotalMassFromSNII", FLOAT, 1, UNIT_CONV_MASS,
-                                 parts, chemistry_data.mass_from_SNII);
-
-  list[9] = io_make_output_field("MetalMassFracFromSNII", FLOAT, 1,
-                                 UNIT_CONV_NO_UNITS, parts,
-                                 chemistry_data.metal_mass_fraction_from_SNII);
-
-  list[10] =
-      io_make_output_field("IronMassFracFromSNIa", FLOAT, 1, UNIT_CONV_NO_UNITS,
-                           parts, chemistry_data.iron_mass_fraction_from_SNIa);
+      "SmoothedMetalMassFractions", FLOAT, 1, UNIT_CONV_NO_UNITS, 0.f, parts,
+      chemistry_data.smoothed_metal_mass_fraction_total,
+      "Smoothed fractions of the particles masses that are in metals");
+
+  list[4] = io_make_output_field(
+      "MassesFromSNIa", FLOAT, 1, UNIT_CONV_MASS, 0.f, parts,
+      chemistry_data.mass_from_SNIa,
+      "Masses of gas that have been produced by SNIa stars");
+
+  list[5] = io_make_output_field("MetalMassFractionsFromSNIa", FLOAT, 1,
+                                 UNIT_CONV_NO_UNITS, 0.f, parts,
+                                 chemistry_data.metal_mass_fraction_from_SNIa,
+                                 "Fractions of the particles' masses that are "
+                                 "in metals produced by SNIa stars");
+
+  list[6] = io_make_output_field(
+      "MassesFromAGB", FLOAT, 1, UNIT_CONV_MASS, 0.f, parts,
+      chemistry_data.mass_from_AGB,
+      "Masses of gas that have been produced by AGN stars");
+
+  list[7] = io_make_output_field("MetalMassFractionsFromAGB", FLOAT, 1,
+                                 UNIT_CONV_NO_UNITS, 0., parts,
+                                 chemistry_data.metal_mass_fraction_from_AGB,
+                                 "Fractions of the particles' masses that are "
+                                 "in metals produced by AGB stars");
+
+  list[8] = io_make_output_field(
+      "MassesFromSNII", FLOAT, 1, UNIT_CONV_MASS, 0.f, parts,
+      chemistry_data.mass_from_SNII,
+      "Masses of gas that have been produced by SNII stars");
+
+  list[9] = io_make_output_field("MetalMassFractionsFromSNII", FLOAT, 1,
+                                 UNIT_CONV_NO_UNITS, 0.f, parts,
+                                 chemistry_data.metal_mass_fraction_from_SNII,
+                                 "Fractions of the particles' masses that are "
+                                 "in metals produced by SNII stars");
+
+  list[10] = io_make_output_field("IronMassFractionsFromSNIa", FLOAT, 1,
+                                  UNIT_CONV_NO_UNITS, 0.f, parts,
+                                  chemistry_data.iron_mass_fraction_from_SNIa,
+                                  "Fractions of the particles' masses that are "
+                                  "in iron produced by SNIa stars");
 
   list[11] = io_make_output_field(
-      "SmoothedIronMassFracFromSNIa", FLOAT, 1, UNIT_CONV_NO_UNITS, parts,
-      chemistry_data.smoothed_iron_mass_fraction_from_SNIa);
+      "SmoothedIronMassFractionsFromSNIa", FLOAT, 1, UNIT_CONV_NO_UNITS, 0.f,
+      parts, chemistry_data.smoothed_iron_mass_fraction_from_SNIa,
+      "Smoothed fractions of the particles' masses that are "
+      "in iron produced by SNIa stars");
 
   return 12;
 }
@@ -118,50 +140,72 @@ INLINE static int chemistry_write_sparticles(const struct spart* sparts,
                                              struct io_props* list) {
 
   /* List what we want to write */
-  list[0] = io_make_output_field("ElementAbundance", FLOAT,
-                                 chemistry_element_count, UNIT_CONV_NO_UNITS,
-                                 sparts, chemistry_data.metal_mass_fraction);
+  list[0] = io_make_output_field(
+      "ElementMassFractions", FLOAT, chemistry_element_count,
+      UNIT_CONV_NO_UNITS, 0.f, sparts, chemistry_data.metal_mass_fraction,
+      "Fractions of the particles' masses that are in the given element");
 
   list[1] = io_make_output_field(
       "SmoothedElementAbundance", FLOAT, chemistry_element_count,
-      UNIT_CONV_NO_UNITS, sparts, chemistry_data.smoothed_metal_mass_fraction);
+      UNIT_CONV_NO_UNITS, 0.f, sparts,
+      chemistry_data.smoothed_metal_mass_fraction,
+      "Smoothed fractions of the particles' masses that are "
+      "in the given element");
 
-  list[2] =
-      io_make_output_field("Metallicity", FLOAT, 1, UNIT_CONV_NO_UNITS, sparts,
-                           chemistry_data.metal_mass_fraction_total);
+  list[2] = io_make_output_field(
+      "MetalMassFractions", FLOAT, 1, UNIT_CONV_NO_UNITS, 0.f, sparts,
+      chemistry_data.metal_mass_fraction_total,
+      "Fractions of the particles' masses that are in metals");
 
   list[3] = io_make_output_field(
-      "SmoothedMetallicity", FLOAT, 1, UNIT_CONV_NO_UNITS, sparts,
-      chemistry_data.smoothed_metal_mass_fraction_total);
-
-  list[4] = io_make_output_field("TotalMassFromSNIa", FLOAT, 1, UNIT_CONV_MASS,
-                                 sparts, chemistry_data.mass_from_SNIa);
-
-  list[5] = io_make_output_field("MetalMassFracFromSNIa", FLOAT, 1,
-                                 UNIT_CONV_NO_UNITS, sparts,
-                                 chemistry_data.metal_mass_fraction_from_SNIa);
-
-  list[6] = io_make_output_field("TotalMassFromAGB", FLOAT, 1, UNIT_CONV_MASS,
-                                 sparts, chemistry_data.mass_from_AGB);
-
-  list[7] =
-      io_make_output_field("MetalMassFracFromAGB", FLOAT, 1, UNIT_CONV_NO_UNITS,
-                           sparts, chemistry_data.metal_mass_fraction_from_AGB);
-
-  list[8] = io_make_output_field("TotalMassFromSNII", FLOAT, 1, UNIT_CONV_MASS,
-                                 sparts, chemistry_data.mass_from_SNII);
-
-  list[9] = io_make_output_field("MetalMassFracFromSNII", FLOAT, 1,
-                                 UNIT_CONV_NO_UNITS, sparts,
-                                 chemistry_data.metal_mass_fraction_from_SNII);
-
-  list[10] =
-      io_make_output_field("IronMassFracFromSNIa", FLOAT, 1, UNIT_CONV_NO_UNITS,
-                           sparts, chemistry_data.iron_mass_fraction_from_SNIa);
+      "SmoothedMetalMassFractions", FLOAT, 1, UNIT_CONV_NO_UNITS, 0.f, sparts,
+      chemistry_data.smoothed_metal_mass_fraction_total,
+      "Smoothed fractions of the particles masses that are in metals");
+
+  list[4] = io_make_output_field(
+      "MassesFromSNIa", FLOAT, 1, UNIT_CONV_MASS, 0.f, sparts,
+      chemistry_data.mass_from_SNIa,
+      "Masses of gas that have been produced by SNIa stars");
+
+  list[5] = io_make_output_field("MetalMassFractionsFromSNIa", FLOAT, 1,
+                                 UNIT_CONV_NO_UNITS, 0.f, sparts,
+                                 chemistry_data.metal_mass_fraction_from_SNIa,
+                                 "Fractions of the particles' masses that are "
+                                 "in metals produced by SNIa stars");
+
+  list[6] = io_make_output_field(
+      "MassesFromAGB", FLOAT, 1, UNIT_CONV_MASS, 0.f, sparts,
+      chemistry_data.mass_from_AGB,
+      "Masses of gas that have been produced by AGN stars");
+
+  list[7] = io_make_output_field("MetalMassFractionsFromAGB", FLOAT, 1,
+                                 UNIT_CONV_NO_UNITS, 0., sparts,
+                                 chemistry_data.metal_mass_fraction_from_AGB,
+                                 "Fractions of the particles' masses that are "
+                                 "in metals produced by AGB stars");
+
+  list[8] = io_make_output_field(
+      "MassesFromSNII", FLOAT, 1, UNIT_CONV_MASS, 0.f, sparts,
+      chemistry_data.mass_from_SNII,
+      "Masses of gas that have been produced by SNII stars");
+
+  list[9] = io_make_output_field("MetalMassFractionsFromSNII", FLOAT, 1,
+                                 UNIT_CONV_NO_UNITS, 0.f, sparts,
+                                 chemistry_data.metal_mass_fraction_from_SNII,
+                                 "Fractions of the particles' masses that are "
+                                 "in metals produced by SNII stars");
+
+  list[10] = io_make_output_field("IronMassFractionsFromSNIa", FLOAT, 1,
+                                  UNIT_CONV_NO_UNITS, 0.f, sparts,
+                                  chemistry_data.iron_mass_fraction_from_SNIa,
+                                  "Fractions of the particles' masses that are "
+                                  "in iron produced by SNIa stars");
 
   list[11] = io_make_output_field(
-      "SmoothedIronMassFracFromSNIa", FLOAT, 1, UNIT_CONV_NO_UNITS, sparts,
-      chemistry_data.smoothed_iron_mass_fraction_from_SNIa);
+      "SmoothedIronMassFractionsFromSNIa", FLOAT, 1, UNIT_CONV_NO_UNITS, 0.f,
+      sparts, chemistry_data.smoothed_iron_mass_fraction_from_SNIa,
+      "Smoothed fractions of the particles' masses that are "
+      "in iron produced by SNIa stars");
 
   return 12;
 }
diff --git a/src/chemistry/GEAR/chemistry_io.h b/src/chemistry/GEAR/chemistry_io.h
index f4084463dd5468bc2a3fe4ad59f9caab168fba23..de70b99340b800fbb801eceb501c743a2a40e0f5 100644
--- a/src/chemistry/GEAR/chemistry_io.h
+++ b/src/chemistry/GEAR/chemistry_io.h
@@ -73,10 +73,11 @@ INLINE static int chemistry_write_particles(const struct part* parts,
                                             struct io_props* list) {
 
   /* List what we want to write */
-  list[0] = io_make_output_field(
-      "SmoothedElementAbundances", FLOAT, chemistry_element_count,
-      UNIT_CONV_NO_UNITS, 0.f, parts, chemistry_data.smoothed_metal_mass_fraction,
-				 "Element abundances smoothed over the neighbors");
+  list[0] =
+      io_make_output_field("SmoothedElementAbundances", FLOAT,
+                           chemistry_element_count, UNIT_CONV_NO_UNITS, 0.f,
+                           parts, chemistry_data.smoothed_metal_mass_fraction,
+                           "Element abundances smoothed over the neighbors");
 
   list[1] = io_make_output_field("Z", FLOAT, 1, UNIT_CONV_NO_UNITS, 0.f, parts,
                                  chemistry_data.Z, "Temporary field");
@@ -84,7 +85,7 @@ INLINE static int chemistry_write_particles(const struct part* parts,
   list[2] = io_make_output_field("ElementAbundances", FLOAT,
                                  chemistry_element_count, UNIT_CONV_NO_UNITS,
                                  0.f, parts, chemistry_data.metal_mass_fraction,
-				 "Mass fraction of each element");
+                                 "Mass fraction of each element");
 
   return 3;
 }
@@ -101,18 +102,19 @@ INLINE static int chemistry_write_sparticles(const struct spart* sparts,
                                              struct io_props* list) {
 
   /* List what we want to write */
-  list[0] = io_make_output_field(
-      "SmoothedElementAbundances", FLOAT, chemistry_element_count,
-      UNIT_CONV_NO_UNITS, 0.f, sparts, chemistry_data.smoothed_metal_mass_fraction,
-				 "Element abundances smoothed over the neighbors");
+  list[0] =
+      io_make_output_field("SmoothedElementAbundances", FLOAT,
+                           chemistry_element_count, UNIT_CONV_NO_UNITS, 0.f,
+                           sparts, chemistry_data.smoothed_metal_mass_fraction,
+                           "Element abundances smoothed over the neighbors");
 
   list[1] = io_make_output_field("Z", FLOAT, 1, UNIT_CONV_NO_UNITS, 0.f, sparts,
                                  chemistry_data.Z, "Temporary field");
 
-  list[2] = io_make_output_field("ElementAbundance", FLOAT,
-                                 chemistry_element_count, UNIT_CONV_NO_UNITS,
-                                 0.f, sparts, chemistry_data.metal_mass_fraction,
-				 "Mass fraction of each element");
+  list[2] = io_make_output_field(
+      "ElementAbundance", FLOAT, chemistry_element_count, UNIT_CONV_NO_UNITS,
+      0.f, sparts, chemistry_data.metal_mass_fraction,
+      "Mass fraction of each element");
 
   return 3;
 }
diff --git a/src/common_io.c b/src/common_io.c
index 1016980740114ef3ea14388d0d5292ea31eac1de..96a06acd19fe4ceadb2648e9b323dd21e17247fd 100644
--- a/src/common_io.c
+++ b/src/common_io.c
@@ -24,10 +24,12 @@
 /* This object's header. */
 #include "common_io.h"
 
+/* First common header */
+#include "engine.h"
+
 /* Local includes. */
 #include "black_holes_io.h"
 #include "chemistry_io.h"
-#include "engine.h"
 #include "error.h"
 #include "gravity_io.h"
 #include "hydro.h"
@@ -1704,6 +1706,7 @@ void io_check_output_fields(const struct swift_params* params,
   struct xpart xp;
   struct spart sp;
   struct gpart gp;
+  struct bpart bp;
 
   /* Copy N_total to array with length == 6 */
   const long long nr_total[swift_type_count] = {N_total[0], N_total[1], 0,
@@ -1731,7 +1734,11 @@ void io_check_output_fields(const struct swift_params* params,
         break;
 
       case swift_type_stars:
-        stars_write_particles(&sp, list, &num_fields);
+        stars_write_particles(&sp, list, &num_fields, 1);
+        break;
+
+      case swift_type_black_hole:
+        black_holes_write_particles(&bp, list, &num_fields, 1);
         break;
 
       default:
@@ -1821,11 +1828,11 @@ void io_write_output_field_parameter(const char* filename) {
         break;
 
       case swift_type_stars:
-        stars_write_particles(NULL, list, &num_fields);
+        stars_write_particles(NULL, list, &num_fields, 1);
         break;
 
       case swift_type_black_hole:
-        black_holes_write_particles(NULL, list, &num_fields);
+        black_holes_write_particles(NULL, list, &num_fields, 1);
         break;
 
       default:
diff --git a/src/cooling/EAGLE/cooling_io.h b/src/cooling/EAGLE/cooling_io.h
index 5508153afc094d84383893f55ac0362a6d427b24..a57e4451d8e350e1f16a64318b510495ed77e811 100644
--- a/src/cooling/EAGLE/cooling_io.h
+++ b/src/cooling/EAGLE/cooling_io.h
@@ -63,9 +63,9 @@ __attribute__((always_inline)) INLINE static int cooling_write_particles(
     const struct part* parts, const struct xpart* xparts, struct io_props* list,
     const struct cooling_function_data* cooling) {
 
-  list[0] = io_make_output_field_convert_part("Temperature", FLOAT, 1,
-                                              UNIT_CONV_TEMPERATURE, parts,
-                                              xparts, convert_part_T);
+  list[0] = io_make_output_field_convert_part(
+      "Temperatures", FLOAT, 1, UNIT_CONV_TEMPERATURE, 0.f, parts, xparts,
+      convert_part_T, "Temperatures of the gas particles");
   return 1;
 }
 
diff --git a/src/cooling/grackle/cooling_io.h b/src/cooling/grackle/cooling_io.h
index c1472c8371f4d67c1c369793bbf8c4dcd859c412..18712f4cf6c74253bfed3bd36f1a4484b9e50943 100644
--- a/src/cooling/grackle/cooling_io.h
+++ b/src/cooling/grackle/cooling_io.h
@@ -63,23 +63,29 @@ __attribute__((always_inline)) INLINE static int cooling_write_particles(
 
 #if COOLING_GRACKLE_MODE >= 1
   /* List what we want to write */
-  list[0] = io_make_output_field("HI", FLOAT, 1, UNIT_CONV_NO_UNITS, 0.f, xparts,
-                                 cooling_data.HI_frac, "HI mass fraction");
+  list[0] =
+      io_make_output_field("HI", FLOAT, 1, UNIT_CONV_NO_UNITS, 0.f, xparts,
+                           cooling_data.HI_frac, "HI mass fraction");
 
-  list[1] = io_make_output_field("HII", FLOAT, 1, UNIT_CONV_NO_UNITS, 0.f, xparts,
-                                 cooling_data.HII_frac, "HII mass fraction");
+  list[1] =
+      io_make_output_field("HII", FLOAT, 1, UNIT_CONV_NO_UNITS, 0.f, xparts,
+                           cooling_data.HII_frac, "HII mass fraction");
 
-  list[2] = io_make_output_field("HeI", FLOAT, 1, UNIT_CONV_NO_UNITS, 0.f, xparts,
-                                 cooling_data.HeI_frac, "HeI mass fraction");
+  list[2] =
+      io_make_output_field("HeI", FLOAT, 1, UNIT_CONV_NO_UNITS, 0.f, xparts,
+                           cooling_data.HeI_frac, "HeI mass fraction");
 
-  list[3] = io_make_output_field("HeII", FLOAT, 1, UNIT_CONV_NO_UNITS, 0.f, xparts,
-                                 cooling_data.HeII_frac, "HeII mass fraction");
+  list[3] =
+      io_make_output_field("HeII", FLOAT, 1, UNIT_CONV_NO_UNITS, 0.f, xparts,
+                           cooling_data.HeII_frac, "HeII mass fraction");
 
-  list[4] = io_make_output_field("HeIII", FLOAT, 1, UNIT_CONV_NO_UNITS, 0.f, xparts,
-                                 cooling_data.HeIII_frac, "HeIII mass fraction");
+  list[4] =
+      io_make_output_field("HeIII", FLOAT, 1, UNIT_CONV_NO_UNITS, 0.f, xparts,
+                           cooling_data.HeIII_frac, "HeIII mass fraction");
 
-  list[5] = io_make_output_field("e", FLOAT, 1, UNIT_CONV_NO_UNITS, 0.f, xparts,
-                                 cooling_data.e_frac, "free electron mass fraction");
+  list[5] =
+      io_make_output_field("e", FLOAT, 1, UNIT_CONV_NO_UNITS, 0.f, xparts,
+                           cooling_data.e_frac, "free electron mass fraction");
 
   num += 6;
 #endif
@@ -87,14 +93,17 @@ __attribute__((always_inline)) INLINE static int cooling_write_particles(
 #if COOLING_GRACKLE_MODE >= 2
   list += num;
 
-  list[0] = io_make_output_field("HM", FLOAT, 1, UNIT_CONV_NO_UNITS, 0.f, xparts,
-                                 cooling_data.HM_frac, "H- mass fraction");
+  list[0] =
+      io_make_output_field("HM", FLOAT, 1, UNIT_CONV_NO_UNITS, 0.f, xparts,
+                           cooling_data.HM_frac, "H- mass fraction");
 
-  list[1] = io_make_output_field("H2I", FLOAT, 1, UNIT_CONV_NO_UNITS, 0.f, xparts,
-                                 cooling_data.H2I_frac, "H2I mass fraction");
+  list[1] =
+      io_make_output_field("H2I", FLOAT, 1, UNIT_CONV_NO_UNITS, 0.f, xparts,
+                           cooling_data.H2I_frac, "H2I mass fraction");
 
-  list[2] = io_make_output_field("H2II", FLOAT, 1, UNIT_CONV_NO_UNITS, 0.f, xparts,
-                                 cooling_data.H2II_frac, "H2II mass fraction");
+  list[2] =
+      io_make_output_field("H2II", FLOAT, 1, UNIT_CONV_NO_UNITS, 0.f, xparts,
+                           cooling_data.H2II_frac, "H2II mass fraction");
 
   num += 3;
 #endif
@@ -102,14 +111,17 @@ __attribute__((always_inline)) INLINE static int cooling_write_particles(
 #if COOLING_GRACKLE_MODE >= 3
   list += num;
 
-  list[0] = io_make_output_field("DI", FLOAT, 1, UNIT_CONV_NO_UNITS, 0.f, xparts,
-                                 cooling_data.DI_frac, "DI mass fraction");
+  list[0] =
+      io_make_output_field("DI", FLOAT, 1, UNIT_CONV_NO_UNITS, 0.f, xparts,
+                           cooling_data.DI_frac, "DI mass fraction");
 
-  list[1] = io_make_output_field("DII", FLOAT, 1, UNIT_CONV_NO_UNITS, 0.f, xparts,
-                                 cooling_data.DII_frac, "DII mass fraction");
+  list[1] =
+      io_make_output_field("DII", FLOAT, 1, UNIT_CONV_NO_UNITS, 0.f, xparts,
+                           cooling_data.DII_frac, "DII mass fraction");
 
-  list[2] = io_make_output_field("HDI", FLOAT, 1, UNIT_CONV_NO_UNITS, 0.f, xparts,
-                                 cooling_data.HDI_frac, "HDI mass fraction");
+  list[2] =
+      io_make_output_field("HDI", FLOAT, 1, UNIT_CONV_NO_UNITS, 0.f, xparts,
+                           cooling_data.HDI_frac, "HDI mass fraction");
 
   num += 3;
 #endif
diff --git a/src/gravity/Default/gravity_io.h b/src/gravity/Default/gravity_io.h
index 47d3e4d2f37a9be9c284a8eb61310fd4084b9749..736cc259730063ba9e2eae98aba5d2c11c5f26f5 100644
--- a/src/gravity/Default/gravity_io.h
+++ b/src/gravity/Default/gravity_io.h
@@ -109,17 +109,23 @@ INLINE static void darkmatter_write_particles(const struct gpart* gparts,
   /* List what we want to write */
   list[0] = io_make_output_field_convert_gpart(
       "Coordinates", DOUBLE, 3, UNIT_CONV_LENGTH, 1.f, gparts,
-      convert_gpart_pos, "Position of the particles");
+      convert_gpart_pos, "Co-moving position of the particles");
+
   list[1] = io_make_output_field_convert_gpart(
       "Velocities", FLOAT, 3, UNIT_CONV_SPEED, 0.f, gparts, convert_gpart_vel,
-      "Peculiar velocites of the particles");
+      "Peculiar velocities of the stars. This is a * dx/dt where x is the "
+      "co-moving position of the particles.");
+
   list[2] = io_make_output_field("Masses", FLOAT, 1, UNIT_CONV_MASS, 0.f,
-                                 gparts, mass, "");
-  list[3] =
-      io_make_output_field("ParticleIDs", ULONGLONG, 1, UNIT_CONV_NO_UNITS, 0.f,
-                           gparts, id_or_neg_offset, "");
-  list[4] = io_make_output_field("GroupIDs", INT, 1, UNIT_CONV_NO_UNITS, 0.f,
-                                 gparts, group_id, "");
+                                 gparts, mass, "Masses of the particles");
+
+  list[3] = io_make_output_field(
+      "ParticleIDs", ULONGLONG, 1, UNIT_CONV_NO_UNITS, 0.f, gparts,
+      id_or_neg_offset, "Unique ID of the particles");
+
+  list[4] = io_make_output_field(
+      "GroupIDs", INT, 1, UNIT_CONV_NO_UNITS, 0.f, gparts, group_id,
+      "Unique ID of the group to which the particles belong");
 }
 
 #endif /* SWIFT_DEFAULT_GRAVITY_IO_H */
diff --git a/src/hydro/Minimal/hydro_io.h b/src/hydro/Minimal/hydro_io.h
index 8137a3e2c8912ed763c6b44267d1770107746eaa..2bebadda0ed65107de8c6fd637c9682af3d5e17c 100644
--- a/src/hydro/Minimal/hydro_io.h
+++ b/src/hydro/Minimal/hydro_io.h
@@ -163,29 +163,40 @@ INLINE static void hydro_write_particles(const struct part* parts,
   list[0] = io_make_output_field_convert_part(
       "Coordinates", DOUBLE, 3, UNIT_CONV_LENGTH, 1.f, parts, xparts,
       convert_part_pos, "Co-moving position of the particles");
+
   list[1] = io_make_output_field_convert_part(
       "Velocities", FLOAT, 3, UNIT_CONV_SPEED, 0.f, parts, xparts,
-      convert_part_vel, "Peculiar velocites of the particles");
+      convert_part_vel,
+      "Peculiar velocities of the stars. This is a * dx/dt where x is the "
+      "co-moving position of the particles.");
+
   list[2] = io_make_output_field("Masses", FLOAT, 1, UNIT_CONV_MASS, 0.f, parts,
                                  mass, "Masses of the particles");
+
   list[3] = io_make_output_field(
-      "SmoothingLength", FLOAT, 1, UNIT_CONV_LENGTH, 1.f, parts, h,
-      "Smoothing lengths (FWHM of the kernel) of the particles");
+      "SmoothingLengths", FLOAT, 1, UNIT_CONV_LENGTH, 1.f, parts, h,
+      "Co-moving smoothing lengths (FWHM of the kernel) of the particles");
+
   list[4] = io_make_output_field(
       "InternalEnergy", FLOAT, 1, UNIT_CONV_ENERGY_PER_UNIT_MASS,
-      3. * hydro_gamma_minus_one, parts, u, "Thermal energy per unit mass.");
+      3. * hydro_gamma_minus_one, parts, u,
+      "Co-moving thermal energies per unit mass of the particles");
+
   list[5] =
       io_make_output_field("ParticleIDs", ULONGLONG, 1, UNIT_CONV_NO_UNITS, 0.f,
                            parts, id, "Unique ID of the particles");
-  list[6] =
-      io_make_output_field("Density", FLOAT, 1, UNIT_CONV_DENSITY, -3.f, parts,
-                           rho, "Physical mass density of the particles");
+
+  list[6] = io_make_output_field("Densities", FLOAT, 1, UNIT_CONV_DENSITY, -3.f,
+                                 parts, rho,
+                                 "Co-moving mass densities of the particles");
+
   list[7] = io_make_output_field_convert_part(
-      "Entropy", FLOAT, 1, UNIT_CONV_ENTROPY_PER_UNIT_MASS, 0.f, parts, xparts,
-      convert_S, "Entropy per unit mass of the particles");
+      "Entropies", FLOAT, 1, UNIT_CONV_ENTROPY_PER_UNIT_MASS, 0.f, parts,
+      xparts, convert_S, "Co-moving Entropies per unit mass of the particles");
+
   list[8] = io_make_output_field_convert_part(
-      "Pressure", FLOAT, 1, UNIT_CONV_PRESSURE, 3 * hydro_gamma, parts, xparts,
-      convert_P, "Pressure of the particles");
+      "Pressures", FLOAT, 1, UNIT_CONV_PRESSURE, 3 * hydro_gamma, parts, xparts,
+      convert_P, "Co-moving pressures of the particles");
 
   list[9] = io_make_output_field_convert_part(
       "Potential", FLOAT, 1, UNIT_CONV_POTENTIAL, 1.f, parts, xparts,
diff --git a/src/parallel_io.c b/src/parallel_io.c
index d58ee80814ed2f2641176f2191b9d9443f23459a..407e4778a291de06e5a355712d1ecda982cb5fe4 100644
--- a/src/parallel_io.c
+++ b/src/parallel_io.c
@@ -1221,7 +1221,7 @@ void prepare_file(struct engine* e, const char* baseName, long long N_total[6],
         break;
 
       case swift_type_stars:
-        stars_write_particles(sparts, list, &num_fields);
+        stars_write_particles(sparts, list, &num_fields, with_cosmology);
         num_fields += chemistry_write_sparticles(sparts, list + num_fields);
         num_fields +=
             tracers_write_sparticles(sparts, list + num_fields, with_cosmology);
@@ -1231,7 +1231,7 @@ void prepare_file(struct engine* e, const char* baseName, long long N_total[6],
         break;
 
       case swift_type_black_hole:
-        black_holes_write_particles(bparts, list, &num_fields);
+        black_holes_write_particles(bparts, list, &num_fields, with_cosmology);
         if (with_stf) {
           num_fields += velociraptor_write_bparts(bparts, list + num_fields);
         }
@@ -1635,7 +1635,7 @@ void write_output_parallel(struct engine* e, const char* baseName,
 
           /* No inhibted particles: easy case */
           Nparticles = Nstars;
-          stars_write_particles(sparts, list, &num_fields);
+          stars_write_particles(sparts, list, &num_fields, with_cosmology);
           num_fields += chemistry_write_sparticles(sparts, list + num_fields);
           num_fields += tracers_write_sparticles(sparts, list + num_fields,
                                                  with_cosmology);
@@ -1658,7 +1658,8 @@ void write_output_parallel(struct engine* e, const char* baseName,
                                      Nstars_written);
 
           /* Select the fields to write */
-          stars_write_particles(sparts_written, list, &num_fields);
+          stars_write_particles(sparts_written, list, &num_fields,
+                                with_cosmology);
           num_fields += chemistry_write_sparticles(sparts, list + num_fields);
           num_fields += tracers_write_sparticles(sparts, list + num_fields,
                                                  with_cosmology);
@@ -1674,7 +1675,8 @@ void write_output_parallel(struct engine* e, const char* baseName,
 
           /* No inhibted particles: easy case */
           Nparticles = Nblackholes;
-          black_holes_write_particles(bparts, list, &num_fields);
+          black_holes_write_particles(bparts, list, &num_fields,
+                                      with_cosmology);
           if (with_stf) {
             num_fields += velociraptor_write_bparts(bparts, list + num_fields);
           }
@@ -1694,7 +1696,8 @@ void write_output_parallel(struct engine* e, const char* baseName,
                                      Nblackholes_written);
 
           /* Select the fields to write */
-          black_holes_write_particles(bparts_written, list, &num_fields);
+          black_holes_write_particles(bparts_written, list, &num_fields,
+                                      with_cosmology);
           if (with_stf) {
             num_fields +=
                 velociraptor_write_bparts(bparts_written, list + num_fields);
diff --git a/src/serial_io.c b/src/serial_io.c
index 80eff5adace188cdc5919a52a8b55ebc6cf1dbae..999425c46a7130e2b8244a266becd18b7c076d0f 100644
--- a/src/serial_io.c
+++ b/src/serial_io.c
@@ -1273,7 +1273,7 @@ void write_output_serial(struct engine* e, const char* baseName,
 
               /* No inhibted particles: easy case */
               Nparticles = Nstars;
-              stars_write_particles(sparts, list, &num_fields);
+              stars_write_particles(sparts, list, &num_fields, with_cosmology);
               num_fields +=
                   chemistry_write_sparticles(sparts, list + num_fields);
               num_fields += tracers_write_sparticles(sparts, list + num_fields,
@@ -1298,7 +1298,8 @@ void write_output_serial(struct engine* e, const char* baseName,
                                          Nstars_written);
 
               /* Select the fields to write */
-              stars_write_particles(sparts_written, list, &num_fields);
+              stars_write_particles(sparts_written, list, &num_fields,
+                                    with_cosmology);
               num_fields +=
                   chemistry_write_sparticles(sparts, list + num_fields);
               num_fields += tracers_write_sparticles(sparts, list + num_fields,
@@ -1315,7 +1316,8 @@ void write_output_serial(struct engine* e, const char* baseName,
 
               /* No inhibted particles: easy case */
               Nparticles = Nblackholes;
-              black_holes_write_particles(bparts, list, &num_fields);
+              black_holes_write_particles(bparts, list, &num_fields,
+                                          with_cosmology);
               if (with_stf) {
                 num_fields +=
                     velociraptor_write_bparts(bparts, list + num_fields);
@@ -1336,7 +1338,8 @@ void write_output_serial(struct engine* e, const char* baseName,
                                          Nblackholes_written);
 
               /* Select the fields to write */
-              black_holes_write_particles(bparts_written, list, &num_fields);
+              black_holes_write_particles(bparts_written, list, &num_fields,
+                                          with_cosmology);
               if (with_stf) {
                 num_fields += velociraptor_write_bparts(bparts_written,
                                                         list + num_fields);
diff --git a/src/single_io.c b/src/single_io.c
index 3e6ea192def788a66a9d8c6badf832cd7f3c1e28..0682480d166e476bdbe0687a0ae20db556e0f4d2 100644
--- a/src/single_io.c
+++ b/src/single_io.c
@@ -1077,7 +1077,7 @@ void write_output_single(struct engine* e, const char* baseName,
 
           /* No inhibted particles: easy case */
           N = Nstars;
-          stars_write_particles(sparts, list, &num_fields);
+          stars_write_particles(sparts, list, &num_fields, with_cosmology);
           num_fields += chemistry_write_sparticles(sparts, list + num_fields);
           num_fields += tracers_write_sparticles(sparts, list + num_fields,
                                                  with_cosmology);
@@ -1100,7 +1100,8 @@ void write_output_single(struct engine* e, const char* baseName,
                                      Nstars_written);
 
           /* Select the fields to write */
-          stars_write_particles(sparts_written, list, &num_fields);
+          stars_write_particles(sparts_written, list, &num_fields,
+                                with_cosmology);
           num_fields +=
               chemistry_write_sparticles(sparts_written, list + num_fields);
           num_fields += tracers_write_sparticles(
@@ -1117,7 +1118,8 @@ void write_output_single(struct engine* e, const char* baseName,
 
           /* No inhibted particles: easy case */
           N = Nblackholes;
-          black_holes_write_particles(bparts, list, &num_fields);
+          black_holes_write_particles(bparts, list, &num_fields,
+                                      with_cosmology);
           if (with_stf) {
             num_fields += velociraptor_write_bparts(bparts, list + num_fields);
           }
@@ -1137,7 +1139,8 @@ void write_output_single(struct engine* e, const char* baseName,
                                      Nblackholes_written);
 
           /* Select the fields to write */
-          black_holes_write_particles(bparts_written, list, &num_fields);
+          black_holes_write_particles(bparts_written, list, &num_fields,
+                                      with_cosmology);
           if (with_stf) {
             num_fields +=
                 velociraptor_write_bparts(bparts_written, list + num_fields);
diff --git a/src/star_formation/EAGLE/star_formation_io.h b/src/star_formation/EAGLE/star_formation_io.h
index cee96326e458d0581af6e62e452ac433dcf407bd..f8bf57145955d41f7ae0ecd4141651a24f1c2727 100644
--- a/src/star_formation/EAGLE/star_formation_io.h
+++ b/src/star_formation/EAGLE/star_formation_io.h
@@ -38,8 +38,11 @@ __attribute__((always_inline)) INLINE static int star_formation_write_particles(
     const struct part* parts, const struct xpart* xparts,
     struct io_props* list) {
 
-  list[0] =
-      io_make_output_field("SFR", FLOAT, 1, UNIT_CONV_SFR, xparts, sf_data.SFR);
+  list[0] = io_make_output_field(
+      "StarFormationRates", FLOAT, 1, UNIT_CONV_SFR, 0.f, xparts, sf_data.SFR,
+      "If positive, star formation rates of the particles. If negative, stores "
+      "the last time/scale-factor at which the gas particle was star-forming. "
+      "If zero, the particle was never star-forming.");
 
   return 1;
 }
diff --git a/src/stars/EAGLE/stars_io.h b/src/stars/EAGLE/stars_io.h
index 0c03bee3007066c7c51c7ce0fb3d88d37a1b2ae3..8e3f6056921936fdcf86337af165e130cc6e02eb 100644
--- a/src/stars/EAGLE/stars_io.h
+++ b/src/stars/EAGLE/stars_io.h
@@ -52,6 +52,53 @@ INLINE static void stars_read_particles(struct spart *sparts,
                                 sparts, mass_init);
 }
 
+INLINE static void convert_spart_pos(const struct engine *e,
+                                     const struct spart *sp, double *ret) {
+
+  if (e->s->periodic) {
+    ret[0] = box_wrap(sp->x[0], 0.0, e->s->dim[0]);
+    ret[1] = box_wrap(sp->x[1], 0.0, e->s->dim[1]);
+    ret[2] = box_wrap(sp->x[2], 0.0, e->s->dim[2]);
+  } else {
+    ret[0] = sp->x[0];
+    ret[1] = sp->x[1];
+    ret[2] = sp->x[2];
+  }
+}
+
+INLINE static void convert_spart_vel(const struct engine *e,
+                                     const struct spart *sp, float *ret) {
+
+  const int with_cosmology = (e->policy & engine_policy_cosmology);
+  const struct cosmology *cosmo = e->cosmology;
+  const integertime_t ti_current = e->ti_current;
+  const double time_base = e->time_base;
+
+  const integertime_t ti_beg = get_integer_time_begin(ti_current, sp->time_bin);
+  const integertime_t ti_end = get_integer_time_end(ti_current, sp->time_bin);
+
+  /* Get time-step since the last kick */
+  float dt_kick_grav;
+  if (with_cosmology) {
+    dt_kick_grav = cosmology_get_grav_kick_factor(cosmo, ti_beg, ti_current);
+    dt_kick_grav -=
+        cosmology_get_grav_kick_factor(cosmo, ti_beg, (ti_beg + ti_end) / 2);
+  } else {
+    dt_kick_grav = (ti_current - ((ti_beg + ti_end) / 2)) * time_base;
+  }
+
+  /* Extrapolate the velocites to the current time */
+  const struct gpart *gp = sp->gpart;
+  ret[0] = gp->v_full[0] + gp->a_grav[0] * dt_kick_grav;
+  ret[1] = gp->v_full[1] + gp->a_grav[1] * dt_kick_grav;
+  ret[2] = gp->v_full[2] + gp->a_grav[2] * dt_kick_grav;
+
+  /* Conversion from internal units to peculiar velocities */
+  ret[0] *= cosmo->a_inv;
+  ret[1] *= cosmo->a_inv;
+  ret[2] *= cosmo->a_inv;
+}
+
 /**
  * @brief Specifies which s-particle fields to write to a dataset
  *
@@ -60,34 +107,65 @@ INLINE static void stars_read_particles(struct spart *sparts,
  * @param num_fields The number of i/o fields to write.
  */
 INLINE static void stars_write_particles(const struct spart *sparts,
-                                         struct io_props *list,
-                                         int *num_fields) {
+                                         struct io_props *list, int *num_fields,
+                                         const int with_cosmology) {
 
   /* Say how much we want to write */
   *num_fields = 10;
 
   /* List what we want to write */
-  list[0] = io_make_output_field("Coordinates", DOUBLE, 3, UNIT_CONV_LENGTH,
-                                 sparts, x);
-  list[1] =
-      io_make_output_field("Velocities", FLOAT, 3, UNIT_CONV_SPEED, sparts, v);
-  list[2] =
-      io_make_output_field("Masses", FLOAT, 1, UNIT_CONV_MASS, sparts, mass);
+  list[0] = io_make_output_field_convert_spart(
+      "Coordinates", DOUBLE, 3, UNIT_CONV_LENGTH, 1.f, sparts,
+      convert_spart_pos, "Co-moving position of the particles");
+
+  list[1] = io_make_output_field_convert_spart(
+      "Velocities", FLOAT, 3, UNIT_CONV_SPEED, 0.f, sparts, convert_spart_vel,
+      "Peculiar velocities of the particles. This is a * dx/dt where x is the "
+      "co-moving position of the particles.");
+
+  list[2] = io_make_output_field("Masses", FLOAT, 1, UNIT_CONV_MASS, 0.f,
+                                 sparts, mass,
+                                 "Masses of the particles at the current point "
+                                 "in time (i.e. after stellar losses");
+
   list[3] = io_make_output_field("ParticleIDs", LONGLONG, 1, UNIT_CONV_NO_UNITS,
-                                 sparts, id);
-  list[4] = io_make_output_field("SmoothingLength", FLOAT, 1, UNIT_CONV_LENGTH,
-                                 sparts, h);
-  list[5] = io_make_output_field("BirthDensity", FLOAT, 1, UNIT_CONV_DENSITY,
-                                 sparts, birth_density);
-  list[6] = io_make_output_field("InitialMasses", FLOAT, 1, UNIT_CONV_MASS,
-                                 sparts, mass_init);
-  list[7] = io_make_output_field("BirthTime", FLOAT, 1, UNIT_CONV_TIME, sparts,
-                                 birth_time);
-  list[8] = io_make_output_field("FeedbackEnergyFraction", FLOAT, 1,
-                                 UNIT_CONV_NO_UNITS, sparts, f_E);
+                                 0.f, sparts, id, "Unique ID of the particles");
+
+  list[4] = io_make_output_field(
+      "SmoothingLengths", FLOAT, 1, UNIT_CONV_LENGTH, 1.f, sparts, h,
+      "Co-moving smoothing lengths (FWHM of the kernel) of the particles");
+
+  list[5] = io_make_output_field(
+      "BirthDensities", FLOAT, 1, UNIT_CONV_DENSITY, 0.f, sparts, birth_density,
+      "Physical densities at the time of birth of the gas particles that "
+      "turned into stars (note that "
+      "we store the physical density at the birth redshift, no conversion is "
+      "needed)");
+
+  list[6] = io_make_output_field("InitialMasses", FLOAT, 1, UNIT_CONV_MASS, 0.f,
+                                 sparts, mass_init,
+                                 "Masses of the star particles at birth time");
+
+  if (with_cosmology) {
+    list[7] = io_make_output_field(
+        "BirthScaleFactors", FLOAT, 1, UNIT_CONV_NO_UNITS, 0.f, sparts,
+        birth_scale_factor, "Scale-factors at which the stars were born");
+  } else {
+    list[7] = io_make_output_field("BirthTimes", FLOAT, 1, UNIT_CONV_TIME, 0.f,
+                                   sparts, birth_time,
+                                   "Times at which the stars were born");
+  }
+
+  list[8] = io_make_output_field(
+      "FeedbackEnergyFractions", FLOAT, 1, UNIT_CONV_NO_UNITS, 0.f, sparts, f_E,
+      "Fractions of the canonical feedback energy that was used for the stars' "
+      "SNII feedback events");
+
   list[9] =
-      io_make_output_field("BirthTemperature", FLOAT, 1, UNIT_CONV_TEMPERATURE,
-                           sparts, birth_temperature);
+      io_make_output_field("BirthTemperatures", FLOAT, 1, UNIT_CONV_TEMPERATURE,
+                           0.f, sparts, birth_temperature,
+                           "Temperatures at the time of birth of the gas "
+                           "particles that turned into stars");
 }
 
 /**
diff --git a/src/tracers/EAGLE/tracers_io.h b/src/tracers/EAGLE/tracers_io.h
index 038cc1c8d3f92c2105d5b5c3ead958f60486ce9f..fe650e9a4ced6f06dae8a5133baf5ce766c22e0c 100644
--- a/src/tracers/EAGLE/tracers_io.h
+++ b/src/tracers/EAGLE/tracers_io.h
@@ -54,20 +54,23 @@ __attribute__((always_inline)) INLINE static int tracers_write_particles(
     const struct part* parts, const struct xpart* xparts, struct io_props* list,
     const int with_cosmology) {
 
-  list[0] = io_make_output_field("Maximal Temperature", FLOAT, 1,
-                                 UNIT_CONV_TEMPERATURE, xparts,
-                                 tracers_data.maximum_temperature);
+  list[0] = io_make_output_field(
+      "MaximalTemperatures", FLOAT, 1, UNIT_CONV_TEMPERATURE, 0.f, xparts,
+      tracers_data.maximum_temperature,
+      "Maximal temperatures ever reached by the particles");
 
   if (with_cosmology) {
     list[1] = io_make_output_field(
-        "Maximal Temperature scale-factor", FLOAT, 1, UNIT_CONV_NO_UNITS,
-        xparts, tracers_data.maximum_temperature_scale_factor);
+        "MaximalTemperatureScaleFactors", FLOAT, 1, UNIT_CONV_NO_UNITS, 0.f,
+        xparts, tracers_data.maximum_temperature_scale_factor,
+        "Scale-factors at which the maximal temperature was reached");
 
   } else {
 
-    list[1] = io_make_output_field("Maximal Temperature time", FLOAT, 1,
-                                   UNIT_CONV_NO_UNITS, xparts,
-                                   tracers_data.maximum_temperature_time);
+    list[1] = io_make_output_field(
+        "MaximalTemperatureTimes", FLOAT, 1, UNIT_CONV_NO_UNITS, 0.f, xparts,
+        tracers_data.maximum_temperature_time,
+        "Times at which the maximal temperature was reached");
   }
 
   return 2;
@@ -77,22 +80,27 @@ __attribute__((always_inline)) INLINE static int tracers_write_sparticles(
     const struct spart* sparts, struct io_props* list,
     const int with_cosmology) {
 
-  list[0] = io_make_output_field("Maximal Temperature", FLOAT, 1,
-                                 UNIT_CONV_TEMPERATURE, sparts,
-                                 tracers_data.maximum_temperature);
+  list[0] = io_make_output_field(
+      "MaximalTemperatures", FLOAT, 1, UNIT_CONV_TEMPERATURE, 0.f, sparts,
+      tracers_data.maximum_temperature,
+      "Maximal temperatures ever reached by the particles before they got "
+      "converted to stars");
 
   if (with_cosmology) {
     list[1] = io_make_output_field(
-        "Maximal Temperature scale-factor", FLOAT, 1, UNIT_CONV_NO_UNITS,
-        sparts, tracers_data.maximum_temperature_scale_factor);
+        "MaximalTemperatureScaleFactors", FLOAT, 1, UNIT_CONV_NO_UNITS, 0.f,
+        sparts, tracers_data.maximum_temperature_scale_factor,
+        "Scale-factors at which the maximal temperature was reached");
 
   } else {
 
-    list[1] = io_make_output_field("Maximal Temperature time", FLOAT, 1,
-                                   UNIT_CONV_NO_UNITS, sparts,
-                                   tracers_data.maximum_temperature_time);
+    list[1] = io_make_output_field(
+        "Maximal Temperature time", FLOAT, 1, UNIT_CONV_NO_UNITS, 0.f, sparts,
+        tracers_data.maximum_temperature_time,
+        "Times at which the maximal temperature was reached");
   }
 
   return 2;
 }
+
 #endif /* SWIFT_TRACERS_EAGLE_IO_H */