diff --git a/src/grid/grid_thermal_spectral.f90 b/src/grid/grid_thermal_spectral.f90
index 193d85c23..cc4b89afb 100644
--- a/src/grid/grid_thermal_spectral.f90
+++ b/src/grid/grid_thermal_spectral.f90
@@ -57,7 +57,7 @@ module grid_thermal_spectral
   integer                     :: totalIter = 0                                                      !< total iteration in current increment
   real(pREAL), dimension(3,3) :: K_ref
   real(pREAL)                 :: mu_ref, Delta_t_
-  logical :: broken
+  integer(kind(STATUS_OK))    :: status
 
   public :: &
     grid_thermal_spectral_init, &
@@ -208,7 +208,7 @@ function grid_thermal_spectral_solution(Delta_t) result(solution)
   call SNESGetConvergedReason(SNES_thermal,reason,err_PETSc)
   CHKERRQ(err_PETSc)
 
-  solution%converged = reason > 0 .and. .not. broken
+  solution%converged = reason > 0 .and. status == STATUS_OK
   solution%iterationsNeeded = merge(totalIter,num%itmax,solution%converged)
 
   call SNESGetDM(SNES_thermal,DM_thermal,err_PETSc)
@@ -323,11 +323,9 @@ subroutine formResidual(residual_subdomain,x_scal,r,dummy,err_PETSc)
 
   integer :: i, j, k, ce
   real(pREAL), dimension(3,cells(1),cells(2),cells3) :: vectorField
-  integer(kind(STATUS_OK)) :: status
 
 
   call homogenization_thermal_response(status,Delta_t_,1,product(cells(1:2))*cells3)
-  broken = status /= STATUS_OK
 
   associate(T => x_scal)
     vectorField = utilities_ScalarGradient(T)
diff --git a/src/mesh/FEM_utilities.f90 b/src/mesh/FEM_utilities.f90
index 64509ed87..ea4e4b87f 100644
--- a/src/mesh/FEM_utilities.f90
+++ b/src/mesh/FEM_utilities.f90
@@ -121,21 +121,19 @@ end subroutine FEM_utilities_init
 !--------------------------------------------------------------------------------------------------
 !> @brief calculates constitutive response
 !--------------------------------------------------------------------------------------------------
-subroutine utilities_constitutiveResponse(broken, Delta_t,P_av,forwardData)
+subroutine utilities_constitutiveResponse(status, Delta_t,P_av,forwardData)
 
-  logical,     intent(out)                :: broken
+  integer(kind(STATUS_OK)),  intent(out)  :: status
   real(pREAL), intent(in)                 :: Delta_t                                                !< loading time
   logical,     intent(in)                 :: forwardData                                            !< age results
   real(pREAL),intent(out), dimension(3,3) :: P_av                                                   !< average PK stress
 
   integer(MPI_INTEGER_KIND) :: err_MPI
-  integer(kind(STATUS_OK)) :: status
 
 
   print'(/,1x,a)', '... evaluating constitutive response ......................................'
 
   call homogenization_mechanical_response(status,Delta_t,1,mesh_maxNips*mesh_NcpElems)              ! calculate P field
-  broken = status /= STATUS_OK
   cutBack = .false.
 
   P_av = sum(homogenization_P,dim=3) * wgt
diff --git a/src/mesh/mesh_mech_FEM.f90 b/src/mesh/mesh_mech_FEM.f90
index 0baebe6c3..6eafe66e3 100644
--- a/src/mesh/mesh_mech_FEM.f90
+++ b/src/mesh/mesh_mech_FEM.f90
@@ -25,6 +25,7 @@ module mesh_mechanical_FEM
   use FEM_quadrature
   use homogenization
   use math
+  use constants
 
 #if (PETSC_VERSION_MAJOR==3 && PETSC_VERSION_MINOR>14) && !defined(PETSC_HAVE_MPI_F90MODULE_VISIBILITY)
   implicit none(type,external)
@@ -68,7 +69,8 @@ module mesh_mechanical_FEM
   character(len=pSTRLEN) :: incInfo
   real(pREAL), dimension(3,3) :: &
     P_av = 0.0_pREAL
-  logical :: ForwardData, broken
+  logical :: ForwardData
+  integer(kind(STATUS_OK)) :: status
   real(pREAL), parameter :: eps = 1.0e-18_pREAL
 
   external :: &                                                                                     ! ToDo: write interfaces
@@ -311,7 +313,7 @@ subroutine FEM_mechanical_init(mechBC,num_mesh)
     call DMPlexVecSetClosure(mechanical_mesh,section,solution_local,cell,px_scal,5,err_PETSc)
     CHKERRQ(err_PETSc)
   end do
-  call utilities_constitutiveResponse(broken,0.0_pREAL,devNull,.true.)
+  call utilities_constitutiveResponse(status,0.0_pREAL,devNull,.true.)
 
 end subroutine FEM_mechanical_init
 
@@ -458,8 +460,8 @@ subroutine FEM_mechanical_formResidual(dm_local,xx_local,f_local,dummy,err_PETSc
 
 !--------------------------------------------------------------------------------------------------
 ! evaluate constitutive response
-  call utilities_constitutiveResponse(broken,params%Delta_t,P_av,ForwardData)
-  call MPI_Allreduce(MPI_IN_PLACE,broken,1_MPI_INTEGER_KIND,MPI_LOGICAL,MPI_LOR,MPI_COMM_WORLD,err_MPI)
+  call utilities_constitutiveResponse(status,params%Delta_t,P_av,ForwardData)
+  call MPI_Allreduce(MPI_IN_PLACE,status,1_MPI_INTEGER_KIND,MPI_INTEGER,MPI_SUM,MPI_COMM_WORLD,err_MPI)
   call parallelization_chkerr(err_MPI)
   ForwardData = .false.
 
@@ -747,7 +749,7 @@ subroutine FEM_mechanical_converged(snes_local,PETScIter,xnorm,snorm,fnorm,reaso
   divTol = max(maxval(abs(P_av(1:dimPlex,1:dimPlex)))*num%eps_struct_rtol,num%eps_struct_atol)
   call SNESConvergedDefault(snes_local,PETScIter,xnorm,snorm,fnorm/divTol,reason,dummy,err_PETSc)
   CHKERRQ(err_PETSc)
-  if (broken) reason = SNES_DIVERGED_FUNCTION_DOMAIN
+  if (status /= STATUS_OK) reason = SNES_DIVERGED_FUNCTION_DOMAIN
   print'(/,1x,a,a,i0,a,f0.3)', trim(incInfo), &
                   ' @ Iteration ',PETScIter,' mechanical residual norm = ',fnorm/divTol
   print'(/,1x,a,/,2(3(2x,f12.4,1x)/),3(2x,f12.4,1x))', &