diff --git a/PRIVATE b/PRIVATE
index 17250a3a2..9c1f83bab 160000
--- a/PRIVATE
+++ b/PRIVATE
@@ -1 +1 @@
-Subproject commit 17250a3a29f07af3c4d4e4083213d46efd534268
+Subproject commit 9c1f83babb7894bfaa16255d6c15a4a438c7f168
diff --git a/src/grid/DAMASK_grid.f90 b/src/grid/DAMASK_grid.f90
index 29eded0a9..39870a8ff 100644
--- a/src/grid/DAMASK_grid.f90
+++ b/src/grid/DAMASK_grid.f90
@@ -107,8 +107,6 @@ program DAMASK_grid
     load_steps, &
     load_step, &
     solver, &
-    initial_conditions, &
-    thermal, &
     step_bc, &
     step_mech, &
     step_discretization
@@ -303,9 +301,7 @@ program DAMASK_grid
     select case (ID(field))
 
       case (FIELD_THERMAL_ID)
-        initial_conditions => config_load%get('initial_conditions',defaultVal=emptyDict)
-        thermal            => initial_conditions%get('thermal',defaultVal=emptyDict)
-        call grid_thermal_spectral_init(thermal%get_asFloat('T'))
+        call grid_thermal_spectral_init()
 
       case (FIELD_DAMAGE_ID)
         call grid_damage_spectral_init()
diff --git a/src/grid/discretization_grid.f90 b/src/grid/discretization_grid.f90
index 57f3b9ac0..3cd779ec0 100644
--- a/src/grid/discretization_grid.f90
+++ b/src/grid/discretization_grid.f90
@@ -37,13 +37,14 @@ module discretization_grid
     size3offset                                                                                     !< (local) size offset in 3rd direction
 
   public :: &
-    discretization_grid_init
+    discretization_grid_init, &
+    discretization_grid_getInitialCondition
 
 contains
 
 
 !--------------------------------------------------------------------------------------------------
-!> @brief reads the geometry file to obtain information on discretization
+!> @brief Read the geometry file to obtain information on discretization.
 !--------------------------------------------------------------------------------------------------
 subroutine discretization_grid_init(restart)
 
@@ -171,7 +172,7 @@ end subroutine discretization_grid_init
 
 
 !---------------------------------------------------------------------------------------------------
-!> @brief Calculate undeformed position of IPs/cell centers (pretend to be an element)
+!> @brief Calculate undeformed position of IPs/cell centers (pretend to be an element).
 !---------------------------------------------------------------------------------------------------
 function IPcoordinates0(cells,geomSize,cells3Offset)
 
@@ -196,7 +197,7 @@ end function IPcoordinates0
 
 
 !---------------------------------------------------------------------------------------------------
-!> @brief Calculate position of undeformed nodes (pretend to be an element)
+!> @brief Calculate position of undeformed nodes (pretend to be an element).
 !---------------------------------------------------------------------------------------------------
 pure function nodes0(cells,geomSize,cells3Offset)
 
@@ -220,7 +221,7 @@ end function nodes0
 
 
 !--------------------------------------------------------------------------------------------------
-!> @brief Calculate IP interface areas
+!> @brief Calculate IP interface areas.
 !--------------------------------------------------------------------------------------------------
 pure function cellSurfaceArea(geomSize,cells)
 
@@ -238,7 +239,7 @@ end function cellSurfaceArea
 
 
 !--------------------------------------------------------------------------------------------------
-!> @brief Calculate IP interface areas normals
+!> @brief Calculate IP interface areas normals.
 !--------------------------------------------------------------------------------------------------
 pure function cellSurfaceNormal(nElems)
 
@@ -257,7 +258,7 @@ end function cellSurfaceNormal
 
 
 !--------------------------------------------------------------------------------------------------
-!> @brief Build IP neighborhood relations
+!> @brief Build IP neighborhood relations.
 !--------------------------------------------------------------------------------------------------
 pure function IPneighborhood(cells)
 
@@ -313,4 +314,40 @@ pure function IPneighborhood(cells)
 end function IPneighborhood
 
 
+!--------------------------------------------------------------------------------------------------
+!> @brief Read initial condition from VTI file.
+!--------------------------------------------------------------------------------------------------
+function discretization_grid_getInitialCondition(label) result(ic)
+
+  character(len=*), intent(in) :: label
+  real(pReal), dimension(cells(1),cells(2),cells3) :: ic
+
+  real(pReal), dimension(:), allocatable :: ic_global, ic_local
+  integer(MPI_INTEGER_KIND) :: err_MPI
+
+  integer, dimension(worldsize) :: &
+    displs, sendcounts
+
+  if (worldrank == 0) then
+    ic_global = VTI_readDataset_real(IO_read(interface_geomFile),label)
+  else
+    allocate(ic_global(0))                                                                          ! needed for IntelMPI
+  endif
+
+  call MPI_Gather(product(cells(1:2))*cells3Offset, 1_MPI_INTEGER_KIND,MPI_INTEGER,displs,&
+                  1_MPI_INTEGER_KIND,MPI_INTEGER,0_MPI_INTEGER_KIND,MPI_COMM_WORLD,err_MPI)
+  if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error'
+  call MPI_Gather(product(cells(1:2))*cells3,      1_MPI_INTEGER_KIND,MPI_INTEGER,sendcounts,&
+                  1_MPI_INTEGER_KIND,MPI_INTEGER,0_MPI_INTEGER_KIND,MPI_COMM_WORLD,err_MPI)
+  if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error'
+
+  allocate(ic_local(product(cells(1:2))*cells3))
+  call MPI_Scatterv(ic_global,sendcounts,displs,MPI_DOUBLE,ic_local,size(ic_local),&
+                    MPI_DOUBLE,0_MPI_INTEGER_KIND,MPI_COMM_WORLD,err_MPI)
+  if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error'
+
+  ic = reshape(ic_local,[cells(1),cells(2),cells3])
+
+end function discretization_grid_getInitialCondition
+
 end module discretization_grid
diff --git a/src/grid/grid_thermal_spectral.f90 b/src/grid/grid_thermal_spectral.f90
index 1066a1f61..6bbeabf00 100644
--- a/src/grid/grid_thermal_spectral.f90
+++ b/src/grid/grid_thermal_spectral.f90
@@ -65,9 +65,7 @@ contains
 !--------------------------------------------------------------------------------------------------
 !> @brief allocates all neccessary fields and fills them with data
 !--------------------------------------------------------------------------------------------------
-subroutine grid_thermal_spectral_init(T_0)
-
-  real(pReal), intent(in) :: T_0
+subroutine grid_thermal_spectral_init()
 
   PetscInt, dimension(0:worldsize-1) :: localK
   integer :: i, j, k, ce
@@ -105,9 +103,9 @@ subroutine grid_thermal_spectral_init(T_0)
 
 !--------------------------------------------------------------------------------------------------
 ! init fields
-  allocate(T_current(cells(1),cells(2),cells3), source=T_0)
-  allocate(T_lastInc(cells(1),cells(2),cells3), source=T_0)
-  allocate(T_stagInc(cells(1),cells(2),cells3), source=T_0)
+  T_current = discretization_grid_getInitialCondition('T')
+  T_lastInc = T_current
+  T_stagInc = T_current
 
 !--------------------------------------------------------------------------------------------------
 ! initialize solver specific parts of PETSc