rdb_ocean_vmix Module

Holds the closure state for the ocean dynamical core’s vertical mixing kernel. Produces 3D kv, kt (and ks once the double-diffusion path lands) diffusivity fields at layer interfaces; the existing rdb_ocean_vdiff Thomas solver consumes them directly via the kv_source argument.

Two closure paths planned: * PP81 (Pacanowski-Philander 1981) — Richardson-number stability function; cheap, no boundary-layer detection, reasonable for stratified-interior dynamics. Live as of this branch. * KPP (Large, McWilliams & Doney 1994) — adds a Ri-bulk boundary-layer detector + shape functions + non-local transport term. Future work; the BL state slots (bl_depth, gamma_t, gamma_s) live here for that.

Interface convention: kv(:, :, k) lives at the bottom interface of layer k (between layers k-1 and k); kv(:, :, 1) is the bed (forced to zero, closed BC), kv(:, :, nz+1) is the surface (forced to zero, closed top). Same convention the rdb_ocean_vdiff impl assumes.

Convective adjustment (vmix_apply_convection) is a Brunt- Vaisala-triggered CONTRIBUTOR: where the interior N^2 < n2_thresh (dense-over-light), it raises kt -> max(kt, kd_conv) and kv -> max(kv, prandtl_conv*kd_conv), applied ONLY below the active surface boundary-layer depth (KPP bl_depth / EPBL mld own the BL response). kd_conv defaults to 1 m^2/s – five orders of magnitude above PP81’s Ri<0 clip ceiling (~1.01e-2 m^2/s) – and is admissible ONLY because vdiff_apply_tracers / vdiff_apply_momentum are backward-Euler (Thomas) tridiagonal solves: unconditionally stable, so kappa*dt/dz^2 >> 1 does not blow up. An explicit vertical diffusion could not use this coefficient. Reference: Brunt-Vaisala convective trigger as implemented in CVMix (Griffies et al., CVMix_convection), wired in MOM6 as USE_CVMix_CONVECTION / KD_CONV / PRANDTL_CONV / BV_SQR_CONV; the underlying “represent convection as a very large diapycnal diffusivity” idea traces to Cox (1984) and Marotzke (1991). Documented divergences from MOM6’s MOM_CVMix_conv: * D1 – the N^2 trigger differences ms%rho_layer (a potential density at the single global eos%p_ref, rdb_eos.F90 eos_wright_impl), not a locally-referenced (interface- pressure) density as MOM6 evaluates. Neglects thermobaricity; identical to what PP81 already assumes for the same expression. * D2 – uses max(kt, kd_conv) / max(kv, prandtl_conv*kd_conv) rather than MOM6’s additive Kd = Kd + kd_col. The two differ by at most the resolved interior Kd (<=1.01e-2, <=1% of kd_conv) – physically immaterial – but max() is an exact, idempotent floor, which matters because this runs every RK2 stage. SAFE ONLY because vmix_compute_pp81 ASSIGNS (not accumulates) kv/kt every stage when interior_closure == VMIX_INTERIOR_PP81 (the only reachable tag today): each stage starts clean so max() cannot ratchet. A future interior closure that does not rewrite kv/kt every stage would let this max() ratchet monotonically – document this dependency for whoever wires VMIX_INTERIOR_LARGE94 / VMIX_INTERIOR_CVMIX. * D3 – MOM6 calls calculate_CVMix_conv BEFORE energetic_PBL_get_MLD fills BLD for that step, so CVMix_conv can mask against a stale/zero BLD under EPBL (MOM6 emits a warning about this). Roundabout does not need the warning: vmix_apply_in_stage runs epbl_compute -> epbl_merge_into_kv_kt -> convection in the SAME stage, so epbl%mld is current when convection reads it. An improvement, not a gap. * D4 – writes kv and kt ONLY, never ks. Unlike kv/kt, ks is not rewritten by any per-stage contributor today (PP81 touches only kv/kt), so ks = max(ks, kd_conv) would ratchet monotonically and never relax. Salt convects for free once the ks <- kt split (a separate PR) lands downstream of this call, because that split copies the already-raised kt. Until then ks stays the constant pp81_kappa_bg background and salt genuinely does not convect – a pre-existing limitation this module does not worsen.

kd_max interaction: vmix_assemble is the single downstream ceiling gate (contract above). With the default kd_max = huge the convective 1 m^2/s passes through untouched; a caller who sets kd_max < kd_conv will silently cap the convective value – that is the gate doing its job, not a convection bug.


Uses

  • module~~rdb_ocean_vmix~~UsesGraph module~rdb_ocean_vmix rdb_ocean_vmix iso_fortran_env iso_fortran_env module~rdb_ocean_vmix->iso_fortran_env module~rdb_constants rdb_constants module~rdb_ocean_vmix->module~rdb_constants module~rdb_eos rdb_eos module~rdb_ocean_vmix->module~rdb_eos module~rdb_grid rdb_grid module~rdb_ocean_vmix->module~rdb_grid module~rdb_mem_report rdb_mem_report module~rdb_ocean_vmix->module~rdb_mem_report module~rdb_multilayer_state rdb_multilayer_state module~rdb_ocean_vmix->module~rdb_multilayer_state module~rdb_ocean_surface_flux rdb_ocean_surface_flux module~rdb_ocean_vmix->module~rdb_ocean_surface_flux module~rdb_ocean_surface_stress rdb_ocean_surface_stress module~rdb_ocean_vmix->module~rdb_ocean_surface_stress pic_types pic_types module~rdb_constants->pic_types module~rdb_eos->module~rdb_constants module~rdb_eos->module~rdb_grid 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_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_tracer rdb_tracer module~rdb_multilayer_state->module~rdb_tracer module~rdb_multilayer_state->pic_logger 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_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_scratch_3d rdb_scratch_3d module~rdb_ocean_surface_stress->module~rdb_scratch_3d module~rdb_efp->iso_fortran_env ieee_arithmetic ieee_arithmetic module~rdb_efp->ieee_arithmetic module~rdb_error_ring->pic_logger 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

Used by

  • module~~rdb_ocean_vmix~~UsedByGraph module~rdb_ocean_vmix rdb_ocean_vmix module~rdb_ocean_dyn rdb_ocean_dyn module~rdb_ocean_dyn->module~rdb_ocean_vmix module~rdb_ocean_setup rdb_ocean_setup module~rdb_ocean_setup->module~rdb_ocean_vmix module~rdb_ocean_setup->module~rdb_ocean_dyn module~rdb_ocean_state rdb_ocean_state module~rdb_ocean_setup->module~rdb_ocean_state module~rdb_ocean_state->module~rdb_ocean_vmix module~rdb_ocean_state->module~rdb_ocean_dyn proc~validate_config validate_config proc~validate_config->module~rdb_ocean_vmix 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_ocean_dyn module~rdb_ocean_engine->module~rdb_ocean_setup 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

Variables

Type Visibility Attributes Name Initial
integer, public, parameter :: BUOY_COEFFS_CONSTANT = 1

constant (DEFAULT): the scalar &ocean_ic_nml alpha_T/beta_S off the EOS handle, whatever the active EOS. Exactly right for eos = "linear" (they ARE that EOS’s coefficients); a constant stand-in under wright/roquet.
integer, public, parameter :: BUOY_COEFFS_EOS = 2

eos: eos_buoyancy_coeffs evaluated from the ACTIVE equation of state at each consumer’s own (T, S, p). Byte-identical to constant under eos = "linear".
integer, public, parameter :: BUOY_COEFFS_INVALID = 0

Unparsed spelling — configure_ocean_vmix fails loud on it rather than falling back (a silently-wrong α is a physics change with no symptom).
real(kind=wp), public, parameter :: HENYEY_KD_MIN_FRAC = 0.01_wp

Fraction of the scalar background tracer diffusivity used as the DEFAULT minimum diffusivity under the Henyey latitude scaling — MOM6 KD_MIN, whose documented default is 0.01*KD. Applied by vmix_resolve_kd_min when bkgnd_kd_min is left at its negative “unset” sentinel. PUBLIC so the unit tests derive the floor oracle from it rather than hard-coding 1e-7.
real(kind=wp), public, parameter :: HENYEY_MIN_SINLAT = 1.0e-10_wp

Floor on |sin(latitude)| used ONLY inside the N0_2Omega / |sin(latitude)| ratio of the Henyey factor (henyey_lat_factor_impl) to avoid the 1/0 singularity exactly at the equator — Henyey, Wright & Flatte (1986) JGR 91:8487; the constant-N0 simplification implemented here is Harrison & Hallberg (2008) JPO 38:1894. A fixed constant, not a namelist knob. The OUTER multiplication in the factor uses the TRUE (unfloored) |sin(latitude)|, so the assembled factor still goes smoothly to zero at the equator rather than blowing up. PUBLIC so the unit tests can DERIVE the poleward-clamp oracle from it instead of hard-coding a magic number that would silently rot if this value ever changed.
integer, public, parameter :: KPP_SW_ALL = 1

all_sw: charge B_0 with the full net heat flux (legacy).
integer, public, parameter :: KPP_SW_LV1 = 3

lv1_sw: subtract the SW that leaks below the top model layer.
integer, public, parameter :: KPP_SW_MXL = 2

mxl_sw (default): subtract the SW that leaks below h_b.
integer, public, parameter :: VMIX_INTERIOR_CVMIX = 3

CVMix-compatible (future).
integer, public, parameter :: VMIX_INTERIOR_LARGE94 = 2

Large et al. 1994 interior closure (future; not yet wired).
integer, public, parameter :: VMIX_INTERIOR_PP81 = 1

Pacanowski-Philander 1981 — implemented; the default.

Derived Types

type, public ::  ocean_vmix_t

Components

Type Visibility Attributes Name Initial
real(kind=wp), public, allocatable :: b0(:,:)

Surface buoyancy flux B_0 (m^2/s^3) the KPP overlay’s convective scale was built from, persisted per column on the SECOND pass (the one that uses the freshly-diagnosed bl_depth). Sign convention matches EPBL’s epbl%b0: > 0 stabilizing (heating / freshening), < 0 destabilizing (cooling / salting), and w_*^3 = max(0, -b0)*bl_depth. Diagnostic only — nothing reads it back into the closure; it exists so the two boundary-layer schemes’ surface forcing can be compared directly (test_ocean_buoyancy_flux). Zero until the first KPP overlay call.
real(kind=wp), public :: bkgnd_delta = 222.0_wp

Transition half-width (m): atan argument is (|z|-z0)/Delta.
logical, public :: bkgnd_henyey = .false.

Henyey, Wright & Flatte (1986) JGR 91:8487 latitude-dependent internal-wave factor, scaling the SCALAR background tracer diffusivities kt_bg/ks_bg by a horizontal-only factor L(phi) computed from geolatT (see henyey_lat_factor_impl for the exact form), with the result floored at bkgnd_kd_min:

Read more…
real(kind=wp), public :: bkgnd_henyey_max_lat = 95.0_wp

Latitude (degN) poleward of which the factor is reset to its equatorial (near-zero) floor. Deliberately > 90 by default so the clamp is INERT for any real latitude out of the box; lower it to activate the optional poleward cutoff.
real(kind=wp), public :: bkgnd_henyey_n0_2omega = 20.0_wp

Ratio of the assumed reference buoyancy frequency N0 to twice the planetary rotation rate (nondim). Physically N0 >> 2*Omega always, so this stays well above 1 (configure-time min=1 guard keeps the internal acosh argument in-domain even under a misconfigured value).
real(kind=wp), public :: bkgnd_kd_deep = 1.3e-4_wp

Deep-asymptote background tracer diffusivity (m^2/s).
real(kind=wp), public :: bkgnd_kd_min = -1.0_wp

Minimum background tracer diffusivity (m^2/s) under the Henyey latitude scaling — MOM6 KD_MIN, applied as max(Kd_min, Kd * L(phi)). Without it kd_bg would collapse toward zero at the equator (L(0 deg) = 0 exactly) and at the poleward max_lat clamp, instead of the documented behaviour “the Henyey profile is returned to the MINIMUM diffusivity”.

Read more…
real(kind=wp), public :: bkgnd_kd_sfc = 1.0e-5_wp

Surface-asymptote background tracer diffusivity (m^2/s). MOM6 BRYAN_LEWIS_C2-side default scale.
real(kind=wp), public :: bkgnd_prandtl = 1.0_wp

Background Prandtl number — Kv_bg = bkgnd_prandtl * Kd_bg (MOM6 ties the background viscosity to the background tracer diffusivity through a Prandtl factor). Default 1.0.
logical, public :: bkgnd_profile = .false.

Master switch for the depth-varying Bryan-Lewis background. Default .false. ⇒ vmix_assemble uses the scalar kv_bg/kt_bg/ks_bg floor exactly as before (bit-identical).
real(kind=wp), public :: bkgnd_z0 = 2500.0_wp

Transition-centre depth (m, positive down) where the profile reaches the (sfc+deep)/2 midpoint.
real(kind=wp), public, allocatable :: bl_depth(:,:)

KPP boundary-layer depth (m, positive down).
integer, public :: buoyancy_coeffs = BUOY_COEFFS_CONSTANT

BUOY_COEFFS_CONSTANT | BUOY_COEFFS_EOS. Read ON-DEVICE (inside vmix_kpp_overlay_impl’s do concurrent bodies), so it rides the same configure-precedes-enter_data contract as rho0 and the pp81_* scalars beside it.
real(kind=wp), public :: c_vt2 = 1.8_wp

Unresolved-turbulence coefficient for the V_t² term in the bulk-Ri denominator (LMD94 eq 23). Folded form of (C_v · √(-β_T) · √(c_s · ε)) / κ. Default 1.8 matches the standard LMD94 tuning. Set to 0 to disable V_t² — recovers the shear-only bulk-Ri sweep bit-identically; useful as a discriminator in unit tests. bl_depth(:, :) from the previous step seeds w_s for V_t², so on the very first call (h_b_lagged = 0) the w_* contribution evaluates to zero and the algorithm self-bootstraps.
logical, public :: conv_enable = .false.

Master switch. Requires use_closure + thermodynamics (validated at configure). Default off => bit-identical.
real(kind=wp), public :: conv_kd = 1.0_wp

Convective tracer diffusivity (m^2/s). MOM6 KD_CONV default 1.0 – ~1e5x the background, admissible only because vdiff is an unconditionally-stable backward-Euler solve.
real(kind=wp), public :: conv_n2_thresh = 0.0_wp

Trigger threshold on N^2 (s^-2). MOM6 BV_SQR_CONV default 0.0. Strict < so an exactly-neutral interface (N^2 == 0) does not trigger.
real(kind=wp), public :: conv_prandtl = 1.0_wp

Kv_conv = conv_prandtl * kd_conv. MOM6 PRANDTL_CONV default 1.0.
real(kind=wp), public :: cs_nonlocal = 6.3_wp

Non-local transport coefficient C_s (LMD94 eq 20). The full form 6.32·(1 - 0.5·exp(-σ/0.1)) is asymptotic; at σ > 0.1 it sits within 1% of the limit value 6.32, and the non-local flux is most important in that range. Using the limit value keeps the kernel branch-free and matches MOM6’s KPP_Cstar default.
logical, public :: ddiff_enable = .false.

Master switch (MOM6 USE_CVMIX_DDIFF). Requires thermodynamics (validated at configure). Default off.
real(kind=wp), public :: ddiff_exp1 = 1.0_wp

Inner (bracket) exponent of the fingering clamped form (CVMix DDIFF_EXP1).
real(kind=wp), public :: ddiff_exp2 = 3.0_wp

Outer exponent of the fingering clamped form (CVMix DDIFF_EXP2); exp1=1, exp2=3 is the Large et al. cubic.
real(kind=wp), public :: ddiff_kappa_s = 1.0e-4_wp

Leading salt-fingering salinity diffusivity K_f (m^2/s, CVMix KAPPA_DDIFF_S). K_T = 0.7*K_S (0.7 hard-wired).
real(kind=wp), public :: ddiff_mol_diff = 1.5e-6_wp

Molecular diffusivity scaling the convection branch (m^2/s, CVMix MOL_DIFF) – the molecular value, NOT a background eddy diffusivity.
real(kind=wp), public :: ddiff_param1 = 0.909_wp

MC76 diffusive-convection exterior coeff (CVMix KAPPA_DDIFF_PARAM1).
real(kind=wp), public :: ddiff_param2 = 4.6_wp

MC76 middle coeff (CVMix KAPPA_DDIFF_PARAM2).
real(kind=wp), public :: ddiff_param3 = -0.54_wp

MC76 interior coeff (CVMix KAPPA_DDIFF_PARAM3).
real(kind=wp), public :: ddiff_strat_param_max = 2.55_wp

R_rho salt-fingering cutoff (CVMix STRAT_PARAM_MAX). Above it fingering diffusivity is zero.
logical, public :: ddiff_use_k90 = .false.

Diffusive-convection form: .false. = Marmorino-Caldwell 1976 (MC76, default), .true. = Kelley 1990 (K90).
type(eos_t), public :: eos
real(kind=wp), public, allocatable :: gamma_s(:,:,:)

Non-local salinity flux (PSU·m/s) at interfaces.
real(kind=wp), public, allocatable :: gamma_t(:,:,:)

Non-local temperature flux (°C·m/s) at interfaces.
real(kind=wp), public :: hmix_fixed = 20.0_wp

Surface-band thickness (m) over which kv_ml_invz2 profile is active. MOM6 production default 20 m.
integer, public :: interior_closure = VMIX_INTERIOR_PP81

Interior mixing closure tag.
logical, public :: is_init = .false.

True between init and destroy. Prefer this to allocated(...) — tracks GPU device attachment too.
real(kind=wp), public, allocatable :: kd_bg(:,:,:)

Bryan-Lewis depth-profile tracer background (m^2/s) at interfaces. kv background floor = bkgnd_prandtl * kd_bg.
real(kind=wp), public :: kd_max = huge(1.0_wp)

Ceiling on tracer diffusivity (kt, ks) (m^2/s). Default huge = no clip. MOM6 Kd_max.
integer, public :: kd_smooth_iterations = 0

Number of 1-2-1 horizontal smoothing passes applied to kv/kt at each interior interface. Default 0 = off (no smoothing kernel runs, bit-identical). MOM6 Kd_smooth.
integer, public :: kpp_sw_method = KPP_SW_MXL

Shortwave-in-boundary-layer method for B_0 (MOM6 KPP_SHORTWAVE_METHOD). Only bites when penetrating SW is active (sf%has_sw); inert at sw_pen_frac = 0 ⇒ the default mxl_sw is bit-identical to the legacy path there.
real(kind=wp), public, allocatable :: ks(:,:,:)

Salinity eddy diffusivity.
real(kind=wp), public :: ks_bg = 1.0e-5_wp

Background floor on salinity diffusivity (m^2/s). Default matches pp81_kappa_bg.
real(kind=wp), public, allocatable :: kt(:,:,:)

Temperature eddy diffusivity.
real(kind=wp), public :: kt_bg = 1.0e-5_wp

Background floor on temperature diffusivity (m^2/s). Default matches pp81_kappa_bg.
real(kind=wp), public, allocatable :: kv(:,:,:)

Momentum (u, v) eddy viscosity at interfaces.
real(kind=wp), public :: kv_bg = 1.0e-4_wp

Background floor on momentum viscosity (m^2/s). Default matches pp81_nu_bg so the floor never raises a PP81 value (PP81 writes pp81_nu_bg + nu0·factor ≥ pp81_nu_bg).
logical, public :: kv_from_restart = .false.

PR-2 (bt-rem-from-av-rem review): set by ocean_state_restart_read (rdb_ocean_state.F90) immediately after a restart read, from the registry’s entry_found("vmix_kv") — .true. iff THIS read actually found vmix_kv in the checkpoint (an older checkpoint without the field, or a cold start, both leave it .false.). configure_ocean_lateral’s mandatory pp81_* config-copy (vmix_seed_backgrounds) reads it to decide whether to skip reseeding kv — see that routine’s docstring. Pure host bookkeeping: never device-mapped, never itself in the restart registry (it describes a read that already happened, not state to carry forward).
real(kind=wp), public :: kv_max = huge(1.0_wp)

Ceiling on momentum viscosity (m^2/s). Default huge = no clip (bit-identical). MOM6 Kd_max momentum analogue.
real(kind=wp), public :: kv_ml_invz2 = 0.0_wp

Extra kinematic viscosity (m²/s) inside the surface band of thickness hmix_fixed. Profile is (hmix_fixed / z)² where z is the distance from the surface to the interface. Mirrors MOM6’s KV_ML_INVZ2 + HMIX_FIXED combo — provides the dissipation that lets the wind drive thin surface layers without exciting a numerical eigenmode. Zero (default) means no addition; the kv field comes purely from PP81 / KPP.
logical, public :: p_top_in_eos = .false.

Mirror of &ocean_psurf_nml in_eos — whether multilayer_state_t%p_top carries a surface load that the in-situ EOS pressure must be measured down from. Seeded at configure from the same knob EPBL’s in_eos takes, so the two boundary-layer schemes build their pressure stacks the same way. Only read when buoyancy_coeffs == BUOY_COEFFS_EOS (p_top is the zero array when the knob is off anyway — this is the belt-and-braces gate EPBL already carries).
real(kind=wp), public :: pp81_alpha = 5.0_wp

Richardson-number multiplier in (1 + α*Ri). Original paper uses 5; some implementations use 4 or 10.
real(kind=wp), public :: pp81_kappa_bg = 1.0e-5_wp

Background tracer diffusivity (m^2/s).
real(kind=wp), public :: pp81_nu0 = 1.0e-2_wp

Numerator viscosity at Ri=0 (m^2/s).
real(kind=wp), public :: pp81_nu_bg = 1.0e-4_wp

Background interior viscosity (m^2/s).
real(kind=wp), public :: rho0 = 1035.0_wp

Boussinesq reference density (kg/m^3). Used in the N² = -g/ρ_0 * dρ/dz expression, in the friction velocity u_* = √(|τ|/ρ_0), and in the kinematic surface fluxes q_T = Q_heat/(ρ_0·cp) / q_S = Q_salt/ρ_0 that set the KPP surface buoyancy flux B_0.

Read more…
real(kind=wp), public :: ri_crit = 0.3_wp

Critical bulk Richardson number for KPP BL depth.
real(kind=wp), public :: shear2_floor = 1.0e-10_wp

Lower bound on |∂u/∂z|² to avoid Ri = N²/0 blow-up in quiescent water columns.
real(kind=wp), public, allocatable :: smooth_scratch(:,:,:)

Scratch (nx, ny, nz+1) for one smoothing pass.
logical, public :: use_closure = .false.

Master switch. When .false. the driver skips the compute step and rdb_ocean_vdiff falls back to its scalar K_v_* constants. Default off so all existing tests stay on the scalar path.
logical, public :: use_kpp = .false.

Enable KPP surface-boundary-layer scheme (future).
logical, public :: vmix_guard = .false.

Debug-gated negative/NaN guard. When .true. the assembly scans kv/kt/ks for a negative or NaN value and error stops (or returns a non-zero status via the testable path). Default off (cheap reduction skipped) so production runs pay nothing.

Type-Bound Procedures

procedure, public, non_overridable :: bytes => ocean_vmix_bytes
procedure, public, non_overridable :: destroy => ocean_vmix_destroy
procedure, public, non_overridable :: enter_data => ocean_vmix_enter_data
procedure, public, non_overridable :: exit_data => ocean_vmix_exit_data
procedure, public, non_overridable :: init => ocean_vmix_init
procedure, public, non_overridable :: seed_backgrounds => vmix_seed_backgrounds

Functions

public pure function bkgnd_henyey_conflicts_profile(bkgnd_henyey, bkgnd_profile) result(conflict)

.true. when both background schemes are selected at once.

Read more…

Arguments

Type IntentOptional Attributes Name
logical, intent(in) :: bkgnd_henyey

&ocean_vmix_nml bkgnd_henyey.
logical, intent(in) :: bkgnd_profile

&ocean_vmix_nml bkgnd_profile (Bryan-Lewis).

Return Value logical

public pure function henyey_lat_factor_impl(lat_deg, n0_2omega, max_lat) result(fac)

Henyey, Wright & Flatte (1986) JGR 91:8487 latitude dependence of the internal-wave-driven mixing rate, in the SIMPLIFIED constant-N0 form of Harrison & Hallberg (2008) JPO 38:1894 — the in-situ column stratification is replaced by a fixed reference N0, so the factor collapses to a pure function of latitude:

Read more…

Arguments

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

T-point latitude, degrees north (may be negative).
real(kind=wp), intent(in) :: n0_2omega

N0/(2Ω) reference-stratification ratio (nondim, ≥ 1).
real(kind=wp), intent(in) :: max_lat

Poleward cutoff latitude (degN, compared against |lat_deg|).

Return Value real(kind=wp)

public pure function kpp_surface_buoyancy_flux(alpha_T, beta_S, rho0, q_T_kin, q_S_kin) result(b0)

Surface buoyancy flux for the KPP overlay’s convective scale:

Read more…

Arguments

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

Dimensional linear-EOS sensitivities (kg/m^3 per degC / psu).
real(kind=wp), intent(in) :: beta_S

Dimensional linear-EOS sensitivities (kg/m^3 per degC / psu).
real(kind=wp), intent(in) :: rho0

Boussinesq reference density (kg/m^3).
real(kind=wp), intent(in) :: q_T_kin

Kinematic surface heat / salt fluxes (K m/s, psu m/s).
real(kind=wp), intent(in) :: q_S_kin

Kinematic surface heat / salt fluxes (K m/s, psu m/s).

Return Value real(kind=wp)

public pure function kpp_sw_method_is_implemented(name) result(ok)

Fail-loud predicate for kpp_sw_method — validate_config aborts on any string this rejects.

Arguments

Type IntentOptional Attributes Name
character(len=*), intent(in) :: name

Return Value logical

public pure function parse_buoyancy_coeffs(name) result(tag)

Map the &ocean_vmix_nml buoyancy_coeffs string to a BUOY_COEFFS_* tag; BUOY_COEFFS_INVALID for an unrecognised string, which configure_ocean_vmix turns into a fail-loud abort (never a silent fallback to the constants). The accepted set must match the nml_enum allowed= list in register_ocean_vmix.

Arguments

Type IntentOptional Attributes Name
character(len=*), intent(in) :: name

Return Value integer

public pure function parse_kpp_sw_method(name) result(tag)

Map the &ocean_thermo_nml kpp_sw_method string to a KPP_SW_* tag; -1 for an unrecognised string (fail-loud at validate_config).

Arguments

Type IntentOptional Attributes Name
character(len=*), intent(in) :: name

Return Value integer

public pure function vmix_interior_closure_is_implemented(code) result(ok)

.true. only for VMIX_INTERIOR_PP81 — the one interior closure with a real kernel today. VMIX_INTERIOR_LARGE94 and VMIX_INTERIOR_CVMIX are declared-but-unimplemented reservations (no kernel); selecting one would leave kv/kt stale (no interior mixing at all). interior_closure has no namelist key yet, so this is defence-in-depth (PR-6) for the next code/config consumer that sets it — the predicate is the gate.

Arguments

Type IntentOptional Attributes Name
integer, intent(in) :: code

Return Value logical

public pure function vmix_resolve_kd_min(kd_min, kt_bg) result(kd_min_eff)

Resolve the bkgnd_kd_min “unset” sentinel to the reference default HENYEY_KD_MIN_FRAC * kt_bg (MOM6 KD_MIN, default 0.01*KD). A NEGATIVE kd_min means “not set by the user”; zero and positive values are taken literally, so kd_min = 0 is a legal way to ask for no floor at all.

Read more…

Arguments

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

The raw bkgnd_kd_min knob; negative = unset sentinel.
real(kind=wp), intent(in) :: kt_bg

Scalar background tracer diffusivity the default is a fraction of.

Return Value real(kind=wp)

private pure function ocean_vmix_bytes(this) result(nbytes)

Counted allocatable footprint of the vertical mixing 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(ocean_vmix_t), intent(in) :: this

Return Value integer(kind=int64)


Subroutines

public pure subroutine vmix_add_kv_ml_invz2(grid, this, ms)

Augment this%kv with an extra near-surface viscosity:

Read more…

Arguments

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

public pure subroutine vmix_apply_convection(grid, this, ms, bl_depth)

Brunt-Vaisala-triggered convective adjustment – a CONTRIBUTOR into kv/kt (see the module docstring for the full physics and the D1-D4 divergences from MOM6’s MOM_CVMix_conv). Early- returns when conv_enable is off (bit-identical). bl_depth (m, positive down) is the caller’s active surface-boundary-layer depth: pass this%bl_depth under KPP / no BL scheme, epbl%mld under EPBL – both are permanently-zero, device-resident fields when their owning scheme is off, so either is a free “no BL scheme -> mask nothing” default (z_int(k) >= 0 for every interior interface).

Arguments

Type IntentOptional Attributes Name
type(hgrid_t), intent(in) :: grid
type(ocean_vmix_t), intent(inout) :: this
type(multilayer_state_t), intent(in) :: ms
real(kind=wp), intent(in) :: bl_depth(:,:)

public pure subroutine vmix_apply_kpp_overlay(grid, this, ms, ss, sf)

Thin host-side shim over vmix_kpp_overlay_impl — same signature as before this PR, so the call site in rdb_ocean_dyn.F90 is unchanged. Selects the shortwave irradiance source HOST-SIDE (sf%sw_from_qsw): the PR-12 q_sw component is allocated only under use_components, so it is only ever passed on the branch guarded by the host flag (validate_config forces enable_components when sw_source="q_sw", making this total). sf%has_sw is passed as the sw_active gate — false ⇒ the _impl’s B_0 reduces to the unmodified legacy source line, bit-for-bit.

Read more…

Arguments

Type IntentOptional Attributes Name
type(hgrid_t), intent(in) :: grid
type(ocean_vmix_t), intent(inout) :: this
type(multilayer_state_t), intent(in) :: ms
type(ocean_surface_stress_t), intent(in) :: ss
type(ocean_surface_flux_t), intent(in) :: sf

public pure subroutine vmix_apply_nonlocal_tendencies(grid, this, ms, dt)

Apply the KPP non-local (counter-gradient) tracer tendency computed by vmix_apply_kpp_overlay. Updates hT, hS from the divergence of gamma_t / gamma_s at interfaces:

Read more…

Arguments

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

public subroutine vmix_assemble(grid, this, ms, geolat, status)

The single downstream gate of the vmix diffusivity assembly — the set_diffusivity-style stage that every interior and overlay closure feeds into before vdiff consumes kv/kt/ks.

Read more…

Arguments

Type IntentOptional Attributes Name
type(hgrid_t), intent(in) :: grid
type(ocean_vmix_t), intent(inout) :: this
type(multilayer_state_t), intent(in) :: ms
real(kind=wp), intent(in), optional :: geolat(grid%nx_total,grid%ny_total)

T-point geographic latitude (degN, ocean_metrics_t%geolatT — the T stagger specifically; geolatBu is corner-shaped (nx+1, ny+1) and is NOT interchangeable here). Explicit-shape so the value flows into vmix_assemble_clip_henyey_impl’s do concurrent without a descriptor walk. Only read when bkgnd_henyey. Required in that case — absent then error stops (a caller forgot to thread metrics through). Every production call site (vmix_apply_in_stage) has metrics in scope and always passes it; test harnesses that never enable bkgnd_henyey may omit it.
integer, intent(out), optional :: status

0 = ok; 1 = guard tripped (negative or NaN K). Only written when vmix_guard is on. When absent and the guard trips, the routine error stops instead.

public pure subroutine vmix_assemble_clip_henyey_impl(nx, ny, nzp1, kv, kt, ks, kv_bg, kt_bg, ks_bg, kv_max, kd_max, n0_2omega, henyey_max_lat, kd_min, geolat)

vmix_assemble_clip_impl with the Henyey latitude factor (henyey_lat_factor_impl) scaling the two SCALAR TRACER background floors, each then floored at the minimum diffusivity:

Read more…

Arguments

Type IntentOptional Attributes Name
integer, intent(in) :: nx
integer, intent(in) :: ny
integer, intent(in) :: nzp1
real(kind=wp), intent(inout) :: kv(nx,ny,nzp1)
real(kind=wp), intent(inout) :: kt(nx,ny,nzp1)
real(kind=wp), intent(inout) :: ks(nx,ny,nzp1)
real(kind=wp), intent(in) :: kv_bg
real(kind=wp), intent(in) :: kt_bg
real(kind=wp), intent(in) :: ks_bg
real(kind=wp), intent(in) :: kv_max
real(kind=wp), intent(in) :: kd_max
real(kind=wp), intent(in) :: n0_2omega

N0/(2Ω) ratio (nondim) and the poleward cutoff latitude (degN).
real(kind=wp), intent(in) :: henyey_max_lat

N0/(2Ω) ratio (nondim) and the poleward cutoff latitude (degN).
real(kind=wp), intent(in) :: kd_min

Minimum background tracer diffusivity (m^2/s), ALREADY resolved through vmix_resolve_kd_min — this kernel never sees the negative sentinel.
real(kind=wp), intent(in) :: geolat(nx,ny)

T-point latitude (degN) — ocean_metrics_t%geolatT, the T stagger. The corner field geolatBu is (nx+1, ny+1) and is NOT a substitute.

public pure subroutine vmix_bkgnd_fill_impl(nx, ny, nzp1, kd_bg, h_layer, kd_sfc, kd_deep, z0, delta)

Fill the per-interface Bryan & Lewis (1979) JGR 84:2503 background tracer-diffusivity profile from the CURRENT column interface depths.

Read more…

Arguments

Type IntentOptional Attributes Name
integer, intent(in) :: nx
integer, intent(in) :: ny
integer, intent(in) :: nzp1
real(kind=wp), intent(inout) :: kd_bg(nx,ny,nzp1)
real(kind=wp), intent(in) :: h_layer(nx,ny,nzp1-1)
real(kind=wp), intent(in) :: kd_sfc

Surface / deep asymptotes (m^2/s), transition centre depth (m), transition half-width (m).
real(kind=wp), intent(in) :: kd_deep

Surface / deep asymptotes (m^2/s), transition centre depth (m), transition half-width (m).
real(kind=wp), intent(in) :: z0

Surface / deep asymptotes (m^2/s), transition centre depth (m), transition half-width (m).
real(kind=wp), intent(in) :: delta

Surface / deep asymptotes (m^2/s), transition centre depth (m), transition half-width (m).

public pure subroutine vmix_compute_pp81(grid, this, ms)

Pacanowski-Philander (1981) Richardson-number closure. Inlined kernel — keeps the do concurrent body adjacent to its derived-type accesses. We tried the outer-shim + _impl pattern but NVHPC’s stdpar codegen produced more descriptor-marshalling memcpys at the shim boundary than the direct-access form generates inside the kernel. Direct ms%foo(i,j,k) access is what other working hot-path kernels (continuity, coriolis_adv, the barotropic substep) use; the shim pattern was an experiment that didn’t help here.

Arguments

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

public pure subroutine vmix_split_kd_heat_salt(grid, this, ms)

Analogue of MOM6’s heat/salt diffusivity split. Derives the per-tracer diffusivities from the assembled interior/boundary- layer diffusivity: Kd_heat = Kd_int + Kd_extra_T -> kt Kd_salt = Kd_int + Kd_extra_S -> ks kt holds Kd_int on entry (every contributor writes it). No double-diffusion contributor exists yet, so Kd_extra_{T,S} = 0 and the split reduces to ks := kt => bit-identical. PR-33 extends this to ks = kt + kd_extra_s ; kt = kt + kd_extra_t BOTH computed from the SAME pre-split kt – read kt into a local before writing it, or the kt update poisons the ks update.

Read more…

Arguments

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

private subroutine ocean_vmix_destroy(this)

Arguments

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

private subroutine ocean_vmix_enter_data(this)

Arguments

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

private subroutine ocean_vmix_enter_data_impl(this)

Arguments

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

private subroutine ocean_vmix_exit_data(this)

Arguments

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

private subroutine ocean_vmix_exit_data_impl(this)

Arguments

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

private subroutine ocean_vmix_init(this, grid, nz_ml)

Allocate the kv / kt / ks diffusivity fields at layer interfaces. Default values: kv = kv_bg (= pp81_nu_bg), kt = ks = kt_bg (= pp81_kappa_bg) at interior interfaces; boundary interfaces k=1 and k=nz+1 are zeroed (closed BC). BL fields (bl_depth, gamma_*) seed at zero. ks is allocated so the assembly stage (vmix_assemble) can floor/clip it; vmix_split_kd_heat_salt derives its live value from kt every stage, and vdiff_apply_tracers consumes it for salinity + every passive tracer.

Read more…

Arguments

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

private pure subroutine vmix_assemble_clip_impl(nx, ny, nzp1, kv, kt, ks, kv_bg, kt_bg, ks_bg, kv_max, kd_max)

Floor + ceiling on interior interfaces k = 2..nzp1-1. Explicit- shape args so the do concurrent stays descriptor-walk free.

Arguments

Type IntentOptional Attributes Name
integer, intent(in) :: nx
integer, intent(in) :: ny
integer, intent(in) :: nzp1
real(kind=wp), intent(inout) :: kv(nx,ny,nzp1)
real(kind=wp), intent(inout) :: kt(nx,ny,nzp1)
real(kind=wp), intent(inout) :: ks(nx,ny,nzp1)
real(kind=wp), intent(in) :: kv_bg
real(kind=wp), intent(in) :: kt_bg
real(kind=wp), intent(in) :: ks_bg
real(kind=wp), intent(in) :: kv_max
real(kind=wp), intent(in) :: kd_max

private pure subroutine vmix_assemble_clip_profile_impl(nx, ny, nzp1, kv, kt, ks, kd_bg, prandtl, kv_max, kd_max)

C7 floor + ceiling: same as vmix_assemble_clip_impl but the tracer floor is the per-interface Bryan-Lewis kd_bg field and the momentum floor is prandtl·kd_bg (MOM6 background Prandtl tie). Interior interfaces k = 2..nzp1-1.

Arguments

Type IntentOptional Attributes Name
integer, intent(in) :: nx
integer, intent(in) :: ny
integer, intent(in) :: nzp1
real(kind=wp), intent(inout) :: kv(nx,ny,nzp1)
real(kind=wp), intent(inout) :: kt(nx,ny,nzp1)
real(kind=wp), intent(inout) :: ks(nx,ny,nzp1)
real(kind=wp), intent(in) :: kd_bg(nx,ny,nzp1)
real(kind=wp), intent(in) :: prandtl
real(kind=wp), intent(in) :: kv_max
real(kind=wp), intent(in) :: kd_max

private pure subroutine vmix_convection_impl(nx, ny, nzp1, kv, kt, rho_layer, h_layer, bl_depth, kd_conv, prandtl_conv, n2_thresh, rho0)

Flat, explicit-shape kernel (the “vmix incident” module – assumed-shape dummies here produced 1.4M per-launch descriptor- walk memcpys; explicit-shape only, no exceptions). Interior interfaces k = 2..nzp1-1 only; boundary interfaces k=1 (bed) and k=nzp1 (surface) stay at the closed-BC zero and are never touched (same contract as PP81 / vmix_add_kv_ml_invz2).

Read more…

Arguments

Type IntentOptional Attributes Name
integer, intent(in) :: nx
integer, intent(in) :: ny
integer, intent(in) :: nzp1
real(kind=wp), intent(inout) :: kv(nx,ny,nzp1)
real(kind=wp), intent(inout) :: kt(nx,ny,nzp1)
real(kind=wp), intent(in) :: rho_layer(nx,ny,nzp1-1)
real(kind=wp), intent(in) :: h_layer(nx,ny,nzp1-1)
real(kind=wp), intent(in) :: bl_depth(nx,ny)
real(kind=wp), intent(in) :: kd_conv
real(kind=wp), intent(in) :: prandtl_conv
real(kind=wp), intent(in) :: n2_thresh
real(kind=wp), intent(in) :: rho0

private pure subroutine vmix_guard_impl(nx, ny, nzp1, kv, kt, ks, bad_count)

Count negative or NaN diffusivities across interior interfaces. A reduction over fresh scratch — uses !$acc parallel loop reduction (the project rule for device reductions; bare count/sum over device data can silently return 0).

Arguments

Type IntentOptional Attributes Name
integer, intent(in) :: nx
integer, intent(in) :: ny
integer, intent(in) :: nzp1
real(kind=wp), intent(in) :: kv(nx,ny,nzp1)
real(kind=wp), intent(in) :: kt(nx,ny,nzp1)
real(kind=wp), intent(in) :: ks(nx,ny,nzp1)
integer, intent(out) :: bad_count

private pure subroutine vmix_kpp_overlay_impl(grid, this, ms, ss, sf, nx_arg, ny_arg, nz_arg, sw_src, temp_h, salt_h, sw_active, sw_pen_frac, sw_R, sw_zeta1, sw_zeta2, sw_method)

KPP boundary-layer overlay on top of the interior closure already in this%kv / this%kt. Phase 1 was shear-driven only; Phase 2 added the convective velocity scale w_* and γ_T/γ_S non-local transport; Phase 3 (this revision) adds the V_t² unresolved-turbulence term in the bulk-Ri denominator (LMD94 eq 23). Surface BL is now feature- complete except for Langmuir / Stokes enhancement.

Read more…

Arguments

Type IntentOptional Attributes Name
type(hgrid_t), intent(in) :: grid
type(ocean_vmix_t), intent(inout) :: this
type(multilayer_state_t), intent(in) :: ms
type(ocean_surface_stress_t), intent(in) :: ss
type(ocean_surface_flux_t), intent(in) :: sf
integer, intent(in) :: nx_arg

Grid extents — declared before sw_src (decl-order hook).
integer, intent(in) :: ny_arg

Grid extents — declared before sw_src (decl-order hook).
integer, intent(in) :: nz_arg

Grid extents — declared before sw_src (decl-order hook).
real(kind=wp), intent(in) :: sw_src(nx_arg,ny_arg)

Caller-selected irradiance source (sf%Q_heat or sf%q_sw), explicit-shape (per-RK2-stage kernel: no assumed-shape waiver).
real(kind=wp), intent(in) :: temp_h(nx_arg,ny_arg,nz_arg)

hTr of the temperature tracer (degC·m), flattened off the registry by the shim. Read ONLY under buoyancy_coeffs == BUOY_COEFFS_EOS, and only at k = nz.
real(kind=wp), intent(in) :: salt_h(nx_arg,ny_arg,nz_arg)

hTr of the salinity tracer (PSU·m). Same contract.
logical, intent(in) :: sw_active

Host-side sf%has_sw gate — false ⇒ B_0 uses the unmodified legacy source line (bit-identity).
real(kind=wp), intent(in) :: sw_pen_frac

Two-band SW parameters (from sf), by value.
real(kind=wp), intent(in) :: sw_R

Two-band SW parameters (from sf), by value.
real(kind=wp), intent(in) :: sw_zeta1

Two-band SW parameters (from sf), by value.
real(kind=wp), intent(in) :: sw_zeta2

Two-band SW parameters (from sf), by value.
integer, intent(in) :: sw_method

KPP_SW_ALL | KPP_SW_MXL | KPP_SW_LV1.

private pure subroutine vmix_seed_backgrounds(this, skip_kv)

Seed kv_bg/kt_bg/ks_bg and the kv/kt/ks/kd_bg arrays from the current pp81_nu_bg/pp81_kappa_bg fields, then zero the closed-BC boundary interfaces on kv/ks. Extracted out of ocean_vmix_init (§2E structural invariant: kv_bg == pp81_nu_bg, kt_bg == ks_bg == pp81_kappa_bg, so vmix_assemble’s background floor is a no-op for the shipped PP81 closure) so BOTH the initial seed and the &ocean_vmix_nml pp81_* config-copy re-derive consistently — a bare field copy without this call would leave the assembly floor clamped against the OLD (type-default) background even after a user sets a new one. Requires kv/kt/ks/kd_bg already allocated (true after init; the config-copy call runs strictly after init_from_config).

Arguments

Type IntentOptional Attributes Name
class(ocean_vmix_t), intent(inout) :: this
logical, intent(in), optional :: skip_kv

PR-2 (bt-rem-from-av-rem, fixed per review): default .false. — the FULL seed always runs (scalars + kv/kt/ks/kd_bg arrays + the kv/ks boundary zero), exactly the historical behaviour. The config-copy call site (configure_ocean_lateral, AFTER engine_setup’s restart read) passes skip_kv = state%vmix%kv_from_restart — .true. ONLY when THIS read actually found vmix_kv in the checkpoint (an older checkpoint without the field, or a cold start, both leave kv_from_restart = .false., so the array still reseeds normally and the run is never left with an uninitialised kv). When skipped, kv is left EXACTLY as the restart read wrote it — no reseed, no boundary re-zero — because a checkpointed kv is a CARRIED field (visc_rem_precompute reads the PREVIOUS stage’s kv before this stage recomputes it) and re-zeroing its boundary rows is not provably idempotent: nothing in this tree asserts every kv-writing closure (PP81/KPP/EPBL/kappa-shear/tidal-mixing/convective adjustment, vmix_assemble) keeps kv(:,:,1) / kv(:,:,nz+1) at exactly 0 throughout a run, so re-asserting it here could diverge a restored run from the continued one it must match bitwise. Restoring the checkpoint verbatim is the only choice that is unconditionally correct.

Read more…

private pure subroutine vmix_smooth_121_impl(nx, ny, nzp1, fld, scratch, wet_mask)

One in-plane 1-2-1 horizontal smoothing pass on interior interfaces k = 2..nzp1-1. Wet-mask aware: contributions from dry neighbours (wet_mask == 0) are excluded and the 9-point stencil weight is renormalised over the wet cells only. A dry centre column (wet_mask(i,j) == 0) is left unchanged — no leakage into or out of dry cells.

Read more…

Arguments

Type IntentOptional Attributes Name
integer, intent(in) :: nx
integer, intent(in) :: ny
integer, intent(in) :: nzp1
real(kind=wp), intent(inout) :: fld(nx,ny,nzp1)
real(kind=wp), intent(inout) :: scratch(nx,ny,nzp1)
real(kind=wp), intent(in) :: wet_mask(nx,ny)

1 = wet, 0 = dry. Shape (nx, ny).

private pure subroutine vmix_split_ddiff_eos_impl(nx, ny, nzp1, kt, ks, temp_h, salt_h, h_layer, p_top, eos, rho0, p_top_in_eos, strat_param_max, kappa_s, exp1, exp2, param1, param2, param3, mol_diff, use_k90)

buoyancy_coeffs = "eos" twin of vmix_split_ddiff_impl — the SAME CVMix closed forms, the same branch structure, the same outputs; the only change is where α and β come from.

Read more…

Arguments

Type IntentOptional Attributes Name
integer, intent(in) :: nx
integer, intent(in) :: ny
integer, intent(in) :: nzp1
real(kind=wp), intent(inout) :: kt(nx,ny,nzp1)
real(kind=wp), intent(inout) :: ks(nx,ny,nzp1)
real(kind=wp), intent(in) :: temp_h(nx,ny,nzp1-1)

hTr of temperature (degC·m).
real(kind=wp), intent(in) :: salt_h(nx,ny,nzp1-1)

hTr of salinity (PSU·m).
real(kind=wp), intent(in) :: h_layer(nx,ny,nzp1-1)
real(kind=wp), intent(in) :: p_top(nx,ny)

Surface load (Pa) — multilayer_state_t%p_top, the zero array unless &ocean_psurf_nml in_eos.
type(eos_t), intent(in) :: eos

Active EOS handle, by value.
real(kind=wp), intent(in) :: rho0

Boussinesq reference density for the hydrostatic accumulation (the one configured ρ₀ of record, via vmix%rho0).
logical, intent(in) :: p_top_in_eos

Whether to seed the stack from p_top — mirrors EPBL’s gate.
real(kind=wp), intent(in) :: strat_param_max
real(kind=wp), intent(in) :: kappa_s
real(kind=wp), intent(in) :: exp1
real(kind=wp), intent(in) :: exp2
real(kind=wp), intent(in) :: param1
real(kind=wp), intent(in) :: param2
real(kind=wp), intent(in) :: param3
real(kind=wp), intent(in) :: mol_diff
logical, intent(in) :: use_k90

private pure subroutine vmix_split_ddiff_impl(nx, ny, nzp1, kt, ks, temp_h, salt_h, h_layer, alpha_T, beta_S, strat_param_max, kappa_s, exp1, exp2, param1, param2, param3, mol_diff, use_k90)

Double-diffusion split. Replaces ks := kt with the asymmetric ks = kt_pre + kd_extra_s ; kt = kt_pre + kd_extra_t both from the SAME pre-split kt (read into kt_pre before either write). CVMix cvmix_coeffs_ddiff algebra (Large et al. 1994 fingering; Marmorino-Caldwell 1976 / Kelley 1990 convection).

Read more…

Arguments

Type IntentOptional Attributes Name
integer, intent(in) :: nx
integer, intent(in) :: ny
integer, intent(in) :: nzp1
real(kind=wp), intent(inout) :: kt(nx,ny,nzp1)
real(kind=wp), intent(inout) :: ks(nx,ny,nzp1)
real(kind=wp), intent(in) :: temp_h(nx,ny,nzp1-1)

Temperature thickness-integral hTr = Th_layer (degCm).
real(kind=wp), intent(in) :: salt_h(nx,ny,nzp1-1)

Salinity thickness-integral hTr = Sh_layer (PSUm).
real(kind=wp), intent(in) :: h_layer(nx,ny,nzp1-1)
real(kind=wp), intent(in) :: alpha_T

|d rho/dT| and d rho/dS (kg/m^3 per degC / per PSU).
real(kind=wp), intent(in) :: beta_S

|d rho/dT| and d rho/dS (kg/m^3 per degC / per PSU).
real(kind=wp), intent(in) :: strat_param_max
real(kind=wp), intent(in) :: kappa_s
real(kind=wp), intent(in) :: exp1
real(kind=wp), intent(in) :: exp2
real(kind=wp), intent(in) :: param1
real(kind=wp), intent(in) :: param2
real(kind=wp), intent(in) :: param3
real(kind=wp), intent(in) :: mol_diff
logical, intent(in) :: use_k90

private pure subroutine vmix_split_kd_heat_salt_impl(nx, ny, nzp1, kt, ks)

Explicit-shape args so the do concurrent stays descriptor-walk free (assumed-shape dummies in a do concurrent make NVHPC walk descriptors per launch – this runs every stage).

Arguments

Type IntentOptional Attributes Name
integer, intent(in) :: nx
integer, intent(in) :: ny
integer, intent(in) :: nzp1
real(kind=wp), intent(in) :: kt(nx,ny,nzp1)
real(kind=wp), intent(inout) :: ks(nx,ny,nzp1)