run_stage_split Subroutine

private subroutine run_stage_split(grid, metrics, dyn, eos, cor, ct, pgf, hv, bd, ss, va, hd, vd, vmix, ms, dt, n_inner, sf, geo, stage, vcoord, bc, sp, t, lateral_mix, epbl, kshear, mle, redi, varmix, vmix_tidal, meke, eta_forcing, td, cav, eta_pf_seam)

One FE stage of the split-explicit step. See the ocean_dyn_step_split header for the design.

Arguments

Type IntentOptional Attributes Name
type(hgrid_t), intent(in) :: grid
type(ocean_metrics_t), intent(in) :: metrics

Curvilinear horizontal metrics — forwarded to the PGF (and, in slice 2, the other geometry-aware kernels).
type(ocean_dyn_t), intent(inout) :: dyn
type(eos_t), intent(in) :: eos
type(coriolis_adv_t), intent(inout) :: cor
type(continuity_t), intent(inout) :: ct
type(ocean_pressure_force_t), intent(inout) :: pgf
type(ocean_horizontal_viscosity_t), intent(inout) :: hv
type(ocean_bottom_drag_t), intent(inout) :: bd
type(ocean_surface_stress_t), intent(inout) :: ss
type(ocean_vertical_advection_t), intent(inout) :: va
type(ocean_hdiff_tracer_t), intent(inout) :: hd
type(ocean_vdiff_t), intent(inout) :: vd
type(ocean_vmix_t), intent(inout) :: vmix
type(multilayer_state_t), intent(inout) :: ms
real(kind=wp), intent(in) :: dt
integer, intent(in) :: n_inner
type(ocean_surface_flux_t), intent(in), optional :: sf
type(ocean_geothermal_t), intent(in), optional :: geo

Geothermal bottom-heat-flux slot. See ocean_dyn_step_split.
integer, intent(in), optional :: stage

Outer SSP-RK2 stage index (1 or 2) — only used by the debug probe so trace lines self-identify which half-step they came from. Production paths leave it unset.
type(ocean_vcoord_t), intent(in), optional :: vcoord

Vertical-coordinate state. When present and coord_type /= VCOORD_EULERIAN_Z the stage runs in Lagrangian mode: vertical advection is skipped (no compute_w_from_continuity / tracer_advect_vertical) and apply_bt_correction’s h-rescale is skipped. The MOM6-constrained slow continuity makes sum_k(h_layer) self-consistent with the barotropic-substep η; the ALE remap at the end of the outer step relayers (h, hTr) onto vcoord%target_h conservatively. Absent or EULERIAN_Z ⇒ the historical Eulerian-z code path.
type(ocean_bc_state_t), intent(inout), optional :: bc

Open-boundary config forwarded straight to the barotropic substep. Absent / all-OBC_WALL keeps the closed-wall path.
type(ocean_sponge_t), intent(in), optional :: sp

Map-driven sponge slot (PR-23). See ocean_dyn_step_split.
real(kind=wp), intent(in), optional :: t

Wall-clock time for OBC_TIDAL constituent evaluation.
type(ocean_lateral_mix_t), intent(inout), optional :: lateral_mix

Flow-aware lateral closure. See the matching arg on ocean_dyn_step_split.
type(ocean_epbl_t), intent(inout), optional :: epbl

EPBL slot. See the matching arg on ocean_dyn_step_split.
type(ocean_kappa_shear_t), intent(inout), optional :: kshear

Kappa-shear slot. See ocean_dyn_step_split.
type(ocean_mle_t), intent(inout), optional :: mle

Fox-Kemper MLE slot (B5). Forwarded to continuity_tracer_step_split to fold the precomputed uhml/vhml into the mass fluxes. See ocean_dyn_step_split.
type(ocean_redi_t), intent(inout), optional :: redi

Redi neutral-diffusion slot (capability [3]). The Phase-A coefficients are precomputed in ocean_dyn_step_split; here Phase B (redi_apply_flux) adds the rotated tracer flux after the along-coordinate tracer_hdiff, at THERMO cadence.
type(ocean_varmix_t), intent(inout), optional :: varmix

VarMix coefficient slot (capability [4]). When enabled, its per-face khtr_u/khtr_v feed the Redi flux (the Visbeck KhTr seam); absent / disabled ⇒ Redi uses its scalar khtr.
type(ocean_tidal_mixing_t), intent(inout), optional :: vmix_tidal

Tidal-mixing slot. See ocean_dyn_step_split.
type(ocean_meke_t), intent(in), optional :: meke

MEKE slot (capability [5]). Read-only here: when meke%backscatter is on, its ku field is injected into the resolved per-face harmonic viscosity right after ocean_lateral_mix_compute (the energy-return seam, Gap 2). Absent / backscatter off ⇒ bit-identical.
real(kind=wp), intent(in), optional :: eta_forcing(grid%nx_total,grid%ny_total)
type(ocean_top_drag_t), intent(inout), optional :: td

Ice-shelf TOP-drag slot (&ocean_tdrag_nml). OPTIONAL so the many direct ocean_dyn_step* / run_stage* call sites in the test suite need no churn; the production driver always passes it. Absent, or present and disabled, => no kernel launch and a bit-identical step.
type(ocean_cavity_flux_t), intent(inout), optional :: cav

Ice-shelf basal-melt slot (&ocean_cavity_melt_nml). OPTIONAL for the same reason td is: the direct ocean_dyn_step* / run_stage* call sites in the test suite need no churn, and the production driver always passes it. Absent, disabled, or freshwater="virtual" => no kernel launch and a bit-identical step. It is threaded down here rather than acted on in engine_step_finalize because the real-freshwater volume must be spent in the SAME stage, at the SAME stage weight and from the SAME melt value as the salt and heat halves the surface-flux apply spends – see ocean_cavity_mass_step’s docstring.
real(kind=wp), intent(in), optional :: eta_pf_seam(grid%nx_total,grid%ny_total)

The part of the eta_forcing seam that ALSO reaches the slow PGF (eta_ib when &ocean_pgf_nml p_top_in_bc puts p_surf in the FV_MOM6 top BC). Absent ⇒ the slow PGF carries no seam load. Read only under &ocean_bt_nml bc_pgf_forcing. Equilibrium-tide elevation (C1), held static across the inner substep loop. Forwarded to barotropic_substep_nonlinear’s PGF; absent ⇒ bit-identical.

Calls

proc~~run_stage_split~~CallsGraph proc~run_stage_split run_stage_split interface~ocean_halo_centre ocean_halo_centre proc~run_stage_split->interface~ocean_halo_centre interface~ocean_halo_face_x ocean_halo_face_x proc~run_stage_split->interface~ocean_halo_face_x interface~ocean_halo_face_y ocean_halo_face_y proc~run_stage_split->interface~ocean_halo_face_y proc~accel_visc_rem_reweight accel_visc_rem_reweight proc~run_stage_split->proc~accel_visc_rem_reweight proc~accel_visc_rem_snapshot accel_visc_rem_snapshot proc~run_stage_split->proc~accel_visc_rem_snapshot proc~add_top_drag_into_f_slow add_top_drag_into_F_slow proc~run_stage_split->proc~add_top_drag_into_f_slow proc~apply_bt_correction apply_bt_correction proc~run_stage_split->proc~apply_bt_correction proc~apply_sw_and_restore apply_sw_and_restore proc~run_stage_split->proc~apply_sw_and_restore proc~barotropic_substep_nonlinear_interior barotropic_substep_nonlinear_interior proc~run_stage_split->proc~barotropic_substep_nonlinear_interior proc~bt_wide_copy_in bt_wide_t%bt_wide_copy_in proc~run_stage_split->proc~bt_wide_copy_in proc~bt_wide_copy_out bt_wide_t%bt_wide_copy_out proc~run_stage_split->proc~bt_wide_copy_out proc~bt_wide_entry_exchange bt_wide_t%bt_wide_entry_exchange proc~run_stage_split->proc~bt_wide_entry_exchange proc~bt_wide_substep bt_wide_substep proc~run_stage_split->proc~bt_wide_substep proc~chksum_bt chksum_bt proc~run_stage_split->proc~chksum_bt proc~chksum_hotface chksum_hotface proc~run_stage_split->proc~chksum_hotface proc~chksum_state chksum_state proc~run_stage_split->proc~chksum_state proc~compute_bt_rem compute_bt_rem proc~run_stage_split->proc~compute_bt_rem proc~compute_bt_rem_from_visc_rem compute_bt_rem_from_visc_rem proc~run_stage_split->proc~compute_bt_rem_from_visc_rem proc~compute_bt_rem_wave_drag compute_bt_rem_wave_drag proc~run_stage_split->proc~compute_bt_rem_wave_drag proc~compute_e_anom compute_e_anom proc~run_stage_split->proc~compute_e_anom proc~compute_gtot_faces compute_gtot_faces proc~run_stage_split->proc~compute_gtot_faces proc~compute_h_face_upstream compute_h_face_upstream proc~run_stage_split->proc~compute_h_face_upstream proc~compute_pbce compute_pbce proc~run_stage_split->proc~compute_pbce proc~coriolis_adv_apply_tendencies coriolis_adv_apply_tendencies proc~run_stage_split->proc~coriolis_adv_apply_tendencies proc~coriolis_adv_compute_tendencies coriolis_adv_compute_tendencies proc~run_stage_split->proc~coriolis_adv_compute_tendencies proc~derive_bt_from_layers derive_bt_from_layers proc~run_stage_split->proc~derive_bt_from_layers proc~face_depth_mean_rem_u face_depth_mean_rem_u proc~run_stage_split->proc~face_depth_mean_rem_u proc~face_depth_mean_rem_v face_depth_mean_rem_v proc~run_stage_split->proc~face_depth_mean_rem_v proc~face_depth_mean_u face_depth_mean_u proc~run_stage_split->proc~face_depth_mean_u proc~face_depth_mean_v face_depth_mean_v proc~run_stage_split->proc~face_depth_mean_v proc~isopycnal_vanish_tol isopycnal_vanish_tol proc~run_stage_split->proc~isopycnal_vanish_tol proc~ke_probe_sample ke_probe_sample proc~run_stage_split->proc~ke_probe_sample proc~lateral_mix_uses_resoln lateral_mix_uses_resoln proc~run_stage_split->proc~lateral_mix_uses_resoln proc~mask_bt_rem mask_bt_rem proc~run_stage_split->proc~mask_bt_rem proc~mask_layer_velocities mask_layer_velocities proc~run_stage_split->proc~mask_layer_velocities proc~meke_backscatter_apply meke_backscatter_apply proc~run_stage_split->proc~meke_backscatter_apply proc~ocean_bc_outer_face_tag ocean_bc_outer_face_tag proc~run_stage_split->proc~ocean_bc_outer_face_tag proc~ocean_bottom_drag_apply_tendencies ocean_bottom_drag_apply_tendencies proc~run_stage_split->proc~ocean_bottom_drag_apply_tendencies proc~ocean_bottom_drag_compute_tendencies ocean_bottom_drag_compute_tendencies proc~run_stage_split->proc~ocean_bottom_drag_compute_tendencies proc~ocean_cavity_mass_step ocean_cavity_mass_step proc~run_stage_split->proc~ocean_cavity_mass_step proc~ocean_channel_drag_apply_tendencies ocean_channel_drag_apply_tendencies proc~run_stage_split->proc~ocean_channel_drag_apply_tendencies proc~ocean_channel_drag_compute_tendencies ocean_channel_drag_compute_tendencies proc~run_stage_split->proc~ocean_channel_drag_compute_tendencies proc~ocean_dyn_is_thermo_step ocean_dyn_t%ocean_dyn_is_thermo_step proc~run_stage_split->proc~ocean_dyn_is_thermo_step proc~ocean_dyn_therm_dt ocean_dyn_t%ocean_dyn_therm_dt proc~run_stage_split->proc~ocean_dyn_therm_dt proc~ocean_eos_compute ocean_eos_compute proc~run_stage_split->proc~ocean_eos_compute proc~ocean_fold_wrap_state ocean_fold_wrap_state proc~run_stage_split->proc~ocean_fold_wrap_state proc~ocean_fold_wrap_time_means ocean_fold_wrap_time_means proc~run_stage_split->proc~ocean_fold_wrap_time_means proc~ocean_geothermal_apply_tracers ocean_geothermal_apply_tracers proc~run_stage_split->proc~ocean_geothermal_apply_tracers proc~ocean_halo_bt_group_2d ocean_halo_bt_group_2d proc~run_stage_split->proc~ocean_halo_bt_group_2d proc~ocean_halo_exchange_ml_state ocean_halo_exchange_ml_state proc~run_stage_split->proc~ocean_halo_exchange_ml_state proc~ocean_halo_is_decomposed_x ocean_halo_is_decomposed_x proc~run_stage_split->proc~ocean_halo_is_decomposed_x proc~ocean_halo_is_decomposed_y ocean_halo_is_decomposed_y proc~run_stage_split->proc~ocean_halo_is_decomposed_y proc~ocean_horizontal_viscosity_apply_tendencies ocean_horizontal_viscosity_apply_tendencies proc~run_stage_split->proc~ocean_horizontal_viscosity_apply_tendencies proc~ocean_horizontal_viscosity_compute_ke_diss ocean_horizontal_viscosity_compute_ke_diss proc~run_stage_split->proc~ocean_horizontal_viscosity_compute_ke_diss proc~ocean_horizontal_viscosity_compute_tendencies ocean_horizontal_viscosity_compute_tendencies proc~run_stage_split->proc~ocean_horizontal_viscosity_compute_tendencies proc~ocean_ideal_age_apply ocean_ideal_age_apply proc~run_stage_split->proc~ocean_ideal_age_apply proc~ocean_lateral_mix_compute ocean_lateral_mix_compute proc~run_stage_split->proc~ocean_lateral_mix_compute proc~ocean_obc_any_open_edge ocean_obc_any_open_edge proc~run_stage_split->proc~ocean_obc_any_open_edge proc~ocean_obc_apply_baroclinic ocean_obc_apply_baroclinic proc~run_stage_split->proc~ocean_obc_apply_baroclinic proc~ocean_periodic_wrap_centre_3d ocean_periodic_wrap_centre_3d proc~run_stage_split->proc~ocean_periodic_wrap_centre_3d proc~ocean_periodic_wrap_face_x_3d ocean_periodic_wrap_face_x_3d proc~run_stage_split->proc~ocean_periodic_wrap_face_x_3d proc~ocean_periodic_wrap_face_y_3d ocean_periodic_wrap_face_y_3d proc~run_stage_split->proc~ocean_periodic_wrap_face_y_3d proc~ocean_periodic_wrap_state ocean_periodic_wrap_state proc~run_stage_split->proc~ocean_periodic_wrap_state proc~ocean_pressure_force_apply ocean_pressure_force_apply proc~run_stage_split->proc~ocean_pressure_force_apply proc~ocean_pressure_force_compute ocean_pressure_force_compute proc~run_stage_split->proc~ocean_pressure_force_compute proc~ocean_sponge_apply ocean_sponge_apply proc~run_stage_split->proc~ocean_sponge_apply proc~ocean_sponge_apply_maps ocean_sponge_apply_maps proc~run_stage_split->proc~ocean_sponge_apply_maps proc~ocean_sponge_apply_tracers ocean_sponge_apply_tracers proc~run_stage_split->proc~ocean_sponge_apply_tracers proc~ocean_surface_flux_apply_tracers ocean_surface_flux_apply_tracers proc~run_stage_split->proc~ocean_surface_flux_apply_tracers proc~ocean_surface_stress_apply_tendencies ocean_surface_stress_apply_tendencies proc~run_stage_split->proc~ocean_surface_stress_apply_tendencies proc~ocean_surface_stress_compute_tendencies ocean_surface_stress_compute_tendencies proc~run_stage_split->proc~ocean_surface_stress_compute_tendencies proc~ocean_top_drag_apply_tendencies ocean_top_drag_apply_tendencies proc~run_stage_split->proc~ocean_top_drag_apply_tendencies proc~ocean_top_drag_compute_tendencies ocean_top_drag_compute_tendencies proc~run_stage_split->proc~ocean_top_drag_compute_tendencies proc~pgf_free_surface_gravity pgf_free_surface_gravity proc~run_stage_split->proc~pgf_free_surface_gravity proc~print_bt_budget print_bt_budget proc~run_stage_split->proc~print_bt_budget proc~probe_ds probe_dS proc~run_stage_split->proc~probe_ds proc~probe_h_vs_eta_residual probe_h_vs_eta_residual proc~run_stage_split->proc~probe_h_vs_eta_residual proc~profiler_start profiler_start proc~run_stage_split->proc~profiler_start proc~profiler_stop profiler_stop proc~run_stage_split->proc~profiler_stop proc~refresh_tracer_ghosts refresh_tracer_ghosts proc~run_stage_split->proc~refresh_tracer_ghosts proc~reset_bt_rem reset_bt_rem proc~run_stage_split->proc~reset_bt_rem proc~reset_vanished_layer_velocities reset_vanished_layer_velocities proc~run_stage_split->proc~reset_vanished_layer_velocities proc~run_continuity_chain run_continuity_chain proc~run_stage_split->proc~run_continuity_chain proc~set_cor_ref_velocity set_cor_ref_velocity proc~run_stage_split->proc~set_cor_ref_velocity proc~set_fast_forcing_eta_pf set_fast_forcing_eta_pf proc~run_stage_split->proc~set_fast_forcing_eta_pf proc~set_local_bt_cont_types set_local_BT_cont_types proc~run_stage_split->proc~set_local_bt_cont_types proc~snapshot_eta_pf snapshot_eta_PF proc~run_stage_split->proc~snapshot_eta_pf proc~subtract_fast_cor_ref subtract_fast_cor_ref proc~run_stage_split->proc~subtract_fast_cor_ref proc~sum_slow_tendencies_into_f_slow sum_slow_tendencies_into_F_slow proc~run_stage_split->proc~sum_slow_tendencies_into_f_slow proc~visc_rem_precompute visc_rem_precompute proc~run_stage_split->proc~visc_rem_precompute proc~vmix_apply_in_stage vmix_apply_in_stage proc~run_stage_split->proc~vmix_apply_in_stage 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~ocean_halo_face_x_2d ocean_halo_face_x_2d interface~ocean_halo_face_x->proc~ocean_halo_face_x_2d proc~ocean_halo_face_x_3d ocean_halo_face_x_3d interface~ocean_halo_face_x->proc~ocean_halo_face_x_3d proc~ocean_halo_face_y_2d ocean_halo_face_y_2d interface~ocean_halo_face_y->proc~ocean_halo_face_y_2d proc~ocean_halo_face_y_3d ocean_halo_face_y_3d interface~ocean_halo_face_y->proc~ocean_halo_face_y_3d local local proc~apply_bt_correction->local reduce reduce proc~apply_bt_correction->reduce proc~ocean_surface_flux_apply_sw_penetration ocean_surface_flux_apply_sw_penetration proc~apply_sw_and_restore->proc~ocean_surface_flux_apply_sw_penetration proc~ocean_surface_restore_apply_tracers ocean_surface_restore_apply_tracers proc~apply_sw_and_restore->proc~ocean_surface_restore_apply_tracers proc~barotropic_substep_nonlinear barotropic_substep_nonlinear proc~barotropic_substep_nonlinear_interior->proc~barotropic_substep_nonlinear proc~bt_wide_copy_in_impl bt_wide_copy_in_impl proc~bt_wide_copy_in->proc~bt_wide_copy_in_impl proc~bt_wide_copy_out_impl bt_wide_copy_out_impl proc~bt_wide_copy_out->proc~bt_wide_copy_out_impl proc~bt_wide_entry_exchange_impl bt_wide_entry_exchange_impl proc~bt_wide_entry_exchange->proc~bt_wide_entry_exchange_impl proc~bt_wide_substep->proc~barotropic_substep_nonlinear interface~rdb_debug_chksum rdb_debug_chksum proc~chksum_bt->interface~rdb_debug_chksum proc~chksum_active chksum_active proc~chksum_bt->proc~chksum_active proc~chksum_hotface->proc~chksum_active proc~chksum_argmax chksum_argmax proc~chksum_hotface->proc~chksum_argmax proc~hotface_row hotface_row proc~chksum_hotface->proc~hotface_row proc~chksum_state->interface~rdb_debug_chksum proc~chksum_state->proc~chksum_active proc~compute_bt_rem->local proc~bt_rem_open_impl bt_rem_open_impl proc~compute_bt_rem->proc~bt_rem_open_impl proc~compute_bt_rem_from_visc_rem->proc~face_depth_mean_u proc~compute_bt_rem_from_visc_rem->proc~face_depth_mean_v proc~compute_bt_rem_wave_drag->local proc~bt_rem_wave_drag_open_impl bt_rem_wave_drag_open_impl proc~compute_bt_rem_wave_drag->proc~bt_rem_wave_drag_open_impl proc~compute_gtot_faces->local proc~compute_h_face_upstream->local proc~h_face_upstream_open_impl h_face_upstream_open_impl proc~compute_h_face_upstream->proc~h_face_upstream_open_impl proc~compute_pbce->local proc~coriolis_adv_compute_tendencies_hk coriolis_adv_compute_tendencies_hk proc~coriolis_adv_compute_tendencies->proc~coriolis_adv_compute_tendencies_hk proc~coriolis_adv_compute_tendencies_sadourny coriolis_adv_compute_tendencies_sadourny proc~coriolis_adv_compute_tendencies->proc~coriolis_adv_compute_tendencies_sadourny proc~coriolis_adv_compute_tendencies_sadourny_energy coriolis_adv_compute_tendencies_sadourny_energy proc~coriolis_adv_compute_tendencies->proc~coriolis_adv_compute_tendencies_sadourny_energy frhat_h_face_step frhat_h_face_step proc~derive_bt_from_layers->frhat_h_face_step proc~derive_bt_from_layers->local proc~face_depth_mean_rem_u->frhat_h_face_step proc~face_depth_mean_rem_u->local proc~face_depth_mean_rem_v->frhat_h_face_step proc~face_depth_mean_rem_v->local proc~face_depth_mean_u->frhat_h_face_step proc~face_depth_mean_u->local proc~face_depth_mean_v->frhat_h_face_step proc~face_depth_mean_v->local proc~ke_probe_active ke_probe_active proc~ke_probe_sample->proc~ke_probe_active proc~ke_regions_impl ke_regions_impl proc~ke_probe_sample->proc~ke_regions_impl proc~meke_backscatter_apply_impl meke_backscatter_apply_impl proc~meke_backscatter_apply->proc~meke_backscatter_apply_impl proc~ocean_bottom_drag_compute_tendencies->local proc~compute_distributed_drag compute_distributed_drag proc~ocean_bottom_drag_compute_tendencies->proc~compute_distributed_drag error error proc~ocean_cavity_mass_step->error proc~cavity_comp_apply_impl cavity_comp_apply_impl proc~ocean_cavity_mass_step->proc~cavity_comp_apply_impl proc~cavity_comp_scale_tracer_impl cavity_comp_scale_tracer_impl proc~ocean_cavity_mass_step->proc~cavity_comp_scale_tracer_impl proc~cavity_comp_withdrawal cavity_comp_withdrawal proc~ocean_cavity_mass_step->proc~cavity_comp_withdrawal proc~cavity_mass_apply_impl cavity_mass_apply_impl proc~ocean_cavity_mass_step->proc~cavity_mass_apply_impl proc~cavity_mass_salt_mirror_impl cavity_mass_salt_mirror_impl proc~ocean_cavity_mass_step->proc~cavity_mass_salt_mirror_impl proc~cavity_mass_thin_is_fatal cavity_mass_thin_is_fatal proc~ocean_cavity_mass_step->proc~cavity_mass_thin_is_fatal proc~cavity_mass_totals_impl cavity_mass_totals_impl proc~ocean_cavity_mass_step->proc~cavity_mass_totals_impl proc~fail fail proc~ocean_cavity_mass_step->proc~fail proc~halo_allreduce_sum halo_allreduce_sum proc~ocean_cavity_mass_step->proc~halo_allreduce_sum to_string to_string proc~ocean_cavity_mass_step->to_string proc~compute_channel_drag_rates compute_channel_drag_rates proc~ocean_channel_drag_compute_tendencies->proc~compute_channel_drag_rates proc~eos_compute_arrays eos_compute_arrays proc~ocean_eos_compute->proc~eos_compute_arrays interface~fold_north_centre fold_north_centre proc~ocean_fold_wrap_state->interface~fold_north_centre interface~fold_north_u_face fold_north_u_face proc~ocean_fold_wrap_state->interface~fold_north_u_face interface~fold_north_v_face fold_north_v_face proc~ocean_fold_wrap_state->interface~fold_north_v_face interface~ocean_fold_pack ocean_fold_pack proc~ocean_fold_wrap_state->interface~ocean_fold_pack interface~ocean_fold_unpack ocean_fold_unpack proc~ocean_fold_wrap_state->interface~ocean_fold_unpack proc~n_tracers n_tracers proc~ocean_fold_wrap_state->proc~n_tracers proc~ocean_fold_begin ocean_fold_begin proc~ocean_fold_wrap_state->proc~ocean_fold_begin proc~ocean_fold_end ocean_fold_end proc~ocean_fold_wrap_state->proc~ocean_fold_end proc~ocean_fold_exchange ocean_fold_exchange proc~ocean_fold_wrap_state->proc~ocean_fold_exchange proc~ocean_fold_is_distributed ocean_fold_is_distributed proc~ocean_fold_wrap_state->proc~ocean_fold_is_distributed proc~ocean_fold_wrap_time_means->interface~fold_north_centre proc~ocean_fold_wrap_time_means->interface~fold_north_u_face proc~ocean_fold_wrap_time_means->interface~fold_north_v_face proc~ocean_fold_wrap_time_means->interface~ocean_fold_pack proc~ocean_fold_wrap_time_means->interface~ocean_fold_unpack proc~ocean_fold_wrap_time_means->proc~ocean_fold_begin proc~ocean_fold_wrap_time_means->proc~ocean_fold_end proc~ocean_fold_wrap_time_means->proc~ocean_fold_exchange proc~ocean_fold_wrap_time_means->proc~ocean_fold_is_distributed proc~apply_geothermal_src_impl apply_geothermal_src_impl proc~ocean_geothermal_apply_tracers->proc~apply_geothermal_src_impl proc~ocean_halo_bt_group_2d->proc~ocean_halo_centre_2d proc~ocean_halo_bt_group_2d->proc~ocean_halo_face_x_2d proc~ocean_halo_bt_group_2d->proc~ocean_halo_face_y_2d proc~oh_count_bt_group oh_count_bt_group proc~ocean_halo_bt_group_2d->proc~oh_count_bt_group proc~oh_count_suppress_off oh_count_suppress_off proc~ocean_halo_bt_group_2d->proc~oh_count_suppress_off proc~oh_count_suppress_on oh_count_suppress_on proc~ocean_halo_bt_group_2d->proc~oh_count_suppress_on proc~ocean_halo_exchange_ml_state->interface~ocean_halo_centre proc~ocean_halo_exchange_ml_state->interface~ocean_halo_face_x proc~ocean_halo_exchange_ml_state->interface~ocean_halo_face_y proc~ocean_halo_exchange_ml_state->proc~profiler_start proc~ocean_halo_exchange_ml_state->proc~profiler_stop proc~oh_count_ml_state oh_count_ml_state proc~ocean_halo_exchange_ml_state->proc~oh_count_ml_state proc~ocean_halo_exchange_ml_state->proc~oh_count_suppress_off proc~ocean_halo_exchange_ml_state->proc~oh_count_suppress_on proc~hvisc_apply_impl hvisc_apply_impl proc~ocean_horizontal_viscosity_apply_tendencies->proc~hvisc_apply_impl proc~hvisc_ke_diss_impl hvisc_ke_diss_impl proc~ocean_horizontal_viscosity_compute_ke_diss->proc~hvisc_ke_diss_impl proc~ocean_horizontal_viscosity_compute_tendencies_on ocean_horizontal_viscosity_compute_tendencies_on proc~ocean_horizontal_viscosity_compute_tendencies->proc~ocean_horizontal_viscosity_compute_tendencies_on proc~ocean_ideal_age_age_step ocean_ideal_age_age_step proc~ocean_ideal_age_apply->proc~ocean_ideal_age_age_step proc~ocean_lateral_mix_compute_leith ocean_lateral_mix_compute_leith proc~ocean_lateral_mix_compute->proc~ocean_lateral_mix_compute_leith proc~ocean_lateral_mix_compute_leith_biharm ocean_lateral_mix_compute_leith_biharm proc~ocean_lateral_mix_compute->proc~ocean_lateral_mix_compute_leith_biharm proc~ocean_lateral_mix_compute_smag ocean_lateral_mix_compute_smag proc~ocean_lateral_mix_compute->proc~ocean_lateral_mix_compute_smag proc~ocean_lateral_mix_compute_smag_ah ocean_lateral_mix_compute_smag_ah proc~ocean_lateral_mix_compute->proc~ocean_lateral_mix_compute_smag_ah proc~ocean_lateral_mix_compute_vel_scale ocean_lateral_mix_compute_vel_scale proc~ocean_lateral_mix_compute->proc~ocean_lateral_mix_compute_vel_scale proc~is_open_ish is_open_ish proc~ocean_obc_any_open_edge->proc~is_open_ish proc~apply_clamped_meridional_north apply_clamped_meridional_north proc~ocean_obc_apply_baroclinic->proc~apply_clamped_meridional_north proc~apply_clamped_meridional_south apply_clamped_meridional_south proc~ocean_obc_apply_baroclinic->proc~apply_clamped_meridional_south proc~apply_clamped_zonal_east apply_clamped_zonal_east proc~ocean_obc_apply_baroclinic->proc~apply_clamped_zonal_east proc~apply_clamped_zonal_west apply_clamped_zonal_west proc~ocean_obc_apply_baroclinic->proc~apply_clamped_zonal_west proc~apply_meridional_baroclinic apply_meridional_baroclinic proc~ocean_obc_apply_baroclinic->proc~apply_meridional_baroclinic proc~apply_nudge_meridional_north apply_nudge_meridional_north proc~ocean_obc_apply_baroclinic->proc~apply_nudge_meridional_north proc~apply_nudge_meridional_south apply_nudge_meridional_south proc~ocean_obc_apply_baroclinic->proc~apply_nudge_meridional_south proc~apply_nudge_zonal_east apply_nudge_zonal_east proc~ocean_obc_apply_baroclinic->proc~apply_nudge_zonal_east proc~apply_nudge_zonal_west apply_nudge_zonal_west proc~ocean_obc_apply_baroclinic->proc~apply_nudge_zonal_west proc~apply_orlanski_east apply_orlanski_east proc~ocean_obc_apply_baroclinic->proc~apply_orlanski_east proc~apply_orlanski_north apply_orlanski_north proc~ocean_obc_apply_baroclinic->proc~apply_orlanski_north proc~apply_orlanski_south apply_orlanski_south proc~ocean_obc_apply_baroclinic->proc~apply_orlanski_south proc~apply_orlanski_west apply_orlanski_west proc~ocean_obc_apply_baroclinic->proc~apply_orlanski_west proc~apply_zonal_baroclinic apply_zonal_baroclinic proc~ocean_obc_apply_baroclinic->proc~apply_zonal_baroclinic proc~fill_uv_layer_ghosts fill_uv_layer_ghosts proc~ocean_obc_apply_baroclinic->proc~fill_uv_layer_ghosts proc~ocean_obc_apply_baroclinic->proc~is_open_ish proc~is_radiating is_radiating proc~ocean_obc_apply_baroclinic->proc~is_radiating proc~snapshot_u_prev_east snapshot_u_prev_east proc~ocean_obc_apply_baroclinic->proc~snapshot_u_prev_east proc~snapshot_u_prev_west snapshot_u_prev_west proc~ocean_obc_apply_baroclinic->proc~snapshot_u_prev_west proc~snapshot_v_prev_north snapshot_v_prev_north proc~ocean_obc_apply_baroclinic->proc~snapshot_v_prev_north proc~snapshot_v_prev_south snapshot_v_prev_south proc~ocean_obc_apply_baroclinic->proc~snapshot_v_prev_south proc~ocean_periodic_wrap_state->proc~ocean_periodic_wrap_centre_3d proc~ocean_periodic_wrap_state->proc~ocean_periodic_wrap_face_x_3d proc~ocean_periodic_wrap_state->proc~ocean_periodic_wrap_face_y_3d proc~ocean_pressure_force_compute->local proc~compute_fv_mom6_impl compute_fv_mom6_impl proc~ocean_pressure_force_compute->proc~compute_fv_mom6_impl proc~compute_fv_mom6_insitu_pcm_impl compute_fv_mom6_insitu_pcm_impl proc~ocean_pressure_force_compute->proc~compute_fv_mom6_insitu_pcm_impl proc~compute_fv_mom6_reconstruct_impl compute_fv_mom6_reconstruct_impl proc~ocean_pressure_force_compute->proc~compute_fv_mom6_reconstruct_impl proc~compute_gprime_impl compute_gprime_impl proc~ocean_pressure_force_compute->proc~compute_gprime_impl proc~eos_wright_pgf_column_sweep_impl eos_wright_pgf_column_sweep_impl proc~ocean_pressure_force_compute->proc~eos_wright_pgf_column_sweep_impl proc~use_insitu_pcm use_insitu_pcm proc~ocean_pressure_force_compute->proc~use_insitu_pcm proc~ocean_sponge_apply->local proc~relax_map_tracer_budget_impl relax_map_tracer_budget_impl proc~ocean_sponge_apply_maps->proc~relax_map_tracer_budget_impl proc~relax_map_tracer_impl relax_map_tracer_impl proc~ocean_sponge_apply_maps->proc~relax_map_tracer_impl proc~relax_map_u_impl relax_map_u_impl proc~ocean_sponge_apply_maps->proc~relax_map_u_impl proc~relax_map_v_impl relax_map_v_impl proc~ocean_sponge_apply_maps->proc~relax_map_v_impl proc~sponge_relax_band_x_tracer sponge_relax_band_x_tracer proc~ocean_sponge_apply_tracers->proc~sponge_relax_band_x_tracer proc~sponge_relax_band_y_tracer sponge_relax_band_y_tracer proc~ocean_sponge_apply_tracers->proc~sponge_relax_band_y_tracer proc~apply_surface_src_2d_dyn_impl apply_surface_src_2d_dyn_impl proc~ocean_surface_flux_apply_tracers->proc~apply_surface_src_2d_dyn_impl proc~apply_surface_src_2d_dyn_nobudget_impl apply_surface_src_2d_dyn_nobudget_impl proc~ocean_surface_flux_apply_tracers->proc~apply_surface_src_2d_dyn_nobudget_impl proc~apply_surface_src_2d_impl apply_surface_src_2d_impl proc~ocean_surface_flux_apply_tracers->proc~apply_surface_src_2d_impl proc~apply_surface_src_2d_nobudget_impl apply_surface_src_2d_nobudget_impl proc~ocean_surface_flux_apply_tracers->proc~apply_surface_src_2d_nobudget_impl proc~surfstress_apply_impl surfstress_apply_impl proc~ocean_surface_stress_apply_tendencies->proc~surfstress_apply_impl proc~surfstress_compute_impl surfstress_compute_impl proc~ocean_surface_stress_compute_tendencies->proc~surfstress_compute_impl proc~surfstress_distributed_impl surfstress_distributed_impl proc~ocean_surface_stress_compute_tendencies->proc~surfstress_distributed_impl proc~top_drag_stress_mag_impl top_drag_stress_mag_impl proc~ocean_top_drag_compute_tendencies->proc~top_drag_stress_mag_impl proc~top_drag_tendencies_impl top_drag_tendencies_impl proc~ocean_top_drag_compute_tendencies->proc~top_drag_tendencies_impl proc~emit_row emit_row proc~print_bt_budget->proc~emit_row proc~region_eta_uv region_eta_uv proc~print_bt_budget->proc~region_eta_uv proc~region_power region_power proc~print_bt_budget->proc~region_power 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~refresh_tracer_ghosts->interface~ocean_halo_centre proc~refresh_tracer_ghosts->proc~profiler_start proc~refresh_tracer_ghosts->proc~profiler_stop proc~ocean_fold_wrap_centre_3d_state ocean_fold_wrap_centre_3d_state proc~refresh_tracer_ghosts->proc~ocean_fold_wrap_centre_3d_state proc~run_continuity_chain->interface~ocean_halo_centre proc~run_continuity_chain->proc~ocean_dyn_is_thermo_step proc~run_continuity_chain->proc~ocean_fold_wrap_state proc~run_continuity_chain->proc~ocean_halo_exchange_ml_state proc~run_continuity_chain->proc~ocean_halo_is_decomposed_x proc~run_continuity_chain->proc~ocean_halo_is_decomposed_y proc~run_continuity_chain->proc~ocean_obc_any_open_edge proc~run_continuity_chain->proc~ocean_periodic_wrap_state proc~run_continuity_chain->proc~probe_ds proc~run_continuity_chain->proc~profiler_start proc~run_continuity_chain->proc~profiler_stop proc~run_continuity_chain->proc~refresh_tracer_ghosts proc~check_h_positive_or_die check_h_positive_or_die proc~run_continuity_chain->proc~check_h_positive_or_die proc~compute_w_from_continuity compute_w_from_continuity proc~run_continuity_chain->proc~compute_w_from_continuity proc~continuity_tracer_step_split continuity_tracer_step_split proc~run_continuity_chain->proc~continuity_tracer_step_split proc~copy_field_3d copy_field_3d proc~run_continuity_chain->proc~copy_field_3d proc~ocean_accumulate_mass_out ocean_accumulate_mass_out proc~run_continuity_chain->proc~ocean_accumulate_mass_out proc~ocean_apply_conservative_min_thickness ocean_apply_conservative_min_thickness proc~run_continuity_chain->proc~ocean_apply_conservative_min_thickness proc~ocean_obc_fill_ghosts ocean_obc_fill_ghosts proc~run_continuity_chain->proc~ocean_obc_fill_ghosts proc~ocean_obc_update_reservoirs ocean_obc_update_reservoirs proc~run_continuity_chain->proc~ocean_obc_update_reservoirs proc~redi_apply_flux redi_apply_flux proc~run_continuity_chain->proc~redi_apply_flux proc~rk2_average_field_3d rk2_average_field_3d proc~run_continuity_chain->proc~rk2_average_field_3d proc~tracer_advect_vertical tracer_advect_vertical proc~run_continuity_chain->proc~tracer_advect_vertical proc~tracer_hdiff tracer_hdiff proc~run_continuity_chain->proc~tracer_hdiff proc~set_cor_ref_velocity->proc~face_depth_mean_rem_u proc~set_cor_ref_velocity->proc~face_depth_mean_rem_v proc~set_cor_ref_velocity->proc~face_depth_mean_u proc~set_cor_ref_velocity->proc~face_depth_mean_v proc~set_fast_forcing_eta_pf->local proc~set_local_bt_cont_types->local proc~subtract_fast_cor_ref->local proc~vdiff_apply_momentum vdiff_apply_momentum proc~visc_rem_precompute->proc~vdiff_apply_momentum proc~visc_rem_halo_refresh visc_rem_halo_refresh proc~visc_rem_precompute->proc~visc_rem_halo_refresh proc~vmix_apply_in_stage->proc~ocean_dyn_is_thermo_step proc~vmix_apply_in_stage->proc~ocean_dyn_therm_dt proc~vmix_apply_in_stage->proc~profiler_start proc~vmix_apply_in_stage->proc~profiler_stop proc~vmix_apply_in_stage->proc~visc_rem_precompute proc~vmix_apply_in_stage->error proc~epbl_compute epbl_compute proc~vmix_apply_in_stage->proc~epbl_compute proc~epbl_merge_into_kv_kt epbl_merge_into_kv_kt proc~vmix_apply_in_stage->proc~epbl_merge_into_kv_kt proc~kappa_shear_compute kappa_shear_compute proc~vmix_apply_in_stage->proc~kappa_shear_compute proc~kappa_shear_merge_into_kv_kt kappa_shear_merge_into_kv_kt proc~vmix_apply_in_stage->proc~kappa_shear_merge_into_kv_kt proc~tidal_mixing_compute tidal_mixing_compute proc~vmix_apply_in_stage->proc~tidal_mixing_compute proc~tidal_mixing_merge_into_kt tidal_mixing_merge_into_kt proc~vmix_apply_in_stage->proc~tidal_mixing_merge_into_kt proc~vmix_apply_in_stage->proc~vdiff_apply_momentum proc~vdiff_apply_tracers vdiff_apply_tracers proc~vmix_apply_in_stage->proc~vdiff_apply_tracers proc~vmix_apply_in_stage->proc~visc_rem_halo_refresh proc~vmix_add_kv_ml_invz2 vmix_add_kv_ml_invz2 proc~vmix_apply_in_stage->proc~vmix_add_kv_ml_invz2 proc~vmix_apply_convection vmix_apply_convection proc~vmix_apply_in_stage->proc~vmix_apply_convection proc~vmix_apply_kpp_overlay vmix_apply_kpp_overlay proc~vmix_apply_in_stage->proc~vmix_apply_kpp_overlay proc~vmix_apply_nonlocal_tendencies vmix_apply_nonlocal_tendencies proc~vmix_apply_in_stage->proc~vmix_apply_nonlocal_tendencies proc~vmix_assemble vmix_assemble proc~vmix_apply_in_stage->proc~vmix_assemble proc~vmix_compute_pp81 vmix_compute_pp81 proc~vmix_apply_in_stage->proc~vmix_compute_pp81 proc~vmix_split_kd_heat_salt vmix_split_kd_heat_salt proc~vmix_apply_in_stage->proc~vmix_split_kd_heat_salt

Called by

proc~~run_stage_split~~CalledByGraph proc~run_stage_split run_stage_split proc~ocean_dyn_step_split ocean_dyn_step_split proc~ocean_dyn_step_split->proc~run_stage_split 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
integer, private :: bc_e_drv

Per-edge BC tags used for the slow-path transport wall-zeroing. Default OBC_WALL; overridden from bc%%bc_type when bc is present.
integer, private :: bc_n_drv

Per-edge BC tags used for the slow-path transport wall-zeroing. Default OBC_WALL; overridden from bc%%bc_type when bc is present.
integer, private :: bc_s_drv

Per-edge BC tags used for the slow-path transport wall-zeroing. Default OBC_WALL; overridden from bc%%bc_type when bc is present.
integer, private :: bc_w_drv

Per-edge BC tags used for the slow-path transport wall-zeroing. Default OBC_WALL; overridden from bc%%bc_type when bc is present.
real(kind=wp), private :: chain_weight
real(kind=wp), private :: dt_inner
real(kind=wp), private :: dt_vel
logical, private :: fold_top

.true. when the ice-shelf top drag is folded into the vdiff k = nz diagonal. Tested instead of present(td) because a DISABLED top-drag slot carries placeholder-sized arrays, and the explicit-shape dummy they would reach in diffuse_velocity_columns_impl is device-mapped unconditionally — see the note in run_stage.
real(kind=wp), private :: h_min_floor
logical, private :: has_e_drv

Physical-edge flags for the transport wall-zeroing. .true. = physical domain edge (WALL zero applies); .false. = MPI seam (halo owns it, O0).
logical, private :: has_n_drv

Physical-edge flags for the transport wall-zeroing. .true. = physical domain edge (WALL zero applies); .false. = MPI seam (halo owns it, O0).
logical, private :: has_s_drv

Physical-edge flags for the transport wall-zeroing. .true. = physical domain edge (WALL zero applies); .false. = MPI seam (halo owns it, O0).
logical, private :: has_w_drv

Physical-edge flags for the transport wall-zeroing. .true. = physical domain edge (WALL zero applies); .false. = MPI seam (halo owns it, O0).
integer, private :: i
integer, private :: i_ss

Loop/extent locals for the inline stress_shelf publish.
logical, private :: is_lagrangian
logical, private :: is_pc

pred_corr stage roles (SPEC §2/§4 S4): is_pred = the predictor (stage 1 under pred_corr) — provisional up = u + BE·dt·accel, h-only continuity into a discarded hp, coefficients/remnant only, no tracer physics. Stage 2 is the corrector: the single full-dt prognostic update. Both .false. under ssp_rk2 ⇒ every gate below is untaken ⇒ bit-identical.
logical, private :: is_pred

pred_corr stage roles (SPEC §2/§4 S4): is_pred = the predictor (stage 1 under pred_corr) — provisional up = u + BE·dt·accel, h-only continuity into a discarded hp, coefficients/remnant only, no tracer physics. Stage 2 is the corrector: the single full-dt prognostic update. Both .false. under ssp_rk2 ⇒ every gate below is untaken ⇒ bit-identical.
integer, private :: j
integer, private :: j_ss

Loop/extent locals for the inline stress_shelf publish.
integer, private :: nx_face
integer, private :: nx_ss

Loop/extent locals for the inline stress_shelf publish.
integer, private :: nx_vface
integer, private :: ny_face
integer, private :: ny_ss

Loop/extent locals for the inline stress_shelf publish.
integer, private :: ny_uface
logical, private :: publish_shelf

.true. when the ice-shelf top-drag slot is live and its stress_top is therefore full-sized and freshly written.
logical, private :: sponge_maps_on
logical, private :: sponge_seam

.true. when the map-driven sponge (PR-23) supersedes the legacy band path. A local logical because Fortran does not guarantee .and. short-circuits past present().
integer, private :: stage_id
integer, private :: step_id
logical, private :: therm_active
real(kind=wp), private :: therm_dt

Source Code

   subroutine run_stage_split(grid, metrics, dyn, eos, cor, ct, pgf, hv, bd, ss, &
                              va, hd, vd, vmix, ms, dt, n_inner, sf, geo, stage, vcoord, bc, sp, t, &
                              lateral_mix, epbl, kshear, mle, redi, varmix, vmix_tidal, meke, &
                              eta_forcing, td, cav, eta_pf_seam)
      !! One FE stage of the split-explicit step.  See the
      !! `ocean_dyn_step_split` header for the design.
      type(hgrid_t), intent(in) :: grid
      type(ocean_metrics_t), intent(in) :: metrics
         !! Curvilinear horizontal metrics — forwarded to the PGF (and,
         !! in slice 2, the other geometry-aware kernels).
      type(ocean_dyn_t), intent(inout) :: dyn
      type(eos_t), intent(in) :: eos
      type(coriolis_adv_t), intent(inout) :: cor
      type(continuity_t), intent(inout) :: ct
      type(ocean_pressure_force_t), intent(inout) :: pgf
      type(ocean_horizontal_viscosity_t), intent(inout) :: hv
      type(ocean_bottom_drag_t), intent(inout) :: bd
      type(ocean_top_drag_t), intent(inout), optional :: td
         !! Ice-shelf TOP-drag slot (`&ocean_tdrag_nml`).  OPTIONAL so the
         !! many direct `ocean_dyn_step*` / `run_stage*` call sites in the
         !! test suite need no churn; the production driver always passes
         !! it.  Absent, or present and disabled, => no kernel launch and
         !! a bit-identical step.
      type(ocean_cavity_flux_t), intent(inout), optional :: cav
         !! Ice-shelf basal-melt slot (`&ocean_cavity_melt_nml`).
         !! OPTIONAL for the same reason `td` is: the direct
         !! `ocean_dyn_step*` / `run_stage*` call sites in the test suite
         !! need no churn, and the production driver always passes it.
         !! Absent, disabled, or `freshwater="virtual"` => no kernel
         !! launch and a bit-identical step.  It is threaded down here
         !! rather than acted on in `engine_step_finalize` because the
         !! real-freshwater volume must be spent in the SAME stage, at
         !! the SAME stage weight and from the SAME `melt` value as the
         !! salt and heat halves the surface-flux apply spends -- see
         !! `ocean_cavity_mass_step`'s docstring.
      type(ocean_surface_stress_t), intent(inout) :: ss
      type(ocean_vertical_advection_t), intent(inout) :: va
      type(ocean_hdiff_tracer_t), intent(inout) :: hd
      type(ocean_vdiff_t), intent(inout) :: vd
      type(ocean_vmix_t), intent(inout) :: vmix
      type(multilayer_state_t), intent(inout) :: ms
      real(wp), intent(in) :: dt
      integer, intent(in) :: n_inner
      type(ocean_surface_flux_t), intent(in), optional :: sf
      type(ocean_geothermal_t), intent(in), optional :: geo
         !! Geothermal bottom-heat-flux slot.  See `ocean_dyn_step_split`.
      integer, intent(in), optional :: stage
         !! Outer SSP-RK2 stage index (1 or 2) — only used by the
         !! debug probe so trace lines self-identify which half-step
         !! they came from.  Production paths leave it unset.
      type(ocean_vcoord_t), intent(in), optional :: vcoord
         !! Vertical-coordinate state.  When present and
         !! `coord_type /= VCOORD_EULERIAN_Z` the stage runs in
         !! Lagrangian mode: vertical advection is skipped (no
         !! `compute_w_from_continuity` / `tracer_advect_vertical`)
         !! and `apply_bt_correction`'s h-rescale is skipped.  The
         !! MOM6-constrained slow continuity makes `sum_k(h_layer)`
         !! self-consistent with the barotropic-substep η; the ALE remap at
         !! the end of the outer step relayers (h, hTr) onto
         !! `vcoord%target_h` conservatively.  Absent or EULERIAN_Z
         !! ⇒ the historical Eulerian-z code path.
      type(ocean_bc_state_t), intent(inout), optional :: bc
         !! Open-boundary config forwarded straight to the barotropic substep.
         !! Absent / all-OBC_WALL keeps the closed-wall path.
      type(ocean_sponge_t), intent(in), optional :: sp
         !! Map-driven sponge slot (PR-23).  See `ocean_dyn_step_split`.
      real(wp), intent(in), optional :: t
         !! Wall-clock time for OBC_TIDAL constituent evaluation.
      type(ocean_lateral_mix_t), intent(inout), optional :: lateral_mix
         !! Flow-aware lateral closure.  See the matching arg on
         !! `ocean_dyn_step_split`.
      type(ocean_epbl_t), intent(inout), optional :: epbl
         !! EPBL slot.  See the matching arg on `ocean_dyn_step_split`.
      type(ocean_kappa_shear_t), intent(inout), optional :: kshear
         !! Kappa-shear slot.  See `ocean_dyn_step_split`.
      type(ocean_tidal_mixing_t), intent(inout), optional :: vmix_tidal
         !! Tidal-mixing slot.  See `ocean_dyn_step_split`.
      type(ocean_mle_t), intent(inout), optional :: mle
         !! Fox-Kemper MLE slot (B5).  Forwarded to
         !! `continuity_tracer_step_split` to fold the precomputed
         !! `uhml`/`vhml` into the mass fluxes.  See `ocean_dyn_step_split`.
      type(ocean_redi_t), intent(inout), optional :: redi
         !! Redi neutral-diffusion slot (capability [3]).  The Phase-A
         !! coefficients are precomputed in `ocean_dyn_step_split`; here
         !! Phase B (`redi_apply_flux`) adds the rotated tracer flux after
         !! the along-coordinate `tracer_hdiff`, at THERMO cadence.
      type(ocean_varmix_t), intent(inout), optional :: varmix
         !! VarMix coefficient slot (capability [4]).  When enabled, its
         !! per-face `khtr_u`/`khtr_v` feed the Redi flux (the Visbeck KhTr
         !! seam); absent / disabled ⇒ Redi uses its scalar `khtr`.
      type(ocean_meke_t), intent(in), optional :: meke
         !! MEKE slot (capability [5]).  Read-only here: when
         !! `meke%backscatter` is on, its `ku` field is injected into the
         !! resolved per-face harmonic viscosity right after
         !! `ocean_lateral_mix_compute` (the energy-return seam, Gap 2).
         !! Absent / `backscatter` off ⇒ bit-identical.
      real(wp), intent(in), optional :: eta_forcing(grid%nx_total, grid%ny_total)
      real(wp), intent(in), optional :: eta_pf_seam(grid%nx_total, grid%ny_total)
         !! The part of the `eta_forcing` seam that ALSO reaches the slow PGF
         !! (`eta_ib` when `&ocean_pgf_nml p_top_in_bc` puts `p_surf` in the
         !! FV_MOM6 top BC).  Absent ⇒ the slow PGF carries no seam load.
         !! Read only under `&ocean_bt_nml bc_pgf_forcing`.
         !! Equilibrium-tide elevation (C1), held static across the inner
         !! substep loop.  Forwarded to `barotropic_substep_nonlinear`'s
         !! PGF; absent ⇒ bit-identical.

      integer :: i, j, nx_face, ny_uface, nx_vface, ny_face
      integer :: stage_id, step_id
      integer :: bc_w_drv, bc_e_drv, bc_s_drv, bc_n_drv
         !! Per-edge BC tags used for the slow-path transport wall-zeroing.
         !! Default OBC_WALL; overridden from bc%<edge>%bc_type when bc is present.
      logical :: has_w_drv, has_e_drv, has_s_drv, has_n_drv
         !! Physical-edge flags for the transport wall-zeroing.  .true. = physical
         !! domain edge (WALL zero applies); .false. = MPI seam (halo owns it, O0).
      logical :: is_lagrangian, therm_active
      logical :: is_pc, is_pred
         !! pred_corr stage roles (SPEC §2/§4 S4): `is_pred` = the predictor
         !! (stage 1 under pred_corr) — provisional `up = u + BE·dt·accel`,
         !! h-only continuity into a discarded hp, coefficients/remnant
         !! only, no tracer physics.  Stage 2 is the corrector: the single
         !! full-dt prognostic update.  Both .false. under ssp_rk2 ⇒ every
         !! gate below is untaken ⇒ bit-identical.
      logical :: sponge_maps_on
      logical :: sponge_seam
         !! .true. when the map-driven sponge (PR-23) supersedes the legacy
         !! band path. A local logical because Fortran does not guarantee
         !! `.and.` short-circuits past `present()`.
      logical :: fold_top
         !! `.true.` when the ice-shelf top drag is folded into the vdiff
         !! `k = nz` diagonal.  Tested instead of `present(td)` because a
         !! DISABLED top-drag slot carries placeholder-sized arrays, and
         !! the explicit-shape dummy they would reach in
         !! `diffuse_velocity_columns_impl` is device-mapped
         !! unconditionally — see the note in `run_stage`.
      logical :: publish_shelf
         !! `.true.` when the ice-shelf top-drag slot is live and its
         !! `stress_top` is therefore full-sized and freshly written.
      integer :: i_ss, j_ss, nx_ss, ny_ss
         !! Loop/extent locals for the inline `stress_shelf` publish.
      real(wp) :: dt_inner, therm_dt
      real(wp) :: h_min_floor
      real(wp) :: chain_weight, dt_vel

      fold_top = .false.
      if (present(td)) fold_top = td%implicit_fold
      stage_id = 0
      if (present(stage)) stage_id = stage
      ! Probes label the OUTER step we're INSIDE — i.e. the one
      ! `outer_step_count` increments to at the end of this call.
      step_id = dyn%outer_step_count + 1
      is_lagrangian = .false.
      if (present(vcoord)) is_lagrangian = vcoord%coord_type /= VCOORD_EULERIAN_Z
      ! Phase-1 Lagrangian h-floor: only non-zero for VCOORD_LAGRANGIAN
      ! specifically (NOT the broad is_lagrangian = /= EULERIAN_Z — R3).
      h_min_floor = 0.0_wp
      if (present(vcoord)) then
         if (vcoord%coord_type == VCOORD_LAGRANGIAN) h_min_floor = ct%angstrom_h
      end if
      is_pc = dyn%split_scheme == SPLIT_SCHEME_PRED_CORR
      is_pred = is_pc .and. stage_id == 1
      ! Predictor: no tracer physics of any kind (MOM6's predictor never
      ! touches thermodynamics — SPEC §2; the corrector owns the full
      ! tracer chain at the thermo cadence).
      therm_active = dyn%enable_thermodynamics .and. dyn%is_thermo_step() &
                     .and. .not. is_pred
      therm_dt = dyn%therm_dt(dt)
      ! Velocity-apply dt: the predictor's provisional velocity advances
      ! to dt_pred = BE·dt (SPEC §2 P8); everything else (substep,
      ! continuity) stays at the full dt as MOM6 does.
      dt_vel = dt
      if (is_pred) dt_vel = dyn%pc_be*dt
      ! Mass-budget weight for the continuity chain: 0.5 per SSP stage;
      ! under pred_corr the predictor's h advance is DISCARDED (weight 0)
      ! and the corrector's is the whole step (weight 1).
      chain_weight = 0.5_wp
      if (is_pc) chain_weight = merge(0.0_wp, 1.0_wp, is_pred)

      call probe_dS(grid, ms, "entry", stage_id, step_id)

      ! ---- 0. Stage-entry periodic wrap ----
      ! For periodic axes: fill ghost cells of h_layer, u/v layer faces,
      ! and all tracers before derive_bt_from_layers, so every slow
      ! tendency kernel (EOS, Coriolis-adv, PGF, hvisc, bdrag, surface
      ! stress) sees physically-correct values at the seam.  No-op when
      ! neither periodic flag is set.
      ! O2: unconditional multi-rank ghost exchange (D0 — no-op on 1 rank).
      call ocean_halo_exchange_ml_state(ms)
      if (present(bc)) call ocean_periodic_wrap_state(grid, bc, ms, &
                                                      skip_x=ocean_halo_is_decomposed_x(), skip_y=ocean_halo_is_decomposed_y())
      ! Tripolar north-fold seam — periodic-FIRST-fold-SECOND (Appendix A):
      ! folds h/u/v/tracers AFTER the periodic wrap so it reads the already
      ! cyclically-wrapped corner columns.  No-op when bc%north_fold is off.
      if (present(bc)) call ocean_fold_wrap_state(grid, bc, ms)

      ! ---- 0b. pred_corr: the SAME seam fill for the step time-means ----
      ! Under `split_scheme = "pred_corr"` the Coriolis-advection and the
      ! horizontal-viscosity tendencies are evaluated on `u_av`/`v_av`/`h_av`,
      ! NOT on the prognostic fields that step 0 just wrapped.  Those means
      ! come out of the continuity solve (`u_cor`) and the h_in/h_out average,
      ! both of which write the INTERIOR only — so without this their ghost
      ! band holds stale values (zeros on the first step) and the seam
      ! tendency is wrong.
      !
      ! Symptom, and why this is not cosmetic: `periodic_shifted_domain_
      ! identity` (tests/test_ocean_periodic.F90) runs the same physical IC
      ! twice with run B's domain circularly shifted half a period and
      ! demands shift(B) == A BIT-FOR-BIT.  With the ghosts unfilled the
      ! answer depends on where the seam falls: max_diff = 8.6e-07 on h ~ 50
      ! after 8 steps, from exactly 0 under ssp_rk2.  ssp_rk2 never reads
      ! these arrays, so the `is_pc` gate keeps it bit-identical.
      if (is_pc) then
         if (allocated(ms%u_av_layer)) then
            ! Multi-rank seam first, exactly as step 0 does for the
            ! prognostics.  GATED on an actually-decomposed axis: the halo
            ! specifics take EXPLICIT-SHAPE dummies sized from the module's
            ! `oh_nx_total`/`oh_ny_total`, which a unit test that never calls
            ! `ocean_halo_init` leaves at 0 -- a mis-shaped device dummy, and
            ! on the GPU build that reads as garbage rather than as an error.
            if (ocean_halo_is_decomposed_x() .or. ocean_halo_is_decomposed_y()) then
               call ocean_halo_face_x(ms%u_av_layer, size(ms%u_av_layer, 3))
               call ocean_halo_face_y(ms%v_av_layer, size(ms%v_av_layer, 3))
               call ocean_halo_centre(ms%h_av_layer, size(ms%h_av_layer, 3))
            end if
            if (present(bc)) then
               if (bc%periodic_x .or. bc%periodic_y) then
                  call ocean_periodic_wrap_face_x_3d( &
                     ms%u_av_layer, size(ms%u_av_layer, 1), size(ms%u_av_layer, 2), &
                     size(ms%u_av_layer, 3), grid%nx_phys, grid%ny_phys, grid%nghost, &
                     bc%periodic_x .and. .not. ocean_halo_is_decomposed_x(), &
                     bc%periodic_y .and. .not. ocean_halo_is_decomposed_y())
                  call ocean_periodic_wrap_face_y_3d( &
                     ms%v_av_layer, size(ms%v_av_layer, 1), size(ms%v_av_layer, 2), &
                     size(ms%v_av_layer, 3), grid%nx_phys, grid%ny_phys, grid%nghost, &
                     bc%periodic_x .and. .not. ocean_halo_is_decomposed_x(), &
                     bc%periodic_y .and. .not. ocean_halo_is_decomposed_y())
                  call ocean_periodic_wrap_centre_3d( &
                     ms%h_av_layer, size(ms%h_av_layer, 1), size(ms%h_av_layer, 2), &
                     size(ms%h_av_layer, 3), grid%nx_phys, grid%ny_phys, grid%nghost, &
                     bc%periodic_x .and. .not. ocean_halo_is_decomposed_x(), &
                     bc%periodic_y .and. .not. ocean_halo_is_decomposed_y())
               end if
               ! Tripolar north fold of the time-means, periodic-FIRST (the
               ! wrap above): the Coriolis/hvisc stencils that read
               ! u_av/v_av/h_av at the seam need the mirrored north ghosts and
               ! the antisymmetric fold-line v_av exactly as step 0 gives the
               ! prognostics.  Without it the pred_corr seam tendencies read
               ! the interior-only continuity output (stale ghost rows).
               ! px > 1: one owner-routed exchange group (u_av, v_av, h_av).
               if (bc%north_fold) call ocean_fold_wrap_time_means(grid, bc, ms)
            end if
         end if
      end if

      ! KE attribution: stage-entry baseline (debug_ke_attr; no-op when off).
      call ke_probe_sample(grid, ms, dyn%ke_probe, "entry", stage_id, step_id)
      call chksum_state(grid, ms, dyn%chksum_probe, "entry", stage_id, step_id)

      ! ---- 1. Snapshot u_bt^n, v_bt^n at start of stage ----
      call derive_bt_from_layers(grid, dyn%bt_work, ms, metrics)
      ! Build the per-face upstream-PPM column-sum thickness on the
      ! same snapshot.  No-op when `use_upstream_h_face = .false.`;
      ! otherwise feeds the BT substep + corrector with the same
      ! face-thickness convention slow continuity uses, eliminating
      ! the centred-vs-upstream mismatch at slopes.
      call compute_h_face_upstream(grid, dyn%bt_work, ms, metrics)
      ! Build the per-face BT_cont_type flux closure from the same ML
      ! snapshot.  No-op when `use_bt_cont_type = .false.`; otherwise
      ! BTCL_u/v feed the BT substep's flux paths instead of the
      ! naive `uh = u·h_face`.
      call set_local_BT_cont_types(grid, metrics, dyn%bt_work, ms, dt)
      nx_face = size(dyn%bt_work%bt_ubt, 1)
      ny_uface = size(dyn%bt_work%bt_ubt, 2)
      nx_vface = size(dyn%bt_work%bt_vbt, 1)
      ny_face = size(dyn%bt_work%bt_vbt, 2)
      do concurrent(j=1:ny_uface, i=1:nx_face)
         dyn%bt_work%ubt_at_n(i, j) = dyn%bt_work%bt_ubt(i, j)
      end do
      do concurrent(j=1:ny_face, i=1:nx_vface)
         dyn%bt_work%vbt_at_n(i, j) = dyn%bt_work%bt_vbt(i, j)
      end do

      ! Snapshot η before the slow PGF runs — used by the bc-PGF
      ! correction below to form `e_anom = (η_avg − eta_PF)`.
      ! No-op when bt_correction_bc_pgf is off; safe to always run.
      if (dyn%bt_work%bt_correction_bc_pgf) then
         call snapshot_eta_PF(dyn%bt_work)
      end if

      ! ---- 2. Compute all slow tendencies ----
      ! All read h^n / u^n / v^n into per-kernel scratch buffers.
      ! No writes to h_layer, hTr, u_face, or v_face yet.
      call profiler_start("ocean_eos")
      ! See the note at the other ocean_eos_compute call site: the EOS is a
      ! dynamics term and must NOT be held across the thermo window. Note the
      ! predictor exclusion (`.not. is_pred`) also drops away here — the
      ! predictor's PGF needs a current rho just as much as the corrector's.
      call ocean_eos_compute(eos, ms, active=dyn%enable_thermodynamics)
      call profiler_stop("ocean_eos")
      call probe_dS(grid, ms, "after eos", stage_id, step_id)
      call profiler_start("ocean_coriolis_adv")
      ! SPEC S3 (`split_scheme = "pred_corr"`): the Coriolis-advection
      ! tendency reads the `u_av`/`h_av` step time-means, never the
      ! prognostic (MOM6 `CorAdCalc(u_av, v_av, h_av, ...)`).
      ! Default ⇒ prognostic, bit-identical.
      !
      ! Mass-consistent CorAdCalc (`&ocean_coriolis_nml use_state_fluxes`):
      ! in the pred_corr CORRECTOR, `ms%mass_flux_*_layer` still hold the
      ! PREDICTOR chain's renormalised transports — the uh/vh from the same
      ! solve that produced this `u_av` evaluation state.  Consuming them
      ! (instead of the kernel's own `u·h_face` recompute) closes the
      ! evaluate-at-u_av / transport-mismatch energy leak the KE-attribution
      ! meter pinned on coriolis_adv (pdc M16, the rim-Kelvin slow
      ! exponential).  Predictor stages keep the recompute.
      if (dyn%split_scheme == SPLIT_SCHEME_PRED_CORR) then
         call coriolis_adv_compute_tendencies(grid, metrics, cor, ms, &
                                              u_src=ms%u_av_layer, &
                                              v_src=ms%v_av_layer, &
                                              h_src=ms%h_av_layer, &
                                              use_state_fluxes=(cor%state_fluxes &
                                                                .and. stage_id == 2))
      else
         call coriolis_adv_compute_tendencies(grid, metrics, cor, ms)
      end if
      call profiler_stop("ocean_coriolis_adv")
      call profiler_start("ocean_pgf")
      call ocean_pressure_force_compute(grid, metrics, pgf, ms, eos=eos)
      call profiler_stop("ocean_pgf")
      ! bc-PGF correction: compute per-layer `pbce` + face-centred
      ! `gtot_*` now (PGF has populated `pgf%e_face`).  Used later
      ! by `apply_bt_correction` with `use_bc_pgf=.true.`.  Gated;
      ! no-op when bt_correction_bc_pgf is off.
      if (dyn%bt_work%bt_correction_bc_pgf) then
         call compute_pbce(grid, dyn%bt_work, pgf, ms)
         call compute_gtot_faces(grid, dyn%bt_work, ms, metrics)
      end if
      call profiler_start("ocean_hvisc")
      ! pred_corr predictor: SKIP the viscous recompute — MOM6's predictor
      ! uses the PREVIOUS step's diffu (SPEC §2 P3, "diffu(u[n-1])"), and
      ! the du_visc/dv_visc buffers persist from the last corrector's C1.
      if (.not. is_pred) then
         ! Flow-aware lateral closure first (Leith / Smagorinsky) so the
         ! hvisc kernel reads the freshly-computed per-face viscosity.
         ! Both calls self-gate on `lateral_mix` absence / `LMIX_NONE`.
         ! Gap 1: when `resoln_scaled_visc` is on AND VarMix is active, forward
         ! the resolution-function face fields so the dynamic coefficients are
         ! scaled by `Res_fn` before the clamps.  Absent / off ⇒ the plain call
         ! ⇒ bit-identical.
         if (lateral_mix_uses_resoln(lateral_mix, varmix)) then
            call ocean_lateral_mix_compute(grid, metrics, lateral_mix, ms, &
                                           res_fn_u=varmix%res_fn_u, &
                                           res_fn_v=varmix%res_fn_v)
         else
            call ocean_lateral_mix_compute(grid, metrics, lateral_mix, ms)
         end if
         ! MEKE backscatter (capability [5], Gap 2): subtract the eddy-energy
         ! harmonic backscatter `ku` (CFL-floored) from the freshly-computed
         ! per-face resolved viscosity, so the hvisc Laplacian returns energy to
         ! the resolved flow.  Reads the prior thermo step's `ku` (one-step lag).
         ! No-op unless `meke%backscatter`; needs the face fields (LMIX active).
         if (present(meke) .and. present(lateral_mix)) then
            if (lateral_mix%is_init) then
               call meke_backscatter_apply(grid, metrics, meke, dt, &
                                           lateral_mix%ah_face_x, lateral_mix%ah_face_y)
            end if
         end if
         ! SPEC S3: under pred_corr the viscous tendency also reads the
         ! time-means (MOM6 `horizontal_viscosity(u_av, v_av, h_av, ...)`).
         ! The lateral-mix COEFFICIENT
         ! computation above still reads the prognostic (documented
         ! deviation: the closure coefficient is one flow-state stale,
         ! second-order in dt; full parity is a follow-on).
         if (dyn%split_scheme == SPLIT_SCHEME_PRED_CORR) then
            call ocean_horizontal_viscosity_compute_tendencies(grid, metrics, hv, ms, &
                                                               lateral_mix=lateral_mix, dt=dt, &
                                                               u_src=ms%u_av_layer, &
                                                               v_src=ms%v_av_layer, &
                                                               h_src=ms%h_av_layer)
         else
            call ocean_horizontal_viscosity_compute_tendencies(grid, metrics, hv, ms, &
                                                               lateral_mix=lateral_mix, dt=dt)
         end if
         ! KE dissipation rate for the MEKE frictional source (no-op unless on).
         call ocean_horizontal_viscosity_compute_ke_diss(hv, ms)
      end if
      call profiler_stop("ocean_hvisc")
      call profiler_start("ocean_bdrag")
      call ocean_bottom_drag_compute_tendencies(grid, bd, ms, dt)
      call ocean_channel_drag_compute_tendencies(grid, metrics, bd, ms)
      ! Ice-shelf top drag (`&ocean_tdrag_nml`) — see `run_stage`.
      if (present(td)) call ocean_top_drag_compute_tendencies(td, ms, dt)
      ! Phase 4b: publish the ice-shelf base stress that BOTH boundary-
      ! layer schemes take `u_*` from.  Under a shelf the wind has been
      ! masked out of `tau` (so `stress_mag` is exactly 0 there) and the
      ! turbulent boundary layer is driven by the ice-ocean stress
      ! instead — `u_*^2 = |tau_top|/rho_0`.
      !
      ! Written INLINE as a `do concurrent`, not as a call handing
      ! `ss%stress_shelf` to an external subroutine: a host-gated call
      ! with a state array as an actual makes nvfortran treat the array
      ! as escaping and pessimises every `do concurrent` in this routine
      ! (CLAUDE.md, measured at +4.8%% on an inert porous pass).
      !
      ! Gated on `td%enable`, not `present(td)`: a DISABLED slot carries
      ! a `(1,1)` placeholder `stress_top`.
      !
      ! Placed here — after the top-drag compute, before
      ! `vmix_apply_in_stage` below — so KPP/EPBL read THIS stage's
      ! stress.  No lag.
      publish_shelf = .false.
      if (present(td)) publish_shelf = td%enable
      if (publish_shelf) then
         nx_ss = size(ss%stress_shelf, 1)
         ny_ss = size(ss%stress_shelf, 2)
         do concurrent(j_ss=1:ny_ss, i_ss=1:nx_ss)
            ss%stress_shelf(i_ss, j_ss) = td%stress_top(i_ss, j_ss)
         end do
      end if
      call profiler_stop("ocean_bdrag")
      call profiler_start("ocean_surfstress")
      call ocean_surface_stress_compute_tendencies(grid, ss, ms)
      call profiler_stop("ocean_surfstress")
      call probe_dS(grid, ms, "after slow computes", stage_id, step_id)

      ! Per-edge OBC tags + physical-edge flags — consumed by
      ! subtract_fast_cor_ref (below) and the post-substep transport
      ! wall-zeroing.  Defaults: WALL tags, .true. flags ⇒ single-rank
      ! closed-basin behaviour.
      bc_w_drv = OBC_WALL
      bc_e_drv = OBC_WALL
      bc_s_drv = OBC_WALL
      bc_n_drv = OBC_WALL
      has_w_drv = .true.
      has_e_drv = .true.
      has_s_drv = .true.
      has_n_drv = .true.
      if (present(bc)) then
         bc_w_drv = ocean_bc_outer_face_tag(bc%west%bc_type)
         bc_e_drv = ocean_bc_outer_face_tag(bc%east%bc_type)
         bc_s_drv = ocean_bc_outer_face_tag(bc%south%bc_type)
         bc_n_drv = ocean_bc_outer_face_tag(bc%north%bc_type)
         has_w_drv = bc%has_west
         has_e_drv = bc%has_east
         has_s_drv = bc%has_south
         has_n_drv = bc%has_north
      end if

      ! ---- 3. Sum slow tendencies into F_slow_u/v + depth-mean to F_bt ----
      ! Builds the constant-per-substep forcing for the barotropic substep.
      ! Uses h^n (which is still untouched) for the depth weighting.
      call profiler_start("ocean_F_slow_assembly")
      ! BT-budget probe: per-region per-term BT power decomposition.
      ! Pull device buffers down to host first (the probe is host-side
      ! Fortran loops + writes to stdout).  Gated on the namelist knob
      ! so production runs pay nothing.
      if (dyn%debug_bt_budget) then
         !$acc update self(ms%h_layer)
         !$acc update self(ms%u_face_x_layer, ms%v_face_y_layer)
         !$acc update self(dyn%bt_work%bt_eta, dyn%bt_work%bt_ubt, dyn%bt_work%bt_vbt)
         !$acc update self(pgf%dpdx_face%data, pgf%dpdy_face%data)
         !$acc update self(cor%pv_flux_x%data, cor%pv_flux_y%data)
         !$acc update self(hv%du_visc%data, hv%dv_visc%data)
         !$acc update self(bd%du_drag%data, bd%dv_drag%data)
         !$acc update self(ss%du_stress%data, ss%dv_stress%data)
         call print_bt_budget(grid, ms, dyn%bt_work, pgf, cor, hv, bd, ss, &
                              real(step_id, wp), "step", "S"//achar(48 + stage_id), &
                              header=(step_id == 1 .and. stage_id == 1))
      end if
      ! Pre-substep visc_rem refresh (MOM6 order parity): the rem-weighted
      ! forcing below, the γ-weighted continuity renormaliser, and the Δu
      ! corrector all read visc_rem_u/v; without this refresh they lag one
      ! stage — fatally at step 1 stage 1, where the init value (≡ 1)
      ! makes F_bt the plain mean of the ballistic spurious tendencies
      ! (PGF_BUG.md §9).  MOM6 computes vertvisc_coef + remnant in the
      ! predictor BEFORE btstep/continuity (SPEC §2 P5).
      if (dyn%bt_work%bt_forcing_visc_rem .or. dyn%bt_work%bt_renorm_visc_rem &
          .or. dyn%bt_work%bt_rem_from_visc_rem .or. is_pc) then
         if (fold_top) then
            call visc_rem_precompute(grid, dyn%bt_work, vmix, vd, ss, bd, ms, dt, kshear=kshear, &
                                     lambda_top_u=td%lambda_top_u, lambda_top_v=td%lambda_top_v, &
                                     cover_u=td%cover_u, cover_v=td%cover_v, bc=bc)
         else
            call visc_rem_precompute(grid, dyn%bt_work, vmix, vd, ss, bd, ms, dt, kshear=kshear, bc=bc)
         end if
      end if
      call sum_slow_tendencies_into_F_slow(dyn%bt_work, pgf, cor, hv, bd, ss, ms)
      ! The top drag MUST reach the barotropic mode the same way the
      ! bottom drag does — through the depth mean of `F_slow`, which the
      ! substep integrates and `apply_bt_correction` then subtracts back
      ! out.  See `add_top_drag_into_F_slow`'s docstring for why a
      ! tendency left out of this sum is both invisible to the fast loop
      ! and mis-corrected on the layers.
      if (present(td)) call add_top_drag_into_F_slow(dyn%bt_work, td, ms)
      ! MOM6 wt_u parity (`&ocean_bt_nml forcing_visc_rem`): weight the
      ! forcing depth-mean by h·visc_rem so layers the implicit friction
      ! will immediately damp (grounded stacks under the vdiff BBL glue)
      ! do not force the fast loop (PGF_BUG.md §9).  The PGF-projection
      ! subtraction below MUST use the
      ! same weights — a plain-h-weighted subtraction would re-inject the
      ! very depth-mean the weighting removed, with the opposite sign.
      ! visc_rem lags one stage (the stage-end vdiff producer), same as
      ! the corrector's documented convention; it initializes to 1 so the
      ! first stage degenerates to the plain h-mean.
      if (dyn%bt_work%bt_forcing_visc_rem) then
         call face_depth_mean_rem_u(grid, dyn%bt_work%F_slow_u, ms%h_layer, &
                                    dyn%bt_work%visc_rem_u, dyn%bt_work%F_bt_u, ms%nz_ml, metrics, &
                                    n_inner, dyn%bt_work%bt_H_ref, dyn%bt_work%frhat_scheme)
         call face_depth_mean_rem_v(grid, dyn%bt_work%F_slow_v, ms%h_layer, &
                                    dyn%bt_work%visc_rem_v, dyn%bt_work%F_bt_v, ms%nz_ml, metrics, &
                                    n_inner, dyn%bt_work%bt_H_ref, dyn%bt_work%frhat_scheme)
      else
         call face_depth_mean_u(grid, dyn%bt_work%F_slow_u, ms%h_layer, dyn%bt_work%F_bt_u, ms%nz_ml, metrics, &
                                dyn%bt_work%bt_H_ref, dyn%bt_work%frhat_scheme)
         call face_depth_mean_v(grid, dyn%bt_work%F_slow_v, ms%h_layer, dyn%bt_work%F_bt_v, ms%nz_ml, metrics, &
                                dyn%bt_work%bt_H_ref, dyn%bt_work%frhat_scheme)
      end if
      ! Remove from the forcing the part of the slow PGF the barotropic
      ! substep re-represents with its own live `-G·∂η/∂x`, or the bt
      ! mode integrates it twice and √(gH) inflates to √(2gH).
      !
      ! `bc_pgf_forcing` (default, MOM6 `BT_force` + `eta_PF`): that part
      ! is the free-surface term the slow PGF carries, `-g_pf·∇η_PF` at the
      ! η it was built on (`g_pf = 0` for the surface-relative MONT/FV_LITE/
      ! FV_WRIGHT forms) — so the depth-mean BAROCLINIC PGF stays in the
      ! forcing.  Legacy (`.false.`): the WHOLE depth-mean PGF is
      ! subtracted, which also throws away its baroclinic part (the JEBAR /
      ! bottom-pressure forcing of the barotropic mode).
      if (dyn%bt_work%bt_bc_pgf_forcing) then
         if (present(eta_pf_seam)) then
            call set_fast_forcing_eta_pf(grid, metrics, dyn%bt_work, grid%nx_total, grid%ny_total, &
                                         pgf_free_surface_gravity(pgf), eta_pf_seam, .true.)
         else
            call set_fast_forcing_eta_pf(grid, metrics, dyn%bt_work, grid%nx_total, grid%ny_total, &
                                         pgf_free_surface_gravity(pgf), dyn%bt_work%bt_eta, .false.)
         end if
      else if (dyn%bt_work%bt_forcing_visc_rem) then
         call face_depth_mean_rem_u(grid, pgf%dpdx_face%data, ms%h_layer, &
                                    dyn%bt_work%visc_rem_u, dyn%bt_work%F_bt_u_fast, ms%nz_ml, metrics, &
                                    n_inner, dyn%bt_work%bt_H_ref, dyn%bt_work%frhat_scheme)
         call face_depth_mean_rem_v(grid, pgf%dpdy_face%data, ms%h_layer, &
                                    dyn%bt_work%visc_rem_v, dyn%bt_work%F_bt_v_fast, ms%nz_ml, metrics, &
                                    n_inner, dyn%bt_work%bt_H_ref, dyn%bt_work%frhat_scheme)
      else
         call face_depth_mean_u(grid, pgf%dpdx_face%data, ms%h_layer, dyn%bt_work%F_bt_u_fast, ms%nz_ml, metrics, &
                                dyn%bt_work%bt_H_ref, dyn%bt_work%frhat_scheme)
         call face_depth_mean_v(grid, pgf%dpdy_face%data, ms%h_layer, dyn%bt_work%F_bt_v_fast, ms%nz_ml, metrics, &
                                dyn%bt_work%bt_H_ref, dyn%bt_work%frhat_scheme)
      end if
      if (.not. dyn%bt_work%bt_bc_pgf_forcing) then
         do concurrent(j=1:ny_uface, i=1:nx_face)
            dyn%bt_work%F_bt_u_fast(i, j) = dyn%bt_work%F_bt_u(i, j) - dyn%bt_work%F_bt_u_fast(i, j)
         end do
         do concurrent(j=1:ny_face, i=1:nx_vface)
            dyn%bt_work%F_bt_v_fast(i, j) = dyn%bt_work%F_bt_v(i, j) - dyn%bt_work%F_bt_v_fast(i, j)
         end do
      end if
      ! Subtract the fast-loop Coriolis + advection evaluated at the
      ! stage-entry bt state — the Coriolis/advection analogue of the PGF
      ! projection subtraction just above (MOM6 `Cor_ref_u/v`).  Without
      ! it the barotropic Coriolis is integrated twice (frozen in F_bt_u
      ! via `cor%pv_flux_*`'s depth mean AND live in the substep), which
      ! pumps an exponential wall-trapped barotropic mode on shelf rims
      ! (the 600² Lagrangian double-gyre h-guard trap).  `bt_ubt/bt_vbt`
      ! still hold the stage-entry state here (`derive_bt_from_layers`
      ! filled them; the substep's copy-in below re-sets them).  The
      ! per-edge tags/flags are derived early (above) for this call.
      !
      ! Wet/dry exclusion (v1 composition deferral, same class as the
      ! wetdry × use_cont_type/upstream_h_face exclusions): the reference
      ! is evaluated on the ungated stage-entry velocities, but the
      ! substep's live terms are zeroed at dry faces by `wd_open_u/v` —
      ! subtracting an ungated reference there injects net forcing at
      ! dry faces and drives `h_layer` negative
      ! (test_ocean_wetdry_driver).  Skipping keeps wetdry runs on the
      ! pre-fix forcing (the double-count is mild at coastal scales).
      if (.not. dyn%bt_work%wetdry_enable) then
         ! Reference velocity first.  Under `ssp_rk2` this is a copy of
         ! the stage-entry `bt_ubt/bt_vbt` the reference used to read
         ! directly (bit-identical); under `pred_corr` it is the depth
         ! mean of `u_av/v_av` — the velocity step 2 actually evaluated
         ! `cor%pv_flux_*` on — under the same weights the forcing
         ! depth-mean above used.  Without this the uncancelled
         ! `f × (v̄_av − v̄^n)` forces every substep and pumps the basin's
         ! gravest Poincaré seiche (see `set_cor_ref_velocity`).
         call set_cor_ref_velocity(grid, dyn%bt_work, ms, is_pc, metrics, n_inner)
         call subtract_fast_cor_ref(grid, metrics, dyn%bt_work, cor%f_corner, &
                                    bc_w_drv, bc_e_drv, bc_s_drv, bc_n_drv, &
                                    has_w_drv, has_e_drv, has_s_drv, has_n_drv)
      end if
      call profiler_stop("ocean_F_slow_assembly")

      ! ---- 4. Run barotropic substep ----
      ! Outputs bt_eta_end, bt_ubt_end, bt_vbt_end (Hallberg end-step
      ! anchors for the Δu correction), bt_eta / bt_ubt / bt_vbt
      ! (time-mean for diagnostics), and — new — bt_uhbt / bt_vhbt
      ! (time-mean depth-integrated transport).  The slow continuity
      ! below renormalises its per-layer mass fluxes to vertically
      ! sum to those, so the layer h evolution is consistent with
      ! the barotropic-substep η evolution by construction (MOM6 split-RK2).
      do concurrent(j=1:ny_uface, i=1:nx_face)
         dyn%bt_work%bt_ubt(i, j) = dyn%bt_work%ubt_at_n(i, j)
      end do
      do concurrent(j=1:ny_face, i=1:nx_vface)
         dyn%bt_work%bt_vbt(i, j) = dyn%bt_work%vbt_at_n(i, j)
      end do
      ! bt_eta is fresh from derive_bt_from_layers — leave it.

      dt_inner = dt/real(n_inner, wp)
      ! MOM6 bt_rem_u: populate the per-face multiplicative damping
      ! factor read by the substep loop.  When the knob is off, leave
      ! `bt_rem_u/v` at their init value of 1 ⇒ multiplication is a
      ! no-op (bit-identical to pre-knob path).
      !
      ! PR-2 (bt-rem-from-av-rem): `bt_rem_from_visc_rem` is a THIRD
      ! resetter, mutually exclusive with `bt_substep_drag` at configure
      ! (D2 — double-counted bed drag) — so this if/else-if chain still
      ! dispatches to exactly one resetter per stage, never more than
      ! one, preserving the multiplicative-accumulator contract (`src/
      ! core/ocean/README.md`).  Built from the SAME visc_rem producer
      ! the BT corrector reads (MOM6's barotropic solver), once
      ! per barotropic call, BEFORE the substeps below.
      if (dyn%bt_work%bt_rem_from_visc_rem) then
         call compute_bt_rem_from_visc_rem(grid, dyn%bt_work, ms, metrics, n_inner)
      else if (dyn%bt_work%bt_substep_drag) then
         call compute_bt_rem(grid, dyn%bt_work, ms, metrics, bd%r_linear, bd%hbbl, dt_inner)
      else if (dyn%bt_work%lwd_enable) then
         ! `bt_rem_u/v` is reset ONLY by `compute_bt_rem`/
         ! `compute_bt_rem_from_visc_rem` above; when both are off but
         ! wave drag is on, nothing else resets it, and
         ! `compute_bt_rem_wave_drag` below MULTIPLIES into it —
         ! without this reset bt_rem would compound geometrically across
         ! outer steps (bt_rem = R^n after n stages), silently annihilating
         ! the barotropic mode. See `src/core/ocean/README.md`.
         call reset_bt_rem(grid, dyn%bt_work)
      end if
      if (dyn%bt_work%lwd_enable) then
         ! Barotropic linear wave drag (Egbert & Ray 2001; Jayne & St
         ! Laurent 2001): MULTIPLIES the static piston-velocity map into
         ! `bt_rem_u/v`, composing with `substep_drag` exactly as MOM6
         ! composes `lin_drag_u` with the viscous remnant.
         ! Must run AFTER the base fill
         ! above and BEFORE `mask_bt_rem` (land masking must be last).
         call compute_bt_rem_wave_drag(grid, dyn%bt_work, ms, metrics, dt_inner)
      end if
      ! Fold the static land face masks into bt_rem (C4 / R4a): zeroes the
      ! BT-substep velocity update across land faces.  No-op for all-wet.
      call mask_bt_rem(grid, metrics, dyn%bt_work)
      call profiler_start("ocean_barotropic_solver")
      ! `eta_forcing` is itself optional here: passing an absent optional
      ! as the actual for the substep's optional dummy propagates absence
      ! (F2018 15.5.2.13) ⇒ bit-identical when the C1 tide is off.
      if (dyn%bt_halo > 0) then
         ! ---- Wide-halo march-in path (Phase 3c) ----
         ! Tides × bt_halo is a configure-time exclusion.
         if (present(eta_forcing)) then
            error stop "run_stage_split: bt_halo > 0 with eta_forcing — excluded at configure time"
         end if
         call dyn%bt_wide%copy_in(grid, &
                                  dyn%bt_work%bt_eta, dyn%bt_work%bt_H_ref, &
                                  dyn%bt_work%bt_ubt, dyn%bt_work%bt_vbt, &
                                  dyn%bt_work%bt_ubt_prev, dyn%bt_work%bt_vbt_prev, &
                                  dyn%bt_work%bt_rem_u, dyn%bt_work%bt_rem_v, &
                                  dyn%bt_work%F_bt_u_fast, dyn%bt_work%F_bt_v_fast)
         call dyn%bt_wide%entry_exchange()
         call bt_wide_substep(dyn%bt_wide, dyn%bt_work, n_inner, dt_inner, bc)
         ! Copy wide outputs back to normal-width bt_work fields.
         ! The copy_out DC loops are synchronous (no acc async); drain async(1)
         ! first so copy_out reads the completed time-mean arrays.
         !$acc wait(1)
         call dyn%bt_wide%copy_out(grid, &
                                   dyn%bt_work%bt_eta, dyn%bt_work%bt_ubt, dyn%bt_work%bt_vbt, &
                                   dyn%bt_work%bt_uhbt, dyn%bt_work%bt_vhbt, &
                                   dyn%bt_work%bt_eta_end, dyn%bt_work%bt_ubt_end, &
                                   dyn%bt_work%bt_vbt_end)
         ! Exit seam freshen: leaves downstream consumers exactly what v1's
         ! last in-loop exchange left (fresh normal-width seam ghosts).
         call profiler_start("ocean_comms_bt")
         call ocean_halo_bt_group_2d(dyn%bt_work%bt_eta, dyn%bt_work%bt_ubt, &
                                     dyn%bt_work%bt_vbt)
         call profiler_stop("ocean_comms_bt")
      else
         call barotropic_substep_nonlinear_interior(grid, metrics, dyn%bt_work, &
                                                    cor%f_corner, n_inner, dt_inner, &
                                                    bc=bc, t=t, eta_forcing=eta_forcing)
      end if
      call profiler_stop("ocean_barotropic_solver")
      call probe_dS(grid, ms, "after barotropic substep", stage_id, step_id)
      ! BT in/out chksums at the substep exit: when the fold's loud count
      ! fires, these rows name which fold INPUT (ubt_end/ubt_at_n/F_bt)
      ! went non-finite — the producer the nan-catch cannot see.
      call chksum_bt(grid, dyn%bt_work, dyn%chksum_probe, "post_bt", stage_id, step_id)

      ! Physical-wall reconciliation for the transport constraint.
      ! The barotropic substep closes `bt_ubt` at array-edge faces (i=1,
      ! i=nx+1) — convenient for its own integration but several
      ! cells away from where the slow continuity closes (the
      ! physical walls at i = nghost+1 and i = nghost+nx_phys+1).
      ! Passing the raw barotropic-substep transport into continuity would
      ! ask for nonzero mass flux through the physical wall, which
      ! the wall-zero step then erases — leaving the constraint
      ! unsatisfied at those faces and a small leak in
      ! `sum_k(h_layer)` vs `H_ref + bt_eta_end`.  Zero the wall
      ! faces of `bt_uhbt / bt_vhbt` so the constraint is
      ! self-consistent with the closed-wall slow continuity.
      !
      ! OBC dispatch: for non-WALL edges, the barotropic substep already
      ! computed a physical transport at the wall face (Flather,
      ! Chapman, clamped data) — DON'T zero it.  The slow continuity
      ! will honour the same `bc` tag and let the matching mass
      ! flux through.  (Tags/flags bc_*_drv / has_*_drv are derived
      ! earlier, before the F_slow assembly, for subtract_fast_cor_ref.)
      ! Seam faces carry real transport; the halo owns them (O0, plan D4).
      do concurrent(j=1:ny_uface)
         if (bc_w_drv == OBC_WALL .and. has_w_drv) dyn%bt_work%bt_uhbt(grid%nghost + 1, j) = 0.0_wp
         if (bc_e_drv == OBC_WALL .and. has_e_drv) dyn%bt_work%bt_uhbt(grid%nghost + grid%nx_phys + 1, j) = 0.0_wp
      end do
      do concurrent(i=1:nx_vface)
         if (bc_s_drv == OBC_WALL .and. has_s_drv) dyn%bt_work%bt_vhbt(i, grid%nghost + 1) = 0.0_wp
         if (bc_n_drv == OBC_WALL .and. has_n_drv) dyn%bt_work%bt_vhbt(i, grid%nghost + grid%ny_phys + 1) = 0.0_wp
      end do

      ! Continuity + tracer chain.  ssp_rk2: here (historical position,
      ! bit-identical).  pred_corr: deferred to AFTER the velocity update +
      ! implicit friction — the forward-backward pairing (SPEC §2 C8).
      if (.not. is_pc) then
         call run_continuity_chain(grid, metrics, dyn, ct, hd, va, redi, varmix, ms, &
                                   dt, therm_dt, therm_active, is_lagrangian, &
                                   h_min_floor, chain_weight, is_pred, &
                                   stage_id, step_id, bc=bc, mle=mle)
      end if

      ! ---- 6. Velocity-tendency applies ----
      ! Async-chained on OpenACC queue 1: these five applies are all
      ! additive forward-Euler accumulations onto u_face/v_face with no
      ! intervening default-queue or host-reading op between them (the
      ! profiler calls are host-side wall-clock timers + NVTX markers; they
      ! neither sync the device nor read device data).  Same queue ⇒ FIFO
      ! ⇒ the additive sequence is order-preserved.  ONE `!$acc wait(1)`
      ! before `apply_bt_correction` — the first device consumer that reads
      ! the freshly-applied u_face/v_face — completes the chain.  This
      ! eliminates the per-launch host sync gap (~half the per-stage apply
      ! cost in nsys).
      ! accel_visc_rem (MOM6 parity): snapshot the pre-apply velocity so
      ! the post-apply reweight below can attenuate the whole explicit
      ! tendency sum by the per-layer viscous remnant.  The `!$acc wait`
      ! orders the snapshot against any still-in-flight async producer of
      ! u_face/v_face (the q1 chain convention).
      if (dyn%accel_visc_rem) then
         !$acc wait
         call accel_visc_rem_snapshot(grid%nx_total + 1, grid%ny_total, ms%nz_ml, &
                                      ms%u_face_x_layer, dyn%avr_u0)
         call accel_visc_rem_snapshot(grid%nx_total, grid%ny_total + 1, ms%nz_ml, &
                                      ms%v_face_y_layer, dyn%avr_v0)
      end if
      call profiler_start("ocean_velocity_apply")
      ! `dt_vel` = dt (historical / corrector) or BE·dt (pred_corr
      ! predictor, SPEC §2 P8).
      ! chksum seams (&ocean_debug_nml chksum, inert when off): each sample
      ! waits the device, so the async queue-1 chain is serialised ONLY inside
      ! the chksum window — the per-tendency attribution that named the
      ! 30 m/s injector.  FIFO order is unaffected (same queue).
      call coriolis_adv_apply_tendencies(cor, ms, dt_vel, no_wait=.true.)
      call chksum_state(grid, ms, dyn%chksum_probe, "post_coradv", stage_id, step_id)
      call chksum_hotface(grid, ms, dyn%bt_work%visc_rem_u, dyn%bt_work%visc_rem_v, &
                          dyn%chksum_probe, "post_coradv", stage_id, step_id)
      call ocean_pressure_force_apply(pgf, ms, dt_vel, no_wait=.true.)
      call chksum_state(grid, ms, dyn%chksum_probe, "post_pgf", stage_id, step_id)
      call chksum_hotface(grid, ms, dyn%bt_work%visc_rem_u, dyn%bt_work%visc_rem_v, &
                          dyn%chksum_probe, "post_pgf", stage_id, step_id)
      call ocean_horizontal_viscosity_apply_tendencies(hv, ms, dt_vel, no_wait=.true.)
      call chksum_state(grid, ms, dyn%chksum_probe, "post_hvisc", stage_id, step_id)
      ! Double-count guard (see run_stage): skip the explicit per-layer
      ! drag/stress apply when it is folded into the implicit vdiff
      ! tridiagonal.  NOTE (split path): the drag/stress tendency buffers
      ! still feed F_slow → the barotropic mode above; the implicit fold
      ! here acts on the post-bt-correction per-layer (baroclinic)
      ! velocity.  PR-19 supplies the visc_rem QUANTITY (produced by
      ! vdiff_apply_momentum, consumed by apply_bt_correction's
      ! visc_rem-weighted fold below); feeding it into F_slow / the
      ! barotropic substep itself — the actual barotropic-coupling flip
      ! — is PR-56's territory (changes the barotropic mode, gated on
      ! the Bleck/Hallberg instability test).  For strict split-path use
      ! today, keep the explicit / &ocean_bdrag_nml implicit split-apply.
      ! Defaults (folds off) ⇒ both applies run ⇒ bit-identical.  The MOM6
      ! BBL glue's piston is the bed drag too (see `run_stage`).
      if (.not. (vd%implicit_drag .or. vd%bbl_glue)) then
         call ocean_bottom_drag_apply_tendencies(bd, ms, dt_vel, no_wait=.true.)
      end if
      call ocean_channel_drag_apply_tendencies(bd, ms, dt_vel, no_wait=.true.)
      ! Top drag: same double-count guard as the bed (see `run_stage`).
      if (present(td)) then
         if (.not. td%implicit_fold) then
            call ocean_top_drag_apply_tendencies(td, ms, dt_vel, no_wait=.true.)
         end if
      end if
      if (.not. vd%implicit_stress) then
         call ocean_surface_stress_apply_tendencies(ss, ms, dt_vel, no_wait=.true.)
      end if
      call chksum_state(grid, ms, dyn%chksum_probe, "post_drag", stage_id, step_id)
      call profiler_stop("ocean_velocity_apply")
      call probe_dS(grid, ms, "after velocity applies", stage_id, step_id)

      ! ---- 7. bt correction ----
      ! In Eulerian-z mode: Δu correction + h-rescale (the rescale
      ! brings any small FP drift in `sum_k(h_layer)` back to
      ! `H_ref + bt_eta_end`, important for gravity-wave temporal
      ! continuity).  In Lagrangian mode: only the Δu correction —
      ! the slow continuity's `h_layer` is authoritative, and any
      ! FP residual gets absorbed by the ALE remap at end of outer
      ! step rather than redistributed across layers (which would
      ! corrupt S = hTr/h).
      ! bc-PGF correction: form `e_anom` from the post-substep
      ! η state and the pre-PGF snapshot.  Then apply_bt_correction
      ! consumes it via the optional `use_bc_pgf` path.  Gated;
      ! no-op when bt_correction_bc_pgf is off.
      if (dyn%bt_work%bt_correction_bc_pgf) then
         call compute_e_anom(dyn%bt_work)
      end if
      ! Close the async velocity-apply chain (queue 1) before the first
      ! device consumer reads the applied u_face/v_face.  apply_bt_correction
      ! reads u_face_x_layer / v_face_y_layer on the default queue, so the
      ! whole coriolis→pgf→hvisc→bdrag→surfstress chain must have landed.
      !$acc wait(1)
      ! accel_visc_rem reweight: `u = u_entry + visc_rem·(u − u_entry)` —
      ! friction-dominated near-massless layers cannot keep a
      ! full-strength explicit dt·F kick (MOM6 MOM_dynamics_split_RK2:
      ! `u = u_init + dt·visc_rem_u·(CAu + PFu + diffu)`; the 2026-07-28
      ! forensics injector-#2 fix).  visc_rem is the LAGGED stage
      ! producer (vdiff_apply_momentum) — ≡ 1 before the first vdiff ⇒
      ! no-op, matching MOM6's semantics.
      if (dyn%accel_visc_rem) then
         call accel_visc_rem_reweight(grid%nx_total + 1, grid%ny_total, ms%nz_ml, &
                                      dyn%avr_u0, dyn%bt_work%visc_rem_u, ms%u_face_x_layer)
         call accel_visc_rem_reweight(grid%nx_total, grid%ny_total + 1, ms%nz_ml, &
                                      dyn%avr_v0, dyn%bt_work%visc_rem_v, ms%v_face_y_layer)
      end if
      call profiler_start("ocean_bt_correction")
      call apply_bt_correction(dyn%bt_work, ms, dt, &
                               skip_h_rescale=is_lagrangian, &
                               grid=grid, &
                               use_bc_pgf=dyn%bt_work%bt_correction_bc_pgf, &
                               use_visc_rem=dyn%bt_work%bt_correction_visc_rem, &
                               metrics=metrics, &
                               scale=merge(dyn%pc_be, 1.0_wp, is_pred), &
                               n_nonfin=dyn%bt_nonfin_step, &
                               n_inner=n_inner)
      if (dyn%bt_nonfin_step > 0) then
         write (output_unit, '("[nan-catch] stage ", i0, " step ", i0, ": ", i0, &
            &" non-finite BT-correction faces (fold skipped — BT loop blown up?)")') &
            stage_id, step_id, dyn%bt_nonfin_step
         flush (output_unit)
      end if
      ! Invariant: ghost face velocities are images of the neighbour's
      ! owned faces before the vertical mixing reads them.  The
      ! visc_rem-weighted fold (and the accel_visc_rem reweight) writes
      ! ghosts with a per-layer weight that is not halo-valid, so refresh
      ! them here.  Configure-time knob, so rank-uniform; a no-op on a
      ! non-periodic single rank.
      if (dyn%bt_work%bt_correction_visc_rem) then
         call ocean_halo_face_x(ms%u_face_x_layer, ms%nz_ml)
         call ocean_halo_face_y(ms%v_face_y_layer, ms%nz_ml)
      end if
      call chksum_state(grid, ms, dyn%chksum_probe, "post_bt_fold", stage_id, step_id)
      ! Land-face velocity reset (C4 / R4b.2): zero re-ingested land-face
      ! layer velocity after the BT correction.  No-op for all-wet.
      call mask_layer_velocities(grid, metrics, ms, bt_work=dyn%bt_work)
      ! Phase-2 vanished-layer velocity reset (R5): runs AFTER apply_bt_correction
      ! so it sees the recombined per-layer velocity.  Only for VCOORD_LAGRANGIAN.
      if (dyn%reset_vanished_u .and. present(vcoord)) then
         if (vcoord%coord_type == VCOORD_LAGRANGIAN) then
            call reset_vanished_layer_velocities(ms, &
                                                 isopycnal_vanish_tol(dyn%angstrom_h, pd_floor=ct%positive_definite))
         end if
      end if
      call profiler_stop("ocean_bt_correction")
      call probe_dS(grid, ms, "after apply_bt_correction", stage_id, step_id)
      call probe_h_vs_eta_residual(grid, ms, dyn%bt_work, stage_id, step_id)

      ! Baroclinic OBC: set per-layer normal velocity at open faces to
      ! Flather mean + zero-gradient baroclinic anomaly (OPEN/TIDAL/
      ! CHAPMAN/NESTED) or uniform clamped_u/v (CLAMPED).  Runs AFTER
      ! apply_bt_correction so it sees the recombined per-layer velocity.
      ! No-op when bc is absent or no edge is open-ish.
      if (present(bc)) then
         call ocean_obc_apply_baroclinic(grid, bc, dyn%bt_work, ms, dt)
         ! The open-edge face values were just rewritten on this tile's
         ! physical rows only; the same faces in a seam ghost (an MPI seam,
         ! OR the local wrap of a single-rank periodic axis) still hold the
         ! pre-OBC velocity, and the boundary-layer scheme below (KPP u*,
         ! shear) reads them before the stage-end exchange.  Refresh
         ! D0-unconditionally, as above.  Collective when decomposed (the
         ! gate is the GLOBAL tags).
         if (ocean_obc_any_open_edge(bc)) then
            call ocean_halo_face_x(ms%u_face_x_layer, ms%nz_ml)
            call ocean_halo_face_y(ms%v_face_y_layer, ms%nz_ml)
         end if
      end if

      ! Sponge.  Map-driven path (&ocean_sponge_nml enable, PR-23)
      ! supersedes the legacy per-edge band; exactly one of the two runs.
      ! Both run AFTER bt-correction so they see the recombined per-layer
      ! velocity. `sponge_maps_on` is a local logical because Fortran does
      ! not guarantee `.and.` short-circuits past `present()`.
      sponge_maps_on = .false.
      if (present(sp)) sponge_maps_on = sp%enable
      if (sponge_maps_on) then
         call ocean_sponge_apply_maps(grid, sp, ms, dt)
      else if (present(bc)) then
         call ocean_sponge_apply(grid, bc, ms, dt)
         call ocean_sponge_apply_tracers(grid, bc, ms, dt)
      end if
      ! Both sponges relax this tile's PHYSICAL cells only; the copies of
      ! those cells in a seam ghost (an MPI seam, OR the local wrap of a
      ! single-rank periodic axis) keep the un-relaxed value, and the
      ! boundary-layer scheme and the corrector's advection read them
      ! before the stage-end exchange.  Refresh D0-unconditionally, NOT
      ! gated on ocean_halo_is_decomposed_x/y(): on a single-rank periodic
      ! run that gate skipped the re-wrap and the stale seam ghost threw
      ! the console `out` term off by orders of magnitude (see
      ! ocean_halo_exchange_ml_state for the no-op / re-wrap / exchange
      ! contract).  Collective when decomposed (the gate is rank-uniform).
      sponge_seam = sponge_maps_on
      if (present(bc)) then
         sponge_seam = sponge_seam .or. &
                       any([bc%west%bc_type, bc%east%bc_type, bc%south%bc_type, &
                            bc%north%bc_type] == OBC_SPONGE)
      end if
      if (sponge_seam) then
         call ocean_halo_face_x(ms%u_face_x_layer, ms%nz_ml)
         call ocean_halo_face_y(ms%v_face_y_layer, ms%nz_ml)
         call refresh_tracer_ghosts(grid, ms, bc=bc)
      end if

      ! ---- 8. Surface tracer fluxes (heat / salt) ----
      ! Same ordering as the unsplit driver: after horizontal +
      ! vertical tracer transport, before vmix/vdiff.  Both `sf`
      ! and `active` are optional — kernel no-ops when either is
      ! absent / false.
      ! Wet/dry: the dynamic cell mask gates the flux (a dry column's
      ! mm-scale sliver must not be heated, plan §4.4).  sw_pen /
      ! restore / geothermal are configure-time excluded with wetdry
      ! (their additive/piston structure mis-composes with a masked
      ! main deposit); knob off ⇒ the original call, byte-identical.
      if (dyn%bt_work%wetdry_enable) then
         call ocean_surface_flux_apply_tracers(grid, sf, ms, therm_dt, &
                                               active=therm_active, &
                                               wet_dyn=dyn%bt_work%wd_wet_dyn)
      else
         call ocean_surface_flux_apply_tracers(grid, sf, ms, therm_dt, active=therm_active)
      end if
      ! Ice-shelf real freshwater MASS -- see the identical block in
      ! `run_stage`.  `chain_weight` is the SAME per-stage weight the
      ! continuity chain's `ocean_accumulate_mass_out` uses (0.5 per
      ! SSP-RK2 stage; 0 / 1 for the pred_corr predictor / corrector), so
      ! the tracked mass source and the tracked boundary outflux are
      ! weighted alike and the console residual closes.
      if (present(cav)) then
         call ocean_cavity_mass_step(grid, metrics, cav, ms, therm_dt, chain_weight, &
                                     active=therm_active)
      end if
      call apply_sw_and_restore(grid, metrics, sf, ms, therm_dt, therm_active)
      call ocean_geothermal_apply_tracers(grid, geo, ms, therm_dt, active=therm_active)
      call probe_dS(grid, ms, "after surface flux", stage_id, step_id)
      ! Ideal-age tracer: interior aging only (1 s/s), thermo-cadence
      ! gated (PR-7).  The surface Dirichlet reset runs once per outer
      ! step, after rk2_average + the ALE remap, in
      ! `ocean_dyn_step_split` — see the ordering comment there.
      !
      ! `.not. is_pred` is what makes the source rate SCHEME-INDEPENDENT.
      ! Under ssp_rk2 this fires on BOTH identical stages and `rk2_average`
      ! turns the two `+therm_dt` into exactly one — the contract
      ! `rdb_ocean_ideal_age`'s header documents.  Under pred_corr there is no
      ! average, so a predictor application would simply survive into the
      ! corrector's and age the ocean at 2x real time (caught by
      ! `split_ideal_age_grows_linearly_via_driver`, which asserts
      ! age == N*DT).  Skipping the predictor leaves ONE application per
      ! outer step, which is the same +therm_dt.  ssp_rk2 is unaffected:
      ! `is_pred` is false on both of its stages.
      call ocean_ideal_age_apply(grid, ms, therm_dt, &
                                 active=dyn%is_thermo_step() .and. .not. is_pred)

      ! ---- 9. Vertical mixing (implicit, stable under any dt) ----
      ! `bt_work=dyn%bt_work` is the visc_rem PRODUCER call: when
      ! `bt_visc_rem_producer` is on (D1 follow-up — decoupled from the
      ! retired `bt_correction_visc_rem`; true whenever forcing/renorm/
      ! bt_rem_from_visc_rem is), this (re)fills `bt_work%visc_rem_u/v`
      ! from THIS stage's momentum solve.  Any consumer that reads it
      ! before the NEXT stage's refresh reads the PREVIOUS stage's γ — a
      ! one-stage (Δt/2) lag, accepted for v1 (see
      ! `vmix_apply_in_stage`'s `bt_work` docstring and PLAN_PR19 §11.1).
      ! At stage 1 of step 1, γ is still at its `source=1.0` init, so the
      ! very first correction is h-only ⇒ benign.
      ! pred_corr predictor: MOM6 applies the implicit friction to the
      ! provisional velocity too — `vertvisc(up, dt_pred)` (SPEC §1 fact 8's
      ! coef-only claim is WRONG) — so the
      ! spurious grounded-layer ballistics are absorbed BEFORE the
      ! predictor continuity forms u_av.  Momentum-only there (tracers
      ! untouched); dt_vel = BE·dt matches MOM6's dt_pred.
      if (is_pred) then
         ! PR-1: thread `dt_remnant=dt` so the visc_rem PRODUCER (when
         ! `bt_visc_rem_producer` is on) is built at the outer step's
         ! full `dt`, NOT the predictor's own `dt_vel = pc_be·dt` —
         ! MOM6's `VISC_REM_TIMESTEP_BUG = .false.` default.  `bc` is
         ! forwarded so the split remnant-only refresh can re-wrap the ghosts.
         if (fold_top) then
            call vmix_apply_in_stage(grid, dyn, vmix, vd, ss, bd, ms, dt_vel, stage, sf, epbl=epbl, &
                                     kshear=kshear, vmix_tidal=vmix_tidal, bt_work=dyn%bt_work, &
                                     apply_tracers=.false., metrics=metrics, &
                                     lambda_top_u=td%lambda_top_u, lambda_top_v=td%lambda_top_v, &
                                     cover_u=td%cover_u, cover_v=td%cover_v, &
                                     dt_remnant=dt, bc=bc)
         else
            call vmix_apply_in_stage(grid, dyn, vmix, vd, ss, bd, ms, dt_vel, stage, sf, epbl=epbl, &
                                     kshear=kshear, vmix_tidal=vmix_tidal, bt_work=dyn%bt_work, &
                                     apply_tracers=.false., metrics=metrics, &
                                     dt_remnant=dt, bc=bc)
         end if
      else
         if (fold_top) then
            call vmix_apply_in_stage(grid, dyn, vmix, vd, ss, bd, ms, dt, stage, sf, epbl=epbl, &
                                     kshear=kshear, vmix_tidal=vmix_tidal, bt_work=dyn%bt_work, &
                                     metrics=metrics, &
                                     lambda_top_u=td%lambda_top_u, lambda_top_v=td%lambda_top_v, &
                                     cover_u=td%cover_u, cover_v=td%cover_v, bc=bc)
         else
            call vmix_apply_in_stage(grid, dyn, vmix, vd, ss, bd, ms, dt, stage, sf, epbl=epbl, &
                                     kshear=kshear, vmix_tidal=vmix_tidal, bt_work=dyn%bt_work, &
                                     metrics=metrics, bc=bc)
         end if
      end if
      ! KE attribution: implicit vertical friction (+ folded drag/stress
      ! when implicit_*) — the stage-close segment (debug_ke_attr).
      call ke_probe_sample(grid, ms, dyn%ke_probe, "vdiff_vmix", stage_id, step_id)
      call chksum_state(grid, ms, dyn%chksum_probe, "stage_close", stage_id, step_id)
      call probe_dS(grid, ms, "after vdiff/vmix (stage end)", stage_id, step_id)

      ! pred_corr: the continuity + tracer chain runs HERE — after the
      ! single prognostic update (corrector) / provisional up (predictor)
      ! and after the implicit friction — so the thickness advances with
      ! the UPDATED velocities via uhbt + u_cor.  This is the
      ! forward-backward gravity-wave pairing that lifts the
      ! internal-wave dt ceiling (SPEC §1 fact 5, §2 C8).
      if (is_pc) then
         ! The surface fluxes and the vertical mixing above updated every
         ! tracer column, ghosts included, but a ghost column's diffusivity
         ! is computed from the TILE's data and is not the neighbour's
         ! interior value wherever the boundary layer / smoothing stencil
         ! reaches past the ghost band — so the ghosts the corrector's
         ! tracer advection reads next are no longer images of the
         ! neighbour.  Serial, that happened only at a periodic wrap; on a
         ! decomposed run it happened at every seam, so the answer
         ! depended on where the seam fell (1 ULP in hTr at a seam cell by
         ! step 12 of the spherical 2x1 case, growing from there).  The
         ! `ssp_rk2` chain runs straight after the stage-entry exchange
         ! and never sees this.  Refresh the tracer ghosts (exchange, local
         ! periodic wrap, north fold) first.
         call refresh_tracer_ghosts(grid, ms, bc)
         call run_continuity_chain(grid, metrics, dyn, ct, hd, va, redi, varmix, ms, &
                                   dt, therm_dt, therm_active, is_lagrangian, &
                                   h_min_floor, chain_weight, is_pred, &
                                   stage_id, step_id, bc=bc, mle=mle)
      end if

      ! ---- 10. Stage-end seam reconciliation (tripolar only) ----
      ! The velocity applies + bt-correction + sponge updated v/u/h/tracers
      ! AFTER the post-continuity fold; the duplicated-DOF v/corner row must
      ! be re-projected and the north ghosts re-folded so the stage output
      ! is fully seam-consistent (the v on-row antisymmetric projection is
      ! "applied after each update of v" — Appendix A).  Periodic-first-
      ! fold-second: re-wrap periodic ghosts, then fold.  No-op when not
      ! folding (and when no periodic axis is set the periodic wrap no-ops too).
      ! O2: unconditional multi-rank ghost exchange (D0 — no-op on 1 rank).
      call ocean_halo_exchange_ml_state(ms)
      if (present(bc)) then
         if (bc%north_fold) 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
      end if
   end subroutine run_stage_split