subroutine run_continuity_chain(grid, metrics, dyn, ct, hd, va, redi, varmix, ms, &
dt, therm_dt, therm_active, is_lagrangian, &
h_min_floor, mass_out_weight, h_only, &
stage_id, step_id, bc, mle)
!! The slow horizontal continuity + tracer chain (ghost fills →
!! constrained continuity+tracer split → reservoirs → halo/wrap →
!! tracer hdiff → Redi → vertical advection), extracted verbatim
!! from `run_stage_split` so the pred_corr path can run it AFTER the
!! velocity update + implicit friction (the forward-backward
!! pairing, SPEC §2 C8/§1 fact 5) while the historical ssp_rk2
!! path keeps it before the applies (bit-identical).
!! `h_only` selects the predictor's TR_MODE_NONE (SPEC §2 P9);
!! `mass_out_weight` is the per-call budget weight (0.5 per SSP
!! stage; 0 for the discarded predictor state, 1 for the
!! corrector).
type(hgrid_t), intent(in) :: grid
type(ocean_metrics_t), intent(in) :: metrics
type(ocean_dyn_t), intent(inout) :: dyn
type(continuity_t), intent(inout) :: ct
type(ocean_hdiff_tracer_t), intent(inout) :: hd
type(ocean_vertical_advection_t), intent(inout) :: va
type(ocean_redi_t), intent(inout), optional :: redi
type(ocean_varmix_t), intent(inout), optional :: varmix
type(multilayer_state_t), intent(inout) :: ms
real(wp), intent(in) :: dt, therm_dt
logical, intent(in) :: therm_active, is_lagrangian
real(wp), intent(in) :: h_min_floor, mass_out_weight
logical, intent(in) :: h_only
integer, intent(in) :: stage_id, step_id
type(ocean_bc_state_t), intent(inout), optional :: bc
type(ocean_mle_t), intent(inout), optional :: mle
integer :: tr_mode
real(wp) :: h_min_pass
! ---- 5. Ghost fills at open edges (before continuity PPM) ----
! Fill h_layer and tracer hTr ghost cells at any open-ish edge
! using zero-gradient (h) and upwind-aware (hTr) logic.
! This generalises the CLAMPED-only ghost fill inside
! rdb_continuity::tracer_advect_zonal/meridional to cover OPEN /
! TIDAL / CHAPMAN / NESTED as well.
!
! Authority split (no conflict): on CLAMPED edges the old fill inside
! rdb_continuity runs LAST (once per Lie-split sub-flux, after h_layer
! has been updated mid-split) and is the authority — it unconditionally
! clamps the ghost to clamped_tracer*h, matching what the BT mode
! imposes. This fill's CLAMPED writes here are simply overwritten by
! it. On the open-class edges (OPEN/TIDAL/CHAPMAN/NESTED) the old fill
! is inert, so THIS fill's upwind-aware values survive and govern.
! No-op when bc is absent or all edges are WALL.
if (present(bc)) then
call ocean_obc_fill_ghosts(grid, bc, ms)
! The fill writes the open-edge ghost rows/columns over this tile's
! PHYSICAL span only; the corner cells beyond a seam (an MPI seam,
! OR the local wrap of a single-rank periodic axis) x an open-edge
! ghost row are someone else's fill, which only a refresh delivers,
! and the Lie-split advection's first pass reads them. So refresh
! D0-unconditionally, NOT gated on ocean_halo_is_decomposed_x/y():
! the halo calls no-op on a single-rank non-periodic axis, re-wrap
! a single-rank periodic one and exchange when decomposed.
! Collective when decomposed (the gate is the GLOBAL tags).
if (ocean_obc_any_open_edge(bc)) then
call ocean_halo_centre(ms%h_layer, ms%nz_ml)
call refresh_tracer_ghosts(grid, ms, bc=bc)
end if
end if
! ---- 5b. Slow horizontal continuity + tracer, constrained ----
! Pass `bt_uhbt, bt_vhbt` so the per-layer mass fluxes are
! renormalised to vertically sum to the barotropic-substep's transport.
! After apply_zonal + apply_meridional, `sum_k(h_layer)` equals
! `H_ref + bt_eta_end` to machine precision — no h rescale
! needed afterwards. Tracer rides the same renormalised
! fluxes ⇒ per-column `T = hTr/h` stays exact.
call profiler_start("ocean_continuity")
! `mle_fold_active = is_thermo_step()` gates the Fox-Kemper fold to
! the thermo cadence: the FK transports are computed once per thermo
! interval (mle_compute_transports, above), so folding them on the
! intervening non-thermo steps would re-apply stale transports. At
! the default dt_therm_ratio = 1 this is .true. every step ⇒ bit-identical.
!
! Phase 2 (6b) horizontal-tracer-advect cadence: ratio == 1 ⇒
! TR_MODE_ADVECT (fused every-step advect, bit-identical); ratio > 1
! ⇒ TR_MODE_ACCUMULATE (advance h + accumulate ½·flux·dt into
! ct%uhtr/vhtr per RK2 stage, hTr frozen). The boundary drain runs
! once per window in ocean_dyn_step_split.
if (h_only) then
! pred_corr PREDICTOR (SPEC §2 P9): continuity advances h (into the
! provisional hp, discarded after the stage) and produces u_av via
! u_cor, but must not touch tracers or the advect window.
tr_mode = TR_MODE_NONE
else if (dyn%dt_tracer_advect_ratio <= 1) then
tr_mode = TR_MODE_ADVECT
else
tr_mode = TR_MODE_ACCUMULATE
end if
! Conservative min-thickness mode: pass h_min=0 so continuity performs
! the RAW (non-injecting) h-update; the sub-floor layers are then repaired
! conservatively by ocean_apply_conservative_min_thickness below. The
! legacy injecting floor stays the default (h_min=h_min_floor).
h_min_pass = h_min_floor
if (ct%conservative_floor) h_min_pass = 0.0_wp
! MOM6 time-mean fields (SPEC §2 C7/C9): stash h_in in h_av, average with
! h_out after continuity. NOTE: under the present two-stage scheme this
! yields a PER-STAGE mean, not the per-step mean MOM6 forms; it becomes
! exact when S4 restructures the outer loop. Nothing reads it until S3.
call copy_field_3d(ms%h_layer, ms%h_av_layer, &
size(ms%h_layer, 1), size(ms%h_layer, 2), size(ms%h_layer, 3))
! SPEC S2b (`&ocean_bt_nml renorm_visc_rem`): forward `visc_rem_u/v`
! so the transport-matching inversion runs in the γ-weighted MOM6
! form (`u_cor = u + du·γ_k`, Jacobian `dy·h_marg·γ_k`) — a
! heavily-frictioned layer receives a smaller share of `du`. Off ⇒
! the historical uniform-`du` renormaliser, bit-identical.
if (dyn%bt_work%bt_renorm_visc_rem) then
if (present(bc)) then
call continuity_tracer_step_split(grid, metrics, ct, ms, dt, &
uhbt=dyn%bt_work%bt_uhbt, &
vhbt=dyn%bt_work%bt_vhbt, bc=bc, mle=mle, &
mle_fold_active=dyn%is_thermo_step(), &
tracer_mode=tr_mode, &
h_min=h_min_pass, &
visc_rem_u=dyn%bt_work%visc_rem_u, &
visc_rem_v=dyn%bt_work%visc_rem_v, &
u_cor=ms%u_av_layer, v_cor=ms%v_av_layer)
else
call continuity_tracer_step_split(grid, metrics, ct, ms, dt, &
uhbt=dyn%bt_work%bt_uhbt, &
vhbt=dyn%bt_work%bt_vhbt, mle=mle, &
mle_fold_active=dyn%is_thermo_step(), &
tracer_mode=tr_mode, &
h_min=h_min_pass, &
visc_rem_u=dyn%bt_work%visc_rem_u, &
visc_rem_v=dyn%bt_work%visc_rem_v, &
u_cor=ms%u_av_layer, v_cor=ms%v_av_layer)
end if
else if (present(bc)) then
call continuity_tracer_step_split(grid, metrics, ct, ms, dt, &
uhbt=dyn%bt_work%bt_uhbt, &
vhbt=dyn%bt_work%bt_vhbt, bc=bc, mle=mle, &
mle_fold_active=dyn%is_thermo_step(), &
tracer_mode=tr_mode, &
h_min=h_min_pass, &
u_cor=ms%u_av_layer, v_cor=ms%v_av_layer)
else
call continuity_tracer_step_split(grid, metrics, ct, ms, dt, &
uhbt=dyn%bt_work%bt_uhbt, &
vhbt=dyn%bt_work%bt_vhbt, mle=mle, &
mle_fold_active=dyn%is_thermo_step(), &
tracer_mode=tr_mode, &
h_min=h_min_pass, &
u_cor=ms%u_av_layer, v_cor=ms%v_av_layer)
end if
! MOM6 C9: h_av = 0.5*(h_in + h_out).
call rk2_average_field_3d(ms%h_layer, ms%h_av_layer, &
size(ms%h_layer, 1), size(ms%h_layer, 2), size(ms%h_layer, 3))
! Conservative minimum-thickness borrow (isopycnal grounding stability).
! Only when the knob is on AND we have a positive floor (⇒ VCOORD_LAGRANGIAN
! by the setup fail-loud guard). No-op on any ungrounded column ⇒ the
! isopycnal interface structure is untouched where all layers meet floor.
! NOT checked under `conservative_floor`: that mode deliberately passes
! `h_min = 0`, and the three clamp sites are gated on `h_min_use > 0`, so
! the h-update here is genuinely RAW — a transiently negative layer is
! legal and is repaired by the borrow below (which flags any
! `h < h_floor`, negatives included, and inflates it back). Checking here
! would abort on the first benign intermediate. With the injecting floor
! the clamp DOES hold, so a negative there is a real defect.
if (dyn%check_h_positive) then
call check_h_positive_or_die(grid, ms, "after continuity_tracer_step_split", &
0, dyn%outer_step_count + 1, &
check_layers=.not. ct%conservative_floor)
end if
if (ct%conservative_floor .and. h_min_floor > 0.0_wp) then
call profiler_start("ocean_min_thickness")
call ocean_apply_conservative_min_thickness(grid, ms, ct%mt_h_new%data, &
ct%mt_grounded%data, h_min_floor)
call profiler_stop("ocean_min_thickness")
if (dyn%check_h_positive) then
call check_h_positive_or_die(grid, ms, "after conservative_min_thickness", &
0, dyn%outer_step_count + 1, check_layers=.true.)
end if
end if
! Mass budget: flux_h_layer now holds the total horizontal divergence;
! accumulate the boundary outflux for this RK2 stage (weight 0.5).
call ocean_accumulate_mass_out(ms, ms%flux_h_layer, metrics%areaT, &
grid%nghost, dt, mass_out_weight)
call profiler_stop("ocean_continuity")
call probe_dS(grid, ms, "after continuity_tracer_split", stage_id, step_id)
! Reservoir update (§1, v2): evolve tres toward T_int / T_data using
! the stage's wall-face mass_flux_*_layer values. Called while those
! arrays still hold this stage's fluxes (before any overwrite).
! No-op when res_lscale_out == 0 and res_lscale_in == 0 (default).
if (present(bc)) call ocean_obc_update_reservoirs(grid, bc, ms, dt)
! Post-continuity periodic wrap (design §1.5): re-wrap h_layer and
! tracers after the Lie-split continuity step and before hdiff, so
! the horizontal diffusion kernel reads correct ghost-zone values.
! No-op when neither periodic axis 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 after the post-continuity periodic wrap, so
! hdiff reads fold-consistent ghost values at the seam.
if (present(bc)) call ocean_fold_wrap_state(grid, bc, ms)
call profiler_start("ocean_tracer_hdiff")
call tracer_hdiff(grid, metrics, hd, ms, therm_dt, active=therm_active, bc=bc)
call profiler_stop("ocean_tracer_hdiff")
call probe_dS(grid, ms, "after tracer_hdiff", stage_id, step_id)
! Redi neutral diffusion (capability [3]) — Phase B. AUGMENTS the
! along-coordinate `tracer_hdiff` above with the rotated (along-
! isopycnal) flux, using the Phase-A coefficients built once this
! outer step. THERMO cadence (self-gated inside on `therm_active`
! via the same `is_thermo_step` window the coeffs were built under).
! No-op when absent / disabled / khtr<=0.
if (present(redi) .and. therm_active) then
call profiler_start("ocean_redi")
! VarMix seam: feed Redi the spatially-varying KhTr when VarMix is
! enabled; otherwise Redi falls back to its scalar khtr.
if (present(varmix)) then
if (varmix%enable) then
call redi_apply_flux(grid, metrics, redi, ms, therm_dt, &
khtr_u_ext=varmix%khtr_u, khtr_v_ext=varmix%khtr_v, bc=bc)
else
call redi_apply_flux(grid, metrics, redi, ms, therm_dt, bc=bc)
end if
else
call redi_apply_flux(grid, metrics, redi, ms, therm_dt, bc=bc)
end if
call profiler_stop("ocean_redi")
end if
! ---- Vertical advection (Eulerian-z only) ----
! In Eulerian-z mode the vertical w-divergence cancels the
! horizontal flux divergence per layer, pinning `h_layer` to
! reference. In Lagrangian mode (`vcoord` present and
! coord_type /= EULERIAN_Z) we skip this — h_layer is allowed
! to evolve under the MOM6-constrained horizontal continuity
! alone, and the ALE remap at end of outer step relayers
! conservatively onto `vcoord%target_h`. Skipping here
! removes the surface F(nz+1) = 0 vs w(nz+1) ≠ 0 CWC
! inconsistency that seeds the salt-baroclinic instability
! under realistic β_S.
if (.not. is_lagrangian) then
call profiler_start("ocean_vertical_advect")
call compute_w_from_continuity(grid, va, ms)
call probe_dS(grid, ms, "after compute_w_from_cont", stage_id, step_id)
call tracer_advect_vertical(grid, va, ms, therm_dt, active=therm_active)
call probe_dS(grid, ms, "after tracer_advect_vertical", stage_id, step_id)
call profiler_stop("ocean_vertical_advect")
end if
end subroutine run_continuity_chain