diff --git a/doc/RTD/source/Snapshots/index.rst b/doc/RTD/source/Snapshots/index.rst
index 19abf8b19ddd8fc08d3c6798f7d134d254c02c0c..d1700edd158ac5d918ebb53a5afe69be32a36d63 100644
--- a/doc/RTD/source/Snapshots/index.rst
+++ b/doc/RTD/source/Snapshots/index.rst
@@ -267,7 +267,9 @@ part designed for users to directly read and in part for machine
 reading of the information. Each field contains the exponent of the
 scale-factor, reduced Hubble constant [#f2]_ and each of the 5 base units
 that is required to convert the field values to physical CGS
-units. These fields are:
+units. The base assumption is that all fields are written in the
+co-moving frame (see below for exceptions).
+These fields are:
 
 +----------------------+---------------------------------------+
 | Meta-data field name | Description                           |
@@ -324,6 +326,12 @@ case of the densities and assuming the usual system of units
 In the case of a non-cosmological simulation, these two expressions
 are identical since :math:`a=1`.
 
+In some special cases, the fields cannot be meaningfully expressed as
+co-moving quantities. In these exceptional circumstances, we set the
+value of the attribute ``Value stored as physical`` to ``1``. And we
+additionally set the attribute ``Property can be converted to
+comoving`` to ``0``.
+
 Particle splitting metadata
 ---------------------------
 
diff --git a/src/black_holes/Default/black_holes_io.h b/src/black_holes/Default/black_holes_io.h
index 43aa06f1fb0f6cb4206d72a2edbd4af47e018547..e74c44e6dd22b377d01db41acaefd84b5cdeaa5f 100644
--- a/src/black_holes/Default/black_holes_io.h
+++ b/src/black_holes/Default/black_holes_io.h
@@ -144,9 +144,9 @@ INLINE static void black_holes_write_particles(const struct bpart* bparts,
   list[2] = io_make_output_field("Masses", FLOAT, 1, UNIT_CONV_MASS, 0.f,
                                  bparts, mass, "Masses of the particles");
 
-  list[3] =
-      io_make_output_field("ParticleIDs", ULONGLONG, 1, UNIT_CONV_NO_UNITS, 0.f,
-                           bparts, id, "Unique ID of the particles");
+  list[3] = io_make_physical_output_field(
+      "ParticleIDs", ULONGLONG, 1, UNIT_CONV_NO_UNITS, 0.f, bparts, id,
+      /*can convert to comoving=*/0, "Unique ID of the particles");
 
   list[4] = io_make_output_field(
       "SmoothingLengths", FLOAT, 1, UNIT_CONV_LENGTH, 1.f, bparts, h,
diff --git a/src/black_holes/EAGLE/black_holes_io.h b/src/black_holes/EAGLE/black_holes_io.h
index 9a588be7eae825f4b9b7894755dfe08574defe31..5ab59fc77d73f8b74c216bccab376748c620324e 100644
--- a/src/black_holes/EAGLE/black_holes_io.h
+++ b/src/black_holes/EAGLE/black_holes_io.h
@@ -185,9 +185,9 @@ INLINE static void black_holes_write_particles(const struct bpart* bparts,
       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", ULONGLONG, 1, UNIT_CONV_NO_UNITS, 0.f,
-                           bparts, id, "Unique ID of the particles");
+  list[3] = io_make_physical_output_field(
+      "ParticleIDs", ULONGLONG, 1, UNIT_CONV_NO_UNITS, 0.f, bparts, id,
+      /*can convert to comoving=*/0, "Unique ID of the particles");
 
   list[4] = io_make_output_field(
       "SmoothingLengths", FLOAT, 1, UNIT_CONV_LENGTH, 1.f, bparts, h,
@@ -198,9 +198,10 @@ INLINE static void black_holes_write_particles(const struct bpart* bparts,
                                  "Subgrid masses of the particles");
 
   if (with_cosmology) {
-    list[6] = io_make_output_field(
+    list[6] = io_make_physical_output_field(
         "FormationScaleFactors", FLOAT, 1, UNIT_CONV_NO_UNITS, 0.f, bparts,
-        formation_scale_factor, "Scale-factors at which the BHs were formed");
+        formation_scale_factor, /*can convert to comoving=*/0,
+        "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,
@@ -233,22 +234,23 @@ INLINE static void black_holes_write_particles(const struct bpart* bparts,
       "birth. This does not include any mass accreted onto any merged black "
       "holes.");
 
-  list[12] = io_make_output_field(
+  list[12] = io_make_physical_output_field(
       "CumulativeNumberOfSeeds", INT, 1, UNIT_CONV_NO_UNITS, 0.f, bparts,
-      cumulative_number_seeds,
+      cumulative_number_seeds, /*can convert to comoving=*/0,
       "Total number of BH seeds that have merged into this black hole");
 
-  list[13] =
-      io_make_output_field("NumberOfMergers", INT, 1, UNIT_CONV_NO_UNITS, 0.f,
-                           bparts, number_of_mergers,
-                           "Number of mergers the black holes went through. "
-                           "This does not include the number of mergers "
-                           "accumulated by any merged black hole.");
+  list[13] = io_make_physical_output_field(
+      "NumberOfMergers", INT, 1, UNIT_CONV_NO_UNITS, 0.f, bparts,
+      number_of_mergers, /*can convert to comoving=*/1,
+      "Number of mergers the black holes went through. "
+      "This does not include the number of mergers "
+      "accumulated by any merged black hole.");
 
   if (with_cosmology) {
-    list[14] = io_make_output_field(
+    list[14] = io_make_physical_output_field(
         "LastHighEddingtonFractionScaleFactors", FLOAT, 1, UNIT_CONV_NO_UNITS,
         0.f, bparts, last_high_Eddington_fraction_scale_factor,
+        /*can convert to comoving=*/0,
         "Scale-factors at which the black holes last reached a large Eddington "
         "ratio. -1 if never reached.");
   } else {
@@ -260,9 +262,9 @@ INLINE static void black_holes_write_particles(const struct bpart* bparts,
   }
 
   if (with_cosmology) {
-    list[15] = io_make_output_field(
+    list[15] = io_make_physical_output_field(
         "LastMinorMergerScaleFactors", FLOAT, 1, UNIT_CONV_NO_UNITS, 0.f,
-        bparts, last_minor_merger_scale_factor,
+        bparts, last_minor_merger_scale_factor, /*can convert to comoving=*/0,
         "Scale-factors at which the black holes last had a minor merger.");
   } else {
     list[15] = io_make_output_field(
@@ -272,9 +274,9 @@ INLINE static void black_holes_write_particles(const struct bpart* bparts,
   }
 
   if (with_cosmology) {
-    list[16] = io_make_output_field(
+    list[16] = io_make_physical_output_field(
         "LastMajorMergerScaleFactors", FLOAT, 1, UNIT_CONV_NO_UNITS, 0.f,
-        bparts, last_major_merger_scale_factor,
+        bparts, last_major_merger_scale_factor, /*can convert to comoving=*/0,
         "Scale-factors at which the black holes last had a major merger.");
   } else {
     list[16] = io_make_output_field(
@@ -308,30 +310,30 @@ INLINE static void black_holes_write_particles(const struct bpart* bparts,
       io_make_output_field("TimeBins", CHAR, 1, UNIT_CONV_NO_UNITS, 0.f, bparts,
                            time_bin, "Time-bins of the particles");
 
-  list[21] = io_make_output_field(
+  list[21] = io_make_physical_output_field(
       "NumberOfSwallows", INT, 1, UNIT_CONV_NO_UNITS, 0.f, bparts,
-      number_of_gas_swallows,
+      number_of_gas_swallows, /*can convert to comoving=*/1,
       "Number of gas particles the black holes have swallowed. "
       "This includes the particles swallowed by any of the black holes that "
       "merged into this one.");
 
-  list[22] = io_make_output_field(
+  list[22] = io_make_physical_output_field(
       "NumberOfDirectSwallows", INT, 1, UNIT_CONV_NO_UNITS, 0.f, bparts,
-      number_of_direct_gas_swallows,
+      number_of_direct_gas_swallows, /*can convert to comoving=*/1,
       "Number of gas particles the black holes have swallowed. "
       "This does not include any particles swallowed by any of the black holes "
       "that merged into this one.");
 
-  list[23] = io_make_output_field(
+  list[23] = io_make_physical_output_field(
       "NumberOfRepositions", INT, 1, UNIT_CONV_NO_UNITS, 0.f, bparts,
-      number_of_repositions,
+      number_of_repositions, /*can convert to comoving=*/1,
       "Number of repositioning events the black holes went through. This does "
       "not include the number of reposition events accumulated by any merged "
       "black holes.");
 
-  list[24] = io_make_output_field(
+  list[24] = io_make_physical_output_field(
       "NumberOfRepositionAttempts", INT, 1, UNIT_CONV_NO_UNITS, 0.f, bparts,
-      number_of_reposition_attempts,
+      number_of_reposition_attempts, /*can convert to comoving=*/1,
       "Number of time steps in which the black holes had an eligible particle "
       "to reposition to. They may or may not have ended up moving there, "
       "depending on their subgrid mass and on whether these particles were at "
@@ -339,9 +341,9 @@ INLINE static void black_holes_write_particles(const struct bpart* bparts,
       "not include attempted repositioning events accumulated by any merged "
       "black holes.");
 
-  list[25] = io_make_output_field(
+  list[25] = io_make_physical_output_field(
       "NumberOfTimeSteps", INT, 1, UNIT_CONV_NO_UNITS, 0.f, bparts,
-      number_of_time_steps,
+      number_of_time_steps, /*can convert to comoving=*/0,
       "Total number of time steps at which the black holes were active.");
 
   list[26] = io_make_output_field(
@@ -360,9 +362,9 @@ INLINE static void black_holes_write_particles(const struct bpart* bparts,
       sound_speed_subgrid_gas,
       "Physical subgrid sound-speeds used in the subgrid-Bondi model.");
 
-  list[29] = io_make_output_field(
+  list[29] = io_make_physical_output_field(
       "BirthGasDensities", FLOAT, 1, UNIT_CONV_DENSITY, 0.f, bparts,
-      formation_gas_density,
+      formation_gas_density, /*can convert to comoving=*/0,
       "Physical densities of the converted part at the time of birth. "
       "We store the physical density at the birth redshift, no conversion is "
       "needed.");
@@ -373,9 +375,9 @@ INLINE static void black_holes_write_particles(const struct bpart* bparts,
       "Physical angular momenta that the black holes have accumulated through "
       "subgrid accretion.");
 
-  list[31] = io_make_output_field(
+  list[31] = io_make_physical_output_field(
       "NumberOfGasNeighbours", INT, 1, UNIT_CONV_NO_UNITS, 0.f, bparts,
-      num_ngbs,
+      num_ngbs, /*can convert to comoving=*/0,
       "Integer number of gas neighbour particles within the black hole "
       "kernels.");
 
@@ -392,23 +394,23 @@ INLINE static void black_holes_write_particles(const struct bpart* bparts,
       "This is 0 for black holes that have never repositioned, or if the "
       "simulation has been run without prescribed repositioning speed.");
 
-  list[34] = io_make_output_field(
+  list[34] = io_make_physical_output_field(
       "NumberOfHeatingEvents", INT, 1, UNIT_CONV_NO_UNITS, 0.f, bparts,
-      AGN_number_of_energy_injections,
+      AGN_number_of_energy_injections, /*can convert to comoving=*/1,
       "Integer number of (thermal) energy injections the black hole has had "
       "so far. This counts each heated gas particle separately, and so can "
       "increase by more than one during a single time step.");
 
-  list[35] = io_make_output_field(
+  list[35] = io_make_physical_output_field(
       "NumberOfAGNEvents", INT, 1, UNIT_CONV_NO_UNITS, 0.f, bparts,
-      AGN_number_of_AGN_events,
+      AGN_number_of_AGN_events, /*can convert to comoving=*/1,
       "Integer number of AGN events the black hole has had so far"
       " (the number of time steps in which the BH did AGN feedback).");
 
   if (with_cosmology) {
-    list[36] = io_make_output_field(
+    list[36] = io_make_physical_output_field(
         "LastAGNFeedbackScaleFactors", FLOAT, 1, UNIT_CONV_NO_UNITS, 0.f,
-        bparts, last_AGN_event_scale_factor,
+        bparts, last_AGN_event_scale_factor, /*can convert to comoving=*/0,
         "Scale-factors at which the black holes last had an AGN event.");
   } else {
     list[36] = io_make_output_field(
diff --git a/src/black_holes/SPIN_JET/black_holes_io.h b/src/black_holes/SPIN_JET/black_holes_io.h
index 825ddddf08c3075e6ee469a66a27859102476115..68273cf59989ebfa53960ee71b77ac6413c37001 100644
--- a/src/black_holes/SPIN_JET/black_holes_io.h
+++ b/src/black_holes/SPIN_JET/black_holes_io.h
@@ -192,9 +192,9 @@ INLINE static void black_holes_write_particles(const struct bpart* bparts,
       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", ULONGLONG, 1, UNIT_CONV_NO_UNITS, 0.f,
-                           bparts, id, "Unique ID of the particles");
+  list[3] = io_make_physical_output_field(
+      "ParticleIDs", ULONGLONG, 1, UNIT_CONV_NO_UNITS, 0.f, bparts, id,
+      /*can convert to comoving=*/0, "Unique ID of the particles");
 
   list[4] = io_make_output_field(
       "SmoothingLengths", FLOAT, 1, UNIT_CONV_LENGTH, 1.f, bparts, h,
@@ -205,9 +205,10 @@ INLINE static void black_holes_write_particles(const struct bpart* bparts,
                                  "Subgrid masses of the particles");
 
   if (with_cosmology) {
-    list[6] = io_make_output_field(
+    list[6] = io_make_physical_output_field(
         "FormationScaleFactors", FLOAT, 1, UNIT_CONV_NO_UNITS, 0.f, bparts,
-        formation_scale_factor, "Scale-factors at which the BHs were formed");
+        formation_scale_factor, /*can convert to comoving=*/0,
+        "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,
@@ -240,22 +241,23 @@ INLINE static void black_holes_write_particles(const struct bpart* bparts,
       "birth. This does not include any mass accreted onto any merged black "
       "holes.");
 
-  list[12] = io_make_output_field(
+  list[12] = io_make_physical_output_field(
       "CumulativeNumberOfSeeds", INT, 1, UNIT_CONV_NO_UNITS, 0.f, bparts,
-      cumulative_number_seeds,
+      cumulative_number_seeds, /*can convert to comoving=*/1,
       "Total number of BH seeds that have merged into this black hole");
 
-  list[13] =
-      io_make_output_field("NumberOfMergers", INT, 1, UNIT_CONV_NO_UNITS, 0.f,
-                           bparts, number_of_mergers,
-                           "Number of mergers the black holes went through. "
-                           "This does not include the number of mergers "
-                           "accumulated by any merged black hole.");
+  list[13] = io_make_physical_output_field(
+      "NumberOfMergers", INT, 1, UNIT_CONV_NO_UNITS, 0.f, bparts,
+      number_of_mergers, /*can convert to comoving=*/1,
+      "Number of mergers the black holes went through. "
+      "This does not include the number of mergers "
+      "accumulated by any merged black hole.");
 
   if (with_cosmology) {
-    list[14] = io_make_output_field(
+    list[14] = io_make_physical_output_field(
         "LastHighEddingtonFractionScaleFactors", FLOAT, 1, UNIT_CONV_NO_UNITS,
         0.f, bparts, last_high_Eddington_fraction_scale_factor,
+        /*can convert to comoving=*/0,
         "Scale-factors at which the black holes last reached a large Eddington "
         "ratio. -1 if never reached.");
   } else {
@@ -267,9 +269,9 @@ INLINE static void black_holes_write_particles(const struct bpart* bparts,
   }
 
   if (with_cosmology) {
-    list[15] = io_make_output_field(
+    list[15] = io_make_physical_output_field(
         "LastMinorMergerScaleFactors", FLOAT, 1, UNIT_CONV_NO_UNITS, 0.f,
-        bparts, last_minor_merger_scale_factor,
+        bparts, last_minor_merger_scale_factor, /*can convert to comoving=*/0,
         "Scale-factors at which the black holes last had a minor merger.");
   } else {
     list[15] = io_make_output_field(
@@ -279,9 +281,9 @@ INLINE static void black_holes_write_particles(const struct bpart* bparts,
   }
 
   if (with_cosmology) {
-    list[16] = io_make_output_field(
+    list[16] = io_make_physical_output_field(
         "LastMajorMergerScaleFactors", FLOAT, 1, UNIT_CONV_NO_UNITS, 0.f,
-        bparts, last_major_merger_scale_factor,
+        bparts, last_major_merger_scale_factor, /*can convert to comoving=*/0,
         "Scale-factors at which the black holes last had a major merger.");
   } else {
     list[16] = io_make_output_field(
@@ -315,30 +317,30 @@ INLINE static void black_holes_write_particles(const struct bpart* bparts,
       io_make_output_field("TimeBins", CHAR, 1, UNIT_CONV_NO_UNITS, 0.f, bparts,
                            time_bin, "Time-bins of the particles");
 
-  list[21] = io_make_output_field(
+  list[21] = io_make_physical_output_field(
       "NumberOfSwallows", INT, 1, UNIT_CONV_NO_UNITS, 0.f, bparts,
-      number_of_gas_swallows,
+      number_of_gas_swallows, /*can convert to comoving=*/1,
       "Number of gas particles the black holes have swallowed. "
       "This includes the particles swallowed by any of the black holes that "
       "merged into this one.");
 
-  list[22] = io_make_output_field(
+  list[22] = io_make_physical_output_field(
       "NumberOfDirectSwallows", INT, 1, UNIT_CONV_NO_UNITS, 0.f, bparts,
-      number_of_direct_gas_swallows,
+      number_of_direct_gas_swallows, /*can convert to comoving=*/1,
       "Number of gas particles the black holes have swallowed. "
       "This does not include any particles swallowed by any of the black holes "
       "that merged into this one.");
 
-  list[23] = io_make_output_field(
+  list[23] = io_make_physical_output_field(
       "NumberOfRepositions", INT, 1, UNIT_CONV_NO_UNITS, 0.f, bparts,
-      number_of_repositions,
+      number_of_repositions, /*can convert to comoving=*/1,
       "Number of repositioning events the black holes went through. This does "
       "not include the number of reposition events accumulated by any merged "
       "black holes.");
 
-  list[24] = io_make_output_field(
+  list[24] = io_make_physical_output_field(
       "NumberOfRepositionAttempts", INT, 1, UNIT_CONV_NO_UNITS, 0.f, bparts,
-      number_of_reposition_attempts,
+      number_of_reposition_attempts, /*can convert to comoving=*/1,
       "Number of time steps in which the black holes had an eligible particle "
       "to reposition to. They may or may not have ended up moving there, "
       "depending on their subgrid mass and on whether these particles were at "
@@ -346,9 +348,9 @@ INLINE static void black_holes_write_particles(const struct bpart* bparts,
       "not include attempted repositioning events accumulated by any merged "
       "black holes.");
 
-  list[25] = io_make_output_field(
+  list[25] = io_make_physical_output_field(
       "NumberOfTimeSteps", INT, 1, UNIT_CONV_NO_UNITS, 0.f, bparts,
-      number_of_time_steps,
+      number_of_time_steps, /*can convert to comoving=*/1,
       "Total number of time steps at which the black holes were active.");
 
   list[26] = io_make_output_field(
@@ -375,29 +377,29 @@ INLINE static void black_holes_write_particles(const struct bpart* bparts,
       "Physical angular momenta that the black holes have accumulated through "
       "subgrid accretion.");
 
-  list[30] = io_make_output_field(
+  list[30] = io_make_physical_output_field(
       "NumberOfGasNeighbours", INT, 1, UNIT_CONV_NO_UNITS, 0.f, bparts,
-      num_ngbs,
+      num_ngbs, /*can convert to comoving=*/1,
       "Integer number of gas neighbour particles within the black hole "
       "kernels.");
 
-  list[31] = io_make_output_field(
+  list[31] = io_make_physical_output_field(
       "NumberOfHeatingEvents", INT, 1, UNIT_CONV_NO_UNITS, 0.f, bparts,
-      AGN_number_of_energy_injections,
+      AGN_number_of_energy_injections, /*can convert to comoving=*/1,
       "Integer number of (thermal) energy injections the black hole has had "
       "so far. This counts each heated gas particle separately, and so can "
       "increase by more than one during a single time step.");
 
-  list[32] = io_make_output_field(
+  list[32] = io_make_physical_output_field(
       "NumberOfAGNEvents", INT, 1, UNIT_CONV_NO_UNITS, 0.f, bparts,
-      AGN_number_of_AGN_events,
+      AGN_number_of_AGN_events, /*can convert to comoving=*/1,
       "Integer number of AGN events the black hole has had so far"
       " (the number of time steps the BH did AGN feedback)");
 
   if (with_cosmology) {
-    list[33] = io_make_output_field(
+    list[33] = io_make_physical_output_field(
         "LastAGNFeedbackScaleFactors", FLOAT, 1, UNIT_CONV_NO_UNITS, 0.f,
-        bparts, last_AGN_event_scale_factor,
+        bparts, last_AGN_event_scale_factor, /*can convert to comoving=*/0,
         "Scale-factors at which the black holes last had an AGN event.");
   } else {
     list[33] = io_make_output_field(
@@ -463,22 +465,22 @@ INLINE static void black_holes_write_particles(const struct bpart* bparts,
       "JetTimeSteps", FLOAT, 1, UNIT_CONV_TIME, 0.f, bparts, dt_jet,
       "Jet-launching-limited time-steps of black holes.");
 
-  list[46] = io_make_output_field(
+  list[46] = io_make_physical_output_field(
       "NumberOfJetParticlesLaunched", INT, 1, UNIT_CONV_NO_UNITS, 0.f, bparts,
-      AGN_number_of_jet_injections,
+      AGN_number_of_jet_injections, /*can convert to comoving=*/1,
       "Integer number of (kinetic) energy injections the black hole has had "
       "so far");
 
-  list[47] = io_make_output_field(
+  list[47] = io_make_physical_output_field(
       "NumberOfAGNJetEvents", INT, 1, UNIT_CONV_NO_UNITS, 0.f, bparts,
-      AGN_number_of_AGN_jet_events,
+      AGN_number_of_AGN_jet_events, /*can convert to comoving=*/1,
       "Integer number of AGN jet launching events the black hole has had"
       " (the number of times the BH did AGN jet feedback)");
 
   if (with_cosmology) {
-    list[48] = io_make_output_field(
+    list[48] = io_make_physical_output_field(
         "LastAGNJetScaleFactors", FLOAT, 1, UNIT_CONV_NO_UNITS, 0.f, bparts,
-        last_AGN_jet_event_scale_factor,
+        last_AGN_jet_event_scale_factor, /*can convert to comoving=*/0,
         "Scale-factors at which the black holes last had an AGN jet event.");
   } else {
     list[48] = io_make_output_field(
diff --git a/src/common_io.c b/src/common_io.c
index 8d3941a246b35d4f490d12af1817bf6bb5ea41ab..8cbdaceb3329cf1b1c7317b42c8183fb5ec862de 100644
--- a/src/common_io.c
+++ b/src/common_io.c
@@ -101,6 +101,8 @@ hid_t io_hdf5_type(enum IO_DATA_TYPE type) {
       return H5T_NATIVE_DOUBLE;
     case CHAR:
       return H5T_NATIVE_CHAR;
+    case BOOL:
+      return H5T_NATIVE_HBOOL;
     default:
       error("Unknown type");
       return 0;
@@ -470,6 +472,16 @@ void io_write_attribute_i(hid_t grp, const char* name, int data) {
   io_write_attribute(grp, name, INT, &data, 1);
 }
 
+/**
+ * @brief Writes a bool value (passed as an int) as an attribute
+ * @param grp The group in which to write
+ * @param name The name of the attribute
+ * @param data The value to write
+ */
+void io_write_attribute_b(hid_t grp, const char* name, int data) {
+  io_write_attribute(grp, name, BOOL, &data, 1);
+}
+
 /**
  * @brief Writes a long value as an attribute
  * @param grp The group in which to write
diff --git a/src/common_io.h b/src/common_io.h
index aced159e04e1bea23c5be84796efda9134896cc7..3e66e574d636ead198e2b82b8eea48051e262bac 100644
--- a/src/common_io.h
+++ b/src/common_io.h
@@ -66,6 +66,7 @@ enum IO_DATA_TYPE {
   FLOAT,
   DOUBLE,
   CHAR,
+  BOOL,
   SIZE_T,
 };
 
@@ -95,6 +96,7 @@ void io_write_attribute(hid_t grp, const char* name, enum IO_DATA_TYPE type,
 void io_write_attribute_d(hid_t grp, const char* name, double data);
 void io_write_attribute_f(hid_t grp, const char* name, float data);
 void io_write_attribute_i(hid_t grp, const char* name, int data);
+void io_write_attribute_b(hid_t grp, const char* name, int data);
 void io_write_attribute_l(hid_t grp, const char* name, long data);
 void io_write_attribute_ll(hid_t grp, const char* name, long long data);
 void io_write_attribute_s(hid_t grp, const char* name, const char* str);
diff --git a/src/common_io_copy.c b/src/common_io_copy.c
index 2be7cb6e4aacfc6796d6592ae1c8f8abfaf36ac2..fdca65cce98fd6c78f06167391100d81b8c878a9 100644
--- a/src/common_io_copy.c
+++ b/src/common_io_copy.c
@@ -488,7 +488,7 @@ void io_copy_temp_buffer(void* temp, const struct engine* e,
 
   } else { /* Converting particle to data */
 
-    if (props.convert_part_f != NULL) {
+    if (props.type == FLOAT && props.parts != NULL) {
 
       /* Prepare some parameters */
       float* temp_f = (float*)temp;
@@ -500,7 +500,7 @@ void io_copy_temp_buffer(void* temp, const struct engine* e,
                      io_convert_part_f_mapper, temp_f, N, copySize,
                      threadpool_auto_chunk_size, (void*)&props);
 
-    } else if (props.convert_part_i != NULL) {
+    } else if (props.type == INT && props.parts != NULL) {
 
       /* Prepare some parameters */
       int* temp_i = (int*)temp;
@@ -512,7 +512,7 @@ void io_copy_temp_buffer(void* temp, const struct engine* e,
                      io_convert_part_i_mapper, temp_i, N, copySize,
                      threadpool_auto_chunk_size, (void*)&props);
 
-    } else if (props.convert_part_d != NULL) {
+    } else if (props.type == DOUBLE && props.parts != NULL) {
 
       /* Prepare some parameters */
       double* temp_d = (double*)temp;
@@ -524,7 +524,7 @@ void io_copy_temp_buffer(void* temp, const struct engine* e,
                      io_convert_part_d_mapper, temp_d, N, copySize,
                      threadpool_auto_chunk_size, (void*)&props);
 
-    } else if (props.convert_part_l != NULL) {
+    } else if (props.type == LONGLONG && props.parts != NULL) {
 
       /* Prepare some parameters */
       long long* temp_l = (long long*)temp;
@@ -536,7 +536,7 @@ void io_copy_temp_buffer(void* temp, const struct engine* e,
                      io_convert_part_l_mapper, temp_l, N, copySize,
                      threadpool_auto_chunk_size, (void*)&props);
 
-    } else if (props.convert_gpart_f != NULL) {
+    } else if (props.type == FLOAT && props.gparts != NULL) {
 
       /* Prepare some parameters */
       float* temp_f = (float*)temp;
@@ -548,7 +548,7 @@ void io_copy_temp_buffer(void* temp, const struct engine* e,
                      io_convert_gpart_f_mapper, temp_f, N, copySize,
                      threadpool_auto_chunk_size, (void*)&props);
 
-    } else if (props.convert_gpart_i != NULL) {
+    } else if (props.type == INT && props.gparts != NULL) {
 
       /* Prepare some parameters */
       int* temp_i = (int*)temp;
@@ -560,7 +560,7 @@ void io_copy_temp_buffer(void* temp, const struct engine* e,
                      io_convert_gpart_i_mapper, temp_i, N, copySize,
                      threadpool_auto_chunk_size, (void*)&props);
 
-    } else if (props.convert_gpart_d != NULL) {
+    } else if (props.type == DOUBLE && props.gparts != NULL) {
 
       /* Prepare some parameters */
       double* temp_d = (double*)temp;
@@ -572,7 +572,7 @@ void io_copy_temp_buffer(void* temp, const struct engine* e,
                      io_convert_gpart_d_mapper, temp_d, N, copySize,
                      threadpool_auto_chunk_size, (void*)&props);
 
-    } else if (props.convert_gpart_l != NULL) {
+    } else if (props.type == LONGLONG && props.gparts != NULL) {
 
       /* Prepare some parameters */
       long long* temp_l = (long long*)temp;
@@ -584,7 +584,7 @@ void io_copy_temp_buffer(void* temp, const struct engine* e,
                      io_convert_gpart_l_mapper, temp_l, N, copySize,
                      threadpool_auto_chunk_size, (void*)&props);
 
-    } else if (props.convert_spart_f != NULL) {
+    } else if (props.type == FLOAT && props.sparts != NULL) {
 
       /* Prepare some parameters */
       float* temp_f = (float*)temp;
@@ -596,7 +596,7 @@ void io_copy_temp_buffer(void* temp, const struct engine* e,
                      io_convert_spart_f_mapper, temp_f, N, copySize,
                      threadpool_auto_chunk_size, (void*)&props);
 
-    } else if (props.convert_spart_i != NULL) {
+    } else if (props.type == INT && props.sparts != NULL) {
 
       /* Prepare some parameters */
       int* temp_i = (int*)temp;
@@ -608,7 +608,7 @@ void io_copy_temp_buffer(void* temp, const struct engine* e,
                      io_convert_spart_i_mapper, temp_i, N, copySize,
                      threadpool_auto_chunk_size, (void*)&props);
 
-    } else if (props.convert_spart_d != NULL) {
+    } else if (props.type == DOUBLE && props.sparts != NULL) {
 
       /* Prepare some parameters */
       double* temp_d = (double*)temp;
@@ -620,7 +620,7 @@ void io_copy_temp_buffer(void* temp, const struct engine* e,
                      io_convert_spart_d_mapper, temp_d, N, copySize,
                      threadpool_auto_chunk_size, (void*)&props);
 
-    } else if (props.convert_spart_l != NULL) {
+    } else if (props.type == LONGLONG && props.sparts != NULL) {
 
       /* Prepare some parameters */
       long long* temp_l = (long long*)temp;
@@ -632,7 +632,7 @@ void io_copy_temp_buffer(void* temp, const struct engine* e,
                      io_convert_spart_l_mapper, temp_l, N, copySize,
                      threadpool_auto_chunk_size, (void*)&props);
 
-    } else if (props.convert_sink_f != NULL) {
+    } else if (props.type == FLOAT && props.sinks != NULL) {
 
       /* Prepare some parameters */
       float* temp_f = (float*)temp;
@@ -644,7 +644,7 @@ void io_copy_temp_buffer(void* temp, const struct engine* e,
                      io_convert_sink_f_mapper, temp_f, N, copySize,
                      threadpool_auto_chunk_size, (void*)&props);
 
-    } else if (props.convert_sink_i != NULL) {
+    } else if (props.type == INT && props.sinks != NULL) {
 
       /* Prepare some parameters */
       int* temp_i = (int*)temp;
@@ -656,7 +656,7 @@ void io_copy_temp_buffer(void* temp, const struct engine* e,
                      io_convert_sink_i_mapper, temp_i, N, copySize,
                      threadpool_auto_chunk_size, (void*)&props);
 
-    } else if (props.convert_sink_d != NULL) {
+    } else if (props.type == DOUBLE && props.sinks != NULL) {
 
       /* Prepare some parameters */
       double* temp_d = (double*)temp;
@@ -668,7 +668,7 @@ void io_copy_temp_buffer(void* temp, const struct engine* e,
                      io_convert_sink_d_mapper, temp_d, N, copySize,
                      threadpool_auto_chunk_size, (void*)&props);
 
-    } else if (props.convert_sink_l != NULL) {
+    } else if (props.type == LONGLONG && props.sinks != NULL) {
 
       /* Prepare some parameters */
       long long* temp_l = (long long*)temp;
@@ -680,7 +680,7 @@ void io_copy_temp_buffer(void* temp, const struct engine* e,
                      io_convert_sink_l_mapper, temp_l, N, copySize,
                      threadpool_auto_chunk_size, (void*)&props);
 
-    } else if (props.convert_bpart_f != NULL) {
+    } else if (props.type == FLOAT && props.bparts != NULL) {
 
       /* Prepare some parameters */
       float* temp_f = (float*)temp;
@@ -692,7 +692,7 @@ void io_copy_temp_buffer(void* temp, const struct engine* e,
                      io_convert_bpart_f_mapper, temp_f, N, copySize,
                      threadpool_auto_chunk_size, (void*)&props);
 
-    } else if (props.convert_bpart_i != NULL) {
+    } else if (props.type == INT && props.bparts != NULL) {
 
       /* Prepare some parameters */
       int* temp_i = (int*)temp;
@@ -704,7 +704,7 @@ void io_copy_temp_buffer(void* temp, const struct engine* e,
                      io_convert_bpart_i_mapper, temp_i, N, copySize,
                      threadpool_auto_chunk_size, (void*)&props);
 
-    } else if (props.convert_bpart_d != NULL) {
+    } else if (props.type == DOUBLE && props.bparts != NULL) {
 
       /* Prepare some parameters */
       double* temp_d = (double*)temp;
@@ -716,7 +716,7 @@ void io_copy_temp_buffer(void* temp, const struct engine* e,
                      io_convert_bpart_d_mapper, temp_d, N, copySize,
                      threadpool_auto_chunk_size, (void*)&props);
 
-    } else if (props.convert_bpart_l != NULL) {
+    } else if (props.type == LONGLONG && props.bparts != NULL) {
 
       /* Prepare some parameters */
       long long* temp_l = (long long*)temp;
@@ -729,7 +729,9 @@ void io_copy_temp_buffer(void* temp, const struct engine* e,
                      threadpool_auto_chunk_size, (void*)&props);
 
     } else {
-      error("Missing conversion function");
+      
+      if (N != 0 && props.ptr_func != NULL)
+	error("Missing conversion function");
     }
   }
 
diff --git a/src/cooling/PS2020/cooling_io.h b/src/cooling/PS2020/cooling_io.h
index dc088df626936a11f879ba2b3cb554d94f70dd8a..242ad0391c95903c055394086bd5efd7b5f116f5 100644
--- a/src/cooling/PS2020/cooling_io.h
+++ b/src/cooling/PS2020/cooling_io.h
@@ -177,18 +177,18 @@ __attribute__((always_inline)) INLINE static int cooling_write_particles(
       "Temperatures", FLOAT, 1, UNIT_CONV_TEMPERATURE, 0.f, parts, xparts,
       convert_part_T, "Temperatures of the gas particles");
 
-  list[1] = io_make_output_field_convert_part(
+  list[1] = io_make_physical_output_field_convert_part(
       "SubgridTemperatures", FLOAT, 1, UNIT_CONV_TEMPERATURE, 0.f, parts,
-      xparts, convert_part_sub_T,
+      xparts, /*can convert to comoving=*/1, convert_part_sub_T,
       "The subgrid temperatures if the particles are within deltaT of the "
       "entropy floor the subgrid temperature is calculated assuming a "
       "pressure equilibrium on the entropy floor, if the particles are "
       "above deltaT of the entropy floor the subgrid temperature is "
       "identical to the SPH temperature.");
 
-  list[2] = io_make_output_field_convert_part(
-      "SubgridPhysicalDensities", FLOAT, 1, UNIT_CONV_DENSITY, 0.f, parts,
-      xparts, convert_part_sub_rho,
+  list[2] = io_make_physical_output_field_convert_part(
+      "SubgridPhysicalDensities", FLOAT, 1, UNIT_CONV_DENSITY, -3.f, parts,
+      xparts, /*can convert to comoving=*/1, convert_part_sub_rho,
       "The subgrid physical density if the particles are within deltaT of the "
       "entropy floor the subgrid density is calculated assuming a pressure "
       "equilibrium on the entropy floor, if the particles are above deltaT "
@@ -222,15 +222,15 @@ __attribute__((always_inline)) INLINE static int cooling_write_particles(
       "floor, by extrapolating to the equilibrium curve assuming constant "
       "pressure.");
 
-  list[6] = io_make_output_field_convert_part(
-      "ElectronNumberDensities", DOUBLE, 1, UNIT_CONV_NUMBER_DENSITY, 0.f,
-      parts, xparts, convert_part_e_density,
+  list[6] = io_make_physical_output_field_convert_part(
+      "ElectronNumberDensities", DOUBLE, 1, UNIT_CONV_NUMBER_DENSITY, -3.f,
+      parts, xparts, /*can convert to comoving=*/1, convert_part_e_density,
       "Electron number densities in the physical frame computed based on the "
       "cooling tables. This is 0 for star-forming particles.");
 
-  list[7] = io_make_output_field_convert_part(
+  list[7] = io_make_physical_output_field_convert_part(
       "ComptonYParameters", DOUBLE, 1, UNIT_CONV_AREA, 0.f, parts, xparts,
-      convert_part_y_compton,
+      /*can convert to comoving=*/0, convert_part_y_compton,
       "Compton y parameters in the physical frame computed based on the "
       "cooling tables. This is 0 for star-forming particles.");
 
diff --git a/src/cooling/const_lambda/cooling_io.h b/src/cooling/const_lambda/cooling_io.h
index 9494731d047b4f1b7b03d03e85d128c286351ae9..120528f3b3087307b10400ec8b8afc4572c1bafb 100644
--- a/src/cooling/const_lambda/cooling_io.h
+++ b/src/cooling/const_lambda/cooling_io.h
@@ -80,9 +80,9 @@ __attribute__((always_inline)) INLINE static int cooling_write_particles(
       "Temperatures", FLOAT, 1, UNIT_CONV_TEMPERATURE, 0.f, parts, xparts,
       convert_part_T, "Temperatures of the gas particles");
 
-  list[1] = io_make_output_field(
+  list[1] = io_make_physical_output_field(
       "RadiatedEnergies", FLOAT, 1, UNIT_CONV_ENERGY, 0.f, xparts,
-      cooling_data.radiated_energy,
+      cooling_data.radiated_energy, /*can convert to comoving=*/0,
       "Thermal energies radiated by the cooling mechanism");
 
   return 2;
diff --git a/src/distributed_io.c b/src/distributed_io.c
index b531a3295712b1d1543ecf2bcc63e1bd1b5b78ef..72bf63302539a2ae3974e9db8660c0247ecbb16e 100644
--- a/src/distributed_io.c
+++ b/src/distributed_io.c
@@ -239,6 +239,9 @@ void write_distributed_array(
   io_write_attribute_f(h_data, "a-scale exponent", props.scale_factor_exponent);
   io_write_attribute_s(h_data, "Expression for physical CGS units", buffer);
   io_write_attribute_s(h_data, "Lossy compression filter", comp_buffer);
+  io_write_attribute_b(h_data, "Value stored as physical", props.is_physical);
+  io_write_attribute_b(h_data, "Property can be converted to comoving",
+                       props.is_convertible_to_comoving);
 
   /* Write the actual number this conversion factor corresponds to */
   const double factor =
@@ -404,6 +407,9 @@ void write_array_virtual(struct engine* e, hid_t grp, const char* fileName_base,
   io_write_attribute_f(h_data, "a-scale exponent", props.scale_factor_exponent);
   io_write_attribute_s(h_data, "Expression for physical CGS units", buffer);
   io_write_attribute_s(h_data, "Lossy compression filter", comp_buffer);
+  io_write_attribute_b(h_data, "Value stored as physical", props.is_physical);
+  io_write_attribute_b(h_data, "Property can be converted to comoving",
+                       props.is_convertible_to_comoving);
 
   /* Write the actual number this conversion factor corresponds to */
   const double factor =
diff --git a/src/extra_io/EAGLE/extra_io.h b/src/extra_io/EAGLE/extra_io.h
index 538e0936068aea9cee49f6646d3209ada15c0482..bdcd2e4b9820d3fedaa57453f6fcf001f388e293 100644
--- a/src/extra_io/EAGLE/extra_io.h
+++ b/src/extra_io/EAGLE/extra_io.h
@@ -63,15 +63,15 @@ INLINE static int extra_io_write_particles(const struct part *parts,
                                            struct io_props *list,
                                            const int with_cosmology) {
 
-  list[0] = io_make_output_field_convert_part(
+  list[0] = io_make_physical_output_field_convert_part(
       "XrayPhotonLuminosities", DOUBLE, 3, UNIT_CONV_PHOTONS_PER_TIME, 0.f,
-      parts, xparts, convert_part_Xray_photons,
+      parts, xparts, /*can convert to comoving=*/0, convert_part_Xray_photons,
       "Intrinsic X-ray photon luminosities in various bands. This is 0 for "
       "star-forming particles.");
 
-  list[1] = io_make_output_field_convert_part(
+  list[1] = io_make_physical_output_field_convert_part(
       "XrayLuminosities", DOUBLE, 3, UNIT_CONV_POWER, 0.f, parts, xparts,
-      convert_part_Xray_energies,
+      /*can convert to comoving=*/0, convert_part_Xray_energies,
       "Intrinsic X-ray luminosities in various bands. This is 0 for "
       "star-forming particles.");
 
diff --git a/src/fof_catalogue_io.c b/src/fof_catalogue_io.c
index bc9060097efa45e890693d73f4d5739d3aaed258..4b44724267f00f642358040be957bc36a2f56c36 100644
--- a/src/fof_catalogue_io.c
+++ b/src/fof_catalogue_io.c
@@ -258,6 +258,9 @@ void write_virtual_fof_hdf5_array(
   io_write_attribute_f(h_data, "a-scale exponent", props.scale_factor_exponent);
   io_write_attribute_s(h_data, "Expression for physical CGS units", buffer);
   io_write_attribute_s(h_data, "Lossy compression filter", comp_buffer);
+  io_write_attribute_b(h_data, "Value stored as physical", props.is_physical);
+  io_write_attribute_b(h_data, "Property can be converted to comoving",
+                       props.is_convertible_to_comoving);
 
   /* Write the actual number this conversion factor corresponds to */
   const double factor =
@@ -316,7 +319,8 @@ void write_fof_virtual_file(const struct fof_props* props,
   struct io_props output_prop;
   output_prop = io_make_output_field_("Masses", DOUBLE, 1, UNIT_CONV_MASS, 0.f,
                                       (char*)props->group_mass, sizeof(double),
-                                      "FOF group masses");
+                                      "FOF group masses", /*physical=*/0,
+                                      /*convertible_to_comoving=*/1);
   write_virtual_fof_hdf5_array(e, h_grp, file_name_base, "Groups", output_prop,
                                num_groups_total, N_counts,
                                compression_write_lossless, e->internal_units,
@@ -324,14 +328,17 @@ void write_fof_virtual_file(const struct fof_props* props,
   output_prop =
       io_make_output_field_("Centres", DOUBLE, 3, UNIT_CONV_LENGTH, 1.f,
                             (char*)props->group_centre_of_mass,
-                            3 * sizeof(double), "FOF group centres of mass");
+                            3 * sizeof(double), "FOF group centres of mass",
+                            /*physical=*/0, /*convertible_to_comoving=*/1);
   write_virtual_fof_hdf5_array(e, h_grp, file_name_base, "Groups", output_prop,
                                num_groups_total, N_counts,
                                compression_write_lossless, e->internal_units,
                                e->snapshot_units);
-  output_prop = io_make_output_field_(
-      "GroupIDs", LONGLONG, 1, UNIT_CONV_NO_UNITS, 0.f,
-      (char*)props->group_index, sizeof(size_t), "FOF group IDs");
+  output_prop =
+      io_make_output_field_("GroupIDs", LONGLONG, 1, UNIT_CONV_NO_UNITS, 0.f,
+                            (char*)props->group_index, sizeof(size_t),
+                            "FOF group IDs", /*physical=*/1,
+                            /*convertible_to_comoving=*/0);
   write_virtual_fof_hdf5_array(e, h_grp, file_name_base, "Groups", output_prop,
                                num_groups_total, N_counts,
                                compression_write_lossless, e->internal_units,
@@ -339,7 +346,8 @@ void write_fof_virtual_file(const struct fof_props* props,
   output_prop =
       io_make_output_field_("Sizes", LONGLONG, 1, UNIT_CONV_NO_UNITS, 0.f,
                             (char*)props->final_group_size, sizeof(long long),
-                            "FOF group length (number of particles)");
+                            "FOF group length (number of particles)",
+                            /*physical=*/1, /*convertible_to_comoving=*/0);
   write_virtual_fof_hdf5_array(e, h_grp, file_name_base, "Groups", output_prop,
                                num_groups_total, N_counts,
                                compression_write_lossless, e->internal_units,
@@ -476,6 +484,9 @@ void write_fof_hdf5_array(
   io_write_attribute_f(h_data, "a-scale exponent", props.scale_factor_exponent);
   io_write_attribute_s(h_data, "Expression for physical CGS units", buffer);
   io_write_attribute_s(h_data, "Lossy compression filter", comp_buffer);
+  io_write_attribute_b(h_data, "Value stored as physical", props.is_physical);
+  io_write_attribute_b(h_data, "Property can be converted to comoving",
+                       props.is_convertible_to_comoving);
 
   /* Write the actual number this conversion factor corresponds to */
   const double factor =
@@ -549,27 +560,32 @@ void write_fof_hdf5_catalogue(const struct fof_props* props,
   struct io_props output_prop;
   output_prop = io_make_output_field_("Masses", DOUBLE, 1, UNIT_CONV_MASS, 0.f,
                                       (char*)props->group_mass, sizeof(double),
-                                      "FOF group masses");
+                                      "FOF group masses", /*physical=*/0,
+                                      /*convertible_to_comoving=*/1);
   write_fof_hdf5_array(e, h_grp, file_name, "Groups", output_prop,
                        num_groups_local, compression_write_lossless,
                        e->internal_units, e->snapshot_units);
   output_prop =
       io_make_output_field_("Centres", DOUBLE, 3, UNIT_CONV_LENGTH, 1.f,
                             (char*)props->group_centre_of_mass,
-                            3 * sizeof(double), "FOF group centres of mass");
+                            3 * sizeof(double), "FOF group centres of mass",
+                            /*physical=*/0, /*convertible_to_comoving=*/1);
   write_fof_hdf5_array(e, h_grp, file_name, "Groups", output_prop,
                        num_groups_local, compression_write_lossless,
                        e->internal_units, e->snapshot_units);
-  output_prop = io_make_output_field_(
-      "GroupIDs", LONGLONG, 1, UNIT_CONV_NO_UNITS, 0.f,
-      (char*)props->group_index, sizeof(size_t), "FOF group IDs");
+  output_prop =
+      io_make_output_field_("GroupIDs", LONGLONG, 1, UNIT_CONV_NO_UNITS, 0.f,
+                            (char*)props->group_index, sizeof(size_t),
+                            "FOF group IDs", /*physical=*/1,
+                            /*convertible_to_comoving=*/0);
   write_fof_hdf5_array(e, h_grp, file_name, "Groups", output_prop,
                        num_groups_local, compression_write_lossless,
                        e->internal_units, e->snapshot_units);
   output_prop =
       io_make_output_field_("Sizes", LONGLONG, 1, UNIT_CONV_NO_UNITS, 0.f,
                             (char*)props->final_group_size, sizeof(long long),
-                            "FOF group length (number of particles)");
+                            "FOF group length (number of particles)",
+                            /*physical=*/1, /*convertible_to_comoving=*/0);
   write_fof_hdf5_array(e, h_grp, file_name, "Groups", output_prop,
                        num_groups_local, compression_write_lossless,
                        e->internal_units, e->snapshot_units);
diff --git a/src/gravity/Default/gravity_io.h b/src/gravity/Default/gravity_io.h
index 9cd4f1e051f9e24698949b98292f6af467626edd..cb1887fe120ca99752e52430db8d9c8af24bf4e9 100644
--- a/src/gravity/Default/gravity_io.h
+++ b/src/gravity/Default/gravity_io.h
@@ -135,9 +135,10 @@ INLINE static void darkmatter_write_particles(const struct gpart* gparts,
   list[2] = io_make_output_field("Masses", FLOAT, 1, UNIT_CONV_MASS, 0.f,
                                  gparts, mass, "Masses of the particles");
 
-  list[3] = io_make_output_field(
+  list[3] = io_make_physical_output_field(
       "ParticleIDs", ULONGLONG, 1, UNIT_CONV_NO_UNITS, 0.f, gparts,
-      id_or_neg_offset, "Unique ID of the particles");
+      id_or_neg_offset, /*can convert to comoving=*/0,
+      "Unique ID of the particles");
 
   list[4] = io_make_output_field_convert_gpart(
       "Potentials", FLOAT, 1, UNIT_CONV_POTENTIAL, -1.f, gparts,
diff --git a/src/gravity/MultiSoftening/gravity_io.h b/src/gravity/MultiSoftening/gravity_io.h
index 81db9762dbb207b0cfbafc83bbabbbc8dc8dbafe..459c15dd34c9bf26ac4f39a4b5d988e34b3a6b15 100644
--- a/src/gravity/MultiSoftening/gravity_io.h
+++ b/src/gravity/MultiSoftening/gravity_io.h
@@ -142,9 +142,10 @@ INLINE static void darkmatter_write_particles(const struct gpart* gparts,
   list[2] = io_make_output_field("Masses", FLOAT, 1, UNIT_CONV_MASS, 0.f,
                                  gparts, mass, "Masses of the particles");
 
-  list[3] = io_make_output_field(
+  list[3] = io_make_physical_output_field(
       "ParticleIDs", ULONGLONG, 1, UNIT_CONV_NO_UNITS, 0.f, gparts,
-      id_or_neg_offset, "Unique ID of the particles");
+      id_or_neg_offset, /*can convert to comoving=*/0,
+      "Unique ID of the particles");
 
   list[4] = io_make_output_field_convert_gpart(
       "Softenings", FLOAT, 1, UNIT_CONV_LENGTH, 1.f, gparts, convert_gpart_soft,
diff --git a/src/hydro/AnarchyPU/hydro_io.h b/src/hydro/AnarchyPU/hydro_io.h
index 3f363780fbe42439b020eb3aa79c448bb76413bb..ded723c9f1fc27613f5f99cfe14cfd203c15ece3 100644
--- a/src/hydro/AnarchyPU/hydro_io.h
+++ b/src/hydro/AnarchyPU/hydro_io.h
@@ -209,9 +209,9 @@ INLINE static void hydro_write_particles(const struct part* parts,
       -3.f * 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 IDs of the particles");
+  list[5] = io_make_physical_output_field(
+      "ParticleIDs", ULONGLONG, 1, UNIT_CONV_NO_UNITS, 0.f, parts, id,
+      /*can convert to comoving=*/0, "Unique IDs of the particles");
 
   list[6] = io_make_output_field("Densities", FLOAT, 1, UNIT_CONV_DENSITY, -3.f,
                                  parts, rho,
diff --git a/src/hydro/Gadget2/hydro_io.h b/src/hydro/Gadget2/hydro_io.h
index 79d86af264dd2172a09451f7795a5f56b1e1066e..2ae73ff7be3db9c0f324ebb1b7b501342a076bd4 100644
--- a/src/hydro/Gadget2/hydro_io.h
+++ b/src/hydro/Gadget2/hydro_io.h
@@ -189,9 +189,9 @@ INLINE static void hydro_write_particles(const struct part* parts,
       "Entropies", FLOAT, 1, UNIT_CONV_ENTROPY_PER_UNIT_MASS, 0.f, parts,
       entropy, "Co-moving entropies per unit mass of the particles");
 
-  list[5] =
-      io_make_output_field("ParticleIDs", ULONGLONG, 1, UNIT_CONV_NO_UNITS, 0.f,
-                           parts, id, "Unique IDs of the particles");
+  list[5] = io_make_physical_output_field(
+      "ParticleIDs", ULONGLONG, 1, UNIT_CONV_NO_UNITS, 0.f, parts, id,
+      /*can convert to comoving=*/0, "Unique IDs of the particles");
 
   list[6] = io_make_output_field("Densities", FLOAT, 1, UNIT_CONV_DENSITY, -3.f,
                                  parts, rho,
diff --git a/src/hydro/Gasoline/hydro_io.h b/src/hydro/Gasoline/hydro_io.h
index b1c69f9af06388dc7dfae8b8f24cc2bcc6ce62de..f066ab5ddd212d9e97677fc44242529b50c1e835 100644
--- a/src/hydro/Gasoline/hydro_io.h
+++ b/src/hydro/Gasoline/hydro_io.h
@@ -204,9 +204,9 @@ INLINE static void hydro_write_particles(const struct part* parts,
       -3.f * 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 IDs of the particles");
+  list[5] = io_make_physical_output_field(
+      "ParticleIDs", ULONGLONG, 1, UNIT_CONV_NO_UNITS, 0.f, parts, id,
+      /*can convert to comoving=*/0, "Unique IDs of the particles");
 
   list[6] = io_make_output_field("Densities", FLOAT, 1, UNIT_CONV_DENSITY, -3.f,
                                  parts, rho,
diff --git a/src/hydro/Gizmo/hydro_io.h b/src/hydro/Gizmo/hydro_io.h
index 25a3d9881dd3d077e67d3d56f5f381a3bd983a8b..1c7a48fa0e79a65e377af3db7a03a0aa81b6b878 100644
--- a/src/hydro/Gizmo/hydro_io.h
+++ b/src/hydro/Gizmo/hydro_io.h
@@ -219,9 +219,9 @@ INLINE static void hydro_write_particles(const struct part* parts,
       -3.f * hydro_gamma_minus_one, parts, xparts, convert_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 IDs of the particles");
+  list[5] = io_make_physical_output_field(
+      "ParticleIDs", ULONGLONG, 1, UNIT_CONV_NO_UNITS, 0.f, parts, id,
+      /*can convert to comoving=*/0, "Unique IDs of the particles");
 
   list[6] = io_make_output_field("Densities", FLOAT, 1, UNIT_CONV_DENSITY, -3.f,
                                  parts, rho,
diff --git a/src/hydro/Minimal/hydro_io.h b/src/hydro/Minimal/hydro_io.h
index 8c32aa7f020fd059fe556e12c147fa397b2d5163..e214a54eb055c9259c6556be8a727e25348c6b9b 100644
--- a/src/hydro/Minimal/hydro_io.h
+++ b/src/hydro/Minimal/hydro_io.h
@@ -204,9 +204,9 @@ INLINE static void hydro_write_particles(const struct part* parts,
       -3.f * 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 IDs of the particles");
+  list[5] = io_make_physical_output_field(
+      "ParticleIDs", ULONGLONG, 1, UNIT_CONV_NO_UNITS, 0.f, parts, id,
+      /*can convert to comoving=*/0, "Unique IDs of the particles");
 
   list[6] = io_make_output_field("Densities", FLOAT, 1, UNIT_CONV_DENSITY, -3.f,
                                  parts, rho,
diff --git a/src/hydro/Phantom/hydro_io.h b/src/hydro/Phantom/hydro_io.h
index 7dbdef3424f02c5a4d3d09d2bcde679269b01455..79a72a835875b30d10850c253895958269413777 100644
--- a/src/hydro/Phantom/hydro_io.h
+++ b/src/hydro/Phantom/hydro_io.h
@@ -210,9 +210,9 @@ INLINE static void hydro_write_particles(const struct part* parts,
       -3.f * 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 IDs of the particles");
+  list[5] = io_make_physical_output_field(
+      "ParticleIDs", ULONGLONG, 1, UNIT_CONV_NO_UNITS, 0.f, parts, id,
+      /*can convert to comoving=*/0, "Unique IDs of the particles");
 
   list[6] = io_make_output_field("Densities", FLOAT, 1, UNIT_CONV_DENSITY, -3.f,
                                  parts, rho,
diff --git a/src/hydro/Planetary/hydro_io.h b/src/hydro/Planetary/hydro_io.h
index f18e1180066513862470e46f99efd5308f0c4495..448e8e0d60ac9065541598d8152bc91f4ece17ab 100644
--- a/src/hydro/Planetary/hydro_io.h
+++ b/src/hydro/Planetary/hydro_io.h
@@ -208,9 +208,9 @@ INLINE static void hydro_write_particles(const struct part* parts,
       "InternalEnergies", FLOAT, 1, UNIT_CONV_ENERGY_PER_UNIT_MASS,
       -3.f * hydro_gamma_minus_one, parts, u,
       "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 IDs of the particles");
+  list[5] = io_make_physical_output_field(
+      "ParticleIDs", ULONGLONG, 1, UNIT_CONV_NO_UNITS, 0.f, parts, id,
+      /*can convert to comoving=*/0, "Unique IDs of the particles");
   list[6] = io_make_output_field("Densities", FLOAT, 1, UNIT_CONV_DENSITY, -3.f,
                                  parts, rho, "Densities of the particles");
   list[7] = io_make_output_field_convert_part(
diff --git a/src/hydro/PressureEnergy/hydro_io.h b/src/hydro/PressureEnergy/hydro_io.h
index 00cc3e4f0c69536a27db5b68374f444ad11098bb..ca8205f70cf928f6e0aa7fb9f3ccb3e9e4272ac0 100644
--- a/src/hydro/PressureEnergy/hydro_io.h
+++ b/src/hydro/PressureEnergy/hydro_io.h
@@ -199,9 +199,9 @@ INLINE static void hydro_write_particles(const struct part* parts,
       -3.f * 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 IDs of the particles");
+  list[5] = io_make_physical_output_field(
+      "ParticleIDs", ULONGLONG, 1, UNIT_CONV_NO_UNITS, 0.f, parts, id,
+      /*can convert to comoving=*/0, "Unique IDs of the particles");
 
   list[6] = io_make_output_field("Densities", FLOAT, 1, UNIT_CONV_DENSITY, -3.f,
                                  parts, rho,
diff --git a/src/hydro/PressureEnergyMorrisMonaghanAV/hydro_io.h b/src/hydro/PressureEnergyMorrisMonaghanAV/hydro_io.h
index 601a99b2d295506f5a7003776add7e2fe33d7d69..467843eeca4d27a5738b1a44bae0cda6e8ddaa27 100644
--- a/src/hydro/PressureEnergyMorrisMonaghanAV/hydro_io.h
+++ b/src/hydro/PressureEnergyMorrisMonaghanAV/hydro_io.h
@@ -200,9 +200,9 @@ INLINE static void hydro_write_particles(const struct part* parts,
       -3.f * 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 IDs of the particles");
+  list[5] = io_make_physical_output_field(
+      "ParticleIDs", ULONGLONG, 1, UNIT_CONV_NO_UNITS, 0.f, parts, id,
+      /*can convert to comoving=*/0, "Unique IDs of the particles");
 
   list[6] = io_make_output_field("Densities", FLOAT, 1, UNIT_CONV_DENSITY, -3.f,
                                  parts, rho,
diff --git a/src/hydro/PressureEntropy/hydro_io.h b/src/hydro/PressureEntropy/hydro_io.h
index 99defbc1f1ce6b86b28c9eadf265d771b20de724..59891d965b1cbc9298e922416c8bb97243515307 100644
--- a/src/hydro/PressureEntropy/hydro_io.h
+++ b/src/hydro/PressureEntropy/hydro_io.h
@@ -200,9 +200,9 @@ INLINE static void hydro_write_particles(const struct part* parts,
       "Entropies", FLOAT, 1, UNIT_CONV_ENTROPY_PER_UNIT_MASS, 0.f, parts,
       entropy, "Co-moving entropies per unit mass of the particles");
 
-  list[5] =
-      io_make_output_field("ParticleIDs", ULONGLONG, 1, UNIT_CONV_NO_UNITS, 0.f,
-                           parts, id, "Unique IDs of the particles");
+  list[5] = io_make_physical_output_field(
+      "ParticleIDs", ULONGLONG, 1, UNIT_CONV_NO_UNITS, 0.f, parts, id,
+      /*can convert to comoving=*/0, "Unique IDs of the particles");
 
   list[6] = io_make_output_field("Densities", FLOAT, 1, UNIT_CONV_DENSITY, -3.f,
                                  parts, rho,
diff --git a/src/hydro/SPHENIX/hydro_io.h b/src/hydro/SPHENIX/hydro_io.h
index 7c4f1d6d85e99c3cd4efb65bdaf7f2a00813b29c..d123f6d1f703bdb45020f04e34b3a6c85b55d417 100644
--- a/src/hydro/SPHENIX/hydro_io.h
+++ b/src/hydro/SPHENIX/hydro_io.h
@@ -217,9 +217,9 @@ INLINE static void hydro_write_particles(const struct part* parts,
       -3.f * 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 IDs of the particles");
+  list[5] = io_make_physical_output_field(
+      "ParticleIDs", ULONGLONG, 1, UNIT_CONV_NO_UNITS, 0.f, parts, id,
+      /*can convert to comoving=*/0, "Unique IDs of the particles");
 
   list[6] = io_make_output_field("Densities", FLOAT, 1, UNIT_CONV_DENSITY, -3.f,
                                  parts, rho,
diff --git a/src/hydro/Shadowswift/hydro_io.h b/src/hydro/Shadowswift/hydro_io.h
index baecc47af05dda85e644186fa9892bd25f052083..3a588d2fd848a37a75c01e84106294e71e52552f 100644
--- a/src/hydro/Shadowswift/hydro_io.h
+++ b/src/hydro/Shadowswift/hydro_io.h
@@ -200,9 +200,9 @@ INLINE static void hydro_write_particles(const struct part* parts,
       -3.f * hydro_gamma_minus_one, parts, xparts, convert_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 IDs of the particles");
+  list[5] = io_make_physical_output_field(
+      "ParticleIDs", ULONGLONG, 1, UNIT_CONV_NO_UNITS, 0.f, parts, id,
+      /*can convert to comoving=*/0, "Unique IDs of the particles");
 
   list[6] = io_make_output_field("Accelerations", FLOAT, 3,
                                  UNIT_CONV_ACCELERATION, 1.f, parts, a_hydro,
diff --git a/src/io_properties.h b/src/io_properties.h
index 35ffef8cb7cc6826d6a627ed41e57fe536dd2564..4f9bc49294a0ab7f2fb567e6656bee402683c840 100644
--- a/src/io_properties.h
+++ b/src/io_properties.h
@@ -110,6 +110,13 @@ struct io_props {
   /* Has this entry been filled */
   int is_used;
 
+  /* Is the entry actually written in the physical frame? */
+  int is_physical;
+
+  /* Is the entry not convertible to comoving (only meaningful if the entry is
+   * physical) */
+  int is_convertible_to_comoving;
+
   /* Is it compulsory ? (input only) */
   enum DATA_IMPORTANCE importance;
 
@@ -152,35 +159,40 @@ struct io_props {
   /* Are we converting? */
   int conversion;
 
-  /* Conversion function for part */
-  conversion_func_part_float convert_part_f;
-  conversion_func_part_int convert_part_i;
-  conversion_func_part_double convert_part_d;
-  conversion_func_part_long_long convert_part_l;
-
-  /* Conversion function for gpart */
-  conversion_func_gpart_float convert_gpart_f;
-  conversion_func_gpart_int convert_gpart_i;
-  conversion_func_gpart_double convert_gpart_d;
-  conversion_func_gpart_long_long convert_gpart_l;
-
-  /* Conversion function for spart */
-  conversion_func_spart_float convert_spart_f;
-  conversion_func_spart_int convert_spart_i;
-  conversion_func_spart_double convert_spart_d;
-  conversion_func_spart_long_long convert_spart_l;
-
-  /* Conversion function for bpart */
-  conversion_func_bpart_float convert_bpart_f;
-  conversion_func_bpart_int convert_bpart_i;
-  conversion_func_bpart_double convert_bpart_d;
-  conversion_func_bpart_long_long convert_bpart_l;
-
-  /* Conversion function for sink */
-  conversion_func_sink_float convert_sink_f;
-  conversion_func_sink_int convert_sink_i;
-  conversion_func_sink_double convert_sink_d;
-  conversion_func_sink_long_long convert_sink_l;
+  union {
+
+    void *ptr_func;
+
+    /* Conversion function for part */
+    conversion_func_part_float convert_part_f;
+    conversion_func_part_int convert_part_i;
+    conversion_func_part_double convert_part_d;
+    conversion_func_part_long_long convert_part_l;
+
+    /* Conversion function for gpart */
+    conversion_func_gpart_float convert_gpart_f;
+    conversion_func_gpart_int convert_gpart_i;
+    conversion_func_gpart_double convert_gpart_d;
+    conversion_func_gpart_long_long convert_gpart_l;
+
+    /* Conversion function for spart */
+    conversion_func_spart_float convert_spart_f;
+    conversion_func_spart_int convert_spart_i;
+    conversion_func_spart_double convert_spart_d;
+    conversion_func_spart_long_long convert_spart_l;
+
+    /* Conversion function for bpart */
+    conversion_func_bpart_float convert_bpart_f;
+    conversion_func_bpart_int convert_bpart_i;
+    conversion_func_bpart_double convert_bpart_d;
+    conversion_func_bpart_long_long convert_bpart_l;
+
+    /* Conversion function for sink */
+    conversion_func_sink_float convert_sink_f;
+    conversion_func_sink_int convert_sink_i;
+    conversion_func_sink_double convert_sink_d;
+    conversion_func_sink_long_long convert_sink_l;
+  };
 };
 
 /**
@@ -276,57 +288,27 @@ INLINE static struct io_props io_make_input_field_(
 #define io_make_output_field(name, type, dim, units, a_exponent, part, field, \
                              desc)                                            \
   io_make_output_field_(name, type, dim, units, a_exponent,                   \
-                        (char *)(&(part[0]).field), sizeof(part[0]), desc)
+                        (char *)(&(part[0]).field), sizeof(part[0]), desc,    \
+                        /*physical=*/0, /*convertible_to_physical=*/1);
 
 /**
- * @brief Construct an #io_props from its parameters
- *
- * @param name Name of the field to read
- * @param type The type of the data
- * @param dimension Dataset dimension (1D, 3D, ...)
- * @param units The units of the dataset
- * @param a_exponent Exponent of the scale-factor to convert to physical units.
- * @param field Pointer to the field of the first particle
- * @param partSize The size in byte of the particle
- * @param description Description of the field added to the meta-data.
+ * @brief Constructs an #io_props from its parameters for a comoving quantity
  *
- * Do not call this function directly. Use the macro defined above.
+ * An alias of io_make_output_field().
  */
-INLINE static struct io_props io_make_output_field_(
-    const char *name, enum IO_DATA_TYPE type, int dimension,
-    enum unit_conversion_factor units, float a_exponent, char *field,
-    size_t partSize, const char *description) {
-
-  struct io_props r;
-  bzero(&r, sizeof(struct io_props));
-
-  safe_strcpy(r.name, name, FIELD_BUFFER_SIZE);
-  if (strlen(description) == 0) {
-    sprintf(r.description, "No description given");
-  } else {
-    safe_strcpy(r.description, description, DESCRIPTION_BUFFER_SIZE);
-  }
-  r.type = type;
-  r.is_used = 1;
-  r.dimension = dimension;
-  r.importance = UNUSED;
-  r.units = units;
-  r.scale_factor_exponent = a_exponent;
-  r.field = field;
-  r.partSize = partSize;
-  r.conversion = 0;
-
-  return r;
-}
+#define io_make_comoving_output_field(name, type, dim, units, a_exponent, \
+                                      part, field, desc)                  \
+  io_make_output_field(name, type, dim, units, a_exponent, part, field, desc)
 
 /**
- * @brief Constructs an #io_props (with conversion) from its parameters
+ * @brief Constructs an #io_props from its parameters for a physical quantity
  */
-#define io_make_output_field_convert_part(name, type, dim, units, a_exponent,  \
-                                          part, xpart, convert, desc)          \
-  io_make_output_field_convert_part_##type(name, type, dim, units, a_exponent, \
-                                           sizeof(part[0]), part, xpart,       \
-                                           convert, desc)
+#define io_make_physical_output_field(name, type, dim, units, a_exponent,  \
+                                      part, field, convertible, desc)      \
+  io_make_output_field_(name, type, dim, units, a_exponent,                \
+                        (char *)(&(part[0]).field), sizeof(part[0]), desc, \
+                        /*physical=*/1,                                    \
+                        /*convertible_to_physical=*/convertible);
 
 /**
  * @brief Construct an #io_props from its parameters
@@ -336,65 +318,22 @@ INLINE static struct io_props io_make_output_field_(
  * @param dimension Dataset dimension (1D, 3D, ...)
  * @param units The units of the dataset
  * @param a_exponent Exponent of the scale-factor to convert to physical units.
+ * @param field Pointer to the field of the first particle
  * @param partSize The size in byte of the particle
- * @param parts The particle array
- * @param xparts The xparticle array
- * @param functionPtr The function used to convert a particle to an int
  * @param description Description of the field added to the meta-data.
+ * @param is_physical Is the quantity written in the physical frame?
+ * @param is_convertible_to_comoving Is the quantity convertible to comoving?
  *
- * Do not call this function directly. Use the macro defined above.
- */
-INLINE static struct io_props io_make_output_field_convert_part_INT(
-    const char *name, enum IO_DATA_TYPE type, int dimension,
-    enum unit_conversion_factor units, float a_exponent, size_t partSize,
-    const struct part *parts, const struct xpart *xparts,
-    conversion_func_part_int functionPtr, const char *description) {
-
-  struct io_props r;
-  bzero(&r, sizeof(struct io_props));
-
-  safe_strcpy(r.name, name, FIELD_BUFFER_SIZE);
-  if (strlen(description) == 0) {
-    sprintf(r.description, "No description given");
-  } else {
-    safe_strcpy(r.description, description, DESCRIPTION_BUFFER_SIZE);
-  }
-  r.type = type;
-  r.is_used = 1;
-  r.dimension = dimension;
-  r.importance = UNUSED;
-  r.units = units;
-  r.scale_factor_exponent = a_exponent;
-  r.partSize = partSize;
-  r.parts = parts;
-  r.xparts = xparts;
-  r.conversion = 1;
-  r.convert_part_i = functionPtr;
-
-  return r;
-}
-
-/**
- * @brief Construct an #io_props from its parameters
- *
- * @param name Name of the field to read
- * @param type The type of the data
- * @param dimension Dataset dimension (1D, 3D, ...)
- * @param units The units of the dataset
- * @param a_exponent Exponent of the scale-factor to convert to physical units.
- * @param partSize The size in byte of the particle
- * @param parts The particle array
- * @param xparts The xparticle array
- * @param functionPtr The function used to convert a particle to a float
- * @param description Description of the field added to the meta-data.
+ * The last argument only makes sense if the quantity is written to
+ * in physical.
  *
  * Do not call this function directly. Use the macro defined above.
  */
-INLINE static struct io_props io_make_output_field_convert_part_FLOAT(
+INLINE static struct io_props io_make_output_field_(
     const char *name, enum IO_DATA_TYPE type, int dimension,
-    enum unit_conversion_factor units, float a_exponent, size_t partSize,
-    const struct part *parts, const struct xpart *xparts,
-    conversion_func_part_float functionPtr, const char *description) {
+    enum unit_conversion_factor units, float a_exponent, char *field,
+    size_t partSize, const char *description, const int is_physical,
+    const int is_convertible_to_comoving) {
 
   struct io_props r;
   bzero(&r, sizeof(struct io_props));
@@ -407,64 +346,39 @@ INLINE static struct io_props io_make_output_field_convert_part_FLOAT(
   }
   r.type = type;
   r.is_used = 1;
+  r.is_physical = is_physical;
+  r.is_convertible_to_comoving = is_convertible_to_comoving;
   r.dimension = dimension;
   r.importance = UNUSED;
   r.units = units;
   r.scale_factor_exponent = a_exponent;
+  r.field = field;
   r.partSize = partSize;
-  r.parts = parts;
-  r.xparts = xparts;
-  r.conversion = 1;
-  r.convert_part_f = functionPtr;
+  r.conversion = 0;
 
   return r;
 }
 
 /**
- * @brief Construct an #io_props from its parameters
- *
- * @param name Name of the field to read
- * @param type The type of the data
- * @param dimension Dataset dimension (1D, 3D, ...)
- * @param units The units of the dataset
- * @param a_exponent Exponent of the scale-factor to convert to physical units.
- * @param partSize The size in byte of the particle
- * @param parts The particle array
- * @param xparts The xparticle array
- * @param functionPtr The function used to convert a particle to a double
- * @param description Description of the field added to the meta-data.
- *
- * Do not call this function directly. Use the macro defined above.
+ * @brief Constructs an #io_props (with conversion) from its parameters
  */
-INLINE static struct io_props io_make_output_field_convert_part_DOUBLE(
-    const char *name, enum IO_DATA_TYPE type, int dimension,
-    enum unit_conversion_factor units, float a_exponent, size_t partSize,
-    const struct part *parts, const struct xpart *xparts,
-    conversion_func_part_double functionPtr, const char *description) {
-
-  struct io_props r;
-  bzero(&r, sizeof(struct io_props));
+#define io_make_output_field_convert_part(name, type, dim, units, a_exponent, \
+                                          part, xpart, convert, desc)         \
+  io_make_output_field_convert_part_(                                         \
+      name, type, dim, units, a_exponent, sizeof(part[0]), part, xpart,       \
+      convert, desc, /*physical=*/0, /*convertible_to_physical=*/1);
 
-  safe_strcpy(r.name, name, FIELD_BUFFER_SIZE);
-  if (strlen(description) == 0) {
-    sprintf(r.description, "No description given");
-  } else {
-    safe_strcpy(r.description, description, DESCRIPTION_BUFFER_SIZE);
-  }
-  r.type = type;
-  r.is_used = 1;
-  r.dimension = dimension;
-  r.importance = UNUSED;
-  r.units = units;
-  r.scale_factor_exponent = a_exponent;
-  r.partSize = partSize;
-  r.parts = parts;
-  r.xparts = xparts;
-  r.conversion = 1;
-  r.convert_part_d = functionPtr;
+#define io_make_comoving_output_field_convert_part(                           \
+    name, type, dim, units, a_exponent, part, xpart, convert, desc)           \
+  io_make_output_field_convert_part(name, type, dim, units, a_exponent, part, \
+                                    xpart, convert, desc)
 
-  return r;
-}
+#define io_make_physical_output_field_convert_part(name, type, dim, units,     \
+                                                   a_exponent, part, xpart,    \
+                                                   convertible, convert, desc) \
+  io_make_output_field_convert_part_(name, type, dim, units, a_exponent,       \
+                                     sizeof(part[0]), part, xpart, convert,    \
+                                     desc, /*physical=*/1, convertible);
 
 /**
  * @brief Construct an #io_props from its parameters
@@ -477,16 +391,17 @@ INLINE static struct io_props io_make_output_field_convert_part_DOUBLE(
  * @param partSize The size in byte of the particle
  * @param parts The particle array
  * @param xparts The xparticle array
- * @param functionPtr The function used to convert a particle to a double
+ * @param functionPtr The function used to convert a particle to an int
  * @param description Description of the field added to the meta-data.
  *
  * Do not call this function directly. Use the macro defined above.
  */
-INLINE static struct io_props io_make_output_field_convert_part_LONGLONG(
+INLINE static struct io_props io_make_output_field_convert_part_(
     const char *name, enum IO_DATA_TYPE type, int dimension,
     enum unit_conversion_factor units, float a_exponent, size_t partSize,
-    const struct part *parts, const struct xpart *xparts,
-    conversion_func_part_long_long functionPtr, const char *description) {
+    const struct part *parts, const struct xpart *xparts, void *functionPtr,
+    const char *description, const int is_physical,
+    const int is_convertible_to_comoving) {
 
   struct io_props r;
   bzero(&r, sizeof(struct io_props));
@@ -499,6 +414,8 @@ INLINE static struct io_props io_make_output_field_convert_part_LONGLONG(
   }
   r.type = type;
   r.is_used = 1;
+  r.is_physical = is_physical;
+  r.is_convertible_to_comoving = is_convertible_to_comoving;
   r.dimension = dimension;
   r.importance = UNUSED;
   r.units = units;
@@ -507,7 +424,7 @@ INLINE static struct io_props io_make_output_field_convert_part_LONGLONG(
   r.parts = parts;
   r.xparts = xparts;
   r.conversion = 1;
-  r.convert_part_l = functionPtr;
+  r.ptr_func = functionPtr;
 
   return r;
 }
@@ -517,53 +434,20 @@ INLINE static struct io_props io_make_output_field_convert_part_LONGLONG(
  */
 #define io_make_output_field_convert_gpart(name, type, dim, units, a_exponent, \
                                            gpart, convert, desc)               \
-  io_make_output_field_convert_gpart_##type(name, type, dim, units,            \
-                                            a_exponent, sizeof(gpart[0]),      \
-                                            gpart, convert, desc)
-
-/**
- * @brief Construct an #io_props from its parameters
- *
- * @param name Name of the field to read
- * @param type The type of the data
- * @param dimension Dataset dimension (1D, 3D, ...)
- * @param units The units of the dataset
- * @param a_exponent Exponent of the scale-factor to convert to physical units.
- * @param gpartSize The size in byte of the particle
- * @param gparts The particle array
- * @param functionPtr The function used to convert a g-particle to a float
- * @param description Description of the field added to the meta-data.
- *
- * Do not call this function directly. Use the macro defined above.
- */
-INLINE static struct io_props io_make_output_field_convert_gpart_INT(
-    const char *name, enum IO_DATA_TYPE type, int dimension,
-    enum unit_conversion_factor units, float a_exponent, size_t gpartSize,
-    const struct gpart *gparts, conversion_func_gpart_int functionPtr,
-    const char *description) {
-
-  struct io_props r;
-  bzero(&r, sizeof(struct io_props));
+  io_make_output_field_convert_gpart_(                                         \
+      name, type, dim, units, a_exponent, sizeof(gpart[0]), gpart, convert,    \
+      desc, /*physical=*/0, /*convertible_to_physical=*/1)
 
-  safe_strcpy(r.name, name, FIELD_BUFFER_SIZE);
-  if (strlen(description) == 0) {
-    sprintf(r.description, "No description given");
-  } else {
-    safe_strcpy(r.description, description, DESCRIPTION_BUFFER_SIZE);
-  }
-  r.type = type;
-  r.is_used = 1;
-  r.dimension = dimension;
-  r.importance = UNUSED;
-  r.units = units;
-  r.scale_factor_exponent = a_exponent;
-  r.partSize = gpartSize;
-  r.gparts = gparts;
-  r.conversion = 1;
-  r.convert_gpart_i = functionPtr;
+#define io_make_comoving_output_field_convert_gpart(                     \
+    name, type, dim, units, a_exponent, gpart, convert, desc)            \
+  io_make_output_field_convert_gpart(name, type, dim, units, a_exponent, \
+                                     gpart, convert, desc)
 
-  return r;
-}
+#define io_make_physical_output_field_convert_gpart(                         \
+    name, type, dim, units, a_exponent, gpart, convertible, convert, desc)   \
+  io_make_output_field_convert_gpart_(name, type, dim, units, a_exponent,    \
+                                      sizeof(gpart[0]), gpart, convert, desc, \
+                                      /*physical=*/1, convertible);
 
 /**
  * @brief Construct an #io_props from its parameters
@@ -580,11 +464,11 @@ INLINE static struct io_props io_make_output_field_convert_gpart_INT(
  *
  * Do not call this function directly. Use the macro defined above.
  */
-INLINE static struct io_props io_make_output_field_convert_gpart_FLOAT(
+INLINE static struct io_props io_make_output_field_convert_gpart_(
     const char *name, enum IO_DATA_TYPE type, int dimension,
     enum unit_conversion_factor units, float a_exponent, size_t gpartSize,
-    const struct gpart *gparts, conversion_func_gpart_float functionPtr,
-    const char *description) {
+    const struct gpart *gparts, void *functionPtr, const char *description,
+    const int is_physical, const int is_convertible_to_comoving) {
 
   struct io_props r;
   bzero(&r, sizeof(struct io_props));
@@ -597,6 +481,8 @@ INLINE static struct io_props io_make_output_field_convert_gpart_FLOAT(
   }
   r.type = type;
   r.is_used = 1;
+  r.is_physical = is_physical;
+  r.is_convertible_to_comoving = is_convertible_to_comoving;
   r.dimension = dimension;
   r.importance = UNUSED;
   r.units = units;
@@ -604,54 +490,30 @@ INLINE static struct io_props io_make_output_field_convert_gpart_FLOAT(
   r.partSize = gpartSize;
   r.gparts = gparts;
   r.conversion = 1;
-  r.convert_gpart_f = functionPtr;
+  r.ptr_func = functionPtr;
 
   return r;
 }
 
 /**
- * @brief Construct an #io_props from its parameters
- *
- * @param name Name of the field to read
- * @param type The type of the data
- * @param dimension Dataset dimension (1D, 3D, ...)
- * @param units The units of the dataset
- * @param a_exponent Exponent of the scale-factor to convert to physical units.
- * @param gpartSize The size in byte of the particle
- * @param gparts The particle array
- * @param functionPtr The function used to convert a g-particle to a double
- * @param description Description of the field added to the meta-data.
- *
- * Do not call this function directly. Use the macro defined above.
+ * @brief Constructs an #io_props (with conversion) from its parameters
  */
-INLINE static struct io_props io_make_output_field_convert_gpart_DOUBLE(
-    const char *name, enum IO_DATA_TYPE type, int dimension,
-    enum unit_conversion_factor units, float a_exponent, size_t gpartSize,
-    const struct gpart *gparts, conversion_func_gpart_double functionPtr,
-    const char *description) {
-
-  struct io_props r;
-  bzero(&r, sizeof(struct io_props));
+#define io_make_output_field_convert_spart(name, type, dim, units, a_exponent, \
+                                           spart, convert, desc)               \
+  io_make_output_field_convert_spart_(                                         \
+      name, type, dim, units, a_exponent, sizeof(spart[0]), spart, convert,    \
+      desc, /*physical=*/0, /*convertible_to_physical=*/1)
 
-  safe_strcpy(r.name, name, FIELD_BUFFER_SIZE);
-  if (strlen(description) == 0) {
-    sprintf(r.description, "No description given");
-  } else {
-    safe_strcpy(r.description, description, DESCRIPTION_BUFFER_SIZE);
-  }
-  r.type = type;
-  r.is_used = 1;
-  r.dimension = dimension;
-  r.importance = UNUSED;
-  r.units = units;
-  r.scale_factor_exponent = a_exponent;
-  r.partSize = gpartSize;
-  r.gparts = gparts;
-  r.conversion = 1;
-  r.convert_gpart_d = functionPtr;
+#define io_make_comoving_output_field_convert_spart(                     \
+    name, type, dim, units, a_exponent, spart, convert, desc)            \
+  io_make_output_field_convert_spart(name, type, dim, units, a_exponent, \
+                                     spart, convert, desc)
 
-  return r;
-}
+#define io_make_physical_output_field_convert_spart(                          \
+    name, type, dim, units, a_exponent, spart, convertible, convert, desc)    \
+  io_make_output_field_convert_spart_(name, type, dim, units, a_exponent,     \
+                                      sizeof(spart[0]), spart, convert, desc, \
+                                      /*physical=*/1, convertible);
 
 /**
  * @brief Construct an #io_props from its parameters
@@ -661,18 +523,18 @@ INLINE static struct io_props io_make_output_field_convert_gpart_DOUBLE(
  * @param dimension Dataset dimension (1D, 3D, ...)
  * @param units The units of the dataset
  * @param a_exponent Exponent of the scale-factor to convert to physical units.
- * @param gpartSize The size in byte of the particle
- * @param gparts The particle array
- * @param functionPtr The function used to convert a g-particle to a double
+ * @param spartSize The size in byte of the particle
+ * @param sparts The particle array
+ * @param functionPtr The function used to convert a s-particle to a float
  * @param description Description of the field added to the meta-data.
  *
  * Do not call this function directly. Use the macro defined above.
  */
-INLINE static struct io_props io_make_output_field_convert_gpart_LONGLONG(
+INLINE static struct io_props io_make_output_field_convert_spart_(
     const char *name, enum IO_DATA_TYPE type, int dimension,
-    enum unit_conversion_factor units, float a_exponent, size_t gpartSize,
-    const struct gpart *gparts, conversion_func_gpart_long_long functionPtr,
-    const char *description) {
+    enum unit_conversion_factor units, float a_exponent, size_t spartSize,
+    const struct spart *sparts, void *functionPtr, const char *description,
+    const int is_physical, const int is_convertible_to_comoving) {
 
   struct io_props r;
   bzero(&r, sizeof(struct io_props));
@@ -685,14 +547,16 @@ INLINE static struct io_props io_make_output_field_convert_gpart_LONGLONG(
   }
   r.type = type;
   r.is_used = 1;
+  r.is_physical = is_physical;
+  r.is_convertible_to_comoving = is_convertible_to_comoving;
   r.dimension = dimension;
   r.importance = UNUSED;
   r.units = units;
   r.scale_factor_exponent = a_exponent;
-  r.partSize = gpartSize;
-  r.gparts = gparts;
+  r.partSize = spartSize;
+  r.sparts = sparts;
   r.conversion = 1;
-  r.convert_gpart_l = functionPtr;
+  r.ptr_func = functionPtr;
 
   return r;
 }
@@ -700,11 +564,22 @@ INLINE static struct io_props io_make_output_field_convert_gpart_LONGLONG(
 /**
  * @brief Constructs an #io_props (with conversion) from its parameters
  */
-#define io_make_output_field_convert_spart(name, type, dim, units, a_exponent, \
-                                           spart, convert, desc)               \
-  io_make_output_field_convert_spart_##type(name, type, dim, units,            \
-                                            a_exponent, sizeof(spart[0]),      \
-                                            spart, convert, desc)
+#define io_make_output_field_convert_bpart(name, type, dim, units, a_exponent, \
+                                           bpart, convert, desc)               \
+  io_make_output_field_convert_bpart_(                                         \
+      name, type, dim, units, a_exponent, sizeof(bpart[0]), bpart, convert,    \
+      desc, /*physical=*/0, /*convertible_to_physical=*/1)
+
+#define io_make_comoving_output_field_convert_bpart(                     \
+    name, type, dim, units, a_exponent, bpart, convert, desc)            \
+  io_make_output_field_convert_bpart(name, type, dim, units, a_exponent, \
+                                     bpart, convert, desc)
+
+#define io_make_physical_output_field_convert_bpart(                          \
+    name, type, dim, units, a_exponent, bpart, convertible, convert, desc)    \
+  io_make_output_field_convert_bpart_(name, type, dim, units, a_exponent,     \
+                                      sizeof(bpart[0]), bpart, convert, desc, \
+                                      /*physical=*/1, convertible);
 
 /**
  * @brief Construct an #io_props from its parameters
@@ -714,18 +589,18 @@ INLINE static struct io_props io_make_output_field_convert_gpart_LONGLONG(
  * @param dimension Dataset dimension (1D, 3D, ...)
  * @param units The units of the dataset
  * @param a_exponent Exponent of the scale-factor to convert to physical units.
- * @param spartSize The size in byte of the particle
- * @param sparts The particle array
- * @param functionPtr The function used to convert a s-particle to a float
+ * @param bpartSize The size in byte of the particle
+ * @param bparts The particle array
+ * @param functionPtr The function used to convert a b-particle to a float
  * @param description Description of the field added to the meta-data.
  *
  * Do not call this function directly. Use the macro defined above.
  */
-INLINE static struct io_props io_make_output_field_convert_spart_INT(
+INLINE static struct io_props io_make_output_field_convert_bpart_(
     const char *name, enum IO_DATA_TYPE type, int dimension,
-    enum unit_conversion_factor units, float a_exponent, size_t spartSize,
-    const struct spart *sparts, conversion_func_spart_int functionPtr,
-    const char *description) {
+    enum unit_conversion_factor units, float a_exponent, size_t bpartSize,
+    const struct bpart *bparts, void *functionPtr, const char *description,
+    const int is_physical, const int is_convertible_to_comoving) {
 
   struct io_props r;
   bzero(&r, sizeof(struct io_props));
@@ -738,344 +613,40 @@ INLINE static struct io_props io_make_output_field_convert_spart_INT(
   }
   r.type = type;
   r.is_used = 1;
+  r.is_physical = is_physical;
+  r.is_convertible_to_comoving = is_convertible_to_comoving;
   r.dimension = dimension;
   r.importance = UNUSED;
   r.units = units;
   r.scale_factor_exponent = a_exponent;
-  r.partSize = spartSize;
-  r.sparts = sparts;
+  r.partSize = bpartSize;
+  r.bparts = bparts;
   r.conversion = 1;
-  r.convert_spart_i = functionPtr;
+  r.ptr_func = functionPtr;
 
   return r;
 }
 
-/**
- * @brief Construct an #io_props from its parameters
- *
- * @param name Name of the field to read
- * @param type The type of the data
- * @param dimension Dataset dimension (1D, 3D, ...)
- * @param units The units of the dataset
- * @param a_exponent Exponent of the scale-factor to convert to physical units.
- * @param spartSize The size in byte of the particle
- * @param sparts The particle array
- * @param functionPtr The function used to convert a g-particle to a float
- * @param description Description of the field added to the meta-data.
- *
- * Do not call this function directly. Use the macro defined above.
- */
-INLINE static struct io_props io_make_output_field_convert_spart_FLOAT(
-    const char *name, enum IO_DATA_TYPE type, int dimension,
-    enum unit_conversion_factor units, float a_exponent, size_t spartSize,
-    const struct spart *sparts, conversion_func_spart_float functionPtr,
-    const char *description) {
-
-  struct io_props r;
-  bzero(&r, sizeof(struct io_props));
-
-  safe_strcpy(r.name, name, FIELD_BUFFER_SIZE);
-  if (strlen(description) == 0) {
-    sprintf(r.description, "No description given");
-  } else {
-    safe_strcpy(r.description, description, DESCRIPTION_BUFFER_SIZE);
-  }
-  r.type = type;
-  r.is_used = 1;
-  r.dimension = dimension;
-  r.importance = UNUSED;
-  r.units = units;
-  r.scale_factor_exponent = a_exponent;
-  r.partSize = spartSize;
-  r.sparts = sparts;
-  r.conversion = 1;
-  r.convert_spart_f = functionPtr;
-
-  return r;
-}
-
-/**
- * @brief Construct an #io_props from its parameters
- *
- * @param name Name of the field to read
- * @param type The type of the data
- * @param dimension Dataset dimension (1D, 3D, ...)
- * @param units The units of the dataset
- * @param a_exponent Exponent of the scale-factor to convert to physical units.
- * @param spartSize The size in byte of the particle
- * @param sparts The particle array
- * @param functionPtr The function used to convert a s-particle to a double
- * @param description Description of the field added to the meta-data.
- *
- * Do not call this function directly. Use the macro defined above.
- */
-INLINE static struct io_props io_make_output_field_convert_spart_DOUBLE(
-    const char *name, enum IO_DATA_TYPE type, int dimension,
-    enum unit_conversion_factor units, float a_exponent, size_t spartSize,
-    const struct spart *sparts, conversion_func_spart_double functionPtr,
-    const char *description) {
-
-  struct io_props r;
-  bzero(&r, sizeof(struct io_props));
-
-  safe_strcpy(r.name, name, FIELD_BUFFER_SIZE);
-  if (strlen(description) == 0) {
-    sprintf(r.description, "No description given");
-  } else {
-    safe_strcpy(r.description, description, DESCRIPTION_BUFFER_SIZE);
-  }
-  r.type = type;
-  r.is_used = 1;
-  r.dimension = dimension;
-  r.importance = UNUSED;
-  r.units = units;
-  r.scale_factor_exponent = a_exponent;
-  r.partSize = spartSize;
-  r.sparts = sparts;
-  r.conversion = 1;
-  r.convert_spart_d = functionPtr;
-
-  return r;
-}
-
-/**
- * @brief Construct an #io_props from its parameters
- *
- * @param name Name of the field to read
- * @param type The type of the data
- * @param dimension Dataset dimension (1D, 3D, ...)
- * @param units The units of the dataset
- * @param a_exponent Exponent of the scale-factor to convert to physical units.
- * @param spartSize The size in byte of the particle
- * @param sparts The particle array
- * @param functionPtr The function used to convert a s-particle to a double
- * @param description Description of the field added to the meta-data.
- *
- * Do not call this function directly. Use the macro defined above.
- */
-INLINE static struct io_props io_make_output_field_convert_spart_LONGLONG(
-    const char *name, enum IO_DATA_TYPE type, int dimension,
-    enum unit_conversion_factor units, float a_exponent, size_t spartSize,
-    const struct spart *sparts, conversion_func_spart_long_long functionPtr,
-    const char *description) {
-
-  struct io_props r;
-  bzero(&r, sizeof(struct io_props));
-
-  safe_strcpy(r.name, name, FIELD_BUFFER_SIZE);
-  if (strlen(description) == 0) {
-    sprintf(r.description, "No description given");
-  } else {
-    safe_strcpy(r.description, description, DESCRIPTION_BUFFER_SIZE);
-  }
-  r.type = type;
-  r.is_used = 1;
-  r.dimension = dimension;
-  r.importance = UNUSED;
-  r.units = units;
-  r.scale_factor_exponent = a_exponent;
-  r.partSize = spartSize;
-  r.sparts = sparts;
-  r.conversion = 1;
-  r.convert_spart_l = functionPtr;
-
-  return r;
-}
-
-/**
- * @brief Constructs an #io_props (with conversion) from its parameters
- */
-#define io_make_output_field_convert_bpart(name, type, dim, units, a_exponent, \
-                                           bpart, convert, desc)               \
-  io_make_output_field_convert_bpart_##type(name, type, dim, units,            \
-                                            a_exponent, sizeof(bpart[0]),      \
-                                            bpart, convert, desc)
-
-/**
- * @brief Construct an #io_props from its parameters
- *
- * @param name Name of the field to read
- * @param type The type of the data
- * @param dimension Dataset dimension (1D, 3D, ...)
- * @param units The units of the dataset
- * @param a_exponent Exponent of the scale-factor to convert to physical units.
- * @param bpartSize The size in byte of the particle
- * @param bparts The particle array
- * @param functionPtr The function used to convert a b-particle to a float
- * @param description Description of the field added to the meta-data.
- *
- * Do not call this function directly. Use the macro defined above.
- */
-INLINE static struct io_props io_make_output_field_convert_bpart_INT(
-    const char *name, enum IO_DATA_TYPE type, int dimension,
-    enum unit_conversion_factor units, float a_exponent, size_t bpartSize,
-    const struct bpart *bparts, conversion_func_bpart_int functionPtr,
-    const char *description) {
-
-  struct io_props r;
-  bzero(&r, sizeof(struct io_props));
-
-  safe_strcpy(r.name, name, FIELD_BUFFER_SIZE);
-  if (strlen(description) == 0) {
-    sprintf(r.description, "No description given");
-  } else {
-    safe_strcpy(r.description, description, DESCRIPTION_BUFFER_SIZE);
-  }
-  r.type = type;
-  r.is_used = 1;
-  r.dimension = dimension;
-  r.importance = UNUSED;
-  r.units = units;
-  r.scale_factor_exponent = a_exponent;
-  r.partSize = bpartSize;
-  r.bparts = bparts;
-  r.conversion = 1;
-  r.convert_bpart_i = functionPtr;
-
-  return r;
-}
-
-/**
- * @brief Construct an #io_props from its parameters
- *
- * @param name Name of the field to read
- * @param type The type of the data
- * @param dimension Dataset dimension (1D, 3D, ...)
- * @param units The units of the dataset
- * @param a_exponent Exponent of the scale-factor to convert to physical units.
- * @param bpartSize The size in byte of the particle
- * @param bparts The particle array
- * @param functionPtr The function used to convert a g-particle to a float
- * @param description Description of the field added to the meta-data.
- *
- * Do not call this function directly. Use the macro defined above.
- */
-INLINE static struct io_props io_make_output_field_convert_bpart_FLOAT(
-    const char *name, enum IO_DATA_TYPE type, int dimension,
-    enum unit_conversion_factor units, float a_exponent, size_t bpartSize,
-    const struct bpart *bparts, conversion_func_bpart_float functionPtr,
-    const char *description) {
-
-  struct io_props r;
-  bzero(&r, sizeof(struct io_props));
-
-  safe_strcpy(r.name, name, FIELD_BUFFER_SIZE);
-  if (strlen(description) == 0) {
-    sprintf(r.description, "No description given");
-  } else {
-    safe_strcpy(r.description, description, DESCRIPTION_BUFFER_SIZE);
-  }
-  r.type = type;
-  r.is_used = 1;
-  r.dimension = dimension;
-  r.importance = UNUSED;
-  r.units = units;
-  r.scale_factor_exponent = a_exponent;
-  r.partSize = bpartSize;
-  r.bparts = bparts;
-  r.conversion = 1;
-  r.convert_bpart_f = functionPtr;
-
-  return r;
-}
-
-/**
- * @brief Construct an #io_props from its parameters
- *
- * @param name Name of the field to read
- * @param type The type of the data
- * @param dimension Dataset dimension (1D, 3D, ...)
- * @param units The units of the dataset
- * @param a_exponent Exponent of the scale-factor to convert to physical units.
- * @param bpartSize The size in byte of the particle
- * @param bparts The particle array
- * @param functionPtr The function used to convert a s-particle to a double
- * @param description Description of the field added to the meta-data.
- *
- * Do not call this function directly. Use the macro defined above.
- */
-INLINE static struct io_props io_make_output_field_convert_bpart_DOUBLE(
-    const char *name, enum IO_DATA_TYPE type, int dimension,
-    enum unit_conversion_factor units, float a_exponent, size_t bpartSize,
-    const struct bpart *bparts, conversion_func_bpart_double functionPtr,
-    const char *description) {
-
-  struct io_props r;
-  bzero(&r, sizeof(struct io_props));
-
-  safe_strcpy(r.name, name, FIELD_BUFFER_SIZE);
-  if (strlen(description) == 0) {
-    sprintf(r.description, "No description given");
-  } else {
-    safe_strcpy(r.description, description, DESCRIPTION_BUFFER_SIZE);
-  }
-  r.type = type;
-  r.is_used = 1;
-  r.dimension = dimension;
-  r.importance = UNUSED;
-  r.units = units;
-  r.scale_factor_exponent = a_exponent;
-  r.partSize = bpartSize;
-  r.bparts = bparts;
-  r.conversion = 1;
-  r.convert_bpart_d = functionPtr;
-
-  return r;
-}
-
-/**
- * @brief Construct an #io_props from its parameters
- *
- * @param name Name of the field to read
- * @param type The type of the data
- * @param dimension Dataset dimension (1D, 3D, ...)
- * @param units The units of the dataset
- * @param a_exponent Exponent of the scale-factor to convert to physical units.
- * @param bpartSize The size in byte of the particle
- * @param bparts The particle array
- * @param functionPtr The function used to convert a s-particle to a double
- * @param description Description of the field added to the meta-data.
- *
- * Do not call this function directly. Use the macro defined above.
- */
-INLINE static struct io_props io_make_output_field_convert_bpart_LONGLONG(
-    const char *name, enum IO_DATA_TYPE type, int dimension,
-    enum unit_conversion_factor units, float a_exponent, size_t bpartSize,
-    const struct bpart *bparts, conversion_func_bpart_long_long functionPtr,
-    const char *description) {
-
-  struct io_props r;
-  bzero(&r, sizeof(struct io_props));
-
-  safe_strcpy(r.name, name, FIELD_BUFFER_SIZE);
-  if (strlen(description) == 0) {
-    sprintf(r.description, "No description given");
-  } else {
-    safe_strcpy(r.description, description, DESCRIPTION_BUFFER_SIZE);
-  }
-  r.type = type;
-  r.is_used = 1;
-  r.dimension = dimension;
-  r.importance = UNUSED;
-  r.units = units;
-  r.scale_factor_exponent = a_exponent;
-  r.partSize = bpartSize;
-  r.bparts = bparts;
-  r.conversion = 1;
-  r.convert_bpart_l = functionPtr;
-
-  return r;
-}
-
-/**
- * @brief Constructs an #io_props (with conversion) from its parameters
- */
-#define io_make_output_field_convert_sink(name, type, dim, units, a_exponent,  \
-                                          sink, convert, desc)                 \
-  io_make_output_field_convert_sink_##type(name, type, dim, units, a_exponent, \
-                                           sizeof(sink[0]), sink, convert,     \
-                                           desc)
-
+/**
+ * @brief Constructs an #io_props (with conversion) from its parameters
+ */
+#define io_make_output_field_convert_sink(name, type, dim, units, a_exponent, \
+                                          sink, convert, desc)                \
+  io_make_output_field_convert_sink_(                                         \
+      name, type, dim, units, a_exponent, sizeof(sink[0]), sink, convert,     \
+      desc, /*physical=*/0, /*convertible_to_physical=*/1)
+
+#define io_make_comoving_output_field_convert_sink(                           \
+    name, type, dim, units, a_exponent, ink, convert, desc)                   \
+  io_make_output_field_convert_sink(name, type, dim, units, a_exponent, sink, \
+                                    convert, desc)
+
+#define io_make_physical_output_field_convert_sink(                        \
+    name, type, dim, units, a_exponent, sink, convertible, convert, desc)  \
+  io_make_output_field_convert_sink_(name, type, dim, units, a_exponent,   \
+                                     sizeof(part[0]), sink, convert, desc, \
+                                     /*physical=*/1, convertible);
+
 /**
  * @brief Construct an #io_props from its parameters
  *
@@ -1091,143 +662,11 @@ INLINE static struct io_props io_make_output_field_convert_bpart_LONGLONG(
  *
  * Do not call this function directly. Use the macro defined above.
  */
-INLINE static struct io_props io_make_output_field_convert_sink_INT(
-    const char *name, enum IO_DATA_TYPE type, int dimension,
-    enum unit_conversion_factor units, float a_exponent, size_t sinkSize,
-    const struct sink *sinks, conversion_func_sink_int functionPtr,
-    const char *description) {
-
-  struct io_props r;
-  bzero(&r, sizeof(struct io_props));
-
-  safe_strcpy(r.name, name, FIELD_BUFFER_SIZE);
-  if (strlen(description) == 0) {
-    sprintf(r.description, "No description given");
-  } else {
-    safe_strcpy(r.description, description, DESCRIPTION_BUFFER_SIZE);
-  }
-  r.type = type;
-  r.is_used = 1;
-  r.dimension = dimension;
-  r.importance = UNUSED;
-  r.units = units;
-  r.scale_factor_exponent = a_exponent;
-  r.partSize = sinkSize;
-  r.sinks = sinks;
-  r.conversion = 1;
-  r.convert_sink_i = functionPtr;
-
-  return r;
-}
-
-/**
- * @brief Construct an #io_props from its parameters
- *
- * @param name Name of the field to read
- * @param type The type of the data
- * @param dimension Dataset dimension (1D, 3D, ...)
- * @param units The units of the dataset
- * @param a_exponent Exponent of the scale-factor to convert to physical units.
- * @param sinkSize The size in byte of the particle
- * @param sinks The particle array
- * @param functionPtr The function used to convert a sink-particle to a float
- * @param description Description of the field added to the meta-data.
- *
- * Do not call this function directly. Use the macro defined above.
- */
-INLINE static struct io_props io_make_output_field_convert_sink_FLOAT(
-    const char *name, enum IO_DATA_TYPE type, int dimension,
-    enum unit_conversion_factor units, float a_exponent, size_t sinkSize,
-    const struct sink *sinks, conversion_func_sink_float functionPtr,
-    const char *description) {
-
-  struct io_props r;
-  bzero(&r, sizeof(struct io_props));
-
-  safe_strcpy(r.name, name, FIELD_BUFFER_SIZE);
-  if (strlen(description) == 0) {
-    sprintf(r.description, "No description given");
-  } else {
-    safe_strcpy(r.description, description, DESCRIPTION_BUFFER_SIZE);
-  }
-  r.type = type;
-  r.is_used = 1;
-  r.dimension = dimension;
-  r.importance = UNUSED;
-  r.units = units;
-  r.scale_factor_exponent = a_exponent;
-  r.partSize = sinkSize;
-  r.sinks = sinks;
-  r.conversion = 1;
-  r.convert_sink_f = functionPtr;
-
-  return r;
-}
-
-/**
- * @brief Construct an #io_props from its parameters
- *
- * @param name Name of the field to read
- * @param type The type of the data
- * @param dimension Dataset dimension (1D, 3D, ...)
- * @param units The units of the dataset
- * @param a_exponent Exponent of the scale-factor to convert to physical units.
- * @param sinkSize The size in byte of the particle
- * @param sinks The particle array
- * @param functionPtr The function used to convert a sink-particle to a double
- * @param description Description of the field added to the meta-data.
- *
- * Do not call this function directly. Use the macro defined above.
- */
-INLINE static struct io_props io_make_output_field_convert_sink_DOUBLE(
-    const char *name, enum IO_DATA_TYPE type, int dimension,
-    enum unit_conversion_factor units, float a_exponent, size_t sinkSize,
-    const struct sink *sinks, conversion_func_sink_double functionPtr,
-    const char *description) {
-
-  struct io_props r;
-  bzero(&r, sizeof(struct io_props));
-
-  safe_strcpy(r.name, name, FIELD_BUFFER_SIZE);
-  if (strlen(description) == 0) {
-    sprintf(r.description, "No description given");
-  } else {
-    safe_strcpy(r.description, description, DESCRIPTION_BUFFER_SIZE);
-  }
-  r.type = type;
-  r.is_used = 1;
-  r.dimension = dimension;
-  r.importance = UNUSED;
-  r.units = units;
-  r.scale_factor_exponent = a_exponent;
-  r.partSize = sinkSize;
-  r.sinks = sinks;
-  r.conversion = 1;
-  r.convert_sink_d = functionPtr;
-
-  return r;
-}
-
-/**
- * @brief Construct an #io_props from its parameters
- *
- * @param name Name of the field to read
- * @param type The type of the data
- * @param dimension Dataset dimension (1D, 3D, ...)
- * @param units The units of the dataset
- * @param a_exponent Exponent of the scale-factor to convert to physical units.
- * @param sinkSize The size in byte of the particle
- * @param sinks The particle array
- * @param functionPtr The function used to convert a sink-particle to a double
- * @param description Description of the field added to the meta-data.
- *
- * Do not call this function directly. Use the macro defined above.
- */
-INLINE static struct io_props io_make_output_field_convert_sink_LONGLONG(
+INLINE static struct io_props io_make_output_field_convert_sink_(
     const char *name, enum IO_DATA_TYPE type, int dimension,
     enum unit_conversion_factor units, float a_exponent, size_t sinkSize,
-    const struct sink *sinks, conversion_func_sink_long_long functionPtr,
-    const char *description) {
+    const struct sink *sinks, void *functionPtr, const char *description,
+    const int is_physical, const int is_convertible_to_comoving) {
 
   struct io_props r;
   bzero(&r, sizeof(struct io_props));
@@ -1240,6 +679,8 @@ INLINE static struct io_props io_make_output_field_convert_sink_LONGLONG(
   }
   r.type = type;
   r.is_used = 1;
+  r.is_physical = is_physical;
+  r.is_convertible_to_comoving = is_convertible_to_comoving;
   r.dimension = dimension;
   r.importance = UNUSED;
   r.units = units;
@@ -1247,7 +688,7 @@ INLINE static struct io_props io_make_output_field_convert_sink_LONGLONG(
   r.partSize = sinkSize;
   r.sinks = sinks;
   r.conversion = 1;
-  r.convert_sink_l = functionPtr;
+  r.ptr_func = functionPtr;
 
   return r;
 }
diff --git a/src/lightcone/lightcone_map.c b/src/lightcone/lightcone_map.c
index 5cf81c21de396972995855c14bc2feb72c70d48f..9bb030e24b3d11700aa2e01c105a2fed84485249 100644
--- a/src/lightcone/lightcone_map.c
+++ b/src/lightcone/lightcone_map.c
@@ -284,6 +284,8 @@ void lightcone_map_write(struct lightcone_map *map, const hid_t loc_id,
   io_write_attribute_f(dset_id, "h-scale exponent", 0.f);
   io_write_attribute_f(dset_id, "a-scale exponent", 0.f);
   io_write_attribute_s(dset_id, "Expression for physical CGS units", buffer);
+  io_write_attribute_b(dset_id, "Value stored as physical", 0);
+  io_write_attribute_b(dset_id, "Property can be converted to comoving", 1);
 
   /* Write the actual number this conversion factor corresponds to */
   const double cgs_factor =
diff --git a/src/line_of_sight.c b/src/line_of_sight.c
index 55dc0222c8eb87ccdfb221923ea929792b93c996..1537ee1a7570b61dbd7ffd886db9aa945e9eaff7 100644
--- a/src/line_of_sight.c
+++ b/src/line_of_sight.c
@@ -410,6 +410,9 @@ void write_los_hdf5_dataset(const struct io_props props, const size_t N,
   io_write_attribute_f(h_data, "a-scale exponent", props.scale_factor_exponent);
   io_write_attribute_s(h_data, "Expression for physical CGS units", buffer);
   io_write_attribute_s(h_data, "Lossy compression filter", comp_buffer);
+  io_write_attribute_b(h_data, "Value stored as physical", props.is_physical);
+  io_write_attribute_b(h_data, "Property can be converted to comoving",
+                       props.is_convertible_to_comoving);
 
   /* Write the actual number this conversion factor corresponds to */
   const double factor =
diff --git a/src/parallel_io.c b/src/parallel_io.c
index 5bdb126342255c3cf903bd9b0f3aa07624345e41..cc024dd6525f4bb4ee92bad053503b1babb6af6f 100644
--- a/src/parallel_io.c
+++ b/src/parallel_io.c
@@ -479,6 +479,9 @@ void prepare_array_parallel(
   io_write_attribute_f(h_data, "a-scale exponent", props.scale_factor_exponent);
   io_write_attribute_s(h_data, "Expression for physical CGS units", buffer);
   io_write_attribute_s(h_data, "Lossy compression filter", comp_buffer);
+  io_write_attribute_b(h_data, "Value stored as physical", props.is_physical);
+  io_write_attribute_b(h_data, "Property can be converted to comoving",
+                       props.is_convertible_to_comoving);
 
   /* Write the actual number this conversion factor corresponds to */
   const double factor =
diff --git a/src/rt/GEAR/rt_io.h b/src/rt/GEAR/rt_io.h
index 49d1096b61d29549612e345655aa4483de81b4a5..cc457bd3787fd2c7405d1dcd3c0938bb9bdb0ac0 100644
--- a/src/rt/GEAR/rt_io.h
+++ b/src/rt/GEAR/rt_io.h
@@ -161,13 +161,14 @@ INLINE static int rt_write_particles(const struct part* parts,
 
   int num_elements = 3;
 
-  list[0] = io_make_output_field_convert_part(
+  list[0] = io_make_physical_output_field_convert_part(
       "PhotonEnergies", FLOAT, RT_NGROUPS, UNIT_CONV_ENERGY, 0, parts,
       /*xparts=*/NULL, rt_convert_radiation_energies,
-      "Photon Energies (all groups)");
-  list[1] = io_make_output_field_convert_part(
+      /*convertible to comoving=*/1, "Photon Energies (all groups)");
+  list[1] = io_make_physical_output_field_convert_part(
       "PhotonFluxes", FLOAT, 3 * RT_NGROUPS, UNIT_CONV_RADIATION_FLUX, 0, parts,
       /*xparts=*/NULL, rt_convert_radiation_fluxes,
+      /*convertible to comoving=*/1,
       "Photon Fluxes (all groups; x, y, and z coordinates)");
   list[2] = io_make_output_field_convert_part(
       "IonMassFractions", FLOAT, 5, UNIT_CONV_NO_UNITS, 0, parts,
@@ -312,6 +313,8 @@ INLINE static void rt_write_flavour(hid_t h_grp, hid_t h_grp_columns,
   io_write_attribute_f(dset, "h-scale exponent", 0.f);
   io_write_attribute_f(dset, "a-scale exponent", 0.f);
   io_write_attribute_s(dset, "Expression for physical CGS units", buffer);
+  io_write_attribute_b(dset, "Value stored as physical", 1);
+  io_write_attribute_b(dset, "Property can be converted to comoving", 0);
 
   /* Write the actual number this conversion factor corresponds to */
   const double factor =
@@ -413,6 +416,8 @@ INLINE static void rt_write_flavour(hid_t h_grp, hid_t h_grp_columns,
   io_write_attribute_f(dset_cred, "a-scale exponent", 0.f);
   io_write_attribute_s(dset_cred, "Expression for physical CGS units",
                        buffer_cred);
+  io_write_attribute_b(dset_cred, "Value stored as physical", 1);
+  io_write_attribute_b(dset_cred, "Property can be converted to comoving", 0);
 
   /* Write the actual number this conversion factor corresponds to */
   /* TODO Mladen: check cosmology. reduced_speed_of_light is physical only for
diff --git a/src/rt/SPHM1RT/rt_io.h b/src/rt/SPHM1RT/rt_io.h
index c1cf3c41f521cfd29b6e8bc70c8167251aea0f1d..e19d062cd82836210c9fe20e4ac3c310396edeb5 100644
--- a/src/rt/SPHM1RT/rt_io.h
+++ b/src/rt/SPHM1RT/rt_io.h
@@ -126,10 +126,8 @@ INLINE static void rt_convert_conserved_photon_fluxes(
 INLINE static int rt_write_particles(const struct part* parts,
                                      struct io_props* list) {
 
-  /* Note that in the output, we write radiation energy and flux
-   * then we convert these quantities from radiation energy per mass and flux
-   * per mass
-   * */
+  /* Note that in the output, we write radiation energy and flux then we convert
+   * these quantities from radiation energy per mass and flux per mass */
   int num_elements = 4;
 
   list[0] = io_make_output_field_convert_part(
@@ -233,6 +231,8 @@ INLINE static void rt_write_flavour(hid_t h_grp, hid_t h_grp_columns,
   io_write_attribute_f(dset, "h-scale exponent", 0.f);
   io_write_attribute_f(dset, "a-scale exponent", 0.f);
   io_write_attribute_s(dset, "Expression for physical CGS units", buffer);
+  io_write_attribute_b(dset, "Value stored as physical", 1);
+  io_write_attribute_b(dset, "Property can be converted to comoving", 0);
 
   /* Write the actual number this conversion factor corresponds to */
   const double factor =
@@ -334,6 +334,8 @@ INLINE static void rt_write_flavour(hid_t h_grp, hid_t h_grp_columns,
   io_write_attribute_f(dset_cred, "a-scale exponent", 0.f);
   io_write_attribute_s(dset_cred, "Expression for physical CGS units",
                        buffer_cred);
+  io_write_attribute_b(dset_cred, "Value stored as physical", 1);
+  io_write_attribute_b(dset_cred, "Property can be converted to comoving", 0);
 
   /* Write the actual number this conversion factor corresponds to */
   /* TODO Mladen: check cosmology. reduced_speed_of_light is physical only for
diff --git a/src/serial_io.c b/src/serial_io.c
index f68d32b5ce406b10820b656d3b7998b462dc5880..4bd6098e31c6d5076cf105d48e912379a1cf185a 100644
--- a/src/serial_io.c
+++ b/src/serial_io.c
@@ -356,6 +356,9 @@ void prepare_array_serial(
   io_write_attribute_f(h_data, "a-scale exponent", props.scale_factor_exponent);
   io_write_attribute_s(h_data, "Expression for physical CGS units", buffer);
   io_write_attribute_s(h_data, "Lossy compression filter", comp_buffer);
+  io_write_attribute_b(h_data, "Value stored as physical", props.is_physical);
+  io_write_attribute_b(h_data, "Property can be converted to comoving",
+                       props.is_convertible_to_comoving);
 
   /* Write the actual number this conversion factor corresponds to */
   const double factor =
diff --git a/src/single_io.c b/src/single_io.c
index 8b1892af199931bbf882014c4d35413893a03905..c63be580b02106b0b9d21a70ae9a592e63ed2c89 100644
--- a/src/single_io.c
+++ b/src/single_io.c
@@ -365,6 +365,9 @@ void write_array_single(const struct engine* e, hid_t grp, const char* fileName,
   io_write_attribute_f(h_data, "a-scale exponent", props.scale_factor_exponent);
   io_write_attribute_s(h_data, "Expression for physical CGS units", buffer);
   io_write_attribute_s(h_data, "Lossy compression filter", comp_buffer);
+  io_write_attribute_b(h_data, "Value stored as physical", props.is_physical);
+  io_write_attribute_b(h_data, "Property can be converted to comoving",
+                       props.is_convertible_to_comoving);
 
   /* Write the actual number this conversion factor corresponds to */
   const double factor =
diff --git a/src/sink/Default/sink_io.h b/src/sink/Default/sink_io.h
index fd47e65d8b94f573835f0f44ceb293c8d4f4c6c5..734dbb3309b3948352d600775f6c98489657d84b 100644
--- a/src/sink/Default/sink_io.h
+++ b/src/sink/Default/sink_io.h
@@ -122,9 +122,9 @@ INLINE static void sink_write_particles(const struct sink* sinks,
   list[2] = io_make_output_field("Masses", FLOAT, 1, UNIT_CONV_MASS, 0.f, sinks,
                                  mass, "Masses of the particles");
 
-  list[3] =
-      io_make_output_field("ParticleIDs", ULONGLONG, 1, UNIT_CONV_NO_UNITS, 0.f,
-                           sinks, id, "Unique ID of the particles");
+  list[3] = io_make_physical_output_field(
+      "ParticleIDs", ULONGLONG, 1, UNIT_CONV_NO_UNITS, 0.f, sinks, id,
+      /*can convert to comoving=*/0, "Unique ID of the particles");
 
 #ifdef DEBUG_INTERACTIONS_SINKS
 
diff --git a/src/sink/GEAR/sink_io.h b/src/sink/GEAR/sink_io.h
index f5fdabcc329cf17619baf799a5f3e32c1fd69017..8ce13dfbc1b2464c1f6c47636635b5431a7f7eb8 100644
--- a/src/sink/GEAR/sink_io.h
+++ b/src/sink/GEAR/sink_io.h
@@ -122,9 +122,9 @@ INLINE static void sink_write_particles(const struct sink* sinks,
   list[2] = io_make_output_field("Masses", FLOAT, 1, UNIT_CONV_MASS, 0.f, sinks,
                                  mass, "Masses of the particles");
 
-  list[3] =
-      io_make_output_field("ParticleIDs", ULONGLONG, 1, UNIT_CONV_NO_UNITS, 0.f,
-                           sinks, id, "Unique ID of the particles");
+  list[3] = io_make_physical_output_field(
+      "ParticleIDs", ULONGLONG, 1, UNIT_CONV_NO_UNITS, 0.f, sinks, id,
+      /*can convert to comoving=*/0, "Unique ID of the particles");
 
   /* chemistry : this should be in chemistry_io.h */
   list[4] = io_make_output_field(
@@ -132,13 +132,15 @@ INLINE static void sink_write_particles(const struct sink* sinks,
       UNIT_CONV_NO_UNITS, 0.f, sinks, chemistry_data.metal_mass_fraction,
       "Mass fraction of each element");
 
-  list[5] = io_make_output_field(
+  list[5] = io_make_physical_output_field(
       "NumberOfSinkSwallows", INT, 1, UNIT_CONV_NO_UNITS, 0.f, sinks,
-      number_of_sink_swallows, "Total number of sink merger events");
+      number_of_sink_swallows, /*can convert to comoving=*/0,
+      "Total number of sink merger events");
 
-  list[6] = io_make_output_field(
+  list[6] = io_make_physical_output_field(
       "NumberOfGasSwallows", INT, 1, UNIT_CONV_NO_UNITS, 0.f, sinks,
-      number_of_gas_swallows, "Total number of gas merger events");
+      number_of_gas_swallows, /*can convert to comoving=*/0,
+      "Total number of gas merger events");
 
 #ifdef DEBUG_INTERACTIONS_SINKS
 
diff --git a/src/star_formation/GEAR/star_formation_io.h b/src/star_formation/GEAR/star_formation_io.h
index 710c085cbcbbe3aa0dce223f031fd47d7c2b9dfe..bc9bff061a7490849b1b8861b454e06d49f1f502 100644
--- a/src/star_formation/GEAR/star_formation_io.h
+++ b/src/star_formation/GEAR/star_formation_io.h
@@ -72,17 +72,17 @@ __attribute__((always_inline)) INLINE static int
 star_formation_write_sparticles(const struct spart* sparts,
                                 struct io_props* list) {
 
-  list[0] = io_make_output_field(
+  list[0] = io_make_physical_output_field(
       "BirthDensities", FLOAT, 1, UNIT_CONV_DENSITY, 0.f, sparts,
-      sf_data.birth_density,
+      sf_data.birth_density, /*can convert to comoving=*/0
       "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[1] =
-      io_make_output_field("BirthTemperatures", FLOAT, 1, UNIT_CONV_TEMPERATURE,
-                           0.f, sparts, sf_data.birth_temperature,
+    io_make_physical_output_field("BirthTemperatures", FLOAT, 1, UNIT_CONV_TEMPERATURE,
+                           0.f, sparts, sf_data.birth_temperature, /*can convert to comoving=*/0
                            "Temperatures at the time of birth of the gas "
                            "particles that turned into stars");
 
diff --git a/src/stars/Basic/stars_io.h b/src/stars/Basic/stars_io.h
index 54ad764954bd42d199ae0d67550cbedb3f855b86..1b018459a964991dd0cd2e5e3e0643cdb70f1f6c 100644
--- a/src/stars/Basic/stars_io.h
+++ b/src/stars/Basic/stars_io.h
@@ -142,9 +142,9 @@ INLINE static void stars_write_particles(const struct spart *sparts,
   list[2] = io_make_output_field("Masses", FLOAT, 1, UNIT_CONV_MASS, 0.f,
                                  sparts, mass, "Masses of the particles");
 
-  list[3] =
-      io_make_output_field("ParticleIDs", ULONGLONG, 1, UNIT_CONV_NO_UNITS, 0.f,
-                           sparts, id, "Unique ID of the particles");
+  list[3] = io_make_physical_output_field(
+      "ParticleIDs", ULONGLONG, 1, UNIT_CONV_NO_UNITS, 0.f, sparts, id,
+      /*can convert to comoving=*/0, "Unique ID of the particles");
 
   list[4] = io_make_output_field(
       "SmoothingLengths", FLOAT, 1, UNIT_CONV_LENGTH, 1.f, sparts, h,
diff --git a/src/stars/EAGLE/stars_io.h b/src/stars/EAGLE/stars_io.h
index 1c499624a67e995e106d4c89225a1e076d8ec21c..cb93e32736b38cc2ceb9542c55cd8ca420896d1c 100644
--- a/src/stars/EAGLE/stars_io.h
+++ b/src/stars/EAGLE/stars_io.h
@@ -160,9 +160,9 @@ INLINE static void stars_write_particles(const struct spart *sparts,
                                  "Masses of the particles at the current point "
                                  "in time (i.e. after stellar losses");
 
-  list[3] =
-      io_make_output_field("ParticleIDs", ULONGLONG, 1, UNIT_CONV_NO_UNITS, 0.f,
-                           sparts, id, "Unique ID of the particles");
+  list[3] = io_make_physical_output_field(
+      "ParticleIDs", ULONGLONG, 1, UNIT_CONV_NO_UNITS, 0.f, sparts, id,
+      /*can convert to comoving=*/0, "Unique ID of the particles");
 
   list[4] = io_make_output_field(
       "SmoothingLengths", FLOAT, 1, UNIT_CONV_LENGTH, 1.f, sparts, h,
@@ -173,9 +173,10 @@ INLINE static void stars_write_particles(const struct spart *sparts,
                                  "Masses of the star particles at birth time");
 
   if (with_cosmology) {
-    list[6] = io_make_output_field(
+    list[6] = io_make_physical_output_field(
         "BirthScaleFactors", FLOAT, 1, UNIT_CONV_NO_UNITS, 0.f, sparts,
-        birth_scale_factor, "Scale-factors at which the stars were born");
+        birth_scale_factor, /*can convert to comoving=*/0,
+        "Scale-factors at which the stars were born");
   } else {
     list[6] = io_make_output_field("BirthTimes", FLOAT, 1, UNIT_CONV_TIME, 0.f,
                                    sparts, birth_time,
@@ -192,18 +193,18 @@ INLINE static void stars_write_particles(const struct spart *sparts,
       number_of_SNII_events,
       "Number of SNII energy injection events the stars went through.");
 
-  list[9] = io_make_output_field(
-      "BirthDensities", FLOAT, 1, UNIT_CONV_DENSITY, 0.f, sparts, birth_density,
+  list[9] = io_make_physical_output_field(
+      "BirthDensities", FLOAT, 1, UNIT_CONV_DENSITY, -3.f, sparts,
+      birth_density, /*can convert to comoving=*/0,
       "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[10] =
-      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");
+      "turned into stars (note that we store the physical density at the birth "
+      "redshift, no conversion is needed)");
+
+  list[10] = io_make_physical_output_field(
+      "BirthTemperatures", FLOAT, 1, UNIT_CONV_TEMPERATURE, 0.f, sparts,
+      birth_temperature, /*can convert to comoving=*/0,
+      "Temperatures at the time of birth of the gas "
+      "particles that turned into stars");
 
   list[11] = io_make_output_field(
       "FeedbackNumberOfHeatingEvents", FLOAT, 1, UNIT_CONV_NO_UNITS, 0.f,
diff --git a/src/stars/GEAR/stars_io.h b/src/stars/GEAR/stars_io.h
index 041551d96c2c545ce8bb51eb2fd89e0e16e27362..3d9c5ce10da6cfc03ec9ec10bd4dafe78d6c4efe 100644
--- a/src/stars/GEAR/stars_io.h
+++ b/src/stars/GEAR/stars_io.h
@@ -144,18 +144,19 @@ INLINE static void stars_write_particles(const struct spart *sparts,
   list[2] = io_make_output_field("Masses", FLOAT, 1, UNIT_CONV_MASS, 0.f,
                                  sparts, mass, "Masses of the particles");
 
-  list[3] =
-      io_make_output_field("ParticleIDs", LONGLONG, 1, UNIT_CONV_NO_UNITS, 0.f,
-                           sparts, id, "Unique IDs of the particles");
+  list[3] = io_make_physical_output_field(
+      "ParticleIDs", LONGLONG, 1, UNIT_CONV_NO_UNITS, 0.f, sparts, id,
+      /*can convert to comoving=*/0, "Unique IDs 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");
 
   if (with_cosmology) {
-    list[5] = io_make_output_field(
+    list[5] = io_make_physical_output_field(
         "BirthScaleFactors", FLOAT, 1, UNIT_CONV_NO_UNITS, 0.f, sparts,
-        birth_scale_factor, "Scale-factors at which the stars were born");
+        birth_scale_factor, /*can convert to comoving=*/0,
+        "Scale-factors at which the stars were born");
   } else {
     list[5] = io_make_output_field("BirthTimes", FLOAT, 1, UNIT_CONV_TIME, 0.f,
                                    sparts, birth_time,
diff --git a/src/tracers/EAGLE/tracers_io.h b/src/tracers/EAGLE/tracers_io.h
index 5971ed99e81e2b721d7ff6ba118a0b2e004b0d06..cc792109239608c06b1548d48c56815ab1608cab 100644
--- a/src/tracers/EAGLE/tracers_io.h
+++ b/src/tracers/EAGLE/tracers_io.h
@@ -104,9 +104,10 @@ __attribute__((always_inline)) INLINE static int tracers_write_particles(
       "Maximal temperatures ever reached by the particles");
 
   if (with_cosmology) {
-    list[1] = io_make_output_field(
+    list[1] = io_make_physical_output_field(
         "MaximalTemperatureScaleFactors", FLOAT, 1, UNIT_CONV_NO_UNITS, 0.f,
         xparts, tracers_data.maximum_temperature_scale_factor,
+        /*can convert to comoving=*/0,
         "Scale-factors at which the maximal temperature was reached");
 
   } else {
@@ -137,39 +138,44 @@ __attribute__((always_inline)) INLINE static int tracers_write_particles(
                                  "Total amount of thermal energy from AGN "
                                  "feedback events received by the particles.");
 
-  list[5] = io_make_output_field(
+  list[5] = io_make_physical_output_field(
       "DensitiesBeforeLastAGNEvent", FLOAT, 1, UNIT_CONV_DENSITY, 0.f, xparts,
       tracers_data.density_before_last_AGN_feedback_event,
+      /*can convert to comoving=*/0,
       "Physical density (not subgrid) of the gas fetched before the last AGN "
       "feedback event that hit the particles. -1 if the particles have never "
       "been heated.");
 
-  list[6] = io_make_output_field(
+  list[6] = io_make_physical_output_field(
       "EntropiesBeforeLastAGNEvent", FLOAT, 1, UNIT_CONV_ENTROPY_PER_UNIT_MASS,
       0.f, xparts, tracers_data.entropy_before_last_AGN_feedback_event,
+      /*can convert to comoving=*/0,
       "Physical entropy (not subgrid) per unit mass of the gas fetched before "
       "the last AGN feedback event that hit the particles. -1 if the particles "
       "have never been heated.");
 
-  list[7] = io_make_output_field(
+  list[7] = io_make_physical_output_field(
       "DensitiesAtLastAGNEvent", FLOAT, 1, UNIT_CONV_DENSITY, 0.f, xparts,
       tracers_data.density_at_last_AGN_feedback_event,
+      /*can convert to comoving=*/0,
       "Physical density (not subgrid) of the gas at the last AGN feedback "
       "event that hit the particles. -1 if the particles have never been "
       "heated.");
 
-  list[8] = io_make_output_field(
+  list[8] = io_make_physical_output_field(
       "EntropiesAtLastAGNEvent", FLOAT, 1, UNIT_CONV_ENTROPY_PER_UNIT_MASS, 0.f,
       xparts, tracers_data.entropy_at_last_AGN_feedback_event,
+      /*can convert to comoving=*/0,
       "Physical entropy (not subgrid) per unit mass of the gas at the last AGN "
       "feedback event that hit the particles. -1 if the particles have never "
       "been heated.");
 
   if (with_cosmology) {
 
-    list[9] = io_make_output_field(
+    list[9] = io_make_physical_output_field(
         "LastAGNFeedbackScaleFactors", FLOAT, 1, UNIT_CONV_NO_UNITS, 0.f,
         xparts, tracers_data.last_AGN_injection_scale_factor,
+        /*can convert to comoving=*/0,
         "Scale-factors at which the particles were last hit by AGN feedback. "
         "-1 if a particle has never been hit by feedback");
 
@@ -196,19 +202,20 @@ __attribute__((always_inline)) INLINE static int tracers_write_particles(
         "an AGN jet feedback event at some point in the past. "
         "If > 0, contains the number of individual events.");
 
-    list[12] = io_make_output_field("EnergiesReceivedFromJetFeedback", FLOAT, 1,
-                                    UNIT_CONV_ENERGY, 0.f, xparts,
-                                    tracers_data.jet_feedback_energy,
-                                    "Total amount of kinetic energy from AGN "
-                                    "jet feedback events received by the "
-                                    "particles while they were still gas "
-                                    "particles.");
+    list[12] = io_make_physical_output_field(
+        "EnergiesReceivedFromJetFeedback", FLOAT, 1, UNIT_CONV_ENERGY, 0.f,
+        xparts, tracers_data.jet_feedback_energy, /*can convert to comoving=*/0,
+        "Total amount of kinetic energy from AGN "
+        "jet feedback events received by the "
+        "particles while they were still gas "
+        "particles.");
 
     if (with_cosmology) {
 
-      list[13] = io_make_output_field(
+      list[13] = io_make_physical_output_field(
           "LastAGNJetFeedbackScaleFactors", FLOAT, 1, UNIT_CONV_NO_UNITS, 0.f,
           xparts, tracers_data.last_AGN_jet_feedback_scale_factor,
+          /*can convert to comoving=*/0,
           "Scale-factors at which the particles were last hit by jet "
           "feedback while they were still gas particles. "
           "-1 if a particle has never been hit by feedback");
@@ -259,9 +266,10 @@ __attribute__((always_inline)) INLINE static int tracers_write_sparticles(
       "converted to stars");
 
   if (with_cosmology) {
-    list[1] = io_make_output_field(
+    list[1] = io_make_physical_output_field(
         "MaximalTemperatureScaleFactors", FLOAT, 1, UNIT_CONV_NO_UNITS, 0.f,
         sparts, tracers_data.maximum_temperature_scale_factor,
+        /*can convert to comoving=*/0,
         "Scale-factors at which the maximal temperature was reached");
 
   } else {
@@ -286,39 +294,43 @@ __attribute__((always_inline)) INLINE static int tracers_write_sparticles(
                            "an AGN feedback event at some point in the past "
                            "when the particle was still a gas particle.");
 
-  list[4] = io_make_output_field(
+  list[4] = io_make_physical_output_field(
       "EnergiesReceivedFromAGNFeedback", FLOAT, 1, UNIT_CONV_ENERGY, 0.f,
-      sparts, tracers_data.AGN_feedback_energy,
+      sparts, tracers_data.AGN_feedback_energy, /*can convert to comoving=*/0,
       "Total amount of thermal energy from AGN feedback events received by the "
       "particles when the particle was still a gas particle.");
 
-  list[5] = io_make_output_field(
+  list[5] = io_make_physical_output_field(
       "DensitiesBeforeLastAGNEvent", FLOAT, 1, UNIT_CONV_DENSITY, 0.f, sparts,
       tracers_data.density_before_last_AGN_feedback_event,
+      /*can convert to comoving=*/0,
       "Physical density (not subgrid) of the gas fetched before the last AGN "
       "feedback "
       "event that hit the particles when they were still gas particles. -1 if "
       "the particles have never been heated.");
 
-  list[6] = io_make_output_field(
+  list[6] = io_make_physical_output_field(
       "EntropiesBeforeLastAGNEvent", FLOAT, 1, UNIT_CONV_ENTROPY_PER_UNIT_MASS,
       0.f, sparts, tracers_data.entropy_before_last_AGN_feedback_event,
+      /*can convert to comoving=*/0,
       "Physical entropy (not subgrid) per unit mass of the gas fetched before "
       "the last AGN "
       "feedback event that hit the particles when they were still gas "
       "particles."
       " -1 if the particles have never been heated.");
 
-  list[7] = io_make_output_field(
+  list[7] = io_make_physical_output_field(
       "DensitiesAtLastAGNEvent", FLOAT, 1, UNIT_CONV_DENSITY, 0.f, sparts,
       tracers_data.density_at_last_AGN_feedback_event,
+      /*can convert to comoving=*/0,
       "Physical density (not subgrid) of the gas at the last AGN feedback "
       "event that hit the particles when they were still gas particles. -1 if "
       "the particles have never been heated.");
 
-  list[8] = io_make_output_field(
+  list[8] = io_make_physical_output_field(
       "EntropiesAtLastAGNEvent", FLOAT, 1, UNIT_CONV_ENTROPY_PER_UNIT_MASS, 0.f,
       sparts, tracers_data.entropy_at_last_AGN_feedback_event,
+      /*can convert to comoving=*/0,
       "Physical entropy (not subgrid) per unit mass of the gas at the last AGN "
       "feedback event that hit the particles when they were still gas "
       "particles."
@@ -326,9 +338,10 @@ __attribute__((always_inline)) INLINE static int tracers_write_sparticles(
 
   if (with_cosmology) {
 
-    list[9] = io_make_output_field(
+    list[9] = io_make_physical_output_field(
         "LastAGNFeedbackScaleFactors", FLOAT, 1, UNIT_CONV_NO_UNITS, 0.f,
         sparts, tracers_data.last_AGN_injection_scale_factor,
+        /*can convert to comoving=*/0,
         "Scale-factors at which the particles were last hit by AGN feedback "
         "when they were still gas particles. -1 if a particle has never been "
         "hit by feedback");
@@ -368,9 +381,10 @@ __attribute__((always_inline)) INLINE static int tracers_write_sparticles(
 
     if (with_cosmology) {
 
-      list[13] = io_make_output_field(
+      list[13] = io_make_physical_output_field(
           "LastAGNJetFeedbackScaleFactors", FLOAT, 1, UNIT_CONV_NO_UNITS, 0.f,
           sparts, tracers_data.last_AGN_jet_feedback_scale_factor,
+          /*can convert to comoving=*/0,
           "Scale-factors at which the particles were last hit by jet "
           "feedback while they were still gas particles. "
           "-1 if a particle has never been hit by feedback");
diff --git a/src/velociraptor_interface.c b/src/velociraptor_interface.c
index 3bd1e32252c70766893c2822f42e9476c394046e..4248bb002b07192cb6c72417107058ccfb45928f 100644
--- a/src/velociraptor_interface.c
+++ b/src/velociraptor_interface.c
@@ -504,6 +504,9 @@ void write_orphan_particle_array(hid_t h_file, const char *name,
   io_write_attribute_f(dset_id, "a-scale exponent",
                        props.scale_factor_exponent);
   io_write_attribute_s(dset_id, "Expression for physical CGS units", buffer);
+  io_write_attribute_b(h_data, "Value stored as physical", props.is_physical);
+  io_write_attribute_b(h_data, "Property can be converted to comoving",
+                       props.is_convertible_to_comoving);
 
   /* Write data, if there is any */
   if (nr_flagged_all > 0) {