Split StokesDrift: function-mode StokesDrift + Field-mode FieldStokesDrift

[glwagner]

Summary

Splits the Stokes-drift abstraction into two concrete subtypes of a new AbstractStokesDrift supertype:

• StokesDrift (function mode, fully backward-compatible): all derivative slots are callables f(x, y, z, t[, params]). Restored to its original 10-parameter form; no uˢ, vˢ, wˢ slots, no Field dispatch, no validation. Existing user code is unchanged.
• FieldStokesDrift (Field mode): six prognostic Field slots (uˢ, vˢ, wˢ, ∂t_uˢ, ∂t_vˢ, ∂t_wˢ) at the C-grid velocity locations for couplings where the Stokes drift comes from an external wave model.

```julia
sd = FieldStokesDrift(grid) # all six Fields default-allocated
sd = FieldStokesDrift(grid; uˢ=my_uˢ_field) # override individual slots
NonhydrostaticModel(grid; stokes_drift=sd, ...)
```

The 7-parameter FieldStokesDrift struct stores only the Field slots — no dead zerofunction placeholders. The per-slot kwarg constructor defaults each missing argument to a freshly-allocated Field at the matching staggered location, so callers can share Fields with other code by overriding individual kwargs.

Behavior

• Spatial derivatives in x_curl_Uˢ_cross_U, y_curl_Uˢ_cross_U, z_curl_Uˢ_cross_U are computed inline from the stored uˢ, vˢ, wˢ Fields via the staggered FD operators (∂xᶠᶠᶜ, ∂yᶠᶠᶜ, ∂zᶠᶜᶠ, ∂zᶜᶠᶠ, ∂xᶠᶜᶠ, ∂yᶜᶠᶠ + ℑ interpolators). Both ∂z uˢ, ∂z vˢ and the cross-derivative terms ∂y uˢ, ∂x vˢ participate — no curl approximation drops them.
• ∂t_uˢ, ∂t_vˢ, ∂t_wˢ are read directly via getindex at the velocity location.
• compute_stokes_drift!(::FieldStokesDrift, grid) fills wˢ and ∂t_wˢ by vertical integration of `∂_z wˢ = -∂_x uˢ - ∂_y vˢ` (same continuity pattern as the hydrostatic model's `compute_w_from_continuity!`), starting from zero at the bottom. No-op on StokesDrift and UniformStokesDrift.
• NonhydrostaticModel.update_state! invokes `compute_stokes_drift!` between the velocity halo fill and auxiliary-field computation, so wˢ and ∂t_wˢ are always in sync with uˢ, vˢ, ∂t_uˢ, ∂t_vˢ.

Shared structure

The three *_curl_Uˢ_cross_U methods share a single body (the (∇×uˢ) × uᴱ algebra) and dispatch on Union{StokesDrift, FieldStokesDrift}; per-derivative _∂{x,y,z}_{u,v,w}ˢ_* helpers branch on the concrete type. UniformStokesDrift retains its own specialized curl-cross methods because it has a simpler structure.

Backward compatibility

Function-mode user code is completely unchanged. StokesDrift(; ∂z_uˢ, ...) returns the same struct as before, with the same type-parameter layout.

Motivation

Developed against the wind-drift instability example in NumericalEarth/Ripple.jl, where the wave-model pseudomomentum p(z) = Q(z)·A·K equals the deep-water Stokes drift when the wave action is calibrated to the target steepness, and the analytic pseudomomentum tendency Q(z)·∂t(A·K) equals the Stokes acceleration the ocean must see for the wave-mean energy budget to close.

Test plan

`test/test_stokes_drift.jl` gains a `FieldStokesDrift` testset:

• All six Field slots allocated by `FieldStokesDrift(grid)` at the expected staggered locations.
• Per-slot kwarg override (e.g., `FieldStokesDrift(grid; uˢ=my_uˢ)`).
• Zero state gives zero contributions.
• Linear vertical shear `uˢ = α·z`, uniform `w = W` → `x_curl_Uˢ_cross_U = W·α` at every interior point.
• Cross-derivative term `uˢ = β·y, v = V` → `V·β` (the term `UniformStokesDrift` drops).
• `∂t_uˢ, ∂t_vˢ` read directly with the right sign.
• `compute_stokes_drift!` fills `wˢ` from continuity when `uˢ, vˢ` have nonzero divergence.
• `compute_stokes_drift!` is a no-op on function-mode `StokesDrift`.

All 326 Stokes-drift tests pass.

🤖 Generated with Claude Code