diff --git a/src/CPFEM.f90 b/src/CPFEM.f90
index 60dc78d2e..d8f4a0c3c 100644
--- a/src/CPFEM.f90
+++ b/src/CPFEM.f90
@@ -150,15 +150,14 @@ subroutine CPFEM_general(mode, ffn, ffn1, temperature_inp, dt, elFE, ip, cauchyS
                                                       H
 
   integer(pInt)                                       elCP, &                                       ! crystal plasticity element number
-                                                      i, j, k, l, m, n, ph, homog, mySource,ma
+                                                      i, j, k, l, m, n, ph, homog, mySource,ce
 
   real(pReal), parameter ::                          ODD_STRESS    = 1e15_pReal, &                  !< return value for stress if terminallyIll
                                                      ODD_JACOBIAN  = 1e50_pReal                     !< return value for jacobian if terminallyIll
 
 
   elCP = discretization_Marc_FEM2DAMASK_elem(elFE)
-
-  ma = (elCP-1) * discretization_nIPs + ip
+  ce   = discretization_Marc_FEM2DAMASK_cell(ip,elFE)
 
   if (debugCPFEM%basic .and. elCP == debugCPFEM%element .and. ip == debugCPFEM%ip) then
     print'(/,a)', '#############################################'
@@ -183,8 +182,8 @@ subroutine CPFEM_general(mode, ffn, ffn1, temperature_inp, dt, elFE, ip, cauchyS
     !    temperature(material_homogenizationAt(elCP))%p(material_homogenizationMemberAt(ip,elCP)) = &
     !      temperature_inp
     !end select chosenThermal1
-    homogenization_F0(1:3,1:3,ma) = ffn
-    homogenization_F(1:3,1:3,ma) = ffn1
+    homogenization_F0(1:3,1:3,ce) = ffn
+    homogenization_F(1:3,1:3,ce) = ffn1
 
   if (iand(mode, CPFEM_CALCRESULTS) /= 0_pInt) then
 
@@ -209,17 +208,17 @@ subroutine CPFEM_general(mode, ffn, ffn1, temperature_inp, dt, elFE, ip, cauchyS
       else terminalIllness
 
         ! translate from P to sigma
-        Kirchhoff = matmul(homogenization_P(1:3,1:3,ma), transpose(homogenization_F(1:3,1:3,ma)))
-        J_inverse  = 1.0_pReal / math_det33(homogenization_F(1:3,1:3,ma))
+        Kirchhoff = matmul(homogenization_P(1:3,1:3,ce), transpose(homogenization_F(1:3,1:3,ce)))
+        J_inverse  = 1.0_pReal / math_det33(homogenization_F(1:3,1:3,ce))
         CPFEM_cs(1:6,ip,elCP) = math_sym33to6(J_inverse * Kirchhoff,weighted=.false.)
 
         !  translate from dP/dF to dCS/dE
         H = 0.0_pReal
         do i=1,3; do j=1,3; do k=1,3; do l=1,3; do m=1,3; do n=1,3
           H(i,j,k,l) = H(i,j,k,l) &
-                     +  homogenization_F(j,m,ma) * homogenization_F(l,n,ma) &
-                                                 * homogenization_dPdF(i,m,k,n,ma) &
-                     -  math_delta(j,l) * homogenization_F(i,m,ma) * homogenization_P(k,m,ma) &
+                     +  homogenization_F(j,m,ce) * homogenization_F(l,n,ce) &
+                                                 * homogenization_dPdF(i,m,k,n,ce) &
+                     -  math_delta(j,l) * homogenization_F(i,m,ce) * homogenization_P(k,m,ce) &
                      +  0.5_pReal * (  Kirchhoff(j,l)*math_delta(i,k) + Kirchhoff(i,k)*math_delta(j,l) &
                                      + Kirchhoff(j,k)*math_delta(i,l) + Kirchhoff(i,l)*math_delta(j,k))
         enddo; enddo; enddo; enddo; enddo; enddo