mom_neutral_diffusion module reference

A column-wise toolbox for implementing neutral diffusion.

More…

Data Types

neutral_diffusion_cs

The control structure for the MOM_neutral_diffusion module.

Functions/Subroutines

neutral_diffusion_init()

Read parameters and allocate control structure for neutral_diffusion module.

neutral_diffusion_calc_coeffs()

Calculate remapping factors for u/v columns used to map adjoining columns to a shared coordinate space.

neutral_diffusion()

Update tracer concentration due to neutral diffusion; layer thickness unchanged by this update.

interface_scalar()

Returns interface scalar, Si, for a column of layer values, S.

ppm_edge()

Returns the PPM quasi-fourth order edge value at k+1/2 following equation 1.6 in Colella & Woodward, 1984: JCP 54, 174-201.

ppm_ave()

Returns the average of a PPM reconstruction between two fractional positions.

signum()

A true signum function that returns either -abs(a), when x<0; or abs(a) when x>0; or 0 when x=0.

plm_diff()

Returns PLM slopes for a column where the slopes are the difference in value across each cell.

fv_diff()

Returns the cell-centered second-order finite volume (unlimited PLM) slope using three consecutive cell widths and average values.

fvlsq_slope()

Returns the cell-centered second-order weighted least squares slope using three consecutive cell widths and average values.

find_neutral_surface_positions_continuous()

Returns positions within left/right columns of combined interfaces using continuous reconstructions of T/S.

interpolate_for_nondim_position()

Returns the non-dimensional position between Pneg and Ppos where the interpolated density difference equals zero.

find_neutral_surface_positions_discontinuous()

Higher order version of find_neutral_surface_positions.

mark_unstable_cells()

Sweep down through the column and mark as stable if the bottom interface of a cell is denser than the top.

search_other_column()

Searches the “other” (searched) column for the position of the neutral surface.

increment_interface()

Increments the interface which was just connected and also set flags if the bottom is reached.

find_neutral_pos_linear()

Search a layer to find where delta_rho = 0 based on a linear interpolation of alpha and beta of the top and bottom being searched and polynomial reconstructions of T and S.

find_neutral_pos_full()

Use the full equation of state to calculate the difference in locally referenced potential density.

calc_delta_rho_and_derivs()

Calculate the difference in density between two points in a variety of ways.

delta_rho_from_derivs()

Calculate delta rho from derivatives and gradients of properties \(\Delta \rho = \frac{1}{2}\left[ (\alpha_1 + \alpha_2)*(T_1-T_2) + (\beta_1 + \beta_2)*(S_1-S_2) + (\gamma^{-1}_1 + \gamma^{-1}_2)*(P_1-P_2) \right]\).

absolute_position()

Converts non-dimensional position within a layer to absolute position (for debugging)

absolute_positions()

Converts non-dimensional positions within layers to absolute positions (for debugging)

neutral_surface_flux()

Returns a single column of neutral diffusion fluxes of a tracer.

neutral_surface_t_eval()

Evaluate various parts of the reconstructions to calculate gradient-based flux limter.

ppm_left_right_edge_values()

Discontinuous PPM reconstructions of the left/right edge values within a cell.

neutral_diffusion_unit_tests()

Returns true if unit tests of neutral_diffusion functions fail.

ndiff_unit_tests_continuous()

Returns true if unit tests of neutral_diffusion functions fail.

ndiff_unit_tests_discontinuous()

test_fv_diff()

Returns true if a test of fv_diff() fails, and conditionally writes results to stream.

test_fvlsq_slope()

Returns true if a test of fvlsq_slope() fails, and conditionally writes results to stream.

test_ifndp()

Returns true if a test of interpolate_for_nondim_position() fails, and conditionally writes results to stream.

test_data1d()

Returns true if comparison of Po and Ptrue fails, and conditionally writes results to stream.

test_data1di()

Returns true if comparison of Po and Ptrue fails, and conditionally writes results to stream.

test_nsp()

Returns true if output of find_neutral_surface_positions() does not match correct values, and conditionally writes results to stream.

compare_nsp_row()

Compares a single row, k, of output from find_neutral_surface_positions()

test_rnp()

Compares output position from refine_nondim_position with an expected value.

neutral_diffusion_end()

Deallocates neutral_diffusion control structure.

Detailed Description

A column-wise toolbox for implementing neutral diffusion.

Type Documentation

type mom_neutral_diffusion/neutral_diffusion_cs

The control structure for the MOM_neutral_diffusion module.

Type fields
  • % nkp1 [integer] :: Number of interfaces for a column = nk + 1.

  • % nsurf [integer] :: Number of neutral surfaces.

  • % deg [integer] :: Degree of polynomial used for reconstructions.

  • % continuous_reconstruction [logical] :: True if using continuous PPM reconstruction at interfaces.

  • % debug [logical] :: If true, write verbose debugging messages.

  • % hard_fail_heff [logical] :: Bring down the model if a problem with heff is detected.

  • % max_iter [integer] :: Maximum number of iterations if refine_position is defined.

  • % drho_tol [real] :: Convergence criterion representing density difference from true neutrality [R ~> kg m-3].

  • % x_tol [real] :: Convergence criterion for how small an update of the position can be.

  • % ref_pres [real] :: Reference pressure, negative if using locally referenced neutral density [R L2 T-2 ~> Pa].

  • % interior_only [logical] :: If true, only applies neutral diffusion in the ocean interior. That is, the algorithm will exclude the surface and bottom boundary layers.

  • % use_unmasked_transport_bug [logical] :: If true, use an older form for the accumulation of neutral-diffusion transports that were unmasked, as used prior to Jan 2018.

  • % upol [real(:,:,:),allocatable] :: Non-dimensional position with left layer uKoL-1, u-point.

  • % upor [real(:,:,:),allocatable] :: Non-dimensional position with right layer uKoR-1, u-point.

  • % ukol [integer(:,:,:),allocatable] :: Index of left interface corresponding to neutral surface, at a u-point.

  • % ukor [integer(:,:,:),allocatable] :: Index of right interface corresponding to neutral surface, at a u-point.

  • % uheff [real(:,:,:),allocatable] :: Effective thickness at u-point [H ~> m or kg m-2].

  • % vpol [real(:,:,:),allocatable] :: Non-dimensional position with left layer uKoL-1, v-point.

  • % vpor [real(:,:,:),allocatable] :: Non-dimensional position with right layer uKoR-1, v-point.

  • % vkol [integer(:,:,:),allocatable] :: Index of left interface corresponding to neutral surface, at a v-point.

  • % vkor [integer(:,:,:),allocatable] :: Index of right interface corresponding to neutral surface, at a v-point.

  • % vheff [real(:,:,:),allocatable] :: Effective thickness at v-point [H ~> m or kg m-2].

  • % ppoly_coeffs_t [real(:,:,:,:),allocatable] :: Polynomial coefficients for temperature.

  • % ppoly_coeffs_s [real(:,:,:,:),allocatable] :: Polynomial coefficients for salinity.

  • % drdt [real(:,:,:),allocatable] :: dRho/dT [R degC-1 ~> kg m-3 degC-1] at interfaces

  • % drds [real(:,:,:),allocatable] :: dRho/dS [R ppt-1 ~> kg m-3 ppt-1] at interfaces

  • % tint [real(:,:,:),allocatable] :: Interface T [degC].

  • % sint [real(:,:,:),allocatable] :: Interface S [ppt].

  • % pint [real(:,:,:),allocatable] :: Interface pressure [R L2 T-2 ~> Pa].

  • % t_i [real(:,:,:,:),allocatable] :: Top edge reconstruction of temperature [degC].

  • % s_i [real(:,:,:,:),allocatable] :: Top edge reconstruction of salinity [ppt].

  • % p_i [real(:,:,:,:),allocatable] :: Interface pressures [R L2 T-2 ~> Pa].

  • % drdt_i [real(:,:,:,:),allocatable] :: dRho/dT [R degC-1 ~> kg m-3 degC-1] at top edge

  • % drds_i [real(:,:,:,:),allocatable] :: dRho/dS [R ppt-1 ~> kg m-3 ppt-1] at top edge

  • % ns [integer(:,:),allocatable] :: Number of interfacs in a column.

  • % stable_cell [logical(:,:,:),allocatable] :: True if the cell is stably stratified wrt to the next cell.

  • % r_to_kg_m3 [real] :: A rescaling factor translating density to kg m-3 for use in diagnostic messages [kg m-3 R-1 ~> 1].

  • % diag [type(diag_ctrl),pointer] :: A structure that is used to regulate the timing of diagnostic output.

  • % neutral_pos_method [integer] :: Method to find the position of a neutral surface within the layer.

  • % delta_rho_form [character (len=40)] :: Determine which (if any) approximation is made to the equation describing the difference in density.

  • % id_uheff_2d [integer] :: Diagnostic IDs.

  • % id_vheff_2d [integer] :: Diagnostic IDs.

  • % eos [type(eos_type),pointer] :: Equation of state parameters.

  • % remap_cs [type(remapping_cs)] :: Remapping control structure used to create sublayers.

  • % remap_answers_2018 [logical] :: If true, use the order of arithmetic and expressions that recover the answers for remapping from the end of 2018. Otherwise, use more robust forms of the same expressions.

  • % kpp_csp [type(kpp_cs),pointer] :: KPP control structure needed to get BLD.

  • % energetic_pbl_csp [type(energetic_pbl_cs),pointer] :: ePBL control structure needed to get MLD

Function/Subroutine Documentation

function mom_neutral_diffusion/neutral_diffusion_init(Time, G, GV, US, param_file, diag, EOS, diabatic_CSp, CS) [logical]

Read parameters and allocate control structure for neutral_diffusion module.

Parameters
  • time :: [in] Time structure

  • g :: [in] Grid structure

  • gv :: [in] The ocean’s vertical grid structure

  • us :: [in] A dimensional unit scaling type

  • diag :: [inout] Diagnostics control structure

  • param_file :: [in] Parameter file structure

  • eos :: [in] Equation of state

  • diabatic_csp :: KPP control structure needed to get BLD

  • cs :: Neutral diffusion control structure

Call to

mom_diabatic_driver::extract_diabatic_member mdl mom_error_handler::mom_error mom_remapping::remappingdefaultscheme mom_remapping::remappingschemesdoc

subroutine mom_neutral_diffusion/neutral_diffusion_calc_coeffs(G, GV, US, h, T, S, CS, p_surf)

Calculate remapping factors for u/v columns used to map adjoining columns to a shared coordinate space.

Parameters
  • g :: [in] Ocean grid structure

  • gv :: [in] ocean vertical grid structure

  • us :: [in] A dimensional unit scaling type

  • h :: [in] Layer thickness [H ~> m or kg m-2]

  • t :: [in] Potential temperature [degC]

  • s :: [in] Salinity [ppt]

  • cs :: Neutral diffusion control structure

  • p_surf :: [in] Surface pressure to include in pressures used for equation of state calculations [R L2 T-2 ~> Pa]

Call to

mom_lateral_boundary_diffusion::boundary_k_range mom_energetic_pbl::energetic_pbl_get_mld polynomial_functions::evaluation_polynomial find_neutral_surface_positions_continuous find_neutral_surface_positions_discontinuous interface_scalar mom_cvmix_kpp::kpp_get_bld mark_unstable_cells mom_lateral_boundary_diffusion::surface

Called from

mom_tracer_hor_diff::tracer_hordiff

subroutine mom_neutral_diffusion/neutral_diffusion(G, GV, h, Coef_x, Coef_y, dt, Reg, US, CS)

Update tracer concentration due to neutral diffusion; layer thickness unchanged by this update.

Parameters
  • g :: [in] Ocean grid structure

  • gv :: [in] ocean vertical grid structure

  • h :: [in] Layer thickness [H ~> m or kg m-2]

  • coef_x :: [in] dt * Kh * dy / dx at u-points [L2 ~> m2]

  • coef_y :: [in] dt * Kh * dx / dy at v-points [L2 ~> m2]

  • dt :: [in] Tracer time step * I_numitts [T ~> s] (I_numitts in tracer_hordiff)

  • reg :: Tracer registry

  • us :: [in] A dimensional unit scaling type

  • cs :: Neutral diffusion control structure

Call to

neutral_surface_flux

Called from

mom_tracer_hor_diff::tracer_hordiff

subroutine mom_neutral_diffusion/interface_scalar(nk, h, S, Si, i_method, h_neglect)

Returns interface scalar, Si, for a column of layer values, S.

Parameters
  • nk :: [in] Number of levels

  • h :: [in] Layer thickness [H ~> m or kg m-2]

  • s :: [in] Layer scalar (conc, e.g. ppt)

  • si :: [inout] Interface scalar (conc, e.g. ppt)

  • i_method :: [in] =1 use average of PLM edges =2 use continuous PPM edge interpolation

  • h_neglect :: [in] A negligibly small thickness [H ~> m or kg m-2]

Call to

plm_diff ppm_edge

Called from

ndiff_unit_tests_continuous neutral_diffusion_calc_coeffs neutral_surface_flux

function mom_neutral_diffusion/ppm_edge(hkm1, hk, hkp1, hkp2, Ak, Akp1, Pk, Pkp1, h_neglect) [real]

Returns the PPM quasi-fourth order edge value at k+1/2 following equation 1.6 in Colella & Woodward, 1984: JCP 54, 174-201.

Parameters
  • hkm1 :: [in] Width of cell k-1

  • hk :: [in] Width of cell k

  • hkp1 :: [in] Width of cell k+1

  • hkp2 :: [in] Width of cell k+2

  • ak :: [in] Average scalar value of cell k

  • akp1 :: [in] Average scalar value of cell k+1

  • pk :: [in] PLM slope for cell k

  • pkp1 :: [in] PLM slope for cell k+1

  • h_neglect :: [in] A negligibly small thickness [H ~> m or kg m-2]

Called from

interface_scalar

function mom_neutral_diffusion/ppm_ave(xL, xR, aL, aR, aMean) [real]

Returns the average of a PPM reconstruction between two fractional positions.

Parameters
  • xl :: [in] Fraction position of left bound (0,1)

  • xr :: [in] Fraction position of right bound (0,1)

  • al :: [in] Left edge scalar value, at x=0

  • ar :: [in] Right edge scalar value, at x=1

  • amean :: [in] Average scalar value of cell

Called from

neutral_surface_flux

function mom_neutral_diffusion/signum(a, x) [real]

A true signum function that returns either -abs(a), when x<0; or abs(a) when x>0; or 0 when x=0.

Parameters
  • a :: [in] The magnitude argument

  • x :: [in] The sign (or zero) argument

Called from

neutral_surface_flux plm_diff ppm_left_right_edge_values

subroutine mom_neutral_diffusion/plm_diff(nk, h, S, c_method, b_method, diff)

Returns PLM slopes for a column where the slopes are the difference in value across each cell. The limiting follows equation 1.8 in Colella & Woodward, 1984: JCP 54, 174-201.

Parameters
  • nk :: [in] Number of levels

  • h :: [in] Layer thickness [H ~> m or kg m-2]

  • s :: [in] Layer salinity (conc, e.g. ppt)

  • c_method :: [in] Method to use for the centered difference

  • b_method :: [in] =1, use PCM in first/last cell, =2 uses linear extrapolation

  • diff :: [inout] Scalar difference across layer (conc, e.g. ppt) determined by the following values for c_method:

System Message: WARNING/2 (docstring of plm_diff, line 12)

Field list ends without a blank line; unexpected unindent.

  1. Second order finite difference (not recommended)

  2. Second order finite volume (used in original PPM)

  3. Finite-volume weighted least squares linear fit

Todo

The use of c_method to choose a scheme is inefficient and should eventually be moved up the call tree.

Call to

fv_diff fvlsq_slope signum

Called from

interface_scalar

function mom_neutral_diffusion/fv_diff(hkm1, hk, hkp1, Skm1, Sk, Skp1) [real]

Returns the cell-centered second-order finite volume (unlimited PLM) slope using three consecutive cell widths and average values. Slope is returned as a difference across the central cell (i.e. units of scalar S). Discretization follows equation 1.7 in Colella & Woodward, 1984: JCP 54, 174-201.

Parameters
  • hkm1 :: [in] Left cell width

  • hk :: [in] Center cell width

  • hkp1 :: [in] Right cell width

  • skm1 :: [in] Left cell average value

  • sk :: [in] Center cell average value

  • skp1 :: [in] Right cell average value

Called from

plm_diff test_fv_diff

function mom_neutral_diffusion/fvlsq_slope(hkm1, hk, hkp1, Skm1, Sk, Skp1) [real]

Returns the cell-centered second-order weighted least squares slope using three consecutive cell widths and average values. Slope is returned as a gradient (i.e. units of scalar S over width units).

Parameters
  • hkm1 :: [in] Left cell width

  • hk :: [in] Center cell width

  • hkp1 :: [in] Right cell width

  • skm1 :: [in] Left cell average value

  • sk :: [in] Center cell average value

  • skp1 :: [in] Right cell average value

Called from

plm_diff test_fvlsq_slope

subroutine mom_neutral_diffusion/find_neutral_surface_positions_continuous(nk, Pl, Tl, Sl, dRdTl, dRdSl, Pr, Tr, Sr, dRdTr, dRdSr, PoL, PoR, KoL, KoR, hEff, bl_kl, bl_kr, bl_zl, bl_zr)

Returns positions within left/right columns of combined interfaces using continuous reconstructions of T/S.

Parameters
  • nk :: [in] Number of levels

  • pl :: [in] Left-column interface pressure [R L2 T-2 ~> Pa] or other units

  • tl :: [in] Left-column interface potential temperature [degC]

  • sl :: [in] Left-column interface salinity [ppt]

  • drdtl :: [in] Left-column dRho/dT [R degC-1 ~> kg m-3 degC-1]

  • drdsl :: [in] Left-column dRho/dS [R ppt-1 ~> kg m-3 ppt-1]

  • pr :: [in] Right-column interface pressure [R L2 T-2 ~> Pa] or other units

  • tr :: [in] Right-column interface potential temperature [degC]

  • sr :: [in] Right-column interface salinity [ppt]

  • drdtr :: [in] Left-column dRho/dT [R degC-1 ~> kg m-3 degC-1]

  • drdsr :: [in] Left-column dRho/dS [R ppt-1 ~> kg m-3 ppt-1]

  • pol :: [inout] Fractional position of neutral surface within layer KoL of left column

  • por :: [inout] Fractional position of neutral surface within layer KoR of right column

  • kol :: [inout] Index of first left interface above neutral surface

  • kor :: [inout] Index of first right interface above neutral surface

  • heff :: [inout] Effective thickness between two neutral surfaces [R L2 T-2 ~> Pa] or other units following Pl and Pr.

  • bl_kl :: [in] Layer index of the boundary layer (left)

  • bl_kr :: [in] Layer index of the boundary layer (right)

  • bl_zl :: [in] Nondimensional position of the boundary layer (left)

  • bl_zr :: [in] Nondimensional position of the boundary layer (right)

Call to

absolute_position interpolate_for_nondim_position

Called from

ndiff_unit_tests_continuous neutral_diffusion_calc_coeffs

function mom_neutral_diffusion/interpolate_for_nondim_position(dRhoNeg, Pneg, dRhoPos, Ppos) [real]

Returns the non-dimensional position between Pneg and Ppos where the interpolated density difference equals zero. The result is always bounded to be between 0 and 1.

Parameters
  • drhoneg :: [in] Negative density difference [R ~> kg m-3]

  • pneg :: [in] Position of negative density difference [R L2 T-2 ~> Pa] or [nondim]

  • drhopos :: [in] Positive density difference [R ~> kg m-3]

  • ppos :: [in] Position of positive density difference [R L2 T-2 ~> Pa] or [nondim]

Called from

find_neutral_surface_positions_continuous search_other_column test_ifndp

subroutine mom_neutral_diffusion/find_neutral_surface_positions_discontinuous(CS, nk, Pres_l, hcol_l, Tl, Sl, ppoly_T_l, ppoly_S_l, stable_l, Pres_r, hcol_r, Tr, Sr, ppoly_T_r, ppoly_S_r, stable_r, PoL, PoR, KoL, KoR, hEff, zeta_bot_L, zeta_bot_R, k_bot_L, k_bot_R, hard_fail_heff)

Higher order version of find_neutral_surface_positions. Returns positions within left/right columns of combined interfaces using intracell reconstructions of T/S. Note that the polynomial reconstrcutions of T and S are optional to aid with unit testing, but will always be passed otherwise.

Parameters
  • cs :: [inout] Neutral diffusion control structure

  • nk :: [in] Number of levels

  • pres_l :: [in] Left-column interface pressure [R L2 T-2 ~> Pa]

  • hcol_l :: [in] Left-column layer thicknesses [H ~> m or kg m-2] or other units

  • tl :: [in] Left-column top interface potential temperature [degC]

  • sl :: [in] Left-column top interface salinity [ppt]

  • ppoly_t_l :: [in] Left-column coefficients of T reconstruction [degC]

  • ppoly_s_l :: [in] Left-column coefficients of S reconstruction [ppt]

  • stable_l :: [in] True where the left-column is stable

  • pres_r :: [in] Right-column interface pressure [R L2 T-2 ~> Pa]

  • hcol_r :: [in] Left-column layer thicknesses [H ~> m or kg m-2] or other units

  • tr :: [in] Right-column top interface potential temperature [degC]

  • sr :: [in] Right-column top interface salinity [ppt]

  • ppoly_t_r :: [in] Right-column coefficients of T reconstruction [degC]

  • ppoly_s_r :: [in] Right-column coefficients of S reconstruction [ppt]

  • stable_r :: [in] True where the right-column is stable

  • pol :: [inout] Fractional position of neutral surface within layer KoL of left column [nondim]

  • por :: [inout] Fractional position of neutral surface within layer KoR of right column [nondim]

  • kol :: [inout] Index of first left interface above neutral surface

  • kor :: [inout] Index of first right interface above neutral surface

  • heff :: [inout] Effective thickness between two neutral surfaces [H ~> m or kg m-2] or other units taken from hcol_l

  • zeta_bot_l :: [in] Non-dimensional distance to where the boundary layer intersetcs the cell (left) [nondim]

  • zeta_bot_r :: [in] Non-dimensional distance to where the boundary layer intersetcs the cell (right) [nondim]

  • k_bot_l :: [in] k-index for the boundary layer (left) [nondim]

  • k_bot_r :: [in] k-index for the boundary layer (right) [nondim]

  • hard_fail_heff :: [in] If true (default) bring down the model if the neutral surfaces ever cross [logical]

Call to

calc_delta_rho_and_derivs increment_interface search_other_column

Called from

ndiff_unit_tests_discontinuous neutral_diffusion_calc_coeffs

subroutine mom_neutral_diffusion/mark_unstable_cells(CS, nk, T, S, P, stable_cell)

Sweep down through the column and mark as stable if the bottom interface of a cell is denser than the top.

Parameters
  • cs :: [inout] Neutral diffusion control structure

  • nk :: [in] Number of levels in a column

  • t :: [in] Temperature at interfaces [degC]

  • s :: [in] Salinity at interfaces [ppt]

  • p :: [in] Pressure at interfaces [R L2 T-2 ~> Pa]

  • stable_cell :: [out] True if this cell is unstably stratified

Call to

calc_delta_rho_and_derivs

Called from

ndiff_unit_tests_discontinuous neutral_diffusion_calc_coeffs

function mom_neutral_diffusion/search_other_column(CS, ksurf, pos_last, T_from, S_from, P_from, T_top, S_top, P_top, T_bot, S_bot, P_bot, T_poly, S_poly) [real]

Searches the “other” (searched) column for the position of the neutral surface.

Parameters
  • cs :: [in] Neutral diffusion control structure

  • ksurf :: [in] Current index of neutral surface

  • pos_last :: [in] Last position within the current layer, used as the lower bound in the root finding algorithm [nondim]

  • t_from :: [in] Temperature at the searched from interface [degC]

  • s_from :: [in] Salinity at the searched from interface [ppt]

  • p_from :: [in] Pressure at the searched from interface [R L2 T-2 ~> Pa]

  • t_top :: [in] Temperature at the searched to top interface [degC]

  • s_top :: [in] Salinity at the searched to top interface [ppt]

  • p_top :: [in] Pressure at the searched to top interface [R L2 T-2 ~> Pa] interface [R L2 T-2 ~> Pa]

  • t_bot :: [in] Temperature at the searched to bottom interface [degC]

  • s_bot :: [in] Salinity at the searched to bottom interface [ppt]

  • p_bot :: [in] Pressure at the searched to bottom interface [R L2 T-2 ~> Pa]

  • t_poly :: [in] Temperature polynomial reconstruction coefficients [degC]

  • s_poly :: [in] Salinity polynomial reconstruction coefficients [ppt]

Call to

calc_delta_rho_and_derivs find_neutral_pos_full find_neutral_pos_linear interpolate_for_nondim_position

Called from

find_neutral_surface_positions_discontinuous

subroutine mom_neutral_diffusion/increment_interface(nk, kl, ki, reached_bottom, searching_this_column, searching_other_column)

Increments the interface which was just connected and also set flags if the bottom is reached.

Parameters
  • nk :: [in] Number of vertical levels

  • kl :: [inout] Current layer (potentially updated)

  • ki :: [inout] Current interface

  • reached_bottom :: [inout] Updated when kl == nk and ki == 2

  • searching_this_column :: [inout] Updated when kl == nk and ki == 2

  • searching_other_column :: [inout] Updated when kl == nk and ki == 2

Called from

find_neutral_surface_positions_discontinuous

function mom_neutral_diffusion/find_neutral_pos_linear(CS, z0, T_ref, S_ref, dRdT_ref, dRdS_ref, dRdT_top, dRdS_top, dRdT_bot, dRdS_bot, ppoly_T, ppoly_S) [real]

Search a layer to find where delta_rho = 0 based on a linear interpolation of alpha and beta of the top and bottom being searched and polynomial reconstructions of T and S. Compressibility is not needed because either, we are assuming incompressibility in the equation of state for this module or alpha and beta are calculated having been displaced to the average pressures of the two pressures We need Newton’s method because the T and S reconstructions make delta_rho a polynomial function of z if using PPM or higher. If Newton’s method would search fall out of the interval [0,1], a bisection step would be taken instead. Also this linearization of alpha, beta means that second derivatives of the EOS are not needed. Note that delta in variable names below refers to horizontal differences and ‘d’ refers to vertical differences.

Parameters
  • cs :: [in] Control structure with parameters for this module

  • z0 :: [in] Lower bound of position, also serves as the initial guess [nondim]

  • t_ref :: [in] Temperature at the searched from interface [degC]

  • s_ref :: [in] Salinity at the searched from interface [ppt]

  • drdt_ref :: [in] dRho/dT at the searched from interface [R degC-1 ~> kg m-3 degC-1]

  • drds_ref :: [in] dRho/dS at the searched from interface [R ppt-1 ~> kg m-3 ppt-1]

  • drdt_top :: [in] dRho/dT at top of layer being searched [R degC-1 ~> kg m-3 degC-1]

  • drds_top :: [in] dRho/dS at top of layer being searched [R ppt-1 ~> kg m-3 ppt-1]

  • drdt_bot :: [in] dRho/dT at bottom of layer being searched [R degC-1 ~> kg m-3 degC-1]

  • drds_bot :: [in] dRho/dS at bottom of layer being searched [R ppt-1 ~> kg m-3 ppt-1]

  • ppoly_t :: [in] Coefficients of the polynomial reconstruction of T within the layer to be searched [degC].

  • ppoly_s :: [in] Coefficients of the polynomial reconstruction of S within the layer to be searched [ppt].

Return

undefined :: Position where drho = 0 [nondim]

Called from

ndiff_unit_tests_discontinuous search_other_column

function mom_neutral_diffusion/find_neutral_pos_full(CS, z0, T_ref, S_ref, P_ref, P_top, P_bot, ppoly_T, ppoly_S) [real]

Use the full equation of state to calculate the difference in locally referenced potential density. The derivatives in this case are not trivial to calculate, so instead we use a regula falsi method.

Parameters
  • cs :: [in] Control structure with parameters for this module

  • z0 :: [in] Lower bound of position, also serves as the initial guess [nondim]

  • t_ref :: [in] Temperature at the searched from interface [degC]

  • s_ref :: [in] Salinity at the searched from interface [ppt]

  • p_ref :: [in] Pressure at the searched from interface [R L2 T-2 ~> Pa]

  • p_top :: [in] Pressure at top of layer being searched [R L2 T-2 ~> Pa]

  • p_bot :: [in] Pressure at bottom of layer being searched [R L2 T-2 ~> Pa]

  • ppoly_t :: [in] Coefficients of the polynomial reconstruction of T within the layer to be searched [degC]

  • ppoly_s :: [in] Coefficients of the polynomial reconstruction of T within the layer to be searched [ppt]

Return

undefined :: Position where drho = 0 [nondim]

Call to

calc_delta_rho_and_derivs

Called from

search_other_column

subroutine mom_neutral_diffusion/calc_delta_rho_and_derivs(CS, T1, S1, p1_in, T2, S2, p2_in, drho, drdt1_out, drds1_out, drdt2_out, drds2_out)

Calculate the difference in density between two points in a variety of ways.

Parameters
  • cs :: Neutral diffusion control structure

  • t1 :: [in] Temperature at point 1 [degC]

  • s1 :: [in] Salinity at point 1 [ppt]

  • p1_in :: [in] Pressure at point 1 [R L2 T-2 ~> Pa]

  • t2 :: [in] Temperature at point 2 [degC]

  • s2 :: [in] Salinity at point 2 [ppt]

  • p2_in :: [in] Pressure at point 2 [R L2 T-2 ~> Pa]

  • drho :: [out] Difference in density between the two points [R ~> kg m-3]

  • drdt1_out :: [out] drho_dt at point 1 [R degC-1 ~> kg m-3 degC-1]

  • drds1_out :: [out] drho_ds at point 1 [R ppt-1 ~> kg m-3 ppt-1]

  • drdt2_out :: [out] drho_dt at point 2 [R degC-1 ~> kg m-3 degC-1]

  • drds2_out :: [out] drho_ds at point 2 [R ppt-1 ~> kg m-3 ppt-1]

Call to

delta_rho_from_derivs

Called from

find_neutral_pos_full find_neutral_surface_positions_discontinuous mark_unstable_cells search_other_column

function mom_neutral_diffusion/delta_rho_from_derivs(T1, S1, P1, dRdT1, dRdS1, T2, S2, P2, dRdT2, dRdS2) [real]

Calculate delta rho from derivatives and gradients of properties \(\Delta \rho = \frac{1}{2}\left[ (\alpha_1 + \alpha_2)*(T_1-T_2) + (\beta_1 + \beta_2)*(S_1-S_2) + (\gamma^{-1}_1 + \gamma^{-1}_2)*(P_1-P_2) \right]\).

Parameters
  • t1 :: Temperature at point 1 [degC]

  • s1 :: Salinity at point 1 [ppt]

  • p1 :: Pressure at point 1 [R L2 T-2 ~> Pa]

  • drdt1 :: The partial derivative of density with temperature at point 1 [R degC-1 ~> kg m-3 degC-1]

  • drds1 :: The partial derivative of density with salinity at point 1 [R ppt-1 ~> kg m-3 ppt-1]

  • t2 :: Temperature at point 2 [degC]

  • s2 :: Salinity at point 2 [ppt]

  • p2 :: Pressure at point 2 [R L2 T-2 ~> Pa]

  • drdt2 :: The partial derivative of density with temperature at point 2 [R degC-1 ~> kg m-3 degC-1]

  • drds2 :: The partial derivative of density with salinity at point 2 [R ppt-1 ~> kg m-3 ppt-1]

Called from

calc_delta_rho_and_derivs

function mom_neutral_diffusion/absolute_position(n, ns, Pint, Karr, NParr, k_surface) [real]

Converts non-dimensional position within a layer to absolute position (for debugging)

Parameters
  • n :: [in] Number of levels

  • ns :: [in] Number of neutral surfaces

  • pint :: [in] Position of interfaces [R L2 T-2 ~> Pa] or other units

  • karr :: [in] Index of interface above position

  • nparr :: [in] Non-dimensional position within layer Karr(:) [nondim]

  • k_surface :: [in] k-interface to query

Return

undefined :: The absolute position of a location [R L2 T-2 ~> Pa] or other units following Pint

Called from

absolute_positions find_neutral_surface_positions_continuous

function mom_neutral_diffusion/absolute_positions(n, ns, Pint, Karr, NParr) [real]

Converts non-dimensional positions within layers to absolute positions (for debugging)

Parameters
  • n :: [in] Number of levels

  • ns :: [in] Number of neutral surfaces

  • pint :: [in] Position of interface [R L2 T-2 ~> Pa] or other units

  • karr :: [in] Indexes of interfaces about positions

  • nparr :: [in] Non-dimensional positions within layers Karr(:)

Return

undefined :: Absolute positions [R L2 T-2 ~> Pa] or other units following Pint

Call to

absolute_position

Called from

ndiff_unit_tests_continuous

subroutine mom_neutral_diffusion/neutral_surface_flux(nk, nsurf, deg, hl, hr, Tl, Tr, PiL, PiR, KoL, KoR, hEff, Flx, continuous, h_neglect, remap_CS, h_neglect_edge)

Returns a single column of neutral diffusion fluxes of a tracer.

Parameters
  • nk :: [in] Number of levels

  • nsurf :: [in] Number of neutral surfaces

  • deg :: [in] Degree of polynomial reconstructions

  • hl :: [in] Left-column layer thickness [H ~> m or kg m-2]

  • hr :: [in] Right-column layer thickness [H ~> m or kg m-2]

  • tl :: [in] Left-column layer tracer (conc, e.g. degC)

  • tr :: [in] Right-column layer tracer (conc, e.g. degC)

  • pil :: [in] Fractional position of neutral surface within layer KoL of left column

  • pir :: [in] Fractional position of neutral surface within layer KoR of right column

  • kol :: [in] Index of first left interface above neutral surface

  • kor :: [in] Index of first right interface above neutral surface

  • heff :: [in] Effective thickness between two neutral surfaces [H ~> m or kg m-2]

  • flx :: [inout] Flux of tracer between pairs of neutral layers (conc H)

  • continuous :: [in] True if using continuous reconstruction

  • h_neglect :: [in] A negligibly small width for the purpose of cell reconstructions [H ~> m or kg m-2]

  • remap_cs :: [in] Remapping control structure used to create sublayers

  • h_neglect_edge :: [in] A negligibly small width used for edge value calculations if continuous is false [H ~> m or kg m-2]

Call to

interface_scalar neutral_surface_t_eval ppm_ave ppm_left_right_edge_values signum

Called from

ndiff_unit_tests_continuous neutral_diffusion

subroutine mom_neutral_diffusion/neutral_surface_t_eval(nk, ns, k_sub, Ks, Ps, T_mean, T_int, deg, iMethod, T_poly, T_top, T_bot, T_sub, T_top_int, T_bot_int, T_layer)

Evaluate various parts of the reconstructions to calculate gradient-based flux limter.

Parameters
  • nk :: [in] Number of cell everages

  • ns :: [in] Number of neutral surfaces

  • k_sub :: [in] Index of current neutral layer

  • ks :: [in] List of the layers associated with each neutral surface

  • ps :: [in] List of the positions within a layer of each surface

  • t_mean :: [in] Cell average of tracer

  • t_int :: [in] Cell interface values of tracer from reconstruction

  • deg :: [in] Degree of reconstruction polynomial (e.g. 1 is linear)

  • imethod :: [in] Method of integration to use

  • t_poly :: [in] Coefficients of polynomial reconstructions

  • t_top :: [out] Tracer value at top (across discontinuity if necessary)

  • t_bot :: [out] Tracer value at bottom (across discontinuity if necessary)

  • t_sub :: [out] Average of the tracer value over the sublayer

  • t_top_int :: [out] Tracer value at top interface of neutral layer

  • t_bot_int :: [out] Tracer value at bottom interface of neutral layer

  • t_layer :: [out] Cell-average that the the reconstruction belongs to

Called from

neutral_surface_flux

subroutine mom_neutral_diffusion/ppm_left_right_edge_values(nk, Tl, Ti, aL, aR)

Discontinuous PPM reconstructions of the left/right edge values within a cell.

Parameters
  • nk :: [in] Number of levels

  • tl :: [in] Layer tracer (conc, e.g. degC)

  • ti :: [in] Interface tracer (conc, e.g. degC)

  • al :: [inout] Left edge value of tracer (conc, e.g. degC)

  • ar :: [inout] Right edge value of tracer (conc, e.g. degC)

Call to

signum

Called from

neutral_surface_flux

function mom_neutral_diffusion/neutral_diffusion_unit_tests(verbose) [logical]

Returns true if unit tests of neutral_diffusion functions fail. Otherwise returns false.

Parameters

verbose :: [in] If true, write results to stdout

Call to

ndiff_unit_tests_continuous ndiff_unit_tests_discontinuous

Called from

mom_unit_tests::unit_tests

function mom_neutral_diffusion/ndiff_unit_tests_continuous(verbose) [logical]

Returns true if unit tests of neutral_diffusion functions fail. Otherwise returns false.

Parameters

verbose :: [in] If true, write results to stdout

Call to

absolute_positions find_neutral_surface_positions_continuous interface_scalar neutral_surface_flux test_data1d test_fv_diff test_fvlsq_slope test_ifndp test_nsp

Called from

neutral_diffusion_unit_tests

function mom_neutral_diffusion/ndiff_unit_tests_discontinuous(verbose) [logical]
Parameters

verbose :: [in] It true, write results to stdout

Call to

find_neutral_pos_linear find_neutral_surface_positions_discontinuous mark_unstable_cells test_nsp test_rnp

Called from

neutral_diffusion_unit_tests

function mom_neutral_diffusion/test_fv_diff(verbose, hkm1, hk, hkp1, Skm1, Sk, Skp1, Ptrue, title) [logical]

Returns true if a test of fv_diff() fails, and conditionally writes results to stream. fails, and conditionally writes results to stream.

Parameters
  • verbose :: [in] If true, write results to stdout

  • hkm1 :: [in] Left cell width [nondim]

  • hk :: [in] Center cell width [nondim]

  • hkp1 :: [in] Right cell width [nondim]

  • skm1 :: [in] Left cell average value

  • sk :: [in] Center cell average value

  • skp1 :: [in] Right cell average value

  • ptrue :: [in] True answer [nondim]

  • title :: [in] Title for messages

Call to

fv_diff

Called from

ndiff_unit_tests_continuous

function mom_neutral_diffusion/test_fvlsq_slope(verbose, hkm1, hk, hkp1, Skm1, Sk, Skp1, Ptrue, title) [logical]

Returns true if a test of fvlsq_slope() fails, and conditionally writes results to stream. fails, and conditionally writes results to stream.

Parameters
  • verbose :: [in] If true, write results to stdout

  • hkm1 :: [in] Left cell width

  • hk :: [in] Center cell width

  • hkp1 :: [in] Right cell width

  • skm1 :: [in] Left cell average value

  • sk :: [in] Center cell average value

  • skp1 :: [in] Right cell average value

  • ptrue :: [in] True answer

  • title :: [in] Title for messages

Call to

fvlsq_slope

Called from

ndiff_unit_tests_continuous

function mom_neutral_diffusion/test_ifndp(verbose, rhoNeg, Pneg, rhoPos, Ppos, Ptrue, title) [logical]

Returns true if a test of interpolate_for_nondim_position() fails, and conditionally writes results to stream. fails, and conditionally writes results to stream.

Parameters
  • verbose :: [in] If true, write results to stdout

  • rhoneg :: [in] Lighter density [R ~> kg m-3]

  • pneg :: [in] Interface position of lighter density [nondim]

  • rhopos :: [in] Heavier density [R ~> kg m-3]

  • ppos :: [in] Interface position of heavier density [nondim]

  • ptrue :: [in] True answer [nondim]

  • title :: [in] Title for messages

Call to

interpolate_for_nondim_position

Called from

ndiff_unit_tests_continuous

function mom_neutral_diffusion/test_data1d(verbose, nk, Po, Ptrue, title) [logical]

Returns true if comparison of Po and Ptrue fails, and conditionally writes results to stream.

Parameters
  • verbose :: [in] If true, write results to stdout

  • nk :: [in] Number of layers

  • po :: [in] Calculated answer

  • ptrue :: [in] True answer

  • title :: [in] Title for messages

Called from

ndiff_unit_tests_continuous

function mom_neutral_diffusion/test_data1di(verbose, nk, Po, Ptrue, title) [logical]

Returns true if comparison of Po and Ptrue fails, and conditionally writes results to stream.

Parameters
  • verbose :: [in] If true, write results to stdout

  • nk :: [in] Number of layers

  • po :: [in] Calculated answer

  • ptrue :: [in] True answer

  • title :: [in] Title for messages

function mom_neutral_diffusion/test_nsp(verbose, ns, KoL, KoR, pL, pR, hEff, KoL0, KoR0, pL0, pR0, hEff0, title) [logical]

Returns true if output of find_neutral_surface_positions() does not match correct values, and conditionally writes results to stream.

Parameters
  • verbose :: [in] If true, write results to stdout

  • ns :: [in] Number of surfaces

  • kol :: [in] Index of first left interface above neutral surface

  • kor :: [in] Index of first right interface above neutral surface

  • pl :: [in] Fractional position of neutral surface within layer KoL of left column

  • pr :: [in] Fractional position of neutral surface within layer KoR of right column

  • heff :: [in] Effective thickness between two neutral surfaces [R L2 T-2 ~> Pa]

  • kol0 :: [in] Correct value for KoL

  • kor0 :: [in] Correct value for KoR

  • pl0 :: [in] Correct value for pL

  • pr0 :: [in] Correct value for pR

  • heff0 :: [in] Correct value for hEff

  • title :: [in] Title for messages

Call to

compare_nsp_row

Called from

ndiff_unit_tests_continuous ndiff_unit_tests_discontinuous

function mom_neutral_diffusion/compare_nsp_row(KoL, KoR, pL, pR, KoL0, KoR0, pL0, pR0) [logical]

Compares a single row, k, of output from find_neutral_surface_positions()

Parameters
  • kol :: [in] Index of first left interface above neutral surface

  • kor :: [in] Index of first right interface above neutral surface

  • pl :: [in] Fractional position of neutral surface within layer KoL of left column

  • pr :: [in] Fractional position of neutral surface within layer KoR of right column

  • kol0 :: [in] Correct value for KoL

  • kor0 :: [in] Correct value for KoR

  • pl0 :: [in] Correct value for pL

  • pr0 :: [in] Correct value for pR

Called from

test_nsp

function mom_neutral_diffusion/test_rnp(expected_pos, test_pos, title) [logical]

Compares output position from refine_nondim_position with an expected value.

Parameters
  • expected_pos :: [in] The expected position

  • test_pos :: [in] The position returned by the code

  • title :: [in] A label for this test

Called from

ndiff_unit_tests_discontinuous

subroutine mom_neutral_diffusion/neutral_diffusion_end(CS)

Deallocates neutral_diffusion control structure.

Parameters

cs :: Neutral diffusion control structure