diff --git a/Template-Config.sh b/Template-Config.sh
index d2f1be5102e9cead3b0f79348bed9c4d45e7f1c8..0a37a2e777a36c7bd649f64113f6b1b8dd231e7a 100644
--- a/Template-Config.sh
+++ b/Template-Config.sh
@@ -109,6 +109,7 @@ DOUBLEPRECISION=1                             # if activated and set to 1, use d
 
 #---------------------------------------- Output/Input options
 
+INITIAL_CONDITIONS_CONTAIN_ENTROPY
 #OUTPUT_VELOCITY_GRADIENT                     # output velocity gradients
 #OUTPUT_PRESSURE                              # output gas pressure   
 #OUTPUT_ENTROPY                               # output gas entropy
diff --git a/documentation/04_config-options.md b/documentation/04_config-options.md
index 8e42f0a8f3f3360fe11ae912c60c0469fe00d677..243fe57e2c9c901e91b020f3c3e1e13387b499a0 100644
--- a/documentation/04_config-options.md
+++ b/documentation/04_config-options.md
@@ -598,6 +598,12 @@ formulation. This is only useful if `PRESSURE_ENTROPY_SPH` is used.
 
 -------
 
+**INITIAL_CONDITIONS_CONTAIN_ENTROPY**
+
+The intial conditions file contains entropy instead of the thermal energy.
+
+------
+
 **GAMMA** = 1.4
 
 Sets the equation of state index in the ideal gas law that is normally
@@ -871,7 +877,6 @@ you need to increase `NUMBER_OF_MPI_LISTENERS_PER_NODE`.
 Output/Input options                                       {#io}
 ====================
 
-
 **POWERSPEC_ON_OUTPUT**
 
 Creates a power spectrum measurement for every output time, i.e. for
diff --git a/src/data/allvars.h b/src/data/allvars.h
index 435be921ebcc574ac3cf1e621fd87aae3e16d254..65fbc2492735281ff8d48007728424f51ba6ee60 100644
--- a/src/data/allvars.h
+++ b/src/data/allvars.h
@@ -97,7 +97,6 @@ struct global_data_all_processes : public parameters
   double InitGasU;         /**< the same, but converted to thermal energy per unit mass */
   double MinEgySpec;       /**< the minimum allowed temperature expressed as energy per unit mass */
 
-  int FlagICsContainedEntropy;
 
   /* some force counters  */
 
diff --git a/src/io/snap_io.cc b/src/io/snap_io.cc
index cae63b05609712c44a132b7680bc409ede45684f..be030887546353e72acaab7d408985b39c7f4de7 100644
--- a/src/io/snap_io.cc
+++ b/src/io/snap_io.cc
@@ -415,11 +415,14 @@ void snap_io::read_ic(const char *fname)
   Sp->NumPart += add_numpart;
 #endif
 
-  All.FlagICsContainedEntropy = 0;
+
 
 #ifdef GADGET2_HEADER
-  if(header.flag_entropy_instead_u)
-    All.FlagICsContainedEntropy = 1;
+#ifndef INITIAL_CONDITIONS_CONTAIN_ENTROPY
+  if(header.flag_entropy_instead_u) Terminate("Initial condition file contains entropy, but INITIAL_CONDITIONS_CONTAIN_ENTROPY is not set\n");
+#else    
+  if(! header.flag_entropy_instead_u)Terminate("Initial condition file contains uthermal, but INITIAL_CONDITIONS_CONTAIN_ENTROPY is set\n");
+#endif
 #endif
 
   TIMER_STOP(CPU_SNAPSHOT);
diff --git a/src/main/init.cc b/src/main/init.cc
index addf4ea8f5f2816d25e096fd538686533e28fff5..198f9fef983e9733052e6b210de04adca52ddeec 100644
--- a/src/main/init.cc
+++ b/src/main/init.cc
@@ -129,7 +129,6 @@ void sim::init(int RestartSnapNum)
 
       Mem.myfree(tmp);
 
-      All.FlagICsContainedEntropy = 0;
 
       int count = 0;
       for(int i = 0; i < Sp.NumPart; i++)
@@ -477,15 +476,16 @@ void sim::init(int RestartSnapNum)
    * Once the density has been computed, we can convert to entropy.
    */
 #ifdef PRESSURE_ENTROPY_SPH
-  if(All.FlagICsContainedEntropy == 0)
+#ifndef INITIAL_CONDITIONS_CONTAIN_ENTROPY
     NgbTree.setup_entropy_to_invgamma();
+#endif
 #endif
 
   double mass = 0;
   for(int i = 0; i < Sp.NumGas; i++)
     {
-      if(All.FlagICsContainedEntropy == 0)
-        {
+#ifndef INITIAL_CONDITIONS_CONTAIN_ENTROPY
+        
           if(ThisTask == 0 && i == 0)
             printf("INIT: Converting u -> entropy\n");
 
@@ -493,8 +493,8 @@ void sim::init(int RestartSnapNum)
           Sp.SphP[i].Entropy = GAMMA_MINUS1 * Sp.SphP[i].Entropy / pow(Sp.SphP[i].Density * All.cf_a3inv, GAMMA_MINUS1);
 #endif
           Sp.SphP[i].EntropyPred = Sp.SphP[i].Entropy;
-        }
-
+        
+#endif
       /* The predicted entropy values have been already set for all SPH formulation */
       /* so it should be ok computing pressure and csound now */
       Sp.SphP[i].set_thermodynamic_variables();