Update auto-generated code as generated by the current version of Kranc

11 files changed, 2913 insertions, 0 deletions

diff --git a/ML_ADM/src/Boundaries.cc b/ML_ADM/src/Boundaries.cc
new file mode 100644
index 0000000..4e20959
--- /dev/null
+++ b/ML_ADM/src/Boundaries.cc
@@ -0,0 +1,1081 @@
+/*  File produced by Kranc */
+
+#include "cctk.h"
+#include "cctk_Arguments.h"
+#include "cctk_Parameters.h"
+#include "cctk_Faces.h"
+#include "util_Table.h"
+#include "Symmetry.h"
+
+
+/* the boundary treatment is split into 3 steps:    */
+/* 1. excision                                      */
+/* 2. symmetries                                    */
+/* 3. "other" boundary conditions, e.g. radiative */
+
+/* to simplify scheduling and testing, the 3 steps  */
+/* are currently applied in separate functions      */
+
+
+extern "C" void ML_ADM_CheckBoundaries(CCTK_ARGUMENTS)
+{
+  DECLARE_CCTK_ARGUMENTS;
+  DECLARE_CCTK_PARAMETERS;
+  
+  return;
+}
+
+extern "C" void ML_ADM_SelectBoundConds(CCTK_ARGUMENTS)
+{
+  DECLARE_CCTK_ARGUMENTS;
+  DECLARE_CCTK_PARAMETERS;
+  
+  CCTK_INT ierr = 0;
+  
+  if (CCTK_EQUALS(ML_curv_bound, "none"  ) ||
+      CCTK_EQUALS(ML_curv_bound, "static") ||
+      CCTK_EQUALS(ML_curv_bound, "flat"  ) ||
+      CCTK_EQUALS(ML_curv_bound, "zero"  ) )
+  {
+    ierr = Boundary_SelectGroupForBC(cctkGH, CCTK_ALL_FACES, 1, -1,
+                      "ML_ADM::ML_curv", ML_curv_bound);
+    if (ierr < 0)
+       CCTK_WARN(0, "Failed to register ML_curv_bound BC for ML_ADM::ML_curv!");
+  }
+  
+  if (CCTK_EQUALS(ML_lapse_bound, "none"  ) ||
+      CCTK_EQUALS(ML_lapse_bound, "static") ||
+      CCTK_EQUALS(ML_lapse_bound, "flat"  ) ||
+      CCTK_EQUALS(ML_lapse_bound, "zero"  ) )
+  {
+    ierr = Boundary_SelectGroupForBC(cctkGH, CCTK_ALL_FACES, 1, -1,
+                      "ML_ADM::ML_lapse", ML_lapse_bound);
+    if (ierr < 0)
+       CCTK_WARN(0, "Failed to register ML_lapse_bound BC for ML_ADM::ML_lapse!");
+  }
+  
+  if (CCTK_EQUALS(ML_metric_bound, "none"  ) ||
+      CCTK_EQUALS(ML_metric_bound, "static") ||
+      CCTK_EQUALS(ML_metric_bound, "flat"  ) ||
+      CCTK_EQUALS(ML_metric_bound, "zero"  ) )
+  {
+    ierr = Boundary_SelectGroupForBC(cctkGH, CCTK_ALL_FACES, 1, -1,
+                      "ML_ADM::ML_metric", ML_metric_bound);
+    if (ierr < 0)
+       CCTK_WARN(0, "Failed to register ML_metric_bound BC for ML_ADM::ML_metric!");
+  }
+  
+  if (CCTK_EQUALS(ML_shift_bound, "none"  ) ||
+      CCTK_EQUALS(ML_shift_bound, "static") ||
+      CCTK_EQUALS(ML_shift_bound, "flat"  ) ||
+      CCTK_EQUALS(ML_shift_bound, "zero"  ) )
+  {
+    ierr = Boundary_SelectGroupForBC(cctkGH, CCTK_ALL_FACES, 1, -1,
+                      "ML_ADM::ML_shift", ML_shift_bound);
+    if (ierr < 0)
+       CCTK_WARN(0, "Failed to register ML_shift_bound BC for ML_ADM::ML_shift!");
+  }
+  
+  if (CCTK_EQUALS(K11_bound, "none"  ) ||
+      CCTK_EQUALS(K11_bound, "static") ||
+      CCTK_EQUALS(K11_bound, "flat"  ) ||
+      CCTK_EQUALS(K11_bound, "zero"  ) )
+  {
+    ierr = Boundary_SelectVarForBC(cctkGH, CCTK_ALL_FACES, 1, -1,
+                      "ML_ADM::K11", K11_bound);
+    if (ierr < 0)
+       CCTK_WARN(0, "Failed to register K11_bound BC for ML_ADM::K11!");
+  }
+  
+  if (CCTK_EQUALS(K12_bound, "none"  ) ||
+      CCTK_EQUALS(K12_bound, "static") ||
+      CCTK_EQUALS(K12_bound, "flat"  ) ||
+      CCTK_EQUALS(K12_bound, "zero"  ) )
+  {
+    ierr = Boundary_SelectVarForBC(cctkGH, CCTK_ALL_FACES, 1, -1,
+                      "ML_ADM::K12", K12_bound);
+    if (ierr < 0)
+       CCTK_WARN(0, "Failed to register K12_bound BC for ML_ADM::K12!");
+  }
+  
+  if (CCTK_EQUALS(K13_bound, "none"  ) ||
+      CCTK_EQUALS(K13_bound, "static") ||
+      CCTK_EQUALS(K13_bound, "flat"  ) ||
+      CCTK_EQUALS(K13_bound, "zero"  ) )
+  {
+    ierr = Boundary_SelectVarForBC(cctkGH, CCTK_ALL_FACES, 1, -1,
+                      "ML_ADM::K13", K13_bound);
+    if (ierr < 0)
+       CCTK_WARN(0, "Failed to register K13_bound BC for ML_ADM::K13!");
+  }
+  
+  if (CCTK_EQUALS(K22_bound, "none"  ) ||
+      CCTK_EQUALS(K22_bound, "static") ||
+      CCTK_EQUALS(K22_bound, "flat"  ) ||
+      CCTK_EQUALS(K22_bound, "zero"  ) )
+  {
+    ierr = Boundary_SelectVarForBC(cctkGH, CCTK_ALL_FACES, 1, -1,
+                      "ML_ADM::K22", K22_bound);
+    if (ierr < 0)
+       CCTK_WARN(0, "Failed to register K22_bound BC for ML_ADM::K22!");
+  }
+  
+  if (CCTK_EQUALS(K23_bound, "none"  ) ||
+      CCTK_EQUALS(K23_bound, "static") ||
+      CCTK_EQUALS(K23_bound, "flat"  ) ||
+      CCTK_EQUALS(K23_bound, "zero"  ) )
+  {
+    ierr = Boundary_SelectVarForBC(cctkGH, CCTK_ALL_FACES, 1, -1,
+                      "ML_ADM::K23", K23_bound);
+    if (ierr < 0)
+       CCTK_WARN(0, "Failed to register K23_bound BC for ML_ADM::K23!");
+  }
+  
+  if (CCTK_EQUALS(K33_bound, "none"  ) ||
+      CCTK_EQUALS(K33_bound, "static") ||
+      CCTK_EQUALS(K33_bound, "flat"  ) ||
+      CCTK_EQUALS(K33_bound, "zero"  ) )
+  {
+    ierr = Boundary_SelectVarForBC(cctkGH, CCTK_ALL_FACES, 1, -1,
+                      "ML_ADM::K33", K33_bound);
+    if (ierr < 0)
+       CCTK_WARN(0, "Failed to register K33_bound BC for ML_ADM::K33!");
+  }
+  
+  if (CCTK_EQUALS(alpha_bound, "none"  ) ||
+      CCTK_EQUALS(alpha_bound, "static") ||
+      CCTK_EQUALS(alpha_bound, "flat"  ) ||
+      CCTK_EQUALS(alpha_bound, "zero"  ) )
+  {
+    ierr = Boundary_SelectVarForBC(cctkGH, CCTK_ALL_FACES, 1, -1,
+                      "ML_ADM::alpha", alpha_bound);
+    if (ierr < 0)
+       CCTK_WARN(0, "Failed to register alpha_bound BC for ML_ADM::alpha!");
+  }
+  
+  if (CCTK_EQUALS(g11_bound, "none"  ) ||
+      CCTK_EQUALS(g11_bound, "static") ||
+      CCTK_EQUALS(g11_bound, "flat"  ) ||
+      CCTK_EQUALS(g11_bound, "zero"  ) )
+  {
+    ierr = Boundary_SelectVarForBC(cctkGH, CCTK_ALL_FACES, 1, -1,
+                      "ML_ADM::g11", g11_bound);
+    if (ierr < 0)
+       CCTK_WARN(0, "Failed to register g11_bound BC for ML_ADM::g11!");
+  }
+  
+  if (CCTK_EQUALS(g12_bound, "none"  ) ||
+      CCTK_EQUALS(g12_bound, "static") ||
+      CCTK_EQUALS(g12_bound, "flat"  ) ||
+      CCTK_EQUALS(g12_bound, "zero"  ) )
+  {
+    ierr = Boundary_SelectVarForBC(cctkGH, CCTK_ALL_FACES, 1, -1,
+                      "ML_ADM::g12", g12_bound);
+    if (ierr < 0)
+       CCTK_WARN(0, "Failed to register g12_bound BC for ML_ADM::g12!");
+  }
+  
+  if (CCTK_EQUALS(g13_bound, "none"  ) ||
+      CCTK_EQUALS(g13_bound, "static") ||
+      CCTK_EQUALS(g13_bound, "flat"  ) ||
+      CCTK_EQUALS(g13_bound, "zero"  ) )
+  {
+    ierr = Boundary_SelectVarForBC(cctkGH, CCTK_ALL_FACES, 1, -1,
+                      "ML_ADM::g13", g13_bound);
+    if (ierr < 0)
+       CCTK_WARN(0, "Failed to register g13_bound BC for ML_ADM::g13!");
+  }
+  
+  if (CCTK_EQUALS(g22_bound, "none"  ) ||
+      CCTK_EQUALS(g22_bound, "static") ||
+      CCTK_EQUALS(g22_bound, "flat"  ) ||
+      CCTK_EQUALS(g22_bound, "zero"  ) )
+  {
+    ierr = Boundary_SelectVarForBC(cctkGH, CCTK_ALL_FACES, 1, -1,
+                      "ML_ADM::g22", g22_bound);
+    if (ierr < 0)
+       CCTK_WARN(0, "Failed to register g22_bound BC for ML_ADM::g22!");
+  }
+  
+  if (CCTK_EQUALS(g23_bound, "none"  ) ||
+      CCTK_EQUALS(g23_bound, "static") ||
+      CCTK_EQUALS(g23_bound, "flat"  ) ||
+      CCTK_EQUALS(g23_bound, "zero"  ) )
+  {
+    ierr = Boundary_SelectVarForBC(cctkGH, CCTK_ALL_FACES, 1, -1,
+                      "ML_ADM::g23", g23_bound);
+    if (ierr < 0)
+       CCTK_WARN(0, "Failed to register g23_bound BC for ML_ADM::g23!");
+  }
+  
+  if (CCTK_EQUALS(g33_bound, "none"  ) ||
+      CCTK_EQUALS(g33_bound, "static") ||
+      CCTK_EQUALS(g33_bound, "flat"  ) ||
+      CCTK_EQUALS(g33_bound, "zero"  ) )
+  {
+    ierr = Boundary_SelectVarForBC(cctkGH, CCTK_ALL_FACES, 1, -1,
+                      "ML_ADM::g33", g33_bound);
+    if (ierr < 0)
+       CCTK_WARN(0, "Failed to register g33_bound BC for ML_ADM::g33!");
+  }
+  
+  if (CCTK_EQUALS(beta1_bound, "none"  ) ||
+      CCTK_EQUALS(beta1_bound, "static") ||
+      CCTK_EQUALS(beta1_bound, "flat"  ) ||
+      CCTK_EQUALS(beta1_bound, "zero"  ) )
+  {
+    ierr = Boundary_SelectVarForBC(cctkGH, CCTK_ALL_FACES, 1, -1,
+                      "ML_ADM::beta1", beta1_bound);
+    if (ierr < 0)
+       CCTK_WARN(0, "Failed to register beta1_bound BC for ML_ADM::beta1!");
+  }
+  
+  if (CCTK_EQUALS(beta2_bound, "none"  ) ||
+      CCTK_EQUALS(beta2_bound, "static") ||
+      CCTK_EQUALS(beta2_bound, "flat"  ) ||
+      CCTK_EQUALS(beta2_bound, "zero"  ) )
+  {
+    ierr = Boundary_SelectVarForBC(cctkGH, CCTK_ALL_FACES, 1, -1,
+                      "ML_ADM::beta2", beta2_bound);
+    if (ierr < 0)
+       CCTK_WARN(0, "Failed to register beta2_bound BC for ML_ADM::beta2!");
+  }
+  
+  if (CCTK_EQUALS(beta3_bound, "none"  ) ||
+      CCTK_EQUALS(beta3_bound, "static") ||
+      CCTK_EQUALS(beta3_bound, "flat"  ) ||
+      CCTK_EQUALS(beta3_bound, "zero"  ) )
+  {
+    ierr = Boundary_SelectVarForBC(cctkGH, CCTK_ALL_FACES, 1, -1,
+                      "ML_ADM::beta3", beta3_bound);
+    if (ierr < 0)
+       CCTK_WARN(0, "Failed to register beta3_bound BC for ML_ADM::beta3!");
+  }
+  
+  if (CCTK_EQUALS(ML_curv_bound, "radiative"))
+  {
+   /* select radiation boundary condition */
+    static CCTK_INT handle_ML_curv_bound = -1;
+    if (handle_ML_curv_bound < 0) handle_ML_curv_bound = Util_TableCreate(UTIL_TABLE_FLAGS_CASE_INSENSITIVE);
+    if (handle_ML_curv_bound < 0) CCTK_WARN(0, "could not create table!");
+    if (Util_TableSetReal(handle_ML_curv_bound , ML_curv_bound_limit, "LIMIT") < 0)
+       CCTK_WARN(0, "could not set LIMIT value in table!");
+    if (Util_TableSetReal(handle_ML_curv_bound ,ML_curv_bound_speed, "SPEED") < 0)
+       CCTK_WARN(0, "could not set SPEED value in table!");
+  
+    ierr = Boundary_SelectGroupForBC(cctkGH, CCTK_ALL_FACES, 1, handle_ML_curv_bound, 
+                      "ML_ADM::ML_curv", "Radiation");
+  
+    if (ierr < 0)
+       CCTK_WARN(0, "Failed to register Radiation BC for ML_ADM::ML_curv!");
+  
+  }
+  
+  if (CCTK_EQUALS(ML_lapse_bound, "radiative"))
+  {
+   /* select radiation boundary condition */
+    static CCTK_INT handle_ML_lapse_bound = -1;
+    if (handle_ML_lapse_bound < 0) handle_ML_lapse_bound = Util_TableCreate(UTIL_TABLE_FLAGS_CASE_INSENSITIVE);
+    if (handle_ML_lapse_bound < 0) CCTK_WARN(0, "could not create table!");
+    if (Util_TableSetReal(handle_ML_lapse_bound , ML_lapse_bound_limit, "LIMIT") < 0)
+       CCTK_WARN(0, "could not set LIMIT value in table!");
+    if (Util_TableSetReal(handle_ML_lapse_bound ,ML_lapse_bound_speed, "SPEED") < 0)
+       CCTK_WARN(0, "could not set SPEED value in table!");
+  
+    ierr = Boundary_SelectGroupForBC(cctkGH, CCTK_ALL_FACES, 1, handle_ML_lapse_bound, 
+                      "ML_ADM::ML_lapse", "Radiation");
+  
+    if (ierr < 0)
+       CCTK_WARN(0, "Failed to register Radiation BC for ML_ADM::ML_lapse!");
+  
+  }
+  
+  if (CCTK_EQUALS(ML_metric_bound, "radiative"))
+  {
+   /* select radiation boundary condition */
+    static CCTK_INT handle_ML_metric_bound = -1;
+    if (handle_ML_metric_bound < 0) handle_ML_metric_bound = Util_TableCreate(UTIL_TABLE_FLAGS_CASE_INSENSITIVE);
+    if (handle_ML_metric_bound < 0) CCTK_WARN(0, "could not create table!");
+    if (Util_TableSetReal(handle_ML_metric_bound , ML_metric_bound_limit, "LIMIT") < 0)
+       CCTK_WARN(0, "could not set LIMIT value in table!");
+    if (Util_TableSetReal(handle_ML_metric_bound ,ML_metric_bound_speed, "SPEED") < 0)
+       CCTK_WARN(0, "could not set SPEED value in table!");
+  
+    ierr = Boundary_SelectGroupForBC(cctkGH, CCTK_ALL_FACES, 1, handle_ML_metric_bound, 
+                      "ML_ADM::ML_metric", "Radiation");
+  
+    if (ierr < 0)
+       CCTK_WARN(0, "Failed to register Radiation BC for ML_ADM::ML_metric!");
+  
+  }
+  
+  if (CCTK_EQUALS(ML_shift_bound, "radiative"))
+  {
+   /* select radiation boundary condition */
+    static CCTK_INT handle_ML_shift_bound = -1;
+    if (handle_ML_shift_bound < 0) handle_ML_shift_bound = Util_TableCreate(UTIL_TABLE_FLAGS_CASE_INSENSITIVE);
+    if (handle_ML_shift_bound < 0) CCTK_WARN(0, "could not create table!");
+    if (Util_TableSetReal(handle_ML_shift_bound , ML_shift_bound_limit, "LIMIT") < 0)
+       CCTK_WARN(0, "could not set LIMIT value in table!");
+    if (Util_TableSetReal(handle_ML_shift_bound ,ML_shift_bound_speed, "SPEED") < 0)
+       CCTK_WARN(0, "could not set SPEED value in table!");
+  
+    ierr = Boundary_SelectGroupForBC(cctkGH, CCTK_ALL_FACES, 1, handle_ML_shift_bound, 
+                      "ML_ADM::ML_shift", "Radiation");
+  
+    if (ierr < 0)
+       CCTK_WARN(0, "Failed to register Radiation BC for ML_ADM::ML_shift!");
+  
+  }
+  
+  if (CCTK_EQUALS(K11_bound, "radiative"))
+  {
+   /* select radiation boundary condition */
+    static CCTK_INT handle_K11_bound = -1;
+    if (handle_K11_bound < 0) handle_K11_bound = Util_TableCreate(UTIL_TABLE_FLAGS_CASE_INSENSITIVE);
+    if (handle_K11_bound < 0) CCTK_WARN(0, "could not create table!");
+    if (Util_TableSetReal(handle_K11_bound , K11_bound_limit, "LIMIT") < 0)
+       CCTK_WARN(0, "could not set LIMIT value in table!");
+    if (Util_TableSetReal(handle_K11_bound ,K11_bound_speed, "SPEED") < 0)
+        CCTK_WARN(0, "could not set SPEED value in table!");
+  
+    ierr = Boundary_SelectVarForBC(cctkGH, CCTK_ALL_FACES, 1, handle_K11_bound, 
+                      "ML_ADM::K11", "Radiation");
+  
+    if (ierr < 0)
+       CCTK_WARN(0, "Failed to register Radiation BC for ML_ADM::K11!");
+  
+  }
+  
+  if (CCTK_EQUALS(K12_bound, "radiative"))
+  {
+   /* select radiation boundary condition */
+    static CCTK_INT handle_K12_bound = -1;
+    if (handle_K12_bound < 0) handle_K12_bound = Util_TableCreate(UTIL_TABLE_FLAGS_CASE_INSENSITIVE);
+    if (handle_K12_bound < 0) CCTK_WARN(0, "could not create table!");
+    if (Util_TableSetReal(handle_K12_bound , K12_bound_limit, "LIMIT") < 0)
+       CCTK_WARN(0, "could not set LIMIT value in table!");
+    if (Util_TableSetReal(handle_K12_bound ,K12_bound_speed, "SPEED") < 0)
+        CCTK_WARN(0, "could not set SPEED value in table!");
+  
+    ierr = Boundary_SelectVarForBC(cctkGH, CCTK_ALL_FACES, 1, handle_K12_bound, 
+                      "ML_ADM::K12", "Radiation");
+  
+    if (ierr < 0)
+       CCTK_WARN(0, "Failed to register Radiation BC for ML_ADM::K12!");
+  
+  }
+  
+  if (CCTK_EQUALS(K13_bound, "radiative"))
+  {
+   /* select radiation boundary condition */
+    static CCTK_INT handle_K13_bound = -1;
+    if (handle_K13_bound < 0) handle_K13_bound = Util_TableCreate(UTIL_TABLE_FLAGS_CASE_INSENSITIVE);
+    if (handle_K13_bound < 0) CCTK_WARN(0, "could not create table!");
+    if (Util_TableSetReal(handle_K13_bound , K13_bound_limit, "LIMIT") < 0)
+       CCTK_WARN(0, "could not set LIMIT value in table!");
+    if (Util_TableSetReal(handle_K13_bound ,K13_bound_speed, "SPEED") < 0)
+        CCTK_WARN(0, "could not set SPEED value in table!");
+  
+    ierr = Boundary_SelectVarForBC(cctkGH, CCTK_ALL_FACES, 1, handle_K13_bound, 
+                      "ML_ADM::K13", "Radiation");
+  
+    if (ierr < 0)
+       CCTK_WARN(0, "Failed to register Radiation BC for ML_ADM::K13!");
+  
+  }
+  
+  if (CCTK_EQUALS(K22_bound, "radiative"))
+  {
+   /* select radiation boundary condition */
+    static CCTK_INT handle_K22_bound = -1;
+    if (handle_K22_bound < 0) handle_K22_bound = Util_TableCreate(UTIL_TABLE_FLAGS_CASE_INSENSITIVE);
+    if (handle_K22_bound < 0) CCTK_WARN(0, "could not create table!");
+    if (Util_TableSetReal(handle_K22_bound , K22_bound_limit, "LIMIT") < 0)
+       CCTK_WARN(0, "could not set LIMIT value in table!");
+    if (Util_TableSetReal(handle_K22_bound ,K22_bound_speed, "SPEED") < 0)
+        CCTK_WARN(0, "could not set SPEED value in table!");
+  
+    ierr = Boundary_SelectVarForBC(cctkGH, CCTK_ALL_FACES, 1, handle_K22_bound, 
+                      "ML_ADM::K22", "Radiation");
+  
+    if (ierr < 0)
+       CCTK_WARN(0, "Failed to register Radiation BC for ML_ADM::K22!");
+  
+  }
+  
+  if (CCTK_EQUALS(K23_bound, "radiative"))
+  {
+   /* select radiation boundary condition */
+    static CCTK_INT handle_K23_bound = -1;
+    if (handle_K23_bound < 0) handle_K23_bound = Util_TableCreate(UTIL_TABLE_FLAGS_CASE_INSENSITIVE);
+    if (handle_K23_bound < 0) CCTK_WARN(0, "could not create table!");
+    if (Util_TableSetReal(handle_K23_bound , K23_bound_limit, "LIMIT") < 0)
+       CCTK_WARN(0, "could not set LIMIT value in table!");
+    if (Util_TableSetReal(handle_K23_bound ,K23_bound_speed, "SPEED") < 0)
+        CCTK_WARN(0, "could not set SPEED value in table!");
+  
+    ierr = Boundary_SelectVarForBC(cctkGH, CCTK_ALL_FACES, 1, handle_K23_bound, 
+                      "ML_ADM::K23", "Radiation");
+  
+    if (ierr < 0)
+       CCTK_WARN(0, "Failed to register Radiation BC for ML_ADM::K23!");
+  
+  }
+  
+  if (CCTK_EQUALS(K33_bound, "radiative"))
+  {
+   /* select radiation boundary condition */
+    static CCTK_INT handle_K33_bound = -1;
+    if (handle_K33_bound < 0) handle_K33_bound = Util_TableCreate(UTIL_TABLE_FLAGS_CASE_INSENSITIVE);
+    if (handle_K33_bound < 0) CCTK_WARN(0, "could not create table!");
+    if (Util_TableSetReal(handle_K33_bound , K33_bound_limit, "LIMIT") < 0)
+       CCTK_WARN(0, "could not set LIMIT value in table!");
+    if (Util_TableSetReal(handle_K33_bound ,K33_bound_speed, "SPEED") < 0)
+        CCTK_WARN(0, "could not set SPEED value in table!");
+  
+    ierr = Boundary_SelectVarForBC(cctkGH, CCTK_ALL_FACES, 1, handle_K33_bound, 
+                      "ML_ADM::K33", "Radiation");
+  
+    if (ierr < 0)
+       CCTK_WARN(0, "Failed to register Radiation BC for ML_ADM::K33!");
+  
+  }
+  
+  if (CCTK_EQUALS(alpha_bound, "radiative"))
+  {
+   /* select radiation boundary condition */
+    static CCTK_INT handle_alpha_bound = -1;
+    if (handle_alpha_bound < 0) handle_alpha_bound = Util_TableCreate(UTIL_TABLE_FLAGS_CASE_INSENSITIVE);
+    if (handle_alpha_bound < 0) CCTK_WARN(0, "could not create table!");
+    if (Util_TableSetReal(handle_alpha_bound , alpha_bound_limit, "LIMIT") < 0)
+       CCTK_WARN(0, "could not set LIMIT value in table!");
+    if (Util_TableSetReal(handle_alpha_bound ,alpha_bound_speed, "SPEED") < 0)
+        CCTK_WARN(0, "could not set SPEED value in table!");
+  
+    ierr = Boundary_SelectVarForBC(cctkGH, CCTK_ALL_FACES, 1, handle_alpha_bound, 
+                      "ML_ADM::alpha", "Radiation");
+  
+    if (ierr < 0)
+       CCTK_WARN(0, "Failed to register Radiation BC for ML_ADM::alpha!");
+  
+  }
+  
+  if (CCTK_EQUALS(g11_bound, "radiative"))
+  {
+   /* select radiation boundary condition */
+    static CCTK_INT handle_g11_bound = -1;
+    if (handle_g11_bound < 0) handle_g11_bound = Util_TableCreate(UTIL_TABLE_FLAGS_CASE_INSENSITIVE);
+    if (handle_g11_bound < 0) CCTK_WARN(0, "could not create table!");
+    if (Util_TableSetReal(handle_g11_bound , g11_bound_limit, "LIMIT") < 0)
+       CCTK_WARN(0, "could not set LIMIT value in table!");
+    if (Util_TableSetReal(handle_g11_bound ,g11_bound_speed, "SPEED") < 0)
+        CCTK_WARN(0, "could not set SPEED value in table!");
+  
+    ierr = Boundary_SelectVarForBC(cctkGH, CCTK_ALL_FACES, 1, handle_g11_bound, 
+                      "ML_ADM::g11", "Radiation");
+  
+    if (ierr < 0)
+       CCTK_WARN(0, "Failed to register Radiation BC for ML_ADM::g11!");
+  
+  }
+  
+  if (CCTK_EQUALS(g12_bound, "radiative"))
+  {
+   /* select radiation boundary condition */
+    static CCTK_INT handle_g12_bound = -1;
+    if (handle_g12_bound < 0) handle_g12_bound = Util_TableCreate(UTIL_TABLE_FLAGS_CASE_INSENSITIVE);
+    if (handle_g12_bound < 0) CCTK_WARN(0, "could not create table!");
+    if (Util_TableSetReal(handle_g12_bound , g12_bound_limit, "LIMIT") < 0)
+       CCTK_WARN(0, "could not set LIMIT value in table!");
+    if (Util_TableSetReal(handle_g12_bound ,g12_bound_speed, "SPEED") < 0)
+        CCTK_WARN(0, "could not set SPEED value in table!");
+  
+    ierr = Boundary_SelectVarForBC(cctkGH, CCTK_ALL_FACES, 1, handle_g12_bound, 
+                      "ML_ADM::g12", "Radiation");
+  
+    if (ierr < 0)
+       CCTK_WARN(0, "Failed to register Radiation BC for ML_ADM::g12!");
+  
+  }
+  
+  if (CCTK_EQUALS(g13_bound, "radiative"))
+  {
+   /* select radiation boundary condition */
+    static CCTK_INT handle_g13_bound = -1;
+    if (handle_g13_bound < 0) handle_g13_bound = Util_TableCreate(UTIL_TABLE_FLAGS_CASE_INSENSITIVE);
+    if (handle_g13_bound < 0) CCTK_WARN(0, "could not create table!");
+    if (Util_TableSetReal(handle_g13_bound , g13_bound_limit, "LIMIT") < 0)
+       CCTK_WARN(0, "could not set LIMIT value in table!");
+    if (Util_TableSetReal(handle_g13_bound ,g13_bound_speed, "SPEED") < 0)
+        CCTK_WARN(0, "could not set SPEED value in table!");
+  
+    ierr = Boundary_SelectVarForBC(cctkGH, CCTK_ALL_FACES, 1, handle_g13_bound, 
+                      "ML_ADM::g13", "Radiation");
+  
+    if (ierr < 0)
+       CCTK_WARN(0, "Failed to register Radiation BC for ML_ADM::g13!");
+  
+  }
+  
+  if (CCTK_EQUALS(g22_bound, "radiative"))
+  {
+   /* select radiation boundary condition */
+    static CCTK_INT handle_g22_bound = -1;
+    if (handle_g22_bound < 0) handle_g22_bound = Util_TableCreate(UTIL_TABLE_FLAGS_CASE_INSENSITIVE);
+    if (handle_g22_bound < 0) CCTK_WARN(0, "could not create table!");
+    if (Util_TableSetReal(handle_g22_bound , g22_bound_limit, "LIMIT") < 0)
+       CCTK_WARN(0, "could not set LIMIT value in table!");
+    if (Util_TableSetReal(handle_g22_bound ,g22_bound_speed, "SPEED") < 0)
+        CCTK_WARN(0, "could not set SPEED value in table!");
+  
+    ierr = Boundary_SelectVarForBC(cctkGH, CCTK_ALL_FACES, 1, handle_g22_bound, 
+                      "ML_ADM::g22", "Radiation");
+  
+    if (ierr < 0)
+       CCTK_WARN(0, "Failed to register Radiation BC for ML_ADM::g22!");
+  
+  }
+  
+  if (CCTK_EQUALS(g23_bound, "radiative"))
+  {
+   /* select radiation boundary condition */
+    static CCTK_INT handle_g23_bound = -1;
+    if (handle_g23_bound < 0) handle_g23_bound = Util_TableCreate(UTIL_TABLE_FLAGS_CASE_INSENSITIVE);
+    if (handle_g23_bound < 0) CCTK_WARN(0, "could not create table!");
+    if (Util_TableSetReal(handle_g23_bound , g23_bound_limit, "LIMIT") < 0)
+       CCTK_WARN(0, "could not set LIMIT value in table!");
+    if (Util_TableSetReal(handle_g23_bound ,g23_bound_speed, "SPEED") < 0)
+        CCTK_WARN(0, "could not set SPEED value in table!");
+  
+    ierr = Boundary_SelectVarForBC(cctkGH, CCTK_ALL_FACES, 1, handle_g23_bound, 
+                      "ML_ADM::g23", "Radiation");
+  
+    if (ierr < 0)
+       CCTK_WARN(0, "Failed to register Radiation BC for ML_ADM::g23!");
+  
+  }
+  
+  if (CCTK_EQUALS(g33_bound, "radiative"))
+  {
+   /* select radiation boundary condition */
+    static CCTK_INT handle_g33_bound = -1;
+    if (handle_g33_bound < 0) handle_g33_bound = Util_TableCreate(UTIL_TABLE_FLAGS_CASE_INSENSITIVE);
+    if (handle_g33_bound < 0) CCTK_WARN(0, "could not create table!");
+    if (Util_TableSetReal(handle_g33_bound , g33_bound_limit, "LIMIT") < 0)
+       CCTK_WARN(0, "could not set LIMIT value in table!");
+    if (Util_TableSetReal(handle_g33_bound ,g33_bound_speed, "SPEED") < 0)
+        CCTK_WARN(0, "could not set SPEED value in table!");
+  
+    ierr = Boundary_SelectVarForBC(cctkGH, CCTK_ALL_FACES, 1, handle_g33_bound, 
+                      "ML_ADM::g33", "Radiation");
+  
+    if (ierr < 0)
+       CCTK_WARN(0, "Failed to register Radiation BC for ML_ADM::g33!");
+  
+  }
+  
+  if (CCTK_EQUALS(beta1_bound, "radiative"))
+  {
+   /* select radiation boundary condition */
+    static CCTK_INT handle_beta1_bound = -1;
+    if (handle_beta1_bound < 0) handle_beta1_bound = Util_TableCreate(UTIL_TABLE_FLAGS_CASE_INSENSITIVE);
+    if (handle_beta1_bound < 0) CCTK_WARN(0, "could not create table!");
+    if (Util_TableSetReal(handle_beta1_bound , beta1_bound_limit, "LIMIT") < 0)
+       CCTK_WARN(0, "could not set LIMIT value in table!");
+    if (Util_TableSetReal(handle_beta1_bound ,beta1_bound_speed, "SPEED") < 0)
+        CCTK_WARN(0, "could not set SPEED value in table!");
+  
+    ierr = Boundary_SelectVarForBC(cctkGH, CCTK_ALL_FACES, 1, handle_beta1_bound, 
+                      "ML_ADM::beta1", "Radiation");
+  
+    if (ierr < 0)
+       CCTK_WARN(0, "Failed to register Radiation BC for ML_ADM::beta1!");
+  
+  }
+  
+  if (CCTK_EQUALS(beta2_bound, "radiative"))
+  {
+   /* select radiation boundary condition */
+    static CCTK_INT handle_beta2_bound = -1;
+    if (handle_beta2_bound < 0) handle_beta2_bound = Util_TableCreate(UTIL_TABLE_FLAGS_CASE_INSENSITIVE);
+    if (handle_beta2_bound < 0) CCTK_WARN(0, "could not create table!");
+    if (Util_TableSetReal(handle_beta2_bound , beta2_bound_limit, "LIMIT") < 0)
+       CCTK_WARN(0, "could not set LIMIT value in table!");
+    if (Util_TableSetReal(handle_beta2_bound ,beta2_bound_speed, "SPEED") < 0)
+        CCTK_WARN(0, "could not set SPEED value in table!");
+  
+    ierr = Boundary_SelectVarForBC(cctkGH, CCTK_ALL_FACES, 1, handle_beta2_bound, 
+                      "ML_ADM::beta2", "Radiation");
+  
+    if (ierr < 0)
+       CCTK_WARN(0, "Failed to register Radiation BC for ML_ADM::beta2!");
+  
+  }
+  
+  if (CCTK_EQUALS(beta3_bound, "radiative"))
+  {
+   /* select radiation boundary condition */
+    static CCTK_INT handle_beta3_bound = -1;
+    if (handle_beta3_bound < 0) handle_beta3_bound = Util_TableCreate(UTIL_TABLE_FLAGS_CASE_INSENSITIVE);
+    if (handle_beta3_bound < 0) CCTK_WARN(0, "could not create table!");
+    if (Util_TableSetReal(handle_beta3_bound , beta3_bound_limit, "LIMIT") < 0)
+       CCTK_WARN(0, "could not set LIMIT value in table!");
+    if (Util_TableSetReal(handle_beta3_bound ,beta3_bound_speed, "SPEED") < 0)
+        CCTK_WARN(0, "could not set SPEED value in table!");
+  
+    ierr = Boundary_SelectVarForBC(cctkGH, CCTK_ALL_FACES, 1, handle_beta3_bound, 
+                      "ML_ADM::beta3", "Radiation");
+  
+    if (ierr < 0)
+       CCTK_WARN(0, "Failed to register Radiation BC for ML_ADM::beta3!");
+  
+  }
+  
+  if (CCTK_EQUALS(ML_curv_bound, "scalar"))
+  {
+   /* select scalar boundary condition */
+    static CCTK_INT handle_ML_curv_bound = -1;
+    if (handle_ML_curv_bound < 0) handle_ML_curv_bound = Util_TableCreate(UTIL_TABLE_FLAGS_CASE_INSENSITIVE);
+    if (handle_ML_curv_bound < 0) CCTK_WARN(0, "could not create table!");
+    if (Util_TableSetReal(handle_ML_curv_bound ,ML_curv_bound_scalar, "SCALAR") < 0)
+        CCTK_WARN(0, "could not set SCALAR value in table!");
+  
+    ierr = Boundary_SelectGroupForBC(cctkGH, CCTK_ALL_FACES, 1, handle_ML_curv_bound, 
+                      "ML_ADM::ML_curv", "scalar");
+  
+    if (ierr < 0)
+       CCTK_WARN(0, "Failed to register Scalar BC for ML_ADM::ML_curv!");
+  
+  }
+  
+  if (CCTK_EQUALS(ML_lapse_bound, "scalar"))
+  {
+   /* select scalar boundary condition */
+    static CCTK_INT handle_ML_lapse_bound = -1;
+    if (handle_ML_lapse_bound < 0) handle_ML_lapse_bound = Util_TableCreate(UTIL_TABLE_FLAGS_CASE_INSENSITIVE);
+    if (handle_ML_lapse_bound < 0) CCTK_WARN(0, "could not create table!");
+    if (Util_TableSetReal(handle_ML_lapse_bound ,ML_lapse_bound_scalar, "SCALAR") < 0)
+        CCTK_WARN(0, "could not set SCALAR value in table!");
+  
+    ierr = Boundary_SelectGroupForBC(cctkGH, CCTK_ALL_FACES, 1, handle_ML_lapse_bound, 
+                      "ML_ADM::ML_lapse", "scalar");
+  
+    if (ierr < 0)
+       CCTK_WARN(0, "Failed to register Scalar BC for ML_ADM::ML_lapse!");
+  
+  }
+  
+  if (CCTK_EQUALS(ML_metric_bound, "scalar"))
+  {
+   /* select scalar boundary condition */
+    static CCTK_INT handle_ML_metric_bound = -1;
+    if (handle_ML_metric_bound < 0) handle_ML_metric_bound = Util_TableCreate(UTIL_TABLE_FLAGS_CASE_INSENSITIVE);
+    if (handle_ML_metric_bound < 0) CCTK_WARN(0, "could not create table!");
+    if (Util_TableSetReal(handle_ML_metric_bound ,ML_metric_bound_scalar, "SCALAR") < 0)
+        CCTK_WARN(0, "could not set SCALAR value in table!");
+  
+    ierr = Boundary_SelectGroupForBC(cctkGH, CCTK_ALL_FACES, 1, handle_ML_metric_bound, 
+                      "ML_ADM::ML_metric", "scalar");
+  
+    if (ierr < 0)
+       CCTK_WARN(0, "Failed to register Scalar BC for ML_ADM::ML_metric!");
+  
+  }
+  
+  if (CCTK_EQUALS(ML_shift_bound, "scalar"))
+  {
+   /* select scalar boundary condition */
+    static CCTK_INT handle_ML_shift_bound = -1;
+    if (handle_ML_shift_bound < 0) handle_ML_shift_bound = Util_TableCreate(UTIL_TABLE_FLAGS_CASE_INSENSITIVE);
+    if (handle_ML_shift_bound < 0) CCTK_WARN(0, "could not create table!");
+    if (Util_TableSetReal(handle_ML_shift_bound ,ML_shift_bound_scalar, "SCALAR") < 0)
+        CCTK_WARN(0, "could not set SCALAR value in table!");
+  
+    ierr = Boundary_SelectGroupForBC(cctkGH, CCTK_ALL_FACES, 1, handle_ML_shift_bound, 
+                      "ML_ADM::ML_shift", "scalar");
+  
+    if (ierr < 0)
+       CCTK_WARN(0, "Failed to register Scalar BC for ML_ADM::ML_shift!");
+  
+  }
+  
+  if (CCTK_EQUALS(K11_bound, "scalar"))
+  {
+   /* select scalar boundary condition */
+    static CCTK_INT handle_K11_bound = -1;
+    if (handle_K11_bound < 0) handle_K11_bound = Util_TableCreate(UTIL_TABLE_FLAGS_CASE_INSENSITIVE);
+    if (handle_K11_bound < 0) CCTK_WARN(0, "could not create table!");
+    if (Util_TableSetReal(handle_K11_bound ,K11_bound_scalar, "SCALAR") < 0)
+      CCTK_WARN(0, "could not set SCALAR value in table!");
+  
+    ierr = Boundary_SelectVarForBC(cctkGH, CCTK_ALL_FACES, 1, handle_K11_bound, 
+                      "ML_ADM::K11", "scalar");
+  
+    if (ierr < 0)
+       CCTK_WARN(0, "Error in registering Scalar BC for ML_ADM::K11!");
+  
+  }
+  
+  if (CCTK_EQUALS(K12_bound, "scalar"))
+  {
+   /* select scalar boundary condition */
+    static CCTK_INT handle_K12_bound = -1;
+    if (handle_K12_bound < 0) handle_K12_bound = Util_TableCreate(UTIL_TABLE_FLAGS_CASE_INSENSITIVE);
+    if (handle_K12_bound < 0) CCTK_WARN(0, "could not create table!");
+    if (Util_TableSetReal(handle_K12_bound ,K12_bound_scalar, "SCALAR") < 0)
+      CCTK_WARN(0, "could not set SCALAR value in table!");
+  
+    ierr = Boundary_SelectVarForBC(cctkGH, CCTK_ALL_FACES, 1, handle_K12_bound, 
+                      "ML_ADM::K12", "scalar");
+  
+    if (ierr < 0)
+       CCTK_WARN(0, "Error in registering Scalar BC for ML_ADM::K12!");
+  
+  }
+  
+  if (CCTK_EQUALS(K13_bound, "scalar"))
+  {
+   /* select scalar boundary condition */
+    static CCTK_INT handle_K13_bound = -1;
+    if (handle_K13_bound < 0) handle_K13_bound = Util_TableCreate(UTIL_TABLE_FLAGS_CASE_INSENSITIVE);
+    if (handle_K13_bound < 0) CCTK_WARN(0, "could not create table!");
+    if (Util_TableSetReal(handle_K13_bound ,K13_bound_scalar, "SCALAR") < 0)
+      CCTK_WARN(0, "could not set SCALAR value in table!");
+  
+    ierr = Boundary_SelectVarForBC(cctkGH, CCTK_ALL_FACES, 1, handle_K13_bound, 
+                      "ML_ADM::K13", "scalar");
+  
+    if (ierr < 0)
+       CCTK_WARN(0, "Error in registering Scalar BC for ML_ADM::K13!");
+  
+  }
+  
+  if (CCTK_EQUALS(K22_bound, "scalar"))
+  {
+   /* select scalar boundary condition */
+    static CCTK_INT handle_K22_bound = -1;
+    if (handle_K22_bound < 0) handle_K22_bound = Util_TableCreate(UTIL_TABLE_FLAGS_CASE_INSENSITIVE);
+    if (handle_K22_bound < 0) CCTK_WARN(0, "could not create table!");
+    if (Util_TableSetReal(handle_K22_bound ,K22_bound_scalar, "SCALAR") < 0)
+      CCTK_WARN(0, "could not set SCALAR value in table!");
+  
+    ierr = Boundary_SelectVarForBC(cctkGH, CCTK_ALL_FACES, 1, handle_K22_bound, 
+                      "ML_ADM::K22", "scalar");
+  
+    if (ierr < 0)
+       CCTK_WARN(0, "Error in registering Scalar BC for ML_ADM::K22!");
+  
+  }
+  
+  if (CCTK_EQUALS(K23_bound, "scalar"))
+  {
+   /* select scalar boundary condition */
+    static CCTK_INT handle_K23_bound = -1;
+    if (handle_K23_bound < 0) handle_K23_bound = Util_TableCreate(UTIL_TABLE_FLAGS_CASE_INSENSITIVE);
+    if (handle_K23_bound < 0) CCTK_WARN(0, "could not create table!");
+    if (Util_TableSetReal(handle_K23_bound ,K23_bound_scalar, "SCALAR") < 0)
+      CCTK_WARN(0, "could not set SCALAR value in table!");
+  
+    ierr = Boundary_SelectVarForBC(cctkGH, CCTK_ALL_FACES, 1, handle_K23_bound, 
+                      "ML_ADM::K23", "scalar");
+  
+    if (ierr < 0)
+       CCTK_WARN(0, "Error in registering Scalar BC for ML_ADM::K23!");
+  
+  }
+  
+  if (CCTK_EQUALS(K33_bound, "scalar"))
+  {
+   /* select scalar boundary condition */
+    static CCTK_INT handle_K33_bound = -1;
+    if (handle_K33_bound < 0) handle_K33_bound = Util_TableCreate(UTIL_TABLE_FLAGS_CASE_INSENSITIVE);
+    if (handle_K33_bound < 0) CCTK_WARN(0, "could not create table!");
+    if (Util_TableSetReal(handle_K33_bound ,K33_bound_scalar, "SCALAR") < 0)
+      CCTK_WARN(0, "could not set SCALAR value in table!");
+  
+    ierr = Boundary_SelectVarForBC(cctkGH, CCTK_ALL_FACES, 1, handle_K33_bound, 
+                      "ML_ADM::K33", "scalar");
+  
+    if (ierr < 0)
+       CCTK_WARN(0, "Error in registering Scalar BC for ML_ADM::K33!");
+  
+  }
+  
+  if (CCTK_EQUALS(alpha_bound, "scalar"))
+  {
+   /* select scalar boundary condition */
+    static CCTK_INT handle_alpha_bound = -1;
+    if (handle_alpha_bound < 0) handle_alpha_bound = Util_TableCreate(UTIL_TABLE_FLAGS_CASE_INSENSITIVE);
+    if (handle_alpha_bound < 0) CCTK_WARN(0, "could not create table!");
+    if (Util_TableSetReal(handle_alpha_bound ,alpha_bound_scalar, "SCALAR") < 0)
+      CCTK_WARN(0, "could not set SCALAR value in table!");
+  
+    ierr = Boundary_SelectVarForBC(cctkGH, CCTK_ALL_FACES, 1, handle_alpha_bound, 
+                      "ML_ADM::alpha", "scalar");
+  
+    if (ierr < 0)
+       CCTK_WARN(0, "Error in registering Scalar BC for ML_ADM::alpha!");
+  
+  }
+  
+  if (CCTK_EQUALS(g11_bound, "scalar"))
+  {
+   /* select scalar boundary condition */
+    static CCTK_INT handle_g11_bound = -1;
+    if (handle_g11_bound < 0) handle_g11_bound = Util_TableCreate(UTIL_TABLE_FLAGS_CASE_INSENSITIVE);
+    if (handle_g11_bound < 0) CCTK_WARN(0, "could not create table!");
+    if (Util_TableSetReal(handle_g11_bound ,g11_bound_scalar, "SCALAR") < 0)
+      CCTK_WARN(0, "could not set SCALAR value in table!");
+  
+    ierr = Boundary_SelectVarForBC(cctkGH, CCTK_ALL_FACES, 1, handle_g11_bound, 
+                      "ML_ADM::g11", "scalar");
+  
+    if (ierr < 0)
+       CCTK_WARN(0, "Error in registering Scalar BC for ML_ADM::g11!");
+  
+  }
+  
+  if (CCTK_EQUALS(g12_bound, "scalar"))
+  {
+   /* select scalar boundary condition */
+    static CCTK_INT handle_g12_bound = -1;
+    if (handle_g12_bound < 0) handle_g12_bound = Util_TableCreate(UTIL_TABLE_FLAGS_CASE_INSENSITIVE);
+    if (handle_g12_bound < 0) CCTK_WARN(0, "could not create table!");
+    if (Util_TableSetReal(handle_g12_bound ,g12_bound_scalar, "SCALAR") < 0)
+      CCTK_WARN(0, "could not set SCALAR value in table!");
+  
+    ierr = Boundary_SelectVarForBC(cctkGH, CCTK_ALL_FACES, 1, handle_g12_bound, 
+                      "ML_ADM::g12", "scalar");
+  
+    if (ierr < 0)
+       CCTK_WARN(0, "Error in registering Scalar BC for ML_ADM::g12!");
+  
+  }
+  
+  if (CCTK_EQUALS(g13_bound, "scalar"))
+  {
+   /* select scalar boundary condition */
+    static CCTK_INT handle_g13_bound = -1;
+    if (handle_g13_bound < 0) handle_g13_bound = Util_TableCreate(UTIL_TABLE_FLAGS_CASE_INSENSITIVE);
+    if (handle_g13_bound < 0) CCTK_WARN(0, "could not create table!");
+    if (Util_TableSetReal(handle_g13_bound ,g13_bound_scalar, "SCALAR") < 0)
+      CCTK_WARN(0, "could not set SCALAR value in table!");
+  
+    ierr = Boundary_SelectVarForBC(cctkGH, CCTK_ALL_FACES, 1, handle_g13_bound, 
+                      "ML_ADM::g13", "scalar");
+  
+    if (ierr < 0)
+       CCTK_WARN(0, "Error in registering Scalar BC for ML_ADM::g13!");
+  
+  }
+  
+  if (CCTK_EQUALS(g22_bound, "scalar"))
+  {
+   /* select scalar boundary condition */
+    static CCTK_INT handle_g22_bound = -1;
+    if (handle_g22_bound < 0) handle_g22_bound = Util_TableCreate(UTIL_TABLE_FLAGS_CASE_INSENSITIVE);
+    if (handle_g22_bound < 0) CCTK_WARN(0, "could not create table!");
+    if (Util_TableSetReal(handle_g22_bound ,g22_bound_scalar, "SCALAR") < 0)
+      CCTK_WARN(0, "could not set SCALAR value in table!");
+  
+    ierr = Boundary_SelectVarForBC(cctkGH, CCTK_ALL_FACES, 1, handle_g22_bound, 
+                      "ML_ADM::g22", "scalar");
+  
+    if (ierr < 0)
+       CCTK_WARN(0, "Error in registering Scalar BC for ML_ADM::g22!");
+  
+  }
+  
+  if (CCTK_EQUALS(g23_bound, "scalar"))
+  {
+   /* select scalar boundary condition */
+    static CCTK_INT handle_g23_bound = -1;
+    if (handle_g23_bound < 0) handle_g23_bound = Util_TableCreate(UTIL_TABLE_FLAGS_CASE_INSENSITIVE);
+    if (handle_g23_bound < 0) CCTK_WARN(0, "could not create table!");
+    if (Util_TableSetReal(handle_g23_bound ,g23_bound_scalar, "SCALAR") < 0)
+      CCTK_WARN(0, "could not set SCALAR value in table!");
+  
+    ierr = Boundary_SelectVarForBC(cctkGH, CCTK_ALL_FACES, 1, handle_g23_bound, 
+                      "ML_ADM::g23", "scalar");
+  
+    if (ierr < 0)
+       CCTK_WARN(0, "Error in registering Scalar BC for ML_ADM::g23!");
+  
+  }
+  
+  if (CCTK_EQUALS(g33_bound, "scalar"))
+  {
+   /* select scalar boundary condition */
+    static CCTK_INT handle_g33_bound = -1;
+    if (handle_g33_bound < 0) handle_g33_bound = Util_TableCreate(UTIL_TABLE_FLAGS_CASE_INSENSITIVE);
+    if (handle_g33_bound < 0) CCTK_WARN(0, "could not create table!");
+    if (Util_TableSetReal(handle_g33_bound ,g33_bound_scalar, "SCALAR") < 0)
+      CCTK_WARN(0, "could not set SCALAR value in table!");
+  
+    ierr = Boundary_SelectVarForBC(cctkGH, CCTK_ALL_FACES, 1, handle_g33_bound, 
+                      "ML_ADM::g33", "scalar");
+  
+    if (ierr < 0)
+       CCTK_WARN(0, "Error in registering Scalar BC for ML_ADM::g33!");
+  
+  }
+  
+  if (CCTK_EQUALS(beta1_bound, "scalar"))
+  {
+   /* select scalar boundary condition */
+    static CCTK_INT handle_beta1_bound = -1;
+    if (handle_beta1_bound < 0) handle_beta1_bound = Util_TableCreate(UTIL_TABLE_FLAGS_CASE_INSENSITIVE);
+    if (handle_beta1_bound < 0) CCTK_WARN(0, "could not create table!");
+    if (Util_TableSetReal(handle_beta1_bound ,beta1_bound_scalar, "SCALAR") < 0)
+      CCTK_WARN(0, "could not set SCALAR value in table!");
+  
+    ierr = Boundary_SelectVarForBC(cctkGH, CCTK_ALL_FACES, 1, handle_beta1_bound, 
+                      "ML_ADM::beta1", "scalar");
+  
+    if (ierr < 0)
+       CCTK_WARN(0, "Error in registering Scalar BC for ML_ADM::beta1!");
+  
+  }
+  
+  if (CCTK_EQUALS(beta2_bound, "scalar"))
+  {
+   /* select scalar boundary condition */
+    static CCTK_INT handle_beta2_bound = -1;
+    if (handle_beta2_bound < 0) handle_beta2_bound = Util_TableCreate(UTIL_TABLE_FLAGS_CASE_INSENSITIVE);
+    if (handle_beta2_bound < 0) CCTK_WARN(0, "could not create table!");
+    if (Util_TableSetReal(handle_beta2_bound ,beta2_bound_scalar, "SCALAR") < 0)
+      CCTK_WARN(0, "could not set SCALAR value in table!");
+  
+    ierr = Boundary_SelectVarForBC(cctkGH, CCTK_ALL_FACES, 1, handle_beta2_bound, 
+                      "ML_ADM::beta2", "scalar");
+  
+    if (ierr < 0)
+       CCTK_WARN(0, "Error in registering Scalar BC for ML_ADM::beta2!");
+  
+  }
+  
+  if (CCTK_EQUALS(beta3_bound, "scalar"))
+  {
+   /* select scalar boundary condition */
+    static CCTK_INT handle_beta3_bound = -1;
+    if (handle_beta3_bound < 0) handle_beta3_bound = Util_TableCreate(UTIL_TABLE_FLAGS_CASE_INSENSITIVE);
+    if (handle_beta3_bound < 0) CCTK_WARN(0, "could not create table!");
+    if (Util_TableSetReal(handle_beta3_bound ,beta3_bound_scalar, "SCALAR") < 0)
+      CCTK_WARN(0, "could not set SCALAR value in table!");
+  
+    ierr = Boundary_SelectVarForBC(cctkGH, CCTK_ALL_FACES, 1, handle_beta3_bound, 
+                      "ML_ADM::beta3", "scalar");
+  
+    if (ierr < 0)
+       CCTK_WARN(0, "Error in registering Scalar BC for ML_ADM::beta3!");
+  
+  }
+  return;
+}
+
+
+
+/* template for entries in parameter file:
+#$bound$#ML_ADM::ML_curv_bound       = "skip"
+#$bound$#ML_ADM::ML_curv_bound_speed = 1.0
+#$bound$#ML_ADM::ML_curv_bound_limit = 0.0
+#$bound$#ML_ADM::ML_curv_bound_scalar = 0.0
+
+#$bound$#ML_ADM::ML_lapse_bound       = "skip"
+#$bound$#ML_ADM::ML_lapse_bound_speed = 1.0
+#$bound$#ML_ADM::ML_lapse_bound_limit = 0.0
+#$bound$#ML_ADM::ML_lapse_bound_scalar = 0.0
+
+#$bound$#ML_ADM::ML_metric_bound       = "skip"
+#$bound$#ML_ADM::ML_metric_bound_speed = 1.0
+#$bound$#ML_ADM::ML_metric_bound_limit = 0.0
+#$bound$#ML_ADM::ML_metric_bound_scalar = 0.0
+
+#$bound$#ML_ADM::ML_shift_bound       = "skip"
+#$bound$#ML_ADM::ML_shift_bound_speed = 1.0
+#$bound$#ML_ADM::ML_shift_bound_limit = 0.0
+#$bound$#ML_ADM::ML_shift_bound_scalar = 0.0
+
+#$bound$#ML_ADM::K11_bound       = "skip"
+#$bound$#ML_ADM::K11_bound_speed = 1.0
+#$bound$#ML_ADM::K11_bound_limit = 0.0
+#$bound$#ML_ADM::K11_bound_scalar = 0.0
+
+#$bound$#ML_ADM::K12_bound       = "skip"
+#$bound$#ML_ADM::K12_bound_speed = 1.0
+#$bound$#ML_ADM::K12_bound_limit = 0.0
+#$bound$#ML_ADM::K12_bound_scalar = 0.0
+
+#$bound$#ML_ADM::K13_bound       = "skip"
+#$bound$#ML_ADM::K13_bound_speed = 1.0
+#$bound$#ML_ADM::K13_bound_limit = 0.0
+#$bound$#ML_ADM::K13_bound_scalar = 0.0
+
+#$bound$#ML_ADM::K22_bound       = "skip"
+#$bound$#ML_ADM::K22_bound_speed = 1.0
+#$bound$#ML_ADM::K22_bound_limit = 0.0
+#$bound$#ML_ADM::K22_bound_scalar = 0.0
+
+#$bound$#ML_ADM::K23_bound       = "skip"
+#$bound$#ML_ADM::K23_bound_speed = 1.0
+#$bound$#ML_ADM::K23_bound_limit = 0.0
+#$bound$#ML_ADM::K23_bound_scalar = 0.0
+
+#$bound$#ML_ADM::K33_bound       = "skip"
+#$bound$#ML_ADM::K33_bound_speed = 1.0
+#$bound$#ML_ADM::K33_bound_limit = 0.0
+#$bound$#ML_ADM::K33_bound_scalar = 0.0
+
+#$bound$#ML_ADM::alpha_bound       = "skip"
+#$bound$#ML_ADM::alpha_bound_speed = 1.0
+#$bound$#ML_ADM::alpha_bound_limit = 0.0
+#$bound$#ML_ADM::alpha_bound_scalar = 0.0
+
+#$bound$#ML_ADM::g11_bound       = "skip"
+#$bound$#ML_ADM::g11_bound_speed = 1.0
+#$bound$#ML_ADM::g11_bound_limit = 0.0
+#$bound$#ML_ADM::g11_bound_scalar = 0.0
+
+#$bound$#ML_ADM::g12_bound       = "skip"
+#$bound$#ML_ADM::g12_bound_speed = 1.0
+#$bound$#ML_ADM::g12_bound_limit = 0.0
+#$bound$#ML_ADM::g12_bound_scalar = 0.0
+
+#$bound$#ML_ADM::g13_bound       = "skip"
+#$bound$#ML_ADM::g13_bound_speed = 1.0
+#$bound$#ML_ADM::g13_bound_limit = 0.0
+#$bound$#ML_ADM::g13_bound_scalar = 0.0
+
+#$bound$#ML_ADM::g22_bound       = "skip"
+#$bound$#ML_ADM::g22_bound_speed = 1.0
+#$bound$#ML_ADM::g22_bound_limit = 0.0
+#$bound$#ML_ADM::g22_bound_scalar = 0.0
+
+#$bound$#ML_ADM::g23_bound       = "skip"
+#$bound$#ML_ADM::g23_bound_speed = 1.0
+#$bound$#ML_ADM::g23_bound_limit = 0.0
+#$bound$#ML_ADM::g23_bound_scalar = 0.0
+
+#$bound$#ML_ADM::g33_bound       = "skip"
+#$bound$#ML_ADM::g33_bound_speed = 1.0
+#$bound$#ML_ADM::g33_bound_limit = 0.0
+#$bound$#ML_ADM::g33_bound_scalar = 0.0
+
+#$bound$#ML_ADM::beta1_bound       = "skip"
+#$bound$#ML_ADM::beta1_bound_speed = 1.0
+#$bound$#ML_ADM::beta1_bound_limit = 0.0
+#$bound$#ML_ADM::beta1_bound_scalar = 0.0
+
+#$bound$#ML_ADM::beta2_bound       = "skip"
+#$bound$#ML_ADM::beta2_bound_speed = 1.0
+#$bound$#ML_ADM::beta2_bound_limit = 0.0
+#$bound$#ML_ADM::beta2_bound_scalar = 0.0
+
+#$bound$#ML_ADM::beta3_bound       = "skip"
+#$bound$#ML_ADM::beta3_bound_speed = 1.0
+#$bound$#ML_ADM::beta3_bound_limit = 0.0
+#$bound$#ML_ADM::beta3_bound_scalar = 0.0
+
+*/
+
diff --git a/ML_ADM/src/ML_ADM_Minkowski.cc b/ML_ADM/src/ML_ADM_Minkowski.cc
new file mode 100644
index 0000000..59b9661
--- /dev/null
+++ b/ML_ADM/src/ML_ADM_Minkowski.cc
@@ -0,0 +1,151 @@
+/*  File produced by Kranc */
+
+#define KRANC_C
+
+#include <assert.h>
+#include <math.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include "cctk.h"
+#include "cctk_Arguments.h"
+#include "cctk_Parameters.h"
+#include "GenericFD.h"
+#include "Differencing.h"
+#include "loopcontrol.h"
+
+/* Define macros used in calculations */
+#define INITVALUE (42)
+#define QAD(x) (SQR(SQR(x)))
+#define INV(x) ((1.0) / (x))
+#define SQR(x) ((x) * (x))
+#define CUB(x) ((x) * (x) * (x))
+
+static void ML_ADM_Minkowski_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[])
+{
+  DECLARE_CCTK_ARGUMENTS;
+  DECLARE_CCTK_PARAMETERS;
+  
+  
+  /* Declare finite differencing variables */
+  
+  if (verbose > 1)
+  {
+    CCTK_VInfo(CCTK_THORNSTRING,"Entering ML_ADM_Minkowski_Body");
+  }
+  
+  if (cctk_iteration % ML_ADM_Minkowski_calc_every != ML_ADM_Minkowski_calc_offset)
+  {
+    return;
+  }
+  
+  const char *groups[] = {"ML_ADM::ML_curv","ML_ADM::ML_lapse","ML_ADM::ML_metric","ML_ADM::ML_shift"};
+  GenericFD_AssertGroupStorage(cctkGH, "ML_ADM_Minkowski", 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);
+  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 dxi = INV(dx);
+  CCTK_REAL const dyi = INV(dy);
+  CCTK_REAL const dzi = INV(dz);
+  CCTK_REAL const khalf = 0.5;
+  CCTK_REAL const kthird = 1/3.0;
+  CCTK_REAL const ktwothird = 2.0/3.0;
+  CCTK_REAL const kfourthird = 4.0/3.0;
+  CCTK_REAL const keightthird = 8.0/3.0;
+  CCTK_REAL const hdxi = 0.5 * dxi;
+  CCTK_REAL const hdyi = 0.5 * dyi;
+  CCTK_REAL const hdzi = 0.5 * dzi;
+  
+  /* Initialize predefined quantities */
+  CCTK_REAL const p1o12dx = 0.0833333333333333333333333333333*INV(dx);
+  CCTK_REAL const p1o12dy = 0.0833333333333333333333333333333*INV(dy);
+  CCTK_REAL const p1o12dz = 0.0833333333333333333333333333333*INV(dz);
+  CCTK_REAL const p1o144dxdy = 0.00694444444444444444444444444444*INV(dx)*INV(dy);
+  CCTK_REAL const p1o144dxdz = 0.00694444444444444444444444444444*INV(dx)*INV(dz);
+  CCTK_REAL const p1o144dydz = 0.00694444444444444444444444444444*INV(dy)*INV(dz);
+  CCTK_REAL const pm1o12dx2 = -0.0833333333333333333333333333333*INV(SQR(dx));
+  CCTK_REAL const pm1o12dy2 = -0.0833333333333333333333333333333*INV(SQR(dy));
+  CCTK_REAL const pm1o12dz2 = -0.0833333333333333333333333333333*INV(SQR(dz));
+  
+  /* Loop over the grid points */
+  #pragma omp parallel
+  LC_LOOP3 (ML_ADM_Minkowski,
+    i,j,k, min[0],min[1],min[2], max[0],max[1],max[2],
+    cctk_lsh[0],cctk_lsh[1],cctk_lsh[2])
+  {
+    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 g11L = 1;
+    
+    CCTK_REAL g12L = 0;
+    
+    CCTK_REAL g13L = 0;
+    
+    CCTK_REAL g22L = 1;
+    
+    CCTK_REAL g23L = 0;
+    
+    CCTK_REAL g33L = 1;
+    
+    CCTK_REAL K11L = 0;
+    
+    CCTK_REAL K12L = 0;
+    
+    CCTK_REAL K13L = 0;
+    
+    CCTK_REAL K22L = 0;
+    
+    CCTK_REAL K23L = 0;
+    
+    CCTK_REAL K33L = 0;
+    
+    CCTK_REAL alphaL = 1;
+    
+    CCTK_REAL beta1L = 0;
+    
+    CCTK_REAL beta2L = 0;
+    
+    CCTK_REAL beta3L = 0;
+    
+    
+    /* Copy local copies back to grid functions */
+    alpha[index] = alphaL;
+    beta1[index] = beta1L;
+    beta2[index] = beta2L;
+    beta3[index] = beta3L;
+    g11[index] = g11L;
+    g12[index] = g12L;
+    g13[index] = g13L;
+    g22[index] = g22L;
+    g23[index] = g23L;
+    g33[index] = g33L;
+    K11[index] = K11L;
+    K12[index] = K12L;
+    K13[index] = K13L;
+    K22[index] = K22L;
+    K23[index] = K23L;
+    K33[index] = K33L;
+  }
+  LC_ENDLOOP3 (ML_ADM_Minkowski);
+}
+
+extern "C" void ML_ADM_Minkowski(CCTK_ARGUMENTS)
+{
+  DECLARE_CCTK_ARGUMENTS;
+  DECLARE_CCTK_PARAMETERS;
+  
+  GenericFD_LoopOverEverything(cctkGH, &ML_ADM_Minkowski_Body);
+}
diff --git a/ML_ADM/src/ML_ADM_RHS.cc b/ML_ADM/src/ML_ADM_RHS.cc
new file mode 100644
index 0000000..a03953e
--- /dev/null
+++ b/ML_ADM/src/ML_ADM_RHS.cc
@@ -0,0 +1,482 @@
+/*  File produced by Kranc */
+
+#define KRANC_C
+
+#include <assert.h>
+#include <math.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include "cctk.h"
+#include "cctk_Arguments.h"
+#include "cctk_Parameters.h"
+#include "GenericFD.h"
+#include "Differencing.h"
+#include "loopcontrol.h"
+
+/* Define macros used in calculations */
+#define INITVALUE (42)
+#define QAD(x) (SQR(SQR(x)))
+#define INV(x) ((1.0) / (x))
+#define SQR(x) ((x) * (x))
+#define CUB(x) ((x) * (x) * (x))
+
+extern "C" void ML_ADM_RHS_SelectBCs(CCTK_ARGUMENTS)
+{
+  DECLARE_CCTK_ARGUMENTS;
+  DECLARE_CCTK_PARAMETERS;
+  
+  CCTK_INT ierr = 0;
+  ierr = Boundary_SelectGroupForBC(cctkGH, CCTK_ALL_FACES, GenericFD_GetBoundaryWidth(cctkGH), -1 /* no table */, "ML_ADM::ML_curvrhs","flat");
+  if (ierr < 0)
+    CCTK_WARN(1, "Failed to register flat BC for ML_ADM::ML_curvrhs.");
+  ierr = Boundary_SelectGroupForBC(cctkGH, CCTK_ALL_FACES, GenericFD_GetBoundaryWidth(cctkGH), -1 /* no table */, "ML_ADM::ML_lapserhs","flat");
+  if (ierr < 0)
+    CCTK_WARN(1, "Failed to register flat BC for ML_ADM::ML_lapserhs.");
+  ierr = Boundary_SelectGroupForBC(cctkGH, CCTK_ALL_FACES, GenericFD_GetBoundaryWidth(cctkGH), -1 /* no table */, "ML_ADM::ML_metricrhs","flat");
+  if (ierr < 0)
+    CCTK_WARN(1, "Failed to register flat BC for ML_ADM::ML_metricrhs.");
+  ierr = Boundary_SelectGroupForBC(cctkGH, CCTK_ALL_FACES, GenericFD_GetBoundaryWidth(cctkGH), -1 /* no table */, "ML_ADM::ML_shiftrhs","flat");
+  if (ierr < 0)
+    CCTK_WARN(1, "Failed to register flat BC for ML_ADM::ML_shiftrhs.");
+  return;
+}
+
+static void ML_ADM_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[])
+{
+  DECLARE_CCTK_ARGUMENTS;
+  DECLARE_CCTK_PARAMETERS;
+  
+  
+  /* Declare finite differencing variables */
+  
+  if (verbose > 1)
+  {
+    CCTK_VInfo(CCTK_THORNSTRING,"Entering ML_ADM_RHS_Body");
+  }
+  
+  if (cctk_iteration % ML_ADM_RHS_calc_every != ML_ADM_RHS_calc_offset)
+  {
+    return;
+  }
+  
+  const char *groups[] = {"ML_ADM::ML_curv","ML_ADM::ML_curvrhs","ML_ADM::ML_lapse","ML_ADM::ML_lapserhs","ML_ADM::ML_metric","ML_ADM::ML_metricrhs","ML_ADM::ML_shift","ML_ADM::ML_shiftrhs"};
+  GenericFD_AssertGroupStorage(cctkGH, "ML_ADM_RHS", 8, 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);
+  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 dxi = INV(dx);
+  CCTK_REAL const dyi = INV(dy);
+  CCTK_REAL const dzi = INV(dz);
+  CCTK_REAL const khalf = 0.5;
+  CCTK_REAL const kthird = 1/3.0;
+  CCTK_REAL const ktwothird = 2.0/3.0;
+  CCTK_REAL const kfourthird = 4.0/3.0;
+  CCTK_REAL const keightthird = 8.0/3.0;
+  CCTK_REAL const hdxi = 0.5 * dxi;
+  CCTK_REAL const hdyi = 0.5 * dyi;
+  CCTK_REAL const hdzi = 0.5 * dzi;
+  
+  /* Initialize predefined quantities */
+  CCTK_REAL const p1o12dx = 0.0833333333333333333333333333333*INV(dx);
+  CCTK_REAL const p1o12dy = 0.0833333333333333333333333333333*INV(dy);
+  CCTK_REAL const p1o12dz = 0.0833333333333333333333333333333*INV(dz);
+  CCTK_REAL const p1o144dxdy = 0.00694444444444444444444444444444*INV(dx)*INV(dy);
+  CCTK_REAL const p1o144dxdz = 0.00694444444444444444444444444444*INV(dx)*INV(dz);
+  CCTK_REAL const p1o144dydz = 0.00694444444444444444444444444444*INV(dy)*INV(dz);
+  CCTK_REAL const pm1o12dx2 = -0.0833333333333333333333333333333*INV(SQR(dx));
+  CCTK_REAL const pm1o12dy2 = -0.0833333333333333333333333333333*INV(SQR(dy));
+  CCTK_REAL const pm1o12dz2 = -0.0833333333333333333333333333333*INV(SQR(dz));
+  
+  /* Loop over the grid points */
+  #pragma omp parallel
+  LC_LOOP3 (ML_ADM_RHS,
+    i,j,k, min[0],min[1],min[2], max[0],max[1],max[2],
+    cctk_lsh[0],cctk_lsh[1],cctk_lsh[2])
+  {
+    ptrdiff_t const index = di*i + dj*j + dk*k;
+    
+    /* Assign local copies of grid functions */
+    CCTK_REAL alphaL = alpha[index];
+    CCTK_REAL beta1L = beta1[index];
+    CCTK_REAL beta2L = beta2[index];
+    CCTK_REAL beta3L = beta3[index];
+    CCTK_REAL g11L = g11[index];
+    CCTK_REAL g12L = g12[index];
+    CCTK_REAL g13L = g13[index];
+    CCTK_REAL g22L = g22[index];
+    CCTK_REAL g23L = g23[index];
+    CCTK_REAL g33L = g33[index];
+    CCTK_REAL K11L = K11[index];
+    CCTK_REAL K12L = K12[index];
+    CCTK_REAL K13L = K13[index];
+    CCTK_REAL K22L = K22[index];
+    CCTK_REAL K23L = K23[index];
+    CCTK_REAL K33L = K33[index];
+    
+    /* Include user supplied include files */
+    
+    /* Precompute derivatives */
+    CCTK_REAL const PDstandardNth1alpha = PDstandardNth1(&alpha[index]);
+    CCTK_REAL const PDstandardNth2alpha = PDstandardNth2(&alpha[index]);
+    CCTK_REAL const PDstandardNth3alpha = PDstandardNth3(&alpha[index]);
+    CCTK_REAL const PDstandardNth11alpha = PDstandardNth11(&alpha[index]);
+    CCTK_REAL const PDstandardNth22alpha = PDstandardNth22(&alpha[index]);
+    CCTK_REAL const PDstandardNth33alpha = PDstandardNth33(&alpha[index]);
+    CCTK_REAL const PDstandardNth12alpha = PDstandardNth12(&alpha[index]);
+    CCTK_REAL const PDstandardNth13alpha = PDstandardNth13(&alpha[index]);
+    CCTK_REAL const PDstandardNth23alpha = PDstandardNth23(&alpha[index]);
+    CCTK_REAL const PDstandardNth1beta1 = PDstandardNth1(&beta1[index]);
+    CCTK_REAL const PDstandardNth2beta1 = PDstandardNth2(&beta1[index]);
+    CCTK_REAL const PDstandardNth3beta1 = PDstandardNth3(&beta1[index]);
+    CCTK_REAL const PDstandardNth1beta2 = PDstandardNth1(&beta2[index]);
+    CCTK_REAL const PDstandardNth2beta2 = PDstandardNth2(&beta2[index]);
+    CCTK_REAL const PDstandardNth3beta2 = PDstandardNth3(&beta2[index]);
+    CCTK_REAL const PDstandardNth1beta3 = PDstandardNth1(&beta3[index]);
+    CCTK_REAL const PDstandardNth2beta3 = PDstandardNth2(&beta3[index]);
+    CCTK_REAL const PDstandardNth3beta3 = PDstandardNth3(&beta3[index]);
+    CCTK_REAL const PDstandardNth1g11 = PDstandardNth1(&g11[index]);
+    CCTK_REAL const PDstandardNth2g11 = PDstandardNth2(&g11[index]);
+    CCTK_REAL const PDstandardNth3g11 = PDstandardNth3(&g11[index]);
+    CCTK_REAL const PDstandardNth22g11 = PDstandardNth22(&g11[index]);
+    CCTK_REAL const PDstandardNth33g11 = PDstandardNth33(&g11[index]);
+    CCTK_REAL const PDstandardNth23g11 = PDstandardNth23(&g11[index]);
+    CCTK_REAL const PDstandardNth1g12 = PDstandardNth1(&g12[index]);
+    CCTK_REAL const PDstandardNth2g12 = PDstandardNth2(&g12[index]);
+    CCTK_REAL const PDstandardNth3g12 = PDstandardNth3(&g12[index]);
+    CCTK_REAL const PDstandardNth33g12 = PDstandardNth33(&g12[index]);
+    CCTK_REAL const PDstandardNth12g12 = PDstandardNth12(&g12[index]);
+    CCTK_REAL const PDstandardNth13g12 = PDstandardNth13(&g12[index]);
+    CCTK_REAL const PDstandardNth23g12 = PDstandardNth23(&g12[index]);
+    CCTK_REAL const PDstandardNth1g13 = PDstandardNth1(&g13[index]);
+    CCTK_REAL const PDstandardNth2g13 = PDstandardNth2(&g13[index]);
+    CCTK_REAL const PDstandardNth3g13 = PDstandardNth3(&g13[index]);
+    CCTK_REAL const PDstandardNth22g13 = PDstandardNth22(&g13[index]);
+    CCTK_REAL const PDstandardNth12g13 = PDstandardNth12(&g13[index]);
+    CCTK_REAL const PDstandardNth13g13 = PDstandardNth13(&g13[index]);
+    CCTK_REAL const PDstandardNth23g13 = PDstandardNth23(&g13[index]);
+    CCTK_REAL const PDstandardNth1g22 = PDstandardNth1(&g22[index]);
+    CCTK_REAL const PDstandardNth2g22 = PDstandardNth2(&g22[index]);
+    CCTK_REAL const PDstandardNth3g22 = PDstandardNth3(&g22[index]);
+    CCTK_REAL const PDstandardNth11g22 = PDstandardNth11(&g22[index]);
+    CCTK_REAL const PDstandardNth33g22 = PDstandardNth33(&g22[index]);
+    CCTK_REAL const PDstandardNth13g22 = PDstandardNth13(&g22[index]);
+    CCTK_REAL const PDstandardNth1g23 = PDstandardNth1(&g23[index]);
+    CCTK_REAL const PDstandardNth2g23 = PDstandardNth2(&g23[index]);
+    CCTK_REAL const PDstandardNth3g23 = PDstandardNth3(&g23[index]);
+    CCTK_REAL const PDstandardNth11g23 = PDstandardNth11(&g23[index]);
+    CCTK_REAL const PDstandardNth12g23 = PDstandardNth12(&g23[index]);
+    CCTK_REAL const PDstandardNth13g23 = PDstandardNth13(&g23[index]);
+    CCTK_REAL const PDstandardNth23g23 = PDstandardNth23(&g23[index]);
+    CCTK_REAL const PDstandardNth1g33 = PDstandardNth1(&g33[index]);
+    CCTK_REAL const PDstandardNth2g33 = PDstandardNth2(&g33[index]);
+    CCTK_REAL const PDstandardNth3g33 = PDstandardNth3(&g33[index]);
+    CCTK_REAL const PDstandardNth11g33 = PDstandardNth11(&g33[index]);
+    CCTK_REAL const PDstandardNth22g33 = PDstandardNth22(&g33[index]);
+    CCTK_REAL const PDstandardNth12g33 = PDstandardNth12(&g33[index]);
+    CCTK_REAL const PDstandardNth1K11 = PDstandardNth1(&K11[index]);
+    CCTK_REAL const PDstandardNth2K11 = PDstandardNth2(&K11[index]);
+    CCTK_REAL const PDstandardNth3K11 = PDstandardNth3(&K11[index]);
+    CCTK_REAL const PDstandardNth1K12 = PDstandardNth1(&K12[index]);
+    CCTK_REAL const PDstandardNth2K12 = PDstandardNth2(&K12[index]);
+    CCTK_REAL const PDstandardNth3K12 = PDstandardNth3(&K12[index]);
+    CCTK_REAL const PDstandardNth1K13 = PDstandardNth1(&K13[index]);
+    CCTK_REAL const PDstandardNth2K13 = PDstandardNth2(&K13[index]);
+    CCTK_REAL const PDstandardNth3K13 = PDstandardNth3(&K13[index]);
+    CCTK_REAL const PDstandardNth1K22 = PDstandardNth1(&K22[index]);
+    CCTK_REAL const PDstandardNth2K22 = PDstandardNth2(&K22[index]);
+    CCTK_REAL const PDstandardNth3K22 = PDstandardNth3(&K22[index]);
+    CCTK_REAL const PDstandardNth1K23 = PDstandardNth1(&K23[index]);
+    CCTK_REAL const PDstandardNth2K23 = PDstandardNth2(&K23[index]);
+    CCTK_REAL const PDstandardNth3K23 = PDstandardNth3(&K23[index]);
+    CCTK_REAL const PDstandardNth1K33 = PDstandardNth1(&K33[index]);
+    CCTK_REAL const PDstandardNth2K33 = PDstandardNth2(&K33[index]);
+    CCTK_REAL const PDstandardNth3K33 = PDstandardNth3(&K33[index]);
+    
+    /* Calculate temporaries and grid functions */
+    CCTK_REAL detg = 2*g12L*g13L*g23L + g33L*(g11L*g22L - SQR(g12L)) - 
+      g22L*SQR(g13L) - g11L*SQR(g23L);
+    
+    CCTK_REAL gu11 = INV(detg)*(g22L*g33L - SQR(g23L));
+    
+    CCTK_REAL gu12 = (g13L*g23L - g12L*g33L)*INV(detg);
+    
+    CCTK_REAL gu13 = (-(g13L*g22L) + g12L*g23L)*INV(detg);
+    
+    CCTK_REAL gu21 = (g13L*g23L - g12L*g33L)*INV(detg);
+    
+    CCTK_REAL gu22 = INV(detg)*(g11L*g33L - SQR(g13L));
+    
+    CCTK_REAL gu23 = (g12L*g13L - g11L*g23L)*INV(detg);
+    
+    CCTK_REAL gu31 = (-(g13L*g22L) + g12L*g23L)*INV(detg);
+    
+    CCTK_REAL gu32 = (g12L*g13L - g11L*g23L)*INV(detg);
+    
+    CCTK_REAL gu33 = INV(detg)*(g11L*g22L - SQR(g12L));
+    
+    CCTK_REAL G111 = 0.5*(gu11*PDstandardNth1g11 + 
+      2*(gu12*PDstandardNth1g12 + gu13*PDstandardNth1g13) - 
+      gu12*PDstandardNth2g11 - gu13*PDstandardNth3g11);
+    
+    CCTK_REAL G211 = 0.5*(gu21*PDstandardNth1g11 + 
+      2*(gu22*PDstandardNth1g12 + gu23*PDstandardNth1g13) - 
+      gu22*PDstandardNth2g11 - gu23*PDstandardNth3g11);
+    
+    CCTK_REAL G311 = 0.5*(gu31*PDstandardNth1g11 + 
+      2*(gu32*PDstandardNth1g12 + gu33*PDstandardNth1g13) - 
+      gu32*PDstandardNth2g11 - gu33*PDstandardNth3g11);
+    
+    CCTK_REAL G112 = 0.5*(gu12*PDstandardNth1g22 + gu11*PDstandardNth2g11 
+      + gu13*(PDstandardNth1g23 + PDstandardNth2g13 - PDstandardNth3g12));
+    
+    CCTK_REAL G212 = 0.5*(gu22*PDstandardNth1g22 + gu21*PDstandardNth2g11 
+      + gu23*(PDstandardNth1g23 + PDstandardNth2g13 - PDstandardNth3g12));
+    
+    CCTK_REAL G312 = 0.5*(gu32*PDstandardNth1g22 + gu31*PDstandardNth2g11 
+      + gu33*(PDstandardNth1g23 + PDstandardNth2g13 - PDstandardNth3g12));
+    
+    CCTK_REAL G113 = 0.5*(gu13*PDstandardNth1g33 + gu11*PDstandardNth3g11 
+      + gu12*(PDstandardNth1g23 - PDstandardNth2g13 + PDstandardNth3g12));
+    
+    CCTK_REAL G213 = 0.5*(gu23*PDstandardNth1g33 + gu21*PDstandardNth3g11 
+      + gu22*(PDstandardNth1g23 - PDstandardNth2g13 + PDstandardNth3g12));
+    
+    CCTK_REAL G313 = 0.5*(gu33*PDstandardNth1g33 + gu31*PDstandardNth3g11 
+      + gu32*(PDstandardNth1g23 - PDstandardNth2g13 + PDstandardNth3g12));
+    
+    CCTK_REAL G122 = 0.5*(gu11*(-PDstandardNth1g22 + 2*PDstandardNth2g12) 
+      + gu12*PDstandardNth2g22 + gu13*(2*PDstandardNth2g23 - 
+      PDstandardNth3g22));
+    
+    CCTK_REAL G222 = 0.5*(gu21*(-PDstandardNth1g22 + 2*PDstandardNth2g12) 
+      + gu22*PDstandardNth2g22 + gu23*(2*PDstandardNth2g23 - 
+      PDstandardNth3g22));
+    
+    CCTK_REAL G322 = 0.5*(gu31*(-PDstandardNth1g22 + 2*PDstandardNth2g12) 
+      + gu32*PDstandardNth2g22 + gu33*(2*PDstandardNth2g23 - 
+      PDstandardNth3g22));
+    
+    CCTK_REAL G123 = 0.5*(gu13*PDstandardNth2g33 + 
+      gu11*(-PDstandardNth1g23 + PDstandardNth2g13 + PDstandardNth3g12) + 
+      gu12*PDstandardNth3g22);
+    
+    CCTK_REAL G223 = 0.5*(gu23*PDstandardNth2g33 + 
+      gu21*(-PDstandardNth1g23 + PDstandardNth2g13 + PDstandardNth3g12) + 
+      gu22*PDstandardNth3g22);
+    
+    CCTK_REAL G323 = 0.5*(gu33*PDstandardNth2g33 + 
+      gu31*(-PDstandardNth1g23 + PDstandardNth2g13 + PDstandardNth3g12) + 
+      gu32*PDstandardNth3g22);
+    
+    CCTK_REAL G133 = 0.5*(gu11*(-PDstandardNth1g33 + 2*PDstandardNth3g13) 
+      + gu12*(-PDstandardNth2g33 + 2*PDstandardNth3g23) + 
+      gu13*PDstandardNth3g33);
+    
+    CCTK_REAL G233 = 0.5*(gu21*(-PDstandardNth1g33 + 2*PDstandardNth3g13) 
+      + gu22*(-PDstandardNth2g33 + 2*PDstandardNth3g23) + 
+      gu23*PDstandardNth3g33);
+    
+    CCTK_REAL G333 = 0.5*(gu31*(-PDstandardNth1g33 + 2*PDstandardNth3g13) 
+      + gu32*(-PDstandardNth2g33 + 2*PDstandardNth3g23) + 
+      gu33*PDstandardNth3g33);
+    
+    CCTK_REAL R11 = 2*(G112*G211 + G113*G311 + G213*G312) - G111*(G111 + 
+      G212 + G313) - G211*(G112 + G222 + G323) - G311*(G113 + G223 + G333) + 
+      0.5*(-(gu22*(PDstandardNth11g22 - 2*PDstandardNth12g12 + 
+      PDstandardNth22g11)) + gu23*(-PDstandardNth11g23 + PDstandardNth12g13 + 
+      PDstandardNth13g12 - PDstandardNth23g11) + gu32*(-PDstandardNth11g23 + 
+      PDstandardNth12g13 + PDstandardNth13g12 - PDstandardNth23g11) - 
+      gu33*(PDstandardNth11g33 - 2*PDstandardNth13g13 + PDstandardNth33g11)) 
+      + SQR(G111) + SQR(G212) + SQR(G313);
+    
+    CCTK_REAL R12 = 0.5*(2*(G122*G211 + G123*G311 + G213*G322 + G313*G323) 
+      - 2*(G112*G212 + G112*G313 + G212*G323 + G312*G333 + 
+      gu12*PDstandardNth12g12) + gu13*(PDstandardNth11g23 - 
+      PDstandardNth12g13 - PDstandardNth13g12) + gu12*(PDstandardNth11g22 + 
+      PDstandardNth22g11) + gu32*PDstandardNth22g13 + gu13*PDstandardNth23g11 
+      + gu32*(-PDstandardNth12g23 + PDstandardNth13g22 - PDstandardNth23g12) 
+      + gu33*PDstandardNth23g13 + gu33*(-PDstandardNth12g33 + 
+      PDstandardNth13g23 - PDstandardNth33g12));
+    
+    CCTK_REAL R13 = 0.5*(2*(G123*G211 + G212*G223 + G133*G311 + G233*G312) 
+      - 2*(G213*G222 + G223*G313 + G113*(G212 + G313) + 
+      gu13*PDstandardNth13g13) + gu12*(PDstandardNth11g23 - 
+      PDstandardNth12g13 - PDstandardNth13g12 + PDstandardNth23g11) + 
+      gu22*(PDstandardNth12g23 - PDstandardNth13g22 - PDstandardNth22g13 + 
+      PDstandardNth23g12) + gu13*(PDstandardNth11g33 + PDstandardNth33g11) + 
+      gu23*(PDstandardNth12g33 - PDstandardNth13g23 - PDstandardNth23g13 + 
+      PDstandardNth33g12));
+    
+    CCTK_REAL R22 = -(G122*(G111 + G212 + G313)) + 2*(G122*G212 + 
+      G123*G312 + G223*G322) - G222*(G112 + G222 + G323) - G322*(G113 + G223 
+      + G333) + 0.5*(-(gu11*(PDstandardNth11g22 - 2*PDstandardNth12g12 + 
+      PDstandardNth22g11)) + gu13*(PDstandardNth12g23 - PDstandardNth13g22 - 
+      PDstandardNth22g13 + PDstandardNth23g12) + gu31*(PDstandardNth12g23 - 
+      PDstandardNth13g22 - PDstandardNth22g13 + PDstandardNth23g12) - 
+      gu33*(PDstandardNth22g33 - 2*PDstandardNth23g23 + PDstandardNth33g22)) 
+      + SQR(G112) + SQR(G222) + SQR(G323);
+    
+    CCTK_REAL R23 = 0.5*(2*(G112*G113 + G122*G213 + G133*G312 + G233*G322) 
+      + gu11*(-PDstandardNth11g23 + PDstandardNth12g13 + PDstandardNth13g12 - 
+      PDstandardNth23g11) + gu21*(-PDstandardNth12g23 + PDstandardNth13g22 + 
+      PDstandardNth22g13 - PDstandardNth23g12) - 2*(G111*G123 + G113*G323 + 
+      G223*(G112 + G323) + gu23*PDstandardNth23g23) + 
+      gu13*(PDstandardNth12g33 - PDstandardNth13g23 - PDstandardNth23g13 + 
+      PDstandardNth33g12) + gu23*(PDstandardNth22g33 + PDstandardNth33g22));
+    
+    CCTK_REAL R33 = -(G133*(G111 + G212 + G313)) + 2*(G123*G213 + 
+      G133*G313) + G233*(-G112 - G222 + G323) - G333*(G113 + G223 + G333) + 
+      0.5*(-(gu11*(PDstandardNth11g33 - 2*PDstandardNth13g13 + 
+      PDstandardNth33g11)) + gu12*(-PDstandardNth12g33 + PDstandardNth13g23 + 
+      PDstandardNth23g13 - PDstandardNth33g12) + gu21*(-PDstandardNth12g33 + 
+      PDstandardNth13g23 + PDstandardNth23g13 - PDstandardNth33g12) - 
+      gu22*(PDstandardNth22g33 - 2*PDstandardNth23g23 + PDstandardNth33g22)) 
+      + SQR(G113) + SQR(G223) + SQR(G333);
+    
+    CCTK_REAL Km11 = gu11*K11L + gu12*K12L + gu13*K13L;
+    
+    CCTK_REAL Km21 = gu21*K11L + gu22*K12L + gu23*K13L;
+    
+    CCTK_REAL Km31 = gu31*K11L + gu32*K12L + gu33*K13L;
+    
+    CCTK_REAL Km12 = gu11*K12L + gu12*K22L + gu13*K23L;
+    
+    CCTK_REAL Km22 = gu21*K12L + gu22*K22L + gu23*K23L;
+    
+    CCTK_REAL Km32 = gu31*K12L + gu32*K22L + gu33*K23L;
+    
+    CCTK_REAL Km13 = gu11*K13L + gu12*K23L + gu13*K33L;
+    
+    CCTK_REAL Km23 = gu21*K13L + gu22*K23L + gu23*K33L;
+    
+    CCTK_REAL Km33 = gu31*K13L + gu32*K23L + gu33*K33L;
+    
+    CCTK_REAL trK = Km11 + Km22 + Km33;
+    
+    CCTK_REAL g11rhsL = -2*alphaL*K11L + 2*(g11L*PDstandardNth1beta1 + 
+      g12L*PDstandardNth1beta2 + g13L*PDstandardNth1beta3) + 
+      beta1L*PDstandardNth1g11 + beta2L*PDstandardNth2g11 + 
+      beta3L*PDstandardNth3g11;
+    
+    CCTK_REAL g12rhsL = -2*alphaL*K12L + g22L*PDstandardNth1beta2 + 
+      g23L*PDstandardNth1beta3 + beta1L*PDstandardNth1g12 + 
+      g11L*PDstandardNth2beta1 + g12L*(PDstandardNth1beta1 + 
+      PDstandardNth2beta2) + g13L*PDstandardNth2beta3 + 
+      beta2L*PDstandardNth2g12 + beta3L*PDstandardNth3g12;
+    
+    CCTK_REAL g13rhsL = -2*alphaL*K13L + g23L*PDstandardNth1beta2 + 
+      g33L*PDstandardNth1beta3 + beta1L*PDstandardNth1g13 + 
+      beta2L*PDstandardNth2g13 + g11L*PDstandardNth3beta1 + 
+      g12L*PDstandardNth3beta2 + g13L*(PDstandardNth1beta1 + 
+      PDstandardNth3beta3) + beta3L*PDstandardNth3g13;
+    
+    CCTK_REAL g22rhsL = -2*alphaL*K22L + beta1L*PDstandardNth1g22 + 
+      2*(g12L*PDstandardNth2beta1 + g22L*PDstandardNth2beta2 + 
+      g23L*PDstandardNth2beta3) + beta2L*PDstandardNth2g22 + 
+      beta3L*PDstandardNth3g22;
+    
+    CCTK_REAL g23rhsL = -2*alphaL*K23L + beta1L*PDstandardNth1g23 + 
+      g13L*PDstandardNth2beta1 + g33L*PDstandardNth2beta3 + 
+      beta2L*PDstandardNth2g23 + g12L*PDstandardNth3beta1 + 
+      g22L*PDstandardNth3beta2 + g23L*(PDstandardNth2beta2 + 
+      PDstandardNth3beta3) + beta3L*PDstandardNth3g23;
+    
+    CCTK_REAL g33rhsL = -2*alphaL*K33L + beta1L*PDstandardNth1g33 + 
+      beta2L*PDstandardNth2g33 + 2*(g13L*PDstandardNth3beta1 + 
+      g23L*PDstandardNth3beta2 + g33L*PDstandardNth3beta3) + 
+      beta3L*PDstandardNth3g33;
+    
+    CCTK_REAL K11rhsL = -PDstandardNth11alpha + G111*PDstandardNth1alpha + 
+      2*(K11L*PDstandardNth1beta1 + K12L*PDstandardNth1beta2 + 
+      K13L*PDstandardNth1beta3) + beta1L*PDstandardNth1K11 + 
+      G211*PDstandardNth2alpha + beta2L*PDstandardNth2K11 + 
+      G311*PDstandardNth3alpha + beta3L*PDstandardNth3K11 + 
+      alphaL*(-2*(K11L*Km11 + K12L*Km21 + K13L*Km31) + R11 + K11L*trK);
+    
+    CCTK_REAL K12rhsL = -PDstandardNth12alpha + G112*PDstandardNth1alpha + 
+      K22L*PDstandardNth1beta2 + K23L*PDstandardNth1beta3 + 
+      beta1L*PDstandardNth1K12 + G212*PDstandardNth2alpha + 
+      K11L*PDstandardNth2beta1 + K12L*(PDstandardNth1beta1 + 
+      PDstandardNth2beta2) + K13L*PDstandardNth2beta3 + 
+      beta2L*PDstandardNth2K12 + G312*PDstandardNth3alpha + 
+      beta3L*PDstandardNth3K12 + alphaL*(-2*(K11L*Km12 + K12L*Km22 + 
+      K13L*Km32) + R12 + K12L*trK);
+    
+    CCTK_REAL K13rhsL = -PDstandardNth13alpha + G113*PDstandardNth1alpha + 
+      K23L*PDstandardNth1beta2 + K33L*PDstandardNth1beta3 + 
+      beta1L*PDstandardNth1K13 + G213*PDstandardNth2alpha + 
+      beta2L*PDstandardNth2K13 + G313*PDstandardNth3alpha + 
+      K11L*PDstandardNth3beta1 + K12L*PDstandardNth3beta2 + 
+      K13L*(PDstandardNth1beta1 + PDstandardNth3beta3) + 
+      beta3L*PDstandardNth3K13 + alphaL*(-2*(K11L*Km13 + K12L*Km23 + 
+      K13L*Km33) + R13 + K13L*trK);
+    
+    CCTK_REAL K22rhsL = G122*PDstandardNth1alpha + 
+      beta1L*PDstandardNth1K22 - PDstandardNth22alpha + 
+      G222*PDstandardNth2alpha + 2*(K12L*PDstandardNth2beta1 + 
+      K22L*PDstandardNth2beta2 + K23L*PDstandardNth2beta3) + 
+      beta2L*PDstandardNth2K22 + G322*PDstandardNth3alpha + 
+      beta3L*PDstandardNth3K22 + alphaL*(-2*(K12L*Km12 + K22L*Km22 + 
+      K23L*Km32) + R22 + K22L*trK);
+    
+    CCTK_REAL K23rhsL = G123*PDstandardNth1alpha + 
+      beta1L*PDstandardNth1K23 - PDstandardNth23alpha + 
+      G223*PDstandardNth2alpha + K13L*PDstandardNth2beta1 + 
+      K33L*PDstandardNth2beta3 + beta2L*PDstandardNth2K23 + 
+      G323*PDstandardNth3alpha + K12L*PDstandardNth3beta1 + 
+      K22L*PDstandardNth3beta2 + K23L*(PDstandardNth2beta2 + 
+      PDstandardNth3beta3) + beta3L*PDstandardNth3K23 + alphaL*(-2*(K12L*Km13 
+      + K22L*Km23 + K23L*Km33) + R23 + K23L*trK);
+    
+    CCTK_REAL K33rhsL = G133*PDstandardNth1alpha + 
+      beta1L*PDstandardNth1K33 + G233*PDstandardNth2alpha + 
+      beta2L*PDstandardNth2K33 - PDstandardNth33alpha + 
+      G333*PDstandardNth3alpha + 2*(K13L*PDstandardNth3beta1 + 
+      K23L*PDstandardNth3beta2 + K33L*PDstandardNth3beta3) + 
+      beta3L*PDstandardNth3K33 + alphaL*(-2*(K13L*Km13 + K23L*Km23 + 
+      K33L*Km33) + R33 + K33L*trK);
+    
+    CCTK_REAL alpharhsL = 0;
+    
+    CCTK_REAL beta1rhsL = 0;
+    
+    CCTK_REAL beta2rhsL = 0;
+    
+    CCTK_REAL beta3rhsL = 0;
+    
+    
+    /* Copy local copies back to grid functions */
+    alpharhs[index] = alpharhsL;
+    beta1rhs[index] = beta1rhsL;
+    beta2rhs[index] = beta2rhsL;
+    beta3rhs[index] = beta3rhsL;
+    g11rhs[index] = g11rhsL;
+    g12rhs[index] = g12rhsL;
+    g13rhs[index] = g13rhsL;
+    g22rhs[index] = g22rhsL;
+    g23rhs[index] = g23rhsL;
+    g33rhs[index] = g33rhsL;
+    K11rhs[index] = K11rhsL;
+    K12rhs[index] = K12rhsL;
+    K13rhs[index] = K13rhsL;
+    K22rhs[index] = K22rhsL;
+    K23rhs[index] = K23rhsL;
+    K33rhs[index] = K33rhsL;
+  }
+  LC_ENDLOOP3 (ML_ADM_RHS);
+}
+
+extern "C" void ML_ADM_RHS(CCTK_ARGUMENTS)
+{
+  DECLARE_CCTK_ARGUMENTS;
+  DECLARE_CCTK_PARAMETERS;
+  
+  GenericFD_LoopOverInterior(cctkGH, &ML_ADM_RHS_Body);
+}
diff --git a/ML_ADM/src/ML_ADM_boundary.cc b/ML_ADM/src/ML_ADM_boundary.cc
new file mode 100644
index 0000000..790928b
--- /dev/null
+++ b/ML_ADM/src/ML_ADM_boundary.cc
@@ -0,0 +1,172 @@
+/*  File produced by Kranc */
+
+#define KRANC_C
+
+#include <assert.h>
+#include <math.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include "cctk.h"
+#include "cctk_Arguments.h"
+#include "cctk_Parameters.h"
+#include "GenericFD.h"
+#include "Differencing.h"
+#include "loopcontrol.h"
+
+/* Define macros used in calculations */
+#define INITVALUE (42)
+#define QAD(x) (SQR(SQR(x)))
+#define INV(x) ((1.0) / (x))
+#define SQR(x) ((x) * (x))
+#define CUB(x) ((x) * (x) * (x))
+
+extern "C" void ML_ADM_boundary_SelectBCs(CCTK_ARGUMENTS)
+{
+  DECLARE_CCTK_ARGUMENTS;
+  DECLARE_CCTK_PARAMETERS;
+  
+  CCTK_INT ierr = 0;
+  ierr = Boundary_SelectGroupForBC(cctkGH, CCTK_ALL_FACES, GenericFD_GetBoundaryWidth(cctkGH), -1 /* no table */, "ML_ADM::ML_curv","flat");
+  if (ierr < 0)
+    CCTK_WARN(1, "Failed to register flat BC for ML_ADM::ML_curv.");
+  ierr = Boundary_SelectGroupForBC(cctkGH, CCTK_ALL_FACES, GenericFD_GetBoundaryWidth(cctkGH), -1 /* no table */, "ML_ADM::ML_lapse","flat");
+  if (ierr < 0)
+    CCTK_WARN(1, "Failed to register flat BC for ML_ADM::ML_lapse.");
+  ierr = Boundary_SelectGroupForBC(cctkGH, CCTK_ALL_FACES, GenericFD_GetBoundaryWidth(cctkGH), -1 /* no table */, "ML_ADM::ML_metric","flat");
+  if (ierr < 0)
+    CCTK_WARN(1, "Failed to register flat BC for ML_ADM::ML_metric.");
+  ierr = Boundary_SelectGroupForBC(cctkGH, CCTK_ALL_FACES, GenericFD_GetBoundaryWidth(cctkGH), -1 /* no table */, "ML_ADM::ML_shift","flat");
+  if (ierr < 0)
+    CCTK_WARN(1, "Failed to register flat BC for ML_ADM::ML_shift.");
+  return;
+}
+
+static void ML_ADM_boundary_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[])
+{
+  DECLARE_CCTK_ARGUMENTS;
+  DECLARE_CCTK_PARAMETERS;
+  
+  
+  /* Declare finite differencing variables */
+  
+  if (verbose > 1)
+  {
+    CCTK_VInfo(CCTK_THORNSTRING,"Entering ML_ADM_boundary_Body");
+  }
+  
+  if (cctk_iteration % ML_ADM_boundary_calc_every != ML_ADM_boundary_calc_offset)
+  {
+    return;
+  }
+  
+  const char *groups[] = {"ML_ADM::ML_curv","ML_ADM::ML_lapse","ML_ADM::ML_metric","ML_ADM::ML_shift"};
+  GenericFD_AssertGroupStorage(cctkGH, "ML_ADM_boundary", 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);
+  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 dxi = INV(dx);
+  CCTK_REAL const dyi = INV(dy);
+  CCTK_REAL const dzi = INV(dz);
+  CCTK_REAL const khalf = 0.5;
+  CCTK_REAL const kthird = 1/3.0;
+  CCTK_REAL const ktwothird = 2.0/3.0;
+  CCTK_REAL const kfourthird = 4.0/3.0;
+  CCTK_REAL const keightthird = 8.0/3.0;
+  CCTK_REAL const hdxi = 0.5 * dxi;
+  CCTK_REAL const hdyi = 0.5 * dyi;
+  CCTK_REAL const hdzi = 0.5 * dzi;
+  
+  /* Initialize predefined quantities */
+  CCTK_REAL const p1o12dx = 0.0833333333333333333333333333333*INV(dx);
+  CCTK_REAL const p1o12dy = 0.0833333333333333333333333333333*INV(dy);
+  CCTK_REAL const p1o12dz = 0.0833333333333333333333333333333*INV(dz);
+  CCTK_REAL const p1o144dxdy = 0.00694444444444444444444444444444*INV(dx)*INV(dy);
+  CCTK_REAL const p1o144dxdz = 0.00694444444444444444444444444444*INV(dx)*INV(dz);
+  CCTK_REAL const p1o144dydz = 0.00694444444444444444444444444444*INV(dy)*INV(dz);
+  CCTK_REAL const pm1o12dx2 = -0.0833333333333333333333333333333*INV(SQR(dx));
+  CCTK_REAL const pm1o12dy2 = -0.0833333333333333333333333333333*INV(SQR(dy));
+  CCTK_REAL const pm1o12dz2 = -0.0833333333333333333333333333333*INV(SQR(dz));
+  
+  /* Loop over the grid points */
+  #pragma omp parallel
+  LC_LOOP3 (ML_ADM_boundary,
+    i,j,k, min[0],min[1],min[2], max[0],max[1],max[2],
+    cctk_lsh[0],cctk_lsh[1],cctk_lsh[2])
+  {
+    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 g11L = 1;
+    
+    CCTK_REAL g12L = 0;
+    
+    CCTK_REAL g13L = 0;
+    
+    CCTK_REAL g22L = 1;
+    
+    CCTK_REAL g23L = 0;
+    
+    CCTK_REAL g33L = 1;
+    
+    CCTK_REAL K11L = 0;
+    
+    CCTK_REAL K12L = 0;
+    
+    CCTK_REAL K13L = 0;
+    
+    CCTK_REAL K22L = 0;
+    
+    CCTK_REAL K23L = 0;
+    
+    CCTK_REAL K33L = 0;
+    
+    CCTK_REAL alphaL = 1;
+    
+    CCTK_REAL beta1L = 0;
+    
+    CCTK_REAL beta2L = 0;
+    
+    CCTK_REAL beta3L = 0;
+    
+    
+    /* Copy local copies back to grid functions */
+    alpha[index] = alphaL;
+    beta1[index] = beta1L;
+    beta2[index] = beta2L;
+    beta3[index] = beta3L;
+    g11[index] = g11L;
+    g12[index] = g12L;
+    g13[index] = g13L;
+    g22[index] = g22L;
+    g23[index] = g23L;
+    g33[index] = g33L;
+    K11[index] = K11L;
+    K12[index] = K12L;
+    K13[index] = K13L;
+    K22[index] = K22L;
+    K23[index] = K23L;
+    K33[index] = K33L;
+  }
+  LC_ENDLOOP3 (ML_ADM_boundary);
+}
+
+extern "C" void ML_ADM_boundary(CCTK_ARGUMENTS)
+{
+  DECLARE_CCTK_ARGUMENTS;
+  DECLARE_CCTK_PARAMETERS;
+  
+  GenericFD_LoopOverBoundaryWithGhosts(cctkGH, &ML_ADM_boundary_Body);
+}
diff --git a/ML_ADM/src/ML_ADM_constraints.cc b/ML_ADM/src/ML_ADM_constraints.cc
new file mode 100644
index 0000000..132f6e6
--- /dev/null
+++ b/ML_ADM/src/ML_ADM_constraints.cc
@@ -0,0 +1,390 @@
+/*  File produced by Kranc */
+
+#define KRANC_C
+
+#include <assert.h>
+#include <math.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include "cctk.h"
+#include "cctk_Arguments.h"
+#include "cctk_Parameters.h"
+#include "GenericFD.h"
+#include "Differencing.h"
+#include "loopcontrol.h"
+
+/* Define macros used in calculations */
+#define INITVALUE (42)
+#define QAD(x) (SQR(SQR(x)))
+#define INV(x) ((1.0) / (x))
+#define SQR(x) ((x) * (x))
+#define CUB(x) ((x) * (x) * (x))
+
+extern "C" void ML_ADM_constraints_SelectBCs(CCTK_ARGUMENTS)
+{
+  DECLARE_CCTK_ARGUMENTS;
+  DECLARE_CCTK_PARAMETERS;
+  
+  CCTK_INT ierr = 0;
+  ierr = Boundary_SelectGroupForBC(cctkGH, CCTK_ALL_FACES, GenericFD_GetBoundaryWidth(cctkGH), -1 /* no table */, "ML_ADM::ML_Ham","flat");
+  if (ierr < 0)
+    CCTK_WARN(1, "Failed to register flat BC for ML_ADM::ML_Ham.");
+  ierr = Boundary_SelectGroupForBC(cctkGH, CCTK_ALL_FACES, GenericFD_GetBoundaryWidth(cctkGH), -1 /* no table */, "ML_ADM::ML_mom","flat");
+  if (ierr < 0)
+    CCTK_WARN(1, "Failed to register flat BC for ML_ADM::ML_mom.");
+  return;
+}
+
+static void ML_ADM_constraints_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[])
+{
+  DECLARE_CCTK_ARGUMENTS;
+  DECLARE_CCTK_PARAMETERS;
+  
+  
+  /* Declare finite differencing variables */
+  
+  if (verbose > 1)
+  {
+    CCTK_VInfo(CCTK_THORNSTRING,"Entering ML_ADM_constraints_Body");
+  }
+  
+  if (cctk_iteration % ML_ADM_constraints_calc_every != ML_ADM_constraints_calc_offset)
+  {
+    return;
+  }
+  
+  const char *groups[] = {"ML_ADM::ML_curv","ML_ADM::ML_Ham","ML_ADM::ML_metric","ML_ADM::ML_mom"};
+  GenericFD_AssertGroupStorage(cctkGH, "ML_ADM_constraints", 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);
+  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 dxi = INV(dx);
+  CCTK_REAL const dyi = INV(dy);
+  CCTK_REAL const dzi = INV(dz);
+  CCTK_REAL const khalf = 0.5;
+  CCTK_REAL const kthird = 1/3.0;
+  CCTK_REAL const ktwothird = 2.0/3.0;
+  CCTK_REAL const kfourthird = 4.0/3.0;
+  CCTK_REAL const keightthird = 8.0/3.0;
+  CCTK_REAL const hdxi = 0.5 * dxi;
+  CCTK_REAL const hdyi = 0.5 * dyi;
+  CCTK_REAL const hdzi = 0.5 * dzi;
+  
+  /* Initialize predefined quantities */
+  CCTK_REAL const p1o12dx = 0.0833333333333333333333333333333*INV(dx);
+  CCTK_REAL const p1o12dy = 0.0833333333333333333333333333333*INV(dy);
+  CCTK_REAL const p1o12dz = 0.0833333333333333333333333333333*INV(dz);
+  CCTK_REAL const p1o144dxdy = 0.00694444444444444444444444444444*INV(dx)*INV(dy);
+  CCTK_REAL const p1o144dxdz = 0.00694444444444444444444444444444*INV(dx)*INV(dz);
+  CCTK_REAL const p1o144dydz = 0.00694444444444444444444444444444*INV(dy)*INV(dz);
+  CCTK_REAL const pm1o12dx2 = -0.0833333333333333333333333333333*INV(SQR(dx));
+  CCTK_REAL const pm1o12dy2 = -0.0833333333333333333333333333333*INV(SQR(dy));
+  CCTK_REAL const pm1o12dz2 = -0.0833333333333333333333333333333*INV(SQR(dz));
+  
+  /* Loop over the grid points */
+  #pragma omp parallel
+  LC_LOOP3 (ML_ADM_constraints,
+    i,j,k, min[0],min[1],min[2], max[0],max[1],max[2],
+    cctk_lsh[0],cctk_lsh[1],cctk_lsh[2])
+  {
+    ptrdiff_t const index = di*i + dj*j + dk*k;
+    
+    /* Assign local copies of grid functions */
+    CCTK_REAL g11L = g11[index];
+    CCTK_REAL g12L = g12[index];
+    CCTK_REAL g13L = g13[index];
+    CCTK_REAL g22L = g22[index];
+    CCTK_REAL g23L = g23[index];
+    CCTK_REAL g33L = g33[index];
+    CCTK_REAL K11L = K11[index];
+    CCTK_REAL K12L = K12[index];
+    CCTK_REAL K13L = K13[index];
+    CCTK_REAL K22L = K22[index];
+    CCTK_REAL K23L = K23[index];
+    CCTK_REAL K33L = K33[index];
+    
+    /* Include user supplied include files */
+    
+    /* Precompute derivatives */
+    CCTK_REAL const PDstandardNth1g11 = PDstandardNth1(&g11[index]);
+    CCTK_REAL const PDstandardNth2g11 = PDstandardNth2(&g11[index]);
+    CCTK_REAL const PDstandardNth3g11 = PDstandardNth3(&g11[index]);
+    CCTK_REAL const PDstandardNth22g11 = PDstandardNth22(&g11[index]);
+    CCTK_REAL const PDstandardNth33g11 = PDstandardNth33(&g11[index]);
+    CCTK_REAL const PDstandardNth23g11 = PDstandardNth23(&g11[index]);
+    CCTK_REAL const PDstandardNth1g12 = PDstandardNth1(&g12[index]);
+    CCTK_REAL const PDstandardNth2g12 = PDstandardNth2(&g12[index]);
+    CCTK_REAL const PDstandardNth3g12 = PDstandardNth3(&g12[index]);
+    CCTK_REAL const PDstandardNth33g12 = PDstandardNth33(&g12[index]);
+    CCTK_REAL const PDstandardNth12g12 = PDstandardNth12(&g12[index]);
+    CCTK_REAL const PDstandardNth13g12 = PDstandardNth13(&g12[index]);
+    CCTK_REAL const PDstandardNth23g12 = PDstandardNth23(&g12[index]);
+    CCTK_REAL const PDstandardNth1g13 = PDstandardNth1(&g13[index]);
+    CCTK_REAL const PDstandardNth2g13 = PDstandardNth2(&g13[index]);
+    CCTK_REAL const PDstandardNth3g13 = PDstandardNth3(&g13[index]);
+    CCTK_REAL const PDstandardNth22g13 = PDstandardNth22(&g13[index]);
+    CCTK_REAL const PDstandardNth12g13 = PDstandardNth12(&g13[index]);
+    CCTK_REAL const PDstandardNth13g13 = PDstandardNth13(&g13[index]);
+    CCTK_REAL const PDstandardNth23g13 = PDstandardNth23(&g13[index]);
+    CCTK_REAL const PDstandardNth1g22 = PDstandardNth1(&g22[index]);
+    CCTK_REAL const PDstandardNth2g22 = PDstandardNth2(&g22[index]);
+    CCTK_REAL const PDstandardNth3g22 = PDstandardNth3(&g22[index]);
+    CCTK_REAL const PDstandardNth11g22 = PDstandardNth11(&g22[index]);
+    CCTK_REAL const PDstandardNth33g22 = PDstandardNth33(&g22[index]);
+    CCTK_REAL const PDstandardNth13g22 = PDstandardNth13(&g22[index]);
+    CCTK_REAL const PDstandardNth1g23 = PDstandardNth1(&g23[index]);
+    CCTK_REAL const PDstandardNth2g23 = PDstandardNth2(&g23[index]);
+    CCTK_REAL const PDstandardNth3g23 = PDstandardNth3(&g23[index]);
+    CCTK_REAL const PDstandardNth11g23 = PDstandardNth11(&g23[index]);
+    CCTK_REAL const PDstandardNth12g23 = PDstandardNth12(&g23[index]);
+    CCTK_REAL const PDstandardNth13g23 = PDstandardNth13(&g23[index]);
+    CCTK_REAL const PDstandardNth23g23 = PDstandardNth23(&g23[index]);
+    CCTK_REAL const PDstandardNth1g33 = PDstandardNth1(&g33[index]);
+    CCTK_REAL const PDstandardNth2g33 = PDstandardNth2(&g33[index]);
+    CCTK_REAL const PDstandardNth3g33 = PDstandardNth3(&g33[index]);
+    CCTK_REAL const PDstandardNth11g33 = PDstandardNth11(&g33[index]);
+    CCTK_REAL const PDstandardNth22g33 = PDstandardNth22(&g33[index]);
+    CCTK_REAL const PDstandardNth12g33 = PDstandardNth12(&g33[index]);
+    CCTK_REAL const PDstandardNth2K11 = PDstandardNth2(&K11[index]);
+    CCTK_REAL const PDstandardNth3K11 = PDstandardNth3(&K11[index]);
+    CCTK_REAL const PDstandardNth1K12 = PDstandardNth1(&K12[index]);
+    CCTK_REAL const PDstandardNth2K12 = PDstandardNth2(&K12[index]);
+    CCTK_REAL const PDstandardNth3K12 = PDstandardNth3(&K12[index]);
+    CCTK_REAL const PDstandardNth1K13 = PDstandardNth1(&K13[index]);
+    CCTK_REAL const PDstandardNth2K13 = PDstandardNth2(&K13[index]);
+    CCTK_REAL const PDstandardNth3K13 = PDstandardNth3(&K13[index]);
+    CCTK_REAL const PDstandardNth1K22 = PDstandardNth1(&K22[index]);
+    CCTK_REAL const PDstandardNth3K22 = PDstandardNth3(&K22[index]);
+    CCTK_REAL const PDstandardNth1K23 = PDstandardNth1(&K23[index]);
+    CCTK_REAL const PDstandardNth2K23 = PDstandardNth2(&K23[index]);
+    CCTK_REAL const PDstandardNth3K23 = PDstandardNth3(&K23[index]);
+    CCTK_REAL const PDstandardNth1K33 = PDstandardNth1(&K33[index]);
+    CCTK_REAL const PDstandardNth2K33 = PDstandardNth2(&K33[index]);
+    
+    /* Calculate temporaries and grid functions */
+    CCTK_REAL detg = 2*g12L*g13L*g23L + g33L*(g11L*g22L - SQR(g12L)) - 
+      g22L*SQR(g13L) - g11L*SQR(g23L);
+    
+    CCTK_REAL gu11 = INV(detg)*(g22L*g33L - SQR(g23L));
+    
+    CCTK_REAL gu12 = (g13L*g23L - g12L*g33L)*INV(detg);
+    
+    CCTK_REAL gu13 = (-(g13L*g22L) + g12L*g23L)*INV(detg);
+    
+    CCTK_REAL gu21 = (g13L*g23L - g12L*g33L)*INV(detg);
+    
+    CCTK_REAL gu22 = INV(detg)*(g11L*g33L - SQR(g13L));
+    
+    CCTK_REAL gu23 = (g12L*g13L - g11L*g23L)*INV(detg);
+    
+    CCTK_REAL gu31 = (-(g13L*g22L) + g12L*g23L)*INV(detg);
+    
+    CCTK_REAL gu32 = (g12L*g13L - g11L*g23L)*INV(detg);
+    
+    CCTK_REAL gu33 = INV(detg)*(g11L*g22L - SQR(g12L));
+    
+    CCTK_REAL G111 = 0.5*(gu11*PDstandardNth1g11 + 
+      2*(gu12*PDstandardNth1g12 + gu13*PDstandardNth1g13) - 
+      gu12*PDstandardNth2g11 - gu13*PDstandardNth3g11);
+    
+    CCTK_REAL G211 = 0.5*(gu21*PDstandardNth1g11 + 
+      2*(gu22*PDstandardNth1g12 + gu23*PDstandardNth1g13) - 
+      gu22*PDstandardNth2g11 - gu23*PDstandardNth3g11);
+    
+    CCTK_REAL G311 = 0.5*(gu31*PDstandardNth1g11 + 
+      2*(gu32*PDstandardNth1g12 + gu33*PDstandardNth1g13) - 
+      gu32*PDstandardNth2g11 - gu33*PDstandardNth3g11);
+    
+    CCTK_REAL G112 = 0.5*(gu12*PDstandardNth1g22 + gu11*PDstandardNth2g11 
+      + gu13*(PDstandardNth1g23 + PDstandardNth2g13 - PDstandardNth3g12));
+    
+    CCTK_REAL G212 = 0.5*(gu22*PDstandardNth1g22 + gu21*PDstandardNth2g11 
+      + gu23*(PDstandardNth1g23 + PDstandardNth2g13 - PDstandardNth3g12));
+    
+    CCTK_REAL G312 = 0.5*(gu32*PDstandardNth1g22 + gu31*PDstandardNth2g11 
+      + gu33*(PDstandardNth1g23 + PDstandardNth2g13 - PDstandardNth3g12));
+    
+    CCTK_REAL G113 = 0.5*(gu13*PDstandardNth1g33 + gu11*PDstandardNth3g11 
+      + gu12*(PDstandardNth1g23 - PDstandardNth2g13 + PDstandardNth3g12));
+    
+    CCTK_REAL G213 = 0.5*(gu23*PDstandardNth1g33 + gu21*PDstandardNth3g11 
+      + gu22*(PDstandardNth1g23 - PDstandardNth2g13 + PDstandardNth3g12));
+    
+    CCTK_REAL G313 = 0.5*(gu33*PDstandardNth1g33 + gu31*PDstandardNth3g11 
+      + gu32*(PDstandardNth1g23 - PDstandardNth2g13 + PDstandardNth3g12));
+    
+    CCTK_REAL G122 = 0.5*(gu11*(-PDstandardNth1g22 + 2*PDstandardNth2g12) 
+      + gu12*PDstandardNth2g22 + gu13*(2*PDstandardNth2g23 - 
+      PDstandardNth3g22));
+    
+    CCTK_REAL G222 = 0.5*(gu21*(-PDstandardNth1g22 + 2*PDstandardNth2g12) 
+      + gu22*PDstandardNth2g22 + gu23*(2*PDstandardNth2g23 - 
+      PDstandardNth3g22));
+    
+    CCTK_REAL G322 = 0.5*(gu31*(-PDstandardNth1g22 + 2*PDstandardNth2g12) 
+      + gu32*PDstandardNth2g22 + gu33*(2*PDstandardNth2g23 - 
+      PDstandardNth3g22));
+    
+    CCTK_REAL G123 = 0.5*(gu13*PDstandardNth2g33 + 
+      gu11*(-PDstandardNth1g23 + PDstandardNth2g13 + PDstandardNth3g12) + 
+      gu12*PDstandardNth3g22);
+    
+    CCTK_REAL G223 = 0.5*(gu23*PDstandardNth2g33 + 
+      gu21*(-PDstandardNth1g23 + PDstandardNth2g13 + PDstandardNth3g12) + 
+      gu22*PDstandardNth3g22);
+    
+    CCTK_REAL G323 = 0.5*(gu33*PDstandardNth2g33 + 
+      gu31*(-PDstandardNth1g23 + PDstandardNth2g13 + PDstandardNth3g12) + 
+      gu32*PDstandardNth3g22);
+    
+    CCTK_REAL G133 = 0.5*(gu11*(-PDstandardNth1g33 + 2*PDstandardNth3g13) 
+      + gu12*(-PDstandardNth2g33 + 2*PDstandardNth3g23) + 
+      gu13*PDstandardNth3g33);
+    
+    CCTK_REAL G233 = 0.5*(gu21*(-PDstandardNth1g33 + 2*PDstandardNth3g13) 
+      + gu22*(-PDstandardNth2g33 + 2*PDstandardNth3g23) + 
+      gu23*PDstandardNth3g33);
+    
+    CCTK_REAL G333 = 0.5*(gu31*(-PDstandardNth1g33 + 2*PDstandardNth3g13) 
+      + gu32*(-PDstandardNth2g33 + 2*PDstandardNth3g23) + 
+      gu33*PDstandardNth3g33);
+    
+    CCTK_REAL R11 = -(G111*(G111 + G122 + G133)) - G211*(G211 + G222 + 
+      G233) - G311*(G311 + G322 + G333) + 0.5*(-(gu22*(PDstandardNth11g22 - 
+      2*PDstandardNth12g12 + PDstandardNth22g11)) + gu23*(-PDstandardNth11g23 
+      + PDstandardNth12g13 + PDstandardNth13g12 - PDstandardNth23g11) + 
+      gu32*(-PDstandardNth11g23 + PDstandardNth12g13 + PDstandardNth13g12 - 
+      PDstandardNth23g11) - gu33*(PDstandardNth11g33 - 2*PDstandardNth13g13 + 
+      PDstandardNth33g11)) + SQR(G111) + SQR(G112) + SQR(G113) + SQR(G211) + 
+      SQR(G212) + SQR(G213) + SQR(G311) + SQR(G312) + SQR(G313);
+    
+    CCTK_REAL R12 = 0.5*(2*(G113*G123 + G213*G223 + G313*G323) - 
+      2*(G112*G133 + G212*G233 + G312*G333 + gu12*PDstandardNth12g12) + 
+      gu13*(PDstandardNth11g23 - PDstandardNth12g13 - PDstandardNth13g12) + 
+      gu12*(PDstandardNth11g22 + PDstandardNth22g11) + 
+      gu32*PDstandardNth22g13 + gu13*PDstandardNth23g11 + 
+      gu32*(-PDstandardNth12g23 + PDstandardNth13g22 - PDstandardNth23g12) + 
+      gu33*PDstandardNth23g13 + gu33*(-PDstandardNth12g33 + 
+      PDstandardNth13g23 - PDstandardNth33g12));
+    
+    CCTK_REAL R13 = 0.5*(2*(G112*G123 + G212*G223 + G312*G323) - 
+      2*(G113*G122 + G213*G222 + G313*G322 + gu13*PDstandardNth13g13) + 
+      gu12*(PDstandardNth11g23 - PDstandardNth12g13 - PDstandardNth13g12 + 
+      PDstandardNth23g11) + gu22*(PDstandardNth12g23 - PDstandardNth13g22 - 
+      PDstandardNth22g13 + PDstandardNth23g12) + gu13*(PDstandardNth11g33 + 
+      PDstandardNth33g11) + gu23*(PDstandardNth12g33 - PDstandardNth13g23 - 
+      PDstandardNth23g13 + PDstandardNth33g12));
+    
+    CCTK_REAL R22 = -(G122*(G111 + G122 + G133)) - G222*(G211 + G222 + 
+      G233) - G322*(G311 + G322 + G333) + 0.5*(-(gu11*(PDstandardNth11g22 - 
+      2*PDstandardNth12g12 + PDstandardNth22g11)) + gu13*(PDstandardNth12g23 
+      - PDstandardNth13g22 - PDstandardNth22g13 + PDstandardNth23g12) + 
+      gu31*(PDstandardNth12g23 - PDstandardNth13g22 - PDstandardNth22g13 + 
+      PDstandardNth23g12) - gu33*(PDstandardNth22g33 - 2*PDstandardNth23g23 + 
+      PDstandardNth33g22)) + SQR(G112) + SQR(G122) + SQR(G123) + SQR(G212) + 
+      SQR(G222) + SQR(G223) + SQR(G312) + SQR(G322) + SQR(G323);
+    
+    CCTK_REAL R23 = 0.5*(2*(G112*G113 + G212*G213 + G312*G313) + 
+      gu11*(-PDstandardNth11g23 + PDstandardNth12g13 + PDstandardNth13g12 - 
+      PDstandardNth23g11) + gu21*(-PDstandardNth12g23 + PDstandardNth13g22 + 
+      PDstandardNth22g13 - PDstandardNth23g12) - 2*(G111*G123 + G211*G223 + 
+      G311*G323 + gu23*PDstandardNth23g23) + gu13*(PDstandardNth12g33 - 
+      PDstandardNth13g23 - PDstandardNth23g13 + PDstandardNth33g12) + 
+      gu23*(PDstandardNth22g33 + PDstandardNth33g22));
+    
+    CCTK_REAL R33 = -(G133*(G111 + G122 + G133)) - G233*(G211 + G222 + 
+      G233) - G333*(G311 + G322 + G333) + 0.5*(-(gu11*(PDstandardNth11g33 - 
+      2*PDstandardNth13g13 + PDstandardNth33g11)) + gu12*(-PDstandardNth12g33 
+      + PDstandardNth13g23 + PDstandardNth23g13 - PDstandardNth33g12) + 
+      gu21*(-PDstandardNth12g33 + PDstandardNth13g23 + PDstandardNth23g13 - 
+      PDstandardNth33g12) - gu22*(PDstandardNth22g33 - 2*PDstandardNth23g23 + 
+      PDstandardNth33g22)) + SQR(G113) + SQR(G123) + SQR(G133) + SQR(G213) + 
+      SQR(G223) + SQR(G233) + SQR(G313) + SQR(G323) + SQR(G333);
+    
+    CCTK_REAL trR = gu11*R11 + (gu12 + gu21)*R12 + (gu13 + gu31)*R13 + 
+      gu22*R22 + (gu23 + gu32)*R23 + gu33*R33;
+    
+    CCTK_REAL Km11 = gu11*K11L + gu12*K12L + gu13*K13L;
+    
+    CCTK_REAL Km21 = gu21*K11L + gu22*K12L + gu23*K13L;
+    
+    CCTK_REAL Km31 = gu31*K11L + gu32*K12L + gu33*K13L;
+    
+    CCTK_REAL Km12 = gu11*K12L + gu12*K22L + gu13*K23L;
+    
+    CCTK_REAL Km22 = gu21*K12L + gu22*K22L + gu23*K23L;
+    
+    CCTK_REAL Km32 = gu31*K12L + gu32*K22L + gu33*K23L;
+    
+    CCTK_REAL Km13 = gu11*K13L + gu12*K23L + gu13*K33L;
+    
+    CCTK_REAL Km23 = gu21*K13L + gu22*K23L + gu23*K33L;
+    
+    CCTK_REAL Km33 = gu31*K13L + gu32*K23L + gu33*K33L;
+    
+    CCTK_REAL trK = Km11 + Km22 + Km33;
+    
+    CCTK_REAL HL = -2*(Km12*Km21 + Km13*Km31 + Km23*Km32) + trR - 
+      SQR(Km11) - SQR(Km22) - SQR(Km33) + SQR(trK);
+    
+    CCTK_REAL M1L = gu21*(-(G112*K11L) + (G111 - G212)*K12L - G312*K13L + 
+      G211*K22L + G311*K23L - PDstandardNth1K12 + PDstandardNth2K11) + 
+      gu22*(-(G122*K11L) + (G112 - G222)*K12L - G322*K13L + G212*K22L + 
+      G312*K23L - PDstandardNth1K22 + PDstandardNth2K12) + gu23*(-(G123*K11L) 
+      + (G113 - G223)*K12L - G323*K13L + G213*K22L + G313*K23L - 
+      PDstandardNth1K23 + PDstandardNth2K13) + gu31*(-(G113*K11L) - G213*K12L 
+      + (G111 - G313)*K13L + G211*K23L + G311*K33L - PDstandardNth1K13 + 
+      PDstandardNth3K11) + gu32*(-(G123*K11L) - G223*K12L + (G112 - 
+      G323)*K13L + G212*K23L + G312*K33L - PDstandardNth1K23 + 
+      PDstandardNth3K12) + gu33*(-(G133*K11L) - G233*K12L + (G113 - 
+      G333)*K13L + G213*K23L + G313*K33L - PDstandardNth1K33 + 
+      PDstandardNth3K13);
+    
+    CCTK_REAL M2L = gu11*(G112*K11L + (-G111 + G212)*K12L + G312*K13L - 
+      G211*K22L - G311*K23L + PDstandardNth1K12 - PDstandardNth2K11) + 
+      gu12*(G122*K11L + (-G112 + G222)*K12L + G322*K13L - G212*K22L - 
+      G312*K23L + PDstandardNth1K22 - PDstandardNth2K12) + gu13*(G123*K11L + 
+      (-G113 + G223)*K12L + G323*K13L - G213*K22L - G313*K23L + 
+      PDstandardNth1K23 - PDstandardNth2K13) + gu31*(-(G113*K12L) + G112*K13L 
+      - G213*K22L + (G212 - G313)*K23L + G312*K33L - PDstandardNth2K13 + 
+      PDstandardNth3K12) + gu32*(-(G123*K12L) + G122*K13L - G223*K22L + (G222 
+      - G323)*K23L + G322*K33L - PDstandardNth2K23 + PDstandardNth3K22) + 
+      gu33*(-(G133*K12L) + G123*K13L - G233*K22L + (G223 - G333)*K23L + 
+      G323*K33L - PDstandardNth2K33 + PDstandardNth3K23);
+    
+    CCTK_REAL M3L = gu11*(G113*K11L + G213*K12L + (-G111 + G313)*K13L - 
+      G211*K23L - G311*K33L + PDstandardNth1K13 - PDstandardNth3K11) + 
+      gu12*(G123*K11L + G223*K12L + (-G112 + G323)*K13L - G212*K23L - 
+      G312*K33L + PDstandardNth1K23 - PDstandardNth3K12) + gu21*(G113*K12L - 
+      G112*K13L + G213*K22L + (-G212 + G313)*K23L - G312*K33L + 
+      PDstandardNth2K13 - PDstandardNth3K12) + gu13*(G133*K11L + G233*K12L + 
+      (-G113 + G333)*K13L - G213*K23L - G313*K33L + PDstandardNth1K33 - 
+      PDstandardNth3K13) + gu22*(G123*K12L - G122*K13L + G223*K22L + (-G222 + 
+      G323)*K23L - G322*K33L + PDstandardNth2K23 - PDstandardNth3K22) + 
+      gu23*(G133*K12L - G123*K13L + G233*K22L + (-G223 + G333)*K23L - 
+      G323*K33L + PDstandardNth2K33 - PDstandardNth3K23);
+    
+    
+    /* Copy local copies back to grid functions */
+    H[index] = HL;
+    M1[index] = M1L;
+    M2[index] = M2L;
+    M3[index] = M3L;
+  }
+  LC_ENDLOOP3 (ML_ADM_constraints);
+}
+
+extern "C" void ML_ADM_constraints(CCTK_ARGUMENTS)
+{
+  DECLARE_CCTK_ARGUMENTS;
+  DECLARE_CCTK_PARAMETERS;
+  
+  GenericFD_LoopOverInterior(cctkGH, &ML_ADM_constraints_Body);
+}
diff --git a/ML_ADM/src/ML_ADM_constraints_boundary.cc b/ML_ADM/src/ML_ADM_constraints_boundary.cc
new file mode 100644
index 0000000..d607362
--- /dev/null
+++ b/ML_ADM/src/ML_ADM_constraints_boundary.cc
@@ -0,0 +1,130 @@
+/*  File produced by Kranc */
+
+#define KRANC_C
+
+#include <assert.h>
+#include <math.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include "cctk.h"
+#include "cctk_Arguments.h"
+#include "cctk_Parameters.h"
+#include "GenericFD.h"
+#include "Differencing.h"
+#include "loopcontrol.h"
+
+/* Define macros used in calculations */
+#define INITVALUE (42)
+#define QAD(x) (SQR(SQR(x)))
+#define INV(x) ((1.0) / (x))
+#define SQR(x) ((x) * (x))
+#define CUB(x) ((x) * (x) * (x))
+
+extern "C" void ML_ADM_constraints_boundary_SelectBCs(CCTK_ARGUMENTS)
+{
+  DECLARE_CCTK_ARGUMENTS;
+  DECLARE_CCTK_PARAMETERS;
+  
+  CCTK_INT ierr = 0;
+  ierr = Boundary_SelectGroupForBC(cctkGH, CCTK_ALL_FACES, GenericFD_GetBoundaryWidth(cctkGH), -1 /* no table */, "ML_ADM::ML_Ham","flat");
+  if (ierr < 0)
+    CCTK_WARN(1, "Failed to register flat BC for ML_ADM::ML_Ham.");
+  ierr = Boundary_SelectGroupForBC(cctkGH, CCTK_ALL_FACES, GenericFD_GetBoundaryWidth(cctkGH), -1 /* no table */, "ML_ADM::ML_mom","flat");
+  if (ierr < 0)
+    CCTK_WARN(1, "Failed to register flat BC for ML_ADM::ML_mom.");
+  return;
+}
+
+static void ML_ADM_constraints_boundary_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[])
+{
+  DECLARE_CCTK_ARGUMENTS;
+  DECLARE_CCTK_PARAMETERS;
+  
+  
+  /* Declare finite differencing variables */
+  
+  if (verbose > 1)
+  {
+    CCTK_VInfo(CCTK_THORNSTRING,"Entering ML_ADM_constraints_boundary_Body");
+  }
+  
+  if (cctk_iteration % ML_ADM_constraints_boundary_calc_every != ML_ADM_constraints_boundary_calc_offset)
+  {
+    return;
+  }
+  
+  const char *groups[] = {"ML_ADM::ML_Ham","ML_ADM::ML_mom"};
+  GenericFD_AssertGroupStorage(cctkGH, "ML_ADM_constraints_boundary", 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);
+  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 dxi = INV(dx);
+  CCTK_REAL const dyi = INV(dy);
+  CCTK_REAL const dzi = INV(dz);
+  CCTK_REAL const khalf = 0.5;
+  CCTK_REAL const kthird = 1/3.0;
+  CCTK_REAL const ktwothird = 2.0/3.0;
+  CCTK_REAL const kfourthird = 4.0/3.0;
+  CCTK_REAL const keightthird = 8.0/3.0;
+  CCTK_REAL const hdxi = 0.5 * dxi;
+  CCTK_REAL const hdyi = 0.5 * dyi;
+  CCTK_REAL const hdzi = 0.5 * dzi;
+  
+  /* Initialize predefined quantities */
+  CCTK_REAL const p1o12dx = 0.0833333333333333333333333333333*INV(dx);
+  CCTK_REAL const p1o12dy = 0.0833333333333333333333333333333*INV(dy);
+  CCTK_REAL const p1o12dz = 0.0833333333333333333333333333333*INV(dz);
+  CCTK_REAL const p1o144dxdy = 0.00694444444444444444444444444444*INV(dx)*INV(dy);
+  CCTK_REAL const p1o144dxdz = 0.00694444444444444444444444444444*INV(dx)*INV(dz);
+  CCTK_REAL const p1o144dydz = 0.00694444444444444444444444444444*INV(dy)*INV(dz);
+  CCTK_REAL const pm1o12dx2 = -0.0833333333333333333333333333333*INV(SQR(dx));
+  CCTK_REAL const pm1o12dy2 = -0.0833333333333333333333333333333*INV(SQR(dy));
+  CCTK_REAL const pm1o12dz2 = -0.0833333333333333333333333333333*INV(SQR(dz));
+  
+  /* Loop over the grid points */
+  #pragma omp parallel
+  LC_LOOP3 (ML_ADM_constraints_boundary,
+    i,j,k, min[0],min[1],min[2], max[0],max[1],max[2],
+    cctk_lsh[0],cctk_lsh[1],cctk_lsh[2])
+  {
+    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 HL = 0;
+    
+    CCTK_REAL M1L = 0;
+    
+    CCTK_REAL M2L = 0;
+    
+    CCTK_REAL M3L = 0;
+    
+    
+    /* Copy local copies back to grid functions */
+    H[index] = HL;
+    M1[index] = M1L;
+    M2[index] = M2L;
+    M3[index] = M3L;
+  }
+  LC_ENDLOOP3 (ML_ADM_constraints_boundary);
+}
+
+extern "C" void ML_ADM_constraints_boundary(CCTK_ARGUMENTS)
+{
+  DECLARE_CCTK_ARGUMENTS;
+  DECLARE_CCTK_PARAMETERS;
+  
+  GenericFD_LoopOverBoundaryWithGhosts(cctkGH, &ML_ADM_constraints_boundary_Body);
+}
diff --git a/ML_ADM/src/ML_ADM_convertFromADMBase.cc b/ML_ADM/src/ML_ADM_convertFromADMBase.cc
new file mode 100644
index 0000000..85a7aa9
--- /dev/null
+++ b/ML_ADM/src/ML_ADM_convertFromADMBase.cc
@@ -0,0 +1,167 @@
+/*  File produced by Kranc */
+
+#define KRANC_C
+
+#include <assert.h>
+#include <math.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include "cctk.h"
+#include "cctk_Arguments.h"
+#include "cctk_Parameters.h"
+#include "GenericFD.h"
+#include "Differencing.h"
+#include "loopcontrol.h"
+
+/* Define macros used in calculations */
+#define INITVALUE (42)
+#define QAD(x) (SQR(SQR(x)))
+#define INV(x) ((1.0) / (x))
+#define SQR(x) ((x) * (x))
+#define CUB(x) ((x) * (x) * (x))
+
+static void ML_ADM_convertFromADMBase_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[])
+{
+  DECLARE_CCTK_ARGUMENTS;
+  DECLARE_CCTK_PARAMETERS;
+  
+  
+  /* Declare finite differencing variables */
+  
+  if (verbose > 1)
+  {
+    CCTK_VInfo(CCTK_THORNSTRING,"Entering ML_ADM_convertFromADMBase_Body");
+  }
+  
+  if (cctk_iteration % ML_ADM_convertFromADMBase_calc_every != ML_ADM_convertFromADMBase_calc_offset)
+  {
+    return;
+  }
+  
+  const char *groups[] = {"ADMBase::curv","ADMBase::lapse","ADMBase::metric","ADMBase::shift","ML_ADM::ML_curv","ML_ADM::ML_lapse","ML_ADM::ML_metric","ML_ADM::ML_shift"};
+  GenericFD_AssertGroupStorage(cctkGH, "ML_ADM_convertFromADMBase", 8, 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);
+  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 dxi = INV(dx);
+  CCTK_REAL const dyi = INV(dy);
+  CCTK_REAL const dzi = INV(dz);
+  CCTK_REAL const khalf = 0.5;
+  CCTK_REAL const kthird = 1/3.0;
+  CCTK_REAL const ktwothird = 2.0/3.0;
+  CCTK_REAL const kfourthird = 4.0/3.0;
+  CCTK_REAL const keightthird = 8.0/3.0;
+  CCTK_REAL const hdxi = 0.5 * dxi;
+  CCTK_REAL const hdyi = 0.5 * dyi;
+  CCTK_REAL const hdzi = 0.5 * dzi;
+  
+  /* Initialize predefined quantities */
+  CCTK_REAL const p1o12dx = 0.0833333333333333333333333333333*INV(dx);
+  CCTK_REAL const p1o12dy = 0.0833333333333333333333333333333*INV(dy);
+  CCTK_REAL const p1o12dz = 0.0833333333333333333333333333333*INV(dz);
+  CCTK_REAL const p1o144dxdy = 0.00694444444444444444444444444444*INV(dx)*INV(dy);
+  CCTK_REAL const p1o144dxdz = 0.00694444444444444444444444444444*INV(dx)*INV(dz);
+  CCTK_REAL const p1o144dydz = 0.00694444444444444444444444444444*INV(dy)*INV(dz);
+  CCTK_REAL const pm1o12dx2 = -0.0833333333333333333333333333333*INV(SQR(dx));
+  CCTK_REAL const pm1o12dy2 = -0.0833333333333333333333333333333*INV(SQR(dy));
+  CCTK_REAL const pm1o12dz2 = -0.0833333333333333333333333333333*INV(SQR(dz));
+  
+  /* Loop over the grid points */
+  #pragma omp parallel
+  LC_LOOP3 (ML_ADM_convertFromADMBase,
+    i,j,k, min[0],min[1],min[2], max[0],max[1],max[2],
+    cctk_lsh[0],cctk_lsh[1],cctk_lsh[2])
+  {
+    ptrdiff_t const index = di*i + dj*j + dk*k;
+    
+    /* Assign local copies of grid functions */
+    CCTK_REAL alpL = alp[index];
+    CCTK_REAL betaxL = betax[index];
+    CCTK_REAL betayL = betay[index];
+    CCTK_REAL betazL = betaz[index];
+    CCTK_REAL gxxL = gxx[index];
+    CCTK_REAL gxyL = gxy[index];
+    CCTK_REAL gxzL = gxz[index];
+    CCTK_REAL gyyL = gyy[index];
+    CCTK_REAL gyzL = gyz[index];
+    CCTK_REAL gzzL = gzz[index];
+    CCTK_REAL kxxL = kxx[index];
+    CCTK_REAL kxyL = kxy[index];
+    CCTK_REAL kxzL = kxz[index];
+    CCTK_REAL kyyL = kyy[index];
+    CCTK_REAL kyzL = kyz[index];
+    CCTK_REAL kzzL = kzz[index];
+    
+    /* Include user supplied include files */
+    
+    /* Precompute derivatives */
+    
+    /* Calculate temporaries and grid functions */
+    CCTK_REAL g11L = gxxL;
+    
+    CCTK_REAL g12L = gxyL;
+    
+    CCTK_REAL g13L = gxzL;
+    
+    CCTK_REAL g22L = gyyL;
+    
+    CCTK_REAL g23L = gyzL;
+    
+    CCTK_REAL g33L = gzzL;
+    
+    CCTK_REAL K11L = kxxL;
+    
+    CCTK_REAL K12L = kxyL;
+    
+    CCTK_REAL K13L = kxzL;
+    
+    CCTK_REAL K22L = kyyL;
+    
+    CCTK_REAL K23L = kyzL;
+    
+    CCTK_REAL K33L = kzzL;
+    
+    CCTK_REAL alphaL = alpL;
+    
+    CCTK_REAL beta1L = betaxL;
+    
+    CCTK_REAL beta2L = betayL;
+    
+    CCTK_REAL beta3L = betazL;
+    
+    
+    /* Copy local copies back to grid functions */
+    alpha[index] = alphaL;
+    beta1[index] = beta1L;
+    beta2[index] = beta2L;
+    beta3[index] = beta3L;
+    g11[index] = g11L;
+    g12[index] = g12L;
+    g13[index] = g13L;
+    g22[index] = g22L;
+    g23[index] = g23L;
+    g33[index] = g33L;
+    K11[index] = K11L;
+    K12[index] = K12L;
+    K13[index] = K13L;
+    K22[index] = K22L;
+    K23[index] = K23L;
+    K33[index] = K33L;
+  }
+  LC_ENDLOOP3 (ML_ADM_convertFromADMBase);
+}
+
+extern "C" void ML_ADM_convertFromADMBase(CCTK_ARGUMENTS)
+{
+  DECLARE_CCTK_ARGUMENTS;
+  DECLARE_CCTK_PARAMETERS;
+  
+  GenericFD_LoopOverEverything(cctkGH, &ML_ADM_convertFromADMBase_Body);
+}
diff --git a/ML_ADM/src/ML_ADM_convertToADMBase.cc b/ML_ADM/src/ML_ADM_convertToADMBase.cc
new file mode 100644
index 0000000..6df8ec1
--- /dev/null
+++ b/ML_ADM/src/ML_ADM_convertToADMBase.cc
@@ -0,0 +1,179 @@
+/*  File produced by Kranc */
+
+#define KRANC_C
+
+#include <assert.h>
+#include <math.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include "cctk.h"
+#include "cctk_Arguments.h"
+#include "cctk_Parameters.h"
+#include "GenericFD.h"
+#include "Differencing.h"
+#include "loopcontrol.h"
+
+/* Define macros used in calculations */
+#define INITVALUE (42)
+#define QAD(x) (SQR(SQR(x)))
+#define INV(x) ((1.0) / (x))
+#define SQR(x) ((x) * (x))
+#define CUB(x) ((x) * (x) * (x))
+
+static void ML_ADM_convertToADMBase_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[])
+{
+  DECLARE_CCTK_ARGUMENTS;
+  DECLARE_CCTK_PARAMETERS;
+  
+  
+  /* Declare finite differencing variables */
+  
+  if (verbose > 1)
+  {
+    CCTK_VInfo(CCTK_THORNSTRING,"Entering ML_ADM_convertToADMBase_Body");
+  }
+  
+  if (cctk_iteration % ML_ADM_convertToADMBase_calc_every != ML_ADM_convertToADMBase_calc_offset)
+  {
+    return;
+  }
+  
+  const char *groups[] = {"ADMBase::curv","ADMBase::dtlapse","ADMBase::dtshift","ADMBase::lapse","ADMBase::metric","ADMBase::shift","ML_ADM::ML_curv","ML_ADM::ML_lapse","ML_ADM::ML_metric","ML_ADM::ML_shift"};
+  GenericFD_AssertGroupStorage(cctkGH, "ML_ADM_convertToADMBase", 10, groups);
+  
+  /* Include user-supplied include files */
+  
+  /* Initialise finite differencing variables */
+  ptrdiff_t const di = 1;
+  ptrdiff_t const dj = CCTK_GFINDEX3D(cctkGH,0,1,0) - CCTK_GFINDEX3D(cctkGH,0,0,0);
+  ptrdiff_t const dk = CCTK_GFINDEX3D(cctkGH,0,0,1) - CCTK_GFINDEX3D(cctkGH,0,0,0);
+  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 dxi = INV(dx);
+  CCTK_REAL const dyi = INV(dy);
+  CCTK_REAL const dzi = INV(dz);
+  CCTK_REAL const khalf = 0.5;
+  CCTK_REAL const kthird = 1/3.0;
+  CCTK_REAL const ktwothird = 2.0/3.0;
+  CCTK_REAL const kfourthird = 4.0/3.0;
+  CCTK_REAL const keightthird = 8.0/3.0;
+  CCTK_REAL const hdxi = 0.5 * dxi;
+  CCTK_REAL const hdyi = 0.5 * dyi;
+  CCTK_REAL const hdzi = 0.5 * dzi;
+  
+  /* Initialize predefined quantities */
+  CCTK_REAL const p1o12dx = 0.0833333333333333333333333333333*INV(dx);
+  CCTK_REAL const p1o12dy = 0.0833333333333333333333333333333*INV(dy);
+  CCTK_REAL const p1o12dz = 0.0833333333333333333333333333333*INV(dz);
+  CCTK_REAL const p1o144dxdy = 0.00694444444444444444444444444444*INV(dx)*INV(dy);
+  CCTK_REAL const p1o144dxdz = 0.00694444444444444444444444444444*INV(dx)*INV(dz);
+  CCTK_REAL const p1o144dydz = 0.00694444444444444444444444444444*INV(dy)*INV(dz);
+  CCTK_REAL const pm1o12dx2 = -0.0833333333333333333333333333333*INV(SQR(dx));
+  CCTK_REAL const pm1o12dy2 = -0.0833333333333333333333333333333*INV(SQR(dy));
+  CCTK_REAL const pm1o12dz2 = -0.0833333333333333333333333333333*INV(SQR(dz));
+  
+  /* Loop over the grid points */
+  #pragma omp parallel
+  LC_LOOP3 (ML_ADM_convertToADMBase,
+    i,j,k, min[0],min[1],min[2], max[0],max[1],max[2],
+    cctk_lsh[0],cctk_lsh[1],cctk_lsh[2])
+  {
+    ptrdiff_t const index = di*i + dj*j + dk*k;
+    
+    /* Assign local copies of grid functions */
+    CCTK_REAL alphaL = alpha[index];
+    CCTK_REAL beta1L = beta1[index];
+    CCTK_REAL beta2L = beta2[index];
+    CCTK_REAL beta3L = beta3[index];
+    CCTK_REAL g11L = g11[index];
+    CCTK_REAL g12L = g12[index];
+    CCTK_REAL g13L = g13[index];
+    CCTK_REAL g22L = g22[index];
+    CCTK_REAL g23L = g23[index];
+    CCTK_REAL g33L = g33[index];
+    CCTK_REAL K11L = K11[index];
+    CCTK_REAL K12L = K12[index];
+    CCTK_REAL K13L = K13[index];
+    CCTK_REAL K22L = K22[index];
+    CCTK_REAL K23L = K23[index];
+    CCTK_REAL K33L = K33[index];
+    
+    /* Include user supplied include files */
+    
+    /* Precompute derivatives */
+    
+    /* Calculate temporaries and grid functions */
+    CCTK_REAL gxxL = g11L;
+    
+    CCTK_REAL gxyL = g12L;
+    
+    CCTK_REAL gxzL = g13L;
+    
+    CCTK_REAL gyyL = g22L;
+    
+    CCTK_REAL gyzL = g23L;
+    
+    CCTK_REAL gzzL = g33L;
+    
+    CCTK_REAL kxxL = K11L;
+    
+    CCTK_REAL kxyL = K12L;
+    
+    CCTK_REAL kxzL = K13L;
+    
+    CCTK_REAL kyyL = K22L;
+    
+    CCTK_REAL kyzL = K23L;
+    
+    CCTK_REAL kzzL = K33L;
+    
+    CCTK_REAL alpL = alphaL;
+    
+    CCTK_REAL dtalpL = 0;
+    
+    CCTK_REAL betaxL = beta1L;
+    
+    CCTK_REAL betayL = beta2L;
+    
+    CCTK_REAL betazL = beta3L;
+    
+    CCTK_REAL dtbetaxL = 0;
+    
+    CCTK_REAL dtbetayL = 0;
+    
+    CCTK_REAL dtbetazL = 0;
+    
+    
+    /* Copy local copies back to grid functions */
+    alp[index] = alpL;
+    betax[index] = betaxL;
+    betay[index] = betayL;
+    betaz[index] = betazL;
+    dtalp[index] = dtalpL;
+    dtbetax[index] = dtbetaxL;
+    dtbetay[index] = dtbetayL;
+    dtbetaz[index] = dtbetazL;
+    gxx[index] = gxxL;
+    gxy[index] = gxyL;
+    gxz[index] = gxzL;
+    gyy[index] = gyyL;
+    gyz[index] = gyzL;
+    gzz[index] = gzzL;
+    kxx[index] = kxxL;
+    kxy[index] = kxyL;
+    kxz[index] = kxzL;
+    kyy[index] = kyyL;
+    kyz[index] = kyzL;
+    kzz[index] = kzzL;
+  }
+  LC_ENDLOOP3 (ML_ADM_convertToADMBase);
+}
+
+extern "C" void ML_ADM_convertToADMBase(CCTK_ARGUMENTS)
+{
+  DECLARE_CCTK_ARGUMENTS;
+  DECLARE_CCTK_PARAMETERS;
+  
+  GenericFD_LoopOverEverything(cctkGH, &ML_ADM_convertToADMBase_Body);
+}
diff --git a/ML_ADM/src/RegisterMoL.cc b/ML_ADM/src/RegisterMoL.cc
new file mode 100644
index 0000000..55605d7
--- /dev/null
+++ b/ML_ADM/src/RegisterMoL.cc
@@ -0,0 +1,32 @@
+/*  File produced by Kranc */
+
+#include "cctk.h"
+#include "cctk_Arguments.h"
+#include "cctk_Parameters.h"
+
+extern "C" void ML_ADM_RegisterVars(CCTK_ARGUMENTS)
+{
+  DECLARE_CCTK_ARGUMENTS;
+  DECLARE_CCTK_PARAMETERS;
+  
+  CCTK_INT ierr = 0;
+  
+  /* Register all the evolved grid functions with MoL */
+  ierr += MoLRegisterEvolved(CCTK_VarIndex("ML_ADM::K11"),  CCTK_VarIndex("ML_ADM::K11rhs"));
+  ierr += MoLRegisterEvolved(CCTK_VarIndex("ML_ADM::K12"),  CCTK_VarIndex("ML_ADM::K12rhs"));
+  ierr += MoLRegisterEvolved(CCTK_VarIndex("ML_ADM::K13"),  CCTK_VarIndex("ML_ADM::K13rhs"));
+  ierr += MoLRegisterEvolved(CCTK_VarIndex("ML_ADM::K22"),  CCTK_VarIndex("ML_ADM::K22rhs"));
+  ierr += MoLRegisterEvolved(CCTK_VarIndex("ML_ADM::K23"),  CCTK_VarIndex("ML_ADM::K23rhs"));
+  ierr += MoLRegisterEvolved(CCTK_VarIndex("ML_ADM::K33"),  CCTK_VarIndex("ML_ADM::K33rhs"));
+  ierr += MoLRegisterEvolved(CCTK_VarIndex("ML_ADM::alpha"),  CCTK_VarIndex("ML_ADM::alpharhs"));
+  ierr += MoLRegisterEvolved(CCTK_VarIndex("ML_ADM::g11"),  CCTK_VarIndex("ML_ADM::g11rhs"));
+  ierr += MoLRegisterEvolved(CCTK_VarIndex("ML_ADM::g12"),  CCTK_VarIndex("ML_ADM::g12rhs"));
+  ierr += MoLRegisterEvolved(CCTK_VarIndex("ML_ADM::g13"),  CCTK_VarIndex("ML_ADM::g13rhs"));
+  ierr += MoLRegisterEvolved(CCTK_VarIndex("ML_ADM::g22"),  CCTK_VarIndex("ML_ADM::g22rhs"));
+  ierr += MoLRegisterEvolved(CCTK_VarIndex("ML_ADM::g23"),  CCTK_VarIndex("ML_ADM::g23rhs"));
+  ierr += MoLRegisterEvolved(CCTK_VarIndex("ML_ADM::g33"),  CCTK_VarIndex("ML_ADM::g33rhs"));
+  ierr += MoLRegisterEvolved(CCTK_VarIndex("ML_ADM::beta1"),  CCTK_VarIndex("ML_ADM::beta1rhs"));
+  ierr += MoLRegisterEvolved(CCTK_VarIndex("ML_ADM::beta2"),  CCTK_VarIndex("ML_ADM::beta2rhs"));
+  ierr += MoLRegisterEvolved(CCTK_VarIndex("ML_ADM::beta3"),  CCTK_VarIndex("ML_ADM::beta3rhs"));
+  return;
+}
diff --git a/ML_ADM/src/RegisterSymmetries.cc b/ML_ADM/src/RegisterSymmetries.cc
new file mode 100644
index 0000000..0d05f97
--- /dev/null
+++ b/ML_ADM/src/RegisterSymmetries.cc
@@ -0,0 +1,119 @@
+/*  File produced by Kranc */
+
+#include "cctk.h"
+#include "cctk_Arguments.h"
+#include "cctk_Parameters.h"
+#include "Symmetry.h"
+
+extern "C" void ML_ADM_RegisterSymmetries(CCTK_ARGUMENTS)
+{
+  DECLARE_CCTK_ARGUMENTS;
+  DECLARE_CCTK_PARAMETERS;
+  
+  
+  /* array holding symmetry definitions */
+  CCTK_INT sym[3];
+  
+  
+  /* Register symmetries of grid functions */
+  sym[0] = 1;
+  sym[1] = 1;
+  sym[2] = 1;
+  SetCartSymVN(cctkGH, sym, "ML_ADM::K11");
+  
+  sym[0] = -1;
+  sym[1] = -1;
+  sym[2] = 1;
+  SetCartSymVN(cctkGH, sym, "ML_ADM::K12");
+  
+  sym[0] = -1;
+  sym[1] = 1;
+  sym[2] = -1;
+  SetCartSymVN(cctkGH, sym, "ML_ADM::K13");
+  
+  sym[0] = 1;
+  sym[1] = 1;
+  sym[2] = 1;
+  SetCartSymVN(cctkGH, sym, "ML_ADM::K22");
+  
+  sym[0] = 1;
+  sym[1] = -1;
+  sym[2] = -1;
+  SetCartSymVN(cctkGH, sym, "ML_ADM::K23");
+  
+  sym[0] = 1;
+  sym[1] = 1;
+  sym[2] = 1;
+  SetCartSymVN(cctkGH, sym, "ML_ADM::K33");
+  
+  sym[0] = 1;
+  sym[1] = 1;
+  sym[2] = 1;
+  SetCartSymVN(cctkGH, sym, "ML_ADM::alpha");
+  
+  sym[0] = 1;
+  sym[1] = 1;
+  sym[2] = 1;
+  SetCartSymVN(cctkGH, sym, "ML_ADM::g11");
+  
+  sym[0] = -1;
+  sym[1] = -1;
+  sym[2] = 1;
+  SetCartSymVN(cctkGH, sym, "ML_ADM::g12");
+  
+  sym[0] = -1;
+  sym[1] = 1;
+  sym[2] = -1;
+  SetCartSymVN(cctkGH, sym, "ML_ADM::g13");
+  
+  sym[0] = 1;
+  sym[1] = 1;
+  sym[2] = 1;
+  SetCartSymVN(cctkGH, sym, "ML_ADM::g22");
+  
+  sym[0] = 1;
+  sym[1] = -1;
+  sym[2] = -1;
+  SetCartSymVN(cctkGH, sym, "ML_ADM::g23");
+  
+  sym[0] = 1;
+  sym[1] = 1;
+  sym[2] = 1;
+  SetCartSymVN(cctkGH, sym, "ML_ADM::g33");
+  
+  sym[0] = -1;
+  sym[1] = 1;
+  sym[2] = 1;
+  SetCartSymVN(cctkGH, sym, "ML_ADM::beta1");
+  
+  sym[0] = 1;
+  sym[1] = -1;
+  sym[2] = 1;
+  SetCartSymVN(cctkGH, sym, "ML_ADM::beta2");
+  
+  sym[0] = 1;
+  sym[1] = 1;
+  sym[2] = -1;
+  SetCartSymVN(cctkGH, sym, "ML_ADM::beta3");
+  
+  sym[0] = 1;
+  sym[1] = 1;
+  sym[2] = 1;
+  SetCartSymVN(cctkGH, sym, "ML_ADM::H");
+  
+  sym[0] = -1;
+  sym[1] = 1;
+  sym[2] = 1;
+  SetCartSymVN(cctkGH, sym, "ML_ADM::M1");
+  
+  sym[0] = 1;
+  sym[1] = -1;
+  sym[2] = 1;
+  SetCartSymVN(cctkGH, sym, "ML_ADM::M2");
+  
+  sym[0] = 1;
+  sym[1] = 1;
+  sym[2] = -1;
+  SetCartSymVN(cctkGH, sym, "ML_ADM::M3");
+  
+}
diff --git a/ML_ADM/src/Startup.cc b/ML_ADM/src/Startup.cc
new file mode 100644
index 0000000..97927b7
--- /dev/null
+++ b/ML_ADM/src/Startup.cc
@@ -0,0 +1,10 @@
+/*  File produced by Kranc */
+
+#include "cctk.h"
+
+extern "C" int ML_ADM_Startup(void)
+{
+  const char * banner = "ML_ADM";
+  CCTK_RegisterBanner(banner);
+  return 0;
+}