metrics_apply_land_mask Subroutine

   subroutine metrics_apply_land_mask(this, wet_mask, grid, &
                                      periodic_x, periodic_y, north_fold, &
                                      mask_wall_velocity, &
                                      wall_west, wall_east, wall_south, wall_north)
      !! Derive the static C-grid face / corner masks from the T-cell
      !! `wet_mask` and zero the face metrics at land faces, so every
      !! transport / gradient / circulation operator that rides those
      !! metrics couples across NO land face (MOM6 pre-masks the face
      !! LENGTHS; Adcroft & Hallberg 2006).
      !!
      !! Must run at SETUP, AFTER `metrics_finalize` (the inverses
      !! `idxCu`/`idyCv` are masked here, so they must already exist) and
      !! AFTER `wet_mask` is seeded, but BEFORE `ocean_state_enter_data`
      !! (the host edit is what the GPU copyin captures).  Plain host
      !! loops — `do concurrent` before `enter_data` would round-trip the
      !! unmapped arrays through the device per loop.
      !!
      !! Halo-aware (R5a): a working copy of `wet_mask` is first filled in
      !! the ghost columns/rows by the SAME periodic wrap / north fold the
      !! metric ghosts use, so the seam u-faces (e.g. a continent that
      !! straddles `x=0≡x=1`) mask correctly.  Wall ghosts already carry
      !! the constant-extrapolated `wet_mask` from the bathymetry fill.
      !!
      !! Masks: `wet_u(i,j) = wet_T(i-1,j)*wet_T(i,j)` (Cu),
      !! `wet_v(i,j) = wet_T(i,j-1)*wet_T(i,j)` (Cv), each forced to 0 on
      !! a face of zero width (`dy_cu`/`dx_cv = 0`, see below),
      !! `wet_q(i,j) = product of the 4 T-cells around corner (i,j)` (Bu,
      !! free-slip).  Zeroed metrics (the EXACT 6 — spec §14 C3):
      !!   `dy_cu, idxCu, dxCu` at `wet_u==0`;
      !!   `dx_cv, idyCv, dyCv` at `wet_v==0`.
      !! NOT touched: `iareaT, areaT, areaCu, areaCv, iareaBu` (zeroing
      !! them would break wet-cell divergence / KE / Coriolis
      !! corner-area normalization).
      !!
      !! Bit-identity: all-wet ⇒ every `wet_*≡1` ⇒ the 6 metrics are
      !! multiplied by 1 (byte-unchanged) and the masks stay inert.
      !!
      !! Solid-wall velocity masking (`mask_wall_velocity`, opt-in): a flat
      !! all-wet channel has `wet_mask≡1` in the WALL-edge ghosts too, so
      !! `wet_v`/`wet_u` at the wall face = 1·1 = 1 (unmasked) — the wall
      !! flux is masked but the raw wall-normal velocity drifts to garbage
      !! (spurious vorticity band).  When enabled, the ghost `wm` beyond
      !! each SOLID WALL edge is zeroed (MOM6: the halo beyond a wall is
      !! land), so the derived face masks are 0 at the wall and the existing
      !! per-stage `mask_layer_velocities` clears the velocity — no bespoke
      !! velocity BC.  Only WALL edges are touched: periodic edges keep
      !! their wrapped (wet) ghosts, open/OBC edges keep the interior value.
      type(ocean_metrics_t), intent(inout) :: this
      real(wp), intent(in) :: wet_mask(:, :)
         !! T-cell wet (1) / land (0) mask, `(nx_total, ny_total)`.
      type(hgrid_t), intent(in) :: grid
      logical, intent(in) :: periodic_x, periodic_y, north_fold
         !! Boundary topology of `wet_mask`'s ghost halo (from the bc
         !! state) — selects the ghost wrap before deriving the masks.
      logical, intent(in), optional :: mask_wall_velocity
         !! Opt-in solid-wall velocity masking (default absent ⇒ .false. ⇒
         !! wall ghosts untouched ⇒ bit-identical to the legacy path).
      logical, intent(in), optional :: wall_west, wall_east, wall_south, wall_north
         !! Per-edge solid-WALL flags (an edge that is NOT periodic, NOT
         !! north-fold, NOT open/OBC).  Only consulted when
         !! `mask_wall_velocity` is .true.; each defaults to "wall" on any
         !! non-periodic / non-fold edge (the closed-default assumption).

      integer :: nx, ny, ni, nj, ng, i, j
      real(wp), allocatable :: wm(:, :)
      logical :: do_wall_mask, w_wall, e_wall, s_wall, n_wall

      nx = grid%nx_total
      ny = grid%ny_total
      ni = grid%nx_phys
      nj = grid%ny_phys
      ng = grid%nghost

      ! ---- Halo-valid working copy of wet_mask (R5a) ----
      ! The incoming `wet_mask` already carries multi-rank seam ghosts: the
      ! caller (configure_ocean_land_mask) exchanges it via ocean_halo_centre
      ! BEFORE this routine so the seam ghost columns hold the neighbour rank's
      ! real wet_T (O3 land x decomp).  Here we only add the PHYSICAL periodic /
      ! fold ghost wraps (disjoint from MPI seams).
      allocate (wm(nx, ny))
      wm = wet_mask
      if (periodic_x) call metrics_periodic_x_2d(wm, grid)
      if (periodic_y) call metrics_periodic_y_2d(wm, grid)
      if (north_fold) call fold_north_centre(wm, nx, ny, ni, nj, ng)

      ! ---- Solid-wall land fill (opt-in; MOM6 mask-in-the-update) ----
      ! Zero the ghost rows/columns beyond each SOLID WALL edge AFTER the
      ! periodic/fold wraps and BEFORE deriving the face masks, so the wall
      ! face products (wet_v/wet_u = wm(interior)*wm(ghost) = *0) vanish and
      ! `mask_layer_velocities` clears the wall-normal velocity each stage.
      ! Untouched when disabled ⇒ bit-identical.  Only WALL edges: periodic
      ! edges keep their wrapped ghosts, open/OBC edges the interior value.
      do_wall_mask = .false.
      if (present(mask_wall_velocity)) do_wall_mask = mask_wall_velocity
      if (do_wall_mask) then
         ! Default: any non-periodic, non-fold edge is a wall (closed default);
         ! callers thread the true per-edge WALL/OPEN flags to override.
         w_wall = .not. periodic_x
         e_wall = .not. periodic_x
         s_wall = .not. periodic_y
         n_wall = (.not. periodic_y) .and. (.not. north_fold)
         if (present(wall_west)) w_wall = wall_west
         if (present(wall_east)) e_wall = wall_east
         if (present(wall_south)) s_wall = wall_south
         if (present(wall_north)) n_wall = wall_north
         ! Ghost cells: west i=1..ng, east i=ng+ni+1..nx,
         !              south j=1..ng, north j=ng+nj+1..ny.
         if (w_wall) then
            do j = 1, ny
               do i = 1, ng
                  wm(i, j) = 0.0_wp
               end do
            end do
         end if
         if (e_wall) then
            do j = 1, ny
               do i = ng + ni + 1, nx
                  wm(i, j) = 0.0_wp
               end do
            end do
         end if
         if (s_wall) then
            do j = 1, ng
               do i = 1, nx
                  wm(i, j) = 0.0_wp
               end do
            end do
         end if
         if (n_wall) then
            do j = ng + nj + 1, ny
               do i = 1, nx
                  wm(i, j) = 0.0_wp
               end do
            end do
         end if
      end if

      ! ---- Store the halo-valid T-cell mask (consumed by PPM mirror-h) ----
      this%wet_T = wm

      ! ---- Derive face / corner masks (plain host loops) ----
      ! wet_u(i,j): u-face i = west face of T-cell (i,j); pairs (i-1,i).
      ! Array outer ring (i=1, nx+1) is outside the ghost band => always
      ! land, decomposition-invariant; physical/seam interface faces at
      ! nghost+1 are bathymetry-masked (wet_T product below), not edge-
      ! position-masked, so no has_west/has_east gate is needed here
      ! (O0 verified: seam face wet_T product = 1*1 = 1, mask stays open).
      !
      ! A face of ZERO width is a wall whatever its neighbours are.  The
      ! tripolar cap's node-aligned pole columns (`tripolar_node_latlon`
      ! places every cap node of a pole column on the pole) have
      ! `dy_cu = 0` between two WET cells: no transport, zero `areaCu`
      ! (no kinetic energy), zero circulation weight at the pole corners
      ! (`iareaBu = 0`).  Left open, its velocity is still a prognostic:
      ! the free-surface gradient across the pole drives it, the
      ! barotropic fast loop's Coriolis couples it to the neighbouring
      ! v faces with a plain 1/4 weight, and under `pred_corr` its time
      ! mean `u_av` is never written (the renormaliser's `u_cor` skips a
      ! face with `sum h*dy_cu = 0`), so `set_cor_ref_velocity` subtracted
      ! a frozen reference while the fast loop integrated the live one.
      ! Measured on the compatibility matrix's tripolar domain (wind +
      ! cooling, 30 d): `En` 0.66 m2/s2 and saturating (pred_corr) vs
      ! 3.2e-3 with the face closed; `ssp_rk2` 4.7e-3 -> 3.2e-3.  Closing
      ! it here makes it a coast to every consumer at once (the six face
      ! metrics below, `mask_layer_velocities`, the `u_av` seed mask, the
      ! barotropic `mask_bt_rem`).  `<= 0` rather than `== 0`: a metric
      ! length is never negative, so this only names the exact zero.
      ! Bit-identical wherever every zero-width face already touches land:
      ! every generator but a node-aligned tripolar cap has none, and the
      ! OM4 1-degree supergrid's 231 (its pole columns) all sit next to a
      ! land cell of `bathy_om1deg.nc` (global_1deg: day-10 En unchanged).
      do j = 1, ny
         this%wet_u(1, j) = 0.0_wp      ! west outer wall (no T-cell i=0)
         do i = 2, nx
            this%wet_u(i, j) = wm(i - 1, j)*wm(i, j)
            if (this%dy_cu(i, j) <= 0.0_wp) this%wet_u(i, j) = 0.0_wp
         end do
         this%wet_u(nx + 1, j) = 0.0_wp  ! east outer wall (no T-cell nx+1)
      end do
      ! wet_v(i,j): v-face j = south face of T-cell (i,j); pairs (j-1,j).
      do i = 1, nx
         this%wet_v(i, 1) = 0.0_wp
         do j = 2, ny
            this%wet_v(i, j) = wm(i, j - 1)*wm(i, j)
            if (this%dx_cv(i, j) <= 0.0_wp) this%wet_v(i, j) = 0.0_wp
         end do
         this%wet_v(i, ny + 1) = 0.0_wp
      end do
      ! wet_q(i,j): SW corner of T-cell (i,j); product of the 4 T-cells
      ! (i-1,j-1),(i,j-1),(i-1,j),(i,j).  Outer ring (i=1/nx+1, j=1/ny+1)
      ! has a missing T-neighbour ⇒ land (matches the domain-wall corner).
      do j = 1, ny + 1
         do i = 1, nx + 1
            if (i >= 2 .and. i <= nx .and. j >= 2 .and. j <= ny) then
               this%wet_q(i, j) = wm(i - 1, j - 1)*wm(i, j - 1)* &
                                  wm(i - 1, j)*wm(i, j)
            else
               this%wet_q(i, j) = 0.0_wp
            end if
         end do
      end do

      ! ---- Zero the 6 face metrics at land faces (spec §14 C3) ----
      ! u-faces (Cu): dy_cu, idxCu, dxCu.
      do j = 1, ny
         do i = 1, nx + 1
            this%dy_cu(i, j) = this%dy_cu(i, j)*this%wet_u(i, j)
            this%idxCu(i, j) = this%idxCu(i, j)*this%wet_u(i, j)
            this%dxCu(i, j) = this%dxCu(i, j)*this%wet_u(i, j)
         end do
      end do
      ! v-faces (Cv): dx_cv, idyCv, dyCv.
      do j = 1, ny + 1
         do i = 1, nx
            this%dx_cv(i, j) = this%dx_cv(i, j)*this%wet_v(i, j)
            this%idyCv(i, j) = this%idyCv(i, j)*this%wet_v(i, j)
            this%dyCv(i, j) = this%dyCv(i, j)*this%wet_v(i, j)
         end do
      end do

      ! Re-sync the BT transport widths with the freshly masked slow-path
      ! widths (this runs AFTER metrics_finalize, which set them equal).
      this%dy_cu_bt = this%dy_cu
      this%dx_cv_bt = this%dx_cv

      deallocate (wm)
   end subroutine metrics_apply_land_mask