rdb_ocean_isopycnal_slopes Module

Neutral-density slope S = -∇ρ/∂_zρ and interface stratification N² at C-grid layer interfaces (u-faces → slope_x/n2_u, v-faces → slope_y/n2_v). Purely diagnostic — consumed by GM/Redi/VarMix/ MLE; no flux consumer here. Harmonic-thickness-weighted FD form (Griffies 1998).

Bottom-up convention (k=1 bed, k=nz_ml surface). Interface K 1..nz_ml+1; K=1 bed and K=nz_ml+1 surface forced to zero slope. Interior K straddles layer k=K (above, surface side) and k=K-1 (below, bed side). A vert-fill pre-pass diffuses T/S into massless layers so T=S=0 ghosts don’t corrupt gradients.

Geopotential interface heights (the bed datum)

The along-layer density gradient is rotated to the horizontal by the interface-tilt term −∂zρ·(e_W − e_E), so e_int must be the TRUE geopotential height of each interface: it is built bed-up from e_int(:,:,1) = −D, with D the slot’s own copy of the bathymetry (barotropic%b, m, positive down below the z = 0 datum — the same datum the FV-MOM6 PGF builds e_face from). The column top is then Σh − D, i.e. η in open ocean and −z_draft + η under an ice shelf (Σh = bt_H_ref + η, bt_H_ref = b − z_draft), and a horizontally uniform stratification over ANY bathymetry reads zero slope on every coordinate. (The pre-fix zero bed datum differenced heights above the LOCAL bed and read a bathymetry step as an isopycnal slope ~ΔD/Δx.) set_bathymetry fills the copy at setup from the wrapped + halo-exchanged b (ghost-correct at periodic and MPI seams), before enter_data; an enabled ocean_slopes_compute fails loud if it was never set.

Partial-step z-level faces (&vcoord_nml zfixed_closed_faces)

The z_fixed target hangs every nominal interface at a fixed depth below z = 0, so with geopotential e_int two columns’ common interior interface differs in height by O(η_W − η_E) only and the tilt term is the (small, correct) free-surface tilt of the coordinate — the general formula, no special case. One addition, host-gated on metrics%use_closed_faces (configure admits it under z_fixed only), knob off ⇒ not taken:

  • Open-column mask. slope / N² are zeroed at every interface that is not strictly inside the face’s open column (open(ka) .and. open(kb)), so no consumer (GM, its gm_src, VarMix) reads a gradient formed against a filler.

Default off (&ocean_slopes_nml enable=.false.) ⇒ slot allocated but ocean_slopes_compute no-ops ⇒ bit-identical.

exposed for the interface-pressure unit test


Uses

  • module~~rdb_ocean_isopycnal_slopes~~UsesGraph module~rdb_ocean_isopycnal_slopes rdb_ocean_isopycnal_slopes iso_fortran_env iso_fortran_env module~rdb_ocean_isopycnal_slopes->iso_fortran_env module~rdb_constants rdb_constants module~rdb_ocean_isopycnal_slopes->module~rdb_constants module~rdb_eos rdb_eos module~rdb_ocean_isopycnal_slopes->module~rdb_eos module~rdb_grid rdb_grid module~rdb_ocean_isopycnal_slopes->module~rdb_grid module~rdb_mem_report rdb_mem_report module~rdb_ocean_isopycnal_slopes->module~rdb_mem_report module~rdb_multilayer_state rdb_multilayer_state module~rdb_ocean_isopycnal_slopes->module~rdb_multilayer_state module~rdb_ocean_metrics rdb_ocean_metrics module~rdb_ocean_isopycnal_slopes->module~rdb_ocean_metrics 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_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_isopycnal_slopes~~UsedByGraph module~rdb_ocean_isopycnal_slopes rdb_ocean_isopycnal_slopes module~rdb_ocean_dyn rdb_ocean_dyn module~rdb_ocean_dyn->module~rdb_ocean_isopycnal_slopes module~rdb_ocean_gm rdb_ocean_gm module~rdb_ocean_dyn->module~rdb_ocean_gm module~rdb_ocean_varmix rdb_ocean_varmix module~rdb_ocean_dyn->module~rdb_ocean_varmix module~rdb_continuity rdb_continuity module~rdb_ocean_dyn->module~rdb_continuity module~rdb_ocean_meke rdb_ocean_meke module~rdb_ocean_dyn->module~rdb_ocean_meke module~rdb_ocean_gm->module~rdb_ocean_isopycnal_slopes module~rdb_ocean_state rdb_ocean_state module~rdb_ocean_state->module~rdb_ocean_isopycnal_slopes module~rdb_ocean_state->module~rdb_ocean_dyn module~rdb_ocean_state->module~rdb_ocean_gm module~rdb_ocean_state->module~rdb_ocean_varmix module~rdb_ocean_state->module~rdb_continuity module~rdb_ocean_state->module~rdb_ocean_meke module~rdb_ocean_varmix->module~rdb_ocean_isopycnal_slopes module~rdb_continuity->module~rdb_ocean_gm 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_ice_transport rdb_ice_transport module~rdb_ocean_engine->module~rdb_ice_transport module~rdb_ocean_meke->module~rdb_ocean_gm module~rdb_ocean_meke->module~rdb_ocean_varmix module~rdb_ocean_setup->module~rdb_ocean_dyn module~rdb_ocean_setup->module~rdb_ocean_state module~rdb_ice_transport->module~rdb_continuity

Derived Types

type, public ::  ocean_slopes_t

Components

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

Bed depth D below the z = 0 datum (m, positive down), (nx, ny) INCLUDING ghosts — the slot’s own copy of barotropic%b, taken by set_bathymetry after the periodic wrap + halo exchange and before enter_data. The bed datum of e_int; see the module docstring.
logical, public :: bathy_set = .false.

True once set_bathymetry has filled bathy. An enabled ocean_slopes_compute fails loud without it: a silent zero datum is exactly the bathymetry-as-slope defect.
real(kind=wp), public, allocatable :: e_int(:,:,:)

Geopotential interface height (m, positive up from z = 0), (nx, ny, nz+1): e_int(:,:,1) = −bathy (bed), then bottom-up cumulative + h_layer; (:,:,nz+1) = column top.
logical, public :: enable = .false.

Master switch (&ocean_slopes_nml enable). Default off ⇒ ocean_slopes_compute no-ops ⇒ bit-identical.
logical, public :: is_init = .false.

True between init and destroy. Guard on this, never on allocated(...) (host pointer only; misses GPU mapping).
real(kind=wp), public :: kd_smooth = 1.0e-6_wp

Vertical diffusivity (m²/s) used by vert_fill_TS to fill massless layers. Multiplied by dt for the smoothing kappa·dt.
real(kind=wp), public :: min_dz_for_n2 = 1.0_wp

Minimum layer thickness (m) used to floor h in the N² vertical-difference denominator, so vanished layers don’t produce a spurious N² spike.
real(kind=wp), public, allocatable :: n2_u(:,:,:)

Brunt-Väisälä N² at u-faces (s⁻²), shape (nx+1, ny, nz+1).
real(kind=wp), public, allocatable :: n2_v(:,:,:)

Brunt-Väisälä N² at v-faces (s⁻²), shape (nx, ny+1, nz+1).
integer, public :: nx_total = 0
integer, public :: ny_total = 0
integer, public :: nz_ml = 0
real(kind=wp), public :: rho0 = 1035.0_wp

Reference density (kg/m³) for the N² scaling g/ρ₀.
real(kind=wp), public, allocatable :: s_fill(:,:,:)

Massless-layer-filled salinity scratch, (nx, ny, nz).
real(kind=wp), public, allocatable :: slope_x(:,:,:)

Isopycnal slope at u-faces, shape (nx+1, ny, nz+1).
real(kind=wp), public, allocatable :: slope_y(:,:,:)

Isopycnal slope at v-faces, shape (nx, ny+1, nz+1).
real(kind=wp), public, allocatable :: t_fill(:,:,:)

Massless-layer-filled temperature scratch, (nx, ny, nz).

Type-Bound Procedures

procedure, public, non_overridable :: bytes => ocean_slopes_bytes
procedure, public, non_overridable :: destroy => ocean_slopes_destroy
procedure, public, non_overridable :: enter_data => ocean_slopes_enter_data
procedure, public, non_overridable :: exit_data => ocean_slopes_exit_data
procedure, public, non_overridable :: init => ocean_slopes_init
procedure, public, non_overridable :: set_bathymetry => ocean_slopes_set_bathymetry

Functions

public pure function pressure_above_x(nx, ny, nz, h_layer, ic, jc, ka, rho0) result(p)

Surface-relative hydrostatic pressure at the interface K straddled by layer ka (above, surface-side) and ka-1 (below): the interface sits at the BOTTOM of layer ka, so the water column above it is layers ka..nz (bottom-up, k=nz the surface). p = g·ρ₀·Σ_{k’=ka}^{nz} h(k’) — the sum INCLUDES ka (the layer directly above the interface); omitting it shorts the pressure by one layer (~5e5 Pa) and biases pressure-dependent EOS derivatives.

Arguments

Type IntentOptional Attributes Name
integer, intent(in) :: nx
integer, intent(in) :: ny
integer, intent(in) :: nz
real(kind=wp), intent(in) :: h_layer(nx,ny,nz)
integer, intent(in) :: ic
integer, intent(in) :: jc
integer, intent(in) :: ka
real(kind=wp), intent(in) :: rho0

Return Value real(kind=wp)

private pure function ocean_slopes_bytes(this) result(nbytes)

Counted allocatable footprint of the isopycnal slopes 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_slopes_t), intent(in) :: this

Return Value integer(kind=int64)


Subroutines

public subroutine ocean_slopes_compute(grid, metrics, eos, slopes, ms, dt)

Public entry point — fill slope_x/slope_y + n2_u/n2_v at all interfaces. No-op if absent / uninitialised / disabled, so the driver can call it unconditionally. Pipeline: vert-fill T/S → build interface heights e_int → u-face pass → v-face pass.

Arguments

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

public subroutine ocean_slopes_vert_fill_ts(nx, ny, nz, h_layer, t_htr, s_htr, kd_smooth, dt, t_fill, s_fill)

Fill massless layers in T/S with sensible values via one pass of constant-kappa·dt vertical diffusion — a SINGLE forward-elim + back-sub Thomas sweep per column (no iteration). Operates on the tracer-from-hTr conversion (T = hTr/h) and writes the scratch t_fill/s_fill; the prognostic tracers are untouched.

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_layer(nx,ny,nz)
real(kind=wp), intent(in) :: t_htr(nx,ny,nz)
real(kind=wp), intent(in) :: s_htr(nx,ny,nz)
real(kind=wp), intent(in) :: kd_smooth
real(kind=wp), intent(in) :: dt
real(kind=wp), intent(out) :: t_fill(nx,ny,nz)
real(kind=wp), intent(out) :: s_fill(nx,ny,nz)

private pure subroutine ocean_slopes_build_e(nx, ny, nz, bathy, h_layer, e_int)

Build GEOPOTENTIAL interface heights bottom-up: e_int(:,:,1) = −bathy (the bed, below the z = 0 datum), e_int(:,:,K+1) = e_int(:,:,K) + h_layer(:,:,K). A per-column serial cumulative sum (parallel over i,j). The across-face difference e_W − e_E feeds the interface-tilt term, so the bed datum is NOT irrelevant: it must be the true bed depth, or a bathymetry step reads as an isopycnal slope.

Arguments

Type IntentOptional Attributes Name
integer, intent(in) :: nx
integer, intent(in) :: ny
integer, intent(in) :: nz
real(kind=wp), intent(in) :: bathy(nx,ny)
real(kind=wp), intent(in) :: h_layer(nx,ny,nz)
real(kind=wp), intent(out) :: e_int(nx,ny,nz+1)

private subroutine ocean_slopes_compute_impl(grid, metrics, eos, slopes, ms, h_layer, t_htr, s_htr, dt)

NVHPC doesn’t descriptor-walk per launch.

Arguments

Type IntentOptional Attributes Name
type(hgrid_t), intent(in) :: grid
type(ocean_metrics_t), intent(in) :: metrics
type(eos_t), intent(in) :: eos
type(ocean_slopes_t), intent(inout) :: slopes
type(multilayer_state_t), intent(in) :: ms
real(kind=wp), intent(in) :: h_layer(slopes%nx_total,slopes%ny_total,slopes%nz_ml)
real(kind=wp), intent(in) :: t_htr(slopes%nx_total,slopes%ny_total,slopes%nz_ml)
real(kind=wp), intent(in) :: s_htr(slopes%nx_total,slopes%ny_total,slopes%nz_ml)
real(kind=wp), intent(in) :: dt

private subroutine ocean_slopes_destroy(this)

Arguments

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

private subroutine ocean_slopes_enter_data(this)

Poly TBP delegating to a type(...)-arg _impl (AMD-crash rule: bare polymorphic copyin(this) maps the stack descriptor → AMD libomptarget cross-slot overlap crash).

Arguments

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

private subroutine ocean_slopes_enter_data_impl(this)

Arguments

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

private subroutine ocean_slopes_exit_data(this)

Arguments

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

private subroutine ocean_slopes_exit_data_impl(this)

Arguments

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

private subroutine ocean_slopes_init(this, grid, nz_ml)

Allocate the slope / N² outputs + the vert-fill T/S scratch + the interface-height buffer. Default nz_ml = 1 preserves the barotropic-only constructor; pass nz_ml = ms%nz_ml for the multilayer driver. Setup code uses plain host allocation (no do concurrent before enter_data).

Arguments

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

private pure subroutine ocean_slopes_mask_open_column(nx, ny, nz, open_u, open_v, slope_x, slope_y, n2_u, n2_v)

Zero slope / N² at every interior interface K that is NOT strictly inside its face’s open column, i.e. unless both layers it separates (K above, K-1 below) are open at that face (&vcoord_nml zfixed_closed_faces). Assigned under a test, never multiplied by the 0/1 mask, so a non-finite value formed against a filler cannot survive as NaN·0.

Arguments

Type IntentOptional Attributes Name
integer, intent(in) :: nx
integer, intent(in) :: ny
integer, intent(in) :: nz
real(kind=wp), intent(in) :: open_u(nx+1,ny,nz)
real(kind=wp), intent(in) :: open_v(nx,ny+1,nz)
real(kind=wp), intent(inout) :: slope_x(nx+1,ny,nz+1)
real(kind=wp), intent(inout) :: slope_y(nx,ny+1,nz+1)
real(kind=wp), intent(inout) :: n2_u(nx+1,ny,nz+1)
real(kind=wp), intent(inout) :: n2_v(nx,ny+1,nz+1)

private pure subroutine ocean_slopes_pass_x(nx, ny, nz, eos, rho0, min_dz, h_layer, t_fill, s_fill, e_int, idxCu, wet_u, slope_x, n2_u)

u-face slope + N² pass. Interface K (interior 2..nz) straddles layer k=K (above, surface side) and k=K-1 (below, bed side). Bed (K=1) + surface (K=nz+1) are forced to zero. The u-face at (i,j) sits between cells (i-1,j) and (i,j); pairs columns iw=i-1 (west) and i (east), so loop i=2:nx.

Arguments

Type IntentOptional Attributes Name
integer, intent(in) :: nx
integer, intent(in) :: ny
integer, intent(in) :: nz
type(eos_t), intent(in) :: eos
real(kind=wp), intent(in) :: rho0
real(kind=wp), intent(in) :: min_dz
real(kind=wp), intent(in) :: h_layer(nx,ny,nz)
real(kind=wp), intent(in) :: t_fill(nx,ny,nz)
real(kind=wp), intent(in) :: s_fill(nx,ny,nz)
real(kind=wp), intent(in) :: e_int(nx,ny,nz+1)
real(kind=wp), intent(in) :: idxCu(nx+1,ny)
real(kind=wp), intent(in) :: wet_u(nx+1,ny)
real(kind=wp), intent(out) :: slope_x(nx+1,ny,nz+1)
real(kind=wp), intent(out) :: n2_u(nx+1,ny,nz+1)

private pure subroutine ocean_slopes_pass_y(nx, ny, nz, eos, rho0, min_dz, h_layer, t_fill, s_fill, e_int, idyCv, wet_v, slope_y, n2_v)

v-face slope + N² pass — mirror of pass_x with v-staggering. The v-face at (i,j) sits between cells (i,j-1) and (i,j); pairs columns js=j-1 (south) and j (north), loop j=2:ny.

Arguments

Type IntentOptional Attributes Name
integer, intent(in) :: nx
integer, intent(in) :: ny
integer, intent(in) :: nz
type(eos_t), intent(in) :: eos
real(kind=wp), intent(in) :: rho0
real(kind=wp), intent(in) :: min_dz
real(kind=wp), intent(in) :: h_layer(nx,ny,nz)
real(kind=wp), intent(in) :: t_fill(nx,ny,nz)
real(kind=wp), intent(in) :: s_fill(nx,ny,nz)
real(kind=wp), intent(in) :: e_int(nx,ny,nz+1)
real(kind=wp), intent(in) :: idyCv(nx,ny+1)
real(kind=wp), intent(in) :: wet_v(nx,ny+1)
real(kind=wp), intent(out) :: slope_y(nx,ny+1,nz+1)
real(kind=wp), intent(out) :: n2_v(nx,ny+1,nz+1)

private subroutine ocean_slopes_set_bathymetry(this, b)

Copy the bed depth b (m, positive down, (nx_total, ny_total) incl. ghosts) into this%bathy on the host. Call it with the WRAPPED + halo-exchanged barotropic%b (the seam faces read the ghosts) and BEFORE enter_data (the device copy is taken from the host values); to refresh after enter_data the caller issues !$acc update device(this%bathy) itself. No-op on an uninitialised slot.

Arguments

Type IntentOptional Attributes Name
class(ocean_slopes_t), intent(inout) :: this
real(kind=wp), intent(in) :: b(:,:)