Ocean vertical-coordinate + ALE remap state.
!! Ocean vertical-coordinate + ALE remap state. module rdb_ocean_vcoord !! Holds the vertical-coordinate configuration and the per-step !! target grid that the ALE remap step relamps `multilayer.h_layer` !! + every `multilayer.tracers(t)%hTr` onto. The coastal path has !! `VCOORD_SIGMA`, `VCOORD_ZSIGMA`, `VCOORD_ZSTAR`, `VCOORD_ZSTAR_SIGMA`, !! `VCOORD_ZSTAR_FULL` plus `rdb_remap_column`; this module is the !! C-grid counterpart for the ocean dynamical core. !! !! Phase 5g build-out status (this commit, Layer 2): !! !! - `target_h(nx, ny, nz_ml)` is allocated on init + bound to the !! device via `enter_data` / `exit_data`. !! - `compute_target_h(this, total_h, eta)` populates `target_h` !! per the current `coord_type`. Working bodies: !! VCOORD_EULERIAN_Z — H · dsig(k) (η ignored) !! VCOORD_SIGMA — (H + η) · dsig(k) !! VCOORD_ZSTAR — MOM6 z*: the z_fixed nominal profile !! dilated by (H + η)/H over bed fillers !! VCOORD_ZSIGMA — smoothstep blend(sigma, fixed z-levels) !! VCOORD_ZSTAR_SIGMA — smoothstep blend(sigma, z*-lite) !! VCOORD_ZSTAR_FULL — per-column z_ref + vanishing-layer floors !! - `build_zref_full(this, h_bed)` populates `z_ref(:, :, 0:nz_ml)` !! per column from local bathymetry. Call once at init (or any !! time the bathymetry changes); the per-step `compute_target_h` !! then walks the cached `z_ref` table. !! - Driver call site (remap kernel invocation between slow stages !! in `ocean_dyn_step_split`) is NOT wired yet — !! simulation state is untouched. External callers can invoke !! `compute_target_h` / `build_zref_full` to inspect the target !! grid without advancing dynamics. !! !! Collaborator hand-off — remaining Phase 5g work (Layer 3): !! 1. Driver wiring in `ocean_dyn_step_split`: call !! `compute_target_h(total_h_2d, dyn%bt_work%bt_eta)` then invoke the !! remap kernel on `multilayer.h_layer` + each `tracers(t)%hTr` !! + face velocities. Recompute `bt_eta` from the new sum of !! `h_layer - bt_H_ref`. !! 2. Face-velocity remap adapter — MOM6 reconstructs u at centres, !! remaps, then projects back to faces with a divergence-free !! correction. Donor-cell on `(h·u)_face` is the simpler !! fallback. !! 3. `VANISHING_LAYER_TOL`-gated CWC for `VCOORD_ZSTAR_FULL`, !! lifted from the coastal multilayer kernels. !! !! See `docs/ROADMAP_OCEAN.md` Phase 5g for the full scope. #ifdef LFORTRAN_PASSING use rdb_constants, only: wp, REMAP_PPM, H_VANISHED, & VCOORD_LAGRANGIAN, VCOORD_EULERIAN_Z, & VCOORD_SIGMA, VCOORD_ZSIGMA, VCOORD_ZSTAR, & VCOORD_ZSTAR_FULL, VCOORD_ZSTAR_SIGMA, & VCOORD_Z_FIXED, VCOORD_RHO, VCOORD_HYCOM #else use rdb_constants, only: NZ_STACK_MAX, wp, REMAP_PPM, H_VANISHED, & VCOORD_LAGRANGIAN, VCOORD_EULERIAN_Z, & VCOORD_SIGMA, VCOORD_ZSIGMA, VCOORD_ZSTAR, & VCOORD_ZSTAR_FULL, VCOORD_ZSTAR_SIGMA, & VCOORD_Z_FIXED, VCOORD_RHO, VCOORD_HYCOM #endif use rdb_grid, only: hgrid_t use rdb_vcoord, only: STRETCH_UNIFORM, STRETCH_LOG, parse_vcoord_type use rdb_eos, only: eos_t, eos_density_point use, intrinsic :: iso_fortran_env, only: int64 use rdb_mem_report, only: arr_bytes implicit none private #ifdef LFORTRAN_PASSING integer, parameter :: NZ_STACK_MAX = 64 !! LFortran 0.64 workaround: module-local copy of the rdb_constants value !! (an imported parameter used as an explicit-shape dummy bound inside a !! PURE call becomes an impure getter under LFortran). Keep in sync (=64). #endif public :: ocean_vcoord_t public :: parse_ocean_vcoord_type public :: invert_density_targets public :: ocean_vcoord_z_fixed_target public :: ocean_vcoord_z_fixed_target_uniform public :: ocean_vcoord_zstar_target public :: ocean_vcoord_set_z_fixed_profile public :: ocean_vcoord_closed_face_masks public :: ocean_vcoord_eta0_target public :: ocean_vcoord_k_top_from_target public :: ocean_vcoord_k_bot_from_target public :: ocean_vcoord_count_ledges public :: ocean_vcoord_count_bed_steps public :: VCOORD_EULERIAN_Z public :: VCOORD_LAGRANGIAN public :: VCOORD_Z_FIXED public :: VCOORD_RHO public :: VCOORD_HYCOM public :: STRETCH_UNIFORM, STRETCH_LOG ! ---- RHO (isopycnal) inversion parameters ---- integer, parameter :: NR_ITERS = 8 !! Fixed (GPU-uniform) Newton iteration budget for the !! density→depth inversion. Unrolled, no data-dependent while. real(wp), parameter :: NR_TOL = 1.0e-12_wp !! Newton convergence tolerance — tested on |delta| AFTER xi += delta. real(wp), parameter :: NR_OFFSET = 1.0e-6_wp !! Out-of-range nudge applied only when the boundary gradient ≈ 0. ! ---- VCOORD_HYCOM z* nominal-floor source ---- integer, parameter :: HYCOM_FLOOR_SIGMA = 0 !! Floor at `Σ dsig·(H+η)` — a column FRACTION, i.e. a sigma floor. !! Only the fallback for a slot no setup path configured !! (`z_fixed_h_ref <= 0`); it was the only floor before 2026-10-02, !! when it set 95 % of the 1-degree Southern Ocean's interfaces and !! made `hycom` terrain-following there (audit finding H1). integer, parameter :: HYCOM_FLOOR_UNIFORM = 1 !! Floor at `Σ (z_fixed_h_ref/nz)·(H+η)/H` — uniform z* in metres !! (`&vcoord_nml z_fixed_profile = "uniform"`, `max_depth/nz`). integer, parameter :: HYCOM_FLOOR_PROFILE = 2 !! Floor at `Σ z_fixed_dz(k)·(H+η)/H` — a stretched z* resolution in !! metres (`z_fixed_profile = "list" | "tanh"`). ! ---- VCOORD_Z_FIXED rigid-top partial cell ---- real(wp), parameter :: Z_FIXED_TOP_PARTIAL_FRAC = 0.1_wp !! Minimum PARTIAL TOP CELL thickness, as a fraction of the nominal !! layer spacing `h_nominal`, for `VCOORD_Z_FIXED` under a rigid top !! (`z_top > 0`, i.e. an ice-shelf cavity). The ice base cuts the !! first live layer wherever the draft crosses a nominal level; left !! unguarded that cut can leave a sliver of arbitrarily small !! thickness (`0 < h << h_nominal`) sitting against the ice. A layer !! that thin is live (it clears `H_VANISHED`) but violates every !! thin-layer CFL the column has, so when the cut would leave less !! than `Z_FIXED_TOP_PARTIAL_FRAC*h_nominal` the sliver is merged !! into the layer BELOW and the vacated index becomes an inert !! filler — exactly mirroring how the BED side treats its own !! slivers today (a would-be sub-`h_min` bed cell collapses to !! `h_min` and hands its water to the layer above). !! !! The two thresholds differ on purpose and the asymmetry is the !! honest one: the bed's threshold is `Z_FIXED_BED_PARTIAL_MIN` !! (below), the top's is a fraction of the spacing. `0.1` is !! MITgcm's `hFacMin` default, which is the minimum partial-cell !! fraction Losch (2008, JGR 113 C08043, §2.1) used for exactly !! this ice-shelf partial-top-cell problem. Not a namelist knob: !! it is inert unless `z_top > 0`, and a cavity on a z-like !! coordinate is itself a new, fenced configuration. real(wp), parameter :: Z_FIXED_BED_PARTIAL_MIN = H_VANISHED !! Minimum LIVE partial BOTTOM cell thickness for `VCOORD_Z_FIXED` !! — the bed mirror of `Z_FIXED_TOP_PARTIAL_FRAC`, expressed as an !! absolute thickness rather than a fraction of the spacing. !! !! **What it closes.** The bed branch used to floor its partial !! cell at `h_min` (`zstar_h_min`) alone, and the contract for this !! family is `zstar_h_min <= H_VANISHED`. That leaves the half-open !! band `(zstar_h_min, H_VANISHED] = (1e-4, 1.5e-4]` with the !! default `h_min`, in which a layer is **live to the coordinate !! and vanished to every consumer** — the EOS substitutes reference !! T/S, the remap reads its concentration as zero, `k_top`-style !! scans walk past it. A millimetre of `η` is enough to flip a bed !! remainder across it, and on !! `validation_examples/ocean/isomip_plus/ocean0_idealised_zfixed.nml` !! five columns did so 165 times in five days. With this floor the !! bed partial cell is either a live layer strictly ABOVE the marker !! or an inert filler at exactly `h_min`; nothing lands in between. !! !! **Why not `Z_FIXED_TOP_PARTIAL_FRAC*h_nominal` at the bed too.** !! That is the full mirror, and it is a different, larger change: !! with a 20 m spacing it would merge away every bed cell thinner !! than 2 m, moving the level structure — and the effective !! bathymetry — of every shipped `z_fixed` configuration. The !! thin-layer CFL argument that motivates the top's 10 % applies at !! the bed as well and that mirror may still be worth taking, but it !! is a bathymetry change and belongs in its own slice with its own !! validation. This constant fixes the defect that is actually !! diagnosed: a 0.05 mm window, so in practice only a column whose !! bed remainder lands inside it moves at all. type :: ocean_vcoord_t logical :: is_init = .false. !! True between `init` and `destroy`. Prefer this to !! `allocated(...)` — tracks GPU device attachment too. ! ---- Active coordinate ---- integer :: coord_type = VCOORD_EULERIAN_Z !! Selected vertical-coordinate variant. ! ---- Per-layer fractional thickness ---- ! Sums to 1.0 across the column. For `VCOORD_SIGMA` this is the ! target σ stencil (and `VCOORD_Z_FIXED` / `VCOORD_ZSTAR` fall back ! to it when no nominal profile was resolved, `z_fixed_h_ref = 0`): ! target_h(i,j,k) = (H(i,j) + eta(i,j)) * dsig(k) ! For `VCOORD_ZSTAR_FULL` dsig is a fallback used when the per- ! column `z_ref` table is not populated. real(wp), allocatable :: dsig(:) !! Per-layer σ-fraction. Sums to 1.0; size `nz_ml`. ! ---- Global z-level reference profile ---- ! Reference z-interfaces in metres (positive-down), `z_ref_global(0) ! = 0` is the surface, `z_ref_global(nz_ml)` is the deepest ! reference interface. Drives the z-level branch of `VCOORD_ZSIGMA` ! and the z*-lite branch of `VCOORD_ZSTAR_SIGMA`. Default at init ! is uniform 0..1 (normalised) — the namelist parser populates it ! with absolute depths when those cases are activated. real(wp), allocatable :: z_ref_global(:) !! Global reference z-interfaces (m, positive-down), shape !! `0:nz_ml`. Used by ZSIGMA / ZSTAR_SIGMA (NOT by ZSTAR, whose !! nominal profile is the `z_fixed` one, `z_fixed_zi`). ! ---- Per-column target thickness ---- ! Recomputed every outer step from (H, eta) per the coord_type ! case. Consumed by the remap kernel (Layer 3) to advance ! `multilayer.h_layer` and every `tracers(t)%hTr`. real(wp), allocatable :: target_h(:, :, :) !! Target layer thickness (m), shape `(nx, ny, nz_ml)`. ! ---- Per-column z* reference profile ---- ! Anchored to local bathymetry for `VCOORD_ZSTAR_FULL`. Indexed ! from `k=0` (surface) to `k=nz_ml` (bed) — opposite of the ! bottom-up state convention so the surface anchor is at index 0 ! (matches the coastal convention in `vcoord_target_dz_column_zstar_full`). ! Populated by `build_zref_full(h_bed)`; consumed by the ! ZSTAR_FULL branch of `compute_target_h`. real(wp), allocatable :: z_ref(:, :, :) !! Per-column z* reference (m), shape `(nx, ny, 0:nz_ml)`. ! ---- Geopotential depth of the column TOP ---- ! The rigid-lid seam of the z-like families. `z_top(i,j)` is the ! depth (m, positive down, `>= 0`) of the top of the WATER column ! below the `z = 0` datum, i.e. `metrics%z_draft` under an ! ice-shelf cavity and identically `0` everywhere else (open ! ocean, and every non-cavity run). Filled once at configure ! (`configure_ocean_cavity`) — the draft is static — and never ! touched again, so the remap driver's signature is unchanged and ! no second static 2-D array is threaded through the `pure` call ! chain. ! ! Allocated UNCONDITIONALLY with `source = 0.0_wp`, so it is always ! safe to hand to an explicit-shape device dummy: unlike ! `metrics%z_draft` (a `(1,1)` placeholder when no cavity is ! configured) this is always `(nx_total, ny_total)`. ! ! Consumed by the `VCOORD_Z_FIXED` branch of `compute_target_h`, ! which measures its nominal interface depths from `z = 0` and ! clips the stack against `z_top`. `z_top ≡ 0` ⇒ every geometric ! branch reproduces its pre-cavity arithmetic bit-for-bit. real(wp), allocatable :: z_top(:, :) !! Geopotential depth of the column top (m, positive down), !! shape `(nx, ny)`. `0` = the free surface at `z = 0`. ! ---- Isopycnal (VCOORD_RHO) target densities ---- ! Monotone-increasing nominal interface potential densities ! (kg/m³) referenced to `rho_ref_pressure`. Indexed `0:nz_ml`: ! `rho_target(0)` is the lightest (surface, maps to the k=nz ! interface in bottom-up state); `rho_target(nz_ml)` the densest ! (bed, k=1 interface). The `compute_target_h_rho` inversion ! places interior interfaces where the reconstructed column ! density equals each interior target. Sized + populated only ! when `coord_type == VCOORD_RHO`; ignored otherwise. real(wp), allocatable :: rho_target(:) !! Target interface potential densities (kg/m³), shape `0:nz_ml`. ! ---- ALE remap workspaces ---- ! Allocated once at init + enter_data'd to device. Previously ! these were allocated per call inside `ocean_apply_ale_remap_*` ! which generated thousands of host-allocated buffers per day ! that the device-side DCs in the remap step had to implicitly ! transfer back and forth. Persistent device-mapped scratch ! eliminates that overhead entirely. real(wp), allocatable :: remap_total_h(:, :) !! Column-total h_layer scratch, shape `(nx, ny)`. real(wp), allocatable :: remap_h_ref(:, :) !! H reference (total_h − bt_eta) scratch, shape `(nx, ny)`. real(wp), allocatable :: remap_h_old(:, :, :) !! Snapshot of `h_layer` before the remap, shape `(nx, ny, nz)`. real(wp), allocatable :: remap_conc_t(:, :, :) !! Layer-mean T concentration scratch for the `VCOORD_RHO` !! density inversion, shape `(nx, ny, nz)`. Built from !! `hTr / remap_h_old` (vanishing-layer-guarded) once per remap. real(wp), allocatable :: remap_conc_s(:, :, :) !! Layer-mean S concentration scratch for `VCOORD_RHO`, shape !! `(nx, ny, nz)`. ! ---- Tuning knobs ---- integer :: remap_method = REMAP_PPM !! ALE remap reconstruction order (REMAP_PCM/PLM/PPM/PPM_H4/PQM). !! PQM falls back to PPM for nz < 5 (see `remap_column_pqm`). real(wp) :: zstar_h_surf_target = 5.0_wp !! Surface-layer thickness anchor for `VCOORD_ZSTAR_FULL` (m). real(wp) :: zstar_h_min = 1.0e-4_wp !! Bed-side vanishing-layer floor (m). **Two contracts, picked by !! the coordinate family — see `rdb_vcoord :: vcoord_h_min_role`.** !! !! On the GEOMETRIC families (`VCOORD_ZSTAR_FULL`, `VCOORD_Z_FIXED`) !! this is the thickness handed to filler layers that lie BELOW the !! local bed. They hold no water; the floor exists ONLY so !! `target_h` is never exactly zero and no h-dividing kernel can !! 1/0. They are MEANT to be classified vanished downstream, so the !! default sits deliberately BELOW the D4 skip/merge marker !! `H_VANISHED = 1.5e-4` — not by accident, and not a floor in the !! `angstrom_h` sense (the D4 taxonomy forbids using `H_VANISHED` !! as a positivity floor). Thinner is also better physics here: !! each filler interface carries the full topographic slope, so the !! spurious rest PGF transport it drives scales WITH the floor (the !! same argument that took `seed_h_layer_uniform_z_impl` off !! `2*H_VANISHED`). `validate_config` warns on a value above !! `H_VANISHED` under these families, and refuses a non-positive one. !! !! On the DENSITY families (`VCOORD_RHO`, `VCOORD_HYCOM`) the !! collapsed layers are real layers the inversion squeezed shut !! anywhere in the column; they carry tracer mass, so !! `compute_target_h_rho_impl` inflates them to !! `max(zstar_h_min, 2*H_VANISHED)` to keep them above the remap !! drain. There `zstar_h_min` is additionally the pre-compaction !! strip threshold, so a large value is meaningful rather than wrong. integer :: zstar_n_surf = 0 !! Number of fine near-surface layers for ZSTAR_FULL. ≤ 0 = !! auto-pick (max(1, nz_ml/3)). integer :: zstar_stretching = STRETCH_UNIFORM !! Stretching mode for ZSTAR_FULL. `STRETCH_UNIFORM` (default) !! or `STRETCH_LOG` for a geometric near-surface fine zone. real(wp) :: zsigma_depth_transition = 200.0_wp !! Sigma → z* transition depth (m) for `VCOORD_ZSTAR_SIGMA`. real(wp) :: zsigma_blend_width = 100.0_wp !! Smoothstep blend width (m) above the transition depth. real(wp) :: rho_ref_pressure = 2.0e7_wp !! Reference pressure (Pa, default 2e7 = 2000 dbar) for the !! potential density that defines the `VCOORD_RHO` coordinate. !! A rdb convention (not MOM6-inherited). real(wp) :: z_fixed_h_ref = 0.0_wp !! Total reference depth (m) for `VCOORD_Z_FIXED` and !! `VCOORD_ZSTAR` (the MOM6 z* nominal profile is the `z_fixed` !! one, dilated per column — `ocean_vcoord_zstar_target`). Layer !! interfaces sit at `z = k · h_ref / nz_ml` from the surface, !! same as MOM6's `COORD_CONFIG = "gprime"` with `MAXIMUM_DEPTH !! = h_ref`. Driver writes from `cfg%ocean%topo%max_depth` at init. !! When 0 (default) the `compute_target_h` Z_FIXED branch falls !! back to a uniform `H · dsig(k)` target so the path stays !! sane in tests that don't explicitly set this knob. !! Under a stretched profile (`z_fixed_use_profile`) it is the !! profile's total depth, `z_fixed_zi(0)`, and the nominal !! interfaces come from `z_fixed_zi` instead of `h_ref/nz`. logical :: z_fixed_use_profile = .false. !! `&vcoord_nml z_fixed_profile /= "uniform"`: the `VCOORD_Z_FIXED` !! (and `VCOORD_ZSTAR`) nominal interfaces come from `z_fixed_zi` / !! `z_fixed_dz` (set by `ocean_vcoord_set_z_fixed_profile`) rather than from the uniform !! `z_fixed_h_ref/nz_ml`. Scalar, rides `copyin(this)`. Default !! `.false.` ⇒ the uniform arithmetic, byte-identical. real(wp), allocatable :: z_fixed_zi(:) !! `VCOORD_Z_FIXED` nominal interface depths (m, positive down, !! below `z = 0`), shape `0:nz_ml`, BOTTOM-UP like the state: !! `z_fixed_zi(k)` is the TOP interface of layer `k`, so !! `z_fixed_zi(nz_ml) = 0` (the surface) and `z_fixed_zi(0)` is the !! profile's total depth. Allocated at init (zeros), so it is never !! a placeholder; read only when `z_fixed_use_profile`. real(wp), allocatable :: z_fixed_dz(:) !! `VCOORD_Z_FIXED` nominal layer thicknesses (m), shape `nz_ml`, !! bottom-up: `z_fixed_zi(k-1) - z_fixed_zi(k)`, stored separately !! so the partial-top-cell threshold uses the exact namelist !! value. Read only when `z_fixed_use_profile`. real(wp) :: regrid_time_scale = 0.0_wp !! Grid time-filter timescale τ (s) for the ALE regrid. After !! `compute_target_h` builds the new target grid, the remap step !! relaxes the coordinate a fraction `dt/(τ+dt)` toward that !! target each outer step rather than jumping to it — damping the !! per-step grid-motion shock that drives the σ/z* PGE !! (White & Adcroft 2008, the grid time-filter). Scalar on the !! type, reaches the device through the existing `copyin(this)`; !! no new device array. Default `0.0` ⇒ `wtd = 1` ⇒ jump to !! target ⇒ bit-identical to the no-filter remap. logical :: remap_boundary_extrap = .false. !! Close the ALE remap's reconstruction at the two boundary cells !! (`k=1`, `k=nz`) with the linear-exact one-sided edge pair !! instead of the PCM flatten (MOM6 `BOUNDARY_EXTRAPOLATION`). !! !! The default closure makes PLM/PPM/PPM_H4/PQM first-order in !! exactly the two cells adjacent to the bed and the surface, so !! a column whose tracer is linear in z is remapped with an O(h) !! error there every thermo step. Under a terrain-following !! coordinate over a slope that error differs between neighbouring !! columns, which is a horizontal density gradient, which is a !! spurious pressure-gradient force — and with rotation it feeds a !! growing grid mode trapped in those same layers (see !! `docs/CAPABILITIES_AND_LIMITATIONS.md`). Scalar on the type, !! reaches the device through the existing `copyin(this)`; no new !! device array. Default `.false.` ⇒ bit-identical. logical :: remap_nonuniform_weights = .false. !! Use the non-uniform-grid reconstruction weights in the ALE !! remap's PLM slope and PPM edge estimate — Colella & Woodward !! (1984) eqs (1.6)-(1.8) — instead of their equal-thickness !! specialisations (`0.5·minmod` and `(7/12, -1/12)`). !! !! The shipped formulae are linear-exact only when the SOURCE !! column is uniform, which under every geometric family but !! `sigma`-on-flat-bed it is not: a stretched column carries an !! O(Δh/h) reconstruction error on a profile linear in z, in the !! whole interior rather than only at the two boundary cells !! `remap_boundary_extrap` addresses. The two knobs are !! complementary — the interior needs this one, the outermost two !! cells need that one, and a column is exact only with BOTH. !! PPM_H4 and PQM carry thickness-weighted stencils already and are !! unaffected (their small-`nz` fallbacks excepted). Scalar on the !! type, reaches the device through the existing `copyin(this)`; no !! new device array. Default `.false.` ⇒ bit-identical. logical :: remap_check_preconditions = .false. !! Assert the ALE remap's column preconditions once per remap and !! fail loud on a violation (audit findings V5, V6). !! !! The overlap sweep every reconstruction shares assumes both !! `dz >= 0` (a negative source thickness makes the cumulative !! interface stack NON-MONOTONE, and the sweep then integrates the !! reversed interval twice — creating mass with no NaN and no bounds !! hit) and `sum(dz_old) == sum(dz_new)` (a short target silently !! deletes the non-overlapping tail; a long one integrates it as !! `q = 0`). Neither has ever been checked, and the target builders !! break the second one on degenerate columns. Diagnostic — a !! per-column reduction at the THERMO cadence, two scalars back to !! the host. Default `.false.` ⇒ the check never runs. logical :: remap_vel_conserve_ke = .false. !! Enable the KE-conserving rescale of the remapped layer !! velocities. After the per-face column remap (which already !! conserves `u·h`, i.e. momentum), rescale the BAROCLINIC !! velocity anomaly per column so column KE `Σ ½ h·u²` is !! preserved (Adcroft & Hallberg 2006 layer-velocity remap), !! capped at a 1.25× rescale factor. The barotropic/depth-mean !! component is never touched (mode-split consistency). Default !! `.false.` ⇒ velocities unchanged ⇒ bit-identical. logical :: check_vanished_content = .false. !! `&vcoord_nml check_vanished_content` — the I1′ tripwire. Carried !! on this slot (rather than on `ocean_dyn_t`) because the vertical !! coordinate is what MAKES vanished layers, so the knob that !! polices them belongs beside `zstar_h_min` and the filler !! contract. A plain scalar on the type: it rides the existing !! `copyin(this)` and adds no device array. Read by !! `check_vanished_invariant_or_die` in `rdb_ocean_dyn`. Default !! `.false.` ⇒ no scan, no cost. logical :: zfixed_closed_faces = .false. !! `&vcoord_nml zfixed_closed_faces` — partial-step z-level face !! closure. Meaningful on the three GEOMETRIC families that !! vanish bed-side layers, `VCOORD_Z_FIXED`, `VCOORD_ZSTAR` (MOM6 !! z*) and `VCOORD_ZSTAR_FULL`, where a layer whose reference range lies !! inside the bed (or, under `z_fixed`, the ice draft) is an !! inert FILLER; a velocity face at which that layer is a filler !! on EITHER side is a z-level WALL, not a thin passage (Adcroft, !! Hill & Marshall 1997; Losch 2008). !! !! The per-layer 0/1 face mask itself lives on `ocean_metrics_t` !! (`open_u`/`open_v`, built once at configure by !! `ocean_vcoord_closed_face_masks` from THIS module's target at !! `eta = 0`, `ocean_vcoord_eta0_target`). The flag is carried here so the ALE !! remap driver — which never sees `ocean_metrics_t` — can build !! its FACE columns as `min(h_L, h_R)` and drop the closed !! layers, instead of pouring momentum into water that is not !! there. Scalar on the type, reaches the device through the !! existing `copyin(this)`. Default `.false.` => bit-identical. ! ---- Cached extents (for kernel loops + sanity checks) ---- integer :: nx_total = 0 !! Total i-extent of `target_h` (incl. halos). integer :: ny_total = 0 !! Total j-extent of `target_h` (incl. halos). integer :: nz_ml = 0 !! Number of active layers. contains procedure, non_overridable :: init => ocean_vcoord_init procedure, non_overridable :: destroy => ocean_vcoord_destroy procedure, non_overridable :: enter_data => ocean_vcoord_enter_data procedure, non_overridable :: exit_data => ocean_vcoord_exit_data procedure, non_overridable :: compute_target_h => ocean_vcoord_compute_target_h procedure, non_overridable :: compute_target_h_rho => ocean_vcoord_compute_target_h_rho procedure, non_overridable :: build_zref_full => ocean_vcoord_build_zref_full procedure, non_overridable :: bytes => ocean_vcoord_bytes end type ocean_vcoord_t contains subroutine ocean_vcoord_init(this, grid, nz_ml) !! Allocate every per-column array the slot owns: `dsig`, !! `z_ref_global`, `target_h`, `z_ref`. All sized once at init — !! grids don't resize. Host allocations only; `enter_data` ships !! them to the device. class(ocean_vcoord_t), intent(inout) :: this type(hgrid_t), intent(in) :: grid integer, intent(in), optional :: nz_ml integer :: nz_local, k nz_local = 1 if (present(nz_ml)) nz_local = nz_ml if (nz_local < 1) nz_local = 1 this%nx_total = grid%nx_total this%ny_total = grid%ny_total this%nz_ml = nz_local allocate (this%dsig(nz_local)) do k = 1, nz_local this%dsig(k) = 1.0_wp/real(nz_local, wp) end do ! Default reference z-interfaces: uniform 0..1 normalised. The ! ZSIGMA / ZSTAR_SIGMA branches that consume this expect absolute ! metre values from the namelist parser; the normalised default ! is only useful for VCOORD_SIGMA (which ignores it) ! and for unit tests that pre-populate before running. allocate (this%z_ref_global(0:nz_local)) do k = 0, nz_local this%z_ref_global(k) = real(k, wp)/real(nz_local, wp) end do allocate (this%target_h(grid%nx_total, grid%ny_total, nz_local), source=0.0_wp) ! Per-column z* reference table — sized but not populated. Callers ! invoke `build_zref_full(h_bed)` once at setup to fill it. Until ! then, ZSTAR_FULL's `compute_target_h` walks a column of zeros, ! which the formula safely degrades to an all-h_min vanishing- ! layer column (the dry / shallow degenerate case). allocate (this%z_ref(grid%nx_total, grid%ny_total, 0:nz_local), source=0.0_wp) ! Geopotential depth of the column top. Zero = the `z = 0` datum; ! `configure_ocean_cavity` overwrites it with the static ice draft ! when a cavity is configured. Allocated unconditionally so the ! Z_FIXED kernel never sees a placeholder-sized array. allocate (this%z_top(grid%nx_total, grid%ny_total), source=0.0_wp) ! `VCOORD_Z_FIXED` stretched nominal profile — allocated ! unconditionally (zeros) so the target kernel's explicit-shape ! dummies never see a placeholder; filled by ! `ocean_vcoord_set_z_fixed_profile` only when a profile is set. allocate (this%z_fixed_zi(0:nz_local), source=0.0_wp) allocate (this%z_fixed_dz(nz_local), source=0.0_wp) ! Isopycnal target densities — sized `0:nz_ml`, populated by the ! setup wiring only when `coord_type == VCOORD_RHO`. Default is a ! benign monotone ramp (1020..1030 kg/m³) so the slot is always ! valid; the setup path overwrites it from the namelist. allocate (this%rho_target(0:nz_local)) do k = 0, nz_local this%rho_target(k) = 1020.0_wp + 10.0_wp*real(k, wp)/real(nz_local, wp) end do ! ALE remap scratch — persistent so the remap step doesn't ! allocate fresh host buffers per call (the DCs there were ! reading device-resident `ms%h_layer` and writing into freshly- ! allocated host arrays, forcing per-iteration H↔D transfers ! that dominated wallclock in production runs). allocate (this%remap_total_h(grid%nx_total, grid%ny_total), source=0.0_wp) allocate (this%remap_h_ref(grid%nx_total, grid%ny_total), source=0.0_wp) allocate (this%remap_h_old(grid%nx_total, grid%ny_total, nz_local), source=0.0_wp) allocate (this%remap_conc_t(grid%nx_total, grid%ny_total, nz_local), source=0.0_wp) allocate (this%remap_conc_s(grid%nx_total, grid%ny_total, nz_local), source=0.0_wp) this%is_init = .true. end subroutine ocean_vcoord_init subroutine ocean_vcoord_destroy(this) class(ocean_vcoord_t), intent(inout) :: this this%is_init = .false. if (allocated(this%dsig)) deallocate (this%dsig) if (allocated(this%z_ref_global)) deallocate (this%z_ref_global) if (allocated(this%target_h)) deallocate (this%target_h) if (allocated(this%z_ref)) deallocate (this%z_ref) if (allocated(this%z_top)) deallocate (this%z_top) if (allocated(this%z_fixed_zi)) deallocate (this%z_fixed_zi) if (allocated(this%z_fixed_dz)) deallocate (this%z_fixed_dz) this%z_fixed_use_profile = .false. if (allocated(this%rho_target)) deallocate (this%rho_target) if (allocated(this%remap_total_h)) deallocate (this%remap_total_h) if (allocated(this%remap_h_ref)) deallocate (this%remap_h_ref) if (allocated(this%remap_h_old)) deallocate (this%remap_h_old) if (allocated(this%remap_conc_t)) deallocate (this%remap_conc_t) if (allocated(this%remap_conc_s)) deallocate (this%remap_conc_s) this%nx_total = 0 this%ny_total = 0 this%nz_ml = 0 end subroutine ocean_vcoord_destroy subroutine ocean_vcoord_enter_data(this) !! Map every host allocatable onto the device. Idempotent guard !! via `is_init`. class(ocean_vcoord_t), intent(inout) :: this select type (this) type is (ocean_vcoord_t) call ocean_vcoord_enter_data_impl(this) end select end subroutine ocean_vcoord_enter_data subroutine ocean_vcoord_enter_data_impl(this) type(ocean_vcoord_t), intent(inout) :: this if (.not. this%is_init) return !$acc enter data copyin(this%dsig, this%z_ref_global, this%target_h, this%z_ref) !$acc enter data copyin(this%z_top) !$acc enter data copyin(this%z_fixed_zi, this%z_fixed_dz) !$acc enter data copyin(this%rho_target) !$acc enter data copyin(this%remap_total_h, this%remap_h_ref, this%remap_h_old) !$acc enter data copyin(this%remap_conc_t, this%remap_conc_s) end subroutine ocean_vcoord_enter_data_impl subroutine ocean_vcoord_exit_data(this) class(ocean_vcoord_t), intent(inout) :: this select type (this) type is (ocean_vcoord_t) call ocean_vcoord_exit_data_impl(this) end select end subroutine ocean_vcoord_exit_data subroutine ocean_vcoord_exit_data_impl(this) type(ocean_vcoord_t), intent(inout) :: this if (.not. this%is_init) return !$acc exit data delete(this%remap_conc_t, this%remap_conc_s) !$acc exit data delete(this%remap_h_old, this%remap_h_ref, this%remap_total_h) !$acc exit data delete(this%rho_target) !$acc exit data delete(this%z_fixed_zi, this%z_fixed_dz) !$acc exit data delete(this%z_top) !$acc exit data delete(this%z_ref, this%target_h, this%z_ref_global, this%dsig) end subroutine ocean_vcoord_exit_data_impl pure subroutine ocean_vcoord_build_zref_full(this, h_bed) !! Populate `z_ref(i, j, 0:nz_ml)` per column from the local !! bathymetry `h_bed(i, j)`. Mirrors `zstar_full_build_column` !! from `src/ALE/rdb_vcoord.F90` but as a 2D loop owned by this !! slot — keeps the coastal helper untouched while letting the !! ocean path own its z_ref lifecycle. !! !! Top-down indexing inside the column: `z_ref(:, :, 0) = 0` is !! the surface, `z_ref(:, :, nz_ml) = h_bed(:, :)` is the bed. !! `compute_target_h(VCOORD_ZSTAR_FULL)` walks this table and !! emits ROMS-ordered `target_h(:, :, 1..nz_ml)`. !! !! Degenerate columns (`h_bed ≤ 0`) get a column of zeros — the !! ZSTAR_FULL branch then produces an all-h_min vanishing-layer !! result which the downstream dry-cell guards already handle. !! !! Call sites: once at setup, then any time bathymetry changes !! (which today is "never" — the ocean path doesn't move the bed). class(ocean_vcoord_t), intent(inout) :: this real(wp), intent(in) :: h_bed(:, :) !! Bed depth at cell centres (m, positive-down). integer :: n_surf_eff, n_coarse, nz, i, j, k real(wp) :: h_surf_eff, h_fine, h_coarse, dz_uniform, r, base, w real(wp) :: r_lo, r_hi, r_mid, f_mid, target_ratio integer :: it if (.not. this%is_init) return nz = this%nz_ml h_surf_eff = this%zstar_h_surf_target if (this%zstar_n_surf <= 0) then n_surf_eff = max(1, nz/3) else n_surf_eff = max(1, min(nz - 1, this%zstar_n_surf)) end if n_coarse = nz - n_surf_eff ! The per-column work is column-local — different (i, j) cells ! don't read each other. The do-concurrent `local()` clause keeps ! the scalar scratch private per thread; the bisection branch in ! the LOG stretching case requires a sequential inner loop so we ! keep it as a non-concurrent block. ! One-time z_ref setup runs on the HOST (plain do, not do concurrent): ! it writes this%z_ref, a vcoord allocatable component that is not yet ! device-mapped at seed time (build_zref_full runs before ! ocean_state_enter_data). A device kernel writing an un-present ! derived-type component faults under OpenMP-target offload (the ! component can't be implicitly mapped); stdpar tolerates it but this ! is one-time init, so host is correct and costs nothing. do j = 1, this%ny_total do i = 1, this%nx_total if (h_bed(i, j) <= 0.0_wp) then ! Degenerate / dry column — leave z_ref at zero. do k = 0, nz this%z_ref(i, j, k) = 0.0_wp end do else if (h_surf_eff <= 0.0_wp .or. nz == 1) then ! Uniform spacing fallback. this%z_ref(i, j, 0) = 0.0_wp do k = 1, nz this%z_ref(i, j, k) = h_bed(i, j)*real(k, wp)/real(nz, wp) end do else this%z_ref(i, j, 0) = 0.0_wp if (this%zstar_stretching == STRETCH_LOG & .and. n_surf_eff >= 2 .and. n_coarse >= 1) then ! Geometric fine zone: layer k thickness = h_surf · r^(k-1). ! Bisect for r so total = h_bed. target_ratio = h_bed(i, j)/h_surf_eff r_lo = 1.000001_wp r_hi = 10.0_wp do it = 1, 60 r_mid = 0.5_wp*(r_lo + r_hi) f_mid = (r_mid**n_surf_eff - 1.0_wp)/(r_mid - 1.0_wp) & + r_mid**(n_surf_eff - 1)*real(n_coarse, wp) if (f_mid > target_ratio) then r_hi = r_mid else r_lo = r_mid end if if (r_hi - r_lo < 1.0e-9_wp) exit end do r = 0.5_wp*(r_lo + r_hi) w = 1.0_wp base = 0.0_wp do k = 1, n_surf_eff base = base + h_surf_eff*w this%z_ref(i, j, k) = base w = w*r end do dz_uniform = h_surf_eff*r**(n_surf_eff - 1) do k = n_surf_eff + 1, nz this%z_ref(i, j, k) = this%z_ref(i, j, k - 1) + dz_uniform end do else ! Uniform fine zone + uniform coarse fill. h_fine = h_surf_eff*real(n_surf_eff, wp) if (h_fine >= h_bed(i, j)) then ! Column too shallow to honour h_surf for all fine layers. ! Match MOM6 isopycnal NK=2: surface absorbs the water, ! deeper layers vanish. Reserve `zstar_h_min` for each ! vanishing layer so target_h ≥ h_min downstream — kernels ! that divide by h_layer can't see exact zero (would NaN). ! Fine layers fill top-down at h_surf_eff each; bottom-most ! fine grabs the residual; coarse vanishes to h_min. base = 0.0_wp do k = 1, n_surf_eff base = base + h_surf_eff ! Leave room for h_min of every layer below this one. if (base > h_bed(i, j) - real(nz - k, wp)*this%zstar_h_min) then base = h_bed(i, j) - real(nz - k, wp)*this%zstar_h_min end if this%z_ref(i, j, k) = base end do do k = n_surf_eff + 1, nz this%z_ref(i, j, k) = this%z_ref(i, j, k - 1) + this%zstar_h_min end do else h_coarse = h_bed(i, j) - h_fine do k = 1, n_surf_eff this%z_ref(i, j, k) = h_surf_eff*real(k, wp) end do if (n_coarse > 0) then dz_uniform = h_coarse/real(n_coarse, wp) do k = n_surf_eff + 1, nz this%z_ref(i, j, k) = h_fine + dz_uniform*real(k - n_surf_eff, wp) end do end if end if end if ! Pin the bed interface to h_bed exactly (round-off guard). this%z_ref(i, j, nz) = h_bed(i, j) end if end do end do end subroutine ocean_vcoord_build_zref_full pure subroutine ocean_vcoord_compute_target_h(this, total_h, eta) !! Thin polymorphic wrapper. A type-bound procedure's passed object !! must be `class(...)`, but mapping a polymorphic list item into a !! `target` / offload region is unspecified behaviour (gfortran !! `-Wopenmp`; see FORTRAN_STYLE.md) — so the `do concurrent` kernels !! live in the `type(ocean_vcoord_t)` `_impl` and this wrapper only !! resolves the concrete type. `ocean_vcoord_t` is never extended, !! so the dynamic type is always the declared type. Mirrors the !! enter_data/exit_data split. class(ocean_vcoord_t), intent(inout) :: this real(wp), intent(in) :: total_h(:, :) ! assumed-shape-ok: thin TBP wrapper, no kernel real(wp), intent(in) :: eta(:, :) ! assumed-shape-ok: thin TBP wrapper, no kernel select type (this) type is (ocean_vcoord_t) call ocean_vcoord_compute_target_h_impl(this, total_h, eta) end select end subroutine ocean_vcoord_compute_target_h pure subroutine ocean_vcoord_compute_target_h_impl(this, total_h, eta) !! Populate `target_h(i,j,k)` from the column-total depth (`H`, !! constant per column for the ocean path; coastal would pass !! the bathymetry) and the free-surface anomaly `η`. !! !! Cases: !! !! VCOORD_EULERIAN_Z — `target_h(:,:,k) = H(:,:) · dsig(k)`. !! The reference grid that `rdb_ocean_vertical_advection` !! pins `h_layer` to via its cancellation trick. No η !! dependence — that's the whole point of the mode. !! !! VCOORD_SIGMA — `target_h(:,:,k) = (H + η) · dsig(k)`. !! Pure terrain-following. !! !! VCOORD_ZSTAR — MOM6 z* (`ocean_vcoord_zstar_target`): the !! `z_fixed` nominal profile dilated per column by the free-surface !! stretching, over a partial bed cell and inert bed fillers. !! !! VCOORD_ZSIGMA / VCOORD_ZSTAR_SIGMA — smoothstep blends, see !! module head comment. !! !! VCOORD_ZSTAR_FULL — walks the per-column `z_ref(:, :, 0:nz)` !! table populated by `build_zref_full`. Surface layer absorbs !! η; subsurface layers keep their reference thicknesses !! when the column is at or above reference depth. When the !! column is shallower than reference (η < 0), bed-side !! layers vanish to `zstar_h_min` and the surface trim makes !! sum(target_h) = H exactly. Mirrors coastal !! `vcoord_target_dz_column_zstar_full`. !! !! All cases preserve `sum_k target_h(i,j,k) = H + η` (or `= H` !! for EULERIAN_Z). The remap kernel will rely on this. type(ocean_vcoord_t), intent(inout) :: this ! assumed-shape-ok: cadence-bounded (once per outer ALE step); flat impl ! behind the thin ocean_vcoord_compute_target_h class wrapper below. real(wp), intent(in) :: total_h(:, :) !! Column-total depth H(i, j) (m). real(wp), intent(in) :: eta(:, :) ! assumed-shape-ok: same reason as total_h above !! Free-surface anomaly η(i, j) (m). real(wp) :: h_nominal if (.not. this%is_init) return ! Lagrangian / isopycnal: the target IS the current h_layer — ! nothing to compute, the ALE remap step is a no-op (see ! `ocean_apply_ale_remap_step`). Return before touching ! `target_h` so the caller keeps the live `h_layer`. if (this%coord_type == VCOORD_LAGRANGIAN) return ! VCOORD_LAGRANGIAN: target is whatever `h_layer` already is. The ! ALE remap is a no-op for this case (see `ocean_apply_ale_remap_step`), ! so `target_h` is never read; we can skip the compute entirely to ! avoid burning a kernel launch. if (this%coord_type == VCOORD_LAGRANGIAN) return ! Every loop runs in a flat kernel on explicit-shape / scalar dummies: ! no `this%` component and no `associate`-name reaches a ! `do concurrent` (ifx's do-concurrent -> OpenMP-target lowering ICEs ! on the former; the latter is the shape CLAUDE.md forbids — ifx ! evaluates such names as zero, and nvfortran -stdpar=gpu hands a ! by-reference one to a device callee as a HOST address, see ! `ocean_vcoord_rho_target`). if (this%coord_type == VCOORD_Z_FIXED) then ! Fixed-z (quasi-geopotential) interfaces with vanishing ! layers at BOTH ends: `h_min` fillers below the bed and — ! under a rigid top (`z_top > 0`, an ice-shelf cavity) — ! `h_min` fillers inside the ice, with a partial cell at each ! live end. See `ocean_vcoord_z_fixed_target`, which owns the ! algorithm and is shared with the initial-thickness seed. ! ! When `z_fixed_h_ref = 0` (knob unset) fall back to uniform ! `(H + η) · dsig(k)` so tests that omit the knob still get ! something sensible: the SIGMA branch of the geometric kernel ! evaluates exactly that expression. ! ! A stretched nominal profile (`z_fixed_use_profile`) replaces ! `h_nominal` with the per-layer `z_fixed_zi` / `z_fixed_dz` ! tables; the uniform path below is untouched. h_nominal = 0.0_wp if (this%z_fixed_h_ref > 0.0_wp) then h_nominal = this%z_fixed_h_ref/real(this%nz_ml, wp) end if if (this%z_fixed_use_profile .or. h_nominal > 0.0_wp) then call ocean_vcoord_z_fixed_target(this%target_h, total_h, eta, this%z_top, & this%nx_total, this%ny_total, this%nz_ml, & h_nominal, this%z_fixed_use_profile, & this%z_fixed_zi, this%z_fixed_dz, & this%zstar_h_min) return end if call ocean_vcoord_geometric_target(VCOORD_SIGMA, this%nx_total, this%ny_total, & this%nz_ml, this%target_h, total_h, eta, & this%dsig, this%z_ref_global, this%z_ref, & this%zsigma_depth_transition, & this%zsigma_blend_width, this%zstar_h_min) return end if if (this%coord_type == VCOORD_ZSTAR) then ! MOM6 z*: the `z_fixed` nominal profile (uniform `z_fixed_h_ref/nz` ! or the stretched `z_fixed_zi` table) dilated by the column's ! free-surface stretching. Same `z_fixed_h_ref = 0` fallback as ! `z_fixed` (uniform sigma) for a slot no profile was resolved on. h_nominal = 0.0_wp if (this%z_fixed_h_ref > 0.0_wp) then h_nominal = this%z_fixed_h_ref/real(this%nz_ml, wp) end if if (this%z_fixed_use_profile .or. h_nominal > 0.0_wp) then call ocean_vcoord_zstar_target(this%target_h, total_h, eta, & this%nx_total, this%ny_total, this%nz_ml, & h_nominal, this%z_fixed_use_profile, & this%z_fixed_zi, this%zstar_h_min) return end if call ocean_vcoord_geometric_target(VCOORD_SIGMA, this%nx_total, this%ny_total, & this%nz_ml, this%target_h, total_h, eta, & this%dsig, this%z_ref_global, this%z_ref, & this%zsigma_depth_transition, & this%zsigma_blend_width, this%zstar_h_min) return end if call ocean_vcoord_geometric_target(this%coord_type, this%nx_total, this%ny_total, & this%nz_ml, this%target_h, total_h, eta, & this%dsig, this%z_ref_global, this%z_ref, & this%zsigma_depth_transition, & this%zsigma_blend_width, this%zstar_h_min) end subroutine ocean_vcoord_compute_target_h_impl pure subroutine ocean_vcoord_geometric_target(coord_type, nx, ny, nz, target_h, & total_h, eta, dsig, z_ref_global, & z_ref, zsigma_depth_transition, & zsigma_blend_width, zstar_h_min) !! Geometric target-grid kernels (EULERIAN_Z, SIGMA, ZSIGMA, !! ZSTAR_SIGMA, ZSTAR_FULL; formulae documented on !! `ocean_vcoord_compute_target_h_impl`). Flat on purpose — every !! array an explicit-shape dummy, every knob a scalar dummy — so no !! derived-type component and no `associate`-name reaches a !! `do concurrent` (see `ocean_vcoord_rho_target` for the GPU fault !! that shape caused there). integer, intent(in), value :: coord_type !! `VCOORD_*` family (not LAGRANGIAN / Z_FIXED / ZSTAR: the !! dispatcher owns those). integer, intent(in), value :: nx !! i-extent of every horizontal array (total, incl. halos). integer, intent(in), value :: ny !! j-extent of every horizontal array (total, incl. halos). integer, intent(in), value :: nz !! Number of layers; `k = 1` is the bed, `k = nz` the surface. real(wp), intent(inout) :: target_h(nx, ny, nz) !! Target layer thickness (m), bottom-up. real(wp), intent(in) :: total_h(nx, ny) !! Column-total depth H (m). real(wp), intent(in) :: eta(nx, ny) !! Free-surface anomaly η (m). real(wp), intent(in) :: dsig(nz) !! Nominal layer fractions, bottom-up. real(wp), intent(in) :: z_ref_global(0:nz) !! Global reference interface depths (ZSIGMA / ZSTAR_SIGMA). real(wp), intent(in) :: z_ref(nx, ny, 0:nz) !! Per-column reference interface depths (ZSTAR_FULL). real(wp), intent(in), value :: zsigma_depth_transition !! Sigma → z* transition depth (m). real(wp), intent(in), value :: zsigma_blend_width !! Smoothstep blend width (m). real(wp), intent(in), value :: zstar_h_min !! Vanished-layer thickness (m). integer :: i, j, k real(wp) :: column_total, alpha, x, z_top_k, z_bot_k, dz_z, dz_sum, deficit real(wp) :: z_ref_nz_inv real(wp) :: h_bed_ref, eta_loc, H_eff, z_upper, z_lower, sum_dz select case (coord_type) case (VCOORD_EULERIAN_Z) do concurrent(k=1:nz, j=1:ny, i=1:nx) target_h(i, j, k) = total_h(i, j)*dsig(k) end do case (VCOORD_SIGMA) do concurrent(k=1:nz, j=1:ny, i=1:nx) column_total = total_h(i, j) + eta(i, j) target_h(i, j, k) = column_total*dsig(k) end do case (VCOORD_ZSIGMA) ! Smoothstep blend: sigma in shallow, fixed z-levels in deep. ! Mirrors `vcoord_target_dz_column` in `src/ALE/rdb_vcoord.F90` ! but built directly on the 2D (H, η) fields. The deep-branch ! z-level intervals are clipped to the local column total so ! sum_k target_h = H + η exactly even when the column is shallower ! than the deepest reference interface; any residual deficit is ! deposited in the bed-side layer (k=1) to preserve the sum. do concurrent(j=1:ny, i=1:nx) & local(column_total, alpha, x, k, z_top_k, z_bot_k, dz_z, dz_sum, deficit) column_total = total_h(i, j) + eta(i, j) if (column_total <= zsigma_depth_transition) then do k = 1, nz target_h(i, j, k) = dsig(k)*column_total end do else if (zsigma_blend_width > 0.0_wp) then x = (column_total - zsigma_depth_transition)/zsigma_blend_width x = max(0.0_wp, min(1.0_wp, x)) alpha = x*x*(3.0_wp - 2.0_wp*x) else alpha = 1.0_wp end if dz_sum = 0.0_wp do k = 1, nz z_top_k = min(z_ref_global(nz - k), column_total) z_bot_k = min(z_ref_global(nz - k + 1), column_total) dz_z = max(z_bot_k - z_top_k, 0.0_wp) target_h(i, j, k) = (1.0_wp - alpha)*dsig(k)*column_total & + alpha*dz_z dz_sum = dz_sum + target_h(i, j, k) end do deficit = column_total - dz_sum target_h(i, j, 1) = target_h(i, j, 1) + deficit end if end do case (VCOORD_ZSTAR_SIGMA) z_ref_nz_inv = 0.0_wp if (z_ref_global(nz) > 0.0_wp) then z_ref_nz_inv = 1.0_wp/z_ref_global(nz) end if do concurrent(j=1:ny, i=1:nx) & local(column_total, alpha, x, k, dz_z) column_total = total_h(i, j) + eta(i, j) if (column_total <= zsigma_depth_transition .or. z_ref_nz_inv == 0.0_wp) then do k = 1, nz target_h(i, j, k) = dsig(k)*column_total end do else if (zsigma_blend_width > 0.0_wp) then x = (column_total - zsigma_depth_transition)/zsigma_blend_width x = max(0.0_wp, min(1.0_wp, x)) alpha = x*x*(3.0_wp - 2.0_wp*x) else alpha = 1.0_wp end if do k = 1, nz dz_z = (z_ref_global(nz - k + 1) - z_ref_global(nz - k)) & *column_total*z_ref_nz_inv target_h(i, j, k) = (1.0_wp - alpha)*dsig(k)*column_total & + alpha*dz_z end do end if end do case (VCOORD_ZSTAR_FULL) ! Per-column z*-full: walk the cached `z_ref(i, j, 0:nz)` from ! `build_zref_full`. Surface layer absorbs η when η ≥ 0; ! bed-side layers vanish to `zstar_h_min` and the surface gets ! trimmed for exact conservation when η < 0. Mirrors ! `vcoord_target_dz_column_zstar_full` (coastal) but emits ! ROMS order (k=1 bed, k=nz surface) directly. do concurrent(j=1:ny, i=1:nx) & local(k, h_bed_ref, eta_loc, H_eff, z_upper, z_lower, sum_dz, deficit) h_bed_ref = z_ref(i, j, nz) eta_loc = (total_h(i, j) + eta(i, j)) - h_bed_ref H_eff = max(total_h(i, j) + eta(i, j), 0.0_wp) if (h_bed_ref <= 0.0_wp) then ! Degenerate column: emit a single vanishing-layer stack. do k = 1, nz target_h(i, j, k) = zstar_h_min end do else if (eta_loc >= 0.0_wp) then ! Column at or above reference: subsurface = z_ref intervals, ! surface (k=nz) gets the +η. do k = 1, nz ! k_top = nz - k + 1 in the top-down z_ref convention. target_h(i, j, k) = max( & z_ref(i, j, nz - k + 1) - z_ref(i, j, nz - k), & 0.0_wp) end do target_h(i, j, nz) = target_h(i, j, nz) + eta_loc else ! Column shallower than reference (η < 0). Walk top-down, ! clip layers to H_eff, vanish below. do k = 1, nz z_upper = z_ref(i, j, nz - k) z_lower = z_ref(i, j, nz - k + 1) if (z_lower <= H_eff) then target_h(i, j, k) = z_lower - z_upper else if (z_upper < H_eff) then target_h(i, j, k) = H_eff - z_upper else target_h(i, j, k) = zstar_h_min end if end do ! Surface trim: drop the vanishing-layer overhead from the ! surface to make sum = H exactly. If the surface would ! itself fall below h_min, leave it at h_min and let the ! downstream dry-cell guards handle the deficit. sum_dz = 0.0_wp do k = 1, nz sum_dz = sum_dz + target_h(i, j, k) end do deficit = sum_dz - H_eff if (deficit > 0.0_wp) then if (target_h(i, j, nz) - deficit >= zstar_h_min) then target_h(i, j, nz) = target_h(i, j, nz) - deficit else target_h(i, j, nz) = zstar_h_min end if end if end if end do case default error stop "ocean_vcoord_geometric_target: unknown coord_type." end select end subroutine ocean_vcoord_geometric_target pure subroutine ocean_vcoord_set_z_fixed_profile(this, dz_surface_first) !! Install a stretched `VCOORD_Z_FIXED` nominal profile: flip the !! surface-first thicknesses into the bottom-up `z_fixed_dz`, build !! the interface table `z_fixed_zi` by accumulating from the surface !! (`z_fixed_zi(nz) = 0` exactly), set `z_fixed_h_ref` to the total !! and raise `z_fixed_use_profile`. Setup-time host code; must run !! BEFORE `enter_data` (the copyin captures the tables). type(ocean_vcoord_t), intent(inout) :: this real(wp), intent(in) :: dz_surface_first(:) !! Nominal thicknesses (m), `dz_surface_first(1)` = top layer; !! size must be `nz_ml`. integer :: k, nz if (.not. this%is_init) return nz = this%nz_ml if (size(dz_surface_first) /= nz) return do k = 1, nz this%z_fixed_dz(k) = dz_surface_first(nz - k + 1) end do this%z_fixed_zi(nz) = 0.0_wp do k = nz, 1, -1 this%z_fixed_zi(k - 1) = this%z_fixed_zi(k) + this%z_fixed_dz(k) end do this%z_fixed_h_ref = this%z_fixed_zi(0) this%z_fixed_use_profile = .true. end subroutine ocean_vcoord_set_z_fixed_profile pure subroutine ocean_vcoord_z_fixed_target_uniform(target_h, total_h, eta, z_top, & nx, ny, nz, h_nominal, h_min) !! `ocean_vcoord_z_fixed_target` on the UNIFORM nominal spacing !! `h_nominal` — the historical signature, for callers (tests, !! setup code without a vcoord slot) that have no profile tables. !! Same kernel, `use_profile = .false.`, so the arithmetic is the !! uniform branch's exactly; the two tables are never read. integer, intent(in) :: nx !! i-extent of every array (total, incl. halos). integer, intent(in) :: ny !! j-extent of every array (total, incl. halos). integer, intent(in) :: nz !! Number of layers; `k = 1` is the bed, `k = nz` the top. real(wp), intent(out) :: target_h(nx, ny, nz) !! Target layer thickness (m). real(wp), intent(in) :: total_h(nx, ny) !! Column reference thickness `H` (m). real(wp), intent(in) :: eta(nx, ny) !! Free-surface anomaly `η` (m). real(wp), intent(in) :: z_top(nx, ny) !! Geopotential depth of the column top (m, positive down). real(wp), intent(in) :: h_nominal !! Nominal layer spacing (m), `> 0`. real(wp), intent(in) :: h_min !! Inert-filler thickness. real(wp) :: zi_unused(0:nz), dz_unused(nz) zi_unused = 0.0_wp dz_unused = 0.0_wp call ocean_vcoord_z_fixed_target(target_h, total_h, eta, z_top, nx, ny, nz, & h_nominal, .false., zi_unused, dz_unused, h_min) end subroutine ocean_vcoord_z_fixed_target_uniform pure subroutine ocean_vcoord_z_fixed_target(target_h, total_h, eta, z_top, & nx, ny, nz, h_nominal, use_profile, & zi, dz_nom, h_min) !! `VCOORD_Z_FIXED` target grid — quasi-geopotential interfaces !! under a rigid top, with inert fillers and a partial cell at BOTH !! ends (Yung, Hallberg, Adcroft & Morrison 2026, JAMES, Fig. 1b: !! quasi-z layers are geopotential and VANISH where they outcrop !! into the ice base; Asay-Davis et al. 2016 §3.1.5: z-level models !! use both partial top and bottom cells). !! !! ### The rule !! !! Nominal interface depths are GEOPOTENTIAL and unchanged by the !! rigid top: the interface above layer `k` sits at depth !! `(nz − k)·h_nominal − η` below `z = 0`, i.e. at !! `(nz − k)·h_nominal − z_top` below the column TOP, which is !! itself at depth `z_top − η`. The walk is bed-up (`k = 1` is the !! bed) in "depth below the column top", `z_below_loc` tracking the !! bottom interface of the layer being laid: !! !! * **bed side.** A layer whose nominal top interface is deeper !! than the remaining column collapses to the inert filler !! `h_min` and hands its water UP; the lowest live layer is the !! partial BOTTOM cell and absorbs `η` plus the bed fillers' !! `h_min` debt. A cut that would leave the partial bottom !! cell at or below `Z_FIXED_BED_PARTIAL_MIN` collapses it too !! — see that constant: no LIVE thickness may land in the !! `(h_min, H_VANISHED]` band. !! * **top side, new.** A layer whose nominal range lies entirely !! above the column top — i.e. inside the ice — collapses to !! `h_min` and stacks immediately under the ice base. The !! layer whose nominal range STRADDLES the column top is the !! partial TOP cell: its thickness is the part of its nominal !! range below the top, less the top fillers' `h_min` debt. !! If that cut would leave less than !! `max(h_min, Z_FIXED_TOP_PARTIAL_FRAC*h_nominal)`, the sliver !! is merged into the layer BELOW and the index becomes another !! filler (`k_live_top` moves down one) — the mirror of the !! bed's own sliver rule. !! !! ### Stretched nominal profile !! !! With `use_profile` the nominal interface above layer `k` sits at !! `zi(k)` (bottom-up table, `zi(nz) = 0`) instead of !! `(nz − k)·h_nominal`, its bottom interface at `zi(k − 1)`, and the !! partial-top threshold is the straddling layer's OWN nominal !! thickness, `max(h_min, Z_FIXED_TOP_PARTIAL_FRAC*dz_nom(k))`. !! Everything else — the bed rule, the top filler debt, the closing !! rule — is shared. Each branch keeps its whole expression !! (`a*b − c` stays one expression in the uniform branch) so an !! FMA-contracting build contracts exactly what it did before. !! !! `Σ_k target_h = total_h + eta` is closed by construction: every !! branch assigns exactly what it subtracts from `z_below_loc`, so !! each END pays for its own fillers and there is ONE closing rule !! that works when both ends vanish. (A degenerate column thinner !! than `nz·h_min` overshoots to `nz·h_min`, exactly as the !! pre-cavity code did.) !! !! ### Bit-identity !! !! `z_top ≡ 0` ⇒ `k_live_top = nz`, so the top branch is taken ONLY !! at `k = nz`, where it evaluates `real(0, wp)*h_min = +0.0` — the !! same value as the pre-cavity `real(nz − nz, wp)*h_nominal`. Every !! other layer evaluates `real(nz − k, wp)*h_nominal − 0.0_wp`, which !! is bit-for-bit `real(nz − k, wp)*h_nominal` in IEEE round-to- !! nearest (and under FMA contraction, since `fma(a, b, −0.0)` !! rounds `a·b` once). Pinned by !! `tests/test_ocean_vcoord_zfixed_cavity.F90`. integer, intent(in) :: nx !! i-extent of every array (total, incl. halos). integer, intent(in) :: ny !! j-extent of every array (total, incl. halos). integer, intent(in) :: nz !! Number of layers; `k = 1` is the bed, `k = nz` the top. real(wp), intent(out) :: target_h(nx, ny, nz) !! Target layer thickness (m). real(wp), intent(in) :: total_h(nx, ny) !! Column reference thickness `H` (m) — `Σ h_layer − η`. real(wp), intent(in) :: eta(nx, ny) !! Free-surface anomaly `η` (m). real(wp), intent(in) :: z_top(nx, ny) !! Geopotential depth of the column top (m, positive down, !! `>= 0`). `0` ⇒ the pre-cavity arithmetic, bit-for-bit. real(wp), intent(in) :: h_nominal !! Nominal layer spacing `z_fixed_h_ref/nz` (m), `> 0`. Unused !! when `use_profile`. logical, intent(in) :: use_profile !! Take the nominal interfaces from `zi` / `dz_nom`. real(wp), intent(in) :: zi(0:nz) !! Nominal interface depths (m), bottom-up, `zi(k)` = top of !! layer `k`, `zi(nz) = 0`. Read only when `use_profile`. real(wp), intent(in) :: dz_nom(nz) !! Nominal layer thicknesses (m), bottom-up. Read only when !! `use_profile`. real(wp), intent(in) :: h_min !! Inert-filler thickness (`zstar_h_min`, `<= H_VANISHED`). integer :: i, j, k, k_live_top real(wp) :: z_below_loc, z_above_nominal_loc, z_top_loc, partial_min real(wp) :: bed_partial_min logical :: vanish_loc partial_min = max(h_min, Z_FIXED_TOP_PARTIAL_FRAC*h_nominal) bed_partial_min = max(h_min, Z_FIXED_BED_PARTIAL_MIN) do concurrent(j=1:ny, i=1:nx) & local(k, k_live_top, z_below_loc, z_above_nominal_loc, z_top_loc, vanish_loc) z_top_loc = z_top(i, j) ! Index of the shallowest layer the rigid top leaves live: the ! largest k whose nominal BOTTOM interface, at depth ! `(nz − k + 1)*h_nominal` below z = 0, clears the top by more ! than the minimum partial-cell thickness. Monotone in k, so ! the last k that passes is the largest. Skipped entirely when ! there is no rigid top (the bit-identity gate). k_live_top = nz if (z_top_loc > 0.0_wp) then k_live_top = 1 do k = 1, nz if (use_profile) then if (zi(k - 1) - z_top_loc > & max(h_min, Z_FIXED_TOP_PARTIAL_FRAC*dz_nom(k))) then k_live_top = k end if else if (real(nz - k + 1, wp)*h_nominal - z_top_loc > partial_min) then k_live_top = k end if end if end do end if z_below_loc = total_h(i, j) + eta(i, j) do k = 1, nz if (k >= k_live_top) then ! At and above the first live layer the stack is measured ! from the column TOP: `(nz − k)*h_min` is the debt the ! top fillers above this layer still owe, so the partial ! top cell (k = k_live_top) is cut at the ice base and ! pays for them, and every filler above lands on h_min. z_above_nominal_loc = real(nz - k, wp)*h_min vanish_loc = z_above_nominal_loc > z_below_loc - h_min else if (use_profile) then z_above_nominal_loc = zi(k) - z_top_loc else z_above_nominal_loc = real(nz - k, wp)*h_nominal - z_top_loc end if ! Bed side: a LIVE partial bottom cell must clear the ! REMAP's vanish marker, not merely `h_min`. `>=` (not ! `>`) because the marker itself reads as vanished — the ! strict-`>` convention of `H_VANISHED`. vanish_loc = z_above_nominal_loc >= z_below_loc - bed_partial_min end if if (vanish_loc) then ! Below the bed / would be a sub-marker sliver — vanish. target_h(i, j, k) = h_min z_below_loc = z_below_loc - h_min else ! Layer fits — nominal thickness, or the end residual. target_h(i, j, k) = z_below_loc - z_above_nominal_loc z_below_loc = z_above_nominal_loc end if end do end do end subroutine ocean_vcoord_z_fixed_target pure subroutine ocean_vcoord_zstar_target(target_h, total_h, eta, nx, ny, nz, & h_nominal, use_profile, zi, h_min) !! `VCOORD_ZSTAR` target grid — MOM6 z* (`REGRIDDING_COORDINATE_MODE !! = "Z*"`, `build_zstar_column`, MOM6 `src/ALE/coord_zlike.F90` !! lines 65-146): the FIXED nominal z profile of `z_fixed` !! (`&vcoord_nml z_fixed_profile` uniform / list / tanh — the same !! `z_fixed_zi` table, not a parallel one) DILATED per column by the !! free-surface stretching, over a partial bed cell and inert bed !! fillers. !! !! ### MOM6 !! !! Without a rigid top MOM6 sets `stretching = (H+η)/H` !! (coord_zlike.F90:109), lays the interfaces top-down from the free !! surface at `z_k = η − stretching·Z_k` (`Z_k` the nominal depth, !! l.130-134), pins the bottom interface at `−H` and clamps upward so !! no layer is thinner than `min_thickness` (l.138-143). The height !! of nominal interface `k` above the bed is therefore !! `stretching·(H − Z_k)`: a layer lies below the bed iff its nominal !! top `Z_k ≥ H`, WHATEVER `η` is (`stretching > 0`). The live/filler !! pattern is exactly static in `η`, of either sign. !! !! ### This kernel !! !! Two passes per column, both the `z_fixed` bed walk !! (`ocean_vcoord_z_fixed_target` with `z_top = 0`): !! !! 1. at `η = 0`, count the bed fillers `n_f` — the layers whose !! nominal top is at or below the bed less the partial-cell !! floor (`Z_FIXED_BED_PARTIAL_MIN`), exactly `z_fixed`'s rule; !! 2. lay the column again from `H + η`: layers `k <= n_f` at !! `h_min`, every live interface at the DILATED nominal depth !! `s·Z_k` below the free surface, with !! `s = (H + η − n_f·h_min)/(H − n_f·h_min)`. !! !! `s` is MOM6's `stretching` with the filler stack taken out of the !! dilation (MOM6 dilates all of `H` and then clamps the fillers back !! to `min_thickness`, which leaves `(s−1)·n_f·h_min` in the partial !! cell — under 1 mm per metre of `η` on the 1-degree Southern Ocean, !! `python_prototypes/mom6_zstar/mom6_zstar.py`). Consequences: !! !! * the liveness of every layer is decided at `η = 0`, so the !! pattern is static BY CONSTRUCTION and the static face mask !! (`zfixed_closed_faces`) applies with every consumer unchanged; !! * every live layer, the partial bed cell included, is its !! `η = 0` thickness times `s` (to round-off) — the dilation keeps !! ratios, so a partial cell (`p0 > Z_FIXED_BED_PARTIAL_MIN` at !! `η = 0` by construction) stays above `H_VANISHED` for every !! `s > H_VANISHED/p0`, i.e. for every `η` that does not all but !! dry the column; !! * at `η = 0`, `s == 1` exactly (numerator and denominator are the !! same expression), `s·Z_k == Z_k`, and the walk is `z_fixed`'s !! operation for operation: the `η = 0` z* target IS the `η = 0` !! `z_fixed` target, bit for bit. !! !! `Σ_k target_h = H + η` by construction (each step assigns what it !! takes from `z_below`, and the surface layer closes the column). !! A degenerate column (`H − n_f·h_min <= h_min`: land, or thinner !! than its filler stack) is laid with `s = 1`, `z_fixed`'s own !! degenerate overshoot. `s` is floored at 0: a column drained past !! its fillers is a dry column, which this coordinate does not !! support (`validate_config` refuses it with wet/dry). integer, intent(in) :: nx !! i-extent of every array (total, incl. halos). integer, intent(in) :: ny !! j-extent of every array (total, incl. halos). integer, intent(in) :: nz !! Number of layers; `k = 1` is the bed, `k = nz` the surface. real(wp), intent(out) :: target_h(nx, ny, nz) !! Target layer thickness (m), bottom-up. real(wp), intent(in) :: total_h(nx, ny) !! Column reference thickness `H` (m) — `Σ h_layer − η`. real(wp), intent(in) :: eta(nx, ny) !! Free-surface anomaly `η` (m). real(wp), intent(in) :: h_nominal !! Uniform nominal spacing `z_fixed_h_ref/nz` (m). Unused when !! `use_profile`. logical, intent(in) :: use_profile !! Take the nominal interfaces from `zi`. real(wp), intent(in) :: zi(0:nz) !! Nominal interface depths (m), bottom-up, `zi(k)` = top of layer !! `k`, `zi(nz) = 0`. Read only when `use_profile`. real(wp), intent(in) :: h_min !! Inert-filler thickness (`zstar_h_min`, `<= H_VANISHED`). integer :: i, j, k, n_f real(wp) :: z_below_loc, z_nom_loc, z_above_loc, stretch_loc, l0_loc real(wp) :: bed_partial_min logical :: walking_loc bed_partial_min = max(h_min, Z_FIXED_BED_PARTIAL_MIN) do concurrent(j=1:ny, i=1:nx) & local(k, n_f, z_below_loc, z_nom_loc, z_above_loc, stretch_loc, l0_loc, & walking_loc) ! Pass 1 — the eta = 0 bed fillers (z_fixed's bed rule, z_top = 0). n_f = 0 z_below_loc = total_h(i, j) walking_loc = .true. do k = 1, nz - 1 if (walking_loc) then if (use_profile) then z_nom_loc = zi(k) else z_nom_loc = real(nz - k, wp)*h_nominal end if if (z_nom_loc >= z_below_loc - bed_partial_min) then n_f = n_f + 1 z_below_loc = z_below_loc - h_min else walking_loc = .false. end if end if end do ! The dilation of the live column. l0_loc = total_h(i, j) - real(n_f, wp)*h_min if (l0_loc > h_min) then stretch_loc = max(((total_h(i, j) + eta(i, j)) - real(n_f, wp)*h_min)/l0_loc, & 0.0_wp) else stretch_loc = 1.0_wp end if ! Pass 2 — the z_fixed walk on the dilated nominal interfaces. z_below_loc = total_h(i, j) + eta(i, j) do k = 1, nz if (k <= n_f) then target_h(i, j, k) = h_min z_below_loc = z_below_loc - h_min else if (k == nz) then ! The surface layer closes the column (z_fixed's top rule). if (0.0_wp > z_below_loc - h_min) then target_h(i, j, k) = h_min else target_h(i, j, k) = z_below_loc end if else if (use_profile) then z_nom_loc = zi(k) else z_nom_loc = real(nz - k, wp)*h_nominal end if z_above_loc = stretch_loc*z_nom_loc target_h(i, j, k) = z_below_loc - z_above_loc z_below_loc = z_above_loc end if end do end do end subroutine ocean_vcoord_zstar_target pure subroutine ocean_vcoord_eta0_target(vc, target_h, total_h, nx, ny, nz) !! The target layer thickness at `eta = 0` of a GEOMETRIC family !! that vanishes layers — the ONE definition of "live" that the !! partial-step face mask (`configure_ocean_closed_faces`) and the !! on-target initial seed share with the running ALE regrid. !! !! * `VCOORD_Z_FIXED` — `ocean_vcoord_z_fixed_target` with the !! slot's nominal layering (uniform `z_fixed_h_ref/nz` or the !! stretched `z_fixed_zi`/`z_fixed_dz` profile) and its rigid top !! `z_top`: exactly the call `configure_ocean_closed_faces` has !! always made, argument for argument. !! * `VCOORD_ZSTAR_FULL` — the `ZSTAR_FULL` branch of !! `ocean_vcoord_geometric_target`, the kernel !! `compute_target_h` dispatches to every regrid, walking the !! per-column `z_ref` table `build_zref_full` laid from the !! bathymetry. The caller must have built that table first. !! * `VCOORD_ZSTAR` — `ocean_vcoord_zstar_target`, the kernel the !! regrid dispatches to, at `eta = 0` (where it is the `z_fixed` !! target with `z_top = 0`, bit for bit). Its live/filler pattern !! is decided at `eta = 0` inside the kernel, so it is EXACTLY the !! pattern of every later regrid, of either sign of `eta`. !! !! Any other family leaves `target_h` untouched (the caller refuses !! it before getting here). Configure / seed time only: the host !! arrays are local, so the `do concurrent` kernels underneath get !! their own implicit data regions. integer, intent(in) :: nx !! i-extent (total, incl. halos). integer, intent(in) :: ny !! j-extent (total, incl. halos). integer, intent(in) :: nz !! Number of layers; `k = 1` is the bed. type(ocean_vcoord_t), intent(in) :: vc !! The vertical-coordinate slot (coord_type + every table). real(wp), intent(inout) :: target_h(nx, ny, nz) !! Target thickness at `eta = 0` (m), bottom-up. real(wp), intent(in) :: total_h(nx, ny) !! Reference column thickness (m) — `bt_H_ref` at configure, !! the seed's water column at IC time. real(wp), allocatable :: eta0(:, :) real(wp) :: h_nominal allocate (eta0(nx, ny), source=0.0_wp) select case (vc%coord_type) case (VCOORD_Z_FIXED) h_nominal = vc%z_fixed_h_ref/real(nz, wp) call ocean_vcoord_z_fixed_target(target_h, total_h, eta0, vc%z_top, & nx, ny, nz, h_nominal, & vc%z_fixed_use_profile, vc%z_fixed_zi, & vc%z_fixed_dz, vc%zstar_h_min) case (VCOORD_ZSTAR_FULL) call ocean_vcoord_geometric_target(VCOORD_ZSTAR_FULL, nx, ny, nz, target_h, & total_h, eta0, vc%dsig, vc%z_ref_global, & vc%z_ref, vc%zsigma_depth_transition, & vc%zsigma_blend_width, vc%zstar_h_min) case (VCOORD_ZSTAR) h_nominal = vc%z_fixed_h_ref/real(nz, wp) call ocean_vcoord_zstar_target(target_h, total_h, eta0, nx, ny, nz, h_nominal, & vc%z_fixed_use_profile, vc%z_fixed_zi, & vc%zstar_h_min) case default end select deallocate (eta0) end subroutine ocean_vcoord_eta0_target pure subroutine ocean_vcoord_closed_face_masks(open_u, open_v, target_h, & nx, ny, nz, h_vanished) !! Partial-step z-level FACE CLOSURE mask for `VCOORD_Z_FIXED` !! (`&vcoord_nml zfixed_closed_faces`; Adcroft, Hill & Marshall !! 1997; Losch 2008 §2.1 for the ice-shelf cavity). !! !! Under `z_fixed` a layer whose nominal geopotential range lies !! inside the bed — or inside the ice draft — is an inert FILLER of !! thickness `zstar_h_min` (`<= h_vanished`). A velocity face at !! which layer `k` is a filler on EITHER side is not a thin !! passage: geometrically there is no water there on one side, so !! it is a **WALL for that layer** — no normal velocity, no mass or !! tracer flux, free-slip on the tangential component. This marks !! those faces. !! !! ### The rule, in one line !! !! ``` !! open_u(I,j,k) = 1 iff target_h(I-1,j,k) > h_vanished !! .and. target_h(I ,j,k) > h_vanished !! ``` !! and the v-face mirror. The mask is **STATIC**: the bed and the !! draft are static, and under `z_fixed` `η` is absorbed by the !! first LIVE layer (the partial cell) — a filler's target is !! `zstar_h_min` whatever `η` does — so the live/filler pattern !! does not move. Build it once, from !! `ocean_vcoord_z_fixed_target` at `η = 0`, so there is exactly !! ONE definition of "live" shared with the ALE regrid and the IC !! seed. !! !! ### Composition !! !! The mask is a THIRD, independent factor on the face width, not a !! replacement for either of the other two: !! ``` !! dy_eff(I,j,k) = dy_cu(I,j) · por_face_area_u(I,j,k) · open_u(I,j,k) !! land (2-D) porous (subgrid) z-level (per layer) !! ``` !! !! ### Array-edge faces !! !! `I = 1` and `I = nx+1` are left fully OPEN (1), exactly as the !! porous kernel leaves them: their `dy_cu` is already zero and the !! continuity wall zeroing owns them. The PHYSICAL seam of a !! periodic axis is an interior index (`nghost+1`), so it is !! covered by the `2:nx` sweep — provided the caller built !! `target_h` from GHOST-FILLED inputs (`bt_H_ref` and `z_top` are !! periodic-wrapped + halo-exchanged before this runs). integer, intent(in) :: nx !! i-extent of the CENTRE arrays (total, incl. halos). integer, intent(in) :: ny !! j-extent of the CENTRE arrays (total, incl. halos). integer, intent(in) :: nz !! Number of layers; `k = 1` is the bed, `k = nz` the top. real(wp), intent(out) :: open_u(nx + 1, ny, nz) !! u-face 0/1 open mask. real(wp), intent(out) :: open_v(nx, ny + 1, nz) !! v-face 0/1 open mask. real(wp), intent(in) :: target_h(nx, ny, nz) !! The `z_fixed` target thickness at `η = 0`. real(wp), intent(in) :: h_vanished !! Inert-filler marker (`H_VANISHED`). A layer is LIVE iff its !! target thickness is strictly greater than this. integer :: i, j, k do concurrent(k=1:nz, j=1:ny, i=1:nx + 1) open_u(i, j, k) = 1.0_wp end do do concurrent(k=1:nz, j=1:ny + 1, i=1:nx) open_v(i, j, k) = 1.0_wp end do do concurrent(k=1:nz, j=1:ny, i=2:nx) if (target_h(i - 1, j, k) > h_vanished .and. & target_h(i, j, k) > h_vanished) then open_u(i, j, k) = 1.0_wp else open_u(i, j, k) = 0.0_wp end if end do do concurrent(k=1:nz, j=2:ny, i=1:nx) if (target_h(i, j - 1, k) > h_vanished .and. & target_h(i, j, k) > h_vanished) then open_v(i, j, k) = 1.0_wp else open_v(i, j, k) = 0.0_wp end if end do end subroutine ocean_vcoord_closed_face_masks pure subroutine ocean_vcoord_k_top_from_target(k_top, k_top_u, k_top_v, & target_h, nx, ny, nz, h_vanished) !! The shared FIRST-LIVE-LAYER index, counting down from the top — !! `multilayer_state_t%k_top` and its two face twins — built from !! a layer-thickness field. !! !! ### The rule !! !! ``` !! k_top(i,j) = the largest k with target_h(i,j,k) > h_vanished, !! or nz when the column has none !! ``` !! !! **Strict `>`**, matching the remap drain's `H_FLOOR` !! (`rdb_ocean_remap.F90`), the melt far-field sampler's !! `<= H_VANISHED ⇒ cycle`, and `ocean_vcoord_closed_face_masks` !! above: a layer sitting exactly ON the marker is dead on every !! side of the contract. `h_vanished` is `H_VANISHED`, never a !! slot-local `h_min = 1.0e-3` anti-zero floor — those are a !! different thing under the D4 taxonomy. !! !! **The `nz` fallback is what makes the whole indirection free.** !! On sigma, z*-lite, and every geometric family that does not !! vanish a layer against the top, no wet column has !! `h(:,:,nz) <= h_vanished`, so `k_top ≡ nz` and every consumer !! that reads `k_top(i,j)` instead of `nz` reads the same memory !! with the same arithmetic. A land column (every layer at the !! marker under the land-state contract) also lands on `nz`, which !! is what those consumers index today, and is then masked out by !! `wet_mask` exactly as before. !! !! ### The face rule is `min`, not `max` !! !! `k_top_u(I,j) = min(k_top(I-1,j), k_top(I,j))` and the v-face !! mirror. A velocity face carries water in layer `k` only where !! BOTH abutting columns are live there — which is precisely the !! statement `ocean_vcoord_closed_face_masks` makes with !! `open_u = 1 iff target_h > h_vanished on both sides` — so the !! shallowest layer the FACE has is the DEEPER of the two column !! tops, i.e. the SMALLER index. `max` would hand the ice-ocean !! top drag and the implicit stress/drag fold a row that is a !! filler on one side, which is the bug this field exists to stop. !! With `k_top ≡ nz` everywhere, `min(nz, nz) = nz` ⇒ the face !! twins are bit-identical too. !! !! Array-edge faces (`I = 1`, `I = nx+1`) take the one column they !! have, mirroring the mask builder leaving them fully open: their !! `dy_cu` is already zero and the wall zeroing owns them. integer, intent(in) :: nx !! i-extent of the CENTRE arrays (total, incl. halos). integer, intent(in) :: ny !! j-extent of the CENTRE arrays (total, incl. halos). integer, intent(in) :: nz !! Number of layers; `k = 1` is the bed, `k = nz` the top. integer, intent(out) :: k_top(nx, ny) !! Cell-centred first live layer. integer, intent(out) :: k_top_u(nx + 1, ny) !! u-face twin. integer, intent(out) :: k_top_v(nx, ny + 1) !! v-face twin. real(wp), intent(in) :: target_h(nx, ny, nz) !! Layer thickness (m) the live/filler pattern is read from — !! the `z_fixed` target at `eta = 0` at configure time. real(wp), intent(in) :: h_vanished !! Inert-filler marker (`H_VANISHED`). A layer is LIVE iff its !! thickness is strictly greater than this. integer :: i, j, k, ia, ib, ka, kb ! Each loop is SELF-CONTAINED — it reads `target_h` and writes ONE ! output array, and no loop reads what another wrote. That is ! deliberate: under `-stdpar=gpu` with `mem:separate` each ! `do concurrent` gets its OWN implicit data region for the host ! arrays it touches, so a face loop that read the `k_top` a ! previous loop had just written saw the array as write-only and ! came back with the wrong indices on device while being correct on ! the host. Re-scanning the column costs two extra passes over ! `target_h` ONCE, at configure. do concurrent(j=1:ny, i=1:nx) local(k) k_top(i, j) = nz do k = nz, 1, -1 if (target_h(i, j, k) > h_vanished) then k_top(i, j) = k exit end if end do end do ! `ia`/`ib` clamp to the one column an ARRAY-EDGE face has, so ! `min(ka, kb)` degenerates to that column's own index there — the ! mask builder leaves those faces fully open for the same reason. do concurrent(j=1:ny, i=1:nx + 1) local(k, ia, ib, ka, kb) ia = max(1, i - 1) ib = min(nx, i) ka = nz do k = nz, 1, -1 if (target_h(ia, j, k) > h_vanished) then ka = k exit end if end do kb = nz do k = nz, 1, -1 if (target_h(ib, j, k) > h_vanished) then kb = k exit end if end do k_top_u(i, j) = min(ka, kb) end do do concurrent(j=1:ny + 1, i=1:nx) local(k, ia, ib, ka, kb) ia = max(1, j - 1) ib = min(ny, j) ka = nz do k = nz, 1, -1 if (target_h(i, ia, k) > h_vanished) then ka = k exit end if end do kb = nz do k = nz, 1, -1 if (target_h(i, ib, k) > h_vanished) then kb = k exit end if end do k_top_v(i, j) = min(ka, kb) end do end subroutine ocean_vcoord_k_top_from_target pure subroutine ocean_vcoord_k_bot_from_target(k_bot, k_bot_u, k_bot_v, & target_h, nx, ny, nz, h_vanished) !! The shared FIRST-LIVE-LAYER index counting UP from the bed — !! `multilayer_state_t%k_bot` and its two face twins. The bed-side !! mirror of `ocean_vcoord_k_top_from_target`, built from the same !! layer-thickness field by the same strict `> h_vanished` test. !! !! ### The rule !! !! ``` !! k_bot(i,j) = the smallest k with target_h(i,j,k) > h_vanished, !! or 1 when the column has none !! ``` !! !! **The `1` fallback is what makes the indirection free**: on every !! family without a static bed filler `h(:,:,1) > h_vanished` on a !! wet column, so `k_bot ≡ 1` and a consumer reading `k_bot(i,j)` !! instead of `1` reads the same memory. A land column (every layer !! AT the marker) also lands on `1`, and `wet_mask` zeroes it as !! before. !! !! ### The face rule is `max`, not `min` !! !! `k_bot_u(I,j) = max(k_bot(I-1,j), k_bot(I,j))`. A face carries !! water in layer `k` only where BOTH abutting columns are live !! there (`ocean_vcoord_closed_face_masks`), so the deepest layer !! the FACE has is the SHALLOWER of the two column bottoms — the !! LARGER index. `min` would hand the bottom drag and the vdiff bed !! row a layer that is a filler on one side, i.e. a closed face. !! With `k_bot ≡ 1`, `max(1, 1) = 1` ⇒ bit-identical. !! !! Array-edge faces (`I = 1`, `I = nx+1`) take the one column they !! have, as `k_top`'s do; the configure driver face-halo-exchanges !! the twins afterwards so a tile seam / periodic wrap carries the !! owner's value there. integer, intent(in) :: nx !! i-extent of the CENTRE arrays (total, incl. halos). integer, intent(in) :: ny !! j-extent of the CENTRE arrays (total, incl. halos). integer, intent(in) :: nz !! Number of layers; `k = 1` is the bed, `k = nz` the top. integer, intent(out) :: k_bot(nx, ny) !! Cell-centred first live layer counting up from the bed. integer, intent(out) :: k_bot_u(nx + 1, ny) !! u-face twin. integer, intent(out) :: k_bot_v(nx, ny + 1) !! v-face twin. real(wp), intent(in) :: target_h(nx, ny, nz) !! Layer thickness (m) the live/filler pattern is read from — !! the `z_fixed` target at `eta = 0` at configure time. real(wp), intent(in) :: h_vanished !! Inert-filler marker (`H_VANISHED`). LIVE iff strictly greater. integer :: i, j, k, ia, ib, ka, kb ! Each loop is SELF-CONTAINED (reads `target_h`, writes ONE output) ! for the reason `ocean_vcoord_k_top_from_target` documents: under ! `-stdpar=gpu` + `mem:separate` a face loop that read the centre ! index a previous loop had just written saw it as write-only. do concurrent(j=1:ny, i=1:nx) local(k) k_bot(i, j) = 1 do k = 1, nz if (target_h(i, j, k) > h_vanished) then k_bot(i, j) = k exit end if end do end do do concurrent(j=1:ny, i=1:nx + 1) local(k, ia, ib, ka, kb) ia = max(1, i - 1) ib = min(nx, i) ka = 1 do k = 1, nz if (target_h(ia, j, k) > h_vanished) then ka = k exit end if end do kb = 1 do k = 1, nz if (target_h(ib, j, k) > h_vanished) then kb = k exit end if end do k_bot_u(i, j) = max(ka, kb) end do do concurrent(j=1:ny + 1, i=1:nx) local(k, ia, ib, ka, kb) ia = max(1, j - 1) ib = min(ny, j) ka = 1 do k = 1, nz if (target_h(i, ia, k) > h_vanished) then ka = k exit end if end do kb = 1 do k = 1, nz if (target_h(i, ib, k) > h_vanished) then kb = k exit end if end do k_bot_v(i, j) = max(ka, kb) end do end subroutine ocean_vcoord_k_bot_from_target pure function ocean_vcoord_count_ledges(open_u, open_v, target_h, & nx, ny, nz, h_vanished) result(n_ledge) !! Count LEDGE cells: a cell that is LIVE at layer `k` but all four !! of whose own-layer faces are closed, i.e. water the mask has !! isolated. A ledge needs a one-cell-wide spike in the bed or the !! draft; it is inert by construction (no flux in or out, and its !! velocity is zeroed every stage), but a non-zero count is worth !! saying out loud once at configure, because it means the mask is !! walling off real water. !! !! Interior cells only (`2:nx-1`, `2:ny-1`) — the ghost ring has no !! four-face neighbourhood of its own. integer, intent(in) :: nx, ny, nz real(wp), intent(in) :: open_u(nx + 1, ny, nz) real(wp), intent(in) :: open_v(nx, ny + 1, nz) real(wp), intent(in) :: target_h(nx, ny, nz) real(wp), intent(in) :: h_vanished integer :: n_ledge integer :: i, j, k n_ledge = 0 do k = 1, nz do j = 2, ny - 1 do i = 2, nx - 1 if (target_h(i, j, k) <= h_vanished) cycle if (open_u(i, j, k) == 0.0_wp .and. & open_u(i + 1, j, k) == 0.0_wp .and. & open_v(i, j, k) == 0.0_wp .and. & open_v(i, j + 1, k) == 0.0_wp) then n_ledge = n_ledge + 1 end if end do end do end do end function ocean_vcoord_count_ledges pure function ocean_vcoord_count_bed_steps(target_h, total_h, dy_cu, dx_cv, & nx, ny, nz, i0, i1, j0, j1, & h_vanished) result(n_step) !! Count the wet velocity faces at which the two columns' BED falls !! in different nominal `z_fixed` layers — a bed staircase step that !! crosses a nominal interface, so that an OPEN face pairs a live !! layer on one side with a bed FILLER on the other. That is the !! face `&vcoord_nml zfixed_closed_faces` closes; left open, the FV !! pressure gradient across the step drives the flow from rest !! (`configure_ocean_closed_faces` refuses the configuration). !! !! The bed layer of a column is its LOWEST live layer under the !! `z_fixed` target at `η = 0` (`target_h > h_vanished`) — the same !! single definition of "live" the closed-face mask is built from, !! so the count is exactly the set of faces whose bed-side layers !! the mask would close. A column with no live layer, or with !! `total_h <= 0`, is dry and pairs with nothing; a face whose !! land-masked width is zero is a wall and is skipped. Top-side !! (ice-draft) fillers are deliberately not counted: only the bed !! side is refused. !! !! Only the OWNED faces are visited — `i0:i1` / `j0:j1` are the !! owned CELL ranges and each owned cell contributes its WEST !! (`I = i`) and SOUTH (`J = j`) face, which pairs it with the !! ghost-filled neighbour — so a step at a tile seam or a periodic !! seam is counted exactly once across the decomposition and the !! per-rank counts sum to the global one. integer, intent(in) :: nx !! i-extent of the CENTRE arrays (total, incl. halos). integer, intent(in) :: ny !! j-extent of the CENTRE arrays (total, incl. halos). integer, intent(in) :: nz !! Number of layers; `k = 1` is the bed, `k = nz` the top. real(wp), intent(in) :: target_h(nx, ny, nz) !! The `z_fixed` target thickness at `η = 0`. real(wp), intent(in) :: total_h(nx, ny) !! Column reference thickness (`bt_H_ref`, ghost-filled). real(wp), intent(in) :: dy_cu(nx + 1, ny) !! Land-masked u-face width. real(wp), intent(in) :: dx_cv(nx, ny + 1) !! Land-masked v-face width. integer, intent(in) :: i0, i1, j0, j1 !! Owned cell range (`nghost+1 : nghost+n_phys`). real(wp), intent(in) :: h_vanished !! Inert-filler marker (`H_VANISHED`). integer :: n_step integer :: i, j n_step = 0 do j = j0, j1 do i = i0, i1 if (dy_cu(i, j) > 0.0_wp) then if (bed_layer(i - 1, j) /= bed_layer(i, j) .and. & bed_layer(i - 1, j) > 0 .and. bed_layer(i, j) > 0) then n_step = n_step + 1 end if end if if (dx_cv(i, j) > 0.0_wp) then if (bed_layer(i, j - 1) /= bed_layer(i, j) .and. & bed_layer(i, j - 1) > 0 .and. bed_layer(i, j) > 0) then n_step = n_step + 1 end if end if end do end do contains pure function bed_layer(ii, jj) result(kb) !! Lowest live layer of column `(ii, jj)`; `0` if dry. integer, intent(in) :: ii, jj integer :: kb integer :: kk kb = 0 if (total_h(ii, jj) <= 0.0_wp) return do kk = 1, nz if (target_h(ii, jj, kk) > h_vanished) then kb = kk return end if end do end function bed_layer end function ocean_vcoord_count_bed_steps pure subroutine ocean_vcoord_compute_target_h_rho(this, total_h, eta, T, S, eos, hybrid) !! Thin polymorphic wrapper for the isopycnal (`VCOORD_RHO`) and !! hybrid z*/isopycnal (`VCOORD_HYCOM`) target-grid build. Mirrors !! `ocean_vcoord_compute_target_h` (resolves the concrete type so !! the `do concurrent` kernel runs on `type(ocean_vcoord_t)`, never !! a polymorphic list item). Separate from `compute_target_h` !! because the RHO/HYCOM branch needs the per-layer T/S !! concentrations and the device-resident EOS coefficients, which !! the shared `pure (total_h, eta)` TBP cannot carry. !! !! `hybrid` selects the variant: `.false.` (default) = pure !! `VCOORD_RHO` (bit-identical with the P2 kernel); `.true.` = !! `VCOORD_HYCOM` (adds the bottom-up density monotonize and the !! z* nominal-floor sweep inside the same column kernel). class(ocean_vcoord_t), intent(inout) :: this real(wp), intent(in) :: total_h(:, :) ! assumed-shape-ok: thin TBP wrapper, no kernel real(wp), intent(in) :: eta(:, :) ! assumed-shape-ok: thin TBP wrapper, no kernel real(wp), intent(in) :: T(:, :, :) ! assumed-shape-ok: thin TBP wrapper, no kernel real(wp), intent(in) :: S(:, :, :) ! assumed-shape-ok: thin TBP wrapper, no kernel type(eos_t), intent(in) :: eos logical, intent(in), optional :: hybrid !! Enable the HYCOM hybrid deltas (default `.false.` = pure RHO). logical :: hybrid_loc hybrid_loc = .false. if (present(hybrid)) hybrid_loc = hybrid select type (this) type is (ocean_vcoord_t) call ocean_vcoord_compute_target_h_rho_impl(this, total_h, eta, T, S, eos, hybrid_loc) end select end subroutine ocean_vcoord_compute_target_h_rho pure subroutine ocean_vcoord_compute_target_h_rho_impl(this, total_h, eta, T, S, eos, hybrid) !! Isopycnal regrid: place layer interfaces on the prescribed !! `rho_target(0:nz)` potential-density surfaces. ONE !! `do concurrent(j,i)` over columns, each running the full !! density-space inversion on `NZ_STACK_MAX` fixed-size locals — !! no host loop, no per-call allocate. !! !! Algorithm (clean-room from Bleck 2002 / White & Adcroft 2008; !! MOM6 `coord_rho` is the behavioural oracle), per the spec: !! !! 0. Pre-compaction — strip source layers `h <= H_MIN`, donate !! their volume to the thickest survivor (volume-conserving), !! build the survivor count. Fast path: `<= 1` survivor → !! `h_new = h_old` (no inversion). !! 1. Layer potential densities via `eos_density_point` at !! `rho_ref_pressure` on the compacted column. !! 2. PPM (Colella-Woodward monotone) reconstruction of the !! density profile over the compacted thicknesses. !! 3. Per interior target, bracket-ordered inversion: light !! boundary → surface; discontinuous-jump sweep; dense !! boundary → bed; else fixed-8-iter Newton on `xi in [0,1]` !! (convergence on `|delta| < NR_TOL` AFTER `xi += delta`; !! zero-gradient `NR_OFFSET` escape at both ends; masked !! fallback to the previous interface on no-bracket). !! 4. Monotone non-decreasing interfaces → `h_new`. !! 5. MOM6 min-thickness inflation, floor = `max(zstar_h_min, !! H_VANISHED)` (MUST-HAVE #2: keeps RHO-collapsed layers !! above the remap-drain `H_FLOOR` so the next regrid does !! not zero their tracer mass); debit the single thickest !! layer once — or, when that would take it below the floor, !! every above-floor layer in proportion to its excess. !! (guard) A column thinner than `nz·h_floor_eff` cannot hold !! every layer at the floor (land columns hold `nz·H_VANISHED`): !! it keeps `h_new = h_old` (the remap is the identity on it). !! Without the guard step 5 wrote a negative thickness or !! minted mass on every such column (audit finding H3). !! !! Internal working frame is TOP-DOWN (index 1 = surface, +down), !! matching the validated prototype and the `rho_target(0)` = !! lightest = surface convention. The final assignment FLIPS to !! the bottom-up state (MUST-HAVE #3): `target_h(:,:,k) = !! h_new_td(nz - k + 1)`, so the lightest target lands at the !! surface (k=nz) and the densest at the bed (k=1). Sum is !! conserved exactly so `sum_k target_h = H` (η is implicit in !! `total_h` here — the caller passes the live column total as !! `total_h`, see `ocean_apply_ale_remap_step`). !! !! HYCOM hybrid (`hybrid = .true.`, Bleck 2002 / MOM6 !! `build_hycom1_column`): two deltas around the unchanged RHO !! inversion, both inside this same column kernel. !! (1b) BOTTOM-UP density monotonize before the PPM reconstruction: !! in the top-down work frame, cap each cell by the one below !! it (toward the bed) — `do k=nk-1,1,-1: rhoc(k)=min(rhoc(k), !! rhoc(k+1))`. Pure RHO omits this (assumes a monotone !! profile + leans on the PPM limiter); HYCOM enforces it so !! the inversion always sees a monotone column. !! (4b) z* NOMINAL-FLOOR sweep after the inversion, before the !! monotone/inflation: walk interior+bottom interfaces from !! the surface down, accumulating the z* NOMINAL THICKNESS IN !! METRES times `stretching = (H+η)/H`, and push each interface !! DOWN to at least that depth (clamped to the column bottom). !! The nominal thicknesses are the z* coordinate resolution — !! the same `&vcoord_nml z_fixed_profile` table `z_fixed` uses !! (`z_fixed_dz` for "list"/"tanh", `max_depth/nz` for !! "uniform"), as MOM6 HYCOM1 takes its `coordinateResolution` !! from the ALE_COORDINATE_CONFIG that would define a z* grid. !! So the band is a fixed depth range in every column (a 2 m !! surface layer stays 2 m over the shelf and the abyss alike) !! and a shallow column's deeper interfaces clamp onto its bed; !! deep isopycnal interfaces already below the floor are !! untouched. Until 2026-10-02 the increment was the column !! FRACTION `dsig·(H+η)` with `dsig ≡ 1/nz` — a sigma floor that !! set 95 % of the 1-degree Southern Ocean's interfaces (audit !! finding H1); it survives only as the fallback for a slot no !! setup path configured (`z_fixed_h_ref <= 0`). !! CRITICAL: `stretching` multiplies a nominal thickness that is !! referenced to `H` (= `total_h`), not `H+η`; scaling by !! `(H+η)` twice over-stretches by `(H+η)/H` — invisible at !! η=0, wrong with a free surface. No renormalize after the !! floor sweep (MOM6 pins the bottom interface + uses the !! debit-thickest inflation instead). !! `hybrid = .false.` (the `VCOORD_RHO` path) skips BOTH deltas and !! is bit-identical to the P2 kernel. !! !! Host dispatcher: resolves the scalar knobs and hands every !! component to the flat `ocean_vcoord_rho_target` kernel as an !! explicit-shape / scalar dummy, so no `this%` reference and no !! `associate`-name reaches the `do concurrent` (see that kernel's !! docstring for why that is load-bearing on the GPU build). type(ocean_vcoord_t), intent(inout) :: this ! assumed-shape-ok: cadence-bounded (once per outer ALE step); forwarded ! to the explicit-shape kernel below, which is where the loop runs. real(wp), intent(in) :: total_h(:, :) !! Column reference depth H(i, j) (m). Caller passes the live !! column total (sum of h_layer) so the new grid spans it exactly. real(wp), intent(in) :: eta(:, :) ! assumed-shape-ok: see total_h !! Free-surface anomaly η(i, j) (m). Added to `total_h` to form !! the column extent the new interfaces span. real(wp), intent(in) :: T(:, :, :) ! assumed-shape-ok: see total_h !! Layer-mean potential temperature concentration (°C), `(nx,ny,nz)`. real(wp), intent(in) :: S(:, :, :) ! assumed-shape-ok: see total_h !! Layer-mean salinity concentration (PSU), `(nx,ny,nz)`. type(eos_t), intent(in) :: eos !! Shared device-resident EOS handle (flat POD, by value). logical, intent(in) :: hybrid !! `.true.` = HYCOM (apply the monotonize + z*-floor deltas); !! `.false.` = pure RHO (bit-identical with the P2 kernel). real(wp) :: h_floor_eff, h_nominal integer :: floor_mode if (.not. this%is_init) return ! Inflation floor must be STRICTLY above H_VANISHED: the remap drain ! (`ocean_remap_tracer_field`) gates on `h_old > H_FLOOR` (== H_VANISHED) ! with a strict `>`, so a layer sitting exactly at H_VANISHED has its ! tracer concentration zeroed when this column is fed back as `h_old` ! on the next regrid. Floor at 2·H_VANISHED so inflated layers always ! survive the drain (closes the multi-regrid mass-loss footgun). h_floor_eff = max(this%zstar_h_min, 2.0_wp*H_VANISHED) ! HYCOM z* nominal floor: METRES from the z* coordinate resolution ! (the `z_fixed` nominal profile — MOM6 HYCOM1 reads its ! `coordinateResolution` from the same ALE_COORDINATE_CONFIG that ! sets a z* grid), stretched by (H+η)/H. A stretched profile ! (`z_fixed_use_profile`) gives per-layer `z_fixed_dz`; otherwise the ! uniform `z_fixed_h_ref/nz` (setup always writes `max_depth` there). ! Only a slot nobody configured (`z_fixed_h_ref <= 0`: unit tests that ! build the vcoord by hand) falls back to the historical column- ! fraction floor `dsig·(H+η)`. if (this%z_fixed_use_profile) then floor_mode = HYCOM_FLOOR_PROFILE h_nominal = 0.0_wp else if (this%z_fixed_h_ref > 0.0_wp) then floor_mode = HYCOM_FLOOR_UNIFORM h_nominal = this%z_fixed_h_ref/real(this%nz_ml, wp) else floor_mode = HYCOM_FLOOR_SIGMA h_nominal = 0.0_wp end if call ocean_vcoord_rho_target(this%nx_total, this%ny_total, this%nz_ml, & this%target_h, this%remap_h_old, total_h, eta, & T, S, this%dsig, this%z_fixed_dz, this%rho_target, eos, & this%rho_ref_pressure, this%zstar_h_min, & h_floor_eff, hybrid, floor_mode, h_nominal) end subroutine ocean_vcoord_compute_target_h_rho_impl pure subroutine ocean_vcoord_rho_target(nx, ny, nz, target_h, remap_h_old, & total_h, eta, t_conc, s_conc, dsig, & floor_dz, rho_target, eos, p_ref, h_min, & h_floor_eff, hybrid, floor_mode, h_nominal) !! Column kernel of the RHO / HYCOM regrid (algorithm: see !! `ocean_vcoord_compute_target_h_rho_impl`). Flat on purpose: !! every array is an explicit-shape dummy and every knob a scalar !! dummy — no derived-type component and no `associate` reaches the !! `do concurrent`. !! !! This is load-bearing on nvfortran `-stdpar=gpu -gpu=mem:separate`. !! The kernel used to run inside `associate (rho_ref_pressure => !! this%rho_ref_pressure, ...)` and hand that name BY REFERENCE to the !! out-of-module `!$acc routine seq` `eos_density_point`; the device !! callee received the HOST address of the component and faulted !! (`CUDA_ERROR_ILLEGAL_ADDRESS`, compute-sanitizer: invalid !! `__global__` read at `rdb_eos.F90` `p_plus_p0 = p + p_0`) on the !! first regrid of every RHO/HYCOM × Wright (or Roquet) run. The !! linear branch never reads `p`, which is why it hid. Gated by !! `tests/test_ocean_vcoord_wright_device.F90`. It also retires the !! `associate`-over-`do concurrent` shape CLAUDE.md forbids for ifx. integer, intent(in), value :: nx !! i-extent of every horizontal array (total, incl. halos). integer, intent(in), value :: ny !! j-extent of every horizontal array (total, incl. halos). integer, intent(in), value :: nz !! Number of layers; `k = 1` is the bed, `k = nz` the surface. real(wp), intent(inout) :: target_h(nx, ny, nz) !! Target layer thickness (m), bottom-up. real(wp), intent(in) :: remap_h_old(nx, ny, nz) !! Pre-remap layer thickness snapshot (m), bottom-up. real(wp), intent(in) :: total_h(nx, ny) !! Column reference depth H (m). real(wp), intent(in) :: eta(nx, ny) !! Free-surface anomaly η (m). real(wp), intent(in) :: t_conc(nx, ny, nz) !! Layer-mean potential temperature (°C), bottom-up. real(wp), intent(in) :: s_conc(nx, ny, nz) !! Layer-mean salinity (PSU), bottom-up. real(wp), intent(in) :: dsig(nz) !! Nominal layer fractions, bottom-up (`dsig(nz)` = surface) — !! the HYCOM floor only under `HYCOM_FLOOR_SIGMA`. real(wp), intent(in) :: floor_dz(nz) !! HYCOM z* nominal layer thicknesses (m), bottom-up !! (`floor_dz(nz)` = surface) — read only under `HYCOM_FLOOR_PROFILE`. real(wp), intent(in) :: rho_target(0:nz) !! Target potential densities (kg/m³), `0` = lightest = surface. type(eos_t), intent(in) :: eos !! Shared EOS handle (flat POD). real(wp), intent(in), value :: p_ref !! Coordinate reference pressure (Pa) — `rho_ref_pressure`. real(wp), intent(in), value :: h_min !! Pre-compaction strip threshold (m) — `zstar_h_min`. real(wp), intent(in), value :: h_floor_eff !! Min-thickness inflation floor (m), `> H_VANISHED`. logical, intent(in), value :: hybrid !! `.true.` = HYCOM deltas; `.false.` = pure RHO. integer, intent(in), value :: floor_mode !! HYCOM z* floor source: `HYCOM_FLOOR_PROFILE` (`floor_dz`), !! `HYCOM_FLOOR_UNIFORM` (`h_nominal`) or `HYCOM_FLOOR_SIGMA` !! (`dsig·H`, the unconfigured-slot fallback). real(wp), intent(in), value :: h_nominal !! Uniform z* nominal thickness (m) for `HYCOM_FLOOR_UNIFORM`. integer :: i, j ! One column per thread: the whole column walk is the same-module ! `!$acc routine seq` `ocean_vcoord_rho_target_column`, so its ! NZ_STACK_MAX work arrays are thread-private and no inner loop is ! spread across threads. (Inlined, nvfortran 26.5 auto-collapsed the ! (j,i) nest, promoted the work arrays AND `nk` to shared memory and ! vectorised the inner loops with shared-memory reductions whose ! scratch fell outside the kernel's shared allocation — ! compute-sanitizer: invalid `__shared__` write in the `donate` ! reduction, every RHO/HYCOM run, both EOS.) do concurrent(j=1:ny, i=1:nx) call ocean_vcoord_rho_target_column(i, j, nx, ny, nz, target_h, remap_h_old, & total_h, eta, t_conc, s_conc, dsig, & floor_dz, rho_target, eos, p_ref, h_min, & h_floor_eff, hybrid, floor_mode, h_nominal) end do end subroutine ocean_vcoord_rho_target pure subroutine ocean_vcoord_rho_target_column(i, j, nx, ny, nz, target_h, remap_h_old, & total_h, eta, t_conc, s_conc, dsig, & floor_dz, rho_target, eos, p_ref, h_min, & h_floor_eff, hybrid, floor_mode, h_nominal) !! One column of the RHO / HYCOM regrid (steps 0-5 of !! `ocean_vcoord_compute_target_h_rho_impl`) — the per-thread body of !! `ocean_vcoord_rho_target`. Same module as its caller (the !! project rule for `!$acc routine seq` callees of a `do concurrent`). !$acc routine seq integer, intent(in), value :: i !! Column i-index. integer, intent(in), value :: j !! Column j-index. integer, intent(in), value :: nx !! i-extent of every horizontal array (total, incl. halos). integer, intent(in), value :: ny !! j-extent of every horizontal array (total, incl. halos). integer, intent(in), value :: nz !! Number of layers; `k = 1` is the bed, `k = nz` the surface. real(wp), intent(inout) :: target_h(nx, ny, nz) !! Target layer thickness (m), bottom-up. real(wp), intent(in) :: remap_h_old(nx, ny, nz) !! Pre-remap layer thickness snapshot (m), bottom-up. real(wp), intent(in) :: total_h(nx, ny) !! Column reference depth H (m). real(wp), intent(in) :: eta(nx, ny) !! Free-surface anomaly η (m). real(wp), intent(in) :: t_conc(nx, ny, nz) !! Layer-mean potential temperature (°C), bottom-up. real(wp), intent(in) :: s_conc(nx, ny, nz) !! Layer-mean salinity (PSU), bottom-up. real(wp), intent(in) :: dsig(nz) !! Nominal layer fractions, bottom-up (`dsig(nz)` = surface) — !! the HYCOM floor only under `HYCOM_FLOOR_SIGMA`. real(wp), intent(in) :: floor_dz(nz) !! HYCOM z* nominal layer thicknesses (m), bottom-up !! (`floor_dz(nz)` = surface) — read only under `HYCOM_FLOOR_PROFILE`. real(wp), intent(in) :: rho_target(0:nz) !! Target potential densities (kg/m³), `0` = lightest = surface. type(eos_t), intent(in) :: eos !! Shared EOS handle (flat POD). real(wp), intent(in), value :: p_ref !! Coordinate reference pressure (Pa) — `rho_ref_pressure`. real(wp), intent(in), value :: h_min !! Pre-compaction strip threshold (m) — `zstar_h_min`. real(wp), intent(in), value :: h_floor_eff !! Min-thickness inflation floor (m), `> H_VANISHED`. logical, intent(in), value :: hybrid !! `.true.` = HYCOM deltas; `.false.` = pure RHO. integer, intent(in), value :: floor_mode !! HYCOM z* floor source: `HYCOM_FLOOR_PROFILE` (`floor_dz`), !! `HYCOM_FLOOR_UNIFORM` (`h_nominal`) or `HYCOM_FLOOR_SIGMA` !! (`dsig·H`, the unconfigured-slot fallback). real(wp), intent(in), value :: h_nominal !! Uniform z* nominal thickness (m) for `HYCOM_FLOOR_UNIFORM`. integer :: k, kk, nk, ns, idx_thick, src, ii integer :: mapping(NZ_STACK_MAX) real(wp) :: h_col(NZ_STACK_MAX), t_col(NZ_STACK_MAX), s_col(NZ_STACK_MAX) real(wp) :: hc(NZ_STACK_MAX), rhoc(NZ_STACK_MAX), rtgt(NZ_STACK_MAX) real(wp) :: z_new(NZ_STACK_MAX + 1), h_new(NZ_STACK_MAX) real(wp) :: col_extent, donate real(wp) :: total_need, thick_max, excess, frac real(wp) :: nominal_z, stretching, h_ref_col ! One column per (j,i). NZ_STACK_MAX fixed-size locals; no name ! shadows a Fortran intrinsic; cross-module pure EOS helper carries ! its own `!$acc routine seq`. ! --- gather TOP-DOWN: working index 1 = surface = state k=nz --- ! Source thicknesses come from the remap snapshot the ! orchestrator placed in `remap_h_old` (the live, pre-remap ! `h_layer`); T/S are the layer-mean concentrations. Flip the ! bottom-up state index (k=nz surface) into the top-down work ! frame (work index 1 = surface). col_extent = max(total_h(i, j) + eta(i, j), 0.0_wp) do k = 1, nz ii = nz - k + 1 ! state (bottom-up) index h_col(k) = remap_h_old(i, j, ii) t_col(k) = t_conc(i, j, ii) s_col(k) = s_conc(i, j, ii) end do ! --- too-thin column: cannot carry the coordinate, keep h_old --- ! A column thinner than `nz*h_floor_eff` (land / dry columns hold ! `nz*H_VANISHED`; a wet/dry or ice-cavity sliver can be any size) ! cannot have every layer at the inflation floor. Step 5 then either ! minted mass (`ns == 0`: every layer to the floor) or debited the one ! surviving layer below zero (`ns > 0`: a 50 x 1.5e-4 m land column ! collapsed into one layer gave 0.0075 - 49*3e-4 = -0.0072 m). Leave ! such a column exactly where it is — the same no-motion answer as the ! `nk <= 1` fast path below — so the remap is the identity on it: no ! negative thickness, no created mass, `sum(h)` preserved bit-for-bit. ! Every column the guard catches was mis-handled by step 5, so it is ! inert on every column that step handled correctly. if (col_extent < real(nz, wp)*h_floor_eff) then do k = 1, nz target_h(i, j, k) = remap_h_old(i, j, k) end do return end if ! --- step 0: pre-compaction (strip h <= h_min, donate) --- nk = 0 do k = 1, nz if (h_col(k) > h_min) then nk = nk + 1 mapping(nk) = k hc(nk) = h_col(k) end if end do if (nk <= 1) then ! Fast path: <= 1 finite layer. nz == nk_state here, so ! keep the source thicknesses unchanged (h_new = h_old), ! flipped back into the bottom-up state. do k = 1, nz target_h(i, j, k) = remap_h_old(i, j, k) end do return end if ! Donate the stripped volume to the thickest survivor. donate = col_extent do kk = 1, nk donate = donate - hc(kk) end do if (donate > 0.0_wp) then idx_thick = 1 do kk = 2, nk if (hc(kk) > hc(idx_thick)) idx_thick = kk end do hc(idx_thick) = hc(idx_thick) + donate end if ! --- step 1: layer potential densities on the compacted column --- do kk = 1, nk src = mapping(kk) rhoc(kk) = eos_density_point(eos, t_col(src), s_col(src), p_ref) end do ! --- step 1b (HYCOM only): bottom-up density monotonize --- ! Work frame is top-down (index 1 = surface, nk = bed): cap each ! cell by the one below it sweeping bed-up so density is ! non-decreasing downward. Pure RHO (hybrid=.false.) skips this ! and stays bit-identical to the merged RHO regrid. if (hybrid) then do kk = nk - 1, 1, -1 rhoc(kk) = min(rhoc(kk), rhoc(kk + 1)) end do end if ! --- steps 2-4: PPM reconstruct + invert each interior target ! density to an interface depth + monotone interfaces. Shared ! density-space inversion (also used by the DENSITY diagnostic ! remap, `rdb_ocean_diag_fills`) — single source of truth for ! the bracket + fixed-iter-Newton solve. Copy the interior ! targets into a stack array so the device call passes a whole ! fixed-size local (no derived-type section descriptor in the ! hot per-column kernel). For HYCOM the rhoc fed in was ! monotonized above; the z* floor below then lifts the result. do kk = 1, nz - 1 rtgt(kk) = rho_target(kk) end do call invert_density_targets(nk, hc, rhoc, nz - 1, rtgt, z_new) ! --- step 4b (HYCOM only): z* nominal-floor sweep --- ! Surface-side minimum-depth floor on the isopycnal interfaces ! (MOM6 `build_hycom1_column`). Walk interfaces from the surface ! down, accumulating the nominal z* depth and pushing any ! too-shallow interface DOWN to it (clamped to the column bottom); ! deep interfaces already below the floor are untouched. ! stretching = col_extent/total_h = (H+η)/H, the z* stretch. The ! nominal increment is a THICKNESS IN METRES — the z* coordinate ! resolution (`floor_dz`, or uniform `h_nominal`) — so the band is ! the same depth range in every column and a shallow column clamps ! its deeper interfaces onto the bed, exactly like a z* grid. Only ! the `HYCOM_FLOOR_SIGMA` fallback accumulates the column FRACTION ! `dsig·H` (a sigma floor: k/nz of every column's depth). Both ! tables are bottom-up (index nz = surface layer); the work layer ! above interface kk maps to bottom-up index nz-kk+2. The floor can ! break monotonicity, so re-monotonize after it. Pure RHO ! (hybrid=.false.) skips this and is bit-identical. if (hybrid) then h_ref_col = total_h(i, j) if (h_ref_col > 0.0_wp) then stretching = col_extent/h_ref_col else stretching = 1.0_wp end if nominal_z = 0.0_wp do kk = 2, nz + 1 select case (floor_mode) case (HYCOM_FLOOR_PROFILE) nominal_z = nominal_z + floor_dz(nz - kk + 2)*stretching case (HYCOM_FLOOR_UNIFORM) nominal_z = nominal_z + h_nominal*stretching case default nominal_z = nominal_z + dsig(nz - kk + 2)*h_ref_col*stretching end select if (z_new(kk) < nominal_z) z_new(kk) = nominal_z if (z_new(kk) > col_extent) z_new(kk) = col_extent end do do kk = 2, nz + 1 if (z_new(kk) < z_new(kk - 1)) z_new(kk) = z_new(kk - 1) end do end if do kk = 1, nz h_new(kk) = z_new(kk + 1) - z_new(kk) end do ! --- step 5: MOM6 min-thickness inflation (floor h_floor_eff) --- ns = 0 do kk = 1, nz if (h_new(kk) > h_floor_eff) ns = ns + 1 end do if (ns == nz) then ! all OK else if (ns == 0) then do kk = 1, nz h_new(kk) = h_floor_eff end do else total_need = 0.0_wp do kk = 1, nz if (h_new(kk) <= h_floor_eff) then total_need = total_need + (h_floor_eff - h_new(kk)) h_new(kk) = h_floor_eff end if end do ! debit the single thickest layer once idx_thick = 1 thick_max = h_new(1) do kk = 2, nz if (h_new(kk) > thick_max) then thick_max = h_new(kk) idx_thick = kk end if end do if (thick_max - total_need >= h_floor_eff) then h_new(idx_thick) = h_new(idx_thick) - total_need else ! The thickest layer alone cannot pay (several comparably thin ! survivors): debit EVERY above-floor layer in proportion to its ! excess over the floor. `col_extent >= nz*h_floor_eff` (guard ! above) makes the total excess >= total_need, so every layer ! stays >= h_floor_eff and the sum is unchanged to round-off. ! The single-layer debit used to drive the thickest one below ! the floor, or negative. excess = 0.0_wp do kk = 1, nz if (h_new(kk) > h_floor_eff) excess = excess + (h_new(kk) - h_floor_eff) end do if (excess > 0.0_wp) then frac = total_need/excess do kk = 1, nz if (h_new(kk) > h_floor_eff) then h_new(kk) = h_new(kk) - frac*(h_new(kk) - h_floor_eff) end if end do end if end if end if ! --- assignment: FLIP top-down working -> bottom-up state --- do k = 1, nz target_h(i, j, k) = h_new(nz - k + 1) end do end subroutine ocean_vcoord_rho_target_column pure function parse_ocean_vcoord_type(name) result(code) !! Ocean-path wrapper around the canonical `parse_vcoord_type` !! in `rdb_vcoord`. Pins the unrecognised-string fallback to !! `VCOORD_EULERIAN_Z` — the ocean path's "do nothing" default, !! distinct from the coastal path's `VCOORD_SIGMA` fallback. !! Kept as a thin name-preserving wrapper so the ocean-only !! semantic (fallback choice) is visible at the call site. character(len=*), intent(in) :: name integer :: code code = parse_vcoord_type(name, default_code=VCOORD_EULERIAN_Z) end function parse_ocean_vcoord_type pure subroutine invert_density_targets(nk, hc, rhoc, n_int, rho_tgt, z_new) !! Density-space interface inversion — the single source of truth !! for the RHO vcoord regrid (`ocean_vcoord_compute_target_h_rho_impl`) !! AND the DENSITY diagnostic remap (`rdb_ocean_diag_fills`). !! !! Given a TOP-DOWN compacted column of `nk` layers with thicknesses !! `hc(1:nk)` (index 1 = surface, depth positive-down) and layer-mean !! potential densities `rhoc(1:nk)`, plus `n_int` monotone-increasing !! interior target densities `rho_tgt(1:n_int)`, return the `n_int + 2` !! interface depths `z_new(1:n_int+2)` (top-down, 0 .. column total), !! monotone non-decreasing. `z_new(1) = 0` (surface), `z_new(n_int+2)` !! = Σ hc (bed); interior interfaces 2..n_int+1 invert each target. !! !! Algorithm (clean-room from Bleck 2002 / White & Adcroft 2008; !! MOM6 `coord_rho` is the behavioural oracle): !! 1. PPM (Colella-Woodward monotone) reconstruction of the density !! profile over `hc`. !! 2. Per interior target: light boundary → surface; discontinuous- !! jump sweep; dense boundary → bed; else fixed-`NR_ITERS`-iter !! Newton on `xi in [0,1]` (convergence on `|delta| < NR_TOL` !! AFTER `xi += delta`; zero-gradient `NR_OFFSET` escape at both !! ends; masked fallback to the previous interface on no-bracket). !! 3. Monotone non-decreasing clamp on the interfaces. !! !! Caller supplies `nk >= 2`. The RHO regrid kernel pre-compacts !! vanished layers (so `nk` is the surviving count) and fast-paths !! `nk <= 1` upstream; the DENSITY diagnostic remap passes the full !! `nz` column, with every vanished layer given zero thickness and the !! density of its nearest live neighbour, so the PPM edges it touches !! are the live layer's own value (no compaction, same effect on the !! inversion). Lightest !! target maps to the surface (index 2), densest to the bed (the !! surface→bed ordering the callers FLIP into the bottom-up state). !$acc routine seq integer, intent(in) :: nk integer, intent(in) :: n_int real(wp), intent(in) :: hc(nk) real(wp), intent(in) :: rhoc(nk) real(wp), intent(in) :: rho_tgt(n_int) real(wp), intent(out) :: z_new(n_int + 2) integer :: kk, ii, k real(wp) :: rhoL(NZ_STACK_MAX), rhoR(NZ_STACK_MAX) real(wp) :: edge(NZ_STACK_MAX + 1), z_old(NZ_STACK_MAX + 1) real(wp) :: tgt, lo, hi, q6, xi, fval, df, delta, grad, ww real(wp) :: rho_light, rho_dense, dd, six logical :: placed ! --- step 1: PPM (Colella-Woodward) edges on the compacted column --- edge(1) = rhoc(1) edge(nk + 1) = rhoc(nk) do kk = 2, nk ww = hc(kk - 1) + hc(kk) ! Floor guards an uncompacted caller (the DENSITY diagnostic remap ! passes the raw column) where a vanished layer pair sums to ~0; ! the RHO regrid pre-compacts so ww is always >> the floor there ! (this branch is inert for it — bit-identical). if (ww < 1.0e-30_wp) ww = 1.0e-30_wp edge(kk) = (hc(kk)*rhoc(kk - 1) + hc(kk - 1)*rhoc(kk))/ww end do do kk = 1, nk rhoL(kk) = edge(kk) rhoR(kk) = edge(kk + 1) ! CW monotonic limiter if ((rhoR(kk) - rhoc(kk))*(rhoc(kk) - rhoL(kk)) <= 0.0_wp) then rhoL(kk) = rhoc(kk) rhoR(kk) = rhoc(kk) else dd = rhoR(kk) - rhoL(kk) six = 6.0_wp*(rhoc(kk) - 0.5_wp*(rhoL(kk) + rhoR(kk))) if (dd*six > dd*dd) then rhoL(kk) = 3.0_wp*rhoc(kk) - 2.0_wp*rhoR(kk) else if (dd*six < -dd*dd) then rhoR(kk) = 3.0_wp*rhoc(kk) - 2.0_wp*rhoL(kk) end if end if end do ! Cumulative compacted-grid interface depths (top-down, 0..H). z_old(1) = 0.0_wp do kk = 1, nk z_old(kk + 1) = z_old(kk) + hc(kk) end do rho_light = rhoL(1) rho_dense = rhoR(nk) ! --- step 2: invert each interior target interface (top-down) --- z_new(1) = 0.0_wp z_new(n_int + 2) = z_old(nk + 1) ! total compacted depth (= col extent) do kk = 2, n_int + 1 ! interior interfaces tgt = rho_tgt(kk - 1) if (tgt <= rho_light) then z_new(kk) = 0.0_wp ! lighter than column -> surface else if (tgt >= rho_dense) then z_new(kk) = z_old(nk + 1) ! denser than column -> bed else placed = .false. do ii = 1, nk lo = rhoL(ii) hi = rhoR(ii) ! Discontinuous jump at the TOP interface of cell ii ! (between cell ii-1's right edge and cell ii's left ! edge). Checked FIRST and independent of whether the ! cells are limiter-flattened — a 2-layer column flattens ! both cells, and a target inside the jump must still ! land on the interface (MOM6 coord_rho behaviour). if (ii > 1) then if (rhoR(ii - 1) <= tgt .and. tgt <= rhoL(ii)) then z_new(kk) = z_old(ii) placed = .true. exit end if end if if (lo == hi) then ! Flat (limiter-collapsed) cell: only an exact match ! places here; otherwise advance to the next cell. if (abs(tgt - lo) < NR_TOL) then z_new(kk) = z_old(ii) placed = .true. exit end if cycle end if if ((lo - tgt)*(hi - tgt) <= 0.0_wp) then ! Newton on xi in [0,1]: ppm(xi) - tgt = 0 q6 = 6.0_wp*rhoc(ii) - 3.0_wp*(lo + hi) xi = 0.5_wp do k = 1, NR_ITERS fval = lo + xi*((hi - lo) + q6*(1.0_wp - xi)) - tgt df = (hi - lo) + q6*(1.0_wp - 2.0_wp*xi) if (abs(df) > 1.0e-30_wp) then delta = -fval/df else delta = 0.0_wp end if xi = xi + delta ! clamp inside the iteration; zero-gradient nudge if (xi < 0.0_wp) then xi = 0.0_wp grad = (hi - lo) + q6 ! d(ppm)/dxi at xi=0 if (abs(grad) < 1.0e-30_wp) xi = NR_OFFSET else if (xi > 1.0_wp) then xi = 1.0_wp grad = (hi - lo) - q6 ! d(ppm)/dxi at xi=1 if (abs(grad) < 1.0e-30_wp) xi = 1.0_wp - NR_OFFSET end if if (abs(delta) < NR_TOL) exit end do z_new(kk) = z_old(ii) + xi*hc(ii) placed = .true. exit end if end do if (.not. placed) then ! masked fallback: previous interface (no FATAL in a DC) z_new(kk) = z_new(kk - 1) end if end if end do ! --- step 3: monotone non-decreasing interfaces --- do kk = 2, n_int + 2 if (z_new(kk) < z_new(kk - 1)) z_new(kk) = z_new(kk - 1) end do end subroutine invert_density_targets pure function ocean_vcoord_bytes(this) result(nbytes) !! Counted allocatable footprint of the vertical coordinate slot !! (0 when unallocated). One arr_bytes term per array — add a !! term here when a new allocatable joins the type. class(ocean_vcoord_t), intent(in) :: this integer(int64) :: nbytes nbytes = arr_bytes(this%dsig) & + arr_bytes(this%z_ref_global) & + arr_bytes(this%target_h) & + arr_bytes(this%z_ref) & + arr_bytes(this%z_top) & + arr_bytes(this%z_fixed_zi) & + arr_bytes(this%z_fixed_dz) & + arr_bytes(this%rho_target) & + arr_bytes(this%remap_total_h) & + arr_bytes(this%remap_h_ref) & + arr_bytes(this%remap_h_old) & + arr_bytes(this%remap_conc_t) & + arr_bytes(this%remap_conc_s) end function ocean_vcoord_bytes end module rdb_ocean_vcoord