Regenerate examples (except Advect)

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

12 files changed, 533 insertions, 299 deletions

diff --git a/Examples/Burgers/configuration.ccl b/Examples/Burgers/configuration.ccl
index 023aac7..0a66ec2 100644
--- a/Examples/Burgers/configuration.ccl
+++ b/Examples/Burgers/configuration.ccl
@@ -1,3 +1,6 @@
 # File produced by Kranc
 
 REQUIRES GenericFD
+OPTIONAL LoopControl
+{
+}
diff --git a/Examples/Burgers/interface.ccl b/Examples/Burgers/interface.ccl
index 873ed69..fadd7cc 100644
--- a/Examples/Burgers/interface.ccl
+++ b/Examples/Burgers/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/Burgers/param.ccl b/Examples/Burgers/param.ccl
index a94fb3d..0b3226d 100644
--- a/Examples/Burgers/param.ccl
+++ b/Examples/Burgers/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
@@ -71,6 +72,12 @@ CCTK_INT Burgers_MaxNumEvolvedVars "Number of evolved variables used by this tho
 } 1
 
 restricted:
+CCTK_INT Burgers_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 :: ""
@@ -83,6 +90,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 burgers_initial_sine_calc_every "burgers_initial_sine_calc_every" STEERABLE=ALWAYS
 {
   *:* :: ""
diff --git a/Examples/Burgers/schedule.ccl b/Examples/Burgers/schedule.ccl
index 73893b1..2384e1a 100644
--- a/Examples/Burgers/schedule.ccl
+++ b/Examples/Burgers/schedule.ccl
@@ -1,11 +1,20 @@
 # File produced by Kranc
 
 
-STORAGE: uF_group[1]
+if (other_timelevels == 1)
+{
+  STORAGE: uF_group[1]
+}
 
-STORAGE: uLeft_group[1]
+if (other_timelevels == 1)
+{
+  STORAGE: uLeft_group[1]
+}
 
-STORAGE: uR_group[1]
+if (other_timelevels == 1)
+{
+  STORAGE: uR_group[1]
+}
 
 if (timelevels == 1)
 {
@@ -102,12 +111,7 @@ schedule group burgers_reconstruct_1_bc_group in burgers_reconstruct_1_group aft
   # no language specified
 } "burgers_reconstruct_1"
 
-schedule group burgers_reconstruct_1_bc_group at CCTK_POSTRESTRICT
-{
-  # no language specified
-} "burgers_reconstruct_1"
-
-schedule group burgers_reconstruct_1_bc_group at CCTK_POSTRESTRICTINITIAL
+schedule group burgers_reconstruct_1_bc_group in MoL_PseudoEvolutionBoundaries after MoL_PostStep
 {
   # no language specified
 } "burgers_reconstruct_1"
@@ -139,12 +143,7 @@ schedule group burgers_flux_1_bc_group in burgers_flux_1_group after burgers_flu
   # no language specified
 } "burgers_flux_1"
 
-schedule group burgers_flux_1_bc_group at CCTK_POSTRESTRICT
-{
-  # no language specified
-} "burgers_flux_1"
-
-schedule group burgers_flux_1_bc_group at CCTK_POSTRESTRICTINITIAL
+schedule group burgers_flux_1_bc_group in MoL_PseudoEvolutionBoundaries after MoL_PostStep
 {
   # no language specified
 } "burgers_flux_1"
diff --git a/Examples/Burgers/src/Differencing.h b/Examples/Burgers/src/Differencing.h
index 94f5300..dec7ce7 100644
--- a/Examples/Burgers/src/Differencing.h
+++ b/Examples/Burgers/src/Differencing.h
@@ -1,12 +1,144 @@
-#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)]))
-#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)])
+#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
+
+#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
+
diff --git a/Examples/Burgers/src/RegisterMoL.cc b/Examples/Burgers/src/RegisterMoL.cc
index d29d0eb..3863790 100644
--- a/Examples/Burgers/src/RegisterMoL.cc
+++ b/Examples/Burgers/src/RegisterMoL.cc
@@ -13,5 +13,7 @@ extern "C" void Burgers_RegisterVars(CCTK_ARGUMENTS)
   
   /* Register all the evolved grid functions with MoL */
   ierr += MoLRegisterEvolved(CCTK_VarIndex("Burgers::u"),  CCTK_VarIndex("Burgers::urhs"));
+  
+  /* Register all the evolved Array functions with MoL */
   return;
 }
diff --git a/Examples/Burgers/src/burgers_flux_1.cc b/Examples/Burgers/src/burgers_flux_1.cc
index 827bbd6..947110f 100644
--- a/Examples/Burgers/src/burgers_flux_1.cc
+++ b/Examples/Burgers/src/burgers_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)
@@ -32,43 +34,28 @@ extern "C" void burgers_flux_1_SelectBCs(CCTK_ARGUMENTS)
   return;
 }
 
-static void burgers_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 burgers_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 burgers_flux_1_Body");
-  }
-  
-  if (cctk_iteration % burgers_flux_1_calc_every != burgers_flux_1_calc_offset)
-  {
-    return;
-  }
-  
-  const char *groups[] = {"Burgers::uF_group","Burgers::uLeft_group","Burgers::uR_group"};
-  GenericFD_AssertGroupStorage(cctkGH, "burgers_flux_1", 3, groups);
-  
-  GenericFD_EnsureStencilFits(cctkGH, "burgers_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);
@@ -87,35 +74,41 @@ static void burgers_flux_1_Body(cGH const * restrict const cctkGH, int const dir
   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 (burgers_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 uLeftL = uLeft[index];
-        CCTK_REAL uRL = uR[index];
-        
-        
-        /* Include user supplied include files */
-        
-        /* Precompute derivatives */
-        CCTK_REAL const ShiftMinus1uR = ShiftMinus1(&uR[index]);
-        
-        /* Calculate temporaries and grid functions */
-        CCTK_REAL uFL = 0.25*(SQR(ShiftMinus1uR) + uLeftL*(uLeftL - 
-          2*ToReal(alpha)) + 2*ShiftMinus1uR*ToReal(alpha));
-        
-        /* Copy local copies back to grid functions */
-        uF[index] = uFL;
-      }
-    }
+    ptrdiff_t const index = di*i + dj*j + dk*k;
+    
+    /* Assign local copies of grid functions */
+    
+    CCTK_REAL uLeftL = uLeft[index];
+    CCTK_REAL uRL = uR[index];
+    
+    
+    /* Include user supplied include files */
+    
+    /* Precompute derivatives */
+    CCTK_REAL const ShiftMinus1uR = ShiftMinus1(&uR[index]);
+    
+    /* Calculate temporaries and grid functions */
+    CCTK_REAL uFL = 0.25*(SQR(ShiftMinus1uR) + uLeftL*(uLeftL - 
+      2*ToReal(alpha)) + 2*ShiftMinus1uR*ToReal(alpha));
+    
+    /* Copy local copies back to grid functions */
+    uF[index] = uFL;
   }
+  CCTK_ENDLOOP3 (burgers_flux_1);
 }
 
 extern "C" void burgers_flux_1(CCTK_ARGUMENTS)
@@ -123,5 +116,26 @@ extern "C" void burgers_flux_1(CCTK_ARGUMENTS)
   DECLARE_CCTK_ARGUMENTS;
   DECLARE_CCTK_PARAMETERS;
   
+  
+  if (verbose > 1)
+  {
+    CCTK_VInfo(CCTK_THORNSTRING,"Entering burgers_flux_1_Body");
+  }
+  
+  if (cctk_iteration % burgers_flux_1_calc_every != burgers_flux_1_calc_offset)
+  {
+    return;
+  }
+  
+  const char *groups[] = {"Burgers::uF_group","Burgers::uLeft_group","Burgers::uR_group"};
+  GenericFD_AssertGroupStorage(cctkGH, "burgers_flux_1", 3, groups);
+  
+  GenericFD_EnsureStencilFits(cctkGH, "burgers_flux_1", 1, 1, 1);
+  
   GenericFD_LoopOverInterior(cctkGH, &burgers_flux_1_Body);
+  
+  if (verbose > 1)
+  {
+    CCTK_VInfo(CCTK_THORNSTRING,"Leaving burgers_flux_1_Body");
+  }
 }
diff --git a/Examples/Burgers/src/burgers_initial_shock.cc b/Examples/Burgers/src/burgers_initial_shock.cc
index 3baf23f..6aaa0fc 100644
--- a/Examples/Burgers/src/burgers_initial_shock.cc
+++ b/Examples/Burgers/src/burgers_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 burgers_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 burgers_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 burgers_initial_shock_Body");
-  }
-  
-  if (cctk_iteration % burgers_initial_shock_calc_every != burgers_initial_shock_calc_offset)
-  {
-    return;
-  }
-  
-  const char *groups[] = {"grid::coordinates","Burgers::u_group"};
-  GenericFD_AssertGroupStorage(cctkGH, "burgers_initial_shock", 2, 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,33 +62,39 @@ static void burgers_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 (burgers_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 uL = ToReal(uL0) + StepFunction(-0.5 + xL)*(-ToReal(uL0) + 
-          ToReal(uR0));
-        
-        /* Copy local copies back to grid functions */
-        u[index] = uL;
-      }
-    }
+    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 uL = ToReal(uL0) + StepFunction(-0.5 + xL)*(-ToReal(uL0) + 
+      ToReal(uR0));
+    
+    /* Copy local copies back to grid functions */
+    u[index] = uL;
   }
+  CCTK_ENDLOOP3 (burgers_initial_shock);
 }
 
 extern "C" void burgers_initial_shock(CCTK_ARGUMENTS)
@@ -108,5 +102,25 @@ extern "C" void burgers_initial_shock(CCTK_ARGUMENTS)
   DECLARE_CCTK_ARGUMENTS;
   DECLARE_CCTK_PARAMETERS;
   
+  
+  if (verbose > 1)
+  {
+    CCTK_VInfo(CCTK_THORNSTRING,"Entering burgers_initial_shock_Body");
+  }
+  
+  if (cctk_iteration % burgers_initial_shock_calc_every != burgers_initial_shock_calc_offset)
+  {
+    return;
+  }
+  
+  const char *groups[] = {"grid::coordinates","Burgers::u_group"};
+  GenericFD_AssertGroupStorage(cctkGH, "burgers_initial_shock", 2, groups);
+  
+  
   GenericFD_LoopOverEverything(cctkGH, &burgers_initial_shock_Body);
+  
+  if (verbose > 1)
+  {
+    CCTK_VInfo(CCTK_THORNSTRING,"Leaving burgers_initial_shock_Body");
+  }
 }
diff --git a/Examples/Burgers/src/burgers_initial_sine.cc b/Examples/Burgers/src/burgers_initial_sine.cc
index 2cc003f..d6aa092 100644
--- a/Examples/Burgers/src/burgers_initial_sine.cc
+++ b/Examples/Burgers/src/burgers_initial_sine.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 burgers_initial_sine_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 burgers_initial_sine_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 burgers_initial_sine_Body");
-  }
-  
-  if (cctk_iteration % burgers_initial_sine_calc_every != burgers_initial_sine_calc_offset)
-  {
-    return;
-  }
-  
-  const char *groups[] = {"grid::coordinates","Burgers::u_group"};
-  GenericFD_AssertGroupStorage(cctkGH, "burgers_initial_sine", 2, 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,32 +62,38 @@ static void burgers_initial_sine_Body(cGH const * restrict const cctkGH, int con
   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 (burgers_initial_sine,
+    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 uL = 1 + Sin(2*Pi*xL)*ToReal(amp);
-        
-        /* Copy local copies back to grid functions */
-        u[index] = uL;
-      }
-    }
+    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 uL = 1 + Sin(2*xL*Pi)*ToReal(amp);
+    
+    /* Copy local copies back to grid functions */
+    u[index] = uL;
   }
+  CCTK_ENDLOOP3 (burgers_initial_sine);
 }
 
 extern "C" void burgers_initial_sine(CCTK_ARGUMENTS)
@@ -107,5 +101,25 @@ extern "C" void burgers_initial_sine(CCTK_ARGUMENTS)
   DECLARE_CCTK_ARGUMENTS;
   DECLARE_CCTK_PARAMETERS;
   
+  
+  if (verbose > 1)
+  {
+    CCTK_VInfo(CCTK_THORNSTRING,"Entering burgers_initial_sine_Body");
+  }
+  
+  if (cctk_iteration % burgers_initial_sine_calc_every != burgers_initial_sine_calc_offset)
+  {
+    return;
+  }
+  
+  const char *groups[] = {"grid::coordinates","Burgers::u_group"};
+  GenericFD_AssertGroupStorage(cctkGH, "burgers_initial_sine", 2, groups);
+  
+  
   GenericFD_LoopOverEverything(cctkGH, &burgers_initial_sine_Body);
+  
+  if (verbose > 1)
+  {
+    CCTK_VInfo(CCTK_THORNSTRING,"Leaving burgers_initial_sine_Body");
+  }
 }
diff --git a/Examples/Burgers/src/burgers_reconstruct_1.cc b/Examples/Burgers/src/burgers_reconstruct_1.cc
index 59622d8..4a79c1e 100644
--- a/Examples/Burgers/src/burgers_reconstruct_1.cc
+++ b/Examples/Burgers/src/burgers_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)
@@ -35,43 +37,28 @@ extern "C" void burgers_reconstruct_1_SelectBCs(CCTK_ARGUMENTS)
   return;
 }
 
-static void burgers_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 burgers_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 burgers_reconstruct_1_Body");
-  }
-  
-  if (cctk_iteration % burgers_reconstruct_1_calc_every != burgers_reconstruct_1_calc_offset)
-  {
-    return;
-  }
-  
-  const char *groups[] = {"Burgers::u_group","Burgers::uLeft_group","Burgers::uR_group"};
-  GenericFD_AssertGroupStorage(cctkGH, "burgers_reconstruct_1", 3, groups);
-  
-  GenericFD_EnsureStencilFits(cctkGH, "burgers_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);
@@ -90,44 +77,50 @@ static void burgers_reconstruct_1_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 (burgers_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 uL = u[index];
-        
-        
-        /* Include user supplied include files */
-        
-        /* Precompute derivatives */
-        CCTK_REAL const DiffPlus1u = DiffPlus1(&u[index]);
-        CCTK_REAL const DiffMinus1u = DiffMinus1(&u[index]);
-        
-        /* Calculate temporaries and grid functions */
-        CCTK_REAL slopeL = DiffMinus1u;
-        
-        CCTK_REAL slopeR = DiffPlus1u;
-        
-        CCTK_REAL slope = IfThen(slopeL*slopeR < 0,0,IfThen(Abs(slopeL) < 
-          Abs(slopeR),slopeL,slopeR));
-        
-        CCTK_REAL uLeftL = -0.5*slope + uL;
-        
-        CCTK_REAL uRL = 0.5*slope + uL;
-        
-        /* Copy local copies back to grid functions */
-        uLeft[index] = uLeftL;
-        uR[index] = uRL;
-      }
-    }
+    ptrdiff_t const index = di*i + dj*j + dk*k;
+    
+    /* Assign local copies of grid functions */
+    
+    CCTK_REAL uL = u[index];
+    
+    
+    /* Include user supplied include files */
+    
+    /* Precompute derivatives */
+    CCTK_REAL const DiffPlus1u = DiffPlus1(&u[index]);
+    CCTK_REAL const DiffMinus1u = DiffMinus1(&u[index]);
+    
+    /* Calculate temporaries and grid functions */
+    CCTK_REAL slopeL = DiffMinus1u;
+    
+    CCTK_REAL slopeR = DiffPlus1u;
+    
+    CCTK_REAL slope = IfThen(slopeL*slopeR < 0,0,IfThen(Abs(slopeL) < 
+      Abs(slopeR),slopeL,slopeR));
+    
+    CCTK_REAL uLeftL = uL - 0.5*slope;
+    
+    CCTK_REAL uRL = uL + 0.5*slope;
+    
+    /* Copy local copies back to grid functions */
+    uLeft[index] = uLeftL;
+    uR[index] = uRL;
   }
+  CCTK_ENDLOOP3 (burgers_reconstruct_1);
 }
 
 extern "C" void burgers_reconstruct_1(CCTK_ARGUMENTS)
@@ -135,5 +128,26 @@ extern "C" void burgers_reconstruct_1(CCTK_ARGUMENTS)
   DECLARE_CCTK_ARGUMENTS;
   DECLARE_CCTK_PARAMETERS;
   
+  
+  if (verbose > 1)
+  {
+    CCTK_VInfo(CCTK_THORNSTRING,"Entering burgers_reconstruct_1_Body");
+  }
+  
+  if (cctk_iteration % burgers_reconstruct_1_calc_every != burgers_reconstruct_1_calc_offset)
+  {
+    return;
+  }
+  
+  const char *groups[] = {"Burgers::u_group","Burgers::uLeft_group","Burgers::uR_group"};
+  GenericFD_AssertGroupStorage(cctkGH, "burgers_reconstruct_1", 3, groups);
+  
+  GenericFD_EnsureStencilFits(cctkGH, "burgers_reconstruct_1", 1, 1, 1);
+  
   GenericFD_LoopOverInterior(cctkGH, &burgers_reconstruct_1_Body);
+  
+  if (verbose > 1)
+  {
+    CCTK_VInfo(CCTK_THORNSTRING,"Leaving burgers_reconstruct_1_Body");
+  }
 }
diff --git a/Examples/Burgers/src/burgers_rhs_1.cc b/Examples/Burgers/src/burgers_rhs_1.cc
index 375ad35..557e7a0 100644
--- a/Examples/Burgers/src/burgers_rhs_1.cc
+++ b/Examples/Burgers/src/burgers_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)
@@ -32,43 +34,28 @@ extern "C" void burgers_rhs_1_SelectBCs(CCTK_ARGUMENTS)
   return;
 }
 
-static void burgers_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 burgers_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 burgers_rhs_1_Body");
-  }
-  
-  if (cctk_iteration % burgers_rhs_1_calc_every != burgers_rhs_1_calc_offset)
-  {
-    return;
-  }
-  
-  const char *groups[] = {"Burgers::uF_group","Burgers::u_grouprhs"};
-  GenericFD_AssertGroupStorage(cctkGH, "burgers_rhs_1", 2, groups);
-  
-  GenericFD_EnsureStencilFits(cctkGH, "burgers_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);
@@ -87,34 +74,40 @@ static void burgers_rhs_1_Body(cGH const * restrict const cctkGH, int const dir,
   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 (burgers_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 uFL = uF[index];
-        CCTK_REAL urhsL = urhs[index];
-        
-        
-        /* Include user supplied include files */
-        
-        /* Precompute derivatives */
-        CCTK_REAL const PDplus1uF = PDplus1(&uF[index]);
-        
-        /* Calculate temporaries and grid functions */
-        urhsL = -PDplus1uF + urhsL;
-        
-        /* Copy local copies back to grid functions */
-        urhs[index] = urhsL;
-      }
-    }
+    ptrdiff_t const index = di*i + dj*j + dk*k;
+    
+    /* Assign local copies of grid functions */
+    
+    CCTK_REAL uFL = uF[index];
+    CCTK_REAL urhsL = urhs[index];
+    
+    
+    /* Include user supplied include files */
+    
+    /* Precompute derivatives */
+    CCTK_REAL const PDplus1uF = PDplus1(&uF[index]);
+    
+    /* Calculate temporaries and grid functions */
+    urhsL = urhsL - PDplus1uF;
+    
+    /* Copy local copies back to grid functions */
+    urhs[index] = urhsL;
   }
+  CCTK_ENDLOOP3 (burgers_rhs_1);
 }
 
 extern "C" void burgers_rhs_1(CCTK_ARGUMENTS)
@@ -122,5 +115,26 @@ extern "C" void burgers_rhs_1(CCTK_ARGUMENTS)
   DECLARE_CCTK_ARGUMENTS;
   DECLARE_CCTK_PARAMETERS;
   
+  
+  if (verbose > 1)
+  {
+    CCTK_VInfo(CCTK_THORNSTRING,"Entering burgers_rhs_1_Body");
+  }
+  
+  if (cctk_iteration % burgers_rhs_1_calc_every != burgers_rhs_1_calc_offset)
+  {
+    return;
+  }
+  
+  const char *groups[] = {"Burgers::uF_group","Burgers::u_grouprhs"};
+  GenericFD_AssertGroupStorage(cctkGH, "burgers_rhs_1", 2, groups);
+  
+  GenericFD_EnsureStencilFits(cctkGH, "burgers_rhs_1", 1, 1, 1);
+  
   GenericFD_LoopOverInterior(cctkGH, &burgers_rhs_1_Body);
+  
+  if (verbose > 1)
+  {
+    CCTK_VInfo(CCTK_THORNSTRING,"Leaving burgers_rhs_1_Body");
+  }
 }
diff --git a/Examples/Burgers/src/burgers_zero_rhs.cc b/Examples/Burgers/src/burgers_zero_rhs.cc
index 538d321..10f020f 100644
--- a/Examples/Burgers/src/burgers_zero_rhs.cc
+++ b/Examples/Burgers/src/burgers_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 burgers_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 burgers_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 burgers_zero_rhs_Body");
-  }
-  
-  if (cctk_iteration % burgers_zero_rhs_calc_every != burgers_zero_rhs_calc_offset)
-  {
-    return;
-  }
-  
-  const char *groups[] = {"Burgers::u_grouprhs"};
-  GenericFD_AssertGroupStorage(cctkGH, "burgers_zero_rhs", 1, 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,31 +62,37 @@ static void burgers_zero_rhs_Body(cGH const * restrict const cctkGH, int const d
   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 (burgers_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 urhsL = 0;
-        
-        /* Copy local copies back to grid functions */
-        urhs[index] = urhsL;
-      }
-    }
+    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 urhsL = 0;
+    
+    /* Copy local copies back to grid functions */
+    urhs[index] = urhsL;
   }
+  CCTK_ENDLOOP3 (burgers_zero_rhs);
 }
 
 extern "C" void burgers_zero_rhs(CCTK_ARGUMENTS)
@@ -106,5 +100,25 @@ extern "C" void burgers_zero_rhs(CCTK_ARGUMENTS)
   DECLARE_CCTK_ARGUMENTS;
   DECLARE_CCTK_PARAMETERS;
   
+  
+  if (verbose > 1)
+  {
+    CCTK_VInfo(CCTK_THORNSTRING,"Entering burgers_zero_rhs_Body");
+  }
+  
+  if (cctk_iteration % burgers_zero_rhs_calc_every != burgers_zero_rhs_calc_offset)
+  {
+    return;
+  }
+  
+  const char *groups[] = {"Burgers::u_grouprhs"};
+  GenericFD_AssertGroupStorage(cctkGH, "burgers_zero_rhs", 1, groups);
+  
+  
   GenericFD_LoopOverEverything(cctkGH, &burgers_zero_rhs_Body);
+  
+  if (verbose > 1)
+  {
+    CCTK_VInfo(CCTK_THORNSTRING,"Leaving burgers_zero_rhs_Body");
+  }
 }