diff --git a/.github/workflows/pydep.yml b/.github/workflows/pydep.yml
index 3d1ede4d..9b577a8d 100644
--- a/.github/workflows/pydep.yml
+++ b/.github/workflows/pydep.yml
@@ -16,7 +16,7 @@ jobs:
     runs-on: ubuntu-latest
     strategy:
       matrix:
-        PYTHON_VERSION: ["3.9", "3.10", "3.11"]
+        PYTHON_VERSION: ["3.9", "3.11", "3.12"]
     env:
       PYTHON_VERSION: ${{ matrix.PYTHON_VERSION }}
     steps:
@@ -26,6 +26,12 @@ jobs:
       run: |
         cat .github/workflows/etchosts.txt | sudo tee -a /etc/hosts
 
+    # daryl is a hack
+    - name: Create symlink
+      run: |
+        set -e
+        sudo ln -s `pwd` /opt/dep
+
     # setup conda-forge with micromamba
     - name: Setup Python
       uses: mamba-org/setup-micromamba@v1
@@ -56,12 +62,13 @@ jobs:
       run: |
         set -e
         python -m pip install . --upgrade --no-deps
-        coverage run --source=pydep setup.py test
-        coverage xml
+        python -m pytest --cov=pydep tests
+        python -m coverage xml
+
 
     - name: Upload codecov
-      if: ${{ env.PYTHON_VERSION == '3.11' }}
-      uses: codecov/codecov-action@v3
+      if: ${{ env.PYTHON_VERSION == '3.12' }}
+      uses: codecov/codecov-action@v4
       with:
         token: ${{ secrets.CODECOV_TOKEN }}
         files: coverage.xml
diff --git a/.gitignore b/.gitignore
index dc5b4025..d484b0ab 100644
--- a/.gitignore
+++ b/.gitignore
@@ -11,6 +11,7 @@
 **/*.xlsx
 **/*.txt
 **/*.pptx
+**/*.feather
 /.vscode
 prj2wepp/wepp/hill.slp
 prj2wepp/wepp/data/managements
diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml
index 83903a88..7fc4355c 100644
--- a/.pre-commit-config.yaml
+++ b/.pre-commit-config.yaml
@@ -1,6 +1,6 @@
 repos:
   - repo: https://github.com/astral-sh/ruff-pre-commit
-    rev: "v0.3.4"
+    rev: "v0.3.5"
     hooks:
       - id: ruff
         args: [--fix, --exit-non-zero-on-fix]
diff --git a/pyproject.toml b/pyproject.toml
index 5bdb6380..882e1aea 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -1,5 +1,11 @@
 [tool.ruff]
 line-length = 79
-lint.select = ["E", "F", "I"]
+lint.select = [
+    "B",  # bugbear
+    "E",
+    "F",
+    "I",
+    # "PD",  # pandas
+]
 target-version = "py39"
 
diff --git a/scripts/cligen/plot_qc_mrms_precip.py b/scripts/cligen/plot_qc_mrms_precip.py
index 36a9be7f..068a3934 100644
--- a/scripts/cligen/plot_qc_mrms_precip.py
+++ b/scripts/cligen/plot_qc_mrms_precip.py
@@ -10,6 +10,7 @@
 import numpy as np
 from geopandas import read_postgis
 from matplotlib.patches import Polygon
+from pyiem.database import get_dbconn
 from pyiem.iemre import (
     EAST,
     NORTH,
@@ -19,7 +20,6 @@
     get_daily_mrms_ncname,
 )
 from pyiem.plot import MapPlot
-from pyiem.util import get_dbconn
 
 YS = np.arange(SOUTH, NORTH, 0.01)
 XS = np.arange(WEST, EAST, 0.01)
@@ -92,7 +92,7 @@ def do(valid):
         geom_col="geom",
         index_col="fpath",
     )
-    for fpath, row in df.iterrows():
+    for _fpath, row in df.iterrows():
         mp.ax.plot(row["geom"].xy[0], row["geom"].xy[1], c="k")
 
     df = read_postgis(
diff --git a/scripts/convergence/plot_avgs.py b/scripts/convergence/plot_avgs.py
index 91d681a2..ba77dd69 100644
--- a/scripts/convergence/plot_avgs.py
+++ b/scripts/convergence/plot_avgs.py
@@ -5,7 +5,7 @@
 import matplotlib.pyplot as plt
 import numpy as np
 import pandas as pd
-from pyiem.util import get_dbconn
+from pyiem.database import get_dbconn
 from scipy import stats
 
 YEAR = int(sys.argv[1])
@@ -62,7 +62,7 @@
     (fig, ax) = plt.subplots(1, 1, figsize=(8, 6))
 
     averages = []
-    for i in range(10):
+    for _ in range(10):
         averages.append([])
 
     for i in range(10):
diff --git a/scripts/dynamic_tillage/accum_acres_plot.py b/scripts/dynamic_tillage/accum_acres_plot.py
index 44649315..c9836d24 100644
--- a/scripts/dynamic_tillage/accum_acres_plot.py
+++ b/scripts/dynamic_tillage/accum_acres_plot.py
@@ -1,8 +1,116 @@
 """Plot accumulated acres diagnostic."""
 
 import pandas as pd
-from pyiem.plot import figure_axes
-from pyiem.util import get_sqlalchemy_conn
+from matplotlib.patches import Rectangle
+from pyiem.database import get_sqlalchemy_conn
+from pyiem.plot import MapPlot, figure_axes
+
+
+def plot_map(i, dt, huc12df, fields):
+    """Make a map diagnostic"""
+    minx, miny, maxx, maxy = huc12df["geom"].to_crs(4326).total_bounds
+    buffer = 0.01
+    huc12 = huc12df.index.values[0]
+
+    mp = MapPlot(
+        title=f"DEP Planting Progress {dt:%Y %b %d} for HUC12: {huc12}",
+        logo="dep",
+        sector="custom",
+        west=minx - buffer,
+        north=maxy + buffer,
+        south=miny - buffer,
+        east=maxx + buffer,
+        caption="Daily Erosion Project",
+        continentalcolor="white",
+    )
+    huc12df.to_crs(mp.panels[0].crs).plot(
+        ax=mp.panels[0].ax,
+        aspect=None,
+        facecolor="None",
+        edgecolor="b",
+        linewidth=2,
+        zorder=100,
+    )
+    fields["color"] = "white"
+    fields.loc[fields["till1"].notna(), "color"] = "tan"
+    fields.loc[fields["plant"].notna(), "color"] = "g"
+    fields.to_crs(mp.panels[0].crs).plot(
+        ax=mp.panels[0].ax,
+        aspect=None,
+        facecolor=fields["color"],
+        edgecolor="k",
+        linewidth=1,
+        zorder=100,
+    )
+    p0 = Rectangle((0, 0), 1, 1, fc="white", ec="k")
+    p1 = Rectangle((0, 0), 1, 1, fc="tan", ec="k")
+    p2 = Rectangle((0, 0), 1, 1, fc="g", ec="k")
+    mp.panels[0].ax.legend(
+        (p0, p1, p2),
+        ("Awaiting", "Tilled", "Planted"),
+        ncol=3,
+        fontsize=11,
+        loc=2,
+    )
+
+    mp.fig.savefig(f"{i:04.0f}.png")
+    mp.close()
+
+
+def plot_timeseries(year, df, huc12):
+    """Make a diagnostic."""
+    (fig, ax) = figure_axes(
+        logo="dep",
+        title=f"DEP {year} Tillage/Plant Operation Timing for HUC12: {huc12}",
+        subtitle="<=10% Daily Rate, All Field OFEs below 0.9 Plastic Limit",
+    )
+    ax2 = ax.twinx()
+    x = pd.date_range(f"{year}/04/15", f"{year}/06/03")
+    ax2.bar(
+        x - pd.Timedelta(hours=8),
+        df["acres_planted"],
+        color="#00ff00",
+        alpha=0.5,
+        width=0.3,
+        align="edge",
+        zorder=2,
+    )
+    ax2.bar(
+        x,
+        df["acres_tilled"],
+        color="#ff0000",
+        alpha=0.5,
+        width=0.3,
+        align="edge",
+        zorder=2,
+    )
+
+    ax.plot(
+        x,
+        df["acres_not_planted"],
+        color="k",
+        label="Not Planted",
+        zorder=3,
+        lw=3,
+    )
+    ax.plot(x, df["acres_to_till"], color="r", label="Tillage", zorder=3, lw=3)
+    ax.plot(
+        x,
+        df["acres_to_plant"],
+        c="g",
+        label="Plant",
+        zorder=3,
+        lw=3,
+    )
+    ax.set_ylabel("Acres Available for Operations")
+    ax.legend(loc=1)
+    ax.grid(True)
+    ax.set_ylim(bottom=-10)
+    ax2.set_ylim(bottom=-10)
+    ax2.set_ylabel("Acreas Worked per Day (bars)")
+    ax.set_zorder(ax2.get_zorder() + 1)
+    ax.patch.set_visible(False)
+    fig.savefig(f"{year}_timeseries.png")
 
 
 def main():
diff --git a/scripts/dynamic_tillage/bootstrap_year.py b/scripts/dynamic_tillage/bootstrap_year.py
new file mode 100644
index 00000000..da7a74c8
--- /dev/null
+++ b/scripts/dynamic_tillage/bootstrap_year.py
@@ -0,0 +1,111 @@
+"""Bootstrap database for dynamic tillage."""
+
+from datetime import date
+
+import click
+from pandas.io.sql import read_sql
+from pydep.tillage import make_tillage
+from pyiem.database import get_dbconnc, get_sqlalchemy_conn
+from pyiem.util import logger
+from sqlalchemy import text
+
+LOG = logger()
+
+
+@click.command()
+@click.option("--year", type=int, help="Year to bootstrap", required=True)
+@click.option("--scenario", type=int, default=0)
+def main(year, scenario):
+    """Run for the given year."""
+    with get_sqlalchemy_conn("idep") as conn:
+        # Delete current year's data
+        res = conn.execute(
+            text("""
+                 delete from field_operations o USING fields f
+                 WHERE o.field_id = f.field_id and o.year = :year
+                 and f.scenario = :scenario
+                 """),
+            {"year": year, "scenario": scenario},
+        )
+        LOG.info("deleted %s field operation rows", res.rowcount)
+        conn.commit()
+        # figure out the fields we need to compute
+        fields = read_sql(
+            text("""
+                select field_id, landuse, management,
+                st_y(st_centroid(st_transform(geom, 4326))) as lat
+                from fields where scenario = :scenario
+                and isag and landuse is not null and management is not null
+            """),
+            conn,
+            params={"scenario": scenario},
+            index_col="field_id",
+        )
+    LOG.info("Processing %s fields", len(fields.index))
+    conn, cursor = get_dbconnc("idep")
+    colidx = year - 2007
+    inserts = 0
+    for field_id, row in fields.iterrows():
+        zone = "KS_NORTH"
+        if row["lat"] >= 42.5:
+            zone = "IA_NORTH"
+        elif row["lat"] >= 41.5:
+            zone = "IA_CENTRAL"
+        elif row["lat"] >= 40.5:
+            zone = "IA_SOUTH"
+        res = make_tillage(
+            scenario,
+            zone,
+            row["landuse"][colidx - 1],
+            row["landuse"][colidx],
+            row["landuse"][colidx - 2],  # HACK
+            row["management"][colidx],
+            year - 2006,
+        )
+        events = 0
+        for line in res.split("\n"):
+            tokens = line.split()
+            if len(tokens) < 5 or tokens[4] != "Tillage" or int(tokens[0]) > 6:
+                continue
+            events += 1
+        # till1, till2, till3, plant
+        dates = [None, None, None, None]
+        if events == 1:
+            dates[3] = date(year, 6, 1)
+        elif events == 2:
+            dates[0] = date(year, 6, 1)
+            dates[3] = date(year, 6, 2)
+        elif events == 3:
+            dates[0] = date(year, 6, 1)
+            dates[1] = date(year, 6, 2)
+            dates[3] = date(year, 6, 3)
+        elif events == 4:
+            dates[0] = date(year, 6, 1)
+            dates[1] = date(year, 6, 2)
+            dates[2] = date(year, 6, 3)
+            dates[3] = date(year, 6, 4)
+        elif events > 4:
+            LOG.info("events: %s", events)
+            LOG.info("res: %s", res)
+            raise ValueError("too many events")
+        cursor.execute(
+            """
+            INSERT into field_operations
+            (field_id, year, till1, till2, till3, plant)
+            VALUES (%s, %s, %s, %s, %s, %s)
+        """,
+            (field_id, year, dates[0], dates[1], dates[2], dates[3]),
+        )
+        inserts += 1
+        if inserts % 10_000 == 0:
+            cursor.close()
+            conn.commit()
+            cursor = conn.cursor()
+
+    cursor.close()
+    conn.commit()
+    conn.close()
+
+
+if __name__ == "__main__":
+    main()
diff --git a/scripts/dynamic_tillage/compute_smstate.py b/scripts/dynamic_tillage/compute_smstate.py
index a47e0946..887c04d0 100644
--- a/scripts/dynamic_tillage/compute_smstate.py
+++ b/scripts/dynamic_tillage/compute_smstate.py
@@ -6,21 +6,22 @@
 """
 
 import os
-import sys
 from multiprocessing import cpu_count
 from multiprocessing.pool import ThreadPool
 
+import click
 import geopandas as gpd
 import pandas as pd
 from pydep.io.dep import read_wb
-from pyiem.util import get_sqlalchemy_conn, logger
+from pyiem.database import get_sqlalchemy_conn
+from pyiem.util import logger
 from sqlalchemy import text
 from tqdm import tqdm
 
 LOG = logger()
 
 
-def job(huc12, year):
+def job(huc12, dt):
     """Do work for a HUC12"""
     with get_sqlalchemy_conn("idep") as conn:
         # Build up cross reference of fields and flowpath/OFEs
@@ -38,8 +39,6 @@ def job(huc12, year):
             index_col=None,
         )
     dfs = []
-    sts = pd.Timestamp(f"{year}/04/14")  # Need "yesterday"
-    ets = pd.Timestamp(f"{year}/06/01")
     for flowpath in df["fpath"].unique():
         flowpath = int(flowpath)
         wbfn = f"/i/0/wb/{huc12[:8]}/{huc12[8:]}/{huc12}_{flowpath}.wb"
@@ -49,7 +48,7 @@ def job(huc12, year):
             smdf = read_wb(wbfn)
         except Exception as exp:
             raise ValueError(f"Read {wbfn} failed") from exp
-        smdf = smdf[(smdf["date"] >= sts) & (smdf["date"] <= ets)]
+        smdf = smdf[smdf["date"] == pd.Timestamp(dt)]
         # Each flowpath + ofe should be associated with a fbndid
         for ofe, gdf in smdf.groupby("ofe"):
             fbndid = df[(df["fpath"] == flowpath) & (df["ofe"] == ofe)].iloc[
@@ -61,20 +60,26 @@ def job(huc12, year):
         df = pd.concat(dfs)
         df["huc12"] = huc12
         return df
+    return None
 
 
-def main(argv):
+@click.command()
+@click.option("--date", "dt", type=click.DateTime(), help="Date to process")
+@click.option("--huc12", type=str, help="HUC12 to process")
+def main(dt, huc12):
     """Go Main Go."""
-    year = int(argv[1])
+    dt = dt.date()
+    huc12limiter = "" if huc12 is None else " and huc_12 = :huc12"
     with get_sqlalchemy_conn("idep") as conn:
         huc12df = gpd.read_postgis(
-            """
+            text(f"""
             SELECT geom, huc_12,
             ST_x(st_transform(st_centroid(geom), 4326)) as lon,
             ST_y(st_transform(st_centroid(geom), 4326)) as lat
-            from huc12 where scenario = 0
-            """,
+            from huc12 where scenario = 0 {huc12limiter}
+            """),
             conn,
+            params={"huc12": huc12},
             geom_col="geom",
             index_col="huc_12",
         )
@@ -94,16 +99,16 @@ def _eb(error):
             print(error)
             pool.terminate()
 
-        for huc12 in huc12df.index.values:
+        for _huc12 in huc12df.index.values:
             pool.apply_async(
-                job, (huc12, year), callback=_cb, error_callback=_eb
+                job, (_huc12, dt), callback=_cb, error_callback=_eb
             )
         pool.close()
         pool.join()
 
     df = pd.concat(dfs).reset_index().drop(columns=["index"])
-    df.to_feather(f"smstate{year}.feather")
+    df.to_feather(f"smstate{dt:%Y%m%d}.feather")
 
 
 if __name__ == "__main__":
-    main(sys.argv)
+    main()
diff --git a/scripts/dynamic_tillage/workflow.py b/scripts/dynamic_tillage/workflow.py
index fbda4ce5..def01fc1 100644
--- a/scripts/dynamic_tillage/workflow.py
+++ b/scripts/dynamic_tillage/workflow.py
@@ -1,194 +1,110 @@
-"""
-Our Dynamic Tillage Workflow!
+"""Dynamic Tillage Workflow.
+
+On a given date between 15 April and 4 June, we will attempt to plant or till
+fields within a HUC12.  We run this script at 9 PM and needs to be done by
+about 12:30 AM.
+
+The `bootstrap_year.py` script creates database entries with default June
+dates and dates that may not be needed given the actual tillage in play.
 
 """
 
-import sys
-from datetime import date, timedelta
+import os
+import shutil
+import subprocess
+import threading
+from datetime import timedelta
 from multiprocessing import cpu_count
 from multiprocessing.pool import ThreadPool
 
+import click
 import geopandas as gpd
 import numpy as np
 import pandas as pd
-from matplotlib.patches import Rectangle
-
-# The MRMS netcdf file use the IEMRE grid domain
+from pyiem.database import get_dbconnc, get_sqlalchemy_conn
 from pyiem.iemre import SOUTH, WEST, daily_offset
-from pyiem.plot import MapPlot, figure_axes
-from pyiem.util import get_dbconn, get_sqlalchemy_conn, logger, ncopen
+from pyiem.util import archive_fetch, logger, ncopen
 from sqlalchemy import text
 from tqdm import tqdm
 
 LOG = logger()
 CPU_COUNT = min([4, cpu_count() / 2])
-MAKE_PLOTS = False
-DBSET_DATE_THRESHOLD = pd.Timestamp(date.today())
-
-
-def plot_map(i, dt, huc12df, fields):
-    """Make a map diagnostic"""
-    minx, miny, maxx, maxy = huc12df["geom"].to_crs(4326).total_bounds
-    buffer = 0.01
-    huc12 = huc12df.index.values[0]
-
-    mp = MapPlot(
-        title=f"DEP Planting Progress {dt:%Y %b %d} for HUC12: {huc12}",
-        logo="dep",
-        sector="custom",
-        west=minx - buffer,
-        north=maxy + buffer,
-        south=miny - buffer,
-        east=maxx + buffer,
-        caption="Daily Erosion Project",
-        continentalcolor="white",
-    )
-    huc12df.to_crs(mp.panels[0].crs).plot(
-        ax=mp.panels[0].ax,
-        aspect=None,
-        facecolor="None",
-        edgecolor="b",
-        linewidth=2,
-        zorder=100,
-    )
-    fields["color"] = "white"
-    fields.loc[fields["till1"].notna(), "color"] = "tan"
-    fields.loc[fields["plant"].notna(), "color"] = "g"
-    fields.to_crs(mp.panels[0].crs).plot(
-        ax=mp.panels[0].ax,
-        aspect=None,
-        facecolor=fields["color"],
-        edgecolor="k",
-        linewidth=1,
-        zorder=100,
-    )
-    p0 = Rectangle((0, 0), 1, 1, fc="white", ec="k")
-    p1 = Rectangle((0, 0), 1, 1, fc="tan", ec="k")
-    p2 = Rectangle((0, 0), 1, 1, fc="g", ec="k")
-    mp.panels[0].ax.legend(
-        (p0, p1, p2),
-        ("Awaiting", "Tilled", "Planted"),
-        ncol=3,
-        fontsize=11,
-        loc=2,
-    )
-
-    mp.fig.savefig(f"{i:04.0f}.png")
-    mp.close()
-
+PROJDIR = "/opt/dep/prj2wepp"
+EXE = f"{PROJDIR}/prj2wepp"
+
+
+def setup_thread():
+    """Ensure that we have temp folders and sym links constructed."""
+    mydir = f"tmp{threading.get_native_id()}"
+    os.makedirs(f"{mydir}/wepp/runs", exist_ok=True)
+    for dn in ["data", "output", "wepp", "weppwin"]:
+        subprocess.call(["ln", "-s", f"../../wepp/{dn}"], cwd=f"{mydir}/wepp")
+    subprocess.call(["ln", "-s", "../userdb"], cwd=mydir)
+
+
+def prj2wepp(huc12, fpath):
+    """Run prj2wepp."""
+    mydir = f"tmp{threading.get_native_id()}"
+    fn = f"/i/0/prj/{huc12[:8]}/{huc12[8:]}/{huc12}_{fpath:.0f}.prj"
+    testfn = os.path.basename(fn)[:-4]
+    cmd = [EXE, fn, testfn, f"{PROJDIR}/{mydir}/wepp", "no"]
+    with subprocess.Popen(
+        cmd, stderr=subprocess.PIPE, stdout=subprocess.PIPE, cwd=mydir
+    ) as proc:
+        # Need to block the above
+        stdout = proc.stdout.read()
+        stderr = proc.stderr.read()
+        if not os.path.isfile(f"{mydir}/{testfn}.man"):
+            print("bzzzz", testfn)
+            print(" ".join(cmd))
+            print(stdout.decode("ascii"))
+            print(stderr.decode("ascii"))
+            return False
+        # This generates .cli, .man, .run, .slp, .sol
+        # We need the .man , .slp , .sol from this process
+        # slp is generated in flowpath2prj
+        for suffix in ["man", "sol"]:
+            shutil.copyfile(
+                f"{mydir}/{testfn}.{suffix}", fn.replace("prj", suffix)
+            )
+        for suffix in ["cli", "man", "run", "slp", "sol"]:
+            os.unlink(f"{mydir}/{testfn}.{suffix}")
+    return True
 
-def plot_timeseries(year, df, huc12):
-    """Make a diagnostic."""
-    (fig, ax) = figure_axes(
-        logo="dep",
-        title=f"DEP {year} Tillage/Plant Operation Timing for HUC12: {huc12}",
-        subtitle="<=10% Daily Rate, All Field OFEs below 0.9 Plastic Limit",
-    )
-    ax2 = ax.twinx()
-    x = pd.date_range(f"{year}/04/15", f"{year}/06/03")
-    print(df)
-    ax2.bar(
-        x - pd.Timedelta(hours=8),
-        df["acres_planted"],
-        color="#00ff00",
-        alpha=0.5,
-        width=0.3,
-        align="edge",
-        zorder=2,
-    )
-    ax2.bar(
-        x,
-        df["acres_tilled"],
-        color="#ff0000",
-        alpha=0.5,
-        width=0.3,
-        align="edge",
-        zorder=2,
-    )
 
-    ax.plot(
-        x,
-        df["acres_not_planted"],
-        color="k",
-        label="Not Planted",
-        zorder=3,
-        lw=3,
-    )
-    ax.plot(x, df["acres_to_till"], color="r", label="Tillage", zorder=3, lw=3)
-    ax.plot(
-        x,
-        df["acres_to_plant"],
-        c="g",
-        label="Plant",
-        zorder=3,
-        lw=3,
-    )
-    ax.set_ylabel("Acres Available for Operations")
-    ax.legend(loc=1)
-    ax.grid(True)
-    ax.set_ylim(bottom=-10)
-    ax2.set_ylim(bottom=-10)
-    ax2.set_ylabel("Acreas Worked per Day (bars)")
-    ax.set_zorder(ax2.get_zorder() + 1)
-    ax.patch.set_visible(False)
-    fig.savefig(f"{year}_timeseries.png")
-
-
-def dbset(cursor, field_id, dt, col):
-    """Update the database."""
-    if dt >= DBSET_DATE_THRESHOLD:
-        return
-    cursor.execute(
-        f"UPDATE field_operations SET {col} = %s WHERE field_id = %s and "
-        "year = %s",
-        (dt, field_id, dt.year),
+def edit_rotfile(year, huc12, row):
+    """Edit up our .rot files."""
+    yearindex = str(year - 2007 + 1)
+    rotfn = (
+        f"/i/0/rot/{huc12[:8]}/{huc12[8:]}/"
+        f"{huc12}_{row['fpath']:.0f}_{row['ofe']:.0f}.rot"
     )
-    if cursor.rowcount == 0:
-        cursor.execute(
-            f"INSERT into field_operations(field_id, year, {col}) "
-            "VALUES(%s, %s, %s)",
-            (field_id, dt.year, dt),
-        )
+    with open(rotfn, encoding="ascii") as fh:
+        lines = fh.readlines()
+    dates = [row["till1"], row["till2"], row["till3"], row["plant"]]
+    for i, line in enumerate(lines):
+        pos = line.find(f" {yearindex} ")  # false positives
+        if pos == -1:
+            continue
+        tokens = line.split()
+        if (
+            tokens[2] == yearindex
+            and int(tokens[0]) < 7
+            and tokens[4].startswith(("Tillage", "Plant"))
+        ):
+            thisdt = row["plant"]
+            if tokens[4] == "Tillage":
+                thisdt = dates.pop(0)
+            if thisdt is None:
+                thisdt = row["plant"]
+            lines[i] = (
+                f"{thisdt:%-m %-d} {yearindex} " f"{' '.join(tokens[3:])}\n"
+            )
+    with open(rotfn, "w", encoding="ascii") as fh:
+        fh.write("".join(lines))
 
 
-def edit_rotfile(year, huc12, ofedf):
-    """Edit up our .rot files."""
-    yearindex = str(year - 2007 + 1)
-    for _, row in ofedf.iterrows():
-        dates = []
-        for col in ["till1", "till2", "plant", "plant"]:  # Hack dupe plant
-            if not pd.isna(row[col]):
-                dates.append(row[col])
-        rotfn = (
-            f"/i/0/rot/{huc12[:8]}/{huc12[8:]}/"
-            f"{huc12}_{row['fpath']}_{row['ofe']}.rot"
-        )
-        with open(rotfn, encoding="ascii") as fh:
-            lines = fh.readlines()
-        for i, line in enumerate(lines):
-            pos = line.find(f" {yearindex} ")  # false positives
-            if pos == -1:
-                continue
-            tokens = line.split()
-            if (
-                tokens[2] == yearindex
-                and int(tokens[0]) < 6
-                and (
-                    tokens[4].startswith("Tillage")
-                    or tokens[4].startswith("Plant")
-                )
-            ):
-                lines[i] = (
-                    f"{dates.pop(0):%m %d} {yearindex} "
-                    f"{' '.join(tokens[3:])}\n"
-                )
-                if not dates:
-                    break
-        with open(rotfn, "w", encoding="ascii") as fh:
-            fh.write("".join(lines))
-
-
-def do_huc12(year, huc12, smdf):
+def do_huc12(dt, huc12, smdf):
     """Go Main Go."""
     with get_sqlalchemy_conn("idep") as conn:
         # build up the cross reference of everyhing we need to know
@@ -196,163 +112,130 @@ def do_huc12(year, huc12, smdf):
             text(
                 """
             with myofes as (
-                select o.ofe, p.fpath, o.fbndid,
-                g.plastic_limit
+                select o.ofe, p.fpath, o.fbndid, g.plastic_limit
                 from flowpaths p, flowpath_ofes o, gssurgo g
                 WHERE o.flowpath = p.fid and p.huc_12 = :huc12
                 and p.scenario = 0 and o.gssurgo_id = g.id),
             agg as (
                 SELECT ofe, fpath, f.fbndid, plastic_limit, f.geom, f.landuse,
                 f.management, f.acres, f.field_id
-                from fields f JOIN myofes m on (f.fbndid = m.fbndid)
-                WHERE f.huc12 = :huc12 ORDER by fpath, ofe)
-            select a.*, o.till1, o.till2, o.plant
-            from agg a LEFT JOIN field_operations o on
-            (a.field_id = o.field_id and o.year = :year)
+                from fields f LEFT JOIN myofes m on (f.fbndid = m.fbndid)
+                WHERE f.huc12 = :huc12 and f.isag ORDER by fpath, ofe)
+            select a.*, o.till1, o.till2, o.till3, o.plant, o.year,
+            plant >= :dt as plant_needed,
+            coalesce(
+                greatest(till1, till2, till3) >= :dt, false) as till_needed,
+            false as operation_done
+            from agg a, field_operations o
+            where a.field_id = o.field_id and o.year = :year
+            and o.plant is not null
             """
             ),
             conn,
-            params={"huc12": huc12, "year": year},
+            params={"huc12": huc12, "year": dt.year, "dt": dt},
             index_col=None,
             geom_col="geom",
         )
+    if df.empty:
+        LOG.info("No fields found for %s", huc12)
+        return
 
-    # Compute things for year
-    char_at = (year - 2007) + 1
     fields = df.groupby("fbndid").first().copy()
+
+    # Soil Moisture check.  We are not modeling every field, so we are a
+    # bit generous here, are 50% of **model** fields below plastic_limit?
+    sw1 = fields[fields["plastic_limit"].notnull()].join(
+        smdf[["fbndid", "sw1"]].groupby("fbndid").min(), on="fbndid"
+    )[["sw1", "plastic_limit"]]
+    hits = (sw1["sw1"] < sw1["plastic_limit"]).sum()
+    if hits < len(sw1.index) / 2:
+        LOG.info(
+            "HUC12: %s has %s/%s fields below plastic limit",
+            huc12,
+            hits,
+            len(sw1.index),
+        )
+        return
+
+    # Compute things for year
+    char_at = (dt.year - 2007) + 1
     fields["geom"].crs = 5070
     fields["crop"] = fields["landuse"].str.slice(char_at - 1, char_at)
     fields["tillage"] = fields["management"].str.slice(char_at - 1, char_at)
-    # augh, Pasture
-    fields = fields[fields["crop"] != "P"]
-
-    fields["status"] = 0
-    # Update status based on what has been done already
-    fields.loc[fields["plant"].notna(), "status"] += 1
-    fields.loc[fields["till1"].notna(), "status"] += 1
-    fields.loc[fields["till2"].notna(), "status"] += 1
-    fields["operations"] = 1
-    fields.loc[fields["tillage"] == "2", "operations"] = 2
-    fields.loc[fields["tillage"].isin(["3", "4", "5", "6"]), "operations"] = 3
+
     total_acres = fields["acres"].sum()
-    operation_acres = (fields["operations"] * fields["acres"]).sum()
-    # 10% acres can get planted and 10% of acres can get tilled
-    LOG.info(
-        "Theoretical min days %s", operation_acres / (total_acres / 10.0) / 2.0
-    )
     limit = total_acres / 10.0
 
-    acres_not_planted = []
-    acres_to_plant = []
-    acres_planted = []
-    acres_to_till = []
-    acres_tilled = []
-
-    pgconn = get_dbconn("idep")
-    cursor = pgconn.cursor()
-    if year == DBSET_DATE_THRESHOLD.year:
-        startdate = max([pd.Timestamp(f"{year}/04/15"), DBSET_DATE_THRESHOLD])
-    else:
-        startdate = pd.Timestamp(f"{year}/04/15")
-    # NOTE once we get to 1 June, everything goes into the ground, so we need
-    # three days to clear the backlog...
-    for i, dt in enumerate(pd.date_range(startdate, f"{year}/06/03")):
-        acres_not_planted.append(fields[fields["plant"].isna()]["acres"].sum())
-        # Can we plant upto 10% of the acreage today?
-        running = 0
-        pop = fields[
-            fields["plant"].isna()
-            & ((fields["operations"] - fields["status"]) == 1)
-        ]
-        acres_to_plant.append(pop["acres"].sum())
-        LOG.info("%s acres to plant: %s", dt, pop["acres"].sum())
-        yesterday = dt - timedelta(days=1)
-        # random iter
-        for fbndid, row in pop.sample(frac=1).iterrows():
-            # Go get the soil moisture state for yesterday as it would be
-            # what is valid at 12 AM today.
-            sm = smdf[(smdf["fbndid"] == fbndid) & (smdf["date"] == yesterday)]
-            # Are we all below?
-            if dt.month < 6 and (sm["sw1"] > row["plastic_limit"]).any():
-                continue
-            fields.at[fbndid, "plant"] = dt
-            dbset(cursor, row["field_id"], dt, "plant")
-            running += row["acres"]
-            if dt.month < 6 and running > limit:
+    acres_tilled = 0
+    # Work on tillage first, so to avoid planting on tilled fields
+    for fbndid, row in fields[fields["till_needed"]].iterrows():
+        acres_tilled += row["acres"]
+        fields.at[fbndid, "operation_done"] = True
+        for i in [1, 2, 3]:
+            if row[f"till{i}"] >= dt:
+                fields.at[fbndid, f"till{i}"] = dt
                 break
-        acres_planted.append(running)
-        LOG.info(
-            "%s acres to plant: %s, did %s", dt, pop["acres"].sum(), running
-        )
-        # Can we do 10% of tillage operations
-        running = 0
-        pop = fields[
-            fields["plant"].isna()
-            & ((fields["operations"] - fields["status"]) > 1)
-        ]
-        acres_to_till.append(pop["acres"].sum())
-        # random iter
-        for fbndid, row in pop.sample(frac=1).iterrows():
-            # Go get the soil moisture state
-            sm = smdf[(smdf["fbndid"] == fbndid) & (smdf["date"] == yesterday)]
-            # Are we all below?
-            if dt.month < 6 and (sm["sw1"] > row["plastic_limit"]).any():
-                continue
-            fields.at[fbndid, f"till{row['status'] + 1}"] = dt
-            fields.at[fbndid, "status"] += 1
-            dbset(
-                cursor,
-                row["field_id"],
-                dt,
-                f"till{fields.at[fbndid, 'status']}",
-            )
-            running += row["acres"]
-            if dt.month < 6 and running > limit:
-                break
-        acres_tilled.append(running)
-        LOG.info(
-            "%s acres to till: %s, did %s", dt, pop["acres"].sum(), running
-        )
-    # Merge field info back into main dataframe
-    df = df[["ofe", "fpath", "fbndid"]].merge(
-        fields[["till1", "till2", "plant", "crop", "tillage"]],
-        left_on="fbndid",
-        right_index=True,
+        if acres_tilled > limit:
+            break
+
+    # Now we need to plant
+    acres_planted = 0
+    for fbndid, row in fields[fields["plant_needed"]].iterrows():
+        if row["operation_done"]:
+            continue
+        acres_planted += row["acres"]
+        fields.at[fbndid, "plant"] = dt
+        fields.at[fbndid, "operation_done"] = True
+        if acres_planted > limit:
+            break
+
+    LOG.info(
+        "plant: %.1f[%.1f%%] tilled: %.1f[%.1f%%] total: %.1f",
+        acres_planted,
+        acres_planted / total_acres * 100.0,
+        acres_tilled,
+        acres_tilled / total_acres * 100.0,
+        total_acres,
     )
-    edit_rotfile(year, huc12, df[df["crop"].isin(["C", "B"])])
-    cursor.close()
-    pgconn.commit()
-    if MAKE_PLOTS:
-        plot_timeseries(
-            year,
-            pd.DataFrame(
-                {
-                    "acres_to_plant": acres_to_plant,
-                    "acres_to_till": acres_to_till,
-                    "acres_not_planted": acres_not_planted,
-                    "acres_planted": acres_planted,
-                    "acres_tilled": acres_tilled,
-                }
-            ),
-            huc12,
-        )
+
+    # Update the database
+    df2 = fields[fields["operation_done"]].reset_index()
+    if df2.empty:
+        return
+    conn, cursor = get_dbconnc("idep")
+    cursor.executemany(
+        "update field_operations set plant = %s, till1 = %s, "
+        "till2 = %s, till3 = %s where field_id = %s and year = %s",
+        df2[
+            ["plant", "till1", "till2", "till3", "field_id", "year"]
+        ].itertuples(index=False),
+    )
+    # FIXME cursor.close()
+    # FIXME conn.commit()
+
+    # Update the .rot files, when necessary
+    df2 = df[(df["field_id"].isin(df2["field_id"])) & (df["ofe"].notna())]
+    for _, row in df2.iterrows():
+        edit_rotfile(dt.year, huc12, row)
+    for fpath in df2["fpath"].unique():
+        prj2wepp(huc12, fpath)
 
 
-def estimate_soiltemp(huc12df):
-    """GFS soil temperature."""
-    with ncopen("/mesonet/data/iemre/gfs_current.nc") as nc:
+def estimate_soiltemp(huc12df, dt):
+    """Write mean GFS soil temperature C into huc12df."""
+    ppath = f"{dt:%Y/%m/%d}/model/gfs/gfs_{dt:%Y%m%d}00_iemre.nc"
+    with archive_fetch(ppath) as ncfn, ncopen(ncfn) as nc:
         tsoil = np.mean(nc.variables["tsoil"][1:7], axis=0) - 273.15
         y = np.digitize(huc12df["lat"].values, nc.variables["lat"][:])
         x = np.digitize(huc12df["lon"].values, nc.variables["lon"][:])
-    for i, idx in enumerate(huc12df.index.values):
-        huc12df.at[idx, "tsoil"] = tsoil[y[i], x[i]]
+        for i, idx in enumerate(huc12df.index.values):
+            huc12df.at[idx, "tsoil"] = tsoil[y[i], x[i]]
 
 
-def estimate_rainfall(huc12df):
+def estimate_rainfall(huc12df, dt):
     """Figure out our rainfall estimate."""
-    today = date.today()
-    tidx = daily_offset(today)
-    with ncopen(f"/mesonet/data/mrms/{today:%Y}_mrms_daily.nc") as nc:
+    tidx = daily_offset(dt)
+    with ncopen(f"/mesonet/data/mrms/{dt:%Y}_mrms_daily.nc") as nc:
         p01d = nc.variables["p01d"][tidx].filled(0)
     for idx, row in huc12df.iterrows():
         y = int((row["lat"] - SOUTH) * 100.0)
@@ -360,42 +243,51 @@ def estimate_rainfall(huc12df):
         huc12df.at[idx, "precip_mm"] = p01d[y, x]
 
 
-def main(argv):
+@click.command()
+@click.option("--scenario", type=int, default=0)
+@click.option("--date", "dt", type=click.DateTime(), required=True)
+@click.option("--huc12", type=str, required=False)
+def main(scenario, dt, huc12):
     """Go Main Go."""
-    scenario = int(argv[1])
-    year = int(argv[2])
+    # Deal with dates not datetime
+    dt = dt.date()
+    LOG.info("Processing %s for scenario %s", dt, scenario)
+    huc12_filter = "" if huc12 is None else " and huc_12 = :huc12"
     with get_sqlalchemy_conn("idep") as conn:
         huc12df = gpd.read_postgis(
-            """
+            text(f"""
             SELECT geom, huc_12,
+            false as limited_by_precip,
+            false as limited_by_soiltemp,
+            false as limited_by_soilmoisture,
+            false as limited,
             ST_x(st_transform(st_centroid(geom), 4326)) as lon,
             ST_y(st_transform(st_centroid(geom), 4326)) as lat
-            from huc12 where scenario = %s
-            """,
+            from huc12 where scenario = :scenario {huc12_filter}
+            """),
             conn,
-            params=(scenario,),
+            params={"scenario": scenario, "huc12": huc12},
             geom_col="geom",
             index_col="huc_12",
         )
-    # Allows for isolated testing.
-    if len(argv) == 4:
-        huc12df = huc12df.loc[[argv[3]]]
     # Estimate today's rainfall so to delay triggers
-    estimate_rainfall(huc12df)
+    estimate_rainfall(huc12df, dt)
     # Any HUC12 over 10mm gets skipped
-    huc12df = huc12df[huc12df["precip_mm"] < 10]
-    if year == date.today().year:
-        # Estimate soil temperature via GFS forecast
-        estimate_soiltemp(huc12df)
-        # Any HUC12 under 6C gets skipped
-        huc12df = huc12df[huc12df["tsoil"] > 5.9]
+    huc12df.loc[huc12df["precip_mm"] > 9.9]["limited_by_precip"] = True
+    # Estimate soil temperature via GFS forecast
+    estimate_soiltemp(huc12df, dt)
+    # Any HUC12 under 6C gets skipped
+    huc12df.loc[huc12df["tsoil"] > 5.9]["limited_by_soiltemp"] = True
+    # The soil moisture state is the day before
+    yesterday = dt - timedelta(days=1)
+    smdf = pd.read_feather(f"smstate{yesterday:%Y%m%d}.feather")
+
+    # Need to run from project dir so local userdb is found
+    os.chdir(PROJDIR)
 
-    LOG.warning("%s threads for %s huc12s", CPU_COUNT, len(huc12df.index))
     progress = tqdm(total=len(huc12df.index))
-
-    smdf = pd.read_feather(f"smstate{year}.feather")
-
-    with ThreadPool(CPU_COUNT) as pool:
+    LOG.warning("%s threads for %s huc12s", CPU_COUNT, len(huc12df.index))
+    with ThreadPool(CPU_COUNT, initializer=setup_thread) as pool:
 
         def _error(exp):
             """Uh oh."""
@@ -406,7 +298,7 @@ def _error(exp):
         def _job(huc12):
             """Do the job."""
             progress.set_description(huc12)
-            do_huc12(year, huc12, smdf[smdf["huc12"] == huc12])
+            do_huc12(dt, huc12, smdf[smdf["huc12"] == huc12])
             progress.update()
 
         for huc12 in huc12df.index.values:
@@ -416,4 +308,4 @@ def _job(huc12):
 
 
 if __name__ == "__main__":
-    main(sys.argv)
+    main()
diff --git a/scripts/gridorder/flowpathlength_totals.py b/scripts/gridorder/flowpathlength_totals.py
index 4dd2848c..a5d0c960 100644
--- a/scripts/gridorder/flowpathlength_totals.py
+++ b/scripts/gridorder/flowpathlength_totals.py
@@ -102,7 +102,7 @@ def load_lengths():
 
     # Begin the processing work now!
     pool = multiprocessing.Pool()
-    for df, huc12, slopes in tqdm(
+    for df, huc12, _slopes in tqdm(
         pool.imap_unordered(do_huc12, huc12s),
         total=len(huc12s),
         disable=(not sys.stdout.isatty()),
diff --git a/scripts/gridorder/single_huc12_map.py b/scripts/gridorder/single_huc12_map.py
index dfaddf8f..26add536 100644
--- a/scripts/gridorder/single_huc12_map.py
+++ b/scripts/gridorder/single_huc12_map.py
@@ -67,7 +67,7 @@ def main(argv):
         geom_col="geo",
         index_col="fpath",
     )
-    for fpath, row in df.iterrows():
+    for _fpath, row in df.iterrows():
         points = row["geo"].xy
         ax.plot(
             points[0],
diff --git a/scripts/import/flowpath2prj.py b/scripts/import/flowpath2prj.py
index ebfdbb02..7cc021dd 100644
--- a/scripts/import/flowpath2prj.py
+++ b/scripts/import/flowpath2prj.py
@@ -16,7 +16,6 @@
   L - Double Crop  (started previous year) (winter wheat + soybea)
   W - Wheat
   S - Durum/Spring Wheat
-  N - Small Grain
   E - Sugarbeets
   J - Rice
   N - Cotton
@@ -39,8 +38,10 @@
 from math import atan2, degrees, pi
 
 import pandas as pd
+from pydep.tillage import make_tillage
 from pydep.util import load_scenarios
-from pyiem.util import get_dbconn, get_sqlalchemy_conn, logger
+from pyiem.database import get_dbconn, get_sqlalchemy_conn
+from pyiem.util import logger
 from sqlalchemy import text
 from tqdm import tqdm
 
@@ -57,31 +58,6 @@
     "P": "IniCropDef.gra_3425",
     "R": "IniCropDef.Aft_12889",
 }
-WHEAT_PLANT = {
-    "KS_SOUTH": datetime.date(2000, 5, 25),
-    "KS_CENTRAL": datetime.date(2000, 5, 25),
-    "KS_NORTH": datetime.date(2000, 5, 25),
-    "IA_SOUTH": datetime.date(2000, 5, 12),
-    "IA_CENTRAL": datetime.date(2000, 5, 17),
-    "IA_NORTH": datetime.date(2000, 5, 23),
-}
-
-SOYBEAN_PLANT = {
-    "KS_SOUTH": datetime.date(2000, 5, 25),
-    "KS_CENTRAL": datetime.date(2000, 5, 25),
-    "KS_NORTH": datetime.date(2000, 5, 25),
-    "IA_SOUTH": datetime.date(2000, 5, 12),
-    "IA_CENTRAL": datetime.date(2000, 5, 17),
-    "IA_NORTH": datetime.date(2000, 5, 23),
-}
-CORN_PLANT = {
-    "KS_SOUTH": datetime.date(2000, 4, 20),
-    "KS_CENTRAL": datetime.date(2000, 4, 25),
-    "KS_NORTH": datetime.date(2000, 4, 30),
-    "IA_SOUTH": datetime.date(2000, 4, 30),
-    "IA_CENTRAL": datetime.date(2000, 5, 5),
-    "IA_NORTH": datetime.date(2000, 5, 10),
-}
 # Management code maps back to a forest crop
 FOREST = {
     "J": "IniCropDef.DEP_forest_95",
@@ -96,124 +72,6 @@
     "A": "IniCropDef.DEP_forest_5",
     "0": "IniCropDef.DEP_forest_45",  # TODO
 }
-CORN = {
-    "KS_SOUTH": "CropDef.Cor_0964",
-    "KS_CENTRAL": "CropDef.Cor_0964",
-    "KS_NORTH": "CropDef.Cor_0964",
-    "IA_SOUTH": "CropDef.Cor_0965",
-    "IA_CENTRAL": "CropDef.Cor_0966",
-    "IA_NORTH": "CropDef.Cor_0967",
-}
-SOYBEAN = {
-    "KS_SOUTH": "CropDef.Soy_2191",
-    "KS_CENTRAL": "CropDef.Soy_2191",
-    "KS_NORTH": "CropDef.Soy_2191",
-    "IA_SOUTH": "CropDef.Soy_2192",
-    "IA_CENTRAL": "CropDef.Soy_2193",
-    "IA_NORTH": "CropDef.Soy_2194",
-}
-WHEAT = {
-    "KS_SOUTH": "CropDef.spwheat1",
-    "KS_CENTRAL": "CropDef.spwheat1",
-    "KS_NORTH": "CropDef.spwheat1",
-    "IA_SOUTH": "CropDef.spwheat1",
-    "IA_CENTRAL": "CropDef.spwheat1",
-    "IA_NORTH": "CropDef.spwheat1",
-}
-
-
-def read_file(scenario, zone, prevcode, code, nextcode, cfactor, year):
-    """Read a block file and do replacements
-
-    Args:
-      scenario (int): The DEP scenario
-      zone (str): The DEP cropping zone
-      prevcode (str): the previous crop code
-      code (str): the crop code
-      nextcode (str): The next year's crop code.
-      cfactor (int): the c-factor for tillage
-      year (int): the year of this crop
-
-    Returns:
-      str with the raw data used for the .rot file
-    """
-    blockfn = f"blocks/{code}{cfactor}.txt"
-    if not os.path.isfile(blockfn):
-        return ""
-    with open(blockfn, "r", encoding="utf8") as fh:
-        data = fh.read()
-    # Special consideration for planting alfalfa
-    if code == "P" and prevcode != "P":
-        # Best we can do now is plant it on Apr 15, sigh
-        data = (
-            f"4  15 {year} 1 Plant-Perennial CropDef.ALFALFA  {{0.000000}}\n"
-            f"{data}"
-        )
-    # Remove fall chisel after corn when going into soybeans for 2
-    if cfactor == 2 and nextcode == "B":
-        pos = data.find("11  1")
-        if pos > 0:
-            data = data[:pos]
-
-    # The fall tillage operation is governed by the next year crop
-    if cfactor == 5 and nextcode == "C":
-        # Replace 11  1 operation with plow
-        pos = data.find("11  1")
-        if pos > 0:
-            data = (
-                f"{data[:pos]}"
-                "11  1  %(yr)s  1 Tillage      OpCropDef.MOPL  {0.203200, 1}\n"
-            )
-
-    # Anhydrous ammonia application when we are going into Corn
-    if nextcode == "C":
-        # HACK: look for a present 1 Nov operation and insert this before it
-        pos = data.find("11  1")
-        extra = ""
-        if pos > 0:
-            extra = data[pos:].replace("11  1", "11  8")
-            data = data[:pos]
-        data = (
-            f"{data}"
-            f"11  1  {year}  1 Tillage   OpCropDef.ANHYDROS  {{0.101600, 2}}\n"
-            f"{extra}"
-        )
-
-    pdate = ""
-    pdatem5 = ""
-    pdatem10 = ""
-    pdatem15 = ""
-    plant = ""
-    # We currently only have zone specific files for Corn and Soybean
-    if code == "C":
-        date = CORN_PLANT.get(scenario, CORN_PLANT[zone])
-        pdate = date.strftime("%m    %d")
-        pdatem5 = (date - datetime.timedelta(days=5)).strftime("%m    %d")
-        pdatem10 = (date - datetime.timedelta(days=10)).strftime("%m    %d")
-        pdatem15 = (date - datetime.timedelta(days=15)).strftime("%m    %d")
-        plant = CORN[zone]
-    elif code in ["B", "L"]:  # TODO support double crop
-        date = SOYBEAN_PLANT.get(scenario, SOYBEAN_PLANT[zone])
-        pdate = date.strftime("%m    %d")
-        pdatem5 = (date - datetime.timedelta(days=5)).strftime("%m    %d")
-        pdatem10 = (date - datetime.timedelta(days=10)).strftime("%m    %d")
-        pdatem15 = (date - datetime.timedelta(days=15)).strftime("%m    %d")
-        plant = SOYBEAN[zone]
-    elif code == "W":
-        date = SOYBEAN_PLANT.get(scenario, SOYBEAN_PLANT[zone])  # TODO
-        plant = WHEAT_PLANT[zone]
-        pdate = date.strftime("%m    %d")
-        pdatem5 = (date - datetime.timedelta(days=5)).strftime("%m    %d")
-        pdatem10 = (date - datetime.timedelta(days=10)).strftime("%m    %d")
-        pdatem15 = (date - datetime.timedelta(days=15)).strftime("%m    %d")
-    return data % {
-        "yr": year,
-        "pdate": pdate,
-        "pdatem5": pdatem5,
-        "pdatem10": pdatem10,
-        "pdatem15": pdatem15,
-        "plant": plant,
-    }
 
 
 def do_rotation(scenario, zone, rotfn, landuse, management):
@@ -243,7 +101,7 @@ def do_rotation(scenario, zone, rotfn, landuse, management):
     else:
         data["initcond"] = INITIAL_COND.get(landuse[0], INITIAL_COND_DEFAULT)
         prevcode = "C" if landuse[0] not in INITIAL_COND else landuse[0]
-        func = partial(read_file, scenario, zone)
+        func = partial(make_tillage, scenario, zone)
         # 2007
         data["yearly"] += func(
             prevcode, landuse[0], landuse[1], int(management[0]), 1
@@ -531,7 +389,7 @@ def generate_combos():
     """Create rotation files for Eduardo inspection."""
     for cfactor in range(1, 7):
         for rot in "CC CB BC BB".split():
-            res = read_file(0, "IA_NORTH", "M", rot[0], rot[1], cfactor, 1)
+            res = make_tillage(0, "IA_NORTH", "M", rot[0], rot[1], cfactor, 1)
             with open(f"ex_{rot}_{cfactor}.rot", "w", encoding="utf8") as fh:
                 fh.write(res)
 
diff --git a/scripts/plots/yearly_summary.py b/scripts/plots/yearly_summary.py
index fef29dce..017b2392 100644
--- a/scripts/plots/yearly_summary.py
+++ b/scripts/plots/yearly_summary.py
@@ -7,8 +7,8 @@
 import numpy as np
 from geopandas import read_postgis
 from matplotlib.patches import Polygon
+from pyiem.database import get_dbconn
 from pyiem.plot import MapPlot
-from pyiem.util import get_dbconn
 
 V2NAME = {
     "avg_loss": "Detachment",
@@ -104,7 +104,7 @@ def main():
     norm = mpcolors.BoundaryNorm(bins, cmap.N)
 
     # m.ax.add_geometries(df['geo'], ccrs.PlateCarree())
-    for i, row in df.iterrows():
+    for _, row in df.iterrows():
         c = cmap(norm([row["data"]]))[0]
         arr = np.asarray(row["geo"].exterior)
         points = m.ax.projection.transform_points(
diff --git a/scripts/switchgrass/mlra_percent_slopes_conv.py b/scripts/switchgrass/mlra_percent_slopes_conv.py
index bd78dc7f..2f28c758 100644
--- a/scripts/switchgrass/mlra_percent_slopes_conv.py
+++ b/scripts/switchgrass/mlra_percent_slopes_conv.py
@@ -1,7 +1,7 @@
 """Print out the percentage of slopes converted."""
 
 from pandas.io.sql import read_sql
-from pyiem.util import get_dbconn
+from pyiem.database import get_dbconn
 
 LABELS = {
     36: "Slopes > 3% to Switchgrass",
@@ -24,7 +24,7 @@ def main(argv):
     )
 
     print("%s," % (mlraxref.at[mlra_id, "mlra_name"],), end="")
-    for col, scenario in enumerate(range(36, 39)):
+    for _, scenario in enumerate(range(36, 39)):
         df = read_sql(
             """
         with myhucs as (
diff --git a/scripts/tillage/flowpath2prj.py b/scripts/tillage/flowpath2prj.py
index 59f2ad48..56b5ef1b 100644
--- a/scripts/tillage/flowpath2prj.py
+++ b/scripts/tillage/flowpath2prj.py
@@ -8,7 +8,7 @@
 
 from psycopg.rows import dict_row
 from pydep.util import get_cli_fname
-from pyiem.util import get_dbconn
+from pyiem.database import get_dbconn
 from tqdm import tqdm
 
 SCENARIO = int(sys.argv[1])
@@ -52,7 +52,7 @@ def simplify(rows):
     newrows = []
     lrow = rows[0]
     newrows.append(lrow)
-    for i, row in enumerate(rows):
+    for _, row in enumerate(rows):
         # If they are the 'same', continue
         if row["landuse"] == "UUUUUUUUUU":
             continue
diff --git a/scripts/ucs/isolatehuc12.py b/scripts/ucs/isolatehuc12.py
index f87e8b22..04f2fb1b 100644
--- a/scripts/ucs/isolatehuc12.py
+++ b/scripts/ucs/isolatehuc12.py
@@ -1,6 +1,6 @@
 import matplotlib.pyplot as plt
 from pandas.io.sql import read_sql
-from pyiem.util import get_dbconn
+from pyiem.database import get_dbconn
 
 pgconn = get_dbconn("idep")
 df = read_sql(
@@ -21,7 +21,7 @@
     x = []
     y = []
     accum = 0
-    for i, row in df2.iterrows():
+    for _, row in df2.iterrows():
         x.append(row["valid"])
         accum += row["avg_loss"] * 4.463
         y.append(accum)
diff --git a/scripts/ucs/top25counties.py b/scripts/ucs/top25counties.py
index 8ec9fc4d..5137765a 100644
--- a/scripts/ucs/top25counties.py
+++ b/scripts/ucs/top25counties.py
@@ -3,7 +3,7 @@
 import sys
 
 from pandas.io.sql import read_sql
-from pyiem.util import get_dbconn
+from pyiem.database import get_dbconn
 
 pgconn = get_dbconn("idep")
 
@@ -29,7 +29,7 @@
 cursor = pgconn.cursor()
 
 DONE = []
-for i, row in df.iterrows():
+for _, row in df.iterrows():
     cursor.execute(
         """SELECT name from ugcs where end_ts is null and
     state = 'IA' and substr(ugc, 3, 1) = 'C' and
diff --git a/src/pydep/io/wepp.py b/src/pydep/io/wepp.py
index aab7a9fa..4241dd06 100644
--- a/src/pydep/io/wepp.py
+++ b/src/pydep/io/wepp.py
@@ -100,7 +100,7 @@ def read_cli(
     filename,
     compute_rfactor=False,
     return_rfactor_metric=True,
-    compute_intensity_over=[],
+    compute_intensity_over=None,
 ):
     """Read WEPP CLI File, Return DataFrame
 
@@ -162,7 +162,7 @@ def read_cli(
                 )
             ),
         }
-        if compute_intensity_over:
+        if compute_intensity_over is not None:
             row.update(intensity(times, points, compute_intensity_over))
         rows.append(row)
 
diff --git a/src/pydep/tillage.py b/src/pydep/tillage.py
new file mode 100644
index 00000000..23723539
--- /dev/null
+++ b/src/pydep/tillage.py
@@ -0,0 +1,149 @@
+"""DEP Tillage Logic and Life Choices."""
+
+import os
+from datetime import date, timedelta
+
+WHEAT_PLANT = {
+    "KS_SOUTH": date(2000, 5, 25),
+    "KS_CENTRAL": date(2000, 5, 25),
+    "KS_NORTH": date(2000, 5, 25),
+    "IA_SOUTH": date(2000, 5, 12),
+    "IA_CENTRAL": date(2000, 5, 17),
+    "IA_NORTH": date(2000, 5, 23),
+}
+
+SOYBEAN_PLANT = {
+    "KS_SOUTH": date(2000, 5, 25),
+    "KS_CENTRAL": date(2000, 5, 25),
+    "KS_NORTH": date(2000, 5, 25),
+    "IA_SOUTH": date(2000, 5, 12),
+    "IA_CENTRAL": date(2000, 5, 17),
+    "IA_NORTH": date(2000, 5, 23),
+}
+CORN_PLANT = {
+    "KS_SOUTH": date(2000, 4, 20),
+    "KS_CENTRAL": date(2000, 4, 25),
+    "KS_NORTH": date(2000, 4, 30),
+    "IA_SOUTH": date(2000, 4, 30),
+    "IA_CENTRAL": date(2000, 5, 5),
+    "IA_NORTH": date(2000, 5, 10),
+}
+CORN = {
+    "KS_SOUTH": "CropDef.Cor_0964",
+    "KS_CENTRAL": "CropDef.Cor_0964",
+    "KS_NORTH": "CropDef.Cor_0964",
+    "IA_SOUTH": "CropDef.Cor_0965",
+    "IA_CENTRAL": "CropDef.Cor_0966",
+    "IA_NORTH": "CropDef.Cor_0967",
+}
+SOYBEAN = {
+    "KS_SOUTH": "CropDef.Soy_2191",
+    "KS_CENTRAL": "CropDef.Soy_2191",
+    "KS_NORTH": "CropDef.Soy_2191",
+    "IA_SOUTH": "CropDef.Soy_2192",
+    "IA_CENTRAL": "CropDef.Soy_2193",
+    "IA_NORTH": "CropDef.Soy_2194",
+}
+WHEAT = {
+    "KS_SOUTH": "CropDef.spwheat1",
+    "KS_CENTRAL": "CropDef.spwheat1",
+    "KS_NORTH": "CropDef.spwheat1",
+    "IA_SOUTH": "CropDef.spwheat1",
+    "IA_CENTRAL": "CropDef.spwheat1",
+    "IA_NORTH": "CropDef.spwheat1",
+}
+
+
+def make_tillage(scenario, zone, prevcode, code, nextcode, cfactor, year):
+    """Read a block file and do replacements
+
+    Args:
+      scenario (int): The DEP scenario
+      zone (str): The DEP cropping zone
+      prevcode (str): the previous crop code
+      code (str): the crop code
+      nextcode (str): The next year's crop code.
+      cfactor (int): the c-factor for tillage
+      year (int): the year of this crop
+
+    Returns:
+      str with the raw data used for the .rot file
+    """
+    # FIXME
+    blockfn = f"/opt/dep/scripts/import/blocks/{code}{cfactor}.txt"
+    if not os.path.isfile(blockfn):
+        return ""
+    with open(blockfn, "r", encoding="utf8") as fh:
+        data = fh.read()
+    # Special consideration for planting alfalfa
+    if code == "P" and prevcode != "P":
+        # Best we can do now is plant it on Apr 15, sigh
+        data = (
+            f"4  15 {year} 1 Plant-Perennial CropDef.ALFALFA  {{0.000000}}\n"
+            f"{data}"
+        )
+    # Remove fall chisel after corn when going into soybeans for 2
+    if cfactor == 2 and nextcode == "B":
+        pos = data.find("11  1")
+        if pos > 0:
+            data = data[:pos]
+
+    # The fall tillage operation is governed by the next year crop
+    if cfactor == 5 and nextcode == "C":
+        # Replace 11  1 operation with plow
+        pos = data.find("11  1")
+        if pos > 0:
+            data = (
+                f"{data[:pos]}"
+                "11  1  %(yr)s  1 Tillage      OpCropDef.MOPL  {0.203200, 1}\n"
+            )
+
+    # Anhydrous ammonia application when we are going into Corn
+    if nextcode == "C":
+        # HACK: look for a present 1 Nov operation and insert this before it
+        pos = data.find("11  1")
+        extra = ""
+        if pos > 0:
+            extra = data[pos:].replace("11  1", "11  8")
+            data = data[:pos]
+        data = (
+            f"{data}"
+            f"11  1  {year}  1 Tillage   OpCropDef.ANHYDROS  {{0.101600, 2}}\n"
+            f"{extra}"
+        )
+
+    pdate = ""
+    pdatem5 = ""
+    pdatem10 = ""
+    pdatem15 = ""
+    plant = ""
+    # We currently only have zone specific files for Corn and Soybean
+    if code == "C":
+        date = CORN_PLANT.get(scenario, CORN_PLANT[zone])
+        pdate = date.strftime("%m    %d")
+        pdatem5 = (date - timedelta(days=5)).strftime("%m    %d")
+        pdatem10 = (date - timedelta(days=10)).strftime("%m    %d")
+        pdatem15 = (date - timedelta(days=15)).strftime("%m    %d")
+        plant = CORN[zone]
+    elif code in ["B", "L"]:  # TODO support double crop
+        date = SOYBEAN_PLANT.get(scenario, SOYBEAN_PLANT[zone])
+        pdate = date.strftime("%m    %d")
+        pdatem5 = (date - timedelta(days=5)).strftime("%m    %d")
+        pdatem10 = (date - timedelta(days=10)).strftime("%m    %d")
+        pdatem15 = (date - timedelta(days=15)).strftime("%m    %d")
+        plant = SOYBEAN[zone]
+    elif code == "W":
+        date = SOYBEAN_PLANT.get(scenario, SOYBEAN_PLANT[zone])  # TODO
+        plant = WHEAT_PLANT[zone]
+        pdate = date.strftime("%m    %d")
+        pdatem5 = (date - timedelta(days=5)).strftime("%m    %d")
+        pdatem10 = (date - timedelta(days=10)).strftime("%m    %d")
+        pdatem15 = (date - timedelta(days=15)).strftime("%m    %d")
+    return data % {
+        "yr": year,
+        "pdate": pdate,
+        "pdatem5": pdatem5,
+        "pdatem10": pdatem10,
+        "pdatem15": pdatem15,
+        "plant": plant,
+    }
diff --git a/tests/test_tillage.py b/tests/test_tillage.py
new file mode 100644
index 00000000..3d19897d
--- /dev/null
+++ b/tests/test_tillage.py
@@ -0,0 +1,25 @@
+"""Test pydep.tillage."""
+
+from pydep import tillage
+
+
+def test_simple():
+    """Test import of API."""
+    plants = ["C", "B", "W", "P"]
+    cfactors = list(range(1, 7))
+    zones = ["IA_NORTH", "IA_SOUTH", "KS_NORTH"]
+    for bp in plants:
+        for pp in plants:
+            for ap in plants:
+                for cfactor in cfactors:
+                    for zone in zones:
+                        res = tillage.make_tillage(
+                            0, zone, bp, pp, ap, cfactor, 1
+                        )
+                        assert res is not None
+
+
+def test_unknown_file():
+    """Test unknown file."""
+    res = tillage.make_tillage(0, "IA_NORTH", "q", "q", "q", 1, 1)
+    assert res == ""