diff --git a/.travis.yml b/.travis.yml
index 8d7303f6..4bfea1b0 100644
--- a/.travis.yml
+++ b/.travis.yml
@@ -15,7 +15,7 @@ before_install:
   - pip install sphinx_rtd_theme
   - pip install recommonmark
   - pip install matplotlib
-  - python setup.py install
+  - pip install -e .
 
 install: 
 
diff --git a/aronnax.f90 b/aronnax.f90
index ff3cdcc9..13521870 100644
--- a/aronnax.f90
+++ b/aronnax.f90
@@ -116,7 +116,7 @@ program aronnax
 
   namelist /GRID/ nx, ny, layers, dx, dy, fUfile, fVfile, wetMaskFile
 
-  namelist /INITIAL_CONDITONS/ initUfile, initVfile, initHfile, initEtaFile
+  namelist /INITIAL_CONDITIONS/ initUfile, initVfile, initHfile, initEtaFile
 
   namelist /EXTERNAL_FORCING/ zonalWindFile, meridionalWindFile, &
       DumpWind, wind_mag_time_series_file
@@ -127,7 +127,7 @@ program aronnax
   read(unit=8, nml=SPONGE)
   read(unit=8, nml=PHYSICS)
   read(unit=8, nml=GRID)
-  read(unit=8, nml=INITIAL_CONDITONS)
+  read(unit=8, nml=INITIAL_CONDITIONS)
   read(unit=8, nml=EXTERNAL_FORCING)
   close(unit=8)
 
diff --git a/aronnax/core.py b/aronnax/core.py
index 2015026f..039aed94 100644
--- a/aronnax/core.py
+++ b/aronnax/core.py
@@ -9,6 +9,7 @@
 
 from contextlib import contextmanager
 import os.path as p
+import re
 
 import numpy as np
 from scipy.io import FortranFile
@@ -94,7 +95,7 @@ def interpret_raw_file(name, nx, ny, layers):
 
 ### General input construction helpers
 
-def write_initial_heights(grid, h_funcs):
+def depths(grid, *h_funcs):
     X,Y = np.meshgrid(grid.x, grid.y)
     initH = np.ones((len(h_funcs), grid.ny, grid.nx))
     for i, f in enumerate(h_funcs):
@@ -102,53 +103,124 @@ def write_initial_heights(grid, h_funcs):
             initH[i,:,:] = f
         else:
             initH[i,:,:] = f(X, Y)
-    with fortran_file('initH.bin', 'w') as f:
-        f.write_record(initH.astype(np.float64))
+    return initH
 
-def write_wind_x(grid, func):
+def wind_x(grid, func):
     X,Y = np.meshgrid(grid.xp1, grid.y)
     if isinstance(func, (int, long, float)):
         wind_x = np.ones(grid.ny, grid.nx+1) * func
     else:
         wind_x = func(X, Y)
-    with fortran_file('wind_x.bin', 'w') as f:
-        f.write_record(wind_x.astype(np.float64))
+    return wind_x
 
-def write_wind_y(grid, func):
+def wind_y(grid, func):
     X,Y = np.meshgrid(grid.y, grid.xp1)
     if isinstance(func, (int, long, float)):
         wind_y = np.ones(grid.ny+1, grid.nx) * func
     else:
         wind_y = func(X, Y)
-    with fortran_file('wind_y.bin', 'w') as f:
-        f.write_record(wind_y.astype(np.float64))
+    return wind_y
 
 ### Specific construction helpers
 
-def write_f_plane(nx, ny, coeff):
-    """Write files defining an f-plane approximation to the Coriolis force."""
-    with fortran_file('fu.bin', 'w') as f:
-        f.write_record(np.ones((nx+1, ny), dtype=np.float64) * coeff)
-    with fortran_file('fv.bin', 'w') as f:
-        f.write_record(np.ones((nx, ny+1), dtype=np.float64) * coeff)
-
-def write_beta_plane(grid, f0, beta):
-    """Write files defining a beta-plane approximation to the Coriolis force."""
-    with fortran_file('fu.bin', 'w') as f:
-        _, Y = np.meshgrid(grid.xp1, grid.y)
-        fu = f0 + Y*beta
-        f.write_record(fu.astype(np.float64))
-    with fortran_file('fv.bin', 'w') as f:
-        _, Y = np.meshgrid(grid.x, grid.yp1)
-        fv = f0 + Y*beta
-        f.write_record(fv.astype(np.float64))
-
-def write_rectangular_pool(nx, ny):
-    """Write the wet mask file for a maximal rectangular pool."""
-    with fortran_file('wetmask.bin', 'w') as f:
-        wetmask = np.ones((nx, ny), dtype=np.float64)
-        wetmask[ 0, :] = 0
-        wetmask[-1, :] = 0
-        wetmask[ :, 0] = 0
-        wetmask[ :,-1] = 0
-        f.write_record(wetmask)
+def f_plane_f_u(grid, coeff):
+    """Define an f-plane approximation to the Coriolis force (u location)."""
+    return np.ones((grid.nx+1, grid.ny), dtype=np.float64) * coeff
+
+def f_plane_f_v(grid, coeff):
+    """Define an f-plane approximation to the Coriolis force (v location)."""
+    return np.ones((grid.nx, grid.ny+1), dtype=np.float64) * coeff
+
+def beta_plane_f_u(grid, f0, beta):
+    """Define a beta-plane approximation to the Coriolis force (u location)."""
+    _, Y = np.meshgrid(grid.xp1, grid.y)
+    fu = f0 + Y*beta
+    return fu
+
+def beta_plane_f_v(grid, f0, beta):
+    """Define a beta-plane approximation to the Coriolis force (v location)."""
+    _, Y = np.meshgrid(grid.x, grid.yp1)
+    fv = f0 + Y*beta
+    return fv
+
+def rectangular_pool(grid):
+    """The wet mask file for a maximal rectangular pool."""
+    nx = grid.nx; ny = grid.ny
+    wetmask = np.ones((nx, ny), dtype=np.float64)
+    wetmask[ 0, :] = 0
+    wetmask[-1, :] = 0
+    wetmask[ :, 0] = 0
+    wetmask[ :,-1] = 0
+    return wetmask
+
+specifier_rx = re.compile(r':(.*):(.*)')
+
+ok_generators = {
+    'depths': depths,
+    'beta_plane_f_u': beta_plane_f_u,
+    'beta_plane_f_v': beta_plane_f_v,
+    'f_plane_f_u': f_plane_f_u,
+    'f_plane_f_v': f_plane_f_v,
+    'rectangular_pool': rectangular_pool,
+    'wind_x': wind_x,
+    'wind_y': wind_y,
+}
+
+def interpret_data_specifier(string):
+    m = re.match(specifier_rx, string)
+    if m:
+        name = m.group(1)
+        arg_str = m.group(2)
+        if len(arg_str) > 0:
+            args = [float(a) for a in arg_str.split(',')]
+        else:
+            args = []
+        return (ok_generators[name], args)
+    else:
+        return None
+
+def interpret_requested_data(requested_data, shape, config):
+    """Interpret a flexible input data specification.
+
+    The requested_data can be one of
+
+    - TODO A string giving the path to a NetCDF file, whose content
+      will be interpolated to match the desired grid specification;
+
+    - A string giving the path to a raw Fortran array file, whose
+      content will be used as-is;
+
+    - TODO A numpy array in memory, whose content will be used as-is,
+      or TODO interpolated; or
+
+    - A string specifying auto-generation of the required data, in this format:
+      :<generator_func_name>:arg1,arg2,...argn
+
+    - Python objects specifying auto-generation of the required data.
+      In this case, `interpret_requested_data` will construct the
+      appropriate `Grid` instance and pass it, together with the
+      `requested_data`, to an appropriate meta-generator for the array
+      shape of the needful datum (determined by the `shape` argument).
+      The exact API varies with the meta-generator, but they typically
+      interpret numbers as that constant and functions as an analytic
+      definition of the field, which is evaluated on appropriate numpy
+      arrays to produce the needed numerical values.
+    """
+    grid = Grid(config.getint("grid", "nx"), config.getint("grid", "ny"),
+                config.getfloat("grid", "dx"), config.getfloat("grid", "dy"))
+    if isinstance(requested_data, basestring):
+        candidate = interpret_data_specifier(requested_data)
+        if candidate is not None:
+            (func, args) = candidate
+            return func(grid, *args)
+        else:
+            # Assume Fortran file name
+            with fortran_file(requested_data, 'r') as f:
+                return f.read_reals(dtype=np.float64)
+    else:
+        if shape == "3d":
+            return depths(grid, *requested_data)
+        if shape == "2dx":
+            return wind_x(grid, requested_data)
+        else:
+            raise Exception("TODO implement custom generation for other input shapes")
diff --git a/aronnax/driver.py b/aronnax/driver.py
new file mode 100644
index 00000000..2353e0e8
--- /dev/null
+++ b/aronnax/driver.py
@@ -0,0 +1,257 @@
+import ConfigParser as par
+import os
+import os.path as p
+import subprocess as sub
+import time
+
+from aronnax.core import fortran_file
+from aronnax.core import interpret_requested_data
+from aronnax.utils import working_directory
+
+self_path = p.dirname(p.abspath(__file__))
+root_path = p.dirname(self_path)
+
+def simulate(work_dir=".", config_path="aronnax.conf", **options):
+    """Main entry point for running an Aronnax simulation.
+
+    A simulation occurs in the working directory given by the
+    `work_dir` parameter, which defaults to the current directory when
+    `simulate` is invoked.  The default arrangement of the working
+    directory is as follows:
+
+        - aronnax.conf - configuration file for that run
+        - aronnax-merged.conf - file to save effective configuration, including
+          effects of options passed to `simulate`. This file is automatically
+          generated by merging the aronnax.conf file with the options passed to
+          this function
+        - parameters.in - relevant portions of aronnax-merged.conf in Fortran
+          namelist format. Also generated automatically
+        - input/ - subdirectory where Aronnax will save input field files
+          in Fortran raw array format
+        - output/ - subdirectory where Aronnax will save output field files
+          in Fortran raw array format
+
+    The process for a simulation is to
+    
+        1. Compute the configuration
+        2. Recompile the Fortran core if necessary
+        3. Save the computed configuration in aronnax-merged.conf
+        4. Write raw-format input fields into input/
+        5. Write parameters.in
+        6. Execute the Fortran core, which writes progress messages
+           to standard output and raw-format output fields into output/
+
+    All the simulation parameters can be controlled from the
+    configuration file aronnax.conf, and additionally can be
+    overridden by passing them as optional arguments to `simulate`.
+
+    Calling `simulate` directly provides one capability that cannot be
+    accessed from the configuration file: custom idealized input
+    generators.
+
+    """
+    config_file = p.join(work_dir, config_path)
+    config = default_configuration()
+    config.read(config_file)
+    merge_config(config, options)
+    with working_directory(work_dir):
+        compile_core(config)
+        # XXX Try to avoid overwriting the input configuration
+        with open('aronnax-merged.conf', 'w') as f:
+            config.write(f)
+        sub.check_call(["rm", "-rf", "input/"])
+        sub.check_call(["rm", "-rf", "output/"])
+        sub.check_call(["mkdir", "-p", "output/"])
+        with working_directory("input"):
+            generate_input_data_files(config)
+        generate_parameters_file(config)
+        then = time.time()
+        run_executable(config)
+        core_run_time = time.time() - then
+        sub.check_call(["rm", "-rf", "netcdf-output/"])
+        sub.check_call(["mkdir", "-p", "netcdf-output/"])
+        convert_output_to_netcdf(config)
+        return core_run_time
+
+def default_configuration():
+    """Configuration defaults before parsing aronnax.conf.
+
+    The configuration is represented as a ConfigParser.RawConfigParser instance."""
+    config = par.RawConfigParser()
+    for section in sections:
+        config.add_section(section)
+    config.set("executable", "valgrind", "False")
+    config.set("executable", "perf", "False")
+    config.optionxform = str
+    return config
+
+sections = ["executable", "numerics", "model", "sponge",
+            "physics", "grid", "initial_conditions", "external_forcing"]
+
+section_map = {
+    "au"                   : "numerics",
+    "ah"                   : "numerics",
+    "ar"                   : "numerics",
+    "botDrag"              : "numerics",
+    "dt"                   : "numerics",
+    "slip"                 : "numerics",
+    "nTimeSteps"           : "numerics",
+    "dumpFreq"             : "numerics",
+    "avFreq"               : "numerics",
+    "hmin"                 : "numerics",
+    "maxits"               : "numerics",
+    "eps"                  : "numerics",
+    "freesurfFac"          : "numerics",
+    "thickness_error"      : "numerics",
+    "hmean"                : "model",
+    "depthFile"            : "model",
+    "H0"                   : "model",
+    "RedGrav"              : "model",
+    "spongeHTimeScaleFile" : "sponge",
+    "spongeUTimeScaleFile" : "sponge",
+    "spongeVTimeScaleFile" : "sponge",
+    "spongeHFile"          : "sponge",
+    "spongeUFile"          : "sponge",
+    "spongeVFile"          : "sponge",
+    "g_vec"                : "physics",
+    "rho0"                 : "physics",
+    "nx"                   : "grid",
+    "ny"                   : "grid",
+    "layers"               : "grid",
+    "dx"                   : "grid",
+    "dy"                   : "grid",
+    "fUfile"               : "grid",
+    "fVfile"               : "grid",
+    "wetMaskFile"          : "grid",
+    "initUfile"            : "initial_conditions",
+    "initVfile"            : "initial_conditions",
+    "initHfile"            : "initial_conditions",
+    "initEtaFile"          : "initial_conditions",
+    "zonalWindFile"        : "external_forcing",
+    "meridionalWindFile"   : "external_forcing",
+    "DumpWind"             : "external_forcing",
+    "wind_mag_time_series_file" : "external_forcing",
+    "exe"                  : "executable",
+    "valgrind"             : "executable",
+    "perf"                 : "executable",
+}
+
+def merge_config(config, options):
+    """Merge the options given in the `options` dict into the RawConfigParser instance `config`.
+
+    Mutates the given config instance."""
+    for k, v in options.iteritems():
+        if k in section_map:
+            section = section_map[k]
+            if not config.has_section(section):
+                config.add_section(section)
+            if v == True: v = "yes"
+            if v == False: v = "no"
+            config.set(section, k, v)
+        else:
+            raise Exception("Unrecognized option", k)
+
+def compile_core(config):
+    """Compile the Aronnax core, if needed."""
+    core_name = config.get("executable", "exe")
+    with working_directory(root_path):
+        sub.check_call(["make", core_name])
+
+data_types = {
+    "depthFile"            : "2d",
+    "spongeHTimeScaleFile" : "3d",
+    "spongeUTimeScaleFile" : "3d",
+    "spongeVTimeScaleFile" : "3d",
+    "spongeHFile"          : "3d",
+    "spongeUFile"          : "3d",
+    "spongeVFile"          : "3d",
+    "fUfile"               : "2dx",
+    "fVfile"               : "2dy",
+    "wetMaskFile"          : "2d",
+    "initUfile"            : "3d",
+    "initVfile"            : "3d",
+    "initHfile"            : "3d",
+    "initEtaFile"          : "2d",
+    "zonalWindFile"        : "2dx",
+    "meridionalWindFile"   : "2dy",
+    "wind_mag_time_series_file" : "time",
+}
+
+def is_file_name_option(name):
+    return name.endswith("File") or name.endswith("file")
+
+def generate_input_data_files(config):
+    for name, section in section_map.iteritems():
+        if not is_file_name_option(name):
+            continue
+        if not config.has_option(section, name):
+            continue
+        requested_data = config.get(section, name)
+        generated_data = interpret_requested_data(
+            requested_data, data_types[name], config)
+        if generated_data is not None:
+            with fortran_file(name + '.bin', 'w') as f:
+                f.write_record(generated_data)
+
+def fortran_option_string(section, name, config):
+    """Convert option values to strings that Fortran namelists will understand correctly.
+
+    Two conversions are of interest: Booleans are rendered as
+    .TRUE. or .FALSE., and options that are input fields are rendered
+    as the file names where `generate_input_data_files` has written
+    those data, or the Fortran empty string `''` if no file is written
+    (which means the core should use its internal default).
+    """
+    if is_file_name_option(name):
+        if config.has_option(section, name):
+            # A file was generated
+            return "'%s'" % (p.join("input", name + '.bin'),)
+        else:
+            return "''"
+    if name in ["RedGrav", "DumpWind"]:
+        if config.getboolean(section, name):
+            return ".TRUE."
+        else:
+            return ".FALSE."
+    else:
+        if config.has_option(section, name):
+            return config.get(section, name)
+        else:
+            return None
+
+def generate_parameters_file(config):
+    with open('parameters.in', 'w') as f:
+        for section in config.sections():
+            f.write(' &')
+            f.write(section.upper())
+            f.write('\n')
+            for (name, section1) in section_map.iteritems():
+                if section1 != section: continue
+                val = fortran_option_string(section, name, config)
+                if val is not None:
+                    f.write(' %s = %s,\n' % (name, val))
+            f.write(' /\n')
+
+def run_executable(config):
+    """Run the compiled Fortran core, possibly in a test or debug regime."""
+    core_name = config.get("executable", "exe")
+    env = dict(os.environ, GFORTRAN_STDERR_UNIT="17")
+    if config.getboolean("executable", "valgrind") \
+       or 'ARONNAX_TEST_VALGRIND_ALL' in os.environ:
+        assert not config.getboolean("executable", "perf")
+        sub.check_call(["valgrind", "--error-exitcode=5", p.join(root_path, core_name)],
+            env=env)
+    elif config.getboolean("executable", "perf"):
+        perf_cmds = ["perf", "stat", "-e", "r530010", # "flops", on my CPU.
+            "-e", "L1-dcache-loads", "-e", "L1-dcache-load-misses",
+            "-e", "L1-dcache-stores", "-e", "L1-dcache-store-misses",
+            "-e", "L1-icache-loads", "-e", "L1-icache-misses",
+            "-e", "L1-dcache-prefetches",
+            "-e", "branch-instructions", "-e", "branch-misses"]
+        sub.check_call(perf_cmds + [p.join(root_path, core_name)], env=env)
+    else:
+        sub.check_call([p.join(root_path, core_name)], env=env)
+
+def convert_output_to_netcdf(config):
+    # TODO Issue #30
+    pass
diff --git a/aronnax/utils.py b/aronnax/utils.py
new file mode 100644
index 00000000..ec3c78df
--- /dev/null
+++ b/aronnax/utils.py
@@ -0,0 +1,14 @@
+from contextlib import contextmanager
+import os
+import subprocess as sub
+
+@contextmanager
+def working_directory(path):
+    old_path = os.getcwd()
+    sub.check_call(["mkdir", "-p", path])
+    os.chdir(path)
+    try:
+        yield
+    finally:
+        os.chdir(old_path)
+
diff --git a/benchmarks/.gitignore b/benchmarks/.gitignore
index 7dce5900..9b68e0e0 100644
--- a/benchmarks/.gitignore
+++ b/benchmarks/.gitignore
@@ -4,3 +4,5 @@ run_finished.txt
 layer thickness dropped below hmin.txt
 *scaling.png
 **/times.pkl
+parameters.in
+aronnax-merged.conf
diff --git a/benchmarks/benchmark.py b/benchmarks/benchmark.py
index ceea7aef..1aadb6e1 100644
--- a/benchmarks/benchmark.py
+++ b/benchmarks/benchmark.py
@@ -7,36 +7,37 @@
 matplotlib.use("Agg")
 import matplotlib.pyplot as plt
 
+from aronnax.utils import working_directory
+import aronnax.driver as aro
+
 self_path = p.dirname(p.abspath(__file__))
 root_path = p.dirname(self_path)
 
-import sys
-sys.path.append(p.join(root_path, 'test'))
-import output_preservation_test as opt
-
 n_time_steps = 502.0
 scale_factor = 1000 / n_time_steps # Show times in ms
 
 def benchmark_gaussian_bump_red_grav_save(grid_points):
     run_time_O1 = np.zeros(len(grid_points))
     run_time_Ofast = np.zeros(len(grid_points))
+    def bump(X, Y):
+        return 500. + 20*np.exp(-((6e5-X)**2 + (5e5-Y)**2)/(2*1e5**2))
 
-    with opt.working_directory(p.join(self_path, "beta_plane_bump_red_grav")):
+    with working_directory(p.join(self_path, "beta_plane_bump_red_grav")):
         aro_exec = "aronnax_test"
         for counter, nx in enumerate(grid_points):
-            run_time_O1[counter] = opt.run_experiment(
-                opt.write_input_beta_plane_bump_red_grav, nx, nx, 1, aro_exec)
+            run_time_O1[counter] = aro.simulate(
+                exe=aro_exec, initHfile=[bump], nx=nx, ny=nx)
 
         aro_exec = "aronnax_core"
         for counter, nx in enumerate(grid_points):
-            run_time_Ofast[counter] = opt.run_experiment(
-                opt.write_input_beta_plane_bump_red_grav, nx, nx, 1, aro_exec)
+            run_time_Ofast[counter] = aro.simulate(
+                exe=aro_exec, initHfile=[bump], nx=nx, ny=nx)
 
         with open("times.pkl", "w") as f:
             pkl.dump((grid_points, run_time_O1, run_time_Ofast), f)
 
 def benchmark_gaussian_bump_red_grav_plot():
-    with opt.working_directory(p.join(self_path, "beta_plane_bump_red_grav")):
+    with working_directory(p.join(self_path, "beta_plane_bump_red_grav")):
         with open("times.pkl", "r") as f:
             (grid_points, run_time_O1, run_time_Ofast) = pkl.load(f)
 
@@ -64,21 +65,23 @@ def benchmark_gaussian_bump_red_grav(grid_points):
 def benchmark_gaussian_bump_save(grid_points):
     run_time_O1 = np.zeros(len(grid_points))
     run_time_Ofast = np.zeros(len(grid_points))
+    def bump(X, Y):
+        return 500. + 20*np.exp(-((6e5-X)**2 + (5e5-Y)**2)/(2*1e5**2))
 
-    with opt.working_directory(p.join(self_path, "beta_plane_bump")):
+    with working_directory(p.join(self_path, "beta_plane_bump")):
         aro_exec = "aronnax_test"
         for counter, nx in enumerate(grid_points):
-            run_time_O1[counter] = opt.run_experiment(
-                  opt.write_input_beta_plane_bump, nx, nx, 2, aro_exec)
+            run_time_O1[counter] = aro.simulate(
+                exe=aro_exec, initHfile=[bump, lambda X, Y: 2000. - bump(X, Y)], nx=nx, ny=nx)
         aro_exec = "aronnax_core"
         for counter, nx in enumerate(grid_points):
-            run_time_Ofast[counter] = opt.run_experiment(
-                  opt.write_input_beta_plane_bump, nx, nx, 2, aro_exec)
+            run_time_Ofast[counter] = aro.simulate(
+                exe=aro_exec, initHfile=[bump, lambda X, Y: 2000. - bump(X, Y)], nx=nx, ny=nx)
         with open("times.pkl", "w") as f:
             pkl.dump((grid_points, run_time_O1, run_time_Ofast), f)
 
 def benchmark_gaussian_bump_plot():
-    with opt.working_directory(p.join(self_path, "beta_plane_bump")):
+    with working_directory(p.join(self_path, "beta_plane_bump")):
         with open("times.pkl", "r") as f:
             (grid_points, run_time_O1, run_time_Ofast) = pkl.load(f)
 
diff --git a/benchmarks/beta_plane_bump/aronnax.conf b/benchmarks/beta_plane_bump/aronnax.conf
new file mode 100644
index 00000000..22bb90d8
--- /dev/null
+++ b/benchmarks/beta_plane_bump/aronnax.conf
@@ -0,0 +1,54 @@
+# Aronnax configuration file. Change the values, but not the names.
+# 
+# au is viscosity
+# ah is thickness diffusivity
+# ar is linear drag between layers
+# dt is time step
+# slip is free-slip (=0), no-slip (=1), or partial slip (something in between)
+# nTimeSteps: number of timesteps before stopping
+# dumpFreq: frequency of snapshot output
+# avFreq: frequency of averaged output
+# hmin: minimum layer thickness allowed by model (for stability)
+# maxits: maximum iterations for the successive over relaxation algorithm. Should be at least max(nx,ny), and probably nx*ny
+# eps: convergence tolerance for SOR solver
+# freesurfFac: 1. = linear implicit free surface, 0. = rigid lid. So far all tests using freesurfFac = 1. have failed 
+# g is the gravity at interfaces (including surface). must have as many entries as there are layers
+# input files are where to look for the various inputs
+
+[numerics]
+au = 500.
+ah = 0.0
+ar = 0.0
+botDrag = 1e-6
+dt = 600.
+slip = 0.0
+nTimeSteps = 502
+dumpFreq = 1.2e5
+avFreq = 1.2e5
+hmin = 100
+maxits = 1000
+eps = 1e-2
+freesurfFac = 0.
+thickness_error = 1e-2
+
+[model]
+hmean = 400.,1600.
+H0 = 2000.
+RedGrav = no
+
+[physics]
+g_vec = 9.8, 0.01
+rho0 = 1035.
+
+[grid]
+nx = 10
+ny = 10
+layers = 2
+dx = 2e4
+dy = 2e4
+fUfile = :beta_plane_f_u:1e-5,2e-11
+fVfile = :beta_plane_f_v:1e-5,2e-11
+wetMaskFile = :rectangular_pool:
+
+[external_forcing]
+DumpWind = no
diff --git a/benchmarks/beta_plane_bump/parameters.in b/benchmarks/beta_plane_bump/parameters.in
deleted file mode 100644
index 45b75577..00000000
--- a/benchmarks/beta_plane_bump/parameters.in
+++ /dev/null
@@ -1,80 +0,0 @@
-! Parameter file. Change the values, but not the names.
-! 
-! Each namelist ends with a /
-! This is important and should not be deleted.
-!
-! au is viscosity
-! ah is thickness diffusivity
-! ar is linear drag between layers
-! dt is time step
-! slip is free-slip (=0), no-slip (=1), or partial slip (something in between)
-! nTimeSteps: number of timesteps before stopping
-! dumpFreq: frequency of snapshot output
-! avFreq: frequency of averaged output
-! hmin: minimum layer thickness allowed by model (for stability)
-! maxits: maximum iterations for the successive over relaxation algorithm. Should be at least max(nx,ny), and probably nx*ny
-! eps: convergence tolerance for SOR solver
-! freesurfFac: 1. = linear implicit free surface, 0. = rigid lid. So far all tests using freesurfFac = 1. have failed 
-! g is the gravity at interfaces (including surface). must have as many entries as there are layers
-! input files are where to look for the various inputs
-!
- &NUMERICS
- au = 500.,
- ah = 0.0,
- ar = 0.0,
- botDrag = 1e-6,
- dt = 600.,
- slip = 0.0,
- nTimeSteps = 502,
- dumpFreq = 1.2e5,
- avFreq = 1.2e5,
- hmin = 100,
- maxits = 1000,
- eps = 1e-2,
- freesurfFac = 0.,
- thickness_error = 1e-2,
- /
- &MODEL
- hmean = 400.,1600.,
- depthFile = '',
- H0 = 2000.,
- RedGrav = .FALSE.,
- /
- &SPONGE
- spongeHTimeScaleFile = '',
- spongeUTimeScaleFile = '',
- spongeVTimeScaleFile = '',
- spongeHFile = '',
- spongeUfile = '',
- spongeVfile = '',
- /
- &PHYSICS
- g_vec = 9.8, 0.01,
- rho0 = 1035.,
- /
- &GRID
- nx = 100,
- ny = 100,
- layers = 2,
- dx = 2e4,
- dy = 2e4,
- fUfile = 'input/fu.bin',
- fVfile = 'input/fv.bin',
- wetMaskFile = 'input/wetmask.bin',
- /
-! Inital conditions for u, v, and h
- &INITIAL_CONDITONS
- initUfile = '',
- initVfile = '',
- initHfile = 'input/initH.bin',
- initEtaFile = '',
- /
- &EXTERNAL_FORCING
- zonalWindFile = '',
- meridionalWindFile = '', 
- DumpWind = .FALSE.,
- wind_mag_time_series_file = '',
- /
-! input/initH.bin
-! 
- &
diff --git a/benchmarks/beta_plane_bump_red_grav/aronnax.conf b/benchmarks/beta_plane_bump_red_grav/aronnax.conf
new file mode 100644
index 00000000..476a767b
--- /dev/null
+++ b/benchmarks/beta_plane_bump_red_grav/aronnax.conf
@@ -0,0 +1,53 @@
+# Aronnax configuration file. Change the values, but not the names.
+# 
+# au is viscosity
+# ah is thickness diffusivity
+# ar is linear drag between layers
+# dt is time step
+# slip is free-slip (=0), no-slip (=1), or partial slip (something in between)
+# nTimeSteps: number of timesteps before stopping
+# dumpFreq: frequency of snapshot output
+# avFreq: frequency of averaged output
+# hmin: minimum layer thickness allowed by model (for stability)
+# maxits: maximum iterations for the successive over relaxation algorithm. Should be at least max(nx,ny), and probably nx*ny
+# eps: convergence tolerance for SOR solver
+# freesurfFac: 1. = linear implicit free surface, 0. = rigid lid. So far all tests using freesurfFac = 1. have failed 
+# g is the gravity at interfaces (including surface). must have as many entries as there are layers
+# input files are where to look for the various inputs
+
+[numerics]
+au = 500.
+ah = 0.0
+ar = 1e-8
+dt = 600.
+slip = 0.0
+nTimeSteps = 502
+dumpFreq = 6e4
+avFreq = 3e5
+hmin = 100
+maxits = 1000
+eps = 1e-2
+freesurfFac = 0.
+thickness_error = 1e-2
+
+[model]
+hmean = 400.
+H0 = 2000.
+RedGrav = yes
+
+[physics]
+g_vec = 0.01
+rho0 = 1035.
+
+[grid]
+nx = 500
+ny = 500
+layers = 1
+dx = 2e4
+dy = 2e4
+fUfile = :beta_plane_f_u:1e-5,2e-11
+fVfile = :beta_plane_f_v:1e-5,2e-11
+wetMaskFile = :rectangular_pool:
+
+[external_forcing]
+DumpWind = no
diff --git a/benchmarks/beta_plane_bump_red_grav/parameters.in b/benchmarks/beta_plane_bump_red_grav/parameters.in
deleted file mode 100644
index 6791797f..00000000
--- a/benchmarks/beta_plane_bump_red_grav/parameters.in
+++ /dev/null
@@ -1,79 +0,0 @@
-! Parameter file. Change the values, but not the names.
-! 
-! Each namelist ends with a /
-! This is important and should not be deleted.
-!
-! au is viscosity
-! ah is thickness diffusivity
-! ar is linear drag between layers
-! dt is time step
-! slip is free-slip (=0), no-slip (=1), or partial slip (something in between)
-! nTimeSteps: number of timesteps before stopping
-! dumpFreq: frequency of snapshot output
-! avFreq: frequency of averaged output
-! hmin: minimum layer thickness allowed by model (for stability)
-! maxits: maximum iterations for the successive over relaxation algorithm. Should be at least max(nx,ny), and probably nx*ny
-! eps: convergence tolerance for SOR solver
-! freesurfFac: 1. = linear implicit free surface, 0. = rigid lid. So far all tests using freesurfFac = 1. have failed 
-! g is the gravity at interfaces (including surface). must have as many entries as there are layers
-! input files are where to look for the various inputs
-!
- &NUMERICS
- au = 500.,
- ah = 0.0,
- ar = 1e-8,
- dt = 600.,
- slip = 0.0,
- nTimeSteps = 502,
- dumpFreq = 6e4,
- avFreq = 3e5,
- hmin = 100,
- maxits = 1000,
- eps = 1e-2,
- freesurfFac = 0.,
- thickness_error = 1e-2,
- /
- &MODEL
- hmean = 400.,
- depthFile = '',
- H0 = 2000.,
- RedGrav = .TRUE.,
- /
- &SPONGE
- spongeHTimeScaleFile = '',
- spongeUTimeScaleFile = '',
- spongeVTimeScaleFile = '',
- spongeHFile = '',
- spongeUfile = '',
- spongeVfile = '',
- /
- &PHYSICS
- g_vec = 0.01,
- rho0 = 1035.,
- /
- &GRID
- nx = 500,
- ny = 500,
- layers = 1,
- dx = 2e4,
- dy = 2e4,
- fUfile = 'input/fu.bin',
- fVfile = 'input/fv.bin',
- wetMaskFile = 'input/wetmask.bin',
- /
-! Inital conditions for u, v, and h
- &INITIAL_CONDITONS
- initUfile = '',
- initVfile = '',
- initHfile = 'input/initH.bin',
- initEtaFile = '',
- /
- &EXTERNAL_FORCING
- zonalWindFile = '',
- meridionalWindFile = '', 
- DumpWind = .FALSE.,
- wind_mag_time_series_file = '',
- /
-! input/initH.bin
-! 
- &
diff --git a/benchmarks/profile.py b/benchmarks/profile.py
index f6d3808d..f3872d01 100644
--- a/benchmarks/profile.py
+++ b/benchmarks/profile.py
@@ -1,21 +1,24 @@
 import os.path as p
+import sys
+
+import numpy as np
+
+from aronnax.utils import working_directory
+import aronnax.driver as aro
 
 self_path = p.dirname(p.abspath(__file__))
-root_path = p.dirname(self_path)
 
-import sys
-sys.path.append(p.join(root_path, 'test'))
-import output_preservation_test as opt
+def bump(X, Y):
+    return 500. + 20*np.exp(-((6e5-X)**2 + (5e5-Y)**2)/(2*1e5**2))
 
 def do_red_grav(nx, aro_build, perf):
-    with opt.working_directory(p.join(self_path, "beta_plane_bump_red_grav")):
-        opt.run_experiment(
-            opt.write_input_beta_plane_bump_red_grav, nx, nx, 1, aro_build, perf=perf)
+    with working_directory(p.join(self_path, "beta_plane_bump_red_grav")):
+        aro.simulate(exe=aro_build, initHfile=[bump], nx=nx, ny=nx, perf=perf)
 
 def do_n_layer(nx, aro_build, perf):
-    with opt.working_directory(p.join(self_path, "beta_plane_bump")):
-        opt.run_experiment(
-            opt.write_input_beta_plane_bump, nx, nx, 2, aro_build, perf=perf)
+    with working_directory(p.join(self_path, "beta_plane_bump")):
+        aro.simulate(exe=aro_build, nx=nx, ny=nx, perf=perf,
+                     initHfile=[bump, lambda X, Y: 2000. - bump(X, Y)])
 
 def main():
     aro_build = "aronnax_prof"
diff --git a/docs/index.rst b/docs/index.rst
index 243ae5a8..405b71ad 100644
--- a/docs/index.rst
+++ b/docs/index.rst
@@ -44,6 +44,7 @@ Aronnax is
    examples
    aronnax_model
    input_generation
+   running_aronnax
    output_helpers
    verification
    benchmarks
diff --git a/docs/input_generation.rst b/docs/input_generation.rst
index 4d67d404..f2ff14ed 100644
--- a/docs/input_generation.rst
+++ b/docs/input_generation.rst
@@ -1,7 +1,9 @@
 Input generators
 *******************
 
-Aronnax includes a number helper functions for generating input fields.
+Aronnax includes a number of helper functions for generating input fields.
+These can be called with numerical arguments (generating very simple input fields)
+directly from the aronnax.conf file, using the syntax `:generator_name:arg1,arg2,...,argn`.
 
 Grid
 =====
@@ -14,22 +16,30 @@ Forcings
 ==========
 
 
-.. autofunction:: aronnax.write_wind_x
+.. autofunction:: aronnax.wind_x
 
-.. autofunction:: aronnax.write_wind_y
+.. autofunction:: aronnax.wind_y
 
 
 Initial conditions
 ===================
 
-.. autofunction:: aronnax.write_initial_heights
+.. autofunction:: aronnax.depths
 
 
-Domain fields
+Coriolis fields
 ===============
 
-.. autofunction:: aronnax.write_f_plane
+.. autofunction:: aronnax.f_plane_f_u
 
-.. autofunction:: aronnax.write_beta_plane
+.. autofunction:: aronnax.f_plane_f_v
 
-.. autofunction:: aronnax.write_rectangular_pool
\ No newline at end of file
+.. autofunction:: aronnax.beta_plane_f_u
+
+.. autofunction:: aronnax.beta_plane_f_v
+
+
+Domain shape
+============
+
+.. autofunction:: aronnax.rectangular_pool
diff --git a/docs/running_aronnax.rst b/docs/running_aronnax.rst
new file mode 100644
index 00000000..7fddf217
--- /dev/null
+++ b/docs/running_aronnax.rst
@@ -0,0 +1,5 @@
+Running Aronnax
+*******************
+
+
+.. autofunction:: aronnax.driver.simulate
diff --git a/test/.gitignore b/test/.gitignore
index a8617d13..abd91f7e 100644
--- a/test/.gitignore
+++ b/test/.gitignore
@@ -2,3 +2,5 @@ input/
 output/
 run_finished.txt
 layer thickness dropped below hmin.txt
+parameters.in
+aronnax-merged.conf
diff --git a/test/beta_plane_bump/aronnax.conf b/test/beta_plane_bump/aronnax.conf
new file mode 100644
index 00000000..22bb90d8
--- /dev/null
+++ b/test/beta_plane_bump/aronnax.conf
@@ -0,0 +1,54 @@
+# Aronnax configuration file. Change the values, but not the names.
+# 
+# au is viscosity
+# ah is thickness diffusivity
+# ar is linear drag between layers
+# dt is time step
+# slip is free-slip (=0), no-slip (=1), or partial slip (something in between)
+# nTimeSteps: number of timesteps before stopping
+# dumpFreq: frequency of snapshot output
+# avFreq: frequency of averaged output
+# hmin: minimum layer thickness allowed by model (for stability)
+# maxits: maximum iterations for the successive over relaxation algorithm. Should be at least max(nx,ny), and probably nx*ny
+# eps: convergence tolerance for SOR solver
+# freesurfFac: 1. = linear implicit free surface, 0. = rigid lid. So far all tests using freesurfFac = 1. have failed 
+# g is the gravity at interfaces (including surface). must have as many entries as there are layers
+# input files are where to look for the various inputs
+
+[numerics]
+au = 500.
+ah = 0.0
+ar = 0.0
+botDrag = 1e-6
+dt = 600.
+slip = 0.0
+nTimeSteps = 502
+dumpFreq = 1.2e5
+avFreq = 1.2e5
+hmin = 100
+maxits = 1000
+eps = 1e-2
+freesurfFac = 0.
+thickness_error = 1e-2
+
+[model]
+hmean = 400.,1600.
+H0 = 2000.
+RedGrav = no
+
+[physics]
+g_vec = 9.8, 0.01
+rho0 = 1035.
+
+[grid]
+nx = 10
+ny = 10
+layers = 2
+dx = 2e4
+dy = 2e4
+fUfile = :beta_plane_f_u:1e-5,2e-11
+fVfile = :beta_plane_f_v:1e-5,2e-11
+wetMaskFile = :rectangular_pool:
+
+[external_forcing]
+DumpWind = no
diff --git a/test/beta_plane_bump/good-output/av.eta.0000000001 b/test/beta_plane_bump/good-output/av.eta.0000000001
index 3ddfc970..7add0c6b 100644
Binary files a/test/beta_plane_bump/good-output/av.eta.0000000001 and b/test/beta_plane_bump/good-output/av.eta.0000000001 differ
diff --git a/test/beta_plane_bump/good-output/av.eta.0000000201 b/test/beta_plane_bump/good-output/av.eta.0000000201
index fab17472..4b99c0aa 100644
Binary files a/test/beta_plane_bump/good-output/av.eta.0000000201 and b/test/beta_plane_bump/good-output/av.eta.0000000201 differ
diff --git a/test/beta_plane_bump/good-output/av.eta.0000000401 b/test/beta_plane_bump/good-output/av.eta.0000000401
index 0362fc7f..8a5462ee 100644
Binary files a/test/beta_plane_bump/good-output/av.eta.0000000401 and b/test/beta_plane_bump/good-output/av.eta.0000000401 differ
diff --git a/test/beta_plane_bump/good-output/av.h.0000000001 b/test/beta_plane_bump/good-output/av.h.0000000001
index 1c9ab437..96ab4617 100644
Binary files a/test/beta_plane_bump/good-output/av.h.0000000001 and b/test/beta_plane_bump/good-output/av.h.0000000001 differ
diff --git a/test/beta_plane_bump/good-output/av.h.0000000201 b/test/beta_plane_bump/good-output/av.h.0000000201
index 13879006..93fd8e4a 100644
Binary files a/test/beta_plane_bump/good-output/av.h.0000000201 and b/test/beta_plane_bump/good-output/av.h.0000000201 differ
diff --git a/test/beta_plane_bump/good-output/av.h.0000000401 b/test/beta_plane_bump/good-output/av.h.0000000401
index c678e0d1..07413229 100644
Binary files a/test/beta_plane_bump/good-output/av.h.0000000401 and b/test/beta_plane_bump/good-output/av.h.0000000401 differ
diff --git a/test/beta_plane_bump/good-output/av.u.0000000001 b/test/beta_plane_bump/good-output/av.u.0000000001
index e2027ed8..7d9a4732 100644
Binary files a/test/beta_plane_bump/good-output/av.u.0000000001 and b/test/beta_plane_bump/good-output/av.u.0000000001 differ
diff --git a/test/beta_plane_bump/good-output/av.u.0000000201 b/test/beta_plane_bump/good-output/av.u.0000000201
index 55632b7f..cd413b97 100644
Binary files a/test/beta_plane_bump/good-output/av.u.0000000201 and b/test/beta_plane_bump/good-output/av.u.0000000201 differ
diff --git a/test/beta_plane_bump/good-output/av.u.0000000401 b/test/beta_plane_bump/good-output/av.u.0000000401
index 9cfc75a8..88d83b92 100644
Binary files a/test/beta_plane_bump/good-output/av.u.0000000401 and b/test/beta_plane_bump/good-output/av.u.0000000401 differ
diff --git a/test/beta_plane_bump/good-output/av.v.0000000001 b/test/beta_plane_bump/good-output/av.v.0000000001
index d16bf26b..fe4ca24c 100644
Binary files a/test/beta_plane_bump/good-output/av.v.0000000001 and b/test/beta_plane_bump/good-output/av.v.0000000001 differ
diff --git a/test/beta_plane_bump/good-output/av.v.0000000201 b/test/beta_plane_bump/good-output/av.v.0000000201
index 5a2c9446..fc935dab 100644
Binary files a/test/beta_plane_bump/good-output/av.v.0000000201 and b/test/beta_plane_bump/good-output/av.v.0000000201 differ
diff --git a/test/beta_plane_bump/good-output/av.v.0000000401 b/test/beta_plane_bump/good-output/av.v.0000000401
index c3155807..a1416ca0 100644
Binary files a/test/beta_plane_bump/good-output/av.v.0000000401 and b/test/beta_plane_bump/good-output/av.v.0000000401 differ
diff --git a/test/beta_plane_bump/good-output/snap.eta.0000000001 b/test/beta_plane_bump/good-output/snap.eta.0000000001
index b9545983..b6baa0a4 100644
Binary files a/test/beta_plane_bump/good-output/snap.eta.0000000001 and b/test/beta_plane_bump/good-output/snap.eta.0000000001 differ
diff --git a/test/beta_plane_bump/good-output/snap.eta.0000000201 b/test/beta_plane_bump/good-output/snap.eta.0000000201
index c7dc102f..64fe68ce 100644
Binary files a/test/beta_plane_bump/good-output/snap.eta.0000000201 and b/test/beta_plane_bump/good-output/snap.eta.0000000201 differ
diff --git a/test/beta_plane_bump/good-output/snap.eta.0000000401 b/test/beta_plane_bump/good-output/snap.eta.0000000401
index f9db8cb0..d727d9ec 100644
Binary files a/test/beta_plane_bump/good-output/snap.eta.0000000401 and b/test/beta_plane_bump/good-output/snap.eta.0000000401 differ
diff --git a/test/beta_plane_bump/good-output/snap.h.0000000001 b/test/beta_plane_bump/good-output/snap.h.0000000001
index 77d00e23..5ae422cf 100644
Binary files a/test/beta_plane_bump/good-output/snap.h.0000000001 and b/test/beta_plane_bump/good-output/snap.h.0000000001 differ
diff --git a/test/beta_plane_bump/good-output/snap.h.0000000201 b/test/beta_plane_bump/good-output/snap.h.0000000201
index 19b35322..1a744717 100644
Binary files a/test/beta_plane_bump/good-output/snap.h.0000000201 and b/test/beta_plane_bump/good-output/snap.h.0000000201 differ
diff --git a/test/beta_plane_bump/good-output/snap.h.0000000401 b/test/beta_plane_bump/good-output/snap.h.0000000401
index d61bbfc8..a71bddd6 100644
Binary files a/test/beta_plane_bump/good-output/snap.h.0000000401 and b/test/beta_plane_bump/good-output/snap.h.0000000401 differ
diff --git a/test/beta_plane_bump/good-output/snap.u.0000000001 b/test/beta_plane_bump/good-output/snap.u.0000000001
index e1ca9e51..07114b13 100644
Binary files a/test/beta_plane_bump/good-output/snap.u.0000000001 and b/test/beta_plane_bump/good-output/snap.u.0000000001 differ
diff --git a/test/beta_plane_bump/good-output/snap.u.0000000201 b/test/beta_plane_bump/good-output/snap.u.0000000201
index ceb90433..402516a4 100644
Binary files a/test/beta_plane_bump/good-output/snap.u.0000000201 and b/test/beta_plane_bump/good-output/snap.u.0000000201 differ
diff --git a/test/beta_plane_bump/good-output/snap.u.0000000401 b/test/beta_plane_bump/good-output/snap.u.0000000401
index 538d9592..f95455c0 100644
Binary files a/test/beta_plane_bump/good-output/snap.u.0000000401 and b/test/beta_plane_bump/good-output/snap.u.0000000401 differ
diff --git a/test/beta_plane_bump/good-output/snap.v.0000000001 b/test/beta_plane_bump/good-output/snap.v.0000000001
index 5e23bd1b..77710e09 100644
Binary files a/test/beta_plane_bump/good-output/snap.v.0000000001 and b/test/beta_plane_bump/good-output/snap.v.0000000001 differ
diff --git a/test/beta_plane_bump/good-output/snap.v.0000000201 b/test/beta_plane_bump/good-output/snap.v.0000000201
index 37d6c27c..3350b76d 100644
Binary files a/test/beta_plane_bump/good-output/snap.v.0000000201 and b/test/beta_plane_bump/good-output/snap.v.0000000201 differ
diff --git a/test/beta_plane_bump/good-output/snap.v.0000000401 b/test/beta_plane_bump/good-output/snap.v.0000000401
index a010fb3c..00f06d3c 100644
Binary files a/test/beta_plane_bump/good-output/snap.v.0000000401 and b/test/beta_plane_bump/good-output/snap.v.0000000401 differ
diff --git a/test/beta_plane_bump/parameters.in b/test/beta_plane_bump/parameters.in
deleted file mode 100644
index 7cfb9e8d..00000000
--- a/test/beta_plane_bump/parameters.in
+++ /dev/null
@@ -1,80 +0,0 @@
-! Parameter file. Change the values, but not the names.
-! 
-! Each namelist ends with a /
-! This is important and should not be deleted.
-!
-! au is viscosity
-! ah is thickness diffusivity
-! ar is linear drag between layers
-! dt is time step
-! slip is free-slip (=0), no-slip (=1), or partial slip (something in between)
-! nTimeSteps: number of timesteps before stopping
-! dumpFreq: frequency of snapshot output
-! avFreq: frequency of averaged output
-! hmin: minimum layer thickness allowed by model (for stability)
-! maxits: maximum iterations for the successive over relaxation algorithm. Should be at least max(nx,ny), and probably nx*ny
-! eps: convergence tolerance for SOR solver
-! freesurfFac: 1. = linear implicit free surface, 0. = rigid lid. So far all tests using freesurfFac = 1. have failed 
-! g is the gravity at interfaces (including surface). must have as many entries as there are layers
-! input files are where to look for the various inputs
-!
- &NUMERICS
- au = 500.,
- ah = 0.0,
- ar = 0.0,
- botDrag = 1e-6,
- dt = 600.,
- slip = 0.0,
- nTimeSteps = 502,
- dumpFreq = 1.2e5,
- avFreq = 1.2e5,
- hmin = 100,
- maxits = 1000,
- eps = 1e-2,
- freesurfFac = 0.,
- thickness_error = 1e-2,
- /
- &MODEL
- hmean = 400.,1600.,
- depthFile = '',
- H0 = 2000.,
- RedGrav = .FALSE.,
- /
- &SPONGE
- spongeHTimeScaleFile = '',
- spongeUTimeScaleFile = '',
- spongeVTimeScaleFile = '',
- spongeHFile = '',
- spongeUfile = '',
- spongeVfile = '',
- /
- &PHYSICS
- g_vec = 9.8, 0.01,
- rho0 = 1035.,
- /
- &GRID
- nx = 10,
- ny = 10,
- layers = 2,
- dx = 2e4,
- dy = 2e4,
- fUfile = 'input/fu.bin',
- fVfile = 'input/fv.bin',
- wetMaskFile = 'input/wetmask.bin',
- /
-! Inital conditions for u, v, and h
- &INITIAL_CONDITONS
- initUfile = '',
- initVfile = '',
- initHfile = 'input/initH.bin',
- initEtaFile = '',
- /
- &EXTERNAL_FORCING
- zonalWindFile = '',
- meridionalWindFile = '', 
- DumpWind = .FALSE.,
- wind_mag_time_series_file = '',
- /
-! input/initH.bin
-! 
- &
diff --git a/test/beta_plane_bump_red_grav/aronnax.conf b/test/beta_plane_bump_red_grav/aronnax.conf
new file mode 100644
index 00000000..88343d76
--- /dev/null
+++ b/test/beta_plane_bump_red_grav/aronnax.conf
@@ -0,0 +1,53 @@
+# Aronnax configuration file. Change the values, but not the names.
+#
+# au is viscosity
+# ah is thickness diffusivity
+# ar is linear drag between layers
+# dt is time step
+# slip is free-slip (=0), no-slip (=1), or partial slip (something in between)
+# nTimeSteps: number of timesteps before stopping
+# dumpFreq: frequency of snapshot output
+# avFreq: frequency of averaged output
+# hmin: minimum layer thickness allowed by model (for stability)
+# maxits: maximum iterations for the successive over relaxation algorithm. Should be at least max(nx,ny), and probably nx*ny
+# eps: convergence tolerance for SOR solver
+# freesurfFac: 1. = linear implicit free surface, 0. = rigid lid. So far all tests using freesurfFac = 1. have failed
+# g is the gravity at interfaces (including surface). must have as many entries as there are layers
+# input files are where to look for the various inputs
+
+[numerics]
+au = 500.
+ah = 0.0
+ar = 1e-8
+dt = 600.
+slip = 0.0
+nTimeSteps = 502
+dumpFreq = 6e4
+avFreq = 3e5
+hmin = 100
+maxits = 1000
+eps = 1e-2
+freesurfFac = 0.
+thickness_error = 1e-2
+
+[model]
+hmean = 400.
+H0 = 2000.
+RedGrav = yes
+
+[physics]
+g_vec = 0.01
+rho0 = 1035.
+
+[grid]
+nx = 10
+ny = 10
+layers = 1
+dx = 2e4
+dy = 2e4
+fUfile = :beta_plane_f_u:1e-5,2e-11
+fVfile = :beta_plane_f_v:1e-5,2e-11
+wetMaskFile = :rectangular_pool:
+
+[external_forcing]
+DumpWind = no
diff --git a/test/beta_plane_bump_red_grav/good-output/av.h.0000000001 b/test/beta_plane_bump_red_grav/good-output/av.h.0000000001
index 32897ab1..2e09f4df 100644
Binary files a/test/beta_plane_bump_red_grav/good-output/av.h.0000000001 and b/test/beta_plane_bump_red_grav/good-output/av.h.0000000001 differ
diff --git a/test/beta_plane_bump_red_grav/good-output/av.h.0000000501 b/test/beta_plane_bump_red_grav/good-output/av.h.0000000501
index c4fe5b27..85d6f74f 100644
Binary files a/test/beta_plane_bump_red_grav/good-output/av.h.0000000501 and b/test/beta_plane_bump_red_grav/good-output/av.h.0000000501 differ
diff --git a/test/beta_plane_bump_red_grav/good-output/av.u.0000000001 b/test/beta_plane_bump_red_grav/good-output/av.u.0000000001
index 953e7c2a..77e7a735 100644
Binary files a/test/beta_plane_bump_red_grav/good-output/av.u.0000000001 and b/test/beta_plane_bump_red_grav/good-output/av.u.0000000001 differ
diff --git a/test/beta_plane_bump_red_grav/good-output/av.u.0000000501 b/test/beta_plane_bump_red_grav/good-output/av.u.0000000501
index 07b7d908..36ce8b51 100644
Binary files a/test/beta_plane_bump_red_grav/good-output/av.u.0000000501 and b/test/beta_plane_bump_red_grav/good-output/av.u.0000000501 differ
diff --git a/test/beta_plane_bump_red_grav/good-output/av.v.0000000001 b/test/beta_plane_bump_red_grav/good-output/av.v.0000000001
index 54fe41da..0a13660c 100644
Binary files a/test/beta_plane_bump_red_grav/good-output/av.v.0000000001 and b/test/beta_plane_bump_red_grav/good-output/av.v.0000000001 differ
diff --git a/test/beta_plane_bump_red_grav/good-output/av.v.0000000501 b/test/beta_plane_bump_red_grav/good-output/av.v.0000000501
index f01071ae..c0b2a62a 100644
Binary files a/test/beta_plane_bump_red_grav/good-output/av.v.0000000501 and b/test/beta_plane_bump_red_grav/good-output/av.v.0000000501 differ
diff --git a/test/beta_plane_bump_red_grav/good-output/snap.h.0000000001 b/test/beta_plane_bump_red_grav/good-output/snap.h.0000000001
index 880dd493..88bcf26b 100644
Binary files a/test/beta_plane_bump_red_grav/good-output/snap.h.0000000001 and b/test/beta_plane_bump_red_grav/good-output/snap.h.0000000001 differ
diff --git a/test/beta_plane_bump_red_grav/good-output/snap.h.0000000101 b/test/beta_plane_bump_red_grav/good-output/snap.h.0000000101
index 4d7fee2b..198263e8 100644
Binary files a/test/beta_plane_bump_red_grav/good-output/snap.h.0000000101 and b/test/beta_plane_bump_red_grav/good-output/snap.h.0000000101 differ
diff --git a/test/beta_plane_bump_red_grav/good-output/snap.h.0000000201 b/test/beta_plane_bump_red_grav/good-output/snap.h.0000000201
index 05cc370b..ed6a5435 100644
Binary files a/test/beta_plane_bump_red_grav/good-output/snap.h.0000000201 and b/test/beta_plane_bump_red_grav/good-output/snap.h.0000000201 differ
diff --git a/test/beta_plane_bump_red_grav/good-output/snap.h.0000000301 b/test/beta_plane_bump_red_grav/good-output/snap.h.0000000301
index f70e60d8..f4ba0791 100644
Binary files a/test/beta_plane_bump_red_grav/good-output/snap.h.0000000301 and b/test/beta_plane_bump_red_grav/good-output/snap.h.0000000301 differ
diff --git a/test/beta_plane_bump_red_grav/good-output/snap.h.0000000401 b/test/beta_plane_bump_red_grav/good-output/snap.h.0000000401
index 7fe1a2e9..62919cea 100644
Binary files a/test/beta_plane_bump_red_grav/good-output/snap.h.0000000401 and b/test/beta_plane_bump_red_grav/good-output/snap.h.0000000401 differ
diff --git a/test/beta_plane_bump_red_grav/good-output/snap.h.0000000501 b/test/beta_plane_bump_red_grav/good-output/snap.h.0000000501
index c6ef3dbb..99548779 100644
Binary files a/test/beta_plane_bump_red_grav/good-output/snap.h.0000000501 and b/test/beta_plane_bump_red_grav/good-output/snap.h.0000000501 differ
diff --git a/test/beta_plane_bump_red_grav/good-output/snap.u.0000000001 b/test/beta_plane_bump_red_grav/good-output/snap.u.0000000001
index 06c43bce..12604f34 100644
Binary files a/test/beta_plane_bump_red_grav/good-output/snap.u.0000000001 and b/test/beta_plane_bump_red_grav/good-output/snap.u.0000000001 differ
diff --git a/test/beta_plane_bump_red_grav/good-output/snap.u.0000000101 b/test/beta_plane_bump_red_grav/good-output/snap.u.0000000101
index 8fc7b9c8..12bb80b8 100644
Binary files a/test/beta_plane_bump_red_grav/good-output/snap.u.0000000101 and b/test/beta_plane_bump_red_grav/good-output/snap.u.0000000101 differ
diff --git a/test/beta_plane_bump_red_grav/good-output/snap.u.0000000201 b/test/beta_plane_bump_red_grav/good-output/snap.u.0000000201
index 04f23f32..57d56564 100644
Binary files a/test/beta_plane_bump_red_grav/good-output/snap.u.0000000201 and b/test/beta_plane_bump_red_grav/good-output/snap.u.0000000201 differ
diff --git a/test/beta_plane_bump_red_grav/good-output/snap.u.0000000301 b/test/beta_plane_bump_red_grav/good-output/snap.u.0000000301
index 76a181c3..9988a5fd 100644
Binary files a/test/beta_plane_bump_red_grav/good-output/snap.u.0000000301 and b/test/beta_plane_bump_red_grav/good-output/snap.u.0000000301 differ
diff --git a/test/beta_plane_bump_red_grav/good-output/snap.u.0000000401 b/test/beta_plane_bump_red_grav/good-output/snap.u.0000000401
index 244fbe87..b383af26 100644
Binary files a/test/beta_plane_bump_red_grav/good-output/snap.u.0000000401 and b/test/beta_plane_bump_red_grav/good-output/snap.u.0000000401 differ
diff --git a/test/beta_plane_bump_red_grav/good-output/snap.u.0000000501 b/test/beta_plane_bump_red_grav/good-output/snap.u.0000000501
index 28054205..5abcb1cb 100644
Binary files a/test/beta_plane_bump_red_grav/good-output/snap.u.0000000501 and b/test/beta_plane_bump_red_grav/good-output/snap.u.0000000501 differ
diff --git a/test/beta_plane_bump_red_grav/good-output/snap.v.0000000001 b/test/beta_plane_bump_red_grav/good-output/snap.v.0000000001
index bb1ae846..09311e8d 100644
Binary files a/test/beta_plane_bump_red_grav/good-output/snap.v.0000000001 and b/test/beta_plane_bump_red_grav/good-output/snap.v.0000000001 differ
diff --git a/test/beta_plane_bump_red_grav/good-output/snap.v.0000000101 b/test/beta_plane_bump_red_grav/good-output/snap.v.0000000101
index 2d630f3e..4047a80b 100644
Binary files a/test/beta_plane_bump_red_grav/good-output/snap.v.0000000101 and b/test/beta_plane_bump_red_grav/good-output/snap.v.0000000101 differ
diff --git a/test/beta_plane_bump_red_grav/good-output/snap.v.0000000201 b/test/beta_plane_bump_red_grav/good-output/snap.v.0000000201
index 5c81c1bf..4434bbc4 100644
Binary files a/test/beta_plane_bump_red_grav/good-output/snap.v.0000000201 and b/test/beta_plane_bump_red_grav/good-output/snap.v.0000000201 differ
diff --git a/test/beta_plane_bump_red_grav/good-output/snap.v.0000000301 b/test/beta_plane_bump_red_grav/good-output/snap.v.0000000301
index 16ff9616..51200d42 100644
Binary files a/test/beta_plane_bump_red_grav/good-output/snap.v.0000000301 and b/test/beta_plane_bump_red_grav/good-output/snap.v.0000000301 differ
diff --git a/test/beta_plane_bump_red_grav/good-output/snap.v.0000000401 b/test/beta_plane_bump_red_grav/good-output/snap.v.0000000401
index 23241830..a69bd593 100644
Binary files a/test/beta_plane_bump_red_grav/good-output/snap.v.0000000401 and b/test/beta_plane_bump_red_grav/good-output/snap.v.0000000401 differ
diff --git a/test/beta_plane_bump_red_grav/good-output/snap.v.0000000501 b/test/beta_plane_bump_red_grav/good-output/snap.v.0000000501
index 51e89e42..0154000d 100644
Binary files a/test/beta_plane_bump_red_grav/good-output/snap.v.0000000501 and b/test/beta_plane_bump_red_grav/good-output/snap.v.0000000501 differ
diff --git a/test/beta_plane_bump_red_grav/parameters.in b/test/beta_plane_bump_red_grav/parameters.in
deleted file mode 100644
index b452b7ec..00000000
--- a/test/beta_plane_bump_red_grav/parameters.in
+++ /dev/null
@@ -1,79 +0,0 @@
-! Parameter file. Change the values, but not the names.
-! 
-! Each namelist ends with a /
-! This is important and should not be deleted.
-!
-! au is viscosity
-! ah is thickness diffusivity
-! ar is linear drag between layers
-! dt is time step
-! slip is free-slip (=0), no-slip (=1), or partial slip (something in between)
-! nTimeSteps: number of timesteps before stopping
-! dumpFreq: frequency of snapshot output
-! avFreq: frequency of averaged output
-! hmin: minimum layer thickness allowed by model (for stability)
-! maxits: maximum iterations for the successive over relaxation algorithm. Should be at least max(nx,ny), and probably nx*ny
-! eps: convergence tolerance for SOR solver
-! freesurfFac: 1. = linear implicit free surface, 0. = rigid lid. So far all tests using freesurfFac = 1. have failed 
-! g is the gravity at interfaces (including surface). must have as many entries as there are layers
-! input files are where to look for the various inputs
-!
- &NUMERICS
- au = 500.,
- ah = 0.0,
- ar = 1e-8,
- dt = 600.,
- slip = 0.0,
- nTimeSteps = 502,
- dumpFreq = 6e4,
- avFreq = 3e5,
- hmin = 100,
- maxits = 1000,
- eps = 1e-2,
- freesurfFac = 0.,
- thickness_error = 1e-2,
- /
- &MODEL
- hmean = 400.,
- depthFile = '',
- H0 = 2000.,
- RedGrav = .TRUE.,
- /
- &SPONGE
- spongeHTimeScaleFile = '',
- spongeUTimeScaleFile = '',
- spongeVTimeScaleFile = '',
- spongeHFile = '',
- spongeUfile = '',
- spongeVfile = '',
- /
- &PHYSICS
- g_vec = 0.01,
- rho0 = 1035.,
- /
- &GRID
- nx = 10,
- ny = 10,
- layers = 1,
- dx = 2e4,
- dy = 2e4,
- fUfile = 'input/fu.bin',
- fVfile = 'input/fv.bin',
- wetMaskFile = 'input/wetmask.bin',
- /
-! Inital conditions for u, v, and h
- &INITIAL_CONDITONS
- initUfile = '',
- initVfile = '',
- initHfile = 'input/initH.bin',
- initEtaFile = '',
- /
- &EXTERNAL_FORCING
- zonalWindFile = '',
- meridionalWindFile = '', 
- DumpWind = .FALSE.,
- wind_mag_time_series_file = '',
- /
-! input/initH.bin
-! 
- &
diff --git a/test/beta_plane_gyre/aronnax.conf b/test/beta_plane_gyre/aronnax.conf
new file mode 100644
index 00000000..0150a13c
--- /dev/null
+++ b/test/beta_plane_gyre/aronnax.conf
@@ -0,0 +1,59 @@
+# Aronnax configuration file. Change the values, but not the names.
+# 
+# au is viscosity
+# ah is thickness diffusivity
+# ar is linear drag between layers
+# dt is time step
+# slip is free-slip (=0), no-slip (=1), or partial slip (something in between)
+# nTimeSteps: number of timesteps before stopping
+# dumpFreq: frequency of snapshot output
+# avFreq: frequency of averaged output
+# hmin: minimum layer thickness allowed by model (for stability)
+# maxits: maximum iterations for the successive over relaxation algorithm. Should be at least max(nx,ny), and probably nx*ny
+# eps: convergence tolerance for SOR solver
+# freesurfFac: 1. = linear implicit free surface, 0. = rigid lid. So far all tests using freesurfFac = 1. have failed 
+# g is the gravity at interfaces (including surface). must have as many entries as there are layers
+# input files are where to look for the various inputs
+
+[numerics]
+au = 500.
+ah = 500.0,500.
+ar = 1e-8
+botDrag = 1e-6
+dt = 600.
+slip = 1.0
+nTimeSteps = 801
+dumpFreq = 12e4
+avFreq = 48e4
+hmin = 100
+maxits = 1000
+eps = 1e-2
+freesurfFac = 0.
+thickness_error = 1e-2
+
+[model]
+hmean = 600.,1400.
+H0 = 2000.
+RedGrav = no
+
+[physics]
+g_vec = 9.8, 0.01
+rho0 = 1035.
+
+[grid]
+nx = 10
+ny = 10
+layers = 2
+dx = 2e4
+dy = 2e4
+fUfile = :beta_plane_f_u:1e-5,2e-11
+fVfile = :beta_plane_f_v:1e-5,2e-11
+wetMaskFile = :rectangular_pool:
+
+# Inital conditions for h
+[initial_conditons]
+initHfile = :depths:600.0,1400.0
+
+[external_forcing]
+DumpWind = yes
+
diff --git a/test/beta_plane_gyre/good-output/av.eta.0000000001 b/test/beta_plane_gyre/good-output/av.eta.0000000001
index 4b5852dd..19ee1b1c 100644
Binary files a/test/beta_plane_gyre/good-output/av.eta.0000000001 and b/test/beta_plane_gyre/good-output/av.eta.0000000001 differ
diff --git a/test/beta_plane_gyre/good-output/av.eta.0000000801 b/test/beta_plane_gyre/good-output/av.eta.0000000801
index 8f9e58e8..8ae443d0 100644
Binary files a/test/beta_plane_gyre/good-output/av.eta.0000000801 and b/test/beta_plane_gyre/good-output/av.eta.0000000801 differ
diff --git a/test/beta_plane_gyre/good-output/av.h.0000000001 b/test/beta_plane_gyre/good-output/av.h.0000000001
index 19f8c67c..09aec92a 100644
Binary files a/test/beta_plane_gyre/good-output/av.h.0000000001 and b/test/beta_plane_gyre/good-output/av.h.0000000001 differ
diff --git a/test/beta_plane_gyre/good-output/av.h.0000000801 b/test/beta_plane_gyre/good-output/av.h.0000000801
index 8cf34b8f..56437ae4 100644
Binary files a/test/beta_plane_gyre/good-output/av.h.0000000801 and b/test/beta_plane_gyre/good-output/av.h.0000000801 differ
diff --git a/test/beta_plane_gyre/good-output/av.u.0000000001 b/test/beta_plane_gyre/good-output/av.u.0000000001
index 3ca79140..423967fe 100644
Binary files a/test/beta_plane_gyre/good-output/av.u.0000000001 and b/test/beta_plane_gyre/good-output/av.u.0000000001 differ
diff --git a/test/beta_plane_gyre/good-output/av.u.0000000801 b/test/beta_plane_gyre/good-output/av.u.0000000801
index 582c937a..46e4a498 100644
Binary files a/test/beta_plane_gyre/good-output/av.u.0000000801 and b/test/beta_plane_gyre/good-output/av.u.0000000801 differ
diff --git a/test/beta_plane_gyre/good-output/av.v.0000000001 b/test/beta_plane_gyre/good-output/av.v.0000000001
index 0c9ea215..493ef307 100644
Binary files a/test/beta_plane_gyre/good-output/av.v.0000000001 and b/test/beta_plane_gyre/good-output/av.v.0000000001 differ
diff --git a/test/beta_plane_gyre/good-output/av.v.0000000801 b/test/beta_plane_gyre/good-output/av.v.0000000801
index b547e6eb..654d7fa9 100644
Binary files a/test/beta_plane_gyre/good-output/av.v.0000000801 and b/test/beta_plane_gyre/good-output/av.v.0000000801 differ
diff --git a/test/beta_plane_gyre/good-output/snap.eta.0000000001 b/test/beta_plane_gyre/good-output/snap.eta.0000000001
index 71d2e1c0..1cf9e46a 100644
Binary files a/test/beta_plane_gyre/good-output/snap.eta.0000000001 and b/test/beta_plane_gyre/good-output/snap.eta.0000000001 differ
diff --git a/test/beta_plane_gyre/good-output/snap.eta.0000000201 b/test/beta_plane_gyre/good-output/snap.eta.0000000201
index c6406c14..6a1e2899 100644
Binary files a/test/beta_plane_gyre/good-output/snap.eta.0000000201 and b/test/beta_plane_gyre/good-output/snap.eta.0000000201 differ
diff --git a/test/beta_plane_gyre/good-output/snap.eta.0000000401 b/test/beta_plane_gyre/good-output/snap.eta.0000000401
index fc2f1fb0..73a442ac 100644
Binary files a/test/beta_plane_gyre/good-output/snap.eta.0000000401 and b/test/beta_plane_gyre/good-output/snap.eta.0000000401 differ
diff --git a/test/beta_plane_gyre/good-output/snap.eta.0000000601 b/test/beta_plane_gyre/good-output/snap.eta.0000000601
index 580a97e0..aea50835 100644
Binary files a/test/beta_plane_gyre/good-output/snap.eta.0000000601 and b/test/beta_plane_gyre/good-output/snap.eta.0000000601 differ
diff --git a/test/beta_plane_gyre/good-output/snap.eta.0000000801 b/test/beta_plane_gyre/good-output/snap.eta.0000000801
index 701e440e..7d4b46b0 100644
Binary files a/test/beta_plane_gyre/good-output/snap.eta.0000000801 and b/test/beta_plane_gyre/good-output/snap.eta.0000000801 differ
diff --git a/test/beta_plane_gyre/good-output/snap.h.0000000001 b/test/beta_plane_gyre/good-output/snap.h.0000000001
index f5a32e7d..183afec7 100644
Binary files a/test/beta_plane_gyre/good-output/snap.h.0000000001 and b/test/beta_plane_gyre/good-output/snap.h.0000000001 differ
diff --git a/test/beta_plane_gyre/good-output/snap.h.0000000201 b/test/beta_plane_gyre/good-output/snap.h.0000000201
index b763c465..ac1628ad 100644
Binary files a/test/beta_plane_gyre/good-output/snap.h.0000000201 and b/test/beta_plane_gyre/good-output/snap.h.0000000201 differ
diff --git a/test/beta_plane_gyre/good-output/snap.h.0000000401 b/test/beta_plane_gyre/good-output/snap.h.0000000401
index a734e09b..c11cc2c1 100644
Binary files a/test/beta_plane_gyre/good-output/snap.h.0000000401 and b/test/beta_plane_gyre/good-output/snap.h.0000000401 differ
diff --git a/test/beta_plane_gyre/good-output/snap.h.0000000601 b/test/beta_plane_gyre/good-output/snap.h.0000000601
index 8f4bfc6c..884cbb2f 100644
Binary files a/test/beta_plane_gyre/good-output/snap.h.0000000601 and b/test/beta_plane_gyre/good-output/snap.h.0000000601 differ
diff --git a/test/beta_plane_gyre/good-output/snap.h.0000000801 b/test/beta_plane_gyre/good-output/snap.h.0000000801
index 4a672ff1..bc8ce2bf 100644
Binary files a/test/beta_plane_gyre/good-output/snap.h.0000000801 and b/test/beta_plane_gyre/good-output/snap.h.0000000801 differ
diff --git a/test/beta_plane_gyre/good-output/snap.u.0000000001 b/test/beta_plane_gyre/good-output/snap.u.0000000001
index b82b9b87..925a4b43 100644
Binary files a/test/beta_plane_gyre/good-output/snap.u.0000000001 and b/test/beta_plane_gyre/good-output/snap.u.0000000001 differ
diff --git a/test/beta_plane_gyre/good-output/snap.u.0000000201 b/test/beta_plane_gyre/good-output/snap.u.0000000201
index 50205893..3682b9bf 100644
Binary files a/test/beta_plane_gyre/good-output/snap.u.0000000201 and b/test/beta_plane_gyre/good-output/snap.u.0000000201 differ
diff --git a/test/beta_plane_gyre/good-output/snap.u.0000000401 b/test/beta_plane_gyre/good-output/snap.u.0000000401
index 4aa873ef..31a0be44 100644
Binary files a/test/beta_plane_gyre/good-output/snap.u.0000000401 and b/test/beta_plane_gyre/good-output/snap.u.0000000401 differ
diff --git a/test/beta_plane_gyre/good-output/snap.u.0000000601 b/test/beta_plane_gyre/good-output/snap.u.0000000601
index b18e5b8a..815f3488 100644
Binary files a/test/beta_plane_gyre/good-output/snap.u.0000000601 and b/test/beta_plane_gyre/good-output/snap.u.0000000601 differ
diff --git a/test/beta_plane_gyre/good-output/snap.u.0000000801 b/test/beta_plane_gyre/good-output/snap.u.0000000801
index 8e763afb..2510b650 100644
Binary files a/test/beta_plane_gyre/good-output/snap.u.0000000801 and b/test/beta_plane_gyre/good-output/snap.u.0000000801 differ
diff --git a/test/beta_plane_gyre/good-output/snap.v.0000000001 b/test/beta_plane_gyre/good-output/snap.v.0000000001
index 2c662ffd..20d551d2 100644
Binary files a/test/beta_plane_gyre/good-output/snap.v.0000000001 and b/test/beta_plane_gyre/good-output/snap.v.0000000001 differ
diff --git a/test/beta_plane_gyre/good-output/snap.v.0000000201 b/test/beta_plane_gyre/good-output/snap.v.0000000201
index a4d99c02..ba2ba7e3 100644
Binary files a/test/beta_plane_gyre/good-output/snap.v.0000000201 and b/test/beta_plane_gyre/good-output/snap.v.0000000201 differ
diff --git a/test/beta_plane_gyre/good-output/snap.v.0000000401 b/test/beta_plane_gyre/good-output/snap.v.0000000401
index 11bd2ca8..e19b103d 100644
Binary files a/test/beta_plane_gyre/good-output/snap.v.0000000401 and b/test/beta_plane_gyre/good-output/snap.v.0000000401 differ
diff --git a/test/beta_plane_gyre/good-output/snap.v.0000000601 b/test/beta_plane_gyre/good-output/snap.v.0000000601
index e04ec273..a682c5e6 100644
Binary files a/test/beta_plane_gyre/good-output/snap.v.0000000601 and b/test/beta_plane_gyre/good-output/snap.v.0000000601 differ
diff --git a/test/beta_plane_gyre/good-output/snap.v.0000000801 b/test/beta_plane_gyre/good-output/snap.v.0000000801
index 0438824b..a52d1184 100644
Binary files a/test/beta_plane_gyre/good-output/snap.v.0000000801 and b/test/beta_plane_gyre/good-output/snap.v.0000000801 differ
diff --git a/test/beta_plane_gyre/parameters.in b/test/beta_plane_gyre/parameters.in
deleted file mode 100644
index e55102fe..00000000
--- a/test/beta_plane_gyre/parameters.in
+++ /dev/null
@@ -1,81 +0,0 @@
-! Parameter file. Change the values, but not the names.
-! 
-! Each namelist ends with a /
-! This is important and should not be deleted.
-!
-! au is viscosity
-! ah is thickness diffusivity
-! ar is linear drag between layers
-! dt is time step
-! slip is free-slip (=0), no-slip (=1), or partial slip (something in between)
-! nTimeSteps: number of timesteps before stopping
-! dumpFreq: frequency of snapshot output
-! avFreq: frequency of averaged output
-! hmin: minimum layer thickness allowed by model (for stability)
-! maxits: maximum iterations for the successive over relaxation algorithm. Should be at least max(nx,ny), and probably nx*ny
-! eps: convergence tolerance for SOR solver
-! freesurfFac: 1. = linear implicit free surface, 0. = rigid lid. So far all tests using freesurfFac = 1. have failed 
-! g is the gravity at interfaces (including surface). must have as many entries as there are layers
-! input files are where to look for the various inputs
-!
- &NUMERICS
- au = 500.,
- ah = 500.0,500.,
- ar = 1e-8,
- botDrag = 1e-6,
- dt = 600.,
- slip = 1.0,
- nTimeSteps = 801,
- dumpFreq = 12e4,
- avFreq = 48e4,
- hmin = 100,
- maxits = 1000,
- eps = 1e-2,
- freesurfFac = 0.,
- thickness_error = 1e-2,
- /
- &MODEL
- hmean = 600.,1400.,
- depthFile = '',
- H0 = 2000.,
- RedGrav = .FALSE.,
- /
- &SPONGE
- spongeHTimeScaleFile = '',
- spongeUTimeScaleFile = '',
- spongeVTimeScaleFile = '',
- spongeHFile = '',
- spongeUfile = '',
- spongeVfile = '',
- /
- &PHYSICS
- g_vec = 9.8, 0.01,
- rho0 = 1035.,
- /
- &GRID
- nx = 10,
- ny = 10,
- layers = 2,
- dx = 2e4,
- dy = 2e4,
- fUfile = 'input/fu.bin',
- fVfile = 'input/fv.bin',
- wetMaskFile = 'input/wetmask.bin',
- /
-! Inital conditions for u, v, and h
- &INITIAL_CONDITONS
- initUfile = '',
- initVfile = '',
- initHfile = '',
-!'input/initH.bin',
- initEtaFile = '',
- /
- &EXTERNAL_FORCING
- zonalWindFile = 'input/wind_x.bin',
- meridionalWindFile = '', 
- DumpWind = .TRUE.,
- wind_mag_time_series_file = '',
- /
-! input/initH.bin
-! 
- &
diff --git a/test/beta_plane_gyre_red_grav/aronnax.conf b/test/beta_plane_gyre_red_grav/aronnax.conf
new file mode 100644
index 00000000..7cb6d806
--- /dev/null
+++ b/test/beta_plane_gyre_red_grav/aronnax.conf
@@ -0,0 +1,58 @@
+# Parameter file. Change the values, but not the names.
+# 
+# au is viscosity
+# ah is thickness diffusivity
+# ar is linear drag between layers
+# dt is time step
+# slip is free-slip (=0), no-slip (=1), or partial slip (something in between)
+# nTimeSteps: number of timesteps before stopping
+# dumpFreq: frequency of snapshot output
+# avFreq: frequency of averaged output
+# hmin: minimum layer thickness allowed by model (for stability)
+# maxits: maximum iterations for the successive over relaxation algorithm. Should be at least max(nx,ny), and probably nx*ny
+# eps: convergence tolerance for SOR solver
+# freesurfFac: 1. = linear implicit free surface, 0. = rigid lid. So far all tests using freesurfFac = 1. have failed 
+# g is the gravity at interfaces (including surface). must have as many entries as there are layers
+# input files are where to look for the various inputs
+
+[numerics]
+au = 500.
+ah = 0.0
+ar = 1e-8
+dt = 600.
+slip = 1.0
+nTimeSteps = 8001
+dumpFreq = 12e5
+avFreq = 48e5
+hmin = 100
+maxits = 1000
+eps = 1e-2
+freesurfFac = 0.
+thickness_error = 1e-2
+
+[model]
+hmean = 400.
+H0 = 2000.
+RedGrav = yes
+
+[physics]
+g_vec = 0.01
+rho0 = 1035.
+
+[grid]
+nx = 10
+ny = 10
+layers = 1
+dx = 2e4
+dy = 2e4
+fUfile = :beta_plane_f_u:1e-5,2e-11
+fVfile = :beta_plane_f_v:1e-5,2e-11
+wetMaskFile = :rectangular_pool:
+
+# Inital conditions h
+[initial_conditons]
+initHfile = :depths:400.0
+
+[external_forcing]
+DumpWind = no
+
diff --git a/test/beta_plane_gyre_red_grav/good-output/av.h.0000000001 b/test/beta_plane_gyre_red_grav/good-output/av.h.0000000001
index 0487e5ce..9c7c3438 100644
Binary files a/test/beta_plane_gyre_red_grav/good-output/av.h.0000000001 and b/test/beta_plane_gyre_red_grav/good-output/av.h.0000000001 differ
diff --git a/test/beta_plane_gyre_red_grav/good-output/av.h.0000008001 b/test/beta_plane_gyre_red_grav/good-output/av.h.0000008001
index f70246d6..1ac62080 100644
Binary files a/test/beta_plane_gyre_red_grav/good-output/av.h.0000008001 and b/test/beta_plane_gyre_red_grav/good-output/av.h.0000008001 differ
diff --git a/test/beta_plane_gyre_red_grav/good-output/av.u.0000000001 b/test/beta_plane_gyre_red_grav/good-output/av.u.0000000001
index 7384fbc4..42c984c0 100644
Binary files a/test/beta_plane_gyre_red_grav/good-output/av.u.0000000001 and b/test/beta_plane_gyre_red_grav/good-output/av.u.0000000001 differ
diff --git a/test/beta_plane_gyre_red_grav/good-output/av.u.0000008001 b/test/beta_plane_gyre_red_grav/good-output/av.u.0000008001
index 6614796e..4865b64a 100644
Binary files a/test/beta_plane_gyre_red_grav/good-output/av.u.0000008001 and b/test/beta_plane_gyre_red_grav/good-output/av.u.0000008001 differ
diff --git a/test/beta_plane_gyre_red_grav/good-output/av.v.0000000001 b/test/beta_plane_gyre_red_grav/good-output/av.v.0000000001
index 3e87eea7..41ce1b04 100644
Binary files a/test/beta_plane_gyre_red_grav/good-output/av.v.0000000001 and b/test/beta_plane_gyre_red_grav/good-output/av.v.0000000001 differ
diff --git a/test/beta_plane_gyre_red_grav/good-output/av.v.0000008001 b/test/beta_plane_gyre_red_grav/good-output/av.v.0000008001
index 25498c48..aa82ba4f 100644
Binary files a/test/beta_plane_gyre_red_grav/good-output/av.v.0000008001 and b/test/beta_plane_gyre_red_grav/good-output/av.v.0000008001 differ
diff --git a/test/beta_plane_gyre_red_grav/good-output/snap.h.0000000001 b/test/beta_plane_gyre_red_grav/good-output/snap.h.0000000001
index 083ce22f..e86dbd0c 100644
Binary files a/test/beta_plane_gyre_red_grav/good-output/snap.h.0000000001 and b/test/beta_plane_gyre_red_grav/good-output/snap.h.0000000001 differ
diff --git a/test/beta_plane_gyre_red_grav/good-output/snap.h.0000002001 b/test/beta_plane_gyre_red_grav/good-output/snap.h.0000002001
index 5b960b0c..48f92298 100644
Binary files a/test/beta_plane_gyre_red_grav/good-output/snap.h.0000002001 and b/test/beta_plane_gyre_red_grav/good-output/snap.h.0000002001 differ
diff --git a/test/beta_plane_gyre_red_grav/good-output/snap.h.0000004001 b/test/beta_plane_gyre_red_grav/good-output/snap.h.0000004001
index 9c86635c..3bed314d 100644
Binary files a/test/beta_plane_gyre_red_grav/good-output/snap.h.0000004001 and b/test/beta_plane_gyre_red_grav/good-output/snap.h.0000004001 differ
diff --git a/test/beta_plane_gyre_red_grav/good-output/snap.h.0000006001 b/test/beta_plane_gyre_red_grav/good-output/snap.h.0000006001
index b1ec13c2..158b9a11 100644
Binary files a/test/beta_plane_gyre_red_grav/good-output/snap.h.0000006001 and b/test/beta_plane_gyre_red_grav/good-output/snap.h.0000006001 differ
diff --git a/test/beta_plane_gyre_red_grav/good-output/snap.h.0000008001 b/test/beta_plane_gyre_red_grav/good-output/snap.h.0000008001
index 620a3122..8d978e05 100644
Binary files a/test/beta_plane_gyre_red_grav/good-output/snap.h.0000008001 and b/test/beta_plane_gyre_red_grav/good-output/snap.h.0000008001 differ
diff --git a/test/beta_plane_gyre_red_grav/good-output/snap.u.0000000001 b/test/beta_plane_gyre_red_grav/good-output/snap.u.0000000001
index e40ecf97..9421e6b4 100644
Binary files a/test/beta_plane_gyre_red_grav/good-output/snap.u.0000000001 and b/test/beta_plane_gyre_red_grav/good-output/snap.u.0000000001 differ
diff --git a/test/beta_plane_gyre_red_grav/good-output/snap.u.0000002001 b/test/beta_plane_gyre_red_grav/good-output/snap.u.0000002001
index cae56c39..2a89dcb7 100644
Binary files a/test/beta_plane_gyre_red_grav/good-output/snap.u.0000002001 and b/test/beta_plane_gyre_red_grav/good-output/snap.u.0000002001 differ
diff --git a/test/beta_plane_gyre_red_grav/good-output/snap.u.0000004001 b/test/beta_plane_gyre_red_grav/good-output/snap.u.0000004001
index ec7f1c8c..d4f2a61c 100644
Binary files a/test/beta_plane_gyre_red_grav/good-output/snap.u.0000004001 and b/test/beta_plane_gyre_red_grav/good-output/snap.u.0000004001 differ
diff --git a/test/beta_plane_gyre_red_grav/good-output/snap.u.0000006001 b/test/beta_plane_gyre_red_grav/good-output/snap.u.0000006001
index 950a75fd..aad9cf19 100644
Binary files a/test/beta_plane_gyre_red_grav/good-output/snap.u.0000006001 and b/test/beta_plane_gyre_red_grav/good-output/snap.u.0000006001 differ
diff --git a/test/beta_plane_gyre_red_grav/good-output/snap.u.0000008001 b/test/beta_plane_gyre_red_grav/good-output/snap.u.0000008001
index 2c1d4958..0a77b19b 100644
Binary files a/test/beta_plane_gyre_red_grav/good-output/snap.u.0000008001 and b/test/beta_plane_gyre_red_grav/good-output/snap.u.0000008001 differ
diff --git a/test/beta_plane_gyre_red_grav/good-output/snap.v.0000000001 b/test/beta_plane_gyre_red_grav/good-output/snap.v.0000000001
index cbd8c4b6..d8f591c9 100644
Binary files a/test/beta_plane_gyre_red_grav/good-output/snap.v.0000000001 and b/test/beta_plane_gyre_red_grav/good-output/snap.v.0000000001 differ
diff --git a/test/beta_plane_gyre_red_grav/good-output/snap.v.0000002001 b/test/beta_plane_gyre_red_grav/good-output/snap.v.0000002001
index a8148221..80f4758a 100644
Binary files a/test/beta_plane_gyre_red_grav/good-output/snap.v.0000002001 and b/test/beta_plane_gyre_red_grav/good-output/snap.v.0000002001 differ
diff --git a/test/beta_plane_gyre_red_grav/good-output/snap.v.0000004001 b/test/beta_plane_gyre_red_grav/good-output/snap.v.0000004001
index 0ef4d553..afb776fb 100644
Binary files a/test/beta_plane_gyre_red_grav/good-output/snap.v.0000004001 and b/test/beta_plane_gyre_red_grav/good-output/snap.v.0000004001 differ
diff --git a/test/beta_plane_gyre_red_grav/good-output/snap.v.0000006001 b/test/beta_plane_gyre_red_grav/good-output/snap.v.0000006001
index 86222fdd..6f561580 100644
Binary files a/test/beta_plane_gyre_red_grav/good-output/snap.v.0000006001 and b/test/beta_plane_gyre_red_grav/good-output/snap.v.0000006001 differ
diff --git a/test/beta_plane_gyre_red_grav/good-output/snap.v.0000008001 b/test/beta_plane_gyre_red_grav/good-output/snap.v.0000008001
index a3adfd59..081f62f2 100644
Binary files a/test/beta_plane_gyre_red_grav/good-output/snap.v.0000008001 and b/test/beta_plane_gyre_red_grav/good-output/snap.v.0000008001 differ
diff --git a/test/beta_plane_gyre_red_grav/parameters.in b/test/beta_plane_gyre_red_grav/parameters.in
deleted file mode 100644
index 95ee246a..00000000
--- a/test/beta_plane_gyre_red_grav/parameters.in
+++ /dev/null
@@ -1,79 +0,0 @@
-! Parameter file. Change the values, but not the names.
-! 
-! Each namelist ends with a /
-! This is important and should not be deleted.
-!
-! au is viscosity
-! ah is thickness diffusivity
-! ar is linear drag between layers
-! dt is time step
-! slip is free-slip (=0), no-slip (=1), or partial slip (something in between)
-! nTimeSteps: number of timesteps before stopping
-! dumpFreq: frequency of snapshot output
-! avFreq: frequency of averaged output
-! hmin: minimum layer thickness allowed by model (for stability)
-! maxits: maximum iterations for the successive over relaxation algorithm. Should be at least max(nx,ny), and probably nx*ny
-! eps: convergence tolerance for SOR solver
-! freesurfFac: 1. = linear implicit free surface, 0. = rigid lid. So far all tests using freesurfFac = 1. have failed 
-! g is the gravity at interfaces (including surface). must have as many entries as there are layers
-! input files are where to look for the various inputs
-!
- &NUMERICS
- au = 500.,
- ah = 0.0,
- ar = 1e-8,
- dt = 600.,
- slip = 1.0,
- nTimeSteps = 8001,
- dumpFreq = 12e5,
- avFreq = 48e5,
- hmin = 100,
- maxits = 1000,
- eps = 1e-2,
- freesurfFac = 0.,
- thickness_error = 1e-2,
- /
- &MODEL
- hmean = 400.,
- depthFile = '',
- H0 = 2000.,
- RedGrav = .TRUE.,
- /
- &SPONGE
- spongeHTimeScaleFile = '',
- spongeUTimeScaleFile = '',
- spongeVTimeScaleFile = '',
- spongeHFile = '',
- spongeUfile = '',
- spongeVfile = '',
- /
- &PHYSICS
- g_vec = 0.01,
- rho0 = 1035.,
- /
- &GRID
- nx = 10,
- ny = 10,
- layers = 1,
- dx = 2e4,
- dy = 2e4,
- fUfile = 'input/fu.bin',
- fVfile = 'input/fv.bin',
- wetMaskFile = 'input/wetmask.bin',
- /
-! Inital conditions for u, v, and h
- &INITIAL_CONDITONS
- initUfile = '',
- initVfile = '',
- initHfile = 'input/initH.bin',
- initEtaFile = '',
- /
- &EXTERNAL_FORCING
- zonalWindFile = 'input/wind_x.bin',
- meridionalWindFile = '', 
- DumpWind = .FALSE.,
- wind_mag_time_series_file = '',
- /
-! input/initH.bin
-! 
- &
diff --git a/test/f_plane/aronnax.conf b/test/f_plane/aronnax.conf
new file mode 100644
index 00000000..ac46d35e
--- /dev/null
+++ b/test/f_plane/aronnax.conf
@@ -0,0 +1,59 @@
+# Aronnax configuration file. Change the values, but not the names.
+# 
+# au is viscosity
+# ah is thickness diffusivity
+# ar is linear drag between layers
+# dt is time step
+# slip is free-slip (=0), no-slip (=1), or partial slip (something in between)
+# nTimeSteps: number of timesteps before stopping
+# dumpFreq: frequency of snapshot output
+# avFreq: frequency of averaged output
+# hmin: minimum layer thickness allowed by model (for stability)
+# maxits: maximum iterations for the successive over relaxation algorithm. Should be at least max(nx,ny), and probably nx*ny
+# eps: convergence tolerance for SOR solver
+# freesurfFac: 1. = linear implicit free surface, 0. = rigid lid. So far all tests using freesurfFac = 1. have failed 
+# g is the gravity at interfaces (including surface). must have as many entries as there are layers
+# input files are where to look for the various inputs
+
+[numerics]
+au = 500.
+ah = 0.0
+ar = 1e-8
+botDrag = 1e-6
+dt = 100.
+slip = 0.0
+nTimeSteps = 61
+dumpFreq = 2e3
+avFreq = 6e3
+hmin = 100
+maxits = 500
+eps = 1e-2
+freesurfFac = 0.
+thickness_error = 1e-2
+
+[model]
+hmean = 400.,1600.
+H0 = 2000.
+RedGrav = no
+
+[physics]
+g_vec = 9.8, 0.01
+rho0 = 1035.
+
+[grid]
+nx = 10
+ny = 10
+layers = 2
+dx = 2e4
+dy = 2e4
+fUfile = :f_plane_f_u:10e-4
+fVfile = :f_plane_f_v:10e-4
+wetMaskFile = :rectangular_pool:
+
+# Inital conditions for h
+[initial_conditons]
+initHfile = :depths:400.0,1600.0
+
+[external_forcing]
+DumpWind = no
+
diff --git a/test/f_plane/parameters.in b/test/f_plane/parameters.in
deleted file mode 100644
index 7a771650..00000000
--- a/test/f_plane/parameters.in
+++ /dev/null
@@ -1,80 +0,0 @@
-! Parameter file. Change the values, but not the names.
-! 
-! Each namelist ends with a /
-! This is important and should not be deleted.
-!
-! au is viscosity
-! ah is thickness diffusivity
-! ar is linear drag between layers
-! dt is time step
-! slip is free-slip (=0), no-slip (=1), or partial slip (something in between)
-! nTimeSteps: number of timesteps before stopping
-! dumpFreq: frequency of snapshot output
-! avFreq: frequency of averaged output
-! hmin: minimum layer thickness allowed by model (for stability)
-! maxits: maximum iterations for the successive over relaxation algorithm. Should be at least max(nx,ny), and probably nx*ny
-! eps: convergence tolerance for SOR solver
-! freesurfFac: 1. = linear implicit free surface, 0. = rigid lid. So far all tests using freesurfFac = 1. have failed 
-! g is the gravity at interfaces (including surface). must have as many entries as there are layers
-! input files are where to look for the various inputs
-!
- &NUMERICS
- au = 500.,
- ah = 0.0,
- ar = 1e-8,
- botDrag = 1e-6,
- dt = 100.,
- slip = 0.0,
- nTimeSteps = 61,
- dumpFreq = 2e3,
- avFreq = 6e3,
- hmin = 100,
- maxits = 500,
- eps = 1e-2,
- freesurfFac = 0.,
- thickness_error = 1e-2,
- /
- &MODEL
- hmean = 400.,1600.,
- depthFile = '',
- H0 = 2000.,
- RedGrav = .FALSE.,
- /
- &SPONGE
- spongeHTimeScaleFile = '',
- spongeUTimeScaleFile = '',
- spongeVTimeScaleFile = '',
- spongeHFile = '',
- spongeUfile = '',
- spongeVfile = '',
- /
- &PHYSICS
- g_vec = 9.8, 0.01,
- rho0 = 1035.,
- /
- &GRID
- nx = 10,
- ny = 10,
- layers = 2,
- dx = 2e4,
- dy = 2e4,
- fUfile = 'input/fu.bin',
- fVfile = 'input/fv.bin',
- wetMaskFile = 'input/wetmask.bin',
- /
-! Inital conditions for u, v, and h
- &INITIAL_CONDITONS
- initUfile = '',
- initVfile = '',
- initHfile = 'input/initH.bin',
- initEtaFile = '',
- /
- &EXTERNAL_FORCING
- zonalWindFile = '',
- meridionalWindFile = '', 
- DumpWind = .FALSE.,
- wind_mag_time_series_file = '',
- /
-! input/initH.bin
-! 
- &
diff --git a/test/f_plane_red_grav/aronnax.conf b/test/f_plane_red_grav/aronnax.conf
new file mode 100644
index 00000000..e58e04a3
--- /dev/null
+++ b/test/f_plane_red_grav/aronnax.conf
@@ -0,0 +1,58 @@
+# Aronnax configuration file. Change the values, but not the names.
+# 
+# au is viscosity
+# ah is thickness diffusivity
+# ar is linear drag between layers
+# dt is time step
+# slip is free-slip (=0), no-slip (=1), or partial slip (something in between)
+# nTimeSteps: number of timesteps before stopping
+# dumpFreq: frequency of snapshot output
+# avFreq: frequency of averaged output
+# hmin: minimum layer thickness allowed by model (for stability)
+# maxits: maximum iterations for the successive over relaxation algorithm. Should be at least max(nx,ny), and probably nx*ny
+# eps: convergence tolerance for SOR solver
+# freesurfFac: 1. = linear implicit free surface, 0. = rigid lid. So far all tests using freesurfFac = 1. have failed 
+# g is the gravity at interfaces (including surface). must have as many entries as there are layers
+# input files are where to look for the various inputs
+
+[numerics]
+au = 500.
+ah = 0.0
+ar = 1e-8
+dt = 600.
+slip = 0.0
+nTimeSteps = 10001
+dumpFreq = 3e6
+avFreq = 3e6
+hmin = 100
+maxits = 1000
+eps = 1e-2
+freesurfFac = 0.
+thickness_error = 1e-2
+
+[model]
+hmean = 400.
+H0 = 2000.
+RedGrav = yes
+
+[physics]
+g_vec = 0.01
+rho0 = 1035.
+
+[grid]
+nx = 10
+ny = 10
+layers = 1
+dx = 2e4
+dy = 2e4
+fUfile = :f_plane_f_u:10e-4
+fVfile = :f_plane_f_v:10e-4
+wetMaskFile = :rectangular_pool:
+
+# Inital conditions for h
+[initial_conditons]
+initHfile = :depths:400.0
+
+[external_forcing]
+DumpWind = no
+
diff --git a/test/f_plane_red_grav/parameters.in b/test/f_plane_red_grav/parameters.in
deleted file mode 100644
index 8663c323..00000000
--- a/test/f_plane_red_grav/parameters.in
+++ /dev/null
@@ -1,79 +0,0 @@
-! Parameter file. Change the values, but not the names.
-! 
-! Each namelist ends with a /
-! This is important and should not be deleted.
-!
-! au is viscosity
-! ah is thickness diffusivity
-! ar is linear drag between layers
-! dt is time step
-! slip is free-slip (=0), no-slip (=1), or partial slip (something in between)
-! nTimeSteps: number of timesteps before stopping
-! dumpFreq: frequency of snapshot output
-! avFreq: frequency of averaged output
-! hmin: minimum layer thickness allowed by model (for stability)
-! maxits: maximum iterations for the successive over relaxation algorithm. Should be at least max(nx,ny), and probably nx*ny
-! eps: convergence tolerance for SOR solver
-! freesurfFac: 1. = linear implicit free surface, 0. = rigid lid. So far all tests using freesurfFac = 1. have failed 
-! g is the gravity at interfaces (including surface). must have as many entries as there are layers
-! input files are where to look for the various inputs
-!
- &NUMERICS
- au = 500.,
- ah = 0.0,
- ar = 1e-8,
- dt = 600.,
- slip = 0.0,
- nTimeSteps = 10001,
- dumpFreq = 3e6,
- avFreq = 3e6,
- hmin = 100,
- maxits = 1000,
- eps = 1e-2,
- freesurfFac = 0.,
- thickness_error = 1e-2,
- /
- &MODEL
- hmean = 400.,
- depthFile = '',
- H0 = 2000.,
- RedGrav = .TRUE.,
- /
- &SPONGE
- spongeHTimeScaleFile = '',
- spongeUTimeScaleFile = '',
- spongeVTimeScaleFile = '',
- spongeHFile = '',
- spongeUfile = '',
- spongeVfile = '',
- /
- &PHYSICS
- g_vec = 0.01,
- rho0 = 1035.,
- /
- &GRID
- nx = 10,
- ny = 10,
- layers = 1,
- dx = 2e4,
- dy = 2e4,
- fUfile = 'input/fu.bin',
- fVfile = 'input/fv.bin',
- wetMaskFile = 'input/wetmask.bin',
- /
-! Inital conditions for u, v, and h
- &INITIAL_CONDITONS
- initUfile = '',
- initVfile = '',
- initHfile = 'input/initH.bin',
- initEtaFile = '',
- /
- &EXTERNAL_FORCING
- zonalWindFile = '',
- meridionalWindFile = '', 
- DumpWind = .FALSE.,
- wind_mag_time_series_file = '',
- /
-! input/initH.bin
-! 
- &
diff --git a/test/output_preservation_test.py b/test/output_preservation_test.py
index b0be41ba..16dc5381 100644
--- a/test/output_preservation_test.py
+++ b/test/output_preservation_test.py
@@ -1,4 +1,3 @@
-from contextlib import contextmanager
 import os
 import os.path as p
 import subprocess as sub
@@ -9,61 +8,13 @@
 import numpy as np
 
 import aronnax as aro
+import aronnax.driver as drv
+from aronnax.utils import working_directory
 
 self_path = p.dirname(p.abspath(__file__))
-root_path = p.dirname(self_path)
 
 ### General helpers
 
-@contextmanager
-def working_directory(path):
-    old_path = os.getcwd()
-    sub.check_call(["mkdir", "-p", path])
-    os.chdir(path)
-    try:
-        yield
-    finally:
-        os.chdir(old_path)
-
-def tweak_parameters(nx, ny, layers):
-    sub.check_call(
-        "cat parameters.in " +
-        "| sed 's/^ nx =.*,$/ nx = %d,/'" % (nx,) +
-        "| sed 's/^ ny =.*,$/ ny = %d,/'" % (ny,) +
-        "| sed 's/^ layers =.*,$/ layers = %d,/'" % (layers,) +
-        "> parameters.new", shell=True)
-    sub.check_call(["mv", "parameters.new", "parameters.in"])
-
-def run_experiment(write_input, nx, ny, layers, aro_exec=None, valgrind=False, perf=False):
-    if aro_exec is None:
-        aro_exec = "aronnax_test"
-    sub.check_call(["rm", "-rf", "input/"])
-    sub.check_call(["rm", "-rf", "output/"])
-    sub.check_call(["mkdir", "-p", "output/"])
-    with working_directory(root_path):
-        sub.check_call(["make", aro_exec])
-    with working_directory("input"):
-        write_input(nx, ny, layers)
-    tweak_parameters(nx, ny, layers)
-    then = time.time()
-    env = dict(os.environ, GFORTRAN_STDERR_UNIT="17")
-    if valgrind or 'ARONNAX_TEST_VALGRIND_ALL' in os.environ:
-        assert not perf
-        sub.check_call(["valgrind", "--error-exitcode=5", p.join(root_path, aro_exec)],
-            env=env)
-    elif perf:
-        perf_cmds = ["perf", "stat", "-e", "r530010", # "flops", on my CPU.
-            "-e", "L1-dcache-loads", "-e", "L1-dcache-load-misses",
-            "-e", "L1-dcache-stores", "-e", "L1-dcache-store-misses",
-            "-e", "L1-icache-loads", "-e", "L1-icache-misses",
-            "-e", "L1-dcache-prefetches",
-            "-e", "branch-instructions", "-e", "branch-misses"]
-        sub.check_call(perf_cmds + [p.join(root_path, aro_exec)], env=env)
-    else:
-        sub.check_call([p.join(root_path, aro_exec)], env=env)
-    run_time = time.time() - then
-    return run_time
-
 def array_relative_error(a1, a2):
     """Return the elementwise absolute difference between the inputs,
 scaled by the maximum value that occurs in the input."""
@@ -80,7 +31,12 @@ def assert_outputs_close(nx, ny, layers, rtol):
     for outfile in outfiles:
         ans = aro.interpret_raw_file(p.join("output/", outfile), nx, ny, layers)
         good_ans = aro.interpret_raw_file(p.join("good-output/", outfile), nx, ny, layers)
-        assert np.amax(array_relative_error(ans, good_ans)) < rtol
+        relerr = np.amax(array_relative_error(ans, good_ans))
+        if relerr >= rtol:
+            print outfile
+            print ans
+            print good_ans
+        assert relerr < rtol
 
 def assert_volume_conservation(nx,ny,layers,rtol):
     hfiles = sorted(glob.glob("output/snap.h.*"))
@@ -100,106 +56,67 @@ def assert_volume_conservation(nx,ny,layers,rtol):
 
 ### The test cases themselves
 
-def write_input_f_plane_red_grav(nx, ny, layers):
-    assert layers == 1
+def test_f_plane_red_grav():
     xlen = 1e6
     ylen = 1e6
-    grid = aro.Grid(nx, ny, xlen / nx, ylen / ny)
-    aro.write_f_plane(nx, ny, 10e-4)
-    aro.write_rectangular_pool(nx, ny)
-    aro.write_initial_heights(grid, [400.0])
-
-def test_f_plane_red_grav():
     with working_directory(p.join(self_path, "f_plane_red_grav")):
-        run_experiment(write_input_f_plane_red_grav, 10, 10, 1)
+        drv.simulate(exe="aronnax_test",
+            nx=10, ny=10, dx=xlen/10, dy=ylen/10)
         assert_outputs_close(10, 10, 1, 1e-15)
         assert_volume_conservation(10, 10, 1, 1e-5)
 
-def write_input_f_plane(nx, ny, layers):
-    assert layers == 2
+def test_f_plane():
     xlen = 1e6
     ylen = 1e6
-    grid = aro.Grid(nx, ny, xlen / nx, ylen / ny)
-    aro.write_f_plane(nx, ny, 10e-4)
-    aro.write_rectangular_pool(nx, ny)
-    aro.write_initial_heights(grid, [400.0, 1600.0])
-
-def test_f_plane():
     with working_directory(p.join(self_path, "f_plane")):
-        run_experiment(write_input_f_plane, 10, 10, 2)
+        drv.simulate(exe="aronnax_test",
+            nx=10, ny=10, dx=xlen/10, dy=ylen/10)
         assert_outputs_close(10, 10, 2, 1e-15)
         assert_volume_conservation(10, 10, 2, 1e-5)
 
-def write_input_beta_plane_bump_red_grav(nx, ny, layers):
-    assert layers == 1
-    xlen = 1e6
-    ylen = 1e6
-    grid = aro.Grid(nx, ny, xlen / nx, ylen / ny)
-
-    aro.write_beta_plane(grid, 1e-5, 2e-11)
-    aro.write_rectangular_pool(nx, ny)
-    def bump(X, Y):
-        return 500. + 20*np.exp(-((6e5-X)**2 + (5e5-Y)**2)/(2*1e5**2))
-    aro.write_initial_heights(grid, [bump])
+def bump(X, Y):
+    return 500. + 20*np.exp(-((6e5-X)**2 + (5e5-Y)**2)/(2*1e5**2))
 
 def test_gaussian_bump_red_grav():
+    xlen = 1e6
+    ylen = 1e6
     with working_directory(p.join(self_path, "beta_plane_bump_red_grav")):
-        run_experiment(write_input_beta_plane_bump_red_grav, 10, 10, 1)
+        drv.simulate(initHfile=[bump], exe="aronnax_test",
+                     nx=10, ny=10, dx=xlen/10, dy=ylen/10)
         assert_outputs_close(10, 10, 1, 1.5e-13)
         assert_volume_conservation(10, 10, 1, 1e-5)
 
-def write_input_beta_plane_bump(nx, ny, layers):
-    assert layers == 2
+def test_gaussian_bump():
     xlen = 1e6
     ylen = 1e6
-    grid = aro.Grid(nx, ny, xlen / nx, ylen / ny)
-
-    aro.write_beta_plane(grid, 1e-5, 2e-11)
-    aro.write_rectangular_pool(nx, ny)
-    def bump(X, Y):
-        return 500. + 20*np.exp(-((6e5-X)**2 + (5e5-Y)**2)/(2*1e5**2))
-    aro.write_initial_heights(grid, [bump, lambda X, Y: 2000. - bump(X, Y)])
-
-def test_gaussian_bump():
     with working_directory(p.join(self_path, "beta_plane_bump")):
-        run_experiment(write_input_beta_plane_bump, 10, 10, 2)
+        drv.simulate(initHfile=[bump, lambda X, Y: 2000. - bump(X, Y)],
+                     nx=10, ny=10, exe="aronnax_test", dx=xlen/10, dy=ylen/10)
         assert_outputs_close(10, 10, 2, 2e-13)
         assert_volume_conservation(10, 10, 2, 1e-5)
 
-def write_input_beta_plane_gyre_red_grav(nx, ny, layers):
-    assert layers == 1
+def test_beta_plane_gyre_red_grav():
     xlen = 1e6
     ylen = 2e6
+    nx = 10; ny = 10
     grid = aro.Grid(nx, ny, xlen / nx, ylen / ny)
-
-    aro.write_beta_plane(grid, 1e-5, 2e-11)
-    aro.write_rectangular_pool(nx, ny)
-    aro.write_initial_heights(grid, [400.0])
     def wind(_, Y):
         return 0.05 * (1 - np.cos(2*np.pi * Y/np.max(grid.y)))
-    aro.write_wind_x(grid, wind)
-
-def test_beta_plane_gyre_red_grav():
     with working_directory(p.join(self_path, "beta_plane_gyre_red_grav")):
-        run_experiment(write_input_beta_plane_gyre_red_grav, 10, 10, 1, valgrind=True)
+        drv.simulate(zonalWindFile=wind, valgrind=True,
+                     nx=10, ny=10, exe="aronnax_test", dx=xlen/10, dy=ylen/10)
         assert_outputs_close(10, 10, 1, 2e-13)
         assert_volume_conservation(10, 10, 1, 1e-5)
 
-def write_input_beta_plane_gyre(nx, ny, layers):
-    assert layers == 2
+def test_beta_plane_gyre():
     xlen = 1e6
     ylen = 2e6
+    nx = 10; ny = 10
     grid = aro.Grid(nx, ny, xlen / nx, ylen / ny)
-
-    aro.write_beta_plane(grid, 1e-5, 2e-11)
-    aro.write_rectangular_pool(nx, ny)
-    aro.write_initial_heights(grid, [600.0, 1400.0])
     def wind(_, Y):
         return 0.05 * (1 - np.cos(2*np.pi * Y/np.max(grid.y)))
-    aro.write_wind_x(grid, wind)
-
-def test_beta_plane_gyre():
     with working_directory(p.join(self_path, "beta_plane_gyre")):
-        run_experiment(write_input_beta_plane_gyre, 10, 10, 2, valgrind=True)
+        drv.simulate(zonalWindFile=wind, valgrind=True,
+                     nx=10, ny=10, exe="aronnax_test", dx=xlen/10, dy=ylen/10)
         assert_outputs_close(10, 10, 2, 3e-12)
         assert_volume_conservation(10, 10, 2, 1e-5)