diff --git a/src/evolve.h b/src/evolve.h
index 6eec7a40..f6a0dee4 100644
--- a/src/evolve.h
+++ b/src/evolve.h
@@ -14,9 +14,13 @@
 #define EULER 1
 #define RUNGEKUTTA 2
 
+#define NO_INTEGRATION 0
+#define FORWARD_INTEGRATION 1
+#define BACKWARD_INTEGRATION -1
+
 /* @cond DOXYGEN_OVERRIDE */
 
-double fdTrapezoidalArea(double,double,double);
+double fdTrapezoidalArea(double, double, double);
 void PropertiesAuxiliary(BODY *, CONTROL *, SYSTEM *, UPDATE *);
 void fdGetUpdateInfo(BODY *, CONTROL *, SYSTEM *, UPDATE *,
                      fnUpdateVariable ***);
diff --git a/src/options.c b/src/options.c
index 09ccee1a..6e3b30e6 100644
--- a/src/options.c
+++ b/src/options.c
@@ -1078,7 +1078,6 @@ void ReadUnitTemp(CONTROL *control, FILES *files, OPTIONS *options, int iFile) {
 
 void ReadSystemName(CONTROL *control, FILES *files, OPTIONS *options,
                     SYSTEM *system, int iFile) {
-  /* System Name */
   int lTmp = -1;
   char cTmp[OPTLEN];
 
@@ -1092,6 +1091,28 @@ void ReadSystemName(CONTROL *control, FILES *files, OPTIONS *options,
   }
 }
 
+void ReadSystemAge(BODY *body, CONTROL *control, FILES *files, OPTIONS *options,
+                   SYSTEM *system, int iFile) {
+  int lTmp = -1;
+  double dTmp;
+
+  AddOptionDouble(files->Infile[iFile].cIn, options->cName, &dTmp, &lTmp,
+                  control->Io.iVerbose);
+  if (lTmp >= 0) {
+    CheckDuplication(files, options, files->Infile[iFile].cIn, lTmp,
+                     control->Io.iVerbose);
+    if (dTmp < 0) {
+      system->dAge = dTmp * dNegativeDouble(*options, files->Infile[iFile].cIn,
+                                            control->Io.iVerbose);
+    } else {
+      system->dAge = dTmp * fdUnitsTime(control->Units[iFile].iTime);
+    }
+    UpdateFoundOption(&files->Infile[iFile], options, lTmp, iFile);
+  } else if (iFile > 0) {
+    AssignDefaultDouble(options, &system->dAge, files->iNumInputs);
+  }
+}
+
 void ReadBodyFileNames(BODY **body, CONTROL *control, FILES *files,
                        OPTIONS *options, char *cFile, char ***saBodyFiles,
                        int *iNumLines, int *iaLines) {
@@ -1442,8 +1463,6 @@ void ReadOutputTime(BODY *body, CONTROL *control, FILES *files,
   }
 }
 
-/* Backward integration stop time */
-
 void ReadStopTime(BODY *body, CONTROL *control, FILES *files, OPTIONS *options,
                   SYSTEM *system, int iFile) {
   /* This parameter can exist in any file, but only once */
@@ -1470,6 +1489,32 @@ void ReadStopTime(BODY *body, CONTROL *control, FILES *files, OPTIONS *options,
   }
 }
 
+void ReadStopAge(BODY *body, CONTROL *control, FILES *files, OPTIONS *options,
+                 SYSTEM *system, int iFile) {
+  /* This parameter can exist in any file, but only once */
+  int lTmp = -1;
+  double dTmp;
+
+  AddOptionDouble(files->Infile[iFile].cIn, options->cName, &dTmp, &lTmp,
+                  control->Io.iVerbose);
+  if (lTmp >= 0) {
+    /* Option was found */
+    CheckDuplication(files, options, files->Infile[iFile].cIn, lTmp,
+                     control->Io.iVerbose);
+    if (dTmp < 0) {
+      if (control->Io.iVerbose >= VERBERR) {
+        fprintf(stderr, "ERROR: %s must be greater than 0.\n", options->cName);
+      }
+      LineExit(files->Infile[iFile].cIn, lTmp);
+    }
+    /* Convert stop time to cgs */
+    control->Evolve.dStopAge = dTmp * fdUnitsTime(control->Units[iFile].iTime);
+    UpdateFoundOption(&files->Infile[iFile], options, lTmp, iFile);
+  } else {
+    AssignDefaultDouble(options, &control->Evolve.dStopAge, files->iNumInputs);
+  }
+}
+
 /* Integration timestep */
 
 void ReadTimeStep(BODY *body, CONTROL *control, FILES *files, OPTIONS *options,
@@ -2617,7 +2662,7 @@ void fvAllocateOutputArrays(char ****saMatch, char ***saOutput, int **baNeg,
 
   for (iIndex = 0; iIndex < iNumArgs; iIndex++) {
     (*saMatch)[iIndex] =
-          malloc(MAXARRAY * sizeof(char*)); // Could be this many matches
+          malloc(MAXARRAY * sizeof(char *)); // Could be this many matches
     for (iMatch = 0; iMatch < MAXARRAY; iMatch++) {
       (*saMatch)[iIndex][iMatch] = NULL;
     }
@@ -3712,7 +3757,7 @@ void InitializeOptionsGeneral(OPTIONS *options, fnReadOption fnRead[]) {
    */
 
   fvFormattedString(&options[OPT_AGE].cName, "dAge");
-  fvFormattedString(&options[OPT_AGE].cDescr, "System Age");
+  fvFormattedString(&options[OPT_AGE].cDescr, "Body's Age");
   fvFormattedString(&options[OPT_AGE].cDefault, "0");
   fvFormattedString(&options[OPT_AGE].cNeg, "Gyr");
   fvFormattedString(&options[OPT_AGE].cDimension, "time");
@@ -3724,6 +3769,20 @@ void InitializeOptionsGeneral(OPTIONS *options, fnReadOption fnRead[]) {
   options[OPT_AGE].iFileType  = 2;
   fnRead[OPT_AGE]             = &ReadAge;
 
+  int iOpt = OPT_SYSTEMAGE;
+  fvFormattedString(&options[iOpt].cName, "dSystemAge");
+  fvFormattedString(&options[iOpt].cDescr, "System Age");
+  fvFormattedString(&options[iOpt].cDefault, "0");
+  fvFormattedString(&options[iOpt].cNeg, "1 Gyr");
+  fvFormattedString(&options[iOpt].cDimension, "time");
+  options[iOpt].dDefault   = 0;
+  options[iOpt].iType      = 2;
+  options[iOpt].iModuleBit = 0;
+  options[iOpt].bNeg       = 1;
+  options[iOpt].dNeg       = 1e9 * YEARSEC;
+  options[iOpt].iFileType  = 2;
+  fnRead[iOpt]             = &ReadSystemAge;
+
   fvFormattedString(&options[OPT_ALBEDOGLOBAL].cName, "dAlbedoGlobal");
   fvFormattedString(&options[OPT_ALBEDOGLOBAL].cDescr,
                     "Globally averaged albedo");
@@ -3805,6 +3864,20 @@ void InitializeOptionsGeneral(OPTIONS *options, fnReadOption fnRead[]) {
   options[OPT_STOPTIME].iFileType  = 2;
   fnRead[OPT_STOPTIME]             = &ReadStopTime;
 
+  iOpt = OPT_STOPAGE;
+  fvFormattedString(&options[iOpt].cName, "dStopAge");
+  fvFormattedString(&options[iOpt].cDescr, "Age to stop integration");
+  fvFormattedString(&options[iOpt].cDefault, "10 Gigayears");
+  fvFormattedString(&options[iOpt].cNeg, "Years");
+  fvFormattedString(&options[iOpt].cDimension, "time");
+  options[iOpt].dDefault   = 1e10 * YEARSEC;
+  options[iOpt].iType      = 2;
+  options[iOpt].iModuleBit = 0;
+  options[iOpt].bNeg       = 1;
+  options[iOpt].dNeg       = YEARSEC;
+  options[iOpt].iFileType  = 2;
+  fnRead[iOpt]             = &ReadStopAge;
+
   fvFormattedString(&options[OPT_TIMESTEP].cName, "dTimeStep");
   fvFormattedString(&options[OPT_TIMESTEP].cDescr, "Integration Timestep");
   fvFormattedString(&options[OPT_TIMESTEP].cDefault, "1 year");
diff --git a/src/options.h b/src/options.h
index 18f67ffd..f11c1693 100644
--- a/src/options.h
+++ b/src/options.h
@@ -30,6 +30,7 @@
 // Regular Options
 
 #define OPT_AGE 100
+#define OPT_SYSTEMAGE 102
 #define OPT_ALBEDOGLOBAL 105
 
 #define OPT_BACK 110
@@ -37,6 +38,7 @@
 #define OPT_ETA 130
 #define OPT_OUTPUTTIME 140
 #define OPT_STOPTIME 150
+#define OPT_STOPAGE 155
 #define OPT_TIMESTEP 160
 #define OPT_VARDT 170
 #define OPT_BODYNAME 180
@@ -128,8 +130,7 @@ void ReadOptions(BODY **, CONTROL *, FILES *, MODULE *, OPTIONS *, OUTPUT *,
                  SYSTEM *, UPDATE **, fnReadOption *, char[]);
 
 double dNegativeDouble(OPTIONS, char[], int);
-void AddOptionStringArray(char[], char[], char***, int *, int *,
-                          int *, int);
+void AddOptionStringArray(char[], char[], char ***, int *, int *, int *, int);
 void AddOptionDoubleArray(char[], char[], double *, int *, int *, int *, int);
 void NotPrimaryInput(int, char[], char[], int, int);
 void AddOptionDouble(char[], char[], double *, int *, int);
diff --git a/src/update.c b/src/update.c
index 34e221aa..195d0f18 100644
--- a/src/update.c
+++ b/src/update.c
@@ -20,10 +20,10 @@ void InitializeUpdateBodyPerts(CONTROL *control, UPDATE *update, int iBody) {
 
 void InitializeUpdateTmpBody(BODY *body, CONTROL *control, MODULE *module,
                              UPDATE *update, int iBody) {
-  int jBody,iModule;
+  int jBody, iModule;
 
-  for (jBody=0;jBody<control->Evolve.iNumBodies;jBody++) {
-   control->Evolve.tmpBody[jBody].cName = NULL;           
+  for (jBody = 0; jBody < control->Evolve.iNumBodies; jBody++) {
+    control->Evolve.tmpBody[jBody].cName = NULL;
   }
 
   for (iModule = 0; iModule < module->iNumModules[iBody]; iModule++) {
@@ -65,21 +65,22 @@ void UpdateCopy(UPDATE *dest, UPDATE *src, int iNumBodies) {
   }
 }
 
-/* Someday we'll fix it so that the each primary variable just calls this function. The current impediment
-is the fnFinalizeUpdateModule -- it doesn't exist yet. Currently each Primary Variable (PV) has its own
-FinalizeUpdate function that makes passing it around more difficult. The "right" way to solve this is to 
-create a generic fnFinalizeUpdateModule function that is 3 dimensions: iBody, iModule, iPV. We'd have to 
-figure out the accounting for iPV, which is why I never did it. In the meantime, it sucks to have to
-make any changes to InitializeUpdate because you have to go through ALL the PVs. -- RB, 21 Feb 22 
-
-void AddPrimaryVariable(BODY *body,CONTROL *control,MODULE *module,UPDATE *update,fnUpdateVariable ***fnUpdate,
-                  double *dPrimaryVariable, double *dTmpPrimaryVariable,fnFinalizeUpdateModule *finalize,
-                   int *iPrimaryVariable,int iNumPrimaryVariable,int iBody,int *iVar,int iID) {
-    *iPrimaryVariable = -1;
-    if (iNumPrimaryVariable) {
-      *iPrimaryVariable    = *iVar;
-      update[iBody].iaVar[*iVar]    = iID;
-      update[iBody].iNumEqns[*iVar] = iNumPrimaryVariable;
+/* Someday we'll fix it so that the each primary variable just calls this
+function. The current impediment is the fnFinalizeUpdateModule -- it doesn't
+exist yet. Currently each Primary Variable (PV) has its own FinalizeUpdate
+function that makes passing it around more difficult. The "right" way to solve
+this is to create a generic fnFinalizeUpdateModule function that is 3
+dimensions: iBody, iModule, iPV. We'd have to figure out the accounting for iPV,
+which is why I never did it. In the meantime, it sucks to have to make any
+changes to InitializeUpdate because you have to go through ALL the PVs. -- RB,
+21 Feb 22
+
+void AddPrimaryVariable(BODY *body,CONTROL *control,MODULE *module,UPDATE
+*update,fnUpdateVariable ***fnUpdate, double *dPrimaryVariable, double
+*dTmpPrimaryVariable,fnFinalizeUpdateModule *finalize, int *iPrimaryVariable,int
+iNumPrimaryVariable,int iBody,int *iVar,int iID) { *iPrimaryVariable = -1; if
+(iNumPrimaryVariable) { *iPrimaryVariable    = *iVar; update[iBody].iaVar[*iVar]
+= iID; update[iBody].iNumEqns[*iVar] = iNumPrimaryVariable;
       update[iBody].pdVar[*iVar]    = dPrimaryVariable;
       update[iBody].iNumBodies[*iVar] =
             malloc(iNumPrimaryVariable * sizeof(int));
@@ -115,8 +116,9 @@ void AddPrimaryVariable(BODY *body,CONTROL *control,MODULE *module,UPDATE *updat
         finalize(body, update, &iEqn,iVar, iBody, iFoo);
       }
 
-      (*fnUpdate)[iBody][*iVar]        = malloc(iEqn * sizeof(fnUpdateVariable));
-      update[iBody].daDerivProc[*iVar] = malloc(iEqn * sizeof(double));
+      (*fnUpdate)[iBody][*iVar]        = malloc(iEqn *
+sizeof(fnUpdateVariable)); update[iBody].daDerivProc[*iVar] = malloc(iEqn *
+sizeof(double));
       (*iVar)++;
     }
 }
@@ -237,33 +239,34 @@ void InitializeUpdate(BODY *body, CONTROL *control, MODULE *module,
     update[iBody].iaBody     = malloc(update[iBody].iNumVars * sizeof(int **));
 
     // May also have to allocate space for the temp UPDATE
+    control->Evolve.tmpUpdate[iBody].iaVar =
+          malloc(update[iBody].iNumVars * sizeof(int));
+    control->Evolve.tmpUpdate[iBody].iNumEqns =
+          malloc(update[iBody].iNumVars * sizeof(double));
+    control->Evolve.tmpUpdate[iBody].iaType =
+          malloc(update[iBody].iNumVars * sizeof(int *));
+    control->Evolve.tmpUpdate[iBody].iaModule =
+          malloc(update[iBody].iNumVars * sizeof(int *));
+    control->Evolve.tmpUpdate[iBody].pdVar =
+          malloc(update[iBody].iNumVars * sizeof(double *));
+    control->Evolve.tmpUpdate[iBody].daDeriv =
+          malloc(update[iBody].iNumVars * sizeof(double *));
+    control->Evolve.tmpUpdate[iBody].daDerivProc =
+          malloc(update[iBody].iNumVars * sizeof(double *));
+    control->Evolve.tmpUpdate[iBody].iNumBodies =
+          malloc(update[iBody].iNumVars * sizeof(int *));
+    control->Evolve.tmpUpdate[iBody].iaBody =
+          malloc(update[iBody].iNumVars * sizeof(int **));
     if (control->Evolve.iOneStep == RUNGEKUTTA) {
-      control->Evolve.tmpUpdate[iBody].iaVar =
-            malloc(update[iBody].iNumVars * sizeof(int));
-      control->Evolve.tmpUpdate[iBody].iNumEqns =
-            malloc(update[iBody].iNumVars * sizeof(double));
-      control->Evolve.tmpUpdate[iBody].iaType =
-            malloc(update[iBody].iNumVars * sizeof(int *));
-      control->Evolve.tmpUpdate[iBody].iaModule =
-            malloc(update[iBody].iNumVars * sizeof(int *));
-      control->Evolve.tmpUpdate[iBody].pdVar =
-            malloc(update[iBody].iNumVars * sizeof(double *));
-      control->Evolve.tmpUpdate[iBody].daDeriv =
-            malloc(update[iBody].iNumVars * sizeof(double *));
-      control->Evolve.tmpUpdate[iBody].daDerivProc =
-            malloc(update[iBody].iNumVars * sizeof(double *));
-      control->Evolve.tmpUpdate[iBody].iNumBodies =
-            malloc(update[iBody].iNumVars * sizeof(int *));
-      control->Evolve.tmpUpdate[iBody].iaBody =
-            malloc(update[iBody].iNumVars * sizeof(int **));
-    }
-    for (iSubStep = 0; iSubStep < 4; iSubStep++) {
-      control->Evolve.daDeriv[iSubStep][iBody] =
-            malloc(update[iBody].iNumVars * sizeof(double));
-      control->Evolve.daDerivProc[iSubStep][iBody] =
-            malloc(update[iBody].iNumVars * sizeof(double*));
+      for (iSubStep = 0; iSubStep < 4; iSubStep++) {
+        control->Evolve.daDeriv[iSubStep][iBody] =
+              malloc(update[iBody].iNumVars * sizeof(double));
+        control->Evolve.daDerivProc[iSubStep][iBody] =
+              malloc(update[iBody].iNumVars * sizeof(double *));
+      }
     }
 
+
     /* Now we malloc some pointers, and perform some initializations for the
        UPDATE struct based on the primary variables required for each's
        planet's assigned modules. */
@@ -297,9 +300,9 @@ void InitializeUpdate(BODY *body, CONTROL *control, MODULE *module,
               malloc(update[iBody].iNumVelX * sizeof(int));
         control->Evolve.tmpUpdate[iBody].iaBody[iVar] =
               malloc(update[iBody].iNumVelX * sizeof(int *));
-        for (iSubStep=0; iSubStep < 4; iSubStep++) {
-            control->Evolve.daDerivProc[iSubStep][iBody][iVar] =
-                  malloc(update[iBody].iNumEqns[iVar] * sizeof(double));
+        for (iSubStep = 0; iSubStep < 4; iSubStep++) {
+          control->Evolve.daDerivProc[iSubStep][iBody][iVar] =
+                malloc(update[iBody].iNumEqns[iVar] * sizeof(double));
         }
       }
 
@@ -341,9 +344,9 @@ void InitializeUpdate(BODY *body, CONTROL *control, MODULE *module,
               malloc(update[iBody].iNumVelY * sizeof(int));
         control->Evolve.tmpUpdate[iBody].iaBody[iVar] =
               malloc(update[iBody].iNumVelY * sizeof(int *));
-        for (iSubStep=0; iSubStep < 4; iSubStep++) {
-            control->Evolve.daDerivProc[iSubStep][iBody][iVar] =
-                  malloc(update[iBody].iNumEqns[iVar] * sizeof(double));
+        for (iSubStep = 0; iSubStep < 4; iSubStep++) {
+          control->Evolve.daDerivProc[iSubStep][iBody][iVar] =
+                malloc(update[iBody].iNumEqns[iVar] * sizeof(double));
         }
       }
 
@@ -385,9 +388,9 @@ void InitializeUpdate(BODY *body, CONTROL *control, MODULE *module,
               malloc(update[iBody].iNumVelZ * sizeof(int));
         control->Evolve.tmpUpdate[iBody].iaBody[iVar] =
               malloc(update[iBody].iNumVelZ * sizeof(int *));
-        for (iSubStep=0; iSubStep < 4; iSubStep++) {
-            control->Evolve.daDerivProc[iSubStep][iBody][iVar] =
-                  malloc(update[iBody].iNumEqns[iVar] * sizeof(double));
+        for (iSubStep = 0; iSubStep < 4; iSubStep++) {
+          control->Evolve.daDerivProc[iSubStep][iBody][iVar] =
+                malloc(update[iBody].iNumEqns[iVar] * sizeof(double));
         }
       }
 
@@ -399,7 +402,7 @@ void InitializeUpdate(BODY *body, CONTROL *control, MODULE *module,
 
       (*fnUpdate)[iBody][iVar]        = malloc(iEqn * sizeof(fnUpdateVariable));
       update[iBody].daDerivProc[iVar] = malloc(iEqn * sizeof(double));
-       iVar++;
+      iVar++;
     }
 
     update[iBody].iPositionX = -1;
@@ -430,9 +433,9 @@ void InitializeUpdate(BODY *body, CONTROL *control, MODULE *module,
               malloc(update[iBody].iNumPositionX * sizeof(int));
         control->Evolve.tmpUpdate[iBody].iaBody[iVar] =
               malloc(update[iBody].iNumPositionX * sizeof(int *));
-        for (iSubStep=0; iSubStep < 4; iSubStep++) {
-            control->Evolve.daDerivProc[iSubStep][iBody][iVar] =
-                  malloc(update[iBody].iNumEqns[iVar] * sizeof(double));
+        for (iSubStep = 0; iSubStep < 4; iSubStep++) {
+          control->Evolve.daDerivProc[iSubStep][iBody][iVar] =
+                malloc(update[iBody].iNumEqns[iVar] * sizeof(double));
         }
       }
 
@@ -475,9 +478,9 @@ void InitializeUpdate(BODY *body, CONTROL *control, MODULE *module,
               malloc(update[iBody].iNumPositionY * sizeof(int));
         control->Evolve.tmpUpdate[iBody].iaBody[iVar] =
               malloc(update[iBody].iNumPositionY * sizeof(int *));
-        for (iSubStep=0; iSubStep < 4; iSubStep++) {
-            control->Evolve.daDerivProc[iSubStep][iBody][iVar] =
-                  malloc(update[iBody].iNumEqns[iVar] * sizeof(double));
+        for (iSubStep = 0; iSubStep < 4; iSubStep++) {
+          control->Evolve.daDerivProc[iSubStep][iBody][iVar] =
+                malloc(update[iBody].iNumEqns[iVar] * sizeof(double));
         }
       }
 
@@ -520,9 +523,9 @@ void InitializeUpdate(BODY *body, CONTROL *control, MODULE *module,
               malloc(update[iBody].iNumPositionZ * sizeof(int));
         control->Evolve.tmpUpdate[iBody].iaBody[iVar] =
               malloc(update[iBody].iNumPositionZ * sizeof(int *));
-        for (iSubStep=0; iSubStep < 4; iSubStep++) {
-            control->Evolve.daDerivProc[iSubStep][iBody][iVar] =
-                  malloc(update[iBody].iNumEqns[iVar] * sizeof(double));
+        for (iSubStep = 0; iSubStep < 4; iSubStep++) {
+          control->Evolve.daDerivProc[iSubStep][iBody][iVar] =
+                malloc(update[iBody].iNumEqns[iVar] * sizeof(double));
         }
       }
 
@@ -567,9 +570,9 @@ void InitializeUpdate(BODY *body, CONTROL *control, MODULE *module,
               malloc(update[iBody].iNumWaterMassMOAtm * sizeof(int));
         control->Evolve.tmpUpdate[iBody].iaBody[iVar] =
               malloc(update[iBody].iNumWaterMassMOAtm * sizeof(int *));
-        for (iSubStep=0; iSubStep < 4; iSubStep++) {
-            control->Evolve.daDerivProc[iSubStep][iBody][iVar] =
-                  malloc(update[iBody].iNumEqns[iVar] * sizeof(double));
+        for (iSubStep = 0; iSubStep < 4; iSubStep++) {
+          control->Evolve.daDerivProc[iSubStep][iBody][iVar] =
+                malloc(update[iBody].iNumEqns[iVar] * sizeof(double));
         }
       }
 
@@ -612,9 +615,9 @@ void InitializeUpdate(BODY *body, CONTROL *control, MODULE *module,
               malloc(update[iBody].iNumWaterMassSol * sizeof(int));
         control->Evolve.tmpUpdate[iBody].iaBody[iVar] =
               malloc(update[iBody].iNumWaterMassSol * sizeof(int *));
-        for (iSubStep=0; iSubStep < 4; iSubStep++) {
-            control->Evolve.daDerivProc[iSubStep][iBody][iVar] =
-                  malloc(update[iBody].iNumEqns[iVar] * sizeof(double));
+        for (iSubStep = 0; iSubStep < 4; iSubStep++) {
+          control->Evolve.daDerivProc[iSubStep][iBody][iVar] =
+                malloc(update[iBody].iNumEqns[iVar] * sizeof(double));
         }
       }
 
@@ -657,9 +660,9 @@ void InitializeUpdate(BODY *body, CONTROL *control, MODULE *module,
               malloc(update[iBody].iNumSurfTemp * sizeof(int));
         control->Evolve.tmpUpdate[iBody].iaBody[iVar] =
               malloc(update[iBody].iNumSurfTemp * sizeof(int *));
-        for (iSubStep=0; iSubStep < 4; iSubStep++) {
-            control->Evolve.daDerivProc[iSubStep][iBody][iVar] =
-                  malloc(update[iBody].iNumEqns[iVar] * sizeof(double));
+        for (iSubStep = 0; iSubStep < 4; iSubStep++) {
+          control->Evolve.daDerivProc[iSubStep][iBody][iVar] =
+                malloc(update[iBody].iNumEqns[iVar] * sizeof(double));
         }
       }
 
@@ -702,9 +705,9 @@ void InitializeUpdate(BODY *body, CONTROL *control, MODULE *module,
               malloc(update[iBody].iNumSolidRadius * sizeof(int));
         control->Evolve.tmpUpdate[iBody].iaBody[iVar] =
               malloc(update[iBody].iNumSolidRadius * sizeof(int *));
-        for (iSubStep=0; iSubStep < 4; iSubStep++) {
-            control->Evolve.daDerivProc[iSubStep][iBody][iVar] =
-                  malloc(update[iBody].iNumEqns[iVar] * sizeof(double));
+        for (iSubStep = 0; iSubStep < 4; iSubStep++) {
+          control->Evolve.daDerivProc[iSubStep][iBody][iVar] =
+                malloc(update[iBody].iNumEqns[iVar] * sizeof(double));
         }
       }
 
@@ -747,9 +750,9 @@ void InitializeUpdate(BODY *body, CONTROL *control, MODULE *module,
               malloc(update[iBody].iNumPotTemp * sizeof(int));
         control->Evolve.tmpUpdate[iBody].iaBody[iVar] =
               malloc(update[iBody].iNumPotTemp * sizeof(int *));
-        for (iSubStep=0; iSubStep < 4; iSubStep++) {
-            control->Evolve.daDerivProc[iSubStep][iBody][iVar] =
-                  malloc(update[iBody].iNumEqns[iVar] * sizeof(double));
+        for (iSubStep = 0; iSubStep < 4; iSubStep++) {
+          control->Evolve.daDerivProc[iSubStep][iBody][iVar] =
+                malloc(update[iBody].iNumEqns[iVar] * sizeof(double));
         }
       }
 
@@ -792,9 +795,9 @@ void InitializeUpdate(BODY *body, CONTROL *control, MODULE *module,
               malloc(update[iBody].iNumOxygenMassMOAtm * sizeof(int));
         control->Evolve.tmpUpdate[iBody].iaBody[iVar] =
               malloc(update[iBody].iNumOxygenMassMOAtm * sizeof(int *));
-        for (iSubStep=0; iSubStep < 4; iSubStep++) {
-            control->Evolve.daDerivProc[iSubStep][iBody][iVar] =
-                  malloc(update[iBody].iNumEqns[iVar] * sizeof(double));
+        for (iSubStep = 0; iSubStep < 4; iSubStep++) {
+          control->Evolve.daDerivProc[iSubStep][iBody][iVar] =
+                malloc(update[iBody].iNumEqns[iVar] * sizeof(double));
         }
       }
 
@@ -806,7 +809,7 @@ void InitializeUpdate(BODY *body, CONTROL *control, MODULE *module,
 
       (*fnUpdate)[iBody][iVar]        = malloc(iEqn * sizeof(fnUpdateVariable));
       update[iBody].daDerivProc[iVar] = malloc(iEqn * sizeof(double));
-       iVar++;
+      iVar++;
     }
     /* until HERE */
     update[iBody].iOxygenMassSol = -1;
@@ -837,9 +840,9 @@ void InitializeUpdate(BODY *body, CONTROL *control, MODULE *module,
               malloc(update[iBody].iNumOxygenMassSol * sizeof(int));
         control->Evolve.tmpUpdate[iBody].iaBody[iVar] =
               malloc(update[iBody].iNumOxygenMassSol * sizeof(int *));
-        for (iSubStep=0; iSubStep < 4; iSubStep++) {
-            control->Evolve.daDerivProc[iSubStep][iBody][iVar] =
-                  malloc(update[iBody].iNumEqns[iVar] * sizeof(double));
+        for (iSubStep = 0; iSubStep < 4; iSubStep++) {
+          control->Evolve.daDerivProc[iSubStep][iBody][iVar] =
+                malloc(update[iBody].iNumEqns[iVar] * sizeof(double));
         }
       }
 
@@ -882,9 +885,9 @@ void InitializeUpdate(BODY *body, CONTROL *control, MODULE *module,
               malloc(update[iBody].iNumHydrogenMassSpace * sizeof(int));
         control->Evolve.tmpUpdate[iBody].iaBody[iVar] =
               malloc(update[iBody].iNumHydrogenMassSpace * sizeof(int *));
-        for (iSubStep=0; iSubStep < 4; iSubStep++) {
-            control->Evolve.daDerivProc[iSubStep][iBody][iVar] =
-                  malloc(update[iBody].iNumEqns[iVar] * sizeof(double));
+        for (iSubStep = 0; iSubStep < 4; iSubStep++) {
+          control->Evolve.daDerivProc[iSubStep][iBody][iVar] =
+                malloc(update[iBody].iNumEqns[iVar] * sizeof(double));
         }
       }
 
@@ -927,9 +930,9 @@ void InitializeUpdate(BODY *body, CONTROL *control, MODULE *module,
               malloc(update[iBody].iNumOxygenMassSpace * sizeof(int));
         control->Evolve.tmpUpdate[iBody].iaBody[iVar] =
               malloc(update[iBody].iNumOxygenMassSpace * sizeof(int *));
-        for (iSubStep=0; iSubStep < 4; iSubStep++) {
-            control->Evolve.daDerivProc[iSubStep][iBody][iVar] =
-                  malloc(update[iBody].iNumEqns[iVar] * sizeof(double));
+        for (iSubStep = 0; iSubStep < 4; iSubStep++) {
+          control->Evolve.daDerivProc[iSubStep][iBody][iVar] =
+                malloc(update[iBody].iNumEqns[iVar] * sizeof(double));
         }
       }
 
@@ -972,9 +975,9 @@ void InitializeUpdate(BODY *body, CONTROL *control, MODULE *module,
               malloc(update[iBody].iNumCO2MassMOAtm * sizeof(int));
         control->Evolve.tmpUpdate[iBody].iaBody[iVar] =
               malloc(update[iBody].iNumCO2MassMOAtm * sizeof(int *));
-        for (iSubStep=0; iSubStep < 4; iSubStep++) {
-            control->Evolve.daDerivProc[iSubStep][iBody][iVar] =
-                  malloc(update[iBody].iNumEqns[iVar] * sizeof(double));
+        for (iSubStep = 0; iSubStep < 4; iSubStep++) {
+          control->Evolve.daDerivProc[iSubStep][iBody][iVar] =
+                malloc(update[iBody].iNumEqns[iVar] * sizeof(double));
         }
       }
 
@@ -1017,9 +1020,9 @@ void InitializeUpdate(BODY *body, CONTROL *control, MODULE *module,
               malloc(update[iBody].iNumCO2MassSol * sizeof(int));
         control->Evolve.tmpUpdate[iBody].iaBody[iVar] =
               malloc(update[iBody].iNumCO2MassSol * sizeof(int *));
-        for (iSubStep=0; iSubStep < 4; iSubStep++) {
-            control->Evolve.daDerivProc[iSubStep][iBody][iVar] =
-                  malloc(update[iBody].iNumEqns[iVar] * sizeof(double));
+        for (iSubStep = 0; iSubStep < 4; iSubStep++) {
+          control->Evolve.daDerivProc[iSubStep][iBody][iVar] =
+                malloc(update[iBody].iNumEqns[iVar] * sizeof(double));
         }
       }
 
@@ -1031,7 +1034,7 @@ void InitializeUpdate(BODY *body, CONTROL *control, MODULE *module,
 
       (*fnUpdate)[iBody][iVar]        = malloc(iEqn * sizeof(fnUpdateVariable));
       update[iBody].daDerivProc[iVar] = malloc(iEqn * sizeof(double));
-       iVar++;
+      iVar++;
     }
     /* HERE
      * End of MagmOc variables
@@ -1066,9 +1069,9 @@ void InitializeUpdate(BODY *body, CONTROL *control, MODULE *module,
               malloc(update[iBody].iNum26AlCore * sizeof(int));
         control->Evolve.tmpUpdate[iBody].iaBody[iVar] =
               malloc(update[iBody].iNum26AlCore * sizeof(int *));
-        for (iSubStep=0; iSubStep < 4; iSubStep++) {
-            control->Evolve.daDerivProc[iSubStep][iBody][iVar] =
-                  malloc(update[iBody].iNumEqns[iVar] * sizeof(double));
+        for (iSubStep = 0; iSubStep < 4; iSubStep++) {
+          control->Evolve.daDerivProc[iSubStep][iBody][iVar] =
+                malloc(update[iBody].iNumEqns[iVar] * sizeof(double));
         }
       }
 
@@ -1112,9 +1115,9 @@ void InitializeUpdate(BODY *body, CONTROL *control, MODULE *module,
               malloc(update[iBody].iNum26AlMan * sizeof(int));
         control->Evolve.tmpUpdate[iBody].iaBody[iVar] =
               malloc(update[iBody].iNum26AlMan * sizeof(int *));
-        for (iSubStep=0; iSubStep < 4; iSubStep++) {
-            control->Evolve.daDerivProc[iSubStep][iBody][iVar] =
-                  malloc(update[iBody].iNumEqns[iVar] * sizeof(double));
+        for (iSubStep = 0; iSubStep < 4; iSubStep++) {
+          control->Evolve.daDerivProc[iSubStep][iBody][iVar] =
+                malloc(update[iBody].iNumEqns[iVar] * sizeof(double));
         }
       }
 
@@ -1158,9 +1161,9 @@ void InitializeUpdate(BODY *body, CONTROL *control, MODULE *module,
               malloc(update[iBody].iNum40KCore * sizeof(int));
         control->Evolve.tmpUpdate[iBody].iaBody[iVar] =
               malloc(update[iBody].iNum40KCore * sizeof(int *));
-        for (iSubStep=0; iSubStep < 4; iSubStep++) {
-            control->Evolve.daDerivProc[iSubStep][iBody][iVar] =
-                  malloc(update[iBody].iNumEqns[iVar] * sizeof(double));
+        for (iSubStep = 0; iSubStep < 4; iSubStep++) {
+          control->Evolve.daDerivProc[iSubStep][iBody][iVar] =
+                malloc(update[iBody].iNumEqns[iVar] * sizeof(double));
         }
       }
 
@@ -1204,9 +1207,9 @@ void InitializeUpdate(BODY *body, CONTROL *control, MODULE *module,
               malloc(update[iBody].iNum40KMan * sizeof(int));
         control->Evolve.tmpUpdate[iBody].iaBody[iVar] =
               malloc(update[iBody].iNum40KMan * sizeof(int *));
-        for (iSubStep=0; iSubStep < 4; iSubStep++) {
-            control->Evolve.daDerivProc[iSubStep][iBody][iVar] =
-                  malloc(update[iBody].iNumEqns[iVar] * sizeof(double));
+        for (iSubStep = 0; iSubStep < 4; iSubStep++) {
+          control->Evolve.daDerivProc[iSubStep][iBody][iVar] =
+                malloc(update[iBody].iNumEqns[iVar] * sizeof(double));
         }
       }
 
@@ -1250,9 +1253,9 @@ void InitializeUpdate(BODY *body, CONTROL *control, MODULE *module,
               malloc(update[iBody].iNum40KCrust * sizeof(int));
         control->Evolve.tmpUpdate[iBody].iaBody[iVar] =
               malloc(update[iBody].iNum40KCrust * sizeof(int *));
-        for (iSubStep=0; iSubStep < 4; iSubStep++) {
-            control->Evolve.daDerivProc[iSubStep][iBody][iVar] =
-                  malloc(update[iBody].iNumEqns[iVar] * sizeof(double));
+        for (iSubStep = 0; iSubStep < 4; iSubStep++) {
+          control->Evolve.daDerivProc[iSubStep][iBody][iVar] =
+                malloc(update[iBody].iNumEqns[iVar] * sizeof(double));
         }
       }
 
@@ -1296,9 +1299,9 @@ void InitializeUpdate(BODY *body, CONTROL *control, MODULE *module,
               malloc(update[iBody].iNum232ThCore * sizeof(int));
         control->Evolve.tmpUpdate[iBody].iaBody[iVar] =
               malloc(update[iBody].iNum232ThCore * sizeof(int *));
-        for (iSubStep=0; iSubStep < 4; iSubStep++) {
-            control->Evolve.daDerivProc[iSubStep][iBody][iVar] =
-                  malloc(update[iBody].iNumEqns[iVar] * sizeof(double));
+        for (iSubStep = 0; iSubStep < 4; iSubStep++) {
+          control->Evolve.daDerivProc[iSubStep][iBody][iVar] =
+                malloc(update[iBody].iNumEqns[iVar] * sizeof(double));
         }
       }
 
@@ -1342,9 +1345,9 @@ void InitializeUpdate(BODY *body, CONTROL *control, MODULE *module,
               malloc(update[iBody].iNum232ThMan * sizeof(int));
         control->Evolve.tmpUpdate[iBody].iaBody[iVar] =
               malloc(update[iBody].iNum232ThMan * sizeof(int *));
-        for (iSubStep=0; iSubStep < 4; iSubStep++) {
-            control->Evolve.daDerivProc[iSubStep][iBody][iVar] =
-                  malloc(update[iBody].iNumEqns[iVar] * sizeof(double));
+        for (iSubStep = 0; iSubStep < 4; iSubStep++) {
+          control->Evolve.daDerivProc[iSubStep][iBody][iVar] =
+                malloc(update[iBody].iNumEqns[iVar] * sizeof(double));
         }
       }
 
@@ -1388,9 +1391,9 @@ void InitializeUpdate(BODY *body, CONTROL *control, MODULE *module,
               malloc(update[iBody].iNum232ThCrust * sizeof(int));
         control->Evolve.tmpUpdate[iBody].iaBody[iVar] =
               malloc(update[iBody].iNum232ThCrust * sizeof(int *));
-        for (iSubStep=0; iSubStep < 4; iSubStep++) {
-            control->Evolve.daDerivProc[iSubStep][iBody][iVar] =
-                  malloc(update[iBody].iNumEqns[iVar] * sizeof(double));
+        for (iSubStep = 0; iSubStep < 4; iSubStep++) {
+          control->Evolve.daDerivProc[iSubStep][iBody][iVar] =
+                malloc(update[iBody].iNumEqns[iVar] * sizeof(double));
         }
       }
 
@@ -1434,9 +1437,9 @@ void InitializeUpdate(BODY *body, CONTROL *control, MODULE *module,
               malloc(update[iBody].iNum235UCore * sizeof(int));
         control->Evolve.tmpUpdate[iBody].iaBody[iVar] =
               malloc(update[iBody].iNum235UCore * sizeof(int *));
-        for (iSubStep=0; iSubStep < 4; iSubStep++) {
-            control->Evolve.daDerivProc[iSubStep][iBody][iVar] =
-                  malloc(update[iBody].iNumEqns[iVar] * sizeof(double));
+        for (iSubStep = 0; iSubStep < 4; iSubStep++) {
+          control->Evolve.daDerivProc[iSubStep][iBody][iVar] =
+                malloc(update[iBody].iNumEqns[iVar] * sizeof(double));
         }
       }
 
@@ -1480,9 +1483,9 @@ void InitializeUpdate(BODY *body, CONTROL *control, MODULE *module,
               malloc(update[iBody].iNum235UMan * sizeof(int));
         control->Evolve.tmpUpdate[iBody].iaBody[iVar] =
               malloc(update[iBody].iNum235UMan * sizeof(int *));
-        for (iSubStep=0; iSubStep < 4; iSubStep++) {
-            control->Evolve.daDerivProc[iSubStep][iBody][iVar] =
-                  malloc(update[iBody].iNumEqns[iVar] * sizeof(double));
+        for (iSubStep = 0; iSubStep < 4; iSubStep++) {
+          control->Evolve.daDerivProc[iSubStep][iBody][iVar] =
+                malloc(update[iBody].iNumEqns[iVar] * sizeof(double));
         }
       }
 
@@ -1526,9 +1529,9 @@ void InitializeUpdate(BODY *body, CONTROL *control, MODULE *module,
               malloc(update[iBody].iNum235UCrust * sizeof(int));
         control->Evolve.tmpUpdate[iBody].iaBody[iVar] =
               malloc(update[iBody].iNum235UCrust * sizeof(int *));
-        for (iSubStep=0; iSubStep < 4; iSubStep++) {
-            control->Evolve.daDerivProc[iSubStep][iBody][iVar] =
-                  malloc(update[iBody].iNumEqns[iVar] * sizeof(double));
+        for (iSubStep = 0; iSubStep < 4; iSubStep++) {
+          control->Evolve.daDerivProc[iSubStep][iBody][iVar] =
+                malloc(update[iBody].iNumEqns[iVar] * sizeof(double));
         }
       }
 
@@ -1572,9 +1575,9 @@ void InitializeUpdate(BODY *body, CONTROL *control, MODULE *module,
               malloc(update[iBody].iNum238UCore * sizeof(int));
         control->Evolve.tmpUpdate[iBody].iaBody[iVar] =
               malloc(update[iBody].iNum238UCore * sizeof(int *));
-        for (iSubStep=0; iSubStep < 4; iSubStep++) {
-            control->Evolve.daDerivProc[iSubStep][iBody][iVar] =
-                  malloc(update[iBody].iNumEqns[iVar] * sizeof(double));
+        for (iSubStep = 0; iSubStep < 4; iSubStep++) {
+          control->Evolve.daDerivProc[iSubStep][iBody][iVar] =
+                malloc(update[iBody].iNumEqns[iVar] * sizeof(double));
         }
       }
 
@@ -1618,9 +1621,9 @@ void InitializeUpdate(BODY *body, CONTROL *control, MODULE *module,
               malloc(update[iBody].iNum238UMan * sizeof(int));
         control->Evolve.tmpUpdate[iBody].iaBody[iVar] =
               malloc(update[iBody].iNum238UMan * sizeof(int *));
-        for (iSubStep=0; iSubStep < 4; iSubStep++) {
-            control->Evolve.daDerivProc[iSubStep][iBody][iVar] =
-                  malloc(update[iBody].iNumEqns[iVar] * sizeof(double));
+        for (iSubStep = 0; iSubStep < 4; iSubStep++) {
+          control->Evolve.daDerivProc[iSubStep][iBody][iVar] =
+                malloc(update[iBody].iNumEqns[iVar] * sizeof(double));
         }
       }
 
@@ -1664,9 +1667,9 @@ void InitializeUpdate(BODY *body, CONTROL *control, MODULE *module,
               malloc(update[iBody].iNum238UCrust * sizeof(int));
         control->Evolve.tmpUpdate[iBody].iaBody[iVar] =
               malloc(update[iBody].iNum238UCrust * sizeof(int *));
-        for (iSubStep=0; iSubStep < 4; iSubStep++) {
-            control->Evolve.daDerivProc[iSubStep][iBody][iVar] =
-                  malloc(update[iBody].iNumEqns[iVar] * sizeof(double));
+        for (iSubStep = 0; iSubStep < 4; iSubStep++) {
+          control->Evolve.daDerivProc[iSubStep][iBody][iVar] =
+                malloc(update[iBody].iNumEqns[iVar] * sizeof(double));
         }
       }
 
@@ -1710,9 +1713,9 @@ void InitializeUpdate(BODY *body, CONTROL *control, MODULE *module,
               malloc(update[iBody].iNumEnvelopeMass * sizeof(int));
         control->Evolve.tmpUpdate[iBody].iaBody[iVar] =
               malloc(update[iBody].iNumEnvelopeMass * sizeof(int *));
-        for (iSubStep=0; iSubStep < 4; iSubStep++) {
-            control->Evolve.daDerivProc[iSubStep][iBody][iVar] =
-                  malloc(update[iBody].iNumEqns[iVar] * sizeof(double));
+        for (iSubStep = 0; iSubStep < 4; iSubStep++) {
+          control->Evolve.daDerivProc[iSubStep][iBody][iVar] =
+                malloc(update[iBody].iNumEqns[iVar] * sizeof(double));
         }
       }
 
@@ -1762,9 +1765,9 @@ void InitializeUpdate(BODY *body, CONTROL *control, MODULE *module,
               malloc(update[iBody].iNumDynEllip * sizeof(int));
         control->Evolve.tmpUpdate[iBody].iaBody[iVar] =
               malloc(update[iBody].iNumDynEllip * sizeof(int *));
-        for (iSubStep=0; iSubStep < 4; iSubStep++) {
-            control->Evolve.daDerivProc[iSubStep][iBody][iVar] =
-                  malloc(update[iBody].iNumEqns[iVar] * sizeof(double));
+        for (iSubStep = 0; iSubStep < 4; iSubStep++) {
+          control->Evolve.daDerivProc[iSubStep][iBody][iVar] =
+                malloc(update[iBody].iNumEqns[iVar] * sizeof(double));
         }
       }
 
@@ -1791,7 +1794,8 @@ void InitializeUpdate(BODY *body, CONTROL *control, MODULE *module,
       update[iBody].iaBody[iVar] =
             malloc(update[iBody].iNumHecc * sizeof(int *));
       update[iBody].iaType[iVar] = malloc(update[iBody].iNumHecc * sizeof(int));
-      // XXX +1 allows for GR correction -- better to set iNumKecc based on user input!
+      // XXX +1 allows for GR correction -- better to set iNumKecc based on user
+      // input!
       update[iBody].iaModule[iVar] =
             malloc(update[iBody].iNumHecc * sizeof(int));
 
@@ -1808,9 +1812,9 @@ void InitializeUpdate(BODY *body, CONTROL *control, MODULE *module,
               malloc(update[iBody].iNumHecc * sizeof(int));
         control->Evolve.tmpUpdate[iBody].iaBody[iVar] =
               malloc(update[iBody].iNumHecc * sizeof(int *));
-        for (iSubStep=0; iSubStep < 4; iSubStep++) {
-            control->Evolve.daDerivProc[iSubStep][iBody][iVar] =
-                  malloc(update[iBody].iNumEqns[iVar] * sizeof(double));
+        for (iSubStep = 0; iSubStep < 4; iSubStep++) {
+          control->Evolve.daDerivProc[iSubStep][iBody][iVar] =
+                malloc(update[iBody].iNumEqns[iVar] * sizeof(double));
         }
       }
 
@@ -1837,7 +1841,8 @@ void InitializeUpdate(BODY *body, CONTROL *control, MODULE *module,
       update[iBody].iaBody[iVar] =
             malloc(update[iBody].iNumKecc * sizeof(int *));
       update[iBody].iaType[iVar] = malloc(update[iBody].iNumKecc * sizeof(int));
-      // XXX +1 allows for GR correction -- better to set iNumKecc based on user input!
+      // XXX +1 allows for GR correction -- better to set iNumKecc based on user
+      // input!
       update[iBody].iaModule[iVar] =
             malloc(update[iBody].iNumKecc * sizeof(int));
 
@@ -1854,9 +1859,9 @@ void InitializeUpdate(BODY *body, CONTROL *control, MODULE *module,
               malloc(update[iBody].iNumKecc * sizeof(int));
         control->Evolve.tmpUpdate[iBody].iaBody[iVar] =
               malloc(update[iBody].iNumKecc * sizeof(int *));
-        for (iSubStep=0; iSubStep < 4; iSubStep++) {
-            control->Evolve.daDerivProc[iSubStep][iBody][iVar] =
-                  malloc(update[iBody].iNumEqns[iVar] * sizeof(double));
+        for (iSubStep = 0; iSubStep < 4; iSubStep++) {
+          control->Evolve.daDerivProc[iSubStep][iBody][iVar] =
+                malloc(update[iBody].iNumEqns[iVar] * sizeof(double));
         }
       }
 
@@ -1874,10 +1879,10 @@ void InitializeUpdate(BODY *body, CONTROL *control, MODULE *module,
     // Luminosity:
     /* Someday we'll fix this, here's a template
     AddPrimaryVariable(body,control,module,update,fnUpdate,&body[iBody].dLuminosity,&control->Evolve.tmpBody[iBody].dLuminosity,
-                       &update[iBody].iLuminosity,&module->fnFinalizeUpdate[iBody][iModule][iVar | VLUMINOSITY?],
-                       update[iBody].iNumLuminosity,iBody,iVar,VLUMINOSITY);
+                       &update[iBody].iLuminosity,&module->fnFinalizeUpdate[iBody][iModule][iVar
+    | VLUMINOSITY?], update[iBody].iNumLuminosity,iBody,iVar,VLUMINOSITY);
     */
-    
+
     update[iBody].iLuminosity = -1;
     if (update[iBody].iNumLuminosity) {
       update[iBody].iLuminosity    = iVar;
@@ -1906,9 +1911,9 @@ void InitializeUpdate(BODY *body, CONTROL *control, MODULE *module,
               malloc(update[iBody].iNumLuminosity * sizeof(int));
         control->Evolve.tmpUpdate[iBody].iaBody[iVar] =
               malloc(update[iBody].iNumLuminosity * sizeof(int *));
-        for (iSubStep=0; iSubStep < 4; iSubStep++) {
-            control->Evolve.daDerivProc[iSubStep][iBody][iVar] =
-                  malloc(update[iBody].iNumEqns[iVar] * sizeof(double));
+        for (iSubStep = 0; iSubStep < 4; iSubStep++) {
+          control->Evolve.daDerivProc[iSubStep][iBody][iVar] =
+                malloc(update[iBody].iNumEqns[iVar] * sizeof(double));
         }
       }
 
@@ -1990,9 +1995,9 @@ void InitializeUpdate(BODY *body, CONTROL *control, MODULE *module,
               malloc(update[iBody].iNumPinc * sizeof(int));
         control->Evolve.tmpUpdate[iBody].iaBody[iVar] =
               malloc(update[iBody].iNumPinc * sizeof(int *));
-        for (iSubStep=0; iSubStep < 4; iSubStep++) {
-            control->Evolve.daDerivProc[iSubStep][iBody][iVar] =
-                  malloc(update[iBody].iNumEqns[iVar] * sizeof(double));
+        for (iSubStep = 0; iSubStep < 4; iSubStep++) {
+          control->Evolve.daDerivProc[iSubStep][iBody][iVar] =
+                malloc(update[iBody].iNumEqns[iVar] * sizeof(double));
         }
       }
 
@@ -2035,9 +2040,9 @@ void InitializeUpdate(BODY *body, CONTROL *control, MODULE *module,
               malloc(update[iBody].iNumQinc * sizeof(int));
         control->Evolve.tmpUpdate[iBody].iaBody[iVar] =
               malloc(update[iBody].iNumQinc * sizeof(int *));
-        for (iSubStep=0; iSubStep < 4; iSubStep++) {
-            control->Evolve.daDerivProc[iSubStep][iBody][iVar] =
-                  malloc(update[iBody].iNumEqns[iVar] * sizeof(double));
+        for (iSubStep = 0; iSubStep < 4; iSubStep++) {
+          control->Evolve.daDerivProc[iSubStep][iBody][iVar] =
+                malloc(update[iBody].iNumEqns[iVar] * sizeof(double));
         }
       }
 
@@ -2081,9 +2086,9 @@ void InitializeUpdate(BODY *body, CONTROL *control, MODULE *module,
               malloc(update[iBody].iNumRadius * sizeof(int));
         control->Evolve.tmpUpdate[iBody].iaBody[iVar] =
               malloc(update[iBody].iNumRadius * sizeof(int *));
-      for (iSubStep=0; iSubStep < 4; iSubStep++) {
-            control->Evolve.daDerivProc[iSubStep][iBody][iVar] =
-                  malloc(update[iBody].iNumEqns[iVar] * sizeof(double));
+        for (iSubStep = 0; iSubStep < 4; iSubStep++) {
+          control->Evolve.daDerivProc[iSubStep][iBody][iVar] =
+                malloc(update[iBody].iNumEqns[iVar] * sizeof(double));
         }
       }
 
@@ -2126,9 +2131,9 @@ void InitializeUpdate(BODY *body, CONTROL *control, MODULE *module,
               malloc(update[iBody].iNumMass * sizeof(int));
         control->Evolve.tmpUpdate[iBody].iaBody[iVar] =
               malloc(update[iBody].iNumMass * sizeof(int *));
-        for (iSubStep=0; iSubStep < 4; iSubStep++) {
-            control->Evolve.daDerivProc[iSubStep][iBody][iVar] =
-                  malloc(update[iBody].iNumEqns[iVar] * sizeof(double));
+        for (iSubStep = 0; iSubStep < 4; iSubStep++) {
+          control->Evolve.daDerivProc[iSubStep][iBody][iVar] =
+                malloc(update[iBody].iNumEqns[iVar] * sizeof(double));
         }
       }
 
@@ -2171,9 +2176,9 @@ void InitializeUpdate(BODY *body, CONTROL *control, MODULE *module,
               malloc(update[iBody].iNumRot * sizeof(int));
         control->Evolve.tmpUpdate[iBody].iaBody[iVar] =
               malloc(update[iBody].iNumRot * sizeof(int *));
-        for (iSubStep=0; iSubStep < 4; iSubStep++) {
-            control->Evolve.daDerivProc[iSubStep][iBody][iVar] =
-                  malloc(update[iBody].iNumEqns[iVar] * sizeof(double));
+        for (iSubStep = 0; iSubStep < 4; iSubStep++) {
+          control->Evolve.daDerivProc[iSubStep][iBody][iVar] =
+                malloc(update[iBody].iNumEqns[iVar] * sizeof(double));
         }
       }
 
@@ -2216,9 +2221,9 @@ void InitializeUpdate(BODY *body, CONTROL *control, MODULE *module,
               malloc(update[iBody].iNumSemi * sizeof(int));
         control->Evolve.tmpUpdate[iBody].iaBody[iVar] =
               malloc(update[iBody].iNumSemi * sizeof(int *));
-        for (iSubStep=0; iSubStep < 4; iSubStep++) {
-            control->Evolve.daDerivProc[iSubStep][iBody][iVar] =
-                  malloc(update[iBody].iNumEqns[iVar] * sizeof(double));
+        for (iSubStep = 0; iSubStep < 4; iSubStep++) {
+          control->Evolve.daDerivProc[iSubStep][iBody][iVar] =
+                malloc(update[iBody].iNumEqns[iVar] * sizeof(double));
         }
       }
 
@@ -2266,9 +2271,9 @@ void InitializeUpdate(BODY *body, CONTROL *control, MODULE *module,
               malloc(update[iBody].iNumSurfaceWaterMass * sizeof(int));
         control->Evolve.tmpUpdate[iBody].iaBody[iVar] =
               malloc(update[iBody].iNumSurfaceWaterMass * sizeof(int *));
-        for (iSubStep=0; iSubStep < 4; iSubStep++) {
-            control->Evolve.daDerivProc[iSubStep][iBody][iVar] =
-                  malloc(update[iBody].iNumEqns[iVar] * sizeof(double));
+        for (iSubStep = 0; iSubStep < 4; iSubStep++) {
+          control->Evolve.daDerivProc[iSubStep][iBody][iVar] =
+                malloc(update[iBody].iNumEqns[iVar] * sizeof(double));
         }
       }
 
@@ -2312,10 +2317,10 @@ void InitializeUpdate(BODY *body, CONTROL *control, MODULE *module,
               malloc(update[iBody].iNumOxygenMass * sizeof(int));
         control->Evolve.tmpUpdate[iBody].iaBody[iVar] =
               malloc(update[iBody].iNumOxygenMass * sizeof(int *));
-        for (iSubStep=0; iSubStep < 4; iSubStep++) {
-            control->Evolve.daDerivProc[iSubStep][iBody][iVar] =
-                  malloc(update[iBody].iNumEqns[iVar] * sizeof(double));
-        }              
+        for (iSubStep = 0; iSubStep < 4; iSubStep++) {
+          control->Evolve.daDerivProc[iSubStep][iBody][iVar] =
+                malloc(update[iBody].iNumEqns[iVar] * sizeof(double));
+        }
       }
 
       iEqn = 0;
@@ -2358,9 +2363,9 @@ void InitializeUpdate(BODY *body, CONTROL *control, MODULE *module,
               malloc(update[iBody].iNumOxygenMantleMass * sizeof(int));
         control->Evolve.tmpUpdate[iBody].iaBody[iVar] =
               malloc(update[iBody].iNumOxygenMantleMass * sizeof(int *));
-        for (iSubStep=0; iSubStep < 4; iSubStep++) {
-            control->Evolve.daDerivProc[iSubStep][iBody][iVar] =
-                  malloc(update[iBody].iNumEqns[iVar] * sizeof(double));
+        for (iSubStep = 0; iSubStep < 4; iSubStep++) {
+          control->Evolve.daDerivProc[iSubStep][iBody][iVar] =
+                malloc(update[iBody].iNumEqns[iVar] * sizeof(double));
         }
       }
 
@@ -2404,9 +2409,9 @@ void InitializeUpdate(BODY *body, CONTROL *control, MODULE *module,
               malloc(update[iBody].iNumTemperature * sizeof(int));
         control->Evolve.tmpUpdate[iBody].iaBody[iVar] =
               malloc(update[iBody].iNumTemperature * sizeof(int *));
-        for (iSubStep=0; iSubStep < 4; iSubStep++) {
-            control->Evolve.daDerivProc[iSubStep][iBody][iVar] =
-                  malloc(update[iBody].iNumEqns[iVar] * sizeof(double));
+        for (iSubStep = 0; iSubStep < 4; iSubStep++) {
+          control->Evolve.daDerivProc[iSubStep][iBody][iVar] =
+                malloc(update[iBody].iNumEqns[iVar] * sizeof(double));
         }
       }
 
@@ -2450,9 +2455,9 @@ void InitializeUpdate(BODY *body, CONTROL *control, MODULE *module,
               malloc(update[iBody].iNumRadGyra * sizeof(int));
         control->Evolve.tmpUpdate[iBody].iaBody[iVar] =
               malloc(update[iBody].iNumRadGyra * sizeof(int *));
-        for (iSubStep=0; iSubStep < 4; iSubStep++) {
-            control->Evolve.daDerivProc[iSubStep][iBody][iVar] =
-                  malloc(update[iBody].iNumEqns[iVar] * sizeof(double));
+        for (iSubStep = 0; iSubStep < 4; iSubStep++) {
+          control->Evolve.daDerivProc[iSubStep][iBody][iVar] =
+                malloc(update[iBody].iNumEqns[iVar] * sizeof(double));
         }
       }
 
@@ -2496,9 +2501,9 @@ void InitializeUpdate(BODY *body, CONTROL *control, MODULE *module,
               malloc(update[iBody].iNumTCore * sizeof(int));
         control->Evolve.tmpUpdate[iBody].iaBody[iVar] =
               malloc(update[iBody].iNumTCore * sizeof(int *));
-        for (iSubStep=0; iSubStep < 4; iSubStep++) {
-            control->Evolve.daDerivProc[iSubStep][iBody][iVar] =
-                  malloc(update[iBody].iNumEqns[iVar] * sizeof(double));
+        for (iSubStep = 0; iSubStep < 4; iSubStep++) {
+          control->Evolve.daDerivProc[iSubStep][iBody][iVar] =
+                malloc(update[iBody].iNumEqns[iVar] * sizeof(double));
         }
       }
 
@@ -2541,9 +2546,9 @@ void InitializeUpdate(BODY *body, CONTROL *control, MODULE *module,
               malloc(update[iBody].iNumTMan * sizeof(int));
         control->Evolve.tmpUpdate[iBody].iaBody[iVar] =
               malloc(update[iBody].iNumTMan * sizeof(int *));
-        for (iSubStep=0; iSubStep < 4; iSubStep++) {
-            control->Evolve.daDerivProc[iSubStep][iBody][iVar] =
-                  malloc(update[iBody].iNumEqns[iVar] * sizeof(double));
+        for (iSubStep = 0; iSubStep < 4; iSubStep++) {
+          control->Evolve.daDerivProc[iSubStep][iBody][iVar] =
+                malloc(update[iBody].iNumEqns[iVar] * sizeof(double));
         }
       }
 
@@ -2587,9 +2592,9 @@ void InitializeUpdate(BODY *body, CONTROL *control, MODULE *module,
         control->Evolve.tmpUpdate[iBody].iaBody[iVar] =
               malloc(update[iBody].iNumXobl * sizeof(int *));
         control->Evolve.tmpUpdate[iBody].iXobl = iVar;
-        for (iSubStep=0; iSubStep < 4; iSubStep++) {
-            control->Evolve.daDerivProc[iSubStep][iBody][iVar] =
-                  malloc(update[iBody].iNumEqns[iVar] * sizeof(double));
+        for (iSubStep = 0; iSubStep < 4; iSubStep++) {
+          control->Evolve.daDerivProc[iSubStep][iBody][iVar] =
+                malloc(update[iBody].iNumEqns[iVar] * sizeof(double));
         }
       }
 
@@ -2633,9 +2638,9 @@ void InitializeUpdate(BODY *body, CONTROL *control, MODULE *module,
         control->Evolve.tmpUpdate[iBody].iaBody[iVar] =
               malloc(update[iBody].iNumYobl * sizeof(int *));
         control->Evolve.tmpUpdate[iBody].iYobl = iVar;
-        for (iSubStep=0; iSubStep < 4; iSubStep++) {
-            control->Evolve.daDerivProc[iSubStep][iBody][iVar] =
-                  malloc(update[iBody].iNumEqns[iVar] * sizeof(double));
+        for (iSubStep = 0; iSubStep < 4; iSubStep++) {
+          control->Evolve.daDerivProc[iSubStep][iBody][iVar] =
+                malloc(update[iBody].iNumEqns[iVar] * sizeof(double));
         }
       }
 
@@ -2679,9 +2684,9 @@ void InitializeUpdate(BODY *body, CONTROL *control, MODULE *module,
         control->Evolve.tmpUpdate[iBody].iaBody[iVar] =
               malloc(update[iBody].iNumZobl * sizeof(int *));
         control->Evolve.tmpUpdate[iBody].iZobl = iVar;
-        for (iSubStep=0; iSubStep < 4; iSubStep++) {
-            control->Evolve.daDerivProc[iSubStep][iBody][iVar] =
-                  malloc(update[iBody].iNumEqns[iVar] * sizeof(double));
+        for (iSubStep = 0; iSubStep < 4; iSubStep++) {
+          control->Evolve.daDerivProc[iSubStep][iBody][iVar] =
+                malloc(update[iBody].iNumEqns[iVar] * sizeof(double));
         }
       }
 
@@ -2726,9 +2731,9 @@ void InitializeUpdate(BODY *body, CONTROL *control, MODULE *module,
               malloc(update[iBody].iNumCBPR * sizeof(int));
         control->Evolve.tmpUpdate[iBody].iaBody[iVar] =
               malloc(update[iBody].iNumCBPR * sizeof(int *));
-        for (iSubStep=0; iSubStep < 4; iSubStep++) {
-            control->Evolve.daDerivProc[iSubStep][iBody][iVar] =
-                  malloc(update[iBody].iNumEqns[iVar] * sizeof(double));
+        for (iSubStep = 0; iSubStep < 4; iSubStep++) {
+          control->Evolve.daDerivProc[iSubStep][iBody][iVar] =
+                malloc(update[iBody].iNumEqns[iVar] * sizeof(double));
         }
       }
 
@@ -2772,9 +2777,9 @@ void InitializeUpdate(BODY *body, CONTROL *control, MODULE *module,
               malloc(update[iBody].iNumCBPZ * sizeof(int));
         control->Evolve.tmpUpdate[iBody].iaBody[iVar] =
               malloc(update[iBody].iNumCBPZ * sizeof(int *));
-        for (iSubStep=0; iSubStep < 4; iSubStep++) {
-            control->Evolve.daDerivProc[iSubStep][iBody][iVar] =
-                  malloc(update[iBody].iNumEqns[iVar] * sizeof(double));
+        for (iSubStep = 0; iSubStep < 4; iSubStep++) {
+          control->Evolve.daDerivProc[iSubStep][iBody][iVar] =
+                malloc(update[iBody].iNumEqns[iVar] * sizeof(double));
         }
       }
 
@@ -2819,9 +2824,9 @@ void InitializeUpdate(BODY *body, CONTROL *control, MODULE *module,
               malloc(update[iBody].iNumCBPPhi * sizeof(int));
         control->Evolve.tmpUpdate[iBody].iaBody[iVar] =
               malloc(update[iBody].iNumCBPPhi * sizeof(int *));
-        for (iSubStep=0; iSubStep < 4; iSubStep++) {
-            control->Evolve.daDerivProc[iSubStep][iBody][iVar] =
-                  malloc(update[iBody].iNumEqns[iVar] * sizeof(double));
+        for (iSubStep = 0; iSubStep < 4; iSubStep++) {
+          control->Evolve.daDerivProc[iSubStep][iBody][iVar] =
+                malloc(update[iBody].iNumEqns[iVar] * sizeof(double));
         }
       }
 
@@ -2866,9 +2871,9 @@ void InitializeUpdate(BODY *body, CONTROL *control, MODULE *module,
               malloc(update[iBody].iNumCBPRDot * sizeof(int));
         control->Evolve.tmpUpdate[iBody].iaBody[iVar] =
               malloc(update[iBody].iNumCBPRDot * sizeof(int *));
-        for (iSubStep=0; iSubStep < 4; iSubStep++) {
-            control->Evolve.daDerivProc[iSubStep][iBody][iVar] =
-                  malloc(update[iBody].iNumEqns[iVar] * sizeof(double));
+        for (iSubStep = 0; iSubStep < 4; iSubStep++) {
+          control->Evolve.daDerivProc[iSubStep][iBody][iVar] =
+                malloc(update[iBody].iNumEqns[iVar] * sizeof(double));
         }
       }
 
@@ -2913,9 +2918,9 @@ void InitializeUpdate(BODY *body, CONTROL *control, MODULE *module,
               malloc(update[iBody].iNumCBPZDot * sizeof(int));
         control->Evolve.tmpUpdate[iBody].iaBody[iVar] =
               malloc(update[iBody].iNumCBPZDot * sizeof(int *));
-        for (iSubStep=0; iSubStep < 4; iSubStep++) {
-            control->Evolve.daDerivProc[iSubStep][iBody][iVar] =
-                  malloc(update[iBody].iNumEqns[iVar] * sizeof(double));
+        for (iSubStep = 0; iSubStep < 4; iSubStep++) {
+          control->Evolve.daDerivProc[iSubStep][iBody][iVar] =
+                malloc(update[iBody].iNumEqns[iVar] * sizeof(double));
         }
       }
 
@@ -2960,9 +2965,9 @@ void InitializeUpdate(BODY *body, CONTROL *control, MODULE *module,
               malloc(update[iBody].iNumCBPPhiDot * sizeof(int));
         control->Evolve.tmpUpdate[iBody].iaBody[iVar] =
               malloc(update[iBody].iNumCBPPhiDot * sizeof(int *));
-        for (iSubStep=0; iSubStep < 4; iSubStep++) {
-            control->Evolve.daDerivProc[iSubStep][iBody][iVar] =
-                  malloc(update[iBody].iNumEqns[iVar] * sizeof(double));
+        for (iSubStep = 0; iSubStep < 4; iSubStep++) {
+          control->Evolve.daDerivProc[iSubStep][iBody][iVar] =
+                malloc(update[iBody].iNumEqns[iVar] * sizeof(double));
         }
       }
 
@@ -3056,9 +3061,9 @@ void InitializeUpdate(BODY *body, CONTROL *control, MODULE *module,
               malloc(update[iBody].iNumEccX * sizeof(int));
         control->Evolve.tmpUpdate[iBody].iaBody[iVar] =
               malloc(update[iBody].iNumEccX * sizeof(int *));
-        for (iSubStep=0; iSubStep < 4; iSubStep++) {
-            control->Evolve.daDerivProc[iSubStep][iBody][iVar] =
-                  malloc(update[iBody].iNumEqns[iVar] * sizeof(double));
+        for (iSubStep = 0; iSubStep < 4; iSubStep++) {
+          control->Evolve.daDerivProc[iSubStep][iBody][iVar] =
+                malloc(update[iBody].iNumEqns[iVar] * sizeof(double));
         }
       }
 
@@ -3101,9 +3106,9 @@ void InitializeUpdate(BODY *body, CONTROL *control, MODULE *module,
               malloc(update[iBody].iNumEccY * sizeof(int));
         control->Evolve.tmpUpdate[iBody].iaBody[iVar] =
               malloc(update[iBody].iNumEccY * sizeof(int *));
-        for (iSubStep=0; iSubStep < 4; iSubStep++) {
-            control->Evolve.daDerivProc[iSubStep][iBody][iVar] =
-                  malloc(update[iBody].iNumEqns[iVar] * sizeof(double));
+        for (iSubStep = 0; iSubStep < 4; iSubStep++) {
+          control->Evolve.daDerivProc[iSubStep][iBody][iVar] =
+                malloc(update[iBody].iNumEqns[iVar] * sizeof(double));
         }
       }
 
@@ -3146,9 +3151,9 @@ void InitializeUpdate(BODY *body, CONTROL *control, MODULE *module,
               malloc(update[iBody].iNumEccZ * sizeof(int));
         control->Evolve.tmpUpdate[iBody].iaBody[iVar] =
               malloc(update[iBody].iNumEccZ * sizeof(int *));
-        for (iSubStep=0; iSubStep < 4; iSubStep++) {
-            control->Evolve.daDerivProc[iSubStep][iBody][iVar] =
-                  malloc(update[iBody].iNumEqns[iVar] * sizeof(double));
+        for (iSubStep = 0; iSubStep < 4; iSubStep++) {
+          control->Evolve.daDerivProc[iSubStep][iBody][iVar] =
+                malloc(update[iBody].iNumEqns[iVar] * sizeof(double));
         }
       }
 
@@ -3192,9 +3197,9 @@ void InitializeUpdate(BODY *body, CONTROL *control, MODULE *module,
               malloc(update[iBody].iNumAngMX * sizeof(int));
         control->Evolve.tmpUpdate[iBody].iaBody[iVar] =
               malloc(update[iBody].iNumAngMX * sizeof(int *));
-        for (iSubStep=0; iSubStep < 4; iSubStep++) {
-            control->Evolve.daDerivProc[iSubStep][iBody][iVar] =
-                  malloc(update[iBody].iNumEqns[iVar] * sizeof(double));
+        for (iSubStep = 0; iSubStep < 4; iSubStep++) {
+          control->Evolve.daDerivProc[iSubStep][iBody][iVar] =
+                malloc(update[iBody].iNumEqns[iVar] * sizeof(double));
         }
       }
 
@@ -3238,9 +3243,9 @@ void InitializeUpdate(BODY *body, CONTROL *control, MODULE *module,
               malloc(update[iBody].iNumAngMY * sizeof(int));
         control->Evolve.tmpUpdate[iBody].iaBody[iVar] =
               malloc(update[iBody].iNumAngMY * sizeof(int *));
-        for (iSubStep=0; iSubStep < 4; iSubStep++) {
-            control->Evolve.daDerivProc[iSubStep][iBody][iVar] =
-                  malloc(update[iBody].iNumEqns[iVar] * sizeof(double));
+        for (iSubStep = 0; iSubStep < 4; iSubStep++) {
+          control->Evolve.daDerivProc[iSubStep][iBody][iVar] =
+                malloc(update[iBody].iNumEqns[iVar] * sizeof(double));
         }
       }
 
@@ -3284,9 +3289,9 @@ void InitializeUpdate(BODY *body, CONTROL *control, MODULE *module,
               malloc(update[iBody].iNumAngMZ * sizeof(int));
         control->Evolve.tmpUpdate[iBody].iaBody[iVar] =
               malloc(update[iBody].iNumAngMZ * sizeof(int *));
-        for (iSubStep=0; iSubStep < 4; iSubStep++) {
-            control->Evolve.daDerivProc[iSubStep][iBody][iVar] =
-                  malloc(update[iBody].iNumEqns[iVar] * sizeof(double));
+        for (iSubStep = 0; iSubStep < 4; iSubStep++) {
+          control->Evolve.daDerivProc[iSubStep][iBody][iVar] =
+                malloc(update[iBody].iNumEqns[iVar] * sizeof(double));
         }
       }
 
@@ -3336,9 +3341,9 @@ void InitializeUpdate(BODY *body, CONTROL *control, MODULE *module,
               malloc(update[iBody].iNumLXUV * sizeof(int));
         control->Evolve.tmpUpdate[iBody].iaBody[iVar] =
               malloc(update[iBody].iNumLXUV * sizeof(int *));
-        for (iSubStep=0; iSubStep < 4; iSubStep++) {
-            control->Evolve.daDerivProc[iSubStep][iBody][iVar] =
-                  malloc(update[iBody].iNumEqns[iVar] * sizeof(double));
+        for (iSubStep = 0; iSubStep < 4; iSubStep++) {
+          control->Evolve.daDerivProc[iSubStep][iBody][iVar] =
+                malloc(update[iBody].iNumEqns[iVar] * sizeof(double));
         }
       }
 
@@ -3382,9 +3387,9 @@ void InitializeUpdate(BODY *body, CONTROL *control, MODULE *module,
               malloc(update[iBody].iNumLostAngMom * sizeof(int));
         control->Evolve.tmpUpdate[iBody].iaBody[iVar] =
               malloc(update[iBody].iNumLostAngMom * sizeof(int *));
-        for (iSubStep=0; iSubStep < 4; iSubStep++) {
-            control->Evolve.daDerivProc[iSubStep][iBody][iVar] =
-                  malloc(update[iBody].iNumEqns[iVar] * sizeof(double));
+        for (iSubStep = 0; iSubStep < 4; iSubStep++) {
+          control->Evolve.daDerivProc[iSubStep][iBody][iVar] =
+                malloc(update[iBody].iNumEqns[iVar] * sizeof(double));
         }
       }
 
@@ -3428,9 +3433,9 @@ void InitializeUpdate(BODY *body, CONTROL *control, MODULE *module,
               malloc(update[iBody].iNumLostEng * sizeof(int));
         control->Evolve.tmpUpdate[iBody].iaBody[iVar] =
               malloc(update[iBody].iNumLostEng * sizeof(int *));
-        for (iSubStep=0; iSubStep < 4; iSubStep++) {
-            control->Evolve.daDerivProc[iSubStep][iBody][iVar] =
-                  malloc(update[iBody].iNumEqns[iVar] * sizeof(double));
+        for (iSubStep = 0; iSubStep < 4; iSubStep++) {
+          control->Evolve.daDerivProc[iSubStep][iBody][iVar] =
+                malloc(update[iBody].iNumEqns[iVar] * sizeof(double));
         }
       }
 
diff --git a/src/verify.c b/src/verify.c
index d480dbb1..29b76192 100644
--- a/src/verify.c
+++ b/src/verify.c
@@ -354,17 +354,20 @@ void IntegrationWarning(char cName1[], char cName2[], char cName3[],
   /* Backward file name */
 }
 
-void VerifyIntegration(BODY *body, CONTROL *control, FILES *files,
-                       OPTIONS *options, SYSTEM *system,
-                       fnIntegrate *fnOneStep) {
-  int iFile, iFile1 = 0, iFile2 = 0;
-  char *cTmp=NULL;
+void AssignIntegrationDirection(CONTROL *control) {
 
+  if (control->Evolve.bDoBackward) {
+    control->Evolve.iDir = BACKWARD_INTEGRATION;
+  } else if (control->Evolve.bDoForward) {
+    control->Evolve.iDir = FORWARD_INTEGRATION;
+  } else {
+    control->Evolve.iDir = NO_INTEGRATION;
+  }
+}
 
-  // Initialize iDir to 0, i.e. assume no integrations requested to start
-  control->Evolve.iDir = 0;
+void VerifyOneDirectionOnly(CONTROL *control, FILES *files, OPTIONS *options) {
+  int iFile, iFile1, iFile2;
 
-  /* Were both Forward and Backward Set? */
   if (control->Evolve.bDoBackward && control->Evolve.bDoForward) {
     fprintf(stderr, "ERROR: Both %s and %s set. Only one is allowed.\n",
             options[OPT_BACK].cName, options[OPT_FORW].cName);
@@ -381,38 +384,55 @@ void VerifyIntegration(BODY *body, CONTROL *control, FILES *files,
           options[OPT_BACK].cFile[iFile1], options[OPT_FORW].cFile[iFile2],
           options[OPT_BACK].iLine[iFile1], options[OPT_FORW].iLine[iFile2]);
   }
+}
+
+void VerifyIntegrationDirection(CONTROL *control, FILES *files,
+                                OPTIONS *options) {
+
+  VerifyOneDirectionOnly(control, files, options);
+  AssignIntegrationDirection(control);
+}
+
+void SetBackwardFile(BODY *body, CONTROL *control, FILES *files,
+                     OPTIONS *options, SYSTEM *system) {
+  int iFile;
 
-  /* Fix backward output file */
   if (control->Evolve.bDoBackward) {
     for (iFile = 1; iFile < files->iNumInputs; iFile++) {
       if (options[OPT_OUTFILE].iLine[iFile] == -1) {
-        fvFormattedString(&files->Outfile[iFile - 1].cOut, "%s.%s.backward", system->cName,
-                body[iFile - 1].cName);
+        fvFormattedString(&files->Outfile[iFile - 1].cOut, "%s.%s.backward",
+                          system->cName, body[iFile - 1].cName);
         if (control->Io.iVerbose >= VERBINPUT) {
           fprintf(stderr, "INFO: %s not set, defaulting to %s.\n",
                   options[OPT_OUTFILE].cName, files->Outfile[iFile - 1].cOut);
         }
       }
     }
-    control->Evolve.iDir = -1;
   }
+}
+
+void SetForwardFile(BODY *body, CONTROL *control, FILES *files,
+                    OPTIONS *options, SYSTEM *system) {
+  int iFile;
 
-  /* Fix forward output file */
   if (control->Evolve.bDoForward) {
     for (iFile = 1; iFile < files->iNumInputs; iFile++) {
       if (options[OPT_OUTFILE].iLine[iFile] == -1) {
-        fvFormattedString(&files->Outfile[iFile - 1].cOut, "%s.%s.forward", system->cName,
-                body[iFile - 1].cName);
+        fvFormattedString(&files->Outfile[iFile - 1].cOut, "%s.%s.forward",
+                          system->cName, body[iFile - 1].cName);
         if (control->Io.iVerbose >= VERBINPUT) {
           fprintf(stderr, "INFO: %s not set, defaulting to %s.\n",
                   options[OPT_OUTFILE].cName, files->Outfile[iFile - 1].cOut);
         }
       }
     }
-    control->Evolve.iDir = 1;
   }
+}
+
+void CheckIntegrationFileExistence(CONTROL *control, FILES *files,
+                                   OPTIONS *options) {
+  int iFile;
 
-  /* Check for file existence */
   for (iFile = 0; iFile < files->iNumInputs - 1; iFile++) {
     if (bFileExists(files->Outfile[iFile].cOut)) {
       if (!control->Io.bOverwrite) {
@@ -424,8 +444,20 @@ void VerifyIntegration(BODY *body, CONTROL *control, FILES *files,
       unlink(files->Outfile[iFile].cOut);
     }
   }
+}
+
+void InitializeIntegrationFiles(BODY *body, CONTROL *control, FILES *files,
+                                OPTIONS *options, SYSTEM *system) {
+  SetBackwardFile(body, control, files, options, system);
+  SetForwardFile(body, control, files, options, system);
+  CheckIntegrationFileExistence(control, files, options);
+}
+
+
+void PrintIntegrationWarnings(CONTROL *control, FILES *files,
+                              OPTIONS *options) {
+  int iFile;
 
-  /* Was DoBackward or DoForward NOT set? */
   if (!control->Evolve.bDoBackward && !control->Evolve.bDoForward) {
     for (iFile = 0; iFile < files->iNumInputs; iFile++) {
       if (options[OPT_ETA].iLine[iFile] > -1) {
@@ -471,13 +503,17 @@ void VerifyIntegration(BODY *body, CONTROL *control, FILES *files,
       }
     }
   }
+}
+
+void AssignIntegrationMethod(CONTROL *control, OPTIONS *options,
+                             fnIntegrate *fnOneStep) {
+  char *cTmp = NULL;
 
   if (control->Evolve.iOneStep == EULER) {
     *fnOneStep = &EulerStep;
   } else if (control->Evolve.iOneStep == RUNGEKUTTA) {
     *fnOneStep = &RungeKutta4Step;
   } else {
-    /* Assign Default */
     fvFormattedString(&cTmp, options[OPT_INTEGRATIONMETHOD].cDefault);
     if (control->Io.iVerbose >= VERBINPUT) {
       fprintf(stderr, "INFO: %s not set, defaulting to %s.\n",
@@ -493,8 +529,11 @@ void VerifyIntegration(BODY *body, CONTROL *control, FILES *files,
       *fnOneStep               = &RungeKutta4Step;
     }
   }
+}
+
+void VerifyOutputTime(CONTROL *control, FILES *files, OPTIONS *options) {
+  int iFile, iFile1, iFile2;
 
-  /* Make sure output interval is less than stop time */
   if (control->Evolve.dStopTime < control->Io.dOutputTime) {
     fprintf(stderr, "ERROR: %s < %s is not allowed.\n",
             options[OPT_STOPTIME].cName, options[OPT_OUTPUTTIME].cName);
@@ -514,6 +553,74 @@ void VerifyIntegration(BODY *body, CONTROL *control, FILES *files,
   }
 }
 
+void StopTimeFromAge(CONTROL *control, FILES *files, OPTIONS *options,
+                     SYSTEM *system, int iFileStopAge) {
+  int iFile, iFileAge = -1;
+
+  for (iFile = 0; iFile < files->iNumInputs; iFile++) {
+    if (options[OPT_AGE].iLine[iFile] > 0) {
+      iFileAge = iFile;
+    }
+  }
+  if (iFileAge == -1) {
+    DoubleLineExit(options[OPT_STOPAGE].cFile[iFileStopAge],
+                   options[OPT_AGE].cFile[iFileAge],
+                   options[OPT_STOPAGE].iLine[iFileStopAge],
+                   options[OPT_AGE].iLine[iFileAge]);
+  }
+
+  control->Evolve.dStopTime = fabs(control->Evolve.dStopAge - system->dAge);
+}
+
+void AssignStopTime(CONTROL *control, FILES *files, OPTIONS *options,
+                    SYSTEM *system) {
+  int iFile, iFileStopTime = -1, iFileStopAge = -1;
+
+  for (iFile = 0; iFile < files->iNumInputs; iFile++) {
+    if (options[OPT_STOPTIME].iLine[iFile] > 0) {
+      iFileStopTime = iFile;
+    }
+    if (options[OPT_STOPAGE].iLine[iFile] > 0) {
+      iFileStopAge = iFile;
+    }
+  }
+
+  if (iFileStopTime >= 0 && iFileStopAge == -1) {
+    return; // control->Evolve.dStopTime assigned in ReadStopTime
+  }
+
+  if (iFileStopTime >= 0 && iFileStopAge >= 0) {
+    DoubleLineExit(options[OPT_STOPTIME].cFile[iFileStopTime],
+                   options[OPT_STOPAGE].cFile[iFileStopAge],
+                   options[OPT_STOPTIME].iLine[iFileStopTime],
+                   options[OPT_STOPAGE].iLine[iFileStopAge]);
+  }
+
+  if (iFileStopTime == -1 && iFileStopAge == -1) {
+    control->Evolve.dStopTime = options->dDefault;
+  }
+
+  if (iFileStopTime == -1 && iFileStopAge >= 0) {
+    StopTimeFromAge(control, files, options, system, iFileStopAge);
+  }
+}
+
+void VerifyIntegration(BODY *body, CONTROL *control, FILES *files,
+                       OPTIONS *options, SYSTEM *system,
+                       fnIntegrate *fnOneStep) {
+
+  VerifyIntegrationDirection(control, files, options);
+  if (control->Evolve.iDir == NO_INTEGRATION) {
+    return;
+  }
+
+  PrintIntegrationWarnings(control, files, options);
+  AssignStopTime(control, files, options, system);
+  VerifyOutputTime(control, files, options);
+  AssignIntegrationMethod(control, options, fnOneStep);
+  InitializeIntegrationFiles(body, control, files, options, system);
+}
+
 /*
  *
  * Mass-Radius Relationships
diff --git a/src/vplanet.h b/src/vplanet.h
index ff52fb16..11c4a80b 100644
--- a/src/vplanet.h
+++ b/src/vplanet.h
@@ -1027,6 +1027,7 @@ typedef double (*fnLaplaceFunction)(double, int);
 
 struct SYSTEM {
   char *cName; /**< System's Name */
+  double dAge;
 
   int iNumBodies; /** Number of bodies in the system; redundant with Evolve! */
 
@@ -1741,6 +1742,7 @@ struct EVOLVE {
   double dTime;     /**< Integration Time */
   double dEta;      /**< Variable Timestep Coefficient */
   double dStopTime; /**< Integration Stop Time */
+  double dStopAge;  /**< Age to stop at */
   double dTimeStep; /**< Integration Time step */
   int bVarDt;       /**< Use Variable Timestep? */
   int iTotalSteps;  /**< Total Number of Steps */
diff --git a/tests/AtmescStellar/Backward/test_Backward.py b/tests/AtmescStellar/Backward/test_Backward.py
index 28827145..7afd986b 100644
--- a/tests/AtmescStellar/Backward/test_Backward.py
+++ b/tests/AtmescStellar/Backward/test_Backward.py
@@ -1,180 +1,312 @@
-from benchmark import Benchmark, benchmark 
-import astropy.units as u 
- 
-@benchmark( 
-   { 
-       "log.initial.system.Age": {"value": 3.155760e+16, "unit": u.sec}, 
-       "log.initial.system.Time": {"value": -0.000000, "unit": u.sec}, 
-       "log.initial.system.TotAngMom": {"value": 5.058770e+42, "unit": (u.kg * u.m ** 2) / u.sec}, 
-       "log.initial.system.TotEnergy": {"value": -2.463825e+41, "unit": u.erg}, 
-       "log.initial.system.PotEnergy": {"value": -2.465658e+41, "unit": u.Joule}, 
-       "log.initial.system.KinEnergy": {"value": 1.833857e+38, "unit": u.Joule}, 
-       "log.initial.system.DeltaTime": {"value": 0.000000, "unit": u.sec}, 
-       "log.initial.a.Mass": {"value": 1.988416e+30, "unit": u.kg}, 
-       "log.initial.a.Radius": {"value": 100.680639, "unit": u.Rearth}, 
-       "log.initial.a.RadGyra": {"value": 0.290831}, 
-       "log.initial.a.RotAngMom": {"value": 5.043469e+42, "unit": (u.kg * u.m ** 2) / u.sec}, 
-       "log.initial.a.RotVel": {"value": 4.669855e+04, "unit": u.m / u.sec}, 
-       "log.initial.a.BodyType": {"value": 0.000000}, 
-       "log.initial.a.RotRate": {"value": 7.272205e-05, "unit": 1 / u.sec}, 
-       "log.initial.a.RotPer": {"value": 1.000000, "unit": u.day}, 
-       "log.initial.a.Density": {"value": 1792.696544, "unit": u.kg / u.m ** 3}, 
-       "log.initial.a.HZLimitDryRunaway": {"value": 1.175525e+11, "unit": u.m}, 
-       "log.initial.a.HZLimRecVenus": {"value": 9.760054e+10, "unit": u.m}, 
-       "log.initial.a.HZLimRunaway": {"value": 1.279108e+11, "unit": u.m}, 
-       "log.initial.a.HZLimMoistGreenhouse": {"value": 1.291407e+11, "unit": u.m}, 
-       "log.initial.a.HZLimMaxGreenhouse": {"value": 2.207286e+11, "unit": u.m}, 
-       "log.initial.a.HZLimEarlyMars": {"value": 2.408255e+11, "unit": u.m}, 
-       "log.initial.a.Instellation": {"value": -1.000000, "unit": u.kg / u.sec ** 3}, 
-       "log.initial.a.CriticalSemiMajorAxis": {"value": -1.000000, "unit": u.m}, 
-       "log.initial.a.LXUVTot": {"value": 1.697389e+22, "unit": u.kg / u.sec ** 3}, 
-       "log.initial.a.LostEnergy": {"value": 5.562685e-309, "unit": u.Joule}, 
-       "log.initial.a.LostAngMom": {"value": 5.562685e-309, "unit": (u.kg * u.m ** 2) / u.sec}, 
-       "log.initial.a.EscapeVelocity": {"value": 6.429130e+05, "unit": u.m / u.sec}, 
-       "log.initial.a.Luminosity": {"value": 0.749501, "unit": u.LSUN}, 
-       "log.initial.a.LXUVStellar": {"value": 1.697389e+22, "unit": u.W}, 
-       "log.initial.a.Temperature": {"value": 5591.770652, "unit": u.K}, 
-       "log.initial.a.LXUVFrac": {"value": 5.888437e-05}, 
-       "log.initial.a.RossbyNumber": {"value": 0.063834}, 
-       "log.initial.a.DRotPerDtStellar": {"value": 4.087843e-11}, 
-       "log.initial.a.WindTorque": {"value": 2.390187e+27}, 
-       "log.initial.b.Mass": {"value": 1.818180, "unit": u.Mearth}, 
-       "log.initial.b.Radius": {"value": 6.987673e+07, "unit": u.m}, 
-       "log.initial.b.RadGyra": {"value": 0.500000}, 
-       "log.initial.b.BodyType": {"value": 0.000000}, 
-       "log.initial.b.Density": {"value": 7.597731, "unit": u.kg / u.m ** 3}, 
-       "log.initial.b.HZLimitDryRunaway": {"value": 1.175525e+11, "unit": u.m}, 
-       "log.initial.b.HZLimRecVenus": {"value": 9.760054e+10, "unit": u.m}, 
-       "log.initial.b.HZLimRunaway": {"value": 1.279108e+11, "unit": u.m}, 
-       "log.initial.b.HZLimMoistGreenhouse": {"value": 1.291407e+11, "unit": u.m}, 
-       "log.initial.b.HZLimMaxGreenhouse": {"value": 2.207286e+11, "unit": u.m}, 
-       "log.initial.b.HZLimEarlyMars": {"value": 2.408255e+11, "unit": u.m}, 
-       "log.initial.b.Instellation": {"value": 1.024993e+05, "unit": u.kg / u.sec ** 3}, 
-       "log.initial.b.MeanMotion": {"value": 6.296061e-06, "unit": 1 / u.sec}, 
-       "log.initial.b.OrbPeriod": {"value": 9.979550e+05, "unit": u.sec}, 
-       "log.initial.b.SemiMajorAxis": {"value": 1.495979e+10, "unit": u.m}, 
-       "log.initial.b.LXUVTot": {"value": -1.000000, "unit": u.kg / u.sec ** 3}, 
-       "log.initial.b.EscapeVelocity": {"value": 4554.450923, "unit": u.m / u.sec}, 
-       "log.initial.b.SurfWaterMass": {"value": 0.000000, "unit": u.kg}, 
-       "log.initial.b.EnvelopeMass": {"value": 0.818180, "unit": u.Mearth}, 
-       "log.initial.b.OxygenMass": {"value": 0.000000, "unit": u.kg}, 
-       "log.initial.b.RGLimit": {"value": 1.228298e+11, "unit": u.m}, 
-       "log.initial.b.XO": {"value": 0.000000}, 
-       "log.initial.b.EtaO": {"value": 0.000000}, 
-       "log.initial.b.PlanetRadius": {"value": 6.987673e+07, "unit": u.m}, 
-       "log.initial.b.OxygenMantleMass": {"value": 0.000000, "unit": u.kg}, 
-       "log.initial.b.RadXUV": {"value": -1.000000, "unit": u.m}, 
-       "log.initial.b.RadSolid": {"value": -1.000000, "unit": u.m}, 
-       "log.initial.b.PresXUV": {"value": 5.000000, "unit": (null)}, 
-       "log.initial.b.ScaleHeight": {"value": -1.000000, "unit": u.m}, 
-       "log.initial.b.ThermTemp": {"value": 400.000000, "unit": u.K}, 
-       "log.initial.b.AtmGasConst": {"value": 4124.000000, "unit": (null)}, 
-       "log.initial.b.PresSurf": {"value": -1.000000, "unit": u.Pa}, 
-       "log.initial.b.DEnvMassDt": {"value": -1.965548e+09, "unit": u.kg / u.sec}, 
-       "log.initial.b.FXUV": {"value": 6.035605, "unit": u.W / u.m ** 2}, 
-       "log.initial.b.AtmXAbsEffH2O": {"value": 0.300000}, 
-       "log.initial.b.RocheRadius": {"value": 1.826583e+08, "unit": u.m}, 
-       "log.initial.b.BondiRadius": {"value": 7.395371e+07, "unit": u.m}, 
-       "log.initial.b.HEscapeRegime": {"value": 3.000000}, 
-       "log.initial.b.RRCriticalFlux": {"value": 0.062660, "unit": u.W / u.m ** 2}, 
-       "log.initial.b.CrossoverMass": {"value": 0.000000, "unit": u.kg}, 
-       "log.initial.b.WaterEscapeRegime": {"value": 8.000000}, 
-       "log.initial.b.FXUVCRITDRAG": {"value": 0.000670, "unit": u.W / u.m ** 2}, 
-       "log.initial.b.HREFFLUX": {"value": 5.787364e+19, "unit": 1 / u.m ** 2 / u.sec}, 
-       "log.initial.b.XO2": {"value": 0.000000}, 
-       "log.initial.b.XH2O": {"value": 0.000000}, 
-       "log.initial.b.HDiffFlux": {"value": 2.874041e+15, "unit": 1 / u.m ** 2 / u.sec}, 
-       "log.initial.b.HRefODragMod": {"value": 1.000000}, 
-       "log.initial.b.KTide": {"value": 0.454161}, 
-       "log.initial.b.CumulativeXUVFlux": {"value": 0.000000, "unit": u.kg / u.sec ** 3}, 
-       "log.initial.b.RunawayGreenhouseFlux": {"value": 1520.422432, "unit": u.kg / u.sec ** 3}, 
-       "log.initial.b.RGDuration": {"value": 0.00000e+00, "unit": u.yr}, 
-       "log.final.system.Age": {"value": 1.577880e+14, "unit": u.sec, "rtol": 1e-4}, 
-       "log.final.system.Time": {"value": -3.139981e+16, "unit": u.sec, "rtol": 1e-4}, 
-       "log.final.system.TotAngMom": {"value": 5.072047e+42, "unit": (u.kg * u.m ** 2) / u.sec, "rtol": 1e-4}, 
-       "log.final.system.TotEnergy": {"value": -2.464150e+41, "unit": u.erg, "rtol": 1e-4}, 
-       "log.final.system.PotEnergy": {"value": -1.587322e+41, "unit": u.Joule, "rtol": 1e-4}, 
-       "log.final.system.KinEnergy": {"value": 5.860165e+40, "unit": u.Joule, "rtol": 1e-4}, 
-       "log.final.a.Mass": {"value": 1.988416e+30, "unit": u.kg, "rtol": 1e-4}, 
-       "log.final.a.Radius": {"value": 156.391718, "unit": u.Rearth, "rtol": 1e-4}, 
-       "log.final.a.RadGyra": {"value": 0.449294, "rtol": 1e-4}, 
-       "log.final.a.RotAngMom": {"value": 2.163512e+44, "unit": (u.kg * u.m ** 2) / u.sec, "rtol": 1e-4}, 
-       "log.final.a.RotVel": {"value": 5.403630e+05, "unit": u.m / u.sec, "rtol": 1e-4}, 
-       "log.final.a.BodyType": {"value": 0.000000, "rtol": 1e-4}, 
-       "log.final.a.RotRate": {"value": 0.000542, "unit": 1 / u.sec, "rtol": 1e-4}, 
-       "log.final.a.RotPer": {"value": 0.134241, "unit": u.day, "rtol": 1e-4}, 
-       "log.final.a.Density": {"value": 478.303342, "unit": u.kg / u.m ** 3, "rtol": 1e-4}, 
-       "log.final.a.HZLimitDryRunaway": {"value": 1.084460e+11, "unit": u.m, "rtol": 1e-4}, 
-       "log.final.a.HZLimRecVenus": {"value": 9.446714e+10, "unit": u.m, "rtol": 1e-4}, 
-       "log.final.a.HZLimRunaway": {"value": 1.254168e+11, "unit": u.m, "rtol": 1e-4}, 
-       "log.final.a.HZLimMoistGreenhouse": {"value": 1.249951e+11, "unit": u.m, "rtol": 1e-4}, 
-       "log.final.a.HZLimMaxGreenhouse": {"value": 2.269129e+11, "unit": u.m, "rtol": 1e-4}, 
-       "log.final.a.HZLimEarlyMars": {"value": 2.475151e+11, "unit": u.m, "rtol": 1e-4}, 
-       "log.final.a.Instellation": {"value": -1.000000, "unit": u.kg / u.sec ** 3, "rtol": 1e-4}, 
-       "log.final.a.CriticalSemiMajorAxis": {"value": -1.000000, "unit": u.m, "rtol": 1e-4}, 
-       "log.final.a.LXUVTot": {"value": 2.453266e+23, "unit": u.kg / u.sec ** 3, "rtol": 1e-4}, 
-       "log.final.a.LostEnergy": {"value": -1.462843e+41, "unit": u.Joule, "rtol": 1e-4}, 
-       "log.final.a.LostAngMom": {"value": -2.113232e+44, "unit": (u.kg * u.m ** 2) / u.sec, "rtol": 1e-4}, 
-       "log.final.a.EscapeVelocity": {"value": 5.158439e+05, "unit": u.m / u.sec, "rtol": 1e-4}, 
-       "log.final.a.Luminosity": {"value": 0.637875, "unit": u.LSUN, "rtol": 1e-4}, 
-       "log.final.a.LXUVStellar": {"value": 2.453266e+23, "unit": u.W, "rtol": 1e-4}, 
-       "log.final.a.Temperature": {"value": 4309.979229, "unit": u.K, "rtol": 1e-4}, 
-       "log.final.a.LXUVFrac": {"value": 0.001000, "rtol": 1e-4}, 
-       "log.final.a.RossbyNumber": {"value": 0.003889, "rtol": 1e-4}, 
-       "log.final.a.DRotPerDtStellar": {"value": -4.682593e-11, "rtol": 1e-4}, 
-       "log.final.a.WindTorque": {"value": 1.342986e+28, "rtol": 1e-4}, 
-       "log.final.b.Mass": {"value": 5.236487, "unit": u.Mearth, "rtol": 1e-4}, 
-       "log.final.b.Radius": {"value": 1.902090e+08, "unit": u.m, "rtol": 1e-4}, 
-       "log.final.b.RadGyra": {"value": 0.500000, "rtol": 1e-4}, 
-       "log.final.b.BodyType": {"value": 0.000000, "rtol": 1e-4}, 
-       "log.final.b.Density": {"value": 1.084905, "unit": u.kg / u.m ** 3, "rtol": 1e-4}, 
-       "log.final.b.HZLimitDryRunaway": {"value": 1.084460e+11, "unit": u.m, "rtol": 1e-4}, 
-       "log.final.b.HZLimRecVenus": {"value": 9.446714e+10, "unit": u.m, "rtol": 1e-4}, 
-       "log.final.b.HZLimRunaway": {"value": 1.254168e+11, "unit": u.m, "rtol": 1e-4}, 
-       "log.final.b.HZLimMoistGreenhouse": {"value": 1.249951e+11, "unit": u.m, "rtol": 1e-4}, 
-       "log.final.b.HZLimMaxGreenhouse": {"value": 2.269129e+11, "unit": u.m, "rtol": 1e-4}, 
-       "log.final.b.HZLimEarlyMars": {"value": 2.475151e+11, "unit": u.m, "rtol": 1e-4}, 
-       "log.final.b.Instellation": {"value": 8.723364e+04, "unit": u.kg / u.sec ** 3, "rtol": 1e-4}, 
-       "log.final.b.MeanMotion": {"value": 6.296093e-06, "unit": 1 / u.sec, "rtol": 1e-4}, 
-       "log.final.b.OrbPeriod": {"value": 9.979499e+05, "unit": u.sec, "rtol": 1e-4}, 
-       "log.final.b.SemiMajorAxis": {"value": 1.495979e+10, "unit": u.m, "rtol": 1e-4}, 
-       "log.final.b.LXUVTot": {"value": -1.000000, "unit": u.kg / u.sec ** 3, "rtol": 1e-4}, 
-       "log.final.b.EscapeVelocity": {"value": 4684.770900, "unit": u.m / u.sec, "rtol": 1e-4}, 
-       "log.final.b.SurfWaterMass": {"value": 0.000000, "unit": u.kg, "rtol": 1e-4}, 
-       "log.final.b.EnvelopeMass": {"value": 4.236487, "unit": u.Mearth, "rtol": 1e-4}, 
-       "log.final.b.OxygenMass": {"value": 0.000000, "unit": u.kg, "rtol": 1e-4}, 
-       "log.final.b.RGLimit": {"value": 1.174025e+11, "unit": u.m, "rtol": 1e-4}, 
-       "log.final.b.XO": {"value": 0.000000, "rtol": 1e-4}, 
-       "log.final.b.EtaO": {"value": 0.000000, "rtol": 1e-4}, 
-       "log.final.b.PlanetRadius": {"value": 1.902090e+08, "unit": u.m, "rtol": 1e-4}, 
-       "log.final.b.OxygenMantleMass": {"value": 0.000000, "unit": u.kg, "rtol": 1e-4}, 
-       "log.final.b.RadXUV": {"value": -1.000000, "unit": u.m, "rtol": 1e-4}, 
-       "log.final.b.RadSolid": {"value": -1.000000, "unit": u.m, "rtol": 1e-4}, 
-       "log.final.b.PresXUV": {"value": 5.000000, "unit": (null), "rtol": 1e-4}, 
-       "log.final.b.ScaleHeight": {"value": -1.000000, "unit": u.m, "rtol": 1e-4}, 
-       "log.final.b.ThermTemp": {"value": 400.000000, "unit": u.K, "rtol": 1e-4}, 
-       "log.final.b.AtmGasConst": {"value": 4124.000000, "unit": (null), "rtol": 1e-4}, 
-       "log.final.b.PresSurf": {"value": -1.000000, "unit": u.Pa, "rtol": 1e-4}, 
-       "log.final.b.DEnvMassDt": {"value": 3.419958e+09, "unit": u.kg / u.sec, "rtol": 1e-4}, 
-       "log.final.b.FXUV": {"value": 87.233642, "unit": u.W / u.m ** 2, "rtol": 1e-4}, 
-       "log.final.b.AtmXAbsEffH2O": {"value": 0.300000, "rtol": 1e-4}, 
-       "log.final.b.RocheRadius": {"value": 2.598805e+08, "unit": u.m, "rtol": 1e-4}, 
-       "log.final.b.BondiRadius": {"value": 2.217192e+08, "unit": u.m, "rtol": 1e-4}, 
-       "log.final.b.HEscapeRegime": {"value": 6.000000, "rtol": 1e-4}, 
-       "log.final.b.RRCriticalFlux": {"value": 0.001249, "unit": u.W / u.m ** 2, "rtol": 1e-4}, 
-       "log.final.b.CrossoverMass": {"value": 0.000000, "unit": u.kg, "rtol": 1e-4}, 
-       "log.final.b.WaterEscapeRegime": {"value": 8.000000, "rtol": 1e-4}, 
-       "log.final.b.FXUVCRITDRAG": {"value": 0.000275, "unit": u.W / u.m ** 2, "rtol": 1e-4}, 
-       "log.final.b.HREFFLUX": {"value": 3.590455e+21, "unit": 1 / u.m ** 2 / u.sec, "rtol": 1e-4}, 
-       "log.final.b.XO2": {"value": 0.000000, "rtol": 1e-4}, 
-       "log.final.b.XH2O": {"value": 0.000000, "rtol": 1e-4}, 
-       "log.final.b.HDiffFlux": {"value": 1.117118e+15, "unit": 1 / u.m ** 2 / u.sec, "rtol": 1e-4}, 
-       "log.final.b.HRefODragMod": {"value": 1.000000, "rtol": 1e-4}, 
-       "log.final.b.KTide": {"value": 0.100000, "rtol": 1e-4}, 
-       "log.final.b.CumulativeXUVFlux": {"value": 8.269339e+17, "unit": u.kg / u.sec ** 3, "rtol": 1e-4}, 
-       "log.final.b.RunawayGreenhouseFlux": {"value": 1416.380837, "unit": u.kg / u.sec ** 3, "rtol": 1e-4}, 
-       "log.final.b.RGDuration": {"value": 0.00000e+00, "unit": u.yr, "rtol": 1e-4}, 
-   } 
+import astropy.units as u
+from benchmark import Benchmark, benchmark
+
+
+@benchmark(
+    {
+        "log.initial.system.Age": {"value": 3.155760e16, "unit": u.sec},
+        "log.initial.system.Time": {"value": -0.000000, "unit": u.sec},
+        "log.initial.system.TotAngMom": {
+            "value": 5.058770e42,
+            "unit": (u.kg * u.m**2) / u.sec,
+        },
+        "log.initial.system.TotEnergy": {"value": -2.463825e41, "unit": u.erg},
+        "log.initial.system.PotEnergy": {"value": -2.465658e41, "unit": u.Joule},
+        "log.initial.system.KinEnergy": {"value": 1.833857e38, "unit": u.Joule},
+        "log.initial.system.DeltaTime": {"value": 0.000000, "unit": u.sec},
+        "log.initial.a.Mass": {"value": 1.988416e30, "unit": u.kg},
+        "log.initial.a.Radius": {"value": 100.680639, "unit": u.Rearth},
+        "log.initial.a.RadGyra": {"value": 0.290831},
+        "log.initial.a.RotAngMom": {
+            "value": 5.043469e42,
+            "unit": (u.kg * u.m**2) / u.sec,
+        },
+        "log.initial.a.RotVel": {"value": 4.669855e04, "unit": u.m / u.sec},
+        "log.initial.a.BodyType": {"value": 0.000000},
+        "log.initial.a.RotRate": {"value": 7.272205e-05, "unit": 1 / u.sec},
+        "log.initial.a.RotPer": {"value": 1.000000, "unit": u.day},
+        "log.initial.a.Density": {"value": 1792.696544, "unit": u.kg / u.m**3},
+        "log.initial.a.HZLimitDryRunaway": {"value": 1.175525e11, "unit": u.m},
+        "log.initial.a.HZLimRecVenus": {"value": 9.760054e10, "unit": u.m},
+        "log.initial.a.HZLimRunaway": {"value": 1.279108e11, "unit": u.m},
+        "log.initial.a.HZLimMoistGreenhouse": {"value": 1.291407e11, "unit": u.m},
+        "log.initial.a.HZLimMaxGreenhouse": {"value": 2.207286e11, "unit": u.m},
+        "log.initial.a.HZLimEarlyMars": {"value": 2.408255e11, "unit": u.m},
+        "log.initial.a.Instellation": {"value": -1.000000, "unit": u.kg / u.sec**3},
+        "log.initial.a.CriticalSemiMajorAxis": {"value": -1.000000, "unit": u.m},
+        "log.initial.a.LXUVTot": {"value": 1.697389e22, "unit": u.kg / u.sec**3},
+        "log.initial.a.LostEnergy": {"value": 5.562685e-309, "unit": u.Joule},
+        "log.initial.a.LostAngMom": {
+            "value": 5.562685e-309,
+            "unit": (u.kg * u.m**2) / u.sec,
+        },
+        "log.initial.a.EscapeVelocity": {"value": 6.429130e05, "unit": u.m / u.sec},
+        "log.initial.a.Luminosity": {"value": 0.749501, "unit": u.LSUN},
+        "log.initial.a.LXUVStellar": {"value": 1.697389e22, "unit": u.W},
+        "log.initial.a.Temperature": {"value": 5591.770652, "unit": u.K},
+        "log.initial.a.LXUVFrac": {"value": 5.888437e-05},
+        "log.initial.a.RossbyNumber": {"value": 0.063834},
+        "log.initial.a.DRotPerDtStellar": {"value": 4.087843e-11},
+        "log.initial.a.WindTorque": {"value": 2.390187e27},
+        "log.initial.b.Mass": {"value": 1.818180, "unit": u.Mearth},
+        "log.initial.b.Radius": {"value": 6.987673e07, "unit": u.m},
+        "log.initial.b.RadGyra": {"value": 0.500000},
+        "log.initial.b.BodyType": {"value": 0.000000},
+        "log.initial.b.Density": {"value": 7.597731, "unit": u.kg / u.m**3},
+        "log.initial.b.HZLimitDryRunaway": {"value": 1.175525e11, "unit": u.m},
+        "log.initial.b.HZLimRecVenus": {"value": 9.760054e10, "unit": u.m},
+        "log.initial.b.HZLimRunaway": {"value": 1.279108e11, "unit": u.m},
+        "log.initial.b.HZLimMoistGreenhouse": {"value": 1.291407e11, "unit": u.m},
+        "log.initial.b.HZLimMaxGreenhouse": {"value": 2.207286e11, "unit": u.m},
+        "log.initial.b.HZLimEarlyMars": {"value": 2.408255e11, "unit": u.m},
+        "log.initial.b.Instellation": {"value": 1.024993e05, "unit": u.kg / u.sec**3},
+        "log.initial.b.MeanMotion": {"value": 6.296061e-06, "unit": 1 / u.sec},
+        "log.initial.b.OrbPeriod": {"value": 9.979550e05, "unit": u.sec},
+        "log.initial.b.SemiMajorAxis": {"value": 1.495979e10, "unit": u.m},
+        "log.initial.b.LXUVTot": {"value": -1.000000, "unit": u.kg / u.sec**3},
+        "log.initial.b.EscapeVelocity": {"value": 4554.450923, "unit": u.m / u.sec},
+        "log.initial.b.SurfWaterMass": {"value": 0.000000, "unit": u.kg},
+        "log.initial.b.EnvelopeMass": {"value": 0.818180, "unit": u.Mearth},
+        "log.initial.b.OxygenMass": {"value": 0.000000, "unit": u.kg},
+        "log.initial.b.RGLimit": {"value": 1.228298e11, "unit": u.m},
+        "log.initial.b.XO": {"value": 0.000000},
+        "log.initial.b.EtaO": {"value": 0.000000},
+        "log.initial.b.PlanetRadius": {"value": 6.987673e07, "unit": u.m},
+        "log.initial.b.OxygenMantleMass": {"value": 0.000000, "unit": u.kg},
+        "log.initial.b.RadXUV": {"value": -1.000000, "unit": u.m},
+        "log.initial.b.RadSolid": {"value": -1.000000, "unit": u.m},
+        "log.initial.b.ScaleHeight": {"value": -1.000000, "unit": u.m},
+        "log.initial.b.ThermTemp": {"value": 400.000000, "unit": u.K},
+        "log.initial.b.PresSurf": {"value": -1.000000, "unit": u.Pa},
+        "log.initial.b.DEnvMassDt": {"value": -1.965548e09, "unit": u.kg / u.sec},
+        "log.initial.b.FXUV": {"value": 6.035605, "unit": u.W / u.m**2},
+        "log.initial.b.AtmXAbsEffH2O": {"value": 0.300000},
+        "log.initial.b.RocheRadius": {"value": 1.826583e08, "unit": u.m},
+        "log.initial.b.BondiRadius": {"value": 7.395371e07, "unit": u.m},
+        "log.initial.b.HEscapeRegime": {"value": 3.000000},
+        "log.initial.b.RRCriticalFlux": {"value": 0.062660, "unit": u.W / u.m**2},
+        "log.initial.b.CrossoverMass": {"value": 0.000000, "unit": u.kg},
+        "log.initial.b.WaterEscapeRegime": {"value": 8.000000},
+        "log.initial.b.FXUVCRITDRAG": {"value": 0.000670, "unit": u.W / u.m**2},
+        "log.initial.b.HREFFLUX": {"value": 5.787364e19, "unit": 1 / u.m**2 / u.sec},
+        "log.initial.b.XO2": {"value": 0.000000},
+        "log.initial.b.XH2O": {"value": 0.000000},
+        "log.initial.b.HDiffFlux": {"value": 2.874041e15, "unit": 1 / u.m**2 / u.sec},
+        "log.initial.b.HRefODragMod": {"value": 1.000000},
+        "log.initial.b.KTide": {"value": 0.454161},
+        "log.initial.b.CumulativeXUVFlux": {
+            "value": 0.000000,
+            "unit": u.kg / u.sec**3,
+        },
+        "log.initial.b.RunawayGreenhouseFlux": {
+            "value": 1520.422432,
+            "unit": u.kg / u.sec**3,
+        },
+        "log.initial.b.RGDuration": {"value": 0.00000e00, "unit": u.yr},
+        "log.final.system.Age": {"value": 1.577880e14, "unit": u.sec, "rtol": 1e-4},
+        "log.final.system.Time": {"value": -3.139981e16, "unit": u.sec, "rtol": 1e-4},
+        "log.final.system.TotAngMom": {
+            "value": 5.072047e42,
+            "unit": (u.kg * u.m**2) / u.sec,
+            "rtol": 1e-4,
+        },
+        "log.final.system.TotEnergy": {
+            "value": -2.464150e41,
+            "unit": u.erg,
+            "rtol": 1e-4,
+        },
+        "log.final.system.PotEnergy": {
+            "value": -1.587322e41,
+            "unit": u.Joule,
+            "rtol": 1e-4,
+        },
+        "log.final.system.KinEnergy": {
+            "value": 5.860165e40,
+            "unit": u.Joule,
+            "rtol": 1e-4,
+        },
+        "log.final.a.Mass": {"value": 1.988416e30, "unit": u.kg, "rtol": 1e-4},
+        "log.final.a.Radius": {"value": 156.391718, "unit": u.Rearth, "rtol": 1e-4},
+        "log.final.a.RadGyra": {"value": 0.449294, "rtol": 1e-4},
+        "log.final.a.RotAngMom": {
+            "value": 2.163512e44,
+            "unit": (u.kg * u.m**2) / u.sec,
+            "rtol": 1e-4,
+        },
+        "log.final.a.RotVel": {"value": 5.403630e05, "unit": u.m / u.sec, "rtol": 1e-4},
+        "log.final.a.BodyType": {"value": 0.000000, "rtol": 1e-4},
+        "log.final.a.RotRate": {"value": 0.000542, "unit": 1 / u.sec, "rtol": 1e-4},
+        "log.final.a.RotPer": {"value": 0.134241, "unit": u.day, "rtol": 1e-4},
+        "log.final.a.Density": {
+            "value": 478.303342,
+            "unit": u.kg / u.m**3,
+            "rtol": 1e-4,
+        },
+        "log.final.a.HZLimitDryRunaway": {
+            "value": 1.084460e11,
+            "unit": u.m,
+            "rtol": 1e-4,
+        },
+        "log.final.a.HZLimRecVenus": {"value": 9.446714e10, "unit": u.m, "rtol": 1e-4},
+        "log.final.a.HZLimRunaway": {"value": 1.254168e11, "unit": u.m, "rtol": 1e-4},
+        "log.final.a.HZLimMoistGreenhouse": {
+            "value": 1.249951e11,
+            "unit": u.m,
+            "rtol": 1e-4,
+        },
+        "log.final.a.HZLimMaxGreenhouse": {
+            "value": 2.269129e11,
+            "unit": u.m,
+            "rtol": 1e-4,
+        },
+        "log.final.a.HZLimEarlyMars": {"value": 2.475151e11, "unit": u.m, "rtol": 1e-4},
+        "log.final.a.Instellation": {
+            "value": -1.000000,
+            "unit": u.kg / u.sec**3,
+            "rtol": 1e-4,
+        },
+        "log.final.a.CriticalSemiMajorAxis": {
+            "value": -1.000000,
+            "unit": u.m,
+            "rtol": 1e-4,
+        },
+        "log.final.a.LXUVTot": {
+            "value": 2.453266e23,
+            "unit": u.kg / u.sec**3,
+            "rtol": 1e-4,
+        },
+        "log.final.a.LostEnergy": {
+            "value": -1.462843e41,
+            "unit": u.Joule,
+            "rtol": 1e-4,
+        },
+        "log.final.a.LostAngMom": {
+            "value": -2.113232e44,
+            "unit": (u.kg * u.m**2) / u.sec,
+            "rtol": 1e-4,
+        },
+        "log.final.a.EscapeVelocity": {
+            "value": 5.158439e05,
+            "unit": u.m / u.sec,
+            "rtol": 1e-4,
+        },
+        "log.final.a.Luminosity": {"value": 0.637875, "unit": u.LSUN, "rtol": 1e-4},
+        "log.final.a.LXUVStellar": {"value": 2.453266e23, "unit": u.W, "rtol": 1e-4},
+        "log.final.a.Temperature": {"value": 4309.979229, "unit": u.K, "rtol": 1e-4},
+        "log.final.a.LXUVFrac": {"value": 0.001000, "rtol": 1e-4},
+        "log.final.a.RossbyNumber": {"value": 0.003889, "rtol": 1e-4},
+        "log.final.a.DRotPerDtStellar": {"value": -4.682593e-11, "rtol": 1e-4},
+        "log.final.a.WindTorque": {"value": 1.342986e28, "rtol": 1e-4},
+        "log.final.b.Mass": {"value": 5.236487, "unit": u.Mearth, "rtol": 1e-4},
+        "log.final.b.Radius": {"value": 1.902090e08, "unit": u.m, "rtol": 1e-4},
+        "log.final.b.RadGyra": {"value": 0.500000, "rtol": 1e-4},
+        "log.final.b.BodyType": {"value": 0.000000, "rtol": 1e-4},
+        "log.final.b.Density": {
+            "value": 1.084905,
+            "unit": u.kg / u.m**3,
+            "rtol": 1e-4,
+        },
+        "log.final.b.HZLimitDryRunaway": {
+            "value": 1.084460e11,
+            "unit": u.m,
+            "rtol": 1e-4,
+        },
+        "log.final.b.HZLimRecVenus": {"value": 9.446714e10, "unit": u.m, "rtol": 1e-4},
+        "log.final.b.HZLimRunaway": {"value": 1.254168e11, "unit": u.m, "rtol": 1e-4},
+        "log.final.b.HZLimMoistGreenhouse": {
+            "value": 1.249951e11,
+            "unit": u.m,
+            "rtol": 1e-4,
+        },
+        "log.final.b.HZLimMaxGreenhouse": {
+            "value": 2.269129e11,
+            "unit": u.m,
+            "rtol": 1e-4,
+        },
+        "log.final.b.HZLimEarlyMars": {"value": 2.475151e11, "unit": u.m, "rtol": 1e-4},
+        "log.final.b.Instellation": {
+            "value": 8.723364e04,
+            "unit": u.kg / u.sec**3,
+            "rtol": 1e-4,
+        },
+        "log.final.b.MeanMotion": {
+            "value": 6.296093e-06,
+            "unit": 1 / u.sec,
+            "rtol": 1e-4,
+        },
+        "log.final.b.OrbPeriod": {"value": 9.979499e05, "unit": u.sec, "rtol": 1e-4},
+        "log.final.b.SemiMajorAxis": {"value": 1.495979e10, "unit": u.m, "rtol": 1e-4},
+        "log.final.b.LXUVTot": {
+            "value": -1.000000,
+            "unit": u.kg / u.sec**3,
+            "rtol": 1e-4,
+        },
+        "log.final.b.EscapeVelocity": {
+            "value": 4684.770900,
+            "unit": u.m / u.sec,
+            "rtol": 1e-4,
+        },
+        "log.final.b.SurfWaterMass": {"value": 0.000000, "unit": u.kg, "rtol": 1e-4},
+        "log.final.b.EnvelopeMass": {"value": 4.236487, "unit": u.Mearth, "rtol": 1e-4},
+        "log.final.b.OxygenMass": {"value": 0.000000, "unit": u.kg, "rtol": 1e-4},
+        "log.final.b.RGLimit": {"value": 1.174025e11, "unit": u.m, "rtol": 1e-4},
+        "log.final.b.XO": {"value": 0.000000, "rtol": 1e-4},
+        "log.final.b.EtaO": {"value": 0.000000, "rtol": 1e-4},
+        "log.final.b.PlanetRadius": {"value": 1.902090e08, "unit": u.m, "rtol": 1e-4},
+        "log.final.b.OxygenMantleMass": {"value": 0.000000, "unit": u.kg, "rtol": 1e-4},
+        "log.final.b.RadXUV": {"value": -1.000000, "unit": u.m, "rtol": 1e-4},
+        "log.final.b.RadSolid": {"value": -1.000000, "unit": u.m, "rtol": 1e-4},
+        "log.final.b.ScaleHeight": {"value": -1.000000, "unit": u.m, "rtol": 1e-4},
+        "log.final.b.ThermTemp": {"value": 400.000000, "unit": u.K, "rtol": 1e-4},
+        "log.final.b.PresSurf": {"value": -1.000000, "unit": u.Pa, "rtol": 1e-4},
+        "log.final.b.DEnvMassDt": {
+            "value": 3.419958e09,
+            "unit": u.kg / u.sec,
+            "rtol": 1e-4,
+        },
+        "log.final.b.FXUV": {"value": 87.233642, "unit": u.W / u.m**2, "rtol": 1e-4},
+        "log.final.b.AtmXAbsEffH2O": {"value": 0.300000, "rtol": 1e-4},
+        "log.final.b.RocheRadius": {"value": 2.598805e08, "unit": u.m, "rtol": 1e-4},
+        "log.final.b.BondiRadius": {"value": 2.217192e08, "unit": u.m, "rtol": 1e-4},
+        "log.final.b.HEscapeRegime": {"value": 6.000000, "rtol": 1e-4},
+        "log.final.b.RRCriticalFlux": {
+            "value": 0.001249,
+            "unit": u.W / u.m**2,
+            "rtol": 1e-4,
+        },
+        "log.final.b.CrossoverMass": {"value": 0.000000, "unit": u.kg, "rtol": 1e-4},
+        "log.final.b.WaterEscapeRegime": {"value": 8.000000, "rtol": 1e-4},
+        "log.final.b.FXUVCRITDRAG": {
+            "value": 0.000275,
+            "unit": u.W / u.m**2,
+            "rtol": 1e-4,
+        },
+        "log.final.b.HREFFLUX": {
+            "value": 3.590455e21,
+            "unit": 1 / u.m**2 / u.sec,
+            "rtol": 1e-4,
+        },
+        "log.final.b.XO2": {"value": 0.000000, "rtol": 1e-4},
+        "log.final.b.XH2O": {"value": 0.000000, "rtol": 1e-4},
+        "log.final.b.HDiffFlux": {
+            "value": 1.117118e15,
+            "unit": 1 / u.m**2 / u.sec,
+            "rtol": 1e-4,
+        },
+        "log.final.b.HRefODragMod": {"value": 1.000000, "rtol": 1e-4},
+        "log.final.b.KTide": {"value": 0.100000, "rtol": 1e-4},
+        "log.final.b.CumulativeXUVFlux": {
+            "value": 8.269339e17,
+            "unit": u.kg / u.sec**3,
+            "rtol": 1e-4,
+        },
+        "log.final.b.RunawayGreenhouseFlux": {
+            "value": 1416.380837,
+            "unit": u.kg / u.sec**3,
+            "rtol": 1e-4,
+        },
+        "log.final.b.RGDuration": {"value": 0.00000e00, "unit": u.yr, "rtol": 1e-4},
+    }
 )
-class Test_Backward(Benchmark): 
-   pass 
+class Test_Backward(Benchmark):
+    pass
diff --git a/tests/AtmescStellar/Backward/vpl.in b/tests/AtmescStellar/Backward/vpl.in
index 4796939f..aea8c374 100644
--- a/tests/AtmescStellar/Backward/vpl.in
+++ b/tests/AtmescStellar/Backward/vpl.in
@@ -21,5 +21,6 @@ dMinValue      1e-10                        # Minimum value of eccentricity/obli
 bDoBackward    1                             # Perform a forward evolution?
 bVarDt        1                             # Use variable timestepping?
 dEta          0.01                         # Coefficient for variable timestepping
-dStopTime     9.95e8                           # Stop time for evolution
 dOutputTime   1e6                           # Output timesteps (assuming in body files)
+dSystemAge    1e9
+dStopAge      5e6