ocean_cavity_mass_step Subroutine

   subroutine ocean_cavity_mass_step(grid, metrics, cav, ms, dt, weight, active)
      !! The real-freshwater MASS update — one call, at the THERMO
      !! cadence, from inside the RK2 stage immediately after
      !! `ocean_surface_flux_apply_tracers`.
      !!
      !! WHY THERE, and not in `engine_step_finalize` next to the melt
      !! solve.  Three things have to line up:
      !!
      !!   1. **The same melt value.**  The salt and heat halves of the
      !!      melt ride `Q_salt`/`Q_heat`, which the stage's surface-flux
      !!      apply spends.  Putting the volume anywhere else would spend
      !!      `melt(n)` for mass and `melt(n-1)` for salt.
      !!   2. **The same stage weight, by construction.**  Under
      !!      `ssp_rk2` both stages apply and `rk2_average` halves the
      !!      pair; under `pred_corr` `therm_active` is false in the
      !!      predictor and the corrector's single application IS the
      !!      step.  The caller passes the matching `weight` for the
      !!      budget accumulator, the same one `ocean_accumulate_mass_out`
      !!      takes.
      !!   3. **The barotropic mode sees it.**  `derive_bt_from_layers`
      !!      rebuilds `bt_eta = sum_k h - bt_H_ref` at the TOP of every
      !!      stage, so a thickness source applied here is in `bt_eta`
      !!      one stage later with no separate barotropic forcing term —
      !!      the free surface under the cavity datum simply rises.  The
      !!      BT substep's transport renormalisation (`bt_uhbt`) is a
      !!      constraint on the layer TRANSPORTS within a stage and never
      !!      reads a thickness source, so it is undisturbed.  This is
      !!      the documented `F_slow`-style operator split, and it is what
      !!      MOM6 does with its own surface mass fluxes.
      !!
      !! It runs BEFORE the vertical-mixing block in the same stage, so
      !! the implicit vdiff/drag solves see the thickened top layer —
      !! which is the right order: the added volume is part of the column
      !! before the column is mixed.  The ALE remap then redistributes it
      !! to the coordinate's target after `rk2_average`.
      !!
      !! NON-`pure`: it logs, it fails loud, it mutates `ms%mass_src`,
      !! and (with compensation on) it issues one collective.  Its six
      !! kernels are all `pure`.
      !!
      !! `mem:separate`: every array is mapped by its owning slot's
      !! `enter_data`; the tracer-registry derefs happen on the HOST
      !! before each kernel (the outer-shim rule).
      use pic_logger, only: global_logger
      use pic_strings, only: to_string
      use rdb_error_ring, only: fail
      use rdb_ocean_status, only: OCEAN_STATUS_ERR_SETUP
      use rdb_halo, only: halo_allreduce_sum
      type(hgrid_t), intent(in) :: grid
         !! Horizontal grid (`nx_total`/`ny_total`/`nghost`).
      type(ocean_metrics_t), intent(in) :: metrics
         !! Reads `cover_frac` and `areaT`.
      type(ocean_cavity_flux_t), intent(inout) :: cav
         !! The cavity-melt slot; reads `melt`/`active`/`t_b`, writes
         !! `comp_scale` and the thin counters.
      type(multilayer_state_t), intent(inout) :: ms
         !! Writes `h_layer(:,:,nz)`, the S/T/passive tracer loads, the
         !! two surface budget contributors and `mass_src`.
      real(wp), intent(in) :: dt
         !! Thermo timestep (s) — the same `therm_dt` the surface-flux
         !! apply was given.
      real(wp), intent(in) :: weight
         !! Per-stage weight for the `mass_src` accumulator (0.5 per
         !! SSP-RK2 stage; 0 / 1 for the pred_corr predictor /
         !! corrector), matching `ocean_accumulate_mass_out`.
      logical, intent(in), optional :: active
         !! Thermo-cadence gate.  Present-and-false ⇒ early return;
         !! absent ⇒ run.
      integer :: nx, ny, nz, idx_t, idx_s, idx_ps, it, n_thin_src, n_thin_comp
      real(wp) :: dt_over_rho0, vol_melt, area_open, dw, tmp

      if (.not. cav%is_init) return
      if (.not. cav%enable) return
      if (cav%freshwater /= CAVITY_FW_MASS) return
      if (present(active)) then
         if (.not. active) return
      end if
      if (.not. allocated(ms%tracers)) return
      idx_t = ms%idx_temperature
      idx_s = ms%idx_salinity
      if (idx_t <= 0 .or. idx_s <= 0) return

      nx = grid%nx_total
      ny = grid%ny_total
      nz = ms%nz_ml
      idx_ps = ms%idx_pseudo_salt
      dt_over_rho0 = dt/cav%rho0

      ! (1) Interior integrals FIRST: the melt volume the budget tracks
      ! and the open-ocean area the sink would spend it over.  Combined
      ! across ranks so a decomposed domain removes the same total the
      ! whole domain gained.
      call cavity_mass_totals_impl(nx, ny, grid%nghost, dt_over_rho0, cav%active, &
                                   ms%wet_mask, metrics%cover_frac, cav%melt, &
                                   metrics%areaT, vol_melt, area_open)
      tmp = vol_melt
      call halo_allreduce_sum(tmp, vol_melt)
      tmp = area_open
      call halo_allreduce_sum(tmp, area_open)
      cav%melt_volume_step = vol_melt
      cav%open_area = area_open

      ! (2) The source.
      call cavity_mass_apply_impl(nx, ny, nz, dt_over_rho0, dt_over_rho0, cav%s_ice, &
                                  cav%active, cav%melt, cav%s_far, cav%t_b, &
                                  ms%h_layer, ms%tracers(idx_s)%hTr, ms%tracers(idx_t)%hTr, &
                                  ms%salt_budget_surface, ms%heat_budget_surface, &
                                  ms%k_top, n_thin_src)
      if (idx_ps > 0) then
         call cavity_mass_salt_mirror_impl(nx, ny, nz, dt_over_rho0, dt_over_rho0, &
                                           cav%s_ice, cav%active, cav%melt, &
                                           cav%s_far, ms%tracers(idx_ps)%hTr, ms%k_top)
      end if
      ms%mass_src = ms%mass_src + weight*RHO_WATER*vol_melt

      ! (3) The sink, when the knob asks for it.
      n_thin_comp = 0
      cav%comp_withdrawal_step = 0.0_wp
      if (cav%volume_comp == CAVITY_VC_UNIFORM_OPEN) then
         dw = cavity_comp_withdrawal(vol_melt, area_open)
         cav%comp_withdrawal_step = dw
         call cavity_comp_apply_impl(nx, ny, nz, dw, ms%wet_mask, metrics%cover_frac, &
                                     ms%h_layer, ms%tracers(idx_s)%hTr, &
                                     ms%tracers(idx_t)%hTr, ms%salt_budget_surface, &
                                     ms%heat_budget_surface, cav%comp_scale, n_thin_comp)
         ! NOTE the sink stays on `k = nz` and is NOT routed through
         ! `k_top`, deliberately: its own gate is `cover_frac < 0.5`, so
         ! it only ever acts on OPEN-OCEAN columns, and an open-ocean
         ! column has `z_top = 0` ⇒ no top-side filler ⇒ `k_top ≡ nz`.
         ! Routing it would be a provable no-op; leaving it spells out
         ! that "the top layer" and "the first live layer" are the same
         ! row wherever this kernel runs.  `cavity_comp_scale_tracer_impl`
         ! below rides the same argument (it is unconditional, but
         ! `comp_scale` is exactly 1 off the sink).
         do it = 1, size(ms%tracers)
            if (it == idx_s .or. it == idx_t) cycle
            call cavity_comp_scale_tracer_impl(nx, ny, nz, cav%comp_scale, &
                                               ms%tracers(it)%hTr)
         end do
         ! The removed interior volume is exactly `dw*area_open` — the
         ! same product the withdrawal was derived from — so with no
         ! clamping this cancels the source term to the last bit and the
         ! domain mass is constant.
         ms%mass_src = ms%mass_src - weight*RHO_WATER*dw*area_open
      end if

      ! (4) Accounting, then fail loud.
      cav%n_thin_step = n_thin_src + n_thin_comp
      cav%n_thin_total = cav%n_thin_total + int(cav%n_thin_step, int64)
      if (cavity_mass_thin_is_fatal(cav%n_thin_step)) then
         call global_logger%error("cavity real freshwater: "// &
                                  to_string(n_thin_src)//" melt column(s) and "// &
                                  to_string(n_thin_comp)//" compensation column(s) "// &
                                  "could not give up the requested thickness")
         call fail("&ocean_cavity_melt_nml freshwater='mass': "// &
                   to_string(cav%n_thin_step)//" column(s) would have been driven "// &
                   "below H_CAVITY_FLOOR (2*H_VANISHED) by the top-layer "// &
                   "withdrawal (a freezing "// &
                   "column thinner than |m|*dt/rho_0, or a compensation sink "// &
                   "deeper than the open-ocean top layer).  The withdrawal was "// &
                   "clamped so the state stays finite, but a clamped withdrawal "// &
                   "no longer matches the tracked mass source, so the console's "// &
                   "Mass Error would stop being a leak measurement while still "// &
                   "being printed as one.  Use a thicker top layer (a coordinate "// &
                   "with a larger surface target), a shorter dt, or "// &
                   "volume_compensation='none'.", code=OCEAN_STATUS_ERR_SETUP)
      end if
   end subroutine ocean_cavity_mass_step