Repository: EC-USGS/pywatershed (language: python)

Issue:
update docs for release
<!-- Feel free to remove check-list items aren't relevant to your change -->

- [x] User visible changes (including notable bug fixes) are documented in `whats-new.rst`
- [x] Refresh/polish notebooks.
- [x] General improvement of the RDT pages.


Discussion / hints:

Fix patch:
diff --git a/.github/workflows/ci.yaml b/.github/workflows/ci.yaml
index ed25a84c..a2af53e8 100644
--- a/.github/workflows/ci.yaml
+++ b/.github/workflows/ci.yaml
@@ -175,10 +175,11 @@ jobs:
           pip list
 
       - name: Run available domains with PRMS and convert csv output to NetCDF
-        working-directory: test_data/scripts
+        working-directory: test_data/generate
         run: |
-          pytest -v -n=auto --durations=0 test_run_domains.py
-          pytest -v -n=auto --durations=0 test_nc_domains.py
+          pytest -v -n=auto --durations=0 run_prms_domains.py
+          pytest -v -n=auto --durations=0 convert_prms_output_to_nc.py
+          pytest -v -n=auto --durations=0 remove_prms_csvs.py
 
       - name: List all NetCDF files in test_data directory
         working-directory: test_data
@@ -198,13 +199,13 @@ jobs:
 
       - name: Upload test results
         if: always()
-        uses: actions/upload-artifact@v2
+        uses: actions/upload-artifact@v3
         with:
           name: Test results for ${{ runner.os }}-${{ matrix.python-version }}
           path: ./autotest/pytest.xml
 
       - name: Upload code coverage to Codecov
-        uses: codecov/codecov-action@v2.1.0
+        uses: codecov/codecov-action@v3
         with:
           file: ./autotest/coverage.xml
           # flags: unittests
diff --git a/.gitignore b/.gitignore
index 8f0e88e1..4c7f1c3f 100644
--- a/.gitignore
+++ b/.gitignore
@@ -61,12 +61,15 @@ generated/
 
 # example output data and gis data
 examples/0*/
+examples/snow_errors
+examples/runoff_errors
 examples/pynhm_gis
 examples/pywatershed_gis
 examples/model_loop_custom_output
 pywatershed/data/pynhm_gis
 pywatershed/data/pywatershed_gis
 
+
 # graphics
 *.png
 *.svg
diff --git a/bin/prms_linux_gfort_dbl_prec b/bin/prms_linux_gfort_dbl_prec
new file mode 100755
index 00000000..00b2ab53
Binary files /dev/null and b/bin/prms_linux_gfort_dbl_prec differ
diff --git a/bin/prms_linux b/bin/prms_linux_gfort_mixed_prec
similarity index 100%
rename from bin/prms_linux
rename to bin/prms_linux_gfort_mixed_prec
diff --git a/bin/prms_mac_intel_gfort_dbl_prec b/bin/prms_mac_intel_gfort_dbl_prec
new file mode 100755
index 00000000..9abb32c4
Binary files /dev/null and b/bin/prms_mac_intel_gfort_dbl_prec differ
diff --git a/bin/prms_mac b/bin/prms_mac_intel_gfort_mixed_prec
similarity index 100%
rename from bin/prms_mac
rename to bin/prms_mac_intel_gfort_mixed_prec
diff --git a/bin/prms_mac_m1_ifort_dbl_prec b/bin/prms_mac_m1_ifort_dbl_prec
new file mode 100755
index 00000000..a9e1ee79
Binary files /dev/null and b/bin/prms_mac_m1_ifort_dbl_prec differ
diff --git a/bin/prms_mac_m1_intel b/bin/prms_mac_m1_ifort_mixed_prec
similarity index 99%
rename from bin/prms_mac_m1_intel
rename to bin/prms_mac_m1_ifort_mixed_prec
index 7485b20d..56403fae 100755
Binary files a/bin/prms_mac_m1_intel and b/bin/prms_mac_m1_ifort_mixed_prec differ
diff --git a/bin/prms_win.exe b/bin/prms_win.exe
deleted file mode 100755
index dfd2476e..00000000
Binary files a/bin/prms_win.exe and /dev/null differ
diff --git a/bin/prms_win_gfort_dbl_prec.exe b/bin/prms_win_gfort_dbl_prec.exe
new file mode 100644
index 00000000..9ac58623
Binary files /dev/null and b/bin/prms_win_gfort_dbl_prec.exe differ
diff --git a/examples/runoff_errors.ipynb b/examples/runoff_errors.ipynb
new file mode 100644
index 00000000..c49d5172
--- /dev/null
+++ b/examples/runoff_errors.ipynb
@@ -0,0 +1,824 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "id": "4bf60102-e977-4b2c-82b1-86780c283c71",
+   "metadata": {},
+   "source": [
+    "# Runoff mass balance errors\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "6936fb48-fb5c-463e-8642-11287a578ab1",
+   "metadata": {
+    "editable": true,
+    "slideshow": {
+     "slide_type": ""
+    }
+   },
+   "outputs": [],
+   "source": [
+    "# auto-format the code in this notebook\n",
+    "%load_ext jupyter_black"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "c7b84519-933d-4284-838d-f5a8d4ec01f7",
+   "metadata": {},
+   "source": [
+    "## Setup"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "d316a5de-19ef-4008-8c68-1e192c3f89d9",
+   "metadata": {
+    "editable": true,
+    "slideshow": {
+     "slide_type": ""
+    }
+   },
+   "outputs": [],
+   "source": [
+    "import pathlib as pl\n",
+    "from pprint import pprint\n",
+    "from shutil import rmtree, copy2\n",
+    "\n",
+    "import hvplot.xarray  # noqa\n",
+    "from icecream import ic\n",
+    "from IPython.display import display\n",
+    "import numpy as np\n",
+    "import pywatershed as pws\n",
+    "import xarray as xr"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "9c188b5b-fc91-480c-a1a2-ce1121f92915",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Configuration\n",
+    "domain_name = \"drb_2yr\"\n",
+    "\n",
+    "nb_output_dir = pl.Path(\"./runoff_errors\")\n",
+    "\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "b9e616a3-03a3-4a13-9a33-181af59f1c71",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "pws_root = pws.constants.__pywatershed_root__\n",
+    "domain_dir = pws_root / f\"../test_data/{domain_name}\"\n",
+    "nb_output_dir.mkdir(exist_ok=True)\n",
+    "\n",
+    "zero = np.zeros([1])[0]\n",
+    "epsilon64 = np.finfo(zero).eps\n",
+    "epsilon32 = np.finfo(zero.astype(\"float32\")).eps"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "0a05c33a-c86b-433a-a400-54c7f2df8eb9",
+   "metadata": {},
+   "source": [
+    "## Run PRMS double precision runs and convert to netcdf"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "daa94738-1c8b-4113-a529-c6dd0b4a0284",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "bin_dir = pws_root / \"../bin/\"\n",
+    "# bin_mixed = bin_dir / \"prms_mac_m1_ifort_mixed_prec\"\n",
+    "bin_double = bin_dir / \"prms_mac_m1_ifort_dbl_prec\""
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "611c0d7a-7bfd-422d-ae53-3978b0c2eb83",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def run_prms(binary: pl.Path, run_dir: pl.Path, skip_if_exists=False):\n",
+    "    import shlex\n",
+    "    import subprocess\n",
+    "\n",
+    "    from pywatershed import CsvFile, Soltab\n",
+    "    from pywatershed.parameters import PrmsParameters\n",
+    "\n",
+    "    if skip_if_exists and run_dir.exists():\n",
+    "        print(\n",
+    "            f\"Run ({run_dir}) already exists and skip_if_exists=True. Using existing run.\"\n",
+    "        )\n",
+    "        return None\n",
+    "    \n",
+    "    run_dir.mkdir()  # must not exist, on user to delete\n",
+    "    copy2(binary, run_dir / binary.name)\n",
+    "    for ff in [\n",
+    "        \"control.test\",\n",
+    "        \"myparam.param\",\n",
+    "        \"tmax.cbh\",\n",
+    "        \"tmin.cbh\",\n",
+    "        \"prcp.cbh\",\n",
+    "        \"sf_data\",\n",
+    "    ]:\n",
+    "        copy2(domain_dir / ff, run_dir / ff)\n",
+    "\n",
+    "    output_dir = run_dir / \"output\"\n",
+    "    output_dir.mkdir()\n",
+    "\n",
+    "    exe_command = f\"time ./{binary.name} control.test -MAXDATALNLEN 60000 2>&1 | tee run.log\"\n",
+    "    result = subprocess.run(\n",
+    "        exe_command,\n",
+    "        shell=True,\n",
+    "        # stdout = subprocess.PIPE,\n",
+    "        stderr=subprocess.STDOUT,\n",
+    "        universal_newlines=True,\n",
+    "        cwd=run_dir,\n",
+    "    )\n",
+    "\n",
+    "    # these will be useful in what follows    \n",
+    "    params = pws.parameters.PrmsParameters.load(\n",
+    "        domain_dir / \"myparam.param\"\n",
+    "    ).parameters\n",
+    "    \n",
+    "    # convert to netcdf\n",
+    "    # could make these arguments\n",
+    "    chunking = {\n",
+    "        \"time\": 0,\n",
+    "        \"doy\": 0,\n",
+    "        \"nhm_id\": 0,\n",
+    "        \"nhm_seg\": 0,\n",
+    "    }\n",
+    "\n",
+    "    output_csvs = output_dir.glob(\"*.csv\")\n",
+    "    for cc in output_csvs:\n",
+    "        if cc.name in [\"stats.csv\"]:\n",
+    "            continue\n",
+    "        nc_path = cc.with_suffix(\".nc\")\n",
+    "        CsvFile(cc).to_netcdf(nc_path, chunk_sizes=chunking)\n",
+    "\n",
+    "    # previous and change variables\n",
+    "    for vv in [\n",
+    "        \"pk_ice\",\n",
+    "        \"freeh2o\",\n",
+    "        \"soil_moist\",\n",
+    "        \"hru_impervstor\",\n",
+    "        \"dprst_stor_hru\",\n",
+    "    ]:\n",
+    "        data = xr.open_dataset(output_dir / f\"{vv}.nc\")[vv]\n",
+    "        prev_da = data.copy()\n",
+    "        prev_da[:] = np.roll(prev_da.values, 1, axis=0)\n",
+    "        assert (prev_da[1:, :].values == data[0:-1, :].values).all()\n",
+    "        prev_da[0, :] = np.zeros(1)[\n",
+    "            0\n",
+    "        ]  # np.nan better but causes plotting to fail\n",
+    "        change_da = data - prev_da\n",
+    "        prev_da.rename(f\"{vv}_prev\").to_dataset().to_netcdf(\n",
+    "            output_dir / f\"{vv}_prev.nc\"\n",
+    "        )\n",
+    "        data[f\"{vv}_prev\"] = xr.open_dataset(output_dir / f\"{vv}_prev.nc\")[\n",
+    "            f\"{vv}_prev\"\n",
+    "        ]\n",
+    "\n",
+    "        change_da.rename(f\"{vv}_change\").to_dataset().to_netcdf(\n",
+    "            output_dir / f\"{vv}_change.nc\"\n",
+    "        )\n",
+    "        data[f\"{vv}_change\"] = xr.open_dataset(output_dir / f\"{vv}_change.nc\")[\n",
+    "            f\"{vv}_change\"\n",
+    "        ]\n",
+    "\n",
+    "    # through_rain\n",
+    "    dep_vars = [\n",
+    "        \"net_ppt\",\n",
+    "        \"pptmix_nopack\",\n",
+    "        \"snowmelt\",\n",
+    "        \"pkwater_equiv\",\n",
+    "        \"snow_evap\",\n",
+    "        \"net_snow\",\n",
+    "        \"net_rain\",\n",
+    "    ]\n",
+    "    data = {}\n",
+    "    for vv in dep_vars:\n",
+    "        data[vv] = xr.open_dataset(output_dir / f\"{vv}.nc\")[vv]\n",
+    "\n",
+    "    nearzero = 1.0e-6\n",
+    "\n",
+    "    cond1 = data[\"net_ppt\"] > zero\n",
+    "    cond2 = data[\"pptmix_nopack\"] != 0\n",
+    "    cond3 = data[\"snowmelt\"] < nearzero\n",
+    "    cond4 = data[\"pkwater_equiv\"] < epsilon32\n",
+    "    cond5 = data[\"snow_evap\"] < nearzero\n",
+    "    cond6 = data[\"net_snow\"] < nearzero\n",
+    "\n",
+    "    through_rain = data[\"net_rain\"] * zero\n",
+    "    # these are in reverse order\n",
+    "    through_rain[:] = np.where(\n",
+    "        cond1 & cond3 & cond4 & cond6, data[\"net_rain\"], zero\n",
+    "    )\n",
+    "    through_rain[:] = np.where(\n",
+    "        cond1 & cond3 & cond4 & cond5, data[\"net_ppt\"], through_rain\n",
+    "    )\n",
+    "    through_rain[:] = np.where(cond1 & cond2, data[\"net_rain\"], through_rain)\n",
+    "\n",
+    "    through_rain.to_dataset(name=\"through_rain\").to_netcdf(\n",
+    "        output_dir / \"through_rain.nc\"\n",
+    "    )\n",
+    "    through_rain.close()\n",
+    "\n",
+    "    # infil_hru\n",
+    "    imperv_frac = params[\"hru_percent_imperv\"]\n",
+    "    dprst_frac = params[\"dprst_frac\"]\n",
+    "    perv_frac = 1.0 - imperv_frac - dprst_frac\n",
+    "    da = xr.open_dataset(output_dir / \"infil.nc\")[\"infil\"].rename(\"infil_hru\")\n",
+    "    da *= perv_frac\n",
+    "    da.to_dataset().to_netcdf(output_dir / \"infil_hru.nc\")\n",
+    "    da.close()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "e591db15-82cc-4843-a46b-87d037823700",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# run_prms(\n",
+    "#     bin_mixed,\n",
+    "#     nb_output_dir / f\"{domain_name}_prms_mixed_run\",\n",
+    "#     skip_if_exists=skip_if_exists_prms_mixed,\n",
+    "# )"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "791f0f08-06fb-445f-b6bb-9262b9f9a6f0",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# %debug"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "155fc8ac-da82-4df7-a33a-a744818d46d7",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "prms_dbl_run_dir = nb_output_dir / f\"{domain_name}_prms_double_run\"\n",
+    "run_prms(\n",
+    "    bin_double, prms_dbl_run_dir, skip_if_exists=skip_if_exists_prms_double\n",
+    ")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "34bcdb09-c0c2-45c4-ab5f-fc22f94e77ef",
+   "metadata": {},
+   "source": [
+    "## Run pywatershed run forced with output from PRMS double precision run"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "01d0f43f-30cd-480e-9200-7c9507d4b263",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "process = [pws.PRMSRunoff]\n",
+    "\n",
+    "pws_run_dir = nb_output_dir / f\"{domain_name}_pws_run\"\n",
+    "input_dir = pws_run_dir / \"pws_input\""
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "80bc300c-2779-497f-877f-6a502531ca4a",
+   "metadata": {
+    "editable": true,
+    "slideshow": {
+     "slide_type": ""
+    }
+   },
+   "outputs": [],
+   "source": [
+    "control = pws.Control.load(domain_dir / \"control.test\")\n",
+    "output_dir = pws_run_dir / \"output\"\n",
+    "control.options = control.options | {\n",
+    "    \"input_dir\": input_dir,\n",
+    "    \"budget_type\": \"error\",\n",
+    "    \"calc_method\": \"numpy\",\n",
+    "    \"netcdf_output_dir\": output_dir,\n",
+    "}\n",
+    "params = pws.parameters.PrmsParameters.load(domain_dir / \"myparam.param\")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "4a0312cf-a286-4247-b4f0-ca8f424a09d2",
+   "metadata": {
+    "editable": true,
+    "slideshow": {
+     "slide_type": ""
+    }
+   },
+   "outputs": [],
+   "source": [
+    "if output_dir.exists() and skip_if_exists_pws:\n",
+    "    print(\n",
+    "        f\"Output ({output_dir}) already exists and skip_if_exists=True. Using existing run.\"\n",
+    "    )\n",
+    "\n",
+    "else:\n",
+    "    input_dir.mkdir(exist_ok=True, parents=True)\n",
+    "    for ff in prms_dbl_run_dir.glob(\"*.nc\"):\n",
+    "        copy2(ff, input_dir / ff.name)\n",
+    "    for ff in (prms_dbl_run_dir / \"output\").glob(\"*.nc\"):\n",
+    "        copy2(ff, input_dir / ff.name)\n",
+    "\n",
+    "    submodel = pws.Model(\n",
+    "        process,\n",
+    "        control=control,\n",
+    "        parameters=params,\n",
+    "    )\n",
+    "\n",
+    "    submodel.run(finalize=True)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "36d3958b-7c37-4e3e-9b46-cb29b9f75f0d",
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "e0e7d6ae-6502-4815-b534-1bcf9c2b948a",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "for vv in process[0].get_variables():\n",
+    "    print(vv)\n",
+    "    assert (output_dir / f\"{vv}.nc\").exists()\n",
+    "    try:\n",
+    "        assert (input_dir / f\"{vv}.nc\").exists()\n",
+    "    except:\n",
+    "        print(f\"********** {vv} not in input_dir\")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "fa4ca04f-c639-479c-8832-8cea51fc79be",
+   "metadata": {},
+   "source": [
+    "## Start by comparing the budget variables"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "0f149cac-3343-4425-bddd-9e641245f3c9",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "budget_terms = process[0].get_mass_budget_terms()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "d2cb23a3-38c3-488b-ace3-118d5fd1ae64",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# additional variables\n",
+    "budget_terms[\"outputs\"] += [\n",
+    "    \"dprst_insroff_hru\",\n",
+    "    \"dprst_stor_hru\",\n",
+    "    \"contrib_fraction\",\n",
+    "    \"infil\",\n",
+    "    \"infil_hru\",\n",
+    "    \"sroff\",\n",
+    "    \"hru_sroffp\",\n",
+    "    \"hru_sroffi\",\n",
+    "    # \"imperv_stor\",\n",
+    "    # \"imperv_evap\",\n",
+    "    \"hru_impervevap\",\n",
+    "    \"hru_impervstor\",\n",
+    "    # \"hru_impervstor_old\",\n",
+    "    \"hru_impervstor_change\",\n",
+    "    # \"dprst_vol_frac\",\n",
+    "    # \"dprst_vol_clos\",\n",
+    "    # \"dprst_vol_open\",\n",
+    "    # \"dprst_vol_clos_frac\",\n",
+    "    # \"dprst_vol_open_frac\",\n",
+    "    # \"dprst_area_clos\",\n",
+    "    # \"dprst_area_open\",\n",
+    "    # \"dprst_area_clos_max\",\n",
+    "    # \"dprst_area_open_max\",\n",
+    "    \"dprst_sroff_hru\",\n",
+    "    \"dprst_seep_hru\",\n",
+    "    \"dprst_evap_hru\",\n",
+    "    \"dprst_insroff_hru\",\n",
+    "    \"dprst_stor_hru\",\n",
+    "    # \"dprst_stor_hru_old\",\n",
+    "    \"dprst_stor_hru_change\",\n",
+    "]\n",
+    "budget_terms[\"outputs\"] = list(set(budget_terms[\"outputs\"]))"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "9b823fdf-e855-4263-b694-621c6b74005a",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "comparisons = {}\n",
+    "for term, vars in budget_terms.items():\n",
+    "    if term == \"inputs\":\n",
+    "        continue\n",
+    "    print(term)\n",
+    "    for vv in vars:\n",
+    "        print(\"    \", vv)\n",
+    "\n",
+    "        pws_file = output_dir / f\"{vv}.nc\"\n",
+    "        assert (pws_file).exists()\n",
+    "        pws_ds = xr.open_dataset(pws_file)[vv].rename(\"pws\")\n",
+    "\n",
+    "        prms_file = input_dir / f\"{vv}.nc\"\n",
+    "        assert (prms_file).exists()\n",
+    "        prms_ds = xr.open_dataset(prms_file)[vv].rename(\"prms\")\n",
+    "\n",
+    "        comparisons[vv] = xr.merge([pws_ds, prms_ds])"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "ec5e0551-9d22-4522-b80e-4c2f54eedd28",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# comparisons"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "66f80562-8679-4ef3-b007-470692de064a",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def plot_var(var_name, diff=False, nhm_id: list = None):\n",
+    "    from textwrap import fill\n",
+    "\n",
+    "    # lines = textwrap.wrap(text, width, break_long_words=False)\n",
+    "    meta = pws.meta.find_variables(var_name)[var_name]\n",
+    "    ylabel = f\"{fill(meta['desc'], 40)}\\n({meta['units']})\"\n",
+    "    title = var_name\n",
+    "    ds = comparisons[var_name]\n",
+    "\n",
+    "    if diff:\n",
+    "        ds = ds.copy()\n",
+    "        ds[\"error\"] = ds[\"pws\"] - ds[\"prms\"]\n",
+    "        ds[\"relative_error\"] = ds[\"error\"] / ds[\"prms\"]\n",
+    "        del ds[\"pws\"], ds[\"prms\"]\n",
+    "        ylabel = \"Difference PWS - PRMS\\n\" + ylabel\n",
+    "        title = \"ERRORS: Difference in \" + title\n",
+    "\n",
+    "    if (nhm_id is not None) and (len(nhm_id) > 0):\n",
+    "        ds = ds.where(ds.nhm_id.isin(nhm_id), drop=True)\n",
+    "\n",
+    "    display(\n",
+    "        ds.hvplot(\n",
+    "            frame_width=700,\n",
+    "            groupby=\"nhm_id\",\n",
+    "            title=title,\n",
+    "            ylabel=ylabel,\n",
+    "            # fontsize={\"ylabel\": \"9px\"},\n",
+    "        )\n",
+    "    )"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "38a84f40-5684-4c73-a36f-027e3dc1bdbb",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def var_close(var_name):\n",
+    "    print(var_name)\n",
+    "    var_ds = comparisons[var_name]\n",
+    "    abs_diff = abs(var_ds[\"pws\"] - var_ds[\"prms\"])\n",
+    "    rel_abs_diff = abs_diff / var_ds[\"prms\"]\n",
+    "    rtol = atol = 1.0e-7\n",
+    "    close = (abs_diff < atol) | (rel_abs_diff < rtol)\n",
+    "    if close.all():\n",
+    "        plot_var(var_name, diff=False)\n",
+    "\n",
+    "    else:\n",
+    "        wh_not_close = np.where(~close)\n",
+    "        nhm_ids = abs_diff.nhm_id[wh_not_close[1]]\n",
+    "        plot_var(var_name, diff=True, nhm_id=nhm_ids)\n",
+    "\n",
+    "    return"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "68c77573-6c17-45a1-8730-6cba7e832958",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "var_close(\"hru_impervstor_change\")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "b3b9d5ff-05e1-4468-bef6-54b3da7f5b0c",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "for var_name in comparisons.keys():\n",
+    "    var_close(var_name)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "62a7b9e8-5e95-4aa6-867c-1ad48f04b5a6",
+   "metadata": {},
+   "source": [
+    "## Look at specific time and budget errors\n",
+    "\n",
+    "> UserWarning: The flux unit balance not equal to the change in unit storageat time 1979-04-26T00:00:00 \n",
+    "> and at the following locations for PRMSRunoff: (array([341, 509]),)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "6ce47dd4-fc1d-4ba3-ba59-b69d8c692409",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "budget_terms = process[0].get_mass_budget_terms()\n",
+    "budget_terms[\"inputs\"] += [\n",
+    "    \"net_ppt\",\n",
+    "    \"net_rain\",\n",
+    "    \"net_snow\",\n",
+    "    \"pptmix_nopack\",\n",
+    "    \"pk_ice_prev\",\n",
+    "    \"freeh2o_prev\",\n",
+    "    \"newsnow\",\n",
+    "    \"snow_evap\",\n",
+    "]"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "5f52f1fb-585e-49bb-9db9-b2e87d62a973",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "budget_cases = [\n",
+    "    (\n",
+    "        \"1979-01-11T00:00:00\",\n",
+    "        [0, 1, 2, 3, 4, 5, 6, 9, 10, 11, 12, 13],\n",
+    "    ),\n",
+    "]\n",
+    "\n",
+    "case_ind = 0\n",
+    "budget_time = np.datetime64(budget_cases[case_ind][0])\n",
+    "budget_location_inds = budget_cases[case_ind][1]\n",
+    "\n",
+    "budget_comps = {}\n",
+    "for term, vars in budget_terms.items():\n",
+    "    print(term)\n",
+    "\n",
+    "    for vv in vars:\n",
+    "        print(\"    \", vv)\n",
+    "\n",
+    "        if term == \"inputs\":\n",
+    "            pws_file = input_dir / f\"{vv}.nc\"\n",
+    "        else:\n",
+    "            pws_file = output_dir / f\"{vv}.nc\"\n",
+    "\n",
+    "        assert (pws_file).exists()\n",
+    "        pws_ds = xr.open_dataset(pws_file)[vv].rename(\"pws\")\n",
+    "\n",
+    "        prms_file = input_dir / f\"{vv}.nc\"\n",
+    "        assert (prms_file).exists()\n",
+    "        prms_ds = xr.open_dataset(prms_file)[vv].rename(\"prms\")\n",
+    "\n",
+    "        budget_comps[vv] = (\n",
+    "            xr.merge([pws_ds, prms_ds])\n",
+    "            .sel(time=budget_time)\n",
+    "            .isel(nhm_id=budget_location_inds)\n",
+    "        )"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "0dad7836-6157-4982-a77e-f0f46210a27f",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "bc = budget_comps\n",
+    "inputs = bc[\"through_rain\"] + bc[\"snowmelt\"] + bc[\"intcp_changeover\"]\n",
+    "outputs = (\n",
+    "    bc[\"hru_sroffi\"]\n",
+    "    + bc[\"hru_sroffp\"]\n",
+    "    + bc[\"dprst_sroff_hru\"]\n",
+    "    + bc[\"infil_hru\"]\n",
+    "    + bc[\"hru_impervevap\"]\n",
+    "    + bc[\"dprst_seep_hru\"]\n",
+    "    + bc[\"dprst_evap_hru\"]\n",
+    ")\n",
+    "storage_changes = bc[\"hru_impervstor_change\"] + bc[\"dprst_stor_hru_change\"]\n",
+    "balance = inputs - outputs - storage_changes"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "127606d6-f28f-43c1-90f0-4f3c9bad3a5f",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "ic(budget_location_inds)\n",
+    "ic(inputs.prms.values)\n",
+    "ic(outputs.prms.values)\n",
+    "ic(storage_changes.prms.values)\n",
+    "\n",
+    "ic(\"-----------\")\n",
+    "\n",
+    "ic(bc[\"through_rain\"].pws.values)\n",
+    "\n",
+    "ic(bc[\"snow_evap\"].prms.values)\n",
+    "ic(bc[\"hru_impervstor_change\"].prms.values)\n",
+    "ic(bc[\"hru_impervstor_change\"].pws.values)\n",
+    "ic(bc[\"dprst_stor_hru_change\"].prms.values)\n",
+    "ic(bc[\"dprst_stor_hru_change\"].pws.values)\n",
+    "ic(balance.prms.values)\n",
+    "\n",
+    "# ic(bc[\"hru_sroffi\"].prms.sum().values)\n",
+    "# ic(bc[\"hru_sroffp\"].prms.sum().values)\n",
+    "# ic(bc[\"dprst_sroff_hru\"].prms.sum().values)\n",
+    "# ic(bc[\"infil_hru\"].prms.sum().values)\n",
+    "# ic(bc[\"hru_impervevap\"].prms.sum().values)\n",
+    "# ic(bc[\"dprst_seep_hru\"].prms.sum().values)\n",
+    "# ic(bc[\"dprst_evap_hru\"].prms.sum().values)\n",
+    "\n",
+    "# ic(storage_changes.prms.values)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "63297394-ec68-4645-babb-a46cd6b680bf",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "ic((balance - bc[\"through_rain\"]).pws.values)\n",
+    "ic((balance - bc[\"through_rain\"]).prms.values)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "52d6a2ed-d900-4f66-bf0b-e75c1d11fce9",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "((balance.pws - balance.prms) < 1.0e-8).all().values"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "c3b27604-a410-4e8c-ac77-a136aac554df",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "balance.pws.sum()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "3f993bb1-7590-4e47-a2a7-b0bc10350ef3",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "ic(bc[\"through_rain\"].pws.values)\n",
+    "ic(bc[\"net_rain\"].pws.values)\n",
+    "ic(bc[\"net_snow\"].pws.values)\n",
+    "ic(bc[\"net_ppt\"].pws.values)\n",
+    "ic(bc[\"pptmix_nopack\"].pws.values)\n",
+    "ic(bc[\"newsnow\"].pws.values)\n",
+    "ic((bc[\"pk_ice_prev\"].pws.values + bc[\"freeh2o_prev\"].pws.values) < epsilon32)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "968ffe3e-a984-4b36-bc9e-7f6a360b245c",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "input_max = max(\n",
+    "    abs(bc[\"through_rain\"]), abs(bc[\"snowmelt\"]), +abs(bc[\"intcp_changeover\"])\n",
+    ")\n",
+    "output_max = max(\n",
+    "    abs(bc[\"hru_sroffi\"]),\n",
+    "    abs(bc[\"hru_sroffp\"]),\n",
+    "    abs(bc[\"dprst_sroff_hru\"]),\n",
+    "    abs(bc[\"infil_hru\"]),\n",
+    "    abs(bc[\"hru_impervevap\"]),\n",
+    "    abs(bc[\"dprst_seep_hru\"]),\n",
+    "    abs(bc[\"dprst_evap_hru\"]),\n",
+    ")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "203b6241-601f-4a4d-b2d7-d0358435be9c",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "input_max.pws.values.tolist()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "7bfbe75c-a67e-4d7a-b29d-3c48d04cc53b",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "output_max.prms.values.tolist()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "9a0e4a51-6794-4aae-b379-c658f7c73bce",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "print((balance.pws / output_max.pws.max()).values)\n",
+    "print((balance.prms / output_max.prms.max()).values)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "f1e738e0-2761-4d07-8b86-a2c8a5403c5c",
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  }
+ ],
+ "metadata": {
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython"
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}
diff --git a/examples/snow_errors.ipynb b/examples/snow_errors.ipynb
new file mode 100644
index 00000000..b8243a21
--- /dev/null
+++ b/examples/snow_errors.ipynb
@@ -0,0 +1,773 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "id": "4bf60102-e977-4b2c-82b1-86780c283c71",
+   "metadata": {},
+   "source": [
+    "# Snow mass balance errors\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "6936fb48-fb5c-463e-8642-11287a578ab1",
+   "metadata": {
+    "editable": true,
+    "slideshow": {
+     "slide_type": ""
+    }
+   },
+   "outputs": [],
+   "source": [
+    "# auto-format the code in this notebook\n",
+    "%load_ext jupyter_black"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "c7b84519-933d-4284-838d-f5a8d4ec01f7",
+   "metadata": {},
+   "source": [
+    "## Setup"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "d316a5de-19ef-4008-8c68-1e192c3f89d9",
+   "metadata": {
+    "editable": true,
+    "slideshow": {
+     "slide_type": ""
+    }
+   },
+   "outputs": [],
+   "source": [
+    "import pathlib as pl\n",
+    "from pprint import pprint\n",
+    "from shutil import rmtree, copy2\n",
+    "\n",
+    "import hvplot.xarray  # noqa\n",
+    "from icecream import ic\n",
+    "from IPython.display import display\n",
+    "import numpy as np\n",
+    "import pywatershed as pws\n",
+    "import xarray as xr"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "9c188b5b-fc91-480c-a1a2-ce1121f92915",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "domain_name = \"ucb_2yr\"\n",
+    "\n",
+    "nb_output_dir = pl.Path(\"./snow_errors\")\n",
+    "\n",
+    "skip_if_exists_prms_mixed = True\n",
+    "skip_if_exists_prms_double = True\n",
+    "skip_if_exists_pws = True"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "6f807cf6-f982-4b5a-91eb-89294bf1c08e",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "pws_root = pws.constants.__pywatershed_root__\n",
+    "domain_dir = pws_root / f\"../test_data/{domain_name}\"\n",
+    "nb_output_dir.mkdir(exist_ok=True)\n",
+    "\n",
+    "zero = pws.constants.zero\n",
+    "epsilon64 = pws.constants.epsilon64\n",
+    "epsilon32 = pws.constants.epsilon32"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "0a05c33a-c86b-433a-a400-54c7f2df8eb9",
+   "metadata": {},
+   "source": [
+    "## Run PRMS mixed and double precision runs and convert to netcdf"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "daa94738-1c8b-4113-a529-c6dd0b4a0284",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "bin_dir = pws_root / \"../bin/\"\n",
+    "# bin_mixed = bin_dir / \"prms_521_mixed_mac_m1_intel\"\n",
+    "bin_double = bin_dir / \"prms_mac_m1_ifort_dbl_prec\""
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "611c0d7a-7bfd-422d-ae53-3978b0c2eb83",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def run_prms(binary: pl.Path, run_dir: pl.Path, skip_if_exists=False):\n",
+    "    import shlex\n",
+    "    import subprocess\n",
+    "\n",
+    "    from pywatershed import CsvFile, Soltab\n",
+    "\n",
+    "    from pywatershed.parameters import PrmsParameters\n",
+    "\n",
+    "    if skip_if_exists and run_dir.exists():\n",
+    "        print(\n",
+    "            f\"Run ({run_dir}) already exists and skip_if_exists=True. Using existing run.\"\n",
+    "        )\n",
+    "        return None\n",
+    "\n",
+    "    run_dir.mkdir()  # must not exist, on user to delete\n",
+    "    copy2(binary, run_dir / binary.name)\n",
+    "    for ff in [\n",
+    "        \"control.test\",\n",
+    "        \"myparam.param\",\n",
+    "        \"tmax.cbh\",\n",
+    "        \"tmin.cbh\",\n",
+    "        \"prcp.cbh\",\n",
+    "        \"sf_data\",\n",
+    "    ]:\n",
+    "        copy2(domain_dir / ff, run_dir / ff)\n",
+    "\n",
+    "    output_dir = run_dir / \"output\"\n",
+    "    output_dir.mkdir()\n",
+    "\n",
+    "    exe_command = f\"time ./{binary.name} control.test -MAXDATALNLEN 60000 2>&1 | tee run.log\"\n",
+    "    result = subprocess.run(\n",
+    "        exe_command,\n",
+    "        shell=True,\n",
+    "        # stdout = subprocess.PIPE,\n",
+    "        stderr=subprocess.STDOUT,\n",
+    "        universal_newlines=True,\n",
+    "        cwd=run_dir,\n",
+    "    )\n",
+    "\n",
+    "    # these will be useful in what follows\n",
+    "    params = pws.parameters.PrmsParameters.load(\n",
+    "        domain_dir / \"myparam.param\"\n",
+    "    ).parameters\n",
+    "\n",
+    "    # convert to netcdf\n",
+    "    # could make these arguments\n",
+    "    chunking = {\n",
+    "        \"time\": 0,\n",
+    "        \"doy\": 0,\n",
+    "        \"nhm_id\": 100,\n",
+    "        \"nhm_seg\": 100,\n",
+    "    }\n",
+    "\n",
+    "    output_csvs = output_dir.glob(\"*.csv\")\n",
+    "    for cc in output_csvs:\n",
+    "        if cc.name in [\"stats.csv\"]:\n",
+    "            continue\n",
+    "        nc_path = cc.with_suffix(\".nc\")\n",
+    "        CsvFile(cc).to_netcdf(nc_path, chunk_sizes=chunking)\n",
+    "\n",
+    "    # solar tables\n",
+    "    soltab_file = run_dir / \"soltab_debug\"\n",
+    "    # the nhm_ids are not available in the solta_debug file currently, so get\n",
+    "    # them from the domain parameters\n",
+    "    params = PrmsParameters.load(run_dir / \"myparam.param\")\n",
+    "    nhm_ids = params.parameters[\"nhm_id\"]\n",
+    "\n",
+    "    soltab = Soltab(\n",
+    "        soltab_file,\n",
+    "        output_dir=output_dir,\n",
+    "        nhm_ids=nhm_ids,\n",
+    "        chunk_sizes=chunking,\n",
+    "    )\n",
+    "\n",
+    "    for var in soltab.variables:\n",
+    "        assert (output_dir / f\"{var}.nc\").exists()\n",
+    "\n",
+    "    # previous and change variables\n",
+    "    for vv in [\"pk_ice\", \"freeh2o\", \"soil_moist\"]:\n",
+    "        data = xr.open_dataset(output_dir / f\"{vv}.nc\")[vv]\n",
+    "        prev_da = data.copy()\n",
+    "        prev_da[:] = np.roll(prev_da.values, 1, axis=0)\n",
+    "        assert (prev_da[1:, :].values == data[0:-1, :].values).all()\n",
+    "        prev_da[0, :] = np.zeros(1)[\n",
+    "            0\n",
+    "        ]  # np.nan better but causes plotting to fail\n",
+    "        change_da = data - prev_da\n",
+    "        prev_da.rename(f\"{vv}_prev\").to_dataset().to_netcdf(\n",
+    "            output_dir / f\"{vv}_prev.nc\"\n",
+    "        )\n",
+    "        data[f\"{vv}_prev\"] = xr.open_dataset(output_dir / f\"{vv}_prev.nc\")[\n",
+    "            f\"{vv}_prev\"\n",
+    "        ]\n",
+    "\n",
+    "        change_da.rename(f\"{vv}_change\").to_dataset().to_netcdf(\n",
+    "            output_dir / f\"{vv}_change.nc\"\n",
+    "        )\n",
+    "        data[f\"{vv}_change\"] = xr.open_dataset(output_dir / f\"{vv}_change.nc\")[\n",
+    "            f\"{vv}_change\"\n",
+    "        ]\n",
+    "    # through_rain\n",
+    "    dep_vars = [\n",
+    "        \"net_ppt\",\n",
+    "        \"pptmix_nopack\",\n",
+    "        \"snowmelt\",\n",
+    "        \"pkwater_equiv\",\n",
+    "        \"snow_evap\",\n",
+    "        \"net_snow\",\n",
+    "        \"net_rain\",\n",
+    "    ]\n",
+    "    data = {}\n",
+    "    for vv in dep_vars:\n",
+    "        data[vv] = xr.open_dataset(output_dir / f\"{vv}.nc\")[vv]\n",
+    "\n",
+    "    nearzero = 1.0e-6\n",
+    "\n",
+    "    cond1 = data[\"net_ppt\"] > zero\n",
+    "    cond2 = data[\"pptmix_nopack\"] != 0\n",
+    "    cond3 = data[\"snowmelt\"] < nearzero\n",
+    "    cond4 = data[\"pkwater_equiv\"] < epsilon32\n",
+    "    cond5 = data[\"snow_evap\"] < nearzero\n",
+    "    cond6 = data[\"net_snow\"] < nearzero\n",
+    "\n",
+    "    through_rain = data[\"net_rain\"] * zero\n",
+    "    # these are in reverse order\n",
+    "    through_rain[:] = np.where(\n",
+    "        cond1 & cond3 & cond4 & cond6, data[\"net_rain\"], zero\n",
+    "    )\n",
+    "    through_rain[:] = np.where(\n",
+    "        cond1 & cond3 & cond4 & cond5, data[\"net_ppt\"], through_rain\n",
+    "    )\n",
+    "    through_rain[:] = np.where(cond1 & cond2, data[\"net_rain\"], through_rain)\n",
+    "\n",
+    "    through_rain.to_dataset(name=\"through_rain\").to_netcdf(\n",
+    "        output_dir / \"through_rain.nc\"\n",
+    "    )\n",
+    "    through_rain.close()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "e591db15-82cc-4843-a46b-87d037823700",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# run_prms(bin_mixed, nb_output_dir / \"prms_mixed_run\")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "155fc8ac-da82-4df7-a33a-a744818d46d7",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "prms_dbl_run_dir = nb_output_dir / f\"{domain_name}_prms_double_run\"\n",
+    "run_prms(\n",
+    "    bin_double, prms_dbl_run_dir, skip_if_exists=skip_if_exists_prms_double\n",
+    ")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "34bcdb09-c0c2-45c4-ab5f-fc22f94e77ef",
+   "metadata": {},
+   "source": [
+    "## pywatershed run"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "01d0f43f-30cd-480e-9200-7c9507d4b263",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "process = [pws.PRMSSnow]\n",
+    "pws_run_dir = nb_output_dir / f\"{domain_name}_pws_run\"\n",
+    "input_dir = pws_run_dir / \"pws_input\""
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "80bc300c-2779-497f-877f-6a502531ca4a",
+   "metadata": {
+    "editable": true,
+    "slideshow": {
+     "slide_type": ""
+    }
+   },
+   "outputs": [],
+   "source": [
+    "control = pws.Control.load(domain_dir / \"control.test\")\n",
+    "output_dir = pws_run_dir / \"output\"\n",
+    "control.options = control.options | {\n",
+    "    \"input_dir\": input_dir,\n",
+    "    \"budget_type\": \"warn\",\n",
+    "    \"calc_method\": \"numpy\",\n",
+    "    \"netcdf_output_dir\": output_dir,\n",
+    "}\n",
+    "params = pws.parameters.PrmsParameters.load(domain_dir / \"myparam.param\")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "65276857-f940-46df-8385-a9791bc79894",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "if output_dir.exists() and skip_if_exists_pws:\n",
+    "    print(\n",
+    "        f\"Output ({output_dir}) already exists and skip_if_exists=True. Using existing run.\"\n",
+    "    )\n",
+    "\n",
+    "else:\n",
+    "    input_dir.mkdir(exist_ok=True, parents=True)\n",
+    "    for ff in prms_dbl_run_dir.glob(\"*.nc\"):\n",
+    "        copy2(ff, input_dir / ff.name)\n",
+    "    for ff in (prms_dbl_run_dir / \"output\").glob(\"*.nc\"):\n",
+    "        copy2(ff, input_dir / ff.name)\n",
+    "\n",
+    "    submodel = pws.Model(\n",
+    "        process,\n",
+    "        control=control,\n",
+    "        parameters=params,\n",
+    "    )\n",
+    "\n",
+    "    submodel.run(finalize=True)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "9a056b6a-d31a-492b-a761-3a33f42ba08c",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "%debug"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "e0e7d6ae-6502-4815-b534-1bcf9c2b948a",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "for vv in process[0].get_variables():\n",
+    "    print(vv)\n",
+    "    assert (output_dir / f\"{vv}.nc\").exists()\n",
+    "    try:\n",
+    "        assert (input_dir / f\"{vv}.nc\").exists()\n",
+    "    except:\n",
+    "        print(f\"********** {vv} not in input_dir\")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "fa4ca04f-c639-479c-8832-8cea51fc79be",
+   "metadata": {},
+   "source": [
+    "## Start by comparing the budget variables"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "0f149cac-3343-4425-bddd-9e641245f3c9",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "budget_terms = pws.PRMSSnow.get_mass_budget_terms()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "d2cb23a3-38c3-488b-ace3-118d5fd1ae64",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# additional variables\n",
+    "budget_terms[\"outputs\"] += [\n",
+    "    \"pk_ice_prev\",\n",
+    "    \"freeh2o_prev\",\n",
+    "    \"newsnow\",\n",
+    "    \"pptmix_nopack\",\n",
+    "    # \"ai\",\n",
+    "    \"albedo\",\n",
+    "    # 'frac_swe',\n",
+    "    \"freeh2o\",\n",
+    "    \"freeh2o_change\",\n",
+    "    \"freeh2o_prev\",\n",
+    "    #' iasw',\n",
+    "    # 'int_alb',\n",
+    "    \"iso\",\n",
+    "    # 'lso',\n",
+    "    # 'lst',\n",
+    "    # \"mso\",\n",
+    "    \"newsnow\",\n",
+    "    \"pk_def\",\n",
+    "    \"pk_den\",\n",
+    "    \"pk_depth\",\n",
+    "    \"pk_ice\",\n",
+    "    \"pk_ice_change\",\n",
+    "    \"pk_ice_prev\",\n",
+    "    #  'pk_precip',\n",
+    "    \"pk_temp\",\n",
+    "    # 'pksv',\n",
+    "    \"pkwater_ante\",\n",
+    "    \"pkwater_equiv\",\n",
+    "    # 'pkwater_equiv_change',\n",
+    "    \"pptmix_nopack\",\n",
+    "    # 'pss',\n",
+    "    \"pst\",\n",
+    "    # \"salb\",\n",
+    "    #' scrv',\n",
+    "    #' slst',\n",
+    "    \"snow_evap\",\n",
+    "    \"snowcov_area\",\n",
+    "    # 'snowcov_areasv',\n",
+    "    \"snowmelt\",\n",
+    "    # 'snsv',\n",
+    "    \"tcal\",\n",
+    "    \"through_rain\",\n",
+    "]"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "9b823fdf-e855-4263-b694-621c6b74005a",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "comparisons = {}\n",
+    "for term, vars in budget_terms.items():\n",
+    "    if term == \"inputs\":\n",
+    "        continue\n",
+    "    print(term)\n",
+    "    for vv in sorted(vars):\n",
+    "        print(\"    \", vv)\n",
+    "\n",
+    "        pws_file = output_dir / f\"{vv}.nc\"\n",
+    "        assert (pws_file).exists()\n",
+    "        pws_ds = xr.open_dataset(pws_file)[vv].rename(\"pws\")\n",
+    "\n",
+    "        prms_file = input_dir / f\"{vv}.nc\"\n",
+    "        assert prms_file.exists()\n",
+    "        prms_ds = xr.open_dataset(prms_file)[vv].rename(\"prms\")\n",
+    "\n",
+    "        comparisons[vv] = xr.merge([pws_ds, prms_ds])"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "ec5e0551-9d22-4522-b80e-4c2f54eedd28",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# comparisons"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "66f80562-8679-4ef3-b007-470692de064a",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def plot_var(var_name, diff=False, nhm_id: list = None):\n",
+    "    from textwrap import fill\n",
+    "\n",
+    "    # lines = textwrap.wrap(text, width, break_long_words=False)\n",
+    "    meta = pws.meta.find_variables(var_name)[var_name]\n",
+    "    ylabel = f\"{fill(meta['desc'], 40)}\\n({meta['units']})\"\n",
+    "    title = var_name\n",
+    "    ds = comparisons[var_name]\n",
+    "\n",
+    "    if diff:\n",
+    "        ds = ds.copy()\n",
+    "        ds[\"error\"] = ds[\"pws\"] - ds[\"prms\"]\n",
+    "        # ds[\"relative_error\"] = ds[\"error\"] / ds[\"prms\"]\n",
+    "        # ds[\"relative_error\"] = xr.where(\n",
+    "        #     abs(ds[\"prms\"]) < 1.0e-7, np.nan, ds[\"relative_error\"]\n",
+    "        # )\n",
+    "        del ds[\"pws\"], ds[\"prms\"]\n",
+    "        ylabel = \"Difference PWS - PRMS\\n\" + ylabel\n",
+    "        title = \"ERRORS: Difference in \" + title\n",
+    "\n",
+    "    if (nhm_id is not None) and (len(nhm_id) > 0):\n",
+    "        ds = ds.where(ds.nhm_id.isin(nhm_id), drop=True)\n",
+    "\n",
+    "    display(\n",
+    "        ds.hvplot(\n",
+    "            frame_width=700,\n",
+    "            groupby=\"nhm_id\",\n",
+    "            title=title,\n",
+    "            ylabel=ylabel,\n",
+    "            # fontsize={\"ylabel\": \"9px\"},\n",
+    "        )\n",
+    "    )"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "38a84f40-5684-4c73-a36f-027e3dc1bdbb",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def var_close(var_name):\n",
+    "    print(var_name)\n",
+    "    var_ds = comparisons[var_name]\n",
+    "    abs_diff = abs(var_ds[\"pws\"] - var_ds[\"prms\"])\n",
+    "    rel_abs_diff = abs_diff / var_ds[\"prms\"]\n",
+    "    rtol = 0.02\n",
+    "    atol = 1.0e-2\n",
+    "    close = (abs_diff < atol) | (rel_abs_diff < rtol)\n",
+    "    if close.all():\n",
+    "        plot_var(var_name, diff=False)\n",
+    "\n",
+    "    else:\n",
+    "        wh_not_close = np.where(~close)\n",
+    "        nhm_ids = abs_diff.nhm_id[wh_not_close[1]]\n",
+    "        plot_var(var_name, diff=False, nhm_id=nhm_ids)\n",
+    "\n",
+    "    return"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "b3b9d5ff-05e1-4468-bef6-54b3da7f5b0c",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "if False:\n",
+    "    for var_name in comparisons.keys():\n",
+    "        var_close(var_name)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "da70848a-3187-4ae7-b045-125c83e451c9",
+   "metadata": {},
+   "source": [
+    "## Mass-balance errors\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "b5d5e323-a739-43a6-a372-5bcb7dc7b775",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "budget_terms = pws.PRMSSnow.get_mass_budget_terms()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "ff49bb08-0ae4-41e9-9e52-adeb0d48cbbe",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "budget_terms"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "bacee4a6-0952-4c96-aa50-b05f5664e74b",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "prms_mass = {}\n",
+    "for term, vars in budget_terms.items():\n",
+    "    prms_mass[term] = {}\n",
+    "    for vv in vars:\n",
+    "        prms_mass[term][vv] = xr.open_dataarray(input_dir / f\"{vv}.nc\")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "3deaaa12-3094-4283-8edf-1c1ad0acdfe7",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "prms_mass_sums = {}\n",
+    "for term, vars in prms_mass.items():\n",
+    "    prms_mass_sums[term] = None\n",
+    "    for vv in vars:\n",
+    "        data = xr.open_dataarray(input_dir / f\"{vv}.nc\")\n",
+    "        if prms_mass_sums[term] is None:\n",
+    "            prms_mass_sums[term] = data\n",
+    "        else:\n",
+    "            prms_mass_sums[term] += xr.open_dataarray(input_dir / f\"{vv}.nc\")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "552ff261-c959-4d89-8f7f-53160c928f75",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "prms_mass_balance = (\n",
+    "    prms_mass_sums[\"inputs\"]\n",
+    "    - prms_mass_sums[\"outputs\"]\n",
+    "    - prms_mass_sums[\"storage_changes\"]\n",
+    ")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "bf5228fd-0a6c-4c9e-bc29-adf8373743f2",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "prms_mass_balance"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "bbd6c480-fc75-4546-9246-2c5f5b2aee46",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "wh_imbalance = np.where(prms_mass_balance.values > 1e-5)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "2be69429-2ab7-4307-90c6-dedb4262635a",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "wh_imbalance"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "7bb8a452-ae77-48e9-9120-02a8c1ebcd75",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "prms_mass_balance.values[wh_imbalance]"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "0c979a40-5740-4971-a0aa-d5e0f5dd88b7",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "freeh2o = xr.open_dataarray(input_dir / \"freeh2o.nc\")\n",
+    "pk_ice = xr.open_dataarray(input_dir / \"pk_ice.nc\")\n",
+    "pkwater_equiv = xr.open_dataarray(input_dir / \"pkwater_equiv.nc\")\n",
+    "freeh2o_prev = xr.open_dataarray(input_dir / \"freeh2o_prev.nc\")\n",
+    "pk_ice_prev = xr.open_dataarray(input_dir / \"pk_ice_prev.nc\")\n",
+    "pkwater_ante = xr.open_dataarray(input_dir / \"pkwater_ante.nc\")\n",
+    "freeh2o_change = xr.open_dataarray(input_dir / \"freeh2o_change.nc\")\n",
+    "pk_ice_change = xr.open_dataarray(input_dir / \"pk_ice_change.nc\")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "be45ba5c-75ce-49ac-89d1-8f91aaaadec3",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "ic(pk_ice.values[wh_imbalance])\n",
+    "ic(pk_ice_prev.values[wh_imbalance])\n",
+    "ic(pk_ice_change.values[wh_imbalance])"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "3220c630-52c8-4074-8794-54a9f462a7a3",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "ic(freeh2o.values[wh_imbalance])\n",
+    "ic(freeh2o_prev.values[wh_imbalance])\n",
+    "ic(freeh2o_change.values[wh_imbalance])"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "d7ae9821-c6b6-422b-89d8-04f1037c83eb",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import numpy as np\n",
+    "\n",
+    "np.set_printoptions(edgeitems=30, linewidth=100000)\n",
+    "\n",
+    "for term, vars in prms_mass.items():\n",
+    "    print(term)\n",
+    "    for vv in vars:\n",
+    "        print(f\"    {vv}\")\n",
+    "        print(f\"    {prms_mass[term][vv].values[wh_imbalance]}\")\n",
+    "print(\"imbalance\")\n",
+    "print(f\"    {prms_mass_balance.values[wh_imbalance]}\")\n",
+    "print(\"freeh2o prev\")\n",
+    "print(f\"    {freeh2o_prev.values[wh_imbalance]}\")\n",
+    "print(\"freeh2o\")\n",
+    "print(f\"    {freeh2o.values[wh_imbalance]}\")\n",
+    "print(\"pk_ice prev\")\n",
+    "print(f\"    {pk_ice_prev.values[wh_imbalance]}\")\n",
+    "print(\"pk_ice\")\n",
+    "print(f\"    {pk_ice.values[wh_imbalance]}\")\n",
+    "print(\"pkwater_ante\")\n",
+    "print(f\"    {pkwater_ante.values[wh_imbalance]}\")\n",
+    "print(\"pkwater_equiv\")\n",
+    "print(f\"    {pkwater_equiv.values[wh_imbalance]}\")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "a882b844-6fa3-4005-8660-36d45deb51f6",
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  }
+ ],
+ "metadata": {
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython"
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}
diff --git a/prms_src/prms5.2.1/makelist b/prms_src/prms5.2.1/makelist
index 6f9ffa48..b4b0e1e6 100644
--- a/prms_src/prms5.2.1/makelist
+++ b/prms_src/prms5.2.1/makelist
@@ -37,7 +37,7 @@ endif
 
 $(info    ----------------)
 $(info    FC is $(FC))
-$(info    FC is $(CC))
+$(info    CC is $(CC))
 $(info    DBL_PREC is false)
 $(info    ----------------)
 $(info    )
diff --git a/prms_src/prms5.2.1/makelist_double_precision b/prms_src/prms5.2.1/makelist_double_precision
index b45d4f4f..223f23cc 100644
--- a/prms_src/prms5.2.1/makelist_double_precision
+++ b/prms_src/prms5.2.1/makelist_double_precision
@@ -37,7 +37,7 @@ endif
 
 $(info    ----------------)
 $(info    FC is $(FC))
-$(info    FC is $(CC))
+$(info    CC is $(CC))
 $(info    DBL_PREC is true)
 $(info    ----------------)
 $(info    )
@@ -78,7 +78,7 @@ endif
 # This is a shortcut to alleviate errors when parsing parameter inputs in soltab.
 ifeq ($(detected_OS), Windows)
 	ifeq ($(FC), gfortran)
-		FFLAGS = $(OPTLEVEL) -fno-second-underscore -fallow-argument-mismatch
+		FFLAGS = -freal-4-real-8 $(OPTLEVEL) -fno-second-underscore -fallow-argument-mismatch
 	endif
 else
 	ifeq ($(FC), gfortran)
@@ -99,7 +99,7 @@ endif
 
 ifeq ($(detected_OS), Windows)
 	ifeq ($(CC), gcc)
-		CFLAGS = $(OPTLEVEL) -D$(ARC) -D_UF -Wall
+		CFLAGS = -Dfloat=double $(OPTLEVEL) -D$(ARC) -D_UF -Wall
 	endif
 else
 	ifeq ($(CC), gcc)
diff --git a/prms_src/prms5.2.1/prms/snowcomp.f90 b/prms_src/prms5.2.1/prms/snowcomp.f90
index fad3cf19..ef2bd151 100644
--- a/prms_src/prms5.2.1/prms/snowcomp.f90
+++ b/prms_src/prms5.2.1/prms/snowcomp.f90
@@ -1435,7 +1435,8 @@ INTEGER FUNCTION snorun()
             ENDIF
           ENDIF
 
-        ENDIF
+        ENDIF  ! ends iff statement between steps 3 and 4 about 265 lines above
+
 
 ! LAST check to clear out all arrays if packwater is gone
         IF ( Pkwater_equiv(i)<=0.0D0 ) THEN
@@ -2598,6 +2599,9 @@ SUBROUTINE snowevap(Potet_sublim, Potet, Snowcov_area, Snow_evap, &
 !RAPCOMMENT - CHANGED TO CHECK FOR NEGATIVE PACK ICE
           ! If all pack ice is removed, then there cannot be a
           ! heat deficit
+          ! JLM mass balance fix, if we take pkwater_equiv - ez below
+          ! JLM we need to take ez from both pk_ice and freeh2o here
+          freeh2o = freeh2o + pk_ice
           Pk_ice = 0.0
           Pk_def = 0.0
           Pk_temp = 0.0
@@ -2613,17 +2617,22 @@ SUBROUTINE snowevap(Potet_sublim, Potet, Snowcov_area, Snow_evap, &
         Pkwater_equiv = Pkwater_equiv - ez
         Snow_evap = ez
       ENDIF
+
       IF ( Snow_evap<0.0 ) THEN
         Pkwater_equiv = Pkwater_equiv - DBLE(Snow_evap)
+
         IF ( Pkwater_equiv<0.0D0 ) THEN
           IF ( Print_debug>DEBUG_less ) THEN
             IF ( Pkwater_equiv<-DNEARZERO ) &
      &           PRINT *, 'snowpack issue, negative pkwater_equiv in snowevap', Pkwater_equiv
-            Pkwater_equiv = 0.0D0
+            Pkwater_equiv = 0.0D0  ! JLM: this is INSIDE a debug statement? will change the answers
+            !  is this in the originial source
           ENDIF
         ENDIF
+
         Snow_evap = 0.0
       ENDIF
+
       avail_et = Potet - Hru_intcpevap - Snow_evap
       IF ( avail_et<0.0 ) THEN
 !        PRINT *, 'snow evap', snow_evap, avail_et, pkwater_equiv
diff --git a/pywatershed/base/budget.py b/pywatershed/base/budget.py
index a36e99d7..d8a428e2 100644
--- a/pywatershed/base/budget.py
+++ b/pywatershed/base/budget.py
@@ -8,7 +8,7 @@
 
 from pywatershed.base.control import Control
 
-from ..constants import zero
+from ..constants import epsilon64, one, zero
 from ..utils.formatting import pretty_print
 from ..utils.netcdf_utils import NetCdfWrite
 from .accessor import Accessor
@@ -298,18 +298,27 @@ def _sum_storage_changes(self):
         return self._sum("storage_changes")
 
     def _calc_unit_balance(self):
-        unit_balance = self._inputs_sum - self._outputs_sum
         self._zero_sum = True
-        if not np.allclose(
-            unit_balance,
-            self._storage_changes_sum,
-            rtol=self.rtol,
-            atol=self.atol,
-        ):
+
+        # roll our own np.allclose so we can diagnose the not close points
+        unit_balance = self._inputs_sum - self._outputs_sum
+        abs_diff = np.abs(unit_balance - self._storage_changes_sum)
+        mask_div_zero = self._storage_changes_sum < epsilon64
+        rel_abs_diff = abs_diff / np.where(
+            mask_div_zero, one, self._storage_changes_sum
+        )
+        abs_close = abs_diff < self.atol
+        rel_close = rel_abs_diff < self.rtol
+        either_close = abs_close | rel_close
+        cond = np.where(mask_div_zero, abs_close, either_close)
+
+        if not cond.all():
             self._zero_sum = False
+            wh_not_cond = np.where(~cond)
             msg = (
-                "The flux unit balance not equal to the change in unit storage"
-                f": {self.description}"
+                "The flux unit balance not equal to the change in unit "
+                f"storage at time {self.control.current_time} and at the "
+                f"following locations for {self.description}: {wh_not_cond}"
             )
             if self.imbalance_fatal:
                 raise ValueError(msg)
diff --git a/pywatershed/base/control.py b/pywatershed/base/control.py
index d99f80b4..3f3730cf 100644
--- a/pywatershed/base/control.py
+++ b/pywatershed/base/control.py
@@ -108,19 +108,19 @@ def current_year(self):
         return datetime_year(self._current_time)
 
     @property
-    def current_month(self):
-        """Get the current month."""
-        return datetime_month(self._current_time)
+    def current_month(self, zero_based: bool = False):
+        """Get the current month in 1-12 (unless zero based)."""
+        return datetime_month(self._current_time, zero_based=zero_based)
 
     @property
-    def current_doy(self):
-        """Get the current day of year."""
-        return datetime_doy(self._current_time)
+    def current_doy(self, zero_based: bool = False):
+        """Get the current day of year in 1-366 (unless zero based)."""
+        return datetime_doy(self._current_time, zero_based=zero_based)
 
     @property
-    def current_dowy(self):
-        """Get the current day of water year."""
-        return datetime_dowy(self._current_time)
+    def current_dowy(self, zero_based: bool = False):
+        """Get the current day of water year in 1-366 (unless zero-based)."""
+        return datetime_dowy(self._current_time, zero_based=zero_based)
 
     @property
     def current_epiweek(self):
diff --git a/pywatershed/hydrology/prms_runoff.py b/pywatershed/hydrology/prms_runoff.py
index a26f9f06..85304316 100644
--- a/pywatershed/hydrology/prms_runoff.py
+++ b/pywatershed/hydrology/prms_runoff.py
@@ -7,7 +7,14 @@
 from ..base.adapter import adaptable
 from ..base.conservative_process import ConservativeProcess
 from ..base.control import Control
-from ..constants import HruType, nan, numba_num_threads, zero
+from ..constants import (
+    HruType,
+    epsilon32,
+    epsilon64,
+    nan,
+    numba_num_threads,
+    zero,
+)
 from ..parameters import Parameters
 
 NEARZERO = 1.0e-6
@@ -25,6 +32,8 @@
 LAND = HruType.LAND.value
 LAKE = HruType.LAKE.value
 
+# TODO: using through_rain and not net_rain and net_ppt is a WIP
+
 
 class PRMSRunoff(ConservativeProcess):
     """PRMS surface runoff.
@@ -188,7 +197,7 @@ def get_init_values() -> dict:
             "hru_impervstor_change": zero,
             "dprst_vol_frac": zero,
             "dprst_vol_clos": zero,
-            "dprst_vol_open": zero,
+            "dprst_vol_open": nan,
             "dprst_vol_clos_frac": zero,
             "dprst_vol_open_frac": zero,
             "dprst_area_clos": zero,
@@ -209,6 +218,7 @@ def get_mass_budget_terms():
         return {
             "inputs": [
                 # "net_rain",
+                # "net_snow",
                 "through_rain",
                 "snowmelt",
                 "intcp_changeover",
@@ -500,6 +510,7 @@ def _calculate(self, time_length, vectorized=False):
             compute_infil=self.compute_infil,
             dprst_comp=self.dprst_comp,
             imperv_et=self.imperv_et,
+            through_rain=self.through_rain,
         )
 
         self.infil_hru[:] = self.infil * self.hru_frac_perv
@@ -586,10 +597,8 @@ def _calculate_numpy(
         compute_infil,
         dprst_comp,
         imperv_et,
+        through_rain,
     ):
-        # move towards replacing net_rain and other variables.
-        # self.net_rain[:] = self.through_rain
-
         dprst_chk = 0
         infil[:] = 0.0
         for k in prange(nhru):
@@ -615,8 +624,7 @@ def _calculate_numpy(
                 hru_impervevap[i] = 0.0
 
             avail_et = potet[i] - snow_evap[i] - hru_intcpevap[i]
-            # apparently not used
-            # availh2o = intcp_changeover[i] + net_rain[i]
+            availh2o = intcp_changeover[i] + net_rain[i]
 
             (
                 sri,
@@ -648,6 +656,7 @@ def _calculate_numpy(
                 srp=srp,
                 check_capacity=check_capacity,
                 perv_comp=perv_comp,
+                through_rain=through_rain[i],
             )
 
             dprst_flag = ACTIVE  # cdl todo: hardwired
@@ -709,7 +718,7 @@ def _calculate_numpy(
                             dprst_flow_coef=dprst_flow_coef[i],
                             dprst_seep_rate_clos=dprst_seep_rate_clos[i],
                             avail_et=avail_et,
-                            net_rain=net_rain[i],
+                            net_rain=availh2o,
                             dprst_in=dprst_in[i],
                             srp=srp,
                             sri=sri,
@@ -816,12 +825,14 @@ def compute_infil(
         srp,
         check_capacity,
         perv_comp,
+        through_rain,
     ):
         isglacier = False  # todo -- hardwired
-        cascade_active = False
         hru_flag = 0
         if hru_type == LAND or isglacier:
             hru_flag = 1
+
+        cascade_active = False
         if cascade_active:
             raise Exception("bad bad bad")
         else:
@@ -846,17 +857,21 @@ def compute_infil(
         # if rain/snow event with no antecedent snowpack,
         # compute the runoff from the rain first and then proceed with the
         # snowmelt computations
-        if pptmix_nopack == ACTIVE:
-            avail_water = avail_water + net_rain
-            infil = infil + net_rain
+        # double_counting = 0
+        cond2 = pptmix_nopack != 0
+        # if pptmix_nopack == ACTIVE:
+        if cond2:
+            # double_counting += 1
+            avail_water = avail_water + through_rain
+            infil = infil + through_rain
             if hru_flag == 1:
                 infil, srp, contrib_fraction = perv_comp(
                     soil_moist_prev,
                     carea_max,
                     smidx_coef,
                     smidx_exp,
-                    net_rain,
-                    net_rain,
+                    through_rain,
+                    through_rain,
                     infil,
                     srp,
                 )
@@ -866,13 +881,19 @@ def compute_infil(
         # procedure is used.  If there is no snowpack and no precip,
         # then check for melt from last of snowpack.  If rain/snow mix
         # with no antecedent snowpack, compute snowmelt portion of runoff.
+        cond3 = snowmelt < NEARZERO
+        cond4 = pkwater_equiv < epsilon64
+        cond1 = net_ppt > zero
+        cond6 = net_snow < NEARZERO
+
         if snowmelt > 0.0:
+            # if not cond3:  # this results in less precision accuracy
             # There was no snowmelt but a snowpack may exist.  If there is
             # no snowpack then check for rain on a snowfree HRU.
             avail_water = avail_water + snowmelt
             infil = infil + snowmelt
             if hru_flag == 1:
-                if pkwater_equiv > 0.0 or net_ppt - net_snow < NEARZERO:
+                if (pkwater_equiv > 0.0) or (net_rain < NEARZERO):
                     # Pervious area computations
                     infil, srp = check_capacity(
                         soil_moist_prev,
@@ -883,6 +904,11 @@ def compute_infil(
                     )
                 else:
                     # Snowmelt occurred and depleted the snowpack
+                    # this frequently gets counted along with pptmix_nopack
+                    # double_counting += 1
+                    # if double_counting > 1:
+                    #     print("snowmelt")
+
                     infil, srp, contrib_fraction = perv_comp(
                         soil_moist_prev,
                         carea_max,
@@ -894,20 +920,27 @@ def compute_infil(
                         srp,
                     )
 
-        elif pkwater_equiv < DNEARZERO:
+        # elif pkwater_equiv < DNEARZERO:
+        elif cond4:
             # If no snowmelt and no snowpack but there was net snow then
             # snowpack was small and was lost to sublimation.
-            if net_snow < NEARZERO and net_rain > 0.0:
-                avail_water = avail_water + net_rain
-                infil = infil + net_rain
+            # if net_snow < NEARZERO and net_rain > 0.0:
+            if cond6 and through_rain > 0.0:
+                # cond3 & cond4 & cond6 & cond1
+                # this is through_rain's top/most narrow case
+                avail_water = avail_water + through_rain
+                infil = infil + through_rain
+                # double_counting += 1
+                # if double_counting > 1:
+                #     print("cond4")
                 if hru_flag == 1:
                     infil, srp, contrib_fraction = perv_comp(
                         soil_moist_prev,
                         carea_max,
                         smidx_coef,
                         smidx_exp,
-                        net_rain,
-                        net_rain,
+                        through_rain,
+                        through_rain,
                         infil,
                         srp,
                     )
@@ -925,6 +958,9 @@ def compute_infil(
                     srp,
                 )
 
+        # if double_counting > 1:
+        #     assert False
+
         if hruarea_imperv > 0.0:
             imperv_stor = imperv_stor + avail_water
             if hru_flag == 1:
@@ -1048,7 +1084,9 @@ def dprst_comp(
             if sri < 0.0:
                 if sri < -NEARZERO:
                     sri = 0.0
-            dprst_insroff_hru = dprst_srp + dprst_sri
+
+        # <<<
+        dprst_insroff_hru = dprst_srp + dprst_sri
 
         dprst_area_open = 0.0
         if dprst_vol_open > 0.0:
diff --git a/pywatershed/hydrology/prms_snow.py b/pywatershed/hydrology/prms_snow.py
index 86694cde..f46f970d 100644
--- a/pywatershed/hydrology/prms_snow.py
+++ b/pywatershed/hydrology/prms_snow.py
@@ -64,6 +64,7 @@
 
 maxalb = 15
 ONETHIRD = 1.0 / 3.0
+tiny_snowpack = 1.0e-12
 
 tcind = 0
 
@@ -697,22 +698,24 @@ def _calculate_numpy(
 
         # cals = zero  # JLM this is unnecessary.
 
-        # newsnow
-        newsnow[:] = False
-        net_snow_gt_zero = net_snow > zero
-        wh_net_snow_gt_zero = np.where(net_snow_gt_zero)
-        newsnow[wh_net_snow_gt_zero] = True
+        # newsnow is a doganostic for prms_snow, so it lives here
+        newsnow = np.where(net_snow > zero, True, False)
 
-        # default assumption
-        pptmix_nopack[:] = False
+        # JLM TODO: there's a conditional here we dont have
+        #  in fotran trd is scalar and the RHS terms are vector?
+        trd = orad_hru / soltab_horad_potsw
+
+        frac_swe[:] = zero
+        pk_precip[:] = zero  # [inches]
+        snowmelt[:] = zero  # [inches]
+        snow_evap[:] = zero  # [inches]
+        tcal[:] = zero
+        ai[:] = zero
 
         for jj in prange(nhru):
             if hru_type[jj] == HruType.LAKE.value:
                 continue
 
-            # <
-            trd = orad_hru[jj] / soltab_horad_potsw[jj]
-
             # If it's the first julian day of the water year, several
             # variables need to be reset:
             # - reset the previous snow water eqivalent plus new snow to 0
@@ -733,15 +736,8 @@ def _calculate_numpy(
             # TIME PERIOD
             # **************************************************************
 
-            # By default, the precipitation added to snowpack, snowmelt,
-            # and snow evaporation are 0.
-            # JLM: this could happen outside the loop
-            pk_precip[jj] = zero  # [inches]
-            snowmelt[jj] = zero  # [inches]
-            snow_evap[jj] = zero  # [inches]
-            frac_swe[jj] = zero
-            ai[jj] = zero
-            tcal[jj] = zero
+            # default assumption
+            pptmix_nopack[jj] = False
 
             # If the day of the water year is beyond the forced melt day
             # indicated by the parameter, then set the flag indicating melt
@@ -764,20 +760,15 @@ def _calculate_numpy(
             #     print(f"pk_ice 0 : {pk_ice[dbgind]}")
             #     print(f"tcal 0 : {tcal[dbgind]}")
 
-            if pkwater_equiv[jj] < epsilon64:
-                # No existing snowpack
-                if not newsnow[jj]:
-                    # Skip the HRU if there is no snowpack and no new snow
-                    # Reset to be sure it is zero if snowpack melted on last
-                    # timestep.
-                    snowcov_area[jj] = zero
-                    continue
-                else:
-                    # We ahave new snow; the initial snow-covered area is
-                    # complete (1)
-                    # JLM: why set this here? just for the case of no existing
-                    # snow? This might be removable.
-                    snowcov_area[jj] = one  # [fraction of area]
+            if (pkwater_equiv[jj] < epsilon64) and (newsnow[jj] == 0):
+                # Skip the HRU if there is no snowpack and no new snow
+                # Reset to be sure it is zero if snowpack melted on last
+                # timestep.
+                snowcov_area[jj] = zero
+                continue
+
+            if newsnow[jj] and (pkwater_equiv[jj] < epsilon64):
+                snowcov_area[jj] = one
 
             # <<
             # HRU STEP 1 - DEAL WITH PRECIPITATION AND ITS EFFECT ON THE WATER
@@ -785,6 +776,7 @@ def _calculate_numpy(
             # ****************************************************************
             month_ind = current_month - 1
 
+            # PRMS conditonal moved inside function
             (
                 freeh2o[jj],
                 iasw[jj],
@@ -906,13 +898,16 @@ def _calculate_numpy(
                 #     )
                 #     print(f"tcal 3 : {tcal[dbgind]}")
 
+            if pkwater_equiv[jj] > zero:
                 # HRU STEP 4 - DETERMINE RADIATION FLUXES AND SNOWPACK
                 #              STATES NECESSARY FOR ENERGY BALANCE
                 # **********************************************************
                 (
                     freeh2o[jj],
+                    iasw[jj],
                     iso[jj],
                     lso[jj],
+                    mso[jj],
                     pk_def[jj],
                     pk_den[jj],
                     pk_depth[jj],
@@ -921,8 +916,9 @@ def _calculate_numpy(
                     pkwater_equiv[jj],
                     pss[jj],
                     tcal[jj],
+                    snowmelt[jj],
                 ) = calc_step_4(
-                    trd,
+                    trd[jj],
                     calc_calin=calc_calin,
                     calc_caloss=calc_caloss,
                     calc_snowbal=calc_snowbal,
@@ -1035,13 +1031,13 @@ def _calculate_numpy(
                     # insufficient to reset albedo, then reduce the cumulative
                     # new snow by the amount melted during the period (but
                     # don't let it be negative).
-                    if lst[jj]:
+                    if lst[jj] > 0:
                         snsv[jj] = snsv[jj] - snowmelt[jj]
 
                         if snsv[jj] < zero:
                             snsv[jj] = zero
 
-            # <<<<
+            # <<<<  The if starting between steps 3 & 4 ends here
             # LAST check to clear out all arrays if packwater is gone
             if pkwater_equiv[jj] <= zero:
                 # if verbose:
@@ -1073,17 +1069,30 @@ def _calculate_numpy(
                 snowcov_areasv[jj] = zero  # rsr, not in original code
                 ai[jj] = zero
                 frac_swe[jj] = zero
+                scrv[jj] = zero
+                pksv[jj] = zero
 
-        # <<
+        # << end of space loop and previous if
 
         pkwater_equiv_change[:] = pkwater_equiv - pkwater_ante
         freeh2o_change[:] = freeh2o - freeh2o_prev
         pk_ice_change[:] = pk_ice - pk_ice_prev
 
-        wh_through = (
-            ((pk_ice_prev + freeh2o_prev) <= epsilon64) & ~newsnow
-        ) | (pptmix_nopack == 1)
-        through_rain[:] = np.where(wh_through, net_rain, zero)
+        nearzero = 1.0e-6
+        cond1 = net_ppt > zero
+        cond2 = pptmix_nopack != 0
+        cond3 = snowmelt < nearzero
+        cond4 = pkwater_equiv < epsilon32
+        cond5 = snow_evap < nearzero
+        cond6 = net_snow < nearzero
+        # reverse order from the if statements
+        through_rain[:] = np.where(
+            cond1 & cond3 & cond4 & cond6, net_rain, zero
+        )
+        through_rain[:] = np.where(
+            cond1 & cond3 & cond4 & cond5, net_ppt, through_rain
+        )
+        through_rain[:] = np.where(cond1 & cond2, net_rain, through_rain)
 
         return (
             ai,
@@ -1203,10 +1212,10 @@ def _calc_ppt_to_pack(
                 snowmelt,
             )
 
-        caln = zero
-        calpr = zero
-        calps = zero
-        pndz = zero
+        # caln = zero
+        # calpr = zero
+        # calps = zero
+        # pndz = zero
 
         tsnow = tavgc  # [degrees C]
 
@@ -1575,16 +1584,6 @@ def _calc_calin(
             pk_def = pk_def - cal  # [cal/cm^2]
             pk_temp = -1 * pk_def / (pkwater_equiv * 1.27)  # [degrees C]
 
-        elif abs(dif) < epsilon32:
-            # JLM: I moved this from an else at the bottom to this
-            #      point in the conditional chain.
-            # JLM: The test had been equality with zero, changed to
-            #      less than epsilon32.
-            # (2) Just enough heat to overcome heat deficit
-            # Set temperature and heat deficit to zero. the pack is "ripe"
-            pk_temp = zero  # [degrees C]
-            pk_def = zero  # [cal/cm^2]
-
         elif dif > zero:
             # (3) More than enough heat to overcome heat deficit (melt ice)...
             # Calculate the potential amount of snowmelt from excess heat in
@@ -1701,6 +1700,17 @@ def _calc_calin(
                     pss = pkwater_equiv  # [inches]
 
         # <<<
+        else:  # abs(dif) < epsilon64:
+            # JLM: The test had been equality with zero, changed to
+            #      less than epsilon64.
+            # (2) Just enough heat to overcome heat deficit
+            # Set temperature and heat deficit to zero. the pack is "ripe"
+            pk_temp = zero  # [degrees C]
+            pk_def = zero  # [cal/cm^2]
+
+        if not (pkwater_equiv > zero):
+            pk_den = zero
+
         return (
             freeh2o,
             iasw,
@@ -1854,10 +1864,11 @@ def _calc_snowcov(
         # <
         # Calculate the ratio of the current packwater equivalent to the
         # maximum packwater equivalent for the given snowpack.
-        if ai == zero:
-            frac_swe = zero
-        else:
+        if ai > epsilon64:
             frac_swe = pkwater_equiv / ai  # [fraction]
+            frac_swe = min(one, frac_swe)
+        else:
+            frac_swe = zero
 
         # <
         # There are 3 potential conditions for the snow area curve:
@@ -2204,7 +2215,6 @@ def _calc_snalbedo(
                         # reset the albedo.
                         # Reset the albedo states.
                         slst = zero  # [days]
-                        # lst = 0  # [flag]
                         lst = False  # [flag]
                         snsv = zero  # [inches]
 
@@ -2235,7 +2245,7 @@ def _calc_snalbedo(
             # threshold.
             # 4 options below (if-then, elseif, elseif, else)
 
-            if pptmix < one:
+            if pptmix == 0:
                 # (3.1) This is not a mixed event...
                 # During the accumulation season, the threshold for resetting
                 # the albedo does not apply if there is a snow-only event.
@@ -2439,12 +2449,6 @@ def _calc_step_4(
         # Save the current value of emissivity
         esv = emis  # [fraction of radiation]
 
-        # The incoming shortwave radiation is the HRU radiation adjusted by the
-        # albedo (some is reflected back into the atmoshphere) and the
-        # transmission coefficient (some is intercepted by the winter
-        # vegetative canopy)
-        swn = swrad * (one - albedo) * rad_trncf  # [cal/cm^2] or [Langleys]
-
         # Set the convection-condensation for a half-day interval
         cec = cecn_coef * 0.5  # [cal/(cm^2 degC)] or [Langleys/degC]
 
@@ -2454,49 +2458,6 @@ def _calc_step_4(
             # RSR: cov_type==4 is valid for trees (coniferous)
             cec = cec * 0.5  # [cal/(cm^2 degC)] or [Langleys/degC]
 
-        # <
-        # Calculate the new snow depth (Riley et al. 1973)
-        # RSR: the following 3 lines of code were developed by Rob Payn,
-        # on 7/10/2013
-        # The snow depth depends on the previous snow pack water equivalent
-        # plus the current net snow.
-        pss = pss + net_snow  # [inches]
-        dpt_before_settle = pk_depth + net_snow * deninv
-        dpt1 = dpt_before_settle + settle_const * (
-            (pss * denmaxinv) - dpt_before_settle
-        )
-
-        # RAPCOMMENT - CHANGED TO THE APPROPRIATE FINITE DIFFERENCE
-        # APPROXIMATION OF SNOW DEPTH
-        # JLM: pk_depth is prognostic here
-        pk_depth = dpt1  # [inches]
-
-        # Calculate the snowpack density
-        if dpt1 > zero:
-            pk_den = pkwater_equiv / dpt1
-        else:
-            pk_den = zero  # [inch water equiv / inch depth]
-
-        # <
-        # The effective thermal conductivity is approximated (empirically)
-        # as zero077 times (snowpack density)^2 [cal / (sec g degC)] Therefore,
-        # the effective conductivity term (inside the square root) in the
-        # equation for conductive heat exchange can be calculated as follows:
-        #   (zero077 * pk_den^2) / (pk_den * 0.5)
-        # where 0.5 is the specific heat of ice [cal / (g degC)]
-        # this simplifies to the following
-        effk = 0.0154 * pk_den  # [unitless]
-
-        # 13751 is the number of seconds in 12 hours over pi
-        # So for a half day, to calculate the conductive heat exchange per cm
-        # of snow per cm^2 area per degree temperature difference is the
-        # following
-        # In effect, multiplying cst times the temperature gradient gives the
-        # heatexchange by heat conducted (calories) per square cm of snowpack
-        cst = pk_den * (
-            np.sqrt(effk * 13751.0)
-        )  # [cal/(cm^2 degC)] or [Langleys / degC]
-
         # Check whether to force spring melt
         # Spring melt is forced if time is before the melt-force day and after
         # the melt-look day (parameters). If between these dates, the spring
@@ -2539,75 +2500,139 @@ def _calc_step_4(
         # Set the flag indicating night time
         niteda = 1  # [flag]
 
-        # No shortwave (solar) radiation at night.
-        sw = zero  # [cal / cm^2] or [Langleys]
-
         # Temperature is halfway between the minimum and average temperature
         # for the day.
         temp = (tminc + tavgc) * 0.5
 
-        # Track total heat flux from both night and day periods
+        if pkwater_equiv > zero:
+            # The incoming shortwave radiation is the HRU radiation adjusted by
+            # the albedo (some is reflected back into the atmoshphere) and the
+            # transmission coefficient (some is intercepted by the winter
+            # vegetative canopy)
+            swn = (
+                swrad * (one - albedo) * rad_trncf
+            )  # [cal/cm^2] or [Langleys]
+
+            # Calculate the new snow depth (Riley et al. 1973)
+            # RSR: the following 3 lines of code were developed by Rob Payn,
+            # on 7/10/2013
+            # The snow depth depends on the previous snow pack water equivalent
+            # plus the current net snow.
+            pss = pss + net_snow  # [inches]
+            # deninv is the inverse of den_init
+            dpt_before_settle = pk_depth + net_snow * deninv
+            dpt1 = dpt_before_settle + settle_const * (
+                (pss * denmaxinv) - dpt_before_settle
+            )
 
-        # Calculate the night time energy balance
-        (
-            tcal,
-            freeh2o,
-            pk_def,
-            pk_ice,
-            pk_temp,
-            pkwater_equiv,
-        ) = calc_snowbal(
-            niteda=niteda,
-            cec=cec,
-            cst=cst,
-            esv=esv,
-            sw=sw,
-            temp=temp,
-            trd=trd,
-            calc_calin=calc_calin,
-            calc_caloss=calc_caloss,
-            canopy_covden=canopy_covden,
-            den_max=den_max,
-            denmaxinv=denmaxinv,
-            emis_noppt=emis_noppt,
-            freeh2o=freeh2o,
-            freeh2o_cap=freeh2o_cap,
-            hru_ppt=hru_ppt,
-            iasw=iasw,
-            pk_def=pk_def,
-            pk_den=pk_den,
-            pk_depth=pk_depth,
-            pk_ice=pk_ice,
-            pk_temp=pk_temp,
-            pkwater_equiv=pkwater_equiv,
-            pss=pss,
-            pst=pst,
-            snowcov_area=snowcov_area,
-            snowmelt=snowmelt,
-            tcal=tcal,
-            tstorm_mo=tstorm_mo,
-        )
+            # RAPCOMMENT - CHANGED TO THE APPROPRIATE FINITE DIFFERENCE
+            # APPROXIMATION OF SNOW DEPTH
+            # JLM: pk_depth is prognostic here
+            pk_depth = dpt1  # [inches]
+
+            # Calculate the snowpack density
+            if dpt1 > zero:
+                pk_den = pkwater_equiv / dpt1
+            else:
+                pk_den = zero  # [inch water equiv / inch depth]
+
+            # <
+            # The effective thermal conductivity is approximated (empirically)
+            # as zero077 times (snowpack density)^2 [cal / (sec g degC)] Therefore,
+            # the effective conductivity term (inside the square root) in the
+            # equation for conductive heat exchange can be calculated as follows:
+            #   (zero077 * pk_den^2) / (pk_den * 0.5)
+            # where 0.5 is the specific heat of ice [cal / (g degC)]
+            # this simplifies to the following
+            effk = 0.0154 * pk_den  # [unitless]
+
+            # 13751 is the number of seconds in 12 hours over pi
+            # So for a half day, to calculate the conductive heat exchange per cm
+            # of snow per cm^2 area per degree temperature difference is the
+            # following
+            # In effect, multiplying cst times the temperature gradient gives the
+            # heatexchange by heat conducted (calories) per square cm of snowpack
+            cst = pk_den * (
+                np.sqrt(effk * 13751.0)
+            )  # [cal/(cm^2 degC)] or [Langleys / degC]
+
+            # No shortwave (solar) radiation at night.
+            sw = zero  # [cal / cm^2] or [Langleys]
+
+            # Calculate the night time energy balance
+            (
+                tcal,
+                freeh2o,
+                iasw,
+                pk_def,
+                pk_den,
+                pk_ice,
+                pk_depth,
+                pk_temp,
+                pss,
+                pst,
+                snowmelt,
+                pkwater_equiv,
+            ) = calc_snowbal(
+                niteda=niteda,
+                cec=cec,
+                cst=cst,
+                esv=esv,
+                sw=sw,
+                temp=temp,
+                trd=trd,
+                calc_calin=calc_calin,
+                calc_caloss=calc_caloss,
+                canopy_covden=canopy_covden,
+                den_max=den_max,
+                denmaxinv=denmaxinv,
+                emis_noppt=emis_noppt,
+                freeh2o=freeh2o,
+                freeh2o_cap=freeh2o_cap,
+                hru_ppt=hru_ppt,
+                iasw=iasw,
+                pk_def=pk_def,
+                pk_den=pk_den,
+                pk_depth=pk_depth,
+                pk_ice=pk_ice,
+                pk_temp=pk_temp,
+                pkwater_equiv=pkwater_equiv,
+                pss=pss,
+                pst=pst,
+                snowcov_area=snowcov_area,
+                snowmelt=snowmelt,
+                tcal=tcal,
+                tstorm_mo=tstorm_mo,
+            )
+            # tcal set directly by calc_snowbal above [cal/cm^2] or [Langleys]
 
-        # [cal/cm^2] or [Langleys]
+        # <
+        # iswn = zero on ly a glacier variable
 
         # Compute energy balance for day period (if the snowpack still exists)
         # THIS SHOULD HAPPEN IN SNOBAL
-        if pkwater_equiv > zero:
-            # Set the flag indicating daytime
-            niteda = 2  # [flag]
+        # Set the flag indicating daytime
+        niteda = 2  # [flag]
+        # Temperature is halfway between the maximum and average
+        # temperature for the day.
+        temp = (tmaxc + tavgc) * 0.5  # [degrees C]
 
+        if pkwater_equiv > zero:
             # Set shortwave radiation as calculated earlier
             sw = swn  # [cal/cm^2] or [Langleys]
 
-            # Temperature is halfway between the maximum and average
-            # temperature for the day.
-            temp = (tmaxc + tavgc) * 0.5  # [degrees C]
             (
                 cals,
                 freeh2o,
+                iasw,
                 pk_def,
+                pk_den,
                 pk_ice,
+                pk_depth,
                 pk_temp,
+                pss,
+                pst,
+                snowmelt,
                 pkwater_equiv,
             ) = calc_snowbal(
                 niteda=niteda,
@@ -2647,8 +2672,10 @@ def _calc_step_4(
         # <
         return (
             freeh2o,
+            iasw,
             iso,
             lso,
+            mso,
             pk_def,
             pk_den,
             pk_depth,
@@ -2657,6 +2684,7 @@ def _calc_step_4(
             pkwater_equiv,
             pss,
             tcal,
+            snowmelt,
         )
 
     @staticmethod
@@ -2832,9 +2860,15 @@ def _calc_snowbal(
             return (
                 cal,
                 freeh2o,
+                iasw,
                 pk_def,
+                pk_den,
                 pk_ice,
+                pk_depth,
                 pk_temp,
+                pss,
+                pst,
+                snowmelt,
                 pkwater_equiv,
             )
 
@@ -3003,9 +3037,15 @@ def _calc_snowbal(
         return (
             cal,
             freeh2o,
+            iasw,
             pk_def,
+            pk_den,
             pk_ice,
+            pk_depth,
             pk_temp,
+            pss,
+            pst,
+            snowmelt,
             pkwater_equiv,
         )
 
@@ -3048,8 +3088,6 @@ def _calc_snowevap(
             # Set the evaporation to the pack water equivalent and set all
             # snowpack
             # variables to no-snowpack values.
-            snow_evap = pkwater_equiv  # [inches]
-
             snow_evap = pkwater_equiv  # [inches]
             pkwater_equiv = zero  # [inches]
             pk_ice = zero  # [inches]
@@ -3068,6 +3106,10 @@ def _calc_snowevap(
                 # RAPCOMMENT - CHANGED TO CHECK FOR NEGATIVE PACK ICE
                 # If all pack ice is removed, then there cannot be a heat
                 # deficit.
+
+                # JLM: mass balance fix only in our 5.2.1 prms version
+                freeh2o = freeh2o + pk_ice
+
                 pk_ice = zero
                 pk_def = zero
                 pk_temp = zero
@@ -3095,8 +3137,10 @@ def _calc_snowevap(
                 #             f"snowevap: {pkwater_equiv}"
                 #         )
                 # #  <<
-
-                pkwater_equiv = zero
+                # zero of pkwater_equiv is inside a debug statement that
+                # is dosent seem to be triggered in test runs
+                # pkwater_equiv = zero
+                pass
 
             # <
             snow_evap = zero
@@ -3114,12 +3158,11 @@ def _calc_snowevap(
                     # if verbose:
                     #     if pkwater_equiv < -epsilon64:
                     #         print(
-                    #             "snowpack issue 2, negative pkwater_equiv in "
+                    #             "snowpack issue 2, negative pkwater_equiv in"
                     #             f"snowevap: {pkwater_equiv}"
                     #         )
                     # # <<
-
-                    # To be sure negative snowpack is ignored
+                    # This zero of pkwater IS NOT in the debug statement
                     pkwater_equiv = zero
 
                 # <
diff --git a/pywatershed/utils/csv_utils.py b/pywatershed/utils/csv_utils.py
index 73a8dcf7..a0c9fffc 100644
--- a/pywatershed/utils/csv_utils.py
+++ b/pywatershed/utils/csv_utils.py
@@ -124,7 +124,7 @@ def to_netcdf(
         clobber: bool = True,
         zlib: bool = True,
         complevel: int = 4,
-        chunk_sizes: dict = {"time": 30, "hruid": 0},
+        chunk_sizes: dict = {"time": 30, "nhm_id": 0, "nhm_seg": 0},
     ) -> None:
         """Output the csv output data to a netcdf file
 
@@ -192,14 +192,17 @@ def to_netcdf(
             ids = self._coordinates[dim_name]
             nids = len(ids)
 
+            dims = ("time", dim_name)
+            chunk_sizes_var = [chunk_sizes[vv] for vv in dims]
+
             var = ds.createVariable(
                 variable_name,
                 variable_type,
-                ("time", dim_name),
+                dims,
                 fill_value=nc4.default_fillvals[variable_type],  # JLM: sus
                 zlib=zlib,
                 complevel=complevel,
-                chunksizes=tuple(chunk_sizes.values()),
+                chunksizes=tuple(chunk_sizes_var),
             )
             # add additional meta data
             if self.meta.is_available(variable_name):
diff --git a/pywatershed/utils/prms5util.py b/pywatershed/utils/prms5util.py
index b3133125..a732651b 100644
--- a/pywatershed/utils/prms5util.py
+++ b/pywatershed/utils/prms5util.py
@@ -249,6 +249,7 @@ def __init__(
         soltab_file: fileish,
         output_dir: fileish = None,
         nhm_ids: np.ndarray = None,
+        chunk_sizes: dict = {"doy": 0, "nhm_id": 0},
     ):
         self.soltab_file = soltab_file
         self.output_dir = output_dir
@@ -273,7 +274,7 @@ def __init__(
         ) = load_soltab_debug(self.soltab_file)
 
         if self.output_dir:
-            self.to_netcdf()
+            self.to_netcdf(chunk_sizes=chunk_sizes)
 
         return
 
@@ -283,7 +284,7 @@ def to_netcdf(
         nhm_ids: np.ndarray = None,
         zlib: bool = True,
         complevel: int = 4,
-        chunk_sizes: dict = {"doy": 0, "hruid": 0},
+        chunk_sizes: dict = {"doy": 0, "nhm_id": 0},
     ):
         # This is just different enough to make it it's own thing compared
         # to the CSV outputs of PRMS.
@@ -320,14 +321,17 @@ def to_netcdf(
             )
             hruid[:] = self.spatial_coord
 
+            dims = ("doy", self.spatial_coord_name)
+            chunk_sizes_var = [chunk_sizes[vv] for vv in dims]
+
             variables[vv] = ds.createVariable(
                 vv,
                 "f4",
-                ("doy", self.spatial_coord_name),
+                dims,
                 fill_value=nc4.default_fillvals["f4"],
                 zlib=zlib,
                 complevel=complevel,
-                chunksizes=tuple(chunk_sizes.values()),
+                chunksizes=tuple(chunk_sizes_var),
             )
             variables[vv][:] = self[vv]
             ds.close()
diff --git a/pywatershed/utils/time_utils.py b/pywatershed/utils/time_utils.py
index 87581b90..dda65c45 100644
--- a/pywatershed/utils/time_utils.py
+++ b/pywatershed/utils/time_utils.py
@@ -18,27 +18,34 @@ def datetime_year(dt64: np.datetime64) -> int:
     return dt64.astype("datetime64[Y]").astype(int) + 1970
 
 
-def datetime_month(dt64: np.datetime64) -> int:
+def datetime_month(dt64: np.datetime64, zero_based=False) -> int:
     """Get the month from np.datetime64"""
-    return dt64.astype("datetime64[M]").astype(int) % 12 + 1
+    return dt64.astype("datetime64[M]").astype(int) % 12 + _offset(zero_based)
 
 
-def datetime_doy(dt64: np.datetime64) -> int:
+def datetime_doy(dt64: np.datetime64, zero_based=False) -> int:
     """Get day of year from np.datetime64"""
     return (dt64 - dt64.astype("datetime64[Y]")).astype(
         "timedelta64[D]"
-    ).astype(int) + 1
+    ).astype(int) + _offset(zero_based)
 
 
-def datetime_dowy(dt64: np.datetime64) -> int:
+def datetime_dowy(dt64: np.datetime64, zero_based=False) -> int:
     """Get day of water year from np.datetime64"""
     year_start = datetime_year(dt64)
     if datetime_month(dt64) < 10:
         year_start -= 1
     diff = dt64 - np.datetime64(f"{year_start}-10-01")
-    return diff.astype("timedelta64[D]").astype(int)
+    return diff.astype("timedelta64[D]").astype(int) + _offset(zero_based)
 
 
 def datetime_epiweek(dt64: np.datetime64) -> int:
     """Get CDC eipweek [1, 53] from np.datetime64"""
     return ew.Week.fromdate(dt64_to_dt(dt64)).week
+
+
+def _offset(zero_based: bool):
+    if zero_based:
+        return 0
+    else:
+        return 1


Test patch:
diff --git a/autotest/test_control.py b/autotest/test_control.py
index 4b215141..946425c1 100644
--- a/autotest/test_control.py
+++ b/autotest/test_control.py
@@ -80,7 +80,7 @@ def test_control_simple(control_simple):
         year = year if month >= 10 else year - 1
         wy_start = np.datetime64(f"{year}-10-01")
         dowy = (current_time - wy_start).astype("timedelta64[D]")
-        assert dowy == control_simple.current_dowy
+        assert dowy == (control_simple.current_dowy - 1)
 
         prev_time = control_simple.current_time
 
diff --git a/autotest/test_model.py b/autotest/test_model.py
index 0f15ad56..8cc1bdd9 100644
--- a/autotest/test_model.py
+++ b/autotest/test_model.py
@@ -287,18 +287,18 @@ def test_model(domain, model_args, tmp_path):
         9: {
             "PRMSChannel": {
                 "seg_outflow": {
-                    "drb_2yr": 1430.6364027142613,
-                    "hru_1": 13.416914151483681,
-                    "ucb_2yr": 1694.5412856707849,
+                    "drb_2yr": 1517.232887980279,
+                    "hru_1": 13.696918669514927,
+                    "ucb_2yr": 1694.5697712423928,
                 },
             },
         },
         99: {
             "PRMSChannel": {
                 "seg_outflow": {
-                    "drb_2yr": 1588.1444684289775,
-                    "hru_1": 19.596412903692578,
-                    "ucb_2yr": 407.2200022510677,
+                    "drb_2yr": 2350.499659332901,
+                    "hru_1": 22.874414994530095,
+                    "ucb_2yr": 733.2293013532435,
                 },
             },
         },
diff --git a/autotest/test_nhm_self_drive.py b/autotest/test_nhm_self_drive.py
index c0a2ecc5..f8b734ba 100644
--- a/autotest/test_nhm_self_drive.py
+++ b/autotest/test_nhm_self_drive.py
@@ -22,10 +22,14 @@
 
 
 def test_drive_indiv_process(domain, tmp_path):
-    """Use output from a full NHM run to drive each of the indiv processes
-    separately: self-driving
+    """Output of a full pywatershed NHM drives indiv process models separately
+
+    The results from the full model should be consistent with the results from
+    the individual models, else there is likely something wrong with the
+    full model.
     """
-    # Full NHM output
+
+    # Run a full pws NHM to use its output to drive individual processes
     nhm_output_dir = pl.Path(tmp_path) / "nhm_output"
 
     params = pws.parameters.PrmsParameters.load(domain["param_file"])
@@ -44,10 +48,9 @@ def test_drive_indiv_process(domain, tmp_path):
     nhm.run(finalize=True)
     del nhm, params, control
 
-    # individual process models
+    # run individual process models
     for proc in nhm_processes:
-        # proc = pws.PRMSRunoff  # TODO: fix this one ASAP
-        if proc in [pws.PRMSSolarGeometry, pws.PRMSAtmosphere, pws.PRMSRunoff]:
+        if proc in [pws.PRMSSolarGeometry, pws.PRMSAtmosphere]:
             # These are not driven by outputs of above, only external outputs
             # or known/static inputs
             continue
@@ -82,10 +85,11 @@ def test_drive_indiv_process(domain, tmp_path):
             ans = xr.open_dataset(nhm_output_dir / f"{vv}.nc")
 
             # Leaving the commented to diagnose what PRMSRunoff later.
-            # try:
-            xr.testing.assert_allclose(res, ans)
-            # except:
-            #     print(vv)
+            try:
+                xr.testing.assert_allclose(res, ans)
+            except:
+                print(vv, abs(res - ans).max())
+                print(vv, (abs(res - ans) / ans).max())
 
             del res, ans
 
diff --git a/autotest/test_preprocess_csv.py b/autotest/test_preprocess_csv.py
index 0364411c..8ce8bfb1 100644
--- a/autotest/test_preprocess_csv.py
+++ b/autotest/test_preprocess_csv.py
@@ -7,6 +7,8 @@
 
 from pywatershed import CsvFile
 
+# these CSV files are protected from deletion in CI by
+# test_data/scripts/test_remove_csvs.py
 csv_test_vars = ["hru_ppt", "intcp_stor", "potet", "gwres_stor"]
 
 
diff --git a/autotest/test_prms_channel.py b/autotest/test_prms_channel.py
index 299db4a7..5160e816 100644
--- a/autotest/test_prms_channel.py
+++ b/autotest/test_prms_channel.py
@@ -2,11 +2,16 @@
 
 import pytest
 
+from pywatershed.base.adapter import adapter_factory
 from pywatershed.base.control import Control
 from pywatershed.base.parameters import Parameters
 from pywatershed.hydrology.prms_channel import PRMSChannel, has_prmschannel_f
 from pywatershed.parameters import PrmsParameters
-from pywatershed.utils.netcdf_utils import NetCdfCompare
+from utils_compare import compare_in_memory, compare_netcdfs
+
+# compare in memory (faster) or full output files?
+compare_output_files = False
+rtol = atol = 1.0e-7
 
 fail_fast = False
 
@@ -67,52 +72,39 @@ def test_compare_prms(
         budget_type="error",
         calc_method=calc_method,
     )
-    nc_parent = tmp_path / domain["domain_name"]
-    channel.initialize_netcdf(nc_parent)
-    # test that init netcdf twice raises a warning
-    with pytest.warns(UserWarning):
+
+    if compare_output_files:
+        nc_parent = tmp_path / domain["domain_name"]
         channel.initialize_netcdf(nc_parent)
+        # test that init netcdf twice raises a warning
+        with pytest.warns(UserWarning):
+            channel.initialize_netcdf(nc_parent)
+
+    else:
+        answers = {}
+        for var in PRMSChannel.get_variables():
+            var_pth = output_dir / f"{var}.nc"
+            answers[var] = adapter_factory(
+                var_pth, variable_name=var, control=control
+            )
 
     for istep in range(control.n_times):
         control.advance()
-
         channel.advance()
-
         channel.calculate(float(istep))
-
         channel.output()
+        if not compare_output_files:
+            compare_in_memory(channel, answers, atol=atol, rtol=rtol)
 
     channel.finalize()
 
-    output_compare = {}
-
-    for key in PRMSChannel.get_variables():
-        base_nc_path = output_dir / f"{key}.nc"
-        compare_nc_path = tmp_path / domain["domain_name"] / f"{key}.nc"
-        # PRMS does not output the storage change in the channel
-        if not base_nc_path.exists():
-            continue
-        output_compare[key] = (base_nc_path, compare_nc_path)
-
-    assert_error = False
-    for key, (base, compare) in output_compare.items():
-        print(f"\nbase_nc_path: {base}")
-        print(f"compare_nc_path: {compare}")
-        success, diff = NetCdfCompare(base, compare).compare()
-        if not success:
-            print(
-                f"comparison for {key} failed: "
-                + f"maximum error {diff[key][0]} "
-                + f"(maximum allowed error {diff[key][1]}) "
-                + f"in column {diff[key][2]}"
-            )
-            assert_error = True
-            if fail_fast:
-                assert False
-
-        else:
-            print(f"comparison for {key} passed")
-
-    assert not assert_error, "comparison failed"
+    if compare_output_files:
+        compare_netcdfs(
+            PRMSChannel.get_variables(),
+            tmp_path / domain["domain_name"],
+            output_dir,
+            atol=atol,
+            rtol=rtol,
+        )
 
     return
diff --git a/autotest/test_prms_groundwater.py b/autotest/test_prms_groundwater.py
index a3381fe0..f0209d8d 100644
--- a/autotest/test_prms_groundwater.py
+++ b/autotest/test_prms_groundwater.py
@@ -2,14 +2,15 @@
 
 import pytest
 
-from pywatershed.base.control import Control
-from pywatershed.base.parameters import Parameters
-from pywatershed.hydrology.prms_groundwater import (
-    PRMSGroundwater,
-    has_prmsgroundwater_f,
-)
+from pywatershed import Control, Parameters, PRMSGroundwater
+from pywatershed.base.adapter import adapter_factory
+from pywatershed.hydrology.prms_groundwater import has_prmsgroundwater_f
 from pywatershed.parameters import PrmsParameters
-from pywatershed.utils.netcdf_utils import NetCdfCompare
+from utils_compare import compare_in_memory, compare_netcdfs
+
+# compare in memory (faster) or full output files?
+compare_output_files = False
+rtol = atol = 1.0e-13
 
 calc_methods = ("numpy", "numba", "fortran")
 params = ("params_sep", "params_one")
@@ -48,7 +49,6 @@ def test_compare_prms(
 
     tmp_path = pl.Path(tmp_path)
 
-    # load csv files into dataframes
     output_dir = domain["prms_output_dir"]
     input_variables = {}
     for key in PRMSGroundwater.get_inputs():
@@ -63,49 +63,36 @@ def test_compare_prms(
         budget_type="error",
         calc_method=calc_method,
     )
-    nc_parent = tmp_path / domain["domain_name"]
-    gw.initialize_netcdf(nc_parent)
-
-    output_compare = {}
-    vars_compare = (
-        "gwres_flow",
-        "gwres_sink",
-        "gwres_stor",
-        "ssr_to_gw",
-        "soil_to_gw",
-    )
-    for key in PRMSGroundwater.get_variables():
-        if key not in vars_compare:
-            continue
-        base_nc_path = output_dir / f"{key}.nc"
-        compare_nc_path = tmp_path / domain["domain_name"] / f"{key}.nc"
-        output_compare[key] = (base_nc_path, compare_nc_path)
 
-    print(f"base_nc_path: {base_nc_path}")
-    print(f"compare_nc_path: {compare_nc_path}")
+    if compare_output_files:
+        nc_parent = tmp_path / domain["domain_name"]
+        gw.initialize_netcdf(nc_parent)
+    else:
+        answers = {}
+        for var in PRMSGroundwater.get_variables():
+            var_pth = output_dir / f"{var}.nc"
+            answers[var] = adapter_factory(
+                var_pth, variable_name=var, control=control
+            )
 
     for istep in range(control.n_times):
         control.advance()
-
         gw.advance()
-
         gw.calculate(float(istep))
-
         gw.output()
 
+        if not compare_output_files:
+            compare_in_memory(gw, answers, atol=atol, rtol=rtol)
+
     gw.finalize()
 
-    assert_error = False
-    for key, (base, compare) in output_compare.items():
-        success, diff = NetCdfCompare(base, compare).compare()
-        if not success:
-            print(
-                f"comparison for {key} failed: "
-                + f"maximum error {diff[key][0]} "
-                + f"(maximum allowed error {diff[key][1]}) "
-                + f"in column {diff[key][2]}"
-            )
-            assert_error = True
-    assert not assert_error, "comparison failed"
+    if compare_output_files:
+        compare_netcdfs(
+            PRMSGroundwater.get_variables(),
+            tmp_path / domain["domain_name"],
+            output_dir,
+            atol=atol,
+            rtol=rtol,
+        )
 
     return
diff --git a/autotest/test_prms_runoff.py b/autotest/test_prms_runoff.py
index 88a25517..f0fe0de6 100644
--- a/autotest/test_prms_runoff.py
+++ b/autotest/test_prms_runoff.py
@@ -51,6 +51,15 @@ def test_compare_prms(
         "sroff",
         "dprst_evap_hru",
         "hru_impervevap",
+        "dprst_insroff_hru",
+        "dprst_stor_hru",
+        "contrib_fraction",
+        "hru_sroffp",
+        "hru_sroffi",
+        # "hru_impervstor_change",
+        "dprst_sroff_hru",
+        "dprst_insroff_hru",
+        # "dprst_stor_hru_change",
     ]
     output_dir = domain["prms_output_dir"]
 
@@ -72,7 +81,7 @@ def test_compare_prms(
         parameters=parameters,
         calc_method=calc_method,
         **input_variables,
-        budget_type="warn",  # intermittent errors currently
+        budget_type="error",
     )
 
     all_success = True
@@ -90,7 +99,7 @@ def test_compare_prms(
         failfast = True
         detailed = True
         if check:
-            atol = 1.0e-5
+            atol = 1.0e-10
             success = check_timestep_results(
                 runoff, istep, ans, atol, detailed
             )
@@ -113,7 +122,7 @@ def check_timestep_results(storageunit, istep, ans, atol, detailed=False):
     for key in ans.keys():
         a1 = ans[key].current
         a2 = storageunit[key]
-        success = np.isclose(a1, a2, atol=atol).all()
+        success = np.isclose(a1, a2, atol=atol, rtol=atol).all()
         if not success:
             all_success = False
             diff = a1 - a2
diff --git a/autotest/test_prms_solar_geom.py b/autotest/test_prms_solar_geom.py
index 4d40cd9a..f2a4560c 100644
--- a/autotest/test_prms_solar_geom.py
+++ b/autotest/test_prms_solar_geom.py
@@ -31,7 +31,6 @@ def parameters(domain, request):
     return params
 
 
-@pytest.mark.xfail
 @pytest.mark.parametrize(
     "from_prms_file", (True, False), ids=("from_prms_file", "compute")
 )
diff --git a/autotest/utils_compare.py b/autotest/utils_compare.py
new file mode 100644
index 00000000..dd379184
--- /dev/null
+++ b/autotest/utils_compare.py
@@ -0,0 +1,121 @@
+import pathlib as pl
+
+import numpy as np
+import xarray as xr
+
+import pywatershed as pws
+
+
+def assert_allclose(
+    actual: np.ndarray,
+    desired: np.ndarray,
+    rtol: float = 1.0e-15,
+    atol: float = 1.0e-15,
+    equal_nan: bool = True,
+    strict: bool = False,
+    also_check_w_np: bool = True,
+    error_message: str = "Comparison unsuccessful (default message)",
+):
+    """Reinvent np.testing.assert_allclose to get useful diagnostincs in debug
+
+    Args:
+        actual: Array obtained.
+        desired: Array desired.
+        rtol: Relative tolerance.
+        atol: Absolute tolerance.
+        equal_nan: If True, NaNs will compare equal.
+        strict: If True, raise an ``AssertionError`` when either the shape or
+            the data type of the arguments does not match. The special
+            handling of scalars mentioned in the Notes section is disabled.
+        also_check_w_np: first check using np.testing.assert_allclose using
+            the same options.
+    """
+
+    if also_check_w_np:
+        np.testing.assert_allclose(
+            actual,
+            desired,
+            rtol=rtol,
+            atol=atol,
+            equal_nan=equal_nan,
+            # strict=strict,  # to add to newer versions of numpy
+        )
+
+    if strict:
+        assert actual.shape == desired.shape
+        assert isinstance(actual, type(desired))
+        assert isinstance(desired, type(actual))
+
+    if equal_nan:
+        actual_nan = np.where(np.isnan(actual), True, False)
+        desired_nan = np.where(np.isnan(desired), True, False)
+        assert (actual_nan == desired_nan).all()
+
+    abs_diff = abs(actual - desired)
+    rel_abs_diff = abs_diff / desired
+
+    abs_close = abs_diff < atol
+    rel_close = rel_abs_diff < rtol
+    rel_close = np.where(np.isnan(rel_close), False, rel_close)
+
+    close = abs_close | rel_close
+
+    assert close.all()
+
+
+def compare_in_memory(
+    process: pws.base.Process,
+    answers: dict[pws.base.adapter.AdapterNetcdf],
+    rtol: float = 1.0e-15,
+    atol: float = 1.0e-15,
+    equal_nan: bool = True,
+    strict: bool = False,
+    also_check_w_np: bool = True,
+    error_message: str = None,
+):
+    # TODO: docstring
+    for var in process.get_variables():
+        answers[var].advance()
+        assert_allclose(
+            process[var],
+            answers[var].current.data,
+            atol=atol,
+            rtol=rtol,
+            equal_nan=equal_nan,
+            strict=strict,
+            also_check_w_np=also_check_w_np,
+            error_message=error_message,
+        )
+
+
+def compare_netcdfs(
+    var_list: list,
+    results_dir: pl.Path,
+    answers_dir: pl.Path,
+    rtol: float = 1.0e-15,
+    atol: float = 1.0e-15,
+    equal_nan: bool = True,
+    strict: bool = False,
+    also_check_w_np: bool = True,
+    error_message: str = None,
+):
+    # TODO: docstring
+    # TODO: improve error message
+    # TODO: collect failures in a try and report at end
+    for var in var_list:
+        answer = xr.open_dataarray(answers_dir / f"{var}.nc")
+        result = xr.open_dataarray(results_dir / f"{var}.nc")
+
+        if error_message is None:
+            error_message = f"Comparison of variable '{var}' was unsuccessful"
+
+        assert_allclose(
+            actual=result.values,
+            desired=answer.values,
+            rtol=rtol,
+            atol=atol,
+            equal_nan=equal_nan,
+            strict=strict,
+            also_check_w_np=also_check_w_np,
+            error_message=error_message,
+        )
diff --git a/test_data/common/control.multi_hru b/test_data/common/control.multi_hru
index df624841..effd9667 100644
--- a/test_data/common/control.multi_hru
+++ b/test_data/common/control.multi_hru
@@ -238,7 +238,7 @@ nhruOutON_OFF
 nhruOutVars
 1
 1
-78
+83
 ####
 nhruOut_freq
 1
@@ -251,13 +251,14 @@ nhruOut_format
 4
 ####
 nhruOutVar_names
-78
+83
 4
 albedo
 cap_infil_tot
 contrib_fraction
 dprst_seep_hru
 dprst_evap_hru
+dprst_insroff_hru
 dprst_sroff_hru
 dprst_stor_hru
 freeh2o
@@ -288,7 +289,11 @@ net_snow
 newsnow
 orad_hru
 perv_actet
+pk_def
+pk_den
+pk_depth
 pk_ice
+pk_temp
 pkwater_ante
 pkwater_equiv
 potet
@@ -345,7 +350,7 @@ nsegmentOutON_OFF
 nsegmentOutVars
 1
 1
-9
+10
 ####
 nsegmentOut_freq
 1
@@ -353,8 +358,9 @@ nsegmentOut_freq
 1
 ####
 nsegmentOutVar_names
-9
+10
 4
+seg_inflow
 seg_outflow
 seg_gwflow
 seg_sroff
diff --git a/test_data/common/control.single_hru b/test_data/common/control.single_hru
index e692a8b6..74ea87b3 100644
--- a/test_data/common/control.single_hru
+++ b/test_data/common/control.single_hru
@@ -238,7 +238,7 @@ nhruOutON_OFF
 nhruOutVars
 1
 1
-76
+81
 ####
 nhruOut_freq
 1
@@ -251,12 +251,13 @@ nhruOut_format
 5
 ####
 nhruOutVar_names
-76
+81
 4
 albedo
 cap_infil_tot
 contrib_fraction
 dprst_evap_hru
+dprst_insroff_hru
 dprst_seep_hru
 dprst_stor_hru
 dprst_sroff_hru
@@ -288,7 +289,11 @@ net_snow
 newsnow
 orad_hru
 perv_actet
+pk_def
+pk_den
+pk_depth
 pk_ice
+pk_temp
 pkwater_ante
 pkwater_equiv
 potet
@@ -343,7 +348,7 @@ nsegmentOutON_OFF
 nsegmentOutVars
 1
 1
-9
+10
 ####
 nsegmentOut_freq
 1
@@ -351,8 +356,9 @@ nsegmentOut_freq
 1
 ####
 nsegmentOutVar_names
-9
+10
 4
+seg_inflow
 seg_outflow
 seg_gwflow
 seg_sroff
diff --git a/test_data/drb_2yr/control.test b/test_data/drb_2yr/control.test
index df624841..effd9667 100644
--- a/test_data/drb_2yr/control.test
+++ b/test_data/drb_2yr/control.test
@@ -238,7 +238,7 @@ nhruOutON_OFF
 nhruOutVars
 1
 1
-78
+83
 ####
 nhruOut_freq
 1
@@ -251,13 +251,14 @@ nhruOut_format
 4
 ####
 nhruOutVar_names
-78
+83
 4
 albedo
 cap_infil_tot
 contrib_fraction
 dprst_seep_hru
 dprst_evap_hru
+dprst_insroff_hru
 dprst_sroff_hru
 dprst_stor_hru
 freeh2o
@@ -288,7 +289,11 @@ net_snow
 newsnow
 orad_hru
 perv_actet
+pk_def
+pk_den
+pk_depth
 pk_ice
+pk_temp
 pkwater_ante
 pkwater_equiv
 potet
@@ -345,7 +350,7 @@ nsegmentOutON_OFF
 nsegmentOutVars
 1
 1
-9
+10
 ####
 nsegmentOut_freq
 1
@@ -353,8 +358,9 @@ nsegmentOut_freq
 1
 ####
 nsegmentOutVar_names
-9
+10
 4
+seg_inflow
 seg_outflow
 seg_gwflow
 seg_sroff
diff --git a/test_data/scripts/cbh_nc_indivual.py b/test_data/generate/cbh_nc_indivual.py
similarity index 100%
rename from test_data/scripts/cbh_nc_indivual.py
rename to test_data/generate/cbh_nc_indivual.py
diff --git a/test_data/scripts/cbh_to_nc.py b/test_data/generate/cbh_to_nc.py
similarity index 100%
rename from test_data/scripts/cbh_to_nc.py
rename to test_data/generate/cbh_to_nc.py
diff --git a/test_data/scripts/conftest.py b/test_data/generate/conftest.py
similarity index 86%
rename from test_data/scripts/conftest.py
rename to test_data/generate/conftest.py
index 2ad89d89..31468dd7 100644
--- a/test_data/scripts/conftest.py
+++ b/test_data/generate/conftest.py
@@ -32,17 +32,16 @@ def pytest_addoption(parser):
 
 @pytest.fixture()
 def exe():
-    exe_name = "prms"
     platform = sys.platform.lower()
     if platform == "win32":
-        exe_name += "_win.exe"
+        exe_name = "prms_win_gfort_dbl_prec.exe"
     elif platform == "darwin":
         if processor() == "arm":
-            exe_name += "_mac_m1_intel"
+            exe_name = "prms_mac_m1_ifort_dbl_prec"
         else:
-            exe_name += "_mac"
+            exe_name = "prms_mac_intel_gfort_dbl_prec"
     elif platform == "linux":
-        exe_name += "_linux"
+        exe_name = "prms_linux_gfort_dbl_prec"
     exe_pth = pl.Path(f"../../bin/{exe_name}").resolve()
     return exe_pth
 
@@ -161,3 +160,16 @@ def pytest_generate_tests(metafunc):
         ]
         ids = [pl.Path(kk).name for kk in simulations.keys()]
         metafunc.parametrize("simulation", sims, ids=ids, scope="session")
+
+    if "final_file" in metafunc.fixturenames:
+        final_var_names = ["through_rain", "seg_lateral_inflow"]
+        final_files = [
+            pl.Path(kk) / var
+            for kk in simulations.keys()
+            for var in final_var_names
+        ]
+        ids = [ff.parent.name + ":" + ff.name for ff in final_files]
+        # these are not really file names they are domain/key_var
+        metafunc.parametrize(
+            "final_file", final_files, ids=ids, scope="session"
+        )
diff --git a/test_data/scripts/conus_2yr_subset.sh b/test_data/generate/conus_2yr_subset.sh
similarity index 100%
rename from test_data/scripts/conus_2yr_subset.sh
rename to test_data/generate/conus_2yr_subset.sh
diff --git a/test_data/generate/convert_prms_output_to_nc.py b/test_data/generate/convert_prms_output_to_nc.py
new file mode 100644
index 00000000..af560a6d
--- /dev/null
+++ b/test_data/generate/convert_prms_output_to_nc.py
@@ -0,0 +1,58 @@
+from pathlib import Path
+from filelock import FileLock
+
+import pytest
+
+from prms_convert_to_netcdf import convert_csv_to_nc, convert_soltab_to_nc
+from prms_diagnostic_variables import (
+    diagnose_simple_vars_to_nc,
+    diagnose_final_vars_to_nc,
+)
+
+
+@pytest.fixture
+def netcdf_file(csv_file):
+    """Convert CSV files from model output to NetCDF"""
+
+    var_name = csv_file.stem
+    data_dir = csv_file.parent
+    convert_csv_to_nc(var_name, data_dir)
+
+    diagnose_simple_vars_to_nc(var_name, data_dir)
+
+    return
+
+
+def make_netcdf_files(netcdf_file):
+    print(f"Created NetCDF from CSV: {netcdf_file}")
+
+
+@pytest.fixture()
+def soltab_netcdf_file(tmp_path_factory, soltab_file):
+    """Convert soltab files to NetCDF, one file for each variable"""
+    domain_dir = soltab_file.parent
+    output_dir = domain_dir / "output"
+    convert_soltab_to_nc(soltab_file, output_dir, domain_dir)
+
+
+def make_soltab_netcdf_files(soltab_netcdf_file):
+    print(f"Creating NetCDF files for soltab file {soltab_netcdf_file}")
+
+
+@pytest.fixture(scope="session")
+def final_netcdf_file(tmp_path_factory, final_file):
+    """Create NetCDF files that depend on multiple other NetCDFs"""
+
+    domain_dir = final_file.parent
+    var_name = final_file.name
+    data_dir = domain_dir / "output"
+
+    root_tmpdir = tmp_path_factory.getbasetemp().parent
+    with FileLock(root_tmpdir / "final_nc.lock"):
+        yield  # do this in session cleanup
+
+        diagnose_final_vars_to_nc(var_name, data_dir, domain_dir)
+
+
+def make_final_netcdf_files(final_netcdf_file):
+    print(f"Creating final NetCDF file {final_netcdf_file}")
diff --git a/test_data/scripts/drb_2yr_check.sh b/test_data/generate/drb_2yr_check.sh
similarity index 100%
rename from test_data/scripts/drb_2yr_check.sh
rename to test_data/generate/drb_2yr_check.sh
diff --git a/test_data/scripts/drb_2yr_subset.sh b/test_data/generate/drb_2yr_subset.sh
similarity index 100%
rename from test_data/scripts/drb_2yr_subset.sh
rename to test_data/generate/drb_2yr_subset.sh
diff --git a/test_data/generate/prms_convert_to_netcdf.py b/test_data/generate/prms_convert_to_netcdf.py
new file mode 100644
index 00000000..e3c526c3
--- /dev/null
+++ b/test_data/generate/prms_convert_to_netcdf.py
@@ -0,0 +1,54 @@
+import pathlib as pl
+
+import pywatershed as pws
+from pywatershed import CsvFile, Soltab
+
+"""This module is intended to aid/make consistent conversions of PRMS files"""
+
+
+def convert_csv_to_nc(
+    var_name: str, data_dir: pl.Path, output_dir: pl.Path = None
+):
+    """Convert PRMS CSV files to netcdf.
+
+    Args:
+        var_name: str name of the variable to create
+        data_dir: where the csv file is found and the netcdf file will be
+            written (could add argument to output to a differnt dir)
+    """
+    if output_dir is None:
+        output_dir = data_dir
+    csv_path = data_dir / f"{var_name}.csv"
+    nc_path = output_dir / f"{var_name}.nc"
+    CsvFile(csv_path).to_netcdf(nc_path)
+    assert nc_path.exists()
+
+
+def convert_soltab_to_nc(
+    soltab_file: pl.Path, output_dir: pl.Path, domain_dir: pl.Path
+):
+    """Convert soltab files to NetCDF, one file for each variable
+
+    The inputs files are "soltab_debug" specifically written by the PRMS5.2.1
+    in the pywatershed repository.
+
+    Args:
+        soltab_file: the pl.Path of the  soltab_debug file.
+        output_dir: pl.Path where the netcdf file output will be written
+        domain_dir: defaults to the parent dir of soltab_file, the pl.Path
+            where domain files (parameter & control) for the domain can be
+            found
+    """
+    if domain_dir is None:
+        domain_dir = soltab_file.parent
+
+    params = pws.parameters.PrmsParameters.load(domain_dir / "myparam.param")
+    nhm_ids = params.parameters["nhm_id"]
+
+    soltab = Soltab(soltab_file, nhm_ids=nhm_ids)
+    soltab.to_netcdf(output_dir=output_dir)
+    print(f"Created NetCDF files from soltab file {soltab_file}:")
+
+    for var in soltab.variables:
+        nc_path = output_dir / f"{var}.nc"
+        assert nc_path.exists()
diff --git a/test_data/generate/prms_diagnostic_variables.py b/test_data/generate/prms_diagnostic_variables.py
new file mode 100644
index 00000000..12610eb3
--- /dev/null
+++ b/test_data/generate/prms_diagnostic_variables.py
@@ -0,0 +1,332 @@
+import pathlib as pl
+
+import pywatershed as pws
+from pywatershed.constants import epsilon32, nan, zero
+
+import numpy as np
+import xarray as xr
+
+"""This module is for generating PRMS diagnostic variables"""
+
+# For generating timeseries of previous states
+previous_vars = {
+    "gwres_stor": pws.PRMSGroundwater,
+    "dprst_stor_hru": pws.PRMSRunoff,
+    "hru_impervstor": pws.PRMSRunoff,
+    "pref_flow_stor": pws.PRMSSoilzone,
+    "slow_stor": pws.PRMSSoilzone,
+    "soil_lower": pws.PRMSSoilzone,
+    "soil_moist": pws.PRMSSoilzone,
+    "soil_rechr": pws.PRMSSoilzone,
+    "ssres_stor": pws.PRMSSoilzone,
+    "freeh2o": pws.PRMSSnow,
+    "pk_ice": pws.PRMSSnow,
+}
+
+prev_rename = {
+    "gwres_stor": "gwres_stor_old",
+}
+
+change_rename = {}
+
+
+def diagnose_simple_vars_to_nc(
+    var_name: str,
+    data_dir: pl.Path,
+    domain_dir: pl.Path = None,
+    output_dir: pl.Path = None,
+):
+    """Diagnose variables from PRMS output with only dependencies on var_name
+
+    These variables are "simple" because they do not have any dependencies
+    besides the "var_name" variable which is a PRMS output variable.
+
+    This is not to say that PRMS dosent output diagnostic variables, it
+    just dosent output these natively.
+
+    Args:
+        var_name: str name of the variable to create
+        data_dir: where the netcdf file will be written, also where to look
+            for necessary inputs to create the diagnostic variable.
+        domain_dir: defaults to the parent dir of output_dir, this is where
+            domain files (parameter & control) for the domain can be found
+        output_dir: defaults to data_dir unless otherwise specified
+
+    """
+
+    if domain_dir is None:
+        domain_dir = data_dir.parent
+
+    if output_dir is None:
+        output_dir = data_dir
+
+    nc_path = data_dir / f"{var_name}.nc"
+
+    if var_name in previous_vars.keys():
+        # the _prev is the desired suffix but PRMS legacy is inconsistent
+        current = xr.open_dataarray(nc_path)
+        prev = current * nan
+        prev[:] = np.roll(current.values, 1, axis=0)
+
+        # the final value (-1) was wrapped to the zeroth position
+        # get the initial conditions for the first time by initializing the
+        # model This works based on the control file, so could handle restart.
+        control = pws.Control.load(domain_dir / "control.test")
+        control.options = control.options | {
+            "input_dir": domain_dir / "output",
+        }
+        params = pws.parameters.PrmsParameters.load(
+            domain_dir / "myparam.param"
+        )
+        proc_class = previous_vars[var_name]
+        inputs = {}
+        for input_name in proc_class.get_inputs():
+            # taken from process._initialize_var
+            meta = pws.meta.find_variables([input_name])
+            dims = [params.dims[vv] for vv in meta[input_name]["dims"]]
+            init_type = meta[input_name]["type"]
+
+            if len(dims) == 1:
+                inputs[input_name] = np.full(dims, nan, dtype=init_type)
+            else:
+                inputs[input_name] = pws.base.timeseries.TimeseriesArray(
+                    var_name=input_name,
+                    control=control,
+                    array=np.full(dims, nan, dtype=init_type),
+                    time=np.arange(
+                        np.datetime64("2017-01-01"),
+                        np.datetime64("2017-01-03"),
+                    ),
+                )
+
+        proc = proc_class(
+            control=control,
+            discretization=None,
+            parameters=params,
+            **inputs,
+        )
+        # these are the init_values
+        prev[0, :] = proc[var_name]
+        del proc
+
+        change = current - prev
+
+        # write the previous file
+        if var_name in prev_rename.keys():
+            nc_name = f"{prev_rename[var_name]}.nc"
+        else:
+            nc_name = f"{var_name}_prev.nc"
+        out_nc_path = output_dir / nc_name
+        prev.rename(out_nc_path.stem).to_netcdf(out_nc_path)
+        assert nc_path.exists()
+
+        # write the change file
+        if var_name in change_rename.keys():
+            nc_name = f"{change_rename[var_name]}.nc"
+        else:
+            nc_name = f"{var_name}_change.nc"
+        out_nc_path = output_dir / nc_name
+        change.rename(out_nc_path.stem).to_netcdf(out_nc_path)
+        assert nc_path.exists()
+        return nc_path
+
+    if var_name in ["sroff", "ssres_flow", "gwres_flow"]:
+        params = pws.parameters.PrmsParameters.load(
+            domain_dir / "myparam.param"
+        )
+        ds = xr.open_dataset(nc_path)
+        ds = ds.rename({var_name: f"{var_name}_vol"})
+        ds = ds * params.data_vars["hru_in_to_cf"]
+        ds.to_netcdf(output_dir / f"{var_name}_vol.nc")
+        ds.close()
+
+    if var_name == "infil":
+        params = pws.parameters.PrmsParameters.load(
+            domain_dir / "myparam.param"
+        ).parameters
+        imperv_frac = params["hru_percent_imperv"]
+        dprst_frac = params["dprst_frac"]
+        perv_frac = 1.0 - imperv_frac - dprst_frac
+        ds = xr.open_dataset(nc_path.with_suffix(".nc"))
+        ds = ds.rename(infil="infil_hru")
+        ds["infil_hru"] = ds["infil_hru"] * perv_frac
+        ds.to_netcdf(output_dir / "infil_hru.nc")
+        ds.close()
+
+
+def diagnose_final_vars_to_nc(
+    var_name: str,
+    data_dir: pl.Path,
+    domain_dir: pl.Path = None,
+    output_dir: pl.Path = None,
+):
+    """Diagnose variables from multiple PRMS outputs
+
+    In this case "var_name" is not a PRMS output variable, the variable is a
+    new variable which requires much/all of the PRMS output to be already
+    present in netcdf output format. This is why this is a final diagnostic
+
+    Currently only: "through_rain"
+
+    Args:
+        var_name: str name of the variable to create, not a PRMS variable.
+        output_dir: where the netcdf file will be written, also where to look
+            for necessary inputs to create the diagnostic variable.
+        domain_dir: defaults to the parent dir of output_dir, this is where
+            domain files (parameter & control) for the domain can be found
+
+    """
+    if domain_dir is None:
+        domain_dir = data_dir.parent
+
+    if output_dir is None:
+        output_dir = data_dir
+
+    if var_name == "through_rain":
+        out_file = output_dir / "through_rain.nc"
+
+        data_vars = [
+            "net_ppt",
+            "pptmix_nopack",
+            "snowmelt",
+            "pkwater_equiv",
+            "snow_evap",
+            "net_snow",
+            "net_rain",
+        ]
+
+        data = {}
+        for vv in data_vars:
+            data_file = data_dir / f"{vv}.nc"
+            data[vv] = xr.open_dataarray(data_file)
+
+        nearzero = 1.0e-6
+
+        cond1 = data["net_ppt"] > zero
+        cond2 = data["pptmix_nopack"] != 0
+        cond3 = data["snowmelt"] < nearzero
+        cond4 = data["pkwater_equiv"] < epsilon32
+        cond5 = data["snow_evap"] < nearzero
+        cond6 = data["net_snow"] < nearzero
+
+        through_rain = data["net_rain"] * zero
+        # these are in reverse order
+        through_rain[:] = np.where(
+            cond1 & cond3 & cond4 & cond6, data["net_rain"], zero
+        )
+        through_rain[:] = np.where(
+            cond1 & cond3 & cond4 & cond5, data["net_ppt"], through_rain
+        )
+        through_rain[:] = np.where(
+            cond1 & cond2, data["net_rain"], through_rain
+        )
+
+        through_rain.to_dataset(name="through_rain").to_netcdf(out_file)
+        through_rain.close()
+
+        for vv in data_vars:
+            data[vv].close()
+
+        assert out_file.exists()
+
+    if var_name in [
+        "channel_sroff_vol",
+        "channel_ssres_flow_vol",
+        "channel_gwres_flow_vol",
+        "seg_lateral_inflow",
+    ]:
+        if var_name != "seg_lateral_inflow":
+            return
+
+        data_vars = [
+            "sroff_vol",
+            "ssres_flow_vol",
+            "gwres_flow_vol",
+            "seg_outflow",
+            "seg_inflow",
+        ]
+        data = {}
+        for vv in data_vars:
+            data_file = data_dir / f"{vv}.nc"
+            data[vv] = xr.open_dataarray(data_file)
+
+        control = pws.Control.load(domain_dir / "control.test")
+        s_per_time = control.time_step_seconds
+        params = pws.parameters.PrmsParameters.load(
+            domain_dir / "myparam.param"
+        )
+
+        ntimes = control.n_times
+        nseg = params.dims["nsegment"]
+        nhru = params.dims["nhru"]
+        hru_segment = params.parameters["hru_segment"] - 1
+        to_segment = (params.parameters["tosegment"] - 1).astype("int64")
+
+        outflow_mask = np.full((len(to_segment)), False)
+        for iseg in range(nseg):
+            if to_segment[iseg] < 0:
+                outflow_mask[iseg] = True
+                continue
+
+        seg_stor_change = (
+            data["seg_inflow"] - data["seg_outflow"]
+        ) * s_per_time
+
+        channel_outflow_vol = (
+            np.where(outflow_mask, data["seg_outflow"], zero)
+        ) * s_per_time
+
+        seg_lateral_inflow = np.zeros((ntimes, nseg))
+        channel_sroff_vol = np.zeros((ntimes, nhru))
+        channel_ssres_flow_vol = np.zeros((ntimes, nhru))
+        channel_gwres_flow_vol = np.zeros((ntimes, nhru))
+
+        for ihru in range(nhru):
+            iseg = hru_segment[ihru]
+            if iseg < 0:
+                # This is bad, selective handling of fluxes is not cool,
+                # mass is being discarded in a way that has to be coordinated
+                # with other parts of the code.
+                # This code shuold be removed evenutally.
+                continue
+
+            else:
+                channel_sroff_vol[:, ihru] = data["sroff_vol"].values[:, ihru]
+                channel_ssres_flow_vol[:, ihru] = data[
+                    "ssres_flow_vol"
+                ].values[:, ihru]
+                channel_gwres_flow_vol[:, ihru] = data[
+                    "gwres_flow_vol"
+                ].values[:, ihru]
+
+            # cubicfeet to cfs
+            lateral_inflow = (
+                channel_sroff_vol[:, ihru]
+                + channel_ssres_flow_vol[:, ihru]
+                + channel_gwres_flow_vol[:, ihru]
+            ) / (s_per_time)
+
+            seg_lateral_inflow[:, iseg] += lateral_inflow
+
+        for hvar in [
+            "channel_sroff_vol",
+            "channel_ssres_flow_vol",
+            "channel_gwres_flow_vol",
+            "seg_lateral_inflow",
+            "seg_stor_change",
+            "channel_outflow_vol",
+        ]:
+            if hvar in [
+                "seg_lateral_inflow",
+                "seg_stor_change",
+                "channel_outflow_vol",
+            ]:
+                dum = data["seg_outflow"].rename(hvar) * np.nan
+            else:
+                dum = data["sroff_vol"].rename(hvar) * np.nan
+
+            dum[:, :] = locals()[hvar]
+            out_file = data_dir / f"{hvar}.nc"
+            dum.to_netcdf(out_file)
+            print(out_file)
+            assert out_file.exists()
diff --git a/test_data/scripts/pytest.ini b/test_data/generate/pytest.ini
similarity index 89%
rename from test_data/scripts/pytest.ini
rename to test_data/generate/pytest.ini
index 51b9641a..b99a9aef 100644
--- a/test_data/scripts/pytest.ini
+++ b/test_data/generate/pytest.ini
@@ -3,6 +3,5 @@ addopts = --order-dependencies
 python_files =
     test_*.py
 python_functions =
-    create_*
     make_*
     test_*
\ No newline at end of file
diff --git a/test_data/generate/remove_prms_csvs.py b/test_data/generate/remove_prms_csvs.py
new file mode 100644
index 00000000..86a139e9
--- /dev/null
+++ b/test_data/generate/remove_prms_csvs.py
@@ -0,0 +1,6 @@
+csvs_to_keep = ["hru_ppt", "intcp_stor", "potet", "gwres_stor"]
+
+
+def test_csv_to_netcdf(csv_file):
+    if csv_file.with_suffix("").name not in (csvs_to_keep):
+        csv_file.unlink()
diff --git a/test_data/scripts/test_run_domains.py b/test_data/generate/run_prms_domains.py
similarity index 100%
rename from test_data/scripts/test_run_domains.py
rename to test_data/generate/run_prms_domains.py
diff --git a/test_data/scripts/submit_domains_tg.sh b/test_data/generate/submit_domains_tg.sh
similarity index 100%
rename from test_data/scripts/submit_domains_tg.sh
rename to test_data/generate/submit_domains_tg.sh
diff --git a/test_data/hru_1/control.test b/test_data/hru_1/control.test
index e692a8b6..74ea87b3 100644
--- a/test_data/hru_1/control.test
+++ b/test_data/hru_1/control.test
@@ -238,7 +238,7 @@ nhruOutON_OFF
 nhruOutVars
 1
 1
-76
+81
 ####
 nhruOut_freq
 1
@@ -251,12 +251,13 @@ nhruOut_format
 5
 ####
 nhruOutVar_names
-76
+81
 4
 albedo
 cap_infil_tot
 contrib_fraction
 dprst_evap_hru
+dprst_insroff_hru
 dprst_seep_hru
 dprst_stor_hru
 dprst_sroff_hru
@@ -288,7 +289,11 @@ net_snow
 newsnow
 orad_hru
 perv_actet
+pk_def
+pk_den
+pk_depth
 pk_ice
+pk_temp
 pkwater_ante
 pkwater_equiv
 potet
@@ -343,7 +348,7 @@ nsegmentOutON_OFF
 nsegmentOutVars
 1
 1
-9
+10
 ####
 nsegmentOut_freq
 1
@@ -351,8 +356,9 @@ nsegmentOut_freq
 1
 ####
 nsegmentOutVar_names
-9
+10
 4
+seg_inflow
 seg_outflow
 seg_gwflow
 seg_sroff
diff --git a/test_data/scripts/test_nc_domains.py b/test_data/scripts/test_nc_domains.py
deleted file mode 100644
index 38cf3686..00000000
--- a/test_data/scripts/test_nc_domains.py
+++ /dev/null
@@ -1,165 +0,0 @@
-from pathlib import Path
-from filelock import FileLock
-
-import numpy as np
-
-from pywatershed import CsvFile, Soltab
-from pywatershed.parameters import PrmsParameters
-from pywatershed.constants import epsilon64, zero
-
-import pytest
-import xarray as xr
-
-
-# For generating timeseries of previous states
-previous_vars = [
-    "dprst_stor_hru",
-    "freeh2o",
-    "hru_impervstor",
-    "pk_ice",
-    "pref_flow_stor",
-    "slow_stor",
-    "soil_lower",
-    "soil_moist",
-    "soil_rechr",
-    "ssres_stor",
-]
-
-
-@pytest.fixture
-def netcdf_file(csv_file) -> Path:
-    """Convert CSV files from model output to NetCDF"""
-
-    nc_path = csv_file.with_suffix(".nc")
-    data_dir = csv_file.parent
-    domain_dir = data_dir.parent
-    CsvFile(csv_file).to_netcdf(nc_path)
-    assert nc_path.exists()
-
-    if csv_file.stem in previous_vars:
-        nc_name = csv_file.stem + "_prev.nc"
-        nc_path = data_dir / nc_name
-        csv = CsvFile({nc_path.stem: csv_file})
-        csv._get_data()  # why so private?
-        orig_dates = csv._data["date"].copy()
-        csv._data = np.roll(csv._data, 1, axis=0)
-        csv._data["date"] = orig_dates
-        # Here we will eventually want to supply the desired initial conditions
-        for hh in range(len(csv._data[0])):
-            if hh == 0:
-                continue
-            csv._data[0][hh] = np.zeros([1])[0]
-        csv.to_netcdf(nc_path)
-        assert nc_path.exists()
-
-    if nc_path.stem in ["sroff", "ssres_flow", "gwres_flow"]:
-        params = PrmsParameters.load(domain_dir / "myparam.param")
-        var = nc_path.stem
-        ds = xr.open_dataset(nc_path)
-        ds = ds.rename({var: f"{var}_vol"})
-        ds = ds * params.data_vars["hru_in_to_cf"]
-        ds.to_netcdf(data_dir / f"{var}_vol.nc")
-        ds.close()
-
-    if nc_path.stem == "infil":
-        params = PrmsParameters.load(domain_dir / "myparam.param").parameters
-        imperv_frac = params["hru_percent_imperv"]
-        dprst_frac = params["dprst_frac"]
-        perv_frac = 1.0 - imperv_frac - dprst_frac
-        ds = xr.open_dataset(nc_path.with_suffix(".nc"))
-        ds = ds.rename(infil="infil_hru")
-        ds["infil_hru"] = ds["infil_hru"] * perv_frac
-        ds.to_netcdf(data_dir / "infil_hru.nc")
-        ds.close()
-
-    return nc_path
-
-
-def make_netcdf_files(netcdf_file):
-    print(f"Created NetCDF from CSV: {netcdf_file}")
-
-
-@pytest.fixture(scope="session")
-def soltab_netcdf_file(tmp_path_factory, soltab_file) -> Path:
-    """Convert soltab files to NetCDF, one file for each variable"""
-
-    # the nhm_ids are not available in the solta_debug file currently, so get
-    # them from the domain parameters
-    domain_dir = soltab_file.parent
-    output_dir = domain_dir / "output"
-    output_dir.mkdir(exist_ok=True)
-    root_tmpdir = tmp_path_factory.getbasetemp().parent
-    nc_paths = []
-
-    with FileLock(root_tmpdir / "soltab_nc.lock"):
-        yield  # postpone the work until session cleanup
-
-        # this is a hack that should probably rely on the yaml if/when this
-        # fails
-        params = PrmsParameters.load(domain_dir / "myparam.param")
-        nhm_ids = params.parameters["nhm_id"]
-        soltab = Soltab(soltab_file, nhm_ids=nhm_ids)
-
-        already_created = (output_dir / "soltab_potsw.nc").is_file()
-        if not already_created:
-            soltab.to_netcdf(output_dir=output_dir)
-            print(f"Created NetCDF files from soltab file {soltab_file}:")
-
-        for var in soltab.variables:
-            nc_path = output_dir / f"{var}.nc"
-            nc_paths.append(nc_path)
-            assert nc_path.exists()
-            print(f"\t{nc_path}")
-
-
-def make_soltab_netcdf_files(soltab_netcdf_file):
-    print(f"Creating NetCDF files for soltab file {soltab_netcdf_file}")
-
-
-@pytest.fixture(scope="session")
-def final_netcdf_file(tmp_path_factory, simulation) -> Path:
-    """Create the final NetCDF file (through_rain.nc) from other NetCDFs"""
-
-    root_tmpdir = tmp_path_factory.getbasetemp().parent
-    data_dir = Path(simulation["ws"]) / "output"
-    data_dir.mkdir(exist_ok=True)
-    nc_path = data_dir / "through_rain.nc"
-
-    with FileLock(root_tmpdir / "final_nc.lock"):
-        yield  # do this in session cleanup
-
-        if nc_path.is_file():
-            return
-
-        data_vars = [
-            "net_rain",
-            "pk_ice_prev",
-            "freeh2o_prev",
-            "newsnow",
-            "pptmix_nopack",
-        ]
-
-        data = {}
-        for vv in data_vars:
-            data_file = data_dir / f"{vv}.nc"
-            result = xr.open_dataset(data_file)[vv]
-            data[vv] = result
-
-        try:
-            wh_through = (
-                ((data["pk_ice_prev"] + data["freeh2o_prev"]) <= epsilon64)
-                & ~(data["newsnow"] == 1)
-            ) | (data["pptmix_nopack"] == 1)
-
-            through_rain = data["net_rain"].copy()
-            through_rain[:] = np.where(wh_through, data["net_rain"], zero)
-            through_rain.to_dataset(name="through_rain").to_netcdf(nc_path)
-            through_rain.close()
-            assert nc_path.exists()
-        finally:
-            for vv in data_vars:
-                data[vv].close()
-
-
-def make_final_netcdf_files(final_netcdf_file):
-    print(f"Creating final NetCDF file {final_netcdf_file}")
diff --git a/test_data/ucb_2yr/control.test b/test_data/ucb_2yr/control.test
index df624841..effd9667 100644
--- a/test_data/ucb_2yr/control.test
+++ b/test_data/ucb_2yr/control.test
@@ -238,7 +238,7 @@ nhruOutON_OFF
 nhruOutVars
 1
 1
-78
+83
 ####
 nhruOut_freq
 1
@@ -251,13 +251,14 @@ nhruOut_format
 4
 ####
 nhruOutVar_names
-78
+83
 4
 albedo
 cap_infil_tot
 contrib_fraction
 dprst_seep_hru
 dprst_evap_hru
+dprst_insroff_hru
 dprst_sroff_hru
 dprst_stor_hru
 freeh2o
@@ -288,7 +289,11 @@ net_snow
 newsnow
 orad_hru
 perv_actet
+pk_def
+pk_den
+pk_depth
 pk_ice
+pk_temp
 pkwater_ante
 pkwater_equiv
 potet
@@ -345,7 +350,7 @@ nsegmentOutON_OFF
 nsegmentOutVars
 1
 1
-9
+10
 ####
 nsegmentOut_freq
 1
@@ -353,8 +358,9 @@ nsegmentOut_freq
 1
 ####
 nsegmentOutVar_names
-9
+10
 4
+seg_inflow
 seg_outflow
 seg_gwflow
 seg_sroff