diff --git a/examples/EAGLE_ICs/EAGLE_100/eagle_100.yml b/examples/EAGLE_ICs/EAGLE_100/eagle_100.yml
index db95e5de66829fd5e29e358f2b96ff00caebbb7f..d81764896aaed82d8562fc759658a10359af5945 100644
--- a/examples/EAGLE_ICs/EAGLE_100/eagle_100.yml
+++ b/examples/EAGLE_ICs/EAGLE_100/eagle_100.yml
@@ -138,20 +138,21 @@ COLIBRECooling:
 
 # EAGLE star formation parameters
 EAGLEStarFormation:
-  SF_model:                          PressureLaw
-  EOS_density_norm_H_p_cm3:          0.1       # Physical density used for the normalisation of the EOS assumed for the star-forming gas in Hydrogen atoms per cm^3.
-  EOS_temperature_norm_K:            8000      # Temperature om the polytropic EOS assumed for star-forming gas at the density normalisation in Kelvin.
-  EOS_gamma_effective:               1.3333333 # Slope the of the polytropic EOS assumed for the star-forming gas.
-  KS_normalisation:                  1.515e-4  # The normalization of the Kennicutt-Schmidt law in Msun / kpc^2 / yr.
-  KS_exponent:                       1.4       # The exponent of the Kennicutt-Schmidt law.
-  min_over_density:                  100.0     # The over-density above which star-formation is allowed.
-  KS_high_density_threshold_H_p_cm3: 1e3       # Hydrogen number density above which the Kennicut-Schmidt law changes slope in Hydrogen atoms per cm^3.
-  KS_high_density_exponent:          2.0       # Slope of the Kennicut-Schmidt law above the high-density threshold.
-  EOS_entropy_margin_dex:            0.3       # Logarithm base 10 of the maximal entropy above the EOS at which stars can form.
-  threshold_norm_H_p_cm3:            0.1       # Normalisation of the metal-dependant density threshold for star formation in Hydrogen atoms per cm^3.
-  threshold_Z0:                      0.002     # Reference metallicity (metal mass fraction) for the metal-dependant threshold for star formation.
-  threshold_slope:                   -0.64     # Slope of the metal-dependant star formation threshold
-  threshold_max_density_H_p_cm3:     10.0      # Maximal density of the metal-dependant density threshold for star formation in Hydrogen atoms per cm^3.
+  SF_threshold:                      Subgrid      # Zdep (Schaye 2004) or Subgrid
+  SF_model:                          PressureLaw  # PressureLaw (Schaye et al. 2008) or SchmidtLaw
+  KS_normalisation:                  1.515e-4     # The normalization of the Kennicutt-Schmidt law in Msun / kpc^2 / yr.
+  KS_exponent:                       1.4          # The exponent of the Kennicutt-Schmidt law.
+  min_over_density:                  100.0        # The over-density above which star-formation is allowed.
+  KS_high_density_threshold_H_p_cm3: 1e3          # Hydrogen number density above which the Kennicut-Schmidt law changes slope in Hydrogen atoms per cm^3.
+  KS_high_density_exponent:          2.0          # Slope of the Kennicut-Schmidt law above the high-density threshold.
+  EOS_entropy_margin_dex:            0.3          # When using Z-based SF threshold, logarithm base 10 of the maximal entropy above the EOS at which stars can form.
+  threshold_norm_H_p_cm3:            0.1          # When using Z-based SF threshold, normalisation of the metal-dependant density threshold for star formation in Hydrogen atoms per cm^3.
+  threshold_Z0:                      0.002        # When using Z-based SF threshold, reference metallicity (metal mass fraction) for the metal-dependant threshold for star formation.
+  threshold_slope:                   -0.64        # When using Z-based SF threshold, slope of the metal-dependant star formation threshold
+  threshold_max_density_H_p_cm3:     10.0         # When using Z-based SF threshold, maximal density of the metal-dependant density threshold for star formation in Hydrogen atoms per cm^3.
+  threshold_temperature1_K:          1000         # When using subgrid-based SF threshold, subgrid temperature below which gas is star-forming.
+  threshold_temperature2_K:          31622        # When using subgrid-based SF threshold, subgrid temperature below which gas is star-forming if also above the density limit.
+  threshold_number_density_H_p_cm3:  10           # When using subgrid-based SF threshold, subgrid number density above which gas is star-forming if also below the second temperature limit.
   
 # Parameters for the EAGLE "equation of state"
 EAGLEEntropyFloor:
diff --git a/examples/EAGLE_ICs/EAGLE_12/eagle_12.yml b/examples/EAGLE_ICs/EAGLE_12/eagle_12.yml
index 84392328cc2e1f9e9fbd4e282e808e395dfd5020..4029cd32263e1000c2ad02a0d67432f31f9cbba6 100644
--- a/examples/EAGLE_ICs/EAGLE_12/eagle_12.yml
+++ b/examples/EAGLE_ICs/EAGLE_12/eagle_12.yml
@@ -137,20 +137,21 @@ COLIBRECooling:
 
 # EAGLE star formation parameters
 EAGLEStarFormation:
-  SF_model:                          PressureLaw
-  EOS_density_norm_H_p_cm3:          0.1       # Physical density used for the normalisation of the EOS assumed for the star-forming gas in Hydrogen atoms per cm^3.
-  EOS_temperature_norm_K:            8000      # Temperature om the polytropic EOS assumed for star-forming gas at the density normalisation in Kelvin.
-  EOS_gamma_effective:               1.3333333 # Slope the of the polytropic EOS assumed for the star-forming gas.
-  KS_normalisation:                  1.515e-4  # The normalization of the Kennicutt-Schmidt law in Msun / kpc^2 / yr.
-  KS_exponent:                       1.4       # The exponent of the Kennicutt-Schmidt law.
-  min_over_density:                  100.0     # The over-density above which star-formation is allowed.
-  KS_high_density_threshold_H_p_cm3: 1e3       # Hydrogen number density above which the Kennicut-Schmidt law changes slope in Hydrogen atoms per cm^3.
-  KS_high_density_exponent:          2.0       # Slope of the Kennicut-Schmidt law above the high-density threshold.
-  EOS_entropy_margin_dex:            0.3       # Logarithm base 10 of the maximal entropy above the EOS at which stars can form.
-  threshold_norm_H_p_cm3:            0.1       # Normalisation of the metal-dependant density threshold for star formation in Hydrogen atoms per cm^3.
-  threshold_Z0:                      0.002     # Reference metallicity (metal mass fraction) for the metal-dependant threshold for star formation.
-  threshold_slope:                   -0.64     # Slope of the metal-dependant star formation threshold
-  threshold_max_density_H_p_cm3:     10.0      # Maximal density of the metal-dependant density threshold for star formation in Hydrogen atoms per cm^3.
+  SF_threshold:                      Subgrid      # Zdep (Schaye 2004) or Subgrid
+  SF_model:                          PressureLaw  # PressureLaw (Schaye et al. 2008) or SchmidtLaw
+  KS_normalisation:                  1.515e-4     # The normalization of the Kennicutt-Schmidt law in Msun / kpc^2 / yr.
+  KS_exponent:                       1.4          # The exponent of the Kennicutt-Schmidt law.
+  min_over_density:                  100.0        # The over-density above which star-formation is allowed.
+  KS_high_density_threshold_H_p_cm3: 1e3          # Hydrogen number density above which the Kennicut-Schmidt law changes slope in Hydrogen atoms per cm^3.
+  KS_high_density_exponent:          2.0          # Slope of the Kennicut-Schmidt law above the high-density threshold.
+  EOS_entropy_margin_dex:            0.3          # When using Z-based SF threshold, logarithm base 10 of the maximal entropy above the EOS at which stars can form.
+  threshold_norm_H_p_cm3:            0.1          # When using Z-based SF threshold, normalisation of the metal-dependant density threshold for star formation in Hydrogen atoms per cm^3.
+  threshold_Z0:                      0.002        # When using Z-based SF threshold, reference metallicity (metal mass fraction) for the metal-dependant threshold for star formation.
+  threshold_slope:                   -0.64        # When using Z-based SF threshold, slope of the metal-dependant star formation threshold
+  threshold_max_density_H_p_cm3:     10.0         # When using Z-based SF threshold, maximal density of the metal-dependant density threshold for star formation in Hydrogen atoms per cm^3.
+  threshold_temperature1_K:          1000         # When using subgrid-based SF threshold, subgrid temperature below which gas is star-forming.
+  threshold_temperature2_K:          31622        # When using subgrid-based SF threshold, subgrid temperature below which gas is star-forming if also above the density limit.
+  threshold_number_density_H_p_cm3:  10           # When using subgrid-based SF threshold, subgrid number density above which gas is star-forming if also below the second temperature limit.
   
 # Parameters for the EAGLE "equation of state"
 EAGLEEntropyFloor:
diff --git a/examples/EAGLE_ICs/EAGLE_25/eagle_25.yml b/examples/EAGLE_ICs/EAGLE_25/eagle_25.yml
index 7cd0437ef7020b0c9ea081bd4465512896742c9d..a52590c9d761da9d8a167e73cbc6f0a8e8a5979f 100644
--- a/examples/EAGLE_ICs/EAGLE_25/eagle_25.yml
+++ b/examples/EAGLE_ICs/EAGLE_25/eagle_25.yml
@@ -137,20 +137,21 @@ COLIBRECooling:
 
 # EAGLE star formation parameters
 EAGLEStarFormation:
-  SF_model:                          PressureLaw
-  EOS_density_norm_H_p_cm3:          0.1       # Physical density used for the normalisation of the EOS assumed for the star-forming gas in Hydrogen atoms per cm^3.
-  EOS_temperature_norm_K:            8000      # Temperature om the polytropic EOS assumed for star-forming gas at the density normalisation in Kelvin.
-  EOS_gamma_effective:               1.3333333 # Slope the of the polytropic EOS assumed for the star-forming gas.
-  KS_normalisation:                  1.515e-4  # The normalization of the Kennicutt-Schmidt law in Msun / kpc^2 / yr.
-  KS_exponent:                       1.4       # The exponent of the Kennicutt-Schmidt law.
-  min_over_density:                  100.0     # The over-density above which star-formation is allowed.
-  KS_high_density_threshold_H_p_cm3: 1e3       # Hydrogen number density above which the Kennicut-Schmidt law changes slope in Hydrogen atoms per cm^3.
-  KS_high_density_exponent:          2.0       # Slope of the Kennicut-Schmidt law above the high-density threshold.
-  EOS_entropy_margin_dex:            0.3       # Logarithm base 10 of the maximal entropy above the EOS at which stars can form.
-  threshold_norm_H_p_cm3:            0.1       # Normalisation of the metal-dependant density threshold for star formation in Hydrogen atoms per cm^3.
-  threshold_Z0:                      0.002     # Reference metallicity (metal mass fraction) for the metal-dependant threshold for star formation.
-  threshold_slope:                   -0.64     # Slope of the metal-dependant star formation threshold
-  threshold_max_density_H_p_cm3:     10.0      # Maximal density of the metal-dependant density threshold for star formation in Hydrogen atoms per cm^3.
+  SF_threshold:                      Subgrid      # Zdep (Schaye 2004) or Subgrid
+  SF_model:                          PressureLaw  # PressureLaw (Schaye et al. 2008) or SchmidtLaw
+  KS_normalisation:                  1.515e-4     # The normalization of the Kennicutt-Schmidt law in Msun / kpc^2 / yr.
+  KS_exponent:                       1.4          # The exponent of the Kennicutt-Schmidt law.
+  min_over_density:                  100.0        # The over-density above which star-formation is allowed.
+  KS_high_density_threshold_H_p_cm3: 1e3          # Hydrogen number density above which the Kennicut-Schmidt law changes slope in Hydrogen atoms per cm^3.
+  KS_high_density_exponent:          2.0          # Slope of the Kennicut-Schmidt law above the high-density threshold.
+  EOS_entropy_margin_dex:            0.3          # When using Z-based SF threshold, logarithm base 10 of the maximal entropy above the EOS at which stars can form.
+  threshold_norm_H_p_cm3:            0.1          # When using Z-based SF threshold, normalisation of the metal-dependant density threshold for star formation in Hydrogen atoms per cm^3.
+  threshold_Z0:                      0.002        # When using Z-based SF threshold, reference metallicity (metal mass fraction) for the metal-dependant threshold for star formation.
+  threshold_slope:                   -0.64        # When using Z-based SF threshold, slope of the metal-dependant star formation threshold
+  threshold_max_density_H_p_cm3:     10.0         # When using Z-based SF threshold, maximal density of the metal-dependant density threshold for star formation in Hydrogen atoms per cm^3.
+  threshold_temperature1_K:          1000         # When using subgrid-based SF threshold, subgrid temperature below which gas is star-forming.
+  threshold_temperature2_K:          31622        # When using subgrid-based SF threshold, subgrid temperature below which gas is star-forming if also above the density limit.
+  threshold_number_density_H_p_cm3:  10           # When using subgrid-based SF threshold, subgrid number density above which gas is star-forming if also below the second temperature limit.
   
 # Parameters for the EAGLE "equation of state"
 EAGLEEntropyFloor:
diff --git a/examples/EAGLE_ICs/EAGLE_25_low_res/eagle_25.yml b/examples/EAGLE_ICs/EAGLE_25_low_res/eagle_25.yml
index 06b80444506ff17c4c645a9b59de041c723427a0..e37450795a9847fbb1a89d0ae92220372873094a 100644
--- a/examples/EAGLE_ICs/EAGLE_25_low_res/eagle_25.yml
+++ b/examples/EAGLE_ICs/EAGLE_25_low_res/eagle_25.yml
@@ -134,20 +134,21 @@ COLIBRECooling:
 
 # EAGLE star formation parameters
 EAGLEStarFormation:
-  SF_model:                          PressureLaw
-  EOS_density_norm_H_p_cm3:          0.1       # Physical density used for the normalisation of the EOS assumed for the star-forming gas in Hydrogen atoms per cm^3.
-  EOS_temperature_norm_K:            8000      # Temperature om the polytropic EOS assumed for star-forming gas at the density normalisation in Kelvin.
-  EOS_gamma_effective:               1.3333333 # Slope the of the polytropic EOS assumed for the star-forming gas.
-  KS_normalisation:                  1.515e-4  # The normalization of the Kennicutt-Schmidt law in Msun / kpc^2 / yr.
-  KS_exponent:                       1.4       # The exponent of the Kennicutt-Schmidt law.
-  min_over_density:                  100.0     # The over-density above which star-formation is allowed.
-  KS_high_density_threshold_H_p_cm3: 1e3       # Hydrogen number density above which the Kennicut-Schmidt law changes slope in Hydrogen atoms per cm^3.
-  KS_high_density_exponent:          2.0       # Slope of the Kennicut-Schmidt law above the high-density threshold.
-  EOS_entropy_margin_dex:            0.3       # Logarithm base 10 of the maximal entropy above the EOS at which stars can form.
-  threshold_norm_H_p_cm3:            0.1       # Normalisation of the metal-dependant density threshold for star formation in Hydrogen atoms per cm^3.
-  threshold_Z0:                      0.002     # Reference metallicity (metal mass fraction) for the metal-dependant threshold for star formation.
-  threshold_slope:                   -0.64     # Slope of the metal-dependant star formation threshold
-  threshold_max_density_H_p_cm3:     10.0      # Maximal density of the metal-dependant density threshold for star formation in Hydrogen atoms per cm^3.
+  SF_threshold:                      Subgrid      # Zdep (Schaye 2004) or Subgrid
+  SF_model:                          PressureLaw  # PressureLaw (Schaye et al. 2008) or SchmidtLaw
+  KS_normalisation:                  1.515e-4     # The normalization of the Kennicutt-Schmidt law in Msun / kpc^2 / yr.
+  KS_exponent:                       1.4          # The exponent of the Kennicutt-Schmidt law.
+  min_over_density:                  100.0        # The over-density above which star-formation is allowed.
+  KS_high_density_threshold_H_p_cm3: 1e3          # Hydrogen number density above which the Kennicut-Schmidt law changes slope in Hydrogen atoms per cm^3.
+  KS_high_density_exponent:          2.0          # Slope of the Kennicut-Schmidt law above the high-density threshold.
+  EOS_entropy_margin_dex:            0.3          # When using Z-based SF threshold, logarithm base 10 of the maximal entropy above the EOS at which stars can form.
+  threshold_norm_H_p_cm3:            0.1          # When using Z-based SF threshold, normalisation of the metal-dependant density threshold for star formation in Hydrogen atoms per cm^3.
+  threshold_Z0:                      0.002        # When using Z-based SF threshold, reference metallicity (metal mass fraction) for the metal-dependant threshold for star formation.
+  threshold_slope:                   -0.64        # When using Z-based SF threshold, slope of the metal-dependant star formation threshold
+  threshold_max_density_H_p_cm3:     10.0         # When using Z-based SF threshold, maximal density of the metal-dependant density threshold for star formation in Hydrogen atoms per cm^3.
+  threshold_temperature1_K:          1000         # When using subgrid-based SF threshold, subgrid temperature below which gas is star-forming.
+  threshold_temperature2_K:          31622        # When using subgrid-based SF threshold, subgrid temperature below which gas is star-forming if also above the density limit.
+  threshold_number_density_H_p_cm3:  10           # When using subgrid-based SF threshold, subgrid number density above which gas is star-forming if also below the second temperature limit.
   
 # Parameters for the EAGLE "equation of state"
 EAGLEEntropyFloor:
diff --git a/examples/EAGLE_ICs/EAGLE_50/eagle_50.yml b/examples/EAGLE_ICs/EAGLE_50/eagle_50.yml
index 2613e3b49d339b9fb9f63ef7dd0b07992930d85e..b7ec98015f93782e315c8ba3d0d7c53f721ccf75 100644
--- a/examples/EAGLE_ICs/EAGLE_50/eagle_50.yml
+++ b/examples/EAGLE_ICs/EAGLE_50/eagle_50.yml
@@ -135,21 +135,22 @@ COLIBRECooling:
 
 # EAGLE star formation parameters
 EAGLEStarFormation:
-  SF_model:                          PressureLaw
-  EOS_density_norm_H_p_cm3:          0.1       # Physical density used for the normalisation of the EOS assumed for the star-forming gas in Hydrogen atoms per cm^3.
-  EOS_temperature_norm_K:            8000      # Temperature om the polytropic EOS assumed for star-forming gas at the density normalisation in Kelvin.
-  EOS_gamma_effective:               1.3333333 # Slope the of the polytropic EOS assumed for the star-forming gas.
-  KS_normalisation:                  1.515e-4  # The normalization of the Kennicutt-Schmidt law in Msun / kpc^2 / yr.
-  KS_exponent:                       1.4       # The exponent of the Kennicutt-Schmidt law.
-  min_over_density:                  100.0     # The over-density above which star-formation is allowed.
-  KS_high_density_threshold_H_p_cm3: 1e3       # Hydrogen number density above which the Kennicut-Schmidt law changes slope in Hydrogen atoms per cm^3.
-  KS_high_density_exponent:          2.0       # Slope of the Kennicut-Schmidt law above the high-density threshold.
-  EOS_entropy_margin_dex:            0.3       # Logarithm base 10 of the maximal entropy above the EOS at which stars can form.
-  threshold_norm_H_p_cm3:            0.1       # Normalisation of the metal-dependant density threshold for star formation in Hydrogen atoms per cm^3.
-  threshold_Z0:                      0.002     # Reference metallicity (metal mass fraction) for the metal-dependant threshold for star formation.
-  threshold_slope:                   -0.64     # Slope of the metal-dependant star formation threshold
-  threshold_max_density_H_p_cm3:     10.0      # Maximal density of the metal-dependant density threshold for star formation in Hydrogen atoms per cm^3.
-  
+  SF_threshold:                      Subgrid      # Zdep (Schaye 2004) or Subgrid
+  SF_model:                          PressureLaw  # PressureLaw (Schaye et al. 2008) or SchmidtLaw
+  KS_normalisation:                  1.515e-4     # The normalization of the Kennicutt-Schmidt law in Msun / kpc^2 / yr.
+  KS_exponent:                       1.4          # The exponent of the Kennicutt-Schmidt law.
+  min_over_density:                  100.0        # The over-density above which star-formation is allowed.
+  KS_high_density_threshold_H_p_cm3: 1e3          # Hydrogen number density above which the Kennicut-Schmidt law changes slope in Hydrogen atoms per cm^3.
+  KS_high_density_exponent:          2.0          # Slope of the Kennicut-Schmidt law above the high-density threshold.
+  EOS_entropy_margin_dex:            0.3          # When using Z-based SF threshold, logarithm base 10 of the maximal entropy above the EOS at which stars can form.
+  threshold_norm_H_p_cm3:            0.1          # When using Z-based SF threshold, normalisation of the metal-dependant density threshold for star formation in Hydrogen atoms per cm^3.
+  threshold_Z0:                      0.002        # When using Z-based SF threshold, reference metallicity (metal mass fraction) for the metal-dependant threshold for star formation.
+  threshold_slope:                   -0.64        # When using Z-based SF threshold, slope of the metal-dependant star formation threshold
+  threshold_max_density_H_p_cm3:     10.0         # When using Z-based SF threshold, maximal density of the metal-dependant density threshold for star formation in Hydrogen atoms per cm^3.
+  threshold_temperature1_K:          1000         # When using subgrid-based SF threshold, subgrid temperature below which gas is star-forming.
+  threshold_temperature2_K:          31622        # When using subgrid-based SF threshold, subgrid temperature below which gas is star-forming if also above the density limit.
+  threshold_number_density_H_p_cm3:  10           # When using subgrid-based SF threshold, subgrid number density above which gas is star-forming if also below the second temperature limit.
+
 # Parameters for the EAGLE "equation of state"
 EAGLEEntropyFloor:
   Jeans_density_threshold_H_p_cm3: 1e-4      # Physical density above which the EAGLE Jeans limiter entropy floor kicks in expressed in Hydrogen atoms per cm^3.
diff --git a/examples/EAGLE_ICs/EAGLE_50_low_res/eagle_50.yml b/examples/EAGLE_ICs/EAGLE_50_low_res/eagle_50.yml
index 6e97297e7b075ff8dc502546cdd9acbbb029b5e3..f91f7aabe9f92f0c5392897171ecd427ae41f7d4 100644
--- a/examples/EAGLE_ICs/EAGLE_50_low_res/eagle_50.yml
+++ b/examples/EAGLE_ICs/EAGLE_50_low_res/eagle_50.yml
@@ -134,20 +134,21 @@ COLIBRECooling:
 
 # EAGLE star formation parameters
 EAGLEStarFormation:
-  SF_model:                          PressureLaw
-  EOS_density_norm_H_p_cm3:          0.1       # Physical density used for the normalisation of the EOS assumed for the star-forming gas in Hydrogen atoms per cm^3.
-  EOS_temperature_norm_K:            8000      # Temperature om the polytropic EOS assumed for star-forming gas at the density normalisation in Kelvin.
-  EOS_gamma_effective:               1.3333333 # Slope the of the polytropic EOS assumed for the star-forming gas.
-  KS_normalisation:                  1.515e-4  # The normalization of the Kennicutt-Schmidt law in Msun / kpc^2 / yr.
-  KS_exponent:                       1.4       # The exponent of the Kennicutt-Schmidt law.
-  min_over_density:                  100.0     # The over-density above which star-formation is allowed.
-  KS_high_density_threshold_H_p_cm3: 1e3       # Hydrogen number density above which the Kennicut-Schmidt law changes slope in Hydrogen atoms per cm^3.
-  KS_high_density_exponent:          2.0       # Slope of the Kennicut-Schmidt law above the high-density threshold.
-  EOS_entropy_margin_dex:            0.3       # Logarithm base 10 of the maximal entropy above the EOS at which stars can form.
-  threshold_norm_H_p_cm3:            0.1       # Normalisation of the metal-dependant density threshold for star formation in Hydrogen atoms per cm^3.
-  threshold_Z0:                      0.002     # Reference metallicity (metal mass fraction) for the metal-dependant threshold for star formation.
-  threshold_slope:                   -0.64     # Slope of the metal-dependant star formation threshold
-  threshold_max_density_H_p_cm3:     10.0      # Maximal density of the metal-dependant density threshold for star formation in Hydrogen atoms per cm^3.
+  SF_threshold:                      Subgrid      # Zdep (Schaye 2004) or Subgrid
+  SF_model:                          PressureLaw  # PressureLaw (Schaye et al. 2008) or SchmidtLaw
+  KS_normalisation:                  1.515e-4     # The normalization of the Kennicutt-Schmidt law in Msun / kpc^2 / yr.
+  KS_exponent:                       1.4          # The exponent of the Kennicutt-Schmidt law.
+  min_over_density:                  100.0        # The over-density above which star-formation is allowed.
+  KS_high_density_threshold_H_p_cm3: 1e3          # Hydrogen number density above which the Kennicut-Schmidt law changes slope in Hydrogen atoms per cm^3.
+  KS_high_density_exponent:          2.0          # Slope of the Kennicut-Schmidt law above the high-density threshold.
+  EOS_entropy_margin_dex:            0.3          # When using Z-based SF threshold, logarithm base 10 of the maximal entropy above the EOS at which stars can form.
+  threshold_norm_H_p_cm3:            0.1          # When using Z-based SF threshold, normalisation of the metal-dependant density threshold for star formation in Hydrogen atoms per cm^3.
+  threshold_Z0:                      0.002        # When using Z-based SF threshold, reference metallicity (metal mass fraction) for the metal-dependant threshold for star formation.
+  threshold_slope:                   -0.64        # When using Z-based SF threshold, slope of the metal-dependant star formation threshold
+  threshold_max_density_H_p_cm3:     10.0         # When using Z-based SF threshold, maximal density of the metal-dependant density threshold for star formation in Hydrogen atoms per cm^3.
+  threshold_temperature1_K:          1000         # When using subgrid-based SF threshold, subgrid temperature below which gas is star-forming.
+  threshold_temperature2_K:          31622        # When using subgrid-based SF threshold, subgrid temperature below which gas is star-forming if also above the density limit.
+  threshold_number_density_H_p_cm3:  10           # When using subgrid-based SF threshold, subgrid number density above which gas is star-forming if also below the second temperature limit.
   
 # Parameters for the EAGLE "equation of state"
 EAGLEEntropyFloor:
diff --git a/examples/EAGLE_ICs/EAGLE_6/eagle_6.yml b/examples/EAGLE_ICs/EAGLE_6/eagle_6.yml
index 10076c0f376047c551420a12880e36ec3662261e..00b68c896d4cf503d2dc61b850fde06e6b87bdb7 100644
--- a/examples/EAGLE_ICs/EAGLE_6/eagle_6.yml
+++ b/examples/EAGLE_ICs/EAGLE_6/eagle_6.yml
@@ -137,20 +137,21 @@ COLIBRECooling:
 
 # EAGLE star formation parameters
 EAGLEStarFormation:
-  SF_model:                          PressureLaw
-  EOS_density_norm_H_p_cm3:          0.1       # Physical density used for the normalisation of the EOS assumed for the star-forming gas in Hydrogen atoms per cm^3.
-  EOS_temperature_norm_K:            8000      # Temperature om the polytropic EOS assumed for star-forming gas at the density normalisation in Kelvin.
-  EOS_gamma_effective:               1.3333333 # Slope the of the polytropic EOS assumed for the star-forming gas.
-  KS_normalisation:                  1.515e-4  # The normalization of the Kennicutt-Schmidt law in Msun / kpc^2 / yr.
-  KS_exponent:                       1.4       # The exponent of the Kennicutt-Schmidt law.
-  min_over_density:                  100.0     # The over-density above which star-formation is allowed.
-  KS_high_density_threshold_H_p_cm3: 1e3       # Hydrogen number density above which the Kennicut-Schmidt law changes slope in Hydrogen atoms per cm^3.
-  KS_high_density_exponent:          2.0       # Slope of the Kennicut-Schmidt law above the high-density threshold.
-  EOS_entropy_margin_dex:            0.3       # Logarithm base 10 of the maximal entropy above the EOS at which stars can form.
-  threshold_norm_H_p_cm3:            0.1       # Normalisation of the metal-dependant density threshold for star formation in Hydrogen atoms per cm^3.
-  threshold_Z0:                      0.002     # Reference metallicity (metal mass fraction) for the metal-dependant threshold for star formation.
-  threshold_slope:                   -0.64     # Slope of the metal-dependant star formation threshold
-  threshold_max_density_H_p_cm3:     10.0      # Maximal density of the metal-dependant density threshold for star formation in Hydrogen atoms per cm^3.
+  SF_threshold:                      Subgrid      # Zdep (Schaye 2004) or Subgrid
+  SF_model:                          PressureLaw  # PressureLaw (Schaye et al. 2008) or SchmidtLaw
+  KS_normalisation:                  1.515e-4     # The normalization of the Kennicutt-Schmidt law in Msun / kpc^2 / yr.
+  KS_exponent:                       1.4          # The exponent of the Kennicutt-Schmidt law.
+  min_over_density:                  100.0        # The over-density above which star-formation is allowed.
+  KS_high_density_threshold_H_p_cm3: 1e3          # Hydrogen number density above which the Kennicut-Schmidt law changes slope in Hydrogen atoms per cm^3.
+  KS_high_density_exponent:          2.0          # Slope of the Kennicut-Schmidt law above the high-density threshold.
+  EOS_entropy_margin_dex:            0.3          # When using Z-based SF threshold, logarithm base 10 of the maximal entropy above the EOS at which stars can form.
+  threshold_norm_H_p_cm3:            0.1          # When using Z-based SF threshold, normalisation of the metal-dependant density threshold for star formation in Hydrogen atoms per cm^3.
+  threshold_Z0:                      0.002        # When using Z-based SF threshold, reference metallicity (metal mass fraction) for the metal-dependant threshold for star formation.
+  threshold_slope:                   -0.64        # When using Z-based SF threshold, slope of the metal-dependant star formation threshold
+  threshold_max_density_H_p_cm3:     10.0         # When using Z-based SF threshold, maximal density of the metal-dependant density threshold for star formation in Hydrogen atoms per cm^3.
+  threshold_temperature1_K:          1000         # When using subgrid-based SF threshold, subgrid temperature below which gas is star-forming.
+  threshold_temperature2_K:          31622        # When using subgrid-based SF threshold, subgrid temperature below which gas is star-forming if also above the density limit.
+  threshold_number_density_H_p_cm3:  10           # When using subgrid-based SF threshold, subgrid number density above which gas is star-forming if also below the second temperature limit.
   
 # Parameters for the EAGLE "equation of state"
 EAGLEEntropyFloor:
diff --git a/examples/EAGLE_ICs/README b/examples/EAGLE_ICs/README
index 2ff5032ffce2882406fbb3238165516190b53761..77a900d7a8c1a09e9dda4509bc77cae146dfb6ad 100644
--- a/examples/EAGLE_ICs/README
+++ b/examples/EAGLE_ICs/README
@@ -35,9 +35,12 @@ the following changes have been made (at standard resolution):
    been removed as the new cooling tables handle this correctly. The
    gamma=4/3 floor has been extended to lower densities
    (i.e. 800K at n_H = 10^-4 cm^-3) as a fail-safe.
- - Particles can be star-forming if they are within 0.3 dex of the
-   entropy floor (was 0.5 dex). These particles also get their
-   subgrid properties (rho, T as well as the HI and H_2 frac) computed.
+ - Particles within 0.3 dex of the entropy floor get their
+   subgrid properties (rho, T as well as the HI and H_2 frac) computed
+   assuming pressure equilibrium.
+ - Particles are star-forming if their subgrid temperature is below
+   10^3 K OR if they are both colder than 10^4.5 K and have a subgrid
+   number density above 10 cm^-3.
  - Particles can be star-forming if they are in an over-density of
    at least 100 (was 57.7).
  - Particles with a density above 10^5 cm^-3 are not turned into
diff --git a/examples/EAGLE_low_z/EAGLE_12/eagle_12.yml b/examples/EAGLE_low_z/EAGLE_12/eagle_12.yml
index d5f6089503fe3df53fa46c1d34c3b16729d0cdf2..3a8e258e1e9cfa63c317133d36e47d24088eb68e 100644
--- a/examples/EAGLE_low_z/EAGLE_12/eagle_12.yml
+++ b/examples/EAGLE_low_z/EAGLE_12/eagle_12.yml
@@ -121,20 +121,21 @@ COLIBRECooling:
 
 # EAGLE star formation parameters
 EAGLEStarFormation:
-  SF_model:                          PressureLaw
-  EOS_density_norm_H_p_cm3:          0.1       # Physical density used for the normalisation of the EOS assumed for the star-forming gas in Hydrogen atoms per cm^3.
-  EOS_temperature_norm_K:            8000      # Temperature om the polytropic EOS assumed for star-forming gas at the density normalisation in Kelvin.
-  EOS_gamma_effective:               1.3333333 # Slope the of the polytropic EOS assumed for the star-forming gas.
-  KS_normalisation:                  1.515e-4  # The normalization of the Kennicutt-Schmidt law in Msun / kpc^2 / yr.
-  KS_exponent:                       1.4       # The exponent of the Kennicutt-Schmidt law.
-  min_over_density:                  100.0     # The over-density above which star-formation is allowed.
-  KS_high_density_threshold_H_p_cm3: 1e3       # Hydrogen number density above which the Kennicut-Schmidt law changes slope in Hydrogen atoms per cm^3.
-  KS_high_density_exponent:          2.0       # Slope of the Kennicut-Schmidt law above the high-density threshold.
-  EOS_entropy_margin_dex:            0.3       # Logarithm base 10 of the maximal entropy above the EOS at which stars can form.
-  threshold_norm_H_p_cm3:            0.1       # Normalisation of the metal-dependant density threshold for star formation in Hydrogen atoms per cm^3.
-  threshold_Z0:                      0.002     # Reference metallicity (metal mass fraction) for the metal-dependant threshold for star formation.
-  threshold_slope:                   -0.64     # Slope of the metal-dependant star formation threshold
-  threshold_max_density_H_p_cm3:     10.0      # Maximal density of the metal-dependant density threshold for star formation in Hydrogen atoms per cm^3.
+  SF_threshold:                      Zdep         # Zdep (Schaye 2004) or Subgrid
+  SF_model:                          PressureLaw  # PressureLaw (Schaye et al. 2008) or SchmidtLaw
+  KS_normalisation:                  1.515e-4     # The normalization of the Kennicutt-Schmidt law in Msun / kpc^2 / yr.
+  KS_exponent:                       1.4          # The exponent of the Kennicutt-Schmidt law.
+  min_over_density:                  100.0        # The over-density above which star-formation is allowed.
+  KS_high_density_threshold_H_p_cm3: 1e3          # Hydrogen number density above which the Kennicut-Schmidt law changes slope in Hydrogen atoms per cm^3.
+  KS_high_density_exponent:          2.0          # Slope of the Kennicut-Schmidt law above the high-density threshold.
+  EOS_entropy_margin_dex:            0.3          # When using Z-based SF threshold, logarithm base 10 of the maximal entropy above the EOS at which stars can form.
+  threshold_norm_H_p_cm3:            0.1          # When using Z-based SF threshold, normalisation of the metal-dependant density threshold for star formation in Hydrogen atoms per cm^3.
+  threshold_Z0:                      0.002        # When using Z-based SF threshold, reference metallicity (metal mass fraction) for the metal-dependant threshold for star formation.
+  threshold_slope:                   -0.64        # When using Z-based SF threshold, slope of the metal-dependant star formation threshold
+  threshold_max_density_H_p_cm3:     10.0         # When using Z-based SF threshold, maximal density of the metal-dependant density threshold for star formation in Hydrogen atoms per cm^3.
+  threshold_temperature1_K:          1000         # When using subgrid-based SF threshold, subgrid temperature below which gas is star-forming.
+  threshold_temperature2_K:          31622        # When using subgrid-based SF threshold, subgrid temperature below which gas is star-forming if also above the density limit.
+  threshold_number_density_H_p_cm3:  10           # When using subgrid-based SF threshold, subgrid number density above which gas is star-forming if also below the second temperature limit.
   
 # Parameters for the EAGLE "equation of state"
 EAGLEEntropyFloor:
diff --git a/examples/EAGLE_low_z/EAGLE_25/eagle_25.yml b/examples/EAGLE_low_z/EAGLE_25/eagle_25.yml
index df7992c172e0d74e99d798cd481c3c82ed8eba21..f03cb7960197339459405ab3ac1f21499c1fbf28 100644
--- a/examples/EAGLE_low_z/EAGLE_25/eagle_25.yml
+++ b/examples/EAGLE_low_z/EAGLE_25/eagle_25.yml
@@ -129,20 +129,21 @@ COLIBRECooling:
 
 # EAGLE star formation parameters
 EAGLEStarFormation:
-  SF_model:                          PressureLaw
-  EOS_density_norm_H_p_cm3:          0.1       # Physical density used for the normalisation of the EOS assumed for the star-forming gas in Hydrogen atoms per cm^3.
-  EOS_temperature_norm_K:            8000      # Temperature om the polytropic EOS assumed for star-forming gas at the density normalisation in Kelvin.
-  EOS_gamma_effective:               1.3333333 # Slope the of the polytropic EOS assumed for the star-forming gas.
-  KS_normalisation:                  1.515e-4  # The normalization of the Kennicutt-Schmidt law in Msun / kpc^2 / yr.
-  KS_exponent:                       1.4       # The exponent of the Kennicutt-Schmidt law.
-  min_over_density:                  100.0     # The over-density above which star-formation is allowed.
-  KS_high_density_threshold_H_p_cm3: 1e3       # Hydrogen number density above which the Kennicut-Schmidt law changes slope in Hydrogen atoms per cm^3.
-  KS_high_density_exponent:          2.0       # Slope of the Kennicut-Schmidt law above the high-density threshold.
-  EOS_entropy_margin_dex:            0.3       # Logarithm base 10 of the maximal entropy above the EOS at which stars can form.
-  threshold_norm_H_p_cm3:            0.1       # Normalisation of the metal-dependant density threshold for star formation in Hydrogen atoms per cm^3.
-  threshold_Z0:                      0.002     # Reference metallicity (metal mass fraction) for the metal-dependant threshold for star formation.
-  threshold_slope:                   -0.64     # Slope of the metal-dependant star formation threshold
-  threshold_max_density_H_p_cm3:     10.0      # Maximal density of the metal-dependant density threshold for star formation in Hydrogen atoms per cm^3.
+  SF_threshold:                      Zdep         # Zdep (Schaye 2004) or Subgrid
+  SF_model:                          PressureLaw  # PressureLaw (Schaye et al. 2008) or SchmidtLaw
+  KS_normalisation:                  1.515e-4     # The normalization of the Kennicutt-Schmidt law in Msun / kpc^2 / yr.
+  KS_exponent:                       1.4          # The exponent of the Kennicutt-Schmidt law.
+  min_over_density:                  100.0        # The over-density above which star-formation is allowed.
+  KS_high_density_threshold_H_p_cm3: 1e3          # Hydrogen number density above which the Kennicut-Schmidt law changes slope in Hydrogen atoms per cm^3.
+  KS_high_density_exponent:          2.0          # Slope of the Kennicut-Schmidt law above the high-density threshold.
+  EOS_entropy_margin_dex:            0.3          # When using Z-based SF threshold, logarithm base 10 of the maximal entropy above the EOS at which stars can form.
+  threshold_norm_H_p_cm3:            0.1          # When using Z-based SF threshold, normalisation of the metal-dependant density threshold for star formation in Hydrogen atoms per cm^3.
+  threshold_Z0:                      0.002        # When using Z-based SF threshold, reference metallicity (metal mass fraction) for the metal-dependant threshold for star formation.
+  threshold_slope:                   -0.64        # When using Z-based SF threshold, slope of the metal-dependant star formation threshold
+  threshold_max_density_H_p_cm3:     10.0         # When using Z-based SF threshold, maximal density of the metal-dependant density threshold for star formation in Hydrogen atoms per cm^3.
+  threshold_temperature1_K:          1000         # When using subgrid-based SF threshold, subgrid temperature below which gas is star-forming.
+  threshold_temperature2_K:          31622        # When using subgrid-based SF threshold, subgrid temperature below which gas is star-forming if also above the density limit.
+  threshold_number_density_H_p_cm3:  10           # When using subgrid-based SF threshold, subgrid number density above which gas is star-forming if also below the second temperature limit.
   
 # Parameters for the EAGLE "equation of state"
 EAGLEEntropyFloor:
diff --git a/examples/EAGLE_low_z/EAGLE_50/eagle_50.yml b/examples/EAGLE_low_z/EAGLE_50/eagle_50.yml
index e8d1a4bd8d9e26d4d3ec6d7f3e6564764875d615..3aac5d4f188405712b733e83b34c595c1b3f97fe 100644
--- a/examples/EAGLE_low_z/EAGLE_50/eagle_50.yml
+++ b/examples/EAGLE_low_z/EAGLE_50/eagle_50.yml
@@ -120,20 +120,21 @@ COLIBRECooling:
 
 # EAGLE star formation parameters
 EAGLEStarFormation:
-  SF_model:                          PressureLaw
-  EOS_density_norm_H_p_cm3:          0.1       # Physical density used for the normalisation of the EOS assumed for the star-forming gas in Hydrogen atoms per cm^3.
-  EOS_temperature_norm_K:            8000      # Temperature om the polytropic EOS assumed for star-forming gas at the density normalisation in Kelvin.
-  EOS_gamma_effective:               1.3333333 # Slope the of the polytropic EOS assumed for the star-forming gas.
-  KS_normalisation:                  1.515e-4  # The normalization of the Kennicutt-Schmidt law in Msun / kpc^2 / yr.
-  KS_exponent:                       1.4       # The exponent of the Kennicutt-Schmidt law.
-  min_over_density:                  100.0     # The over-density above which star-formation is allowed.
-  KS_high_density_threshold_H_p_cm3: 1e3       # Hydrogen number density above which the Kennicut-Schmidt law changes slope in Hydrogen atoms per cm^3.
-  KS_high_density_exponent:          2.0       # Slope of the Kennicut-Schmidt law above the high-density threshold.
-  EOS_entropy_margin_dex:            0.3       # Logarithm base 10 of the maximal entropy above the EOS at which stars can form.
-  threshold_norm_H_p_cm3:            0.1       # Normalisation of the metal-dependant density threshold for star formation in Hydrogen atoms per cm^3.
-  threshold_Z0:                      0.002     # Reference metallicity (metal mass fraction) for the metal-dependant threshold for star formation.
-  threshold_slope:                   -0.64     # Slope of the metal-dependant star formation threshold
-  threshold_max_density_H_p_cm3:     10.0      # Maximal density of the metal-dependant density threshold for star formation in Hydrogen atoms per cm^3.
+  SF_threshold:                      Zdep         # Zdep (Schaye 2004) or Subgrid
+  SF_model:                          PressureLaw  # PressureLaw (Schaye et al. 2008) or SchmidtLaw
+  KS_normalisation:                  1.515e-4     # The normalization of the Kennicutt-Schmidt law in Msun / kpc^2 / yr.
+  KS_exponent:                       1.4          # The exponent of the Kennicutt-Schmidt law.
+  min_over_density:                  100.0        # The over-density above which star-formation is allowed.
+  KS_high_density_threshold_H_p_cm3: 1e3          # Hydrogen number density above which the Kennicut-Schmidt law changes slope in Hydrogen atoms per cm^3.
+  KS_high_density_exponent:          2.0          # Slope of the Kennicut-Schmidt law above the high-density threshold.
+  EOS_entropy_margin_dex:            0.3          # When using Z-based SF threshold, logarithm base 10 of the maximal entropy above the EOS at which stars can form.
+  threshold_norm_H_p_cm3:            0.1          # When using Z-based SF threshold, normalisation of the metal-dependant density threshold for star formation in Hydrogen atoms per cm^3.
+  threshold_Z0:                      0.002        # When using Z-based SF threshold, reference metallicity (metal mass fraction) for the metal-dependant threshold for star formation.
+  threshold_slope:                   -0.64        # When using Z-based SF threshold, slope of the metal-dependant star formation threshold
+  threshold_max_density_H_p_cm3:     10.0         # When using Z-based SF threshold, maximal density of the metal-dependant density threshold for star formation in Hydrogen atoms per cm^3.
+  threshold_temperature1_K:          1000         # When using subgrid-based SF threshold, subgrid temperature below which gas is star-forming.
+  threshold_temperature2_K:          31622        # When using subgrid-based SF threshold, subgrid temperature below which gas is star-forming if also above the density limit.
+  threshold_number_density_H_p_cm3:  10           # When using subgrid-based SF threshold, subgrid number density above which gas is star-forming if also below the second temperature limit.
   
 # Parameters for the EAGLE "equation of state"
 EAGLEEntropyFloor:
diff --git a/examples/EAGLE_low_z/EAGLE_6/eagle_6.yml b/examples/EAGLE_low_z/EAGLE_6/eagle_6.yml
index 68682d92a7be4eb801f48819d0bd5644239dd611..b678dd11a18b340d1945a9bdf28af3b1816f7f1f 100644
--- a/examples/EAGLE_low_z/EAGLE_6/eagle_6.yml
+++ b/examples/EAGLE_low_z/EAGLE_6/eagle_6.yml
@@ -132,20 +132,21 @@ COLIBRECooling:
   
 # EAGLE star formation parameters
 EAGLEStarFormation:
-  SF_model:                          PressureLaw
-  EOS_density_norm_H_p_cm3:          0.1       # Physical density used for the normalisation of the EOS assumed for the star-forming gas in Hydrogen atoms per cm^3.
-  EOS_temperature_norm_K:            8000      # Temperature om the polytropic EOS assumed for star-forming gas at the density normalisation in Kelvin.
-  EOS_gamma_effective:               1.3333333 # Slope the of the polytropic EOS assumed for the star-forming gas.
-  KS_normalisation:                  1.515e-4  # The normalization of the Kennicutt-Schmidt law in Msun / kpc^2 / yr.
-  KS_exponent:                       1.4       # The exponent of the Kennicutt-Schmidt law.
-  min_over_density:                  100.0     # The over-density above which star-formation is allowed.
-  KS_high_density_threshold_H_p_cm3: 1e3       # Hydrogen number density above which the Kennicut-Schmidt law changes slope in Hydrogen atoms per cm^3.
-  KS_high_density_exponent:          2.0       # Slope of the Kennicut-Schmidt law above the high-density threshold.
-  EOS_entropy_margin_dex:            0.3       # Logarithm base 10 of the maximal entropy above the EOS at which stars can form.
-  threshold_norm_H_p_cm3:            0.1       # Normalisation of the metal-dependant density threshold for star formation in Hydrogen atoms per cm^3.
-  threshold_Z0:                      0.002     # Reference metallicity (metal mass fraction) for the metal-dependant threshold for star formation.
-  threshold_slope:                   -0.64     # Slope of the metal-dependant star formation threshold
-  threshold_max_density_H_p_cm3:     10.0      # Maximal density of the metal-dependant density threshold for star formation in Hydrogen atoms per cm^3.
+  SF_threshold:                      Zdep         # Zdep (Schaye 2004) or Subgrid
+  SF_model:                          PressureLaw  # PressureLaw (Schaye et al. 2008) or SchmidtLaw
+  KS_normalisation:                  1.515e-4     # The normalization of the Kennicutt-Schmidt law in Msun / kpc^2 / yr.
+  KS_exponent:                       1.4          # The exponent of the Kennicutt-Schmidt law.
+  min_over_density:                  100.0        # The over-density above which star-formation is allowed.
+  KS_high_density_threshold_H_p_cm3: 1e3          # Hydrogen number density above which the Kennicut-Schmidt law changes slope in Hydrogen atoms per cm^3.
+  KS_high_density_exponent:          2.0          # Slope of the Kennicut-Schmidt law above the high-density threshold.
+  EOS_entropy_margin_dex:            0.3          # When using Z-based SF threshold, logarithm base 10 of the maximal entropy above the EOS at which stars can form.
+  threshold_norm_H_p_cm3:            0.1          # When using Z-based SF threshold, normalisation of the metal-dependant density threshold for star formation in Hydrogen atoms per cm^3.
+  threshold_Z0:                      0.002        # When using Z-based SF threshold, reference metallicity (metal mass fraction) for the metal-dependant threshold for star formation.
+  threshold_slope:                   -0.64        # When using Z-based SF threshold, slope of the metal-dependant star formation threshold
+  threshold_max_density_H_p_cm3:     10.0         # When using Z-based SF threshold, maximal density of the metal-dependant density threshold for star formation in Hydrogen atoms per cm^3.
+  threshold_temperature1_K:          1000         # When using subgrid-based SF threshold, subgrid temperature below which gas is star-forming.
+  threshold_temperature2_K:          31622        # When using subgrid-based SF threshold, subgrid temperature below which gas is star-forming if also above the density limit.
+  threshold_number_density_H_p_cm3:  10           # When using subgrid-based SF threshold, subgrid number density above which gas is star-forming if also below the second temperature limit.
   
 # Parameters for the EAGLE "equation of state"
 EAGLEEntropyFloor:
diff --git a/examples/IsolatedGalaxy/IsolatedGalaxy_feedback/isolated_galaxy.yml b/examples/IsolatedGalaxy/IsolatedGalaxy_feedback/isolated_galaxy.yml
index 91df03828740b7e2cfc1d4e368ce51730a228a58..a0e0317e87f532b1169d13322735de9fb78e5356 100644
--- a/examples/IsolatedGalaxy/IsolatedGalaxy_feedback/isolated_galaxy.yml
+++ b/examples/IsolatedGalaxy/IsolatedGalaxy_feedback/isolated_galaxy.yml
@@ -111,20 +111,21 @@ HernquistPotential:
  
 # EAGLE star formation parameters
 EAGLEStarFormation:
-  SF_model:                          PressureLaw
-  EOS_density_norm_H_p_cm3:          0.1       # Physical density used for the normalisation of the EOS assumed for the star-forming gas in Hydrogen atoms per cm^3.
-  EOS_temperature_norm_K:            8000      # Temperature om the polytropic EOS assumed for star-forming gas at the density normalisation in Kelvin.
-  EOS_gamma_effective:               1.3333333 # Slope the of the polytropic EOS assumed for the star-forming gas.
-  KS_normalisation:                  1.515e-4  # The normalization of the Kennicutt-Schmidt law in Msun / kpc^2 / yr.
-  KS_exponent:                       1.4       # The exponent of the Kennicutt-Schmidt law.
-  min_over_density:                  100.0     # The over-density above which star-formation is allowed.
-  KS_high_density_threshold_H_p_cm3: 1e3       # Hydrogen number density above which the Kennicut-Schmidt law changes slope in Hydrogen atoms per cm^3.
-  KS_high_density_exponent:          2.0       # Slope of the Kennicut-Schmidt law above the high-density threshold.
-  EOS_entropy_margin_dex:            0.3       # Logarithm base 10 of the maximal entropy above the EOS at which stars can form.
-  threshold_norm_H_p_cm3:            0.1       # Normalisation of the metal-dependant density threshold for star formation in Hydrogen atoms per cm^3.
-  threshold_Z0:                      0.002     # Reference metallicity (metal mass fraction) for the metal-dependant threshold for star formation.
-  threshold_slope:                   -0.64     # Slope of the metal-dependant star formation threshold
-  threshold_max_density_H_p_cm3:     10.0      # Maximal density of the metal-dependant density threshold for star formation in Hydrogen atoms per cm^3.
+  SF_threshold:                      Zdep         # Zdep (Schaye 2004) or Subgrid
+  SF_model:                          PressureLaw  # PressureLaw (Schaye et al. 2008) or SchmidtLaw
+  KS_normalisation:                  1.515e-4     # The normalization of the Kennicutt-Schmidt law in Msun / kpc^2 / yr.
+  KS_exponent:                       1.4          # The exponent of the Kennicutt-Schmidt law.
+  min_over_density:                  100.0        # The over-density above which star-formation is allowed.
+  KS_high_density_threshold_H_p_cm3: 1e3          # Hydrogen number density above which the Kennicut-Schmidt law changes slope in Hydrogen atoms per cm^3.
+  KS_high_density_exponent:          2.0          # Slope of the Kennicut-Schmidt law above the high-density threshold.
+  EOS_entropy_margin_dex:            0.3          # When using Z-based SF threshold, logarithm base 10 of the maximal entropy above the EOS at which stars can form.
+  threshold_norm_H_p_cm3:            0.1          # When using Z-based SF threshold, normalisation of the metal-dependant density threshold for star formation in Hydrogen atoms per cm^3.
+  threshold_Z0:                      0.002        # When using Z-based SF threshold, reference metallicity (metal mass fraction) for the metal-dependant threshold for star formation.
+  threshold_slope:                   -0.64        # When using Z-based SF threshold, slope of the metal-dependant star formation threshold
+  threshold_max_density_H_p_cm3:     10.0         # When using Z-based SF threshold, maximal density of the metal-dependant density threshold for star formation in Hydrogen atoms per cm^3.
+  threshold_temperature1_K:          1000         # When using subgrid-based SF threshold, subgrid temperature below which gas is star-forming.
+  threshold_temperature2_K:          31622        # When using subgrid-based SF threshold, subgrid temperature below which gas is star-forming if also above the density limit.
+  threshold_number_density_H_p_cm3:  10           # When using subgrid-based SF threshold, subgrid number density above which gas is star-forming if also below the second temperature limit.
   
 # Parameters for the EAGLE "equation of state"
 EAGLEEntropyFloor:
diff --git a/examples/IsolatedGalaxy/IsolatedGalaxy_starformation/isolated_galaxy.yml b/examples/IsolatedGalaxy/IsolatedGalaxy_starformation/isolated_galaxy.yml
index 983d185c9089fae93253d2859dba3050497d2d51..af334ff7d4f2834c18c3af97d8b01259e559ebed 100644
--- a/examples/IsolatedGalaxy/IsolatedGalaxy_starformation/isolated_galaxy.yml
+++ b/examples/IsolatedGalaxy/IsolatedGalaxy_starformation/isolated_galaxy.yml
@@ -103,21 +103,22 @@ HernquistPotential:
  
 # EAGLE star formation parameters
 EAGLEStarFormation:
-  SF_model:                          PressureLaw
-  EOS_density_norm_H_p_cm3:          0.1       # Physical density used for the normalisation of the EOS assumed for the star-forming gas in Hydrogen atoms per cm^3.
-  EOS_temperature_norm_K:            8000      # Temperature om the polytropic EOS assumed for star-forming gas at the density normalisation in Kelvin.
-  EOS_gamma_effective:               1.3333333 # Slope the of the polytropic EOS assumed for the star-forming gas.
-  gas_fraction:                      0.3       # The gas fraction used internally by the model.
-  KS_normalisation:                  1.515e-4  # The normalization of the Kennicutt-Schmidt law in Msun / kpc^2 / yr.
-  KS_exponent:                       1.4       # The exponent of the Kennicutt-Schmidt law.
-  min_over_density:                  100.0     # The over-density above which star-formation is allowed.
-  KS_high_density_threshold_H_p_cm3: 1e3       # Hydrogen number density above which the Kennicut-Schmidt law changes slope in Hydrogen atoms per cm^3.
-  KS_high_density_exponent:          2.0       # Slope of the Kennicut-Schmidt law above the high-density threshold.
-  EOS_entropy_margin_dex:            0.3       # Logarithm base 10 of the maximal entropy above the EOS at which stars can form.
-  threshold_norm_H_p_cm3:            0.1       # Normalisation of the metal-dependant density threshold for star formation in Hydrogen atoms per cm^3.
-  threshold_Z0:                      0.002     # Reference metallicity (metal mass fraction) for the metal-dependant threshold for star formation.
-  threshold_slope:                   -0.64     # Slope of the metal-dependant star formation threshold
-  threshold_max_density_H_p_cm3:     10.0      # Maximal density of the metal-dependant density threshold for star formation in Hydrogen atoms per cm^3.
+  SF_threshold:                      Zdep         # Zdep (Schaye 2004) or Subgrid
+  SF_model:                          PressureLaw  # PressureLaw (Schaye et al. 2008) or SchmidtLaw
+  gas_fraction:                      0.3          # The gas fraction used internally by the model.
+  KS_normalisation:                  1.515e-4     # The normalization of the Kennicutt-Schmidt law in Msun / kpc^2 / yr.
+  KS_exponent:                       1.4          # The exponent of the Kennicutt-Schmidt law.
+  min_over_density:                  100.0        # The over-density above which star-formation is allowed.
+  KS_high_density_threshold_H_p_cm3: 1e3          # Hydrogen number density above which the Kennicut-Schmidt law changes slope in Hydrogen atoms per cm^3.
+  KS_high_density_exponent:          2.0          # Slope of the Kennicut-Schmidt law above the high-density threshold.
+  EOS_entropy_margin_dex:            0.3          # When using Z-based SF threshold, logarithm base 10 of the maximal entropy above the EOS at which stars can form.
+  threshold_norm_H_p_cm3:            0.1          # When using Z-based SF threshold, normalisation of the metal-dependant density threshold for star formation in Hydrogen atoms per cm^3.
+  threshold_Z0:                      0.002        # When using Z-based SF threshold, reference metallicity (metal mass fraction) for the metal-dependant threshold for star formation.
+  threshold_slope:                   -0.64        # When using Z-based SF threshold, slope of the metal-dependant star formation threshold
+  threshold_max_density_H_p_cm3:     10.0         # When using Z-based SF threshold, maximal density of the metal-dependant density threshold for star formation in Hydrogen atoms per cm^3.
+  threshold_temperature1_K:          1000         # When using subgrid-based SF threshold, subgrid temperature below which gas is star-forming.
+  threshold_temperature2_K:          31622        # When using subgrid-based SF threshold, subgrid temperature below which gas is star-forming if also above the density limit.
+  threshold_number_density_H_p_cm3:  10           # When using subgrid-based SF threshold, subgrid number density above which gas is star-forming if also below the second temperature limit.
   
 # Parameters for the EAGLE "equation of state"
 EAGLEEntropyFloor:
diff --git a/examples/main.c b/examples/main.c
index edb6c6461f13c6171ff389b5f859df0a4157bb4d..0510211985882f05e118dd7681ea89b21494c733 100644
--- a/examples/main.c
+++ b/examples/main.c
@@ -1043,8 +1043,8 @@ int main(int argc, char *argv[]) {
 #ifdef STAR_FORMATION_NONE
       error("ERROR: Running with star formation but compiled without it!");
 #endif
-      starformation_init(params, &prog_const, &us, &hydro_properties,
-                         &starform);
+      starformation_init(params, &prog_const, &us, &hydro_properties, &cosmo,
+                         &entropy_floor, &starform);
     }
     if (with_star_formation && myrank == 0) starformation_print(&starform);
 
diff --git a/examples/parameter_example.yml b/examples/parameter_example.yml
index a1a6e2c44a57f094049e977d55528e6bd0d590ab..6a881d5a4409ffe64542c540a90116f8b62c7390 100644
--- a/examples/parameter_example.yml
+++ b/examples/parameter_example.yml
@@ -458,7 +458,8 @@ GEARChemistry:
 
 # EAGLE star formation model (Schaye and Dalla Vecchia 2008)
 EAGLEStarFormation:
-  SF_model:                     SchmidtLaw     # "SchmidtLaw" or "PressureLaw".
+  SF_threshold:                   Subgrid      # Zdep (Schaye 2004) or Subgrid
+  SF_model:                       PressureLaw  # PressureLaw (Schaye et al. 2008) or SchmidtLaw
   star_formation_efficiency:         0.01      # In the SchmidtLaw model this regulates the star formation efficiency per free-fall time.
   KS_normalisation:                  1.515e-4  # In the PressureLaw model, normalization of the Kennicutt-Schmidt law in Msun / kpc^2 / yr.
   KS_exponent:                       1.4       # In the PressureLaw model, exponent of the Kennicutt-Schmidt law.
@@ -466,15 +467,14 @@ EAGLEStarFormation:
   KS_high_density_exponent:          2.0       # In the PressureLaw model, Slope of the Kennicut-Schmidt law above the high-density threshold.
   gas_fraction:                      0.25      # (Optional) In the PressureLaw model, The gas fraction used internally by the model (Defaults to 1).
   density_direct_H_p_cm3:            1e5       # (Optional) Density above which a gas particle gets automatically turned into a star in Hydrogen atoms per cm^3 (Defaults to FLT_MAX).
-  EOS_density_norm_H_p_cm3:          0.1       # Physical density used for the normalisation of the EoS assumed for the star-forming gas in Hydrogen atoms per cm^3.
-  EOS_temperature_norm_K:            8000      # Temperature om the polytropic EoS assumed for star-forming gas at the density normalisation in Kelvin.
-  EOS_gamma_effective:               1.3333333 # Slope the of the polytropic EoS assumed for the star-forming gas.
-  EOS_entropy_margin_dex:            0.5       # (Optional) Logarithm base 10 of the maximal entropy above the EoS at which stars can form.
-  min_over_density:                  57.7      # Over-density above which star-formation is allowed.
-  threshold_norm_H_p_cm3:            0.1       # Normalisation of the metal-dependant density threshold for star formation in Hydrogen atoms per cm^3.
-  threshold_Z0:                      0.002     # Reference metallicity (metal mass fraction) for the metal-dependant threshold for star formation.
-  threshold_slope:                   -0.64     # Slope of the metal-dependant star formation threshold
-  threshold_max_density_H_p_cm3:     10.0      # Maximal density of the metal-dependant density threshold for star formation in Hydrogen atoms per cm^3.  
+  EOS_entropy_margin_dex:            0.3       # When using Z-based SF threshold, logarithm base 10 of the maximal entropy above the EOS at which stars can form.
+  threshold_norm_H_p_cm3:            0.1       # When using Z-based SF threshold, normalisation of the metal-dependant density threshold for star formation in Hydrogen atoms per cm^3.
+  threshold_Z0:                      0.002     # When using Z-based SF threshold, reference metallicity (metal mass fraction) for the metal-dependant threshold for star formation.
+  threshold_slope:                   -0.64     # When using Z-based SF threshold, slope of the metal-dependant star formation threshold
+  threshold_max_density_H_p_cm3:     10.0      # When using Z-based SF threshold, maximal density of the metal-dependant density threshold for star formation in Hydrogen atoms per cm^3.
+  threshold_temperature1_K:          1000      # When using subgrid-based SF threshold, subgrid temperature below which gas is star-forming.
+  threshold_temperature2_K:          31622     # When using subgrid-based SF threshold, subgrid temperature below which gas is star-forming if also above the density limit.
+  threshold_number_density_H_p_cm3:  10        # When using subgrid-based SF threshold, subgrid number density above which gas is star-forming if also below the second temperature limit.
 
 # Quick Lyman-alpha star formation parameters
 QLAStarFormation:
diff --git a/src/chemistry/EAGLE/chemistry.h b/src/chemistry/EAGLE/chemistry.h
index 4d7c37800d7d578c0ce7501e9c8f9a5ae447ff21..cab8916bce17fb5beda5f01ef679bf2fdbdaf0ad 100644
--- a/src/chemistry/EAGLE/chemistry.h
+++ b/src/chemistry/EAGLE/chemistry.h
@@ -603,7 +603,7 @@ chemistry_get_total_metal_mass_for_stats(const struct part* restrict p) {
  * @brief Returns the total metallicity (metal mass fraction) of the
  * star particle to be used in the stats related routines.
  *
- * @param p Pointer to the particle data.
+ * @param sp Pointer to the star particle data.
  */
 __attribute__((always_inline)) INLINE static float
 chemistry_get_star_total_metal_mass_for_stats(const struct spart* restrict sp) {
@@ -615,7 +615,7 @@ chemistry_get_star_total_metal_mass_for_stats(const struct spart* restrict sp) {
  * @brief Returns the total metallicity (metal mass fraction) of the
  * black hole particle to be used in the stats related routines.
  *
- * @param p Pointer to the particle data.
+ * @param bp Pointer to the BH particle data.
  */
 __attribute__((always_inline)) INLINE static float
 chemistry_get_bh_total_metal_mass_for_stats(const struct bpart* restrict bp) {
diff --git a/src/chemistry/GEAR/chemistry.h b/src/chemistry/GEAR/chemistry.h
index ef3f0e999c59efbffaa9f03867ef461b189a6943..754d1636ddf9ecf30931e3ecf2db7ad3025822b0 100644
--- a/src/chemistry/GEAR/chemistry.h
+++ b/src/chemistry/GEAR/chemistry.h
@@ -415,7 +415,7 @@ chemistry_get_metal_mass_fraction_for_feedback(const struct part* restrict p) {
  *
  * This is unused in GEAR. --> return 0
  *
- * @param sp Pointer to the particle data.
+ * @param p Pointer to the particle data.
  */
 __attribute__((always_inline)) INLINE static float
 chemistry_get_total_metal_mass_fraction_for_feedback(
@@ -519,7 +519,7 @@ chemistry_get_total_metal_mass_for_stats(const struct part* restrict p) {
  * @brief Returns the total metallicity (metal mass fraction) of the
  * star particle to be used in the stats related routines.
  *
- * @param p Pointer to the particle data.
+ * @param sp Pointer to the star particle data.
  */
 __attribute__((always_inline)) INLINE static float
 chemistry_get_star_total_metal_mass_for_stats(const struct spart* restrict sp) {
@@ -533,7 +533,7 @@ chemistry_get_star_total_metal_mass_for_stats(const struct spart* restrict sp) {
  * @brief Returns the total metallicity (metal mass fraction) of the
  * black hole particle to be used in the stats related routines.
  *
- * @param p Pointer to the particle data.
+ * @param bp Pointer to the BH particle data.
  */
 __attribute__((always_inline)) INLINE static float
 chemistry_get_bh_total_metal_mass_for_stats(const struct bpart* restrict bp) {
diff --git a/src/cooling/COLIBRE/cooling.c b/src/cooling/COLIBRE/cooling.c
index b4bb935d62cd6c972c90c6c0d4d93e29244fbb2b..0109baf3ea51e08a79a1a29d9ab0709b8f57a074 100644
--- a/src/cooling/COLIBRE/cooling.c
+++ b/src/cooling/COLIBRE/cooling.c
@@ -1072,6 +1072,30 @@ void cooling_set_particle_subgrid_properties(
       cooling_compute_subgrid_density);
 }
 
+/**
+ * @brief Returns the subgrid temperature of a particle.
+ *
+ * @param p The particle.
+ * @param xp The extended particle data.
+ * @return The subgrid temperature in internal units.
+ */
+float cooling_get_subgrid_temperature(const struct part *p,
+                                      const struct xpart *xp) {
+  return p->cooling_data.subgrid_temp;
+}
+
+/**
+ * @brief Returns the subgrid density of a particle.
+ *
+ * @param p The particle.
+ * @param xp The extended particle data.
+ * @return The subgrid density in physical internal units.
+ */
+float cooling_get_subgrid_density(const struct part *p,
+                                  const struct xpart *xp) {
+  return p->cooling_data.subgrid_dens;
+}
+
 /**
  * @brief Returns the total radiated energy by this particle.
  *
diff --git a/src/cooling/COLIBRE/cooling.h b/src/cooling/COLIBRE/cooling.h
index 94d3e2bdb225737f392e1f154fc7c207581d650a..52208601fdb7c855ee3ad2b9c345970b80fa3de1 100644
--- a/src/cooling/COLIBRE/cooling.h
+++ b/src/cooling/COLIBRE/cooling.h
@@ -134,6 +134,11 @@ double compute_subgrid_property(
     const float abundance_ratio[colibre_cooling_N_elementtypes],
     const double log_u_cgs, const enum cooling_subgrid_properties isub);
 
+float cooling_get_subgrid_temperature(const struct part *p,
+                                      const struct xpart *xp);
+
+float cooling_get_subgrid_density(const struct part *p, const struct xpart *xp);
+
 float cooling_get_radiated_energy(const struct xpart *xp);
 
 void cooling_split_part(struct part *p, struct xpart *xp, double n);
diff --git a/src/cooling/Compton/cooling.h b/src/cooling/Compton/cooling.h
index ca8b428d94324703e3c253ff92faa61fdb3082b5..7b2530c04590e9257f96e8b15c78e0c2915c72a0 100644
--- a/src/cooling/Compton/cooling.h
+++ b/src/cooling/Compton/cooling.h
@@ -302,6 +302,34 @@ INLINE static float cooling_get_temperature(
     return T_transition;
 }
 
+/**
+ * @brief Returns the subgrid temperature of a particle.
+ *
+ * This model has no subgrid quantity. We return an error.
+ *
+ * @param p The particle.
+ * @param xp The extended particle data.
+ */
+INLINE static float cooling_get_subgrid_temperature(const struct part* p,
+                                                    const struct xpart* xp) {
+  error("This cooling model does not use subgrid quantities!");
+  return -1.f;
+}
+
+/**
+ * @brief Returns the subgrid density of a particle.
+ *
+ * This model has no subgrid quantity. We return an error.
+ *
+ * @param p The particle.
+ * @param xp The extended particle data.
+ */
+INLINE static float cooling_get_subgrid_density(const struct part* p,
+                                                const struct xpart* xp) {
+  error("This cooling model does not use subgrid quantities!");
+  return -1.f;
+}
+
 /**
  * @brief Returns the total radiated energy by this particle.
  *
diff --git a/src/cooling/EAGLE/cooling.c b/src/cooling/EAGLE/cooling.c
index 22ff25442496ee4e74b72fd5c46e078def370190..481eca2842ad3056551efe632fe0488b743a1152 100644
--- a/src/cooling/EAGLE/cooling.c
+++ b/src/cooling/EAGLE/cooling.c
@@ -709,6 +709,34 @@ double compute_subgrid_density(
   return rho_phys;
 }
 
+/**
+ * @brief Returns the subgrid temperature of a particle.
+ *
+ * This model has no subgrid quantity. We return an error.
+ *
+ * @param p The particle.
+ * @param xp The extended particle data.
+ */
+float cooling_get_subgrid_temperature(const struct part *p,
+                                      const struct xpart *xp) {
+  error("This cooling model does not use subgrid quantities!");
+  return -1.f;
+}
+
+/**
+ * @brief Returns the subgrid density of a particle.
+ *
+ * This model has no subgrid quantity. We return an error.
+ *
+ * @param p The particle.
+ * @param xp The extended particle data.
+ */
+float cooling_get_subgrid_density(const struct part *p,
+                                  const struct xpart *xp) {
+  error("This cooling model does not use subgrid quantities!");
+  return -1.f;
+}
+
 /**
  * @brief Returns the total radiated energy by this particle.
  *
diff --git a/src/cooling/EAGLE/cooling.h b/src/cooling/EAGLE/cooling.h
index a851019aeaa23520f42ca060bfedb82840d88e55..960d22d25ba07ddb390804aab1103f9b05cdb238 100644
--- a/src/cooling/EAGLE/cooling.h
+++ b/src/cooling/EAGLE/cooling.h
@@ -26,6 +26,7 @@
 
 /* Local includes. */
 #include "cooling_properties.h"
+#include "error.h"
 
 struct part;
 struct xpart;
@@ -87,6 +88,11 @@ double compute_subgrid_property(
     const float *abundance_ratio, const double log_u_cgs,
     const enum cooling_subgrid_properties isub);
 
+float cooling_get_subgrid_temperature(const struct part *p,
+                                      const struct xpart *xp);
+
+float cooling_get_subgrid_density(const struct part *p, const struct xpart *xp);
+
 float cooling_get_radiated_energy(const struct xpart *restrict xp);
 
 void cooling_split_part(struct part *p, struct xpart *xp, double n);
diff --git a/src/cooling/QLA/cooling.c b/src/cooling/QLA/cooling.c
index 7f321b4d270f8b8395d8daa9ac35e5d7e18fe5de..ebf9d32135e09c006a4db9a9d6bdf290f5b9483d 100644
--- a/src/cooling/QLA/cooling.c
+++ b/src/cooling/QLA/cooling.c
@@ -685,6 +685,34 @@ __attribute__((always_inline)) INLINE void cooling_first_init_part(
     const struct cooling_function_data *cooling, struct part *p,
     struct xpart *xp) {}
 
+/**
+ * @param Returns the subgrid temperature of a particle.
+ *
+ * This model has no subgrid quantity. We return an error.
+ *
+ * @param p The particle.
+ * @param xp The extended particle data.
+ */
+float cooling_get_subgrid_temperature(const struct part *p,
+                                      const struct xpart *xp) {
+  error("This cooling model does not use subgrid quantities!");
+  return -1.f;
+}
+
+/**
+ * @param Returns the subgrid density of a particle.
+ *
+ * This model has no subgrid quantity. We return an error.
+ *
+ * @param p The particle.
+ * @param xp The extended particle data.
+ */
+float cooling_get_subgrid_density(const struct part *p,
+                                  const struct xpart *xp) {
+  error("This cooling model does not use subgrid quantities!");
+  return -1.f;
+}
+
 /**
  * @brief Returns the total radiated energy by this particle.
  *
diff --git a/src/cooling/QLA/cooling.h b/src/cooling/QLA/cooling.h
index b11acc51f8530d1e03a6bbb4f245ba0af1e68c5b..f62756064ed98214fb586d51c6fd94253f905541 100644
--- a/src/cooling/QLA/cooling.h
+++ b/src/cooling/QLA/cooling.h
@@ -76,6 +76,11 @@ float cooling_get_temperature(const struct phys_const *phys_const,
                               const struct cooling_function_data *cooling,
                               const struct part *p, const struct xpart *xp);
 
+float cooling_get_subgrid_temperature(const struct part *p,
+                                      const struct xpart *xp);
+
+float cooling_get_subgrid_density(const struct part *p, const struct xpart *xp);
+
 float cooling_get_radiated_energy(const struct xpart *xp);
 
 void cooling_split_part(struct part *p, struct xpart *xp, double n);
diff --git a/src/cooling/const_du/cooling.h b/src/cooling/const_du/cooling.h
index 399db01cf290c3511e0fac63079c401622b1152f..c5b4f4ee28cd12c1267b4298c670e4ae0ad8bef6 100644
--- a/src/cooling/const_du/cooling.h
+++ b/src/cooling/const_du/cooling.h
@@ -218,6 +218,34 @@ INLINE static float cooling_get_temperature(
     return T_transition;
 }
 
+/**
+ * @brief Returns the subgrid temperature of a particle.
+ *
+ * This model has no subgrid quantity. We return an error.
+ *
+ * @param p The particle.
+ * @param xp The extended particle data.
+ */
+INLINE static float cooling_get_subgrid_temperature(const struct part* p,
+                                                    const struct xpart* xp) {
+  error("This cooling model does not use subgrid quantities!");
+  return -1.f;
+}
+
+/**
+ * @brief Returns the subgrid density of a particle.
+ *
+ * This model has no subgrid quantity. We return an error.
+ *
+ * @param p The particle.
+ * @param xp The extended particle data.
+ */
+INLINE static float cooling_get_subgrid_density(const struct part* p,
+                                                const struct xpart* xp) {
+  error("This cooling model does not use subgrid quantities!");
+  return -1.f;
+}
+
 /**
  * @brief Returns the total radiated energy by this particle.
  *
diff --git a/src/cooling/const_lambda/cooling.h b/src/cooling/const_lambda/cooling.h
index 92d11067c5dac5ca2b40b9e68eadc18fed0dae3c..7659177ef797bc2ff3294c8c3af338e647daaea8 100644
--- a/src/cooling/const_lambda/cooling.h
+++ b/src/cooling/const_lambda/cooling.h
@@ -296,6 +296,34 @@ INLINE static float cooling_get_temperature(
     return T_transition;
 }
 
+/**
+ * @brief Returns the subgrid temperature of a particle.
+ *
+ * This model has no subgrid quantity. We return an error.
+ *
+ * @param p The particle.
+ * @param xp The extended particle data.
+ */
+INLINE static float cooling_get_subgrid_temperature(const struct part* p,
+                                                    const struct xpart* xp) {
+  error("This cooling model does not use subgrid quantities!");
+  return -1.f;
+}
+
+/**
+ * @brief Returns the subgrid density of a particle.
+ *
+ * This model has no subgrid quantity. We return an error.
+ *
+ * @param p The particle.
+ * @param xp The extended particle data.
+ */
+INLINE static float cooling_get_subgrid_density(const struct part* p,
+                                                const struct xpart* xp) {
+  error("This cooling model does not use subgrid quantities!");
+  return -1.f;
+}
+
 /**
  * @brief Returns the total radiated energy by this particle.
  *
diff --git a/src/cooling/grackle/cooling.c b/src/cooling/grackle/cooling.c
index 84b9fa4457f5cce15ea0cdb61f530d8548c43ae0..7d4f88a013a0dd7172d126faa2f8e215675a977b 100644
--- a/src/cooling/grackle/cooling.c
+++ b/src/cooling/grackle/cooling.c
@@ -256,6 +256,34 @@ void cooling_first_init_part(const struct phys_const* restrict phys_const,
 #endif
 }
 
+/**
+ * @brief Returns the subgrid temperature of a particle.
+ *
+ * This model has no subgrid quantity. We return an error.
+ *
+ * @param p The particle.
+ * @param xp The extended particle data.
+ */
+INLINE static float cooling_get_subgrid_temperature(const struct part* p,
+                                                    const struct xpart* xp) {
+  error("This cooling model does not use subgrid quantities!");
+  return -1.f;
+}
+
+/**
+ * @brief Returns the subgrid density of a particle.
+ *
+ * This model has no subgrid quantity. We return an error.
+ *
+ * @param p The particle.
+ * @param xp The extended particle data.
+ */
+INLINE static float cooling_get_subgrid_density(const struct part* p,
+                                                const struct xpart* xp) {
+  error("This cooling model does not use subgrid quantities!");
+  return -1.f;
+}
+
 /**
  * @brief Returns the total radiated energy by this particle.
  *
diff --git a/src/cooling/grackle/cooling.h b/src/cooling/grackle/cooling.h
index bdaf0654dd7cd63c9a0e7e320adc07c086a0a3e1..4c58ef8fb1b2cddcfd965cf7a6c90aed2cbbe82d 100644
--- a/src/cooling/grackle/cooling.h
+++ b/src/cooling/grackle/cooling.h
@@ -68,6 +68,11 @@ void cooling_first_init_part(const struct phys_const* restrict phys_const,
                              const struct part* restrict p,
                              struct xpart* restrict xp);
 
+float cooling_get_subgrid_temperature(const struct part* p,
+                                      const struct xpart* xp);
+
+float cooling_get_subgrid_density(const struct part* p, const struct xpart* xp);
+
 float cooling_get_radiated_energy(const struct xpart* restrict xp);
 void cooling_print_backend(const struct cooling_function_data* cooling);
 
diff --git a/src/cooling/none/cooling.h b/src/cooling/none/cooling.h
index 92cf1fd75855a938b22368ccde2777907b1253fa..518e5f435b49e751a7bd6f0635921479248412ab 100644
--- a/src/cooling/none/cooling.h
+++ b/src/cooling/none/cooling.h
@@ -167,6 +167,32 @@ INLINE static float cooling_get_temperature(
     return T_transition;
 }
 
+/**
+ * @param Returns the subgrid temperature of a particle.
+ *
+ * This model has no subgrid quantity. We return -1.
+ *
+ * @param p The particle.
+ * @param xp The extended particle data.
+ */
+INLINE static float cooling_get_subgrid_temperature(const struct part* p,
+                                                    const struct xpart* xp) {
+  return -1.f;
+}
+
+/**
+ * @param Returns the subgrid density of a particle.
+ *
+ * This model has no subgrid quantity. We return -1.
+ *
+ * @param p The particle.
+ * @param xp The extended particle data.
+ */
+INLINE static float cooling_get_subgrid_density(const struct part* p,
+                                                const struct xpart* xp) {
+  return -1.f;
+}
+
 /**
  * @brief Returns the total radiated energy by this particle.
  *
diff --git a/src/distributed_io.c b/src/distributed_io.c
index 3ab6d39a9a519e3fa9440fb275fbc6c868692dd8..f5cf8b3d3ce64a2c7a085fca544b31cb129cdccf 100644
--- a/src/distributed_io.c
+++ b/src/distributed_io.c
@@ -74,6 +74,7 @@
  * the HDF5 file.
  * @param props The #io_props of the field to read
  * @param N The number of particles to write.
+ * @param lossy_compression Level of lossy compression to use for this field.
  * @param internal_units The #unit_system used internally
  * @param snapshot_units The #unit_system used in the snapshots
  *
diff --git a/src/entropy_floor/EAGLE/entropy_floor.h b/src/entropy_floor/EAGLE/entropy_floor.h
index 32e3828a7fa20b70630906fdb2c2a9c37f8f4c80..739f18ccdf0c3bb9a095c8583f37a1c78f617a9e 100644
--- a/src/entropy_floor/EAGLE/entropy_floor.h
+++ b/src/entropy_floor/EAGLE/entropy_floor.h
@@ -87,24 +87,21 @@ struct entropy_floor_properties {
 };
 
 /**
- * @brief Compute the entropy floor of a given #part.
+ * @brief Compute the pressure from the entropy floor at a given density
  *
- * Note that the particle is not updated!!
+ * This is the pressure exactly corresponding to the imposed EoS shape.
+ * It only matches the entropy returned by the entropy_floor() function
+ * for a neutral gas with primoridal abundance.
  *
- * @param p The #part.
+ * @param rho_phys The physical density (internal units).
+ * @param rho_com The comoving density (internal units).
  * @param cosmo The cosmological model.
  * @param props The properties of the entropy floor.
  */
-static INLINE float entropy_floor(
-    const struct part *p, const struct cosmology *cosmo,
+static INLINE float entropy_floor_gas_pressure(
+    const float rho_phys, const float rho_com, const struct cosmology *cosmo,
     const struct entropy_floor_properties *props) {
 
-  /* Comoving density in internal units */
-  const float rho_com = hydro_get_comoving_density(p);
-
-  /* Physical density in internal units */
-  const float rho_phys = hydro_get_physical_density(p, cosmo);
-
   /* Mean baryon density in co-moving internal units for over-density condition
    * (Recall cosmo->critical_density_0 is 0 in a non-cosmological run,
    * making the over-density condition a no-op) */
@@ -138,8 +135,33 @@ static INLINE float entropy_floor(
     pressure = max(pressure, pressure_Cool);
   }
 
+  return pressure;
+}
+
+/**
+ * @brief Compute the entropy floor of a given #part.
+ *
+ * Note that the particle is not updated!!
+ *
+ * @param p The #part.
+ * @param cosmo The cosmological model.
+ * @param props The properties of the entropy floor.
+ */
+static INLINE float entropy_floor(
+    const struct part *p, const struct cosmology *cosmo,
+    const struct entropy_floor_properties *props) {
+
+  /* Comoving density in internal units */
+  const float rho_com = hydro_get_comoving_density(p);
+
+  /* Physical density in internal units */
+  const float rho_phys = hydro_get_physical_density(p, cosmo);
+
+  const float pressure =
+      entropy_floor_gas_pressure(rho_phys, rho_com, cosmo, props);
+
   /* Convert to an entropy.
-   * (Recall that the entropy is the same in co-moving and phycial frames) */
+   * (Recall that the entropy is the same in co-moving and physical frames) */
   return gas_entropy_from_pressure(rho_phys, pressure);
 }
 
diff --git a/src/entropy_floor/none/entropy_floor.h b/src/entropy_floor/none/entropy_floor.h
index 7269ee59f490ed154c6b3c1ca40df4601b64f076..80eef1b3f010899274028163a184efad2d30f25c 100644
--- a/src/entropy_floor/none/entropy_floor.h
+++ b/src/entropy_floor/none/entropy_floor.h
@@ -37,6 +37,22 @@ struct phys_const;
  */
 struct entropy_floor_properties {};
 
+/**
+ * @brief Compute the pressure from the entropy floor at a given density
+ *
+ * Simply return 0 (no floor).
+ *
+ * @param rho_phys The physical density (internal units).
+ * @param rho_com The comoving density (internal units).
+ * @param cosmo The cosmological model.
+ * @param props The properties of the entropy floor.
+ */
+static INLINE float entropy_floor_gas_pressure(
+    const float rho_phys, const float rho_com, const struct cosmology *cosmo,
+    const struct entropy_floor_properties *props) {
+  return 0.f;
+}
+
 /**
  * @brief Compute the entropy floor of a given #part.
  *
@@ -53,6 +69,23 @@ static INLINE float entropy_floor(
   return 0.f;
 }
 
+/**
+ * @brief Compute the temperature from the entropy floor at a given density
+ *
+ * Simply return 0 (no floor).
+ *
+ * @param rho_phys The physical density (internal units).
+ * @param rho_com The comoving density (internal units).
+ * @param cosmo The cosmological model.
+ * @param props The properties of the entropy floor.
+ */
+static INLINE float entropy_floor_gas_temperature(
+    const float rho_phys, const float rho_com, const struct cosmology *cosmo,
+    const struct entropy_floor_properties *props) {
+
+  return 0.f;
+}
+
 /**
  * @brief Compute the temperature from the entropy floor for a given #part
  *
diff --git a/src/feedback/EAGLE/feedback.c b/src/feedback/EAGLE/feedback.c
index d6fd5a6f168fddf20fdbd385c27b786b3413ed6d..8046de77262a49d3b250359d3861dc1e8f09e58f 100644
--- a/src/feedback/EAGLE/feedback.c
+++ b/src/feedback/EAGLE/feedback.c
@@ -878,6 +878,7 @@ INLINE static void evolve_AGB(const double log10_min_mass,
  * functions to calculate feedback due to SNIa, SNII and AGB
  *
  * @param feedback_props feedback_props data structure
+ * @param phys_const The physical constants in internal units.
  * @param cosmo The cosmological model.
  * @param sp spart that we're evolving
  * @param us unit_system data structure
diff --git a/src/output_options.c b/src/output_options.c
index 50cd088c05aac1ccfb3a10abc954832153944dbd..f5e1a10575a016d9f4e845f8856fdcb715191fda 100644
--- a/src/output_options.c
+++ b/src/output_options.c
@@ -134,7 +134,7 @@ void output_options_struct_restore(struct output_options* output_options,
  *        output (i.e. top level section in the yaml).
  * @param field_name pointer to a char array containing the name of the
  *        relevant field.
- * @param part_type integer particle type
+ * @param ptype integer particle type
  * @param compression_level_current_default The default output strategy
  *.       based on the snapshot_type and part_type.
  *
@@ -143,13 +143,13 @@ void output_options_struct_restore(struct output_options* output_options,
  **/
 enum lossy_compression_schemes output_options_get_field_compression(
     const struct output_options* output_options, const char* snapshot_type,
-    const char* field_name, const enum part_type part_type,
+    const char* field_name, const enum part_type ptype,
     const enum lossy_compression_schemes compression_level_current_default) {
 
   /* Full name for the field path */
   char field[PARSER_MAX_LINE_SIZE];
   sprintf(field, "%.*s:%.*s_%s", FIELD_BUFFER_SIZE, snapshot_type,
-          FIELD_BUFFER_SIZE, field_name, part_type_names[part_type]);
+          FIELD_BUFFER_SIZE, field_name, part_type_names[ptype]);
 
   char compression_level[FIELD_BUFFER_SIZE];
   parser_get_opt_param_string(
@@ -178,16 +178,16 @@ enum lossy_compression_schemes output_options_get_field_compression(
  *
  * @param output_params The parsed select output file.
  * @param snapshot_type The type of snapshot we are writing
- * @param part_type The #part_type we are considering.
+ * @param ptype The #part_type we are considering.
  */
 enum lossy_compression_schemes output_options_get_ptype_default_compression(
     struct swift_params* output_params, const char* snapshot_type,
-    const enum part_type part_type) {
+    const enum part_type ptype) {
 
   /* Full name for the default path */
   char field[PARSER_MAX_LINE_SIZE];
   sprintf(field, "%.*s:Standard_%s", FIELD_BUFFER_SIZE, snapshot_type,
-          part_type_names[part_type]);
+          part_type_names[ptype]);
 
   char compression_level[FIELD_BUFFER_SIZE];
   parser_get_opt_param_string(
@@ -210,7 +210,7 @@ enum lossy_compression_schemes output_options_get_ptype_default_compression(
         "A lossy default compression strategy was specified for snapshot "
         "type %s and particle type %d. This is not allowed, lossy "
         "compression must be set on a field-by-field basis.",
-        snapshot_type, part_type);
+        snapshot_type, ptype);
 
 #ifdef SWIFT_DEBUG_CHECKS
   /* Check whether we could translate the level string to a known entry. */
@@ -218,13 +218,12 @@ enum lossy_compression_schemes output_options_get_ptype_default_compression(
     error(
         "Could not resolve compression level \"%s\" as default compression "
         "level of particle type %s in snapshot type %s.",
-        compression_level, part_type_names[part_type], snapshot_type);
+        compression_level, part_type_names[ptype], snapshot_type);
 
   message(
       "Determined default compression level of %s in snapshot type %s "
       "as \"%s\", corresponding to level code %d",
-      part_type_names[part_type], snapshot_type, compression_level,
-      level_index);
+      part_type_names[ptype], snapshot_type, compression_level, level_index);
 #endif
 
   return (enum lossy_compression_schemes)level_index;
diff --git a/src/star_formation.c b/src/star_formation.c
index 60cff1e2e68feaf7e71705b5079294ec478fad42..c0e8dfc2049c62ff827c7caa945bb08ee46e92e8 100644
--- a/src/star_formation.c
+++ b/src/star_formation.c
@@ -35,16 +35,21 @@
  * @param phys_const Physical constants in internal units
  * @param us the current internal system of units
  * @param hydro_props The propoerties of the hydro scheme.
+ * @param cosmo The cosmology model.
+ * @param entropy_floor The properties of the entropy floor used in this
+ * simulation.
  * @param starform the properties of the star formation law
  */
 void starformation_init(struct swift_params* parameter_file,
                         const struct phys_const* phys_const,
                         const struct unit_system* us,
                         const struct hydro_props* hydro_props,
+                        const struct cosmology* cosmo,
+                        const struct entropy_floor_properties* entropy_floor,
                         struct star_formation* starform) {
 
-  starformation_init_backend(parameter_file, phys_const, us, hydro_props,
-                             starform);
+  starformation_init_backend(parameter_file, phys_const, us, hydro_props, cosmo,
+                             entropy_floor, starform);
 }
 
 /**
diff --git a/src/star_formation.h b/src/star_formation.h
index 3ac28ad02b557aa06ba03e9c43da0d6f511a56f7..cd21c8aa43c5c6a0a6747c3bd1e7c3544e3cc4e9 100644
--- a/src/star_formation.h
+++ b/src/star_formation.h
@@ -52,6 +52,8 @@ void starformation_init(struct swift_params* parameter_file,
                         const struct phys_const* phys_const,
                         const struct unit_system* us,
                         const struct hydro_props* hydro_props,
+                        const struct cosmology* cosmo,
+                        const struct entropy_floor_properties* entropy_floor,
                         struct star_formation* starform);
 
 void starformation_print(const struct star_formation* starform);
diff --git a/src/star_formation/EAGLE/star_formation.h b/src/star_formation/EAGLE/star_formation.h
index 57f97089289d089404f13effd876ea616ed4ad8f..318c2610e12c8d4c6561b8f81de22aa5e27f4e9a 100644
--- a/src/star_formation/EAGLE/star_formation.h
+++ b/src/star_formation/EAGLE/star_formation.h
@@ -51,11 +51,23 @@ enum star_formation_law {
   eagle_star_formation_pressure_law /*<! Pressure law */
 };
 
+/**
+ * @brief Choice of star formation threshold
+ */
+enum star_formation_threshold {
+  eagle_star_formation_threshold_Z_dep,  /*<! SF threshold based on Schaye+2004
+                                            Z-dependence */
+  eagle_star_formation_threshold_subgrid /*<! SF threshold based on subgrid
+                                            properties */
+};
+
 /**
  * @brief Properties of the EAGLE star formation model.
  */
 struct star_formation {
 
+  /* SF law ------------------------------------------------------------*/
+
   /*! Which SF law are we using? */
   enum star_formation_law SF_law;
 
@@ -112,141 +124,177 @@ struct star_formation {
 
   } pressure_law;
 
-  /* SF threshold properties -------------------------------------------- */
-
-  /*! Critical overdensity */
-  double min_over_den;
+  /* Density for direct conversion to star -------------------------------- */
 
-  /*! Density threshold to form stars (internal units) */
-  double density_threshold;
+  /*! Max physical density (H atoms per cm^3)*/
+  double gas_density_direct_HpCM3;
 
-  /*! Density threshold to form stars in user units */
-  double density_threshold_HpCM3;
+  /*! Max physical density (internal units) */
+  double gas_density_direct;
 
-  /*! Maximum density threshold to form stars (internal units) */
-  double density_threshold_max;
+  /* SF threshold --------------------------------------------------------- */
 
-  /*! Maximum density threshold to form stars (H atoms per cm^3) */
-  double density_threshold_max_HpCM3;
+  enum star_formation_threshold SF_threshold;
 
-  /*! Reference metallicity for metal-dependant threshold */
-  double Z0;
+  /*! Critical overdensity above which SF is allowed */
+  double min_over_den;
 
-  /*! Inverse of reference metallicity */
-  double Z0_inv;
+  struct {
 
-  /*! critical density Metallicity power law (internal units) */
-  double n_Z0;
+    /*! Density threshold to form stars (internal units) */
+    double density_threshold;
 
-  /* Internal EoS properties ----------------------------------------------- */
+    /*! Density threshold to form stars in user units */
+    double density_threshold_HpCM3;
 
-  /*! Polytropic index */
-  double EOS_polytropic_index;
+    /*! Maximum density threshold to form stars (internal units) */
+    double density_threshold_max;
 
-  /*! EOS density norm (H atoms per cm^3) */
-  double EOS_density_norm_HpCM3;
+    /*! Maximum density threshold to form stars (H atoms per cm^3) */
+    double density_threshold_max_HpCM3;
 
-  /*! EOS Temperature norm (Kelvin)  */
-  double EOS_temperature_norm_K;
+    /*! Reference metallicity for metal-dependant threshold */
+    double Z0;
 
-  /*! EOS pressure norm, eq. 13 of Schaye & Dalla Vecchia 2008 (internal units)
-   */
-  double EOS_pressure_c;
+    /*! Inverse of reference metallicity */
+    double Z0_inv;
 
-  /*! EOS Temperarure norm, eq. 13 of Schaye & Dalla Vecchia 2008 (internal
-   * units) */
-  double EOS_temperature_c;
+    /*! critical density Metallicity power law (internal units) */
+    double n_Z0;
 
-  /*! EOS density norm, eq. 13 of Schaye & Dalla Vecchia 2008 (internal units)
-   */
-  double EOS_density_c;
+    /*! Dalla Vecchia & Schaye entropy differnce criterion */
+    double entropy_margin_threshold_dex;
 
-  /*! Inverse of EOS density norm (internal units) */
-  double EOS_density_c_inv;
+    /*! 10^Tdex of Dalla Vecchia & Schaye entropy difference criterion */
+    double ten_to_entropy_margin_threshold_dex;
 
-  /* Maximal offset from entropy floor allowed for SF --------------------- */
+  } Z_dep_thresh;
 
-  /*! Dalla Vecchia & Schaye entropy differnce criterion */
-  double entropy_margin_threshold_dex;
+  struct {
 
-  /*! 10^Tdex of Dalla Vecchia & Schaye entropy difference criterion */
-  double ten_to_entropy_margin_threshold_dex;
+    /*! (Subgrid) temperature threshold for SF to use on its own */
+    double T_threshold1;
 
-  /* Density for direct conversion to star -------------------------------- */
+    /*! (Subgrid) temperature threshold for SF to use combined with the density
+     * threshold */
+    double T_threshold2;
 
-  /*! Max physical density (H atoms per cm^3)*/
-  double gas_density_direct_HpCM3;
+    /*! (Subgrid) Hydrogen number density threshold for SF */
+    double nH_threshold;
 
-  /*! Max physical density (internal units) */
-  double gas_density_direct;
+  } subgrid_thresh;
 };
 
 /**
- * @brief Computes the density threshold for star-formation fo a given total
- * metallicity.
- *
- * Follows Schaye (2004) eq. 19 and 24 (see also Schaye et al. 2015, eq. 2).
+ * @brief Calculate if the satisfies the conditions for star formation.
  *
- * @param Z The metallicity (metal mass fraction).
- * @param starform The properties of the star formation model.
- * @param phys_const The physical constants.
- * @return The physical density threshold for star formation in internal units.
+ * @param starform the star formation law properties to use.
+ * @param p the gas particles.
+ * @param xp the additional properties of the gas particles.
+ * @param phys_const the physical constants in internal units.
+ * @param cosmo the cosmological parameters and properties.
+ * @param hydro_props The properties of the hydro scheme.
+ * @param us The internal system of units.
+ * @param cooling The cooling data struct.
+ * @param entropy_floor_props The entropy floor assumed in this run.
  */
-INLINE static double star_formation_threshold(
-    const double Z, const struct star_formation* starform,
-    const struct phys_const* phys_const) {
+INLINE static int star_formation_is_star_forming_Z_dep(
+    const struct part* p, const struct xpart* xp,
+    const struct star_formation* starform, const struct phys_const* phys_const,
+    const struct cosmology* cosmo, const struct hydro_props* hydro_props,
+    const struct unit_system* us, const struct cooling_function_data* cooling,
+    const struct entropy_floor_properties* entropy_floor_props) {
+
+  /* Physical density of the particle */
+  const double physical_density = hydro_get_physical_density(p, cosmo);
+
+  /* Get the Hydrogen number density (assuming primordial H abundance) */
+  const double n_H = physical_density * hydro_props->hydrogen_mass_fraction;
+
+  /* Get the density threshold for star formation */
+  const double Z =
+      chemistry_get_total_metal_mass_fraction_for_star_formation(p);
 
   double density_threshold;
 
   /* Schaye (2004), eq. 19 and 24 */
   if (Z > 0.) {
-    density_threshold = starform->density_threshold *
-                        powf(Z * starform->Z0_inv, starform->n_Z0);
-    density_threshold = min(density_threshold, starform->density_threshold_max);
+    density_threshold =
+        starform->Z_dep_thresh.density_threshold *
+        powf(Z * starform->Z_dep_thresh.Z0_inv, starform->Z_dep_thresh.n_Z0);
+    density_threshold =
+        min(density_threshold, starform->Z_dep_thresh.density_threshold_max);
   } else {
-    density_threshold = starform->density_threshold_max;
+    density_threshold = starform->Z_dep_thresh.density_threshold_max;
   }
 
   /* Convert to mass density */
-  return density_threshold * phys_const->const_proton_mass;
-}
+  density_threshold *= phys_const->const_proton_mass;
 
-/**
- * @brief Compute the pressure on the polytropic equation of state for a given
- * Hydrogen number density.
- *
- * Schaye & Dalla Vecchia 2008, eq. 13.
- *
- * @param n_H The Hydrogen number density in internal units.
- * @param starform The properties of the star formation model.
- * @return The pressure on the equation of state in internal units.
- */
-INLINE static double EOS_pressure(const double n_H,
-                                  const struct star_formation* starform) {
+  /* Check if it exceeded the minimum density */
+  if (n_H < density_threshold) return 0;
 
-  return starform->EOS_pressure_c *
-         pow(n_H * starform->EOS_density_c_inv, starform->EOS_polytropic_index);
+  /* Calculate the entropy of the particle */
+  const double entropy = hydro_get_physical_entropy(p, xp, cosmo);
+
+  /* Calculate the entropy that will be used to calculate
+   * the off-set from the EoS */
+  const double entropy_eos = entropy_floor(p, cosmo, entropy_floor_props);
+
+  /* Check the Schaye & Dalla Vecchia 2012 EOS-based temperature criterion */
+  return (entropy <
+          entropy_eos *
+              starform->Z_dep_thresh.ten_to_entropy_margin_threshold_dex);
 }
 
 /**
- * @brief Compute the entropy of the polytropic equation of state for a given
- * Hydrogen number density.
+ * @brief Calculate if the satisfies the conditions for star formation.
  *
- * @param n_H The Hydrogen number density in internal units.
- * @param starform The properties of the star formation model.
- * @param rho The physical density
- * @return The pressure on the equation of state in internal units.
+ * @param starform the star formation law properties to use.
+ * @param p the gas particles.
+ * @param xp the additional properties of the gas particles.
+ * @param phys_const the physical constants in internal units.
+ * @param cosmo the cosmological parameters and properties.
+ * @param hydro_props The properties of the hydro scheme.
+ * @param us The internal system of units.
+ * @param cooling The cooling data struct.
+ * @param entropy_floor_props The entropy floor assumed in this run.
  */
-INLINE static double EOS_entropy(const double n_H,
-                                 const struct star_formation* starform,
-                                 const double rho) {
-
-  return gas_entropy_from_pressure(rho, EOS_pressure(n_H, starform));
+INLINE static int star_formation_is_star_forming_subgrid(
+    const struct part* p, const struct xpart* xp,
+    const struct star_formation* starform, const struct phys_const* phys_const,
+    const struct cosmology* cosmo, const struct hydro_props* hydro_props,
+    const struct unit_system* us, const struct cooling_function_data* cooling,
+    const struct entropy_floor_properties* entropy_floor_props) {
+
+  const double number_density_to_cgs =
+      units_cgs_conversion_factor(us, UNIT_CONV_NUMBER_DENSITY);
+
+  /* Get the Hydrogen mass fraction */
+  const double XH = chemistry_get_metal_mass_fraction_for_star_formation(
+      p)[chemistry_element_H];
+
+  /* Get the subgrid properties
+   * Note these are both in physical frame already */
+  const double subgrid_T_cgs = cooling_get_subgrid_temperature(p, xp);
+  const double subgrid_rho = cooling_get_subgrid_density(p, xp);
+  const double subgrid_n_H = subgrid_rho * XH / phys_const->const_proton_mass;
+  const double subgrid_n_H_cgs = subgrid_n_H * number_density_to_cgs;
+
+  /* Now, determine whether we are very cold or (cold and dense) enough
+   *
+   * This would typically be (T < 10^3 OR (T < 10^4.5 AND n_H > 10))
+   * with T and n_H subgrid properties.
+   *
+   * Recall that particles above the EoS have T_sub = T and rho_sub = rho.
+   */
+  return ((subgrid_T_cgs < starform->subgrid_thresh.T_threshold1) ||
+          (subgrid_T_cgs < starform->subgrid_thresh.T_threshold2 &&
+           subgrid_n_H_cgs > starform->subgrid_thresh.nH_threshold));
 }
 
 /**
- * @brief Calculate if the gas has the potential of becoming
- * a star.
+ * @brief Calculate if the satisfies the conditions for star formation.
  *
  * @param starform the star formation law properties to use.
  * @param p the gas particles.
@@ -259,58 +307,44 @@ INLINE static double EOS_entropy(const double n_H,
  * @param entropy_floor_props The entropy floor assumed in this run.
  */
 INLINE static int star_formation_is_star_forming(
-    const struct part* restrict p, const struct xpart* restrict xp,
+    const struct part* p, const struct xpart* xp,
     const struct star_formation* starform, const struct phys_const* phys_const,
-    const struct cosmology* cosmo,
-    const struct hydro_props* restrict hydro_props,
-    const struct unit_system* restrict us,
-    const struct cooling_function_data* restrict cooling,
-    const struct entropy_floor_properties* restrict entropy_floor_props) {
+    const struct cosmology* cosmo, const struct hydro_props* hydro_props,
+    const struct unit_system* us, const struct cooling_function_data* cooling,
+    const struct entropy_floor_properties* entropy_floor_props) {
 
-  /* Minimal density (converted from mean baryonic density) for star formation
-   */
+  /* Minimal density (converted from mean baryonic density)
+   * for star formation */
   const double rho_mean_b_times_min_over_den =
       cosmo->mean_density_Omega_b * starform->min_over_den;
 
   /* Physical density of the particle */
   const double physical_density = hydro_get_physical_density(p, cosmo);
 
-  /* Deside whether we should form stars or not,
-   * first we deterime if we have the correct over density
-   * if that is true we calculate if either the maximum density
-   * threshold is reached or if the metallicity dependent
-   * threshold is reached, after this we calculate if the
-   * temperature is appropriate */
-  if (physical_density < rho_mean_b_times_min_over_den) return 0;
-
-  /* In this case there are actually multiple possibilities
-   * because we also need to check if the physical density exceeded
-   * the appropriate limit */
-
-  /* Get the Hydrogen number density (assuming primordial H abundance) */
-  const double n_H = physical_density * hydro_props->hydrogen_mass_fraction;
+  /* Deside whether we should form stars or not */
 
-  /* Get the density threshold for star formation */
-  const double Z =
-      chemistry_get_total_metal_mass_fraction_for_star_formation(p);
-  const double density_threshold =
-      star_formation_threshold(Z, starform, phys_const);
+  /* First, deterime if we have the correct over density */
+  if (physical_density < rho_mean_b_times_min_over_den) return 0;
 
-  /* Check if it exceeded the minimum density */
-  if (n_H < density_threshold) return 0;
+  /* Determine which star formation model to use */
+  switch (starform->SF_threshold) {
 
-  /* Calculate the entropy of the particle */
-  const double entropy = hydro_get_physical_entropy(p, xp, cosmo);
+    case eagle_star_formation_threshold_Z_dep:
+      return star_formation_is_star_forming_Z_dep(p, xp, starform, phys_const,
+                                                  cosmo, hydro_props, us,
+                                                  cooling, entropy_floor_props);
+      break;
 
-  /* Calculate the entropy that will be used to calculate
-   * the off-set, this is the maximum between the entropy
-   * floor and the star formation polytropic EOS. */
-  const double entropy_eos = max(entropy_floor(p, cosmo, entropy_floor_props),
-                                 EOS_entropy(n_H, starform, physical_density));
+    case eagle_star_formation_threshold_subgrid:
+      return star_formation_is_star_forming_subgrid(
+          p, xp, starform, phys_const, cosmo, hydro_props, us, cooling,
+          entropy_floor_props);
+      break;
 
-  /* Check the Schaye & Dalla Vecchia 2012 EOS-based temperature criterion */
-  return (entropy <
-          entropy_eos * starform->ten_to_entropy_margin_threshold_dex);
+    default:
+      error("Invalid star formation threshold model!!!");
+      return 0;
+  }
 }
 
 /**
@@ -370,15 +404,12 @@ INLINE static void star_formation_compute_SFR_pressure_law(
   /* Hydrogen number density of this particle (assuming primordial H abundance)
    */
   const double physical_density = hydro_get_physical_density(p, cosmo);
-  const double n_H = physical_density * hydro_props->hydrogen_mass_fraction;
 
   /* Get the pressure used for the star formation, this is
-   * the maximum of the star formation EOS pressure,
-   * the physical pressure of the particle and the
+   * the maximum the physical pressure of the particle and the
    * floor pressure. The floor pressure is used implicitly
    * when getting the physical pressure. */
-  const double pressure =
-      max(EOS_pressure(n_H, starform), hydro_get_physical_pressure(p, cosmo));
+  const double pressure = hydro_get_physical_pressure(p, cosmo);
 
   /* Calculate the specific star formation rate */
   double SFRpergasmass;
@@ -411,7 +442,7 @@ INLINE static void star_formation_compute_SFR_pressure_law(
  * @param dt_star The time-step of this particle.
  */
 INLINE static void star_formation_compute_SFR(
-    const struct part* restrict p, struct xpart* restrict xp,
+    const struct part* p, struct xpart* xp,
     const struct star_formation* starform, const struct phys_const* phys_const,
     const struct hydro_props* hydro_props, const struct cosmology* cosmo,
     const double dt_star) {
@@ -542,10 +573,8 @@ INLINE static void star_formation_copy_properties(
     const struct part* p, const struct xpart* xp, struct spart* sp,
     const struct engine* e, const struct star_formation* starform,
     const struct cosmology* cosmo, const int with_cosmology,
-    const struct phys_const* phys_const,
-    const struct hydro_props* restrict hydro_props,
-    const struct unit_system* restrict us,
-    const struct cooling_function_data* restrict cooling,
+    const struct phys_const* phys_const, const struct hydro_props* hydro_props,
+    const struct unit_system* us, const struct cooling_function_data* cooling,
     const int convert_part) {
 
   /* Store the current mass */
@@ -589,22 +618,21 @@ INLINE static void star_formation_copy_properties(
  * @param phys_const Physical constants in internal units
  * @param us The current internal system of units.
  * @param hydro_props The propertis of the hydro model.
+ * @param cosmo The current cosmological model.
+ * @param entropy_floor The properties of the entropy floor used in this
+ * simulation.
  * @param starform the star formation law properties to initialize
  */
 INLINE static void starformation_init_backend(
     struct swift_params* parameter_file, const struct phys_const* phys_const,
     const struct unit_system* us, const struct hydro_props* hydro_props,
+    const struct cosmology* cosmo,
+    const struct entropy_floor_properties* entropy_floor,
     struct star_formation* starform) {
 
   /* Get the Gravitational constant */
   const double G_newton = phys_const->const_newton_G;
 
-  /* Initial Hydrogen abundance (mass fraction) */
-  const double X_H = hydro_props->hydrogen_mass_fraction;
-
-  /* Mean molecular weight assuming neutral gas */
-  const double mean_molecular_weight = hydro_props->mu_neutral;
-
   /* Get the surface density unit Msun / pc^2 in internal units */
   const double Msun_per_pc2 =
       phys_const->const_solar_mass /
@@ -619,28 +647,6 @@ INLINE static void starformation_init_backend(
   const double number_density_from_cgs =
       1. / units_cgs_conversion_factor(us, UNIT_CONV_NUMBER_DENSITY);
 
-  /* Load the equation of state for this model */
-  starform->EOS_polytropic_index = parser_get_param_double(
-      parameter_file, "EAGLEStarFormation:EOS_gamma_effective");
-  starform->EOS_temperature_norm_K = parser_get_param_double(
-      parameter_file, "EAGLEStarFormation:EOS_temperature_norm_K");
-  starform->EOS_density_norm_HpCM3 = parser_get_param_double(
-      parameter_file, "EAGLEStarFormation:EOS_density_norm_H_p_cm3");
-  starform->EOS_density_c =
-      starform->EOS_density_norm_HpCM3 * number_density_from_cgs;
-  starform->EOS_density_c_inv = 1. / starform->EOS_density_c;
-
-  /* Calculate the EOS pressure normalization */
-  starform->EOS_pressure_c =
-      starform->EOS_density_c * starform->EOS_temperature_norm_K *
-      phys_const->const_boltzmann_k / mean_molecular_weight / X_H;
-
-  /* Normalisation of the temperature in the EOS calculatio */
-  starform->EOS_temperature_c =
-      starform->EOS_pressure_c / phys_const->const_boltzmann_k;
-  starform->EOS_temperature_c *=
-      pow(starform->EOS_density_c, starform->EOS_polytropic_index);
-
   /* Check if we are using the Schmidt law for the star formation rate,
    * defaults to pressure law if is not explicitely set to a Schmidt law */
   char temp[32];
@@ -715,9 +721,15 @@ INLINE static void starformation_init_backend(
         starform->pressure_law.KS_high_den_thresh_HpCM3 *
         number_density_from_cgs;
 
-    /* Pressure at the high-density threshold */
+    /* Pressure on the entropy floor at the high-density threshold
+     *
+     * Note that we use FLT_MAX as the comoving density to make sure
+     * the floor is applied no matter what redshift we are at. This will
+     * always be a density above the comoving density threashold for the floor
+     * to be used.*/
     const double EOS_high_den_pressure =
-        EOS_pressure(starform->pressure_law.KS_high_den_thresh, starform);
+        entropy_floor_gas_pressure(starform->pressure_law.KS_high_den_thresh,
+                                   FLT_MAX, cosmo, entropy_floor);
 
     /* Calculate the KS high density normalization
      * We want the SF law to be continous so the normalisation of the second
@@ -756,37 +768,70 @@ INLINE static void starformation_init_backend(
   starform->gas_density_direct =
       starform->gas_density_direct_HpCM3 * number_density_from_cgs;
 
-  starform->entropy_margin_threshold_dex = parser_get_opt_param_double(
-      parameter_file, "EAGLEStarFormation:EOS_entropy_margin_dex", FLT_MAX);
+  /* Check if we are using the Schmidt law for the star formation rate,
+   * defaults to pressure law if is not explicitely set to a Schmidt law */
+  char temp_SF[32];
+  parser_get_param_string(parameter_file, "EAGLEStarFormation:SF_threshold",
+                          temp_SF);
+
+  if (strcmp(temp_SF, "Zdep") == 0) {
+
+    /* Z-dep (Schaye+2004) model */
+    starform->SF_threshold = eagle_star_formation_threshold_Z_dep;
+
+    starform->Z_dep_thresh.entropy_margin_threshold_dex =
+        parser_get_opt_param_double(parameter_file,
+                                    "EAGLEStarFormation:EOS_entropy_margin_dex",
+                                    FLT_MAX);
 
-  starform->ten_to_entropy_margin_threshold_dex =
-      exp10(starform->entropy_margin_threshold_dex);
+    starform->Z_dep_thresh.ten_to_entropy_margin_threshold_dex =
+        exp10(starform->Z_dep_thresh.entropy_margin_threshold_dex);
 
-  /* Read the normalization of the metallicity dependent critical
-   * density*/
-  starform->density_threshold_HpCM3 = parser_get_param_double(
-      parameter_file, "EAGLEStarFormation:threshold_norm_H_p_cm3");
+    /* Read the normalization of the metallicity dependent critical
+     * density*/
+    starform->Z_dep_thresh.density_threshold_HpCM3 = parser_get_param_double(
+        parameter_file, "EAGLEStarFormation:threshold_norm_H_p_cm3");
 
-  /* Convert to internal units */
-  starform->density_threshold =
-      starform->density_threshold_HpCM3 * number_density_from_cgs;
+    /* Convert to internal units */
+    starform->Z_dep_thresh.density_threshold =
+        starform->Z_dep_thresh.density_threshold_HpCM3 *
+        number_density_from_cgs;
+
+    /* Read the scale metallicity Z0 */
+    starform->Z_dep_thresh.Z0 = parser_get_param_double(
+        parameter_file, "EAGLEStarFormation:threshold_Z0");
+    starform->Z_dep_thresh.Z0_inv = 1. / starform->Z_dep_thresh.Z0;
+
+    /* Read the power law of the critical density scaling */
+    starform->Z_dep_thresh.n_Z0 = parser_get_param_double(
+        parameter_file, "EAGLEStarFormation:threshold_slope");
+
+    /* Read the maximum allowed density for star formation */
+    starform->Z_dep_thresh.density_threshold_max_HpCM3 =
+        parser_get_param_double(
+            parameter_file, "EAGLEStarFormation:threshold_max_density_H_p_cm3");
+
+    /* Convert to internal units */
+    starform->Z_dep_thresh.density_threshold_max =
+        starform->Z_dep_thresh.density_threshold_max_HpCM3 *
+        number_density_from_cgs;
 
-  /* Read the scale metallicity Z0 */
-  starform->Z0 = parser_get_param_double(parameter_file,
-                                         "EAGLEStarFormation:threshold_Z0");
-  starform->Z0_inv = 1. / starform->Z0;
+  } else if (strcmp(temp_SF, "Subgrid") == 0) {
 
-  /* Read the power law of the critical density scaling */
-  starform->n_Z0 = parser_get_param_double(
-      parameter_file, "EAGLEStarFormation:threshold_slope");
+    /* Subgrid quantities based model */
+    starform->SF_threshold = eagle_star_formation_threshold_subgrid;
 
-  /* Read the maximum allowed density for star formation */
-  starform->density_threshold_max_HpCM3 = parser_get_param_double(
-      parameter_file, "EAGLEStarFormation:threshold_max_density_H_p_cm3");
+    /* Read threshold properties */
+    starform->subgrid_thresh.T_threshold1 = parser_get_param_double(
+        parameter_file, "EAGLEStarFormation:threshold_temperature1_K");
+    starform->subgrid_thresh.T_threshold2 = parser_get_param_double(
+        parameter_file, "EAGLEStarFormation:threshold_temperature2_K");
+    starform->subgrid_thresh.nH_threshold = parser_get_param_double(
+        parameter_file, "EAGLEStarFormation:threshold_number_density_H_p_cm3");
 
-  /* Convert to internal units */
-  starform->density_threshold_max =
-      starform->density_threshold_max_HpCM3 * number_density_from_cgs;
+  } else {
+    error("Invalid SF threshold model: '%s'", temp_SF);
+  }
 }
 
 /**
@@ -797,7 +842,41 @@ INLINE static void starformation_init_backend(
 INLINE static void starformation_print_backend(
     const struct star_formation* starform) {
 
-  message("Star formation law is EAGLE");
+  message("Star formation model is EAGLE");
+
+  switch (starform->SF_threshold) {
+    case eagle_star_formation_threshold_Z_dep:
+
+      message("Density threshold follows Schaye (2004)");
+      message(
+          "the normalization of the density threshold is given by"
+          " %e #/cm^3, with metallicity slope of %e, and metallicity "
+          "normalization"
+          " of %e, the maximum density threshold is given by %e #/cm^3",
+          starform->Z_dep_thresh.density_threshold_HpCM3,
+          starform->Z_dep_thresh.n_Z0, starform->Z_dep_thresh.Z0,
+          starform->Z_dep_thresh.density_threshold_max_HpCM3);
+      message(
+          "Temperature threshold is given by Dalla Vecchia and Schaye (2012)");
+      message(
+          "The temperature threshold offset from the EOS is given by: %e dex",
+          starform->Z_dep_thresh.entropy_margin_threshold_dex);
+
+      break;
+    case eagle_star_formation_threshold_subgrid:
+
+      message("Density threshold uses subgrid quantities");
+      message(
+          "Particles are star-forming if their properties obey (T_sub < %e K "
+          "OR (T_sub < %e K AND n_H,sub > %e cm^-3))",
+          starform->subgrid_thresh.T_threshold1,
+          starform->subgrid_thresh.T_threshold2,
+          starform->subgrid_thresh.nH_threshold);
+
+      break;
+    default:
+      error("Invalid star formation threshold!!!");
+  }
 
   switch (starform->SF_law) {
     case eagle_star_formation_schmidt_law:
@@ -819,25 +898,9 @@ INLINE static void starformation_print_backend(
               starform->pressure_law.KS_high_den_power_law);
       break;
     default:
-      error("Invalid star formation model!!!");
+      error("Invalid star formation law!!!");
   }
 
-  message(
-      "The effective equation of state is given by: polytropic "
-      "index = %e , normalization density = %e #/cm^3 and normalization "
-      "temperature = %e K",
-      starform->EOS_polytropic_index, starform->EOS_density_norm_HpCM3,
-      starform->EOS_temperature_norm_K);
-  message("Density threshold follows Schaye (2004)");
-  message(
-      "the normalization of the density threshold is given by"
-      " %e #/cm^3, with metallicity slope of %e, and metallicity normalization"
-      " of %e, the maximum density threshold is given by %e #/cm^3",
-      starform->density_threshold_HpCM3, starform->n_Z0, starform->Z0,
-      starform->density_threshold_max_HpCM3);
-  message("Temperature threshold is given by Dalla Vecchia and Schaye (2012)");
-  message("The temperature threshold offset from the EOS is given by: %e dex",
-          starform->entropy_margin_threshold_dex);
   message("Running with a direct conversion density of: %e #/cm^3",
           starform->gas_density_direct_HpCM3);
 }
@@ -854,8 +917,8 @@ INLINE static void starformation_print_backend(
  * @param cosmo The current cosmological model.
  */
 __attribute__((always_inline)) INLINE static void star_formation_end_density(
-    struct part* restrict p, struct xpart* restrict xp,
-    const struct star_formation* cd, const struct cosmology* cosmo) {}
+    struct part* p, struct xpart* xp, const struct star_formation* cd,
+    const struct cosmology* cosmo) {}
 
 /**
  * @brief Sets all particle fields to sensible values when the #part has 0 ngbs.
@@ -869,8 +932,7 @@ __attribute__((always_inline)) INLINE static void star_formation_end_density(
  * @param cosmo The current cosmological model.
  */
 __attribute__((always_inline)) INLINE static void
-star_formation_part_has_no_neighbours(struct part* restrict p,
-                                      struct xpart* restrict xp,
+star_formation_part_has_no_neighbours(struct part* p, struct xpart* xp,
                                       const struct star_formation* cd,
                                       const struct cosmology* cosmo) {}
 
@@ -885,7 +947,7 @@ star_formation_part_has_no_neighbours(struct part* restrict p,
  * @param p Pointer to the particle data.
  */
 __attribute__((always_inline)) INLINE static void star_formation_init_part(
-    struct part* restrict p, const struct star_formation* data) {}
+    struct part* p, const struct star_formation* data) {}
 
 /**
  * @brief Sets the star_formation properties of the (x-)particles to a valid
@@ -901,12 +963,11 @@ __attribute__((always_inline)) INLINE static void star_formation_init_part(
  * @param xp Pointer to the extended particle data.
  */
 __attribute__((always_inline)) INLINE static void
-star_formation_first_init_part(const struct phys_const* restrict phys_const,
-                               const struct unit_system* restrict us,
-                               const struct cosmology* restrict cosmo,
+star_formation_first_init_part(const struct phys_const* phys_const,
+                               const struct unit_system* us,
+                               const struct cosmology* cosmo,
                                const struct star_formation* data,
-                               const struct part* restrict p,
-                               struct xpart* restrict xp) {}
+                               const struct part* p, struct xpart* xp) {}
 
 /**
  * @brief Split the star formation content of a particle into n pieces
diff --git a/src/star_formation/GEAR/star_formation_io.h b/src/star_formation/GEAR/star_formation_io.h
index 46ac7c6e3de567fab148b08a62d8d88e1e8113da..e5a7231d9678e0b5b8396e87d85789551719231b 100644
--- a/src/star_formation/GEAR/star_formation_io.h
+++ b/src/star_formation/GEAR/star_formation_io.h
@@ -104,12 +104,16 @@ star_formation_write_sparticles(const struct spart* sparts,
  * @param phys_const Physical constants in internal units
  * @param us The current internal system of units
  * @param hydro_props The #hydro_props.
+ * @param cosmo The current cosmological model.
+ * @param entropy_floor The properties of the entropy floor used in this
+ * simulation.
  * @param starform the star formation law properties to initialize
- *
  */
 INLINE static void starformation_init_backend(
     struct swift_params* parameter_file, const struct phys_const* phys_const,
     const struct unit_system* us, const struct hydro_props* hydro_props,
+    const struct cosmology* cosmo,
+    const struct entropy_floor_properties* entropy_floor,
     struct star_formation* starform) {
 
   /* Star formation efficiency */
diff --git a/src/star_formation/QLA/star_formation.h b/src/star_formation/QLA/star_formation.h
index 5f8a36f2103256464d7a3f469da09fc43a4b5b2a..dd0244bf7731c7809a3003502e8b6e9c75f6b370 100644
--- a/src/star_formation/QLA/star_formation.h
+++ b/src/star_formation/QLA/star_formation.h
@@ -94,7 +94,7 @@ INLINE static int star_formation_is_star_forming(
  * @param dt_star The time-step of this particle.
  */
 INLINE static void star_formation_compute_SFR(
-    const struct part* restrict p, struct xpart* restrict xp,
+    const struct part* p, struct xpart* xp,
     const struct star_formation* starform, const struct phys_const* phys_const,
     const struct hydro_props* hydro_props, const struct cosmology* cosmo,
     const double dt_star) {
@@ -159,10 +159,8 @@ INLINE static void star_formation_copy_properties(
     const struct part* p, const struct xpart* xp, struct spart* sp,
     const struct engine* e, const struct star_formation* starform,
     const struct cosmology* cosmo, const int with_cosmology,
-    const struct phys_const* phys_const,
-    const struct hydro_props* restrict hydro_props,
-    const struct unit_system* restrict us,
-    const struct cooling_function_data* restrict cooling,
+    const struct phys_const* phys_const, const struct hydro_props* hydro_props,
+    const struct unit_system* us, const struct cooling_function_data* cooling,
     const int convert_part) {
 
   /* Store the current mass */
@@ -183,11 +181,16 @@ INLINE static void star_formation_copy_properties(
  * @param phys_const Physical constants in internal units
  * @param us The current internal system of units.
  * @param hydro_props The propertis of the hydro model.
+ * @param cosmo The current cosmological model.
+ * @param entropy_floor The properties of the entropy floor used in this
+ * simulation.
  * @param starform the star formation law properties to initialize
  */
 INLINE static void starformation_init_backend(
     struct swift_params* parameter_file, const struct phys_const* phys_const,
     const struct unit_system* us, const struct hydro_props* hydro_props,
+    const struct cosmology* cosmo,
+    const struct entropy_floor_properties* entropy_floor,
     struct star_formation* starform) {
 
   /* Read the critical density contrast from the parameter file*/
@@ -219,8 +222,8 @@ INLINE static void starformation_print_backend(
  * @param cosmo The current cosmological model.
  */
 __attribute__((always_inline)) INLINE static void star_formation_end_density(
-    struct part* restrict p, struct xpart* restrict xp,
-    const struct star_formation* cd, const struct cosmology* cosmo) {}
+    struct part* p, struct xpart* xp, const struct star_formation* cd,
+    const struct cosmology* cosmo) {}
 
 /**
  * @brief Sets all particle fields to sensible values when the #part has 0 ngbs.
@@ -234,8 +237,7 @@ __attribute__((always_inline)) INLINE static void star_formation_end_density(
  * @param cosmo The current cosmological model.
  */
 __attribute__((always_inline)) INLINE static void
-star_formation_part_has_no_neighbours(struct part* restrict p,
-                                      struct xpart* restrict xp,
+star_formation_part_has_no_neighbours(struct part* p, struct xpart* xp,
                                       const struct star_formation* cd,
                                       const struct cosmology* cosmo) {}
 
@@ -250,7 +252,7 @@ star_formation_part_has_no_neighbours(struct part* restrict p,
  * @param p Pointer to the particle data.
  */
 __attribute__((always_inline)) INLINE static void star_formation_init_part(
-    struct part* restrict p, const struct star_formation* data) {}
+    struct part* p, const struct star_formation* data) {}
 
 /**
  * @brief Sets the star_formation properties of the (x-)particles to a valid
@@ -266,12 +268,11 @@ __attribute__((always_inline)) INLINE static void star_formation_init_part(
  * @param xp Pointer to the extended particle data.
  */
 __attribute__((always_inline)) INLINE static void
-star_formation_first_init_part(const struct phys_const* restrict phys_const,
-                               const struct unit_system* restrict us,
-                               const struct cosmology* restrict cosmo,
+star_formation_first_init_part(const struct phys_const* phys_const,
+                               const struct unit_system* us,
+                               const struct cosmology* cosmo,
                                const struct star_formation* data,
-                               const struct part* restrict p,
-                               struct xpart* restrict xp) {
+                               const struct part* p, struct xpart* xp) {
 
   xp->sf_data.convert_to_star = 0;
 }
diff --git a/src/star_formation/none/star_formation.h b/src/star_formation/none/star_formation.h
index e88dfefe4c65c7258e90aa0dd812a49a45236b18..57544c837368cb0c4290f2dbdbaedebd369d971d 100644
--- a/src/star_formation/none/star_formation.h
+++ b/src/star_formation/none/star_formation.h
@@ -55,13 +55,11 @@ struct star_formation {};
  *
  */
 INLINE static int star_formation_is_star_forming(
-    const struct part* restrict p, const struct xpart* restrict xp,
+    const struct part* p, const struct xpart* xp,
     const struct star_formation* starform, const struct phys_const* phys_const,
-    const struct cosmology* cosmo,
-    const struct hydro_props* restrict hydro_props,
-    const struct unit_system* restrict us,
-    const struct cooling_function_data* restrict cooling,
-    const struct entropy_floor_properties* restrict entropy_floor) {
+    const struct cosmology* cosmo, const struct hydro_props* hydro_props,
+    const struct unit_system* us, const struct cooling_function_data* cooling,
+    const struct entropy_floor_properties* entropy_floor) {
 
   return 0;
 }
@@ -81,7 +79,7 @@ INLINE static int star_formation_is_star_forming(
  * @param dt_star The time-step of this particle.
  */
 INLINE static void star_formation_compute_SFR(
-    const struct part* restrict p, struct xpart* restrict xp,
+    const struct part* p, struct xpart* xp,
     const struct star_formation* starform, const struct phys_const* phys_const,
     const struct hydro_props* hydro_props, const struct cosmology* cosmo,
     const double dt_star) {}
@@ -157,10 +155,8 @@ INLINE static void star_formation_copy_properties(
     const struct part* p, const struct xpart* xp, struct spart* sp,
     const struct engine* e, const struct star_formation* starform,
     const struct cosmology* cosmo, const int with_cosmology,
-    const struct phys_const* phys_const,
-    const struct hydro_props* restrict hydro_props,
-    const struct unit_system* restrict us,
-    const struct cooling_function_data* restrict cooling,
+    const struct phys_const* phys_const, const struct hydro_props* hydro_props,
+    const struct unit_system* us, const struct cooling_function_data* cooling,
     const int convert_part) {}
 
 /**
@@ -169,12 +165,17 @@ INLINE static void star_formation_copy_properties(
  * @param parameter_file The parsed parameter file
  * @param phys_const Physical constants in internal units
  * @param us The current internal system of units
+ * @param hydro_props The propertis of the hydro model.
+ * @param cosmo The current cosmological model.
+ * @param entropy_floor The properties of the entropy floor used in this
+ * simulation.
  * @param starform the star formation law properties to initialize
- *
  */
 INLINE static void starformation_init_backend(
     struct swift_params* parameter_file, const struct phys_const* phys_const,
     const struct unit_system* us, const struct hydro_props* hydro_props,
+    const struct cosmology* cosmo,
+    const struct entropy_floor_properties* entropy_floor,
     const struct star_formation* starform) {}
 
 /**
@@ -199,8 +200,8 @@ INLINE static void starformation_print_backend(
  * @param cosmo The current cosmological model.
  */
 __attribute__((always_inline)) INLINE static void star_formation_end_density(
-    struct part* restrict p, struct xpart* restrict xp,
-    const struct star_formation* cd, const struct cosmology* cosmo) {}
+    struct part* p, struct xpart* xp, const struct star_formation* cd,
+    const struct cosmology* cosmo) {}
 
 /**
  * @brief Sets all particle fields to sensible values when the #part has 0 ngbs.
@@ -213,8 +214,7 @@ __attribute__((always_inline)) INLINE static void star_formation_end_density(
  * @param cosmo The current cosmological model.
  */
 __attribute__((always_inline)) INLINE static void
-star_formation_part_has_no_neighbours(struct part* restrict p,
-                                      struct xpart* restrict xp,
+star_formation_part_has_no_neighbours(struct part* p, struct xpart* xp,
                                       const struct star_formation* cd,
                                       const struct cosmology* cosmo) {}
 
@@ -228,7 +228,7 @@ star_formation_part_has_no_neighbours(struct part* restrict p,
  * @param p Pointer to the particle data.
  */
 __attribute__((always_inline)) INLINE static void star_formation_init_part(
-    struct part* restrict p, const struct star_formation* data) {}
+    struct part* p, const struct star_formation* data) {}
 
 /**
  * @brief Sets the star_formation properties of the (x-)particles to a valid
@@ -244,12 +244,11 @@ __attribute__((always_inline)) INLINE static void star_formation_init_part(
  * @param xp Pointer to the extended particle data.
  */
 __attribute__((always_inline)) INLINE static void
-star_formation_first_init_part(const struct phys_const* restrict phys_const,
-                               const struct unit_system* restrict us,
-                               const struct cosmology* restrict cosmo,
+star_formation_first_init_part(const struct phys_const* phys_const,
+                               const struct unit_system* us,
+                               const struct cosmology* cosmo,
                                const struct star_formation* data,
-                               const struct part* restrict p,
-                               struct xpart* restrict xp) {}
+                               const struct part* p, struct xpart* xp) {}
 
 /**
  * @brief Split the star formation content of a particle into n pieces
diff --git a/src/statistics.c b/src/statistics.c
index 711782a84f45610215f2a6fbfa50af2fe2c57053..edaaf94002115442f2bfafe46a0b95d8e93a7641 100644
--- a/src/statistics.c
+++ b/src/statistics.c
@@ -209,7 +209,7 @@ void stats_collect_part_mapper(void *map_data, int nr_parts, void *extra_data) {
  * @brief The #threadpool mapper function used to collect statistics for #spart.
  *
  * @param map_data Pointer to the particles.
- * @param nr_parts The number of particles in this chunk
+ * @param nr_sparts The number of particles in this chunk
  * @param extra_data The #statistics aggregator.
  */
 void stats_collect_spart_mapper(void *map_data, int nr_sparts,
@@ -298,7 +298,7 @@ void stats_collect_spart_mapper(void *map_data, int nr_sparts,
  * @brief The #threadpool mapper function used to collect statistics for #sink.
  *
  * @param map_data Pointer to the particles.
- * @param nr_parts The number of particles in this chunk
+ * @param nr_sinks The number of particles in this chunk
  * @param extra_data The #statistics aggregator.
  */
 void stats_collect_sink_mapper(void *map_data, int nr_sinks, void *extra_data) {
@@ -383,7 +383,7 @@ void stats_collect_sink_mapper(void *map_data, int nr_sinks, void *extra_data) {
  * @brief The #threadpool mapper function used to collect statistics for #bpart.
  *
  * @param map_data Pointer to the particles.
- * @param nr_parts The number of particles in this chunk
+ * @param nr_bparts The number of particles in this chunk
  * @param extra_data The #statistics aggregator.
  */
 void stats_collect_bpart_mapper(void *map_data, int nr_bparts,