diff --git a/docs/llm/architecture.md b/docs/llm/architecture.md
index e49d697..210e185 100644
--- a/docs/llm/architecture.md
+++ b/docs/llm/architecture.md
@@ -24,7 +24,8 @@ HFCollection {path: netCDFMeta}                    gents/meta.py
 TSCollection — list of order dicts                 gents/timeseries.py
    │  modifiers: include/exclude, add_args, apply_compression,
    │             apply_chunk_target_bytes, apply_path_swap, apply_overwrite,
-   │             append_timestep_dirs, add_attrs
+   │             append_timestep_dirs, add_attrs, skip_existing (resume a
+   │             continued run -- narrows/drops orders against existing output)
    │  .execute() — batches orders sharing source files (optimize=True),
    │               pool (or in-process) → generate_time_series() per batch
    ▼
@@ -46,9 +47,9 @@ bundled YAML config (`gents/configs/*.yaml`) merged with command-line flags.
 |---|---|---|
 | `datastore.py` | `GenTSDataStore` — the *only* place `netCDF4.Dataset` is constructed in the pipeline; context manager + delegation | netCDF4 |
 | `meta.py` | `netCDFMeta` (cached header metadata, with opt-out load flags), `is_var_secondary` (primary/secondary rules), `get_meta_from_path` (picklable factory), `get_attributes`, `get_time_variables_names` | datastore |
-| `hfcollection.py` | `HFCollection` + module-level helpers for discovery, grouping, year math, fragmentation merging | meta, utils |
+| `hfcollection.py` | `HFCollection` + module-level helpers for discovery, grouping, date math (`get_time_boundary_num`, `get_year_boundary_num`), fragmentation merging | meta, utils |
 | `mhfdataset.py` | `MHFDataset` — presents a file group as one virtual dataset; time→file mapping; memory-bounded variable cache; tile reassembly | datastore, meta |
-| `timeseries.py` | `TSCollection`, order construction/execution, `write_timeseries_file`, `compute_chunksizes`, `CHUNK_TARGET_BYTES`, `check_timeseries_integrity`, `check_timeseries_conform`, `get_timestamp_format`, `get_timestep_label` | mhfdataset, datastore, meta, utils |
+| `timeseries.py` | `TSCollection` (including `skip_existing`, the resume modifier), order construction/execution, `write_timeseries_file`, `compute_chunksizes`, `CHUNK_TARGET_BYTES`, `check_timeseries_integrity`, `check_timeseries_conform`, `get_timestamp_format`, `get_timestep_label` | mhfdataset, datastore, meta, utils |
 | `cli.py` | argparse, YAML config loading/merging, `check_config`, `main()` | hfcollection, timeseries, utils |
 | `utils.py` | `enable_logging`, `ProgressBar`, `get_version`, `get_time_stamp`, `log_hfcollection_info`, `log_tscollection_info`, `LOG_LEVEL_IO_WARNING = 5` | — |
 
@@ -160,13 +161,15 @@ number of simultaneously open handles is therefore 1, independent of group size.
 - **Time mapping:** `{float_time: [(file_index, sub_time_index), ...]}`. `sub_time_index`
   is the step's position inside its own file, precomputed so reads never rescan a time
   array. More than one entry per time ⇒ fragmented group.
-- **Cache:** `open()` caches every secondary variable plus as many primaries from
-  `preload_var_list` as fit under `memory_limit_bytes` (`__plan_cacheable_vars`, sized
+- **Cache:** `open()` caches as many secondaries, then as many primaries from
+  `preload_var_list`, as fit under `memory_limit_bytes` (`__plan_cacheable_vars`, sized
   from one file's shape × file count). All-or-nothing per variable — a partial entry
   would serve the wrong file's data. Anything left over is cached on first use if it
   fits, otherwise only the requested time slice is read from disk per access.
   `preload_primaries=False` (used by `no_data` runs) skips all primary preloading —
-  primaries stay readable on demand.
+  primaries stay readable on demand. Every read path uses raw values
+  (`set_auto_maskandscale(False)`), so cached data is plain `ndarray` whose size the
+  budget counts exactly.
 - **Eviction:** reads move through one variable at a time, so switching variables frees
   the previous one's cache. A single file's entry is *not* freed after one read, since
   write chunks don't align with source file boundaries.
@@ -205,16 +208,22 @@ workers are logged (`raise_errors=False` default) rather than raised, so a bad f
 doesn't kill a long batch job. Dask was removed; `dask_client` kwargs remain only as
 deprecation shims.
 
+Neither pool sets a start method, so workers start with the platform default: `fork`
+on Linux before Python 3.14 (including the Docker images), `forkserver` from 3.14.
+
 Memory is bounded per worker, not globally: `execute(memory_limit_bytes=...)` (CLI
-`--memory-limit`, in GB) is forwarded to every `MHFDataset` a worker opens, so the
-process-wide ceiling is roughly `tscores × memory_limit`. The default is
-`DEFAULT_MEMORY_LIMIT_BYTES` (4 GiB per worker); constructing an `MHFDataset` directly
-is still unbounded by default.
+`--memory-limit`, in GB) is forwarded to every `MHFDataset` a worker opens and bounds
+that worker's cache. The parent is not bounded: it holds every file's `netCDFMeta`
+for the whole run (5.8 GiB on a 40k-file CESM3 case), and forked workers can copy part
+of it. See the `MHFDataset` cache invariants and gotchas in
+[conventions.md](conventions.md). The default is `DEFAULT_MEMORY_LIMIT_BYTES` (4 GiB
+per worker); constructing an `MHFDataset` directly is still unbounded by default.
 
 ## Error-handling philosophy
 
-Skip-and-warn, not fail-fast: invalid files are dropped by `check_validity`, minority
-variable sets are dropped by `check_groups_by_variables`, worker exceptions are logged.
+Skip-and-warn, not fail-fast: invalid files are dropped by `check_validity` and worker
+exceptions are logged. Variable-set consistency within a group is *not* checked
+(`check_groups_by_variables` exists but is never called).
 Logging uses the `"gents"` logger hierarchy; `enable_logging(verbose=True)` turns on
 per-file I/O traces at custom level `LOG_LEVEL_IO_WARNING = 5`.
 
diff --git a/docs/llm/concepts.md b/docs/llm/concepts.md
index 61cfe72..e288bcc 100644
--- a/docs/llm/concepts.md
+++ b/docs/llm/concepts.md
@@ -16,7 +16,8 @@
   `h0`, `h1`, `h2`, ... (CESM) or `eam.h0`... (E3SM). Different streams have different
   variable sets and time step sizes, so they must never be mixed in one output file.
 - **Model component:** Subdirectory per model subsystem in a case's output tree:
-  `atm`, `ocn`, `lnd`, `ice`, `glc`, `rof`. Conventionally raw files live under
+  `atm`, `ocn`, `lnd`, `ice`, `glc`, `rof` (CESM3 cases also have `cpl` and `wav`,
+  which the bundled config's `"*.nc"` include pulls in). Conventionally raw files live under
   `<component>/hist/` and processed output under `<component>/proc/tseries/`.
 - **Case:** One model simulation/run and its output directory tree.
 
@@ -61,6 +62,25 @@
 - **Multistep slicing:** A single HF can contain many time steps that straddle a slice
   boundary; `HFCollection` records per-file `(start_index, end_index)` tuples
   (`get_multistep_slices`) so `TSCollection` can split within a file.
+- **Resume / incremental generation (`TSCollection.skip_existing`):** Narrows each
+  already-built order to the HF source files not yet covered by existing output, for a
+  CESM/E3SM run that was continued from a restart after time series were already generated
+  for the earlier portion. Scans the output directory (its own directory only, not
+  recursively) for files matching `{ts_path_template}.{primary_var}.*.nc`, and takes the
+  latest raw time any *complete* (`check_timeseries_integrity`-passing) match covers as the
+  cutoff. Drops each `hf_paths` entry *whole* if its own last raw time step does not exceed
+  that cutoff, keeps it whole otherwise — **file granularity, not step granularity**: a
+  continued run always starts a fresh HF file at the restart boundary (a stream's step count
+  can't change mid-run without a namelist edit and a fresh run), so no single file can
+  straddle "already covered" and "genuinely new," and there is nothing to split within one.
+  `ts_start_index`/`ts_end_index` are therefore left untouched by this modifier; the
+  trimmed order's `ts_string` is recomputed from the first surviving file's own first
+  timestep so the output name still describes what it actually covers. An order with no
+  existing output is untouched; one already covered in full is dropped. Per-order
+  granularity means two variables in the same group can have independently different resume
+  points. Does not detect or fill a gap between two disjoint existing outputs, and does not
+  warn about one either — see `conventions.md`. Library-only for now: nothing in `cli.py`
+  or the bundled YAML configs calls it.
 - **Fragmented (tiled) files:** Some models split one time step *spatially* across
   several files (e.g. per latitude band); detected by paths not ending in `.nc` and
   merged by matching non-time dimension bounds (`merge_fragmented_groups`).
@@ -69,15 +89,17 @@
 - **Order:** A plain dict describing one TS output file to generate — the unit of work
   of `TSCollection`. See [architecture.md](architecture.md) for the exact schema.
 - **Variable cache / memory limit:** `MHFDataset` reads a group's files one at a time and
-  keeps *data* in memory rather than keeping *handles* open. On `open()` it caches every
-  secondary variable plus as many primaries as fit under `memory_limit_bytes`
+  keeps *data* in memory rather than keeping *handles* open. On `open()` it caches as
+  many secondaries, then primaries, as fit under `memory_limit_bytes`
   (`execute(memory_limit_bytes=...)`, CLI `--memory-limit`, in GB); the rest are read on
   demand and cached if they fit. Caching is all-or-nothing per variable and per group, and
   a variable's cache is released when reads move on to the next variable. The limit is
-  per worker process, so the pipeline-wide ceiling is roughly `tscores × limit`; the
-  default is `timeseries.DEFAULT_MEMORY_LIMIT_BYTES` (4 GiB per worker) — a bare
-  `MHFDataset` constructed directly is still unbounded. Data too big to cache is read
-  from disk per run of consecutive steps, fetching only the needed slice.
+  per worker process and covers only the cache: the pipeline-wide peak is roughly
+  `parent metadata + tscores × (limit + overhead)`. The default is
+  `timeseries.DEFAULT_MEMORY_LIMIT_BYTES` (4 GiB per worker) — a bare `MHFDataset`
+  constructed directly is still unbounded. Data too big to cache is read from disk per
+  run of consecutive steps, fetching only the needed slice. All reads are raw (no
+  masking, no `scale_factor` applied), so packed data is copied as stored.
 - **Lazy date decoding:** `cftime.num2date` is measurably expensive per file, and the
   pipeline works almost entirely on raw float time values, which order identically to
   their decoded dates within one `(units, calendar)` reference. `get_meta_from_path`
@@ -85,15 +107,18 @@
   grouping, slicing, timestamp ranges) stays in the float domain and only endpoint
   values are decoded, via `netCDFMeta.decode_time_values` /
   `decode_time_bounds_values`. Year-window membership is tested against per-reference
-  float boundaries from `hfcollection.get_year_boundary_num`. Cross-file comparisons use
-  the decoded CFTime endpoints, so mixed time references between files stay correct.
+  float boundaries from `hfcollection.get_year_boundary_num`. Cross-file comparisons in
+  `HFCollection`'s year/slice/timestep math use decoded CFTime endpoints, so mixed time
+  references stay correct there. Within a group, `sort_along_time`, `MHFDataset` and
+  the written time axis still assume one reference (see `conventions.md`).
   Sibling flags (`load_time_bounds`, `load_variable_attrs`, `compute_dim_bounds`) skip
   other per-file reads; each of those getters raises rather than lying when its data was
   skipped, while `get_cftimes()`/`get_cftime_bounds()` simply decode on first call.
 - **Integrity stamp:** Every completed TS file gets a `gents_version` global attribute.
   Dual use: (1) output files lacking it are considered corrupt/partial and are
-  regenerated; (2) *input* files carrying it are recognized as GenTS output and
-  excluded from processing (`netCDFMeta.is_valid`).
+  regenerated, while stamped ones are not rewritten (the check happens at write time,
+  after the group's input has already been read); (2) *input* files carrying it are
+  recognized as GenTS output and excluded from processing (`netCDFMeta.is_valid`).
 - **Conforming chunking:** Output chunking convention checked by
   `check_timeseries_conform`: `time` stored contiguously; large variables chunked along
   time so each chunk is ≈ 4 MiB (CMOR-friendly). Files < 4 MiB are stored contiguously.
diff --git a/docs/llm/conventions.md b/docs/llm/conventions.md
index c1439b7..68e56c9 100644
--- a/docs/llm/conventions.md
+++ b/docs/llm/conventions.md
@@ -2,7 +2,7 @@
 
 > Rules the codebase relies on that are easy to violate, plus known stale spots.
 > Verify a gotcha still exists before acting on it — this file describes the code as of
-> August 2026 (branch `validation`).
+> September 2026 (branch `main`).
 
 ## API invariants — do not break
 
@@ -17,7 +17,9 @@
   filter first, read headers second). Methods needing metadata call `check_pulled()`,
   which auto-pulls.
 - **Copies inherit metadata.** Filtering a pulled collection yields pulled
-  sub-collections; never re-read headers that a parent already read.
+  sub-collections; never re-read headers that a parent already read. (They also inherit
+  the parent's *groups*, which is a bug — see "`include()`/`exclude()` after a pull"
+  below.)
 - **Worker picklability.** Anything submitted to `ProcessPoolExecutor` must be a
   module-level function with picklable args. `get_meta_from_path` exists precisely as
   the picklable factory for `netCDFMeta`; keep that pattern.
@@ -46,6 +48,41 @@
 - **`[sorting_pivot]` group-key suffix** is the string protocol between
   `HFCollection.slice_groups` and `TSCollection.update_ts_orders`. Both sides parse it
   literally; change it in both places or nowhere.
+- **`include_time`/`include_years` half-open range.** `include_time(start_year, end_year, ...)`
+  keeps `start_date <= representative_time < end_date` (end **exclusive**), matching the
+  `[lower_num, upper_num)` convention `slice_groups` already used. `include_years` is a
+  thin wrapper (`include_time(start_year, end_year + 1)`); do not "fix" it to inclusive
+  bounds without checking both together. Only calendar-agnostic `year`/`month`/`day`
+  components are ever taken from a caller — there is no finer granularity (hour/minute/
+  second), and no date-like object is accepted, specifically so the boundary can be
+  re-expressed fresh in any file's own calendar via `get_time_boundary_num` without ever
+  comparing two `cftime.datetime` objects that might carry different calendars (that raises
+  `TypeError`; see the float-domain time handling section below).
+- **`TSCollection.skip_existing()` strict inequality, whole-file granularity.** The resume
+  cutoff compares with `>`, not `>=` (a continuation run can re-emit the prior run's exact
+  last timestep as its own first timestep, and `>` is what drops that duplicate rather than
+  keeping or double-counting it), and a source file is kept or dropped *whole* rather than
+  split. That's deliberate, not a missed case: an output stream's step count cannot change
+  mid-run without a namelist edit and a fresh run, so a continuation always starts a new HF
+  file at the restart boundary and no single file can straddle "already covered" and
+  "genuinely new." Do not "simplify" the inequality to `>=`, and do not reintroduce
+  `ts_start_index`/`ts_end_index` narrowing to handle a straddling file — it cannot happen
+  in this workflow.
+- **`skip_existing()` compares raw floats across files.** It compares an existing TS
+  output's raw time value directly against a candidate HF file's raw time value, with no
+  `decode_time_values()` call — a scoped exception to the float-domain rule below, which
+  otherwise treats an undecoded cross-file comparison as unsafe. It relies on both sides
+  (the order's own output and its own input) sharing one `(units, calendar)` reference.
+  That is usual but **not guaranteed** within a stream: WW3 writes a new `units` string per
+  restart segment (see the WW3 gotcha below). There the raw values happen to stay
+  continuous, so the comparison still works, but only by accident. Do not generalize this
+  pattern to a comparison spanning two different streams.
+- **`skip_existing()` does not fill a gap, and does not warn about one either.** If more
+  than one existing, valid time series file is found for one order, only the *latest*
+  covered end is used as the cutoff — a gap between two disjoint existing outputs (e.g.
+  years 0-3 and 6-7 exist, 4-5 does not) is neither detected, backfilled, nor logged. This
+  is a known, accepted limitation, not a bug to silently "fix" by unioning ranges — see the
+  resume concept in `concepts.md`.
 - **All user-facing filters are `fnmatch` globs** applied to *absolute path strings*
   (or variable names for `var_glob`). Not regex, not `pathlib.match`. Every glob
   parameter takes **either a single string or a list of them**, normalised at the top
@@ -72,22 +109,35 @@
   array per file in the group, or the variable is absent from the cache entirely. A
   partial (prefix-only) list makes `__get_hf_data` return a *different file's* data by
   index, silently.
+- **Secondaries are budgeted too.** `open()` plans secondaries before primaries, and
+  both go through `__plan_cacheable_vars`, so a secondary that doesn't fit (e.g. CICE's
+  static grid × 3,650 daily files) is read on demand, not cached. Its size estimate is
+  one file's shape × file count even for variables with no time dimension.
 - **One file open at a time.** `open()` and `__cache_variable` both use
   `with GenTSDataStore(path)` per file. Nothing in the class holds a handle open across
   files, and nothing should — that property is what retired the `EMFILE` gotcha.
+- **Every read is raw.** Each handle `MHFDataset` opens calls
+  `set_auto_maskandscale(False)` before reading. Masked arrays carry a mask that
+  `__get_cache_dsize` doesn't count (+25% over `memory_limit_bytes` for float32), and
+  scaled reads of packed ints get truncated back to the on-disk dtype by `get_var_vals`.
+  A new read site must do the same.
 - **Eviction is per-variable, on switch.** Do not free a file's entry after reading it:
   write chunks don't align with source file boundaries, so the same entry is legitimately
   read more than once per variable.
-- **`memory_limit_bytes` is per `MHFDataset`**, i.e. per worker process. The pipeline-wide
-  ceiling is roughly `tscores × memory_limit`. `execute()` and the CLI default to
-  `timeseries.DEFAULT_MEMORY_LIMIT_BYTES` (4 GiB); only a directly-constructed
+- **`memory_limit_bytes` is per `MHFDataset`**, i.e. per worker process, and bounds only
+  the cache. The pipeline-wide peak is roughly
+  `parent metadata + tscores × (memory_limit + worker overhead)`. The parent holds every
+  file's `netCDFMeta` (measured at 5.8 GiB for a 40k-file CESM3 case), and under `fork`
+  workers can end up copying part of it (see the gotchas below). `execute()` and the CLI
+  default to `timeseries.DEFAULT_MEMORY_LIMIT_BYTES` (4 GiB); only a directly-constructed
   `MHFDataset` is unbounded by default.
 
 ### `netCDFMeta` load flags
 
 `decode_dates`, `load_time_bounds`, `load_variable_attrs`, `compute_dim_bounds` exist to
 skip measured per-file cost for callers that don't need the result (`MHFDataset` needs
-none of them; `get_meta_from_path` — the whole pipeline — opts out of `decode_dates`).
+only the first file's variable attrs; `get_meta_from_path` — the whole pipeline — opts
+out of `decode_dates` only, so it still loads attrs nobody reads — see gotchas).
 Opting out of the latter three must stay **strict**: the corresponding getter raises
 `RuntimeError` rather than returning empty or stale data. Don't "helpfully" soften that —
 silent empties would surface as wrong output, not an error. `decode_dates=False` is the
@@ -108,6 +158,12 @@ did). Year-window membership is tested against per-reference boundaries from
 `num2date` — the regression is guarded by size-capped `num2date` assertions in
 `test_hfcollection.py` and `test_tscollection.py`.
 
+Reference-safety covers the `HFCollection` year/slice/timestep math only. Several
+places still compare raw floats across files and assume **one reference per group**:
+`sort_along_time` (sorts all files by raw first time), `MHFDataset`'s time mapping, and
+the written time axis (raw values from every file, labelled with the first file's
+`units`). A group that mixes references gets a wrong output time axis.
+
 ## Conformity conventions (`gents/conformity/`) — do not break
 
 `gents/conformity/README.md` is authoritative; these are the rules most likely to be
@@ -154,7 +210,45 @@ violated by a well-meaning refactor.
 - Private attributes use double-underscore name mangling (`self.__hf_to_meta_map`).
 - Loggers are `logging.getLogger(__name__)` under the `"gents"` hierarchy.
 
-## Known gotchas / stale spots (verified 2026-08)
+## Known gotchas / stale spots (verified 2026-09)
+
+- **Static secondaries that don't fit the cache are truncated to their first row.**
+  `__cache_variable` stores a no-time variable as a bare array
+  (`__data_var_cache[var] = hf_ds[var][:]`) instead of a one-item list, so
+  `__get_hf_data(0, var)` returns `cache[var][0]`, its first row. `write_timeseries_file`
+  then broadcasts that row across the whole output variable. It triggers whenever a static
+  secondary misses the `open()` budget (every daily CICE stream: `TLAT`, `tarea`,
+  `*_bounds`, …). Missing-value clones hide it because every value is fill. Repro: 3
+  files, a static `(4, 5)` variable, `memory_limit_bytes=200` →
+  `get_var_vals` returns shape `(5,)`.
+- **`include()`/`exclude()` after a pull keep stale groups.** `copy()` inherits
+  `self.get_groups()` whenever the collection is pulled, so a filtered copy still lists
+  every group of its parent. `slice_groups()` then raises `KeyError` on the first dropped
+  path (`hfc.pull_metadata(); hfc.include("*ww3.hi*").slice_groups()`). Filtering *before*
+  the pull (as `cli.main` does) avoids it.
+- **Worker pools use the platform's default start method.** Neither pool passes
+  `mp_context`. That's `fork` on Linux before Python 3.14, including the Docker images
+  (`DEFAULT_PYTHON=3.12`). Forked workers inherit the parent's `netCDFMeta` heap, and
+  CPython's refcount/GC writes gradually copy those pages. Measured with a 5.8 GiB parent:
+  one `--no-data` execute worker reached 2.4 GiB private memory under `fork` vs 440 MiB
+  under `forkserver`.
+- **The parent keeps per-variable attrs it never reads.** `get_meta_from_path` builds
+  metas with `load_variable_attrs=True`, but nothing reads attrs from an `HFCollection`'s
+  metas (`MHFDataset` builds its own). It is 53–76% of each meta's size on CESM3 streams
+  (`cpl.hi`: 58.7 → 27.4 KiB pickled).
+- **Re-running a finished case still reads every group.** The integrity-stamp skip lives
+  inside `write_timeseries_file`, *after* `generate_time_series` has opened the
+  `MHFDataset` (one pass over every file, preloading primaries up to the limit) and read
+  the secondaries. Output is not rewritten, but the input is read again. Only
+  `skip_existing()` drops work before any file is opened.
+- **Dead code in `hfcollection.py`.** `check_groups_by_variables`, `filter_by_variables`,
+  `sort_metas_by_time` and `is_ds_within_years` are defined but never called. Nothing
+  enforces a consistent variable set within a group.
+- **WW3 history files have self-inconsistent time units.** Each restart segment's `time`
+  counts seconds from model start (0001), but its `units` say `seconds since <restart
+  date>`. A file named `ww3.hi.0005-01-02` therefore decodes to `0009-01-02`, and a
+  24-year run slices as years 1–50. GenTS decodes each file faithfully; the input is
+  wrong. Other CESM3 components keep one `days since 0001-01-01` throughout.
 
 - **`gents/conformity/` is not packaged.** `pyproject.toml` sets
   `packages = ["gents", "gents.configs"]`, an explicit list, so a built wheel/sdist
@@ -180,13 +274,6 @@ violated by a well-meaning refactor.
   returns a collection containing only `TREFHT` — every other order is silently discarded
   rather than left un-prefixed. Harmless at the default `"*"`, which is the only way it is
   currently called.
-- **`append_timestep_dirs` is dead config in the CLI.** `gents_cesm3.yaml` and
-  `gents_example.yaml` both set `output_ts.append_timestep_dirs: true`, but `cli.main`
-  only reads `path_swaps` and `compression` from `output_ts` — `TSCollection.append_timestep_dirs()`
-  is never called anywhere outside its own definition and its tests. So CLI output has no
-  `month_1/`-style frequency directories despite the config asking for them, and any
-  CLI-driven output fails the conformity check for it. Either wire the key up in
-  `cli.main` or drop it from the YAML.
 - **`calculate_year_slices` quirks:** the guard's error message is inverted
   ("Maximum year cannot exceed minimum year" fires when max < min), and the early
   return triggers when `slice_size_years >= max_year - min_year`, so a span exactly
@@ -244,9 +331,10 @@ violated by a well-meaning refactor.
   lookup over a `GenTSDataStore` but returns empty-handed where `netCDFMeta` raises. Keep
   that split: the clone must describe every file it copies, not just the GenTS-legible
   ones.
-- **`MHFDataset` trusts the first file** of a group for variable dims/dtype/attrs;
-  variable-set consistency is enforced earlier by `check_groups_by_variables`
-  (majority wins, minority files dropped with a warning).
+- **`MHFDataset` trusts the first file** of a group for variable dims/dtype/attrs, and
+  nothing upstream enforces a consistent variable set (`check_groups_by_variables` is
+  dead code, see above). A secondary missing from a later file is skipped at `open()`;
+  a primary missing from a later file fails that worker call when it is read.
 - **Time-bounds variable names are matched case-insensitively** against exactly:
   `time_bnds`, `time_bnd`, `time_bounds`, `time_bound`; the time variable must have
   `units` and `calendar` attributes or the file is rejected.
@@ -277,7 +365,15 @@ violated by a well-meaning refactor.
 - **`EMFILE` on wide streams.** `MHFDataset` used to open every file of a group at once,
   so a ~1800-file daily stream died with `OSError: [Errno 24] Too many open files`. It now
   opens one file at a time and caches data instead of handles, so the ulimit workaround is
-  no longer needed. `gents/conformity/README.md` still carries the old warning.
+  no longer needed.
+- **`append_timestep_dirs` was dead config in the CLI.** Both bundled YAMLs set
+  `output_ts.append_timestep_dirs: true` but `cli.main` never read it. It is now wired up
+  (`cli.main` calls `tsc.append_timestep_dirs()` when the key is true).
+- **Lazy `MHFDataset` reads were masked and scaled.** `__cache_variable` and the
+  direct-read fallback opened handles without `set_auto_maskandscale(False)`, unlike
+  `open()`. The uncounted masks pushed the cache ~25% over `memory_limit_bytes`, and
+  scaled reads of packed ints were truncated back to the on-disk dtype. Fixed in
+  `ea8f36f`; see "Every read is raw" above.
 - **`hfcollection.check_config` / `get_default_config`.** Removed (2026-08): they asserted a
   `{name, include, exclude}` config shape that no version of GenTS still uses, were called
   from nowhere, and collided by name with the live `gents.cli.check_config`. The YAML
diff --git a/docs/llm/index.md b/docs/llm/index.md
index 5f48dcd..b02b368 100644
--- a/docs/llm/index.md
+++ b/docs/llm/index.md
@@ -41,7 +41,7 @@ verified against the source; when docs and code disagree, trust the code and fix
 ## Repository map
 
 ```
-gents/                  Package source (7 pipeline modules ~3,100 lines, + conformity/)
+gents/                  Package source (7 pipeline modules ~3,700 lines, + conformity/)
   hfcollection.py       HFCollection: discover/filter/group/slice history files
   timeseries.py         TSCollection: build + execute time-series "orders"; file writing
   meta.py               netCDFMeta: cached per-file metadata; primary/secondary classification
@@ -56,13 +56,14 @@ gents/                  Package source (7 pipeline modules ~3,100 lines, + confo
     check.py              gents_conform — run a model spec against generated output
     report.py             pass/fail/skip collector; text + JSON rendering
     models/cesm3.py       what correct CESM3 output looks like (the file that matters)
-  tests/                pytest suite (242 tests); test_cases.py generates synthetic fixtures
+  tests/                pytest suite (284 tests); test_cases.py generates synthetic fixtures
 docs/                   Sphinx docs (index/install/user/dev/api .rst)
 benchmarks/             ASV performance benchmarks + fixtures.build_bench_case()
 pipeline_bench.py       Standalone py-spy profiling driver (repo root, not an ASV suite)
 Dockerfile              Multi-stage: runtime / test / bench / dev / deptest-floor / deptest-latest
 .github/workflows/      tests.yml (docker --target test), dependency_checks, release
 build/lib/              STALE build artifact copy of the package — never read or edit
+downloads/              gitignored; local copies of user-reported cases (e.g. clone tarballs)
 ```
 
 ## The pipeline in one paragraph
@@ -81,7 +82,7 @@ group once into a memory-bounded cache) and writing one netCDF file per variable
 
 ```bash
 pip install -e ".[dev]"          # editable install with test/docs/bench extras
-pytest gents/tests/              # full test suite (242 tests, ~25 s)
+pytest gents/tests/              # full test suite (284 tests, ~35 s)
 docker build --target test -t gents-tests . && docker run --rm gents-tests   # CI-equivalent
 run_gents <hf_dir> --model CESM3 --dryrun   # CLI dry run (read-only validation)
 ```
diff --git a/docs/llm/workflows.md b/docs/llm/workflows.md
index 2c3eb54..0d1a409 100644
--- a/docs/llm/workflows.md
+++ b/docs/llm/workflows.md
@@ -1,7 +1,7 @@
 # Workflows & Recipes
 
 > Commands and step-by-step recipes for common tasks. Commands verified against the
-> Dockerfile, CI workflow, and pyproject as of August 2026.
+> Dockerfile, CI workflow, and pyproject as of September 2026.
 
 ## Environment setup
 
@@ -15,7 +15,7 @@ pip install -r requirements.txt && pip install -e .
 ## Testing
 
 ```bash
-pytest gents/tests/                         # full suite (242 tests, ~25 s)
+pytest gents/tests/                         # full suite (284 tests, ~35 s)
 pytest gents/tests/test_workflow.py         # one file
 ```
 
@@ -93,17 +93,16 @@ asv run --machine gents-container HEAD^!
 
 - `--machine <fixed name>` — asv defaults the machine name to the hostname, which is a
   fresh random hex string in every container, so results never accumulate or compare.
-- `asv.conf.json` lists `branches: ["main", "validation"]` — a working/feature branch
-  needs to be listed (or use an explicit range like `HEAD^!`) for `asv run` to discover
-  its commits.
-
-Two things that used to require workarounds are fixed in the tree, not just documented
-around: `asv.conf.json` no longer pins `pythons` (asv defaults to the interpreter it is
-running under, so any image works out of the box), and the `bench` Dockerfile stage sets
-`ENV HOME=/usr/local/gents`, so `asv machine` no longer fails with `PermissionError` on
-`/.asv-machine.json`. Without a bind mount the image's baked-in `/usr/local/gents` is
-root-owned from `COPY` and still isn't writable by uid 1000 — a separate, pre-existing
-issue.
+- `asv.conf.json` lists `branches: ["main"]` — a working/feature branch needs an
+  explicit range like `HEAD^!` for `asv run` to discover its commits.
+- `asv.conf.json` pins `pythons: ["3.14"]`, but the `bench`/`dev` images default to
+  `DEFAULT_PYTHON=3.12`. Build them with `--build-arg DEFAULT_PYTHON=3.14` so asv finds
+  the pinned interpreter.
+
+The `bench` Dockerfile stage sets `ENV HOME=/usr/local/gents`, so `asv machine` no
+longer fails with `PermissionError` on `/.asv-machine.json`. Without a bind mount the
+image's baked-in `/usr/local/gents` is root-owned from `COPY` and still isn't writable
+by uid 1000 — a separate, pre-existing issue.
 
 ## Performance profiling (py-spy + `pipeline_bench.py`)
 
@@ -231,8 +230,10 @@ mark a group twice). `--slice_start_year` overrides `start_year` in every batch
 way. `--compression` requires `--level`.
 `--model` is case-insensitive; omitted → `gents_example.yaml`. Output dir defaults to
 the input dir (path swaps like `/hist/` → `/proc/tseries/` come from the YAML config).
-`--memory-limit` is in GB and applies **per TS worker**, so the process-wide ceiling is
-about `tscores ×` that; defaults to 4 GB per worker (`DEFAULT_MEMORY_LIMIT_BYTES`). `-nd/--no-data` (→ `execute(no_data=True)`)
+`--memory-limit` is in GB and caps each **TS worker's cache**; defaults to 4 GB per worker
+(`DEFAULT_MEMORY_LIMIT_BYTES`). Size a job for the parent's metadata (it grows with file
+count; 5.8 GiB for a 40k-file CESM3 case, `--dryrun` alone shows it) plus
+`--tscores ×` a bit over the limit. `-nd/--no-data` (→ `execute(no_data=True)`)
 builds the full directory/file structure but skips reading/writing primary-variable data
 — primaries read back as their fill value; the fast path for conformity runs over
 missing-value clones (see the `--no-data` concept in [concepts.md](concepts.md)). Don't
@@ -357,8 +358,7 @@ first — several obvious-looking refactors are explicitly unwanted there.
 Required top-level keys (asserted by `cli.check_config`): `version`, `model`,
 `input_hf`, `output_ts`. Under `input_hf`: `match` (discovery glob), `include`,
 `exclude`, `slicing` (list of `slice_groups` kwarg dicts). Under `output_ts`:
-`append_timestep_dirs` (currently ignored by the CLI — see
-[conventions.md](conventions.md)), `path_swaps` (list of `apply_path_swap` kwargs),
+`append_timestep_dirs` (bool → `append_timestep_dirs()`), `path_swaps` (list of `apply_path_swap` kwargs),
 `compression` (list of `apply_compression` kwargs). To add support for a new model, add a
 YAML here, register it in `model_config_files` in `cli.main`, include it in
 `[tool.setuptools.package-data]` (already `*.yaml`), and add tests mirroring
@@ -384,10 +384,13 @@ tsc = tsc.append_timestep_dirs()
 output_paths = tsc.execute(memory_limit_bytes=8 * 1024**3)     # returns written file paths
 ```
 
-Filter *before* `pull_metadata()`/`TSCollection` construction to avoid needless header
-reads. Re-running is cheap: existing complete outputs (integrity-stamped) are skipped
-unless `apply_overwrite` / `--overwrite` is set. Pass `show_progress=False` to
-`pull_metadata`/`execute` when driving GenTS from a script that owns stdout.
+Filter *before* `pull_metadata()`/`TSCollection` construction, both to avoid needless
+header reads and because filtering a pulled collection keeps its stale groups (see
+[conventions.md](conventions.md)). Re-running doesn't rewrite complete
+(integrity-stamped) outputs unless `apply_overwrite` / `--overwrite` is set, but every
+group's input is still read; use `skip_existing()` to drop finished work up front. Pass
+`show_progress=False` to `pull_metadata`/`execute` when driving GenTS from a script that
+owns stdout.
 
 ### Inspecting state interactively
 
diff --git a/gents/cli.py b/gents/cli.py
index ec923e6..0afa546 100644
--- a/gents/cli.py
+++ b/gents/cli.py
@@ -86,7 +86,7 @@ def parse_arguments():
         "--slice_start_year",
         type=int,
         default=None,
-        help="Year to start slice windows at. (Default is start year for history files)"
+        help="Year to align slice windows to; earlier data gets windows on the same alignment. (Default is start year for history files)"
     )
     parser.add_argument(
         "-hc", "--hfcores",
@@ -134,8 +134,10 @@ def parse_arguments():
     parser.add_argument(
         "--align_method",
         type=str,
-        default="midpoint",
-        help="Method to use when aligning the history files by time. ('midpoint', 'direct_time', 'start_bound', 'end_bound')"
+        default=None,
+        choices=["midpoint", "direct_time", "start_bound", "end_bound"],
+        help="Method to use when aligning the history files by time, for slicing and "
+             "output file names. Overrides the model configuration. (Default midpoint)"
     )
     parser.add_argument(
         "--compression",
@@ -247,12 +249,16 @@ def main():
     for slice_batch in slice_batches:
         if args.slice_start_year is not None:
             slice_batch["start_year"] = args.slice_start_year
+        if args.align_method is not None:
+            slice_batch["time_alignment_method"] = args.align_method
         hfc = hfc.slice_groups(**slice_batch)
 
     tsc = TSCollection(hfc, args.outputdir, num_processes=args.tscores)
-
-    if args.compression is not None:
-        tsc = tsc.apply_compression(alg=args.compression, level=args.level, path_glob="*")
+    if args.align_method is not None:
+        # Output file names use the same alignment as the slicing.
+        tsc = tsc.update_ts_orders(time_alignment_method=args.align_method)
+    if args.overwrite:
+        tsc = tsc.apply_overwrite(path_glob="*")
 
     if "path_swaps" in yaml_config["output_ts"]:
         for swap_batch in yaml_config["output_ts"]["path_swaps"]:
@@ -266,6 +272,10 @@ def main():
         for comp_batch in yaml_config["output_ts"]["compression"]:
             tsc = tsc.apply_compression(**comp_batch)
 
+    # After the YAML's, so an explicit --compression wins over the model default.
+    if args.compression is not None:
+        tsc = tsc.apply_compression(alg=args.compression, level=args.level, path_glob="*")
+
     if args.verbose:
         log_hfcollection_info(hfc)
         log_tscollection_info(tsc)
diff --git a/gents/hfcollection.py b/gents/hfcollection.py
index b97bd93..e5ae0de 100644
--- a/gents/hfcollection.py
+++ b/gents/hfcollection.py
@@ -63,7 +63,7 @@ def calculate_year_slices(slice_size_years, min_year, max_year):
     """
     if max_year < min_year:
         raise ValueError("Maximum year cannot exceed minimum year.")
-    if slice_size_years >= max_year - min_year:
+    if max_year - min_year < slice_size_years:
         return [(min_year, max_year)]
 
     start_year = min_year
@@ -162,34 +162,61 @@ def get_year_boundary_num(year, units, calendar):
     return get_time_boundary_num(units, calendar, year)
 
 
-def get_year_bounds(hf_to_meta_map):
+def get_aligned_times(meta, time_alignment_method="midpoint"):
+    """
+    Returns one file's representative time for each step, with the reference
+    needed to decode them.
+
+    ``'midpoint'``, ``'start_bound'`` and ``'end_bound'`` pick from the time bounds;
+    ``'direct_time'``, or a file without bounds, uses the time values themselves.
+
+    :param meta: Metadata of the history file.
+    :type meta: gents.meta.netCDFMeta
+    :param time_alignment_method: ``'midpoint'``, ``'start_bound'``,
+        ``'end_bound'`` or ``'direct_time'``.
+    :type time_alignment_method: str
+    :returns: ``(times, units, calendar)``.
+    :rtype: tuple[numpy.ndarray, str, str]
+    :raises ValueError: If ``time_alignment_method`` is not one of those four.
+    """
+    float_bnds = meta.get_float_time_bounds()
+    if float_bnds is None or time_alignment_method == "direct_time":
+        times = np.ma.getdata(np.atleast_1d(meta.get_float_times()))
+        return times, meta.get_time_units(), meta.get_time_calendar()
+
+    bnds = np.ma.getdata(float_bnds)
+    if time_alignment_method == "midpoint":
+        times = bnds[:, 0] + (bnds[:, 1] - bnds[:, 0]) / 2
+    elif time_alignment_method == "start_bound":
+        times = bnds[:, 0]
+    elif time_alignment_method == "end_bound":
+        times = bnds[:, 1]
+    else:
+        raise ValueError(f"'{time_alignment_method}' is an invalid time-alignment method. Valid methods are ['direct_time', 'midpoint', 'start_bound', 'end_bound']")
+    return times, meta.get_time_bounds_units(), meta.get_time_bounds_calendar()
+
+
+def get_year_bounds(hf_to_meta_map, time_alignment_method="midpoint"):
     """
     Returns the ``(min_year, max_year)`` covered by a set of history files.
 
-    A file's year comes from the midpoint of each time bound, or from the time
-    value itself when the file has no bounds.
+    Years come from each step's aligned time (see :func:`get_aligned_times`), the
+    same times :meth:`HFCollection.slice_groups` assigns files to windows by.
 
     :param hf_to_meta_map: ``{path: netCDFMeta}`` mapping to inspect.
     :type hf_to_meta_map: dict
+    :param time_alignment_method: See :func:`get_aligned_times`.
+    :type time_alignment_method: str
     :rtype: tuple[int, int]
     """
     min_year = np.inf
     max_year = -np.inf
 
-    for path in list(hf_to_meta_map.keys()):
-        meta = hf_to_meta_map[path]
-        float_bounds = meta.get_float_time_bounds()
-        # Midpoints are computed on the raw values, which order identically to
-        # their decoded dates, so only the two extremes need decoding per file.
-        if float_bounds is None:
-            midpoints = np.ma.getdata(np.atleast_1d(meta.get_float_times()))
-            decode = meta.decode_time_values
-        else:
-            bounds = np.ma.getdata(float_bounds)
-            midpoints = bounds[:, 0] + (bounds[:, 1] - bounds[:, 0]) / 2
-            decode = meta.decode_time_bounds_values
-
-        extremes = np.atleast_1d(decode(np.array([np.min(midpoints), np.max(midpoints)])))
+    for meta in hf_to_meta_map.values():
+        # Raw values order identically to their decoded dates, so only the two
+        # extremes need decoding per file.
+        times, units, calendar = get_aligned_times(meta, time_alignment_method)
+        extremes = np.atleast_1d(num2date(np.array([np.min(times), np.max(times)]), units=units, calendar=calendar))
         if extremes[-1].year > max_year:
             max_year = extremes[-1].year
         if extremes[0].year < min_year:
@@ -390,8 +417,9 @@ def merge_fragmented_groups(hf_groups, hf_meta_map):
 
     dim_hashes = {}
     for pattern in fragmented_groups:
-        dims = hf_meta_map[fragmented_groups[pattern][0]].get_dim_bounds()
-        dims = {variable: dims[variable] for variable in dims if variable != "time"}
+        init_meta = hf_meta_map[fragmented_groups[pattern][0]]
+        dims = init_meta.get_dim_bounds()
+        dims = {variable: dims[variable] for variable in dims if variable != init_meta.get_time_var_name()}
         dims_hash = str(dims)
 
         if dims_hash not in dim_hashes:
@@ -426,14 +454,16 @@ class HFCollection:
     tree.
     """
 
-    def __init__(self, hf_dir, num_processes=1, meta_map=None, hf_groups=None, step_map=None, hf_glob_pattern="*.nc*", dask_client=None, multistep_slice_map={}):
+    def __init__(self, hf_dir, num_processes=1, meta_map=None, hf_groups=None, step_map=None, hf_glob_pattern="*.nc*", dask_client=None, multistep_slice_map=None):
         """
         Discovers history files under ``hf_dir``, without reading their metadata.
 
         The pre-computed arguments below are how :meth:`copy` hands state to a
         derived collection; callers normally pass only the first few.
 
-        :param hf_dir: Root directory to search for history files.
+        :param hf_dir: Root directory to search for history files. Made absolute
+            (without resolving symlinks), so relative paths are relative to the
+            current directory at construction.
         :type hf_dir: str
         :param num_processes: Worker processes used for parallel metadata reads.
         :type num_processes: int
@@ -443,38 +473,44 @@ def __init__(self, hf_dir, num_processes=1, meta_map=None, hf_groups=None, step_
         :type hf_groups: dict or None
         :param step_map: Pre-computed ``{path: timedelta}`` mapping.
         :type step_map: dict or None
-        :param hf_glob_pattern: ``fnmatch`` pattern used to discover files.
+        :param hf_glob_pattern: ``fnmatch`` pattern used to discover files. Unused
+            when ``meta_map`` is given: its keys are the collection's paths.
         :type hf_glob_pattern: str
         :param multistep_slice_map: Pre-computed slice indices for multi-timestep
             files (see :meth:`get_multistep_slices`).
-        :type multistep_slice_map: dict
+        :type multistep_slice_map: dict or None
         :param dask_client: Deprecated. Pass ``num_processes`` instead.
-        :raises FileNotFoundError: If no file under ``hf_dir`` matches the pattern.
+        :raises FileNotFoundError: If ``meta_map`` is not given and ``hf_dir`` is not
+            an existing directory, or no file under it matches the pattern.
         """
         if dask_client is not None:
             warnings.warn("Dask is no longer implemented in GenTS. Use the 'num_processes' argument to enable parallelism or reference the ReadTheDocs for using Dask..", DeprecationWarning, stacklevel=2)
 
-        self.__raw_paths = find_files(hf_dir, hf_glob_pattern)
+        # Checked before abspath, which would turn "" into the current directory.
+        if meta_map is None and not os.path.isdir(hf_dir):
+            raise FileNotFoundError(f"History file directory '{hf_dir}' does not exist.")
+        # Absolute (symlinks kept), so every path and group key starts with it and
+        # output templates can strip it as a prefix.
+        hf_dir = os.path.abspath(hf_dir)
         self.__num_processes = num_processes
 
-        if len(self.__raw_paths) == 0:
-            raise FileNotFoundError(f"No files matching '{hf_glob_pattern}' found in '{hf_dir}'")
-
-        self.__hf_to_meta_map = {}
-        self.__hf_multistep_slices = multistep_slice_map
         if meta_map is None:
-            for path in self.__raw_paths:
-                self.__hf_to_meta_map[path] = None
-        else:
-            self.__hf_to_meta_map = meta_map
-        
+            meta_map = {path: None for path in find_files(hf_dir, hf_glob_pattern)}
+            if len(meta_map) == 0:
+                raise FileNotFoundError(f"No files matching '{hf_glob_pattern}' found in '{hf_dir}'")
+            if hf_groups is None:
+                logger.info(f"Initialized HFCollection at '{hf_dir}'")
+                logger.info(f"{len(meta_map)} netCDF files found.")
+
+        self.__hf_to_meta_map = meta_map
+        self.__pulled = False
+        self.__hf_multistep_slices = multistep_slice_map
+        if self.__hf_multistep_slices is None:
+            self.__hf_multistep_slices = {}
+
         self.__hf_groups = hf_groups
         self.__hf_dir = Path(hf_dir)
 
-        if meta_map is None and hf_groups is None:
-            logger.info(f"Initialized HFCollection at '{hf_dir}'")
-            logger.info(f"{len(self.__raw_paths)} netCDF files found.")
-
         self.__hf_to_timestep_delta_map = step_map
 
     def __getitem__(self, key):
@@ -502,12 +538,14 @@ def is_pulled(self):
         """
         Returns whether metadata has been loaded for every file in the collection.
 
+        Cached once true: entries only ever go from ``None`` to loaded, so a pulled
+        collection stays pulled. The per-file getters call this once per file.
+
         :rtype: bool
         """
-        for path in self.__hf_to_meta_map:
-            if self.__hf_to_meta_map[path] is None:
-                return False
-        return True
+        if not self.__pulled:
+            self.__pulled = all(meta is not None for meta in self.__hf_to_meta_map.values())
+        return self.__pulled
 
     def get_multistep_slices(self, hf_path):
         """
@@ -587,6 +625,9 @@ def copy(self, num_processes=None, meta_map=None, hf_groups=None, step_map=None,
             meta_map = self.__hf_to_meta_map
         if hf_groups is None and self.is_pulled():
             hf_groups = self.get_groups()
+            if meta_map is not self.__hf_to_meta_map:
+                hf_groups = {group: [path for path in paths if path in meta_map] for group, paths in hf_groups.items()}
+                hf_groups = {group: paths for group, paths in hf_groups.items() if paths}
         if step_map is None:
             step_map = self.__hf_to_timestep_delta_map
         if multistep_slice_map is None:
@@ -841,17 +882,14 @@ def include_time(self, start_year, end_year, start_month=1, start_day=1, end_mon
                 filtered_path_map[path] = self.__hf_to_meta_map[path]
 
         logger.debug(f"Filtered from {start_year}-{start_month:02d}-{start_day:02d} to {end_year}-{end_month:02d}-{end_day:02d} applied to following glob patterns: '{glob_patterns}'")
-        hf_groups = None
-        if self.__hf_groups is not None:
-            hf_groups = sort_hf_groups(list(filtered_path_map.keys()))
-
-        return self.copy(meta_map=filtered_path_map, hf_groups=hf_groups)
+        return self.copy(meta_map=filtered_path_map)
 
     def get_groups(self, check_fragmented=True):
         """
         Returns the collection's ``{group ID: [paths]}`` mapping.
 
-        Groups are built by :func:`sort_hf_groups` on the first call and cached.
+        Groups are built by :func:`sort_hf_groups` on the first call and cached,
+        with each group's paths in time order (by decoded first time value).
 
         :param check_fragmented: Also merge spatially tiled groups via
             :func:`merge_fragmented_groups`, which requires metadata.
@@ -864,9 +902,17 @@ def get_groups(self, check_fragmented=True):
             if check_fragmented:
                 self.check_pulled()
                 self.__hf_groups = merge_fragmented_groups(self.__hf_groups, self.__hf_to_meta_map)
+                self.__hf_groups = {
+                    group: sorted(paths, key=self.__first_time)
+                    for group, paths in self.__hf_groups.items()
+                }
 
         return self.__hf_groups
 
+    def __first_time(self, hf_path):
+        meta = self.__hf_to_meta_map[hf_path]
+        return meta.decode_time_values(np.ma.getdata(np.atleast_1d(meta.get_float_times()))[0])
+
     def slice_groups(self, slice_size_years=10, start_year=0, pattern="*", time_alignment_method="midpoint"):
         """
         Returns a new collection with its groups partitioned into year windows.
@@ -880,7 +926,9 @@ def slice_groups(self, slice_size_years=10, start_year=0, pattern="*", time_alig
         :param slice_size_years: Maximum width of each window in years.
         :type slice_size_years: int
         :param start_year: Year to align windows to; ``None`` uses the collection's
-            own earliest year.
+            own earliest year. Data before it is still covered, by windows on the
+            same alignment (e.g. ``start_year=1851`` with 2-year slices and data
+            from 1850 gives ``1849-1850``, ``1851-1852``).
         :type start_year: int or None
         :param pattern: One or more ``fnmatch`` globs restricting which groups are
             sliced; a group matching none of them passes through unsliced. A
@@ -920,33 +968,19 @@ def slice_groups(self, slice_size_years=10, start_year=0, pattern="*", time_alig
 
             group_meta_map = {path: self.__hf_to_meta_map[path] for path in hf_paths}
             
-            min_year, max_year = get_year_bounds(group_meta_map)
+            min_year, max_year = get_year_bounds(group_meta_map, time_alignment_method)
             if start_year is not None:
-                min_year = start_year
+                if start_year > min_year:
+                    min_year -= (min_year - start_year) % slice_size_years
+                else:
+                    min_year = start_year
             
             time_slices = calculate_year_slices(slice_size_years, min_year, max_year)
 
             hf_slices = {}
             boundary_cache = {}
             for hf_path in hf_paths:
-                meta_ds = self.__hf_to_meta_map[hf_path]
-                float_bnds = meta_ds.get_float_time_bounds()
-                if float_bnds is None or time_alignment_method == "direct_time":
-                    times = np.ma.getdata(np.atleast_1d(meta_ds.get_float_times()))
-                    ref_units = meta_ds.get_time_units()
-                    ref_calendar = meta_ds.get_time_calendar()
-                else:
-                    bnds = np.ma.getdata(float_bnds)
-                    if time_alignment_method == "midpoint":
-                        times = bnds[:, 0] + (bnds[:, 1] - bnds[:, 0]) / 2
-                    elif time_alignment_method == "start_bound":
-                        times = bnds[:, 0]
-                    elif time_alignment_method == "end_bound":
-                        times = bnds[:, 1]
-                    else:
-                        raise ValueError(f"'{time_alignment_method}' is an invalid time-alignment method. Valid methods are ['direct_time', 'midpoint', 'start_bound', 'end_bound']")
-                    ref_units = meta_ds.get_time_bounds_units()
-                    ref_calendar = meta_ds.get_time_bounds_calendar()
+                times, ref_units, ref_calendar = get_aligned_times(self.__hf_to_meta_map[hf_path], time_alignment_method)
 
                 for time_slice in time_slices:
                     cache_key = (time_slice, ref_units, ref_calendar)
@@ -968,7 +1002,7 @@ def slice_groups(self, slice_size_years=10, start_year=0, pattern="*", time_alig
 
                     if len(times) > 1:
                         past_window = times >= upper_num
-                        end_index = int(np.argmax(past_window)) if past_window.any() else len(times) - 1
+                        end_index = int(np.argmax(past_window)) if past_window.any() else len(times)
                         if start_index != 0 or times[-1] >= upper_num:
                             if hf_path in self.__hf_multistep_slices:
                                 assert f"{time_slice[0]}-{time_slice[1]}" not in self.__hf_multistep_slices[hf_path]
diff --git a/gents/meta.py b/gents/meta.py
index 7057a1b..2288dfc 100644
--- a/gents/meta.py
+++ b/gents/meta.py
@@ -498,7 +498,8 @@ def get_meta_from_path(path: str):
     Opens a netCDF file and returns a :class:`netCDFMeta` built from it.
 
     Picklable factory, so metadata can be read inside ``ProcessPoolExecutor``
-    workers.
+    workers. Dates are left undecoded and per-variable attributes unread; the
+    pipeline needs neither from collection metadata.
 
     :param path: Path to the netCDF history file.
     :type path: str
@@ -509,7 +510,7 @@ def get_meta_from_path(path: str):
     ds_meta = None
     try:
         with GenTSDataStore(path, 'r') as ds:
-            ds_meta = netCDFMeta(ds, path, decode_dates=False)
+            ds_meta = netCDFMeta(ds, path, decode_dates=False, load_variable_attrs=False)
     except Exception as e:
         raise type(e)(f"{e} Path: {path}") from e
 
diff --git a/gents/mhfdataset.py b/gents/mhfdataset.py
index 9c8c535..13d8a81 100644
--- a/gents/mhfdataset.py
+++ b/gents/mhfdataset.py
@@ -54,7 +54,7 @@ class MHFDataset:
     variable.
     """
 
-    def __init__(self, hf_paths, preload_var_list=None, memory_limit_bytes=np.inf, load_secondaries=True, preload_primaries=True):
+    def __init__(self, hf_paths, preload_var_list=None, memory_limit_bytes=np.inf, preload_primaries=True):
         """
         Stores the group's paths; opens nothing until :meth:`open` is called.
 
@@ -66,8 +66,6 @@ def __init__(self, hf_paths, preload_var_list=None, memory_limit_bytes=np.inf, l
         :param memory_limit_bytes: Ceiling on cached variable data. Unbounded by
             default.
         :type memory_limit_bytes: float
-        :param load_secondaries: Cache secondary variable data during :meth:`open`.
-        :type load_secondaries: bool
         :param preload_primaries: Cache primary variable data during :meth:`open`.
             ``False`` skips all primary reads up front (``preload_var_list`` is
             ignored); primaries are still read on demand. Used by ``no_data``
@@ -85,16 +83,25 @@ def __init__(self, hf_paths, preload_var_list=None, memory_limit_bytes=np.inf, l
         self.__last_var_read = None
         self.__memory_limit_bytes = memory_limit_bytes
         self.__past_vars_read = []
-        self.__load_secondaries = load_secondaries
         self.__preload_primaries = preload_primaries
         self.__preload_var_list = list(preload_var_list) if preload_var_list is not None else []
         self.__sorted_time_vals = None
         self.__fragmented = None
 
+    def __is_read_once(self, var_name, hf_meta):
+        """
+        Whether a variable is cached once, from the first file: it has no time
+        dimension and isn't a coordinate variable. Coordinate variables stay per
+        file because fragmented groups place each tile by its own coordinates.
+        """
+        dims = hf_meta.get_variable_dims(var_name)
+        return hf_meta.get_time_var_name() not in dims and var_name not in dims
+
     def __estimate_var_bytes(self, var_name, hf_meta):
         """
         Estimates a variable's aggregated size across the group by scaling one
-        file's shape by the file count.
+        file's shape by the file count, or one file's size for a variable cached
+        once (see :meth:`__is_read_once`).
 
         Exact when every file holds the same number of time steps; approximate
         otherwise. Used only to plan the preload cache before the group has been
@@ -102,6 +109,8 @@ def __estimate_var_bytes(self, var_name, hf_meta):
         """
         shape = hf_meta.get_variable_shapes(var_name)
         itemsize = hf_meta.get_variable_dtype(var_name).itemsize
+        if self.__is_read_once(var_name, hf_meta):
+            return int(np.prod(shape)) * itemsize
         return int(np.prod(shape)) * itemsize * len(self.__hf_files)
 
     def __plan_cacheable_vars(self, candidate_names, hf_meta, running_total=0):
@@ -131,10 +140,12 @@ def open(self):
         ``sub_time_index`` is that value's position within its own file's time
         array -- precomputed so :meth:`get_var_vals` never rescans a time array.
 
-        Secondary variables are cached per file rather than read once from the
-        first file, since time, bounds and per-tile coordinates differ between
-        files; they are small enough for that to be cheap. Preloaded primaries
-        ride along in the same pass so no file has to be reopened for them.
+        Time-varying secondaries (time, bounds) and coordinate variables are
+        cached per file: fragmented groups place each tile by its own
+        coordinates. Other variables without a time dimension (grids, areas)
+        are cached once, from the first file, the only copy :meth:`get_var_vals`
+        reads. Preloaded primaries ride along in the same pass so no file has to
+        be reopened for them.
 
         :raises Exception: If the spatial fragmentation is not consistent over time.
         """
@@ -153,7 +164,7 @@ def open(self):
                     # group does not actually have.
                     checked_preload_list = []
                     for var_name in self.__preload_var_list:
-                        if self.__load_secondaries and var_name in hf_meta.get_secondary_variables():
+                        if var_name in hf_meta.get_secondary_variables():
                             continue
                         if var_name not in hf_meta.get_variables():
                             raise KeyError(f"Attempted to cache non-existent variable '{var_name}'.")
@@ -169,8 +180,7 @@ def open(self):
                     else:
                         self.__preload_var_list = []
 
-                    secondary_names = list(hf_meta.get_secondary_variables()) if self.__load_secondaries else []
-                    vars_to_cache_secondary = self.__plan_cacheable_vars(secondary_names, hf_meta)
+                    vars_to_cache_secondary = self.__plan_cacheable_vars(list(hf_meta.get_secondary_variables()), hf_meta)
                     secondary_bytes = sum(self.__estimate_var_bytes(name, hf_meta) for name in vars_to_cache_secondary)
                     vars_to_cache_primary = self.__plan_cacheable_vars(self.__preload_var_list, hf_meta, secondary_bytes)
 
@@ -191,18 +201,21 @@ def open(self):
                     else:
                         self.__time_mapping[time] = [(hf_index, sub_t_index)]
 
-                if self.__load_secondaries:
-                    for var_name in vars_to_cache_secondary:
-                        if var_name not in hf_meta.get_secondary_variables():
-                            continue
-                        if var_name in self.__data_secondary_var_cache:
-                            self.__data_secondary_var_cache[var_name].append(hf_ds[var_name][:])
-                        else:
-                            self.__data_secondary_var_cache[var_name] = [hf_ds[var_name][:]]
+                for var_name in vars_to_cache_secondary:
+                    if var_name not in hf_meta.get_secondary_variables():
+                        continue
+                    if hf_index > 0 and self.__is_read_once(var_name, hf_meta):
+                        continue
+                    if var_name in self.__data_secondary_var_cache:
+                        self.__data_secondary_var_cache[var_name].append(hf_ds[var_name][:])
+                    else:
+                        self.__data_secondary_var_cache[var_name] = [hf_ds[var_name][:]]
 
                 for var_name in vars_to_cache_primary:
                     if var_name not in hf_meta.get_primary_variables():
                         continue
+                    if hf_index > 0 and self.__is_read_once(var_name, hf_meta):
+                        continue
                     if var_name in self.__data_var_cache:
                         self.__data_var_cache[var_name].append(hf_ds[var_name][:])
                     else:
@@ -263,23 +276,22 @@ def __cache_variable(self, target_var_name, cache_ahead=True):
                     cache_size_b += self.get_var_dsize(nvar_name)
                     vars_to_cache.append(nvar_name)
 
-        # A variable without a time dimension is identical in every file, so it
-        # is read once from the first rather than once per file.
-        notime_vars_to_cache = []
-        time_vars_to_cache = []
+        once_vars_to_cache = []
+        per_file_vars_to_cache = []
         for var_name in vars_to_cache:
-            if self.__time_name in self.get_var_dimensions(var_name):
-                time_vars_to_cache.append(var_name)
+            if self.__is_read_once(var_name, self.__hf_metas[0]):
+                once_vars_to_cache.append(var_name)
             else:
-                notime_vars_to_cache.append(var_name)
+                per_file_vars_to_cache.append(var_name)
 
-        for index, path in enumerate(self.__hf_files):
+        hf_files = self.__hf_files if per_file_vars_to_cache else self.__hf_files[:1]
+        for index, path in enumerate(hf_files):
             with GenTSDataStore(path, 'r') as hf_ds:
+                hf_ds.set_auto_maskandscale(False)
                 if index == 0:
-                    for var_name in notime_vars_to_cache:
-                        self.__data_var_cache[var_name] = hf_ds[var_name][:]
-                
-                for var_name in time_vars_to_cache:
+                    for var_name in once_vars_to_cache:
+                        self.__data_var_cache[var_name] = [hf_ds[var_name][:]]
+                for var_name in per_file_vars_to_cache:
                     if var_name in self.__data_var_cache:
                         self.__data_var_cache[var_name].append(hf_ds[var_name][:])
                     else:
@@ -320,6 +332,7 @@ def __get_hf_data(self, index, var_name, time_slice=None):
         else:
             # Too big to cache whole; reread just the requested region.
             with GenTSDataStore(self.__hf_files[index], 'r') as hf_ds:
+                hf_ds.set_auto_maskandscale(False)
                 if time_slice is None:
                     return hf_ds[var_name][:]
                 return hf_ds[var_name][time_slice]
@@ -334,6 +347,15 @@ def get_var_dsize(self, var_name):
         """
         return np.prod(self.get_var_data_shape(var_name))*self.get_var_dtype(var_name).itemsize
 
+    def get_time_var_name(self):
+        """
+        Returns the name of the group's time variable, as spelled in its files
+        (``time``, ``Time``, ...).
+
+        :rtype: str
+        """
+        return self.__time_name
+
     def get_time_vals(self):
         """
         Returns the sorted, unique float time values across the group.
@@ -449,7 +471,7 @@ def get_var_vals(self, var_name, time_index_start=0, time_index_end=None):
         if var_name in self.__data_coords:
             return self.__data_coords[var_name]
 
-        if "time" not in self.get_var_dimensions(var_name):
+        if self.__time_name not in self.get_var_dimensions(var_name):
             return self.__get_hf_data(0, var_name)
 
         time_vals = self.get_time_vals()[time_index_start:time_index_end]
diff --git a/gents/tests/test_cases.py b/gents/tests/test_cases.py
index 88b141b..a3b1b63 100644
--- a/gents/tests/test_cases.py
+++ b/gents/tests/test_cases.py
@@ -492,6 +492,28 @@ def extraneous_file_case(tmp_path_factory):
     return head_hf_dir, head_ts_dir
 
 
+@pytest.fixture(scope="function")
+def straddling_case(tmp_path_factory):
+    """
+    A callable writing 6-step monthly history files that start in April, so each
+    one straddles a year boundary. File ``k`` holds months ``3+6k .. 8+6k``
+    (month 0 = January of ``CASE_START_YEAR``), and ``VAR0`` holds each step's
+    month index so tests can see which months were written.
+    """
+    head_hf_dir = tmp_path_factory.mktemp("straddling_history_files")
+    head_ts_dir = tmp_path_factory.mktemp("straddling_timeseries_files")
+
+    def write(file_indices):
+        for k in file_indices:
+            months = np.arange(3 + 6*k, 9 + 6*k)
+            path = f"{head_hf_dir}/testing.hf.{str(k).zfill(5)}.nc"
+            generate_history_file(path, (months + 0.5)*30, [[m*30, (m+1)*30] for m in months], num_vars=1)
+            with GenTSDataStore(path, "a") as ds:
+                ds["VAR0"][:] = months[:, None, None] * np.ones((len(months), 3, 4))
+
+    return head_hf_dir, head_ts_dir, write
+
+
 @pytest.fixture(scope="function")
 def continued_case(tmp_path_factory):
     """
diff --git a/gents/tests/test_cli.py b/gents/tests/test_cli.py
index e926783..8603881 100644
--- a/gents/tests/test_cli.py
+++ b/gents/tests/test_cli.py
@@ -290,3 +290,63 @@ def test_cli_append_keeps_config_filters(structured_case):
     assert len(ts_paths) > 0
     for path in ts_paths:
         assert "/0_dir/" not in str(path)
+
+
+def test_cli_overwrite_rewrites_complete_output(simple_case):
+    """-w/--overwrite rewrites existing complete output; without it, it is kept."""
+    input_head_dir, output_head_dir = simple_case
+    argv = ["run_gents", str(input_head_dir), "-o", str(output_head_dir)]
+    with patch.object(sys, "argv", argv):
+        main()
+    ts_path = find_files(output_head_dir, "*.nc")[0]
+
+    for extra_args, marker_survives in (([], True), (["-w"], False)):
+        with GenTSDataStore(ts_path, "a") as ts_ds:
+            ts_ds.setncattr("marker", "from the previous run")
+        with patch.object(sys, "argv", argv + extra_args):
+            main()
+        with GenTSDataStore(ts_path, "r") as ts_ds:
+            assert ("marker" in ts_ds.ncattrs()) == marker_survives
+
+
+def test_cli_compression_overrides_model_config(tmp_path):
+    """--compression/--level win over the model YAML's compression (CESM3: zlib 2)."""
+    hist_dir = tmp_path / "case" / "atm" / "hist"
+    hist_dir.mkdir(parents=True)
+    for index in range(3):
+        generate_history_file(f"{hist_dir}/case.cam.h0.{index:05d}.nc", [(index+0.5)*30], [[index*30, (index+1)*30]], num_vars=1)
+    output_dir = tmp_path / "out"
+
+    with patch.object(sys, "argv", [
+        "run_gents", str(tmp_path / "case"), "-o", str(output_dir), "--model", "CESM3",
+        "--compression", "zlib", "--level", "5"
+    ]):
+        main()
+
+    ts_paths = find_files(output_dir, "*.nc")
+    assert len(ts_paths) == 1
+    with GenTSDataStore(ts_paths[0], "r") as ts_ds:
+        assert ts_ds["VAR0"].filters()["complevel"] == 5
+
+
+def test_parse_align_method_rejects_unknown_methods():
+    """--align_method only accepts the four methods slice_groups/update_ts_orders know."""
+    with patch.object(sys, "argv", ["run_gents", "hf_dir", "--align_method", "bogus_method"]):
+        with pytest.raises(SystemExit):
+            parse_arguments()
+
+
+def test_cli_align_method_used_for_slicing_and_names(tmp_path):
+    """--align_method reaches both slicing and output naming. With end_bound, each
+    December (ending on January 1st) moves into the next year's slice."""
+    for index in range(24):  # monthly, 1850-1851
+        generate_history_file(f"{tmp_path}/testing.hf.{index:05d}.nc", [(index+0.5)*30], [[index*30, (index+1)*30]], num_vars=1)
+    output_dir = tmp_path / "out"
+
+    with patch.object(sys, "argv", [
+        "run_gents", str(tmp_path), "-o", str(output_dir), "--slice", "1", "--align_method", "end_bound"
+    ]):
+        main()
+
+    ranges = sorted(path.name.split(".")[-2] for path in find_files(output_dir, "*.nc"))
+    assert ranges == ["185002-185012", "185101-185112", "185201-185201"]
diff --git a/gents/tests/test_hfcollection.py b/gents/tests/test_hfcollection.py
index bc14ce6..dd57914 100644
--- a/gents/tests/test_hfcollection.py
+++ b/gents/tests/test_hfcollection.py
@@ -1,3 +1,4 @@
+from unittest.mock import patch
 from gents.tests.test_cases import *
 from gents.hfcollection import *
 from gents.meta import netCDFMeta
@@ -27,6 +28,9 @@ def test_calculate_year_slices():
     assert calculate_year_slices(10, 1, 30) == [(1, 10), (11, 20), (21, 30)]
     assert calculate_year_slices(1, 0, 3) == [(0, 0), (1, 1), (2, 2), (3, 3)]
     assert calculate_year_slices(5, 3, 12) == [(3, 7), (8, 12)]
+    assert calculate_year_slices(10, 1, 10) == [(1, 10)]
+    assert calculate_year_slices(10, 1, 11) == [(1, 10), (11, 20)]
+    assert calculate_year_slices(1, 1850, 1851) == [(1850, 1850), (1851, 1851)]
 
 
 def test_hf_sorting(structured_case):
@@ -185,7 +189,7 @@ def test_simple_hfcollection(simple_case, caplog):
         repeat_years = []
         year_slices = calculate_year_slices(slice_size, min_year, max_year)
         for lower, upper in year_slices:
-            assert upper - lower <= slice_size
+            assert upper - lower < slice_size
             assert lower <= upper
             assert lower not in repeat_years
             assert upper not in repeat_years
@@ -500,10 +504,15 @@ def test_get_timestep_delta_matches_reference(simple_case, multistep_large_case,
         assert hf_collection.get_timestep_delta(path) == timedelta(0)
 
 
-def test_no_history_files():
-    """No history files found should raise an error."""
-    with pytest.raises(FileNotFoundError) as exc:
-        empty_hfcollection = HFCollection("")
+def test_no_history_files(tmp_path, monkeypatch):
+    """An empty or missing input path raises, rather than falling back to the current
+    directory (run from one holding .nc files, as the Docker test stages are)."""
+    generate_history_file(f"{tmp_path}/stray.nc", [15], [[0, 30]])
+    monkeypatch.chdir(tmp_path)
+
+    for missing in ("", "no_such_dir"):
+        with pytest.raises(FileNotFoundError):
+            HFCollection(missing)
 
 
 def test_extraneous_hfcollection(extraneous_file_case):
@@ -758,3 +767,103 @@ def test_include_time_month_granularity(long_case):
     # Only the second half.
     filtered = hf_collection.include_time(CASE_START_YEAR, CASE_START_YEAR + 1, start_month=7, end_month=1)
     assert len(filtered) == 6
+
+
+def test_collections_do_not_share_multistep_slices(straddling_case):
+    """Each new HFCollection starts with its own multistep slice map, so slicing the
+    same files again in one process (e.g. a continued run in a notebook) works."""
+    input_head_dir, output_head_dir, write = straddling_case
+    write(range(5))
+
+    for _ in range(2):
+        hf_collection = HFCollection(input_head_dir)
+        hf_collection.pull_metadata()
+        assert hf_collection.get_multistep_slices(next(iter(hf_collection))) is None
+        hf_collection.slice_groups(slice_size_years=1, start_year=CASE_START_YEAR)
+
+
+def test_filters_after_pull_drop_filtered_paths_from_groups(structured_case):
+    """include()/exclude() after pull_metadata() keep only surviving paths in the
+    inherited groups, so slice_groups() doesn't look up a dropped file."""
+    input_head_dir, output_head_dir = structured_case
+    hf_collection = HFCollection(input_head_dir)
+    hf_collection.pull_metadata()
+
+    for filtered in (hf_collection.include("*/0_dir/*"), hf_collection.exclude("*/0_dir/*")):
+        grouped_paths = [path for paths in filtered.get_groups().values() for path in paths]
+        assert sorted(grouped_paths) == sorted(filtered)
+        filtered.slice_groups()
+
+
+def test_include_years_keeps_merged_and_sliced_groups(spatial_fragment_case):
+    """A time filter that drops nothing leaves fragment merging and slicing intact."""
+    input_head_dir, output_head_dir = spatial_fragment_case
+    hf_collection = HFCollection(input_head_dir)
+    hf_collection.pull_metadata()
+    sliced = hf_collection.slice_groups(slice_size_years=1)
+
+    assert hf_collection.include_years(0, 99999).get_groups() == hf_collection.get_groups()
+    assert sliced.include_years(0, 99999).get_groups() == sliced.get_groups()
+
+
+def test_slice_start_year_after_data_extends_windows_backwards(tmp_path):
+    """start_year aligns windows; data before it gets windows on the same alignment
+    instead of being dropped."""
+    for index in range(36):  # monthly, 1850-1852
+        generate_history_file(f"{tmp_path}/testing.hf.{index:05d}.nc", [(index+0.5)*30], [[index*30, (index+1)*30]], num_vars=1)
+    hf_collection = HFCollection(tmp_path)
+    hf_collection.pull_metadata()
+
+    groups = hf_collection.slice_groups(slice_size_years=2, start_year=CASE_START_YEAR + 1).get_groups()
+
+    assert [group.split("[sorting_pivot]")[1] for group in groups] == ["1849-1850", "1851-1852"]
+    assert sum(len(paths) for paths in groups.values()) == 36
+
+
+@pytest.mark.parametrize("time_name", ["time", "Time"])
+def test_fragmented_tiles_merge_whatever_the_time_name(tmp_path, time_name):
+    """Tiled files merge into one group across timesteps whether the time
+    coordinate is spelled 'time' or 'Time' (it differs per timestep, so it must be
+    left out of the tiles' coordinate key)."""
+    dim_shapes = {time_name: None, "bnds": 2, "lat": 1, "lon": 2}
+    for step in range(4):
+        for tile, lat in enumerate([-45.0, 45.0]):
+            generate_history_file(f"{tmp_path}/testing.hf.{step:05d}.nc.{tile}", [(step+0.5)*30], [[step*30, (step+1)*30]], num_vars=1,
+                                  time_name=time_name, time_bounds_name=f"{time_name}_bounds", dim_shapes=dim_shapes,
+                                  dim_vals={"lat": [lat], "lon": [0.0, 90.0]})
+    hf_collection = HFCollection(tmp_path)
+    hf_collection.pull_metadata()
+
+    assert [len(paths) for paths in hf_collection.get_groups().values()] == [8]
+
+
+def test_copies_reuse_paths_instead_of_rewalking(tmp_path):
+    """Filters build copies from the parent's paths rather than re-walking the input
+    tree, which also used the default glob and so broke collections built with
+    another (and walked a 41k-file tree 10 times per CLI run)."""
+    for index in range(3):
+        generate_history_file(f"{tmp_path}/testing.hf.{index:05d}.hist", [(index+0.5)*30], [[index*30, (index+1)*30]], num_vars=1)
+
+    with patch("gents.hfcollection.find_files", wraps=find_files) as walks:
+        hf_collection = HFCollection(tmp_path, hf_glob_pattern="*.hist")
+        filtered = hf_collection.include("*").exclude("*00002*")
+        emptied = hf_collection.include([])
+
+    assert walks.call_count == 1  # construction only; the filters reuse its paths
+    assert len(filtered) == 2
+    assert len(emptied) == 0
+
+
+@pytest.mark.parametrize("method", ["midpoint", "start_bound", "end_bound", "direct_time"])
+def test_slice_groups_keeps_every_file_for_every_alignment(tmp_path, method):
+    """Slice windows span the years of the same aligned times the files are assigned
+    by. With end_bound, the last December ends on January 1st of the next year and
+    used to fall past the last (midpoint-based) window."""
+    for index in range(24):  # monthly, 1850-1851
+        generate_history_file(f"{tmp_path}/testing.hf.{index:05d}.nc", [(index+0.5)*30], [[index*30, (index+1)*30]], num_vars=1)
+    hf_collection = HFCollection(tmp_path)
+    hf_collection.pull_metadata()
+
+    groups = hf_collection.slice_groups(slice_size_years=1, start_year=None, time_alignment_method=method).get_groups()
+
+    assert sum(len(paths) for paths in groups.values()) == 24
diff --git a/gents/tests/test_meta.py b/gents/tests/test_meta.py
index 136fc09..8d47cd6 100644
--- a/gents/tests/test_meta.py
+++ b/gents/tests/test_meta.py
@@ -395,6 +395,16 @@ def test_get_meta_from_path_defers_decoding(tmp_path):
     assert len(meta.get_cftimes()) == 2
 
 
+def test_get_meta_from_path_skips_variable_attrs(tmp_path):
+    """get_meta_from_path() leaves per-variable attrs unread: nothing reads them from
+    collection metadata, and holding them for every file dominates the parent's memory."""
+    path = str(tmp_path / "test.nc")
+    generate_history_file(path, [15.0], [[0.0, 30.0]])
+    meta = get_meta_from_path(path)
+    with pytest.raises(RuntimeError, match="load_variable_attrs"):
+        meta.get_variable_attrs("VAR0")
+
+
 def test_netcdfmeta_is_valid_does_not_decode(tmp_path):
     """is_valid() decides from raw values alone, never triggering a cftime decode."""
     from unittest.mock import patch
diff --git a/gents/tests/test_mhfdataset.py b/gents/tests/test_mhfdataset.py
index f333bae..91855d4 100644
--- a/gents/tests/test_mhfdataset.py
+++ b/gents/tests/test_mhfdataset.py
@@ -405,3 +405,73 @@ def test_MHFDataset_memory_limited_reads_correct_slices(tmp_path):
         vals = ds.get_var_vals("VAR0", time_index_start=start, time_index_end=end)
         assert np.array_equal(vals, source[start:end])
     ds.close()
+
+
+def test_MHFDataset_reads_raw_values_on_every_path(tmp_path):
+    """Preloaded, lazily cached and uncacheable direct reads all return the raw on-disk
+    values, as open() does: no scaling (which corrupts packed integers) and no masking
+    (whose uncounted mask array lets the cache overrun memory_limit_bytes)."""
+    raw = np.arange(36, dtype="i2").reshape(3, 3, 4)
+    paths = []
+    for findex in range(3):
+        path = str(tmp_path / f"testing.hf.{findex:05d}.nc")
+        generate_history_file(path, [findex * 30.0], [[findex * 30.0, (findex + 1) * 30.0]], num_vars=1, dtype="i2")
+        with GenTSDataStore(path, "a") as ds:
+            ds["VAR0"].scale_factor = 0.5
+            ds["VAR0"].set_auto_scale(False)
+            ds["VAR0"][:] = raw[findex:findex + 1]
+        paths.append(path)
+
+    preloaded = MHFDataset(paths)
+    lazy = MHFDataset(paths, preload_primaries=False)
+    direct = MHFDataset(paths, memory_limit_bytes=1)
+    for ds in (preloaded, lazy, direct):
+        with ds:
+            assert np.array_equal(ds.get_var_vals("VAR0"), raw)
+            for cached in ds._MHFDataset__data_var_cache.get("VAR0", []):
+                assert not np.ma.isMaskedArray(cached)
+
+
+def test_MHFDataset_uncached_static_variable_reads_whole_array_from_first_file(tmp_path):
+    """A no-time variable too big for the cache is read directly from the first file
+    alone, and whole: not truncated to its first row, and without opening the rest of
+    the group (41 statics x 3,650 daily CICE files was ~150k opens)."""
+    grid = np.arange(20, dtype=float).reshape(4, 5)
+    paths = []
+    for findex in range(3):
+        path = str(tmp_path / f"testing.hf.{findex:05d}.nc")
+        generate_history_file(path, [findex * 30.0], [[findex * 30.0, (findex + 1) * 30.0]], num_vars=1)
+        with GenTSDataStore(path, "a") as ds:
+            ds.createDimension("ny", 4)
+            ds.createDimension("nx", 5)
+            ds.createVariable("grid", float, ("ny", "nx"))[:] = grid
+        paths.append(path)
+
+    # 150 B: under one grid copy (160 B), so it is never cached.
+    with MHFDataset(paths, memory_limit_bytes=150, preload_primaries=False) as ds:
+        assert "grid" not in ds._MHFDataset__data_secondary_var_cache
+        with patch("gents.mhfdataset.GenTSDataStore", wraps=GenTSDataStore) as opens:
+            assert np.array_equal(ds.get_var_vals("grid"), grid)
+        assert opens.call_count == 1
+
+
+def test_MHFDataset_open_caches_static_secondary_once(tmp_path):
+    """open() budgets and caches a no-time secondary as one copy from the first file,
+    so reading it later opens no files."""
+    grid = np.arange(20, dtype=float).reshape(4, 5)
+    paths = []
+    for findex in range(3):
+        path = str(tmp_path / f"testing.hf.{findex:05d}.nc")
+        generate_history_file(path, [findex * 30.0], [[findex * 30.0, (findex + 1) * 30.0]], num_vars=1)
+        with GenTSDataStore(path, "a") as ds:
+            ds.createDimension("ny", 4)
+            ds.createDimension("nx", 5)
+            ds.createVariable("grid", float, ("ny", "nx"))[:] = grid
+        paths.append(path)
+
+    # 300 B fits one grid copy (160 B) plus time and bounds, not one per file (480 B).
+    with MHFDataset(paths, memory_limit_bytes=300, preload_primaries=False) as ds:
+        assert len(ds._MHFDataset__data_secondary_var_cache["grid"]) == 1
+        with patch("gents.mhfdataset.GenTSDataStore", wraps=GenTSDataStore) as opens:
+            assert np.array_equal(ds.get_var_vals("grid"), grid)
+        assert opens.call_count == 0
diff --git a/gents/tests/test_parallel.py b/gents/tests/test_parallel.py
index ea61e49..0bb0ba9 100644
--- a/gents/tests/test_parallel.py
+++ b/gents/tests/test_parallel.py
@@ -121,19 +121,40 @@ def test_execute_memory_limit_bytes_threaded_to_MHFDataset(simple_case):
             assert call.kwargs["memory_limit_bytes"] == DEFAULT_MEMORY_LIMIT_BYTES
 
 
-def test_execute_memory_limit_bytes_not_yet_enforced_by_preload(simple_case):
-    """
-    Known current limitation, pinned down so a future change doesn't silently
-    alter it either way: memory_limit_bytes reaches MHFDataset correctly (see
-    test_execute_memory_limit_bytes_threaded_to_MHFDataset), but MHFDataset's
-    preload_var_list path -- what every order taken through TSCollection.execute()
-    uses -- does not itself check the limit while preloading. The limit is only
-    enforced by the on-demand __cache_variable fallback, which preload bypasses.
-    An absurdly small limit should therefore NOT prevent a normal run from
-    succeeding today. If this test starts failing because execute() now raises
-    or drops data under a tiny limit, that's preload becoming memory-aware --
-    update this test to assert the new, real enforcement instead of removing it.
-    """
+def test_execute_peak_memory_stays_near_memory_limit(tmp_path):
+    """Regression test for a user-reported OOM on a missing-value clone: the peak of
+    everything allocated while generating time series stays within memory_limit_bytes
+    plus a few write chunks. Lazily cached variables used to carry an uncounted
+    mask array (+25% for float32, full for all-fill data) that pushed the cache past
+    the limit."""
+    import tracemalloc
+    from gents.timeseries import CHUNK_TARGET_BYTES
+
+    hf_dir = tmp_path / "hf"
+    makedirs(hf_dir)
+    dim_shapes = {"time": None, "bnds": 2, "lat": 1000, "lon": 1000, "lev": 1}
+    for findex in range(24):
+        path = str(hf_dir / f"case.cam.h0.{findex:04d}.nc")
+        generate_history_file(path, [findex * 30.0 + 15], [[findex * 30.0, (findex + 1) * 30.0]],
+                              num_vars=0, dim_shapes=dim_shapes)
+        with GenTSDataStore(path, "a") as ds:
+            # Never written, so every value reads back as the NaN fill, as in a clone.
+            for name in ("VAR0", "VAR1"):
+                ds.createVariable(name, "f4", ("time", "lat", "lon"), fill_value=np.nan)
+
+    # Room for one variable (24 x 4 MB = 91.6 MiB) but not both, so VAR1 is cached lazily.
+    limit = 96 * 1024**2
+    ts_collection = TSCollection(HFCollection(hf_dir), tmp_path / "ts", num_processes=1)
+    tracemalloc.start()
+    ts_collection.execute(memory_limit_bytes=limit, show_progress=False)
+    peak = tracemalloc.get_traced_memory()[1]
+    tracemalloc.stop()
+    assert peak < limit + 3 * CHUNK_TARGET_BYTES, f"peak {peak / 1024**2:.1f} MiB"
+
+
+def test_execute_tiny_memory_limit_still_writes_every_file(simple_case):
+    """A memory_limit_bytes too small to cache anything still produces every file:
+    open() preloads nothing and each read streams from disk instead."""
     input_head_dir, output_head_dir = simple_case
     hf_collection = HFCollection(input_head_dir, num_processes=1)
     ts_collection = TSCollection(hf_collection, output_head_dir, num_processes=1)
diff --git a/gents/tests/test_tscollection.py b/gents/tests/test_tscollection.py
index 5e7a2b4..fec691e 100644
--- a/gents/tests/test_tscollection.py
+++ b/gents/tests/test_tscollection.py
@@ -790,3 +790,145 @@ def test_skip_existing_ts_string_reflects_trimmed_range(continued_case):
 
     for order in resumed:
         assert order["ts_string"] == "185101-185112"
+
+
+def test_skip_existing_recomputes_cut_indices_for_trimmed_order(straddling_case):
+    """Trimming leading files from an order whose files straddle slice boundaries
+    shifts its time indices onto the trimmed group."""
+    input_head_dir, output_head_dir, write = straddling_case
+
+    def build():
+        hf_collection = HFCollection(input_head_dir)
+        hf_collection.pull_metadata()
+        return TSCollection(hf_collection.slice_groups(slice_size_years=1, start_year=CASE_START_YEAR), str(output_head_dir))
+
+    write(range(5))    # Apr 1850 .. Sep 1852
+    build().execute(raise_errors=True, show_progress=False)
+    write(range(5, 7))  # Oct 1852 .. Sep 1853
+
+    resumed = build().skip_existing()
+
+    # 1852's order was files 3-5 (steps 3..15); files 3-4 are covered, leaving
+    # file 5's first three steps (Oct-Dec 1852).
+    trimmed = [order for order in resumed if order["ts_string"].startswith("1852")][0]
+    assert [Path(path).name for path in trimmed["hf_paths"]] == ["testing.hf.00005.nc"]
+    assert (trimmed["ts_start_index"], trimmed["ts_end_index"]) == (0, 3)
+
+    resumed.execute(raise_errors=True, show_progress=False)
+    with GenTSDataStore(f"{output_head_dir}/testing.hf.VAR0.185210-185212.nc", "r") as ds:
+        assert ds["VAR0"][:, 0, 0].tolist() == [33, 34, 35]
+
+
+def test_trailing_single_file_slice_keeps_last_step(straddling_case):
+    """A slice made of one multi-step file that starts before the window and ends
+    inside it keeps every in-window step, including the file's last."""
+    input_head_dir, output_head_dir, write = straddling_case
+    write(range(4))  # Apr 1850 .. Mar 1852; 1852's slice is file 3's Jan-Mar alone
+
+    hf_collection = HFCollection(input_head_dir)
+    hf_collection.pull_metadata()
+    ts_collection = TSCollection(hf_collection.slice_groups(slice_size_years=1, start_year=CASE_START_YEAR), str(output_head_dir))
+    ts_collection.execute(raise_errors=True, show_progress=False)
+
+    with GenTSDataStore(f"{output_head_dir}/testing.hf.VAR0.185201-185203.nc", "r") as ds:
+        assert ds["VAR0"][:, 0, 0].tolist() == [24, 25, 26]
+
+
+def test_relative_input_dir_keeps_output_names(tmp_path, monkeypatch):
+    """A relative input directory whose name repeats in the file prefix (CESM's
+    archive/<case>/atm/hist/<case>.cam.h0.* run as `run_gents <case>`), or ".",
+    gives the same output template as an absolute one."""
+    hist_dir = tmp_path / "case" / "atm" / "hist"
+    hist_dir.mkdir(parents=True)
+    for index in range(2):
+        generate_history_file(f"{hist_dir}/case.cam.h0.{index:05d}.nc", [(index+0.5)*30], [[index*30, (index+1)*30]], num_vars=1)
+    output_dir = str(tmp_path / "out")
+    expected = f"{output_dir}/atm/hist/case.cam.h0"
+
+    for cwd, input_dir in ((tmp_path, "case"), (tmp_path / "case", ".")):
+        monkeypatch.chdir(cwd)
+        assert TSCollection(HFCollection(input_dir), output_dir)[0]["ts_path_template"] == expected
+
+
+def test_auxiliary_order_keeps_every_step_when_first_variable_is_static(tmp_path):
+    """A group with no primaries writes every time step even when its first variable
+    is static (its first dimension used to become the step count)."""
+    input_dir, output_dir = tmp_path / "hist", tmp_path / "out"
+    input_dir.mkdir()
+    for index in range(12):
+        with GenTSDataStore(f"{input_dir}/testing.hf.{index:05d}.nc", "w") as ds:
+            ds.createDimension("ny", 3)
+            ds.createDimension("time", None)
+            ds.createDimension("nbnd", 2)
+            ds.createVariable("area", float, ("ny",))[:] = [1.0, 2.0, 3.0]
+            time = ds.createVariable("time", float, ("time",))
+            time.units, time.calendar = "days since 1850-01-01", "noleap"
+            time[:] = [(index + 0.5) * 30]
+            ds.createVariable("time_bnds", float, ("time", "nbnd"))[:] = [[index * 30, (index + 1) * 30]]
+
+    ts_collection = TSCollection(HFCollection(input_dir), str(output_dir))
+    assert [order["primary_var"] for order in ts_collection] == ["auxiliary"]
+    ts_collection.execute(raise_errors=True, show_progress=False)
+
+    with GenTSDataStore(f"{output_dir}/testing.hf.auxiliary.185001-185012.nc", "r") as ds:
+        assert len(ds["time"]) == 12
+        assert ds["area"][:].tolist() == [1.0, 2.0, 3.0]
+
+
+def test_groups_follow_time_order_when_file_names_do_not(tmp_path):
+    """Cut indices assume a group's first path is its earliest file, so groups are
+    ordered by time even when file names sort the other way."""
+    input_dir, output_dir = tmp_path / "hist", tmp_path / "out"
+    input_dir.mkdir()
+    for k in range(5):  # 6-step files from April 1850, named in reverse time order
+        months = np.arange(3 + 6*k, 9 + 6*k)
+        path = f"{input_dir}/testing.hf.{9 - k:05d}.nc"
+        generate_history_file(path, (months + 0.5)*30, [[m*30, (m+1)*30] for m in months], num_vars=1)
+        with GenTSDataStore(path, "a") as ds:
+            ds["VAR0"][:] = months[:, None, None] * np.ones((len(months), 3, 4))
+
+    hf_collection = HFCollection(input_dir)
+    hf_collection.pull_metadata()
+    for paths in hf_collection.get_groups().values():
+        first_times = [hf_collection[path].get_float_times()[0] for path in paths]
+        assert first_times == sorted(first_times)
+
+    ts_collection = TSCollection(hf_collection.slice_groups(slice_size_years=1, start_year=CASE_START_YEAR), str(output_dir))
+    ts_collection.execute(raise_errors=True, show_progress=False)
+    with GenTSDataStore(f"{output_dir}/testing.hf.VAR0.185101-185112.nc", "r") as ds:
+        assert ds["VAR0"][:, 0, 0].tolist() == list(range(12, 24))
+
+
+@pytest.mark.parametrize("time_name", ["time", "Time"])
+def test_time_dimension_name_is_not_assumed_lowercase(tmp_path, time_name):
+    """A time dimension spelled 'Time' (MOM6 standalone) is written per step like
+    'time', not treated as static with the first file's data repeated."""
+    for index in range(3):
+        path = f"{tmp_path}/testing.hf.{index:05d}.nc"
+        generate_history_file(path, [(index+0.5)*30], [[index*30, (index+1)*30]], num_vars=1,
+                              time_name=time_name, time_bounds_name=f"{time_name}_bounds")
+        with GenTSDataStore(path, "a") as ds:
+            ds["VAR0"][:] = index
+    output_dir = tmp_path / "out"
+
+    TSCollection(HFCollection(tmp_path), str(output_dir)).execute(raise_errors=True, show_progress=False)
+
+    with GenTSDataStore(f"{output_dir}/testing.hf.VAR0.185001-185003.nc", "r") as ds:
+        assert ds.dimensions[time_name].isunlimited()
+        assert ds["VAR0"][:, 0, 0].tolist() == [0, 1, 2]
+        assert ds[f"{time_name}_bounds"][:].tolist() == [[0, 30], [30, 60], [60, 90]]
+
+
+def test_batched_execute_keeps_per_variable_output_paths(tmp_path):
+    """execute(optimize=True) batches variables that share source files, but each
+    still goes to its own ts_path_template (e.g. a var_glob-limited path swap)."""
+    (tmp_path / "hist").mkdir()
+    for index in range(3):
+        generate_history_file(f"{tmp_path}/hist/testing.hf.{index:05d}.nc", [(index+0.5)*30], [[index*30, (index+1)*30]], num_vars=2)
+    output_dir = tmp_path / "out"
+
+    ts_collection = TSCollection(HFCollection(tmp_path), str(output_dir))
+    ts_collection.apply_path_swap("/hist/", "/swapped/", var_glob="VAR1").execute(raise_errors=True, show_progress=False)
+
+    assert isfile(f"{output_dir}/hist/testing.hf.VAR0.185001-185003.nc")
+    assert isfile(f"{output_dir}/swapped/testing.hf.VAR1.185001-185003.nc")
diff --git a/gents/tests/test_utils.py b/gents/tests/test_utils.py
index 4771d2f..e6aa2e0 100644
--- a/gents/tests/test_utils.py
+++ b/gents/tests/test_utils.py
@@ -65,6 +65,18 @@ def test_progress_bar_silent_without_a_terminal(monkeypatch):
     assert bar.count == 2
 
 
+def test_progress_bar_with_nothing_to_do_draws_complete(monkeypatch):
+    """A zero-total bar (an empty execute(), e.g. after skip_existing() finds
+    everything done) draws as complete on a terminal instead of dividing by zero."""
+    stdout = FakeStdout(tty=True)
+    monkeypatch.setattr(sys, "stdout", stdout)
+
+    ProgressBar(total=0, label="Testing")
+
+    assert "100.00%" in stdout.written
+    assert stdout.written.endswith("\n")
+
+
 @pytest.fixture(scope="session")
 def log_output_dir(tmp_path_factory):
     """Session-scoped temp directory for log file output."""
diff --git a/gents/timeseries.py b/gents/timeseries.py
index 332fdd5..9424580 100644
--- a/gents/timeseries.py
+++ b/gents/timeseries.py
@@ -18,6 +18,7 @@
 import traceback
 import logging
 import copy
+import gc
 import warnings
 
 logger = logging.getLogger(__name__)
@@ -186,18 +187,19 @@ def write_timeseries_file(agg_hf_ds, ts_out_path, primary_var, secondary_vars_da
 
     if ts_start_index is None:
         ts_start_index = 0
+    if ts_end_index is None:
+        ts_end_index = len(agg_hf_ds.get_time_vals())
+    time_name = agg_hf_ds.get_time_var_name()
 
     with GenTSDataStore(ts_out_path, mode="w") as ts_ds:
         if primary_var != "auxiliary":
             var_shape = agg_hf_ds.get_var_data_shape(primary_var)
             var_dims = agg_hf_ds.get_var_dimensions(primary_var)
             
-            if ts_end_index is None:
-                ts_end_index = var_shape[0]
             var_shape[0] = ts_end_index - ts_start_index
 
             for index, dim in enumerate(var_dims):
-                if dim == "time":
+                if dim == time_name:
                     ts_ds.createDimension(dim, None)
                 else:
                     ts_ds.createDimension(dim, var_shape[index])
@@ -228,7 +230,7 @@ def write_timeseries_file(agg_hf_ds, ts_out_path, primary_var, secondary_vars_da
 
             if no_data:
                 pass
-            elif len(var_shape) > 0 and "time" in var_dims:
+            elif len(var_shape) > 0 and time_name in var_dims:
                 for i in range(0, var_shape[0], chunksizes[0]):
                     end = min(i + chunksizes[0], var_shape[0])
                     chunk = agg_hf_ds.get_var_vals(
@@ -245,14 +247,12 @@ def write_timeseries_file(agg_hf_ds, ts_out_path, primary_var, secondary_vars_da
             var_shape = agg_hf_ds.get_var_data_shape(secondary_var)
             var_dims = agg_hf_ds.get_var_dimensions(secondary_var)
 
-            if ts_end_index is None:
-                ts_end_index = var_shape[0]
-            if "time" in var_dims:
+            if time_name in var_dims:
                 var_shape[0] = ts_end_index - ts_start_index
 
             for index, dim in enumerate(var_dims):
                 if dim not in ts_ds.dimensions:
-                    if dim == "time":
+                    if dim == time_name:
                         ts_ds.createDimension(dim, None)
                     else:
                         ts_ds.createDimension(dim, var_shape[index])
@@ -275,7 +275,7 @@ def write_timeseries_file(agg_hf_ds, ts_out_path, primary_var, secondary_vars_da
             ts_ds[secondary_var].setncatts(
                 {key: val for key, val in svar_attrs.items() if key != "_FillValue"}
             )
-            if "time" in var_dims:
+            if time_name in var_dims:
                 svar_vals = secondary_vars_data[secondary_var][ts_start_index:ts_end_index]
             else:
                 svar_vals = secondary_vars_data[secondary_var]
@@ -526,7 +526,7 @@ def update_ts_orders(self, strfrmt_kwargs={}, time_alignment_method="midpoint"):
         orders = []
         for index, glob_template in enumerate(self.__groups):
             hf_paths = self.__groups[glob_template]
-            output_template = glob_template.split(str(self.__hf_collection.get_input_dir()))[1]
+            output_template = glob_template.removeprefix(str(self.__hf_collection.get_input_dir()))
             if "[sorting_pivot]" in output_template:
                 output_template, slice_years = output_template.split("[sorting_pivot]")
                 logger.debug(f"Group [{index+1}/{len(self.__groups)}] {len(hf_paths)} files: {output_template}, sliced to [{slice_years}]")
@@ -847,12 +847,18 @@ def skip_existing(self):
             if latest_existing_time is not None:
                 original_hf_paths = order_dict["hf_paths"]
                 new_hf_paths = []
+                dropped_steps = 0
                 for hf_path in original_hf_paths:
                     hf_times = self.__hf_collection[hf_path].get_float_times()
                     hf_times = np.ma.getdata(np.atleast_1d(hf_times))
 
                     if hf_times[-1] > latest_existing_time:
                         new_hf_paths.append(hf_path)
+                    else:
+                        dropped_steps += hf_times.shape[0]
+                if dropped_steps > 0 and order_dict["ts_end_index"] is not None:
+                    order_dict["ts_start_index"] = 0
+                    order_dict["ts_end_index"] -= dropped_steps
 
                 if new_hf_paths and new_hf_paths != original_hf_paths:
                     first_meta = self.__hf_collection[new_hf_paths[0]]
@@ -923,8 +929,9 @@ def execute(self, optimize=True, optimize_batch_n=200, raise_errors=False, no_da
         """
         Runs every order, writing the time series files.
 
-        Orders sharing a first source file and time slice are batched together so
-        each group of history files is opened once rather than once per variable.
+        Orders sharing source files, time slice, output path template and
+        secondary variables are batched together so each group of history files is
+        opened once rather than once per variable.
         Work runs over a process pool when ``num_processes > 1`` and in-process
         otherwise; per-order failures are logged and the rest of the run continues.
 
@@ -954,10 +961,10 @@ def execute(self, optimize=True, optimize_batch_n=200, raise_errors=False, no_da
         if optimize:
             order_index_merge_map = {}
             for index, order in enumerate(self.__orders):
-                first_hf_path = order["hf_paths"][0]
-                start_index = order["ts_start_index"]
-                end_index = order["ts_end_index"]
-                key = f"{first_hf_path}.{start_index}.{end_index}"
+                key = (
+                    tuple(order["hf_paths"]), order["ts_start_index"], order["ts_end_index"],
+                    order["ts_path_template"], tuple(order["secondary_vars"]),
+                )
                 if key in order_index_merge_map:
                     order_index_merge_map[key].append(index)
                 else:
@@ -1009,18 +1016,25 @@ def execute(self, optimize=True, optimize_batch_n=200, raise_errors=False, no_da
                 })
         prog_bar = ProgressBar(total=len(optimized_orders), label="Generating Timeseries", quiet=not show_progress)
         if self.__num_processes > 1:
-            with ProcessPoolExecutor(max_workers=self.__num_processes) as executor:
-                futures = {executor.submit(generate_time_series, **args): args for args in optimized_orders}
-                for future in as_completed(futures):
-                    try:
-                        results.append(future.result())
-                    except Exception as exc:
-                        order = futures[future]
-                        logger.warning(f"Failed to generate time series for {order['ts_path_template']}: {exc}", exc_info=True)
-                        if raise_errors:
-                            raise
-                    finally:
-                        prog_bar.step()
+            # Forked workers share the parent's heap copy-on-write. Freezing the GC
+            # keeps collections from writing to those pages, so workers don't end up
+            # copying the collection's metadata.
+            gc.freeze()
+            try:
+                with ProcessPoolExecutor(max_workers=self.__num_processes) as executor:
+                    futures = {executor.submit(generate_time_series, **args): args for args in optimized_orders}
+                    for future in as_completed(futures):
+                        try:
+                            results.append(future.result())
+                        except Exception as exc:
+                            order = futures[future]
+                            logger.warning(f"Failed to generate time series for {order['ts_path_template']}: {exc}", exc_info=True)
+                            if raise_errors:
+                                raise
+                        finally:
+                            prog_bar.step()
+            finally:
+                gc.unfreeze()
         else:
             for args in optimized_orders:
                 try:
diff --git a/gents/utils.py b/gents/utils.py
index a2a1c7a..d2b01f7 100644
--- a/gents/utils.py
+++ b/gents/utils.py
@@ -119,7 +119,7 @@ def step(self):
         self.count += 1
         if self.quiet:
             return
-        percent = self.count / self.total
+        percent = self.count / self.total if self.total else 1.0
         filled_length = int(self.length * percent)
         bar = "█" * filled_length + "-" * (self.length - filled_length)
         elapsed = time() - self.start_time
diff --git a/pyproject.toml b/pyproject.toml
index 90200b4..ab255db 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -41,7 +41,7 @@ Homepage = "https://github.com/AgentOxygen/GenTS"
 Issues = "https://github.com/AgentOxygen/GenTS"
 
 [tool.setuptools]
-packages = ["gents", "gents.configs"]
+packages = ["gents", "gents.configs", "gents.conformity", "gents.conformity.models"]
 
 [tool.setuptools.package-data]
 "gents.configs" = ["*.yaml"]