Merge commit 'origin/master'

Conflicts:

	m/Makefile

11 files changed, 687 insertions, 458 deletions

diff --git a/ML_ADM/interface.ccl b/ML_ADM/interface.ccl
index 305c740..cdbd29e 100644
--- a/ML_ADM/interface.ccl
+++ b/ML_ADM/interface.ccl
@@ -5,7 +5,7 @@
 
 implements: ML_ADM
 
-inherits:   ADMBase TmunuBase Grid GenericFD Boundary
+inherits:   ADMBase Grid GenericFD Boundary
 
 
 
@@ -31,21 +31,21 @@ CCTK_INT FUNCTION Boundary_SelectVarForBC(CCTK_POINTER_TO_CONST IN GH, CCTK_INT
 USES FUNCTION Boundary_SelectVarForBC
 
 public:
-CCTK_REAL Ham type=GF timelevels=1 tags='tensortypealias="Scalar" tensorweight=1.0000000000000000000'
+CCTK_REAL ML_Ham type=GF timelevels=1 tags='tensortypealias="Scalar" tensorweight=1.0000000000000000000'
 {
   H
-} "Ham"
+} "ML_Ham"
 
 public:
-CCTK_REAL mom type=GF timelevels=1 tags='tensortypealias="D" tensorweight=1.0000000000000000000'
+CCTK_REAL ML_mom type=GF timelevels=1 tags='tensortypealias="D" tensorweight=1.0000000000000000000'
 {
   M1,
   M2,
   M3
-} "mom"
+} "ML_mom"
 
 public:
-CCTK_REAL ml_curvrhs type=GF timelevels=1 tags='tensortypealias="DD_sym" tensorweight=1.0000000000000000000'
+CCTK_REAL ML_curvrhs type=GF timelevels=1 tags='tensortypealias="DD_sym" tensorweight=1.0000000000000000000'
 {
   K11rhs,
   K12rhs,
@@ -53,16 +53,16 @@ CCTK_REAL ml_curvrhs type=GF timelevels=1 tags='tensortypealias="DD_sym" tensorw
   K22rhs,
   K23rhs,
   K33rhs
-} "ml_curvrhs"
+} "ML_curvrhs"
 
 public:
-CCTK_REAL ml_lapserhs type=GF timelevels=1 tags='tensortypealias="Scalar" tensorweight=1.0000000000000000000'
+CCTK_REAL ML_lapserhs type=GF timelevels=1 tags='tensortypealias="Scalar" tensorweight=1.0000000000000000000'
 {
   alpharhs
-} "ml_lapserhs"
+} "ML_lapserhs"
 
 public:
-CCTK_REAL ml_metricrhs type=GF timelevels=1 tags='tensortypealias="DD_sym" tensorweight=1.0000000000000000000'
+CCTK_REAL ML_metricrhs type=GF timelevels=1 tags='tensortypealias="DD_sym" tensorweight=1.0000000000000000000'
 {
   g11rhs,
   g12rhs,
@@ -70,18 +70,18 @@ CCTK_REAL ml_metricrhs type=GF timelevels=1 tags='tensortypealias="DD_sym" tenso
   g22rhs,
   g23rhs,
   g33rhs
-} "ml_metricrhs"
+} "ML_metricrhs"
 
 public:
-CCTK_REAL ml_shiftrhs type=GF timelevels=1 tags='tensortypealias="U" tensorweight=1.0000000000000000000'
+CCTK_REAL ML_shiftrhs type=GF timelevels=1 tags='tensortypealias="U" tensorweight=1.0000000000000000000'
 {
   beta1rhs,
   beta2rhs,
   beta3rhs
-} "ml_shiftrhs"
+} "ML_shiftrhs"
 
 public:
-CCTK_REAL ml_curv type=GF timelevels=3 tags='tensortypealias="DD_sym" tensorweight=1.0000000000000000000'
+CCTK_REAL ML_curv type=GF timelevels=3 tags='tensortypealias="DD_sym" tensorweight=1.0000000000000000000'
 {
   K11,
   K12,
@@ -89,16 +89,16 @@ CCTK_REAL ml_curv type=GF timelevels=3 tags='tensortypealias="DD_sym" tensorweig
   K22,
   K23,
   K33
-} "ml_curv"
+} "ML_curv"
 
 public:
-CCTK_REAL ml_lapse type=GF timelevels=3 tags='tensortypealias="Scalar" tensorweight=1.0000000000000000000'
+CCTK_REAL ML_lapse type=GF timelevels=3 tags='tensortypealias="Scalar" tensorweight=1.0000000000000000000'
 {
   alpha
-} "ml_lapse"
+} "ML_lapse"
 
 public:
-CCTK_REAL ml_metric type=GF timelevels=3 tags='tensortypealias="DD_sym" tensorweight=1.0000000000000000000'
+CCTK_REAL ML_metric type=GF timelevels=3 tags='tensortypealias="DD_sym" tensorweight=1.0000000000000000000'
 {
   g11,
   g12,
@@ -106,12 +106,12 @@ CCTK_REAL ml_metric type=GF timelevels=3 tags='tensortypealias="DD_sym" tensorwe
   g22,
   g23,
   g33
-} "ml_metric"
+} "ML_metric"
 
 public:
-CCTK_REAL ml_shift type=GF timelevels=3 tags='tensortypealias="U" tensorweight=1.0000000000000000000'
+CCTK_REAL ML_shift type=GF timelevels=3 tags='tensortypealias="U" tensorweight=1.0000000000000000000'
 {
   beta1,
   beta2,
   beta3
-} "ml_shift"
+} "ML_shift"
diff --git a/ML_ADM/param.ccl b/ML_ADM/param.ccl
index db99057..b22e4ba 100644
--- a/ML_ADM/param.ccl
+++ b/ML_ADM/param.ccl
@@ -19,9 +19,9 @@ USES CCTK_INT MoL_Num_Evolved_Vars
 USES CCTK_INT MoL_Num_Constrained_Vars
 
 restricted:
-CCTK_REAL verbose "verbose"
+CCTK_INT verbose "verbose"
 {
-  "*:*" :: ""
+  *:* :: ""
 } 0
 
 private:
@@ -47,8 +47,8 @@ CCTK_INT ML_ADM_MaxNumEvolvedVars "Number of evolved variables used by this thor
 restricted:
 CCTK_INT ML_ADM_MaxNumConstrainedVars "Number of constrained variables used by this thorn" ACCUMULATOR-BASE=MethodofLines::MoL_Num_Constrained_Vars
 {
-  38:38 :: "Number of constrained variables used by this thorn"
-} 38
+  65:65 :: "Number of constrained variables used by this thorn"
+} 65
 
 restricted:
 CCTK_INT timelevels "Number of active timelevels"
@@ -333,7 +333,7 @@ KEYWORD beta3_bound "Boundary condition to implement"
 } "skip"
 
 private:
-KEYWORD ml_curv_bound "Boundary condition to implement"
+KEYWORD ML_curv_bound "Boundary condition to implement"
 {
   "flat" :: "Flat boundary condition"
   "none" :: "No boundary condition"
@@ -345,7 +345,7 @@ KEYWORD ml_curv_bound "Boundary condition to implement"
 } "skip"
 
 private:
-KEYWORD ml_lapse_bound "Boundary condition to implement"
+KEYWORD ML_lapse_bound "Boundary condition to implement"
 {
   "flat" :: "Flat boundary condition"
   "none" :: "No boundary condition"
@@ -357,7 +357,7 @@ KEYWORD ml_lapse_bound "Boundary condition to implement"
 } "skip"
 
 private:
-KEYWORD ml_metric_bound "Boundary condition to implement"
+KEYWORD ML_metric_bound "Boundary condition to implement"
 {
   "flat" :: "Flat boundary condition"
   "none" :: "No boundary condition"
@@ -369,7 +369,7 @@ KEYWORD ml_metric_bound "Boundary condition to implement"
 } "skip"
 
 private:
-KEYWORD ml_shift_bound "Boundary condition to implement"
+KEYWORD ML_shift_bound "Boundary condition to implement"
 {
   "flat" :: "Flat boundary condition"
   "none" :: "No boundary condition"
@@ -477,25 +477,25 @@ CCTK_REAL beta3_bound_speed "characteristic speed at boundary"
 } 1.
 
 private:
-CCTK_REAL ml_curv_bound_speed "characteristic speed at boundary"
+CCTK_REAL ML_curv_bound_speed "characteristic speed at boundary"
 {
   "0:*" :: "outgoing characteristic speed > 0"
 } 1.
 
 private:
-CCTK_REAL ml_lapse_bound_speed "characteristic speed at boundary"
+CCTK_REAL ML_lapse_bound_speed "characteristic speed at boundary"
 {
   "0:*" :: "outgoing characteristic speed > 0"
 } 1.
 
 private:
-CCTK_REAL ml_metric_bound_speed "characteristic speed at boundary"
+CCTK_REAL ML_metric_bound_speed "characteristic speed at boundary"
 {
   "0:*" :: "outgoing characteristic speed > 0"
 } 1.
 
 private:
-CCTK_REAL ml_shift_bound_speed "characteristic speed at boundary"
+CCTK_REAL ML_shift_bound_speed "characteristic speed at boundary"
 {
   "0:*" :: "outgoing characteristic speed > 0"
 } 1.
@@ -597,25 +597,25 @@ CCTK_REAL beta3_bound_limit "limit value for r -> infinity"
 } 0.
 
 private:
-CCTK_REAL ml_curv_bound_limit "limit value for r -> infinity"
+CCTK_REAL ML_curv_bound_limit "limit value for r -> infinity"
 {
   "*:*" :: "value of limit value is unrestricted"
 } 0.
 
 private:
-CCTK_REAL ml_lapse_bound_limit "limit value for r -> infinity"
+CCTK_REAL ML_lapse_bound_limit "limit value for r -> infinity"
 {
   "*:*" :: "value of limit value is unrestricted"
 } 0.
 
 private:
-CCTK_REAL ml_metric_bound_limit "limit value for r -> infinity"
+CCTK_REAL ML_metric_bound_limit "limit value for r -> infinity"
 {
   "*:*" :: "value of limit value is unrestricted"
 } 0.
 
 private:
-CCTK_REAL ml_shift_bound_limit "limit value for r -> infinity"
+CCTK_REAL ML_shift_bound_limit "limit value for r -> infinity"
 {
   "*:*" :: "value of limit value is unrestricted"
 } 0.
@@ -717,25 +717,25 @@ CCTK_REAL beta3_bound_scalar "Dirichlet boundary value"
 } 0.
 
 private:
-CCTK_REAL ml_curv_bound_scalar "Dirichlet boundary value"
+CCTK_REAL ML_curv_bound_scalar "Dirichlet boundary value"
 {
   "*:*" :: "unrestricted"
 } 0.
 
 private:
-CCTK_REAL ml_lapse_bound_scalar "Dirichlet boundary value"
+CCTK_REAL ML_lapse_bound_scalar "Dirichlet boundary value"
 {
   "*:*" :: "unrestricted"
 } 0.
 
 private:
-CCTK_REAL ml_metric_bound_scalar "Dirichlet boundary value"
+CCTK_REAL ML_metric_bound_scalar "Dirichlet boundary value"
 {
   "*:*" :: "unrestricted"
 } 0.
 
 private:
-CCTK_REAL ml_shift_bound_scalar "Dirichlet boundary value"
+CCTK_REAL ML_shift_bound_scalar "Dirichlet boundary value"
 {
   "*:*" :: "unrestricted"
 } 0.
diff --git a/ML_ADM/schedule.ccl b/ML_ADM/schedule.ccl
index e66f883..88c07e5 100644
--- a/ML_ADM/schedule.ccl
+++ b/ML_ADM/schedule.ccl
@@ -4,68 +4,68 @@
 # Mathematica script written by Ian Hinder and Sascha Husa
 
 
-STORAGE: Ham[1]
+STORAGE: ML_Ham[1]
 
-STORAGE: mom[1]
+STORAGE: ML_mom[1]
 
-STORAGE: ml_curvrhs[1]
+STORAGE: ML_curvrhs[1]
 
-STORAGE: ml_lapserhs[1]
+STORAGE: ML_lapserhs[1]
 
-STORAGE: ml_metricrhs[1]
+STORAGE: ML_metricrhs[1]
 
-STORAGE: ml_shiftrhs[1]
+STORAGE: ML_shiftrhs[1]
 
 if (timelevels == 1)
 {
-  STORAGE: ml_curv[1]
+  STORAGE: ML_curv[1]
 }
 if (timelevels == 2)
 {
-  STORAGE: ml_curv[2]
+  STORAGE: ML_curv[2]
 }
 if (timelevels == 3)
 {
-  STORAGE: ml_curv[3]
+  STORAGE: ML_curv[3]
 }
 
 if (timelevels == 1)
 {
-  STORAGE: ml_lapse[1]
+  STORAGE: ML_lapse[1]
 }
 if (timelevels == 2)
 {
-  STORAGE: ml_lapse[2]
+  STORAGE: ML_lapse[2]
 }
 if (timelevels == 3)
 {
-  STORAGE: ml_lapse[3]
+  STORAGE: ML_lapse[3]
 }
 
 if (timelevels == 1)
 {
-  STORAGE: ml_metric[1]
+  STORAGE: ML_metric[1]
 }
 if (timelevels == 2)
 {
-  STORAGE: ml_metric[2]
+  STORAGE: ML_metric[2]
 }
 if (timelevels == 3)
 {
-  STORAGE: ml_metric[3]
+  STORAGE: ML_metric[3]
 }
 
 if (timelevels == 1)
 {
-  STORAGE: ml_shift[1]
+  STORAGE: ML_shift[1]
 }
 if (timelevels == 2)
 {
-  STORAGE: ml_shift[2]
+  STORAGE: ML_shift[2]
 }
 if (timelevels == 3)
 {
-  STORAGE: ml_shift[3]
+  STORAGE: ML_shift[3]
 }
 
 schedule ML_ADM_Startup at STARTUP
@@ -112,10 +112,10 @@ schedule ML_ADM_RHS IN MoL_CalcRHS
 schedule ML_ADM_RHS AT analysis
 {
   LANG: C
-  SYNC: ml_curvrhs
-  SYNC: ml_lapserhs
-  SYNC: ml_metricrhs
-  SYNC: ml_shiftrhs
+  SYNC: ML_curvrhs
+  SYNC: ML_lapserhs
+  SYNC: ML_metricrhs
+  SYNC: ML_shiftrhs
 } "ML_ADM_RHS"
 
 
@@ -127,7 +127,7 @@ if (CCTK_EQUALS(my_boundary_condition, "Minkowski"))
   } "ML_ADM_boundary"
 }
 
-schedule ML_ADM_convertToADMBase IN MoL_PostStep AFTER (ML_ADM_ApplyBCs ML_ADM_boundary)
+schedule ML_ADM_convertToADMBase IN MoL_PostStep AFTER ML_ADM_ApplyBCs
 {
   LANG: C
 } "ML_ADM_convertToADMBase"
@@ -135,10 +135,10 @@ schedule ML_ADM_convertToADMBase IN MoL_PostStep AFTER (ML_ADM_ApplyBCs ML_ADM_b
 schedule ML_ADM_constraints AT analysis
 {
   LANG: C
-  SYNC: Ham
-  SYNC: mom
-  TRIGGERS: Ham
-  TRIGGERS: mom
+  SYNC: ML_Ham
+  SYNC: ML_mom
+  TRIGGERS: ML_Ham
+  TRIGGERS: ML_mom
 } "ML_ADM_constraints"
 
 schedule ML_ADM_constraints_boundary AT analysis AFTER ML_ADM_constraints
@@ -150,10 +150,10 @@ schedule ML_ADM_ApplyBoundConds in MoL_PostStep
 {
   LANG: C
   OPTIONS: level
-  SYNC: ml_curv
-  SYNC: ml_lapse
-  SYNC: ml_metric
-  SYNC: ml_shift
+  SYNC: ML_curv
+  SYNC: ML_lapse
+  SYNC: ML_metric
+  SYNC: ML_shift
 } "apply boundary conditions"
 
 schedule ML_ADM_CheckBoundaries at BASEGRID
diff --git a/ML_ADM/src/Boundaries.c b/ML_ADM/src/Boundaries.c
index 8c6c288..de39fc6 100644
--- a/ML_ADM/src/Boundaries.c
+++ b/ML_ADM/src/Boundaries.c
@@ -35,48 +35,48 @@ void ML_ADM_ApplyBoundConds(CCTK_ARGUMENTS)
   
   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"  ) ) 
+  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);
+                      "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!");
+       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"  ) ) 
+  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);
+                      "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!");
+       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"  ) ) 
+  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);
+                      "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!");
+       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"  ) ) 
+  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);
+                      "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!");
+       CCTK_WARN(0, "Failed to register ML_shift_bound BC for ML_ADM::ML_shift!");
   }
   
   if (CCTK_EQUALS(K11_bound, "none"  ) ||
@@ -255,79 +255,79 @@ void ML_ADM_ApplyBoundConds(CCTK_ARGUMENTS)
        CCTK_WARN(0, "Failed to register beta3_bound BC for ML_ADM::beta3!");
   }
   
-  if (CCTK_EQUALS(ml_curv_bound, "radiative"))
+  if (CCTK_EQUALS(ML_curv_bound, "radiative"))
   {
    /* apply 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)
+    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)
+    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");
+    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!");
+       CCTK_WARN(0, "Failed to register Radiation BC for ML_ADM::ML_curv!");
   
   }
   
-  if (CCTK_EQUALS(ml_lapse_bound, "radiative"))
+  if (CCTK_EQUALS(ML_lapse_bound, "radiative"))
   {
    /* apply 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)
+    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)
+    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");
+    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!");
+       CCTK_WARN(0, "Failed to register Radiation BC for ML_ADM::ML_lapse!");
   
   }
   
-  if (CCTK_EQUALS(ml_metric_bound, "radiative"))
+  if (CCTK_EQUALS(ML_metric_bound, "radiative"))
   {
    /* apply 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)
+    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)
+    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");
+    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!");
+       CCTK_WARN(0, "Failed to register Radiation BC for ML_ADM::ML_metric!");
   
   }
   
-  if (CCTK_EQUALS(ml_shift_bound, "radiative"))
+  if (CCTK_EQUALS(ML_shift_bound, "radiative"))
   {
    /* apply 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)
+    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)
+    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");
+    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!");
+       CCTK_WARN(0, "Failed to register Radiation BC for ML_ADM::ML_shift!");
   
   }
   
@@ -635,71 +635,71 @@ void ML_ADM_ApplyBoundConds(CCTK_ARGUMENTS)
   
   }
   
-  if (CCTK_EQUALS(ml_curv_bound, "scalar"))
+  if (CCTK_EQUALS(ML_curv_bound, "scalar"))
   {
    /* apply 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)
+    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");
+    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!");
+       CCTK_WARN(0, "Failed to register Scalar BC for ML_ADM::ML_curv!");
   
   }
   
-  if (CCTK_EQUALS(ml_lapse_bound, "scalar"))
+  if (CCTK_EQUALS(ML_lapse_bound, "scalar"))
   {
    /* apply 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)
+    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");
+    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!");
+       CCTK_WARN(0, "Failed to register Scalar BC for ML_ADM::ML_lapse!");
   
   }
   
-  if (CCTK_EQUALS(ml_metric_bound, "scalar"))
+  if (CCTK_EQUALS(ML_metric_bound, "scalar"))
   {
    /* apply 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)
+    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");
+    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!");
+       CCTK_WARN(0, "Failed to register Scalar BC for ML_ADM::ML_metric!");
   
   }
   
-  if (CCTK_EQUALS(ml_shift_bound, "scalar"))
+  if (CCTK_EQUALS(ML_shift_bound, "scalar"))
   {
    /* apply 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)
+    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");
+    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!");
+       CCTK_WARN(0, "Failed to register Scalar BC for ML_ADM::ML_shift!");
   
   }
   
@@ -980,25 +980,25 @@ void ML_ADM_ApplyBoundConds(CCTK_ARGUMENTS)
 
 
 /* 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_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_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_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::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
diff --git a/ML_ADM/src/ML_ADM_Minkowski.c b/ML_ADM/src/ML_ADM_Minkowski.c
index dd11e04..8912c16 100644
--- a/ML_ADM/src/ML_ADM_Minkowski.c
+++ b/ML_ADM/src/ML_ADM_Minkowski.c
@@ -5,7 +5,10 @@
 
 #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"
@@ -84,7 +87,7 @@ void ML_ADM_Minkowski_Body(cGH *cctkGH, CCTK_INT dir, CCTK_INT face, CCTK_REAL n
   pm1o12dz2 = -pow(dz,-2)/12.;
   
   /* Loop over the grid points */
-  _Pragma ("omp parallel")
+  #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])
diff --git a/ML_ADM/src/ML_ADM_RHS.c b/ML_ADM/src/ML_ADM_RHS.c
index 424275d..936ed3f 100644
--- a/ML_ADM/src/ML_ADM_RHS.c
+++ b/ML_ADM/src/ML_ADM_RHS.c
@@ -5,7 +5,10 @@
 
 #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"
@@ -84,7 +87,7 @@ void ML_ADM_RHS_Body(cGH *cctkGH, CCTK_INT dir, CCTK_INT face, CCTK_REAL normal[
   pm1o12dz2 = -pow(dz,-2)/12.;
   
   /* Loop over the grid points */
-  _Pragma ("omp parallel")
+  #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])
@@ -99,7 +102,8 @@ void ML_ADM_RHS_Body(cGH *cctkGH, CCTK_INT dir, CCTK_INT face, CCTK_REAL normal[
     CCTK_REAL G111 = INITVALUE, G112 = INITVALUE, G113 = INITVALUE, G122 = INITVALUE, G123 = INITVALUE, G133 = INITVALUE;
     CCTK_REAL G211 = INITVALUE, G212 = INITVALUE, G213 = INITVALUE, G222 = INITVALUE, G223 = INITVALUE, G233 = INITVALUE;
     CCTK_REAL G311 = INITVALUE, G312 = INITVALUE, G313 = INITVALUE, G322 = INITVALUE, G323 = INITVALUE, G333 = INITVALUE;
-    CCTK_REAL gu11 = INITVALUE, gu21 = INITVALUE, gu22 = INITVALUE, gu31 = INITVALUE, gu32 = INITVALUE, gu33 = INITVALUE;
+    CCTK_REAL gu11 = INITVALUE, gu12 = INITVALUE, gu13 = INITVALUE, gu21 = INITVALUE, gu22 = INITVALUE, gu23 = INITVALUE;
+    CCTK_REAL gu31 = INITVALUE, gu32 = INITVALUE, gu33 = INITVALUE;
     CCTK_REAL Km11 = INITVALUE, Km12 = INITVALUE, Km13 = INITVALUE, Km21 = INITVALUE, Km22 = INITVALUE, Km23 = INITVALUE;
     CCTK_REAL Km31 = INITVALUE, Km32 = INITVALUE, Km33 = INITVALUE;
     CCTK_REAL R11 = INITVALUE, R12 = INITVALUE, R13 = INITVALUE, R22 = INITVALUE, R23 = INITVALUE, R33 = INITVALUE;
@@ -136,31 +140,30 @@ void ML_ADM_RHS_Body(cGH *cctkGH, CCTK_INT dir, CCTK_INT face, CCTK_REAL normal[
     CCTK_REAL PDstandardNth1g11 = INITVALUE;
     CCTK_REAL PDstandardNth2g11 = INITVALUE;
     CCTK_REAL PDstandardNth3g11 = INITVALUE;
+    CCTK_REAL PDstandardNth11g11 = INITVALUE;
     CCTK_REAL PDstandardNth22g11 = INITVALUE;
     CCTK_REAL PDstandardNth33g11 = INITVALUE;
     CCTK_REAL PDstandardNth12g11 = INITVALUE;
     CCTK_REAL PDstandardNth13g11 = INITVALUE;
-    CCTK_REAL PDstandardNth21g11 = INITVALUE;
     CCTK_REAL PDstandardNth23g11 = INITVALUE;
-    CCTK_REAL PDstandardNth31g11 = INITVALUE;
-    CCTK_REAL PDstandardNth32g11 = INITVALUE;
     CCTK_REAL PDstandardNth1g12 = INITVALUE;
     CCTK_REAL PDstandardNth2g12 = INITVALUE;
     CCTK_REAL PDstandardNth3g12 = INITVALUE;
+    CCTK_REAL PDstandardNth11g12 = INITVALUE;
+    CCTK_REAL PDstandardNth22g12 = INITVALUE;
     CCTK_REAL PDstandardNth33g12 = INITVALUE;
     CCTK_REAL PDstandardNth12g12 = INITVALUE;
     CCTK_REAL PDstandardNth13g12 = INITVALUE;
     CCTK_REAL PDstandardNth21g12 = INITVALUE;
     CCTK_REAL PDstandardNth23g12 = INITVALUE;
-    CCTK_REAL PDstandardNth31g12 = INITVALUE;
-    CCTK_REAL PDstandardNth32g12 = INITVALUE;
     CCTK_REAL PDstandardNth1g13 = INITVALUE;
     CCTK_REAL PDstandardNth2g13 = INITVALUE;
     CCTK_REAL PDstandardNth3g13 = INITVALUE;
+    CCTK_REAL PDstandardNth11g13 = INITVALUE;
     CCTK_REAL PDstandardNth22g13 = INITVALUE;
+    CCTK_REAL PDstandardNth33g13 = INITVALUE;
     CCTK_REAL PDstandardNth12g13 = INITVALUE;
     CCTK_REAL PDstandardNth13g13 = INITVALUE;
-    CCTK_REAL PDstandardNth21g13 = INITVALUE;
     CCTK_REAL PDstandardNth23g13 = INITVALUE;
     CCTK_REAL PDstandardNth31g13 = INITVALUE;
     CCTK_REAL PDstandardNth32g13 = INITVALUE;
@@ -168,17 +171,18 @@ void ML_ADM_RHS_Body(cGH *cctkGH, CCTK_INT dir, CCTK_INT face, CCTK_REAL normal[
     CCTK_REAL PDstandardNth2g22 = INITVALUE;
     CCTK_REAL PDstandardNth3g22 = INITVALUE;
     CCTK_REAL PDstandardNth11g22 = INITVALUE;
+    CCTK_REAL PDstandardNth22g22 = INITVALUE;
     CCTK_REAL PDstandardNth33g22 = INITVALUE;
     CCTK_REAL PDstandardNth12g22 = INITVALUE;
     CCTK_REAL PDstandardNth13g22 = INITVALUE;
     CCTK_REAL PDstandardNth21g22 = INITVALUE;
     CCTK_REAL PDstandardNth23g22 = INITVALUE;
-    CCTK_REAL PDstandardNth31g22 = INITVALUE;
-    CCTK_REAL PDstandardNth32g22 = INITVALUE;
     CCTK_REAL PDstandardNth1g23 = INITVALUE;
     CCTK_REAL PDstandardNth2g23 = INITVALUE;
     CCTK_REAL PDstandardNth3g23 = INITVALUE;
     CCTK_REAL PDstandardNth11g23 = INITVALUE;
+    CCTK_REAL PDstandardNth22g23 = INITVALUE;
+    CCTK_REAL PDstandardNth33g23 = INITVALUE;
     CCTK_REAL PDstandardNth12g23 = INITVALUE;
     CCTK_REAL PDstandardNth13g23 = INITVALUE;
     CCTK_REAL PDstandardNth21g23 = INITVALUE;
@@ -190,29 +194,20 @@ void ML_ADM_RHS_Body(cGH *cctkGH, CCTK_INT dir, CCTK_INT face, CCTK_REAL normal[
     CCTK_REAL PDstandardNth3g33 = INITVALUE;
     CCTK_REAL PDstandardNth11g33 = INITVALUE;
     CCTK_REAL PDstandardNth22g33 = INITVALUE;
+    CCTK_REAL PDstandardNth33g33 = INITVALUE;
     CCTK_REAL PDstandardNth12g33 = INITVALUE;
     CCTK_REAL PDstandardNth13g33 = INITVALUE;
-    CCTK_REAL PDstandardNth21g33 = INITVALUE;
     CCTK_REAL PDstandardNth23g33 = INITVALUE;
     CCTK_REAL PDstandardNth31g33 = INITVALUE;
     CCTK_REAL PDstandardNth32g33 = INITVALUE;
     CCTK_REAL PDstandardNth1K11 = INITVALUE;
-    CCTK_REAL PDstandardNth2K11 = INITVALUE;
-    CCTK_REAL PDstandardNth3K11 = INITVALUE;
     CCTK_REAL PDstandardNth1K12 = INITVALUE;
     CCTK_REAL PDstandardNth2K12 = INITVALUE;
-    CCTK_REAL PDstandardNth3K12 = INITVALUE;
     CCTK_REAL PDstandardNth1K13 = INITVALUE;
-    CCTK_REAL PDstandardNth2K13 = INITVALUE;
     CCTK_REAL PDstandardNth3K13 = INITVALUE;
-    CCTK_REAL PDstandardNth1K22 = INITVALUE;
     CCTK_REAL PDstandardNth2K22 = INITVALUE;
-    CCTK_REAL PDstandardNth3K22 = INITVALUE;
-    CCTK_REAL PDstandardNth1K23 = INITVALUE;
     CCTK_REAL PDstandardNth2K23 = INITVALUE;
     CCTK_REAL PDstandardNth3K23 = INITVALUE;
-    CCTK_REAL PDstandardNth1K33 = INITVALUE;
-    CCTK_REAL PDstandardNth2K33 = INITVALUE;
     CCTK_REAL PDstandardNth3K33 = INITVALUE;
     
     /* Assign local copies of grid functions */
@@ -259,12 +254,17 @@ void ML_ADM_RHS_Body(cGH *cctkGH, CCTK_INT dir, CCTK_INT face, CCTK_REAL normal[
     PDstandardNth1g11 = PDstandardNth1(g11, i, j, k);
     PDstandardNth2g11 = PDstandardNth2(g11, i, j, k);
     PDstandardNth3g11 = PDstandardNth3(g11, i, j, k);
+    PDstandardNth11g11 = PDstandardNth11(g11, i, j, k);
     PDstandardNth22g11 = PDstandardNth22(g11, i, j, k);
     PDstandardNth33g11 = PDstandardNth33(g11, i, j, k);
+    PDstandardNth12g11 = PDstandardNth12(g11, i, j, k);
+    PDstandardNth13g11 = PDstandardNth13(g11, i, j, k);
     PDstandardNth23g11 = PDstandardNth23(g11, i, j, k);
     PDstandardNth1g12 = PDstandardNth1(g12, i, j, k);
     PDstandardNth2g12 = PDstandardNth2(g12, i, j, k);
     PDstandardNth3g12 = PDstandardNth3(g12, i, j, k);
+    PDstandardNth11g12 = PDstandardNth11(g12, i, j, k);
+    PDstandardNth22g12 = PDstandardNth22(g12, i, j, k);
     PDstandardNth33g12 = PDstandardNth33(g12, i, j, k);
     PDstandardNth12g12 = PDstandardNth12(g12, i, j, k);
     PDstandardNth13g12 = PDstandardNth13(g12, i, j, k);
@@ -272,7 +272,9 @@ void ML_ADM_RHS_Body(cGH *cctkGH, CCTK_INT dir, CCTK_INT face, CCTK_REAL normal[
     PDstandardNth1g13 = PDstandardNth1(g13, i, j, k);
     PDstandardNth2g13 = PDstandardNth2(g13, i, j, k);
     PDstandardNth3g13 = PDstandardNth3(g13, i, j, k);
+    PDstandardNth11g13 = PDstandardNth11(g13, i, j, k);
     PDstandardNth22g13 = PDstandardNth22(g13, i, j, k);
+    PDstandardNth33g13 = PDstandardNth33(g13, i, j, k);
     PDstandardNth12g13 = PDstandardNth12(g13, i, j, k);
     PDstandardNth13g13 = PDstandardNth13(g13, i, j, k);
     PDstandardNth23g13 = PDstandardNth23(g13, i, j, k);
@@ -280,12 +282,17 @@ void ML_ADM_RHS_Body(cGH *cctkGH, CCTK_INT dir, CCTK_INT face, CCTK_REAL normal[
     PDstandardNth2g22 = PDstandardNth2(g22, i, j, k);
     PDstandardNth3g22 = PDstandardNth3(g22, i, j, k);
     PDstandardNth11g22 = PDstandardNth11(g22, i, j, k);
+    PDstandardNth22g22 = PDstandardNth22(g22, i, j, k);
     PDstandardNth33g22 = PDstandardNth33(g22, i, j, k);
+    PDstandardNth12g22 = PDstandardNth12(g22, i, j, k);
     PDstandardNth13g22 = PDstandardNth13(g22, i, j, k);
+    PDstandardNth23g22 = PDstandardNth23(g22, i, j, k);
     PDstandardNth1g23 = PDstandardNth1(g23, i, j, k);
     PDstandardNth2g23 = PDstandardNth2(g23, i, j, k);
     PDstandardNth3g23 = PDstandardNth3(g23, i, j, k);
     PDstandardNth11g23 = PDstandardNth11(g23, i, j, k);
+    PDstandardNth22g23 = PDstandardNth22(g23, i, j, k);
+    PDstandardNth33g23 = PDstandardNth33(g23, i, j, k);
     PDstandardNth12g23 = PDstandardNth12(g23, i, j, k);
     PDstandardNth13g23 = PDstandardNth13(g23, i, j, k);
     PDstandardNth23g23 = PDstandardNth23(g23, i, j, k);
@@ -294,206 +301,297 @@ void ML_ADM_RHS_Body(cGH *cctkGH, CCTK_INT dir, CCTK_INT face, CCTK_REAL normal[
     PDstandardNth3g33 = PDstandardNth3(g33, i, j, k);
     PDstandardNth11g33 = PDstandardNth11(g33, i, j, k);
     PDstandardNth22g33 = PDstandardNth22(g33, i, j, k);
+    PDstandardNth33g33 = PDstandardNth33(g33, i, j, k);
     PDstandardNth12g33 = PDstandardNth12(g33, i, j, k);
+    PDstandardNth13g33 = PDstandardNth13(g33, i, j, k);
+    PDstandardNth23g33 = PDstandardNth23(g33, i, j, k);
     PDstandardNth1K11 = PDstandardNth1(K11, i, j, k);
-    PDstandardNth2K11 = PDstandardNth2(K11, i, j, k);
-    PDstandardNth3K11 = PDstandardNth3(K11, i, j, k);
     PDstandardNth1K12 = PDstandardNth1(K12, i, j, k);
     PDstandardNth2K12 = PDstandardNth2(K12, i, j, k);
-    PDstandardNth3K12 = PDstandardNth3(K12, i, j, k);
     PDstandardNth1K13 = PDstandardNth1(K13, i, j, k);
-    PDstandardNth2K13 = PDstandardNth2(K13, i, j, k);
     PDstandardNth3K13 = PDstandardNth3(K13, i, j, k);
-    PDstandardNth1K22 = PDstandardNth1(K22, i, j, k);
     PDstandardNth2K22 = PDstandardNth2(K22, i, j, k);
-    PDstandardNth3K22 = PDstandardNth3(K22, i, j, k);
-    PDstandardNth1K23 = PDstandardNth1(K23, i, j, k);
     PDstandardNth2K23 = PDstandardNth2(K23, i, j, k);
     PDstandardNth3K23 = PDstandardNth3(K23, i, j, k);
-    PDstandardNth1K33 = PDstandardNth1(K33, i, j, k);
-    PDstandardNth2K33 = PDstandardNth2(K33, i, j, k);
     PDstandardNth3K33 = PDstandardNth3(K33, i, j, k);
     
     /* Precompute derivatives (old style) */
     
     /* Calculate temporaries and grid functions */
-    detg  =  2*g12L*g13L*g23L + g33L*(g11L*g22L - SQR(g12L)) - g22L*SQR(g13L) - g11L*SQR(g23L);
+    detg  =  2*g12L*g13L*g23L + g33L*(g11L*g22L - SQR(g12L)) - g22L*SQR(g13L) - 
+        g11L*SQR(g23L);
     
     gu11  =  INV(detg)*(g22L*g33L - SQR(g23L));
     
-    gu21  =  (g13L*g23L - g12L*g33L)*INV(detg);
+    gu12  =  (g13L*g23L - g12L*g33L)*INV(detg);
     
-    gu31  =  (-(g13L*g22L) + g12L*g23L)*INV(detg);
+    gu13  =  (-(g13L*g22L) + g12L*g23L)*INV(detg);
+    
+    gu21  =  (g13L*g23L - g12L*g33L)*INV(detg);
     
     gu22  =  INV(detg)*(g11L*g33L - SQR(g13L));
     
+    gu23  =  (g12L*g13L - g11L*g23L)*INV(detg);
+    
+    gu31  =  (-(g13L*g22L) + g12L*g23L)*INV(detg);
+    
     gu32  =  (g12L*g13L - g11L*g23L)*INV(detg);
     
     gu33  =  INV(detg)*(g11L*g22L - SQR(g12L));
     
-    G111  =  khalf*(gu11*PDstandardNth1g11 + 2*(gu21*PDstandardNth1g12 + gu31*PDstandardNth1g13) - gu21*PDstandardNth2g11 - 
-          gu31*PDstandardNth3g11);
+    G111  =  khalf*(gu11*PDstandardNth1g11 + 
+          2*(gu12*PDstandardNth1g12 + gu13*PDstandardNth1g13) - 
+          gu12*PDstandardNth2g11 - gu13*PDstandardNth3g11);
     
-    G211  =  khalf*(gu21*PDstandardNth1g11 + 2*(gu22*PDstandardNth1g12 + gu32*PDstandardNth1g13) - gu22*PDstandardNth2g11 - 
-          gu32*PDstandardNth3g11);
+    G211  =  khalf*(gu21*PDstandardNth1g11 + 
+          2*(gu22*PDstandardNth1g12 + gu23*PDstandardNth1g13) - 
+          gu22*PDstandardNth2g11 - gu23*PDstandardNth3g11);
     
-    G311  =  khalf*(gu31*PDstandardNth1g11 + 2*(gu32*PDstandardNth1g12 + gu33*PDstandardNth1g13) - gu32*PDstandardNth2g11 - 
-          gu33*PDstandardNth3g11);
+    G311  =  khalf*(gu31*PDstandardNth1g11 + 
+          2*(gu32*PDstandardNth1g12 + gu33*PDstandardNth1g13) - 
+          gu32*PDstandardNth2g11 - gu33*PDstandardNth3g11);
     
-    G112  =  khalf*(gu21*PDstandardNth1g22 + gu11*PDstandardNth2g11 + 
-          gu31*(PDstandardNth1g23 + PDstandardNth2g13 - PDstandardNth3g12));
+    G112  =  khalf*(gu12*PDstandardNth1g22 + gu11*PDstandardNth2g11 + 
+          gu13*(PDstandardNth1g23 + PDstandardNth2g13 - PDstandardNth3g12));
     
     G212  =  khalf*(gu22*PDstandardNth1g22 + gu21*PDstandardNth2g11 + 
-          gu32*(PDstandardNth1g23 + PDstandardNth2g13 - PDstandardNth3g12));
+          gu23*(PDstandardNth1g23 + PDstandardNth2g13 - PDstandardNth3g12));
     
     G312  =  khalf*(gu32*PDstandardNth1g22 + gu31*PDstandardNth2g11 + 
           gu33*(PDstandardNth1g23 + PDstandardNth2g13 - PDstandardNth3g12));
     
-    G113  =  khalf*(gu31*PDstandardNth1g33 + gu11*PDstandardNth3g11 + 
-          gu21*(PDstandardNth1g23 - PDstandardNth2g13 + PDstandardNth3g12));
+    G113  =  khalf*(gu13*PDstandardNth1g33 + gu11*PDstandardNth3g11 + 
+          gu12*(PDstandardNth1g23 - PDstandardNth2g13 + PDstandardNth3g12));
     
-    G213  =  khalf*(gu32*PDstandardNth1g33 + gu21*PDstandardNth3g11 + 
+    G213  =  khalf*(gu23*PDstandardNth1g33 + gu21*PDstandardNth3g11 + 
           gu22*(PDstandardNth1g23 - PDstandardNth2g13 + PDstandardNth3g12));
     
     G313  =  khalf*(gu33*PDstandardNth1g33 + gu31*PDstandardNth3g11 + 
           gu32*(PDstandardNth1g23 - PDstandardNth2g13 + PDstandardNth3g12));
     
-    G122  =  khalf*(gu11*(-PDstandardNth1g22 + 2*PDstandardNth2g12) + gu21*PDstandardNth2g22 + 
-          gu31*(2*PDstandardNth2g23 - PDstandardNth3g22));
+    G122  =  khalf*(gu11*(-PDstandardNth1g22 + 2*PDstandardNth2g12) + 
+          gu12*PDstandardNth2g22 + gu13*(2*PDstandardNth2g23 - PDstandardNth3g22));
     
-    G222  =  khalf*(gu21*(-PDstandardNth1g22 + 2*PDstandardNth2g12) + gu22*PDstandardNth2g22 + 
-          gu32*(2*PDstandardNth2g23 - PDstandardNth3g22));
+    G222  =  khalf*(gu21*(-PDstandardNth1g22 + 2*PDstandardNth2g12) + 
+          gu22*PDstandardNth2g22 + gu23*(2*PDstandardNth2g23 - PDstandardNth3g22));
     
-    G322  =  khalf*(gu31*(-PDstandardNth1g22 + 2*PDstandardNth2g12) + gu32*PDstandardNth2g22 + 
-          gu33*(2*PDstandardNth2g23 - PDstandardNth3g22));
+    G322  =  khalf*(gu31*(-PDstandardNth1g22 + 2*PDstandardNth2g12) + 
+          gu32*PDstandardNth2g22 + gu33*(2*PDstandardNth2g23 - PDstandardNth3g22));
     
-    G123  =  khalf*(gu31*PDstandardNth2g33 + gu11*(-PDstandardNth1g23 + PDstandardNth2g13 + PDstandardNth3g12) + 
-          gu21*PDstandardNth3g22);
+    G123  =  khalf*(gu13*PDstandardNth2g33 + 
+          gu11*(-PDstandardNth1g23 + PDstandardNth2g13 + PDstandardNth3g12) + 
+          gu12*PDstandardNth3g22);
     
-    G223  =  khalf*(gu32*PDstandardNth2g33 + gu21*(-PDstandardNth1g23 + PDstandardNth2g13 + PDstandardNth3g12) + 
+    G223  =  khalf*(gu23*PDstandardNth2g33 + 
+          gu21*(-PDstandardNth1g23 + PDstandardNth2g13 + PDstandardNth3g12) + 
           gu22*PDstandardNth3g22);
     
-    G323  =  khalf*(gu33*PDstandardNth2g33 + gu31*(-PDstandardNth1g23 + PDstandardNth2g13 + PDstandardNth3g12) + 
+    G323  =  khalf*(gu33*PDstandardNth2g33 + 
+          gu31*(-PDstandardNth1g23 + PDstandardNth2g13 + PDstandardNth3g12) + 
           gu32*PDstandardNth3g22);
     
-    G133  =  khalf*(-(gu11*PDstandardNth1g33) - gu21*PDstandardNth2g33 + 2*gu11*PDstandardNth3g13 + 
-          2*gu21*PDstandardNth3g23 + gu31*PDstandardNth3g33);
-    
-    G233  =  khalf*(-(gu21*PDstandardNth1g33) - gu22*PDstandardNth2g33 + 2*gu21*PDstandardNth3g13 + 
-          2*gu22*PDstandardNth3g23 + gu32*PDstandardNth3g33);
-    
-    G333  =  khalf*(-(gu31*PDstandardNth1g33) - gu32*PDstandardNth2g33 + 2*gu31*PDstandardNth3g13 + 
-          2*gu32*PDstandardNth3g23 + gu33*PDstandardNth3g33);
-    
-    R11  =  khalf*(4*G213*G312 + G211*(2*G112 - 2*G222 - 2*G323) + G311*(2*G113 - 2*G333) - gu22*PDstandardNth11g22 - 
-          2*(G111*G212 + G223*G311 + G111*G313 + gu32*PDstandardNth11g23) - gu33*PDstandardNth11g33 + 
-          2*gu22*PDstandardNth12g12 + 2*gu32*PDstandardNth12g13 + 2*gu32*PDstandardNth13g12 + 2*gu33*PDstandardNth13g13 - 
-          gu22*PDstandardNth22g11 - 2*gu32*PDstandardNth23g11 - gu33*PDstandardNth33g11 + 2*SQR(G212) + 2*SQR(G313));
-    
-    R12  =  khalf*(2*(G122*G211 + G123*G311 + G213*G322 + G313*G323) - 
-          2*(G112*G212 + G112*G313 + G212*G323 + G312*G333 + gu21*PDstandardNth12g12) - gu32*PDstandardNth12g23 - 
-          gu33*PDstandardNth12g33 + gu31*(PDstandardNth11g23 - PDstandardNth12g13 - PDstandardNth13g12) + 
-          gu32*PDstandardNth13g22 + gu33*PDstandardNth13g23 + gu21*(PDstandardNth11g22 + PDstandardNth22g11) + 
-          gu32*PDstandardNth22g13 + gu31*PDstandardNth23g11 - gu32*PDstandardNth23g12 + gu33*PDstandardNth23g13 - 
-          gu33*PDstandardNth33g12);
-    
-    R13  =  khalf*(2*(G123*G211 + G212*G223 + G133*G311 + G233*G312) - 
-          2*(G213*G222 + G223*G313 + G113*(G212 + G313) + gu31*PDstandardNth13g13) + 
-          gu21*(PDstandardNth11g23 - PDstandardNth12g13 - PDstandardNth13g12 + PDstandardNth23g11) + 
-          gu22*(PDstandardNth12g23 - PDstandardNth13g22 - PDstandardNth22g13 + PDstandardNth23g12) + 
-          gu31*(PDstandardNth11g33 + PDstandardNth33g11) + 
-          gu32*(PDstandardNth12g33 - PDstandardNth13g23 - PDstandardNth23g13 + PDstandardNth33g12));
-    
-    R22  =  khalf*(4*G123*G312 + G122*(2*G212 - 2*(G111 + G313)) + G322*(2*G223 - 2*G333) - 
-          2*(G113*G322 + G222*(G112 + G323) + gu31*PDstandardNth13g22) + 
-          gu11*(-PDstandardNth11g22 + 2*PDstandardNth12g12 - PDstandardNth22g11) + 
-          gu31*(-2*PDstandardNth22g13 + 2*(PDstandardNth12g23 + PDstandardNth23g12)) + 
-          gu33*(-PDstandardNth22g33 + 2*PDstandardNth23g23 - PDstandardNth33g22) + 2*(SQR(G112) + SQR(G323)));
-    
-    R23  =  khalf*(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) + gu32*PDstandardNth23g23) + 
-          gu31*(PDstandardNth12g33 - PDstandardNth13g23 - PDstandardNth23g13 + PDstandardNth33g12) + 
-          gu32*(PDstandardNth22g33 + PDstandardNth33g22));
-    
-    R33  =  khalf*(4*G123*G213 - gu11*PDstandardNth11g33 - 
-          2*(G111*G133 + G133*G212 + G112*G233 + G222*G233 + G113*G333 + G223*G333 + gu21*PDstandardNth12g33) + 
-          2*(G133*G313 + G233*G323 + gu11*PDstandardNth13g13) + 2*gu21*PDstandardNth13g23 - gu22*PDstandardNth22g33 + 
-          2*gu21*PDstandardNth23g13 + 2*gu22*PDstandardNth23g23 - gu11*PDstandardNth33g11 - 2*gu21*PDstandardNth33g12 - 
-          gu22*PDstandardNth33g22 + 2*SQR(G113) + 2*SQR(G223));
-    
-    Km11  =  gu11*K11L + gu21*K12L + gu31*K13L;
-    
-    Km21  =  gu21*K11L + gu22*K12L + gu32*K13L;
+    G133  =  khalf*(-(gu11*PDstandardNth1g33) - gu12*PDstandardNth2g33 + 
+          2*gu11*PDstandardNth3g13 + 2*gu12*PDstandardNth3g23 + 
+          gu13*PDstandardNth3g33);
+    
+    G233  =  khalf*(-(gu21*PDstandardNth1g33) - gu22*PDstandardNth2g33 + 
+          2*gu21*PDstandardNth3g13 + 2*gu22*PDstandardNth3g23 + 
+          gu23*PDstandardNth3g33);
+    
+    G333  =  khalf*(-(gu31*PDstandardNth1g33) - gu32*PDstandardNth2g33 + 
+          2*gu31*PDstandardNth3g13 + 2*gu32*PDstandardNth3g23 + 
+          gu33*PDstandardNth3g33);
+    
+    R11  =  2*(G112*G211 + G113*G311 + G213*G312) - G111*(G111 + G212 + G313) - 
+        G211*(G112 + G222 + G323) - G311*(G113 + G223 + G333) + 
+        khalf*(-6*gu11*PDstandardNth11g11 + 
+           gu12*(-8*PDstandardNth11g12 + 2*PDstandardNth12g12) + 
+           gu22*(-9*PDstandardNth11g22 + 3*PDstandardNth12g22) + 
+           gu13*(-8*PDstandardNth11g13 + 2*PDstandardNth13g13) + 
+           gu33*(-9*PDstandardNth11g33 + 3*PDstandardNth13g33) + 
+           gu21*(-6*PDstandardNth11g12 - PDstandardNth12g12 + PDstandardNth22g12) + 
+           gu23*(-9*PDstandardNth11g23 + 3*PDstandardNth13g23 + 
+              PDstandardNth22g23 - PDstandardNth23g23) + 
+           gu31*(-6*PDstandardNth11g13 - PDstandardNth13g13 + PDstandardNth33g13) + 
+           gu32*(-9*PDstandardNth11g23 + 3*PDstandardNth12g23 - 
+              PDstandardNth23g23 + PDstandardNth33g23)) + SQR(G111) + SQR(G212) + 
+        SQR(G313);
+    
+    R12  =  khalf*(-2*(G112*G212 + G112*G313 + G212*G323 + G312*G333) - 
+          6*(gu11*PDstandardNth12g11 + gu31*PDstandardNth12g13) - 
+          9*(gu22*PDstandardNth12g22 + gu23*PDstandardNth12g23 + 
+             gu32*PDstandardNth12g23 + gu33*PDstandardNth12g33) - 
+          gu31*PDstandardNth13g13 + gu21*
+           (-7*PDstandardNth12g12 + PDstandardNth22g12) + 
+          gu12*(PDstandardNth11g12 - 10*PDstandardNth12g12 + 
+             3*PDstandardNth22g12) + 3*gu22*PDstandardNth22g22 + 
+          gu23*PDstandardNth22g23 + 3*gu32*PDstandardNth22g23 + 
+          gu13*(PDstandardNth11g13 - 9*PDstandardNth12g13 - PDstandardNth13g13 + 
+             3*PDstandardNth23g13) - gu32*PDstandardNth23g23 + 
+          2*(G122*G211 + G123*G311 + G213*G322 + G313*G323 + 
+             gu23*PDstandardNth23g23) + 3*gu33*PDstandardNth23g33 + 
+          gu31*PDstandardNth33g13 + gu32*PDstandardNth33g23);
+    
+    R13  =  khalf*(-2*(G213*G222 + G223*G313 + G113*(G212 + G313)) - 
+          6*gu11*PDstandardNth13g11 - 
+          9*(gu22*PDstandardNth13g22 + (gu23 + gu32)*PDstandardNth13g23 + 
+             gu33*PDstandardNth13g33) + 
+          gu21*(-PDstandardNth12g12 - 6*PDstandardNth13g12 + PDstandardNth22g12) + 
+          gu12*(PDstandardNth11g12 - PDstandardNth12g12 - 9*PDstandardNth13g12 + 
+             3*PDstandardNth23g12) + 
+          gu23*(PDstandardNth22g23 - PDstandardNth23g23) + 
+          2*(G123*G211 + G212*G223 + G133*G311 + G233*G312 + 
+             gu32*PDstandardNth23g23) + 
+          gu31*(-7*PDstandardNth13g13 + PDstandardNth33g13) + 
+          gu13*(PDstandardNth11g13 - 10*PDstandardNth13g13 + 
+             3*PDstandardNth33g13) + gu32*PDstandardNth33g23 + 
+          3*(gu22*PDstandardNth23g22 + gu23*PDstandardNth33g23 + 
+             gu33*PDstandardNth33g33));
+    
+    R22  =  -(G122*(G111 + G212 + G313)) + 2*(G122*G212 + G123*G312 + G223*G322) - 
+        G222*(G112 + G222 + G323) - G322*(G113 + G223 + G333) + 
+        khalf*(3*gu11*(PDstandardNth12g11 - 3*PDstandardNth22g11) + 
+           gu21*(2*PDstandardNth12g12 - 8*PDstandardNth22g12) + 
+           gu12*(PDstandardNth11g12 - PDstandardNth12g12 - 6*PDstandardNth22g12) - 
+           6*gu22*PDstandardNth22g22 + 
+           gu13*(PDstandardNth11g13 - PDstandardNth13g13 - 9*PDstandardNth22g13 + 
+              3*PDstandardNth23g13) - 
+           gu23*(8*PDstandardNth22g23 - 2*PDstandardNth23g23) + 
+           gu33*(-9*PDstandardNth22g33 + 3*PDstandardNth23g33) + 
+           gu31*(3*PDstandardNth12g13 - PDstandardNth13g13 - 9*PDstandardNth22g13 + 
+              PDstandardNth33g13) + gu32*
+            (-6*PDstandardNth22g23 - PDstandardNth23g23 + PDstandardNth33g23)) + 
+        SQR(G112) + SQR(G222) + SQR(G323);
+    
+    R23  =  khalf*(-2*(G111*G123 + G113*G323 + G223*(G112 + G323)) + 
+          2*(G112*G113 + G122*G213 + G133*G312 + G233*G322 + 
+             gu31*PDstandardNth13g13) + 
+          gu11*(3*PDstandardNth13g11 - 9*PDstandardNth23g11) + 
+          gu21*(-PDstandardNth12g12 + 3*PDstandardNth13g12 + PDstandardNth22g12 - 
+             9*PDstandardNth23g12) + 
+          gu12*(PDstandardNth11g12 - PDstandardNth12g12 - 6*PDstandardNth23g12) - 
+          6*gu22*PDstandardNth23g22 - 
+          9*(gu31*PDstandardNth23g13 + gu33*PDstandardNth23g33) + 
+          gu31*PDstandardNth33g13 + gu13*
+           (PDstandardNth11g13 - PDstandardNth13g13 - 9*PDstandardNth23g13 + 
+             3*PDstandardNth33g13) + 
+          gu32*(-7*PDstandardNth23g23 + PDstandardNth33g23) + 
+          gu23*(PDstandardNth22g23 - 10*PDstandardNth23g23 + 
+             3*PDstandardNth33g23) + 3*gu33*PDstandardNth33g33);
+    
+    R33  =  -(G133*(G111 + G212 + G313)) + 2*(G123*G213 + G133*G313) + 2*G233*G323 - 
+        G233*(G112 + G222 + G323) - G333*(G113 + G223 + G333) + 
+        khalf*(3*gu11*(PDstandardNth13g11 - 3*PDstandardNth33g11) + 
+           gu21*(-PDstandardNth12g12 + 3*PDstandardNth13g12 + PDstandardNth22g12 - 
+              9*PDstandardNth33g12) + 
+           gu12*(PDstandardNth11g12 - PDstandardNth12g12 + 3*PDstandardNth23g12 - 
+              9*PDstandardNth33g12) + 
+           gu31*(2*PDstandardNth13g13 - 8*PDstandardNth33g13) + 
+           gu13*(PDstandardNth11g13 - PDstandardNth13g13 - 6*PDstandardNth33g13) + 
+           3*gu22*(PDstandardNth23g22 - 3*PDstandardNth33g22) + 
+           gu32*(2*PDstandardNth23g23 - 8*PDstandardNth33g23) + 
+           gu23*(PDstandardNth22g23 - PDstandardNth23g23 - 6*PDstandardNth33g23) - 
+           6*gu33*PDstandardNth33g33) + SQR(G113) + SQR(G223) + SQR(G333);
+    
+    Km11  =  gu11*K11L + gu12*K12L + gu13*K13L;
+    
+    Km21  =  gu21*K11L + gu22*K12L + gu23*K13L;
     
     Km31  =  gu31*K11L + gu32*K12L + gu33*K13L;
     
-    Km12  =  gu11*K12L + gu21*K22L + gu31*K23L;
+    Km12  =  gu11*K12L + gu12*K22L + gu13*K23L;
     
-    Km22  =  gu21*K12L + gu22*K22L + gu32*K23L;
+    Km22  =  gu21*K12L + gu22*K22L + gu23*K23L;
     
     Km32  =  gu31*K12L + gu32*K22L + gu33*K23L;
     
-    Km13  =  gu11*K13L + gu21*K23L + gu31*K33L;
+    Km13  =  gu11*K13L + gu12*K23L + gu13*K33L;
     
-    Km23  =  gu21*K13L + gu22*K23L + gu32*K33L;
+    Km23  =  gu21*K13L + gu22*K23L + gu23*K33L;
     
     Km33  =  gu31*K13L + gu32*K23L + gu33*K33L;
     
     trK  =  Km11 + Km22 + Km33;
     
-    g11rhsL  =  -2*alphaL*K11L + 2*(g11L*PDstandardNth1beta1 + g12L*PDstandardNth1beta2 + g13L*PDstandardNth1beta3) + 
-        beta1L*PDstandardNth1g11 + beta2L*PDstandardNth2g11 + beta3L*PDstandardNth3g11;
-    
-    g12rhsL  =  -2*alphaL*K12L + g22L*PDstandardNth1beta2 + g23L*PDstandardNth1beta3 + beta1L*PDstandardNth1g12 + 
-        g11L*PDstandardNth2beta1 + g12L*(PDstandardNth1beta1 + PDstandardNth2beta2) + g13L*PDstandardNth2beta3 + 
-        beta2L*PDstandardNth2g12 + beta3L*PDstandardNth3g12;
-    
-    g13rhsL  =  -2*alphaL*K13L + g23L*PDstandardNth1beta2 + g33L*PDstandardNth1beta3 + beta1L*PDstandardNth1g13 + 
-        beta2L*PDstandardNth2g13 + g11L*PDstandardNth3beta1 + g12L*PDstandardNth3beta2 + 
-        g13L*(PDstandardNth1beta1 + PDstandardNth3beta3) + beta3L*PDstandardNth3g13;
-    
-    g22rhsL  =  -2*alphaL*K22L + beta1L*PDstandardNth1g22 + 
-        2*(g12L*PDstandardNth2beta1 + g22L*PDstandardNth2beta2 + g23L*PDstandardNth2beta3) + beta2L*PDstandardNth2g22 + 
-        beta3L*PDstandardNth3g22;
-    
-    g23rhsL  =  -2*alphaL*K23L + beta1L*PDstandardNth1g23 + g13L*PDstandardNth2beta1 + g33L*PDstandardNth2beta3 + 
-        beta2L*PDstandardNth2g23 + g12L*PDstandardNth3beta1 + g22L*PDstandardNth3beta2 + 
-        g23L*(PDstandardNth2beta2 + PDstandardNth3beta3) + beta3L*PDstandardNth3g23;
-    
-    g33rhsL  =  -2*alphaL*K33L + beta1L*PDstandardNth1g33 + beta2L*PDstandardNth2g33 + 
-        2*(g13L*PDstandardNth3beta1 + g23L*PDstandardNth3beta2 + g33L*PDstandardNth3beta3) + beta3L*PDstandardNth3g33;
-    
-    K11rhsL  =  -PDstandardNth11alpha + G111*PDstandardNth1alpha + 
-        2*(K11L*PDstandardNth1beta1 + K12L*PDstandardNth1beta2 + K13L*PDstandardNth1beta3) + beta1L*PDstandardNth1K11 + 
-        G211*PDstandardNth2alpha + beta2L*PDstandardNth2K11 + G311*PDstandardNth3alpha + beta3L*PDstandardNth3K11 + 
+    g11rhsL  =  -2*alphaL*K11L + 4*(g11L*PDstandardNth1beta1 + 
+           g12L*PDstandardNth1beta2 + g13L*PDstandardNth1beta3) + 
+        beta1L*PDstandardNth1g11 + beta2L*PDstandardNth2g12 + 
+        beta3L*PDstandardNth3g13;
+    
+    g12rhsL  =  -2*alphaL*K12L + 3*(g12L*PDstandardNth1beta1 + 
+           g22L*PDstandardNth1beta2 + g23L*PDstandardNth1beta3) + 
+        beta1L*PDstandardNth1g12 + g12L*PDstandardNth2beta1 + 
+        g22L*PDstandardNth2beta2 + g23L*PDstandardNth2beta3 + 
+        beta2L*PDstandardNth2g22 + beta3L*PDstandardNth3g23;
+    
+    g13rhsL  =  -2*alphaL*K13L + 3*(g13L*PDstandardNth1beta1 + 
+           g23L*PDstandardNth1beta2 + g33L*PDstandardNth1beta3) + 
+        beta1L*PDstandardNth1g13 + beta2L*PDstandardNth2g23 + 
+        g13L*PDstandardNth3beta1 + g23L*PDstandardNth3beta2 + 
+        g33L*PDstandardNth3beta3 + beta3L*PDstandardNth3g33;
+    
+    g22rhsL  =  -2*alphaL*K22L + beta1L*PDstandardNth1g12 + 
+        4*(g12L*PDstandardNth2beta1 + g22L*PDstandardNth2beta2 + 
+           g23L*PDstandardNth2beta3) + beta2L*PDstandardNth2g22 + 
+        beta3L*PDstandardNth3g23;
+    
+    g23rhsL  =  -2*alphaL*K23L + beta1L*PDstandardNth1g13 + 
+        3*(g13L*PDstandardNth2beta1 + g23L*PDstandardNth2beta2 + 
+           g33L*PDstandardNth2beta3) + beta2L*PDstandardNth2g23 + 
+        g13L*PDstandardNth3beta1 + g23L*PDstandardNth3beta2 + 
+        g33L*PDstandardNth3beta3 + beta3L*PDstandardNth3g33;
+    
+    g33rhsL  =  -2*alphaL*K33L + beta1L*PDstandardNth1g13 + 
+        beta2L*PDstandardNth2g23 + 4*
+         (g13L*PDstandardNth3beta1 + g23L*PDstandardNth3beta2 + 
+           g33L*PDstandardNth3beta3) + beta3L*PDstandardNth3g33;
+    
+    K11rhsL  =  -9*PDstandardNth11alpha + G111*PDstandardNth1alpha + 
+        4*(K11L*PDstandardNth1beta1 + K12L*PDstandardNth1beta2 + 
+           K13L*PDstandardNth1beta3) + beta1L*PDstandardNth1K11 + 
+        G212*PDstandardNth2alpha + beta2L*PDstandardNth2K12 + 
+        G313*PDstandardNth3alpha + beta3L*PDstandardNth3K13 + 
         alphaL*(-2*(K11L*Km11 + K12L*Km21 + K13L*Km31) + R11 + K11L*trK);
     
-    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);
-    
-    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);
-    
-    K22rhsL  =  G122*PDstandardNth1alpha + beta1L*PDstandardNth1K22 - PDstandardNth22alpha + G222*PDstandardNth2alpha + 
-        2*(K12L*PDstandardNth2beta1 + K22L*PDstandardNth2beta2 + K23L*PDstandardNth2beta3) + beta2L*PDstandardNth2K22 + 
-        G322*PDstandardNth3alpha + beta3L*PDstandardNth3K22 + 
+    K12rhsL  =  -9*PDstandardNth12alpha + G112*PDstandardNth1alpha + 
+        3*(K12L*PDstandardNth1beta1 + K22L*PDstandardNth1beta2 + 
+           K23L*PDstandardNth1beta3) + beta1L*PDstandardNth1K12 + 
+        G222*PDstandardNth2alpha + K12L*PDstandardNth2beta1 + 
+        K22L*PDstandardNth2beta2 + K23L*PDstandardNth2beta3 + 
+        beta2L*PDstandardNth2K22 + G323*PDstandardNth3alpha + 
+        beta3L*PDstandardNth3K23 + alphaL*
+         (-2*(K11L*Km12 + K12L*Km22 + K13L*Km32) + R12 + K12L*trK);
+    
+    K13rhsL  =  -9*PDstandardNth13alpha + G113*PDstandardNth1alpha + 
+        3*(K13L*PDstandardNth1beta1 + K23L*PDstandardNth1beta2 + 
+           K33L*PDstandardNth1beta3) + beta1L*PDstandardNth1K13 + 
+        G223*PDstandardNth2alpha + beta2L*PDstandardNth2K23 + 
+        G333*PDstandardNth3alpha + K13L*PDstandardNth3beta1 + 
+        K23L*PDstandardNth3beta2 + K33L*PDstandardNth3beta3 + 
+        beta3L*PDstandardNth3K33 + alphaL*
+         (-2*(K11L*Km13 + K12L*Km23 + K13L*Km33) + R13 + K13L*trK);
+    
+    K22rhsL  =  G112*PDstandardNth1alpha + beta1L*PDstandardNth1K12 - 
+        9*PDstandardNth22alpha + G222*PDstandardNth2alpha + 
+        4*(K12L*PDstandardNth2beta1 + K22L*PDstandardNth2beta2 + 
+           K23L*PDstandardNth2beta3) + beta2L*PDstandardNth2K22 + 
+        G323*PDstandardNth3alpha + beta3L*PDstandardNth3K23 + 
         alphaL*(-2*(K12L*Km12 + K22L*Km22 + K23L*Km32) + R22 + K22L*trK);
     
-    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);
-    
-    K33rhsL  =  G133*PDstandardNth1alpha + beta1L*PDstandardNth1K33 + G233*PDstandardNth2alpha + beta2L*PDstandardNth2K33 - 
-        PDstandardNth33alpha + G333*PDstandardNth3alpha + 
-        2*(K13L*PDstandardNth3beta1 + K23L*PDstandardNth3beta2 + K33L*PDstandardNth3beta3) + beta3L*PDstandardNth3K33 + 
+    K23rhsL  =  G113*PDstandardNth1alpha + beta1L*PDstandardNth1K13 - 
+        9*PDstandardNth23alpha + G223*PDstandardNth2alpha + 
+        3*(K13L*PDstandardNth2beta1 + K23L*PDstandardNth2beta2 + 
+           K33L*PDstandardNth2beta3) + beta2L*PDstandardNth2K23 + 
+        G333*PDstandardNth3alpha + K13L*PDstandardNth3beta1 + 
+        K23L*PDstandardNth3beta2 + K33L*PDstandardNth3beta3 + 
+        beta3L*PDstandardNth3K33 + alphaL*
+         (-2*(K12L*Km13 + K22L*Km23 + K23L*Km33) + R23 + K23L*trK);
+    
+    K33rhsL  =  G113*PDstandardNth1alpha + beta1L*PDstandardNth1K13 + 
+        G223*PDstandardNth2alpha + beta2L*PDstandardNth2K23 - 
+        9*PDstandardNth33alpha + G333*PDstandardNth3alpha + 
+        4*(K13L*PDstandardNth3beta1 + K23L*PDstandardNth3beta2 + 
+           K33L*PDstandardNth3beta3) + beta3L*PDstandardNth3K33 + 
         alphaL*(-2*(K13L*Km13 + K23L*Km23 + K33L*Km33) + R33 + K33L*trK);
     
     alpharhsL  =  0;
diff --git a/ML_ADM/src/ML_ADM_boundary.c b/ML_ADM/src/ML_ADM_boundary.c
index 99b4b47..14f9375 100644
--- a/ML_ADM/src/ML_ADM_boundary.c
+++ b/ML_ADM/src/ML_ADM_boundary.c
@@ -5,7 +5,10 @@
 
 #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"
@@ -84,7 +87,7 @@ void ML_ADM_boundary_Body(cGH *cctkGH, CCTK_INT dir, CCTK_INT face, CCTK_REAL no
   pm1o12dz2 = -pow(dz,-2)/12.;
   
   /* Loop over the grid points */
-  _Pragma ("omp parallel")
+  #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])
diff --git a/ML_ADM/src/ML_ADM_constraints.c b/ML_ADM/src/ML_ADM_constraints.c
index cacb53f..6d8c808 100644
--- a/ML_ADM/src/ML_ADM_constraints.c
+++ b/ML_ADM/src/ML_ADM_constraints.c
@@ -5,7 +5,10 @@
 
 #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"
@@ -84,7 +87,7 @@ void ML_ADM_constraints_Body(cGH *cctkGH, CCTK_INT dir, CCTK_INT face, CCTK_REAL
   pm1o12dz2 = -pow(dz,-2)/12.;
   
   /* Loop over the grid points */
-  _Pragma ("omp parallel")
+  #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])
@@ -99,7 +102,8 @@ void ML_ADM_constraints_Body(cGH *cctkGH, CCTK_INT dir, CCTK_INT face, CCTK_REAL
     CCTK_REAL G111 = INITVALUE, G112 = INITVALUE, G113 = INITVALUE, G122 = INITVALUE, G123 = INITVALUE, G133 = INITVALUE;
     CCTK_REAL G211 = INITVALUE, G212 = INITVALUE, G213 = INITVALUE, G222 = INITVALUE, G223 = INITVALUE, G233 = INITVALUE;
     CCTK_REAL G311 = INITVALUE, G312 = INITVALUE, G313 = INITVALUE, G322 = INITVALUE, G323 = INITVALUE, G333 = INITVALUE;
-    CCTK_REAL gu11 = INITVALUE, gu21 = INITVALUE, gu22 = INITVALUE, gu31 = INITVALUE, gu32 = INITVALUE, gu33 = INITVALUE;
+    CCTK_REAL gu11 = INITVALUE, gu12 = INITVALUE, gu13 = INITVALUE, gu21 = INITVALUE, gu22 = INITVALUE, gu23 = INITVALUE;
+    CCTK_REAL gu31 = INITVALUE, gu32 = INITVALUE, gu33 = INITVALUE;
     CCTK_REAL Km11 = INITVALUE, Km12 = INITVALUE, Km13 = INITVALUE, Km21 = INITVALUE, Km22 = INITVALUE, Km23 = INITVALUE;
     CCTK_REAL Km31 = INITVALUE, Km32 = INITVALUE, Km33 = INITVALUE;
     CCTK_REAL R11 = INITVALUE, R12 = INITVALUE, R13 = INITVALUE, R22 = INITVALUE, R23 = INITVALUE, R33 = INITVALUE;
@@ -117,31 +121,30 @@ void ML_ADM_constraints_Body(cGH *cctkGH, CCTK_INT dir, CCTK_INT face, CCTK_REAL
     CCTK_REAL PDstandardNth1g11 = INITVALUE;
     CCTK_REAL PDstandardNth2g11 = INITVALUE;
     CCTK_REAL PDstandardNth3g11 = INITVALUE;
+    CCTK_REAL PDstandardNth11g11 = INITVALUE;
     CCTK_REAL PDstandardNth22g11 = INITVALUE;
     CCTK_REAL PDstandardNth33g11 = INITVALUE;
     CCTK_REAL PDstandardNth12g11 = INITVALUE;
     CCTK_REAL PDstandardNth13g11 = INITVALUE;
-    CCTK_REAL PDstandardNth21g11 = INITVALUE;
     CCTK_REAL PDstandardNth23g11 = INITVALUE;
-    CCTK_REAL PDstandardNth31g11 = INITVALUE;
-    CCTK_REAL PDstandardNth32g11 = INITVALUE;
     CCTK_REAL PDstandardNth1g12 = INITVALUE;
     CCTK_REAL PDstandardNth2g12 = INITVALUE;
     CCTK_REAL PDstandardNth3g12 = INITVALUE;
+    CCTK_REAL PDstandardNth11g12 = INITVALUE;
+    CCTK_REAL PDstandardNth22g12 = INITVALUE;
     CCTK_REAL PDstandardNth33g12 = INITVALUE;
     CCTK_REAL PDstandardNth12g12 = INITVALUE;
     CCTK_REAL PDstandardNth13g12 = INITVALUE;
     CCTK_REAL PDstandardNth21g12 = INITVALUE;
     CCTK_REAL PDstandardNth23g12 = INITVALUE;
-    CCTK_REAL PDstandardNth31g12 = INITVALUE;
-    CCTK_REAL PDstandardNth32g12 = INITVALUE;
     CCTK_REAL PDstandardNth1g13 = INITVALUE;
     CCTK_REAL PDstandardNth2g13 = INITVALUE;
     CCTK_REAL PDstandardNth3g13 = INITVALUE;
+    CCTK_REAL PDstandardNth11g13 = INITVALUE;
     CCTK_REAL PDstandardNth22g13 = INITVALUE;
+    CCTK_REAL PDstandardNth33g13 = INITVALUE;
     CCTK_REAL PDstandardNth12g13 = INITVALUE;
     CCTK_REAL PDstandardNth13g13 = INITVALUE;
-    CCTK_REAL PDstandardNth21g13 = INITVALUE;
     CCTK_REAL PDstandardNth23g13 = INITVALUE;
     CCTK_REAL PDstandardNth31g13 = INITVALUE;
     CCTK_REAL PDstandardNth32g13 = INITVALUE;
@@ -149,17 +152,18 @@ void ML_ADM_constraints_Body(cGH *cctkGH, CCTK_INT dir, CCTK_INT face, CCTK_REAL
     CCTK_REAL PDstandardNth2g22 = INITVALUE;
     CCTK_REAL PDstandardNth3g22 = INITVALUE;
     CCTK_REAL PDstandardNth11g22 = INITVALUE;
+    CCTK_REAL PDstandardNth22g22 = INITVALUE;
     CCTK_REAL PDstandardNth33g22 = INITVALUE;
     CCTK_REAL PDstandardNth12g22 = INITVALUE;
     CCTK_REAL PDstandardNth13g22 = INITVALUE;
     CCTK_REAL PDstandardNth21g22 = INITVALUE;
     CCTK_REAL PDstandardNth23g22 = INITVALUE;
-    CCTK_REAL PDstandardNth31g22 = INITVALUE;
-    CCTK_REAL PDstandardNth32g22 = INITVALUE;
     CCTK_REAL PDstandardNth1g23 = INITVALUE;
     CCTK_REAL PDstandardNth2g23 = INITVALUE;
     CCTK_REAL PDstandardNth3g23 = INITVALUE;
     CCTK_REAL PDstandardNth11g23 = INITVALUE;
+    CCTK_REAL PDstandardNth22g23 = INITVALUE;
+    CCTK_REAL PDstandardNth33g23 = INITVALUE;
     CCTK_REAL PDstandardNth12g23 = INITVALUE;
     CCTK_REAL PDstandardNth13g23 = INITVALUE;
     CCTK_REAL PDstandardNth21g23 = INITVALUE;
@@ -171,12 +175,13 @@ void ML_ADM_constraints_Body(cGH *cctkGH, CCTK_INT dir, CCTK_INT face, CCTK_REAL
     CCTK_REAL PDstandardNth3g33 = INITVALUE;
     CCTK_REAL PDstandardNth11g33 = INITVALUE;
     CCTK_REAL PDstandardNth22g33 = INITVALUE;
+    CCTK_REAL PDstandardNth33g33 = INITVALUE;
     CCTK_REAL PDstandardNth12g33 = INITVALUE;
     CCTK_REAL PDstandardNth13g33 = INITVALUE;
-    CCTK_REAL PDstandardNth21g33 = INITVALUE;
     CCTK_REAL PDstandardNth23g33 = INITVALUE;
     CCTK_REAL PDstandardNth31g33 = INITVALUE;
     CCTK_REAL PDstandardNth32g33 = INITVALUE;
+    CCTK_REAL PDstandardNth1K11 = INITVALUE;
     CCTK_REAL PDstandardNth2K11 = INITVALUE;
     CCTK_REAL PDstandardNth3K11 = INITVALUE;
     CCTK_REAL PDstandardNth1K12 = INITVALUE;
@@ -186,12 +191,14 @@ void ML_ADM_constraints_Body(cGH *cctkGH, CCTK_INT dir, CCTK_INT face, CCTK_REAL
     CCTK_REAL PDstandardNth2K13 = INITVALUE;
     CCTK_REAL PDstandardNth3K13 = INITVALUE;
     CCTK_REAL PDstandardNth1K22 = INITVALUE;
+    CCTK_REAL PDstandardNth2K22 = INITVALUE;
     CCTK_REAL PDstandardNth3K22 = INITVALUE;
     CCTK_REAL PDstandardNth1K23 = INITVALUE;
     CCTK_REAL PDstandardNth2K23 = INITVALUE;
     CCTK_REAL PDstandardNth3K23 = INITVALUE;
     CCTK_REAL PDstandardNth1K33 = INITVALUE;
     CCTK_REAL PDstandardNth2K33 = INITVALUE;
+    CCTK_REAL PDstandardNth3K33 = INITVALUE;
     
     /* Assign local copies of grid functions */
     g11L = g11[index];
@@ -215,12 +222,17 @@ void ML_ADM_constraints_Body(cGH *cctkGH, CCTK_INT dir, CCTK_INT face, CCTK_REAL
     PDstandardNth1g11 = PDstandardNth1(g11, i, j, k);
     PDstandardNth2g11 = PDstandardNth2(g11, i, j, k);
     PDstandardNth3g11 = PDstandardNth3(g11, i, j, k);
+    PDstandardNth11g11 = PDstandardNth11(g11, i, j, k);
     PDstandardNth22g11 = PDstandardNth22(g11, i, j, k);
     PDstandardNth33g11 = PDstandardNth33(g11, i, j, k);
+    PDstandardNth12g11 = PDstandardNth12(g11, i, j, k);
+    PDstandardNth13g11 = PDstandardNth13(g11, i, j, k);
     PDstandardNth23g11 = PDstandardNth23(g11, i, j, k);
     PDstandardNth1g12 = PDstandardNth1(g12, i, j, k);
     PDstandardNth2g12 = PDstandardNth2(g12, i, j, k);
     PDstandardNth3g12 = PDstandardNth3(g12, i, j, k);
+    PDstandardNth11g12 = PDstandardNth11(g12, i, j, k);
+    PDstandardNth22g12 = PDstandardNth22(g12, i, j, k);
     PDstandardNth33g12 = PDstandardNth33(g12, i, j, k);
     PDstandardNth12g12 = PDstandardNth12(g12, i, j, k);
     PDstandardNth13g12 = PDstandardNth13(g12, i, j, k);
@@ -228,7 +240,9 @@ void ML_ADM_constraints_Body(cGH *cctkGH, CCTK_INT dir, CCTK_INT face, CCTK_REAL
     PDstandardNth1g13 = PDstandardNth1(g13, i, j, k);
     PDstandardNth2g13 = PDstandardNth2(g13, i, j, k);
     PDstandardNth3g13 = PDstandardNth3(g13, i, j, k);
+    PDstandardNth11g13 = PDstandardNth11(g13, i, j, k);
     PDstandardNth22g13 = PDstandardNth22(g13, i, j, k);
+    PDstandardNth33g13 = PDstandardNth33(g13, i, j, k);
     PDstandardNth12g13 = PDstandardNth12(g13, i, j, k);
     PDstandardNth13g13 = PDstandardNth13(g13, i, j, k);
     PDstandardNth23g13 = PDstandardNth23(g13, i, j, k);
@@ -236,12 +250,17 @@ void ML_ADM_constraints_Body(cGH *cctkGH, CCTK_INT dir, CCTK_INT face, CCTK_REAL
     PDstandardNth2g22 = PDstandardNth2(g22, i, j, k);
     PDstandardNth3g22 = PDstandardNth3(g22, i, j, k);
     PDstandardNth11g22 = PDstandardNth11(g22, i, j, k);
+    PDstandardNth22g22 = PDstandardNth22(g22, i, j, k);
     PDstandardNth33g22 = PDstandardNth33(g22, i, j, k);
+    PDstandardNth12g22 = PDstandardNth12(g22, i, j, k);
     PDstandardNth13g22 = PDstandardNth13(g22, i, j, k);
+    PDstandardNth23g22 = PDstandardNth23(g22, i, j, k);
     PDstandardNth1g23 = PDstandardNth1(g23, i, j, k);
     PDstandardNth2g23 = PDstandardNth2(g23, i, j, k);
     PDstandardNth3g23 = PDstandardNth3(g23, i, j, k);
     PDstandardNth11g23 = PDstandardNth11(g23, i, j, k);
+    PDstandardNth22g23 = PDstandardNth22(g23, i, j, k);
+    PDstandardNth33g23 = PDstandardNth33(g23, i, j, k);
     PDstandardNth12g23 = PDstandardNth12(g23, i, j, k);
     PDstandardNth13g23 = PDstandardNth13(g23, i, j, k);
     PDstandardNth23g23 = PDstandardNth23(g23, i, j, k);
@@ -250,7 +269,11 @@ void ML_ADM_constraints_Body(cGH *cctkGH, CCTK_INT dir, CCTK_INT face, CCTK_REAL
     PDstandardNth3g33 = PDstandardNth3(g33, i, j, k);
     PDstandardNth11g33 = PDstandardNth11(g33, i, j, k);
     PDstandardNth22g33 = PDstandardNth22(g33, i, j, k);
+    PDstandardNth33g33 = PDstandardNth33(g33, i, j, k);
     PDstandardNth12g33 = PDstandardNth12(g33, i, j, k);
+    PDstandardNth13g33 = PDstandardNth13(g33, i, j, k);
+    PDstandardNth23g33 = PDstandardNth23(g33, i, j, k);
+    PDstandardNth1K11 = PDstandardNth1(K11, i, j, k);
     PDstandardNth2K11 = PDstandardNth2(K11, i, j, k);
     PDstandardNth3K11 = PDstandardNth3(K11, i, j, k);
     PDstandardNth1K12 = PDstandardNth1(K12, i, j, k);
@@ -260,170 +283,263 @@ void ML_ADM_constraints_Body(cGH *cctkGH, CCTK_INT dir, CCTK_INT face, CCTK_REAL
     PDstandardNth2K13 = PDstandardNth2(K13, i, j, k);
     PDstandardNth3K13 = PDstandardNth3(K13, i, j, k);
     PDstandardNth1K22 = PDstandardNth1(K22, i, j, k);
+    PDstandardNth2K22 = PDstandardNth2(K22, i, j, k);
     PDstandardNth3K22 = PDstandardNth3(K22, i, j, k);
     PDstandardNth1K23 = PDstandardNth1(K23, i, j, k);
     PDstandardNth2K23 = PDstandardNth2(K23, i, j, k);
     PDstandardNth3K23 = PDstandardNth3(K23, i, j, k);
     PDstandardNth1K33 = PDstandardNth1(K33, i, j, k);
     PDstandardNth2K33 = PDstandardNth2(K33, i, j, k);
+    PDstandardNth3K33 = PDstandardNth3(K33, i, j, k);
     
     /* Precompute derivatives (old style) */
     
     /* Calculate temporaries and grid functions */
-    detg  =  2*g12L*g13L*g23L + g33L*(g11L*g22L - SQR(g12L)) - g22L*SQR(g13L) - g11L*SQR(g23L);
+    detg  =  2*g12L*g13L*g23L + g33L*(g11L*g22L - SQR(g12L)) - g22L*SQR(g13L) - 
+        g11L*SQR(g23L);
     
     gu11  =  INV(detg)*(g22L*g33L - SQR(g23L));
     
-    gu21  =  (g13L*g23L - g12L*g33L)*INV(detg);
+    gu12  =  (g13L*g23L - g12L*g33L)*INV(detg);
     
-    gu31  =  (-(g13L*g22L) + g12L*g23L)*INV(detg);
+    gu13  =  (-(g13L*g22L) + g12L*g23L)*INV(detg);
+    
+    gu21  =  (g13L*g23L - g12L*g33L)*INV(detg);
     
     gu22  =  INV(detg)*(g11L*g33L - SQR(g13L));
     
+    gu23  =  (g12L*g13L - g11L*g23L)*INV(detg);
+    
+    gu31  =  (-(g13L*g22L) + g12L*g23L)*INV(detg);
+    
     gu32  =  (g12L*g13L - g11L*g23L)*INV(detg);
     
     gu33  =  INV(detg)*(g11L*g22L - SQR(g12L));
     
-    G111  =  khalf*(gu11*PDstandardNth1g11 + 2*(gu21*PDstandardNth1g12 + gu31*PDstandardNth1g13) - gu21*PDstandardNth2g11 - 
-          gu31*PDstandardNth3g11);
+    G111  =  khalf*(gu11*PDstandardNth1g11 + 
+          2*(gu12*PDstandardNth1g12 + gu13*PDstandardNth1g13) - 
+          gu12*PDstandardNth2g11 - gu13*PDstandardNth3g11);
     
-    G211  =  khalf*(gu21*PDstandardNth1g11 + 2*(gu22*PDstandardNth1g12 + gu32*PDstandardNth1g13) - gu22*PDstandardNth2g11 - 
-          gu32*PDstandardNth3g11);
+    G211  =  khalf*(gu21*PDstandardNth1g11 + 
+          2*(gu22*PDstandardNth1g12 + gu23*PDstandardNth1g13) - 
+          gu22*PDstandardNth2g11 - gu23*PDstandardNth3g11);
     
-    G311  =  khalf*(gu31*PDstandardNth1g11 + 2*(gu32*PDstandardNth1g12 + gu33*PDstandardNth1g13) - gu32*PDstandardNth2g11 - 
-          gu33*PDstandardNth3g11);
+    G311  =  khalf*(gu31*PDstandardNth1g11 + 
+          2*(gu32*PDstandardNth1g12 + gu33*PDstandardNth1g13) - 
+          gu32*PDstandardNth2g11 - gu33*PDstandardNth3g11);
     
-    G112  =  khalf*(gu21*PDstandardNth1g22 + gu11*PDstandardNth2g11 + 
-          gu31*(PDstandardNth1g23 + PDstandardNth2g13 - PDstandardNth3g12));
+    G112  =  khalf*(gu12*PDstandardNth1g22 + gu11*PDstandardNth2g11 + 
+          gu13*(PDstandardNth1g23 + PDstandardNth2g13 - PDstandardNth3g12));
     
     G212  =  khalf*(gu22*PDstandardNth1g22 + gu21*PDstandardNth2g11 + 
-          gu32*(PDstandardNth1g23 + PDstandardNth2g13 - PDstandardNth3g12));
+          gu23*(PDstandardNth1g23 + PDstandardNth2g13 - PDstandardNth3g12));
     
     G312  =  khalf*(gu32*PDstandardNth1g22 + gu31*PDstandardNth2g11 + 
           gu33*(PDstandardNth1g23 + PDstandardNth2g13 - PDstandardNth3g12));
     
-    G113  =  khalf*(gu31*PDstandardNth1g33 + gu11*PDstandardNth3g11 + 
-          gu21*(PDstandardNth1g23 - PDstandardNth2g13 + PDstandardNth3g12));
+    G113  =  khalf*(gu13*PDstandardNth1g33 + gu11*PDstandardNth3g11 + 
+          gu12*(PDstandardNth1g23 - PDstandardNth2g13 + PDstandardNth3g12));
     
-    G213  =  khalf*(gu32*PDstandardNth1g33 + gu21*PDstandardNth3g11 + 
+    G213  =  khalf*(gu23*PDstandardNth1g33 + gu21*PDstandardNth3g11 + 
           gu22*(PDstandardNth1g23 - PDstandardNth2g13 + PDstandardNth3g12));
     
     G313  =  khalf*(gu33*PDstandardNth1g33 + gu31*PDstandardNth3g11 + 
           gu32*(PDstandardNth1g23 - PDstandardNth2g13 + PDstandardNth3g12));
     
-    G122  =  khalf*(gu11*(-PDstandardNth1g22 + 2*PDstandardNth2g12) + gu21*PDstandardNth2g22 + 
-          gu31*(2*PDstandardNth2g23 - PDstandardNth3g22));
+    G122  =  khalf*(gu11*(-PDstandardNth1g22 + 2*PDstandardNth2g12) + 
+          gu12*PDstandardNth2g22 + gu13*(2*PDstandardNth2g23 - PDstandardNth3g22));
     
-    G222  =  khalf*(gu21*(-PDstandardNth1g22 + 2*PDstandardNth2g12) + gu22*PDstandardNth2g22 + 
-          gu32*(2*PDstandardNth2g23 - PDstandardNth3g22));
+    G222  =  khalf*(gu21*(-PDstandardNth1g22 + 2*PDstandardNth2g12) + 
+          gu22*PDstandardNth2g22 + gu23*(2*PDstandardNth2g23 - PDstandardNth3g22));
     
-    G322  =  khalf*(gu31*(-PDstandardNth1g22 + 2*PDstandardNth2g12) + gu32*PDstandardNth2g22 + 
-          gu33*(2*PDstandardNth2g23 - PDstandardNth3g22));
+    G322  =  khalf*(gu31*(-PDstandardNth1g22 + 2*PDstandardNth2g12) + 
+          gu32*PDstandardNth2g22 + gu33*(2*PDstandardNth2g23 - PDstandardNth3g22));
     
-    G123  =  khalf*(gu31*PDstandardNth2g33 + gu11*(-PDstandardNth1g23 + PDstandardNth2g13 + PDstandardNth3g12) + 
-          gu21*PDstandardNth3g22);
+    G123  =  khalf*(gu13*PDstandardNth2g33 + 
+          gu11*(-PDstandardNth1g23 + PDstandardNth2g13 + PDstandardNth3g12) + 
+          gu12*PDstandardNth3g22);
     
-    G223  =  khalf*(gu32*PDstandardNth2g33 + gu21*(-PDstandardNth1g23 + PDstandardNth2g13 + PDstandardNth3g12) + 
+    G223  =  khalf*(gu23*PDstandardNth2g33 + 
+          gu21*(-PDstandardNth1g23 + PDstandardNth2g13 + PDstandardNth3g12) + 
           gu22*PDstandardNth3g22);
     
-    G323  =  khalf*(gu33*PDstandardNth2g33 + gu31*(-PDstandardNth1g23 + PDstandardNth2g13 + PDstandardNth3g12) + 
+    G323  =  khalf*(gu33*PDstandardNth2g33 + 
+          gu31*(-PDstandardNth1g23 + PDstandardNth2g13 + PDstandardNth3g12) + 
           gu32*PDstandardNth3g22);
     
-    G133  =  khalf*(-(gu11*PDstandardNth1g33) - gu21*PDstandardNth2g33 + 2*gu11*PDstandardNth3g13 + 
-          2*gu21*PDstandardNth3g23 + gu31*PDstandardNth3g33);
-    
-    G233  =  khalf*(-(gu21*PDstandardNth1g33) - gu22*PDstandardNth2g33 + 2*gu21*PDstandardNth3g13 + 
-          2*gu22*PDstandardNth3g23 + gu32*PDstandardNth3g33);
-    
-    G333  =  khalf*(-(gu31*PDstandardNth1g33) - gu32*PDstandardNth2g33 + 2*gu31*PDstandardNth3g13 + 
-          2*gu32*PDstandardNth3g23 + gu33*PDstandardNth3g33);
-    
-    R11  =  khalf*(-(gu22*PDstandardNth11g22) - 2*(G111*G122 + G111*G133 + G211*G222 + G211*G233 + G311*G322 + G311*G333 + 
-             gu32*PDstandardNth11g23) - gu33*PDstandardNth11g33 + 2*gu22*PDstandardNth12g12 + 2*gu32*PDstandardNth12g13 + 
-          2*gu32*PDstandardNth13g12 + 2*gu33*PDstandardNth13g13 - gu22*PDstandardNth22g11 - 2*gu32*PDstandardNth23g11 - 
-          gu33*PDstandardNth33g11 + 2*SQR(G112) + 2*SQR(G113) + 2*SQR(G212) + 2*SQR(G213) + 2*SQR(G312) + 2*SQR(G313));
-    
-    R12  =  khalf*(2*(G113*G123 + G213*G223 + G313*G323) - 2*(G112*G133 + G212*G233 + G312*G333 + gu21*PDstandardNth12g12) - 
-          gu32*PDstandardNth12g23 - gu33*PDstandardNth12g33 + 
-          gu31*(PDstandardNth11g23 - PDstandardNth12g13 - PDstandardNth13g12) + gu32*PDstandardNth13g22 + 
-          gu33*PDstandardNth13g23 + gu21*(PDstandardNth11g22 + PDstandardNth22g11) + gu32*PDstandardNth22g13 + 
-          gu31*PDstandardNth23g11 - gu32*PDstandardNth23g12 + gu33*PDstandardNth23g13 - gu33*PDstandardNth33g12);
-    
-    R13  =  khalf*(2*(G112*G123 + G212*G223 + G312*G323) - 2*(G113*G122 + G213*G222 + G313*G322 + gu31*PDstandardNth13g13) + 
-          gu21*(PDstandardNth11g23 - PDstandardNth12g13 - PDstandardNth13g12 + PDstandardNth23g11) + 
-          gu22*(PDstandardNth12g23 - PDstandardNth13g22 - PDstandardNth22g13 + PDstandardNth23g12) + 
-          gu31*(PDstandardNth11g33 + PDstandardNth33g11) + 
-          gu32*(PDstandardNth12g33 - PDstandardNth13g23 - PDstandardNth23g13 + PDstandardNth33g12));
-    
-    R22  =  khalf*(-2*(G122*(G111 + G133) + G222*(G211 + G233) + G322*(G311 + G333) + gu31*PDstandardNth13g22) + 
-          gu11*(-PDstandardNth11g22 + 2*PDstandardNth12g12 - PDstandardNth22g11) + 
-          gu31*(-2*PDstandardNth22g13 + 2*(PDstandardNth12g23 + PDstandardNth23g12)) + 
-          gu33*(-PDstandardNth22g33 + 2*PDstandardNth23g23 - PDstandardNth33g22) + 
-          2*(SQR(G112) + SQR(G123) + SQR(G212) + SQR(G223) + SQR(G312) + SQR(G323)));
-    
-    R23  =  khalf*(2*(G112*G113 + G212*G213 + G312*G313) + 
-          gu11*(-PDstandardNth11g23 + PDstandardNth12g13 + PDstandardNth13g12 - PDstandardNth23g11) + 
-          gu21*(-PDstandardNth12g23 + PDstandardNth13g22 + PDstandardNth22g13 - PDstandardNth23g12) - 
-          2*(G111*G123 + G211*G223 + G311*G323 + gu32*PDstandardNth23g23) + 
-          gu31*(PDstandardNth12g33 - PDstandardNth13g23 - PDstandardNth23g13 + PDstandardNth33g12) + 
-          gu32*(PDstandardNth22g33 + PDstandardNth33g22));
-    
-    R33  =  khalf*(gu11*(-PDstandardNth11g33 + 2*PDstandardNth13g13 - PDstandardNth33g11) - 
-          2*((G111 + G122)*G133 + (G211 + G222)*G233 + (G311 + G322)*G333 + 
-             gu21*(PDstandardNth12g33 + PDstandardNth33g12)) + 
-          gu22*(-PDstandardNth22g33 + 2*PDstandardNth23g23 - PDstandardNth33g22) + 
-          2*(gu21*(PDstandardNth13g23 + PDstandardNth23g13) + SQR(G113) + SQR(G123) + SQR(G213) + SQR(G223) + SQR(G313) + 
-             SQR(G323)));
-    
-    trR  =  gu11*R11 + gu22*R22 + 2*(gu21*R12 + gu31*R13 + gu32*R23) + gu33*R33;
-    
-    Km11  =  gu11*K11L + gu21*K12L + gu31*K13L;
-    
-    Km21  =  gu21*K11L + gu22*K12L + gu32*K13L;
+    G133  =  khalf*(-(gu11*PDstandardNth1g33) - gu12*PDstandardNth2g33 + 
+          2*gu11*PDstandardNth3g13 + 2*gu12*PDstandardNth3g23 + 
+          gu13*PDstandardNth3g33);
+    
+    G233  =  khalf*(-(gu21*PDstandardNth1g33) - gu22*PDstandardNth2g33 + 
+          2*gu21*PDstandardNth3g13 + 2*gu22*PDstandardNth3g23 + 
+          gu23*PDstandardNth3g33);
+    
+    G333  =  khalf*(-(gu31*PDstandardNth1g33) - gu32*PDstandardNth2g33 + 
+          2*gu31*PDstandardNth3g13 + 2*gu32*PDstandardNth3g23 + 
+          gu33*PDstandardNth3g33);
+    
+    R11  =  -(G111*(G111 + G122 + G133)) - G211*(G211 + G222 + G233) - 
+        G311*(G311 + G322 + G333) + khalf*
+         (-6*gu11*PDstandardNth11g11 + 
+           gu12*(-8*PDstandardNth11g12 + 2*PDstandardNth12g12) + 
+           gu22*(-9*PDstandardNth11g22 + 3*PDstandardNth12g22) + 
+           gu13*(-8*PDstandardNth11g13 + 2*PDstandardNth13g13) + 
+           gu33*(-9*PDstandardNth11g33 + 3*PDstandardNth13g33) + 
+           gu21*(-6*PDstandardNth11g12 - PDstandardNth12g12 + PDstandardNth22g12) + 
+           gu23*(-9*PDstandardNth11g23 + 3*PDstandardNth13g23 + 
+              PDstandardNth22g23 - PDstandardNth23g23) + 
+           gu31*(-6*PDstandardNth11g13 - PDstandardNth13g13 + PDstandardNth33g13) + 
+           gu32*(-9*PDstandardNth11g23 + 3*PDstandardNth12g23 - 
+              PDstandardNth23g23 + PDstandardNth33g23)) + SQR(G111) + SQR(G112) + 
+        SQR(G113) + SQR(G211) + SQR(G212) + SQR(G213) + SQR(G311) + SQR(G312) + 
+        SQR(G313);
+    
+    R12  =  khalf*(-2*(G112*G133 + G212*G233 + G312*G333) - 
+          6*(gu11*PDstandardNth12g11 + gu31*PDstandardNth12g13) - 
+          9*(gu22*PDstandardNth12g22 + gu23*PDstandardNth12g23 + 
+             gu32*PDstandardNth12g23 + gu33*PDstandardNth12g33) - 
+          gu31*PDstandardNth13g13 + gu21*
+           (-7*PDstandardNth12g12 + PDstandardNth22g12) + 
+          gu12*(PDstandardNth11g12 - 10*PDstandardNth12g12 + 
+             3*PDstandardNth22g12) + 3*gu22*PDstandardNth22g22 + 
+          gu23*PDstandardNth22g23 + 3*gu32*PDstandardNth22g23 + 
+          gu13*(PDstandardNth11g13 - 9*PDstandardNth12g13 - PDstandardNth13g13 + 
+             3*PDstandardNth23g13) - gu32*PDstandardNth23g23 + 
+          2*(G113*G123 + G213*G223 + G313*G323 + gu23*PDstandardNth23g23) + 
+          3*gu33*PDstandardNth23g33 + gu31*PDstandardNth33g13 + 
+          gu32*PDstandardNth33g23);
+    
+    R13  =  khalf*(-2*(G113*G122 + G213*G222 + G313*G322) - 
+          6*gu11*PDstandardNth13g11 - 
+          9*(gu22*PDstandardNth13g22 + (gu23 + gu32)*PDstandardNth13g23 + 
+             gu33*PDstandardNth13g33) + 
+          gu21*(-PDstandardNth12g12 - 6*PDstandardNth13g12 + PDstandardNth22g12) + 
+          gu12*(PDstandardNth11g12 - PDstandardNth12g12 - 9*PDstandardNth13g12 + 
+             3*PDstandardNth23g12) + 
+          gu23*(PDstandardNth22g23 - PDstandardNth23g23) + 
+          2*(G112*G123 + G212*G223 + G312*G323 + gu32*PDstandardNth23g23) + 
+          gu31*(-7*PDstandardNth13g13 + PDstandardNth33g13) + 
+          gu13*(PDstandardNth11g13 - 10*PDstandardNth13g13 + 
+             3*PDstandardNth33g13) + gu32*PDstandardNth33g23 + 
+          3*(gu22*PDstandardNth23g22 + gu23*PDstandardNth33g23 + 
+             gu33*PDstandardNth33g33));
+    
+    R22  =  -(G122*(G111 + G122 + G133)) - G222*(G211 + G222 + G233) - 
+        G322*(G311 + G322 + G333) + khalf*
+         (3*gu11*(PDstandardNth12g11 - 3*PDstandardNth22g11) + 
+           gu21*(2*PDstandardNth12g12 - 8*PDstandardNth22g12) + 
+           gu12*(PDstandardNth11g12 - PDstandardNth12g12 - 6*PDstandardNth22g12) - 
+           6*gu22*PDstandardNth22g22 + 
+           gu13*(PDstandardNth11g13 - PDstandardNth13g13 - 9*PDstandardNth22g13 + 
+              3*PDstandardNth23g13) - 
+           gu23*(8*PDstandardNth22g23 - 2*PDstandardNth23g23) + 
+           gu33*(-9*PDstandardNth22g33 + 3*PDstandardNth23g33) + 
+           gu31*(3*PDstandardNth12g13 - PDstandardNth13g13 - 9*PDstandardNth22g13 + 
+              PDstandardNth33g13) + gu32*
+            (-6*PDstandardNth22g23 - PDstandardNth23g23 + PDstandardNth33g23)) + 
+        SQR(G112) + SQR(G122) + SQR(G123) + SQR(G212) + SQR(G222) + SQR(G223) + 
+        SQR(G312) + SQR(G322) + SQR(G323);
+    
+    R23  =  khalf*(-2*(G111*G123 + G211*G223 + G311*G323) + 
+          2*(G112*G113 + G212*G213 + G312*G313 + gu31*PDstandardNth13g13) + 
+          gu11*(3*PDstandardNth13g11 - 9*PDstandardNth23g11) + 
+          gu21*(-PDstandardNth12g12 + 3*PDstandardNth13g12 + PDstandardNth22g12 - 
+             9*PDstandardNth23g12) + 
+          gu12*(PDstandardNth11g12 - PDstandardNth12g12 - 6*PDstandardNth23g12) - 
+          6*gu22*PDstandardNth23g22 - 
+          9*(gu31*PDstandardNth23g13 + gu33*PDstandardNth23g33) + 
+          gu31*PDstandardNth33g13 + gu13*
+           (PDstandardNth11g13 - PDstandardNth13g13 - 9*PDstandardNth23g13 + 
+             3*PDstandardNth33g13) + 
+          gu32*(-7*PDstandardNth23g23 + PDstandardNth33g23) + 
+          gu23*(PDstandardNth22g23 - 10*PDstandardNth23g23 + 
+             3*PDstandardNth33g23) + 3*gu33*PDstandardNth33g33);
+    
+    R33  =  -(G133*(G111 + G122 + G133)) - G233*(G211 + G222 + G233) - 
+        G333*(G311 + G322 + G333) + khalf*
+         (3*gu11*(PDstandardNth13g11 - 3*PDstandardNth33g11) + 
+           gu21*(-PDstandardNth12g12 + 3*PDstandardNth13g12 + PDstandardNth22g12 - 
+              9*PDstandardNth33g12) + 
+           gu12*(PDstandardNth11g12 - PDstandardNth12g12 + 3*PDstandardNth23g12 - 
+              9*PDstandardNth33g12) + 
+           gu31*(2*PDstandardNth13g13 - 8*PDstandardNth33g13) + 
+           gu13*(PDstandardNth11g13 - PDstandardNth13g13 - 6*PDstandardNth33g13) + 
+           3*gu22*(PDstandardNth23g22 - 3*PDstandardNth33g22) + 
+           gu32*(2*PDstandardNth23g23 - 8*PDstandardNth33g23) + 
+           gu23*(PDstandardNth22g23 - PDstandardNth23g23 - 6*PDstandardNth33g23) - 
+           6*gu33*PDstandardNth33g33) + SQR(G113) + SQR(G123) + SQR(G133) + 
+        SQR(G213) + SQR(G223) + SQR(G233) + SQR(G313) + SQR(G323) + SQR(G333);
+    
+    trR  =  gu11*R11 + (gu12 + gu21)*R12 + (gu13 + gu31)*R13 + gu22*R22 + 
+        (gu23 + gu32)*R23 + gu33*R33;
+    
+    Km11  =  gu11*K11L + gu12*K12L + gu13*K13L;
+    
+    Km21  =  gu21*K11L + gu22*K12L + gu23*K13L;
     
     Km31  =  gu31*K11L + gu32*K12L + gu33*K13L;
     
-    Km12  =  gu11*K12L + gu21*K22L + gu31*K23L;
+    Km12  =  gu11*K12L + gu12*K22L + gu13*K23L;
     
-    Km22  =  gu21*K12L + gu22*K22L + gu32*K23L;
+    Km22  =  gu21*K12L + gu22*K22L + gu23*K23L;
     
     Km32  =  gu31*K12L + gu32*K22L + gu33*K23L;
     
-    Km13  =  gu11*K13L + gu21*K23L + gu31*K33L;
+    Km13  =  gu11*K13L + gu12*K23L + gu13*K33L;
     
-    Km23  =  gu21*K13L + gu22*K23L + gu32*K33L;
+    Km23  =  gu21*K13L + gu22*K23L + gu23*K33L;
     
     Km33  =  gu31*K13L + gu32*K23L + gu33*K33L;
     
     trK  =  Km11 + Km22 + Km33;
     
-    HL  =  -2*(Km12*Km21 + Km13*Km31 + Km23*Km32) + trR - SQR(Km11) - SQR(Km22) - SQR(Km33) + SQR(trK);
-    
-    M1L  =  gu21*(-(G112*K11L) + G111*K12L - G212*K12L - G312*K13L + G211*K22L + G311*K23L - PDstandardNth1K12 + 
-           PDstandardNth2K11) + gu22*(-(G122*K11L) + G112*K12L - G222*K12L - G322*K13L + G212*K22L + G312*K23L - 
-           PDstandardNth1K22 + PDstandardNth2K12) + gu31*
-         (-(G113*K11L) - G213*K12L + G111*K13L - G313*K13L + G211*K23L + G311*K33L - PDstandardNth1K13 + PDstandardNth3K11)\
-         + gu32*(G113*K12L + G112*K13L + G213*K22L + (G212 + G313)*K23L + G312*K33L - 
-           2*(G123*K11L + G223*K12L + G323*K13L + PDstandardNth1K23) + PDstandardNth2K13 + PDstandardNth3K12) + 
-        gu33*(-(G133*K11L) - G233*K12L + G113*K13L - G333*K13L + G213*K23L + G313*K33L - PDstandardNth1K33 + 
-           PDstandardNth3K13);
-    
-    M2L  =  gu11*(G112*K11L + (-G111 + G212)*K12L + G312*K13L - G211*K22L - G311*K23L + PDstandardNth1K12 - 
-           PDstandardNth2K11) + gu21*(G122*K11L + (-G112 + G222)*K12L + G322*K13L - G212*K22L - G312*K23L + 
-           PDstandardNth1K22 - PDstandardNth2K12) + gu31*
-         (G123*K11L + (-2*G113 + G223)*K12L + (G112 + G323)*K13L + G212*K23L + G312*K33L + PDstandardNth1K23 - 
-           2*(G213*K22L + G313*K23L + PDstandardNth2K13) + PDstandardNth3K12) + 
-        gu32*(-(G123*K12L) + G122*K13L - G223*K22L + G222*K23L - G323*K23L + G322*K33L - PDstandardNth2K23 + 
-           PDstandardNth3K22) + gu33*(-(G133*K12L) + G123*K13L - G233*K22L + G223*K23L - G333*K23L + G323*K33L - 
-           PDstandardNth2K33 + PDstandardNth3K23);
-    
-    M3L  =  gu11*(G113*K11L + G213*K12L + (-G111 + G313)*K13L - G211*K23L - G311*K33L + PDstandardNth1K13 - 
-           PDstandardNth3K11) + gu21*(G123*K11L + (G113 + G223)*K12L + (-2*G112 + G323)*K13L + G213*K22L + 
-           (-2*G212 + G313)*K23L + PDstandardNth1K23 + PDstandardNth2K13 - 2*(G312*K33L + PDstandardNth3K12)) + 
-        gu31*(G133*K11L + G233*K12L + (-G113 + G333)*K13L - G213*K23L - G313*K33L + PDstandardNth1K33 - 
-           PDstandardNth3K13) + gu22*(G123*K12L - G122*K13L + G223*K22L - G222*K23L + G323*K23L - G322*K33L + 
-           PDstandardNth2K23 - PDstandardNth3K22) + gu32*
-         (G133*K12L - G123*K13L + G233*K22L - G223*K23L + G333*K23L - G323*K33L + PDstandardNth2K33 - PDstandardNth3K23);
+    HL  =  -2*(Km12*Km21 + Km13*Km31 + Km23*Km32) + trR - SQR(Km11) - SQR(Km22) - 
+        SQR(Km33) + SQR(trK);
+    
+    M1L  =  -2*(gu11*PDstandardNth1K11 + gu12*PDstandardNth1K12 + 
+           gu13*PDstandardNth1K13) + 
+        gu21*(-(G112*K11L) + G111*K12L - G212*K12L - G312*K13L + G211*K22L + 
+           G311*K23L - 3*PDstandardNth1K12 + PDstandardNth2K12) + 
+        gu22*(-(G122*K11L) + G112*K12L - G222*K12L - G322*K13L + G212*K22L + 
+           G312*K23L - 3*PDstandardNth1K22 + PDstandardNth2K22) + 
+        gu23*(-(G123*K11L) + G113*K12L - G223*K12L - G323*K13L + G213*K22L + 
+           G313*K23L - 3*PDstandardNth1K23 + PDstandardNth2K23) + 
+        gu31*(-(G113*K11L) - G213*K12L + G111*K13L - G313*K13L + G211*K23L + 
+           G311*K33L - 3*PDstandardNth1K13 + PDstandardNth3K13) + 
+        gu32*(-(G123*K11L) - G223*K12L + G112*K13L - G323*K13L + G212*K23L + 
+           G312*K33L - 3*PDstandardNth1K23 + PDstandardNth3K23) + 
+        gu33*(-(G133*K11L) - G233*K12L + G113*K13L - G333*K13L + G213*K23L + 
+           G313*K33L - 3*PDstandardNth1K33 + PDstandardNth3K33);
+    
+    M2L  =  gu11*(G112*K11L + (-G111 + G212)*K12L + G312*K13L - G211*K22L - 
+           G311*K23L + PDstandardNth1K11 - 3*PDstandardNth2K11) + 
+        gu12*(G122*K11L + (-G112 + G222)*K12L + G322*K13L - G212*K22L - G312*K23L + 
+           PDstandardNth1K12 - 3*PDstandardNth2K12) + 
+        gu13*(G123*K11L + (-G113 + G223)*K12L + G323*K13L - G213*K22L - G313*K23L + 
+           PDstandardNth1K13 - 3*PDstandardNth2K13) - 
+        2*(gu21*PDstandardNth2K12 + gu22*PDstandardNth2K22 + 
+           gu23*PDstandardNth2K23) + 
+        gu31*(-(G113*K12L) + G112*K13L - G213*K22L + G212*K23L - G313*K23L + 
+           G312*K33L - 3*PDstandardNth2K13 + PDstandardNth3K13) + 
+        gu32*(-(G123*K12L) + G122*K13L - G223*K22L + G222*K23L - G323*K23L + 
+           G322*K33L - 3*PDstandardNth2K23 + PDstandardNth3K23) + 
+        gu33*(-(G133*K12L) + G123*K13L - G233*K22L + G223*K23L - G333*K23L + 
+           G323*K33L - 3*PDstandardNth2K33 + PDstandardNth3K33);
+    
+    M3L  =  gu11*(G113*K11L + G213*K12L + (-G111 + G313)*K13L - G211*K23L - 
+           G311*K33L + PDstandardNth1K11 - 3*PDstandardNth3K11) + 
+        gu12*(G123*K11L + G223*K12L + (-G112 + G323)*K13L - G212*K23L - G312*K33L + 
+           PDstandardNth1K12 - 3*PDstandardNth3K12) + 
+        gu21*(G113*K12L - G112*K13L + G213*K22L - G212*K23L + G313*K23L - 
+           G312*K33L + PDstandardNth2K12 - 3*PDstandardNth3K12) + 
+        gu13*(G133*K11L + G233*K12L + (-G113 + G333)*K13L - G213*K23L - G313*K33L + 
+           PDstandardNth1K13 - 3*PDstandardNth3K13) + 
+        gu22*(G123*K12L - G122*K13L + G223*K22L - G222*K23L + G323*K23L - 
+           G322*K33L + PDstandardNth2K22 - 3*PDstandardNth3K22) + 
+        gu23*(G133*K12L - G123*K13L + G233*K22L - G223*K23L + G333*K23L - 
+           G323*K33L + PDstandardNth2K23 - 3*PDstandardNth3K23) - 
+        2*(gu31*PDstandardNth3K13 + gu32*PDstandardNth3K23 + gu33*PDstandardNth3K33);
     
     
     /* Copy local copies back to grid functions */
diff --git a/ML_ADM/src/ML_ADM_constraints_boundary.c b/ML_ADM/src/ML_ADM_constraints_boundary.c
index 33e504d..4108bff 100644
--- a/ML_ADM/src/ML_ADM_constraints_boundary.c
+++ b/ML_ADM/src/ML_ADM_constraints_boundary.c
@@ -5,7 +5,10 @@
 
 #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"
@@ -84,7 +87,7 @@ void ML_ADM_constraints_boundary_Body(cGH *cctkGH, CCTK_INT dir, CCTK_INT face,
   pm1o12dz2 = -pow(dz,-2)/12.;
   
   /* Loop over the grid points */
-  _Pragma ("omp parallel")
+  #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])
diff --git a/ML_ADM/src/ML_ADM_convertFromADMBase.c b/ML_ADM/src/ML_ADM_convertFromADMBase.c
index a78d3de..8ece38e 100644
--- a/ML_ADM/src/ML_ADM_convertFromADMBase.c
+++ b/ML_ADM/src/ML_ADM_convertFromADMBase.c
@@ -5,7 +5,10 @@
 
 #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"
@@ -84,7 +87,7 @@ void ML_ADM_convertFromADMBase_Body(cGH *cctkGH, CCTK_INT dir, CCTK_INT face, CC
   pm1o12dz2 = -pow(dz,-2)/12.;
   
   /* Loop over the grid points */
-  _Pragma ("omp parallel")
+  #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])
diff --git a/ML_ADM/src/ML_ADM_convertToADMBase.c b/ML_ADM/src/ML_ADM_convertToADMBase.c
index 50150f1..a7047e3 100644
--- a/ML_ADM/src/ML_ADM_convertToADMBase.c
+++ b/ML_ADM/src/ML_ADM_convertToADMBase.c
@@ -5,7 +5,10 @@
 
 #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"
@@ -84,7 +87,7 @@ void ML_ADM_convertToADMBase_Body(cGH *cctkGH, CCTK_INT dir, CCTK_INT face, CCTK
   pm1o12dz2 = -pow(dz,-2)/12.;
   
   /* Loop over the grid points */
-  _Pragma ("omp parallel")
+  #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])