continuity_gm_apply Subroutine

public subroutine continuity_gm_apply(grid, metrics, this, ms, gm, dt, budget_w, tracer_mode, bc, set_flux_h)

Gent-McWilliams thickness diffusion as its OWN sequential operator: move h_layer AND every tracer by the bolus transport gm%uhD/ gm%vhD, which gm_compute_transports has JUST filled from this same, untouched h_layer with this same dt.

MOM6 parity: thickness_diffuse runs after step_MOM_dyn_split_RK2 and updates h in place, h -= dt·IareaT·(div uhD), while adding uhD·dt to uhtr so the SAME tracer advection that carries the resolved transport carries the bolus one. Here:

  • every-step tracers (TR_MODE_ADVECT): the bolus flux goes through the resolved path’s own PPM tracer kernels (tracer_advect_{zonal,meridional}_one_impl), interleaved with the two thickness applies exactly as in continuity_tracer_step_split (x: advect at h, apply; seam refresh; y: advect at h*, apply) — so a uniform tracer stays uniform (CWC) and the content is conserved to round-off;
  • windowed tracers (TR_MODE_ACCUMULATE): uhD·dt/vhD·dt join the window accumulator uhtr/vhtr (MOM6’s uhtr += uhD·dt) and the concentration hold is re-weighted onto the new h, so the drain spends the bolus transport with the resolved one.

Positivity: uhD is capped per face by A·(h − H_VANISHED)/(4·dt) of the DONOR at the h passed in, so the four faces of a cell remove at most h − H_VANISHED over dt — no layer is taken below min(h, H_VANISHED), whatever the dynamics left. (Folded into the resolved sweeps, as until 2026-10, the cap bounded the stage-entry h and the resolved outflow came on top — see rdb_ocean_gm.)

Edges: the bolus transport is ZEROED on every physical edge face that is not periodic or a tripolar fold — walls, sponges and every open-boundary type — as MOM6 masks the GM slopes / KhTh by OBCmaskCu/Cv: no GM flux leaves the domain, so no budget term and no ghost fill is needed. An MPI seam (has_* false) is interior and keeps its transport.

I1′: fillers keep their donor’s concentration through this operator — the PPM kernel reads a filler’s hTr/h, which IS c_live under I1′, and a filler cannot DONATE (its availability is 0); what it receives arrives at its live neighbour’s concentration. The pool (multilayer_state_t%enforce_vanished_content) at the tail of the outer step restores I1′ exactly, as after the resolved continuity.

budget_w multiplies the heat/salt horizontal-advection budget increments: this operator runs AFTER the RK2 stage average, so it records 1/ocean_budget_stage_weight (2 under ssp_rk2, 1 under pred_corr) for the console’s per-step weight to recover it 1:1.

set_flux_h (optional, default false): leave the bolus thickness divergence in ms%flux_h_layer for the eulerian_z vertical advection (compute_w_from_continuity), which then cancels it per layer exactly as it cancels the resolved divergence.

The caller refreshes the h / tracer ghosts afterwards (ocean_halo_exchange_ml_state + periodic wrap + fold). No-op when GM is uninitialised or disabled.

Arguments

Type IntentOptional Attributes Name
type(hgrid_t), intent(in) :: grid
type(ocean_metrics_t), intent(in) :: metrics
type(continuity_t), intent(inout) :: this
type(multilayer_state_t), intent(inout) :: ms
type(ocean_gm_t), intent(inout) :: gm

uhD/vhD (inout: the edge closure and the fold-line projection are applied to them in place).
real(kind=wp), intent(in) :: dt

The step the transports were capped for (s).
real(kind=wp), intent(in) :: budget_w

Budget bookkeeping weight (see above).
integer, intent(in) :: tracer_mode

TR_MODE_ADVECT or TR_MODE_ACCUMULATE.
type(ocean_bc_state_t), intent(in), optional :: bc

Edge tags (absent ⇒ every edge a wall).
logical, intent(in), optional :: set_flux_h

Leave the bolus divergence in ms%flux_h_layer (eulerian_z).

Calls

proc~~continuity_gm_apply~~CallsGraph proc~continuity_gm_apply continuity_gm_apply interface~ocean_fold_north_v_face ocean_fold_north_v_face proc~continuity_gm_apply->interface~ocean_fold_north_v_face interface~ocean_halo_centre ocean_halo_centre proc~continuity_gm_apply->interface~ocean_halo_centre local local proc~continuity_gm_apply->local proc~drain_copy_3d drain_copy_3d proc~continuity_gm_apply->proc~drain_copy_3d proc~drain_rescale_htr drain_rescale_hTr proc~continuity_gm_apply->proc~drain_rescale_htr proc~drain_rescale_htr_budget drain_rescale_hTr_budget proc~continuity_gm_apply->proc~drain_rescale_htr_budget proc~gm_tracer_advect_x gm_tracer_advect_x proc~continuity_gm_apply->proc~gm_tracer_advect_x proc~gm_tracer_advect_y gm_tracer_advect_y proc~continuity_gm_apply->proc~gm_tracer_advect_y proc~ocean_bc_outer_face_tag ocean_bc_outer_face_tag proc~continuity_gm_apply->proc~ocean_bc_outer_face_tag proc~ocean_fold_wrap_centre_3d_state ocean_fold_wrap_centre_3d_state proc~continuity_gm_apply->proc~ocean_fold_wrap_centre_3d_state proc~ocean_halo_is_decomposed_x ocean_halo_is_decomposed_x proc~continuity_gm_apply->proc~ocean_halo_is_decomposed_x proc~ocean_halo_is_decomposed_y ocean_halo_is_decomposed_y proc~continuity_gm_apply->proc~ocean_halo_is_decomposed_y proc~ocean_periodic_wrap_centre_3d ocean_periodic_wrap_centre_3d proc~continuity_gm_apply->proc~ocean_periodic_wrap_centre_3d proc~profiler_start profiler_start proc~continuity_gm_apply->proc~profiler_start proc~profiler_stop profiler_stop proc~continuity_gm_apply->proc~profiler_stop proc~fold_v_2d fold_v_2d interface~ocean_fold_north_v_face->proc~fold_v_2d proc~fold_v_3d fold_v_3d interface~ocean_fold_north_v_face->proc~fold_v_3d proc~ocean_halo_centre_2d ocean_halo_centre_2d interface~ocean_halo_centre->proc~ocean_halo_centre_2d proc~ocean_halo_centre_3d ocean_halo_centre_3d interface~ocean_halo_centre->proc~ocean_halo_centre_3d proc~drain_rescale_htr_budget->local proc~tracer_advect_zonal_one_impl tracer_advect_zonal_one_impl proc~gm_tracer_advect_x->proc~tracer_advect_zonal_one_impl proc~tracer_advect_meridional_one_impl tracer_advect_meridional_one_impl proc~gm_tracer_advect_y->proc~tracer_advect_meridional_one_impl interface~fold_north_centre fold_north_centre proc~ocean_fold_wrap_centre_3d_state->interface~fold_north_centre interface~ocean_fold_pack ocean_fold_pack proc~ocean_fold_wrap_centre_3d_state->interface~ocean_fold_pack interface~ocean_fold_unpack ocean_fold_unpack proc~ocean_fold_wrap_centre_3d_state->interface~ocean_fold_unpack proc~n_tracers n_tracers proc~ocean_fold_wrap_centre_3d_state->proc~n_tracers proc~ocean_fold_begin ocean_fold_begin proc~ocean_fold_wrap_centre_3d_state->proc~ocean_fold_begin proc~ocean_fold_end ocean_fold_end proc~ocean_fold_wrap_centre_3d_state->proc~ocean_fold_end proc~ocean_fold_exchange ocean_fold_exchange proc~ocean_fold_wrap_centre_3d_state->proc~ocean_fold_exchange proc~ocean_fold_is_distributed ocean_fold_is_distributed proc~ocean_fold_wrap_centre_3d_state->proc~ocean_fold_is_distributed proc~find_or_create_region find_or_create_region proc~profiler_start->proc~find_or_create_region proc~get_wall_time get_wall_time proc~profiler_start->proc~get_wall_time proc~nvtx_range_push nvtx_range_push proc~profiler_start->proc~nvtx_range_push proc~profiler_stop->proc~get_wall_time proc~nvtx_range_pop nvtx_range_pop proc~profiler_stop->proc~nvtx_range_pop proc~fold_north_centre_2d fold_north_centre_2d interface~fold_north_centre->proc~fold_north_centre_2d proc~fold_north_centre_3d fold_north_centre_3d interface~fold_north_centre->proc~fold_north_centre_3d proc~ocean_fold_pack_2d ocean_fold_pack_2d interface~ocean_fold_pack->proc~ocean_fold_pack_2d proc~ocean_fold_pack_3d ocean_fold_pack_3d interface~ocean_fold_pack->proc~ocean_fold_pack_3d proc~ocean_fold_unpack_2d ocean_fold_unpack_2d interface~ocean_fold_unpack->proc~ocean_fold_unpack_2d proc~ocean_fold_unpack_3d ocean_fold_unpack_3d interface~ocean_fold_unpack->proc~ocean_fold_unpack_3d proc~fold_v_2d->interface~ocean_fold_pack proc~fold_v_2d->interface~ocean_fold_unpack proc~fold_v_2d->proc~ocean_fold_begin proc~fold_v_2d->proc~ocean_fold_end proc~fold_v_2d->proc~ocean_fold_exchange proc~fold_v_2d->proc~ocean_fold_is_distributed interface~fold_north_v_face fold_north_v_face proc~fold_v_2d->interface~fold_north_v_face proc~fold_stagger_nrows fold_stagger_nrows proc~fold_v_2d->proc~fold_stagger_nrows proc~fold_v_3d->interface~ocean_fold_pack proc~fold_v_3d->interface~ocean_fold_unpack proc~fold_v_3d->proc~ocean_fold_begin proc~fold_v_3d->proc~ocean_fold_end proc~fold_v_3d->proc~ocean_fold_exchange proc~fold_v_3d->proc~ocean_fold_is_distributed proc~fold_v_3d->interface~fold_north_v_face proc~fold_v_3d->proc~fold_stagger_nrows proc~grow_buffers grow_buffers proc~ocean_fold_begin->proc~grow_buffers to_string to_string proc~ocean_fold_begin->to_string warning warning proc~ocean_fold_begin->warning comm_irecv_real_sp_array_n comm_irecv_real_sp_array_n proc~ocean_fold_exchange->comm_irecv_real_sp_array_n comm_isend_real_sp_array_n comm_isend_real_sp_array_n proc~ocean_fold_exchange->comm_isend_real_sp_array_n proc~comm_env_compute_comm comm_env_compute_comm proc~ocean_fold_exchange->proc~comm_env_compute_comm waitall waitall proc~ocean_fold_exchange->waitall proc~ocean_halo_centre_2d_impl ocean_halo_centre_2d_impl proc~ocean_halo_centre_2d->proc~ocean_halo_centre_2d_impl proc~oh_count_centre_2d oh_count_centre_2d proc~ocean_halo_centre_2d->proc~oh_count_centre_2d proc~ocean_halo_centre_3d->proc~ocean_periodic_wrap_centre_3d proc~ocean_halo_centre_3d->comm_irecv_real_sp_array_n proc~ocean_halo_centre_3d->comm_isend_real_sp_array_n proc~ocean_halo_centre_3d->proc~comm_env_compute_comm proc~ew_rank_east ew_rank_east proc~ocean_halo_centre_3d->proc~ew_rank_east proc~ew_rank_west ew_rank_west proc~ocean_halo_centre_3d->proc~ew_rank_west proc~needs_flags needs_flags proc~ocean_halo_centre_3d->proc~needs_flags proc~ns_rank_north ns_rank_north proc~ocean_halo_centre_3d->proc~ns_rank_north proc~ns_rank_south ns_rank_south proc~ocean_halo_centre_3d->proc~ns_rank_south proc~ocean_halo_buffers_ensure_nz ocean_halo_buffers_ensure_nz proc~ocean_halo_centre_3d->proc~ocean_halo_buffers_ensure_nz proc~oh_count_centre_3d oh_count_centre_3d proc~ocean_halo_centre_3d->proc~oh_count_centre_3d proc~oh_count_msgs oh_count_msgs proc~ocean_halo_centre_3d->proc~oh_count_msgs proc~ocean_halo_centre_3d->waitall proc~tracer_advect_meridional_one_impl->local proc~ppm_cell_limiter ppm_cell_limiter proc~tracer_advect_meridional_one_impl->proc~ppm_cell_limiter proc~ppm_limited_slope ppm_limited_slope proc~tracer_advect_meridional_one_impl->proc~ppm_limited_slope proc~ppm_mirror_h ppm_mirror_h proc~tracer_advect_meridional_one_impl->proc~ppm_mirror_h proc~tracer_advect_zonal_one_impl->local proc~tracer_advect_zonal_one_impl->proc~ppm_cell_limiter proc~tracer_advect_zonal_one_impl->proc~ppm_limited_slope proc~tracer_advect_zonal_one_impl->proc~ppm_mirror_h proc~fold_north_v_face_2d fold_north_v_face_2d interface~fold_north_v_face->proc~fold_north_v_face_2d proc~fold_north_v_face_3d fold_north_v_face_3d interface~fold_north_v_face->proc~fold_north_v_face_3d comm_world comm_world proc~comm_env_compute_comm->comm_world proc~decomp_rank_from_coords decomp_rank_from_coords proc~ew_rank_east->proc~decomp_rank_from_coords proc~ew_rank_west->proc~decomp_rank_from_coords proc~ns_rank_north->proc~decomp_rank_from_coords proc~ns_rank_south->proc~decomp_rank_from_coords proc~ocean_fold_pack_2d->proc~ocean_fold_pack_3d proc~ocean_fold_pack_3d->proc~fold_stagger_nrows proc~fold_stagger_family fold_stagger_family proc~ocean_fold_pack_3d->proc~fold_stagger_family proc~ocean_fold_unpack_2d->proc~ocean_fold_unpack_3d proc~ocean_fold_unpack_3d->proc~fold_stagger_nrows proc~ocean_fold_unpack_3d->proc~fold_stagger_family proc~ocean_halo_buffers_ensure_nz->to_string proc~ocean_halo_buffers_ensure_nz->warning proc~ocean_halo_centre_2d_impl->comm_irecv_real_sp_array_n proc~ocean_halo_centre_2d_impl->comm_isend_real_sp_array_n proc~ocean_halo_centre_2d_impl->proc~comm_env_compute_comm proc~ocean_halo_centre_2d_impl->proc~ew_rank_east proc~ocean_halo_centre_2d_impl->proc~ew_rank_west proc~ocean_halo_centre_2d_impl->proc~needs_flags proc~ocean_halo_centre_2d_impl->proc~ns_rank_north proc~ocean_halo_centre_2d_impl->proc~ns_rank_south proc~ocean_halo_centre_2d_impl->proc~oh_count_msgs proc~ocean_halo_centre_2d_impl->waitall proc~ocean_periodic_wrap_centre_2d ocean_periodic_wrap_centre_2d proc~ocean_halo_centre_2d_impl->proc~ocean_periodic_wrap_centre_2d proc~fold_north_v_face_2d->local proc~fold_north_v_face_3d->local

Called by

proc~~continuity_gm_apply~~CalledByGraph proc~continuity_gm_apply continuity_gm_apply proc~run_gm_step run_gm_step proc~run_gm_step->proc~continuity_gm_apply proc~ocean_dyn_step_split ocean_dyn_step_split proc~ocean_dyn_step_split->proc~run_gm_step proc~engine_step engine_step proc~engine_step->proc~ocean_dyn_step_split proc~driver_run_ocean driver_run_ocean proc~driver_run_ocean->proc~engine_step proc~rdb_ocean_step rdb_ocean_step proc~rdb_ocean_step->proc~engine_step proc~driver_run driver_run proc~driver_run->proc~driver_run_ocean

Variables

Type Visibility Attributes Name Initial
real(kind=wp), private :: div
integer, private :: i
integer, private :: it
integer, private :: j
integer, private :: k
logical, private :: keep_e
logical, private :: keep_n
logical, private :: keep_s
logical, private :: keep_w
integer, private :: nghost
integer, private :: nx
integer, private :: nx_phys
integer, private :: ny
integer, private :: ny_phys
integer, private :: nz
logical, private :: per_x
logical, private :: per_y
integer, private :: tag_e
integer, private :: tag_n
integer, private :: tag_s
integer, private :: tag_w
logical, private :: want_flux_h

Source Code

   subroutine continuity_gm_apply(grid, metrics, this, ms, gm, dt, budget_w, tracer_mode, bc, &
                                  set_flux_h)
      !! Gent-McWilliams thickness diffusion as its OWN sequential operator:
      !! move `h_layer` AND every tracer by the bolus transport `gm%uhD`/
      !! `gm%vhD`, which `gm_compute_transports` has JUST filled from this
      !! same, untouched `h_layer` with this same `dt`.
      !!
      !! MOM6 parity: `thickness_diffuse` runs after `step_MOM_dyn_split_RK2`
      !! and updates `h` in place, `h -= dt·IareaT·(div uhD)`, while
      !! adding `uhD·dt` to `uhtr` so the SAME tracer advection that carries
      !! the resolved transport carries the bolus one.  Here:
      !!
      !!   * every-step tracers (`TR_MODE_ADVECT`): the bolus flux goes
      !!     through the resolved path's own PPM tracer kernels
      !!     (`tracer_advect_{zonal,meridional}_one_impl`), interleaved with
      !!     the two thickness applies exactly as in
      !!     `continuity_tracer_step_split` (x: advect at h, apply; seam
      !!     refresh; y: advect at h*, apply) — so a uniform tracer stays
      !!     uniform (CWC) and the content is conserved to round-off;
      !!   * windowed tracers (`TR_MODE_ACCUMULATE`): `uhD·dt`/`vhD·dt` join
      !!     the window accumulator `uhtr`/`vhtr` (MOM6's `uhtr += uhD·dt`)
      !!     and the concentration hold is re-weighted onto the new `h`, so
      !!     the drain spends the bolus transport with the resolved one.
      !!
      !! Positivity: `uhD` is capped per face by `A·(h − H_VANISHED)/(4·dt)` of
      !! the DONOR at the `h` passed in, so the four faces of a cell remove at
      !! most `h − H_VANISHED` over `dt` — no layer is taken below
      !! `min(h, H_VANISHED)`, whatever the dynamics left.  (Folded into the
      !! resolved sweeps, as until 2026-10, the cap bounded the stage-entry
      !! `h` and the resolved outflow came on top — see `rdb_ocean_gm`.)
      !!
      !! Edges: the bolus transport is ZEROED on every physical edge face
      !! that is not periodic or a tripolar fold — walls, sponges and every
      !! open-boundary type — as MOM6 masks the GM slopes / KhTh by
      !! `OBCmaskCu/Cv`: no GM flux leaves the domain, so no budget term
      !! and no ghost fill is needed.  An MPI seam (`has_*` false) is
      !! interior and keeps its transport.
      !!
      !! I1′: fillers keep their donor's concentration through this operator
      !! — the PPM kernel reads a filler's `hTr/h`, which IS `c_live` under
      !! I1′, and a filler cannot DONATE (its availability is 0); what it
      !! receives arrives at its live neighbour's concentration.  The pool
      !! (`multilayer_state_t%enforce_vanished_content`) at the tail of the
      !! outer step restores I1′ exactly, as after the resolved continuity.
      !!
      !! `budget_w` multiplies the heat/salt horizontal-advection budget
      !! increments: this operator runs AFTER the RK2 stage average, so it
      !! records `1/ocean_budget_stage_weight` (2 under ssp_rk2, 1 under
      !! pred_corr) for the console's per-step weight to recover it 1:1.
      !!
      !! `set_flux_h` (optional, default false): leave the bolus thickness
      !! divergence in `ms%flux_h_layer` for the `eulerian_z` vertical
      !! advection (`compute_w_from_continuity`), which then cancels it per
      !! layer exactly as it cancels the resolved divergence.
      !!
      !! The caller refreshes the h / tracer ghosts afterwards
      !! (`ocean_halo_exchange_ml_state` + periodic wrap + fold).  No-op when
      !! GM is uninitialised or disabled.
      type(hgrid_t), intent(in) :: grid
      type(ocean_metrics_t), intent(in) :: metrics
      type(continuity_t), intent(inout) :: this
      type(multilayer_state_t), intent(inout) :: ms
      type(ocean_gm_t), intent(inout) :: gm
         !! `uhD`/`vhD` (inout: the edge closure and the fold-line
         !! projection are applied to them in place).
      real(wp), intent(in) :: dt
         !! The step the transports were capped for (s).
      real(wp), intent(in) :: budget_w
         !! Budget bookkeeping weight (see above).
      integer, intent(in) :: tracer_mode
         !! `TR_MODE_ADVECT` or `TR_MODE_ACCUMULATE`.
      type(ocean_bc_state_t), intent(in), optional :: bc
         !! Edge tags (absent ⇒ every edge a wall).
      logical, intent(in), optional :: set_flux_h
         !! Leave the bolus divergence in `ms%flux_h_layer` (eulerian_z).

      integer :: nx, ny, nz, nx_phys, ny_phys, nghost, it
      integer :: i, j, k
      integer :: tag_w, tag_e, tag_s, tag_n
      logical :: keep_w, keep_e, keep_s, keep_n, per_x, per_y, want_flux_h
      real(wp) :: div

      if (.not. gm%is_init) return
      if (.not. gm%enable) return
      if (.not. allocated(gm%uhD)) return

      nx = grid%nx_total
      ny = grid%ny_total
      nz = ms%nz_ml
      nx_phys = grid%nx_phys
      ny_phys = grid%ny_phys
      nghost = grid%nghost
      want_flux_h = .false.
      if (present(set_flux_h)) want_flux_h = set_flux_h

      ! ---- Edge closure (MOM6 OBCmaskCu/Cv): keep only periodic / fold
      ! edges and MPI seams.
      tag_w = OBC_WALL
      tag_e = OBC_WALL
      tag_s = OBC_WALL
      tag_n = OBC_WALL
      per_x = .false.
      per_y = .false.
      if (present(bc)) then
         tag_w = ocean_bc_outer_face_tag(bc%west%bc_type)
         tag_e = ocean_bc_outer_face_tag(bc%east%bc_type)
         tag_s = ocean_bc_outer_face_tag(bc%south%bc_type)
         tag_n = ocean_bc_outer_face_tag(bc%north%bc_type)
         if (.not. bc%has_west) tag_w = OBC_PERIODIC
         if (.not. bc%has_east) tag_e = OBC_PERIODIC
         if (.not. bc%has_south) tag_s = OBC_PERIODIC
         if (.not. bc%has_north) tag_n = OBC_PERIODIC
         per_x = bc%periodic_x .and. .not. ocean_halo_is_decomposed_x()
         per_y = bc%periodic_y .and. .not. ocean_halo_is_decomposed_y()
      end if
      keep_w = tag_w == OBC_PERIODIC .or. tag_w == OBC_TRIPOLAR_FOLD
      keep_e = tag_e == OBC_PERIODIC .or. tag_e == OBC_TRIPOLAR_FOLD
      keep_s = tag_s == OBC_PERIODIC .or. tag_s == OBC_TRIPOLAR_FOLD
      keep_n = tag_n == OBC_PERIODIC .or. tag_n == OBC_TRIPOLAR_FOLD
      do concurrent(k=1:nz, j=1:ny)
         if (.not. keep_w) gm%uhD(nghost + 1, j, k) = 0.0_wp
         if (.not. keep_e) gm%uhD(nghost + nx_phys + 1, j, k) = 0.0_wp
      end do
      do concurrent(k=1:nz, i=1:nx)
         if (.not. keep_s) gm%vhD(i, nghost + 1, k) = 0.0_wp
         if (.not. keep_n) gm%vhD(i, nghost + ny_phys + 1, k) = 0.0_wp
      end do
      ! Tripolar fold-line projection (see `continuity_tracer_step_split`):
      ! the duplicated fold-line face must carry ONE antisymmetric flux so
      ! what leaves (i,nj) north is exactly what enters (ni+1-i,nj).
      if (present(bc)) then
         if (bc%north_fold) call ocean_fold_north_v_face(gm%vhD, nx, ny + 1, nz, &
                                                         nx_phys, ny_phys, nghost)
      end if

      ! ---- x half: tracers at h^n, then h^n -> h*.
      if (tracer_mode == TR_MODE_ACCUMULATE) then
         call drain_copy_3d(nx, ny, nz, ms%h_layer, this%hprev_work)
         do concurrent(k=1:nz, j=1:ny, i=1:nx + 1)
            this%uhtr(i, j, k) = this%uhtr(i, j, k) + gm%uhD(i, j, k)*dt
         end do
      else
         call gm_tracer_advect_x(grid, metrics, this, ms, gm%uhD, dt, budget_w)
      end if
      do concurrent(k=1:nz, j=1:ny, i=1:nx) local(div)
         div = (gm%uhD(i + 1, j, k) - gm%uhD(i, j, k))*metrics%iareaT(i, j)
         ms%h_layer(i, j, k) = ms%h_layer(i, j, k) - dt*div
      end do
      if (want_flux_h) then
         do concurrent(k=1:nz, j=1:ny, i=1:nx)
            ms%flux_h_layer(i, j, k) = (gm%uhD(i + 1, j, k) - gm%uhD(i, j, k))*metrics%iareaT(i, j)
         end do
      end if

      ! ---- Mid-split seam refresh (the meridional PPM reads h / tracer
      ! ghosts the neighbour's x half just moved) — the same three steps as
      ! `continuity_tracer_step_split`.
      call profiler_start("ocean_comms_ml")
      call ocean_halo_centre(ms%h_layer, nz)
      if (allocated(ms%tracers)) then
         do it = 1, size(ms%tracers)
            if (.not. allocated(ms%tracers(it)%hTr)) cycle
            call ocean_halo_centre(ms%tracers(it)%hTr, nz)
         end do
      end if
      call profiler_stop("ocean_comms_ml")
      if (per_x .or. per_y) then
         call ocean_periodic_wrap_centre_3d(ms%h_layer, nx, ny, nz, &
                                            nx_phys, ny_phys, nghost, per_x, per_y, no_wait=.true.)
         if (allocated(ms%tracers)) then
            do it = 1, size(ms%tracers)
               if (.not. allocated(ms%tracers(it)%hTr)) cycle
               call ocean_periodic_wrap_centre_3d(ms%tracers(it)%hTr, nx, ny, nz, &
                                                  nx_phys, ny_phys, nghost, per_x, per_y, no_wait=.true.)
            end do
         end if
         !$acc wait(1)
      end if
      if (present(bc)) call ocean_fold_wrap_centre_3d_state(grid, bc, ms)

      ! ---- y half: tracers at h*, then h* -> h^{n+1}.
      if (tracer_mode == TR_MODE_ACCUMULATE) then
         do concurrent(k=1:nz, j=1:ny + 1, i=1:nx)
            this%vhtr(i, j, k) = this%vhtr(i, j, k) + gm%vhD(i, j, k)*dt
         end do
      else
         call gm_tracer_advect_y(grid, metrics, this, ms, gm%vhD, dt, budget_w)
      end if
      do concurrent(k=1:nz, j=1:ny, i=1:nx) local(div)
         div = (gm%vhD(i, j + 1, k) - gm%vhD(i, j, k))*metrics%iareaT(i, j)
         ms%h_layer(i, j, k) = ms%h_layer(i, j, k) - dt*div
      end do
      if (want_flux_h) then
         do concurrent(k=1:nz, j=1:ny, i=1:nx)
            ms%flux_h_layer(i, j, k) = ms%flux_h_layer(i, j, k) + &
                                       (gm%vhD(i, j + 1, k) - gm%vhD(i, j, k))*metrics%iareaT(i, j)
         end do
      end if

      ! ---- Windowed mode: re-weight the held concentration onto the new h
      ! (the same hold `continuity_tracer_step_split` keeps per stage), with
      ! the post-average budget weight.
      if (tracer_mode == TR_MODE_ACCUMULATE .and. allocated(ms%tracers)) then
         do it = 1, size(ms%tracers)
            if (.not. allocated(ms%tracers(it)%hTr)) cycle
            if (.not. ms%tracers(it)%do_horizontal_advection) cycle
            select case (ms%tracers(it)%budget_id)
            case (TRACER_BUDGET_HEAT)
               call drain_rescale_hTr_budget(nx, ny, nz, ms%h_layer, this%hprev_work, &
                                             budget_w, ms%tracers(it)%hTr, ms%heat_budget_horiz_adv)
            case (TRACER_BUDGET_SALT)
               call drain_rescale_hTr_budget(nx, ny, nz, ms%h_layer, this%hprev_work, &
                                             budget_w, ms%tracers(it)%hTr, ms%salt_budget_horiz_adv)
            case default
               call drain_rescale_hTr(nx, ny, nz, ms%h_layer, this%hprev_work, &
                                      ms%tracers(it)%hTr)
            end select
         end do
      end if
   end subroutine continuity_gm_apply