rdb_continuity Module

Holds scheme-variant flags and reusable workspace for the C-grid continuity-PPM (Lin & Rood) thickness-flux kernel, plus the public step routines that the split-explicit driver calls per barotropic substep. The kernel is the primary mass-flux producer for the ocean path: it consumes face velocities + cell-centred thickness and emits per-face mass fluxes that the tracer-advection kernels then reuse under CWC for free.

Phase 2b: method-of-lines PPM face reconstruction (Colella & Woodward 1984, eq 1.6 face value + eq 1.10 monotonic limiter) on the barotropic C-grid state. Closed-wall BC. Cells too close to a wall (< 2 cells from the boundary) fall back to first-order (h_L = h_R = h_centre); the resulting kernel preserves uniform fields bit-for-bit (the constancy-preservation property the lake-at-rest tests guard) and propagates Gaussian humps with < 5% peak diffusion over their own width.


Uses

  • module~~rdb_continuity~~UsesGraph module~rdb_continuity rdb_continuity iso_fortran_env iso_fortran_env module~rdb_continuity->iso_fortran_env module~rdb_barotropic_state rdb_barotropic_state module~rdb_continuity->module~rdb_barotropic_state module~rdb_constants rdb_constants module~rdb_continuity->module~rdb_constants module~rdb_grid rdb_grid module~rdb_continuity->module~rdb_grid module~rdb_mem_report rdb_mem_report module~rdb_continuity->module~rdb_mem_report module~rdb_multilayer_state rdb_multilayer_state module~rdb_continuity->module~rdb_multilayer_state module~rdb_ocean_boundary_types rdb_ocean_boundary_types module~rdb_continuity->module~rdb_ocean_boundary_types module~rdb_ocean_fold rdb_ocean_fold module~rdb_continuity->module~rdb_ocean_fold module~rdb_ocean_fold_apply rdb_ocean_fold_apply module~rdb_continuity->module~rdb_ocean_fold_apply module~rdb_ocean_fold_exchange rdb_ocean_fold_exchange module~rdb_continuity->module~rdb_ocean_fold_exchange module~rdb_ocean_gm rdb_ocean_gm module~rdb_continuity->module~rdb_ocean_gm module~rdb_ocean_halo rdb_ocean_halo module~rdb_continuity->module~rdb_ocean_halo module~rdb_ocean_metrics rdb_ocean_metrics module~rdb_continuity->module~rdb_ocean_metrics module~rdb_ocean_mle rdb_ocean_mle module~rdb_continuity->module~rdb_ocean_mle module~rdb_ocean_periodic rdb_ocean_periodic module~rdb_continuity->module~rdb_ocean_periodic module~rdb_profiler rdb_profiler module~rdb_continuity->module~rdb_profiler module~rdb_recon_weno rdb_recon_weno module~rdb_continuity->module~rdb_recon_weno module~rdb_scratch_3d rdb_scratch_3d module~rdb_continuity->module~rdb_scratch_3d module~rdb_tracer rdb_tracer module~rdb_continuity->module~rdb_tracer module~rdb_barotropic_state->iso_fortran_env module~rdb_barotropic_state->module~rdb_constants module~rdb_barotropic_state->module~rdb_grid module~rdb_barotropic_state->module~rdb_mem_report pic_types pic_types module~rdb_constants->pic_types module~rdb_grid->module~rdb_constants module~rdb_mem_report->iso_fortran_env module~rdb_mem_report->module~rdb_constants pic_logger pic_logger module~rdb_mem_report->pic_logger pic_strings pic_strings module~rdb_mem_report->pic_strings module~rdb_multilayer_state->iso_fortran_env module~rdb_multilayer_state->module~rdb_constants module~rdb_multilayer_state->module~rdb_grid module~rdb_multilayer_state->module~rdb_mem_report module~rdb_multilayer_state->module~rdb_tracer module~rdb_efp rdb_efp module~rdb_multilayer_state->module~rdb_efp module~rdb_error_ring rdb_error_ring module~rdb_multilayer_state->module~rdb_error_ring module~rdb_multilayer_state->pic_logger module~rdb_ocean_boundary_types->iso_fortran_env module~rdb_ocean_boundary_types->module~rdb_constants module~rdb_ocean_boundary_types->module~rdb_grid module~rdb_ocean_boundary_types->module~rdb_mem_report module~rdb_ocean_status rdb_ocean_status module~rdb_ocean_boundary_types->module~rdb_ocean_status module~rdb_ocean_tide_astro rdb_ocean_tide_astro module~rdb_ocean_boundary_types->module~rdb_ocean_tide_astro pic_ascii pic_ascii module~rdb_ocean_boundary_types->pic_ascii module~rdb_ocean_boundary_types->pic_logger module~rdb_ocean_fold->module~rdb_constants module~rdb_ocean_fold_apply->module~rdb_constants module~rdb_ocean_fold_apply->module~rdb_grid module~rdb_ocean_fold_apply->module~rdb_multilayer_state module~rdb_ocean_fold_apply->module~rdb_ocean_boundary_types module~rdb_ocean_fold_apply->module~rdb_ocean_fold module~rdb_ocean_fold_apply->module~rdb_ocean_fold_exchange module~rdb_ocean_fold_exchange->module~rdb_constants module~rdb_ocean_fold_exchange->module~rdb_ocean_fold module~rdb_comm_env rdb_comm_env module~rdb_ocean_fold_exchange->module~rdb_comm_env module~rdb_decomp rdb_decomp module~rdb_ocean_fold_exchange->module~rdb_decomp module~rdb_ocean_fold_exchange->module~rdb_error_ring module~rdb_ocean_fold_plan rdb_ocean_fold_plan module~rdb_ocean_fold_exchange->module~rdb_ocean_fold_plan module~rdb_ocean_fold_exchange->module~rdb_ocean_status module~rdb_ocean_fold_exchange->pic_logger pic_mpi_lib pic_mpi_lib module~rdb_ocean_fold_exchange->pic_mpi_lib module~rdb_ocean_fold_exchange->pic_strings module~rdb_ocean_gm->iso_fortran_env module~rdb_ocean_gm->module~rdb_constants module~rdb_ocean_gm->module~rdb_grid module~rdb_ocean_gm->module~rdb_mem_report module~rdb_ocean_gm->module~rdb_multilayer_state module~rdb_ocean_gm->module~rdb_ocean_metrics ieee_arithmetic ieee_arithmetic module~rdb_ocean_gm->ieee_arithmetic module~rdb_ocean_isopycnal_slopes rdb_ocean_isopycnal_slopes module~rdb_ocean_gm->module~rdb_ocean_isopycnal_slopes module~rdb_ocean_halo->module~rdb_constants module~rdb_ocean_halo->module~rdb_ocean_periodic module~rdb_ocean_halo->module~rdb_comm_env module~rdb_ocean_halo->module~rdb_decomp module~rdb_ocean_halo->module~rdb_error_ring module~rdb_ocean_halo_counters rdb_ocean_halo_counters module~rdb_ocean_halo->module~rdb_ocean_halo_counters module~rdb_ocean_halo->module~rdb_ocean_status module~rdb_ocean_halo->pic_logger module~rdb_ocean_halo->pic_mpi_lib module~rdb_ocean_halo->pic_strings module~rdb_ocean_metrics->iso_fortran_env module~rdb_ocean_metrics->module~rdb_constants module~rdb_ocean_metrics->module~rdb_grid module~rdb_ocean_metrics->module~rdb_mem_report module~rdb_ocean_metrics->module~rdb_ocean_fold module~rdb_ocean_metrics->module~rdb_error_ring module~rdb_io_netcdf rdb_io_netcdf module~rdb_ocean_metrics->module~rdb_io_netcdf module~rdb_ocean_bipolar rdb_ocean_bipolar module~rdb_ocean_metrics->module~rdb_ocean_bipolar module~rdb_ocean_metrics->module~rdb_ocean_status netcdf netcdf module~rdb_ocean_metrics->netcdf module~rdb_ocean_metrics->pic_logger module~rdb_ocean_metrics->pic_strings module~rdb_ocean_mle->iso_fortran_env module~rdb_ocean_mle->module~rdb_constants module~rdb_ocean_mle->module~rdb_grid module~rdb_ocean_mle->module~rdb_mem_report module~rdb_ocean_mle->module~rdb_multilayer_state module~rdb_ocean_mle->module~rdb_ocean_boundary_types module~rdb_ocean_mle->module~rdb_ocean_metrics module~rdb_ocean_epbl rdb_ocean_epbl module~rdb_ocean_mle->module~rdb_ocean_epbl module~rdb_ocean_surface_stress rdb_ocean_surface_stress module~rdb_ocean_mle->module~rdb_ocean_surface_stress module~rdb_ocean_periodic->module~rdb_constants module~rdb_ocean_periodic->module~rdb_grid module~rdb_ocean_periodic->module~rdb_multilayer_state module~rdb_ocean_periodic->module~rdb_ocean_boundary_types module~rdb_profiler->iso_fortran_env module~rdb_profiler->pic_logger module~rdb_recon_weno->module~rdb_constants module~rdb_recon_weno->module~rdb_grid module~rdb_scratch_3d->iso_fortran_env module~rdb_scratch_3d->module~rdb_constants module~rdb_scratch_3d->module~rdb_mem_report module~rdb_tracer->iso_fortran_env module~rdb_tracer->module~rdb_constants module~rdb_tracer->module~rdb_grid module~rdb_tracer->module~rdb_mem_report module~rdb_comm_env->iso_fortran_env module~rdb_comm_env->module~rdb_constants module~rdb_comm_env->pic_mpi_lib module~rdb_config rdb_config module~rdb_decomp->module~rdb_config module~rdb_efp->iso_fortran_env module~rdb_efp->ieee_arithmetic module~rdb_error_ring->pic_logger module~rdb_io_netcdf->iso_fortran_env module~rdb_io_netcdf->module~rdb_constants module~rdb_io_netcdf->module~rdb_error_ring module~rdb_io_netcdf->netcdf module~rdb_io_netcdf->pic_logger module~rdb_io_netcdf->pic_strings iso_c_binding iso_c_binding module~rdb_io_netcdf->iso_c_binding module~rdb_ocean_bipolar->module~rdb_constants module~rdb_ocean_epbl->iso_fortran_env module~rdb_ocean_epbl->module~rdb_constants module~rdb_ocean_epbl->module~rdb_grid module~rdb_ocean_epbl->module~rdb_mem_report module~rdb_ocean_epbl->module~rdb_multilayer_state module~rdb_ocean_epbl->module~rdb_scratch_3d module~rdb_ocean_epbl->module~rdb_ocean_surface_stress module~rdb_eos rdb_eos module~rdb_ocean_epbl->module~rdb_eos module~rdb_ocean_surface_flux rdb_ocean_surface_flux module~rdb_ocean_epbl->module~rdb_ocean_surface_flux module~rdb_ocean_halo_counters->iso_fortran_env module~rdb_ocean_halo_counters->pic_strings module~rdb_ocean_isopycnal_slopes->iso_fortran_env module~rdb_ocean_isopycnal_slopes->module~rdb_constants module~rdb_ocean_isopycnal_slopes->module~rdb_grid module~rdb_ocean_isopycnal_slopes->module~rdb_mem_report module~rdb_ocean_isopycnal_slopes->module~rdb_multilayer_state module~rdb_ocean_isopycnal_slopes->module~rdb_ocean_metrics module~rdb_ocean_isopycnal_slopes->module~rdb_eos module~rdb_ocean_surface_stress->iso_fortran_env module~rdb_ocean_surface_stress->module~rdb_constants module~rdb_ocean_surface_stress->module~rdb_grid module~rdb_ocean_surface_stress->module~rdb_mem_report module~rdb_ocean_surface_stress->module~rdb_multilayer_state module~rdb_ocean_surface_stress->module~rdb_scratch_3d module~rdb_ocean_tide_astro->iso_fortran_env module~rdb_ocean_tide_astro->module~rdb_constants module~rdb_config->module~rdb_constants module~rdb_config->module~rdb_error_ring module~rdb_config->module~rdb_ocean_status module~rdb_config->pic_ascii module~rdb_config->pic_logger module~rdb_config->pic_strings module~rdb_ice_enthalpy rdb_ice_enthalpy module~rdb_config->module~rdb_ice_enthalpy module~rdb_ice_init rdb_ice_init module~rdb_config->module~rdb_ice_init module~rdb_nml_schema rdb_nml_schema module~rdb_config->module~rdb_nml_schema module~rdb_eos->module~rdb_constants module~rdb_eos->module~rdb_grid module~rdb_ocean_surface_flux->iso_fortran_env module~rdb_ocean_surface_flux->module~rdb_constants module~rdb_ocean_surface_flux->module~rdb_grid module~rdb_ocean_surface_flux->module~rdb_mem_report module~rdb_ocean_surface_flux->module~rdb_multilayer_state module~rdb_ice_enthalpy->module~rdb_constants module~rdb_ice_init->module~rdb_constants module~rdb_ice_init->module~rdb_grid module~rdb_ice_init->module~rdb_multilayer_state module~rdb_ice_init->module~rdb_ocean_metrics module~rdb_ice_init->module~rdb_ice_enthalpy module~rdb_ice_column rdb_ice_column module~rdb_ice_init->module~rdb_ice_column module~rdb_ice_state rdb_ice_state module~rdb_ice_init->module~rdb_ice_state module~rdb_nml_schema->module~rdb_constants module~rdb_nml_schema->module~rdb_error_ring module~rdb_nml_schema->pic_logger

Used by

  • module~~rdb_continuity~~UsedByGraph module~rdb_continuity rdb_continuity module~rdb_ice_transport rdb_ice_transport module~rdb_ice_transport->module~rdb_continuity module~rdb_ocean_dyn rdb_ocean_dyn module~rdb_ocean_dyn->module~rdb_continuity module~rdb_ocean_state rdb_ocean_state module~rdb_ocean_state->module~rdb_continuity module~rdb_ocean_state->module~rdb_ocean_dyn module~rdb_driver rdb_driver module~rdb_driver->module~rdb_ocean_dyn module~rdb_driver->module~rdb_ocean_state module~rdb_ocean_engine rdb_ocean_engine module~rdb_driver->module~rdb_ocean_engine module~rdb_handle rdb_handle module~rdb_handle->module~rdb_ocean_state module~rdb_handle->module~rdb_ocean_engine module~rdb_ocean_api rdb_ocean_api module~rdb_ocean_api->module~rdb_ocean_dyn module~rdb_ocean_api->module~rdb_handle module~rdb_ocean_diag_derived rdb_ocean_diag_derived module~rdb_ocean_api->module~rdb_ocean_diag_derived module~rdb_ocean_diag_fills rdb_ocean_diag_fills module~rdb_ocean_api->module~rdb_ocean_diag_fills module~rdb_ocean_api->module~rdb_ocean_engine module~rdb_ocean_diag_derived->module~rdb_ocean_state module~rdb_ocean_diag_derived->module~rdb_ocean_diag_fills module~rdb_ocean_diag_fills->module~rdb_ocean_state module~rdb_ocean_engine->module~rdb_ice_transport module~rdb_ocean_engine->module~rdb_ocean_dyn module~rdb_ocean_engine->module~rdb_ocean_state module~rdb_ocean_engine->module~rdb_ocean_diag_derived module~rdb_ocean_engine->module~rdb_ocean_diag_fills module~rdb_ocean_setup rdb_ocean_setup module~rdb_ocean_engine->module~rdb_ocean_setup module~rdb_ocean_setup->module~rdb_ocean_dyn module~rdb_ocean_setup->module~rdb_ocean_state

Variables

Type Visibility Attributes Name Initial
integer, public, parameter :: TR_MODE_ACCUMULATE = 1
integer, public, parameter :: TR_MODE_ADVECT = 0

Every-step fused mode: advance h AND advect tracers each call (the historical default; ratio = 1 bypass).
integer, public, parameter :: TR_MODE_NONE = 2

h-only continuity: skip tracer advection AND window accumulation entirely. Used by the pred_corr predictor (SPEC §2 P9 — MOM6’s predictor continuity advances hp but never touches tracers; the predictor state is discarded except for u_av / h_av). Windowed mode (ratio > 1): advance h, accumulate 0.5·mass_flux·dt into uhtr/vhtr, SKIP the per-step tracer advect (hTr held frozen until the boundary drain).
integer, private, parameter :: AVAIL_LIMIT_PASS = 8
real(kind=wp), private, parameter :: DRAIN_BUDGET_IN_STAGE_WEIGHT = 1.0_wp
real(kind=wp), private, parameter :: DRAIN_BUDGET_POST_AVERAGE_WEIGHT = 2.0_wp
real(kind=wp), private, parameter :: DRAIN_MIN_H = 1.0e-10_wp
real(kind=wp), private, parameter :: DRAIN_MIN_VOL = 1.0e-10_wp
real(kind=wp), private, parameter :: RENORM_CFL = 0.25_wp

CFL cap on the reconciliation correction (MOM6 CONTINUITY_CFL_LIMIT, default 0.5). MOM6 brackets its Newton solve by this and will ACCEPT a residual uhbt mismatch rather than hand any layer a super-CFL velocity. Without the bracket the solve can satisfy the transport constraint by assigning a huge du to a near-massless layer – which is exactly the grounded sliver case here.
integer, private, parameter :: RENORM_MAXIT = 8

Newton iterations for the Sum_k uh_k == uhbt reconciliation. The flux is a NONLINEAR function of the correction du because the donor cell is picked by sign(u): adding du can flip a layer’s velocity sign, changing its h_face, so a single linear step does not actually land on uhbt. MOM6 solves the same equation with Newton + bisection (zonal_flux_adjust, up to 20 its). When no donor flips – the overwhelmingly common case – iteration 2 sees zero residual and exits, so the result is BIT-IDENTICAL to the previous single-step form.
integer, private, parameter :: RENORM_MAXIT_CONSISTENT = 20

Iteration cap when renorm_consistent_flux is on (MOM6 zonal_flux_adjust also allows 20). Newton lands on a single-kink root in <= 3 iterations; the head-room is for the bisection fallback on a many-layer face with several donor flips.
real(kind=wp), private, parameter :: RENORM_TOL = 1.0e-12_wp

Relative convergence tolerance on |uhbt - Sum_k uh_k|.
real(kind=wp), private, parameter :: RENORM_VR_MIN = 1.0e-12_wp

Floor below which a layer’s viscous remnant γ_k is treated as zero when bracketing du — such a layer receives no barotropic increment, so it constrains nothing and must not divide the CFL bound.

Derived Types

type, public ::  continuity_t

Components

Type Visibility Attributes Name Initial
real(kind=wp), public :: angstrom_h = 0.0_wp

Phase-1 Lagrangian minimum-thickness floor (m). Set from cfg%ocean%isopycnal%angstrom_h at setup, but only PASSED to the h-update kernels when the active vcoord is VCOORD_LAGRANGIAN (gated at the dyn call site). 0 ⇒ off ⇒ bit-identical. R7: floor lifts h but NOT hTr, so Tr=hTr/h shifts on a floored layer. Harmless for the adiabatic isopycnal config; thermo-on isopycnal correctness is OUT OF SCOPE for v1.
real(kind=wp), public :: cfl_max = 0.5_wp

Soft cap on per-face CFL before falling back to upwind.
logical, public :: conservative_floor = .false.

When .true. (set from cfg%ocean%isopycnal%conservative_floor at setup; requires angstrom_h > 0 + VCOORD_LAGRANGIAN) the injecting max(h, angstrom_h) floor is skipped in the h-update and replaced by the conservative per-column borrow in rdb_ocean_min_thickness, invoked from the dyn continuity site. Uses mt_h_new as scratch. Default .false. ⇒ legacy injecting floor ⇒ bit-identical.
logical, public :: hTr_holds_conc = .false.

True while the frozen tracer content hTr has been re-weighted onto the CURRENT h_layer so that the concentration T = hTr/h_layer is the (invariant) window-start value. Set by continuity_tracer_step_split in TR_MODE_ACCUMULATE, cleared by continuity_tracer_drain once it has converted hTr back onto the reconstructed window-start thickness hprev.

Read more…
type(scratch_3d_buffer_t), public :: h_face_left_x

Left-state thickness at east faces.
type(scratch_3d_buffer_t), public :: h_face_left_y

Left-state thickness at north faces.
type(scratch_3d_buffer_t), public :: h_face_right_x

Right-state thickness at east faces.
type(scratch_3d_buffer_t), public :: h_face_right_y

Right-state thickness at north faces.
real(kind=wp), public :: h_lim = 0.0_wp

Positive-definite thickness floor (m). Derived at setup: cfg%ocean%isopycnal%angstrom_h on VCOORD_LAGRANGIAN, else 0. Only consumed under positive_definite = .true.; at h_lim = 0 (every non-Lagrangian coord) the P1 max(edge, 2·h_lim) floor is inert (edges are already >= 0), but the P2 outflux limiter still engages to keep every layer >= 0 (avail = max(h, 0)).
real(kind=wp), public :: h_min = 1.0e-6_wp

Lower clip on cell-centred thickness during the update. Also used as the floor for ppm_limit_pos when use_ppm_limit_pos = .true. — same semantic role as MOM6’s GV%Angstrom_H for the vanishing-layer limiter.
real(kind=wp), public, allocatable :: h_win_start(:,:,:)

Layer thickness captured when an accumulation window OPENS. Paired with hTr_holds_conc: the per-stage concentration hold re-weights hTr onto the evolving h_layer, and the drain undoes it with THIS array, which returns hTr bit-for-bit to the frozen window-start content the drain has always consumed. (Undoing via the reconstructed hprev instead would be equivalent only where hprev == h_win_start exactly, and would silently convert any reconstruction residual — e.g. from a thin-layer h_min clip — into a tracer-mass drift.) Evolving (drained) layer thickness during the sub-cycle (m).
real(kind=wp), public, allocatable :: hprev_work(:,:,:)
logical, public :: is_init = .false.

True between init and destroy. Prefer this to allocated(...) — tracks GPU device attachment too.
type(scratch_3d_buffer_t), public :: mt_grounded

(nx,ny,1) grounded-column mask (1.0 = any layer below floor) for the borrow’s early-exit restructure: built in one coalesced pass, consumed by every borrow kernel in place of per-face column re-scans. Same lifetime/mapping as mt_h_new.
type(scratch_3d_buffer_t), public :: mt_h_new

Cell-centred (nx,ny,nz) scratch for the conservative minimum-thickness borrow (conservative_floor). Holds the floor-only target thickness field between the h-update and the h_layer overwrite. Unused (but allocated + mapped) when the knob is off — the DC kernels never touch it in that case.
integer, public :: n_limited_step = 0

P2 diagnostic counter: number of interior faces whose mass flux the positive-definite limiter scaled (θ_donor < 1) this outer-step call, summed over the zonal + meridional passes. Host-side scalar (the !$acc parallel loop reduction returns to the host); zeroed at the top of every continuity_tracer_step_split call. P3 drains it to the console stats line. 0 ⇒ no limiting fired (the healthy case; MOM6-style graceful degradation must be loud).
integer(kind=int64), public :: n_limited_total = 0_int64

P3 running total of n_limited_step across every continuity_tracer_step_split call (accumulated once per call, at the end — one add per split call, mirroring dyn%ntrunc_total). The driver drains its per-report DELTA to the console next to the CFL-truncation line. Grows ONLY when positive_definite = .true. (else the passes are skipped), so a nonzero total is itself the signal the limiter is active. Loud-by-design caveat: a uniform_z isopycnal (VCOORD_LAGRANGIAN) stack shows a permanently LARGE, steadily-growing count because ~most layers sit AT the floor and θ=0 correctly freezes their (massless) outflow every call — that is not a pathology. A healthy zstar/sigma run sits at 0. int64 (not int32 like dyn%ntrunc_total) precisely because that floored-stack case can accumulate > 2·10⁹ over a long run.
real(kind=wp), public, allocatable :: pa6(:,:,:)

CW parabola curvature a6 = 6·Tr − 3·(aL+aR) (rebuilt per pass).
real(kind=wp), public, allocatable :: pal(:,:,:)

CW parabola left-edge value per cell (rebuilt per pass).
real(kind=wp), public, allocatable :: par(:,:,:)

CW parabola right-edge value per cell (rebuilt per pass).
type(scratch_3d_buffer_t), public :: pd_theta

(nx,ny,nz) cell-centred per-donor availability factor θ(i,j,k) for the P2 positive-definite outflux limiter. Built once per direction pass (over the FULL range incl. ghosts, so any cell that can donate to a swept face has a current θ), then each interior face is scaled by its upwind donor’s θ. Allocated + device-mapped UNCONDITIONALLY (like mt_h_new/mt_grounded) — the enter_data contract has no config visibility, and at nx·ny·nz reals the footprint matches mt_h_new already carried; unused (but present) when positive_definite = .false..
logical, public :: positive_definite = .false.

Positive-definite split continuity master switch. When .true. the split layer continuity (continuity_tracer_step_split via continuity_zonal_flux/continuity_meridional_flux) keeps every layer >= h_lim: P1 floors the PPM reconstruction edges at 2·h_lim (MOM6’s positive-definite reconstruction); P2 (later) scales down per-donor outfluxes so no mass is created. Set from cfg%ocean%continuity%positive_definite at setup. Host-side control scalar — read into a local before each DC loop (never dereferenced inside a device loop). Default .false. ⇒ untaken branches only ⇒ bit-identical.
logical, public :: renorm_consistent_flux = .true.

&ocean_continuity_nml renorm_consistent_flux. When .true. the uhbt/vhbt renormalisation evaluates a layer whose upwind donor FLIPS under the correction as (u0 + du)·h_face(new donor), i.e. as the flux of its corrected velocity, and brackets the Newton solve with bisection (MOM6 zonal_flux_adjust). The historical model flux0 + du·h_face(new donor) keeps the OLD donor’s u0·h_old and is DISCONTINUOUS at the flip, by u0·(h_new − h_old)·w: whenever uhbt falls in that gap (a face where the corrected velocity must change sign across a thickness jump — a sigma layer over a bathymetric step, where h_old ≠ h_new after the PPM limiter flattens both edges), Newton has NO root, cycles for RENORM_MAXIT iterations and hands continuity a layer transport of the wrong SIGN. The layer η then departs from the barotropic η_end by O(η) at the step every such step — a spurious η dipole that pumps the (undamped) barotropic grid-scale mode. Default .true. (MOM6 behaviour); faces where no donor flips are bit-identical to the historical model, which .false. restores.
logical, public :: renorm_legacy_single_step = .false.

Use the pre-Newton single-linear-step uhbt renormalisation (donors picked at the UNCORRECTED velocity, no CFL bracket, no iteration) — wet/dry composition, set from cfg%ocean%wetdry%enable at setup. The Newton form’s donor RE-PICK + CFL bracket are what stabilise Lagrangian grounding (LAGRANGIAN_PGF_BUG.md §0), but under wet/dry they interact with the drying-front face gating and drive a drying column’s h_layer negative (test_ocean_wetdry_driver); wet/dry owns its positivity via ppm_limit_pos + the BT limiter and ran validated on the single-step form, so it keeps it (bit-identical there).
real(kind=wp), public :: t_dyn_rel_adv = 0.0_wp

Elapsed dynamics time (s) accumulated since the last tracer advect / accumulator reset. Adds dt once per outer step.
real(kind=wp), public, allocatable :: tr_flux_x(:,:,:)

Per-pass zonal tracer flux F (m^3 · concentration).
real(kind=wp), public, allocatable :: tr_flux_y(:,:,:)

Per-pass meridional tracer flux F (m^3 · concentration).
real(kind=wp), public, allocatable :: tr_work(:,:,:)

Current concentration Tr = hTr/hprev_work (rebuilt per pass).
integer, public :: tracer_recon = TRACER_RECON_PPM

Face-reconstruction scheme for the WINDOWED tracer-advection drain (Q6; dt_tracer_advect_ratio > 1 path only). 0 = CW-PPM (default, bit-identical), 1/2/3 = WENO5/7/9-Z swept-average (rung-adaptive, degrades near land/walls). Set at configure from &ocean_vmix_nml tracer_recon. Host-side control knob: the drain reads it to pick a face kernel — never dereferenced inside a device loop, so it rides the struct copyin. The every-step (ratio = 1) advect path is unaffected (stays CW-PPM) — WENO is a drain-only scheme in this wave.
real(kind=wp), public, allocatable :: uhh_x(:,:,:)

Per-pass limited zonal transport portion (m^3).
real(kind=wp), public, allocatable :: uhh_y(:,:,:)

Per-pass limited meridional transport portion (m^3).
real(kind=wp), public, allocatable :: uhr_x(:,:,:)

Remaining unspent zonal transport this window (m^3).
real(kind=wp), public, allocatable :: uhr_y(:,:,:)

Remaining unspent meridional transport this window (m^3).
real(kind=wp), public, allocatable :: uhtr(:,:,:)

Accumulated zonal face transport (m^3, area-weighted ·dt).
logical, public :: use_ppm_limit_pos = .false.

MOM6 PPM_limit_pos analogue. When .true., the PPM face-thickness reconstruction in continuity adds a positivity-preserving limiter that shrinks h_left / h_right toward h_centre whenever the parabolic fit would produce an interior minimum below h_min. At the h_centre ≤ h_min limit the reconstruction collapses to a constant (= upwind for that cell), bounding the mass flux by the actual layer thickness. Default off keeps the pre-knob behaviour bit-identical. Driver writes from cfg%ocean%continuity%ppm_limit_pos.
real(kind=wp), public, allocatable :: vhtr(:,:,:)

Accumulated meridional face transport (m^3, area-weighted ·dt).
logical, public :: vol_cfl = .false.

MOM6 vol_CFL analogue. When .false. (default) the continuity flux uses the PPM downwind EDGE value (the CFL→0 limit), bit-identical to the pre-knob behaviour. When .true. the donor-side face thickness is the swept-volume integral of the reconstructed parabola (volcfl_face), adding the missing O(CFL) term. Fixes the dt-sensitive near-bed residual at steep shelf breaks. Driver writes from cfg%ocean%continuity%vol_cfl.
logical, public :: windowed_advection = .false.

Gates the ALLOCATION of the Phase-2 windowed-advection state (13 3D arrays: uhtr/vhtr accumulators + the (6b) drain workspace, ~4.2 GB at 1000x800x50) — consumed only when dt_tracer_advect_ratio > 1. Latched from cfg BEFORE init(grid) by ocean_state_init_from_config (the same conditional-allocation contract as the default-off closures); the enter_data/exit_data/drain paths already guard on allocated().

Read more…

Type-Bound Procedures

procedure, public, non_overridable :: bytes => continuity_bytes
procedure, public, non_overridable :: destroy => continuity_destroy
procedure, public, non_overridable :: enter_data => continuity_enter_data
procedure, public, non_overridable :: exit_data => continuity_exit_data
procedure, public, non_overridable :: init => continuity_init

Functions

public pure elemental function ppm_mirror_h(h_nbr, h_loc, w_nbr) result(h_out)

Mirror-h at a land neighbour (spec §14 C2 / MOM6’s reflected-coast PPM): substitute the LOCAL cell’s thickness (or tracer) for a LAND neighbour’s held floor value so the PPM parabola sees a flat, reflected coast and the wet-side face value is not biased by the dry column. Branchless: h_out = w_nbr·h_nbr + (1-w_nbr)·h_loc. Wet neighbour (w_nbr=1) ⇒ h_out = h_nbr (literal no-op); land neighbour (w_nbr=0) ⇒ h_out = h_loc.

Read more…

Arguments

Type IntentOptional Attributes Name
real(kind=wp), intent(in) :: h_nbr

neighbour value
real(kind=wp), intent(in) :: h_loc

local-cell value (the mirror target)
real(kind=wp), intent(in) :: w_nbr

neighbour wet mask (0/1)

Return Value real(kind=wp)

public pure elemental function volcfl_face(h_edge, dh, curv3, cfl) result(h_face)

MOM6 swept-volume continuity-PPM face thickness (Lin & Rood / MOM_continuity_PPM flux_elem). Integrates the donor cell’s reconstructed parabola over the swept volume rather than sampling the edge value, adding the O(CFL) correction:

Read more…

Arguments

Type IntentOptional Attributes Name
real(kind=wp), intent(in) :: h_edge

downwind-facing donor PPM edge value
real(kind=wp), intent(in) :: dh

swept-oriented donor edge difference
real(kind=wp), intent(in) :: curv3

donor parabola curvature (h_L+h_R−2h)
real(kind=wp), intent(in) :: cfl

donor-cell Courant number (>= 0)

Return Value real(kind=wp)

private pure function continuity_bytes(this) result(nbytes)

Counted allocatable footprint of the continuity-PPM slot (0 when unallocated). One arr_bytes term per array — add a term here when a new allocatable joins the type.

Arguments

Type IntentOptional Attributes Name
class(continuity_t), intent(in) :: this

Return Value integer(kind=int64)

private pure function weno_face_conc_x(nx, ny, nz, tr, wet_T, cc, jj, kk, d, avail_up, avail_down, rung_max, cfl) result(conc)

Swept-average WENO donor concentration at a zonal face. cc is the donor cell (i-index), d = +1 (u>0, downwind toward +i) or d = -1 (u<0, downwind toward -i). Gathers the mirrored/clamped stencil in downwind-positive order and dispatches to the coastal swept-average face helper at the highest feasible rung.

Arguments

Type IntentOptional Attributes Name
integer, intent(in) :: nx
integer, intent(in) :: ny
integer, intent(in) :: nz
real(kind=wp), intent(in) :: tr(nx,ny,nz)
real(kind=wp), intent(in) :: wet_T(nx,ny)
integer, intent(in) :: cc
integer, intent(in) :: jj
integer, intent(in) :: kk
integer, intent(in) :: d
integer, intent(in) :: avail_up
integer, intent(in) :: avail_down
integer, intent(in) :: rung_max
real(kind=wp), intent(in) :: cfl

Return Value real(kind=wp)

private pure function weno_face_conc_y(nx, ny, nz, tr, wet_T, ii, cc, kk, d, avail_up, avail_down, rung_max, cfl) result(conc)

Meridional analogue of weno_face_conc_x. cc is the donor cell (j-index); the stencil steps along j with d = +1 (u>0) / -1 (u<0).

Arguments

Type IntentOptional Attributes Name
integer, intent(in) :: nx
integer, intent(in) :: ny
integer, intent(in) :: nz
real(kind=wp), intent(in) :: tr(nx,ny,nz)
real(kind=wp), intent(in) :: wet_T(nx,ny)
integer, intent(in) :: ii
integer, intent(in) :: cc
integer, intent(in) :: kk
integer, intent(in) :: d
integer, intent(in) :: avail_up
integer, intent(in) :: avail_down
integer, intent(in) :: rung_max
real(kind=wp), intent(in) :: cfl

Return Value real(kind=wp)


Subroutines

public pure subroutine continuity_apply_fluxes_barotropic(bs, dt)

Forward-Euler step: h <- h - dt * flux_h. The full split-explicit RK2 scheme (Phase 4) wraps two of these calls around an RK2 averaging pass; for Phase 2 this single-stage step is enough to exercise the kernel under the lake-at-rest, Gaussian-hump, and mass-conservation tests.

Arguments

Type IntentOptional Attributes Name
type(barotropic_state_t), intent(inout) :: bs
real(kind=wp), intent(in) :: dt

public pure subroutine continuity_compute_fluxes_barotropic(grid, metrics, this, bs)

PPM face reconstruction + per-face mass flux + cell-centred flux divergence for the barotropic C-grid state.

Read more…

Arguments

Type IntentOptional Attributes Name
type(hgrid_t), intent(in) :: grid
type(ocean_metrics_t), intent(in) :: metrics
type(continuity_t), intent(inout) :: this
type(barotropic_state_t), intent(inout) :: bs

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

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

Read more…

Arguments

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

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

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

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

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

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

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

public subroutine continuity_tracer_drain(grid, metrics, this, ms, ratio, bc)

Phase-2 (6b) windowed horizontal tracer-advection drain.

Read more…

Arguments

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

dt_tracer_advect_ratio for this window; sets the fixed-budget pass count max_iter = 2·ratio+1 (ratio is an integer ⇒ ceil(ratio) = ratio).
type(ocean_bc_state_t), intent(in), optional :: bc

public subroutine continuity_tracer_step_split(grid, metrics, this, ms, dt, uhbt, vhbt, bc, mle, mle_fold_active, tracer_mode, h_min, visc_rem_u, visc_rem_v, u_cor, v_cor)

Production entry point for the directionally-split continuity + tracer step. Interleaves the two so the CWC discrete theorem holds in the split form:

Read more…

Arguments

Type IntentOptional Attributes Name
type(hgrid_t), intent(in) :: grid
type(ocean_metrics_t), intent(in) :: metrics
type(continuity_t), intent(inout) :: this
type(multilayer_state_t), intent(inout) :: ms
real(kind=wp), intent(in) :: dt
real(kind=wp), intent(in), optional :: uhbt(:,:)
real(kind=wp), intent(in), optional :: vhbt(:,:)
type(ocean_bc_state_t), intent(in), optional :: bc

When present, per-edge OBC tags gate the wall-zero step inside the flux kernels. OBC_WALL keeps the Phase 3 closure; OBC_OPEN (and other non-wall tags) leaves the computed mass flux at the wall face for the downstream transport. Absent ⇒ closed-wall everywhere.
type(ocean_mle_t), intent(in), optional :: mle

Fox-Kemper mixed-layer-eddy transports (B5). When present and enabled, mle%uhml/vhml are folded into the per-layer mass fluxes AFTER each direction’s flux fill and BEFORE the matching tracer advect + divergence — so the augmented flux transports both h and tracers (conservative; velocity untouched). Absent / disabled ⇒ bit-identical no-op.
logical, intent(in), optional :: mle_fold_active

Gates the Fox-Kemper fold to the THERMO cadence. Absent or .true. ⇒ the fold applies (bit-identical default — the case at dt_therm_ratio = 1, where every step is a thermo step). .false. skips the fold so the stale FK transports (computed once per thermo interval) are NOT re-applied on the intervening non-thermo outer steps when dt_therm_ratio > 1.
integer, intent(in), optional :: tracer_mode

Phase 2 (6b) windowed-advection mode. TR_MODE_ADVECT (default, absent) ⇒ the historical fused path: advance h AND advect tracers each call (bit-identical to pre-6b). TR_MODE_ACCUMULATE ⇒ advance h, accumulate 0.5·mass_flux·dt into this%uhtr/vhtr (one += per RK2 stage, weight 0.5 baked in — closes the reconstruction against the RK2-averaged h), and SKIP the per-step tracer advect so hTr stays frozen until the boundary drain.
real(kind=wp), intent(in), optional :: h_min

Phase-1 Lagrangian minimum-thickness floor (m). When > 0, passed to continuity_apply_zonal/_meridional to clamp h_new >= h_min. Absent or 0 ⇒ off ⇒ bit-identical.
real(kind=wp), intent(in), optional :: visc_rem_u(:,:,:)

Per-layer viscous remnant gamma_k on east / north faces. Forwarded to the flux renormalisers, where it weights the barotropic increment (MOM6 u_cor = u + du*visc_rem). Absent => gamma == 1, bit-identical.
real(kind=wp), intent(in), optional :: visc_rem_v(:,:,:)

Per-layer viscous remnant gamma_k on east / north faces. Forwarded to the flux renormalisers, where it weights the barotropic increment (MOM6 u_cor = u + du*visc_rem). Absent => gamma == 1, bit-identical.
real(kind=wp), intent(inout), optional :: u_cor(:,:,:)

MOM6 u_cor/v_cor destinations — the step TIME-MEAN velocity (u_av/v_av), never the prognostic. Absent => flux-only.
real(kind=wp), intent(inout), optional :: v_cor(:,:,:)

MOM6 u_cor/v_cor destinations — the step TIME-MEAN velocity (u_av/v_av), never the prognostic. Absent => flux-only.

public pure subroutine ppm_cell_limiter(h_centre, h_left, h_right)

Colella-Woodward 1984 eq 1.10 monotonic limiter on the parabolic profile in a single cell. Three branches:

Read more…

Arguments

Type IntentOptional Attributes Name
real(kind=wp), intent(in) :: h_centre
real(kind=wp), intent(inout) :: h_left
real(kind=wp), intent(inout) :: h_right

public pure subroutine ppm_limit_pos(h_centre, h_left, h_right, h_min)

Positivity-preserving limiter on the PPM reconstruction. Mirrors MOM6’s PPM_limit_pos: when the parabolic fit predicts a minimum interior to the cell that dips below h_min, shrink h_left / h_right toward h_centre so the minimum sits at exactly h_min. Pure scalar form per cell; runs after ppm_cell_limiter so the monotonic-limited reconstruction is the input.

Read more…

Arguments

Type IntentOptional Attributes Name
real(kind=wp), intent(in) :: h_centre
real(kind=wp), intent(inout) :: h_left
real(kind=wp), intent(inout) :: h_right
real(kind=wp), intent(in) :: h_min

public pure subroutine ppm_limited_slope(h_im1, h_i, h_ip1, dh)

Van Leer monotonized centred slope for cell i. Returns 0 at local extrema (sign change between left and right differences) and the slope-limited centred derivative otherwise. Standard PPM convention; see Colella-Woodward 1984.

Read more…

Arguments

Type IntentOptional Attributes Name
real(kind=wp), intent(in) :: h_im1
real(kind=wp), intent(in) :: h_i
real(kind=wp), intent(in) :: h_ip1
real(kind=wp), intent(out) :: dh

private pure subroutine accumulate_flux_x(nx_face, ny, nz, dt, mass_flux_x, uhtr)

Accumulate one RK2 stage’s zonal mass flux into the windowed accumulator with the ½ RK2 weight baked in: uhtr += 0.5·mass_flux_x·dt. mass_flux_x is already area-weighted (m³/s = u·h_face·dy_cu); the product is m³.

Arguments

Type IntentOptional Attributes Name
integer, intent(in) :: nx_face
integer, intent(in) :: ny
integer, intent(in) :: nz
real(kind=wp), intent(in) :: dt
real(kind=wp), intent(in) :: mass_flux_x(nx_face,ny,nz)
real(kind=wp), intent(inout) :: uhtr(nx_face,ny,nz)

private pure subroutine accumulate_flux_y(nx, ny_face, nz, dt, mass_flux_y, vhtr)

Meridional analogue of accumulate_flux_x.

Arguments

Type IntentOptional Attributes Name
integer, intent(in) :: nx
integer, intent(in) :: ny_face
integer, intent(in) :: nz
real(kind=wp), intent(in) :: dt
real(kind=wp), intent(in) :: mass_flux_y(nx,ny_face,nz)
real(kind=wp), intent(inout) :: vhtr(nx,ny_face,nz)

private pure subroutine continuity_apply_fluxes(ms, dt, h_min)

Test-only (no production caller): unsplit apply, paired with continuity_compute_fluxes as the split path’s reference oracle. Per-layer forward-Euler thickness update.

Read more…

Arguments

Type IntentOptional Attributes Name
type(multilayer_state_t), intent(inout) :: ms
real(kind=wp), intent(in) :: dt
real(kind=wp), intent(in), optional :: h_min

Minimum-thickness floor (m). 0 or absent ⇒ off ⇒ bit-identical.

private pure subroutine continuity_apply_meridional(grid, metrics, ms, dt, h_min)

Apply the meridional (y-flux) thickness update on top of the zonally-updated state: h(i,j,k) ← h(i,j,k) - dt · (Φy(i,j+1,k) - Φy(i,j,k)) · iareaT Adds the y-divergence to flux_h_layer so the field ends the split step holding the total horizontal divergence that the vertical-advection kernel consumes (w_interface(k+1) = w(k) - flux_h_layer(k)). Φy carries dx_cv; iareaT = inv_dy on uniform metrics.

Read more…

Arguments

Type IntentOptional Attributes Name
type(hgrid_t), intent(in) :: grid
type(ocean_metrics_t), intent(in) :: metrics
type(multilayer_state_t), intent(inout) :: ms
real(kind=wp), intent(in) :: dt
real(kind=wp), intent(in), optional :: h_min

Minimum-thickness floor (m). 0 or absent ⇒ off ⇒ bit-identical.

private pure subroutine continuity_apply_zonal(grid, metrics, ms, dt, h_min)

Apply the zonal (x-flux) thickness update: h(i,j,k) ← h(i,j,k) - dt · (Φx(i+1,j,k) - Φx(i,j,k)) · iareaT Overwrites flux_h_layer with the x-divergence so the companion meridional apply can accumulate the total. Φx is the width-weighted transport (m³/s); iareaT closes the divergence to a per-area rate (= inv_dx on uniform metrics).

Read more…

Arguments

Type IntentOptional Attributes Name
type(hgrid_t), intent(in) :: grid
type(ocean_metrics_t), intent(in) :: metrics
type(multilayer_state_t), intent(inout) :: ms
real(kind=wp), intent(in) :: dt
real(kind=wp), intent(in), optional :: h_min

Minimum-thickness floor (m). 0 or absent ⇒ off ⇒ bit-identical.

private pure subroutine continuity_compute_fluxes(grid, metrics, this, ms)

Test-only (no production caller): the unsplit reference path, kept as the oracle the split production path is checked against. Multilayer counterpart to continuity_compute_fluxes_barotropic: identical PPM reconstruction + upwind face pick + flux divergence, lifted per-layer. Each k-slice is independent (the PPM stencil reads only the same k), so the do-concurrent kernels parallelize over (k, j, i) simultaneously for GPU occupancy.

Read more…

Arguments

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

private subroutine continuity_destroy(this)

Arguments

Type IntentOptional Attributes Name
class(continuity_t), intent(inout) :: this

private subroutine continuity_enter_data(this)

Bare copyin(this) removed (stack-descriptor map → AMD cross-slot overlap; see ocean_surfstress_enter_data). The face buffers attach below; ct-descriptor presence (so DCs touching ct%h_face_left_x%data don’t per-launch memcpy) comes from the root copyin(state) in ocean_state_enter_data. A V100 A/B with copyin(this) gone is bit-identical and faster overall, so the root copy fully covers it.

Arguments

Type IntentOptional Attributes Name
class(continuity_t), intent(inout) :: this

private subroutine continuity_enter_data_impl(this)

Arguments

Type IntentOptional Attributes Name
type(continuity_t), intent(inout) :: this

private subroutine continuity_exit_data(this)

Arguments

Type IntentOptional Attributes Name
class(continuity_t), intent(inout) :: this

private subroutine continuity_exit_data_impl(this)

Arguments

Type IntentOptional Attributes Name
type(continuity_t), intent(inout) :: this

private subroutine continuity_init(this, grid, nz_ml)

Allocate the 4 face-reconstruction scratch buffers sized at (nx_face, ny_face, nz). Default nz=1 covers the barotropic kernel; passing nz_ml sizes them for the multilayer kernel without forcing a separate init routine. Ocean init passes state%multilayer%nz_ml when the multilayer state is in play.

Arguments

Type IntentOptional Attributes Name
class(continuity_t), intent(inout) :: this
type(hgrid_t), intent(in) :: grid
integer, intent(in), optional :: nz_ml

private pure subroutine continuity_meridional_flux(grid, metrics, this, ms, dt, vhbt, bc, visc_rem, v_cor)

Meridional (y-only) PPM reconstruction + per-face mass flux. Mirror of continuity_zonal_flux, with the same optional vhbt transport-constraint renormalisation. Writes ms%mass_flux_y_layer. Walls at j=1 and j=ny+1 zeroed. In the Lie split this runs after the zonal apply, so it reconstructs against the already-updated ms%h_layer.

Arguments

Type IntentOptional Attributes Name
type(hgrid_t), intent(in) :: grid
type(ocean_metrics_t), intent(in) :: metrics
type(continuity_t), intent(inout) :: this
type(multilayer_state_t), intent(inout) :: ms
real(kind=wp), intent(in) :: dt

Time increment (s). Used only for the swept-volume CFL when this%vol_cfl = .true.; ignored (bit-identical) otherwise.
real(kind=wp), intent(in), optional :: vhbt(:,:)

Time-mean north-face transport from barotropic substep (m³/s), shape (nx, ny+1). Width-weighted (carries dx_cv).
type(ocean_bc_state_t), intent(in), optional :: bc
real(kind=wp), intent(in), optional :: visc_rem(:,:,:)

Per-layer viscous remnant gamma_k. Absent => 1, bit-identical.
real(kind=wp), intent(inout), optional :: v_cor(:,:,:)

Forwarded MOM6 v_cor destination — separate time-mean field.

private subroutine continuity_step_split(grid, metrics, this, ms, dt)

Test-only (no production caller): continuity-only split wrapper; production runs the tracer-interleaved continuity_tracer_step_split. Directionally-split (Lie) PPM continuity step over dt:

Read more…

Arguments

Type IntentOptional Attributes Name
type(hgrid_t), intent(in) :: grid
type(ocean_metrics_t), intent(in) :: metrics
type(continuity_t), intent(inout) :: this
type(multilayer_state_t), intent(inout) :: ms
real(kind=wp), intent(in) :: dt

private pure subroutine continuity_zonal_flux(grid, metrics, this, ms, dt, uhbt, bc, visc_rem, u_cor)

Zonal (x-only) PPM reconstruction + per-face mass flux on the multilayer C-grid. Companion to continuity_meridional_flux for the directionally-split (Lie) continuity step. Writes ms%mass_flux_x_layer and leaves mass_flux_y_layer / flux_h_layer untouched. Wall faces at i=1 and i=nx+1 are zeroed (closed-wall BC).

Read more…

Arguments

Type IntentOptional Attributes Name
type(hgrid_t), intent(in) :: grid
type(ocean_metrics_t), intent(in) :: metrics
type(continuity_t), intent(inout) :: this
type(multilayer_state_t), intent(inout) :: ms
real(kind=wp), intent(in) :: dt

Time increment (s). Used only for the swept-volume CFL when this%vol_cfl = .true.; ignored (bit-identical) otherwise.
real(kind=wp), intent(in), optional :: uhbt(:,:)

Time-mean east-face transport from barotropic substep (m³/s), shape (nx+1, ny). Width-weighted (carries dy_cu) so it constrains the same transport mass_flux_x_layer now holds.
type(ocean_bc_state_t), intent(in), optional :: bc

Per-edge OBC tags. Default (absent) -> OBC_WALL on both ends. Non-WALL tags skip the wall-zero step at that edge.
real(kind=wp), intent(in), optional :: visc_rem(:,:,:)

Per-layer viscous remnant γ_k, forwarded to renormalise_zonal_flux_to_uhbt. Absent ⇒ γ ≡ 1, bit-identical.
real(kind=wp), intent(inout), optional :: u_cor(:,:,:)

Forwarded MOM6 u_cor destination — a SEPARATE time-mean field, never the prognostic. Absent ⇒ flux-only, bit-identical.

private pure subroutine drain_avail_limit(nx, ny, nz, areaT, h_min, h_end, uhtr, vhtr, scratch)

Conservative upfront availability limiter on the accumulated window transports uhtr/vhtr, applied BEFORE drain_reconstruct_hprev.

Read more…

Arguments

Type IntentOptional Attributes Name
integer, intent(in) :: nx
integer, intent(in) :: ny
integer, intent(in) :: nz
real(kind=wp), intent(in) :: areaT(nx,ny)
real(kind=wp), intent(in) :: h_min
real(kind=wp), intent(in) :: h_end(nx,ny,nz)
real(kind=wp), intent(inout) :: uhtr(nx+1,ny,nz)
real(kind=wp), intent(inout) :: vhtr(nx,ny+1,nz)
real(kind=wp), intent(inout) :: scratch(nx,ny,nz)

Per-cell inflow scale factor in [0,1] (centre-shaped slot).

private pure subroutine drain_avail_scale_x(nx, ny, nz, scratch, uhtr)

Scale each interior x-face transport by its INFLOW-receiving cell’s factor (drain_avail_limit step 2). Face i between cell (i-1) and cell (i): uhtr(i)>0 ⇒ receiver i, uhtr(i)<0 ⇒ receiver i-1.

Arguments

Type IntentOptional Attributes Name
integer, intent(in) :: nx
integer, intent(in) :: ny
integer, intent(in) :: nz
real(kind=wp), intent(in) :: scratch(nx,ny,nz)
real(kind=wp), intent(inout) :: uhtr(nx+1,ny,nz)

private pure subroutine drain_avail_scale_y(nx, ny, nz, scratch, vhtr)

Meridional analogue of drain_avail_scale_x. Face j between cell (j-1) and cell (j): vhtr(j)>0 ⇒ receiver j, vhtr(j)<0 ⇒ receiver j-1.

Arguments

Type IntentOptional Attributes Name
integer, intent(in) :: nx
integer, intent(in) :: ny
integer, intent(in) :: nz
real(kind=wp), intent(in) :: scratch(nx,ny,nz)
real(kind=wp), intent(inout) :: vhtr(nx,ny+1,nz)

private pure subroutine drain_copy_3d(n1, n2, n3, src, dst)

dst = src (explicit-shape device copy).

Arguments

Type IntentOptional Attributes Name
integer, intent(in) :: n1
integer, intent(in) :: n2
integer, intent(in) :: n3
real(kind=wp), intent(in) :: src(n1,n2,n3)
real(kind=wp), intent(inout) :: dst(n1,n2,n3)

private pure subroutine drain_fill_conc(nx, ny, nz, hTr, hprev, tr)

Concentration field Tr = hTr / max(hprev, DRAIN_MIN_H) for the WENO drain (the CW path builds the same field as the first loop of drain_parabola_*; factored out so the WENO path can reuse it without the parabola coefficients).

Arguments

Type IntentOptional Attributes Name
integer, intent(in) :: nx
integer, intent(in) :: ny
integer, intent(in) :: nz
real(kind=wp), intent(in) :: hTr(nx,ny,nz)
real(kind=wp), intent(in) :: hprev(nx,ny,nz)
real(kind=wp), intent(inout) :: tr(nx,ny,nz)

private pure subroutine drain_limit_x(nx, ny, nz, areaT, h_min, uhr_x, hprev, uhh_x)

MOM6 hup/hlos/min_h two-test limiter on the zonal face transport (volume units). Face i between cell (i-1) and cell (i). Positive flow (uhr_x(i) > 0), donor = cell (i-1): hup = areaT(i-1)·hprev(i-1) − areaT(i-1)·min_h hlos = max(0, −uhr_x(i-1)) (already-committed outflow via the donor’s OTHER (west) face) cap when (hup−hlos)−uhr < 0 AND 0.5·hup−uhr < 0. Negative flow mirror, donor = cell (i).

Arguments

Type IntentOptional Attributes Name
integer, intent(in) :: nx
integer, intent(in) :: ny
integer, intent(in) :: nz
real(kind=wp), intent(in) :: areaT(nx,ny)
real(kind=wp), intent(in) :: h_min
real(kind=wp), intent(in) :: uhr_x(nx+1,ny,nz)
real(kind=wp), intent(in) :: hprev(nx,ny,nz)
real(kind=wp), intent(inout) :: uhh_x(nx+1,ny,nz)

private pure subroutine drain_limit_y(nx, ny, nz, areaT, h_min, uhr_y, hprev, uhh_y)

Meridional analogue of drain_limit_x. Face j between cell (i,j-1) and cell (i,j); positive donor = cell (i,j-1).

Arguments

Type IntentOptional Attributes Name
integer, intent(in) :: nx
integer, intent(in) :: ny
integer, intent(in) :: nz
real(kind=wp), intent(in) :: areaT(nx,ny)
real(kind=wp), intent(in) :: h_min
real(kind=wp), intent(in) :: uhr_y(nx,ny+1,nz)
real(kind=wp), intent(in) :: hprev(nx,ny,nz)
real(kind=wp), intent(inout) :: uhh_y(nx,ny+1,nz)

private pure subroutine drain_parabola_x(nx, ny, nz, wet_T, hTr, hprev, tr, aL, aR, a6)

Rebuild the per-cell zonal CW PPM parabola from the CURRENT Tr = hTr/hprev (V2 — per pass). Interior cells (3..nx-2) use the limited PPM edges; the 2-cell boundary band falls back to PCM (aL=aR=Tr ⇒ swept reduces to the donor value, 1st order), matching tracer_advect_zonal_one_impl’s near-wall band. TODO(MOM6-fidelity): MOM6 advect_tracer keeps full PPM up to the wall (dropping to PCM only at genuine local extrema / zero mask2dCu faces), so we are 1st-order in the 2 cells nearest a true WALL where MOM6 is PPM-with-mask (periodic seams are fine — the wrap restores full PPM). Tracked divergence; revisit if near-wall tracer diffusion matters. a6 = 6·Tr − 3·(aL+aR). Mirror-T at land neighbours (C2); bit-identical for all-wet (wet_T≡1).

Arguments

Type IntentOptional Attributes Name
integer, intent(in) :: nx
integer, intent(in) :: ny
integer, intent(in) :: nz
real(kind=wp), intent(in) :: wet_T(nx,ny)
real(kind=wp), intent(in) :: hTr(nx,ny,nz)
real(kind=wp), intent(in) :: hprev(nx,ny,nz)
real(kind=wp), intent(inout) :: tr(nx,ny,nz)
real(kind=wp), intent(inout) :: aL(nx,ny,nz)
real(kind=wp), intent(inout) :: aR(nx,ny,nz)
real(kind=wp), intent(inout) :: a6(nx,ny,nz)

private pure subroutine drain_parabola_y(nx, ny, nz, wet_T, hTr, hprev, tr, aL, aR, a6)

Meridional analogue of drain_parabola_x. aL = south-edge, aR = north-edge value of each cell. Mirror-T at land neighbours (C2); bit-identical for all-wet.

Arguments

Type IntentOptional Attributes Name
integer, intent(in) :: nx
integer, intent(in) :: ny
integer, intent(in) :: nz
real(kind=wp), intent(in) :: wet_T(nx,ny)
real(kind=wp), intent(in) :: hTr(nx,ny,nz)
real(kind=wp), intent(in) :: hprev(nx,ny,nz)
real(kind=wp), intent(inout) :: tr(nx,ny,nz)
real(kind=wp), intent(inout) :: aL(nx,ny,nz)
real(kind=wp), intent(inout) :: aR(nx,ny,nz)
real(kind=wp), intent(inout) :: a6(nx,ny,nz)

private pure subroutine drain_reconstruct_hprev(nx, ny, nz, areaT, iareaT, h_end, uhtr, vhtr, hprev)

hprev = max(0, areaT·h_end + div(uhtr,vhtr)) · iareaT, then the vanishing-layer hatch hprev += max(0, 1e-13·hprev − h_end) (Adcroft & Hallberg 2006; reuse of VANISHING_LAYER_TOL thinking).

Arguments

Type IntentOptional Attributes Name
integer, intent(in) :: nx
integer, intent(in) :: ny
integer, intent(in) :: nz
real(kind=wp), intent(in) :: areaT(nx,ny)
real(kind=wp), intent(in) :: iareaT(nx,ny)
real(kind=wp), intent(in) :: h_end(nx,ny,nz)
real(kind=wp), intent(in) :: uhtr(nx+1,ny,nz)
real(kind=wp), intent(in) :: vhtr(nx,ny+1,nz)
real(kind=wp), intent(inout) :: hprev(nx,ny,nz)

private pure subroutine drain_rescale_hTr(nx, ny, nz, h_new, h_old, hTr)

Re-weight a tracer’s thickness-weighted content onto a new layer thickness, holding the CONCENTRATION fixed:

Read more…

Arguments

Type IntentOptional Attributes Name
integer, intent(in) :: nx
integer, intent(in) :: ny
integer, intent(in) :: nz
real(kind=wp), intent(in) :: h_new(nx,ny,nz)
real(kind=wp), intent(in) :: h_old(nx,ny,nz)
real(kind=wp), intent(inout) :: hTr(nx,ny,nz)

private pure subroutine drain_rescale_hTr_budget(nx, ny, nz, h_new, h_old, w, hTr, budget_adv)

drain_rescale_hTr + the closed-budget fill. The concentration hold / un-hold pair is NOT content-conserving cell by cell – hTr := Tr·h moves Σ areaT·hTr by Σ areaT·Tr·δh, which is only zero when Tr is uniform – so both halves have to be recorded or the closed budget is only valid on window boundaries and wobbles at every mid-window report. Recording both makes them cancel exactly, since the un-hold is the arithmetic inverse of the accumulated hold.

Read more…

Arguments

Type IntentOptional Attributes Name
integer, intent(in) :: nx
integer, intent(in) :: ny
integer, intent(in) :: nz
real(kind=wp), intent(in) :: h_new(nx,ny,nz)
real(kind=wp), intent(in) :: h_old(nx,ny,nz)
real(kind=wp), intent(in) :: w
real(kind=wp), intent(inout) :: hTr(nx,ny,nz)
real(kind=wp), intent(inout) :: budget_adv(nx,ny,nz)

private pure subroutine drain_subtract_3d(n1, n2, n3, sub, fld)

fld = fld - sub (uhr -= uhh).

Arguments

Type IntentOptional Attributes Name
integer, intent(in) :: n1
integer, intent(in) :: n2
integer, intent(in) :: n3
real(kind=wp), intent(in) :: sub(n1,n2,n3)
real(kind=wp), intent(inout) :: fld(n1,n2,n3)

private pure subroutine drain_swept_flux_x(nx, ny, nz, areaT, uhh, hprev, aL, aR, a6, F)

MOM6 swept-average CW parabola flux for the zonal faces. Face i, donor = cell (i-1) if uhh>0 else cell (i). Per-pass Courant CFL = |uhh| / (areaT·hprev) on the donor, clamped [0,1]. uhh >= 0: F = uhh·( aR − 0.5·CFL·((aR−aL) − a6·(1 − ⅔·CFL)) ) uhh < 0: F = uhh·( aL + 0.5·CFL·((aR−aL) + a6·(1 − ⅔·CFL)) )

Arguments

Type IntentOptional Attributes Name
integer, intent(in) :: nx
integer, intent(in) :: ny
integer, intent(in) :: nz
real(kind=wp), intent(in) :: areaT(nx,ny)
real(kind=wp), intent(in) :: uhh(nx+1,ny,nz)
real(kind=wp), intent(in) :: hprev(nx,ny,nz)
real(kind=wp), intent(in) :: aL(nx,ny,nz)
real(kind=wp), intent(in) :: aR(nx,ny,nz)
real(kind=wp), intent(in) :: a6(nx,ny,nz)
real(kind=wp), intent(inout) :: F(nx+1,ny,nz)

private pure subroutine drain_swept_flux_x_weno(nx, ny, nz, nghost, periodic, areaT, uhh, hprev, tr, wet_T, recon, F)

WENO analogue of drain_swept_flux_x: swept-average zonal face flux from the concentration field tr using the WENO rung ladder. recon is TRACER_RECON_WENO5/7/9 (1/2/3); the internal rung_max is recon + 1 (WENO5→2, WENO7→3, WENO9→4).

Read more…

Arguments

Type IntentOptional Attributes Name
integer, intent(in) :: nx
integer, intent(in) :: ny
integer, intent(in) :: nz
integer, intent(in) :: nghost
logical, intent(in) :: periodic
real(kind=wp), intent(in) :: areaT(nx,ny)
real(kind=wp), intent(in) :: uhh(nx+1,ny,nz)
real(kind=wp), intent(in) :: hprev(nx,ny,nz)
real(kind=wp), intent(in) :: tr(nx,ny,nz)
real(kind=wp), intent(in) :: wet_T(nx,ny)
integer, intent(in) :: recon
real(kind=wp), intent(inout) :: F(nx+1,ny,nz)

private pure subroutine drain_swept_flux_y(nx, ny, nz, areaT, uhh, hprev, aL, aR, a6, F)

Meridional analogue of drain_swept_flux_x.

Arguments

Type IntentOptional Attributes Name
integer, intent(in) :: nx
integer, intent(in) :: ny
integer, intent(in) :: nz
real(kind=wp), intent(in) :: areaT(nx,ny)
real(kind=wp), intent(in) :: uhh(nx,ny+1,nz)
real(kind=wp), intent(in) :: hprev(nx,ny,nz)
real(kind=wp), intent(in) :: aL(nx,ny,nz)
real(kind=wp), intent(in) :: aR(nx,ny,nz)
real(kind=wp), intent(in) :: a6(nx,ny,nz)
real(kind=wp), intent(inout) :: F(nx,ny+1,nz)

private pure subroutine drain_swept_flux_y_weno(nx, ny, nz, nghost, periodic, areaT, uhh, hprev, tr, wet_T, recon, F)

Meridional analogue of drain_swept_flux_x_weno.

Arguments

Type IntentOptional Attributes Name
integer, intent(in) :: nx
integer, intent(in) :: ny
integer, intent(in) :: nz
integer, intent(in) :: nghost
logical, intent(in) :: periodic
real(kind=wp), intent(in) :: areaT(nx,ny)
real(kind=wp), intent(in) :: uhh(nx,ny+1,nz)
real(kind=wp), intent(in) :: hprev(nx,ny,nz)
real(kind=wp), intent(in) :: tr(nx,ny,nz)
real(kind=wp), intent(in) :: wet_T(nx,ny)
integer, intent(in) :: recon
real(kind=wp), intent(inout) :: F(nx,ny+1,nz)

private pure subroutine drain_update_h_x(nx, ny, nz, iareaT, uhh, hprev)

hprev(i) -= (uhh(i+1) - uhh(i))·iareaT (volume div → thickness).

Arguments

Type IntentOptional Attributes Name
integer, intent(in) :: nx
integer, intent(in) :: ny
integer, intent(in) :: nz
real(kind=wp), intent(in) :: iareaT(nx,ny)
real(kind=wp), intent(in) :: uhh(nx+1,ny,nz)
real(kind=wp), intent(inout) :: hprev(nx,ny,nz)

private pure subroutine drain_update_h_y(nx, ny, nz, iareaT, uhh, hprev)

Meridional analogue of drain_update_h_x.

Arguments

Type IntentOptional Attributes Name
integer, intent(in) :: nx
integer, intent(in) :: ny
integer, intent(in) :: nz
real(kind=wp), intent(in) :: iareaT(nx,ny)
real(kind=wp), intent(in) :: uhh(nx,ny+1,nz)
real(kind=wp), intent(inout) :: hprev(nx,ny,nz)

private pure subroutine drain_update_tracer_x(nx, ny, nz, iareaT, F, hTr)

hTr(i) -= (F(i+1) - F(i))·iareaT (zonal flux divergence; dt is already baked into uhh⊂uhtr, so no dt here).

Arguments

Type IntentOptional Attributes Name
integer, intent(in) :: nx
integer, intent(in) :: ny
integer, intent(in) :: nz
real(kind=wp), intent(in) :: iareaT(nx,ny)
real(kind=wp), intent(in) :: F(nx+1,ny,nz)
real(kind=wp), intent(inout) :: hTr(nx,ny,nz)

private pure subroutine drain_update_tracer_x_budget(nx, ny, nz, iareaT, F, w, hTr, budget_adv)

drain_update_tracer_x + the closed-budget fill: the SAME increment written to hTr is accumulated (times w) into budget_adv, so the console out term sees the horizontal tracer transport the windowed drain performs. Without this the drain moves tracer that ms%*_budget_horiz_adv never records, and the Heat/Salt Error column has to fall back to raw drift (which then reports a live surface flux as a “leak”). w is DRAIN_BUDGET_POST_AVERAGE_WEIGHT for every drain call site.

Read more…

Arguments

Type IntentOptional Attributes Name
integer, intent(in) :: nx
integer, intent(in) :: ny
integer, intent(in) :: nz
real(kind=wp), intent(in) :: iareaT(nx,ny)
real(kind=wp), intent(in) :: F(nx+1,ny,nz)
real(kind=wp), intent(in) :: w
real(kind=wp), intent(inout) :: hTr(nx,ny,nz)
real(kind=wp), intent(inout) :: budget_adv(nx,ny,nz)

private pure subroutine drain_update_tracer_y(nx, ny, nz, iareaT, F, hTr)

Meridional analogue of drain_update_tracer_x.

Arguments

Type IntentOptional Attributes Name
integer, intent(in) :: nx
integer, intent(in) :: ny
integer, intent(in) :: nz
real(kind=wp), intent(in) :: iareaT(nx,ny)
real(kind=wp), intent(in) :: F(nx,ny+1,nz)
real(kind=wp), intent(inout) :: hTr(nx,ny,nz)

private pure subroutine drain_update_tracer_y_budget(nx, ny, nz, iareaT, F, w, hTr, budget_adv)

Meridional analogue of drain_update_tracer_x_budget.

Arguments

Type IntentOptional Attributes Name
integer, intent(in) :: nx
integer, intent(in) :: ny
integer, intent(in) :: nz
real(kind=wp), intent(in) :: iareaT(nx,ny)
real(kind=wp), intent(in) :: F(nx,ny+1,nz)
real(kind=wp), intent(in) :: w
real(kind=wp), intent(inout) :: hTr(nx,ny,nz)
real(kind=wp), intent(inout) :: budget_adv(nx,ny,nz)

private subroutine drain_wrap_centre(fld, nx, ny, nz, nx_phys, ny_phys, nghost, per_x, per_y, fold_n, no_wait)

Periodic wrap (+ north fold) of a cell-centred drain field. no_wait (optional): when .true. AND not folding, the periodic wrap is issued async on queue 1 without syncing, so a caller can batch several independent wraps (e.g. the pal/par/pa6 parabola triple) and !$acc wait(1) once. Ignored when fold_n (the fold reads the wrapped field, so the periodic wrap must complete first).

Arguments

Type IntentOptional Attributes Name
real(kind=wp), intent(inout) :: fld(nx,ny,nz)
integer, intent(in) :: nx
integer, intent(in) :: ny
integer, intent(in) :: nz
integer, intent(in) :: nx_phys
integer, intent(in) :: ny_phys
integer, intent(in) :: nghost
logical, intent(in) :: per_x
logical, intent(in) :: per_y
logical, intent(in) :: fold_n
logical, intent(in), optional :: no_wait

private subroutine drain_wrap_face_x(fld, nx, ny, nz, nx_phys, ny_phys, nghost, per_x, per_y, fold_n)

Periodic wrap (+ north fold) of an x-face drain field (nx+1,ny,nz).

Arguments

Type IntentOptional Attributes Name
real(kind=wp), intent(inout) :: fld(nx+1,ny,nz)
integer, intent(in) :: nx
integer, intent(in) :: ny
integer, intent(in) :: nz
integer, intent(in) :: nx_phys
integer, intent(in) :: ny_phys
integer, intent(in) :: nghost
logical, intent(in) :: per_x
logical, intent(in) :: per_y
logical, intent(in) :: fold_n

private subroutine drain_wrap_face_y(fld, nx, ny, nz, nx_phys, ny_phys, nghost, per_x, per_y, fold_n)

Periodic wrap (+ north fold) of a y-face drain field (nx,ny+1,nz).

Arguments

Type IntentOptional Attributes Name
real(kind=wp), intent(inout) :: fld(nx,ny+1,nz)
integer, intent(in) :: nx
integer, intent(in) :: ny
integer, intent(in) :: nz
integer, intent(in) :: nx_phys
integer, intent(in) :: ny_phys
integer, intent(in) :: nghost
logical, intent(in) :: per_x
logical, intent(in) :: per_y
logical, intent(in) :: fold_n

private pure subroutine drain_zero_3d(n1, n2, n3, fld)

fld = 0.

Arguments

Type IntentOptional Attributes Name
integer, intent(in) :: n1
integer, intent(in) :: n2
integer, intent(in) :: n3
real(kind=wp), intent(inout) :: fld(n1,n2,n3)

private pure subroutine gm_tracer_advect_x(grid, metrics, this, ms, uflux, dt, budget_w)

Zonal PPM tracer advection of every horizontally-advected tracer by the GM bolus flux uflux (the resolved path’s own kernel, tracer_advect_zonal_one_impl; heat/salt budgets weighted by budget_w). No OBC ghost override: continuity_gm_apply has closed every non-periodic edge face.

Arguments

Type IntentOptional Attributes Name
type(hgrid_t), intent(in) :: grid
type(ocean_metrics_t), intent(in) :: metrics
type(continuity_t), intent(inout) :: this
type(multilayer_state_t), intent(inout) :: ms
real(kind=wp), intent(in) :: uflux(grid%nx_total+1,grid%ny_total,ms%nz_ml)

GM zonal bolus transport (m^3/s).
real(kind=wp), intent(in) :: dt
real(kind=wp), intent(in) :: budget_w

private pure subroutine gm_tracer_advect_y(grid, metrics, this, ms, vflux, dt, budget_w)

Meridional twin of gm_tracer_advect_x.

Arguments

Type IntentOptional Attributes Name
type(hgrid_t), intent(in) :: grid
type(ocean_metrics_t), intent(in) :: metrics
type(continuity_t), intent(inout) :: this
type(multilayer_state_t), intent(inout) :: ms
real(kind=wp), intent(in) :: vflux(grid%nx_total,grid%ny_total+1,ms%nz_ml)

GM meridional bolus transport (m^3/s).
real(kind=wp), intent(in) :: dt
real(kind=wp), intent(in) :: budget_w

private subroutine pd_limit_meridional_impl(nx, ny, nz, dt, h_lim, iareaT, h_layer, mass_flux_y, theta, n_limited, v_cor)

Positive-definite per-donor outflux limiter — meridional (y) pass (P2). Mirror of pd_limit_zonal_impl; reads the post-zonal-apply h_layer (= h*), which is exactly the availability the second Lie pass must respect, and scales mass_flux_y so h_layer >= h_lim holds after continuity_apply_meridional. Cell (i,j,k) outflow = (north face j+1 when positive) + (south face j when negative); interior face j ∈ 2..ny scaled by its upwind donor’s θ (south cell j−1 when the face flux ≥ 0, else north cell j). See the zonal twin for the θ construction, ghost range, MPI-seam determinism, and the v1.1 v_cor re-matching rationale (MOM6 v_cor scaled by the same per-face θ so it stays consistent with the limited flux).

Arguments

Type IntentOptional Attributes Name
integer, intent(in) :: nx
integer, intent(in) :: ny
integer, intent(in) :: nz
real(kind=wp), intent(in) :: dt
real(kind=wp), intent(in) :: h_lim
real(kind=wp), intent(in) :: iareaT(nx,ny)
real(kind=wp), intent(in) :: h_layer(nx,ny,nz)
real(kind=wp), intent(inout) :: mass_flux_y(nx,ny+1,nz)
real(kind=wp), intent(inout) :: theta(nx,ny,nz)
integer, intent(inout) :: n_limited
real(kind=wp), intent(inout), optional :: v_cor(nx,ny+1,nz)

MOM6 v_cor capture; when present, re-scaled by the SAME per-face θ as mass_flux_y so it stays consistent with the limited flux the mom6-scheme corrector reads.

private subroutine pd_limit_zonal_impl(nx, ny, nz, dt, h_lim, iareaT, h_layer, mass_flux_x, theta, n_limited, u_cor)

Positive-definite per-donor outflux limiter — zonal (x) pass (P2, plan §P2 / decision D2). Scales the per-layer east-face mass fluxes DOWN so no donor cell loses more thickness than it holds above the floor h_lim: guarantees h_layer >= h_lim after continuity_apply_zonal, with ZERO mass created — outfluxes shrink, thickness is never inflated (the deliberate contrast with MOM6’s max(h, Angstrom) injection; the conservative borrow stays the backstop). Two device passes:

Read more…

Arguments

Type IntentOptional Attributes Name
integer, intent(in) :: nx
integer, intent(in) :: ny
integer, intent(in) :: nz
real(kind=wp), intent(in) :: dt
real(kind=wp), intent(in) :: h_lim
real(kind=wp), intent(in) :: iareaT(nx,ny)
real(kind=wp), intent(in) :: h_layer(nx,ny,nz)
real(kind=wp), intent(inout) :: mass_flux_x(nx+1,ny,nz)
real(kind=wp), intent(inout) :: theta(nx,ny,nz)
integer, intent(inout) :: n_limited
real(kind=wp), intent(inout), optional :: u_cor(nx+1,ny,nz)

MOM6 u_cor capture (transport-matched velocity); when present, re-scaled by the SAME per-face θ as mass_flux_x so it stays consistent with the limited flux the mom6-scheme corrector reads.

private pure subroutine renormalise_meridional_flux_to_vhbt(grid, metrics, this, ms, vhbt, dt, skip_walls, has_south, has_north, visc_rem, v_cor, use_por, por, use_open, open_f)

Apply a uniform per-face velocity correction so Σ_k mass_flux_y_layer(i, j, k) = vhbt(i, j) at every face. Mirror of renormalise_zonal_flux_to_uhbt. See that routine for the skip_walls and has_* (MPI seam, O4 fix) semantics.

Arguments

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

Read-only here (only the PPM edge buffers are consulted). intent(in) is load-bearing, not tidiness: the knob-off call sites pass one of THESE buffers as the inert por stand-in, and intent(in) turns “the callee never defines it” from a comment into a compiler-enforced invariant, so the argument association can never become aliasing.
type(multilayer_state_t), intent(inout) :: ms
real(kind=wp), intent(in) :: vhbt(:,:)
real(kind=wp), intent(in) :: dt

Outer-step dt (s), for the CFL bracket on dv.
logical, intent(in), optional :: skip_walls
logical, intent(in), optional :: has_south

Physical-edge flags (default .true. = single-rank behaviour, bit-identical); .false. at a y-decomposition seam forces the local wall-position face to be renormalised as interior.
logical, intent(in), optional :: has_north

Physical-edge flags (default .true. = single-rank behaviour, bit-identical); .false. at a y-decomposition seam forces the local wall-position face to be renormalised as interior.
real(kind=wp), intent(in), optional :: visc_rem(:,:,:)

Per-layer viscous remnant γ_k. Absent ⇒ γ ≡ 1, bit-identical. See renormalise_zonal_flux_to_uhbt for the full rationale.
real(kind=wp), intent(inout), optional :: v_cor(:,:,:)

MOM6 v_cor. Separate time-mean field, NEVER the prognostic — see the zonal twin and docs/MOM6_SPLIT_RK2_SPEC.md §5 trap 1.
logical, intent(in) :: use_por

Porous barriers active. .false. => por is never indexed and the arithmetic below stays byte-identical to the un-narrowed form.
real(kind=wp), intent(in) :: por(grid%nx_total,grid%ny_total+1,ms%nz_ml)

Layer-averaged open-area fraction at this stagger (nondim), read ONLY when use_por.

Read more…
logical, intent(in) :: use_open

z-level closed faces active. .false. => open_f is never indexed; byte-identical to the un-masked form. See the zonal twin for the full rationale.
real(kind=wp), intent(in) :: open_f(grid%nx_total,grid%ny_total+1,ms%nz_ml)

Per-layer 0/1 face-open mask at this stagger, read ONLY when use_open; enters the SAME weight wk = dx_cv * por * open. Knob-off callers pass h_face_right_y as the inert stand-in (never the (1,1,1) placeholder).

private pure subroutine renormalise_zonal_flux_to_uhbt(grid, metrics, this, ms, uhbt, dt, skip_walls, has_west, has_east, visc_rem, u_cor, use_por, por, use_open, open_f)

Apply a uniform per-face velocity correction so Σ_k mass_flux_x_layer(i, j, k) = uhbt(i, j) at every face. Helper for continuity_zonal_flux. skip_walls (default true) bypasses the physical-wall faces where mass_flux was zeroed; pass .false. for periodic axes so the wall faces (which carry real transport) are renormalised.

Arguments

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

Read-only here (only the PPM edge buffers are consulted). intent(in) is load-bearing, not tidiness: the knob-off call sites pass one of THESE buffers as the inert por stand-in, and intent(in) turns “the callee never defines it” from a comment into a compiler-enforced invariant, so the argument association can never become aliasing.
type(multilayer_state_t), intent(inout) :: ms
real(kind=wp), intent(in) :: uhbt(:,:)
real(kind=wp), intent(in) :: dt

Outer-step dt (s), for the CFL bracket on du.
logical, intent(in), optional :: skip_walls

When .true. (default), cycle the physical-wall faces. When .false. (periodic axis), include them in the renorm.
logical, intent(in), optional :: has_west

Physical-edge flags (default .true. = single-rank behaviour, bit-identical). .false. at an MPI seam: the local “wall position” face i = nghost+1 (west) / i = nghost+nx_phys+1 (east) is a REAL interior face carrying transport, so it MUST be renormalised like any interior face (O4 seam fix — the un-renormalised seam face left per-layer fluxes inconsistent with uhbt; the Eulerian-z h-rescale hid it in h while tracers rode the raw fluxes, breaking hTr/h at the seam).
logical, intent(in), optional :: has_east

Physical-edge flags (default .true. = single-rank behaviour, bit-identical). .false. at an MPI seam: the local “wall position” face i = nghost+1 (west) / i = nghost+nx_phys+1 (east) is a REAL interior face carrying transport, so it MUST be renormalised like any interior face (O4 seam fix — the un-renormalised seam face left per-layer fluxes inconsistent with uhbt; the Eulerian-z h-rescale hid it in h while tracers rode the raw fluxes, breaking hTr/h at the seam).
real(kind=wp), intent(in), optional :: visc_rem(:,:,:)

Per-layer viscous remnant γ_k weighting the barotropic increment (MOM6 u_cor = u + du·visc_rem; Jacobian duhdu = dy·h_marg·visc_rem). ABSENT ⇒ γ ≡ 1 ⇒ the historical uniform-du form, bit-identical.
real(kind=wp), intent(inout), optional :: u_cor(:,:,:)

MOM6’s u_cor return: the transport-consistent velocity u0 + du*gamma_k, i.e. the velocity that yields uhbt as the depth-integrated transport. This is a SEPARATE time-mean field — it must NEVER be the prognostic velocity. MOM6 writes it to u_av and evaluates the slow tendencies on it; writing it back into the prognostic kills the run by step 5. See docs/MOM6_SPLIT_RK2_SPEC.md §5 trap 1. Absent => flux-only, bit-identical.
logical, intent(in) :: use_por

Porous barriers active. .false. => por is never indexed and the arithmetic below stays byte-identical to the un-narrowed form.
real(kind=wp), intent(in) :: por(grid%nx_total+1,grid%ny_total,ms%nz_ml)

Layer-averaged open-area fraction at this stagger (nondim), read ONLY when use_por.

Read more…
logical, intent(in) :: use_open

z-level closed faces active (&vcoord_nml zfixed_closed_faces). .false. => open_f is never indexed and the arithmetic below stays byte-identical to the un-masked form.
real(kind=wp), intent(in) :: open_f(grid%nx_total+1,grid%ny_total,ms%nz_ml)

Per-layer 0/1 face-open mask at this stagger, read ONLY when use_open. It enters the SAME weight wk the porous fraction does – wk = dy_cu * por * open – which is the whole reason the barotropic transport is distributed over the OPEN layers only: a closed layer gets wk = 0, so it receives no du and contributes nothing to sum_h, and sum_k mass_flux = uhbt stays the exact fixed point the Newton solve iterates to.

Read more…

private subroutine tracer_advect(grid, metrics, this, ms, dt)

Test-only (no production caller): unsplit 2-D tracer advection, the reference oracle for split tracer_advect_zonal/_meridional. Per-layer PPM tracer advection — iterates over the tracer registry on the multilayer C-grid state and forwards each tracer’s hTr array to the flat-impl below. Extendable by construction: appending a new entry to ms%tracers(:) (BGC, sediment, passive scalar) drops it in without touching this routine. Per-tracer behaviour gates on tracer_t%do_horizontal_advection — set to .false. for tracers that should be diagnostic / forced externally.

Read more…

Arguments

Type IntentOptional Attributes Name
type(hgrid_t), intent(in) :: grid
type(ocean_metrics_t), intent(in) :: metrics
type(continuity_t), intent(inout) :: this
type(multilayer_state_t), intent(inout) :: ms
real(kind=wp), intent(in) :: dt

private pure subroutine tracer_advect_meridional(grid, metrics, this, ms, dt, bc)

Meridional half of the direction-split tracer advection. Reads the post-zonal h_layer (since continuity_apply_zonal has already updated h in the split flow) and the just-computed mass_flux_y_layer from continuity_meridional_flux.

Read more…

Arguments

Type IntentOptional Attributes Name
type(hgrid_t), intent(in) :: grid
type(ocean_metrics_t), intent(in) :: metrics
type(continuity_t), intent(inout) :: this
type(multilayer_state_t), intent(inout) :: ms
real(kind=wp), intent(in) :: dt
type(ocean_bc_state_t), intent(in), optional :: bc

private pure subroutine tracer_advect_meridional_one_impl(nx, ny, nz, dt, iareaT, wet_T, h, hTr, mass_flux_y, Tr_face_left_y, Tr_face_right_y, budget_adv, budget_w)

Meridional half of tracer_advect_one_impl. Same shape as the zonal impl, applied to y. In the split flow, h here is the post-zonal-apply thickness so the Tr = hTr/h reconstruction stays consistent with what continuity used in continuity_meridional_flux. iareaT = inv_dy on uniform metrics; mass_flux_y carries dx_cv. Mirror-T at land neighbours (C2); bit-identical for all-wet.

Read more…

Arguments

Type IntentOptional Attributes Name
integer, intent(in) :: nx
integer, intent(in) :: ny
integer, intent(in) :: nz
real(kind=wp), intent(in) :: dt
real(kind=wp), intent(in) :: iareaT(nx,ny)
real(kind=wp), intent(in) :: wet_T(nx,ny)
real(kind=wp), intent(in) :: h(nx,ny,nz)
real(kind=wp), intent(inout) :: hTr(nx,ny,nz)
real(kind=wp), intent(in) :: mass_flux_y(nx,ny+1,nz)
real(kind=wp), intent(inout) :: Tr_face_left_y(nx,ny+1,nz)
real(kind=wp), intent(inout) :: Tr_face_right_y(nx,ny+1,nz)
real(kind=wp), intent(inout), optional :: budget_adv(nx,ny,nz)

Per-cell meridional-advection budget accumulator (PSU·m or °C·m per cell). Added to the same array as the zonal half so the net entry covers both directions. Absent ⇒ inert.
real(kind=wp), intent(in), optional :: budget_w

Bookkeeping weight on the budget_adv increment; see the zonal impl. Absent ⇒ 1, bit-identical.

private pure subroutine tracer_advect_one_impl(nx, ny, nz, dt, iareaT, wet_T, h, hTr, mass_flux_x, mass_flux_y, Tr_face_left_x, Tr_face_right_x, Tr_face_left_y, Tr_face_right_y)

One-tracer PPM advection. Flat-impl: takes bare 3D arrays (no derived-type deref inside do-concurrent), so NVHPC stdpar handles it cleanly even for tracers stored in an array-of-derived-types registry.

Read more…

Arguments

Type IntentOptional Attributes Name
integer, intent(in) :: nx
integer, intent(in) :: ny
integer, intent(in) :: nz
real(kind=wp), intent(in) :: dt
real(kind=wp), intent(in) :: iareaT(nx,ny)
real(kind=wp), intent(in) :: wet_T(nx,ny)
real(kind=wp), intent(in) :: h(nx,ny,nz)
real(kind=wp), intent(inout) :: hTr(nx,ny,nz)
real(kind=wp), intent(in) :: mass_flux_x(nx+1,ny,nz)
real(kind=wp), intent(in) :: mass_flux_y(nx,ny+1,nz)
real(kind=wp), intent(inout) :: Tr_face_left_x(nx+1,ny,nz)
real(kind=wp), intent(inout) :: Tr_face_right_x(nx+1,ny,nz)
real(kind=wp), intent(inout) :: Tr_face_left_y(nx,ny+1,nz)
real(kind=wp), intent(inout) :: Tr_face_right_y(nx,ny+1,nz)

private pure subroutine tracer_advect_zonal(grid, metrics, this, ms, dt, bc)

Zonal half of the direction-split tracer advection. Mirrors tracer_advect but updates hTr using only the x-direction tracer mass flux. Companion to tracer_advect_meridional. Both are called interleaved with the continuity substeps by continuity_tracer_step_split to preserve CWC.

Read more…

Arguments

Type IntentOptional Attributes Name
type(hgrid_t), intent(in) :: grid
type(ocean_metrics_t), intent(in) :: metrics
type(continuity_t), intent(inout) :: this
type(multilayer_state_t), intent(inout) :: ms
real(kind=wp), intent(in) :: dt
type(ocean_bc_state_t), intent(in), optional :: bc

private pure subroutine tracer_advect_zonal_one_impl(nx, ny, nz, dt, iareaT, wet_T, h, hTr, mass_flux_x, Tr_face_left_x, Tr_face_right_x, budget_adv, budget_w)

Zonal half of tracer_advect_one_impl. Same three-pass pattern (PPM reconstruction → upwind pick into tracer mass flux → forward-Euler update) but only the x-direction half. Reads the input h for the Tr = hTr/h reconstruction. The tracer transport mass_flux_x·Tr_face inherits dy_cu, and the divergence closes with iareaT (= inv_dx on uniform). Mirror-T at land neighbours (C2): a held land column’s tracer is reflected to the local cell so the wet-side face value is unbiased; bit-identical for all-wet (wet_T≡1).

Read more…

Arguments

Type IntentOptional Attributes Name
integer, intent(in) :: nx
integer, intent(in) :: ny
integer, intent(in) :: nz
real(kind=wp), intent(in) :: dt
real(kind=wp), intent(in) :: iareaT(nx,ny)
real(kind=wp), intent(in) :: wet_T(nx,ny)
real(kind=wp), intent(in) :: h(nx,ny,nz)
real(kind=wp), intent(inout) :: hTr(nx,ny,nz)
real(kind=wp), intent(in) :: mass_flux_x(nx+1,ny,nz)
real(kind=wp), intent(inout) :: Tr_face_left_x(nx+1,ny,nz)
real(kind=wp), intent(inout) :: Tr_face_right_x(nx+1,ny,nz)
real(kind=wp), intent(inout), optional :: budget_adv(nx,ny,nz)

Per-cell accumulator for the horizontal-advection budget (same sign/units as hTr). When present, the zonal-flux divergence is added (+=) here after the prognostic update. Absent ⇒ inert (byte-identical to the pre-feature build).
real(kind=wp), intent(in), optional :: budget_w

Bookkeeping weight on the budget_adv increment (absent ⇒ 1, bit-identical). A caller writing OUTSIDE an RK2 stage (the GM operator, after the stage average) passes the reciprocal of the console’s per-step weight, 1/ocean_budget_stage_weight.