Velocity Solver

coriolis_parameter(lat_deg)

f = 2Ω sin(φ) with the equatorial floor |f| ≥ F_MIN applied.

solve_baroclinic!(u, v, p, tau_x, tau_y, grid)

Compute baroclinic horizontal velocities from the frictional-geostrophic balance.

Arguments

• u : Field{Face, Center, Center} — zonal velocity output
• v : Field{Center, Face, Center} — meridional velocity output
• p : Field{Center, Center, Center} — hydrostatic pressure, bar (halos filled)
• tau_x: Field{Face, Center, Nothing} or Matrix — zonal wind stress, N m⁻² (at U-face)
• tau_y: Field{Center, Face, Nothing} or Matrix — meridional wind stress (at V-face)
• grid : the shared ImmersedBoundaryGrid

BarotropicSolver

Pre-assembled sparse system for the barotropic streamfunction. Assembled once per land-sea mask and reused every timestep.

Fields:

• A : sparse coefficient matrix
• rhs : right-hand side vector (updated each timestep)
• psi_vec : solution vector (nlon*nlat unknowns, land cells zeroed)
• cell_index: (nlon, nlat) → linear index into the ocean-cell DOF ordering
• islands : detected island list

build_barotropic_solver(ocean_mask, f, H, dx, dy, lat_deg) -> BarotropicSolver

Assemble the sparse linear system for the barotropic streamfunction.

Arguments

• ocean_mask: Bool (nlon, nlat)
• f : Coriolis (nlat,), s⁻¹, with F_MIN floor applied
• H : ocean depth (nlon, nlat), m, positive downward
• dx : zonal grid spacing (nlat,), m
• dy : meridional grid spacing, m (uniform assumed)
• r_bt : barotropic friction field (nlon, nlat), s⁻¹

build_barotropic_solver(grid::AbstractGrid) -> BarotropicSolver

Convenience constructor: extracts ocean mask, bathymetric depth, Coriolis, and grid spacings from the Oceananigans ImmersedBoundaryGrid, then calls the array-based assembler.

A column (i,j) is ocean when bottom_height[i,j] < 0; H = −bottom_height.

compute_jebar_forcing(T, S, grid) -> Matrix{Float64}

Compute the JEBAR (Joint Effect of Baroclinicity And Relief) forcing for the barotropic-vorticity RHS, in units of kg m⁻³ s⁻² (= Pa m⁻², equivalent to the wind-stress curl term before division by ρ₀).

The CLIMBER-X formulation (jbar.f90):

bp(i, j, k) = − ∫{zbottom}^{z_k} ρ g dz (Pa, zero at the bottom)

For each face F between two T-cells, integrate bp · dz from a common face-bottom level k_face = max(k_bot_left, k_bot_right) — using the deeper of the two cells' active ranges so any sub-shelf contributions that would differ between the two cells are excluded. Then form

EF = Σ{k≥kface} ½(bpleft(k) + bpright(k)) dz(k) (Pa·m) invHF = 1 / Σ{k≥kface} dz(k) (m⁻¹)

The T-cell JEBAR is the central difference on those face values:

ubar_jbar[i,j] = ∂(1/H)/∂x · ∂E/∂y − ∂(1/H)/∂y · ∂E/∂x

with ∂(1/H)/∂x = (invHE − invHW)/Δx (and similarly for y), and likewise for ∂E. At ocean–land faces the face-bottom level lies above any active layer, so both EF and invHF vanish identically — the closed-boundary treatment falls out of the staggering rather than a special case.

Periodic in λ; bounded north/south.

compute_rbt(ocean_mask, H, dx, dy) -> Matrix{Float64}

Compute spatially variable barotropic friction: RBTBASE × RBTFAC near continental boundaries and in shallow water (H < 500 m), RBTBASE elsewhere.

solve_barotropic!(solver, tau_x, tau_y, H, f, dx, dy, ocean_mask;
                   cos_phi_T=nothing, ubar_jbar=nothing) -> Matrix{Float64}

Solve for the barotropic streamfunction ψ in m³ s⁻¹ (1 Sv = 10⁶ m³ s⁻¹). Returns ψ on the full (nlon, nlat) grid (land cells = 0).

Wind stress tau_x, tau_y in N m⁻²; the assembler divides by RHO_0 so the equation is dimensionally consistent.

Spherical curl: when cos_phi_T is supplied (length nlat, cosines of T-row latitudes), the meridional piece of curl(τ/H) is computed in the CLIMBER-X form (1/(R cosφ_T))·∂(τx·cosφ/H)/∂φ rather than the flat-earth ∂(τx/H)/∂y. At high latitudes the difference (τx·tanφ/(R·H) metric correction) is order-unity. Without cos_phi_T the flat-earth form is used (kept for tests that build a test grid without lat metadata).

diagnose_w!(w, u, v, grid)

Integrate continuity upward from the seafloor to diagnose w.

The surface value w[:, :, Nz+1] is the rigid-lid residual; it is ≈ 0 only if the depth-integrated horizontal divergence of (u, v) is ≈ 0.

IslandInfo

Holds connectivity data for one island.

Fields:

• id : island index (1-based)
• perimeter : vector of (i,j) tuples marking ocean cells adjacent to this island
• interior_mask: Bool matrix marking the island's land interior cells

detect_islands(ocean_mask) -> Vector{IslandInfo}

Detect islands (connected land regions not touching the domain boundary) from the binary ocean mask. Returns one IslandInfo per island.

The main continent + domain boundary are excluded — only interior land cells enclosed by ocean on all sides count as islands.

ocean_mask: Bool (nlon, nlat), true = ocean cell.