Regenerate examples (except Advect)

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

22 files changed, 2052 insertions, 1640 deletions

diff --git a/Examples/EulerSR/configuration.ccl b/Examples/EulerSR/configuration.ccl
index 023aac7..0a66ec2 100644
--- a/Examples/EulerSR/configuration.ccl
+++ b/Examples/EulerSR/configuration.ccl
@@ -1,3 +1,6 @@
 # File produced by Kranc
 
 REQUIRES GenericFD
+OPTIONAL LoopControl
+{
+}
diff --git a/Examples/EulerSR/interface.ccl b/Examples/EulerSR/interface.ccl
index eb300ee..41c6ab6 100644
--- a/Examples/EulerSR/interface.ccl
+++ b/Examples/EulerSR/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/EulerSR/param.ccl b/Examples/EulerSR/param.ccl
index 211a7e7..87f07c5 100644
--- a/Examples/EulerSR/param.ccl
+++ b/Examples/EulerSR/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 EulerSR_MaxNumEvolvedVars "Number of evolved variables used by this tho
 } 5
 
 restricted:
+CCTK_INT EulerSR_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 eulersr_initial_shock_calc_every "eulersr_initial_shock_calc_every" STEERABLE=ALWAYS
 {
   *:* :: ""
diff --git a/Examples/EulerSR/schedule.ccl b/Examples/EulerSR/schedule.ccl
index ffeaf07..d6b29e0 100644
--- a/Examples/EulerSR/schedule.ccl
+++ b/Examples/EulerSR/schedule.ccl
@@ -1,47 +1,110 @@
 # 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: epsi_group[1]
+if (other_timelevels == 1)
+{
+  STORAGE: epsi_group[1]
+}
 
-STORAGE: epsi_lr_group[1]
+if (other_timelevels == 1)
+{
+  STORAGE: epsi_lr_group[1]
+}
 
-STORAGE: h_group[1]
+if (other_timelevels == 1)
+{
+  STORAGE: h_group[1]
+}
 
-STORAGE: p_group[1]
+if (other_timelevels == 1)
+{
+  STORAGE: p_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: tau_flux_group[1]
+if (other_timelevels == 1)
+{
+  STORAGE: tau_flux_group[1]
+}
 
-STORAGE: tau_lr_group[1]
+if (other_timelevels == 1)
+{
+  STORAGE: tau_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]
+}
 
-STORAGE: W_group[1]
+if (other_timelevels == 1)
+{
+  STORAGE: W_group[1]
+}
 
 if (timelevels == 1)
 {
@@ -194,12 +257,7 @@ schedule group eulersr_cons_calc_reconstruct_1_bc_group in eulersr_cons_calc_rec
   # no language specified
 } "eulersr_cons_calc_reconstruct_1"
 
-schedule group eulersr_cons_calc_reconstruct_1_bc_group at CCTK_POSTRESTRICT
-{
-  # no language specified
-} "eulersr_cons_calc_reconstruct_1"
-
-schedule group eulersr_cons_calc_reconstruct_1_bc_group at CCTK_POSTRESTRICTINITIAL
+schedule group eulersr_cons_calc_reconstruct_1_bc_group in MoL_PseudoEvolutionBoundaries after MoL_PostStep
 {
   # no language specified
 } "eulersr_cons_calc_reconstruct_1"
@@ -240,12 +298,7 @@ schedule group eulersr_cons_calc_flux_1_bc_group in eulersr_cons_calc_flux_1_gro
   # no language specified
 } "eulersr_cons_calc_flux_1"
 
-schedule group eulersr_cons_calc_flux_1_bc_group at CCTK_POSTRESTRICT
-{
-  # no language specified
-} "eulersr_cons_calc_flux_1"
-
-schedule group eulersr_cons_calc_flux_1_bc_group at CCTK_POSTRESTRICTINITIAL
+schedule group eulersr_cons_calc_flux_1_bc_group in MoL_PseudoEvolutionBoundaries after MoL_PostStep
 {
   # no language specified
 } "eulersr_cons_calc_flux_1"
@@ -286,12 +339,7 @@ schedule group eulersr_cons_calc_reconstruct_2_bc_group in eulersr_cons_calc_rec
   # no language specified
 } "eulersr_cons_calc_reconstruct_2"
 
-schedule group eulersr_cons_calc_reconstruct_2_bc_group at CCTK_POSTRESTRICT
-{
-  # no language specified
-} "eulersr_cons_calc_reconstruct_2"
-
-schedule group eulersr_cons_calc_reconstruct_2_bc_group at CCTK_POSTRESTRICTINITIAL
+schedule group eulersr_cons_calc_reconstruct_2_bc_group in MoL_PseudoEvolutionBoundaries after MoL_PostStep
 {
   # no language specified
 } "eulersr_cons_calc_reconstruct_2"
@@ -332,12 +380,7 @@ schedule group eulersr_cons_calc_flux_2_bc_group in eulersr_cons_calc_flux_2_gro
   # no language specified
 } "eulersr_cons_calc_flux_2"
 
-schedule group eulersr_cons_calc_flux_2_bc_group at CCTK_POSTRESTRICT
-{
-  # no language specified
-} "eulersr_cons_calc_flux_2"
-
-schedule group eulersr_cons_calc_flux_2_bc_group at CCTK_POSTRESTRICTINITIAL
+schedule group eulersr_cons_calc_flux_2_bc_group in MoL_PseudoEvolutionBoundaries after MoL_PostStep
 {
   # no language specified
 } "eulersr_cons_calc_flux_2"
@@ -378,12 +421,7 @@ schedule group eulersr_cons_calc_reconstruct_3_bc_group in eulersr_cons_calc_rec
   # no language specified
 } "eulersr_cons_calc_reconstruct_3"
 
-schedule group eulersr_cons_calc_reconstruct_3_bc_group at CCTK_POSTRESTRICT
-{
-  # no language specified
-} "eulersr_cons_calc_reconstruct_3"
-
-schedule group eulersr_cons_calc_reconstruct_3_bc_group at CCTK_POSTRESTRICTINITIAL
+schedule group eulersr_cons_calc_reconstruct_3_bc_group in MoL_PseudoEvolutionBoundaries after MoL_PostStep
 {
   # no language specified
 } "eulersr_cons_calc_reconstruct_3"
@@ -424,12 +462,7 @@ schedule group eulersr_cons_calc_flux_3_bc_group in eulersr_cons_calc_flux_3_gro
   # no language specified
 } "eulersr_cons_calc_flux_3"
 
-schedule group eulersr_cons_calc_flux_3_bc_group at CCTK_POSTRESTRICT
-{
-  # no language specified
-} "eulersr_cons_calc_flux_3"
-
-schedule group eulersr_cons_calc_flux_3_bc_group at CCTK_POSTRESTRICTINITIAL
+schedule group eulersr_cons_calc_flux_3_bc_group in MoL_PseudoEvolutionBoundaries after MoL_PostStep
 {
   # no language specified
 } "eulersr_cons_calc_flux_3"
diff --git a/Examples/EulerSR/src/Differencing.h b/Examples/EulerSR/src/Differencing.h
index 0908172..913cfa8 100644
--- a/Examples/EulerSR/src/Differencing.h
+++ b/Examples/EulerSR/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/EulerSR/src/RegisterMoL.cc b/Examples/EulerSR/src/RegisterMoL.cc
index 5f77b35..9ef9a4d 100644
--- a/Examples/EulerSR/src/RegisterMoL.cc
+++ b/Examples/EulerSR/src/RegisterMoL.cc
@@ -17,5 +17,7 @@ extern "C" void EulerSR_RegisterVars(CCTK_ARGUMENTS)
   ierr += MoLRegisterEvolved(CCTK_VarIndex("EulerSR::S2"),  CCTK_VarIndex("EulerSR::S2rhs"));
   ierr += MoLRegisterEvolved(CCTK_VarIndex("EulerSR::S3"),  CCTK_VarIndex("EulerSR::S3rhs"));
   ierr += MoLRegisterEvolved(CCTK_VarIndex("EulerSR::tau"),  CCTK_VarIndex("EulerSR::taurhs"));
+  
+  /* Register all the evolved Array functions with MoL */
   return;
 }
diff --git a/Examples/EulerSR/src/eulersr_cons_calc_conserved.cc b/Examples/EulerSR/src/eulersr_cons_calc_conserved.cc
index 788478c..673ff99 100644
--- a/Examples/EulerSR/src/eulersr_cons_calc_conserved.cc
+++ b/Examples/EulerSR/src/eulersr_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 eulersr_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 eulersr_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 eulersr_cons_calc_conserved_Body");
-  }
-  
-  if (cctk_iteration % eulersr_cons_calc_conserved_calc_every != eulersr_cons_calc_conserved_calc_offset)
-  {
-    return;
-  }
-  
-  const char *groups[] = {"EulerSR::Den_group","EulerSR::epsi_group","EulerSR::h_group","EulerSR::p_group","EulerSR::rho_group","EulerSR::S_group","EulerSR::tau_group","EulerSR::v_group","EulerSR::W_group"};
-  GenericFD_AssertGroupStorage(cctkGH, "eulersr_cons_calc_conserved", 9, 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 eulersr_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 (eulersr_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 DenL = Den[index];
-        CCTK_REAL epsiL = epsi[index];
-        CCTK_REAL hL = h[index];
-        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];
-        CCTK_REAL WL = W[index];
-        
-        
-        /* Include user supplied include files */
-        
-        /* Precompute derivatives */
-        
-        /* Calculate temporaries and grid functions */
-        CCTK_REAL Wx = 1 - SQR(v1L) - SQR(v2L) - SQR(v3L);
-        
-        WL = INV(sqrt(Wx));
-        
-        pL = epsiL*rhoL*(-1 + ToReal(gamma));
-        
-        hL = 1 + epsiL + pL*INV(rhoL);
-        
-        DenL = rhoL*WL;
-        
-        CCTK_REAL S1L = hL*rhoL*v1L*SQR(WL);
-        
-        CCTK_REAL S2L = hL*rhoL*v2L*SQR(WL);
-        
-        CCTK_REAL S3L = hL*rhoL*v3L*SQR(WL);
-        
-        CCTK_REAL tauL = -DenL - pL + hL*rhoL*SQR(WL);
-        
-        /* Copy local copies back to grid functions */
-        Den[index] = DenL;
-        h[index] = hL;
-        p[index] = pL;
-        S1[index] = S1L;
-        S2[index] = S2L;
-        S3[index] = S3L;
-        tau[index] = tauL;
-        W[index] = WL;
-      }
-    }
+    ptrdiff_t const index = di*i + dj*j + dk*k;
+    
+    /* Assign local copies of grid functions */
+    
+    CCTK_REAL DenL = Den[index];
+    CCTK_REAL epsiL = epsi[index];
+    CCTK_REAL hL = h[index];
+    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];
+    CCTK_REAL WL = W[index];
+    
+    
+    /* Include user supplied include files */
+    
+    /* Precompute derivatives */
+    
+    /* Calculate temporaries and grid functions */
+    CCTK_REAL Wx = 1 - SQR(v1L) - SQR(v2L) - SQR(v3L);
+    
+    WL = INV(sqrt(Wx));
+    
+    pL = epsiL*rhoL*(-1 + ToReal(gamma));
+    
+    hL = 1 + epsiL + pL*INV(rhoL);
+    
+    DenL = rhoL*WL;
+    
+    CCTK_REAL S1L = hL*rhoL*v1L*SQR(WL);
+    
+    CCTK_REAL S2L = hL*rhoL*v2L*SQR(WL);
+    
+    CCTK_REAL S3L = hL*rhoL*v3L*SQR(WL);
+    
+    CCTK_REAL tauL = -DenL - pL + hL*rhoL*SQR(WL);
+    
+    /* Copy local copies back to grid functions */
+    Den[index] = DenL;
+    h[index] = hL;
+    p[index] = pL;
+    S1[index] = S1L;
+    S2[index] = S2L;
+    S3[index] = S3L;
+    tau[index] = tauL;
+    W[index] = WL;
   }
+  CCTK_ENDLOOP3 (eulersr_cons_calc_conserved);
 }
 
 extern "C" void eulersr_cons_calc_conserved(CCTK_ARGUMENTS)
@@ -138,5 +132,25 @@ extern "C" void eulersr_cons_calc_conserved(CCTK_ARGUMENTS)
   DECLARE_CCTK_ARGUMENTS;
   DECLARE_CCTK_PARAMETERS;
   
+  
+  if (verbose > 1)
+  {
+    CCTK_VInfo(CCTK_THORNSTRING,"Entering eulersr_cons_calc_conserved_Body");
+  }
+  
+  if (cctk_iteration % eulersr_cons_calc_conserved_calc_every != eulersr_cons_calc_conserved_calc_offset)
+  {
+    return;
+  }
+  
+  const char *groups[] = {"EulerSR::Den_group","EulerSR::epsi_group","EulerSR::h_group","EulerSR::p_group","EulerSR::rho_group","EulerSR::S_group","EulerSR::tau_group","EulerSR::v_group","EulerSR::W_group"};
+  GenericFD_AssertGroupStorage(cctkGH, "eulersr_cons_calc_conserved", 9, groups);
+  
+  
   GenericFD_LoopOverEverything(cctkGH, &eulersr_cons_calc_conserved_Body);
+  
+  if (verbose > 1)
+  {
+    CCTK_VInfo(CCTK_THORNSTRING,"Leaving eulersr_cons_calc_conserved_Body");
+  }
 }
diff --git a/Examples/EulerSR/src/eulersr_cons_calc_flux_1.cc b/Examples/EulerSR/src/eulersr_cons_calc_flux_1.cc
index eefd139..843de50 100644
--- a/Examples/EulerSR/src/eulersr_cons_calc_flux_1.cc
+++ b/Examples/EulerSR/src/eulersr_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 eulersr_cons_calc_flux_1_SelectBCs(CCTK_ARGUMENTS)
   return;
 }
 
-static void eulersr_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 eulersr_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 eulersr_cons_calc_flux_1_Body");
-  }
-  
-  if (cctk_iteration % eulersr_cons_calc_flux_1_calc_every != eulersr_cons_calc_flux_1_calc_offset)
-  {
-    return;
-  }
-  
-  const char *groups[] = {"EulerSR::Den_flux_group","EulerSR::Den_lr_group","EulerSR::epsi_lr_group","EulerSR::rho_lr_group","EulerSR::S1_flux_group","EulerSR::S1_lr_group","EulerSR::S2_flux_group","EulerSR::S2_lr_group","EulerSR::S3_flux_group","EulerSR::S3_lr_group","EulerSR::tau_flux_group","EulerSR::tau_lr_group","EulerSR::v1_lr_group"};
-  GenericFD_AssertGroupStorage(cctkGH, "eulersr_cons_calc_flux_1", 13, groups);
-  
-  GenericFD_EnsureStencilFits(cctkGH, "eulersr_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,96 +86,102 @@ static void eulersr_cons_calc_flux_1_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 (eulersr_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 epsiLeftL = epsiLeft[index];
-        CCTK_REAL epsiRightL = epsiRight[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 tauLeftL = tauLeft[index];
-        CCTK_REAL tauRightL = tauRight[index];
-        CCTK_REAL v1LeftL = v1Left[index];
-        CCTK_REAL v1RightL = v1Right[index];
-        
-        
-        /* Include user supplied include files */
-        
-        /* Precompute derivatives */
-        CCTK_REAL const ShiftMinus1DenRight = ShiftMinus1(&DenRight[index]);
-        CCTK_REAL const ShiftMinus1epsiRight = ShiftMinus1(&epsiRight[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 ShiftMinus1tauRight = ShiftMinus1(&tauRight[index]);
-        CCTK_REAL const ShiftMinus1v1Right = ShiftMinus1(&v1Right[index]);
-        
-        /* Calculate temporaries and grid functions */
-        CCTK_REAL DenFluxLeft = DenLeftL*v1LeftL;
-        
-        CCTK_REAL DenFluxRight = ShiftMinus1DenRight*ShiftMinus1v1Right;
-        
-        CCTK_REAL DenFluxL = 0.5*(DenFluxLeft + DenFluxRight + (-DenLeftL + 
-          ShiftMinus1DenRight)*ToReal(hlleAlpha));
-        
-        CCTK_REAL S1FluxLeft = S1LeftL*v1LeftL + epsiLeftL*rhoLeftL*(-1 + 
-          ToReal(gamma));
-        
-        CCTK_REAL S1FluxRight = ShiftMinus1S1Right*ShiftMinus1v1Right + 
-          ShiftMinus1epsiRight*ShiftMinus1rhoRight*(-1 + ToReal(gamma));
-        
-        CCTK_REAL S1FluxL = 0.5*(S1FluxLeft + S1FluxRight + (-S1LeftL + 
-          ShiftMinus1S1Right)*ToReal(hlleAlpha));
-        
-        CCTK_REAL S2FluxLeft = S2LeftL*v1LeftL;
-        
-        CCTK_REAL S2FluxRight = ShiftMinus1S2Right*ShiftMinus1v1Right;
-        
-        CCTK_REAL S2FluxL = 0.5*(S2FluxLeft + S2FluxRight + (-S2LeftL + 
-          ShiftMinus1S2Right)*ToReal(hlleAlpha));
-        
-        CCTK_REAL S3FluxLeft = S3LeftL*v1LeftL;
-        
-        CCTK_REAL S3FluxRight = ShiftMinus1S3Right*ShiftMinus1v1Right;
-        
-        CCTK_REAL S3FluxL = 0.5*(S3FluxLeft + S3FluxRight + (-S3LeftL + 
-          ShiftMinus1S3Right)*ToReal(hlleAlpha));
-        
-        CCTK_REAL tauFluxLeft = v1LeftL*(tauLeftL + epsiLeftL*rhoLeftL*(-1 + 
-          ToReal(gamma)));
-        
-        CCTK_REAL tauFluxRight = ShiftMinus1v1Right*(ShiftMinus1tauRight + 
-          ShiftMinus1epsiRight*ShiftMinus1rhoRight*(-1 + ToReal(gamma)));
-        
-        CCTK_REAL tauFluxL = 0.5*(tauFluxLeft + tauFluxRight + 
-          (ShiftMinus1tauRight - tauLeftL)*ToReal(hlleAlpha));
-        
-        /* Copy local copies back to grid functions */
-        DenFlux[index] = DenFluxL;
-        S1Flux[index] = S1FluxL;
-        S2Flux[index] = S2FluxL;
-        S3Flux[index] = S3FluxL;
-        tauFlux[index] = tauFluxL;
-      }
-    }
+    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 epsiLeftL = epsiLeft[index];
+    CCTK_REAL epsiRightL = epsiRight[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 tauLeftL = tauLeft[index];
+    CCTK_REAL tauRightL = tauRight[index];
+    CCTK_REAL v1LeftL = v1Left[index];
+    CCTK_REAL v1RightL = v1Right[index];
+    
+    
+    /* Include user supplied include files */
+    
+    /* Precompute derivatives */
+    CCTK_REAL const ShiftMinus1DenRight = ShiftMinus1(&DenRight[index]);
+    CCTK_REAL const ShiftMinus1epsiRight = ShiftMinus1(&epsiRight[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 ShiftMinus1tauRight = ShiftMinus1(&tauRight[index]);
+    CCTK_REAL const ShiftMinus1v1Right = ShiftMinus1(&v1Right[index]);
+    
+    /* Calculate temporaries and grid functions */
+    CCTK_REAL DenFluxLeft = DenLeftL*v1LeftL;
+    
+    CCTK_REAL DenFluxRight = ShiftMinus1DenRight*ShiftMinus1v1Right;
+    
+    CCTK_REAL DenFluxL = 0.5*(DenFluxLeft + DenFluxRight + (-DenLeftL + 
+      ShiftMinus1DenRight)*ToReal(hlleAlpha));
+    
+    CCTK_REAL S1FluxLeft = S1LeftL*v1LeftL + 
+      epsiLeftL*rhoLeftL*(-1 + ToReal(gamma));
+    
+    CCTK_REAL S1FluxRight = ShiftMinus1S1Right*ShiftMinus1v1Right + 
+      ShiftMinus1epsiRight*ShiftMinus1rhoRight*(-1 + ToReal(gamma));
+    
+    CCTK_REAL S1FluxL = 0.5*(S1FluxLeft + S1FluxRight + (-S1LeftL + 
+      ShiftMinus1S1Right)*ToReal(hlleAlpha));
+    
+    CCTK_REAL S2FluxLeft = S2LeftL*v1LeftL;
+    
+    CCTK_REAL S2FluxRight = ShiftMinus1S2Right*ShiftMinus1v1Right;
+    
+    CCTK_REAL S2FluxL = 0.5*(S2FluxLeft + S2FluxRight + (-S2LeftL + 
+      ShiftMinus1S2Right)*ToReal(hlleAlpha));
+    
+    CCTK_REAL S3FluxLeft = S3LeftL*v1LeftL;
+    
+    CCTK_REAL S3FluxRight = ShiftMinus1S3Right*ShiftMinus1v1Right;
+    
+    CCTK_REAL S3FluxL = 0.5*(S3FluxLeft + S3FluxRight + (-S3LeftL + 
+      ShiftMinus1S3Right)*ToReal(hlleAlpha));
+    
+    CCTK_REAL tauFluxLeft = v1LeftL*(tauLeftL + 
+      epsiLeftL*rhoLeftL*(-1 + ToReal(gamma)));
+    
+    CCTK_REAL tauFluxRight = ShiftMinus1v1Right*(ShiftMinus1tauRight + 
+      ShiftMinus1epsiRight*ShiftMinus1rhoRight*(-1 + ToReal(gamma)));
+    
+    CCTK_REAL tauFluxL = 0.5*(tauFluxLeft + tauFluxRight + (-tauLeftL + 
+      ShiftMinus1tauRight)*ToReal(hlleAlpha));
+    
+    /* Copy local copies back to grid functions */
+    DenFlux[index] = DenFluxL;
+    S1Flux[index] = S1FluxL;
+    S2Flux[index] = S2FluxL;
+    S3Flux[index] = S3FluxL;
+    tauFlux[index] = tauFluxL;
   }
+  CCTK_ENDLOOP3 (eulersr_cons_calc_flux_1);
 }
 
 extern "C" void eulersr_cons_calc_flux_1(CCTK_ARGUMENTS)
@@ -196,5 +189,26 @@ extern "C" void eulersr_cons_calc_flux_1(CCTK_ARGUMENTS)
   DECLARE_CCTK_ARGUMENTS;
   DECLARE_CCTK_PARAMETERS;
   
+  
+  if (verbose > 1)
+  {
+    CCTK_VInfo(CCTK_THORNSTRING,"Entering eulersr_cons_calc_flux_1_Body");
+  }
+  
+  if (cctk_iteration % eulersr_cons_calc_flux_1_calc_every != eulersr_cons_calc_flux_1_calc_offset)
+  {
+    return;
+  }
+  
+  const char *groups[] = {"EulerSR::Den_flux_group","EulerSR::Den_lr_group","EulerSR::epsi_lr_group","EulerSR::rho_lr_group","EulerSR::S1_flux_group","EulerSR::S1_lr_group","EulerSR::S2_flux_group","EulerSR::S2_lr_group","EulerSR::S3_flux_group","EulerSR::S3_lr_group","EulerSR::tau_flux_group","EulerSR::tau_lr_group","EulerSR::v1_lr_group"};
+  GenericFD_AssertGroupStorage(cctkGH, "eulersr_cons_calc_flux_1", 13, groups);
+  
+  GenericFD_EnsureStencilFits(cctkGH, "eulersr_cons_calc_flux_1", 1, 1, 1);
+  
   GenericFD_LoopOverInterior(cctkGH, &eulersr_cons_calc_flux_1_Body);
+  
+  if (verbose > 1)
+  {
+    CCTK_VInfo(CCTK_THORNSTRING,"Leaving eulersr_cons_calc_flux_1_Body");
+  }
 }
diff --git a/Examples/EulerSR/src/eulersr_cons_calc_flux_2.cc b/Examples/EulerSR/src/eulersr_cons_calc_flux_2.cc
index b650410..fb8553c 100644
--- a/Examples/EulerSR/src/eulersr_cons_calc_flux_2.cc
+++ b/Examples/EulerSR/src/eulersr_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 eulersr_cons_calc_flux_2_SelectBCs(CCTK_ARGUMENTS)
   return;
 }
 
-static void eulersr_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 eulersr_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 eulersr_cons_calc_flux_2_Body");
-  }
-  
-  if (cctk_iteration % eulersr_cons_calc_flux_2_calc_every != eulersr_cons_calc_flux_2_calc_offset)
-  {
-    return;
-  }
-  
-  const char *groups[] = {"EulerSR::Den_flux_group","EulerSR::Den_lr_group","EulerSR::epsi_lr_group","EulerSR::rho_lr_group","EulerSR::S1_flux_group","EulerSR::S1_lr_group","EulerSR::S2_flux_group","EulerSR::S2_lr_group","EulerSR::S3_flux_group","EulerSR::S3_lr_group","EulerSR::tau_flux_group","EulerSR::tau_lr_group","EulerSR::v2_lr_group"};
-  GenericFD_AssertGroupStorage(cctkGH, "eulersr_cons_calc_flux_2", 13, groups);
-  
-  GenericFD_EnsureStencilFits(cctkGH, "eulersr_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,96 +86,102 @@ static void eulersr_cons_calc_flux_2_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 (eulersr_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 epsiLeftL = epsiLeft[index];
-        CCTK_REAL epsiRightL = epsiRight[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 tauLeftL = tauLeft[index];
-        CCTK_REAL tauRightL = tauRight[index];
-        CCTK_REAL v2LeftL = v2Left[index];
-        CCTK_REAL v2RightL = v2Right[index];
-        
-        
-        /* Include user supplied include files */
-        
-        /* Precompute derivatives */
-        CCTK_REAL const ShiftMinus2DenRight = ShiftMinus2(&DenRight[index]);
-        CCTK_REAL const ShiftMinus2epsiRight = ShiftMinus2(&epsiRight[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 ShiftMinus2tauRight = ShiftMinus2(&tauRight[index]);
-        CCTK_REAL const ShiftMinus2v2Right = ShiftMinus2(&v2Right[index]);
-        
-        /* Calculate temporaries and grid functions */
-        CCTK_REAL DenFluxLeft = DenLeftL*v2LeftL;
-        
-        CCTK_REAL DenFluxRight = ShiftMinus2DenRight*ShiftMinus2v2Right;
-        
-        CCTK_REAL DenFluxL = 0.5*(DenFluxLeft + DenFluxRight + (-DenLeftL + 
-          ShiftMinus2DenRight)*ToReal(hlleAlpha));
-        
-        CCTK_REAL S1FluxLeft = S1LeftL*v2LeftL;
-        
-        CCTK_REAL S1FluxRight = ShiftMinus2S1Right*ShiftMinus2v2Right;
-        
-        CCTK_REAL S1FluxL = 0.5*(S1FluxLeft + S1FluxRight + (-S1LeftL + 
-          ShiftMinus2S1Right)*ToReal(hlleAlpha));
-        
-        CCTK_REAL S2FluxLeft = S2LeftL*v2LeftL + epsiLeftL*rhoLeftL*(-1 + 
-          ToReal(gamma));
-        
-        CCTK_REAL S2FluxRight = ShiftMinus2S2Right*ShiftMinus2v2Right + 
-          ShiftMinus2epsiRight*ShiftMinus2rhoRight*(-1 + ToReal(gamma));
-        
-        CCTK_REAL S2FluxL = 0.5*(S2FluxLeft + S2FluxRight + (-S2LeftL + 
-          ShiftMinus2S2Right)*ToReal(hlleAlpha));
-        
-        CCTK_REAL S3FluxLeft = S3LeftL*v2LeftL;
-        
-        CCTK_REAL S3FluxRight = ShiftMinus2S3Right*ShiftMinus2v2Right;
-        
-        CCTK_REAL S3FluxL = 0.5*(S3FluxLeft + S3FluxRight + (-S3LeftL + 
-          ShiftMinus2S3Right)*ToReal(hlleAlpha));
-        
-        CCTK_REAL tauFluxLeft = v2LeftL*(tauLeftL + epsiLeftL*rhoLeftL*(-1 + 
-          ToReal(gamma)));
-        
-        CCTK_REAL tauFluxRight = ShiftMinus2v2Right*(ShiftMinus2tauRight + 
-          ShiftMinus2epsiRight*ShiftMinus2rhoRight*(-1 + ToReal(gamma)));
-        
-        CCTK_REAL tauFluxL = 0.5*(tauFluxLeft + tauFluxRight + 
-          (ShiftMinus2tauRight - tauLeftL)*ToReal(hlleAlpha));
-        
-        /* Copy local copies back to grid functions */
-        DenFlux[index] = DenFluxL;
-        S1Flux[index] = S1FluxL;
-        S2Flux[index] = S2FluxL;
-        S3Flux[index] = S3FluxL;
-        tauFlux[index] = tauFluxL;
-      }
-    }
+    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 epsiLeftL = epsiLeft[index];
+    CCTK_REAL epsiRightL = epsiRight[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 tauLeftL = tauLeft[index];
+    CCTK_REAL tauRightL = tauRight[index];
+    CCTK_REAL v2LeftL = v2Left[index];
+    CCTK_REAL v2RightL = v2Right[index];
+    
+    
+    /* Include user supplied include files */
+    
+    /* Precompute derivatives */
+    CCTK_REAL const ShiftMinus2DenRight = ShiftMinus2(&DenRight[index]);
+    CCTK_REAL const ShiftMinus2epsiRight = ShiftMinus2(&epsiRight[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 ShiftMinus2tauRight = ShiftMinus2(&tauRight[index]);
+    CCTK_REAL const ShiftMinus2v2Right = ShiftMinus2(&v2Right[index]);
+    
+    /* Calculate temporaries and grid functions */
+    CCTK_REAL DenFluxLeft = DenLeftL*v2LeftL;
+    
+    CCTK_REAL DenFluxRight = ShiftMinus2DenRight*ShiftMinus2v2Right;
+    
+    CCTK_REAL DenFluxL = 0.5*(DenFluxLeft + DenFluxRight + (-DenLeftL + 
+      ShiftMinus2DenRight)*ToReal(hlleAlpha));
+    
+    CCTK_REAL S1FluxLeft = S1LeftL*v2LeftL;
+    
+    CCTK_REAL S1FluxRight = ShiftMinus2S1Right*ShiftMinus2v2Right;
+    
+    CCTK_REAL S1FluxL = 0.5*(S1FluxLeft + S1FluxRight + (-S1LeftL + 
+      ShiftMinus2S1Right)*ToReal(hlleAlpha));
+    
+    CCTK_REAL S2FluxLeft = S2LeftL*v2LeftL + 
+      epsiLeftL*rhoLeftL*(-1 + ToReal(gamma));
+    
+    CCTK_REAL S2FluxRight = ShiftMinus2S2Right*ShiftMinus2v2Right + 
+      ShiftMinus2epsiRight*ShiftMinus2rhoRight*(-1 + ToReal(gamma));
+    
+    CCTK_REAL S2FluxL = 0.5*(S2FluxLeft + S2FluxRight + (-S2LeftL + 
+      ShiftMinus2S2Right)*ToReal(hlleAlpha));
+    
+    CCTK_REAL S3FluxLeft = S3LeftL*v2LeftL;
+    
+    CCTK_REAL S3FluxRight = ShiftMinus2S3Right*ShiftMinus2v2Right;
+    
+    CCTK_REAL S3FluxL = 0.5*(S3FluxLeft + S3FluxRight + (-S3LeftL + 
+      ShiftMinus2S3Right)*ToReal(hlleAlpha));
+    
+    CCTK_REAL tauFluxLeft = v2LeftL*(tauLeftL + 
+      epsiLeftL*rhoLeftL*(-1 + ToReal(gamma)));
+    
+    CCTK_REAL tauFluxRight = ShiftMinus2v2Right*(ShiftMinus2tauRight + 
+      ShiftMinus2epsiRight*ShiftMinus2rhoRight*(-1 + ToReal(gamma)));
+    
+    CCTK_REAL tauFluxL = 0.5*(tauFluxLeft + tauFluxRight + (-tauLeftL + 
+      ShiftMinus2tauRight)*ToReal(hlleAlpha));
+    
+    /* Copy local copies back to grid functions */
+    DenFlux[index] = DenFluxL;
+    S1Flux[index] = S1FluxL;
+    S2Flux[index] = S2FluxL;
+    S3Flux[index] = S3FluxL;
+    tauFlux[index] = tauFluxL;
   }
+  CCTK_ENDLOOP3 (eulersr_cons_calc_flux_2);
 }
 
 extern "C" void eulersr_cons_calc_flux_2(CCTK_ARGUMENTS)
@@ -196,5 +189,26 @@ extern "C" void eulersr_cons_calc_flux_2(CCTK_ARGUMENTS)
   DECLARE_CCTK_ARGUMENTS;
   DECLARE_CCTK_PARAMETERS;
   
+  
+  if (verbose > 1)
+  {
+    CCTK_VInfo(CCTK_THORNSTRING,"Entering eulersr_cons_calc_flux_2_Body");
+  }
+  
+  if (cctk_iteration % eulersr_cons_calc_flux_2_calc_every != eulersr_cons_calc_flux_2_calc_offset)
+  {
+    return;
+  }
+  
+  const char *groups[] = {"EulerSR::Den_flux_group","EulerSR::Den_lr_group","EulerSR::epsi_lr_group","EulerSR::rho_lr_group","EulerSR::S1_flux_group","EulerSR::S1_lr_group","EulerSR::S2_flux_group","EulerSR::S2_lr_group","EulerSR::S3_flux_group","EulerSR::S3_lr_group","EulerSR::tau_flux_group","EulerSR::tau_lr_group","EulerSR::v2_lr_group"};
+  GenericFD_AssertGroupStorage(cctkGH, "eulersr_cons_calc_flux_2", 13, groups);
+  
+  GenericFD_EnsureStencilFits(cctkGH, "eulersr_cons_calc_flux_2", 1, 1, 1);
+  
   GenericFD_LoopOverInterior(cctkGH, &eulersr_cons_calc_flux_2_Body);
+  
+  if (verbose > 1)
+  {
+    CCTK_VInfo(CCTK_THORNSTRING,"Leaving eulersr_cons_calc_flux_2_Body");
+  }
 }
diff --git a/Examples/EulerSR/src/eulersr_cons_calc_flux_3.cc b/Examples/EulerSR/src/eulersr_cons_calc_flux_3.cc
index 7f4dbea..90fc2b7 100644
--- a/Examples/EulerSR/src/eulersr_cons_calc_flux_3.cc
+++ b/Examples/EulerSR/src/eulersr_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 eulersr_cons_calc_flux_3_SelectBCs(CCTK_ARGUMENTS)
   return;
 }
 
-static void eulersr_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 eulersr_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 eulersr_cons_calc_flux_3_Body");
-  }
-  
-  if (cctk_iteration % eulersr_cons_calc_flux_3_calc_every != eulersr_cons_calc_flux_3_calc_offset)
-  {
-    return;
-  }
-  
-  const char *groups[] = {"EulerSR::Den_flux_group","EulerSR::Den_lr_group","EulerSR::epsi_lr_group","EulerSR::rho_lr_group","EulerSR::S1_flux_group","EulerSR::S1_lr_group","EulerSR::S2_flux_group","EulerSR::S2_lr_group","EulerSR::S3_flux_group","EulerSR::S3_lr_group","EulerSR::tau_flux_group","EulerSR::tau_lr_group","EulerSR::v3_lr_group"};
-  GenericFD_AssertGroupStorage(cctkGH, "eulersr_cons_calc_flux_3", 13, groups);
-  
-  GenericFD_EnsureStencilFits(cctkGH, "eulersr_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,96 +86,102 @@ static void eulersr_cons_calc_flux_3_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 (eulersr_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 epsiLeftL = epsiLeft[index];
-        CCTK_REAL epsiRightL = epsiRight[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 tauLeftL = tauLeft[index];
-        CCTK_REAL tauRightL = tauRight[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 ShiftMinus3epsiRight = ShiftMinus3(&epsiRight[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 ShiftMinus3tauRight = ShiftMinus3(&tauRight[index]);
-        CCTK_REAL const ShiftMinus3v3Right = ShiftMinus3(&v3Right[index]);
-        
-        /* Calculate temporaries and grid functions */
-        CCTK_REAL DenFluxLeft = DenLeftL*v3LeftL;
-        
-        CCTK_REAL DenFluxRight = ShiftMinus3DenRight*ShiftMinus3v3Right;
-        
-        CCTK_REAL DenFluxL = 0.5*(DenFluxLeft + DenFluxRight + (-DenLeftL + 
-          ShiftMinus3DenRight)*ToReal(hlleAlpha));
-        
-        CCTK_REAL S1FluxLeft = S1LeftL*v3LeftL;
-        
-        CCTK_REAL S1FluxRight = ShiftMinus3S1Right*ShiftMinus3v3Right;
-        
-        CCTK_REAL S1FluxL = 0.5*(S1FluxLeft + S1FluxRight + (-S1LeftL + 
-          ShiftMinus3S1Right)*ToReal(hlleAlpha));
-        
-        CCTK_REAL S2FluxLeft = S2LeftL*v3LeftL;
-        
-        CCTK_REAL S2FluxRight = ShiftMinus3S2Right*ShiftMinus3v3Right;
-        
-        CCTK_REAL S2FluxL = 0.5*(S2FluxLeft + S2FluxRight + (-S2LeftL + 
-          ShiftMinus3S2Right)*ToReal(hlleAlpha));
-        
-        CCTK_REAL S3FluxLeft = S3LeftL*v3LeftL + epsiLeftL*rhoLeftL*(-1 + 
-          ToReal(gamma));
-        
-        CCTK_REAL S3FluxRight = ShiftMinus3S3Right*ShiftMinus3v3Right + 
-          ShiftMinus3epsiRight*ShiftMinus3rhoRight*(-1 + ToReal(gamma));
-        
-        CCTK_REAL S3FluxL = 0.5*(S3FluxLeft + S3FluxRight + (-S3LeftL + 
-          ShiftMinus3S3Right)*ToReal(hlleAlpha));
-        
-        CCTK_REAL tauFluxLeft = v3LeftL*(tauLeftL + epsiLeftL*rhoLeftL*(-1 + 
-          ToReal(gamma)));
-        
-        CCTK_REAL tauFluxRight = ShiftMinus3v3Right*(ShiftMinus3tauRight + 
-          ShiftMinus3epsiRight*ShiftMinus3rhoRight*(-1 + ToReal(gamma)));
-        
-        CCTK_REAL tauFluxL = 0.5*(tauFluxLeft + tauFluxRight + 
-          (ShiftMinus3tauRight - tauLeftL)*ToReal(hlleAlpha));
-        
-        /* Copy local copies back to grid functions */
-        DenFlux[index] = DenFluxL;
-        S1Flux[index] = S1FluxL;
-        S2Flux[index] = S2FluxL;
-        S3Flux[index] = S3FluxL;
-        tauFlux[index] = tauFluxL;
-      }
-    }
+    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 epsiLeftL = epsiLeft[index];
+    CCTK_REAL epsiRightL = epsiRight[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 tauLeftL = tauLeft[index];
+    CCTK_REAL tauRightL = tauRight[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 ShiftMinus3epsiRight = ShiftMinus3(&epsiRight[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 ShiftMinus3tauRight = ShiftMinus3(&tauRight[index]);
+    CCTK_REAL const ShiftMinus3v3Right = ShiftMinus3(&v3Right[index]);
+    
+    /* Calculate temporaries and grid functions */
+    CCTK_REAL DenFluxLeft = DenLeftL*v3LeftL;
+    
+    CCTK_REAL DenFluxRight = ShiftMinus3DenRight*ShiftMinus3v3Right;
+    
+    CCTK_REAL DenFluxL = 0.5*(DenFluxLeft + DenFluxRight + (-DenLeftL + 
+      ShiftMinus3DenRight)*ToReal(hlleAlpha));
+    
+    CCTK_REAL S1FluxLeft = S1LeftL*v3LeftL;
+    
+    CCTK_REAL S1FluxRight = ShiftMinus3S1Right*ShiftMinus3v3Right;
+    
+    CCTK_REAL S1FluxL = 0.5*(S1FluxLeft + S1FluxRight + (-S1LeftL + 
+      ShiftMinus3S1Right)*ToReal(hlleAlpha));
+    
+    CCTK_REAL S2FluxLeft = S2LeftL*v3LeftL;
+    
+    CCTK_REAL S2FluxRight = ShiftMinus3S2Right*ShiftMinus3v3Right;
+    
+    CCTK_REAL S2FluxL = 0.5*(S2FluxLeft + S2FluxRight + (-S2LeftL + 
+      ShiftMinus3S2Right)*ToReal(hlleAlpha));
+    
+    CCTK_REAL S3FluxLeft = S3LeftL*v3LeftL + 
+      epsiLeftL*rhoLeftL*(-1 + ToReal(gamma));
+    
+    CCTK_REAL S3FluxRight = ShiftMinus3S3Right*ShiftMinus3v3Right + 
+      ShiftMinus3epsiRight*ShiftMinus3rhoRight*(-1 + ToReal(gamma));
+    
+    CCTK_REAL S3FluxL = 0.5*(S3FluxLeft + S3FluxRight + (-S3LeftL + 
+      ShiftMinus3S3Right)*ToReal(hlleAlpha));
+    
+    CCTK_REAL tauFluxLeft = v3LeftL*(tauLeftL + 
+      epsiLeftL*rhoLeftL*(-1 + ToReal(gamma)));
+    
+    CCTK_REAL tauFluxRight = ShiftMinus3v3Right*(ShiftMinus3tauRight + 
+      ShiftMinus3epsiRight*ShiftMinus3rhoRight*(-1 + ToReal(gamma)));
+    
+    CCTK_REAL tauFluxL = 0.5*(tauFluxLeft + tauFluxRight + (-tauLeftL + 
+      ShiftMinus3tauRight)*ToReal(hlleAlpha));
+    
+    /* Copy local copies back to grid functions */
+    DenFlux[index] = DenFluxL;
+    S1Flux[index] = S1FluxL;
+    S2Flux[index] = S2FluxL;
+    S3Flux[index] = S3FluxL;
+    tauFlux[index] = tauFluxL;
   }
+  CCTK_ENDLOOP3 (eulersr_cons_calc_flux_3);
 }
 
 extern "C" void eulersr_cons_calc_flux_3(CCTK_ARGUMENTS)
@@ -196,5 +189,26 @@ extern "C" void eulersr_cons_calc_flux_3(CCTK_ARGUMENTS)
   DECLARE_CCTK_ARGUMENTS;
   DECLARE_CCTK_PARAMETERS;
   
+  
+  if (verbose > 1)
+  {
+    CCTK_VInfo(CCTK_THORNSTRING,"Entering eulersr_cons_calc_flux_3_Body");
+  }
+  
+  if (cctk_iteration % eulersr_cons_calc_flux_3_calc_every != eulersr_cons_calc_flux_3_calc_offset)
+  {
+    return;
+  }
+  
+  const char *groups[] = {"EulerSR::Den_flux_group","EulerSR::Den_lr_group","EulerSR::epsi_lr_group","EulerSR::rho_lr_group","EulerSR::S1_flux_group","EulerSR::S1_lr_group","EulerSR::S2_flux_group","EulerSR::S2_lr_group","EulerSR::S3_flux_group","EulerSR::S3_lr_group","EulerSR::tau_flux_group","EulerSR::tau_lr_group","EulerSR::v3_lr_group"};
+  GenericFD_AssertGroupStorage(cctkGH, "eulersr_cons_calc_flux_3", 13, groups);
+  
+  GenericFD_EnsureStencilFits(cctkGH, "eulersr_cons_calc_flux_3", 1, 1, 1);
+  
   GenericFD_LoopOverInterior(cctkGH, &eulersr_cons_calc_flux_3_Body);
+  
+  if (verbose > 1)
+  {
+    CCTK_VInfo(CCTK_THORNSTRING,"Leaving eulersr_cons_calc_flux_3_Body");
+  }
 }
diff --git a/Examples/EulerSR/src/eulersr_cons_calc_intercell_conserved_1.cc b/Examples/EulerSR/src/eulersr_cons_calc_intercell_conserved_1.cc
index 6885845..df163c1 100644
--- a/Examples/EulerSR/src/eulersr_cons_calc_intercell_conserved_1.cc
+++ b/Examples/EulerSR/src/eulersr_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 eulersr_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 eulersr_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 eulersr_cons_calc_intercell_conserved_1_Body");
-  }
-  
-  if (cctk_iteration % eulersr_cons_calc_intercell_conserved_1_calc_every != eulersr_cons_calc_intercell_conserved_1_calc_offset)
-  {
-    return;
-  }
-  
-  const char *groups[] = {"EulerSR::Den_lr_group","EulerSR::epsi_lr_group","EulerSR::h_group","EulerSR::p_group","EulerSR::rho_lr_group","EulerSR::S1_lr_group","EulerSR::S2_lr_group","EulerSR::S3_lr_group","EulerSR::tau_lr_group","EulerSR::v1_lr_group","EulerSR::v2_lr_group","EulerSR::v3_lr_group","EulerSR::W_group"};
-  GenericFD_AssertGroupStorage(cctkGH, "eulersr_cons_calc_intercell_conserved_1", 13, 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,92 +62,99 @@ static void eulersr_cons_calc_intercell_conserved_1_Body(cGH const * restrict co
   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 (eulersr_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 DenLeftL = DenLeft[index];
-        CCTK_REAL DenRightL = DenRight[index];
-        CCTK_REAL epsiLeftL = epsiLeft[index];
-        CCTK_REAL epsiRightL = epsiRight[index];
-        CCTK_REAL hL = h[index];
-        CCTK_REAL pL = p[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];
-        CCTK_REAL WL = W[index];
-        
-        
-        /* Include user supplied include files */
-        
-        /* Precompute derivatives */
-        
-        /* Calculate temporaries and grid functions */
-        CCTK_REAL Wx = 1 - SQR(v1LeftL) - SQR(v2LeftL) - SQR(v3LeftL);
-        
-        WL = INV(sqrt(Wx));
-        
-        pL = epsiLeftL*rhoLeftL*(-1 + ToReal(gamma));
-        
-        hL = 1 + epsiLeftL + pL*INV(rhoLeftL);
-        
-        DenLeftL = rhoLeftL*WL;
-        
-        CCTK_REAL S1LeftL = hL*rhoLeftL*v1LeftL*SQR(WL);
-        
-        CCTK_REAL S2LeftL = hL*rhoLeftL*v2LeftL*SQR(WL);
-        
-        CCTK_REAL S3LeftL = hL*rhoLeftL*v3LeftL*SQR(WL);
-        
-        CCTK_REAL tauLeftL = -DenLeftL - pL + hL*rhoLeftL*SQR(WL);
-        
-        Wx = 1 - SQR(v1RightL) - SQR(v2RightL) - SQR(v3RightL);
-        
-        WL = INV(sqrt(Wx));
-        
-        pL = epsiRightL*rhoRightL*(-1 + ToReal(gamma));
-        
-        hL = 1 + epsiRightL + pL*INV(rhoRightL);
-        
-        DenRightL = rhoRightL*WL;
-        
-        CCTK_REAL S1RightL = hL*rhoRightL*v1RightL*SQR(WL);
-        
-        CCTK_REAL S2RightL = hL*rhoRightL*v2RightL*SQR(WL);
-        
-        CCTK_REAL S3RightL = hL*rhoRightL*v3RightL*SQR(WL);
-        
-        CCTK_REAL tauRightL = -DenRightL - pL + hL*rhoRightL*SQR(WL);
-        
-        /* Copy local copies back to grid functions */
-        DenLeft[index] = DenLeftL;
-        DenRight[index] = DenRightL;
-        h[index] = hL;
-        p[index] = pL;
-        S1Left[index] = S1LeftL;
-        S1Right[index] = S1RightL;
-        S2Left[index] = S2LeftL;
-        S2Right[index] = S2RightL;
-        S3Left[index] = S3LeftL;
-        S3Right[index] = S3RightL;
-        tauLeft[index] = tauLeftL;
-        tauRight[index] = tauRightL;
-        W[index] = WL;
-      }
-    }
+    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 epsiLeftL = epsiLeft[index];
+    CCTK_REAL epsiRightL = epsiRight[index];
+    CCTK_REAL hL = h[index];
+    CCTK_REAL pL = p[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];
+    CCTK_REAL WL = W[index];
+    
+    
+    /* Include user supplied include files */
+    
+    /* Precompute derivatives */
+    
+    /* Calculate temporaries and grid functions */
+    CCTK_REAL Wx = 1 - SQR(v1LeftL) - SQR(v2LeftL) - SQR(v3LeftL);
+    
+    WL = INV(sqrt(Wx));
+    
+    pL = epsiLeftL*rhoLeftL*(-1 + ToReal(gamma));
+    
+    hL = 1 + epsiLeftL + pL*INV(rhoLeftL);
+    
+    DenLeftL = rhoLeftL*WL;
+    
+    CCTK_REAL S1LeftL = hL*rhoLeftL*v1LeftL*SQR(WL);
+    
+    CCTK_REAL S2LeftL = hL*rhoLeftL*v2LeftL*SQR(WL);
+    
+    CCTK_REAL S3LeftL = hL*rhoLeftL*v3LeftL*SQR(WL);
+    
+    CCTK_REAL tauLeftL = -DenLeftL - pL + hL*rhoLeftL*SQR(WL);
+    
+    Wx = 1 - SQR(v1RightL) - SQR(v2RightL) - SQR(v3RightL);
+    
+    WL = INV(sqrt(Wx));
+    
+    pL = epsiRightL*rhoRightL*(-1 + ToReal(gamma));
+    
+    hL = 1 + epsiRightL + pL*INV(rhoRightL);
+    
+    DenRightL = rhoRightL*WL;
+    
+    CCTK_REAL S1RightL = hL*rhoRightL*v1RightL*SQR(WL);
+    
+    CCTK_REAL S2RightL = hL*rhoRightL*v2RightL*SQR(WL);
+    
+    CCTK_REAL S3RightL = hL*rhoRightL*v3RightL*SQR(WL);
+    
+    CCTK_REAL tauRightL = -DenRightL - pL + 
+      hL*rhoRightL*SQR(WL);
+    
+    /* Copy local copies back to grid functions */
+    DenLeft[index] = DenLeftL;
+    DenRight[index] = DenRightL;
+    h[index] = hL;
+    p[index] = pL;
+    S1Left[index] = S1LeftL;
+    S1Right[index] = S1RightL;
+    S2Left[index] = S2LeftL;
+    S2Right[index] = S2RightL;
+    S3Left[index] = S3LeftL;
+    S3Right[index] = S3RightL;
+    tauLeft[index] = tauLeftL;
+    tauRight[index] = tauRightL;
+    W[index] = WL;
   }
+  CCTK_ENDLOOP3 (eulersr_cons_calc_intercell_conserved_1);
 }
 
 extern "C" void eulersr_cons_calc_intercell_conserved_1(CCTK_ARGUMENTS)
@@ -167,5 +162,25 @@ extern "C" void eulersr_cons_calc_intercell_conserved_1(CCTK_ARGUMENTS)
   DECLARE_CCTK_ARGUMENTS;
   DECLARE_CCTK_PARAMETERS;
   
+  
+  if (verbose > 1)
+  {
+    CCTK_VInfo(CCTK_THORNSTRING,"Entering eulersr_cons_calc_intercell_conserved_1_Body");
+  }
+  
+  if (cctk_iteration % eulersr_cons_calc_intercell_conserved_1_calc_every != eulersr_cons_calc_intercell_conserved_1_calc_offset)
+  {
+    return;
+  }
+  
+  const char *groups[] = {"EulerSR::Den_lr_group","EulerSR::epsi_lr_group","EulerSR::h_group","EulerSR::p_group","EulerSR::rho_lr_group","EulerSR::S1_lr_group","EulerSR::S2_lr_group","EulerSR::S3_lr_group","EulerSR::tau_lr_group","EulerSR::v1_lr_group","EulerSR::v2_lr_group","EulerSR::v3_lr_group","EulerSR::W_group"};
+  GenericFD_AssertGroupStorage(cctkGH, "eulersr_cons_calc_intercell_conserved_1", 13, groups);
+  
+  
   GenericFD_LoopOverEverything(cctkGH, &eulersr_cons_calc_intercell_conserved_1_Body);
+  
+  if (verbose > 1)
+  {
+    CCTK_VInfo(CCTK_THORNSTRING,"Leaving eulersr_cons_calc_intercell_conserved_1_Body");
+  }
 }
diff --git a/Examples/EulerSR/src/eulersr_cons_calc_intercell_conserved_2.cc b/Examples/EulerSR/src/eulersr_cons_calc_intercell_conserved_2.cc
index 2263c6a..83616ac 100644
--- a/Examples/EulerSR/src/eulersr_cons_calc_intercell_conserved_2.cc
+++ b/Examples/EulerSR/src/eulersr_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 eulersr_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 eulersr_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 eulersr_cons_calc_intercell_conserved_2_Body");
-  }
-  
-  if (cctk_iteration % eulersr_cons_calc_intercell_conserved_2_calc_every != eulersr_cons_calc_intercell_conserved_2_calc_offset)
-  {
-    return;
-  }
-  
-  const char *groups[] = {"EulerSR::Den_lr_group","EulerSR::epsi_lr_group","EulerSR::h_group","EulerSR::p_group","EulerSR::rho_lr_group","EulerSR::S1_lr_group","EulerSR::S2_lr_group","EulerSR::S3_lr_group","EulerSR::tau_lr_group","EulerSR::v1_lr_group","EulerSR::v2_lr_group","EulerSR::v3_lr_group","EulerSR::W_group"};
-  GenericFD_AssertGroupStorage(cctkGH, "eulersr_cons_calc_intercell_conserved_2", 13, 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,92 +62,99 @@ static void eulersr_cons_calc_intercell_conserved_2_Body(cGH const * restrict co
   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 (eulersr_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 DenLeftL = DenLeft[index];
-        CCTK_REAL DenRightL = DenRight[index];
-        CCTK_REAL epsiLeftL = epsiLeft[index];
-        CCTK_REAL epsiRightL = epsiRight[index];
-        CCTK_REAL hL = h[index];
-        CCTK_REAL pL = p[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];
-        CCTK_REAL WL = W[index];
-        
-        
-        /* Include user supplied include files */
-        
-        /* Precompute derivatives */
-        
-        /* Calculate temporaries and grid functions */
-        CCTK_REAL Wx = 1 - SQR(v1LeftL) - SQR(v2LeftL) - SQR(v3LeftL);
-        
-        WL = INV(sqrt(Wx));
-        
-        pL = epsiLeftL*rhoLeftL*(-1 + ToReal(gamma));
-        
-        hL = 1 + epsiLeftL + pL*INV(rhoLeftL);
-        
-        DenLeftL = rhoLeftL*WL;
-        
-        CCTK_REAL S1LeftL = hL*rhoLeftL*v1LeftL*SQR(WL);
-        
-        CCTK_REAL S2LeftL = hL*rhoLeftL*v2LeftL*SQR(WL);
-        
-        CCTK_REAL S3LeftL = hL*rhoLeftL*v3LeftL*SQR(WL);
-        
-        CCTK_REAL tauLeftL = -DenLeftL - pL + hL*rhoLeftL*SQR(WL);
-        
-        Wx = 1 - SQR(v1RightL) - SQR(v2RightL) - SQR(v3RightL);
-        
-        WL = INV(sqrt(Wx));
-        
-        pL = epsiRightL*rhoRightL*(-1 + ToReal(gamma));
-        
-        hL = 1 + epsiRightL + pL*INV(rhoRightL);
-        
-        DenRightL = rhoRightL*WL;
-        
-        CCTK_REAL S1RightL = hL*rhoRightL*v1RightL*SQR(WL);
-        
-        CCTK_REAL S2RightL = hL*rhoRightL*v2RightL*SQR(WL);
-        
-        CCTK_REAL S3RightL = hL*rhoRightL*v3RightL*SQR(WL);
-        
-        CCTK_REAL tauRightL = -DenRightL - pL + hL*rhoRightL*SQR(WL);
-        
-        /* Copy local copies back to grid functions */
-        DenLeft[index] = DenLeftL;
-        DenRight[index] = DenRightL;
-        h[index] = hL;
-        p[index] = pL;
-        S1Left[index] = S1LeftL;
-        S1Right[index] = S1RightL;
-        S2Left[index] = S2LeftL;
-        S2Right[index] = S2RightL;
-        S3Left[index] = S3LeftL;
-        S3Right[index] = S3RightL;
-        tauLeft[index] = tauLeftL;
-        tauRight[index] = tauRightL;
-        W[index] = WL;
-      }
-    }
+    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 epsiLeftL = epsiLeft[index];
+    CCTK_REAL epsiRightL = epsiRight[index];
+    CCTK_REAL hL = h[index];
+    CCTK_REAL pL = p[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];
+    CCTK_REAL WL = W[index];
+    
+    
+    /* Include user supplied include files */
+    
+    /* Precompute derivatives */
+    
+    /* Calculate temporaries and grid functions */
+    CCTK_REAL Wx = 1 - SQR(v1LeftL) - SQR(v2LeftL) - SQR(v3LeftL);
+    
+    WL = INV(sqrt(Wx));
+    
+    pL = epsiLeftL*rhoLeftL*(-1 + ToReal(gamma));
+    
+    hL = 1 + epsiLeftL + pL*INV(rhoLeftL);
+    
+    DenLeftL = rhoLeftL*WL;
+    
+    CCTK_REAL S1LeftL = hL*rhoLeftL*v1LeftL*SQR(WL);
+    
+    CCTK_REAL S2LeftL = hL*rhoLeftL*v2LeftL*SQR(WL);
+    
+    CCTK_REAL S3LeftL = hL*rhoLeftL*v3LeftL*SQR(WL);
+    
+    CCTK_REAL tauLeftL = -DenLeftL - pL + hL*rhoLeftL*SQR(WL);
+    
+    Wx = 1 - SQR(v1RightL) - SQR(v2RightL) - SQR(v3RightL);
+    
+    WL = INV(sqrt(Wx));
+    
+    pL = epsiRightL*rhoRightL*(-1 + ToReal(gamma));
+    
+    hL = 1 + epsiRightL + pL*INV(rhoRightL);
+    
+    DenRightL = rhoRightL*WL;
+    
+    CCTK_REAL S1RightL = hL*rhoRightL*v1RightL*SQR(WL);
+    
+    CCTK_REAL S2RightL = hL*rhoRightL*v2RightL*SQR(WL);
+    
+    CCTK_REAL S3RightL = hL*rhoRightL*v3RightL*SQR(WL);
+    
+    CCTK_REAL tauRightL = -DenRightL - pL + 
+      hL*rhoRightL*SQR(WL);
+    
+    /* Copy local copies back to grid functions */
+    DenLeft[index] = DenLeftL;
+    DenRight[index] = DenRightL;
+    h[index] = hL;
+    p[index] = pL;
+    S1Left[index] = S1LeftL;
+    S1Right[index] = S1RightL;
+    S2Left[index] = S2LeftL;
+    S2Right[index] = S2RightL;
+    S3Left[index] = S3LeftL;
+    S3Right[index] = S3RightL;
+    tauLeft[index] = tauLeftL;
+    tauRight[index] = tauRightL;
+    W[index] = WL;
   }
+  CCTK_ENDLOOP3 (eulersr_cons_calc_intercell_conserved_2);
 }
 
 extern "C" void eulersr_cons_calc_intercell_conserved_2(CCTK_ARGUMENTS)
@@ -167,5 +162,25 @@ extern "C" void eulersr_cons_calc_intercell_conserved_2(CCTK_ARGUMENTS)
   DECLARE_CCTK_ARGUMENTS;
   DECLARE_CCTK_PARAMETERS;
   
+  
+  if (verbose > 1)
+  {
+    CCTK_VInfo(CCTK_THORNSTRING,"Entering eulersr_cons_calc_intercell_conserved_2_Body");
+  }
+  
+  if (cctk_iteration % eulersr_cons_calc_intercell_conserved_2_calc_every != eulersr_cons_calc_intercell_conserved_2_calc_offset)
+  {
+    return;
+  }
+  
+  const char *groups[] = {"EulerSR::Den_lr_group","EulerSR::epsi_lr_group","EulerSR::h_group","EulerSR::p_group","EulerSR::rho_lr_group","EulerSR::S1_lr_group","EulerSR::S2_lr_group","EulerSR::S3_lr_group","EulerSR::tau_lr_group","EulerSR::v1_lr_group","EulerSR::v2_lr_group","EulerSR::v3_lr_group","EulerSR::W_group"};
+  GenericFD_AssertGroupStorage(cctkGH, "eulersr_cons_calc_intercell_conserved_2", 13, groups);
+  
+  
   GenericFD_LoopOverEverything(cctkGH, &eulersr_cons_calc_intercell_conserved_2_Body);
+  
+  if (verbose > 1)
+  {
+    CCTK_VInfo(CCTK_THORNSTRING,"Leaving eulersr_cons_calc_intercell_conserved_2_Body");
+  }
 }
diff --git a/Examples/EulerSR/src/eulersr_cons_calc_intercell_conserved_3.cc b/Examples/EulerSR/src/eulersr_cons_calc_intercell_conserved_3.cc
index f636073..4cc9ef2 100644
--- a/Examples/EulerSR/src/eulersr_cons_calc_intercell_conserved_3.cc
+++ b/Examples/EulerSR/src/eulersr_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 eulersr_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 eulersr_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 eulersr_cons_calc_intercell_conserved_3_Body");
-  }
-  
-  if (cctk_iteration % eulersr_cons_calc_intercell_conserved_3_calc_every != eulersr_cons_calc_intercell_conserved_3_calc_offset)
-  {
-    return;
-  }
-  
-  const char *groups[] = {"EulerSR::Den_lr_group","EulerSR::epsi_lr_group","EulerSR::h_group","EulerSR::p_group","EulerSR::rho_lr_group","EulerSR::S1_lr_group","EulerSR::S2_lr_group","EulerSR::S3_lr_group","EulerSR::tau_lr_group","EulerSR::v1_lr_group","EulerSR::v2_lr_group","EulerSR::v3_lr_group","EulerSR::W_group"};
-  GenericFD_AssertGroupStorage(cctkGH, "eulersr_cons_calc_intercell_conserved_3", 13, 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,92 +62,99 @@ static void eulersr_cons_calc_intercell_conserved_3_Body(cGH const * restrict co
   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 (eulersr_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 DenLeftL = DenLeft[index];
-        CCTK_REAL DenRightL = DenRight[index];
-        CCTK_REAL epsiLeftL = epsiLeft[index];
-        CCTK_REAL epsiRightL = epsiRight[index];
-        CCTK_REAL hL = h[index];
-        CCTK_REAL pL = p[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];
-        CCTK_REAL WL = W[index];
-        
-        
-        /* Include user supplied include files */
-        
-        /* Precompute derivatives */
-        
-        /* Calculate temporaries and grid functions */
-        CCTK_REAL Wx = 1 - SQR(v1LeftL) - SQR(v2LeftL) - SQR(v3LeftL);
-        
-        WL = INV(sqrt(Wx));
-        
-        pL = epsiLeftL*rhoLeftL*(-1 + ToReal(gamma));
-        
-        hL = 1 + epsiLeftL + pL*INV(rhoLeftL);
-        
-        DenLeftL = rhoLeftL*WL;
-        
-        CCTK_REAL S1LeftL = hL*rhoLeftL*v1LeftL*SQR(WL);
-        
-        CCTK_REAL S2LeftL = hL*rhoLeftL*v2LeftL*SQR(WL);
-        
-        CCTK_REAL S3LeftL = hL*rhoLeftL*v3LeftL*SQR(WL);
-        
-        CCTK_REAL tauLeftL = -DenLeftL - pL + hL*rhoLeftL*SQR(WL);
-        
-        Wx = 1 - SQR(v1RightL) - SQR(v2RightL) - SQR(v3RightL);
-        
-        WL = INV(sqrt(Wx));
-        
-        pL = epsiRightL*rhoRightL*(-1 + ToReal(gamma));
-        
-        hL = 1 + epsiRightL + pL*INV(rhoRightL);
-        
-        DenRightL = rhoRightL*WL;
-        
-        CCTK_REAL S1RightL = hL*rhoRightL*v1RightL*SQR(WL);
-        
-        CCTK_REAL S2RightL = hL*rhoRightL*v2RightL*SQR(WL);
-        
-        CCTK_REAL S3RightL = hL*rhoRightL*v3RightL*SQR(WL);
-        
-        CCTK_REAL tauRightL = -DenRightL - pL + hL*rhoRightL*SQR(WL);
-        
-        /* Copy local copies back to grid functions */
-        DenLeft[index] = DenLeftL;
-        DenRight[index] = DenRightL;
-        h[index] = hL;
-        p[index] = pL;
-        S1Left[index] = S1LeftL;
-        S1Right[index] = S1RightL;
-        S2Left[index] = S2LeftL;
-        S2Right[index] = S2RightL;
-        S3Left[index] = S3LeftL;
-        S3Right[index] = S3RightL;
-        tauLeft[index] = tauLeftL;
-        tauRight[index] = tauRightL;
-        W[index] = WL;
-      }
-    }
+    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 epsiLeftL = epsiLeft[index];
+    CCTK_REAL epsiRightL = epsiRight[index];
+    CCTK_REAL hL = h[index];
+    CCTK_REAL pL = p[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];
+    CCTK_REAL WL = W[index];
+    
+    
+    /* Include user supplied include files */
+    
+    /* Precompute derivatives */
+    
+    /* Calculate temporaries and grid functions */
+    CCTK_REAL Wx = 1 - SQR(v1LeftL) - SQR(v2LeftL) - SQR(v3LeftL);
+    
+    WL = INV(sqrt(Wx));
+    
+    pL = epsiLeftL*rhoLeftL*(-1 + ToReal(gamma));
+    
+    hL = 1 + epsiLeftL + pL*INV(rhoLeftL);
+    
+    DenLeftL = rhoLeftL*WL;
+    
+    CCTK_REAL S1LeftL = hL*rhoLeftL*v1LeftL*SQR(WL);
+    
+    CCTK_REAL S2LeftL = hL*rhoLeftL*v2LeftL*SQR(WL);
+    
+    CCTK_REAL S3LeftL = hL*rhoLeftL*v3LeftL*SQR(WL);
+    
+    CCTK_REAL tauLeftL = -DenLeftL - pL + hL*rhoLeftL*SQR(WL);
+    
+    Wx = 1 - SQR(v1RightL) - SQR(v2RightL) - SQR(v3RightL);
+    
+    WL = INV(sqrt(Wx));
+    
+    pL = epsiRightL*rhoRightL*(-1 + ToReal(gamma));
+    
+    hL = 1 + epsiRightL + pL*INV(rhoRightL);
+    
+    DenRightL = rhoRightL*WL;
+    
+    CCTK_REAL S1RightL = hL*rhoRightL*v1RightL*SQR(WL);
+    
+    CCTK_REAL S2RightL = hL*rhoRightL*v2RightL*SQR(WL);
+    
+    CCTK_REAL S3RightL = hL*rhoRightL*v3RightL*SQR(WL);
+    
+    CCTK_REAL tauRightL = -DenRightL - pL + 
+      hL*rhoRightL*SQR(WL);
+    
+    /* Copy local copies back to grid functions */
+    DenLeft[index] = DenLeftL;
+    DenRight[index] = DenRightL;
+    h[index] = hL;
+    p[index] = pL;
+    S1Left[index] = S1LeftL;
+    S1Right[index] = S1RightL;
+    S2Left[index] = S2LeftL;
+    S2Right[index] = S2RightL;
+    S3Left[index] = S3LeftL;
+    S3Right[index] = S3RightL;
+    tauLeft[index] = tauLeftL;
+    tauRight[index] = tauRightL;
+    W[index] = WL;
   }
+  CCTK_ENDLOOP3 (eulersr_cons_calc_intercell_conserved_3);
 }
 
 extern "C" void eulersr_cons_calc_intercell_conserved_3(CCTK_ARGUMENTS)
@@ -167,5 +162,25 @@ extern "C" void eulersr_cons_calc_intercell_conserved_3(CCTK_ARGUMENTS)
   DECLARE_CCTK_ARGUMENTS;
   DECLARE_CCTK_PARAMETERS;
   
+  
+  if (verbose > 1)
+  {
+    CCTK_VInfo(CCTK_THORNSTRING,"Entering eulersr_cons_calc_intercell_conserved_3_Body");
+  }
+  
+  if (cctk_iteration % eulersr_cons_calc_intercell_conserved_3_calc_every != eulersr_cons_calc_intercell_conserved_3_calc_offset)
+  {
+    return;
+  }
+  
+  const char *groups[] = {"EulerSR::Den_lr_group","EulerSR::epsi_lr_group","EulerSR::h_group","EulerSR::p_group","EulerSR::rho_lr_group","EulerSR::S1_lr_group","EulerSR::S2_lr_group","EulerSR::S3_lr_group","EulerSR::tau_lr_group","EulerSR::v1_lr_group","EulerSR::v2_lr_group","EulerSR::v3_lr_group","EulerSR::W_group"};
+  GenericFD_AssertGroupStorage(cctkGH, "eulersr_cons_calc_intercell_conserved_3", 13, groups);
+  
+  
   GenericFD_LoopOverEverything(cctkGH, &eulersr_cons_calc_intercell_conserved_3_Body);
+  
+  if (verbose > 1)
+  {
+    CCTK_VInfo(CCTK_THORNSTRING,"Leaving eulersr_cons_calc_intercell_conserved_3_Body");
+  }
 }
diff --git a/Examples/EulerSR/src/eulersr_cons_calc_primitives.cc b/Examples/EulerSR/src/eulersr_cons_calc_primitives.cc
index c1ca939..b43f41a 100644
--- a/Examples/EulerSR/src/eulersr_cons_calc_primitives.cc
+++ b/Examples/EulerSR/src/eulersr_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 eulersr_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 eulersr_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 eulersr_cons_calc_primitives_Body");
-  }
-  
-  if (cctk_iteration % eulersr_cons_calc_primitives_calc_every != eulersr_cons_calc_primitives_calc_offset)
-  {
-    return;
-  }
-  
-  const char *groups[] = {"EulerSR::Den_group","EulerSR::epsi_group","EulerSR::h_group","EulerSR::p_group","EulerSR::rho_group","EulerSR::S_group","EulerSR::tau_group","EulerSR::v_group","EulerSR::W_group"};
-  GenericFD_AssertGroupStorage(cctkGH, "eulersr_cons_calc_primitives", 9, 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,178 +62,185 @@ static void eulersr_cons_calc_primitives_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 (eulersr_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 epsiL = epsi[index];
-        CCTK_REAL hL = h[index];
-        CCTK_REAL pL = p[index];
-        CCTK_REAL rhoL = rho[index];
-        CCTK_REAL S1L = S1[index];
-        CCTK_REAL S2L = S2[index];
-        CCTK_REAL S3L = S3[index];
-        CCTK_REAL tauL = tau[index];
-        CCTK_REAL WL = W[index];
-        
-        
-        /* Include user supplied include files */
-        
-        /* Precompute derivatives */
-        
-        /* Calculate temporaries and grid functions */
-        CCTK_REAL pBar = pL;
-        
-        CCTK_REAL f = 10;
-        
-        CCTK_REAL Z = DenL + pBar + tauL;
-        
-        CCTK_REAL Ssq = SQR(S1L) + SQR(S2L) + SQR(S3L);
-        
-        CCTK_REAL vsq = Ssq*INV(SQR(Z));
-        
-        WL = INV(sqrt(1 - vsq));
-        
-        rhoL = DenL*INV(WL);
-        
-        hL = Z*INV(rhoL)*INV(SQR(WL));
-        
-        epsiL = hL - (pBar + rhoL)*INV(rhoL);
-        
-        CCTK_REAL pEOS = epsiL*rhoL*(-1 + ToReal(gamma));
-        
-        f = -pBar + pEOS;
-        
-        CCTK_REAL cs = sqrt(epsiL*INV(hL)*(-1 + ToReal(gamma))*ToReal(gamma));
-        
-        CCTK_REAL df = -1 + vsq*SQR(cs);
-        
-        pBar = pBar - f*INV(df);
-        
-        Z = DenL + pBar + tauL;
-        
-        Ssq = SQR(S1L) + SQR(S2L) + SQR(S3L);
-        
-        vsq = Ssq*INV(SQR(Z));
-        
-        WL = INV(sqrt(1 - vsq));
-        
-        rhoL = DenL*INV(WL);
-        
-        hL = Z*INV(rhoL)*INV(SQR(WL));
-        
-        epsiL = hL - (pBar + rhoL)*INV(rhoL);
-        
-        pEOS = epsiL*rhoL*(-1 + ToReal(gamma));
-        
-        f = -pBar + pEOS;
-        
-        cs = sqrt(epsiL*INV(hL)*(-1 + ToReal(gamma))*ToReal(gamma));
-        
-        df = -1 + vsq*SQR(cs);
-        
-        pBar = pBar - f*INV(df);
-        
-        Z = DenL + pBar + tauL;
-        
-        Ssq = SQR(S1L) + SQR(S2L) + SQR(S3L);
-        
-        vsq = Ssq*INV(SQR(Z));
-        
-        WL = INV(sqrt(1 - vsq));
-        
-        rhoL = DenL*INV(WL);
-        
-        hL = Z*INV(rhoL)*INV(SQR(WL));
-        
-        epsiL = hL - (pBar + rhoL)*INV(rhoL);
-        
-        pEOS = epsiL*rhoL*(-1 + ToReal(gamma));
-        
-        f = -pBar + pEOS;
-        
-        cs = sqrt(epsiL*INV(hL)*(-1 + ToReal(gamma))*ToReal(gamma));
-        
-        df = -1 + vsq*SQR(cs);
-        
-        pBar = pBar - f*INV(df);
-        
-        Z = DenL + pBar + tauL;
-        
-        Ssq = SQR(S1L) + SQR(S2L) + SQR(S3L);
-        
-        vsq = Ssq*INV(SQR(Z));
-        
-        WL = INV(sqrt(1 - vsq));
-        
-        rhoL = DenL*INV(WL);
-        
-        hL = Z*INV(rhoL)*INV(SQR(WL));
-        
-        epsiL = hL - (pBar + rhoL)*INV(rhoL);
-        
-        pEOS = epsiL*rhoL*(-1 + ToReal(gamma));
-        
-        f = -pBar + pEOS;
-        
-        cs = sqrt(epsiL*INV(hL)*(-1 + ToReal(gamma))*ToReal(gamma));
-        
-        df = -1 + vsq*SQR(cs);
-        
-        pBar = pBar - f*INV(df);
-        
-        Z = DenL + pBar + tauL;
-        
-        Ssq = SQR(S1L) + SQR(S2L) + SQR(S3L);
-        
-        vsq = Ssq*INV(SQR(Z));
-        
-        WL = INV(sqrt(1 - vsq));
-        
-        rhoL = DenL*INV(WL);
-        
-        hL = Z*INV(rhoL)*INV(SQR(WL));
-        
-        epsiL = hL - (pBar + rhoL)*INV(rhoL);
-        
-        pEOS = epsiL*rhoL*(-1 + ToReal(gamma));
-        
-        f = -pBar + pEOS;
-        
-        cs = sqrt(epsiL*INV(hL)*(-1 + ToReal(gamma))*ToReal(gamma));
-        
-        df = -1 + vsq*SQR(cs);
-        
-        pBar = pBar - f*INV(df);
-        
-        pL = pBar;
-        
-        CCTK_REAL v1L = S1L*INV(hL)*INV(rhoL)*INV(SQR(WL));
-        
-        CCTK_REAL v2L = S2L*INV(hL)*INV(rhoL)*INV(SQR(WL));
-        
-        CCTK_REAL v3L = S3L*INV(hL)*INV(rhoL)*INV(SQR(WL));
-        
-        /* Copy local copies back to grid functions */
-        epsi[index] = epsiL;
-        h[index] = hL;
-        p[index] = pL;
-        rho[index] = rhoL;
-        v1[index] = v1L;
-        v2[index] = v2L;
-        v3[index] = v3L;
-        W[index] = WL;
-      }
-    }
+    ptrdiff_t const index = di*i + dj*j + dk*k;
+    
+    /* Assign local copies of grid functions */
+    
+    CCTK_REAL DenL = Den[index];
+    CCTK_REAL epsiL = epsi[index];
+    CCTK_REAL hL = h[index];
+    CCTK_REAL pL = p[index];
+    CCTK_REAL rhoL = rho[index];
+    CCTK_REAL S1L = S1[index];
+    CCTK_REAL S2L = S2[index];
+    CCTK_REAL S3L = S3[index];
+    CCTK_REAL tauL = tau[index];
+    CCTK_REAL WL = W[index];
+    
+    
+    /* Include user supplied include files */
+    
+    /* Precompute derivatives */
+    
+    /* Calculate temporaries and grid functions */
+    CCTK_REAL pBar = pL;
+    
+    CCTK_REAL f = 10;
+    
+    CCTK_REAL Z = DenL + tauL + pBar;
+    
+    CCTK_REAL Ssq = SQR(S1L) + SQR(S2L) + SQR(S3L);
+    
+    CCTK_REAL vsq = Ssq*INV(SQR(Z));
+    
+    WL = INV(sqrt(1 - vsq));
+    
+    rhoL = DenL*INV(WL);
+    
+    hL = Z*INV(rhoL)*INV(SQR(WL));
+    
+    epsiL = hL - (rhoL + pBar)*INV(rhoL);
+    
+    CCTK_REAL pEOS = epsiL*rhoL*(-1 + ToReal(gamma));
+    
+    f = -pBar + pEOS;
+    
+    CCTK_REAL cs = sqrt(epsiL*INV(hL)*(-1 + 
+      ToReal(gamma))*ToReal(gamma));
+    
+    CCTK_REAL df = -1 + vsq*SQR(cs);
+    
+    pBar = pBar - f*INV(df);
+    
+    Z = DenL + tauL + pBar;
+    
+    Ssq = SQR(S1L) + SQR(S2L) + SQR(S3L);
+    
+    vsq = Ssq*INV(SQR(Z));
+    
+    WL = INV(sqrt(1 - vsq));
+    
+    rhoL = DenL*INV(WL);
+    
+    hL = Z*INV(rhoL)*INV(SQR(WL));
+    
+    epsiL = hL - (rhoL + pBar)*INV(rhoL);
+    
+    pEOS = epsiL*rhoL*(-1 + ToReal(gamma));
+    
+    f = -pBar + pEOS;
+    
+    cs = sqrt(epsiL*INV(hL)*(-1 + ToReal(gamma))*ToReal(gamma));
+    
+    df = -1 + vsq*SQR(cs);
+    
+    pBar = pBar - f*INV(df);
+    
+    Z = DenL + tauL + pBar;
+    
+    Ssq = SQR(S1L) + SQR(S2L) + SQR(S3L);
+    
+    vsq = Ssq*INV(SQR(Z));
+    
+    WL = INV(sqrt(1 - vsq));
+    
+    rhoL = DenL*INV(WL);
+    
+    hL = Z*INV(rhoL)*INV(SQR(WL));
+    
+    epsiL = hL - (rhoL + pBar)*INV(rhoL);
+    
+    pEOS = epsiL*rhoL*(-1 + ToReal(gamma));
+    
+    f = -pBar + pEOS;
+    
+    cs = sqrt(epsiL*INV(hL)*(-1 + ToReal(gamma))*ToReal(gamma));
+    
+    df = -1 + vsq*SQR(cs);
+    
+    pBar = pBar - f*INV(df);
+    
+    Z = DenL + tauL + pBar;
+    
+    Ssq = SQR(S1L) + SQR(S2L) + SQR(S3L);
+    
+    vsq = Ssq*INV(SQR(Z));
+    
+    WL = INV(sqrt(1 - vsq));
+    
+    rhoL = DenL*INV(WL);
+    
+    hL = Z*INV(rhoL)*INV(SQR(WL));
+    
+    epsiL = hL - (rhoL + pBar)*INV(rhoL);
+    
+    pEOS = epsiL*rhoL*(-1 + ToReal(gamma));
+    
+    f = -pBar + pEOS;
+    
+    cs = sqrt(epsiL*INV(hL)*(-1 + ToReal(gamma))*ToReal(gamma));
+    
+    df = -1 + vsq*SQR(cs);
+    
+    pBar = pBar - f*INV(df);
+    
+    Z = DenL + tauL + pBar;
+    
+    Ssq = SQR(S1L) + SQR(S2L) + SQR(S3L);
+    
+    vsq = Ssq*INV(SQR(Z));
+    
+    WL = INV(sqrt(1 - vsq));
+    
+    rhoL = DenL*INV(WL);
+    
+    hL = Z*INV(rhoL)*INV(SQR(WL));
+    
+    epsiL = hL - (rhoL + pBar)*INV(rhoL);
+    
+    pEOS = epsiL*rhoL*(-1 + ToReal(gamma));
+    
+    f = -pBar + pEOS;
+    
+    cs = sqrt(epsiL*INV(hL)*(-1 + ToReal(gamma))*ToReal(gamma));
+    
+    df = -1 + vsq*SQR(cs);
+    
+    pBar = pBar - f*INV(df);
+    
+    pL = pBar;
+    
+    CCTK_REAL v1L = S1L*INV(hL)*INV(rhoL)*INV(SQR(WL));
+    
+    CCTK_REAL v2L = S2L*INV(hL)*INV(rhoL)*INV(SQR(WL));
+    
+    CCTK_REAL v3L = S3L*INV(hL)*INV(rhoL)*INV(SQR(WL));
+    
+    /* Copy local copies back to grid functions */
+    epsi[index] = epsiL;
+    h[index] = hL;
+    p[index] = pL;
+    rho[index] = rhoL;
+    v1[index] = v1L;
+    v2[index] = v2L;
+    v3[index] = v3L;
+    W[index] = WL;
   }
+  CCTK_ENDLOOP3 (eulersr_cons_calc_primitives);
 }
 
 extern "C" void eulersr_cons_calc_primitives(CCTK_ARGUMENTS)
@@ -253,5 +248,25 @@ extern "C" void eulersr_cons_calc_primitives(CCTK_ARGUMENTS)
   DECLARE_CCTK_ARGUMENTS;
   DECLARE_CCTK_PARAMETERS;
   
+  
+  if (verbose > 1)
+  {
+    CCTK_VInfo(CCTK_THORNSTRING,"Entering eulersr_cons_calc_primitives_Body");
+  }
+  
+  if (cctk_iteration % eulersr_cons_calc_primitives_calc_every != eulersr_cons_calc_primitives_calc_offset)
+  {
+    return;
+  }
+  
+  const char *groups[] = {"EulerSR::Den_group","EulerSR::epsi_group","EulerSR::h_group","EulerSR::p_group","EulerSR::rho_group","EulerSR::S_group","EulerSR::tau_group","EulerSR::v_group","EulerSR::W_group"};
+  GenericFD_AssertGroupStorage(cctkGH, "eulersr_cons_calc_primitives", 9, groups);
+  
+  
   GenericFD_LoopOverEverything(cctkGH, &eulersr_cons_calc_primitives_Body);
+  
+  if (verbose > 1)
+  {
+    CCTK_VInfo(CCTK_THORNSTRING,"Leaving eulersr_cons_calc_primitives_Body");
+  }
 }
diff --git a/Examples/EulerSR/src/eulersr_cons_calc_reconstruct_1.cc b/Examples/EulerSR/src/eulersr_cons_calc_reconstruct_1.cc
index fd2111a..347fb53 100644
--- a/Examples/EulerSR/src/eulersr_cons_calc_reconstruct_1.cc
+++ b/Examples/EulerSR/src/eulersr_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 eulersr_cons_calc_reconstruct_1_SelectBCs(CCTK_ARGUMENTS)
   return;
 }
 
-static void eulersr_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 eulersr_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 eulersr_cons_calc_reconstruct_1_Body");
-  }
-  
-  if (cctk_iteration % eulersr_cons_calc_reconstruct_1_calc_every != eulersr_cons_calc_reconstruct_1_calc_offset)
-  {
-    return;
-  }
-  
-  const char *groups[] = {"EulerSR::epsi_group","EulerSR::epsi_lr_group","EulerSR::rho_group","EulerSR::rho_lr_group","EulerSR::v1_lr_group","EulerSR::v2_lr_group","EulerSR::v3_lr_group","EulerSR::v_group"};
-  GenericFD_AssertGroupStorage(cctkGH, "eulersr_cons_calc_reconstruct_1", 8, groups);
-  
-  GenericFD_EnsureStencilFits(cctkGH, "eulersr_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 eulersr_cons_calc_reconstruct_1_Body(cGH const * restrict const cctk
   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 (eulersr_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 epsiL = epsi[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 DiffPlus1epsi = DiffPlus1(&epsi[index]);
-        CCTK_REAL const DiffMinus1epsi = DiffMinus1(&epsi[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 = DiffMinus1epsi;
-        
-        slopeR = DiffPlus1epsi;
-        
-        slope = VanLeer(slopeL,slopeR);
-        
-        CCTK_REAL epsiLeftL = epsiL - 0.5*slope;
-        
-        CCTK_REAL epsiRightL = epsiL + 0.5*slope;
-        
-        /* Copy local copies back to grid functions */
-        epsiLeft[index] = epsiLeftL;
-        epsiRight[index] = epsiRightL;
-        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 epsiL = epsi[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 DiffPlus1epsi = DiffPlus1(&epsi[index]);
+    CCTK_REAL const DiffMinus1epsi = DiffMinus1(&epsi[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 = DiffMinus1epsi;
+    
+    slopeR = DiffPlus1epsi;
+    
+    slope = VanLeer(slopeL,slopeR);
+    
+    CCTK_REAL epsiLeftL = epsiL - 0.5*slope;
+    
+    CCTK_REAL epsiRightL = epsiL + 0.5*slope;
+    
+    /* Copy local copies back to grid functions */
+    epsiLeft[index] = epsiLeftL;
+    epsiRight[index] = epsiRightL;
+    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 (eulersr_cons_calc_reconstruct_1);
 }
 
 extern "C" void eulersr_cons_calc_reconstruct_1(CCTK_ARGUMENTS)
@@ -203,5 +196,26 @@ extern "C" void eulersr_cons_calc_reconstruct_1(CCTK_ARGUMENTS)
   DECLARE_CCTK_ARGUMENTS;
   DECLARE_CCTK_PARAMETERS;
   
+  
+  if (verbose > 1)
+  {
+    CCTK_VInfo(CCTK_THORNSTRING,"Entering eulersr_cons_calc_reconstruct_1_Body");
+  }
+  
+  if (cctk_iteration % eulersr_cons_calc_reconstruct_1_calc_every != eulersr_cons_calc_reconstruct_1_calc_offset)
+  {
+    return;
+  }
+  
+  const char *groups[] = {"EulerSR::epsi_group","EulerSR::epsi_lr_group","EulerSR::rho_group","EulerSR::rho_lr_group","EulerSR::v1_lr_group","EulerSR::v2_lr_group","EulerSR::v3_lr_group","EulerSR::v_group"};
+  GenericFD_AssertGroupStorage(cctkGH, "eulersr_cons_calc_reconstruct_1", 8, groups);
+  
+  GenericFD_EnsureStencilFits(cctkGH, "eulersr_cons_calc_reconstruct_1", 1, 1, 1);
+  
   GenericFD_LoopOverInterior(cctkGH, &eulersr_cons_calc_reconstruct_1_Body);
+  
+  if (verbose > 1)
+  {
+    CCTK_VInfo(CCTK_THORNSTRING,"Leaving eulersr_cons_calc_reconstruct_1_Body");
+  }
 }
diff --git a/Examples/EulerSR/src/eulersr_cons_calc_reconstruct_2.cc b/Examples/EulerSR/src/eulersr_cons_calc_reconstruct_2.cc
index 6725583..420741b 100644
--- a/Examples/EulerSR/src/eulersr_cons_calc_reconstruct_2.cc
+++ b/Examples/EulerSR/src/eulersr_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 eulersr_cons_calc_reconstruct_2_SelectBCs(CCTK_ARGUMENTS)
   return;
 }
 
-static void eulersr_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 eulersr_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 eulersr_cons_calc_reconstruct_2_Body");
-  }
-  
-  if (cctk_iteration % eulersr_cons_calc_reconstruct_2_calc_every != eulersr_cons_calc_reconstruct_2_calc_offset)
-  {
-    return;
-  }
-  
-  const char *groups[] = {"EulerSR::epsi_group","EulerSR::epsi_lr_group","EulerSR::rho_group","EulerSR::rho_lr_group","EulerSR::v1_lr_group","EulerSR::v2_lr_group","EulerSR::v3_lr_group","EulerSR::v_group"};
-  GenericFD_AssertGroupStorage(cctkGH, "eulersr_cons_calc_reconstruct_2", 8, groups);
-  
-  GenericFD_EnsureStencilFits(cctkGH, "eulersr_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 eulersr_cons_calc_reconstruct_2_Body(cGH const * restrict const cctk
   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 (eulersr_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 epsiL = epsi[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 DiffPlus2epsi = DiffPlus2(&epsi[index]);
-        CCTK_REAL const DiffMinus2epsi = DiffMinus2(&epsi[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 = DiffMinus2epsi;
-        
-        slopeR = DiffPlus2epsi;
-        
-        slope = VanLeer(slopeL,slopeR);
-        
-        CCTK_REAL epsiLeftL = epsiL - 0.5*slope;
-        
-        CCTK_REAL epsiRightL = epsiL + 0.5*slope;
-        
-        /* Copy local copies back to grid functions */
-        epsiLeft[index] = epsiLeftL;
-        epsiRight[index] = epsiRightL;
-        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 epsiL = epsi[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 DiffPlus2epsi = DiffPlus2(&epsi[index]);
+    CCTK_REAL const DiffMinus2epsi = DiffMinus2(&epsi[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 = DiffMinus2epsi;
+    
+    slopeR = DiffPlus2epsi;
+    
+    slope = VanLeer(slopeL,slopeR);
+    
+    CCTK_REAL epsiLeftL = epsiL - 0.5*slope;
+    
+    CCTK_REAL epsiRightL = epsiL + 0.5*slope;
+    
+    /* Copy local copies back to grid functions */
+    epsiLeft[index] = epsiLeftL;
+    epsiRight[index] = epsiRightL;
+    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 (eulersr_cons_calc_reconstruct_2);
 }
 
 extern "C" void eulersr_cons_calc_reconstruct_2(CCTK_ARGUMENTS)
@@ -203,5 +196,26 @@ extern "C" void eulersr_cons_calc_reconstruct_2(CCTK_ARGUMENTS)
   DECLARE_CCTK_ARGUMENTS;
   DECLARE_CCTK_PARAMETERS;
   
+  
+  if (verbose > 1)
+  {
+    CCTK_VInfo(CCTK_THORNSTRING,"Entering eulersr_cons_calc_reconstruct_2_Body");
+  }
+  
+  if (cctk_iteration % eulersr_cons_calc_reconstruct_2_calc_every != eulersr_cons_calc_reconstruct_2_calc_offset)
+  {
+    return;
+  }
+  
+  const char *groups[] = {"EulerSR::epsi_group","EulerSR::epsi_lr_group","EulerSR::rho_group","EulerSR::rho_lr_group","EulerSR::v1_lr_group","EulerSR::v2_lr_group","EulerSR::v3_lr_group","EulerSR::v_group"};
+  GenericFD_AssertGroupStorage(cctkGH, "eulersr_cons_calc_reconstruct_2", 8, groups);
+  
+  GenericFD_EnsureStencilFits(cctkGH, "eulersr_cons_calc_reconstruct_2", 1, 1, 1);
+  
   GenericFD_LoopOverInterior(cctkGH, &eulersr_cons_calc_reconstruct_2_Body);
+  
+  if (verbose > 1)
+  {
+    CCTK_VInfo(CCTK_THORNSTRING,"Leaving eulersr_cons_calc_reconstruct_2_Body");
+  }
 }
diff --git a/Examples/EulerSR/src/eulersr_cons_calc_reconstruct_3.cc b/Examples/EulerSR/src/eulersr_cons_calc_reconstruct_3.cc
index 7a913d0..d0d6b42 100644
--- a/Examples/EulerSR/src/eulersr_cons_calc_reconstruct_3.cc
+++ b/Examples/EulerSR/src/eulersr_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 eulersr_cons_calc_reconstruct_3_SelectBCs(CCTK_ARGUMENTS)
   return;
 }
 
-static void eulersr_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 eulersr_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 eulersr_cons_calc_reconstruct_3_Body");
-  }
-  
-  if (cctk_iteration % eulersr_cons_calc_reconstruct_3_calc_every != eulersr_cons_calc_reconstruct_3_calc_offset)
-  {
-    return;
-  }
-  
-  const char *groups[] = {"EulerSR::epsi_group","EulerSR::epsi_lr_group","EulerSR::rho_group","EulerSR::rho_lr_group","EulerSR::v1_lr_group","EulerSR::v2_lr_group","EulerSR::v3_lr_group","EulerSR::v_group"};
-  GenericFD_AssertGroupStorage(cctkGH, "eulersr_cons_calc_reconstruct_3", 8, groups);
-  
-  GenericFD_EnsureStencilFits(cctkGH, "eulersr_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 eulersr_cons_calc_reconstruct_3_Body(cGH const * restrict const cctk
   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 (eulersr_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 epsiL = epsi[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 DiffPlus3epsi = DiffPlus3(&epsi[index]);
-        CCTK_REAL const DiffMinus3epsi = DiffMinus3(&epsi[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 = DiffMinus3epsi;
-        
-        slopeR = DiffPlus3epsi;
-        
-        slope = VanLeer(slopeL,slopeR);
-        
-        CCTK_REAL epsiLeftL = epsiL - 0.5*slope;
-        
-        CCTK_REAL epsiRightL = epsiL + 0.5*slope;
-        
-        /* Copy local copies back to grid functions */
-        epsiLeft[index] = epsiLeftL;
-        epsiRight[index] = epsiRightL;
-        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 epsiL = epsi[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 DiffPlus3epsi = DiffPlus3(&epsi[index]);
+    CCTK_REAL const DiffMinus3epsi = DiffMinus3(&epsi[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 = DiffMinus3epsi;
+    
+    slopeR = DiffPlus3epsi;
+    
+    slope = VanLeer(slopeL,slopeR);
+    
+    CCTK_REAL epsiLeftL = epsiL - 0.5*slope;
+    
+    CCTK_REAL epsiRightL = epsiL + 0.5*slope;
+    
+    /* Copy local copies back to grid functions */
+    epsiLeft[index] = epsiLeftL;
+    epsiRight[index] = epsiRightL;
+    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 (eulersr_cons_calc_reconstruct_3);
 }
 
 extern "C" void eulersr_cons_calc_reconstruct_3(CCTK_ARGUMENTS)
@@ -203,5 +196,26 @@ extern "C" void eulersr_cons_calc_reconstruct_3(CCTK_ARGUMENTS)
   DECLARE_CCTK_ARGUMENTS;
   DECLARE_CCTK_PARAMETERS;
   
+  
+  if (verbose > 1)
+  {
+    CCTK_VInfo(CCTK_THORNSTRING,"Entering eulersr_cons_calc_reconstruct_3_Body");
+  }
+  
+  if (cctk_iteration % eulersr_cons_calc_reconstruct_3_calc_every != eulersr_cons_calc_reconstruct_3_calc_offset)
+  {
+    return;
+  }
+  
+  const char *groups[] = {"EulerSR::epsi_group","EulerSR::epsi_lr_group","EulerSR::rho_group","EulerSR::rho_lr_group","EulerSR::v1_lr_group","EulerSR::v2_lr_group","EulerSR::v3_lr_group","EulerSR::v_group"};
+  GenericFD_AssertGroupStorage(cctkGH, "eulersr_cons_calc_reconstruct_3", 8, groups);
+  
+  GenericFD_EnsureStencilFits(cctkGH, "eulersr_cons_calc_reconstruct_3", 1, 1, 1);
+  
   GenericFD_LoopOverInterior(cctkGH, &eulersr_cons_calc_reconstruct_3_Body);
+  
+  if (verbose > 1)
+  {
+    CCTK_VInfo(CCTK_THORNSTRING,"Leaving eulersr_cons_calc_reconstruct_3_Body");
+  }
 }
diff --git a/Examples/EulerSR/src/eulersr_cons_calc_rhs_1.cc b/Examples/EulerSR/src/eulersr_cons_calc_rhs_1.cc
index 0dee465..f5ca0d1 100644
--- a/Examples/EulerSR/src/eulersr_cons_calc_rhs_1.cc
+++ b/Examples/EulerSR/src/eulersr_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 eulersr_cons_calc_rhs_1_SelectBCs(CCTK_ARGUMENTS)
   return;
 }
 
-static void eulersr_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 eulersr_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 eulersr_cons_calc_rhs_1_Body");
-  }
-  
-  if (cctk_iteration % eulersr_cons_calc_rhs_1_calc_every != eulersr_cons_calc_rhs_1_calc_offset)
-  {
-    return;
-  }
-  
-  const char *groups[] = {"EulerSR::Den_flux_group","EulerSR::Den_grouprhs","EulerSR::S1_flux_group","EulerSR::S2_flux_group","EulerSR::S3_flux_group","EulerSR::S_grouprhs","EulerSR::tau_flux_group","EulerSR::tau_grouprhs"};
-  GenericFD_AssertGroupStorage(cctkGH, "eulersr_cons_calc_rhs_1", 8, groups);
-  
-  GenericFD_EnsureStencilFits(cctkGH, "eulersr_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 eulersr_cons_calc_rhs_1_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 (eulersr_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 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];
-        CCTK_REAL tauFluxL = tauFlux[index];
-        CCTK_REAL taurhsL = taurhs[index];
-        
-        
-        /* Include user supplied include files */
-        
-        /* Precompute derivatives */
-        CCTK_REAL const PDplus1DenFlux = PDplus1(&DenFlux[index]);
-        CCTK_REAL const PDplus1S1Flux = PDplus1(&S1Flux[index]);
-        CCTK_REAL const PDplus1S2Flux = PDplus1(&S2Flux[index]);
-        CCTK_REAL const PDplus1S3Flux = PDplus1(&S3Flux[index]);
-        CCTK_REAL const PDplus1tauFlux = PDplus1(&tauFlux[index]);
-        
-        /* Calculate temporaries and grid functions */
-        DenrhsL = DenrhsL - PDplus1DenFlux;
-        
-        S1rhsL = -PDplus1S1Flux + S1rhsL;
-        
-        S2rhsL = -PDplus1S2Flux + S2rhsL;
-        
-        S3rhsL = -PDplus1S3Flux + S3rhsL;
-        
-        taurhsL = -PDplus1tauFlux + taurhsL;
-        
-        /* Copy local copies back to grid functions */
-        Denrhs[index] = DenrhsL;
-        S1rhs[index] = S1rhsL;
-        S2rhs[index] = S2rhsL;
-        S3rhs[index] = S3rhsL;
-        taurhs[index] = taurhsL;
-      }
-    }
+    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 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];
+    CCTK_REAL tauFluxL = tauFlux[index];
+    CCTK_REAL taurhsL = taurhs[index];
+    
+    
+    /* Include user supplied include files */
+    
+    /* Precompute derivatives */
+    CCTK_REAL const PDplus1DenFlux = PDplus1(&DenFlux[index]);
+    CCTK_REAL const PDplus1S1Flux = PDplus1(&S1Flux[index]);
+    CCTK_REAL const PDplus1S2Flux = PDplus1(&S2Flux[index]);
+    CCTK_REAL const PDplus1S3Flux = PDplus1(&S3Flux[index]);
+    CCTK_REAL const PDplus1tauFlux = PDplus1(&tauFlux[index]);
+    
+    /* Calculate temporaries and grid functions */
+    DenrhsL = DenrhsL - PDplus1DenFlux;
+    
+    S1rhsL = S1rhsL - PDplus1S1Flux;
+    
+    S2rhsL = S2rhsL - PDplus1S2Flux;
+    
+    S3rhsL = S3rhsL - PDplus1S3Flux;
+    
+    taurhsL = taurhsL - PDplus1tauFlux;
+    
+    /* Copy local copies back to grid functions */
+    Denrhs[index] = DenrhsL;
+    S1rhs[index] = S1rhsL;
+    S2rhs[index] = S2rhsL;
+    S3rhs[index] = S3rhsL;
+    taurhs[index] = taurhsL;
   }
+  CCTK_ENDLOOP3 (eulersr_cons_calc_rhs_1);
 }
 
 extern "C" void eulersr_cons_calc_rhs_1(CCTK_ARGUMENTS)
@@ -152,5 +145,26 @@ extern "C" void eulersr_cons_calc_rhs_1(CCTK_ARGUMENTS)
   DECLARE_CCTK_ARGUMENTS;
   DECLARE_CCTK_PARAMETERS;
   
+  
+  if (verbose > 1)
+  {
+    CCTK_VInfo(CCTK_THORNSTRING,"Entering eulersr_cons_calc_rhs_1_Body");
+  }
+  
+  if (cctk_iteration % eulersr_cons_calc_rhs_1_calc_every != eulersr_cons_calc_rhs_1_calc_offset)
+  {
+    return;
+  }
+  
+  const char *groups[] = {"EulerSR::Den_flux_group","EulerSR::Den_grouprhs","EulerSR::S1_flux_group","EulerSR::S2_flux_group","EulerSR::S3_flux_group","EulerSR::S_grouprhs","EulerSR::tau_flux_group","EulerSR::tau_grouprhs"};
+  GenericFD_AssertGroupStorage(cctkGH, "eulersr_cons_calc_rhs_1", 8, groups);
+  
+  GenericFD_EnsureStencilFits(cctkGH, "eulersr_cons_calc_rhs_1", 1, 1, 1);
+  
   GenericFD_LoopOverInterior(cctkGH, &eulersr_cons_calc_rhs_1_Body);
+  
+  if (verbose > 1)
+  {
+    CCTK_VInfo(CCTK_THORNSTRING,"Leaving eulersr_cons_calc_rhs_1_Body");
+  }
 }
diff --git a/Examples/EulerSR/src/eulersr_cons_calc_rhs_2.cc b/Examples/EulerSR/src/eulersr_cons_calc_rhs_2.cc
index f8a470b..33d6d0e 100644
--- a/Examples/EulerSR/src/eulersr_cons_calc_rhs_2.cc
+++ b/Examples/EulerSR/src/eulersr_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 eulersr_cons_calc_rhs_2_SelectBCs(CCTK_ARGUMENTS)
   return;
 }
 
-static void eulersr_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 eulersr_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 eulersr_cons_calc_rhs_2_Body");
-  }
-  
-  if (cctk_iteration % eulersr_cons_calc_rhs_2_calc_every != eulersr_cons_calc_rhs_2_calc_offset)
-  {
-    return;
-  }
-  
-  const char *groups[] = {"EulerSR::Den_flux_group","EulerSR::Den_grouprhs","EulerSR::S1_flux_group","EulerSR::S2_flux_group","EulerSR::S3_flux_group","EulerSR::S_grouprhs","EulerSR::tau_flux_group","EulerSR::tau_grouprhs"};
-  GenericFD_AssertGroupStorage(cctkGH, "eulersr_cons_calc_rhs_2", 8, groups);
-  
-  GenericFD_EnsureStencilFits(cctkGH, "eulersr_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 eulersr_cons_calc_rhs_2_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 (eulersr_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 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];
-        CCTK_REAL tauFluxL = tauFlux[index];
-        CCTK_REAL taurhsL = taurhs[index];
-        
-        
-        /* Include user supplied include files */
-        
-        /* Precompute derivatives */
-        CCTK_REAL const PDplus2DenFlux = PDplus2(&DenFlux[index]);
-        CCTK_REAL const PDplus2S1Flux = PDplus2(&S1Flux[index]);
-        CCTK_REAL const PDplus2S2Flux = PDplus2(&S2Flux[index]);
-        CCTK_REAL const PDplus2S3Flux = PDplus2(&S3Flux[index]);
-        CCTK_REAL const PDplus2tauFlux = PDplus2(&tauFlux[index]);
-        
-        /* Calculate temporaries and grid functions */
-        DenrhsL = DenrhsL - PDplus2DenFlux;
-        
-        S1rhsL = -PDplus2S1Flux + S1rhsL;
-        
-        S2rhsL = -PDplus2S2Flux + S2rhsL;
-        
-        S3rhsL = -PDplus2S3Flux + S3rhsL;
-        
-        taurhsL = -PDplus2tauFlux + taurhsL;
-        
-        /* Copy local copies back to grid functions */
-        Denrhs[index] = DenrhsL;
-        S1rhs[index] = S1rhsL;
-        S2rhs[index] = S2rhsL;
-        S3rhs[index] = S3rhsL;
-        taurhs[index] = taurhsL;
-      }
-    }
+    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 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];
+    CCTK_REAL tauFluxL = tauFlux[index];
+    CCTK_REAL taurhsL = taurhs[index];
+    
+    
+    /* Include user supplied include files */
+    
+    /* Precompute derivatives */
+    CCTK_REAL const PDplus2DenFlux = PDplus2(&DenFlux[index]);
+    CCTK_REAL const PDplus2S1Flux = PDplus2(&S1Flux[index]);
+    CCTK_REAL const PDplus2S2Flux = PDplus2(&S2Flux[index]);
+    CCTK_REAL const PDplus2S3Flux = PDplus2(&S3Flux[index]);
+    CCTK_REAL const PDplus2tauFlux = PDplus2(&tauFlux[index]);
+    
+    /* Calculate temporaries and grid functions */
+    DenrhsL = DenrhsL - PDplus2DenFlux;
+    
+    S1rhsL = S1rhsL - PDplus2S1Flux;
+    
+    S2rhsL = S2rhsL - PDplus2S2Flux;
+    
+    S3rhsL = S3rhsL - PDplus2S3Flux;
+    
+    taurhsL = taurhsL - PDplus2tauFlux;
+    
+    /* Copy local copies back to grid functions */
+    Denrhs[index] = DenrhsL;
+    S1rhs[index] = S1rhsL;
+    S2rhs[index] = S2rhsL;
+    S3rhs[index] = S3rhsL;
+    taurhs[index] = taurhsL;
   }
+  CCTK_ENDLOOP3 (eulersr_cons_calc_rhs_2);
 }
 
 extern "C" void eulersr_cons_calc_rhs_2(CCTK_ARGUMENTS)
@@ -152,5 +145,26 @@ extern "C" void eulersr_cons_calc_rhs_2(CCTK_ARGUMENTS)
   DECLARE_CCTK_ARGUMENTS;
   DECLARE_CCTK_PARAMETERS;
   
+  
+  if (verbose > 1)
+  {
+    CCTK_VInfo(CCTK_THORNSTRING,"Entering eulersr_cons_calc_rhs_2_Body");
+  }
+  
+  if (cctk_iteration % eulersr_cons_calc_rhs_2_calc_every != eulersr_cons_calc_rhs_2_calc_offset)
+  {
+    return;
+  }
+  
+  const char *groups[] = {"EulerSR::Den_flux_group","EulerSR::Den_grouprhs","EulerSR::S1_flux_group","EulerSR::S2_flux_group","EulerSR::S3_flux_group","EulerSR::S_grouprhs","EulerSR::tau_flux_group","EulerSR::tau_grouprhs"};
+  GenericFD_AssertGroupStorage(cctkGH, "eulersr_cons_calc_rhs_2", 8, groups);
+  
+  GenericFD_EnsureStencilFits(cctkGH, "eulersr_cons_calc_rhs_2", 1, 1, 1);
+  
   GenericFD_LoopOverInterior(cctkGH, &eulersr_cons_calc_rhs_2_Body);
+  
+  if (verbose > 1)
+  {
+    CCTK_VInfo(CCTK_THORNSTRING,"Leaving eulersr_cons_calc_rhs_2_Body");
+  }
 }
diff --git a/Examples/EulerSR/src/eulersr_cons_calc_rhs_3.cc b/Examples/EulerSR/src/eulersr_cons_calc_rhs_3.cc
index 2b33c05..6a288fe 100644
--- a/Examples/EulerSR/src/eulersr_cons_calc_rhs_3.cc
+++ b/Examples/EulerSR/src/eulersr_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 eulersr_cons_calc_rhs_3_SelectBCs(CCTK_ARGUMENTS)
   return;
 }
 
-static void eulersr_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 eulersr_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 eulersr_cons_calc_rhs_3_Body");
-  }
-  
-  if (cctk_iteration % eulersr_cons_calc_rhs_3_calc_every != eulersr_cons_calc_rhs_3_calc_offset)
-  {
-    return;
-  }
-  
-  const char *groups[] = {"EulerSR::Den_flux_group","EulerSR::Den_grouprhs","EulerSR::S1_flux_group","EulerSR::S2_flux_group","EulerSR::S3_flux_group","EulerSR::S_grouprhs","EulerSR::tau_flux_group","EulerSR::tau_grouprhs"};
-  GenericFD_AssertGroupStorage(cctkGH, "eulersr_cons_calc_rhs_3", 8, groups);
-  
-  GenericFD_EnsureStencilFits(cctkGH, "eulersr_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 eulersr_cons_calc_rhs_3_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 (eulersr_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 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];
-        CCTK_REAL tauFluxL = tauFlux[index];
-        CCTK_REAL taurhsL = taurhs[index];
-        
-        
-        /* Include user supplied include files */
-        
-        /* Precompute derivatives */
-        CCTK_REAL const PDplus3DenFlux = PDplus3(&DenFlux[index]);
-        CCTK_REAL const PDplus3S1Flux = PDplus3(&S1Flux[index]);
-        CCTK_REAL const PDplus3S2Flux = PDplus3(&S2Flux[index]);
-        CCTK_REAL const PDplus3S3Flux = PDplus3(&S3Flux[index]);
-        CCTK_REAL const PDplus3tauFlux = PDplus3(&tauFlux[index]);
-        
-        /* Calculate temporaries and grid functions */
-        DenrhsL = DenrhsL - PDplus3DenFlux;
-        
-        S1rhsL = -PDplus3S1Flux + S1rhsL;
-        
-        S2rhsL = -PDplus3S2Flux + S2rhsL;
-        
-        S3rhsL = -PDplus3S3Flux + S3rhsL;
-        
-        taurhsL = -PDplus3tauFlux + taurhsL;
-        
-        /* Copy local copies back to grid functions */
-        Denrhs[index] = DenrhsL;
-        S1rhs[index] = S1rhsL;
-        S2rhs[index] = S2rhsL;
-        S3rhs[index] = S3rhsL;
-        taurhs[index] = taurhsL;
-      }
-    }
+    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 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];
+    CCTK_REAL tauFluxL = tauFlux[index];
+    CCTK_REAL taurhsL = taurhs[index];
+    
+    
+    /* Include user supplied include files */
+    
+    /* Precompute derivatives */
+    CCTK_REAL const PDplus3DenFlux = PDplus3(&DenFlux[index]);
+    CCTK_REAL const PDplus3S1Flux = PDplus3(&S1Flux[index]);
+    CCTK_REAL const PDplus3S2Flux = PDplus3(&S2Flux[index]);
+    CCTK_REAL const PDplus3S3Flux = PDplus3(&S3Flux[index]);
+    CCTK_REAL const PDplus3tauFlux = PDplus3(&tauFlux[index]);
+    
+    /* Calculate temporaries and grid functions */
+    DenrhsL = DenrhsL - PDplus3DenFlux;
+    
+    S1rhsL = S1rhsL - PDplus3S1Flux;
+    
+    S2rhsL = S2rhsL - PDplus3S2Flux;
+    
+    S3rhsL = S3rhsL - PDplus3S3Flux;
+    
+    taurhsL = taurhsL - PDplus3tauFlux;
+    
+    /* Copy local copies back to grid functions */
+    Denrhs[index] = DenrhsL;
+    S1rhs[index] = S1rhsL;
+    S2rhs[index] = S2rhsL;
+    S3rhs[index] = S3rhsL;
+    taurhs[index] = taurhsL;
   }
+  CCTK_ENDLOOP3 (eulersr_cons_calc_rhs_3);
 }
 
 extern "C" void eulersr_cons_calc_rhs_3(CCTK_ARGUMENTS)
@@ -152,5 +145,26 @@ extern "C" void eulersr_cons_calc_rhs_3(CCTK_ARGUMENTS)
   DECLARE_CCTK_ARGUMENTS;
   DECLARE_CCTK_PARAMETERS;
   
+  
+  if (verbose > 1)
+  {
+    CCTK_VInfo(CCTK_THORNSTRING,"Entering eulersr_cons_calc_rhs_3_Body");
+  }
+  
+  if (cctk_iteration % eulersr_cons_calc_rhs_3_calc_every != eulersr_cons_calc_rhs_3_calc_offset)
+  {
+    return;
+  }
+  
+  const char *groups[] = {"EulerSR::Den_flux_group","EulerSR::Den_grouprhs","EulerSR::S1_flux_group","EulerSR::S2_flux_group","EulerSR::S3_flux_group","EulerSR::S_grouprhs","EulerSR::tau_flux_group","EulerSR::tau_grouprhs"};
+  GenericFD_AssertGroupStorage(cctkGH, "eulersr_cons_calc_rhs_3", 8, groups);
+  
+  GenericFD_EnsureStencilFits(cctkGH, "eulersr_cons_calc_rhs_3", 1, 1, 1);
+  
   GenericFD_LoopOverInterior(cctkGH, &eulersr_cons_calc_rhs_3_Body);
+  
+  if (verbose > 1)
+  {
+    CCTK_VInfo(CCTK_THORNSTRING,"Leaving eulersr_cons_calc_rhs_3_Body");
+  }
 }
diff --git a/Examples/EulerSR/src/eulersr_cons_calc_zero_rhs.cc b/Examples/EulerSR/src/eulersr_cons_calc_zero_rhs.cc
index 1919515..6a93d61 100644
--- a/Examples/EulerSR/src/eulersr_cons_calc_zero_rhs.cc
+++ b/Examples/EulerSR/src/eulersr_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 eulersr_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 eulersr_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 eulersr_cons_calc_zero_rhs_Body");
-  }
-  
-  if (cctk_iteration % eulersr_cons_calc_zero_rhs_calc_every != eulersr_cons_calc_zero_rhs_calc_offset)
-  {
-    return;
-  }
-  
-  const char *groups[] = {"EulerSR::Den_grouprhs","EulerSR::S_grouprhs","EulerSR::tau_grouprhs"};
-  GenericFD_AssertGroupStorage(cctkGH, "eulersr_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 eulersr_cons_calc_zero_rhs_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 (eulersr_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;
-        
-        S1rhsL = 0;
-        
-        S1rhsL = 0;
-        
-        CCTK_REAL S2rhsL = 0;
-        
-        S2rhsL = 0;
-        
-        S2rhsL = 0;
-        
-        CCTK_REAL S3rhsL = 0;
-        
-        S3rhsL = 0;
-        
-        S3rhsL = 0;
-        
-        CCTK_REAL taurhsL = 0;
-        
-        taurhsL = 0;
-        
-        taurhsL = 0;
-        
-        /* Copy local copies back to grid functions */
-        Denrhs[index] = DenrhsL;
-        S1rhs[index] = S1rhsL;
-        S2rhs[index] = S2rhsL;
-        S3rhs[index] = S3rhsL;
-        taurhs[index] = taurhsL;
-      }
-    }
+    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;
+    
+    S1rhsL = 0;
+    
+    S1rhsL = 0;
+    
+    CCTK_REAL S2rhsL = 0;
+    
+    S2rhsL = 0;
+    
+    S2rhsL = 0;
+    
+    CCTK_REAL S3rhsL = 0;
+    
+    S3rhsL = 0;
+    
+    S3rhsL = 0;
+    
+    CCTK_REAL taurhsL = 0;
+    
+    taurhsL = 0;
+    
+    taurhsL = 0;
+    
+    /* Copy local copies back to grid functions */
+    Denrhs[index] = DenrhsL;
+    S1rhs[index] = S1rhsL;
+    S2rhs[index] = S2rhsL;
+    S3rhs[index] = S3rhsL;
+    taurhs[index] = taurhsL;
   }
+  CCTK_ENDLOOP3 (eulersr_cons_calc_zero_rhs);
 }
 
 extern "C" void eulersr_cons_calc_zero_rhs(CCTK_ARGUMENTS)
@@ -138,5 +132,25 @@ extern "C" void eulersr_cons_calc_zero_rhs(CCTK_ARGUMENTS)
   DECLARE_CCTK_ARGUMENTS;
   DECLARE_CCTK_PARAMETERS;
   
+  
+  if (verbose > 1)
+  {
+    CCTK_VInfo(CCTK_THORNSTRING,"Entering eulersr_cons_calc_zero_rhs_Body");
+  }
+  
+  if (cctk_iteration % eulersr_cons_calc_zero_rhs_calc_every != eulersr_cons_calc_zero_rhs_calc_offset)
+  {
+    return;
+  }
+  
+  const char *groups[] = {"EulerSR::Den_grouprhs","EulerSR::S_grouprhs","EulerSR::tau_grouprhs"};
+  GenericFD_AssertGroupStorage(cctkGH, "eulersr_cons_calc_zero_rhs", 3, groups);
+  
+  
   GenericFD_LoopOverEverything(cctkGH, &eulersr_cons_calc_zero_rhs_Body);
+  
+  if (verbose > 1)
+  {
+    CCTK_VInfo(CCTK_THORNSTRING,"Leaving eulersr_cons_calc_zero_rhs_Body");
+  }
 }
diff --git a/Examples/EulerSR/src/eulersr_initial_shock.cc b/Examples/EulerSR/src/eulersr_initial_shock.cc
index 2881056..73eb475 100644
--- a/Examples/EulerSR/src/eulersr_initial_shock.cc
+++ b/Examples/EulerSR/src/eulersr_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 eulersr_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 eulersr_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 eulersr_initial_shock_Body");
-  }
-  
-  if (cctk_iteration % eulersr_initial_shock_calc_every != eulersr_initial_shock_calc_offset)
-  {
-    return;
-  }
-  
-  const char *groups[] = {"EulerSR::epsi_group","grid::coordinates","EulerSR::rho_group","EulerSR::v_group"};
-  GenericFD_AssertGroupStorage(cctkGH, "eulersr_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,51 +62,57 @@ static void eulersr_initial_shock_Body(cGH const * restrict const cctkGH, int co
   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 (eulersr_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 X = xL;
-        
-        CCTK_REAL rhoL = ToReal(rhoL0) + StepFunction(X)*(-ToReal(rhoL0) + 
-          ToReal(rhoR0));
-        
-        CCTK_REAL v1L = ToReal(vL0) + StepFunction(X)*(-ToReal(vL0) + 
-          ToReal(vR0));
-        
-        CCTK_REAL v2L = ToReal(vL0) + StepFunction(X)*(-ToReal(vL0) + 
-          ToReal(vR0));
-        
-        CCTK_REAL v3L = ToReal(vL0) + StepFunction(X)*(-ToReal(vL0) + 
-          ToReal(vR0));
-        
-        CCTK_REAL epsiL = ToReal(epsiL0) + StepFunction(X)*(-ToReal(epsiL0) + 
-          ToReal(epsiR0));
-        
-        /* Copy local copies back to grid functions */
-        epsi[index] = epsiL;
-        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 X = xL;
+    
+    CCTK_REAL rhoL = ToReal(rhoL0) + StepFunction(X)*(-ToReal(rhoL0) + 
+      ToReal(rhoR0));
+    
+    CCTK_REAL v1L = ToReal(vL0) + StepFunction(X)*(-ToReal(vL0) + 
+      ToReal(vR0));
+    
+    CCTK_REAL v2L = ToReal(vL0) + StepFunction(X)*(-ToReal(vL0) + 
+      ToReal(vR0));
+    
+    CCTK_REAL v3L = ToReal(vL0) + StepFunction(X)*(-ToReal(vL0) + 
+      ToReal(vR0));
+    
+    CCTK_REAL epsiL = ToReal(epsiL0) + StepFunction(X)*(-ToReal(epsiL0) + 
+      ToReal(epsiR0));
+    
+    /* Copy local copies back to grid functions */
+    epsi[index] = epsiL;
+    rho[index] = rhoL;
+    v1[index] = v1L;
+    v2[index] = v2L;
+    v3[index] = v3L;
   }
+  CCTK_ENDLOOP3 (eulersr_initial_shock);
 }
 
 extern "C" void eulersr_initial_shock(CCTK_ARGUMENTS)
@@ -126,5 +120,25 @@ extern "C" void eulersr_initial_shock(CCTK_ARGUMENTS)
   DECLARE_CCTK_ARGUMENTS;
   DECLARE_CCTK_PARAMETERS;
   
+  
+  if (verbose > 1)
+  {
+    CCTK_VInfo(CCTK_THORNSTRING,"Entering eulersr_initial_shock_Body");
+  }
+  
+  if (cctk_iteration % eulersr_initial_shock_calc_every != eulersr_initial_shock_calc_offset)
+  {
+    return;
+  }
+  
+  const char *groups[] = {"EulerSR::epsi_group","grid::coordinates","EulerSR::rho_group","EulerSR::v_group"};
+  GenericFD_AssertGroupStorage(cctkGH, "eulersr_initial_shock", 4, groups);
+  
+  
   GenericFD_LoopOverEverything(cctkGH, &eulersr_initial_shock_Body);
+  
+  if (verbose > 1)
+  {
+    CCTK_VInfo(CCTK_THORNSTRING,"Leaving eulersr_initial_shock_Body");
+  }
 }