Repository: barbagroup/geoclaw-landspill (language: python)

Issue:
Prepare for JOSS submission
- [x] `CODE_OF_CONDUCT.md` (GitHub has some templates)
- [x] `CONTRIBUTING.md` There should be clear guidelines for third-parties wishing to:
    * Contribute to the software
    * Report issues or problems with the software
    * Seek support

- [x] `paper.md` as file name, including a title, summary, author names, affiliations, and key references
    * A list of the authors of the software and their affiliations, using the correct format (see the example below).
    * A summary describing the high-level functionality and purpose of the software for a diverse, non-specialist audience.
    * A clear Statement of Need that illustrates the research purpose of the software.
    * A list of key references, including to other software addressing related needs.
    * Mention (if applicable) a representative set of past or ongoing research projects using the software and recent scholarly publications enabled by it.
    * Acknowledgement of any financial support.

- [x] `paper.bib`
- [x] `codemeta.json`: create a metadata file describing your software and include it in your repository. (use [this script](https://gist.github.com/arfon/478b2ed49e11f984d6fb) to generate it)
- [x] A  markdown to explain examples, how to run examples, and what results are expected .

Reference:
1. https://joss.readthedocs.io/en/latest/submitting.html
2. https://joss.readthedocs.io/en/latest/review_criteria.html


Discussion / hints:

Fix patch:
diff --git a/CODE_OF_CONDUCT.md b/CODE_OF_CONDUCT.md
new file mode 100644
index 0000000..9105775
--- /dev/null
+++ b/CODE_OF_CONDUCT.md
@@ -0,0 +1,134 @@
+
+# Contributor Covenant Code of Conduct
+
+## Our Pledge
+
+We as members, contributors, and leaders pledge to make participation in our
+community a harassment-free experience for everyone, regardless of age, body
+size, visible or invisible disability, ethnicity, sex characteristics, gender
+identity and expression, level of experience, education, socio-economic status,
+nationality, personal appearance, race, religion, or sexual identity
+and orientation.
+
+We pledge to act and interact in ways that contribute to an open, welcoming,
+diverse, inclusive, and healthy community.
+
+## Our Standards
+
+Examples of behavior that contributes to a positive environment for our
+community include:
+
+* Demonstrating empathy and kindness toward other people
+* Being respectful of differing opinions, viewpoints, and experiences
+* Giving and gracefully accepting constructive feedback
+* Accepting responsibility and apologizing to those affected by our mistakes,
+  and learning from the experience
+* Focusing on what is best not just for us as individuals, but for the
+  overall community
+
+Examples of unacceptable behavior include:
+
+* The use of sexualized language or imagery, and sexual attention or
+  advances of any kind
+* Trolling, insulting or derogatory comments, and personal or political attacks
+* Public or private harassment
+* Publishing others' private information, such as a physical or email
+  address, without their explicit permission
+* Other conduct which could reasonably be considered inappropriate in a
+  professional setting
+
+## Enforcement Responsibilities
+
+Community leaders are responsible for clarifying and enforcing our standards of
+acceptable behavior and will take appropriate and fair corrective action in
+response to any behavior that they deem inappropriate, threatening, offensive,
+or harmful.
+
+Community leaders have the right and responsibility to remove, edit, or reject
+comments, commits, code, wiki edits, issues, and other contributions that are
+not aligned to this Code of Conduct, and will communicate reasons for moderation
+decisions when appropriate.
+
+## Scope
+
+This Code of Conduct applies within all community spaces, and also applies when
+an individual is officially representing the community in public spaces.
+Examples of representing our community include using an official e-mail address,
+posting via an official social media account, or acting as an appointed
+representative at an online or offline event.
+
+## Enforcement
+
+Instances of abusive, harassing, or otherwise unacceptable behavior may be
+reported to the community leaders responsible for enforcement at
+pychuang@gwu.edu or labarba@gwu.edu.
+All complaints will be reviewed and investigated promptly and fairly.
+
+All community leaders are obligated to respect the privacy and security of the
+reporter of any incident.
+
+## Enforcement Guidelines
+
+Community leaders will follow these Community Impact Guidelines in determining
+the consequences for any action they deem in violation of this Code of Conduct:
+
+### 1. Correction
+
+**Community Impact**: Use of inappropriate language or other behavior deemed
+unprofessional or unwelcome in the community.
+
+**Consequence**: A private, written warning from community leaders, providing
+clarity around the nature of the violation and an explanation of why the
+behavior was inappropriate. A public apology may be requested.
+
+### 2. Warning
+
+**Community Impact**: A violation through a single incident or series
+of actions.
+
+**Consequence**: A warning with consequences for continued behavior. No
+interaction with the people involved, including unsolicited interaction with
+those enforcing the Code of Conduct, for a specified period of time. This
+includes avoiding interactions in community spaces as well as external channels
+like social media. Violating these terms may lead to a temporary or
+permanent ban.
+
+### 3. Temporary Ban
+
+**Community Impact**: A serious violation of community standards, including
+sustained inappropriate behavior.
+
+**Consequence**: A temporary ban from any sort of interaction or public
+communication with the community for a specified period of time. No public or
+private interaction with the people involved, including unsolicited interaction
+with those enforcing the Code of Conduct, is allowed during this period.
+Violating these terms may lead to a permanent ban.
+
+### 4. Permanent Ban
+
+**Community Impact**: Demonstrating a pattern of violation of community
+standards, including sustained inappropriate behavior,  harassment of an
+individual, or aggression toward or disparagement of classes of individuals.
+
+**Consequence**: A permanent ban from any sort of public interaction within
+the community.
+
+## Attribution
+
+This Code of Conduct is adapted from the [Contributor Covenant][homepage],
+version 2.0, available at
+[https://www.contributor-covenant.org/version/2/0/code_of_conduct.html][v2.0].
+
+Community Impact Guidelines were inspired by 
+[Mozilla's code of conduct enforcement ladder][Mozilla CoC].
+
+For answers to common questions about this code of conduct, see the FAQ at
+[https://www.contributor-covenant.org/faq][FAQ]. Translations are available 
+at [https://www.contributor-covenant.org/translations][translations].
+
+[homepage]: https://www.contributor-covenant.org
+[v2.0]: https://www.contributor-covenant.org/version/2/0/code_of_conduct.html
+[Mozilla CoC]: https://github.com/mozilla/diversity
+[FAQ]: https://www.contributor-covenant.org/faq
+[translations]: https://www.contributor-covenant.org/translations
+
diff --git a/CONTRIBUTING.md b/CONTRIBUTING.md
new file mode 100644
index 0000000..d2a7878
--- /dev/null
+++ b/CONTRIBUTING.md
@@ -0,0 +1,49 @@
+Contributing
+============
+
+Thank you for considering contributing to this project. *geoclaw-landspil* is an
+open-source project hosted on GitHub at
+https://github.com/barbagroup/geoclaw-landspill.
+All contributions are welcome, including but not limited to:
+
+* bug reports,
+* bug fixes,
+* documentation improvement,
+* more examples or tests,
+* new features, and
+* performance enhancement.
+
+Don't hesitate to ask any questions! Questions help us know what can be improved
+and make this project more helpful to people.
+
+--------------------
+## How to contribute
+
+Please open an issue at our
+[GitHub repository](https://github.com/barbagroup/geoclaw-landspill) if there
+are any questions, bugs, or suggestions.
+
+Make a pull request if you want to add new features/modifications, fix bugs, or
+add more examples/tests. As we try to keep the project simple, all pull requests
+should make the `master` branch as the base branch. A pull request against the
+`master` branch triggers Travis CI and tests defined with GitHub Actions. Please
+make sure the newly added code passes all tests.
+
+------------------------------------------------------------------------
+## How to install *geoclaw-landspill* for development and run tests
+
+In addition to the standard installation methods described in
+[`deps_install_tests.md`](doc/deps_install_tests.md), developers should consider
+an editable installation. With an editable installation, code modifications take
+effect immediately without re-installation. Assuming currently under the
+top-level folder, do
+
+```
+$ pip install --editable .
+```
+
+You may need to uninstall previously installed *geoclaw-landspill* with
+`$ pip uninstall geoclaw-landspill`.
+
+To run tests, see the last section in
+[`deps_install_tests.md`](doc/deps_install_tests.md).
diff --git a/README.md b/README.md
index 07f2a43..7946a44 100644
--- a/README.md
+++ b/README.md
@@ -32,7 +32,9 @@ ruptures. These features include:
 ## Documentation
 1. [Dependencies, installation, and tests](doc/deps_install_tests.md)
 2. [Usage](doc/usage.md)
-3. [Containers: Docker and Singularity](doc/container.md)
+3. [Configuration file: `setrun.py`](doc/configuration.md)
+3. [Example cases](cases/README.md)
+4. [Containers: Docker and Singularity](doc/container.md)
 
 ------------------------------------------------------------------------
 ## Quick start
diff --git a/Singularityfiles/Singularity.v1.0.dev3 b/Singularityfiles/Singularity.v1.0.dev3
new file mode 100644
index 0000000..1d5f741
--- /dev/null
+++ b/Singularityfiles/Singularity.v1.0.dev3
@@ -0,0 +1,76 @@
+# multi-stage build: builder
+# ==========================
+Bootstrap: docker
+From: ubuntu:focal
+Stage: builder
+
+%environment
+    export MPLBACKEND=agg
+
+%post
+    # prepare environment
+    ln -snf /usr/share/zoneinfo/UTC /etc/localtime
+    echo UTC > /etc/timezone
+    apt update && apt -y full-upgrade
+    apt install -y --no-install-recommends git ca-certificates gfortran python3 python3-pip
+
+    # requires users to provide variable ${VER}
+    git clone --branch v1.0.dev3 --recurse-submodules https://github.com/barbagroup/geoclaw-landspill.git
+
+    # install dependencies
+    cd /geoclaw-landspill && pip3 install -r requirements-build.txt
+
+    # build binary wheel
+    cd /geoclaw-landspill && python3 setup.py bdist_wheel --build-type RELEASE
+
+
+# multi-stage build: production
+# =============================
+Bootstrap: docker
+From: ubuntu:focal
+Stage: production
+
+%labels
+    Author Pi-Yueh Chuang (pychuang@gwu.edu)
+    Version v1.0.dev3
+
+# copy data from builder
+%files from builder
+    /geoclaw-landspill/dist /root/dist
+    /geoclaw-landspill/cases /root/cases
+
+%post
+    # move example to a proper location
+    mv /root/cases /opt/geoclaw-landspill-cases 
+
+    # create runtime environment
+    ln -snf /usr/share/zoneinfo/UTC /etc/localtime && echo UTC > /etc/timezone
+    apt update && apt -y full-upgrade
+    apt install -y --no-install-recommends libgfortran5 libgomp1 python3 python3-pip
+    rm -rf /var/lib/apt/lists/*
+
+    # install geoclaw-landspill
+    echo $(find /root -name "*.whl")
+    pip3 install $(find /root -name "*.whl") && rm -rf /root/*
+
+%runscript
+    exec geoclaw-landspill $@
+
+%help
+
+    This help only covers the usage of this Singularity image. For the usage of
+    geoclaw-landspill package, see https://github.com/barbagroup/geoclaw-landspill
+
+    A Singularity image can be treated as an executable, which is equivalent to
+    the executable `geoclaw-landspill` from geoclaw-landspill package. Assuming
+    this Singularity image is called "landspill.sif", for example:
+
+        * To run a simulation: `$ landspill.sif run <path to case>`
+        * To create a NetCDF file: `$ landspill.sif createnc <path to case>`
+        * To create quick depth plots: `$ landspill.sif plotdepth <path to case>`
+
+    Alternatively, users can run these subcommands with, for example,
+
+        * To run a simulation: `$ singularity run landspill.sif run <path to case>`
+        * To create a NetCDF file: `$ singularity run landspill.sif createnc <path to case>`
+        * To create quick depth plots: `$ singularity run landspill.sif plotdepth <path to case>`
diff --git a/cases/README.md b/cases/README.md
new file mode 100644
index 0000000..998f1ed
--- /dev/null
+++ b/cases/README.md
@@ -0,0 +1,135 @@
+Examples
+========
+
+Currently, all examples are from the nearby region of Salt Lake City, Utah, USA.
+The naming of examples follows
+**{location}-{terrain type}-{oil type}-{additional setting}**.
+
+There are three files Under each example folder:
+
+* `setrun.py`: it serves as the configuration file of the simulation.
+* `roughness.txt`: all examples use a Darcy-Weisbach model that considers
+  surface roughness. This file provides roughness information to the model,
+  though the values are simply `0.1` everywhere.
+* `animation.gif`: this is an animation of the flow to give users a sense of how
+  their results should look.
+
+For the details of `setrun.py`, please refer to the documentation.
+
+To run an example, do
+
+```
+$ OMP_NUM_THREADS={number of CPU threads to use} geoclaw-landspill run {path to a case folder}
+```
+
+To create the `animation.gif`s shown below, do
+
+```
+$ geoclaw-landspill plotdepth --use-sat {path to a case folder}
+```
+
+This command only creates frames used by the animations. The frames are saved
+in the folder `_plots/sat/level02/` under each case. Finally, we use
+[`ffmpeg`](https://ffmpeg.org/) to combine frames into a GIF animation file:
+
+```
+$ ffmpeg \
+    -i {case folder}/_plots/sat/level02/frame00%03d.png \
+    -vf "fps=20,split[s0][s1];[s0]palettegen[p];[s1][p]paletteuse" \
+    {filename to save the resulting animation}
+```
+
+For a more serious analysis, users should consider creating a temporal NetCDF
+raster file for each case with
+
+```
+$ geoclaw-landspill createnc {path to a case folder}
+```
+
+The resulting NetCDF file will be at `_output/{case name}-depth-lvl02.nc`.
+
+----------
+## List of cases and results
+
+### 1. utah-flat-gasoline-no-evap
+
+Gasoline without evaporation on flat terrain.
+
+![utah-flat-gasoline-no-evap](./utah-flat-gasoline-no-evap/animation.gif)
+
+### 2. utah-flat-gasoline
+
+Gasoline with evaporation on flat terrain. File `_output/evaporated_fluid.dat`
+stores the total volume evaporated at the end of the simulation.
+
+![utah-flat-gasoline](./utah-flat-gasoline/animation.gif)
+
+### 3. utah-flat-maya-no-evap
+
+Maya crude oil without evaporation on flat terrain.
+
+![utah-flat-maya-no-evap](./utah-flat-maya-no-evap/animation.gif)
+
+### 4. utah-flat-maya
+
+Maya crude oil with evaporation on flat terrain. File `_output/evaporated_fluid.dat`
+stores the total volume evaporated at the end of the simulation.
+
+![utah-flat-maya](./utah-flat-maya/animation.gif)
+
+### 5. utah-hill-maya-no-hydro
+
+Maya crude oil on hilly terrain. Evaporation is on, but the in-land waterbody
+(a creek in this case) is turned off to see how the solver works with a
+drainage system.
+
+![utah-hill-maya-no-hydro](./utah-hill-maya-no-hydro/animation.gif)
+
+### 6. utah-hill-maya
+
+Maya crude oil on hilly terrain. Evaporation is on. The in-land waterbody
+(a creek in this case) is on. The creek catches all oil.
+`_output/removed_fluid.csv` records contact points and the fluid volumes flowing
+into waterbodies at the contact points.
+
+
+![utah-hill-maya](./utah-hill-maya/animation.gif)
+
+### 7. utah-structures-maya
+
+Maya crude oil on flat terrain but with some structures in the area. Both
+evaporation and in-land waterbodies are on.
+
+![utah-structures-maya](./utah-structures-maya/animation.gif)
+
+### 8. utah-waterbody-gasoline-no-evap
+
+Gasoline above land surrounded by in-land waterbodies.
+`_output/removed_fluid.csv` records contact points and the fluid volumes flowing
+into waterbodies at the contact points. The evaporation is off.
+
+![utah-waterbody-gasoline-no-evap](./utah-waterbody-gasoline-no-evap/animation.gif)
+
+### 9. utah-waterbody-gasoline
+
+Gasoline above land surrounded by in-land waterbodies.
+`_output/removed_fluid.csv` records contact points and the fluid volumes flowing
+into waterbodies at the contact points. The evaporation is on.
+
+![utah-waterbody-gasoline](./utah-waterbody-gasoline/animation.gif)
+
+### 10. utah-waterbody-maya-no-evap
+
+Maya crude oil above land surrounded by in-land waterbodies.
+`_output/removed_fluid.csv` records contact points and the fluid volumes flowing
+into waterbodies at the contact points. The evaporation is off.
+
+![utah-waterbody-maya-no-evap](./utah-waterbody-maya-no-evap/animation.gif)
+
+### 11. utah-waterbody-maya
+
+Maya crude oil above land surrounded by in-land waterbodies.
+`_output/removed_fluid.csv` records contact points and the fluid volumes flowing
+into waterbodies at the contact points. The evaporation is on.
+
+![utah-waterbody-maya](./utah-waterbody-maya/animation.gif)
diff --git a/cases/utah-flat-gasoline-no-evap/animation.gif b/cases/utah-flat-gasoline-no-evap/animation.gif
new file mode 100644
index 0000000..2910bc0
Binary files /dev/null and b/cases/utah-flat-gasoline-no-evap/animation.gif differ
diff --git a/cases/utah-flat-gasoline-no-evap/setrun.py b/cases/utah-flat-gasoline-no-evap/setrun.py
index f9d7903..00760e5 100644
--- a/cases/utah-flat-gasoline-no-evap/setrun.py
+++ b/cases/utah-flat-gasoline-no-evap/setrun.py
@@ -2,304 +2,63 @@
 
 The values set in the function setrun are then written out to data files
 that will be read in by the Fortran code.
-
 """
 # pylint: disable=no-member, too-many-statements
-import functools
-from gclandspill.clawutil import data
-from gclandspill.data import LandSpillData
+import gclandspill.data
 
 
-def setrun(claw_pkg: str = 'geoclaw'):
+def setrun():
     """Define the parameters used for running Clawpack.
 
-    Arguments
-    ---------
-    claw_pkg : str
-        Expected to be "geoclaw" for this setrun.
-
     Returns
     -------
-    rundata : clawutil.ClawRunData
+    rundata : gclandspill.data.ClawRunData
     """
 
-    assert claw_pkg.lower() == 'geoclaw',  "Expected claw_pkg = 'geoclaw'"
-
-    num_dim = 2
-    rundata = data.ClawRunData(claw_pkg, num_dim)
-
-    # GeoClaw specific parameters:
-    rundata = setgeo(rundata)
+    rundata = gclandspill.data.ClawRunData()
 
-    # AmrClaw specific parameters:
-    rundata = setamr(rundata)
+    # for convenience
+    break_point = [-12459650., 4986000.]
 
-    # an alias to object holding standard Clawpack parameters
-    clawdata = rundata.clawdata
-
-    # number of space dimensions:
-    clawdata.num_dim = num_dim
-
-    # lower and upper edge of computational domain:
-    clawdata.lower[0] = -12459650.-400.
-    clawdata.upper[0] = -12459650.+400.
-    clawdata.lower[1] = 4986000.-400
-    clawdata.upper[1] = 4986000.+400
+    # lower and upper edges of computational domain
+    rundata.clawdata.lower[0] = break_point[0]-400.
+    rundata.clawdata.upper[0] = break_point[0]+400.
+    rundata.clawdata.lower[1] = break_point[1]-400
+    rundata.clawdata.upper[1] = break_point[1]+400
 
     # number of grid cells: coarsest grid
-    clawdata.num_cells[0] = 200
-    clawdata.num_cells[1] = 200
-
-    # number of equations in the system: continuity eq. and x, y momemtums
-    clawdata.num_eqn = 3
-
-    # number of auxiliary variables: elevation and hydrological data
-    clawdata.num_aux = 2
-
-    # index of aux array corresponding to capacity function, if there is one:
-    clawdata.capa_index = 0
-
-    # initial time:
-    clawdata.t0 = 0.0
-
-    # restart from checkpoint file of a previous run?
-    clawdata.restart = False               # True to restart from prior results
-    clawdata.restart_file = 'fort.chk00006'  # File to use for restart data
-
-    # temporal data output style
-    clawdata.output_style = 1
-
-    if clawdata.output_style == 1:  # equally spaced times up to tfinal
-        clawdata.num_output_times = 240
-        clawdata.tfinal = 28800
-        clawdata.output_t0 = True  # output at initial (or restart) time?
-    elif clawdata.output_style == 2:  # a list of output times
-        clawdata.output_times = [0., 1800.]
-    elif clawdata.output_style == 3:  # every iout timesteps with a total of ntot time steps:
-        clawdata.output_step_interval = 80
-        clawdata.total_steps = 4000
-        clawdata.output_t0 = True
-
-    # output data format
-    clawdata.output_format = 'binary'
-
-    # which veriable to output
-    clawdata.output_q_components = 'all'  # could be list such as [True,True]
-
-    # which aux data to output
-    clawdata.output_aux_components = 'all'  # could be list
-
-    # whether to output aux data only at t0
-    clawdata.output_aux_onlyonce = False
-
-    # print information up to this AMR level
-    clawdata.verbosity = 3
-
-    # use veriable time step size based on CFL
-    clawdata.dt_variable = 1
-
-    # data required to calculate a proper initial time step size
-    cell_area = \
-        (clawdata.upper[0] - clawdata.lower[0]) * (clawdata.upper[1] - clawdata.lower[1]) / \
-        (clawdata.num_cells[0] * clawdata.num_cells[1])
-    cell_area /= functools.reduce(
-        lambda x, y: x*y, rundata.amrdata.refinement_ratios_x[:rundata.amrdata.amr_levels_max], 1)
-    cell_area /= functools.reduce(
-        lambda x, y: x*y, rundata.amrdata.refinement_ratios_y[:rundata.amrdata.amr_levels_max], 1)
-    vrate = rundata.landspill_data.point_sources.point_sources[0][-1][0]
-
-    # initial time step for variable dt.
-    clawdata.dt_initial = 0.3 * cell_area / vrate
-
-    # max time step allowed if variable dt
-    clawdata.dt_max = 4.0
-
-    # desired CFL if variable dt, and max CFL allowed without retaking a step with a smaller dt
-    clawdata.cfl_desired = 0.9
-    clawdata.cfl_max = 0.95
-
-    # maximum number of time steps allowed between two output times:
-    clawdata.steps_max = 100000
-
-    # order of accuracy:  1 => Godunov,  2 => Lax-Wendroff plus limiters
-    clawdata.order = 2
-
-    # use dimensional splitting? (not yet available for AMR)
-    clawdata.dimensional_split = 'unsplit'
-
-    # 2 or 'all' => corner transport of 2nd order corrections too
-    clawdata.transverse_waves = 2
-
-    # number of waves in the Riemann solution:
-    clawdata.num_waves = 3
-
-    # list of limiters to use for each wave family:
-    clawdata.limiter = ['mc', 'mc', 'mc']
-
-    # use f-wave version of algorithms
-    clawdata.use_fwaves = True
-
-    # src_split == 1 or 'godunov' ==> Godunov (1st order) splitting used,
-    clawdata.source_split = 'godunov'
-
-    # number of ghost cells (usually 2)
-    clawdata.num_ghost = 2
+    rundata.clawdata.num_cells[0] = 200
+    rundata.clawdata.num_cells[1] = 200
 
-    # BC: 1 => extrapolation (non-reflecting outflow)
-    clawdata.bc_lower[0] = 1
-    clawdata.bc_upper[0] = 1
-    clawdata.bc_lower[1] = 1
-    clawdata.bc_upper[1] = 1
+    # temporal data output
+    rundata.clawdata.output_style = 1
+    rundata.clawdata.num_output_times = 240
+    rundata.clawdata.tfinal = 28800
+    rundata.clawdata.output_t0 = True  # output at initial (or restart) time?
 
-    # checkout mechanism selection
-    clawdata.checkpt_style = 0
-
-    if clawdata.checkpt_style == 0:  # do not checkpoint at all
-        pass
-    elif clawdata.checkpt_style == 1:  # checkpoint only at tfinal.
-        pass
-    elif clawdata.checkpt_style == 2:  # specify a list of checkpoint times.
-        clawdata.checkpt_times = [0.1, 0.15]
-    elif clawdata.checkpt_style == 3:  # checkpoint every checkpt_interval timesteps (on Level 1)
-        clawdata.checkpt_interval = 5
-
-    return rundata
-
-
-def setamr(rundata: data.ClawRunData):
-    """Set AmrClaw specific runtime parameters."""
-
-    try:
-        amrdata = rundata.amrdata
-    except AttributeError as err:
-        raise AttributeError("This rundata has no amrdata attribute") from err
-
-    # max number of refinement levels:
-    amrdata.amr_levels_max = 2
-
-    # List of refinement ratios at each level (length at least mxnest-1)
-    amrdata.refinement_ratios_x = [4]
-    amrdata.refinement_ratios_y = [4]
-    amrdata.refinement_ratios_t = [4]
-
-    # aux variables are defined at cell centers
-    amrdata.aux_type = ['center', 'center']
-
-    # use flag2refine routine to flag refinement (rather than using richardson error)
-    amrdata.flag_richardson = False
-    amrdata.flag2refine = True
-
-    # steps to take on each level L between regriddings of level L+1:
-    amrdata.regrid_interval = 1
-
-    # width of buffer zone around flagged points
-    amrdata.regrid_buffer_width = 1
-
-    # clustering alg. cutoff for (# flagged pts) / (total # of cells refined)
-    amrdata.clustering_cutoff = 0.80000
-
-    # print info about each regridding up to this level:
-    amrdata.verbosity_regrid = 0
-
-    #  ----- For developers -----
-    # Toggle debugging print statements:
-    amrdata.dprint = False      # print domain flags
-    amrdata.eprint = False      # print err est flags
-    amrdata.edebug = False      # even more err est flags
-    amrdata.gprint = False      # grid bisection/clustering
-    amrdata.nprint = False      # proper nesting output
-    amrdata.pprint = False      # proj. of tagged points
-    amrdata.rprint = False      # print regridding summary
-    amrdata.sprint = False      # space/memory output
-    amrdata.tprint = False      # time step reporting each level
-    amrdata.uprint = False      # update/upbnd reporting
-
-    # customize regions of refinement
-    # rundata.regiondata.regions.append([minlevel, maxlevel, t1, t2, x1, x2, y1, y2])
-
-    # gauges
-    # rundata.gaugedata.gauges.append([gaugeno, x, y, t1, t2])
-
-    return rundata
-
-
-def setgeo(rundata: data.ClawRunData):
-    """Set GeoClaw specific runtime parameters."""
-
-    try:
-        geo_data = rundata.geo_data
-    except AttributeError as err:
-        raise AttributeError("This rundata has no geo_data attribute") from err
-
-    # == Physics ==
-    geo_data.gravity = 9.81
-    geo_data.coordinate_system = 1
-    geo_data.earth_radius = 6367.5e3
-
-    # == Forcing Options
-    geo_data.coriolis_forcing = False
-
-    # == Algorithm and Initial Conditions ==
-    geo_data.sea_level = 0.0
-    geo_data.dry_tolerance = 1.e-4
-    geo_data.friction_forcing = False
-    geo_data.manning_coefficient = 0.035
-    geo_data.friction_depth = 1.e6
-
-    # Refinement data
-    refinement_data = rundata.refinement_data
-    refinement_data.wave_tolerance = 1.e-5
-    refinement_data.speed_tolerance = [1e-8]
-    refinement_data.variable_dt_refinement_ratios = True
-
-    # for topography, append lines of the form [topotype, fname]
-    topo_data = rundata.topo_data
-    topo_data.topofiles.append([3, '../common-files/utah-flat.asc'])
-
-    # for moving topography, append lines of the form: [topotype, fname]
-    # dtopo_data = rundata.dtopo_data
-
-    # for qinit perturbations, append lines of the form: [fname]
-    rundata.qinit_data.qinit_type = 0
-    rundata.qinit_data.qinitfiles = []
-
-    # for fixed grids append lines of the form
-    # [t1, t2, noutput, x1, x2, y1, y2, xpoints, ypoints, ioutarrivaltimes, ioutsurfacemax]
-    # fixedgrids = rundata.fixed_grid_data
+    # add a topography file; format [topotype, fname]
+    rundata.topo_data.topofiles.append([3, '../common-files/utah-flat.asc'])
 
     # Land-spill module settings
-    rundata.add_data(LandSpillData(), 'landspill_data')
-    landspill = rundata.landspill_data
-    landspill.ref_mu = 0.6512  # cP @ 15 degree C
-    landspill.ref_temperature = 15.
-    landspill.ambient_temperature = 25.
-    landspill.density = 800.  # kg / m^3 @ 15 degree C; will overwrite rho in GeoClaw
-
-    # extra parameters
-    landspill.update_tol = geo_data.dry_tolerance
-    landspill.refine_tol = 0.0
+    rundata.landspill_data.ref_mu = 0.6512  # cP @ 15 degree C
+    rundata.landspill_data.ref_temperature = 15.
+    rundata.landspill_data.ambient_temperature = 25.
+    rundata.landspill_data.density = 800.  # kg / m^3 @ 15 degree C; will overwrite rho in GeoClaw
 
     # Point sources
-    ptsources_data = landspill.point_sources
-    ptsources_data.n_point_sources = 1
-    ptsources_data.point_sources.append(
-        [[-12459650., 4986000.], 2, [1800., 12600.], [0.5, 0.1]])
+    rundata.landspill_data.point_sources.n_point_sources = 1
+    rundata.landspill_data.point_sources.point_sources.append([break_point, 2, [1800., 12600.], [0.5, 0.1]])
 
     # Darcy-Weisbach friction
-    darcy_weisbach_data = landspill.darcy_weisbach_friction
-    darcy_weisbach_data.type = 4
-    darcy_weisbach_data.dry_tol = 1e-4
-    darcy_weisbach_data.friction_tol = 1e6
-    darcy_weisbach_data.default_roughness = 0.0
-    darcy_weisbach_data.filename = "roughness.txt"
+    rundata.landspill_data.darcy_weisbach_friction.type = 4
+    rundata.landspill_data.darcy_weisbach_friction.default_roughness = 0.0
+    rundata.landspill_data.darcy_weisbach_friction.filename = "roughness.txt"
 
-    # hydrological features
-    hydro_feature_data = landspill.hydro_features
-    hydro_feature_data.files.append("../common-files/utah-flat-hydro.asc")
+    # add a hydrological feature raster file
+    rundata.landspill_data.hydro_features.files.append("../common-files/utah-flat-hydro.asc")
 
     # Evaporation
-    evaporation_data = landspill.evaporation
-    evaporation_data.type = 0
-    evaporation_data.coefficients = [13.2, 0.21]
+    rundata.landspill_data.evaporation.type = 0
+    rundata.landspill_data.evaporation.coefficients = [13.2, 0.21]
 
     return rundata
diff --git a/cases/utah-flat-gasoline/animation.gif b/cases/utah-flat-gasoline/animation.gif
new file mode 100644
index 0000000..ec4e09b
Binary files /dev/null and b/cases/utah-flat-gasoline/animation.gif differ
diff --git a/cases/utah-flat-gasoline/setrun.py b/cases/utah-flat-gasoline/setrun.py
index 39870c5..ef921e9 100644
--- a/cases/utah-flat-gasoline/setrun.py
+++ b/cases/utah-flat-gasoline/setrun.py
@@ -2,304 +2,63 @@
 
 The values set in the function setrun are then written out to data files
 that will be read in by the Fortran code.
-
 """
 # pylint: disable=no-member, too-many-statements
-import functools
-from gclandspill.clawutil import data
-from gclandspill.data import LandSpillData
+import gclandspill.data
 
 
-def setrun(claw_pkg: str = 'geoclaw'):
+def setrun():
     """Define the parameters used for running Clawpack.
 
-    Arguments
-    ---------
-    claw_pkg : str
-        Expected to be "geoclaw" for this setrun.
-
     Returns
     -------
-    rundata : clawutil.ClawRunData
+    rundata : gclandspill.data.ClawRunData
     """
 
-    assert claw_pkg.lower() == 'geoclaw',  "Expected claw_pkg = 'geoclaw'"
-
-    num_dim = 2
-    rundata = data.ClawRunData(claw_pkg, num_dim)
-
-    # GeoClaw specific parameters:
-    rundata = setgeo(rundata)
+    rundata = gclandspill.data.ClawRunData()
 
-    # AmrClaw specific parameters:
-    rundata = setamr(rundata)
+    # for convenience
+    break_point = [-12459650., 4986000.]
 
-    # an alias to object holding standard Clawpack parameters
-    clawdata = rundata.clawdata
-
-    # number of space dimensions:
-    clawdata.num_dim = num_dim
-
-    # lower and upper edge of computational domain:
-    clawdata.lower[0] = -12459650.-400.
-    clawdata.upper[0] = -12459650.+400.
-    clawdata.lower[1] = 4986000.-400
-    clawdata.upper[1] = 4986000.+400
+    # lower and upper edges of computational domain
+    rundata.clawdata.lower[0] = break_point[0]-400.
+    rundata.clawdata.upper[0] = break_point[0]+400.
+    rundata.clawdata.lower[1] = break_point[1]-400
+    rundata.clawdata.upper[1] = break_point[1]+400
 
     # number of grid cells: coarsest grid
-    clawdata.num_cells[0] = 200
-    clawdata.num_cells[1] = 200
-
-    # number of equations in the system: continuity eq. and x, y momemtums
-    clawdata.num_eqn = 3
-
-    # number of auxiliary variables: elevation and hydrological data
-    clawdata.num_aux = 2
-
-    # index of aux array corresponding to capacity function, if there is one:
-    clawdata.capa_index = 0
-
-    # initial time:
-    clawdata.t0 = 0.0
-
-    # restart from checkpoint file of a previous run?
-    clawdata.restart = False               # True to restart from prior results
-    clawdata.restart_file = 'fort.chk00006'  # File to use for restart data
-
-    # temporal data output style
-    clawdata.output_style = 1
-
-    if clawdata.output_style == 1:  # equally spaced times up to tfinal
-        clawdata.num_output_times = 240
-        clawdata.tfinal = 28800
-        clawdata.output_t0 = True  # output at initial (or restart) time?
-    elif clawdata.output_style == 2:  # a list of output times
-        clawdata.output_times = [0., 1800.]
-    elif clawdata.output_style == 3:  # every iout timesteps with a total of ntot time steps:
-        clawdata.output_step_interval = 80
-        clawdata.total_steps = 4000
-        clawdata.output_t0 = True
-
-    # output data format
-    clawdata.output_format = 'binary'
-
-    # which veriable to output
-    clawdata.output_q_components = 'all'  # could be list such as [True,True]
-
-    # which aux data to output
-    clawdata.output_aux_components = 'all'  # could be list
-
-    # whether to output aux data only at t0
-    clawdata.output_aux_onlyonce = False
-
-    # print information up to this AMR level
-    clawdata.verbosity = 3
-
-    # use veriable time step size based on CFL
-    clawdata.dt_variable = 1
-
-    # data required to calculate a proper initial time step size
-    cell_area = \
-        (clawdata.upper[0] - clawdata.lower[0]) * (clawdata.upper[1] - clawdata.lower[1]) / \
-        (clawdata.num_cells[0] * clawdata.num_cells[1])
-    cell_area /= functools.reduce(
-        lambda x, y: x*y, rundata.amrdata.refinement_ratios_x[:rundata.amrdata.amr_levels_max], 1)
-    cell_area /= functools.reduce(
-        lambda x, y: x*y, rundata.amrdata.refinement_ratios_y[:rundata.amrdata.amr_levels_max], 1)
-    vrate = rundata.landspill_data.point_sources.point_sources[0][-1][0]
-
-    # initial time step for variable dt.
-    clawdata.dt_initial = 0.3 * cell_area / vrate
-
-    # max time step allowed if variable dt
-    clawdata.dt_max = 4.0
-
-    # desired CFL if variable dt, and max CFL allowed without retaking a step with a smaller dt
-    clawdata.cfl_desired = 0.9
-    clawdata.cfl_max = 0.95
-
-    # maximum number of time steps allowed between two output times:
-    clawdata.steps_max = 100000
-
-    # order of accuracy:  1 => Godunov,  2 => Lax-Wendroff plus limiters
-    clawdata.order = 2
-
-    # use dimensional splitting? (not yet available for AMR)
-    clawdata.dimensional_split = 'unsplit'
-
-    # 2 or 'all' => corner transport of 2nd order corrections too
-    clawdata.transverse_waves = 2
-
-    # number of waves in the Riemann solution:
-    clawdata.num_waves = 3
-
-    # list of limiters to use for each wave family:
-    clawdata.limiter = ['mc', 'mc', 'mc']
-
-    # use f-wave version of algorithms
-    clawdata.use_fwaves = True
-
-    # src_split == 1 or 'godunov' ==> Godunov (1st order) splitting used,
-    clawdata.source_split = 'godunov'
-
-    # number of ghost cells (usually 2)
-    clawdata.num_ghost = 2
+    rundata.clawdata.num_cells[0] = 200
+    rundata.clawdata.num_cells[1] = 200
 
-    # BC: 1 => extrapolation (non-reflecting outflow)
-    clawdata.bc_lower[0] = 1
-    clawdata.bc_upper[0] = 1
-    clawdata.bc_lower[1] = 1
-    clawdata.bc_upper[1] = 1
+    # temporal data output
+    rundata.clawdata.output_style = 1
+    rundata.clawdata.num_output_times = 240
+    rundata.clawdata.tfinal = 28800
+    rundata.clawdata.output_t0 = True  # output at initial (or restart) time?
 
-    # checkout mechanism selection
-    clawdata.checkpt_style = 0
-
-    if clawdata.checkpt_style == 0:  # do not checkpoint at all
-        pass
-    elif clawdata.checkpt_style == 1:  # checkpoint only at tfinal.
-        pass
-    elif clawdata.checkpt_style == 2:  # specify a list of checkpoint times.
-        clawdata.checkpt_times = [0.1, 0.15]
-    elif clawdata.checkpt_style == 3:  # checkpoint every checkpt_interval timesteps (on Level 1)
-        clawdata.checkpt_interval = 5
-
-    return rundata
-
-
-def setamr(rundata: data.ClawRunData):
-    """Set AmrClaw specific runtime parameters."""
-
-    try:
-        amrdata = rundata.amrdata
-    except AttributeError as err:
-        raise AttributeError("This rundata has no amrdata attribute") from err
-
-    # max number of refinement levels:
-    amrdata.amr_levels_max = 2
-
-    # List of refinement ratios at each level (length at least mxnest-1)
-    amrdata.refinement_ratios_x = [4]
-    amrdata.refinement_ratios_y = [4]
-    amrdata.refinement_ratios_t = [4]
-
-    # aux variables are defined at cell centers
-    amrdata.aux_type = ['center', 'center']
-
-    # use flag2refine routine to flag refinement (rather than using richardson error)
-    amrdata.flag_richardson = False
-    amrdata.flag2refine = True
-
-    # steps to take on each level L between regriddings of level L+1:
-    amrdata.regrid_interval = 1
-
-    # width of buffer zone around flagged points
-    amrdata.regrid_buffer_width = 1
-
-    # clustering alg. cutoff for (# flagged pts) / (total # of cells refined)
-    amrdata.clustering_cutoff = 0.80000
-
-    # print info about each regridding up to this level:
-    amrdata.verbosity_regrid = 0
-
-    #  ----- For developers -----
-    # Toggle debugging print statements:
-    amrdata.dprint = False      # print domain flags
-    amrdata.eprint = False      # print err est flags
-    amrdata.edebug = False      # even more err est flags
-    amrdata.gprint = False      # grid bisection/clustering
-    amrdata.nprint = False      # proper nesting output
-    amrdata.pprint = False      # proj. of tagged points
-    amrdata.rprint = False      # print regridding summary
-    amrdata.sprint = False      # space/memory output
-    amrdata.tprint = False      # time step reporting each level
-    amrdata.uprint = False      # update/upbnd reporting
-
-    # customize regions of refinement
-    # rundata.regiondata.regions.append([minlevel, maxlevel, t1, t2, x1, x2, y1, y2])
-
-    # gauges
-    # rundata.gaugedata.gauges.append([gaugeno, x, y, t1, t2])
-
-    return rundata
-
-
-def setgeo(rundata: data.ClawRunData):
-    """Set GeoClaw specific runtime parameters."""
-
-    try:
-        geo_data = rundata.geo_data
-    except AttributeError as err:
-        raise AttributeError("This rundata has no geo_data attribute") from err
-
-    # == Physics ==
-    geo_data.gravity = 9.81
-    geo_data.coordinate_system = 1
-    geo_data.earth_radius = 6367.5e3
-
-    # == Forcing Options
-    geo_data.coriolis_forcing = False
-
-    # == Algorithm and Initial Conditions ==
-    geo_data.sea_level = 0.0
-    geo_data.dry_tolerance = 1.e-4
-    geo_data.friction_forcing = False
-    geo_data.manning_coefficient = 0.035
-    geo_data.friction_depth = 1.e6
-
-    # Refinement data
-    refinement_data = rundata.refinement_data
-    refinement_data.wave_tolerance = 1.e-5
-    refinement_data.speed_tolerance = [1e-8]
-    refinement_data.variable_dt_refinement_ratios = True
-
-    # for topography, append lines of the form [topotype, fname]
-    topo_data = rundata.topo_data
-    topo_data.topofiles.append([3, '../common-files/utah-flat.asc'])
-
-    # for moving topography, append lines of the form: [topotype, fname]
-    # dtopo_data = rundata.dtopo_data
-
-    # for qinit perturbations, append lines of the form: [fname]
-    rundata.qinit_data.qinit_type = 0
-    rundata.qinit_data.qinitfiles = []
-
-    # for fixed grids append lines of the form
-    # [t1, t2, noutput, x1, x2, y1, y2, xpoints, ypoints, ioutarrivaltimes, ioutsurfacemax]
-    # fixedgrids = rundata.fixed_grid_data
+    # add a topography file; format [topotype, fname]
+    rundata.topo_data.topofiles.append([3, '../common-files/utah-flat.asc'])
 
     # Land-spill module settings
-    rundata.add_data(LandSpillData(), 'landspill_data')
-    landspill = rundata.landspill_data
-    landspill.ref_mu = 0.6512  # cP @ 15 degree C
-    landspill.ref_temperature = 15.
-    landspill.ambient_temperature = 25.
-    landspill.density = 800.  # kg / m^3 @ 15 degree C; will overwrite rho in GeoClaw
-
-    # extra parameters
-    landspill.update_tol = geo_data.dry_tolerance
-    landspill.refine_tol = 0.0
+    rundata.landspill_data.ref_mu = 0.6512  # cP @ 15 degree C
+    rundata.landspill_data.ref_temperature = 15.
+    rundata.landspill_data.ambient_temperature = 25.
+    rundata.landspill_data.density = 800.  # kg / m^3 @ 15 degree C; will overwrite rho in GeoClaw
 
     # Point sources
-    ptsources_data = landspill.point_sources
-    ptsources_data.n_point_sources = 1
-    ptsources_data.point_sources.append(
-        [[-12459650., 4986000.], 2, [1800., 12600.], [0.5, 0.1]])
+    rundata.landspill_data.point_sources.n_point_sources = 1
+    rundata.landspill_data.point_sources.point_sources.append([break_point, 2, [1800., 12600.], [0.5, 0.1]])
 
     # Darcy-Weisbach friction
-    darcy_weisbach_data = landspill.darcy_weisbach_friction
-    darcy_weisbach_data.type = 4
-    darcy_weisbach_data.dry_tol = 1e-4
-    darcy_weisbach_data.friction_tol = 1e6
-    darcy_weisbach_data.default_roughness = 0.0
-    darcy_weisbach_data.filename = "roughness.txt"
+    rundata.landspill_data.darcy_weisbach_friction.type = 4
+    rundata.landspill_data.darcy_weisbach_friction.default_roughness = 0.0
+    rundata.landspill_data.darcy_weisbach_friction.filename = "roughness.txt"
 
-    # hydrological features
-    hydro_feature_data = landspill.hydro_features
-    hydro_feature_data.files.append("../common-files/utah-flat-hydro.asc")
+    # add a hydrological feature raster file
+    rundata.landspill_data.hydro_features.files.append("../common-files/utah-flat-hydro.asc")
 
     # Evaporation
-    evaporation_data = landspill.evaporation
-    evaporation_data.type = 1
-    evaporation_data.coefficients = [13.2, 0.21]
+    rundata.landspill_data.evaporation.type = 1
+    rundata.landspill_data.evaporation.coefficients = [13.2, 0.21]
 
     return rundata
diff --git a/cases/utah-flat-maya-no-evap/animation.gif b/cases/utah-flat-maya-no-evap/animation.gif
new file mode 100644
index 0000000..864b4ed
Binary files /dev/null and b/cases/utah-flat-maya-no-evap/animation.gif differ
diff --git a/cases/utah-flat-maya-no-evap/setrun.py b/cases/utah-flat-maya-no-evap/setrun.py
index 2ba6315..c026dbb 100644
--- a/cases/utah-flat-maya-no-evap/setrun.py
+++ b/cases/utah-flat-maya-no-evap/setrun.py
@@ -2,304 +2,63 @@
 
 The values set in the function setrun are then written out to data files
 that will be read in by the Fortran code.
-
 """
 # pylint: disable=no-member, too-many-statements
-import functools
-from gclandspill.clawutil import data
-from gclandspill.data import LandSpillData
+import gclandspill.data
 
 
-def setrun(claw_pkg: str = 'geoclaw'):
+def setrun():
     """Define the parameters used for running Clawpack.
 
-    Arguments
-    ---------
-    claw_pkg : str
-        Expected to be "geoclaw" for this setrun.
-
     Returns
     -------
-    rundata : clawutil.ClawRunData
+    rundata : gclandspill.data.ClawRunData
     """
 
-    assert claw_pkg.lower() == 'geoclaw',  "Expected claw_pkg = 'geoclaw'"
-
-    num_dim = 2
-    rundata = data.ClawRunData(claw_pkg, num_dim)
-
-    # GeoClaw specific parameters:
-    rundata = setgeo(rundata)
+    rundata = gclandspill.data.ClawRunData()
 
-    # AmrClaw specific parameters:
-    rundata = setamr(rundata)
+    # for convenience
+    break_point = [-12459650., 4986000.]
 
-    # an alias to object holding standard Clawpack parameters
-    clawdata = rundata.clawdata
-
-    # number of space dimensions:
-    clawdata.num_dim = num_dim
-
-    # lower and upper edge of computational domain:
-    clawdata.lower[0] = -12459650.-400.
-    clawdata.upper[0] = -12459650.+400.
-    clawdata.lower[1] = 4986000.-400
-    clawdata.upper[1] = 4986000.+400
+    # lower and upper edges of computational domain
+    rundata.clawdata.lower[0] = break_point[0]-400.
+    rundata.clawdata.upper[0] = break_point[0]+400.
+    rundata.clawdata.lower[1] = break_point[1]-400
+    rundata.clawdata.upper[1] = break_point[1]+400
 
     # number of grid cells: coarsest grid
-    clawdata.num_cells[0] = 200
-    clawdata.num_cells[1] = 200
-
-    # number of equations in the system: continuity eq. and x, y momemtums
-    clawdata.num_eqn = 3
-
-    # number of auxiliary variables: elevation and hydrological data
-    clawdata.num_aux = 2
-
-    # index of aux array corresponding to capacity function, if there is one:
-    clawdata.capa_index = 0
-
-    # initial time:
-    clawdata.t0 = 0.0
-
-    # restart from checkpoint file of a previous run?
-    clawdata.restart = False               # True to restart from prior results
-    clawdata.restart_file = 'fort.chk00006'  # File to use for restart data
-
-    # temporal data output style
-    clawdata.output_style = 1
-
-    if clawdata.output_style == 1:  # equally spaced times up to tfinal
-        clawdata.num_output_times = 240
-        clawdata.tfinal = 28800
-        clawdata.output_t0 = True  # output at initial (or restart) time?
-    elif clawdata.output_style == 2:  # a list of output times
-        clawdata.output_times = [0., 1800.]
-    elif clawdata.output_style == 3:  # every iout timesteps with a total of ntot time steps:
-        clawdata.output_step_interval = 80
-        clawdata.total_steps = 4000
-        clawdata.output_t0 = True
-
-    # output data format
-    clawdata.output_format = 'binary'
-
-    # which veriable to output
-    clawdata.output_q_components = 'all'  # could be list such as [True,True]
-
-    # which aux data to output
-    clawdata.output_aux_components = 'all'  # could be list
-
-    # whether to output aux data only at t0
-    clawdata.output_aux_onlyonce = False
-
-    # print information up to this AMR level
-    clawdata.verbosity = 3
-
-    # use veriable time step size based on CFL
-    clawdata.dt_variable = 1
-
-    # data required to calculate a proper initial time step size
-    cell_area = \
-        (clawdata.upper[0] - clawdata.lower[0]) * (clawdata.upper[1] - clawdata.lower[1]) / \
-        (clawdata.num_cells[0] * clawdata.num_cells[1])
-    cell_area /= functools.reduce(
-        lambda x, y: x*y, rundata.amrdata.refinement_ratios_x[:rundata.amrdata.amr_levels_max], 1)
-    cell_area /= functools.reduce(
-        lambda x, y: x*y, rundata.amrdata.refinement_ratios_y[:rundata.amrdata.amr_levels_max], 1)
-    vrate = rundata.landspill_data.point_sources.point_sources[0][-1][0]
-
-    # initial time step for variable dt.
-    clawdata.dt_initial = 0.3 * cell_area / vrate
-
-    # max time step allowed if variable dt
-    clawdata.dt_max = 4.0
-
-    # desired CFL if variable dt, and max CFL allowed without retaking a step with a smaller dt
-    clawdata.cfl_desired = 0.9
-    clawdata.cfl_max = 0.95
-
-    # maximum number of time steps allowed between two output times:
-    clawdata.steps_max = 100000
-
-    # order of accuracy:  1 => Godunov,  2 => Lax-Wendroff plus limiters
-    clawdata.order = 2
-
-    # use dimensional splitting? (not yet available for AMR)
-    clawdata.dimensional_split = 'unsplit'
-
-    # 2 or 'all' => corner transport of 2nd order corrections too
-    clawdata.transverse_waves = 2
-
-    # number of waves in the Riemann solution:
-    clawdata.num_waves = 3
-
-    # list of limiters to use for each wave family:
-    clawdata.limiter = ['mc', 'mc', 'mc']
-
-    # use f-wave version of algorithms
-    clawdata.use_fwaves = True
-
-    # src_split == 1 or 'godunov' ==> Godunov (1st order) splitting used,
-    clawdata.source_split = 'godunov'
-
-    # number of ghost cells (usually 2)
-    clawdata.num_ghost = 2
+    rundata.clawdata.num_cells[0] = 200
+    rundata.clawdata.num_cells[1] = 200
 
-    # BC: 1 => extrapolation (non-reflecting outflow)
-    clawdata.bc_lower[0] = 1
-    clawdata.bc_upper[0] = 1
-    clawdata.bc_lower[1] = 1
-    clawdata.bc_upper[1] = 1
+    # temporal data output
+    rundata.clawdata.output_style = 1
+    rundata.clawdata.num_output_times = 240
+    rundata.clawdata.tfinal = 28800
+    rundata.clawdata.output_t0 = True  # output at initial (or restart) time?
 
-    # checkout mechanism selection
-    clawdata.checkpt_style = 0
-
-    if clawdata.checkpt_style == 0:  # do not checkpoint at all
-        pass
-    elif clawdata.checkpt_style == 1:  # checkpoint only at tfinal.
-        pass
-    elif clawdata.checkpt_style == 2:  # specify a list of checkpoint times.
-        clawdata.checkpt_times = [0.1, 0.15]
-    elif clawdata.checkpt_style == 3:  # checkpoint every checkpt_interval timesteps (on Level 1)
-        clawdata.checkpt_interval = 5
-
-    return rundata
-
-
-def setamr(rundata: data.ClawRunData):
-    """Set AmrClaw specific runtime parameters."""
-
-    try:
-        amrdata = rundata.amrdata
-    except AttributeError as err:
-        raise AttributeError("This rundata has no amrdata attribute") from err
-
-    # max number of refinement levels:
-    amrdata.amr_levels_max = 2
-
-    # List of refinement ratios at each level (length at least mxnest-1)
-    amrdata.refinement_ratios_x = [4]
-    amrdata.refinement_ratios_y = [4]
-    amrdata.refinement_ratios_t = [4]
-
-    # aux variables are defined at cell centers
-    amrdata.aux_type = ['center', 'center']
-
-    # use flag2refine routine to flag refinement (rather than using richardson error)
-    amrdata.flag_richardson = False
-    amrdata.flag2refine = True
-
-    # steps to take on each level L between regriddings of level L+1:
-    amrdata.regrid_interval = 1
-
-    # width of buffer zone around flagged points
-    amrdata.regrid_buffer_width = 1
-
-    # clustering alg. cutoff for (# flagged pts) / (total # of cells refined)
-    amrdata.clustering_cutoff = 0.80000
-
-    # print info about each regridding up to this level:
-    amrdata.verbosity_regrid = 0
-
-    #  ----- For developers -----
-    # Toggle debugging print statements:
-    amrdata.dprint = False      # print domain flags
-    amrdata.eprint = False      # print err est flags
-    amrdata.edebug = False      # even more err est flags
-    amrdata.gprint = False      # grid bisection/clustering
-    amrdata.nprint = False      # proper nesting output
-    amrdata.pprint = False      # proj. of tagged points
-    amrdata.rprint = False      # print regridding summary
-    amrdata.sprint = False      # space/memory output
-    amrdata.tprint = False      # time step reporting each level
-    amrdata.uprint = False      # update/upbnd reporting
-
-    # customize regions of refinement
-    # rundata.regiondata.regions.append([minlevel, maxlevel, t1, t2, x1, x2, y1, y2])
-
-    # gauges
-    # rundata.gaugedata.gauges.append([gaugeno, x, y, t1, t2])
-
-    return rundata
-
-
-def setgeo(rundata: data.ClawRunData):
-    """Set GeoClaw specific runtime parameters."""
-
-    try:
-        geo_data = rundata.geo_data
-    except AttributeError as err:
-        raise AttributeError("This rundata has no geo_data attribute") from err
-
-    # == Physics ==
-    geo_data.gravity = 9.81
-    geo_data.coordinate_system = 1
-    geo_data.earth_radius = 6367.5e3
-
-    # == Forcing Options
-    geo_data.coriolis_forcing = False
-
-    # == Algorithm and Initial Conditions ==
-    geo_data.sea_level = 0.0
-    geo_data.dry_tolerance = 1.e-4
-    geo_data.friction_forcing = False
-    geo_data.manning_coefficient = 0.035
-    geo_data.friction_depth = 1.e6
-
-    # Refinement data
-    refinement_data = rundata.refinement_data
-    refinement_data.wave_tolerance = 1.e-5
-    refinement_data.speed_tolerance = [1e-8]
-    refinement_data.variable_dt_refinement_ratios = True
-
-    # for topography, append lines of the form [topotype, fname]
-    topo_data = rundata.topo_data
-    topo_data.topofiles.append([3, '../common-files/utah-flat.asc'])
-
-    # for moving topography, append lines of the form: [topotype, fname]
-    # dtopo_data = rundata.dtopo_data
-
-    # for qinit perturbations, append lines of the form: [fname]
-    rundata.qinit_data.qinit_type = 0
-    rundata.qinit_data.qinitfiles = []
-
-    # for fixed grids append lines of the form
-    # [t1, t2, noutput, x1, x2, y1, y2, xpoints, ypoints, ioutarrivaltimes, ioutsurfacemax]
-    # fixedgrids = rundata.fixed_grid_data
+    # add a topography file; format [topotype, fname]
+    rundata.topo_data.topofiles.append([3, '../common-files/utah-flat.asc'])
 
     # Land-spill module settings
-    rundata.add_data(LandSpillData(), 'landspill_data')
-    landspill = rundata.landspill_data
-    landspill.ref_mu = 332.  # cP @ 15 degree C
-    landspill.ref_temperature = 15.
-    landspill.ambient_temperature = 25.
-    landspill.density = 9.266e2  # kg / m^3 @ 15 degree C; will overwrite rho in GeoClaw
-
-    # extra parameters
-    landspill.update_tol = geo_data.dry_tolerance
-    landspill.refine_tol = 0.0
+    rundata.landspill_data.ref_mu = 332.  # cP @ 15 degree C
+    rundata.landspill_data.ref_temperature = 15.
+    rundata.landspill_data.ambient_temperature = 25.
+    rundata.landspill_data.density = 9.266e2  # kg / m^3 @ 15 degree C; will overwrite rho in GeoClaw
 
     # Point sources
-    ptsources_data = landspill.point_sources
-    ptsources_data.n_point_sources = 1
-    ptsources_data.point_sources.append(
-        [[-12459650., 4986000.], 2, [1800., 12600.], [0.5, 0.1]])
+    rundata.landspill_data.point_sources.n_point_sources = 1
+    rundata.landspill_data.point_sources.point_sources.append([break_point, 2, [1800., 12600.], [0.5, 0.1]])
 
     # Darcy-Weisbach friction
-    darcy_weisbach_data = landspill.darcy_weisbach_friction
-    darcy_weisbach_data.type = 4
-    darcy_weisbach_data.dry_tol = 1e-4
-    darcy_weisbach_data.friction_tol = 1e6
-    darcy_weisbach_data.default_roughness = 0.0
-    darcy_weisbach_data.filename = "roughness.txt"
+    rundata.landspill_data.darcy_weisbach_friction.type = 4
+    rundata.landspill_data.darcy_weisbach_friction.default_roughness = 0.0
+    rundata.landspill_data.darcy_weisbach_friction.filename = "roughness.txt"
 
-    # hydrological features
-    hydro_feature_data = landspill.hydro_features
-    hydro_feature_data.files.append("../common-files/utah-flat-hydro.asc")
+    # add a hydrological feature raster file
+    rundata.landspill_data.hydro_features.files.append("../common-files/utah-flat-hydro.asc")
 
     # Evaporation
-    evaporation_data = landspill.evaporation
-    evaporation_data.type = 0
-    evaporation_data.coefficients = [1.38, 0.045]
+    rundata.landspill_data.evaporation.type = 0
+    rundata.landspill_data.evaporation.coefficients = [1.38, 0.045]
 
     return rundata
diff --git a/cases/utah-flat-maya/animation.gif b/cases/utah-flat-maya/animation.gif
new file mode 100644
index 0000000..a134005
Binary files /dev/null and b/cases/utah-flat-maya/animation.gif differ
diff --git a/cases/utah-flat-maya/setrun.py b/cases/utah-flat-maya/setrun.py
index d61c4d2..a6391f2 100644
--- a/cases/utah-flat-maya/setrun.py
+++ b/cases/utah-flat-maya/setrun.py
@@ -2,304 +2,63 @@
 
 The values set in the function setrun are then written out to data files
 that will be read in by the Fortran code.
-
 """
 # pylint: disable=no-member, too-many-statements
-import functools
-from gclandspill.clawutil import data
-from gclandspill.data import LandSpillData
+import gclandspill.data
 
 
-def setrun(claw_pkg: str = 'geoclaw'):
+def setrun():
     """Define the parameters used for running Clawpack.
 
-    Arguments
-    ---------
-    claw_pkg : str
-        Expected to be "geoclaw" for this setrun.
-
     Returns
     -------
-    rundata : clawutil.ClawRunData
+    rundata : gclandspill.data.ClawRunData
     """
 
-    assert claw_pkg.lower() == 'geoclaw',  "Expected claw_pkg = 'geoclaw'"
-
-    num_dim = 2
-    rundata = data.ClawRunData(claw_pkg, num_dim)
-
-    # GeoClaw specific parameters:
-    rundata = setgeo(rundata)
+    rundata = gclandspill.data.ClawRunData()
 
-    # AmrClaw specific parameters:
-    rundata = setamr(rundata)
+    # for convenience
+    break_point = [-12459650., 4986000.]
 
-    # an alias to object holding standard Clawpack parameters
-    clawdata = rundata.clawdata
-
-    # number of space dimensions:
-    clawdata.num_dim = num_dim
-
-    # lower and upper edge of computational domain:
-    clawdata.lower[0] = -12459650.-400.
-    clawdata.upper[0] = -12459650.+400.
-    clawdata.lower[1] = 4986000.-400
-    clawdata.upper[1] = 4986000.+400
+    # lower and upper edges of computational domain
+    rundata.clawdata.lower[0] = break_point[0]-400.
+    rundata.clawdata.upper[0] = break_point[0]+400.
+    rundata.clawdata.lower[1] = break_point[1]-400
+    rundata.clawdata.upper[1] = break_point[1]+400
 
     # number of grid cells: coarsest grid
-    clawdata.num_cells[0] = 200
-    clawdata.num_cells[1] = 200
-
-    # number of equations in the system: continuity eq. and x, y momemtums
-    clawdata.num_eqn = 3
-
-    # number of auxiliary variables: elevation and hydrological data
-    clawdata.num_aux = 2
-
-    # index of aux array corresponding to capacity function, if there is one:
-    clawdata.capa_index = 0
-
-    # initial time:
-    clawdata.t0 = 0.0
-
-    # restart from checkpoint file of a previous run?
-    clawdata.restart = False               # True to restart from prior results
-    clawdata.restart_file = 'fort.chk00006'  # File to use for restart data
-
-    # temporal data output style
-    clawdata.output_style = 1
-
-    if clawdata.output_style == 1:  # equally spaced times up to tfinal
-        clawdata.num_output_times = 240
-        clawdata.tfinal = 28800
-        clawdata.output_t0 = True  # output at initial (or restart) time?
-    elif clawdata.output_style == 2:  # a list of output times
-        clawdata.output_times = [0., 1800.]
-    elif clawdata.output_style == 3:  # every iout timesteps with a total of ntot time steps:
-        clawdata.output_step_interval = 80
-        clawdata.total_steps = 4000
-        clawdata.output_t0 = True
-
-    # output data format
-    clawdata.output_format = 'binary'
-
-    # which veriable to output
-    clawdata.output_q_components = 'all'  # could be list such as [True,True]
-
-    # which aux data to output
-    clawdata.output_aux_components = 'all'  # could be list
-
-    # whether to output aux data only at t0
-    clawdata.output_aux_onlyonce = False
-
-    # print information up to this AMR level
-    clawdata.verbosity = 3
-
-    # use veriable time step size based on CFL
-    clawdata.dt_variable = 1
-
-    # data required to calculate a proper initial time step size
-    cell_area = \
-        (clawdata.upper[0] - clawdata.lower[0]) * (clawdata.upper[1] - clawdata.lower[1]) / \
-        (clawdata.num_cells[0] * clawdata.num_cells[1])
-    cell_area /= functools.reduce(
-        lambda x, y: x*y, rundata.amrdata.refinement_ratios_x[:rundata.amrdata.amr_levels_max], 1)
-    cell_area /= functools.reduce(
-        lambda x, y: x*y, rundata.amrdata.refinement_ratios_y[:rundata.amrdata.amr_levels_max], 1)
-    vrate = rundata.landspill_data.point_sources.point_sources[0][-1][0]
-
-    # initial time step for variable dt.
-    clawdata.dt_initial = 0.3 * cell_area / vrate
-
-    # max time step allowed if variable dt
-    clawdata.dt_max = 4.0
-
-    # desired CFL if variable dt, and max CFL allowed without retaking a step with a smaller dt
-    clawdata.cfl_desired = 0.9
-    clawdata.cfl_max = 0.95
-
-    # maximum number of time steps allowed between two output times:
-    clawdata.steps_max = 100000
-
-    # order of accuracy:  1 => Godunov,  2 => Lax-Wendroff plus limiters
-    clawdata.order = 2
-
-    # use dimensional splitting? (not yet available for AMR)
-    clawdata.dimensional_split = 'unsplit'
-
-    # 2 or 'all' => corner transport of 2nd order corrections too
-    clawdata.transverse_waves = 2
-
-    # number of waves in the Riemann solution:
-    clawdata.num_waves = 3
-
-    # list of limiters to use for each wave family:
-    clawdata.limiter = ['mc', 'mc', 'mc']
-
-    # use f-wave version of algorithms
-    clawdata.use_fwaves = True
-
-    # src_split == 1 or 'godunov' ==> Godunov (1st order) splitting used,
-    clawdata.source_split = 'godunov'
-
-    # number of ghost cells (usually 2)
-    clawdata.num_ghost = 2
+    rundata.clawdata.num_cells[0] = 200
+    rundata.clawdata.num_cells[1] = 200
 
-    # BC: 1 => extrapolation (non-reflecting outflow)
-    clawdata.bc_lower[0] = 1
-    clawdata.bc_upper[0] = 1
-    clawdata.bc_lower[1] = 1
-    clawdata.bc_upper[1] = 1
+    # temporal data output
+    rundata.clawdata.output_style = 1
+    rundata.clawdata.num_output_times = 240
+    rundata.clawdata.tfinal = 28800
+    rundata.clawdata.output_t0 = True  # output at initial (or restart) time?
 
-    # checkout mechanism selection
-    clawdata.checkpt_style = 0
-
-    if clawdata.checkpt_style == 0:  # do not checkpoint at all
-        pass
-    elif clawdata.checkpt_style == 1:  # checkpoint only at tfinal.
-        pass
-    elif clawdata.checkpt_style == 2:  # specify a list of checkpoint times.
-        clawdata.checkpt_times = [0.1, 0.15]
-    elif clawdata.checkpt_style == 3:  # checkpoint every checkpt_interval timesteps (on Level 1)
-        clawdata.checkpt_interval = 5
-
-    return rundata
-
-
-def setamr(rundata: data.ClawRunData):
-    """Set AmrClaw specific runtime parameters."""
-
-    try:
-        amrdata = rundata.amrdata
-    except AttributeError as err:
-        raise AttributeError("This rundata has no amrdata attribute") from err
-
-    # max number of refinement levels:
-    amrdata.amr_levels_max = 2
-
-    # List of refinement ratios at each level (length at least mxnest-1)
-    amrdata.refinement_ratios_x = [4]
-    amrdata.refinement_ratios_y = [4]
-    amrdata.refinement_ratios_t = [4]
-
-    # aux variables are defined at cell centers
-    amrdata.aux_type = ['center', 'center']
-
-    # use flag2refine routine to flag refinement (rather than using richardson error)
-    amrdata.flag_richardson = False
-    amrdata.flag2refine = True
-
-    # steps to take on each level L between regriddings of level L+1:
-    amrdata.regrid_interval = 1
-
-    # width of buffer zone around flagged points
-    amrdata.regrid_buffer_width = 1
-
-    # clustering alg. cutoff for (# flagged pts) / (total # of cells refined)
-    amrdata.clustering_cutoff = 0.80000
-
-    # print info about each regridding up to this level:
-    amrdata.verbosity_regrid = 0
-
-    #  ----- For developers -----
-    # Toggle debugging print statements:
-    amrdata.dprint = False      # print domain flags
-    amrdata.eprint = False      # print err est flags
-    amrdata.edebug = False      # even more err est flags
-    amrdata.gprint = False      # grid bisection/clustering
-    amrdata.nprint = False      # proper nesting output
-    amrdata.pprint = False      # proj. of tagged points
-    amrdata.rprint = False      # print regridding summary
-    amrdata.sprint = False      # space/memory output
-    amrdata.tprint = False      # time step reporting each level
-    amrdata.uprint = False      # update/upbnd reporting
-
-    # customize regions of refinement
-    # rundata.regiondata.regions.append([minlevel, maxlevel, t1, t2, x1, x2, y1, y2])
-
-    # gauges
-    # rundata.gaugedata.gauges.append([gaugeno, x, y, t1, t2])
-
-    return rundata
-
-
-def setgeo(rundata: data.ClawRunData):
-    """Set GeoClaw specific runtime parameters."""
-
-    try:
-        geo_data = rundata.geo_data
-    except AttributeError as err:
-        raise AttributeError("This rundata has no geo_data attribute") from err
-
-    # == Physics ==
-    geo_data.gravity = 9.81
-    geo_data.coordinate_system = 1
-    geo_data.earth_radius = 6367.5e3
-
-    # == Forcing Options
-    geo_data.coriolis_forcing = False
-
-    # == Algorithm and Initial Conditions ==
-    geo_data.sea_level = 0.0
-    geo_data.dry_tolerance = 1.e-4
-    geo_data.friction_forcing = False
-    geo_data.manning_coefficient = 0.035
-    geo_data.friction_depth = 1.e6
-
-    # Refinement data
-    refinement_data = rundata.refinement_data
-    refinement_data.wave_tolerance = 1.e-5
-    refinement_data.speed_tolerance = [1e-8]
-    refinement_data.variable_dt_refinement_ratios = True
-
-    # for topography, append lines of the form [topotype, fname]
-    topo_data = rundata.topo_data
-    topo_data.topofiles.append([3, '../common-files/utah-flat.asc'])
-
-    # for moving topography, append lines of the form: [topotype, fname]
-    # dtopo_data = rundata.dtopo_data
-
-    # for qinit perturbations, append lines of the form: [fname]
-    rundata.qinit_data.qinit_type = 0
-    rundata.qinit_data.qinitfiles = []
-
-    # for fixed grids append lines of the form
-    # [t1, t2, noutput, x1, x2, y1, y2, xpoints, ypoints, ioutarrivaltimes, ioutsurfacemax]
-    # fixedgrids = rundata.fixed_grid_data
+    # add a topography file; format [topotype, fname]
+    rundata.topo_data.topofiles.append([3, '../common-files/utah-flat.asc'])
 
     # Land-spill module settings
-    rundata.add_data(LandSpillData(), 'landspill_data')
-    landspill = rundata.landspill_data
-    landspill.ref_mu = 332.  # cP @ 15 degree C
-    landspill.ref_temperature = 15.
-    landspill.ambient_temperature = 25.
-    landspill.density = 9.266e2  # kg / m^3 @ 15 degree C; will overwrite rho in GeoClaw
-
-    # extra parameters
-    landspill.update_tol = geo_data.dry_tolerance
-    landspill.refine_tol = 0.0
+    rundata.landspill_data.ref_mu = 332.  # cP @ 15 degree C
+    rundata.landspill_data.ref_temperature = 15.
+    rundata.landspill_data.ambient_temperature = 25.
+    rundata.landspill_data.density = 9.266e2  # kg / m^3 @ 15 degree C; will overwrite rho in GeoClaw
 
     # Point sources
-    ptsources_data = landspill.point_sources
-    ptsources_data.n_point_sources = 1
-    ptsources_data.point_sources.append(
-        [[-12459650., 4986000.], 2, [1800., 12600.], [0.5, 0.1]])
+    rundata.landspill_data.point_sources.n_point_sources = 1
+    rundata.landspill_data.point_sources.point_sources.append([break_point, 2, [1800., 12600.], [0.5, 0.1]])
 
     # Darcy-Weisbach friction
-    darcy_weisbach_data = landspill.darcy_weisbach_friction
-    darcy_weisbach_data.type = 4
-    darcy_weisbach_data.dry_tol = 1e-4
-    darcy_weisbach_data.friction_tol = 1e6
-    darcy_weisbach_data.default_roughness = 0.0
-    darcy_weisbach_data.filename = "roughness.txt"
+    rundata.landspill_data.darcy_weisbach_friction.type = 4
+    rundata.landspill_data.darcy_weisbach_friction.default_roughness = 0.0
+    rundata.landspill_data.darcy_weisbach_friction.filename = "roughness.txt"
 
-    # hydrological features
-    hydro_feature_data = landspill.hydro_features
-    hydro_feature_data.files.append("../common-files/utah-flat-hydro.asc")
+    # add a hydrological feature raster file
+    rundata.landspill_data.hydro_features.files.append("../common-files/utah-flat-hydro.asc")
 
     # Evaporation
-    evaporation_data = landspill.evaporation
-    evaporation_data.type = 1
-    evaporation_data.coefficients = [1.38, 0.045]
+    rundata.landspill_data.evaporation.type = 1
+    rundata.landspill_data.evaporation.coefficients = [1.38, 0.045]
 
     return rundata
diff --git a/cases/utah-hill-maya-no-hydro/animation.gif b/cases/utah-hill-maya-no-hydro/animation.gif
new file mode 100644
index 0000000..0ea33e6
Binary files /dev/null and b/cases/utah-hill-maya-no-hydro/animation.gif differ
diff --git a/cases/utah-hill-maya-no-hydro/setrun.py b/cases/utah-hill-maya-no-hydro/setrun.py
index aad8cf3..da72309 100644
--- a/cases/utah-hill-maya-no-hydro/setrun.py
+++ b/cases/utah-hill-maya-no-hydro/setrun.py
@@ -2,300 +2,60 @@
 
 The values set in the function setrun are then written out to data files
 that will be read in by the Fortran code.
-
 """
 # pylint: disable=no-member, too-many-statements
-import functools
-from gclandspill.clawutil import data
-from gclandspill.data import LandSpillData
+import gclandspill.data
 
 
-def setrun(claw_pkg: str = 'geoclaw'):
+def setrun():
     """Define the parameters used for running Clawpack.
 
-    Arguments
-    ---------
-    claw_pkg : str
-        Expected to be "geoclaw" for this setrun.
-
     Returns
     -------
-    rundata : clawutil.ClawRunData
+    rundata : gclandspill.data.ClawRunData
     """
 
-    assert claw_pkg.lower() == 'geoclaw',  "Expected claw_pkg = 'geoclaw'"
-
-    num_dim = 2
-    rundata = data.ClawRunData(claw_pkg, num_dim)
-
-    # GeoClaw specific parameters:
-    rundata = setgeo(rundata)
+    rundata = gclandspill.data.ClawRunData()
 
-    # AmrClaw specific parameters:
-    rundata = setamr(rundata)
+    # for convenience
+    break_point = [-12443619., 4977641.]
 
-    # an alias to object holding standard Clawpack parameters
-    clawdata = rundata.clawdata
-
-    # number of space dimensions:
-    clawdata.num_dim = num_dim
-
-    # lower and upper edge of computational domain:
-    clawdata.lower[0] = -12443619.-100.
-    clawdata.upper[0] = -12443619.+700.
-    clawdata.lower[1] = 4977641.-650
-    clawdata.upper[1] = 4977641.+150
+    # lower and upper edges of computational domain
+    rundata.clawdata.lower[0] = break_point[0]-100.
+    rundata.clawdata.upper[0] = break_point[0]+700.
+    rundata.clawdata.lower[1] = break_point[1]-650
+    rundata.clawdata.upper[1] = break_point[1]+150
 
     # number of grid cells: coarsest grid
-    clawdata.num_cells[0] = 200
-    clawdata.num_cells[1] = 200
-
-    # number of equations in the system: continuity eq. and x, y momemtums
-    clawdata.num_eqn = 3
-
-    # number of auxiliary variables: elevation and hydrological data
-    clawdata.num_aux = 2
-
-    # index of aux array corresponding to capacity function, if there is one:
-    clawdata.capa_index = 0
-
-    # initial time:
-    clawdata.t0 = 0.0
-
-    # restart from checkpoint file of a previous run?
-    clawdata.restart = False               # True to restart from prior results
-    clawdata.restart_file = 'fort.chk00006'  # File to use for restart data
-
-    # temporal data output style
-    clawdata.output_style = 1
-
-    if clawdata.output_style == 1:  # equally spaced times up to tfinal
-        clawdata.num_output_times = 240
-        clawdata.tfinal = 28800
-        clawdata.output_t0 = True  # output at initial (or restart) time?
-    elif clawdata.output_style == 2:  # a list of output times
-        clawdata.output_times = [0., 1800.]
-    elif clawdata.output_style == 3:  # every iout timesteps with a total of ntot time steps:
-        clawdata.output_step_interval = 80
-        clawdata.total_steps = 4000
-        clawdata.output_t0 = True
-
-    # output data format
-    clawdata.output_format = 'binary'
-
-    # which veriable to output
-    clawdata.output_q_components = 'all'  # could be list such as [True,True]
-
-    # which aux data to output
-    clawdata.output_aux_components = 'all'  # could be list
-
-    # whether to output aux data only at t0
-    clawdata.output_aux_onlyonce = False
-
-    # print information up to this AMR level
-    clawdata.verbosity = 3
-
-    # use veriable time step size based on CFL
-    clawdata.dt_variable = 1
-
-    # data required to calculate a proper initial time step size
-    cell_area = \
-        (clawdata.upper[0] - clawdata.lower[0]) * (clawdata.upper[1] - clawdata.lower[1]) / \
-        (clawdata.num_cells[0] * clawdata.num_cells[1])
-    cell_area /= functools.reduce(
-        lambda x, y: x*y, rundata.amrdata.refinement_ratios_x[:rundata.amrdata.amr_levels_max], 1)
-    cell_area /= functools.reduce(
-        lambda x, y: x*y, rundata.amrdata.refinement_ratios_y[:rundata.amrdata.amr_levels_max], 1)
-    vrate = rundata.landspill_data.point_sources.point_sources[0][-1][0]
-
-    # initial time step for variable dt.
-    clawdata.dt_initial = 0.3 * cell_area / vrate
-
-    # max time step allowed if variable dt
-    clawdata.dt_max = 2.5
-
-    # desired CFL if variable dt, and max CFL allowed without retaking a step with a smaller dt
-    clawdata.cfl_desired = 0.9
-    clawdata.cfl_max = 0.95
-
-    # maximum number of time steps allowed between two output times:
-    clawdata.steps_max = 100000
-
-    # order of accuracy:  1 => Godunov,  2 => Lax-Wendroff plus limiters
-    clawdata.order = 2
-
-    # use dimensional splitting? (not yet available for AMR)
-    clawdata.dimensional_split = 'unsplit'
-
-    # 2 or 'all' => corner transport of 2nd order corrections too
-    clawdata.transverse_waves = 2
-
-    # number of waves in the Riemann solution:
-    clawdata.num_waves = 3
-
-    # list of limiters to use for each wave family:
-    clawdata.limiter = ['mc', 'mc', 'mc']
-
-    # use f-wave version of algorithms
-    clawdata.use_fwaves = True
-
-    # src_split == 1 or 'godunov' ==> Godunov (1st order) splitting used,
-    clawdata.source_split = 'godunov'
-
-    # number of ghost cells (usually 2)
-    clawdata.num_ghost = 2
+    rundata.clawdata.num_cells[0] = 200
+    rundata.clawdata.num_cells[1] = 200
 
-    # BC: 1 => extrapolation (non-reflecting outflow)
-    clawdata.bc_lower[0] = 1
-    clawdata.bc_upper[0] = 1
-    clawdata.bc_lower[1] = 1
-    clawdata.bc_upper[1] = 1
+    # temporal data output
+    rundata.clawdata.output_style = 1
+    rundata.clawdata.num_output_times = 240
+    rundata.clawdata.tfinal = 28800
+    rundata.clawdata.output_t0 = True  # output at initial (or restart) time?
 
-    # checkout mechanism selection
-    clawdata.checkpt_style = 0
-
-    if clawdata.checkpt_style == 0:  # do not checkpoint at all
-        pass
-    elif clawdata.checkpt_style == 1:  # checkpoint only at tfinal.
-        pass
-    elif clawdata.checkpt_style == 2:  # specify a list of checkpoint times.
-        clawdata.checkpt_times = [0.1, 0.15]
-    elif clawdata.checkpt_style == 3:  # checkpoint every checkpt_interval timesteps (on Level 1)
-        clawdata.checkpt_interval = 5
-
-    return rundata
-
-
-def setamr(rundata: data.ClawRunData):
-    """Set AmrClaw specific runtime parameters."""
-
-    try:
-        amrdata = rundata.amrdata
-    except AttributeError as err:
-        raise AttributeError("This rundata has no amrdata attribute") from err
-
-    # max number of refinement levels:
-    amrdata.amr_levels_max = 2
-
-    # List of refinement ratios at each level (length at least mxnest-1)
-    amrdata.refinement_ratios_x = [4]
-    amrdata.refinement_ratios_y = [4]
-    amrdata.refinement_ratios_t = [4]
-
-    # aux variables are defined at cell centers
-    amrdata.aux_type = ['center', 'center']
-
-    # use flag2refine routine to flag refinement (rather than using richardson error)
-    amrdata.flag_richardson = False
-    amrdata.flag2refine = True
-
-    # steps to take on each level L between regriddings of level L+1:
-    amrdata.regrid_interval = 1
-
-    # width of buffer zone around flagged points
-    amrdata.regrid_buffer_width = 1
-
-    # clustering alg. cutoff for (# flagged pts) / (total # of cells refined)
-    amrdata.clustering_cutoff = 0.80000
-
-    # print info about each regridding up to this level:
-    amrdata.verbosity_regrid = 0
-
-    #  ----- For developers -----
-    # Toggle debugging print statements:
-    amrdata.dprint = False      # print domain flags
-    amrdata.eprint = False      # print err est flags
-    amrdata.edebug = False      # even more err est flags
-    amrdata.gprint = False      # grid bisection/clustering
-    amrdata.nprint = False      # proper nesting output
-    amrdata.pprint = False      # proj. of tagged points
-    amrdata.rprint = False      # print regridding summary
-    amrdata.sprint = False      # space/memory output
-    amrdata.tprint = False      # time step reporting each level
-    amrdata.uprint = False      # update/upbnd reporting
-
-    # customize regions of refinement
-    # rundata.regiondata.regions.append([minlevel, maxlevel, t1, t2, x1, x2, y1, y2])
-
-    # gauges
-    # rundata.gaugedata.gauges.append([gaugeno, x, y, t1, t2])
-
-    return rundata
-
-
-def setgeo(rundata: data.ClawRunData):
-    """Set GeoClaw specific runtime parameters."""
-
-    try:
-        geo_data = rundata.geo_data
-    except AttributeError as err:
-        raise AttributeError("This rundata has no geo_data attribute") from err
-
-    # == Physics ==
-    geo_data.gravity = 9.81
-    geo_data.coordinate_system = 1
-    geo_data.earth_radius = 6367.5e3
-
-    # == Forcing Options
-    geo_data.coriolis_forcing = False
-
-    # == Algorithm and Initial Conditions ==
-    geo_data.sea_level = 0.0
-    geo_data.dry_tolerance = 1.e-4
-    geo_data.friction_forcing = False
-    geo_data.manning_coefficient = 0.035
-    geo_data.friction_depth = 1.e6
-
-    # Refinement data
-    refinement_data = rundata.refinement_data
-    refinement_data.wave_tolerance = 1.e-5
-    refinement_data.speed_tolerance = [1e-8]
-    refinement_data.variable_dt_refinement_ratios = True
-
-    # for topography, append lines of the form [topotype, fname]
-    topo_data = rundata.topo_data
-    topo_data.topofiles.append([3, '../common-files/utah-hill.asc'])
-
-    # for moving topography, append lines of the form: [topotype, fname]
-    # dtopo_data = rundata.dtopo_data
-
-    # for qinit perturbations, append lines of the form: [fname]
-    rundata.qinit_data.qinit_type = 0
-    rundata.qinit_data.qinitfiles = []
-
-    # for fixed grids append lines of the form
-    # [t1, t2, noutput, x1, x2, y1, y2, xpoints, ypoints, ioutarrivaltimes, ioutsurfacemax]
-    # fixedgrids = rundata.fixed_grid_data
+    # add a topography file; format [topotype, fname]
+    rundata.topo_data.topofiles.append([3, '../common-files/utah-hill.asc'])
 
     # Land-spill module settings
-    rundata.add_data(LandSpillData(), 'landspill_data')
-    landspill = rundata.landspill_data
-    landspill.ref_mu = 332.  # cP @ 15 degree C
-    landspill.ref_temperature = 15.
-    landspill.ambient_temperature = 25.
-    landspill.density = 9.266e2  # kg / m^3 @ 15 degree C; will overwrite rho in GeoClaw
-
-    # extra parameters
-    landspill.update_tol = geo_data.dry_tolerance
-    landspill.refine_tol = 0.0
+    rundata.landspill_data.ref_mu = 332.  # cP @ 15 degree C
+    rundata.landspill_data.ref_temperature = 15.
+    rundata.landspill_data.ambient_temperature = 25.
+    rundata.landspill_data.density = 9.266e2  # kg / m^3 @ 15 degree C; will overwrite rho in GeoClaw
 
     # Point sources
-    ptsources_data = landspill.point_sources
-    ptsources_data.n_point_sources = 1
-    ptsources_data.point_sources.append(
-        [[-12443619., 4977641.], 2, [1800., 12600.], [0.5, 0.1]])
+    rundata.landspill_data.point_sources.n_point_sources = 1
+    rundata.landspill_data.point_sources.point_sources.append([break_point, 2, [1800., 12600.], [0.5, 0.1]])
 
     # Darcy-Weisbach friction
-    darcy_weisbach_data = landspill.darcy_weisbach_friction
-    darcy_weisbach_data.type = 4
-    darcy_weisbach_data.dry_tol = 1e-4
-    darcy_weisbach_data.friction_tol = 1e6
-    darcy_weisbach_data.default_roughness = 0.0
-    darcy_weisbach_data.filename = "roughness.txt"
+    rundata.landspill_data.darcy_weisbach_friction.type = 4
+    rundata.landspill_data.darcy_weisbach_friction.default_roughness = 0.0
+    rundata.landspill_data.darcy_weisbach_friction.filename = "roughness.txt"
 
     # Evaporation
-    evaporation_data = landspill.evaporation
-    evaporation_data.type = 1
-    evaporation_data.coefficients = [1.38, 0.045]
+    rundata.landspill_data.evaporation.type = 1
+    rundata.landspill_data.evaporation.coefficients = [1.38, 0.045]
 
     return rundata
diff --git a/cases/utah-hill-maya/animation.gif b/cases/utah-hill-maya/animation.gif
new file mode 100644
index 0000000..8cf4e6b
Binary files /dev/null and b/cases/utah-hill-maya/animation.gif differ
diff --git a/cases/utah-hill-maya/setrun.py b/cases/utah-hill-maya/setrun.py
index e589a29..3cd958e 100644
--- a/cases/utah-hill-maya/setrun.py
+++ b/cases/utah-hill-maya/setrun.py
@@ -2,304 +2,63 @@
 
 The values set in the function setrun are then written out to data files
 that will be read in by the Fortran code.
-
 """
 # pylint: disable=no-member, too-many-statements
-import functools
-from gclandspill.clawutil import data
-from gclandspill.data import LandSpillData
+import gclandspill.data
 
 
-def setrun(claw_pkg: str = 'geoclaw'):
+def setrun():
     """Define the parameters used for running Clawpack.
 
-    Arguments
-    ---------
-    claw_pkg : str
-        Expected to be "geoclaw" for this setrun.
-
     Returns
     -------
-    rundata : clawutil.ClawRunData
+    rundata : gclandspill.data.ClawRunData
     """
 
-    assert claw_pkg.lower() == 'geoclaw',  "Expected claw_pkg = 'geoclaw'"
-
-    num_dim = 2
-    rundata = data.ClawRunData(claw_pkg, num_dim)
-
-    # GeoClaw specific parameters:
-    rundata = setgeo(rundata)
+    rundata = gclandspill.data.ClawRunData()
 
-    # AmrClaw specific parameters:
-    rundata = setamr(rundata)
+    # for convenience
+    break_point = [-12443619., 4977641.]
 
-    # an alias to object holding standard Clawpack parameters
-    clawdata = rundata.clawdata
-
-    # number of space dimensions:
-    clawdata.num_dim = num_dim
-
-    # lower and upper edge of computational domain:
-    clawdata.lower[0] = -12443619.-100.
-    clawdata.upper[0] = -12443619.+700.
-    clawdata.lower[1] = 4977641.-650
-    clawdata.upper[1] = 4977641.+150
+    # lower and upper edges of computational domain
+    rundata.clawdata.lower[0] = break_point[0]-100.
+    rundata.clawdata.upper[0] = break_point[0]+700.
+    rundata.clawdata.lower[1] = break_point[1]-650
+    rundata.clawdata.upper[1] = break_point[1]+150
 
     # number of grid cells: coarsest grid
-    clawdata.num_cells[0] = 200
-    clawdata.num_cells[1] = 200
-
-    # number of equations in the system: continuity eq. and x, y momemtums
-    clawdata.num_eqn = 3
-
-    # number of auxiliary variables: elevation and hydrological data
-    clawdata.num_aux = 2
-
-    # index of aux array corresponding to capacity function, if there is one:
-    clawdata.capa_index = 0
-
-    # initial time:
-    clawdata.t0 = 0.0
-
-    # restart from checkpoint file of a previous run?
-    clawdata.restart = False               # True to restart from prior results
-    clawdata.restart_file = 'fort.chk00006'  # File to use for restart data
-
-    # temporal data output style
-    clawdata.output_style = 1
-
-    if clawdata.output_style == 1:  # equally spaced times up to tfinal
-        clawdata.num_output_times = 240
-        clawdata.tfinal = 28800
-        clawdata.output_t0 = True  # output at initial (or restart) time?
-    elif clawdata.output_style == 2:  # a list of output times
-        clawdata.output_times = [0., 1800.]
-    elif clawdata.output_style == 3:  # every iout timesteps with a total of ntot time steps:
-        clawdata.output_step_interval = 80
-        clawdata.total_steps = 4000
-        clawdata.output_t0 = True
-
-    # output data format
-    clawdata.output_format = 'binary'
-
-    # which veriable to output
-    clawdata.output_q_components = 'all'  # could be list such as [True,True]
-
-    # which aux data to output
-    clawdata.output_aux_components = 'all'  # could be list
-
-    # whether to output aux data only at t0
-    clawdata.output_aux_onlyonce = False
-
-    # print information up to this AMR level
-    clawdata.verbosity = 3
-
-    # use veriable time step size based on CFL
-    clawdata.dt_variable = 1
-
-    # data required to calculate a proper initial time step size
-    cell_area = \
-        (clawdata.upper[0] - clawdata.lower[0]) * (clawdata.upper[1] - clawdata.lower[1]) / \
-        (clawdata.num_cells[0] * clawdata.num_cells[1])
-    cell_area /= functools.reduce(
-        lambda x, y: x*y, rundata.amrdata.refinement_ratios_x[:rundata.amrdata.amr_levels_max], 1)
-    cell_area /= functools.reduce(
-        lambda x, y: x*y, rundata.amrdata.refinement_ratios_y[:rundata.amrdata.amr_levels_max], 1)
-    vrate = rundata.landspill_data.point_sources.point_sources[0][-1][0]
-
-    # initial time step for variable dt.
-    clawdata.dt_initial = 0.3 * cell_area / vrate
-
-    # max time step allowed if variable dt
-    clawdata.dt_max = 2.5
-
-    # desired CFL if variable dt, and max CFL allowed without retaking a step with a smaller dt
-    clawdata.cfl_desired = 0.9
-    clawdata.cfl_max = 0.95
-
-    # maximum number of time steps allowed between two output times:
-    clawdata.steps_max = 100000
-
-    # order of accuracy:  1 => Godunov,  2 => Lax-Wendroff plus limiters
-    clawdata.order = 2
-
-    # use dimensional splitting? (not yet available for AMR)
-    clawdata.dimensional_split = 'unsplit'
-
-    # 2 or 'all' => corner transport of 2nd order corrections too
-    clawdata.transverse_waves = 2
-
-    # number of waves in the Riemann solution:
-    clawdata.num_waves = 3
-
-    # list of limiters to use for each wave family:
-    clawdata.limiter = ['mc', 'mc', 'mc']
-
-    # use f-wave version of algorithms
-    clawdata.use_fwaves = True
-
-    # src_split == 1 or 'godunov' ==> Godunov (1st order) splitting used,
-    clawdata.source_split = 'godunov'
-
-    # number of ghost cells (usually 2)
-    clawdata.num_ghost = 2
+    rundata.clawdata.num_cells[0] = 200
+    rundata.clawdata.num_cells[1] = 200
 
-    # BC: 1 => extrapolation (non-reflecting outflow)
-    clawdata.bc_lower[0] = 1
-    clawdata.bc_upper[0] = 1
-    clawdata.bc_lower[1] = 1
-    clawdata.bc_upper[1] = 1
+    # temporal data output
+    rundata.clawdata.output_style = 1
+    rundata.clawdata.num_output_times = 240
+    rundata.clawdata.tfinal = 28800
+    rundata.clawdata.output_t0 = True  # output at initial (or restart) time?
 
-    # checkout mechanism selection
-    clawdata.checkpt_style = 0
-
-    if clawdata.checkpt_style == 0:  # do not checkpoint at all
-        pass
-    elif clawdata.checkpt_style == 1:  # checkpoint only at tfinal.
-        pass
-    elif clawdata.checkpt_style == 2:  # specify a list of checkpoint times.
-        clawdata.checkpt_times = [0.1, 0.15]
-    elif clawdata.checkpt_style == 3:  # checkpoint every checkpt_interval timesteps (on Level 1)
-        clawdata.checkpt_interval = 5
-
-    return rundata
-
-
-def setamr(rundata: data.ClawRunData):
-    """Set AmrClaw specific runtime parameters."""
-
-    try:
-        amrdata = rundata.amrdata
-    except AttributeError as err:
-        raise AttributeError("This rundata has no amrdata attribute") from err
-
-    # max number of refinement levels:
-    amrdata.amr_levels_max = 2
-
-    # List of refinement ratios at each level (length at least mxnest-1)
-    amrdata.refinement_ratios_x = [4]
-    amrdata.refinement_ratios_y = [4]
-    amrdata.refinement_ratios_t = [4]
-
-    # aux variables are defined at cell centers
-    amrdata.aux_type = ['center', 'center']
-
-    # use flag2refine routine to flag refinement (rather than using richardson error)
-    amrdata.flag_richardson = False
-    amrdata.flag2refine = True
-
-    # steps to take on each level L between regriddings of level L+1:
-    amrdata.regrid_interval = 1
-
-    # width of buffer zone around flagged points
-    amrdata.regrid_buffer_width = 1
-
-    # clustering alg. cutoff for (# flagged pts) / (total # of cells refined)
-    amrdata.clustering_cutoff = 0.80000
-
-    # print info about each regridding up to this level:
-    amrdata.verbosity_regrid = 0
-
-    #  ----- For developers -----
-    # Toggle debugging print statements:
-    amrdata.dprint = False      # print domain flags
-    amrdata.eprint = False      # print err est flags
-    amrdata.edebug = False      # even more err est flags
-    amrdata.gprint = False      # grid bisection/clustering
-    amrdata.nprint = False      # proper nesting output
-    amrdata.pprint = False      # proj. of tagged points
-    amrdata.rprint = False      # print regridding summary
-    amrdata.sprint = False      # space/memory output
-    amrdata.tprint = False      # time step reporting each level
-    amrdata.uprint = False      # update/upbnd reporting
-
-    # customize regions of refinement
-    # rundata.regiondata.regions.append([minlevel, maxlevel, t1, t2, x1, x2, y1, y2])
-
-    # gauges
-    # rundata.gaugedata.gauges.append([gaugeno, x, y, t1, t2])
-
-    return rundata
-
-
-def setgeo(rundata: data.ClawRunData):
-    """Set GeoClaw specific runtime parameters."""
-
-    try:
-        geo_data = rundata.geo_data
-    except AttributeError as err:
-        raise AttributeError("This rundata has no geo_data attribute") from err
-
-    # == Physics ==
-    geo_data.gravity = 9.81
-    geo_data.coordinate_system = 1
-    geo_data.earth_radius = 6367.5e3
-
-    # == Forcing Options
-    geo_data.coriolis_forcing = False
-
-    # == Algorithm and Initial Conditions ==
-    geo_data.sea_level = 0.0
-    geo_data.dry_tolerance = 1.e-4
-    geo_data.friction_forcing = False
-    geo_data.manning_coefficient = 0.035
-    geo_data.friction_depth = 1.e6
-
-    # Refinement data
-    refinement_data = rundata.refinement_data
-    refinement_data.wave_tolerance = 1.e-5
-    refinement_data.speed_tolerance = [1e-8]
-    refinement_data.variable_dt_refinement_ratios = True
-
-    # for topography, append lines of the form [topotype, fname]
-    topo_data = rundata.topo_data
-    topo_data.topofiles.append([3, '../common-files/utah-hill.asc'])
-
-    # for moving topography, append lines of the form: [topotype, fname]
-    # dtopo_data = rundata.dtopo_data
-
-    # for qinit perturbations, append lines of the form: [fname]
-    rundata.qinit_data.qinit_type = 0
-    rundata.qinit_data.qinitfiles = []
-
-    # for fixed grids append lines of the form
-    # [t1, t2, noutput, x1, x2, y1, y2, xpoints, ypoints, ioutarrivaltimes, ioutsurfacemax]
-    # fixedgrids = rundata.fixed_grid_data
+    # add a topography file; format [topotype, fname]
+    rundata.topo_data.topofiles.append([3, '../common-files/utah-hill.asc'])
 
     # Land-spill module settings
-    rundata.add_data(LandSpillData(), 'landspill_data')
-    landspill = rundata.landspill_data
-    landspill.ref_mu = 332.  # cP @ 15 degree C
-    landspill.ref_temperature = 15.
-    landspill.ambient_temperature = 25.
-    landspill.density = 9.266e2  # kg / m^3 @ 15 degree C; will overwrite rho in GeoClaw
-
-    # extra parameters
-    landspill.update_tol = geo_data.dry_tolerance
-    landspill.refine_tol = 0.0
+    rundata.landspill_data.ref_mu = 332.  # cP @ 15 degree C
+    rundata.landspill_data.ref_temperature = 15.
+    rundata.landspill_data.ambient_temperature = 25.
+    rundata.landspill_data.density = 9.266e2  # kg / m^3 @ 15 degree C; will overwrite rho in GeoClaw
 
     # Point sources
-    ptsources_data = landspill.point_sources
-    ptsources_data.n_point_sources = 1
-    ptsources_data.point_sources.append(
-        [[-12443619., 4977641.], 2, [1800., 12600.], [0.5, 0.1]])
+    rundata.landspill_data.point_sources.n_point_sources = 1
+    rundata.landspill_data.point_sources.point_sources.append([break_point, 2, [1800., 12600.], [0.5, 0.1]])
 
     # Darcy-Weisbach friction
-    darcy_weisbach_data = landspill.darcy_weisbach_friction
-    darcy_weisbach_data.type = 4
-    darcy_weisbach_data.dry_tol = 1e-4
-    darcy_weisbach_data.friction_tol = 1e6
-    darcy_weisbach_data.default_roughness = 0.0
-    darcy_weisbach_data.filename = "roughness.txt"
+    rundata.landspill_data.darcy_weisbach_friction.type = 4
+    rundata.landspill_data.darcy_weisbach_friction.default_roughness = 0.0
+    rundata.landspill_data.darcy_weisbach_friction.filename = "roughness.txt"
 
-    # hydrological features
-    hydro_feature_data = landspill.hydro_features
-    hydro_feature_data.files.append("../common-files/utah-hill-hydro.asc")
+    # add a hydrological feature raster file
+    rundata.landspill_data.hydro_features.files.append("../common-files/utah-hill-hydro.asc")
 
     # Evaporation
-    evaporation_data = landspill.evaporation
-    evaporation_data.type = 1
-    evaporation_data.coefficients = [1.38, 0.045]
+    rundata.landspill_data.evaporation.type = 1
+    rundata.landspill_data.evaporation.coefficients = [1.38, 0.045]
 
     return rundata
diff --git a/cases/utah-structures-maya/animation.gif b/cases/utah-structures-maya/animation.gif
new file mode 100644
index 0000000..2ceb957
Binary files /dev/null and b/cases/utah-structures-maya/animation.gif differ
diff --git a/cases/utah-structures/roughness.txt b/cases/utah-structures-maya/roughness.txt
similarity index 100%
rename from cases/utah-structures/roughness.txt
rename to cases/utah-structures-maya/roughness.txt
diff --git a/cases/utah-structures-maya/setrun.py b/cases/utah-structures-maya/setrun.py
new file mode 100644
index 0000000..76d772a
--- /dev/null
+++ b/cases/utah-structures-maya/setrun.py
@@ -0,0 +1,64 @@
+"""Module to set up run time parameters for Clawpack.
+
+The values set in the function setrun are then written out to data files
+that will be read in by the Fortran code.
+"""
+# pylint: disable=no-member, too-many-statements
+import gclandspill.data
+
+
+def setrun():
+    """Define the parameters used for running Clawpack.
+
+    Returns
+    -------
+    rundata : gclandspill.data.ClawRunData
+    """
+
+    rundata = gclandspill.data.ClawRunData()
+
+    # for convenience
+    break_point = [-12458815.605194421, 4987307.342152054]
+
+    # lower and upper edges of computational domain
+    rundata.clawdata.lower[0] = break_point[0]-1000.
+    rundata.clawdata.upper[0] = break_point[0]+1000.
+    rundata.clawdata.lower[1] = break_point[1]-1000
+    rundata.clawdata.upper[1] = break_point[1]+1000
+
+    # number of grid cells: coarsest grid
+    rundata.clawdata.num_cells[0] = 500
+    rundata.clawdata.num_cells[1] = 500
+
+    # temporal data output
+    rundata.clawdata.output_style = 1
+    rundata.clawdata.num_output_times = 240
+    rundata.clawdata.tfinal = 28800
+    rundata.clawdata.output_t0 = True  # output at initial (or restart) time?
+
+    # add a topography file; format [topotype, fname]
+    rundata.topo_data.topofiles.append([3, '../common-files/utah-structures.asc'])
+
+    # Land-spill module settings
+    rundata.landspill_data.ref_mu = 332.  # cP @ 15 degree C
+    rundata.landspill_data.ref_temperature = 15.
+    rundata.landspill_data.ambient_temperature = 25.
+    rundata.landspill_data.density = 9.266e2  # kg / m^3 @ 15 degree C; will overwrite rho in GeoClaw
+
+    # Point sources
+    rundata.landspill_data.point_sources.n_point_sources = 1
+    rundata.landspill_data.point_sources.point_sources.append([break_point, 2, [1800., 12600.], [0.5, 0.1]])
+
+    # Darcy-Weisbach friction
+    rundata.landspill_data.darcy_weisbach_friction.type = 4
+    rundata.landspill_data.darcy_weisbach_friction.default_roughness = 0.0
+    rundata.landspill_data.darcy_weisbach_friction.filename = "roughness.txt"
+
+    # add a hydrological feature raster file
+    rundata.landspill_data.hydro_features.files.append("../common-files/utah-structures-hydro.asc")
+
+    # Evaporation
+    rundata.landspill_data.evaporation.type = 1
+    rundata.landspill_data.evaporation.coefficients = [1.38, 0.045]
+
+    return rundata
diff --git a/cases/utah-structures/setrun.py b/cases/utah-structures/setrun.py
deleted file mode 100644
index 0072454..0000000
--- a/cases/utah-structures/setrun.py
+++ /dev/null
@@ -1,305 +0,0 @@
-"""Module to set up run time parameters for Clawpack.
-
-The values set in the function setrun are then written out to data files
-that will be read in by the Fortran code.
-
-"""
-# pylint: disable=no-member, too-many-statements
-import functools
-from gclandspill.clawutil import data
-from gclandspill.data import LandSpillData
-
-
-def setrun(claw_pkg: str = 'geoclaw'):
-    """Define the parameters used for running Clawpack.
-
-    Arguments
-    ---------
-    claw_pkg : str
-        Expected to be "geoclaw" for this setrun.
-
-    Returns
-    -------
-    rundata : clawutil.ClawRunData
-    """
-
-    assert claw_pkg.lower() == 'geoclaw',  "Expected claw_pkg = 'geoclaw'"
-
-    num_dim = 2
-    rundata = data.ClawRunData(claw_pkg, num_dim)
-
-    # GeoClaw specific parameters:
-    rundata = setgeo(rundata)
-
-    # AmrClaw specific parameters:
-    rundata = setamr(rundata)
-
-    # an alias to object holding standard Clawpack parameters
-    clawdata = rundata.clawdata
-
-    # number of space dimensions:
-    clawdata.num_dim = num_dim
-
-    # lower and upper edge of computational domain:
-    clawdata.lower[0] = -12458815.605194421-1000.0
-    clawdata.upper[0] = -12458815.605194421+1000.0
-    clawdata.lower[1] = 4987307.342152054-1000.0
-    clawdata.upper[1] = 4987307.342152054+1000.0
-
-    # number of grid cells: coarsest grid
-    clawdata.num_cells[0] = 500
-    clawdata.num_cells[1] = 500
-
-    # number of equations in the system: continuity eq. and x, y momemtums
-    clawdata.num_eqn = 3
-
-    # number of auxiliary variables: elevation and hydrological data
-    clawdata.num_aux = 2
-
-    # index of aux array corresponding to capacity function, if there is one:
-    clawdata.capa_index = 0
-
-    # initial time:
-    clawdata.t0 = 0.0
-
-    # restart from checkpoint file of a previous run?
-    clawdata.restart = False               # True to restart from prior results
-    clawdata.restart_file = 'fort.chk00006'  # File to use for restart data
-
-    # temporal data output style
-    clawdata.output_style = 1
-
-    if clawdata.output_style == 1:  # equally spaced times up to tfinal
-        clawdata.num_output_times = 240
-        clawdata.tfinal = 28800
-        clawdata.output_t0 = True  # output at initial (or restart) time?
-    elif clawdata.output_style == 2:  # a list of output times
-        clawdata.output_times = [0., 1800.]
-    elif clawdata.output_style == 3:  # every iout timesteps with a total of ntot time steps:
-        clawdata.output_step_interval = 80
-        clawdata.total_steps = 4000
-        clawdata.output_t0 = True
-
-    # output data format
-    clawdata.output_format = 'binary'
-
-    # which veriable to output
-    clawdata.output_q_components = 'all'  # could be list such as [True,True]
-
-    # which aux data to output
-    clawdata.output_aux_components = 'all'  # could be list
-
-    # whether to output aux data only at t0
-    clawdata.output_aux_onlyonce = False
-
-    # print information up to this AMR level
-    clawdata.verbosity = 3
-
-    # use veriable time step size based on CFL
-    clawdata.dt_variable = 1
-
-    # data required to calculate a proper initial time step size
-    cell_area = \
-        (clawdata.upper[0] - clawdata.lower[0]) * (clawdata.upper[1] - clawdata.lower[1]) / \
-        (clawdata.num_cells[0] * clawdata.num_cells[1])
-    cell_area /= functools.reduce(
-        lambda x, y: x*y, rundata.amrdata.refinement_ratios_x[:rundata.amrdata.amr_levels_max], 1)
-    cell_area /= functools.reduce(
-        lambda x, y: x*y, rundata.amrdata.refinement_ratios_y[:rundata.amrdata.amr_levels_max], 1)
-    vrate = rundata.landspill_data.point_sources.point_sources[0][-1][0]
-
-    # initial time step for variable dt.
-    clawdata.dt_initial = 0.3 * cell_area / vrate
-
-    # max time step allowed if variable dt
-    clawdata.dt_max = 4.0
-
-    # desired CFL if variable dt, and max CFL allowed without retaking a step with a smaller dt
-    clawdata.cfl_desired = 0.9
-    clawdata.cfl_max = 0.95
-
-    # maximum number of time steps allowed between two output times:
-    clawdata.steps_max = 100000
-
-    # order of accuracy:  1 => Godunov,  2 => Lax-Wendroff plus limiters
-    clawdata.order = 2
-
-    # use dimensional splitting? (not yet available for AMR)
-    clawdata.dimensional_split = 'unsplit'
-
-    # 2 or 'all' => corner transport of 2nd order corrections too
-    clawdata.transverse_waves = 2
-
-    # number of waves in the Riemann solution:
-    clawdata.num_waves = 3
-
-    # list of limiters to use for each wave family:
-    clawdata.limiter = ['mc', 'mc', 'mc']
-
-    # use f-wave version of algorithms
-    clawdata.use_fwaves = True
-
-    # src_split == 1 or 'godunov' ==> Godunov (1st order) splitting used,
-    clawdata.source_split = 'godunov'
-
-    # number of ghost cells (usually 2)
-    clawdata.num_ghost = 2
-
-    # BC: 1 => extrapolation (non-reflecting outflow)
-    clawdata.bc_lower[0] = 1
-    clawdata.bc_upper[0] = 1
-    clawdata.bc_lower[1] = 1
-    clawdata.bc_upper[1] = 1
-
-    # checkout mechanism selection
-    clawdata.checkpt_style = 0
-
-    if clawdata.checkpt_style == 0:  # do not checkpoint at all
-        pass
-    elif clawdata.checkpt_style == 1:  # checkpoint only at tfinal.
-        pass
-    elif clawdata.checkpt_style == 2:  # specify a list of checkpoint times.
-        clawdata.checkpt_times = [0.1, 0.15]
-    elif clawdata.checkpt_style == 3:  # checkpoint every checkpt_interval timesteps (on Level 1)
-        clawdata.checkpt_interval = 5
-
-    return rundata
-
-
-def setamr(rundata: data.ClawRunData):
-    """Set AmrClaw specific runtime parameters."""
-
-    try:
-        amrdata = rundata.amrdata
-    except AttributeError as err:
-        raise AttributeError("This rundata has no amrdata attribute") from err
-
-    # max number of refinement levels:
-    amrdata.amr_levels_max = 2
-
-    # List of refinement ratios at each level (length at least mxnest-1)
-    amrdata.refinement_ratios_x = [4]
-    amrdata.refinement_ratios_y = [4]
-    amrdata.refinement_ratios_t = [4]
-
-    # aux variables are defined at cell centers
-    amrdata.aux_type = ['center', 'center']
-
-    # use flag2refine routine to flag refinement (rather than using richardson error)
-    amrdata.flag_richardson = False
-    amrdata.flag2refine = True
-
-    # steps to take on each level L between regriddings of level L+1:
-    amrdata.regrid_interval = 1
-
-    # width of buffer zone around flagged points
-    amrdata.regrid_buffer_width = 1
-
-    # clustering alg. cutoff for (# flagged pts) / (total # of cells refined)
-    amrdata.clustering_cutoff = 0.80000
-
-    # print info about each regridding up to this level:
-    amrdata.verbosity_regrid = 0
-
-    #  ----- For developers -----
-    # Toggle debugging print statements:
-    amrdata.dprint = False      # print domain flags
-    amrdata.eprint = False      # print err est flags
-    amrdata.edebug = False      # even more err est flags
-    amrdata.gprint = False      # grid bisection/clustering
-    amrdata.nprint = False      # proper nesting output
-    amrdata.pprint = False      # proj. of tagged points
-    amrdata.rprint = False      # print regridding summary
-    amrdata.sprint = False      # space/memory output
-    amrdata.tprint = False      # time step reporting each level
-    amrdata.uprint = False      # update/upbnd reporting
-
-    # customize regions of refinement
-    # rundata.regiondata.regions.append([minlevel, maxlevel, t1, t2, x1, x2, y1, y2])
-
-    # gauges
-    # rundata.gaugedata.gauges.append([gaugeno, x, y, t1, t2])
-
-    return rundata
-
-
-def setgeo(rundata: data.ClawRunData):
-    """Set GeoClaw specific runtime parameters."""
-
-    try:
-        geo_data = rundata.geo_data
-    except AttributeError as err:
-        raise AttributeError("This rundata has no geo_data attribute") from err
-
-    # == Physics ==
-    geo_data.gravity = 9.81
-    geo_data.coordinate_system = 1
-    geo_data.earth_radius = 6367.5e3
-
-    # == Forcing Options
-    geo_data.coriolis_forcing = False
-
-    # == Algorithm and Initial Conditions ==
-    geo_data.sea_level = 0.0
-    geo_data.dry_tolerance = 1.e-4
-    geo_data.friction_forcing = False
-    geo_data.manning_coefficient = 0.035
-    geo_data.friction_depth = 1.e6
-
-    # Refinement data
-    refinement_data = rundata.refinement_data
-    refinement_data.wave_tolerance = 1.e-5
-    refinement_data.speed_tolerance = [1e-8]
-    refinement_data.variable_dt_refinement_ratios = True
-
-    # for topography, append lines of the form [topotype, fname]
-    topo_data = rundata.topo_data
-    topo_data.topofiles.append([3, "../common-files/utah-structure.asc"])
-
-    # for moving topography, append lines of the form: [topotype, fname]
-    # dtopo_data = rundata.dtopo_data
-
-    # for qinit perturbations, append lines of the form: [fname]
-    rundata.qinit_data.qinit_type = 0
-    rundata.qinit_data.qinitfiles = []
-
-    # for fixed grids append lines of the form
-    # [t1, t2, noutput, x1, x2, y1, y2, xpoints, ypoints, ioutarrivaltimes, ioutsurfacemax]
-    # fixedgrids = rundata.fixed_grid_data
-
-    # Land-spill module settings
-    rundata.add_data(LandSpillData(), 'landspill_data')
-    landspill = rundata.landspill_data
-    landspill.ref_mu = 332.  # cP @ 15 degree C
-    landspill.ref_temperature = 15.
-    landspill.ambient_temperature = 25.
-    landspill.density = 9.266e2  # kg / m^3 @ 15 degree C; will overwrite rho in GeoClaw
-
-    # extra parameters
-    landspill.update_tol = geo_data.dry_tolerance
-    landspill.refine_tol = 0.0
-
-    # Point sources
-    ptsources_data = landspill.point_sources
-    ptsources_data.n_point_sources = 1
-    ptsources_data.point_sources.append(
-        [[-12458815.605194421, 4987307.342152054], 2, [ 1800., 12600.], [0.5, 0.1]])
-
-    # Darcy-Weisbach friction
-    darcy_weisbach_data = landspill.darcy_weisbach_friction
-    darcy_weisbach_data.type = 4
-    darcy_weisbach_data.dry_tol = 1e-4
-    darcy_weisbach_data.friction_tol = 1e6
-    darcy_weisbach_data.default_roughness = 0.0
-    darcy_weisbach_data.filename = "roughness.txt"
-
-    # hydrological features
-    hydro_feature_data = landspill.hydro_features
-    hydro_feature_data.files.append("../common-files/utah-structure-hydro.asc")
-
-    # Evaporation
-    evaporation_data = landspill.evaporation
-    evaporation_data.type = 1
-    evaporation_data.coefficients = [1.38, 0.045]
-
-    return rundata
diff --git a/cases/utah-waterbody-gasoline-no-evap/animation.gif b/cases/utah-waterbody-gasoline-no-evap/animation.gif
new file mode 100644
index 0000000..0eedbb4
Binary files /dev/null and b/cases/utah-waterbody-gasoline-no-evap/animation.gif differ
diff --git a/cases/utah-waterbody-gasoline-no-evap/setrun.py b/cases/utah-waterbody-gasoline-no-evap/setrun.py
index b464d02..52c6db6 100644
--- a/cases/utah-waterbody-gasoline-no-evap/setrun.py
+++ b/cases/utah-waterbody-gasoline-no-evap/setrun.py
@@ -2,304 +2,63 @@
 
 The values set in the function setrun are then written out to data files
 that will be read in by the Fortran code.
-
 """
 # pylint: disable=no-member, too-many-statements
-import functools
-from gclandspill.clawutil import data
-from gclandspill.data import LandSpillData
+import gclandspill.data
 
 
-def setrun(claw_pkg: str = 'geoclaw'):
+def setrun():
     """Define the parameters used for running Clawpack.
 
-    Arguments
-    ---------
-    claw_pkg : str
-        Expected to be "geoclaw" for this setrun.
-
     Returns
     -------
-    rundata : clawutil.ClawRunData
+    rundata : gclandspill.data.ClawRunData
     """
 
-    assert claw_pkg.lower() == 'geoclaw',  "Expected claw_pkg = 'geoclaw'"
-
-    num_dim = 2
-    rundata = data.ClawRunData(claw_pkg, num_dim)
-
-    # GeoClaw specific parameters:
-    rundata = setgeo(rundata)
+    rundata = gclandspill.data.ClawRunData()
 
-    # AmrClaw specific parameters:
-    rundata = setamr(rundata)
+    # for convenience
+    break_point = [-12460209.5, 4985137.4]
 
-    # an alias to object holding standard Clawpack parameters
-    clawdata = rundata.clawdata
-
-    # number of space dimensions:
-    clawdata.num_dim = num_dim
-
-    # lower and upper edge of computational domain:
-    clawdata.lower[0] = -12460209.5-400.
-    clawdata.upper[0] = -12460209.5+400.
-    clawdata.lower[1] = 4985137.4-400
-    clawdata.upper[1] = 4985137.4+400
+    # lower and upper edges of computational domain
+    rundata.clawdata.lower[0] = break_point[0]-400.
+    rundata.clawdata.upper[0] = break_point[0]+400.
+    rundata.clawdata.lower[1] = break_point[1]-400
+    rundata.clawdata.upper[1] = break_point[1]+400
 
     # number of grid cells: coarsest grid
-    clawdata.num_cells[0] = 200
-    clawdata.num_cells[1] = 200
-
-    # number of equations in the system: continuity eq. and x, y momemtums
-    clawdata.num_eqn = 3
-
-    # number of auxiliary variables: elevation and hydrological data
-    clawdata.num_aux = 2
-
-    # index of aux array corresponding to capacity function, if there is one:
-    clawdata.capa_index = 0
-
-    # initial time:
-    clawdata.t0 = 0.0
-
-    # restart from checkpoint file of a previous run?
-    clawdata.restart = False               # True to restart from prior results
-    clawdata.restart_file = 'fort.chk00006'  # File to use for restart data
-
-    # temporal data output style
-    clawdata.output_style = 1
-
-    if clawdata.output_style == 1:  # equally spaced times up to tfinal
-        clawdata.num_output_times = 240
-        clawdata.tfinal = 28800
-        clawdata.output_t0 = True  # output at initial (or restart) time?
-    elif clawdata.output_style == 2:  # a list of output times
-        clawdata.output_times = [0., 1800.]
-    elif clawdata.output_style == 3:  # every iout timesteps with a total of ntot time steps:
-        clawdata.output_step_interval = 80
-        clawdata.total_steps = 4000
-        clawdata.output_t0 = True
-
-    # output data format
-    clawdata.output_format = 'binary'
-
-    # which veriable to output
-    clawdata.output_q_components = 'all'  # could be list such as [True,True]
-
-    # which aux data to output
-    clawdata.output_aux_components = 'all'  # could be list
-
-    # whether to output aux data only at t0
-    clawdata.output_aux_onlyonce = False
-
-    # print information up to this AMR level
-    clawdata.verbosity = 3
-
-    # use veriable time step size based on CFL
-    clawdata.dt_variable = 1
-
-    # data required to calculate a proper initial time step size
-    cell_area = \
-        (clawdata.upper[0] - clawdata.lower[0]) * (clawdata.upper[1] - clawdata.lower[1]) / \
-        (clawdata.num_cells[0] * clawdata.num_cells[1])
-    cell_area /= functools.reduce(
-        lambda x, y: x*y, rundata.amrdata.refinement_ratios_x[:rundata.amrdata.amr_levels_max], 1)
-    cell_area /= functools.reduce(
-        lambda x, y: x*y, rundata.amrdata.refinement_ratios_y[:rundata.amrdata.amr_levels_max], 1)
-    vrate = rundata.landspill_data.point_sources.point_sources[0][-1][0]
-
-    # initial time step for variable dt.
-    clawdata.dt_initial = 0.3 * cell_area / vrate
-
-    # max time step allowed if variable dt
-    clawdata.dt_max = 4.0
-
-    # desired CFL if variable dt, and max CFL allowed without retaking a step with a smaller dt
-    clawdata.cfl_desired = 0.9
-    clawdata.cfl_max = 0.95
-
-    # maximum number of time steps allowed between two output times:
-    clawdata.steps_max = 100000
-
-    # order of accuracy:  1 => Godunov,  2 => Lax-Wendroff plus limiters
-    clawdata.order = 2
-
-    # use dimensional splitting? (not yet available for AMR)
-    clawdata.dimensional_split = 'unsplit'
-
-    # 2 or 'all' => corner transport of 2nd order corrections too
-    clawdata.transverse_waves = 2
-
-    # number of waves in the Riemann solution:
-    clawdata.num_waves = 3
-
-    # list of limiters to use for each wave family:
-    clawdata.limiter = ['mc', 'mc', 'mc']
-
-    # use f-wave version of algorithms
-    clawdata.use_fwaves = True
-
-    # src_split == 1 or 'godunov' ==> Godunov (1st order) splitting used,
-    clawdata.source_split = 'godunov'
-
-    # number of ghost cells (usually 2)
-    clawdata.num_ghost = 2
+    rundata.clawdata.num_cells[0] = 200
+    rundata.clawdata.num_cells[1] = 200
 
-    # BC: 1 => extrapolation (non-reflecting outflow)
-    clawdata.bc_lower[0] = 1
-    clawdata.bc_upper[0] = 1
-    clawdata.bc_lower[1] = 1
-    clawdata.bc_upper[1] = 1
+    # temporal data output
+    rundata.clawdata.output_style = 1
+    rundata.clawdata.num_output_times = 240
+    rundata.clawdata.tfinal = 28800
+    rundata.clawdata.output_t0 = True  # output at initial (or restart) time?
 
-    # checkout mechanism selection
-    clawdata.checkpt_style = 0
-
-    if clawdata.checkpt_style == 0:  # do not checkpoint at all
-        pass
-    elif clawdata.checkpt_style == 1:  # checkpoint only at tfinal.
-        pass
-    elif clawdata.checkpt_style == 2:  # specify a list of checkpoint times.
-        clawdata.checkpt_times = [0.1, 0.15]
-    elif clawdata.checkpt_style == 3:  # checkpoint every checkpt_interval timesteps (on Level 1)
-        clawdata.checkpt_interval = 5
-
-    return rundata
-
-
-def setamr(rundata: data.ClawRunData):
-    """Set AmrClaw specific runtime parameters."""
-
-    try:
-        amrdata = rundata.amrdata
-    except AttributeError as err:
-        raise AttributeError("This rundata has no amrdata attribute") from err
-
-    # max number of refinement levels:
-    amrdata.amr_levels_max = 2
-
-    # List of refinement ratios at each level (length at least mxnest-1)
-    amrdata.refinement_ratios_x = [4]
-    amrdata.refinement_ratios_y = [4]
-    amrdata.refinement_ratios_t = [4]
-
-    # aux variables are defined at cell centers
-    amrdata.aux_type = ['center', 'center']
-
-    # use flag2refine routine to flag refinement (rather than using richardson error)
-    amrdata.flag_richardson = False
-    amrdata.flag2refine = True
-
-    # steps to take on each level L between regriddings of level L+1:
-    amrdata.regrid_interval = 1
-
-    # width of buffer zone around flagged points
-    amrdata.regrid_buffer_width = 1
-
-    # clustering alg. cutoff for (# flagged pts) / (total # of cells refined)
-    amrdata.clustering_cutoff = 0.80000
-
-    # print info about each regridding up to this level:
-    amrdata.verbosity_regrid = 0
-
-    #  ----- For developers -----
-    # Toggle debugging print statements:
-    amrdata.dprint = False      # print domain flags
-    amrdata.eprint = False      # print err est flags
-    amrdata.edebug = False      # even more err est flags
-    amrdata.gprint = False      # grid bisection/clustering
-    amrdata.nprint = False      # proper nesting output
-    amrdata.pprint = False      # proj. of tagged points
-    amrdata.rprint = False      # print regridding summary
-    amrdata.sprint = False      # space/memory output
-    amrdata.tprint = False      # time step reporting each level
-    amrdata.uprint = False      # update/upbnd reporting
-
-    # customize regions of refinement
-    # rundata.regiondata.regions.append([minlevel, maxlevel, t1, t2, x1, x2, y1, y2])
-
-    # gauges
-    # rundata.gaugedata.gauges.append([gaugeno, x, y, t1, t2])
-
-    return rundata
-
-
-def setgeo(rundata: data.ClawRunData):
-    """Set GeoClaw specific runtime parameters."""
-
-    try:
-        geo_data = rundata.geo_data
-    except AttributeError as err:
-        raise AttributeError("This rundata has no geo_data attribute") from err
-
-    # == Physics ==
-    geo_data.gravity = 9.81
-    geo_data.coordinate_system = 1
-    geo_data.earth_radius = 6367.5e3
-
-    # == Forcing Options
-    geo_data.coriolis_forcing = False
-
-    # == Algorithm and Initial Conditions ==
-    geo_data.sea_level = 0.0
-    geo_data.dry_tolerance = 1.e-4
-    geo_data.friction_forcing = False
-    geo_data.manning_coefficient = 0.035
-    geo_data.friction_depth = 1.e6
-
-    # Refinement data
-    refinement_data = rundata.refinement_data
-    refinement_data.wave_tolerance = 1.e-5
-    refinement_data.speed_tolerance = [1e-8]
-    refinement_data.variable_dt_refinement_ratios = True
-
-    # for topography, append lines of the form [topotype, fname]
-    topo_data = rundata.topo_data
-    topo_data.topofiles.append([3, '../common-files/utah-waterbody.asc'])
-
-    # for moving topography, append lines of the form: [topotype, fname]
-    # dtopo_data = rundata.dtopo_data
-
-    # for qinit perturbations, append lines of the form: [fname]
-    rundata.qinit_data.qinit_type = 0
-    rundata.qinit_data.qinitfiles = []
-
-    # for fixed grids append lines of the form
-    # [t1, t2, noutput, x1, x2, y1, y2, xpoints, ypoints, ioutarrivaltimes, ioutsurfacemax]
-    # fixedgrids = rundata.fixed_grid_data
+    # add a topography file; format [topotype, fname]
+    rundata.topo_data.topofiles.append([3, '../common-files/utah-waterbody.asc'])
 
     # Land-spill module settings
-    rundata.add_data(LandSpillData(), 'landspill_data')
-    landspill = rundata.landspill_data
-    landspill.ref_mu = 0.6512  # cP @ 15 degree C
-    landspill.ref_temperature = 15.
-    landspill.ambient_temperature = 25.
-    landspill.density = 800.  # kg / m^3 @ 15 degree C; will overwrite rho in GeoClaw
-
-    # extra parameters
-    landspill.update_tol = geo_data.dry_tolerance
-    landspill.refine_tol = 0.0
+    rundata.landspill_data.ref_mu = 0.6512  # cP @ 15 degree C
+    rundata.landspill_data.ref_temperature = 15.
+    rundata.landspill_data.ambient_temperature = 25.
+    rundata.landspill_data.density = 800.  # kg / m^3 @ 15 degree C; will overwrite rho in GeoClaw
 
     # Point sources
-    ptsources_data = landspill.point_sources
-    ptsources_data.n_point_sources = 1
-    ptsources_data.point_sources.append(
-        [[-12460209.5, 4985137.4], 2, [1800., 12600.], [0.5, 0.1]])
+    rundata.landspill_data.point_sources.n_point_sources = 1
+    rundata.landspill_data.point_sources.point_sources.append([break_point, 2, [1800., 12600.], [0.5, 0.1]])
 
     # Darcy-Weisbach friction
-    darcy_weisbach_data = landspill.darcy_weisbach_friction
-    darcy_weisbach_data.type = 4
-    darcy_weisbach_data.dry_tol = 1e-4
-    darcy_weisbach_data.friction_tol = 1e6
-    darcy_weisbach_data.default_roughness = 0.0
-    darcy_weisbach_data.filename = "roughness.txt"
+    rundata.landspill_data.darcy_weisbach_friction.type = 4
+    rundata.landspill_data.darcy_weisbach_friction.default_roughness = 0.0
+    rundata.landspill_data.darcy_weisbach_friction.filename = "roughness.txt"
 
-    # hydrological features
-    hydro_feature_data = landspill.hydro_features
-    hydro_feature_data.files.append("../common-files/utah-waterbody-hydro.asc")
+    # add a hydrological feature raster file
+    rundata.landspill_data.hydro_features.files.append("../common-files/utah-waterbody-hydro.asc")
 
     # Evaporation
-    evaporation_data = landspill.evaporation
-    evaporation_data.type = 0
-    evaporation_data.coefficients = [13.2, 0.21]
+    rundata.landspill_data.evaporation.type = 0
+    rundata.landspill_data.evaporation.coefficients = [13.2, 0.21]
 
     return rundata
diff --git a/cases/utah-waterbody-gasoline/animation.gif b/cases/utah-waterbody-gasoline/animation.gif
new file mode 100644
index 0000000..ea4c37a
Binary files /dev/null and b/cases/utah-waterbody-gasoline/animation.gif differ
diff --git a/cases/utah-waterbody-gasoline/setrun.py b/cases/utah-waterbody-gasoline/setrun.py
index 54d18b8..cd2164c 100644
--- a/cases/utah-waterbody-gasoline/setrun.py
+++ b/cases/utah-waterbody-gasoline/setrun.py
@@ -2,304 +2,63 @@
 
 The values set in the function setrun are then written out to data files
 that will be read in by the Fortran code.
-
 """
 # pylint: disable=no-member, too-many-statements
-import functools
-from gclandspill.clawutil import data
-from gclandspill.data import LandSpillData
+import gclandspill.data
 
 
-def setrun(claw_pkg: str = 'geoclaw'):
+def setrun():
     """Define the parameters used for running Clawpack.
 
-    Arguments
-    ---------
-    claw_pkg : str
-        Expected to be "geoclaw" for this setrun.
-
     Returns
     -------
-    rundata : clawutil.ClawRunData
+    rundata : gclandspill.data.ClawRunData
     """
 
-    assert claw_pkg.lower() == 'geoclaw',  "Expected claw_pkg = 'geoclaw'"
-
-    num_dim = 2
-    rundata = data.ClawRunData(claw_pkg, num_dim)
-
-    # GeoClaw specific parameters:
-    rundata = setgeo(rundata)
+    rundata = gclandspill.data.ClawRunData()
 
-    # AmrClaw specific parameters:
-    rundata = setamr(rundata)
+    # for convenience
+    break_point = [-12460209.5, 4985137.4]
 
-    # an alias to object holding standard Clawpack parameters
-    clawdata = rundata.clawdata
-
-    # number of space dimensions:
-    clawdata.num_dim = num_dim
-
-    # lower and upper edge of computational domain:
-    clawdata.lower[0] = -12460209.5-400.
-    clawdata.upper[0] = -12460209.5+400.
-    clawdata.lower[1] = 4985137.4-400
-    clawdata.upper[1] = 4985137.4+400
+    # lower and upper edges of computational domain
+    rundata.clawdata.lower[0] = break_point[0]-400.
+    rundata.clawdata.upper[0] = break_point[0]+400.
+    rundata.clawdata.lower[1] = break_point[1]-400
+    rundata.clawdata.upper[1] = break_point[1]+400
 
     # number of grid cells: coarsest grid
-    clawdata.num_cells[0] = 200
-    clawdata.num_cells[1] = 200
-
-    # number of equations in the system: continuity eq. and x, y momemtums
-    clawdata.num_eqn = 3
-
-    # number of auxiliary variables: elevation and hydrological data
-    clawdata.num_aux = 2
-
-    # index of aux array corresponding to capacity function, if there is one:
-    clawdata.capa_index = 0
-
-    # initial time:
-    clawdata.t0 = 0.0
-
-    # restart from checkpoint file of a previous run?
-    clawdata.restart = False               # True to restart from prior results
-    clawdata.restart_file = 'fort.chk00006'  # File to use for restart data
-
-    # temporal data output style
-    clawdata.output_style = 1
-
-    if clawdata.output_style == 1:  # equally spaced times up to tfinal
-        clawdata.num_output_times = 240
-        clawdata.tfinal = 28800
-        clawdata.output_t0 = True  # output at initial (or restart) time?
-    elif clawdata.output_style == 2:  # a list of output times
-        clawdata.output_times = [0., 1800.]
-    elif clawdata.output_style == 3:  # every iout timesteps with a total of ntot time steps:
-        clawdata.output_step_interval = 80
-        clawdata.total_steps = 4000
-        clawdata.output_t0 = True
-
-    # output data format
-    clawdata.output_format = 'binary'
-
-    # which veriable to output
-    clawdata.output_q_components = 'all'  # could be list such as [True,True]
-
-    # which aux data to output
-    clawdata.output_aux_components = 'all'  # could be list
-
-    # whether to output aux data only at t0
-    clawdata.output_aux_onlyonce = False
-
-    # print information up to this AMR level
-    clawdata.verbosity = 3
-
-    # use veriable time step size based on CFL
-    clawdata.dt_variable = 1
-
-    # data required to calculate a proper initial time step size
-    cell_area = \
-        (clawdata.upper[0] - clawdata.lower[0]) * (clawdata.upper[1] - clawdata.lower[1]) / \
-        (clawdata.num_cells[0] * clawdata.num_cells[1])
-    cell_area /= functools.reduce(
-        lambda x, y: x*y, rundata.amrdata.refinement_ratios_x[:rundata.amrdata.amr_levels_max], 1)
-    cell_area /= functools.reduce(
-        lambda x, y: x*y, rundata.amrdata.refinement_ratios_y[:rundata.amrdata.amr_levels_max], 1)
-    vrate = rundata.landspill_data.point_sources.point_sources[0][-1][0]
-
-    # initial time step for variable dt.
-    clawdata.dt_initial = 0.3 * cell_area / vrate
-
-    # max time step allowed if variable dt
-    clawdata.dt_max = 4.0
-
-    # desired CFL if variable dt, and max CFL allowed without retaking a step with a smaller dt
-    clawdata.cfl_desired = 0.9
-    clawdata.cfl_max = 0.95
-
-    # maximum number of time steps allowed between two output times:
-    clawdata.steps_max = 100000
-
-    # order of accuracy:  1 => Godunov,  2 => Lax-Wendroff plus limiters
-    clawdata.order = 2
-
-    # use dimensional splitting? (not yet available for AMR)
-    clawdata.dimensional_split = 'unsplit'
-
-    # 2 or 'all' => corner transport of 2nd order corrections too
-    clawdata.transverse_waves = 2
-
-    # number of waves in the Riemann solution:
-    clawdata.num_waves = 3
-
-    # list of limiters to use for each wave family:
-    clawdata.limiter = ['mc', 'mc', 'mc']
-
-    # use f-wave version of algorithms
-    clawdata.use_fwaves = True
-
-    # src_split == 1 or 'godunov' ==> Godunov (1st order) splitting used,
-    clawdata.source_split = 'godunov'
-
-    # number of ghost cells (usually 2)
-    clawdata.num_ghost = 2
+    rundata.clawdata.num_cells[0] = 200
+    rundata.clawdata.num_cells[1] = 200
 
-    # BC: 1 => extrapolation (non-reflecting outflow)
-    clawdata.bc_lower[0] = 1
-    clawdata.bc_upper[0] = 1
-    clawdata.bc_lower[1] = 1
-    clawdata.bc_upper[1] = 1
+    # temporal data output
+    rundata.clawdata.output_style = 1
+    rundata.clawdata.num_output_times = 240
+    rundata.clawdata.tfinal = 28800
+    rundata.clawdata.output_t0 = True  # output at initial (or restart) time?
 
-    # checkout mechanism selection
-    clawdata.checkpt_style = 0
-
-    if clawdata.checkpt_style == 0:  # do not checkpoint at all
-        pass
-    elif clawdata.checkpt_style == 1:  # checkpoint only at tfinal.
-        pass
-    elif clawdata.checkpt_style == 2:  # specify a list of checkpoint times.
-        clawdata.checkpt_times = [0.1, 0.15]
-    elif clawdata.checkpt_style == 3:  # checkpoint every checkpt_interval timesteps (on Level 1)
-        clawdata.checkpt_interval = 5
-
-    return rundata
-
-
-def setamr(rundata: data.ClawRunData):
-    """Set AmrClaw specific runtime parameters."""
-
-    try:
-        amrdata = rundata.amrdata
-    except AttributeError as err:
-        raise AttributeError("This rundata has no amrdata attribute") from err
-
-    # max number of refinement levels:
-    amrdata.amr_levels_max = 2
-
-    # List of refinement ratios at each level (length at least mxnest-1)
-    amrdata.refinement_ratios_x = [4]
-    amrdata.refinement_ratios_y = [4]
-    amrdata.refinement_ratios_t = [4]
-
-    # aux variables are defined at cell centers
-    amrdata.aux_type = ['center', 'center']
-
-    # use flag2refine routine to flag refinement (rather than using richardson error)
-    amrdata.flag_richardson = False
-    amrdata.flag2refine = True
-
-    # steps to take on each level L between regriddings of level L+1:
-    amrdata.regrid_interval = 1
-
-    # width of buffer zone around flagged points
-    amrdata.regrid_buffer_width = 1
-
-    # clustering alg. cutoff for (# flagged pts) / (total # of cells refined)
-    amrdata.clustering_cutoff = 0.80000
-
-    # print info about each regridding up to this level:
-    amrdata.verbosity_regrid = 0
-
-    #  ----- For developers -----
-    # Toggle debugging print statements:
-    amrdata.dprint = False      # print domain flags
-    amrdata.eprint = False      # print err est flags
-    amrdata.edebug = False      # even more err est flags
-    amrdata.gprint = False      # grid bisection/clustering
-    amrdata.nprint = False      # proper nesting output
-    amrdata.pprint = False      # proj. of tagged points
-    amrdata.rprint = False      # print regridding summary
-    amrdata.sprint = False      # space/memory output
-    amrdata.tprint = False      # time step reporting each level
-    amrdata.uprint = False      # update/upbnd reporting
-
-    # customize regions of refinement
-    # rundata.regiondata.regions.append([minlevel, maxlevel, t1, t2, x1, x2, y1, y2])
-
-    # gauges
-    # rundata.gaugedata.gauges.append([gaugeno, x, y, t1, t2])
-
-    return rundata
-
-
-def setgeo(rundata: data.ClawRunData):
-    """Set GeoClaw specific runtime parameters."""
-
-    try:
-        geo_data = rundata.geo_data
-    except AttributeError as err:
-        raise AttributeError("This rundata has no geo_data attribute") from err
-
-    # == Physics ==
-    geo_data.gravity = 9.81
-    geo_data.coordinate_system = 1
-    geo_data.earth_radius = 6367.5e3
-
-    # == Forcing Options
-    geo_data.coriolis_forcing = False
-
-    # == Algorithm and Initial Conditions ==
-    geo_data.sea_level = 0.0
-    geo_data.dry_tolerance = 1.e-4
-    geo_data.friction_forcing = False
-    geo_data.manning_coefficient = 0.035
-    geo_data.friction_depth = 1.e6
-
-    # Refinement data
-    refinement_data = rundata.refinement_data
-    refinement_data.wave_tolerance = 1.e-5
-    refinement_data.speed_tolerance = [1e-8]
-    refinement_data.variable_dt_refinement_ratios = True
-
-    # for topography, append lines of the form [topotype, fname]
-    topo_data = rundata.topo_data
-    topo_data.topofiles.append([3, '../common-files/utah-waterbody.asc'])
-
-    # for moving topography, append lines of the form: [topotype, fname]
-    # dtopo_data = rundata.dtopo_data
-
-    # for qinit perturbations, append lines of the form: [fname]
-    rundata.qinit_data.qinit_type = 0
-    rundata.qinit_data.qinitfiles = []
-
-    # for fixed grids append lines of the form
-    # [t1, t2, noutput, x1, x2, y1, y2, xpoints, ypoints, ioutarrivaltimes, ioutsurfacemax]
-    # fixedgrids = rundata.fixed_grid_data
+    # add a topography file; format [topotype, fname]
+    rundata.topo_data.topofiles.append([3, '../common-files/utah-waterbody.asc'])
 
     # Land-spill module settings
-    rundata.add_data(LandSpillData(), 'landspill_data')
-    landspill = rundata.landspill_data
-    landspill.ref_mu = 0.6512  # cP @ 15 degree C
-    landspill.ref_temperature = 15.
-    landspill.ambient_temperature = 25.
-    landspill.density = 800.  # kg / m^3 @ 15 degree C; will overwrite rho in GeoClaw
-
-    # extra parameters
-    landspill.update_tol = geo_data.dry_tolerance
-    landspill.refine_tol = 0.0
+    rundata.landspill_data.ref_mu = 0.6512  # cP @ 15 degree C
+    rundata.landspill_data.ref_temperature = 15.
+    rundata.landspill_data.ambient_temperature = 25.
+    rundata.landspill_data.density = 800.  # kg / m^3 @ 15 degree C; will overwrite rho in GeoClaw
 
     # Point sources
-    ptsources_data = landspill.point_sources
-    ptsources_data.n_point_sources = 1
-    ptsources_data.point_sources.append(
-        [[-12460209.5, 4985137.4], 2, [1800., 12600.], [0.5, 0.1]])
+    rundata.landspill_data.point_sources.n_point_sources = 1
+    rundata.landspill_data.point_sources.point_sources.append([break_point, 2, [1800., 12600.], [0.5, 0.1]])
 
     # Darcy-Weisbach friction
-    darcy_weisbach_data = landspill.darcy_weisbach_friction
-    darcy_weisbach_data.type = 4
-    darcy_weisbach_data.dry_tol = 1e-4
-    darcy_weisbach_data.friction_tol = 1e6
-    darcy_weisbach_data.default_roughness = 0.0
-    darcy_weisbach_data.filename = "roughness.txt"
+    rundata.landspill_data.darcy_weisbach_friction.type = 4
+    rundata.landspill_data.darcy_weisbach_friction.default_roughness = 0.0
+    rundata.landspill_data.darcy_weisbach_friction.filename = "roughness.txt"
 
-    # hydrological features
-    hydro_feature_data = landspill.hydro_features
-    hydro_feature_data.files.append("../common-files/utah-waterbody-hydro.asc")
+    # add a hydrological feature raster file
+    rundata.landspill_data.hydro_features.files.append("../common-files/utah-waterbody-hydro.asc")
 
     # Evaporation
-    evaporation_data = landspill.evaporation
-    evaporation_data.type = 1
-    evaporation_data.coefficients = [13.2, 0.21]
+    rundata.landspill_data.evaporation.type = 1
+    rundata.landspill_data.evaporation.coefficients = [13.2, 0.21]
 
     return rundata
diff --git a/cases/utah-waterbody-maya-no-evap/animation.gif b/cases/utah-waterbody-maya-no-evap/animation.gif
new file mode 100644
index 0000000..1903626
Binary files /dev/null and b/cases/utah-waterbody-maya-no-evap/animation.gif differ
diff --git a/cases/utah-waterbody-maya-no-evap/setrun.py b/cases/utah-waterbody-maya-no-evap/setrun.py
index aac29c1..1b54acc 100644
--- a/cases/utah-waterbody-maya-no-evap/setrun.py
+++ b/cases/utah-waterbody-maya-no-evap/setrun.py
@@ -2,304 +2,63 @@
 
 The values set in the function setrun are then written out to data files
 that will be read in by the Fortran code.
-
 """
 # pylint: disable=no-member, too-many-statements
-import functools
-from gclandspill.clawutil import data
-from gclandspill.data import LandSpillData
+import gclandspill.data
 
 
-def setrun(claw_pkg: str = 'geoclaw'):
+def setrun():
     """Define the parameters used for running Clawpack.
 
-    Arguments
-    ---------
-    claw_pkg : str
-        Expected to be "geoclaw" for this setrun.
-
     Returns
     -------
-    rundata : clawutil.ClawRunData
+    rundata : gclandspill.data.ClawRunData
     """
 
-    assert claw_pkg.lower() == 'geoclaw',  "Expected claw_pkg = 'geoclaw'"
-
-    num_dim = 2
-    rundata = data.ClawRunData(claw_pkg, num_dim)
-
-    # GeoClaw specific parameters:
-    rundata = setgeo(rundata)
+    rundata = gclandspill.data.ClawRunData()
 
-    # AmrClaw specific parameters:
-    rundata = setamr(rundata)
+    # for convenience
+    break_point = [-12460209.5, 4985137.4]
 
-    # an alias to object holding standard Clawpack parameters
-    clawdata = rundata.clawdata
-
-    # number of space dimensions:
-    clawdata.num_dim = num_dim
-
-    # lower and upper edge of computational domain:
-    clawdata.lower[0] = -12460209.5-400.
-    clawdata.upper[0] = -12460209.5+400.
-    clawdata.lower[1] = 4985137.4-400
-    clawdata.upper[1] = 4985137.4+400
+    # lower and upper edges of computational domain
+    rundata.clawdata.lower[0] = break_point[0]-400.
+    rundata.clawdata.upper[0] = break_point[0]+400.
+    rundata.clawdata.lower[1] = break_point[1]-400
+    rundata.clawdata.upper[1] = break_point[1]+400
 
     # number of grid cells: coarsest grid
-    clawdata.num_cells[0] = 200
-    clawdata.num_cells[1] = 200
-
-    # number of equations in the system: continuity eq. and x, y momemtums
-    clawdata.num_eqn = 3
-
-    # number of auxiliary variables: elevation and hydrological data
-    clawdata.num_aux = 2
-
-    # index of aux array corresponding to capacity function, if there is one:
-    clawdata.capa_index = 0
-
-    # initial time:
-    clawdata.t0 = 0.0
-
-    # restart from checkpoint file of a previous run?
-    clawdata.restart = False               # True to restart from prior results
-    clawdata.restart_file = 'fort.chk00006'  # File to use for restart data
-
-    # temporal data output style
-    clawdata.output_style = 1
-
-    if clawdata.output_style == 1:  # equally spaced times up to tfinal
-        clawdata.num_output_times = 240
-        clawdata.tfinal = 28800
-        clawdata.output_t0 = True  # output at initial (or restart) time?
-    elif clawdata.output_style == 2:  # a list of output times
-        clawdata.output_times = [0., 1800.]
-    elif clawdata.output_style == 3:  # every iout timesteps with a total of ntot time steps:
-        clawdata.output_step_interval = 80
-        clawdata.total_steps = 4000
-        clawdata.output_t0 = True
-
-    # output data format
-    clawdata.output_format = 'binary'
-
-    # which veriable to output
-    clawdata.output_q_components = 'all'  # could be list such as [True,True]
-
-    # which aux data to output
-    clawdata.output_aux_components = 'all'  # could be list
-
-    # whether to output aux data only at t0
-    clawdata.output_aux_onlyonce = False
-
-    # print information up to this AMR level
-    clawdata.verbosity = 3
-
-    # use veriable time step size based on CFL
-    clawdata.dt_variable = 1
-
-    # data required to calculate a proper initial time step size
-    cell_area = \
-        (clawdata.upper[0] - clawdata.lower[0]) * (clawdata.upper[1] - clawdata.lower[1]) / \
-        (clawdata.num_cells[0] * clawdata.num_cells[1])
-    cell_area /= functools.reduce(
-        lambda x, y: x*y, rundata.amrdata.refinement_ratios_x[:rundata.amrdata.amr_levels_max], 1)
-    cell_area /= functools.reduce(
-        lambda x, y: x*y, rundata.amrdata.refinement_ratios_y[:rundata.amrdata.amr_levels_max], 1)
-    vrate = rundata.landspill_data.point_sources.point_sources[0][-1][0]
-
-    # initial time step for variable dt.
-    clawdata.dt_initial = 0.3 * cell_area / vrate
-
-    # max time step allowed if variable dt
-    clawdata.dt_max = 4.0
-
-    # desired CFL if variable dt, and max CFL allowed without retaking a step with a smaller dt
-    clawdata.cfl_desired = 0.9
-    clawdata.cfl_max = 0.95
-
-    # maximum number of time steps allowed between two output times:
-    clawdata.steps_max = 100000
-
-    # order of accuracy:  1 => Godunov,  2 => Lax-Wendroff plus limiters
-    clawdata.order = 2
-
-    # use dimensional splitting? (not yet available for AMR)
-    clawdata.dimensional_split = 'unsplit'
-
-    # 2 or 'all' => corner transport of 2nd order corrections too
-    clawdata.transverse_waves = 2
-
-    # number of waves in the Riemann solution:
-    clawdata.num_waves = 3
-
-    # list of limiters to use for each wave family:
-    clawdata.limiter = ['mc', 'mc', 'mc']
-
-    # use f-wave version of algorithms
-    clawdata.use_fwaves = True
-
-    # src_split == 1 or 'godunov' ==> Godunov (1st order) splitting used,
-    clawdata.source_split = 'godunov'
-
-    # number of ghost cells (usually 2)
-    clawdata.num_ghost = 2
+    rundata.clawdata.num_cells[0] = 200
+    rundata.clawdata.num_cells[1] = 200
 
-    # BC: 1 => extrapolation (non-reflecting outflow)
-    clawdata.bc_lower[0] = 1
-    clawdata.bc_upper[0] = 1
-    clawdata.bc_lower[1] = 1
-    clawdata.bc_upper[1] = 1
+    # temporal data output
+    rundata.clawdata.output_style = 1
+    rundata.clawdata.num_output_times = 240
+    rundata.clawdata.tfinal = 28800
+    rundata.clawdata.output_t0 = True  # output at initial (or restart) time?
 
-    # checkout mechanism selection
-    clawdata.checkpt_style = 0
-
-    if clawdata.checkpt_style == 0:  # do not checkpoint at all
-        pass
-    elif clawdata.checkpt_style == 1:  # checkpoint only at tfinal.
-        pass
-    elif clawdata.checkpt_style == 2:  # specify a list of checkpoint times.
-        clawdata.checkpt_times = [0.1, 0.15]
-    elif clawdata.checkpt_style == 3:  # checkpoint every checkpt_interval timesteps (on Level 1)
-        clawdata.checkpt_interval = 5
-
-    return rundata
-
-
-def setamr(rundata: data.ClawRunData):
-    """Set AmrClaw specific runtime parameters."""
-
-    try:
-        amrdata = rundata.amrdata
-    except AttributeError as err:
-        raise AttributeError("This rundata has no amrdata attribute") from err
-
-    # max number of refinement levels:
-    amrdata.amr_levels_max = 2
-
-    # List of refinement ratios at each level (length at least mxnest-1)
-    amrdata.refinement_ratios_x = [4]
-    amrdata.refinement_ratios_y = [4]
-    amrdata.refinement_ratios_t = [4]
-
-    # aux variables are defined at cell centers
-    amrdata.aux_type = ['center', 'center']
-
-    # use flag2refine routine to flag refinement (rather than using richardson error)
-    amrdata.flag_richardson = False
-    amrdata.flag2refine = True
-
-    # steps to take on each level L between regriddings of level L+1:
-    amrdata.regrid_interval = 1
-
-    # width of buffer zone around flagged points
-    amrdata.regrid_buffer_width = 1
-
-    # clustering alg. cutoff for (# flagged pts) / (total # of cells refined)
-    amrdata.clustering_cutoff = 0.80000
-
-    # print info about each regridding up to this level:
-    amrdata.verbosity_regrid = 0
-
-    #  ----- For developers -----
-    # Toggle debugging print statements:
-    amrdata.dprint = False      # print domain flags
-    amrdata.eprint = False      # print err est flags
-    amrdata.edebug = False      # even more err est flags
-    amrdata.gprint = False      # grid bisection/clustering
-    amrdata.nprint = False      # proper nesting output
-    amrdata.pprint = False      # proj. of tagged points
-    amrdata.rprint = False      # print regridding summary
-    amrdata.sprint = False      # space/memory output
-    amrdata.tprint = False      # time step reporting each level
-    amrdata.uprint = False      # update/upbnd reporting
-
-    # customize regions of refinement
-    # rundata.regiondata.regions.append([minlevel, maxlevel, t1, t2, x1, x2, y1, y2])
-
-    # gauges
-    # rundata.gaugedata.gauges.append([gaugeno, x, y, t1, t2])
-
-    return rundata
-
-
-def setgeo(rundata: data.ClawRunData):
-    """Set GeoClaw specific runtime parameters."""
-
-    try:
-        geo_data = rundata.geo_data
-    except AttributeError as err:
-        raise AttributeError("This rundata has no geo_data attribute") from err
-
-    # == Physics ==
-    geo_data.gravity = 9.81
-    geo_data.coordinate_system = 1
-    geo_data.earth_radius = 6367.5e3
-
-    # == Forcing Options
-    geo_data.coriolis_forcing = False
-
-    # == Algorithm and Initial Conditions ==
-    geo_data.sea_level = 0.0
-    geo_data.dry_tolerance = 1.e-4
-    geo_data.friction_forcing = False
-    geo_data.manning_coefficient = 0.035
-    geo_data.friction_depth = 1.e6
-
-    # Refinement data
-    refinement_data = rundata.refinement_data
-    refinement_data.wave_tolerance = 1.e-5
-    refinement_data.speed_tolerance = [1e-8]
-    refinement_data.variable_dt_refinement_ratios = True
-
-    # for topography, append lines of the form [topotype, fname]
-    topo_data = rundata.topo_data
-    topo_data.topofiles.append([3, '../common-files/utah-waterbody.asc'])
-
-    # for moving topography, append lines of the form: [topotype, fname]
-    # dtopo_data = rundata.dtopo_data
-
-    # for qinit perturbations, append lines of the form: [fname]
-    rundata.qinit_data.qinit_type = 0
-    rundata.qinit_data.qinitfiles = []
-
-    # for fixed grids append lines of the form
-    # [t1, t2, noutput, x1, x2, y1, y2, xpoints, ypoints, ioutarrivaltimes, ioutsurfacemax]
-    # fixedgrids = rundata.fixed_grid_data
+    # add a topography file; format [topotype, fname]
+    rundata.topo_data.topofiles.append([3, '../common-files/utah-waterbody.asc'])
 
     # Land-spill module settings
-    rundata.add_data(LandSpillData(), 'landspill_data')
-    landspill = rundata.landspill_data
-    landspill.ref_mu = 332.  # cP @ 15 degree C
-    landspill.ref_temperature = 15.
-    landspill.ambient_temperature = 25.
-    landspill.density = 9.266e2  # kg / m^3 @ 15 degree C; will overwrite rho in GeoClaw
-
-    # extra parameters
-    landspill.update_tol = geo_data.dry_tolerance
-    landspill.refine_tol = 0.0
+    rundata.landspill_data.ref_mu = 332.  # cP @ 15 degree C
+    rundata.landspill_data.ref_temperature = 15.
+    rundata.landspill_data.ambient_temperature = 25.
+    rundata.landspill_data.density = 9.266e2  # kg / m^3 @ 15 degree C; will overwrite rho in GeoClaw
 
     # Point sources
-    ptsources_data = landspill.point_sources
-    ptsources_data.n_point_sources = 1
-    ptsources_data.point_sources.append(
-        [[-12460209.5, 4985137.4], 2, [1800., 12600.], [0.5, 0.1]])
+    rundata.landspill_data.point_sources.n_point_sources = 1
+    rundata.landspill_data.point_sources.point_sources.append([break_point, 2, [1800., 12600.], [0.5, 0.1]])
 
     # Darcy-Weisbach friction
-    darcy_weisbach_data = landspill.darcy_weisbach_friction
-    darcy_weisbach_data.type = 4
-    darcy_weisbach_data.dry_tol = 1e-4
-    darcy_weisbach_data.friction_tol = 1e6
-    darcy_weisbach_data.default_roughness = 0.0
-    darcy_weisbach_data.filename = "roughness.txt"
+    rundata.landspill_data.darcy_weisbach_friction.type = 4
+    rundata.landspill_data.darcy_weisbach_friction.default_roughness = 0.0
+    rundata.landspill_data.darcy_weisbach_friction.filename = "roughness.txt"
 
-    # hydrological features
-    hydro_feature_data = landspill.hydro_features
-    hydro_feature_data.files.append("../common-files/utah-waterbody-hydro.asc")
+    # add a hydrological feature raster file
+    rundata.landspill_data.hydro_features.files.append("../common-files/utah-waterbody-hydro.asc")
 
     # Evaporation
-    evaporation_data = landspill.evaporation
-    evaporation_data.type = 0
-    evaporation_data.coefficients = [1.38, 0.045]
+    rundata.landspill_data.evaporation.type = 0
+    rundata.landspill_data.evaporation.coefficients = [1.38, 0.045]
 
     return rundata
diff --git a/cases/utah-waterbody-maya/animation.gif b/cases/utah-waterbody-maya/animation.gif
new file mode 100644
index 0000000..9e83611
Binary files /dev/null and b/cases/utah-waterbody-maya/animation.gif differ
diff --git a/cases/utah-waterbody-maya/setrun.py b/cases/utah-waterbody-maya/setrun.py
index 8780c29..f6f39c9 100644
--- a/cases/utah-waterbody-maya/setrun.py
+++ b/cases/utah-waterbody-maya/setrun.py
@@ -2,304 +2,63 @@
 
 The values set in the function setrun are then written out to data files
 that will be read in by the Fortran code.
-
 """
 # pylint: disable=no-member, too-many-statements
-import functools
-from gclandspill.clawutil import data
-from gclandspill.data import LandSpillData
+import gclandspill.data
 
 
-def setrun(claw_pkg: str = 'geoclaw'):
+def setrun():
     """Define the parameters used for running Clawpack.
 
-    Arguments
-    ---------
-    claw_pkg : str
-        Expected to be "geoclaw" for this setrun.
-
     Returns
     -------
-    rundata : clawutil.ClawRunData
+    rundata : gclandspill.data.ClawRunData
     """
 
-    assert claw_pkg.lower() == 'geoclaw',  "Expected claw_pkg = 'geoclaw'"
-
-    num_dim = 2
-    rundata = data.ClawRunData(claw_pkg, num_dim)
-
-    # GeoClaw specific parameters:
-    rundata = setgeo(rundata)
+    rundata = gclandspill.data.ClawRunData()
 
-    # AmrClaw specific parameters:
-    rundata = setamr(rundata)
+    # for convenience
+    break_point = [-12460209.5, 4985137.4]
 
-    # an alias to object holding standard Clawpack parameters
-    clawdata = rundata.clawdata
-
-    # number of space dimensions:
-    clawdata.num_dim = num_dim
-
-    # lower and upper edge of computational domain:
-    clawdata.lower[0] = -12460209.5-400.
-    clawdata.upper[0] = -12460209.5+400.
-    clawdata.lower[1] = 4985137.4-400
-    clawdata.upper[1] = 4985137.4+400
+    # lower and upper edges of computational domain
+    rundata.clawdata.lower[0] = break_point[0]-400.
+    rundata.clawdata.upper[0] = break_point[0]+400.
+    rundata.clawdata.lower[1] = break_point[1]-400
+    rundata.clawdata.upper[1] = break_point[1]+400
 
     # number of grid cells: coarsest grid
-    clawdata.num_cells[0] = 200
-    clawdata.num_cells[1] = 200
-
-    # number of equations in the system: continuity eq. and x, y momemtums
-    clawdata.num_eqn = 3
-
-    # number of auxiliary variables: elevation and hydrological data
-    clawdata.num_aux = 2
-
-    # index of aux array corresponding to capacity function, if there is one:
-    clawdata.capa_index = 0
-
-    # initial time:
-    clawdata.t0 = 0.0
-
-    # restart from checkpoint file of a previous run?
-    clawdata.restart = False               # True to restart from prior results
-    clawdata.restart_file = 'fort.chk00006'  # File to use for restart data
-
-    # temporal data output style
-    clawdata.output_style = 1
-
-    if clawdata.output_style == 1:  # equally spaced times up to tfinal
-        clawdata.num_output_times = 240
-        clawdata.tfinal = 28800
-        clawdata.output_t0 = True  # output at initial (or restart) time?
-    elif clawdata.output_style == 2:  # a list of output times
-        clawdata.output_times = [0., 1800.]
-    elif clawdata.output_style == 3:  # every iout timesteps with a total of ntot time steps:
-        clawdata.output_step_interval = 80
-        clawdata.total_steps = 4000
-        clawdata.output_t0 = True
-
-    # output data format
-    clawdata.output_format = 'binary'
-
-    # which veriable to output
-    clawdata.output_q_components = 'all'  # could be list such as [True,True]
-
-    # which aux data to output
-    clawdata.output_aux_components = 'all'  # could be list
-
-    # whether to output aux data only at t0
-    clawdata.output_aux_onlyonce = False
-
-    # print information up to this AMR level
-    clawdata.verbosity = 3
-
-    # use veriable time step size based on CFL
-    clawdata.dt_variable = 1
-
-    # data required to calculate a proper initial time step size
-    cell_area = \
-        (clawdata.upper[0] - clawdata.lower[0]) * (clawdata.upper[1] - clawdata.lower[1]) / \
-        (clawdata.num_cells[0] * clawdata.num_cells[1])
-    cell_area /= functools.reduce(
-        lambda x, y: x*y, rundata.amrdata.refinement_ratios_x[:rundata.amrdata.amr_levels_max], 1)
-    cell_area /= functools.reduce(
-        lambda x, y: x*y, rundata.amrdata.refinement_ratios_y[:rundata.amrdata.amr_levels_max], 1)
-    vrate = rundata.landspill_data.point_sources.point_sources[0][-1][0]
-
-    # initial time step for variable dt.
-    clawdata.dt_initial = 0.3 * cell_area / vrate
-
-    # max time step allowed if variable dt
-    clawdata.dt_max = 4.0
-
-    # desired CFL if variable dt, and max CFL allowed without retaking a step with a smaller dt
-    clawdata.cfl_desired = 0.9
-    clawdata.cfl_max = 0.95
-
-    # maximum number of time steps allowed between two output times:
-    clawdata.steps_max = 100000
-
-    # order of accuracy:  1 => Godunov,  2 => Lax-Wendroff plus limiters
-    clawdata.order = 2
-
-    # use dimensional splitting? (not yet available for AMR)
-    clawdata.dimensional_split = 'unsplit'
-
-    # 2 or 'all' => corner transport of 2nd order corrections too
-    clawdata.transverse_waves = 2
-
-    # number of waves in the Riemann solution:
-    clawdata.num_waves = 3
-
-    # list of limiters to use for each wave family:
-    clawdata.limiter = ['mc', 'mc', 'mc']
-
-    # use f-wave version of algorithms
-    clawdata.use_fwaves = True
-
-    # src_split == 1 or 'godunov' ==> Godunov (1st order) splitting used,
-    clawdata.source_split = 'godunov'
-
-    # number of ghost cells (usually 2)
-    clawdata.num_ghost = 2
+    rundata.clawdata.num_cells[0] = 200
+    rundata.clawdata.num_cells[1] = 200
 
-    # BC: 1 => extrapolation (non-reflecting outflow)
-    clawdata.bc_lower[0] = 1
-    clawdata.bc_upper[0] = 1
-    clawdata.bc_lower[1] = 1
-    clawdata.bc_upper[1] = 1
+    # temporal data output
+    rundata.clawdata.output_style = 1
+    rundata.clawdata.num_output_times = 240
+    rundata.clawdata.tfinal = 28800
+    rundata.clawdata.output_t0 = True  # output at initial (or restart) time?
 
-    # checkout mechanism selection
-    clawdata.checkpt_style = 0
-
-    if clawdata.checkpt_style == 0:  # do not checkpoint at all
-        pass
-    elif clawdata.checkpt_style == 1:  # checkpoint only at tfinal.
-        pass
-    elif clawdata.checkpt_style == 2:  # specify a list of checkpoint times.
-        clawdata.checkpt_times = [0.1, 0.15]
-    elif clawdata.checkpt_style == 3:  # checkpoint every checkpt_interval timesteps (on Level 1)
-        clawdata.checkpt_interval = 5
-
-    return rundata
-
-
-def setamr(rundata: data.ClawRunData):
-    """Set AmrClaw specific runtime parameters."""
-
-    try:
-        amrdata = rundata.amrdata
-    except AttributeError as err:
-        raise AttributeError("This rundata has no amrdata attribute") from err
-
-    # max number of refinement levels:
-    amrdata.amr_levels_max = 2
-
-    # List of refinement ratios at each level (length at least mxnest-1)
-    amrdata.refinement_ratios_x = [4]
-    amrdata.refinement_ratios_y = [4]
-    amrdata.refinement_ratios_t = [4]
-
-    # aux variables are defined at cell centers
-    amrdata.aux_type = ['center', 'center']
-
-    # use flag2refine routine to flag refinement (rather than using richardson error)
-    amrdata.flag_richardson = False
-    amrdata.flag2refine = True
-
-    # steps to take on each level L between regriddings of level L+1:
-    amrdata.regrid_interval = 1
-
-    # width of buffer zone around flagged points
-    amrdata.regrid_buffer_width = 1
-
-    # clustering alg. cutoff for (# flagged pts) / (total # of cells refined)
-    amrdata.clustering_cutoff = 0.80000
-
-    # print info about each regridding up to this level:
-    amrdata.verbosity_regrid = 0
-
-    #  ----- For developers -----
-    # Toggle debugging print statements:
-    amrdata.dprint = False      # print domain flags
-    amrdata.eprint = False      # print err est flags
-    amrdata.edebug = False      # even more err est flags
-    amrdata.gprint = False      # grid bisection/clustering
-    amrdata.nprint = False      # proper nesting output
-    amrdata.pprint = False      # proj. of tagged points
-    amrdata.rprint = False      # print regridding summary
-    amrdata.sprint = False      # space/memory output
-    amrdata.tprint = False      # time step reporting each level
-    amrdata.uprint = False      # update/upbnd reporting
-
-    # customize regions of refinement
-    # rundata.regiondata.regions.append([minlevel, maxlevel, t1, t2, x1, x2, y1, y2])
-
-    # gauges
-    # rundata.gaugedata.gauges.append([gaugeno, x, y, t1, t2])
-
-    return rundata
-
-
-def setgeo(rundata: data.ClawRunData):
-    """Set GeoClaw specific runtime parameters."""
-
-    try:
-        geo_data = rundata.geo_data
-    except AttributeError as err:
-        raise AttributeError("This rundata has no geo_data attribute") from err
-
-    # == Physics ==
-    geo_data.gravity = 9.81
-    geo_data.coordinate_system = 1
-    geo_data.earth_radius = 6367.5e3
-
-    # == Forcing Options
-    geo_data.coriolis_forcing = False
-
-    # == Algorithm and Initial Conditions ==
-    geo_data.sea_level = 0.0
-    geo_data.dry_tolerance = 1.e-4
-    geo_data.friction_forcing = False
-    geo_data.manning_coefficient = 0.035
-    geo_data.friction_depth = 1.e6
-
-    # Refinement data
-    refinement_data = rundata.refinement_data
-    refinement_data.wave_tolerance = 1.e-5
-    refinement_data.speed_tolerance = [1e-8]
-    refinement_data.variable_dt_refinement_ratios = True
-
-    # for topography, append lines of the form [topotype, fname]
-    topo_data = rundata.topo_data
-    topo_data.topofiles.append([3, '../common-files/utah-waterbody.asc'])
-
-    # for moving topography, append lines of the form: [topotype, fname]
-    # dtopo_data = rundata.dtopo_data
-
-    # for qinit perturbations, append lines of the form: [fname]
-    rundata.qinit_data.qinit_type = 0
-    rundata.qinit_data.qinitfiles = []
-
-    # for fixed grids append lines of the form
-    # [t1, t2, noutput, x1, x2, y1, y2, xpoints, ypoints, ioutarrivaltimes, ioutsurfacemax]
-    # fixedgrids = rundata.fixed_grid_data
+    # add a topography file; format [topotype, fname]
+    rundata.topo_data.topofiles.append([3, '../common-files/utah-waterbody.asc'])
 
     # Land-spill module settings
-    rundata.add_data(LandSpillData(), 'landspill_data')
-    landspill = rundata.landspill_data
-    landspill.ref_mu = 332.  # cP @ 15 degree C
-    landspill.ref_temperature = 15.
-    landspill.ambient_temperature = 25.
-    landspill.density = 9.266e2  # kg / m^3 @ 15 degree C; will overwrite rho in GeoClaw
-
-    # extra parameters
-    landspill.update_tol = geo_data.dry_tolerance
-    landspill.refine_tol = 0.0
+    rundata.landspill_data.ref_mu = 332.  # cP @ 15 degree C
+    rundata.landspill_data.ref_temperature = 15.
+    rundata.landspill_data.ambient_temperature = 25.
+    rundata.landspill_data.density = 9.266e2  # kg / m^3 @ 15 degree C; will overwrite rho in GeoClaw
 
     # Point sources
-    ptsources_data = landspill.point_sources
-    ptsources_data.n_point_sources = 1
-    ptsources_data.point_sources.append(
-        [[-12460209.5, 4985137.4], 2, [1800., 12600.], [0.5, 0.1]])
+    rundata.landspill_data.point_sources.n_point_sources = 1
+    rundata.landspill_data.point_sources.point_sources.append([break_point, 2, [1800., 12600.], [0.5, 0.1]])
 
     # Darcy-Weisbach friction
-    darcy_weisbach_data = landspill.darcy_weisbach_friction
-    darcy_weisbach_data.type = 4
-    darcy_weisbach_data.dry_tol = 1e-4
-    darcy_weisbach_data.friction_tol = 1e6
-    darcy_weisbach_data.default_roughness = 0.0
-    darcy_weisbach_data.filename = "roughness.txt"
+    rundata.landspill_data.darcy_weisbach_friction.type = 4
+    rundata.landspill_data.darcy_weisbach_friction.default_roughness = 0.0
+    rundata.landspill_data.darcy_weisbach_friction.filename = "roughness.txt"
 
-    # hydrological features
-    hydro_feature_data = landspill.hydro_features
-    hydro_feature_data.files.append("../common-files/utah-waterbody-hydro.asc")
+    # add a hydrological feature raster file
+    rundata.landspill_data.hydro_features.files.append("../common-files/utah-waterbody-hydro.asc")
 
     # Evaporation
-    evaporation_data = landspill.evaporation
-    evaporation_data.type = 1
-    evaporation_data.coefficients = [1.38, 0.045]
+    rundata.landspill_data.evaporation.type = 1
+    rundata.landspill_data.evaporation.coefficients = [1.38, 0.045]
 
     return rundata
diff --git a/codemeta.json b/codemeta.json
new file mode 100644
index 0000000..5415f68
--- /dev/null
+++ b/codemeta.json
@@ -0,0 +1,37 @@
+{
+  "@context": "https://raw.githubusercontent.com/codemeta/codemeta/master/codemeta.jsonld",
+  "@type": "Code",
+  "author": [
+    {
+      "@id": "https://orcid.org/0000-0001-6330-2709",
+      "@type": "Person",
+      "email": "pychuang@gwu.edu",
+      "name": "Pi-Yueh Chuang",
+      "affiliation": "Department of Mechanical and Aerospace Engineering, The George Washington University, Washington, DC, USA"
+    },
+    {
+      "@id": "",
+      "@type": "Person",
+      "email": "Tracy.Thorleifson@g2-is.com",
+      "name": "Tracy Thorleifson",
+      "affiliation": "G2 Integrated Solutions, Houston, TX, USA"
+    },
+    {
+      "@id": "http://orcid.org/0000-0001-5812-2711",
+      "@type": "Person",
+      "email": "labarba@gwu.edu",
+      "name": "Lorena A. Barba",
+      "affiliation": "Department of Mechanical and Aerospace Engineering, The George Washington University, Washington, DC, USA"
+    }
+  ],
+  "identifier": "",
+  "codeRepository": "https://github.com/barbagroup/geoclaw-landspill",
+  "datePublished": "2021-01-11",
+  "dateModified": "2021-01-11",
+  "dateCreated": "2021-01-11",
+  "description": "Oil land-spill and overland flow simulator for pipeline rupture events",
+  "keywords": "geoclaw, overland-flow, pipeline, shallow-water-equations, pipeline-ruptures, land-spill",
+  "license": "BSD 3-Clause",
+  "title": "geoclaw-landspill",
+  "version": "v1.0.dev3"
+}
diff --git a/doc/configuration.md b/doc/configuration.md
new file mode 100644
index 0000000..41917f6
--- /dev/null
+++ b/doc/configuration.md
@@ -0,0 +1,328 @@
+Configuration file: `setrun.py`
+===============================
+
+A simulation case folder must at least contains a `setrun.py`, which describes
+the configurations of a simulation. `setrun.py` defines a function called
+`setrun` that accepts no argument and returns an instance of
+`gclandspill.data.ClawRunData`. For those familiar with Clawpack's ecosystem,
+`gclandspill.data.ClawRunData` is a modified version of
+`clawpack.clawutil.data.ClawRunData`. Our modified `ClawRunData` has some
+default settings configured to meet *geoclaw-landspill*'s needs. Also, it has
+additional settings for landspill simulations.
+
+The basic working function `setrun` looks like:
+
+```python
+import gclandspill
+
+def setrun():
+    """Whatever docstring."""
+
+    rundata = gclandspill.data.ClawRunData()
+
+    <...setting up rundata...>
+
+    return rundata
+```
+
+For those unfamiliar with Clawpack/GeoClaw, roughly speaking, there are four
+parts in the `rundata` object that needs to be configured: core solver, AMR
+(adaptive mesh refinement), GeoClaw-specific, and landspill-specific
+configurations. The core solver, AMR, and GeoClaw configurations are described
+in the official Clawpack documentation:
+
+* Core solver: [Specifying classic run-time parameters in setrun.py](http://www.clawpack.org/setrun.html)
+* AMR: [Specifying AMRClaw run-time parameters in setrun.py](http://www.clawpack.org/setrun_amrclaw.html)
+* GeoClaw: [Specifying GeoClaw parameters in setrun.py](http://www.clawpack.org/setrun_geoclaw.html)
+
+This document covers only required settings and also settings specific for
+landspill. To better understand how a `setrun.py` looks like, please see the
+`setrun.py` in examples under the `cases` folder.
+
+--------------------------------------------
+## I. Settings in core solver, AMR, and GeoClaw
+
+Most settings in the core solver, AMR, and GeoClaw are tuned to work with
+*geoclaw-landspill*, so end-users usually don't need to configure them. Here we
+list some settings in these three parts that users must configure:
+
+* `rundata.clawdata.lower[:]`: A list of two `float` representing *xmin* and
+  *ymin*, respectively. In other words, lower boundaries in x and y.
+
+* `rundata.clawdata.upper[:]`: A list of two `float` representing *xmax* and
+  *ymax*, respectively. In other words, upper boundaries in x and y.
+
+* `rundata.clawdata.num_cells[:]`: A list of two `int` representing the numbers
+  of cells at the coarsest AMR level in the x and y direction.
+
+* `rundata.topo_data.topofiles`: A list of two-element lists. For each of the
+  two-element lists, the first element is an `int` indicating the topography
+  file's format. The topography file is specified in the second element. If a
+  relative path is used for the topography file, it must be relative to the
+  case's folder. Usually, the format is `3`, which means it's an
+  [Esri ASCII format raster](https://en.wikipedia.org/wiki/Esri_grid) file.
+  See [Topography data](https://www.clawpack.org/topo.html#topo) for other
+  acceptable raster formats.
+
+* `rundata.clawdata.output_style`: An `int` specifying which output style to
+  use for saving temporal results. Each style comes with different follow-up
+  settings that need to be set:
+
+  * `1`: Output a fixed number of frames at equally spaced times up to a final
+    time.
+    * `rundata.clawdata.num_output_times`: an `int`; total number of frames
+    * `rundata.clawdata.tfinal`: a `float`; time of the last frame
+    * `rundata.clawdata.output_t0`: a `bool`; whether to additionally output the
+      initial state
+  * `2`: Specify a list of times.
+    * `rundata.clawdata.output_times`: a list of `float`; output at these times
+  * `3`: Specify a number of time steps, and output every this number of time
+    steps up to a given number of total steps.
+    * `rundata.clawdata.output_step_interval`: an `int`; output every this
+      number of steps
+    * `rundata.clawdata.total_steps`: an `int`; output until the total number of
+      steps reaches this value
+    * `rundata.clawdata.output_t0`: a `bool` indicating whether to also output
+      the initial state
+
+--------------------------------------
+## II. Settings specific for landspill
+
+### i. Basic
+
+* `rundata.landspill_data.ref_mu`: Reference dynamic viscosity at
+  `ref_temperature`. The default is the viscosity of Maya crude oil. Unit: cP
+  (i.e., mPa-s, or 1e-3 kg/s-m). (default: 332.0)
+
+* `rundata.landspill_data.ref_temperature`: The temperature at which the
+  `ref_mu` is measured. Unit: Celsius. (default: 15.0)
+
+* `rundata.landspill_data.ambient_temperature`: The ambient temperature in
+  simulation. The solver adjusts the viscosity and density based on this
+  temperature. Unit: Celsius. (default: 25.0)
+
+* `rundata.landspill_data.density`: The density measured at `ref_temperature`.
+  The default is the density of Maya crude oil. Unit: kg/m^3. (default: 926.6)
+
+### ii. Point sources
+
+A point source is used to mimic a rupture point along a pipeline. It provides
+fluid inflow into a computational domain. The inflow profile can be a
+multi-stage constant rate.
+
+* `rundata.landspill_data.point_source.n_point_sources`: An `int` indicating how
+  many point sources. ***NOTE: currently only cases with one point source are
+  well tested.*** (default: 0)
+
+* `rundata.landspill_data.point_source.point_sources`: A list of four-element
+  lists. Each four-element list describes a point source. (default: [])
+  * The first element is a length-2 list of x and y coordinates of the point
+    source.
+  * The second element is the number of stages in the multi-stage inflow rate.
+  * The third element is a list indicating the end time of each stage.
+  * The fourth element is a list providing the volumetric rate of each stage.
+
+  For example, if there is one single point source located at x=10, y=11. and if
+  it has three stages of inflow profile (excluding when the rate is zero):
+  * 1 m^3/sec when time &lt; 60 seconds
+  * 0.5 m^3/sec when time &ge; 60 and &lt; 1800 seconds
+  * 0.1 m^3/sec when time &ge; 1800 and &lt; 7200 seconds
+  * 0 m^3/sec after 7200 seconds
+
+  Then the corresponding configuration is
+
+  ```python
+  rundata.landspill_data.point_source.n_point_sources = 1
+  rundata.landspill_data.point_source.point_sources = [
+     [[10., 11.], 3, [60., 1800., 7200.], [1., 0.5, 0.1]]
+  ]
+  ```
+
+### iii. Darcy-Weisbach friction model
+
+In Darcy-Weisbach friction model, the friction coefficient can be calculated
+by several different models. We implemented some of them. The type of models to
+use is controlled by:
+
+* `rundata.landspill_data.darcy_weisbach_friction.type`: (default: 0)
+
+Each type has its own set of follow-up parameters.
+
+* type is `0`: No Darcy-Weisbach friction.
+
+* type is `1`: Use a constant friction coefficient everywhere. The only
+follow-up parameter is
+`rundata.landspill_data.darcy_weisbach_friction.coefficient` (default: 0.25).
+
+* type is `2`: Block-wise constant coefficients. Users specify several blocks in
+  a computational domain, and each block has a constant coefficient.  Available
+  parameters are:
+
+    * `rundata.landspill_data.darcy_weisbach_friction.n_blocks`: Number of blocks.
+      (default: 0) 
+    * `rundata.landspill_data.darcy_weisbach_friction.xlowers`: The lower boundaries
+      in x direction of all blocks. (default: [])
+    * `rundata.landspill_data.darcy_weisbach_friction.ylowers`: The lower boundaries
+      in y direction of all blocks. (default: [])
+    * `rundata.landspill_data.darcy_weisbach_friction.xuppers`: The upper boundaries
+      in x direction of all blocks. (default: [])
+    * `rundata.landspill_data.darcy_weisbach_friction.yuppers`: The upper boundaries
+      in y direction of all blocks. (default: [])
+    * `rundata.landspill_data.darcy_weisbach_friction.coefficients`: The
+      coefficients in all blocks. (default: [])
+    * `rundata.landspill_data.darcy_weisbach_friction.default_coefficient`: Regions
+      that are not covered by blocks will have this value. (default: 0.25)
+
+* type is `3`: Cell-wise coefficients. The coefficients are set through a raster
+  file in Esri ASCII format.
+
+    * `rundata.landspill_data.darcy_weisbach_friction.filename`: the name of the
+      raster file for friction coefficients. (default: "")
+    * `rundata.landspill_data.darcy_weisbach_friction.default_coefficient`: the
+      coefficient for regions not covered by the raster file. (default: 0.25)
+
+* type is `4`: Three-regime coefficient model. The coefficient is determined
+in each cell based on whether the cell is laminar, transient, or turbulent
+flow. See reference [1]. 
+
+    * `rundata.landspill_data.darcy_weisbach_friction.filename`: the name of the
+      raster file for surface roughness. (default: "")
+    * `rundata.landspill_data.darcy_weisbach_friction.default_roughness`: the
+      surface roughness for regions not covered by the raster file. (default:
+      0.0)
+
+* type is `5`: Churchill's coefficient model. See reference [2].
+
+    * `rundata.landspill_data.darcy_weisbach_friction.filename`: the name of the
+      raster file for surface roughness. (default: "")
+    * `rundata.landspill_data.darcy_weisbach_friction.default_roughness`: the
+      surface roughness for regions not covered by the raster file. (default:
+      0.0)
+
+* type is `6`: Two-regime coefficient. Similar to the three-regime model but
+ignore transient flow regime.
+
+    * `rundata.landspill_data.darcy_weisbach_friction.filename`: the name of the
+      raster file for surface roughness. (default: "")
+    * `rundata.landspill_data.darcy_weisbach_friction.default_roughness`: the
+      surface roughness for regions not covered by the raster file. (default:
+      0.0)
+
+Usually, type `4` or `5` is used because they consider the surface roughness.
+
+### iv. In-land waterbodies
+
+*geoclaw-landspill* is able to records how much fluid volume flowing into
+in-land waterbodies and at what locations. Users have to provide rasterized
+hydrology data to enable this feature.
+
+* `rundata.landspill_data.hydro_features.files`: A list of filenames. If
+  relative paths are provided, they are assumed to be relative to the case
+  folder. (default: [])
+
+***Note: using multiple files has not been widely tested. It's recommended to
+combine all data into one single raster file.***
+
+Usually, hydrology data are vector data, such as those obtained from the USGS
+database. Users can burn the vector data into raster files using
+[`gdal_rasterize`](https://gdal.org/programs/gdal_rasterize.html).
+Alternatively, users can specify a filename without actually providing the file.
+*geoclaw-landspill* automatically downloads and rasterize hydrology data in this
+case.
+
+### v. Evaporation model
+
+*geoclaw-landspill* implements Fingas' 1996 model, including the natural log
+model and square root model. See reference [3].
+
+* `rundata.landspill_data.evaporation.type`: The type of evaporation models to
+  use. (default: 0)
+
+  * `0`: No evaporation.
+  * `1`: Fingas' natural log model
+  * `2`: Fingas' square root model
+
+* `rundata.landspill_data.evaporation.coefficients`: a list of `float` for model
+  coefficients. Currently, we only have Fingas' models, so this variable should
+  be a list of two `float`. See reference [3] for coefficients of a variety of
+  oils.
+
+-------------------------------------------------------------------
+## III. Optional settings affecting stability, accuracy, and performance
+
+The following optional settings are commonly used to control numerical
+stability, accuracy, or/and performance:
+
+* `rundata.clawdata.dt_initial`: This determines the size of the first time
+  step. By default, it is automatically determined by the following
+  formula:
+
+  ```
+  dt_initial = 0.3 * (cell size at the finest AMR grid) / (inflow rate)
+  ```
+
+  This makes the cell containing the first point source to have a depth of 0.3
+  meters at the end of the first time step. Currently, the coefficient, 0.3, is
+  hard-coded. This may give some trouble to simulations under small scales.
+  For example, if the inflow rate of a point source is 1e-6 m^3 / sec, and if
+  the finest cell size is 1e-4 m^2, this formula returns a `dt_initial` of 30
+  seconds, which is apparently too big for a problem at this scale. In this
+  case, users may want to manually set `dt_initial`.
+
+* `rundata.clawdata.dt_max`: The maximum time step size allowed during a
+  simulation. By default, the solver adjusts time step sizes based on flow
+  conditions. Users sometimes may want to limit or increase how big a step size
+  can be. (default: 5.0)
+
+* `rundata.clawdata.cfl_desired`: The desired Courant–Friedrichs–Lewy (CFL)
+  number. The solver adjusts time step sizes to make the resulting CFL numbers
+  below this value. (default: 0.9)
+
+* `rundata.clawdata.cfl_max`: The maximum allowed CFL number. If a CFL number is
+  larger than this value, it triggers the solver to adjust time step sizes.
+  (default: 0.95)
+
+* `rundata.refinement_data.variable_dt_refinement_ratios`: Whether to use
+  GeoClaw's adaptive time-step refinement mechanism. When AMR refines a coarse
+  cell into several smaller cells, the smaller cells need smaller time step
+  sizes. The time step sizes' refinement ratios can be fixed values or
+  automatically determined by the GeoClaw solver. (default: True)
+
+* `rundata.amrdata.amr_levels_max`: The number of AMR levels to be used. By
+  default, *geoclaw-landspill* uses two levels: the coarse level for dry
+  regions and the fine level for wet regions. (default: 2)
+
+* `rundata.amrdata.refinement_ratios_x`: The refinement ratio between each two
+  consecutive AMR levels in the x direction. (default: [4])
+
+* `rundata.amrdata.refinement_ratios_y`: The refinement ratio between each two
+  consecutive AMR levels in the y direction. (default: [4])
+
+* `rundata.amrdata.refinement_ratios_t`: The refinement ratio between each two
+  consecutive AMR levels in time step sizes. This only have effect when
+  `rundata.refinement_data.variable_dt_refinement_ratios` is `False`. (default:
+  [4])
+
+* `rundata.geo_data.dry_tolerance`: The solver considers cells with depth below
+  this value to be dry cells. (default: 1.e-4)
+
+* `rundata.landspill_data.update_tol`: Affect how a coarse cell should
+  distribute the depth to its children cells. This setting affects the
+  conservation of mass. (default: the same as the `dry_tolerance`)
+
+* `rundata.landspill_data.refine_tol`: By default, as long as there's fluid in a
+  cell, it is refined by AMR (even if the depth is smaller than `dry_tolerance`.
+  The behavior can be changed with this variable. (default: 0.0)
+
+------------
+## Reference
+
+[1] B. C. Yen, “Open Channel Flow Resistance,” Journal of Hydraulic Engineering,
+vol. 128, no. 1, pp. 20–39, Jan. 2002, doi: 10.1061/(ASCE)0733-9429(2002)128:1(20).
+
+[2] S. W. Churchill, “Friction-factor equation spans all fluid flow regimes,”
+Chem. Eng., vol. 84, no. 24, pp. 91–92, 1977
+
+[3] M. F. Fingas, “Modeling evaporation using models that are not boundary-layer
+regulated,” Journal of Hazardous Materials, vol. 107, no. 1–2, pp. 27–36, Feb.
+2004, doi: 10.1016/j.jhazmat.2003.11.007.
diff --git a/gclandspill/__init__.py b/gclandspill/__init__.py
index 96a021e..e16c7c0 100644
--- a/gclandspill/__init__.py
+++ b/gclandspill/__init__.py
@@ -18,5 +18,5 @@
 from . import data
 
 # meta
-__version__ = "1.0.dev2"
+__version__ = "1.0.dev3"
 __author__ = "Pi-Yueh Chuang <pychuang@gwu.edu>"
diff --git a/gclandspill/_misc.py b/gclandspill/_misc.py
index b797276..a169d8d 100644
--- a/gclandspill/_misc.py
+++ b/gclandspill/_misc.py
@@ -40,6 +40,11 @@ class NoFrameDataError(Exception):
     pass  # pylint: disable=unnecessary-pass
 
 
+class FrictionModelError(Exception):
+    """An error to raise when both Manning and Darcy-Weisbach models are enabled."""
+    pass  # pylint: disable=unnecessary-pass
+
+
 class DatetimeCtrlParams(TypedDict):
     """A type definition for the `dict` controlling the timestamps."""
     apply_datetime_stamp = bool
diff --git a/gclandspill/data.py b/gclandspill/data.py
index 2be069b..42fbad1 100644
--- a/gclandspill/data.py
+++ b/gclandspill/data.py
@@ -7,17 +7,149 @@
 # Distributed under terms of the BSD 3-Clause license.
 
 
-"""Classes representing parameters for geoclaw-landspill runs
+"""Classes representing parameters for geoclaw-landspill runs.
 """
 import os
-import gclandspill.clawutil.data
-# pylint: disable=no-member, attribute-defined-outside-init, fixme
-# TODO: merge the friction models of GeoClaw and landspill
-# TODO: merge dry_tol and dry_tolerance from GeoClaw and landspill
+from functools import reduce as _reduce
+from gclandspill import clawutil as _clawutil
+from gclandspill import _misc
+
+
+class ClawRunData(_clawutil.data.ClawRunData):
+    """A modified version of ClawRunData for geoclaw-landspill.
+
+    This ClawRunData has an additional LandSpillData instance. Also, the following default values are modified:
+
+        # core classical clawpack configuration
+        self.clawdata.num_eqn = 2
+        self.clawdata.num_waves = 3
+        self.clawdata.num_aux = 2
+        self.clawdata.output_format = 3  # binary
+        self.clawdata.output_aux_components = "all"
+        self.clawdata.output_aux_onlyonce = False
+        self.clawdata.dt_initial = None
+        self.clawdata.dt_max = 5.0
+        self.clawdata.cfl_max = 0.95
+        self.clawdata.steps_max = 100000
+        self.clawdata.verbosity = 5
+        self.clawdata.verbosity_regrid = 5
+        self.clawdata.source_split = 1  # godunov
+        self.clawdata.limiter = [4, 4, 4]  # use "mc" limiter for all equations
+        self.clawdata.use_fwaves = True
+        self.clawdata.bc_lower = [1, 1]
+        self.clawdata.bc_upper = [1, 1]
+
+        # AMRClaw configuration
+        self.amrdata.amr_levels_max = 2
+        self.amrdata.refinement_ratios_x = [4]
+        self.amrdata.refinement_ratios_y = [4]
+        self.amrdata.refinement_ratios_t = [4]
+        self.amrdata.aux_type = ["center", "center"]  # aux variables are defined at cell centers
+        self.amrdata.regrid_interval = 1
+        self.amrdata.verbosity_regrid = 5
+
+        # GeoClaw basic configuration
+        self.geo_data.gravity = 9.81
+        self.geo_data.coriolis_forcing = False
+        self.geo_data.sea_level = -1000.
+        self.geo_data.dry_tolerance = 1e-4
+        self.geo_data.friction_forcing = False  # turned off in favor of Darcy-Weisbach friction
+
+        # GeoClaw refinement mechanism
+        self.refinement_data.wave_tolerance = 1.0e-5
+        self.refinement_data.speed_tolerance = [1e-5] * 6
+        self.refinement_data.variable_dt_refinement_ratios = True
+    """
+
+    def __init__(self):
+        # pylint: disable=no-member
+        super().__init__("geoclaw", 2)  # to first get attributes from geoclaw and with dimension 2
+
+        # modify attributes for geoclaw-landspill
+        self.pkg = "geoclaw-landspill"
+        self.xclawcmd = "geoclaw-landspill-bin"
+
+        # append a landspill data instance to the data object
+        self.add_data(LandSpillData(), 'landspill_data')
+
+        # core classical clawpack configuration
+        self.clawdata.num_eqn = 3
+        self.clawdata.num_waves = 3
+        self.clawdata.num_aux = 2
+        self.clawdata.output_format = 3  # binary
+        self.clawdata.output_aux_components = "all"
+        self.clawdata.output_aux_onlyonce = False
+        self.clawdata.dt_initial = None
+        self.clawdata.dt_max = 5.0
+        self.clawdata.cfl_max = 0.95
+        self.clawdata.steps_max = 100000
+        self.clawdata.verbosity = 5
+        self.clawdata.verbosity_regrid = 5
+        self.clawdata.source_split = 1  # godunov
+        self.clawdata.limiter = [4, 4, 4]  # use "mc" limiter for all equations
+        self.clawdata.use_fwaves = True
+        self.clawdata.bc_lower = [1, 1]
+        self.clawdata.bc_upper = [1, 1]
+
+        # AMRClaw configuration
+        self.amrdata.amr_levels_max = 2
+        self.amrdata.refinement_ratios_x = [4]
+        self.amrdata.refinement_ratios_y = [4]
+        self.amrdata.refinement_ratios_t = [4]
+        self.amrdata.aux_type = ["center", "center"]  # aux variables are defined at cell centers
+        self.amrdata.regrid_interval = 1
+        self.amrdata.verbosity_regrid = 5
+
+        # GeoClaw basic configuration
+        self.geo_data.gravity = 9.81
+        self.geo_data.coriolis_forcing = False
+        self.geo_data.sea_level = -1000.
+        self.geo_data.dry_tolerance = 1e-4
+        self.geo_data.friction_forcing = False  # turned off in favor of Darcy-Weisbach friction
+
+        # GeoClaw refinement mechanism
+        self.refinement_data.wave_tolerance = 1.0e-5
+        self.refinement_data.speed_tolerance = [1e-5] * 6
+        self.refinement_data.variable_dt_refinement_ratios = True
+
+    def write(self, out_dir=""):
+
+        # check and set some automatically determined values
+        self._check()
+
+        super().write(out_dir)
+
+    def _check(self):
+        """Check and set some automatically determined values."""
+        # pylint: disable=no-member
+
+        # automatically determine dt_initial
+        if self.clawdata.dt_initial is None:
+            cell_area = \
+                (self.clawdata.upper[0] - self.clawdata.lower[0]) * \
+                (self.clawdata.upper[1] - self.clawdata.lower[1]) / \
+                (self.clawdata.num_cells[0] * self.clawdata.num_cells[1])
+            cell_area /= _reduce(lambda x, y: x*y, self.amrdata.refinement_ratios_x[:self.amrdata.amr_levels_max], 1)
+            cell_area /= _reduce(lambda x, y: x*y, self.amrdata.refinement_ratios_y[:self.amrdata.amr_levels_max], 1)
+            vrate = self.landspill_data.point_sources.point_sources[0][-1][0]
+
+            self.clawdata.dt_initial = 0.3 * cell_area / vrate
+
+        if self.landspill_data.update_tol is None:
+            self.landspill_data.update_tol = self.geo_data.dry_tolerance
+
+        if self.landspill_data.darcy_weisbach_friction.type != 0 and self.geo_data.friction_forcing:
+            raise _misc.FrictionModelError("Both Manning and Darcy-Weisbach models are enabled.")
+
+        if self.landspill_data.darcy_weisbach_friction.dry_tol is None:
+            self.landspill_data.darcy_weisbach_friction.dry_tol = self.geo_data.dry_tolerance
+
+        if self.landspill_data.darcy_weisbach_friction.friction_tol is None:
+            self.landspill_data.darcy_weisbach_friction.friction_tol = self.geo_data.friction_depth
 
 
 # Land-spill data
-class LandSpillData(gclandspill.clawutil.data.ClawData):
+class LandSpillData(_clawutil.data.ClawData):
     """Data object describing land spill simulation configurations"""
 
     def __init__(self):
@@ -38,7 +170,7 @@ def __init__(self):
         self.add_attribute('density', 926.6)  # unit: kg / m^3
 
         # extra parameters to hack the AMR algorithms for landspill applications
-        self.add_attribute("update_tol", 0.0)
+        self.add_attribute("update_tol", None)
         self.add_attribute("refine_tol", 0.0)
 
         # add point source data
@@ -83,29 +215,29 @@ def write(self, out_file='landspill.data', data_source="setrun.py"):
         # output point sources data
         self.data_write('point_sources_file',
                         description='File name of point sources settings')
-        self.point_sources.write(os.path.join(base, self.point_sources_file))
+        self.point_sources.write(os.path.join(base, self.point_sources_file))  # pylint: disable=no-member
 
         # output Darcy-Weisbach data
         self.data_write('darcy_weisbach_file',
                         description='File name of Darcy-Weisbach settings')
-        self.darcy_weisbach_friction.write(os.path.join(base, self.darcy_weisbach_file))
+        self.darcy_weisbach_friction.write(os.path.join(base, self.darcy_weisbach_file))  # pylint: disable=no-member
 
         # output hydroological feature data
         self.data_write('hydro_feature_file',
                         description='File name of hydrological feature settings')
-        self.hydro_features.write(os.path.join(base, self.hydro_feature_file))
+        self.hydro_features.write(os.path.join(base, self.hydro_feature_file))  # pylint: disable=no-member
 
         # output evaporation data
         self.data_write('evaporation_file',
                         description='File name of evaporation settings')
-        self.evaporation.write(os.path.join(base, self.evaporation_file))
+        self.evaporation.write(os.path.join(base, self.evaporation_file))  # pylint: disable=no-member
 
         # close the output file
         self.close_data_file()
 
 
 # Point source data
-class PointSourceData(gclandspill.clawutil.data.ClawData):
+class PointSourceData(_clawutil.data.ClawData):
     """Data object describing point sources"""
 
     def __init__(self):
@@ -132,7 +264,7 @@ def write(self, out_file='point_source.data', data_source="setrun.py"):
         self.data_write('n_point_sources', description='Number of point sources')
 
         # write point sources
-        for i, pts in enumerate(self.point_sources):
+        for i, pts in enumerate(self.point_sources):  # pylint: disable=no-member
             self.data_write()  # a blank line
             self.data_write("id", i, description="ID of this point source")
             self.data_write("coord", pts[0], description="coordinates")
@@ -145,6 +277,7 @@ def write(self, out_file='point_source.data', data_source="setrun.py"):
 
     def _check(self):
         """Check if the data are consistent"""
+        # pylint: disable=no-member
 
         # check if the number of data set match n_point_sources
         assert self.n_point_sources == len(self.point_sources),  \
@@ -174,7 +307,7 @@ def _check(self):
 
 
 # Darcy-Weisbach data
-class DarcyWeisbachData(gclandspill.clawutil.data.ClawData):
+class DarcyWeisbachData(_clawutil.data.ClawData):
     """Data object describing Darcy-Weisbach friction model"""
 
     def __init__(self):
@@ -192,11 +325,11 @@ def __init__(self):
         # 6: two-regime coefficient model: laminar + Colebrook-White
         self.add_attribute('type', 0)
 
-        # friction_tol, the same as friction_tol in GeoClaw.
-        self.add_attribute('friction_tol', 1e6)
+        # friction_tol, the same as the friction_depth in GeoClaw.
+        self.add_attribute('friction_tol', None)
 
         # dry_tol, the same as dry_tolerance in GeoClaw.
-        self.add_attribute('dry_tol', 1e-6)
+        self.add_attribute('dry_tol', None)
 
         # single constant coefficient.
         self.add_attribute('coefficient', 0.25)
@@ -224,6 +357,7 @@ def __init__(self):
 
     def write(self, out_file='darcy_weisbach.data', data_source="setrun.py"):
         """Write out the data file to the path given"""
+        # pylint: disable=access-member-before-definition, attribute-defined-outside-init, no-member
 
         # check data consistency
         self._check()
@@ -285,6 +419,7 @@ def write(self, out_file='darcy_weisbach.data', data_source="setrun.py"):
 
     def _check(self):
         """Check if the data are consistent"""
+        # pylint: disable=no-member
 
         assert isinstance(self.type, int), "Type should be an integer."
 
@@ -342,7 +477,7 @@ def _check(self):
 
 
 # Hydrologic feature data
-class HydroFeatureData(gclandspill.clawutil.data.ClawData):
+class HydroFeatureData(_clawutil.data.ClawData):
     """Data object describing hydrologic features"""
 
     def __init__(self):
@@ -360,11 +495,11 @@ def write(self, out_file='hydro_feature.data', data_source="setrun.py"):
         self.open_data_file(out_file, data_source)
 
         # write number of files
-        self.data_write('n_files', len(self.files),
+        self.data_write('n_files', len(self.files),  # pylint: disable=no-member
                         description='Number of hydro files')
 
         # write file names line by line
-        for i, single_file in enumerate(self.files):
+        for i, single_file in enumerate(self.files):  # pylint: disable=no-member
             # relative path will be relative to the folder containing `out_file` (following how
             # geoclaw handles topo files)
             if not os.path.isabs(single_file):
@@ -376,7 +511,7 @@ def write(self, out_file='hydro_feature.data', data_source="setrun.py"):
 
 
 # Evaporation data
-class EvaporationData(gclandspill.clawutil.data.ClawData):
+class EvaporationData(_clawutil.data.ClawData):
     """Data object describing evaporation"""
 
     def __init__(self):
@@ -406,10 +541,10 @@ def write(self, out_file='evaporation.data', data_source="setrun.py"):
         self.data_write('type', description='Evaporation type')
 
         # number of coefficients
-        self.data_write('n_coefficients', len(self.coefficients),
+        self.data_write('n_coefficients', len(self.coefficients),  # pylint: disable=no-member
                         description='Number of evaporation coefficients.')
 
-        for i, coeff in enumerate(self.coefficients):
+        for i, coeff in enumerate(self.coefficients):  # pylint: disable=no-member
             self.data_write('C{}'.format(i), coeff,
                             description='Coefficient {}'.format(i))
 
@@ -419,5 +554,5 @@ def write(self, out_file='evaporation.data', data_source="setrun.py"):
     def _check(self):
         """Check if the data are consistent"""
 
-        assert isinstance(self.type, int), "Evaporation type should be an integer."
-        assert self.type in [0, 1, 2], "Evaporation type should be in [0, 1, 2]"
+        assert isinstance(self.type, int), "Evaporation type should be an integer."  # pylint: disable=no-member
+        assert self.type in [0, 1, 2], "Evaporation type should be in [0, 1, 2]"  # pylint: disable=no-member
diff --git a/paper/paper.bib b/paper/paper.bib
new file mode 100644
index 0000000..5926b03
--- /dev/null
+++ b/paper/paper.bib
@@ -0,0 +1,202 @@
+
+@article{Berger2011,
+	title = {The {GeoClaw} software for depth-averaged flows with adaptive refinement},
+	volume = {34},
+	issn = {03091708},
+	url = {https://linkinghub.elsevier.com/retrieve/pii/S0309170811000480},
+	doi = {10.1016/j.advwatres.2011.02.016},
+	pages = {1195--1206},
+	number = {9},
+	journaltitle = {Advances in Water Resources},
+	author = {Berger, Marsha J. and George, David L. and {LeVeque}, Randall J. and Mandli, Kyle T.},
+	date = {2011-09},
+	keywords = {{GeoClaw}, proposal},
+}
+
+@article{mehrotra_generalized_1991,
+	title = {Generalized one-parameter viscosity equation for light and medium liquid hydrocarbons},
+	volume = {30},
+	pages = {1367--1372},
+	number = {6},
+	journaltitle = {Industrial \& engineering chemistry research},
+	author = {Mehrotra, Anil K.},
+	date = {1991},
+	langid = {english},
+}
+
+@article{Yen2002,
+	title = {Open Channel Flow Resistance},
+	volume = {128},
+	issn = {0733-9429},
+	url = {http://ascelibrary.org/doi/10.1061/%28ASCE%290733-9429%282002%29128%3A1%2820%29},
+	doi = {10.1061/(ASCE)0733-9429(2002)128:1(20)},
+	pages = {20--39},
+	number = {1},
+	journaltitle = {Journal of Hydraulic Engineering},
+	author = {Yen, Ben Chie},
+	date = {2002-01},
+	note = {{ISBN}: 0733-9429},
+}
+
+@article{churchill_friction-factor_1977,
+	title = {Friction-factor equation spans all fluid flow regimes},
+	volume = {84},
+	pages = {91--92},
+	number = {24},
+	journaltitle = {Chemical Engineering},
+	shortjournal = {Chem. Eng.},
+	author = {Churchill, Stuart W.},
+	date = {1977},
+	note = {{ISBN}: 00092460},
+	keywords = {oil overland flow},
+}
+
+@article{fingas_modeling_2004,
+	title = {Modeling evaporation using models that are not boundary-layer regulated},
+	volume = {107},
+	issn = {03043894},
+	url = {http://www.ncbi.nlm.nih.gov/pubmed/15036640},
+	doi = {10.1016/j.jhazmat.2003.11.007},
+	pages = {27--36},
+	number = {1},
+	journaltitle = {Journal of Hazardous Materials},
+	author = {Fingas, Mervin F.},
+	date = {2004-02-27},
+	pmid = {15036640},
+	note = {{ISBN}: 0304-3894},
+}
+
+@article{belvederesi_statistical_2018,
+	title = {Statistical analysis of environmental consequences of hazardous liquid pipeline accidents},
+	volume = {4},
+	issn = {24058440},
+	url = {https://linkinghub.elsevier.com/retrieve/pii/S2405844018356020},
+	doi = {10.1016/j.heliyon.2018.e00901},
+	pages = {e00901},
+	number = {11},
+	journaltitle = {Heliyon},
+	shortjournal = {Heliyon},
+	author = {Belvederesi, Chiara and Thompson, Megan S. and Komers, Petr E.},
+	urldate = {2021-01-15},
+	date = {2018-11},
+	langid = {english},
+	keywords = {oil overland flow},
+}
+
+@inproceedings{Paige2002,
+	title = {Terrestrial Spill Modeling -- Increasing Confidence in the Estimation of High Consequence Area ({HCA}) Impact},
+	pages = {23--25},
+	booktitle = {{GeoFields} Inc, {GIS} for Oil and Gas Conference},
+	author = {Paige, Don and Park, Nick and Posner, Jonathan},
+	date = {2002},
+}
+
+@inproceedings{Ronnie2004,
+	title = {Prediction of the Consequences of Accidental Releases From Liquid Pipelines},
+	isbn = {0-7918-4176-6},
+	url = {http://proceedings.asmedigitalcollection.asme.org/proceeding.aspx?articleid=1645597},
+	doi = {10.1115/IPC2004-0652},
+	pages = {581--586},
+	booktitle = {2004 International Pipeline Conference, Volumes 1, 2, and 3},
+	publisher = {{ASME}},
+	author = {Ronnie, S E},
+	date = {2004},
+}
+
+@article{Simmons2003,
+	title = {Status of Models for Land Surface Spills of Nonaqueous Liquids},
+	pages = {62},
+	issue = {August},
+	author = {Simmons, C. S. and Keller, J. M.},
+	date = {2003},
+}
+
+@article{Guo2006,
+	title = {Overland flow model for asphalt oil spills},
+	volume = {78},
+	issn = {03014797},
+	url = {http://linkinghub.elsevier.com/retrieve/pii/S0301479705002264},
+	doi = {10.1016/j.jenvman.2005.05.018},
+	pages = {102--105},
+	number = {1},
+	journaltitle = {Journal of Environmental Management},
+	author = {Guo, J.C.Y.},
+	date = {2006-01},
+}
+
+@inproceedings{farrar_gis_2005,
+	location = {San Diego, {CA}, {USA}},
+	title = {{GIS} technology applied to modeling oil spills on land},
+	volume = {20},
+	url = {http://proceedings.esri.com/library/userconf/proc05},
+	eventtitle = {25th Annual Esri International User Conference},
+	booktitle = {25th Annual Esri International User Conference},
+	author = {Farrar, Will and Galagan, Chris and Isaji, Tatsu and Knee, Kelly},
+	date = {2005},
+}
+
+@inproceedings{Zuczek2008,
+	title = {An Overland-Hydrographical Spill Model and Its Application to Pipeline Consequence Modeling},
+	isbn = {978-0-7918-4860-9},
+	url = {http://proceedings.asmedigitalcollection.asme.org/proceeding.aspx?articleid=1641159},
+	doi = {10.1115/IPC2008-64389},
+	pages = {429--437},
+	booktitle = {2008 7th International Pipeline Conference, Volume 4},
+	publisher = {{ASME}},
+	author = {Zuczek, Patricia and Deng, Chuntao and Mihell, Jim and Adams, Keith},
+	date = {2008},
+}
+
+@inproceedings{Su2017,
+	title = {Simulation Model for Fuel Spreading on Complex Topography Based on Hydrodynamics},
+	isbn = {978-981-322-035-5},
+	url = {https://www.worldscientific.com/doi/abs/10.1142/9789813220362_0021},
+	doi = {10.1142/9789813220362_0021},
+	pages = {151--157},
+	booktitle = {Electrical Engineering and Automation},
+	publisher = {{WORLD} {SCIENTIFIC}},
+	author = {Su, Bo-ni and Huang, Hong and Zhang, Nan},
+	date = {2017-06},
+}
+
+@online{RPSGroup,
+	title = {{OILMAPLAND}: Land Spill Model and Analysis System},
+	url = {http://asascience.com/software/oilmap/oilmapland.shtml},
+	author = {{RPS Group}},
+	urldate = {2019-08-26}
+}
+
+@online{Hydronia,
+	title = {{OilFlow}2D},
+	url = {http://www.hydronia.com/oilflow2d/},
+	author = {{Hydronia}},
+	urldate = {2019-08-26}
+}
+
+@article{Hussein2002,
+	title = {Development and verification of a screening model for surface spreading of petroleum},
+	volume = {57},
+	issn = {01697722},
+	url = {http://linkinghub.elsevier.com/retrieve/pii/S0169772201002200},
+	doi = {10.1016/S0169-7722(01)00220-0},
+	pages = {281--302},
+	number = {3},
+	journaltitle = {Journal of Contaminant Hydrology},
+	author = {Hussein, Maged and Jin, Minghui and Weaver, James W.},
+	date = {2002-08},
+	pmid = {12180813},
+	note = {{ISBN}: 0169-7722},
+}
+
+@inproceedings{Gin2012,
+	location = {Calgary, Alberta, Canada},
+	title = {Spill Impact Assessment: Crude Oil Pipeline {ERCB} License Number 3353-1},
+	isbn = {978-0-7918-4512-7},
+	url = {http://proceedings.asmedigitalcollection.asme.org/proceeding.aspx?doi=10.1115/IPC2012-90252},
+	doi = {10.1115/IPC2012-90252},
+	pages = {283},
+	booktitle = {Proceedings of the 2012 9th International Pipeline Conference, {IPC}2012},
+	publisher = {{ASME}},
+	author = {Gin, Gordon E. and Davis, John},
+	date = {2012-09-24},
+}
diff --git a/paper/paper.md b/paper/paper.md
new file mode 100644
index 0000000..7492530
--- /dev/null
+++ b/paper/paper.md
@@ -0,0 +1,100 @@
+---
+title: "geoclaw-landspill: an oil land-spill and overland flow simulator for pipeline rupture events"
+tags:
+  - shallow-water-equations
+  - overland-flow
+  - land-spill
+  - pipeline
+  - geoclaw
+authors:
+  - name: Pi-Yueh Chuang
+    orcid: 0000-0001-6330-2709
+    affiliation: 1
+  - name: Tracy Thorleifson
+    affiliation: 2
+  - name: Lorena A. Barba
+    orcid: 0000-0001-5812-2711
+    affiliation: 1
+affiliations:
+  - name: Department of Mechanical and Aerospace Engineering, the George Washington University, Washington, DC, USA
+    index: 1
+  - name: G2 Integrated Solutions, Houston, TX, USA
+    index: 2
+date: 15 January 2021
+bibliography: paper.bib
+---
+
+# Summary
+
+*geoclaw-landspill* is a package to numerically simulate oil land-spill and overland flow that occurs during pipeline accidents.
+It helps us understand how oil flows above ground after accidental release and to study a pipeline's impact on the environment and economy.
+By understanding the impact of a to-be-constructed pipeline, we can choose a pipeline route with the least loss if accidents happen.
+On the other hand, understanding how oil flows also helps develop rescue teams' deployment strategies and remedy plans for potential accidents.
+
+The package provides a numerical solver that solves full shallow-water equations and post-processing utilities.
+The numerical solver is an expanded version of GeoClaw [@Berger2011].
+GeoClaw is a parallel shallow-water equation solver for tsunami simulations using adaptive mesh refinement (AMR) and finite-volume methods.
+We added several new features and modifications to simulate the overland flow of pipeline rupture events, including (citations refer to the details of the adopted models):
+
+* point sources with multi-stage inflow rates to mimic rupture points along a pipeline,
+* Lewis Squires Correlation for temperature-dependent flow viscosity [@mehrotra_generalized_1991],
+* Darcy-Weisbach friction model with multi-regime coefficient models (laminar, transient, and turbulent regimes) [@Yen2002] and Churchill's model [@churchill_friction-factor_1977],
+* in-land waterbody interactions,
+* Fingas' evaporation models [@fingas_modeling_2004], and
+* optimizations to improve performance in overland flow simulations.
+
+In addition to the numerical solver, *geoclaw-landspill* is also able to:
+
+* automatically download high-resolution topography and hydrology data, and
+* create CF-compliant NetCDF raster files for mainstream GIS software (e.g., QGIS, ArcGIS).
+
+*geoclaw-landspill* includes several examples to showcase its capability to simulate the overland flow on different terrain.
+
+We implemented the core numerical solver and new features using Fortran 2008 to better integrate into the original GeoClaw code.
+Other utilities are in Python.
+Users can find *geoclaw-landspill* as a Python package on PyPI and install it through pip.
+The only dependency that pip does not manage is the Fortran compiler.
+Docker images and Singularity images are also available.
+They ease the deployment of the solver and simulations to cloud-based high-performance clusters.
+
+# Statement of need
+
+In the US, between 2010 and 2017, an average of 388 hazardous liquid pipeline accidents happened per year.
+50% of accidents contaminate soil, and 41% of accidents affect areas with high consequences in either ecology or economy.
+Moreover, 85% on average of the released oil was not recovered and kept damaging the environment.
+[@belvederesi_statistical_2018]
+From the perspective of risk management, while pipelines are unavoidable in modern days, it is necessary to understand how a pipeline may impact the environment if any accidental release happens.
+*geoclaw-landspill* serves this purpose.
+*geoclaw-landspill* provides a free and open-source simulation tool to researchers investigating the danger, risk, and loss posted by potential pipeline accidents.
+
+To our knowledge, *geoclaw-landspill* is the only open-source high-fidelity flow simulator for oil pipeline rupture events.
+High fidelity means the results provide more details and accuracy because of high-resolution digital elevation data, fine spatial discretization, and full shallow-water equations.
+Commercial products with a similar capability to *geoclaw-landspill* are available [@Zuczek2008; @RPSGroup; @Hydronia; @Gin2012].
+Other non-commercial software more or less serving a similar purpose usually relies on simplified models, such as 1D open-channel models, diffusive wave approximation, gravity current models, and gradient-based route selection models [@Hussein2002; @Simmons2003; @Ronnie2004; @farrar_gis_2005; @Guo2006; @Su2017].
+Moreover, these non-commercial codes are either non-exist in modern days or are not open-source.
+
+Another value of *geoclaw-landspill* is to provide a platform for scholars who study oil flow modeling to implement and test their models.
+As the main flow solver is under the BSD 3-Clause License, scholars can add their models to *geoclaw-landspill* freely.
+
+# Past or ongoing research projects using the software
+
+The following conference presentations and posters used previous versions of
+*geoclaw-landspill*:
+
+1. Chuang, P.-Y., Thorleifson, T., & Barba, L. A. (2019a, May). GeoClaw-ArcGIS Integration for Advanced Modeling of Overland Hydrocarbon Plumes. *2019 Petroleum GIS Conference Proceedings. 2019 Esri Petroleum GIS Conference, Houston, TX, USA*.
+
+2. Chuang, P.-Y., Thorleifson, T., & Barba, L. A. (2019b, July). Python Workflow for High-Fidelity Modeling of Overland Hydrocarbon Flows with GeoClaw and Cloud Computing. *Proceedings of the 18th Python in Science Conference. 18th Python in Science Conference (SciPy 2019), Austin, TX, USA*.
+
+
+# Acknowledgments
+
+The project received funding from G2 Integrated Solutions, Houston, TX, USA.
+
+# Author Contributions
+
+Pi-Yueh Chuang developed and is the maintainer of this software.
+Tracy Thorleifson was involved in software design, model selection, and software testing.
+Lorena A. Barba designed the software specification, development roadmap, and software framework.
+She also oversaw the progress of the project.
+
+# Reference


Test patch:
diff --git a/tests/data/regression-1/fort.a0000 b/tests/data/regression-1/fort.a0000
new file mode 100644
index 0000000..488cdb8
Binary files /dev/null and b/tests/data/regression-1/fort.a0000 differ
diff --git a/tests/data/regression-1/fort.a0001 b/tests/data/regression-1/fort.a0001
new file mode 100644
index 0000000..e4cc42a
Binary files /dev/null and b/tests/data/regression-1/fort.a0001 differ
diff --git a/tests/data/regression-1/fort.a0002 b/tests/data/regression-1/fort.a0002
new file mode 100644
index 0000000..79d4e3a
Binary files /dev/null and b/tests/data/regression-1/fort.a0002 differ
diff --git a/tests/data/regression-1/fort.a0003 b/tests/data/regression-1/fort.a0003
new file mode 100644
index 0000000..9742abd
Binary files /dev/null and b/tests/data/regression-1/fort.a0003 differ
diff --git a/tests/data/regression-1/fort.a0004 b/tests/data/regression-1/fort.a0004
new file mode 100644
index 0000000..7c6d26b
Binary files /dev/null and b/tests/data/regression-1/fort.a0004 differ
diff --git a/tests/data/regression-1/fort.a0005 b/tests/data/regression-1/fort.a0005
new file mode 100644
index 0000000..5d7285a
Binary files /dev/null and b/tests/data/regression-1/fort.a0005 differ
diff --git a/tests/data/regression-1/fort.b0001 b/tests/data/regression-1/fort.b0001
index 6cb084d..19d15db 100644
Binary files a/tests/data/regression-1/fort.b0001 and b/tests/data/regression-1/fort.b0001 differ
diff --git a/tests/data/regression-1/fort.b0002 b/tests/data/regression-1/fort.b0002
index 6e13cba..ef2f823 100644
Binary files a/tests/data/regression-1/fort.b0002 and b/tests/data/regression-1/fort.b0002 differ
diff --git a/tests/data/regression-1/fort.b0003 b/tests/data/regression-1/fort.b0003
index acc58e3..15afa10 100644
Binary files a/tests/data/regression-1/fort.b0003 and b/tests/data/regression-1/fort.b0003 differ
diff --git a/tests/data/regression-1/fort.b0004 b/tests/data/regression-1/fort.b0004
index 6fe92f0..c6fd7a1 100644
Binary files a/tests/data/regression-1/fort.b0004 and b/tests/data/regression-1/fort.b0004 differ
diff --git a/tests/data/regression-1/fort.b0005 b/tests/data/regression-1/fort.b0005
index 749f825..de0bff7 100644
Binary files a/tests/data/regression-1/fort.b0005 and b/tests/data/regression-1/fort.b0005 differ
diff --git a/tests/data/regression-1/fort.q0001 b/tests/data/regression-1/fort.q0001
index b4e80f1..b80b758 100644
--- a/tests/data/regression-1/fort.q0001
+++ b/tests/data/regression-1/fort.q0001
@@ -34,7 +34,7 @@
     0.1000000000000000E-02    dx
     0.1000000000000000E-02    dy
 
-     2                 grid_number
+     3                 grid_number
      2                 AMR_level
     32                 mx
     32                 my
diff --git a/tests/data/regression-1/fort.q0002 b/tests/data/regression-1/fort.q0002
index c032473..614cbf9 100644
--- a/tests/data/regression-1/fort.q0002
+++ b/tests/data/regression-1/fort.q0002
@@ -16,7 +16,7 @@
     0.1000000000000000E-02    dx
     0.1000000000000000E-02    dy
 
-     8                 grid_number
+     7                 grid_number
      2                 AMR_level
     40                 mx
     40                 my
@@ -25,7 +25,7 @@
     0.1000000000000000E-02    dx
     0.1000000000000000E-02    dy
 
-     7                 grid_number
+     8                 grid_number
      2                 AMR_level
     40                 mx
     40                 my
@@ -34,7 +34,7 @@
     0.1000000000000000E-02    dx
     0.1000000000000000E-02    dy
 
-     2                 grid_number
+     3                 grid_number
      2                 AMR_level
     40                 mx
     40                 my
diff --git a/tests/data/regression-1/fort.q0003 b/tests/data/regression-1/fort.q0003
index d476f7d..ff3144a 100644
--- a/tests/data/regression-1/fort.q0003
+++ b/tests/data/regression-1/fort.q0003
@@ -7,7 +7,7 @@
     0.4000000000000000E-02    dx
     0.4000000000000000E-02    dy
 
-     9                 grid_number
+     6                 grid_number
      2                 AMR_level
     48                 mx
     48                 my
@@ -16,7 +16,7 @@
     0.1000000000000000E-02    dx
     0.1000000000000000E-02    dy
 
-     7                 grid_number
+     4                 grid_number
      2                 AMR_level
     44                 mx
     48                 my
@@ -25,7 +25,7 @@
     0.1000000000000000E-02    dx
     0.1000000000000000E-02    dy
 
-     8                 grid_number
+     5                 grid_number
      2                 AMR_level
     48                 mx
     44                 my
diff --git a/tests/data/regression-1/fort.q0004 b/tests/data/regression-1/fort.q0004
index 72b6d06..48c09b6 100644
--- a/tests/data/regression-1/fort.q0004
+++ b/tests/data/regression-1/fort.q0004
@@ -7,7 +7,7 @@
     0.4000000000000000E-02    dx
     0.4000000000000000E-02    dy
 
-     9                 grid_number
+     6                 grid_number
      2                 AMR_level
     52                 mx
     52                 my
@@ -16,7 +16,7 @@
     0.1000000000000000E-02    dx
     0.1000000000000000E-02    dy
 
-     8                 grid_number
+     5                 grid_number
      2                 AMR_level
     52                 mx
     52                 my
@@ -25,7 +25,7 @@
     0.1000000000000000E-02    dx
     0.1000000000000000E-02    dy
 
-     7                 grid_number
+     4                 grid_number
      2                 AMR_level
     52                 mx
     52                 my
diff --git a/tests/data/regression-1/fort.q0005 b/tests/data/regression-1/fort.q0005
index 9b6d7aa..665a4d3 100644
--- a/tests/data/regression-1/fort.q0005
+++ b/tests/data/regression-1/fort.q0005
@@ -7,7 +7,7 @@
     0.4000000000000000E-02    dx
     0.4000000000000000E-02    dy
 
-     9                 grid_number
+     6                 grid_number
      2                 AMR_level
     56                 mx
     56                 my
@@ -16,7 +16,7 @@
     0.1000000000000000E-02    dx
     0.1000000000000000E-02    dy
 
-     8                 grid_number
+     5                 grid_number
      2                 AMR_level
     56                 mx
     56                 my
@@ -25,7 +25,7 @@
     0.1000000000000000E-02    dx
     0.1000000000000000E-02    dy
 
-     7                 grid_number
+     4                 grid_number
      2                 AMR_level
     56                 mx
     56                 my
diff --git a/tests/data/regression-1/frame00000.png b/tests/data/regression-1/frame00000.png
index af381e3..90a5474 100644
Binary files a/tests/data/regression-1/frame00000.png and b/tests/data/regression-1/frame00000.png differ
diff --git a/tests/data/regression-1/frame00001.png b/tests/data/regression-1/frame00001.png
index 29826af..90c2563 100644
Binary files a/tests/data/regression-1/frame00001.png and b/tests/data/regression-1/frame00001.png differ
diff --git a/tests/data/regression-1/frame00002.png b/tests/data/regression-1/frame00002.png
index 6450f2b..c551c12 100644
Binary files a/tests/data/regression-1/frame00002.png and b/tests/data/regression-1/frame00002.png differ
diff --git a/tests/data/regression-1/frame00003.png b/tests/data/regression-1/frame00003.png
index 3c94ceb..c054ff0 100644
Binary files a/tests/data/regression-1/frame00003.png and b/tests/data/regression-1/frame00003.png differ
diff --git a/tests/data/regression-1/frame00004.png b/tests/data/regression-1/frame00004.png
index f466a31..2c00d0a 100644
Binary files a/tests/data/regression-1/frame00004.png and b/tests/data/regression-1/frame00004.png differ
diff --git a/tests/data/regression-1/frame00005.png b/tests/data/regression-1/frame00005.png
index 781b81e..3df9e1d 100644
Binary files a/tests/data/regression-1/frame00005.png and b/tests/data/regression-1/frame00005.png differ
diff --git a/tests/data/regression-1/regression-1.nc b/tests/data/regression-1/regression-1.nc
index 07bd0cb..be9b7fa 100644
Binary files a/tests/data/regression-1/regression-1.nc and b/tests/data/regression-1/regression-1.nc differ
diff --git a/tests/data/regression-1/volumes.csv b/tests/data/regression-1/volumes.csv
index 7cec668..2d0141b 100644
--- a/tests/data/regression-1/volumes.csv
+++ b/tests/data/regression-1/volumes.csv
@@ -1,7 +1,7 @@
 frame,level 1,level 2
 0,0.0,0.0
-1,1.5280947247343468e-06,1.4799999999999992e-06
-2,2.9930331957813305e-06,2.959999999999998e-06
-3,4.477314971709564e-06,4.4400000000000295e-06
-4,5.977680238761498e-06,5.920000000000061e-06
-5,7.445095228902317e-06,7.4000000000000935e-06
+1,1.5238739454696168e-06,1.4800000755474338e-06
+2,2.994774725214481e-06,2.9600000755474357e-06
+3,4.472597858500439e-06,4.440000075547438e-06
+4,5.956159041901612e-06,5.920000075547439e-06
+5,7.438311361212e-06,7.400000075547437e-06
diff --git a/tests/data/regression-2/depth/frame00001.png b/tests/data/regression-2/depth/frame00001.png
index 633a3be..0507360 100644
Binary files a/tests/data/regression-2/depth/frame00001.png and b/tests/data/regression-2/depth/frame00001.png differ
diff --git a/tests/data/regression-2/depth/frame00002.png b/tests/data/regression-2/depth/frame00002.png
index 3069a02..ecf6e69 100644
Binary files a/tests/data/regression-2/depth/frame00002.png and b/tests/data/regression-2/depth/frame00002.png differ
diff --git a/tests/data/regression-2/depth/frame00003.png b/tests/data/regression-2/depth/frame00003.png
index ddf4729..ef7f7da 100644
Binary files a/tests/data/regression-2/depth/frame00003.png and b/tests/data/regression-2/depth/frame00003.png differ
diff --git a/tests/data/regression-2/depth/frame00004.png b/tests/data/regression-2/depth/frame00004.png
index eae9d03..9ed1128 100644
Binary files a/tests/data/regression-2/depth/frame00004.png and b/tests/data/regression-2/depth/frame00004.png differ
diff --git a/tests/data/regression-2/depth/frame00005.png b/tests/data/regression-2/depth/frame00005.png
index d851848..b95231f 100644
Binary files a/tests/data/regression-2/depth/frame00005.png and b/tests/data/regression-2/depth/frame00005.png differ
diff --git a/tests/data/regression-2/fort.a0001 b/tests/data/regression-2/fort.a0001
index aed25c6..f79f1ee 100644
Binary files a/tests/data/regression-2/fort.a0001 and b/tests/data/regression-2/fort.a0001 differ
diff --git a/tests/data/regression-2/fort.a0002 b/tests/data/regression-2/fort.a0002
index 8f75e32..75b7a0a 100644
Binary files a/tests/data/regression-2/fort.a0002 and b/tests/data/regression-2/fort.a0002 differ
diff --git a/tests/data/regression-2/fort.b0001 b/tests/data/regression-2/fort.b0001
index a84e39e..1eb9802 100644
Binary files a/tests/data/regression-2/fort.b0001 and b/tests/data/regression-2/fort.b0001 differ
diff --git a/tests/data/regression-2/fort.b0002 b/tests/data/regression-2/fort.b0002
index a41fdc9..c3101c2 100644
Binary files a/tests/data/regression-2/fort.b0002 and b/tests/data/regression-2/fort.b0002 differ
diff --git a/tests/data/regression-2/fort.b0003 b/tests/data/regression-2/fort.b0003
index 6f75cb7..3b37212 100644
Binary files a/tests/data/regression-2/fort.b0003 and b/tests/data/regression-2/fort.b0003 differ
diff --git a/tests/data/regression-2/fort.b0004 b/tests/data/regression-2/fort.b0004
index aea7850..49fa7cf 100644
Binary files a/tests/data/regression-2/fort.b0004 and b/tests/data/regression-2/fort.b0004 differ
diff --git a/tests/data/regression-2/fort.b0005 b/tests/data/regression-2/fort.b0005
index 03bd4ba..7025a5f 100644
Binary files a/tests/data/regression-2/fort.b0005 and b/tests/data/regression-2/fort.b0005 differ
diff --git a/tests/data/regression-2/fort.q0001 b/tests/data/regression-2/fort.q0001
index d55261d..18f3d46 100644
--- a/tests/data/regression-2/fort.q0001
+++ b/tests/data/regression-2/fort.q0001
@@ -7,9 +7,18 @@
     0.4000000000000000E+01    dx
     0.4000000000000000E+01    dy
 
-     2                 grid_number
+     5                 grid_number
      2                 AMR_level
-    56                 mx
+    32                 mx
+    44                 my
+    0.3200000000000000E+02    xlow
+    0.8000000000000000E+01    ylow
+    0.1000000000000000E+01    dx
+    0.1000000000000000E+01    dy
+
+     3                 grid_number
+     2                 AMR_level
+    28                 mx
     44                 my
     0.4000000000000000E+01    xlow
     0.8000000000000000E+01    ylow
diff --git a/tests/data/regression-2/fort.q0002 b/tests/data/regression-2/fort.q0002
index 19523fd..0d1f58d 100644
--- a/tests/data/regression-2/fort.q0002
+++ b/tests/data/regression-2/fort.q0002
@@ -10,16 +10,16 @@
      5                 grid_number
      2                 AMR_level
     44                 mx
-    48                 my
+    44                 my
     0.4400000000000000E+02    xlow
     0.8000000000000000E+01    ylow
     0.1000000000000000E+01    dx
     0.1000000000000000E+01    dy
 
-     2                 grid_number
+     3                 grid_number
      2                 AMR_level
     40                 mx
-    48                 my
+    44                 my
     0.4000000000000000E+01    xlow
     0.8000000000000000E+01    ylow
     0.1000000000000000E+01    dx
diff --git a/tests/data/regression-2/fort.q0003 b/tests/data/regression-2/fort.q0003
index 6fe6be9..608335e 100644
--- a/tests/data/regression-2/fort.q0003
+++ b/tests/data/regression-2/fort.q0003
@@ -7,7 +7,7 @@
     0.4000000000000000E+01    dx
     0.4000000000000000E+01    dy
 
-     4                 grid_number
+     5                 grid_number
      2                 AMR_level
     56                 mx
     48                 my
diff --git a/tests/data/regression-2/fort.q0004 b/tests/data/regression-2/fort.q0004
index 5acc0f2..2c2cce0 100644
--- a/tests/data/regression-2/fort.q0004
+++ b/tests/data/regression-2/fort.q0004
@@ -7,7 +7,7 @@
     0.4000000000000000E+01    dx
     0.4000000000000000E+01    dy
 
-     8                 grid_number
+     9                 grid_number
      2                 AMR_level
     36                 mx
     48                 my
@@ -16,7 +16,7 @@
     0.1000000000000000E+01    dx
     0.1000000000000000E+01    dy
 
-     4                 grid_number
+     7                 grid_number
      2                 AMR_level
     32                 mx
     48                 my
@@ -25,7 +25,7 @@
     0.1000000000000000E+01    dx
     0.1000000000000000E+01    dy
 
-     7                 grid_number
+     3                 grid_number
      2                 AMR_level
     32                 mx
     48                 my
@@ -34,7 +34,7 @@
     0.1000000000000000E+01    dx
     0.1000000000000000E+01    dy
 
-     3                 grid_number
+     5                 grid_number
      2                 AMR_level
     32                 mx
     48                 my
diff --git a/tests/data/regression-2/fort.q0005 b/tests/data/regression-2/fort.q0005
index 6b81ca1..4f4b755 100644
--- a/tests/data/regression-2/fort.q0005
+++ b/tests/data/regression-2/fort.q0005
@@ -16,7 +16,7 @@
     0.1000000000000000E+01    dx
     0.1000000000000000E+01    dy
 
-     2                 grid_number
+     3                 grid_number
      2                 AMR_level
     32                 mx
     52                 my
@@ -25,7 +25,7 @@
     0.1000000000000000E+01    dx
     0.1000000000000000E+01    dy
 
-     6                 grid_number
+     7                 grid_number
      2                 AMR_level
     32                 mx
     52                 my
diff --git a/tests/data/regression-2/fort.t0001 b/tests/data/regression-2/fort.t0001
index a4c1b68..0c6fdfb 100644
--- a/tests/data/regression-2/fort.t0001
+++ b/tests/data/regression-2/fort.t0001
@@ -1,6 +1,6 @@
     0.13000000E+02    time
      4                 meqn
-     2                 ngrids
+     3                 ngrids
      2                 naux
      2                 ndim
      2                 nghost
diff --git a/tests/data/regression-2/regression-2.nc b/tests/data/regression-2/regression-2.nc
index 173cd05..29f5ed0 100644
Binary files a/tests/data/regression-2/regression-2.nc and b/tests/data/regression-2/regression-2.nc differ
diff --git a/tests/data/regression-2/removed_fluid.csv b/tests/data/regression-2/removed_fluid.csv
index 8ac673d..5fe3101 100644
--- a/tests/data/regression-2/removed_fluid.csv
+++ b/tests/data/regression-2/removed_fluid.csv
@@ -1,50 +1,50 @@
-123.5, 18.5, 63.5,  3.668478925962791E-005,  3.668478925962791E-005
-123.5, 19.5, 61.8,  9.961772121586661E-004,  9.961772121586661E-004
-123.5, 20.5, 60.0,  4.431936843537299E-003,  4.431936843537299E-003
-123.5, 21.5, 58.7,  1.064163356293725E-002,  1.064163356293725E-002
-123.5, 22.5, 57.6,  1.927188400688665E-002,  1.927188400688665E-002
-123.5, 23.5, 56.5,  3.009516315905580E-002,  3.009516315905580E-002
-123.5, 24.5, 55.5,  4.308024054375853E-002,  4.308024054375853E-002
-123.5, 25.5, 54.6,  5.804303344140620E-002,  5.804303344140620E-002
-123.5, 26.5, 53.9,  7.457786243877053E-002,  7.457786243877053E-002
-123.5, 27.5, 53.1,  9.271841963928681E-002,  9.271841963928681E-002
-123.5, 28.5, 52.4,  1.127352402468812E-001,  1.127352402468812E-001
-123.5, 29.5, 52.0,  1.349990508091892E-001,  1.349990508091892E-001
-123.5, 30.5, 52.0,  1.450560213027895E-001,  1.450560213027895E-001
-123.5, 31.5, 52.0,  1.320412956013401E-001,  1.320412956013401E-001
-123.5, 32.5, 52.4,  1.105216188357969E-001,  1.105216188357969E-001
-123.5, 33.5, 53.3,  9.103580858396133E-002,  9.103580858396133E-002
-123.5, 34.5, 53.9,  7.316729688034768E-002,  7.316729688034768E-002
-123.5, 35.5, 54.8,  5.684898881537943E-002,  5.684898881537943E-002
-123.5, 36.5, 55.5,  4.208409380158845E-002,  4.208409380158845E-002
-123.5, 37.5, 56.5,  2.930904432315403E-002,  2.930904432315403E-002
-123.5, 38.5, 57.6,  1.871083022469346E-002,  1.871083022469346E-002
-123.5, 39.5, 58.7,  1.024599501088745E-002,  1.024599501088745E-002
-123.5, 40.5, 60.2,  4.167540542684677E-003,  4.167540542684677E-003
-123.5, 41.5, 62.0,  8.967553004447496E-004,  8.967553004447496E-004
-123.5, 42.5, 63.7,  3.011500220659512E-005,  3.011500220659512E-005
-124.5, 18.5, 63.5,  3.668478925962791E-005,  3.668478925962791E-005
-124.5, 19.5, 61.8,  9.961772121586661E-004,  9.961772121586661E-004
-124.5, 20.5, 60.0,  4.431936843537299E-003,  4.431936843537299E-003
-124.5, 21.5, 58.7,  1.064163356293725E-002,  1.064163356293725E-002
-124.5, 22.5, 57.6,  1.927188400688665E-002,  1.927188400688665E-002
-124.5, 23.5, 56.5,  3.009516315905580E-002,  3.009516315905580E-002
-124.5, 24.5, 55.5,  4.308024054375853E-002,  4.308024054375853E-002
-124.5, 25.5, 54.6,  5.804303344140620E-002,  5.804303344140620E-002
-124.5, 26.5, 53.9,  7.457786243877053E-002,  7.457786243877053E-002
-124.5, 27.5, 53.1,  9.271841963928681E-002,  9.271841963928681E-002
-124.5, 28.5, 52.4,  1.127352402468812E-001,  1.127352402468812E-001
-124.5, 29.5, 52.0,  1.349990508091892E-001,  1.349990508091892E-001
-124.5, 30.5, 52.0,  1.450560213027895E-001,  1.450560213027895E-001
-124.5, 31.5, 52.0,  1.320412956013401E-001,  1.320412956013401E-001
-124.5, 32.5, 52.4,  1.105216188357969E-001,  1.105216188357969E-001
-124.5, 33.5, 53.3,  9.103580858396133E-002,  9.103580858396133E-002
-124.5, 34.5, 53.9,  7.316729688034768E-002,  7.316729688034768E-002
-124.5, 35.5, 54.8,  5.684898881537943E-002,  5.684898881537943E-002
-124.5, 36.5, 55.5,  4.208409380158845E-002,  4.208409380158845E-002
-124.5, 37.5, 56.5,  2.930904432315403E-002,  2.930904432315403E-002
-124.5, 38.5, 57.6,  1.871083022469346E-002,  1.871083022469346E-002
-124.5, 39.5, 58.7,  1.024599501088745E-002,  1.024599501088745E-002
-124.5, 40.5, 60.2,  4.167540542684677E-003,  4.167540542684677E-003
-124.5, 41.5, 62.0,  8.967553004447496E-004,  8.967553004447496E-004
-124.5, 42.5, 63.7,  3.011500220659512E-005,  3.011500220659512E-005
+123.5, 18.5, 64.1,  1.256902257230944E-005,  1.256902257230944E-005
+123.5, 19.5, 62.5,  6.747228055041653E-004,  6.747228055041653E-004
+123.5, 20.5, 60.8,  3.624994585755931E-003,  3.624994585755931E-003
+123.5, 21.5, 59.4,  9.601612922003750E-003,  9.601612922003750E-003
+123.5, 22.5, 58.2,  1.805194639732591E-002,  1.805194639732591E-002
+123.5, 23.5, 57.3,  2.852189375620896E-002,  2.852189375620896E-002
+123.5, 24.5, 56.3,  4.090493276212451E-002,  4.090493276212451E-002
+123.5, 25.5, 55.3,  5.520602562905776E-002,  5.520602562905776E-002
+123.5, 26.5, 54.5,  7.131434851715725E-002,  7.131434851715725E-002
+123.5, 27.5, 54.0,  8.935426915935542E-002,  8.935426915935542E-002
+123.5, 28.5, 53.0,  1.099244080518900E-001,  1.099244080518900E-001
+123.5, 29.5, 52.6,  1.318461654343450E-001,  1.318461654343450E-001
+123.5, 30.5, 52.6,  1.405353562160279E-001,  1.405353562160279E-001
+123.5, 31.5, 52.8,  1.282096821652821E-001,  1.282096821652821E-001
+123.5, 32.5, 53.2,  1.074104536846571E-001,  1.074104536846571E-001
+123.5, 33.5, 54.0,  8.738164861312597E-002,  8.738164861312597E-002
+123.5, 34.5, 54.7,  6.968934401172296E-002,  6.968934401172296E-002
+123.5, 35.5, 55.5,  5.390427220879378E-002,  5.390427220879378E-002
+123.5, 36.5, 56.3,  3.991132139262153E-002,  3.991132139262153E-002
+123.5, 37.5, 57.3,  2.777860753596639E-002,  2.777860753596639E-002
+123.5, 38.5, 58.4,  1.750622927131418E-002,  1.750622927131418E-002
+123.5, 39.5, 59.6,  9.215732155515861E-003,  9.215732155515861E-003
+123.5, 40.5, 61.0,  3.392269209949621E-003,  3.392269209949621E-003
+123.5, 41.5, 62.7,  5.919076313246389E-004,  5.919076313246389E-004
+123.5, 42.5, 64.3,  8.366893364293919E-006,  8.366893364293919E-006
+124.5, 18.5, 64.1,  1.256902257230944E-005,  1.256902257230944E-005
+124.5, 19.5, 62.5,  6.747228055041653E-004,  6.747228055041653E-004
+124.5, 20.5, 60.8,  3.624994585755931E-003,  3.624994585755931E-003
+124.5, 21.5, 59.4,  9.601612922003750E-003,  9.601612922003750E-003
+124.5, 22.5, 58.2,  1.805194639732591E-002,  1.805194639732591E-002
+124.5, 23.5, 57.3,  2.852189375620896E-002,  2.852189375620896E-002
+124.5, 24.5, 56.3,  4.090493276212451E-002,  4.090493276212451E-002
+124.5, 25.5, 55.3,  5.520602562905776E-002,  5.520602562905776E-002
+124.5, 26.5, 54.5,  7.131434851715725E-002,  7.131434851715725E-002
+124.5, 27.5, 54.0,  8.935426915935542E-002,  8.935426915935542E-002
+124.5, 28.5, 53.0,  1.099244080518900E-001,  1.099244080518900E-001
+124.5, 29.5, 52.6,  1.318461654343450E-001,  1.318461654343450E-001
+124.5, 30.5, 52.6,  1.405353562160279E-001,  1.405353562160279E-001
+124.5, 31.5, 52.8,  1.282096821652821E-001,  1.282096821652821E-001
+124.5, 32.5, 53.2,  1.074104536846571E-001,  1.074104536846571E-001
+124.5, 33.5, 54.0,  8.738164861312597E-002,  8.738164861312597E-002
+124.5, 34.5, 54.7,  6.968934401172296E-002,  6.968934401172296E-002
+124.5, 35.5, 55.5,  5.390427220879378E-002,  5.390427220879378E-002
+124.5, 36.5, 56.3,  3.991132139262153E-002,  3.991132139262153E-002
+124.5, 37.5, 57.3,  2.777860753596639E-002,  2.777860753596639E-002
+124.5, 38.5, 58.4,  1.750622927131418E-002,  1.750622927131418E-002
+124.5, 39.5, 59.6,  9.215732155515861E-003,  9.215732155515861E-003
+124.5, 40.5, 61.0,  3.392269209949621E-003,  3.392269209949621E-003
+124.5, 41.5, 62.7,  5.919076313246389E-004,  5.919076313246389E-004
+124.5, 42.5, 64.3,  8.366893364293919E-006,  8.366893364293919E-006
diff --git a/tests/data/regression-2/topo/frame00001.png b/tests/data/regression-2/topo/frame00001.png
index 263a1ac..3621219 100644
Binary files a/tests/data/regression-2/topo/frame00001.png and b/tests/data/regression-2/topo/frame00001.png differ
diff --git a/tests/data/regression-2/topo/frame00002.png b/tests/data/regression-2/topo/frame00002.png
index 444e1c7..7231652 100644
Binary files a/tests/data/regression-2/topo/frame00002.png and b/tests/data/regression-2/topo/frame00002.png differ
diff --git a/tests/data/regression-2/volumes.csv b/tests/data/regression-2/volumes.csv
index f982492..43b7f04 100644
--- a/tests/data/regression-2/volumes.csv
+++ b/tests/data/regression-2/volumes.csv
@@ -1,7 +1,7 @@
 frame,level 1,level 2
 0,0.0,0.0
-1,6.1910603304118474,6.782206658433722
-2,7.885807185595178,13.282840527703026
-3,8.137599209763328,19.784740872555574
-4,12.646051835409626,26.287573889036814
-5,15.372926000275601,30.198551957488988
+1,6.16501724037713,6.51215195738637
+2,7.657355617762287,13.013088065081753
+3,8.33010238488509,19.513117559729356
+4,12.87672163498295,26.015447289876143
+5,15.84968233014592,30.027687960401312
diff --git a/tests/test_regression_1.py b/tests/test_regression_1.py
index 21b3b34..83166f8 100644
--- a/tests/test_regression_1.py
+++ b/tests/test_regression_1.py
@@ -8,11 +8,13 @@
 
 """Regression test 1."""
 # pylint: disable=redefined-outer-name
+import os
 import sys
 import csv
 import pathlib
 import pytest
 import rasterio
+import matplotlib
 import matplotlib.pyplot
 import gclandspill.__main__
 import gclandspill.pyclaw
@@ -22,54 +24,30 @@
 def create_case(tmp_path_factory):
     """A pytest fixture to make a temp case available across different tests."""
 
+    # force to use only 4 threads
+    os.environ["OMP_NUM_THREADS"] = "4"
+
+    # to avoid wayland session or none-X environment
+    matplotlib.use("agg")
+
     case_dir = tmp_path_factory.mktemp("regression-test-1")
 
     content = (
         "#! /usr/bin/env python\n"
         "import gclandspill\n"
-        "def setrun(prog='geoclaw'):\n"
-        "\trundata = gclandspill.clawutil.data.ClawRunData('geoclaw', 2)\n"
-        "\trundata.add_data(gclandspill.data.LandSpillData(), 'landspill_data')\n"
-        "\trundata.clawdata.num_dim = 2\n"
+        "def setrun():\n"
+        "\trundata = gclandspill.data.ClawRunData()\n"
         "\trundata.clawdata.lower[:] = [-0.1, -0.1]\n"
         "\trundata.clawdata.upper[:] = [0.1, 0.1]\n"
         "\trundata.clawdata.num_cells[:] = [50, 50]\n"
-        "\trundata.clawdata.num_eqn = 3\n"
-        "\trundata.clawdata.num_aux = 2\n"
         "\trundata.clawdata.output_style = 2\n"
         "\trundata.clawdata.output_times = list(range(6))\n"
-        "\trundata.clawdata.output_format = 'binary'\n"
-        "\trundata.clawdata.dt_variable = True\n"
-        "\trundata.clawdata.dt_initial = 1e-4\n"
-        "\trundata.clawdata.dt_max = 0.01\n"
-        "\trundata.clawdata.transverse_waves = 2\n"
-        "\trundata.clawdata.num_waves = 3\n"
-        "\trundata.clawdata.limiter = ['mc', 'mc', 'mc']\n"
-        "\trundata.clawdata.use_fwaves = True\n"
-        "\trundata.clawdata.source_split = 'godunov'\n"
-        "\trundata.clawdata.num_ghost = 2\n"
-        "\trundata.clawdata.bc_lower[:] = [1, 1]\n"
-        "\trundata.clawdata.bc_upper[:] = [1, 1]\n"
-        "\trundata.amrdata.amr_levels_max = 2\n"
-        "\trundata.amrdata.refinement_ratios_x = [4]\n"
-        "\trundata.amrdata.refinement_ratios_y = [4]\n"
-        "\trundata.amrdata.refinement_ratios_t = [4]\n"
-        "\trundata.amrdata.aux_type = ['center', 'center']\n"
-        "\trundata.geo_data.gravity = 9.81\n"
-        "\trundata.geo_data.coriolis_forcing = False\n"
-        "\trundata.geo_data.sea_level = 0.0\n"
-        "\trundata.geo_data.dry_tolerance = 1.e-4\n"
-        "\trundata.geo_data.friction_forcing = False\n"
-        "\trundata.refinement_data.wave_tolerance = 1.e-5\n"
-        "\trundata.refinement_data.speed_tolerance = [1e-8]\n"
-        "\trundata.refinement_data.variable_dt_refinement_ratios = True\n"
+        "\trundata.clawdata.dt_initial = 0.005\n"
+        "\trundata.clawdata.dt_max = 0.02\n"
         "\trundata.topo_data.topofiles.append([3, 'toy.asc'])\n"
-        "\trundata.landspill_data.update_tol = rundata.geo_data.dry_tolerance\n"
         "\trundata.landspill_data.point_sources.n_point_sources = 1\n"
         "\trundata.landspill_data.point_sources.point_sources.append([[0., 0.], 1, [3600.], [1.48e-6]])\n"
         "\trundata.landspill_data.darcy_weisbach_friction.type = 4\n"
-        "\trundata.landspill_data.darcy_weisbach_friction.dry_tol = 1e-4\n"
-        "\trundata.landspill_data.darcy_weisbach_friction.friction_tol = 1e6\n"
         "\trundata.landspill_data.darcy_weisbach_friction.default_roughness = 0.0\n"
         "\trundata.landspill_data.darcy_weisbach_friction.filename = 'roughness.asc'\n"
         "\trundata.landspill_data.evaporation.type = 1\n"
@@ -136,7 +114,7 @@ def test_raw_result(create_case, frame):
     case_dir = create_case
 
     ref = gclandspill.pyclaw.Solution()
-    ref.read(frame, pathlib.Path(__file__).parent.joinpath("data", "regression-1"), "binary",  read_aux=False)
+    ref.read(frame, pathlib.Path(__file__).parent.joinpath("data", "regression-1"), "binary",  read_aux=True)
 
     soln = gclandspill.pyclaw.Solution()
     soln.read(frame, case_dir.joinpath("_output"), "binary",  read_aux=False)
diff --git a/tests/test_regression_2.py b/tests/test_regression_2.py
index 2cf2809..2ff9294 100644
--- a/tests/test_regression_2.py
+++ b/tests/test_regression_2.py
@@ -8,12 +8,14 @@
 
 """Regression test 2."""
 # pylint: disable=redefined-outer-name
+import os
 import sys
 import csv
 import pathlib
 import numpy
 import pytest
 import rasterio
+import matplotlib
 import matplotlib.pyplot
 import gclandspill.__main__
 import gclandspill.pyclaw
@@ -24,62 +26,34 @@ def create_case(tmp_path_factory):
     """A pytest fixture to make a temp case available across different tests."""
     # pylint: disable=invalid-name
 
+    # force to use only 4 threads
+    os.environ["OMP_NUM_THREADS"] = "4"
+
+    # to avoid wayland session or none-X environment
+    matplotlib.use("agg")
+
     case_dir = tmp_path_factory.mktemp("regression-test-2")
 
     content = (
         "#! /usr/bin/env python\n"
         "import gclandspill\n"
-        "def setrun(prog='geoclaw'):\n"
-        "\trundata = gclandspill.clawutil.data.ClawRunData('geoclaw', 2)\n"
-        "\trundata.add_data(gclandspill.data.LandSpillData(), 'landspill_data')\n"
-        "\trundata.clawdata.num_dim = 2\n"
+        "def setrun():\n"
+        "\trundata = gclandspill.data.ClawRunData()\n"
         "\trundata.clawdata.lower[:] = [0.0, 0.0]\n"
         "\trundata.clawdata.upper[:] = [152.0, 60.0]\n"
         "\trundata.clawdata.num_cells[:] = [38, 15]\n"
-        "\trundata.clawdata.num_eqn = 3\n"
-        "\trundata.clawdata.num_aux = 2\n"
         "\trundata.clawdata.output_style = 1\n"
         "\trundata.clawdata.num_output_times = 5\n"
         "\trundata.clawdata.tfinal = 65\n"
         "\trundata.clawdata.output_t0 = True\n"
-        "\trundata.clawdata.output_format = 'binary'\n"
-        "\trundata.clawdata.output_aux_components = 'all'\n"
-        "\trundata.clawdata.output_aux_onlyonce = False\n"
-        "\trundata.clawdata.dt_variable = True\n"
         "\trundata.clawdata.dt_initial = 0.6\n"
         "\trundata.clawdata.dt_max = 4.0\n"
-        "\trundata.clawdata.transverse_waves = 2\n"
-        "\trundata.clawdata.num_waves = 3\n"
-        "\trundata.clawdata.limiter = ['mc', 'mc', 'mc']\n"
-        "\trundata.clawdata.use_fwaves = True\n"
-        "\trundata.clawdata.source_split = 'godunov'\n"
-        "\trundata.clawdata.num_ghost = 2\n"
-        "\trundata.clawdata.bc_lower[:] = [1, 1]\n"
-        "\trundata.clawdata.bc_upper[:] = [1, 1]\n"
-        "\trundata.amrdata.amr_levels_max = 2\n"
-        "\trundata.amrdata.refinement_ratios_x = [4]\n"
-        "\trundata.amrdata.refinement_ratios_y = [4]\n"
-        "\trundata.amrdata.refinement_ratios_t = [4]\n"
-        "\trundata.amrdata.aux_type = ['center', 'center']\n"
-        "\trundata.geo_data.gravity = 9.81\n"
-        "\trundata.geo_data.coriolis_forcing = False\n"
-        "\trundata.geo_data.sea_level = -10.0\n"
-        "\trundata.geo_data.dry_tolerance = 1.e-4\n"
-        "\trundata.geo_data.friction_forcing = False\n"
-        "\trundata.refinement_data.wave_tolerance = 1.e-5\n"
-        "\trundata.refinement_data.speed_tolerance = [1e-8]\n"
-        "\trundata.refinement_data.variable_dt_refinement_ratios = True\n"
         "\trundata.topo_data.topofiles.append([3, 'topo.asc'])\n"
         "\trundata.landspill_data.ref_mu = 0.6512\n"
-        "\trundata.landspill_data.ref_temperature = 15.\n"
-        "\trundata.landspill_data.ambient_temperature = 25.\n"
         "\trundata.landspill_data.density = 800.\n"
-        "\trundata.landspill_data.update_tol = rundata.geo_data.dry_tolerance\n"
         "\trundata.landspill_data.point_sources.n_point_sources = 1\n"
         "\trundata.landspill_data.point_sources.point_sources.append([[20., 30.], 1, [1800.], [0.5]])\n"
         "\trundata.landspill_data.darcy_weisbach_friction.type = 4\n"
-        "\trundata.landspill_data.darcy_weisbach_friction.dry_tol = 1e-4\n"
-        "\trundata.landspill_data.darcy_weisbach_friction.friction_tol = 1e6\n"
         "\trundata.landspill_data.darcy_weisbach_friction.default_roughness = 0.0\n"
         "\trundata.landspill_data.darcy_weisbach_friction.filename = 'roughness.asc'\n"
         "\trundata.landspill_data.hydro_features.files.append('hydro.asc')\n"
@@ -192,7 +166,7 @@ def test_raw_result(create_case, frame):
     case_dir = create_case
 
     ref = gclandspill.pyclaw.Solution()
-    ref.read(frame, pathlib.Path(__file__).parent.joinpath("data", "regression-2"), "binary",  read_aux=False)
+    ref.read(frame, pathlib.Path(__file__).parent.joinpath("data", "regression-2"), "binary",  read_aux=True)
 
     soln = gclandspill.pyclaw.Solution()
     soln.read(frame, case_dir.joinpath("_output"), "binary",  read_aux=False)