diff --git a/CMakeLists.txt b/CMakeLists.txt
index d337db7..29d48c7 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -1166,7 +1166,7 @@ set_tests_properties(cable_statics_modes PROPERTIES LABELS "fortran;driver;hermi
 # clamped hang-off moment under vessel pitch.
 add_executable(test_vessel_motion tests/test_vessel_motion.f90)
 target_link_libraries(test_vessel_motion PRIVATE cabledyn_core)
-file(MAKE_DIRECTORY ${CMAKE_BINARY_DIR}/vessel_motion)
+file(MAKE_DIRECTORY ${CMAKE_BINARY_DIR}/vessel_motion/vm_wk7)
 add_test(NAME vessel_motion COMMAND test_vessel_motion WORKING_DIRECTORY ${CMAKE_BINARY_DIR}/vessel_motion)
 set_tests_properties(vessel_motion PROPERTIES LABELS "fortran;unit;dynamic;validation")
 
diff --git a/doc/driver_format.md b/doc/driver_format.md
index 922859b..b49f153 100644
--- a/doc/driver_format.md
+++ b/doc/driver_format.md
@@ -1329,7 +1329,7 @@ every node's position and time as MoorDyn-C interpolates them. Standalone decks
 | Option | File | Meaning |
 |---|---|---|
 | `WaveKin 3` | `wave_elevation.txt` | `time elevation` rows from `t = 0` (no header rows), resampled linearly at `dtWave` to `floor(t_last/dtWave) + 1` samples (the last dropped when the count is odd) and reduced to its Fourier components; components above `WvHiCOff` (default 0.5 Hz, MoorDyn-C's fixed cut) are dropped and the mean is a mean-level offset, as MoorDyn-C does. Waves travel along +x |
-| `WaveKin 7` | `wave_frequencies.txt` | `omega Re Im [beta]` rows [rad/s, m, m, rad], the first at `omega = 0`: each row is the component `\|c\| sin(omega t - k d + arg c)`, `d = x cos(beta) + y sin(beta)`; the `omega = 0` row is a mean-level offset. All rows must share one `beta` |
+| `WaveKin 7` | `wave_frequencies.txt` | `omega Re Im [beta]` rows [rad/s, m, m, rad], the first at `omega = 0`: each row is the component `\|c\| sin(omega t - k d + arg c)`, `d = x cos(beta) + y sin(beta)` with the row's own `beta` (0 when the column is left out); the `omega = 0` row is a mean-level offset and its `beta` has no effect. Rows of different directions describe a short-crested sea, summed in ascending frequency |
 | `Currents 1` | `current_profile.txt` | `z ux uy uz` rows after its header rows, strictly increasing `z`; linear in `z` and held beyond the end rows |
 | `Currents 2` | `current_profile_dynamic.txt` | a depth profile in time: lines 1-4 and 6 are free text, line 5 lists the depths `z1 ... zn` (strictly increasing, read in single precision as MoorDyn-C does), and each later row is `t ux1 ... uxn [uy1 ... uyn [uz1 ... uzn]]` (missing components zero), at least two rows with strictly increasing `t` |
 | `Currents 5` | `current_profile_4d.txt` | a grid in x, y, z and time: line 1 is `nx ny nz nt`, lines 2-5 the x, y, z and t coordinates (each strictly increasing, at least two times), then one `x y z t ux uy uz` row per grid point, in any order, matched to the grid by value; a text line 6 is a header |
@@ -1368,16 +1368,25 @@ where it is sampled on `water_grid.txt`; CableDyn adds the current at the nodes,
 Rejected by name: `WaveKin 1` (kinematics set node by node through the MoorDyn-C API; a deck has
 no source for them), `WaveKin 2` (the FFT grid: `wave_frequencies.txt` resampled to an even
 spectrum; `WaveKin 7` evaluates the same file's components exactly), `WaveKin 4`-`6` and
-`Currents 3`-`4` (no kinematics source in MoorDyn-C v2.7.1), and a `wave_frequencies.txt` whose
-rows have different directions. A deck that also declares a `waves` option, a WaterKin file, or a
-second current source is rejected as double counting. On shared inputs the `wavekin_modes` test
-holds `WaveKin 7` and `Currents 1` identical to MoorDyn-C v2.7.1 (`5e-12` relative) and `WaveKin 3`
+`Currents 3`-`4` (no kinematics source in MoorDyn-C v2.7.1). A deck that also declares a `waves`
+option, a WaterKin file, or a second current source is rejected as double counting. On shared
+inputs the `wavekin_modes` test holds `WaveKin 7` and `Currents 1` identical to MoorDyn-C v2.7.1 (`5e-12` relative) and `WaveKin 3`
 within `3e-7` m in elevation and `1e-4` in velocity and acceleration at MoorDyn-C's grid points,
 the remainder being MoorDyn-C's grid interpolation and its approximate wave number. The
 `mdc_currents` test holds `Currents 2` and `5` identical to MoorDyn-C v2.7.1 at 739 points and
 times (including past the last sample), and `mdc_currents_runs` matches the tension change of a
 time-varying and of a sheared current to MoorDyn-C's within 0.9 % on a 70 m rope.
 
+`wave_frequencies.txt` is the file MoorDyn-C v2.7.1 reads, a `beta` per row included: its
+`SUM_COMPONENTS` mode gives each row its own direction too. `WaveRefX` and `WaveRefY`
+([options](options.rst)) are an extension of CableDyn and no part of MoorDyn-C or MoorDyn-F: a
+deck with them is not a plain MoorDyn deck, and a MoorDyn solver would not place the phases there.
+
+A non-zero `omega = 0` row raises the still-water level. CableDyn then takes every wave number at
+the depth `WtrDpth` plus that level, as it does for the mean of an elevation record, where
+MoorDyn-C keeps `WtrDpth` and adds the level to the elevation alone. The two differ by the order
+of the level over the depth in the long components.
+
 (waterkin-file-modes)=
 ### MoorDyn-F `WaterKin` file modes
 
diff --git a/doc/file_formats.rst b/doc/file_formats.rst
index eeb19de..65232a1 100644
--- a/doc/file_formats.rst
+++ b/doc/file_formats.rst
@@ -549,7 +549,8 @@ two columns:
    * - 2
      - ``elevation``
      - m
-     - Free-surface elevation η at the reference location.
+     - Free-surface elevation η at the reference location: the origin, or the point of the
+       ``WaveRefX`` and ``WaveRefY`` options, an extension of CableDyn.
 
 * Blank lines are skipped. Lines before the first data row whose first token is not
   a number, such as a ``time elevation`` header, are skipped. Once data has
diff --git a/doc/options.rst b/doc/options.rst
index 1075098..84cd525 100644
--- a/doc/options.rst
+++ b/doc/options.rst
@@ -761,6 +761,16 @@ Ambient fluid and prescribed motion
        component at or below the cut-off is refused
      - standalone
      - ``6.2832 WvHiCOff``
+   * - ``WaveRefX``, ``WaveRefY``
+     - number [m], at most 10\ :sup:`6` in magnitude
+     - ``0``
+     - the horizontal point that the phases of an elevation record or a component file
+       (WaterKin ``WaveKinMod 1``, ``WaveKin 3`` and ``7``) refer to: the record was measured
+       there, and each phase uses the position relative to it (:ref:`theory-waves`). An
+       extension of CableDyn, which neither MoorDyn-C nor MoorDyn-F reads. A deck with any other
+       sea, or none, is refused by name
+     - standalone
+     - ``250.0 WaveRefX``
    * - ``Currents``
      - ``0``, ``1``, ``2`` or ``5``
      - ``0``
diff --git a/doc/theory.rst b/doc/theory.rst
index 2d00636..8b9f9ed 100644
--- a/doc/theory.rst
+++ b/doc/theory.rst
@@ -1112,6 +1112,19 @@ commensurate, the record does not repeat (an equally spaced comb would repeat af
 its component table (amplitude, frequency, phase), for example a WaterKin ``WaveKinMod`` 1
 elevation history, which is reduced once to its Fourier components on the run time base.
 
+The phase of a component of an elevation history or a ``WaveKin 7`` table is
+
+.. math::
+
+   \theta_i = k_i\big[(x - x_r)\cos\beta_i + (y - y_r)\sin\beta_i\big] - \omega_i t + \phi_i,
+
+with :math:`(x_r, y_r)` the ``WaveRefX``, ``WaveRefY`` point where the record was measured, the
+origin by default. The sea of a record referred to a point is therefore the same sea with the
+deck translated, and a vessel RAO takes the vessel's position relative to that point. A history
+has one direction :math:`\beta_i = \beta`. The rows of a ``WaveKin 7`` table may each have
+their own: a short-crested sea, whose elevation and kinematics are the sums of those of its
+directions, stretched by the one total elevation.
+
 Every point sums every component of the elevation. The sums run in interleaved partial sums
 added up in a fixed order, so a result does not depend on the vector width or the thread
 count, and each point is a function of its position and time alone: nothing carries from one
diff --git a/examples/cabledyn_options_reference.dat b/examples/cabledyn_options_reference.dat
index 85a57e0..b9fc312 100644
--- a/examples/cabledyn_options_reference.dat
+++ b/examples/cabledyn_options_reference.dat
@@ -99,6 +99,8 @@ none      waves             - Wave model and parameters (m, s, deg) {none; airy;
 -- 0       Currents         - MoorDyn-C current mode {0=none; 1=current_profile.txt} [default: 0]
 -- 0.25    dtWave           - MoorDyn-C WaveKin 3 resampling step (s) [default: 0.25]
 -- 3.14159 WvHiCOff         - High cut-off of an elevation record, WaveKinMod 1 or WaveKin 3 (rad/s) {value; none} [default: 3.14159]
+-- 0.0     WaveRefX         - x of the point an elevation record or component file refers to (m) [default: 0]
+-- 0.0     WaveRefY         - y of the point an elevation record or component file refers to (m) [default: 0]
 -- Keyword aliases (equivalent spellings, case-insensitive):
 --   g = gravity; rhoW = WtrDnsty = water_density; WtrDpth = water_depth; kBot = kb; cBot = cb
 --   bathymetryFile = bathymetry_file = seafloorFile = seafloor_file; frictionMu = mu = frictionCoefficient
diff --git a/python/cabledyn/deck_file.py b/python/cabledyn/deck_file.py
index c2a54e1..5ba451b 100644
--- a/python/cabledyn/deck_file.py
+++ b/python/cabledyn/deck_file.py
@@ -12,6 +12,7 @@
 
 from __future__ import annotations
 
+import bisect
 import contextlib
 import copy
 import hashlib
@@ -22,10 +23,11 @@
 import re
 import sys
 import tempfile
-from collections.abc import Callable, Mapping
+from collections.abc import Callable, Iterable, Iterator, Mapping, Sequence
 from dataclasses import dataclass, field
 from pathlib import Path, PureWindowsPath
 from types import MappingProxyType
+from typing import Any
 
 from cabledyn._paths import windows_device_component
 from cabledyn.errors import DeckFormatError
@@ -277,6 +279,8 @@
     "bodysubstep",
     "dtwave",
     "wvhicoff",
+    "waverefx",
+    "waverefy",
     "bodywetting",
     "bodyhydro",
     "rodhydro",
@@ -1229,6 +1233,8 @@ def _scalar_option_value_error(keyword: str, value: str) -> str | None:
         "dtout",
         "dtwave",
         "wvhicoff",
+        "waverefx",
+        "waverefy",
         "inertialf",
         "inertialf_rampt",
     }:
@@ -1441,6 +1447,911 @@ def _scalar_option_value_error(keyword: str, value: str) -> str | None:
     return None
 
 
+# --------------------------------------------------------------------------- ancillary files
+# The native readers of the files a deck names, as far as they refuse a deck. Each function
+# follows the subroutine it names in src/CableDyn_DeckDriver.f90 (or the module given) and raises
+# _AncillaryError with the native message; DeckFile adds the deck's label.
+
+_NATIVE_RECORD_LENGTH = 512  # DECK_RECLEN: the longest non-comment record of next_record
+_NATIVE_TABLE_ROW_LENGTH = 65536  # LONG_ROW of read_moordyn_c_kinematics
+_NATIVE_EPSILON = sys.float_info.epsilon
+_NUMERIC_SHAPE = re.compile(
+    r"[+-]?(?:(?i:inf|infinity|nan)|(?:[0-9]+\.?[0-9]*|\.[0-9]+)(?:[eEdD][+-]?[0-9]+)?)"
+)
+_MAX_GRID_STEPS = float(2**31 - 2)
+_MAX_MOTION_SPEED = 1.0e4
+_MAX_MOTION_ACCEL = 1.0e6
+_MAX_FLUID_SPEED = 1.0e3
+_MAX_WAVE_SAMPLES = 1.0e6
+_MAX_CURRENT_SAMPLES = 2000000
+_SYROPE_CURVE_SAMPLES = 30  # CD_SYROPE_NWC
+_SYROPE_SCAN_SAMPLES = 200  # CD_SYROPE_NSCAN
+_BATHY_CREASE_FRACTION = 0.05
+_BATHY_CREASE_MIN_HEIGHT = 1.0e-6
+_SEA_DEFAULT_DIRECTIONS = 9
+_SPREADING_QUADRATURE = 128
+_STATIC_CHECK_DT = 1.0  # the nominal dtM of a standalone TMax = 0 without one (parse_deck)
+
+
+class _AncillaryError(ValueError):
+    """A refusal of a native ancillary-file reader."""
+
+
+def _token_real(token: str) -> float | None:
+    """Mirror the native ``parse_real_token``: the finite, normal real of a whole token."""
+    if _NUMERIC_SHAPE.fullmatch(token) is None:
+        return None
+    try:
+        value = float(re.sub("[dD]", "e", token))
+    except ValueError:
+        return None
+    if not math.isfinite(value) or _is_subnormal(value):
+        return None
+    return value
+
+
+def _real_row(tokens: list[str], count: int) -> tuple[int, list[float]]:
+    """Mirror the native ``parse_real_row``: status 0, 1 (no number) or 2 (malformed), values."""
+    if not tokens or _NUMERIC_SHAPE.fullmatch(tokens[0]) is None:
+        return 1, []
+    values: list[float] = []
+    for token in tokens[:count]:
+        value = _token_real(token)
+        if value is None:
+            return 2, []
+        values.append(value)
+    return (0, values) if len(values) == count else (2, [])
+
+
+def _data_rows(records: list[str], what: str) -> Iterator[tuple[int, list[str]]]:
+    """Mirror the native ``next_record``: non-blank rows without comments, by file line."""
+    for number, raw in enumerate(records, start=1):
+        visible = _strip_inline_comment(raw)
+        if not visible.strip(_WHITESPACE):
+            continue
+        defect = _character_error(visible)
+        if defect is not None:
+            raise _AncillaryError(f"{what} line {number} {defect}")
+        if len(visible) > _NATIVE_RECORD_LENGTH:
+            raise _AncillaryError(f"{what} line {number} is longer than the record limit")
+        yield number, visible.split()
+
+
+def _increasing(values: Sequence[float]) -> bool:
+    return all(later > earlier for earlier, later in itertools.pairwise(values))
+
+
+def _padded_fft_size(wanted: int) -> int:
+    """Mirror ``padded_fft_size``: twice the next 7-smooth half length."""
+    half = max((wanted + 1) // 2, 1)
+    while True:
+        rest = half
+        for prime in (2, 3, 5, 7):
+            while rest % prime == 0:
+                rest //= prime
+        if rest == 1:
+            return 2 * half
+        half += 1
+
+
+def _check_elevation_record(
+    records: list[str], dt_wave: float, window: float, cutoff: float, depth: float | None
+) -> None:
+    """Mirror ``read_wave_elevation_components`` and the component bounds of validation."""
+    import numpy as np
+
+    times: list[float] = []
+    heights: list[float] = []
+    for raw in records:
+        tokens = raw.split()
+        if not tokens:
+            continue
+        status, values = _real_row(tokens, 2)
+        if status == 2:
+            raise _AncillaryError(f'WaterKin WaveKinFile row "{raw}" is not "time elevation"')
+        if status != 0:
+            if times:
+                raise _AncillaryError("WaterKin WaveKinFile contains a non-data row in its history")
+            continue
+        times.append(values[0])
+        heights.append(values[1])
+    if len(times) < 4:
+        raise _AncillaryError("WaterKin WaveKinFile needs at least four `time elevation` rows")
+    tolerance = 1000.0 * _NATIVE_EPSILON * max(1.0, abs(times[-1]))
+    if abs(times[0]) > tolerance:
+        raise _AncillaryError("WaterKin WaveKinFile time history must start at t = 0")
+    if not _increasing(times):
+        raise _AncillaryError("WaterKin WaveKinFile times must be strictly increasing")
+    span = window if window > 0.0 else times[-1]
+    if not span / dt_wave <= _MAX_WAVE_SAMPLES:
+        raise _AncillaryError("WaterKin TMax/dtWave needs more than 1000000 wave samples")
+    wanted = math.ceil(span / dt_wave)
+    samples = _padded_fft_size(wanted)
+    series = np.zeros(samples)
+    series[:wanted] = np.interp(np.arange(wanted) * dt_wave, times, heights, right=0.0)
+    kept = samples // 2
+    if cutoff >= 0.0:
+        step = 2.0 * math.pi / (samples * dt_wave)
+        kept = min(kept, math.floor(cutoff / step * (1.0 + 1.0e-15)))
+    if kept < 1:
+        raise _AncillaryError("WvHiCOff keeps no component of the WaterKin WaveKinFile record")
+    mean = float(series.sum() / samples)
+    amplitudes = np.abs(np.fft.rfft(series - mean))[1 : kept + 1] * (2.0 / samples)
+    if samples % 2 == 0 and kept == samples // 2:
+        amplitudes[-1] *= 0.5
+    top = 2.0 * math.pi * kept / (samples * dt_wave)
+    if float(amplitudes.max()) > 1.0e3 or top > 1.0e4:
+        raise _AncillaryError(
+            "the wave components (WaterKin WaveKinMod 1 or WaveKin 3/7) need amplitudes of at "
+            "most 1e3 m and frequencies of at most 1e4 rad/s"
+        )
+    if depth is not None and depth + mean <= 0.0:
+        raise _AncillaryError("WaterKin mean elevation leaves a non-positive effective depth")
+
+
+def _check_waterkin_current(records: list[str]) -> None:
+    """Mirror the CurrentMod 1 table of ``read_waterkin_records`` (records after CurrentMod)."""
+    if len(records) < 2:
+        raise _AncillaryError("WaterKin file ended inside the current-table header")
+    malformed = "malformed WaterKin current row (expected z ux uy as plain numbers)"
+    rows: list[list[float]] = []
+    start = 2
+    for offset in range(4):
+        if start + offset >= len(records):
+            break
+        status, values = _real_row(records[start + offset].split(), 3)
+        if status == 2:
+            raise _AncillaryError(malformed)
+        if status == 0:
+            rows.append(values)
+            start += offset + 1
+            break
+    if not rows:
+        raise _AncillaryError("could not find the WaterKin current profile table")
+    for record in records[start:]:
+        status, values = _real_row(record.split(), 3)
+        if status == 2:
+            raise _AncillaryError(malformed)
+        if status != 0:
+            break
+        rows.append(values)
+    if len(rows) < 2:
+        raise _AncillaryError("WaterKin CurrentMod 1 needs at least two `z ux uy` depth rows")
+    depths = [row[0] for row in rows]
+    if any(value > 0.0 for value in depths):
+        if not all(value >= 0.0 for value in depths):
+            raise _AncillaryError("WaterKin current depth column mixes signs")
+        depths = [-value for value in depths]
+    if depths[0] > depths[-1]:
+        depths.reverse()
+    if any(b - a <= 100.0 * _NATIVE_EPSILON for a, b in itertools.pairwise(depths)):
+        raise _AncillaryError("WaterKin current depths must be strictly monotonic")
+    _check_current_bounds(depths, [value for row in rows for value in row[1:]])
+
+
+def _check_current_bounds(depths: Sequence[float], velocities: Iterable[float]) -> None:
+    if any(abs(value) > _NATIVE_MAX_COORDINATE for value in depths) or any(
+        abs(value) > _MAX_FLUID_SPEED for value in velocities
+    ):
+        raise _AncillaryError(
+            "current profile depths must be at most 1e6 m and velocities at most 1e3 m/s in "
+            "magnitude"
+        )
+
+
+def _wave_frequency_rows(records: list[str], label: str) -> list[list[float]]:
+    """Mirror ``read_table``: rows `omega Re Im beta`, the direction zero when left out."""
+    rows: list[list[float]] = []
+    for number, record in enumerate(records, start=1):
+        tokens = record.split()
+        if not tokens:
+            continue
+        status, values = _real_row(tokens, 4)
+        if status != 0:
+            status, values = _real_row(tokens, 3)
+            values.append(0.0)
+        if status != 0:
+            raise DeckFormatError(f"{label} line {number}: row is not 3 to 4 plain numbers")
+        rows.append(values)
+    return rows
+
+
+def _table_tokens(record: str) -> list[str]:
+    """Mirror ``row_tokens``: the tokens up to the first one holding a ``#``."""
+    tokens: list[str] = []
+    for token in record.split():
+        if "#" in token:
+            break
+        tokens.append(token)
+    return tokens
+
+
+def _table_number(token: str, name: str, number: int) -> float:
+    value = _token_real(token)
+    if value is None:
+        raise _AncillaryError(f'{name} line {number}: value "{token}" is not a plain finite number')
+    return value
+
+
+def _check_current_profile(records: list[str], name: str) -> None:
+    """Mirror ``Currents 1``: `z ux uy uz` rows after their header rows."""
+    rows: list[list[float]] = []
+    for number, raw in enumerate(records, start=1):
+        tokens = raw.split()
+        if not tokens:
+            continue
+        status, values = _real_row(tokens, 4)
+        if status != 0:
+            status, values = _real_row(tokens, 3)
+        if status == 1 and not rows:
+            continue
+        if status != 0:
+            raise _AncillaryError(f"{name} line {number}: row is not 3 to 4 plain numbers")
+        rows.append(values)
+    if not rows:
+        raise _AncillaryError(f"Currents 1: {name} has no `z ux uy uz` rows")
+    depths = [row[0] for row in rows]
+    if not _increasing(depths):
+        raise _AncillaryError(f"Currents 1: {name} depths must be strictly increasing")
+    if any(b - a <= 100.0 * _NATIVE_EPSILON for a, b in itertools.pairwise(depths)):
+        raise _AncillaryError("current profile depths must be strictly increasing")
+    _check_current_bounds(depths, [value for row in rows for value in row[1:]])
+
+
+def _check_current_resample(
+    times: list[float], spatial: int, velocities: Any, coordinates: Iterable[float]
+) -> None:
+    """Mirror ``CD_Current_Grid_Resample`` and the grid bounds; velocities is (nt, points)."""
+    import numpy as np
+
+    step = min(b - a for a, b in itertools.pairwise(times))
+    if times[-1] < 0.0:
+        raise _AncillaryError("current grid: the last time must not be negative")
+    if not times[-1] / _MAX_CURRENT_SAMPLES < step:
+        raise _AncillaryError(
+            "current grid: the last time over the smallest time interval exceeds the sample limit"
+        )
+    samples = math.floor(times[-1] / step) + 1
+    if samples * spatial > _MAX_CURRENT_SAMPLES:
+        raise _AncillaryError("current grid: more than 2e6 resampled space-time samples")
+    at = np.arange(samples) * step
+    upper = np.clip(np.searchsorted(times, at, side="right"), 1, len(times) - 1)
+    lower = upper - 1
+    knots = np.asarray(times)
+    weight = np.clip((at - knots[lower]) / (knots[upper] - knots[lower]), 0.0, 1.0)[:, None]
+    resampled = (1.0 - weight) * velocities[lower] + weight * velocities[upper]
+    if any(abs(value) > _NATIVE_MAX_COORDINATE for value in coordinates) or (
+        float(np.abs(resampled).max()) > _MAX_FLUID_SPEED
+    ):
+        raise _AncillaryError(
+            "current grid coordinates must be at most 1e6 m and velocities at most 1e3 m/s in "
+            "magnitude"
+        )
+
+
+def _check_dynamic_current(records: list[str], name: str) -> None:
+    """Mirror ``read_dynamic_profile`` (Currents 2)."""
+    import numpy as np
+
+    depths: list[float] = []
+    times: list[float] = []
+    rows: list[list[float]] = []
+    for number, raw in enumerate(records, start=1):
+        if len(raw) > _NATIVE_TABLE_ROW_LENGTH:
+            raise _AncillaryError(f"{name} line {number}: row is too long")
+        if number <= 4 or number == 6:
+            continue
+        tokens = _table_tokens(raw)
+        if number == 5:
+            if not tokens:
+                raise _AncillaryError(f"Currents 2: line 5 of {name} must list the depths")
+            for token in tokens:
+                depth = float(np.float32(_table_number(token, name, number)))
+                if not math.isfinite(depth):
+                    raise _AncillaryError(f"{name} depth is not a single-precision number")
+                depths.append(depth)
+            continue
+        if not tokens:
+            continue
+        count = len(depths)
+        if len(tokens) not in (1 + count, 1 + 2 * count, 1 + 3 * count):
+            raise _AncillaryError(f"Currents 2: {name} line {number}: wrong number of values")
+        values = [_table_number(token, name, number) for token in tokens]
+        times.append(values[0])
+        rows.append(values[1:] + [0.0] * (3 * count - len(values) + 1))
+    if not depths or len(records) < 7:
+        raise _AncillaryError(f"Currents 2: {name} needs at least 7 lines")
+    if len(times) < 2:
+        raise _AncillaryError(f"Currents 2: {name} needs at least two time rows")
+    if not _increasing(depths):
+        raise _AncillaryError(f"Currents 2: {name} depths must be strictly increasing")
+    if not _increasing(times):
+        raise _AncillaryError(f"Currents 2: {name} times must be strictly increasing")
+    _check_current_resample(times, len(depths), np.asarray(rows), [*depths, 0.0])
+
+
+def _check_grid_current(records: list[str], name: str) -> None:
+    """Mirror ``read_4d_grid`` (Currents 5)."""
+    import numpy as np
+
+    sizes: list[int] = []
+    axes: list[list[float]] = []
+    points: dict[tuple[int, int, int, int], list[float]] = {}
+    for number, raw in enumerate(records, start=1):
+        if len(raw) > _NATIVE_TABLE_ROW_LENGTH:
+            raise _AncillaryError(f"{name} line {number}: row is too long")
+        tokens = _table_tokens(raw)
+        if number == 1:
+            good = len(tokens) == 4 and all(re.fullmatch(r"[+-]?[0-9]+", t) for t in tokens)
+            sizes = [int(token) for token in tokens] if good else []
+            if not good or any(size < 1 for size in sizes):
+                raise _AncillaryError(f"Currents 5: line 1 of {name} must be `nx ny nz nt`")
+            if math.prod(sizes) > _MAX_CURRENT_SAMPLES:
+                raise _AncillaryError(f"Currents 5: {name} declares more than 2e6 grid points")
+        elif number <= 5:
+            if len(tokens) != sizes[number - 2]:
+                raise _AncillaryError(f"Currents 5: line {number} of {name} has the wrong length")
+            axes.append([_table_number(token, name, number) for token in tokens])
+        else:
+            if not tokens:
+                continue
+            if number == 6 and _NUMERIC_SHAPE.fullmatch(tokens[0]) is None:
+                continue
+            if len(tokens) != 7:
+                raise _AncillaryError(f"Currents 5: {name} line {number}: a row is 7 numbers")
+            row = [_table_number(token, name, number) for token in tokens]
+            where: list[int] = []
+            for axis, value in zip(axes, row, strict=False):
+                if value not in axis:
+                    raise _AncillaryError(f"Currents 5: {name} line {number}: off the grid")
+                where.append(axis.index(value))
+            key = (where[0], where[1], where[2], where[3])
+            if key in points:
+                raise _AncillaryError(f"Currents 5: {name} line {number}: repeats a grid point")
+            points[key] = row[4:]
+    if len(axes) < 4:
+        raise _AncillaryError(f"Currents 5: {name} needs the grid sizes and coordinates")
+    if not all(_increasing(axis) for axis in axes):
+        raise _AncillaryError(f"Currents 5: {name} coordinates must be strictly increasing")
+    if len(axes[3]) < 2:
+        raise _AncillaryError(f"Currents 5: {name} needs at least two times")
+    if len(points) < math.prod(sizes):
+        raise _AncillaryError(f"Currents 5: {name} has no row for a grid point")
+    spatial = math.prod(sizes[:3])
+    field = np.zeros((sizes[3], spatial * 3))
+    for (ix, iy, iz, it), velocity in points.items():
+        column = 3 * ((ix * sizes[1] + iy) * sizes[2] + iz)
+        field[it, column : column + 3] = velocity
+    _check_current_resample(axes[3], spatial, field, [*axes[0], *axes[1], *axes[2]])
+
+
+def _motion_time_index(
+    time: float, dt: float, n_time: int, tolerance: float, start_row_only: bool
+) -> int | None:
+    """Mirror ``motion_grid_index``: the 0-based step, -1 past TMax, ``None`` off the grid.
+
+    start_row_only: a static check without dtM, whose grid is the t = 0 row: a later row is not
+    used, whatever its time step.
+    """
+    if start_row_only and time > tolerance:
+        return -1
+    if not abs(time / dt) < _MAX_GRID_STEPS:
+        return None
+    step = round(time / dt)
+    if step < 0:
+        return None
+    if step < n_time:
+        return step if abs(step * dt - time) <= tolerance else None
+    if abs(step * dt - time) <= 100.0 * _NATIVE_EPSILON * max(1.0, abs(time)):
+        return -1
+    return None
+
+
+def _check_motion_rows(
+    records: list[str],
+    dt: float,
+    tmax: float,
+    eligible: set[int] | None,
+    roll_ends: set[int] | None,
+    turbines: dict[int, list[int]] | None,
+    start_row_only: bool,
+) -> None:
+    """Mirror ``read_motion_file``.
+
+    eligible: the ids the file must move at every time (``None``: on a route that takes
+    other points as well, their rows are only read). roll_ends: with a line solved in torsion,
+    the points whose moving end is restrained in torsion. turbines: the point ids of each
+    TURBINES row, on a deck that has them.
+    """
+    n_time = round(tmax / dt) + 1
+    tolerance = 100.0 * _NATIVE_EPSILON * max(1.0, tmax)
+    what = "motionFile"
+    seen: set[tuple[int, int]] = set()
+    roll_rows: dict[int, dict[int, float]] = {}
+    for number, tokens in _data_rows(records, what):
+        where = f"{what} line {number}"
+        if turbines is not None:
+            status, row = _real_row(tokens, 21)
+            if status != 0:
+                raise _AncillaryError(f"{where}: TURBINES motion rows must be 21 numbers")
+            if (
+                any(abs(value) > _NATIVE_MAX_COORDINATE for value in row[2:9])
+                or any(abs(value) > _MAX_MOTION_SPEED for value in row[9:15])
+                or any(abs(value) > _MAX_MOTION_ACCEL for value in row[15:21])
+            ):
+                raise _AncillaryError(f"{where}: value outside the admissible range")
+            turbine = int(row[1]) if row[1] == int(row[1]) else 0
+            if turbine not in turbines:
+                raise _AncillaryError(f"{where}: J must be the integer id of a TURBINES row")
+            if abs(math.sqrt(sum(value * value for value in row[5:9])) - 1.0) > 1.0e-6:
+                raise _AncillaryError(f"{where}: the quaternion must have unit norm (within 1e-6)")
+            step = _motion_time_index(row[0], dt, n_time, tolerance, start_row_only)
+            if step is None:
+                raise _AncillaryError(f"{where}: time is outside the dtM/TMax grid")
+            if step == -1:
+                continue
+            if (-turbine, step) in seen:
+                raise _AncillaryError(f"{where}: duplicate row for turbine {turbine}")
+            seen.add((-turbine, step))
+            seen.update((point, step) for point in turbines[turbine])
+            continue
+        status, row = _real_row(tokens, 11)
+        if status != 0:
+            raise _AncillaryError(f"{where}: rows must be: time point_id x y z vx vy vz ax ay az")
+        if (
+            any(abs(value) > _NATIVE_MAX_COORDINATE for value in row[2:5])
+            or any(abs(value) > _MAX_MOTION_SPEED for value in row[5:8])
+            or any(abs(value) > _MAX_MOTION_ACCEL for value in row[8:11])
+        ):
+            raise _AncillaryError(f"{where}: value outside the admissible range")
+        point = int(row[1]) if abs(row[1]) <= 2**31 - 1 and row[1] == int(row[1]) else 0
+        if point < 1:
+            raise _AncillaryError(f"{where}: point_id must be a positive integer")
+        step = _motion_time_index(row[0], dt, n_time, tolerance, start_row_only)
+        if step is None:
+            raise _AncillaryError(f"{where}: time is outside the dtM/TMax grid")
+        if eligible is not None and point not in eligible:
+            raise _AncillaryError(f"{where}: references point id {point}, which is not eligible")
+        if step == -1:
+            continue
+        if (point, step) in seen:
+            raise _AncillaryError(f"{where}: duplicate row for point {point} at this time")
+        seen.add((point, step))
+        if roll_ends is not None and len(tokens) >= 12 and _NUMERIC_SHAPE.fullmatch(tokens[11]):
+            roll = _token_real(tokens[11])
+            if roll is None:
+                raise _AncillaryError(f"{where}: the roll column must be a finite number")
+            roll_rows.setdefault(point, {})[step] = roll
+    if eligible is not None and any(
+        (point, step) not in seen for point in eligible for step in range(n_time)
+    ):
+        raise _AncillaryError(
+            "motionFile must provide every Coupled/Vessel point at every dtM time"
+        )
+    if roll_ends is None:
+        return
+    for point, rolls in roll_rows.items():
+        if len(rolls) != n_time:
+            raise _AncillaryError(
+                f"motionFile: point {point} gives the roll column on some rows only"
+            )
+        if any(rolls.values()):
+            if point not in roll_ends:
+                raise _AncillaryError(
+                    f"motionFile: point {point} has a non-zero roll column, but no line restrained "
+                    "in torsion has its moving end there"
+                )
+            if rolls[0] != 0.0:
+                raise _AncillaryError(
+                    f"motionFile: the roll column of point {point} must be 0 at t = 0"
+                )
+
+
+def _check_vessel_rows(records: list[str], dt: float, tmax: float, start_row_only: bool) -> None:
+    """Mirror ``read_vessel_record``."""
+    n_time = round(tmax / dt) + 1
+    tolerance = 100.0 * _NATIVE_EPSILON * max(1.0, tmax)
+    what = "vesselMotion file"
+    form = 0
+    seen: set[int] = set()
+    for number, tokens in _data_rows(records, what):
+        where = f"{what} line {number}"
+        if form == 0 and len(tokens) in (19, 20):
+            form = len(tokens)
+        status, row = (1, []) if len(tokens) != form else _real_row(tokens, form)
+        if status != 0:
+            raise _AncillaryError(f"{where}: rows must hold 19 or 20 values, one form for the file")
+        if any(abs(value) > _NATIVE_MAX_COORDINATE for value in row[1:]):
+            raise _AncillaryError(f"{where}: value outside the admissible range")
+        step = _motion_time_index(row[0], dt, n_time, tolerance, start_row_only)
+        if step is None:
+            raise _AncillaryError(f"{where}: time is outside the dtM/TMax grid")
+        if step == -1:
+            continue
+        if step in seen:
+            raise _AncillaryError(f"{where}: duplicate row at this time")
+        if form == 20 and abs(math.sqrt(sum(v * v for v in row[4:8])) - 1.0) > 1.0e-6:
+            raise _AncillaryError(f"{where}: the quaternion must have unit norm (within 1e-6)")
+        seen.add(step)
+    if len(seen) < n_time:
+        raise _AncillaryError("vesselMotion file must provide one row at every dtM time")
+
+
+def _check_rao_table(records: list[str], headings_in_use: Iterable[float]) -> None:
+    """Mirror ``read_rao_table``, ``CD_Vessel_RAO_Set`` and the heading test of ``_RAO_Eval``."""
+    what = "vesselRAO file"
+    heads: list[float] = []
+    blocks: list[list[float]] = []
+    amplitudes: list[float] = []
+    for number, tokens in _data_rows(records, what):
+        where = f"{what} line {number}"
+        if tokens[0].lower() == "heading":
+            value = _token_real(tokens[1]) if len(tokens) == 2 else None
+            if value is None or len(heads) >= 361:
+                raise _AncillaryError(f'{where}: a heading row must be "HEADING <deg>"')
+            heads.append(value)
+            blocks.append([])
+            continue
+        status, row = (1, []) if len(tokens) != 13 else _real_row(tokens, 13)
+        if status != 0 or not heads or len(blocks[-1]) >= 2000:
+            raise _AncillaryError(f'{where}: data rows follow a "HEADING <deg>" row with 13 values')
+        blocks[-1].append(row[0])
+        amplitudes.extend(row[1::2])
+    if not heads:
+        raise _AncillaryError('vesselRAO file holds no "HEADING <deg>" block')
+    if any(not block for block in blocks):
+        raise _AncillaryError("a vesselRAO heading block holds no period rows")
+    if any(block != blocks[0] for block in blocks):
+        raise _AncillaryError("every vesselRAO heading block must list the same periods")
+    periods = sorted(blocks[0])
+    heads = sorted(heads)
+    if periods[0] <= 0.0:
+        raise _AncillaryError("vesselRAO file: RAO periods must be positive")
+    if any(value < 0.0 for value in amplitudes):
+        raise _AncillaryError("vesselRAO file: RAO amplitudes must be non-negative")
+    if not _increasing(periods):
+        raise _AncillaryError("vesselRAO file: RAO periods must be strictly ascending")
+    if not _increasing(heads):
+        raise _AncillaryError("vesselRAO file: RAO headings must be strictly ascending")
+    if heads[-1] - heads[0] >= 360.0:
+        raise _AncillaryError("vesselRAO file: RAO headings must lie within one 360 deg turn")
+    for heading in headings_in_use:
+        turned = heads[0] + (heading - heads[0]) % 360.0
+        if turned > heads[-1] + 1.0e-9 and abs(turned - 360.0 - heads[0]) <= 1.0e-9:
+            turned = heads[0]
+        if len(heads) == 1:
+            if abs(turned - heads[0]) > 1.0e-9:
+                raise _AncillaryError("relative wave heading is not in the one-heading RAO table")
+        elif turned > heads[-1] + 1.0e-9:
+            raise _AncillaryError("relative wave heading lies outside the RAO table heading range")
+
+
+def _floor_axis(axis: list[float], q: float) -> tuple[int, float, bool]:
+    """Mirror ``locate_axis`` (0-based): the cell, the fraction and whether q is off the grid."""
+    if q <= axis[0]:
+        return 0, 0.0, q < axis[0]
+    if q >= axis[-1]:
+        return len(axis) - 2, 1.0, q > axis[-1]
+    cell = bisect.bisect_left(axis, q) - 1
+    return cell, (q - axis[cell]) / (axis[cell + 1] - axis[cell]), False
+
+
+def _crease_band(axis: list[float], cell: int, q: float) -> tuple[int, float, float]:
+    """Mirror ``crease_band`` (0-based; -1: none)."""
+    if cell > 0:
+        band = _BATHY_CREASE_FRACTION * min(
+            axis[cell] - axis[cell - 1], axis[cell + 1] - axis[cell]
+        )
+        if q - axis[cell] < band:
+            return cell, q - axis[cell], band
+    if cell + 2 < len(axis):
+        band = _BATHY_CREASE_FRACTION * min(
+            axis[cell + 1] - axis[cell], axis[cell + 2] - axis[cell + 1]
+        )
+        if axis[cell + 1] - q < band:
+            return cell + 1, q - axis[cell + 1], band
+    return -1, 0.0, 0.0
+
+
+def _bathymetry_floor(
+    xs: list[float], ys: list[float], depth: dict[tuple[int, int], float], x: float, y: float
+) -> tuple[float, float]:
+    """Mirror ``floor_eval``: the rounded and the bilinear seabed elevation at (x, y)."""
+    ix, tx, clamp_x = _floor_axis(xs, x)
+    iy, ty, clamp_y = _floor_axis(ys, y)
+    grid = (
+        -(1.0 - tx) * (1.0 - ty) * depth[ix, iy]
+        - tx * (1.0 - ty) * depth[ix + 1, iy]
+        - (1.0 - tx) * ty * depth[ix, iy + 1]
+        - tx * ty * depth[ix + 1, iy + 1]
+    )
+    if len(xs) <= 2 and len(ys) <= 2:
+        return grid, grid
+
+    def x_jump(k: int, j: int, band: float) -> tuple[float, bool]:
+        jump = (depth[k, j] - depth[k + 1, j]) / (xs[k + 1] - xs[k]) - (
+            depth[k - 1, j] - depth[k, j]
+        ) / (xs[k] - xs[k - 1])
+        active = abs(jump) * band / 6.0 > _BATHY_CREASE_MIN_HEIGHT
+        return (jump if active else 0.0), active
+
+    def y_jump(i: int, k: int, band: float) -> tuple[float, bool]:
+        jump = (depth[i, k] - depth[i, k + 1]) / (ys[k + 1] - ys[k]) - (
+            depth[i, k - 1] - depth[i, k]
+        ) / (ys[k] - ys[k - 1])
+        active = abs(jump) * band / 6.0 > _BATHY_CREASE_MIN_HEIGHT
+        return (jump if active else 0.0), active
+
+    kx, u, band_x = (-1, 0.0, 0.0) if clamp_x else _crease_band(xs, ix, x)
+    ky, v, band_y = (-1, 0.0, 0.0) if clamp_y else _crease_band(ys, iy, y)
+    if kx < 0 and ky < 0:
+        return grid, grid
+    floor = grid
+    hx = 0.0
+    if kx >= 0:
+        hx = (band_x - abs(u)) ** 3 / (6.0 * band_x * band_x)
+        jx0, on0 = x_jump(kx, iy, band_x)
+        jx1, on1 = x_jump(kx, iy + 1, band_x)
+        if on0 or on1:
+            floor += hx * ((1.0 - ty) * jx0 + ty * jx1)
+    if ky >= 0:
+        hy = (band_y - abs(v)) ** 3 / (6.0 * band_y * band_y)
+        jy0, on0 = y_jump(ix, ky, band_y)
+        jy1, on1 = y_jump(ix + 1, ky, band_y)
+        mixed = 0.0
+        on_mixed = False
+        if kx >= 0:
+            below, _ = x_jump(kx, ky - 1, band_x)
+            centre, _ = x_jump(kx, ky, band_x)
+            above, _ = x_jump(kx, ky + 1, band_x)
+            mixed = (above - centre) / (ys[ky + 1] - ys[ky]) - (centre - below) / (
+                ys[ky] - ys[ky - 1]
+            )
+            on_mixed = abs(mixed) * band_x * band_y / 36.0 > _BATHY_CREASE_MIN_HEIGHT
+            if not on_mixed:
+                mixed = 0.0
+        if on0 or on1 or on_mixed:
+            floor += hy * ((1.0 - tx) * jy0 + tx * jy1 + hx * mixed)
+    return floor, grid
+
+
+def _check_bathymetry(
+    records: list[str], fixed_points: list[tuple[int, float, float, float]]
+) -> None:
+    """Mirror ``load_bathymetry_file`` and ``check_fixed_points_on_bathymetry``."""
+    what = "bathymetry file"
+    rows: list[list[float]] = []
+    for number, tokens in _data_rows(records, what):
+        if len(tokens) != 3:
+            raise _AncillaryError(f'{what} line {number}: rows must be "x y depth"')
+        status, row = _real_row(tokens, 3)
+        if status != 0:
+            raise _AncillaryError(f"{what} line {number}: malformed row")
+        rows.append(row)
+    if len(rows) < 4:
+        raise _AncillaryError("bathymetry file needs at least a 2x2 grid")
+    if any(row[2] <= 0.0 for row in rows):
+        raise _AncillaryError("bathymetry file depths must be positive")
+
+    def unique(values: list[float]) -> list[float]:
+        kept: list[float] = []
+        for value in sorted(values):
+            if not kept or abs(value - kept[-1]) > 16.0 * _NATIVE_EPSILON * max(
+                1.0, abs(value), abs(kept[-1])
+            ):
+                kept.append(value)
+        return kept
+
+    def nearest(axis: list[float], value: float) -> int:
+        at = bisect.bisect_left(axis, value)
+        return min(
+            (i for i in (at - 1, at) if 0 <= i < len(axis)), key=lambda i: abs(axis[i] - value)
+        )
+
+    xs = unique([row[0] for row in rows])
+    ys = unique([row[1] for row in rows])
+    if len(xs) < 2 or len(ys) < 2 or len(rows) != len(xs) * len(ys):
+        raise _AncillaryError("bathymetry file must define one complete rectangular grid")
+    depth: dict[tuple[int, int], float] = {}
+    for x, y, value in rows:
+        key = (nearest(xs, x), nearest(ys, y))
+        if key in depth:
+            raise _AncillaryError("bathymetry file has duplicate x/y entries")
+        depth[key] = value
+    for point, x, y, z in fixed_points:
+        floor, grid = _bathymetry_floor(xs, ys, depth, x, y)
+        tolerance = max(1.0e-3, 1.0e-5 * abs(floor)) * (1.0 + 1.0e-9)
+        if abs(z - floor) <= tolerance or abs(z - grid) <= tolerance or z > floor:
+            continue
+        raise _AncillaryError(f"POINT {point} (Fixed) lies below the seabed: z = {z:g}")
+
+
+def _read_syrope_settings(records: list[str]) -> tuple[str, str, float, float]:
+    """Mirror the settings rows of ``read_syrope_settings``: OWC path, WCType, k1, k2."""
+    what = "Syrope settings file"
+    found: dict[str, str] = {}
+    for number, tokens in _data_rows(records, what):
+        if len(tokens) < 2:
+            continue
+        name = tokens[1].lower()
+        value = _strip_quotes(tokens[0])
+        if name not in {"owc", "wctype", "k1", "k2"}:
+            continue
+        if name in found:
+            raise _AncillaryError(f"{what} line {number}: duplicate Syrope {tokens[1]} row")
+        if name == "wctype" and value.lower() not in {"linear", "quadratic", "exp"}:
+            raise _AncillaryError(
+                f"{what} line {number}: Syrope WCType must be LINEAR, QUADRATIC, or EXP"
+            )
+        if name in {"k1", "k2"} and _token_real(value) is None:
+            raise _AncillaryError(f"{what} line {number}: Syrope {tokens[1]} is not a plain number")
+        found[name] = value
+    if len(found) < 4:
+        raise _AncillaryError("the Syrope settings file needs OWC, WCType, k1, and k2 rows")
+    k1 = _token_real(found["k1"])
+    k2 = _token_real(found["k2"])
+    assert k1 is not None and k2 is not None
+    return found["owc"], found["wctype"].lower(), k1, k2
+
+
+def _read_syrope_table(records: list[str]) -> list[list[float]]:
+    """Mirror the OWC table rows of ``read_syrope_settings``."""
+    what = "Syrope OWC table"
+    rows: list[list[float]] = []
+    for number, tokens in _data_rows(records, what):
+        if not any(_NUMERIC_SHAPE.fullmatch(token) for token in tokens):
+            continue
+        status, row = _real_row(tokens, 2)
+        if status != 0:
+            raise _AncillaryError(f"{what} line {number}: malformed row")
+        rows.append(row)
+    if len(rows) < 2:
+        raise _AncillaryError("the Syrope OWC table needs at least two rows")
+    return rows
+
+
+def _syrope_curve_defect(
+    strain: list[float],
+    tension: list[float],
+    shape: str,
+    k1: float,
+    k2: float,
+    alpha: float,
+    beta: float,
+    t_max: float,
+) -> bool:
+    """True when ``CD_Syrope_Working_Curve`` fails at the running maximum ``t_max``."""
+
+    def fast(value: float) -> float:
+        return math.log(1.0 + beta / alpha * value) / beta
+
+    if not t_max > 0.0 or t_max < tension[0] or t_max > tension[-1]:
+        return True
+    at = bisect.bisect_left(tension, t_max)
+    if at == 0:
+        eps_max = strain[0]
+    elif at == len(tension):
+        eps_max = strain[-1]
+    else:
+        eps_max = strain[at - 1] + (strain[at] - strain[at - 1]) * (t_max - tension[at - 1]) / (
+            tension[at] - tension[at - 1]
+        )
+    eps_min = strain[0] + k1 * (eps_max - strain[0])
+    if shape == "linear" and k1 >= 1.0:
+        eps_min = eps_max - t_max / k1
+    if not (eps_min >= 0.0 and eps_max - eps_min > 0.0):
+        return True
+    previous: tuple[float, float] | None = None
+    for index in range(_SYROPE_CURVE_SAMPLES):
+        at_strain = eps_min + (eps_max - eps_min) * index / (_SYROPE_CURVE_SAMPLES - 1)
+        xi = (at_strain - eps_min) / (eps_max - eps_min)
+        if shape == "linear":
+            value = t_max * xi
+        elif shape == "quadratic":
+            value = t_max * xi * (k2 * xi + (1.0 - k2))
+        else:
+            value = t_max * (1.0 - math.exp(k2 * xi)) / (1.0 - math.exp(k2))
+        argument = 1.0 + beta / alpha * value
+        if argument <= 0.0:
+            return True
+        slow = at_strain - fast(value)
+        if previous is not None and (value <= previous[0] or slow <= previous[1]):
+            return True
+        previous = (value, slow)
+    return False
+
+
+def _check_syrope_curve(
+    rows: list[list[float]],
+    shape: str,
+    k1: float,
+    k2: float,
+    alpha: float,
+    beta: float,
+) -> None:
+    """Mirror ``CD_Syrope_Init`` (the constants BA_s and BA_d are checked with the row)."""
+    strain = [row[0] for row in rows]
+    tension = [row[1] for row in rows]
+    if not (_increasing(strain) and _increasing(tension)):
+        raise _AncillaryError("the OWC strain and tension columns must be strictly increasing")
+    if shape == "exp" and not k2 > 0.0:
+        raise _AncillaryError("the EXP working curve needs a positive shape parameter p2 (k2 > 0)")
+    if shape == "linear" and k1 < 0.0:
+        raise _AncillaryError("the LINEAR working curve requires k1 >= 0")
+    if shape == "quadratic" and (k1 < 0.0 or k1 >= 1.0 or k2 <= 0.0 or k2 > 1.0):
+        raise _AncillaryError("the QUADRATIC working curve requires 0 <= k1 < 1 and 0 < k2 <= 1")
+    if shape == "exp" and (k1 < 0.0 or k1 >= 1.0 or k2 < 1.0e-6 or k2 > 700.0):
+        raise _AncillaryError("the EXP working curve requires 0 <= k1 < 1 and 1e-6 <= k2 <= 700")
+    if any(1.0 + beta / alpha * value <= 0.0 for value in tension):
+        raise _AncillaryError("the OWC slow-spring static strain is non-finite")
+    slow = [
+        e - math.log(1.0 + beta / alpha * t) / beta for e, t in zip(strain, tension, strict=True)
+    ]
+    if not _increasing(slow):
+        raise _AncillaryError("the OWC slow-spring static strain must be strictly increasing")
+    low = max(tension[0], 0.0)
+    high = tension[-1]
+    if all(
+        _syrope_curve_defect(
+            strain,
+            tension,
+            shape,
+            k1,
+            k2,
+            alpha,
+            beta,
+            low + (high - low) * k / _SYROPE_SCAN_SAMPLES,
+        )
+        for k in range(1, _SYROPE_SCAN_SAMPLES + 1)
+    ):
+        raise _AncillaryError("no running maximum in the OWC table range admits a working curve")
+
+
+def _train_direction(values: tuple[str, ...]) -> tuple[float, float]:
+    """Mirror ``CD_Sea_Parse_Train``: the direction and spreading of a wave row's values."""
+    kind = values[0].lower()
+    if kind == "airy":
+        return _native_float(values[3]), 0.0
+    index = 4 if kind == "jonswap" else 7 if kind.startswith("ochi") else 3
+    spreading = _native_float(values[index + 1]) if len(values) > index + 1 else 0.0
+    return _native_float(values[index]), spreading
+
+
+def _spreading_weights(spreading: float, bins: int) -> list[float]:
+    """Mirror ``CD_Spreading_Weights``: the cos-2s weight of each equal-angle bin."""
+    constant = math.exp(math.lgamma(spreading + 1.0) - math.lgamma(spreading + 0.5))
+    constant /= math.sqrt(math.pi)
+    weights = []
+    for j in range(bins):
+        start = -0.5 * math.pi + math.pi * j / bins
+        step = (-0.5 * math.pi + math.pi * (j + 1) / bins - start) / _SPREADING_QUADRATURE
+        weight = 0.0
+        for q in range(_SPREADING_QUADRATURE + 1):
+            density = constant * max(0.0, math.cos(start + step * q)) ** (2.0 * spreading)
+            end = q in (0, _SPREADING_QUADRATURE)
+            weight += density * (1.0 if end else 4.0 if q % 2 == 1 else 2.0)
+        weights.append(weight * step / 3.0)
+    total = sum(weights)
+    if total > 0.0:
+        weights = [weight / total for weight in weights]
+    return [weight if weight >= sys.float_info.min else 0.0 for weight in weights]
+
+
+def _spread_headings(direction: float, spreading: float, bins: int) -> list[float]:
+    """The headings of a train (``CD_Sea_Synthesise``) that carry energy.
+
+    A spread train has one heading per bin centre. A narrow spreading leaves the outer bins with
+    no weight, and ``synthesize_vessel_rao`` skips their zero-amplitude components.
+    """
+    if spreading <= 0.0:
+        return [direction]
+    return [
+        direction - 90.0 + 180.0 * (j + 0.5) / bins
+        for j, weight in enumerate(_spreading_weights(spreading, bins))
+        if weight > 0.0
+    ]
+
+
 @dataclass(frozen=True)
 class DeckRecord:
     """One parsed table row and its zero-based source-line location.
@@ -2424,12 +3335,13 @@ def validate(self) -> None:
         between rows. A deck that passes can still be rejected by the solver on
         physical grounds.
 
-        The ancillary files a deck names (motion, bathymetry, WaterKin, Syrope,
-        and the MoorDyn-C kinematics files) are not read here, so the checks
-        that depend on their contents stay with the solver: for example the
-        anchor height on a bathymetry surface, whether a WaterKin file carries
-        waves (``WaveKinMod``) for ``vesselRAO``, and whether its current
-        doubles a ``current`` row or ``Currents 1``, ``2`` or ``5``.
+        The ancillary files a deck names (WaterKin and its wave record, the
+        MoorDyn-C kinematics files, motion and vessel records, the RAO table,
+        bathymetry, and the Syrope settings and OWC table) are read from the
+        deck's folder, with the rules of their native readers. A file that is
+        not there, as for a deck held in memory under its nominal path, is not
+        an error at this stage: the deck passes these checks and the solver,
+        which must open the file, refuses it by name if it is still missing.
 
         Raises
         ------
@@ -2583,6 +3495,8 @@ def effective(key: str, default: float) -> float:
             has_waterkin_file = waterkin_value.lower() not in {"", "none", "seastate"} and any(
                 character in _WATERKIN_FILENAME_LETTERS for character in waterkin_value
             )
+            if has_waterkin_file and motion_source == "vesselrao":
+                has_waterkin_file = self._waterkin_wave_mode(waterkin_value) in {None, 1}
             has_trains = any(o.keyword.lower() in _WAVETRAIN_KEYWORDS for o in options)
             # WaveKin 3/7 installs linear component waves (native has_wave)
             if motion_source == "vesselrao" and not (
@@ -2675,6 +3589,35 @@ def effective(key: str, default: float) -> float:
                 raise DeckFormatError(
                     f"{self._label}: WaveKin {mdc_wave_kin} waves require dtM, TMax and WtrDpth"
                 )
+            if mdc_wave_kin == 7:
+                self._check_wave_frequencies(_native_float(option_by_key["water_depth"].values[0]))
+        # native validation: WaveRefX and WaveRefY place the phases of component waves alone,
+        # which a WaterKin file carries under WaveKinMod 1.
+        # The bound applies to the value in force: a later row of an option replaces an earlier one.
+        wave_ref_rows = {
+            o.keyword.lower(): o for o in options if o.keyword.lower() in {"waverefx", "waverefy"}
+        }
+        if wave_ref_rows:
+            if any(
+                abs(_native_float(o.values[0])) > _NATIVE_MAX_COORDINATE
+                for o in wave_ref_rows.values()
+            ):
+                raise DeckFormatError(
+                    f"{self._label}: WaveRefX and WaveRefY must be at most 1e6 m in magnitude"
+                )
+            waterkin_ref = option_by_key.get("waterkin")
+            waterkin_text = "" if waterkin_ref is None else waterkin_ref.values[0]
+            has_record = bool(mdc_wave_kin) or (
+                waterkin_text.lower() not in {"", "none", "seastate"}
+                and any(character in _WATERKIN_FILENAME_LETTERS for character in waterkin_text)
+                and self._waterkin_wave_mode(waterkin_text) in {None, 1}
+            )
+            if self.caller_driven or has_waves or trains or not has_record:
+                raise DeckFormatError(
+                    f"{self._label}: WaveRefX and WaveRefY place the reference point of an "
+                    "elevation record or a component file (WaterKin WaveKinMod 1, WaveKin 3 or 7); "
+                    "this deck has none"
+                )
         if mdc_currents:
             if current is not None and current.values[0].lower() != "none":
                 raise DeckFormatError(
@@ -4311,6 +5254,431 @@ def is_body_kind(kind: str, pinned: bool) -> bool:
                         )
                     continue
             raise DeckFormatError(f"{self._label}: unsupported OUTPUT channel {channel!r}")
+        try:
+            self._check_ancillary_files(
+                option_by_key,
+                options,
+                motion_source,
+                mdc_wave_kin,
+                mdc_currents,
+                has_waves or bool(trains),
+                current is not None and current.values[0].lower() != "none",
+                (
+                    {
+                        node
+                        for line_id in torsion_lines
+                        for node, end in zip(
+                            normalized_line_endpoints[line_id], ("a", "b"), strict=True
+                        )
+                        if torsion_ends.get(line_id, {}).get(end, False)
+                        and point_types[node] != "fixed"
+                    }
+                    if torsion_lines
+                    else None
+                ),
+                bool(finite_length_rods),
+                has_rigid6,
+                static_check,
+            )
+        except _AncillaryError as exc:
+            raise DeckFormatError(f"{self._label}: {exc}") from None
+
+    def _check_ancillary_files(
+        self,
+        option_by_key: dict[str, OptionRecord],
+        options: tuple[OptionRecord, ...],
+        motion_source: str | None,
+        mdc_wave_kin: int,
+        mdc_currents: int,
+        has_sea: bool,
+        has_current: bool,
+        roll_ends: set[str] | None,
+        has_rods: bool,
+        has_rigid6: bool,
+        static_check: bool,
+    ) -> None:
+        """Run the native rules of the ancillary files the deck names, those that are there.
+
+        A static check runs on the nominal clock the native parser installs before it validates.
+        """
+
+        def number(key: str) -> float | None:
+            row = option_by_key.get(key)
+            return None if row is None else _native_float(row.values[0])
+
+        dt = _STATIC_CHECK_DT if static_check else number("dtm")
+        tmax = number("tmax") or 0.0
+        depth = number("water_depth")
+        record_dir: float | None = None
+        component_headings: list[float] | None = None
+        waterkin_row = option_by_key.get("waterkin")
+        if waterkin_row is not None and not mdc_wave_kin:
+            value = waterkin_row.values[0]
+            if value.lower() not in {"", "none", "seastate"} and any(
+                character in _WATERKIN_FILENAME_LETTERS for character in value
+            ):
+                record_dir = self._check_waterkin_file(
+                    value, tmax, dt, depth, has_sea, has_current, mdc_currents, option_by_key
+                )
+        if mdc_wave_kin == 3:
+            component_headings = [0.0]
+        elif mdc_wave_kin == 7:
+            component_headings = self._wave_frequency_headings()
+        if mdc_currents:
+            name = _MOORDYN_C_CURRENT_FILES[mdc_currents]
+            records = self._ancillary_records(name)
+            if records is not None:
+                check = {
+                    1: _check_current_profile,
+                    2: _check_dynamic_current,
+                    5: _check_grid_current,
+                }[mdc_currents]
+                check(records, name)
+        rows = self._rows("POINTS")
+        for key in ("bathymetry",):
+            row = option_by_key.get(key)
+            if row is None:
+                continue
+            records = self._ancillary_records(row.values[0])
+            if records is not None:
+                _check_bathymetry(
+                    records,
+                    [
+                        (
+                            _native_int(point.tokens[0]),
+                            _native_float(point.tokens[2]),
+                            _native_float(point.tokens[3]),
+                            _native_float(point.tokens[4]),
+                        )
+                        for point in rows
+                        if point.tokens[1].lower() == "fixed"
+                    ],
+                )
+        self._check_syrope_files()
+        if motion_source is None or dt is None or self.caller_driven:
+            return
+        source = next(o for o in reversed(options) if o.keyword.lower() == motion_source)
+        records = self._ancillary_records(source.values[0])
+        if records is None:
+            return
+        if motion_source == "vesselmotion":
+            if has_rods or has_rigid6:
+                raise _AncillaryError(
+                    "vesselMotion/vesselRAO move Coupled/Vessel line points; decks with prescribed "
+                    "rods or Rigid6 bodies take motionFile rows"
+                )
+            _check_vessel_rows(records, dt, tmax, static_check)
+        elif motion_source == "vesselrao":
+            if has_rods or has_rigid6:
+                raise _AncillaryError(
+                    "vesselMotion/vesselRAO move Coupled/Vessel line points; decks with prescribed "
+                    "rods or Rigid6 bodies take motionFile rows"
+                )
+            headings = component_headings
+            if record_dir is not None:
+                headings = [record_dir]
+            if headings is None:
+                headings = self._sea_headings(option_by_key, options)
+            _check_rao_table(records, headings)
+        else:
+            self._check_motion_file(
+                records, dt, tmax, roll_ends, has_rods, has_rigid6, static_check
+            )
+
+    def _check_waterkin_file(
+        self,
+        value: str,
+        tmax: float,
+        dt: float | None,
+        depth: float | None,
+        has_sea: bool,
+        has_current: bool,
+        mdc_currents: int,
+        option_by_key: dict[str, OptionRecord],
+    ) -> float | None:
+        """Mirror ``read_waterkin_records``; returns the WaveDir of a WaveKinMod 1 record."""
+        records = self._ancillary_records(value)
+        mode = self._waterkin_wave_mode(value)
+        if records is None or mode is None:
+            return None
+
+        def line(index: int, name: str) -> str:
+            if index >= len(records):
+                raise _AncillaryError(f"WaterKin file has no {name} line")
+            return records[index]
+
+        tokens = line(4, "WaveKinFile").split()
+        wave_file = _strip_quotes(tokens[0]) if tokens else ""
+        tokens = line(5, "dtWave").split()
+        dt_wave = _token_real(tokens[0]) if tokens else None
+        if dt_wave is None or dt_wave < 0.0:
+            raise _AncillaryError("WaterKin dtWave must be finite and non-negative")
+        tokens = line(6, "WaveDir").split()
+        wave_dir = _token_real(tokens[0]) if tokens else None
+        if wave_dir is None:
+            raise _AncillaryError("WaterKin WaveDir must be finite")
+        if len(records) < 13:
+            raise _AncillaryError("WaterKin file ended inside the waves block")
+        if mode == 1:
+            if self.caller_driven:
+                raise _AncillaryError(
+                    "WaterKin WaveKinMod 1 is supported only by the standalone driver; coupled "
+                    "aggregate runs must use WaveKinMod 2/SEASTATE through external_fluid"
+                )
+            if not wave_file:
+                raise _AncillaryError("WaterKin WaveKinMod 1 requires WaveKinFile")
+            if dt_wave <= 0.0:
+                raise _AncillaryError("WaterKin WaveKinMod 1 requires positive dtWave")
+            if has_sea:
+                raise _AncillaryError(
+                    "the deck declares BOTH a waves OPTION and WaterKin WaveKinMod 1 "
+                    "(double-counting); keep one"
+                )
+            if dt is None or depth is None:
+                raise _AncillaryError("deck waves require dtM, TMax and WtrDpth")
+            history = self._ancillary_records(wave_file)
+            if history is None:
+                history = self._ancillary_records(
+                    str(self._ancillary_path(value).parent / wave_file)
+                )
+            if history is not None:
+                cutoff_row = option_by_key.get("wvhicoff")
+                cutoff = math.pi
+                if cutoff_row is not None:
+                    cutoff = (
+                        -1.0
+                        if cutoff_row.values[0].lower() == "none"
+                        else _native_float(cutoff_row.values[0])
+                    )
+                _check_elevation_record(history, dt_wave, tmax, cutoff, depth)
+        elif mode == 2 and not self.caller_driven:
+            raise _AncillaryError(
+                "WaterKin WaveKinMod 2/SEASTATE is coupled-only; the standalone driver has no host "
+                "SeaState field"
+            )
+        line(13, "current header")
+        tokens = line(14, "CurrentMod").split()
+        if not tokens or re.fullmatch(r"[+-]?[0-9]{1,18}", tokens[0]) is None:
+            raise _AncillaryError("WaterKin CurrentMod must be an integer (0, 1, or 2)")
+        current_mode = int(tokens[0])
+        if not 0 <= current_mode <= 2:
+            raise _AncillaryError(
+                "WaterKin CurrentMod supports 0 (none), 1 (depth profile), or 2 (SeaState hybrid)"
+            )
+        if current_mode == 2:
+            if mode == 1:
+                raise _AncillaryError(
+                    "WaterKin WaveKinMod 1 plus CurrentMod 2 mixes self-driven waves with a host "
+                    "SeaState current, which the nodewise host field cannot separate"
+                )
+            if not self.caller_driven:
+                raise _AncillaryError(
+                    "WaterKin WaveKinMod 2/SEASTATE is coupled-only; the standalone driver has no "
+                    "host SeaState field"
+                )
+        if current_mode == 1:
+            _check_waterkin_current(records[15:])
+            if has_current:
+                raise _AncillaryError(
+                    "the deck declares BOTH a current OPTION and a WaterKin current "
+                    "(double-counting); keep one"
+                )
+            if mdc_currents:
+                raise _AncillaryError(
+                    "the deck declares BOTH a current OPTION (or a WaterKin current) and "
+                    f"Currents {mdc_currents} (double-counting); keep one"
+                )
+        return wave_dir if mode == 1 else None
+
+    def _wave_frequency_headings(self) -> list[float] | None:
+        """The directions [deg] of the components of ``wave_frequencies.txt``, if it is there."""
+        records = self._ancillary_records(_MOORDYN_C_WAVE_FILES[7])
+        if records is None:
+            return None
+        headings = [
+            math.degrees(math.atan2(math.sin(beta), math.cos(beta)))
+            for omega, real_part, imaginary, beta in _wave_frequency_rows(records, "")
+            if omega > 0.0 and math.hypot(real_part, imaginary) > 0.0
+        ]
+        return headings
+
+    def _sea_headings(
+        self, option_by_key: dict[str, OptionRecord], options: tuple[OptionRecord, ...]
+    ) -> list[float]:
+        """The headings of a deck's own sea, as ``CD_Sea_Synthesise`` lays them out."""
+        bins_row = option_by_key.get("wave_directions")
+        bins = _SEA_DEFAULT_DIRECTIONS
+        if bins_row is not None and round(_native_float(bins_row.values[0])) > 0:
+            bins = round(_native_float(bins_row.values[0]))
+        spreading_row = option_by_key.get("wave_spreading")
+        spreading = 0.0 if spreading_row is None else _native_float(spreading_row.values[0])
+        headings: list[float] = []
+        waves = option_by_key.get("waves")
+        if waves is not None and waves.values[0].lower() != "none":
+            direction, _ = _train_direction(waves.values)
+            headings.extend(_spread_headings(direction, spreading, bins))
+        for option in options:
+            if option.keyword.lower() in _WAVETRAIN_KEYWORDS:
+                direction, own = _train_direction(option.values)
+                headings.extend(_spread_headings(direction, own, bins))
+        return headings
+
+    def _check_syrope_files(self) -> None:
+        """Mirror ``read_syrope_settings``, and ``CD_Syrope_Init`` for the types lines use."""
+        used = {row.tokens[1].lower() for row in self.sections}
+        for row in self.line_types:
+            parts = row.tokens[3].split("|")
+            if not parts[0].lower().startswith("syrope:"):
+                continue
+            records = self._ancillary_records(parts[0][len("syrope:") :])
+            if records is None:
+                continue
+            owc, shape, k1, k2 = _read_syrope_settings(records)
+            settings = self._ancillary_path(parts[0][len("syrope:") :])
+            name = _strip_quotes(owc)
+            absolute = name.startswith(("/", "\\")) or name[1:2] == ":"
+            table = self._ancillary_records(name if absolute else str(settings.parent / name))
+            if table is None:
+                continue
+            rows = _read_syrope_table(table)
+            if row.tokens[0].lower() in used:
+                _check_syrope_curve(
+                    rows, shape, k1, k2, _native_float(parts[1]), _native_float(parts[2])
+                )
+
+    def _check_motion_file(
+        self,
+        records: list[str],
+        dt: float,
+        tmax: float,
+        roll_ends: set[str] | None,
+        has_rods: bool,
+        has_rigid6: bool,
+        start_row_only: bool,
+    ) -> None:
+        """Mirror ``read_motion_file`` on the route the native dispatch takes.
+
+        Only the routes with lines read the file, and only those of lines alone its roll column.
+        """
+        if not self.lines and (has_rods or has_rigid6):
+            return
+        moved = {"coupled", "vessel"}
+        ids: dict[int, str] = {
+            _native_int(row.tokens[0]): row.tokens[1].lower() for row in self.points
+        }
+        turbine_rows = self._rows("TURBINES")
+        eligible: set[int] | None = {
+            point for point, kind in ids.items() if kind in moved or _is_turbine_point_type(kind)
+        }
+        if has_rods and has_rigid6:
+            eligible = None
+        elif has_rods:
+            assert eligible is not None
+            eligible |= {point for point, kind in ids.items() if _rod_point(kind) is not None}
+        elif has_rigid6:
+            assert eligible is not None
+            eligible |= {point for point, kind in ids.items() if kind.startswith("body")}
+        if eligible is not None and not eligible:
+            raise _AncillaryError("motionFile requires at least one Coupled/Vessel point")
+        turbines: dict[int, list[int]] | None = None
+        if turbine_rows:
+            turbines = {_native_int(row.tokens[0]): [] for row in turbine_rows}
+            for point, kind in ids.items():
+                if _is_turbine_point_type(kind):
+                    turbines.setdefault(_turbine_id(kind), []).append(point)
+        _check_motion_rows(
+            records,
+            dt,
+            tmax,
+            eligible,
+            None
+            if roll_ends is None or has_rods or has_rigid6
+            else {int(point) for point in roll_ends},
+            turbines,
+            start_row_only,
+        )
+
+    def _ancillary_path(self, value: str) -> Path:
+        """Return the file an ancillary value names, as the native ``path_join``."""
+        name = _strip_quotes(value)
+        absolute = name.startswith(("/", "\\")) or name[1:2] == ":"
+        return Path(name) if absolute else self.path.parent / name
+
+    def _ancillary_records(self, value: str) -> list[str] | None:
+        """Return the records of the ancillary file ``value`` names, ``None`` if it is not there.
+
+        A relative name is taken from the deck's folder.
+        """
+        source = self._ancillary_path(value)
+        try:
+            payload = source.read_bytes()
+        except OSError:
+            return None
+        records, _ = _split_records(payload.decode("utf-8", "surrogateescape"))
+        if records and records[0].startswith(_UTF8_BOM):
+            records[0] = records[0][len(_UTF8_BOM) :]
+        return records
+
+    def _waterkin_wave_mode(self, value: str) -> int | None:
+        """Return the ``WaveKinMod`` of WaterKin file ``value``, ``None`` if it is not there.
+
+        Mirrors the head of the native ``read_waterkin_records``: two header lines, a waves
+        header, then the WaveKinMod line.
+        """
+        records = self._ancillary_records(value)
+        if records is None:
+            return None
+        if len(records) < 3:
+            raise DeckFormatError(f"{self._label}: WaterKin file ended inside its header")
+        if len(records) < 4:
+            raise DeckFormatError(f"{self._label}: WaterKin file has no WaveKinMod line")
+        tokens = records[3].split()
+        if not tokens:
+            raise DeckFormatError(f"{self._label}: WaterKin WaveKinMod line has no value")
+        if _strip_quotes(tokens[0]).lower() == "seastate":
+            return 2
+        if re.fullmatch(r"[+-]?[0-9]{1,18}", tokens[0]) is None:
+            raise DeckFormatError(
+                f"{self._label}: WaterKin WaveKinMod must be 0, 1, 2, or SEASTATE"
+            )
+        mode = int(tokens[0])
+        if not 0 <= mode <= 2:
+            raise DeckFormatError(
+                f"{self._label}: WaterKin WaveKinMod supports 0 (none), 1 (elevation history), "
+                "or 2 (coupled SeaState hybrid)"
+            )
+        return mode
+
+    def _check_wave_frequencies(self, depth: float) -> None:
+        """Check ``wave_frequencies.txt`` in the deck's folder, if it is there.
+
+        Mirrors the ``WaveKin 7`` branch of the native ``read_moordyn_c_kinematics``, the
+        component bounds of its validation and the effective depth of the wave precompute.
+        A row takes four numbers, or three and a direction of zero.
+        """
+        records = self._ancillary_records(_MOORDYN_C_WAVE_FILES[7])
+        if records is None:
+            return
+        label = f"{self._label}: WaveKin 7: wave_frequencies.txt"
+        rows = _wave_frequency_rows(records, label)
+        if len(rows) < 2:
+            raise DeckFormatError(f"{label} needs at least two `omega Re Im` rows")
+        if rows[0][0] != 0.0:
+            raise DeckFormatError(f"{label}: the first row must be at omega = 0")
+        for index, (omega, real_part, imaginary, beta) in enumerate(rows):
+            if abs(beta) > 2.0 * math.pi:
+                raise DeckFormatError(f"{label} directions are radians within +-2 pi")
+            if index == 0:
+                continue
+            if omega <= 0.0:
+                raise DeckFormatError(f"{label} frequencies after the first row must be positive")
+            if math.hypot(real_part, imaginary) > 1.0e3 or omega > 1.0e4:
+                raise DeckFormatError(
+                    f"{label} components need amplitudes of at most 1e3 m and frequencies of "
+                    "at most 1e4 rad/s"
+                )
+        if depth + rows[0][2] <= 0.0:
+            raise DeckFormatError(f"{label}: the mean level leaves a non-positive effective depth")
 
     def _check_torsion_normal(
         self,
diff --git a/python/tests/test_deck_file.py b/python/tests/test_deck_file.py
index 543d24c..25fea08 100644
--- a/python/tests/test_deck_file.py
+++ b/python/tests/test_deck_file.py
@@ -834,7 +834,7 @@ def _vessel_source(rows: str) -> str:
 
 def test_vessel_motion_options_parse_and_rebase(tmp_path):
     source = _vessel_source("vessel.dat vesselMotion\n80.0|0.0|-10.0 vesselRef")
-    path = tmp_path / "vessel.dat"
+    path = tmp_path / "deck.dat"
     path.write_text(source, encoding="utf-8")
 
     deck = DeckFile.read(path)
@@ -2210,7 +2210,10 @@ def test_case_generation_rebases_relative_ancillary_paths(tmp_path):
 
     bathymetry = source_directory / "data" / "bathy_profile.dat"
     bathymetry.parent.mkdir()
-    bathymetry.write_text("placeholder\n", encoding="utf-8")
+    bathymetry.write_text(
+        "".join(f"{x} {y} 100.0\n" for x in (-2000.0, 2000.0) for y in (-2000.0, 2000.0)),
+        encoding="utf-8",
+    )
     deck_text = (ROOT / "examples/chain_catenary_r3_100m.dat").read_text(encoding="utf-8")
     deck_text = deck_text.replace("100.0        WtrDpth", "data/bathy_profile.dat bathymetryFile")
     deck_path = source_directory / "base.dat"
diff --git a/python/tests/test_deck_native_parity.py b/python/tests/test_deck_native_parity.py
index 4f2ad44..f968114 100644
--- a/python/tests/test_deck_native_parity.py
+++ b/python/tests/test_deck_native_parity.py
@@ -12,6 +12,7 @@
 import json
 import os
 import re
+import shutil
 import subprocess
 from pathlib import Path
 
@@ -86,6 +87,671 @@ def _rejects(text: str, message: str, *, caller_driven: bool = False) -> None:
 # --------------------------------------------------------------------------- options
 
 
+_WATERKIN = (
+    "MoorDyn v2 water kinematics file\n(header line 2)\n--- WAVES ---\n"
+    '{mode} WaveKinMod\n"eta.txt" WaveKinFile\n0.25 dtWave\n37.0 WaveDir\n'
+    "0 X\n0.0\n0 Y\n0.0\n0 Z\n0.0\n--- CURRENT ---\n0 CurrentMod\n0 Z\n0.0\n"
+    "--- need this line ---\n"
+)
+_ETA = "time elevation\n" + "".join(
+    f"{0.5 * i:g} {0.4 * (-1.0) ** (i // 4):g}\n" for i in range(41)
+)
+_RAO = (EXAMPLES / "data" / "vessel" / "sample_rao.txt").read_text(encoding="utf-8")
+_ROWS = "0 0 0 0\n0.5 0.5 0 0.3\n0.6 0.3 0.1 2.4\n0.6 0.2 -0.1 -1.0\n"
+
+
+def _waterkin(mode: str) -> dict[str, str]:
+    return {"wk.dat": _WATERKIN.format(mode=mode), "eta.txt": _ETA}
+
+
+# (rows of the deck, files next to it): every case is judged by DeckFile and by the driver
+_WAVE_RECORD_CASES = [
+    # a reference point needs component waves
+    (("12.5 WaveRefX",), {}),
+    (("0 WaveRefY",), {}),
+    (("airy 2 8 0 waves", "1 WaveRefX"), {}),
+    (("airy 2 8 0 wavetrain", "1 WaveRefY"), {}),
+    (("7 WaterKin", "12.5 WaveRefX", "-3 WaveRefY"), {"wave_frequencies.txt": _ROWS}),
+    (("7 WaterKin", "1.0e6 WaveRefX"), {"wave_frequencies.txt": _ROWS}),
+    (("7 WaterKin", "1.0000001e6 WaveRefX"), {"wave_frequencies.txt": _ROWS}),
+    (("7 WaterKin", "-2e6 WaveRefY"), {"wave_frequencies.txt": _ROWS}),
+    (("7 WaterKin", "nan WaveRefY"), {"wave_frequencies.txt": _ROWS}),
+    (("7 WaterKin", "north WaveRefY"), {"wave_frequencies.txt": _ROWS}),
+    (("7 WaterKin", "1e-320 WaveRefX"), {"wave_frequencies.txt": _ROWS}),
+    (("7 WaterKin", "2e6 WaveRefX", "1 WaveRefX"), {"wave_frequencies.txt": _ROWS}),
+    (("7 WaterKin", "1 WAVEREFX"), {"wave_frequencies.txt": _ROWS}),
+    # the WaveKinMod of a WaterKin file decides
+    (("wk.dat WaterKin", "12.5 WaveRefX"), _waterkin("1")),
+    (("wk.dat WaterKin", "12.5 WaveRefX"), _waterkin("0")),
+    (("wk.dat WaterKin", "0 WaveRefY"), _waterkin("0")),
+    (("wk.dat WaterKin", "12.5 WaveRefX"), _waterkin("2")),
+    (("wk.dat WaterKin", "12.5 WaveRefX"), _waterkin("SEASTATE")),
+    (("wk.dat WaterKin", "12.5 WaveRefX"), _waterkin("3")),
+    (("wk.dat WaterKin", "12.5 WaveRefX"), _waterkin("one")),
+    (("wk.dat WaterKin", "12.5 WaveRefX"), {"wk.dat": "header\nheader\n"}),
+    (("wk.dat WaterKin", "12.5 WaveRefX"), {"wk.dat": "header\nheader\nwaves\n"}),
+    (("wk.dat WaterKin", "12.5 WaveRefX"), {"wk.dat": "header\nheader\nwaves\n\n"}),
+    (("wk.dat WaterKin", "rao.txt vesselRAO"), {**_waterkin("0"), "rao.txt": _RAO}),
+    (("wk.dat WaterKin", "rao.txt vesselRAO"), {**_waterkin("1"), "rao.txt": _RAO}),
+    # the rows of wave_frequencies.txt
+    (("7 WaterKin",), {"wave_frequencies.txt": _ROWS}),
+    (("7 WaterKin",), {"wave_frequencies.txt": "0 0 0 1.0\n0.5 0.5 0 0.3\n\n0.6 0.3 0.1 0.3\n"}),
+    (("7 WaterKin",), {"wave_frequencies.txt": "0 0 0\n0.5 0.5 0\n0.6 0.3 0.1 -6.2\n"}),
+    (("7 WaterKin",), {"wave_frequencies.txt": "0 0 0 0\n0.5 0.5 0 6.3\n"}),
+    (("7 WaterKin",), {"wave_frequencies.txt": "0 0 0 -6.3\n0.5 0.5 0 0\n"}),
+    (("7 WaterKin",), {"wave_frequencies.txt": "0 0 0 0\n0.5 0.5 0 east\n"}),
+    (("7 WaterKin",), {"wave_frequencies.txt": "0 0 0 0\n0.5 0.5 0 nan\n"}),
+    (("7 WaterKin",), {"wave_frequencies.txt": "0 0 0 0\n0.5 nan 0 0\n"}),
+    (("7 WaterKin",), {"wave_frequencies.txt": "0 0 0 0\n0.5 0.5\n"}),
+    (("7 WaterKin",), {"wave_frequencies.txt": "omega Re Im\n0 0 0 0\n0.5 0.5 0 0\n"}),
+    (("7 WaterKin",), {"wave_frequencies.txt": "0 0 0 0\n"}),
+    (("7 WaterKin",), {"wave_frequencies.txt": "0.1 0 0 0\n0.5 0.5 0 0\n"}),
+    (("7 WaterKin",), {"wave_frequencies.txt": "0 0 0 0\n0.5 0.5 0 0\n0 0.5 0 0\n"}),
+    (("7 WaterKin",), {"wave_frequencies.txt": "0 0 0 0\n0.5 2000 0 0\n"}),
+    (("7 WaterKin",), {"wave_frequencies.txt": "0 0 0 0\n2e4 0.1 0 0\n"}),
+    (("7 WaterKin",), {"wave_frequencies.txt": "0 0 -60 0\n0.5 0.5 0 0\n"}),
+]
+
+
+_CURRENT_TABLES = {
+    "2": "a\nb\nc\nd\n-40 -10\ne\n0 0.1 0.2\n5 0.2 0.3\n10 0.2 0.3\n",
+    "5": "1 1 2 2\n0\n0\n-40 -10\n0 10\n"
+    + "".join(f"0 0 {z} {t} 0.1 0 0\n" for z in (-40, -10) for t in (0, 10)),
+}
+_PROFILE = "z ux uy uz\n-40 0.1 0 0\n-10 0.2 0.1 0\n"
+_TABLE = "-40 0.1 0.0\n-10 0.2 0.1\n"
+
+
+def _kin(wave="0", current="0", table="", dt="0.25", direction="37.0", record='"eta.txt"') -> str:
+    return (
+        "MoorDyn v2 water kinematics file\n(header line 2)\n--- WAVES ---\n"
+        f"{wave} WaveKinMod\n{record} WaveKinFile\n{dt} dtWave\n{direction} WaveDir\n"
+        "0 X\n0.0\n0 Y\n0.0\n0 Z\n0.0\n--- CURRENT ---\n"
+        f"{current} CurrentMod\nz ux uy\n(m) (m/s) (m/s)\n{table}--- need this line ---\n"
+    )
+
+
+def _hold(point=2, *, tail=(), steps=101, last=None) -> str:
+    rows = [f"{0.1 * i:.1f} {point} 0 0 0 0 0 0 0 0 0\n" for i in range(steps)]
+    return "".join(rows[:last]) + "".join(tail)
+
+
+def _vessel(*, form=19, steps=101, tail=()) -> str:
+    pose = "0 0 0 0 0 0" if form == 19 else "0 0 0 1 0 0 0"
+    rows = [f"{0.1 * i:.1f} {pose} {'0 ' * 12}\n" for i in range(steps)]
+    return "".join(rows) + "".join(tail)
+
+
+def _rao(
+    heads=(-180.0, 179.0), period_row="0.5 0 0 0 0 1 0 0 0 0 0 0 0\n100 0 0 0 0 1 0 0 0 0 0 0 0\n"
+):
+    return "".join(f"HEADING {head}\n{period_row}" for head in heads)
+
+
+def _sea_bed(depth=50.0, **changes) -> str:
+    rows = {(x, y): depth for x in (-1000.0, 1000.0) for y in (-1000.0, 1000.0)}
+    rows.update(changes.get("rows", {}))
+    return "".join(f"{x} {y} {value}\n" for (x, y), value in rows.items()) + changes.get("tail", "")
+
+
+_SETTINGS = (EXAMPLES / "data" / "syrope" / "syrope_settings.dat").read_text(encoding="utf-8")
+_OWC = (EXAMPLES / "data" / "syrope" / "syrope_owc.dat").read_text(encoding="utf-8")
+
+
+def _syrope(settings=_SETTINGS, owc=_OWC) -> dict[str, str]:
+    return {"data/syrope/syrope_settings.dat": settings, "data/syrope/syrope_owc.dat": owc}
+
+
+# a deck of objects and no lines: the only route that does not read a motion file
+_OBJECTS = """Dry rod pinned to a heavy body
+--- ROD TYPES ---
+Name Diam Mass Cd Ca CdEnd CaEnd
+(-) (m) (kg/m) (-) (-) (-) (-)
+pend 0.2 50 0 0 0 0
+--- BODIES ---
+ID Attachment X0 Y0 Z0 r0 p0 y0 Mass CG* I* Volume CdA* Ca*
+(#) (-) (m) (m) (m) (deg) (deg) (deg) (kg) (m) (kg-m^2) (m^3) (m^2) (-)
+1 Free 5 0 20 0 0 0 5.0e10 0 1.0e13 4.87804882926829E+07 0 0
+--- RODS ---
+ID RodType Attachment XA YA ZA XB YB ZB NumSegs Outputs
+(-) (-) (-) (m) (m) (m) (m) (m) (m) (-) (-)
+1 pend Body1Pinned 0 0 0 5.0 0 -8.66025403784439 4 -
+--- OPTIONS ---
+100 WtrDpth
+deck bodyIC
+moordyn bodyWetting
+0.05 dtM
+1.0 TMax
+--- OUTPUTS ---
+Rod1Px Body1Pz
+--- end ---
+"""
+_FREE = "wave_frequencies.txt"
+_ONE_BETA = "0 0 0 0.5\n0.5 0.5 0 0.5\n"
+_WK = ("wk.dat WaterKin",)
+_C1 = "current_profile.txt"
+_C2 = "current_profile_dynamic.txt"
+_C5 = "current_profile_4d.txt"
+_RAO_ROWS = ("rao.txt vesselRAO", "0|0|0 vesselRef")
+# (deck, rows, files): every case is judged by DeckFile and by the driver
+_ANCILLARY_CASES = [
+    # WaterKin: the rest of the file
+    ("held", _WK, {"wk.dat": _kin()}),
+    ("held", _WK, {"wk.dat": _kin(dt="-1")}),
+    ("held", _WK, {"wk.dat": _kin(dt="x")}),
+    ("held", _WK, {"wk.dat": _kin(direction="north")}),
+    ("held", _WK, {"wk.dat": 'a\nb\nc\n0\n""\n0\n0\n0\n'}),
+    ("held", _WK, {"wk.dat": _kin(current="x")}),
+    ("held", _WK, {"wk.dat": _kin(current="3")}),
+    ("held", _WK, {"wk.dat": _kin(current="2")}),
+    ("held", _WK, {"wk.dat": _kin(wave="1", current="2"), "eta.txt": _ETA}),
+    ("held", _WK, {"wk.dat": _kin(wave="1", record='""')}),
+    ("held", _WK, {"wk.dat": _kin(wave="1", dt="0"), "eta.txt": _ETA}),
+    ("held", ("wk.dat WaterKin", "airy 2 8 0 waves"), {"wk.dat": _kin("1"), "eta.txt": _ETA}),
+    # the WaterKin current table, and a current doubled
+    ("held", _WK, {"wk.dat": _kin(current="1", table=_TABLE)}),
+    ("held", _WK, {"wk.dat": _kin(current="1", table="10 0.2 0.1\n40 0.1 0.0\n")}),
+    (
+        "held",
+        ("wk.dat WaterKin", "uniform 0.1 0 0 current"),
+        {"wk.dat": _kin(current="1", table=_TABLE)},
+    ),
+    ("held", ("wk.dat WaterKin", "uniform 0.1 0 0 current"), {"wk.dat": _kin()}),
+    (
+        "held",
+        ("wk.dat WaterKin", "1 Currents"),
+        {"wk.dat": _kin(current="1", table=_TABLE), _C1: _PROFILE},
+    ),
+    ("held", ("wk.dat WaterKin", "1 Currents"), {"wk.dat": _kin(), _C1: _PROFILE}),
+    ("held", _WK, {"wk.dat": _kin(current="1", table="-40 0.1 0\n")}),
+    ("held", _WK, {"wk.dat": _kin(current="1", table="-40 0.1 0\n10 0.2 0\n")}),
+    ("held", _WK, {"wk.dat": _kin(current="1", table="-40 0.1 0\n-40 0.2 0\n")}),
+    ("held", _WK, {"wk.dat": _kin(current="1", table="-40 0.1\n-10 0.2 0\n")}),
+    ("held", _WK, {"wk.dat": _kin(current="1", table="-40 2000 0\n-10 0.2 0\n")}),
+    # the record of WaveKinMod 1
+    ("held", _WK, {"wk.dat": _kin("1"), "eta.txt": "0 0\n1 0\n2 0\n"}),
+    ("held", _WK, {"wk.dat": _kin("1"), "eta.txt": "1 0\n2 0\n3 0\n4 0\n"}),
+    ("held", _WK, {"wk.dat": _kin("1"), "eta.txt": "0 0\n2 0\n1 0\n3 0\n"}),
+    ("held", _WK, {"wk.dat": _kin("1"), "eta.txt": "0 0\n1 0\nend\n2 0\n3 0\n"}),
+    ("held", _WK, {"wk.dat": _kin("1"), "eta.txt": "0 0\n1 x\n2 0\n3 0\n"}),
+    ("held", _WK, {"wk.dat": _kin("1"), "eta.txt": _ETA.replace(" 0.4\n", " 4000\n")}),
+    (
+        "held",
+        _WK,
+        {
+            "wk.dat": _kin("1"),
+            "eta.txt": _ETA.replace(" 0.4\n", " -70\n").replace(" -0.4\n", " -70\n"),
+        },
+    ),
+    ("held", ("wk.dat WaterKin", "0.1 WvHiCOff"), {"wk.dat": _kin("1"), "eta.txt": _ETA}),
+    ("held", ("wk.dat WaterKin", "none WvHiCOff"), {"wk.dat": _kin("1"), "eta.txt": _ETA}),
+    # Currents 1, 2 and 5
+    ("held", ("1 Currents",), {_C1: _PROFILE}),
+    ("held", ("1 Currents",), {_C1: "z ux uy uz\n"}),
+    ("held", ("1 Currents",), {_C1: "-10 0 0 0\n-40 0 0 0\n"}),
+    ("held", ("1 Currents",), {_C1: "-40 2000 0 0\n-10 0 0 0\n"}),
+    ("held", ("1 Currents",), {_C1: "-40 0 0\n-10 0 0 x\n"}),
+    ("held", ("1 Currents",), {_C1: "-40 0\n"}),
+    ("held", ("2 Currents",), {_C2: _CURRENT_TABLES["2"]}),
+    ("held", ("2 Currents",), {_C2: "a\nb\n"}),
+    ("held", ("2 Currents",), {_C2: "a\nb\nc\nd\n-40 -10\ne\n0 0.1 0.2\n"}),
+    ("held", ("2 Currents",), {_C2: "a\nb\nc\nd\n-10 -40\ne\n0 1 1\n5 1 1\n"}),
+    ("held", ("2 Currents",), {_C2: "a\nb\nc\nd\n-40 -10\ne\n0 1 1\n0 1 1\n"}),
+    ("held", ("2 Currents",), {_C2: "a\nb\nc\nd\n-40 -10\ne\n0 1\n5 1 1\n"}),
+    ("held", ("2 Currents",), {_C2: "a\nb\nc\nd\n-40 -10\ne\n0 x 1\n5 1 1\n"}),
+    ("held", ("2 Currents",), {_C2: "a\nb\nc\nd\n-40 -10\ne\n0 2000 1\n5 1 1\n"}),
+    ("held", ("5 Currents",), {_C5: _CURRENT_TABLES["5"]}),
+    ("held", ("5 Currents",), {_C5: "1 1 x 2\n"}),
+    ("held", ("5 Currents",), {_C5: "1 1 2 2\n0\n0\n-40\n0 10\n"}),
+    ("held", ("5 Currents",), {_C5: _CURRENT_TABLES["5"] + "0 0 -40 0 0.1 0 0\n"}),
+    (
+        "held",
+        ("5 Currents",),
+        {_C5: _CURRENT_TABLES["5"].replace("0 0 -40 0 0.1", "0 0 -41 0 0.1")},
+    ),
+    ("held", ("5 Currents",), {_C5: _CURRENT_TABLES["5"].replace("0 0 -40 0 0.1 0 0\n", "")}),
+    ("held", ("5 Currents",), {_C5: "1 1 1 1\n0\n0\n-40\n0\n0 0 -40 0 0.1 0 0\n"}),
+    (
+        "held",
+        ("5 Currents",),
+        {_C5: _CURRENT_TABLES["5"].replace("0 0 -40 0 0.1", "0 0 -40 0 2000")},
+    ),
+    # motionFile
+    ("held", ("m.txt motionFile",), {"m.txt": _hold()}),
+    ("held", ("m.txt motionFile",), {"m.txt": _hold(tail=("10.1 2 0 0 0 0 0 0 0 0 0\n",))}),
+    ("held", ("m.txt motionFile",), {"m.txt": _hold(last=100)}),
+    ("held", ("m.txt motionFile",), {"m.txt": _hold(tail=("0.0 2 0 0 0 0 0 0 0 0 0\n",))}),
+    ("held", ("m.txt motionFile",), {"m.txt": _hold(tail=("0.05 2 0 0 0 0 0 0 0 0 0\n",))}),
+    ("held", ("m.txt motionFile",), {"m.txt": _hold(tail=("0.0 1 0 0 0 0 0 0 0 0 0\n",))}),
+    ("held", ("m.txt motionFile",), {"m.txt": _hold(tail=("0.0 2.5 0 0 0 0 0 0 0 0 0\n",))}),
+    ("held", ("m.txt motionFile",), {"m.txt": _hold(tail=("0.0 2 0 0 0\n",))}),
+    ("held", ("m.txt motionFile",), {"m.txt": _hold(point=1)}),
+    ("held", ("m.txt motionFile",), {"m.txt": "0.0 2 2e6 0 0 0 0 0 0 0 0\n"}),
+    ("held", ("m.txt motionFile",), {"m.txt": "0.0 2 0 0 0 2e4 0 0 0 0 0\n"}),
+    ("held", ("m.txt motionFile",), {"m.txt": "0.0 2 0 0 0 0 0 0 2e6 0 0\n"}),
+    # vesselMotion
+    ("held", ("v.txt vesselMotion", "0|0|0 vesselRef"), {"v.txt": _vessel()}),
+    ("held", ("v.txt vesselMotion", "0|0|0 vesselRef"), {"v.txt": _vessel(form=20)}),
+    ("held", ("v.txt vesselMotion", "0|0|0 vesselRef"), {"v.txt": _vessel(steps=100)}),
+    (
+        "held",
+        ("v.txt vesselMotion", "0|0|0 vesselRef"),
+        {"v.txt": _vessel(form=19, tail=("0.0 " + "0 " * 19 + "\n",))},
+    ),
+    (
+        "held",
+        ("v.txt vesselMotion", "0|0|0 vesselRef"),
+        {"v.txt": _vessel(form=20).replace("0 0 0 1 0 0 0", "0 0 0 2 0 0 0", 1)},
+    ),
+    ("held", ("v.txt vesselMotion", "0|0|0 vesselRef"), {"v.txt": "0.05 " + "0 " * 18 + "\n"}),
+    ("held", ("v.txt vesselMotion", "0|0|0 vesselRef"), {"v.txt": "0.0 2e6 " + "0 " * 17 + "\n"}),
+    # vesselRAO, its rows and the headings of the sea
+    ("held", ("airy 2 8 30 waves", *_RAO_ROWS), {"rao.txt": _rao()}),
+    ("held", ("airy 2 8 30 waves", *_RAO_ROWS), {"rao.txt": _rao(heads=(30.0,))}),
+    ("held", ("airy 2 8 30 waves", *_RAO_ROWS), {"rao.txt": _rao(heads=(0.0,))}),
+    ("held", ("airy 2 8 30 waves", *_RAO_ROWS), {"rao.txt": _rao(heads=(40.0, 90.0))}),
+    ("held", ("airy 2 8 30 waves", *_RAO_ROWS), {"rao.txt": _rao(heads=(0.0, 90.0))}),
+    ("held", ("airy 2 8 30 waves", *_RAO_ROWS), {"rao.txt": _rao(heads=(0.0, 360.0))}),
+    ("held", ("airy 2 8 30 waves", *_RAO_ROWS), {"rao.txt": "1 0 0 0 0 1 0 0 0 0 0 0 0\n"}),
+    ("held", ("airy 2 8 30 waves", *_RAO_ROWS), {"rao.txt": "HEADING 30\n"}),
+    ("held", ("airy 2 8 30 waves", *_RAO_ROWS), {"rao.txt": "HEADING 30 40\n"}),
+    ("held", ("airy 2 8 30 waves", *_RAO_ROWS), {"rao.txt": "HEADING 30\n1 0 0\n"}),
+    ("held", ("airy 2 8 30 waves", *_RAO_ROWS), {"rao.txt": "HEADING 30\n-1 " + "0 " * 12 + "\n"}),
+    (
+        "held",
+        ("airy 2 8 30 waves", *_RAO_ROWS),
+        {"rao.txt": "HEADING 30\n1 -1 " + "0 " * 11 + "\n"},
+    ),
+    (
+        "held",
+        ("airy 2 8 30 waves", *_RAO_ROWS),
+        {"rao.txt": "HEADING 30\n1 " + "0 " * 12 + "\nHEADING 40\n2 " + "0 " * 12 + "\n"},
+    ),
+    ("held", ("jonswap 2 9 3.3 30 waves", *_RAO_ROWS), {"rao.txt": _rao(heads=(30.0,))}),
+    (
+        "held",
+        ("jonswap 2 9 3.3 30 waves", "4 WaveSpreading", *_RAO_ROWS),
+        {"rao.txt": _rao(heads=(30.0,))},
+    ),
+    (
+        "held",
+        ("jonswap 2 9 3.3 30 waves", "4 WaveSpreading", *_RAO_ROWS),
+        {"rao.txt": _rao(heads=(-60.0, 120.0))},
+    ),
+    (
+        "held",
+        ("jonswap 2 9 3.3 30 waves", "4 WaveSpreading", "5 WaveDirections", *_RAO_ROWS),
+        {"rao.txt": _rao(heads=(-42.0, 102.0))},
+    ),
+    # a narrow spreading underflows the outer bins: only the bins with energy need a heading
+    (
+        "held",
+        ("jonswap 2 9 3.3 0 waves", "1000 WaveSpreading", *_RAO_ROWS),
+        {"rao.txt": _rao(heads=(-40.0, 40.0))},
+    ),
+    (
+        "held",
+        ("jonswap 2 9 3.3 0 waves", "1000 WaveSpreading", *_RAO_ROWS),
+        {"rao.txt": _rao(heads=(-20.0, 20.0))},
+    ),
+    (
+        "held",
+        ("jonswap 2 9 3.3 0 waves", "838 WaveSpreading", *_RAO_ROWS),
+        {"rao.txt": _rao(heads=(-40.0, 40.0))},
+    ),
+    (
+        "held",
+        ("jonswap 2 9 3.3 0 waves", "700 WaveSpreading", *_RAO_ROWS),
+        {"rao.txt": _rao(heads=(-40.0, 40.0))},
+    ),
+    (
+        "held",
+        ("jonswap 2 9 3.3 0 1000 wavetrain", *_RAO_ROWS),
+        {"rao.txt": _rao(heads=(-40.0, 40.0))},
+    ),
+    (
+        "held",
+        ("jonswap 2 9 3.3 0 1000 wavetrain", "5 WaveDirections", *_RAO_ROWS),
+        {"rao.txt": _rao(heads=(-36.0, 36.0))},
+    ),
+    ("held", ("7 WaterKin", *_RAO_ROWS), {"rao.txt": _rao(), _FREE: _ROWS}),
+    ("held", ("7 WaterKin", *_RAO_ROWS), {"rao.txt": _rao(heads=(0.0, 100.0)), _FREE: _ROWS}),
+    ("held", ("7 WaterKin", *_RAO_ROWS), {"rao.txt": _rao(heads=(28.0, 29.0)), _FREE: _ONE_BETA}),
+    (
+        "held",
+        ("wk.dat WaterKin", *_RAO_ROWS),
+        {"rao.txt": _rao(heads=(37.0,)), "wk.dat": _kin("1"), "eta.txt": _ETA},
+    ),
+    (
+        "held",
+        ("wk.dat WaterKin", *_RAO_ROWS),
+        {"rao.txt": _rao(heads=(38.0,)), "wk.dat": _kin("1"), "eta.txt": _ETA},
+    ),
+    # bathymetry
+    ("bed", (), {"b.txt": _sea_bed()}),
+    ("bed", (), {"b.txt": _sea_bed(depth=49.0)}),
+    ("bed", (), {"b.txt": _sea_bed(depth=50.0004)}),
+    ("bed", (), {"b.txt": _sea_bed(depth=0.0)}),
+    ("bed", (), {"b.txt": "0 0 50\n1 0 50\n"}),
+    ("bed", (), {"b.txt": _sea_bed(tail="0 0 50\n")}),
+    ("bed", (), {"b.txt": _sea_bed(tail="-1000.0 -1000.0 50\n")}),
+    ("bed", (), {"b.txt": _sea_bed(tail="1 2\n")}),
+    ("bed", (), {"b.txt": _sea_bed(tail="1 x 3\n")}),
+    ("bed", (), {"b.txt": _sea_bed(tail="-1000.0 0.0 50\n1000.0 0.0 50\n")}),
+    # a crease at the anchor: the rounded and the tabulated surfaces differ there
+    (
+        "bed",
+        (),
+        {
+            "b.txt": "".join(
+                f"{x} {y} {d}\n"
+                for x, d in ((-1000.0, 40.0), (400.0, 50.0), (1000.0, 40.0))
+                for y in (-1000.0, 1000.0)
+            )
+        },
+    ),
+    # Syrope settings and OWC table
+    ("syrope", (), _syrope()),
+    ("syrope", (), _syrope(settings=_SETTINGS + "1.0 k1\n")),
+    ("syrope", (), _syrope(settings="data/x.dat OWC\n")),
+    (
+        "syrope",
+        (),
+        _syrope(
+            settings=_SETTINGS.replace("LINEAR", "CUBIC")
+            .replace("QUADRATIC", "CUBIC")
+            .replace("EXP", "CUBIC")
+        ),
+    ),
+    ("syrope", (), _syrope(owc="strain tension\n0.0 0.0\n")),
+    ("syrope", (), _syrope(owc="0.0 0.0\n0.01 x\n")),
+    ("syrope", (), _syrope(owc="0.02 1e6\n0.01 2e6\n0.03 3e6\n")),
+    ("syrope", (), _syrope(owc="0.0 2e6\n0.01 1e6\n0.03 3e6\n")),
+    # the routes of the objects: which of them read a motion file
+    ("rigid6_buoy.dat", ("m.txt motionFile",), {"m.txt": "not a motion row\n"}),
+    ("rod_moored_spar.dat", ("m.txt motionFile",), {"m.txt": "not a motion row\n"}),
+    ("spar_pinned_rods.dat", ("m.txt motionFile",), {"m.txt": "not a motion row\n"}),
+    ("mixed_body_rods_points.dat", ("m.txt motionFile",), {"m.txt": "not a motion row\n"}),
+    ("objects", ("m.txt motionFile",), {"m.txt": "not a motion row\n"}),
+]
+# every spelling the reader takes of the options and kinds these rules read
+_KEYWORD_SPELLINGS = (
+    ("waves", "wave"),
+    ("wavetrain", "wave_train"),
+    ("WaveSpreading", "wave_spreading"),
+    ("WaveDirections", "wave_directions"),
+    ("WaveComponents", "wave_components"),
+    ("WtrDpth", "water_depth"),
+    ("dtM", "dt"),
+    ("WaterKin", "WaveKin"),
+    ("bathymetryFile", "seafloor_file"),
+    ("rampTime", "tramp"),
+)
+_KIND_VALUES = {
+    "airy": ("2 8 30", "2 8 30"),
+    "Airy": ("2 8 30", "2 8 30"),
+    "jonswap": ("2 9 3.3 30", "2 9 3.3 30 0"),
+    "JONSWAP": ("2 9 3.3 30", "2 9 3.3 30 0"),
+    "pm": ("2 9 30", "2 9 30 0"),
+    "issc": ("2 9 30", "2 9 30 0"),
+    "bretschneider": ("2 9 30", "2 9 30 0"),
+    "torsethaugen": ("6 12 30", "6 12 30 0"),
+    "ochihubble": ("2 8 1 1 12 1 30", "2 8 1 1 12 1 30 0"),
+    "ochi-hubble": ("2 8 1 1 12 1 30", "2 8 1 1 12 1 30 0"),
+    "ochi_hubble": ("2 8 1 1 12 1 30", "2 8 1 1 12 1 30 0"),
+    "Ochi_Hubble": ("2 8 1 1 12 1 30", "2 8 1 1 12 1 30 0"),
+}
+# every kind and keyword against a one-heading RAO table: at the sea's heading and off it, on
+# the dynamic clock and on the static check's
+_SPELLED_CASES = [
+    (base, (f"{kind} {values[train]} {keyword}", *_RAO_ROWS), {"rao.txt": _rao(heads=(heading,))})
+    for kind, values in _KIND_VALUES.items()
+    for train, keyword in ((0, "waves"), (0, "wave"), (1, "wavetrain"), (1, "wave_train"))
+    for heading in (30.0, 31.0)
+    for base in ("held", "held_static")
+]
+_SPELLED_CASES += [
+    # the kinds whose rows spread, and the spellings of the options that spread them
+    (base, (f"pm 2 9 30 {keyword}", f"4 {spread}", *_RAO_ROWS), {"rao.txt": _rao(heads=heads)})
+    for base in ("held", "held_static")
+    for keyword in ("waves", "wave")
+    for spread in ("WaveSpreading", "wave_spreading", "WAVESPREADING")
+    for heads in ((30.0,), (-180.0, 179.0))
+] + [
+    (base, (f"ochi_hubble 2 8 1 1 12 1 30 4 {keyword}", *_RAO_ROWS), {"rao.txt": _rao(heads=heads)})
+    for base in ("held", "held_static")
+    for keyword in ("wavetrain", "wave_train")
+    for heads in ((30.0,), (-180.0, 179.0))
+]
+# a record or a kinematics file on the static check's clock, under each spelling of its selector
+_STATIC_CASES = [
+    ("held_static", (f"{value} {keyword}", *rows), files)
+    for keyword in ("WaterKin", "WaveKin", "WATERKIN")
+    for value, rows, files in (
+        ("wk.dat", (), {"wk.dat": _kin("1"), "eta.txt": _ETA}),
+        ("wk.dat", ("12.5 WaveRefX",), {"wk.dat": _kin("1"), "eta.txt": _ETA}),
+        ("wk.dat", ("12.5 WaveRefX",), {"wk.dat": _kin("0"), "eta.txt": _ETA}),
+        ("wk.dat", (), {"wk.dat": _kin("1"), "eta.txt": "0 0\n1 0\n"}),
+        ("wk.dat", (), {"wk.dat": _kin("1", dt="100"), "eta.txt": _ETA}),
+        (
+            "wk.dat",
+            _RAO_ROWS,
+            {"wk.dat": _kin("1"), "eta.txt": _ETA, "rao.txt": _rao(heads=(37.0,))},
+        ),
+        (
+            "wk.dat",
+            _RAO_ROWS,
+            {"wk.dat": _kin("1"), "eta.txt": _ETA, "rao.txt": _rao(heads=(38.0,))},
+        ),
+        ("7", (), {_FREE: _ROWS}),
+        ("7", ("12.5 WaveRefY",), {_FREE: _ROWS}),
+        ("7", (), {_FREE: "0 0 0 0\n0.5 0.5 0 6.3\n"}),
+        ("3", ("0.1 dtWave",), {"wave_elevation.txt": _ETA.replace("time elevation\n", "")}),
+    )
+] + [
+    ("held_static", ("m.txt motionFile",), {"m.txt": "0.0 2 0 0 0 0 0 0 0 0 0\n"}),
+    (
+        "held_static",
+        ("m.txt motionFile",),
+        {"m.txt": "0.0 2 0 0 0 0 0 0 0 0 0\n0.05 2 0 0 0 0 0 0 0 0 0\n"},
+    ),
+    ("held_static", ("m.txt motionFile",), {"m.txt": _hold()}),
+    ("held_static", ("m.txt motionFile",), {"m.txt": "1.0 2 0 0 0 0 0 0 0 0 0\n"}),
+    ("held_static", ("v.txt vesselMotion", "0|0|0 vesselRef"), {"v.txt": _vessel(steps=1)}),
+    ("held_static", ("v.txt vesselMotion", "0|0|0 vesselRef"), {"v.txt": _vessel()}),
+    ("bed_static", (), {"b.txt": _sea_bed()}),
+    ("held", ("airy 2 8 0 waves", "4 WaveSpreading"), {}),
+]
+
+
+def _shipped_decks() -> list[str]:
+    decks = [path.name for path in sorted(EXAMPLES.glob("*.dat"))]
+    decks.remove("cabledyn_options_reference.dat")
+    return [*decks, "moordynC_wavekin/chain_wavekin7_currents1.dat"]
+
+
+def _carries_waves(name: str) -> bool:
+    text = (EXAMPLES / name).read_text(encoding="utf-8")
+    return re.search(r"(?mi)^(?!none\b|--)\S.*\s(?:waves?|wave_?train)\b", text) is not None or (
+        re.search(r"(?mi)^(?!0\s|none\s)\S+\s+(?:WaterKin|WaveKin)\b", text) is not None
+    )
+
+
+# every shipped deck as a static check without dtM; every shipped wave deck under each spelling
+_SHIPPED_CASES = [(f"shipped:static:{name}", (), {}) for name in _shipped_decks()] + [
+    (f"shipped:{variant}:{name}", (), {})
+    for name in _shipped_decks()
+    if _carries_waves(name)
+    for variant in ("spelled", "upper")
+]
+_ALL_FILE_CASES = [
+    *[("held", rows, files) for rows, files in _WAVE_RECORD_CASES],
+    *_ANCILLARY_CASES,
+    *_SPELLED_CASES,
+    *_STATIC_CASES,
+    *_SHIPPED_CASES,
+]
+
+
+def _respell(text: str, old: str, new: str) -> str:
+    """``text`` with the first token ``old`` of each row (any case) spelled ``new``."""
+    token = re.compile(rf"(?i)(?<!\S){re.escape(old)}(?!\S)")
+    return "\n".join(
+        line if line.startswith("--") else token.sub(new, line, count=1)
+        for line in text.split("\n")
+    )
+
+
+def _static_check(text: str, *, keep_dtm: bool = False) -> str:
+    """``text`` as a static check: TMax 0, and without dtM unless it is kept."""
+    text = re.sub(r"(?mi)^\s*\S+(\s+TMax\b)", r"0.0\1", text, count=1)
+    return text if keep_dtm else re.sub(r"(?mi)^\s*\S+\s+dtM\b.*\n", "", text, count=1)
+
+
+def _shipped_text(variant: str, name: str) -> str:
+    text = (EXAMPLES / name).read_text(encoding="utf-8")
+    if variant == "static":
+        return _static_check(text)
+    text = _static_check(text, keep_dtm=True)
+    for old, new in _KEYWORD_SPELLINGS:
+        text = _respell(text, old, new if variant == "spelled" else old.upper())
+    if variant == "upper":
+        for kind in ("airy", "jonswap", "pm", "issc", "torsethaugen", "ochihubble", "stream"):
+            text = _respell(text, kind, kind.upper())
+    return text
+
+
+def _deck_text(deck: str, rows: tuple[str, ...]) -> str:
+    if deck.startswith("shipped:"):
+        return _shipped_text(*deck.split(":")[1:])
+    if deck.endswith("_static"):
+        return _static_check(_deck_text(deck[: -len("_static")], rows))
+    if deck == "bed":
+        return _options(_HELD, "b.txt bathymetryFile", *rows).replace("50.0 WtrDpth\n", "")
+    if deck == "syrope":
+        return (EXAMPLES / "syrope_polyester_mooring.dat").read_text(encoding="utf-8")
+    if deck == "objects":
+        return _OBJECTS.replace("1.0 TMax\n", "1.0 TMax\n" + "".join(f"{r}\n" for r in rows))
+    if deck.endswith(".dat"):
+        text = (EXAMPLES / deck).read_text(encoding="utf-8")
+        return re.sub(r"(?m)^(\S+\s+TMax\b.*)$", "\\1\n" + "\n".join(rows), text, count=1)
+    return _options(_HELD, *rows)
+
+
+def _stage(folder: Path, rows: tuple[str, ...], files: dict[str, str], deck: str = "held") -> Path:
+    folder.mkdir(parents=True, exist_ok=True)
+    for name, text in files.items():
+        (folder / name).parent.mkdir(parents=True, exist_ok=True)
+        (folder / name).write_text(text, encoding="utf-8")
+    if deck.startswith("shipped:"):
+        source = (EXAMPLES / deck.split(":")[2]).parent
+        for item in source.iterdir():
+            if item.is_dir():
+                shutil.copytree(item, folder / item.name, dirs_exist_ok=True)
+            elif item.suffix != ".dat":
+                shutil.copy2(item, folder / item.name)
+    path = folder / "deck.dat"
+    path.write_text(_deck_text(deck, rows), encoding="utf-8")
+    return path
+
+
+def _accepted(deck: Path) -> bool:
+    try:
+        DeckFile.read(deck)
+    except DeckFormatError:
+        return False
+    return True
+
+
+@pytest.mark.skipif(
+    not os.environ.get("CABLEDYN_DRIVER", "").strip(),
+    reason="needs the native driver (CABLEDYN_DRIVER)",
+)
+@pytest.mark.parametrize(("base", "rows", "files"), _ALL_FILE_CASES)
+def test_wave_record_rules_match_the_driver(tmp_path, base, rows, files):
+    deck = _stage(tmp_path, rows, files, base)
+    run = subprocess.run(
+        [os.environ["CABLEDYN_DRIVER"], str(deck), str(tmp_path / "out")],
+        capture_output=True,
+        text=True,
+        errors="replace",
+        env={**os.environ, "OMP_NUM_THREADS": "1", "OPENBLAS_NUM_THREADS": "1"},
+        check=False,
+    )
+    assert _accepted(deck) == (run.returncode == 0), run.stdout[-400:] + run.stderr[-400:]
+
+
+@pytest.mark.parametrize(
+    ("rows", "files", "message"),
+    [
+        (
+            ("wk.dat WaterKin", "1 WaveRefX"),
+            _waterkin("0"),
+            "reference point of an elevation record",
+        ),
+        (
+            ("wk.dat WaterKin", "1 WaveRefX"),
+            _waterkin("SEASTATE"),
+            "reference point of an elevation",
+        ),
+        (("wk.dat WaterKin", "1 WaveRefX"), _waterkin("3"), "WaveKinMod supports 0"),
+        (("wk.dat WaterKin", "1 WaveRefX"), _waterkin("one"), "must be 0, 1, 2, or SEASTATE"),
+        (("wk.dat WaterKin", "1 WaveRefX"), {"wk.dat": "a\nb\n"}, "ended inside its header"),
+        (("wk.dat WaterKin", "1 WaveRefX"), {"wk.dat": "a\nb\nc\n"}, "no WaveKinMod line"),
+        (("wk.dat WaterKin", "1 WaveRefX"), {"wk.dat": "a\nb\nc\n \n"}, "line has no value"),
+        (("7 WaterKin",), {"wave_frequencies.txt": "0 0 0 0\n0.5 0.5 0 6.3\n"}, "within \\+-2 pi"),
+        (("7 WaterKin",), {"wave_frequencies.txt": "0 0 0 0\n"}, "at least two"),
+        (("7 WaterKin",), {"wave_frequencies.txt": "1 0 0 0\n2 1 0 0\n"}, "at omega = 0"),
+        (("7 WaterKin",), {"wave_frequencies.txt": "0 0 0 0\n0 1 0 0\n"}, "must be positive"),
+        (("7 WaterKin",), {"wave_frequencies.txt": "0 0 0 0\n0.5 x 0 0\n"}, "plain numbers"),
+        (("7 WaterKin",), {"wave_frequencies.txt": "0 0 0 0\n0.5 2000 0 0\n"}, "at most 1e3 m"),
+        (("7 WaterKin",), {"wave_frequencies.txt": "0 0 -60 0\n0.5 1 0 0\n"}, "effective depth"),
+    ],
+)
+def test_wave_files_next_to_the_deck_are_read(tmp_path, rows, files, message):
+    with pytest.raises(DeckFormatError, match=message):
+        DeckFile.read(_stage(tmp_path, rows, files))
+
+
+@pytest.mark.parametrize("kind", ["ochihubble", "ochi-hubble", "ochi_hubble", "OCHI_HUBBLE"])
+def test_every_spelling_of_a_kind_takes_its_heading_from_its_own_column(tmp_path, kind):
+    rows = (f"{kind} 2 8 1 1 12 1 30 waves", *_RAO_ROWS)
+    DeckFile.read(_stage(tmp_path, rows, {"rao.txt": _rao(heads=(30.0,))}))
+    with pytest.raises(DeckFormatError, match="one-heading RAO table"):
+        DeckFile.read(_stage(tmp_path, rows, {"rao.txt": _rao(heads=(31.0,))}))
+
+
+def test_a_static_check_reads_a_record_on_its_own_clock(tmp_path):
+    # TMax = 0 without dtM: the record resamples on its dtWave, whatever the nominal step
+    files = {"wk.dat": _kin("1"), "eta.txt": _ETA}
+    DeckFile.read(_stage(tmp_path, _WK, files, "held_static"))
+    with pytest.raises(DeckFormatError, match="four `time elevation` rows"):
+        DeckFile.read(_stage(tmp_path, _WK, {"eta.txt": "0 0\n1 0\n"}, "held_static"))
+
+
+def test_wave_files_that_are_not_there_are_left_to_the_solver(tmp_path):
+    # In memory under the nominal path, or on disk ahead of its files: the deck passes here,
+    # and the driver, which must open the file, refuses it by name while it is missing.
+    for rows in (("wk.dat WaterKin", "1 WaveRefX"), ("7 WaterKin", "1 WaveRefX")):
+        DeckFile.from_text(_options(_HELD, *rows))
+        DeckFile.from_text(_options(_HELD, *rows), path=tmp_path / "deck.dat")
+    # once the file is there the same text is judged by it
+    (tmp_path / "wk.dat").write_text(_WATERKIN.format(mode="0"), encoding="utf-8")
+    with pytest.raises(DeckFormatError, match="reference point of an elevation record"):
+        DeckFile.from_text(
+            _options(_HELD, "wk.dat WaterKin", "1 WaveRefX"), path=tmp_path / "deck.dat"
+        )
+
+
+def test_a_waterkin_file_is_found_as_the_native_reader_finds_it(tmp_path):
+    # quoted, in a subfolder, and by an absolute path
+    (tmp_path / "sub").mkdir()
+    (tmp_path / "sub" / "wk.dat").write_text(_WATERKIN.format(mode="0"), encoding="utf-8")
+    for value in ('"sub/wk.dat"', "sub/wk.dat", (tmp_path / "sub" / "wk.dat").as_posix()):
+        with pytest.raises(DeckFormatError, match="reference point of an elevation record"):
+            DeckFile.read(_stage(tmp_path, (f"{value} WaterKin", "1 WaveRefX"), {}))
+
+
 @pytest.mark.parametrize("value", ["0", "-0.1"])
 def test_dtwave_must_be_positive(value):
     _rejects(_options(_HELD, f"{value} dtWave"), "OPTION dtWave must be positive")
@@ -105,6 +771,16 @@ def test_dtwave_must_be_positive(value):
         (("uniform 0.1 0 0 current", "1 Currents"), "both a current OPTION and Currents 1"),
         (("uniform 0.1 0 0 current", "2 Currents"), "both a current OPTION and Currents 2"),
         (("uniform 0.1 0 0 current", "5 Currents"), "both a current OPTION and Currents 5"),
+        (("7 WaterKin", "12.5 WaveRefX", "-3 WaveRefY"), None),
+        (("3 WaveKin", "0.1 dtWave", "0 WaveRefY"), None),
+        (("wk.dat WaterKin", "12.5 WaveRefX"), None),
+        (("12.5 WaveRefX",), "reference point of an elevation record"),
+        (("0 WaveRefY",), "reference point of an elevation record"),
+        (("airy 2 8 0 waves", "1 WaveRefX"), "reference point of an elevation record"),
+        (("airy 2 8 0 wavetrain", "1 WaveRefY"), "reference point of an elevation record"),
+        (("7 WaterKin", "2e6 WaveRefX"), "at most 1e6 m"),
+        (("7 WaterKin", "nan WaveRefY"), "must be"),
+        (("7 WaterKin", "north WaveRefY"), "must be numeric"),
     ],
 )
 def test_moordyn_c_kinematics_follow_the_native_rules(rows, message):
@@ -633,7 +1309,7 @@ def test_generated_cases_carry_the_fixed_name_kinematics_files(tmp_path):
     ("mode", "name"), [("2", "current_profile_dynamic.txt"), ("5", "current_profile_4d.txt")]
 )
 def test_generated_cases_carry_the_moordyn_c_current_tables(tmp_path, mode, name):
-    table = "a MoorDyn-C current table\n"
+    table = _CURRENT_TABLES[mode]
     base = _kinematic_base(tmp_path / "src", f"{mode} Currents", files={name: table})
     out = tmp_path / "cases"
     generate_deck_cases(base, out, {"a": {"option.dtM": 0.05}})
diff --git a/src/CableDyn_DeckDriver.f90 b/src/CableDyn_DeckDriver.f90
index d4c02c2..9c43c02 100644
--- a/src/CableDyn_DeckDriver.f90
+++ b/src/CableDyn_DeckDriver.f90
@@ -664,7 +664,14 @@ MODULE CableDyn_DeckDriver
     REAL(wp), ALLOCATABLE :: wave_component_amplitude(:), wave_component_phase(:)
     ! exp(-2 k (wtrdpth + wave_eta_mean)) and the tail_sums of the table (CD_Component_Wave_Kinematics)
     REAL(wp), ALLOCATABLE :: wave_component_e2kd(:), wave_component_tail(:, :)
+    ! a short-crested table (WaveKin 7 rows of different directions): each component's direction
+    ! cosine and sine, in ascending frequency; unallocated when all share wave_direction
+    REAL(wp), ALLOCATABLE :: wave_component_cos(:), wave_component_sin(:)
     REAL(wp) :: wave_eta_mean = CD_ZERO
+    ! OPTION WaveRefX, WaveRefY: the point the phases of an elevation record or a component file
+    ! refer to (the record was measured there); the origin by default
+    REAL(wp) :: wave_ref_x = CD_ZERO, wave_ref_y = CD_ZERO
+    LOGICAL :: has_wave_ref = .FALSE.
     ! N-level vertical current profile (ascending z): from the two-point deck
     ! OPTION row or a MoorDyn WaterKin file's CurrentMod-1 depth table.
     REAL(wp), ALLOCATABLE :: current_profile_z(:)
@@ -13802,8 +13809,10 @@ SUBROUTINE read_moordyn_c_kinematics(dir, opts, ErrStat, ErrMsg)
     !!    the mean-level offset, as MoorDyn-C does. Waves travel along +x.
     !!  - WaveKin 7, wave_frequencies.txt (`omega Re Im [beta]` rows, the first at omega = 0):
     !!    each row is the component |c| sin(omega t - k d + arg c) of MoorDyn-C's SUM_COMPONENTS
-    !!    mode, d = x cos(beta) + y sin(beta); the omega = 0 row is the mean-level offset. The rows
-    !!    must share one direction beta.
+    !!    mode, d = x cos(beta) + y sin(beta) with the row's own direction beta; the omega = 0 row
+    !!    is the mean-level offset and its beta has no effect. Rows that share one direction are
+    !!    a long-crested table in file order; rows that differ a short-crested table in ascending
+    !!    frequency, the order the tail bound of the sums needs.
     !!  - Currents 1, current_profile.txt (`z ux uy uz` rows after its header rows): the steady
     !!    profile, linear in z and held beyond its ends, as MoorDyn-C's current grid.
     !!  - Currents 2, current_profile_dynamic.txt (constructDynamicCurrentGrid of MoorDyn-C's
@@ -13823,9 +13832,11 @@ SUBROUTINE read_moordyn_c_kinematics(dir, opts, ErrStat, ErrMsg)
     TYPE(DeckOptions), INTENT(INOUT) :: opts
     INTEGER, INTENT(OUT) :: ErrStat
     CHARACTER(*), INTENT(OUT) :: ErrMsg
-    REAL(wp), ALLOCATABLE :: tab(:, :), eta(:), comp_w(:), comp_a(:), comp_p(:)
+    REAL(wp), ALLOCATABLE :: tab(:, :), eta(:), comp_w(:), comp_a(:), comp_p(:), comp_b(:)
     REAL(wp) :: dt, tl, frac, ca, sa, angle, scale, beta, mean_level, arg
     INTEGER :: n, nt, j, k, nc, bracket
+    INTEGER, ALLOCATABLE :: order(:)
+    LOGICAL :: one_direction
     INTEGER(KIND=8) :: kj
     CHARACTER(32) :: fname
     ! the longest row of a MoorDyn-C current table (a time and 3 nz velocities per row)
@@ -13933,15 +13944,15 @@ SUBROUTINE read_moordyn_c_kinematics(dir, opts, ErrStat, ErrMsg)
       IF (ABS(tab(1, 1)) > CD_ZERO) THEN
         CALL fail(ErrStat, ErrMsg, 'WaveKin 7: the first wave_frequencies.txt row must be at omega = 0'); RETURN
       END IF
+      ! the omega = 0 row has no direction of its own: the first component's stands for the table
       beta = tab(4, 1)
+      IF (ABS(tab(4, 2) - beta) > 1.0e-12_wp) beta = tab(4, 2)
+      one_direction = .TRUE.
       DO j = 1, n
         IF (ABS(tab(4, j)) > TWO_PI) THEN
           CALL fail(ErrStat, ErrMsg, 'WaveKin 7: wave_frequencies.txt directions are radians within +-2 pi'); RETURN
         END IF
-        IF (ABS(tab(4, j) - beta) > 1.0e-12_wp) THEN
-          CALL fail(ErrStat, ErrMsg, 'WaveKin 7: wave_frequencies.txt rows with different directions '// &
-                    '(short-crested seas) are not supported; give every row one direction'); RETURN
-        END IF
+        IF (j > 1 .AND. ABS(tab(4, j) - beta) > 1.0e-12_wp) one_direction = .FALSE.
         IF (j > 1) THEN
           IF (.NOT. tab(1, j) > CD_ZERO) THEN
             CALL fail(ErrStat, ErrMsg, 'WaveKin 7: wave_frequencies.txt frequencies after the first row '// &
@@ -13954,6 +13965,14 @@ SUBROUTINE read_moordyn_c_kinematics(dir, opts, ErrStat, ErrMsg)
       comp_w = tab(1, 2:n)
       comp_a = SQRT(tab(2, 2:n)**2 + tab(3, 2:n)**2)
       comp_p = 0.25_wp*TWO_PI - ATAN2(tab(3, 2:n), tab(2, 2:n))
+      IF (.NOT. one_direction) THEN
+        ALLOCATE (order(n - 1))
+        CALL sort_index(comp_w, order)
+        comp_w = comp_w(order)
+        comp_a = comp_a(order)
+        comp_p = comp_p(order)
+        comp_b = tab(4, order + 1)
+      END IF
       CALL install_components()
       IF (ErrStat /= CD_DECKDRV_OK) RETURN
     END SELECT
@@ -14378,7 +14397,7 @@ END FUNCTION real_to_str
 
     SUBROUTINE install_components()
       !! The component waves (omega, amplitude, phase in CableDyn's cos(k d - omega t + phase)
-      !! convention), their direction beta [rad] and the mean level.
+      !! convention), their direction beta [rad], or each their own comp_b, and the mean level.
       IF (.NOT. (ALL(IEEE_IS_FINITE(comp_w)) .AND. ALL(IEEE_IS_FINITE(comp_a)) .AND. &
                  ALL(IEEE_IS_FINITE(comp_p)) .AND. IEEE_IS_FINITE(mean_level))) THEN
         CALL fail(ErrStat, ErrMsg, 'WaveKin '//TRIM(int_to_str(opts%mdc_wave_kin))//': '//TRIM(fname)// &
@@ -14400,6 +14419,11 @@ SUBROUTINE install_components()
       opts%wave_component_phase = comp_p
       opts%wave_eta_mean = mean_level
       opts%wave_direction = beta*45.0_wp/ATAN(CD_ONE)
+      IF (ALLOCATED(opts%wave_component_cos)) DEALLOCATE (opts%wave_component_cos, opts%wave_component_sin)
+      IF (ALLOCATED(comp_b)) THEN
+        opts%wave_component_cos = COS(comp_b)
+        opts%wave_component_sin = SIN(comp_b)
+      END IF
       opts%wave_kind = WAVE_COMPONENTS
       opts%has_wave = .TRUE.
     END SUBROUTINE install_components
@@ -16620,6 +16644,12 @@ SUBROUTINE apply_option(line, opts, ErrStat, ErrMsg)
       ELSE
         CALL rd(val, opts%wv_hi_cutoff, ErrStat, ErrMsg)
       END IF
+    CASE ('waverefx')
+      CALL rd(val, opts%wave_ref_x, ErrStat, ErrMsg)
+      opts%has_wave_ref = .TRUE.
+    CASE ('waverefy')
+      CALL rd(val, opts%wave_ref_y, ErrStat, ErrMsg)
+      opts%has_wave_ref = .TRUE.
     CASE ('ramptime', 'ramp_time', 'tramp')
       CALL rd(val, opts%ramp_time, ErrStat, ErrMsg)
     CASE ('maxstrain', 'max_strain')
@@ -17124,7 +17154,7 @@ SUBROUTINE check_option_row_value(kw, val, ErrStat, ErrMsg)
       rule = 'non-negative'
     CASE ('aslaidtension', 'as_laid_tension')
       rule = 'non-negative'
-    CASE ('layazimuth', 'lay_azimuth')
+    CASE ('layazimuth', 'lay_azimuth', 'waverefx', 'waverefy')
       rule = 'finite'
     CASE ('rhoinf', 'rho_inf')
       rule = 'unit'
@@ -17758,6 +17788,15 @@ SUBROUTINE finalize_checks(types, bodies, rod_types, rods, points, lines, sectio
         END IF
       END IF
     END IF
+    IF (opts%has_wave_ref) THEN
+      IF (ABS(opts%wave_ref_x) > MAX_COORDINATE .OR. ABS(opts%wave_ref_y) > MAX_COORDINATE) THEN
+        CALL fail(ErrStat, ErrMsg, 'WaveRefX and WaveRefY must be at most 1e6 m in magnitude'); RETURN
+      END IF
+      IF (.NOT. (opts%has_wave .AND. opts%wave_kind == WAVE_COMPONENTS)) THEN
+        CALL fail(ErrStat, ErrMsg, 'WaveRefX and WaveRefY place the reference point of an elevation record or a '// &
+                  'component file (WaterKin WaveKinMod 1, WaveKin 3 or 7); this deck has none'); RETURN
+      END IF
+    END IF
     IF (opts%has_wave) THEN
       IF (.NOT. opts%has_dtm) THEN
         CALL fail(ErrStat, ErrMsg, 'deck waves require dtM and TMax'); RETURN
@@ -32232,21 +32271,23 @@ SUBROUTINE deck_wave_kinematics(opts, x, y, z, time, eta, velocity, acceleration
                                 pdyn)
     CASE (WAVE_COMPONENTS)
       IF (ramp < CD_ONE) THEN
-        CALL CD_Component_Wave_Kinematics(x, y, z - opts%wave_eta_mean, time, &
+        CALL CD_Component_Wave_Kinematics(x - opts%wave_ref_x, y - opts%wave_ref_y, z - opts%wave_eta_mean, time, &
                                           opts%wtrdpth + opts%wave_eta_mean, opts%wave_direction, .TRUE., &
                                           opts%wave_component_omega, opts%wave_component_k, &
                                           opts%wave_component_amplitude, opts%wave_component_phase, &
                                           eta, velocity, acceleration, ErrStat, ErrMsg, inputs_validated=.TRUE., &
                                           amplitude_scale=ramp, pdyn=pdyn, exp_m2kd=opts%wave_component_e2kd, &
-                                          tail_sums=opts%wave_component_tail)
+                                          tail_sums=opts%wave_component_tail, comp_cos=opts%wave_component_cos, &
+                                          comp_sin=opts%wave_component_sin)
       ELSE
-        CALL CD_Component_Wave_Kinematics(x, y, z - opts%wave_eta_mean, time, &
+        CALL CD_Component_Wave_Kinematics(x - opts%wave_ref_x, y - opts%wave_ref_y, z - opts%wave_eta_mean, time, &
                                           opts%wtrdpth + opts%wave_eta_mean, opts%wave_direction, .TRUE., &
                                           opts%wave_component_omega, opts%wave_component_k, &
                                           opts%wave_component_amplitude, opts%wave_component_phase, &
                                           eta, velocity, acceleration, ErrStat, ErrMsg, inputs_validated=.TRUE., &
                                           pdyn=pdyn, exp_m2kd=opts%wave_component_e2kd, &
-                                          tail_sums=opts%wave_component_tail)
+                                          tail_sums=opts%wave_component_tail, comp_cos=opts%wave_component_cos, &
+                                          comp_sin=opts%wave_component_sin)
       END IF
       eta = eta + opts%wave_eta_mean
     CASE DEFAULT
@@ -35097,6 +35138,10 @@ SUBROUTINE synthesize_vessel_rao(path, opts, motion, rref, vref, aref, ErrStat,
     CASE (WAVE_COMPONENTS)
       ALLOCATE (om, SOURCE=opts%wave_component_omega); ALLOCATE (kk, SOURCE=opts%wave_component_k)
       ALLOCATE (amp, SOURCE=opts%wave_component_amplitude); ALLOCATE (phs, SOURCE=opts%wave_component_phase)
+      IF (ALLOCATED(opts%wave_component_cos)) THEN
+        ALLOCATE (dir(SIZE(om)))
+        dir = ATAN2(opts%wave_component_sin, opts%wave_component_cos)*R2D
+      END IF
     CASE (WAVE_SEA)
       nc = opts%sea%n_comp
       ALLOCATE (om, SOURCE=opts%sea%omega(1:nc)); ALLOCATE (kk, SOURCE=opts%sea%k(1:nc))
@@ -35122,7 +35167,7 @@ SUBROUTINE synthesize_vessel_rao(path, opts, motion, rref, vref, aref, ErrStat,
     DO i = 1, nc
       cb = COS(dir(i)*QUARTER_PI/45.0_wp)
       sb = SIN(dir(i)*QUARTER_PI/45.0_wp)
-      eps(i) = kk(i)*(opts%vessel_ref(1)*cb + opts%vessel_ref(2)*sb) + phs(i)
+      eps(i) = kk(i)*((opts%vessel_ref(1) - opts%wave_ref_x)*cb + (opts%vessel_ref(2) - opts%wave_ref_y)*sb) + phs(i)
       IF (.NOT. (amp(i) > CD_ZERO)) CYCLE
       CALL CD_Vessel_RAO_Eval(rao, TWO_PI/om(i), dir(i), cre(:, i), cim(:, i), clamped, es, em)
       IF (es /= CD_VESSEL_OK) THEN
diff --git a/src/CableDyn_Hydro.f90 b/src/CableDyn_Hydro.f90
index 487e5c5..2a6789a 100644
--- a/src/CableDyn_Hydro.f90
+++ b/src/CableDyn_Hydro.f90
@@ -2239,7 +2239,7 @@ END SUBROUTINE CD_JONSWAP_Wave_Kinematics_Precomputed
   SUBROUTINE CD_Component_Wave_Kinematics(x, y, z, t, depth, direction_deg, stretch, &
                                           omega_comp, k_comp, amplitude, phase, &
                                           eta, velocity, acceleration, ErrStat, ErrMsg, inputs_validated, &
-                                          amplitude_scale, pdyn, exp_m2kd, tail_sums)
+                                          amplitude_scale, pdyn, exp_m2kd, tail_sums, comp_cos, comp_sin)
     !! Evaluate linear irregular-wave kinematics from a CALLER-SUPPLIED component table:
     !! the general (arbitrary length, explicit per-component phase) form of the fixed-size
     !! deterministic-phase JONSWAP evaluator above, sharing its conventions exactly --
@@ -2275,6 +2275,9 @@ SUBROUTINE CD_Component_Wave_Kinematics(x, y, z, t, depth, direction_deg, stretc
     !! once that bound is below half an ulp of the magnitude already summed: a deep point drops
     !! the short components that have decayed to nothing there, and the result changes at
     !! round-off only. The elevation always sums every component.
+    !! comp_cos and comp_sin, given together, hold the cosine and sine of each component's own
+    !! direction, in place of the one direction_deg of a long-crested table: a short-crested
+    !! table. The bound of tail_sums holds for any direction, since it bounds the magnitude.
     REAL(wp), INTENT(IN) :: x, y, z, t, depth, direction_deg
     LOGICAL, INTENT(IN) :: stretch
     REAL(wp), INTENT(IN) :: omega_comp(:), k_comp(:), amplitude(:), phase(:)
@@ -2286,7 +2289,7 @@ SUBROUTINE CD_Component_Wave_Kinematics(x, y, z, t, depth, direction_deg, stretc
     ! the linear dynamic pressure per unit density (sum over the components of
     ! (omega^2/k) a cosh(k(z+h))/sinh(kh) cos(theta)), as in CD_Airy_Wave_Kinematics_Precomputed
     REAL(wp), INTENT(OUT), OPTIONAL :: pdyn
-    REAL(wp), INTENT(IN), OPTIONAL :: exp_m2kd(:), tail_sums(:, :)
+    REAL(wp), INTENT(IN), OPTIONAL :: exp_m2kd(:), tail_sums(:, :), comp_cos(:), comp_sin(:)
     INTEGER :: n
     REAL(wp) :: cb, sb, scl
     LOGICAL :: trust_inputs
@@ -2316,6 +2319,22 @@ SUBROUTINE CD_Component_Wave_Kinematics(x, y, z, t, depth, direction_deg, stretc
         RETURN
       END IF
     END IF
+    IF (PRESENT(comp_cos) .NEQV. PRESENT(comp_sin)) THEN
+      CALL fail(ErrStat, ErrMsg, 'comp_cos and comp_sin must be given together')
+      RETURN
+    END IF
+    IF (PRESENT(comp_cos)) THEN
+      IF (SIZE(comp_cos) /= n .OR. SIZE(comp_sin) /= n) THEN
+        CALL fail(ErrStat, ErrMsg, 'comp_cos and comp_sin must match the component-wave arrays')
+        RETURN
+      END IF
+      IF (.NOT. trust_inputs) THEN
+        IF (.NOT. (CD_All_Finite(comp_cos) .AND. CD_All_Finite(comp_sin))) THEN
+          CALL fail(ErrStat, ErrMsg, 'comp_cos and comp_sin must be finite')
+          RETURN
+        END IF
+      END IF
+    END IF
     ! Fortran does not require short-circuit evaluation: keep the O(n) scans inside a
     ! separate branch so the setup-validated hot path portably skips them at -O0 too.
     IF (.NOT. trust_inputs) THEN
@@ -2345,7 +2364,7 @@ SUBROUTINE CD_Component_Wave_Kinematics(x, y, z, t, depth, direction_deg, stretc
     ! a unit scale is exact, so the unscaled table is reproduced bit for bit
     scl = CD_ONE
     IF (PRESENT(amplitude_scale)) scl = amplitude_scale
-    IF (n == 1 .AND. ABS(phase(1)) <= CD_ZERO .AND. scl*amplitude(1) > CD_ZERO) THEN
+    IF (n == 1 .AND. ABS(phase(1)) <= CD_ZERO .AND. scl*amplitude(1) > CD_ZERO .AND. .NOT. PRESENT(comp_cos)) THEN
       CALL CD_Airy_Wave_Kinematics_Precomputed(x, y, z, t, 2.0_wp*(scl*amplitude(1)), omega_comp(1), k_comp(1), &
                                                depth, direction_deg, stretch, eta, velocity, acceleration, &
                                                ErrStat, ErrMsg, pdyn)
@@ -2353,25 +2372,28 @@ SUBROUTINE CD_Component_Wave_Kinematics(x, y, z, t, depth, direction_deg, stretc
     END IF
     cb = COS(direction_deg*DEG2RAD)
     sb = SIN(direction_deg*DEG2RAD)
-    CALL component_sums(n, x*cb + y*sb, z, t, depth, stretch, scl, cb, sb, omega_comp, k_comp, amplitude, phase, eta, &
-                        velocity, acceleration, ErrStat, ErrMsg, pdyn, exp_m2kd, tail_sums)
+    CALL component_sums(n, x, y, z, t, depth, stretch, scl, cb, sb, omega_comp, k_comp, amplitude, phase, eta, &
+                        velocity, acceleration, ErrStat, ErrMsg, pdyn, exp_m2kd, tail_sums, comp_cos, comp_sin)
   END SUBROUTINE CD_Component_Wave_Kinematics
 
-  SUBROUTINE component_sums(n, x_along, z, t, depth, stretch, scl, cb, sb, omega_comp, k_comp, amplitude, phase, eta, &
-                            velocity, acceleration, ErrStat, ErrMsg, pdyn, exp_m2kd, tail_sums)
+  SUBROUTINE component_sums(n, x, y, z, t, depth, stretch, scl, cb, sb, omega_comp, k_comp, amplitude, phase, eta, &
+                            velocity, acceleration, ErrStat, ErrMsg, pdyn, exp_m2kd, tail_sums, comp_cos, comp_sin)
     !! The sums of CD_Component_Wave_Kinematics, on explicit-shape tables: their blocks pass to
-    !! the block kernels without a copy, and the loops run at unit stride.
+    !! the block kernels without a copy, and the loops run at unit stride. A long-crested table
+    !! forms each phase from the one distance along its direction and sums the horizontal
+    !! velocity and acceleration once, for cb and sb to multiply; a short-crested one (comp_cos,
+    !! comp_sin) forms each phase from the component's own direction and sums x and y apart.
     INTEGER, INTENT(IN) :: n
-    REAL(wp), INTENT(IN) :: x_along, z, t, depth, scl, cb, sb
+    REAL(wp), INTENT(IN) :: x, y, z, t, depth, scl, cb, sb
     LOGICAL, INTENT(IN) :: stretch
     REAL(wp), INTENT(IN) :: omega_comp(n), k_comp(n), amplitude(n), phase(n)
     REAL(wp), INTENT(INOUT) :: eta, velocity(3), acceleration(3)
     INTEGER, INTENT(INOUT) :: ErrStat
     CHARACTER(*), INTENT(INOUT) :: ErrMsg
     REAL(wp), INTENT(INOUT), OPTIONAL :: pdyn
-    REAL(wp), INTENT(IN), OPTIONAL :: exp_m2kd(n), tail_sums(3, n)
+    REAL(wp), INTENT(IN), OPTIONAL :: exp_m2kd(n), tail_sums(3, n), comp_cos(n), comp_sin(n)
     INTEGER :: nb, n_small, i0, i1, j, nk, m, lanes
-    REAL(wp) :: z_eval
+    REAL(wp) :: z_eval, x_along
     ! the velocity sums run in blocks of NLANE and the elevation in chunks of NCHUNK (NBUF a
     ! multiple of NCHUNK and NCHUNK of NLANE, so a block lies wholly inside or past the buffer)
     INTEGER, PARAMETER :: NBUF = 4096, NCHUNK = 512
@@ -2379,10 +2401,12 @@ SUBROUTINE component_sums(n, x_along, z, t, depth, stretch, scl, cb, sb, omega_c
     REAL(wp) :: denom, e2, first_ratio
     REAL(wp) :: cth(NBUF), sth(NBUF), th(NCHUNK), bc(NCHUNK), bs(NCHUNK), ch(NLANE), sh(NLANE), e_bed(NLANE)
     REAL(wp) :: lane_eta(NLANE), lane_vc(NLANE), lane_vs(NLANE), lane_as(NLANE), lane_ac(NLANE)
-    REAL(wp) :: lane_p(NLANE), lane_m(NLANE, 3), sum_vc, sum_vs, sum_as, sum_ac
+    REAL(wp) :: lane_p(NLANE), lane_m(NLANE, 3), lane_vy(NLANE), lane_ay(NLANE), sum_vc, sum_vs, sum_as, sum_ac
     REAL(wp) :: mag(3), rest(3)
-    LOGICAL :: at_bed, use_tail, want_p
+    LOGICAL :: at_bed, use_tail, want_p, spread
 
+    spread = PRESENT(comp_cos)
+    x_along = x*cb + y*sb
     nb = MIN(n, NBUF)
     lanes = MIN(n, NLANE)
     ! the elevation sums every component (Wheeler stretching needs the total eta); the phase
@@ -2391,9 +2415,15 @@ SUBROUTINE component_sums(n, x_along, z, t, depth, stretch, scl, cb, sb, omega_c
     DO i0 = 1, n, NCHUNK
       nk = MIN(n - i0 + 1, NCHUNK)
       i1 = i0 + nk - 1
-      DO j = 1, nk
-        th(j) = k_comp(i0 + j - 1)*x_along - omega_comp(i0 + j - 1)*t + phase(i0 + j - 1)
-      END DO
+      IF (spread) THEN
+        DO j = i0, i1
+          th(j - i0 + 1) = k_comp(j)*(x*comp_cos(j) + y*comp_sin(j)) - omega_comp(j)*t + phase(j)
+        END DO
+      ELSE
+        DO j = 1, nk
+          th(j) = k_comp(i0 + j - 1)*x_along - omega_comp(i0 + j - 1)*t + phase(i0 + j - 1)
+        END DO
+      END IF
       IF (i0 <= nb) THEN
         CALL CD_Wave_SinCos(nk, th, sth(i0), cth(i0))
         CALL add_elevation_lanes(nk, amplitude(i0:i1), cth(i0), lane_eta)
@@ -2429,6 +2459,10 @@ SUBROUTINE component_sums(n, x_along, z, t, depth, stretch, scl, cb, sb, omega_c
     lane_ac(1:lanes) = CD_ZERO
     lane_p(1:lanes) = CD_ZERO
     lane_m(1:lanes, :) = CD_ZERO
+    IF (spread) THEN
+      lane_vy(1:lanes) = CD_ZERO
+      lane_ay(1:lanes) = CD_ZERO
+    END IF
     first_ratio = CD_ZERO
     DO i0 = 1, n, NLANE
       nk = MIN(n - i0 + 1, NLANE)
@@ -2479,20 +2513,42 @@ SUBROUTINE component_sums(n, x_along, z, t, depth, stretch, scl, cb, sb, omega_c
           END DO
         END IF
       END IF
-      IF (i0 <= nb) THEN
+      IF (i0 > nb) THEN
+        IF (spread) THEN
+          DO j = i0, i1
+            th(j - i0 + 1) = k_comp(j)*(x*comp_cos(j) + y*comp_sin(j)) - omega_comp(j)*t + phase(j)
+          END DO
+        ELSE
+          DO j = 1, nk
+            th(j) = k_comp(i0 + j - 1)*x_along - omega_comp(i0 + j - 1)*t + phase(i0 + j - 1)
+          END DO
+        END IF
+        CALL CD_Wave_SinCos(nk, th, bs, bc)
+      END IF
+      IF (spread) THEN
+        IF (i0 <= nb) THEN
+          CALL add_spread_lanes(nk, amplitude(i0:i1), omega_comp(i0:i1), comp_cos(i0:i1), comp_sin(i0:i1), ch, sh, &
+                                cth(i0), sth(i0), lane_vc, lane_vy, lane_vs, lane_as, lane_ay, lane_ac, lane_m(1, 1), &
+                                lane_m(1, 2))
+        ELSE
+          CALL add_spread_lanes(nk, amplitude(i0:i1), omega_comp(i0:i1), comp_cos(i0:i1), comp_sin(i0:i1), ch, sh, &
+                                bc, bs, lane_vc, lane_vy, lane_vs, lane_as, lane_ay, lane_ac, lane_m(1, 1), &
+                                lane_m(1, 2))
+        END IF
+      ELSE IF (i0 <= nb) THEN
         CALL add_velocity_lanes(nk, amplitude(i0:i1), omega_comp(i0:i1), ch, sh, cth(i0), sth(i0), lane_vc, &
                                 lane_vs, lane_as, lane_ac, lane_m(1, 1), lane_m(1, 2))
-        IF (want_p) CALL add_pressure_lanes(nk, amplitude(i0:i1), omega_comp(i0:i1), k_comp(i0:i1), ch, cth(i0), &
-                                            lane_p, lane_m(1, 3))
       ELSE
-        DO j = 1, nk
-          th(j) = k_comp(i0 + j - 1)*x_along - omega_comp(i0 + j - 1)*t + phase(i0 + j - 1)
-        END DO
-        CALL CD_Wave_SinCos(nk, th, bs, bc)
         CALL add_velocity_lanes(nk, amplitude(i0:i1), omega_comp(i0:i1), ch, sh, bc, bs, lane_vc, lane_vs, lane_as, &
                                 lane_ac, lane_m(1, 1), lane_m(1, 2))
-        IF (want_p) CALL add_pressure_lanes(nk, amplitude(i0:i1), omega_comp(i0:i1), k_comp(i0:i1), ch, bc, lane_p, &
-                                            lane_m(1, 3))
+      END IF
+      IF (want_p) THEN
+        IF (i0 <= nb) THEN
+          CALL add_pressure_lanes(nk, amplitude(i0:i1), omega_comp(i0:i1), k_comp(i0:i1), ch, cth(i0), lane_p, &
+                                  lane_m(1, 3))
+        ELSE
+          CALL add_pressure_lanes(nk, amplitude(i0:i1), omega_comp(i0:i1), k_comp(i0:i1), ch, bc, lane_p, lane_m(1, 3))
+        END IF
       END IF
       IF (i0 == 1) first_ratio = ch(1)
       IF (use_tail) THEN
@@ -2515,8 +2571,13 @@ SUBROUTINE component_sums(n, x_along, z, t, depth, stretch, scl, cb, sb, omega_c
     sum_vs = scl*CD_Wave_Sum_Lanes(lane_vs, lanes)
     sum_as = scl*CD_Wave_Sum_Lanes(lane_as, lanes)
     sum_ac = scl*CD_Wave_Sum_Lanes(lane_ac, lanes)
-    velocity = [sum_vc*cb, sum_vc*sb, sum_vs]
-    acceleration = [sum_as*cb, sum_as*sb, -sum_ac]
+    IF (spread) THEN
+      velocity = [sum_vc, scl*CD_Wave_Sum_Lanes(lane_vy, lanes), sum_vs]
+      acceleration = [sum_as, scl*CD_Wave_Sum_Lanes(lane_ay, lanes), -sum_ac]
+    ELSE
+      velocity = [sum_vc*cb, sum_vc*sb, sum_vs]
+      acceleration = [sum_as*cb, sum_as*sb, -sum_ac]
+    END IF
     IF (want_p) pdyn = scl*CD_Wave_Sum_Lanes(lane_p, lanes)
   END SUBROUTINE component_sums
 
@@ -2557,6 +2618,34 @@ PURE SUBROUTINE add_velocity_lanes(n, amplitude, omega, ch, sh, c, s, lane_vc, l
     END DO
   END SUBROUTINE add_velocity_lanes
 
+  PURE SUBROUTINE add_spread_lanes(n, amplitude, omega, cosb, sinb, ch, sh, c, s, lane_vx, lane_vy, lane_vs, lane_ax, &
+                                   lane_ay, lane_ac, lane_mv, lane_ma)
+    !! add_velocity_lanes for components each with its own direction: the horizontal sums go to
+    !! x and y apart.
+    INTEGER, INTENT(IN) :: n
+    REAL(wp), INTENT(IN) :: amplitude(n), omega(n), cosb(n), sinb(n), ch(n), sh(n), c(n), s(n)
+    REAL(wp), INTENT(INOUT) :: lane_vx(n), lane_vy(n), lane_vs(n), lane_ax(n), lane_ay(n), lane_ac(n), lane_mv(n), &
+                               lane_ma(n)
+    INTEGER :: j
+    REAL(wp) :: wv, wa, gv, ga, uh
+    DO j = 1, n
+      wv = amplitude(j)*omega(j)
+      wa = wv*omega(j)
+      gv = wv*ch(j)
+      ga = wa*ch(j)
+      uh = gv*c(j)
+      lane_vx(j) = lane_vx(j) + uh*cosb(j)
+      lane_vy(j) = lane_vy(j) + uh*sinb(j)
+      lane_vs(j) = lane_vs(j) + (wv*sh(j))*s(j)
+      uh = ga*s(j)
+      lane_ax(j) = lane_ax(j) + uh*cosb(j)
+      lane_ay(j) = lane_ay(j) + uh*sinb(j)
+      lane_ac(j) = lane_ac(j) + (wa*sh(j))*c(j)
+      lane_mv(j) = lane_mv(j) + gv
+      lane_ma(j) = lane_ma(j) + ga
+    END DO
+  END SUBROUTINE add_spread_lanes
+
   PURE SUBROUTINE add_pressure_lanes(n, amplitude, omega, k, ch, c, lane_p, lane_mp)
     !! One block (n <= NLANE) of the pressure sum of CD_Component_Wave_Kinematics.
     INTEGER, INTENT(IN) :: n
diff --git a/src/CableDyn_WaveSpectra.f90 b/src/CableDyn_WaveSpectra.f90
index f2b45c0..5ff8461 100644
--- a/src/CableDyn_WaveSpectra.f90
+++ b/src/CableDyn_WaveSpectra.f90
@@ -233,7 +233,9 @@ END FUNCTION CD_Spreading_Constant
   PURE SUBROUTINE CD_Spreading_Weights(s, n, offset, weight)
     !! Discretise D(theta) = K(s) cos^2s(theta) on [-pi/2, pi/2] into n equal-angle bins:
     !! offset(j) [rad] is the bin centre and weight(j) the bin integral of D (composite
-    !! Simpson quadrature, 128 intervals per bin), normalised to sum 1.
+    !! Simpson quadrature, 128 intervals per bin), normalised to sum 1. A bin below the smallest
+    !! normal number holds no energy: its weight is exactly zero, so which bins carry the sea
+    !! depends on s and n alone, not on where the spectrum's product with it underflows.
     REAL(wp), INTENT(IN) :: s
     INTEGER, INTENT(IN) :: n
     REAL(wp), INTENT(OUT) :: offset(n), weight(n)
@@ -263,6 +265,7 @@ PURE SUBROUTINE CD_Spreading_Weights(s, n, offset, weight)
     END DO
     total = SUM(weight)
     IF (total > CD_ZERO) weight = weight/total
+    WHERE (weight < TINY(CD_ONE)) weight = CD_ZERO
   END SUBROUTINE CD_Spreading_Weights
 
   SUBROUTINE CD_Sea_Reset(sea)
diff --git a/tests/test_vessel_motion.f90 b/tests/test_vessel_motion.f90
index 2dd92d1..1f2616f 100644
--- a/tests/test_vessel_motion.f90
+++ b/tests/test_vessel_motion.f90
@@ -25,6 +25,9 @@ PROGRAM test_vessel_motion
   !!   7. Fail-closed option and file grammar.
   !!   8. Deck, vesselRAO in a spread sea: each spread component drives the vessel with its own
   !!      heading and phase; a nonlinear stream wave is rejected by name.
+  !!   9. Deck, vesselRAO in a WaveKin 7 sea of rows with different directions referred to a point
+  !!      (WaveRefX, WaveRefY): each row drives the vessel with its own heading, and its phase
+  !!      from the position of the vessel relative to that point.
   USE CableDyn_Precision, ONLY: wp
   USE CableDyn_Vessel, ONLY: CD_VesselRAOType, CD_Vessel_Euler_DCM, CD_Vessel_Euler_Angular, &
                              CD_Vessel_Point_Kinematics, CD_Vessel_RAO_Set, CD_Vessel_RAO_Eval, &
@@ -49,6 +52,7 @@ PROGRAM test_vessel_motion
   CALL case_deck_rao_regular()
   CALL case_deck_rao_irregular()
   CALL case_deck_rao_spread()
+  CALL case_deck_rao_wavekin7()
   CALL case_deck_rao_stream_rejected()
   CALL case_deck_rejections()
 
@@ -701,6 +705,67 @@ SUBROUTINE case_deck_rao_spread()
     CALL require(err < 1.0e-5_wp, 'rao spread: every spread component drives the vessel with its own phase')
   END SUBROUTINE case_deck_rao_spread
 
+  SUBROUTINE case_deck_rao_wavekin7()
+    !! A heave-only RAO over the full circle of headings: the fairlead below the reference point
+    !! heaves with the elevation of the file's rows there.
+    INTEGER, PARAMETER :: NC = 3
+    REAL(wp), PARAMETER :: W(NC) = [0.6_wp, 0.75_wp, 0.9_wp], RE(NC) = [0.5_wp, -0.2_wp, 0.3_wp]
+    REAL(wp), PARAMETER :: IM(NC) = [0.1_wp, 0.4_wp, -0.25_wp], BETA(NC) = [0.65_wp, 2.9_wp, -2.2_wp]
+    REAL(wp), PARAMETER :: X0 = 31.5_wp, Y0 = -12.0_wp
+    CHARACTER(40) :: opts(9)
+    CHARACTER(512) :: msg
+    LOGICAL :: ok
+    REAL(wp), ALLOCATABLE :: d(:, :)
+    REAL(wp) :: kk(NC), t, eta, err, dist, r, rd, rdd
+    INTEGER :: nr, n, i, es, u
+    CHARACTER(200) :: em
+    OPEN (NEWUNIT=u, FILE='vm_wk7/vm_rao.txt', STATUS='REPLACE', ACTION='WRITE')
+    DO i = 1, 2
+      WRITE (u, '(A,F7.1)') 'HEADING ', MERGE(-180.0_wp, 179.0_wp, i == 1)
+      WRITE (u, '(A)') '0.5 0 0 0 0 1 0 0 0 0 0 0 0'
+      WRITE (u, '(A)') '100 0 0 0 0 1 0 0 0 0 0 0 0'
+    END DO
+    CLOSE (u)
+    OPEN (NEWUNIT=u, FILE='vm_wk7/wave_frequencies.txt', STATUS='REPLACE', ACTION='WRITE')
+    WRITE (u, '(A)') '0.0 0.0 0.0 0.0'
+    DO i = 1, NC
+      WRITE (u, '(4ES24.16)') W(i), RE(i), IM(i), BETA(i)
+    END DO
+    CLOSE (u)
+    opts(1) = '0.05 dtM'
+    opts(2) = '20.0 TMax'
+    opts(3) = '60.0 WtrDpth'
+    opts(4) = '7 WaveKin'
+    opts(5) = 'vm_rao.txt vesselRAO'
+    opts(6) = '90.0|0.0|0.0 vesselRef'
+    opts(7) = '31.5 WaveRefX'
+    opts(8) = '-12.0 WaveRefY'
+    opts(9) = '6.0 rampTime'
+    CALL write_deck('vm_wk7/vm_wk7.dat', opts)
+    CALL run('vm_wk7/vm_wk7.dat', 'vm_wk7/out', ok, msg)
+    CALL require(ok, 'rao WaveKin 7: deck runs: '//TRIM(msg))
+    CALL read_out('vm_wk7/out.out', d, nr, 10)
+    CALL require(nr == 401, 'rao WaveKin 7: full record written')
+    IF (nr /= 401) RETURN
+    DO i = 1, NC
+      CALL CD_Solve_Dispersion_Wavenumber(W(i), 60.0_wp, 9.80665_wp, kk(i), es, em)
+      CALL require(es == 0, 'rao WaveKin 7: dispersion solve: '//TRIM(em))
+    END DO
+    err = 0.0_wp
+    DO n = 1, nr
+      t = d(1, n)
+      CALL CD_Vessel_Ramp(6.0_wp, t, r, rd, rdd)
+      eta = 0.0_wp
+      DO i = 1, NC
+        dist = (90.0_wp - X0)*COS(BETA(i)) + (0.0_wp - Y0)*SIN(BETA(i))
+        eta = eta + SQRT(RE(i)**2 + IM(i)**2)*SIN(W(i)*t - kk(i)*dist + ATAN2(IM(i), RE(i)))
+      END DO
+      err = MAX(err, ABS(d(7, n) - (-14.0_wp + r*eta)))
+    END DO
+    WRITE (*, '(A,ES12.4)') '  rao WaveKin 7: max fairlead heave error [m] = ', err
+    CALL require(err < 1.0e-5_wp, 'rao WaveKin 7: each row drives the vessel with its own heading and phase')
+  END SUBROUTINE case_deck_rao_wavekin7
+
   SUBROUTINE case_deck_rao_stream_rejected()
     CHARACTER(40) :: opts(6)
     CHARACTER(512) :: msg
diff --git a/tests/test_wave_component_sum.f90 b/tests/test_wave_component_sum.f90
index 1992e96..eddb3d3 100644
--- a/tests/test_wave_component_sum.f90
+++ b/tests/test_wave_component_sum.f90
@@ -13,7 +13,11 @@ PROGRAM test_wave_component_sum
   !!      sum of squares: they are independent) and, for the stretched depth, of the elevation;
   !!   2. the sum is a function of its arguments alone: the same call after others gives the same
   !!      bits, and so do the tables in any thread count (the sums keep a fixed order);
-  !!   3. amplitude_scale equals a scaled table to round-off.
+  !!   3. amplitude_scale equals a scaled table to round-off;
+  !!   4. a short-crested table (a direction per component) equals the two long-crested tables it
+  !!      is made of to round-off, for directions in every quadrant, one direction everywhere
+  !!      equals the long-crested sum, and the stretched sum equals the unstretched one at the
+  !!      stretched depth; half a table of cosines fails closed.
   USE, INTRINSIC :: ISO_FORTRAN_ENV, ONLY: INT64
   USE CableDyn_Precision, ONLY: wp
   USE CableDyn_Hydro, ONLY: CD_Component_Wave_Kinematics, CD_Solve_Dispersion_Wavenumber, CD_HYDRO_OK
@@ -45,6 +49,7 @@ PROGRAM test_wave_component_sum
   END DO
   WRITE (*, '(A,F6.3)') '    largest error over its round-off bound: ', worst
   CALL check_pure()
+  CALL check_short_crested()
 
   IF (nfail > 0) THEN
     WRITE (*, '(A,I0,A)') 'FAIL: ', nfail, ' assertion(s) failed'
@@ -263,6 +268,85 @@ SUBROUTINE check_pure()
                  'amplitude_scale equals a scaled table to round-off')
   END SUBROUTINE check_pure
 
+  SUBROUTINE check_short_crested()
+    !! A table past the phase buffer in which every other component goes in direction a or b.
+    INTEGER, PARAMETER :: N = 4200, NH = N/2, NPAIR = 5
+    REAL(wp), PARAMETER :: PAIRS(2, NPAIR) = RESHAPE([0.0_wp, 90.0_wp, 37.5_wp, 127.0_wp, 200.0_wp, 311.3_wp, &
+                                                      -45.0_wp, 233.3_wp, 180.0_wp, 271.0_wp], [2, NPAIR])
+    REAL(wp), PARAMETER :: X = 41.0_wp, Y = -17.5_wp, T = 1234.5_wp, D = 60.0_wp
+    REAL(wp) :: om(N), k(N), amp(N), ph(N), e2kd(N), tail(3, N), cb(N), sb(N), half(4, NH, 2)
+    REAL(wp) :: e0, u0(3), a0(3), p0, e1, u1(3), a1(3), p1, e2, u2(3), a2(3), p2, tol, ze
+    REAL(wp), PARAMETER :: ZS(3) = [-0.8_wp, -7.0_wp, -33.0_wp]
+    INTEGER :: ip, iz, es
+    CHARACTER(200) :: em
+
+    CALL make_table(N, D, 4.0_wp, om, k, amp, ph, e2kd, tail)
+    tol = 32.0_wp*EPSILON(1.0_wp)*SUM(amp*MAX(1.0_wp, om*om, om*om/k))
+    DO iz = 1, 2
+      half(1, :, iz) = om(iz:N:2)
+      half(2, :, iz) = k(iz:N:2)
+      half(3, :, iz) = amp(iz:N:2)
+      half(4, :, iz) = ph(iz:N:2)
+    END DO
+    DO ip = 1, NPAIR
+      cb(1:N:2) = COS(PAIRS(1, ip)*PI/180.0_wp)
+      sb(1:N:2) = SIN(PAIRS(1, ip)*PI/180.0_wp)
+      cb(2:N:2) = COS(PAIRS(2, ip)*PI/180.0_wp)
+      sb(2:N:2) = SIN(PAIRS(2, ip)*PI/180.0_wp)
+      DO iz = 1, SIZE(ZS)
+        CALL CD_Component_Wave_Kinematics(X, Y, ZS(iz), T, D, 0.0_wp, .FALSE., om, k, amp, ph, e0, u0, a0, es, em, &
+                                          pdyn=p0, exp_m2kd=e2kd, tail_sums=tail, comp_cos=cb, comp_sin=sb)
+        CALL require(es == CD_HYDRO_OK, 'short-crested sum: '//TRIM(em))
+        CALL long_crested(X, Y, ZS(iz), T, D, PAIRS(1, ip), half(:, :, 1), e1, u1, a1, p1, es, em)
+        CALL require(es == CD_HYDRO_OK, 'first long-crested half: '//TRIM(em))
+        CALL long_crested(X, Y, ZS(iz), T, D, PAIRS(2, ip), half(:, :, 2), e2, u2, a2, p2, es, em)
+        CALL require(es == CD_HYDRO_OK, 'second long-crested half: '//TRIM(em))
+        CALL require(ABS(e0 - e1 - e2) <= tol .AND. MAXVAL(ABS(u0 - u1 - u2)) <= tol .AND. &
+                     MAXVAL(ABS(a0 - a1 - a2)) <= tol .AND. ABS(p0 - p1 - p2) <= tol, &
+                     'a short-crested table is the sum of its long-crested halves')
+        ! stretched: the sum at the stretched depth of the one total elevation
+        CALL CD_Component_Wave_Kinematics(X, Y, ZS(iz), T, D, 0.0_wp, .TRUE., om, k, amp, ph, e1, u1, a1, es, em, &
+                                          exp_m2kd=e2kd, tail_sums=tail, comp_cos=cb, comp_sin=sb)
+        CALL require(es == CD_HYDRO_OK, 'stretched short-crested sum: '//TRIM(em))
+        IF (ZS(iz) > e1) CYCLE   ! above this trough the point is dry
+        ze = (ZS(iz) - e1)*D/(D + e1)
+        CALL CD_Component_Wave_Kinematics(X, Y, ze, T, D, 0.0_wp, .FALSE., om, k, amp, ph, e2, u2, a2, es, em, &
+                                          exp_m2kd=e2kd, comp_cos=cb, comp_sin=sb)
+        CALL require(es == CD_HYDRO_OK, 'at the stretched depth: '//TRIM(em))
+        CALL require(MAXVAL(ABS(u1 - u2)) <= tol .AND. MAXVAL(ABS(a1 - a2)) <= tol, &
+                     'stretched short-crested sum at the depth of the total elevation')
+      END DO
+      ! one direction everywhere
+      cb = COS(PAIRS(2, ip)*PI/180.0_wp)
+      sb = SIN(PAIRS(2, ip)*PI/180.0_wp)
+      CALL CD_Component_Wave_Kinematics(X, Y, -4.0_wp, T, D, 0.0_wp, .TRUE., om, k, amp, ph, e1, u1, a1, es, em, &
+                                        comp_cos=cb, comp_sin=sb)
+      CALL CD_Component_Wave_Kinematics(X, Y, -4.0_wp, T, D, PAIRS(2, ip), .TRUE., om, k, amp, ph, e2, u2, a2, es, em)
+      CALL require(ABS(e1 - e2) <= tol .AND. MAXVAL(ABS(u1 - u2)) <= tol .AND. MAXVAL(ABS(a1 - a2)) <= tol, &
+                   'one direction per component equals the long-crested sum')
+    END DO
+    CALL CD_Component_Wave_Kinematics(X, Y, -4.0_wp, T, D, 0.0_wp, .TRUE., om, k, amp, ph, e1, u1, a1, es, em, &
+                                      comp_cos=cb)
+    CALL require(es /= CD_HYDRO_OK, 'comp_cos without comp_sin fails closed')
+    CALL CD_Component_Wave_Kinematics(X, Y, -4.0_wp, T, D, 0.0_wp, .TRUE., om, k, amp, ph, e1, u1, a1, es, em, &
+                                      comp_cos=cb(1:NH), comp_sin=sb(1:NH))
+    CALL require(es /= CD_HYDRO_OK, 'direction tables of the wrong length fail closed')
+  END SUBROUTINE check_short_crested
+
+  SUBROUTINE long_crested(x, y, z, t, depth, direction, table, eta, u, a, p, es, em)
+    !! The unstretched sum of the table rows omega, k, amplitude and phase along one direction.
+    REAL(wp), INTENT(IN) :: x, y, z, t, depth, direction, table(:, :)
+    REAL(wp), INTENT(OUT) :: eta, u(3), a(3), p
+    INTEGER, INTENT(OUT) :: es
+    CHARACTER(*), INTENT(OUT) :: em
+    REAL(wp) :: om(SIZE(table, 2)), k(SIZE(table, 2)), amp(SIZE(table, 2)), ph(SIZE(table, 2))
+    om = table(1, :)
+    k = table(2, :)
+    amp = table(3, :)
+    ph = table(4, :)
+    CALL CD_Component_Wave_Kinematics(x, y, z, t, depth, direction, .FALSE., om, k, amp, ph, eta, u, a, es, em, pdyn=p)
+  END SUBROUTINE long_crested
+
   SUBROUTINE sample(p, t, om, k, amp, ph, e2kd, tail, out)
     REAL(wp), INTENT(IN) :: p(3), t, om(:), k(:), amp(:), ph(:), e2kd(:), tail(:, :)
     REAL(wp), INTENT(OUT) :: out(7)
diff --git a/tests/test_wavekin_modes.f90 b/tests/test_wavekin_modes.f90
index 896e9b3..2cfcf46 100644
--- a/tests/test_wavekin_modes.f90
+++ b/tests/test_wavekin_modes.f90
@@ -17,7 +17,13 @@
 !>  - WvHiCOff: 1 rad/s keeps only the record's 8 s harmonic, exactly; none keeps both, as the
 !>    0.5 Hz default does (the record has nothing above it);
 !>  - Currents 1: identical to 1e-12 m/s;
-!>  - the MoorDyn-C modes CableDyn does not support fail closed by name.
+!>  - WaveKin 7 rows of different directions, across the quadrants: the elevation is the sum of
+!>    the two single-direction seas and the velocity, acceleration likewise at the depth stretched
+!>    by the total elevation, to round-off; the beta of the omega = 0 row has no effect;
+!>  - WaveRefX, WaveRefY: an elevation record or a component file referred to a point equals the
+!>    same sea with the deck translated by it, for WaveKin 3 and 7 and WaterKin WaveKinMod 1;
+!>  - the MoorDyn-C modes CableDyn does not support, and a reference point with no record to
+!>    place, fail closed by name.
 PROGRAM test_wavekin_modes
   USE CableDyn_Precision, ONLY: wp
   USE CableDyn_DeckDriver, ONLY: CD_Deck_Fluid_Probe, CD_DECKDRV_OK
@@ -32,6 +38,8 @@ PROGRAM test_wavekin_modes
   CALL check_cutoff()
   CALL check_mode('wavekin_mdc/currents1.dat', 'wavekin_mdc/mdc_currents1.csv', 'Currents 1', 1.0e-12_wp, &
                   1.0e-12_wp, 1.0e-12_wp, absolute=.TRUE.)
+  CALL check_short_crested()
+  CALL check_reference_point()
   CALL check_rejections()
   CALL check_still_water_ramp()
   IF (n_fail > 0) THEN
@@ -234,6 +242,205 @@ SUBROUTINE expect_failure(wavekin, currents, needle)
                  ' fails closed naming "'//needle//'": '//TRIM(em))
   END SUBROUTINE expect_failure
 
+  SUBROUTINE expect_row_failure(wavekin, row, needle)
+    CHARACTER(*), INTENT(IN) :: wavekin, row, needle
+    REAL(wp) :: pos(3, 1), t(1), eta(1), vel(3, 1), acc(3, 1)
+    INTEGER :: es
+    CHARACTER(512) :: em
+    pos = 0.0_wp
+    pos(3, 1) = -5.0_wp
+    t = 0.0_wp
+    CALL write_fail_deck('wavekin_fail.dat', wavekin, '0', row)
+    CALL CD_Deck_Fluid_Probe('wavekin_fail.dat', pos, t, eta, vel, acc, es, em)
+    CALL require(es /= CD_DECKDRV_OK .AND. INDEX(em, needle) > 0, 'WaveKin '//wavekin//' with '//row// &
+                 ' fails closed naming "'//needle//'": '//TRIM(em))
+  END SUBROUTINE expect_row_failure
+
+  SUBROUTINE write_components(beta_a, beta_b, which, beta_mean)
+    !! wave_frequencies.txt of two component sets, set a at direction beta_a and set b at beta_b
+    !! (radians), in file order a then b, so that the frequencies do not ascend; which: 0 both,
+    !! 1 a alone, 2 b alone. beta_mean: the direction of the omega = 0 row (default that of the
+    !! first component).
+    REAL(wp), INTENT(IN) :: beta_a, beta_b
+    INTEGER, INTENT(IN) :: which
+    REAL(wp), INTENT(IN), OPTIONAL :: beta_mean
+    INTEGER :: u
+    REAL(wp) :: b0
+    b0 = MERGE(beta_b, beta_a, which == 2)
+    IF (PRESENT(beta_mean)) b0 = beta_mean
+    OPEN (NEWUNIT=u, FILE='wave_frequencies.txt', STATUS='REPLACE', ACTION='WRITE')
+    WRITE (u, '(A,ES24.16)') '0.0 0.0 0.0 ', b0
+    IF (which /= 2) THEN
+      WRITE (u, '(A,ES24.16)') '0.6 0.7 0.2 ', beta_a
+      WRITE (u, '(A,ES24.16)') '1.1 0.25 -0.3 ', beta_a
+    END IF
+    IF (which /= 1) THEN
+      WRITE (u, '(A,ES24.16)') '0.6 -0.3 0.5 ', beta_b
+      WRITE (u, '(A,ES24.16)') '0.9 0.4 0.1 ', beta_b
+    END IF
+    CLOSE (u)
+  END SUBROUTINE write_components
+
+  SUBROUTINE probe_components(beta_a, beta_b, which, pos, t, eta, vel, acc)
+    REAL(wp), INTENT(IN) :: beta_a, beta_b, pos(:, :), t(:)
+    INTEGER, INTENT(IN) :: which
+    REAL(wp), INTENT(OUT) :: eta(:), vel(:, :), acc(:, :)
+    INTEGER :: es
+    CHARACTER(512) :: em
+    CALL write_components(beta_a, beta_b, which)
+    CALL write_fail_deck('wavekin_short.dat', '7', '0')
+    CALL CD_Deck_Fluid_Probe('wavekin_short.dat', pos, t, eta, vel, acc, es, em)
+    CALL require(es == CD_DECKDRV_OK, 'WaveKin 7 component deck evaluated: '//TRIM(em))
+  END SUBROUTINE probe_components
+
+  SUBROUTINE check_short_crested()
+    !! Rows of different directions against the two single-direction seas they are made of. The
+    !! deck stretches each sea by its own elevation, so a single sea is probed at the depth that
+    !! maps to the depth the total elevation stretches the point to.
+    INTEGER, PARAMETER :: NP = 6, NPAIR = 6
+    REAL(wp), PARAMETER :: D = 50.0_wp
+    ! radians, within +-2 pi: the quadrants, the axes, and directions that are no multiple of 90 deg
+    REAL(wp), PARAMETER :: PAIRS(2, NPAIR) = RESHAPE([0.0_wp, 1.5707963267948966_wp, 0.3_wp, 2.4_wp, -1.0_wp, 3.7_wp, &
+                                                      5.1_wp, -2.2_wp, 3.141592653589793_wp, 4.71238898038469_wp, &
+                                                      0.65_wp, -6.0_wp], [2, NPAIR])
+    REAL(wp) :: pos(3, NP), t(NP), eta(NP), vel(3, NP), acc(3, NP), ea(NP), va(3, NP), aa(3, NP)
+    REAL(wp) :: eb(NP), vb(3, NP), ab(3, NP), ze(NP), pa(3, NP), pb(3, NP), err_e, err_v, err_a, tmp(NP)
+    REAL(wp) :: tmpv(3, NP), tmpa(3, NP), e0(NP), v0(3, NP), a0(3, NP)
+    INTEGER :: ip, es
+    CHARACTER(512) :: em
+    pos(1, :) = [0.0_wp, 13.0_wp, -41.5_wp, 77.0_wp, -8.25_wp, 150.0_wp]
+    pos(2, :) = [0.0_wp, -7.0_wp, 22.0_wp, 90.5_wp, -33.0_wp, -120.0_wp]
+    pos(3, :) = [-2.0_wp, -5.0_wp, -11.0_wp, -20.0_wp, -33.0_wp, -46.0_wp]
+    t = [0.0_wp, 1.7_wp, 3.3_wp, 5.9_wp, 7.1_wp, 9.8_wp]
+    err_e = 0.0_wp
+    err_v = 0.0_wp
+    err_a = 0.0_wp
+    DO ip = 1, NPAIR
+      CALL probe_components(PAIRS(1, ip), PAIRS(2, ip), 0, pos, t, eta, vel, acc)
+      CALL probe_components(PAIRS(1, ip), PAIRS(2, ip), 1, pos, t, ea, tmpv, tmpa)
+      CALL probe_components(PAIRS(1, ip), PAIRS(2, ip), 2, pos, t, eb, tmpv, tmpa)
+      err_e = MAX(err_e, MAXVAL(ABS(eta - ea - eb)))
+      ze = (pos(3, :) - eta)*D/(D + eta)
+      pa = pos
+      pb = pos
+      pa(3, :) = ze*(D + ea)/D + ea
+      pb(3, :) = ze*(D + eb)/D + eb
+      CALL probe_components(PAIRS(1, ip), PAIRS(2, ip), 1, pa, t, tmp, va, aa)
+      CALL probe_components(PAIRS(1, ip), PAIRS(2, ip), 2, pb, t, tmp, vb, ab)
+      err_v = MAX(err_v, MAXVAL(ABS(vel - va - vb)))
+      err_a = MAX(err_a, MAXVAL(ABS(acc - aa - ab)))
+    END DO
+    WRITE (*, '(A,3ES10.3)') 'WaveKin 7 two directions vs two seas: elevation, velocity, acceleration ', err_e, &
+      err_v, err_a
+    CALL require(err_e <= 1.0e-13_wp, 'two directions: the elevation is the sum of the two seas')
+    CALL require(err_v <= 1.0e-13_wp, 'two directions: the velocity is the sum of the two seas')
+    CALL require(err_a <= 1.0e-13_wp, 'two directions: the acceleration is the sum of the two seas')
+    ! the mean row carries no direction: a different beta there leaves a one-direction table
+    CALL write_components(0.4_wp, 0.4_wp, 0)
+    CALL write_fail_deck('wavekin_short.dat', '7', '0')
+    CALL CD_Deck_Fluid_Probe('wavekin_short.dat', pos, t, e0, v0, a0, es, em)
+    CALL require(es == CD_DECKDRV_OK, 'one direction evaluated: '//TRIM(em))
+    CALL write_components(0.4_wp, 0.4_wp, 0, beta_mean=1.3_wp)
+    CALL CD_Deck_Fluid_Probe('wavekin_short.dat', pos, t, eta, vel, acc, es, em)
+    CALL require(es == CD_DECKDRV_OK, 'mean-row direction evaluated: '//TRIM(em))
+    CALL require(ALL(eta == e0) .AND. ALL(vel == v0) .AND. ALL(acc == a0), &
+                 'the direction of the omega = 0 row has no effect')
+  END SUBROUTINE check_short_crested
+
+  SUBROUTINE copy_deck_with(source, target, match, row, replace)
+    !! source with row after the line that holds match, or in its place (replace).
+    CHARACTER(*), INTENT(IN) :: source, target, match, row
+    LOGICAL, INTENT(IN) :: replace
+    CHARACTER(1024) :: line
+    INTEGER :: u_in, u_out, ios
+    OPEN (NEWUNIT=u_in, FILE=source, STATUS='OLD', ACTION='READ')
+    OPEN (NEWUNIT=u_out, FILE=target, STATUS='REPLACE', ACTION='WRITE')
+    DO
+      READ (u_in, '(A)', IOSTAT=ios) line
+      IF (ios /= 0) EXIT
+      IF (INDEX(line, match) > 0) THEN
+        IF (.NOT. replace) WRITE (u_out, '(A)') TRIM(line)
+        WRITE (u_out, '(A)') row
+      ELSE
+        WRITE (u_out, '(A)') TRIM(line)
+      END IF
+    END DO
+    CLOSE (u_in)
+    CLOSE (u_out)
+  END SUBROUTINE copy_deck_with
+
+  SUBROUTINE compare_translated(plain, shifted, shift, label)
+    !! The deck shifted, referred to the point shift, at positions p + shift, against the deck
+    !! plain at p.
+    CHARACTER(*), INTENT(IN) :: plain, shifted, label
+    REAL(wp), INTENT(IN) :: shift(2)
+    INTEGER, PARAMETER :: NP = 5
+    REAL(wp) :: pos(3, NP), moved(3, NP), t(NP), e0(NP), v0(3, NP), a0(3, NP), e1(NP), v1(3, NP), a1(3, NP), err
+    INTEGER :: es
+    CHARACTER(512) :: em
+    pos(1, :) = [0.0_wp, 11.0_wp, -37.0_wp, 60.5_wp, -5.0_wp]
+    pos(2, :) = [0.0_wp, -3.5_wp, 19.0_wp, 44.0_wp, -71.0_wp]
+    pos(3, :) = [-1.5_wp, -6.0_wp, -14.0_wp, -27.0_wp, -41.0_wp]
+    t = [0.3_wp, 1.9_wp, 4.4_wp, 6.2_wp, 9.1_wp]
+    moved = pos
+    moved(1, :) = pos(1, :) + shift(1)
+    moved(2, :) = pos(2, :) + shift(2)
+    CALL CD_Deck_Fluid_Probe(plain, pos, t, e0, v0, a0, es, em)
+    CALL require(es == CD_DECKDRV_OK, label//' plain deck evaluated: '//TRIM(em))
+    IF (es /= CD_DECKDRV_OK) RETURN
+    CALL CD_Deck_Fluid_Probe(shifted, moved, t, e1, v1, a1, es, em)
+    CALL require(es == CD_DECKDRV_OK, label//' referred deck evaluated: '//TRIM(em))
+    IF (es /= CD_DECKDRV_OK) RETURN
+    err = MAX(MAXVAL(ABS(e1 - e0)), MAXVAL(ABS(v1 - v0)), MAXVAL(ABS(a1 - a0)))
+    WRITE (*, '(A,A,ES10.3)') label, ' referred to a point vs the deck translated: ', err
+    CALL require(err <= 1.0e-12_wp, label//': a shifted reference point equals the deck translated')
+    ! and the point matters: the same positions without the shift are another sea
+    CALL CD_Deck_Fluid_Probe(shifted, pos, t, e1, v1, a1, es, em)
+    CALL require(MAXVAL(ABS(e1 - e0)) > 1.0e-3_wp, label//': the reference point moves the sea')
+  END SUBROUTINE compare_translated
+
+  SUBROUTINE check_reference_point()
+    REAL(wp), PARAMETER :: SHIFT(2) = [250.37_wp, -70.19_wp]
+    CHARACTER(*), PARAMETER :: ROWS = '250.37 WaveRefX'//NEW_LINE('a')//'-70.19 WaveRefY'
+    INTEGER :: u
+    ! WaveKin 7, one direction and two
+    CALL write_components(0.9_wp, 0.9_wp, 0)
+    CALL write_fail_deck('wavekin_ref0.dat', '7', '0')
+    CALL write_fail_deck('wavekin_ref1.dat', '7', '0', ROWS)
+    CALL compare_translated('wavekin_ref0.dat', 'wavekin_ref1.dat', SHIFT, 'WaveKin 7')
+    CALL write_components(0.65_wp, 2.9_wp, 0)
+    CALL compare_translated('wavekin_ref0.dat', 'wavekin_ref1.dat', SHIFT, 'WaveKin 7, two directions')
+    ! WaveKin 3
+    CALL copy_deck_with('wavekin_mdc/wavekin3.dat', 'wavekin_mdc/wavekin3_ref.dat', 'dtWave', ROWS, .FALSE.)
+    CALL compare_translated('wavekin_mdc/wavekin3.dat', 'wavekin_mdc/wavekin3_ref.dat', SHIFT, 'WaveKin 3')
+    ! WaterKin WaveKinMod 1, along a direction off the axes
+    OPEN (NEWUNIT=u, FILE='wavekin_mdc/waterkin_ref.dat', STATUS='REPLACE', ACTION='WRITE')
+    WRITE (u, '(A)') 'MoorDyn v2 water kinematics file'
+    WRITE (u, '(A)') '(header line 2)'
+    WRITE (u, '(A)') '--------------------------- WAVES -------------------------------------'
+    WRITE (u, '(A)') '1                    WaveKinMod'
+    WRITE (u, '(A)') '"wave_elevation.txt" WaveKinFile'
+    WRITE (u, '(A)') '0.1                  dtWave'
+    WRITE (u, '(A)') '37.0                 WaveDir'
+    WRITE (u, '(A)') '0                    X wave input'
+    WRITE (u, '(A)') '0.0'
+    WRITE (u, '(A)') '0                    Y wave input'
+    WRITE (u, '(A)') '0.0'
+    WRITE (u, '(A)') '0                    Z wave input'
+    WRITE (u, '(A)') '0.0'
+    WRITE (u, '(A)') '--------------------------- CURRENT -----------------------------------'
+    WRITE (u, '(A)') '0                    CurrentMod'
+    WRITE (u, '(A)') '0                    Z current'
+    WRITE (u, '(A)') '0.0'
+    WRITE (u, '(A)') '--------------------- need this line ----------------------------------'
+    CLOSE (u)
+    CALL copy_deck_with('wavekin_mdc/wavekin3.dat', 'wavekin_mdc/waterkin0.dat', 'WaveKin ', &
+                        'waterkin_ref.dat WaterKin', .TRUE.)
+    CALL copy_deck_with('wavekin_mdc/waterkin0.dat', 'wavekin_mdc/waterkin1.dat', 'WaterKin', &
+                        'waterkin_ref.dat WaterKin'//NEW_LINE('a')//ROWS, .TRUE.)
+    CALL compare_translated('wavekin_mdc/waterkin0.dat', 'wavekin_mdc/waterkin1.dat', SHIFT, 'WaterKin WaveKinMod 1')
+  END SUBROUTINE check_reference_point
+
   SUBROUTINE check_still_water_ramp()
     !! rampTime is valid without waves (it also ramps prescribed motion). The still-water
     !! kinematics must then be exactly zero, including at t = 0 where the ramp is closed
@@ -252,19 +459,17 @@ SUBROUTINE check_still_water_ramp()
   END SUBROUTINE check_still_water_ramp
 
   SUBROUTINE check_rejections()
-    INTEGER :: u
     CALL expect_failure('1', '0', 'MoorDyn-C API')
     CALL expect_failure('2', '0', 'FFT_GRID')
     CALL expect_failure('4', '0', 'no kinematics source in MoorDyn-C')
     CALL expect_failure('6', '0', 'no kinematics source in MoorDyn-C')
     CALL expect_failure('0', '3', 'no current source in MoorDyn-C')
     CALL expect_failure('0', '4', 'no current source in MoorDyn-C')
-    ! short-crested components (rows with different directions)
-    OPEN (NEWUNIT=u, FILE='wave_frequencies.txt', STATUS='REPLACE', ACTION='WRITE')
-    WRITE (u, '(A)') '0.0 0.0 0.0 0.0'
-    WRITE (u, '(A)') '0.6 0.7 0.2 0.0'
-    WRITE (u, '(A)') '0.9 -0.3 0.5 1.57'
-    CLOSE (u)
-    CALL expect_failure('7', '0', 'different directions')
+    ! a reference point needs a record or component file to place, and a bounded position
+    CALL write_components(0.4_wp, 0.4_wp, 0)
+    CALL expect_row_failure('0', '1.0 WaveRefX', 'reference point')
+    CALL expect_row_failure('0', '1.0 WaveRefY'//NEW_LINE('a')//'airy 2.0 8.0 0.0 waves', 'reference point')
+    CALL expect_row_failure('7', '2.0e6 WaveRefX', 'at most 1e6 m')
+    CALL expect_row_failure('7', 'nan WaveRefY', 'WaveRefY')
   END SUBROUTINE check_rejections
 END PROGRAM test_wavekin_modes