diff --git a/src/grid/DAMASK_grid.f90 b/src/grid/DAMASK_grid.f90
index b1da0c2a2..a5f441cc0 100644
--- a/src/grid/DAMASK_grid.f90
+++ b/src/grid/DAMASK_grid.f90
@@ -466,7 +466,14 @@ program DAMASK_grid
         call MPI_Allreduce(interface_SIGUSR2,signal,1,MPI_LOGICAL,MPI_LOR,MPI_COMM_WORLD,ierr)
         if (ierr /= 0) error stop 'MPI error'
         if (mod(inc,loadCases(l)%f_restart) == 0 .or. signal) then
-          call mechanical_restartWrite
+          do field = 1, nActiveFields
+            select case (ID(field))
+              case(FIELD_MECH_ID)
+                call mechanical_restartWrite
+              case(FIELD_THERMAL_ID)
+                call grid_thermal_spectral_restartWrite
+            end select
+          end do
           call CPFEM_restartWrite
         endif
         if (signal) call interface_setSIGUSR2(.false.)
diff --git a/src/grid/grid_thermal_spectral.f90 b/src/grid/grid_thermal_spectral.f90
index 67c7ba1c3..0f539b73b 100644
--- a/src/grid/grid_thermal_spectral.f90
+++ b/src/grid/grid_thermal_spectral.f90
@@ -17,6 +17,9 @@ module grid_thermal_spectral
   use parallelization
   use IO
   use spectral_utilities
+  use DAMASK_interface
+  use HDF5
+  use HDF5_utilities
   use discretization_grid
   use homogenization
   use YAML_types
@@ -55,6 +58,7 @@ module grid_thermal_spectral
   public :: &
     grid_thermal_spectral_init, &
     grid_thermal_spectral_solution, &
+    grid_thermal_spectral_restartWrite, &
     grid_thermal_spectral_forward
 
 contains
@@ -70,7 +74,9 @@ subroutine grid_thermal_spectral_init(T_0)
   PetscInt, dimension(0:worldsize-1) :: localK
   integer :: i, j, k, ce
   DM :: thermal_grid
+  real(pReal), dimension(:), allocatable :: T_restart, T_lastInc_restart
   PetscScalar, dimension(:,:,:), pointer :: T_PETSc
+  integer(HID_T) :: fileHandle, groupHandle
   PetscErrorCode :: ierr
   class(tNode), pointer :: &
     num_grid
@@ -130,17 +136,26 @@ subroutine grid_thermal_spectral_init(T_0)
   xend = xstart + xend - 1
   yend = ystart + yend - 1
   zend = zstart + zend - 1
-  allocate(T_current(grid(1),grid(2),grid3), source=0.0_pReal)
-  allocate(T_lastInc(grid(1),grid(2),grid3), source=0.0_pReal)
-  allocate(T_stagInc(grid(1),grid(2),grid3), source=0.0_pReal)
+  allocate(T_current(grid(1),grid(2),grid3), source=T_0)
+  allocate(T_lastInc(grid(1),grid(2),grid3), source=T_0)
+  allocate(T_stagInc(grid(1),grid(2),grid3), source=T_0)
+
+  restartRead: if (interface_restartInc > 0) then
+    print'(/,1x,a,i0,a)', 'reading restart data of increment ', interface_restartInc, ' from file'
+
+    fileHandle  = HDF5_openFile(getSolverJobName()//'_restart.hdf5','r')
+    groupHandle = HDF5_openGroup(fileHandle,'solver')
+
+    call HDF5_read(T_restart,groupHandle,'T',.false.)
+    T_current = reshape(T_restart,[grid(1),grid(2),grid3])
+    call HDF5_read(T_lastInc_restart,groupHandle,'T_lastInc',.false.)
+    T_lastInc = reshape(T_lastInc_restart,[grid(1),grid(2),grid3])
+  end if restartRead 
 
   ce = 0
   do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1)
     ce = ce + 1
-    T_current(i,j,k) = T_0
-    T_lastInc(i,j,k) = T_current(i,j,k)
-    T_stagInc(i,j,k) = T_current(i,j,k)
-    call homogenization_thermal_setField(T_0,0.0_pReal,ce)
+    call homogenization_thermal_setField(T_current(i,j,k),0.0_pReal,ce)
   end do; end do; end do
 
   call DMDAVecGetArrayF90(thermal_grid,solution_vec,T_PETSc,ierr); CHKERRQ(ierr)
@@ -242,6 +257,34 @@ subroutine grid_thermal_spectral_forward(cutBack)
 end subroutine grid_thermal_spectral_forward
 
 
+!--------------------------------------------------------------------------------------------------
+!> @brief Write current solver and constitutive data for restart to file
+!--------------------------------------------------------------------------------------------------
+subroutine grid_thermal_spectral_restartWrite
+
+  PetscErrorCode :: ierr
+  DM :: dm_local
+  integer(HID_T) :: fileHandle, groupHandle
+  PetscScalar, dimension(:,:,:), pointer :: T
+
+  call SNESGetDM(thermal_snes,dm_local,ierr); CHKERRQ(ierr)
+  call DMDAVecGetArrayF90(dm_local,solution_vec,T,ierr); CHKERRQ(ierr)
+
+  print'(1x,a)', 'writing solver data required for restart to file'; flush(IO_STDOUT)
+
+  fileHandle  = HDF5_openFile(getSolverJobName()//'_restart.hdf5','w')
+  groupHandle = HDF5_addGroup(fileHandle,'solver')
+  call HDF5_write(T,groupHandle,'T')
+  call HDF5_write(T_lastInc,groupHandle,'T_lastInc')
+  call HDF5_closeGroup(groupHandle)
+  call HDF5_closeFile(fileHandle)
+
+  call DMDAVecRestoreArrayF90(dm_local,solution_vec,T,ierr); CHKERRQ(ierr)
+
+end subroutine grid_thermal_spectral_restartWrite
+
+
+
 !--------------------------------------------------------------------------------------------------
 !> @brief forms the spectral thermal residual vector
 !--------------------------------------------------------------------------------------------------
diff --git a/src/phase.f90 b/src/phase.f90
index 6035b4491..9ee3d5a4b 100644
--- a/src/phase.f90
+++ b/src/phase.f90
@@ -123,11 +123,20 @@ module phase
       integer, intent(in) :: ph
     end subroutine mechanical_restartWrite
 
+    module subroutine thermal_restartWrite(groupHandle,ph)
+      integer(HID_T), intent(in) :: groupHandle
+      integer, intent(in) :: ph
+    end subroutine thermal_restartWrite
+
     module subroutine mechanical_restartRead(groupHandle,ph)
       integer(HID_T), intent(in) :: groupHandle
       integer, intent(in) :: ph
     end subroutine mechanical_restartRead
 
+    module subroutine thermal_restartRead(groupHandle,ph)
+      integer(HID_T), intent(in) :: groupHandle
+      integer, intent(in) :: ph
+    end subroutine thermal_restartRead
 
     module function mechanical_S(ph,en) result(S)
       integer, intent(in) :: ph,en
@@ -640,6 +649,7 @@ subroutine phase_restartWrite(fileHandle)
     groupHandle(2) = HDF5_addGroup(groupHandle(1),material_name_phase(ph))
 
     call mechanical_restartWrite(groupHandle(2),ph)
+    call thermal_restartWrite(groupHandle(2),ph)
 
     call HDF5_closeGroup(groupHandle(2))
 
@@ -668,6 +678,7 @@ subroutine phase_restartRead(fileHandle)
     groupHandle(2) = HDF5_openGroup(groupHandle(1),material_name_phase(ph))
 
     call mechanical_restartRead(groupHandle(2),ph)
+    call thermal_restartRead(groupHandle(2),ph)
 
     call HDF5_closeGroup(groupHandle(2))
 
diff --git a/src/phase_thermal.f90 b/src/phase_thermal.f90
index a3e4dd628..fb98039fd 100644
--- a/src/phase_thermal.f90
+++ b/src/phase_thermal.f90
@@ -254,6 +254,38 @@ function integrateThermalState(Delta_t, ph,en) result(broken)
 end function integrateThermalState
 
 
+module subroutine thermal_restartWrite(groupHandle,ph)
+
+  integer(HID_T), intent(in) :: groupHandle
+  integer, intent(in) :: ph
+
+  integer(HID_T) :: new_group
+  integer :: so
+
+  do so = 1,thermal_Nsources(ph)
+    !new_group = HDF5_addGroup(groupHandle,
+    call HDF5_write(thermalState(ph)%p(so)%state,groupHandle,'omega2')
+  enddo
+
+end subroutine thermal_restartWrite
+
+
+module subroutine thermal_restartRead(groupHandle,ph)
+
+  integer(HID_T), intent(in) :: groupHandle
+  integer, intent(in) :: ph
+
+  integer(HID_T) :: new_group
+  integer :: so
+
+  do so = 1,thermal_Nsources(ph)
+    !new_group = HDF5_addGroup(groupHandle,
+    call HDF5_read(thermalState(ph)%p(so)%state0,groupHandle,'omega2')
+  enddo
+
+end subroutine thermal_restartRead
+
+
 module subroutine thermal_forward()
 
   integer :: ph, so