diff --git a/code/constitutive_dislotwin.f90 b/code/constitutive_dislotwin.f90
index f87e14c34..daefcf26d 100644
--- a/code/constitutive_dislotwin.f90
+++ b/code/constitutive_dislotwin.f90
@@ -542,12 +542,15 @@ do i = 1_pInt,maxNinstance
 
     
    !* Elasticity matrix and shear modulus according to material.config
-   constitutive_dislotwin_Cslip_66(:,:,i) = lattice_symmetrizeC66(constitutive_dislotwin_structureName(i),&
-                                                                      constitutive_dislotwin_Cslip_66) 
+   constitutive_dislotwin_Cslip_66(1:6,1:6,i) = lattice_symmetrizeC66(constitutive_dislotwin_structureName(i),&
+                                                                     constitutive_dislotwin_Cslip_66) 
    constitutive_dislotwin_Gmod(i) = &
-   0.2_pReal*(constitutive_dislotwin_Cslip_66(1,1,i)-constitutive_dislotwin_Cslip_66(1,2,i))+0.3_pReal*constitutive_dislotwin_Cslip_66(4,4,i)
-   constitutive_dislotwin_Cslip_66(:,:,i) = math_Mandel3333to66(math_Voigt66to3333(constitutive_dislotwin_Cslip_66(:,:,i)))
-   constitutive_dislotwin_Cslip_3333(:,:,:,:,i) = math_Voigt66to3333(constitutive_dislotwin_Cslip_66(:,:,i))
+      0.2_pReal*(constitutive_dislotwin_Cslip_66(1,1,i)-constitutive_dislotwin_Cslip_66(1,2,i)) &
+     +0.3_pReal*constitutive_dislotwin_Cslip_66(4,4,i)
+   constitutive_dislotwin_Cslip_66(1:6,1:6,i) = &
+      math_Mandel3333to66(math_Voigt66to3333(constitutive_dislotwin_Cslip_66(1:6,1:6,i)))
+   constitutive_dislotwin_Cslip_3333(1:3,1:3,1:3,1:3,i) = &
+      math_Voigt66to3333(constitutive_dislotwin_Cslip_66(1:6,1:6,i))
 
 
    !* Process slip related parameters ------------------------------------------------
@@ -621,7 +624,8 @@ do i = 1_pInt,maxNinstance
              lattice_Qtwin(o,s,index_otherFamily+j,myStructure)
          enddo ; enddo ; enddo ; enddo
        enddo ; enddo ; enddo ; enddo
-       constitutive_dislotwin_Ctwin_66(:,:,index_myFamily+j,i) = math_Mandel3333to66(constitutive_dislotwin_Ctwin_3333(:,:,:,:,index_myFamily+j,i))
+       constitutive_dislotwin_Ctwin_66(1:6,1:6,index_myFamily+j,i) = &
+         math_Mandel3333to66(constitutive_dislotwin_Ctwin_3333(1:3,1:3,1:3,1:3,index_myFamily+j,i))
 
     !* Interaction matrices
 
diff --git a/code/constitutive_phenopowerlaw.f90 b/code/constitutive_phenopowerlaw.f90
index 717196968..b71331995 100644
--- a/code/constitutive_phenopowerlaw.f90
+++ b/code/constitutive_phenopowerlaw.f90
@@ -680,9 +680,11 @@ subroutine constitutive_phenopowerlaw_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v,Temp
                                            math_Mandel6to33(lattice_Sslip_v(1:6,2*k+1,index_myFamily+i,structID))
      enddo
      gdot_slip_pos(j) = 0.5_pReal*constitutive_phenopowerlaw_gdot0_slip(matID)* &
-                    ((abs(tau_slip_pos(j))/state(ipc,ip,el)%p(j))**constitutive_phenopowerlaw_n_slip(matID))*sign(1.0_pReal,tau_slip_pos(j))
+                    ((abs(tau_slip_pos(j))/state(ipc,ip,el)%p(j))**constitutive_phenopowerlaw_n_slip(matID))*&
+                                                                    sign(1.0_pReal,tau_slip_pos(j))
      gdot_slip_neg(j) = 0.5_pReal*constitutive_phenopowerlaw_gdot0_slip(matID)* &
-                    ((abs(tau_slip_neg(j))/state(ipc,ip,el)%p(j))**constitutive_phenopowerlaw_n_slip(matID))*sign(1.0_pReal,tau_slip_neg(j))
+                    ((abs(tau_slip_neg(j))/state(ipc,ip,el)%p(j))**constitutive_phenopowerlaw_n_slip(matID))*&
+                                                                    sign(1.0_pReal,tau_slip_neg(j))
      Lp = Lp + (1.0_pReal-state(ipc,ip,el)%p(index_F))*&                     ! 1-F
                (gdot_slip_pos(j)+gdot_slip_neg(j))*lattice_Sslip(1:3,1:3,index_myFamily+i,structID)
 
diff --git a/code/constitutive_titanmod.f90 b/code/constitutive_titanmod.f90
index 71791b6fb..25432d067 100644
--- a/code/constitutive_titanmod.f90
+++ b/code/constitutive_titanmod.f90
@@ -864,15 +864,18 @@ do i = 1_pInt,maxNinstance
        endif
    enddo
    
-write(6,*) 'Determining elasticity matrix'
+   write(6,*) 'Determining elasticity matrix'
 
    !* Elasticity matrix and shear modulus according to material.config
    constitutive_titanmod_Cslip_66(:,:,i) = lattice_symmetrizeC66(constitutive_titanmod_structureName(i),&
                                                                       constitutive_titanmod_Cslip_66)  
    constitutive_titanmod_Gmod(i) = &
-   0.2_pReal*(constitutive_titanmod_Cslip_66(1,1,i)-constitutive_titanmod_Cslip_66(1,2,i))+0.3_pReal*constitutive_titanmod_Cslip_66(4,4,i)
-   constitutive_titanmod_Cslip_66(:,:,i) =     math_Mandel3333to66(math_Voigt66to3333(constitutive_titanmod_Cslip_66(:,:,i)))
-   constitutive_titanmod_Cslip_3333(:,:,:,:,i) = math_Voigt66to3333(constitutive_titanmod_Cslip_66(:,:,i))
+        0.2_pReal*(constitutive_titanmod_Cslip_66(1,1,i)-constitutive_titanmod_Cslip_66(1,2,i))&
+      + 0.3_pReal*constitutive_titanmod_Cslip_66(4,4,i)
+   constitutive_titanmod_Cslip_66(1:6,1:6,i) = &
+      math_Mandel3333to66(math_Voigt66to3333(constitutive_titanmod_Cslip_66(1:6,1:6,i)))
+   constitutive_titanmod_Cslip_3333(1:3,1:3,1:3,1:3,i) = &
+      math_Voigt66to3333(constitutive_titanmod_Cslip_66(1:6,1:6,i))
    
    !* Construction of the twin elasticity matrices
    do j=1_pInt,lattice_maxNtwinFamily
@@ -887,7 +890,8 @@ write(6,*) 'Determining elasticity matrix'
              lattice_Qtwin(n,r,sum(lattice_NslipSystem(1:j-1_pInt,myStructure))+k,myStructure)* &
              lattice_Qtwin(o,s,sum(lattice_NslipSystem(1:j-1_pInt,myStructure))+k,myStructure)
            enddo ; enddo ; enddo ; enddo ; enddo ; enddo ; enddo ; enddo
-         constitutive_titanmod_Ctwin_66(:,:,k,i) = math_Mandel3333to66(constitutive_titanmod_Ctwin_3333(:,:,:,:,k,i))
+         constitutive_titanmod_Ctwin_66(1:6,1:6,k,i) =  &
+            math_Mandel3333to66(constitutive_titanmod_Ctwin_3333(1:3,1:3,1:3,1:3,k,i))
         enddo
    enddo
 
diff --git a/code/homogenization_RGC.f90 b/code/homogenization_RGC.f90
index 78380d073..1e0c32c73 100644
--- a/code/homogenization_RGC.f90
+++ b/code/homogenization_RGC.f90
@@ -845,7 +845,7 @@ endsubroutine
 !********************************************************************
 ! derive average stress and stiffness from constituent quantities
 !********************************************************************
-function homogenization_RGC_averageTemperature(&
+pure function homogenization_RGC_averageTemperature(&
    Temperature, &   ! temperature
    ip, &            ! my integration point
    el  &            ! my element
@@ -1302,7 +1302,7 @@ endfunction
 !********************************************************************
 ! subroutine to map grain ID from in 3D (local position) to in 1D (global array)
 !********************************************************************
-function homogenization_RGC_grain3to1(&
+pure function homogenization_RGC_grain3to1(&
    grain3, &       ! grain ID in 3D array (pos.x,pos.y,pos.z)
    homID &         ! homogenization ID
   )
@@ -1313,7 +1313,7 @@ function homogenization_RGC_grain3to1(&
  integer(pInt), dimension (3), intent(in) :: grain3
  integer(pInt)                :: homogenization_RGC_grain3to1
  integer(pInt), dimension (3) :: nGDim
- integer(pInt) homID
+ integer(pInt), intent(in) :: homID
 
 !* Get the grain ID
  nGDim = homogenization_RGC_Ngrains(:,homID)
@@ -1324,7 +1324,7 @@ endfunction
 !********************************************************************
 ! subroutine to map interface ID from 4D (normal and local position) into 1D (global array)
 !********************************************************************
-function homogenization_RGC_interface4to1(&
+pure function homogenization_RGC_interface4to1(&
    iFace4D, &        ! interface ID in 4D array (n.dir,pos.x,pos.y,pos.z)
    homID &           ! homogenization ID
   )
@@ -1335,7 +1335,7 @@ function homogenization_RGC_interface4to1(&
  integer(pInt), dimension (4), intent(in) :: iFace4D
  integer(pInt)                :: homogenization_RGC_interface4to1
  integer(pInt), dimension (3) :: nGDim,nIntFace
- integer(pInt) homID
+ integer(pInt), intent(in) :: homID
 
  nGDim = homogenization_RGC_Ngrains(:,homID)
 !* Compute the total number of interfaces, which ...
@@ -1374,7 +1374,7 @@ function homogenization_RGC_interface1to4(&
  integer(pInt), dimension (4) :: homogenization_RGC_interface1to4
  integer(pInt), intent(in)    :: iFace1D
  integer(pInt), dimension (3) :: nGDim,nIntFace
- integer(pInt) homID
+ integer(pInt), intent(in) :: homID
 
  nGDim = homogenization_RGC_Ngrains(:,homID)
 !* Compute the total number of interfaces, which ...