rdb_ocean_wave_speed Module

Per-column first-baroclinic internal gravity-wave speed cg1 (m/s) and first-mode Rossby radius Rd (m), plus Rd/dx (the GM/Redi/MEKE resolution ratio). Solves the rigid-lid Sturm-Liouville eigenproblem discretised from layer thicknesses and per-interface reduced gravities gprime = (g/rho0)*max(0,drho) (rho_layer authoritative); largest c^2 via a fixed-budget Sturm-count bisection (no early exit -> warp-divergence-free). Reference: Chelton et al. (1998).

Vertical ordering (load-bearing): rdb is bottom-up (k=1 bed, k=nz surface); the eigensolve is surface-down so the column kernel FLIPS on gather (local k_loc=1 is the surface layer).

Diagnostic, default off (&ocean_wavespeed_nml enable=.false.) ⇒ kernel never called, bit-identical.


Uses

  • module~~rdb_ocean_wave_speed~~UsesGraph module~rdb_ocean_wave_speed rdb_ocean_wave_speed iso_fortran_env iso_fortran_env module~rdb_ocean_wave_speed->iso_fortran_env module~rdb_constants rdb_constants module~rdb_ocean_wave_speed->module~rdb_constants module~rdb_grid rdb_grid module~rdb_ocean_wave_speed->module~rdb_grid module~rdb_mem_report rdb_mem_report module~rdb_ocean_wave_speed->module~rdb_mem_report module~rdb_multilayer_state rdb_multilayer_state module~rdb_ocean_wave_speed->module~rdb_multilayer_state module~rdb_ocean_metrics rdb_ocean_metrics module~rdb_ocean_wave_speed->module~rdb_ocean_metrics 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_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_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_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_fold rdb_ocean_fold module~rdb_ocean_metrics->module~rdb_ocean_fold module~rdb_ocean_status rdb_ocean_status 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_efp->iso_fortran_env ieee_arithmetic ieee_arithmetic 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_fold->module~rdb_constants 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_wave_speed~~UsedByGraph module~rdb_ocean_wave_speed rdb_ocean_wave_speed module~rdb_ocean_dyn rdb_ocean_dyn module~rdb_ocean_dyn->module~rdb_ocean_wave_speed module~rdb_ocean_meke rdb_ocean_meke module~rdb_ocean_dyn->module~rdb_ocean_meke module~rdb_ocean_varmix rdb_ocean_varmix module~rdb_ocean_dyn->module~rdb_ocean_varmix module~rdb_ocean_meke->module~rdb_ocean_wave_speed module~rdb_ocean_meke->module~rdb_ocean_varmix module~rdb_ocean_state rdb_ocean_state module~rdb_ocean_state->module~rdb_ocean_wave_speed module~rdb_ocean_state->module~rdb_ocean_dyn module~rdb_ocean_state->module~rdb_ocean_meke module~rdb_ocean_state->module~rdb_ocean_varmix module~rdb_ocean_varmix->module~rdb_ocean_wave_speed 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_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
real(kind=wp), private, parameter :: C2_SCALE = 1.0_wp/(4096.0_wp*4096.0_wp)

Per-row determinant rescale s (slows det growth between rows).
real(kind=wp), private, parameter :: F_DENOM_FLOOR = 1.0e-10_wp

Floor on the Rd denominator (guards f = beta = 0).
real(kind=wp), private, parameter :: I_RESCALE = 1.0_wp/(1024.0_wp**4)
real(kind=wp), private, parameter :: LAM_SEED = 1.0e-20_wp

Tiny lower seed for the doubling prelude (below any physical eigenvalue).
integer, private, parameter :: MAX_DBL = 128

Cap on the Sturm-count doubling prelude (data-dependent trip count, but O(1) det evals; the inner bisection stays fixed).
integer, private, parameter :: MAX_ITT = 40

Fixed bisection budget — no early exit (GPU-divergence-free).
real(kind=wp), private, parameter :: MIN_SPEED2 = 1.0e-8_wp

(1e-4 m/s)^2 floor: speed2_tot <= MIN_SPEED2 => cg1 = 0 (no resolvable first-baroclinic mode).
real(kind=wp), private, parameter :: RD_GUARD = 1.0e-20_wp

Inside-sqrt guard for the smooth equatorial Rd blend.
real(kind=wp), private, parameter :: RESCALE = 1024.0_wp**4

Dynamic-rescale ceiling to keep det representable for large kc.
real(kind=wp), private, parameter :: TOL_MERGE = 0.001_wp

Relative backtracking-merge threshold (MOM6 wave_speed_tol default). The criterion is relative to the column’s own stratification scale drxh_sum; no fixed absolute floor.

Derived Types

type, public ::  ocean_wave_speed_t

Components

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

|grad f| at cell centres (1/(m*s)); static, filled at configure by build_static from the same Coriolis field (the meke_length_scales centred-difference idiom). Feeds the equatorial branch of wavespeed_rd — replaces the old namelist-scalar beta.
real(kind=wp), public, allocatable :: cg1(:,:)

First-baroclinic gravity-wave speed (m/s).
logical, public :: enable = .false.

Master switch. Default off — bit-identity preserved.
real(kind=wp), public, allocatable :: f_centre(:,:)

|f| at cell centres (1/s); filled by build_static from metrics_fill_coriolis (planetary or beta-plane, per &ocean_grid_nml coriolis_scheme).
logical, public :: is_init = .false.

True between init and destroy.
real(kind=wp), public :: mono_n2_depth = -1.0_wp

depth. < 0 = off.
integer, public :: n_wavespeed = 1

diagnostic). Default every step (cheap when default-off).
real(kind=wp), public, allocatable :: rd(:,:)

First-mode Rossby deformation radius (m).
real(kind=wp), public, allocatable :: rd_over_dx(:,:)

Rd / dx (nondim) — the B2 resolution ratio. dx here is metrics%dxT (metres), NOT grid%dx (degrees on spherical/supergrid/tripolar).
real(kind=wp), public :: rho0 = 1035.0_wp

Boussinesq reference density (kg/m^3) for gprime.
logical, public :: use_ebt = .false.

DEFERRED: pressure-Neumann / equivalent-barotropic variant.

Type-Bound Procedures

procedure, public :: build_static => ocean_wave_speed_build_static
procedure, public, non_overridable :: bytes => ocean_wave_speed_bytes
procedure, public :: destroy => ocean_wave_speed_destroy
procedure, public :: enter_data => ocean_wave_speed_enter_data
procedure, public :: exit_data => ocean_wave_speed_exit_data
procedure, public :: init => ocean_wave_speed_init

Functions

public pure function wavespeed_rd(cg1, fabs, beta) result(rd)

Smooth equatorial Rd blend: Rd = cg1/sqrt(f^2 + 2betacg1). Reduces to cg1/|f| away from the equator and sqrt(cg1/(2*beta)) at f=0; a small inside-sqrt guard + denominator floor handle f = beta = 0.

Arguments

Type IntentOptional Attributes Name
real(kind=wp), intent(in) :: cg1
real(kind=wp), intent(in) :: fabs
real(kind=wp), intent(in) :: beta

Return Value real(kind=wp)

private pure function det_sign(igu, igl, kc, lam) result(sgn)

Sign of det(M(lam)) on rows 2..kc (Sturm/Hallberg recursion).

Arguments

Type IntentOptional Attributes Name
real(kind=wp), intent(in) :: igu(NZ_STACK_MAX+1)
real(kind=wp), intent(in) :: igl(NZ_STACK_MAX+1)
integer, intent(in) :: kc
real(kind=wp), intent(in) :: lam

Return Value integer

private pure function ocean_wave_speed_bytes(this) result(nbytes)

Counted allocatable footprint of the wave speed 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_wave_speed_t), intent(in) :: this

Return Value integer(kind=int64)

private pure function sturm_count(igu, igl, kc, lam) result(n_chg)

Number of Sturm-sequence sign changes (eigenvalues < lam) via the three-term determinant recursion with dynamic rescaling.

Arguments

Type IntentOptional Attributes Name
real(kind=wp), intent(in) :: igu(NZ_STACK_MAX+1)
real(kind=wp), intent(in) :: igl(NZ_STACK_MAX+1)
integer, intent(in) :: kc
real(kind=wp), intent(in) :: lam

Return Value integer


Subroutines

public pure subroutine wavespeed_cg1_column(nz, h_rak, rho_rak, rho0, cg1)

First-baroclinic wave speed for ONE column. h_rak/rho_rak are in Roundabout ordering (k=1 bed, k=nz surface), fixed-size NZ_STACK_MAX arrays; only 1..nz are read. Returns cg1 (m/s), 0 for land / homogeneous / kc<2 / sub-floor columns. Gathers+flips surface-down, backtracking convective merge, symmetric tridiag, fixed-budget Sturm-count bisection.

Arguments

Type IntentOptional Attributes Name
integer, intent(in) :: nz
real(kind=wp), intent(in) :: h_rak(NZ_STACK_MAX)
real(kind=wp), intent(in) :: rho_rak(NZ_STACK_MAX)
real(kind=wp), intent(in) :: rho0
real(kind=wp), intent(out) :: cg1

public pure subroutine wavespeed_compute(grid, metrics, this, ms)

Fill this%cg1, this%rd, this%rd_over_dx over the domain. rho_layer is a top-level allocatable that reaches the device directly, so no outer-shim tracer dereference is needed (unlike EPBL/kappa-shear). Call at the n_wavespeed cadence. Host guards + dereference here; the explicit-shape do concurrent kernel lives in wavespeed_compute_impl (outer-shim + flat-impl pattern — mirrors varmix_compute).

Arguments

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

private subroutine ocean_wave_speed_build_static(this, grid, metrics, f_centre)

Copy a pre-filled cell-centre Coriolis magnitude |f| (1/s) onto the slot (mirror of ocean_meke_set_f_centre — the caller builds f_centre via fill_coriolis_centre / metrics_fill_coriolis, which handles beta-plane AND planetary/spherical), and fill the static beta_centre = |grad f| field with the SAME centred-difference stencil meke_length_scales uses (edge rows/columns left at 0 -> the extratropical Rd = cg1/|f| branch there).

Read more…

Arguments

Type IntentOptional Attributes Name
class(ocean_wave_speed_t), intent(inout) :: this
type(hgrid_t), intent(in) :: grid
type(ocean_metrics_t), intent(in) :: metrics
real(kind=wp), intent(in) :: f_centre(grid%nx_total,grid%ny_total)

private subroutine ocean_wave_speed_destroy(this)

Arguments

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

private subroutine ocean_wave_speed_enter_data(this)

Arguments

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

private subroutine ocean_wave_speed_enter_data_impl(this)

Arguments

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

private subroutine ocean_wave_speed_exit_data(this)

Arguments

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

private subroutine ocean_wave_speed_exit_data_impl(this)

Arguments

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

private subroutine ocean_wave_speed_init(this, grid)

Allocate the persistent (nx, ny) fields. Always allocates (configure runs after init, so enable is not known yet); the off-state footprint is four 2D arrays.

Arguments

Type IntentOptional Attributes Name
class(ocean_wave_speed_t), intent(inout) :: this
type(hgrid_t), intent(in) :: grid

private pure subroutine wavespeed_compute_impl(nx, ny, nz, rho0, h_layer, rho_layer, wet_mask, f_centre, beta_centre, dxT, cg1, rd, rd_over_dx)

Flat-impl wavespeed kernel (explicit-shape; NVHPC descriptor-walk-free). Per-column Sturm-Liouville solve (wavespeed_cg1_column) + the deformation-radius blend (wavespeed_rd), then the metres-denominated resolution ratio rd_over_dx = rd / dxT.

Arguments

Type IntentOptional Attributes Name
integer, intent(in) :: nx
integer, intent(in) :: ny
integer, intent(in) :: nz
real(kind=wp), intent(in) :: rho0
real(kind=wp), intent(in) :: h_layer(nx,ny,nz)
real(kind=wp), intent(in) :: rho_layer(nx,ny,nz)
real(kind=wp), intent(in) :: wet_mask(nx,ny)
real(kind=wp), intent(in) :: f_centre(nx,ny)
real(kind=wp), intent(in) :: beta_centre(nx,ny)
real(kind=wp), intent(in) :: dxT(nx,ny)
real(kind=wp), intent(out) :: cg1(nx,ny)
real(kind=wp), intent(out) :: rd(nx,ny)
real(kind=wp), intent(out) :: rd_over_dx(nx,ny)