rdb_remap_column Module

Pure per-column conservative vertical remap, shared by both backends. Reconstruction methods: PCM — piecewise constant (donor cell) PLM — piecewise linear, minmod limiter PPM — piecewise parabolic (Colella & Woodward 1984) PPM_H4 — PPM with non-uniform 4th-order edge values (White & Adcroft 2008) PQM — piecewise quartic (White & Adcroft 2008)

All routines are pure with NZ_STACK_MAX stack workspace so they run inside do concurrent (one GPU thread per column, sequential O(nz) vertical sweep). Conservation: sum(q_newdz_new) = sum(q_olddz_old) to machine precision when sum(dz_old) = sum(dz_new).

Boundary-cell closure. Every reconstruction above PCM needs a stencil the outermost cells do not have. By default (and matching MOM6 BOUNDARY_EXTRAPOLATION = False) k=1 and k=nz collapse to PCM, so the remap is FIRST-ORDER in the two cells adjacent to the boundary no matter which method is selected — PLM, PPM, PPM_H4 and PQM all share that closure and all remap a linear-in-z profile with the same O(h) error there. The optional bnd_extrap argument selects boundary_half_jump instead: the linear-exact one-sided edge pair, which makes the whole column exact for a profile linear in z. It is the remap-side twin of rdb_ocean_pgf_reconstruct :: boundary_edges_linear. Default .false. everywhere ⇒ bit-identical.

Non-uniform-grid weights. PLM’s slope and PPM’s edge estimate are both written for a UNIFORM source column — the minmod half-difference 0.5·minmod(Δq_l, Δq_r) and the (7/12, -1/12) four-cell average are the equal-thickness specialisations of Colella & Woodward (1984) eqs (1.6)-(1.8). On a STRETCHED source column they are only first-order-consistent, so a profile linear in z is NOT reproduced — the same defect the boundary closure above fixes at k=1/k=nz, but in the interior and driven by the thickness RATIO rather than by the missing stencil. The optional nonunif argument selects the proper thickness-weighted forms (MOM6 PLM_slope_cw / CW84 (1.6)-(1.8)), which reduce to the shipped formulae exactly on a uniform column. PPM_H4 and PQM already carry thickness-weighted stencils and are unaffected except through their small-nz fallbacks. Default .false. everywhere ⇒ bit-identical.


Uses

  • module~~rdb_remap_column~~UsesGraph module~rdb_remap_column rdb_remap_column module~rdb_constants rdb_constants module~rdb_remap_column->module~rdb_constants pic_types pic_types module~rdb_constants->pic_types

Used by

  • module~~rdb_remap_column~~UsedByGraph module~rdb_remap_column rdb_remap_column module~rdb_ocean_diag_fills rdb_ocean_diag_fills module~rdb_ocean_diag_fills->module~rdb_remap_column module~rdb_ocean_state rdb_ocean_state module~rdb_ocean_diag_fills->module~rdb_ocean_state module~rdb_ocean_min_thickness rdb_ocean_min_thickness module~rdb_ocean_min_thickness->module~rdb_remap_column module~rdb_ocean_remap rdb_ocean_remap module~rdb_ocean_remap->module~rdb_remap_column module~rdb_ocean_api rdb_ocean_api module~rdb_ocean_api->module~rdb_ocean_diag_fills module~rdb_ocean_diag_derived rdb_ocean_diag_derived module~rdb_ocean_api->module~rdb_ocean_diag_derived module~rdb_ocean_dyn rdb_ocean_dyn module~rdb_ocean_api->module~rdb_ocean_dyn module~rdb_ocean_engine rdb_ocean_engine module~rdb_ocean_api->module~rdb_ocean_engine module~rdb_handle rdb_handle module~rdb_ocean_api->module~rdb_handle module~rdb_ocean_diag_derived->module~rdb_ocean_diag_fills module~rdb_ocean_diag_derived->module~rdb_ocean_state module~rdb_ocean_dyn->module~rdb_ocean_min_thickness module~rdb_ocean_dyn->module~rdb_ocean_remap module~rdb_ocean_engine->module~rdb_ocean_diag_fills module~rdb_ocean_engine->module~rdb_ocean_diag_derived module~rdb_ocean_engine->module~rdb_ocean_dyn module~rdb_ocean_setup rdb_ocean_setup module~rdb_ocean_engine->module~rdb_ocean_setup module~rdb_ocean_engine->module~rdb_ocean_state module~rdb_driver rdb_driver module~rdb_driver->module~rdb_ocean_dyn module~rdb_driver->module~rdb_ocean_engine module~rdb_driver->module~rdb_ocean_state module~rdb_handle->module~rdb_ocean_engine module~rdb_handle->module~rdb_ocean_state module~rdb_ocean_setup->module~rdb_ocean_dyn module~rdb_ocean_setup->module~rdb_ocean_state module~rdb_ocean_state->module~rdb_ocean_dyn

Variables

Type Visibility Attributes Name Initial
real(kind=wp), private, parameter :: H_MIN_FRAC = 1.0e-5_wp

Relative thickness floor in the PPM_H4 edge stencil when a consecutive thickness pair sums to ~0.
real(kind=wp), private, parameter :: H_NEGLECT = 1.0e-30_wp

Tiny thickness floor for the PPM_H4 edge singularity guard; only fires on degenerate (vanishing-layer) columns. White & Adcroft 2008.
real(kind=wp), private, parameter :: H_RATIO_FLOOR = 1.0e-12_wp

Relative thickness-ratio floor in the PQM implicit-h4 edge-value stencil (White & Adcroft 2008, roundoff-safe).
real(kind=wp), private, parameter :: PQM_MIN_FRAC = 1.0e-6_wp

Boundary-closure guard for the PQM 4th-order one-sided end fit (White & Adcroft 2008).

Functions

public pure function remap_column_preconditions_ok(nz, dz_old, dz_new, rel_tol) result(ok)

Precondition test for one remap column, as a pure predicate so the caller decides what to do about a violation (audit findings V5, V6).

Read more…

Arguments

Type IntentOptional Attributes Name
integer, intent(in) :: nz
real(kind=wp), intent(in) :: dz_old(nz)

Source-column thicknesses.
real(kind=wp), intent(in) :: dz_new(nz)

Target-column thicknesses.
real(kind=wp), intent(in) :: rel_tol

Relative tolerance on the column-total match, applied against the larger of the two totals (so a land column of zero total passes trivially).

Return Value logical

.true. when both preconditions hold.


Subroutines

public pure subroutine remap_column(method, nz, dz_old, dz_new, q_old, q_new, bnd_extrap, nonunif)

Dispatch to the requested remapping method.

Arguments

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

REMAP_PCM, REMAP_PLM, REMAP_PPM, REMAP_PPM_H4, or REMAP_PQM
integer, intent(in) :: nz
real(kind=wp), intent(in) :: dz_old(nz)
real(kind=wp), intent(in) :: dz_new(nz)
real(kind=wp), intent(in) :: q_old(nz)
real(kind=wp), intent(out) :: q_new(nz)
logical, intent(in), optional :: bnd_extrap

Boundary extrapolation (MOM6 BOUNDARY_EXTRAPOLATION). Absent or .false. (the default) ⇒ the boundary cells k=1 and k=nz reconstruct as PCM, which is first-order there. .true. ⇒ boundary_half_jump, the linear-exact one-sided closure. Ignored by PCM (no reconstruction to close).
logical, intent(in), optional :: nonunif

Non-uniform-grid reconstruction weights (&vcoord_nml remap_nonuniform_weights). Absent or .false. (the default) ⇒ PLM’s slope and PPM’s edge estimate use their equal-thickness specialisations, which are linear-exact only on a uniform SOURCE column. .true. ⇒ the Colella & Woodward (1984) (1.6)-(1.8) thickness-weighted forms, linear-exact on any source column. Inert for PCM; reaches PPM_H4/PQM only through their small-nz fallbacks, their own stencils already being thickness-weighted.

public pure subroutine remap_column_pcm(nz, dz_old, dz_new, q_old, q_new)

Piecewise-constant (donor cell) remap. Diffusive, guaranteed monotone.

Arguments

Type IntentOptional Attributes Name
integer, intent(in) :: nz
real(kind=wp), intent(in) :: dz_old(nz)

Old layer thicknesses (must sum to same total as dz_new)
real(kind=wp), intent(in) :: dz_new(nz)

New layer thicknesses
real(kind=wp), intent(in) :: q_old(nz)

Old cell-average scalar values
real(kind=wp), intent(out) :: q_new(nz)

New cell-average scalar values (conservative)

public pure subroutine remap_column_plm(nz, dz_old, dz_new, q_old, q_new, bnd_extrap, nonunif)

Piecewise-linear (minmod-limited) remap. Monotone (no new extrema). Per old layer k: q_hat(xi) = q(k) + slope(k)(2xi - 1), xi in [0,1], slope(k) = 0.5*minmod(q(k+1)-q(k), q(k)-q(k-1)). bnd_extrap (absent/.false. = default) closes the boundary cells with boundary_half_jump instead of the PCM flatten — see remap_column. nonunif (absent/.false. = default) replaces the minmod half-difference — which assumes EQUAL source thicknesses — with the thickness-weighted CW84 (1.7)/(1.8) slope; see plm_slope_nonuniform.

Arguments

Type IntentOptional Attributes Name
integer, intent(in) :: nz
real(kind=wp), intent(in) :: dz_old(nz)

Old layer thicknesses
real(kind=wp), intent(in) :: dz_new(nz)

New layer thicknesses
real(kind=wp), intent(in) :: q_old(nz)

Old cell-average scalar values
real(kind=wp), intent(out) :: q_new(nz)

New cell-average scalar values (conservative)
logical, intent(in), optional :: bnd_extrap

Linear-exact boundary-cell closure (MOM6 BOUNDARY_EXTRAPOLATION).
logical, intent(in), optional :: nonunif

Non-uniform-grid slope weights (CW84 1.7/1.8).

public pure subroutine remap_column_ppm(nz, dz_old, dz_new, q_old, q_new, bnd_extrap, nonunif)

Piecewise-parabolic (Colella & Woodward 1984) remap. Per old layer k, xi in [0,1]: q_hat(xi) = q_L + xi(q_R - q_L + q6(1 - xi)), q6 = 6q_bar - 3(q_L+q_R) Edge values: 4th-order interp + CW monotonicity limiting; boundary layers fall back to PLM-quality edges, or — under bnd_extrap — to the linear-exact one-sided pair (boundary_half_jump), which zeroes q6 there so the boundary cell carries a straight line. nonunif swaps the (7/12, -1/12) edge estimate — an EQUAL- thickness specialisation — for CW84 (1.6) on the true stencil thicknesses, and the 1|2 / (nz-1)|nz edges for the thickness-weighted two-cell value; see ppm_edge_nonuniform.

Arguments

Type IntentOptional Attributes Name
integer, intent(in) :: nz
real(kind=wp), intent(in) :: dz_old(nz)

Old layer thicknesses
real(kind=wp), intent(in) :: dz_new(nz)

New layer thicknesses
real(kind=wp), intent(in) :: q_old(nz)

Old cell-average scalar values
real(kind=wp), intent(out) :: q_new(nz)

New cell-average scalar values (conservative)
logical, intent(in), optional :: bnd_extrap

Linear-exact boundary-cell closure (MOM6 BOUNDARY_EXTRAPOLATION).
logical, intent(in), optional :: nonunif

Non-uniform-grid edge weights (CW84 1.6-1.8).

public pure subroutine remap_column_ppm_h4(nz, dz_old, dz_new, q_old, q_new, bnd_extrap, nonunif)

PPM with non-uniform 4th-order (H4) edge values (White & Adcroft 2008). As remap_column_ppm but the interior edge estimate is the thickness-weighted exactly-4th-order stencil (reduces to PPM’s (7/12,-1/12) on uniform layers), cutting spurious diapycnal mixing per remap (Ilicak et al. 2012). Limiter/reconstruction/integration are identical to PPM. H4 edge algebra is inlined (helper extraction costs 4-6% on this hot kernel). Boundary edges: outermost = PCM, second-from-boundary = non-uniform 3-cell (H3) quadratic; CW limiter clamps all edges to local monotone bounds.

Arguments

Type IntentOptional Attributes Name
integer, intent(in) :: nz
real(kind=wp), intent(in) :: dz_old(nz)

Old layer thicknesses
real(kind=wp), intent(in) :: dz_new(nz)

New layer thicknesses
real(kind=wp), intent(in) :: q_old(nz)

Old cell-average scalar values
real(kind=wp), intent(out) :: q_new(nz)

New cell-average scalar values (conservative)
logical, intent(in), optional :: bnd_extrap

Linear-exact boundary-cell closure (MOM6 BOUNDARY_EXTRAPOLATION).
logical, intent(in), optional :: nonunif

Non-uniform-grid weights. The H4/H3 stencils below are ALREADY thickness-weighted, so this only reaches the nz == 2 PLM fallback; passed through for consistency.

public pure subroutine remap_column_pqm(nz, dz_old, dz_new, q_old, q_new, bnd_extrap, nonunif)

Piecewise-quartic (PQM_IH4IH3) conservative remap (White & Adcroft 2008). Implicit-h4 edge VALUES + implicit-h3 edge SLOPES (each a diagonally-dominant tridiagonal solve with one-sided 4-cell boundary closure), per-cell quartic, W&A monotonicity limiter, conservative quartic overlap integral. Cuts diapycnal mixing per remap vs PPM/PPM_H4 (Ilicak et al. 2012). Per cell k, xi in [0,1]: q_hat(xi) = a + bxi + cxi^2 + dxi^3 + exi^4. Boundary cells reconstruct as PCM. nz < 5 falls back to REMAP_PPM (W&A boundary closure needs >= 4 cells).

Arguments

Type IntentOptional Attributes Name
integer, intent(in) :: nz
real(kind=wp), intent(in) :: dz_old(nz)

Old layer thicknesses (must sum to same total as dz_new)
real(kind=wp), intent(in) :: dz_new(nz)

New layer thicknesses
real(kind=wp), intent(in) :: q_old(nz)

Old cell-average scalar values
real(kind=wp), intent(out) :: q_new(nz)

New cell-average scalar values (conservative)
logical, intent(in), optional :: bnd_extrap

Linear-exact boundary-cell closure (MOM6 BOUNDARY_EXTRAPOLATION).
logical, intent(in), optional :: nonunif

Non-uniform-grid weights. The implicit-h4/h3 stencils below are ALREADY thickness-weighted, so this only reaches the nz < 5 PPM fallback; passed through for consistency.

private pure subroutine boundary_half_jump(h_self, h_nbr, dq_up, d)

Linear-exact half-jump across a BOUNDARY cell (k=1 or k=nz), where a centred stencil has no second neighbour.

Read more…

Arguments

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

Thickness of the boundary cell itself.
real(kind=wp), intent(in) :: h_nbr

Thickness of its single interior neighbour.
real(kind=wp), intent(in) :: dq_up

Cell-mean increment toward the surface (neighbour -> self at the surface cell, self -> neighbour at the bed cell).
real(kind=wp), intent(out) :: d

Half-jump across the boundary cell; edges are q ± d.

private pure subroutine plm_slope_nonuniform(h_l, h_c, h_r, q_l, q_c, q_r, slope)

Thickness-weighted PLM slope — Colella & Woodward (1984) eq (1.7) with the (1.8) bound, the form MOM6 ships as PLM_slope_cw.

Read more…

Arguments

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

Thickness of the cell below (toward the bed).
real(kind=wp), intent(in) :: h_c

Thickness of the cell being reconstructed.
real(kind=wp), intent(in) :: h_r

Thickness of the cell above (toward the surface).
real(kind=wp), intent(in) :: q_l

Cell mean below.
real(kind=wp), intent(in) :: q_c

Cell mean here.
real(kind=wp), intent(in) :: q_r

Cell mean above.
real(kind=wp), intent(out) :: slope

Limited half-jump across the cell.

private pure subroutine ppm_edge_nonuniform(h0, h1, h2, h3, q1, q2, dq1, dq2, edge)

Colella & Woodward (1984) eq (1.6): the fourth-order interface value between cells 1 and 2 on a NON-UNIFORM stencil h0,h1,h2,h3.

Read more…

Arguments

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

Thickness two cells below the interface.
real(kind=wp), intent(in) :: h1

Thickness of the cell just below the interface.
real(kind=wp), intent(in) :: h2

Thickness of the cell just above the interface.
real(kind=wp), intent(in) :: h3

Thickness two cells above the interface.
real(kind=wp), intent(in) :: q1

Cell mean just below the interface.
real(kind=wp), intent(in) :: q2

Cell mean just above the interface.
real(kind=wp), intent(in) :: dq1

Limited CW84 (1.7) jump of the cell below.
real(kind=wp), intent(in) :: dq2

Limited CW84 (1.7) jump of the cell above.
real(kind=wp), intent(out) :: edge

Interface value.

private pure subroutine ppm_edge_two_cell(h_l, h_r, q_l, q_r, edge)

Thickness-weighted two-cell interface value — the non-uniform generalisation of 0.5*(q_l + q_r).

Read more…

Arguments

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

Thickness of the cell below the interface.
real(kind=wp), intent(in) :: h_r

Thickness of the cell above the interface.
real(kind=wp), intent(in) :: q_l

Cell mean below.
real(kind=wp), intent(in) :: q_r

Cell mean above.
real(kind=wp), intent(out) :: edge

Interface value.

private pure subroutine ppm_jump_nonuniform(h_l, h_c, h_r, q_l, q_c, q_r, dq)

Colella & Woodward (1984) eq (1.7) — the thickness-weighted second-order jump delta a across the cell, which eq (1.6) consumes. Returned UNLIMITED, deliberately.

Read more…

Arguments

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

Thickness of the cell below.
real(kind=wp), intent(in) :: h_c

Thickness of this cell.
real(kind=wp), intent(in) :: h_r

Thickness of the cell above.
real(kind=wp), intent(in) :: q_l

Cell mean below.
real(kind=wp), intent(in) :: q_c

Cell mean here.
real(kind=wp), intent(in) :: q_r

Cell mean above.
real(kind=wp), intent(out) :: dq

Full jump across the cell.

private pure subroutine pqm_end_value_h4(dz, u, csys)

One-sided 4th-order polynomial fit of the cell averages u to the four boundary layers dz (thicknesses, must be positive), returning the four coefficients csys of the fit (White & Adcroft 2008, appendix; roundoff-safe closed form). csys(1) is the edge VALUE at the boundary interface and csys(2) is the edge SLOPE there.

Arguments

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

Thicknesses of the 4 boundary layers, starting at the edge
real(kind=wp), intent(in) :: u(4)

Cell averages of the 4 boundary layers, starting at the edge
real(kind=wp), intent(out) :: csys(4)

Coefficients of the 4th-order fit polynomial in z

private pure subroutine pqm_solve_diag_dominant(n, al, ac, au, r, x)

Diagonally-dominant tridiagonal solve; central diagonal supplied as the OFFSET ac from al + au (full pivot = ac + al + au). Never divides by zero for positive-definite ac, al, au (White & Adcroft 2008).

Arguments

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

Number of unknowns (= number of edges = nz+1)
real(kind=wp), intent(in) :: al(n)

Lower diagonal (al(1) unused)
real(kind=wp), intent(in) :: ac(n)

Central-diagonal OFFSET from al+au (full diagonal = ac+al+au)
real(kind=wp), intent(in) :: au(n)

Upper diagonal (au(n) unused)
real(kind=wp), intent(in) :: r(n)

Right-hand side
real(kind=wp), intent(out) :: x(n)

Solution vector