diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
index 709655863..df8991900 100644
--- a/.gitlab-ci.yml
+++ b/.gitlab-ci.yml
@@ -1,3 +1,4 @@
+---
 stages:
   - prepareAll
   - preprocessing
@@ -23,26 +24,30 @@ stages:
 
 ###################################################################################################
 before_script:
-  - if [ $(awk "/$CI_PIPELINE_ID/{print NR}" $TESTROOT/GitLabCI.queue)x == 'x' ]; then echo $CI_PIPELINE_ID >> $TESTROOT/GitLabCI.queue; fi
-  - while [ $(awk "/$CI_PIPELINE_ID/{print NR}" $TESTROOT/GitLabCI.queue) -ne 1 ];do sleep 5m; done
-  - source $DAMASKROOT/DAMASK_env.sh
+  - if [ $(awk "/$CI_PIPELINE_ID/{print NR}" $TESTROOT/GitLabCI.queue)x == 'x' ];
+      then echo $CI_PIPELINE_ID >> $TESTROOT/GitLabCI.queue;
+    fi
+  - while [ $(awk "/$CI_PIPELINE_ID/{print NR}" $TESTROOT/GitLabCI.queue) != 1 ];
+      do sleep 5m;
+    done
+  - source $DAMASKROOT/env/DAMASK.sh
   - cd $DAMASKROOT/PRIVATE/testing
 
 ###################################################################################################
 variables:
-  #================================================================================================
+  # ===============================================================================================
   # GitLab Settings
-  #================================================================================================
+  # ===============================================================================================
   GIT_SUBMODULE_STRATEGY: none
 
-  #================================================================================================
+  # ===============================================================================================
   # Shortcut names
-  #================================================================================================
-  DAMASKROOT:                      "$TESTROOT/GitLabCI_Pipeline_$CI_PIPELINE_ID/DAMASK"
+  # ===============================================================================================
+  DAMASKROOT: "$TESTROOT/GitLabCI_Pipeline_$CI_PIPELINE_ID/DAMASK"
 
-  #================================================================================================
+  # ===============================================================================================
   # Names of module files to load
-  #================================================================================================
+  # ===============================================================================================
   # ++++++++++++ Compiler ++++++++++++++++++++++++++++++++++++++++++++++
   IntelCompiler16_0:               "Compiler/Intel/16.0 Libraries/IMKL/2016"
   IntelCompiler17_0:               "Compiler/Intel/17.0 Libraries/IMKL/2017"
@@ -81,18 +86,21 @@ variables:
 
 
 ###################################################################################################
-checkout: 
+checkout:
   stage: prepareAll
-  before_script: 
+  before_script:
     - echo $CI_PIPELINE_ID >> $TESTROOT/GitLabCI.queue
-    - while [ $(awk "/$CI_PIPELINE_ID/{print NR}" $TESTROOT/GitLabCI.queue) -ne 1 ];do sleep 5m; done
+    - while [ $(awk "/$CI_PIPELINE_ID/{print NR}" $TESTROOT/GitLabCI.queue) != 1 ];
+        do sleep 5m;
+      done
   script:
     - mkdir -p $DAMASKROOT
     - cd $DAMASKROOT
+    - if [ -d DAMASK ]; then rm -rf DAMASK; fi                                                      # there might be some leftovers from a failed clone
     - git clone -q git@magit1.mpie.de:damask/DAMASK.git .
     - git checkout $CI_COMMIT_SHA
     - git submodule update --init
-    - source DAMASK_env.sh
+    - source env/DAMASK.sh
     - make processing
   except:
     - master
@@ -186,7 +194,7 @@ Post_ParaviewRelated:
 
 ###################################################################################################
 Compile_Intel:
-  stage: compileSpectralIntel 
+  stage: compileSpectralIntel
   script:
     - module load $IntelCompiler $MPICH_Intel $PETSc_MPICH_Intel
     - SpectralAll_compile/test.py
@@ -214,7 +222,7 @@ Compile_Intel_Prepare:
   except:
     - master
     - release
- 
+
 ###################################################################################################
 Spectral_PackedGeometry:
   stage: spectral
@@ -442,9 +450,9 @@ SpectralRuntime:
   except:
     - master
     - release
- 
+
 ###################################################################################################
-AbaqusExp: 
+AbaqusExp:
   stage: createDocumentation
   script:
     - module load $IntelCompiler $MPICH_Intel $PETSc_MPICH_Intel $Doxygen
@@ -453,7 +461,7 @@ AbaqusExp:
     - master
     - release
 
-AbaqusStd: 
+AbaqusStd:
   stage: createDocumentation
   script:
     - module load $IntelCompiler $MPICH_Intel $PETSc_MPICH_Intel $Doxygen
@@ -462,7 +470,7 @@ AbaqusStd:
     - master
     - release
 
-Marc: 
+Marc:
   stage: createDocumentation
   script:
     - module load $IntelCompiler $MPICH_Intel $PETSc_MPICH_Intel $Doxygen
@@ -471,7 +479,7 @@ Marc:
     - master
     - release
 
-Spectral: 
+Spectral:
   stage: createDocumentation
   script:
     - module load $IntelCompiler $MPICH_Intel $PETSc_MPICH_Intel $Doxygen
@@ -481,7 +489,7 @@ Spectral:
     - release
 
 ##################################################################################################
-backupData: 
+backupData:
   stage: saveDocumentation
   script:
     - cd $TESTROOT/performance                                                                      # location of new runtime results
@@ -496,7 +504,7 @@ backupData:
     - development
 
 ##################################################################################################
-mergeIntoMaster: 
+mergeIntoMaster:
   stage: updateMaster
   script:
     - cd $DAMASKROOT
@@ -527,9 +535,9 @@ removeData:
     - release
 
 ###################################################################################################
-removeLock: 
+removeLock:
   stage: releaseLock
-  before_script: 
+  before_script:
     - echo 'Do nothing'
   when: always
   script: sed -i "/$CI_PIPELINE_ID/d" $TESTROOT/GitLabCI.queue
diff --git a/VERSION b/VERSION
index 5a76f3ec4..12e0d6df0 100644
--- a/VERSION
+++ b/VERSION
@@ -1 +1 @@
-v2.0.1-728-g66efd90
+v2.0.1-805-gdc3eda3
diff --git a/env/DAMASK.csh b/env/DAMASK.csh
index 81d4de421..819b263ef 100644
--- a/env/DAMASK.csh
+++ b/env/DAMASK.csh
@@ -3,16 +3,21 @@
 
 set CALLED=($_)
 set DIRNAME=`dirname $CALLED[2]`
+
+# transition compatibility (renamed $DAMASK_ROOT/DAMASK_env.csh to $DAMASK_ROOT/env/DAMASK.csh)
+set FILENAME=`basename $CALLED[2]`
+if ($FILENAME == "DAMASK.csh") then
+  set DIRNAME=$DIRNAME"/../"
+endif
+
 set DAMASK_ROOT=`python -c "import os,sys; print(os.path.realpath(os.path.expanduser(sys.argv[1])))" $DIRNAME`
 
+
 source $DAMASK_ROOT/CONFIG
 
-# if DAMASK_BIN is present and not in $PATH, add it
+# if DAMASK_BIN is present
 if ( $?DAMASK_BIN) then
-  set MATCH=`echo :${PATH}: | grep ${DAMASK_BIN}:`
-  if ( "x$MATCH" == "x" ) then
-    set PATH=${DAMASK_BIN}:${PATH}
-  endif
+  set path = ($DAMASK_BIN $path)
 endif
 
 set SOLVER=`which DAMASK_spectral`                                                          
@@ -62,4 +67,8 @@ if ( $?prompt ) then
 endif
 
 setenv DAMASK_NUM_THREADS $DAMASK_NUM_THREADS
-setenv PYTHONPATH $DAMASK_ROOT/lib:$PYTHONPATH
+if ( ! $?PYTHONPATH ) then
+  setenv PYTHONPATH $DAMASK_ROOT/lib
+else
+  setenv PYTHONPATH $DAMASK_ROOT/lib:$PYTHONPATH
+endif
diff --git a/env/DAMASK.sh b/env/DAMASK.sh
index 264ae9c94..3e9889fec 100644
--- a/env/DAMASK.sh
+++ b/env/DAMASK.sh
@@ -1,5 +1,5 @@
 # sets up an environment for DAMASK on bash
-# usage:  source DAMASK_env.sh
+# usage:  source DAMASK.sh
 
 
 if [ "$OSTYPE" == "linux-gnu" ] || [ "$OSTYPE" == 'linux' ]; then
@@ -10,8 +10,13 @@ else
   DAMASK_ROOT=${STAT##* }
 fi
 
+# transition compatibility (renamed $DAMASK_ROOT/DAMASK_env.sh to $DAMASK_ROOT/env/DAMASK.sh)
+if [ ${BASH_SOURCE##*/} == "DAMASK.sh" ]; then
+  DAMASK_ROOT=$DAMASK_ROOT'/..'
+fi
+
 # shorthand command to change to DAMASK_ROOT directory
-eval "function damask() { cd $DAMASK_ROOT; }"
+eval "function DAMASK_root() { cd $DAMASK_ROOT; }"
 
 # defining set() allows to source the same file for tcsh and bash, with and without space around =
 set() {
@@ -20,22 +25,16 @@ set() {
 source $DAMASK_ROOT/CONFIG
 unset -f set
 
-# add DAMASK_BIN if present but not yet in $PATH
-if [[ "x$DAMASK_BIN" != "x" && ! $(echo ":$PATH:" | grep $DAMASK_BIN:) ]]; then
-  export PATH=$DAMASK_BIN:$PATH
-fi
+# add DAMASK_BIN if present
+[ "x$DAMASK_BIN" != "x" ] && PATH=$DAMASK_BIN:$PATH
 
 SOLVER=$(which DAMASK_spectral || true 2>/dev/null)
-if [ "x$SOLVER" == "x" ]; then
-  SOLVER='Not found!'
-fi
+[ "x$SOLVER" == "x" ] && SOLVER='Not found!'
+
 PROCESSING=$(which postResults || true 2>/dev/null)
-if [ "x$PROCESSING" == "x" ]; then
-  PROCESSING='Not found!'
-fi
-if [ "x$DAMASK_NUM_THREADS" == "x" ]; then
-  DAMASK_NUM_THREADS=1
-fi
+[ "x$PROCESSING" == "x" ] && PROCESSING='Not found!'
+
+[ "x$DAMASK_NUM_THREADS" == "x" ] && DAMASK_NUM_THREADS=1
 
 # according to http://software.intel.com/en-us/forums/topic/501500
 # this seems to make sense for the stack size
diff --git a/env/DAMASK.zsh b/env/DAMASK.zsh
index 9098352bf..c662e6aee 100644
--- a/env/DAMASK.zsh
+++ b/env/DAMASK.zsh
@@ -1,10 +1,15 @@
 # sets up an environment for DAMASK on zsh
-# usage:  source DAMASK_env.zsh
+# usage:  source DAMASK.zsh
 
-DAMASK_ROOT=${0:a:h}
+# transition compatibility (renamed $DAMASK_ROOT/DAMASK_env.zsh to $DAMASK_ROOT/env/DAMASK.zsh)
+if [ ${0:t:r} = 'DAMASK' ]; then
+  DAMASK_ROOT=${0:a:h}'/..'
+else
+  DAMASK_ROOT=${0:a:h}
+fi
 
 # shorthand command to change to DAMASK_ROOT directory
-eval "function damask() { cd $DAMASK_ROOT; }"
+eval "function DAMASK_root() { cd $DAMASK_ROOT; }"
 
 # defining set() allows to source the same file for tcsh and zsh, with and without space around =
 set() {
@@ -13,17 +18,12 @@ set() {
 source $DAMASK_ROOT/CONFIG
 unset -f set
 
-# add DAMASK_BIN if present but not yet in $PATH
-MATCH=`echo ":$PATH:" | grep $DAMASK_BIN:`
-if [[ ( "x$DAMASK_BIN" != "x" ) && ( "x$MATCH" = "x" ) ]]; then
-  export PATH=$DAMASK_BIN:$PATH
-fi
+# add DAMASK_BIN if present
+[ "x$DAMASK_BIN" != "x" ] && PATH=$DAMASK_BIN:$PATH
 
 SOLVER=`which DAMASK_spectral || True 2>/dev/null`
 PROCESSING=`which postResults || True 2>/dev/null`
-if [ "x$DAMASK_NUM_THREADS" = "x" ]; then
-  DAMASK_NUM_THREADS=1
-fi
+[ "x$DAMASK_NUM_THREADS" = "x" ] && DAMASK_NUM_THREADS=1
 
 # according to http://software.intel.com/en-us/forums/topic/501500
 # this seems to make sense for the stack size
diff --git a/examples/ConfigFiles/Kinematics_Thermal_Expansion.config b/examples/ConfigFiles/Kinematics_Thermal_Expansion.config
index 6c778136f..0cedff049 100644
--- a/examples/ConfigFiles/Kinematics_Thermal_Expansion.config
+++ b/examples/ConfigFiles/Kinematics_Thermal_Expansion.config
@@ -1,2 +1,2 @@
 (kinematics) thermal_expansion
-thermal_expansion11         0.00231
+thermal_expansion11         0.0000231
diff --git a/examples/ConfigFiles/Phase_Dislotwin_TWIP-Steel-FeMnC.config b/examples/ConfigFiles/Phase_Dislotwin_TWIP-Steel-FeMnC.config
index e7c9d4e19..0330e6f4a 100644
--- a/examples/ConfigFiles/Phase_Dislotwin_TWIP-Steel-FeMnC.config
+++ b/examples/ConfigFiles/Phase_Dislotwin_TWIP-Steel-FeMnC.config
@@ -43,7 +43,7 @@ q_slip              1.0             # q-exponent in glide velocity
 CLambdaSlip         10.0            # Adj. parameter controlling dislocation mean free path
 D0                  4.0e-5          # Vacancy diffusion prefactor [m**2/s]
 Qsd                 4.5e-19         # Activation energy for climb [J]
-Catomicvolume       1.0             # Adj. parameter controlling the atomic volume [in b]
+Catomicvolume       1.0             # Adj. parameter controlling the atomic volume [in b^3]
 Cedgedipmindistance 1.0             # Adj. parameter controlling the minimum dipole distance [in b]
 atol_rho            1.0
 interactionSlipSlip 0.122 0.122 0.625 0.07 0.137 0.122 # Interaction coefficients (Kubin et al. 2008)
diff --git a/installation/mods_MarcMentat/apply_DAMASK_modifications.sh b/installation/mods_MarcMentat/apply_DAMASK_modifications.sh
index d1c1e37dd..7eef9c729 100755
--- a/installation/mods_MarcMentat/apply_DAMASK_modifications.sh
+++ b/installation/mods_MarcMentat/apply_DAMASK_modifications.sh
@@ -102,8 +102,7 @@ for filename in 'edit_window' \
 done
 
 # Mentat scripts
-echo ''
-echo 'adapting Mentat menus...'
+echo -e '\nadapting Mentat menus...'
 theDIR=$INSTALLDIR/mentat$VERSION/menus
 for filename in 'job_run.ms'; do
   cp $SCRIPTLOCATION/$VERSION/Mentat_menus/$filename $theDIR
@@ -163,5 +162,4 @@ echo ''
 echo 'precompiling $VERSION HYPELA2 user subroutine...'
 echo 'not yet implemented..!'
 
-echo ''
-echo 'done.'
+echo -e '\ndone.'
diff --git a/installation/symlink_Processing.py b/installation/symlink_Processing.py
index d10b5af55..50ea86e68 100755
--- a/installation/symlink_Processing.py
+++ b/installation/symlink_Processing.py
@@ -30,7 +30,7 @@ for subDir in processing_subDirs:
 
   for theFile in os.listdir(theDir):
     theName,theExt = os.path.splitext(theFile)
-    if theExt in processing_extensions:                           # only consider files with proper extensions
+    if theExt in processing_extensions:                                                             # only consider files with proper extensions
 
       src      = os.path.abspath(os.path.join(theDir,theFile))
       sym_link = os.path.abspath(os.path.join(binDir,theName))
diff --git a/processing/post/addGaussian.py b/processing/post/addGaussian.py
new file mode 100755
index 000000000..bc5599d49
--- /dev/null
+++ b/processing/post/addGaussian.py
@@ -0,0 +1,137 @@
+#!/usr/bin/env python2.7
+# -*- coding: UTF-8 no BOM -*-
+
+import os,sys
+import numpy as np
+from optparse import OptionParser
+from scipy import ndimage
+import damask
+
+scriptName = os.path.splitext(os.path.basename(__file__))[0]
+scriptID   = ' '.join([scriptName,damask.version])
+
+
+# --------------------------------------------------------------------
+#                                MAIN
+# --------------------------------------------------------------------
+
+parser = OptionParser(option_class=damask.extendableOption, usage='%prog option(s) [ASCIItable(s)]', description = """
+Add column(s) containing Gaussian filtered values of requested column(s).
+Operates on periodic and non-periodic ordered three-dimensional data sets.
+For Details see scipy.ndimage documentation.
+
+""", version = scriptID)
+
+parser.add_option('-p','--pos','--periodiccellcenter',
+                  dest = 'pos',
+                  type = 'string', metavar = 'string',
+                  help = 'label of coordinates [%default]')
+parser.add_option('-s','--scalar',
+                  dest = 'scalar',
+                  action = 'extend', metavar = '<string LIST>',
+                  help = 'label(s) of scalar field values')
+parser.add_option('-o','--order',
+                  dest = 'order',
+                  type = int,
+                  metavar = 'int',
+                  help = 'order of the filter')
+parser.add_option('--sigma',
+                  dest = 'sigma',
+                  type = float,
+                  metavar = 'float',
+                  help = 'standard deviation')
+parser.add_option('--periodic',
+                  dest = 'periodic',
+                  action = 'store_true',
+                  help = 'assume periodic grain structure'
+                 )
+
+
+
+parser.set_defaults(pos = 'pos',
+                    order = 0,
+                    sigma = 1,
+                    periodic = False
+                   )
+
+(options,filenames) = parser.parse_args()
+
+if options.scalar is None:
+  parser.error('no data column specified.')
+
+# --- loop over input files ------------------------------------------------------------------------
+
+if filenames == []: filenames = [None]
+
+for name in filenames:
+  try:    table = damask.ASCIItable(name = name,buffered = False)
+  except: continue
+  damask.util.report(scriptName,name)
+
+# ------------------------------------------ read header ------------------------------------------
+
+  table.head_read()
+
+# ------------------------------------------ sanity checks ----------------------------------------
+
+  items = {
+            'scalar': {'dim': 1, 'shape': [1], 'labels':options.scalar, 'active':[], 'column': []},
+          }
+  errors  = []
+  remarks = []
+  column = {}
+  
+  if table.label_dimension(options.pos) != 3: errors.append('coordinates {} are not a vector.'.format(options.pos))
+  else: colCoord = table.label_index(options.pos)
+
+  for type, data in items.iteritems():
+    for what in (data['labels'] if data['labels'] is not None else []):
+      dim = table.label_dimension(what)
+      if dim != data['dim']: remarks.append('column {} is not a {}.'.format(what,type))
+      else:
+        items[type]['active'].append(what)
+        items[type]['column'].append(table.label_index(what))
+
+  if remarks != []: damask.util.croak(remarks)
+  if errors  != []:
+    damask.util.croak(errors)
+    table.close(dismiss = True)
+    continue
+
+# ------------------------------------------ assemble header --------------------------------------
+
+  table.info_append(scriptID + '\t' + ' '.join(sys.argv[1:]))
+  for type, data in items.iteritems():
+    for label in data['active']:
+      table.labels_append(['Gauss{}({})'.format(options.sigma,label)])                              # extend ASCII header with new labels
+  table.head_write()
+
+# --------------- figure out size and grid ---------------------------------------------------------
+
+  table.data_readArray()
+
+  coords = [np.unique(table.data[:,colCoord+i]) for i in range(3)]
+  mincorner = np.array(map(min,coords))
+  maxcorner = np.array(map(max,coords))
+  grid   = np.array(map(len,coords),'i')
+  size   = grid/np.maximum(np.ones(3,'d'), grid-1.0) * (maxcorner-mincorner)                        # size from edge to edge = dim * n/(n-1) 
+  size   = np.where(grid > 1, size, min(size[grid > 1]/grid[grid > 1]))                             # spacing for grid==1 equal to smallest among other ones
+
+# ------------------------------------------ process value field -----------------------------------
+
+  stack = [table.data]
+  for type, data in items.iteritems():
+    for i,label in enumerate(data['active']):
+      stack.append(ndimage.filters.gaussian_filter(table.data[:,data['column'][i]],
+                                                   options.sigma,options.order,
+                                                   mode = 'wrap' if options.periodic else 'nearest'
+                                                   ).reshape([table.data.shape[0],1])
+                   )
+
+# ------------------------------------------ output result -----------------------------------------
+  if len(stack) > 1: table.data = np.hstack(tuple(stack))
+  table.data_writeArray('%.12g')
+
+# ------------------------------------------ output finalization -----------------------------------
+
+  table.close()                                                                                     # close input ASCII table (works for stdin)
diff --git a/processing/post/addGrainID.py b/processing/post/addGrainID.py
index c39d67983..45034034b 100755
--- a/processing/post/addGrainID.py
+++ b/processing/post/addGrainID.py
@@ -148,7 +148,7 @@ for name in filenames:
 
   bg.set_message('reading positions...')
 
-  table.data_readArray(options.pos)                                                              # read position vectors
+  table.data_readArray(options.pos)                                                                 # read position vectors
   grainID = -np.ones(len(table.data),dtype=int)
 
   start = tick = time.clock()
diff --git a/processing/pre/geom_fromTable.py b/processing/pre/geom_fromTable.py
index 9af92888d..b10bc9f88 100755
--- a/processing/pre/geom_fromTable.py
+++ b/processing/pre/geom_fromTable.py
@@ -23,7 +23,7 @@ Generate geometry description and material configuration from position, phase, a
 parser.add_option('--coordinates',
                   dest = 'pos',
                   type = 'string', metavar = 'string',
-                  help = 'coordinates label')
+                  help = 'coordinates label (%default)')
 parser.add_option('--phase',
                   dest = 'phase',
                   type = 'string', metavar = 'string',
@@ -90,6 +90,7 @@ parser.set_defaults(symmetry       = [damask.Symmetry.lattices[-1]],
                     homogenization = 1,
                     crystallite    = 1,
                     verbose        = False,
+                    pos            = 'pos',
                    )
 
 (options,filenames) = parser.parse_args()
diff --git a/processing/pre/mentat_pbcOnBoxMesh.py b/processing/pre/mentat_pbcOnBoxMesh.py
old mode 100644
new mode 100755
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index 83647d8d0..eb1448028 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -1,16 +1,5 @@
 # special flags for some files
 if (${CMAKE_Fortran_COMPILER_ID} STREQUAL "GNU")
-
-    SET_SOURCE_FILES_PROPERTIES( "prec.f90" PROPERTIES
-                                    COMPILE_FLAGS "-fno-range-check -fall-intrinsics -fno-fast-math")
-    # fno-range-check:  Disable range checking on results of simplification of constant expressions during compilation
-    #      -->   allows the definition of DAMASK_NaN
-    #-fall-intrinsics:  all intrinsic procedures (including the GNU-specific extensions) are accepted. -Wintrinsics-std will be ignored
-    #                   and no user-defined procedure with the same name as any intrinsic will be called except when it is explicitly declared external
-    #      -->   allows the use of 'isnan'
-    #-fno-fast-math:
-    #      -->   otherwise, when setting -ffast-math, isnan always evaluates to false (I would call it a bug)
-
     SET_SOURCE_FILES_PROPERTIES( "lattice.f90" PROPERTIES
                                 COMPILE_FLAGS "-ffree-line-length-240")
     # long lines for interaction matrix
@@ -185,8 +174,10 @@ if ("${PROJECT_NAME}" STREQUAL "DAMASK_spectral")
       list(APPEND OBJECTFILES $<TARGET_OBJECTS:SPECTRAL_SOLVER>)
       if(NOT "${CMAKE_BUILD_TYPE}" STREQUAL "SYNTAXONLY")
         add_executable(DAMASK_spectral "DAMASK_spectral.f90" ${OBJECTFILES})
-        add_dependencies(DAMASK_spectral SPECTRAL_SOLVER)
+      else()
+        add_library(DAMASK_spectral OBJECT "DAMASK_spectral.f90")
       endif()
+      add_dependencies(DAMASK_spectral SPECTRAL_SOLVER)
 elseif ("${PROJECT_NAME}" STREQUAL "DAMASK_FEM")
       add_library(FEM_UTILITIES OBJECT "FEM_utilities.f90")
       add_dependencies(FEM_UTILITIES DAMASK_CPFE)
diff --git a/src/DAMASK_spectral.f90 b/src/DAMASK_spectral.f90
old mode 100644
new mode 100755
index 6f1794e35..dfa1746b2
--- a/src/DAMASK_spectral.f90
+++ b/src/DAMASK_spectral.f90
@@ -364,12 +364,12 @@ program DAMASK_spectral
          case (DAMASK_spectral_SolverBasicPETSc_label)
            call basicPETSc_init
          case (DAMASK_spectral_SolverAL_label)
-           if(iand(debug_level(debug_spectral),debug_levelBasic)/= 0 .and. worldrank == 0_pInt) &
+           if(iand(debug_level(debug_spectral),debug_levelBasic)/= 0) &
            call IO_warning(42_pInt, ext_msg='debug Divergence')
            call AL_init
 
          case (DAMASK_spectral_SolverPolarisation_label)
-           if(iand(debug_level(debug_spectral),debug_levelBasic)/= 0 .and. worldrank == 0_pInt) &
+           if(iand(debug_level(debug_spectral),debug_levelBasic)/= 0) &
            call IO_warning(42_pInt, ext_msg='debug Divergence')
            call Polarisation_init
 
@@ -450,8 +450,7 @@ program DAMASK_spectral
      if (ierr /= 0_pInt) call IO_error(error_ID=894_pInt, ext_msg='MPI_file_write')
    enddo
    fileOffset = fileOffset + sum(outputSize)                                                        ! forward to current file position
-   if (worldrank == 0) &
-     write(6,'(1/,a)') ' ... writing initial configuration to file ........................'
+   write(6,'(1/,a)') ' ... writing initial configuration to file ........................'
  endif
 !--------------------------------------------------------------------------------------------------
 ! loopping over loadcases
@@ -498,22 +497,20 @@ program DAMASK_spectral
 
 !--------------------------------------------------------------------------------------------------
 ! report begin of new increment
-         if (worldrank == 0) then
-           write(6,'(/,a)') ' ###########################################################################'
-           write(6,'(1x,a,es12.5'//&
-                   ',a,'//IO_intOut(inc)//',a,'//IO_intOut(loadCases(currentLoadCase)%incs)//&
-                   ',a,'//IO_intOut(stepFraction)//',a,'//IO_intOut(subStepFactor**cutBackLevel)//&
-                   ',a,'//IO_intOut(currentLoadCase)//',a,'//IO_intOut(size(loadCases))//')') &
-                   'Time', time, &
-                   's: Increment ', inc, '/', loadCases(currentLoadCase)%incs,&
-                   '-', stepFraction, '/', subStepFactor**cutBackLevel,&
-                   ' of load case ', currentLoadCase,'/',size(loadCases)
-           flush(6)
-           write(incInfo,'(a,'//IO_intOut(totalIncsCounter)//',a,'//IO_intOut(sum(loadCases%incs))//&
-                 ',a,'//IO_intOut(stepFraction)//',a,'//IO_intOut(subStepFactor**cutBackLevel)//')') &
-                 'Increment ',totalIncsCounter,'/',sum(loadCases%incs),&
-                 '-',stepFraction, '/', subStepFactor**cutBackLevel
-         endif
+         write(6,'(/,a)') ' ###########################################################################'
+         write(6,'(1x,a,es12.5'//&
+                 ',a,'//IO_intOut(inc)//',a,'//IO_intOut(loadCases(currentLoadCase)%incs)//&
+                 ',a,'//IO_intOut(stepFraction)//',a,'//IO_intOut(subStepFactor**cutBackLevel)//&
+                 ',a,'//IO_intOut(currentLoadCase)//',a,'//IO_intOut(size(loadCases))//')') &
+                 'Time', time, &
+                 's: Increment ', inc, '/', loadCases(currentLoadCase)%incs,&
+                 '-', stepFraction, '/', subStepFactor**cutBackLevel,&
+                 ' of load case ', currentLoadCase,'/',size(loadCases)
+         flush(6)
+         write(incInfo,'(a,'//IO_intOut(totalIncsCounter)//',a,'//IO_intOut(sum(loadCases%incs))//&
+               ',a,'//IO_intOut(stepFraction)//',a,'//IO_intOut(subStepFactor**cutBackLevel)//')') &
+               'Increment ',totalIncsCounter,'/',sum(loadCases%incs),&
+               '-',stepFraction, '/', subStepFactor**cutBackLevel
 
 !--------------------------------------------------------------------------------------------------
 ! forward fields
@@ -595,7 +592,7 @@ program DAMASK_spectral
          cutBack = .False.
          if(solres(1)%termIll .or. .not. all(solres(:)%converged .and. solres(:)%stagConverged)) then ! no solution found
            if (cutBackLevel < maxCutBack) then                                                      ! do cut back
-             if (worldrank == 0) write(6,'(/,a)') ' cut back detected'
+             write(6,'(/,a)') ' cut back detected'
              cutBack = .True.
              stepFraction = (stepFraction - 1_pInt) * subStepFactor                                 ! adjust to new denominator
              cutBackLevel = cutBackLevel + 1_pInt
@@ -603,11 +600,15 @@ program DAMASK_spectral
              timeinc = timeinc/2.0_pReal
            elseif (solres(1)%termIll) then                                                          ! material point model cannot find a solution, exit in any casy
              call IO_warning(850_pInt)
+             call MPI_file_close(resUnit,ierr)
+             close(statUnit)
              call quit(-1_pInt*(lastRestartWritten+1_pInt))                                         ! quit and provide information about last restart inc written
            elseif (continueCalculation == 1_pInt)  then
              guess = .true.                                                                         ! accept non converged BVP solution
            else                                                                                     ! default behavior, exit if spectral solver does not converge
              call IO_warning(850_pInt)
+             call MPI_file_close(resUnit,ierr)
+             close(statUnit)
              call quit(-1_pInt*(lastRestartWritten+1_pInt))                                         ! quit and provide information about last restart inc written
            endif
          else
@@ -624,13 +625,11 @@ program DAMASK_spectral
        cutBackLevel = max(0_pInt, cutBackLevel - 1_pInt)                                            ! try half number of subincs next inc
        if(all(solres(:)%converged)) then                                                            ! report converged inc
          convergedCounter = convergedCounter + 1_pInt
-         if (worldrank == 0) &
-           write(6,'(/,a,'//IO_intOut(totalIncsCounter)//',a)') &
-                                     ' increment ', totalIncsCounter, ' converged'
+         write(6,'(/,a,'//IO_intOut(totalIncsCounter)//',a)') &
+                                   ' increment ', totalIncsCounter, ' converged'
        else
-         if (worldrank == 0) &
-           write(6,'(/,a,'//IO_intOut(totalIncsCounter)//',a)') &                                   ! report non-converged inc
-                                     ' increment ', totalIncsCounter, ' NOT converged'
+         write(6,'(/,a,'//IO_intOut(totalIncsCounter)//',a)') &                                   ! report non-converged inc
+                                   ' increment ', totalIncsCounter, ' NOT converged'
          notConvergedCounter = notConvergedCounter + 1_pInt
        endif; flush(6)
        if (mod(inc,loadCases(currentLoadCase)%outputFrequency) == 0_pInt) then                      ! at output frequency
@@ -665,45 +664,15 @@ program DAMASK_spectral
 
 !--------------------------------------------------------------------------------------------------
 ! report summary of whole calculation
- if (worldrank == 0) then
-   write(6,'(/,a)') ' ###########################################################################'
-   write(6,'(1x,i6.6,a,i6.6,a,f5.1,a)') convergedCounter, ' out of ', &
-                                     notConvergedCounter + convergedCounter, ' (', &
-                                     real(convergedCounter, pReal)/&
-                                     real(notConvergedCounter + convergedCounter,pReal)*100.0_pReal, &
-                                     ' %) increments converged!'
- endif
+ write(6,'(/,a)') ' ###########################################################################'
+ write(6,'(1x,i6.6,a,i6.6,a,f5.1,a)') convergedCounter, ' out of ', &
+                                   notConvergedCounter + convergedCounter, ' (', &
+                                   real(convergedCounter, pReal)/&
+                                   real(notConvergedCounter + convergedCounter,pReal)*100.0_pReal, &
+                                   ' %) increments converged!'
  call MPI_file_close(resUnit,ierr)
  close(statUnit)
 
- do field = 1, nActiveFields
-   select case(loadCases(1)%ID(field))
-     case(FIELD_MECH_ID)
-       select case (spectral_solver)
-         case (DAMASK_spectral_SolverBasicPETSc_label)
-           call BasicPETSC_destroy()
-         case (DAMASK_spectral_SolverAL_label)
-           call AL_destroy()
-         case (DAMASK_spectral_SolverPolarisation_label)
-           call Polarisation_destroy()
-       end select
-     case(FIELD_THERMAL_ID)
-       call spectral_thermal_destroy()
-     case(FIELD_DAMAGE_ID)
-       call spectral_damage_destroy()
-   end select
- enddo
- call utilities_destroy()
-
- call PETScFinalize(ierr); CHKERRQ(ierr)
-
-#ifdef _OPENMP
-  call MPI_finalize(i)
-  if (i /= 0_pInt) then 
-    call IO_error(error_ID=894, el=i, ext_msg="Finalize()")
-  endif
-#endif
-
  if (notConvergedCounter > 0_pInt) call quit(3_pInt)                                                ! error if some are not converged
  call quit(0_pInt)                                                                                  ! no complains ;)
 
@@ -721,11 +690,47 @@ end program DAMASK_spectral
 subroutine quit(stop_id)
  use prec, only: &
    pInt
+ use spectral_mech_Basic, only: &
+   BasicPETSC_destroy
+ use spectral_mech_AL, only: &
+   AL_destroy
+ use spectral_mech_Polarisation, only: &
+   Polarisation_destroy
+ use spectral_damage, only: &
+   spectral_damage_destroy
+ use spectral_thermal, only: &
+   spectral_thermal_destroy
+ use spectral_utilities, only: &
+   utilities_destroy
 
  implicit none
+ 
+#include <petsc/finclude/petscsys.h>
  integer(pInt), intent(in) :: stop_id
  integer, dimension(8) :: dateAndTime                                                               ! type default integer
+ integer(pInt) :: error = 0_pInt
+ PetscErrorCode :: ierr = 0
+ logical :: ErrorInQuit
+ 
+ external :: &
+   PETScFinalize, &
+   MPI_finalize
 
+ call BasicPETSC_destroy()
+ call AL_destroy()
+ call Polarisation_destroy()
+ call spectral_damage_destroy()
+ call spectral_thermal_destroy()
+ call utilities_destroy()
+
+ call PETScFinalize(ierr)
+ 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'
+#endif
+ ErrorInQuit = (ierr /= 0 .or. error /= 0_pInt)
+ 
  call date_and_time(values = dateAndTime)
  write(6,'(/,a)') 'DAMASK terminated on:'
  write(6,'(a,2(i2.2,a),i4.4)') 'Date:               ',dateAndTime(3),'/',&
@@ -735,12 +740,13 @@ subroutine quit(stop_id)
                                                       dateAndTime(6),':',&
                                                       dateAndTime(7)
 
- if (stop_id == 0_pInt) stop 0                                                                      ! normal termination
- if (stop_id <  0_pInt) then                                                                        ! terminally ill, restart might help
+ 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
    write(0,'(a,i6)') 'restart information available at ', stop_id*(-1_pInt)
    stop 2
  endif
- if (stop_id == 3_pInt) stop 3                                                                      ! not all incs converged
+ if (stop_id == 3_pInt .and. .not. ErrorInQuit) stop 3                                              ! not all incs converged
+ 
  stop 1                                                                                             ! error (message from IO_error)
 
 end subroutine quit
diff --git a/src/crystallite.f90 b/src/crystallite.f90
index efeab9dcd..c5bd4d979 100644
--- a/src/crystallite.f90
+++ b/src/crystallite.f90
@@ -1232,8 +1232,8 @@ end subroutine crystallite_stressAndItsTangent
 !> @brief integrate stress, state with 4th order explicit Runge Kutta method
 !--------------------------------------------------------------------------------------------------
 subroutine crystallite_integrateStateRK4()
- use prec, only: &
-   prec_isNaN
+ use, intrinsic :: &
+   IEEE_arithmetic
  use numerics, only: &
    numerics_integrationMode
  use debug, only: &
@@ -1331,9 +1331,9 @@ subroutine crystallite_integrateStateRK4()
      if (crystallite_todo(g,i,e)) then
        c = phasememberAt(g,i,e)
        p = phaseAt(g,i,e)
-       NaN = any(prec_isNaN(plasticState(p)%dotState(:,c)))
+       NaN = any(IEEE_is_NaN(plasticState(p)%dotState(:,c)))
        do mySource = 1_pInt, phase_Nsources(p)
-         NaN = NaN .or. any(prec_isNaN(sourceState(p)%p(mySource)%dotState(:,c)))
+         NaN = NaN .or. any(IEEE_is_NaN(sourceState(p)%p(mySource)%dotState(:,c)))
        enddo
        if (NaN) then                                                                                    ! NaN occured in any dotState
          if (.not. crystallite_localPlasticity(g,i,e)) then                                             ! if broken non-local...
@@ -1475,9 +1475,9 @@ subroutine crystallite_integrateStateRK4()
 
            p = phaseAt(g,i,e)
            c = phasememberAt(g,i,e)
-           NaN = any(prec_isNaN(plasticState(p)%dotState(:,c)))
+           NaN = any(IEEE_is_NaN(plasticState(p)%dotState(:,c)))
            do mySource = 1_pInt, phase_Nsources(p)
-             NaN = NaN .or. any(prec_isNaN(sourceState(p)%p(mySource)%dotState(:,c)))
+             NaN = NaN .or. any(IEEE_is_NaN(sourceState(p)%p(mySource)%dotState(:,c)))
            enddo
            if (NaN) then                                                                                   ! NaN occured in any dotState
              if (.not. crystallite_localPlasticity(g,i,e)) then                                            ! if broken non-local...
@@ -1528,8 +1528,8 @@ end subroutine crystallite_integrateStateRK4
 !> adaptive step size  (use 5th order solution to advance = "local extrapolation")
 !--------------------------------------------------------------------------------------------------
 subroutine crystallite_integrateStateRKCK45()
- use prec, only: &
-   prec_isNaN
+ use, intrinsic :: &
+   IEEE_arithmetic
  use debug, only: &
    debug_level, &
    debug_crystallite, &
@@ -1647,9 +1647,9 @@ subroutine crystallite_integrateStateRKCK45()
      if (crystallite_todo(g,i,e)) then
        cc = phasememberAt(g,i,e)
        p = phaseAt(g,i,e)
-       NaN = any(prec_isNaN(plasticState(p)%dotState(:,cc)))
+       NaN = any(IEEE_is_NaN(plasticState(p)%dotState(:,cc)))
        do mySource = 1_pInt, phase_Nsources(p)
-         NaN = NaN .or. any(prec_isNaN(sourceState(p)%p(mySource)%dotState(:,cc)))
+         NaN = NaN .or. any(IEEE_is_NaN(sourceState(p)%p(mySource)%dotState(:,cc)))
        enddo
        if (NaN) then                                                                                       ! NaN occured in any dotState
          if (.not. crystallite_localPlasticity(g,i,e)) then                                                ! if broken non-local...
@@ -1801,9 +1801,9 @@ subroutine crystallite_integrateStateRKCK45()
 
          p = phaseAt(g,i,e)
          cc = phasememberAt(g,i,e)
-         NaN = any(prec_isNaN(plasticState(p)%dotState(:,cc)))
+         NaN = any(IEEE_is_NaN(plasticState(p)%dotState(:,cc)))
          do mySource = 1_pInt, phase_Nsources(p)
-           NaN = NaN .or. any(prec_isNaN(sourceState(p)%p(mySource)%dotState(:,cc)))
+           NaN = NaN .or. any(IEEE_is_NaN(sourceState(p)%p(mySource)%dotState(:,cc)))
          enddo
          if (NaN) then                                                                                   ! NaN occured in any dotState
            if (.not. crystallite_localPlasticity(g,i,e)) then                                              ! if broken non-local...
@@ -2031,8 +2031,8 @@ end subroutine crystallite_integrateStateRKCK45
 !> @brief integrate stress, state with 1st order Euler method with adaptive step size
 !--------------------------------------------------------------------------------------------------
 subroutine crystallite_integrateStateAdaptiveEuler()
- use prec, only: &
-   prec_isNaN
+ use, intrinsic :: &
+   IEEE_arithmetic
  use debug, only: &
    debug_level, &
    debug_crystallite, &
@@ -2133,9 +2133,9 @@ subroutine crystallite_integrateStateAdaptiveEuler()
        if (crystallite_todo(g,i,e)) then
          p = phaseAt(g,i,e)
          c = phasememberAt(g,i,e)
-         NaN = any(prec_isNaN(plasticState(p)%dotState(:,c)))
+         NaN = any(IEEE_is_NaN(plasticState(p)%dotState(:,c)))
          do mySource = 1_pInt, phase_Nsources(p)
-           NaN = NaN .or. any(prec_isNaN(sourceState(p)%p(mySource)%dotState(:,c)))
+           NaN = NaN .or. any(IEEE_is_NaN(sourceState(p)%p(mySource)%dotState(:,c)))
          enddo
          if (NaN) then                                                                                       ! NaN occured in any dotState
            if (.not. crystallite_localPlasticity(g,i,e)) then                                                ! if broken non-local...
@@ -2254,9 +2254,9 @@ subroutine crystallite_integrateStateAdaptiveEuler()
        if (crystallite_todo(g,i,e)) then
          p = phaseAt(g,i,e)
          c = phasememberAt(g,i,e)
-         NaN = any(prec_isNaN(plasticState(p)%dotState(:,c)))
+         NaN = any(IEEE_is_NaN(plasticState(p)%dotState(:,c)))
          do mySource = 1_pInt, phase_Nsources(p)
-           NaN = NaN .or. any(prec_isNaN(sourceState(p)%p(mySource)%dotState(:,c)))
+           NaN = NaN .or. any(IEEE_is_NaN(sourceState(p)%p(mySource)%dotState(:,c)))
          enddo
          if (NaN) then                                                                                     ! NaN occured in any dotState
            if (.not. crystallite_localPlasticity(g,i,e)) then                                              ! if broken non-local...
@@ -2391,8 +2391,8 @@ end subroutine crystallite_integrateStateAdaptiveEuler
 !> @brief integrate stress, and state with 1st order explicit Euler method
 !--------------------------------------------------------------------------------------------------
 subroutine crystallite_integrateStateEuler()
- use prec, only: &
-   prec_isNaN
+ use, intrinsic :: &
+   IEEE_arithmetic
  use debug, only: &
    debug_level, &
    debug_crystallite, &
@@ -2471,9 +2471,9 @@ eIter = FEsolving_execElem(1:2)
        if (crystallite_todo(g,i,e) .and. .not. crystallite_converged(g,i,e)) then
          c = phasememberAt(g,i,e)
          p = phaseAt(g,i,e)
-         NaN = any(prec_isNaN(plasticState(p)%dotState(:,c)))
+         NaN = any(IEEE_is_NaN(plasticState(p)%dotState(:,c)))
          do mySource = 1_pInt, phase_Nsources(p)
-           NaN = NaN .or. any(prec_isNaN(sourceState(p)%p(mySource)%dotState(:,c)))
+           NaN = NaN .or. any(IEEE_is_NaN(sourceState(p)%p(mySource)%dotState(:,c)))
          enddo
          if (NaN) then                                                                                       ! NaN occured in any dotState
            if (.not. crystallite_localPlasticity(g,i,e) .and. .not. numerics_timeSyncing) then               ! if broken non-local...
@@ -2614,8 +2614,8 @@ end subroutine crystallite_integrateStateEuler
 !> using Fixed Point Iteration to adapt the stepsize
 !--------------------------------------------------------------------------------------------------
 subroutine crystallite_integrateStateFPI()
- use prec, only: &
-   prec_isNaN
+ use, intrinsic :: &
+   IEEE_arithmetic
  use debug, only: &
    debug_e, &
    debug_i, &
@@ -2737,9 +2737,9 @@ subroutine crystallite_integrateStateFPI()
      if (crystallite_todo(g,i,e)) then
        p = phaseAt(g,i,e)
        c = phasememberAt(g,i,e)
-       NaN = any(prec_isNaN(plasticState(p)%dotState(:,c)))
+       NaN = any(IEEE_is_NaN(plasticState(p)%dotState(:,c)))
        do mySource = 1_pInt, phase_Nsources(p)
-         NaN = NaN .or. any(prec_isNaN(sourceState(p)%p(mySource)%dotState(:,c)))
+         NaN = NaN .or. any(IEEE_is_NaN(sourceState(p)%p(mySource)%dotState(:,c)))
        enddo
        if (NaN) then                                                                                       ! NaN occured in any dotState
          if (.not. crystallite_localPlasticity(g,i,e)) then                                                ! if broken is a non-local...
@@ -2851,9 +2851,9 @@ subroutine crystallite_integrateStateFPI()
        if (crystallite_todo(g,i,e) .and. .not. crystallite_converged(g,i,e)) then
          p = phaseAt(g,i,e)
          c = phasememberAt(g,i,e)
-         NaN = any(prec_isNaN(plasticState(p)%dotState(:,c)))
+         NaN = any(IEEE_is_NaN(plasticState(p)%dotState(:,c)))
          do mySource = 1_pInt, phase_Nsources(p)
-           NaN = NaN .or. any(prec_isNaN(sourceState(p)%p(mySource)%dotState(:,c)))
+           NaN = NaN .or. any(IEEE_is_NaN(sourceState(p)%p(mySource)%dotState(:,c)))
          enddo
          if (NaN) then                                                                                       ! NaN occured in any dotState
            crystallite_todo(g,i,e) = .false.                                                                 ! ... skip me next time
@@ -3062,8 +3062,9 @@ end subroutine crystallite_integrateStateFPI
 !> returns true, if state jump was successfull or not needed. false indicates NaN in delta state
 !--------------------------------------------------------------------------------------------------
 logical function crystallite_stateJump(ipc,ip,el)
+ use, intrinsic :: &
+   IEEE_arithmetic
  use prec, only: &
-   prec_isNaN, &
    dNeq0
  use debug, only: &
    debug_level, &
@@ -3098,7 +3099,7 @@ logical function crystallite_stateJump(ipc,ip,el)
  p = phaseAt(ipc,ip,el)
  call constitutive_collectDeltaState(crystallite_Tstar_v(1:6,ipc,ip,el), crystallite_Fe(1:3,1:3,ipc,ip,el), ipc,ip,el)
  mySizePlasticDeltaState = plasticState(p)%sizeDeltaState
- if( any(prec_isNaN(plasticState(p)%deltaState(:,c)))) then                                         ! NaN occured in deltaState
+ if( any(IEEE_is_NaN(plasticState(p)%deltaState(:,c)))) then                                       ! NaN occured in deltaState
    crystallite_stateJump = .false.
    return
  endif
@@ -3106,7 +3107,7 @@ logical function crystallite_stateJump(ipc,ip,el)
                                                       plasticState(p)%deltaState(1:mySizePlasticDeltaState,c)
  do mySource = 1_pInt, phase_Nsources(p)
    mySizeSourceDeltaState = sourceState(p)%p(mySource)%sizeDeltaState
-   if( any(prec_isNaN(sourceState(p)%p(mySource)%deltaState(:,c)))) then                            ! NaN occured in deltaState
+   if( any(IEEE_is_NaN(sourceState(p)%p(mySource)%deltaState(:,c)))) then                          ! NaN occured in deltaState
      crystallite_stateJump = .false.
      return
    endif
@@ -3169,15 +3170,18 @@ logical function crystallite_integrateStress(&
       el,&          ! element number
       timeFraction &
       )
+ use, intrinsic :: &
+   IEEE_arithmetic
  use prec, only:         pLongInt, &
                          tol_math_check, &
-                         prec_isNaN, &
                          dEq0
  use numerics, only:     nStress, &
                          aTol_crystalliteStress, &
                          rTol_crystalliteStress, &
                          iJacoLpresiduum, &
-                         numerics_integrationMode
+                         numerics_integrationMode, &
+                         subStepSizeLp, &
+                         subStepSizeLi
  use debug, only:        debug_level, &
                          debug_crystallite, &
                          debug_levelBasic, &
@@ -3263,9 +3267,7 @@ logical function crystallite_integrateStress(&
                                      dLp_dT3333, &
                                      dLi_dT3333
  real(pReal)                         detInvFi, &                                                     ! determinant of InvFi
-                                     steplengthLp0, &
                                      steplengthLp, &
-                                     steplengthLi0, &
                                      steplengthLi, &
                                      dt, &                                                           ! time increment
                                      aTolLp, &
@@ -3351,8 +3353,7 @@ logical function crystallite_integrateStress(&
 
  NiterationStressLi = 0_pInt
  jacoCounterLi      = 0_pInt
- steplengthLi0      = 1.0_pReal
- steplengthLi       = steplengthLi0
+ steplengthLi       = 1.0_pReal
  residuumLi_old     = 0.0_pReal
 
  LiLoop: do
@@ -3372,8 +3373,7 @@ logical function crystallite_integrateStress(&
 
    NiterationStressLp = 0_pInt
    jacoCounterLp      = 0_pInt
-   steplengthLp0      = 1.0_pReal
-   steplengthLp       = steplengthLp0
+   steplengthLp       = 1.0_pReal
    residuumLp_old     = 0.0_pReal
    Lpguess_old        = Lpguess
 
@@ -3430,7 +3430,7 @@ logical function crystallite_integrateStress(&
                   aTol_crystalliteStress)                                                            ! minimum lower cutoff
      residuumLp = Lpguess - Lp_constitutive
 
-     if (any(prec_isNaN(residuumLp))) then                                                           ! NaN in residuum...
+     if (any(IEEE_is_NaN(residuumLp))) then                                                         ! NaN in residuum...
 #ifndef _OPENMP
        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 ', &
@@ -3445,9 +3445,9 @@ logical function crystallite_integrateStress(&
              .or. norm2(residuumLp) < norm2(residuumLp_old)) then                                     ! not converged, but improved norm of residuum (always proceed in first iteration)...
        residuumLp_old = residuumLp                                                                    ! ...remember old values and...
        Lpguess_old    = Lpguess
-       steplengthLp   = steplengthLp0                                                                 ! ...proceed with normal step length (calculate new search direction)
+       steplengthLp   = 1.0_pReal                                                                     ! ...proceed with normal step length (calculate new search direction)
      else                                                                                             ! not converged and residuum not improved...
-       steplengthLp = 0.5_pReal * steplengthLp                                                        ! ...try with smaller step length in same direction
+       steplengthLp = subStepSizeLp * steplengthLp                                                    ! ...try with smaller step length in same direction
        Lpguess    = Lpguess_old + steplengthLp * deltaLp
        cycle LpLoop
      endif
@@ -3520,7 +3520,7 @@ logical function crystallite_integrateStress(&
    aTolLi = max(rTol_crystalliteStress * max(norm2(Liguess),norm2(Li_constitutive)), &              ! absolute tolerance from largest acceptable relative error
                 aTol_crystalliteStress)                                                             ! minimum lower cutoff
    residuumLi = Liguess - Li_constitutive
-   if (any(prec_isNaN(residuumLi))) then                                                            ! NaN in residuum...
+   if (any(IEEE_is_NaN(residuumLi))) then                                                           ! NaN in residuum...
      return                                                                                         ! ...me = .false. to inform integrator about problem
    elseif (norm2(residuumLi) < aTolLi) then                                                         ! converged if below absolute tolerance
      exit LiLoop                                                                                    ! ...leave iteration loop
@@ -3528,9 +3528,9 @@ logical function crystallite_integrateStress(&
            .or. norm2(residuumLi) < norm2(residuumLi_old)) then                                     ! not converged, but improved norm of residuum (always proceed in first iteration)...
      residuumLi_old = residuumLi                                                                    ! ...remember old values and...
      Liguess_old    = Liguess
-     steplengthLi   = steplengthLi0                                                                 ! ...proceed with normal step length (calculate new search direction)
+     steplengthLi   = 1.0_pReal                                                                     ! ...proceed with normal step length (calculate new search direction)
    else                                                                                             ! not converged and residuum not improved...
-     steplengthLi   = 0.5_pReal * steplengthLi                                                      ! ...try with smaller step length in same direction
+     steplengthLi   = subStepSizeLi * steplengthLi                                                  ! ...try with smaller step length in same direction
      Liguess        = Liguess_old + steplengthLi * deltaLi
      cycle LiLoop
    endif
diff --git a/src/homogenization.f90 b/src/homogenization.f90
index 34bc93af0..c8c5fad01 100644
--- a/src/homogenization.f90
+++ b/src/homogenization.f90
@@ -923,12 +923,11 @@ subroutine materialpoint_postResults
    material_phase, &
    homogenization_Ngrains, &
    microstructure_crystallite
- use constitutive, only: &
 #ifdef FEM
+ use constitutive, only: &
    constitutive_plasticity_maxSizePostResults, &
-   constitutive_source_maxSizePostResults, &
+   constitutive_source_maxSizePostResults
 #endif
-   constitutive_postResults
  use crystallite, only: &
 #ifdef FEM
    crystallite_maxSizePostResults, &
diff --git a/src/lattice.f90 b/src/lattice.f90
index 7726c772d..a970ed85a 100644
--- a/src/lattice.f90
+++ b/src/lattice.f90
@@ -2170,8 +2170,8 @@ pure function lattice_symmetrize33(struct,T33)
 !> @brief figures whether unit quat falls into stereographic standard triangle
 !--------------------------------------------------------------------------------------------------
 logical pure function lattice_qInSST(Q, struct)
- use prec, only: &
-   prec_isNaN
+ use, intrinsic :: &
+   IEEE_arithmetic
  use math, only: &
    math_qToRodrig
 
@@ -2181,7 +2181,7 @@ logical pure function lattice_qInSST(Q, struct)
  real(pReal), dimension(3) ::                  Rodrig                                               ! Rodrigues vector of Q
 
  Rodrig = math_qToRodrig(Q)
- if (any(prec_isNaN(Rodrig))) then
+ if (any(IEEE_is_NaN(Rodrig))) then
    lattice_qInSST = .false.
  else
    select case (struct)
diff --git a/src/math.f90 b/src/math.f90
index 9629401ca..3a0533ddf 100644
--- a/src/math.f90
+++ b/src/math.f90
@@ -1666,15 +1666,16 @@ end function math_qToEulerAxisAngle
 !> @brief Rodrigues vector (x, y, z) from unit quaternion (w+ix+jy+kz)
 !--------------------------------------------------------------------------------------------------
 pure function math_qToRodrig(Q)
+ use, intrinsic :: &
+   IEEE_arithmetic
  use prec, only: &
-   DAMASK_NaN, &
    tol_math_check
 
  implicit none
  real(pReal), dimension(4), intent(in) :: Q
  real(pReal), dimension(3) :: math_qToRodrig
 
- math_qToRodrig = merge(Q(2:4)/Q(1),DAMASK_NaN,abs(Q(1)) > tol_math_check)                          ! NaN for 180 deg since Rodrig is unbound
+ math_qToRodrig = merge(Q(2:4)/Q(1),IEEE_value(1.0_pReal,IEEE_quiet_NaN),abs(Q(1)) > tol_math_check)! NaN for 180 deg since Rodrig is unbound
 
 end function math_qToRodrig
 
@@ -2092,11 +2093,11 @@ end function math_rotationalPart33
 
 !--------------------------------------------------------------------------------------------------
 !> @brief Eigenvalues of symmetric matrix m
-! will return NaN  on error
+! will return NaN on error
 !--------------------------------------------------------------------------------------------------
 function math_eigenvaluesSym(m)
- use prec, only: &
-   DAMASK_NaN
+ use, intrinsic :: &
+   IEEE_arithmetic
  
  implicit none
  real(pReal), dimension(:,:),                  intent(in)  :: m
@@ -2109,7 +2110,7 @@ function math_eigenvaluesSym(m)
 
  vectors = m                                                                                         ! copy matrix to input (doubles as output) array
  call dsyev('N','U',size(m,1),vectors,size(m,1),math_eigenvaluesSym,work,(64+2)*size(m,1),info)
- if (info /= 0_pInt) math_eigenvaluesSym = DAMASK_NaN
+ if (info /= 0_pInt) math_eigenvaluesSym = IEEE_value(1.0_pReal,IEEE_quiet_NaN)
 
 end function math_eigenvaluesSym
 
@@ -2701,29 +2702,13 @@ pure function math_rotate_forward3333(tensor,rot_tensor)
 end function math_rotate_forward3333
 
 
-!--------------------------------------------------------------------------------------------------
-!> @brief calculate average of tensor field
-!--------------------------------------------------------------------------------------------------
-function math_tensorAvg(field)
- 
- implicit none
- real(pReal), dimension(3,3) :: math_tensorAvg
- real(pReal), intent(in), dimension(:,:,:,:,:) :: field
- real(pReal) :: wgt
- 
- wgt = 1.0_pReal/real(size(field,3)*size(field,4)*size(field,5), pReal) 
- math_tensorAvg = sum(sum(sum(field,dim=5),dim=4),dim=3)*wgt
-
-end function math_tensorAvg
-
-
 !--------------------------------------------------------------------------------------------------
 !> @brief limits a scalar value to a certain range (either one or two sided)
 ! Will return NaN if left > right
 !--------------------------------------------------------------------------------------------------
 real(pReal) pure function math_limit(a, left, right)
- use prec, only: &
-   DAMASK_NaN
+ use, intrinsic :: &
+   IEEE_arithmetic
 
  implicit none
  real(pReal), intent(in) :: a
@@ -2735,7 +2720,8 @@ real(pReal) pure function math_limit(a, left, right)
                         merge(right, huge(a), present(right)) &
                   )
 
- if (present(left) .and. present(right)) math_limit = merge (DAMASK_NaN,math_limit, left>right)
+ if (present(left) .and. present(right)) &
+   math_limit = merge (IEEE_value(1.0_pReal,IEEE_quiet_NaN),math_limit, left>right)
 
 end function math_limit
  
diff --git a/src/numerics.f90 b/src/numerics.f90
index eb974b3c4..db7bf0fe4 100644
--- a/src/numerics.f90
+++ b/src/numerics.f90
@@ -42,6 +42,8 @@ module numerics
    subStepMinHomog            =  1.0e-3_pReal, &                                                    !< minimum (relative) size of sub-step allowed during cutback in homogenization
    subStepSizeCryst           =  0.25_pReal, &                                                      !< size of first substep when cutback in crystallite
    subStepSizeHomog           =  0.25_pReal, &                                                      !< size of first substep when cutback in homogenization
+   subStepSizeLp              =  0.5_pReal, &                                                       !< size of first substep when cutback in Lp calculation
+   subStepSizeLi              =  0.5_pReal, &                                                       !< size of first substep when cutback in Li calculation
    stepIncreaseCryst          =  1.5_pReal, &                                                       !< increase of next substep size when previous substep converged in crystallite
    stepIncreaseHomog          =  1.5_pReal, &                                                       !< increase of next substep size when previous substep converged in homogenization
    rTol_crystalliteState      =  1.0e-6_pReal, &                                                    !< relative tolerance in crystallite state loop 
@@ -295,6 +297,10 @@ subroutine numerics_init
          subStepSizeCryst = IO_floatValue(line,chunkPos,2_pInt)
        case ('stepincreasecryst')
          stepIncreaseCryst = IO_floatValue(line,chunkPos,2_pInt)
+       case ('substepsizelp')
+         subStepSizeLp = IO_floatValue(line,chunkPos,2_pInt)
+       case ('substepsizeli')
+         subStepSizeLi = IO_floatValue(line,chunkPos,2_pInt)
        case ('substepminhomog')
          subStepMinHomog = IO_floatValue(line,chunkPos,2_pInt)
        case ('substepsizehomog')
@@ -515,6 +521,8 @@ subroutine numerics_init
  write(6,'(a24,1x,es8.1)')  ' subStepMinCryst:        ',subStepMinCryst
  write(6,'(a24,1x,es8.1)')  ' subStepSizeCryst:       ',subStepSizeCryst
  write(6,'(a24,1x,es8.1)')  ' stepIncreaseCryst:      ',stepIncreaseCryst
+ write(6,'(a24,1x,es8.1)')  ' subStepSizeLp:          ',subStepSizeLp
+ write(6,'(a24,1x,es8.1)')  ' subStepSizeLi:          ',subStepSizeLi
  write(6,'(a24,1x,i8)')     ' nState:                 ',nState
  write(6,'(a24,1x,i8)')     ' nStress:                ',nStress
  write(6,'(a24,1x,es8.1)')  ' rTol_crystalliteState:  ',rTol_crystalliteState
@@ -643,6 +651,8 @@ subroutine numerics_init
  if (subStepMinCryst <= 0.0_pReal)         call IO_error(301_pInt,ext_msg='subStepMinCryst')
  if (subStepSizeCryst <= 0.0_pReal)        call IO_error(301_pInt,ext_msg='subStepSizeCryst')
  if (stepIncreaseCryst <= 0.0_pReal)       call IO_error(301_pInt,ext_msg='stepIncreaseCryst')
+ if (subStepSizeLp <= 0.0_pReal)           call IO_error(301_pInt,ext_msg='subStepSizeLp')
+ if (subStepSizeLi <= 0.0_pReal)           call IO_error(301_pInt,ext_msg='subStepSizeLi')
  if (subStepMinHomog <= 0.0_pReal)         call IO_error(301_pInt,ext_msg='subStepMinHomog')
  if (subStepSizeHomog <= 0.0_pReal)        call IO_error(301_pInt,ext_msg='subStepSizeHomog')
  if (stepIncreaseHomog <= 0.0_pReal)       call IO_error(301_pInt,ext_msg='stepIncreaseHomog')
diff --git a/src/plastic_isotropic.f90 b/src/plastic_isotropic.f90
index a82b9c820..bea9c616e 100644
--- a/src/plastic_isotropic.f90
+++ b/src/plastic_isotropic.f90
@@ -9,8 +9,7 @@
 module plastic_isotropic
  use prec, only: &
    pReal,&
-   pInt, &
-   DAMASK_NaN
+   pInt
  
  implicit none
  private
@@ -36,14 +35,14 @@ module plastic_isotropic
    integer(kind(undefined_ID)), allocatable, dimension(:) :: & 
      outputID
   real(pReal) :: &
-     fTaylor        = DAMASK_NaN, &
-     tau0           = DAMASK_NaN, &
-     gdot0          = DAMASK_NaN, &
-     n              = DAMASK_NaN, &
-     h0             = DAMASK_NaN, &
+     fTaylor, &
+     tau0, &
+     gdot0, &
+     n, &
+     h0, &
      h0_slopeLnRate = 0.0_pReal, &
-     tausat         = DAMASK_NaN, &
-     a              = DAMASK_NaN, &
+     tausat, &
+     a, &
      aTolFlowstress = 1.0_pReal, &
      aTolShear      = 1.0e-6_pReal, &
      tausat_SinhFitA= 0.0_pReal, &
diff --git a/src/plastic_nonlocal.f90 b/src/plastic_nonlocal.f90
index 0c0a48523..0a8c4c3f9 100644
--- a/src/plastic_nonlocal.f90
+++ b/src/plastic_nonlocal.f90
@@ -2377,9 +2377,9 @@ end subroutine plastic_nonlocal_deltaState
 !---------------------------------------------------------------------------------------------------
 subroutine plastic_nonlocal_dotState(Tstar_v, Fe, Fp, Temperature, &
                                      timestep,subfrac, ip,el)
-
-use prec,     only: DAMASK_NaN, &
-                    dNeq0, &
+use, intrinsic :: &
+  IEEE_arithmetic
+use prec,     only: dNeq0, &
                     dNeq, &
                     dEq0
 use numerics, only: numerics_integrationMode, &
@@ -2701,7 +2701,7 @@ if (.not. phase_localPlasticity(material_phase(1_pInt,ip,el))) then
       write(6,'(a)') '<< CONST >> enforcing cutback !!!'
     endif
 #endif
-    plasticState(p)%dotState =  DAMASK_NaN                                                     ! -> return NaN and, hence, enforce cutback
+    plasticState(p)%dotState = IEEE_value(1.0_pReal,IEEE_quiet_NaN)                                 ! -> return NaN and, hence, enforce cutback
     return
   endif
 
@@ -2984,7 +2984,7 @@ if (    any(rhoSglOriginal(1:ns,1:4) + rhoDot(1:ns,1:4) * timestep < -aTolRho(in
     write(6,'(a)') '<< CONST >> enforcing cutback !!!'
   endif
 #endif
-  plasticState(p)%dotState = DAMASK_NaN
+  plasticState(p)%dotState = IEEE_value(1.0_pReal,IEEE_quiet_NaN)
   return
 else
   forall (s = 1:ns, t = 1_pInt:4_pInt) 
diff --git a/src/prec.f90 b/src/prec.f90
index c6ae656b8..671e15990 100644
--- a/src/prec.f90
+++ b/src/prec.f90
@@ -5,22 +5,12 @@
 !> @author   Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH
 !> @author   Luv Sharma, Max-Planck-Institut für Eisenforschung GmbH
 !> @brief    setting precision for real and int type
-!> @details  setting precision for real and int type and for DAMASK_NaN. Definition is made 
-!!           depending on makro "INT" defined during compilation
-!!           for details on NaN see https://software.intel.com/en-us/forums/topic/294680
 !--------------------------------------------------------------------------------------------------
 module prec
-
-#if !(defined(__GFORTRAN__) && __GNUC__ < 5)
- use, intrinsic :: &                                                                                ! unfortunately not avialable in gfortran <= 5
-  IEEE_arithmetic
-#endif
-
  implicit none
  private 
 #if (FLOAT==8)
  integer,     parameter, public :: pReal = 8                                                        !< floating point double precision (was selected_real_kind(15,300), number with 15 significant digits, up to 1e+-300)
- real(pReal), parameter, public :: DAMASK_NaN = real(Z'7FF8000000000000',pReal)                     !< quiet NaN for double precision (from http://www.hpc.unimelb.edu.au/doc/f90lrm/dfum_035.html)
 #else
  NO SUITABLE PRECISION FOR REAL SELECTED, STOPPING COMPILATION
 #endif
@@ -106,7 +96,6 @@ module prec
 
  public :: &
    prec_init, &
-   prec_isNaN, &
    dEq, &
    dEq0, &
    cEq, &
@@ -118,7 +107,7 @@ contains
 
 
 !--------------------------------------------------------------------------------------------------
-!> @brief reporting precision and checking if DAMASK_NaN is set correctly
+!> @brief reporting precision
 !--------------------------------------------------------------------------------------------------
 subroutine prec_init
  use, intrinsic :: &
@@ -133,36 +122,13 @@ subroutine prec_init
  write(6,'(a,i3)')    ' Bytes for pReal:    ',pReal
  write(6,'(a,i3)')    ' Bytes for pInt:     ',pInt
  write(6,'(a,i3)')    ' Bytes for pLongInt: ',pLongInt
- write(6,'(a,e10.3)') ' NaN:           ',     DAMASK_NaN
- write(6,'(a,l3)')    ' NaN != NaN:         ',DAMASK_NaN /= DAMASK_NaN
- write(6,'(a,l3,/)')  ' NaN check passed    ',prec_isNAN(DAMASK_NaN)
 
- if ((.not. prec_isNaN(DAMASK_NaN)) .or. (DAMASK_NaN == DAMASK_NaN)) call quit(9000)
  realloc_lhs_test = [1_pInt,2_pInt]
  if (realloc_lhs_test(2)/=2_pInt) call quit(9000)
  
-
 end subroutine prec_init
 
 
-!--------------------------------------------------------------------------------------------------
-!> @brief figures out if a floating point number is NaN
-! basically just a small wrapper, because gfortran < 5.0 does not have the IEEE module
-!--------------------------------------------------------------------------------------------------
-logical elemental pure function prec_isNaN(a)
-
- implicit none
- real(pReal), intent(in) :: a
-
-#if (defined(__GFORTRAN__) && __GNUC__ < 5)
- intrinsic :: isNaN
- prec_isNaN = isNaN(a)
-#else
- prec_isNaN = IEEE_is_NaN(a)
-#endif
-end function prec_isNaN
-
-
 !--------------------------------------------------------------------------------------------------
 !> @brief equality comparison for float with double precision
 ! replaces "==" but for certain (relative) tolerance. Counterpart to dNeq
diff --git a/src/spectral_utilities.f90 b/src/spectral_utilities.f90
index a33c7a9f5..0773a9065 100644
--- a/src/spectral_utilities.f90
+++ b/src/spectral_utilities.f90
@@ -741,8 +741,8 @@ end function utilities_curlRMS
 !> @brief calculates mask compliance tensor used to adjust F to fullfill stress BC
 !--------------------------------------------------------------------------------------------------
 function utilities_maskedCompliance(rot_BC,mask_stress,C)
- use prec, only: &
-   prec_isNaN
+ use, intrinsic :: &
+   IEEE_arithmetic
  use IO, only: &
    IO_error
  use math, only: &
@@ -794,7 +794,7 @@ function utilities_maskedCompliance(rot_BC,mask_stress,C)
    endif; enddo; endif; enddo
 
    call math_invert(size_reduced, c_reduced, s_reduced, errmatinv)                                  ! invert reduced stiffness
-   if (any(prec_isNaN(s_reduced))) errmatinv = .true.
+   if (any(IEEE_is_NaN(s_reduced))) errmatinv = .true.
    if(errmatinv) call IO_error(error_ID=400_pInt,ext_msg='utilities_maskedCompliance')
    temp99_Real = 0.0_pReal                                                                          ! fill up compliance with zeros
     k = 0_pInt