configure_ocean_closed_faces Subroutine

   subroutine configure_ocean_closed_faces(cfg, ocean_state, grid, compute_rank, ierr)
      !! Build the static partial-step z-level FACE-CLOSURE mask
      !! (`&vcoord_nml zfixed_closed_faces`; Adcroft, Hill & Marshall
      !! 1997; Losch 2008 §2.1 for the ice-shelf cavity).
      !!
      !! Under `vcoord_type = "z_fixed"` a layer whose nominal
      !! geopotential range lies inside the bed — or inside the ice draft
      !! — carries an inert FILLER of thickness `zstar_h_min`
      !! (`<= H_VANISHED`).  A velocity face at which layer `k` is a
      !! filler on EITHER side stays OPEN today, and the FV pressure
      !! gradient evaluated across that staircase step drives
      !! `a = |ρ′|·g·Δz_step/(ρ₀·dx)` out of a resting stratified state —
      !! independent of the filler thickness, so no `h`-gate can reach
      !! it.  A z-LEVEL model treats such a face as a WALL for that
      !! layer: no normal velocity, no mass or tracer flux, free-slip.
      !!
      !! This fills `metrics%open_u/open_v` with that wall, ONCE, from
      !! the coordinate's target at `η = 0` (`ocean_vcoord_eta0_target`)
      !! — the same kernel the ALE regrid and the IC seed use, so there
      !! is no second definition of "live".  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), so the live/filler
      !! pattern does not move to first order in `η/h_partial`.
      !!
      !! **`zstar_full` too.**  `build_zref_full` lays a z-level fine
      !! zone (`h_surf` layers from the surface) over a terrain-following
      !! coarse zone; a column shallower than the fine zone ends in a
      !! partial cell and every layer below it is a `zstar_h_min` filler —
      !! the same staircase as `z_fixed`'s bed, with the same
      !! thickness-independent FV-PGF defect across it.  Its target puts
      !! `η ≥ 0` in the surface layer (pattern EXACTLY static) and clips
      !! `η < 0` from the bed (a bed-most live layer thinner than `|η|`
      !! flips — fewer columns than `z_fixed` flips at the same `|η|` on
      !! the 1-degree Southern Ocean).  So the mask is built from the
      !! `ZSTAR_FULL` target exactly as it is from the `z_fixed` one, and
      !! every consumer is reused unchanged.  It closes FILLER faces
      !! only; the terrain-following coarse zone's open faces keep the
      !! sigma PGF error.
      !!
      !! **`zstar` (MOM6 z*).**  `ocean_vcoord_zstar_target` dilates the
      !! `z_fixed` nominal profile by the column's free-surface stretching
      !! and decides every layer's liveness at `η = 0`, so its pattern is
      !! EXACTLY static for `η` of either sign (MOM6 `build_zstar_column`:
      !! the dilation keeps the ratios) and its `η = 0` target is the
      !! `z_fixed` one bit for bit — the same mask, every consumer reused.
      !!
      !! It also seeds the barotropic face widths `dy_cu_bt`/`dx_cv_bt`
      !! with the OPEN-depth fraction of the face, so the barotropic
      !! solve is not blind to the closed layers.  That seed is refreshed
      !! from the LIVE `h` every outer step by `ocean_porous_refresh`;
      !! this is only the `η = 0` value the first stage reads.
      !!
      !! **Ordering.** MUST run AFTER `configure_ocean_cavity` (which
      !! fills `vcoord%z_top`), after `configure_ocean_bt_split` (which
      !! lays `bt_H_ref`) and after `configure_ocean_land_mask` and the
      !! periodic-wrap / halo pass (the target is built from GHOST-FILLED
      !! `bt_H_ref` and `z_top`, which is what makes the mask correct at
      !! a periodic seam — the physical seam face is an interior index).
      !! And BEFORE `ocean_state_enter_data`: the host fill is what the
      !! `copyin` captures, and `metrics_closed_faces_alloc` reallocs.
      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
      integer, intent(out), optional :: ierr
         !! Non-zero on a configuration conflict when present; absent
         !! behaves as today (`error stop`).

      integer :: nx, ny, nz, i, j, k
      integer :: n_closed_u, n_closed_v, n_open_u, n_open_v, n_ledge
      real(wp) :: h_face, sum_all, sum_open
      real(wp), allocatable :: tgt(:, :, :)
      character(len=16) :: vcoord_label

      if (present(ierr)) ierr = OCEAN_STATUS_OK
      if (.not. cfg%zfixed_closed_faces) then
         call refuse_open_zfixed_staircase(ocean_state, grid, ierr)
         return
      end if

      if (.not. ocean_state%multilayer%is_init) then
         call fail("&vcoord_nml zfixed_closed_faces requires the ocean "// &
                   "multilayer path (the mask is per-layer)", ierr, OCEAN_STATUS_ERR_SETUP)
         return
      end if
      if (ocean_state%vcoord%coord_type /= VCOORD_Z_FIXED .and. &
          ocean_state%vcoord%coord_type /= VCOORD_ZSTAR .and. &
          ocean_state%vcoord%coord_type /= VCOORD_ZSTAR_FULL) then
         call fail("&vcoord_nml zfixed_closed_faces is only defined for "// &
                   "vcoord_type='z_fixed', 'zstar' and 'zstar_full': the "// &
                   "live/filler staircase it closes is made by a GEOMETRIC "// &
                   "coordinate that vanishes bed-side layers at fixed "// &
                   "reference depths. On sigma / zstar_sigma every layer is live on "// &
                   "every wet face, and rho / hycom vanish layers by DENSITY, "// &
                   "so their pattern is not static", &
                   ierr, OCEAN_STATUS_ERR_SETUP)
         return
      end if
      if ((ocean_state%vcoord%coord_type == VCOORD_Z_FIXED .or. &
           ocean_state%vcoord%coord_type == VCOORD_ZSTAR) .and. &
          ocean_state%vcoord%z_fixed_h_ref <= 0.0_wp) then
         call fail("&vcoord_nml zfixed_closed_faces needs a resolved "// &
                   "z_fixed_h_ref (set &ocean_topo_nml max_depth): without "// &
                   "it the z_fixed / zstar target degenerates to uniform sigma, "// &
                   "there are no fillers, and the mask would close nothing", &
                   ierr, OCEAN_STATUS_ERR_SETUP)
         return
      end if
      if (ocean_state%vcoord%coord_type == VCOORD_ZSTAR_FULL .and. &
          ocean_state%vcoord%zstar_h_surf_target <= 0.0_wp) then
         call fail("&vcoord_nml zfixed_closed_faces under vcoord_type="// &
                   "'zstar_full' needs zstar_h_surf_target > 0: without it "// &
                   "build_zref_full lays uniform sigma, there are no fillers, "// &
                   "and the mask would close nothing", &
                   ierr, OCEAN_STATUS_ERR_SETUP)
         return
      end if
      if (ocean_state%dyn%bt_work%is_init) then
         if (maxval(ocean_state%dyn%bt_work%bt_H_ref) <= 0.0_wp) then
            call fail("&vcoord_nml zfixed_closed_faces: bt_H_ref is empty — "// &
                      "the barotropic datum must be built (&ocean_bt_nml "// &
                      "n_inner >= 1) before the static face mask can be laid", &
                      ierr, OCEAN_STATUS_ERR_SETUP)
            return
         end if
      end if
      if (cfg%ocean%wetdry%enable) then
         call fail("&vcoord_nml zfixed_closed_faces is incompatible with "// &
                   "&ocean_wetdry_nml enable: wet/dry moves the live/filler "// &
                   "pattern under the running state and the mask is static", &
                   ierr, OCEAN_STATUS_ERR_SETUP)
         return
      end if
      if (cfg%ocean%bt%bt_halo > 0) then
         call fail("&vcoord_nml zfixed_closed_faces is incompatible with "// &
                   "&ocean_bt_nml bt_halo > 0: the wide-halo barotropic "// &
                   "clone carries its own metrics and no face mask, so the "// &
                   "wide BT loop would transport through closed faces "// &
                   "(the same argument &ocean_porous_nml already makes)", &
                   ierr, OCEAN_STATUS_ERR_SETUP)
         return
      end if
      ! GM composes: `gm_column_x/y` build the streamfunction on each face's
      ! OPEN column (`metrics%open_u/open_v` x live on both sides) and the
      ! slopes slot masks slope / N^2 to it, so `uhD`/`vhD` are zero on
      ! every closed face-layer by construction (`test_ocean_gm_zfixed`).
      ! Redi composes: each face pairs only its OPEN WINDOW (the contiguous
      ! open + live-both-sides run), so the neutral-surface sweep, the PPM
      ! reconstruction and the flux scatter never touch a closed face-layer
      ! or a filler (`test_ocean_redi_zfixed`).
      ! Fox-Kemper MLE composes: the ML walk skips fillers and each face
      ! builds its overturning on its OPEN column (masked face thickness,
      ! `H_vel` clamped to it), so `uhml`/`vhml` — folded into
      ! mass_flux_*_layer after the per-layer mask — are zero on every
      ! closed face-layer and filler and still sum to zero per face
      ! (`test_ocean_mle_zfixed`).
      ! The velocity-form BIHARMONIC (scalar `nu_4` and the flow-aware
      ! `smag_ah` / `leith_biharm` `nu4_face_*`) composes: both chained
      ! Laplacians gate every difference by `open(a)*open(b)` (free-slip,
      ! `hvisc_biharm_lap_closed`).  The stress-tensor assembly (and
      ! `kh_aniso`, which only it reads) does not: its T-cell tension and
      ! corner shear are built from the face velocities with the 2-D
      ! `wet_*` masks only.
      if (cfg%ocean%hvisc%stress_tensor) then
         call fail("&vcoord_nml zfixed_closed_faces does not yet compose "// &
                   "with &ocean_hvisc_nml stress_tensor (nor kh_aniso, which "// &
                   "only the stress-tensor path reads): its T-cell tension and "// &
                   "corner shear are masked by the 2-D wet_* fields only, so a "// &
                   "closed face-layer's zeroed velocity enters the strain as a "// &
                   "Dirichlet value (no-slip at every staircase step).  The "// &
                   "velocity-Laplacian harmonic kernels and both biharmonic "// &
                   "paths (nu_4, smag_ah, leith_biharm) carry the free-slip "// &
                   "closure; use those", &
                   ierr, OCEAN_STATUS_ERR_SETUP)
         return
      end if
      ! ---- Barotropic paths still on FULL-column weights ----
      !
      ! `derive_bt_from_layers`, `face_depth_mean_*`, `set_cor_ref_velocity`,
      ! the `apply_bt_correction` fold, `compute_h_face_upstream`
      ! (`upstream_h_face`), `compute_bt_rem` (`substep_drag`) and
      ! `compute_bt_rem_wave_drag` (`wave_drag`) weight by `h_face·open`
      ! under this knob.  The bc-PGF correction does NOT: it sums
      ! `0.5·(h_L + h_R)` over EVERY layer.  Its answer is not wrong by a
      ! round-off — the closed layers' thickness is O(h_nominal) at a
      ! staircase face — so it is refused until it is ported, rather than
      ! left to run on the wrong column.
      if (cfg%ocean%bt%correction_bc_pgf) then
         call fail("&vcoord_nml zfixed_closed_faces does not yet compose "// &
                   "with &ocean_bt_nml correction_bc_pgf: compute_pbce, "// &
                   "compute_gtot_faces and the bc-PGF du_bc block in "// &
                   "apply_bt_correction all weight by the FULL column "// &
                   "(0.5*(h_L+h_R) over every layer), so the correction's "// &
                   "depth-mean-zero identity holds on the full column, not "// &
                   "on the OPEN column ubt is the mean of — it would leak a "// &
                   "barotropic increment into the open layers and push one "// &
                   "into the closed ones.  Set correction_bc_pgf=.false.", &
                   ierr, OCEAN_STATUS_ERR_SETUP)
         return
      end if

      nx = grid%nx_total
      ny = grid%ny_total
      nz = ocean_state%multilayer%nz_ml
      vcoord_label = "z_fixed"
      if (ocean_state%vcoord%coord_type == VCOORD_ZSTAR_FULL) vcoord_label = "zstar_full"
      if (ocean_state%vcoord%coord_type == VCOORD_ZSTAR) vcoord_label = "zstar"

      call metrics_closed_faces_alloc(ocean_state%metrics, grid, nz)

      ! The eta = 0 target of the running coordinate — `z_fixed`, `zstar`
      ! or `zstar_full` — through the SAME kernel the ALE regrid dispatches
      ! to, so "live" has one definition.  Under `zstar_full` it walks the
      ! per-column `z_ref` table the engine rebuilt from the
      ! halo-exchanged bathymetry just before the static-geometry pass.
      allocate (tgt(nx, ny, nz), source=0.0_wp)
      call ocean_vcoord_eta0_target(ocean_state%vcoord, tgt, &
                                    ocean_state%dyn%bt_work%bt_H_ref, nx, ny, nz)
      call ocean_vcoord_closed_face_masks(ocean_state%metrics%open_u, &
                                          ocean_state%metrics%open_v, &
                                          tgt, nx, ny, nz, H_VANISHED)
      ! The builder cannot evaluate the OUTERMOST face of the array (it
      ! needs a cell beyond it) and leaves it open.  On a tile that face
      ! is a SEAM ghost, three cells from the owned region, and on a
      ! periodic edge it is the wrap of an interior face — so without
      ! this exchange the ghost band of the mask is decomposition-
      ! dependent.  The harmonic kernels never read that deep; the
      ! biharmonic's chained stencil does (its ghost-band tendency at
      ! depth 1 reads the mask at depth 3), which broke 1x4 bit-identity
      ! in `test_ocean_decomp_bitid_mpi` (`file_readers`).  A plain face
      ! exchange (no sign flip: a 0/1 mask, not a vector) makes every
      ! ghost face carry its owner's value.  Single-rank non-periodic:
      ! a no-op.
      if (ocean_halo_is_init()) then
         call ocean_halo_face_x(ocean_state%metrics%open_u, nz, device_resident=.false.)
         call ocean_halo_face_y(ocean_state%metrics%open_v, nz, device_resident=.false.)
      end if

      ! Seed the BAROTROPIC face widths with the eta = 0 open-depth
      ! fraction.  `ocean_porous_refresh` recomputes them from the live
      ! `h` every outer step; this is only what stage 1 of step 1 reads.
      ! The width and the depth are the same number here: the BT
      ! transport is `ubt * FA * dy_cu_bt` with `FA = sum_k h_face`, so
      ! `dy_cu_bt = dy_cu * (sum_open h_face)/(sum_k h_face)` makes that
      ! product identically `ubt * (sum_open h_face) * dy_cu`.
      do j = 1, ny
         do i = 2, nx
            sum_all = 0.0_wp
            sum_open = 0.0_wp
            do k = 1, nz
               h_face = 0.5_wp*(tgt(i - 1, j, k) + tgt(i, j, k))
               sum_all = sum_all + h_face
               if (ocean_state%metrics%open_u(i, j, k) > 0.0_wp) then
                  sum_open = sum_open + h_face
               end if
            end do
            if (sum_all > 0.0_wp) then
               ocean_state%metrics%dy_cu_bt(i, j) = &
                  ocean_state%metrics%dy_cu(i, j)*(sum_open/sum_all)
            end if
         end do
      end do
      do j = 2, ny
         do i = 1, nx
            sum_all = 0.0_wp
            sum_open = 0.0_wp
            do k = 1, nz
               h_face = 0.5_wp*(tgt(i, j - 1, k) + tgt(i, j, k))
               sum_all = sum_all + h_face
               if (ocean_state%metrics%open_v(i, j, k) > 0.0_wp) then
                  sum_open = sum_open + h_face
               end if
            end do
            if (sum_all > 0.0_wp) then
               ocean_state%metrics%dx_cv_bt(i, j) = &
                  ocean_state%metrics%dx_cv(i, j)*(sum_open/sum_all)
            end if
         end do
      end do

      ! Cheap one-shot census for the configure line: how much of the
      ! array the mask actually closes, and whether it isolated any
      ! water.  Host-side, once, O(nx*ny*nz) — not a per-step diagnostic.
      n_closed_u = 0
      n_open_u = 0
      do k = 1, nz
         do j = 1, ny
            do i = 2, nx
               if (ocean_state%metrics%open_u(i, j, k) > 0.0_wp) then
                  n_open_u = n_open_u + 1
               else
                  n_closed_u = n_closed_u + 1
               end if
            end do
         end do
      end do
      n_closed_v = 0
      n_open_v = 0
      do k = 1, nz
         do j = 2, ny
            do i = 1, nx
               if (ocean_state%metrics%open_v(i, j, k) > 0.0_wp) then
                  n_open_v = n_open_v + 1
               else
                  n_closed_v = n_closed_v + 1
               end if
            end do
         end do
      end do
      n_ledge = ocean_vcoord_count_ledges(ocean_state%metrics%open_u, &
                                          ocean_state%metrics%open_v, &
                                          tgt, nx, ny, nz, H_VANISHED)

      deallocate (tgt)

      ! Flip the switches LAST — every consumer branches on them, and the
      ! mask has to be in place before any of them can read a closed face.
      ocean_state%metrics%use_closed_faces = .true.
      ocean_state%vcoord%zfixed_closed_faces = .true.
      ocean_state%vdiff%zlevel_faces = .true.

      if (compute_rank == 0) then
         call logger%info(trim(vcoord_label)//" closed faces: ON (partial steps, "// &
                          "Adcroft/Hill/Marshall 1997) — u closed "// &
                          to_string(n_closed_u)//"/"// &
                          to_string(n_closed_u + n_open_u)//", v closed "// &
                          to_string(n_closed_v)//"/"// &
                          to_string(n_closed_v + n_open_v)// &
                          ", isolated ledge cells "//to_string(n_ledge))
         if (n_ledge > 0) then
            call logger%warning(trim(vcoord_label)//" closed faces: "//to_string(n_ledge)// &
                                " LIVE cells have all four own-layer faces "// &
                                "closed — the mask has isolated water.  They "// &
                                "are inert (no flux in or out, velocity zeroed "// &
                                "every stage) but a one-cell spike in the bed "// &
                                "or the draft is worth looking at.")
         end if
      end if
      if (present(ierr)) ierr = OCEAN_STATUS_OK
   end subroutine configure_ocean_closed_faces