diff --git a/src/Physics/PDE_Coefficients.F90 b/src/Physics/PDE_Coefficients.F90
index eb1e102c..9a778848 100755
--- a/src/Physics/PDE_Coefficients.F90
+++ b/src/Physics/PDE_Coefficients.F90
@@ -689,9 +689,10 @@ Subroutine Augment_Reference()
 
         If (my_rank .eq. 0) Then
             Call stdout%print('Reference state will be augmented.')
-            Call stdout%print('Only heating and buoyancy may be modified.')
+            Call stdout%print('Only heating, buoyancy, or background dTdr/dSdr may be modified.')
             Call stdout%print('Heating requires both c_10 and f_6 to be set.')
             Call stdout%print('Buoyancy requires both c_2 and f_2 to be set.')
+            Call stdout%print('dTdr/dSdr requires f_14 to be set.')            
             Call stdout%print('Reading from: '//Trim(custom_reference_file))
         Endif
 
@@ -703,11 +704,10 @@ Subroutine Augment_Reference()
         Allocate(temp_functions(1:n_r, 1:n_ra_functions))
         Allocate(temp_constants(1:n_ra_constants))
         temp_functions(:,:) = ra_functions(:,:)
-
+        
         ! Note that ra_constants is allocated up to max_ra_constants
         temp_constants(:) = ra_constants(1:n_ra_constants)
 
-
         Call Read_Custom_Reference_File(custom_reference_file)
 
         If (use_custom_constant(10) .and. use_custom_function(6)) Then
@@ -732,6 +732,16 @@ Subroutine Augment_Reference()
             temp_constants(2) = ra_constants(2)
         Endif
 
+        If (use_custom_function(14)) Then
+            If (my_rank .eq. 0) Then
+                Call stdout%print('Background thermal gradient is set to:')
+                Call stdout%print('f_14')
+                Call stdout%print(' ')
+            Endif        
+            ref%dsdr(:) = ra_functions(:,14)
+            temp_functions(:,14) = ra_functions(:,14)
+        Endif
+
         ra_constants(1:n_ra_constants) = temp_constants(:)
         ra_functions(:,:) = temp_functions(:,:)
         DeAllocate(temp_functions, temp_constants)