Regenerate examples (except Advect)

Advect doesn't work at the moment due to a problem with TensorTools.

22 files changed, 1847 insertions, 1442 deletions

diff --git a/Examples/EulerAuto/configuration.ccl b/Examples/EulerAuto/configuration.ccl
index 023aac7..0a66ec2 100644
--- a/Examples/EulerAuto/configuration.ccl
+++ b/Examples/EulerAuto/configuration.ccl
@@ -1,3 +1,6 @@
 # File produced by Kranc
 
 REQUIRES GenericFD
+OPTIONAL LoopControl
+{
+}
diff --git a/Examples/EulerAuto/interface.ccl b/Examples/EulerAuto/interface.ccl
index 409ad9e..0323809 100644
--- a/Examples/EulerAuto/interface.ccl
+++ b/Examples/EulerAuto/interface.ccl
@@ -10,6 +10,7 @@ USES INCLUDE: GenericFD.h
 USES INCLUDE: Symmetry.h
 USES INCLUDE: sbp_calc_coeffs.h
 USES INCLUDE: Boundary.h
+USES INCLUDE: loopcontrol.h
 
 CCTK_INT FUNCTION MoLRegisterEvolved(CCTK_INT IN EvolvedIndex, CCTK_INT IN RHSIndex)
 USES FUNCTION MoLRegisterEvolved
diff --git a/Examples/EulerAuto/param.ccl b/Examples/EulerAuto/param.ccl
index e36bd21..06e8caf 100644
--- a/Examples/EulerAuto/param.ccl
+++ b/Examples/EulerAuto/param.ccl
@@ -8,6 +8,7 @@ shares: GenericFD
 shares: MethodOfLines
 
 USES CCTK_INT MoL_Num_Evolved_Vars
+USES CCTK_INT MoL_Num_ArrayEvolved_Vars
 
 restricted:
 CCTK_INT verbose "verbose" STEERABLE=ALWAYS
@@ -94,6 +95,12 @@ CCTK_INT EulerAuto_MaxNumEvolvedVars "Number of evolved variables used by this t
 } 5
 
 restricted:
+CCTK_INT EulerAuto_MaxNumArrayEvolvedVars "Number of Array evolved variables used by this thorn" ACCUMULATOR-BASE=MethodofLines::MoL_Num_ArrayEvolved_Vars STEERABLE=RECOVER
+{
+  0:0 :: "Number of Array evolved variables used by this thorn"
+} 0
+
+restricted:
 CCTK_INT timelevels "Number of active timelevels" STEERABLE=RECOVER
 {
   0:3 :: ""
@@ -106,6 +113,12 @@ CCTK_INT rhs_timelevels "Number of active RHS timelevels" STEERABLE=RECOVER
 } 1
 
 restricted:
+CCTK_INT other_timelevels "Number of active timelevels for non-evolved grid functions" STEERABLE=RECOVER
+{
+  0:3 :: ""
+} 1
+
+restricted:
 CCTK_INT eulerauto_initial_shock_calc_every "eulerauto_initial_shock_calc_every" STEERABLE=ALWAYS
 {
   *:* :: ""
diff --git a/Examples/EulerAuto/schedule.ccl b/Examples/EulerAuto/schedule.ccl
index 4f5faa0..2911d7b 100644
--- a/Examples/EulerAuto/schedule.ccl
+++ b/Examples/EulerAuto/schedule.ccl
@@ -1,41 +1,95 @@
 # File produced by Kranc
 
 
-STORAGE: Den_flux_group[1]
+if (other_timelevels == 1)
+{
+  STORAGE: Den_flux_group[1]
+}
 
-STORAGE: Den_lr_group[1]
+if (other_timelevels == 1)
+{
+  STORAGE: Den_lr_group[1]
+}
 
-STORAGE: En_flux_group[1]
+if (other_timelevels == 1)
+{
+  STORAGE: En_flux_group[1]
+}
 
-STORAGE: En_lr_group[1]
+if (other_timelevels == 1)
+{
+  STORAGE: En_lr_group[1]
+}
 
-STORAGE: p_group[1]
+if (other_timelevels == 1)
+{
+  STORAGE: p_group[1]
+}
 
-STORAGE: p_lr_group[1]
+if (other_timelevels == 1)
+{
+  STORAGE: p_lr_group[1]
+}
 
-STORAGE: rho_group[1]
+if (other_timelevels == 1)
+{
+  STORAGE: rho_group[1]
+}
 
-STORAGE: rho_lr_group[1]
+if (other_timelevels == 1)
+{
+  STORAGE: rho_lr_group[1]
+}
 
-STORAGE: S1_flux_group[1]
+if (other_timelevels == 1)
+{
+  STORAGE: S1_flux_group[1]
+}
 
-STORAGE: S1_lr_group[1]
+if (other_timelevels == 1)
+{
+  STORAGE: S1_lr_group[1]
+}
 
-STORAGE: S2_flux_group[1]
+if (other_timelevels == 1)
+{
+  STORAGE: S2_flux_group[1]
+}
 
-STORAGE: S2_lr_group[1]
+if (other_timelevels == 1)
+{
+  STORAGE: S2_lr_group[1]
+}
 
-STORAGE: S3_flux_group[1]
+if (other_timelevels == 1)
+{
+  STORAGE: S3_flux_group[1]
+}
 
-STORAGE: S3_lr_group[1]
+if (other_timelevels == 1)
+{
+  STORAGE: S3_lr_group[1]
+}
 
-STORAGE: v1_lr_group[1]
+if (other_timelevels == 1)
+{
+  STORAGE: v1_lr_group[1]
+}
 
-STORAGE: v2_lr_group[1]
+if (other_timelevels == 1)
+{
+  STORAGE: v2_lr_group[1]
+}
 
-STORAGE: v3_lr_group[1]
+if (other_timelevels == 1)
+{
+  STORAGE: v3_lr_group[1]
+}
 
-STORAGE: v_group[1]
+if (other_timelevels == 1)
+{
+  STORAGE: v_group[1]
+}
 
 if (timelevels == 1)
 {
@@ -188,12 +242,7 @@ schedule group eulerauto_cons_calc_reconstruct_1_bc_group in eulerauto_cons_calc
   # no language specified
 } "eulerauto_cons_calc_reconstruct_1"
 
-schedule group eulerauto_cons_calc_reconstruct_1_bc_group at CCTK_POSTRESTRICT
-{
-  # no language specified
-} "eulerauto_cons_calc_reconstruct_1"
-
-schedule group eulerauto_cons_calc_reconstruct_1_bc_group at CCTK_POSTRESTRICTINITIAL
+schedule group eulerauto_cons_calc_reconstruct_1_bc_group in MoL_PseudoEvolutionBoundaries after MoL_PostStep
 {
   # no language specified
 } "eulerauto_cons_calc_reconstruct_1"
@@ -234,12 +283,7 @@ schedule group eulerauto_cons_calc_flux_1_bc_group in eulerauto_cons_calc_flux_1
   # no language specified
 } "eulerauto_cons_calc_flux_1"
 
-schedule group eulerauto_cons_calc_flux_1_bc_group at CCTK_POSTRESTRICT
-{
-  # no language specified
-} "eulerauto_cons_calc_flux_1"
-
-schedule group eulerauto_cons_calc_flux_1_bc_group at CCTK_POSTRESTRICTINITIAL
+schedule group eulerauto_cons_calc_flux_1_bc_group in MoL_PseudoEvolutionBoundaries after MoL_PostStep
 {
   # no language specified
 } "eulerauto_cons_calc_flux_1"
@@ -280,12 +324,7 @@ schedule group eulerauto_cons_calc_reconstruct_2_bc_group in eulerauto_cons_calc
   # no language specified
 } "eulerauto_cons_calc_reconstruct_2"
 
-schedule group eulerauto_cons_calc_reconstruct_2_bc_group at CCTK_POSTRESTRICT
-{
-  # no language specified
-} "eulerauto_cons_calc_reconstruct_2"
-
-schedule group eulerauto_cons_calc_reconstruct_2_bc_group at CCTK_POSTRESTRICTINITIAL
+schedule group eulerauto_cons_calc_reconstruct_2_bc_group in MoL_PseudoEvolutionBoundaries after MoL_PostStep
 {
   # no language specified
 } "eulerauto_cons_calc_reconstruct_2"
@@ -326,12 +365,7 @@ schedule group eulerauto_cons_calc_flux_2_bc_group in eulerauto_cons_calc_flux_2
   # no language specified
 } "eulerauto_cons_calc_flux_2"
 
-schedule group eulerauto_cons_calc_flux_2_bc_group at CCTK_POSTRESTRICT
-{
-  # no language specified
-} "eulerauto_cons_calc_flux_2"
-
-schedule group eulerauto_cons_calc_flux_2_bc_group at CCTK_POSTRESTRICTINITIAL
+schedule group eulerauto_cons_calc_flux_2_bc_group in MoL_PseudoEvolutionBoundaries after MoL_PostStep
 {
   # no language specified
 } "eulerauto_cons_calc_flux_2"
@@ -372,12 +406,7 @@ schedule group eulerauto_cons_calc_reconstruct_3_bc_group in eulerauto_cons_calc
   # no language specified
 } "eulerauto_cons_calc_reconstruct_3"
 
-schedule group eulerauto_cons_calc_reconstruct_3_bc_group at CCTK_POSTRESTRICT
-{
-  # no language specified
-} "eulerauto_cons_calc_reconstruct_3"
-
-schedule group eulerauto_cons_calc_reconstruct_3_bc_group at CCTK_POSTRESTRICTINITIAL
+schedule group eulerauto_cons_calc_reconstruct_3_bc_group in MoL_PseudoEvolutionBoundaries after MoL_PostStep
 {
   # no language specified
 } "eulerauto_cons_calc_reconstruct_3"
@@ -418,12 +447,7 @@ schedule group eulerauto_cons_calc_flux_3_bc_group in eulerauto_cons_calc_flux_3
   # no language specified
 } "eulerauto_cons_calc_flux_3"
 
-schedule group eulerauto_cons_calc_flux_3_bc_group at CCTK_POSTRESTRICT
-{
-  # no language specified
-} "eulerauto_cons_calc_flux_3"
-
-schedule group eulerauto_cons_calc_flux_3_bc_group at CCTK_POSTRESTRICTINITIAL
+schedule group eulerauto_cons_calc_flux_3_bc_group in MoL_PseudoEvolutionBoundaries after MoL_PostStep
 {
   # no language specified
 } "eulerauto_cons_calc_flux_3"
diff --git a/Examples/EulerAuto/src/Differencing.h b/Examples/EulerAuto/src/Differencing.h
index 0908172..913cfa8 100644
--- a/Examples/EulerAuto/src/Differencing.h
+++ b/Examples/EulerAuto/src/Differencing.h
@@ -1,12 +1,144 @@
-#define DiffPlus1(u) (p1o1*(-(u)[di*(0)+dj*(0)+dk*(0)] + (u)[di*(1)+dj*(0)+dk*(0)]))
-#define DiffPlus2(u) (p1o1*(-(u)[di*(0)+dj*(0)+dk*(0)] + (u)[di*(0)+dj*(1)+dk*(0)]))
-#define DiffPlus3(u) (p1o1*(-(u)[di*(0)+dj*(0)+dk*(0)] + (u)[di*(0)+dj*(0)+dk*(1)]))
-#define DiffMinus1(u) (p1o1*((u)[di*(0)+dj*(0)+dk*(0)] - (u)[di*(-1)+dj*(0)+dk*(0)]))
-#define DiffMinus2(u) (p1o1*((u)[di*(0)+dj*(0)+dk*(0)] - (u)[di*(0)+dj*(-1)+dk*(0)]))
-#define DiffMinus3(u) (p1o1*((u)[di*(0)+dj*(0)+dk*(0)] - (u)[di*(0)+dj*(0)+dk*(-1)]))
-#define ShiftMinus1(u) (p1o1*(u)[di*(-1)+dj*(0)+dk*(0)])
-#define ShiftMinus2(u) (p1o1*(u)[di*(0)+dj*(-1)+dk*(0)])
-#define ShiftMinus3(u) (p1o1*(u)[di*(0)+dj*(0)+dk*(-1)])
-#define PDplus1(u) (p1odx*(-(u)[di*(0)+dj*(0)+dk*(0)] + (u)[di*(1)+dj*(0)+dk*(0)]))
-#define PDplus2(u) (p1ody*(-(u)[di*(0)+dj*(0)+dk*(0)] + (u)[di*(0)+dj*(1)+dk*(0)]))
-#define PDplus3(u) (p1odz*(-(u)[di*(0)+dj*(0)+dk*(0)] + (u)[di*(0)+dj*(0)+dk*(1)]))
+#ifndef KRANC_DIFF_FUNCTIONS
+#  define DiffPlus1(u) ((-KRANC_GFOFFSET3D(u,0,0,0) + KRANC_GFOFFSET3D(u,1,0,0))*p1o1)
+#else
+#  define DiffPlus1(u) (DiffPlus1_impl(u,p1o1,cdj,cdk))
+static CCTK_REAL DiffPlus1_impl(CCTK_REAL const* restrict const u, CCTK_REAL const p1o1, ptrdiff_t const cdj, ptrdiff_t const cdk) CCTK_ATTRIBUTE_NOINLINE CCTK_ATTRIBUTE_UNUSED;
+static CCTK_REAL DiffPlus1_impl(CCTK_REAL const* restrict const u, CCTK_REAL const p1o1, ptrdiff_t const cdj, ptrdiff_t const cdk)
+{
+  ptrdiff_t const cdi=sizeof(CCTK_REAL);
+  return (-KRANC_GFOFFSET3D(u,0,0,0) + KRANC_GFOFFSET3D(u,1,0,0))*p1o1;
+}
+#endif
+
+#ifndef KRANC_DIFF_FUNCTIONS
+#  define DiffPlus2(u) ((-KRANC_GFOFFSET3D(u,0,0,0) + KRANC_GFOFFSET3D(u,0,1,0))*p1o1)
+#else
+#  define DiffPlus2(u) (DiffPlus2_impl(u,p1o1,cdj,cdk))
+static CCTK_REAL DiffPlus2_impl(CCTK_REAL const* restrict const u, CCTK_REAL const p1o1, ptrdiff_t const cdj, ptrdiff_t const cdk) CCTK_ATTRIBUTE_NOINLINE CCTK_ATTRIBUTE_UNUSED;
+static CCTK_REAL DiffPlus2_impl(CCTK_REAL const* restrict const u, CCTK_REAL const p1o1, ptrdiff_t const cdj, ptrdiff_t const cdk)
+{
+  ptrdiff_t const cdi=sizeof(CCTK_REAL);
+  return (-KRANC_GFOFFSET3D(u,0,0,0) + KRANC_GFOFFSET3D(u,0,1,0))*p1o1;
+}
+#endif
+
+#ifndef KRANC_DIFF_FUNCTIONS
+#  define DiffPlus3(u) ((-KRANC_GFOFFSET3D(u,0,0,0) + KRANC_GFOFFSET3D(u,0,0,1))*p1o1)
+#else
+#  define DiffPlus3(u) (DiffPlus3_impl(u,p1o1,cdj,cdk))
+static CCTK_REAL DiffPlus3_impl(CCTK_REAL const* restrict const u, CCTK_REAL const p1o1, ptrdiff_t const cdj, ptrdiff_t const cdk) CCTK_ATTRIBUTE_NOINLINE CCTK_ATTRIBUTE_UNUSED;
+static CCTK_REAL DiffPlus3_impl(CCTK_REAL const* restrict const u, CCTK_REAL const p1o1, ptrdiff_t const cdj, ptrdiff_t const cdk)
+{
+  ptrdiff_t const cdi=sizeof(CCTK_REAL);
+  return (-KRANC_GFOFFSET3D(u,0,0,0) + KRANC_GFOFFSET3D(u,0,0,1))*p1o1;
+}
+#endif
+
+#ifndef KRANC_DIFF_FUNCTIONS
+#  define DiffMinus1(u) ((KRANC_GFOFFSET3D(u,0,0,0) - KRANC_GFOFFSET3D(u,-1,0,0))*p1o1)
+#else
+#  define DiffMinus1(u) (DiffMinus1_impl(u,p1o1,cdj,cdk))
+static CCTK_REAL DiffMinus1_impl(CCTK_REAL const* restrict const u, CCTK_REAL const p1o1, ptrdiff_t const cdj, ptrdiff_t const cdk) CCTK_ATTRIBUTE_NOINLINE CCTK_ATTRIBUTE_UNUSED;
+static CCTK_REAL DiffMinus1_impl(CCTK_REAL const* restrict const u, CCTK_REAL const p1o1, ptrdiff_t const cdj, ptrdiff_t const cdk)
+{
+  ptrdiff_t const cdi=sizeof(CCTK_REAL);
+  return (KRANC_GFOFFSET3D(u,0,0,0) - KRANC_GFOFFSET3D(u,-1,0,0))*p1o1;
+}
+#endif
+
+#ifndef KRANC_DIFF_FUNCTIONS
+#  define DiffMinus2(u) ((KRANC_GFOFFSET3D(u,0,0,0) - KRANC_GFOFFSET3D(u,0,-1,0))*p1o1)
+#else
+#  define DiffMinus2(u) (DiffMinus2_impl(u,p1o1,cdj,cdk))
+static CCTK_REAL DiffMinus2_impl(CCTK_REAL const* restrict const u, CCTK_REAL const p1o1, ptrdiff_t const cdj, ptrdiff_t const cdk) CCTK_ATTRIBUTE_NOINLINE CCTK_ATTRIBUTE_UNUSED;
+static CCTK_REAL DiffMinus2_impl(CCTK_REAL const* restrict const u, CCTK_REAL const p1o1, ptrdiff_t const cdj, ptrdiff_t const cdk)
+{
+  ptrdiff_t const cdi=sizeof(CCTK_REAL);
+  return (KRANC_GFOFFSET3D(u,0,0,0) - KRANC_GFOFFSET3D(u,0,-1,0))*p1o1;
+}
+#endif
+
+#ifndef KRANC_DIFF_FUNCTIONS
+#  define DiffMinus3(u) ((KRANC_GFOFFSET3D(u,0,0,0) - KRANC_GFOFFSET3D(u,0,0,-1))*p1o1)
+#else
+#  define DiffMinus3(u) (DiffMinus3_impl(u,p1o1,cdj,cdk))
+static CCTK_REAL DiffMinus3_impl(CCTK_REAL const* restrict const u, CCTK_REAL const p1o1, ptrdiff_t const cdj, ptrdiff_t const cdk) CCTK_ATTRIBUTE_NOINLINE CCTK_ATTRIBUTE_UNUSED;
+static CCTK_REAL DiffMinus3_impl(CCTK_REAL const* restrict const u, CCTK_REAL const p1o1, ptrdiff_t const cdj, ptrdiff_t const cdk)
+{
+  ptrdiff_t const cdi=sizeof(CCTK_REAL);
+  return (KRANC_GFOFFSET3D(u,0,0,0) - KRANC_GFOFFSET3D(u,0,0,-1))*p1o1;
+}
+#endif
+
+#ifndef KRANC_DIFF_FUNCTIONS
+#  define ShiftMinus1(u) (KRANC_GFOFFSET3D(u,-1,0,0)*p1o1)
+#else
+#  define ShiftMinus1(u) (ShiftMinus1_impl(u,p1o1,cdj,cdk))
+static CCTK_REAL ShiftMinus1_impl(CCTK_REAL const* restrict const u, CCTK_REAL const p1o1, ptrdiff_t const cdj, ptrdiff_t const cdk) CCTK_ATTRIBUTE_NOINLINE CCTK_ATTRIBUTE_UNUSED;
+static CCTK_REAL ShiftMinus1_impl(CCTK_REAL const* restrict const u, CCTK_REAL const p1o1, ptrdiff_t const cdj, ptrdiff_t const cdk)
+{
+  ptrdiff_t const cdi=sizeof(CCTK_REAL);
+  return KRANC_GFOFFSET3D(u,-1,0,0)*p1o1;
+}
+#endif
+
+#ifndef KRANC_DIFF_FUNCTIONS
+#  define ShiftMinus2(u) (KRANC_GFOFFSET3D(u,0,-1,0)*p1o1)
+#else
+#  define ShiftMinus2(u) (ShiftMinus2_impl(u,p1o1,cdj,cdk))
+static CCTK_REAL ShiftMinus2_impl(CCTK_REAL const* restrict const u, CCTK_REAL const p1o1, ptrdiff_t const cdj, ptrdiff_t const cdk) CCTK_ATTRIBUTE_NOINLINE CCTK_ATTRIBUTE_UNUSED;
+static CCTK_REAL ShiftMinus2_impl(CCTK_REAL const* restrict const u, CCTK_REAL const p1o1, ptrdiff_t const cdj, ptrdiff_t const cdk)
+{
+  ptrdiff_t const cdi=sizeof(CCTK_REAL);
+  return KRANC_GFOFFSET3D(u,0,-1,0)*p1o1;
+}
+#endif
+
+#ifndef KRANC_DIFF_FUNCTIONS
+#  define ShiftMinus3(u) (KRANC_GFOFFSET3D(u,0,0,-1)*p1o1)
+#else
+#  define ShiftMinus3(u) (ShiftMinus3_impl(u,p1o1,cdj,cdk))
+static CCTK_REAL ShiftMinus3_impl(CCTK_REAL const* restrict const u, CCTK_REAL const p1o1, ptrdiff_t const cdj, ptrdiff_t const cdk) CCTK_ATTRIBUTE_NOINLINE CCTK_ATTRIBUTE_UNUSED;
+static CCTK_REAL ShiftMinus3_impl(CCTK_REAL const* restrict const u, CCTK_REAL const p1o1, ptrdiff_t const cdj, ptrdiff_t const cdk)
+{
+  ptrdiff_t const cdi=sizeof(CCTK_REAL);
+  return KRANC_GFOFFSET3D(u,0,0,-1)*p1o1;
+}
+#endif
+
+#ifndef KRANC_DIFF_FUNCTIONS
+#  define PDplus1(u) ((-KRANC_GFOFFSET3D(u,0,0,0) + KRANC_GFOFFSET3D(u,1,0,0))*p1odx)
+#else
+#  define PDplus1(u) (PDplus1_impl(u,p1odx,cdj,cdk))
+static CCTK_REAL PDplus1_impl(CCTK_REAL const* restrict const u, CCTK_REAL const p1odx, ptrdiff_t const cdj, ptrdiff_t const cdk) CCTK_ATTRIBUTE_NOINLINE CCTK_ATTRIBUTE_UNUSED;
+static CCTK_REAL PDplus1_impl(CCTK_REAL const* restrict const u, CCTK_REAL const p1odx, ptrdiff_t const cdj, ptrdiff_t const cdk)
+{
+  ptrdiff_t const cdi=sizeof(CCTK_REAL);
+  return (-KRANC_GFOFFSET3D(u,0,0,0) + KRANC_GFOFFSET3D(u,1,0,0))*p1odx;
+}
+#endif
+
+#ifndef KRANC_DIFF_FUNCTIONS
+#  define PDplus2(u) ((-KRANC_GFOFFSET3D(u,0,0,0) + KRANC_GFOFFSET3D(u,0,1,0))*p1ody)
+#else
+#  define PDplus2(u) (PDplus2_impl(u,p1ody,cdj,cdk))
+static CCTK_REAL PDplus2_impl(CCTK_REAL const* restrict const u, CCTK_REAL const p1ody, ptrdiff_t const cdj, ptrdiff_t const cdk) CCTK_ATTRIBUTE_NOINLINE CCTK_ATTRIBUTE_UNUSED;
+static CCTK_REAL PDplus2_impl(CCTK_REAL const* restrict const u, CCTK_REAL const p1ody, ptrdiff_t const cdj, ptrdiff_t const cdk)
+{
+  ptrdiff_t const cdi=sizeof(CCTK_REAL);
+  return (-KRANC_GFOFFSET3D(u,0,0,0) + KRANC_GFOFFSET3D(u,0,1,0))*p1ody;
+}
+#endif
+
+#ifndef KRANC_DIFF_FUNCTIONS
+#  define PDplus3(u) ((-KRANC_GFOFFSET3D(u,0,0,0) + KRANC_GFOFFSET3D(u,0,0,1))*p1odz)
+#else
+#  define PDplus3(u) (PDplus3_impl(u,p1odz,cdj,cdk))
+static CCTK_REAL PDplus3_impl(CCTK_REAL const* restrict const u, CCTK_REAL const p1odz, ptrdiff_t const cdj, ptrdiff_t const cdk) CCTK_ATTRIBUTE_NOINLINE CCTK_ATTRIBUTE_UNUSED;
+static CCTK_REAL PDplus3_impl(CCTK_REAL const* restrict const u, CCTK_REAL const p1odz, ptrdiff_t const cdj, ptrdiff_t const cdk)
+{
+  ptrdiff_t const cdi=sizeof(CCTK_REAL);
+  return (-KRANC_GFOFFSET3D(u,0,0,0) + KRANC_GFOFFSET3D(u,0,0,1))*p1odz;
+}
+#endif
+
diff --git a/Examples/EulerAuto/src/RegisterMoL.cc b/Examples/EulerAuto/src/RegisterMoL.cc
index c01d0fd..76554bb 100644
--- a/Examples/EulerAuto/src/RegisterMoL.cc
+++ b/Examples/EulerAuto/src/RegisterMoL.cc
@@ -17,5 +17,7 @@ extern "C" void EulerAuto_RegisterVars(CCTK_ARGUMENTS)
   ierr += MoLRegisterEvolved(CCTK_VarIndex("EulerAuto::S1"),  CCTK_VarIndex("EulerAuto::S1rhs"));
   ierr += MoLRegisterEvolved(CCTK_VarIndex("EulerAuto::S2"),  CCTK_VarIndex("EulerAuto::S2rhs"));
   ierr += MoLRegisterEvolved(CCTK_VarIndex("EulerAuto::S3"),  CCTK_VarIndex("EulerAuto::S3rhs"));
+  
+  /* Register all the evolved Array functions with MoL */
   return;
 }
diff --git a/Examples/EulerAuto/src/eulerauto_cons_calc_conserved.cc b/Examples/EulerAuto/src/eulerauto_cons_calc_conserved.cc
index 9f483da..b41803a 100644
--- a/Examples/EulerAuto/src/eulerauto_cons_calc_conserved.cc
+++ b/Examples/EulerAuto/src/eulerauto_cons_calc_conserved.cc
@@ -12,6 +12,8 @@
 #include "cctk_Parameters.h"
 #include "GenericFD.h"
 #include "Differencing.h"
+#include "cctk_Loop.h"
+#include "loopcontrol.h"
 
 /* Define macros used in calculations */
 #define INITVALUE (42)
@@ -20,42 +22,28 @@
 #define SQR(x) ((x) * (x))
 #define CUB(x) ((x) * (x) * (x))
 
-static void eulerauto_cons_calc_conserved_Body(cGH const * restrict const cctkGH, int const dir, int const face, CCTK_REAL const normal[3], CCTK_REAL const tangentA[3], CCTK_REAL const tangentB[3], int const min[3], int const max[3], int const n_subblock_gfs, CCTK_REAL * restrict const subblock_gfs[])
+static void eulerauto_cons_calc_conserved_Body(cGH const * restrict const cctkGH, int const dir, int const face, CCTK_REAL const normal[3], CCTK_REAL const tangentA[3], CCTK_REAL const tangentB[3], int const imin[3], int const imax[3], int const n_subblock_gfs, CCTK_REAL * restrict const subblock_gfs[])
 {
   DECLARE_CCTK_ARGUMENTS;
   DECLARE_CCTK_PARAMETERS;
   
   
-  /* Declare the variables used for looping over grid points */
-  CCTK_INT i, j, k;
-  // CCTK_INT index = INITVALUE;
-  
   /* Declare finite differencing variables */
   
-  if (verbose > 1)
-  {
-    CCTK_VInfo(CCTK_THORNSTRING,"Entering eulerauto_cons_calc_conserved_Body");
-  }
-  
-  if (cctk_iteration % eulerauto_cons_calc_conserved_calc_every != eulerauto_cons_calc_conserved_calc_offset)
-  {
-    return;
-  }
-  
-  const char *groups[] = {"EulerAuto::Den_group","EulerAuto::En_group","EulerAuto::p_group","EulerAuto::rho_group","EulerAuto::S_group","EulerAuto::v_group"};
-  GenericFD_AssertGroupStorage(cctkGH, "eulerauto_cons_calc_conserved", 6, groups);
-  
-  
   /* Include user-supplied include files */
   
   /* Initialise finite differencing variables */
   ptrdiff_t const di = 1;
   ptrdiff_t const dj = CCTK_GFINDEX3D(cctkGH,0,1,0) - CCTK_GFINDEX3D(cctkGH,0,0,0);
   ptrdiff_t const dk = CCTK_GFINDEX3D(cctkGH,0,0,1) - CCTK_GFINDEX3D(cctkGH,0,0,0);
+  ptrdiff_t const cdi = sizeof(CCTK_REAL) * di;
+  ptrdiff_t const cdj = sizeof(CCTK_REAL) * dj;
+  ptrdiff_t const cdk = sizeof(CCTK_REAL) * dk;
   CCTK_REAL const dx = ToReal(CCTK_DELTA_SPACE(0));
   CCTK_REAL const dy = ToReal(CCTK_DELTA_SPACE(1));
   CCTK_REAL const dz = ToReal(CCTK_DELTA_SPACE(2));
   CCTK_REAL const dt = ToReal(CCTK_DELTA_TIME);
+  CCTK_REAL const t = ToReal(cctk_time);
   CCTK_REAL const dxi = INV(dx);
   CCTK_REAL const dyi = INV(dy);
   CCTK_REAL const dzi = INV(dz);
@@ -74,49 +62,55 @@ static void eulerauto_cons_calc_conserved_Body(cGH const * restrict const cctkGH
   CCTK_REAL const p1ody = INV(dy);
   CCTK_REAL const p1odz = INV(dz);
   
+  /* Assign local copies of arrays functions */
+  
+  
+  
+  /* Calculate temporaries and arrays functions */
+  
+  /* Copy local copies back to grid functions */
+  
   /* Loop over the grid points */
-  for (k = min[2]; k < max[2]; k++)
+  #pragma omp parallel
+  CCTK_LOOP3 (eulerauto_cons_calc_conserved,
+    i,j,k, imin[0],imin[1],imin[2], imax[0],imax[1],imax[2],
+    cctk_lsh[0],cctk_lsh[1],cctk_lsh[2])
   {
-    for (j = min[1]; j < max[1]; j++)
-    {
-      for (i = min[0]; i < max[0]; i++)
-      {
-         int  const  index  =  CCTK_GFINDEX3D(cctkGH,i,j,k) ;
-        
-        /* Assign local copies of grid functions */
-        
-        CCTK_REAL pL = p[index];
-        CCTK_REAL rhoL = rho[index];
-        CCTK_REAL v1L = v1[index];
-        CCTK_REAL v2L = v2[index];
-        CCTK_REAL v3L = v3[index];
-        
-        
-        /* Include user supplied include files */
-        
-        /* Precompute derivatives */
-        
-        /* Calculate temporaries and grid functions */
-        CCTK_REAL DenL = rhoL;
-        
-        CCTK_REAL S1L = rhoL*v1L;
-        
-        CCTK_REAL S2L = rhoL*v2L;
-        
-        CCTK_REAL S3L = rhoL*v3L;
-        
-        CCTK_REAL EnL = pL*INV(-1 + ToReal(gamma)) + 0.5*rhoL*(SQR(v1L) + 
-          SQR(v2L) + SQR(v3L));
-        
-        /* Copy local copies back to grid functions */
-        Den[index] = DenL;
-        En[index] = EnL;
-        S1[index] = S1L;
-        S2[index] = S2L;
-        S3[index] = S3L;
-      }
-    }
+    ptrdiff_t const index = di*i + dj*j + dk*k;
+    
+    /* Assign local copies of grid functions */
+    
+    CCTK_REAL pL = p[index];
+    CCTK_REAL rhoL = rho[index];
+    CCTK_REAL v1L = v1[index];
+    CCTK_REAL v2L = v2[index];
+    CCTK_REAL v3L = v3[index];
+    
+    
+    /* Include user supplied include files */
+    
+    /* Precompute derivatives */
+    
+    /* Calculate temporaries and grid functions */
+    CCTK_REAL DenL = rhoL;
+    
+    CCTK_REAL S1L = rhoL*v1L;
+    
+    CCTK_REAL S2L = rhoL*v2L;
+    
+    CCTK_REAL S3L = rhoL*v3L;
+    
+    CCTK_REAL EnL = 0.5*INV(-1 + ToReal(gamma))*(2*pL + 
+      rhoL*(SQR(v1L) + SQR(v2L) + SQR(v3L))*(-1 + ToReal(gamma)));
+    
+    /* Copy local copies back to grid functions */
+    Den[index] = DenL;
+    En[index] = EnL;
+    S1[index] = S1L;
+    S2[index] = S2L;
+    S3[index] = S3L;
   }
+  CCTK_ENDLOOP3 (eulerauto_cons_calc_conserved);
 }
 
 extern "C" void eulerauto_cons_calc_conserved(CCTK_ARGUMENTS)
@@ -124,5 +118,25 @@ extern "C" void eulerauto_cons_calc_conserved(CCTK_ARGUMENTS)
   DECLARE_CCTK_ARGUMENTS;
   DECLARE_CCTK_PARAMETERS;
   
+  
+  if (verbose > 1)
+  {
+    CCTK_VInfo(CCTK_THORNSTRING,"Entering eulerauto_cons_calc_conserved_Body");
+  }
+  
+  if (cctk_iteration % eulerauto_cons_calc_conserved_calc_every != eulerauto_cons_calc_conserved_calc_offset)
+  {
+    return;
+  }
+  
+  const char *groups[] = {"EulerAuto::Den_group","EulerAuto::En_group","EulerAuto::p_group","EulerAuto::rho_group","EulerAuto::S_group","EulerAuto::v_group"};
+  GenericFD_AssertGroupStorage(cctkGH, "eulerauto_cons_calc_conserved", 6, groups);
+  
+  
   GenericFD_LoopOverEverything(cctkGH, &eulerauto_cons_calc_conserved_Body);
+  
+  if (verbose > 1)
+  {
+    CCTK_VInfo(CCTK_THORNSTRING,"Leaving eulerauto_cons_calc_conserved_Body");
+  }
 }
diff --git a/Examples/EulerAuto/src/eulerauto_cons_calc_flux_1.cc b/Examples/EulerAuto/src/eulerauto_cons_calc_flux_1.cc
index 579c9dc..02eec61 100644
--- a/Examples/EulerAuto/src/eulerauto_cons_calc_flux_1.cc
+++ b/Examples/EulerAuto/src/eulerauto_cons_calc_flux_1.cc
@@ -12,6 +12,8 @@
 #include "cctk_Parameters.h"
 #include "GenericFD.h"
 #include "Differencing.h"
+#include "cctk_Loop.h"
+#include "loopcontrol.h"
 
 /* Define macros used in calculations */
 #define INITVALUE (42)
@@ -44,43 +46,28 @@ extern "C" void eulerauto_cons_calc_flux_1_SelectBCs(CCTK_ARGUMENTS)
   return;
 }
 
-static void eulerauto_cons_calc_flux_1_Body(cGH const * restrict const cctkGH, int const dir, int const face, CCTK_REAL const normal[3], CCTK_REAL const tangentA[3], CCTK_REAL const tangentB[3], int const min[3], int const max[3], int const n_subblock_gfs, CCTK_REAL * restrict const subblock_gfs[])
+static void eulerauto_cons_calc_flux_1_Body(cGH const * restrict const cctkGH, int const dir, int const face, CCTK_REAL const normal[3], CCTK_REAL const tangentA[3], CCTK_REAL const tangentB[3], int const imin[3], int const imax[3], int const n_subblock_gfs, CCTK_REAL * restrict const subblock_gfs[])
 {
   DECLARE_CCTK_ARGUMENTS;
   DECLARE_CCTK_PARAMETERS;
   
   
-  /* Declare the variables used for looping over grid points */
-  CCTK_INT i, j, k;
-  // CCTK_INT index = INITVALUE;
-  
   /* Declare finite differencing variables */
   
-  if (verbose > 1)
-  {
-    CCTK_VInfo(CCTK_THORNSTRING,"Entering eulerauto_cons_calc_flux_1_Body");
-  }
-  
-  if (cctk_iteration % eulerauto_cons_calc_flux_1_calc_every != eulerauto_cons_calc_flux_1_calc_offset)
-  {
-    return;
-  }
-  
-  const char *groups[] = {"EulerAuto::Den_flux_group","EulerAuto::Den_lr_group","EulerAuto::En_flux_group","EulerAuto::En_lr_group","EulerAuto::p_lr_group","EulerAuto::rho_lr_group","EulerAuto::S1_flux_group","EulerAuto::S1_lr_group","EulerAuto::S2_flux_group","EulerAuto::S2_lr_group","EulerAuto::S3_flux_group","EulerAuto::S3_lr_group","EulerAuto::v1_lr_group","EulerAuto::v2_lr_group","EulerAuto::v3_lr_group"};
-  GenericFD_AssertGroupStorage(cctkGH, "eulerauto_cons_calc_flux_1", 15, groups);
-  
-  GenericFD_EnsureStencilFits(cctkGH, "eulerauto_cons_calc_flux_1", 1, 1, 1);
-  
   /* Include user-supplied include files */
   
   /* Initialise finite differencing variables */
   ptrdiff_t const di = 1;
   ptrdiff_t const dj = CCTK_GFINDEX3D(cctkGH,0,1,0) - CCTK_GFINDEX3D(cctkGH,0,0,0);
   ptrdiff_t const dk = CCTK_GFINDEX3D(cctkGH,0,0,1) - CCTK_GFINDEX3D(cctkGH,0,0,0);
+  ptrdiff_t const cdi = sizeof(CCTK_REAL) * di;
+  ptrdiff_t const cdj = sizeof(CCTK_REAL) * dj;
+  ptrdiff_t const cdk = sizeof(CCTK_REAL) * dk;
   CCTK_REAL const dx = ToReal(CCTK_DELTA_SPACE(0));
   CCTK_REAL const dy = ToReal(CCTK_DELTA_SPACE(1));
   CCTK_REAL const dz = ToReal(CCTK_DELTA_SPACE(2));
   CCTK_REAL const dt = ToReal(CCTK_DELTA_TIME);
+  CCTK_REAL const t = ToReal(cctk_time);
   CCTK_REAL const dxi = INV(dx);
   CCTK_REAL const dyi = INV(dy);
   CCTK_REAL const dzi = INV(dz);
@@ -99,102 +86,108 @@ static void eulerauto_cons_calc_flux_1_Body(cGH const * restrict const cctkGH, i
   CCTK_REAL const p1ody = INV(dy);
   CCTK_REAL const p1odz = INV(dz);
   
+  /* Assign local copies of arrays functions */
+  
+  
+  
+  /* Calculate temporaries and arrays functions */
+  
+  /* Copy local copies back to grid functions */
+  
   /* Loop over the grid points */
-  for (k = min[2]; k < max[2]; k++)
+  #pragma omp parallel
+  CCTK_LOOP3 (eulerauto_cons_calc_flux_1,
+    i,j,k, imin[0],imin[1],imin[2], imax[0],imax[1],imax[2],
+    cctk_lsh[0],cctk_lsh[1],cctk_lsh[2])
   {
-    for (j = min[1]; j < max[1]; j++)
-    {
-      for (i = min[0]; i < max[0]; i++)
-      {
-         int  const  index  =  CCTK_GFINDEX3D(cctkGH,i,j,k) ;
-        
-        /* Assign local copies of grid functions */
-        
-        CCTK_REAL DenLeftL = DenLeft[index];
-        CCTK_REAL DenRightL = DenRight[index];
-        CCTK_REAL EnLeftL = EnLeft[index];
-        CCTK_REAL EnRightL = EnRight[index];
-        CCTK_REAL pLeftL = pLeft[index];
-        CCTK_REAL pRightL = pRight[index];
-        CCTK_REAL rhoLeftL = rhoLeft[index];
-        CCTK_REAL rhoRightL = rhoRight[index];
-        CCTK_REAL S1LeftL = S1Left[index];
-        CCTK_REAL S1RightL = S1Right[index];
-        CCTK_REAL S2LeftL = S2Left[index];
-        CCTK_REAL S2RightL = S2Right[index];
-        CCTK_REAL S3LeftL = S3Left[index];
-        CCTK_REAL S3RightL = S3Right[index];
-        CCTK_REAL v1LeftL = v1Left[index];
-        CCTK_REAL v1RightL = v1Right[index];
-        CCTK_REAL v2LeftL = v2Left[index];
-        CCTK_REAL v2RightL = v2Right[index];
-        CCTK_REAL v3LeftL = v3Left[index];
-        CCTK_REAL v3RightL = v3Right[index];
-        
-        
-        /* Include user supplied include files */
-        
-        /* Precompute derivatives */
-        CCTK_REAL const ShiftMinus1DenRight = ShiftMinus1(&DenRight[index]);
-        CCTK_REAL const ShiftMinus1EnRight = ShiftMinus1(&EnRight[index]);
-        CCTK_REAL const ShiftMinus1pRight = ShiftMinus1(&pRight[index]);
-        CCTK_REAL const ShiftMinus1rhoRight = ShiftMinus1(&rhoRight[index]);
-        CCTK_REAL const ShiftMinus1S1Right = ShiftMinus1(&S1Right[index]);
-        CCTK_REAL const ShiftMinus1S2Right = ShiftMinus1(&S2Right[index]);
-        CCTK_REAL const ShiftMinus1S3Right = ShiftMinus1(&S3Right[index]);
-        CCTK_REAL const ShiftMinus1v1Right = ShiftMinus1(&v1Right[index]);
-        CCTK_REAL const ShiftMinus1v2Right = ShiftMinus1(&v2Right[index]);
-        CCTK_REAL const ShiftMinus1v3Right = ShiftMinus1(&v3Right[index]);
-        
-        /* Calculate temporaries and grid functions */
-        CCTK_REAL DenFluxLeft = rhoLeftL*v1LeftL;
-        
-        CCTK_REAL DenFluxRight = ShiftMinus1rhoRight*ShiftMinus1v1Right;
-        
-        CCTK_REAL DenFluxL = 0.5*(DenFluxLeft + DenFluxRight + (-DenLeftL + 
-          ShiftMinus1DenRight)*ToReal(hlleAlpha));
-        
-        CCTK_REAL S1FluxLeft = pLeftL + rhoLeftL*SQR(v1LeftL);
-        
-        CCTK_REAL S1FluxRight = ShiftMinus1pRight + 
-          ShiftMinus1rhoRight*SQR(ShiftMinus1v1Right);
-        
-        CCTK_REAL S1FluxL = 0.5*(S1FluxLeft + S1FluxRight + (-S1LeftL + 
-          ShiftMinus1S1Right)*ToReal(hlleAlpha));
-        
-        CCTK_REAL S2FluxLeft = rhoLeftL*v1LeftL*v2LeftL;
-        
-        CCTK_REAL S2FluxRight = 
-          ShiftMinus1rhoRight*ShiftMinus1v1Right*ShiftMinus1v2Right;
-        
-        CCTK_REAL S2FluxL = 0.5*(S2FluxLeft + S2FluxRight + (-S2LeftL + 
-          ShiftMinus1S2Right)*ToReal(hlleAlpha));
-        
-        CCTK_REAL S3FluxLeft = rhoLeftL*v1LeftL*v3LeftL;
-        
-        CCTK_REAL S3FluxRight = 
-          ShiftMinus1rhoRight*ShiftMinus1v1Right*ShiftMinus1v3Right;
-        
-        CCTK_REAL S3FluxL = 0.5*(S3FluxLeft + S3FluxRight + (-S3LeftL + 
-          ShiftMinus1S3Right)*ToReal(hlleAlpha));
-        
-        CCTK_REAL EnFluxLeft = (EnLeftL + pLeftL)*v1LeftL;
-        
-        CCTK_REAL EnFluxRight = (ShiftMinus1EnRight + 
-          ShiftMinus1pRight)*ShiftMinus1v1Right;
-        
-        CCTK_REAL EnFluxL = 0.5*(EnFluxLeft + EnFluxRight + (-EnLeftL + 
-          ShiftMinus1EnRight)*ToReal(hlleAlpha));
-        
-        /* Copy local copies back to grid functions */
-        DenFlux[index] = DenFluxL;
-        EnFlux[index] = EnFluxL;
-        S1Flux[index] = S1FluxL;
-        S2Flux[index] = S2FluxL;
-        S3Flux[index] = S3FluxL;
-      }
-    }
+    ptrdiff_t const index = di*i + dj*j + dk*k;
+    
+    /* Assign local copies of grid functions */
+    
+    CCTK_REAL DenLeftL = DenLeft[index];
+    CCTK_REAL DenRightL = DenRight[index];
+    CCTK_REAL EnLeftL = EnLeft[index];
+    CCTK_REAL EnRightL = EnRight[index];
+    CCTK_REAL pLeftL = pLeft[index];
+    CCTK_REAL pRightL = pRight[index];
+    CCTK_REAL rhoLeftL = rhoLeft[index];
+    CCTK_REAL rhoRightL = rhoRight[index];
+    CCTK_REAL S1LeftL = S1Left[index];
+    CCTK_REAL S1RightL = S1Right[index];
+    CCTK_REAL S2LeftL = S2Left[index];
+    CCTK_REAL S2RightL = S2Right[index];
+    CCTK_REAL S3LeftL = S3Left[index];
+    CCTK_REAL S3RightL = S3Right[index];
+    CCTK_REAL v1LeftL = v1Left[index];
+    CCTK_REAL v1RightL = v1Right[index];
+    CCTK_REAL v2LeftL = v2Left[index];
+    CCTK_REAL v2RightL = v2Right[index];
+    CCTK_REAL v3LeftL = v3Left[index];
+    CCTK_REAL v3RightL = v3Right[index];
+    
+    
+    /* Include user supplied include files */
+    
+    /* Precompute derivatives */
+    CCTK_REAL const ShiftMinus1DenRight = ShiftMinus1(&DenRight[index]);
+    CCTK_REAL const ShiftMinus1EnRight = ShiftMinus1(&EnRight[index]);
+    CCTK_REAL const ShiftMinus1pRight = ShiftMinus1(&pRight[index]);
+    CCTK_REAL const ShiftMinus1rhoRight = ShiftMinus1(&rhoRight[index]);
+    CCTK_REAL const ShiftMinus1S1Right = ShiftMinus1(&S1Right[index]);
+    CCTK_REAL const ShiftMinus1S2Right = ShiftMinus1(&S2Right[index]);
+    CCTK_REAL const ShiftMinus1S3Right = ShiftMinus1(&S3Right[index]);
+    CCTK_REAL const ShiftMinus1v1Right = ShiftMinus1(&v1Right[index]);
+    CCTK_REAL const ShiftMinus1v2Right = ShiftMinus1(&v2Right[index]);
+    CCTK_REAL const ShiftMinus1v3Right = ShiftMinus1(&v3Right[index]);
+    
+    /* Calculate temporaries and grid functions */
+    CCTK_REAL DenFluxLeft = rhoLeftL*v1LeftL;
+    
+    CCTK_REAL DenFluxRight = ShiftMinus1rhoRight*ShiftMinus1v1Right;
+    
+    CCTK_REAL DenFluxL = 0.5*(DenFluxLeft + DenFluxRight + (-DenLeftL + 
+      ShiftMinus1DenRight)*ToReal(hlleAlpha));
+    
+    CCTK_REAL S1FluxLeft = pLeftL + rhoLeftL*SQR(v1LeftL);
+    
+    CCTK_REAL S1FluxRight = ShiftMinus1pRight + 
+      ShiftMinus1rhoRight*SQR(ShiftMinus1v1Right);
+    
+    CCTK_REAL S1FluxL = 0.5*(S1FluxLeft + S1FluxRight + (-S1LeftL + 
+      ShiftMinus1S1Right)*ToReal(hlleAlpha));
+    
+    CCTK_REAL S2FluxLeft = rhoLeftL*v1LeftL*v2LeftL;
+    
+    CCTK_REAL S2FluxRight = 
+      ShiftMinus1rhoRight*ShiftMinus1v1Right*ShiftMinus1v2Right;
+    
+    CCTK_REAL S2FluxL = 0.5*(S2FluxLeft + S2FluxRight + (-S2LeftL + 
+      ShiftMinus1S2Right)*ToReal(hlleAlpha));
+    
+    CCTK_REAL S3FluxLeft = rhoLeftL*v1LeftL*v3LeftL;
+    
+    CCTK_REAL S3FluxRight = 
+      ShiftMinus1rhoRight*ShiftMinus1v1Right*ShiftMinus1v3Right;
+    
+    CCTK_REAL S3FluxL = 0.5*(S3FluxLeft + S3FluxRight + (-S3LeftL + 
+      ShiftMinus1S3Right)*ToReal(hlleAlpha));
+    
+    CCTK_REAL EnFluxLeft = (EnLeftL + pLeftL)*v1LeftL;
+    
+    CCTK_REAL EnFluxRight = (ShiftMinus1EnRight + 
+      ShiftMinus1pRight)*ShiftMinus1v1Right;
+    
+    CCTK_REAL EnFluxL = 0.5*(EnFluxLeft + EnFluxRight + (-EnLeftL + 
+      ShiftMinus1EnRight)*ToReal(hlleAlpha));
+    
+    /* Copy local copies back to grid functions */
+    DenFlux[index] = DenFluxL;
+    EnFlux[index] = EnFluxL;
+    S1Flux[index] = S1FluxL;
+    S2Flux[index] = S2FluxL;
+    S3Flux[index] = S3FluxL;
   }
+  CCTK_ENDLOOP3 (eulerauto_cons_calc_flux_1);
 }
 
 extern "C" void eulerauto_cons_calc_flux_1(CCTK_ARGUMENTS)
@@ -202,5 +195,26 @@ extern "C" void eulerauto_cons_calc_flux_1(CCTK_ARGUMENTS)
   DECLARE_CCTK_ARGUMENTS;
   DECLARE_CCTK_PARAMETERS;
   
+  
+  if (verbose > 1)
+  {
+    CCTK_VInfo(CCTK_THORNSTRING,"Entering eulerauto_cons_calc_flux_1_Body");
+  }
+  
+  if (cctk_iteration % eulerauto_cons_calc_flux_1_calc_every != eulerauto_cons_calc_flux_1_calc_offset)
+  {
+    return;
+  }
+  
+  const char *groups[] = {"EulerAuto::Den_flux_group","EulerAuto::Den_lr_group","EulerAuto::En_flux_group","EulerAuto::En_lr_group","EulerAuto::p_lr_group","EulerAuto::rho_lr_group","EulerAuto::S1_flux_group","EulerAuto::S1_lr_group","EulerAuto::S2_flux_group","EulerAuto::S2_lr_group","EulerAuto::S3_flux_group","EulerAuto::S3_lr_group","EulerAuto::v1_lr_group","EulerAuto::v2_lr_group","EulerAuto::v3_lr_group"};
+  GenericFD_AssertGroupStorage(cctkGH, "eulerauto_cons_calc_flux_1", 15, groups);
+  
+  GenericFD_EnsureStencilFits(cctkGH, "eulerauto_cons_calc_flux_1", 1, 1, 1);
+  
   GenericFD_LoopOverInterior(cctkGH, &eulerauto_cons_calc_flux_1_Body);
+  
+  if (verbose > 1)
+  {
+    CCTK_VInfo(CCTK_THORNSTRING,"Leaving eulerauto_cons_calc_flux_1_Body");
+  }
 }
diff --git a/Examples/EulerAuto/src/eulerauto_cons_calc_flux_2.cc b/Examples/EulerAuto/src/eulerauto_cons_calc_flux_2.cc
index cf67145..eb6fac5 100644
--- a/Examples/EulerAuto/src/eulerauto_cons_calc_flux_2.cc
+++ b/Examples/EulerAuto/src/eulerauto_cons_calc_flux_2.cc
@@ -12,6 +12,8 @@
 #include "cctk_Parameters.h"
 #include "GenericFD.h"
 #include "Differencing.h"
+#include "cctk_Loop.h"
+#include "loopcontrol.h"
 
 /* Define macros used in calculations */
 #define INITVALUE (42)
@@ -44,43 +46,28 @@ extern "C" void eulerauto_cons_calc_flux_2_SelectBCs(CCTK_ARGUMENTS)
   return;
 }
 
-static void eulerauto_cons_calc_flux_2_Body(cGH const * restrict const cctkGH, int const dir, int const face, CCTK_REAL const normal[3], CCTK_REAL const tangentA[3], CCTK_REAL const tangentB[3], int const min[3], int const max[3], int const n_subblock_gfs, CCTK_REAL * restrict const subblock_gfs[])
+static void eulerauto_cons_calc_flux_2_Body(cGH const * restrict const cctkGH, int const dir, int const face, CCTK_REAL const normal[3], CCTK_REAL const tangentA[3], CCTK_REAL const tangentB[3], int const imin[3], int const imax[3], int const n_subblock_gfs, CCTK_REAL * restrict const subblock_gfs[])
 {
   DECLARE_CCTK_ARGUMENTS;
   DECLARE_CCTK_PARAMETERS;
   
   
-  /* Declare the variables used for looping over grid points */
-  CCTK_INT i, j, k;
-  // CCTK_INT index = INITVALUE;
-  
   /* Declare finite differencing variables */
   
-  if (verbose > 1)
-  {
-    CCTK_VInfo(CCTK_THORNSTRING,"Entering eulerauto_cons_calc_flux_2_Body");
-  }
-  
-  if (cctk_iteration % eulerauto_cons_calc_flux_2_calc_every != eulerauto_cons_calc_flux_2_calc_offset)
-  {
-    return;
-  }
-  
-  const char *groups[] = {"EulerAuto::Den_flux_group","EulerAuto::Den_lr_group","EulerAuto::En_flux_group","EulerAuto::En_lr_group","EulerAuto::p_lr_group","EulerAuto::rho_lr_group","EulerAuto::S1_flux_group","EulerAuto::S1_lr_group","EulerAuto::S2_flux_group","EulerAuto::S2_lr_group","EulerAuto::S3_flux_group","EulerAuto::S3_lr_group","EulerAuto::v1_lr_group","EulerAuto::v2_lr_group","EulerAuto::v3_lr_group"};
-  GenericFD_AssertGroupStorage(cctkGH, "eulerauto_cons_calc_flux_2", 15, groups);
-  
-  GenericFD_EnsureStencilFits(cctkGH, "eulerauto_cons_calc_flux_2", 1, 1, 1);
-  
   /* Include user-supplied include files */
   
   /* Initialise finite differencing variables */
   ptrdiff_t const di = 1;
   ptrdiff_t const dj = CCTK_GFINDEX3D(cctkGH,0,1,0) - CCTK_GFINDEX3D(cctkGH,0,0,0);
   ptrdiff_t const dk = CCTK_GFINDEX3D(cctkGH,0,0,1) - CCTK_GFINDEX3D(cctkGH,0,0,0);
+  ptrdiff_t const cdi = sizeof(CCTK_REAL) * di;
+  ptrdiff_t const cdj = sizeof(CCTK_REAL) * dj;
+  ptrdiff_t const cdk = sizeof(CCTK_REAL) * dk;
   CCTK_REAL const dx = ToReal(CCTK_DELTA_SPACE(0));
   CCTK_REAL const dy = ToReal(CCTK_DELTA_SPACE(1));
   CCTK_REAL const dz = ToReal(CCTK_DELTA_SPACE(2));
   CCTK_REAL const dt = ToReal(CCTK_DELTA_TIME);
+  CCTK_REAL const t = ToReal(cctk_time);
   CCTK_REAL const dxi = INV(dx);
   CCTK_REAL const dyi = INV(dy);
   CCTK_REAL const dzi = INV(dz);
@@ -99,102 +86,108 @@ static void eulerauto_cons_calc_flux_2_Body(cGH const * restrict const cctkGH, i
   CCTK_REAL const p1ody = INV(dy);
   CCTK_REAL const p1odz = INV(dz);
   
+  /* Assign local copies of arrays functions */
+  
+  
+  
+  /* Calculate temporaries and arrays functions */
+  
+  /* Copy local copies back to grid functions */
+  
   /* Loop over the grid points */
-  for (k = min[2]; k < max[2]; k++)
+  #pragma omp parallel
+  CCTK_LOOP3 (eulerauto_cons_calc_flux_2,
+    i,j,k, imin[0],imin[1],imin[2], imax[0],imax[1],imax[2],
+    cctk_lsh[0],cctk_lsh[1],cctk_lsh[2])
   {
-    for (j = min[1]; j < max[1]; j++)
-    {
-      for (i = min[0]; i < max[0]; i++)
-      {
-         int  const  index  =  CCTK_GFINDEX3D(cctkGH,i,j,k) ;
-        
-        /* Assign local copies of grid functions */
-        
-        CCTK_REAL DenLeftL = DenLeft[index];
-        CCTK_REAL DenRightL = DenRight[index];
-        CCTK_REAL EnLeftL = EnLeft[index];
-        CCTK_REAL EnRightL = EnRight[index];
-        CCTK_REAL pLeftL = pLeft[index];
-        CCTK_REAL pRightL = pRight[index];
-        CCTK_REAL rhoLeftL = rhoLeft[index];
-        CCTK_REAL rhoRightL = rhoRight[index];
-        CCTK_REAL S1LeftL = S1Left[index];
-        CCTK_REAL S1RightL = S1Right[index];
-        CCTK_REAL S2LeftL = S2Left[index];
-        CCTK_REAL S2RightL = S2Right[index];
-        CCTK_REAL S3LeftL = S3Left[index];
-        CCTK_REAL S3RightL = S3Right[index];
-        CCTK_REAL v1LeftL = v1Left[index];
-        CCTK_REAL v1RightL = v1Right[index];
-        CCTK_REAL v2LeftL = v2Left[index];
-        CCTK_REAL v2RightL = v2Right[index];
-        CCTK_REAL v3LeftL = v3Left[index];
-        CCTK_REAL v3RightL = v3Right[index];
-        
-        
-        /* Include user supplied include files */
-        
-        /* Precompute derivatives */
-        CCTK_REAL const ShiftMinus2DenRight = ShiftMinus2(&DenRight[index]);
-        CCTK_REAL const ShiftMinus2EnRight = ShiftMinus2(&EnRight[index]);
-        CCTK_REAL const ShiftMinus2pRight = ShiftMinus2(&pRight[index]);
-        CCTK_REAL const ShiftMinus2rhoRight = ShiftMinus2(&rhoRight[index]);
-        CCTK_REAL const ShiftMinus2S1Right = ShiftMinus2(&S1Right[index]);
-        CCTK_REAL const ShiftMinus2S2Right = ShiftMinus2(&S2Right[index]);
-        CCTK_REAL const ShiftMinus2S3Right = ShiftMinus2(&S3Right[index]);
-        CCTK_REAL const ShiftMinus2v1Right = ShiftMinus2(&v1Right[index]);
-        CCTK_REAL const ShiftMinus2v2Right = ShiftMinus2(&v2Right[index]);
-        CCTK_REAL const ShiftMinus2v3Right = ShiftMinus2(&v3Right[index]);
-        
-        /* Calculate temporaries and grid functions */
-        CCTK_REAL DenFluxLeft = rhoLeftL*v2LeftL;
-        
-        CCTK_REAL DenFluxRight = ShiftMinus2rhoRight*ShiftMinus2v2Right;
-        
-        CCTK_REAL DenFluxL = 0.5*(DenFluxLeft + DenFluxRight + (-DenLeftL + 
-          ShiftMinus2DenRight)*ToReal(hlleAlpha));
-        
-        CCTK_REAL S1FluxLeft = rhoLeftL*v1LeftL*v2LeftL;
-        
-        CCTK_REAL S1FluxRight = 
-          ShiftMinus2rhoRight*ShiftMinus2v1Right*ShiftMinus2v2Right;
-        
-        CCTK_REAL S1FluxL = 0.5*(S1FluxLeft + S1FluxRight + (-S1LeftL + 
-          ShiftMinus2S1Right)*ToReal(hlleAlpha));
-        
-        CCTK_REAL S2FluxLeft = pLeftL + rhoLeftL*SQR(v2LeftL);
-        
-        CCTK_REAL S2FluxRight = ShiftMinus2pRight + 
-          ShiftMinus2rhoRight*SQR(ShiftMinus2v2Right);
-        
-        CCTK_REAL S2FluxL = 0.5*(S2FluxLeft + S2FluxRight + (-S2LeftL + 
-          ShiftMinus2S2Right)*ToReal(hlleAlpha));
-        
-        CCTK_REAL S3FluxLeft = rhoLeftL*v2LeftL*v3LeftL;
-        
-        CCTK_REAL S3FluxRight = 
-          ShiftMinus2rhoRight*ShiftMinus2v2Right*ShiftMinus2v3Right;
-        
-        CCTK_REAL S3FluxL = 0.5*(S3FluxLeft + S3FluxRight + (-S3LeftL + 
-          ShiftMinus2S3Right)*ToReal(hlleAlpha));
-        
-        CCTK_REAL EnFluxLeft = (EnLeftL + pLeftL)*v2LeftL;
-        
-        CCTK_REAL EnFluxRight = (ShiftMinus2EnRight + 
-          ShiftMinus2pRight)*ShiftMinus2v2Right;
-        
-        CCTK_REAL EnFluxL = 0.5*(EnFluxLeft + EnFluxRight + (-EnLeftL + 
-          ShiftMinus2EnRight)*ToReal(hlleAlpha));
-        
-        /* Copy local copies back to grid functions */
-        DenFlux[index] = DenFluxL;
-        EnFlux[index] = EnFluxL;
-        S1Flux[index] = S1FluxL;
-        S2Flux[index] = S2FluxL;
-        S3Flux[index] = S3FluxL;
-      }
-    }
+    ptrdiff_t const index = di*i + dj*j + dk*k;
+    
+    /* Assign local copies of grid functions */
+    
+    CCTK_REAL DenLeftL = DenLeft[index];
+    CCTK_REAL DenRightL = DenRight[index];
+    CCTK_REAL EnLeftL = EnLeft[index];
+    CCTK_REAL EnRightL = EnRight[index];
+    CCTK_REAL pLeftL = pLeft[index];
+    CCTK_REAL pRightL = pRight[index];
+    CCTK_REAL rhoLeftL = rhoLeft[index];
+    CCTK_REAL rhoRightL = rhoRight[index];
+    CCTK_REAL S1LeftL = S1Left[index];
+    CCTK_REAL S1RightL = S1Right[index];
+    CCTK_REAL S2LeftL = S2Left[index];
+    CCTK_REAL S2RightL = S2Right[index];
+    CCTK_REAL S3LeftL = S3Left[index];
+    CCTK_REAL S3RightL = S3Right[index];
+    CCTK_REAL v1LeftL = v1Left[index];
+    CCTK_REAL v1RightL = v1Right[index];
+    CCTK_REAL v2LeftL = v2Left[index];
+    CCTK_REAL v2RightL = v2Right[index];
+    CCTK_REAL v3LeftL = v3Left[index];
+    CCTK_REAL v3RightL = v3Right[index];
+    
+    
+    /* Include user supplied include files */
+    
+    /* Precompute derivatives */
+    CCTK_REAL const ShiftMinus2DenRight = ShiftMinus2(&DenRight[index]);
+    CCTK_REAL const ShiftMinus2EnRight = ShiftMinus2(&EnRight[index]);
+    CCTK_REAL const ShiftMinus2pRight = ShiftMinus2(&pRight[index]);
+    CCTK_REAL const ShiftMinus2rhoRight = ShiftMinus2(&rhoRight[index]);
+    CCTK_REAL const ShiftMinus2S1Right = ShiftMinus2(&S1Right[index]);
+    CCTK_REAL const ShiftMinus2S2Right = ShiftMinus2(&S2Right[index]);
+    CCTK_REAL const ShiftMinus2S3Right = ShiftMinus2(&S3Right[index]);
+    CCTK_REAL const ShiftMinus2v1Right = ShiftMinus2(&v1Right[index]);
+    CCTK_REAL const ShiftMinus2v2Right = ShiftMinus2(&v2Right[index]);
+    CCTK_REAL const ShiftMinus2v3Right = ShiftMinus2(&v3Right[index]);
+    
+    /* Calculate temporaries and grid functions */
+    CCTK_REAL DenFluxLeft = rhoLeftL*v2LeftL;
+    
+    CCTK_REAL DenFluxRight = ShiftMinus2rhoRight*ShiftMinus2v2Right;
+    
+    CCTK_REAL DenFluxL = 0.5*(DenFluxLeft + DenFluxRight + (-DenLeftL + 
+      ShiftMinus2DenRight)*ToReal(hlleAlpha));
+    
+    CCTK_REAL S1FluxLeft = rhoLeftL*v1LeftL*v2LeftL;
+    
+    CCTK_REAL S1FluxRight = 
+      ShiftMinus2rhoRight*ShiftMinus2v1Right*ShiftMinus2v2Right;
+    
+    CCTK_REAL S1FluxL = 0.5*(S1FluxLeft + S1FluxRight + (-S1LeftL + 
+      ShiftMinus2S1Right)*ToReal(hlleAlpha));
+    
+    CCTK_REAL S2FluxLeft = pLeftL + rhoLeftL*SQR(v2LeftL);
+    
+    CCTK_REAL S2FluxRight = ShiftMinus2pRight + 
+      ShiftMinus2rhoRight*SQR(ShiftMinus2v2Right);
+    
+    CCTK_REAL S2FluxL = 0.5*(S2FluxLeft + S2FluxRight + (-S2LeftL + 
+      ShiftMinus2S2Right)*ToReal(hlleAlpha));
+    
+    CCTK_REAL S3FluxLeft = rhoLeftL*v2LeftL*v3LeftL;
+    
+    CCTK_REAL S3FluxRight = 
+      ShiftMinus2rhoRight*ShiftMinus2v2Right*ShiftMinus2v3Right;
+    
+    CCTK_REAL S3FluxL = 0.5*(S3FluxLeft + S3FluxRight + (-S3LeftL + 
+      ShiftMinus2S3Right)*ToReal(hlleAlpha));
+    
+    CCTK_REAL EnFluxLeft = (EnLeftL + pLeftL)*v2LeftL;
+    
+    CCTK_REAL EnFluxRight = (ShiftMinus2EnRight + 
+      ShiftMinus2pRight)*ShiftMinus2v2Right;
+    
+    CCTK_REAL EnFluxL = 0.5*(EnFluxLeft + EnFluxRight + (-EnLeftL + 
+      ShiftMinus2EnRight)*ToReal(hlleAlpha));
+    
+    /* Copy local copies back to grid functions */
+    DenFlux[index] = DenFluxL;
+    EnFlux[index] = EnFluxL;
+    S1Flux[index] = S1FluxL;
+    S2Flux[index] = S2FluxL;
+    S3Flux[index] = S3FluxL;
   }
+  CCTK_ENDLOOP3 (eulerauto_cons_calc_flux_2);
 }
 
 extern "C" void eulerauto_cons_calc_flux_2(CCTK_ARGUMENTS)
@@ -202,5 +195,26 @@ extern "C" void eulerauto_cons_calc_flux_2(CCTK_ARGUMENTS)
   DECLARE_CCTK_ARGUMENTS;
   DECLARE_CCTK_PARAMETERS;
   
+  
+  if (verbose > 1)
+  {
+    CCTK_VInfo(CCTK_THORNSTRING,"Entering eulerauto_cons_calc_flux_2_Body");
+  }
+  
+  if (cctk_iteration % eulerauto_cons_calc_flux_2_calc_every != eulerauto_cons_calc_flux_2_calc_offset)
+  {
+    return;
+  }
+  
+  const char *groups[] = {"EulerAuto::Den_flux_group","EulerAuto::Den_lr_group","EulerAuto::En_flux_group","EulerAuto::En_lr_group","EulerAuto::p_lr_group","EulerAuto::rho_lr_group","EulerAuto::S1_flux_group","EulerAuto::S1_lr_group","EulerAuto::S2_flux_group","EulerAuto::S2_lr_group","EulerAuto::S3_flux_group","EulerAuto::S3_lr_group","EulerAuto::v1_lr_group","EulerAuto::v2_lr_group","EulerAuto::v3_lr_group"};
+  GenericFD_AssertGroupStorage(cctkGH, "eulerauto_cons_calc_flux_2", 15, groups);
+  
+  GenericFD_EnsureStencilFits(cctkGH, "eulerauto_cons_calc_flux_2", 1, 1, 1);
+  
   GenericFD_LoopOverInterior(cctkGH, &eulerauto_cons_calc_flux_2_Body);
+  
+  if (verbose > 1)
+  {
+    CCTK_VInfo(CCTK_THORNSTRING,"Leaving eulerauto_cons_calc_flux_2_Body");
+  }
 }
diff --git a/Examples/EulerAuto/src/eulerauto_cons_calc_flux_3.cc b/Examples/EulerAuto/src/eulerauto_cons_calc_flux_3.cc
index 797da55..140cd47 100644
--- a/Examples/EulerAuto/src/eulerauto_cons_calc_flux_3.cc
+++ b/Examples/EulerAuto/src/eulerauto_cons_calc_flux_3.cc
@@ -12,6 +12,8 @@
 #include "cctk_Parameters.h"
 #include "GenericFD.h"
 #include "Differencing.h"
+#include "cctk_Loop.h"
+#include "loopcontrol.h"
 
 /* Define macros used in calculations */
 #define INITVALUE (42)
@@ -44,43 +46,28 @@ extern "C" void eulerauto_cons_calc_flux_3_SelectBCs(CCTK_ARGUMENTS)
   return;
 }
 
-static void eulerauto_cons_calc_flux_3_Body(cGH const * restrict const cctkGH, int const dir, int const face, CCTK_REAL const normal[3], CCTK_REAL const tangentA[3], CCTK_REAL const tangentB[3], int const min[3], int const max[3], int const n_subblock_gfs, CCTK_REAL * restrict const subblock_gfs[])
+static void eulerauto_cons_calc_flux_3_Body(cGH const * restrict const cctkGH, int const dir, int const face, CCTK_REAL const normal[3], CCTK_REAL const tangentA[3], CCTK_REAL const tangentB[3], int const imin[3], int const imax[3], int const n_subblock_gfs, CCTK_REAL * restrict const subblock_gfs[])
 {
   DECLARE_CCTK_ARGUMENTS;
   DECLARE_CCTK_PARAMETERS;
   
   
-  /* Declare the variables used for looping over grid points */
-  CCTK_INT i, j, k;
-  // CCTK_INT index = INITVALUE;
-  
   /* Declare finite differencing variables */
   
-  if (verbose > 1)
-  {
-    CCTK_VInfo(CCTK_THORNSTRING,"Entering eulerauto_cons_calc_flux_3_Body");
-  }
-  
-  if (cctk_iteration % eulerauto_cons_calc_flux_3_calc_every != eulerauto_cons_calc_flux_3_calc_offset)
-  {
-    return;
-  }
-  
-  const char *groups[] = {"EulerAuto::Den_flux_group","EulerAuto::Den_lr_group","EulerAuto::En_flux_group","EulerAuto::En_lr_group","EulerAuto::p_lr_group","EulerAuto::rho_lr_group","EulerAuto::S1_flux_group","EulerAuto::S1_lr_group","EulerAuto::S2_flux_group","EulerAuto::S2_lr_group","EulerAuto::S3_flux_group","EulerAuto::S3_lr_group","EulerAuto::v1_lr_group","EulerAuto::v2_lr_group","EulerAuto::v3_lr_group"};
-  GenericFD_AssertGroupStorage(cctkGH, "eulerauto_cons_calc_flux_3", 15, groups);
-  
-  GenericFD_EnsureStencilFits(cctkGH, "eulerauto_cons_calc_flux_3", 1, 1, 1);
-  
   /* Include user-supplied include files */
   
   /* Initialise finite differencing variables */
   ptrdiff_t const di = 1;
   ptrdiff_t const dj = CCTK_GFINDEX3D(cctkGH,0,1,0) - CCTK_GFINDEX3D(cctkGH,0,0,0);
   ptrdiff_t const dk = CCTK_GFINDEX3D(cctkGH,0,0,1) - CCTK_GFINDEX3D(cctkGH,0,0,0);
+  ptrdiff_t const cdi = sizeof(CCTK_REAL) * di;
+  ptrdiff_t const cdj = sizeof(CCTK_REAL) * dj;
+  ptrdiff_t const cdk = sizeof(CCTK_REAL) * dk;
   CCTK_REAL const dx = ToReal(CCTK_DELTA_SPACE(0));
   CCTK_REAL const dy = ToReal(CCTK_DELTA_SPACE(1));
   CCTK_REAL const dz = ToReal(CCTK_DELTA_SPACE(2));
   CCTK_REAL const dt = ToReal(CCTK_DELTA_TIME);
+  CCTK_REAL const t = ToReal(cctk_time);
   CCTK_REAL const dxi = INV(dx);
   CCTK_REAL const dyi = INV(dy);
   CCTK_REAL const dzi = INV(dz);
@@ -99,102 +86,108 @@ static void eulerauto_cons_calc_flux_3_Body(cGH const * restrict const cctkGH, i
   CCTK_REAL const p1ody = INV(dy);
   CCTK_REAL const p1odz = INV(dz);
   
+  /* Assign local copies of arrays functions */
+  
+  
+  
+  /* Calculate temporaries and arrays functions */
+  
+  /* Copy local copies back to grid functions */
+  
   /* Loop over the grid points */
-  for (k = min[2]; k < max[2]; k++)
+  #pragma omp parallel
+  CCTK_LOOP3 (eulerauto_cons_calc_flux_3,
+    i,j,k, imin[0],imin[1],imin[2], imax[0],imax[1],imax[2],
+    cctk_lsh[0],cctk_lsh[1],cctk_lsh[2])
   {
-    for (j = min[1]; j < max[1]; j++)
-    {
-      for (i = min[0]; i < max[0]; i++)
-      {
-         int  const  index  =  CCTK_GFINDEX3D(cctkGH,i,j,k) ;
-        
-        /* Assign local copies of grid functions */
-        
-        CCTK_REAL DenLeftL = DenLeft[index];
-        CCTK_REAL DenRightL = DenRight[index];
-        CCTK_REAL EnLeftL = EnLeft[index];
-        CCTK_REAL EnRightL = EnRight[index];
-        CCTK_REAL pLeftL = pLeft[index];
-        CCTK_REAL pRightL = pRight[index];
-        CCTK_REAL rhoLeftL = rhoLeft[index];
-        CCTK_REAL rhoRightL = rhoRight[index];
-        CCTK_REAL S1LeftL = S1Left[index];
-        CCTK_REAL S1RightL = S1Right[index];
-        CCTK_REAL S2LeftL = S2Left[index];
-        CCTK_REAL S2RightL = S2Right[index];
-        CCTK_REAL S3LeftL = S3Left[index];
-        CCTK_REAL S3RightL = S3Right[index];
-        CCTK_REAL v1LeftL = v1Left[index];
-        CCTK_REAL v1RightL = v1Right[index];
-        CCTK_REAL v2LeftL = v2Left[index];
-        CCTK_REAL v2RightL = v2Right[index];
-        CCTK_REAL v3LeftL = v3Left[index];
-        CCTK_REAL v3RightL = v3Right[index];
-        
-        
-        /* Include user supplied include files */
-        
-        /* Precompute derivatives */
-        CCTK_REAL const ShiftMinus3DenRight = ShiftMinus3(&DenRight[index]);
-        CCTK_REAL const ShiftMinus3EnRight = ShiftMinus3(&EnRight[index]);
-        CCTK_REAL const ShiftMinus3pRight = ShiftMinus3(&pRight[index]);
-        CCTK_REAL const ShiftMinus3rhoRight = ShiftMinus3(&rhoRight[index]);
-        CCTK_REAL const ShiftMinus3S1Right = ShiftMinus3(&S1Right[index]);
-        CCTK_REAL const ShiftMinus3S2Right = ShiftMinus3(&S2Right[index]);
-        CCTK_REAL const ShiftMinus3S3Right = ShiftMinus3(&S3Right[index]);
-        CCTK_REAL const ShiftMinus3v1Right = ShiftMinus3(&v1Right[index]);
-        CCTK_REAL const ShiftMinus3v2Right = ShiftMinus3(&v2Right[index]);
-        CCTK_REAL const ShiftMinus3v3Right = ShiftMinus3(&v3Right[index]);
-        
-        /* Calculate temporaries and grid functions */
-        CCTK_REAL DenFluxLeft = rhoLeftL*v3LeftL;
-        
-        CCTK_REAL DenFluxRight = ShiftMinus3rhoRight*ShiftMinus3v3Right;
-        
-        CCTK_REAL DenFluxL = 0.5*(DenFluxLeft + DenFluxRight + (-DenLeftL + 
-          ShiftMinus3DenRight)*ToReal(hlleAlpha));
-        
-        CCTK_REAL S1FluxLeft = rhoLeftL*v1LeftL*v3LeftL;
-        
-        CCTK_REAL S1FluxRight = 
-          ShiftMinus3rhoRight*ShiftMinus3v1Right*ShiftMinus3v3Right;
-        
-        CCTK_REAL S1FluxL = 0.5*(S1FluxLeft + S1FluxRight + (-S1LeftL + 
-          ShiftMinus3S1Right)*ToReal(hlleAlpha));
-        
-        CCTK_REAL S2FluxLeft = rhoLeftL*v2LeftL*v3LeftL;
-        
-        CCTK_REAL S2FluxRight = 
-          ShiftMinus3rhoRight*ShiftMinus3v2Right*ShiftMinus3v3Right;
-        
-        CCTK_REAL S2FluxL = 0.5*(S2FluxLeft + S2FluxRight + (-S2LeftL + 
-          ShiftMinus3S2Right)*ToReal(hlleAlpha));
-        
-        CCTK_REAL S3FluxLeft = pLeftL + rhoLeftL*SQR(v3LeftL);
-        
-        CCTK_REAL S3FluxRight = ShiftMinus3pRight + 
-          ShiftMinus3rhoRight*SQR(ShiftMinus3v3Right);
-        
-        CCTK_REAL S3FluxL = 0.5*(S3FluxLeft + S3FluxRight + (-S3LeftL + 
-          ShiftMinus3S3Right)*ToReal(hlleAlpha));
-        
-        CCTK_REAL EnFluxLeft = (EnLeftL + pLeftL)*v3LeftL;
-        
-        CCTK_REAL EnFluxRight = (ShiftMinus3EnRight + 
-          ShiftMinus3pRight)*ShiftMinus3v3Right;
-        
-        CCTK_REAL EnFluxL = 0.5*(EnFluxLeft + EnFluxRight + (-EnLeftL + 
-          ShiftMinus3EnRight)*ToReal(hlleAlpha));
-        
-        /* Copy local copies back to grid functions */
-        DenFlux[index] = DenFluxL;
-        EnFlux[index] = EnFluxL;
-        S1Flux[index] = S1FluxL;
-        S2Flux[index] = S2FluxL;
-        S3Flux[index] = S3FluxL;
-      }
-    }
+    ptrdiff_t const index = di*i + dj*j + dk*k;
+    
+    /* Assign local copies of grid functions */
+    
+    CCTK_REAL DenLeftL = DenLeft[index];
+    CCTK_REAL DenRightL = DenRight[index];
+    CCTK_REAL EnLeftL = EnLeft[index];
+    CCTK_REAL EnRightL = EnRight[index];
+    CCTK_REAL pLeftL = pLeft[index];
+    CCTK_REAL pRightL = pRight[index];
+    CCTK_REAL rhoLeftL = rhoLeft[index];
+    CCTK_REAL rhoRightL = rhoRight[index];
+    CCTK_REAL S1LeftL = S1Left[index];
+    CCTK_REAL S1RightL = S1Right[index];
+    CCTK_REAL S2LeftL = S2Left[index];
+    CCTK_REAL S2RightL = S2Right[index];
+    CCTK_REAL S3LeftL = S3Left[index];
+    CCTK_REAL S3RightL = S3Right[index];
+    CCTK_REAL v1LeftL = v1Left[index];
+    CCTK_REAL v1RightL = v1Right[index];
+    CCTK_REAL v2LeftL = v2Left[index];
+    CCTK_REAL v2RightL = v2Right[index];
+    CCTK_REAL v3LeftL = v3Left[index];
+    CCTK_REAL v3RightL = v3Right[index];
+    
+    
+    /* Include user supplied include files */
+    
+    /* Precompute derivatives */
+    CCTK_REAL const ShiftMinus3DenRight = ShiftMinus3(&DenRight[index]);
+    CCTK_REAL const ShiftMinus3EnRight = ShiftMinus3(&EnRight[index]);
+    CCTK_REAL const ShiftMinus3pRight = ShiftMinus3(&pRight[index]);
+    CCTK_REAL const ShiftMinus3rhoRight = ShiftMinus3(&rhoRight[index]);
+    CCTK_REAL const ShiftMinus3S1Right = ShiftMinus3(&S1Right[index]);
+    CCTK_REAL const ShiftMinus3S2Right = ShiftMinus3(&S2Right[index]);
+    CCTK_REAL const ShiftMinus3S3Right = ShiftMinus3(&S3Right[index]);
+    CCTK_REAL const ShiftMinus3v1Right = ShiftMinus3(&v1Right[index]);
+    CCTK_REAL const ShiftMinus3v2Right = ShiftMinus3(&v2Right[index]);
+    CCTK_REAL const ShiftMinus3v3Right = ShiftMinus3(&v3Right[index]);
+    
+    /* Calculate temporaries and grid functions */
+    CCTK_REAL DenFluxLeft = rhoLeftL*v3LeftL;
+    
+    CCTK_REAL DenFluxRight = ShiftMinus3rhoRight*ShiftMinus3v3Right;
+    
+    CCTK_REAL DenFluxL = 0.5*(DenFluxLeft + DenFluxRight + (-DenLeftL + 
+      ShiftMinus3DenRight)*ToReal(hlleAlpha));
+    
+    CCTK_REAL S1FluxLeft = rhoLeftL*v1LeftL*v3LeftL;
+    
+    CCTK_REAL S1FluxRight = 
+      ShiftMinus3rhoRight*ShiftMinus3v1Right*ShiftMinus3v3Right;
+    
+    CCTK_REAL S1FluxL = 0.5*(S1FluxLeft + S1FluxRight + (-S1LeftL + 
+      ShiftMinus3S1Right)*ToReal(hlleAlpha));
+    
+    CCTK_REAL S2FluxLeft = rhoLeftL*v2LeftL*v3LeftL;
+    
+    CCTK_REAL S2FluxRight = 
+      ShiftMinus3rhoRight*ShiftMinus3v2Right*ShiftMinus3v3Right;
+    
+    CCTK_REAL S2FluxL = 0.5*(S2FluxLeft + S2FluxRight + (-S2LeftL + 
+      ShiftMinus3S2Right)*ToReal(hlleAlpha));
+    
+    CCTK_REAL S3FluxLeft = pLeftL + rhoLeftL*SQR(v3LeftL);
+    
+    CCTK_REAL S3FluxRight = ShiftMinus3pRight + 
+      ShiftMinus3rhoRight*SQR(ShiftMinus3v3Right);
+    
+    CCTK_REAL S3FluxL = 0.5*(S3FluxLeft + S3FluxRight + (-S3LeftL + 
+      ShiftMinus3S3Right)*ToReal(hlleAlpha));
+    
+    CCTK_REAL EnFluxLeft = (EnLeftL + pLeftL)*v3LeftL;
+    
+    CCTK_REAL EnFluxRight = (ShiftMinus3EnRight + 
+      ShiftMinus3pRight)*ShiftMinus3v3Right;
+    
+    CCTK_REAL EnFluxL = 0.5*(EnFluxLeft + EnFluxRight + (-EnLeftL + 
+      ShiftMinus3EnRight)*ToReal(hlleAlpha));
+    
+    /* Copy local copies back to grid functions */
+    DenFlux[index] = DenFluxL;
+    EnFlux[index] = EnFluxL;
+    S1Flux[index] = S1FluxL;
+    S2Flux[index] = S2FluxL;
+    S3Flux[index] = S3FluxL;
   }
+  CCTK_ENDLOOP3 (eulerauto_cons_calc_flux_3);
 }
 
 extern "C" void eulerauto_cons_calc_flux_3(CCTK_ARGUMENTS)
@@ -202,5 +195,26 @@ extern "C" void eulerauto_cons_calc_flux_3(CCTK_ARGUMENTS)
   DECLARE_CCTK_ARGUMENTS;
   DECLARE_CCTK_PARAMETERS;
   
+  
+  if (verbose > 1)
+  {
+    CCTK_VInfo(CCTK_THORNSTRING,"Entering eulerauto_cons_calc_flux_3_Body");
+  }
+  
+  if (cctk_iteration % eulerauto_cons_calc_flux_3_calc_every != eulerauto_cons_calc_flux_3_calc_offset)
+  {
+    return;
+  }
+  
+  const char *groups[] = {"EulerAuto::Den_flux_group","EulerAuto::Den_lr_group","EulerAuto::En_flux_group","EulerAuto::En_lr_group","EulerAuto::p_lr_group","EulerAuto::rho_lr_group","EulerAuto::S1_flux_group","EulerAuto::S1_lr_group","EulerAuto::S2_flux_group","EulerAuto::S2_lr_group","EulerAuto::S3_flux_group","EulerAuto::S3_lr_group","EulerAuto::v1_lr_group","EulerAuto::v2_lr_group","EulerAuto::v3_lr_group"};
+  GenericFD_AssertGroupStorage(cctkGH, "eulerauto_cons_calc_flux_3", 15, groups);
+  
+  GenericFD_EnsureStencilFits(cctkGH, "eulerauto_cons_calc_flux_3", 1, 1, 1);
+  
   GenericFD_LoopOverInterior(cctkGH, &eulerauto_cons_calc_flux_3_Body);
+  
+  if (verbose > 1)
+  {
+    CCTK_VInfo(CCTK_THORNSTRING,"Leaving eulerauto_cons_calc_flux_3_Body");
+  }
 }
diff --git a/Examples/EulerAuto/src/eulerauto_cons_calc_intercell_conserved_1.cc b/Examples/EulerAuto/src/eulerauto_cons_calc_intercell_conserved_1.cc
index 3f91fa6..fdcd6f7 100644
--- a/Examples/EulerAuto/src/eulerauto_cons_calc_intercell_conserved_1.cc
+++ b/Examples/EulerAuto/src/eulerauto_cons_calc_intercell_conserved_1.cc
@@ -12,6 +12,8 @@
 #include "cctk_Parameters.h"
 #include "GenericFD.h"
 #include "Differencing.h"
+#include "cctk_Loop.h"
+#include "loopcontrol.h"
 
 /* Define macros used in calculations */
 #define INITVALUE (42)
@@ -20,42 +22,28 @@
 #define SQR(x) ((x) * (x))
 #define CUB(x) ((x) * (x) * (x))
 
-static void eulerauto_cons_calc_intercell_conserved_1_Body(cGH const * restrict const cctkGH, int const dir, int const face, CCTK_REAL const normal[3], CCTK_REAL const tangentA[3], CCTK_REAL const tangentB[3], int const min[3], int const max[3], int const n_subblock_gfs, CCTK_REAL * restrict const subblock_gfs[])
+static void eulerauto_cons_calc_intercell_conserved_1_Body(cGH const * restrict const cctkGH, int const dir, int const face, CCTK_REAL const normal[3], CCTK_REAL const tangentA[3], CCTK_REAL const tangentB[3], int const imin[3], int const imax[3], int const n_subblock_gfs, CCTK_REAL * restrict const subblock_gfs[])
 {
   DECLARE_CCTK_ARGUMENTS;
   DECLARE_CCTK_PARAMETERS;
   
   
-  /* Declare the variables used for looping over grid points */
-  CCTK_INT i, j, k;
-  // CCTK_INT index = INITVALUE;
-  
   /* Declare finite differencing variables */
   
-  if (verbose > 1)
-  {
-    CCTK_VInfo(CCTK_THORNSTRING,"Entering eulerauto_cons_calc_intercell_conserved_1_Body");
-  }
-  
-  if (cctk_iteration % eulerauto_cons_calc_intercell_conserved_1_calc_every != eulerauto_cons_calc_intercell_conserved_1_calc_offset)
-  {
-    return;
-  }
-  
-  const char *groups[] = {"EulerAuto::Den_lr_group","EulerAuto::En_lr_group","EulerAuto::p_lr_group","EulerAuto::rho_lr_group","EulerAuto::S1_lr_group","EulerAuto::S2_lr_group","EulerAuto::S3_lr_group","EulerAuto::v1_lr_group","EulerAuto::v2_lr_group","EulerAuto::v3_lr_group"};
-  GenericFD_AssertGroupStorage(cctkGH, "eulerauto_cons_calc_intercell_conserved_1", 10, groups);
-  
-  
   /* Include user-supplied include files */
   
   /* Initialise finite differencing variables */
   ptrdiff_t const di = 1;
   ptrdiff_t const dj = CCTK_GFINDEX3D(cctkGH,0,1,0) - CCTK_GFINDEX3D(cctkGH,0,0,0);
   ptrdiff_t const dk = CCTK_GFINDEX3D(cctkGH,0,0,1) - CCTK_GFINDEX3D(cctkGH,0,0,0);
+  ptrdiff_t const cdi = sizeof(CCTK_REAL) * di;
+  ptrdiff_t const cdj = sizeof(CCTK_REAL) * dj;
+  ptrdiff_t const cdk = sizeof(CCTK_REAL) * dk;
   CCTK_REAL const dx = ToReal(CCTK_DELTA_SPACE(0));
   CCTK_REAL const dy = ToReal(CCTK_DELTA_SPACE(1));
   CCTK_REAL const dz = ToReal(CCTK_DELTA_SPACE(2));
   CCTK_REAL const dt = ToReal(CCTK_DELTA_TIME);
+  CCTK_REAL const t = ToReal(cctk_time);
   CCTK_REAL const dxi = INV(dx);
   CCTK_REAL const dyi = INV(dy);
   CCTK_REAL const dzi = INV(dz);
@@ -74,70 +62,78 @@ static void eulerauto_cons_calc_intercell_conserved_1_Body(cGH const * restrict
   CCTK_REAL const p1ody = INV(dy);
   CCTK_REAL const p1odz = INV(dz);
   
+  /* Assign local copies of arrays functions */
+  
+  
+  
+  /* Calculate temporaries and arrays functions */
+  
+  /* Copy local copies back to grid functions */
+  
   /* Loop over the grid points */
-  for (k = min[2]; k < max[2]; k++)
+  #pragma omp parallel
+  CCTK_LOOP3 (eulerauto_cons_calc_intercell_conserved_1,
+    i,j,k, imin[0],imin[1],imin[2], imax[0],imax[1],imax[2],
+    cctk_lsh[0],cctk_lsh[1],cctk_lsh[2])
   {
-    for (j = min[1]; j < max[1]; j++)
-    {
-      for (i = min[0]; i < max[0]; i++)
-      {
-         int  const  index  =  CCTK_GFINDEX3D(cctkGH,i,j,k) ;
-        
-        /* Assign local copies of grid functions */
-        
-        CCTK_REAL pLeftL = pLeft[index];
-        CCTK_REAL pRightL = pRight[index];
-        CCTK_REAL rhoLeftL = rhoLeft[index];
-        CCTK_REAL rhoRightL = rhoRight[index];
-        CCTK_REAL v1LeftL = v1Left[index];
-        CCTK_REAL v1RightL = v1Right[index];
-        CCTK_REAL v2LeftL = v2Left[index];
-        CCTK_REAL v2RightL = v2Right[index];
-        CCTK_REAL v3LeftL = v3Left[index];
-        CCTK_REAL v3RightL = v3Right[index];
-        
-        
-        /* Include user supplied include files */
-        
-        /* Precompute derivatives */
-        
-        /* Calculate temporaries and grid functions */
-        CCTK_REAL DenLeftL = rhoLeftL;
-        
-        CCTK_REAL S1LeftL = rhoLeftL*v1LeftL;
-        
-        CCTK_REAL S2LeftL = rhoLeftL*v2LeftL;
-        
-        CCTK_REAL S3LeftL = rhoLeftL*v3LeftL;
-        
-        CCTK_REAL EnLeftL = pLeftL*INV(-1 + ToReal(gamma)) + 
-          0.5*rhoLeftL*(SQR(v1LeftL) + SQR(v2LeftL) + SQR(v3LeftL));
-        
-        CCTK_REAL DenRightL = rhoRightL;
-        
-        CCTK_REAL S1RightL = rhoRightL*v1RightL;
-        
-        CCTK_REAL S2RightL = rhoRightL*v2RightL;
-        
-        CCTK_REAL S3RightL = rhoRightL*v3RightL;
-        
-        CCTK_REAL EnRightL = pRightL*INV(-1 + ToReal(gamma)) + 
-          0.5*rhoRightL*(SQR(v1RightL) + SQR(v2RightL) + SQR(v3RightL));
-        
-        /* Copy local copies back to grid functions */
-        DenLeft[index] = DenLeftL;
-        DenRight[index] = DenRightL;
-        EnLeft[index] = EnLeftL;
-        EnRight[index] = EnRightL;
-        S1Left[index] = S1LeftL;
-        S1Right[index] = S1RightL;
-        S2Left[index] = S2LeftL;
-        S2Right[index] = S2RightL;
-        S3Left[index] = S3LeftL;
-        S3Right[index] = S3RightL;
-      }
-    }
+    ptrdiff_t const index = di*i + dj*j + dk*k;
+    
+    /* Assign local copies of grid functions */
+    
+    CCTK_REAL pLeftL = pLeft[index];
+    CCTK_REAL pRightL = pRight[index];
+    CCTK_REAL rhoLeftL = rhoLeft[index];
+    CCTK_REAL rhoRightL = rhoRight[index];
+    CCTK_REAL v1LeftL = v1Left[index];
+    CCTK_REAL v1RightL = v1Right[index];
+    CCTK_REAL v2LeftL = v2Left[index];
+    CCTK_REAL v2RightL = v2Right[index];
+    CCTK_REAL v3LeftL = v3Left[index];
+    CCTK_REAL v3RightL = v3Right[index];
+    
+    
+    /* Include user supplied include files */
+    
+    /* Precompute derivatives */
+    
+    /* Calculate temporaries and grid functions */
+    CCTK_REAL DenLeftL = rhoLeftL;
+    
+    CCTK_REAL S1LeftL = rhoLeftL*v1LeftL;
+    
+    CCTK_REAL S2LeftL = rhoLeftL*v2LeftL;
+    
+    CCTK_REAL S3LeftL = rhoLeftL*v3LeftL;
+    
+    CCTK_REAL EnLeftL = 0.5*INV(-1 + ToReal(gamma))*(2*pLeftL + 
+      rhoLeftL*(SQR(v1LeftL) + SQR(v2LeftL) + SQR(v3LeftL))*(-1 + 
+      ToReal(gamma)));
+    
+    CCTK_REAL DenRightL = rhoRightL;
+    
+    CCTK_REAL S1RightL = rhoRightL*v1RightL;
+    
+    CCTK_REAL S2RightL = rhoRightL*v2RightL;
+    
+    CCTK_REAL S3RightL = rhoRightL*v3RightL;
+    
+    CCTK_REAL EnRightL = 0.5*INV(-1 + ToReal(gamma))*(2*pRightL + 
+      rhoRightL*(SQR(v1RightL) + SQR(v2RightL) + SQR(v3RightL))*(-1 + 
+      ToReal(gamma)));
+    
+    /* Copy local copies back to grid functions */
+    DenLeft[index] = DenLeftL;
+    DenRight[index] = DenRightL;
+    EnLeft[index] = EnLeftL;
+    EnRight[index] = EnRightL;
+    S1Left[index] = S1LeftL;
+    S1Right[index] = S1RightL;
+    S2Left[index] = S2LeftL;
+    S2Right[index] = S2RightL;
+    S3Left[index] = S3LeftL;
+    S3Right[index] = S3RightL;
   }
+  CCTK_ENDLOOP3 (eulerauto_cons_calc_intercell_conserved_1);
 }
 
 extern "C" void eulerauto_cons_calc_intercell_conserved_1(CCTK_ARGUMENTS)
@@ -145,5 +141,25 @@ extern "C" void eulerauto_cons_calc_intercell_conserved_1(CCTK_ARGUMENTS)
   DECLARE_CCTK_ARGUMENTS;
   DECLARE_CCTK_PARAMETERS;
   
+  
+  if (verbose > 1)
+  {
+    CCTK_VInfo(CCTK_THORNSTRING,"Entering eulerauto_cons_calc_intercell_conserved_1_Body");
+  }
+  
+  if (cctk_iteration % eulerauto_cons_calc_intercell_conserved_1_calc_every != eulerauto_cons_calc_intercell_conserved_1_calc_offset)
+  {
+    return;
+  }
+  
+  const char *groups[] = {"EulerAuto::Den_lr_group","EulerAuto::En_lr_group","EulerAuto::p_lr_group","EulerAuto::rho_lr_group","EulerAuto::S1_lr_group","EulerAuto::S2_lr_group","EulerAuto::S3_lr_group","EulerAuto::v1_lr_group","EulerAuto::v2_lr_group","EulerAuto::v3_lr_group"};
+  GenericFD_AssertGroupStorage(cctkGH, "eulerauto_cons_calc_intercell_conserved_1", 10, groups);
+  
+  
   GenericFD_LoopOverEverything(cctkGH, &eulerauto_cons_calc_intercell_conserved_1_Body);
+  
+  if (verbose > 1)
+  {
+    CCTK_VInfo(CCTK_THORNSTRING,"Leaving eulerauto_cons_calc_intercell_conserved_1_Body");
+  }
 }
diff --git a/Examples/EulerAuto/src/eulerauto_cons_calc_intercell_conserved_2.cc b/Examples/EulerAuto/src/eulerauto_cons_calc_intercell_conserved_2.cc
index 209949e..9a8709d 100644
--- a/Examples/EulerAuto/src/eulerauto_cons_calc_intercell_conserved_2.cc
+++ b/Examples/EulerAuto/src/eulerauto_cons_calc_intercell_conserved_2.cc
@@ -12,6 +12,8 @@
 #include "cctk_Parameters.h"
 #include "GenericFD.h"
 #include "Differencing.h"
+#include "cctk_Loop.h"
+#include "loopcontrol.h"
 
 /* Define macros used in calculations */
 #define INITVALUE (42)
@@ -20,42 +22,28 @@
 #define SQR(x) ((x) * (x))
 #define CUB(x) ((x) * (x) * (x))
 
-static void eulerauto_cons_calc_intercell_conserved_2_Body(cGH const * restrict const cctkGH, int const dir, int const face, CCTK_REAL const normal[3], CCTK_REAL const tangentA[3], CCTK_REAL const tangentB[3], int const min[3], int const max[3], int const n_subblock_gfs, CCTK_REAL * restrict const subblock_gfs[])
+static void eulerauto_cons_calc_intercell_conserved_2_Body(cGH const * restrict const cctkGH, int const dir, int const face, CCTK_REAL const normal[3], CCTK_REAL const tangentA[3], CCTK_REAL const tangentB[3], int const imin[3], int const imax[3], int const n_subblock_gfs, CCTK_REAL * restrict const subblock_gfs[])
 {
   DECLARE_CCTK_ARGUMENTS;
   DECLARE_CCTK_PARAMETERS;
   
   
-  /* Declare the variables used for looping over grid points */
-  CCTK_INT i, j, k;
-  // CCTK_INT index = INITVALUE;
-  
   /* Declare finite differencing variables */
   
-  if (verbose > 1)
-  {
-    CCTK_VInfo(CCTK_THORNSTRING,"Entering eulerauto_cons_calc_intercell_conserved_2_Body");
-  }
-  
-  if (cctk_iteration % eulerauto_cons_calc_intercell_conserved_2_calc_every != eulerauto_cons_calc_intercell_conserved_2_calc_offset)
-  {
-    return;
-  }
-  
-  const char *groups[] = {"EulerAuto::Den_lr_group","EulerAuto::En_lr_group","EulerAuto::p_lr_group","EulerAuto::rho_lr_group","EulerAuto::S1_lr_group","EulerAuto::S2_lr_group","EulerAuto::S3_lr_group","EulerAuto::v1_lr_group","EulerAuto::v2_lr_group","EulerAuto::v3_lr_group"};
-  GenericFD_AssertGroupStorage(cctkGH, "eulerauto_cons_calc_intercell_conserved_2", 10, groups);
-  
-  
   /* Include user-supplied include files */
   
   /* Initialise finite differencing variables */
   ptrdiff_t const di = 1;
   ptrdiff_t const dj = CCTK_GFINDEX3D(cctkGH,0,1,0) - CCTK_GFINDEX3D(cctkGH,0,0,0);
   ptrdiff_t const dk = CCTK_GFINDEX3D(cctkGH,0,0,1) - CCTK_GFINDEX3D(cctkGH,0,0,0);
+  ptrdiff_t const cdi = sizeof(CCTK_REAL) * di;
+  ptrdiff_t const cdj = sizeof(CCTK_REAL) * dj;
+  ptrdiff_t const cdk = sizeof(CCTK_REAL) * dk;
   CCTK_REAL const dx = ToReal(CCTK_DELTA_SPACE(0));
   CCTK_REAL const dy = ToReal(CCTK_DELTA_SPACE(1));
   CCTK_REAL const dz = ToReal(CCTK_DELTA_SPACE(2));
   CCTK_REAL const dt = ToReal(CCTK_DELTA_TIME);
+  CCTK_REAL const t = ToReal(cctk_time);
   CCTK_REAL const dxi = INV(dx);
   CCTK_REAL const dyi = INV(dy);
   CCTK_REAL const dzi = INV(dz);
@@ -74,70 +62,78 @@ static void eulerauto_cons_calc_intercell_conserved_2_Body(cGH const * restrict
   CCTK_REAL const p1ody = INV(dy);
   CCTK_REAL const p1odz = INV(dz);
   
+  /* Assign local copies of arrays functions */
+  
+  
+  
+  /* Calculate temporaries and arrays functions */
+  
+  /* Copy local copies back to grid functions */
+  
   /* Loop over the grid points */
-  for (k = min[2]; k < max[2]; k++)
+  #pragma omp parallel
+  CCTK_LOOP3 (eulerauto_cons_calc_intercell_conserved_2,
+    i,j,k, imin[0],imin[1],imin[2], imax[0],imax[1],imax[2],
+    cctk_lsh[0],cctk_lsh[1],cctk_lsh[2])
   {
-    for (j = min[1]; j < max[1]; j++)
-    {
-      for (i = min[0]; i < max[0]; i++)
-      {
-         int  const  index  =  CCTK_GFINDEX3D(cctkGH,i,j,k) ;
-        
-        /* Assign local copies of grid functions */
-        
-        CCTK_REAL pLeftL = pLeft[index];
-        CCTK_REAL pRightL = pRight[index];
-        CCTK_REAL rhoLeftL = rhoLeft[index];
-        CCTK_REAL rhoRightL = rhoRight[index];
-        CCTK_REAL v1LeftL = v1Left[index];
-        CCTK_REAL v1RightL = v1Right[index];
-        CCTK_REAL v2LeftL = v2Left[index];
-        CCTK_REAL v2RightL = v2Right[index];
-        CCTK_REAL v3LeftL = v3Left[index];
-        CCTK_REAL v3RightL = v3Right[index];
-        
-        
-        /* Include user supplied include files */
-        
-        /* Precompute derivatives */
-        
-        /* Calculate temporaries and grid functions */
-        CCTK_REAL DenLeftL = rhoLeftL;
-        
-        CCTK_REAL S1LeftL = rhoLeftL*v1LeftL;
-        
-        CCTK_REAL S2LeftL = rhoLeftL*v2LeftL;
-        
-        CCTK_REAL S3LeftL = rhoLeftL*v3LeftL;
-        
-        CCTK_REAL EnLeftL = pLeftL*INV(-1 + ToReal(gamma)) + 
-          0.5*rhoLeftL*(SQR(v1LeftL) + SQR(v2LeftL) + SQR(v3LeftL));
-        
-        CCTK_REAL DenRightL = rhoRightL;
-        
-        CCTK_REAL S1RightL = rhoRightL*v1RightL;
-        
-        CCTK_REAL S2RightL = rhoRightL*v2RightL;
-        
-        CCTK_REAL S3RightL = rhoRightL*v3RightL;
-        
-        CCTK_REAL EnRightL = pRightL*INV(-1 + ToReal(gamma)) + 
-          0.5*rhoRightL*(SQR(v1RightL) + SQR(v2RightL) + SQR(v3RightL));
-        
-        /* Copy local copies back to grid functions */
-        DenLeft[index] = DenLeftL;
-        DenRight[index] = DenRightL;
-        EnLeft[index] = EnLeftL;
-        EnRight[index] = EnRightL;
-        S1Left[index] = S1LeftL;
-        S1Right[index] = S1RightL;
-        S2Left[index] = S2LeftL;
-        S2Right[index] = S2RightL;
-        S3Left[index] = S3LeftL;
-        S3Right[index] = S3RightL;
-      }
-    }
+    ptrdiff_t const index = di*i + dj*j + dk*k;
+    
+    /* Assign local copies of grid functions */
+    
+    CCTK_REAL pLeftL = pLeft[index];
+    CCTK_REAL pRightL = pRight[index];
+    CCTK_REAL rhoLeftL = rhoLeft[index];
+    CCTK_REAL rhoRightL = rhoRight[index];
+    CCTK_REAL v1LeftL = v1Left[index];
+    CCTK_REAL v1RightL = v1Right[index];
+    CCTK_REAL v2LeftL = v2Left[index];
+    CCTK_REAL v2RightL = v2Right[index];
+    CCTK_REAL v3LeftL = v3Left[index];
+    CCTK_REAL v3RightL = v3Right[index];
+    
+    
+    /* Include user supplied include files */
+    
+    /* Precompute derivatives */
+    
+    /* Calculate temporaries and grid functions */
+    CCTK_REAL DenLeftL = rhoLeftL;
+    
+    CCTK_REAL S1LeftL = rhoLeftL*v1LeftL;
+    
+    CCTK_REAL S2LeftL = rhoLeftL*v2LeftL;
+    
+    CCTK_REAL S3LeftL = rhoLeftL*v3LeftL;
+    
+    CCTK_REAL EnLeftL = 0.5*INV(-1 + ToReal(gamma))*(2*pLeftL + 
+      rhoLeftL*(SQR(v1LeftL) + SQR(v2LeftL) + SQR(v3LeftL))*(-1 + 
+      ToReal(gamma)));
+    
+    CCTK_REAL DenRightL = rhoRightL;
+    
+    CCTK_REAL S1RightL = rhoRightL*v1RightL;
+    
+    CCTK_REAL S2RightL = rhoRightL*v2RightL;
+    
+    CCTK_REAL S3RightL = rhoRightL*v3RightL;
+    
+    CCTK_REAL EnRightL = 0.5*INV(-1 + ToReal(gamma))*(2*pRightL + 
+      rhoRightL*(SQR(v1RightL) + SQR(v2RightL) + SQR(v3RightL))*(-1 + 
+      ToReal(gamma)));
+    
+    /* Copy local copies back to grid functions */
+    DenLeft[index] = DenLeftL;
+    DenRight[index] = DenRightL;
+    EnLeft[index] = EnLeftL;
+    EnRight[index] = EnRightL;
+    S1Left[index] = S1LeftL;
+    S1Right[index] = S1RightL;
+    S2Left[index] = S2LeftL;
+    S2Right[index] = S2RightL;
+    S3Left[index] = S3LeftL;
+    S3Right[index] = S3RightL;
   }
+  CCTK_ENDLOOP3 (eulerauto_cons_calc_intercell_conserved_2);
 }
 
 extern "C" void eulerauto_cons_calc_intercell_conserved_2(CCTK_ARGUMENTS)
@@ -145,5 +141,25 @@ extern "C" void eulerauto_cons_calc_intercell_conserved_2(CCTK_ARGUMENTS)
   DECLARE_CCTK_ARGUMENTS;
   DECLARE_CCTK_PARAMETERS;
   
+  
+  if (verbose > 1)
+  {
+    CCTK_VInfo(CCTK_THORNSTRING,"Entering eulerauto_cons_calc_intercell_conserved_2_Body");
+  }
+  
+  if (cctk_iteration % eulerauto_cons_calc_intercell_conserved_2_calc_every != eulerauto_cons_calc_intercell_conserved_2_calc_offset)
+  {
+    return;
+  }
+  
+  const char *groups[] = {"EulerAuto::Den_lr_group","EulerAuto::En_lr_group","EulerAuto::p_lr_group","EulerAuto::rho_lr_group","EulerAuto::S1_lr_group","EulerAuto::S2_lr_group","EulerAuto::S3_lr_group","EulerAuto::v1_lr_group","EulerAuto::v2_lr_group","EulerAuto::v3_lr_group"};
+  GenericFD_AssertGroupStorage(cctkGH, "eulerauto_cons_calc_intercell_conserved_2", 10, groups);
+  
+  
   GenericFD_LoopOverEverything(cctkGH, &eulerauto_cons_calc_intercell_conserved_2_Body);
+  
+  if (verbose > 1)
+  {
+    CCTK_VInfo(CCTK_THORNSTRING,"Leaving eulerauto_cons_calc_intercell_conserved_2_Body");
+  }
 }
diff --git a/Examples/EulerAuto/src/eulerauto_cons_calc_intercell_conserved_3.cc b/Examples/EulerAuto/src/eulerauto_cons_calc_intercell_conserved_3.cc
index 9398982..2c4f3a1 100644
--- a/Examples/EulerAuto/src/eulerauto_cons_calc_intercell_conserved_3.cc
+++ b/Examples/EulerAuto/src/eulerauto_cons_calc_intercell_conserved_3.cc
@@ -12,6 +12,8 @@
 #include "cctk_Parameters.h"
 #include "GenericFD.h"
 #include "Differencing.h"
+#include "cctk_Loop.h"
+#include "loopcontrol.h"
 
 /* Define macros used in calculations */
 #define INITVALUE (42)
@@ -20,42 +22,28 @@
 #define SQR(x) ((x) * (x))
 #define CUB(x) ((x) * (x) * (x))
 
-static void eulerauto_cons_calc_intercell_conserved_3_Body(cGH const * restrict const cctkGH, int const dir, int const face, CCTK_REAL const normal[3], CCTK_REAL const tangentA[3], CCTK_REAL const tangentB[3], int const min[3], int const max[3], int const n_subblock_gfs, CCTK_REAL * restrict const subblock_gfs[])
+static void eulerauto_cons_calc_intercell_conserved_3_Body(cGH const * restrict const cctkGH, int const dir, int const face, CCTK_REAL const normal[3], CCTK_REAL const tangentA[3], CCTK_REAL const tangentB[3], int const imin[3], int const imax[3], int const n_subblock_gfs, CCTK_REAL * restrict const subblock_gfs[])
 {
   DECLARE_CCTK_ARGUMENTS;
   DECLARE_CCTK_PARAMETERS;
   
   
-  /* Declare the variables used for looping over grid points */
-  CCTK_INT i, j, k;
-  // CCTK_INT index = INITVALUE;
-  
   /* Declare finite differencing variables */
   
-  if (verbose > 1)
-  {
-    CCTK_VInfo(CCTK_THORNSTRING,"Entering eulerauto_cons_calc_intercell_conserved_3_Body");
-  }
-  
-  if (cctk_iteration % eulerauto_cons_calc_intercell_conserved_3_calc_every != eulerauto_cons_calc_intercell_conserved_3_calc_offset)
-  {
-    return;
-  }
-  
-  const char *groups[] = {"EulerAuto::Den_lr_group","EulerAuto::En_lr_group","EulerAuto::p_lr_group","EulerAuto::rho_lr_group","EulerAuto::S1_lr_group","EulerAuto::S2_lr_group","EulerAuto::S3_lr_group","EulerAuto::v1_lr_group","EulerAuto::v2_lr_group","EulerAuto::v3_lr_group"};
-  GenericFD_AssertGroupStorage(cctkGH, "eulerauto_cons_calc_intercell_conserved_3", 10, groups);
-  
-  
   /* Include user-supplied include files */
   
   /* Initialise finite differencing variables */
   ptrdiff_t const di = 1;
   ptrdiff_t const dj = CCTK_GFINDEX3D(cctkGH,0,1,0) - CCTK_GFINDEX3D(cctkGH,0,0,0);
   ptrdiff_t const dk = CCTK_GFINDEX3D(cctkGH,0,0,1) - CCTK_GFINDEX3D(cctkGH,0,0,0);
+  ptrdiff_t const cdi = sizeof(CCTK_REAL) * di;
+  ptrdiff_t const cdj = sizeof(CCTK_REAL) * dj;
+  ptrdiff_t const cdk = sizeof(CCTK_REAL) * dk;
   CCTK_REAL const dx = ToReal(CCTK_DELTA_SPACE(0));
   CCTK_REAL const dy = ToReal(CCTK_DELTA_SPACE(1));
   CCTK_REAL const dz = ToReal(CCTK_DELTA_SPACE(2));
   CCTK_REAL const dt = ToReal(CCTK_DELTA_TIME);
+  CCTK_REAL const t = ToReal(cctk_time);
   CCTK_REAL const dxi = INV(dx);
   CCTK_REAL const dyi = INV(dy);
   CCTK_REAL const dzi = INV(dz);
@@ -74,70 +62,78 @@ static void eulerauto_cons_calc_intercell_conserved_3_Body(cGH const * restrict
   CCTK_REAL const p1ody = INV(dy);
   CCTK_REAL const p1odz = INV(dz);
   
+  /* Assign local copies of arrays functions */
+  
+  
+  
+  /* Calculate temporaries and arrays functions */
+  
+  /* Copy local copies back to grid functions */
+  
   /* Loop over the grid points */
-  for (k = min[2]; k < max[2]; k++)
+  #pragma omp parallel
+  CCTK_LOOP3 (eulerauto_cons_calc_intercell_conserved_3,
+    i,j,k, imin[0],imin[1],imin[2], imax[0],imax[1],imax[2],
+    cctk_lsh[0],cctk_lsh[1],cctk_lsh[2])
   {
-    for (j = min[1]; j < max[1]; j++)
-    {
-      for (i = min[0]; i < max[0]; i++)
-      {
-         int  const  index  =  CCTK_GFINDEX3D(cctkGH,i,j,k) ;
-        
-        /* Assign local copies of grid functions */
-        
-        CCTK_REAL pLeftL = pLeft[index];
-        CCTK_REAL pRightL = pRight[index];
-        CCTK_REAL rhoLeftL = rhoLeft[index];
-        CCTK_REAL rhoRightL = rhoRight[index];
-        CCTK_REAL v1LeftL = v1Left[index];
-        CCTK_REAL v1RightL = v1Right[index];
-        CCTK_REAL v2LeftL = v2Left[index];
-        CCTK_REAL v2RightL = v2Right[index];
-        CCTK_REAL v3LeftL = v3Left[index];
-        CCTK_REAL v3RightL = v3Right[index];
-        
-        
-        /* Include user supplied include files */
-        
-        /* Precompute derivatives */
-        
-        /* Calculate temporaries and grid functions */
-        CCTK_REAL DenLeftL = rhoLeftL;
-        
-        CCTK_REAL S1LeftL = rhoLeftL*v1LeftL;
-        
-        CCTK_REAL S2LeftL = rhoLeftL*v2LeftL;
-        
-        CCTK_REAL S3LeftL = rhoLeftL*v3LeftL;
-        
-        CCTK_REAL EnLeftL = pLeftL*INV(-1 + ToReal(gamma)) + 
-          0.5*rhoLeftL*(SQR(v1LeftL) + SQR(v2LeftL) + SQR(v3LeftL));
-        
-        CCTK_REAL DenRightL = rhoRightL;
-        
-        CCTK_REAL S1RightL = rhoRightL*v1RightL;
-        
-        CCTK_REAL S2RightL = rhoRightL*v2RightL;
-        
-        CCTK_REAL S3RightL = rhoRightL*v3RightL;
-        
-        CCTK_REAL EnRightL = pRightL*INV(-1 + ToReal(gamma)) + 
-          0.5*rhoRightL*(SQR(v1RightL) + SQR(v2RightL) + SQR(v3RightL));
-        
-        /* Copy local copies back to grid functions */
-        DenLeft[index] = DenLeftL;
-        DenRight[index] = DenRightL;
-        EnLeft[index] = EnLeftL;
-        EnRight[index] = EnRightL;
-        S1Left[index] = S1LeftL;
-        S1Right[index] = S1RightL;
-        S2Left[index] = S2LeftL;
-        S2Right[index] = S2RightL;
-        S3Left[index] = S3LeftL;
-        S3Right[index] = S3RightL;
-      }
-    }
+    ptrdiff_t const index = di*i + dj*j + dk*k;
+    
+    /* Assign local copies of grid functions */
+    
+    CCTK_REAL pLeftL = pLeft[index];
+    CCTK_REAL pRightL = pRight[index];
+    CCTK_REAL rhoLeftL = rhoLeft[index];
+    CCTK_REAL rhoRightL = rhoRight[index];
+    CCTK_REAL v1LeftL = v1Left[index];
+    CCTK_REAL v1RightL = v1Right[index];
+    CCTK_REAL v2LeftL = v2Left[index];
+    CCTK_REAL v2RightL = v2Right[index];
+    CCTK_REAL v3LeftL = v3Left[index];
+    CCTK_REAL v3RightL = v3Right[index];
+    
+    
+    /* Include user supplied include files */
+    
+    /* Precompute derivatives */
+    
+    /* Calculate temporaries and grid functions */
+    CCTK_REAL DenLeftL = rhoLeftL;
+    
+    CCTK_REAL S1LeftL = rhoLeftL*v1LeftL;
+    
+    CCTK_REAL S2LeftL = rhoLeftL*v2LeftL;
+    
+    CCTK_REAL S3LeftL = rhoLeftL*v3LeftL;
+    
+    CCTK_REAL EnLeftL = 0.5*INV(-1 + ToReal(gamma))*(2*pLeftL + 
+      rhoLeftL*(SQR(v1LeftL) + SQR(v2LeftL) + SQR(v3LeftL))*(-1 + 
+      ToReal(gamma)));
+    
+    CCTK_REAL DenRightL = rhoRightL;
+    
+    CCTK_REAL S1RightL = rhoRightL*v1RightL;
+    
+    CCTK_REAL S2RightL = rhoRightL*v2RightL;
+    
+    CCTK_REAL S3RightL = rhoRightL*v3RightL;
+    
+    CCTK_REAL EnRightL = 0.5*INV(-1 + ToReal(gamma))*(2*pRightL + 
+      rhoRightL*(SQR(v1RightL) + SQR(v2RightL) + SQR(v3RightL))*(-1 + 
+      ToReal(gamma)));
+    
+    /* Copy local copies back to grid functions */
+    DenLeft[index] = DenLeftL;
+    DenRight[index] = DenRightL;
+    EnLeft[index] = EnLeftL;
+    EnRight[index] = EnRightL;
+    S1Left[index] = S1LeftL;
+    S1Right[index] = S1RightL;
+    S2Left[index] = S2LeftL;
+    S2Right[index] = S2RightL;
+    S3Left[index] = S3LeftL;
+    S3Right[index] = S3RightL;
   }
+  CCTK_ENDLOOP3 (eulerauto_cons_calc_intercell_conserved_3);
 }
 
 extern "C" void eulerauto_cons_calc_intercell_conserved_3(CCTK_ARGUMENTS)
@@ -145,5 +141,25 @@ extern "C" void eulerauto_cons_calc_intercell_conserved_3(CCTK_ARGUMENTS)
   DECLARE_CCTK_ARGUMENTS;
   DECLARE_CCTK_PARAMETERS;
   
+  
+  if (verbose > 1)
+  {
+    CCTK_VInfo(CCTK_THORNSTRING,"Entering eulerauto_cons_calc_intercell_conserved_3_Body");
+  }
+  
+  if (cctk_iteration % eulerauto_cons_calc_intercell_conserved_3_calc_every != eulerauto_cons_calc_intercell_conserved_3_calc_offset)
+  {
+    return;
+  }
+  
+  const char *groups[] = {"EulerAuto::Den_lr_group","EulerAuto::En_lr_group","EulerAuto::p_lr_group","EulerAuto::rho_lr_group","EulerAuto::S1_lr_group","EulerAuto::S2_lr_group","EulerAuto::S3_lr_group","EulerAuto::v1_lr_group","EulerAuto::v2_lr_group","EulerAuto::v3_lr_group"};
+  GenericFD_AssertGroupStorage(cctkGH, "eulerauto_cons_calc_intercell_conserved_3", 10, groups);
+  
+  
   GenericFD_LoopOverEverything(cctkGH, &eulerauto_cons_calc_intercell_conserved_3_Body);
+  
+  if (verbose > 1)
+  {
+    CCTK_VInfo(CCTK_THORNSTRING,"Leaving eulerauto_cons_calc_intercell_conserved_3_Body");
+  }
 }
diff --git a/Examples/EulerAuto/src/eulerauto_cons_calc_primitives.cc b/Examples/EulerAuto/src/eulerauto_cons_calc_primitives.cc
index 0a7228a..1b16069 100644
--- a/Examples/EulerAuto/src/eulerauto_cons_calc_primitives.cc
+++ b/Examples/EulerAuto/src/eulerauto_cons_calc_primitives.cc
@@ -12,6 +12,8 @@
 #include "cctk_Parameters.h"
 #include "GenericFD.h"
 #include "Differencing.h"
+#include "cctk_Loop.h"
+#include "loopcontrol.h"
 
 /* Define macros used in calculations */
 #define INITVALUE (42)
@@ -20,42 +22,28 @@
 #define SQR(x) ((x) * (x))
 #define CUB(x) ((x) * (x) * (x))
 
-static void eulerauto_cons_calc_primitives_Body(cGH const * restrict const cctkGH, int const dir, int const face, CCTK_REAL const normal[3], CCTK_REAL const tangentA[3], CCTK_REAL const tangentB[3], int const min[3], int const max[3], int const n_subblock_gfs, CCTK_REAL * restrict const subblock_gfs[])
+static void eulerauto_cons_calc_primitives_Body(cGH const * restrict const cctkGH, int const dir, int const face, CCTK_REAL const normal[3], CCTK_REAL const tangentA[3], CCTK_REAL const tangentB[3], int const imin[3], int const imax[3], int const n_subblock_gfs, CCTK_REAL * restrict const subblock_gfs[])
 {
   DECLARE_CCTK_ARGUMENTS;
   DECLARE_CCTK_PARAMETERS;
   
   
-  /* Declare the variables used for looping over grid points */
-  CCTK_INT i, j, k;
-  // CCTK_INT index = INITVALUE;
-  
   /* Declare finite differencing variables */
   
-  if (verbose > 1)
-  {
-    CCTK_VInfo(CCTK_THORNSTRING,"Entering eulerauto_cons_calc_primitives_Body");
-  }
-  
-  if (cctk_iteration % eulerauto_cons_calc_primitives_calc_every != eulerauto_cons_calc_primitives_calc_offset)
-  {
-    return;
-  }
-  
-  const char *groups[] = {"EulerAuto::Den_group","EulerAuto::En_group","EulerAuto::p_group","EulerAuto::rho_group","EulerAuto::S_group","EulerAuto::v_group"};
-  GenericFD_AssertGroupStorage(cctkGH, "eulerauto_cons_calc_primitives", 6, groups);
-  
-  
   /* Include user-supplied include files */
   
   /* Initialise finite differencing variables */
   ptrdiff_t const di = 1;
   ptrdiff_t const dj = CCTK_GFINDEX3D(cctkGH,0,1,0) - CCTK_GFINDEX3D(cctkGH,0,0,0);
   ptrdiff_t const dk = CCTK_GFINDEX3D(cctkGH,0,0,1) - CCTK_GFINDEX3D(cctkGH,0,0,0);
+  ptrdiff_t const cdi = sizeof(CCTK_REAL) * di;
+  ptrdiff_t const cdj = sizeof(CCTK_REAL) * dj;
+  ptrdiff_t const cdk = sizeof(CCTK_REAL) * dk;
   CCTK_REAL const dx = ToReal(CCTK_DELTA_SPACE(0));
   CCTK_REAL const dy = ToReal(CCTK_DELTA_SPACE(1));
   CCTK_REAL const dz = ToReal(CCTK_DELTA_SPACE(2));
   CCTK_REAL const dt = ToReal(CCTK_DELTA_TIME);
+  CCTK_REAL const t = ToReal(cctk_time);
   CCTK_REAL const dxi = INV(dx);
   CCTK_REAL const dyi = INV(dy);
   CCTK_REAL const dzi = INV(dz);
@@ -74,49 +62,55 @@ static void eulerauto_cons_calc_primitives_Body(cGH const * restrict const cctkG
   CCTK_REAL const p1ody = INV(dy);
   CCTK_REAL const p1odz = INV(dz);
   
+  /* Assign local copies of arrays functions */
+  
+  
+  
+  /* Calculate temporaries and arrays functions */
+  
+  /* Copy local copies back to grid functions */
+  
   /* Loop over the grid points */
-  for (k = min[2]; k < max[2]; k++)
+  #pragma omp parallel
+  CCTK_LOOP3 (eulerauto_cons_calc_primitives,
+    i,j,k, imin[0],imin[1],imin[2], imax[0],imax[1],imax[2],
+    cctk_lsh[0],cctk_lsh[1],cctk_lsh[2])
   {
-    for (j = min[1]; j < max[1]; j++)
-    {
-      for (i = min[0]; i < max[0]; i++)
-      {
-         int  const  index  =  CCTK_GFINDEX3D(cctkGH,i,j,k) ;
-        
-        /* Assign local copies of grid functions */
-        
-        CCTK_REAL DenL = Den[index];
-        CCTK_REAL EnL = En[index];
-        CCTK_REAL S1L = S1[index];
-        CCTK_REAL S2L = S2[index];
-        CCTK_REAL S3L = S3[index];
-        
-        
-        /* Include user supplied include files */
-        
-        /* Precompute derivatives */
-        
-        /* Calculate temporaries and grid functions */
-        CCTK_REAL rhoL = DenL;
-        
-        CCTK_REAL v1L = S1L*INV(DenL);
-        
-        CCTK_REAL v2L = S2L*INV(DenL);
-        
-        CCTK_REAL v3L = S3L*INV(DenL);
-        
-        CCTK_REAL pL = (EnL - 0.5*INV(DenL)*(SQR(S1L) + SQR(S2L) + 
-          SQR(S3L)))*(-1 + ToReal(gamma));
-        
-        /* Copy local copies back to grid functions */
-        p[index] = pL;
-        rho[index] = rhoL;
-        v1[index] = v1L;
-        v2[index] = v2L;
-        v3[index] = v3L;
-      }
-    }
+    ptrdiff_t const index = di*i + dj*j + dk*k;
+    
+    /* Assign local copies of grid functions */
+    
+    CCTK_REAL DenL = Den[index];
+    CCTK_REAL EnL = En[index];
+    CCTK_REAL S1L = S1[index];
+    CCTK_REAL S2L = S2[index];
+    CCTK_REAL S3L = S3[index];
+    
+    
+    /* Include user supplied include files */
+    
+    /* Precompute derivatives */
+    
+    /* Calculate temporaries and grid functions */
+    CCTK_REAL rhoL = DenL;
+    
+    CCTK_REAL v1L = S1L*INV(DenL);
+    
+    CCTK_REAL v2L = S2L*INV(DenL);
+    
+    CCTK_REAL v3L = S3L*INV(DenL);
+    
+    CCTK_REAL pL = 0.5*INV(DenL)*(2*DenL*EnL - SQR(S1L) - 
+      SQR(S2L) - SQR(S3L))*(-1 + ToReal(gamma));
+    
+    /* Copy local copies back to grid functions */
+    p[index] = pL;
+    rho[index] = rhoL;
+    v1[index] = v1L;
+    v2[index] = v2L;
+    v3[index] = v3L;
   }
+  CCTK_ENDLOOP3 (eulerauto_cons_calc_primitives);
 }
 
 extern "C" void eulerauto_cons_calc_primitives(CCTK_ARGUMENTS)
@@ -124,5 +118,25 @@ extern "C" void eulerauto_cons_calc_primitives(CCTK_ARGUMENTS)
   DECLARE_CCTK_ARGUMENTS;
   DECLARE_CCTK_PARAMETERS;
   
+  
+  if (verbose > 1)
+  {
+    CCTK_VInfo(CCTK_THORNSTRING,"Entering eulerauto_cons_calc_primitives_Body");
+  }
+  
+  if (cctk_iteration % eulerauto_cons_calc_primitives_calc_every != eulerauto_cons_calc_primitives_calc_offset)
+  {
+    return;
+  }
+  
+  const char *groups[] = {"EulerAuto::Den_group","EulerAuto::En_group","EulerAuto::p_group","EulerAuto::rho_group","EulerAuto::S_group","EulerAuto::v_group"};
+  GenericFD_AssertGroupStorage(cctkGH, "eulerauto_cons_calc_primitives", 6, groups);
+  
+  
   GenericFD_LoopOverEverything(cctkGH, &eulerauto_cons_calc_primitives_Body);
+  
+  if (verbose > 1)
+  {
+    CCTK_VInfo(CCTK_THORNSTRING,"Leaving eulerauto_cons_calc_primitives_Body");
+  }
 }
diff --git a/Examples/EulerAuto/src/eulerauto_cons_calc_reconstruct_1.cc b/Examples/EulerAuto/src/eulerauto_cons_calc_reconstruct_1.cc
index 0506ef0..be96c9f 100644
--- a/Examples/EulerAuto/src/eulerauto_cons_calc_reconstruct_1.cc
+++ b/Examples/EulerAuto/src/eulerauto_cons_calc_reconstruct_1.cc
@@ -12,6 +12,8 @@
 #include "cctk_Parameters.h"
 #include "GenericFD.h"
 #include "Differencing.h"
+#include "cctk_Loop.h"
+#include "loopcontrol.h"
 
 /* Define macros used in calculations */
 #define INITVALUE (42)
@@ -44,43 +46,28 @@ extern "C" void eulerauto_cons_calc_reconstruct_1_SelectBCs(CCTK_ARGUMENTS)
   return;
 }
 
-static void eulerauto_cons_calc_reconstruct_1_Body(cGH const * restrict const cctkGH, int const dir, int const face, CCTK_REAL const normal[3], CCTK_REAL const tangentA[3], CCTK_REAL const tangentB[3], int const min[3], int const max[3], int const n_subblock_gfs, CCTK_REAL * restrict const subblock_gfs[])
+static void eulerauto_cons_calc_reconstruct_1_Body(cGH const * restrict const cctkGH, int const dir, int const face, CCTK_REAL const normal[3], CCTK_REAL const tangentA[3], CCTK_REAL const tangentB[3], int const imin[3], int const imax[3], int const n_subblock_gfs, CCTK_REAL * restrict const subblock_gfs[])
 {
   DECLARE_CCTK_ARGUMENTS;
   DECLARE_CCTK_PARAMETERS;
   
   
-  /* Declare the variables used for looping over grid points */
-  CCTK_INT i, j, k;
-  // CCTK_INT index = INITVALUE;
-  
   /* Declare finite differencing variables */
   
-  if (verbose > 1)
-  {
-    CCTK_VInfo(CCTK_THORNSTRING,"Entering eulerauto_cons_calc_reconstruct_1_Body");
-  }
-  
-  if (cctk_iteration % eulerauto_cons_calc_reconstruct_1_calc_every != eulerauto_cons_calc_reconstruct_1_calc_offset)
-  {
-    return;
-  }
-  
-  const char *groups[] = {"EulerAuto::p_group","EulerAuto::p_lr_group","EulerAuto::rho_group","EulerAuto::rho_lr_group","EulerAuto::v1_lr_group","EulerAuto::v2_lr_group","EulerAuto::v3_lr_group","EulerAuto::v_group"};
-  GenericFD_AssertGroupStorage(cctkGH, "eulerauto_cons_calc_reconstruct_1", 8, groups);
-  
-  GenericFD_EnsureStencilFits(cctkGH, "eulerauto_cons_calc_reconstruct_1", 1, 1, 1);
-  
   /* Include user-supplied include files */
   
   /* Initialise finite differencing variables */
   ptrdiff_t const di = 1;
   ptrdiff_t const dj = CCTK_GFINDEX3D(cctkGH,0,1,0) - CCTK_GFINDEX3D(cctkGH,0,0,0);
   ptrdiff_t const dk = CCTK_GFINDEX3D(cctkGH,0,0,1) - CCTK_GFINDEX3D(cctkGH,0,0,0);
+  ptrdiff_t const cdi = sizeof(CCTK_REAL) * di;
+  ptrdiff_t const cdj = sizeof(CCTK_REAL) * dj;
+  ptrdiff_t const cdk = sizeof(CCTK_REAL) * dk;
   CCTK_REAL const dx = ToReal(CCTK_DELTA_SPACE(0));
   CCTK_REAL const dy = ToReal(CCTK_DELTA_SPACE(1));
   CCTK_REAL const dz = ToReal(CCTK_DELTA_SPACE(2));
   CCTK_REAL const dt = ToReal(CCTK_DELTA_TIME);
+  CCTK_REAL const t = ToReal(cctk_time);
   CCTK_REAL const dxi = INV(dx);
   CCTK_REAL const dyi = INV(dy);
   CCTK_REAL const dzi = INV(dz);
@@ -99,103 +86,109 @@ static void eulerauto_cons_calc_reconstruct_1_Body(cGH const * restrict const cc
   CCTK_REAL const p1ody = INV(dy);
   CCTK_REAL const p1odz = INV(dz);
   
+  /* Assign local copies of arrays functions */
+  
+  
+  
+  /* Calculate temporaries and arrays functions */
+  
+  /* Copy local copies back to grid functions */
+  
   /* Loop over the grid points */
-  for (k = min[2]; k < max[2]; k++)
+  #pragma omp parallel
+  CCTK_LOOP3 (eulerauto_cons_calc_reconstruct_1,
+    i,j,k, imin[0],imin[1],imin[2], imax[0],imax[1],imax[2],
+    cctk_lsh[0],cctk_lsh[1],cctk_lsh[2])
   {
-    for (j = min[1]; j < max[1]; j++)
-    {
-      for (i = min[0]; i < max[0]; i++)
-      {
-         int  const  index  =  CCTK_GFINDEX3D(cctkGH,i,j,k) ;
-        
-        /* Assign local copies of grid functions */
-        
-        CCTK_REAL pL = p[index];
-        CCTK_REAL rhoL = rho[index];
-        CCTK_REAL v1L = v1[index];
-        CCTK_REAL v2L = v2[index];
-        CCTK_REAL v3L = v3[index];
-        
-        
-        /* Include user supplied include files */
-        
-        /* Precompute derivatives */
-        CCTK_REAL const DiffPlus1p = DiffPlus1(&p[index]);
-        CCTK_REAL const DiffMinus1p = DiffMinus1(&p[index]);
-        CCTK_REAL const DiffPlus1rho = DiffPlus1(&rho[index]);
-        CCTK_REAL const DiffMinus1rho = DiffMinus1(&rho[index]);
-        CCTK_REAL const DiffPlus1v1 = DiffPlus1(&v1[index]);
-        CCTK_REAL const DiffMinus1v1 = DiffMinus1(&v1[index]);
-        CCTK_REAL const DiffPlus1v2 = DiffPlus1(&v2[index]);
-        CCTK_REAL const DiffMinus1v2 = DiffMinus1(&v2[index]);
-        CCTK_REAL const DiffPlus1v3 = DiffPlus1(&v3[index]);
-        CCTK_REAL const DiffMinus1v3 = DiffMinus1(&v3[index]);
-        
-        /* Calculate temporaries and grid functions */
-        CCTK_REAL slopeL = DiffMinus1rho;
-        
-        CCTK_REAL slopeR = DiffPlus1rho;
-        
-        CCTK_REAL slope = VanLeer(slopeL,slopeR);
-        
-        CCTK_REAL rhoLeftL = rhoL - 0.5*slope;
-        
-        CCTK_REAL rhoRightL = rhoL + 0.5*slope;
-        
-        slopeL = DiffMinus1v1;
-        
-        slopeR = DiffPlus1v1;
-        
-        slope = VanLeer(slopeL,slopeR);
-        
-        CCTK_REAL v1LeftL = -0.5*slope + v1L;
-        
-        CCTK_REAL v1RightL = 0.5*slope + v1L;
-        
-        slopeL = DiffMinus1v2;
-        
-        slopeR = DiffPlus1v2;
-        
-        slope = VanLeer(slopeL,slopeR);
-        
-        CCTK_REAL v2LeftL = -0.5*slope + v2L;
-        
-        CCTK_REAL v2RightL = 0.5*slope + v2L;
-        
-        slopeL = DiffMinus1v3;
-        
-        slopeR = DiffPlus1v3;
-        
-        slope = VanLeer(slopeL,slopeR);
-        
-        CCTK_REAL v3LeftL = -0.5*slope + v3L;
-        
-        CCTK_REAL v3RightL = 0.5*slope + v3L;
-        
-        slopeL = DiffMinus1p;
-        
-        slopeR = DiffPlus1p;
-        
-        slope = VanLeer(slopeL,slopeR);
-        
-        CCTK_REAL pLeftL = pL - 0.5*slope;
-        
-        CCTK_REAL pRightL = pL + 0.5*slope;
-        
-        /* Copy local copies back to grid functions */
-        pLeft[index] = pLeftL;
-        pRight[index] = pRightL;
-        rhoLeft[index] = rhoLeftL;
-        rhoRight[index] = rhoRightL;
-        v1Left[index] = v1LeftL;
-        v1Right[index] = v1RightL;
-        v2Left[index] = v2LeftL;
-        v2Right[index] = v2RightL;
-        v3Left[index] = v3LeftL;
-        v3Right[index] = v3RightL;
-      }
-    }
+    ptrdiff_t const index = di*i + dj*j + dk*k;
+    
+    /* Assign local copies of grid functions */
+    
+    CCTK_REAL pL = p[index];
+    CCTK_REAL rhoL = rho[index];
+    CCTK_REAL v1L = v1[index];
+    CCTK_REAL v2L = v2[index];
+    CCTK_REAL v3L = v3[index];
+    
+    
+    /* Include user supplied include files */
+    
+    /* Precompute derivatives */
+    CCTK_REAL const DiffPlus1p = DiffPlus1(&p[index]);
+    CCTK_REAL const DiffMinus1p = DiffMinus1(&p[index]);
+    CCTK_REAL const DiffPlus1rho = DiffPlus1(&rho[index]);
+    CCTK_REAL const DiffMinus1rho = DiffMinus1(&rho[index]);
+    CCTK_REAL const DiffPlus1v1 = DiffPlus1(&v1[index]);
+    CCTK_REAL const DiffMinus1v1 = DiffMinus1(&v1[index]);
+    CCTK_REAL const DiffPlus1v2 = DiffPlus1(&v2[index]);
+    CCTK_REAL const DiffMinus1v2 = DiffMinus1(&v2[index]);
+    CCTK_REAL const DiffPlus1v3 = DiffPlus1(&v3[index]);
+    CCTK_REAL const DiffMinus1v3 = DiffMinus1(&v3[index]);
+    
+    /* Calculate temporaries and grid functions */
+    CCTK_REAL slopeL = DiffMinus1rho;
+    
+    CCTK_REAL slopeR = DiffPlus1rho;
+    
+    CCTK_REAL slope = VanLeer(slopeL,slopeR);
+    
+    CCTK_REAL rhoLeftL = rhoL - 0.5*slope;
+    
+    CCTK_REAL rhoRightL = rhoL + 0.5*slope;
+    
+    slopeL = DiffMinus1v1;
+    
+    slopeR = DiffPlus1v1;
+    
+    slope = VanLeer(slopeL,slopeR);
+    
+    CCTK_REAL v1LeftL = v1L - 0.5*slope;
+    
+    CCTK_REAL v1RightL = v1L + 0.5*slope;
+    
+    slopeL = DiffMinus1v2;
+    
+    slopeR = DiffPlus1v2;
+    
+    slope = VanLeer(slopeL,slopeR);
+    
+    CCTK_REAL v2LeftL = v2L - 0.5*slope;
+    
+    CCTK_REAL v2RightL = v2L + 0.5*slope;
+    
+    slopeL = DiffMinus1v3;
+    
+    slopeR = DiffPlus1v3;
+    
+    slope = VanLeer(slopeL,slopeR);
+    
+    CCTK_REAL v3LeftL = v3L - 0.5*slope;
+    
+    CCTK_REAL v3RightL = v3L + 0.5*slope;
+    
+    slopeL = DiffMinus1p;
+    
+    slopeR = DiffPlus1p;
+    
+    slope = VanLeer(slopeL,slopeR);
+    
+    CCTK_REAL pLeftL = pL - 0.5*slope;
+    
+    CCTK_REAL pRightL = pL + 0.5*slope;
+    
+    /* Copy local copies back to grid functions */
+    pLeft[index] = pLeftL;
+    pRight[index] = pRightL;
+    rhoLeft[index] = rhoLeftL;
+    rhoRight[index] = rhoRightL;
+    v1Left[index] = v1LeftL;
+    v1Right[index] = v1RightL;
+    v2Left[index] = v2LeftL;
+    v2Right[index] = v2RightL;
+    v3Left[index] = v3LeftL;
+    v3Right[index] = v3RightL;
   }
+  CCTK_ENDLOOP3 (eulerauto_cons_calc_reconstruct_1);
 }
 
 extern "C" void eulerauto_cons_calc_reconstruct_1(CCTK_ARGUMENTS)
@@ -203,5 +196,26 @@ extern "C" void eulerauto_cons_calc_reconstruct_1(CCTK_ARGUMENTS)
   DECLARE_CCTK_ARGUMENTS;
   DECLARE_CCTK_PARAMETERS;
   
+  
+  if (verbose > 1)
+  {
+    CCTK_VInfo(CCTK_THORNSTRING,"Entering eulerauto_cons_calc_reconstruct_1_Body");
+  }
+  
+  if (cctk_iteration % eulerauto_cons_calc_reconstruct_1_calc_every != eulerauto_cons_calc_reconstruct_1_calc_offset)
+  {
+    return;
+  }
+  
+  const char *groups[] = {"EulerAuto::p_group","EulerAuto::p_lr_group","EulerAuto::rho_group","EulerAuto::rho_lr_group","EulerAuto::v1_lr_group","EulerAuto::v2_lr_group","EulerAuto::v3_lr_group","EulerAuto::v_group"};
+  GenericFD_AssertGroupStorage(cctkGH, "eulerauto_cons_calc_reconstruct_1", 8, groups);
+  
+  GenericFD_EnsureStencilFits(cctkGH, "eulerauto_cons_calc_reconstruct_1", 1, 1, 1);
+  
   GenericFD_LoopOverInterior(cctkGH, &eulerauto_cons_calc_reconstruct_1_Body);
+  
+  if (verbose > 1)
+  {
+    CCTK_VInfo(CCTK_THORNSTRING,"Leaving eulerauto_cons_calc_reconstruct_1_Body");
+  }
 }
diff --git a/Examples/EulerAuto/src/eulerauto_cons_calc_reconstruct_2.cc b/Examples/EulerAuto/src/eulerauto_cons_calc_reconstruct_2.cc
index f71f106..34b7f55 100644
--- a/Examples/EulerAuto/src/eulerauto_cons_calc_reconstruct_2.cc
+++ b/Examples/EulerAuto/src/eulerauto_cons_calc_reconstruct_2.cc
@@ -12,6 +12,8 @@
 #include "cctk_Parameters.h"
 #include "GenericFD.h"
 #include "Differencing.h"
+#include "cctk_Loop.h"
+#include "loopcontrol.h"
 
 /* Define macros used in calculations */
 #define INITVALUE (42)
@@ -44,43 +46,28 @@ extern "C" void eulerauto_cons_calc_reconstruct_2_SelectBCs(CCTK_ARGUMENTS)
   return;
 }
 
-static void eulerauto_cons_calc_reconstruct_2_Body(cGH const * restrict const cctkGH, int const dir, int const face, CCTK_REAL const normal[3], CCTK_REAL const tangentA[3], CCTK_REAL const tangentB[3], int const min[3], int const max[3], int const n_subblock_gfs, CCTK_REAL * restrict const subblock_gfs[])
+static void eulerauto_cons_calc_reconstruct_2_Body(cGH const * restrict const cctkGH, int const dir, int const face, CCTK_REAL const normal[3], CCTK_REAL const tangentA[3], CCTK_REAL const tangentB[3], int const imin[3], int const imax[3], int const n_subblock_gfs, CCTK_REAL * restrict const subblock_gfs[])
 {
   DECLARE_CCTK_ARGUMENTS;
   DECLARE_CCTK_PARAMETERS;
   
   
-  /* Declare the variables used for looping over grid points */
-  CCTK_INT i, j, k;
-  // CCTK_INT index = INITVALUE;
-  
   /* Declare finite differencing variables */
   
-  if (verbose > 1)
-  {
-    CCTK_VInfo(CCTK_THORNSTRING,"Entering eulerauto_cons_calc_reconstruct_2_Body");
-  }
-  
-  if (cctk_iteration % eulerauto_cons_calc_reconstruct_2_calc_every != eulerauto_cons_calc_reconstruct_2_calc_offset)
-  {
-    return;
-  }
-  
-  const char *groups[] = {"EulerAuto::p_group","EulerAuto::p_lr_group","EulerAuto::rho_group","EulerAuto::rho_lr_group","EulerAuto::v1_lr_group","EulerAuto::v2_lr_group","EulerAuto::v3_lr_group","EulerAuto::v_group"};
-  GenericFD_AssertGroupStorage(cctkGH, "eulerauto_cons_calc_reconstruct_2", 8, groups);
-  
-  GenericFD_EnsureStencilFits(cctkGH, "eulerauto_cons_calc_reconstruct_2", 1, 1, 1);
-  
   /* Include user-supplied include files */
   
   /* Initialise finite differencing variables */
   ptrdiff_t const di = 1;
   ptrdiff_t const dj = CCTK_GFINDEX3D(cctkGH,0,1,0) - CCTK_GFINDEX3D(cctkGH,0,0,0);
   ptrdiff_t const dk = CCTK_GFINDEX3D(cctkGH,0,0,1) - CCTK_GFINDEX3D(cctkGH,0,0,0);
+  ptrdiff_t const cdi = sizeof(CCTK_REAL) * di;
+  ptrdiff_t const cdj = sizeof(CCTK_REAL) * dj;
+  ptrdiff_t const cdk = sizeof(CCTK_REAL) * dk;
   CCTK_REAL const dx = ToReal(CCTK_DELTA_SPACE(0));
   CCTK_REAL const dy = ToReal(CCTK_DELTA_SPACE(1));
   CCTK_REAL const dz = ToReal(CCTK_DELTA_SPACE(2));
   CCTK_REAL const dt = ToReal(CCTK_DELTA_TIME);
+  CCTK_REAL const t = ToReal(cctk_time);
   CCTK_REAL const dxi = INV(dx);
   CCTK_REAL const dyi = INV(dy);
   CCTK_REAL const dzi = INV(dz);
@@ -99,103 +86,109 @@ static void eulerauto_cons_calc_reconstruct_2_Body(cGH const * restrict const cc
   CCTK_REAL const p1ody = INV(dy);
   CCTK_REAL const p1odz = INV(dz);
   
+  /* Assign local copies of arrays functions */
+  
+  
+  
+  /* Calculate temporaries and arrays functions */
+  
+  /* Copy local copies back to grid functions */
+  
   /* Loop over the grid points */
-  for (k = min[2]; k < max[2]; k++)
+  #pragma omp parallel
+  CCTK_LOOP3 (eulerauto_cons_calc_reconstruct_2,
+    i,j,k, imin[0],imin[1],imin[2], imax[0],imax[1],imax[2],
+    cctk_lsh[0],cctk_lsh[1],cctk_lsh[2])
   {
-    for (j = min[1]; j < max[1]; j++)
-    {
-      for (i = min[0]; i < max[0]; i++)
-      {
-         int  const  index  =  CCTK_GFINDEX3D(cctkGH,i,j,k) ;
-        
-        /* Assign local copies of grid functions */
-        
-        CCTK_REAL pL = p[index];
-        CCTK_REAL rhoL = rho[index];
-        CCTK_REAL v1L = v1[index];
-        CCTK_REAL v2L = v2[index];
-        CCTK_REAL v3L = v3[index];
-        
-        
-        /* Include user supplied include files */
-        
-        /* Precompute derivatives */
-        CCTK_REAL const DiffPlus2p = DiffPlus2(&p[index]);
-        CCTK_REAL const DiffMinus2p = DiffMinus2(&p[index]);
-        CCTK_REAL const DiffPlus2rho = DiffPlus2(&rho[index]);
-        CCTK_REAL const DiffMinus2rho = DiffMinus2(&rho[index]);
-        CCTK_REAL const DiffPlus2v1 = DiffPlus2(&v1[index]);
-        CCTK_REAL const DiffMinus2v1 = DiffMinus2(&v1[index]);
-        CCTK_REAL const DiffPlus2v2 = DiffPlus2(&v2[index]);
-        CCTK_REAL const DiffMinus2v2 = DiffMinus2(&v2[index]);
-        CCTK_REAL const DiffPlus2v3 = DiffPlus2(&v3[index]);
-        CCTK_REAL const DiffMinus2v3 = DiffMinus2(&v3[index]);
-        
-        /* Calculate temporaries and grid functions */
-        CCTK_REAL slopeL = DiffMinus2rho;
-        
-        CCTK_REAL slopeR = DiffPlus2rho;
-        
-        CCTK_REAL slope = VanLeer(slopeL,slopeR);
-        
-        CCTK_REAL rhoLeftL = rhoL - 0.5*slope;
-        
-        CCTK_REAL rhoRightL = rhoL + 0.5*slope;
-        
-        slopeL = DiffMinus2v1;
-        
-        slopeR = DiffPlus2v1;
-        
-        slope = VanLeer(slopeL,slopeR);
-        
-        CCTK_REAL v1LeftL = -0.5*slope + v1L;
-        
-        CCTK_REAL v1RightL = 0.5*slope + v1L;
-        
-        slopeL = DiffMinus2v2;
-        
-        slopeR = DiffPlus2v2;
-        
-        slope = VanLeer(slopeL,slopeR);
-        
-        CCTK_REAL v2LeftL = -0.5*slope + v2L;
-        
-        CCTK_REAL v2RightL = 0.5*slope + v2L;
-        
-        slopeL = DiffMinus2v3;
-        
-        slopeR = DiffPlus2v3;
-        
-        slope = VanLeer(slopeL,slopeR);
-        
-        CCTK_REAL v3LeftL = -0.5*slope + v3L;
-        
-        CCTK_REAL v3RightL = 0.5*slope + v3L;
-        
-        slopeL = DiffMinus2p;
-        
-        slopeR = DiffPlus2p;
-        
-        slope = VanLeer(slopeL,slopeR);
-        
-        CCTK_REAL pLeftL = pL - 0.5*slope;
-        
-        CCTK_REAL pRightL = pL + 0.5*slope;
-        
-        /* Copy local copies back to grid functions */
-        pLeft[index] = pLeftL;
-        pRight[index] = pRightL;
-        rhoLeft[index] = rhoLeftL;
-        rhoRight[index] = rhoRightL;
-        v1Left[index] = v1LeftL;
-        v1Right[index] = v1RightL;
-        v2Left[index] = v2LeftL;
-        v2Right[index] = v2RightL;
-        v3Left[index] = v3LeftL;
-        v3Right[index] = v3RightL;
-      }
-    }
+    ptrdiff_t const index = di*i + dj*j + dk*k;
+    
+    /* Assign local copies of grid functions */
+    
+    CCTK_REAL pL = p[index];
+    CCTK_REAL rhoL = rho[index];
+    CCTK_REAL v1L = v1[index];
+    CCTK_REAL v2L = v2[index];
+    CCTK_REAL v3L = v3[index];
+    
+    
+    /* Include user supplied include files */
+    
+    /* Precompute derivatives */
+    CCTK_REAL const DiffPlus2p = DiffPlus2(&p[index]);
+    CCTK_REAL const DiffMinus2p = DiffMinus2(&p[index]);
+    CCTK_REAL const DiffPlus2rho = DiffPlus2(&rho[index]);
+    CCTK_REAL const DiffMinus2rho = DiffMinus2(&rho[index]);
+    CCTK_REAL const DiffPlus2v1 = DiffPlus2(&v1[index]);
+    CCTK_REAL const DiffMinus2v1 = DiffMinus2(&v1[index]);
+    CCTK_REAL const DiffPlus2v2 = DiffPlus2(&v2[index]);
+    CCTK_REAL const DiffMinus2v2 = DiffMinus2(&v2[index]);
+    CCTK_REAL const DiffPlus2v3 = DiffPlus2(&v3[index]);
+    CCTK_REAL const DiffMinus2v3 = DiffMinus2(&v3[index]);
+    
+    /* Calculate temporaries and grid functions */
+    CCTK_REAL slopeL = DiffMinus2rho;
+    
+    CCTK_REAL slopeR = DiffPlus2rho;
+    
+    CCTK_REAL slope = VanLeer(slopeL,slopeR);
+    
+    CCTK_REAL rhoLeftL = rhoL - 0.5*slope;
+    
+    CCTK_REAL rhoRightL = rhoL + 0.5*slope;
+    
+    slopeL = DiffMinus2v1;
+    
+    slopeR = DiffPlus2v1;
+    
+    slope = VanLeer(slopeL,slopeR);
+    
+    CCTK_REAL v1LeftL = v1L - 0.5*slope;
+    
+    CCTK_REAL v1RightL = v1L + 0.5*slope;
+    
+    slopeL = DiffMinus2v2;
+    
+    slopeR = DiffPlus2v2;
+    
+    slope = VanLeer(slopeL,slopeR);
+    
+    CCTK_REAL v2LeftL = v2L - 0.5*slope;
+    
+    CCTK_REAL v2RightL = v2L + 0.5*slope;
+    
+    slopeL = DiffMinus2v3;
+    
+    slopeR = DiffPlus2v3;
+    
+    slope = VanLeer(slopeL,slopeR);
+    
+    CCTK_REAL v3LeftL = v3L - 0.5*slope;
+    
+    CCTK_REAL v3RightL = v3L + 0.5*slope;
+    
+    slopeL = DiffMinus2p;
+    
+    slopeR = DiffPlus2p;
+    
+    slope = VanLeer(slopeL,slopeR);
+    
+    CCTK_REAL pLeftL = pL - 0.5*slope;
+    
+    CCTK_REAL pRightL = pL + 0.5*slope;
+    
+    /* Copy local copies back to grid functions */
+    pLeft[index] = pLeftL;
+    pRight[index] = pRightL;
+    rhoLeft[index] = rhoLeftL;
+    rhoRight[index] = rhoRightL;
+    v1Left[index] = v1LeftL;
+    v1Right[index] = v1RightL;
+    v2Left[index] = v2LeftL;
+    v2Right[index] = v2RightL;
+    v3Left[index] = v3LeftL;
+    v3Right[index] = v3RightL;
   }
+  CCTK_ENDLOOP3 (eulerauto_cons_calc_reconstruct_2);
 }
 
 extern "C" void eulerauto_cons_calc_reconstruct_2(CCTK_ARGUMENTS)
@@ -203,5 +196,26 @@ extern "C" void eulerauto_cons_calc_reconstruct_2(CCTK_ARGUMENTS)
   DECLARE_CCTK_ARGUMENTS;
   DECLARE_CCTK_PARAMETERS;
   
+  
+  if (verbose > 1)
+  {
+    CCTK_VInfo(CCTK_THORNSTRING,"Entering eulerauto_cons_calc_reconstruct_2_Body");
+  }
+  
+  if (cctk_iteration % eulerauto_cons_calc_reconstruct_2_calc_every != eulerauto_cons_calc_reconstruct_2_calc_offset)
+  {
+    return;
+  }
+  
+  const char *groups[] = {"EulerAuto::p_group","EulerAuto::p_lr_group","EulerAuto::rho_group","EulerAuto::rho_lr_group","EulerAuto::v1_lr_group","EulerAuto::v2_lr_group","EulerAuto::v3_lr_group","EulerAuto::v_group"};
+  GenericFD_AssertGroupStorage(cctkGH, "eulerauto_cons_calc_reconstruct_2", 8, groups);
+  
+  GenericFD_EnsureStencilFits(cctkGH, "eulerauto_cons_calc_reconstruct_2", 1, 1, 1);
+  
   GenericFD_LoopOverInterior(cctkGH, &eulerauto_cons_calc_reconstruct_2_Body);
+  
+  if (verbose > 1)
+  {
+    CCTK_VInfo(CCTK_THORNSTRING,"Leaving eulerauto_cons_calc_reconstruct_2_Body");
+  }
 }
diff --git a/Examples/EulerAuto/src/eulerauto_cons_calc_reconstruct_3.cc b/Examples/EulerAuto/src/eulerauto_cons_calc_reconstruct_3.cc
index 5204483..5481fba 100644
--- a/Examples/EulerAuto/src/eulerauto_cons_calc_reconstruct_3.cc
+++ b/Examples/EulerAuto/src/eulerauto_cons_calc_reconstruct_3.cc
@@ -12,6 +12,8 @@
 #include "cctk_Parameters.h"
 #include "GenericFD.h"
 #include "Differencing.h"
+#include "cctk_Loop.h"
+#include "loopcontrol.h"
 
 /* Define macros used in calculations */
 #define INITVALUE (42)
@@ -44,43 +46,28 @@ extern "C" void eulerauto_cons_calc_reconstruct_3_SelectBCs(CCTK_ARGUMENTS)
   return;
 }
 
-static void eulerauto_cons_calc_reconstruct_3_Body(cGH const * restrict const cctkGH, int const dir, int const face, CCTK_REAL const normal[3], CCTK_REAL const tangentA[3], CCTK_REAL const tangentB[3], int const min[3], int const max[3], int const n_subblock_gfs, CCTK_REAL * restrict const subblock_gfs[])
+static void eulerauto_cons_calc_reconstruct_3_Body(cGH const * restrict const cctkGH, int const dir, int const face, CCTK_REAL const normal[3], CCTK_REAL const tangentA[3], CCTK_REAL const tangentB[3], int const imin[3], int const imax[3], int const n_subblock_gfs, CCTK_REAL * restrict const subblock_gfs[])
 {
   DECLARE_CCTK_ARGUMENTS;
   DECLARE_CCTK_PARAMETERS;
   
   
-  /* Declare the variables used for looping over grid points */
-  CCTK_INT i, j, k;
-  // CCTK_INT index = INITVALUE;
-  
   /* Declare finite differencing variables */
   
-  if (verbose > 1)
-  {
-    CCTK_VInfo(CCTK_THORNSTRING,"Entering eulerauto_cons_calc_reconstruct_3_Body");
-  }
-  
-  if (cctk_iteration % eulerauto_cons_calc_reconstruct_3_calc_every != eulerauto_cons_calc_reconstruct_3_calc_offset)
-  {
-    return;
-  }
-  
-  const char *groups[] = {"EulerAuto::p_group","EulerAuto::p_lr_group","EulerAuto::rho_group","EulerAuto::rho_lr_group","EulerAuto::v1_lr_group","EulerAuto::v2_lr_group","EulerAuto::v3_lr_group","EulerAuto::v_group"};
-  GenericFD_AssertGroupStorage(cctkGH, "eulerauto_cons_calc_reconstruct_3", 8, groups);
-  
-  GenericFD_EnsureStencilFits(cctkGH, "eulerauto_cons_calc_reconstruct_3", 1, 1, 1);
-  
   /* Include user-supplied include files */
   
   /* Initialise finite differencing variables */
   ptrdiff_t const di = 1;
   ptrdiff_t const dj = CCTK_GFINDEX3D(cctkGH,0,1,0) - CCTK_GFINDEX3D(cctkGH,0,0,0);
   ptrdiff_t const dk = CCTK_GFINDEX3D(cctkGH,0,0,1) - CCTK_GFINDEX3D(cctkGH,0,0,0);
+  ptrdiff_t const cdi = sizeof(CCTK_REAL) * di;
+  ptrdiff_t const cdj = sizeof(CCTK_REAL) * dj;
+  ptrdiff_t const cdk = sizeof(CCTK_REAL) * dk;
   CCTK_REAL const dx = ToReal(CCTK_DELTA_SPACE(0));
   CCTK_REAL const dy = ToReal(CCTK_DELTA_SPACE(1));
   CCTK_REAL const dz = ToReal(CCTK_DELTA_SPACE(2));
   CCTK_REAL const dt = ToReal(CCTK_DELTA_TIME);
+  CCTK_REAL const t = ToReal(cctk_time);
   CCTK_REAL const dxi = INV(dx);
   CCTK_REAL const dyi = INV(dy);
   CCTK_REAL const dzi = INV(dz);
@@ -99,103 +86,109 @@ static void eulerauto_cons_calc_reconstruct_3_Body(cGH const * restrict const cc
   CCTK_REAL const p1ody = INV(dy);
   CCTK_REAL const p1odz = INV(dz);
   
+  /* Assign local copies of arrays functions */
+  
+  
+  
+  /* Calculate temporaries and arrays functions */
+  
+  /* Copy local copies back to grid functions */
+  
   /* Loop over the grid points */
-  for (k = min[2]; k < max[2]; k++)
+  #pragma omp parallel
+  CCTK_LOOP3 (eulerauto_cons_calc_reconstruct_3,
+    i,j,k, imin[0],imin[1],imin[2], imax[0],imax[1],imax[2],
+    cctk_lsh[0],cctk_lsh[1],cctk_lsh[2])
   {
-    for (j = min[1]; j < max[1]; j++)
-    {
-      for (i = min[0]; i < max[0]; i++)
-      {
-         int  const  index  =  CCTK_GFINDEX3D(cctkGH,i,j,k) ;
-        
-        /* Assign local copies of grid functions */
-        
-        CCTK_REAL pL = p[index];
-        CCTK_REAL rhoL = rho[index];
-        CCTK_REAL v1L = v1[index];
-        CCTK_REAL v2L = v2[index];
-        CCTK_REAL v3L = v3[index];
-        
-        
-        /* Include user supplied include files */
-        
-        /* Precompute derivatives */
-        CCTK_REAL const DiffPlus3p = DiffPlus3(&p[index]);
-        CCTK_REAL const DiffMinus3p = DiffMinus3(&p[index]);
-        CCTK_REAL const DiffPlus3rho = DiffPlus3(&rho[index]);
-        CCTK_REAL const DiffMinus3rho = DiffMinus3(&rho[index]);
-        CCTK_REAL const DiffPlus3v1 = DiffPlus3(&v1[index]);
-        CCTK_REAL const DiffMinus3v1 = DiffMinus3(&v1[index]);
-        CCTK_REAL const DiffPlus3v2 = DiffPlus3(&v2[index]);
-        CCTK_REAL const DiffMinus3v2 = DiffMinus3(&v2[index]);
-        CCTK_REAL const DiffPlus3v3 = DiffPlus3(&v3[index]);
-        CCTK_REAL const DiffMinus3v3 = DiffMinus3(&v3[index]);
-        
-        /* Calculate temporaries and grid functions */
-        CCTK_REAL slopeL = DiffMinus3rho;
-        
-        CCTK_REAL slopeR = DiffPlus3rho;
-        
-        CCTK_REAL slope = VanLeer(slopeL,slopeR);
-        
-        CCTK_REAL rhoLeftL = rhoL - 0.5*slope;
-        
-        CCTK_REAL rhoRightL = rhoL + 0.5*slope;
-        
-        slopeL = DiffMinus3v1;
-        
-        slopeR = DiffPlus3v1;
-        
-        slope = VanLeer(slopeL,slopeR);
-        
-        CCTK_REAL v1LeftL = -0.5*slope + v1L;
-        
-        CCTK_REAL v1RightL = 0.5*slope + v1L;
-        
-        slopeL = DiffMinus3v2;
-        
-        slopeR = DiffPlus3v2;
-        
-        slope = VanLeer(slopeL,slopeR);
-        
-        CCTK_REAL v2LeftL = -0.5*slope + v2L;
-        
-        CCTK_REAL v2RightL = 0.5*slope + v2L;
-        
-        slopeL = DiffMinus3v3;
-        
-        slopeR = DiffPlus3v3;
-        
-        slope = VanLeer(slopeL,slopeR);
-        
-        CCTK_REAL v3LeftL = -0.5*slope + v3L;
-        
-        CCTK_REAL v3RightL = 0.5*slope + v3L;
-        
-        slopeL = DiffMinus3p;
-        
-        slopeR = DiffPlus3p;
-        
-        slope = VanLeer(slopeL,slopeR);
-        
-        CCTK_REAL pLeftL = pL - 0.5*slope;
-        
-        CCTK_REAL pRightL = pL + 0.5*slope;
-        
-        /* Copy local copies back to grid functions */
-        pLeft[index] = pLeftL;
-        pRight[index] = pRightL;
-        rhoLeft[index] = rhoLeftL;
-        rhoRight[index] = rhoRightL;
-        v1Left[index] = v1LeftL;
-        v1Right[index] = v1RightL;
-        v2Left[index] = v2LeftL;
-        v2Right[index] = v2RightL;
-        v3Left[index] = v3LeftL;
-        v3Right[index] = v3RightL;
-      }
-    }
+    ptrdiff_t const index = di*i + dj*j + dk*k;
+    
+    /* Assign local copies of grid functions */
+    
+    CCTK_REAL pL = p[index];
+    CCTK_REAL rhoL = rho[index];
+    CCTK_REAL v1L = v1[index];
+    CCTK_REAL v2L = v2[index];
+    CCTK_REAL v3L = v3[index];
+    
+    
+    /* Include user supplied include files */
+    
+    /* Precompute derivatives */
+    CCTK_REAL const DiffPlus3p = DiffPlus3(&p[index]);
+    CCTK_REAL const DiffMinus3p = DiffMinus3(&p[index]);
+    CCTK_REAL const DiffPlus3rho = DiffPlus3(&rho[index]);
+    CCTK_REAL const DiffMinus3rho = DiffMinus3(&rho[index]);
+    CCTK_REAL const DiffPlus3v1 = DiffPlus3(&v1[index]);
+    CCTK_REAL const DiffMinus3v1 = DiffMinus3(&v1[index]);
+    CCTK_REAL const DiffPlus3v2 = DiffPlus3(&v2[index]);
+    CCTK_REAL const DiffMinus3v2 = DiffMinus3(&v2[index]);
+    CCTK_REAL const DiffPlus3v3 = DiffPlus3(&v3[index]);
+    CCTK_REAL const DiffMinus3v3 = DiffMinus3(&v3[index]);
+    
+    /* Calculate temporaries and grid functions */
+    CCTK_REAL slopeL = DiffMinus3rho;
+    
+    CCTK_REAL slopeR = DiffPlus3rho;
+    
+    CCTK_REAL slope = VanLeer(slopeL,slopeR);
+    
+    CCTK_REAL rhoLeftL = rhoL - 0.5*slope;
+    
+    CCTK_REAL rhoRightL = rhoL + 0.5*slope;
+    
+    slopeL = DiffMinus3v1;
+    
+    slopeR = DiffPlus3v1;
+    
+    slope = VanLeer(slopeL,slopeR);
+    
+    CCTK_REAL v1LeftL = v1L - 0.5*slope;
+    
+    CCTK_REAL v1RightL = v1L + 0.5*slope;
+    
+    slopeL = DiffMinus3v2;
+    
+    slopeR = DiffPlus3v2;
+    
+    slope = VanLeer(slopeL,slopeR);
+    
+    CCTK_REAL v2LeftL = v2L - 0.5*slope;
+    
+    CCTK_REAL v2RightL = v2L + 0.5*slope;
+    
+    slopeL = DiffMinus3v3;
+    
+    slopeR = DiffPlus3v3;
+    
+    slope = VanLeer(slopeL,slopeR);
+    
+    CCTK_REAL v3LeftL = v3L - 0.5*slope;
+    
+    CCTK_REAL v3RightL = v3L + 0.5*slope;
+    
+    slopeL = DiffMinus3p;
+    
+    slopeR = DiffPlus3p;
+    
+    slope = VanLeer(slopeL,slopeR);
+    
+    CCTK_REAL pLeftL = pL - 0.5*slope;
+    
+    CCTK_REAL pRightL = pL + 0.5*slope;
+    
+    /* Copy local copies back to grid functions */
+    pLeft[index] = pLeftL;
+    pRight[index] = pRightL;
+    rhoLeft[index] = rhoLeftL;
+    rhoRight[index] = rhoRightL;
+    v1Left[index] = v1LeftL;
+    v1Right[index] = v1RightL;
+    v2Left[index] = v2LeftL;
+    v2Right[index] = v2RightL;
+    v3Left[index] = v3LeftL;
+    v3Right[index] = v3RightL;
   }
+  CCTK_ENDLOOP3 (eulerauto_cons_calc_reconstruct_3);
 }
 
 extern "C" void eulerauto_cons_calc_reconstruct_3(CCTK_ARGUMENTS)
@@ -203,5 +196,26 @@ extern "C" void eulerauto_cons_calc_reconstruct_3(CCTK_ARGUMENTS)
   DECLARE_CCTK_ARGUMENTS;
   DECLARE_CCTK_PARAMETERS;
   
+  
+  if (verbose > 1)
+  {
+    CCTK_VInfo(CCTK_THORNSTRING,"Entering eulerauto_cons_calc_reconstruct_3_Body");
+  }
+  
+  if (cctk_iteration % eulerauto_cons_calc_reconstruct_3_calc_every != eulerauto_cons_calc_reconstruct_3_calc_offset)
+  {
+    return;
+  }
+  
+  const char *groups[] = {"EulerAuto::p_group","EulerAuto::p_lr_group","EulerAuto::rho_group","EulerAuto::rho_lr_group","EulerAuto::v1_lr_group","EulerAuto::v2_lr_group","EulerAuto::v3_lr_group","EulerAuto::v_group"};
+  GenericFD_AssertGroupStorage(cctkGH, "eulerauto_cons_calc_reconstruct_3", 8, groups);
+  
+  GenericFD_EnsureStencilFits(cctkGH, "eulerauto_cons_calc_reconstruct_3", 1, 1, 1);
+  
   GenericFD_LoopOverInterior(cctkGH, &eulerauto_cons_calc_reconstruct_3_Body);
+  
+  if (verbose > 1)
+  {
+    CCTK_VInfo(CCTK_THORNSTRING,"Leaving eulerauto_cons_calc_reconstruct_3_Body");
+  }
 }
diff --git a/Examples/EulerAuto/src/eulerauto_cons_calc_rhs_1.cc b/Examples/EulerAuto/src/eulerauto_cons_calc_rhs_1.cc
index 24e51d4..5f38d60 100644
--- a/Examples/EulerAuto/src/eulerauto_cons_calc_rhs_1.cc
+++ b/Examples/EulerAuto/src/eulerauto_cons_calc_rhs_1.cc
@@ -12,6 +12,8 @@
 #include "cctk_Parameters.h"
 #include "GenericFD.h"
 #include "Differencing.h"
+#include "cctk_Loop.h"
+#include "loopcontrol.h"
 
 /* Define macros used in calculations */
 #define INITVALUE (42)
@@ -38,43 +40,28 @@ extern "C" void eulerauto_cons_calc_rhs_1_SelectBCs(CCTK_ARGUMENTS)
   return;
 }
 
-static void eulerauto_cons_calc_rhs_1_Body(cGH const * restrict const cctkGH, int const dir, int const face, CCTK_REAL const normal[3], CCTK_REAL const tangentA[3], CCTK_REAL const tangentB[3], int const min[3], int const max[3], int const n_subblock_gfs, CCTK_REAL * restrict const subblock_gfs[])
+static void eulerauto_cons_calc_rhs_1_Body(cGH const * restrict const cctkGH, int const dir, int const face, CCTK_REAL const normal[3], CCTK_REAL const tangentA[3], CCTK_REAL const tangentB[3], int const imin[3], int const imax[3], int const n_subblock_gfs, CCTK_REAL * restrict const subblock_gfs[])
 {
   DECLARE_CCTK_ARGUMENTS;
   DECLARE_CCTK_PARAMETERS;
   
   
-  /* Declare the variables used for looping over grid points */
-  CCTK_INT i, j, k;
-  // CCTK_INT index = INITVALUE;
-  
   /* Declare finite differencing variables */
   
-  if (verbose > 1)
-  {
-    CCTK_VInfo(CCTK_THORNSTRING,"Entering eulerauto_cons_calc_rhs_1_Body");
-  }
-  
-  if (cctk_iteration % eulerauto_cons_calc_rhs_1_calc_every != eulerauto_cons_calc_rhs_1_calc_offset)
-  {
-    return;
-  }
-  
-  const char *groups[] = {"EulerAuto::Den_flux_group","EulerAuto::Den_grouprhs","EulerAuto::En_flux_group","EulerAuto::En_grouprhs","EulerAuto::S1_flux_group","EulerAuto::S2_flux_group","EulerAuto::S3_flux_group","EulerAuto::S_grouprhs"};
-  GenericFD_AssertGroupStorage(cctkGH, "eulerauto_cons_calc_rhs_1", 8, groups);
-  
-  GenericFD_EnsureStencilFits(cctkGH, "eulerauto_cons_calc_rhs_1", 1, 1, 1);
-  
   /* Include user-supplied include files */
   
   /* Initialise finite differencing variables */
   ptrdiff_t const di = 1;
   ptrdiff_t const dj = CCTK_GFINDEX3D(cctkGH,0,1,0) - CCTK_GFINDEX3D(cctkGH,0,0,0);
   ptrdiff_t const dk = CCTK_GFINDEX3D(cctkGH,0,0,1) - CCTK_GFINDEX3D(cctkGH,0,0,0);
+  ptrdiff_t const cdi = sizeof(CCTK_REAL) * di;
+  ptrdiff_t const cdj = sizeof(CCTK_REAL) * dj;
+  ptrdiff_t const cdk = sizeof(CCTK_REAL) * dk;
   CCTK_REAL const dx = ToReal(CCTK_DELTA_SPACE(0));
   CCTK_REAL const dy = ToReal(CCTK_DELTA_SPACE(1));
   CCTK_REAL const dz = ToReal(CCTK_DELTA_SPACE(2));
   CCTK_REAL const dt = ToReal(CCTK_DELTA_TIME);
+  CCTK_REAL const t = ToReal(cctk_time);
   CCTK_REAL const dxi = INV(dx);
   CCTK_REAL const dyi = INV(dy);
   CCTK_REAL const dzi = INV(dz);
@@ -93,58 +80,64 @@ static void eulerauto_cons_calc_rhs_1_Body(cGH const * restrict const cctkGH, in
   CCTK_REAL const p1ody = INV(dy);
   CCTK_REAL const p1odz = INV(dz);
   
+  /* Assign local copies of arrays functions */
+  
+  
+  
+  /* Calculate temporaries and arrays functions */
+  
+  /* Copy local copies back to grid functions */
+  
   /* Loop over the grid points */
-  for (k = min[2]; k < max[2]; k++)
+  #pragma omp parallel
+  CCTK_LOOP3 (eulerauto_cons_calc_rhs_1,
+    i,j,k, imin[0],imin[1],imin[2], imax[0],imax[1],imax[2],
+    cctk_lsh[0],cctk_lsh[1],cctk_lsh[2])
   {
-    for (j = min[1]; j < max[1]; j++)
-    {
-      for (i = min[0]; i < max[0]; i++)
-      {
-         int  const  index  =  CCTK_GFINDEX3D(cctkGH,i,j,k) ;
-        
-        /* Assign local copies of grid functions */
-        
-        CCTK_REAL DenFluxL = DenFlux[index];
-        CCTK_REAL DenrhsL = Denrhs[index];
-        CCTK_REAL EnFluxL = EnFlux[index];
-        CCTK_REAL EnrhsL = Enrhs[index];
-        CCTK_REAL S1FluxL = S1Flux[index];
-        CCTK_REAL S1rhsL = S1rhs[index];
-        CCTK_REAL S2FluxL = S2Flux[index];
-        CCTK_REAL S2rhsL = S2rhs[index];
-        CCTK_REAL S3FluxL = S3Flux[index];
-        CCTK_REAL S3rhsL = S3rhs[index];
-        
-        
-        /* Include user supplied include files */
-        
-        /* Precompute derivatives */
-        CCTK_REAL const PDplus1DenFlux = PDplus1(&DenFlux[index]);
-        CCTK_REAL const PDplus1EnFlux = PDplus1(&EnFlux[index]);
-        CCTK_REAL const PDplus1S1Flux = PDplus1(&S1Flux[index]);
-        CCTK_REAL const PDplus1S2Flux = PDplus1(&S2Flux[index]);
-        CCTK_REAL const PDplus1S3Flux = PDplus1(&S3Flux[index]);
-        
-        /* Calculate temporaries and grid functions */
-        DenrhsL = DenrhsL - PDplus1DenFlux;
-        
-        EnrhsL = EnrhsL - PDplus1EnFlux;
-        
-        S1rhsL = -PDplus1S1Flux + S1rhsL;
-        
-        S2rhsL = -PDplus1S2Flux + S2rhsL;
-        
-        S3rhsL = -PDplus1S3Flux + S3rhsL;
-        
-        /* Copy local copies back to grid functions */
-        Denrhs[index] = DenrhsL;
-        Enrhs[index] = EnrhsL;
-        S1rhs[index] = S1rhsL;
-        S2rhs[index] = S2rhsL;
-        S3rhs[index] = S3rhsL;
-      }
-    }
+    ptrdiff_t const index = di*i + dj*j + dk*k;
+    
+    /* Assign local copies of grid functions */
+    
+    CCTK_REAL DenFluxL = DenFlux[index];
+    CCTK_REAL DenrhsL = Denrhs[index];
+    CCTK_REAL EnFluxL = EnFlux[index];
+    CCTK_REAL EnrhsL = Enrhs[index];
+    CCTK_REAL S1FluxL = S1Flux[index];
+    CCTK_REAL S1rhsL = S1rhs[index];
+    CCTK_REAL S2FluxL = S2Flux[index];
+    CCTK_REAL S2rhsL = S2rhs[index];
+    CCTK_REAL S3FluxL = S3Flux[index];
+    CCTK_REAL S3rhsL = S3rhs[index];
+    
+    
+    /* Include user supplied include files */
+    
+    /* Precompute derivatives */
+    CCTK_REAL const PDplus1DenFlux = PDplus1(&DenFlux[index]);
+    CCTK_REAL const PDplus1EnFlux = PDplus1(&EnFlux[index]);
+    CCTK_REAL const PDplus1S1Flux = PDplus1(&S1Flux[index]);
+    CCTK_REAL const PDplus1S2Flux = PDplus1(&S2Flux[index]);
+    CCTK_REAL const PDplus1S3Flux = PDplus1(&S3Flux[index]);
+    
+    /* Calculate temporaries and grid functions */
+    DenrhsL = DenrhsL - PDplus1DenFlux;
+    
+    EnrhsL = EnrhsL - PDplus1EnFlux;
+    
+    S1rhsL = S1rhsL - PDplus1S1Flux;
+    
+    S2rhsL = S2rhsL - PDplus1S2Flux;
+    
+    S3rhsL = S3rhsL - PDplus1S3Flux;
+    
+    /* Copy local copies back to grid functions */
+    Denrhs[index] = DenrhsL;
+    Enrhs[index] = EnrhsL;
+    S1rhs[index] = S1rhsL;
+    S2rhs[index] = S2rhsL;
+    S3rhs[index] = S3rhsL;
   }
+  CCTK_ENDLOOP3 (eulerauto_cons_calc_rhs_1);
 }
 
 extern "C" void eulerauto_cons_calc_rhs_1(CCTK_ARGUMENTS)
@@ -152,5 +145,26 @@ extern "C" void eulerauto_cons_calc_rhs_1(CCTK_ARGUMENTS)
   DECLARE_CCTK_ARGUMENTS;
   DECLARE_CCTK_PARAMETERS;
   
+  
+  if (verbose > 1)
+  {
+    CCTK_VInfo(CCTK_THORNSTRING,"Entering eulerauto_cons_calc_rhs_1_Body");
+  }
+  
+  if (cctk_iteration % eulerauto_cons_calc_rhs_1_calc_every != eulerauto_cons_calc_rhs_1_calc_offset)
+  {
+    return;
+  }
+  
+  const char *groups[] = {"EulerAuto::Den_flux_group","EulerAuto::Den_grouprhs","EulerAuto::En_flux_group","EulerAuto::En_grouprhs","EulerAuto::S1_flux_group","EulerAuto::S2_flux_group","EulerAuto::S3_flux_group","EulerAuto::S_grouprhs"};
+  GenericFD_AssertGroupStorage(cctkGH, "eulerauto_cons_calc_rhs_1", 8, groups);
+  
+  GenericFD_EnsureStencilFits(cctkGH, "eulerauto_cons_calc_rhs_1", 1, 1, 1);
+  
   GenericFD_LoopOverInterior(cctkGH, &eulerauto_cons_calc_rhs_1_Body);
+  
+  if (verbose > 1)
+  {
+    CCTK_VInfo(CCTK_THORNSTRING,"Leaving eulerauto_cons_calc_rhs_1_Body");
+  }
 }
diff --git a/Examples/EulerAuto/src/eulerauto_cons_calc_rhs_2.cc b/Examples/EulerAuto/src/eulerauto_cons_calc_rhs_2.cc
index caee229..c4859bc 100644
--- a/Examples/EulerAuto/src/eulerauto_cons_calc_rhs_2.cc
+++ b/Examples/EulerAuto/src/eulerauto_cons_calc_rhs_2.cc
@@ -12,6 +12,8 @@
 #include "cctk_Parameters.h"
 #include "GenericFD.h"
 #include "Differencing.h"
+#include "cctk_Loop.h"
+#include "loopcontrol.h"
 
 /* Define macros used in calculations */
 #define INITVALUE (42)
@@ -38,43 +40,28 @@ extern "C" void eulerauto_cons_calc_rhs_2_SelectBCs(CCTK_ARGUMENTS)
   return;
 }
 
-static void eulerauto_cons_calc_rhs_2_Body(cGH const * restrict const cctkGH, int const dir, int const face, CCTK_REAL const normal[3], CCTK_REAL const tangentA[3], CCTK_REAL const tangentB[3], int const min[3], int const max[3], int const n_subblock_gfs, CCTK_REAL * restrict const subblock_gfs[])
+static void eulerauto_cons_calc_rhs_2_Body(cGH const * restrict const cctkGH, int const dir, int const face, CCTK_REAL const normal[3], CCTK_REAL const tangentA[3], CCTK_REAL const tangentB[3], int const imin[3], int const imax[3], int const n_subblock_gfs, CCTK_REAL * restrict const subblock_gfs[])
 {
   DECLARE_CCTK_ARGUMENTS;
   DECLARE_CCTK_PARAMETERS;
   
   
-  /* Declare the variables used for looping over grid points */
-  CCTK_INT i, j, k;
-  // CCTK_INT index = INITVALUE;
-  
   /* Declare finite differencing variables */
   
-  if (verbose > 1)
-  {
-    CCTK_VInfo(CCTK_THORNSTRING,"Entering eulerauto_cons_calc_rhs_2_Body");
-  }
-  
-  if (cctk_iteration % eulerauto_cons_calc_rhs_2_calc_every != eulerauto_cons_calc_rhs_2_calc_offset)
-  {
-    return;
-  }
-  
-  const char *groups[] = {"EulerAuto::Den_flux_group","EulerAuto::Den_grouprhs","EulerAuto::En_flux_group","EulerAuto::En_grouprhs","EulerAuto::S1_flux_group","EulerAuto::S2_flux_group","EulerAuto::S3_flux_group","EulerAuto::S_grouprhs"};
-  GenericFD_AssertGroupStorage(cctkGH, "eulerauto_cons_calc_rhs_2", 8, groups);
-  
-  GenericFD_EnsureStencilFits(cctkGH, "eulerauto_cons_calc_rhs_2", 1, 1, 1);
-  
   /* Include user-supplied include files */
   
   /* Initialise finite differencing variables */
   ptrdiff_t const di = 1;
   ptrdiff_t const dj = CCTK_GFINDEX3D(cctkGH,0,1,0) - CCTK_GFINDEX3D(cctkGH,0,0,0);
   ptrdiff_t const dk = CCTK_GFINDEX3D(cctkGH,0,0,1) - CCTK_GFINDEX3D(cctkGH,0,0,0);
+  ptrdiff_t const cdi = sizeof(CCTK_REAL) * di;
+  ptrdiff_t const cdj = sizeof(CCTK_REAL) * dj;
+  ptrdiff_t const cdk = sizeof(CCTK_REAL) * dk;
   CCTK_REAL const dx = ToReal(CCTK_DELTA_SPACE(0));
   CCTK_REAL const dy = ToReal(CCTK_DELTA_SPACE(1));
   CCTK_REAL const dz = ToReal(CCTK_DELTA_SPACE(2));
   CCTK_REAL const dt = ToReal(CCTK_DELTA_TIME);
+  CCTK_REAL const t = ToReal(cctk_time);
   CCTK_REAL const dxi = INV(dx);
   CCTK_REAL const dyi = INV(dy);
   CCTK_REAL const dzi = INV(dz);
@@ -93,58 +80,64 @@ static void eulerauto_cons_calc_rhs_2_Body(cGH const * restrict const cctkGH, in
   CCTK_REAL const p1ody = INV(dy);
   CCTK_REAL const p1odz = INV(dz);
   
+  /* Assign local copies of arrays functions */
+  
+  
+  
+  /* Calculate temporaries and arrays functions */
+  
+  /* Copy local copies back to grid functions */
+  
   /* Loop over the grid points */
-  for (k = min[2]; k < max[2]; k++)
+  #pragma omp parallel
+  CCTK_LOOP3 (eulerauto_cons_calc_rhs_2,
+    i,j,k, imin[0],imin[1],imin[2], imax[0],imax[1],imax[2],
+    cctk_lsh[0],cctk_lsh[1],cctk_lsh[2])
   {
-    for (j = min[1]; j < max[1]; j++)
-    {
-      for (i = min[0]; i < max[0]; i++)
-      {
-         int  const  index  =  CCTK_GFINDEX3D(cctkGH,i,j,k) ;
-        
-        /* Assign local copies of grid functions */
-        
-        CCTK_REAL DenFluxL = DenFlux[index];
-        CCTK_REAL DenrhsL = Denrhs[index];
-        CCTK_REAL EnFluxL = EnFlux[index];
-        CCTK_REAL EnrhsL = Enrhs[index];
-        CCTK_REAL S1FluxL = S1Flux[index];
-        CCTK_REAL S1rhsL = S1rhs[index];
-        CCTK_REAL S2FluxL = S2Flux[index];
-        CCTK_REAL S2rhsL = S2rhs[index];
-        CCTK_REAL S3FluxL = S3Flux[index];
-        CCTK_REAL S3rhsL = S3rhs[index];
-        
-        
-        /* Include user supplied include files */
-        
-        /* Precompute derivatives */
-        CCTK_REAL const PDplus2DenFlux = PDplus2(&DenFlux[index]);
-        CCTK_REAL const PDplus2EnFlux = PDplus2(&EnFlux[index]);
-        CCTK_REAL const PDplus2S1Flux = PDplus2(&S1Flux[index]);
-        CCTK_REAL const PDplus2S2Flux = PDplus2(&S2Flux[index]);
-        CCTK_REAL const PDplus2S3Flux = PDplus2(&S3Flux[index]);
-        
-        /* Calculate temporaries and grid functions */
-        DenrhsL = DenrhsL - PDplus2DenFlux;
-        
-        EnrhsL = EnrhsL - PDplus2EnFlux;
-        
-        S1rhsL = -PDplus2S1Flux + S1rhsL;
-        
-        S2rhsL = -PDplus2S2Flux + S2rhsL;
-        
-        S3rhsL = -PDplus2S3Flux + S3rhsL;
-        
-        /* Copy local copies back to grid functions */
-        Denrhs[index] = DenrhsL;
-        Enrhs[index] = EnrhsL;
-        S1rhs[index] = S1rhsL;
-        S2rhs[index] = S2rhsL;
-        S3rhs[index] = S3rhsL;
-      }
-    }
+    ptrdiff_t const index = di*i + dj*j + dk*k;
+    
+    /* Assign local copies of grid functions */
+    
+    CCTK_REAL DenFluxL = DenFlux[index];
+    CCTK_REAL DenrhsL = Denrhs[index];
+    CCTK_REAL EnFluxL = EnFlux[index];
+    CCTK_REAL EnrhsL = Enrhs[index];
+    CCTK_REAL S1FluxL = S1Flux[index];
+    CCTK_REAL S1rhsL = S1rhs[index];
+    CCTK_REAL S2FluxL = S2Flux[index];
+    CCTK_REAL S2rhsL = S2rhs[index];
+    CCTK_REAL S3FluxL = S3Flux[index];
+    CCTK_REAL S3rhsL = S3rhs[index];
+    
+    
+    /* Include user supplied include files */
+    
+    /* Precompute derivatives */
+    CCTK_REAL const PDplus2DenFlux = PDplus2(&DenFlux[index]);
+    CCTK_REAL const PDplus2EnFlux = PDplus2(&EnFlux[index]);
+    CCTK_REAL const PDplus2S1Flux = PDplus2(&S1Flux[index]);
+    CCTK_REAL const PDplus2S2Flux = PDplus2(&S2Flux[index]);
+    CCTK_REAL const PDplus2S3Flux = PDplus2(&S3Flux[index]);
+    
+    /* Calculate temporaries and grid functions */
+    DenrhsL = DenrhsL - PDplus2DenFlux;
+    
+    EnrhsL = EnrhsL - PDplus2EnFlux;
+    
+    S1rhsL = S1rhsL - PDplus2S1Flux;
+    
+    S2rhsL = S2rhsL - PDplus2S2Flux;
+    
+    S3rhsL = S3rhsL - PDplus2S3Flux;
+    
+    /* Copy local copies back to grid functions */
+    Denrhs[index] = DenrhsL;
+    Enrhs[index] = EnrhsL;
+    S1rhs[index] = S1rhsL;
+    S2rhs[index] = S2rhsL;
+    S3rhs[index] = S3rhsL;
   }
+  CCTK_ENDLOOP3 (eulerauto_cons_calc_rhs_2);
 }
 
 extern "C" void eulerauto_cons_calc_rhs_2(CCTK_ARGUMENTS)
@@ -152,5 +145,26 @@ extern "C" void eulerauto_cons_calc_rhs_2(CCTK_ARGUMENTS)
   DECLARE_CCTK_ARGUMENTS;
   DECLARE_CCTK_PARAMETERS;
   
+  
+  if (verbose > 1)
+  {
+    CCTK_VInfo(CCTK_THORNSTRING,"Entering eulerauto_cons_calc_rhs_2_Body");
+  }
+  
+  if (cctk_iteration % eulerauto_cons_calc_rhs_2_calc_every != eulerauto_cons_calc_rhs_2_calc_offset)
+  {
+    return;
+  }
+  
+  const char *groups[] = {"EulerAuto::Den_flux_group","EulerAuto::Den_grouprhs","EulerAuto::En_flux_group","EulerAuto::En_grouprhs","EulerAuto::S1_flux_group","EulerAuto::S2_flux_group","EulerAuto::S3_flux_group","EulerAuto::S_grouprhs"};
+  GenericFD_AssertGroupStorage(cctkGH, "eulerauto_cons_calc_rhs_2", 8, groups);
+  
+  GenericFD_EnsureStencilFits(cctkGH, "eulerauto_cons_calc_rhs_2", 1, 1, 1);
+  
   GenericFD_LoopOverInterior(cctkGH, &eulerauto_cons_calc_rhs_2_Body);
+  
+  if (verbose > 1)
+  {
+    CCTK_VInfo(CCTK_THORNSTRING,"Leaving eulerauto_cons_calc_rhs_2_Body");
+  }
 }
diff --git a/Examples/EulerAuto/src/eulerauto_cons_calc_rhs_3.cc b/Examples/EulerAuto/src/eulerauto_cons_calc_rhs_3.cc
index f485670..0ce65f6 100644
--- a/Examples/EulerAuto/src/eulerauto_cons_calc_rhs_3.cc
+++ b/Examples/EulerAuto/src/eulerauto_cons_calc_rhs_3.cc
@@ -12,6 +12,8 @@
 #include "cctk_Parameters.h"
 #include "GenericFD.h"
 #include "Differencing.h"
+#include "cctk_Loop.h"
+#include "loopcontrol.h"
 
 /* Define macros used in calculations */
 #define INITVALUE (42)
@@ -38,43 +40,28 @@ extern "C" void eulerauto_cons_calc_rhs_3_SelectBCs(CCTK_ARGUMENTS)
   return;
 }
 
-static void eulerauto_cons_calc_rhs_3_Body(cGH const * restrict const cctkGH, int const dir, int const face, CCTK_REAL const normal[3], CCTK_REAL const tangentA[3], CCTK_REAL const tangentB[3], int const min[3], int const max[3], int const n_subblock_gfs, CCTK_REAL * restrict const subblock_gfs[])
+static void eulerauto_cons_calc_rhs_3_Body(cGH const * restrict const cctkGH, int const dir, int const face, CCTK_REAL const normal[3], CCTK_REAL const tangentA[3], CCTK_REAL const tangentB[3], int const imin[3], int const imax[3], int const n_subblock_gfs, CCTK_REAL * restrict const subblock_gfs[])
 {
   DECLARE_CCTK_ARGUMENTS;
   DECLARE_CCTK_PARAMETERS;
   
   
-  /* Declare the variables used for looping over grid points */
-  CCTK_INT i, j, k;
-  // CCTK_INT index = INITVALUE;
-  
   /* Declare finite differencing variables */
   
-  if (verbose > 1)
-  {
-    CCTK_VInfo(CCTK_THORNSTRING,"Entering eulerauto_cons_calc_rhs_3_Body");
-  }
-  
-  if (cctk_iteration % eulerauto_cons_calc_rhs_3_calc_every != eulerauto_cons_calc_rhs_3_calc_offset)
-  {
-    return;
-  }
-  
-  const char *groups[] = {"EulerAuto::Den_flux_group","EulerAuto::Den_grouprhs","EulerAuto::En_flux_group","EulerAuto::En_grouprhs","EulerAuto::S1_flux_group","EulerAuto::S2_flux_group","EulerAuto::S3_flux_group","EulerAuto::S_grouprhs"};
-  GenericFD_AssertGroupStorage(cctkGH, "eulerauto_cons_calc_rhs_3", 8, groups);
-  
-  GenericFD_EnsureStencilFits(cctkGH, "eulerauto_cons_calc_rhs_3", 1, 1, 1);
-  
   /* Include user-supplied include files */
   
   /* Initialise finite differencing variables */
   ptrdiff_t const di = 1;
   ptrdiff_t const dj = CCTK_GFINDEX3D(cctkGH,0,1,0) - CCTK_GFINDEX3D(cctkGH,0,0,0);
   ptrdiff_t const dk = CCTK_GFINDEX3D(cctkGH,0,0,1) - CCTK_GFINDEX3D(cctkGH,0,0,0);
+  ptrdiff_t const cdi = sizeof(CCTK_REAL) * di;
+  ptrdiff_t const cdj = sizeof(CCTK_REAL) * dj;
+  ptrdiff_t const cdk = sizeof(CCTK_REAL) * dk;
   CCTK_REAL const dx = ToReal(CCTK_DELTA_SPACE(0));
   CCTK_REAL const dy = ToReal(CCTK_DELTA_SPACE(1));
   CCTK_REAL const dz = ToReal(CCTK_DELTA_SPACE(2));
   CCTK_REAL const dt = ToReal(CCTK_DELTA_TIME);
+  CCTK_REAL const t = ToReal(cctk_time);
   CCTK_REAL const dxi = INV(dx);
   CCTK_REAL const dyi = INV(dy);
   CCTK_REAL const dzi = INV(dz);
@@ -93,58 +80,64 @@ static void eulerauto_cons_calc_rhs_3_Body(cGH const * restrict const cctkGH, in
   CCTK_REAL const p1ody = INV(dy);
   CCTK_REAL const p1odz = INV(dz);
   
+  /* Assign local copies of arrays functions */
+  
+  
+  
+  /* Calculate temporaries and arrays functions */
+  
+  /* Copy local copies back to grid functions */
+  
   /* Loop over the grid points */
-  for (k = min[2]; k < max[2]; k++)
+  #pragma omp parallel
+  CCTK_LOOP3 (eulerauto_cons_calc_rhs_3,
+    i,j,k, imin[0],imin[1],imin[2], imax[0],imax[1],imax[2],
+    cctk_lsh[0],cctk_lsh[1],cctk_lsh[2])
   {
-    for (j = min[1]; j < max[1]; j++)
-    {
-      for (i = min[0]; i < max[0]; i++)
-      {
-         int  const  index  =  CCTK_GFINDEX3D(cctkGH,i,j,k) ;
-        
-        /* Assign local copies of grid functions */
-        
-        CCTK_REAL DenFluxL = DenFlux[index];
-        CCTK_REAL DenrhsL = Denrhs[index];
-        CCTK_REAL EnFluxL = EnFlux[index];
-        CCTK_REAL EnrhsL = Enrhs[index];
-        CCTK_REAL S1FluxL = S1Flux[index];
-        CCTK_REAL S1rhsL = S1rhs[index];
-        CCTK_REAL S2FluxL = S2Flux[index];
-        CCTK_REAL S2rhsL = S2rhs[index];
-        CCTK_REAL S3FluxL = S3Flux[index];
-        CCTK_REAL S3rhsL = S3rhs[index];
-        
-        
-        /* Include user supplied include files */
-        
-        /* Precompute derivatives */
-        CCTK_REAL const PDplus3DenFlux = PDplus3(&DenFlux[index]);
-        CCTK_REAL const PDplus3EnFlux = PDplus3(&EnFlux[index]);
-        CCTK_REAL const PDplus3S1Flux = PDplus3(&S1Flux[index]);
-        CCTK_REAL const PDplus3S2Flux = PDplus3(&S2Flux[index]);
-        CCTK_REAL const PDplus3S3Flux = PDplus3(&S3Flux[index]);
-        
-        /* Calculate temporaries and grid functions */
-        DenrhsL = DenrhsL - PDplus3DenFlux;
-        
-        EnrhsL = EnrhsL - PDplus3EnFlux;
-        
-        S1rhsL = -PDplus3S1Flux + S1rhsL;
-        
-        S2rhsL = -PDplus3S2Flux + S2rhsL;
-        
-        S3rhsL = -PDplus3S3Flux + S3rhsL;
-        
-        /* Copy local copies back to grid functions */
-        Denrhs[index] = DenrhsL;
-        Enrhs[index] = EnrhsL;
-        S1rhs[index] = S1rhsL;
-        S2rhs[index] = S2rhsL;
-        S3rhs[index] = S3rhsL;
-      }
-    }
+    ptrdiff_t const index = di*i + dj*j + dk*k;
+    
+    /* Assign local copies of grid functions */
+    
+    CCTK_REAL DenFluxL = DenFlux[index];
+    CCTK_REAL DenrhsL = Denrhs[index];
+    CCTK_REAL EnFluxL = EnFlux[index];
+    CCTK_REAL EnrhsL = Enrhs[index];
+    CCTK_REAL S1FluxL = S1Flux[index];
+    CCTK_REAL S1rhsL = S1rhs[index];
+    CCTK_REAL S2FluxL = S2Flux[index];
+    CCTK_REAL S2rhsL = S2rhs[index];
+    CCTK_REAL S3FluxL = S3Flux[index];
+    CCTK_REAL S3rhsL = S3rhs[index];
+    
+    
+    /* Include user supplied include files */
+    
+    /* Precompute derivatives */
+    CCTK_REAL const PDplus3DenFlux = PDplus3(&DenFlux[index]);
+    CCTK_REAL const PDplus3EnFlux = PDplus3(&EnFlux[index]);
+    CCTK_REAL const PDplus3S1Flux = PDplus3(&S1Flux[index]);
+    CCTK_REAL const PDplus3S2Flux = PDplus3(&S2Flux[index]);
+    CCTK_REAL const PDplus3S3Flux = PDplus3(&S3Flux[index]);
+    
+    /* Calculate temporaries and grid functions */
+    DenrhsL = DenrhsL - PDplus3DenFlux;
+    
+    EnrhsL = EnrhsL - PDplus3EnFlux;
+    
+    S1rhsL = S1rhsL - PDplus3S1Flux;
+    
+    S2rhsL = S2rhsL - PDplus3S2Flux;
+    
+    S3rhsL = S3rhsL - PDplus3S3Flux;
+    
+    /* Copy local copies back to grid functions */
+    Denrhs[index] = DenrhsL;
+    Enrhs[index] = EnrhsL;
+    S1rhs[index] = S1rhsL;
+    S2rhs[index] = S2rhsL;
+    S3rhs[index] = S3rhsL;
   }
+  CCTK_ENDLOOP3 (eulerauto_cons_calc_rhs_3);
 }
 
 extern "C" void eulerauto_cons_calc_rhs_3(CCTK_ARGUMENTS)
@@ -152,5 +145,26 @@ extern "C" void eulerauto_cons_calc_rhs_3(CCTK_ARGUMENTS)
   DECLARE_CCTK_ARGUMENTS;
   DECLARE_CCTK_PARAMETERS;
   
+  
+  if (verbose > 1)
+  {
+    CCTK_VInfo(CCTK_THORNSTRING,"Entering eulerauto_cons_calc_rhs_3_Body");
+  }
+  
+  if (cctk_iteration % eulerauto_cons_calc_rhs_3_calc_every != eulerauto_cons_calc_rhs_3_calc_offset)
+  {
+    return;
+  }
+  
+  const char *groups[] = {"EulerAuto::Den_flux_group","EulerAuto::Den_grouprhs","EulerAuto::En_flux_group","EulerAuto::En_grouprhs","EulerAuto::S1_flux_group","EulerAuto::S2_flux_group","EulerAuto::S3_flux_group","EulerAuto::S_grouprhs"};
+  GenericFD_AssertGroupStorage(cctkGH, "eulerauto_cons_calc_rhs_3", 8, groups);
+  
+  GenericFD_EnsureStencilFits(cctkGH, "eulerauto_cons_calc_rhs_3", 1, 1, 1);
+  
   GenericFD_LoopOverInterior(cctkGH, &eulerauto_cons_calc_rhs_3_Body);
+  
+  if (verbose > 1)
+  {
+    CCTK_VInfo(CCTK_THORNSTRING,"Leaving eulerauto_cons_calc_rhs_3_Body");
+  }
 }
diff --git a/Examples/EulerAuto/src/eulerauto_cons_calc_zero_rhs.cc b/Examples/EulerAuto/src/eulerauto_cons_calc_zero_rhs.cc
index ba14718..5418435 100644
--- a/Examples/EulerAuto/src/eulerauto_cons_calc_zero_rhs.cc
+++ b/Examples/EulerAuto/src/eulerauto_cons_calc_zero_rhs.cc
@@ -12,6 +12,8 @@
 #include "cctk_Parameters.h"
 #include "GenericFD.h"
 #include "Differencing.h"
+#include "cctk_Loop.h"
+#include "loopcontrol.h"
 
 /* Define macros used in calculations */
 #define INITVALUE (42)
@@ -20,42 +22,28 @@
 #define SQR(x) ((x) * (x))
 #define CUB(x) ((x) * (x) * (x))
 
-static void eulerauto_cons_calc_zero_rhs_Body(cGH const * restrict const cctkGH, int const dir, int const face, CCTK_REAL const normal[3], CCTK_REAL const tangentA[3], CCTK_REAL const tangentB[3], int const min[3], int const max[3], int const n_subblock_gfs, CCTK_REAL * restrict const subblock_gfs[])
+static void eulerauto_cons_calc_zero_rhs_Body(cGH const * restrict const cctkGH, int const dir, int const face, CCTK_REAL const normal[3], CCTK_REAL const tangentA[3], CCTK_REAL const tangentB[3], int const imin[3], int const imax[3], int const n_subblock_gfs, CCTK_REAL * restrict const subblock_gfs[])
 {
   DECLARE_CCTK_ARGUMENTS;
   DECLARE_CCTK_PARAMETERS;
   
   
-  /* Declare the variables used for looping over grid points */
-  CCTK_INT i, j, k;
-  // CCTK_INT index = INITVALUE;
-  
   /* Declare finite differencing variables */
   
-  if (verbose > 1)
-  {
-    CCTK_VInfo(CCTK_THORNSTRING,"Entering eulerauto_cons_calc_zero_rhs_Body");
-  }
-  
-  if (cctk_iteration % eulerauto_cons_calc_zero_rhs_calc_every != eulerauto_cons_calc_zero_rhs_calc_offset)
-  {
-    return;
-  }
-  
-  const char *groups[] = {"EulerAuto::Den_grouprhs","EulerAuto::En_grouprhs","EulerAuto::S_grouprhs"};
-  GenericFD_AssertGroupStorage(cctkGH, "eulerauto_cons_calc_zero_rhs", 3, groups);
-  
-  
   /* Include user-supplied include files */
   
   /* Initialise finite differencing variables */
   ptrdiff_t const di = 1;
   ptrdiff_t const dj = CCTK_GFINDEX3D(cctkGH,0,1,0) - CCTK_GFINDEX3D(cctkGH,0,0,0);
   ptrdiff_t const dk = CCTK_GFINDEX3D(cctkGH,0,0,1) - CCTK_GFINDEX3D(cctkGH,0,0,0);
+  ptrdiff_t const cdi = sizeof(CCTK_REAL) * di;
+  ptrdiff_t const cdj = sizeof(CCTK_REAL) * dj;
+  ptrdiff_t const cdk = sizeof(CCTK_REAL) * dk;
   CCTK_REAL const dx = ToReal(CCTK_DELTA_SPACE(0));
   CCTK_REAL const dy = ToReal(CCTK_DELTA_SPACE(1));
   CCTK_REAL const dz = ToReal(CCTK_DELTA_SPACE(2));
   CCTK_REAL const dt = ToReal(CCTK_DELTA_TIME);
+  CCTK_REAL const t = ToReal(cctk_time);
   CCTK_REAL const dxi = INV(dx);
   CCTK_REAL const dyi = INV(dy);
   CCTK_REAL const dzi = INV(dz);
@@ -74,63 +62,69 @@ static void eulerauto_cons_calc_zero_rhs_Body(cGH const * restrict const cctkGH,
   CCTK_REAL const p1ody = INV(dy);
   CCTK_REAL const p1odz = INV(dz);
   
+  /* Assign local copies of arrays functions */
+  
+  
+  
+  /* Calculate temporaries and arrays functions */
+  
+  /* Copy local copies back to grid functions */
+  
   /* Loop over the grid points */
-  for (k = min[2]; k < max[2]; k++)
+  #pragma omp parallel
+  CCTK_LOOP3 (eulerauto_cons_calc_zero_rhs,
+    i,j,k, imin[0],imin[1],imin[2], imax[0],imax[1],imax[2],
+    cctk_lsh[0],cctk_lsh[1],cctk_lsh[2])
   {
-    for (j = min[1]; j < max[1]; j++)
-    {
-      for (i = min[0]; i < max[0]; i++)
-      {
-         int  const  index  =  CCTK_GFINDEX3D(cctkGH,i,j,k) ;
-        
-        /* Assign local copies of grid functions */
-        
-        
-        
-        /* Include user supplied include files */
-        
-        /* Precompute derivatives */
-        
-        /* Calculate temporaries and grid functions */
-        CCTK_REAL DenrhsL = 0;
-        
-        DenrhsL = 0;
-        
-        DenrhsL = 0;
-        
-        CCTK_REAL S1rhsL = 0;
-        
-        CCTK_REAL S2rhsL = 0;
-        
-        CCTK_REAL S3rhsL = 0;
-        
-        S1rhsL = 0;
-        
-        S2rhsL = 0;
-        
-        S3rhsL = 0;
-        
-        S1rhsL = 0;
-        
-        S2rhsL = 0;
-        
-        S3rhsL = 0;
-        
-        CCTK_REAL EnrhsL = 0;
-        
-        EnrhsL = 0;
-        
-        EnrhsL = 0;
-        
-        /* Copy local copies back to grid functions */
-        Denrhs[index] = DenrhsL;
-        Enrhs[index] = EnrhsL;
-        S1rhs[index] = S1rhsL;
-        S2rhs[index] = S2rhsL;
-        S3rhs[index] = S3rhsL;
-      }
-    }
+    ptrdiff_t const index = di*i + dj*j + dk*k;
+    
+    /* Assign local copies of grid functions */
+    
+    
+    
+    /* Include user supplied include files */
+    
+    /* Precompute derivatives */
+    
+    /* Calculate temporaries and grid functions */
+    CCTK_REAL DenrhsL = 0;
+    
+    DenrhsL = 0;
+    
+    DenrhsL = 0;
+    
+    CCTK_REAL S1rhsL = 0;
+    
+    CCTK_REAL S2rhsL = 0;
+    
+    CCTK_REAL S3rhsL = 0;
+    
+    S1rhsL = 0;
+    
+    S2rhsL = 0;
+    
+    S3rhsL = 0;
+    
+    S1rhsL = 0;
+    
+    S2rhsL = 0;
+    
+    S3rhsL = 0;
+    
+    CCTK_REAL EnrhsL = 0;
+    
+    EnrhsL = 0;
+    
+    EnrhsL = 0;
+    
+    /* Copy local copies back to grid functions */
+    Denrhs[index] = DenrhsL;
+    Enrhs[index] = EnrhsL;
+    S1rhs[index] = S1rhsL;
+    S2rhs[index] = S2rhsL;
+    S3rhs[index] = S3rhsL;
   }
+  CCTK_ENDLOOP3 (eulerauto_cons_calc_zero_rhs);
 }
 
 extern "C" void eulerauto_cons_calc_zero_rhs(CCTK_ARGUMENTS)
@@ -138,5 +132,25 @@ extern "C" void eulerauto_cons_calc_zero_rhs(CCTK_ARGUMENTS)
   DECLARE_CCTK_ARGUMENTS;
   DECLARE_CCTK_PARAMETERS;
   
+  
+  if (verbose > 1)
+  {
+    CCTK_VInfo(CCTK_THORNSTRING,"Entering eulerauto_cons_calc_zero_rhs_Body");
+  }
+  
+  if (cctk_iteration % eulerauto_cons_calc_zero_rhs_calc_every != eulerauto_cons_calc_zero_rhs_calc_offset)
+  {
+    return;
+  }
+  
+  const char *groups[] = {"EulerAuto::Den_grouprhs","EulerAuto::En_grouprhs","EulerAuto::S_grouprhs"};
+  GenericFD_AssertGroupStorage(cctkGH, "eulerauto_cons_calc_zero_rhs", 3, groups);
+  
+  
   GenericFD_LoopOverEverything(cctkGH, &eulerauto_cons_calc_zero_rhs_Body);
+  
+  if (verbose > 1)
+  {
+    CCTK_VInfo(CCTK_THORNSTRING,"Leaving eulerauto_cons_calc_zero_rhs_Body");
+  }
 }
diff --git a/Examples/EulerAuto/src/eulerauto_initial_shock.cc b/Examples/EulerAuto/src/eulerauto_initial_shock.cc
index f6893b7..ea1b342 100644
--- a/Examples/EulerAuto/src/eulerauto_initial_shock.cc
+++ b/Examples/EulerAuto/src/eulerauto_initial_shock.cc
@@ -12,6 +12,8 @@
 #include "cctk_Parameters.h"
 #include "GenericFD.h"
 #include "Differencing.h"
+#include "cctk_Loop.h"
+#include "loopcontrol.h"
 
 /* Define macros used in calculations */
 #define INITVALUE (42)
@@ -20,42 +22,28 @@
 #define SQR(x) ((x) * (x))
 #define CUB(x) ((x) * (x) * (x))
 
-static void eulerauto_initial_shock_Body(cGH const * restrict const cctkGH, int const dir, int const face, CCTK_REAL const normal[3], CCTK_REAL const tangentA[3], CCTK_REAL const tangentB[3], int const min[3], int const max[3], int const n_subblock_gfs, CCTK_REAL * restrict const subblock_gfs[])
+static void eulerauto_initial_shock_Body(cGH const * restrict const cctkGH, int const dir, int const face, CCTK_REAL const normal[3], CCTK_REAL const tangentA[3], CCTK_REAL const tangentB[3], int const imin[3], int const imax[3], int const n_subblock_gfs, CCTK_REAL * restrict const subblock_gfs[])
 {
   DECLARE_CCTK_ARGUMENTS;
   DECLARE_CCTK_PARAMETERS;
   
   
-  /* Declare the variables used for looping over grid points */
-  CCTK_INT i, j, k;
-  // CCTK_INT index = INITVALUE;
-  
   /* Declare finite differencing variables */
   
-  if (verbose > 1)
-  {
-    CCTK_VInfo(CCTK_THORNSTRING,"Entering eulerauto_initial_shock_Body");
-  }
-  
-  if (cctk_iteration % eulerauto_initial_shock_calc_every != eulerauto_initial_shock_calc_offset)
-  {
-    return;
-  }
-  
-  const char *groups[] = {"grid::coordinates","EulerAuto::p_group","EulerAuto::rho_group","EulerAuto::v_group"};
-  GenericFD_AssertGroupStorage(cctkGH, "eulerauto_initial_shock", 4, groups);
-  
-  
   /* Include user-supplied include files */
   
   /* Initialise finite differencing variables */
   ptrdiff_t const di = 1;
   ptrdiff_t const dj = CCTK_GFINDEX3D(cctkGH,0,1,0) - CCTK_GFINDEX3D(cctkGH,0,0,0);
   ptrdiff_t const dk = CCTK_GFINDEX3D(cctkGH,0,0,1) - CCTK_GFINDEX3D(cctkGH,0,0,0);
+  ptrdiff_t const cdi = sizeof(CCTK_REAL) * di;
+  ptrdiff_t const cdj = sizeof(CCTK_REAL) * dj;
+  ptrdiff_t const cdk = sizeof(CCTK_REAL) * dk;
   CCTK_REAL const dx = ToReal(CCTK_DELTA_SPACE(0));
   CCTK_REAL const dy = ToReal(CCTK_DELTA_SPACE(1));
   CCTK_REAL const dz = ToReal(CCTK_DELTA_SPACE(2));
   CCTK_REAL const dt = ToReal(CCTK_DELTA_TIME);
+  CCTK_REAL const t = ToReal(cctk_time);
   CCTK_REAL const dxi = INV(dx);
   CCTK_REAL const dyi = INV(dy);
   CCTK_REAL const dzi = INV(dz);
@@ -74,47 +62,53 @@ static void eulerauto_initial_shock_Body(cGH const * restrict const cctkGH, int
   CCTK_REAL const p1ody = INV(dy);
   CCTK_REAL const p1odz = INV(dz);
   
+  /* Assign local copies of arrays functions */
+  
+  
+  
+  /* Calculate temporaries and arrays functions */
+  
+  /* Copy local copies back to grid functions */
+  
   /* Loop over the grid points */
-  for (k = min[2]; k < max[2]; k++)
+  #pragma omp parallel
+  CCTK_LOOP3 (eulerauto_initial_shock,
+    i,j,k, imin[0],imin[1],imin[2], imax[0],imax[1],imax[2],
+    cctk_lsh[0],cctk_lsh[1],cctk_lsh[2])
   {
-    for (j = min[1]; j < max[1]; j++)
-    {
-      for (i = min[0]; i < max[0]; i++)
-      {
-         int  const  index  =  CCTK_GFINDEX3D(cctkGH,i,j,k) ;
-        
-        /* Assign local copies of grid functions */
-        
-        CCTK_REAL xL = x[index];
-        
-        
-        /* Include user supplied include files */
-        
-        /* Precompute derivatives */
-        
-        /* Calculate temporaries and grid functions */
-        CCTK_REAL rhoL = ToReal(rhoL0) + StepFunction(-0.5 + 
-          xL)*(-ToReal(rhoL0) + ToReal(rhoR0));
-        
-        CCTK_REAL v1L = ToReal(vL0) + StepFunction(-0.5 + xL)*(-ToReal(vL0) + 
-          ToReal(vR0));
-        
-        CCTK_REAL v2L = 0;
-        
-        CCTK_REAL v3L = 0;
-        
-        CCTK_REAL pL = ToReal(pL0) + StepFunction(-0.5 + xL)*(-ToReal(pL0) + 
-          ToReal(pR0));
-        
-        /* Copy local copies back to grid functions */
-        p[index] = pL;
-        rho[index] = rhoL;
-        v1[index] = v1L;
-        v2[index] = v2L;
-        v3[index] = v3L;
-      }
-    }
+    ptrdiff_t const index = di*i + dj*j + dk*k;
+    
+    /* Assign local copies of grid functions */
+    
+    CCTK_REAL xL = x[index];
+    
+    
+    /* Include user supplied include files */
+    
+    /* Precompute derivatives */
+    
+    /* Calculate temporaries and grid functions */
+    CCTK_REAL rhoL = ToReal(rhoL0) + StepFunction(-0.5 + 
+      xL)*(-ToReal(rhoL0) + ToReal(rhoR0));
+    
+    CCTK_REAL v1L = ToReal(vL0) + StepFunction(-0.5 + xL)*(-ToReal(vL0) 
+      + ToReal(vR0));
+    
+    CCTK_REAL v2L = 0;
+    
+    CCTK_REAL v3L = 0;
+    
+    CCTK_REAL pL = ToReal(pL0) + StepFunction(-0.5 + xL)*(-ToReal(pL0) + 
+      ToReal(pR0));
+    
+    /* Copy local copies back to grid functions */
+    p[index] = pL;
+    rho[index] = rhoL;
+    v1[index] = v1L;
+    v2[index] = v2L;
+    v3[index] = v3L;
   }
+  CCTK_ENDLOOP3 (eulerauto_initial_shock);
 }
 
 extern "C" void eulerauto_initial_shock(CCTK_ARGUMENTS)
@@ -122,5 +116,25 @@ extern "C" void eulerauto_initial_shock(CCTK_ARGUMENTS)
   DECLARE_CCTK_ARGUMENTS;
   DECLARE_CCTK_PARAMETERS;
   
+  
+  if (verbose > 1)
+  {
+    CCTK_VInfo(CCTK_THORNSTRING,"Entering eulerauto_initial_shock_Body");
+  }
+  
+  if (cctk_iteration % eulerauto_initial_shock_calc_every != eulerauto_initial_shock_calc_offset)
+  {
+    return;
+  }
+  
+  const char *groups[] = {"grid::coordinates","EulerAuto::p_group","EulerAuto::rho_group","EulerAuto::v_group"};
+  GenericFD_AssertGroupStorage(cctkGH, "eulerauto_initial_shock", 4, groups);
+  
+  
   GenericFD_LoopOverEverything(cctkGH, &eulerauto_initial_shock_Body);
+  
+  if (verbose > 1)
+  {
+    CCTK_VInfo(CCTK_THORNSTRING,"Leaving eulerauto_initial_shock_Body");
+  }
 }