diff --git a/Z4cow/src/adm.cxx b/Z4cow/src/adm.cxx
index ec676c96..1c9d0a19 100644
--- a/Z4cow/src/adm.cxx
+++ b/Z4cow/src/adm.cxx
@@ -26,9 +26,6 @@ extern "C" void Z4cow_ADM(CCTK_ARGUMENTS) {
   DECLARE_CCTK_ARGUMENTSX_Z4cow_ADM;
   DECLARE_CCTK_PARAMETERS;
 
-  const array<int, dim> indextype = {0, 0, 0};
-  const GF3D2layout layout2(cctkGH, indextype);
-
   // Input Z4c grid functions
   const GF3D2<const CCTK_REAL> &gf_W = W;
   const smat<GF3D2<const CCTK_REAL>, 3> gf_gamt{gammatxx, gammatxy, gammatxz,
diff --git a/Z4cow/src/constraint.cxx b/Z4cow/src/constraint.cxx
index f6038f23..1b7ca6da 100644
--- a/Z4cow/src/constraint.cxx
+++ b/Z4cow/src/constraint.cxx
@@ -43,13 +43,11 @@ extern "C" void Z4cow_Constraints(CCTK_ARGUMENTS) {
       CCTK_DELTA_SPACE(2),
   };
 
-  const array<int, dim> indextype = {0, 0, 0};
   const array<int, dim> nghostzones = {cctk_nghostzones[0], cctk_nghostzones[1],
                                        cctk_nghostzones[2]};
   vect<int, dim> imin, imax;
   GridDescBase(cctkGH).box_int<0, 0, 0>(nghostzones, imin, imax);
   // suffix 2: with ghost zones, suffix 5: without ghost zones
-  const GF3D2layout layout2(cctkGH, indextype);
   const GF3D5layout layout5(imin, imax);
 
   // Input grid functions
diff --git a/Z4cow/src/enforce.cxx b/Z4cow/src/enforce.cxx
index 37990cbd..7fcfab52 100644
--- a/Z4cow/src/enforce.cxx
+++ b/Z4cow/src/enforce.cxx
@@ -22,9 +22,6 @@ extern "C" void Z4cow_Enforce(CCTK_ARGUMENTS) {
   DECLARE_CCTK_ARGUMENTSX_Z4cow_Enforce;
   DECLARE_CCTK_PARAMETERS;
 
-  const array<int, dim> indextype = {0, 0, 0};
-  const GF3D2layout layout1(cctkGH, indextype);
-
   const GF3D2<CCTK_REAL> &gf_W = W;
   const GF3D2<CCTK_REAL> &gf_alphaG = alphaG;
   const smat<GF3D2<CCTK_REAL>, 3> gf_gammat{gammatxx, gammatxy, gammatxz,
@@ -41,13 +38,12 @@ extern "C" void Z4cow_Enforce(CCTK_ARGUMENTS) {
   grid.loop_int_device<0, 0, 0, vsize>(
       grid.nghostzones, [=] ARITH_DEVICE(const PointDesc &p) ARITH_INLINE {
         const vbool mask = mask_for_loop_tail<vbool>(p.i, p.imax);
-        const GF3D2index index1(layout1, p.I);
 
         // Load
-        const vreal W_old = gf_W(mask, index1);
-        const vreal alphaG_old = gf_alphaG(mask, index1);
-        const smat<vreal, 3> gammat_old = gf_gammat(mask, index1);
-        const smat<vreal, 3> At_old = gf_At(mask, index1);
+        const vreal W_old = gf_W(mask, p.I);
+        const vreal alphaG_old = gf_alphaG(mask, p.I);
+        const smat<vreal, 3> gammat_old = gf_gammat(mask, p.I);
+        const smat<vreal, 3> At_old = gf_At(mask, p.I);
 
         // Enforce floors
         const vreal W = fmax(vreal(W_floor), W_old);
@@ -102,10 +98,10 @@ extern "C" void Z4cow_Enforce(CCTK_ARGUMENTS) {
 #endif
 
         // Store
-        gf_W.store(mask, index1, W);
-        gf_gammat.store(mask, index1, gammat);
-        gf_At.store(mask, index1, At);
-        gf_alphaG.store(mask, index1, alphaG);
+        gf_W.store(mask, p.I, W);
+        gf_gammat.store(mask, p.I, gammat);
+        gf_At.store(mask, p.I, At);
+        gf_alphaG.store(mask, p.I, alphaG);
       });
 #ifdef __CUDACC__
   nvtxRangeEnd(range);
diff --git a/Z4cow/src/initial1.cxx b/Z4cow/src/initial1.cxx
index d07dc0ef..687c1304 100644
--- a/Z4cow/src/initial1.cxx
+++ b/Z4cow/src/initial1.cxx
@@ -19,9 +19,6 @@ using namespace Loop;
 extern "C" void Z4cow_Initial1(CCTK_ARGUMENTS) {
   DECLARE_CCTK_ARGUMENTSX_Z4cow_Initial1;
 
-  const array<int, dim> indextype = {0, 0, 0};
-  const GF3D2layout layout1(cctkGH, indextype);
-
   // Input grid functions
   const smat<GF3D2<const CCTK_REAL>, 3> gf_g1{gxx, gxy, gxz, gyy, gyz, gzz};
   const smat<GF3D2<const CCTK_REAL>, 3> gf_K1{kxx, kxy, kxz, kyy, kyz, kzz};
@@ -48,13 +45,12 @@ extern "C" void Z4cow_Initial1(CCTK_ARGUMENTS) {
   grid.loop_int_device<0, 0, 0, vsize>(
       grid.nghostzones, [=] ARITH_DEVICE(const PointDesc &p) ARITH_INLINE {
         const vbool mask = mask_for_loop_tail<vbool>(p.i, p.imax);
-        const GF3D2index index1(layout1, p.I);
 
         // Load (initilize those masked grid with more reasonable value)
-        const smat<vreal, 3> g = gf_g1(mask, index1, one<smat<int, 3>>()());
-        const smat<vreal, 3> K = gf_K1(mask, index1);
-        const vreal alp = gf_alp1(mask, index1, 1);
-        const vec<vreal, 3> beta = gf_beta1(mask, index1);
+        const smat<vreal, 3> g = gf_g1(mask, p.I, one<smat<int, 3>>()());
+        const smat<vreal, 3> K = gf_K1(mask, p.I);
+        const vreal alp = gf_alp1(mask, p.I, 1);
+        const vec<vreal, 3> beta = gf_beta1(mask, p.I);
 
         // Calculate Z4c variables (all except Gammat)
         const vreal detg = calc_det(g);
@@ -74,13 +70,13 @@ extern "C" void Z4cow_Initial1(CCTK_ARGUMENTS) {
         const vec<vreal, 3> betaG([&](int a) ARITH_INLINE { return beta(a); });
 
         // Store
-        gf_W1.store(mask, index1, W);
-        gf_gammat1.store(mask, index1, gammat);
-        gf_Kh1.store(mask, index1, Kh);
-        gf_At1.store(mask, index1, At);
-        gf_Theta1.store(mask, index1, Theta);
-        gf_alphaG1.store(mask, index1, alphaG);
-        gf_betaG1.store(mask, index1, betaG);
+        gf_W1.store(mask, p.I, W);
+        gf_gammat1.store(mask, p.I, gammat);
+        gf_Kh1.store(mask, p.I, Kh);
+        gf_At1.store(mask, p.I, At);
+        gf_Theta1.store(mask, p.I, Theta);
+        gf_alphaG1.store(mask, p.I, alphaG);
+        gf_betaG1.store(mask, p.I, betaG);
       });
 #ifdef __CUDACC__
   nvtxRangeEnd(range);
diff --git a/Z4cow/src/initial2.cxx b/Z4cow/src/initial2.cxx
index 21607ed1..ce5406de 100644
--- a/Z4cow/src/initial2.cxx
+++ b/Z4cow/src/initial2.cxx
@@ -46,9 +46,6 @@ extern "C" void Z4cow_Initial2(CCTK_ARGUMENTS) {
       CCTK_DELTA_SPACE(2),
   };
 
-  const array<int, dim> indextype = {0, 0, 0};
-  const GF3D2layout layout1(cctkGH, indextype);
-
   // Input grid functions
   const smat<GF3D2<const CCTK_REAL>, 3> gf_gammat1{
       gammatxx, gammatxy, gammatxz, gammatyy, gammatyz, gammatzz};
@@ -87,10 +84,9 @@ extern "C" void Z4cow_Initial2(CCTK_ARGUMENTS) {
   grid.loop_int_device<0, 0, 0, vsize>(
       grid.nghostzones, [=] ARITH_DEVICE(const PointDesc &p) ARITH_INLINE {
         const vbool mask = mask_for_loop_tail<vbool>(p.i, p.imax);
-        const GF3D2index index1(layout1, p.I);
 
         // Load
-        const smat<vreal, 3> gammat = gf_gammat1(mask, index1);
+        const smat<vreal, 3> gammat = gf_gammat1(mask, p.I);
 
         // Calculate Z4c variables (only Gamt)
         const smat<vreal, 3> gammatu = calc_inv(gammat, vreal(1));
@@ -108,7 +104,7 @@ extern "C" void Z4cow_Initial2(CCTK_ARGUMENTS) {
         });
 
         // Store
-        gf_Gamt1.store(mask, index1, Gamt);
+        gf_Gamt1.store(mask, p.I, Gamt);
       });
 #ifdef __CUDACC__
   nvtxRangeEnd(range);
diff --git a/Z4cow/src/rhs.cxx b/Z4cow/src/rhs.cxx
index 34add147..004a2a46 100644
--- a/Z4cow/src/rhs.cxx
+++ b/Z4cow/src/rhs.cxx
@@ -45,13 +45,11 @@ extern "C" void Z4cow_RHS(CCTK_ARGUMENTS) {
       CCTK_DELTA_SPACE(2),
   };
 
-  const array<int, dim> indextype = {0, 0, 0};
   const array<int, dim> nghostzones = {cctk_nghostzones[0], cctk_nghostzones[1],
                                        cctk_nghostzones[2]};
   vect<int, dim> imin, imax;
   GridDescBase(cctkGH).box_int<0, 0, 0>(nghostzones, imin, imax);
   // suffix 2: with ghost zones, suffix 5: without ghost zones
-  const GF3D2layout layout2(cctkGH, indextype);
   const GF3D5layout layout5(imin, imax);
 
   // Input grid functions
diff --git a/Z4cow/wolfram/Z4cow_set_ADM.hxx b/Z4cow/wolfram/Z4cow_set_ADM.hxx
index b736f557..6abab4ad 100644
--- a/Z4cow/wolfram/Z4cow_set_ADM.hxx
+++ b/Z4cow/wolfram/Z4cow_set_ADM.hxx
@@ -24,15 +24,14 @@ const GF3D2<CCTK_REAL> &local_ADMbeta3 = gf_ADMbeta(2);
 grid.loop_all_device<0, 0, 0, vsize>(
   grid.nghostzones, [=] ARITH_DEVICE(const PointDesc &p) ARITH_INLINE {
   const vbool mask = mask_for_loop_tail<vbool>(p.i, p.imax);
-  const GF3D2index index2(layout2, p.I);
 
-const auto &tmp_W = gf_W(mask, index2);
-const auto &tmp_gamt = gf_gamt(mask, index2);
-const auto &tmp_exKh = gf_exKh(mask, index2);
-const auto &tmp_exAt = gf_exAt(mask, index2);
-const auto &tmp_Theta = gf_Theta(mask, index2);
-const auto &tmp_alpha = gf_alpha(mask, index2);
-const auto &tmp_beta = gf_beta(mask, index2);
+const auto &tmp_W = gf_W(mask, p.I);
+const auto &tmp_gamt = gf_gamt(mask, p.I);
+const auto &tmp_exKh = gf_exKh(mask, p.I);
+const auto &tmp_exAt = gf_exAt(mask, p.I);
+const auto &tmp_Theta = gf_Theta(mask, p.I);
+const auto &tmp_alpha = gf_alpha(mask, p.I);
+const auto &tmp_beta = gf_beta(mask, p.I);
 
 const vreal W = tmp_W;
 const vreal gamt11 = tmp_gamt(0,0);
@@ -54,67 +53,67 @@ const vreal beta1 = tmp_beta(0);
 const vreal beta2 = tmp_beta(1);
 const vreal beta3 = tmp_beta(2);
 
-local_ADMgam11.store(mask, index2, 
+local_ADMgam11.store(mask, p.I, 
 gamt11/Power(W,2)
 );
 
-local_ADMgam12.store(mask, index2, 
+local_ADMgam12.store(mask, p.I, 
 gamt12/Power(W,2)
 );
 
-local_ADMgam13.store(mask, index2, 
+local_ADMgam13.store(mask, p.I, 
 gamt13/Power(W,2)
 );
 
-local_ADMgam22.store(mask, index2, 
+local_ADMgam22.store(mask, p.I, 
 gamt22/Power(W,2)
 );
 
-local_ADMgam23.store(mask, index2, 
+local_ADMgam23.store(mask, p.I, 
 gamt23/Power(W,2)
 );
 
-local_ADMgam33.store(mask, index2, 
+local_ADMgam33.store(mask, p.I, 
 gamt33/Power(W,2)
 );
 
-local_ADMK11.store(mask, index2, 
+local_ADMK11.store(mask, p.I, 
 (exAt11 + (gamt11*(exKh + 2*Theta))/3.)/Power(W,2)
 );
 
-local_ADMK12.store(mask, index2, 
+local_ADMK12.store(mask, p.I, 
 (exAt12 + (gamt12*(exKh + 2*Theta))/3.)/Power(W,2)
 );
 
-local_ADMK13.store(mask, index2, 
+local_ADMK13.store(mask, p.I, 
 (exAt13 + (gamt13*(exKh + 2*Theta))/3.)/Power(W,2)
 );
 
-local_ADMK22.store(mask, index2, 
+local_ADMK22.store(mask, p.I, 
 (exAt22 + (gamt22*(exKh + 2*Theta))/3.)/Power(W,2)
 );
 
-local_ADMK23.store(mask, index2, 
+local_ADMK23.store(mask, p.I, 
 (exAt23 + (gamt23*(exKh + 2*Theta))/3.)/Power(W,2)
 );
 
-local_ADMK33.store(mask, index2, 
+local_ADMK33.store(mask, p.I, 
 (exAt33 + (gamt33*(exKh + 2*Theta))/3.)/Power(W,2)
 );
 
-local_ADMalpha.store(mask, index2, 
+local_ADMalpha.store(mask, p.I, 
 alpha
 );
 
-local_ADMbeta1.store(mask, index2, 
+local_ADMbeta1.store(mask, p.I, 
 beta1
 );
 
-local_ADMbeta2.store(mask, index2, 
+local_ADMbeta2.store(mask, p.I, 
 beta2
 );
 
-local_ADMbeta3.store(mask, index2, 
+local_ADMbeta3.store(mask, p.I, 
 beta3
 );
 
diff --git a/Z4cow/wolfram/Z4cow_set_ADM.wl b/Z4cow/wolfram/Z4cow_set_ADM.wl
index b1b183a4..eb8ee8cf 100644
--- a/Z4cow/wolfram/Z4cow_set_ADM.wl
+++ b/Z4cow/wolfram/Z4cow_set_ADM.wl
@@ -37,10 +37,9 @@ $MainPrint[] :=
     pr["grid.loop_all_device<0, 0, 0, vsize>("];
     pr["  grid.nghostzones, [=] ARITH_DEVICE(const PointDesc &p) ARITH_INLINE {"];
     pr["  const vbool mask = mask_for_loop_tail<vbool>(p.i, p.imax);"];
-    pr["  const GF3D2index index2(layout2, p.I);"];
     pr[];
 
-    PrintListInitializations[Drop[EvolVarlist, {5}], "gf_", "index2"];
+    PrintListInitializations[Drop[EvolVarlist, {5}], "gf_", "p.I"];
     pr[];
 
     PrintInitializations[{Mode -> "MainIn"}, Drop[EvolVarlist, {5}]];
diff --git a/Z4cow/wolfram/Z4cow_set_constraint.hxx b/Z4cow/wolfram/Z4cow_set_constraint.hxx
index 3b01441a..cbc575a8 100644
--- a/Z4cow/wolfram/Z4cow_set_constraint.hxx
+++ b/Z4cow/wolfram/Z4cow_set_constraint.hxx
@@ -16,12 +16,11 @@ noinline([&]() __attribute__((__flatten__, __hot__)) {
   grid.loop_int_device<0, 0, 0, vsize>(
     grid.nghostzones, [=] ARITH_DEVICE(const PointDesc &p) ARITH_INLINE {
     const vbool mask = mask_for_loop_tail<vbool>(p.i, p.imax);
-    const GF3D2index index2(layout2, p.I);
     const GF3D5index index5(layout5, p.I);
 
-const auto &tmp_eTtt = gf_eTtt(mask, index2);
-const auto &tmp_eTt = gf_eTt(mask, index2);
-const auto &tmp_eT = gf_eT(mask, index2);
+const auto &tmp_eTtt = gf_eTtt(mask, p.I);
+const auto &tmp_eTt = gf_eTt(mask, p.I);
+const auto &tmp_eT = gf_eT(mask, p.I);
 const auto &tmp_W = tl_W(mask, index5);
 const auto &tmp_gamt = tl_gamt(mask, index5);
 const auto &tmp_exKh = tl_exKh(mask, index5);
@@ -980,25 +979,25 @@ vreal trdexAtUU3
 ;
 
 
-local_ZtC1.store(mask, index2, 
+local_ZtC1.store(mask, p.I, 
 (trGt1 - trGtd1)/2.
 );
 
-local_ZtC2.store(mask, index2, 
+local_ZtC2.store(mask, p.I, 
 (trGt2 - trGtd2)/2.
 );
 
-local_ZtC3.store(mask, index2, 
+local_ZtC3.store(mask, p.I, 
 (trGt3 - trGtd3)/2.
 );
 
-local_HC.store(mask, index2, 
+local_HC.store(mask, p.I, 
 exAt11*exAtUU11 + 2*exAt12*exAtUU12 + 2*exAt13*exAtUU13 + exAt22*exAtUU22 + 
   2*exAt23*exAtUU23 + exAt33*exAtUU33 - (2*Power(exKh,2))/3. - 16*cpi*rho - 
   (8*exKh*Theta)/3. - (8*Power(Theta,2))/3. + trR
 );
 
-local_MtC1.store(mask, index2, 
+local_MtC1.store(mask, p.I, 
 (-4*dlnW1*exAtUU11 - 4*dlnW2*exAtUU12 - 4*dlnW3*exAtUU13 + 
     3*exAtUU11*Gt111 + 6*exAtUU12*Gt112 + 6*exAtUU13*Gt113 + 
     3*exAtUU22*Gt122 + 6*exAtUU23*Gt123 + 3*exAtUU33*Gt133 - 
@@ -1008,7 +1007,7 @@ local_MtC1.store(mask, index2,
     3*trdexAtUU1)/3.
 );
 
-local_MtC2.store(mask, index2, 
+local_MtC2.store(mask, p.I, 
 (-4*dlnW1*exAtUU12 - 4*dlnW2*exAtUU22 - 4*dlnW3*exAtUU23 + 
     3*exAtUU11*Gt211 + 6*exAtUU12*Gt212 + 6*exAtUU13*Gt213 + 
     3*exAtUU22*Gt222 + 6*exAtUU23*Gt223 + 3*exAtUU33*Gt233 - 
@@ -1018,7 +1017,7 @@ local_MtC2.store(mask, index2,
     3*trdexAtUU2)/3.
 );
 
-local_MtC3.store(mask, index2, 
+local_MtC3.store(mask, p.I, 
 (-4*dlnW1*exAtUU13 - 4*dlnW2*exAtUU23 - 4*dlnW3*exAtUU33 + 
     3*exAtUU11*Gt311 + 6*exAtUU12*Gt312 + 6*exAtUU13*Gt313 + 
     3*exAtUU22*Gt322 + 6*exAtUU23*Gt323 + 3*exAtUU33*Gt333 - 
diff --git a/Z4cow/wolfram/Z4cow_set_constraint.wl b/Z4cow/wolfram/Z4cow_set_constraint.wl
index 7ebd97aa..e76c28a5 100644
--- a/Z4cow/wolfram/Z4cow_set_constraint.wl
+++ b/Z4cow/wolfram/Z4cow_set_constraint.wl
@@ -43,11 +43,10 @@ $MainPrint[] :=
     pr["  grid.loop_int_device<0, 0, 0, vsize>("];
     pr["    grid.nghostzones, [=] ARITH_DEVICE(const PointDesc &p) ARITH_INLINE {"];
     pr["    const vbool mask = mask_for_loop_tail<vbool>(p.i, p.imax);"];
-    pr["    const GF3D2index index2(layout2, p.I);"];
     pr["    const GF3D5index index5(layout5, p.I);"];
     pr[];
 
-    PrintListInitializations[TmunuVarlist, "gf_", "index2"];
+    PrintListInitializations[TmunuVarlist, "gf_", "p.I"];
     PrintListInitializations[EvolVarlist, "tl_", "index5"];
     PrintListInitializations[Drop[dEvolVarlist, -2], "tl_", "index5"];
     PrintListInitializations[Drop[ddEvolVarlist, -2], "tl_", "index5"];
diff --git a/Z4cow/wolfram/Z4cow_set_rhs.hxx b/Z4cow/wolfram/Z4cow_set_rhs.hxx
index fcbf5d52..4aa15dfe 100644
--- a/Z4cow/wolfram/Z4cow_set_rhs.hxx
+++ b/Z4cow/wolfram/Z4cow_set_rhs.hxx
@@ -31,15 +31,14 @@ noinline([&]() __attribute__((__flatten__, __hot__)) {
   grid.loop_int_device<0, 0, 0, vsize>(
     grid.nghostzones, [=] ARITH_DEVICE(const PointDesc &p) ARITH_INLINE {
     const vbool mask = mask_for_loop_tail<vbool>(p.i, p.imax);
-    const GF3D2index index2(layout2, p.I);
     const GF3D5index index5(layout5, p.I);
 
 const vreal x0 = p.x + Arith::iota<vreal>() * p.dx;
 const vreal ceta = calceta(x0, p.y, p.z);
 
-const auto &tmp_eTtt = gf_eTtt(mask, index2);
-const auto &tmp_eTt = gf_eTt(mask, index2);
-const auto &tmp_eT = gf_eT(mask, index2);
+const auto &tmp_eTtt = gf_eTtt(mask, p.I);
+const auto &tmp_eTt = gf_eTt(mask, p.I);
+const auto &tmp_eT = gf_eT(mask, p.I);
 const auto &tmp_W = tl_W(mask, index5);
 const auto &tmp_gamt = tl_gamt(mask, index5);
 const auto &tmp_exKh = tl_exKh(mask, index5);
@@ -1290,51 +1289,51 @@ invgam11*Ss11 + 2*invgam12*Ss12 + 2*invgam13*Ss13 + invgam22*Ss22 +
 ;
 
 
-local_dtW.store(mask, index2, 
+local_dtW.store(mask, p.I, 
 beta1*dW1 + beta2*dW2 + beta3*dW3 - ((dbeta11 + dbeta22 + dbeta33)*W)/3. + 
   (alpha*exKh*W)/3. + (2*alpha*Theta*W)/3.
 );
 
-local_dtgamt11.store(mask, index2, 
+local_dtgamt11.store(mask, p.I, 
 beta1*dgamt111 + beta2*dgamt211 + beta3*dgamt311 - 2*alpha*exAt11 + 
   2*dbeta11*gamt11 - (2*(dbeta11 + dbeta22 + dbeta33)*gamt11)/3. + 
   2*dbeta12*gamt12 + 2*dbeta13*gamt13
 );
 
-local_dtgamt12.store(mask, index2, 
+local_dtgamt12.store(mask, p.I, 
 beta1*dgamt112 + beta2*dgamt212 + beta3*dgamt312 - 2*alpha*exAt12 + 
   dbeta21*gamt11 + dbeta11*gamt12 + dbeta22*gamt12 - 
   (2*(dbeta11 + dbeta22 + dbeta33)*gamt12)/3. + dbeta23*gamt13 + 
   dbeta12*gamt22 + dbeta13*gamt23
 );
 
-local_dtgamt13.store(mask, index2, 
+local_dtgamt13.store(mask, p.I, 
 beta1*dgamt113 + beta2*dgamt213 + beta3*dgamt313 - 2*alpha*exAt13 + 
   dbeta31*gamt11 + dbeta32*gamt12 + dbeta11*gamt13 + dbeta33*gamt13 - 
   (2*(dbeta11 + dbeta22 + dbeta33)*gamt13)/3. + dbeta12*gamt23 + 
   dbeta13*gamt33
 );
 
-local_dtgamt22.store(mask, index2, 
+local_dtgamt22.store(mask, p.I, 
 beta1*dgamt122 + beta2*dgamt222 + beta3*dgamt322 - 2*alpha*exAt22 + 
   2*dbeta21*gamt12 + 2*dbeta22*gamt22 - 
   (2*(dbeta11 + dbeta22 + dbeta33)*gamt22)/3. + 2*dbeta23*gamt23
 );
 
-local_dtgamt23.store(mask, index2, 
+local_dtgamt23.store(mask, p.I, 
 beta1*dgamt123 + beta2*dgamt223 + beta3*dgamt323 - 2*alpha*exAt23 + 
   dbeta31*gamt12 + dbeta21*gamt13 + dbeta32*gamt22 + dbeta22*gamt23 + 
   dbeta33*gamt23 - (2*(dbeta11 + dbeta22 + dbeta33)*gamt23)/3. + 
   dbeta23*gamt33
 );
 
-local_dtgamt33.store(mask, index2, 
+local_dtgamt33.store(mask, p.I, 
 beta1*dgamt133 + beta2*dgamt233 + beta3*dgamt333 - 2*alpha*exAt33 + 
   2*dbeta31*gamt13 + 2*dbeta32*gamt23 + 2*dbeta33*gamt33 - 
   (2*(dbeta11 + dbeta22 + dbeta33)*gamt33)/3.
 );
 
-local_dtexKh.store(mask, index2, 
+local_dtexKh.store(mask, p.I, 
 beta1*dexKh1 + beta2*dexKh2 + beta3*dexKh3 + alpha*exAt11*exAtUU11 + 
   2*alpha*exAt12*exAtUU12 + 2*alpha*exAt13*exAtUU13 + 
   alpha*exAt22*exAtUU22 + 2*alpha*exAt23*exAtUU23 + alpha*exAt33*exAtUU33 + 
@@ -1345,7 +1344,7 @@ beta1*dexKh1 + beta2*dexKh2 + beta3*dexKh3 + alpha*exAt11*exAtUU11 +
   (4*alpha*Power(Theta,2))/3. + 4*alpha*cpi*trSs
 );
 
-local_dtexAt11.store(mask, index2, 
+local_dtexAt11.store(mask, p.I, 
 (3*beta1*dexAt111 + 3*beta2*dexAt211 + 3*beta3*dexAt311 + 
     4*dbeta11*exAt11 - 2*dbeta22*exAt11 - 2*dbeta33*exAt11 + 
     6*dbeta12*exAt12 + 6*dbeta13*exAt13 + 3*alpha*exAt11*exKh - 
@@ -1373,7 +1372,7 @@ local_dtexAt11.store(mask, index2,
     8*alpha*cpi*gam11*invgam33*Ss33*Power(W,2))/3.
 );
 
-local_dtexAt12.store(mask, index2, 
+local_dtexAt12.store(mask, p.I, 
 (3*beta1*dexAt112 + 3*beta2*dexAt212 + 3*beta3*dexAt312 + 
     3*dbeta21*exAt11 + dbeta11*exAt12 + dbeta22*exAt12 - 2*dbeta33*exAt12 + 
     3*dbeta23*exAt13 + 3*dbeta12*exAt22 + 3*dbeta13*exAt23 + 
@@ -1404,7 +1403,7 @@ local_dtexAt12.store(mask, index2,
     8*alpha*cpi*gam12*invgam33*Ss33*Power(W,2))/3.
 );
 
-local_dtexAt13.store(mask, index2, 
+local_dtexAt13.store(mask, p.I, 
 (3*beta1*dexAt113 + 3*beta2*dexAt213 + 3*beta3*dexAt313 + 
     3*dbeta31*exAt11 + 3*dbeta32*exAt12 + dbeta11*exAt13 - 
     2*dbeta22*exAt13 + dbeta33*exAt13 + 3*dbeta12*exAt23 + 
@@ -1435,7 +1434,7 @@ local_dtexAt13.store(mask, index2,
     8*alpha*cpi*gam13*invgam33*Ss33*Power(W,2))/3.
 );
 
-local_dtexAt22.store(mask, index2, 
+local_dtexAt22.store(mask, p.I, 
 (3*beta1*dexAt122 + 3*beta2*dexAt222 + 3*beta3*dexAt322 + 
     6*dbeta21*exAt12 - 2*dbeta11*exAt22 + 4*dbeta22*exAt22 - 
     2*dbeta33*exAt22 + 6*dbeta23*exAt23 + 3*alpha*exAt22*exKh - 
@@ -1464,7 +1463,7 @@ local_dtexAt22.store(mask, index2,
     8*alpha*cpi*gam22*invgam33*Ss33*Power(W,2))/3.
 );
 
-local_dtexAt23.store(mask, index2, 
+local_dtexAt23.store(mask, p.I, 
 (3*beta1*dexAt123 + 3*beta2*dexAt223 + 3*beta3*dexAt323 + 
     3*dbeta31*exAt12 + 3*dbeta21*exAt13 + 3*dbeta32*exAt22 - 
     2*dbeta11*exAt23 + dbeta22*exAt23 + dbeta33*exAt23 + 3*dbeta23*exAt33 + 
@@ -1495,7 +1494,7 @@ local_dtexAt23.store(mask, index2,
     8*alpha*cpi*gam23*invgam33*Ss33*Power(W,2))/3.
 );
 
-local_dtexAt33.store(mask, index2, 
+local_dtexAt33.store(mask, p.I, 
 (3*beta1*dexAt133 + 3*beta2*dexAt233 + 3*beta3*dexAt333 + 
     6*dbeta31*exAt13 + 6*dbeta32*exAt23 - 2*dbeta11*exAt33 - 
     2*dbeta22*exAt33 + 4*dbeta33*exAt33 + 3*alpha*exAt33*exKh - 
@@ -1524,7 +1523,7 @@ local_dtexAt33.store(mask, index2,
     8*alpha*cpi*gam33*invgam33*Ss33*Power(W,2))/3.
 );
 
-local_dttrGt1.store(mask, index2, 
+local_dttrGt1.store(mask, p.I, 
 (3*beta1*dtrGt11 + 3*beta2*dtrGt21 + 3*beta3*dtrGt31 - 6*dalpha1*exAtUU11 - 
     18*alpha*dlnW1*exAtUU11 - 6*dalpha2*exAtUU12 - 
     18*alpha*dlnW2*exAtUU12 - 6*dalpha3*exAtUU13 - 
@@ -1545,7 +1544,7 @@ local_dttrGt1.store(mask, index2,
     2*dbeta33*trGtd1 - 3*dbeta21*trGtd2 - 3*dbeta31*trGtd3)/3.
 );
 
-local_dttrGt2.store(mask, index2, 
+local_dttrGt2.store(mask, p.I, 
 (3*beta1*dtrGt12 + 3*beta2*dtrGt22 + 3*beta3*dtrGt32 - 6*dalpha1*exAtUU12 - 
     18*alpha*dlnW1*exAtUU12 - 6*dalpha2*exAtUU22 - 
     18*alpha*dlnW2*exAtUU22 - 6*dalpha3*exAtUU23 - 
@@ -1567,7 +1566,7 @@ local_dttrGt2.store(mask, index2,
   3.
 );
 
-local_dttrGt3.store(mask, index2, 
+local_dttrGt3.store(mask, p.I, 
 (3*beta1*dtrGt13 + 3*beta2*dtrGt23 + 3*beta3*dtrGt33 - 6*dalpha1*exAtUU13 - 
     18*alpha*dlnW1*exAtUU13 - 6*dalpha2*exAtUU23 - 
     18*alpha*dlnW2*exAtUU23 - 6*dalpha3*exAtUU33 - 
@@ -1588,7 +1587,7 @@ local_dttrGt3.store(mask, index2,
     2*dbeta22*trGtd3 - dbeta33*trGtd3)/3.
 );
 
-local_dtTheta.store(mask, index2, 
+local_dtTheta.store(mask, p.I, 
 (6*beta1*dTheta1 + 6*beta2*dTheta2 + 6*beta3*dTheta3 - 
     3*alpha*exAt11*exAtUU11 - 6*alpha*exAt12*exAtUU12 - 
     6*alpha*exAt13*exAtUU13 - 3*alpha*exAt22*exAtUU22 - 
@@ -1598,19 +1597,19 @@ local_dtTheta.store(mask, index2,
     8*alpha*Power(Theta,2) + 3*alpha*trR)/6.
 );
 
-local_dtalpha.store(mask, index2, 
+local_dtalpha.store(mask, p.I, 
 beta1*dalpha1 + beta2*dalpha2 + beta3*dalpha3 - alpha*cmuL*exKh
 );
 
-local_dtbeta1.store(mask, index2, 
+local_dtbeta1.store(mask, p.I, 
 beta1*(-ceta + dbeta11) + beta2*dbeta21 + beta3*dbeta31 + cmuS*trGt1
 );
 
-local_dtbeta2.store(mask, index2, 
+local_dtbeta2.store(mask, p.I, 
 beta1*dbeta12 + beta2*(-ceta + dbeta22) + beta3*dbeta32 + cmuS*trGt2
 );
 
-local_dtbeta3.store(mask, index2, 
+local_dtbeta3.store(mask, p.I, 
 beta1*dbeta13 + beta2*dbeta23 + beta3*(-ceta + dbeta33) + cmuS*trGt3
 );
 
diff --git a/Z4cow/wolfram/Z4cow_set_rhs.wl b/Z4cow/wolfram/Z4cow_set_rhs.wl
index f672e627..58f8853b 100644
--- a/Z4cow/wolfram/Z4cow_set_rhs.wl
+++ b/Z4cow/wolfram/Z4cow_set_rhs.wl
@@ -43,7 +43,6 @@ $MainPrint[] :=
     pr["  grid.loop_int_device<0, 0, 0, vsize>("];
     pr["    grid.nghostzones, [=] ARITH_DEVICE(const PointDesc &p) ARITH_INLINE {"];
     pr["    const vbool mask = mask_for_loop_tail<vbool>(p.i, p.imax);"];
-    pr["    const GF3D2index index2(layout2, p.I);"];
     pr["    const GF3D5index index5(layout5, p.I);"];
     pr[];
 
@@ -51,7 +50,7 @@ $MainPrint[] :=
     pr["const vreal ceta = calceta(x0, p.y, p.z);"];
     pr[];
 
-    PrintListInitializations[TmunuVarlist, "gf_", "index2"];
+    PrintListInitializations[TmunuVarlist, "gf_", "p.I"];
     PrintListInitializations[EvolVarlist, "tl_", "index5"];
     PrintListInitializations[dEvolVarlist, "tl_", "index5"];
     PrintListInitializations[ddEvolVarlist, "tl_", "index5"];