diff --git a/CMakeLists.txt b/CMakeLists.txt
index 9770996b1..f5d6546a9 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -124,6 +124,7 @@ endif ()
 
 # Predefined sets for OPTIMIZATION/OPENMP based on BUILD_TYPE
 if ("${CMAKE_BUILD_TYPE}" STREQUAL "DEBUG" OR "${CMAKE_BUILD_TYPE}" STREQUAL "SYNTAXONLY" )
+    set (DEBUG_FLAGS "${DEBUG_FLAGS} -DDEBUG")
     set (PARALLEL "OFF")
     set (OPTI "OFF")
 elseif ("${CMAKE_BUILD_TYPE}" STREQUAL "RELEASE")
diff --git a/src/DAMASK_spectral.f90 b/src/DAMASK_spectral.f90
index dc529b2e8..f32bfb7b3 100644
--- a/src/DAMASK_spectral.f90
+++ b/src/DAMASK_spectral.f90
@@ -447,8 +447,8 @@ program DAMASK_spectral
                              min(i*((maxRealOut)/materialpoint_sizeResults),size(materialpoint_results,3))],pLongInt)
      call MPI_file_write(resUnit, &
                          reshape(materialpoint_results(:,:,outputIndex(1):outputIndex(2)), &
-                                [(outputIndex(2)-outputIndex(1)+1)*materialpoint_sizeResults]), &
-                         (outputIndex(2)-outputIndex(1)+1)*materialpoint_sizeResults, &
+                                [(outputIndex(2)-outputIndex(1)+1)*int(materialpoint_sizeResults,pLongInt)]), &
+                         (outputIndex(2)-outputIndex(1)+1)*int(materialpoint_sizeResults,pLongInt), &
                          MPI_DOUBLE, MPI_STATUS_IGNORE, ierr)
      if (ierr /= 0_pInt) call IO_error(error_ID=894_pInt, ext_msg='MPI_file_write')
    enddo
@@ -645,8 +645,8 @@ program DAMASK_spectral
            outputIndex=int([(i-1_pInt)*((maxRealOut)/materialpoint_sizeResults)+1_pInt, &
                       min(i*((maxRealOut)/materialpoint_sizeResults),size(materialpoint_results,3))],pLongInt)
            call MPI_file_write(resUnit,reshape(materialpoint_results(:,:,outputIndex(1):outputIndex(2)),&
-                                         [(outputIndex(2)-outputIndex(1)+1)*materialpoint_sizeResults]), &
-                               (outputIndex(2)-outputIndex(1)+1)*materialpoint_sizeResults,&
+                                         [(outputIndex(2)-outputIndex(1)+1)*int(materialpoint_sizeResults,pLongInt)]), &
+                               (outputIndex(2)-outputIndex(1)+1)*int(materialpoint_sizeResults,pLongInt),&
                                MPI_DOUBLE, MPI_STATUS_IGNORE, ierr)
            if(ierr /=0_pInt) call IO_error(894_pInt, ext_msg='MPI_file_write')
          enddo
@@ -727,10 +727,10 @@ subroutine quit(stop_id)
  call utilities_destroy()
 
  call PETScFinalize(ierr)
- if(ierr /= 0) write(6,'(a)') ' Error in PETScFinalize'
+ if (ierr /= 0) write(6,'(a)') ' Error in PETScFinalize'
 #ifdef _OPENMP
  call MPI_finalize(error)
- if(error /= 0) write(6,'(a)') ' Error in MPI_finalize'
+ if (error /= 0) write(6,'(a)') ' Error in MPI_finalize'
 #endif
  ErrorInQuit = (ierr /= 0 .or. error /= 0_pInt)
  
@@ -744,7 +744,7 @@ subroutine quit(stop_id)
                                                       dateAndTime(7)
 
  if (stop_id == 0_pInt .and. .not. ErrorInQuit) stop 0                                              ! normal termination
- if (stop_id <  0_pInt .and. .not. ErrorInQuit) then                                                                        ! terminally ill, restart might help
+ if (stop_id <  0_pInt .and. .not. ErrorInQuit) then                                                ! terminally ill, restart might help
    write(0,'(a,i6)') 'restart information available at ', stop_id*(-1_pInt)
    stop 2
  endif
diff --git a/src/crystallite.f90 b/src/crystallite.f90
index c1af4424b..7f4f6fa44 100644
--- a/src/crystallite.f90
+++ b/src/crystallite.f90
@@ -554,7 +554,6 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
    FEsolving_execIP
  use mesh, only: &
    mesh_element, &
-   mesh_NcpElems, &
    mesh_maxNips, &
    mesh_ipNeighborhood, &
    FE_NipNeighbors, &
@@ -565,8 +564,7 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
    plasticState, &
    sourceState, &
    phase_Nsources, &
-   phaseAt, phasememberAt, &
-   homogenization_maxNgrains
+   phaseAt, phasememberAt
  use constitutive, only:  &
    constitutive_TandItsTangent, &
    constitutive_LpAndItsTangent, &
@@ -794,7 +792,7 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
                        if (.not. crystallite_localPlasticity(1,neighboring_i,neighboring_e) &
                            .and. .not. crystallite_converged(1,neighboring_i,neighboring_e)) then
                          crystallite_neighborEnforcedCutback(i,e) = .true.
-#ifndef _OPENMP
+#ifdef DEBUG
                          if (iand(debug_level(debug_crystallite),debug_levelExtensive) /= 0_pInt) &
                             write(6,'(a12,i5,1x,i2,a,i5,1x,i2)') '<< CRYST >> ', neighboring_e,neighboring_i, &
                                                                  ' enforced cutback at ',e,i
@@ -829,7 +827,7 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
                      if (.not. crystallite_localPlasticity(1,neighboring_i,neighboring_e) &
                          .and. .not. crystallite_converged(1,neighboring_i,neighboring_e)) then
                        crystallite_syncSubFrac(i,e) = .true.
-#ifndef _OPENMP
+#ifdef DEBUG
                        if (iand(debug_level(debug_crystallite),debug_levelExtensive) /= 0_pInt) &
                          write(6,'(a12,i5,1x,i2,a,i5,1x,i2)') '<< CRYST >> ',neighboring_e,neighboring_i, &
                                                              ' enforced time synchronization at ',e,i
@@ -937,7 +935,7 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
                  crystallite_todo(c,i,e) = .true.
                endif
                !$OMP FLUSH(crystallite_todo)
-#ifndef _OPENMP
+#ifdef DEBUG
                if (iand(debug_level(debug_crystallite),debug_levelBasic) /= 0_pInt &
                    .and. ((e == debug_e .and. i == debug_i .and. c == debug_g) &
                           .or. .not. iand(debug_level(debug_crystallite), debug_levelSelective) /= 0_pInt)) &
@@ -987,7 +985,7 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
                                                                                                      ! cant restore dotState here, since not yet calculated in first cutback after initialization
              crystallite_todo(c,i,e) = crystallite_subStep(c,i,e) > subStepMinCryst                  ! still on track or already done (beyond repair)
              !$OMP FLUSH(crystallite_todo)
-#ifndef _OPENMP
+#ifdef DEBUG
              if (iand(debug_level(debug_crystallite),debug_levelBasic) /= 0_pInt) then
                if (crystallite_todo(c,i,e)) then
                  write(6,'(a,f12.8,a,i8,1x,i2,1x,i3,/)') '<< CRYST >> cutback step in crystallite_stressAndItsTangent &
@@ -1393,7 +1391,7 @@ subroutine crystallite_integrateStateRK4()
          * crystallite_subdt(g,i,e) * timeStepFraction(n)
          enddo
 
-#ifndef _OPENMP
+#ifdef DEBUG
          if (n == 4 &
              .and. iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt &
              .and. ((e == debug_e .and. i == debug_i .and. g == debug_g) &
@@ -1784,7 +1782,7 @@ subroutine crystallite_integrateStateRKCK45()
 
 
    ! --- dot state and RK dot state---
-#ifndef _OPENMP
+#ifdef DEBUG
    if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt) &
      write(6,'(a,1x,i1)') '<< CRYST >> Runge--Kutta step',stage+1_pInt
 #endif
@@ -1933,7 +1931,7 @@ subroutine crystallite_integrateStateRKCK45()
                                        sourceState(p)%p(mySource)%aTolState(1:mySizeSourceDotState))
        enddo
 
-#ifndef _OPENMP
+#ifdef DEBUG
        if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt&
            .and. ((e == debug_e .and. i == debug_i .and. g == debug_g)&
                   .or. .not. iand(debug_level(debug_crystallite), debug_levelSelective) /= 0_pInt)) then
@@ -2317,7 +2315,7 @@ subroutine crystallite_integrateStateAdaptiveEuler()
          !$OMP FLUSH(relPlasticStateResiduum)
          !$OMP FLUSH(relSourceStateResiduum)
 
-#ifndef _OPENMP
+#ifdef DEBUG
 
          if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt &
              .and. ((e == debug_e .and. i == debug_i .and. g == debug_g)&
@@ -2513,7 +2511,7 @@ eIter = FEsolving_execElem(1:2)
          * crystallite_subdt(g,i,e)
          enddo
 
-#ifndef _OPENMP
+#ifdef DEBUG
          if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt &
              .and. ((e == debug_e .and. i == debug_i .and. g == debug_g) &
                      .or. .not. iand(debug_level(debug_crystallite), debug_levelSelective) /= 0_pInt)) then
@@ -2962,7 +2960,7 @@ subroutine crystallite_integrateStateFPI()
            * (1.0_pReal - sourceStateDamper)
          enddo
 
-#ifndef _OPENMP
+#ifdef DEBUG
          if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt &
              .and. ((e == debug_e .and. i == debug_i .and. g == debug_g) &
                     .or. .not. iand(debug_level(debug_crystallite), debug_levelSelective) /= 0_pInt)) then
@@ -3134,7 +3132,7 @@ logical function crystallite_stateJump(ipc,ip,el)
      sourceState(p)%p(mySource)%deltaState(1:mySizeSourceDeltaState,c)
  enddo
 
-#ifndef _OPENMP
+#ifdef DEBUG
  if (any(dNeq0(plasticState(p)%deltaState(1:mySizePlasticDeltaState,c))) &
      .and. iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt &
      .and. ((el == debug_e .and. ip == debug_i .and. ipc == debug_g) &
@@ -3309,7 +3307,7 @@ logical function crystallite_integrateStress(&
 
  !* be pessimistic
  crystallite_integrateStress = .false.
-#ifndef _OPENMP
+#ifdef DEBUG
  if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt &
      .and. ((el == debug_e .and. ip == debug_i .and. ipc == debug_g) &
             .or. .not. iand(debug_level(debug_crystallite), debug_levelSelective) /= 0_pInt)) &
@@ -3342,9 +3340,9 @@ logical function crystallite_integrateStress(&
 
  invFp_current = math_inv33(Fp_current)
  failedInversionFp: if (all(dEq0(invFp_current))) then
-#ifndef _OPENMP
+#ifdef DEBUG
    if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0_pInt) then
-     write(6,'(a,i8,1x,a,i8,a,1x,i2,1x,i3)') '<< CRYST >> integrateStress failed on inversion of Fp_current at el (elFE) ip g ',&
+     write(6,'(a,i8,1x,a,i8,a,1x,i2,1x,i3)') '<< CRYST >> integrateStress failed on inversion of Fp_current at el (elFE) ip ipc ',&
        el,'(',mesh_element(1,el),')',ip,ipc
      if (iand(debug_level(debug_crystallite), debug_levelExtensive) > 0_pInt) &
        write(6,'(/,a,/,3(12x,3(f12.7,1x)/))') '<< CRYST >> Fp_current',math_transpose33(Fp_current(1:3,1:3))
@@ -3358,7 +3356,7 @@ logical function crystallite_integrateStress(&
 
  invFi_current = math_inv33(Fi_current)
  failedInversionFi: if (all(dEq0(invFi_current))) then
-#ifndef _OPENMP
+#ifdef DEBUG
    if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0_pInt) then
      write(6,'(a,i8,1x,a,i8,a,1x,i2,1x,i3)') '<< CRYST >> integrateStress failed on inversion of Fi_current at el (elFE) ip ipc ',&
        el,'(',mesh_element(1,el),')',ip,ipc
@@ -3379,10 +3377,10 @@ logical function crystallite_integrateStress(&
  LiLoop: do
    NiterationStressLi = NiterationStressLi + 1_pInt
    IloopsExeced: if (NiterationStressLi > nStress) then
-#ifndef _OPENMP
+#ifdef DEBUG
      if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0_pInt) &
        write(6,'(a,i3,a,i8,1x,a,i8,a,1x,i2,1x,i3,/)') '<< CRYST >> integrateStress reached inelastic loop limit',nStress, &
-       ' at el (elFE) ip ipc ', el,mesh_element(1,el),ip,ipc
+       ' at el (elFE) ip ipc ', el,'(',mesh_element(1,el),')',ip,ipc
 #endif
      return
    endif IloopsExeced
@@ -3400,7 +3398,7 @@ logical function crystallite_integrateStress(&
    LpLoop: do                                    ! inner stress integration loop for consistency with Fi
      NiterationStressLp = NiterationStressLp + 1_pInt
      loopsExeced: if (NiterationStressLp > nStress) then
-#ifndef _OPENMP
+#ifdef DEBUG
      if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0_pInt) &
          write(6,'(a,i3,a,i8,1x,a,i8,a,1x,i2,1x,i3,/)') '<< CRYST >> integrateStress reached loop limit',nStress, &
          ' at el (elFE) ip ipc ', el,mesh_element(1,el),ip,ipc
@@ -3433,7 +3431,7 @@ logical function crystallite_integrateStress(&
        !$OMP END CRITICAL (debugTimingLpTangent)
      endif
 
-#ifndef _OPENMP
+#ifdef DEBUG
      if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt &
          .and. ((el == debug_e .and. ip == debug_i .and. ipc == debug_g) &
                 .or. .not. iand(debug_level(debug_crystallite), debug_levelSelective) /= 0_pInt)) then
@@ -3450,11 +3448,11 @@ logical function crystallite_integrateStress(&
                   aTol_crystalliteStress)                                                            ! minimum lower cutoff
      residuumLp = Lpguess - Lp_constitutive
 
-     if (any(IEEE_is_NaN(residuumLp))) then                                                         ! NaN in residuum...
-#ifndef _OPENMP
+     if (any(IEEE_is_NaN(residuumLp))) then                                                          ! NaN in residuum...
+#ifdef DEBUG
        if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0_pInt) &
          write(6,'(a,i8,1x,a,i8,a,1x,i2,1x,i3,a,i3,a)') '<< CRYST >> integrateStress encountered NaN at el (elFE) ip ipc ', &
-           el,mesh_element(1,el),ip,ipc, &
+           el,'(',mesh_element(1,el),')',ip,ipc, &
                                                         ' ; iteration ', NiterationStressLp,&
                                                         ' >> returning..!'
 #endif
@@ -3486,10 +3484,10 @@ logical function crystallite_integrateStress(&
        work = math_plain33to9(residuumLp)
        call dgesv(9,1,dRLp_dLp2,9,ipiv,work,9,ierr)                                                   ! solve dRLp/dLp * delta Lp = -res for delta Lp
        if (ierr /= 0_pInt) then
-#ifndef _OPENMP
+#ifdef DEBUG
          if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0_pInt) then
-           write(6,'(a,i8,1x,a,i8,a,1x,i2,1x,i3,a,i3)') '<< CRYST >> integrateStress failed on dR/dLp inversion at el ip ipc ', &
-             el,mesh_element(1,el),ip,ipc
+           write(6,'(a,i8,1x,a,i8,a,1x,i2,1x,i3)') '<< CRYST >> integrateStress failed on dR/dLp inversion at el (elFE) ip ipc ', &
+             el,'(',mesh_element(1,el),')',ip,ipc
            if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt &
                .and. ((el == debug_e .and. ip == debug_i .and. ipc == debug_g)&
                       .or. .not. iand(debug_level(debug_crystallite), debug_levelSelective) /= 0_pInt)) then
@@ -3527,7 +3525,7 @@ logical function crystallite_integrateStress(&
    call constitutive_LiAndItsTangent(Li_constitutive, dLi_dT3333, dLi_dFi3333, &
                                      Tstar_v, Fi_new, ipc, ip, el)
 
-#ifndef _OPENMP
+#ifdef DEBUG
      if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt &
          .and. ((el == debug_e .and. ip == debug_i .and. ipc == debug_g) &
                 .or. .not. iand(debug_level(debug_crystallite), debug_levelSelective) /= 0_pInt)) then
@@ -3575,10 +3573,10 @@ logical function crystallite_integrateStress(&
      work = math_plain33to9(residuumLi)
      call dgesv(9,1,dRLi_dLi,9,ipiv,work,9,ierr)                                                    ! solve dRLi/dLp * delta Li = -res for delta Li
      if (ierr /= 0_pInt) then
-#ifndef _OPENMP
+#ifdef DEBUG
        if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0_pInt) then
-         write(6,'(a,i8,1x,a,i8,a,1x,i2,1x,i3,a,i3)') '<< CRYST >> integrateStress failed on dR/dLi inversion at el ip ipc ', &
-               el,mesh_element(1,el),ip,ipc
+         write(6,'(a,i8,1x,a,i8,a,1x,i2,1x,i3)') '<< CRYST >> integrateStress failed on dR/dLi inversion at el (elFE) ip ipc ', &
+               el,'(',mesh_element(1,el),')',ip,ipc
          if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt &
              .and. ((el == debug_e .and. ip == debug_i .and. ipc == debug_g)&
                     .or. .not. iand(debug_level(debug_crystallite), debug_levelSelective) /= 0_pInt)) then
@@ -3615,10 +3613,10 @@ logical function crystallite_integrateStress(&
  invFp_new = invFp_new / math_det33(invFp_new)**(1.0_pReal/3.0_pReal)                               ! regularize by det
  Fp_new = math_inv33(invFp_new)
  failedInversionInvFp: if (all(dEq0(Fp_new))) then
-#ifndef _OPENMP
+#ifdef DEBUG
    if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0_pInt) then
-     write(6,'(a,i8,1x,a,i8,a,1x,i2,1x,i3,a,i3)') '<< CRYST >> integrateStress failed on invFp_new inversion at el ip ipc ',&
-       el,mesh_element(1,el),ip,ipc, ' ; iteration ', NiterationStressLp
+     write(6,'(a,i8,1x,a,i8,a,1x,i2,1x,i3,a,i3)') '<< CRYST >> integrateStress failed on invFp_new inversion at el (elFE) ip ipc ',&
+       el,'(',mesh_element(1,el),')',ip,ipc, ' ; iteration ', NiterationStressLp
      if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt &
          .and. ((el == debug_e .and. ip == debug_i .and. ipc == debug_g) &
                 .or. .not. iand(debug_level(debug_crystallite), debug_levelSelective) /= 0_pInt)) &
@@ -3649,7 +3647,7 @@ logical function crystallite_integrateStress(&
  !* set return flag to true
 
  crystallite_integrateStress = .true.
-#ifndef _OPENMP
+#ifdef DEBUG
  if (iand(debug_level(debug_crystallite),debug_levelExtensive) /= 0_pInt &
      .and. ((el == debug_e .and. ip == debug_i .and. ipc == debug_g) &
              .or. .not. iand(debug_level(debug_crystallite), debug_levelSelective) /= 0_pInt)) then
diff --git a/src/homogenization.f90 b/src/homogenization.f90
index 504f68e8c..5a30a72c8 100644
--- a/src/homogenization.f90
+++ b/src/homogenization.f90
@@ -542,6 +542,7 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt)
    debug_level, &
    debug_homogenization, &
    debug_levelBasic, &
+   debug_levelExtensive, &
    debug_levelSelective, &
    debug_e, &
    debug_i, &
@@ -638,8 +639,8 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt)
      IpLooping1: do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e)
 
        converged: if ( materialpoint_converged(i,e) ) then
-#ifndef _OPENMP
-         if (iand(debug_level(debug_homogenization), debug_levelBasic) /= 0_pInt &
+#ifdef DEBUG
+         if (iand(debug_level(debug_homogenization), debug_levelExtensive) /= 0_pInt &
             .and. ((e == debug_e .and. i == debug_i) &
                    .or. .not. iand(debug_level(debug_homogenization),debug_levelSelective) /= 0_pInt)) then
            write(6,'(a,1x,f12.8,1x,a,1x,f12.8,1x,a,i8,1x,i2/)') '<< HOMOG >> winding forward from', &
@@ -741,8 +742,8 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt)
            materialpoint_subStep(i,e) = subStepSizeHomog * materialpoint_subStep(i,e)               ! crystallite had severe trouble, so do a significant cutback
            !$OMP FLUSH(materialpoint_subStep)
 
-#ifndef _OPENMP
-           if (iand(debug_level(debug_homogenization), debug_levelBasic) /= 0_pInt &
+#ifdef DEBUG
+           if (iand(debug_level(debug_homogenization), debug_levelExtensive) /= 0_pInt &
               .and. ((e == debug_e .and. i == debug_i) &
                     .or. .not. iand(debug_level(debug_homogenization), debug_levelSelective) /= 0_pInt)) then
              write(6,'(a,1x,f12.8,a,i8,1x,i2/)') &
diff --git a/src/mesh.f90 b/src/mesh.f90
index 87160f2c7..666fe1e33 100644
--- a/src/mesh.f90
+++ b/src/mesh.f90
@@ -519,6 +519,8 @@ subroutine mesh_init(ip,el)
  integer(pInt) :: j
  logical :: myDebug
  
+ external :: MPI_comm_size
+
  write(6,'(/,a)')   ' <<<+-  mesh init  -+>>>'
  write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
 #include "compilation_info.f90"
diff --git a/src/plastic_nonlocal.f90 b/src/plastic_nonlocal.f90
index 0a8c4c3f9..55871737d 100644
--- a/src/plastic_nonlocal.f90
+++ b/src/plastic_nonlocal.f90
@@ -1823,14 +1823,14 @@ plasticState(ph)%state(iRhoF(1:ns,instance),of) = rhoForest
 plasticState(ph)%state(iTauF(1:ns,instance),of) = tauThreshold
 plasticState(ph)%state(iTauB(1:ns,instance),of) = tauBack
 
-#ifndef _OPENMP
+#ifdef DEBUG
   if (iand(debug_level(debug_constitutive),debug_levelExtensive) /= 0_pInt &
       .and. ((debug_e == el .and. debug_i == ip)&
              .or. .not. iand(debug_level(debug_constitutive),debug_levelSelective) /= 0_pInt)) then
     write(6,'(/,a,i8,1x,i2,1x,i1,/)') '<< CONST >> nonlocal_microstructure at el ip ',el,ip
     write(6,'(a,/,12x,12(e10.3,1x))') '<< CONST >> rhoForest', rhoForest
-    write(6,'(a,/,12x,12(f10.5,1x))') '<< CONST >> tauThreshold / MPa', tauThreshold/1e6
-    write(6,'(a,/,12x,12(f10.5,1x),/)') '<< CONST >> tauBack / MPa', tauBack/1e6
+    write(6,'(a,/,12x,12(f10.5,1x))') '<< CONST >> tauThreshold / MPa', tauThreshold*1e-6
+    write(6,'(a,/,12x,12(f10.5,1x),/)') '<< CONST >> tauBack / MPa', tauBack*1e-6
   endif
 #endif
 
@@ -1978,15 +1978,15 @@ if (Temperature > 0.0_pReal) then
 endif
     
 
-#ifndef _OPENMP
+#ifdef DEBUG
   if (iand(debug_level(debug_constitutive),debug_levelExtensive) /= 0_pInt &
       .and. ((debug_e == el .and. debug_i == ip)&
              .or. .not. iand(debug_level(debug_constitutive),debug_levelSelective) /= 0_pInt)) then
     write(6,'(/,a,i8,1x,i2,1x,i1,/)') '<< CONST >> nonlocal_kinetics at el ip',el,ip
-    write(6,'(a,/,12x,12(f12.5,1x))') '<< CONST >> tauThreshold / MPa', tauThreshold / 1e6_pReal
-    write(6,'(a,/,12x,12(f12.5,1x))') '<< CONST >> tau / MPa', tau / 1e6_pReal
-    write(6,'(a,/,12x,12(f12.5,1x))') '<< CONST >> tauNS / MPa', tauNS / 1e6_pReal
-    write(6,'(a,/,12x,12(f12.5,1x))') '<< CONST >> v / 1e-3m/s', v * 1e3
+    write(6,'(a,/,12x,12(f12.5,1x))') '<< CONST >> tauThreshold / MPa', tauThreshold * 1e-6_pReal
+    write(6,'(a,/,12x,12(f12.5,1x))') '<< CONST >> tau / MPa', tau * 1e-6_pReal
+    write(6,'(a,/,12x,12(f12.5,1x))') '<< CONST >> tauNS / MPa', tauNS * 1e-6_pReal
+    write(6,'(a,/,12x,12(f12.5,1x))') '<< CONST >> v / mm/s', v * 1e3
     write(6,'(a,/,12x,12(e12.5,1x))') '<< CONST >> dv_dtau', dv_dtau
     write(6,'(a,/,12x,12(e12.5,1x))') '<< CONST >> dv_dtauNS', dv_dtauNS
   endif
@@ -2176,12 +2176,12 @@ enddo
 dLp_dTstar99 = math_Plain3333to99(dLp_dTstar3333)
 
 
-#ifndef _OPENMP
+#ifdef DEBUG
   if (iand(debug_level(debug_constitutive),debug_levelExtensive) /= 0_pInt &
       .and. ((debug_e == el .and. debug_i == ip)&
              .or. .not. iand(debug_level(debug_constitutive),debug_levelSelective) /= 0_pInt )) then
     write(6,'(/,a,i8,1x,i2,1x,i1,/)') '<< CONST >> nonlocal_LpandItsTangent at el ip',el,ip
-    write(6,'(a,/,12x,12(f12.5,1x))') '<< CONST >> gdot total / 1e-3',gdotTotal*1e3_pReal
+    write(6,'(a,/,12x,12(f12.5,1x))') '<< CONST >> gdot total',gdotTotal
     write(6,'(a,/,3(12x,3(f12.7,1x),/))') '<< CONST >> Lp',transpose(Lp)
   endif
 #endif
@@ -2248,7 +2248,7 @@ real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(1,ip,e
                                             dUpperOld, &                  ! old maximum stable dipole distance for edges and screws
                                             deltaDUpper                   ! change in maximum stable dipole distance for edges and screws
 
-#ifndef _OPENMP
+#ifdef DEBUG
   if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0_pInt &
       .and. ((debug_e == el .and. debug_i == ip)&
              .or. .not. iand(debug_level(debug_constitutive),debug_levelSelective) /= 0_pInt)) &
@@ -2361,7 +2361,7 @@ forall (s = 1:ns, c = 1_pInt:2_pInt) &
   plasticState(ph)%deltaState(iRhoD(s,c,instance),of) = deltaRho(s,c+8_pInt)
 
 
-#ifndef _OPENMP
+#ifdef DEBUG
   if (iand(debug_level(debug_constitutive),debug_levelExtensive) /= 0_pInt &
       .and. ((debug_e == el .and. debug_i == ip)&
              .or. .not. iand(debug_level(debug_constitutive),debug_levelSelective) /= 0_pInt )) then
@@ -2522,11 +2522,11 @@ logical                                     considerEnteringFlux, &
 
 
 
-#ifndef _OPENMP
+#ifdef DEBUG
   if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0_pInt &
       .and. ((debug_e == el .and. debug_i == ip)&
              .or. .not. iand(debug_level(debug_constitutive),debug_levelSelective) /= 0_pInt)) &
-    write(6,'(/,a,i8,1x,i2,1x,i1,/)') '<< CONST >> nonlocal_dotState at el ip ',el,ip
+    write(6,'(/,a,i8,1x,i2,/)') '<< CONST >> nonlocal_dotState at el ip ',el,ip
 #endif
 
 ph = material_phase(1_pInt,ip,el)
@@ -2589,7 +2589,7 @@ endif
 forall (t = 1_pInt:4_pInt) &
   gdot(1_pInt:ns,t) = rhoSgl(1_pInt:ns,t) * burgers(1:ns,instance) * v(1:ns,t)
 
-#ifndef _OPENMP
+#ifdef DEBUG
   if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0_pInt &
       .and. ((debug_e == el .and. debug_i == ip)&
              .or. .not. iand(debug_level(debug_constitutive),debug_levelSelective) /= 0_pInt )) then
@@ -2663,7 +2663,7 @@ else
           / burgers(s,instance) * sqrt(rhoForest(s)) / lambda0(s,instance)
       endif
     enddo
-#ifndef _OPENMP
+#ifdef DEBUG
     if (iand(debug_level(debug_constitutive),debug_levelExtensive) /= 0_pInt &
         .and. ((debug_e == el .and. debug_i == ip)&
              .or. .not. iand(debug_level(debug_constitutive),debug_levelSelective) /= 0_pInt )) &
@@ -2690,7 +2690,7 @@ if (.not. phase_localPlasticity(material_phase(1_pInt,ip,el))) then
   if (any( abs(gdot) > 0.0_pReal &                                                                  ! any active slip system ...
           .and. CFLfactor(instance) * abs(v) * timestep &
               > mesh_ipVolume(ip,el) / maxval(mesh_ipArea(:,ip,el)))) then                          ! ...with velocity above critical value (we use the reference volume and area for simplicity here)
-#ifndef _OPENMP
+#ifdef DEBUG
     if (iand(debug_level(debug_constitutive),debug_levelExtensive) /= 0_pInt) then 
       write(6,'(a,i5,a,i2)') '<< CONST >> CFL condition not fullfilled at el ',el,' ip ',ip
       write(6,'(a,e10.3,a,e10.3)') '<< CONST >> velocity is at  ', &
@@ -2952,7 +2952,7 @@ if (numerics_integrationMode == 1_pInt) then
 endif
 
 
-#ifndef _OPENMP
+#ifdef DEBUG
   if (iand(debug_level(debug_constitutive),debug_levelExtensive) /= 0_pInt &
       .and. ((debug_e == el .and. debug_i == ip .and. debug_g == 1_pInt)&
              .or. .not. iand(debug_level(debug_constitutive),debug_levelSelective) /= 0_pInt )) then
@@ -2978,7 +2978,7 @@ endif
 
 if (    any(rhoSglOriginal(1:ns,1:4) + rhoDot(1:ns,1:4) * timestep < -aTolRho(instance)) &
    .or. any(rhoDipOriginal(1:ns,1:2) + rhoDot(1:ns,9:10) * timestep < -aTolRho(instance))) then
-#ifndef _OPENMP
+#ifdef DEBUG
   if (iand(debug_level(debug_constitutive),debug_levelExtensive) /= 0_pInt) then 
     write(6,'(a,i5,a,i2)') '<< CONST >> evolution rate leads to negative density at el ',el,' ip ',ip
     write(6,'(a)') '<< CONST >> enforcing cutback !!!'
diff --git a/src/spectral_mech_Basic.f90 b/src/spectral_mech_Basic.f90
index cfb727129..55403ee7c 100644
--- a/src/spectral_mech_Basic.f90
+++ b/src/spectral_mech_Basic.f90
@@ -310,7 +310,6 @@ subroutine BasicPETSC_formResidual(in,x_scal,f_scal,dummy,ierr)
    debug_spectral, &
    debug_spectralRotation
  use spectral_utilities, only: &
-   wgt, &
    tensorField_real, &
    utilities_FFTtensorForward, &
    utilities_fourierGammaConvolution, &
diff --git a/src/spectral_utilities.f90 b/src/spectral_utilities.f90
index 1ad251747..1bbf2e608 100644
--- a/src/spectral_utilities.f90
+++ b/src/spectral_utilities.f90
@@ -970,6 +970,9 @@ subroutine utilities_constitutiveResponse(F_lastInc,F,timeinc, &
  real(pReal), dimension(3,3,3,3) :: max_dPdF, min_dPdF
  real(pReal)   :: max_dPdF_norm, min_dPdF_norm, defgradDetMin, defgradDetMax, defgradDet
 
+ external :: &
+  MPI_Allreduce
+
  write(6,'(/,a)') ' ... evaluating constitutive response ......................................'
  flush(6)
  age = .False.