diff --git a/post_processing/reference_tools.py b/post_processing/reference_tools.py
index 8fd34839..7c192b10 100644
--- a/post_processing/reference_tools.py
+++ b/post_processing/reference_tools.py
@@ -13,8 +13,20 @@ class equation_coefficients:
                  'eta':7, 'd_ln_rho':8, 'd2_ln_rho':9, 'd_ln_T':10, 'd_ln_nu':11, 'd_ln_kappa':12,
                  'd_ln_eta':13, 'ds_dr':14}
 
-    def __init__(self,radius=[], file=None):
-        if (len(radius) != 0):
+    def __init__(self, radius=[], file=None, override_nconst=False):
+        """Reads/writes equation_coefficient and custom reference state files.
+
+        Keyword arguments:
+        radius:          If provided, a new state will be created with these
+                         radii and zero coefficients.
+        file:            If provided, the reference state will be read from the
+                         provided file.
+        override_nconst: Fixes file reading for equation_coefficients created
+                         between versions that already had 11 constants but
+                         had not updated the version number to 2 yet.
+                         See https://github.com/geodynamics/Rayleigh/issues/560
+                         for details.
+        """
         if len(radius) != 0:
             if file is not None:
                 raise RuntimeError("Cannot provide radius and file at the same time.")
@@ -26,8 +38,8 @@ def __init__(self,radius=[], file=None):
             self.constants = numpy.zeros(self.nconst     , dtype='float64' )
             self.cset      = numpy.zeros(self.nconst     , dtype='int32'   )
             self.fset      = numpy.zeros(self.nfunc      , dtype='int32'   )
-            self.read(filename=file)
         elif file is not None:
+            self.read(filename=file, override_nconst=override_nconst)
 
     def __getattr__(self, name):
         if name in self.f_dict:
@@ -82,7 +94,7 @@ def write(self, filename='ecoefs.dat'):
         self.functions.tofile(fd)
         fd.close() 
         
-    def read(self, filename='equation_coefficients'):
+    def read(self, filename='equation_coefficients', override_nconst=False):
         class_version = self.version
         fd = open(filename,'rb')
         picheck = numpy.fromfile(fd,dtype='int32',count=1)[0]
@@ -91,8 +103,14 @@ def read(self, filename='equation_coefficients'):
         if self.version > 1:
             self.nconst = numpy.fromfile(fd,dtype='int32', count=1)[0]
             self.nfunc = numpy.fromfile(fd,dtype='int32', count=1)[0]
-        else: # if the version is 1, nconst was 10
-            self.nconst = 10
+        else:
+            if override_nconst:
+                # Override for the versions of Rayleigh that already had 11 constants
+                # but had not updated the version number of the file yet.
+                self.nconst = 11
+            else:
+                # if the version is 1, nconst was 10
+                self.nconst = 10
         self.cset = numpy.fromfile(fd,dtype='int32',count=self.nconst)
         self.fset = numpy.fromfile(fd,dtype='int32',count=self.nfunc)
         self.constants = numpy.fromfile(fd,dtype='float64',count=self.nconst)