diff --git a/post_processing/reference_tools.py b/post_processing/reference_tools.py
index fb3dbdaf..e1e90495 100644
--- a/post_processing/reference_tools.py
+++ b/post_processing/reference_tools.py
@@ -3,12 +3,12 @@
 
 class equation_coefficients:
     """ equation coeff class  """
-    nconst = 10
+    nconst = 11 # new "version 2" default
     nfunc = 14
-    version=1
+    version=2
     c_dict    = {'two_omega':1, 'buoy_fact':2, 'p_fact':3, 'lorentz_fact':4, 'visc_fact':5,
                  'diff_fact':6, 'resist_fact':7, 'nu_heat_fact':8, 'ohm_heat_fact':9,
-                 'luminosity':10}
+                 'luminosity':10, 'dsdr_scale': 11}
     f_dict    = {'density':1, 'buoy':2, 'nu':3, 'temperature':4, 'kappa':5, 'heating':6,
                  '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}
@@ -65,9 +65,13 @@ def write(self, filename='ecoefs.dat'):
         pi = numpy.array([314],dtype='int32')
         nr = numpy.array([self.nr],dtype='int32')
         version = numpy.array([self.version],dtype='int32')
+        nconst = numpy.array([self.nconst],dtype='int32')
+        nfunc = numpy.array([self.nfunc],dtype='int32')
         fd = open(filename,'wb')
         pi.tofile(fd)
         version.tofile(fd)
+        nconst.tofile(fd)
+        nfunc.tofile(fd)
         self.cset.tofile(fd)
         self.fset.tofile(fd)
         self.constants.tofile(fd)
@@ -77,20 +81,33 @@ def write(self, filename='ecoefs.dat'):
         fd.close() 
         
     def read(self, filename='equation_coefficients'):
-
+        class_version = self.version
         fd = open(filename,'rb')
         picheck = numpy.fromfile(fd,dtype='int32',count=1)[0]
        
         self.version = numpy.fromfile(fd,dtype='int32', count=1)[0]
+        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
         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)
+        self.constants = numpy.fromfile(fd,dtype='float64',count=self.nconst)        
         self.nr = numpy.fromfile(fd,dtype='int32',count=1)[0]
         self.radius = numpy.fromfile(fd,dtype='float64',count=self.nr)
         functions=numpy.fromfile(fd,dtype='float64',count=self.nr*self.nfunc)
         self.functions = numpy.reshape(functions, (self.nfunc,self.nr))
         fd.close()
 
+        # now if we read in v. 1, convert current instance to v. 2
+        if self.version == 1:
+            self.nconst = 11
+            self.version = class_version
+            self.constants[10] = 1.
+            print("Version of input file was 1.")
+            print("Converting current equation_coefficients instance's version to %i." %class_version)
+
 class background_state:
     nr = None
     radius = None
diff --git a/src/Physics/PDE_Coefficients.F90 b/src/Physics/PDE_Coefficients.F90
index 4153bdfe..8e4bf992 100644
--- a/src/Physics/PDE_Coefficients.F90
+++ b/src/Physics/PDE_Coefficients.F90
@@ -75,7 +75,7 @@ Module PDE_Coefficients
 
     ! Version number for the "equation_coefficients" container that is output to the simulation directory
     ! (i.e., the human-obtainable version of "ref")
-    Integer, Parameter  :: eqn_coeff_version = 1
+    Integer, Parameter  :: eqn_coeff_version = 2
 
     ! Which background state to use; default 1 (non-dimensional Boussinesq)
     Integer :: reference_type = 1 
@@ -1224,6 +1224,8 @@ Subroutine Write_Equation_Coefficients_File(filename)
             Write(15) eqn_coeff_version
             ra_constant_set(:) = 1
             ra_function_set(:) = 1
+            Write(15) n_ra_constants
+            Write(15) n_ra_functions
             Write(15) (ra_constant_set(i), i=1,n_ra_constants)
             Write(15) (ra_function_set(i), i=1,n_ra_functions)
             Write(15) (ra_constants(i), i=1,n_ra_constants)
@@ -1242,7 +1244,7 @@ Subroutine Read_Custom_Reference_File(filename)
         Character*120, Intent(In) :: filename
         Character*120 :: ref_file
         Integer :: pi_integer,nr_ref, eqversion
-        Integer :: i, k, j, n_scalars
+        Integer :: i, k, j, n_scalars, dummy
         Integer :: cset(1:n_ra_constants), fset(1:n_ra_functions)
         Real*8  :: input_constants(1:n_ra_constants)
         Real*8, Allocatable :: ref_arr_old(:,:), rtmp(:), rtmp2(:)
@@ -1290,6 +1292,9 @@ Subroutine Read_Custom_Reference_File(filename)
                 cset(11) = 1 ! treat this as if c_11 = 1 was specified in the custom reference file
                 input_constants(11) = 1.0d0 
             Else
+                Read(15) dummy ! n_ra_constants, but its up to the user to make sure this is the same 
+                               ! same number specified in Rayleigh via "n_active_scalars"
+                Read(15) dummy ! n_ra_functions, which Rayleigh should already know           
                 Read(15) cset(1:n_ra_constants)
                 Read(15) fset(1:n_ra_functions) 
                 Read(15) input_constants(1:n_ra_constants)