`PerturbationAdvection`: support 3D / laterally-varying density at open boundaries

## Current behavior The `density` kwarg of `PerturbationAdvection` ([`src/BoundaryConditions/perturbation_advection.jl#L121-L124`](https://github.com/CliMA/Oceananigans.jl/blob/main/src/BoundaryConditions/perturbation_advection.jl#L121-L124)) is used to convert density-weighted prognostic variables to/from intensive form via: ```julia @inline to_intensive(::Nothing, ψ, k) = ψ @inline to_intensive(ρ, ψ, k) = ψ / @inbounds ρ[1, 1, k] @inline to_extensive(::Nothing, ψ, k) = ψ @inline to_extensive(ρ, ψ, k) = @inbounds ρ[1, 1, k] * ψ ``` The hardcoded `[1, 1, k]` indexing assumes density is **column-uniform** (no horizontal variation). This is appropriate for anelastic-style reference profiles `ρᵣ(z)` stored as `Field{Nothing, Nothing, Center}(grid)` — the documented use case in the docstring. ## Use case requiring laterally-varying density For **regional atmospheric hindcasts** driven by reanalysis (ERA5, GFS, MERRA-2, …), the prescribed boundary state varies along the lateral boundary in `(x, y)`. Density at the boundary face is diagnosed from the driver's `(T, q_v, p)` via the equation of state: ``` ρ_boundary(i, j, k, t) = p_FTS(i, j, k, t) / (R_d · T_v_FTS(i, j, k, t)) ``` For `PerturbationAdvection` to correctly convert `ρu → u → ρu_new` at the boundary face, it should index `ρ_boundary` at the actual face position `[i, j, k]` rather than a single column `[1, 1, k]`. The current code silently uses the wrong density at every boundary point except the centerline. Similarly, fully-compressible models with a 3D prognostic density would want to pass `model.dynamics.density` directly and have PA index it correctly per-face. ## Proposed fix Thread `(i, j)` through `to_intensive` / `to_extensive`, with location-aware dispatch to keep the existing anelastic reference-profile case backwards-compatible: ```julia # Existing — Nothing-located horizontal axes (anelastic reference profile) @inline to_intensive(ρ::AbstractField{Nothing, Nothing, Center}, ψ, i, j, k) = ψ / @inbounds ρ[1, 1, k] # New — full 3D density (compressible / regional hindcast) @inline to_intensive(ρ::AbstractField{Center, Center, Center}, ψ, i, j, k) = ψ / @inbounds ρ[i, j, k] # Same for to_extensive ``` Plus pass `(i, j)` through `step_right_open_boundary!` / `step_left_open_boundary!` (and the corresponding south/north/bottom/top variants). ## Scope Small patch — roughly 10–15 lines. No API change beyond the addition of the `[i, j, k]` indexing path; users passing a `Field{Nothing, Nothing, Center}` density (the current documented anelastic use case) continue to work unchanged. ## Context Came up while planning downstream Breeze (`NumericalEarth/Breeze.jl`) open-boundary support for compressible regional hindcasts using `OpenBoundaryCondition(...; scheme=PerturbationAdvection(...))` on density-weighted prognostic variables (ρu, ρv, ρθ, ρqᵗ). With the current indexing, the boundary mass flux at off-centerline faces is biased by the density mismatch. cc @glwagner