MoL: add Multirate capabilities. This add three new multirate RK schemes to MoL.

Flags indicate whether it is time to execute slow RHS computation.
For instance, in the RK4-RK2 scheme, there are 4 substeps in total, but the RK2 RHS are only evaluated in the very first and in the very last step of the four substeps.

From: Christian Reisswig, minor changes by Roland Haas



git-svn-id: http://svn.cactuscode.org/arrangements/CactusNumerical/MoL/trunk@175 578cdeb0-5ea1-4b81-8215-5a3b8777ee0b

18 files changed, 2336 insertions, 5 deletions

diff --git a/interface.ccl b/interface.ccl
index fa72e0c..d2d14f9 100644
--- a/interface.ccl
+++ b/interface.ccl
@@ -34,14 +34,21 @@ USES FUNCTION EnableProlongating
 
 CCTK_INT FUNCTION MoLRegisterEvolved(CCTK_INT IN EvolvedIndex, \
                                      CCTK_INT IN RHSIndex)
+CCTK_INT FUNCTION MoLRegisterEvolvedSlow(CCTK_INT IN EvolvedIndex, \
+                                         CCTK_INT IN RHSIndexSlow)
+
 CCTK_INT FUNCTION MoLRegisterConstrained(CCTK_INT IN ConstrainedIndex)
 CCTK_INT FUNCTION MoLRegisterSaveAndRestore(CCTK_INT IN SandRIndex)
 CCTK_INT FUNCTION MoLRegisterEvolvedGroup(CCTK_INT IN EvolvedIndex, \
                                           CCTK_INT IN RHSIndex)
+CCTK_INT FUNCTION MoLRegisterEvolvedGroupSlow(CCTK_INT IN EvolvedIndex, \
+                                              CCTK_INT IN RHSIndexSlow)
 CCTK_INT FUNCTION MoLRegisterConstrainedGroup(CCTK_INT IN ConstrainedIndex)
 CCTK_INT FUNCTION MoLRegisterSaveAndRestoreGroup(CCTK_INT IN SandRIndex)
 CCTK_INT FUNCTION MoLChangeToEvolved(CCTK_INT IN EvolvedIndex, \
                                      CCTK_INT IN RHSIndex)
+CCTK_INT FUNCTION MoLChangeToEvolvedSlow(CCTK_INT IN EvolvedIndex, \
+                                         CCTK_INT IN RHSIndexSlow)
 CCTK_INT FUNCTION MoLChangeToConstrained(CCTK_INT IN ConstrainedIndex)
 CCTK_INT FUNCTION MoLChangeToSaveAndRestore(CCTK_INT IN SandRIndex)
 CCTK_INT FUNCTION MoLChangeToNone(CCTK_INT IN RemoveIndex)
@@ -49,17 +56,21 @@ CCTK_INT FUNCTION MoLQueryEvolvedRHS(CCTK_INT IN EvolvedIndex)
 CCTK_INT FUNCTION MoLNumIntegratorSubsteps()
 
 PROVIDES FUNCTION MoLRegisterEvolved WITH MoL_RegisterEvolved LANGUAGE C
+PROVIDES FUNCTION MoLRegisterEvolvedSlow WITH MoL_RegisterEvolvedSlow LANGUAGE C
 PROVIDES FUNCTION MoLRegisterConstrained WITH MoL_RegisterConstrained \
                   LANGUAGE C
 PROVIDES FUNCTION MoLRegisterSaveAndRestore WITH MoL_RegisterSaveAndRestore \
                   LANGUAGE C
 PROVIDES FUNCTION MoLRegisterEvolvedGroup WITH MoL_RegisterEvolvedGroup \
                   LANGUAGE C
+PROVIDES FUNCTION MoLRegisterEvolvedGroupSlow WITH MoL_RegisterEvolvedGroupSlow \
+                  LANGUAGE C
 PROVIDES FUNCTION MoLRegisterConstrainedGroup WITH \
                   MoL_RegisterConstrainedGroup LANGUAGE C
 PROVIDES FUNCTION MoLRegisterSaveAndRestoreGroup WITH \
                   MoL_RegisterSaveAndRestoreGroup LANGUAGE C
 PROVIDES FUNCTION MoLChangeToEvolved WITH MoL_ChangeToEvolved LANGUAGE C
+PROVIDES FUNCTION MoLChangeToEvolvedSlow WITH MoL_ChangeToEvolvedSlow LANGUAGE C
 PROVIDES FUNCTION MoLChangeToConstrained WITH MoL_ChangeToConstrained \
                   LANGUAGE C
 PROVIDES FUNCTION MoLChangeToSaveAndRestore WITH MoL_ChangeToSaveAndRestore \
@@ -226,6 +237,17 @@ CCTK_INT MoL_Counters \
 {
   MoL_Intermediate_Step
   MoL_Stepsize_Bad
+  
+  # A flag indicating whether it is time for slow RHS evaluation. 
+  # Oustide the MoL loop, it is guaranteed to be 1.
+  # It is only zero for certain MoL substeps when multirate methods are used.
+  MoL_SlowStep
+  
+  # A flag indicating whether it is time for slow post step computations (e.g. applying BCs)
+  # Oustide the MoL loop, it is guaranteed to be 1.
+  # It is only zero for certain MoL substeps when multirate methods are used.
+  MoL_SlowPostStep
+  
 } "The counter for the time integration method"
 
 CCTK_REAL MoL_Original_Time \
@@ -241,6 +263,11 @@ CCTK_REAL ScratchSpace[MoL_Num_Evolved_Vars*MoL_Num_Scratch_Levels] \
             Timelevels = 1 \
             TAGS       = 'Prolongation="None" Checkpoint="no"'
 
+CCTK_REAL ScratchSpaceSlow[MoL_Num_Evolved_Vars_Slow*MoL_Num_Scratch_Levels] \
+            TYPE       = GF \
+            Timelevels = 1 \
+            TAGS       = 'Prolongation="None" Checkpoint="no"'
+
 CCTK_REAL SandRScratchSpace[MoL_Num_SaveAndRestore_Vars] \
             TYPE       = GF \
             Timelevels = 1 \
diff --git a/param.ccl b/param.ccl
index ba43b51..e82dd07 100644
--- a/param.ccl
+++ b/param.ccl
@@ -8,6 +8,12 @@ CCTK_INT MoL_Num_Evolved_Vars "The maximum number of variables to be evolved by
   0:*		:: "Anything non negative. Added to by other thorns."
 } 0
 
+CCTK_INT MoL_Num_Evolved_Vars_Slow "The maximum number of 'slow' variables to be evolved by MoL" ACCUMULATOR = (x+y)
+{
+  0:*		:: "Anything non negative. Added to by other thorns."
+} 0
+
+
 CCTK_INT MoL_Num_Constrained_Vars "The maximum number of constrained variables with timelevels that MoL needs to know about" ACCUMULATOR = (x+y)
 {
   0:*		:: "Anything non negative. Added to by other thorns."
@@ -92,6 +98,9 @@ KEYWORD ODE_Method "The ODE method use by MoL to do time integration"
   "RK45CK"      :: "RK45CK (Cash-Karp) with error estimation"
   "RK65"        :: "RK65 with error estimation"
   "RK87"        :: "RK87 with error estimation"
+  "RK2-MR-2:1"  :: "2nd order 2:1 multirate RK scheme based on RK2 due to Schlegel et al 2009"
+  "RK4-MR-2:1"  :: "3rd order 2:1 multirate RK scheme based on RK43 due to Schlegel et al 2009"
+  "RK4-RK2"     :: "RK4 as fast method and RK2 as slow method"
 } "ICN"
 
 KEYWORD Generic_Type "If using the generic method, which sort"
diff --git a/schedule.ccl b/schedule.ccl
index f0aed52..8dc65f9 100644
--- a/schedule.ccl
+++ b/schedule.ccl
@@ -399,7 +399,13 @@ else
     }
   }
 }
-    
+
+
+if (MoL_Num_Evolved_Vars_Slow)
+{
+   STORAGE: ScratchSpaceSlow
+}
+
 ######################################################
 ###  StartStep contains the routines that just do  ###
 ###  internal MoL stuff; setting the counter, and  ###
@@ -607,6 +613,28 @@ else if (CCTK_Equals(ODE_Method,"ICN-avg"))
     LANG: C
   } "Updates calculated with the averaging ICN method"
 }
+else if (CCTK_Equals(ODE_Method,"RK2-MR-2:1")) 
+{
+  schedule MoL_RK2_MR_2_1_Add AS MoL_Add IN MoL_Step AFTER MoL_CalcRHS BEFORE MoL_PostStep
+  {
+    LANG: C
+  } "Updates calculated with the multirate Runge-Kutta 2 method"
+}
+else if (CCTK_Equals(ODE_Method,"RK4-MR-2:1")) 
+{
+  schedule MoL_RK4_MR_2_1_Add AS MoL_Add IN MoL_Step AFTER MoL_CalcRHS BEFORE MoL_PostStep
+  {
+    LANG: C
+  } "Updates calculated with the multirate Runge-Kutta 4 method"
+}
+else if (CCTK_Equals(ODE_Method,"RK4-RK2")) 
+{
+  schedule MoL_RK4_RK2_Add AS MoL_Add IN MoL_Step AFTER MoL_CalcRHS BEFORE MoL_PostStep
+  {
+    LANG: C
+  } "Updates calculated with the multirate RK4/RK2 method"
+}
+
 
 ##################################################
 ###  Physics thorns can apply boundaries and   ###
diff --git a/src/ChangeType.c b/src/ChangeType.c
index e0bc333..9c1db4a 100644
--- a/src/ChangeType.c
+++ b/src/ChangeType.c
@@ -27,6 +27,7 @@ CCTK_FILEVERSION(CactusBase_MoL_ChangeType_c);
 #define MOL_EVOLVED_VARTYPE 1
 #define MOL_CONSTRAINED_VARTYPE 2
 #define MOL_SANDR_VARTYPE 3
+#define MOL_EVOLVEDSLOW_VARTYPE 4
 
 /********************************************************************
  ********************* Local Routine Prototypes *********************
@@ -42,6 +43,8 @@ CCTK_FILEVERSION(CactusBase_MoL_ChangeType_c);
 
 CCTK_INT MoL_ChangeToEvolved(CCTK_INT EvolvedIndex, CCTK_INT RHSIndex);
 
+CCTK_INT MoL_ChangeToEvolvedSlow(CCTK_INT EvolvedIndex, CCTK_INT RHSIndexSlow);
+
 CCTK_INT MoL_ChangeToConstrained(CCTK_INT ConstrainedIndex);
 
 CCTK_INT MoL_ChangeToSaveAndRestore(CCTK_INT SandRIndex);
@@ -95,6 +98,16 @@ CCTK_INT MoL_ChangeToEvolved(CCTK_INT EvolvedIndex, CCTK_INT RHSIndex)
     }
   }
   
+  for (index = 0; (index < MoLNumEvolvedVariablesSlow)&&(!vartype); index++)
+  {
+    vartype = MOL_EVOLVEDSLOW_VARTYPE * 
+      (EvolvedVariableIndexSlow[index] == EvolvedIndex);
+    if (vartype)
+    {
+      usedindex = index;
+    }
+  }
+  
   for (index = 0; (index < MoLNumConstrainedVariables)&&(!vartype); index++)
   {
     vartype = MOL_CONSTRAINED_VARTYPE * 
@@ -174,6 +187,46 @@ CCTK_INT MoL_ChangeToEvolved(CCTK_INT EvolvedIndex, CCTK_INT RHSIndex)
         break;
       }
 
+    case MOL_EVOLVEDSLOW_VARTYPE:
+      
+      {
+        timelevs = CCTK_MaxTimeLevelsVI(RHSIndex);
+        if (timelevs < 1) {
+          CCTK_VWarn(1,__LINE__,__FILE__,CCTK_THORNSTRING,
+                     "Warning whilst trying to change variable "
+                     "index %i to evolved type from constrained.", RHSIndex);
+          CCTK_WARN(0, "The RHS index passed does not correspond to a GF.");
+        }
+
+        if (!(MoLNumEvolvedVariables < MoL_Num_Evolved_Vars))
+        {
+          CCTK_VWarn(1,__LINE__,__FILE__,CCTK_THORNSTRING,
+                     "Warning whilst trying to change variable "
+                     "index %i (%s) to evolved.", 
+                     EvolvedIndex, CCTK_VarName(EvolvedIndex));
+          CCTK_WARN(0, "When changing type there are more variables "
+                    "than the accumulator parameter "
+                    "MoL_Num_Evolved_Vars allows. Check that "
+                    "you are accumulating onto this parameter correctly");
+        }
+
+        for (index = usedindex; index < MoLNumEvolvedVariablesSlow - 1; 
+             index++)
+        {
+          EvolvedVariableIndexSlow[index] = EvolvedVariableIndexSlow[index+1];
+        }
+        MoLNumEvolvedVariablesSlow--;
+        EvolvedVariableIndex[MoLNumEvolvedVariables] = EvolvedIndex;
+        RHSVariableIndex[MoLNumEvolvedVariables] = RHSIndex;
+        MoLNumEvolvedVariables++;
+#ifdef MOLDEBUG
+        printf("Changing (constrained): vars %d var %d (%s).\n",
+               MoLNumEvolvedVariables, EvolvedIndex,
+               CCTK_VarName(EvolvedIndex));
+#endif
+        break;
+      }
+      
     case MOL_CONSTRAINED_VARTYPE:
       
       {
@@ -268,6 +321,266 @@ CCTK_INT MoL_ChangeToEvolved(CCTK_INT EvolvedIndex, CCTK_INT RHSIndex)
 }
 
  /*@@
+   @routine    MoL_ChangeToEvolvedSlow
+   @date       Thu May 30 16:45:30 2002
+   @author     Ian Hawke, Roland Haas
+   @desc 
+   Changes a variable to evolved slow type. Checks to see which type it was
+   before and does the bookkeeping on the index arrays.
+   @enddesc 
+   @calls     
+   @calledby   
+   @history 
+ 
+   @endhistory 
+
+@@*/
+
+CCTK_INT MoL_ChangeToEvolvedSlow(CCTK_INT EvolvedIndex, CCTK_INT RHSIndex)
+{
+
+  DECLARE_CCTK_PARAMETERS;
+
+  CCTK_INT index, usedindex;
+  CCTK_INT vartype; /* See the defines at the top of file */
+  CCTK_INT timelevs;
+  
+  vartype = 0;
+  usedindex = -1;
+  
+  for (index = 0; (index < MoLNumEvolvedVariables)&&(!vartype); index++)
+  {
+    vartype = MOL_EVOLVED_VARTYPE * 
+      (EvolvedVariableIndex[index] == EvolvedIndex);
+    if (vartype)
+    {
+      usedindex = index;
+    }
+  }
+  
+  for (index = 0; (index < MoLNumEvolvedVariablesSlow)&&(!vartype); index++)
+  {
+    vartype = MOL_EVOLVEDSLOW_VARTYPE * 
+      (EvolvedVariableIndexSlow[index] == EvolvedIndex);
+    if (vartype)
+    {
+      usedindex = index;
+    }
+  }
+  
+  for (index = 0; (index < MoLNumConstrainedVariables)&&(!vartype); index++)
+  {
+    vartype = MOL_CONSTRAINED_VARTYPE * 
+      (ConstrainedVariableIndex[index] == EvolvedIndex);
+    if (vartype)
+    {
+      usedindex = index;
+    }
+  }
+  
+  for (index = 0; (index < MoLNumSandRVariables)&&(!vartype); index++)
+  {
+    vartype = MOL_SANDR_VARTYPE * 
+      (SandRVariableIndex[index] == EvolvedIndex);
+    if (vartype)
+    {
+      usedindex = index;
+    }
+  }
+
+  switch (vartype)
+  {
+    
+    case MOL_UNKNOWN_VARTYPE:
+
+      {
+        timelevs = CCTK_MaxTimeLevelsVI(EvolvedIndex);
+        if (timelevs < 1)
+        {
+          CCTK_VWarn(1,__LINE__,__FILE__,CCTK_THORNSTRING,
+                     "Warning whilst trying to change variable "
+                     "index %i to slow evolved.", EvolvedIndex);
+          CCTK_WARN(0, "The index passed does not correspond to a GF.");
+        }
+        else if (timelevs == 1) 
+        {
+          CCTK_VWarn(1,__LINE__,__FILE__,CCTK_THORNSTRING,
+                     "Warning whilst trying to change variable "
+                     "index %i to slow evolved.", EvolvedIndex);
+          CCTK_WARN(0, "The index passed only has a single timelevel.");
+        }
+        timelevs = CCTK_MaxTimeLevelsVI(RHSIndex);
+        if (timelevs < 1) {
+          CCTK_VWarn(1,__LINE__,__FILE__,CCTK_THORNSTRING,
+                     "Warning whilst trying to change variable "
+                     "index %i to slow evolved (RHS GF).", RHSIndex);
+          CCTK_WARN(0, "The index passed does not correspond to a GF.");
+        }
+
+        if (!(MoLNumEvolvedVariablesSlow < MoL_Num_Evolved_Vars_Slow))
+        {
+          CCTK_VWarn(1,__LINE__,__FILE__,CCTK_THORNSTRING,
+                     "Warning whilst trying to change variable "
+                     "index %i (%s) to slow evolved.", 
+                     EvolvedIndex, CCTK_VarName(EvolvedIndex));
+          CCTK_WARN(0, "When changing type there are more variables "
+                    "than the accumulator parameter "
+                    "MoL_Num_Evolved_Vars_Slow allows. Check that "
+                    "you are accumulating onto this parameter correctly");
+        }
+
+        EvolvedVariableIndexSlow[MoLNumEvolvedVariablesSlow] = EvolvedIndex;
+        RHSVariableIndexSlow[MoLNumEvolvedVariablesSlow] = RHSIndex;
+        MoLNumEvolvedVariablesSlow++;
+#ifdef MOLDEBUG
+        printf("Changing (unknown): vars %d var %d (%s).\n",
+               MoLNumEvolvedVariablesSlow, EvolvedIndex,
+               CCTK_VarName(EvolvedIndex));
+#endif
+        break;
+      }
+      
+    case MOL_EVOLVED_VARTYPE:
+      
+      {
+        timelevs = CCTK_MaxTimeLevelsVI(RHSIndex);
+        if (timelevs < 1) {
+          CCTK_VWarn(1,__LINE__,__FILE__,CCTK_THORNSTRING,
+                     "Warning whilst trying to change variable "
+                     "index %i to slow evolved type from constrained.", RHSIndex);
+          CCTK_WARN(0, "The RHS index passed does not correspond to a GF.");
+        }
+
+        if (!(MoLNumEvolvedVariablesSlow < MoL_Num_Evolved_Vars_Slow))
+        {
+          CCTK_VWarn(1,__LINE__,__FILE__,CCTK_THORNSTRING,
+                     "Warning whilst trying to change variable "
+                     "index %i (%s) to slow evolved.", 
+                     EvolvedIndex, CCTK_VarName(EvolvedIndex));
+          CCTK_WARN(0, "When changing type there are more variables "
+                    "than the accumulator parameter "
+                    "MoL_Num_Evolved_Vars_Slow allows. Check that "
+                    "you are accumulating onto this parameter correctly");
+        }
+
+        for (index = usedindex; index < MoLNumEvolvedVariables - 1; 
+             index++)
+        {
+          EvolvedVariableIndex[index] = EvolvedVariableIndex[index+1];
+        }
+        MoLNumEvolvedVariables--;
+        EvolvedVariableIndexSlow[MoLNumEvolvedVariablesSlow] = EvolvedIndex;
+        RHSVariableIndexSlow[MoLNumEvolvedVariablesSlow] = RHSIndex;
+        MoLNumEvolvedVariablesSlow++;
+#ifdef MOLDEBUG
+        printf("Changing (constrained): vars %d var %d (%s).\n",
+               MoLNumEvolvedVariablesSlow, EvolvedIndex,
+               CCTK_VarName(EvolvedIndex));
+#endif
+        break;
+      }
+      
+    case MOL_EVOLVEDSLOW_VARTYPE:
+      
+      {
+        RHSVariableIndex[usedindex] = RHSIndex;
+        break;
+      }
+
+    case MOL_CONSTRAINED_VARTYPE:
+      
+      {
+        timelevs = CCTK_MaxTimeLevelsVI(RHSIndex);
+        if (timelevs < 1) {
+          CCTK_VWarn(1,__LINE__,__FILE__,CCTK_THORNSTRING,
+                     "Warning whilst trying to change variable "
+                     "index %i to slow evolved type from constrained.", RHSIndex);
+          CCTK_WARN(0, "The RHS index passed does not correspond to a GF.");
+        }
+
+        if (!(MoLNumEvolvedVariablesSlow < MoL_Num_Evolved_Vars_Slow))
+        {
+          CCTK_VWarn(1,__LINE__,__FILE__,CCTK_THORNSTRING,
+                     "Warning whilst trying to change variable "
+                     "index %i (%s) to slow evolved.", 
+                     EvolvedIndex, CCTK_VarName(EvolvedIndex));
+          CCTK_WARN(0, "When changing type there are more variables "
+                    "than the accumulator parameter "
+                    "MoL_Num_Evolved_Vars_Slow allows. Check that "
+                    "you are accumulating onto this parameter correctly");
+        }
+
+        for (index = usedindex; index < MoLNumConstrainedVariables - 1; 
+             index++)
+        {
+          ConstrainedVariableIndex[index] = ConstrainedVariableIndex[index+1];
+        }
+        MoLNumConstrainedVariables--;
+        EvolvedVariableIndexSlow[MoLNumEvolvedVariablesSlow] = EvolvedIndex;
+        RHSVariableIndexSlow[MoLNumEvolvedVariablesSlow] = RHSIndex;
+        MoLNumEvolvedVariablesSlow++;
+#ifdef MOLDEBUG
+        printf("Changing (constrained): vars %d var %d (%s).\n",
+               MoLNumEvolvedVariables, EvolvedIndex,
+               CCTK_VarName(EvolvedIndex));
+#endif
+        break;
+      }
+      
+    case MOL_SANDR_VARTYPE:
+      
+      {
+        timelevs = CCTK_MaxTimeLevelsVI(RHSIndex);
+        if (timelevs < 1) {
+          CCTK_VWarn(1,__LINE__,__FILE__,CCTK_THORNSTRING,
+                     "Warning whilst trying to change variable "
+                     "index %i to slow evolved type from save and "
+                     "restore.", RHSIndex);
+          CCTK_WARN(0, "The RHS index passed does not correspond to a GF.");
+        }
+
+        if (!(MoLNumEvolvedVariablesSlow < MoL_Num_Evolved_Vars_Slow))
+        {
+          CCTK_VWarn(1,__LINE__,__FILE__,CCTK_THORNSTRING,
+                     "Warning whilst trying to change variable "
+                     "index %i (%s) to slow evolved.", 
+                     EvolvedIndex, CCTK_VarName(EvolvedIndex));
+          CCTK_WARN(0, "When changing type there are more variables "
+                    "than the accumulator parameter "
+                    "MoL_Num_Evolved_Vars_Slow allows. Check that "
+                    "you are accumulating onto this parameter correctly");
+        }
+
+        for (index = usedindex; index < MoLNumSandRVariables - 1; index++)
+        {
+          SandRVariableIndex[index] = SandRVariableIndex[index+1];
+        }
+        MoLNumSandRVariables--;
+        EvolvedVariableIndexSlow[MoLNumEvolvedVariablesSlow] = EvolvedIndex;
+        RHSVariableIndexSlow[MoLNumEvolvedVariablesSlow] = RHSIndex;
+        MoLNumEvolvedVariablesSlow++;
+#ifdef MOLDEBUG
+        printf("Changing (SandR): vars %d var %d (%s).\n",
+               MoLNumEvolvedVariablesSlow, EvolvedIndex,
+               CCTK_VarName(EvolvedIndex));
+#endif
+        break;
+      }
+      
+    default:
+      
+      {
+        CCTK_WARN(0, "Something is seriously wrong in ChangeType.c! "
+                  "Case out of range in switch statement.");
+      }
+      
+  }
+      
+  return 0;
+  
+}
+
+ /*@@
    @routine    MoL_ChangeToConstrained
    @date       Thu May 30 16:47:08 2002
    @author     Ian Hawke
@@ -321,6 +634,16 @@ CCTK_INT MoL_ChangeToConstrained(CCTK_INT ConstrainedIndex)
     }
   }
   
+  for (index = 0; (index < MoLNumEvolvedVariablesSlow)&&(!vartype); index++)
+  {
+    vartype = MOL_EVOLVEDSLOW_VARTYPE * 
+      (EvolvedVariableIndexSlow[index] == ConstrainedIndex);
+    if (vartype)
+    {
+      usedindex = index;
+    }
+  }
+  
   for (index = 0; (index < MoLNumConstrainedVariables)&&(!vartype); index++)
   {
     vartype = MOL_CONSTRAINED_VARTYPE * 
@@ -397,6 +720,34 @@ CCTK_INT MoL_ChangeToConstrained(CCTK_INT ConstrainedIndex)
         break;
       }
 
+    case MOL_EVOLVEDSLOW_VARTYPE:
+      
+      {
+        for (index = usedindex; index < MoLNumEvolvedVariablesSlow - 1; index++)
+        {
+          EvolvedVariableIndex[index] = EvolvedVariableIndex[index+1];
+          RHSVariableIndex[index] = RHSVariableIndex[index+1];
+        }
+        MoLNumEvolvedVariablesSlow--;
+        if (timelevs > 1)
+        {
+          if ( !(MoLNumConstrainedVariables < MoL_Num_Constrained_Vars))
+          {
+            CCTK_VWarn(1,__LINE__,__FILE__,CCTK_THORNSTRING,
+                       "Warning whilst trying to change variable "
+                       "index %i (%s) to constrained.", 
+                       ConstrainedIndex, CCTK_VarName(ConstrainedIndex));
+            CCTK_WARN(0, "When changing type there are more variables "
+                      "than the accumulator parameter "
+                      "MoL_Num_Constrained_Vars allows. Check that "
+                      "you are accumulating onto this parameter correctly");
+          }
+          ConstrainedVariableIndex[MoLNumConstrainedVariables] = ConstrainedIndex;
+          MoLNumConstrainedVariables++;
+        }
+        break;
+      }
+
     case MOL_CONSTRAINED_VARTYPE:
       
       {
@@ -489,6 +840,16 @@ CCTK_INT MoL_ChangeToSaveAndRestore(CCTK_INT SandRIndex)
     }
   }
   
+  for (index = 0; (index < MoLNumEvolvedVariablesSlow)&&(!vartype); index++)
+  {
+    vartype = MOL_EVOLVEDSLOW_VARTYPE * 
+      (EvolvedVariableIndexSlow[index] == SandRIndex);
+    if (vartype)
+    {
+      usedindex = index;
+    }
+  }
+  
   for (index = 0; (index < MoLNumConstrainedVariables)&&(!vartype); index++)
   {
     vartype = MOL_CONSTRAINED_VARTYPE * 
@@ -558,6 +919,32 @@ CCTK_INT MoL_ChangeToSaveAndRestore(CCTK_INT SandRIndex)
         break;
       }
 
+    case MOL_EVOLVEDSLOW_VARTYPE:
+      
+      {
+
+        if (!(MoLNumSandRVariables < MoL_Num_SaveAndRestore_Vars))
+        {
+          CCTK_VWarn(1,__LINE__,__FILE__,CCTK_THORNSTRING,
+                     "Warning whilst trying to change variable "
+                     "index %i (%s) to save and restore.", 
+                     SandRIndex, CCTK_VarName(SandRIndex));
+          CCTK_WARN(0, "When changing type there are more variables "
+                    "than the accumulator parameter "
+                    "MoL_Num_SaveAndRestore_Vars allows. Check that "
+                    "you are accumulating onto this parameter correctly");
+        }
+        for (index = usedindex; index < MoLNumEvolvedVariablesSlow - 1; index++)
+        {
+          EvolvedVariableIndex[index] = EvolvedVariableIndex[index+1];
+          RHSVariableIndex[index] = RHSVariableIndex[index+1];
+        }
+        MoLNumEvolvedVariablesSlow--;
+        SandRVariableIndex[MoLNumSandRVariables] = SandRIndex;
+        MoLNumSandRVariables++;
+        break;
+      }
+
     case MOL_CONSTRAINED_VARTYPE:
       
       {
@@ -639,6 +1026,16 @@ CCTK_INT MoL_ChangeToNone(CCTK_INT RemoveIndex)
     }
   }
   
+  for (index = 0; (index < MoLNumEvolvedVariablesSlow)&&(!vartype); index++)
+  {
+    vartype = MOL_EVOLVEDSLOW_VARTYPE * 
+      (EvolvedVariableIndexSlow[index] == RemoveIndex);
+    if (vartype)
+    {
+      usedindex = index;
+    }
+  }
+  
   for (index = 0; (index < MoLNumConstrainedVariables)&&(!vartype); index++)
   {
     vartype = MOL_CONSTRAINED_VARTYPE * 
@@ -680,6 +1077,18 @@ CCTK_INT MoL_ChangeToNone(CCTK_INT RemoveIndex)
         break;
       }
 
+    case MOL_EVOLVEDSLOW_VARTYPE:
+      
+      {
+        for (index = usedindex; index < MoLNumEvolvedVariablesSlow - 1; index++)
+        {
+          EvolvedVariableIndex[index] = EvolvedVariableIndex[index+1];
+          RHSVariableIndex[index] = RHSVariableIndex[index+1];
+        }
+        MoLNumEvolvedVariablesSlow--;
+        break;
+      }
+
     case MOL_CONSTRAINED_VARTYPE:
       
       {
diff --git a/src/Counter.c b/src/Counter.c
index 44515de..0093f1f 100644
--- a/src/Counter.c
+++ b/src/Counter.c
@@ -69,7 +69,13 @@ void MoL_SetCounter(CCTK_ARGUMENTS)
 
   *MoL_Intermediate_Step = MoL_Intermediate_Steps;
 
-   
+  // Always indicate execution of slow steps.
+  // If multirate methods are used, they are set to zero
+  // for certain substeps.
+  // In any case, outside the MoL loop, these scalars are guaranteed to be 1
+  *MoL_SlowStep = 1;
+  *MoL_SlowPostStep = 1;
+  
 /* #ifdef HAVE_CARPET */
   if ((*MoL_Intermediate_Step))
   {
@@ -114,6 +120,71 @@ void MoL_DecrementCounter(CCTK_ARGUMENTS)
   
   (*MoL_Intermediate_Step) --;
 
+
+  if (CCTK_EQUALS(ODE_Method, "RK2-MR-2:1"))
+  {
+     switch (MoL_Intermediate_Steps - *MoL_Intermediate_Step)
+     {
+        case 2:
+        case 4:
+           if (*MoL_SlowStep)
+              *MoL_SlowPostStep = 1;
+           else
+              *MoL_SlowPostStep = 0;
+           // don't do slow step
+           *MoL_SlowStep = 0;
+           break;
+        default:
+           if (!(*MoL_SlowStep))
+              *MoL_SlowPostStep = 0;
+           else
+              *MoL_SlowPostStep = 1;
+           // do a slow step!
+           *MoL_SlowStep = 1;
+     }
+  }
+
+  if (CCTK_EQUALS(ODE_Method, "RK4-MR-2:1"))
+  {
+     switch (MoL_Intermediate_Steps - *MoL_Intermediate_Step)
+     {
+        case 2:
+        case 4:
+        case 7:
+        case 9:
+           // don't do slow step
+           *MoL_SlowStep = 0;
+           break;
+        default:
+           // do a slow step!
+           *MoL_SlowStep = 1;
+     }
+  }
+
+  if (CCTK_EQUALS(ODE_Method, "RK4-RK2"))
+  {
+     switch (MoL_Intermediate_Steps - *MoL_Intermediate_Step)
+     {
+        case 0:
+           // do a slow step!
+           *MoL_SlowStep = 1;
+           // don't sync!
+           *MoL_SlowPostStep = 0;
+           break;
+        case 1:
+        case 2:
+           // don't do slow step
+           *MoL_SlowStep = 0;
+           *MoL_SlowPostStep = 0;
+           break;
+        case 3:
+           // do a slow step!
+           *MoL_SlowStep = 1;
+           // sync!
+           *MoL_SlowPostStep = 1;
+     }
+  }
+
 /* #ifdef HAVE_CARPET */
   if (! (*MoL_Intermediate_Step))
   {
diff --git a/src/ExternalVariables.h b/src/ExternalVariables.h
index 54a5129..5f8c1e6 100644
--- a/src/ExternalVariables.h
+++ b/src/ExternalVariables.h
@@ -16,12 +16,15 @@
 
 
 extern CCTK_INT *EvolvedVariableIndex;
+extern CCTK_INT *EvolvedVariableIndexSlow;
 extern CCTK_INT *RHSVariableIndex;
+extern CCTK_INT *RHSVariableIndexSlow;
 extern CCTK_INT *ConstrainedVariableIndex;
 extern CCTK_INT *SandRVariableIndex;
 
 
 extern CCTK_INT MoLNumEvolvedVariables;
+extern CCTK_INT MoLNumEvolvedVariablesSlow;
 extern CCTK_INT MoLNumConstrainedVariables;
 extern CCTK_INT MoLNumSandRVariables;
 
diff --git a/src/IndexArrays.c b/src/IndexArrays.c
index d3ceccb..485a4bd 100644
--- a/src/IndexArrays.c
+++ b/src/IndexArrays.c
@@ -110,6 +110,23 @@ void MoL_SetupIndexArrays(CCTK_ARGUMENTS)
       }
   }
 
+  if (MoL_Num_Evolved_Vars_Slow)
+  {
+    EvolvedVariableIndexSlow = (CCTK_INT *)malloc(MoL_Num_Evolved_Vars_Slow * 
+					      sizeof(CCTK_INT));
+    if (!EvolvedVariableIndexSlow)
+      {
+	CCTK_WARN(0,"Failed to allocate the slow evolved variable index array");
+      }
+
+    RHSVariableIndexSlow = (CCTK_INT *)malloc(MoL_Num_Evolved_Vars_Slow * 
+					  sizeof(CCTK_INT));
+    if (!RHSVariableIndexSlow)
+      {
+	CCTK_WARN(0,"Failed to allocate the slow RHS variable index array");
+      }
+  }
+
   if (MoL_Num_Constrained_Vars)
   {
     ConstrainedVariableIndex = (CCTK_INT *)malloc(MoL_Num_Constrained_Vars * 
@@ -335,6 +352,18 @@ void MoL_SetupIndexArrays(CCTK_ARGUMENTS)
             "Averaging iterative Crank Nicholson with %i iterations", 
             MoL_Intermediate_Steps);
   }
+  else if (CCTK_EQUALS(ODE_Method,"RK2-MR-2:1")) 
+  {
+    sprintf(infoline, "Multi-rate 2:1 Runge-Kutta 2");
+  }
+  else if (CCTK_EQUALS(ODE_Method,"RK4-MR-2:1")) 
+  {
+    sprintf(infoline, "Multi-rate 2:1 Runge-Kutta 4");
+  }
+  else if (CCTK_EQUALS(ODE_Method,"RK4-RK2")) 
+  {
+    sprintf(infoline, "Multi-rate 2:1 Runge-Kutta 4 and Runge-Kutta 2");
+  }
   else 
   {
     CCTK_WARN(0, "ODE_Method not recognized!");
@@ -345,6 +374,12 @@ void MoL_SetupIndexArrays(CCTK_ARGUMENTS)
   free(infoline);
   infoline = NULL;
 
+  // These scalars must be 1 oustide of the MoL loop.
+  // They will only be zero for certain substeps when multirate methods are used.
+  // Otherwise, they are guaranteed to always be ONE.
+  *MoL_SlowPostStep = 1;
+  *MoL_SlowStep = 1;
+
   return;
 
 }
@@ -381,6 +416,18 @@ void MoL_FreeIndexArrays(CCTK_ARGUMENTS)
     RHSVariableIndex = NULL;
   }
 
+  if (EvolvedVariableIndexSlow)
+  {
+    free(EvolvedVariableIndexSlow);
+    EvolvedVariableIndexSlow = NULL;
+  }
+
+  if (RHSVariableIndexSlow)
+  {
+    free(RHSVariableIndexSlow);  
+    RHSVariableIndexSlow = NULL;
+  }
+
   if (ConstrainedVariableIndex)
   {
     free(ConstrainedVariableIndex);
diff --git a/src/MoL.h b/src/MoL.h
index 2135870..707c7e2 100644
--- a/src/MoL.h
+++ b/src/MoL.h
@@ -20,6 +20,7 @@ extern "C" {
 #ifdef CCODE
 
 CCTK_INT MoL_RegisterEvolved(CCTK_INT EvolvedIndex, CCTK_INT RHSIndex);
+CCTK_INT MoL_RegisterEvolvedSlow(CCTK_INT EvolvedIndex, CCTK_INT RHSIndexSlow);
 CCTK_INT MoL_RegisterConstrained(CCTK_INT ConstrainedIndex);
 CCTK_INT MoL_RegisterSaveAndRestore(CCTK_INT SandRIndex);
 CCTK_INT MoL_ChangeToEvolved(CCTK_INT EvolvedIndex, CCTK_INT RHSIndex);
@@ -28,6 +29,8 @@ CCTK_INT MoL_ChangeToSaveAndRestore(CCTK_INT SandRIndex);
 CCTK_INT MoL_ChangeToNone(CCTK_INT RemoveIndex);
 CCTK_INT MoL_RegisterEvolvedGroup(CCTK_INT EvolvedGroupIndex, 
                                   CCTK_INT RHSGroupIndex);
+CCTK_INT MoL_RegisterEvolvedGroupSlow(CCTK_INT EvolvedGroupIndexSlow,
+                                      CCTK_INT RHSGroupIndexSlow);
 CCTK_INT MoL_RegisterConstrainedGroup(CCTK_INT ConstrainedGroupIndex);
 CCTK_INT MoL_RegisterSaveAndRestoreGroup(CCTK_INT SandRGroupIndex);
 
diff --git a/src/MoLFunctions.h b/src/MoLFunctions.h
index 9f99e73..6d92bf0 100644
--- a/src/MoLFunctions.h
+++ b/src/MoLFunctions.h
@@ -13,6 +13,7 @@
 #define MOL_FUNCTIONS_H
 
 CCTK_INT MoL_RegisterEvolved(CCTK_INT EvolvedIndex, CCTK_INT RHSIndex);
+CCTK_INT MoL_RegisterEvolvedSlow(CCTK_INT EvolvedIndex, CCTK_INT RHSIndexSlow);
 CCTK_INT MoL_RegisterConstrained(CCTK_INT ConstrainedIndex);
 CCTK_INT MoL_RegisterSaveAndRestore(CCTK_INT SandRIndex);
 CCTK_INT MoL_ChangeToEvolved(CCTK_INT EvolvedIndex, CCTK_INT RHSIndex);
@@ -21,6 +22,8 @@ CCTK_INT MoL_ChangeToSaveAndRestore(CCTK_INT SandRIndex);
 CCTK_INT MoL_ChangeToNone(CCTK_INT RemoveIndex);
 CCTK_INT MoL_RegisterEvolvedGroup(CCTK_INT EvolvedGroupIndex, 
                                   CCTK_INT RHSGroupIndex);
+CCTK_INT MoL_RegisterEvolvedGroupSlow(CCTK_INT EvolvedGroupIndexSlow,
+                                      CCTK_INT RHSGroupIndexSlow);
 CCTK_INT MoL_RegisterConstrainedGroup(CCTK_INT ConstrainedGroupIndex);
 CCTK_INT MoL_RegisterSaveAndRestoreGroup(CCTK_INT SandRGroupIndex);
 
diff --git a/src/ParamCheck.c b/src/ParamCheck.c
index 0257eb2..5ae51d6 100644
--- a/src/ParamCheck.c
+++ b/src/ParamCheck.c
@@ -159,6 +159,24 @@ void MoL_ParamCheck(CCTK_ARGUMENTS)
                    "number of intermediate steps must be 13 " 
                    "and the number of scratch levels at least 13.");
   }
+  
+  if ( (CCTK_Equals(ODE_Method, "RK2-MR-2:1"))&&( !((MoL_Intermediate_Steps == 5) || (MoL_Num_Scratch_Levels > 4))) )
+  {
+    CCTK_PARAMWARN("When using the multirate 2-1 RK2 evolver the "
+                   "number of intermediate steps must be 5 and the number of scratch levels at least 5");
+  }
+
+  if ( (CCTK_Equals(ODE_Method, "RK4-MR-2:1"))&&( !((MoL_Intermediate_Steps == 10) || (MoL_Num_Scratch_Levels > 9))) )
+  {
+    CCTK_PARAMWARN("When using the multirate 2-1 RK4 evolver the "
+                   "number of intermediate steps must be 10 and the number of scratch levels at least 10");
+  }
+
+  if ( (CCTK_Equals(ODE_Method, "RK4-RK2"))&&( !((MoL_Intermediate_Steps == 4) || (MoL_Num_Scratch_Levels > 3))) )
+  {
+    CCTK_PARAMWARN("When using the multirate RK4-RK2 evolver the "
+                   "number of intermediate steps must be 4 and the number of scratch levels at least 3");
+  }
 
   if (adaptive_stepsize)
   {
diff --git a/src/RK2-MR-2_1.c b/src/RK2-MR-2_1.c
new file mode 100644
index 0000000..8eaabd2
--- /dev/null
+++ b/src/RK2-MR-2_1.c
@@ -0,0 +1,358 @@
+ /*@@
+   @file      RK2-MR-2_1.c
+   @date      2012-03-25
+   @author    Christian Reisswig
+   @desc 
+   A routine to perform RK2 2:1 MR evolution. Mostly copied from
+   genericRK.c
+   @enddesc 
+   @version   $Header$
+ @@*/
+
+#include "cctk.h"
+#include "cctk_Arguments.h"
+#include "cctk_Parameters.h"
+
+#include <stdio.h>
+#include "ExternalVariables.h"
+
+//#define MOLDEBUG
+
+static const char *rcsid = "$Header$";
+
+CCTK_FILEVERSION(CactusBase_MoL_RK2_MR_2_1_c);
+
+/********************************************************************
+ *********************     Local Data Types   ***********************
+ ********************************************************************/
+
+/********************************************************************
+ ********************* Local Routine Prototypes *********************
+ ********************************************************************/
+
+/********************************************************************
+ ***************** Scheduled Routine Prototypes *********************
+ ********************************************************************/
+
+void MoL_RK2_MR_2_1_Add(CCTK_ARGUMENTS);
+
+/********************************************************************
+ ********************* Other Routine Prototypes *********************
+ ********************************************************************/
+
+/********************************************************************
+ *********************     Local Data   *****************************
+ ********************************************************************/
+
+/********************************************************************
+ *********************     External Routines   **********************
+ ********************************************************************/
+
+ /*@@
+   @routine    MoL_RK2_MR_2_1_Add
+   @date       
+   @author     
+   @desc 
+   Performs a single step of a RK2_MR_2_1 type time
+   integration.
+   @enddesc 
+   @calls     
+   @calledby   
+   @history 
+ 
+   @endhistory 
+
+@@*/
+
+void MoL_RK2_MR_2_1_Add(CCTK_ARGUMENTS)
+{
+
+  DECLARE_CCTK_ARGUMENTS;
+  DECLARE_CCTK_PARAMETERS;
+    
+  CCTK_INT arraydim;
+
+  /* FIXME */
+
+#ifdef MOLDOESCOMPLEX
+
+CCTK_WARN(0, "not implemented");
+  CCTK_COMPLEX Complex_alpha, Complex_beta, Complex_Delta_Time;
+  CCTK_COMPLEX * restrict OldComplexVar;
+  CCTK_COMPLEX * restrict UpdateComplexVar;
+  CCTK_COMPLEX const * restrict RHSComplexVar;
+  CCTK_COMPLEX * restrict ScratchComplexVar;
+
+#endif
+
+  static CCTK_INT scratchspace_firstindex = -99;
+  static CCTK_INT scratchspace_firstindex_slow = -99;
+  CCTK_INT index, var, scratchstep;
+  CCTK_INT totalsize;
+  CCTK_REAL alpha[4], beta[5];
+  CCTK_REAL alpha_slow[4], beta_slow[5];
+  CCTK_REAL * restrict UpdateVar;
+  CCTK_REAL * restrict OldVar;
+  CCTK_REAL const * restrict RHSVar;
+  CCTK_REAL * restrict ScratchVar;
+
+  totalsize = 1;
+  for (arraydim = 0; arraydim < cctk_dim; arraydim++)
+  {
+     totalsize *= cctk_lsh[arraydim];
+  }
+
+  if (scratchspace_firstindex == -99)
+  {
+    scratchspace_firstindex = CCTK_FirstVarIndex("MOL::SCRATCHSPACE");
+  }
+  
+  if (scratchspace_firstindex_slow == -99)
+  {
+    scratchspace_firstindex_slow = CCTK_FirstVarIndex("MOL::SCRATCHSPACESLOW");
+  }
+
+  switch (MoL_Intermediate_Steps - (*MoL_Intermediate_Step))
+  {
+    case 0:
+      alpha[0] = 1.0/2.0; 
+      alpha[1] = 0;
+      alpha[2] = 0;
+      alpha[3] = 0;
+      alpha_slow[0] = 1.0/2.0; 
+      alpha_slow[1] = 0;
+      alpha_slow[2] = 0;
+      alpha_slow[3] = 0;
+      break;
+    case 1:
+      alpha[0] = 1.0/4.0;
+      alpha[1] = 1.0/4.0;
+      alpha[2] = 0;
+      alpha[3] = 0;
+      alpha_slow[0] = 1.0/2.0; 
+      alpha_slow[1] = 0;
+      alpha_slow[2] = 0;
+      alpha_slow[3] = 0;
+      break;
+    case 2:
+      alpha[0] = 1.0/4.0; 
+      alpha[1] = 1.0/4.0;
+      alpha[2] = 1.0/2.0;
+      alpha[3] = 0.0;
+      alpha_slow[0] = 1.0; 
+      alpha_slow[1] = 0;
+      alpha_slow[2] = 0;
+      alpha_slow[3] = 0;
+      break;
+    case 3:
+      alpha[0] = 1.0/4.0;
+      alpha[1] = 1.0/4.0; 
+      alpha[2] = 1.0/4.0; 
+      alpha[3] = 1.0/4.0;
+      alpha_slow[0] = 1.0; 
+      alpha_slow[1] = 0;
+      alpha_slow[2] = 0;
+      alpha_slow[3] = 0;
+  }
+
+  beta[0] = 1.0/4.0; 
+  beta[1] = 1.0/4.0; 
+  beta[2] = 1.0/4.0; 
+  beta[3] = 1.0/4.0;
+  beta[4] = 0.0;
+
+  beta_slow[0] = 1.0/2.0; 
+  beta_slow[1] = 0.0; 
+  beta_slow[2] = 0.0; 
+  beta_slow[3] = 0.0;
+  beta_slow[4] = 1.0/2.0;
+
+
+  /* FIXME */
+
+
+  /* Real GFs, the "fast" part */
+
+  for (var = 0; var < MoLNumEvolvedVariables; var++)
+  {
+    
+    UpdateVar = (CCTK_REAL *)CCTK_VarDataPtrI(cctkGH, 0, 
+                                              EvolvedVariableIndex[var]);
+    OldVar = (CCTK_REAL *)CCTK_VarDataPtrI(cctkGH, 1, 
+                                              EvolvedVariableIndex[var]);
+    RHSVar = (CCTK_REAL const *)CCTK_VarDataPtrI(cctkGH, 0, 
+                                                 RHSVariableIndex[var]);
+/* #define MOLDEBUG 1 */
+#ifdef MOLDEBUG
+    printf("In multirate RK. Variable %d (%s). RHS %d (%s). beta %g.\n",
+           EvolvedVariableIndex[var],
+           CCTK_VarName(EvolvedVariableIndex[var]),
+           RHSVariableIndex[var],
+           CCTK_VarName(RHSVariableIndex[var]),
+           beta[MoL_Intermediate_Steps - (*MoL_Intermediate_Step)]);
+#endif
+
+    ScratchVar = CCTK_VarDataPtrI(cctkGH, 0, 
+                                      scratchspace_firstindex
+                                      + var
+                                      + MoLNumEvolvedVariables * (MoL_Intermediate_Steps - (*MoL_Intermediate_Step)));
+
+#pragma omp parallel for
+    for (index = 0; index < totalsize; index++)
+    {
+      ScratchVar[index] = (*Original_Delta_Time) / cctkGH->cctk_timefac * RHSVar[index];
+    
+#ifdef MOLDEBUG
+      if (CCTK_EQUALS(verbose,"extreme"))
+      {
+        printf("Variable: %d. Index: %d. dt: %f. beta %f. RHS: %f. q: %f.\n",
+               var, index, (*Original_Delta_Time) / cctkGH->cctk_timefac, beta, RHSVar[index], 
+               UpdateVar[index]);
+      }
+#endif
+    }
+    
+    
+#pragma omp parallel for
+    for (index = 0; index < totalsize; index++)
+    {
+      UpdateVar[index] = OldVar[index];
+    }
+
+    if ((*MoL_Intermediate_Step)>1)
+    {
+      //printf("Step %d \n", MoL_Intermediate_Steps - (*MoL_Intermediate_Step));
+      for (scratchstep = 0; scratchstep <= MoL_Intermediate_Steps - (*MoL_Intermediate_Step); scratchstep++)
+      {
+         
+         //printf("Scratch Step %d, alpha %g \n", scratchstep, alpha[scratchstep]);
+         
+         ScratchVar = CCTK_VarDataPtrI(cctkGH, 0, 
+                                       scratchspace_firstindex
+                                       + var
+                                       + MoLNumEvolvedVariables * scratchstep);
+      
+#pragma omp parallel for
+         for (index = 0; index < totalsize; index++)
+         {
+            UpdateVar[index] += alpha[scratchstep] * ScratchVar[index];
+         }
+      }
+    }
+    else
+    {
+      //printf("Final Step!\n");
+    
+      for (scratchstep = 0; scratchstep < MoL_Intermediate_Steps; scratchstep++)
+      {
+         
+         //printf("Scratch Step %d, beta %g \n", scratchstep, beta[scratchstep]);
+         
+         ScratchVar = CCTK_VarDataPtrI(cctkGH, 0, 
+                                       scratchspace_firstindex
+                                       + var
+                                       + MoLNumEvolvedVariables * scratchstep);
+      
+#pragma omp parallel for
+         for (index = 0; index < totalsize; index++)
+         {
+            UpdateVar[index] += beta[scratchstep] * ScratchVar[index];
+         }
+      }
+    }
+
+  }
+
+
+  for (var = 0; var < MoLNumEvolvedVariablesSlow; var++)
+  {
+    
+    UpdateVar = (CCTK_REAL *)CCTK_VarDataPtrI(cctkGH, 0, 
+                                              EvolvedVariableIndexSlow[var]);
+    OldVar = (CCTK_REAL *)CCTK_VarDataPtrI(cctkGH, 1, 
+                                              EvolvedVariableIndexSlow[var]);
+    RHSVar = (CCTK_REAL const *)CCTK_VarDataPtrI(cctkGH, 0, 
+                                                 RHSVariableIndexSlow[var]);
+/* #define MOLDEBUG 1 */
+#ifdef MOLDEBUG
+    printf("In multirate RK. SLOW Variable %d (%s). RHS %d (%s). beta %g.\n",
+           EvolvedVariableIndexSlow[var],
+           CCTK_VarName(EvolvedVariableIndexSlow[var]),
+           RHSVariableIndexSlow[var],
+           CCTK_VarName(RHSVariableIndexSlow[var]),
+           beta[MoL_Intermediate_Steps - (*MoL_Intermediate_Step)]);
+#endif
+
+    ScratchVar = CCTK_VarDataPtrI(cctkGH, 0, 
+                                      scratchspace_firstindex_slow
+                                      + var
+                                      + MoLNumEvolvedVariablesSlow * (MoL_Intermediate_Steps - (*MoL_Intermediate_Step)));
+
+#pragma omp parallel for
+    for (index = 0; index < totalsize; index++)
+    {
+      ScratchVar[index] = (*Original_Delta_Time) / cctkGH->cctk_timefac * RHSVar[index];
+    
+#ifdef MOLDEBUG
+      if (CCTK_EQUALS(verbose,"extreme"))
+      {
+        printf("SLOW Variable: %d. Index: %d. dt: %f. beta %f. RHS: %f. q: %f.\n",
+               var, index, (*Original_Delta_Time) / cctkGH->cctk_timefac, beta, RHSVar[index], 
+               UpdateVar[index]);
+      }
+#endif
+    }
+    
+    
+#pragma omp parallel for
+    for (index = 0; index < totalsize; index++)
+    {
+      UpdateVar[index] = OldVar[index];
+    }
+
+    if ((*MoL_Intermediate_Step)>1)
+    {
+      //printf("Step %d \n", MoL_Intermediate_Steps - (*MoL_Intermediate_Step));
+      for (scratchstep = 0; scratchstep <= /*MoL_Intermediate_Steps - (*MoL_Intermediate_Step)*/ 0; scratchstep++)
+      {
+         
+         //printf("Scratch Step %d, alpha %g \n", scratchstep, alpha_slow[scratchstep]);
+         
+         ScratchVar = CCTK_VarDataPtrI(cctkGH, 0, 
+                                       scratchspace_firstindex_slow
+                                       + var
+                                       + MoLNumEvolvedVariablesSlow * scratchstep);
+      
+#pragma omp parallel for
+         for (index = 0; index < totalsize; index++)
+         {
+            UpdateVar[index] += alpha_slow[scratchstep] * ScratchVar[index];
+         }
+      }
+    }
+    else
+    {
+      //printf("Final Step!\n");
+    
+      for (scratchstep = 0; scratchstep < MoL_Intermediate_Steps; scratchstep+=4)
+      {
+         
+         //printf("Scratch Step %d, beta %g \n", scratchstep, beta_slow[scratchstep]);
+         
+         ScratchVar = CCTK_VarDataPtrI(cctkGH, 0, 
+                                       scratchspace_firstindex_slow
+                                       + var
+                                       + MoLNumEvolvedVariablesSlow * scratchstep);
+      
+#pragma omp parallel for
+         for (index = 0; index < totalsize; index++)
+         {
+            UpdateVar[index] += beta_slow[scratchstep] * ScratchVar[index];
+         }
+      }
+    }
+
+  }
+
+  return;
+}
diff --git a/src/RK4-MR-2_1.c b/src/RK4-MR-2_1.c
new file mode 100644
index 0000000..815a41e
--- /dev/null
+++ b/src/RK4-MR-2_1.c
@@ -0,0 +1,502 @@
+ /*@@
+   @file      RK4-MR-2_1.c
+   @date      2012-03-25
+   @author    Christian Reisswig
+   @desc 
+   A routine to perform 3rd order 2:1 multirate RK evolution. Mostly copied
+   from genericRK.c
+   @enddesc 
+   @version   $Header$
+ @@*/
+
+#include "cctk.h"
+#include "cctk_Arguments.h"
+#include "cctk_Parameters.h"
+
+#include <stdio.h>
+#include "ExternalVariables.h"
+
+//#define MOLDEBUG
+
+static const char *rcsid = "$Header$";
+
+CCTK_FILEVERSION(CactusBase_MoL_RK4_MR_2_1_c);
+
+/********************************************************************
+ *********************     Local Data Types   ***********************
+ ********************************************************************/
+
+/********************************************************************
+ ********************* Local Routine Prototypes *********************
+ ********************************************************************/
+
+/********************************************************************
+ ***************** Scheduled Routine Prototypes *********************
+ ********************************************************************/
+
+void MoL_RK4_MR_2_1_Add(CCTK_ARGUMENTS);
+
+/********************************************************************
+ ********************* Other Routine Prototypes *********************
+ ********************************************************************/
+
+/********************************************************************
+ *********************     Local Data   *****************************
+ ********************************************************************/
+
+/********************************************************************
+ *********************     External Routines   **********************
+ ********************************************************************/
+
+ /*@@
+   @routine    MoL_RK4_MR_2_1_Add
+   @date       
+   @author     
+   @desc 
+   Performs a single step of a RK43 2:1 type time
+   integration.
+   @enddesc 
+   @calls     
+   @calledby   
+   @history 
+ 
+   @endhistory 
+
+@@*/
+
+void MoL_RK4_MR_2_1_Add(CCTK_ARGUMENTS)
+{
+
+  DECLARE_CCTK_ARGUMENTS;
+  DECLARE_CCTK_PARAMETERS;
+    
+  CCTK_INT arraydim;
+
+  /* FIXME */
+
+#ifdef MOLDOESCOMPLEX
+
+CCTK_WARN(0, "not implemented");
+  CCTK_COMPLEX Complex_alpha, Complex_beta, Complex_Delta_Time;
+  CCTK_COMPLEX * restrict OldComplexVar;
+  CCTK_COMPLEX * restrict UpdateComplexVar;
+  CCTK_COMPLEX const * restrict RHSComplexVar;
+  CCTK_COMPLEX * restrict ScratchComplexVar;
+
+#endif
+
+  static CCTK_INT scratchspace_firstindex = -99;
+  static CCTK_INT scratchspace_firstindex_slow = -99;
+  CCTK_INT index, var, scratchstep;
+  CCTK_INT totalsize;
+  CCTK_REAL alpha[9], beta[10];
+  CCTK_REAL alpha_slow[9], beta_slow[10];
+  CCTK_REAL * restrict UpdateVar;
+  CCTK_REAL * restrict OldVar;
+  CCTK_REAL const * restrict RHSVar;
+  CCTK_REAL * restrict ScratchVar;
+
+  totalsize = 1;
+  for (arraydim = 0; arraydim < cctk_dim; arraydim++)
+  {
+     totalsize *= cctk_lsh[arraydim];
+  }
+
+  if (scratchspace_firstindex == -99)
+  {
+    scratchspace_firstindex = CCTK_FirstVarIndex("MOL::SCRATCHSPACE");
+  }
+  
+  if (scratchspace_firstindex_slow == -99)
+  {
+    scratchspace_firstindex_slow = CCTK_FirstVarIndex("MOL::SCRATCHSPACESLOW");
+  }
+
+  switch (MoL_Intermediate_Steps - (*MoL_Intermediate_Step))
+  {
+    case 0:
+      alpha[0] = 1.0/4.0; 
+      alpha[1] = 0;
+      alpha[2] = 0;
+      alpha[3] = 0;
+      alpha[4] = 0;
+      alpha[5] = 0;
+      alpha[6] = 0;
+      alpha[7] = 0;
+      alpha[8] = 0;
+      alpha_slow[0] = 1.0/4.0; 
+      alpha_slow[1] = 0;
+      alpha_slow[2] = 0;
+      alpha_slow[3] = 0;
+      alpha_slow[4] = 0.0; 
+      alpha_slow[5] = 0;
+      alpha_slow[6] = 0;
+      alpha_slow[7] = 0;
+      alpha_slow[8] = 0;
+      break;
+    case 1:
+      alpha[0] = -1.0/12.0;
+      alpha[1] = 1.0/3.0;
+      alpha[2] = 0;
+      alpha[3] = 0;
+      alpha[4] = 0.0;
+      alpha[5] = 0.0;
+      alpha[6] = 0;
+      alpha[7] = 0;
+      alpha[8] = 0;
+      alpha_slow[0] = 1.0/4.0; 
+      alpha_slow[1] = 0;
+      alpha_slow[2] = 0;
+      alpha_slow[3] = 0;
+      alpha_slow[4] = 0.0; 
+      alpha_slow[5] = 0;
+      alpha_slow[6] = 0;
+      alpha_slow[7] = 0;
+      alpha_slow[8] = 0;
+      break;
+    case 2:
+      alpha[0] = 1.0/6.0; 
+      alpha[1] = -1.0/6.0;
+      alpha[2] = 1.0/2.0;
+      alpha[3] = 0.0;
+      alpha[4] = 0.0; 
+      alpha[5] = 0.0;
+      alpha[6] = 0.0;
+      alpha[7] = 0.0;
+      alpha[8] = 0.0;
+      alpha_slow[0] = 1.0/2.0; 
+      alpha_slow[1] = 0;
+      alpha_slow[2] = 0;
+      alpha_slow[3] = 0;
+      alpha_slow[4] = 0.0; 
+      alpha_slow[5] = 0;
+      alpha_slow[6] = 0;
+      alpha_slow[7] = 0;
+      alpha_slow[8] = 0;
+      break;
+    case 3:
+      alpha[0] = 1.0/12.0;
+      alpha[1] = 1.0/6.0; 
+      alpha[2] = 1.0/6.0; 
+      alpha[3] = 1.0/12.0;
+      alpha[4] = 0.0; 
+      alpha[5] = 0.0;
+      alpha[6] = 0.0;
+      alpha[7] = 0.0;
+      alpha[8] = 0.0;
+      alpha_slow[0] = 1.0/2.0; 
+      alpha_slow[1] = 0;
+      alpha_slow[2] = 0;
+      alpha_slow[3] = 0;
+      alpha_slow[4] = 0.0; 
+      alpha_slow[5] = 0;
+      alpha_slow[6] = 0;
+      alpha_slow[7] = 0;
+      alpha_slow[8] = 0;
+    case 4:
+      alpha[0] = 1.0/12.0;
+      alpha[1] = 1.0/6.0; 
+      alpha[2] = 1.0/6.0; 
+      alpha[3] = 1.0/6.0;
+      alpha[4] = 1.0/12.0; 
+      alpha[5] = 0.0;
+      alpha[6] = 0.0;
+      alpha[7] = 0.0;
+      alpha[8] = 0.0;
+      alpha_slow[0] = -1.0/6.0; 
+      alpha_slow[1] = 0;
+      alpha_slow[2] = 0;
+      alpha_slow[3] = 0;
+      alpha_slow[4] = 2.0/3.0; 
+      alpha_slow[5] = 0;
+      alpha_slow[6] = 0;
+      alpha_slow[7] = 0;
+      alpha_slow[7] = 0;
+    case 5:
+      alpha[0] = 1.0/12.0;
+      alpha[1] = 1.0/6.0; 
+      alpha[2] = 1.0/6.0; 
+      alpha[3] = 1.0/12.0;
+      alpha[4] = 0.0; 
+      alpha[5] = 1.0/4.0;
+      alpha[6] = 0.0;
+      alpha[7] = 0.0;
+      alpha[8] = 0.0;
+      alpha_slow[0] = 1.0/12.0;
+      alpha_slow[1] = 0;
+      alpha_slow[2] = 0;
+      alpha_slow[3] = 0;
+      alpha_slow[4] = 1.0/6.0; 
+      alpha_slow[5] = 1.0/2.0;
+      alpha_slow[6] = 0;
+      alpha_slow[7] = 0;
+      alpha_slow[7] = 0;
+    case 6:
+      alpha[0] = 1.0/12.0;
+      alpha[1] = 1.0/6.0; 
+      alpha[2] = 1.0/6.0; 
+      alpha[3] = 1.0/12.0;
+      alpha[4] = 0.0; 
+      alpha[5] = -1.0/12.0;
+      alpha[6] = 1.0/3.0;
+      alpha[7] = 0.0;
+      alpha[8] = 0.0;
+      alpha_slow[0] = 1.0/12.0; 
+      alpha_slow[1] = 0;
+      alpha_slow[2] = 0;
+      alpha_slow[3] = 0;
+      alpha_slow[4] = 1.0/6.0; 
+      alpha_slow[5] = 1.0/2.0;
+      alpha_slow[6] = 0;
+      alpha_slow[7] = 0;
+      alpha_slow[7] = 0;
+    case 7:
+      alpha[0] = 1.0/12.0;
+      alpha[1] = 1.0/6.0; 
+      alpha[2] = 1.0/6.0; 
+      alpha[3] = 1.0/12.0;
+      alpha[4] = 0.0; 
+      alpha[5] = 1.0/6.0;
+      alpha[6] = -1.0/6.0;
+      alpha[7] = 1.0/2.0;
+      alpha[8] = 0.0;
+      alpha_slow[0] = 1.0/3.0; 
+      alpha_slow[1] = 0;
+      alpha_slow[2] = 0;
+      alpha_slow[3] = 0;
+      alpha_slow[4] = -1.0/3.0; 
+      alpha_slow[5] = 1.0;
+      alpha_slow[6] = 0;
+      alpha_slow[7] = 0;
+      alpha_slow[7] = 0;
+    case 8:
+      alpha[0] = 1.0/12.0;
+      alpha[1] = 1.0/6.0; 
+      alpha[2] = 1.0/6.0; 
+      alpha[3] = 1.0/12.0;
+      alpha[4] = 0.0; 
+      alpha[5] = 1.0/12.0;
+      alpha[6] = 1.0/6.0;
+      alpha[7] = 1.0/6.0;
+      alpha[8] = 1.0/12.0;
+      alpha_slow[0] = 1.0/3.0; 
+      alpha_slow[1] = 0;
+      alpha_slow[2] = 0;
+      alpha_slow[3] = 0;
+      alpha_slow[4] = -1.0/3.0; 
+      alpha_slow[5] = 1.0;
+      alpha_slow[6] = 0;
+      alpha_slow[7] = 0;
+      alpha_slow[7] = 0;
+  }
+
+  beta[0] = 1.0/12.0;
+  beta[1] = 1.0/6.0; 
+  beta[2] = 1.0/6.0; 
+  beta[3] = 1.0/12.0;
+  beta[4] = 0.0;
+  beta[5] = 1.0/12.0;
+  beta[6] = 1.0/6.0; 
+  beta[7] = 1.0/6.0; 
+  beta[8] = 1.0/12.0;
+  beta[9] = 0.0;
+
+  beta_slow[0] = 1.0/6.0; 
+  beta_slow[1] = 0.0; 
+  beta_slow[2] = 0.0; 
+  beta_slow[3] = 0.0;
+  beta_slow[4] = 1.0/3.0;
+  beta_slow[5] = 1.0/3.0; 
+  beta_slow[6] = 0.0; 
+  beta_slow[7] = 0.0; 
+  beta_slow[8] = 0.0;
+  beta_slow[9] = 1.0/6.0;
+
+  /* FIXME */
+
+
+  /* Real GFs, the "fast" part */
+
+  for (var = 0; var < MoLNumEvolvedVariables; var++)
+  {
+    
+    UpdateVar = (CCTK_REAL *)CCTK_VarDataPtrI(cctkGH, 0, 
+                                              EvolvedVariableIndex[var]);
+    OldVar = (CCTK_REAL *)CCTK_VarDataPtrI(cctkGH, 1, 
+                                              EvolvedVariableIndex[var]);
+    RHSVar = (CCTK_REAL const *)CCTK_VarDataPtrI(cctkGH, 0, 
+                                                 RHSVariableIndex[var]);
+/* #define MOLDEBUG 1 */
+#ifdef MOLDEBUG
+    printf("In multirate RK. Variable %d (%s). RHS %d (%s). beta %g.\n",
+           EvolvedVariableIndex[var],
+           CCTK_VarName(EvolvedVariableIndex[var]),
+           RHSVariableIndex[var],
+           CCTK_VarName(RHSVariableIndex[var]),
+           beta[MoL_Intermediate_Steps - (*MoL_Intermediate_Step)]);
+#endif
+
+    ScratchVar = CCTK_VarDataPtrI(cctkGH, 0, 
+                                      scratchspace_firstindex
+                                      + var
+                                      + MoLNumEvolvedVariables * (MoL_Intermediate_Steps - (*MoL_Intermediate_Step)));
+
+#pragma omp parallel for
+    for (index = 0; index < totalsize; index++)
+    {
+      ScratchVar[index] = (*Original_Delta_Time) / cctkGH->cctk_timefac * RHSVar[index];
+    
+#ifdef MOLDEBUG
+      if (CCTK_EQUALS(verbose,"extreme"))
+      {
+        printf("Variable: %d. Index: %d. dt: %f. beta %f. RHS: %f. q: %f.\n",
+               var, index, (*Original_Delta_Time) / cctkGH->cctk_timefac, beta, RHSVar[index], 
+               UpdateVar[index]);
+      }
+#endif
+    }
+    
+    
+#pragma omp parallel for
+    for (index = 0; index < totalsize; index++)
+    {
+      UpdateVar[index] = OldVar[index];
+    }
+
+    if ((*MoL_Intermediate_Step)>1)
+    {
+      //printf("Step %d \n", MoL_Intermediate_Steps - (*MoL_Intermediate_Step));
+      for (scratchstep = 0; scratchstep <= MoL_Intermediate_Steps - (*MoL_Intermediate_Step); scratchstep++)
+      {
+         
+         //printf("Scratch Step %d, alpha %g \n", scratchstep, alpha[scratchstep]);
+         
+         ScratchVar = CCTK_VarDataPtrI(cctkGH, 0, 
+                                       scratchspace_firstindex
+                                       + var
+                                       + MoLNumEvolvedVariables * scratchstep);
+      
+#pragma omp parallel for
+         for (index = 0; index < totalsize; index++)
+         {
+            UpdateVar[index] += alpha[scratchstep] * ScratchVar[index];
+         }
+      }
+    }
+    else
+    {
+      //printf("Final Step!\n");
+    
+      for (scratchstep = 0; scratchstep < MoL_Intermediate_Steps; scratchstep++)
+      {
+         
+         //printf("Scratch Step %d, beta %g \n", scratchstep, beta[scratchstep]);
+         
+         ScratchVar = CCTK_VarDataPtrI(cctkGH, 0, 
+                                       scratchspace_firstindex
+                                       + var
+                                       + MoLNumEvolvedVariables * scratchstep);
+      
+#pragma omp parallel for
+         for (index = 0; index < totalsize; index++)
+         {
+            UpdateVar[index] += beta[scratchstep] * ScratchVar[index];
+         }
+      }
+    }
+
+  }
+
+
+  for (var = 0; var < MoLNumEvolvedVariablesSlow; var++)
+  {
+    
+    UpdateVar = (CCTK_REAL *)CCTK_VarDataPtrI(cctkGH, 0, 
+                                              EvolvedVariableIndexSlow[var]);
+    OldVar = (CCTK_REAL *)CCTK_VarDataPtrI(cctkGH, 1, 
+                                              EvolvedVariableIndexSlow[var]);
+    RHSVar = (CCTK_REAL const *)CCTK_VarDataPtrI(cctkGH, 0, 
+                                                 RHSVariableIndexSlow[var]);
+/* #define MOLDEBUG 1 */
+#ifdef MOLDEBUG
+    printf("In multirate RK. SLOW Variable %d (%s). RHS %d (%s). beta %g.\n",
+           EvolvedVariableIndexSlow[var],
+           CCTK_VarName(EvolvedVariableIndexSlow[var]),
+           RHSVariableIndexSlow[var],
+           CCTK_VarName(RHSVariableIndexSlow[var]),
+           beta[MoL_Intermediate_Steps - (*MoL_Intermediate_Step)]);
+#endif
+
+    ScratchVar = CCTK_VarDataPtrI(cctkGH, 0, 
+                                      scratchspace_firstindex_slow
+                                      + var
+                                      + MoLNumEvolvedVariablesSlow * (MoL_Intermediate_Steps - (*MoL_Intermediate_Step)));
+
+#pragma omp parallel for
+    for (index = 0; index < totalsize; index++)
+    {
+      ScratchVar[index] = (*Original_Delta_Time) / cctkGH->cctk_timefac * RHSVar[index];
+    
+#ifdef MOLDEBUG
+      if (CCTK_EQUALS(verbose,"extreme"))
+      {
+        printf("SLOW Variable: %d. Index: %d. dt: %f. beta %f. RHS: %f. q: %f.\n",
+               var, index, (*Original_Delta_Time) / cctkGH->cctk_timefac, beta, RHSVar[index], 
+               UpdateVar[index]);
+      }
+#endif
+    }
+    
+    
+#pragma omp parallel for
+    for (index = 0; index < totalsize; index++)
+    {
+      UpdateVar[index] = OldVar[index];
+    }
+
+    if ((*MoL_Intermediate_Step)>1)
+    {
+      //printf("Step %d \n", MoL_Intermediate_Steps - (*MoL_Intermediate_Step));
+      for (scratchstep = 0; scratchstep <= MoL_Intermediate_Steps - (*MoL_Intermediate_Step); scratchstep++)
+      {
+         
+         //printf("Scratch Step %d, alpha %g \n", scratchstep, alpha[scratchstep]);
+         
+         ScratchVar = CCTK_VarDataPtrI(cctkGH, 0, 
+                                       scratchspace_firstindex_slow
+                                       + var
+                                       + MoLNumEvolvedVariablesSlow * scratchstep);
+      
+#pragma omp parallel for
+         for (index = 0; index < totalsize; index++)
+         {
+            UpdateVar[index] += alpha_slow[scratchstep] * ScratchVar[index];
+         }
+      }
+    }
+    else
+    {
+      //printf("Final Step!\n");
+    
+      for (scratchstep = 0; scratchstep < MoL_Intermediate_Steps; scratchstep++)
+      {
+         
+         //printf("Scratch Step %d, beta %g \n", scratchstep, beta[scratchstep]);
+         
+         ScratchVar = CCTK_VarDataPtrI(cctkGH, 0, 
+                                       scratchspace_firstindex_slow
+                                       + var
+                                       + MoLNumEvolvedVariablesSlow * scratchstep);
+      
+#pragma omp parallel for
+         for (index = 0; index < totalsize; index++)
+         {
+            UpdateVar[index] += beta_slow[scratchstep] * ScratchVar[index];
+         }
+      }
+    }
+
+  }
+
+  return;
+}
diff --git a/src/RK4-RK2.c b/src/RK4-RK2.c
new file mode 100644
index 0000000..40e84fb
--- /dev/null
+++ b/src/RK4-RK2.c
@@ -0,0 +1,342 @@
+ /*@@
+   @file      RK4-RK2.c
+   @date      2012-03-25
+   @author    Christian Reisswig
+   @desc 
+   A routine to perform homegrown RK4RK2 evolution. Mostly copied from
+   genericRK.c
+   @enddesc 
+   @version   $Header$
+ @@*/
+
+#include "cctk.h"
+#include "cctk_Arguments.h"
+#include "cctk_Parameters.h"
+
+#include <stdio.h>
+#include "ExternalVariables.h"
+
+//#define MOLDEBUG
+
+static const char *rcsid = "$Header$";
+
+CCTK_FILEVERSION(CactusBase_MoL_RK4_RK2_c);
+
+/********************************************************************
+ *********************     Local Data Types   ***********************
+ ********************************************************************/
+
+/********************************************************************
+ ********************* Local Routine Prototypes *********************
+ ********************************************************************/
+
+/********************************************************************
+ ***************** Scheduled Routine Prototypes *********************
+ ********************************************************************/
+
+void MoL_RK4_RK2_Add(CCTK_ARGUMENTS);
+
+/********************************************************************
+ ********************* Other Routine Prototypes *********************
+ ********************************************************************/
+
+/********************************************************************
+ *********************     Local Data   *****************************
+ ********************************************************************/
+
+/********************************************************************
+ *********************     External Routines   **********************
+ ********************************************************************/
+
+ /*@@
+   @routine    MoL_RK4_RK2_Add
+   @date       
+   @author     
+   @desc 
+   Performs a single step of a RK4_RK2 type time
+   integration.
+   @enddesc 
+   @calls     
+   @calledby   
+   @history 
+ 
+   @endhistory 
+
+@@*/
+
+void MoL_RK4_RK2_Add(CCTK_ARGUMENTS)
+{
+
+  DECLARE_CCTK_ARGUMENTS;
+  DECLARE_CCTK_PARAMETERS;
+    
+  CCTK_INT arraydim;
+
+  /* FIXME */
+
+#ifdef MOLDOESCOMPLEX
+
+CCTK_WARN(0, "not implemented");
+  CCTK_COMPLEX Complex_alpha, Complex_beta, Complex_Delta_Time;
+  CCTK_COMPLEX * restrict OldComplexVar;
+  CCTK_COMPLEX * restrict UpdateComplexVar;
+  CCTK_COMPLEX const * restrict RHSComplexVar;
+  CCTK_COMPLEX * restrict ScratchComplexVar;
+
+#endif
+
+  static CCTK_INT scratchspace_firstindex = -99;
+  static CCTK_INT scratchspace_firstindex_slow = -99;
+  CCTK_INT index, var, scratchstep;
+  CCTK_INT totalsize;
+  CCTK_REAL alpha[3], beta[4];
+  CCTK_REAL alpha_slow[3], beta_slow[4];
+  CCTK_REAL * restrict UpdateVar;
+  CCTK_REAL * restrict OldVar;
+  CCTK_REAL const * restrict RHSVar;
+  CCTK_REAL * restrict ScratchVar;
+
+  totalsize = 1;
+  for (arraydim = 0; arraydim < cctk_dim; arraydim++)
+  {
+     totalsize *= cctk_lsh[arraydim];
+  }
+
+  if (scratchspace_firstindex == -99)
+  {
+    scratchspace_firstindex = CCTK_FirstVarIndex("MOL::SCRATCHSPACE");
+  }
+  
+  if (scratchspace_firstindex_slow == -99)
+  {
+    scratchspace_firstindex_slow = CCTK_FirstVarIndex("MOL::SCRATCHSPACESLOW");
+  }
+
+  switch (MoL_Intermediate_Steps - (*MoL_Intermediate_Step))
+  {
+    case 0:
+      alpha[0] = 1.0/2.0;
+      alpha[1] = 0;
+      alpha[2] = 0;
+      alpha_slow[0] = 0;
+      alpha_slow[1] = 0;
+      alpha_slow[2] = 0;
+      break;
+    case 1:
+      alpha[0] = 0.0;
+      alpha[1] = 1.0/2.0;
+      alpha[2] = 0;
+      alpha_slow[0] = 0;
+      alpha_slow[1] = 0;
+      alpha_slow[2] = 0;
+      break;
+    case 2:
+      alpha[0] = 0.0;
+      alpha[1] = 0.0;
+      alpha[2] = 1.0;
+      alpha_slow[0] = 1.0;
+      alpha_slow[1] = 0;
+      alpha_slow[2] = 0;
+      break;
+  }
+
+  beta[0] = 1.0/3.0; 
+  beta[1] = 1.0/6.0; 
+  beta[2] = 1.0/6.0; 
+  beta[3] = 1.0/3.0;
+
+  beta_slow[0] = 1.0/2.0; 
+  beta_slow[1] = 0.0; 
+  beta_slow[2] = 0.0; 
+  beta_slow[3] = 1.0/2.0;
+
+
+  /* FIXME */
+
+
+  /* Real GFs, the "fast" part */
+
+  for (var = 0; var < MoLNumEvolvedVariables; var++)
+  {
+    
+    UpdateVar = (CCTK_REAL *)CCTK_VarDataPtrI(cctkGH, 0, 
+                                              EvolvedVariableIndex[var]);
+    OldVar = (CCTK_REAL *)CCTK_VarDataPtrI(cctkGH, 1, 
+                                              EvolvedVariableIndex[var]);
+    RHSVar = (CCTK_REAL const *)CCTK_VarDataPtrI(cctkGH, 0, 
+                                                 RHSVariableIndex[var]);
+/* #define MOLDEBUG 1 */
+#ifdef MOLDEBUG
+    printf("In multirate RK. Variable %d (%s). RHS %d (%s). beta %g.\n",
+           EvolvedVariableIndex[var],
+           CCTK_VarName(EvolvedVariableIndex[var]),
+           RHSVariableIndex[var],
+           CCTK_VarName(RHSVariableIndex[var]),
+           beta[MoL_Intermediate_Steps - (*MoL_Intermediate_Step)]);
+#endif
+
+    ScratchVar = CCTK_VarDataPtrI(cctkGH, 0, 
+                                      scratchspace_firstindex
+                                      + var
+                                      + MoLNumEvolvedVariables * (MoL_Intermediate_Steps - (*MoL_Intermediate_Step)));
+
+#pragma omp parallel for
+    for (index = 0; index < totalsize; index++)
+    {
+      ScratchVar[index] = (*Original_Delta_Time) / cctkGH->cctk_timefac * RHSVar[index];
+    
+#ifdef MOLDEBUG
+      if (CCTK_EQUALS(verbose,"extreme"))
+      {
+        printf("Variable: %d. Index: %d. dt: %f. beta %f. RHS: %f. q: %f.\n",
+               var, index, (*Original_Delta_Time) / cctkGH->cctk_timefac, beta, RHSVar[index], 
+               UpdateVar[index]);
+      }
+#endif
+    }
+    
+    
+#pragma omp parallel for
+    for (index = 0; index < totalsize; index++)
+    {
+      UpdateVar[index] = OldVar[index];
+    }
+    
+
+    if ((*MoL_Intermediate_Step)>1)
+    {
+      //printf("Step %d \n", MoL_Intermediate_Steps - (*MoL_Intermediate_Step));
+      for (scratchstep = 0; scratchstep <= MoL_Intermediate_Steps - (*MoL_Intermediate_Step); scratchstep++)
+      {
+         
+         //printf("Scratch Step %d, alpha %g \n", scratchstep, alpha[scratchstep]);
+         
+         ScratchVar = CCTK_VarDataPtrI(cctkGH, 0, 
+                                       scratchspace_firstindex
+                                       + var
+                                       + MoLNumEvolvedVariables * scratchstep);
+      
+#pragma omp parallel for
+         for (index = 0; index < totalsize; index++)
+         {
+            UpdateVar[index] += alpha[scratchstep] * ScratchVar[index];
+         }
+      }
+    }
+    else
+    {
+      //printf("Final Step!\n");
+    
+      for (scratchstep = 0; scratchstep < MoL_Intermediate_Steps; scratchstep++)
+      {
+         
+         //printf("Scratch Step %d, beta %g \n", scratchstep, beta[scratchstep]);
+         
+         ScratchVar = CCTK_VarDataPtrI(cctkGH, 0, 
+                                       scratchspace_firstindex
+                                       + var
+                                       + MoLNumEvolvedVariables * scratchstep);
+      
+#pragma omp parallel for
+         for (index = 0; index < totalsize; index++)
+         {
+            UpdateVar[index] += beta[scratchstep] * ScratchVar[index];
+         }
+      }
+    }
+
+  }
+
+
+  for (var = 0; var < MoLNumEvolvedVariablesSlow; var++)
+  {
+    
+    UpdateVar = (CCTK_REAL *)CCTK_VarDataPtrI(cctkGH, 0, 
+                                              EvolvedVariableIndexSlow[var]);
+    OldVar = (CCTK_REAL *)CCTK_VarDataPtrI(cctkGH, 1, 
+                                              EvolvedVariableIndexSlow[var]);
+    RHSVar = (CCTK_REAL const *)CCTK_VarDataPtrI(cctkGH, 0, 
+                                                 RHSVariableIndexSlow[var]);
+/* #define MOLDEBUG 1 */
+#ifdef MOLDEBUG
+    printf("In multirate RK. SLOW Variable %d (%s). RHS %d (%s). beta %g.\n",
+           EvolvedVariableIndexSlow[var],
+           CCTK_VarName(EvolvedVariableIndexSlow[var]),
+           RHSVariableIndexSlow[var],
+           CCTK_VarName(RHSVariableIndexSlow[var]),
+           beta[MoL_Intermediate_Steps - (*MoL_Intermediate_Step)]);
+#endif
+
+    ScratchVar = CCTK_VarDataPtrI(cctkGH, 0, 
+                                      scratchspace_firstindex_slow
+                                      + var
+                                      + MoLNumEvolvedVariablesSlow * (MoL_Intermediate_Steps - (*MoL_Intermediate_Step)));
+
+#pragma omp parallel for
+    for (index = 0; index < totalsize; index++)
+    {
+      ScratchVar[index] = (*Original_Delta_Time) / cctkGH->cctk_timefac * RHSVar[index];
+    
+#ifdef MOLDEBUG
+      if (CCTK_EQUALS(verbose,"extreme"))
+      {
+        printf("SLOW Variable: %d. Index: %d. dt: %f. beta %f. RHS: %f. q: %f.\n",
+               var, index, (*Original_Delta_Time) / cctkGH->cctk_timefac, beta, RHSVar[index], 
+               UpdateVar[index]);
+      }
+#endif
+    }
+    
+    
+#pragma omp parallel for
+    for (index = 0; index < totalsize; index++)
+    {
+      UpdateVar[index] = OldVar[index];
+    }
+
+    if ((*MoL_Intermediate_Step)>1)
+    {
+      //printf("Step %d \n", MoL_Intermediate_Steps - (*MoL_Intermediate_Step));
+      for (scratchstep = 0; scratchstep <= /*MoL_Intermediate_Steps - (*MoL_Intermediate_Step)*/ 0; scratchstep++)
+      {
+         
+         //printf("Scratch Step %d, alpha %g \n", scratchstep, alpha_slow[scratchstep]);
+         
+         ScratchVar = CCTK_VarDataPtrI(cctkGH, 0, 
+                                       scratchspace_firstindex_slow
+                                       + var
+                                       + MoLNumEvolvedVariablesSlow * scratchstep);
+      
+#pragma omp parallel for
+         for (index = 0; index < totalsize; index++)
+         {
+            UpdateVar[index] += alpha_slow[scratchstep] * ScratchVar[index];
+         }
+      }
+    }
+    else
+    {
+      //printf("Final Step!\n");
+    
+      for (scratchstep = 0; scratchstep < MoL_Intermediate_Steps; scratchstep++)
+      {
+         
+         //printf("Scratch Step %d, beta %g \n", scratchstep, beta_slow[scratchstep]);
+         
+         ScratchVar = CCTK_VarDataPtrI(cctkGH, 0, 
+                                       scratchspace_firstindex_slow
+                                       + var
+                                       + MoLNumEvolvedVariablesSlow * scratchstep);
+      
+#pragma omp parallel for
+         for (index = 0; index < totalsize; index++)
+         {
+            UpdateVar[index] += beta_slow[scratchstep] * ScratchVar[index];
+         }
+      }
+    }
+
+  }
+
+  return;
+}
diff --git a/src/Registration.c b/src/Registration.c
index 06cf953..5f58d36 100644
--- a/src/Registration.c
+++ b/src/Registration.c
@@ -40,6 +40,8 @@ void MoL_SetScheduleStatus(CCTK_ARGUMENTS);
 
 CCTK_INT MoL_RegisterEvolved(CCTK_INT EvolvedIndex, CCTK_INT RHSIndex);
 
+CCTK_INT MoL_RegisterEvolvedSlow(CCTK_INT EvolvedIndex, CCTK_INT RHSIndexSlow);
+
 CCTK_INT MoL_RegisterConstrained(CCTK_INT ConstrainedIndex);
 
 CCTK_INT MoL_RegisterSaveAndRestore(CCTK_INT SandRIndex);
@@ -52,9 +54,10 @@ CCTK_INT MoL_RegisterConstrainedGroup(CCTK_INT ConstrainedGroupIndex);
 CCTK_INT MoL_RegisterSaveAndRestoreGroup(CCTK_INT SandRGroupIndex);
 
 
-
 CCTK_INT MoL_RegisterEvolvedReal(CCTK_INT EvolvedIndex, CCTK_INT RHSIndex);
 
+CCTK_INT MoL_RegisterEvolvedRealSlow(CCTK_INT EvolvedIndex, CCTK_INT RHSIndexSlow);
+
 CCTK_INT MoL_RegisterConstrainedReal(CCTK_INT ConstrainedIndex);
 
 CCTK_INT MoL_RegisterSaveAndRestoreReal(CCTK_INT SandRIndex);
@@ -62,6 +65,10 @@ CCTK_INT MoL_RegisterSaveAndRestoreReal(CCTK_INT SandRIndex);
 CCTK_INT MoL_RegisterEvolvedRealGroup(CCTK_INT EvolvedGroupIndex, 
                                   CCTK_INT RHSGroupIndex);
 
+CCTK_INT MoL_RegisterEvolvedRealGroupSlow(CCTK_INT EvolvedGroupIndex, 
+                                          CCTK_INT RHSGroupIndexSlow);
+
+
 CCTK_INT MoL_RegisterConstrainedRealGroup(CCTK_INT ConstrainedGroupIndex);
 
 CCTK_INT MoL_RegisterSaveAndRestoreRealGroup(CCTK_INT SandRGroupIndex);
@@ -269,6 +276,109 @@ CCTK_INT MoL_RegisterEvolved(CCTK_INT EvolvedIndex, CCTK_INT RHSIndex)
 
 }
 
+
+ /*@@
+   @routine    MoL_RegisterEvolvedSlow
+   @date       Thu May 30 11:36:59 2002
+   @author     Ian Hawke
+   @desc 
+   Given the index of the GF to be evolved and the RHS GF, it stores 
+   the indexes for later use together with various error checking.
+   @enddesc 
+   @calls     
+   @calledby   
+   @history 
+ 
+   @endhistory 
+
+@@*/
+
+CCTK_INT MoL_RegisterEvolvedSlow(CCTK_INT EvolvedIndex, CCTK_INT RHSIndexSlow)
+{
+
+  CCTK_INT retval, GroupType;
+  
+  retval = 0;
+
+  if (!ScheduleStatus)
+  {
+    CCTK_VWarn(0, __LINE__, __FILE__, CCTK_THORNSTRING,
+               "MoL registration routine called too early!\n"
+               "Trying to register variable '%s',"
+               "Please ensure that all calls to MoL registration routines "
+               "occur within the \"MoL_Register\" timebin.",
+               CCTK_VarName(EvolvedIndex));
+    retval++;
+  }
+
+  GroupType = CCTK_GroupTypeFromVarI(EvolvedIndex);
+  if (GroupType < 0)
+  {
+    CCTK_VWarn(0, __LINE__, __FILE__, CCTK_THORNSTRING,
+               "Evolved index %i is not a real variable index.",
+               EvolvedIndex);
+    retval++;
+  }
+
+  if (!retval)
+  {
+    switch (CCTK_VarTypeI(EvolvedIndex))
+    {
+      case CCTK_VARIABLE_REAL:
+        {
+          switch (GroupType)
+          {
+            case CCTK_GF:
+              {
+                retval +=
+                  MoL_RegisterEvolvedRealSlow(EvolvedIndex,
+                                              RHSIndexSlow);
+                break;
+              }
+            case CCTK_ARRAY:
+            case CCTK_SCALAR:
+              {
+                CCTK_VWarn(0,__LINE__,__FILE__,CCTK_THORNSTRING,
+                           "The variable '%s' is not supported for multirate methods", 
+                           CCTK_VarName(EvolvedIndex));
+                retval++;
+                break;
+              }
+            default:
+              {
+                CCTK_VWarn(0,__LINE__,__FILE__,CCTK_THORNSTRING,
+                           "The variable '%s' is not supported for multirate methods", 
+                           CCTK_VarName(EvolvedIndex));
+                retval++;
+                break;
+              }
+          }
+          break;
+        }
+      case CCTK_VARIABLE_COMPLEX:
+        {
+          CCTK_VWarn(0,__LINE__,__FILE__,CCTK_THORNSTRING,
+                           "The COMPLEX variable '%s' is not supported for multirate methods", 
+                           CCTK_VarName(EvolvedIndex));
+          break;
+        }
+      default:
+        {
+          CCTK_VWarn(0,__LINE__,__FILE__,CCTK_THORNSTRING,
+                     "The variable '%s' is neither REAL nor COMPLEX.", 
+                     CCTK_VarName(EvolvedIndex));
+          retval++;
+          break;
+        }
+    }  
+  }
+  
+  return retval;
+
+}
+
+
+
  /*@@
    @routine    MoL_RegisterConstrained
    @date       Thu May 30 12:35:58 2002
@@ -606,6 +716,91 @@ CCTK_INT MoL_RegisterEvolvedGroup(CCTK_INT EvolvedGroupIndex,
 }
 
 
+CCTK_INT MoL_RegisterEvolvedGroupSlow(CCTK_INT EvolvedGroupIndex, 
+                                      CCTK_INT RHSGroupIndex2)
+{
+
+  CCTK_INT retval, GroupFirstVar;
+  
+  retval = 0;
+
+  if (!ScheduleStatus)
+  {
+    CCTK_VWarn(0, __LINE__, __FILE__, CCTK_THORNSTRING,
+               "MoL registration routine called too early!\n"
+               "Trying to register group '%s',"
+               "Please ensure that all calls to MoL registration routines "
+               "occur within the \"MoL_Register\" timebin.",
+               CCTK_GroupName(EvolvedGroupIndex));
+    retval++;
+  }
+
+  GroupFirstVar = CCTK_FirstVarIndexI(EvolvedGroupIndex);
+  if (GroupFirstVar < 0)
+  {
+    CCTK_VWarn(0, __LINE__, __FILE__, CCTK_THORNSTRING,
+               "Evolved group index %i is not a real group index.",
+               EvolvedGroupIndex);
+    retval++;
+  }
+  
+  switch (CCTK_VarTypeI(CCTK_FirstVarIndexI(EvolvedGroupIndex)))
+  {
+    case CCTK_VARIABLE_REAL:
+      {
+        switch (CCTK_GroupTypeI(EvolvedGroupIndex))
+        {
+          case CCTK_GF:
+            {
+              retval +=
+                MoL_RegisterEvolvedRealGroupSlow(EvolvedGroupIndex, 
+                                                 RHSGroupIndex2);
+              break;
+            }
+          case CCTK_ARRAY:
+          case CCTK_SCALAR:
+            {
+              CCTK_VWarn(0,__LINE__,__FILE__,CCTK_THORNSTRING,
+                         "The group '%s' is not supported by multirate RK", 
+                         CCTK_GroupName(EvolvedGroupIndex));
+              retval++;
+              break;
+            }
+          default:
+            {
+              CCTK_VWarn(0,__LINE__,__FILE__,CCTK_THORNSTRING,
+                         "The group '%s' is not supported by multirate RK", 
+                         CCTK_GroupName(EvolvedGroupIndex));
+              retval++;
+              break;
+            }
+        }
+        break;
+      }
+    case CCTK_VARIABLE_COMPLEX:
+      {
+        CCTK_VWarn(0,__LINE__,__FILE__,CCTK_THORNSTRING,
+                         "The COMPLEX group '%s' is not supported by multirate RK", 
+                         CCTK_GroupName(EvolvedGroupIndex));
+              retval++;
+        break;
+      }
+    default:
+      {
+        CCTK_VWarn(0,__LINE__,__FILE__,CCTK_THORNSTRING,
+                   "The group '%s' is neither REAL nor COMPLEX.", 
+                   CCTK_GroupName(EvolvedGroupIndex));
+        retval++;
+        break;
+      }
+  }  
+  
+  return retval;
+
+}
+
+
+
 CCTK_INT MoL_RegisterConstrainedGroup(CCTK_INT ConstrainedGroupIndex)
 {
 
@@ -992,6 +1187,198 @@ CCTK_INT MoL_RegisterEvolvedReal(CCTK_INT EvolvedIndex, CCTK_INT RHSIndex)
 
 }
 
+
+ /*@@
+   @routine    MoL_RegisterEvolvedSlow
+   @date       Thu May 30 11:36:59 2002
+   @author     Ian Hawke
+   @desc 
+   Given the index of the GF to be evolved and the RHS GF, it stores 
+   the indexes for later use together with various error checking.
+   @enddesc 
+   @calls     
+   @calledby   
+   @history 
+ 
+   @endhistory 
+
+@@*/
+
+CCTK_INT MoL_RegisterEvolvedRealSlow(CCTK_INT EvolvedIndex, CCTK_INT RHSIndexSlow)
+{
+
+  DECLARE_CCTK_PARAMETERS;
+
+  CCTK_INT /*  ierr, */ index, varused, numtimelevs1, numtimelevs2;
+  
+#ifdef MOLDEBUG
+  printf("Arrived in MoLRegisterEvolvedSlow \n");
+  printf("The indexes are %d and %d.\n",EvolvedIndex, RHSIndexSlow);
+  printf("These correspond to variables %s and %s.\n",
+         CCTK_VarName(EvolvedIndex),CCTK_VarName(RHSIndexSlow));
+  printf("The pointer to EvolvedVariableIndex: %p\n",
+         EvolvedVariableIndex);
+#endif
+
+  if (!(CCTK_GroupTypeFromVarI(EvolvedIndex)==CCTK_GF))
+  {
+    CCTK_VWarn(0,__LINE__,__FILE__,CCTK_THORNSTRING,
+               "The variable %s is not a GF and so should "
+               "not be registered with MoLRegisterEvolved.", 
+               CCTK_VarName(EvolvedIndex)); 
+  }
+
+  if (!(CCTK_GroupTypeFromVarI(RHSIndexSlow)==CCTK_GF))
+  {
+    CCTK_VWarn(0,__LINE__,__FILE__,CCTK_THORNSTRING,
+               "The rhs variable %s is not a GF and so should "
+               "not be registered with MoLRegisterEvolved.", 
+               CCTK_VarName(RHSIndexSlow)); 
+  }    
+
+  if (!(CCTK_VarTypeI(EvolvedIndex)==CCTK_VARIABLE_REAL))
+  {
+    CCTK_VWarn(0,__LINE__,__FILE__,CCTK_THORNSTRING,
+               "The variable %s is not of type CCTK_REAL and so "
+               "should not be registered with MoLRegisterEvolved.", 
+               CCTK_VarName(EvolvedIndex)); 
+  }
+
+  if (!(CCTK_VarTypeI(RHSIndexSlow)==CCTK_VARIABLE_REAL))
+  {
+    CCTK_VWarn(0,__LINE__,__FILE__,CCTK_THORNSTRING,
+               "The rhs variable %s is not of type CCTK_REAL and "
+               "so should not be registered with MoLRegisterEvolved.", 
+               CCTK_VarName(RHSIndexSlow));
+  }
+
+  numtimelevs1 = CCTK_MaxTimeLevelsVI(EvolvedIndex);
+  numtimelevs2 = CCTK_MaxTimeLevelsVI(RHSIndexSlow);
+  
+  if ( (numtimelevs1 < 0) || (numtimelevs2 < 0) ) 
+  {
+    CCTK_VWarn(1,__LINE__,__FILE__,CCTK_THORNSTRING,
+               "Warning for variable index %i", EvolvedIndex); 
+    CCTK_WARN(0, "The index passed does not correspond to a GF.");
+  }
+
+  if (numtimelevs1 < 2)
+  {
+    CCTK_VWarn(1,__LINE__,__FILE__,CCTK_THORNSTRING,
+               "Warning for variable index %i name %s", 
+               EvolvedIndex, CCTK_VarName(EvolvedIndex));
+    CCTK_WARN(0, "The GF passed only has one timelevel. "
+              "It must have at least two.");
+  }
+
+  varused = 0;
+
+  for (index = 0; (index < MoLNumEvolvedVariablesSlow)&&(!varused); index++)
+  {
+    varused = (EvolvedIndex == EvolvedVariableIndexSlow[index]);
+#ifdef MOLDEBUG
+    printf("Registering %d. Checking index %d which is %d\n",EvolvedIndex,
+           index,EvolvedVariableIndexSlow[index]);
+#endif
+  }
+
+  if (varused)
+  {
+
+    CCTK_VWarn(2,__LINE__,__FILE__,CCTK_THORNSTRING,
+               "The GF %s has already been registered "
+               "as an evolved variable with RHS GF %s. "
+               "The attempt to register with RHS GF %s will be ignored",
+               CCTK_VarName(EvolvedIndex),
+               CCTK_VarName(RHSVariableIndexSlow[index-1]),
+               CCTK_VarName(RHSIndexSlow));
+
+  }
+  else
+  {
+
+    if (MoLNumEvolvedVariablesSlow+1 > MoL_Num_Evolved_Vars_Slow)
+    {
+      CCTK_WARN(0,"You have tried to register more evolved "
+                "variables than the accumulator parameter "
+                "MoL_Num_Evolved_Variables allows. Check that "
+                "you are accumulating onto this parameter correctly");
+    }
+  
+    EvolvedVariableIndexSlow[MoLNumEvolvedVariablesSlow] = EvolvedIndex;
+    RHSVariableIndexSlow[MoLNumEvolvedVariablesSlow] = RHSIndexSlow;
+  
+    MoLNumEvolvedVariablesSlow++;
+
+#ifdef MOLDEBUG
+    printf("The max number is now %d. Just added %d (%s).\n",
+           MoLNumEvolvedVariablesSlow, EvolvedIndex,
+           CCTK_VarName(EvolvedIndex));
+#endif
+
+  }
+
+  varused = -1;
+
+  for (index = 0; (index < MoLNumConstrainedVariables)&&(!(varused+1)); 
+       index++)
+  {
+    if (EvolvedIndex == ConstrainedVariableIndex[index])
+    {
+      varused = index;
+    }
+    
+  }
+  
+  if ((varused+1))
+  {
+    for (index = varused; index < MoLNumConstrainedVariables-1; index++)
+    {
+      ConstrainedVariableIndex[index] = ConstrainedVariableIndex[index+1];
+    }
+    MoLNumConstrainedVariables--;
+  }
+    
+  varused = -1;
+
+  for (index = 0; (index < MoLNumSandRVariables)&&(!(varused+1)); index++)
+  {
+    if (EvolvedIndex == SandRVariableIndex[index])
+    {
+      varused = index;
+    }
+    
+#ifdef MOLDEBUG
+    printf("Checking SandR var %d. Index %d (evolvedindex %d).\n",
+           index, SandRVariableIndex[index], EvolvedIndex);    
+#endif
+
+  }
+  
+  if ((varused+1))
+  {
+    for (index = varused; index < MoLNumSandRVariables-1; index++)
+    {
+      SandRVariableIndex[index] = SandRVariableIndex[index+1];
+
+#ifdef MOLDEBUG
+    printf("The registered evolved variable was SandR."
+           " Now index %d is %d (%s).\n",
+           index, SandRVariableIndex[index], 
+           CCTK_VarName(SandRVariableIndex[index]));    
+#endif
+
+    }
+    MoLNumSandRVariables--;
+  }
+  
+  return 0;
+
+}
+
+
+
+
  /*@@
    @routine    MoL_RegisterConstrained
    @date       Thu May 30 12:35:58 2002
@@ -1256,6 +1643,51 @@ CCTK_INT MoL_RegisterEvolvedRealGroup(CCTK_INT EvolvedGroupIndex,
 }
 
 
+CCTK_INT MoL_RegisterEvolvedRealGroupSlow(CCTK_INT EvolvedGroupIndex, 
+                                          CCTK_INT RHSGroupIndexSlow)
+{
+
+  CCTK_INT EvolvedGroupFirstVar, RHSGroupFirstVar, GroupNumVars, retval;
+  CCTK_INT EvolvedVar, RHSVar;
+  
+  EvolvedGroupFirstVar = CCTK_FirstVarIndexI(EvolvedGroupIndex);
+  if (EvolvedGroupFirstVar < 0)
+  {
+    CCTK_VWarn(0, __LINE__, __FILE__, CCTK_THORNSTRING,
+               "Evolved group index %i is not a real group index.",
+              EvolvedGroupIndex);
+  }
+  
+  RHSGroupFirstVar = CCTK_FirstVarIndexI(RHSGroupIndexSlow);
+  if (RHSGroupFirstVar < 0)
+  {
+    CCTK_VWarn(0, __LINE__, __FILE__, CCTK_THORNSTRING,
+               "RHS group index %d is not a real group index.",
+              RHSGroupIndexSlow);
+  }
+
+  GroupNumVars = CCTK_NumVarsInGroupI(EvolvedGroupIndex);
+  if (CCTK_NumVarsInGroupI(RHSGroupIndexSlow) != GroupNumVars)
+  {
+    CCTK_VWarn(0, __LINE__, __FILE__, CCTK_THORNSTRING,
+               "There are a different number of variables in "
+               "evolved group %d and RHS group %d.", 
+               EvolvedGroupIndex, RHSGroupIndexSlow);
+  }
+  
+  retval = 0;
+  
+  for (EvolvedVar = EvolvedGroupFirstVar, RHSVar = RHSGroupFirstVar;
+       EvolvedVar < EvolvedGroupFirstVar + GroupNumVars;
+       EvolvedVar++, RHSVar++)
+  {
+    retval += MoL_RegisterEvolvedRealSlow(EvolvedVar, RHSVar);
+  }
+  
+  return retval;
+}
+
+
 CCTK_INT MoL_RegisterConstrainedRealGroup(CCTK_INT ConstrainedGroupIndex)
 {
   
@@ -2741,6 +3173,7 @@ CCTK_INT MoL_RegisterSaveAndRestoreComplexArrayGroup(CCTK_INT SandRGroupIndex)
    @desc 
    Given the index of a registered evolved GF,
    returns the corresponding RHS GF.
+   Checks both regular and slow evolved variables.
    @enddesc 
    @calls     
    @calledby   
@@ -2769,6 +3202,17 @@ CCTK_INT MoL_QueryEvolvedRHS(CCTK_INT EvolvedIndex)
     }
   }
   
+  // If we haven't found any registered variable that matches,
+  // we also try to check whether this has been registered with the slow multirate sector!
+  
+  for (index = 0; index < MoLNumEvolvedVariablesSlow; index++)
+  {
+    if (EvolvedVariableIndexSlow[index] == EvolvedIndex)
+    {
+      return RHSVariableIndexSlow[index];
+    }
+  }
+  
   return -1;
 
 }
diff --git a/src/SetTime.c b/src/SetTime.c
index c04cc0b..32e893c 100644
--- a/src/SetTime.c
+++ b/src/SetTime.c
@@ -298,6 +298,63 @@ void MoL_ResetTime(CCTK_ARGUMENTS)
          + ((alpha_array87[substep] - 1)
             * (* Original_Delta_Time) / cctkGH->cctk_timefac));
   }
+  else if (CCTK_EQUALS(ODE_Method,"RK2-MR-2:1"))
+  {
+    const int substep = MoL_Intermediate_Steps - (* MoL_Intermediate_Step);
+
+    CCTK_REAL dt = (*Original_Delta_Time)/cctkGH->cctk_timefac;
+    switch (substep)
+    {
+      case 1:
+      case 2:
+        dt *= 0.5;
+        break;
+      default:
+        dt = 0;
+    }
+    cctkGH->cctk_time = (*Original_Time) - dt;
+  }
+  else if (CCTK_EQUALS(ODE_Method,"RK4-MR-2:1"))
+  {
+    const int substep = MoL_Intermediate_Steps - (* MoL_Intermediate_Step);
+
+    CCTK_REAL dt = (*Original_Delta_Time)/cctkGH->cctk_timefac;
+    switch (substep)
+    {
+      case 1:
+      case 2:
+        dt *= 3.0/4.0;
+        break;
+      case 3:
+      case 4:
+      case 5:
+        dt *= 1.0/2.0;
+        break;
+      case 6:
+      case 7:
+        dt *= 1.0/4.0;
+        break;
+      default:
+        dt = 0;
+    }
+    cctkGH->cctk_time = (*Original_Time) - dt;
+  }
+  else if (CCTK_EQUALS(ODE_Method,"RK4-RK2"))
+  {
+    const int substep = MoL_Intermediate_Steps - (* MoL_Intermediate_Step);
+
+    CCTK_REAL dt = (*Original_Delta_Time)/cctkGH->cctk_timefac;
+    switch (substep)
+    {
+      case 1:
+      case 2:
+        dt *= 0.5;
+        break;
+      default:
+        dt = 0;
+    }
+    cctkGH->cctk_time = (*Original_Time) - dt;
+  }
 #ifdef MOLDEBUG
   printf("MoL has once more reset t (%d): %f.\n", 
          *MoL_Intermediate_Step, cctkGH->cctk_time);
diff --git a/src/Startup.c b/src/Startup.c
index 06ae44d..e5a2d5e 100644
--- a/src/Startup.c
+++ b/src/Startup.c
@@ -26,11 +26,14 @@ CCTK_FILEVERSION(CactusBase_MoL_Startup_c);
  ********************************************************************/
 
 CCTK_INT *EvolvedVariableIndex = NULL;
+CCTK_INT *EvolvedVariableIndexSlow = NULL;
 CCTK_INT *RHSVariableIndex = NULL;
+CCTK_INT *RHSVariableIndexSlow = NULL;
 CCTK_INT *ConstrainedVariableIndex = NULL;
 CCTK_INT *SandRVariableIndex = NULL;
 
 CCTK_INT MoLNumEvolvedVariables = 0;
+CCTK_INT MoLNumEvolvedVariablesSlow = 0;
 CCTK_INT MoLNumConstrainedVariables = 0;
 CCTK_INT MoLNumSandRVariables = 0;
 
diff --git a/src/StepSize.c b/src/StepSize.c
index cb5499e..b20c752 100644
--- a/src/StepSize.c
+++ b/src/StepSize.c
@@ -358,6 +358,10 @@ void MoL_FinishLoop(CCTK_ARGUMENTS)
 
   /* Keep time step size unchanged */
   *MoL_Stepsize_Bad = 0;
+  
+  // Set these flags to ONE outside of MoL-loop!
+  *MoL_SlowPostStep = 1;
+  *MoL_SlowStep = 1;
 }
 
 /********************************************************************
diff --git a/src/make.code.defn b/src/make.code.defn
index 926dbbb..80e1b56 100644
--- a/src/make.code.defn
+++ b/src/make.code.defn
@@ -21,8 +21,11 @@ SRCS = 	ChangeType.c \
 	SandR.c \
 	SetTime.c \
 	Startup.c \
-	StepSize.c
-
+	StepSize.c \
+	RK2-MR-2_1.c \
+	RK4-MR-2_1.c \
+	RK4-RK2.c
+	
 # Subdirectories containing source files
 SUBDIRS =