configure_ocean_land_mask Subroutine

   subroutine configure_ocean_land_mask(cfg, ocean_state, grid, compute_rank, warm_restart)
      !! Derive the static C-grid land masks from the seeded T-cell
      !! `wet_mask` and zero the 6 face metrics at land faces
      !! (`metrics_apply_land_mask`).  Run AFTER `configure_ocean_metrics`
      !! (the metrics + inverses must exist) AND `configure_ocean_bc` (the
      !! periodic / fold flags drive the halo-aware mask derivation), but
      !! BEFORE `ocean_state_enter_data` (the host edit is what the GPU
      !! copyin captures).
      !!
      !! For a domain with no land (`wet_mask≡1`, the flat-bottom /
      !! analytical path) every mask is 1.0 and the metric multiply is a
      !! literal no-op ⇒ bit-identical to a no-mask build — except that a
      !! face of ZERO width (`dy_cu`/`dx_cv = 0`, a node-aligned tripolar
      !! cap's pole columns) is closed even between wet cells.  This call
      !! is the only place that closure happens, so every geometry path
      !! must run it.
      type(config_t), intent(in) :: cfg
      type(ocean_state_t), intent(inout) :: ocean_state
      type(hgrid_t), intent(in) :: grid
      integer, intent(in) :: compute_rank
      logical, intent(in), optional :: warm_restart
         !! `.true.` when the prognostic state was just read from a
         !! checkpoint: the land-state seed below is then SKIPPED (the
         !! masks and metrics are still derived).  Absent ⇒ cold start.

      logical :: seed_land

      ! Multi-rank seam ghost fill of wet_mask (O3 land x decomp): a subdomain
      ! seam that bisects a continent needs the NEIGHBOUR rank's wet_T in the
      ! seam ghost columns so the C-grid face/corner masks derived inside
      ! metrics_apply_land_mask close the seam faces correctly.  Done HERE (the
      ! caller) rather than inside metrics_apply_land_mask because importing the
      ! comm-layer rdb_ocean_halo into the low-level rdb_ocean_metrics leaf
      ! creates an NVFORTRAN USE cycle (config -> lateral_mix -> metrics).
      ! ocean_halo_init has already run (driver orders it before this call).
      ! Gated on ocean_halo_is_DECOMPOSED (px>1 or py>1) — a genuine
      ! cross-rank seam is the ONLY case wet_mask needs a message: the
      ! periodic ghost wrap is done locally inside metrics_apply_land_mask,
      ! and a single-rank non-periodic run's ghosts are domain-boundary
      ! walls.  Gating on the seam (not periodicity / not merely is_init)
      ! also keeps the MPI runtime untouched in non-MPI ocean unit tests
      ! that init the halo via this setup path but never call MPI_Init
      ! (those are px=py=1 ⇒ skip).  Runs host-side before enter_data; a
      ! GPU+MPI init-time host-halo variant is an O3+ follow-up.
      if (ocean_halo_is_decomposed()) then
         call ocean_halo_centre(ocean_state%multilayer%wet_mask, device_resident=.false.)
      end if
      ! Distributed tripolar fold (px > 1, plan site S2):
      ! `metrics_apply_land_mask` folds its working copy with the LOCAL
      ! kernel, which is exact only on a tile holding the whole fold row.
      ! Fold the STORED wet mask here instead, through the owner-routed
      ! exchange (collective over the north rank row, host mode), and tell
      ! the land mask not to fold again.  Folding the stored mask, not a
      ! copy, also closes the seed-time gap for it: the seed derives
      ! `wet_mask` elementwise from the water column, whose north ghosts
      ! the single-rank seed folded first (`seed_wrap_static_2d`) and the
      ! px > 1 seed could not, so the stored ghosts now equal the serial
      ! run's bit for bit (an elementwise map commutes with the fold copy).
      if (ocean_fold_is_distributed()) then
         call ocean_fold_wrap_eta_2d(grid, ocean_state%bc, ocean_state%multilayer%wet_mask, &
                                     device_resident=.false.)
      end if

      ! Solid-wall velocity masking (&ocean_bc_nml mask_wall_velocity):
      ! thread the per-edge WALL flags so only genuine solid walls get their
      ! ghost wet_T zeroed (periodic + open/OBC edges keep their ghosts).
      !
      ! The `has_*` conjunct is load-bearing under MPI.  `bc_type` is the
      ! GLOBAL edge tag, identical on every rank; `has_*` is the per-rank
      ! "I actually sit on that domain edge" flag (set from decomp).
      ! Without the AND, every rank zeroed its own west/east/south/north
      ! ghost band, so an interior MPI SEAM was masked as land: wet_u at
      ! the seam face went to 0 and the six masked face metrics with it.
      ! Measured on double_gyre_mom6 (44x40, 10 days, wall basin): px=1
      ! gave En 3.036E-03 / MaxCFL 0.01000, px=2 gave 2.998E-03 / 0.01197
      ! — and MaxCFL is a max reduction, which is exactly order-invariant,
      ! so that spread could not have been reduction roundoff.  With the
      ! conjunct, px=1 and px=2 agree.
      ! The periodic flags select a LOCAL wrap of the working wet mask's
      ! ghosts, which is right only on an axis this rank holds whole: on a
      ! decomposed periodic axis the wrap ghosts are an MPI seam the
      ! exchange above already filled, and a local wrap overwrote them with
      ! the tile's own opposite edge (every rank of a px >= 3 periodic run
      ! with land near a seam then masked the wrong ghost faces; the
      ! decomposed-vs-serial bit-identity test caught it on 4x1).  The wall
      ! flags below carry the edge tags, so nothing else reads these two.
      call metrics_apply_land_mask(ocean_state%metrics, &
                                   ocean_state%multilayer%wet_mask, grid, &
                                   ocean_state%bc%periodic_x .and. .not. ocean_halo_is_decomposed_x(), &
                                   ocean_state%bc%periodic_y .and. .not. ocean_halo_is_decomposed_y(), &
                                   ocean_state%bc%north_fold .and. .not. ocean_fold_is_distributed(), &
                                   mask_wall_velocity=cfg%ocean%bc%mask_wall_velocity, &
                                   wall_west=(ocean_bc_outer_face_tag(ocean_state%bc%west%bc_type) == OBC_WALL &
                                              .and. ocean_state%bc%has_west), &
                                   wall_east=(ocean_bc_outer_face_tag(ocean_state%bc%east%bc_type) == OBC_WALL &
                                              .and. ocean_state%bc%has_east), &
                                   wall_south=(ocean_bc_outer_face_tag(ocean_state%bc%south%bc_type) == OBC_WALL &
                                               .and. ocean_state%bc%has_south), &
                                   wall_north=(ocean_bc_outer_face_tag(ocean_state%bc%north%bc_type) == OBC_WALL &
                                               .and. ocean_state%bc%has_north))

      ! Hold land T-cells at finite reference values so masked arithmetic
      ! never multiplies 0 by a NaN (0*NaN = NaN).  Land h_layer is floored
      ! to H_VANISHED (never 0), tracers held at the IC reference, layer +
      ! face velocities zeroed on land.
      !
      ! Cold start only.  The seed is the land state's INITIAL value, not an
      ! invariant the step maintains: the ALE remap regrids a land column
      ! like any other (a `z_fixed` land column leaves step 1 as `nz-1`
      ! `zstar_h_min` fillers over a bed cell holding the rest of its
      ! `nz·H_VANISHED`), and the checkpoint saves that state verbatim.
      ! Re-seeding after a restart read rewrote every land column back to
      ! uniform `H_VANISHED` -- column total unchanged, layout not -- so the
      ! resumed run was not the run that wrote the file (1/4-degree Southern
      ! Ocean, 2026-10-01: 56 046 land columns x 50 layers per rank).
      seed_land = .true.
      if (present(warm_restart)) seed_land = .not. warm_restart
      if (seed_land) call ocean_state_seed_land_cells(ocean_state, grid)

      if (compute_rank == 0) then
         associate (unused => cfg%ocean%grid%grid_config)
         end associate
      end if
   end subroutine configure_ocean_land_mask