configure_ocean_cavity Subroutine

   subroutine configure_ocean_cavity(cfg, ocean_state, grid, compute_rank, ierr)
      !! Build the static ice-shelf cavity LOAD field
      !! `metrics%p_ice_ref = (rho_ref*GRAVITY)*z_draft` (Pa), ASSEMBLE it
      !! into the top-of-column pressure `multilayer_state_t%p_top`, and
      !! assert the counted-once datum invariant.
      !!
      !! The GEOMETRY (`z_draft`, `cover_frac`) is filled much earlier, in
      !! `ocean_state_seed_from_cfg`, because the wet mask and the layer
      !! split are seeded from `b − z_draft`.  What is left for configure
      !! is the part that needs the PGF's reference density, which
      !! `configure_ocean_pgf` / `configure_ocean_reference_density` only
      !! settle later — hence this runs after them and before
      !! `ocean_state_enter_data`, like every other static field the
      !! device map has to capture.
      !!
      !! `rho_ref*GRAVITY` is formed as ONE product, the same one the
      !! FV_MOM6 Pass-1 surface BC forms, so that
      !! `pa(nz+1) = rho_ref*g*(−z_draft) + p_ice_ref` cancels to bit-zero
      !! at rest (exactly when the toolchain rounds the product before the
      !! add; under FMA contraction, to the rounding of `p_ice_ref`).
      !!
      !! ### The partition (P5.2)
      !!
      !! ```
      !! ms%p_top = metrics%p_ice_ref  +  sf%p_surf
      !!            (static ice load)     (atmospheric / anomaly load)
      !! ```
      !!
      !! `p_ice_ref` goes HERE and to the datum (`bt_H_ref = b − z_draft`)
      !! and NOWHERE else — in particular it is never added into
      !! `sf%p_surf`, because `eta_ib = −p_surf/(ρ₀ g_bt)` is built from
      !! the assembled total and the datum already carries exactly this
      !! much.  Only the load ANOMALY reaches the `eta_forcing` seam, and
      !! for the Boussinesq-isostatic default that anomaly IS `sf%p_surf`
      !! (a cavity with no atmospheric load sends the seam nothing at
      !! all).  See `src/core/ocean/README.md`'s `p_top` / `eta_forcing`
      !! seam contracts.
      !!
      !! This seed is the FINAL value for a cavity without the psurf seam
      !! (the draft is static, so there is nothing to refresh); with psurf
      !! enabled, `ocean_dyn_step_split` rebuilds the same sum once per
      !! outer step from the live `sf%p_surf`.  Host-side plain `do`
      !! loops, before `ocean_state_enter_data` maps the result.
      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 cavity configuration conflict when present;
         !! absent behaves as today (`error stop`).

      integer :: nx, ny, i, j
      real(wp) :: resid, datum_tol

      if (.not. ocean_state%metrics%use_cavity) then
         if (present(ierr)) ierr = OCEAN_STATUS_OK
         return
      end if

      nx = size(ocean_state%metrics%z_draft, 1)
      ny = size(ocean_state%metrics%z_draft, 2)

      call cavity_fill_p_ice_ref(ocean_state%metrics%p_ice_ref, &
                                 ocean_state%metrics%z_draft, &
                                 ocean_state%pressure_force%rho_ref*GRAVITY, nx, ny)

      ! P5.2 — THE LOAD MUST HAVE A CONSUMER WHEN IT HAS A GRADIENT.
      !
      ! `&ocean_pgf_nml p_top_in_bc` is the only route by which a cavity
      ! load reaches the pressure stack, and a cavity whose draft VARIES
      ! is not in hydrostatic balance without it: `pa(nz+1)` then carries
      ! the uncancelled `-rho_ref*g*z_draft`, so the stack sits ~5e6 Pa
      ! off its anomaly scale (every `h_neglect` face-divisor leak grows
      ! by the same factor), the UNSPLIT driver — which has no
      ! depth-mean replacement — feels a raw `g*grad(z_draft)` ~ 0.1 m/s^2,
      ! and a non-uniform barotropic-correction weight turns the
      ! uncancelled depth-uniform force into a real per-layer shear.
      ! Refused, not auto-enabled: the namelist should say what the run does.
      !
      ! EXEMPTION, and it is a theorem rather than a courtesy: a draft
      ! that is UNIFORM over the whole array has a load with no gradient,
      ! and a gradient-free `p_top` is bit-identically inert in the top BC
      ! (`test_ocean_pgf_p_top_bc::uniform_p_top_bit_identical`).  Such a
      ! run — the flat-lid datum-equivalence case — is refused by nothing.
      !
      ! Tested on the FILLED array, not on the namelist shape: land
      ! exclusion can zero `z_draft` under a formula that looks uniform.
      ! `validate_config` carries the same refusal on the config-level
      ! predicate so the user meets it before any state is built.
      if (.not. ocean_state%pressure_force%p_top_in_bc) then
         if (maxval(ocean_state%metrics%z_draft) /= &
             minval(ocean_state%metrics%z_draft)) then
            call fail("&ocean_cavity_dyn_nml enable=.true. with a NON-UNIFORM "// &
                      "draft requires &ocean_pgf_nml p_top_in_bc=.true.  The "// &
                      "isostatic load rho_ref*g*z_draft would otherwise never "// &
                      "reach the FV_MOM6 pa(nz+1) surface boundary condition, "// &
                      "leaving the pressure stack ~"// &
                      to_string(maxval(ocean_state%metrics%p_ice_ref))// &
                      " Pa off its anomaly scale and the column out of "// &
                      "hydrostatic balance.  (A UNIFORM draft is exempt: a "// &
                      "load with no gradient is provably inert in the top BC.)", &
                      ierr, OCEAN_STATUS_ERR_SETUP)
            return
         end if
      end if

      ! P5.2: assemble the top-of-column load.  Rebuilt from scratch (not
      ! `+=` onto whatever the earlier psurf seed left) so the result does
      ! not depend on which of the two seeds ran first, and over the WHOLE
      ! array including ghosts — `p_top` owes no halo exchange of its own
      ! precisely because both of its sources are already ghost-valid.
      if (.not. (size(ocean_state%multilayer%p_top, 1) == nx .and. &
                 size(ocean_state%multilayer%p_top, 2) == ny)) then
         call fail("configure_ocean_cavity: ms%p_top and metrics%p_ice_ref "// &
                   "have different shapes — the top-of-column load cannot be "// &
                   "assembled.", ierr, OCEAN_STATUS_ERR_SETUP)
         return
      end if
      do j = 1, ny
         do i = 1, nx
            ocean_state%multilayer%p_top(i, j) = ocean_state%metrics%p_ice_ref(i, j)
         end do
      end do
      if (allocated(ocean_state%surface_flux%p_surf)) then
         do j = 1, ny
            do i = 1, nx
               ocean_state%multilayer%p_top(i, j) = &
                  ocean_state%multilayer%p_top(i, j) + &
                  ocean_state%surface_flux%p_surf(i, j)
            end do
         end do
      end if

      ! ---- The vertical coordinate's rigid-top seam (P6.2) ----
      ! `vcoord%z_top(i,j)` is the geopotential depth of the top of the
      ! WATER column — the ice base under a shelf, `z = 0` elsewhere.
      ! The draft is static, so this is a configure-time copy and the
      ! vcoord slot stays self-contained at run time: the ALE remap
      ! driver never sees `metrics` and the target builder's signature is
      ! unchanged.  Same pattern (and the same reason) as the sponge's
      ! `sp%z_top` in `configure_ocean_sponge`.
      !
      ! Consumed by the `VCOORD_Z_FIXED` branch of `compute_target_h`.
      ! Without a cavity this routine has already returned, so `z_top`
      ! keeps its init-time zero fill and every geometric family
      ! reproduces its pre-cavity arithmetic bit-for-bit.  BEFORE
      ! `ocean_state_enter_data` maps it.
      if (allocated(ocean_state%vcoord%z_top)) then
         if (size(ocean_state%vcoord%z_top, 1) == nx .and. &
             size(ocean_state%vcoord%z_top, 2) == ny) then
            do j = 1, ny
               do i = 1, nx
                  ocean_state%vcoord%z_top(i, j) = ocean_state%metrics%z_draft(i, j)
               end do
            end do
         else
            call fail("configure_ocean_cavity: vcoord%z_top and metrics%z_draft "// &
                      "have different shapes — the vertical coordinate cannot see "// &
                      "the ice base.", ierr, OCEAN_STATUS_ERR_SETUP)
            return
         end if
      end if

      ! The counted-once invariant (I), in metres of reference depth, over
      ! the WET columns:
      !     rho*g*z_draft + (bt_H_ref - b)*rho*g == 0   <=>   bt_H_ref == b - z_draft.
      ! GROUNDED columns are excluded because they are LAND: every face
      ! metric on them is zero, so they carry no barotropic momentum
      ! equation and there is no load on them to count once or twice.
      ! Their datum is deliberately 0 (`cavity_datum_impl`).
      ! Asserted on the common positive factor divided out — scale-free,
      ! and it does not fabricate a product the code never forms.  The
      ! bound is a pure round-off allowance on the ONE subtraction
      ! `b - z_draft`: both operands are O(max depth), so the result
      ! carries at most a few ulp of it.  (It is NOT asserted as bit-zero:
      ! the latch and this check evaluate the same difference in two
      ! places, and an FMA-contracting build is free to round them
      ! differently.)
      if (cfg%ocean%bt%n_inner >= 1) then
         datum_tol = 8.0_wp*epsilon(1.0_wp)* &
                     max(maxval(abs(ocean_state%barotropic%b)), 1.0_wp)
         resid = cavity_datum_residual(ocean_state%dyn%bt_work%bt_H_ref, &
                                       ocean_state%barotropic%b, &
                                       ocean_state%metrics%z_draft, &
                                       cfg%ocean%cavity_dyn%h_min_cavity, nx, ny)
         if (.not. (resid <= datum_tol)) then
            call fail("&ocean_cavity_dyn_nml: the barotropic datum and the ice "// &
                      "draft disagree (max |bt_H_ref - (b - z_draft)| = "// &
                      to_string(resid)//" m > "//to_string(datum_tol)//" m).  The "// &
                      "ice load would then be counted twice, or not at all — check "// &
                      "that z_draft went through the same periodic/fold re-wrap and "// &
                      "halo exchange as b.", ierr, OCEAN_STATUS_ERR_SETUP)
            return
         end if
      end if

      ! ---- Cover mask on the atmospheric forcing (P2c) ----
      ! There is no atmosphere under an ice shelf.  Two static forcing
      ! fields are masked ONCE, here, because this is the first point at
      ! which `cover_frac` exists AND the forcing has been seeded
      ! (`configure_ocean_forcing` runs much earlier, before the cavity
      ! geometry):
      !
      !   * the wind-stress PAIR, masked on every face touching a covered
      !     cell and followed by the `stress_mag` refresh in the same
      !     call — see `ocean_surface_stress_apply_cover`.  Masking the
      !     source rather than the derived views is what also silences
      !     the implicit vdiff stress fold and the MLE front sampler,
      !     which read `ss%tau_x` raw;
      !   * the SCALAR `&ocean_thermo_nml q_heat` / `q_salt` fill, but
      !     ONLY when the component set is off.  With components on those
      !     two are assembler outputs and the mask belongs in
      !     `ocean_surface_flux_assemble` (a second writer here would
      !     break the fill contract); the call below no-ops itself in
      !     that case.
      !
      ! Both are host-side and both run BEFORE `ocean_state_enter_data`,
      ! so the masked values are what the device map captures.  Both are
      ! idempotent.  The time-varying twin of the wind mask lives in
      ! `ocean_seam_refresh_surface_stress` (per data-forcing bracket).
      call ocean_surface_stress_apply_cover(ocean_state%surface_stress, &
                                            ocean_state%metrics%cover_frac)
      call ocean_surface_flux_apply_cover_const(ocean_state%surface_flux, &
                                                ocean_state%metrics%cover_frac)

      if (compute_rank == 0) then
         call logger%info("Ice-shelf cavity: isostatic load p_ice_ref = "// &
                          "rho_ref*g*z_draft, max = "// &
                          to_string(maxval(ocean_state%metrics%p_ice_ref))// &
                          " Pa (rho_ref = "// &
                          to_string(ocean_state%pressure_force%rho_ref)//" kg/m^3)")
         call logger%info("                  assembled into ms%p_top (max = "// &
                          to_string(maxval(ocean_state%multilayer%p_top))// &
                          " Pa = p_ice_ref + sf%p_surf); NOT into sf%p_surf, "// &
                          "which the eta_forcing seam is built from")
         if (ocean_state%pressure_force%p_top_in_bc) then
            call logger%info("                  consumed by the FV_MOM6 pa(nz+1) "// &
                             "surface BC (&ocean_pgf_nml p_top_in_bc)")
         else
            call logger%info("                  the draft is UNIFORM here, so the "// &
                             "load has no gradient and the PGF top BC is provably "// &
                             "inert; &ocean_pgf_nml p_top_in_bc is not required")
         end if
      end if
      if (present(ierr)) ierr = OCEAN_STATUS_OK
   end subroutine configure_ocean_cavity