diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml
index 4631571..01f714c 100644
--- a/.github/workflows/CI.yml
+++ b/.github/workflows/CI.yml
@@ -62,6 +62,7 @@ jobs:
           - r3x
           - bc
           - multipatch
+          - nesting_support
           - tile_multipatch
           - basis_cache
           - relocation
diff --git a/CHANGELOG.md b/CHANGELOG.md
index 04b529e..82700b9 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -9,6 +9,37 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 ### Added
 
+- **Temporal-nesting support for two-way nested models** (DeMaria et al. 1992;
+  Ooyama 2001):
+  - **`lerp_payload!(dest, p0, p1, θ)`** linearly interpolates two
+    [`InterfacePayload`]s in time, bitwise-exact at the endpoints. A subcycling
+    child patch applies the interpolated parent trio at each of its substeps
+    between the parent's bracketing steps.
+  - **`evaluate_grid_ipoints(grid, xq)`** (and the in-place `!` form) evaluates a
+    grid's spectral representation at arbitrary i-direction points with the same
+    variable/derivative-slice layout as `grid.physical` (1-D spline `R` and
+    `RiRk` grids). This is the fine→coarse feedback primitive: the coarse patch's
+    collar quadrature points are evaluated on the fine grid and injected into the
+    coarse Galerkin loads, per DeMaria et al. (1992) eq. 2.22 — feedback through
+    the tendencies, not through a boundary condition.
+  - **Collar interfaces** — a parent patch extended one cell past the nominal
+    junction into the child's domain, so the trio the child's R3X boundary reads
+    consists of interior, freely-fitted parent amplitudes — are covered by the
+    existing `PatchInterface(...; is_stacked=true)` path; the new
+    `test/nesting_support.jl` suite documents the construction and adds an
+    anti-freeze regression against BC-based (dual-R3X) feedback, which is
+    degenerate (a rank-3 spline's border trio is identically its `ahat`).
+  - **`evaluate_grid_points(grid::RL_Grid, pts; kmax)`** — values plus
+    ∂r/∂²r/∂λ/∂²λ at arbitrary `(r, λ)` points for radially-nested RL models,
+    with an optional per-point wavenumber truncation (the azimuthal analogue of
+    ring transmissibility: injected values must not exceed the target ring's
+    supported wavenumbers). Registry-aware via `_get_ahat_cache_rl`.
+  - **Nest-annulus grids**: an RL patch may set `patchOffsetL` explicitly (with
+    `spectralIndexL` left patch-relative) so its ring point counts and
+    wavenumber support follow the global ring numbering, and
+    `_create_tile_from_patch` now composes the patch's own `patchOffsetL` into
+    its tiles (identity for ordinary patches).
+
 - **`Thermodynamics.potential_temperature(p_Pa, rho_d)`** — a new two-argument method giving the
   dry potential temperature from the `(p, ρ_d)` pair, `θ_d ≡ (p₀^κ/R_d)·p^(1−κ)/ρ_d`. It is the
   same quantity as the existing `potential_temperature(s, rho_d, q_v)`, but needs no entropy
@@ -84,6 +115,32 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 ### Fixed
 
+- **Unstructured RL/RLZ evaluation flipped every sine (odd-in-λ) component.**
+  `_eval_unstructured_rl`/`_eval_unstructured_rlz` (behind `evaluate_unstructured`,
+  `interpolate_to_grid`, and grid relocation) synthesized the Fourier series with
+  `+ aI·sin(kλ)`, but the FFTW halfcomplex "imag" slots hold the **negative**
+  sine sums — so any azimuthally-asymmetric field was evaluated as its mirror
+  image. The synthesis now subtracts the sine term, matching the `HC2R` inverse
+  used by `FItransform!`; a new mish-point test pins the convention against
+  `gridTransform!` (the existing interpolation/relocation tests passed under
+  either sign and never caught this). **This changes results** for any consumer
+  that evaluated asymmetric RL/RLZ fields at unstructured points.
+- **The tiled RL b→a solves now reload the coupled per-wavenumber border.**
+  The 2- and 3-argument `splineTransform!(…, ::RL_Grid)` reused the three
+  k0/real/imag spline objects across all wavenumbers without reloading their
+  `ahat` from the multi-patch registry, so an R3X-coupled RL patch carried the
+  *last-applied* wavenumber's border on every mode (only k=0 was right). Both
+  methods — and the `_get_ahat_cache_rl` coefficient cache behind the
+  unstructured evaluators — now mirror `gridTransform`'s per-wavenumber registry
+  loads.
+- **The out-of-place `SAtransform(spline, b)` now honors the R3X `ahat`.** This
+  allocating form is what every 3-argument (tiled) `splineTransform!` method uses
+  for the b→a solve on the patch splines, and it ignored `spline.ahat` — so a
+  rank-3-coupled (R3X) patch lost its donated border trio on the distributed tile
+  path, with the border silently pinned to zero while the in-place `SAtransform!`
+  honored the coupling. The rank-3 inhomogeneous path now matches `SAtransform!`
+  exactly; a regression test asserts allocating == in-place for an R3X spline
+  with nonzero `ahat`.
 - **`grid_from_regular_data` (and `grid_from_netcdf`, which forwards to it) never accepted
   `BoundaryConditions` structs.** Its BC keyword arguments were annotated `::Dict`, so a bare
   `DirichletBC()` — which is not a `Dict` — raised a `MethodError`. The struct-BC API simply
diff --git a/src/CubicBSpline.jl b/src/CubicBSpline.jl
index 4294145..c1e0a6d 100644
--- a/src/CubicBSpline.jl
+++ b/src/CubicBSpline.jl
@@ -1354,6 +1354,14 @@ end
 
 function SAtransform(spline::Spline1D, b::AbstractVector)
 
+    # R3X (inhomogeneous rank-3) splines carry coupled border coefficients in
+    # spline.ahat; honor them exactly as the in-place SAtransform! does. The
+    # tiled b→a path (3-arg splineTransform!) relies on this so that patch
+    # splines coupled to a nest neighbor reconstruct with the donated border.
+    if _has_r3x(spline.params)
+        btilde = spline.gammaBC * (b .- (spline.pq * spline.ahat))
+        return (spline.gammaBC' * (spline.pqFactor \ btilde)) .+ spline.ahat
+    end
     a = spline.gammaBC' * (spline.pqFactor \ (spline.gammaBC * b))
     return a
 end
diff --git a/src/Springsteel.jl b/src/Springsteel.jl
index 171843e..2ae23c4 100644
--- a/src/Springsteel.jl
+++ b/src/Springsteel.jl
@@ -285,6 +285,7 @@ export compute_interface_payload, compute_interface_payload!
 export apply_interface_payload!
 export PatchChain, PatchEmbedded
 export AbstractMultiGrid, SpringsteelMultiGrid, createMultiGrid
+export lerp_payload!, evaluate_grid_ipoints, evaluate_grid_ipoints!, evaluate_grid_points
 
 # Boundary condition type system
 export BoundaryConditions, BCSpec, bc_rank, is_periodic, is_inhomogeneous
diff --git a/src/interpolation.jl b/src/interpolation.jl
index af98d21..e82949f 100644
--- a/src/interpolation.jl
+++ b/src/interpolation.jl
@@ -1567,9 +1567,19 @@ function _get_ahat_cache_rl(source::SpringsteelGrid, sv::Int)
         return entry.a_coeffs
     end
 
+    # Multi-patch coupling: reload the per-wavenumber border coefficients into
+    # the reused splines before each solve (see gridTransform(_RLGrid)). The
+    # cache key hashes the spectral b only; that is sufficient in practice
+    # because a coupled patch's b changes whenever its neighbors' data does
+    # (every step of a time integration).
+    has_wn_ahat = _has_wavenumber_ahat(source)
+
     a_cache = zeros(Float64, b_iDim, n_kslots)
     sp0 = source.ibasis.data[1, sv]
     sp0.b .= view(source.spectral, 1:b_iDim, sv)
+    if has_wn_ahat
+        sp0.ahat .= _get_wavenumber_ahat(source, sv, 0)
+    end
     SAtransform!(sp0)
     a_cache[:, 1] .= sp0.a
 
@@ -1580,11 +1590,17 @@ function _get_ahat_cache_rl(source::SpringsteelGrid, sv::Int)
 
         spc = source.ibasis.data[2, sv]
         spc.b .= view(source.spectral, p1c:p2c, sv)
+        if has_wn_ahat
+            spc.ahat .= _get_wavenumber_ahat(source, sv, p)
+        end
         SAtransform!(spc)
         a_cache[:, 2k] .= spc.a
 
         sps = source.ibasis.data[3, sv]
         sps.b .= view(source.spectral, p1s:p2s, sv)
+        if has_wn_ahat
+            sps.ahat .= _get_wavenumber_ahat(source, sv, p + 1)
+        end
         SAtransform!(sps)
         a_cache[:, 2k + 1] .= sps.a
     end
@@ -1732,8 +1748,12 @@ function _eval_unstructured_rl(source::SpringsteelGrid, pts::AbstractMatrix{Floa
     for n in 1:npts
         f = sc.ak[n, 1]
         λ = pts[n, 2]
+        # FFTW halfcomplex convention: the "imag" slot holds the NEGATIVE sine
+        # sum, so synthesis is a0 + Σ 2(aR cos kλ − aI sin kλ). The previous
+        # `+ sin` sign silently flipped every sine (odd-in-λ) component of the
+        # evaluated field.
         for k in 1:kDim
-            f += 2.0 * (sc.ak[n, 2k] * cos(k * λ) + sc.ak[n, 2k + 1] * sin(k * λ))
+            f += 2.0 * (sc.ak[n, 2k] * cos(k * λ) - sc.ak[n, 2k + 1] * sin(k * λ))
         end
         result[n] = f
     end
@@ -1783,8 +1803,10 @@ function _eval_unstructured_rlz(source::SpringsteelGrid, pts::AbstractMatrix{Flo
 
         for z_b in 1:b_kDim
             val = sc.spline_vals[n, z_b, 1]
+            # FFTW halfcomplex: "imag" slots hold the NEGATIVE sine sum (see
+            # _eval_unstructured_rl) — synthesis subtracts the sine term.
             for k in 1:kDim_wn
-                val += 2.0 * (sc.spline_vals[n, z_b, 2k] * cos(k * λ) +
+                val += 2.0 * (sc.spline_vals[n, z_b, 2k] * cos(k * λ) -
                               sc.spline_vals[n, z_b, 2k + 1] * sin(k * λ))
             end
             cheb_col.b[z_b] = val
diff --git a/src/multipatch.jl b/src/multipatch.jl
index 287bc4e..da0480e 100644
--- a/src/multipatch.jl
+++ b/src/multipatch.jl
@@ -1573,3 +1573,273 @@ function createMultiGrid(config::Dict{Symbol, Any};
     end
     return SpringsteelMultiGrid(cfg, mpg)
 end
+
+# ────────────────────────────────────────────────────────────────────────────
+# Temporal nesting support: payload time-interpolation + collar evaluation
+# ────────────────────────────────────────────────────────────────────────────
+#
+# Two-way nesting à la DeMaria et al. (1992, MWR) / Ooyama (2001, JAS):
+#
+# - Coarse→fine: the fine patch's R3X border trio comes from the coarse
+#   amplitudes at the interface node (COUPLING_MATRIX_2X).  When the fine
+#   patch subcycles with a smaller timestep, the coarse trio is linearly
+#   interpolated in time between the bracketing coarse steps
+#   (`lerp_payload!`).
+# - Fine→coarse: NOT a boundary condition.  The coarse patch extends one
+#   coarse cell past the interface into the fine domain (the "collar"); its
+#   Galerkin load integrals there are evaluated from the fine-mesh
+#   representation.  `evaluate_grid_ipoints` provides those fine-mesh values
+#   (and derivative slices, matching the `physical` array layout) at the
+#   coarse patch's collar quadrature points.
+#
+# A collar interface is constructed with the existing
+# `PatchInterface(...; is_stacked=true)` path: the primary (coarse) domain
+# extends one cell past the secondary (fine) boundary, so the extraction trio
+# sits at interior, freely-fitted primary nodes.
+
+"""
+    lerp_payload!(dest, p0, p1, θ) -> InterfacePayload
+
+Linearly interpolate two [`InterfacePayload`](@ref)s in time:
+`dest.border = (1-θ)·p0.border + θ·p1.border`.
+
+Used by subcycling drivers: the primary computes payloads at its step
+endpoints `t_n` (`p0`) and `t_{n+1}` (`p1`); each secondary substep at
+fraction `θ ∈ [0,1]` of the primary step applies the interpolated payload.
+`θ == 0.0` and `θ == 1.0` reproduce `p0`/`p1` bitwise.
+
+All three payloads must agree in scheme, side, nvars, and n_slots.
+"""
+function lerp_payload!(dest::InterfacePayload, p0::InterfacePayload,
+                       p1::InterfacePayload, θ::Float64)
+    for p in (p0, p1)
+        p.scheme === dest.scheme || throw(ArgumentError(
+            "Payload scheme :$(p.scheme) ≠ dest scheme :$(dest.scheme)"))
+        p.side === dest.side || throw(ArgumentError(
+            "Payload side :$(p.side) ≠ dest side :$(dest.side)"))
+        p.nvars == dest.nvars || throw(ArgumentError(
+            "Payload nvars $(p.nvars) ≠ dest nvars $(dest.nvars)"))
+        p.n_slots == dest.n_slots || throw(ArgumentError(
+            "Payload n_slots $(p.n_slots) ≠ dest n_slots $(dest.n_slots)"))
+    end
+    (0.0 <= θ <= 1.0) || throw(ArgumentError("θ must be in [0,1], got $θ"))
+    if θ == 0.0
+        copyto!(dest.border, p0.border)
+    elseif θ == 1.0
+        copyto!(dest.border, p1.border)
+    else
+        @. dest.border = (1.0 - θ) * p0.border + θ * p1.border
+    end
+    return dest
+end
+
+"""
+    evaluate_grid_ipoints(grid, xq) -> Array{Float64,3}
+
+Evaluate `grid`'s spectral representation at arbitrary i-direction points
+`xq` (which must lie inside the grid's i-domain), returning an array with the
+same variable and derivative-slice layout as `grid.physical`, but with the
+i-mish points replaced by `xq`.
+
+- 1-D spline grid (`R`): output is `(length(xq), nvars, 3)` with slices
+  1=value, 2=∂i, 3=∂²i.
+- 2-D spline×spline grid (`RiRk`): output is `(length(xq)·kDim, nvars, 5)`
+  with rows ordered `(q-1)·kDim + z` (i-outer, k-inner, matching `physical`)
+  and slices 1=value, 2=∂i, 3=∂²i, 4=∂k, 5=∂²k.  The k-direction is
+  evaluated at the grid's own vertical mish points.
+
+This is the fine→coarse feedback primitive for two-way nesting: the coarse
+patch's collar quadrature points are evaluated on the fine grid, and the
+result overwrites the coarse `physical` rows there before the coarse
+tendency computation (DeMaria et al. 1992, eq. 2.22 discussion).
+
+Reads `grid.spectral`; mutates the grid's basis spline scratch (`.b`, `.a`)
+exactly as `gridTransform!` does, so it must not run concurrently with other
+transforms on the same grid.
+"""
+function evaluate_grid_ipoints(grid::_1DCartesianGrid, xq::AbstractVector{Float64})
+    nvars = size(grid.spectral, 2)
+    out = zeros(Float64, length(xq), nvars, 3)
+    return evaluate_grid_ipoints!(out, grid, xq)
+end
+
+function evaluate_grid_ipoints!(out::AbstractArray{Float64,3},
+                                grid::_1DCartesianGrid,
+                                xq::AbstractVector{Float64})
+    nvars = size(grid.spectral, 2)
+    size(out) == (length(xq), nvars, 3) || throw(DimensionMismatch(
+        "out must be ($(length(xq)), $nvars, 3), got $(size(out))"))
+    _check_ipoints_in_domain(grid, xq)
+    for v in 1:nvars
+        spline = grid.ibasis.data[1, v]
+        spline.b .= view(grid.spectral, :, v)
+        CubicBSpline.SAtransform!(spline)
+        for dr in 0:2
+            CubicBSpline.SItransform(spline.params, spline.a, xq,
+                                     view(out, :, v, dr + 1), dr)
+        end
+    end
+    return out
+end
+
+function evaluate_grid_ipoints(grid::_2DCartesianRiRk, xq::AbstractVector{Float64})
+    kDim  = grid.params.kDim
+    nvars = size(grid.spectral, 2)
+    out = zeros(Float64, length(xq) * kDim, nvars, 5)
+    return evaluate_grid_ipoints!(out, grid, xq)
+end
+
+function evaluate_grid_ipoints!(out::AbstractArray{Float64,3},
+                                grid::_2DCartesianRiRk,
+                                xq::AbstractVector{Float64})
+    nq     = length(xq)
+    kDim   = grid.params.kDim
+    b_iDim = grid.params.b_iDim
+    b_kDim = grid.params.b_kDim
+    nvars  = size(grid.spectral, 2)
+    size(out) == (nq * kDim, nvars, 5) || throw(DimensionMismatch(
+        "out must be ($(nq * kDim), $nvars, 5), got $(size(out))"))
+    _check_ipoints_in_domain(grid, xq)
+
+    buf = Matrix{Float64}(undef, nq, b_kDim)   # i-evaluated k-mode coefficients
+    kscratch = Vector{Float64}(undef, kDim)
+
+    for v in 1:nvars
+        kcol = grid.kbasis.data[v]
+        for dr in 0:2
+            # i-direction: evaluate each k-spectral mode's i-spline at xq
+            for z in 1:b_kDim
+                r1 = (z - 1) * b_iDim + 1
+                isp = grid.ibasis.data[z, v]
+                copyto!(isp.b, 1, grid.spectral, r1 + (v - 1) * size(grid.spectral, 1), b_iDim)
+                CubicBSpline.SAtransform!(isp)
+                CubicBSpline.SItransform(isp.params, isp.a, xq,
+                                         view(buf, :, z), dr)
+            end
+            # k-direction: inverse transform per evaluation point
+            for q in 1:nq
+                @inbounds for z in 1:b_kDim
+                    kcol.b[z] = buf[q, z]
+                end
+                CubicBSpline.SAtransform!(kcol)
+                z1 = (q - 1) * kDim + 1
+                z2 = z1 + kDim - 1
+                if dr == 0
+                    CubicBSpline.SItransform!(kcol)
+                    copyto!(view(out, z1:z2, v, 1), kcol.uMish)
+                    CubicBSpline.SIxtransform(kcol, kscratch)
+                    copyto!(view(out, z1:z2, v, 4), kscratch)
+                    CubicBSpline.SIxxtransform(kcol, kscratch)
+                    copyto!(view(out, z1:z2, v, 5), kscratch)
+                elseif dr == 1
+                    CubicBSpline.SItransform!(kcol)
+                    copyto!(view(out, z1:z2, v, 2), kcol.uMish)
+                else
+                    CubicBSpline.SItransform!(kcol)
+                    copyto!(view(out, z1:z2, v, 3), kcol.uMish)
+                end
+            end
+        end
+    end
+    return out
+end
+
+function _check_ipoints_in_domain(grid::SpringsteelGrid, xq::AbstractVector{Float64})
+    lo, hi = _get_domain_bounds(grid, :i)
+    tol = 1e-9 * max(abs(lo), abs(hi), 1.0)
+    for x in xq
+        (lo - tol <= x <= hi + tol) || throw(ArgumentError(
+            "Evaluation point $x outside grid i-domain [$lo, $hi]"))
+    end
+    return nothing
+end
+
+"""
+    evaluate_grid_points(grid::RL_Grid, pts; kmax=nothing) -> Array{Float64,3}
+
+Evaluate an RL (spline-radius × Fourier-azimuth) grid's spectral
+representation at arbitrary `(r, λ)` points — `pts` is `(N, 2)` with radius in
+column 1 and azimuth in column 2 — returning `(N, nvars, 5)` with the RL
+`physical` slice layout: 1=value, 2=∂r, 3=∂²r, 4=∂λ, 5=∂²λ.
+
+`kmax` (optional, `Vector{Int}` of length N) truncates the wavenumber sum per
+point. RL rings only support wavenumbers up to their global ring index, so a
+consumer injecting these values into another RL grid's rings should pass the
+TARGET ring's supported maximum — the azimuthal analogue of the radial
+transmissibility truncation (higher-k content would alias on the target
+ring's coarser azimuth points).
+
+This is the fine→coarse collar-feedback primitive for radially-nested RL
+models: the coarse annulus's collar ring points are evaluated on the fine
+patch and injected into the coarse Galerkin loads. Reads `grid.spectral` (b)
+and honors the per-wavenumber coupled border registry via
+`_get_ahat_cache_rl`. The λ-derivatives are analytic (FFTW halfcomplex
+convention: "imag" slots are negative sine sums).
+"""
+function evaluate_grid_points(grid::_RLGrid, pts::AbstractMatrix{Float64};
+                              kmax::Union{Nothing, Vector{Int}} = nothing)
+    gp = grid.params
+    b_iDim = gp.b_iDim
+    kDim = gp.iDim + gp.patchOffsetL
+    npts = size(pts, 1)
+    kmax === nothing || length(kmax) == npts || throw(ArgumentError(
+        "kmax must have one entry per point ($(npts)), got $(length(kmax))"))
+    size(pts, 2) == 2 || throw(ArgumentError("pts must be (N, 2) [r λ], got $(size(pts))"))
+    n_kslots = 1 + 2 * kDim
+    nvars = length(gp.vars)
+
+    lo, hi = gp.iMin, gp.iMax
+    tol = 1e-9 * max(abs(lo), abs(hi), 1.0)
+    for n in 1:npts
+        (lo - tol <= pts[n, 1] <= hi + tol) || throw(ArgumentError(
+            "Evaluation radius $(pts[n, 1]) outside grid domain [$lo, $hi]"))
+    end
+
+    out = zeros(Float64, npts, nvars, 5)
+    ak0 = Matrix{Float64}(undef, npts, n_kslots)   # radial value
+    ak1 = Matrix{Float64}(undef, npts, n_kslots)   # ∂r
+    ak2 = Matrix{Float64}(undef, npts, n_kslots)   # ∂²r
+    r_view = view(pts, :, 1)
+
+    for v in 1:nvars
+        a_cache = _get_ahat_cache_rl(grid, v)
+        sp0 = grid.ibasis.data[1, v]
+        for slot in 1:n_kslots
+            a_slice = view(a_cache, :, slot)
+            CubicBSpline.SItransform(sp0.params, a_slice, r_view,
+                                     view(ak0, :, slot), 0)
+            CubicBSpline.SItransform(sp0.params, a_slice, r_view,
+                                     view(ak1, :, slot), 1)
+            CubicBSpline.SItransform(sp0.params, a_slice, r_view,
+                                     view(ak2, :, slot), 2)
+        end
+        @inbounds for n in 1:npts
+            λ = pts[n, 2]
+            f   = ak0[n, 1]
+            fr  = ak1[n, 1]
+            frr = ak2[n, 1]
+            fl  = 0.0
+            fll = 0.0
+            # FFTW halfcomplex convention: the "imag" slot holds the NEGATIVE
+            # sine sum, so synthesis is a0 + Σ 2(aR cos kλ − aI sin kλ)
+            # (matching the HC2R inverse used by FItransform!).
+            kn = kmax === nothing ? kDim : min(kDim, kmax[n])
+            for k in 1:kn
+                ck = cos(k * λ)
+                sk = sin(k * λ)
+                aR0 = ak0[n, 2k]; aI0 = ak0[n, 2k + 1]
+                f   += 2.0 * (aR0 * ck - aI0 * sk)
+                fr  += 2.0 * (ak1[n, 2k] * ck - ak1[n, 2k + 1] * sk)
+                frr += 2.0 * (ak2[n, 2k] * ck - ak2[n, 2k + 1] * sk)
+                fl  += 2.0 * k * (-aR0 * sk - aI0 * ck)
+                fll -= 2.0 * k * k * (aR0 * ck - aI0 * sk)
+            end
+            out[n, v, 1] = f
+            out[n, v, 2] = fr
+            out[n, v, 3] = frr
+            out[n, v, 4] = fl
+            out[n, v, 5] = fll
+        end
+    end
+    return out
+end
diff --git a/src/tiling.jl b/src/tiling.jl
index 23c8f83..89cd260 100644
--- a/src/tiling.jl
+++ b/src/tiling.jl
@@ -82,6 +82,14 @@ function _create_tile_from_patch(patch::SpringsteelGrid,
         mubar          = patch.params.mubar,
         quadrature     = patch.params.quadrature,
         spectralIndexL = spectralIndexL,
+        # Compose the patch's own physical offset (nonzero for nest annulus
+        # patches whose ring numbering / wavenumber support follows the GLOBAL
+        # grid) with the tile's position inside the patch. spectralIndexL stays
+        # patch-relative — it windows the tile into the patch spectral array —
+        # while patchOffsetL carries the global mish offset. For ordinary
+        # patches (patchOffsetL == 0) this reduces to the previous derivation.
+        patchOffsetL   = patch.params.patchOffsetL +
+                         (spectralIndexL - 1) * patch.params.mubar,
         tile_num       = tile_num)
     return createGrid(tile_gp)
 end
@@ -1514,10 +1522,17 @@ function splineTransform!(sharedSpectral::SharedArray{real}, tile::_RLGrid)
     kDim   = tile.params.iDim + tile.params.patchOffsetL
     nvars  = length(tile.params.vars)
 
+    # Reload the coupled per-wavenumber border coefficients before every solve
+    # on the reused splines (see the 3-arg method and gridTransform).
+    has_wn_ahat = _has_wavenumber_ahat(tile)
+
     for v in 1:nvars
         # k = 0 block
         k1 = 1
         k2 = b_iDim
+        if has_wn_ahat
+            tile.ibasis.data[1, v].ahat .= _get_wavenumber_ahat(tile, v, 0)
+        end
         tile.spectral[k1:k2, v] .= SAtransform(tile.ibasis.data[1, v],
                                                  view(sharedSpectral, k1:k2, v))
 
@@ -1526,10 +1541,16 @@ function splineTransform!(sharedSpectral::SharedArray{real}, tile::_RLGrid)
             p  = k * 2
             p1 = ((p - 1) * b_iDim) + 1
             p2 = p * b_iDim
+            if has_wn_ahat
+                tile.ibasis.data[2, v].ahat .= _get_wavenumber_ahat(tile, v, p)
+            end
             tile.spectral[p1:p2, v] .= SAtransform(tile.ibasis.data[2, v],
                                                      view(sharedSpectral, p1:p2, v))
             p1 = p * b_iDim + 1
             p2 = (p + 1) * b_iDim
+            if has_wn_ahat
+                tile.ibasis.data[3, v].ahat .= _get_wavenumber_ahat(tile, v, p + 1)
+            end
             tile.spectral[p1:p2, v] .= SAtransform(tile.ibasis.data[3, v],
                                                      view(sharedSpectral, p1:p2, v))
         end
@@ -3511,8 +3532,18 @@ function splineTransform!(sharedSpectral::SharedArray{real},
     kDim         = tile.params.iDim + tile.params.patchOffsetL
     nvars        = length(tile.params.vars)
 
+    # Multi-patch coupling: the reused splines (k0/real/imag) are shared across
+    # wavenumbers, so the coupled border coefficients must be reloaded from the
+    # per-wavenumber registry before EVERY SAtransform, exactly as the non-tiled
+    # gridTransform does — otherwise every k >= 1 solve reuses whatever ahat the
+    # last apply_interface_payload! left behind.
+    has_wn_ahat = _has_wavenumber_ahat(patch)
+
     for v in 1:nvars
         # k = 0 block
+        if has_wn_ahat
+            patch.ibasis.data[1, v].ahat .= _get_wavenumber_ahat(patch, v, 0)
+        end
         patch.spectral[1:b_iDim_patch, v] .=
             SAtransform(patch.ibasis.data[1, v],
                         view(sharedSpectral, 1:b_iDim_patch, v))
@@ -3526,6 +3557,9 @@ function splineTransform!(sharedSpectral::SharedArray{real},
             # Real part
             pp1 = (p - 1) * b_iDim_patch + 1
             tp1 = (p - 1) * b_iDim_tile + 1
+            if has_wn_ahat
+                patch.ibasis.data[2, v].ahat .= _get_wavenumber_ahat(patch, v, p)
+            end
             patch.spectral[pp1:pp1+b_iDim_patch-1, v] .=
                 SAtransform(patch.ibasis.data[2, v],
                             view(sharedSpectral, pp1:pp1+b_iDim_patch-1, v))
@@ -3535,6 +3569,9 @@ function splineTransform!(sharedSpectral::SharedArray{real},
             # Imaginary part
             pp1 = p * b_iDim_patch + 1
             tp1 = p * b_iDim_tile + 1
+            if has_wn_ahat
+                patch.ibasis.data[3, v].ahat .= _get_wavenumber_ahat(patch, v, p + 1)
+            end
             patch.spectral[pp1:pp1+b_iDim_patch-1, v] .=
                 SAtransform(patch.ibasis.data[3, v],
                             view(sharedSpectral, pp1:pp1+b_iDim_patch-1, v))
diff --git a/test/nesting_support.jl b/test/nesting_support.jl
new file mode 100644
index 0000000..98d12a3
--- /dev/null
+++ b/test/nesting_support.jl
@@ -0,0 +1,361 @@
+using Serialization
+
+# Tests for temporal-nesting support: collar interfaces, payload time
+# interpolation, and fine-mesh evaluation at arbitrary i-points.
+#
+# Two-way nesting scheme (DeMaria et al. 1992; Ooyama 2001): coarse→fine via
+# R3X border trio at the interface node (interior to the collar-extended
+# coarse patch), fine→coarse via collar-zone evaluation feeding the coarse
+# Galerkin loads.
+
+@testset "Nesting support (collar / lerp / point evaluation)" begin
+
+    # ── lerp_payload! ────────────────────────────────────────────────────────
+
+    @testset "lerp_payload! endpoints and midpoint" begin
+        meta_dummy = nothing  # payloads constructed directly
+        mk(vals) = InterfacePayload(:per_mode, :left, 1, 2,
+                                    reshape(collect(Float64, vals), 3, 2, 1))
+        p0 = mk(1:6)
+        p1 = mk(7:12)
+        dest = mk(zeros(6))
+
+        lerp_payload!(dest, p0, p1, 0.0)
+        @test dest.border == p0.border            # bitwise
+        lerp_payload!(dest, p0, p1, 1.0)
+        @test dest.border == p1.border            # bitwise
+        lerp_payload!(dest, p0, p1, 0.5)
+        @test dest.border ≈ 0.5 .* (p0.border .+ p1.border)
+        lerp_payload!(dest, p0, p1, 0.25)
+        @test dest.border ≈ 0.75 .* p0.border .+ 0.25 .* p1.border
+
+        @test_throws ArgumentError lerp_payload!(dest, p0, p1, 1.5)
+        bad = InterfacePayload(:per_mode, :right, 1, 2, zeros(3, 2, 1))
+        @test_throws ArgumentError lerp_payload!(bad, p0, p1, 0.5)
+    end
+
+    # ── evaluate_grid_ipoints, 1D R grid ─────────────────────────────────────
+
+    @testset "evaluate_grid_ipoints matches gridTransform! on mish points (1D)" begin
+        gp = SpringsteelGridParameters(
+            geometry="R", iMin=-5.0, iMax=5.0, num_cells=20,
+            BCL=Dict("u" => CubicBSpline.R0),
+            BCR=Dict("u" => CubicBSpline.R0),
+            vars=Dict("u" => 1))
+        g = createGrid(gp)
+        pts = getGridpoints(g)
+        for i in eachindex(pts)
+            g.physical[i, 1, 1] = sin(0.7 * pts[i]) + 0.3 * pts[i]^2
+        end
+        spectralTransform!(g)
+        gridTransform!(g)
+
+        xq = vec(pts)
+        out = evaluate_grid_ipoints(g, xq)
+        for s in 1:3
+            @test maximum(abs.(out[:, 1, s] .- g.physical[:, 1, s])) < 1e-12
+        end
+    end
+
+    @testset "evaluate_grid_ipoints linear exactness at off-mish points (1D)" begin
+        f(x)  = 3x + 5
+        gp = SpringsteelGridParameters(
+            geometry="R", iMin=0.0, iMax=10.0, num_cells=10,
+            BCL=Dict("u" => CubicBSpline.R0),
+            BCR=Dict("u" => CubicBSpline.R0),
+            vars=Dict("u" => 1))
+        g = createGrid(gp)
+        pts = getGridpoints(g)
+        for i in eachindex(pts)
+            g.physical[i, 1, 1] = f(pts[i])
+        end
+        spectralTransform!(g)
+
+        xq = [0.123, 3.456, 6.789, 9.876]
+        out = evaluate_grid_ipoints(g, xq)
+        @test out[:, 1, 1] ≈ f.(xq) atol=1e-10
+        @test out[:, 1, 2] ≈ fill(3.0, 4) atol=1e-10
+        @test out[:, 1, 3] ≈ zeros(4) atol=1e-9
+
+        @test_throws ArgumentError evaluate_grid_ipoints(g, [10.5])
+    end
+
+    # ── Collar interface: interior-node extraction, coarse→fine ─────────────
+
+    @testset "Collar interface transfers coarse trio at interior node" begin
+        f(x) = 3x + 5
+        # Coarse patch: nominal [0,10] + one-cell collar → [0,11], DX=1
+        gp_c = SpringsteelGridParameters(
+            geometry="R", iMin=0.0, iMax=11.0, num_cells=11,
+            BCL=Dict("u" => CubicBSpline.R0),
+            BCR=Dict("u" => CubicBSpline.R0),
+            vars=Dict("u" => 1))
+        # Fine patch: [10, 20], DX=0.5, R3X on the interface (left) side
+        gp_f = SpringsteelGridParameters(
+            geometry="R", iMin=10.0, iMax=20.0, num_cells=20,
+            BCL=Dict("u" => CubicBSpline.R3X),
+            BCR=Dict("u" => CubicBSpline.R0),
+            vars=Dict("u" => 1))
+        gc = createGrid(gp_c)
+        gf = createGrid(gp_f)
+
+        iface = PatchInterface(gc, gf, :right, :left, :i; is_stacked=true)
+        # Interface node x=10 is m=10 → array index 12 = b_iDim-2 (interior)
+        @test iface.primary_node_indices == (11, 12, 13)
+        @test iface.is_stacked
+        @test iface.coupling_matrix == Springsteel.COUPLING_MATRIX_2X
+
+        pts_c = getGridpoints(gc)
+        pts_f = getGridpoints(gf)
+        for i in eachindex(pts_c); gc.physical[i, 1, 1] = f(pts_c[i]); end
+        for i in eachindex(pts_f); gf.physical[i, 1, 1] = f(pts_f[i]); end
+        spectralTransform!(gc)
+        spectralTransform!(gf)
+
+        gridTransform!(gc)          # sets coarse spline .a
+        update_interface!(iface)    # coarse trio → fine ahat
+        gridTransform!(gf)
+
+        for i in eachindex(pts_f)
+            @test gf.physical[i, 1, 1] ≈ f(pts_f[i]) atol=1e-9
+        end
+    end
+
+    # ── Anti-freeze regression ───────────────────────────────────────────────
+    # The payload a patch donates must track its own interior dynamics.  A
+    # BC-based fine→coarse exchange (dual R3X) fails this: an R3X spline's
+    # border trio is identically its ahat, so the donated data would be the
+    # neighbor's own stale values and the interface would freeze.
+
+    @testset "Anti-freeze: donated payload tracks donor interior" begin
+        gp_c = SpringsteelGridParameters(
+            geometry="R", iMin=0.0, iMax=11.0, num_cells=11,
+            BCL=Dict("u" => CubicBSpline.R0),
+            BCR=Dict("u" => CubicBSpline.R0),
+            vars=Dict("u" => 1))
+        gp_f = SpringsteelGridParameters(
+            geometry="R", iMin=10.0, iMax=20.0, num_cells=20,
+            BCL=Dict("u" => CubicBSpline.R3X),
+            BCR=Dict("u" => CubicBSpline.R0),
+            vars=Dict("u" => 1))
+        gc = createGrid(gp_c)
+        gf = createGrid(gp_f)
+        iface = PatchInterface(gc, gf, :right, :left, :i; is_stacked=true)
+
+        pts_c = getGridpoints(gc)
+        for i in eachindex(pts_c); gc.physical[i, 1, 1] = 1.0; end
+        spectralTransform!(gc)
+        gridTransform!(gc)
+        p_before = compute_interface_payload(iface.metadata, gc)
+
+        # Perturb the coarse field near the interface (interior dynamics)
+        for i in eachindex(pts_c)
+            gc.physical[i, 1, 1] = 1.0 + 0.5 * exp(-(pts_c[i] - 9.5)^2)
+        end
+        spectralTransform!(gc)
+        gridTransform!(gc)
+        p_after = compute_interface_payload(iface.metadata, gc)
+
+        @test maximum(abs.(p_after.border .- p_before.border)) > 0.05
+    end
+
+    # ── Payload serialization round-trip (collar metadata) ──────────────────
+
+    @testset "Collar payload serializes byte-identically" begin
+        gp_c = SpringsteelGridParameters(
+            geometry="R", iMin=0.0, iMax=11.0, num_cells=11,
+            BCL=Dict("u" => CubicBSpline.R0),
+            BCR=Dict("u" => CubicBSpline.R0),
+            vars=Dict("u" => 1))
+        gp_f = SpringsteelGridParameters(
+            geometry="R", iMin=10.0, iMax=20.0, num_cells=20,
+            BCL=Dict("u" => CubicBSpline.R3X),
+            BCR=Dict("u" => CubicBSpline.R0),
+            vars=Dict("u" => 1))
+        gc = createGrid(gp_c)
+        gf = createGrid(gp_f)
+        iface = PatchInterface(gc, gf, :right, :left, :i; is_stacked=true)
+
+        pts_c = getGridpoints(gc)
+        for i in eachindex(pts_c); gc.physical[i, 1, 1] = sin(pts_c[i]); end
+        spectralTransform!(gc)
+        gridTransform!(gc)
+        p = compute_interface_payload(iface.metadata, gc)
+
+        buf = IOBuffer()
+        Serialization.serialize(buf, p)
+        seekstart(buf)
+        p2 = Serialization.deserialize(buf)
+        @test p2.border == p.border
+        @test p2.scheme === p.scheme && p2.side === p.side
+    end
+
+    # ── Out-of-place SAtransform honors R3X ahat ─────────────────────────────
+    # The tiled b→a path (3-arg splineTransform!) uses the allocating
+    # SAtransform(spline, b); it must reproduce SAtransform!'s ahat handling,
+    # otherwise nest patches lose their coupled borders on the tiled path.
+
+    @testset "SAtransform(spline, b) matches SAtransform! for R3X" begin
+        gp = SpringsteelGridParameters(
+            geometry="R", iMin=0.0, iMax=10.0, num_cells=10,
+            BCL=Dict("u" => CubicBSpline.R3X),
+            BCR=Dict("u" => CubicBSpline.R0),
+            vars=Dict("u" => 1))
+        g = createGrid(gp)
+        pts = getGridpoints(g)
+        for i in eachindex(pts)
+            g.physical[i, 1, 1] = sin(pts[i])
+        end
+        spectralTransform!(g)
+        spline = g.ibasis.data[1, 1]
+        Springsteel.CubicBSpline.set_ahat_r3x!(spline, 0.3, 0.7, 1.1, :left)
+
+        spline.b .= view(g.spectral, :, 1)
+        CubicBSpline.SAtransform!(spline)
+        a_inplace = copy(spline.a)
+        a_alloc = CubicBSpline.SAtransform(spline, view(g.spectral, :, 1))
+        @test a_alloc ≈ a_inplace atol=1e-14
+        # Rank-3 invariant: the constrained border trio IS the ahat
+        @test a_alloc[1:3] ≈ spline.ahat[1:3] atol=1e-14
+        @test any(!iszero, spline.ahat[1:3])
+    end
+
+    # ── evaluate_grid_ipoints, RiRk grid ─────────────────────────────────────
+
+    @testset "evaluate_grid_ipoints matches gridTransform! on mish points (RiRk)" begin
+        gp = SpringsteelGridParameters(
+            geometry="RiRk",
+            iMin=0.0, iMax=10.0, num_cells=10,
+            kMin=0.0, kMax=5.0, num_cells_k=5,
+            BCL=Dict("u" => CubicBSpline.R0, "w" => CubicBSpline.R0),
+            BCR=Dict("u" => CubicBSpline.R0, "w" => CubicBSpline.R0),
+            BCB=Dict("u" => CubicBSpline.R0, "w" => CubicBSpline.R0),
+            BCT=Dict("u" => CubicBSpline.R0, "w" => CubicBSpline.R0),
+            vars=Dict("u" => 1, "w" => 2))
+        g = createGrid(gp)
+        pts = getGridpoints(g)          # (iDim*kDim, 2), i-outer k-inner
+        for i in axes(pts, 1)
+            x, z = pts[i, 1], pts[i, 2]
+            g.physical[i, 1, 1] = sin(0.5x) * cos(0.8z)
+            g.physical[i, 2, 1] = 0.1 * x * z + z^2
+        end
+        spectralTransform!(g)
+        gridTransform!(g)
+
+        kDim = g.params.kDim
+        iDim = g.params.iDim
+        xq = [pts[(q - 1) * kDim + 1, 1] for q in 1:iDim]   # unique i-mish points
+        out = evaluate_grid_ipoints(g, xq)
+        for v in 1:2, s in 1:5
+            @test maximum(abs.(out[:, v, s] .- g.physical[:, v, s])) < 1e-10
+        end
+    end
+
+end
+
+# ── RL (cylindrical) nesting support ─────────────────────────────────────────
+
+@testset "RL nesting support" begin
+
+    @testset "evaluate_grid_points matches gridTransform! on mish points (RL)" begin
+        gp = SpringsteelGridParameters(
+            geometry="RL", iMin=0.0, iMax=50.0, num_cells=10,
+            BCL=Dict("u" => NaturalBC()), BCR=Dict("u" => NaturalBC()),
+            vars=Dict("u" => 1))
+        g = createGrid(gp)
+        pts = getGridpoints(g)
+        for i in 1:size(pts, 1)
+            r, λ = pts[i, 1], pts[i, 2]
+            g.physical[i, 1, 1] = 0.1 * r + r * cos(λ) + 0.5 * r * sin(2λ)
+        end
+        spectralTransform!(g)
+        gridTransform!(g)
+
+        # gridTransform reconstructs ring ri with wavenumbers k ≤ ri only;
+        # pass the same per-point truncation for an exact comparison.
+        kmax = Int[]
+        for r in 1:g.params.iDim
+            ri = r + g.params.patchOffsetL
+            append!(kmax, fill(ri, 4 + 4 * ri))
+        end
+        out = evaluate_grid_points(g, pts; kmax = kmax)
+        for s in 1:5
+            err = maximum(abs.(out[:, 1, s] .- g.physical[:, 1, s]))
+            @test err < 1e-8
+        end
+    end
+
+    @testset "tiled RL splineTransform! honors the per-wavenumber registry" begin
+        # Disc-in-annulus fixture with distinct k=0 / k=1-real / k=1-imag borders
+        gp_annulus = SpringsteelGridParameters(
+            geometry="RL", iMin=20.0, iMax=100.0, num_cells=10,
+            BCL=Dict("u" => NaturalBC()), BCR=Dict("u" => NaturalBC()),
+            vars=Dict("u" => 1))
+        gp_disc = SpringsteelGridParameters(
+            geometry="RL", iMin=0.0, iMax=20.0, num_cells=5,
+            BCL=Dict("u" => NaturalBC()), BCR=Dict("u" => FixedBC()),
+            vars=Dict("u" => 1))
+        g_annulus = createGrid(gp_annulus)
+        g_disc = createGrid(gp_disc)
+        f(r, λ) = (2r + 5) + 0.3r * cos(λ) + 0.7r * sin(λ)
+        pts_a = getGridpoints(g_annulus)
+        pts_d = getGridpoints(g_disc)
+        for i in 1:size(pts_a, 1); g_annulus.physical[i, 1, 1] = f(pts_a[i, 1], pts_a[i, 2]); end
+        for i in 1:size(pts_d, 1); g_disc.physical[i, 1, 1] = f(pts_d[i, 1], pts_d[i, 2]); end
+        spectralTransform!(g_annulus)
+        spectralTransform!(g_disc)
+
+        iface = PatchInterface(g_annulus, g_disc, :left, :right, :i)
+        gridTransform!(g_annulus)
+        update_interface!(iface)
+
+        # Known-correct registry-aware path
+        gridTransform!(g_disc)
+        want = copy(g_disc.physical)
+
+        # Tiled path: 3-arg splineTransform! (patch splines) + tileTransform!
+        shared = SharedArray{Float64,2}(size(g_disc.spectral))
+        shared .= g_disc.spectral
+        tile = createGrid(gp_disc)
+        splineTransform!(shared, g_disc, tile)
+        tileTransform!(shared, tile, tile.physical, tile.spectral)
+        for s in 1:5
+            err = maximum(abs.(tile.physical[:, 1, s] .- want[:, 1, s]))
+            @test err < 1e-10
+        end
+
+        # Single-tile 2-arg path on the coupled grid object itself
+        splineTransform!(shared, g_disc)
+        tileTransform!(shared, g_disc, g_disc.physical, g_disc.spectral)
+        for s in 1:5
+            err = maximum(abs.(g_disc.physical[:, 1, s] .- want[:, 1, s]))
+            @test err < 1e-10
+        end
+    end
+
+    @testset "offset RL annulus follows global ring numbering" begin
+        # Nest-annulus convention: patchOffsetL carries the GLOBAL ring offset
+        # explicitly, while spectralIndexL stays 1 (the annulus is its own
+        # patch, so tile spectral windows are annulus-relative).
+        mubar = 3
+        gp = SpringsteelGridParameters(
+            geometry="RL", iMin=150.0e3, iMax=300.0e3, num_cells=50,
+            patchOffsetL = 50 * mubar,    # 50 inner cells of the same DX
+            BCL=Dict("u" => NaturalBC()), BCR=Dict("u" => NaturalBC()),
+            vars=Dict("u" => 1))
+        g = createGrid(gp)
+        @test g.params.patchOffsetL == 50 * mubar
+        @test g.params.spectralIndexL == 1
+        pts = getGridpoints(g)
+        # First ring's lpoints follows the global index: 4 + 4*(1 + offset)
+        lp1 = 4 + 4 * (1 + g.params.patchOffsetL)
+        @test count(x -> x ≈ pts[1, 1], pts[1:lp1 + 8, 1]) == lp1
+
+        # A single tile of the annulus composes the offset (rings match the
+        # patch) while keeping its spectral window patch-relative.
+        tiles = calcTileSizes(g, 1)
+        @test tiles[1].params.spectralIndexL == 1
+        @test tiles[1].params.patchOffsetL == g.params.patchOffsetL
+        @test size(getGridpoints(tiles[1]), 1) == size(pts, 1)
+    end
+end
diff --git a/test/runtests.jl b/test/runtests.jl
index fb6e94d..a947d69 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -26,6 +26,7 @@ const TEST_GROUP = get(ENV, "TEST_GROUP", "all")
     TEST_GROUP in ("all", "r3x")            && include("r3x.jl")
     TEST_GROUP in ("all", "bc")             && include("boundary_conditions.jl")
     TEST_GROUP in ("all", "multipatch")     && include("multipatch.jl")
+    TEST_GROUP in ("all", "nesting_support") && include("nesting_support.jl")
     TEST_GROUP in ("all", "tile_multipatch") && include("tile_multipatch.jl")
     TEST_GROUP in ("all", "operator_algebra") && include("operator_algebra.jl")
     TEST_GROUP in ("all", "solver_problem")   && include("solver_problem.jl")