diff --git a/src/constitutive.f90 b/src/constitutive.f90
index 58442283a..f76b94579 100644
--- a/src/constitutive.f90
+++ b/src/constitutive.f90
@@ -47,7 +47,7 @@ module constitutive
     !
     crystallite_S0, &                                                                               !< 2nd Piola-Kirchhoff stress vector at start of FE inc
     crystallite_partitionedS0                                                                       !< 2nd Piola-Kirchhoff stress vector at start of homog inc
-  real(pReal),               dimension(:,:,:,:,:),    allocatable, public, protected :: &
+  real(pReal),               dimension(:,:,:,:,:),    allocatable, public :: &
     crystallite_P, &                                                                                !< 1st Piola-Kirchhoff stress per grain
     crystallite_Lp, &                                                                               !< current plastic velocitiy grad (end of converged time step)
     crystallite_S, &                                                                                !< current 2nd Piola-Kirchhoff stress vector (end of converged time step)
diff --git a/src/constitutive_mech.f90 b/src/constitutive_mech.f90
index dea5ed647..07e48b537 100644
--- a/src/constitutive_mech.f90
+++ b/src/constitutive_mech.f90
@@ -922,7 +922,7 @@ end function integrateStress
 !> @brief integrate stress, state with adaptive 1st order explicit Euler method
 !> using Fixed Point Iteration to adapt the stepsize
 !--------------------------------------------------------------------------------------------------
-subroutine integrateStateFPI(g,i,e)
+module subroutine integrateStateFPI(g,i,e)
 
   integer, intent(in) :: &
     e, &                                                                                            !< element index in element loop
@@ -1027,7 +1027,7 @@ end subroutine integrateStateFPI
 !--------------------------------------------------------------------------------------------------
 !> @brief integrate state with 1st order explicit Euler method
 !--------------------------------------------------------------------------------------------------
-subroutine integrateStateEuler(g,i,e)
+module subroutine integrateStateEuler(g,i,e)
 
   integer, intent(in) :: &
     e, &                                                                                            !< element index in element loop
@@ -1068,7 +1068,7 @@ end subroutine integrateStateEuler
 !--------------------------------------------------------------------------------------------------
 !> @brief integrate stress, state with 1st order Euler method with adaptive step size
 !--------------------------------------------------------------------------------------------------
-subroutine integrateStateAdaptiveEuler(g,i,e)
+module subroutine integrateStateAdaptiveEuler(g,i,e)
 
   integer, intent(in) :: &
     e, &                                                                                            !< element index in element loop
@@ -1127,7 +1127,7 @@ end subroutine integrateStateAdaptiveEuler
 !---------------------------------------------------------------------------------------------------
 !> @brief Integrate state (including stress integration) with the classic Runge Kutta method
 !---------------------------------------------------------------------------------------------------
-subroutine integrateStateRK4(g,i,e)
+module subroutine integrateStateRK4(g,i,e)
 
   integer, intent(in) :: g,i,e
 
@@ -1150,7 +1150,7 @@ end subroutine integrateStateRK4
 !---------------------------------------------------------------------------------------------------
 !> @brief Integrate state (including stress integration) with the Cash-Carp method
 !---------------------------------------------------------------------------------------------------
-subroutine integrateStateRKCK45(g,i,e)
+module subroutine integrateStateRKCK45(g,i,e)
 
   integer, intent(in) :: g,i,e