diff --git a/src/constitutive.f90 b/src/constitutive.f90
index 36400fca7..e68d90de6 100644
--- a/src/constitutive.f90
+++ b/src/constitutive.f90
@@ -112,13 +112,15 @@ module constitutive
   interface
 
 ! == cleaned:begin =================================================================================
-    module subroutine mech_init
+    module subroutine mech_init(phases)
+      class(tNode), pointer :: phases
     end subroutine mech_init
 
     module subroutine damage_init
     end subroutine damage_init
 
-    module subroutine thermal_init
+    module subroutine thermal_init(phases)
+      class(tNode), pointer :: phases
     end subroutine thermal_init
 
 
@@ -197,10 +199,10 @@ module constitutive
       real(pReal), dimension(3,3) :: P
     end function constitutive_mech_getP
 
-    module function constitutive_thermal_T(co,ip,el) result(T)
-      integer, intent(in) :: co, ip, el
+    module function thermal_T(ph,me) result(T)
+      integer, intent(in) :: ph,me
       real(pReal) :: T
-    end function constitutive_thermal_T
+    end function thermal_T
 
 
     module subroutine constitutive_mech_setF(F,co,ip,el)
@@ -463,6 +465,8 @@ subroutine constitutive_init
     phases
 
 
+  print'(/,a)', ' <<<+-  constitutive init  -+>>>'; flush(IO_STDOUT)
+
   debug_constitutive => config_debug%get('constitutive', defaultVal=emptyList)
   debugConstitutive%basic      =  debug_constitutive%contains('basic')
   debugConstitutive%extensive  =  debug_constitutive%contains('extensive')
@@ -471,15 +475,14 @@ subroutine constitutive_init
   debugConstitutive%ip         =  config_debug%get_asInt('integrationpoint',defaultVal = 1)
   debugConstitutive%grain      =  config_debug%get_asInt('grain',defaultVal = 1)
 
-!--------------------------------------------------------------------------------------------------
-! initialize constitutive laws
-  print'(/,a)', ' <<<+-  constitutive init  -+>>>'; flush(IO_STDOUT)
-  call mech_init
-  call damage_init
-  call thermal_init
-
 
   phases => config_material%get('phase')
+
+  call mech_init(phases)
+  call damage_init
+  call thermal_init(phases)
+
+
   constitutive_source_maxSizeDotState = 0
   PhaseLoop2:do ph = 1,phases%length
 !--------------------------------------------------------------------------------------------------
diff --git a/src/constitutive_mech.f90 b/src/constitutive_mech.f90
index f5a5cd0a2..97bf7d853 100644
--- a/src/constitutive_mech.f90
+++ b/src/constitutive_mech.f90
@@ -319,7 +319,10 @@ contains
 !> @brief Initialize mechanical field related constitutive models
 !> @details Initialize elasticity, plasticity and stiffness degradation models.
 !--------------------------------------------------------------------------------------------------
-module subroutine mech_init
+module subroutine mech_init(phases)
+
+  class(tNode), pointer :: &
+    phases
 
   integer :: &
     el, &
@@ -331,7 +334,6 @@ module subroutine mech_init
     Nconstituents
   class(tNode), pointer :: &
     num_crystallite, &
-    phases, &
     phase, &
     mech, &
     elastic, &
@@ -341,7 +343,6 @@ module subroutine mech_init
 
 !-------------------------------------------------------------------------------------------------
 ! initialize elasticity (hooke)                         !ToDO: Maybe move to elastic submodule along with function homogenizedC?
-  phases => config_material%get('phase')
   allocate(phase_elasticity(phases%length), source = ELASTICITY_undefined_ID)
   allocate(phase_elasticityInstance(phases%length), source = 0)
   allocate(phase_NstiffnessDegradations(phases%length),source=0)
diff --git a/src/constitutive_thermal.f90 b/src/constitutive_thermal.f90
index 01d517124..bdcd0bc26 100644
--- a/src/constitutive_thermal.f90
+++ b/src/constitutive_thermal.f90
@@ -2,6 +2,12 @@
 !> @brief internal microstructure state for all thermal sources and kinematics constitutive models
 !----------------------------------------------------------------------------------------------------
 submodule(constitutive) constitutive_thermal
+  
+  type :: tDataContainer
+    real(pReal), dimension(:), allocatable :: T
+  end type tDataContainer
+ 
+  type(tDataContainer), dimension(:), allocatable :: current
 
   interface
 
@@ -49,8 +55,29 @@ contains
 !----------------------------------------------------------------------------------------------
 !< @brief initializes thermal sources and kinematics mechanism
 !----------------------------------------------------------------------------------------------
-module subroutine thermal_init
+module subroutine thermal_init(phases)
+  
+  class(tNode), pointer :: &
+    phases
+  
+  integer :: &
+    ph, &
+    Nconstituents
 
+
+  print'(/,a)', ' <<<+-  constitutive_mech init  -+>>>'
+
+  allocate(current(phases%length))
+
+
+  do ph = 1, phases%length
+    
+    Nconstituents = count(material_phaseAt == ph) * discretization_nIPs
+
+    allocate(current(ph)%T(Nconstituents))
+
+  enddo
+ 
 ! initialize source mechanisms
   if(maxval(phase_Nsources) /= 0) then
     where(source_thermal_dissipation_init (maxval(phase_Nsources))) phase_source = SOURCE_thermal_dissipation_ID
@@ -122,21 +149,16 @@ end subroutine constitutive_thermal_getRateAndItsTangents
 
 
 
-
 ! getter for non-thermal (e.g. mech)
-module function constitutive_thermal_T(co,ip,el) result(T)
+module function thermal_T(ph,me) result(T)
 
-  integer, intent(in) :: co, ip, el
+  integer, intent(in) :: ph, me
   real(pReal) :: T
 
-  integer :: ho, tme
 
-  ho = material_homogenizationAt(el)
-  tme = material_homogenizationMemberAt(ip,el)
+  T = current(ph)%T(me)
 
-  T = temperature(ho)%p(tme)
-
-end function constitutive_thermal_T
+end function thermal_T
 
 
 ! setter for homogenization