ocean_state_seed_from_cfg Subroutine

   subroutine ocean_state_seed_from_cfg(state, grid, cfg, ierr, injected_b, injected_b_convention, &
                                        periodic_x, periodic_y)
      !! Populate the ocean prognostic state with an analytical IC
      !! derived from cfg scalars.  Bathymetry is set per `cfg%ocean%topo%topo_config`
      !! (`"flat"` → uniform `ocean_max_depth`; `"spoon"` → MOM6 spoon
      !! shape) — UNLESS `injected_b` is present, in which case it
      !! overrides `topo_config` entirely (P2.5 pre-create geometry
      !! injection).  Water column thickness `h = b` so the free surface
      !! starts at SSH = 0.  Layers split the local depth evenly
      !! (`h_layer(i,j,k) = b(i,j) / nz_ml`).  Tracers carry per-layer
      !! `h * Tr` at the configured initial T/S.  Velocities zeroed.
      !!
      !! Caller must have run `state%init(grid)` first so every slot's
      !! allocations are in place; this routine only writes into them
      !! on the host.  Run before `ocean_state_enter_data` so the GPU
      !! mapping captures the seeded values.
      !!
      !! Every position-aware formula fill below reads its global offsets
      !! and global extents from `grid` (`i_offset_global` /
      !! `j_offset_global`, `nx_global` / `ny_global`), so each rank seeds
      !! its own window on ONE global analytical IC.  There is deliberately
      !! no `decomp` argument and no optional offset quartet: an optional
      !! that defaults to the LOCAL extents is silently wrong under MPI, and
      !! forgetting to pass it shipped three separate bugs.
      type(ocean_state_t), intent(inout) :: state
      type(hgrid_t), intent(in) :: grid
      type(config_t), intent(in) :: cfg
      integer, intent(out), optional :: ierr
         !! Non-zero on an initial-condition configuration/data failure
         !! when present; absent behaves as today (`error stop`).
      real(wp), intent(in), optional :: injected_b(:, :)
         !! P2.5 pre-create geometry injection: interior-sized
         !! `(grid%nx_phys, grid%ny_phys)` bathymetry, RAW in the caller's
         !! own sign convention (see `injected_b_convention`) — REQUIRES
         !! `injected_b_convention`.  When present, this OVERRIDES
         !! `cfg%ocean%topo%topo_config` entirely: sign-normalised to
         !! positive-down depth + wet-fraction validated
         !! (`bathymetry_normalise_sign`), written into the interior of
         !! `state%barotropic%b`, and ghost-filled by constant
         !! extrapolation (`bathymetry_fill_ghosts_array`) — mirroring
         !! `topo_config='file'` rather than trusting the caller
         !! (CLAUDE.md formula-bathymetry ghost-fill gotcha).  Everything
         !! downstream in this routine (h=b, layer split, wet mask, T/S
         !! seed, IC overlays, the ALE z_ref table) then runs unchanged
         !! against the injected bathymetry.
      integer, intent(in), optional :: injected_b_convention
         !! REQUIRED alongside `injected_b` — no default (D6.2: sign is
         !! the single most dangerous argument in the geometry API). One
         !! of `BATHY_CONVENTION_DEPTH_POSITIVE_DOWN` /
         !! `_HEIGHT_POSITIVE_UP` (`rdb_ocean_bathymetry_inject`).
      logical, intent(in), optional :: periodic_x
         !! Grid topology the static geometry is wrapped with (see
         !! `seed_wrap_static_2d`).  Absent ⇒ derived from the
         !! `&ocean_bc_nml` west/east tags, exactly as `configure_ocean_bc`
         !! derives `bc%periodic_x`.  The engine passes it when a staged
         !! topology (`rdb_ocean_stage_topology`) overrides the tags.
      logical, intent(in), optional :: periodic_y
         !! As `periodic_x`, for south/north.

      integer :: nz_ml, idx_S, idx_T, i, j, k, nx, ny, local_ierr
      logical :: per_x, per_y, north_fold
      logical :: seeded_on_eta0_target
         !! The FOURTH or FIFTH `h_layer` seed branch (`zstar_full` under
         !! `zfixed_closed_faces`; MOM6 `zstar` under the knob or zinit) was
         !! taken ⇒ establish I1′ after the tracer IC.
      real(wp), allocatable :: water(:, :)
         !! Reference water-column thickness the IC seeds work on:
         !! `b − z_draft` under an ice shelf, a byte copy of `b`
         !! otherwise.  Host-only setup scratch, released on return.
      real(wp), allocatable :: eta_trim(:, :)
         !! Initial free-surface anomaly of the TRIMMED cavity IC
         !! (`&ocean_cavity_dyn_nml trim_ic_for_p_surf`); 0 when off.
      real(wp), allocatable :: h_col(:, :)
         !! Initial water column `water + eta_trim` the layer split seeds
         !! from; a byte copy of `water` when the trim is off.
      logical :: trim_ic
         !! `&ocean_cavity_dyn_nml trim_ic_for_p_surf` under an active cavity.

      nz_ml = state%multilayer%nz_ml
      idx_S = state%multilayer%idx_salinity
      idx_T = state%multilayer%idx_temperature
      nx = size(state%barotropic%b, 1)
      ny = size(state%barotropic%b, 2)

      ! Grid topology of the static geometry (the same rule
      ! `ocean_bc_state_init` applies to the tags, read here because the
      ! seed runs BEFORE `configure_ocean_bc`).
      per_x = ocean_bc_type_from_string(cfg%ocean%bc%west) == OBC_PERIODIC .and. &
              ocean_bc_type_from_string(cfg%ocean%bc%east) == OBC_PERIODIC
      per_y = ocean_bc_type_from_string(cfg%ocean%bc%south) == OBC_PERIODIC .and. &
              ocean_bc_type_from_string(cfg%ocean%bc%north) == OBC_PERIODIC
      if (present(periodic_x)) per_x = periodic_x
      if (present(periodic_y)) per_y = periodic_y
      north_fold = ocean_bc_type_from_string(cfg%ocean%bc%north) == OBC_TRIPOLAR_FOLD

      ! Bathymetry.  `slope_scale` (the spoon/seamount length scale) is a
      ! metres knob; `set_bathymetry_*` works in GRID coordinate units,
      ! which are metres on Cartesian grids but DEGREES on spherical/
      ! curvilinear ones.  Convert here so the length scale is consistent
      ! with the grid coordinate (without it, a metres scale against a
      ! degrees position collapses the seamount/spoon to a flat basin).
      if (present(injected_b)) then
         ! P2.5 pre-create geometry injection — bypasses topo_config
         ! entirely.  `injected_b_convention` has no default (D6.2).
         if (.not. present(injected_b_convention)) then
            call fail("ocean_state_seed_from_cfg: injected_b requires "// &
                      "injected_b_convention (no default — sign is the single most "// &
                      "dangerous argument; see D6.2 in the Python API design notes)", &
                      ierr, OCEAN_STATUS_ERR_BAD_SHAPE)
            return
         end if
         if (size(injected_b, 1) /= grid%nx_phys .or. size(injected_b, 2) /= grid%ny_phys) then
            call fail("ocean_state_seed_from_cfg: injected_b shape mismatch — got ("// &
                      to_string(size(injected_b, 1))//","//to_string(size(injected_b, 2))// &
                      "), expected the physical interior ("//to_string(grid%nx_phys)//","// &
                      to_string(grid%ny_phys)//")", ierr, OCEAN_STATUS_ERR_BAD_SHAPE)
            return
         end if
         block
            real(wp) :: b_local(grid%nx_phys, grid%ny_phys)
            integer :: ng
            b_local = injected_b
            call bathymetry_normalise_sign(b_local, injected_b_convention, ierr=local_ierr)
            if (local_ierr /= OCEAN_STATUS_OK) then
               if (present(ierr)) then
                  ierr = local_ierr
                  return
               end if
               error stop "ocean_state_seed_from_cfg: injected_b sign normalisation failed"
            end if
            ng = grid%nghost
            state%barotropic%b(ng + 1:ng + grid%nx_phys, ng + 1:ng + grid%ny_phys) = b_local
            call bathymetry_fill_ghosts_array(state%barotropic%b, grid)
         end block
      else
         select case (trim(cfg%ocean%topo%topo_config))
         case ("spoon")
            call set_bathymetry_spoon(state%barotropic%b, grid, &
                                      cfg%ocean%topo%max_depth, cfg%ocean%topo%edge_depth, &
                                      topo_length_to_grid_units(cfg%ocean%topo%slope_scale, &
                                                                cfg%ocean%grid%grid_config, &
                                                                cfg%ocean%grid%rad_earth))
         case ("seamount")
            call set_bathymetry_seamount(state%barotropic%b, grid, &
                                         cfg%ocean%topo%max_depth, cfg%ocean%topo%edge_depth, &
                                         topo_length_to_grid_units(cfg%ocean%topo%slope_scale, &
                                                                   cfg%ocean%grid%grid_config, &
                                                                   cfg%ocean%grid%rad_earth))
         case ("neverworld2")
            call set_bathymetry_neverworld2(state%barotropic%b, grid, &
                                            cfg%ocean%topo%max_depth, &
                                            cfg%ocean%topo%nl_continent_amp, &
                                            cfg%ocean%topo%nl_roughness_amp, &
                                            cfg%ocean%topo%nl_min_depth)
         case ("island")
            ! Flat basin at `max_depth` with a central LAND square (depth
            ! forced below LAND_DEPTH_THRESHOLD).  Produces the interior
            ! land block the static land-mask validation needs.  The square
            ! spans the central `island_frac` fraction of the physical
            ! domain in each direction (default via edge_depth/slope_scale
            ! re-use; see set_bathymetry_island).
            call set_bathymetry_island(state%barotropic%b, grid, &
                                       cfg%ocean%topo%max_depth, &
                                       cfg%ocean%topo%slope_scale)
         case ("double_drake")
            ! Ferreira et al. (2010) two-wall supercontinent with a
            ! reentrant southern channel — the land-mask showcase.
            call set_bathymetry_double_drake(state%barotropic%b, grid, &
                                             cfg%ocean%topo%max_depth, &
                                             cfg%ocean%topo%slope_scale)
         case ("isomip_plus")
            ! MISMIP+/ISOMIP+ analytic bedrock (Asay-Davis et al. 2016,
            ! Eqs. 1-4 + Table 1).  `max_depth` IS the deep clip
            ! (`-z_b,deep`, protocol 720 m) and `x_origin` places the
            ! model's west edge on the paper's absolute x axis (protocol
            ! 320 km).  The formula is written in METRES, so the setter
            ! is handed metres-per-grid-unit rather than a converted
            ! length: `topo_length_to_grid_units(1, ...)` is grid units
            ! per metre, and this is its reciprocal (exactly 1 on a
            ! Cartesian grid).
            call set_bathymetry_isomip_plus(state%barotropic%b, grid, &
                                            cfg%ocean%topo%max_depth, &
                                            cfg%ocean%topo%x_origin, &
                                            1.0_wp/topo_length_to_grid_units(1.0_wp, &
                                                                             cfg%ocean%grid%grid_config, &
                                                                             cfg%ocean%grid%rad_earth))
         case ("file")
            ! Real bathymetry from NetCDF.  File must be pre-projected onto
            ! the model's Cartesian grid (matching nx_phys × ny_phys); the
            ! loader fills ghost rows by constant extrapolation.  Sign
            ! convention: `b` is bottom depth positive-down (consistent with
            ! the spoon + flat branches).  GEBCO + ETOPO ship elevation
            ! positive-up — the preprocessing script should flip the sign
            ! before writing the NetCDF.
            if (len_trim(cfg%bathymetry_file) == 0) then
               call fail("ocean_state_seed_from_cfg: topo_config='file' requires "// &
                         "bathymetry_file", ierr, OCEAN_STATUS_ERR_IC_SEED)
               return
            end if
#ifndef RDB_NO_NETCDF
            if (present(ierr)) then
               call load_bathymetry_into_array(trim(cfg%bathymetry_file), &
                                               state%barotropic%b, grid, ierr=local_ierr)
               if (local_ierr /= OCEAN_STATUS_OK) then
                  ierr = local_ierr
                  return
               end if
            else
               call load_bathymetry_into_array(trim(cfg%bathymetry_file), &
                                               state%barotropic%b, grid)
            end if
#else
            ! File-loaded bathymetry needs the NetCDF reader in `rdb_bathymetry`,
            ! which isn't compiled in when RDB_ENABLE_NETCDF=OFF.  Fail fast
            ! with a descriptive error so the user knows to rebuild with NetCDF
            ! enabled rather than hitting a less-obvious runtime issue.
            call fail("ocean_state_seed_from_cfg: topo_config='file' requires "// &
                      "RDB_ENABLE_NETCDF=ON at build time (the NetCDF-backed "// &
                      "rdb_bathymetry reader is needed to load the file).", ierr, OCEAN_STATUS_ERR_IO)
            return
#endif
         case default  ! "flat"
            state%barotropic%b = cfg%ocean%topo%max_depth
         end select
      end if

      ! Make the bathymetry SEAM-CONSISTENT before anything reads it.  The
      ! file loader and the injected-array path fill the ghosts by constant
      ! extrapolation, and a formula setter evaluates its formula OUTSIDE
      ! the domain; on a periodic or folded edge neither is the value the
      ! seam needs.  Every field seeded below (the water column, h_layer,
      ! the wet mask, the tracers, the zstar_full z_ref table) and every
      ! setup-time consumer before the engine's first halo pass (the PGF's
      ! own bathymetry copy, bt_H_ref, the wet/dry and sponge setup) reads
      ! these ghosts — a stale seam ghost is a 1.9 m/s jet on the seam face
      ! of the 1-degree global grid within three hours.  The engine still
      ! re-wraps + halo-exchanges `b` later (the only fill a DECOMPOSED
      ! axis can get); on the local axes that is now a no-op.
      call seed_wrap_static_2d(state%barotropic%b, grid, per_x, per_y, north_fold)

      ! ---- Static ice-shelf cavity geometry (&ocean_cavity_dyn_nml, P5.1) ----
      ! ORDERING IS LOAD-BEARING and this is the only place it can go: the
      ! draft must exist before the water column `b − z_draft` seeds the
      ! layer split and the wet mask (grounding), both of which happen a
      ! few lines below, and long before any `configure_ocean_*` pass.
      ! The formula setters fill the FULL array including ghosts; the
      ! periodic/fold re-wrap + halo exchange that `barotropic%b` gets are
      ! applied to `z_draft` alongside it in `rdb_ocean_engine`.
      !
      ! `water` is the reference water-column thickness every seed below
      ! works on.  With the knob off it is a byte copy of `b`, so there is
      ! ONE code path and the default run is bit-identical.
      allocate (water(nx, ny))
      if (state%metrics%use_cavity) then
         call seed_cavity_draft(state, grid, cfg, ierr=local_ierr)
         if (local_ierr /= OCEAN_STATUS_OK) then
            if (present(ierr)) then
               ierr = local_ierr
               return
            end if
            error stop "ocean_state_seed_from_cfg: ice-shelf cavity geometry failed"
         end if
         ! Same seam rule as `b` just above, for the same reason: the
         ! draft and its cover fraction are bathymetry-class geometry, and
         ! the water column below reads their ghosts.
         call seed_wrap_static_2d(state%metrics%z_draft, grid, per_x, per_y, north_fold)
         call seed_wrap_static_2d(state%metrics%cover_frac, grid, per_x, per_y, north_fold)
         call cavity_water_column_impl(water, state%barotropic%b, &
                                       state%metrics%z_draft, nx, ny)
      else
         water = state%barotropic%b
      end if

      ! MOM6 TRIM_IC_FOR_P_SURF (`&ocean_cavity_dyn_nml
      ! trim_ic_for_p_surf`, default off).  The ice LOAD stays the
      ! Boussinesq-isostatic `rho_ref*g*z_draft`; each loaded column's
      ! initial TOP moves to the depth where the displaced water's own
      ! weight equals it, so the MOM6 barotropic split starts at rest (see
      ! `cavity_trim_eta_linear_impl`).  `water` itself — the datum and
      ! the grounding decision — is untouched: the trim is an initial
      ! `bt_eta`, not a geometry change.  The density is the linear EOS
      ! over the analytic zinit profile, both enforced by
      ! `validate_config`, and taken from `cfg` because they are exactly
      ! what `ocean_state_init` copied onto `state%eos`.
      allocate (eta_trim(nx, ny), source=0.0_wp)
      trim_ic = state%metrics%use_cavity .and. cfg%ocean%cavity_dyn%trim_ic_for_p_surf
      if (trim_ic) then
         block
            real(wp) :: rho_surf, drho_dz
            logical :: trim_ok
            rho_surf = cfg%ocean%ic%rho_0 + &
                       cfg%ocean%ic%beta_S*(cfg%ocean%zinit%lin_s_ref - cfg%ocean%ic%S_ref) - &
                       cfg%ocean%ic%alpha_T*(cfg%ocean%zinit%lin_t_ref - cfg%ocean%ic%T_ref)
            drho_dz = cfg%ocean%ic%beta_S*cfg%ocean%zinit%lin_ds_dz - &
                      cfg%ocean%ic%alpha_T*cfg%ocean%zinit%lin_dt_dz
            call cavity_trim_eta_linear_impl(eta_trim, trim_ok, state%metrics%z_draft, &
                                             water, cfg%ocean%cavity_dyn%h_min_cavity, &
                                             cfg%ocean%ic%rho_0, rho_surf, drho_dz, nx, ny)
            if (.not. trim_ok) then
               call fail("ocean_state_seed_from_cfg: &ocean_cavity_dyn_nml "// &
                         "trim_ic_for_p_surf found no admissible trim depth (the "// &
                         "initial density must be positive and stably stratified, "// &
                         "and the trimmed column must stay non-empty)", &
                         ierr, OCEAN_STATUS_ERR_IC_SEED)
               return
            end if
            call logger%info("ocean_cavity: trimmed the initial column under the ice "// &
                             "to the load (MOM6 TRIM_IC_FOR_P_SURF): eta in ["// &
                             to_string(minval(eta_trim))//", "// &
                             to_string(maxval(eta_trim))//"] m")
         end block
      end if
      allocate (h_col(nx, ny))
      if (trim_ic) then
         h_col = water + eta_trim
      else
         h_col = water
      end if

      ! Water column thickness h = b − z_draft → the free-surface anomaly
      ! `bt_eta = Σ h_layer − bt_H_ref` starts at zero (SSH = 0 without a
      ! cavity; the loaded equilibrium under one) — or at `eta_trim` under
      ! a trimmed cavity IC.
      state%barotropic%h = h_col
      state%barotropic%u_face_x = 0.0_wp
      state%barotropic%v_face_y = 0.0_wp
      state%barotropic%hu_face_x = 0.0_wp
      state%barotropic%hv_face_y = 0.0_wp

      ! Per-column layer thickness = depth / nz_ml.  Pass flat
      ! allocatables into a `_impl` helper so the `do concurrent` body
      ! reads/writes plain arrays — registry / array-of-DT indirection
      ! crashes NVHPC's device codegen with CUDA_ERROR_ILLEGAL_ADDRESS.
      ! With wet/dry ON, floor the seed to 2·H_VANISHED so an emerged
      ! intertidal rest-bed (b < 0 ⇒ now a LIVE wet_mask=1 column) never
      ! seeds a negative layer / negative barotropic depth D = Σ h_layer,
      ! and each seeded layer clears the strict `h_old > H_VANISHED`
      ! remap-drain gate (so the seeded S/T survives the first regrid).
      ! Tracers are seeded from THIS floored h_layer just below, so hTr
      ! stays consistent (const·2·H_VANISHED, finite T/S, no negative mass).
      ! `thickness_config = "uniform_z"` swaps the local-depth even split for
      ! MOM6's uniform-z-interface seed (flat resting isopycnals under a
      ! horizontally-uniform density stack).  Validated as mutually exclusive
      ! with wet/dry, so the two branches never both need the emerged-column
      ! floor.  Default "sigma" ⇒ byte-identical to the pre-knob path.
      !
      ! THIRD BRANCH — `VCOORD_Z_FIXED` under a cavity.  The running
      ! coordinate there is quasi-geopotential with inert fillers inside
      ! the ice, so a sigma-style seed is NOT on the coordinate: the very
      ! first ALE remap would relamp the whole column in one step, and a
      ! T/S profile that `&ocean_zinit_nml source="linear"` made exactly
      ! linear in geopotential z would come back through the PPM boundary
      ! closure NOT exactly linear — column by column, because the draft
      ! (and so the cut) differs column to column.  That difference IS a
      ! horizontal density gradient, i.e. exactly the spurious rest
      ! current this coordinate exists to remove.  So seed `h_layer`
      ! directly FROM the target (`η = 0`, which is the cavity datum's
      ! own resting state) and let the zinit overlay evaluate T/S at
      ! those layer centres: exact by construction, first remap an
      ! identity, no step-1 regrid shock.
      !
      ! The same holds WITHOUT a cavity whenever the T/S come from the
      ! GEOPOTENTIAL `&ocean_zinit_nml` overlay: a sigma-style seed puts
      ! layer `k` of a column of depth `H` at `(nz-k+1/2)*H/nz`, the zinit
      ! profile is evaluated THERE, and only the step-1 regrid moves it onto
      ! the `z_fixed` layers.  Anything snapshotted from the seed in between
      ! is then indexed on the wrong layers -- the `&ocean_sponge_nml
      ! target_source = "ic"` reference is (`ocean_sponge_snapshot_reference`
      ! runs right after this seed), so the sponge relaxed layer `k` of the
      ! 50-level tanh grid toward the zinit value hundreds of metres deeper,
      ! and toward a DIFFERENT depth in every column of a different `H`: a
      ! grid-scale, bathymetry-following horizontal density forcing across
      ! the whole band (measured: the band's surface layers pulled 5-6 degC
      ! cold, domain-mean T -0.55 degC in 5 days on the coastal-noise box).
      ! Seeding on the target makes the overlay exact and the snapshot
      ! consistent.  `z_fixed` WITHOUT zinit keeps the sigma-style seed: its
      ! `&tracer_nml` IC is defined PER LAYER INDEX, so moving the layers
      ! would change what that IC means.
      !
      ! FOURTH BRANCH — `VCOORD_ZSTAR_FULL` under `&vcoord_nml
      ! zfixed_closed_faces`.  The closed-face mask is built from the
      ! `ZSTAR_FULL` target at `η = 0`; a sigma-style seed is not on that
      ! coordinate, so step 1 would run the FULL sigma pressure gradient on
      ! the `b/nz` stack (on the 1-degree Southern Ocean, 0.67 m/s from
      ! rest in one step, at faces the mask does not even see) before the
      ! first regrid moved the layers, and a `target_source = "ic"` sponge
      ! would snapshot T/S on the wrong layers — the `z_fixed` seed bug
      ! above.  So seed `h_layer` from the SAME target the mask is built
      ! from (`ocean_vcoord_eta0_target`, after laying the `z_ref` table
      ! it walks), with or without zinit: a per-index `&tracer_nml` IC
      ! then means "per coordinate layer", which is the only reading under
      ! which the mask's live/filler pattern is the IC's.  Knob-gated, so
      ! every existing `zstar_full` namelist keeps its sigma-style seed
      ! byte for byte.
      !
      ! FIFTH BRANCH — `VCOORD_ZSTAR` (MOM6 z*) under `zfixed_closed_faces`
      ! OR `&ocean_zinit_nml`.  The same reasons as the fourth (the mask's
      ! pattern is the `η = 0` target's) and the third (a z-level IC is
      ! exact only at the running coordinate's own layer centres; a
      ! `target_source = "ic"` sponge snapshots what is seeded), on the
      ! `η = 0` z* target, which is the `z_fixed` one with `z_top = 0`.  A
      ! per-layer-index `&tracer_nml` IC without the knob keeps the
      ! sigma-style seed and is moved onto z* by the first regrid, as on
      ! `z_fixed`.  z* is refused under a cavity, so `eta_trim ≡ 0` here.
      seeded_on_eta0_target = .false.
      if (trim(cfg%thickness_config) == "uniform_z") then
         call seed_h_layer_uniform_z_impl(state%multilayer%h_layer, &
                                          water, nz_ml, &
                                          cfg%ocean%topo%max_depth, &
                                          cfg%ocean%isopycnal%angstrom_h)
      else if ((state%metrics%use_cavity .or. cfg%ocean%zinit%enable) .and. &
               parse_ocean_vcoord_type(cfg%vcoord_type) == VCOORD_Z_FIXED .and. &
               cfg%ocean%topo%max_depth > 0.0_wp) then
         block
            real(wp) :: h_min_seed
            real(wp), allocatable :: z_top_seed(:, :)
            ! Column-top depth: the ice draft under a cavity, else `z = 0`
            ! (`metrics%z_draft` is only a `(1,1)` placeholder then).
            if (state%metrics%use_cavity) then
               z_top_seed = state%metrics%z_draft
            else
               allocate (z_top_seed(nx, ny), source=0.0_wp)
            end if
            ! `zstar_h_min` comes off the SLOT, not off `cfg`: there is
            ! one source of truth for the filler thickness and it is the
            ! one the running target builder will use.  `engine_setup`
            ! copies both `zstar_*` knobs onto the slot immediately BEFORE
            ! this seed (it has to — the seed's own tail calls
            ! `vcoord%build_zref_full`, which reads them), so the value is
            ! already the namelist's.  Fall back to `cfg` only for a
            ! caller that seeds a state whose vcoord slot was never
            ! initialised.
            h_min_seed = cfg%zstar_h_min
            if (state%vcoord%is_init) h_min_seed = state%vcoord%zstar_h_min
            ! A stretched nominal profile (`&vcoord_nml z_fixed_profile`)
            ! is installed on the slot by `engine_setup` before this seed,
            ! for the same one-source-of-truth reason as `zstar_h_min`.
            if (state%vcoord%is_init .and. state%vcoord%z_fixed_use_profile) then
               call ocean_vcoord_z_fixed_target(state%multilayer%h_layer, water, eta_trim, &
                                                z_top_seed, nx, ny, nz_ml, &
                                                cfg%ocean%topo%max_depth/real(nz_ml, wp), &
                                                .true., state%vcoord%z_fixed_zi, &
                                                state%vcoord%z_fixed_dz, h_min_seed)
            else
               call ocean_vcoord_z_fixed_target_uniform(state%multilayer%h_layer, water, &
                                                        eta_trim, z_top_seed, &
                                                        nx, ny, nz_ml, &
                                                        cfg%ocean%topo%max_depth/real(nz_ml, wp), &
                                                        h_min_seed)
            end if
         end block
      else if (cfg%zfixed_closed_faces .and. state%vcoord%is_init .and. &
               state%vcoord%coord_type == VCOORD_ZSTAR_FULL .and. &
               parse_ocean_vcoord_type(cfg%vcoord_type) == VCOORD_ZSTAR_FULL) then
         call state%vcoord%build_zref_full(state%barotropic%b)
         call ocean_vcoord_eta0_target(state%vcoord, state%multilayer%h_layer, &
                                       water, nx, ny, nz_ml)
         seeded_on_eta0_target = .true.
      else if ((cfg%zfixed_closed_faces .or. cfg%ocean%zinit%enable) .and. &
               state%vcoord%is_init .and. &
               state%vcoord%coord_type == VCOORD_ZSTAR .and. &
               parse_ocean_vcoord_type(cfg%vcoord_type) == VCOORD_ZSTAR .and. &
               state%vcoord%z_fixed_h_ref > 0.0_wp) then
         call ocean_vcoord_eta0_target(state%vcoord, state%multilayer%h_layer, &
                                       water, nx, ny, nz_ml)
         seeded_on_eta0_target = .true.
      else
         call seed_h_layer_uniform_impl(state%multilayer%h_layer, &
                                        h_col, nz_ml, &
                                        apply_wetdry_floor=cfg%ocean%wetdry%enable)
      end if
      ! Keep the barotropic water-column prognostic non-negative on the
      ! same emerged band (it is inert in the multilayer dyn-core — SSH is
      ! Σ h_layer − b — but is restart-registered as `bt_h`; floor it so a
      ! checkpoint never carries a negative rest depth).  Match the layer
      ! floor exactly: on an emerged column Σ h_layer = nz·2·H_VANISHED, so
      ! floor bt_h to the SAME value (not 0) — otherwise a restart that
      ! re-derives D from bt_h would disagree with Σ h_layer by nz·2·H_VANISHED
      ! on every emerged column.  Knob-off ⇒ h = b (byte-identical).
      if (cfg%ocean%wetdry%enable) then
         state%barotropic%h = max(water, &
                                  real(nz_ml, wp)*2.0_wp*H_VANISHED)
      end if
      state%multilayer%u_face_x_layer = 0.0_wp
      state%multilayer%v_face_y_layer = 0.0_wp
      state%multilayer%hu_face_x_layer = 0.0_wp
      state%multilayer%hv_face_y_layer = 0.0_wp

      ! Wet/dry mask: 1.0 where `b >= cutoff`, 0.0 elsewhere (land).
      ! Default cutoff = LAND_DEPTH_THRESHOLD (byte-identical for all
      ! non-wetdry configs).  When wetdry is enabled, pass `land_cutoff =
      ! -land_margin` so intertidal columns (bed above rest MSL but below the
      ! flood headroom) stay wet_mask=1 and the dynamic wd_wet_dyn gate
      ! handles their wetting/drying instead of the static land mask.
      if (cfg%ocean%wetdry%enable) then
         call seed_wet_mask_impl(state%multilayer%wet_mask, water, &
                                 land_cutoff=-cfg%ocean%wetdry%land_margin)
      else if (state%metrics%use_cavity) then
         ! GROUNDING.  A column with less than `h_min_cavity` of water
         ! under the ice is LAND — routed through the SAME wet-mask seed
         ! the bathymetry uses, so the static metric-zeroing land mask
         ! (`configure_ocean_land_mask`) and the finite land-state hold
         ! (`ocean_state_seed_land_cells`) follow for free.  Never a thin
         ! film of water under grounded ice.  Note the cutoff is applied
         ! to `b − z_draft`, so an ordinary land column (`b` below
         ! LAND_DEPTH_THRESHOLD, draft already zeroed there) is land for
         ! the same reason it always was.
         call seed_wet_mask_impl(state%multilayer%wet_mask, water, &
                                 land_cutoff=cfg%ocean%cavity_dyn%h_min_cavity)
      else
         call seed_wet_mask_impl(state%multilayer%wet_mask, water)
      end if

      ! Tracers carry per-layer h*Tr.  Multiply the (now spatially-
      ! varying) h_layer by the configured uniform scalar.
      if (idx_S > 0) then
         block
            real(wp) :: s_layer(nz_ml)
            real(wp) :: dS_dlayer
            logical :: stratify_s
            ! Linear S(z) when both surface + bottom are specified —
            ! EXACT mirror of the temperature branch below, including
            ! the gate (`both /= 0`), the vertical convention (k=1 is
            ! the bed = S_init_bottom, k=nz_ml the surface =
            ! S_init_surface), the single-layer fall-through, and the
            ! seed helper (so ghosts and land columns are filled the
            ! same way: hTr = S(k)·h_layer everywhere, land included,
            ! since h_layer is already 0/floored there).  The profile is
            ! linear in LAYER INDEX, which under the sigma-style
            ! `h_layer = b/nz_ml` seed is linear in layer-centre depth
            ! on every column — so a sloping bed gets the same endpoint
            ! values with a depth-proportional gradient.
            ! Note the stable polarity is the INVERSE of temperature:
            ! dense/salty water belongs at the bed, so a stable haline
            ! column has `S_init_bottom > S_init_surface`.
            ! Both-zero (the default) ⇒ uniform `initial_salinity`,
            ! byte-identical to the pre-knob path.
            stratify_s = (cfg%S_init_surface /= 0.0_wp) .and. &
                         (cfg%S_init_bottom /= 0.0_wp) .and. &
                         (nz_ml > 1)
            if (stratify_s) then
               dS_dlayer = (cfg%S_init_surface - cfg%S_init_bottom)/ &
                           real(nz_ml - 1, wp)
               do k = 1, nz_ml
                  s_layer(k) = cfg%S_init_bottom + dS_dlayer*real(k - 1, wp)
               end do
               call seed_tracer_stratified_impl( &
                  state%multilayer%tracers(idx_S)%hTr, &
                  state%multilayer%h_layer, s_layer, nz_ml)
            else
               call seed_tracer_uniform_impl( &
                  state%multilayer%tracers(idx_S)%hTr, &
                  state%multilayer%h_layer, &
                  cfg%initial_salinity, nz_ml)
            end if
         end block
      end if
      if (idx_T > 0) then
         block
            real(wp) :: t_layer(nz_ml)
            real(wp) :: dT_dlayer
            logical :: stratify
            ! Linear T(z) when both surface + bottom are specified.
            ! Convention: k=1 is the bed (T = T_init_bottom), k=nz_ml is
            ! the surface (T = T_init_surface).  Single-layer case
            ! falls through to T_init_surface (no profile to build).
            stratify = (cfg%T_init_surface /= 0.0_wp) .and. &
                       (cfg%T_init_bottom /= 0.0_wp) .and. &
                       (nz_ml > 1)
            if (stratify) then
               dT_dlayer = (cfg%T_init_surface - cfg%T_init_bottom)/ &
                           real(nz_ml - 1, wp)
               do k = 1, nz_ml
                  t_layer(k) = cfg%T_init_bottom + dT_dlayer*real(k - 1, wp)
               end do
            else
               t_layer = cfg%initial_temperature
            end if
            call seed_tracer_stratified_impl( &
               state%multilayer%tracers(idx_T)%hTr, &
               state%multilayer%h_layer, t_layer, nz_ml)
         end block
      end if

      ! Optional IC overlay applied after bathymetry + uniform h_layer
      ! + analytical T/S are in place.
      ! ierr threaded down ONLY when THIS routine's own ierr is present:
      ! otherwise each seed_*_ic helper must keep reaching its own
      ! `error stop` (specific text) rather than the generic wrapper
      ! message below (P0.1 review F2).
      if (present(ierr)) then
         select case (trim(cfg%ocean%ic%ic_config))
         case ("eady")
            call seed_eady_ic(state, grid, cfg, ierr=local_ierr)
         case ("geostrophic_adjustment")
            call seed_geostrophic_adjustment_ic(state, grid, cfg, ierr=local_ierr)
         case ("baroclinic_jet")
            call seed_baroclinic_jet_ic(state, grid, cfg, ierr=local_ierr)
         case default
            ! "" — keep the default analytical IC unchanged.
            local_ierr = 0
         end select
         if (local_ierr /= 0) then
            ierr = local_ierr
            return
         end if
      else
         select case (trim(cfg%ocean%ic%ic_config))
         case ("eady")
            call seed_eady_ic(state, grid, cfg)
         case ("geostrophic_adjustment")
            call seed_geostrophic_adjustment_ic(state, grid, cfg)
         case ("baroclinic_jet")
            call seed_baroclinic_jet_ic(state, grid, cfg)
         case default
            ! "" — keep the default analytical IC unchanged.
         end select
      end if

      ! Z-level T/S IC overlay (capability A2).  Runs after bathymetry +
      ! uniform h_layer + wet_mask + analytical T/S are in place (it reads
      ! the seeded column depths + wet_mask) and intentionally OVERWRITES
      ! any analytical T/S.  Default off (`enable = .false.`) preserves
      ! bit-identity.  NetCDF-only: the reader lives in rdb_ocean_z_init,
      ! which only compiles with RDB_ENABLE_NETCDF=ON.
      !
      ! Under a cavity the overlay is handed `metrics%z_draft` so every
      ! layer centre's depth is measured from `z = 0` rather than from the
      ! ice base; see `seed_zinit_overlay`.
      if (cfg%ocean%zinit%enable) then
#ifndef RDB_NO_NETCDF
         ! A trimmed cavity IC moves the column top to `z_draft - eta_trim`.
         if (trim_ic) then
            if (present(ierr)) then
               call seed_zinit_overlay(state, grid, cfg, ierr=local_ierr, &
                                       z_top=state%metrics%z_draft - eta_trim)
               if (local_ierr /= 0) then
                  ierr = local_ierr
                  return
               end if
            else
               call seed_zinit_overlay(state, grid, cfg, &
                                       z_top=state%metrics%z_draft - eta_trim)
            end if
         else if (present(ierr)) then
            call seed_zinit_overlay(state, grid, cfg, ierr=local_ierr)
            if (local_ierr /= 0) then
               ierr = local_ierr
               return
            end if
         else
            call seed_zinit_overlay(state, grid, cfg)
         end if
#else
         call fail("ocean_state_seed_from_cfg: ocean_zinit requires "// &
                   "RDB_ENABLE_NETCDF=ON at build time (the NetCDF-backed "// &
                   "rdb_ocean_z_init reader is needed to load the z-level T/S file).", ierr, OCEAN_STATUS_ERR_IO)
         return
#endif
      end if

      ! The on-target `zstar_full` / `zstar` seed (FOURTH / FIFTH BRANCH
      ! above) laid inert `zstar_h_min` fillers, and every IC writer
      ! since (`&tracer_nml` per layer index, the zinit overlay at the
      ! filler's own depth) gave them a concentration that is not their
      ! donor's.  Establish invariant I1′ NOW, with the one definition (host
      ! twin — before `enter_data`), so the sponge snapshot and the budget
      ! latch see the state every later step holds.  Left to the first
      ! in-step enforcement, the pooling moves that foreign content into the
      ! live partial cell above at step 1: a real horizontal density
      ! difference at a face the mask keeps open (measured 4.6e-6 m/s after
      ! two hours on `test_ocean_zstar_full_closed_faces`' resting
      ! staircase).  Column-conservative.  Before the pseudo-salt seed, so
      ! that copies the settled S.
      if (seeded_on_eta0_target) then
         call state%multilayer%enforce_vanished_content_host(nx, ny)
      end if

      ! Pseudo-salt seed — MUST run after every write to salinity's initial
      ! condition above (analytical IC, the eady/geostrophic/baroclinic-jet
      ! overlays, and critically the z-file overlay, which OVERWRITES any
      ! analytical S) and before ocean_state_enter_data (driver-ordered).
      ! Seeding before the z-file overlay would leave D(t=0) /= 0 in exactly
      ! the configuration (z-file IC) where the diagnostic matters most.
      ! Self-gates on idx_pseudo_salt <= 0 (knob off).
      call ocean_pseudo_salt_seed(state%multilayer)

      ! Direct per-layer density init — MOM6 `COORD_CONFIG="gprime"` analogue.
      ! When `ocean_layer_rho_init` is set (any non-sentinel value), write
      ! `ms%rho_layer(:,:,k) = ocean_layer_rho_init(k)` directly, bypassing
      ! the EOS path entirely.  Intended use case: `ocean_enable_thermodynamics
      ! = .false.` + this knob ⇒ static reduced-gravity stratification, mirroring
      ! MOM6's adiabatic gprime IC.  Without this, `rho_layer` stays at its
      ! alloc-time zero when thermo is off (the dyn step's EOS call is gated
      ! on `therm_active`), and FV-LITE's pressure-stack collapses to zero
      ! ⇒ no PGF, no Sverdrup balance, no WBC dynamics — only the wind +
      ! Coriolis side of the momentum equation produces anything.  Setting
      ! the per-layer ρ here is what turns the PGF back on under the
      ! adiabatic-stratified semantics MOM6 uses for `double_gyre`.
      !
      ! Convention: `k=1` bed → `k=nz_ml` surface (Roundabout bottom-up).  Pass
      ! the heavier value first.  Count of non-sentinel entries must match
      ! `nz_ml`; otherwise we abort with a descriptive error.
      !
      ! `rho_lightest >= 0` selects the linear density-range generator
      ! (MOM6 `COORD_CONFIG="linear"`): build `nz_ml` linearly-spaced
      ! densities and reuse the same write path, so the column scales by
      ! just bumping `nz_layers`.  Mutually exclusive with the explicit
      ! `layer_rho_init` list — set only one.
      if (cfg%ocean%ic%rho_lightest >= 0.0_wp) then
         if (count(cfg%ocean%ic%layer_rho_init >= 0.0_wp) > 0) then
            call fail("ocean_ic: rho_lightest (linear density-range) and "// &
                      "layer_rho_init (explicit list) are mutually exclusive — "// &
                      "set only one.", ierr, OCEAN_STATUS_ERR_IC_SEED)
            return
         end if
         if (present(ierr)) then
            call apply_layer_rho_init(state%multilayer, &
                                      ocean_linear_layer_density(cfg%ocean%ic%rho_lightest, &
                                                                 cfg%ocean%ic%rho_range, nz_ml), &
                                      nz_ml, ierr=local_ierr)
            if (local_ierr /= 0) then
               ierr = local_ierr
               return
            end if
         else
            call apply_layer_rho_init(state%multilayer, &
                                      ocean_linear_layer_density(cfg%ocean%ic%rho_lightest, &
                                                                 cfg%ocean%ic%rho_range, nz_ml), &
                                      nz_ml)
         end if
      else
         if (present(ierr)) then
            call apply_layer_rho_init(state%multilayer, cfg%ocean%ic%layer_rho_init, nz_ml, &
                                      ierr=local_ierr)
            if (local_ierr /= 0) then
               ierr = local_ierr
               return
            end if
         else
            call apply_layer_rho_init(state%multilayer, cfg%ocean%ic%layer_rho_init, nz_ml)
         end if
      end if

      ! Populate the per-column z_ref table for VCOORD_ZSTAR_FULL.  Without
      ! this, `compute_target_h(ZSTAR_FULL)` walks an all-zero table and
      ! the first ALE remap collapses every layer to `zstar_h_min` — bug
      ! caught by the double_gyre adiabatic NK=2 setup (2026-05-26).
      ! Other vcoord types ignore z_ref so the call is safe regardless.
      call state%vcoord%build_zref_full(state%barotropic%b)
      if (present(ierr)) ierr = OCEAN_STATUS_OK
   end subroutine ocean_state_seed_from_cfg