diff --git a/.github/workflows/Fortran.yml b/.github/workflows/Fortran.yml
index e9175bd53..cd361f248 100644
--- a/.github/workflows/Fortran.yml
+++ b/.github/workflows/Fortran.yml
@@ -2,7 +2,7 @@ name: Grid and Mesh Solver
 on: [push]
 
 env:
-  PETSC_VERSION: '3.17.1'
+  PETSC_VERSION: '3.18.0'
   HOMEBREW_NO_ANALYTICS: 'ON'     # Make Homebrew installation a little quicker
   HOMEBREW_NO_AUTO_UPDATE: 'ON'
   HOMEBREW_NO_BOTTLE_SOURCE_FALLBACK: 'ON'
@@ -158,7 +158,10 @@ jobs:
         if: contains( matrix.intel_v, 'classic')
         run: |
           cd petsc-${PETSC_VERSION}
-          ./configure --with-fc=mpiifort --with-cc=mpiicc --with-cxx=mpiicpc \
+          ./configure \
+          --with-fc='mpiifort -fc=ifort -diag-disable=10441' \
+          --with-cc='mpiicc -cc=icc -diag-disable=10441' \
+          --with-cxx='mpiicpc -cxx=icpc -diag-disable=10441' \
           --download-fftw --download-hdf5 --download-hdf5-fortran-bindings=1 --download-zlib
           make all
 
@@ -166,7 +169,10 @@ jobs:
         if: contains( matrix.intel_v, 'llvm')
         run: |
           cd petsc-${PETSC_VERSION}
-          ./configure --with-fc=mpiifort --with-cc="mpiicc -cc=icx" --with-cxx="mpiicpc -cxx=icpx" \
+          ./configure \
+          --with-fc='mpiifort -fc=ifx' \
+          --with-cc='mpiicc -cc=icx' \
+          --with-cxx='mpiicpc -cxx=icpx' \
           --download-fftw --download-hdf5 --download-hdf5-fortran-bindings=1 --download-zlib
           make all
 
diff --git a/.github/workflows/Python.yml b/.github/workflows/Python.yml
index 2f2a5d111..1e0fa4a76 100644
--- a/.github/workflows/Python.yml
+++ b/.github/workflows/Python.yml
@@ -24,7 +24,7 @@ jobs:
       - name: Install dependencies
         run: |
           python -m pip install --upgrade pip
-          pip install pytest pandas scipy h5py vtk matplotlib pyyaml
+          pip install pytest pandas scipy h5py vtk matplotlib pyyaml build
 
       - name: Strip git hash (Unix)
         if: runner.os != 'Windows'
@@ -39,6 +39,13 @@ jobs:
           $VERSION,$_ = $VERSION -Split '-g',2,"simplematch"
           $VERSION | Out-File VERSION
 
+      - name: Build and Install
+        run: |
+          cd python
+          python -m build
+          python -m pip install dist/*.whl
+          python -c 'import damask;print(damask.__version__)'
+
       - name: Install and run unit tests (Unix)
         if: runner.os != 'Windows'
         run: |
diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
index faaf66564..735335391 100644
--- a/.gitlab-ci.yml
+++ b/.gitlab-ci.yml
@@ -63,6 +63,13 @@ create_testroot:
 
 
 ###################################################################################################
+setuptools:
+  stage: python
+  script:
+    - sed -i 's/-[[:digit:]]*-.*//' VERSION
+    - cd python
+    - python3 -m build --wheel --no-isolation
+
 pytest:
   stage: python
   script:
diff --git a/CMakeLists.txt b/CMakeLists.txt
index 6827d31b8..7fb3f2d1b 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -11,7 +11,7 @@ endif()
 project(Prerequisites LANGUAGES)
 set(ENV{PKG_CONFIG_PATH} "$ENV{PETSC_DIR}/$ENV{PETSC_ARCH}/lib/pkgconfig:$ENV{PKG_CONFIG_PATH}")
 pkg_check_modules(PETSC_MIN REQUIRED PETSc>=3.12.0 QUIET) #CMake does not support version range
-pkg_check_modules(PETSC REQUIRED PETSc<3.18.0)
+pkg_check_modules(PETSC REQUIRED PETSc<3.19.0)
 
 pkg_get_variable(CMAKE_Fortran_COMPILER PETSc fcompiler)
 pkg_get_variable(CMAKE_C_COMPILER PETSc ccompiler)
diff --git a/PRIVATE b/PRIVATE
index 8c0a24006..36c2fe3cd 160000
--- a/PRIVATE
+++ b/PRIVATE
@@ -1 +1 @@
-Subproject commit 8c0a24006f355ce8d7602daa77a33348183a2bc7
+Subproject commit 36c2fe3cde7fc1397ddd0fef71b2f4b40b72bba8
diff --git a/VERSION b/VERSION
index 246a3966a..f28804687 100644
--- a/VERSION
+++ b/VERSION
@@ -1 +1 @@
-3.0.0-alpha6-561-geb562ac41
+3.0.0-alpha7-26-gb9f354a49
diff --git a/install/MarcMentat/2022.1/Marc_tools/comp_damask_hmp.patch b/install/MarcMentat/2022.1/Marc_tools/comp_damask_hmp.patch
index 8dda8e1ce..886ebf008 100644
--- a/install/MarcMentat/2022.1/Marc_tools/comp_damask_hmp.patch
+++ b/install/MarcMentat/2022.1/Marc_tools/comp_damask_hmp.patch
@@ -32,7 +32,7 @@
  
  
      $LOAD ${program}  $DIR/lib/main.o\
-@@ -33,9 +42,13 @@ echo "program: $program"
+@@ -33,9 +42,13 @@
      $TKLIBS \
      $MRCLIBS \
      $METISLIBS \
@@ -46,4 +46,4 @@
      /bin/rm $userobj
 +	/bin/rm $DIRJOB/*.mod
 +	/bin/rm $DIRJOB/*.smod
-+    /bin/rm $DIRJOB/*_genmod.f90
++   /bin/rm $DIRJOB/*_genmod.f90
diff --git a/install/MarcMentat/2022.1/Marc_tools/comp_damask_lmp.patch b/install/MarcMentat/2022.1/Marc_tools/comp_damask_lmp.patch
index 6ea9c098f..191cb1a53 100644
--- a/install/MarcMentat/2022.1/Marc_tools/comp_damask_lmp.patch
+++ b/install/MarcMentat/2022.1/Marc_tools/comp_damask_lmp.patch
@@ -32,7 +32,7 @@
  
  
      $LOAD ${program}  $DIR/lib/main.o\
-@@ -33,9 +42,13 @@ echo "program: $program"
+@@ -33,9 +42,13 @@
      $TKLIBS \
      $MRCLIBS \
      $METISLIBS \
@@ -46,4 +46,4 @@
      /bin/rm $userobj
 +	/bin/rm $DIRJOB/*.mod
 +	/bin/rm $DIRJOB/*.smod
-+    /bin/rm $DIRJOB/*_genmod.f90
++   /bin/rm $DIRJOB/*_genmod.f90
diff --git a/install/MarcMentat/2022.1/Marc_tools/comp_damask_mp.patch b/install/MarcMentat/2022.1/Marc_tools/comp_damask_mp.patch
index 3316227d6..7c9cf7ba7 100644
--- a/install/MarcMentat/2022.1/Marc_tools/comp_damask_mp.patch
+++ b/install/MarcMentat/2022.1/Marc_tools/comp_damask_mp.patch
@@ -32,7 +32,7 @@
  
  
      $LOAD ${program}  $DIR/lib/main.o\
-@@ -33,9 +42,13 @@ echo "program: $program"
+@@ -33,9 +42,13 @@
      $TKLIBS \
      $MRCLIBS \
      $METISLIBS \
@@ -46,4 +46,4 @@
      /bin/rm $userobj
 +	/bin/rm $DIRJOB/*.mod
 +	/bin/rm $DIRJOB/*.smod
-+    /bin/rm $DIRJOB/*_genmod.f90
++   /bin/rm $DIRJOB/*_genmod.f90
diff --git a/install/MarcMentat/2022.1/Marc_tools/include_linux64.patch b/install/MarcMentat/2022.1/Marc_tools/include_linux64.patch
index d11e84f1a..18c92edd5 100644
--- a/install/MarcMentat/2022.1/Marc_tools/include_linux64.patch
+++ b/install/MarcMentat/2022.1/Marc_tools/include_linux64.patch
@@ -1,6 +1,6 @@
 ---
 +++
-@@ -166,6 +166,11 @@ if test -n "$MSCCOSIM_HOME"; then
+@@ -166,6 +166,15 @@
   MARC_COSIM_LIB="$MSCCOSIM_HOME/Cosim$MSCCOSIM_VERSION/Dcosim$MSCCOSIM_VERSION/lib"
  fi
  
@@ -16,17 +16,7 @@
  # AEM
  if test "$MARCDLLOUTDIR" = ""; then
   DLLOUTDIR="$MARC_LIB"
-@@ -477,8 +486,8 @@ if test "$MARC_INTEGER_SIZE" = "i4" ; then
-   I8DEFINES=
-   I8CDEFINES=
- else
--  I8FFLAGS="-i8"
-+  I8FFLAGS="-i8 -integer-size 64"
-   I8DEFINES="-DI64"
-   I8CDEFINES="-U_DOUBLE -D_SINGLE"
- fi
- 
-@@ -594,7 +605,7 @@ then
+@@ -594,7 +603,7 @@
  PROFILE=" $PROFILE -pg"
  fi
  
@@ -35,7 +25,7 @@
  if test "$MTHREAD" = "OPENMP"
  then
    FORT_OPT=" $FORT_OPT -qopenmp"
-@@ -607,7 +616,7 @@ else
+@@ -607,7 +616,7 @@
    FORT_OPT=" $FORT_OPT -save -zero"
  fi
  if test "$MARCHDF_HDF" = "HDF"; then
@@ -44,10 +34,10 @@
  fi
  
  FORTLOW="$FCOMP $FORT_OPT $PROFILE -O0 $I8FFLAGS -I$MARC_SOURCE/common \
-@@ -621,6 +630,29 @@ FORTNA="$FCOMP $FORT_OPT -fno-alias -O3 $I8FFLAGS -I$MARC_SOURCE/common \
+@@ -621,6 +630,29 @@
  #  for compiling free form f90 files. high opt, integer(4)
  FORTF90="$FCOMP -c -O3"
-
+ 
 +# determine DAMASK version
 +if test -n "$DAMASK_USER"; then
 +  DAMASKROOT=`dirname $DAMASK_USER`/..
@@ -74,30 +64,12 @@
  if test "$MARCDEBUG" = "ON"
  then
    FORTLOW="$FCOMP $FORT_OPT $PROFILE $I8FFLAGS -I$MARC_SOURCE/common \
-@@ -778,7 +806,7 @@ SECLIBS="-L$MARC_LIB -llapi"
-
+@@ -778,7 +810,7 @@
+ 
  SOLVERLIBS="${BCSSOLVERLIBS} ${VKISOLVERLIBS} ${CASISOLVERLIBS} ${MF2SOLVERLIBS} \
    -L$MARC_MKL \
 -  $MARC_LIB/blas_src.a ${ACSI_LIB}/ACSI_MarcLib.a $KDTREE2_LIB/libkdtree2.a $MARC_LIB/libtetmeshinterface.a   $MARC_LIB/libcaefatigueinterface.a -L$MARC_LIB -lmkl_blacs_intelmpi_ilp64 -lmkl_scalapack_ilp64 -lmkl_intel_ilp64 -lmkl_intel_thread -lmkl_core -liomp5 -ltetmesh -lmeshgems -lmg-tetra -lmeshgems_stubs $HDF_LIBS $SOLVER2LIBS"
 +  $MARC_LIB/blas_src.a ${ACSI_LIB}/ACSI_MarcLib.a $KDTREE2_LIB/libkdtree2.a $MARC_LIB/libtetmeshinterface.a   $MARC_LIB/libcaefatigueinterface.a -L$MARC_LIB -lmkl_blacs_intelmpi_ilp64 -lmkl_scalapack_ilp64 -lmkl_intel_ilp64 -lmkl_intel_thread -lmkl_core -liomp5 -ltetmesh -lmeshgems -lmg-tetra -lmeshgems_stubs $HDF5_LIB $SOLVER2LIBS"
-
+ 
  SOLVERLIBS_DLL=${SOLVERLIBS}
  if test "$AEM_DLL" -eq 1
-@@ -802,7 +830,7 @@ then
-  OPENSSL=NONE
- fi
- 
--SYSLIBS=" $OPENMP -lpthread -shared-intel -cxxlib $MARC_RPC_LIB"
-+SYSLIBS=" $OPENMP -lpthread -cxxlib $MARC_RPC_LIB"
-
- # Uncomment the following lines to turn on the trace and comment out the next 4 lines
- # if test $MPITYPE = intelmpi
-@@ -812,7 +840,7 @@ SYSLIBS=" $OPENMP -lpthread -shared-intel -cxxlib $MARC_RPC_LIB"
- # fi
- if test $MPITYPE = intelmpi
- then
-- SYSLIBS="-L${MPI_ROOT}/lib/release -lmpi -L${MPI_ROOT}/lib -lmpifort -lrt $OPENMP -threads -lpthread -shared-intel -cxxlib $MARC_RPC_LIB"
-+ SYSLIBS="-L${MPI_ROOT}/lib/release -lmpi -L${MPI_ROOT}/lib -lmpifort -lrt $OPENMP -threads -lpthread -cxxlib $MARC_RPC_LIB"
- fi
- 
- if test "$ZLIB" = "ZLIB"; then
diff --git a/install/MarcMentat/2022.1/Marc_tools/run_damask_hmp.patch b/install/MarcMentat/2022.1/Marc_tools/run_damask_hmp.patch
index e3a364e00..cdeeecf4b 100644
--- a/install/MarcMentat/2022.1/Marc_tools/run_damask_hmp.patch
+++ b/install/MarcMentat/2022.1/Marc_tools/run_damask_hmp.patch
@@ -1,6 +1,6 @@
 ---
 +++
-@@ -302,7 +302,23 @@ fi
+@@ -302,7 +302,23 @@
  
  . "$DIR/getarch"
  
@@ -24,7 +24,7 @@
  #
  
  #
-@@ -405,7 +421,7 @@ sid=
+@@ -405,7 +421,7 @@
  did=
  vid=
  user=
@@ -33,16 +33,7 @@
  objs=
  qid=background
  cpu=
-@@ -573,7 +589,7 @@ do
- 		    justlist=yes
- 		;;
-                 -fe* | -FE*)
--                    feature=$value
-+		    feature=$value
- 
- 		;;
- 		-pr* | -PR*)
-@@ -676,50 +692,19 @@ do
+@@ -676,50 +692,19 @@
                          esac        
                  ;;
  		-u* | -U*)
@@ -99,39 +90,7 @@
  		;;
  		-obj | -OBJ)
  			objs="$value"
-@@ -739,19 +724,19 @@ do
- 				;;
- 			esac
- 		;;
--                -dl | -DL)
--                        case $value in
--                                y* | Y*)
--                                deletelog=yes
--                                ;;
--                                n* | N*)
--                                deletelog=no
--                                ;;
--                                *)
--                                ;;
--                        esac
-+        -dl | -DL)
-+            case $value in
-+                    y* | Y*)
-+                    deletelog=yes
-+                    ;;
-+                    n* | N*)
-+                    deletelog=no
-+                    ;;
-+                    *)
-+                    ;;
-+            esac
- 
--                ;;
-+        ;;
- 		-at | -AT)
- 			att=$value
- 		;;
-@@ -1207,12 +1192,12 @@ post file $DIRPID/$pid.t16 or $DIRPID/$pid.t19 or $DIRPID/$pid.h5 not accessible
+@@ -1207,12 +1192,12 @@
  				fi
  			fi
  		fi
@@ -147,16 +106,7 @@
  			fi
  		fi
  		if test "$objs"
-@@ -1386,7 +1371,7 @@ else
-     else
- 	error="$error
- job id required"
--    fi
-+fi
- fi
- 
- case $qid in
-@@ -1531,7 +1516,7 @@ Program name         : $prog
+@@ -1531,7 +1516,7 @@
  Marc shared lib      : $progdll
  Version type         : $mode
  Job ID               : $DIRJID/$jid$extra_job_info
@@ -165,7 +115,7 @@
  User objects/libs    : $objs
  Restart file job ID  : $rid
  Substructure file ID : $sid
-@@ -1564,7 +1549,7 @@ Program name         : $prog
+@@ -1564,7 +1549,7 @@
  Marc shared lib      : $progdll
  Version type         : $mode
  Job ID               : $DIRJID/$jid$extra_job_info
@@ -174,7 +124,7 @@
  User objects/libs    : $objs
  Restart file job ID  : $rid
  Substructure file ID : $sid
-@@ -1687,7 +1672,7 @@ Program name ($prog)? $ECHOTXT"
+@@ -1687,7 +1672,7 @@
                      ;;
                    esac        
                  fi
@@ -183,7 +133,7 @@
  		read value
  		if test "$value"
  		then
-@@ -1696,50 +1681,19 @@ Program name ($prog)? $ECHOTXT"
+@@ -1696,50 +1681,19 @@
                        user=
                      ;;
                      *)
@@ -209,16 +159,13 @@
 -                      fi
 +                      user=$value
                        case $user in
--                         \/*)
--                           ;;
--                         *)
-+                        \/*)
-+                        ;;
-+                        *)
+                          \/*)
+                            ;;
+                          *)
                            user=`pwd`/$user
 -                          usersubname=`pwd`/$usersubname
--                           ;;
--                      esac
+                            ;;
+                       esac
 -                      if test ! -f $usersubname
 -                      then
 -                       if test -f $usersubname.f 
@@ -235,8 +182,6 @@
 -                        usersubname=$usersubname.F90
 -                       fi
 -                      fi
-+                        ;;
-+                      esac        
 +					  usernoext=$user
 +					  usernoext=`dirname $usernoext`/`$BASENAME $usernoext .f`
 +					  usernoext=`dirname $usernoext`/`$BASENAME $usernoext .F`
@@ -245,7 +190,7 @@
                      ;;
                    esac
  		fi
-@@ -2274,11 +2228,12 @@ fi
+@@ -2274,11 +2228,12 @@
  #
  # user subroutine used
  #
@@ -260,7 +205,7 @@
      case $program in
        \/* | \.\/*)
          bd=
-@@ -2391,7 +2346,7 @@ fi
+@@ -2391,7 +2346,7 @@
      fi
      if test "$user"
      then
@@ -269,7 +214,7 @@
       usersub=1
      fi
      export execpath      
-@@ -3274,44 +3229,27 @@ then
+@@ -3274,44 +3229,27 @@
      echo
      if test "$user"
      then
@@ -319,21 +264,15 @@
      fi
  
  
-@@ -3331,10 +3269,11 @@ then
+@@ -3331,6 +3269,7 @@
      $TKLIBS  \
      $MRCLIBS     \
      $METISLIBS \
--    $SFLIB \
 +    $DAMASK   \
-+    $SFLIB  \
+     $SFLIB \
      $OPENSSL_LIB \
      $SYSLIBS  \
--    $SECLIBS  || \
-+    $SECLIBS      || \
-     {
-     echo "$0: link failed for ${user:+$userobj }$objs"
-     exit 1
-@@ -3344,6 +3283,9 @@ else
+@@ -3344,6 +3283,9 @@
  	prgsav=yes
  fi
  /bin/rm $userobj 2>/dev/null
@@ -343,7 +282,7 @@
  
  #
  # run marc
-@@ -3390,7 +3332,7 @@ if test $dllrun -eq 0; then
+@@ -3390,7 +3332,7 @@
    fi
  else
    if test $cpdll = yes; then
@@ -352,7 +291,7 @@
       /bin/cp $DIRJOB/$marcdll $DIRJOB/${filename}_$marcdll 2>/dev/null
    fi
    if test $rmdll = yes
-@@ -3556,7 +3498,7 @@ then
+@@ -3556,7 +3498,7 @@
                # first copy over the user sub if local directories
                if test ${dirstatus[$counter]} = "local"
                then
@@ -361,7 +300,7 @@
                fi
                # do the compilation on the other machine
                if test ${dirstatus[$counter]} = "shared"
-@@ -3569,21 +3511,21 @@ then
+@@ -3569,21 +3511,21 @@
                remoteuser=$DIR1/`$BASENAME $user`
                $RSH $i /bin/rm $remoteprog 2> /dev/null
                echo
@@ -386,7 +325,7 @@
                fi
              fi
            fi
-@@ -3593,39 +3535,27 @@ then
+@@ -3593,39 +3535,27 @@
      if test "$userhost"
      then
        echo
@@ -432,21 +371,15 @@
    fi # if test $user
  
  
-@@ -3645,10 +3575,11 @@ then
+@@ -3645,6 +3575,7 @@
    $TKLIBS  \
    $MRCLIBS     \
    $METISLIBS \
--  $SFLIB \
 +  $DAMASK \
-+  $SFLIB  \
+   $SFLIB \
    $OPENSSL_LIB \
    $SYSLIBS  \
--  $SECLIBS  || \
-+  $SECLIBS    || \
-   {
-     echo "$0: link failed for ${user:+$userobj }$objs"
-     echo "         $PRODUCT Exit number 3"
-@@ -3686,6 +3617,9 @@ else # if test $link
+@@ -3686,6 +3617,9 @@
  	prgsav=yes
  fi   # if test $link
  /bin/rm $userobj 2>/dev/null
@@ -456,69 +389,7 @@
  
  #
  # run marc
-@@ -3744,42 +3678,42 @@ then
-       counter=0
-       if test -f "$host_filt"
-       then
--        for i in `$AWK -v n=$numfield '{print $n}' $host_filt`
--        do
--          ibase=${i%%.*}
--          if test $ibase != $thishost
-+      for i in `$AWK -v n=$numfield '{print $n}' $host_filt`
-+      do
-+        ibase=${i%%.*}
-+        if test $ibase != $thishost
-+        then
-+          counter=$((counter+1))
-+          DIR1=$DIRJOB
-+          line=`grep -v '^#' $userhost | grep "^$ibase "`
-+          workdir=`echo $line | $AWK '{print $3}'`
-+          if test -n "$workdir"
-           then
--            counter=$((counter+1))
--            DIR1=$DIRJOB
--            line=`grep -v '^#' $userhost | grep "^$ibase "`
--            workdir=`echo $line | $AWK '{print $3}'`
--            if test -n "$workdir"
--            then
--              DIR1=$workdir
--            fi
--            #  if an incompatible host uses shared directory,
--            #  then the root machine deletes the executable
--            if test ${dirstatus[$counter]} = "shared" -a  ${compstatus[$counter]} = "no"
--            then
--              hname=_$ibase
--              /bin/rm ${execname}$hname
--            fi
--            # if local directory used, the remote machine
--            # deletes the executable
--            if test ${dirstatus[$counter]} = "local"
--            then
--              $RSH $i /bin/rm $DIR1/${execname} 2>/dev/null
--            fi
-+            DIR1=$workdir
-           fi
--        done
--      fi
-+          #  if an incompatible host uses shared directory,
-+          #  then the root machine deletes the executable
-+          if test ${dirstatus[$counter]} = "shared" -a  ${compstatus[$counter]} = "no"
-+          then
-+            hname=_$ibase
-+            /bin/rm ${execname}$hname
-+          fi
-+          # if local directory used, the remote machine
-+          # deletes the executable
-+          if test ${dirstatus[$counter]} = "local"
-+          then
-+            $RSH $i /bin/rm $DIR1/${execname} 2>/dev/null
-+          fi
-+        fi
-+      done
-     fi
-   fi
- fi
-+fi
+@@ -3779,7 +3713,7 @@
  else
  #dllrun >0
    if test $cpdll = yes; then
@@ -527,7 +398,7 @@
       /bin/cp $DIRJOB/$marcdll $DIRJOB/${filename}_$marcdll 2>/dev/null
    fi
    if test $rmdll = yes;then
-@@ -3904,7 +3838,7 @@ then
+@@ -3904,7 +3838,7 @@
                # first copy over the user sub if local directories
                if test ${dirstatus[$counter]} = "local"
                then
@@ -536,7 +407,7 @@
                fi
                # do the compilation on the other machine
                if test ${dirstatus[$counter]} = "shared"
-@@ -3917,20 +3851,20 @@ then
+@@ -3917,20 +3851,20 @@
                remoteuser=$DIR1/`$BASENAME $user`
                $RSH $i /bin/rm $remoteprog 2> /dev/null
                echo
@@ -560,7 +431,7 @@
                fi
              fi
            fi
-@@ -3940,37 +3874,25 @@ then
+@@ -3940,37 +3874,25 @@
      if test "$userhost"
      then
        echo
@@ -604,7 +475,7 @@
    fi # if test $user
  
  
-@@ -3990,10 +3912,11 @@ then
+@@ -3990,6 +3912,7 @@
    $TKLIBS  \
    $MRCLIBS     \
    $METISLIBS \
@@ -612,12 +483,7 @@
    $SFLIB \
    $OPENSSL_LIB \
    $SYSLIBS \
--  $SECLIBS  || \
-+  $SECLIBS         || \
-   {
-     echo "$0: link failed for ${user:+$userobj }$objs"
-     exit 1
-@@ -4030,7 +3953,9 @@ else # if test $link
+@@ -4030,7 +3953,9 @@
  	prgsav=yes
  fi   # if test $link
  /bin/rm $userobj 2>/dev/null
@@ -628,78 +494,7 @@
  # done if no job id given
  if test -z "$jid"
  then
-@@ -4070,7 +3995,7 @@ if test $ddm_arc -gt 0; then
-   RUN_JOB="$BINDIR/exeddm $RUN_JOB -ddm $ddm_arc "
- fi
- 
--$RUN_JOB 
-+    $RUN_JOB 
- 
- if test $nprocd -gt 1
- then
-@@ -4114,42 +4039,42 @@ then
-       counter=0
-       if test -f "$host_filt"
-       then
--        for i in `$AWK -v n=$numfield '{print $n}' $host_filt`
--        do
--          ibase=${i%%.*}
--          if test $ibase != $thishost
-+      for i in `$AWK -v n=$numfield '{print $n}' $host_filt`
-+      do
-+        ibase=${i%%.*}
-+        if test $ibase != $thishost
-+        then
-+          counter=$((counter+1))
-+          DIR1=$DIRJOB
-+          line=`grep -v '^#' $userhost | grep "^$ibase "`
-+          workdir=`echo $line | $AWK '{print $3}'`
-+          if test -n "$workdir"
-           then
--            counter=$((counter+1))
--            DIR1=$DIRJOB
--            line=`grep -v '^#' $userhost | grep "^$ibase "`
--            workdir=`echo $line | $AWK '{print $3}'`
--            if test -n "$workdir"
--            then
--              DIR1=$workdir
--            fi
--            #  if an incompatible host uses shared directory,
--            #  then the root machine deletes the executable
--            if test ${dirstatus[$counter]} = "shared" -a  ${compstatus[$counter]} = "no"
--            then
--              hname=_$ibase
--              /bin/rm ${execname}$hname
--            fi
--            # if local directory used, the remote machine
--            # deletes the executable
--            if test ${dirstatus[$counter]} = "local"
--            then
--              $RSH $i /bin/rm $DIR1/${execname} 2>/dev/null
--            fi
-+            DIR1=$workdir
-           fi
--        done
--      fi
-+          #  if an incompatible host uses shared directory,
-+          #  then the root machine deletes the executable
-+          if test ${dirstatus[$counter]} = "shared" -a  ${compstatus[$counter]} = "no"
-+          then
-+            hname=_$ibase
-+            /bin/rm ${execname}$hname
-+          fi
-+          # if local directory used, the remote machine
-+          # deletes the executable
-+          if test ${dirstatus[$counter]} = "local"
-+          then
-+            $RSH $i /bin/rm $DIR1/${execname} 2>/dev/null
-+          fi
-+        fi
-+      done
-     fi
-   fi
- fi
-+fi
+@@ -4149,7 +4074,7 @@
  else
  #dllrun >0
    if test $cpdll = yes; then
diff --git a/install/MarcMentat/2022.1/Marc_tools/run_damask_lmp.patch b/install/MarcMentat/2022.1/Marc_tools/run_damask_lmp.patch
index a5989bfca..21ed967f6 100644
--- a/install/MarcMentat/2022.1/Marc_tools/run_damask_lmp.patch
+++ b/install/MarcMentat/2022.1/Marc_tools/run_damask_lmp.patch
@@ -1,6 +1,6 @@
 ---
 +++
-@@ -302,7 +302,23 @@ fi
+@@ -302,7 +302,23 @@
  
  . "$DIR/getarch"
  
@@ -24,7 +24,7 @@
  #
  
  #
-@@ -405,7 +421,7 @@ sid=
+@@ -405,7 +421,7 @@
  did=
  vid=
  user=
@@ -33,16 +33,7 @@
  objs=
  qid=background
  cpu=
-@@ -573,7 +589,7 @@ do
- 		    justlist=yes
- 		;;
-                 -fe* | -FE*)
--                    feature=$value
-+		    feature=$value
- 
- 		;;
- 		-pr* | -PR*)
-@@ -676,50 +692,19 @@ do
+@@ -676,50 +692,19 @@
                          esac        
                  ;;
  		-u* | -U*)
@@ -99,39 +90,7 @@
  		;;
  		-obj | -OBJ)
  			objs="$value"
-@@ -739,19 +724,19 @@ do
- 				;;
- 			esac
- 		;;
--                -dl | -DL)
--                        case $value in
--                                y* | Y*)
--                                deletelog=yes
--                                ;;
--                                n* | N*)
--                                deletelog=no
--                                ;;
--                                *)
--                                ;;
--                        esac
-+        -dl | -DL)
-+            case $value in
-+                    y* | Y*)
-+                    deletelog=yes
-+                    ;;
-+                    n* | N*)
-+                    deletelog=no
-+                    ;;
-+                    *)
-+                    ;;
-+            esac
- 
--                ;;
-+        ;;
- 		-at | -AT)
- 			att=$value
- 		;;
-@@ -1207,12 +1192,12 @@ post file $DIRPID/$pid.t16 or $DIRPID/$pid.t19 or $DIRPID/$pid.h5 not accessible
+@@ -1207,12 +1192,12 @@
  				fi
  			fi
  		fi
@@ -147,16 +106,7 @@
  			fi
  		fi
  		if test "$objs"
-@@ -1386,7 +1371,7 @@ else
-     else
- 	error="$error
- job id required"
--    fi
-+fi
- fi
- 
- case $qid in
-@@ -1531,7 +1516,7 @@ Program name         : $prog
+@@ -1531,7 +1516,7 @@
  Marc shared lib      : $progdll
  Version type         : $mode
  Job ID               : $DIRJID/$jid$extra_job_info
@@ -165,7 +115,7 @@
  User objects/libs    : $objs
  Restart file job ID  : $rid
  Substructure file ID : $sid
-@@ -1564,7 +1549,7 @@ Program name         : $prog
+@@ -1564,7 +1549,7 @@
  Marc shared lib      : $progdll
  Version type         : $mode
  Job ID               : $DIRJID/$jid$extra_job_info
@@ -174,7 +124,7 @@
  User objects/libs    : $objs
  Restart file job ID  : $rid
  Substructure file ID : $sid
-@@ -1687,7 +1672,7 @@ Program name ($prog)? $ECHOTXT"
+@@ -1687,7 +1672,7 @@
                      ;;
                    esac        
                  fi
@@ -183,7 +133,7 @@
  		read value
  		if test "$value"
  		then
-@@ -1696,50 +1681,19 @@ Program name ($prog)? $ECHOTXT"
+@@ -1696,50 +1681,19 @@
                        user=
                      ;;
                      *)
@@ -209,16 +159,13 @@
 -                      fi
 +                      user=$value
                        case $user in
--                         \/*)
--                           ;;
--                         *)
-+                        \/*)
-+                        ;;
-+                        *)
+                          \/*)
+                            ;;
+                          *)
                            user=`pwd`/$user
 -                          usersubname=`pwd`/$usersubname
--                           ;;
--                      esac
+                            ;;
+                       esac
 -                      if test ! -f $usersubname
 -                      then
 -                       if test -f $usersubname.f 
@@ -235,8 +182,6 @@
 -                        usersubname=$usersubname.F90
 -                       fi
 -                      fi
-+                        ;;
-+                      esac        
 +					  usernoext=$user
 +					  usernoext=`dirname $usernoext`/`$BASENAME $usernoext .f`
 +					  usernoext=`dirname $usernoext`/`$BASENAME $usernoext .F`
@@ -245,7 +190,7 @@
                      ;;
                    esac
  		fi
-@@ -2274,11 +2228,12 @@ fi
+@@ -2274,11 +2228,12 @@
  #
  # user subroutine used
  #
@@ -260,7 +205,7 @@
      case $program in
        \/* | \.\/*)
          bd=
-@@ -2391,7 +2346,7 @@ fi
+@@ -2391,7 +2346,7 @@
      fi
      if test "$user"
      then
@@ -269,7 +214,7 @@
       usersub=1
      fi
      export execpath      
-@@ -3274,44 +3229,27 @@ then
+@@ -3274,44 +3229,27 @@
      echo
      if test "$user"
      then
@@ -319,21 +264,15 @@
      fi
  
  
-@@ -3331,10 +3269,11 @@ then
+@@ -3331,6 +3269,7 @@
      $TKLIBS  \
      $MRCLIBS     \
      $METISLIBS \
--    $SFLIB \
 +    $DAMASK   \
-+    $SFLIB  \
+     $SFLIB \
      $OPENSSL_LIB \
      $SYSLIBS  \
--    $SECLIBS  || \
-+    $SECLIBS      || \
-     {
-     echo "$0: link failed for ${user:+$userobj }$objs"
-     exit 1
-@@ -3344,6 +3283,9 @@ else
+@@ -3344,6 +3283,9 @@
  	prgsav=yes
  fi
  /bin/rm $userobj 2>/dev/null
@@ -343,7 +282,7 @@
  
  #
  # run marc
-@@ -3390,7 +3332,7 @@ if test $dllrun -eq 0; then
+@@ -3390,7 +3332,7 @@
    fi
  else
    if test $cpdll = yes; then
@@ -352,7 +291,7 @@
       /bin/cp $DIRJOB/$marcdll $DIRJOB/${filename}_$marcdll 2>/dev/null
    fi
    if test $rmdll = yes
-@@ -3556,7 +3498,7 @@ then
+@@ -3556,7 +3498,7 @@
                # first copy over the user sub if local directories
                if test ${dirstatus[$counter]} = "local"
                then
@@ -361,7 +300,7 @@
                fi
                # do the compilation on the other machine
                if test ${dirstatus[$counter]} = "shared"
-@@ -3569,21 +3511,21 @@ then
+@@ -3569,21 +3511,21 @@
                remoteuser=$DIR1/`$BASENAME $user`
                $RSH $i /bin/rm $remoteprog 2> /dev/null
                echo
@@ -386,7 +325,7 @@
                fi
              fi
            fi
-@@ -3593,39 +3535,27 @@ then
+@@ -3593,39 +3535,27 @@
      if test "$userhost"
      then
        echo
@@ -432,21 +371,15 @@
    fi # if test $user
  
  
-@@ -3645,10 +3575,11 @@ then
+@@ -3645,6 +3575,7 @@
    $TKLIBS  \
    $MRCLIBS     \
    $METISLIBS \
--  $SFLIB \
 +  $DAMASK \
-+  $SFLIB  \
+   $SFLIB \
    $OPENSSL_LIB \
    $SYSLIBS  \
--  $SECLIBS  || \
-+  $SECLIBS    || \
-   {
-     echo "$0: link failed for ${user:+$userobj }$objs"
-     echo "         $PRODUCT Exit number 3"
-@@ -3686,6 +3617,9 @@ else # if test $link
+@@ -3686,6 +3617,9 @@
  	prgsav=yes
  fi   # if test $link
  /bin/rm $userobj 2>/dev/null
@@ -456,69 +389,7 @@
  
  #
  # run marc
-@@ -3744,42 +3678,42 @@ then
-       counter=0
-       if test -f "$host_filt"
-       then
--        for i in `$AWK -v n=$numfield '{print $n}' $host_filt`
--        do
--          ibase=${i%%.*}
--          if test $ibase != $thishost
-+      for i in `$AWK -v n=$numfield '{print $n}' $host_filt`
-+      do
-+        ibase=${i%%.*}
-+        if test $ibase != $thishost
-+        then
-+          counter=$((counter+1))
-+          DIR1=$DIRJOB
-+          line=`grep -v '^#' $userhost | grep "^$ibase "`
-+          workdir=`echo $line | $AWK '{print $3}'`
-+          if test -n "$workdir"
-           then
--            counter=$((counter+1))
--            DIR1=$DIRJOB
--            line=`grep -v '^#' $userhost | grep "^$ibase "`
--            workdir=`echo $line | $AWK '{print $3}'`
--            if test -n "$workdir"
--            then
--              DIR1=$workdir
--            fi
--            #  if an incompatible host uses shared directory,
--            #  then the root machine deletes the executable
--            if test ${dirstatus[$counter]} = "shared" -a  ${compstatus[$counter]} = "no"
--            then
--              hname=_$ibase
--              /bin/rm ${execname}$hname
--            fi
--            # if local directory used, the remote machine
--            # deletes the executable
--            if test ${dirstatus[$counter]} = "local"
--            then
--              $RSH $i /bin/rm $DIR1/${execname} 2>/dev/null
--            fi
-+            DIR1=$workdir
-           fi
--        done
--      fi
-+          #  if an incompatible host uses shared directory,
-+          #  then the root machine deletes the executable
-+          if test ${dirstatus[$counter]} = "shared" -a  ${compstatus[$counter]} = "no"
-+          then
-+            hname=_$ibase
-+            /bin/rm ${execname}$hname
-+          fi
-+          # if local directory used, the remote machine
-+          # deletes the executable
-+          if test ${dirstatus[$counter]} = "local"
-+          then
-+            $RSH $i /bin/rm $DIR1/${execname} 2>/dev/null
-+          fi
-+        fi
-+      done
-     fi
-   fi
- fi
-+fi
+@@ -3779,7 +3713,7 @@
  else
  #dllrun >0
    if test $cpdll = yes; then
@@ -527,7 +398,7 @@
       /bin/cp $DIRJOB/$marcdll $DIRJOB/${filename}_$marcdll 2>/dev/null
    fi
    if test $rmdll = yes;then
-@@ -3904,7 +3838,7 @@ then
+@@ -3904,7 +3838,7 @@
                # first copy over the user sub if local directories
                if test ${dirstatus[$counter]} = "local"
                then
@@ -536,7 +407,7 @@
                fi
                # do the compilation on the other machine
                if test ${dirstatus[$counter]} = "shared"
-@@ -3917,20 +3851,20 @@ then
+@@ -3917,20 +3851,20 @@
                remoteuser=$DIR1/`$BASENAME $user`
                $RSH $i /bin/rm $remoteprog 2> /dev/null
                echo
@@ -560,7 +431,7 @@
                fi
              fi
            fi
-@@ -3940,37 +3874,25 @@ then
+@@ -3940,37 +3874,25 @@
      if test "$userhost"
      then
        echo
@@ -604,7 +475,7 @@
    fi # if test $user
  
  
-@@ -3990,10 +3912,11 @@ then
+@@ -3990,6 +3912,7 @@
    $TKLIBS  \
    $MRCLIBS     \
    $METISLIBS \
@@ -612,12 +483,7 @@
    $SFLIB \
    $OPENSSL_LIB \
    $SYSLIBS \
--  $SECLIBS  || \
-+  $SECLIBS         || \
-   {
-     echo "$0: link failed for ${user:+$userobj }$objs"
-     exit 1
-@@ -4030,7 +3953,9 @@ else # if test $link
+@@ -4030,7 +3953,9 @@
  	prgsav=yes
  fi   # if test $link
  /bin/rm $userobj 2>/dev/null
@@ -628,78 +494,7 @@
  # done if no job id given
  if test -z "$jid"
  then
-@@ -4070,7 +3995,7 @@ if test $ddm_arc -gt 0; then
-   RUN_JOB="$BINDIR/exeddm $RUN_JOB -ddm $ddm_arc "
- fi
- 
--$RUN_JOB 
-+    $RUN_JOB 
- 
- if test $nprocd -gt 1
- then
-@@ -4114,42 +4039,42 @@ then
-       counter=0
-       if test -f "$host_filt"
-       then
--        for i in `$AWK -v n=$numfield '{print $n}' $host_filt`
--        do
--          ibase=${i%%.*}
--          if test $ibase != $thishost
-+      for i in `$AWK -v n=$numfield '{print $n}' $host_filt`
-+      do
-+        ibase=${i%%.*}
-+        if test $ibase != $thishost
-+        then
-+          counter=$((counter+1))
-+          DIR1=$DIRJOB
-+          line=`grep -v '^#' $userhost | grep "^$ibase "`
-+          workdir=`echo $line | $AWK '{print $3}'`
-+          if test -n "$workdir"
-           then
--            counter=$((counter+1))
--            DIR1=$DIRJOB
--            line=`grep -v '^#' $userhost | grep "^$ibase "`
--            workdir=`echo $line | $AWK '{print $3}'`
--            if test -n "$workdir"
--            then
--              DIR1=$workdir
--            fi
--            #  if an incompatible host uses shared directory,
--            #  then the root machine deletes the executable
--            if test ${dirstatus[$counter]} = "shared" -a  ${compstatus[$counter]} = "no"
--            then
--              hname=_$ibase
--              /bin/rm ${execname}$hname
--            fi
--            # if local directory used, the remote machine
--            # deletes the executable
--            if test ${dirstatus[$counter]} = "local"
--            then
--              $RSH $i /bin/rm $DIR1/${execname} 2>/dev/null
--            fi
-+            DIR1=$workdir
-           fi
--        done
--      fi
-+          #  if an incompatible host uses shared directory,
-+          #  then the root machine deletes the executable
-+          if test ${dirstatus[$counter]} = "shared" -a  ${compstatus[$counter]} = "no"
-+          then
-+            hname=_$ibase
-+            /bin/rm ${execname}$hname
-+          fi
-+          # if local directory used, the remote machine
-+          # deletes the executable
-+          if test ${dirstatus[$counter]} = "local"
-+          then
-+            $RSH $i /bin/rm $DIR1/${execname} 2>/dev/null
-+          fi
-+        fi
-+      done
-     fi
-   fi
- fi
-+fi
+@@ -4149,7 +4074,7 @@
  else
  #dllrun >0
    if test $cpdll = yes; then
diff --git a/install/MarcMentat/2022.1/Marc_tools/run_damask_mp.patch b/install/MarcMentat/2022.1/Marc_tools/run_damask_mp.patch
index 924af6fbd..cdeeecf4b 100644
--- a/install/MarcMentat/2022.1/Marc_tools/run_damask_mp.patch
+++ b/install/MarcMentat/2022.1/Marc_tools/run_damask_mp.patch
@@ -1,8 +1,6 @@
-diff --git "a/C:\\Users\\f.roters\\Documents\\Forschung\\FEM\\Kristallplastizit\303\244t\\Programme\\Userroutinen\\Marc\\Marc_mods\\2022.1\\Marc_tools\\run_marc.original" "b/C:\\Users\\f.roters\\Documents\\Forschung\\FEM\\Kristallplastizit\303\244t\\Programme\\Userroutinen\\Marc\\Marc_mods\\2022.1\\Marc_tools\\run_damask_mp"
-index 88a3027..85beabe 100644
---- "a/C:\\Users\\f.roters\\Documents\\Forschung\\FEM\\Kristallplastizit\303\244t\\Programme\\Userroutinen\\Marc\\Marc_mods\\2022.1\\Marc_tools\\run_marc.original"
-+++ "b/C:\\Users\\f.roters\\Documents\\Forschung\\FEM\\Kristallplastizit\303\244t\\Programme\\Userroutinen\\Marc\\Marc_mods\\2022.1\\Marc_tools\\run_damask_mp"
-@@ -302,7 +302,23 @@ fi
+---
++++
+@@ -302,7 +302,23 @@
  
  . "$DIR/getarch"
  
@@ -26,7 +24,7 @@ index 88a3027..85beabe 100644
  #
  
  #
-@@ -405,7 +421,7 @@ sid=
+@@ -405,7 +421,7 @@
  did=
  vid=
  user=
@@ -35,16 +33,7 @@ index 88a3027..85beabe 100644
  objs=
  qid=background
  cpu=
-@@ -573,7 +589,7 @@ do
- 		    justlist=yes
- 		;;
-                 -fe* | -FE*)
--                    feature=$value
-+		    feature=$value
- 
- 		;;
- 		-pr* | -PR*)
-@@ -676,50 +692,19 @@ do
+@@ -676,50 +692,19 @@
                          esac        
                  ;;
  		-u* | -U*)
@@ -101,39 +90,7 @@ index 88a3027..85beabe 100644
  		;;
  		-obj | -OBJ)
  			objs="$value"
-@@ -739,19 +724,19 @@ do
- 				;;
- 			esac
- 		;;
--                -dl | -DL)
--                        case $value in
--                                y* | Y*)
--                                deletelog=yes
--                                ;;
--                                n* | N*)
--                                deletelog=no
--                                ;;
--                                *)
--                                ;;
--                        esac
-+        -dl | -DL)
-+            case $value in
-+                    y* | Y*)
-+                    deletelog=yes
-+                    ;;
-+                    n* | N*)
-+                    deletelog=no
-+                    ;;
-+                    *)
-+                    ;;
-+            esac
- 
--                ;;
-+        ;;
- 		-at | -AT)
- 			att=$value
- 		;;
-@@ -1207,12 +1192,12 @@ post file $DIRPID/$pid.t16 or $DIRPID/$pid.t19 or $DIRPID/$pid.h5 not accessible
+@@ -1207,12 +1192,12 @@
  				fi
  			fi
  		fi
@@ -149,16 +106,7 @@ index 88a3027..85beabe 100644
  			fi
  		fi
  		if test "$objs"
-@@ -1386,7 +1371,7 @@ else
-     else
- 	error="$error
- job id required"
--    fi
-+fi
- fi
- 
- case $qid in
-@@ -1531,7 +1516,7 @@ Program name         : $prog
+@@ -1531,7 +1516,7 @@
  Marc shared lib      : $progdll
  Version type         : $mode
  Job ID               : $DIRJID/$jid$extra_job_info
@@ -167,16 +115,7 @@ index 88a3027..85beabe 100644
  User objects/libs    : $objs
  Restart file job ID  : $rid
  Substructure file ID : $sid
-@@ -1548,8 +1533,6 @@ GPGPU option         : $gpuids
- Host file name       : $host" > $jid.log
- if test "$iprintsimufact" = true ; then
-  echo "DDM with ARC Mapper  : $ddm_arc" >> $jid.log
-- echo "Mapper directory     : $MESHERDIR" >> $jid.log
-- echo "Meshing directory    : $MESHERDIR" >> $jid.log
- fi
- echo \
- "Message passing type : $itree
-@@ -1564,7 +1547,7 @@ Program name         : $prog
+@@ -1564,7 +1549,7 @@
  Marc shared lib      : $progdll
  Version type         : $mode
  Job ID               : $DIRJID/$jid$extra_job_info
@@ -185,16 +124,7 @@ index 88a3027..85beabe 100644
  User objects/libs    : $objs
  Restart file job ID  : $rid
  Substructure file ID : $sid
-@@ -1579,6 +1562,8 @@ Solver processes     : $nsolverprint
- Solver threads       : $ntsprint"
- if test "$iprintsimufact" = true ; then
-  echo "DDM with ARC Mapper  : $ddm_arc"
-+ echo "Mapper directory     : $MESHERDIR" >> $jid.log
-+ echo "Meshing directory    : $MESHERDIR" >> $jid.log
- fi
- echo \
- "GPGPU option         : $gpuids
-@@ -1687,7 +1672,7 @@ Program name ($prog)? $ECHOTXT"
+@@ -1687,7 +1672,7 @@
                      ;;
                    esac        
                  fi
@@ -203,7 +133,7 @@ index 88a3027..85beabe 100644
  		read value
  		if test "$value"
  		then
-@@ -1696,50 +1681,19 @@ Program name ($prog)? $ECHOTXT"
+@@ -1696,50 +1681,19 @@
                        user=
                      ;;
                      *)
@@ -229,16 +159,13 @@ index 88a3027..85beabe 100644
 -                      fi
 +                      user=$value
                        case $user in
--                         \/*)
--                           ;;
--                         *)
-+                        \/*)
-+                        ;;
-+                        *)
+                          \/*)
+                            ;;
+                          *)
                            user=`pwd`/$user
 -                          usersubname=`pwd`/$usersubname
--                           ;;
--                      esac
+                            ;;
+                       esac
 -                      if test ! -f $usersubname
 -                      then
 -                       if test -f $usersubname.f 
@@ -255,8 +182,6 @@ index 88a3027..85beabe 100644
 -                        usersubname=$usersubname.F90
 -                       fi
 -                      fi
-+                        ;;
-+                      esac        
 +					  usernoext=$user
 +					  usernoext=`dirname $usernoext`/`$BASENAME $usernoext .f`
 +					  usernoext=`dirname $usernoext`/`$BASENAME $usernoext .F`
@@ -265,7 +190,7 @@ index 88a3027..85beabe 100644
                      ;;
                    esac
  		fi
-@@ -2274,11 +2228,12 @@ fi
+@@ -2274,11 +2228,12 @@
  #
  # user subroutine used
  #
@@ -280,7 +205,7 @@ index 88a3027..85beabe 100644
      case $program in
        \/* | \.\/*)
          bd=
-@@ -2391,7 +2346,7 @@ fi
+@@ -2391,7 +2346,7 @@
      fi
      if test "$user"
      then
@@ -289,7 +214,7 @@ index 88a3027..85beabe 100644
       usersub=1
      fi
      export execpath      
-@@ -3274,44 +3229,27 @@ then
+@@ -3274,44 +3229,27 @@
      echo
      if test "$user"
      then
@@ -315,7 +240,7 @@ index 88a3027..85beabe 100644
          if test $MACHINENAME = "CRAY"
          then
 -            $FORTRAN $usersub || \
-+            $DFORTRANMP $user || \
++            $DFORTHIGHMP $user || \
                  {
 -                echo "$0: compile failed for $user.f"
 +                echo "$0: compile failed for $user"
@@ -324,7 +249,7 @@ index 88a3027..85beabe 100644
              /bin/rm $program 2>/dev/null
          else
 -            $FORTRAN $usersub -o $userobj || \
-+            $DFORTRANMP $user -o $userobj || \
++            $DFORTHIGHMP $user -o $userobj || \
                  {
 -                echo "$0: compile failed for $user.f"
 +                echo "$0: compile failed for $user"
@@ -339,21 +264,15 @@ index 88a3027..85beabe 100644
      fi
  
  
-@@ -3331,10 +3269,11 @@ then
+@@ -3331,6 +3269,7 @@
      $TKLIBS  \
      $MRCLIBS     \
      $METISLIBS \
--    $SFLIB \
 +    $DAMASK   \
-+    $SFLIB  \
+     $SFLIB \
      $OPENSSL_LIB \
      $SYSLIBS  \
--    $SECLIBS  || \
-+    $SECLIBS      || \
-     {
-     echo "$0: link failed for ${user:+$userobj }$objs"
-     exit 1
-@@ -3344,6 +3283,9 @@ else
+@@ -3344,6 +3283,9 @@
  	prgsav=yes
  fi
  /bin/rm $userobj 2>/dev/null
@@ -363,7 +282,7 @@ index 88a3027..85beabe 100644
  
  #
  # run marc
-@@ -3390,7 +3332,7 @@ if test $dllrun -eq 0; then
+@@ -3390,7 +3332,7 @@
    fi
  else
    if test $cpdll = yes; then
@@ -372,7 +291,7 @@ index 88a3027..85beabe 100644
       /bin/cp $DIRJOB/$marcdll $DIRJOB/${filename}_$marcdll 2>/dev/null
    fi
    if test $rmdll = yes
-@@ -3556,7 +3498,7 @@ then
+@@ -3556,7 +3498,7 @@
                # first copy over the user sub if local directories
                if test ${dirstatus[$counter]} = "local"
                then
@@ -381,12 +300,12 @@ index 88a3027..85beabe 100644
                fi
                # do the compilation on the other machine
                if test ${dirstatus[$counter]} = "shared"
-@@ -3569,21 +3511,21 @@ then
+@@ -3569,21 +3511,21 @@
                remoteuser=$DIR1/`$BASENAME $user`
                $RSH $i /bin/rm $remoteprog 2> /dev/null
                echo
 -              $RSH $i $DIR2/tools/comp_user $DIR2 $DIR1 $remoteuser $remoteprog
-+              $RSH $i $DIR2/tools/comp_damask_mp $DIR2 $DIR1 $remoteuser $remoteprog
++              $RSH $i $DIR2/tools/comp_damask_hmp $DIR2 $DIR1 $remoteuser $remoteprog
                # check if successful, the new executable should be there
                line=`$RSH $i /bin/ls $remoteprog 2> /dev/null`
                if test "$line"
@@ -406,7 +325,7 @@ index 88a3027..85beabe 100644
                fi
              fi
            fi
-@@ -3593,39 +3535,27 @@ then
+@@ -3593,39 +3535,27 @@
      if test "$userhost"
      then
        echo
@@ -426,7 +345,7 @@ index 88a3027..85beabe 100644
          if test $MACHINENAME = "CRAY"
          then
 -            $FORTRAN $usersub || \
-+            $DFORTRANMP $user || \
++            $DFORTHIGHMP $user || \
                  {
 -                echo "$0: compile failed for $user.f"
 +                echo "$0: compile failed for $user"
@@ -436,7 +355,7 @@ index 88a3027..85beabe 100644
              /bin/rm $program 2>/dev/null
          else
 -            $FORTRAN $usersub -o $userobj || \
-+            $DFORTRANMP $user -o $userobj || \
++            $DFORTHIGHMP $user -o $userobj || \
                  {
 -                echo "$0: compile failed for $user.f"
 +                echo "$0: compile failed for $user"
@@ -452,21 +371,15 @@ index 88a3027..85beabe 100644
    fi # if test $user
  
  
-@@ -3645,10 +3575,11 @@ then
+@@ -3645,6 +3575,7 @@
    $TKLIBS  \
    $MRCLIBS     \
    $METISLIBS \
--  $SFLIB \
 +  $DAMASK \
-+  $SFLIB  \
+   $SFLIB \
    $OPENSSL_LIB \
    $SYSLIBS  \
--  $SECLIBS  || \
-+  $SECLIBS    || \
-   {
-     echo "$0: link failed for ${user:+$userobj }$objs"
-     echo "         $PRODUCT Exit number 3"
-@@ -3686,6 +3617,9 @@ else # if test $link
+@@ -3686,6 +3617,9 @@
  	prgsav=yes
  fi   # if test $link
  /bin/rm $userobj 2>/dev/null
@@ -476,69 +389,7 @@ index 88a3027..85beabe 100644
  
  #
  # run marc
-@@ -3744,42 +3678,42 @@ then
-       counter=0
-       if test -f "$host_filt"
-       then
--        for i in `$AWK -v n=$numfield '{print $n}' $host_filt`
--        do
--          ibase=${i%%.*}
--          if test $ibase != $thishost
-+      for i in `$AWK -v n=$numfield '{print $n}' $host_filt`
-+      do
-+        ibase=${i%%.*}
-+        if test $ibase != $thishost
-+        then
-+          counter=$((counter+1))
-+          DIR1=$DIRJOB
-+          line=`grep -v '^#' $userhost | grep "^$ibase "`
-+          workdir=`echo $line | $AWK '{print $3}'`
-+          if test -n "$workdir"
-           then
--            counter=$((counter+1))
--            DIR1=$DIRJOB
--            line=`grep -v '^#' $userhost | grep "^$ibase "`
--            workdir=`echo $line | $AWK '{print $3}'`
--            if test -n "$workdir"
--            then
--              DIR1=$workdir
--            fi
--            #  if an incompatible host uses shared directory,
--            #  then the root machine deletes the executable
--            if test ${dirstatus[$counter]} = "shared" -a  ${compstatus[$counter]} = "no"
--            then
--              hname=_$ibase
--              /bin/rm ${execname}$hname
--            fi
--            # if local directory used, the remote machine
--            # deletes the executable
--            if test ${dirstatus[$counter]} = "local"
--            then
--              $RSH $i /bin/rm $DIR1/${execname} 2>/dev/null
--            fi
-+            DIR1=$workdir
-           fi
--        done
--      fi
-+          #  if an incompatible host uses shared directory,
-+          #  then the root machine deletes the executable
-+          if test ${dirstatus[$counter]} = "shared" -a  ${compstatus[$counter]} = "no"
-+          then
-+            hname=_$ibase
-+            /bin/rm ${execname}$hname
-+          fi
-+          # if local directory used, the remote machine
-+          # deletes the executable
-+          if test ${dirstatus[$counter]} = "local"
-+          then
-+            $RSH $i /bin/rm $DIR1/${execname} 2>/dev/null
-+          fi
-+        fi
-+      done
-     fi
-   fi
- fi
-+fi
+@@ -3779,7 +3713,7 @@
  else
  #dllrun >0
    if test $cpdll = yes; then
@@ -547,7 +398,7 @@ index 88a3027..85beabe 100644
       /bin/cp $DIRJOB/$marcdll $DIRJOB/${filename}_$marcdll 2>/dev/null
    fi
    if test $rmdll = yes;then
-@@ -3904,7 +3838,7 @@ then
+@@ -3904,7 +3838,7 @@
                # first copy over the user sub if local directories
                if test ${dirstatus[$counter]} = "local"
                then
@@ -556,12 +407,12 @@ index 88a3027..85beabe 100644
                fi
                # do the compilation on the other machine
                if test ${dirstatus[$counter]} = "shared"
-@@ -3917,20 +3851,20 @@ then
+@@ -3917,20 +3851,20 @@
                remoteuser=$DIR1/`$BASENAME $user`
                $RSH $i /bin/rm $remoteprog 2> /dev/null
                echo
 -              $RSH $i $DIR2/tools/comp_user $DIR2 $DIR1 $remoteuser $remoteprog
-+              $RSH $i $DIR2/tools/comp_damask_mp $DIR2 $DIR1 $remoteuser $remoteprog
++              $RSH $i $DIR2/tools/comp_damask_hmp $DIR2 $DIR1 $remoteuser $remoteprog
                # check if successful, the new executable should be there
                line=`$RSH $i /bin/ls $remoteprog 2> /dev/null`
                if test "$line"
@@ -580,7 +431,7 @@ index 88a3027..85beabe 100644
                fi
              fi
            fi
-@@ -3940,37 +3874,25 @@ then
+@@ -3940,37 +3874,25 @@
      if test "$userhost"
      then
        echo
@@ -600,7 +451,7 @@ index 88a3027..85beabe 100644
          if test $MACHINENAME = "CRAY"
          then
 -            $FORTRAN $usersub || \
-+            $DFORTRANMP $user || \
++            $DFORTHIGHMP $user || \
                  {
 -                echo "$0: compile failed for $user.f"
 +                echo "$0: compile failed for $user"
@@ -609,7 +460,7 @@ index 88a3027..85beabe 100644
              /bin/rm $program 2>/dev/null
          else
 -            $FORTRAN $usersub -o $userobj || \
-+            $DFORTRANMP $user -o $userobj || \
++            $DFORTHIGHMP $user -o $userobj || \
                  {
 -                echo "$0: compile failed for $user.f"
 +                echo "$0: compile failed for $user"
@@ -624,7 +475,7 @@ index 88a3027..85beabe 100644
    fi # if test $user
  
  
-@@ -3990,10 +3912,11 @@ then
+@@ -3990,6 +3912,7 @@
    $TKLIBS  \
    $MRCLIBS     \
    $METISLIBS \
@@ -632,12 +483,7 @@ index 88a3027..85beabe 100644
    $SFLIB \
    $OPENSSL_LIB \
    $SYSLIBS \
--  $SECLIBS  || \
-+  $SECLIBS         || \
-   {
-     echo "$0: link failed for ${user:+$userobj }$objs"
-     exit 1
-@@ -4030,7 +3953,9 @@ else # if test $link
+@@ -4030,7 +3953,9 @@
  	prgsav=yes
  fi   # if test $link
  /bin/rm $userobj 2>/dev/null
@@ -648,78 +494,7 @@ index 88a3027..85beabe 100644
  # done if no job id given
  if test -z "$jid"
  then
-@@ -4070,7 +3995,7 @@ if test $ddm_arc -gt 0; then
-   RUN_JOB="$BINDIR/exeddm $RUN_JOB -ddm $ddm_arc "
- fi
- 
--$RUN_JOB 
-+    $RUN_JOB 
- 
- if test $nprocd -gt 1
- then
-@@ -4114,42 +4039,42 @@ then
-       counter=0
-       if test -f "$host_filt"
-       then
--        for i in `$AWK -v n=$numfield '{print $n}' $host_filt`
--        do
--          ibase=${i%%.*}
--          if test $ibase != $thishost
-+      for i in `$AWK -v n=$numfield '{print $n}' $host_filt`
-+      do
-+        ibase=${i%%.*}
-+        if test $ibase != $thishost
-+        then
-+          counter=$((counter+1))
-+          DIR1=$DIRJOB
-+          line=`grep -v '^#' $userhost | grep "^$ibase "`
-+          workdir=`echo $line | $AWK '{print $3}'`
-+          if test -n "$workdir"
-           then
--            counter=$((counter+1))
--            DIR1=$DIRJOB
--            line=`grep -v '^#' $userhost | grep "^$ibase "`
--            workdir=`echo $line | $AWK '{print $3}'`
--            if test -n "$workdir"
--            then
--              DIR1=$workdir
--            fi
--            #  if an incompatible host uses shared directory,
--            #  then the root machine deletes the executable
--            if test ${dirstatus[$counter]} = "shared" -a  ${compstatus[$counter]} = "no"
--            then
--              hname=_$ibase
--              /bin/rm ${execname}$hname
--            fi
--            # if local directory used, the remote machine
--            # deletes the executable
--            if test ${dirstatus[$counter]} = "local"
--            then
--              $RSH $i /bin/rm $DIR1/${execname} 2>/dev/null
--            fi
-+            DIR1=$workdir
-           fi
--        done
--      fi
-+          #  if an incompatible host uses shared directory,
-+          #  then the root machine deletes the executable
-+          if test ${dirstatus[$counter]} = "shared" -a  ${compstatus[$counter]} = "no"
-+          then
-+            hname=_$ibase
-+            /bin/rm ${execname}$hname
-+          fi
-+          # if local directory used, the remote machine
-+          # deletes the executable
-+          if test ${dirstatus[$counter]} = "local"
-+          then
-+            $RSH $i /bin/rm $DIR1/${execname} 2>/dev/null
-+          fi
-+        fi
-+      done
-     fi
-   fi
- fi
-+fi
+@@ -4149,7 +4074,7 @@
  else
  #dllrun >0
    if test $cpdll = yes; then
diff --git a/install/MarcMentat/2022.1/Mentat_bin/submit4.patch b/install/MarcMentat/2022.1/Mentat_bin/submit4.patch
index 3f1371fdb..98c51e76d 100644
--- a/install/MarcMentat/2022.1/Mentat_bin/submit4.patch
+++ b/install/MarcMentat/2022.1/Mentat_bin/submit4.patch
@@ -1,6 +1,6 @@
 ---
 +++
-@@ -63,10 +64,10 @@ doe_nparallel=$6
+@@ -63,10 +63,10 @@
  if [ "$slv" != "" -a "$slv" != "marc" -a "$slv" != "datfit" ]; then
    slv="-iam sfm"
  fi
@@ -13,7 +13,7 @@
    slv="-iam datfit"
  fi
  
-@@ -91,6 +92,7 @@ if [ "$srcfile" != "" -a "$srcfile" != "-" ]; then
+@@ -91,6 +91,7 @@
        srcfile="-u $srcfile -save y"
        ;;
      runsaved)
@@ -21,7 +21,7 @@
        srcfile="-prog $srcfile"
        ;;
    esac
-@@ -189,12 +191,12 @@ unset PYTHONHOME
+@@ -189,12 +190,12 @@
  unset PYTHONPATH
  
  if [ "$doe_first" = "-" ]; then # submit of regular Marc job
diff --git a/install/MarcMentat/2022.1/Mentat_bin/submit5.patch b/install/MarcMentat/2022.1/Mentat_bin/submit5.patch
index 9614d8c69..ab32b1058 100644
--- a/install/MarcMentat/2022.1/Mentat_bin/submit5.patch
+++ b/install/MarcMentat/2022.1/Mentat_bin/submit5.patch
@@ -1,6 +1,6 @@
 ---
 +++
-@@ -63,10 +64,10 @@ doe_nparallel=$6
+@@ -63,10 +63,10 @@
  if [ "$slv" != "" -a "$slv" != "marc" -a "$slv" != "datfit" ]; then
    slv="-iam sfm"
  fi
@@ -13,7 +13,7 @@
    slv="-iam datfit"
  fi
  
-@@ -91,6 +92,7 @@ if [ "$srcfile" != "" -a "$srcfile" != "-" ]; then
+@@ -91,6 +91,7 @@
        srcfile="-u $srcfile -save y"
        ;;
      runsaved)
@@ -21,7 +21,7 @@
        srcfile="-prog $srcfile"
        ;;
    esac
-@@ -189,12 +191,12 @@ unset PYTHONHOME
+@@ -189,12 +190,12 @@
  unset PYTHONPATH
  
  if [ "$doe_first" = "-" ]; then # submit of regular Marc job
diff --git a/install/MarcMentat/2022.1/Mentat_bin/submit6.patch b/install/MarcMentat/2022.1/Mentat_bin/submit6.patch
index a3ed16135..d5ea3cfde 100644
--- a/install/MarcMentat/2022.1/Mentat_bin/submit6.patch
+++ b/install/MarcMentat/2022.1/Mentat_bin/submit6.patch
@@ -1,6 +1,6 @@
 ---
 +++
-@@ -63,10 +64,10 @@ doe_nparallel=$6
+@@ -63,10 +63,10 @@
  if [ "$slv" != "" -a "$slv" != "marc" -a "$slv" != "datfit" ]; then
    slv="-iam sfm"
  fi
@@ -13,7 +13,7 @@
    slv="-iam datfit"
  fi
  
-@@ -91,6 +92,7 @@ if [ "$srcfile" != "" -a "$srcfile" != "-" ]; then
+@@ -91,6 +91,7 @@
        srcfile="-u $srcfile -save y"
        ;;
      runsaved)
@@ -21,7 +21,7 @@
        srcfile="-prog $srcfile"
        ;;
    esac
-@@ -189,12 +191,12 @@ unset PYTHONHOME
+@@ -189,12 +190,12 @@
  unset PYTHONPATH
  
  if [ "$doe_first" = "-" ]; then # submit of regular Marc job
diff --git a/install/MarcMentat/2022.1/Mentat_menus/job_run.ms.patch b/install/MarcMentat/2022.1/Mentat_menus/job_run.ms.patch
index 6cbcf895c..9ff6b8ae0 100644
--- a/install/MarcMentat/2022.1/Mentat_menus/job_run.ms.patch
+++ b/install/MarcMentat/2022.1/Mentat_menus/job_run.ms.patch
@@ -1,6 +1,6 @@
 ---
 +++
-@@ -261,11 +261,18 @@ popmenu job_run_popmenu {
+@@ -261,12 +261,19 @@
        }
        button {
         position +25 =
@@ -10,21 +10,21 @@
         help "job_run#Job Submission And Control"
         popmenu job_submit_adv_pm
        }
-+      button {
+       button {
 +       position +18 =
 +       size 7 4
 +       text "DAMASK"
 +       help "damask_run#Job Submission And Control"
 +       popmenu damask
 +      }
-       button {
++      button {
         position 0 +4
         size 12 4
-@@ -1189,6 +1196,135 @@ popmenu job_submit_adv_pm {
- }
+        text "UPDATE"
+@@ -1190,6 +1197,135 @@
  
  
-+#--------------------------------------------------------------------------------------------------
+ #--------------------------------------------------------------------------------------------------
 +popmenu damask {
 +
 +#ifdef QT_MENTAT
@@ -153,6 +153,7 @@
 + mode permanent
 +}
 +
- #--------------------------------------------------------------------------------------------------
++#--------------------------------------------------------------------------------------------------
  popmenu job_exit_msg_pm {
  
+  text "EXIT MESSAGE"
diff --git a/python/setup.cfg b/python/setup.cfg
index 321040658..e183e6b3f 100644
--- a/python/setup.cfg
+++ b/python/setup.cfg
@@ -6,7 +6,7 @@ author_email = damask@mpie.de
 url = https://damask.mpie.de
 description = DAMASK processing tools
 long_description = Pre- and post-processing tools for DAMASK
-license: AGPL3
+license = AGPL3
 classifiers =
     Intended Audience :: Science/Research
     Topic :: Scientific/Engineering
@@ -19,12 +19,11 @@ packages = find:
 include_package_data = true
 python_requires = >= 3.8
 install_requires =
-    importlib-metadata; python_version<"3.8"
-    pandas; python_version<="0.24"                                                                  # requires numpy
-    numpy; python_version<="1.17"                                                                   # needed for default_rng
-    scipy; python_version<="1.2"
-    h5py; python_version<="2.9"                                                                     # requires numpy
-    vtk; python_version<="8.1"
-    matplotlib; python_version<="3.0"                                                               # requires numpy, pillow
-    pyyaml; python_version<="3.12"
+    pandas>=0.24                                                                                    # requires numpy
+    numpy>=1.17                                                                                     # needed for default_rng
+    scipy>=1.2
+    h5py>=2.9                                                                                       # requires numpy
+    vtk>=8.1
+    matplotlib>=3.0                                                                                 # requires numpy, pillow
+    pyyaml>=3.12
 setup_requires = setuptools
diff --git a/src/CLI.f90 b/src/CLI.f90
index ee1555cb2..4fef460ff 100644
--- a/src/CLI.f90
+++ b/src/CLI.f90
@@ -7,7 +7,7 @@
 !--------------------------------------------------------------------------------------------------
 #define PETSC_MAJOR 3
 #define PETSC_MINOR_MIN 12
-#define PETSC_MINOR_MAX 17
+#define PETSC_MINOR_MAX 18
 
 module CLI
   use, intrinsic :: ISO_fortran_env
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index a473069b9..d0a95fa23 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -1,7 +1,8 @@
 # special flags for some files
 if(CMAKE_Fortran_COMPILER_ID STREQUAL "GNU")
     # long lines for interaction matrix
-    SET_SOURCE_FILES_PROPERTIES("lattice.f90" PROPERTIES COMPILE_FLAGS "-ffree-line-length-240")
+    set_source_files_properties("lattice.f90" PROPERTIES COMPILE_FLAGS "-ffree-line-length-240")
+    set_source_files_properties("parallelization.f90" PROPERTIES COMPILE_FLAGS "-ffree-line-length-none")
 endif()
 
 file(GLOB damask-sources CONFIGURE_DEPENDS *.f90 *.c)
@@ -14,6 +15,15 @@ elseif(PROJECT_NAME STREQUAL "damask-mesh")
     file(GLOB solver-sources CONFIGURE_DEPENDS mesh/*.f90)
 endif()
 
+foreach(solver-source ${solver-sources})
+    file(READ ${solver-source} content)
+    string(FIND "${content}" "CHKERR" found)
+    if(NOT ${found} EQUAL -1)
+        set_source_files_properties(${solver-source} PROPERTIES COMPILE_FLAGS "-ffree-line-length-none")
+    endif()
+endforeach()
+
+
 if(NOT CMAKE_BUILD_TYPE STREQUAL "SYNTAXONLY")
     add_executable(${executable-name} ${damask-sources} ${solver-sources})
     install(TARGETS ${executable-name} RUNTIME DESTINATION bin)
diff --git a/src/HDF5_utilities.f90 b/src/HDF5_utilities.f90
index ee1e3cf8c..4829d0729 100644
--- a/src/HDF5_utilities.f90
+++ b/src/HDF5_utilities.f90
@@ -673,13 +673,9 @@ subroutine HDF5_read_real1(dataset,loc_id,datasetName,parallel)
     totalShape                                                                                      !< shape of the dataset (all processes)
   integer :: hdferr
 
-!---------------------------------------------------------------------------------------------------
-! determine shape of dataset
-  myShape = int(shape(dataset),HSIZE_T)
-  if (any(myShape(1:size(myShape)-1) == 0)) return                                                  !< empty dataset (last dimension can be empty)
 
-!---------------------------------------------------------------------------------------------------
-! initialize HDF5 data structures
+  myShape = int(shape(dataset),HSIZE_T)
+
   if (present(parallel)) then
     call initialize_read(dset_id, filespace_id, memspace_id, plist_id, aplist_id, &
                          myStart, totalShape, loc_id,myShape,datasetName,parallel)
@@ -687,6 +683,7 @@ subroutine HDF5_read_real1(dataset,loc_id,datasetName,parallel)
     call initialize_read(dset_id, filespace_id, memspace_id, plist_id, aplist_id, &
                          myStart, totalShape, loc_id,myShape,datasetName,parallel_default)
   end if
+  if (any(totalShape == 0)) return
 
   call H5Dread_f(dset_id, H5T_NATIVE_DOUBLE,dataset,totalShape, hdferr,&
                  file_space_id = filespace_id, xfer_prp = plist_id, mem_space_id = memspace_id)
@@ -713,13 +710,9 @@ subroutine HDF5_read_real2(dataset,loc_id,datasetName,parallel)
     totalShape                                                                                      !< shape of the dataset (all processes)
   integer :: hdferr
 
-!---------------------------------------------------------------------------------------------------
-! determine shape of dataset
-  myShape = int(shape(dataset),HSIZE_T)
-  if (any(myShape(1:size(myShape)-1) == 0)) return                                                  !< empty dataset (last dimension can be empty)
 
-!---------------------------------------------------------------------------------------------------
-! initialize HDF5 data structures
+  myShape = int(shape(dataset),HSIZE_T)
+
   if (present(parallel)) then
     call initialize_read(dset_id, filespace_id, memspace_id, plist_id, aplist_id, &
                          myStart, totalShape, loc_id,myShape,datasetName,parallel)
@@ -727,6 +720,7 @@ subroutine HDF5_read_real2(dataset,loc_id,datasetName,parallel)
     call initialize_read(dset_id, filespace_id, memspace_id, plist_id, aplist_id, &
                          myStart, totalShape, loc_id,myShape,datasetName,parallel_default)
   end if
+  if (any(totalShape == 0)) return
 
   call H5Dread_f(dset_id, H5T_NATIVE_DOUBLE,dataset,totalShape, hdferr,&
                  file_space_id = filespace_id, xfer_prp = plist_id, mem_space_id = memspace_id)
@@ -753,13 +747,9 @@ subroutine HDF5_read_real3(dataset,loc_id,datasetName,parallel)
     totalShape                                                                                      !< shape of the dataset (all processes)
   integer :: hdferr
 
-!---------------------------------------------------------------------------------------------------
-! determine shape of dataset
-  myShape = int(shape(dataset),HSIZE_T)
-  if (any(myShape(1:size(myShape)-1) == 0)) return                                                  !< empty dataset (last dimension can be empty)
 
-!---------------------------------------------------------------------------------------------------
-! initialize HDF5 data structures
+  myShape = int(shape(dataset),HSIZE_T)
+
   if (present(parallel)) then
     call initialize_read(dset_id, filespace_id, memspace_id, plist_id, aplist_id, &
                          myStart, totalShape, loc_id,myShape,datasetName,parallel)
@@ -767,6 +757,7 @@ subroutine HDF5_read_real3(dataset,loc_id,datasetName,parallel)
     call initialize_read(dset_id, filespace_id, memspace_id, plist_id, aplist_id, &
                          myStart, totalShape, loc_id,myShape,datasetName,parallel_default)
   end if
+  if (any(totalShape == 0)) return
 
   call H5Dread_f(dset_id, H5T_NATIVE_DOUBLE,dataset,totalShape, hdferr,&
                  file_space_id = filespace_id, xfer_prp = plist_id, mem_space_id = memspace_id)
@@ -793,13 +784,10 @@ subroutine HDF5_read_real4(dataset,loc_id,datasetName,parallel)
     totalShape                                                                                      !< shape of the dataset (all processes)
   integer :: hdferr
 
-!---------------------------------------------------------------------------------------------------
-! determine shape of dataset
+
   myShape = int(shape(dataset),HSIZE_T)
   if (any(myShape(1:size(myShape)-1) == 0)) return                                                  !< empty dataset (last dimension can be empty)
 
-!---------------------------------------------------------------------------------------------------
-! initialize HDF5 data structures
   if (present(parallel)) then
     call initialize_read(dset_id, filespace_id, memspace_id, plist_id, aplist_id, &
                          myStart, totalShape, loc_id,myShape,datasetName,parallel)
@@ -807,6 +795,7 @@ subroutine HDF5_read_real4(dataset,loc_id,datasetName,parallel)
     call initialize_read(dset_id, filespace_id, memspace_id, plist_id, aplist_id, &
                          myStart, totalShape, loc_id,myShape,datasetName,parallel_default)
   end if
+  if (any(totalShape == 0)) return
 
   call H5Dread_f(dset_id, H5T_NATIVE_DOUBLE,dataset,totalShape, hdferr,&
                  file_space_id = filespace_id, xfer_prp = plist_id, mem_space_id = memspace_id)
@@ -833,13 +822,10 @@ subroutine HDF5_read_real5(dataset,loc_id,datasetName,parallel)
     totalShape                                                                                      !< shape of the dataset (all processes)
   integer :: hdferr
 
-!---------------------------------------------------------------------------------------------------
-! determine shape of dataset
+
   myShape = int(shape(dataset),HSIZE_T)
   if (any(myShape(1:size(myShape)-1) == 0)) return                                                  !< empty dataset (last dimension can be empty)
 
-!---------------------------------------------------------------------------------------------------
-! initialize HDF5 data structures
   if (present(parallel)) then
     call initialize_read(dset_id, filespace_id, memspace_id, plist_id, aplist_id, &
                          myStart, totalShape, loc_id,myShape,datasetName,parallel)
@@ -847,6 +833,7 @@ subroutine HDF5_read_real5(dataset,loc_id,datasetName,parallel)
     call initialize_read(dset_id, filespace_id, memspace_id, plist_id, aplist_id, &
                          myStart, totalShape, loc_id,myShape,datasetName,parallel_default)
   end if
+  if (any(totalShape == 0)) return
 
   call H5Dread_f(dset_id, H5T_NATIVE_DOUBLE,dataset,totalShape, hdferr,&
                  file_space_id = filespace_id, xfer_prp = plist_id, mem_space_id = memspace_id)
@@ -873,13 +860,10 @@ subroutine HDF5_read_real6(dataset,loc_id,datasetName,parallel)
     totalShape                                                                                      !< shape of the dataset (all processes)
   integer :: hdferr
 
-!---------------------------------------------------------------------------------------------------
-! determine shape of dataset
+
   myShape = int(shape(dataset),HSIZE_T)
   if (any(myShape(1:size(myShape)-1) == 0)) return                                                  !< empty dataset (last dimension can be empty)
 
-!---------------------------------------------------------------------------------------------------
-! initialize HDF5 data structures
   if (present(parallel)) then
     call initialize_read(dset_id, filespace_id, memspace_id, plist_id, aplist_id, &
                          myStart, totalShape, loc_id,myShape,datasetName,parallel)
@@ -887,6 +871,7 @@ subroutine HDF5_read_real6(dataset,loc_id,datasetName,parallel)
     call initialize_read(dset_id, filespace_id, memspace_id, plist_id, aplist_id, &
                          myStart, totalShape, loc_id,myShape,datasetName,parallel_default)
   end if
+  if (any(totalShape == 0)) return
 
   call H5Dread_f(dset_id, H5T_NATIVE_DOUBLE,dataset,totalShape, hdferr,&
                  file_space_id = filespace_id, xfer_prp = plist_id, mem_space_id = memspace_id)
@@ -913,13 +898,10 @@ subroutine HDF5_read_real7(dataset,loc_id,datasetName,parallel)
     totalShape                                                                                      !< shape of the dataset (all processes)
   integer :: hdferr
 
-!---------------------------------------------------------------------------------------------------
-! determine shape of dataset
+
   myShape = int(shape(dataset),HSIZE_T)
   if (any(myShape(1:size(myShape)-1) == 0)) return                                                  !< empty dataset (last dimension can be empty)
 
-!---------------------------------------------------------------------------------------------------
-! initialize HDF5 data structures
   if (present(parallel)) then
     call initialize_read(dset_id, filespace_id, memspace_id, plist_id, aplist_id, &
                          myStart, totalShape, loc_id,myShape,datasetName,parallel)
@@ -927,6 +909,7 @@ subroutine HDF5_read_real7(dataset,loc_id,datasetName,parallel)
     call initialize_read(dset_id, filespace_id, memspace_id, plist_id, aplist_id, &
                          myStart, totalShape, loc_id,myShape,datasetName,parallel_default)
   end if
+  if (any(totalShape == 0)) return
 
   call H5Dread_f(dset_id, H5T_NATIVE_DOUBLE,dataset,totalShape, hdferr,&
                  file_space_id = filespace_id, xfer_prp = plist_id, mem_space_id = memspace_id)
@@ -955,13 +938,8 @@ subroutine HDF5_read_int1(dataset,loc_id,datasetName,parallel)
     totalShape                                                                                      !< shape of the dataset (all processes)
   integer :: hdferr
 
-!---------------------------------------------------------------------------------------------------
-! determine shape of dataset
   myShape = int(shape(dataset),HSIZE_T)
-  if (any(myShape(1:size(myShape)-1) == 0)) return                                                  !< empty dataset (last dimension can be empty)
 
-!---------------------------------------------------------------------------------------------------
-! initialize HDF5 data structures
   if (present(parallel)) then
     call initialize_read(dset_id, filespace_id, memspace_id, plist_id, aplist_id, &
                          myStart, totalShape, loc_id,myShape,datasetName,parallel)
@@ -969,6 +947,7 @@ subroutine HDF5_read_int1(dataset,loc_id,datasetName,parallel)
     call initialize_read(dset_id, filespace_id, memspace_id, plist_id, aplist_id, &
                          myStart, totalShape, loc_id,myShape,datasetName,parallel_default)
   end if
+  if (any(totalShape == 0)) return
 
   call H5Dread_f(dset_id, H5T_NATIVE_INTEGER,dataset,totalShape, hdferr,&
                  file_space_id = filespace_id, xfer_prp = plist_id, mem_space_id = memspace_id)
@@ -995,13 +974,10 @@ subroutine HDF5_read_int2(dataset,loc_id,datasetName,parallel)
     totalShape                                                                                      !< shape of the dataset (all processes)
   integer :: hdferr
 
-!---------------------------------------------------------------------------------------------------
-! determine shape of dataset
+
   myShape = int(shape(dataset),HSIZE_T)
   if (any(myShape(1:size(myShape)-1) == 0)) return                                                  !< empty dataset (last dimension can be empty)
 
-!---------------------------------------------------------------------------------------------------
-! initialize HDF5 data structures
   if (present(parallel)) then
     call initialize_read(dset_id, filespace_id, memspace_id, plist_id, aplist_id, &
                          myStart, totalShape, loc_id,myShape,datasetName,parallel)
@@ -1009,6 +985,7 @@ subroutine HDF5_read_int2(dataset,loc_id,datasetName,parallel)
     call initialize_read(dset_id, filespace_id, memspace_id, plist_id, aplist_id, &
                          myStart, totalShape, loc_id,myShape,datasetName,parallel_default)
   end if
+  if (any(totalShape == 0)) return
 
   call H5Dread_f(dset_id, H5T_NATIVE_INTEGER,dataset,totalShape, hdferr,&
                  file_space_id = filespace_id, xfer_prp = plist_id, mem_space_id = memspace_id)
@@ -1035,13 +1012,9 @@ subroutine HDF5_read_int3(dataset,loc_id,datasetName,parallel)
     totalShape                                                                                      !< shape of the dataset (all processes)
   integer :: hdferr
 
-!---------------------------------------------------------------------------------------------------
-! determine shape of dataset
-  myShape = int(shape(dataset),HSIZE_T)
-  if (any(myShape(1:size(myShape)-1) == 0)) return                                                  !< empty dataset (last dimension can be empty)
 
-!---------------------------------------------------------------------------------------------------
-! initialize HDF5 data structures
+  myShape = int(shape(dataset),HSIZE_T)
+
   if (present(parallel)) then
     call initialize_read(dset_id, filespace_id, memspace_id, plist_id, aplist_id, &
                          myStart, totalShape, loc_id,myShape,datasetName,parallel)
@@ -1049,6 +1022,7 @@ subroutine HDF5_read_int3(dataset,loc_id,datasetName,parallel)
    call initialize_read(dset_id, filespace_id, memspace_id, plist_id, aplist_id, &
                         myStart, totalShape, loc_id,myShape,datasetName,parallel_default)
   end if
+  if (any(totalShape == 0)) return
 
   call H5Dread_f(dset_id, H5T_NATIVE_INTEGER,dataset,totalShape, hdferr,&
                  file_space_id = filespace_id, xfer_prp = plist_id, mem_space_id = memspace_id)
@@ -1075,13 +1049,9 @@ subroutine HDF5_read_int4(dataset,loc_id,datasetName,parallel)
     totalShape                                                                                      !< shape of the dataset (all processes)
   integer :: hdferr
 
-!---------------------------------------------------------------------------------------------------
-! determine shape of dataset
-  myShape = int(shape(dataset),HSIZE_T)
-  if (any(myShape(1:size(myShape)-1) == 0)) return                                                  !< empty dataset (last dimension can be empty)
 
-!---------------------------------------------------------------------------------------------------
-! initialize HDF5 data structures
+  myShape = int(shape(dataset),HSIZE_T)
+
   if (present(parallel)) then
     call initialize_read(dset_id, filespace_id, memspace_id, plist_id, aplist_id, &
                          myStart, totalShape, loc_id,myShape,datasetName,parallel)
@@ -1089,6 +1059,7 @@ subroutine HDF5_read_int4(dataset,loc_id,datasetName,parallel)
     call initialize_read(dset_id, filespace_id, memspace_id, plist_id, aplist_id, &
                          myStart, totalShape, loc_id,myShape,datasetName,parallel_default)
   end if
+  if (any(totalShape == 0)) return
 
   call H5Dread_f(dset_id, H5T_NATIVE_INTEGER,dataset,totalShape, hdferr,&
                  file_space_id = filespace_id, xfer_prp = plist_id, mem_space_id = memspace_id)
@@ -1115,13 +1086,9 @@ subroutine HDF5_read_int5(dataset,loc_id,datasetName,parallel)
     totalShape                                                                                      !< shape of the dataset (all processes)
   integer :: hdferr
 
-!---------------------------------------------------------------------------------------------------
-! determine shape of dataset
-  myShape = int(shape(dataset),HSIZE_T)
-  if (any(myShape(1:size(myShape)-1) == 0)) return                                                  !< empty dataset (last dimension can be empty)
 
-!---------------------------------------------------------------------------------------------------
-! initialize HDF5 data structures
+  myShape = int(shape(dataset),HSIZE_T)
+
   if (present(parallel)) then
     call initialize_read(dset_id, filespace_id, memspace_id, plist_id, aplist_id, &
                          myStart, totalShape, loc_id,myShape,datasetName,parallel)
@@ -1129,6 +1096,7 @@ subroutine HDF5_read_int5(dataset,loc_id,datasetName,parallel)
     call initialize_read(dset_id, filespace_id, memspace_id, plist_id, aplist_id, &
                          myStart, totalShape, loc_id,myShape,datasetName,parallel_default)
   end if
+  if (any(totalShape == 0)) return
 
   call H5Dread_f(dset_id, H5T_NATIVE_INTEGER,dataset,totalShape, hdferr,&
                  file_space_id = filespace_id, xfer_prp = plist_id, mem_space_id = memspace_id)
@@ -1155,13 +1123,10 @@ subroutine HDF5_read_int6(dataset,loc_id,datasetName,parallel)
     totalShape                                                                                      !< shape of the dataset (all processes)
   integer :: hdferr
 
-!---------------------------------------------------------------------------------------------------
-! determine shape of dataset
+
   myShape = int(shape(dataset),HSIZE_T)
   if (any(myShape(1:size(myShape)-1) == 0)) return                                                  !< empty dataset (last dimension can be empty)
 
-!---------------------------------------------------------------------------------------------------
-! initialize HDF5 data structures
   if (present(parallel)) then
     call initialize_read(dset_id, filespace_id, memspace_id, plist_id, aplist_id, &
                          myStart, totalShape, loc_id,myShape,datasetName,parallel)
@@ -1169,6 +1134,7 @@ subroutine HDF5_read_int6(dataset,loc_id,datasetName,parallel)
     call initialize_read(dset_id, filespace_id, memspace_id, plist_id, aplist_id, &
                          myStart, totalShape, loc_id,myShape,datasetName,parallel_default)
   end if
+  if (any(totalShape == 0)) return
 
   call H5Dread_f(dset_id, H5T_NATIVE_INTEGER,dataset,totalShape, hdferr,&
                  file_space_id = filespace_id, xfer_prp = plist_id, mem_space_id = memspace_id)
@@ -1195,13 +1161,10 @@ subroutine HDF5_read_int7(dataset,loc_id,datasetName,parallel)
     totalShape                                                                                      !< shape of the dataset (all processes)
   integer :: hdferr
 
-!---------------------------------------------------------------------------------------------------
-! determine shape of dataset
+
   myShape = int(shape(dataset),HSIZE_T)
   if (any(myShape(1:size(myShape)-1) == 0)) return                                                  !< empty dataset (last dimension can be empty)
 
-!---------------------------------------------------------------------------------------------------
-! initialize HDF5 data structures
   if (present(parallel)) then
     call initialize_read(dset_id, filespace_id, memspace_id, plist_id, aplist_id, &
                          myStart, totalShape, loc_id,myShape,datasetName,parallel)
@@ -1209,6 +1172,7 @@ subroutine HDF5_read_int7(dataset,loc_id,datasetName,parallel)
     call initialize_read(dset_id, filespace_id, memspace_id, plist_id, aplist_id, &
                          myStart, totalShape, loc_id,myShape,datasetName,parallel_default)
   end if
+  if (any(totalShape == 0)) return
 
   call H5Dread_f(dset_id, H5T_NATIVE_INTEGER,dataset,totalShape, hdferr,&
                  file_space_id = filespace_id, xfer_prp = plist_id, mem_space_id = memspace_id)
@@ -2031,6 +1995,12 @@ subroutine initialize_read(dset_id, filespace_id, memspace_id, plist_id, aplist_
   myStart(ubound(myStart))  = int(sum(readSize(1:worldrank)),HSIZE_T)
   globalShape = [localShape(1:ubound(localShape,1)-1),int(sum(readSize),HSIZE_T)]
 
+  if (any(globalShape == 0)) then
+    call H5Pclose_f(plist_id, hdferr)
+    if(hdferr < 0) error stop 'HDF5 error'
+    return
+  end if
+
 !--------------------------------------------------------------------------------------------------
 ! create dataspace in memory (local shape)
   call H5Screate_simple_f(size(localShape), localShape, memspace_id, hdferr, localShape)
diff --git a/src/grid/grid_thermal_spectral.f90 b/src/grid/grid_thermal_spectral.f90
index ff0c5b403..6192420df 100644
--- a/src/grid/grid_thermal_spectral.f90
+++ b/src/grid/grid_thermal_spectral.f90
@@ -50,8 +50,8 @@ module grid_thermal_spectral
   real(pReal), dimension(:,:,:), allocatable :: &
     T_current, &                                                                                    !< field of current temperature
     T_lastInc, &                                                                                    !< field of previous temperature
-    T_stagInc                                                                                       !< field of staggered temperature
-
+    T_stagInc, &                                                                                    !< field of staggered temperature
+    dotT_lastInc
 !--------------------------------------------------------------------------------------------------
 ! reference diffusion tensor, mobility etc.
   integer                     :: totalIter = 0                                                      !< total iteration in current increment
@@ -110,6 +110,7 @@ subroutine grid_thermal_spectral_init()
   T_current = discretization_grid_getInitialCondition('T')
   T_lastInc = T_current
   T_stagInc = T_current
+  dotT_lastInc = 0.0_pReal * T_current
 
 !--------------------------------------------------------------------------------------------------
 ! initialize solver specific parts of PETSc
@@ -152,6 +153,7 @@ subroutine grid_thermal_spectral_init()
 
     call HDF5_read(T_current,groupHandle,'T',.false.)
     call HDF5_read(T_lastInc,groupHandle,'T_lastInc',.false.)
+    call HDF5_read(dotT_lastInc,groupHandle,'dotT_lastInc',.false.)
   end if restartRead
 
   ce = 0
@@ -261,9 +263,10 @@ subroutine grid_thermal_spectral_forward(cutBack)
     ce = 0
     do k = 1, cells3;  do j = 1, cells(2);  do i = 1,cells(1)
       ce = ce + 1
-      call homogenization_thermal_setField(T_current(i,j,k),(T_current(i,j,k)-T_lastInc(i,j,k))/params%Delta_t,ce)
+      call homogenization_thermal_setField(T_current(i,j,k),dotT_lastInc(i,j,k),ce)
     end do; end do; end do
   else
+    dotT_lastInc = (T_current - T_lastInc)/params%Delta_t
     T_lastInc = T_current
     call updateReference
   end if
@@ -292,6 +295,7 @@ subroutine grid_thermal_spectral_restartWrite
   groupHandle = HDF5_openGroup(fileHandle,'solver')
   call HDF5_write(T,groupHandle,'T')
   call HDF5_write(T_lastInc,groupHandle,'T_lastInc')
+  call HDF5_write(dotT_lastInc,groupHandle,'dotT_lastInc')
   call HDF5_closeGroup(groupHandle)
   call HDF5_closeFile(fileHandle)
 
diff --git a/src/mesh/FEM_utilities.f90 b/src/mesh/FEM_utilities.f90
index c97c2bbe2..c880bd4cf 100644
--- a/src/mesh/FEM_utilities.f90
+++ b/src/mesh/FEM_utilities.f90
@@ -192,7 +192,8 @@ subroutine utilities_projectBCValues(localVec,section,field,comp,bcPointsIS,BCVa
   call ISGetSize(bcPointsIS,nBcPoints,err_PETSc)
   CHKERRQ(err_PETSc)
   if (nBcPoints > 0) call ISGetIndicesF90(bcPointsIS,bcPoints,err_PETSc)
-  call VecGetArrayF90(localVec,localArray,err_PETSc); CHKERRQ(err_PETSc)
+  call VecGetArrayF90(localVec,localArray,err_PETSc)
+  CHKERRQ(err_PETSc)
   do point = 1, nBcPoints
     call PetscSectionGetFieldDof(section,bcPoints(point),field,numDof,err_PETSc)
     CHKERRQ(err_PETSc)
@@ -202,9 +203,12 @@ subroutine utilities_projectBCValues(localVec,section,field,comp,bcPointsIS,BCVa
       localArray(dof) = localArray(dof) + BCValue + BCDotValue*timeinc
     end do
   end do
-  call VecRestoreArrayF90(localVec,localArray,err_PETSc); CHKERRQ(err_PETSc)
-  call VecAssemblyBegin(localVec, err_PETSc); CHKERRQ(err_PETSc)
-  call VecAssemblyEnd  (localVec, err_PETSc); CHKERRQ(err_PETSc)
+  call VecRestoreArrayF90(localVec,localArray,err_PETSc)
+  CHKERRQ(err_PETSc)
+  call VecAssemblyBegin(localVec, err_PETSc)
+  CHKERRQ(err_PETSc)
+  call VecAssemblyEnd  (localVec, err_PETSc)
+  CHKERRQ(err_PETSc)
   if (nBcPoints > 0) call ISRestoreIndicesF90(bcPointsIS,bcPoints,err_PETSc)
   CHKERRQ(err_PETSc)
 
diff --git a/src/mesh/discretization_mesh.f90 b/src/mesh/discretization_mesh.f90
index a53c50655..d8151b0a4 100644
--- a/src/mesh/discretization_mesh.f90
+++ b/src/mesh/discretization_mesh.f90
@@ -158,7 +158,8 @@ subroutine discretization_mesh_init(restart)
   call MPI_Bcast(mesh_boundaries,mesh_Nboundaries,MPI_INTEGER,0_MPI_INTEGER_KIND,MPI_COMM_WORLD,err_MPI)
   if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error'
 
-  call DMDestroy(globalMesh,err_PETSc); CHKERRQ(err_PETSc)
+  call DMDestroy(globalMesh,err_PETSc)
+  CHKERRQ(err_PETSc)
 
   call DMGetStratumSize(geomMesh,'depth',dimPlex,Nelems,err_PETSc)
   CHKERRQ(err_PETSc)
diff --git a/src/mesh/mesh_mech_FEM.f90 b/src/mesh/mesh_mech_FEM.f90
index 3b162d97b..ef20ed350 100644
--- a/src/mesh/mesh_mech_FEM.f90
+++ b/src/mesh/mesh_mech_FEM.f90
@@ -146,8 +146,10 @@ subroutine FEM_mechanical_init(fieldBC)
 
 !--------------------------------------------------------------------------------------------------
 ! Setup FEM mech mesh
-  call DMClone(geomMesh,mechanical_mesh,err_PETSc); CHKERRQ(err_PETSc)
-  call DMGetDimension(mechanical_mesh,dimPlex,err_PETSc); CHKERRQ(err_PETSc)
+  call DMClone(geomMesh,mechanical_mesh,err_PETSc)
+  CHKERRQ(err_PETSc)
+  call DMGetDimension(mechanical_mesh,dimPlex,err_PETSc)
+  CHKERRQ(err_PETSc)
 
 !--------------------------------------------------------------------------------------------------
 ! Setup FEM mech discretization
@@ -162,24 +164,34 @@ subroutine FEM_mechanical_init(fieldBC)
   call PetscQuadratureSetData(mechQuad,dimPlex,nc,int(nQuadrature,pPETSCINT),qPointsP,qWeightsP,err_PETSc)
   CHKERRQ(err_PETSc)
   call PetscFECreateDefault(PETSC_COMM_SELF,dimPlex,nc,PETSC_TRUE,prefix, &
-                            num%p_i,mechFE,err_PETSc); CHKERRQ(err_PETSc)
-  call PetscFESetQuadrature(mechFE,mechQuad,err_PETSc); CHKERRQ(err_PETSc)
-  call PetscFEGetDimension(mechFE,nBasis,err_PETSc); CHKERRQ(err_PETSc)
+                            num%p_i,mechFE,err_PETSc)
+  CHKERRQ(err_PETSc)
+  call PetscFESetQuadrature(mechFE,mechQuad,err_PETSc)
+  CHKERRQ(err_PETSc)
+  call PetscFEGetDimension(mechFE,nBasis,err_PETSc)
+  CHKERRQ(err_PETSc)
   nBasis = nBasis/nc
   call DMAddField(mechanical_mesh,PETSC_NULL_DMLABEL,mechFE,err_PETSc)
   CHKERRQ(err_PETSc)
-  call DMCreateDS(mechanical_mesh,err_PETSc); CHKERRQ(err_PETSc)
-  call DMGetDS(mechanical_mesh,mechDS,err_PETSc); CHKERRQ(err_PETSc)
-  call PetscDSGetTotalDimension(mechDS,cellDof,err_PETSc); CHKERRQ(err_PETSc)
-  call PetscFEDestroy(mechFE,err_PETSc); CHKERRQ(err_PETSc)
-  call PetscQuadratureDestroy(mechQuad,err_PETSc); CHKERRQ(err_PETSc)
+  call DMCreateDS(mechanical_mesh,err_PETSc)
+  CHKERRQ(err_PETSc)
+  call DMGetDS(mechanical_mesh,mechDS,err_PETSc)
+  CHKERRQ(err_PETSc)
+  call PetscDSGetTotalDimension(mechDS,cellDof,err_PETSc)
+  CHKERRQ(err_PETSc)
+  call PetscFEDestroy(mechFE,err_PETSc)
+  CHKERRQ(err_PETSc)
+  call PetscQuadratureDestroy(mechQuad,err_PETSc)
+  CHKERRQ(err_PETSc)
 
 !--------------------------------------------------------------------------------------------------
 ! Setup FEM mech boundary conditions
   call DMGetLabel(mechanical_mesh,'Face Sets',BCLabel,err_PETSc)
   CHKERRQ(err_PETSc)
-  call DMPlexLabelComplete(mechanical_mesh,BCLabel,err_PETSc); CHKERRQ(err_PETSc)
-  call DMGetLocalSection(mechanical_mesh,section,err_PETSc); CHKERRQ(err_PETSc)
+  call DMPlexLabelComplete(mechanical_mesh,BCLabel,err_PETSc)
+  CHKERRQ(err_PETSc)
+  call DMGetLocalSection(mechanical_mesh,section,err_PETSc)
+  CHKERRQ(err_PETSc)
   allocate(pnumComp(1), source=dimPlex)
   allocate(pnumDof(0:dimPlex), source = 0_pPETSCINT)
   do topologDim = 0, dimPlex
@@ -206,11 +218,14 @@ subroutine FEM_mechanical_init(fieldBC)
       if (bcSize > 0) then
         call DMGetStratumIS(mechanical_mesh,'Face Sets',mesh_boundaries(faceSet),bcPoint,err_PETSc)
         CHKERRQ(err_PETSc)
-        call ISGetIndicesF90(bcPoint,pBcPoint,err_PETSc); CHKERRQ(err_PETSc)
+        call ISGetIndicesF90(bcPoint,pBcPoint,err_PETSc)
+        CHKERRQ(err_PETSc)
         call ISCreateGeneral(PETSC_COMM_WORLD,bcSize,pBcPoint,PETSC_COPY_VALUES,pbcPoints(numBC),err_PETSc)
         CHKERRQ(err_PETSc)
-        call ISRestoreIndicesF90(bcPoint,pBcPoint,err_PETSc); CHKERRQ(err_PETSc)
-        call ISDestroy(bcPoint,err_PETSc); CHKERRQ(err_PETSc)
+        call ISRestoreIndicesF90(bcPoint,pBcPoint,err_PETSc)
+        CHKERRQ(err_PETSc)
+        call ISDestroy(bcPoint,err_PETSc)
+        CHKERRQ(err_PETSc)
       else
         call ISCreateGeneral(PETSC_COMM_WORLD,0_pPETSCINT,[0_pPETSCINT],PETSC_COPY_VALUES,pbcPoints(numBC),err_PETSc)
         CHKERRQ(err_PETSc)
@@ -220,14 +235,17 @@ subroutine FEM_mechanical_init(fieldBC)
   call DMPlexCreateSection(mechanical_mesh,nolabel,pNumComp,pNumDof, &
                            numBC,pBcField,pBcComps,pBcPoints,PETSC_NULL_IS,section,err_PETSc)
   CHKERRQ(err_PETSc)
-  call DMSetSection(mechanical_mesh,section,err_PETSc); CHKERRQ(err_PETSc)
+  call DMSetSection(mechanical_mesh,section,err_PETSc)
+  CHKERRQ(err_PETSc)
   do faceSet = 1, numBC
-    call ISDestroy(pbcPoints(faceSet),err_PETSc); CHKERRQ(err_PETSc)
+    call ISDestroy(pbcPoints(faceSet),err_PETSc)
+    CHKERRQ(err_PETSc)
   enddo
 
 !--------------------------------------------------------------------------------------------------
 ! initialize solver specific parts of PETSc
-  call SNESCreate(PETSC_COMM_WORLD,mechanical_snes,err_PETSc);CHKERRQ(err_PETSc)
+  call SNESCreate(PETSC_COMM_WORLD,mechanical_snes,err_PETSc)
+  CHKERRQ(err_PETSc)
   call SNESSetOptionsPrefix(mechanical_snes,'mechanical_',err_PETSc)
   CHKERRQ(err_PETSc)
   call SNESSetDM(mechanical_snes,mechanical_mesh,err_PETSc)                                         ! set the mesh for non-linear solver
@@ -248,12 +266,15 @@ subroutine FEM_mechanical_init(fieldBC)
   CHKERRQ(err_PETSc)
   call SNESSetTolerances(mechanical_snes,1.0_pReal,0.0_pReal,0.0_pReal,num%itmax,num%itmax,err_PETSc)
   CHKERRQ(err_PETSc)
-  call SNESSetFromOptions(mechanical_snes,err_PETSc); CHKERRQ(err_PETSc)
+  call SNESSetFromOptions(mechanical_snes,err_PETSc)
+  CHKERRQ(err_PETSc)
 
 !--------------------------------------------------------------------------------------------------
 ! init fields
-  call VecSet(solution     ,0.0_pReal,err_PETSc); CHKERRQ(err_PETSc)
-  call VecSet(solution_rate,0.0_pReal,err_PETSc); CHKERRQ(err_PETSc)
+  call VecSet(solution     ,0.0_pReal,err_PETSc)
+  CHKERRQ(err_PETSc)
+  call VecSet(solution_rate,0.0_pReal,err_PETSc)
+  CHKERRQ(err_PETSc)
   allocate(x_scal(cellDof))
   allocate(nodalWeightsP(1))
   allocate(nodalPointsP(dimPlex))
@@ -263,7 +284,8 @@ subroutine FEM_mechanical_init(fieldBC)
   allocate(cellJMat(dimPlex,dimPlex))
   call PetscDSGetDiscretization(mechDS,0_pPETSCINT,mechFE,err_PETSc)
   CHKERRQ(err_PETSc)
-  call PetscFEGetDualSpace(mechFE,mechDualSpace,err_PETSc); CHKERRQ(err_PETSc)
+  call PetscFEGetDualSpace(mechFE,mechDualSpace,err_PETSc)
+  CHKERRQ(err_PETSc)
   call DMPlexGetHeightStratum(mechanical_mesh,0_pPETSCINT,cellStart,cellEnd,err_PETSc)
   CHKERRQ(err_PETSc)
   do cell = cellStart, cellEnd-1                                                                    !< loop over all elements
@@ -365,8 +387,10 @@ subroutine FEM_mechanical_formResidual(dm_local,xx_local,f_local,dummy,err_PETSc
   allocate(pinvcellJ(dimPlex**2))
   allocate(x_scal(cellDof))
 
-  call DMGetLocalSection(dm_local,section,err_PETSc); CHKERRQ(err_PETSc)
-  call DMGetDS(dm_local,prob,err_PETSc); CHKERRQ(err_PETSc)
+  call DMGetLocalSection(dm_local,section,err_PETSc)
+  CHKERRQ(err_PETSc)
+  call DMGetDS(dm_local,prob,err_PETSc)
+  CHKERRQ(err_PETSc)
   call PetscDSGetTabulation(prob,0_pPETSCINT,basisField,basisFieldDer,err_PETSc)
   CHKERRQ(err_PETSc)
   call DMPlexGetHeightStratum(dm_local,0_pPETSCINT,cellStart,cellEnd,err_PETSc)
@@ -383,7 +407,8 @@ subroutine FEM_mechanical_formResidual(dm_local,xx_local,f_local,dummy,err_PETSc
         CHKERRQ(err_PETSc)
         call utilities_projectBCValues(x_local,section,0_pPETSCINT,field-1,bcPoints, &
                                        0.0_pReal,params%fieldBC%componentBC(field)%Value(face),params%timeinc)
-        call ISDestroy(bcPoints,err_PETSc); CHKERRQ(err_PETSc)
+        call ISDestroy(bcPoints,err_PETSc)
+        CHKERRQ(err_PETSc)
       endif
     endif
   enddo; enddo
@@ -464,7 +489,8 @@ subroutine FEM_mechanical_formResidual(dm_local,xx_local,f_local,dummy,err_PETSc
     call DMPlexVecRestoreClosure(dm_local,section,x_local,cell,x_scal,err_PETSc)
     CHKERRQ(err_PETSc)
   enddo
-  call DMRestoreLocalVector(dm_local,x_local,err_PETSc); CHKERRQ(err_PETSc)
+  call DMRestoreLocalVector(dm_local,x_local,err_PETSc)
+  CHKERRQ(err_PETSc)
 
 end subroutine FEM_mechanical_formResidual
 
@@ -531,7 +557,8 @@ subroutine FEM_mechanical_formJacobian(dm_local,xx_local,Jac_pre,Jac,dummy,err_P
         CHKERRQ(err_PETSc)
         call utilities_projectBCValues(x_local,section,0_pPETSCINT,field-1,bcPoints, &
                                        0.0_pReal,params%fieldBC%componentBC(field)%Value(face),params%timeinc)
-        call ISDestroy(bcPoints,err_PETSc); CHKERRQ(err_PETSc)
+        call ISDestroy(bcPoints,err_PETSc)
+        CHKERRQ(err_PETSc)
       endif
     endif
   enddo; enddo
@@ -598,23 +625,32 @@ subroutine FEM_mechanical_formJacobian(dm_local,xx_local,Jac_pre,Jac,dummy,err_P
     call DMPlexVecRestoreClosure(dm_local,section,x_local,cell,x_scal,err_PETSc)
     CHKERRQ(err_PETSc)
   enddo
-  call MatAssemblyBegin(Jac,MAT_FINAL_ASSEMBLY,err_PETSc); CHKERRQ(err_PETSc)
-  call MatAssemblyEnd(Jac,MAT_FINAL_ASSEMBLY,err_PETSc); CHKERRQ(err_PETSc)
-  call MatAssemblyBegin(Jac_pre,MAT_FINAL_ASSEMBLY,err_PETSc); CHKERRQ(err_PETSc)
-  call MatAssemblyEnd(Jac_pre,MAT_FINAL_ASSEMBLY,err_PETSc); CHKERRQ(err_PETSc)
-  call DMRestoreLocalVector(dm_local,x_local,err_PETSc); CHKERRQ(err_PETSc)
+  call MatAssemblyBegin(Jac,MAT_FINAL_ASSEMBLY,err_PETSc)
+  CHKERRQ(err_PETSc)
+  call MatAssemblyEnd(Jac,MAT_FINAL_ASSEMBLY,err_PETSc)
+  CHKERRQ(err_PETSc)
+  call MatAssemblyBegin(Jac_pre,MAT_FINAL_ASSEMBLY,err_PETSc)
+  CHKERRQ(err_PETSc)
+  call MatAssemblyEnd(Jac_pre,MAT_FINAL_ASSEMBLY,err_PETSc)
+  CHKERRQ(err_PETSc)
+  call DMRestoreLocalVector(dm_local,x_local,err_PETSc)
+  CHKERRQ(err_PETSc)
 
 !--------------------------------------------------------------------------------------------------
 ! apply boundary conditions
 #if (PETSC_VERSION_MINOR < 14)
-  call DMPlexCreateRigidBody(dm_local,matnull,err_PETSc); CHKERRQ(err_PETSc)
+  call DMPlexCreateRigidBody(dm_local,matnull,err_PETSc)
+  CHKERRQ(err_PETSc)
 #else
   call DMPlexCreateRigidBody(dm_local,0_pPETSCINT,matnull,err_PETSc)
   CHKERRQ(err_PETSc)
 #endif
-  call MatSetNullSpace(Jac,matnull,err_PETSc); CHKERRQ(err_PETSc)
-  call MatSetNearNullSpace(Jac,matnull,err_PETSc); CHKERRQ(err_PETSc)
-  call MatNullSpaceDestroy(matnull,err_PETSc); CHKERRQ(err_PETSc)
+  call MatSetNullSpace(Jac,matnull,err_PETSc)
+  CHKERRQ(err_PETSc)
+  call MatSetNearNullSpace(Jac,matnull,err_PETSc)
+  CHKERRQ(err_PETSc)
+  call MatNullSpaceDestroy(matnull,err_PETSc)
+  CHKERRQ(err_PETSc)
 
 end subroutine FEM_mechanical_formJacobian
 
@@ -644,15 +680,20 @@ subroutine FEM_mechanical_forward(guess,timeinc,timeinc_old,fieldBC)
   if (guess .and. .not. cutBack) then
     ForwardData = .True.
     homogenization_F0 = homogenization_F
-    call SNESGetDM(mechanical_snes,dm_local,err_PETSc); CHKERRQ(err_PETSc)                               !< retrieve mesh info from mechanical_snes into dm_local
-    call DMGetSection(dm_local,section,err_PETSc); CHKERRQ(err_PETSc)
-    call DMGetLocalVector(dm_local,x_local,err_PETSc); CHKERRQ(err_PETSc)
-    call VecSet(x_local,0.0_pReal,err_PETSc); CHKERRQ(err_PETSc)
+    call SNESGetDM(mechanical_snes,dm_local,err_PETSc)                                              !< retrieve mesh info from mechanical_snes into dm_local
+    CHKERRQ(err_PETSc)
+    call DMGetSection(dm_local,section,err_PETSc)
+    CHKERRQ(err_PETSc)
+    call DMGetLocalVector(dm_local,x_local,err_PETSc)
+    CHKERRQ(err_PETSc)
+    call VecSet(x_local,0.0_pReal,err_PETSc)
+    CHKERRQ(err_PETSc)
     call DMGlobalToLocalBegin(dm_local,solution,INSERT_VALUES,x_local,err_PETSc)                         !< retrieve my partition of global solution vector
     CHKERRQ(err_PETSc)
     call DMGlobalToLocalEnd(dm_local,solution,INSERT_VALUES,x_local,err_PETSc)
     CHKERRQ(err_PETSc)
-    call VecAXPY(solution_local,1.0_pReal,x_local,err_PETSc); CHKERRQ(err_PETSc)
+    call VecAXPY(solution_local,1.0_pReal,x_local,err_PETSc)
+    CHKERRQ(err_PETSc)
     do field = 1, dimPlex; do face = 1, mesh_Nboundaries
       if (fieldBC%componentBC(field)%Mask(face)) then
         call DMGetStratumSize(dm_local,'Face Sets',mesh_boundaries(face),bcSize,err_PETSc)
@@ -661,19 +702,25 @@ subroutine FEM_mechanical_forward(guess,timeinc,timeinc_old,fieldBC)
           CHKERRQ(err_PETSc)
           call utilities_projectBCValues(solution_local,section,0_pPETSCINT,field-1,bcPoints, &
                                          0.0_pReal,fieldBC%componentBC(field)%Value(face),timeinc_old)
-          call ISDestroy(bcPoints,err_PETSc); CHKERRQ(err_PETSc)
+          call ISDestroy(bcPoints,err_PETSc)
+          CHKERRQ(err_PETSc)
         endif
       endif
     enddo; enddo
-    call DMRestoreLocalVector(dm_local,x_local,err_PETSc); CHKERRQ(err_PETSc)
+    call DMRestoreLocalVector(dm_local,x_local,err_PETSc)
+    CHKERRQ(err_PETSc)
 
 !--------------------------------------------------------------------------------------------------
 ! update rate and forward last inc
-    call VecCopy(solution,solution_rate,err_PETSc); CHKERRQ(err_PETSc)
-    call VecScale(solution_rate,timeinc_old**(-1),err_PETSc); CHKERRQ(err_PETSc)
+    call VecCopy(solution,solution_rate,err_PETSc)
+    CHKERRQ(err_PETSc)
+    call VecScale(solution_rate,timeinc_old**(-1),err_PETSc)
+    CHKERRQ(err_PETSc)
   endif
-  call VecCopy(solution_rate,solution,err_PETSc); CHKERRQ(err_PETSc)
-  call VecScale(solution,timeinc,err_PETSc); CHKERRQ(err_PETSc)
+  call VecCopy(solution_rate,solution,err_PETSc)
+  CHKERRQ(err_PETSc)
+  call VecScale(solution,timeinc,err_PETSc)
+  CHKERRQ(err_PETSc)
 
 end subroutine FEM_mechanical_forward
 
@@ -732,24 +779,32 @@ subroutine FEM_mechanical_updateCoords()
     nodeCoords_linear                                                                               !< nodal coordinates (dimPlex*Nnodes)
   PetscScalar, dimension(:), pointer :: x_scal
 
-  call SNESGetDM(mechanical_snes,dm_local,err_PETSc); CHKERRQ(err_PETSc)
-  call DMGetDS(dm_local,mechQuad,err_PETSc); CHKERRQ(err_PETSc)
-  call DMGetLocalSection(dm_local,section,err_PETSc); CHKERRQ(err_PETSc)
-  call DMGetLocalVector(dm_local,x_local,err_PETSc); CHKERRQ(err_PETSc)
-  call DMGetDimension(dm_local,dimPlex,err_PETSc); CHKERRQ(err_PETSc)
+  call SNESGetDM(mechanical_snes,dm_local,err_PETSc)
+  CHKERRQ(err_PETSc)
+  call DMGetDS(dm_local,mechQuad,err_PETSc)
+  CHKERRQ(err_PETSc)
+  call DMGetLocalSection(dm_local,section,err_PETSc)
+  CHKERRQ(err_PETSc)
+  call DMGetLocalVector(dm_local,x_local,err_PETSc)
+  CHKERRQ(err_PETSc)
+  call DMGetDimension(dm_local,dimPlex,err_PETSc)
+  CHKERRQ(err_PETSc)
 
   ! write cell vertex displacements
   call DMPlexGetDepthStratum(dm_local,0_pPETSCINT,pStart,pEnd,err_PETSc)
   CHKERRQ(err_PETSc)
   allocate(nodeCoords(3,pStart:pEnd-1),source=0.0_pReal)
-  call VecGetArrayF90(x_local,nodeCoords_linear,err_PETSc); CHKERRQ(err_PETSc)
+  call VecGetArrayF90(x_local,nodeCoords_linear,err_PETSc)
+  CHKERRQ(err_PETSc)
   do p=pStart, pEnd-1
-    call DMPlexGetPointLocal(dm_local, p, s, e, err_PETSc); CHKERRQ(err_PETSc)
+    call DMPlexGetPointLocal(dm_local, p, s, e, err_PETSc)
+    CHKERRQ(err_PETSc)
     nodeCoords(1:dimPlex,p)=nodeCoords_linear(s+1:e)
   enddo
 
   call discretization_setNodeCoords(nodeCoords)
-  call VecRestoreArrayF90(x_local,nodeCoords_linear,err_PETSc); CHKERRQ(err_PETSc)
+  call VecRestoreArrayF90(x_local,nodeCoords_linear,err_PETSc)
+  CHKERRQ(err_PETSc)
 
   ! write ip displacements
   call DMPlexGetHeightStratum(dm_local,0_pPETSCINT,cellStart,cellEnd,err_PETSc)
@@ -779,7 +834,8 @@ subroutine FEM_mechanical_updateCoords()
     CHKERRQ(err_PETSc)
   end do
   call discretization_setIPcoords(reshape(ipCoords,[3,mesh_NcpElems*nQuadrature]))
-  call DMRestoreLocalVector(dm_local,x_local,err_PETSc); CHKERRQ(err_PETSc)
+  call DMRestoreLocalVector(dm_local,x_local,err_PETSc)
+  CHKERRQ(err_PETSc)
 
 end subroutine FEM_mechanical_updateCoords
 
diff --git a/src/phase_mechanical_plastic_phenopowerlaw.f90 b/src/phase_mechanical_plastic_phenopowerlaw.f90
index f0dc04869..325e4cd64 100644
--- a/src/phase_mechanical_plastic_phenopowerlaw.f90
+++ b/src/phase_mechanical_plastic_phenopowerlaw.f90
@@ -347,8 +347,7 @@ module function phenopowerlaw_dotState(Mp,ph,en) result(dotState)
     sumF
   real(pReal), dimension(param(ph)%sum_N_sl) :: &
     dot_gamma_sl_pos,dot_gamma_sl_neg, &
-    right_SlipSlip
-
+    left_SlipSlip
 
   associate(prm => param(ph), stt => state(ph), &
             dot_xi_sl => dotState(indexDotState(ph)%xi_sl(1):indexDotState(ph)%xi_sl(2)), &
@@ -356,17 +355,17 @@ module function phenopowerlaw_dotState(Mp,ph,en) result(dotState)
             dot_gamma_sl => dotState(indexDotState(ph)%gamma_sl(1):indexDotState(ph)%gamma_sl(2)), &
             dot_gamma_tw => dotState(indexDotState(ph)%gamma_tw(1):indexDotState(ph)%gamma_tw(2)))
 
-    call kinetics_sl(Mp,ph,en,dot_gamma_sl_pos,dot_gamma_sl_neg)
+    call kinetics_sl(Mp,ph,en, dot_gamma_sl_pos,dot_gamma_sl_neg)
     dot_gamma_sl = abs(dot_gamma_sl_pos+dot_gamma_sl_neg)
-    call kinetics_tw(Mp,ph,en,dot_gamma_tw)
-
+    call kinetics_tw(Mp,ph,en, dot_gamma_tw)
     sumF = sum(stt%gamma_tw(:,en)/prm%gamma_char)
-    xi_sl_sat_offset = prm%f_sat_sl_tw*sqrt(sumF)
-    right_SlipSlip = sign(abs(1.0_pReal-stt%xi_sl(:,en) / (prm%xi_inf_sl+xi_sl_sat_offset))**prm%a_sl, &
-                          1.0_pReal-stt%xi_sl(:,en) / (prm%xi_inf_sl+xi_sl_sat_offset))
 
-    dot_xi_sl = prm%h_0_sl_sl * (1.0_pReal + prm%c_1*sumF** prm%c_2) * (1.0_pReal + prm%h_int) &
-                * matmul(prm%h_sl_sl,dot_gamma_sl*right_SlipSlip) &
+    xi_sl_sat_offset = prm%f_sat_sl_tw*sqrt(sumF)
+    left_SlipSlip = sign(abs(1.0_pReal-stt%xi_sl(:,en) / (prm%xi_inf_sl+xi_sl_sat_offset))**prm%a_sl, &
+                             1.0_pReal-stt%xi_sl(:,en) / (prm%xi_inf_sl+xi_sl_sat_offset))
+
+    dot_xi_sl = prm%h_0_sl_sl * (1.0_pReal + prm%c_1 * sumF**prm%c_2) * (1.0_pReal + prm%h_int) &
+                * left_SlipSlip * matmul(prm%h_sl_sl,dot_gamma_sl) &
               + matmul(prm%h_sl_tw,dot_gamma_tw)
 
     dot_xi_tw = prm%h_0_tw_sl * sum(stt%gamma_sl(:,en))**prm%c_3 &