run_gm_step Subroutine

   subroutine run_gm_step(grid, metrics, dyn, ct, va, ms, dt, gm, slopes, varmix, &
                          wavespeed, vcoord, bc)
      !! Gent-McWilliams thickness diffusion as its OWN sequential operator,
      !! run once per outer step AFTER the dynamics (the stage loop and, under
      !! ssp_rk2, the stage average) — where MOM6 calls `thickness_diffuse`
      !! after `step_MOM_dyn_split_RK2`, every dynamics step:
      !!
      !!   1. `gm_compute_transports` fills `uhD`/`vhD` + `gm_src` from the
      !!      thickness the dynamics LEFT, with the stored slopes and the
      !!      VarMix/MEKE base KhTh refreshed at the top of the step;
      !!   2. `continuity_gm_apply` moves `h_layer` and every tracer by them
      !!      with the same `dt` — so the per-face availability cap
      !!      `A·(h − H_VANISHED)/(4·dt)` bounds what is actually there;
      !!   3. `eulerian_z` only: the bolus divergence is cancelled per layer by
      !!      the vertical advection, exactly as the resolved one is;
      !!   4. the h / tracer ghosts are refreshed (exchange, periodic wrap,
      !!      fold), as after the resolved continuity.
      !!
      !! Until 2026-10 the transports were computed at the top of the step
      !! from the stage-entry thickness and FOLDED into the resolved
      !! continuity sweeps; the cap then bounded the wrong `h` and a partial
      !! bed cell on an open z* step was driven negative (`rdb_ocean_gm`).
      !!
      !! Cadence: every outer step, with the outer `dt` (MOM6: every
      !! dynamics step).  The fold path applied GM on thermo steps only, so
      !! under `dt_therm_ratio > 1` it ran at 1/ratio of its strength.
      !! No-op when GM is absent / disabled, the slopes slot is absent, or
      !! thermodynamics is off.
      type(hgrid_t), intent(in) :: grid
      type(ocean_metrics_t), intent(in) :: metrics
      type(ocean_dyn_t), intent(inout) :: dyn
      type(continuity_t), intent(inout) :: ct
      type(ocean_vertical_advection_t), intent(inout) :: va
      type(multilayer_state_t), intent(inout) :: ms
      real(wp), intent(in) :: dt
      type(ocean_gm_t), intent(inout), optional :: gm
      type(ocean_slopes_t), intent(inout), optional :: slopes
      type(ocean_varmix_t), intent(inout), optional :: varmix
      type(ocean_wave_speed_t), intent(inout), optional :: wavespeed
      type(ocean_vcoord_t), intent(in), optional :: vcoord
      type(ocean_bc_state_t), intent(inout), optional :: bc

      integer :: tr_mode
      logical :: use_ext, is_eulerian
      real(wp) :: budget_w

      if (.not. present(gm)) return
      if (.not. present(slopes)) return
      if (.not. gm%enable) return
      if (.not. dyn%enable_thermodynamics) return

      call profiler_start("ocean_gm")
      use_ext = .false.
      if (present(varmix) .and. present(wavespeed)) use_ext = varmix%enable
      if (use_ext) then
         call gm_compute_transports(grid, metrics, gm, slopes, ms, dt, &
                                    khth_ext_u=varmix%khth_u, khth_ext_v=varmix%khth_v)
      else
         call gm_compute_transports(grid, metrics, gm, slopes, ms, dt)
      end if

      tr_mode = TR_MODE_ADVECT
      if (dyn%dt_tracer_advect_ratio > 1) tr_mode = TR_MODE_ACCUMULATE
      ! Post-average budget weight: the console multiplies the heat/salt
      ! accumulators by `ocean_budget_stage_weight` (0.5 ssp_rk2, 1
      ! pred_corr); this operator runs once, after the average.
      budget_w = 1.0_wp/ocean_budget_stage_weight(dyn%split_scheme == SPLIT_SCHEME_PRED_CORR)
      is_eulerian = .false.
      if (present(vcoord)) is_eulerian = vcoord%coord_type == VCOORD_EULERIAN_Z
      if (present(bc)) then
         call continuity_gm_apply(grid, metrics, ct, ms, gm, dt, budget_w, tr_mode, bc=bc, &
                                  set_flux_h=is_eulerian)
      else
         call continuity_gm_apply(grid, metrics, ct, ms, gm, dt, budget_w, tr_mode, &
                                  set_flux_h=is_eulerian)
      end if
      if (is_eulerian) then
         ! Pin the Eulerian layers: the vertical w-divergence cancels the
         ! bolus divergence `continuity_gm_apply` left in `flux_h_layer`,
         ! carrying the tracers with it.  Column-local, so it runs BEFORE
         ! the ghost refresh below: a ghost column's `flux_h_layer` is the
         ! tile's own (the outermost ghost face carries no bolus), not its
         ! owner's, so a ghost advanced here and NOT refreshed afterwards
         ! made the step depend on where the seam fell (eulerian_z cells of
         ! the compatibility matrix's DECOMP leg).
         call compute_w_from_continuity(grid, va, ms)
         call tracer_advect_vertical(grid, va, ms, dt)
      end if
      call ocean_halo_exchange_ml_state(ms)
      if (present(bc)) then
         call ocean_periodic_wrap_state(grid, bc, ms, skip_x=ocean_halo_is_decomposed_x(), &
                                        skip_y=ocean_halo_is_decomposed_y())
         call ocean_fold_wrap_state(grid, bc, ms)
      end if
      if (dyn%check_h_positive) then
         call check_h_positive_or_die(grid, ms, "after the GM operator", 3, &
                                      dyn%outer_step_count + 1, check_layers=.true.)
      end if
      call profiler_stop("ocean_gm")
      call probe_dS(grid, ms, "after GM", 3, dyn%outer_step_count + 1)
   end subroutine run_gm_step