diff --git a/src/math.f90 b/src/math.f90
index 07b657cab..a875741b3 100644
--- a/src/math.f90
+++ b/src/math.f90
@@ -486,22 +486,25 @@ end function math_inv33
 !> @details Direct Cramer inversion of matrix A. Also returns determinant
 !  Returns an error if not possible, i.e. if determinant is close to zero
 !--------------------------------------------------------------------------------------------------
-pure subroutine math_invert33(InvA, DetA, error, A)
+pure subroutine math_invert33(InvA,DetA,error, A)
 
   real(pReal), dimension(3,3), intent(out) :: InvA
-  real(pReal),                 intent(out) :: DetA
+  real(pReal),                 intent(out), optional :: DetA
   logical,                     intent(out) :: error
   real(pReal), dimension(3,3), intent(in)  :: A
 
+  real(pReal) :: Det
+
 
   InvA(1,1) =  A(2,2) * A(3,3) - A(2,3) * A(3,2)
   InvA(2,1) = -A(2,1) * A(3,3) + A(2,3) * A(3,1)
   InvA(3,1) =  A(2,1) * A(3,2) - A(2,2) * A(3,1)
 
-  DetA = A(1,1) * InvA(1,1) + A(1,2) * InvA(2,1) + A(1,3) * InvA(3,1)
+  Det = A(1,1) * InvA(1,1) + A(1,2) * InvA(2,1) + A(1,3) * InvA(3,1)
 
-  if (dEq0(DetA)) then
+  if (dEq0(Det)) then
     InvA = 0.0_pReal
+    if (present(DetA)) DetA = 0.0_pReal
     error = .true.
   else
     InvA(1,2) = -A(1,2) * A(3,3) + A(1,3) * A(3,2)
@@ -512,7 +515,8 @@ pure subroutine math_invert33(InvA, DetA, error, A)
     InvA(2,3) = -A(1,1) * A(2,3) + A(1,3) * A(2,1)
     InvA(3,3) =  A(1,1) * A(2,2) - A(1,2) * A(2,1)
 
-    InvA = InvA/DetA
+    InvA = InvA/Det
+    if (present(DetA)) DetA = Det
     error = .false.
   end if
 
diff --git a/src/phase_mechanical.f90 b/src/phase_mechanical.f90
index 3a66bdf95..440e196bc 100644
--- a/src/phase_mechanical.f90
+++ b/src/phase_mechanical.f90
@@ -399,8 +399,7 @@ function integrateStress(F,subFp0,subFi0,Delta_t,ph,en) result(broken)
   real(pReal)                         steplengthLp, &
                                       steplengthLi, &
                                       atol_Lp, &
-                                      atol_Li, &
-                                      devNull
+                                      atol_Li
   integer                             NiterationStressLp, &                                         ! number of stress integrations
                                       NiterationStressLi, &                                         ! number of inner stress integrations
                                       ierr, &                                                       ! error indicator for LAPACK
@@ -417,9 +416,9 @@ function integrateStress(F,subFp0,subFi0,Delta_t,ph,en) result(broken)
   Lpguess = phase_mechanical_Lp(ph)%data(1:3,1:3,en)                                              ! take as first guess
   Liguess = phase_mechanical_Li(ph)%data(1:3,1:3,en)                                              ! take as first guess
 
-  call math_invert33(invFp_current,devNull,error,subFp0)
+  call math_invert33(invFp_current,error=error,A=subFp0)
   if (error) return ! error
-  call math_invert33(invFi_current,devNull,error,subFi0)
+  call math_invert33(invFi_current,error=error,A=subFi0)
   if (error) return ! error
 
   A = matmul(F,invFp_current)                                                                       ! intermediate tensor needed later to calculate dFe_dLp
@@ -541,7 +540,7 @@ function integrateStress(F,subFp0,subFi0,Delta_t,ph,en) result(broken)
   end do LiLoop
 
   invFp_new = matmul(invFp_current,B)
-  call math_invert33(Fp_new,devNull,error,invFp_new)
+  call math_invert33(Fp_new,error=error,A=invFp_new)
   if (error) return ! error
 
   phase_mechanical_P(ph)%data(1:3,1:3,en)  = matmul(matmul(F,invFp_new),matmul(S,transpose(invFp_new)))