diff --git a/PRIVATE b/PRIVATE
index 7cae6f6fb..e5c897781 160000
--- a/PRIVATE
+++ b/PRIVATE
@@ -1 +1 @@
-Subproject commit 7cae6f6fb6daca513fefbce21e4f07ab36fd3c78
+Subproject commit e5c8977814ea8afca021fad2181dbb2a986bea70
diff --git a/src/phase_mechanical_eigen_thermalexpansion.f90 b/src/phase_mechanical_eigen_thermalexpansion.f90
index b62db88ef..d03ea0eb2 100644
--- a/src/phase_mechanical_eigen_thermalexpansion.f90
+++ b/src/phase_mechanical_eigen_thermalexpansion.f90
@@ -98,7 +98,7 @@ module subroutine thermalexpansion_LiAndItsTangent(Li, dLi_dTstar, ph,me)
 
   associate(prm => param(kinematics_thermal_expansion_instance(ph)))
     Li = dot_T * ( &
-                  prm%A(1:3,1:3,1)*(T - prm%T_ref)**0 &                                             ! constant  coefficient
+                  prm%A(1:3,1:3,1)*1.0_pReal &                                                      ! constant  coefficient
                 + prm%A(1:3,1:3,2)*(T - prm%T_ref)**1 &                                             ! linear    coefficient
                 + prm%A(1:3,1:3,3)*(T - prm%T_ref)**2 &                                             ! quadratic coefficient
                 ) / &
diff --git a/src/phase_mechanical_elastic.f90 b/src/phase_mechanical_elastic.f90
index e2999d96f..9e5c2653c 100644
--- a/src/phase_mechanical_elastic.f90
+++ b/src/phase_mechanical_elastic.f90
@@ -8,6 +8,8 @@ submodule(phase:mechanical) elastic
       C_33 = 0.0_pReal, &
       C_44 = 0.0_pReal, &
       C_66 = 0.0_pReal
+    real(pReal) :: &
+      T_ref
   end type tParameters
 
   type(tParameters), allocatable, dimension(:) :: param
@@ -45,6 +47,8 @@ module subroutine elastic_init(phases)
 
     associate(prm => param(ph))
 
+      prm%T_ref = elastic%get_asFloat('T_ref', defaultVal=T_ROOM)
+
       prm%C_11(1) = elastic%get_asFloat('C_11')
       prm%C_11(2) = elastic%get_asFloat('C_11,T',  defaultVal=0.0_pReal)
       prm%C_11(3) = elastic%get_asFloat('C_11,T^2',defaultVal=0.0_pReal)
@@ -93,16 +97,36 @@ module function elastic_C66(ph,en) result(C66)
   associate(prm => param(ph))
     C66 = 0.0_pReal
     T = thermal_T(ph,en)
-    C66(1,1) = prm%C_11(1) + prm%C_11(2)*T + prm%C_11(3)*T**2
-    C66(1,2) = prm%C_12(1) + prm%C_12(2)*T + prm%C_12(3)*T**2
-    C66(4,4) = prm%C_44(1) + prm%C_44(2)*T + prm%C_44(3)*T**2
+
+    C66(1,1) =   prm%C_11(1)* 1.0_pReal &
+               + prm%C_11(2)*(T - prm%T_ref)**1 &
+               + prm%C_11(3)*(T - prm%T_ref)**2
+
+    C66(1,2) =   prm%C_12(1)* 1.0_pReal &
+               + prm%C_12(2)*(T - prm%T_ref)**1 &
+               + prm%C_12(3)*(T - prm%T_ref)**2
+
+    C66(4,4) =   prm%C_44(1)*1.0_pReal &
+               + prm%C_44(2)*(T - prm%T_ref)**1 &
+               + prm%C_44(3)*(T - prm%T_ref)**2
+
 
     if (any(phase_lattice(ph) == ['hP','tI'])) then
-      C66(1,3) = prm%C_13(1) + prm%C_13(2)*T + prm%C_13(3)*T**2
-      C66(3,3) = prm%C_33(1) + prm%C_33(2)*T + prm%C_33(3)*T**2
+      C66(1,3) =   prm%C_13(1)*1.0_pReal &
+                 + prm%C_13(2)*(T - prm%T_ref)**1 &
+                 + prm%C_13(3)*(T - prm%T_ref)**2
+
+      C66(3,3) =   prm%C_33(1)*1.0_pReal &
+                 + prm%C_33(2)*(T - prm%T_ref)**1 &
+                 + prm%C_33(3)*(T - prm%T_ref)**2
+
     end if
 
-    if (phase_lattice(ph) == 'tI') C66(6,6) = prm%C_66(1) + prm%C_66(2)*T + prm%C_66(3)*T**2
+    if (phase_lattice(ph) == 'tI') then
+      C66(6,6) =   prm%C_66(1)*1.0_pReal &
+                 + prm%C_66(2)*(T - prm%T_ref)**1 &
+                 + prm%C_66(3)*(T - prm%T_ref)**2
+    end if
 
     C66 = lattice_symmetrize_C66(C66,phase_lattice(ph))