Repository: VirtualPlanetaryLaboratory/vplanet (language: python)

Issue:
Removing stupid new feature
Hi, the new feature you implemented is terrible. This PR removes it.


Discussion / hints:

Fix patch:
diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml
index be1e768e1..850aeaa92 100644
--- a/.pre-commit-config.yaml
+++ b/.pre-commit-config.yaml
@@ -4,7 +4,7 @@ exclude: "^\
 
 repos:
   - repo: https://github.com/pre-commit/pre-commit-hooks
-    rev: v4.1.0
+    rev: v4.4.0
     hooks:
       - id: check-added-large-files # prevents giant files from being committed.
       - id: check-case-conflict # checks for files that would conflict in case-insensitive filesystems.
@@ -18,17 +18,17 @@ repos:
       - id: trailing-whitespace # trims trailing whitespace.
 
   - repo: https://github.com/psf/black
-    rev: 21.12b0
+    rev: 23.1.0
     hooks:
       - id: black
 
   - repo: https://github.com/PyCQA/isort
-    rev: 5.10.1
+    rev: 5.12.0
     hooks:
       - id: isort
         args: [--profile=black]
 
   - repo: https://github.com/pre-commit/mirrors-clang-format
-    rev: v13.0.0
+    rev: v15.0.7
     hooks:
       - id: clang-format
diff --git a/docs/parametersweep.rst b/docs/parametersweep.rst
index 0ce90c0ee..7bd08b032 100644
--- a/docs/parametersweep.rst
+++ b/docs/parametersweep.rst
@@ -118,6 +118,14 @@ This command returns output like the following to the terminal:
 
 But with the proper numbers shown.
 
+.. warning::
+
+    If you decide to *rerun* a parameter sweep, you must delete the checkpoint file! 
+    The name of this file is ``.sDestFolder``, where sDestFolder is the option in :code:`VSPACE`
+    the provdes the name for the directory that contains the simulations. If you do not delete
+    this file, the :code:`MultiPlanet` will conclude your sweep has finished and will not restart
+    the simulations.
+
 After :code:`MultiPlanet` completes, you may have a large number of directories with gigabytes
 of data. Storing, analyzing, and plotting these data can be tedious as each output file
 from each directory must be opened and read in sequentially. To streamline this process,
diff --git a/examples/FfdProxCen/FfdProxCen.png b/examples/FfdProxCen/FfdProxCen.png
index 514a16523..7f660fb70 100644
Binary files a/examples/FfdProxCen/FfdProxCen.png and b/examples/FfdProxCen/FfdProxCen.png differ
diff --git a/examples/FfdProxCen/README.rst b/examples/FfdProxCen/README.rst
index 2c9e9aa8e..46f326edf 100644
--- a/examples/FfdProxCen/README.rst
+++ b/examples/FfdProxCen/README.rst
@@ -1,4 +1,4 @@
-Flare frequency distribution of Proxima Centauri with Davenport et. al (2019) and Lacy et. al (1976) models .
+Flare frequency distribution of Proxima Centauri
 ====================================================================================================
 
 Overview
@@ -26,11 +26,11 @@ To run this example
 Expected output
 ---------------
 
-.. figure:: FfdReproduced.png
+.. figure:: FfdProxCen.png
    :width:  7200px
    :align: center
 
-   Flare frequency distribution (FFD) of energies between 1e33 and 1e36 ergs for Proxima Centauri .
-   The blue, purple, brown and red lines represent the FFD using Davenport et. al (2019) model , at 1 Myr,
-   10 Myr, 100 Myr and 1Gyr. The orange line represents the Lacy et. al (1976) model, where the FFD is
-   constant along the stellar age.
+   Flare frequency distribution (FFD) for Proxima Centauri between 10^34 and 10^36 ergs.
+   The red, orange, dark blue, and pale blue lines represent the FFD from Davenport et. al (2019) at 1 Myr,
+   10 Myr, 100 Myr and 1Gyr, respectively. The black line is the Lacy et. al (1976) model, where the FFD is
+   constant for all the stellar ages.
diff --git a/examples/FfdProxCen/makeplot.py b/examples/FfdProxCen/makeplot.py
index d90a5a544..3146c8586 100644
--- a/examples/FfdProxCen/makeplot.py
+++ b/examples/FfdProxCen/makeplot.py
@@ -5,21 +5,13 @@
 Laura N.  R. do Amaral, Universidad Nacional Autónoma de México, 2021
 date: Sep. 2021
 """
-import os
 import pathlib
-import subprocess
 import sys
 from math import exp, expm1, log, log10
 
-import matplotlib as mpl
-import matplotlib.colors as colors
 import matplotlib.pyplot as plt
-import matplotlib.style
 import numpy as np
-import seaborn as sns
 from matplotlib.lines import Line2D
-from matplotlib.patches import Patch
-from scipy.interpolate import interp1d
 
 import vplanet
 
@@ -37,23 +29,14 @@
     "./lacy/vpl.in",
 ]
 
-directory = os.getcwd()
 
-print(directory)
 # Run the simulations
 for i in range(0, 2):
     vplanet.run(path / run[i], units=False)
 
 # Loading the data
-directory = os.getcwd()
-
-print(directory)
-data_daven = np.loadtxt("./davenport/davenport.star.forward")
-data_lacy = np.loadtxt("./lacy/lacy.star.forward")
-
-directory = os.getcwd()
-
-print(directory)
+data_daven = np.loadtxt(path / "./davenport/davenport.star.forward")
+data_lacy = np.loadtxt(path / "./lacy/lacy.star.forward")
 
 age_daven = data_daven[:, 0]
 
@@ -79,25 +62,25 @@
 
 
 # Plot
-fig = plt.figure(figsize=(4, 3))
+fig = plt.figure(figsize=(6.5, 6.5))
 
-a = 0.5
+a = 2
 
 fp = [1, 1, 1, 1]
 np.interp(lacy[5:9], lacy[1:5], fp)
-plt.plot(lacy[5:9], lacy[1:5], "-", color="Orange", linewidth=a)
+plt.plot(lacy[5:9], lacy[1:5], "-", color="k", linewidth=a)
 
 np.interp(Mdata[5:9], Mdata[1:5], fp)
-plt.plot(Mdata[5:9], Mdata[1:5], "-", color="cornflowerblue", linewidth=a)
+plt.plot(Mdata[5:9], Mdata[1:5], "-", color=vpl.colors.red, linewidth=a)
 
 np.interp(Ndata[5:9], Ndata[1:5], fp)
-plt.plot(Ndata[5:9], Ndata[1:5], "-", color="darkviolet", linewidth=a)
+plt.plot(Ndata[5:9], Ndata[1:5], "-", color=vpl.colors.orange, linewidth=a)
 
 np.interp(Odata[5:9], Odata[1:5], fp)
-plt.plot(Odata[5:9], Odata[1:5], "-", color="darkred", linewidth=a)
+plt.plot(Odata[5:9], Odata[1:5], "-", color=vpl.colors.dark_blue, linewidth=a)
 
 np.interp(Pdata[5:9], Pdata[1:5], fp)
-plt.plot(Pdata[5:9], Pdata[1:5], "-", color="red", linewidth=a)
+plt.plot(Pdata[5:9], Pdata[1:5], "-", color=vpl.colors.pale_blue, linewidth=a)
 
 
 for i in range(1, 5):
@@ -111,27 +94,55 @@
 # Limits and scale
 plt.yscale("log")
 plt.xscale("log")
-plt.xlim(1e33, 1e36)
+plt.xlim(1e34, 1e36)
 
 
 # Legend
 legend_elements = [
     Line2D(
-        [0], [0], ls="-", color="cornflowerblue", lw=1, label="Davenport mode, 1 Myr"
+        [0],
+        [0],
+        ls="-",
+        color=vpl.colors.red,
+        lw=a,
+        label="Davenport et al. (2019), 1 Myr",
+    ),
+    Line2D(
+        [0],
+        [0],
+        ls="-",
+        color=vpl.colors.orange,
+        lw=a,
+        label="Davenport et al. (2019), 10 Myr",
+    ),
+    Line2D(
+        [0],
+        [0],
+        ls="-",
+        color=vpl.colors.dark_blue,
+        lw=a,
+        label="Davenport et al. (2019), 100 Myr",
+    ),
+    Line2D(
+        [0],
+        [0],
+        ls="-",
+        color=vpl.colors.pale_blue,
+        lw=a,
+        label="Davenport et al. (2019), 1 Gyr",
     ),
-    Line2D([0], [0], ls="-", color="darkviolet", lw=1, label="Davenport mode, 10 Myr"),
-    Line2D([0], [0], ls="-", color="darkred", lw=1, label="Davenport mode, 100 Myr"),
-    Line2D([0], [0], ls="-", color="red", lw=1, label="Davenport mode, 1 Gyr"),
-    Line2D([0], [0], color="orange", lw=1, label="Lacy mode"),
+    Line2D([0], [0], color="k", lw=a, label="Lacy et al. (1976)"),
 ]
 
-plt.legend(handles=legend_elements, ncol=1, loc="lower left", fontsize=7)
+plt.legend(handles=legend_elements, ncol=1, loc="lower left", fontsize=12)
 
 
 # Label
-plt.title("Flare Frequency Distribution")
-plt.xlabel("log Flare Energy (erg)")
-plt.ylabel("Cumulative Flare Freq (#/day)")
+plt.title("Flare Frequency Distribution", fontsize=15)
+plt.xlabel("log(Flare Energy/ergs)", fontsize=15)
+plt.ylabel("Cumulative Flare Freq (#/day)", fontsize=15)
+plt.xticks(fontsize=12)
+plt.yticks(fontsize=12)
 
 # Saving figure
 ext = get_args().ext
diff --git a/src/body.c b/src/body.c
index c5ef74f76..a3f0dd35e 100644
--- a/src/body.c
+++ b/src/body.c
@@ -460,6 +460,7 @@ void BodyCopy(BODY *dest, BODY *src, EVOLVE *evolve) {
     dest[iBody].dLostAngMom   = src[iBody].dLostAngMom;
     dest[iBody].dLostEng      = src[iBody].dLostEng;
     dest[iBody].dAlbedoGlobal = src[iBody].dAlbedoGlobal;
+    dest[iBody].bCalcDynEllip = src[iBody].bCalcDynEllip;
 
     dest[iBody].bBinary   = src[iBody].bBinary;
     dest[iBody].bDistOrb  = src[iBody].bDistOrb;
diff --git a/src/eqtide.c b/src/eqtide.c
index 278b9c84d..04daa4335 100644
--- a/src/eqtide.c
+++ b/src/eqtide.c
@@ -1268,12 +1268,13 @@ void VerifyCTL(BODY *body, CONTROL *control, FILES *files, OPTIONS *options,
     if (options[OPT_TIDALTAU].iLine[iBody + 1] >= 0) {
       /* Tidal Tau was also set */
       if (control->Io.iVerbose >= VERBINPUT) {
+        // System units (SI) by here
         fprintf(stderr,
                 "INFO: Time lag model selected, %s and %s set in file %s. "
-                "Using %s = %lf and ignoring %s.\n",
+                "Using %s = %lf seconds and ignoring %s.\n",
                 options[OPT_TIDALTAU].cName, options[OPT_TIDALQ].cName,
                 options[OPT_TIDALTAU].cFile[iBody + 1],
-                options[OPT_TIDALTAU].cName, body[iBody + 1].dTidalTau,
+                options[OPT_TIDALTAU].cName, body[iBody].dTidalTau,
                 options[OPT_TIDALQ].cName);
       }
     } else {
@@ -1489,11 +1490,11 @@ void VerifyCPL(BODY *body, CONTROL *control, FILES *files, OPTIONS *options,
       /* Tidal Q was also set */
       if (control->Io.iVerbose >= VERBINPUT) {
         fprintf(stderr,
-                "INFO: Phase lag model selected, %s and %s set in file %s. "
+                "WARNING: Phase lag model selected, but both %s and %s set in file %s. "
                 "Using %s = %lf and ignoring %s.\n",
                 options[OPT_TIDALTAU].cName, options[OPT_TIDALQ].cName,
                 options[OPT_TIDALTAU].cFile[iBody + 1],
-                options[OPT_TIDALQ].cName, body[iBody + 1].dTidalQ,
+                options[OPT_TIDALQ].cName, body[iBody].dTidalQ,
                 options[OPT_TIDALTAU].cName);
       }
     } else {
@@ -1590,10 +1591,10 @@ void VerifyCPL(BODY *body, CONTROL *control, FILES *files, OPTIONS *options,
 void VerifyPerturbersEqtide(BODY *body, FILES *files, OPTIONS *options,
                             UPDATE *update, int iNumBodies, int iBody) {
   int iPert, iBodyPert, iVar, ok;
-  int *bFound = malloc(iNumBodies);
+  int *bFound = malloc(iNumBodies*sizeof(int));
 
   for (iBody = 0; iBody < iNumBodies; iBody++) {
-
+    fprintf(stderr,"Body: %s\n",body[iBody].cName);
     if (body[iBody].bEqtide) {
       if (body[iBody].iTidePerts > 0) {
         for (iPert = 0; iPert < body[iBody].iTidePerts; iPert++) {
@@ -1610,7 +1611,7 @@ void VerifyPerturbersEqtide(BODY *body, FILES *files, OPTIONS *options,
                 // Was eqtide set for the perturbed body?
                 if (!body[iBodyPert].bEqtide) {
                   fprintf(stderr,
-                          "ERROR: %s tidally perturbs %s, but module EQTIDE "
+                          "ERROR: %s tidally perturbs %s, but module EqTide "
                           "was not selected for body %s.\n",
                           body[iBody].cName, body[iBodyPert].cName,
                           body[iBodyPert].cName);
@@ -2663,15 +2664,8 @@ void WriteEccTimescaleEqtide(BODY *body, CONTROL *control, OUTPUT *output,
 void WriteEqRotPer(BODY *body, CONTROL *control, OUTPUT *output, SYSTEM *system,
                    UNITS *units, UPDATE *update, int iBody, double *dTmp,
                    char cUnit[]) {
-  int iOrbiter;
-
-  if (!bPrimary(body, iBody)) {
-    iOrbiter = iBody;
-  } else {
-    // Only 1 pertuber allowed -- Maybe check in VerifyOutputEqtide?
-    iOrbiter = body[iBody].iaTidePerts[0];
-  }
 
+  int iOrbiter = fiAssignTidalOrbiter(body,iBody);
   *dTmp = fdFreqToPer(fdEqRotRate(
         body, iBody, body[iOrbiter].dMeanMotion, body[iOrbiter].dEccSq,
         control->Evolve.iEqtideModel, control->Evolve.bDiscreteRot));
@@ -2688,14 +2682,8 @@ void WriteEqRotPer(BODY *body, CONTROL *control, OUTPUT *output, SYSTEM *system,
 void WriteEqRotPerCont(BODY *body, CONTROL *control, OUTPUT *output,
                        SYSTEM *system, UNITS *units, UPDATE *update, int iBody,
                        double *dTmp, char cUnit[]) {
-  int iOrbiter;
 
-  if (!bPrimary(body, iBody)) {
-    iOrbiter = iBody;
-  } else {
-    // Only 1 pertuber allowed -- Maybe check in VerifyOutputEqtide?
-    iOrbiter = body[iBody].iaTidePerts[0];
-  }
+  int iOrbiter = fiAssignTidalOrbiter(body,iBody);
 
   // To CPL, or to CTL? That is the question
   if (control->Evolve.iEqtideModel == CPL) {
@@ -2719,15 +2707,8 @@ void WriteEqRotPerCont(BODY *body, CONTROL *control, OUTPUT *output,
 void WriteEqRotPerDiscrete(BODY *body, CONTROL *control, OUTPUT *output,
                            SYSTEM *system, UNITS *units, UPDATE *update,
                            int iBody, double *dTmp, char cUnit[]) {
-  int iOrbiter;
-
-  if (!bPrimary(body, iBody)) {
-    iOrbiter = iBody;
-  } else {
-    // Only 1 pertuber allowed -- Maybe check in VerifyOutputEqtide?
-    iOrbiter = body[iBody].iaTidePerts[0];
-  }
 
+  int iOrbiter = fiAssignTidalOrbiter(body,iBody);
   if (control->Evolve.iEqtideModel == CPL) {
     *dTmp = fdFreqToPer(fdCPLEqRotRateDiscrete(body[iOrbiter].dMeanMotion,
                                                body[iOrbiter].dEccSq));
@@ -2747,15 +2728,8 @@ void WriteEqRotPerDiscrete(BODY *body, CONTROL *control, OUTPUT *output,
 void WriteEqRotRate(BODY *body, CONTROL *control, OUTPUT *output,
                     SYSTEM *system, UNITS *units, UPDATE *update, int iBody,
                     double *dTmp, char cUnit[]) {
-  int iOrbiter;
-
-  if (!bPrimary(body, iBody)) {
-    iOrbiter = iBody;
-  } else {
-    // Only 1 pertuber allowed -- Maybe check in VerifyOutputEqtide?
-    iOrbiter = body[iBody].iaTidePerts[0];
-  }
 
+  int iOrbiter = fiAssignTidalOrbiter(body,iBody);
   *dTmp = fdEqRotRate(body, iBody, body[iOrbiter].dMeanMotion,
                       body[iOrbiter].dEccSq, control->Evolve.iEqtideModel,
                       control->Evolve.bDiscreteRot);
@@ -2772,14 +2746,8 @@ void WriteEqRotRate(BODY *body, CONTROL *control, OUTPUT *output,
 void WriteEqRotRateCont(BODY *body, CONTROL *control, OUTPUT *output,
                         SYSTEM *system, UNITS *units, UPDATE *update, int iBody,
                         double *dTmp, char cUnit[]) {
-  int iOrbiter;
 
-  if (!bPrimary(body, iBody)) {
-    iOrbiter = iBody;
-  } else {
-    // Only 1 pertuber allowed -- Maybe check in VerifyOutputEqtide?
-    iOrbiter = body[iBody].iaTidePerts[0];
-  }
+  int iOrbiter = fiAssignTidalOrbiter(body,iBody);
 
   // To CPL, or to CTL? That is the question
   if (control->Evolve.iEqtideModel == CPL) {
@@ -2802,14 +2770,7 @@ void WriteEqRotRateCont(BODY *body, CONTROL *control, OUTPUT *output,
 void WriteEqRotRateDiscrete(BODY *body, CONTROL *control, OUTPUT *output,
                             SYSTEM *system, UNITS *units, UPDATE *update,
                             int iBody, double *dTmp, char cUnit[]) {
-  int iOrbiter;
-
-  if (!bPrimary(body, iBody)) {
-    iOrbiter = iBody;
-  } else {
-    // Only 1 pertuber allowed -- Maybe check in VerifyOutputEqtide?
-    iOrbiter = body[iBody].iaTidePerts[0];
-  }
+  int iOrbiter = fiAssignTidalOrbiter(body,iBody);
 
   if (control->Evolve.iEqtideModel == CPL) {
     *dTmp = fdCPLEqRotRateDiscrete(body[iOrbiter].dMeanMotion,
@@ -2830,15 +2791,9 @@ void WriteEqRotRateDiscrete(BODY *body, CONTROL *control, OUTPUT *output,
 void WriteEqTidePower(BODY *body, CONTROL *control, OUTPUT *output,
                       SYSTEM *system, UNITS *units, UPDATE *update, int iBody,
                       double *dTmp, char cUnit[]) {
-  int iOrbiter;
-
-  if (!bPrimary(body, iBody)) {
-    iOrbiter = iBody;
-  } else {
-    // Only 1 pertuber allowed -- Maybe check in VerifyOutputEqtide?
-    iOrbiter = body[iBody].iaTidePerts[0];
-  }
 
+  int iOrbiter = fiAssignTidalOrbiter(body,iBody);
+  // Why is this Eq??????? XXX
   if (control->Evolve.iEqtideModel == CPL) {
     *dTmp = fdCPLTidePowerEq(body[iBody].dTidalZ[iOrbiter], body[iBody].dEccSq,
                              body[iBody].dMeanMotion, body[iBody].dObliquity,
@@ -2928,6 +2883,7 @@ void WriteOblTimescaleEqtide(BODY *body, CONTROL *control, OUTPUT *output,
   }
   */
   *dTmp = -1;
+  strcpy(cUnit, "");
 }
 
 void WriteRotTimescaleEqtide(BODY *body, CONTROL *control, OUTPUT *output,
@@ -2975,6 +2931,16 @@ void WritePowerEqtide(BODY *body, CONTROL *control, OUTPUT *output,
   *dTmp = fdTidePower(body, iBody, control->Evolve.iEqtideModel);
   //*dTmp = body[iBody].dTidalPowMan;
 
+
+  int jBody=0;
+  if (iBody==0)
+    jBody++;
+    /*
+  fprintf(stderr,"\niBody: %d\n",iBody);
+  fprintf(stderr,"\njBody: %d\n",jBody);
+  fprintf(stderr,"TidalZ: %lf\n",body[iBody].dTidalZ[jBody]);
+*/
+
   if (output->bDoNeg[iBody]) {
     *dTmp *= output->dNeg;
     strcpy(cUnit, output->cNeg);
@@ -3467,6 +3433,7 @@ void LogBodyEqtide(BODY *body, CONTROL *control, OUTPUT *output, SYSTEM *system,
   fprintf(fp, "----- EQTIDE PARAMETERS (%s)------\n", body[iBody].cName);
   for (iOut = iStart; iOut < OUTENDEQTIDE; iOut++) {
     if (output[iOut].iNum > 0) {
+      //fprintf(stderr,"iOut = %d.\n",iOut);
       WriteLogEntry(body, control, &output[iOut], system, update, fnWrite[iOut],
                     fp, iBody);
     }
@@ -3937,6 +3904,40 @@ void ForceBehaviorEqtide(BODY *body, MODULE *module, EVOLVE *evolve, IO *io,
   }
 }
 
+/**
+ * Identify and return the index for the perturbing body.
+ * 
+ * 
+*/
+
+int fiAssignTidalPerturber(BODY *body,int iBody) {
+  int iPerturber;
+
+  if (!bPrimary(body, iBody)) {
+    iPerturber = iBody;
+  } else {
+    // Without more info, assume closest body dominates
+    // XXX Should either be calculable from these variables or be rewritten
+    iPerturber = 0;
+  }
+  return iPerturber;
+}
+
+/* Non-central body contains the orbital parameters */
+
+int fiAssignTidalOrbiter(BODY *body,int iBody) {
+  int iOrbiter;
+
+  if (!bPrimary(body, iBody)) {
+    iOrbiter = iBody;
+  } else {
+    // Without more info, assume closest body dominates
+    // XXX Should either be calculable from these variables or be rewritten
+    iOrbiter = 1;
+  }
+  return iOrbiter;
+}
+
 /*
  ************************ CPL Functions ******************
  */
@@ -3957,10 +3958,14 @@ double fdCPLTidePower(BODY *body, int iBody) {
     } else {
       iOrbiter = iBody;
     }
+    // iIndex should be removed by optimization
     iIndex = body[iBody].iaTidePerts[iPert];
 
     // Does this work with DF's changes to da/dt with the synchronous case?
     // See Fleming et al., 2018
+    //fprintf(stderr,"\niBody: %d\n",iBody);
+    //fprintf(stderr,"TidalZ[%d]: %lf\n",iIndex,body[iBody].dTidalZ[iIndex]);
+
     dOrbPow += -body[iBody].dTidalZ[iIndex] / 8 *
                (4 * body[iBody].iTidalEpsilon[iIndex][0] +
                 body[iOrbiter].dEccSq *
@@ -4064,6 +4069,7 @@ double fdCPLTidePowerEq(double dTidalZ, double dEccSq, double dMeanMotion,
   double dGammaOrb, dGammaRot, dRotRateEq;
   int *iEpsilon;
 
+  // XXX Yikes! This shouldn't be getting mallocked every output!
   iEpsilon = malloc(10 * sizeof(int));
 
   /* Must reset spin rate in order to get the phase lags correctly */
diff --git a/src/eqtide.h b/src/eqtide.h
index 9d0e2a91e..bef4fb018 100644
--- a/src/eqtide.h
+++ b/src/eqtide.h
@@ -276,6 +276,9 @@ void ForceBehaviorEqtide(BODY *, MODULE *, EVOLVE *, IO *, SYSTEM *, UPDATE *,
 
 void PropsAuxEqtide(BODY *, EVOLVE *, IO *, UPDATE *, int);
 
+int fiAssignTidalPerturber(BODY*,int);
+int fiAssignTidalOrbiter(BODY*,int);
+
 /******************** CPL Functions ********************/
 
 /* Auxiliary parameters */
diff --git a/src/flare.c b/src/flare.c
index 618fafc10..3f126ef17 100644
--- a/src/flare.c
+++ b/src/flare.c
@@ -614,31 +614,31 @@ void InitializeOptionsFlare(OPTIONS *options, fnReadOption fnRead[]) {
   options[OPT_FLAREBANDPASS].bMultiFile = 1;
   fnRead[OPT_FLAREBANDPASS]             = &ReadFlareBandPass;
   sprintf(options[OPT_FLAREBANDPASS].cLongDescr,
-          /*
-                  "If UV or GOES is selected, the code will convert \n"
-                  "the input energy of flares from the UV band \n"
-                  "(3000-4300 Å) or GOES band (1-8 Å) to the Kepler band \n"
-                  "(4000-9000 Å) to calculate the FFD and the SXR \n"
-                  "band (1.24 - 1239.85 Å) to calculate the luminosity.\n"
-                  " If SXR is selected, the code will convert the \n"
-                  "input energy of flares from the SXR band to the Kepler\n"
-                  "band to calculate the FFD and will use the same \n"
-                  "value of energy to calculate the luminosity .\n"
-                  "If KEPLER is selected, the code will use the same\n"
-                  "input energy to calculate the FFD and will convert \n"
-                  "the input energy to SXR band to calculate the \n"
-                  "luminosity. The conversion values are taken from \n"
-                  "Osten and Wolk (2015) (doi:10.1088/0004-637X/809/1/79).\n"
-                  "If TESSUV is selected, then the code will convert\n"
-                  "the input energy of flares using the conversion value\n"
-                  "to the band U (2000-2800 Å) to the TESS data from \n"
-                  "Gunther et al 2020
-             (https://doi.org/10.3847/1538-3881/ab5d3a). \n" "If the BOLOMETRIC
-             its selected, the code will convert the input \n" "energy of flares
-             using the conversion values are taken \n" "from Osten and Wolk
-             (2015) (doi:10.1088/0004-637X/809/1/79).\n"
-                  */
-          "Currently suppressed due to execution errors.");
+
+          "If UV or GOES is selected, the code will convert \n"
+          "the input energy of flares from the UV band \n"
+          "(3000-4300 Å) or GOES band (1-8 Å) to the Kepler band \n"
+          "(4000-9000 Å) to calculate the FFD and the SXR \n"
+          "band (1.24 - 1239.85 Å) to calculate the luminosity.\n"
+          " If SXR is selected, the code will convert the \n"
+          "input energy of flares from the SXR band to the Kepler\n"
+          "band to calculate the FFD and will use the same \n"
+          "value of energy to calculate the luminosity .\n"
+          "If KEPLER is selected, the code will use the same\n"
+          "input energy to calculate the FFD and will convert \n"
+          "the input energy to SXR band to calculate the \n"
+          "luminosity. The conversion values are taken from \n"
+          "Osten and Wolk (2015) (doi:10.1088/0004-637X/809/1/79).\n"
+          "If TESSUV is selected, then the code will convert\n"
+          "the input energy of flares using the conversion value\n"
+          "to the band U (2000-2800 Å) to the TESS data from \n"
+          "Gunther et al 2020 (https://doi.org/10.3847/1538-3881/ab5d3a). \n"
+          "If the BOLOMETRIC its selected, the code will convert the input \n"
+          "energy of flares using the conversion values are taken \n"
+          "from Osten and Wolk (2015) (doi:10.1088/0004-637X/809/1/79).\n"
+
+          //"Currently suppressed due to execution errors."
+  );
   /*
     sprintf(options[OPT_FLARESLOPEUNITS].cName, "sFlareSlopeUnits");
     sprintf(options[OPT_FLARESLOPEUNITS].cDescr,
diff --git a/src/options.c b/src/options.c
index 148742e0a..10901f10e 100644
--- a/src/options.c
+++ b/src/options.c
@@ -436,6 +436,8 @@ int iGetNumLines(char cFile[]) {
   FILE *fp;
   char cLine[LINE];
 
+  fprintf(stderr,"File: %s\n",cFile);
+
   fp = fopen(cFile, "r");
   if (fp == NULL) {
     fprintf(stderr, "Unable to open %s.\n", cFile);
@@ -487,6 +489,8 @@ void InitializeInput(INFILE *input) {
   FILE *fp;
   char cLine[LINE];
 
+fprintf(stderr,"File: %s\n",input->cIn);
+
   fp = fopen(input->cIn, "r");
   if (fp == NULL) {
     fprintf(stderr, "Unable to open %s.\n", input->cIn);
diff --git a/src/output.c b/src/output.c
index 073c831e4..487df799e 100644
--- a/src/output.c
+++ b/src/output.c
@@ -2453,9 +2453,11 @@ void WriteOutput(BODY *body, CONTROL *control, FILES *files, OUTPUT *output,
   }
 }
 
-void InitializeOutput(OUTPUT *output, fnWriteOutput fnWrite[]) {
+void InitializeOutput(FILES *files, OUTPUT *output, fnWriteOutput fnWrite[]) {
   int iOut, iBody, iModule;
 
+  memset(files->cLog, '\0', NAMELEN);
+
   for (iOut = 0; iOut < MODULEOUTEND; iOut++) {
     memset(output[iOut].cName, '\0', OPTLEN);
     sprintf(output[iOut].cName, "null");
diff --git a/src/output.h b/src/output.h
index 487047c89..428b325db 100644
--- a/src/output.h
+++ b/src/output.h
@@ -121,7 +121,7 @@ void WriteOutput(BODY *, CONTROL *, FILES *, OUTPUT *, SYSTEM *, UPDATE *,
                  fnWriteOutput *, double, double);
 void WriteLog(BODY *, CONTROL *, FILES *, MODULE *, OPTIONS *, OUTPUT *,
               SYSTEM *, UPDATE *, fnUpdateVariable ***, fnWriteOutput *, int);
-void InitializeOutput(OUTPUT *, fnWriteOutput *);
+void InitializeOutput(FILES*,OUTPUT *, fnWriteOutput *);
 void WriteTime(BODY *, CONTROL *, OUTPUT *, SYSTEM *, UNITS *, UPDATE *, int,
                double *, char[]);
 
diff --git a/src/poise.c b/src/poise.c
index e02da311b..f0cc4e5da 100644
--- a/src/poise.c
+++ b/src/poise.c
@@ -16,26 +16,26 @@ void BodyCopyPoise(BODY *dest, BODY *src, int iTideModel, int iNumBodies,
   dest[iBody].bReadOrbitOblData = src[iBody].bReadOrbitOblData;
 }
 
-void InitializeUpdateTmpBodyPoise(BODY *body,CONTROL *control,UPDATE *update,\
+void InitializeUpdateTmpBodyPoise(BODY *body, CONTROL *control, UPDATE *update,
                                   int iBody) {
   if (body[iBody].bReadOrbitOblData) {
     int iLine;
 
-    control->Evolve.tmpBody[iBody].daSemiSeries = \
-    malloc(body[iBody].iNLines*sizeof(double));
-    control->Evolve.tmpBody[iBody].daHeccSeries = \
-    malloc(body[iBody].iNLines*sizeof(double));
-    control->Evolve.tmpBody[iBody].daKeccSeries = \
-    malloc(body[iBody].iNLines*sizeof(double));
+    control->Evolve.tmpBody[iBody].daSemiSeries =
+          malloc(body[iBody].iNLines * sizeof(double));
+    control->Evolve.tmpBody[iBody].daHeccSeries =
+          malloc(body[iBody].iNLines * sizeof(double));
+    control->Evolve.tmpBody[iBody].daKeccSeries =
+          malloc(body[iBody].iNLines * sizeof(double));
 
 
-    for (iLine=0;iLine<body[iBody].iNLines;iLine++) {
+    for (iLine = 0; iLine < body[iBody].iNLines; iLine++) {
       control->Evolve.tmpBody[iBody].daSemiSeries[iLine] =
-      body[iBody].daSemiSeries[iLine];
-      control->Evolve.tmpBody[iBody].daHeccSeries[iLine] = \
-      body[iBody].daHeccSeries[iLine];
-      control->Evolve.tmpBody[iBody].daKeccSeries[iLine] = \
-      body[iBody].daKeccSeries[iLine];
+            body[iBody].daSemiSeries[iLine];
+      control->Evolve.tmpBody[iBody].daHeccSeries[iLine] =
+            body[iBody].daHeccSeries[iLine];
+      control->Evolve.tmpBody[iBody].daKeccSeries[iLine] =
+            body[iBody].daKeccSeries[iLine];
     }
   }
 }
@@ -118,37 +118,36 @@ void ReadSurfAlbedo(BODY *body, CONTROL *control, FILES *files,
   }
 }
 
-void ReadOrbitOblData(BODY *body,CONTROL *control,FILES *files,\
-     OPTIONS *options, SYSTEM *system,int iFile) {
-  int lTmp=-1,bTmp;
-  AddOptionBool(files->Infile[iFile].cIn,options->cName,&bTmp,\
-     &lTmp,control->Io.iVerbose);
+void ReadOrbitOblData(BODY *body, CONTROL *control, FILES *files,
+                      OPTIONS *options, SYSTEM *system, int iFile) {
+  int lTmp = -1, bTmp;
+  AddOptionBool(files->Infile[iFile].cIn, options->cName, &bTmp, &lTmp,
+                control->Io.iVerbose);
   if (lTmp >= 0) {
-    NotPrimaryInput(iFile,options->cName,files->Infile[iFile].cIn,lTmp,\
-        control->Io.iVerbose);
+    NotPrimaryInput(iFile, options->cName, files->Infile[iFile].cIn, lTmp,
+                    control->Io.iVerbose);
     /* Option was found */
-    body[iFile-1].bReadOrbitOblData = bTmp;
-    UpdateFoundOption(&files->Infile[iFile],options,lTmp,iFile);
+    body[iFile - 1].bReadOrbitOblData = bTmp;
+    UpdateFoundOption(&files->Infile[iFile], options, lTmp, iFile);
   } else
-    body[iFile-1].bReadOrbitOblData = options->dDefault;
+    body[iFile - 1].bReadOrbitOblData = options->dDefault;
 }
 
-void ReadFileOrbitOblData(BODY *body,CONTROL *control,FILES *files,\
-     OPTIONS *options,SYSTEM *system,int iFile) {
-  int lTmp=-1;
+void ReadFileOrbitOblData(BODY *body, CONTROL *control, FILES *files,
+                          OPTIONS *options, SYSTEM *system, int iFile) {
+  int lTmp = -1;
   char cTmp[OPTLEN];
 
-  AddOptionString(files->Infile[iFile].cIn,options->cName,cTmp,&lTmp,\
-     control->Io.iVerbose);
+  AddOptionString(files->Infile[iFile].cIn, options->cName, cTmp, &lTmp,
+                  control->Io.iVerbose);
   if (lTmp >= 0) {
     /* Cannot exist in primary input file -- Each body has an output file */
-    NotPrimaryInput(iFile,options->cName,files->Infile[iFile].cIn,lTmp,\
-        control->Io.iVerbose);
-    strcpy(body[iFile-1].cFileOrbitOblData,cTmp);
-    UpdateFoundOption(&files->Infile[iFile],options,lTmp,iFile);
-  } else
-    if (iFile > 0)
-      strcpy(body[iFile-1].cFileOrbitOblData,options->cDefault);
+    NotPrimaryInput(iFile, options->cName, files->Infile[iFile].cIn, lTmp,
+                    control->Io.iVerbose);
+    strcpy(body[iFile - 1].cFileOrbitOblData, cTmp);
+    UpdateFoundOption(&files->Infile[iFile], options, lTmp, iFile);
+  } else if (iFile > 0)
+    strcpy(body[iFile - 1].cFileOrbitOblData, options->cDefault);
 }
 
 void ReadIceAlbedo(BODY *body, CONTROL *control, FILES *files, OPTIONS *options,
@@ -1132,29 +1131,29 @@ void InitializeOptionsPoise(OPTIONS *options, fnReadOption fnRead[]) {
   options[OPT_LATCELLNUM].dDefault   = 50;
   options[OPT_LATCELLNUM].iType      = 1;
   options[OPT_LATCELLNUM].bMultiFile = 1;
-  fnRead[OPT_LATCELLNUM] = &ReadLatCellNum;
+  fnRead[OPT_LATCELLNUM]             = &ReadLatCellNum;
 
-  sprintf(options[OPT_READORBITOBLDATA].cName,"bReadOrbitOblData");
-  sprintf(options[OPT_READORBITOBLDATA].cDescr,"Read in orbital and obliquity \
+  sprintf(options[OPT_READORBITOBLDATA].cName, "bReadOrbitOblData");
+  sprintf(options[OPT_READORBITOBLDATA].cDescr, "Read in orbital and obliquity \
     data and use with poise");
-  sprintf(options[OPT_READORBITOBLDATA].cDefault,"0");
-  options[OPT_READORBITOBLDATA].dDefault = 0;
-  options[OPT_READORBITOBLDATA].iType = 0;
+  sprintf(options[OPT_READORBITOBLDATA].cDefault, "0");
+  options[OPT_READORBITOBLDATA].dDefault   = 0;
+  options[OPT_READORBITOBLDATA].iType      = 0;
   options[OPT_READORBITOBLDATA].bMultiFile = 1;
-  fnRead[OPT_READORBITOBLDATA] = &ReadOrbitOblData;
+  fnRead[OPT_READORBITOBLDATA]             = &ReadOrbitOblData;
 
-  sprintf(options[OPT_FILEORBITOBLDATA].cName,"sFileOrbitOblData");
-  sprintf(options[OPT_FILEORBITOBLDATA].cDescr,"Name of file containing orbit\
+  sprintf(options[OPT_FILEORBITOBLDATA].cName, "sFileOrbitOblData");
+  sprintf(options[OPT_FILEORBITOBLDATA].cDescr, "Name of file containing orbit\
     and obliquity time series");
-  sprintf(options[OPT_FILEORBITOBLDATA].cDefault,"Obl_data.txt");
+  sprintf(options[OPT_FILEORBITOBLDATA].cDefault, "Obl_data.txt");
   options[OPT_FILEORBITOBLDATA].iType = 3;
-  fnRead[OPT_FILEORBITOBLDATA] = &ReadFileOrbitOblData;
+  fnRead[OPT_FILEORBITOBLDATA]        = &ReadFileOrbitOblData;
 
-  sprintf(options[OPT_PLANCKA].cName,"dPlanckA");
-  sprintf(options[OPT_PLANCKA].cDescr,"Constant 'A' used in OLR calculation");
-  sprintf(options[OPT_PLANCKA].cDefault,"203.3");
-  options[OPT_PLANCKA].dDefault = 203.3;
-  options[OPT_PLANCKA].iType = 2;
+  sprintf(options[OPT_PLANCKA].cName, "dPlanckA");
+  sprintf(options[OPT_PLANCKA].cDescr, "Constant 'A' used in OLR calculation");
+  sprintf(options[OPT_PLANCKA].cDefault, "203.3");
+  options[OPT_PLANCKA].dDefault   = 203.3;
+  options[OPT_PLANCKA].iType      = 2;
   options[OPT_PLANCKA].bMultiFile = 1;
   fnRead[OPT_PLANCKA]             = &ReadPlanckA;
 
@@ -1623,12 +1622,14 @@ void InitializeOptionsPoise(OPTIONS *options, fnReadOption fnRead[]) {
   sprintf(options[OPT_MINICEHEIGHT].cNeg, "meters");
   fnRead[OPT_MINICEHEIGHT] = &ReadMinIceSheetHeight;
   sprintf(options[OPT_MINICEHEIGHT].cLongDescr,
-    "The minimum thickness of permanent ice in a latitude bin for it to be\n"
-    "labeled ice-covered. In some cases, such as rapid thawing, a latituden\n"
-    "can have a very low value of ice height, primarily for numerical\n"
-    "reasons. Parameter forces unphysically small values of the ice height\n"
-    "to be ignored.\n"
-  );
+          "The minimum thickness of permanent ice in a latitude bin for it to "
+          "be\n"
+          "labeled ice-covered. In some cases, such as rapid thawing, a "
+          "latituden\n"
+          "can have a very low value of ice height, primarily for numerical\n"
+          "reasons. Parameter forces unphysically small values of the ice "
+          "height\n"
+          "to be ignored.\n");
 }
 
 void ReadOptionsPoise(BODY *body, CONTROL *control, FILES *files,
@@ -1782,99 +1783,109 @@ void VerifyOLR(BODY *body, OPTIONS *options, char cFile[], int iBody,
   }
 }
 
-void VerifyOrbitOblData(BODY *body,CONTROL *control,OPTIONS *options,\
-       int iBody){
+void VerifyOrbitOblData(BODY *body, CONTROL *control, OPTIONS *options,
+                        int iBody) {
   int iNLines, iLine, c;
-  double dttmp, datmp, detmp, ditmp, daptmp, dlatmp, dmatmp, dobltmp,dprecatmp;
+  double dttmp, datmp, detmp, ditmp, daptmp, dlatmp, dmatmp, dobltmp, dprecatmp;
   FILE *fileorb;
 
   if (body[iBody].bReadOrbitOblData) {
-    if (options[OPT_FILEORBITOBLDATA].iLine[iBody+1] == -1) {
-      fprintf(stderr,"ERROR: Must set %s if using %s for file %s\n",\
-        options[OPT_FILEORBITOBLDATA].cName,\
-        options[OPT_READORBITOBLDATA].cName,body[iBody].cName);
+    if (options[OPT_FILEORBITOBLDATA].iLine[iBody + 1] == -1) {
+      fprintf(stderr, "ERROR: Must set %s if using %s for file %s\n",
+              options[OPT_FILEORBITOBLDATA].cName,
+              options[OPT_READORBITOBLDATA].cName, body[iBody].cName);
       exit(EXIT_INPUT);
     } else {
-      fileorb = fopen(body[iBody].cFileOrbitOblData,"r");
+      fileorb = fopen(body[iBody].cFileOrbitOblData, "r");
       if (fileorb == NULL) {
         printf("ERROR: File %s not found.\n", body[iBody].cFileOrbitOblData);
         exit(EXIT_INPUT);
       }
       iNLines = 0;
       while ((c = getc(fileorb)) != EOF) {
-        if (c == '\n') iNLines++;              //add 1 for each new line
+        if (c == '\n')
+          iNLines++; // add 1 for each new line
       }
       rewind(fileorb);
 
-      body[iBody].iNLines = iNLines;
-      body[iBody].daTimeSeries = malloc(iNLines*sizeof(double));
-      body[iBody].daSemiSeries = malloc(iNLines*sizeof(double));
-      body[iBody].daEccSeries = malloc(iNLines*sizeof(double));
-      body[iBody].daArgPSeries = malloc(iNLines*sizeof(double));
-      body[iBody].daLongASeries = malloc(iNLines*sizeof(double));
-      body[iBody].daOblSeries = malloc(iNLines*sizeof(double));
-      body[iBody].daPrecASeries = malloc(iNLines*sizeof(double));
+      body[iBody].iNLines       = iNLines;
+      body[iBody].daTimeSeries  = malloc(iNLines * sizeof(double));
+      body[iBody].daSemiSeries  = malloc(iNLines * sizeof(double));
+      body[iBody].daEccSeries   = malloc(iNLines * sizeof(double));
+      body[iBody].daArgPSeries  = malloc(iNLines * sizeof(double));
+      body[iBody].daLongASeries = malloc(iNLines * sizeof(double));
+      body[iBody].daOblSeries   = malloc(iNLines * sizeof(double));
+      body[iBody].daPrecASeries = malloc(iNLines * sizeof(double));
 
-      body[iBody].daHeccSeries = malloc(iNLines*sizeof(double));
-      body[iBody].daKeccSeries = malloc(iNLines*sizeof(double));
+      body[iBody].daHeccSeries = malloc(iNLines * sizeof(double));
+      body[iBody].daKeccSeries = malloc(iNLines * sizeof(double));
 
       printf("file open\n");
 
       iLine = 0;
       while (feof(fileorb) == 0) {
-        fscanf(fileorb, "%lf %lf %lf %lf %lf %lf %lf", &dttmp, &datmp, &detmp,\
-          &daptmp, &dlatmp, &dobltmp, &dprecatmp);
+        fscanf(fileorb, "%lf %lf %lf %lf %lf %lf %lf", &dttmp, &datmp, &detmp,
+               &daptmp, &dlatmp, &dobltmp, &dprecatmp);
 
-        body[iBody].daTimeSeries[iLine] = \
-          dttmp*fdUnitsTime(control->Units[iBody+1].iTime);
-        body[iBody].daSemiSeries[iLine] = \
-          datmp*fdUnitsLength(control->Units[iBody+1].iLength);
+        body[iBody].daTimeSeries[iLine] =
+              dttmp * fdUnitsTime(control->Units[iBody + 1].iTime);
+        body[iBody].daSemiSeries[iLine] =
+              datmp * fdUnitsLength(control->Units[iBody + 1].iLength);
         body[iBody].daEccSeries[iLine] = detmp;
 
-        if (control->Units[iBody+1].iAngle == 0) {
-          body[iBody].daArgPSeries[iLine] = daptmp;
+        if (control->Units[iBody + 1].iAngle == 0) {
+          body[iBody].daArgPSeries[iLine]  = daptmp;
           body[iBody].daLongASeries[iLine] = dlatmp;
-          body[iBody].daOblSeries[iLine] = dobltmp;
+          body[iBody].daOblSeries[iLine]   = dobltmp;
           body[iBody].daPrecASeries[iLine] = dprecatmp;
         } else {
-          body[iBody].daArgPSeries[iLine] = daptmp*DEGRAD;
-          body[iBody].daLongASeries[iLine] = dlatmp*DEGRAD;
-          body[iBody].daOblSeries[iLine] = dobltmp*DEGRAD;
-          body[iBody].daPrecASeries[iLine] = dprecatmp*DEGRAD;
+          body[iBody].daArgPSeries[iLine]  = daptmp * DEGRAD;
+          body[iBody].daLongASeries[iLine] = dlatmp * DEGRAD;
+          body[iBody].daOblSeries[iLine]   = dobltmp * DEGRAD;
+          body[iBody].daPrecASeries[iLine] = dprecatmp * DEGRAD;
         }
-        body[iBody].daHeccSeries[iLine] = body[iBody].daEccSeries[iLine]*\
-            sin(body[iBody].daArgPSeries[iLine]+\
-            body[iBody].daLongASeries[iLine]);
-        body[iBody].daKeccSeries[iLine] = body[iBody].daEccSeries[iLine]*\
-            cos(body[iBody].daArgPSeries[iLine]+\
-            body[iBody].daLongASeries[iLine]);
+        body[iBody].daHeccSeries[iLine] = body[iBody].daEccSeries[iLine] *
+                                          sin(body[iBody].daArgPSeries[iLine] +
+                                              body[iBody].daLongASeries[iLine]);
+        body[iBody].daKeccSeries[iLine] = body[iBody].daEccSeries[iLine] *
+                                          cos(body[iBody].daArgPSeries[iLine] +
+                                              body[iBody].daLongASeries[iLine]);
 
         iLine++;
       }
       fclose(fileorb);
-
     }
     body[iBody].iCurrentStep = 0;
     if (control->Evolve.bVarDt) {
-      fprintf(stderr,"ERROR: Cannot use variable time step (%s = 1) if %s = 1\n",options[OPT_VARDT].cName,options[OPT_READORBITDATA].cName);
+      fprintf(stderr,
+              "ERROR: Cannot use variable time step (%s = 1) if %s = 1\n",
+              options[OPT_VARDT].cName, options[OPT_READORBITDATA].cName);
       exit(EXIT_INPUT);
     }
     if (control->Evolve.bDoForward) {
       if (body[iBody].daTimeSeries[1] != control->Evolve.dTimeStep) {
-        fprintf(stderr,"ERROR: Time step size (%s = 1) must match orbital data if %s = 1\n",options[OPT_TIMESTEP].cName,options[OPT_READORBITDATA].cName);
+        fprintf(stderr,
+                "ERROR: Time step size (%s = 1) must match orbital data if %s "
+                "= 1\n",
+                options[OPT_TIMESTEP].cName, options[OPT_READORBITDATA].cName);
         exit(EXIT_INPUT);
       }
     } else if (control->Evolve.bDoBackward) {
-      if (body[iBody].daTimeSeries[1] != -1*control->Evolve.dTimeStep) {
-        fprintf(stderr,"ERROR: Time step size (%s = 1) must match orbital data if %s = 1\n",options[OPT_TIMESTEP].cName,options[OPT_READORBITDATA].cName);
+      if (body[iBody].daTimeSeries[1] != -1 * control->Evolve.dTimeStep) {
+        fprintf(stderr,
+                "ERROR: Time step size (%s = 1) must match orbital data if %s "
+                "= 1\n",
+                options[OPT_TIMESTEP].cName, options[OPT_READORBITDATA].cName);
         exit(EXIT_INPUT);
       }
     }
-    if (iNLines < (control->Evolve.dStopTime/control->Evolve.dTimeStep+1) ) {
-      fprintf(stderr,"ERROR: Input orbit data must at least as long as vplanet integration (%f years)\n",control->Evolve.dStopTime/YEARSEC);
+    if (iNLines < (control->Evolve.dStopTime / control->Evolve.dTimeStep + 1)) {
+      fprintf(stderr,
+              "ERROR: Input orbit data must at least as long as vplanet "
+              "integration (%f years)\n",
+              control->Evolve.dStopTime / YEARSEC);
       exit(EXIT_INPUT);
     }
-
   }
 }
 
@@ -4500,15 +4511,15 @@ void AddModulePoise(CONTROL *control, MODULE *module, int iBody, int iModule) {
 
 /************* POISE Functions ***********/
 void UpdateOrbitOblData(BODY *body, EVOLVE *evolve, int iBody) {
-  body[iBody].dSemi = body[iBody].daSemiSeries[body[iBody].iCurrentStep];
-  body[iBody].dEcc = body[iBody].daEccSeries[body[iBody].iCurrentStep];
-  body[iBody].dArgP = body[iBody].daArgPSeries[body[iBody].iCurrentStep];
-  body[iBody].dLongA = body[iBody].daLongASeries[body[iBody].iCurrentStep];
+  body[iBody].dSemi      = body[iBody].daSemiSeries[body[iBody].iCurrentStep];
+  body[iBody].dEcc       = body[iBody].daEccSeries[body[iBody].iCurrentStep];
+  body[iBody].dArgP      = body[iBody].daArgPSeries[body[iBody].iCurrentStep];
+  body[iBody].dLongA     = body[iBody].daLongASeries[body[iBody].iCurrentStep];
   body[iBody].dObliquity = body[iBody].daOblSeries[body[iBody].iCurrentStep];
-  body[iBody].dPrecA = body[iBody].daPrecASeries[body[iBody].iCurrentStep];
+  body[iBody].dPrecA     = body[iBody].daPrecASeries[body[iBody].iCurrentStep];
 
-  body[iBody].dXobl = sin(body[iBody].dObliquity)*cos(body[iBody].dPrecA);
-  body[iBody].dYobl = sin(body[iBody].dObliquity)*sin(body[iBody].dPrecA);
+  body[iBody].dXobl = sin(body[iBody].dObliquity) * cos(body[iBody].dPrecA);
+  body[iBody].dYobl = sin(body[iBody].dObliquity) * sin(body[iBody].dPrecA);
   body[iBody].dZobl = cos(body[iBody].dObliquity);
 
   body[iBody].dHecc = body[iBody].daHeccSeries[body[iBody].iCurrentStep];
@@ -6964,6 +6975,323 @@ void EnergyResiduals(BODY *body, int iBody, int iNday) {
   }
 }
 
+void fvPoiseSeasonalInitialize(BODY *body, int iBody, int iYear) {
+  body[iBody].dTGlobal       = 0.0;
+  body[iBody].dFluxInGlobal  = 0.0;
+  body[iBody].dFluxOutGlobal = 0.0;
+
+  for (int iLat = 0; iLat < body[iBody].iNumLats; iLat++) {
+    // start of year, reset annual averages to zero
+    body[iBody].daTempAvg[iLat]          = 0.0;
+    body[iBody].daPlanckBAvg[iLat]       = 0.0;
+    body[iBody].daAlbedoAvg[iLat]        = 0.0;
+    body[iBody].daTempAvgL[iLat]         = 0.0;
+    body[iBody].daTempAvgW[iLat]         = 0.0;
+    body[iBody].daAlbedoAvgL[iLat]       = 0.0;
+    body[iBody].daAlbedoAvgW[iLat]       = 0.0;
+    body[iBody].daFluxAvg[iLat]          = 0.0;
+    body[iBody].daFluxInAvg[iLat]        = 0.0;
+    body[iBody].daDivFluxAvg[iLat]       = 0.0;
+    body[iBody].daFluxOutAvg[iLat]       = 0.0;
+    body[iBody].daIceBalanceAnnual[iLat] = 0.0;
+    body[iBody].daEnerResLAnn[iLat]      = 0.0;
+    body[iBody].daEnerResWAnn[iLat]      = 0.0;
+    body[iBody].daTempMinLW[iLat]        = dHUGE;
+    body[iBody].daTempMaxLW[iLat]        = -1 * dHUGE;
+    body[iBody].daTempMaxLand[iLat]      = -1 * dHUGE;
+    body[iBody].daTempMaxWater[iLat]     = -1 * dHUGE;
+
+    if (iYear == 0) {
+      // reset ice sheet stuff only on first year
+      if (body[iBody].bIceSheets) {
+        body[iBody].daIceMassTmp[iLat]   = body[iBody].daIceMass[iLat];
+        body[iBody].daIceAccumTot[iLat]  = 0;
+        body[iBody].daIceAblateTot[iLat] = 0;
+      }
+    }
+  }
+}
+
+void fvCalcPlanckAB(BODY *body, int iBody, int iLat) {
+  body[iBody].daTempLW[iLat] =
+        body[iBody].daLandFrac[iLat] * body[iBody].daTempLand[iLat] +
+        body[iBody].daWaterFrac[iLat] * body[iBody].daTempWater[iLat];
+  body[iBody].dTGlobalTmp += body[iBody].daTempLW[iLat] / body[iBody].iNumLats;
+
+  if (body[iBody].bCalcAB) {
+    if (body[iBody].iOLRModel == WK97) {
+      /* Calculate A and B from williams and kasting 97 result */
+      body[iBody].daPlanckBSea[iLat] = fdOLRdTwk97(body, iBody, iLat, SEA);
+      body[iBody].daPlanckASea[iLat] =
+            fdOLRwk97(body, iBody, iLat, SEA) -
+            body[iBody].daPlanckBSea[iLat] * (body[iBody].daTempLW[iLat]);
+    } else if (body[iBody].iOLRModel == HM16) {
+      /* Calculate A and B from haqq-misra+ 2016 result */
+      body[iBody].daPlanckBSea[iLat] = fdOLRdThm16(body, iBody, iLat, SEA);
+      body[iBody].daPlanckASea[iLat] =
+            fdOLRhm16(body, iBody, iLat, SEA) -
+            body[iBody].daPlanckBSea[iLat] * (body[iBody].daTempLW[iLat]);
+    } else {
+      /* Calculate A and B from spiegel+ 2009 model */
+      body[iBody].daPlanckBSea[iLat] = fdOLRdTsms09(body, iBody, iLat, SEA);
+      body[iBody].daPlanckASea[iLat] =
+            fdOLRsms09(body, iBody, iLat, SEA) -
+            body[iBody].daPlanckBSea[iLat] * (body[iBody].daTempLW[iLat]);
+    }
+
+    // XXX This section doesn't seem releated to calculating A and B,
+    // but was part of this if-else originally
+
+    if (body[iBody].bMEPDiff) {
+      if (iLat == 0) {
+        body[iBody].daDiffusionSea[iLat] = body[iBody].daPlanckBSea[iLat] / 4.0;
+
+      } else if (iLat == (body[iBody].iNumLats - 1)) {
+        body[iBody].daDiffusionSea[iLat] =
+              (body[iBody].daPlanckBSea[iLat] +
+               body[iBody].daPlanckBSea[iLat - 1]) /
+              8.0;
+
+        body[iBody].daDiffusionSea[iLat + 1] =
+              body[iBody].daPlanckBSea[iLat] / 4.0;
+
+      } else {
+        body[iBody].daDiffusionSea[iLat] =
+              (body[iBody].daPlanckBSea[iLat] +
+               body[iBody].daPlanckBSea[iLat - 1]) /
+              8.0;
+      }
+    }
+  }
+}
+
+void fvFluxesByLatitude(BODY *body, int iBody, int iNyear, int iNstep) {
+  for (int iLat = 0; iLat < body[iBody].iNumLats; iLat++) {
+    body[iBody].daDMidPt[iLat] = 0.5 * (body[iBody].daDiffusionSea[iLat + 1] +
+                                        body[iBody].daDiffusionSea[iLat]);
+    body[iBody].daFlux[iLat] =
+          -2. * PI * (body[iBody].dRadius * body[iBody].dRadius) *
+          sqrt(1.0 - (sin(body[iBody].daLats[iLat]) *
+                      sin(body[iBody].daLats[iLat]))) *
+          body[iBody].daDMidPt[iLat] * body[iBody].daTGrad[iLat];
+    body[iBody].daFluxAvg[iLat] +=
+          body[iBody].daFlux[iLat] / body[iBody].iNStepInYear;
+    body[iBody].daFluxDaily[iLat][iNyear * body[iBody].iNStepInYear + iNstep] =
+          body[iBody].daFlux[iLat];
+
+    body[iBody].daDivFlux[iLat] = 0.0;
+    for (int jLat = 0; jLat < body[iBody].iNumLats; jLat++) {
+      body[iBody].daDivFlux[iLat] +=
+            -body[iBody].daMDiffSea[iLat][jLat] * body[iBody].daTempLW[jLat];
+    }
+    body[iBody].daDivFluxAvg[iLat] +=
+          body[iBody].daDivFlux[iLat] / body[iBody].iNStepInYear;
+    body[iBody]
+          .daDivFluxDaily[iLat][iNyear * body[iBody].iNStepInYear + iNstep] =
+          body[iBody].daDivFlux[iLat];
+  }
+}
+
+// annual averages by latitude
+void fvPoiseAnnualAveragesByLatitude(BODY *body, int iBody, int iLat) {
+  body[iBody].daTempAvg[iLat] +=
+        body[iBody].daTempLW[iLat] / body[iBody].iNStepInYear;
+  body[iBody].daPlanckBAvg[iLat] +=
+        body[iBody].daPlanckBSea[iLat] / body[iBody].iNStepInYear;
+  body[iBody].daTempAvgL[iLat] +=
+        body[iBody].daTempLand[iLat] / body[iBody].iNStepInYear;
+  body[iBody].daTempAvgW[iLat] +=
+        body[iBody].daTempWater[iLat] / body[iBody].iNStepInYear;
+  body[iBody].daAlbedoAvgL[iLat] +=
+        body[iBody].daAlbedoLand[iLat] / body[iBody].iNStepInYear;
+  body[iBody].daAlbedoAvgW[iLat] +=
+        body[iBody].daAlbedoWater[iLat] / body[iBody].iNStepInYear;
+  body[iBody].daAlbedoAvg[iLat] +=
+        body[iBody].daAlbedoLW[iLat] / body[iBody].iNStepInYear;
+  body[iBody].daFluxInAvg[iLat] +=
+        body[iBody].daFluxIn[iLat] / body[iBody].iNStepInYear;
+  body[iBody].daFluxOutAvg[iLat] +=
+        body[iBody].daFluxOut[iLat] / body[iBody].iNStepInYear;
+}
+
+void fvPoiseDailyProps(BODY *body, int iBody, int iLat, int iNyear,
+                       int iNstep) {
+  body[iBody].daTempDaily[iLat][iNyear * body[iBody].iNStepInYear + iNstep] =
+        body[iBody].daTempLW[iLat];
+  body[iBody].daFluxInDaily[iLat][iNyear * body[iBody].iNStepInYear + iNstep] =
+        body[iBody].daFluxIn[iLat];
+  body[iBody].daFluxOutDaily[iLat][iNyear * body[iBody].iNStepInYear + iNstep] =
+        body[iBody].daFluxOut[iLat];
+  body[iBody].daPlanckBDaily[iLat][iNyear * body[iBody].iNStepInYear + iNstep] =
+        body[iBody].daPlanckBSea[iLat];
+}
+
+void fvPoiseSeasonalMinMax(BODY *body, int iBody, int iLat) {
+  if (body[iBody].daTempLW[iLat] < body[iBody].daTempMinLW[iLat]) {
+    body[iBody].daTempMinLW[iLat] = body[iBody].daTempLW[iLat];
+  }
+  if (body[iBody].daTempLW[iLat] > body[iBody].daTempMaxLW[iLat]) {
+    body[iBody].daTempMaxLW[iLat] = body[iBody].daTempLW[iLat];
+  }
+  if (body[iBody].daTempLand[iLat] > body[iBody].daTempMaxLand[iLat]) {
+    body[iBody].daTempMaxLand[iLat] = body[iBody].daTempLand[iLat];
+  }
+  if (body[iBody].daTempWater[iLat] > body[iBody].daTempMaxWater[iLat]) {
+    body[iBody].daTempMaxWater[iLat] = body[iBody].daTempWater[iLat];
+  }
+}
+
+void fvCalculateIceSheets(BODY *body, double dStepsize, int iBody, int iLat,
+                          int iNstep) {
+  // calculate derivative of ice mass density and take an euler step
+  body[iBody].daIceBalance[iLat][iNstep] = fdIceMassBalance(body, iBody, iLat);
+
+  body[iBody].daIceMassTmp[iLat] +=
+        dStepsize * body[iBody].daIceBalance[iLat][iNstep];
+
+  if (body[iBody].daIceBalance[iLat][iNstep] >= 0) {
+    body[iBody].daIceAccumTot[iLat] += dStepsize *
+                                       body[iBody].daIceBalance[iLat][iNstep] /
+                                       body[iBody].iNumYears;
+
+  } else {
+    if (body[iBody].daIceMassTmp[iLat] > 0) {
+      if (body[iBody].daIceMassTmp[iLat] >=
+          dStepsize * body[iBody].daIceBalance[iLat][iNstep]) {
+
+        body[iBody].daIceAblateTot[iLat] +=
+              dStepsize * body[iBody].daIceBalance[iLat][iNstep] /
+              body[iBody].iNumYears;
+
+      } else {
+        body[iBody].daIceAblateTot[iLat] += body[iBody].daIceMassTmp[iLat];
+      }
+    }
+  }
+  if (body[iBody].daIceMassTmp[iLat] < 0.0) {
+    body[iBody].daIceMassTmp[iLat] = 0.0;
+  } // don't let ice mass become negative
+  if (body[iBody].daIceBalance[iLat][iNstep] < 0 &&
+      body[iBody].daIceMassTmp[iLat] != 0) {
+
+    if (body[iBody].daIceMassTmp[iLat] <=
+        fabs(dStepsize * body[iBody].daIceBalance[iLat][iNstep])) {
+
+      // adjust temperature
+      body[iBody].daTempLand[iLat] +=
+            -body[iBody].daIceMassTmp[iLat] * LFICE / body[iBody].dHeatCapLand;
+
+    } else {
+      body[iBody].daTempLand[iLat] +=
+            dStepsize * body[iBody].daIceBalance[iLat][iNstep] * LFICE /
+            body[iBody].dHeatCapLand; // adjust temperature
+    }
+  } else if (body[iBody].daIceBalance[iLat][iNstep] > 0) {
+    body[iBody].daTempLand[iLat] +=
+          dStepsize * body[iBody].daIceBalance[iLat][iNstep] * LFICE /
+          body[iBody].dHeatCapLand; // adjust temperature
+  }
+}
+
+void fvCalculateFluxesByLatitude(BODY *body, int iBody, int iLat, int iNday) {
+  body[iBody].daFluxOutLand[iLat] =
+        body[iBody].daPlanckASea[iLat] +
+        body[iBody].daPlanckBSea[iLat] * body[iBody].daTempLand[iLat];
+  body[iBody].daFluxOutWater[iLat] =
+        body[iBody].daPlanckASea[iLat] +
+        body[iBody].daPlanckBSea[iLat] * body[iBody].daTempWater[iLat];
+  body[iBody].daFluxOut[iLat] =
+        body[iBody].daLandFrac[iLat] * body[iBody].daFluxOutLand[iLat] +
+        body[iBody].daWaterFrac[iLat] * body[iBody].daFluxOutWater[iLat];
+  body[iBody].dFluxOutGlobalTmp +=
+        body[iBody].daFluxOut[iLat] / body[iBody].iNumLats;
+
+  body[iBody].daFluxInLand[iLat] = (1.0 - body[iBody].daAlbedoLand[iLat]) *
+                                   body[iBody].daInsol[iLat][iNday];
+  body[iBody].daFluxInWater[iLat] = (1.0 - body[iBody].daAlbedoWater[iLat]) *
+                                    body[iBody].daInsol[iLat][iNday];
+  body[iBody].daFluxIn[iLat] =
+        body[iBody].daLandFrac[iLat] * body[iBody].daFluxInLand[iLat] +
+        body[iBody].daWaterFrac[iLat] * body[iBody].daFluxInWater[iLat];
+  body[iBody].dFluxInGlobalTmp +=
+        body[iBody].daFluxIn[iLat] / body[iBody].iNumLats;
+}
+
+void fvPoiseGlobalProps(BODY *body, int iBody) {
+  body[iBody].dTGlobal += body[iBody].dTGlobalTmp / (body[iBody].iNStepInYear);
+  body[iBody].dFluxOutGlobal +=
+        body[iBody].dFluxOutGlobalTmp / (body[iBody].iNStepInYear);
+  body[iBody].dFluxInGlobal +=
+        body[iBody].dFluxInGlobalTmp / (body[iBody].iNStepInYear);
+}
+
+void fvFinishSeaIce(BODY *body, double dStepsize, int iBody, int iLat,
+                    int iNyear, int iNday, int iNstep) {
+  fvSnowball(body, iBody);
+  // ice growth/ablation
+  if (body[iBody].bIceSheets) {
+    fvCalculateIceSheets(body, dStepsize, iBody, iLat, iNstep);
+  }
+
+  fvCalcPlanckAB(body, iBody, iLat);
+  fvCalculateFluxesByLatitude(body, iBody, iLat, iNday);
+  fvPoiseAnnualAveragesByLatitude(body, iBody, iLat);
+  fvPoiseDailyProps(body, iBody, iLat, iNyear, iNstep);
+  fvPoiseSeasonalMinMax(body, iBody, iLat);
+}
+
+void fvCalculateSeaIce(BODY *body, double dStepsize, int iBody, int iNyear,
+                       int iNday, int iNstep) {
+  fvSeaIce(body, iBody);
+  for (int iLat = 0; iLat < body[iBody].iNumLats; iLat++) {
+    if (body[iBody].daSeaIceHeight[iLat] > 0) {
+      // adjust height of present sea ice
+      body[iBody].daSeaIceHeight[iLat] -= body[iBody].dSeasDeltat /
+                                          body[iBody].dLatentHeatIce *
+                                          body[iBody].daFluxSeaIce[iLat];
+
+      if (body[iBody].daSeaIceHeight[iLat] < 0) {
+        body[iBody].daSeaIceHeight[iLat] = 0;
+      }
+    } else {
+      if (body[iBody].daTempWater[iLat] < body[iBody].dFrzTSeaIce) {
+        // create new sea ice
+        body[iBody].daSeaIceHeight[iLat] =
+              -body[iBody].dHeatCapWater * body[iBody].dMeanMotion / (2 * PI) /
+              body[iBody].dLatentHeatIce *
+              (body[iBody].daTempWater[iLat] - body[iBody].dFrzTSeaIce);
+      }
+    }
+
+    fvFinishSeaIce(body, dStepsize, iBody, iLat, iNyear, iNday, iNstep);
+  }
+}
+
+void fvStaticSeaIce(BODY *body, double dStepsize, int iBody, int iNyear,
+                    int iNday, int iNstep) {
+  int iLat, jLat;
+
+  for (iLat = 0; iLat < 2 * body[iBody].iNumLats; iLat++) {
+    body[iBody].daTmpTempSea[iLat] = 0.0;
+    for (jLat = 0; jLat < 2 * body[iBody].iNumLats; jLat++) {
+      body[iBody].daTmpTempSea[iLat] +=
+            body[iBody].daInvMSea[iLat][jLat] * body[iBody].daSourceLW[jLat];
+    }
+  }
+  for (iLat = 0; iLat < body[iBody].iNumLats; iLat++) {
+    // temp change this time step
+    body[iBody].daDeltaTempL[iLat] =
+          body[iBody].daTmpTempSea[2 * iLat] - body[iBody].daTempLand[iLat];
+    body[iBody].daDeltaTempW[iLat] = body[iBody].daTmpTempSea[2 * iLat + 1] -
+                                     body[iBody].daTempWater[iLat];
+    // calculate temperature and fluxes by latitude and global average
+    body[iBody].daTempLand[iLat]  = body[iBody].daTmpTempSea[2 * iLat];
+    body[iBody].daTempWater[iLat] = body[iBody].daTmpTempSea[2 * iLat + 1];
+
+    fvFinishSeaIce(body, dStepsize, iBody, iLat, iNyear, iNday, iNstep);
+  }
+}
+
 /**
 Runs the seasonal EBM, called from ForceBehavior. The seasonal EBM runs on a
 fixed time-step (of order days) to resolve the seasonal cycle. It can be said
@@ -6979,43 +7307,12 @@ void PoiseSeasonal(BODY *body, int iBody) {
   dStepsize = 2 * PI / body[iBody].dMeanMotion / body[iBody].iNStepInYear;
 
   // fvMatrixSeasonal(body,iBody);
-  /* main loop */
   // XXX confusing! Should distinguish b/w years/orbits
-  for (iNyear = 0; iNyear < body[iBody].iNumYears; iNyear++) {
-    body[iBody].dTGlobal       = 0.0;
-    body[iBody].dFluxInGlobal  = 0.0;
-    body[iBody].dFluxOutGlobal = 0.0;
+  // XXX Needs to be broken down into functions!
 
-    for (iLat = 0; iLat < body[iBody].iNumLats; iLat++) {
-      // start of year, reset annual averages to zero
-      body[iBody].daTempAvg[iLat]          = 0.0;
-      body[iBody].daPlanckBAvg[iLat]       = 0.0;
-      body[iBody].daAlbedoAvg[iLat]        = 0.0;
-      body[iBody].daTempAvgL[iLat]         = 0.0;
-      body[iBody].daTempAvgW[iLat]         = 0.0;
-      body[iBody].daAlbedoAvgL[iLat]       = 0.0;
-      body[iBody].daAlbedoAvgW[iLat]       = 0.0;
-      body[iBody].daFluxAvg[iLat]          = 0.0;
-      body[iBody].daFluxInAvg[iLat]        = 0.0;
-      body[iBody].daDivFluxAvg[iLat]       = 0.0;
-      body[iBody].daFluxOutAvg[iLat]       = 0.0;
-      body[iBody].daIceBalanceAnnual[iLat] = 0.0;
-      body[iBody].daEnerResLAnn[iLat]      = 0.0;
-      body[iBody].daEnerResWAnn[iLat]      = 0.0;
-      body[iBody].daTempMinLW[iLat]        = dHUGE;
-      body[iBody].daTempMaxLW[iLat]        = -1 * dHUGE;
-      body[iBody].daTempMaxLand[iLat]      = -1 * dHUGE;
-      body[iBody].daTempMaxWater[iLat]     = -1 * dHUGE;
-
-      if (iNyear == 0) {
-        // reset ice sheet stuff only on first year
-        if (body[iBody].bIceSheets) {
-          body[iBody].daIceMassTmp[iLat]   = body[iBody].daIceMass[iLat];
-          body[iBody].daIceAccumTot[iLat]  = 0;
-          body[iBody].daIceAblateTot[iLat] = 0;
-        }
-      }
-    }
+  /* main loop */
+  for (iNyear = 0; iNyear < body[iBody].iNumYears; iNyear++) {
+    fvPoiseSeasonalInitialize(body, iBody, iNyear);
 
     fvAlbedoSeasonal(body, iBody, 0);
     body[iBody].dAlbedoGlobal = 0.0;
@@ -7029,499 +7326,15 @@ void PoiseSeasonal(BODY *body, int iBody) {
       fvSourceFSeas(body, iBody, iNday);
 
       if (body[iBody].bSeaIceModel == 1) {
-        fvSeaIce(body, iBody);
-        for (iLat = 0; iLat < body[iBody].iNumLats; iLat++) {
-          if (body[iBody].daSeaIceHeight[iLat] > 0) {
-            // adjust height of present sea ice
-            body[iBody].daSeaIceHeight[iLat] -= body[iBody].dSeasDeltat /
-                                                body[iBody].dLatentHeatIce *
-                                                body[iBody].daFluxSeaIce[iLat];
-
-            if (body[iBody].daSeaIceHeight[iLat] < 0) {
-              body[iBody].daSeaIceHeight[iLat] = 0;
-            }
-
-          } else {
-
-            if (body[iBody].daTempWater[iLat] < body[iBody].dFrzTSeaIce) {
-              // create new sea ice
-              body[iBody].daSeaIceHeight[iLat] =
-                    -body[iBody].dHeatCapWater * body[iBody].dMeanMotion /
-                    (2 * PI) / body[iBody].dLatentHeatIce *
-                    (body[iBody].daTempWater[iLat] - body[iBody].dFrzTSeaIce);
-            }
-          }
-
-          fvSnowball(body, iBody);
-          // ice growth/ablation
-          if (body[iBody].bIceSheets) {
-            // calculate derivative of ice mass density and take an euler step
-            body[iBody].daIceBalance[iLat][iNstep] =
-                  fdIceMassBalance(body, iBody, iLat);
-
-            body[iBody].daIceMassTmp[iLat] +=
-                  dStepsize * body[iBody].daIceBalance[iLat][iNstep];
-
-            if (body[iBody].daIceBalance[iLat][iNstep] >= 0) {
-              body[iBody].daIceAccumTot[iLat] +=
-                    dStepsize * body[iBody].daIceBalance[iLat][iNstep] /
-                    body[iBody].iNumYears;
-
-            } else {
-              if (body[iBody].daIceMassTmp[iLat] > 0) {
-                if (body[iBody].daIceMassTmp[iLat] >=
-                    dStepsize * body[iBody].daIceBalance[iLat][iNstep]) {
-
-                  body[iBody].daIceAblateTot[iLat] +=
-                        dStepsize * body[iBody].daIceBalance[iLat][iNstep] /
-                        body[iBody].iNumYears;
-
-                } else {
-                  body[iBody].daIceAblateTot[iLat] +=
-                        body[iBody].daIceMassTmp[iLat];
-                }
-              }
-            }
-            if (body[iBody].daIceMassTmp[iLat] < 0.0) {
-              body[iBody].daIceMassTmp[iLat] = 0.0;
-            } // don't let ice mass become negative
-            if (body[iBody].daIceBalance[iLat][iNstep] < 0 &&
-                body[iBody].daIceMassTmp[iLat] != 0) {
-
-              if (body[iBody].daIceMassTmp[iLat] <=
-                  fabs(dStepsize * body[iBody].daIceBalance[iLat][iNstep])) {
-
-                // adjust temperature
-                body[iBody].daTempLand[iLat] +=
-                      -body[iBody].daIceMassTmp[iLat] * LFICE /
-                      body[iBody].dHeatCapLand;
-
-              } else {
-                body[iBody].daTempLand[iLat] +=
-                      dStepsize * body[iBody].daIceBalance[iLat][iNstep] *
-                      LFICE / body[iBody].dHeatCapLand; // adjust temperature
-              }
-            } else if (body[iBody].daIceBalance[iLat][iNstep] > 0) {
-              body[iBody].daTempLand[iLat] +=
-                    dStepsize * body[iBody].daIceBalance[iLat][iNstep] * LFICE /
-                    body[iBody].dHeatCapLand; // adjust temperature
-            }
-          }
-
-          body[iBody].daTempLW[iLat] =
-                body[iBody].daLandFrac[iLat] * body[iBody].daTempLand[iLat] +
-                body[iBody].daWaterFrac[iLat] * body[iBody].daTempWater[iLat];
-          body[iBody].dTGlobalTmp +=
-                body[iBody].daTempLW[iLat] / body[iBody].iNumLats;
-
-          if (body[iBody].bCalcAB) {
-            if (body[iBody].iOLRModel == WK97) {
-              /* Calculate A and B from williams and kasting 97 result */
-              body[iBody].daPlanckBSea[iLat] =
-                    fdOLRdTwk97(body, iBody, iLat, SEA);
-              body[iBody].daPlanckASea[iLat] =
-                    fdOLRwk97(body, iBody, iLat, SEA) -
-                    body[iBody].daPlanckBSea[iLat] *
-                          (body[iBody].daTempLW[iLat]);
-            } else if (body[iBody].iOLRModel == HM16) {
-              /* Calculate A and B from haqq-misra+ 2016 result */
-              body[iBody].daPlanckBSea[iLat] =
-                    fdOLRdThm16(body, iBody, iLat, SEA);
-              body[iBody].daPlanckASea[iLat] =
-                    fdOLRhm16(body, iBody, iLat, SEA) -
-                    body[iBody].daPlanckBSea[iLat] *
-                          (body[iBody].daTempLW[iLat]);
-            } else {
-              /* Calculate A and B from spiegel+ 2009 model */
-              body[iBody].daPlanckBSea[iLat] =
-                    fdOLRdTsms09(body, iBody, iLat, SEA);
-              body[iBody].daPlanckASea[iLat] =
-                    fdOLRsms09(body, iBody, iLat, SEA) -
-                    body[iBody].daPlanckBSea[iLat] *
-                          (body[iBody].daTempLW[iLat]);
-            }
-
-
-            if (body[iBody].bMEPDiff) {
-              if (iLat == 0) {
-                body[iBody].daDiffusionSea[iLat] =
-                      body[iBody].daPlanckBSea[iLat] / 4.0;
-
-              } else if (iLat == (body[iBody].iNumLats - 1)) {
-                body[iBody].daDiffusionSea[iLat] =
-                      (body[iBody].daPlanckBSea[iLat] +
-                       body[iBody].daPlanckBSea[iLat - 1]) /
-                      8.0;
-
-                body[iBody].daDiffusionSea[iLat + 1] =
-                      body[iBody].daPlanckBSea[iLat] / 4.0;
-
-              } else {
-                body[iBody].daDiffusionSea[iLat] =
-                      (body[iBody].daPlanckBSea[iLat] +
-                       body[iBody].daPlanckBSea[iLat - 1]) /
-                      8.0;
-              }
-            }
-          }
-
-          // calculate fluxes by latitude and global average
-          body[iBody].daFluxOutLand[iLat] =
-                body[iBody].daPlanckASea[iLat] +
-                body[iBody].daPlanckBSea[iLat] * body[iBody].daTempLand[iLat];
-          body[iBody].daFluxOutWater[iLat] =
-                body[iBody].daPlanckASea[iLat] +
-                body[iBody].daPlanckBSea[iLat] * body[iBody].daTempWater[iLat];
-          body[iBody].daFluxOut[iLat] =
-                body[iBody].daLandFrac[iLat] * body[iBody].daFluxOutLand[iLat] +
-                body[iBody].daWaterFrac[iLat] *
-                      body[iBody].daFluxOutWater[iLat];
-          body[iBody].dFluxOutGlobalTmp +=
-                body[iBody].daFluxOut[iLat] / body[iBody].iNumLats;
-          body[iBody].daFluxInLand[iLat] =
-                (1.0 - body[iBody].daAlbedoLand[iLat]) *
-                body[iBody].daInsol[iLat][iNday];
-          body[iBody].daFluxInWater[iLat] =
-                (1.0 - body[iBody].daAlbedoWater[iLat]) *
-                body[iBody].daInsol[iLat][iNday];
-          body[iBody].daFluxIn[iLat] =
-                body[iBody].daLandFrac[iLat] * body[iBody].daFluxInLand[iLat] +
-                body[iBody].daWaterFrac[iLat] * body[iBody].daFluxInWater[iLat];
-          body[iBody].dFluxInGlobalTmp +=
-                body[iBody].daFluxIn[iLat] / body[iBody].iNumLats;
-
-          // annual averages by latitude
-          body[iBody].daTempAvg[iLat] +=
-                body[iBody].daTempLW[iLat] / body[iBody].iNStepInYear;
-          body[iBody].daPlanckBAvg[iLat] +=
-                body[iBody].daPlanckBSea[iLat] / body[iBody].iNStepInYear;
-          body[iBody].daAlbedoAvg[iLat] +=
-                body[iBody].daAlbedoLW[iLat] / body[iBody].iNStepInYear;
-          body[iBody].daFluxInAvg[iLat] +=
-                body[iBody].daFluxIn[iLat] / body[iBody].iNStepInYear;
-          body[iBody].daFluxOutAvg[iLat] +=
-                body[iBody].daFluxOut[iLat] / body[iBody].iNStepInYear;
-
-          body[iBody]
-                .daTempDaily[iLat][iNyear * body[iBody].iNStepInYear + iNstep] =
-                body[iBody].daTempLW[iLat];
-          body[iBody].daTempAvgL[iLat] +=
-                body[iBody].daTempLand[iLat] / body[iBody].iNStepInYear;
-          body[iBody].daTempAvgW[iLat] +=
-                body[iBody].daTempWater[iLat] / body[iBody].iNStepInYear;
-          body[iBody].daAlbedoAvgL[iLat] +=
-                body[iBody].daAlbedoLand[iLat] / body[iBody].iNStepInYear;
-          body[iBody].daAlbedoAvgW[iLat] +=
-                body[iBody].daAlbedoWater[iLat] / body[iBody].iNStepInYear;
-          body[iBody].daAlbedoAvg[iLat] +=
-                body[iBody].daAlbedoLW[iLat] / body[iBody].iNStepInYear;
-          body[iBody].daFluxInAvg[iLat] +=
-                body[iBody].daFluxIn[iLat] / body[iBody].iNStepInYear;
-          body[iBody].daFluxOutAvg[iLat] +=
-                body[iBody].daFluxOut[iLat] / body[iBody].iNStepInYear;
-
-
-          // daily values for output over year
-          body[iBody]
-                .daTempDaily[iLat][iNyear * body[iBody].iNStepInYear + iNstep] =
-                body[iBody].daTempLW[iLat];
-
-          body[iBody].daFluxInDaily[iLat][iNyear * body[iBody].iNStepInYear +
-                                          iNstep] = body[iBody].daFluxIn[iLat];
-
-          body[iBody].daFluxOutDaily[iLat][iNyear * body[iBody].iNStepInYear +
-                                           iNstep] =
-                body[iBody].daFluxOut[iLat];
-
-          body[iBody].daPlanckBDaily[iLat][iNyear * body[iBody].iNStepInYear +
-                                           iNstep] =
-                body[iBody].daPlanckBSea[iLat];
-
-          if (body[iBody].daTempLW[iLat] < body[iBody].daTempMinLW[iLat]) {
-            body[iBody].daTempMinLW[iLat] = body[iBody].daTempLW[iLat];
-          }
-          if (body[iBody].daTempLW[iLat] > body[iBody].daTempMaxLW[iLat]) {
-            body[iBody].daTempMaxLW[iLat] = body[iBody].daTempLW[iLat];
-          }
-          if (body[iBody].daTempLand[iLat] > body[iBody].daTempMaxLand[iLat]) {
-            body[iBody].daTempMaxLand[iLat] = body[iBody].daTempLand[iLat];
-          }
-          if (body[iBody].daTempWater[iLat] >
-              body[iBody].daTempMaxWater[iLat]) {
-            body[iBody].daTempMaxWater[iLat] = body[iBody].daTempWater[iLat];
-          }
-        }
-
-        fvTempGradientSea(body, body[iBody].dSeasDeltax, iBody);
-        for (iLat = 0; iLat < body[iBody].iNumLats; iLat++) {
-          body[iBody].daDMidPt[iLat] =
-                0.5 * (body[iBody].daDiffusionSea[iLat + 1] +
-                       body[iBody].daDiffusionSea[iLat]);
-          body[iBody].daFlux[iLat] =
-                -2. * PI * (body[iBody].dRadius * body[iBody].dRadius) *
-                sqrt(1.0 - (sin(body[iBody].daLats[iLat]) *
-                            sin(body[iBody].daLats[iLat]))) *
-                body[iBody].daDMidPt[iLat] * body[iBody].daTGrad[iLat];
-          body[iBody].daFluxAvg[iLat] +=
-                body[iBody].daFlux[iLat] / body[iBody].iNStepInYear;
-          body[iBody]
-                .daFluxDaily[iLat][iNyear * body[iBody].iNStepInYear + iNstep] =
-                body[iBody].daFlux[iLat];
-
-          body[iBody].daDivFlux[iLat] = 0.0;
-          for (jLat = 0; jLat < body[iBody].iNumLats; jLat++) {
-            body[iBody].daDivFlux[iLat] += -body[iBody].daMDiffSea[iLat][jLat] *
-                                           body[iBody].daTempLW[jLat];
-          }
-          body[iBody].daDivFluxAvg[iLat] +=
-                body[iBody].daDivFlux[iLat] / body[iBody].iNStepInYear;
-          body[iBody].daDivFluxDaily[iLat][iNyear * body[iBody].iNStepInYear +
-                                           iNstep] =
-                body[iBody].daDivFlux[iLat];
-        }
-
-        body[iBody].dTGlobal +=
-              body[iBody].dTGlobalTmp / (body[iBody].iNStepInYear);
-        body[iBody].dFluxOutGlobal +=
-              body[iBody].dFluxOutGlobalTmp / (body[iBody].iNStepInYear);
-        body[iBody].dFluxInGlobal +=
-              body[iBody].dFluxInGlobalTmp / (body[iBody].iNStepInYear);
-
+        fvCalculateSeaIce(body, dStepsize, iBody, iNyear, iNday, iNstep);
       } else {
-        for (iLat = 0; iLat < 2 * body[iBody].iNumLats; iLat++) {
-          body[iBody].daTmpTempSea[iLat] = 0.0;
-          for (jLat = 0; jLat < 2 * body[iBody].iNumLats; jLat++) {
-            body[iBody].daTmpTempSea[iLat] +=
-                  body[iBody].daInvMSea[iLat][jLat] *
-                  body[iBody].daSourceLW[jLat];
-          }
-        }
-        for (iLat = 0; iLat < body[iBody].iNumLats; iLat++) {
-          // temp change this time step
-          body[iBody].daDeltaTempL[iLat] = body[iBody].daTmpTempSea[2 * iLat] -
-                                           body[iBody].daTempLand[iLat];
-          body[iBody].daDeltaTempW[iLat] =
-                body[iBody].daTmpTempSea[2 * iLat + 1] -
-                body[iBody].daTempWater[iLat];
-          // calculate temperature and fluxes by latitude and global average
-          body[iBody].daTempLand[iLat] = body[iBody].daTmpTempSea[2 * iLat];
-          body[iBody].daTempWater[iLat] =
-                body[iBody].daTmpTempSea[2 * iLat + 1];
-
-          fvSnowball(body, iBody);
-          // ice growth/ablation
-          if (body[iBody].bIceSheets) {
-            // calculate derivative of ice mass density and take an euler step
-            body[iBody].daIceBalance[iLat][iNstep] =
-                  fdIceMassBalance(body, iBody, iLat);
-            body[iBody].daIceMassTmp[iLat] +=
-                  dStepsize * body[iBody].daIceBalance[iLat][iNstep];
-            //             if (iLat==0) printf("iNstep = %d,
-            //             ice mass 0 =
-            //             %f\n",iNstep,body[iBody].daIceMassTmp[iLat]);
-            if (body[iBody].daIceBalance[iLat][iNstep] >= 0) {
-              body[iBody].daIceAccumTot[iLat] +=
-                    dStepsize * body[iBody].daIceBalance[iLat][iNstep] /
-                    body[iBody].iNumYears;
-            } else {
-              if (body[iBody].daIceMassTmp[iLat] > 0) {
-                if (body[iBody].daIceMassTmp[iLat] >=
-                    dStepsize * body[iBody].daIceBalance[iLat][iNstep]) {
-
-                  body[iBody].daIceAblateTot[iLat] +=
-                        dStepsize * body[iBody].daIceBalance[iLat][iNstep] /
-                        body[iBody].iNumYears;
-                } else {
-                  body[iBody].daIceAblateTot[iLat] +=
-                        body[iBody].daIceMassTmp[iLat];
-                }
-              }
-            }
-            if (body[iBody].daIceMassTmp[iLat] < 0.0) {
-              body[iBody].daIceMassTmp[iLat] = 0.0;
-            } // don't let ice mass become negative
-            if (body[iBody].daIceBalance[iLat][iNstep] < 0 &&
-                body[iBody].daIceMassTmp[iLat] != 0) {
-
-              if (body[iBody].daIceMassTmp[iLat] <=
-                  fabs(dStepsize * body[iBody].daIceBalance[iLat][iNstep])) {
-                // adjust temperature
-                body[iBody].daTempLand[iLat] +=
-                      -body[iBody].daIceMassTmp[iLat] * LFICE /
-                      body[iBody].dHeatCapLand;
-              } else {
-                // adjust temperature
-                body[iBody].daTempLand[iLat] +=
-                      dStepsize * body[iBody].daIceBalance[iLat][iNstep] *
-                      LFICE / body[iBody].dHeatCapLand;
-              }
-            } else if (body[iBody].daIceBalance[iLat][iNstep] > 0) {
-              // adjust temperature
-              body[iBody].daTempLand[iLat] +=
-                    dStepsize * body[iBody].daIceBalance[iLat][iNstep] * LFICE /
-                    body[iBody].dHeatCapLand;
-            }
-          }
-
-          body[iBody].daTempLW[iLat] =
-                body[iBody].daLandFrac[iLat] * body[iBody].daTempLand[iLat] +
-                body[iBody].daWaterFrac[iLat] * body[iBody].daTempWater[iLat];
-          body[iBody].dTGlobalTmp +=
-                body[iBody].daTempLW[iLat] / body[iBody].iNumLats;
-
-          if (body[iBody].bCalcAB) {
-            if (body[iBody].iOLRModel == WK97) {
-              body[iBody].daPlanckBSea[iLat] =
-                    fdOLRdTwk97(body, iBody, iLat, SEA);
-              body[iBody].daPlanckASea[iLat] =
-                    fdOLRwk97(body, iBody, iLat, SEA) -
-                    body[iBody].daPlanckBSea[iLat] *
-                          (body[iBody].daTempLW[iLat]);
-            } else if (body[iBody].iOLRModel == HM16) {
-              /* Calculate A and B from haqq-misra+ 2016 result */
-              body[iBody].daPlanckBSea[iLat] =
-                    fdOLRdThm16(body, iBody, iLat, SEA);
-              body[iBody].daPlanckASea[iLat] =
-                    fdOLRhm16(body, iBody, iLat, SEA) -
-                    body[iBody].daPlanckBSea[iLat] *
-                          (body[iBody].daTempLW[iLat]);
-            } else {
-              /* Calculate A and B from spiegel+ 2009 model */
-              body[iBody].daPlanckBSea[iLat] =
-                    fdOLRdTsms09(body, iBody, iLat, SEA);
-              body[iBody].daPlanckASea[iLat] =
-                    fdOLRsms09(body, iBody, iLat, SEA) -
-                    body[iBody].daPlanckBSea[iLat] *
-                          (body[iBody].daTempLW[iLat]);
-            }
-
-            if (body[iBody].bMEPDiff) {
-              if (iLat == 0) {
-                body[iBody].daDiffusionSea[iLat] =
-                      body[iBody].daPlanckBSea[iLat] / 4.0;
-              } else if (iLat == (body[iBody].iNumLats - 1)) {
-                body[iBody].daDiffusionSea[iLat] =
-                      (body[iBody].daPlanckBSea[iLat] +
-                       body[iBody].daPlanckBSea[iLat - 1]) /
-                      8.0;
-                body[iBody].daDiffusionSea[iLat + 1] =
-                      body[iBody].daPlanckBSea[iLat] / 4.0;
-              } else {
-                body[iBody].daDiffusionSea[iLat] =
-                      (body[iBody].daPlanckBSea[iLat] +
-                       body[iBody].daPlanckBSea[iLat - 1]) /
-                      8.0;
-              }
-            }
-          }
-
-          body[iBody].daFluxOutLand[iLat] =
-                body[iBody].daPlanckASea[iLat] +
-                body[iBody].daPlanckBSea[iLat] * body[iBody].daTempLand[iLat];
-          body[iBody].daFluxOutWater[iLat] =
-                body[iBody].daPlanckASea[iLat] +
-                body[iBody].daPlanckBSea[iLat] * body[iBody].daTempWater[iLat];
-          body[iBody].daFluxOut[iLat] =
-                body[iBody].daLandFrac[iLat] * body[iBody].daFluxOutLand[iLat] +
-                body[iBody].daWaterFrac[iLat] *
-                      body[iBody].daFluxOutWater[iLat];
-          body[iBody].dFluxOutGlobalTmp +=
-                body[iBody].daFluxOut[iLat] / body[iBody].iNumLats;
-
-          body[iBody].daFluxInLand[iLat] =
-                (1.0 - body[iBody].daAlbedoLand[iLat]) *
-                body[iBody].daInsol[iLat][iNday];
-          body[iBody].daFluxInWater[iLat] =
-                (1.0 - body[iBody].daAlbedoWater[iLat]) *
-                body[iBody].daInsol[iLat][iNday];
-          body[iBody].daFluxIn[iLat] =
-                body[iBody].daLandFrac[iLat] * body[iBody].daFluxInLand[iLat] +
-                body[iBody].daWaterFrac[iLat] * body[iBody].daFluxInWater[iLat];
-          body[iBody].dFluxInGlobalTmp +=
-                body[iBody].daFluxIn[iLat] / body[iBody].iNumLats;
-
-          // annual averages by latitude
-          body[iBody].daTempAvg[iLat] +=
-                body[iBody].daTempLW[iLat] / body[iBody].iNStepInYear;
-          body[iBody].daPlanckBAvg[iLat] +=
-                body[iBody].daPlanckBSea[iLat] / body[iBody].iNStepInYear;
-          body[iBody].daTempAvgL[iLat] +=
-                body[iBody].daTempLand[iLat] / body[iBody].iNStepInYear;
-          body[iBody].daTempAvgW[iLat] +=
-                body[iBody].daTempWater[iLat] / body[iBody].iNStepInYear;
-          body[iBody].daAlbedoAvgL[iLat] +=
-                body[iBody].daAlbedoLand[iLat] / body[iBody].iNStepInYear;
-          body[iBody].daAlbedoAvgW[iLat] +=
-                body[iBody].daAlbedoWater[iLat] / body[iBody].iNStepInYear;
-          body[iBody].daAlbedoAvg[iLat] +=
-                body[iBody].daAlbedoLW[iLat] / body[iBody].iNStepInYear;
-          body[iBody].daFluxInAvg[iLat] +=
-                body[iBody].daFluxIn[iLat] / body[iBody].iNStepInYear;
-          body[iBody].daFluxOutAvg[iLat] +=
-                body[iBody].daFluxOut[iLat] / body[iBody].iNStepInYear;
-
-          body[iBody]
-                .daTempDaily[iLat][iNyear * body[iBody].iNStepInYear + iNstep] =
-                body[iBody].daTempLW[iLat];
-          body[iBody].daFluxInDaily[iLat][iNyear * body[iBody].iNStepInYear +
-                                          iNstep] = body[iBody].daFluxIn[iLat];
-          body[iBody].daFluxOutDaily[iLat][iNyear * body[iBody].iNStepInYear +
-                                           iNstep] =
-                body[iBody].daFluxOut[iLat];
-          body[iBody].daPlanckBDaily[iLat][iNyear * body[iBody].iNStepInYear +
-                                           iNstep] =
-                body[iBody].daPlanckBSea[iLat];
-          if (body[iBody].daTempLW[iLat] < body[iBody].daTempMinLW[iLat]) {
-            body[iBody].daTempMinLW[iLat] = body[iBody].daTempLW[iLat];
-          }
-          if (body[iBody].daTempLW[iLat] > body[iBody].daTempMaxLW[iLat]) {
-            body[iBody].daTempMaxLW[iLat] = body[iBody].daTempLW[iLat];
-          }
-          if (body[iBody].daTempLand[iLat] > body[iBody].daTempMaxLand[iLat]) {
-            body[iBody].daTempMaxLand[iLat] = body[iBody].daTempLand[iLat];
-          }
-          if (body[iBody].daTempWater[iLat] >
-              body[iBody].daTempMaxWater[iLat]) {
-            body[iBody].daTempMaxWater[iLat] = body[iBody].daTempWater[iLat];
-          }
-        }
-
-        fvTempGradientSea(body, body[iBody].dSeasDeltax, iBody);
-        for (iLat = 0; iLat < body[iBody].iNumLats; iLat++) {
-          body[iBody].daDMidPt[iLat] =
-                0.5 * (body[iBody].daDiffusionSea[iLat + 1] +
-                       body[iBody].daDiffusionSea[iLat]);
-          body[iBody].daFlux[iLat] =
-                -2. * PI * (body[iBody].dRadius * body[iBody].dRadius) *
-                sqrt(1.0 - (sin(body[iBody].daLats[iLat]) *
-                            sin(body[iBody].daLats[iLat]))) *
-                body[iBody].daDMidPt[iLat] * body[iBody].daTGrad[iLat];
-          body[iBody].daFluxAvg[iLat] +=
-                body[iBody].daFlux[iLat] / body[iBody].iNStepInYear;
-          body[iBody]
-                .daFluxDaily[iLat][iNyear * body[iBody].iNStepInYear + iNstep] =
-                body[iBody].daFlux[iLat];
-
-          body[iBody].daDivFlux[iLat] = 0.0;
-          for (jLat = 0; jLat < body[iBody].iNumLats; jLat++) {
-            body[iBody].daDivFlux[iLat] += -body[iBody].daMDiffSea[iLat][jLat] *
-                                           body[iBody].daTempLW[jLat];
-          }
-          body[iBody].daDivFluxAvg[iLat] +=
-                body[iBody].daDivFlux[iLat] / body[iBody].iNStepInYear;
-          body[iBody].daDivFluxDaily[iLat][iNyear * body[iBody].iNStepInYear +
-                                           iNstep] =
-                body[iBody].daDivFlux[iLat];
-        }
-
-        body[iBody].dTGlobal +=
-              body[iBody].dTGlobalTmp / (body[iBody].iNStepInYear);
-        body[iBody].dFluxOutGlobal +=
-              body[iBody].dFluxOutGlobalTmp / (body[iBody].iNStepInYear);
-        body[iBody].dFluxInGlobal +=
-              body[iBody].dFluxInGlobalTmp / (body[iBody].iNStepInYear);
+        fvStaticSeaIce(body, dStepsize, iBody, iNyear, iNday, iNstep);
       }
 
+      fvTempGradientSea(body, body[iBody].dSeasDeltax, iBody);
+      fvFluxesByLatitude(body, iBody, iNyear, iNstep);
+      fvPoiseGlobalProps(body, iBody);
+
       EnergyResiduals(body, iBody, iNday);
       fvAlbedoSeasonal(body, iBody, iNday);
       if (body[iBody].bAccuracyMode) {
diff --git a/src/system.c b/src/system.c
index fc780d7b4..bc1431639 100644
--- a/src/system.c
+++ b/src/system.c
@@ -30,6 +30,8 @@ double fdSemiToMeanMotion(double dSemi, double dMass) {
   return pow(BIGG * dMass / (dSemi * dSemi * dSemi), 0.5);
 }
 
+
+
 /*
  * Angular Momentum
  */
diff --git a/src/vplanet.c b/src/vplanet.c
index 976d79204..99d977f3e 100644
--- a/src/vplanet.c
+++ b/src/vplanet.c
@@ -33,9 +33,9 @@ int main_impl(int argc, char *argv[]) {
   _MM_SET_EXCEPTION_MASK(_MM_GET_EXCEPTION_MASK() & ~_MM_MASK_INVALID);
   _MM_SET_EXCEPTION_MASK(_MM_GET_EXCEPTION_MASK() & ~_MM_MASK_OVERFLOW);
 #else
-  fprintf(
-        stderr,
-        "WARNING: Floating point trapping only enabled for x86 architectures.");
+  fprintf(stderr,
+          "WARNING: Floating point trapping only enabled for x86 "
+          "architectures.\n");
 #endif
 #endif
 
@@ -79,7 +79,7 @@ int main_impl(int argc, char *argv[]) {
   options = malloc(MODULEOPTEND * sizeof(OPTIONS));
   InitializeOptions(options, fnRead);
   output = malloc(MODULEOUTEND * sizeof(OUTPUT));
-  InitializeOutput(output, fnWrite);
+  InitializeOutput(&files, output, fnWrite);
 
   /* Set to IntegrationMethod to 0, so default can be
      assigned if necessary */


Test patch:
diff --git a/tests/CPL_POISE/gl514.in b/tests/CPL_POISE/gl514.in
new file mode 100644
index 000000000..f340e2f4b
--- /dev/null
+++ b/tests/CPL_POISE/gl514.in
@@ -0,0 +1,23 @@
+# Gl 514 parameters
+sName        gl514
+
+# Physical Parameters
+dMass         0.510
+#dSemi         0
+#dEcc          0
+dObliquity    0
+dRotPeriod    -30.55
+dRadGyra      0.5		# Radius of gyration (moment of inertia constant)
+dRadius       -54.64919165
+dLuminosity   1.65378e25
+
+# EQTIDE Parameters
+#dTidalTau    -1		# Tidal time lag, negative -> seconds
+dTidalQ	      1e6	# Tidal phase lag
+dK2           0.5	# Love number of degree 2
+sTideModel    p2		# Tidal model, p2=CPL, t8=CTL
+#sTideModel   t8		# Uncomment sTideModel, dTidalTau && dTidalQ to compare
+saTidePerts   gl514b		# Body name(s) of tidal perturberssStellarModel BARAFFE
+
+saModules     eqtide stellar
+saOutputOrder Time -RotPer
diff --git a/tests/CPL_POISE/gl514b.in b/tests/CPL_POISE/gl514b.in
new file mode 100644
index 000000000..c705debce
--- /dev/null
+++ b/tests/CPL_POISE/gl514b.in
@@ -0,0 +1,79 @@
+sName             gl514b
+saModules         poise eqtide
+dMass             1.56169977e-05
+dRadius           -2.10 # From Chen & Kippin 2017 forecaster calc
+dRotPeriod        -1.00000
+dObliquity        23.5
+#bForceObliq      1
+#dObliqAmp        180.
+#dObliqPer        1000
+#dOrbPeriod       0.384739726 # Either the semi or orbital period should be given, not both
+dSemi             0.422
+dEcc              0.45
+#bForceEcc        1
+#dEccPer          1000
+dPrecA            0.0
+#dDynEllip        0.0     # Shape of planet (0 = sphere)
+#dPrecRate        0.0
+dRefHeight        1000
+
+#_______eqtide parameters (have fun with these!)_________________________________________
+dRadGyra          0.5     # Radius of gyration (moment of inertia constant)
+bDiscreteRot      0       # Setting to 0 since ecc is greater than 0.324
+dTidalQ           12      # Tidal Q Earth value
+dTidalTau        -638     # Earth value
+dK2               0.3     # Love Number of degree 2
+saTidePerts       gl514   # Name of tidal perturber (only central body may have >1)
+dMaxLockDiff      0.1
+
+
+#_______addition disorb/distrot parameters (leave these alone for now)__________________
+dInc              0
+dLongA            0
+#bForcePrecRate   0
+#bOverrideMaxEcc  1
+#dHaltMaxEcc      0.4
+
+#_______poise parameters (have fun with these!)_________________________________________
+iLatCellNum       150
+sClimateModel     sea
+dTGlobalInit      14.85
+iNumYears         10
+iNStepInYear      60
+#dSurfAlbedo      0.4
+
+#__ice params_________
+bIceSheets        1
+dInitIceLat       90.
+dInitIceHeight    0.
+dIceDepRate       2.9e-5
+dIceAlbedo        0.6
+iIceDt            5
+iReRunSeas        500
+bSeaIceModel      0
+#iAlbedoType      fix
+bElevFB           0
+
+#__heat diffusion______
+#bMEPDiff         1
+#bHadley          1
+dDiffusion        0.44
+dNuLandWater      0.81
+
+#__outgoing flux_______
+dPlanckA          203.3
+dPlanckB          2.09
+bCalcAB           0
+#dpCO2            0.00028
+
+#__surface properties__
+dAlbedoLand       0.363
+dAlbedoWater      0.263
+dHeatCapLand      1.42e7
+dHeatCapWater     4.2e6
+dMixingDepth      70
+
+#________output options!_____________________________________________
+dSeasOutputTime   0
+
+
diff --git a/tests/CPL_POISE/test_CPL_POISE.py b/tests/CPL_POISE/test_CPL_POISE.py
new file mode 100644
index 000000000..187ed5d35
--- /dev/null
+++ b/tests/CPL_POISE/test_CPL_POISE.py
@@ -0,0 +1,436 @@
+from benchmark import Benchmark, benchmark 
+import astropy.units as u 
+import pytest 
+ 
+@benchmark( 
+   { 
+       "log.initial.system.Age": {"value": 0.0000000000000000, "unit": u.sec}, 
+       "log.initial.system.Time": {"value": 0.0000000000000000, "unit": u.sec}, 
+       "log.initial.system.TotAngMom": {"value": 8.9534036048945070e+41, "unit": (u.kg * u.m ** 2) / u.sec}, 
+       "log.initial.system.TotEnergy": {"value": -3.1347047276701209e+40, "unit": u.Joule}, 
+       "log.initial.system.PotEnergy": {"value": -3.1348025171406187e+40, "unit": u.Joule}, 
+       "log.initial.system.KinEnergy": {"value": 9.9742042956183368e+35, "unit": u.Joule}, 
+       "log.initial.system.DeltaTime": {"value": 0.0000000000000000, "unit": u.sec}, 
+       "log.initial.gl514.Mass": {"value": 1.0140921600000001e+30, "unit": u.kg}, 
+       "log.initial.gl514.Obliquity": {"value": 0.0000000000000000, "unit": u.rad}, 
+       "log.initial.gl514.PrecA": {"value": 0.0000000000000000, "unit": u.rad}, 
+       "log.initial.gl514.Xobl": {"value": 0.0000000000000000}, 
+       "log.initial.gl514.Yobl": {"value": 0.0000000000000000}, 
+       "log.initial.gl514.Zobl": {"value": 1.0000000000000000}, 
+       "log.initial.gl514.Radius": {"value": 1.3137125118741000e+09, "unit": u.m}, 
+       "log.initial.gl514.RadGyra": {"value": 0.4484976455000000}, 
+       "log.initial.gl514.RotAngMom": {"value": 8.3801799358956994e+41, "unit": (u.kg * u.m ** 2) / u.sec}, 
+       "log.initial.gl514.RotKinEnergy": {"value": 9.9742042956183368e+35, "unit": u.Joule}, 
+       "log.initial.gl514.RotVel": {"value": 3127.1970481244061375, "unit": u.m / u.sec}, 
+       "log.initial.gl514.BodyType": {"value": 0.0000000000000000}, 
+       "log.initial.gl514.RotRate": {"value": 2.3804272394903566e-06, "unit": 1 / u.sec}, 
+       "log.initial.gl514.RotPer": {"value": 30.5499999999999972, "unit": u.day}, 
+       "log.initial.gl514.Density": {"value": 106.7794814169220103, "unit": u.kg / u.m ** 3}, 
+       "log.initial.gl514.SurfEnFluxTotal": {"value": 5.6927005838186184e-10, "unit": u.kg / u.sec ** 3}, 
+       "log.initial.gl514.TidalQ": {"value": 1.0000000000000000e+06}, 
+       "log.initial.gl514.ImK2": {"value": -4.9999999999999998e-07}, 
+       "log.initial.gl514.K2": {"value": 0.5000000000000000}, 
+       "log.initial.gl514.K2Man": {"value": 0.0100000000000000}, 
+       "log.initial.gl514.Imk2Man": {"value": 0.0000000000000000}, 
+       "log.initial.gl514.TidalQMantle": {"value": 100.0000000000000000}, 
+       "log.initial.gl514.HEcc": {"value": 0.0000000000000000}, 
+       "log.initial.gl514.HZLimitDryRunaway": {"value": 1.1215301622470532e+11, "unit": u.m}, 
+       "log.initial.gl514.HZLimRecVenus": {"value": 9.9104353063518631e+10, "unit": u.m}, 
+       "log.initial.gl514.HZLimRunaway": {"value": 1.3136272247075568e+11, "unit": u.m}, 
+       "log.initial.gl514.HZLimMoistGreenhouse": {"value": 1.3113105818748537e+11, "unit": u.m}, 
+       "log.initial.gl514.HZLimMaxGreenhouse": {"value": 2.4324580025875378e+11, "unit": u.m}, 
+       "log.initial.gl514.HZLimEarlyMars": {"value": 2.6530630489685046e+11, "unit": u.m}, 
+       "log.initial.gl514.Instellation": {"value": -1.0000000000000000, "unit": u.kg / u.sec ** 3}, 
+       "log.initial.gl514.KEcc": {"value": 0.0000000000000000}, 
+       "log.initial.gl514.Eccentricity": {"value": -1.0000000000000000}, 
+       "log.initial.gl514.OrbEnergy": {"value": 0.0000000000000000, "unit": u.Joule}, 
+       "log.initial.gl514.MeanMotion": {"value": -1.0000000000000000, "unit": 1 / u.sec}, 
+       "log.initial.gl514.OrbPeriod": {"value": -1.0000000000000000, "unit": u.sec}, 
+       "log.initial.gl514.SemiMajorAxis": {"value": -1.0000000000000000, "unit": u.m}, 
+       "log.initial.gl514.CriticalSemiMajorAxis": {"value": -1.0000000000000000, "unit": u.m}, 
+       "log.initial.gl514.COPP": {"value": 0.0000000000000000}, 
+       "log.initial.gl514.OrbAngMom": {"value": 0.0000000000000000, "unit": (u.kg * u.m ** 2) / u.sec}, 
+       "log.initial.gl514.LongP": {"value": 0.0000000000000000, "unit": u.rad}, 
+       "log.initial.gl514.LXUVTot": {"value": 2.6238552211501671e+23, "unit": u.kg / u.sec ** 3}, 
+       "log.initial.gl514.TotOrbEnergy": {"value": -1.6646340895690071e+34, "unit": u.Joule}, 
+       "log.initial.gl514.OrbPotEnergy": {"value": -1.0000000000000000, "unit": u.Joule}, 
+       "log.initial.gl514.LostEnergy": {"value": 5.5626846462680035e-309, "unit": u.Joule}, 
+       "log.initial.gl514.LostAngMom": {"value": 5.5626846462680035e-309, "unit": (u.kg * u.m ** 2) / u.sec}, 
+       "log.initial.gl514.LockTime": {"value": -1.0000000000000000, "unit": u.sec}, 
+       "log.initial.gl514.BodyDsemiDtEqtide": {"value": -1.0000000000000000}, 
+       "log.initial.gl514.BodyDeccDt": {"value": -1.0000000000000000}, 
+       "log.initial.gl514.DOblDtEqtide": {"value": 0.0000000000000000, "unit": u.rad / u.sec}, 
+       "log.initial.gl514.DRotPerDtEqtide": {"value": 2.3670446824643342e-20}, 
+       "log.initial.gl514.DRotRateDtEqtide": {"value": -2.1346978387088952e-32, "unit": 1 / u.sec ** 2}, 
+       "log.initial.gl514.EqRotRateDiscrete": {"value": 7.7800437205935984e-07, "unit": 1 / u.sec}, 
+       "log.initial.gl514.EqRotPerDiscrete": {"value": 8.0760282754557505e+06, "unit": u.sec}, 
+       "log.initial.gl514.EqRotRateCont": {"value": 1.5164601885390357e-06, "unit": 1 / u.sec}, 
+       "log.initial.gl514.EqRotPerCont": {"value": 4.1433236128887986e+06, "unit": u.sec}, 
+       "log.initial.gl514.EqRotPer": {"value": 4.1433236128887986e+06, "unit": u.sec}, 
+       "log.initial.gl514.EqTidePower": {"value": 0.0000000000000000, "unit": 1 / u.sec}, 
+       "log.initial.gl514.GammaRot": {"value": -1.0000000000000000, "unit": u.sec}, 
+       "log.initial.gl514.GammaOrb": {"value": -1.0000000000000000, "unit": u.sec}, 
+       "log.initial.gl514.OceanK2": {"value": 0.0100000000000000}, 
+       "log.initial.gl514.EnvTidalQ": {"value": -1.0000000000000000}, 
+       "log.initial.gl514.OceanTidalQ": {"value": -1.0000000000000000}, 
+       "log.initial.gl514.TideLock": {"value": 0.0000000000000000}, 
+       "log.initial.gl514.RotTimeEqtide": {"value": 1.1151120295929485e+26, "unit": u.sec}, 
+       "log.initial.gl514.EnvK2": {"value": 0.0100000000000000}, 
+       "log.initial.gl514.OblTimeEqtide": {"value": -1.0000000000000000}, 
+       "log.initial.gl514.PowerEqtide": {"value": 1.2346074076706444e+10, "unit": u.W}, 
+       "log.initial.gl514.SurfEnFluxEqtide": {"value": 5.6927005838186184e-10, "unit": u.kg / u.sec ** 3}, 
+       "log.initial.gl514.Luminosity": {"value": 2.6238552211501671e+26, "unit": u.W}, 
+       "log.initial.gl514.LXUVStellar": {"value": 2.6238552211501671e+23, "unit": u.W}, 
+       "log.initial.gl514.Temperature": {"value": 3818.8020013755422042, "unit": u.K}, 
+       "log.initial.gl514.LXUVFrac": {"value": 0.0010000000000000}, 
+       "log.initial.gl514.RossbyNumber": {"value": 0.6873965284599931}, 
+       "log.initial.gl514.DRotPerDtStellar": {"value": 2.3504479341501690e-10}, 
+       "log.initial.gl514b.Mass": {"value": 3.1053088098643201e+25, "unit": u.kg}, 
+       "log.initial.gl514b.Obliquity": {"value": 0.4101523742069445, "unit": u.rad}, 
+       "log.initial.gl514b.PrecA": {"value": 0.0000000000000000, "unit": u.rad}, 
+       "log.initial.gl514b.Xobl": {"value": 0.3987490689144955}, 
+       "log.initial.gl514b.Yobl": {"value": 0.0000000000000000}, 
+       "log.initial.gl514b.Zobl": {"value": 0.9170600743897985}, 
+       "log.initial.gl514b.Radius": {"value": 1.3394010000000000e+07, "unit": u.m}, 
+       "log.initial.gl514b.RadGyra": {"value": 0.5000000000000000}, 
+       "log.initial.gl514b.RotAngMom": {"value": 1.0128197643484414e+35, "unit": (u.kg * u.m ** 2) / u.sec}, 
+       "log.initial.gl514b.RotKinEnergy": {"value": 3.6827165869069548e+30, "unit": u.Joule}, 
+       "log.initial.gl514b.RotVel": {"value": 974.0398939376904082, "unit": u.m / u.sec}, 
+       "log.initial.gl514b.BodyType": {"value": 0.0000000000000000}, 
+       "log.initial.gl514b.RotRate": {"value": 7.2722052166430395e-05, "unit": 1 / u.sec}, 
+       "log.initial.gl514b.RotPer": {"value": 8.6400000000000000e+04, "unit": u.sec}, 
+       "log.initial.gl514b.Density": {"value": 3085.2071528551509800, "unit": u.kg / u.m ** 3}, 
+       "log.initial.gl514b.SurfEnFluxTotal": {"value": 1.4096431035620620, "unit": u.kg / u.sec ** 3}, 
+       "log.initial.gl514b.TidalQ": {"value": 12.0000000000000000}, 
+       "log.initial.gl514b.ImK2": {"value": -0.0250000000000000}, 
+       "log.initial.gl514b.K2": {"value": 0.3000000000000000}, 
+       "log.initial.gl514b.K2Man": {"value": 0.0100000000000000}, 
+       "log.initial.gl514b.Imk2Man": {"value": 0.0000000000000000}, 
+       "log.initial.gl514b.TidalQMantle": {"value": 100.0000000000000000}, 
+       "log.initial.gl514b.HEcc": {"value": 0.0000000000000000}, 
+       "log.initial.gl514b.HZLimitDryRunaway": {"value": 1.0509802225988144e+11, "unit": u.m}, 
+       "log.initial.gl514b.HZLimRecVenus": {"value": 9.9104353063518631e+10, "unit": u.m}, 
+       "log.initial.gl514b.HZLimRunaway": {"value": 1.3136272247075568e+11, "unit": u.m}, 
+       "log.initial.gl514b.HZLimMoistGreenhouse": {"value": 1.3113105818748537e+11, "unit": u.m}, 
+       "log.initial.gl514b.HZLimMaxGreenhouse": {"value": 2.4324580025875378e+11, "unit": u.m}, 
+       "log.initial.gl514b.HZLimEarlyMars": {"value": 2.6530630489685046e+11, "unit": u.m}, 
+       "log.initial.gl514b.Instellation": {"value": 5866.6331509136580280, "unit": u.kg / u.sec ** 3}, 
+       "log.initial.gl514b.KEcc": {"value": 0.4500000000000000}, 
+       "log.initial.gl514b.Eccentricity": {"value": 0.4500000000000000}, 
+       "log.initial.gl514b.OrbEnergy": {"value": -1.6646340895690071e+34, "unit": u.Joule}, 
+       "log.initial.gl514b.MeanMotion": {"value": 5.1866958137290656e-07, "unit": 1 / u.sec}, 
+       "log.initial.gl514b.OrbPeriod": {"value": 1.2114042413183626e+07, "unit": u.sec}, 
+       "log.initial.gl514b.SemiMajorAxis": {"value": 6.3130301435400002e+10, "unit": u.m}, 
+       "log.initial.gl514b.CriticalSemiMajorAxis": {"value": -1.0000000000000000, "unit": u.m}, 
+       "log.initial.gl514b.COPP": {"value": 0.0000000000000000}, 
+       "log.initial.gl514b.OrbAngMom": {"value": 5.7322265617904285e+40, "unit": (u.kg * u.m ** 2) / u.sec}, 
+       "log.initial.gl514b.LongP": {"value": 0.0000000000000000, "unit": u.rad}, 
+       "log.initial.gl514b.TotOrbEnergy": {"value": -1.6646340895690071e+34, "unit": u.Joule}, 
+       "log.initial.gl514b.OrbPotEnergy": {"value": -3.3292681791380142e+34, "unit": u.Joule}, 
+       "log.initial.gl514b.LostEnergy": {"value": 5.5626846462680035e-309, "unit": u.Joule}, 
+       "log.initial.gl514b.TidalRadius": {"value": 1.3394010000000000e+07, "unit": u.m}, 
+       "log.initial.gl514b.DsemiDtEqtide": {"value": 1.2351242905116021e-10, "unit": u.m / u.sec}, 
+       "log.initial.gl514b.DeccDtEqtide": {"value": 5.8376560316556904e-22, "unit": 1 / u.sec}, 
+       "log.initial.gl514b.DMeanMotionDtEqtide": {"value": -1.5221408360263325e-27, "unit": 1 / u.sec ** 2}, 
+       "log.initial.gl514b.DOrbPerDtEqtide": {"value": 3.5551108660842267e-14}, 
+       "log.initial.gl514b.EccTimeEqtide": {"value": 7.7085733993198282e+20, "unit": u.sec}, 
+       "log.initial.gl514b.SemiTimeEqtide": {"value": 5.1112509016603289e+20, "unit": u.sec}, 
+       "log.initial.gl514b.DHEccDtEqtide": {"value": 0.0000000000000000, "unit": 1 / u.sec}, 
+       "log.initial.gl514b.DKEccDtEqtide": {"value": 5.8376560316556904e-22, "unit": 1 / u.sec}, 
+       "log.initial.gl514b.DXoblDtEqtide": {"value": 3.1641063776408136e-17, "unit": 1 / u.sec}, 
+       "log.initial.gl514b.DYoblDtEqtide": {"value": 0.0000000000000000, "unit": 1 / u.sec}, 
+       "log.initial.gl514b.DZoblDtEqtide": {"value": -1.3757926086469334e-17, "unit": 1 / u.sec}, 
+       "log.initial.gl514b.LockTime": {"value": -1.0000000000000000, "unit": u.sec}, 
+       "log.initial.gl514b.BodyDsemiDtEqtide": {"value": -1.0000000000000000}, 
+       "log.initial.gl514b.BodyDeccDt": {"value": -1.0000000000000000}, 
+       "log.initial.gl514b.DOblDtEqtide": {"value": 3.4502716517739303e-17, "unit": u.rad / u.sec}, 
+       "log.initial.gl514b.DRotPerDtEqtide": {"value": 3.7660240995856781e-11}, 
+       "log.initial.gl514b.DRotRateDtEqtide": {"value": -3.1698264008113854e-20, "unit": 1 / u.sec ** 2}, 
+       "log.initial.gl514b.EqRotRateDiscrete": {"value": 7.7800437205935984e-07, "unit": 1 / u.sec}, 
+       "log.initial.gl514b.EqRotPerDiscrete": {"value": 8.0760282754557505e+06, "unit": u.sec}, 
+       "log.initial.gl514b.EqRotRateCont": {"value": 1.5164601885390357e-06, "unit": 1 / u.sec}, 
+       "log.initial.gl514b.EqRotPerCont": {"value": 4.1433236128887986e+06, "unit": u.sec}, 
+       "log.initial.gl514b.EqRotPer": {"value": 4.1433236128887986e+06, "unit": u.sec}, 
+       "log.initial.gl514b.EqTidePower": {"value": 0.0000000000000000, "unit": 1 / u.sec}, 
+       "log.initial.gl514b.GammaRot": {"value": -1.0000000000000000, "unit": u.sec}, 
+       "log.initial.gl514b.GammaOrb": {"value": -1.0000000000000000, "unit": u.sec}, 
+       "log.initial.gl514b.OceanK2": {"value": 0.0100000000000000}, 
+       "log.initial.gl514b.EnvTidalQ": {"value": -1.0000000000000000}, 
+       "log.initial.gl514b.OceanTidalQ": {"value": -1.0000000000000000}, 
+       "log.initial.gl514b.TideLock": {"value": 0.0000000000000000}, 
+       "log.initial.gl514b.RotTimeEqtide": {"value": 2.2941966836990070e+15, "unit": u.sec}, 
+       "log.initial.gl514b.EnvK2": {"value": 0.0100000000000000}, 
+       "log.initial.gl514b.OblTimeEqtide": {"value": -1.0000000000000000}, 
+       "log.initial.gl514b.PowerEqtide": {"value": 3.1779003342801765e+15, "unit": u.W}, 
+       "log.initial.gl514b.SurfEnFluxEqtide": {"value": 1.4096431035620620, "unit": u.kg / u.sec ** 3}, 
+       "log.initial.gl514b.TGlobal": {"value": 676.2679353212223532, "unit": u.sec}, 
+       "log.initial.gl514b.AlbedoGlobal": {"value": 0.2996778493802214}, 
+       "log.initial.gl514b.FluxInGlobal": {"value": 1046.1316028565233864, "unit": u.kg / u.sec ** 3}, 
+       "log.initial.gl514b.FluxOutGlobal": {"value": 1046.1299848213550376, "unit": u.kg / u.sec ** 3}, 
+       "log.initial.gl514b.TotIceMass": {"value": 0.0000000000000000, "unit": u.kg}, 
+       "log.initial.gl514b.TotIceFlow": {"value": 0.0000000000000000, "unit": u.kg}, 
+       "log.initial.gl514b.TotIceBalance": {"value": 0.0000000000000000, "unit": u.kg}, 
+       "log.initial.gl514b.SkipSeas": {"value": 0.0000000000000000}, 
+       "log.initial.gl514b.AreaIceCov": {"value": 0.0000000000000000}, 
+       "log.initial.gl514b.Latitude": {"value": -1.4552620265106593, "unit": u.rad}, 
+       "log.initial.gl514b.TempLat": {"value": 557.2247797713517912, "unit": u.sec}, 
+       "log.initial.gl514b.AlbedoLat": {"value": 0.3675649878041329}, 
+       "log.initial.gl514b.AnnInsol": {"value": 755.3407450595735781, "unit": u.kg / u.sec ** 3}, 
+       "log.initial.gl514b.FluxMerid": {"value": -2.0352457824211096e+16, "unit": u.Joule}, 
+       "log.initial.gl514b.FluxIn": {"value": 486.1623336259239636, "unit": u.kg / u.sec ** 3}, 
+       "log.initial.gl514b.FluxOut": {"value": 797.3297897221252697, "unit": u.kg / u.sec ** 3}, 
+       "log.initial.gl514b.DivFlux": {"value": -311.1690738881308107, "unit": u.kg / u.sec ** 3}, 
+       "log.initial.gl514b.IceMass": {"value": 0.0000000000000000}, 
+       "log.initial.gl514b.IceHeight": {"value": 0.0000000000000000, "unit": u.m}, 
+       "log.initial.gl514b.DIceMassDt": {"value": 0.0000000000000000, "unit": u.m}, 
+       "log.initial.gl514b.IceFlow": {"value": 0.0000000000000000}, 
+       "log.initial.gl514b.EnergyResL": {"value": -4.7634785005357116e-11, "unit": u.kg / u.sec ** 3}, 
+       "log.initial.gl514b.EnergyResW": {"value": 4.0836312109604478e-10, "unit": u.kg / u.sec ** 3}, 
+       "log.initial.gl514b.BedrockH": {"value": 0.0000000000000000, "unit": u.m}, 
+       "log.initial.gl514b.TempLandLat": {"value": 545.8662422200939091, "unit": u.sec}, 
+       "log.initial.gl514b.TempWaterLat": {"value": 563.0761476007876354, "unit": u.sec}, 
+       "log.initial.gl514b.AlbedoLandLat": {"value": 0.4335649878041328}, 
+       "log.initial.gl514b.AlbedoWaterLat": {"value": 0.3335649878041328}, 
+       "log.initial.gl514b.TempMinLat": {"value": 520.6583069986143073, "unit": u.sec}, 
+       "log.initial.gl514b.TempMaxLat": {"value": 598.1727520370602633, "unit": u.sec}, 
+       "log.initial.gl514b.Snowball": {"value": 0.0000000000000000}, 
+       "log.initial.gl514b.PlanckBAvg": {"value": 2.0899999999999990}, 
+       "log.initial.gl514b.IceAccum": {"value": 0.0000000000000000}, 
+       "log.initial.gl514b.IceAblate": {"value": 0.0000000000000000}, 
+       "log.initial.gl514b.TempMaxLand": {"value": 657.9976852023412448, "unit": u.sec}, 
+       "log.initial.gl514b.TempMaxWater": {"value": 568.6086486017817379, "unit": u.sec}, 
+       "log.initial.gl514b.PeakInsol": {"value": 4243.8792145061070187, "unit": u.kg / u.sec ** 3}, 
+       "log.initial.gl514b.IceCapNorthLand": {"value": 0.0000000000000000}, 
+       "log.initial.gl514b.IceCapNorthSea": {"value": 0.0000000000000000}, 
+       "log.initial.gl514b.IceCapSouthLand": {"value": 0.0000000000000000}, 
+       "log.initial.gl514b.IceCapSouthSea": {"value": 0.0000000000000000}, 
+       "log.initial.gl514b.IceBeltLand": {"value": 0.0000000000000000}, 
+       "log.initial.gl514b.IceBeltSea": {"value": 0.0000000000000000}, 
+       "log.initial.gl514b.SnowballLand": {"value": 0.0000000000000000}, 
+       "log.initial.gl514b.SnowballSea": {"value": 0.0000000000000000}, 
+       "log.initial.gl514b.IceFree": {"value": 1.0000000000000000}, 
+       "log.initial.gl514b.IceCapNorthLatLand": {"value": 100.0000000000000000, "unit": u.rad}, 
+       "log.initial.gl514b.IceCapNorthLatSea": {"value": 100.0000000000000000, "unit": u.rad}, 
+       "log.initial.gl514b.IceCapSouthLatLand": {"value": 100.0000000000000000, "unit": u.rad}, 
+       "log.initial.gl514b.IceCapSouthLatSea": {"value": 100.0000000000000000, "unit": u.rad}, 
+       "log.initial.gl514b.IceBeltNorthLatLand": {"value": 100.0000000000000000, "unit": u.rad}, 
+       "log.initial.gl514b.IceBeltNorthLatSea": {"value": 100.0000000000000000, "unit": u.rad}, 
+       "log.initial.gl514b.IceBeltSouthLatLand": {"value": 100.0000000000000000, "unit": u.rad}, 
+       "log.initial.gl514b.IceBeltSouthLatSea": {"value": 100.0000000000000000, "unit": u.rad}, 
+       "log.final.system.Age": {"value": 3.1557600000000000e+07, "unit": u.sec, "rtol": 1e-4}, 
+       "log.final.system.Time": {"value": 3.1557600000000000e+07, "unit": u.sec, "rtol": 1e-4}, 
+       "log.final.system.TotAngMom": {"value": 8.9534036048945039e+41, "unit": (u.kg * u.m ** 2) / u.sec, "rtol": 1e-4}, 
+       "log.final.system.TotEnergy": {"value": -3.1347047276701209e+40, "unit": u.Joule, "rtol": 1e-4}, 
+       "log.final.system.PotEnergy": {"value": -3.1348025171406187e+40, "unit": u.Joule, "rtol": 1e-4}, 
+       "log.final.system.KinEnergy": {"value": 9.9742042395602989e+35, "unit": u.Joule, "rtol": 1e-4}, 
+       "log.final.system.DeltaTime": {"value": 3.1557600000000000e+07, "unit": u.sec, "rtol": 1e-4}, 
+       "log.final.gl514.Mass": {"value": 1.0140921600000001e+30, "unit": u.kg, "rtol": 1e-4}, 
+       "log.final.gl514.Obliquity": {"value": 0.0000000000000000, "unit": u.rad, "rtol": 1e-4}, 
+       "log.final.gl514.PrecA": {"value": 0.0000000000000000, "unit": u.rad, "rtol": 1e-4}, 
+       "log.final.gl514.Xobl": {"value": 0.0000000000000000, "rtol": 1e-4}, 
+       "log.final.gl514.Yobl": {"value": 0.0000000000000000, "rtol": 1e-4}, 
+       "log.final.gl514.Zobl": {"value": 1.0000000000000000, "rtol": 1e-4}, 
+       "log.final.gl514.Radius": {"value": 1.3137125118741000e+09, "unit": u.m, "rtol": 1e-4}, 
+       "log.final.gl514.RadGyra": {"value": 0.4484976455000000, "rtol": 1e-4}, 
+       "log.final.gl514.RotAngMom": {"value": 8.3801799123461292e+41, "unit": (u.kg * u.m ** 2) / u.sec, "rtol": 1e-4}, 
+       "log.final.gl514.RotKinEnergy": {"value": 9.9742042395602989e+35, "unit": u.Joule, "rtol": 1e-4}, 
+       "log.final.gl514.RotVel": {"value": 3127.1970393365104428, "unit": u.m / u.sec, "rtol": 1e-4}, 
+       "log.final.gl514.BodyType": {"value": 0.0000000000000000, "rtol": 1e-4}, 
+       "log.final.gl514.RotRate": {"value": 2.3804272328009969e-06, "unit": 1 / u.sec, "rtol": 1e-4}, 
+       "log.final.gl514.RotPer": {"value": 30.5500000858501117, "unit": u.day, "rtol": 1e-4}, 
+       "log.final.gl514.Density": {"value": 106.7794814169220103, "unit": u.kg / u.m ** 3, "rtol": 1e-4}, 
+       "log.final.gl514.SurfEnFluxTotal": {"value": 5.6927005606372488e-10, "unit": u.kg / u.sec ** 3, "rtol": 1e-4}, 
+       "log.final.gl514.TidalQ": {"value": 1.0000000000000000e+06, "rtol": 1e-4}, 
+       "log.final.gl514.ImK2": {"value": -4.9999999999999998e-07, "rtol": 1e-4}, 
+       "log.final.gl514.K2": {"value": 0.5000000000000000, "rtol": 1e-4}, 
+       "log.final.gl514.K2Man": {"value": 0.0100000000000000, "rtol": 1e-4}, 
+       "log.final.gl514.Imk2Man": {"value": 0.0000000000000000, "rtol": 1e-4}, 
+       "log.final.gl514.TidalQMantle": {"value": 100.0000000000000000, "rtol": 1e-4}, 
+       "log.final.gl514.HEcc": {"value": 0.0000000000000000, "rtol": 1e-4}, 
+       "log.final.gl514.HZLimitDryRunaway": {"value": 1.1215301622470532e+11, "unit": u.m, "rtol": 1e-4}, 
+       "log.final.gl514.HZLimRecVenus": {"value": 9.9104353063518631e+10, "unit": u.m, "rtol": 1e-4}, 
+       "log.final.gl514.HZLimRunaway": {"value": 1.3136272247075568e+11, "unit": u.m, "rtol": 1e-4}, 
+       "log.final.gl514.HZLimMoistGreenhouse": {"value": 1.3113105818748537e+11, "unit": u.m, "rtol": 1e-4}, 
+       "log.final.gl514.HZLimMaxGreenhouse": {"value": 2.4324580025875378e+11, "unit": u.m, "rtol": 1e-4}, 
+       "log.final.gl514.HZLimEarlyMars": {"value": 2.6530630489685046e+11, "unit": u.m, "rtol": 1e-4}, 
+       "log.final.gl514.Instellation": {"value": -1.0000000000000000, "unit": u.kg / u.sec ** 3, "rtol": 1e-4}, 
+       "log.final.gl514.KEcc": {"value": 0.0000000000000000, "rtol": 1e-4}, 
+       "log.final.gl514.Eccentricity": {"value": -1.0000000000000000, "rtol": 1e-4}, 
+       "log.final.gl514.OrbEnergy": {"value": 0.0000000000000000, "unit": u.Joule, "rtol": 1e-4}, 
+       "log.final.gl514.MeanMotion": {"value": -1.0000000000000000, "unit": 1 / u.sec, "rtol": 1e-4}, 
+       "log.final.gl514.OrbPeriod": {"value": -1.0000000000000000, "unit": u.sec, "rtol": 1e-4}, 
+       "log.final.gl514.SemiMajorAxis": {"value": -1.0000000000000000, "unit": u.m, "rtol": 1e-4}, 
+       "log.final.gl514.CriticalSemiMajorAxis": {"value": -1.0000000000000000, "unit": u.m, "rtol": 1e-4}, 
+       "log.final.gl514.COPP": {"value": 0.0000000000000000, "rtol": 1e-4}, 
+       "log.final.gl514.OrbAngMom": {"value": 0.0000000000000000, "unit": (u.kg * u.m ** 2) / u.sec, "rtol": 1e-4}, 
+       "log.final.gl514.LongP": {"value": 0.0000000000000000, "unit": u.rad, "rtol": 1e-4}, 
+       "log.final.gl514.LXUVTot": {"value": 2.6238552211501671e+23, "unit": u.kg / u.sec ** 3, "rtol": 1e-4}, 
+       "log.final.gl514.TotOrbEnergy": {"value": -1.6646340895689042e+34, "unit": u.Joule, "rtol": 1e-4}, 
+       "log.final.gl514.OrbPotEnergy": {"value": -1.0000000000000000, "unit": u.Joule, "rtol": 1e-4}, 
+       "log.final.gl514.LostEnergy": {"value": 5.6058038409805720e+27, "unit": u.Joule, "rtol": 1e-4}, 
+       "log.final.gl514.LostAngMom": {"value": 2.3549570242977578e+33, "unit": (u.kg * u.m ** 2) / u.sec, "rtol": 1e-4}, 
+       "log.final.gl514.LockTime": {"value": -1.0000000000000000, "unit": u.sec, "rtol": 1e-4}, 
+       "log.final.gl514.BodyDsemiDtEqtide": {"value": -1.0000000000000000, "rtol": 1e-4}, 
+       "log.final.gl514.BodyDeccDt": {"value": -1.0000000000000000, "rtol": 1e-4}, 
+       "log.final.gl514.DOblDtEqtide": {"value": 0.0000000000000000, "unit": u.rad / u.sec, "rtol": 1e-4}, 
+       "log.final.gl514.DRotPerDtEqtide": {"value": 2.3670446957670796e-20, "rtol": 1e-4}, 
+       "log.final.gl514.DRotRateDtEqtide": {"value": -2.1346978387082101e-32, "unit": 1 / u.sec ** 2, "rtol": 1e-4}, 
+       "log.final.gl514.EqRotRateDiscrete": {"value": 7.7800437205928774e-07, "unit": 1 / u.sec, "rtol": 1e-4}, 
+       "log.final.gl514.EqRotPerDiscrete": {"value": 8.0760282754564993e+06, "unit": u.sec, "rtol": 1e-4}, 
+       "log.final.gl514.EqRotRateCont": {"value": 1.5164601885389766e-06, "unit": 1 / u.sec, "rtol": 1e-4}, 
+       "log.final.gl514.EqRotPerCont": {"value": 4.1433236128889602e+06, "unit": u.sec, "rtol": 1e-4}, 
+       "log.final.gl514.EqRotPer": {"value": 4.1433236128889602e+06, "unit": u.sec, "rtol": 1e-4}, 
+       "log.final.gl514.EqTidePower": {"value": 0.0000000000000000, "unit": 1 / u.sec, "rtol": 1e-4}, 
+       "log.final.gl514.GammaRot": {"value": -1.0000000000000000, "unit": u.sec, "rtol": 1e-4}, 
+       "log.final.gl514.GammaOrb": {"value": -1.0000000000000000, "unit": u.sec, "rtol": 1e-4}, 
+       "log.final.gl514.OceanK2": {"value": 0.0100000000000000, "rtol": 1e-4}, 
+       "log.final.gl514.EnvTidalQ": {"value": -1.0000000000000000, "rtol": 1e-4}, 
+       "log.final.gl514.OceanTidalQ": {"value": -1.0000000000000000, "rtol": 1e-4}, 
+       "log.final.gl514.TideLock": {"value": 0.0000000000000000, "rtol": 1e-4}, 
+       "log.final.gl514.RotTimeEqtide": {"value": 1.1151120264596733e+26, "unit": u.sec, "rtol": 1e-4}, 
+       "log.final.gl514.EnvK2": {"value": 0.0100000000000000, "rtol": 1e-4}, 
+       "log.final.gl514.OblTimeEqtide": {"value": -1.0000000000000000, "rtol": 1e-4}, 
+       "log.final.gl514.PowerEqtide": {"value": 1.2346074026431728e+10, "unit": u.W, "rtol": 1e-4}, 
+       "log.final.gl514.SurfEnFluxEqtide": {"value": 5.6927005606372488e-10, "unit": u.kg / u.sec ** 3, "rtol": 1e-4}, 
+       "log.final.gl514.Luminosity": {"value": 2.6238552211501671e+26, "unit": u.W, "rtol": 1e-4}, 
+       "log.final.gl514.LXUVStellar": {"value": 2.6238552211501671e+23, "unit": u.W, "rtol": 1e-4}, 
+       "log.final.gl514.Temperature": {"value": 3818.8020013755422042, "unit": u.K, "rtol": 1e-4}, 
+       "log.final.gl514.LXUVFrac": {"value": 0.0010000000000000, "rtol": 1e-4}, 
+       "log.final.gl514.RossbyNumber": {"value": 0.6873965303916811, "rtol": 1e-4}, 
+       "log.final.gl514.DRotPerDtStellar": {"value": 2.3504479143348279e-10, "rtol": 1e-4}, 
+       "log.final.gl514b.Mass": {"value": 3.1053088098643201e+25, "unit": u.kg, "rtol": 1e-4}, 
+       "log.final.gl514b.Obliquity": {"value": 0.4101523752957674, "unit": u.rad, "rtol": 1e-4}, 
+       "log.final.gl514b.PrecA": {"value": 0.0000000000000000, "unit": u.rad, "rtol": 1e-4}, 
+       "log.final.gl514b.Xobl": {"value": 0.3987490699130116, "rtol": 1e-4}, 
+       "log.final.gl514b.Yobl": {"value": 0.0000000000000000, "rtol": 1e-4}, 
+       "log.final.gl514b.Zobl": {"value": 0.9170600739556314, "rtol": 1e-4}, 
+       "log.final.gl514b.Radius": {"value": 1.3394010000000000e+07, "unit": u.m, "rtol": 1e-4}, 
+       "log.final.gl514b.RadGyra": {"value": 0.5000000000000000, "rtol": 1e-4}, 
+       "log.final.gl514b.RotAngMom": {"value": 1.0128197504166973e+35, "unit": (u.kg * u.m ** 2) / u.sec, "rtol": 1e-4}, 
+       "log.final.gl514b.RotKinEnergy": {"value": 3.6827164855924537e+30, "unit": u.Joule, "rtol": 1e-4}, 
+       "log.final.gl514b.RotVel": {"value": 974.0398805393790553, "unit": u.m / u.sec, "rtol": 1e-4}, 
+       "log.final.gl514b.BodyType": {"value": 0.0000000000000000, "rtol": 1e-4}, 
+       "log.final.gl514b.RotRate": {"value": 7.2722051166109258e-05, "unit": 1 / u.sec, "rtol": 1e-4}, 
+       "log.final.gl514b.RotPer": {"value": 8.6400001188466835e+04, "unit": u.sec, "rtol": 1e-4}, 
+       "log.final.gl514b.Density": {"value": 3085.2071528551509800, "unit": u.kg / u.m ** 3, "rtol": 1e-4}, 
+       "log.final.gl514b.SurfEnFluxTotal": {"value": 1.4096430839727192, "unit": u.kg / u.sec ** 3, "rtol": 1e-4}, 
+       "log.final.gl514b.TidalQ": {"value": 12.0000000000000000, "rtol": 1e-4}, 
+       "log.final.gl514b.ImK2": {"value": -0.0250000000000000, "rtol": 1e-4}, 
+       "log.final.gl514b.K2": {"value": 0.3000000000000000, "rtol": 1e-4}, 
+       "log.final.gl514b.K2Man": {"value": 0.0100000000000000, "rtol": 1e-4}, 
+       "log.final.gl514b.Imk2Man": {"value": 0.0000000000000000, "rtol": 1e-4}, 
+       "log.final.gl514b.TidalQMantle": {"value": 100.0000000000000000, "rtol": 1e-4}, 
+       "log.final.gl514b.HEcc": {"value": 0.0000000000000000, "rtol": 1e-4}, 
+       "log.final.gl514b.HZLimitDryRunaway": {"value": 1.0509802225887555e+11, "unit": u.m, "rtol": 1e-4}, 
+       "log.final.gl514b.HZLimRecVenus": {"value": 9.9104353063518631e+10, "unit": u.m, "rtol": 1e-4}, 
+       "log.final.gl514b.HZLimRunaway": {"value": 1.3136272247075568e+11, "unit": u.m, "rtol": 1e-4}, 
+       "log.final.gl514b.HZLimMoistGreenhouse": {"value": 1.3113105818748537e+11, "unit": u.m, "rtol": 1e-4}, 
+       "log.final.gl514b.HZLimMaxGreenhouse": {"value": 2.4324580025875378e+11, "unit": u.m, "rtol": 1e-4}, 
+       "log.final.gl514b.HZLimEarlyMars": {"value": 2.6530630489685046e+11, "unit": u.m, "rtol": 1e-4}, 
+       "log.final.gl514b.Instellation": {"value": 5866.6331509129950064, "unit": u.kg / u.sec ** 3, "rtol": 1e-4}, 
+       "log.final.gl514b.KEcc": {"value": 0.4500000000000184, "rtol": 1e-4}, 
+       "log.final.gl514b.Eccentricity": {"value": 0.4500000000000184, "rtol": 1e-4}, 
+       "log.final.gl514b.OrbEnergy": {"value": -1.6646340895689042e+34, "unit": u.Joule, "rtol": 1e-4}, 
+       "log.final.gl514b.MeanMotion": {"value": 5.1866958137285849e-07, "unit": 1 / u.sec, "rtol": 1e-4}, 
+       "log.final.gl514b.OrbPeriod": {"value": 1.2114042413184749e+07, "unit": u.sec, "rtol": 1e-4}, 
+       "log.final.gl514b.SemiMajorAxis": {"value": 6.3130301435403900e+10, "unit": u.m, "rtol": 1e-4}, 
+       "log.final.gl514b.CriticalSemiMajorAxis": {"value": -1.0000000000000000, "unit": u.m, "rtol": 1e-4}, 
+       "log.final.gl514b.COPP": {"value": 0.0000000000000000, "rtol": 1e-4}, 
+       "log.final.gl514b.OrbAngMom": {"value": 5.7322265617905455e+40, "unit": (u.kg * u.m ** 2) / u.sec, "rtol": 1e-4}, 
+       "log.final.gl514b.LongP": {"value": 0.0000000000000000, "unit": u.rad, "rtol": 1e-4}, 
+       "log.final.gl514b.TotOrbEnergy": {"value": -1.6646340895689042e+34, "unit": u.Joule, "rtol": 1e-4}, 
+       "log.final.gl514b.OrbPotEnergy": {"value": -3.3292681791378085e+34, "unit": u.Joule, "rtol": 1e-4}, 
+       "log.final.gl514b.LostEnergy": {"value": 1.0028690689225319e+23, "unit": u.Joule, "rtol": 1e-4}, 
+       "log.final.gl514b.TidalRadius": {"value": 1.3394010000000000e+07, "unit": u.m, "rtol": 1e-4}, 
+       "log.final.gl514b.DsemiDtEqtide": {"value": 1.2351242905112557e-10, "unit": u.m / u.sec, "rtol": 1e-4}, 
+       "log.final.gl514b.DeccDtEqtide": {"value": 5.8376560316535877e-22, "unit": 1 / u.sec, "rtol": 1e-4}, 
+       "log.final.gl514b.DMeanMotionDtEqtide": {"value": -1.5221408360256703e-27, "unit": 1 / u.sec ** 2, "rtol": 1e-4}, 
+       "log.final.gl514b.DOrbPerDtEqtide": {"value": 3.5551108660833406e-14, "rtol": 1e-4}, 
+       "log.final.gl514b.EccTimeEqtide": {"value": 7.7085733993229202e+20, "unit": u.sec, "rtol": 1e-4}, 
+       "log.final.gl514b.SemiTimeEqtide": {"value": 5.1112509016620781e+20, "unit": u.sec, "rtol": 1e-4}, 
+       "log.final.gl514b.DHEccDtEqtide": {"value": 0.0000000000000000, "unit": 1 / u.sec, "rtol": 1e-4}, 
+       "log.final.gl514b.DKEccDtEqtide": {"value": 5.8376560316535877e-22, "unit": 1 / u.sec, "rtol": 1e-4}, 
+       "log.final.gl514b.DXoblDtEqtide": {"value": 3.1641064275885809e-17, "unit": 1 / u.sec, "rtol": 1e-4}, 
+       "log.final.gl514b.DYoblDtEqtide": {"value": 0.0000000000000000, "unit": 1 / u.sec, "rtol": 1e-4}, 
+       "log.final.gl514b.DZoblDtEqtide": {"value": -1.3757926344613390e-17, "unit": 1 / u.sec, "rtol": 1e-4}, 
+       "log.final.gl514b.LockTime": {"value": -1.0000000000000000, "unit": u.sec, "rtol": 1e-4}, 
+       "log.final.gl514b.BodyDsemiDtEqtide": {"value": -1.0000000000000000, "rtol": 1e-4}, 
+       "log.final.gl514b.BodyDeccDt": {"value": -1.0000000000000000, "rtol": 1e-4}, 
+       "log.final.gl514b.DOblDtEqtide": {"value": 3.4502717078725040e-17, "unit": u.rad / u.sec, "rtol": 1e-4}, 
+       "log.final.gl514b.DRotPerDtEqtide": {"value": 3.7660242031908298e-11, "rtol": 1e-4}, 
+       "log.final.gl514b.DRotRateDtEqtide": {"value": -3.1698264008103682e-20, "unit": 1 / u.sec ** 2, "rtol": 1e-4}, 
+       "log.final.gl514b.EqRotRateDiscrete": {"value": 7.7800437205928774e-07, "unit": 1 / u.sec, "rtol": 1e-4}, 
+       "log.final.gl514b.EqRotPerDiscrete": {"value": 8.0760282754564993e+06, "unit": u.sec, "rtol": 1e-4}, 
+       "log.final.gl514b.EqRotRateCont": {"value": 1.5164601885389766e-06, "unit": 1 / u.sec, "rtol": 1e-4}, 
+       "log.final.gl514b.EqRotPerCont": {"value": 4.1433236128889602e+06, "unit": u.sec, "rtol": 1e-4}, 
+       "log.final.gl514b.EqRotPer": {"value": 4.1433236128889602e+06, "unit": u.sec, "rtol": 1e-4}, 
+       "log.final.gl514b.EqTidePower": {"value": 0.0000000000000000, "unit": 1 / u.sec, "rtol": 1e-4}, 
+       "log.final.gl514b.GammaRot": {"value": -1.0000000000000000, "unit": u.sec, "rtol": 1e-4}, 
+       "log.final.gl514b.GammaOrb": {"value": -1.0000000000000000, "unit": u.sec, "rtol": 1e-4}, 
+       "log.final.gl514b.OceanK2": {"value": 0.0100000000000000, "rtol": 1e-4}, 
+       "log.final.gl514b.EnvTidalQ": {"value": -1.0000000000000000, "rtol": 1e-4}, 
+       "log.final.gl514b.OceanTidalQ": {"value": -1.0000000000000000, "rtol": 1e-4}, 
+       "log.final.gl514b.TideLock": {"value": 0.0000000000000000, "rtol": 1e-4}, 
+       "log.final.gl514b.RotTimeEqtide": {"value": 2.2941966521421430e+15, "unit": u.sec, "rtol": 1e-4}, 
+       "log.final.gl514b.EnvK2": {"value": 0.0100000000000000, "rtol": 1e-4}, 
+       "log.final.gl514b.OblTimeEqtide": {"value": -1.0000000000000000, "rtol": 1e-4}, 
+       "log.final.gl514b.PowerEqtide": {"value": 3.1779002901179495e+15, "unit": u.W, "rtol": 1e-4}, 
+       "log.final.gl514b.SurfEnFluxEqtide": {"value": 1.4096430839727192, "unit": u.kg / u.sec ** 3, "rtol": 1e-4}, 
+       "log.final.gl514b.TGlobal": {"value": 676.2684488390327715, "unit": u.sec, "rtol": 1e-4}, 
+       "log.final.gl514b.AlbedoGlobal": {"value": 0.2996778493936415, "rtol": 1e-4}, 
+       "log.final.gl514b.FluxInGlobal": {"value": 1046.1316028588294103, "unit": u.kg / u.sec ** 3, "rtol": 1e-4}, 
+       "log.final.gl514b.FluxOutGlobal": {"value": 1046.1310580735782878, "unit": u.kg / u.sec ** 3, "rtol": 1e-4}, 
+       "log.final.gl514b.TotIceMass": {"value": 0.0000000000000000, "unit": u.kg, "rtol": 1e-4}, 
+       "log.final.gl514b.TotIceFlow": {"value": 0.0000000000000000, "unit": u.kg, "rtol": 1e-4}, 
+       "log.final.gl514b.TotIceBalance": {"value": 0.0000000000000000, "unit": u.kg, "rtol": 1e-4}, 
+       "log.final.gl514b.SkipSeas": {"value": 0.0000000000000000, "rtol": 1e-4}, 
+       "log.final.gl514b.AreaIceCov": {"value": 0.0000000000000000, "rtol": 1e-4}, 
+       "log.final.gl514b.Latitude": {"value": -1.4552620265106593, "unit": u.rad, "rtol": 1e-4}, 
+       "log.final.gl514b.TempLat": {"value": 557.2252934226596608, "unit": u.sec, "rtol": 1e-4}, 
+       "log.final.gl514b.AlbedoLat": {"value": 0.3675649877648788, "rtol": 1e-4}, 
+       "log.final.gl514b.AnnInsol": {"value": 755.3407468728647700, "unit": u.kg / u.sec ** 3, "rtol": 1e-4}, 
+       "log.final.gl514b.FluxMerid": {"value": -2.0352457779997064e+16, "unit": u.Joule, "rtol": 1e-4}, 
+       "log.final.gl514b.FluxIn": {"value": 486.1623348490650187, "unit": u.kg / u.sec ** 3, "rtol": 1e-4}, 
+       "log.final.gl514b.FluxOut": {"value": 797.3308632533586433, "unit": u.kg / u.sec ** 3, "rtol": 1e-4}, 
+       "log.final.gl514b.DivFlux": {"value": -311.1690732121418819, "unit": u.kg / u.sec ** 3, "rtol": 1e-4}, 
+       "log.final.gl514b.IceMass": {"value": 0.0000000000000000, "rtol": 1e-4}, 
+       "log.final.gl514b.IceHeight": {"value": 0.0000000000000000, "unit": u.m, "rtol": 1e-4}, 
+       "log.final.gl514b.DIceMassDt": {"value": 0.0000000000000000, "unit": u.m, "rtol": 1e-4}, 
+       "log.final.gl514b.IceFlow": {"value": 0.0000000000000000, "rtol": 1e-4}, 
+       "log.final.gl514b.EnergyResL": {"value": 7.9239725891966373e-11, "unit": u.kg / u.sec ** 3, "rtol": 1e-4}, 
+       "log.final.gl514b.EnergyResW": {"value": 6.8553163146134466e-11, "unit": u.kg / u.sec ** 3, "rtol": 1e-4}, 
+       "log.final.gl514b.BedrockH": {"value": 0.0000000000000000, "unit": u.m, "rtol": 1e-4}, 
+       "log.final.gl514b.TempLandLat": {"value": 545.8665750487834885, "unit": u.sec, "rtol": 1e-4}, 
+       "log.final.gl514b.TempWaterLat": {"value": 563.0767544031414218, "unit": u.sec, "rtol": 1e-4}, 
+       "log.final.gl514b.AlbedoLandLat": {"value": 0.4335649877648788, "rtol": 1e-4}, 
+       "log.final.gl514b.AlbedoWaterLat": {"value": 0.3335649877648789, "rtol": 1e-4}, 
+       "log.final.gl514b.TempMinLat": {"value": 520.6587997956662548, "unit": u.sec, "rtol": 1e-4}, 
+       "log.final.gl514b.TempMaxLat": {"value": 598.1732888468328611, "unit": u.sec, "rtol": 1e-4}, 
+       "log.final.gl514b.Snowball": {"value": 0.0000000000000000, "rtol": 1e-4}, 
+       "log.final.gl514b.PlanckBAvg": {"value": 2.0899999999999990, "rtol": 1e-4}, 
+       "log.final.gl514b.IceAccum": {"value": 0.0000000000000000, "rtol": 1e-4}, 
+       "log.final.gl514b.IceAblate": {"value": 0.0000000000000000, "rtol": 1e-4}, 
+       "log.final.gl514b.TempMaxLand": {"value": 657.9980332140321480, "unit": u.sec, "rtol": 1e-4}, 
+       "log.final.gl514b.TempMaxWater": {"value": 568.6092758126177387, "unit": u.sec, "rtol": 1e-4}, 
+       "log.final.gl514b.PeakInsol": {"value": 4243.8792251330132785, "unit": u.kg / u.sec ** 3, "rtol": 1e-4}, 
+       "log.final.gl514b.IceCapNorthLand": {"value": 0.0000000000000000, "rtol": 1e-4}, 
+       "log.final.gl514b.IceCapNorthSea": {"value": 0.0000000000000000, "rtol": 1e-4}, 
+       "log.final.gl514b.IceCapSouthLand": {"value": 0.0000000000000000, "rtol": 1e-4}, 
+       "log.final.gl514b.IceCapSouthSea": {"value": 0.0000000000000000, "rtol": 1e-4}, 
+       "log.final.gl514b.IceBeltLand": {"value": 0.0000000000000000, "rtol": 1e-4}, 
+       "log.final.gl514b.IceBeltSea": {"value": 0.0000000000000000, "rtol": 1e-4}, 
+       "log.final.gl514b.SnowballLand": {"value": 0.0000000000000000, "rtol": 1e-4}, 
+       "log.final.gl514b.SnowballSea": {"value": 0.0000000000000000, "rtol": 1e-4}, 
+       "log.final.gl514b.IceFree": {"value": 1.0000000000000000, "rtol": 1e-4}, 
+       "log.final.gl514b.IceCapNorthLatLand": {"value": 100.0000000000000000, "unit": u.rad, "rtol": 1e-4}, 
+       "log.final.gl514b.IceCapNorthLatSea": {"value": 100.0000000000000000, "unit": u.rad, "rtol": 1e-4}, 
+       "log.final.gl514b.IceCapSouthLatLand": {"value": 100.0000000000000000, "unit": u.rad, "rtol": 1e-4}, 
+       "log.final.gl514b.IceCapSouthLatSea": {"value": 100.0000000000000000, "unit": u.rad, "rtol": 1e-4}, 
+       "log.final.gl514b.IceBeltNorthLatLand": {"value": 100.0000000000000000, "unit": u.rad, "rtol": 1e-4}, 
+       "log.final.gl514b.IceBeltNorthLatSea": {"value": 100.0000000000000000, "unit": u.rad, "rtol": 1e-4}, 
+       "log.final.gl514b.IceBeltSouthLatLand": {"value": 100.0000000000000000, "unit": u.rad, "rtol": 1e-4}, 
+       "log.final.gl514b.IceBeltSouthLatSea": {"value": 100.0000000000000000, "unit": u.rad, "rtol": 1e-4}, 
+   } 
+)
+class TestCPL_POISE(Benchmark): 
+   pass 
diff --git a/tests/CPL_POISE/vpl.in b/tests/CPL_POISE/vpl.in
new file mode 100644
index 000000000..d152c18b6
--- /dev/null
+++ b/tests/CPL_POISE/vpl.in
@@ -0,0 +1,16 @@
+sSystemName   gl514
+iVerbose      5
+bOverwrite    1
+sUnitMass     solar
+sUnitLength   au
+sUnitTime     y
+sUnitAngle    d
+sUnitTemp     K
+bDoLog        1
+iDigits       16
+saBodyFiles   gl514.in gl514b.in
+bDoForward    1
+bVarDt        1
+dEta          0.1
+dStopTime     1
+dOutputTime   1
diff --git a/tests/CTL_POISE/gl514.in b/tests/CTL_POISE/gl514.in
new file mode 100644
index 000000000..49f4acf9c
--- /dev/null
+++ b/tests/CTL_POISE/gl514.in
@@ -0,0 +1,23 @@
+# Gl 514 parameters
+sName        gl514
+
+# Physical Parameters
+dMass         0.510
+#dSemi         0
+#dEcc          0
+dObliquity    0
+dRotPeriod    -30.55
+dRadGyra      0.5		# Radius of gyration (moment of inertia constant)
+dRadius       -54.64919165
+dLuminosity   1.65378e25
+
+# EQTIDE Parameters
+dTidalTau    -0.1		# Tidal time lag, negative -> seconds
+dTidalQ	      1e6	# Tidal phase lag
+dK2           0.5	# Love number of degree 2
+#sTideModel    p2		# Tidal model, p2=CPL, t8=CTL
+sTideModel   t8		# Uncomment sTideModel, dTidalTau && dTidalQ to compare
+saTidePerts   gl514b		# Body name(s) of tidal perturberssStellarModel BARAFFE
+
+saModules     eqtide stellar
+saOutputOrder Time -RotPer
diff --git a/tests/CTL_POISE/gl514b.in b/tests/CTL_POISE/gl514b.in
new file mode 100644
index 000000000..08184f0cd
--- /dev/null
+++ b/tests/CTL_POISE/gl514b.in
@@ -0,0 +1,79 @@
+sName             gl514b
+saModules         poise eqtide
+dMass             1.56169977e-05
+dRadius           -2.10 # From Chen & Kippin 2017 forecaster calc
+dRotPeriod        -1.00000
+dObliquity        23.5
+#bForceObliq      1
+#dObliqAmp        180.
+#dObliqPer        1000
+#dOrbPeriod       0.384739726 # Either the semi or orbital period should be given, not both
+dSemi             0.422
+dEcc              0.45
+#bForceEcc        1
+#dEccPer          1000
+dPrecA            0.0
+#dDynEllip        0.0     # Shape of planet (0 = sphere)
+#dPrecRate        0.0
+dRefHeight        1000
+
+#_______eqtide parameters (have fun with these!)_________________________________________
+dRadGyra          0.5     # Radius of gyration (moment of inertia constant)
+#bDiscreteRot      0       # Setting to 0 since ecc is greater than 0.324
+dTidalQ           12      # Tidal Q Earth value
+dTidalTau        -638     # Earth value
+dK2               0.3     # Love Number of degree 2
+saTidePerts       gl514   # Name of tidal perturber (only central body may have >1)
+dMaxLockDiff      0.1
+
+
+#_______addition disorb/distrot parameters (leave these alone for now)__________________
+dInc              0
+dLongA            0
+#bForcePrecRate   0
+#bOverrideMaxEcc  1
+#dHaltMaxEcc      0.4
+
+#_______poise parameters (have fun with these!)_________________________________________
+iLatCellNum       150
+sClimateModel     sea
+dTGlobalInit      14.85
+iNumYears         10
+iNStepInYear      60
+#dSurfAlbedo      0.4
+
+#__ice params_________
+bIceSheets        1
+dInitIceLat       90.
+dInitIceHeight    0.
+dIceDepRate       2.9e-5
+dIceAlbedo        0.6
+iIceDt            5
+iReRunSeas        500
+bSeaIceModel      0
+#iAlbedoType      fix
+bElevFB           0
+
+#__heat diffusion______
+#bMEPDiff         1
+#bHadley          1
+dDiffusion        0.44
+dNuLandWater      0.81
+
+#__outgoing flux_______
+dPlanckA          203.3
+dPlanckB          2.09
+bCalcAB           0
+#dpCO2            0.00028
+
+#__surface properties__
+dAlbedoLand       0.363
+dAlbedoWater      0.263
+dHeatCapLand      1.42e7
+dHeatCapWater     4.2e6
+dMixingDepth      70
+
+#________output options!_____________________________________________
+dSeasOutputTime   0
+
+
diff --git a/tests/CTL_POISE/test_CTL_POISE.py b/tests/CTL_POISE/test_CTL_POISE.py
new file mode 100644
index 000000000..c52d19b99
--- /dev/null
+++ b/tests/CTL_POISE/test_CTL_POISE.py
@@ -0,0 +1,412 @@
+from benchmark import Benchmark, benchmark 
+import astropy.units as u 
+import pytest 
+ 
+@benchmark( 
+   { 
+       "log.initial.system.Age": {"value": 0.0000000000000000, "unit": u.sec}, 
+       "log.initial.system.Time": {"value": 0.0000000000000000, "unit": u.sec}, 
+       "log.initial.system.TotAngMom": {"value": 8.9534036048945070e+41, "unit": (u.kg * u.m ** 2) / u.sec}, 
+       "log.initial.system.TotEnergy": {"value": -3.1347047276701209e+40, "unit": u.Joule}, 
+       "log.initial.system.PotEnergy": {"value": -3.1348025171406187e+40, "unit": u.Joule}, 
+       "log.initial.system.KinEnergy": {"value": 9.9742042956183368e+35, "unit": u.Joule}, 
+       "log.initial.system.DeltaTime": {"value": 0.0000000000000000, "unit": u.sec}, 
+       "log.initial.gl514.Mass": {"value": 1.0140921600000001e+30, "unit": u.kg}, 
+       "log.initial.gl514.Obliquity": {"value": 0.0000000000000000, "unit": u.rad}, 
+       "log.initial.gl514.PrecA": {"value": 0.0000000000000000, "unit": u.rad}, 
+       "log.initial.gl514.Xobl": {"value": 0.0000000000000000}, 
+       "log.initial.gl514.Yobl": {"value": 0.0000000000000000}, 
+       "log.initial.gl514.Zobl": {"value": 1.0000000000000000}, 
+       "log.initial.gl514.Radius": {"value": 1.3137125118741000e+09, "unit": u.m}, 
+       "log.initial.gl514.RadGyra": {"value": 0.4484976455000000}, 
+       "log.initial.gl514.RotAngMom": {"value": 8.3801799358956994e+41, "unit": (u.kg * u.m ** 2) / u.sec}, 
+       "log.initial.gl514.RotKinEnergy": {"value": 9.9742042956183368e+35, "unit": u.Joule}, 
+       "log.initial.gl514.RotVel": {"value": 3127.1970481244061375, "unit": u.m / u.sec}, 
+       "log.initial.gl514.BodyType": {"value": 0.0000000000000000}, 
+       "log.initial.gl514.RotRate": {"value": 2.3804272394903566e-06, "unit": 1 / u.sec}, 
+       "log.initial.gl514.RotPer": {"value": 30.5499999999999972, "unit": u.day}, 
+       "log.initial.gl514.Density": {"value": 106.7794814169220103, "unit": u.kg / u.m ** 3}, 
+       "log.initial.gl514.SurfEnFluxTotal": {"value": 1.9756861074797136e-10, "unit": u.kg / u.sec ** 3}, 
+       "log.initial.gl514.ImK2": {"value": -4.9999999999999998e-07}, 
+       "log.initial.gl514.K2": {"value": 0.5000000000000000}, 
+       "log.initial.gl514.K2Man": {"value": 0.0100000000000000}, 
+       "log.initial.gl514.Imk2Man": {"value": 0.0000000000000000}, 
+       "log.initial.gl514.TidalQMantle": {"value": 100.0000000000000000}, 
+       "log.initial.gl514.HEcc": {"value": 0.0000000000000000}, 
+       "log.initial.gl514.HZLimitDryRunaway": {"value": 1.1215301622470532e+11, "unit": u.m}, 
+       "log.initial.gl514.HZLimRecVenus": {"value": 9.9104353063518631e+10, "unit": u.m}, 
+       "log.initial.gl514.HZLimRunaway": {"value": 1.3136272247075568e+11, "unit": u.m}, 
+       "log.initial.gl514.HZLimMoistGreenhouse": {"value": 1.3113105818748537e+11, "unit": u.m}, 
+       "log.initial.gl514.HZLimMaxGreenhouse": {"value": 2.4324580025875378e+11, "unit": u.m}, 
+       "log.initial.gl514.HZLimEarlyMars": {"value": 2.6530630489685046e+11, "unit": u.m}, 
+       "log.initial.gl514.Instellation": {"value": -1.0000000000000000, "unit": u.kg / u.sec ** 3}, 
+       "log.initial.gl514.KEcc": {"value": 0.0000000000000000}, 
+       "log.initial.gl514.Eccentricity": {"value": -1.0000000000000000}, 
+       "log.initial.gl514.OrbEnergy": {"value": 0.0000000000000000, "unit": u.Joule}, 
+       "log.initial.gl514.MeanMotion": {"value": -1.0000000000000000, "unit": 1 / u.sec}, 
+       "log.initial.gl514.OrbPeriod": {"value": -1.0000000000000000, "unit": u.sec}, 
+       "log.initial.gl514.SemiMajorAxis": {"value": -1.0000000000000000, "unit": u.m}, 
+       "log.initial.gl514.CriticalSemiMajorAxis": {"value": -1.0000000000000000, "unit": u.m}, 
+       "log.initial.gl514.COPP": {"value": 0.0000000000000000}, 
+       "log.initial.gl514.OrbAngMom": {"value": 0.0000000000000000, "unit": (u.kg * u.m ** 2) / u.sec}, 
+       "log.initial.gl514.LongP": {"value": 0.0000000000000000, "unit": u.rad}, 
+       "log.initial.gl514.LXUVTot": {"value": 2.6238552211501671e+23, "unit": u.kg / u.sec ** 3}, 
+       "log.initial.gl514.TotOrbEnergy": {"value": -1.6646340895690071e+34, "unit": u.Joule}, 
+       "log.initial.gl514.OrbPotEnergy": {"value": -1.0000000000000000, "unit": u.Joule}, 
+       "log.initial.gl514.LostEnergy": {"value": 5.5626846462680035e-309, "unit": u.Joule}, 
+       "log.initial.gl514.LostAngMom": {"value": 5.5626846462680035e-309, "unit": (u.kg * u.m ** 2) / u.sec}, 
+       "log.initial.gl514.LockTime": {"value": -1.0000000000000000, "unit": u.sec}, 
+       "log.initial.gl514.BodyDsemiDtEqtide": {"value": -1.0000000000000000}, 
+       "log.initial.gl514.BodyDeccDt": {"value": -1.0000000000000000}, 
+       "log.initial.gl514.DOblDtEqtide": {"value": 0.0000000000000000, "unit": u.rad / u.sec}, 
+       "log.initial.gl514.DRotPerDtEqtide": {"value": 9.6722127519728598e-21}, 
+       "log.initial.gl514.DRotRateDtEqtide": {"value": -8.7227975923433734e-33, "unit": 1 / u.sec ** 2}, 
+       "log.initial.gl514.EqRotPer": {"value": 5.0877011498424802e+06, "unit": u.sec}, 
+       "log.initial.gl514.EqTidePower": {"value": -1.0000000000000000, "unit": 1 / u.sec}, 
+       "log.initial.gl514.OceanK2": {"value": 0.0100000000000000}, 
+       "log.initial.gl514.EnvTidalQ": {"value": -1.0000000000000000}, 
+       "log.initial.gl514.OceanTidalQ": {"value": -1.0000000000000000}, 
+       "log.initial.gl514.TidalTau": {"value": 0.1000000000000000, "unit": u.sec}, 
+       "log.initial.gl514.TideLock": {"value": 0.0000000000000000}, 
+       "log.initial.gl514.RotTimeEqtide": {"value": 2.7289722297119776e+26, "unit": u.sec}, 
+       "log.initial.gl514.EnvK2": {"value": 0.0100000000000000}, 
+       "log.initial.gl514.OblTimeEqtide": {"value": -1.0000000000000000}, 
+       "log.initial.gl514.PowerEqtide": {"value": 4.2847795481459198e+09, "unit": u.W}, 
+       "log.initial.gl514.SurfEnFluxEqtide": {"value": 1.9756861074797136e-10, "unit": u.kg / u.sec ** 3}, 
+       "log.initial.gl514.Luminosity": {"value": 2.6238552211501671e+26, "unit": u.W}, 
+       "log.initial.gl514.LXUVStellar": {"value": 2.6238552211501671e+23, "unit": u.W}, 
+       "log.initial.gl514.Temperature": {"value": 3818.8020013755422042, "unit": u.K}, 
+       "log.initial.gl514.LXUVFrac": {"value": 0.0010000000000000}, 
+       "log.initial.gl514.RossbyNumber": {"value": 0.6873965284599931}, 
+       "log.initial.gl514.DRotPerDtStellar": {"value": 2.3504479341501690e-10}, 
+       "log.initial.gl514b.Mass": {"value": 3.1053088098643201e+25, "unit": u.kg}, 
+       "log.initial.gl514b.Obliquity": {"value": 0.4101523742069445, "unit": u.rad}, 
+       "log.initial.gl514b.PrecA": {"value": 0.0000000000000000, "unit": u.rad}, 
+       "log.initial.gl514b.Xobl": {"value": 0.3987490689144955}, 
+       "log.initial.gl514b.Yobl": {"value": 0.0000000000000000}, 
+       "log.initial.gl514b.Zobl": {"value": 0.9170600743897985}, 
+       "log.initial.gl514b.Radius": {"value": 1.3394010000000000e+07, "unit": u.m}, 
+       "log.initial.gl514b.RadGyra": {"value": 0.5000000000000000}, 
+       "log.initial.gl514b.RotAngMom": {"value": 1.0128197643484414e+35, "unit": (u.kg * u.m ** 2) / u.sec}, 
+       "log.initial.gl514b.RotKinEnergy": {"value": 3.6827165869069548e+30, "unit": u.Joule}, 
+       "log.initial.gl514b.RotVel": {"value": 974.0398939376904082, "unit": u.m / u.sec}, 
+       "log.initial.gl514b.BodyType": {"value": 0.0000000000000000}, 
+       "log.initial.gl514b.RotRate": {"value": 7.2722052166430395e-05, "unit": 1 / u.sec}, 
+       "log.initial.gl514b.RotPer": {"value": 8.6400000000000000e+04, "unit": u.sec}, 
+       "log.initial.gl514b.Density": {"value": 3085.2071528551509800, "unit": u.kg / u.m ** 3}, 
+       "log.initial.gl514b.SurfEnFluxTotal": {"value": 2.5164492816849116, "unit": u.kg / u.sec ** 3}, 
+       "log.initial.gl514b.ImK2": {"value": -0.0250000000000000}, 
+       "log.initial.gl514b.K2": {"value": 0.3000000000000000}, 
+       "log.initial.gl514b.K2Man": {"value": 0.0100000000000000}, 
+       "log.initial.gl514b.Imk2Man": {"value": 0.0000000000000000}, 
+       "log.initial.gl514b.TidalQMantle": {"value": 100.0000000000000000}, 
+       "log.initial.gl514b.HEcc": {"value": 0.0000000000000000}, 
+       "log.initial.gl514b.HZLimitDryRunaway": {"value": 1.0509802225988144e+11, "unit": u.m}, 
+       "log.initial.gl514b.HZLimRecVenus": {"value": 9.9104353063518631e+10, "unit": u.m}, 
+       "log.initial.gl514b.HZLimRunaway": {"value": 1.3136272247075568e+11, "unit": u.m}, 
+       "log.initial.gl514b.HZLimMoistGreenhouse": {"value": 1.3113105818748537e+11, "unit": u.m}, 
+       "log.initial.gl514b.HZLimMaxGreenhouse": {"value": 2.4324580025875378e+11, "unit": u.m}, 
+       "log.initial.gl514b.HZLimEarlyMars": {"value": 2.6530630489685046e+11, "unit": u.m}, 
+       "log.initial.gl514b.Instellation": {"value": 5866.6331509136580280, "unit": u.kg / u.sec ** 3}, 
+       "log.initial.gl514b.KEcc": {"value": 0.4500000000000000}, 
+       "log.initial.gl514b.Eccentricity": {"value": 0.4500000000000000}, 
+       "log.initial.gl514b.OrbEnergy": {"value": -1.6646340895690071e+34, "unit": u.Joule}, 
+       "log.initial.gl514b.MeanMotion": {"value": 5.1866958137290656e-07, "unit": 1 / u.sec}, 
+       "log.initial.gl514b.OrbPeriod": {"value": 1.2114042413183626e+07, "unit": u.sec}, 
+       "log.initial.gl514b.SemiMajorAxis": {"value": 6.3130301435400002e+10, "unit": u.m}, 
+       "log.initial.gl514b.CriticalSemiMajorAxis": {"value": -1.0000000000000000, "unit": u.m}, 
+       "log.initial.gl514b.COPP": {"value": 0.0000000000000000}, 
+       "log.initial.gl514b.OrbAngMom": {"value": 5.7322265617904285e+40, "unit": (u.kg * u.m ** 2) / u.sec}, 
+       "log.initial.gl514b.LongP": {"value": 0.0000000000000000, "unit": u.rad}, 
+       "log.initial.gl514b.TotOrbEnergy": {"value": -1.6646340895690071e+34, "unit": u.Joule}, 
+       "log.initial.gl514b.OrbPotEnergy": {"value": -3.3292681791380142e+34, "unit": u.Joule}, 
+       "log.initial.gl514b.LostEnergy": {"value": 5.5626846462680035e-309, "unit": u.Joule}, 
+       "log.initial.gl514b.TidalRadius": {"value": 1.3394010000000000e+07, "unit": u.m}, 
+       "log.initial.gl514b.DsemiDtEqtide": {"value": 3.6833472670881283e-10, "unit": u.m / u.sec}, 
+       "log.initial.gl514b.DeccDtEqtide": {"value": 2.7299498000248231e-21, "unit": 1 / u.sec}, 
+       "log.initial.gl514b.DMeanMotionDtEqtide": {"value": -4.5392786228651742e-27, "unit": 1 / u.sec ** 2}, 
+       "log.initial.gl514b.DOrbPerDtEqtide": {"value": 1.0601935362604418e-13}, 
+       "log.initial.gl514b.EccTimeEqtide": {"value": 1.6483819592430168e+20, "unit": u.sec}, 
+       "log.initial.gl514b.SemiTimeEqtide": {"value": 1.7139383516588075e+20, "unit": u.sec}, 
+       "log.initial.gl514b.DHEccDtEqtide": {"value": 0.0000000000000000, "unit": 1 / u.sec}, 
+       "log.initial.gl514b.DKEccDtEqtide": {"value": 2.7299498000248231e-21, "unit": 1 / u.sec}, 
+       "log.initial.gl514b.DXoblDtEqtide": {"value": 1.3978299777162946e-16, "unit": 1 / u.sec}, 
+       "log.initial.gl514b.DYoblDtEqtide": {"value": 0.0000000000000000, "unit": 1 / u.sec}, 
+       "log.initial.gl514b.DZoblDtEqtide": {"value": -6.0779377238292621e-17, "unit": 1 / u.sec}, 
+       "log.initial.gl514b.LockTime": {"value": -1.0000000000000000, "unit": u.sec}, 
+       "log.initial.gl514b.BodyDsemiDtEqtide": {"value": -1.0000000000000000}, 
+       "log.initial.gl514b.BodyDeccDt": {"value": -1.0000000000000000}, 
+       "log.initial.gl514b.DOblDtEqtide": {"value": 1.5242512641785164e-16, "unit": u.rad / u.sec}, 
+       "log.initial.gl514b.DRotPerDtEqtide": {"value": 6.7687224755314529e-11}, 
+       "log.initial.gl514b.DRotRateDtEqtide": {"value": -5.6971688537695393e-20, "unit": 1 / u.sec ** 2}, 
+       "log.initial.gl514b.EqRotPer": {"value": 5.1067830269352868e+06, "unit": u.sec}, 
+       "log.initial.gl514b.EqTidePower": {"value": -1.0000000000000000, "unit": 1 / u.sec}, 
+       "log.initial.gl514b.OceanK2": {"value": 0.0100000000000000}, 
+       "log.initial.gl514b.EnvTidalQ": {"value": -1.0000000000000000}, 
+       "log.initial.gl514b.OceanTidalQ": {"value": -1.0000000000000000}, 
+       "log.initial.gl514b.TidalTau": {"value": 638.0000000000000000, "unit": u.sec}, 
+       "log.initial.gl514b.TideLock": {"value": 0.0000000000000000}, 
+       "log.initial.gl514b.RotTimeEqtide": {"value": 1.2764594842280960e+15, "unit": u.sec}, 
+       "log.initial.gl514b.EnvK2": {"value": 0.0100000000000000}, 
+       "log.initial.gl514b.OblTimeEqtide": {"value": -1.0000000000000000}, 
+       "log.initial.gl514b.PowerEqtide": {"value": 5.6730849058585900e+15, "unit": u.W}, 
+       "log.initial.gl514b.SurfEnFluxEqtide": {"value": 2.5164492816849116, "unit": u.kg / u.sec ** 3}, 
+       "log.initial.gl514b.TGlobal": {"value": 676.2679353212223532, "unit": u.sec}, 
+       "log.initial.gl514b.AlbedoGlobal": {"value": 0.2996778493802214}, 
+       "log.initial.gl514b.FluxInGlobal": {"value": 1046.1316028565233864, "unit": u.kg / u.sec ** 3}, 
+       "log.initial.gl514b.FluxOutGlobal": {"value": 1046.1299848213550376, "unit": u.kg / u.sec ** 3}, 
+       "log.initial.gl514b.TotIceMass": {"value": 0.0000000000000000, "unit": u.kg}, 
+       "log.initial.gl514b.TotIceFlow": {"value": 0.0000000000000000, "unit": u.kg}, 
+       "log.initial.gl514b.TotIceBalance": {"value": 0.0000000000000000, "unit": u.kg}, 
+       "log.initial.gl514b.SkipSeas": {"value": 0.0000000000000000}, 
+       "log.initial.gl514b.AreaIceCov": {"value": 0.0000000000000000}, 
+       "log.initial.gl514b.Latitude": {"value": -1.4552620265106593, "unit": u.rad}, 
+       "log.initial.gl514b.TempLat": {"value": 557.2247797713517912, "unit": u.sec}, 
+       "log.initial.gl514b.AlbedoLat": {"value": 0.3675649878041329}, 
+       "log.initial.gl514b.AnnInsol": {"value": 755.3407450595735781, "unit": u.kg / u.sec ** 3}, 
+       "log.initial.gl514b.FluxMerid": {"value": -2.0352457824211096e+16, "unit": u.Joule}, 
+       "log.initial.gl514b.FluxIn": {"value": 486.1623336259239636, "unit": u.kg / u.sec ** 3}, 
+       "log.initial.gl514b.FluxOut": {"value": 797.3297897221252697, "unit": u.kg / u.sec ** 3}, 
+       "log.initial.gl514b.DivFlux": {"value": -311.1690738881308107, "unit": u.kg / u.sec ** 3}, 
+       "log.initial.gl514b.IceMass": {"value": 0.0000000000000000}, 
+       "log.initial.gl514b.IceHeight": {"value": 0.0000000000000000, "unit": u.m}, 
+       "log.initial.gl514b.DIceMassDt": {"value": 0.0000000000000000, "unit": u.m}, 
+       "log.initial.gl514b.IceFlow": {"value": 0.0000000000000000}, 
+       "log.initial.gl514b.EnergyResL": {"value": -4.7634785005357116e-11, "unit": u.kg / u.sec ** 3}, 
+       "log.initial.gl514b.EnergyResW": {"value": 4.0836312109604478e-10, "unit": u.kg / u.sec ** 3}, 
+       "log.initial.gl514b.BedrockH": {"value": 0.0000000000000000, "unit": u.m}, 
+       "log.initial.gl514b.TempLandLat": {"value": 545.8662422200939091, "unit": u.sec}, 
+       "log.initial.gl514b.TempWaterLat": {"value": 563.0761476007876354, "unit": u.sec}, 
+       "log.initial.gl514b.AlbedoLandLat": {"value": 0.4335649878041328}, 
+       "log.initial.gl514b.AlbedoWaterLat": {"value": 0.3335649878041328}, 
+       "log.initial.gl514b.TempMinLat": {"value": 520.6583069986143073, "unit": u.sec}, 
+       "log.initial.gl514b.TempMaxLat": {"value": 598.1727520370602633, "unit": u.sec}, 
+       "log.initial.gl514b.Snowball": {"value": 0.0000000000000000}, 
+       "log.initial.gl514b.PlanckBAvg": {"value": 2.0899999999999990}, 
+       "log.initial.gl514b.IceAccum": {"value": 0.0000000000000000}, 
+       "log.initial.gl514b.IceAblate": {"value": 0.0000000000000000}, 
+       "log.initial.gl514b.TempMaxLand": {"value": 657.9976852023412448, "unit": u.sec}, 
+       "log.initial.gl514b.TempMaxWater": {"value": 568.6086486017817379, "unit": u.sec}, 
+       "log.initial.gl514b.PeakInsol": {"value": 4243.8792145061070187, "unit": u.kg / u.sec ** 3}, 
+       "log.initial.gl514b.IceCapNorthLand": {"value": 0.0000000000000000}, 
+       "log.initial.gl514b.IceCapNorthSea": {"value": 0.0000000000000000}, 
+       "log.initial.gl514b.IceCapSouthLand": {"value": 0.0000000000000000}, 
+       "log.initial.gl514b.IceCapSouthSea": {"value": 0.0000000000000000}, 
+       "log.initial.gl514b.IceBeltLand": {"value": 0.0000000000000000}, 
+       "log.initial.gl514b.IceBeltSea": {"value": 0.0000000000000000}, 
+       "log.initial.gl514b.SnowballLand": {"value": 0.0000000000000000}, 
+       "log.initial.gl514b.SnowballSea": {"value": 0.0000000000000000}, 
+       "log.initial.gl514b.IceFree": {"value": 1.0000000000000000}, 
+       "log.initial.gl514b.IceCapNorthLatLand": {"value": 100.0000000000000000, "unit": u.rad}, 
+       "log.initial.gl514b.IceCapNorthLatSea": {"value": 100.0000000000000000, "unit": u.rad}, 
+       "log.initial.gl514b.IceCapSouthLatLand": {"value": 100.0000000000000000, "unit": u.rad}, 
+       "log.initial.gl514b.IceCapSouthLatSea": {"value": 100.0000000000000000, "unit": u.rad}, 
+       "log.initial.gl514b.IceBeltNorthLatLand": {"value": 100.0000000000000000, "unit": u.rad}, 
+       "log.initial.gl514b.IceBeltNorthLatSea": {"value": 100.0000000000000000, "unit": u.rad}, 
+       "log.initial.gl514b.IceBeltSouthLatLand": {"value": 100.0000000000000000, "unit": u.rad}, 
+       "log.initial.gl514b.IceBeltSouthLatSea": {"value": 100.0000000000000000, "unit": u.rad}, 
+       "log.final.system.Age": {"value": 3.1557600000000000e+07, "unit": u.sec, "rtol": 1e-4}, 
+       "log.final.system.Time": {"value": 3.1557600000000000e+07, "unit": u.sec, "rtol": 1e-4}, 
+       "log.final.system.TotAngMom": {"value": 8.9534036048945070e+41, "unit": (u.kg * u.m ** 2) / u.sec, "rtol": 1e-4}, 
+       "log.final.system.TotEnergy": {"value": -3.1347047276701209e+40, "unit": u.Joule, "rtol": 1e-4}, 
+       "log.final.system.PotEnergy": {"value": -3.1348025171406187e+40, "unit": u.Joule, "rtol": 1e-4}, 
+       "log.final.system.KinEnergy": {"value": 9.9742042395602989e+35, "unit": u.Joule, "rtol": 1e-4}, 
+       "log.final.system.DeltaTime": {"value": 3.1557600000000000e+07, "unit": u.sec, "rtol": 1e-4}, 
+       "log.final.gl514.Mass": {"value": 1.0140921600000001e+30, "unit": u.kg, "rtol": 1e-4}, 
+       "log.final.gl514.Obliquity": {"value": 0.0000000000000000, "unit": u.rad, "rtol": 1e-4}, 
+       "log.final.gl514.PrecA": {"value": 0.0000000000000000, "unit": u.rad, "rtol": 1e-4}, 
+       "log.final.gl514.Xobl": {"value": 0.0000000000000000, "rtol": 1e-4}, 
+       "log.final.gl514.Yobl": {"value": 0.0000000000000000, "rtol": 1e-4}, 
+       "log.final.gl514.Zobl": {"value": 1.0000000000000000, "rtol": 1e-4}, 
+       "log.final.gl514.Radius": {"value": 1.3137125118741000e+09, "unit": u.m, "rtol": 1e-4}, 
+       "log.final.gl514.RadGyra": {"value": 0.4484976455000000, "rtol": 1e-4}, 
+       "log.final.gl514.RotAngMom": {"value": 8.3801799123461292e+41, "unit": (u.kg * u.m ** 2) / u.sec, "rtol": 1e-4}, 
+       "log.final.gl514.RotKinEnergy": {"value": 9.9742042395602989e+35, "unit": u.Joule, "rtol": 1e-4}, 
+       "log.final.gl514.RotVel": {"value": 3127.1970393365104428, "unit": u.m / u.sec, "rtol": 1e-4}, 
+       "log.final.gl514.BodyType": {"value": 0.0000000000000000, "rtol": 1e-4}, 
+       "log.final.gl514.RotRate": {"value": 2.3804272328009969e-06, "unit": 1 / u.sec, "rtol": 1e-4}, 
+       "log.final.gl514.RotPer": {"value": 30.5500000858501117, "unit": u.day, "rtol": 1e-4}, 
+       "log.final.gl514.Density": {"value": 106.7794814169220103, "unit": u.kg / u.m ** 3, "rtol": 1e-4}, 
+       "log.final.gl514.SurfEnFluxTotal": {"value": 1.9756860885353259e-10, "unit": u.kg / u.sec ** 3, "rtol": 1e-4}, 
+       "log.final.gl514.ImK2": {"value": -4.9999999999999998e-07, "rtol": 1e-4}, 
+       "log.final.gl514.K2": {"value": 0.5000000000000000, "rtol": 1e-4}, 
+       "log.final.gl514.K2Man": {"value": 0.0100000000000000, "rtol": 1e-4}, 
+       "log.final.gl514.Imk2Man": {"value": 0.0000000000000000, "rtol": 1e-4}, 
+       "log.final.gl514.TidalQMantle": {"value": 100.0000000000000000, "rtol": 1e-4}, 
+       "log.final.gl514.HEcc": {"value": 0.0000000000000000, "rtol": 1e-4}, 
+       "log.final.gl514.HZLimitDryRunaway": {"value": 1.1215301622470532e+11, "unit": u.m, "rtol": 1e-4}, 
+       "log.final.gl514.HZLimRecVenus": {"value": 9.9104353063518631e+10, "unit": u.m, "rtol": 1e-4}, 
+       "log.final.gl514.HZLimRunaway": {"value": 1.3136272247075568e+11, "unit": u.m, "rtol": 1e-4}, 
+       "log.final.gl514.HZLimMoistGreenhouse": {"value": 1.3113105818748537e+11, "unit": u.m, "rtol": 1e-4}, 
+       "log.final.gl514.HZLimMaxGreenhouse": {"value": 2.4324580025875378e+11, "unit": u.m, "rtol": 1e-4}, 
+       "log.final.gl514.HZLimEarlyMars": {"value": 2.6530630489685046e+11, "unit": u.m, "rtol": 1e-4}, 
+       "log.final.gl514.Instellation": {"value": -1.0000000000000000, "unit": u.kg / u.sec ** 3, "rtol": 1e-4}, 
+       "log.final.gl514.KEcc": {"value": 0.0000000000000000, "rtol": 1e-4}, 
+       "log.final.gl514.Eccentricity": {"value": -1.0000000000000000, "rtol": 1e-4}, 
+       "log.final.gl514.OrbEnergy": {"value": 0.0000000000000000, "unit": u.Joule, "rtol": 1e-4}, 
+       "log.final.gl514.MeanMotion": {"value": -1.0000000000000000, "unit": 1 / u.sec, "rtol": 1e-4}, 
+       "log.final.gl514.OrbPeriod": {"value": -1.0000000000000000, "unit": u.sec, "rtol": 1e-4}, 
+       "log.final.gl514.SemiMajorAxis": {"value": -1.0000000000000000, "unit": u.m, "rtol": 1e-4}, 
+       "log.final.gl514.CriticalSemiMajorAxis": {"value": -1.0000000000000000, "unit": u.m, "rtol": 1e-4}, 
+       "log.final.gl514.COPP": {"value": 0.0000000000000000, "rtol": 1e-4}, 
+       "log.final.gl514.OrbAngMom": {"value": 0.0000000000000000, "unit": (u.kg * u.m ** 2) / u.sec, "rtol": 1e-4}, 
+       "log.final.gl514.LongP": {"value": 0.0000000000000000, "unit": u.rad, "rtol": 1e-4}, 
+       "log.final.gl514.LXUVTot": {"value": 2.6238552211501671e+23, "unit": u.kg / u.sec ** 3, "rtol": 1e-4}, 
+       "log.final.gl514.TotOrbEnergy": {"value": -1.6646340895687004e+34, "unit": u.Joule, "rtol": 1e-4}, 
+       "log.final.gl514.OrbPotEnergy": {"value": -1.0000000000000000, "unit": u.Joule, "rtol": 1e-4}, 
+       "log.final.gl514.LostEnergy": {"value": 5.6058038407261769e+27, "unit": u.Joule, "rtol": 1e-4}, 
+       "log.final.gl514.LostAngMom": {"value": 2.3549570242977578e+33, "unit": (u.kg * u.m ** 2) / u.sec, "rtol": 1e-4}, 
+       "log.final.gl514.LockTime": {"value": -1.0000000000000000, "unit": u.sec, "rtol": 1e-4}, 
+       "log.final.gl514.BodyDsemiDtEqtide": {"value": -1.0000000000000000, "rtol": 1e-4}, 
+       "log.final.gl514.BodyDeccDt": {"value": -1.0000000000000000, "rtol": 1e-4}, 
+       "log.final.gl514.DOblDtEqtide": {"value": 0.0000000000000000, "unit": u.rad / u.sec, "rtol": 1e-4}, 
+       "log.final.gl514.DRotPerDtEqtide": {"value": 9.6722127498436088e-21, "rtol": 1e-4}, 
+       "log.final.gl514.DRotRateDtEqtide": {"value": -8.7227975413983740e-33, "unit": 1 / u.sec ** 2, "rtol": 1e-4}, 
+       "log.final.gl514.EqRotPer": {"value": 5.0877011498425016e+06, "unit": u.sec, "rtol": 1e-4}, 
+       "log.final.gl514.EqTidePower": {"value": -1.0000000000000000, "unit": 1 / u.sec, "rtol": 1e-4}, 
+       "log.final.gl514.OceanK2": {"value": 0.0100000000000000, "rtol": 1e-4}, 
+       "log.final.gl514.EnvTidalQ": {"value": -1.0000000000000000, "rtol": 1e-4}, 
+       "log.final.gl514.OceanTidalQ": {"value": -1.0000000000000000, "rtol": 1e-4}, 
+       "log.final.gl514.TidalTau": {"value": 0.1000000000000000, "unit": u.sec, "rtol": 1e-4}, 
+       "log.final.gl514.TideLock": {"value": 0.0000000000000000, "rtol": 1e-4}, 
+       "log.final.gl514.RotTimeEqtide": {"value": 2.7289722379815599e+26, "unit": u.sec, "rtol": 1e-4}, 
+       "log.final.gl514.EnvK2": {"value": 0.0100000000000000, "rtol": 1e-4}, 
+       "log.final.gl514.OblTimeEqtide": {"value": -1.0000000000000000, "rtol": 1e-4}, 
+       "log.final.gl514.PowerEqtide": {"value": 4.2847795070601802e+09, "unit": u.W, "rtol": 1e-4}, 
+       "log.final.gl514.SurfEnFluxEqtide": {"value": 1.9756860885353259e-10, "unit": u.kg / u.sec ** 3, "rtol": 1e-4}, 
+       "log.final.gl514.Luminosity": {"value": 2.6238552211501671e+26, "unit": u.W, "rtol": 1e-4}, 
+       "log.final.gl514.LXUVStellar": {"value": 2.6238552211501671e+23, "unit": u.W, "rtol": 1e-4}, 
+       "log.final.gl514.Temperature": {"value": 3818.8020013755422042, "unit": u.K, "rtol": 1e-4}, 
+       "log.final.gl514.LXUVFrac": {"value": 0.0010000000000000, "rtol": 1e-4}, 
+       "log.final.gl514.RossbyNumber": {"value": 0.6873965303916811, "rtol": 1e-4}, 
+       "log.final.gl514.DRotPerDtStellar": {"value": 2.3504479143348279e-10, "rtol": 1e-4}, 
+       "log.final.gl514b.Mass": {"value": 3.1053088098643201e+25, "unit": u.kg, "rtol": 1e-4}, 
+       "log.final.gl514b.Obliquity": {"value": 0.4101523790171157, "unit": u.rad, "rtol": 1e-4}, 
+       "log.final.gl514b.PrecA": {"value": 0.0000000000000000, "unit": u.rad, "rtol": 1e-4}, 
+       "log.final.gl514b.Xobl": {"value": 0.3987490733257115, "rtol": 1e-4}, 
+       "log.final.gl514b.Yobl": {"value": 0.0000000000000000, "rtol": 1e-4}, 
+       "log.final.gl514b.Zobl": {"value": 0.9170600724717473, "rtol": 1e-4}, 
+       "log.final.gl514b.Radius": {"value": 1.3394010000000000e+07, "unit": u.m, "rtol": 1e-4}, 
+       "log.final.gl514b.RadGyra": {"value": 0.5000000000000000, "rtol": 1e-4}, 
+       "log.final.gl514b.RotAngMom": {"value": 1.0128197393087430e+35, "unit": (u.kg * u.m ** 2) / u.sec, "rtol": 1e-4}, 
+       "log.final.gl514b.RotKinEnergy": {"value": 3.6827164048131317e+30, "unit": u.Joule, "rtol": 1e-4}, 
+       "log.final.gl514b.RotVel": {"value": 974.0398698567373685, "unit": u.m / u.sec, "rtol": 1e-4}, 
+       "log.final.gl514b.BodyType": {"value": 0.0000000000000000, "rtol": 1e-4}, 
+       "log.final.gl514b.RotRate": {"value": 7.2722050368540664e-05, "unit": 1 / u.sec, "rtol": 1e-4}, 
+       "log.final.gl514b.RotPer": {"value": 8.6400002136046387e+04, "unit": u.sec, "rtol": 1e-4}, 
+       "log.final.gl514b.Density": {"value": 3085.2071528551509800, "unit": u.kg / u.m ** 3, "rtol": 1e-4}, 
+       "log.final.gl514b.SurfEnFluxTotal": {"value": 2.5164491503354709, "unit": u.kg / u.sec ** 3, "rtol": 1e-4}, 
+       "log.final.gl514b.ImK2": {"value": -0.0250000000000000, "rtol": 1e-4}, 
+       "log.final.gl514b.K2": {"value": 0.3000000000000000, "rtol": 1e-4}, 
+       "log.final.gl514b.K2Man": {"value": 0.0100000000000000, "rtol": 1e-4}, 
+       "log.final.gl514b.Imk2Man": {"value": 0.0000000000000000, "rtol": 1e-4}, 
+       "log.final.gl514b.TidalQMantle": {"value": 100.0000000000000000, "rtol": 1e-4}, 
+       "log.final.gl514b.HEcc": {"value": 0.0000000000000000, "rtol": 1e-4}, 
+       "log.final.gl514b.HZLimitDryRunaway": {"value": 1.0509802225543790e+11, "unit": u.m, "rtol": 1e-4}, 
+       "log.final.gl514b.HZLimRecVenus": {"value": 9.9104353063518631e+10, "unit": u.m, "rtol": 1e-4}, 
+       "log.final.gl514b.HZLimRunaway": {"value": 1.3136272247075568e+11, "unit": u.m, "rtol": 1e-4}, 
+       "log.final.gl514b.HZLimMoistGreenhouse": {"value": 1.3113105818748537e+11, "unit": u.m, "rtol": 1e-4}, 
+       "log.final.gl514b.HZLimMaxGreenhouse": {"value": 2.4324580025875378e+11, "unit": u.m, "rtol": 1e-4}, 
+       "log.final.gl514b.HZLimEarlyMars": {"value": 2.6530630489685046e+11, "unit": u.m, "rtol": 1e-4}, 
+       "log.final.gl514b.Instellation": {"value": 5866.6331509117835594, "unit": u.kg / u.sec ** 3, "rtol": 1e-4}, 
+       "log.final.gl514b.KEcc": {"value": 0.4500000000000862, "rtol": 1e-4}, 
+       "log.final.gl514b.Eccentricity": {"value": 0.4500000000000862, "rtol": 1e-4}, 
+       "log.final.gl514b.OrbEnergy": {"value": -1.6646340895687004e+34, "unit": u.Joule, "rtol": 1e-4}, 
+       "log.final.gl514b.MeanMotion": {"value": 5.1866958137276320e-07, "unit": 1 / u.sec, "rtol": 1e-4}, 
+       "log.final.gl514b.OrbPeriod": {"value": 1.2114042413186973e+07, "unit": u.sec, "rtol": 1e-4}, 
+       "log.final.gl514b.SemiMajorAxis": {"value": 6.3130301435411629e+10, "unit": u.m, "rtol": 1e-4}, 
+       "log.final.gl514b.CriticalSemiMajorAxis": {"value": -1.0000000000000000, "unit": u.m, "rtol": 1e-4}, 
+       "log.final.gl514b.COPP": {"value": 0.0000000000000000, "rtol": 1e-4}, 
+       "log.final.gl514b.OrbAngMom": {"value": 5.7322265617906780e+40, "unit": (u.kg * u.m ** 2) / u.sec, "rtol": 1e-4}, 
+       "log.final.gl514b.LongP": {"value": 0.0000000000000000, "unit": u.rad, "rtol": 1e-4}, 
+       "log.final.gl514b.TotOrbEnergy": {"value": -1.6646340895687004e+34, "unit": u.Joule, "rtol": 1e-4}, 
+       "log.final.gl514b.OrbPotEnergy": {"value": -3.3292681791374008e+34, "unit": u.Joule, "rtol": 1e-4}, 
+       "log.final.gl514b.LostEnergy": {"value": 1.7902893955276926e+23, "unit": u.Joule, "rtol": 1e-4}, 
+       "log.final.gl514b.TidalRadius": {"value": 1.3394010000000000e+07, "unit": u.m, "rtol": 1e-4}, 
+       "log.final.gl514b.DsemiDtEqtide": {"value": 3.6833471661025376e-10, "unit": u.m / u.sec, "rtol": 1e-4}, 
+       "log.final.gl514b.DeccDtEqtide": {"value": 2.7299497249402080e-21, "unit": 1 / u.sec, "rtol": 1e-4}, 
+       "log.final.gl514b.DMeanMotionDtEqtide": {"value": -4.5392784984105951e-27, "unit": 1 / u.sec ** 2, "rtol": 1e-4}, 
+       "log.final.gl514b.DOrbPerDtEqtide": {"value": 1.0601935071934267e-13, "rtol": 1e-4}, 
+       "log.final.gl514b.EccTimeEqtide": {"value": 1.6483820045804768e+20, "unit": u.sec, "rtol": 1e-4}, 
+       "log.final.gl514b.SemiTimeEqtide": {"value": 1.7139383986498273e+20, "unit": u.sec, "rtol": 1e-4}, 
+       "log.final.gl514b.DHEccDtEqtide": {"value": 0.0000000000000000, "unit": 1 / u.sec, "rtol": 1e-4}, 
+       "log.final.gl514b.DKEccDtEqtide": {"value": 2.7299497249402080e-21, "unit": 1 / u.sec, "rtol": 1e-4}, 
+       "log.final.gl514b.DXoblDtEqtide": {"value": 1.3978299858905663e-16, "unit": 1 / u.sec, "rtol": 1e-4}, 
+       "log.final.gl514b.DYoblDtEqtide": {"value": 0.0000000000000000, "unit": 1 / u.sec, "rtol": 1e-4}, 
+       "log.final.gl514b.DZoblDtEqtide": {"value": -6.0779378393221619e-17, "unit": 1 / u.sec, "rtol": 1e-4}, 
+       "log.final.gl514b.LockTime": {"value": -1.0000000000000000, "unit": u.sec, "rtol": 1e-4}, 
+       "log.final.gl514b.BodyDsemiDtEqtide": {"value": -1.0000000000000000, "rtol": 1e-4}, 
+       "log.final.gl514b.BodyDeccDt": {"value": -1.0000000000000000, "rtol": 1e-4}, 
+       "log.final.gl514b.DOblDtEqtide": {"value": 1.5242512762800832e-16, "unit": u.rad / u.sec, "rtol": 1e-4}, 
+       "log.final.gl514b.DRotPerDtEqtide": {"value": 6.7687226270763294e-11, "rtol": 1e-4}, 
+       "log.final.gl514b.DRotRateDtEqtide": {"value": -5.6971686996239482e-20, "unit": 1 / u.sec ** 2, "rtol": 1e-4}, 
+       "log.final.gl514b.EqRotPer": {"value": 5.1067830278577860e+06, "unit": u.sec, "rtol": 1e-4}, 
+       "log.final.gl514b.EqTidePower": {"value": -1.0000000000000000, "unit": 1 / u.sec, "rtol": 1e-4}, 
+       "log.final.gl514b.OceanK2": {"value": 0.0100000000000000, "rtol": 1e-4}, 
+       "log.final.gl514b.EnvTidalQ": {"value": -1.0000000000000000, "rtol": 1e-4}, 
+       "log.final.gl514b.OceanTidalQ": {"value": -1.0000000000000000, "rtol": 1e-4}, 
+       "log.final.gl514b.TidalTau": {"value": 638.0000000000000000, "unit": u.sec, "rtol": 1e-4}, 
+       "log.final.gl514b.TideLock": {"value": 0.0000000000000000, "rtol": 1e-4}, 
+       "log.final.gl514b.RotTimeEqtide": {"value": 1.2764594872070545e+15, "unit": u.sec, "rtol": 1e-4}, 
+       "log.final.gl514b.EnvK2": {"value": 0.0100000000000000, "rtol": 1e-4}, 
+       "log.final.gl514b.OblTimeEqtide": {"value": -1.0000000000000000, "rtol": 1e-4}, 
+       "log.final.gl514b.PowerEqtide": {"value": 5.6730846097443250e+15, "unit": u.W, "rtol": 1e-4}, 
+       "log.final.gl514b.SurfEnFluxEqtide": {"value": 2.5164491503354709, "unit": u.kg / u.sec ** 3, "rtol": 1e-4}, 
+       "log.final.gl514b.TGlobal": {"value": 676.2684488417830835, "unit": u.sec, "rtol": 1e-4}, 
+       "log.final.gl514b.AlbedoGlobal": {"value": 0.2996778494395087, "rtol": 1e-4}, 
+       "log.final.gl514b.FluxInGlobal": {"value": 1046.1316028669016305, "unit": u.kg / u.sec ** 3, "rtol": 1e-4}, 
+       "log.final.gl514b.FluxOutGlobal": {"value": 1046.1310580793265217, "unit": u.kg / u.sec ** 3, "rtol": 1e-4}, 
+       "log.final.gl514b.TotIceMass": {"value": 0.0000000000000000, "unit": u.kg, "rtol": 1e-4}, 
+       "log.final.gl514b.TotIceFlow": {"value": 0.0000000000000000, "unit": u.kg, "rtol": 1e-4}, 
+       "log.final.gl514b.TotIceBalance": {"value": 0.0000000000000000, "unit": u.kg, "rtol": 1e-4}, 
+       "log.final.gl514b.SkipSeas": {"value": 0.0000000000000000, "rtol": 1e-4}, 
+       "log.final.gl514b.AreaIceCov": {"value": 0.0000000000000000, "rtol": 1e-4}, 
+       "log.final.gl514b.Latitude": {"value": -1.4552620265106593, "unit": u.rad, "rtol": 1e-4}, 
+       "log.final.gl514b.TempLat": {"value": 557.2252940986816157, "unit": u.sec, "rtol": 1e-4}, 
+       "log.final.gl514b.AlbedoLat": {"value": 0.3675649876307175, "rtol": 1e-4}, 
+       "log.final.gl514b.AnnInsol": {"value": 755.3407530704208739, "unit": u.kg / u.sec ** 3, "rtol": 1e-4}, 
+       "log.final.gl514b.FluxMerid": {"value": -2.0352457616350744e+16, "unit": u.Joule, "rtol": 1e-4}, 
+       "log.final.gl514b.FluxIn": {"value": 486.1623390295726495, "unit": u.kg / u.sec ** 3, "rtol": 1e-4}, 
+       "log.final.gl514b.FluxOut": {"value": 797.3308646662439969, "unit": u.kg / u.sec ** 3, "rtol": 1e-4}, 
+       "log.final.gl514b.DivFlux": {"value": -311.1690707101504927, "unit": u.kg / u.sec ** 3, "rtol": 1e-4}, 
+       "log.final.gl514b.IceMass": {"value": 0.0000000000000000, "rtol": 1e-4}, 
+       "log.final.gl514b.IceHeight": {"value": 0.0000000000000000, "unit": u.m, "rtol": 1e-4}, 
+       "log.final.gl514b.DIceMassDt": {"value": 0.0000000000000000, "unit": u.m, "rtol": 1e-4}, 
+       "log.final.gl514b.IceFlow": {"value": 0.0000000000000000, "rtol": 1e-4}, 
+       "log.final.gl514b.EnergyResL": {"value": 4.1836756281554699e-11, "unit": u.kg / u.sec ** 3, "rtol": 1e-4}, 
+       "log.final.gl514b.EnergyResW": {"value": -3.7482550396816805e-10, "unit": u.kg / u.sec ** 3, "rtol": 1e-4}, 
+       "log.final.gl514b.BedrockH": {"value": 0.0000000000000000, "unit": u.m, "rtol": 1e-4}, 
+       "log.final.gl514b.TempLandLat": {"value": 545.8665757148477269, "unit": u.sec, "rtol": 1e-4}, 
+       "log.final.gl514b.TempWaterLat": {"value": 563.0767550842926994, "unit": u.sec, "rtol": 1e-4}, 
+       "log.final.gl514b.AlbedoLandLat": {"value": 0.4335649876307174, "rtol": 1e-4}, 
+       "log.final.gl514b.AlbedoWaterLat": {"value": 0.3335649876307175, "rtol": 1e-4}, 
+       "log.final.gl514b.TempMinLat": {"value": 520.6588002187510256, "unit": u.sec, "rtol": 1e-4}, 
+       "log.final.gl514b.TempMaxLat": {"value": 598.1732898182863210, "unit": u.sec, "rtol": 1e-4}, 
+       "log.final.gl514b.Snowball": {"value": 0.0000000000000000, "rtol": 1e-4}, 
+       "log.final.gl514b.PlanckBAvg": {"value": 2.0899999999999990, "rtol": 1e-4}, 
+       "log.final.gl514b.IceAccum": {"value": 0.0000000000000000, "rtol": 1e-4}, 
+       "log.final.gl514b.IceAblate": {"value": 0.0000000000000000, "rtol": 1e-4}, 
+       "log.final.gl514b.TempMaxLand": {"value": 657.9980346867481558, "unit": u.sec, "rtol": 1e-4}, 
+       "log.final.gl514b.TempMaxWater": {"value": 568.6092765386927113, "unit": u.sec, "rtol": 1e-4}, 
+       "log.final.gl514b.PeakInsol": {"value": 4243.8792614541798685, "unit": u.kg / u.sec ** 3, "rtol": 1e-4}, 
+       "log.final.gl514b.IceCapNorthLand": {"value": 0.0000000000000000, "rtol": 1e-4}, 
+       "log.final.gl514b.IceCapNorthSea": {"value": 0.0000000000000000, "rtol": 1e-4}, 
+       "log.final.gl514b.IceCapSouthLand": {"value": 0.0000000000000000, "rtol": 1e-4}, 
+       "log.final.gl514b.IceCapSouthSea": {"value": 0.0000000000000000, "rtol": 1e-4}, 
+       "log.final.gl514b.IceBeltLand": {"value": 0.0000000000000000, "rtol": 1e-4}, 
+       "log.final.gl514b.IceBeltSea": {"value": 0.0000000000000000, "rtol": 1e-4}, 
+       "log.final.gl514b.SnowballLand": {"value": 0.0000000000000000, "rtol": 1e-4}, 
+       "log.final.gl514b.SnowballSea": {"value": 0.0000000000000000, "rtol": 1e-4}, 
+       "log.final.gl514b.IceFree": {"value": 1.0000000000000000, "rtol": 1e-4}, 
+       "log.final.gl514b.IceCapNorthLatLand": {"value": 100.0000000000000000, "unit": u.rad, "rtol": 1e-4}, 
+       "log.final.gl514b.IceCapNorthLatSea": {"value": 100.0000000000000000, "unit": u.rad, "rtol": 1e-4}, 
+       "log.final.gl514b.IceCapSouthLatLand": {"value": 100.0000000000000000, "unit": u.rad, "rtol": 1e-4}, 
+       "log.final.gl514b.IceCapSouthLatSea": {"value": 100.0000000000000000, "unit": u.rad, "rtol": 1e-4}, 
+       "log.final.gl514b.IceBeltNorthLatLand": {"value": 100.0000000000000000, "unit": u.rad, "rtol": 1e-4}, 
+       "log.final.gl514b.IceBeltNorthLatSea": {"value": 100.0000000000000000, "unit": u.rad, "rtol": 1e-4}, 
+       "log.final.gl514b.IceBeltSouthLatLand": {"value": 100.0000000000000000, "unit": u.rad, "rtol": 1e-4}, 
+       "log.final.gl514b.IceBeltSouthLatSea": {"value": 100.0000000000000000, "unit": u.rad, "rtol": 1e-4}, 
+   } 
+)
+class TestCTL_POISE(Benchmark): 
+   pass 
diff --git a/tests/CTL_POISE/vpl.in b/tests/CTL_POISE/vpl.in
new file mode 100644
index 000000000..d152c18b6
--- /dev/null
+++ b/tests/CTL_POISE/vpl.in
@@ -0,0 +1,16 @@
+sSystemName   gl514
+iVerbose      5
+bOverwrite    1
+sUnitMass     solar
+sUnitLength   au
+sUnitTime     y
+sUnitAngle    d
+sUnitTemp     K
+bDoLog        1
+iDigits       16
+saBodyFiles   gl514.in gl514b.in
+bDoForward    1
+bVarDt        1
+dEta          0.1
+dStopTime     1
+dOutputTime   1
diff --git a/tests/maketest.py b/tests/maketest.py
index fd66aaf8e..2423914eb 100644
--- a/tests/maketest.py
+++ b/tests/maketest.py
@@ -41,6 +41,7 @@ def Main(dir, initial=False):
         if dirname not in initial_list:
             initial = False
 
+        # Should check for trailing "/" and remove it XXX
         # grab infiles
         infiles = []
         for file in os.listdir(dirname):
@@ -438,10 +439,7 @@ def WriteTest(data, dirname, stellar):
     for i in badchars:
         dirname = dirname.replace(i, "")
 
-    print(stellar)
-
     test_file = "test_" + dirname + ".py"
-    print(test_file)
     with open(test_file, "w") as t:
         t.write("from benchmark import Benchmark, benchmark \n")
         t.write("import astropy.units as u \n")
@@ -513,6 +511,9 @@ def WriteTest(data, dirname, stellar):
         t.write("   pass")
         t.write(" \n")
 
+    print("Successfuly created new test file: "+dirname+"/"+test_file)
+
+
 
 def Arguments():
     parser = argparse.ArgumentParser(