MOM6
mom_neutral_diffusion Module Reference

Detailed Description

A column-wise toolbox for implementing neutral diffusion.

Data Types

type  neutral_diffusion_cs
 The control structure for the MOM_neutral_diffusion module. More...
 

Functions/Subroutines

logical function, public neutral_diffusion_init (Time, G, param_file, diag, EOS, diabatic_CSp, CS)
 Read parameters and allocate control structure for neutral_diffusion module. More...
 
subroutine, public neutral_diffusion_calc_coeffs (G, GV, US, h, T, S, CS)
 Calculate remapping factors for u/v columns used to map adjoining columns to a shared coordinate space. More...
 
subroutine, public 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. More...
 
subroutine interface_scalar (nk, h, S, Si, i_method, h_neglect)
 Returns interface scalar, Si, for a column of layer values, S. More...
 
real function ppm_edge (hkm1, hk, hkp1, hkp2, Ak, Akp1, Pk, Pkp1, h_neglect)
 Returns the PPM quasi-fourth order edge value at k+1/2 following equation 1.6 in Colella & Woodward, 1984: JCP 54, 174-201. More...
 
real function ppm_ave (xL, xR, aL, aR, aMean)
 Returns the average of a PPM reconstruction between two fractional positions. More...
 
real function signum (a, x)
 A true signum function that returns either -abs(a), when x<0; or abs(a) when x>0; or 0 when x=0. More...
 
subroutine 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. More...
 
real function fv_diff (hkm1, hk, hkp1, Skm1, Sk, Skp1)
 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. More...
 
real function fvlsq_slope (hkm1, hk, hkp1, Skm1, Sk, Skp1)
 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). More...
 
subroutine 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. More...
 
real function interpolate_for_nondim_position (dRhoNeg, Pneg, dRhoPos, Ppos)
 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. More...
 
subroutine 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. More...
 
subroutine 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. More...
 
real function 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)
 Searches the "other" (searched) column for the position of the neutral surface. More...
 
subroutine 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. More...
 
real function find_neutral_pos_linear (CS, z0, T_ref, S_ref, P_ref, dRdT_ref, dRdS_ref, P_top, dRdT_top, dRdS_top, P_bot, dRdT_bot, dRdS_bot, ppoly_T, ppoly_S)
 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. More...
 
real function find_neutral_pos_full (CS, z0, T_ref, S_ref, P_ref, P_top, P_bot, ppoly_T, ppoly_S)
 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. More...
 
subroutine 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. More...
 
real function delta_rho_from_derivs (T1, S1, P1, dRdT1, dRdS1, T2, S2, P2, dRdT2, dRdS2)
 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] \). More...
 
real function absolute_position (n, ns, Pint, Karr, NParr, k_surface)
 Converts non-dimensional position within a layer to absolute position (for debugging) More...
 
real function, dimension(ns) absolute_positions (n, ns, Pint, Karr, NParr)
 Converts non-dimensional positions within layers to absolute positions (for debugging) More...
 
subroutine 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. More...
 
subroutine 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. More...
 
subroutine ppm_left_right_edge_values (nk, Tl, Ti, aL, aR)
 Discontinuous PPM reconstructions of the left/right edge values within a cell. More...
 
logical function, public neutral_diffusion_unit_tests (verbose)
 Returns true if unit tests of neutral_diffusion functions fail. Otherwise returns false. More...
 
logical function ndiff_unit_tests_continuous (verbose)
 Returns true if unit tests of neutral_diffusion functions fail. Otherwise returns false. More...
 
logical function ndiff_unit_tests_discontinuous (verbose)
 
logical function test_fv_diff (verbose, hkm1, hk, hkp1, Skm1, Sk, Skp1, Ptrue, title)
 Returns true if a test of fv_diff() fails, and conditionally writes results to stream. More...
 
logical function test_fvlsq_slope (verbose, hkm1, hk, hkp1, Skm1, Sk, Skp1, Ptrue, title)
 Returns true if a test of fvlsq_slope() fails, and conditionally writes results to stream. More...
 
logical function test_ifndp (verbose, rhoNeg, Pneg, rhoPos, Ppos, Ptrue, title)
 Returns true if a test of interpolate_for_nondim_position() fails, and conditionally writes results to stream. More...
 
logical function test_data1d (verbose, nk, Po, Ptrue, title)
 Returns true if comparison of Po and Ptrue fails, and conditionally writes results to stream. More...
 
logical function test_data1di (verbose, nk, Po, Ptrue, title)
 Returns true if comparison of Po and Ptrue fails, and conditionally writes results to stream. More...
 
logical function test_nsp (verbose, ns, KoL, KoR, pL, pR, hEff, KoL0, KoR0, pL0, pR0, hEff0, title)
 Returns true if output of find_neutral_surface_positions() does not match correct values, and conditionally writes results to stream. More...
 
logical function compare_nsp_row (KoL, KoR, pL, pR, KoL0, KoR0, pL0, pR0)
 Compares a single row, k, of output from find_neutral_surface_positions() More...
 
logical function test_rnp (expected_pos, test_pos, title)
 Compares output position from refine_nondim_position with an expected value. More...
 
subroutine, public neutral_diffusion_end (CS)
 Deallocates neutral_diffusion control structure. More...
 

Variables

character(len=40) mdl = "MOM_neutral_diffusion"
 module name More...
 

Function/Subroutine Documentation

◆ absolute_position()

real function mom_neutral_diffusion::absolute_position ( integer, intent(in)  n,
integer, intent(in)  ns,
real, dimension(n+1), intent(in)  Pint,
integer, dimension(ns), intent(in)  Karr,
real, dimension(ns), intent(in)  NParr,
integer, intent(in)  k_surface 
)
private

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

Parameters
[in]nNumber of levels
[in]nsNumber of neutral surfaces
[in]pintPosition of interfaces [Pa]
[in]karrIndex of interface above position
[in]nparrNon-dimensional position within layer Karr(:)
[in]k_surfacek-interface to query

Definition at line 1783 of file MOM_neutral_diffusion.F90.

1783  integer, intent(in) :: n !< Number of levels
1784  integer, intent(in) :: ns !< Number of neutral surfaces
1785  real, intent(in) :: Pint(n+1) !< Position of interfaces [Pa]
1786  integer, intent(in) :: Karr(ns) !< Index of interface above position
1787  real, intent(in) :: NParr(ns) !< Non-dimensional position within layer Karr(:)
1788  integer, intent(in) :: k_surface !< k-interface to query
1789  ! Local variables
1790  integer :: k
1791 
1792  k = karr(k_surface)
1793  if (k>n) stop 'absolute_position: k>nk is out of bounds!'
1794  absolute_position = pint(k) + nparr(k_surface) * ( pint(k+1) - pint(k) )
1795 

Referenced by absolute_positions(), and find_neutral_surface_positions_continuous().

Here is the caller graph for this function:

◆ absolute_positions()

real function, dimension(ns) mom_neutral_diffusion::absolute_positions ( integer, intent(in)  n,
integer, intent(in)  ns,
real, dimension(n+1), intent(in)  Pint,
integer, dimension(ns), intent(in)  Karr,
real, dimension(ns), intent(in)  NParr 
)
private

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

Parameters
[in]nNumber of levels
[in]nsNumber of neutral surfaces
[in]pintPosition of interface [Pa]
[in]karrIndexes of interfaces about positions
[in]nparrNon-dimensional positions within layers Karr(:)

Definition at line 1800 of file MOM_neutral_diffusion.F90.

1800  integer, intent(in) :: n !< Number of levels
1801  integer, intent(in) :: ns !< Number of neutral surfaces
1802  real, intent(in) :: Pint(n+1) !< Position of interface [Pa]
1803  integer, intent(in) :: Karr(ns) !< Indexes of interfaces about positions
1804  real, intent(in) :: NParr(ns) !< Non-dimensional positions within layers Karr(:)
1805 
1806  real, dimension(ns) :: absolute_positions ! Absolute positions [Pa]
1807 
1808  ! Local variables
1809  integer :: k_surface, k
1810 
1811  do k_surface = 1, ns
1812  absolute_positions(k_surface) = absolute_position(n,ns,pint,karr,nparr,k_surface)
1813  enddo
1814 

References absolute_position().

Referenced by ndiff_unit_tests_continuous().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ calc_delta_rho_and_derivs()

subroutine mom_neutral_diffusion::calc_delta_rho_and_derivs ( type(neutral_diffusion_cs CS,
real, intent(in)  T1,
real, intent(in)  S1,
real, intent(in)  p1_in,
real, intent(in)  T2,
real, intent(in)  S2,
real, intent(in)  p2_in,
real, intent(out)  drho,
real, intent(out), optional  drdt1_out,
real, intent(out), optional  drds1_out,
real, intent(out), optional  drdt2_out,
real, intent(out), optional  drds2_out 
)
private

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

Parameters
csNeutral diffusion control structure
[in]t1Temperature at point 1
[in]s1Salinity at point 1
[in]p1_inPressure at point 1
[in]t2Temperature at point 2
[in]s2Salinity at point 2
[in]p2_inPressure at point 2
[out]drhoDifference in density between the two points
[out]drdt1_outdrho_dt at point 1
[out]drds1_outdrho_ds at point 1
[out]drdt2_outdrho_dt at point 2
[out]drds2_outdrho_ds at point 2

Definition at line 1706 of file MOM_neutral_diffusion.F90.

1706  type(neutral_diffusion_CS) :: CS !< Neutral diffusion control structure
1707  real, intent(in ) :: T1 !< Temperature at point 1
1708  real, intent(in ) :: S1 !< Salinity at point 1
1709  real, intent(in ) :: p1_in !< Pressure at point 1
1710  real, intent(in ) :: T2 !< Temperature at point 2
1711  real, intent(in ) :: S2 !< Salinity at point 2
1712  real, intent(in ) :: p2_in !< Pressure at point 2
1713  real, intent( out) :: drho !< Difference in density between the two points
1714  real, optional, intent( out) :: dRdT1_out !< drho_dt at point 1
1715  real, optional, intent( out) :: dRdS1_out !< drho_ds at point 1
1716  real, optional, intent( out) :: dRdT2_out !< drho_dt at point 2
1717  real, optional, intent( out) :: dRdS2_out !< drho_ds at point 2
1718  ! Local variables
1719  real :: rho1, rho2, p1, p2, pmid
1720  real :: drdt1, drdt2, drds1, drds2, drdp1, drdp2, rho_dummy
1721 
1722  ! Use the same reference pressure or the in-situ pressure
1723  if (cs%ref_pres > 0.) then
1724  p1 = cs%ref_pres
1725  p2 = cs%ref_pres
1726  else
1727  p1 = p1_in
1728  p2 = p2_in
1729  endif
1730 
1731  ! Use the full linear equation of state to calculate the difference in density (expensive!)
1732  if (trim(cs%delta_rho_form) == 'full') then
1733  pmid = 0.5 * (p1 + p2)
1734  call calculate_density( t1, s1, pmid, rho1, cs%EOS )
1735  call calculate_density( t2, s2, pmid, rho2, cs%EOS )
1736  drho = rho1 - rho2
1737  ! Use the density derivatives at the average of pressures and the differentces int temperature
1738  elseif (trim(cs%delta_rho_form) == 'mid_pressure') then
1739  pmid = 0.5 * (p1 + p2)
1740  if (cs%ref_pres>=0) pmid = cs%ref_pres
1741  call calculate_density_derivs(t1, s1, pmid, drdt1, drds1, cs%EOS)
1742  call calculate_density_derivs(t2, s2, pmid, drdt2, drds2, cs%EOS)
1743  drho = delta_rho_from_derivs( t1, s1, p1, drdt1, drds1, t2, s2, p2, drdt2, drds2)
1744  elseif (trim(cs%delta_rho_form) == 'local_pressure') then
1745  call calculate_density_derivs(t1, s1, p1, drdt1, drds1, cs%EOS)
1746  call calculate_density_derivs(t2, s2, p2, drdt2, drds2, cs%EOS)
1747  drho = delta_rho_from_derivs( t1, s1, p1, drdt1, drds1, t2, s2, p2, drdt2, drds2)
1748  else
1749  call mom_error(fatal, "delta_rho_form is not recognized")
1750  endif
1751 
1752  if (PRESENT(drdt1_out)) drdt1_out = drdt1
1753  if (PRESENT(drds1_out)) drds1_out = drds1
1754  if (PRESENT(drdt2_out)) drdt2_out = drdt2
1755  if (PRESENT(drds2_out)) drds2_out = drds2
1756 

References delta_rho_from_derivs(), and mom_error_handler::mom_error().

Referenced by find_neutral_pos_full(), find_neutral_surface_positions_discontinuous(), mark_unstable_cells(), and search_other_column().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ compare_nsp_row()

logical function mom_neutral_diffusion::compare_nsp_row ( integer, intent(in)  KoL,
integer, intent(in)  KoR,
real, intent(in)  pL,
real, intent(in)  pR,
integer, intent(in)  KoL0,
integer, intent(in)  KoR0,
real, intent(in)  pL0,
real, intent(in)  pR0 
)
private

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

Parameters
[in]kolIndex of first left interface above neutral surface
[in]korIndex of first right interface above neutral surface
[in]plFractional position of neutral surface within layer KoL of left column
[in]prFractional position of neutral surface within layer KoR of right column
[in]kol0Correct value for KoL
[in]kor0Correct value for KoR
[in]pl0Correct value for pL
[in]pr0Correct value for pR

Definition at line 2775 of file MOM_neutral_diffusion.F90.

2775  integer, intent(in) :: KoL !< Index of first left interface above neutral surface
2776  integer, intent(in) :: KoR !< Index of first right interface above neutral surface
2777  real, intent(in) :: pL !< Fractional position of neutral surface within layer KoL of left column
2778  real, intent(in) :: pR !< Fractional position of neutral surface within layer KoR of right column
2779  integer, intent(in) :: KoL0 !< Correct value for KoL
2780  integer, intent(in) :: KoR0 !< Correct value for KoR
2781  real, intent(in) :: pL0 !< Correct value for pL
2782  real, intent(in) :: pR0 !< Correct value for pR
2783 
2784  compare_nsp_row = .false.
2785  if (kol /= kol0) compare_nsp_row = .true.
2786  if (kor /= kor0) compare_nsp_row = .true.
2787  if (pl /= pl0) compare_nsp_row = .true.
2788  if (pr /= pr0) compare_nsp_row = .true.

Referenced by test_nsp().

Here is the caller graph for this function:

◆ delta_rho_from_derivs()

real function mom_neutral_diffusion::delta_rho_from_derivs ( real  T1,
real  S1,
real  P1,
real  dRdT1,
real  dRdS1,
real  T2,
real  S2,
real  P2,
real  dRdT2,
real  dRdS2 
)
private

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
t1Temperature at point 1
s1Salinity at point 1
p1Pressure at point 1
drdt1Pressure at point 1
drds1Pressure at point 1
t2Temperature at point 2
s2Salinity at point 2
p2Pressure at point 2
drdt2Pressure at point 2
drds2Pressure at point 2

Definition at line 1765 of file MOM_neutral_diffusion.F90.

1765  real :: T1 !< Temperature at point 1
1766  real :: S1 !< Salinity at point 1
1767  real :: P1 !< Pressure at point 1
1768  real :: dRdT1 !< Pressure at point 1
1769  real :: dRdS1 !< Pressure at point 1
1770  real :: T2 !< Temperature at point 2
1771  real :: S2 !< Salinity at point 2
1772  real :: P2 !< Pressure at point 2
1773  real :: dRdT2 !< Pressure at point 2
1774  real :: dRdS2 !< Pressure at point 2
1775  ! Local variables
1776  real :: drho
1777 
1778  drho = 0.5 * ( (drdt1+drdt2)*(t1-t2) + (drds1+drds2)*(s1-s2))
1779 

Referenced by calc_delta_rho_and_derivs(), and find_neutral_pos_linear().

Here is the caller graph for this function:

◆ find_neutral_pos_full()

real function mom_neutral_diffusion::find_neutral_pos_full ( type(neutral_diffusion_cs), intent(in)  CS,
real, intent(in)  z0,
real, intent(in)  T_ref,
real, intent(in)  S_ref,
real, intent(in)  P_ref,
real, intent(in)  P_top,
real, intent(in)  P_bot,
real, dimension(:), intent(in)  ppoly_T,
real, dimension(:), intent(in)  ppoly_S 
)
private

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
[in]csControl structure with parameters for this module
[in]z0Lower bound of position, also serves as the initial guess
[in]t_refTemperature at the searched from interface
[in]s_refSalinity at the searched from interface
[in]p_refPressure at the searched from interface
[in]p_topPressure at top of layer being searched
[in]p_botPressure at bottom of layer being searched
[in]ppoly_tCoefficients of the polynomial reconstruction of T within the layer to be searched.
[in]ppoly_sCoefficients of the polynomial reconstruction of T within the layer to be searched.
Returns
Position where drho = 0

Definition at line 1616 of file MOM_neutral_diffusion.F90.

1616  type(neutral_diffusion_CS),intent(in) :: CS !< Control structure with parameters for this module
1617  real, intent(in) :: z0 !< Lower bound of position, also serves as the initial guess
1618  real, intent(in) :: T_ref !< Temperature at the searched from interface
1619  real, intent(in) :: S_ref !< Salinity at the searched from interface
1620  real, intent(in) :: P_ref !< Pressure at the searched from interface
1621  real, intent(in) :: P_top !< Pressure at top of layer being searched
1622  real, intent(in) :: P_bot !< Pressure at bottom of layer being searched
1623  real, dimension(:), intent(in) :: ppoly_T !< Coefficients of the polynomial reconstruction of T within
1624  !! the layer to be searched.
1625  real, dimension(:), intent(in) :: ppoly_S !< Coefficients of the polynomial reconstruction of T within
1626  !! the layer to be searched.
1627  real :: z !< Position where drho = 0
1628  ! Local variables
1629  integer :: iter
1630  integer :: nterm
1631 
1632  real :: drho_a, drho_b, drho_c
1633  real :: a, b, c, Ta, Tb, Tc, Sa, Sb, Sc, Pa, Pb, Pc
1634  integer :: side
1635 
1636  side = 0
1637  ! Set the first two evaluation to the endpoints of the interval
1638  b = z0; c = 1
1639  nterm = SIZE(ppoly_t)
1640 
1641  ! Calculate drho at the minimum bound
1642  tb = evaluation_polynomial( ppoly_t, nterm, b )
1643  sb = evaluation_polynomial( ppoly_s, nterm, b )
1644  pb = p_top*(1.-b) + p_bot*b
1645  call calc_delta_rho_and_derivs(cs, tb, sb, pb, t_ref, s_ref, p_ref, drho_b)
1646 
1647  ! Calculate drho at the maximum bound
1648  tc = evaluation_polynomial( ppoly_t, nterm, 1. )
1649  sc = evaluation_polynomial( ppoly_s, nterm, 1. )
1650  pc = p_bot
1651  call calc_delta_rho_and_derivs(cs, tc, sc, pc, t_ref, s_ref, p_ref, drho_c)
1652 
1653  if (drho_b >= 0.) then
1654  z = z0
1655  return
1656  elseif (drho_c == 0.) then
1657  z = 1.
1658  return
1659  endif
1660  if ( sign(1.,drho_b) == sign(1.,drho_c) ) then
1661  z = z0
1662  return
1663  endif
1664 
1665  do iter = 1, cs%max_iter
1666  ! Calculate new position and evaluate if we have converged
1667  a = (drho_b*c - drho_c*b)/(drho_b-drho_c)
1668  ta = evaluation_polynomial( ppoly_t, nterm, a )
1669  sa = evaluation_polynomial( ppoly_s, nterm, a )
1670  pa = p_top*(1.-a) + p_bot*a
1671  call calc_delta_rho_and_derivs(cs, ta, sa, pa, t_ref, s_ref, p_ref, drho_a)
1672  if (abs(drho_a) < cs%drho_tol) then
1673  z = a
1674  return
1675  endif
1676 
1677  if (drho_a*drho_c > 0.) then
1678  if ( abs(a-c)<cs%x_tol) then
1679  z = a
1680  return
1681  endif
1682  c = a ; drho_c = drho_a;
1683  if (side == -1) drho_b = 0.5*drho_b
1684  side = -1
1685  elseif ( drho_b*drho_a > 0 ) then
1686  if ( abs(a-b)<cs%x_tol) then
1687  z = a
1688  return
1689  endif
1690  b = a ; drho_b = drho_a
1691  if (side == 1) drho_c = 0.5*drho_c
1692  side = 1
1693  else
1694  z = a
1695  return
1696  endif
1697  enddo
1698 
1699  z = a
1700 

References calc_delta_rho_and_derivs(), and polynomial_functions::evaluation_polynomial().

Referenced by search_other_column().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ find_neutral_pos_linear()

real function mom_neutral_diffusion::find_neutral_pos_linear ( type(neutral_diffusion_cs), intent(in)  CS,
real, intent(in)  z0,
real, intent(in)  T_ref,
real, intent(in)  S_ref,
real, intent(in)  P_ref,
real, intent(in)  dRdT_ref,
real, intent(in)  dRdS_ref,
real, intent(in)  P_top,
real, intent(in)  dRdT_top,
real, intent(in)  dRdS_top,
real, intent(in)  P_bot,
real, intent(in)  dRdT_bot,
real, intent(in)  dRdS_bot,
real, dimension(:), intent(in)  ppoly_T,
real, dimension(:), intent(in)  ppoly_S 
)
private

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
[in]csControl structure with parameters for this module
[in]z0Lower bound of position, also serves as the initial guess
[in]t_refTemperature at the searched from interface
[in]s_refSalinity at the searched from interface
[in]p_refPressure at the searched from interface
[in]drdt_refdRho/dT at the searched from interface
[in]drds_refdRho/dS at the searched from interface
[in]p_topPressure at top of layer being searched
[in]drdt_topdRho/dT at top of layer being searched
[in]drds_topdRho/dS at top of layer being searched
[in]p_botPressure at bottom of layer being searched
[in]drdt_botdRho/dT at bottom of layer being searched
[in]drds_botdRho/dS at bottom of layer being searched
[in]ppoly_tCoefficients of the polynomial reconstruction of T within the layer to be searched.
[in]ppoly_sCoefficients of the polynomial reconstruction of T within the layer to be searched.
Returns
Position where drho = 0

Definition at line 1497 of file MOM_neutral_diffusion.F90.

1497  type(neutral_diffusion_CS),intent(in) :: CS !< Control structure with parameters for this module
1498  real, intent(in) :: z0 !< Lower bound of position, also serves as the initial guess
1499  real, intent(in) :: T_ref !< Temperature at the searched from interface
1500  real, intent(in) :: S_ref !< Salinity at the searched from interface
1501  real, intent(in) :: P_ref !< Pressure at the searched from interface
1502  real, intent(in) :: dRdT_ref !< dRho/dT at the searched from interface
1503  real, intent(in) :: dRdS_ref !< dRho/dS at the searched from interface
1504  real, intent(in) :: P_top !< Pressure at top of layer being searched
1505  real, intent(in) :: dRdT_top !< dRho/dT at top of layer being searched
1506  real, intent(in) :: dRdS_top !< dRho/dS at top of layer being searched
1507  real, intent(in) :: P_bot !< Pressure at bottom of layer being searched
1508  real, intent(in) :: dRdT_bot !< dRho/dT at bottom of layer being searched
1509  real, intent(in) :: dRdS_bot !< dRho/dS at bottom of layer being searched
1510  real, dimension(:), intent(in) :: ppoly_T !< Coefficients of the polynomial reconstruction of T within
1511  !! the layer to be searched.
1512  real, dimension(:), intent(in) :: ppoly_S !< Coefficients of the polynomial reconstruction of T within
1513  !! the layer to be searched.
1514  real :: z !< Position where drho = 0
1515  ! Local variables
1516  real :: dRdT_diff, dRdS_diff
1517  real :: drho, drho_dz, dRdT_z, dRdS_z, T_z, S_z, deltaT, deltaS, deltaP, dT_dz, dS_dz
1518  real :: drho_min, drho_max, ztest, zmin, zmax, dRdT_sum, dRdS_sum, dz, P_z, dP_dz
1519  real :: a1, a2
1520  integer :: iter
1521  integer :: nterm
1522  real :: T_top, T_bot, S_top, S_bot
1523 
1524  nterm = SIZE(ppoly_t)
1525 
1526  ! Position independent quantities
1527  drdt_diff = drdt_bot - drdt_top
1528  drds_diff = drds_bot - drds_top
1529  ! Assume a linear increase in pressure from top and bottom of the cell
1530  dp_dz = p_bot - p_top
1531  ! Initial starting drho (used for bisection)
1532  zmin = z0 ! Lower bounding interval
1533  zmax = 1. ! Maximum bounding interval (bottom of layer)
1534  a1 = 1. - zmin
1535  a2 = zmin
1536  t_z = evaluation_polynomial( ppoly_t, nterm, zmin )
1537  s_z = evaluation_polynomial( ppoly_s, nterm, zmin )
1538  drdt_z = a1*drdt_top + a2*drdt_bot
1539  drds_z = a1*drds_top + a2*drds_bot
1540  p_z = a1*p_top + a2*p_bot
1541  drho_min = delta_rho_from_derivs(t_z, s_z, p_z, drdt_z, drds_z, &
1542  t_ref, s_ref, p_ref, drdt_ref, drds_ref)
1543 
1544  t_z = evaluation_polynomial( ppoly_t, nterm, 1. )
1545  s_z = evaluation_polynomial( ppoly_s, nterm, 1. )
1546  drho_max = delta_rho_from_derivs(t_z, s_z, p_bot, drdt_bot, drds_bot, &
1547  t_ref, s_ref, p_ref, drdt_ref, drds_ref)
1548 
1549  if (drho_min >= 0.) then
1550  z = z0
1551  return
1552  elseif (drho_max == 0.) then
1553  z = 1.
1554  return
1555  endif
1556  if ( sign(1.,drho_min) == sign(1.,drho_max) ) then
1557  call mom_error(fatal, "drho_min is the same sign as dhro_max")
1558  endif
1559 
1560  z = z0
1561  ztest = z0
1562  do iter = 1, cs%max_iter
1563  ! Calculate quantities at the current nondimensional position
1564  a1 = 1.-z
1565  a2 = z
1566  drdt_z = a1*drdt_top + a2*drdt_bot
1567  drds_z = a1*drds_top + a2*drds_bot
1568  t_z = evaluation_polynomial( ppoly_t, nterm, z )
1569  s_z = evaluation_polynomial( ppoly_s, nterm, z )
1570  p_z = a1*p_top + a2*p_bot
1571  deltat = t_z - t_ref
1572  deltas = s_z - s_ref
1573  deltap = p_z - p_ref
1574  drdt_sum = drdt_ref + drdt_z
1575  drds_sum = drds_ref + drds_z
1576  drho = delta_rho_from_derivs(t_z, s_z, p_z, drdt_z, drds_z, &
1577  t_ref, s_ref, p_ref, drdt_ref, drds_ref)
1578 
1579  ! Check for convergence
1580  if (abs(drho) <= cs%drho_tol) exit
1581  ! Update bisection bracketing intervals
1582  if (drho < 0. .and. drho > drho_min) then
1583  drho_min = drho
1584  zmin = z
1585  elseif (drho > 0. .and. drho < drho_max) then
1586  drho_max = drho
1587  zmax = z
1588  endif
1589 
1590  ! Calculate a Newton step
1591  dt_dz = first_derivative_polynomial( ppoly_t, nterm, z )
1592  ds_dz = first_derivative_polynomial( ppoly_s, nterm, z )
1593  drho_dz = 0.5*( (drdt_diff*deltat + drdt_sum*dt_dz) + (drds_diff*deltas + drds_sum*ds_dz) )
1594 
1595  ztest = z - drho/drho_dz
1596  ! Take a bisection if z falls out of [zmin,zmax]
1597  if (ztest < zmin .or. ztest > zmax) then
1598  if ( drho < 0. ) then
1599  ztest = 0.5*(z + zmax)
1600  else
1601  ztest = 0.5*(zmin + z)
1602  endif
1603  endif
1604 
1605  ! Test to ensure we haven't stalled out
1606  if ( abs(z-ztest) <= cs%x_tol ) exit
1607  ! Reset for next iteration
1608  z = ztest
1609  enddo
1610 

References delta_rho_from_derivs(), polynomial_functions::evaluation_polynomial(), polynomial_functions::first_derivative_polynomial(), and mom_error_handler::mom_error().

Referenced by ndiff_unit_tests_discontinuous(), and search_other_column().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ find_neutral_surface_positions_continuous()

subroutine mom_neutral_diffusion::find_neutral_surface_positions_continuous ( integer, intent(in)  nk,
real, dimension(nk+1), intent(in)  Pl,
real, dimension(nk+1), intent(in)  Tl,
real, dimension(nk+1), intent(in)  Sl,
real, dimension(nk+1), intent(in)  dRdTl,
real, dimension(nk+1), intent(in)  dRdSl,
real, dimension(nk+1), intent(in)  Pr,
real, dimension(nk+1), intent(in)  Tr,
real, dimension(nk+1), intent(in)  Sr,
real, dimension(nk+1), intent(in)  dRdTr,
real, dimension(nk+1), intent(in)  dRdSr,
real, dimension(2*nk+2), intent(inout)  PoL,
real, dimension(2*nk+2), intent(inout)  PoR,
integer, dimension(2*nk+2), intent(inout)  KoL,
integer, dimension(2*nk+2), intent(inout)  KoR,
real, dimension(2*nk+1), intent(inout)  hEff,
integer, intent(in), optional  bl_kl,
integer, intent(in), optional  bl_kr,
real, intent(in), optional  bl_zl,
real, intent(in), optional  bl_zr 
)
private

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

Parameters
[in]nkNumber of levels
[in]plLeft-column interface pressure [Pa]
[in]tlLeft-column interface potential temperature [degC]
[in]slLeft-column interface salinity [ppt]
[in]drdtlLeft-column dRho/dT [kg m-3 degC-1]
[in]drdslLeft-column dRho/dS [kg m-3 ppt-1]
[in]prRight-column interface pressure [Pa]
[in]trRight-column interface potential temperature [degC]
[in]srRight-column interface salinity [ppt]
[in]drdtrLeft-column dRho/dT [kg m-3 degC-1]
[in]drdsrLeft-column dRho/dS [kg m-3 ppt-1]
[in,out]polFractional position of neutral surface within layer KoL of left column
[in,out]porFractional position of neutral surface within layer KoR of right column
[in,out]kolIndex of first left interface above neutral surface
[in,out]korIndex of first right interface above neutral surface
[in,out]heffEffective thickness between two neutral surfaces [Pa]
[in]bl_klLayer index of the boundary layer (left)
[in]bl_krLayer index of the boundary layer (right)
[in]bl_zlNondimensional position of the boundary layer (left)
[in]bl_zrNondimensional position of the boundary layer (right)

Definition at line 894 of file MOM_neutral_diffusion.F90.

894  integer, intent(in) :: nk !< Number of levels
895  real, dimension(nk+1), intent(in) :: Pl !< Left-column interface pressure [Pa]
896  real, dimension(nk+1), intent(in) :: Tl !< Left-column interface potential temperature [degC]
897  real, dimension(nk+1), intent(in) :: Sl !< Left-column interface salinity [ppt]
898  real, dimension(nk+1), intent(in) :: dRdTl !< Left-column dRho/dT [kg m-3 degC-1]
899  real, dimension(nk+1), intent(in) :: dRdSl !< Left-column dRho/dS [kg m-3 ppt-1]
900  real, dimension(nk+1), intent(in) :: Pr !< Right-column interface pressure [Pa]
901  real, dimension(nk+1), intent(in) :: Tr !< Right-column interface potential temperature [degC]
902  real, dimension(nk+1), intent(in) :: Sr !< Right-column interface salinity [ppt]
903  real, dimension(nk+1), intent(in) :: dRdTr !< Left-column dRho/dT [kg m-3 degC-1]
904  real, dimension(nk+1), intent(in) :: dRdSr !< Left-column dRho/dS [kg m-3 ppt-1]
905  real, dimension(2*nk+2), intent(inout) :: PoL !< Fractional position of neutral surface within
906  !! layer KoL of left column
907  real, dimension(2*nk+2), intent(inout) :: PoR !< Fractional position of neutral surface within
908  !! layer KoR of right column
909  integer, dimension(2*nk+2), intent(inout) :: KoL !< Index of first left interface above neutral surface
910  integer, dimension(2*nk+2), intent(inout) :: KoR !< Index of first right interface above neutral surface
911  real, dimension(2*nk+1), intent(inout) :: hEff !< Effective thickness between two neutral surfaces [Pa]
912  integer, optional, intent(in) :: bl_kl !< Layer index of the boundary layer (left)
913  integer, optional, intent(in) :: bl_kr !< Layer index of the boundary layer (right)
914  real, optional, intent(in) :: bl_zl !< Nondimensional position of the boundary layer (left)
915  real, optional, intent(in) :: bl_zr !< Nondimensional position of the boundary layer (right)
916 
917  ! Local variables
918  integer :: ns ! Number of neutral surfaces
919  integer :: k_surface ! Index of neutral surface
920  integer :: kl ! Index of left interface
921  integer :: kr ! Index of right interface
922  real :: dRdT, dRdS ! dRho/dT and dRho/dS for the neutral surface
923  logical :: searching_left_column ! True if searching for the position of a right interface in the left column
924  logical :: searching_right_column ! True if searching for the position of a left interface in the right column
925  logical :: reached_bottom ! True if one of the bottom-most interfaces has been used as the target
926  integer :: krm1, klm1
927  real :: dRho, dRhoTop, dRhoBot, hL, hR
928  integer :: lastK_left, lastK_right
929  real :: lastP_left, lastP_right
930  logical :: interior_limit
931 
932  ns = 2*nk+2
933 
934  ! Initialize variables for the search
935  kr = 1 ;
936  kl = 1 ;
937  lastp_right = 0.
938  lastp_left = 0.
939  lastk_right = 1
940  lastk_left = 1
941  reached_bottom = .false.
942 
943  ! Check to see if we should limit the diffusion to the interior
944  interior_limit = PRESENT(bl_kl) .and. PRESENT(bl_kr) .and. PRESENT(bl_zr) .and. PRESENT(bl_zl)
945 
946  ! Loop over each neutral surface, working from top to bottom
947  neutral_surfaces: do k_surface = 1, ns
948  klm1 = max(kl-1, 1)
949  if (klm1>nk) stop 'find_neutral_surface_positions(): klm1 went out of bounds!'
950  krm1 = max(kr-1, 1)
951  if (krm1>nk) stop 'find_neutral_surface_positions(): krm1 went out of bounds!'
952 
953  ! Potential density difference, rho(kr) - rho(kl)
954  drho = 0.5 * ( ( drdtr(kr) + drdtl(kl) ) * ( tr(kr) - tl(kl) ) &
955  + ( drdsr(kr) + drdsl(kl) ) * ( sr(kr) - sl(kl) ) )
956  ! Which column has the lighter surface for the current indexes, kr and kl
957  if (.not. reached_bottom) then
958  if (drho < 0.) then
959  searching_left_column = .true.
960  searching_right_column = .false.
961  elseif (drho > 0.) then
962  searching_right_column = .true.
963  searching_left_column = .false.
964  else ! dRho == 0.
965  if (kl + kr == 2) then ! Still at surface
966  searching_left_column = .true.
967  searching_right_column = .false.
968  else ! Not the surface so we simply change direction
969  searching_left_column = .not. searching_left_column
970  searching_right_column = .not. searching_right_column
971  endif
972  endif
973  endif
974 
975  if (searching_left_column) then
976  ! Interpolate for the neutral surface position within the left column, layer klm1
977  ! Potential density difference, rho(kl-1) - rho(kr) (should be negative)
978  drhotop = 0.5 * ( ( drdtl(klm1) + drdtr(kr) ) * ( tl(klm1) - tr(kr) ) &
979  + ( drdsl(klm1) + drdsr(kr) ) * ( sl(klm1) - sr(kr) ) )
980  ! Potential density difference, rho(kl) - rho(kr) (will be positive)
981  drhobot = 0.5 * ( ( drdtl(klm1+1) + drdtr(kr) ) * ( tl(klm1+1) - tr(kr) ) &
982  + ( drdsl(klm1+1) + drdsr(kr) ) * ( sl(klm1+1) - sr(kr) ) )
983 
984  ! Because we are looking left, the right surface, kr, is lighter than klm1+1 and should be denser than klm1
985  ! unless we are still at the top of the left column (kl=1)
986  if (drhotop > 0. .or. kr+kl==2) then
987  pol(k_surface) = 0. ! The right surface is lighter than anything in layer klm1
988  elseif (drhotop >= drhobot) then ! Left layer is unstratified
989  pol(k_surface) = 1.
990  else
991  ! Linearly interpolate for the position between Pl(kl-1) and Pl(kl) where the density difference
992  ! between right and left is zero.
993  pol(k_surface) = interpolate_for_nondim_position( drhotop, pl(klm1), drhobot, pl(klm1+1) )
994  endif
995  if (pol(k_surface)>=1. .and. klm1<nk) then ! >= is really ==, when PoL==1 we point to the bottom of the cell
996  klm1 = klm1 + 1
997  pol(k_surface) = pol(k_surface) - 1.
998  endif
999  if (real(klm1-lastk_left)+(pol(k_surface)-lastp_left)<0.) then
1000  pol(k_surface) = lastp_left
1001  klm1 = lastk_left
1002  endif
1003  kol(k_surface) = klm1
1004  if (kr <= nk) then
1005  por(k_surface) = 0.
1006  kor(k_surface) = kr
1007  else
1008  por(k_surface) = 1.
1009  kor(k_surface) = nk
1010  endif
1011  if (kr <= nk) then
1012  kr = kr + 1
1013  else
1014  reached_bottom = .true.
1015  searching_right_column = .true.
1016  searching_left_column = .false.
1017  endif
1018  elseif (searching_right_column) then
1019  ! Interpolate for the neutral surface position within the right column, layer krm1
1020  ! Potential density difference, rho(kr-1) - rho(kl) (should be negative)
1021  drhotop = 0.5 * ( ( drdtr(krm1) + drdtl(kl) ) * ( tr(krm1) - tl(kl) ) &
1022  + ( drdsr(krm1) + drdsl(kl) ) * ( sr(krm1) - sl(kl) ) )
1023  ! Potential density difference, rho(kr) - rho(kl) (will be positive)
1024  drhobot = 0.5 * ( ( drdtr(krm1+1) + drdtl(kl) ) * ( tr(krm1+1) - tl(kl) ) &
1025  + ( drdsr(krm1+1) + drdsl(kl) ) * ( sr(krm1+1) - sl(kl) ) )
1026 
1027  ! Because we are looking right, the left surface, kl, is lighter than krm1+1 and should be denser than krm1
1028  ! unless we are still at the top of the right column (kr=1)
1029  if (drhotop >= 0. .or. kr+kl==2) then
1030  por(k_surface) = 0. ! The left surface is lighter than anything in layer krm1
1031  elseif (drhotop >= drhobot) then ! Right layer is unstratified
1032  por(k_surface) = 1.
1033  else
1034  ! Linearly interpolate for the position between Pr(kr-1) and Pr(kr) where the density difference
1035  ! between right and left is zero.
1036  por(k_surface) = interpolate_for_nondim_position( drhotop, pr(krm1), drhobot, pr(krm1+1) )
1037  endif
1038  if (por(k_surface)>=1. .and. krm1<nk) then ! >= is really ==, when PoR==1 we point to the bottom of the cell
1039  krm1 = krm1 + 1
1040  por(k_surface) = por(k_surface) - 1.
1041  endif
1042  if (real(krm1-lastk_right)+(por(k_surface)-lastp_right)<0.) then
1043  por(k_surface) = lastp_right
1044  krm1 = lastk_right
1045  endif
1046  kor(k_surface) = krm1
1047  if (kl <= nk) then
1048  pol(k_surface) = 0.
1049  kol(k_surface) = kl
1050  else
1051  pol(k_surface) = 1.
1052  kol(k_surface) = nk
1053  endif
1054  if (kl <= nk) then
1055  kl = kl + 1
1056  else
1057  reached_bottom = .true.
1058  searching_right_column = .false.
1059  searching_left_column = .true.
1060  endif
1061  else
1062  stop 'Else what?'
1063  endif
1064  if (interior_limit) then
1065  if (kol(k_surface)<=bl_kl) then
1066  kol(k_surface) = bl_kl
1067  if (pol(k_surface)<bl_zl) then
1068  pol(k_surface) = bl_zl
1069  endif
1070  endif
1071  if (kor(k_surface)<=bl_kr) then
1072  kor(k_surface) = bl_kr
1073  if (por(k_surface)<bl_zr) then
1074  por(k_surface) = bl_zr
1075  endif
1076  endif
1077  endif
1078 
1079  lastk_left = kol(k_surface) ; lastp_left = pol(k_surface)
1080  lastk_right = kor(k_surface) ; lastp_right = por(k_surface)
1081  ! Effective thickness
1082  ! NOTE: This would be better expressed in terms of the layers thicknesses rather
1083  ! than as differences of position - AJA
1084  if (k_surface>1) then
1085  hl = absolute_position(nk,ns,pl,kol,pol,k_surface) - absolute_position(nk,ns,pl,kol,pol,k_surface-1)
1086  hr = absolute_position(nk,ns,pr,kor,por,k_surface) - absolute_position(nk,ns,pr,kor,por,k_surface-1)
1087  if ( hl + hr > 0.) then
1088  heff(k_surface-1) = 2. * hl * hr / ( hl + hr ) ! Harmonic mean of layer thicknesses
1089  else
1090  heff(k_surface-1) = 0.
1091  endif
1092  endif
1093 
1094  enddo neutral_surfaces
1095 

References absolute_position(), and interpolate_for_nondim_position().

Referenced by ndiff_unit_tests_continuous(), and neutral_diffusion_calc_coeffs().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ find_neutral_surface_positions_discontinuous()

subroutine mom_neutral_diffusion::find_neutral_surface_positions_discontinuous ( type(neutral_diffusion_cs), intent(inout)  CS,
integer, intent(in)  nk,
real, dimension(nk,2), intent(in)  Pres_l,
real, dimension(nk), intent(in)  hcol_l,
real, dimension(nk,2), intent(in)  Tl,
real, dimension(nk,2), intent(in)  Sl,
real, dimension(:,:), intent(in)  ppoly_T_l,
real, dimension(:,:), intent(in)  ppoly_S_l,
logical, dimension(nk), intent(in)  stable_l,
real, dimension(nk,2), intent(in)  Pres_r,
real, dimension(nk), intent(in)  hcol_r,
real, dimension(nk,2), intent(in)  Tr,
real, dimension(nk,2), intent(in)  Sr,
real, dimension(:,:), intent(in)  ppoly_T_r,
real, dimension(:,:), intent(in)  ppoly_S_r,
logical, dimension(nk), intent(in)  stable_r,
real, dimension(4*nk), intent(inout)  PoL,
real, dimension(4*nk), intent(inout)  PoR,
integer, dimension(4*nk), intent(inout)  KoL,
integer, dimension(4*nk), intent(inout)  KoR,
real, dimension(4*nk-1), intent(inout)  hEff,
real, intent(in), optional  zeta_bot_L,
real, intent(in), optional  zeta_bot_R,
integer, intent(in), optional  k_bot_L,
integer, intent(in), optional  k_bot_R,
logical, intent(in), optional  hard_fail_heff 
)
private

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
[in,out]csNeutral diffusion control structure
[in]nkNumber of levels
[in]pres_lLeft-column interface pressure (Pa)
[in]hcol_lLeft-column layer thicknesses
[in]tlLeft-column top interface potential temperature (degC)
[in]slLeft-column top interface salinity (ppt)
[in]ppoly_t_lLeft-column coefficients of T reconstruction
[in]ppoly_s_lLeft-column coefficients of S reconstruction
[in]stable_lLeft-column, top interface dRho/dS (kg/m3/ppt)
[in]pres_rRight-column interface pressure (Pa)
[in]hcol_rLeft-column layer thicknesses
[in]trRight-column top interface potential temperature (degC)
[in]srRight-column top interface salinity (ppt)
[in]ppoly_t_rRight-column coefficients of T reconstruction
[in]ppoly_s_rRight-column coefficients of S reconstruction
[in]stable_rLeft-column, top interface dRho/dS (kg/m3/ppt)
[in,out]polFractional position of neutral surface within layer KoL of left column
[in,out]porFractional position of neutral surface within layer KoR of right column
[in,out]kolIndex of first left interface above neutral surface
[in,out]korIndex of first right interface above neutral surface
[in,out]heffEffective thickness between two neutral surfaces (Pa)
[in]zeta_bot_lNon-dimensional distance to where the boundary layer intersetcs the cell (left) [nondim]
[in]zeta_bot_rNon-dimensional distance to where the boundary layer intersetcs the cell (right) [nondim]
[in]k_bot_lk-index for the boundary layer (left) [nondim]
[in]k_bot_rk-index for the boundary layer (right) [nondim]
[in]hard_fail_heffIf true (default) bring down the model if the neutral surfaces ever cross [logical]

Definition at line 1141 of file MOM_neutral_diffusion.F90.

1141 
1142  type(neutral_diffusion_CS), intent(inout) :: CS !< Neutral diffusion control structure
1143  integer, intent(in) :: nk !< Number of levels
1144  real, dimension(nk,2), intent(in) :: Pres_l !< Left-column interface pressure (Pa)
1145  real, dimension(nk), intent(in) :: hcol_l !< Left-column layer thicknesses
1146  real, dimension(nk,2), intent(in) :: Tl !< Left-column top interface potential temperature (degC)
1147  real, dimension(nk,2), intent(in) :: Sl !< Left-column top interface salinity (ppt)
1148  real, dimension(:,:), intent(in) :: ppoly_T_l !< Left-column coefficients of T reconstruction
1149  real, dimension(:,:), intent(in) :: ppoly_S_l !< Left-column coefficients of S reconstruction
1150  logical, dimension(nk), intent(in) :: stable_l !< Left-column, top interface dRho/dS (kg/m3/ppt)
1151  real, dimension(nk,2), intent(in) :: Pres_r !< Right-column interface pressure (Pa)
1152  real, dimension(nk), intent(in) :: hcol_r !< Left-column layer thicknesses
1153  real, dimension(nk,2), intent(in) :: Tr !< Right-column top interface potential temperature (degC)
1154  real, dimension(nk,2), intent(in) :: Sr !< Right-column top interface salinity (ppt)
1155  real, dimension(:,:), intent(in) :: ppoly_T_r !< Right-column coefficients of T reconstruction
1156  real, dimension(:,:), intent(in) :: ppoly_S_r !< Right-column coefficients of S reconstruction
1157  logical, dimension(nk), intent(in) :: stable_r !< Left-column, top interface dRho/dS (kg/m3/ppt)
1158  real, dimension(4*nk), intent(inout) :: PoL !< Fractional position of neutral surface within
1159  !! layer KoL of left column
1160  real, dimension(4*nk), intent(inout) :: PoR !< Fractional position of neutral surface within
1161  !! layer KoR of right column
1162  integer, dimension(4*nk), intent(inout) :: KoL !< Index of first left interface above neutral surface
1163  integer, dimension(4*nk), intent(inout) :: KoR !< Index of first right interface above neutral surface
1164  real, dimension(4*nk-1), intent(inout) :: hEff !< Effective thickness between two neutral surfaces (Pa)
1165  real, optional, intent(in) :: zeta_bot_L!< Non-dimensional distance to where the boundary layer
1166  !! intersetcs the cell (left) [nondim]
1167  real, optional, intent(in) :: zeta_bot_R!< Non-dimensional distance to where the boundary layer
1168  !! intersetcs the cell (right) [nondim]
1169 
1170  integer, optional, intent(in) :: k_bot_L !< k-index for the boundary layer (left) [nondim]
1171  integer, optional, intent(in) :: k_bot_R !< k-index for the boundary layer (right) [nondim]
1172  logical, optional, intent(in) :: hard_fail_heff !< If true (default) bring down the model if the
1173  !! neutral surfaces ever cross [logical]
1174  ! Local variables
1175  integer :: ns ! Number of neutral surfaces
1176  integer :: k_surface ! Index of neutral surface
1177  integer :: kl_left, kl_right ! Index of layers on the left/right
1178  integer :: ki_left, ki_right ! Index of interfaces on the left/right
1179  logical :: searching_left_column ! True if searching for the position of a right interface in the left column
1180  logical :: searching_right_column ! True if searching for the position of a left interface in the right column
1181  logical :: reached_bottom ! True if one of the bottom-most interfaces has been used as the target
1182  logical :: search_layer
1183  logical :: fail_heff ! By default,
1184  real :: dRho, dRhoTop, dRhoBot, hL, hR
1185  real :: z0, pos
1186  real :: dRdT_from_top, dRdS_from_top ! Density derivatives at the searched from interface
1187  real :: dRdT_from_bot, dRdS_from_bot ! Density derivatives at the searched from interface
1188  real :: dRdT_to_top, dRdS_to_top ! Density derivatives at the interfaces being searched
1189  real :: dRdT_to_bot, dRdS_to_bot ! Density derivatives at the interfaces being searched
1190  real :: T_ref, S_ref, P_ref, P_top, P_bot
1191  real :: lastP_left, lastP_right
1192  integer :: k_init_L, k_init_R ! Starting indices layers for left and right
1193  real :: p_init_L, p_init_R ! Starting positions for left and right
1194  ! Initialize variables for the search
1195  ns = 4*nk
1196  ki_right = 1
1197  ki_left = 1
1198  kl_left = 1
1199  kl_right = 1
1200  lastp_left = 0.
1201  lastp_right = 0.
1202  reached_bottom = .false.
1203  searching_left_column = .false.
1204  searching_right_column = .false.
1205 
1206  fail_heff = .true.
1207  if (PRESENT(hard_fail_heff)) fail_heff = hard_fail_heff
1208 
1209  if (PRESENT(k_bot_l) .and. PRESENT(k_bot_r) .and. PRESENT(zeta_bot_l) .and. PRESENT(zeta_bot_r)) then
1210  k_init_l = k_bot_l; k_init_r = k_bot_r
1211  p_init_l = zeta_bot_l; p_init_r = zeta_bot_r
1212  lastp_left = zeta_bot_l; lastp_right = zeta_bot_r
1213  kl_left = k_bot_l; kl_right = k_bot_r
1214  else
1215  k_init_l = 1 ; k_init_r = 1
1216  p_init_l = 0. ; p_init_r = 0.
1217  endif
1218  ! Loop over each neutral surface, working from top to bottom
1219  neutral_surfaces: do k_surface = 1, ns
1220 
1221  if (k_surface == ns) then
1222  pol(k_surface) = 1.
1223  por(k_surface) = 1.
1224  kol(k_surface) = nk
1225  kor(k_surface) = nk
1226  ! If the layers are unstable, then simply point the surface to the previous location
1227  elseif (.not. stable_l(kl_left)) then
1228  if (k_surface > 1) then
1229  pol(k_surface) = ki_left - 1 ! Top interface is at position = 0., Bottom is at position = 1
1230  kol(k_surface) = kl_left
1231  por(k_surface) = por(k_surface-1)
1232  kor(k_surface) = kor(k_surface-1)
1233  else
1234  por(k_surface) = p_init_r
1235  kor(k_surface) = k_init_r
1236  pol(k_surface) = p_init_l
1237  kol(k_surface) = k_init_l
1238  endif
1239  call increment_interface(nk, kl_left, ki_left, reached_bottom, searching_left_column, searching_right_column)
1240  searching_left_column = .true.
1241  searching_right_column = .false.
1242  elseif (.not. stable_r(kl_right)) then ! Check the right layer for stability
1243  if (k_surface > 1) then
1244  por(k_surface) = ki_right - 1 ! Top interface is at position = 0., Bottom is at position = 1
1245  kor(k_surface) = kl_right
1246  pol(k_surface) = pol(k_surface-1)
1247  kol(k_surface) = kol(k_surface-1)
1248  else
1249  por(k_surface) = 0.
1250  kor(k_surface) = 1
1251  pol(k_surface) = 0.
1252  kol(k_surface) = 1
1253  endif
1254  call increment_interface(nk, kl_right, ki_right, reached_bottom, searching_right_column, searching_left_column)
1255  searching_left_column = .false.
1256  searching_right_column = .true.
1257  else ! Layers are stable so need to figure out whether we need to search right or left
1258  ! For convenience, the left column uses the searched "from" interface variables, and the right column
1259  ! uses the searched 'to'. These will get reset in subsequent calc_delta_rho calls
1260 
1261  call calc_delta_rho_and_derivs(cs, &
1262  tr(kl_right, ki_right), sr(kl_right, ki_right), pres_r(kl_right,ki_right), &
1263  tl(kl_left, ki_left), sl(kl_left, ki_left) , pres_l(kl_left,ki_left), &
1264  drho)
1265  if (cs%debug) write(*,'(A,I2,A,E12.4,A,I2,A,I2,A,I2,A,I2)') "k_surface=",k_surface," dRho=",drho, &
1266  "kl_left=",kl_left, " ki_left=",ki_left," kl_right=",kl_right, " ki_right=",ki_right
1267  ! Which column has the lighter surface for the current indexes, kr and kl
1268  if (.not. reached_bottom) then
1269  if (drho < 0.) then
1270  searching_left_column = .true.
1271  searching_right_column = .false.
1272  elseif (drho > 0.) then
1273  searching_left_column = .false.
1274  searching_right_column = .true.
1275  else ! dRho == 0.
1276  if ( ( kl_left + kl_right == 2 ) .and. (ki_left + ki_right == 2) ) then ! Still at surface
1277  searching_left_column = .true.
1278  searching_right_column = .false.
1279  else ! Not the surface so we simply change direction
1280  searching_left_column = .not. searching_left_column
1281  searching_right_column = .not. searching_right_column
1282  endif
1283  endif
1284  endif
1285  if (searching_left_column) then
1286  ! Position of the right interface is known and all quantities are fixed
1287  por(k_surface) = ki_right - 1.
1288  kor(k_surface) = kl_right
1289  pol(k_surface) = search_other_column(cs, k_surface, lastp_left, &
1290  tr(kl_right, ki_right), sr(kl_right, ki_right), pres_r(kl_right, ki_right), &
1291  tl(kl_left,1), sl(kl_left,1), pres_l(kl_left,1), &
1292  tl(kl_left,2), sl(kl_left,2), pres_l(kl_left,2), &
1293  ppoly_t_l(kl_left,:), ppoly_s_l(kl_left,:))
1294  kol(k_surface) = kl_left
1295 
1296  if (cs%debug) then
1297  write(*,'(A,I2)') "Searching left layer ", kl_left
1298  write(*,'(A,I2,X,I2)') "Searching from right: ", kl_right, ki_right
1299  write(*,*) "Temp/Salt Reference: ", tr(kl_right,ki_right), sr(kl_right,ki_right)
1300  write(*,*) "Temp/Salt Top L: ", tl(kl_left,1), sl(kl_left,1)
1301  write(*,*) "Temp/Salt Bot L: ", tl(kl_left,2), sl(kl_left,2)
1302  endif
1303  call increment_interface(nk, kl_right, ki_right, reached_bottom, searching_right_column, searching_left_column)
1304  lastp_left = pol(k_surface)
1305  ! If the right layer increments, then we need to reset the last position on the right
1306  if ( kl_right == (kor(k_surface) + 1) ) lastp_right = 0.
1307 
1308  elseif (searching_right_column) then
1309  ! Position of the right interface is known and all quantities are fixed
1310  pol(k_surface) = ki_left - 1.
1311  kol(k_surface) = kl_left
1312  por(k_surface) = search_other_column(cs, k_surface, lastp_right, &
1313  tl(kl_left, ki_left), sl(kl_left, ki_left), pres_l(kl_left, ki_left), &
1314  tr(kl_right,1), sr(kl_right,1), pres_r(kl_right,1), &
1315  tr(kl_right,2), sr(kl_right,2), pres_r(kl_right,2), &
1316  ppoly_t_r(kl_right,:), ppoly_s_r(kl_right,:))
1317  kor(k_surface) = kl_right
1318 
1319  if (cs%debug) then
1320  write(*,'(A,I2)') "Searching right layer ", kl_right
1321  write(*,'(A,I2,X,I2)') "Searching from left: ", kl_left, ki_left
1322  write(*,*) "Temp/Salt Reference: ", tl(kl_left,ki_left), sl(kl_left,ki_left)
1323  write(*,*) "Temp/Salt Top L: ", tr(kl_right,1), sr(kl_right,1)
1324  write(*,*) "Temp/Salt Bot L: ", tr(kl_right,2), sr(kl_right,2)
1325  endif
1326  call increment_interface(nk, kl_left, ki_left, reached_bottom, searching_left_column, searching_right_column)
1327  lastp_right = por(k_surface)
1328  ! If the right layer increments, then we need to reset the last position on the right
1329  if ( kl_left == (kol(k_surface) + 1) ) lastp_left = 0.
1330  else
1331  stop 'Else what?'
1332  endif
1333  if (cs%debug) write(*,'(A,I3,A,ES16.6,A,I2,A,ES16.6)') "KoL:", kol(k_surface), " PoL:", pol(k_surface), &
1334  " KoR:", kor(k_surface), " PoR:", por(k_surface)
1335  endif
1336  ! Effective thickness
1337  if (k_surface>1) then
1338  if ( kol(k_surface) == kol(k_surface-1) .and. kor(k_surface) == kor(k_surface-1) ) then
1339  hl = (pol(k_surface) - pol(k_surface-1))*hcol_l(kol(k_surface))
1340  hr = (por(k_surface) - por(k_surface-1))*hcol_r(kor(k_surface))
1341  if (hl < 0. .or. hr < 0.) then
1342  if (fail_heff) then
1343  call mom_error(fatal,"Negative thicknesses in neutral diffusion")
1344  else
1345  if (searching_left_column) then
1346  pol(k_surface) = pol(k_surface-1)
1347  kol(k_surface) = kol(k_surface-1)
1348  elseif (searching_right_column) then
1349  por(k_surface) = por(k_surface-1)
1350  kor(k_surface) = kor(k_surface-1)
1351  endif
1352  endif
1353  elseif ( hl + hr == 0. ) then
1354  heff(k_surface-1) = 0.
1355  else
1356  heff(k_surface-1) = 2. * ( (hl * hr) / ( hl + hr ) )! Harmonic mean
1357  if ( kol(k_surface) /= kol(k_surface-1) ) then
1358  call mom_error(fatal,"Neutral sublayer spans multiple layers")
1359  endif
1360  if ( kor(k_surface) /= kor(k_surface-1) ) then
1361  call mom_error(fatal,"Neutral sublayer spans multiple layers")
1362  endif
1363  endif
1364  else
1365  heff(k_surface-1) = 0.
1366  endif
1367  endif
1368  enddo neutral_surfaces

References calc_delta_rho_and_derivs(), increment_interface(), mom_error_handler::mom_error(), and search_other_column().

Referenced by ndiff_unit_tests_discontinuous(), and neutral_diffusion_calc_coeffs().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ fv_diff()

real function mom_neutral_diffusion::fv_diff ( real, intent(in)  hkm1,
real, intent(in)  hk,
real, intent(in)  hkp1,
real, intent(in)  Skm1,
real, intent(in)  Sk,
real, intent(in)  Skp1 
)
private

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
[in]hkm1Left cell width
[in]hkCenter cell width
[in]hkp1Right cell width
[in]skm1Left cell average value
[in]skCenter cell average value
[in]skp1Right cell average value

Definition at line 837 of file MOM_neutral_diffusion.F90.

837  real, intent(in) :: hkm1 !< Left cell width
838  real, intent(in) :: hk !< Center cell width
839  real, intent(in) :: hkp1 !< Right cell width
840  real, intent(in) :: Skm1 !< Left cell average value
841  real, intent(in) :: Sk !< Center cell average value
842  real, intent(in) :: Skp1 !< Right cell average value
843 
844  ! Local variables
845  real :: h_sum, hp, hm
846 
847  h_sum = ( hkm1 + hkp1 ) + hk
848  if (h_sum /= 0.) h_sum = 1./ h_sum
849  hm = hkm1 + hk
850  if (hm /= 0.) hm = 1./ hm
851  hp = hkp1 + hk
852  if (hp /= 0.) hp = 1./ hp
853  fv_diff = ( hk * h_sum ) * &
854  ( ( 2. * hkm1 + hk ) * hp * ( skp1 - sk ) &
855  + ( 2. * hkp1 + hk ) * hm * ( sk - skm1 ) )

Referenced by plm_diff(), and test_fv_diff().

Here is the caller graph for this function:

◆ fvlsq_slope()

real function mom_neutral_diffusion::fvlsq_slope ( real, intent(in)  hkm1,
real, intent(in)  hk,
real, intent(in)  hkp1,
real, intent(in)  Skm1,
real, intent(in)  Sk,
real, intent(in)  Skp1 
)
private

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
[in]hkm1Left cell width
[in]hkCenter cell width
[in]hkp1Right cell width
[in]skm1Left cell average value
[in]skCenter cell average value
[in]skp1Right cell average value

Definition at line 863 of file MOM_neutral_diffusion.F90.

863  real, intent(in) :: hkm1 !< Left cell width
864  real, intent(in) :: hk !< Center cell width
865  real, intent(in) :: hkp1 !< Right cell width
866  real, intent(in) :: Skm1 !< Left cell average value
867  real, intent(in) :: Sk !< Center cell average value
868  real, intent(in) :: Skp1 !< Right cell average value
869 
870  ! Local variables
871  real :: xkm1, xkp1
872  real :: h_sum, hx_sum, hxsq_sum, hxy_sum, hy_sum, det
873 
874  xkm1 = -0.5 * ( hk + hkm1 )
875  xkp1 = 0.5 * ( hk + hkp1 )
876  h_sum = ( hkm1 + hkp1 ) + hk
877  hx_sum = hkm1*xkm1 + hkp1*xkp1
878  hxsq_sum = hkm1*(xkm1**2) + hkp1*(xkp1**2)
879  hxy_sum = hkm1*xkm1*skm1 + hkp1*xkp1*skp1
880  hy_sum = ( hkm1*skm1 + hkp1*skp1 ) + hk*sk
881  det = h_sum * hxsq_sum - hx_sum**2
882  if (det /= 0.) then
883  !a = ( hxsq_sum * hy_sum - hx_sum*hxy_sum ) / det ! a would be mean of straight line fit
884  fvlsq_slope = ( h_sum * hxy_sum - hx_sum*hy_sum ) / det ! Gradient of straight line fit
885  else
886  fvlsq_slope = 0. ! Adcroft's reciprocal rule
887  endif

Referenced by plm_diff(), and test_fvlsq_slope().

Here is the caller graph for this function:

◆ increment_interface()

subroutine mom_neutral_diffusion::increment_interface ( integer, intent(in)  nk,
integer, intent(inout)  kl,
integer, intent(inout)  ki,
logical, intent(inout)  reached_bottom,
logical, intent(inout)  searching_this_column,
logical, intent(inout)  searching_other_column 
)
private

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

Parameters
[in]nkNumber of vertical levels
[in,out]klCurrent layer (potentially updated)
[in,out]kiCurrent interface
[in,out]reached_bottomUpdated when kl == nk and ki == 2
[in,out]searching_this_columnUpdated when kl == nk and ki == 2
[in,out]searching_other_columnUpdated when kl == nk and ki == 2

Definition at line 1461 of file MOM_neutral_diffusion.F90.

1461  integer, intent(in ) :: nk !< Number of vertical levels
1462  integer, intent(inout) :: kl !< Current layer (potentially updated)
1463  integer, intent(inout) :: ki !< Current interface
1464  logical, intent(inout) :: reached_bottom !< Updated when kl == nk and ki == 2
1465  logical, intent(inout) :: searching_this_column !< Updated when kl == nk and ki == 2
1466  logical, intent(inout) :: searching_other_column !< Updated when kl == nk and ki == 2
1467  integer :: k
1468 
1469  reached_bottom = .false.
1470  if (ki == 2) then ! At the bottom interface
1471  if ((ki == 2) .and. (kl < nk) ) then ! Not at the bottom so just go to the next layer
1472  kl = kl+1
1473  ki = 1
1474  elseif ((kl == nk) .and. (ki==2)) then
1475  reached_bottom = .true.
1476  searching_this_column = .false.
1477  searching_other_column = .true.
1478  endif
1479  elseif (ki==1) then ! At the top interface
1480  ki = 2 ! Next interface is same layer, but bottom interface
1481  else
1482  call mom_error(fatal,"Unanticipated eventuality in increment_interface")
1483  endif

References mom_error_handler::mom_error().

Referenced by find_neutral_surface_positions_discontinuous().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ interface_scalar()

subroutine mom_neutral_diffusion::interface_scalar ( integer, intent(in)  nk,
real, dimension(nk), intent(in)  h,
real, dimension(nk), intent(in)  S,
real, dimension(nk+1), intent(inout)  Si,
integer, intent(in)  i_method,
real, intent(in)  h_neglect 
)
private

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

Parameters
[in]nkNumber of levels
[in]hLayer thickness [H ~> m or kg m-2]
[in]sLayer scalar (conc, e.g. ppt)
[in,out]siInterface scalar (conc, e.g. ppt)
[in]i_method=1 use average of PLM edges =2 use continuous PPM edge interpolation
[in]h_neglectA negligibly small thickness [H ~> m or kg m-2]

Definition at line 654 of file MOM_neutral_diffusion.F90.

654  integer, intent(in) :: nk !< Number of levels
655  real, dimension(nk), intent(in) :: h !< Layer thickness [H ~> m or kg m-2]
656  real, dimension(nk), intent(in) :: S !< Layer scalar (conc, e.g. ppt)
657  real, dimension(nk+1), intent(inout) :: Si !< Interface scalar (conc, e.g. ppt)
658  integer, intent(in) :: i_method !< =1 use average of PLM edges
659  !! =2 use continuous PPM edge interpolation
660  real, intent(in) :: h_neglect !< A negligibly small thickness [H ~> m or kg m-2]
661  ! Local variables
662  integer :: k, km2, kp1
663  real, dimension(nk) :: diff
664  real :: Sb, Sa
665 
666  call plm_diff(nk, h, s, 2, 1, diff)
667  si(1) = s(1) - 0.5 * diff(1)
668  if (i_method==1) then
669  do k = 2, nk
670  ! Average of the two edge values (will be bounded and,
671  ! when slopes are unlimited, notionally second-order accurate)
672  sa = s(k-1) + 0.5 * diff(k-1) ! Lower edge value of a PLM reconstruction for layer above
673  sb = s(k) - 0.5 * diff(k) ! Upper edge value of a PLM reconstruction for layer below
674  si(k) = 0.5 * ( sa + sb )
675  enddo
676  elseif (i_method==2) then
677  do k = 2, nk
678  ! PPM quasi-fourth order interpolation for edge values following
679  ! equation 1.6 in Colella & Woodward, 1984: JCP 54, 174-201.
680  km2 = max(1, k-2)
681  kp1 = min(nk, k+1)
682  si(k) = ppm_edge(h(km2), h(k-1), h(k), h(kp1), s(k-1), s(k), diff(k-1), diff(k), h_neglect)
683  enddo
684  endif
685  si(nk+1) = s(nk) + 0.5 * diff(nk)
686 

References plm_diff(), and ppm_edge().

Referenced by ndiff_unit_tests_continuous(), neutral_diffusion_calc_coeffs(), and neutral_surface_flux().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ interpolate_for_nondim_position()

real function mom_neutral_diffusion::interpolate_for_nondim_position ( real, intent(in)  dRhoNeg,
real, intent(in)  Pneg,
real, intent(in)  dRhoPos,
real, intent(in)  Ppos 
)
private

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
[in]drhonegNegative density difference
[in]pnegPosition of negative density difference
[in]drhoposPositive density difference
[in]pposPosition of positive density difference

Definition at line 1101 of file MOM_neutral_diffusion.F90.

1101  real, intent(in) :: dRhoNeg !< Negative density difference
1102  real, intent(in) :: Pneg !< Position of negative density difference
1103  real, intent(in) :: dRhoPos !< Positive density difference
1104  real, intent(in) :: Ppos !< Position of positive density difference
1105 
1106  if (ppos<pneg) then
1107  stop 'interpolate_for_nondim_position: Houston, we have a problem! Ppos<Pneg'
1108  elseif (drhoneg>drhopos) then
1109  write(0,*) 'dRhoNeg, Pneg, dRhoPos, Ppos=',drhoneg, pneg, drhopos, ppos
1110  elseif (drhoneg>drhopos) then
1111  stop 'interpolate_for_nondim_position: Houston, we have a problem! dRhoNeg>dRhoPos'
1112  endif
1113  if (ppos<=pneg) then ! Handle vanished or inverted layers
1114  interpolate_for_nondim_position = 0.5
1115  elseif ( drhopos - drhoneg > 0. ) then
1116  interpolate_for_nondim_position = min( 1., max( 0., -drhoneg / ( drhopos - drhoneg ) ) )
1117  elseif ( drhopos - drhoneg == 0) then
1118  if (drhoneg>0.) then
1119  interpolate_for_nondim_position = 0.
1120  elseif (drhoneg<0.) then
1121  interpolate_for_nondim_position = 1.
1122  else ! dRhoPos = dRhoNeg = 0
1123  interpolate_for_nondim_position = 0.5
1124  endif
1125  else ! dRhoPos - dRhoNeg < 0
1126  interpolate_for_nondim_position = 0.5
1127  endif
1128  if ( interpolate_for_nondim_position < 0. ) &
1129  stop 'interpolate_for_nondim_position: Houston, we have a problem! Pint < Pneg'
1130  if ( interpolate_for_nondim_position > 1. ) &
1131  stop 'interpolate_for_nondim_position: Houston, we have a problem! Pint > Ppos'

Referenced by find_neutral_surface_positions_continuous(), search_other_column(), and test_ifndp().

Here is the caller graph for this function:

◆ mark_unstable_cells()

subroutine mom_neutral_diffusion::mark_unstable_cells ( type(neutral_diffusion_cs), intent(inout)  CS,
integer, intent(in)  nk,
real, dimension(nk,2), intent(in)  T,
real, dimension(nk,2), intent(in)  S,
real, dimension(nk,2), intent(in)  P,
logical, dimension(nk), intent(out)  stable_cell 
)
private

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

Parameters
[in,out]csNeutral diffusion control structure
[in]nkNumber of levels in a column
[in]tTemperature at interfaces
[in]sSalinity at interfaces
[in]pPressure at interfaces
[out]stable_cellTrue if this cell is unstably stratified

Definition at line 1373 of file MOM_neutral_diffusion.F90.

1373  type(neutral_diffusion_CS), intent(inout) :: CS !< Neutral diffusion control structure
1374  integer, intent(in) :: nk !< Number of levels in a column
1375  real, dimension(nk,2), intent(in) :: T !< Temperature at interfaces
1376  real, dimension(nk,2), intent(in) :: S !< Salinity at interfaces
1377  real, dimension(nk,2), intent(in) :: P !< Pressure at interfaces
1378  logical, dimension(nk), intent( out) :: stable_cell !< True if this cell is unstably stratified
1379 
1380  integer :: k, first_stable, prev_stable
1381  real :: delta_rho
1382 
1383  do k = 1,nk
1384  call calc_delta_rho_and_derivs( cs, t(k,2), s(k,2), max(p(k,2),cs%ref_pres), &
1385  t(k,1), s(k,1), max(p(k,1),cs%ref_pres), delta_rho )
1386  stable_cell(k) = delta_rho > 0.
1387  enddo

References calc_delta_rho_and_derivs().

Referenced by ndiff_unit_tests_discontinuous(), and neutral_diffusion_calc_coeffs().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ ndiff_unit_tests_continuous()

logical function mom_neutral_diffusion::ndiff_unit_tests_continuous ( logical, intent(in)  verbose)
private

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

Parameters
[in]verboseIf true, write results to stdout

Definition at line 2034 of file MOM_neutral_diffusion.F90.

2034  logical, intent(in) :: verbose !< If true, write results to stdout
2035  ! Local variables
2036  integer, parameter :: nk = 4
2037  real, dimension(nk+1) :: TiL, TiR1, TiR2, TiR4, Tio ! Test interface temperatures
2038  real, dimension(nk) :: TL ! Test layer temperatures
2039  real, dimension(nk+1) :: SiL ! Test interface salinities
2040  real, dimension(nk+1) :: PiL, PiR4 ! Test interface positions
2041  real, dimension(2*nk+2) :: PiLRo, PiRLo ! Test positions
2042  integer, dimension(2*nk+2) :: KoL, KoR ! Test indexes
2043  real, dimension(2*nk+1) :: hEff ! Test positions
2044  real, dimension(2*nk+1) :: Flx ! Test flux
2045  integer :: k
2046  logical :: v
2047  real :: h_neglect, h_neglect_edge
2048 
2049  h_neglect_edge = 1.0e-10 ; h_neglect = 1.0e-30
2050 
2051  v = verbose
2052 
2053  ndiff_unit_tests_continuous = .false. ! Normally return false
2054  write(*,*) '==== MOM_neutral_diffusion: ndiff_unit_tests_continuous ='
2055 
2056  ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. &
2057  test_fv_diff(v,1.,1.,1., 0.,1.,2., 1., 'FV: Straight line on uniform grid')
2058  ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. &
2059  test_fv_diff(v,1.,1.,0., 0.,4.,8., 7., 'FV: Vanished right cell')
2060  ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. &
2061  test_fv_diff(v,0.,1.,1., 0.,4.,8., 7., 'FV: Vanished left cell')
2062  ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. &
2063  test_fv_diff(v,1.,2.,4., 0.,3.,9., 4., 'FV: Stretched grid')
2064  ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. &
2065  test_fv_diff(v,2.,0.,2., 0.,1.,2., 0., 'FV: Vanished middle cell')
2066  ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. &
2067  test_fv_diff(v,0.,1.,0., 0.,1.,2., 2., 'FV: Vanished on both sides')
2068  ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. &
2069  test_fv_diff(v,1.,0.,0., 0.,1.,2., 0., 'FV: Two vanished cell sides')
2070  ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. &
2071  test_fv_diff(v,0.,0.,0., 0.,1.,2., 0., 'FV: All vanished cells')
2072 
2073  ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. &
2074  test_fvlsq_slope(v,1.,1.,1., 0.,1.,2., 1., 'LSQ: Straight line on uniform grid')
2075  ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. &
2076  test_fvlsq_slope(v,1.,1.,0., 0.,1.,2., 1., 'LSQ: Vanished right cell')
2077  ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. &
2078  test_fvlsq_slope(v,0.,1.,1., 0.,1.,2., 1., 'LSQ: Vanished left cell')
2079  ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. &
2080  test_fvlsq_slope(v,1.,2.,4., 0.,3.,9., 2., 'LSQ: Stretched grid')
2081  ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. &
2082  test_fvlsq_slope(v,1.,0.,1., 0.,1.,2., 2., 'LSQ: Vanished middle cell')
2083  ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. &
2084  test_fvlsq_slope(v,0.,1.,0., 0.,1.,2., 0., 'LSQ: Vanished on both sides')
2085  ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. &
2086  test_fvlsq_slope(v,1.,0.,0., 0.,1.,2., 0., 'LSQ: Two vanished cell sides')
2087  ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. &
2088  test_fvlsq_slope(v,0.,0.,0., 0.,1.,2., 0., 'LSQ: All vanished cells')
2089 
2090  call interface_scalar(4, (/10.,10.,10.,10./), (/24.,18.,12.,6./), tio, 1, h_neglect)
2091  !ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. &
2092  ! test_data1d(5, Tio, (/27.,21.,15.,9.,3./), 'Linear profile, interface temperatures')
2093  ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. &
2094  test_data1d(v,5, tio, (/24.,22.5,15.,7.5,6./), 'Linear profile, linear interface temperatures')
2095  call interface_scalar(4, (/10.,10.,10.,10./), (/24.,18.,12.,6./), tio, 2, h_neglect)
2096  ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. &
2097  test_data1d(v,5, tio, (/24.,22.,15.,8.,6./), 'Linear profile, PPM interface temperatures')
2098 
2099  ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. &
2100  test_ifndp(v,-1.0, 0., 1.0, 1.0, 0.5, 'Check mid-point')
2101  ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. &
2102  test_ifndp(v, 0.0, 0., 1.0, 1.0, 0.0, 'Check bottom')
2103  ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. &
2104  test_ifndp(v, 0.1, 0., 1.1, 1.0, 0.0, 'Check below')
2105  ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. &
2106  test_ifndp(v,-1.0, 0., 0.0, 1.0, 1.0, 'Check top')
2107  ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. &
2108  test_ifndp(v,-1.0, 0., -0.1, 1.0, 1.0, 'Check above')
2109  ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. &
2110  test_ifndp(v,-1.0, 0., 3.0, 1.0, 0.25, 'Check 1/4')
2111  ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. &
2112  test_ifndp(v,-3.0, 0., 1.0, 1.0, 0.75, 'Check 3/4')
2113  ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. &
2114  test_ifndp(v, 1.0, 0., 1.0, 1.0, 0.0, 'Check dRho=0 below')
2115  ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. &
2116  test_ifndp(v,-1.0, 0., -1.0, 1.0, 1.0, 'Check dRho=0 above')
2117  ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. &
2118  test_ifndp(v, 0.0, 0., 0.0, 1.0, 0.5, 'Check dRho=0 mid')
2119  ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. &
2120  test_ifndp(v,-2.0, .5, 5.0, 0.5, 0.5, 'Check dP=0')
2121 
2122  ! Identical columns
2123  call find_neutral_surface_positions_continuous(3, &
2124  (/0.,10.,20.,30./), (/22.,18.,14.,10./), (/0.,0.,0.,0./), & ! Left positions, T and S
2125  (/-1.,-1.,-1.,-1./), (/1.,1.,1.,1./), &! Left dRdT and dRdS
2126  (/0.,10.,20.,30./), (/22.,18.,14.,10./), (/0.,0.,0.,0./), & ! Right positions, T and S
2127  (/-1.,-1.,-1.,-1./), (/1.,1.,1.,1./), &! Right dRdT and dRdS
2128  pilro, pirlo, kol, kor, heff)
2129  ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. test_nsp(v, 8, kol, kor, pilro, pirlo, heff, &
2130  (/1,1,2,2,3,3,3,3/), & ! KoL
2131  (/1,1,2,2,3,3,3,3/), & ! KoR
2132  (/0.,0.,0.,0.,0.,0.,1.,1./), & ! pL
2133  (/0.,0.,0.,0.,0.,0.,1.,1./), & ! pR
2134  (/0.,10.,0.,10.,0.,10.,0./), & ! hEff
2135  'Identical columns')
2136  ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. test_data1d(v, 8, &
2137  absolute_positions(3, 8, (/0.,10.,20.,30./), kol, pilro), &
2138  (/0.,0.,10.,10.,20.,20.,30.,30./), '... left positions')
2139  ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. test_data1d(v, 8, &
2140  absolute_positions(3, 8, (/0.,10.,20.,30./), kor, pirlo), &
2141  (/0.,0.,10.,10.,20.,20.,30.,30./), '... right positions')
2142  call neutral_surface_flux(3, 2*3+2, 2, (/10.,10.,10./), (/10.,10.,10./), & ! nk, hL, hR
2143  (/20.,16.,12./), (/20.,16.,12./), & ! Tl, Tr
2144  pilro, pirlo, kol, kor, heff, flx, .true., &
2145  h_neglect, h_neglect_edge=h_neglect_edge)
2146  ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. test_data1d(v, 7, flx, &
2147  (/0.,0.,0.,0.,0.,0.,0./), 'Identical columns, rho flux (=0)')
2148  call neutral_surface_flux(3, 2*3+2, 2, (/10.,10.,10./), (/10.,10.,10./), & ! nk, hL, hR
2149  (/-1.,-1.,-1./), (/1.,1.,1./), & ! Sl, Sr
2150  pilro, pirlo, kol, kor, heff, flx, .true., &
2151  h_neglect, h_neglect_edge=h_neglect_edge)
2152  ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. test_data1d(v, 7, flx, &
2153  (/0.,20.,0.,20.,0.,20.,0./), 'Identical columns, S flux')
2154 
2155  ! Right column slightly cooler than left
2156  call find_neutral_surface_positions_continuous(3, &
2157  (/0.,10.,20.,30./), (/22.,18.,14.,10./), (/0.,0.,0.,0./), & ! Left positions, T and S
2158  (/-1.,-1.,-1.,-1./), (/1.,1.,1.,1./), &! Left dRdT and dRdS
2159  (/0.,10.,20.,30./), (/20.,16.,12.,8./), (/0.,0.,0.,0./), & ! Right positions, T and S
2160  (/-1.,-1.,-1.,-1./), (/1.,1.,1.,1./), &! Right dRdT and dRdS
2161  pilro, pirlo, kol, kor, heff)
2162  ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. test_nsp(v, 8, kol, kor, pilro, pirlo, heff, &
2163  (/1,1,2,2,3,3,3,3/), & ! kL
2164  (/1,1,1,2,2,3,3,3/), & ! kR
2165  (/0.,0.5,0.,0.5,0.,0.5,1.,1./), & ! pL
2166  (/0.,0.,0.5,0.,0.5,0.,0.5,1./), & ! pR
2167  (/0.,5.,5.,5.,5.,5.,0./), & ! hEff
2168  'Right column slightly cooler')
2169  ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. test_data1d(v, 8, &
2170  absolute_positions(3, 8, (/0.,10.,20.,30./), kol, pilro), &
2171  (/0.,5.,10.,15.,20.,25.,30.,30./), '... left positions')
2172  ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. test_data1d(v, 8, &
2173  absolute_positions(3, 8, (/0.,10.,20.,30./), kor, pirlo), &
2174  (/0.,0.,5.,10.,15.,20.,25.,30./), '... right positions')
2175 
2176  ! Right column slightly warmer than left
2177  call find_neutral_surface_positions_continuous(3, &
2178  (/0.,10.,20.,30./), (/22.,18.,14.,10./), (/0.,0.,0.,0./), & ! Left positions, T and S
2179  (/-1.,-1.,-1.,-1./), (/1.,1.,1.,1./), &! Left dRdT and dRdS
2180  (/0.,10.,20.,30./), (/24.,20.,16.,12./), (/0.,0.,0.,0./), & ! Right positions, T and S
2181  (/-1.,-1.,-1.,-1./), (/1.,1.,1.,1./), &! Right dRdT and dRdS
2182  pilro, pirlo, kol, kor, heff)
2183  ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. test_nsp(v, 8, kol, kor, pilro, pirlo, heff, &
2184  (/1,1,1,2,2,3,3,3/), & ! kL
2185  (/1,1,2,2,3,3,3,3/), & ! kR
2186  (/0.,0.,0.5,0.,0.5,0.,0.5,1./), & ! pL
2187  (/0.,0.5,0.,0.5,0.,0.5,1.,1./), & ! pR
2188  (/0.,5.,5.,5.,5.,5.,0./), & ! hEff
2189  'Right column slightly warmer')
2190 
2191  ! Right column somewhat cooler than left
2192  call find_neutral_surface_positions_continuous(3, &
2193  (/0.,10.,20.,30./), (/22.,18.,14.,10./), (/0.,0.,0.,0./), & ! Left positions, T and S
2194  (/-1.,-1.,-1.,-1./), (/1.,1.,1.,1./), &! Left dRdT and dRdS
2195  (/0.,10.,20.,30./), (/16.,12.,8.,4./), (/0.,0.,0.,0./), & ! Right positions, T and S
2196  (/-1.,-1.,-1.,-1./), (/1.,1.,1.,1./), &! Right dRdT and dRdS
2197  pilro, pirlo, kol, kor, heff)
2198  ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. test_nsp(v, 8, kol, kor, pilro, pirlo, heff, &
2199  (/1,2,2,3,3,3,3,3/), & ! kL
2200  (/1,1,1,1,2,2,3,3/), & ! kR
2201  (/0.,0.,0.5,0.,0.5,1.,1.,1./), & ! pL
2202  (/0.,0.,0.,0.5,0.,0.5,0.,1./), & ! pR
2203  (/0.,0.,5.,5.,5.,0.,0./), & ! hEff
2204  'Right column somewhat cooler')
2205 
2206  ! Right column much colder than left with no overlap
2207  call find_neutral_surface_positions_continuous(3, &
2208  (/0.,10.,20.,30./), (/22.,18.,14.,10./), (/0.,0.,0.,0./), & ! Left positions, T and S
2209  (/-1.,-1.,-1.,-1./), (/1.,1.,1.,1./), &! Left dRdT and dRdS
2210  (/0.,10.,20.,30./), (/9.,7.,5.,3./), (/0.,0.,0.,0./), & ! Right positions, T and S
2211  (/-1.,-1.,-1.,-1./), (/1.,1.,1.,1./), &! Right dRdT and dRdS
2212  pilro, pirlo, kol, kor, heff)
2213  ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. test_nsp(v, 8, kol, kor, pilro, pirlo, heff, &
2214  (/1,2,3,3,3,3,3,3/), & ! kL
2215  (/1,1,1,1,1,2,3,3/), & ! kR
2216  (/0.,0.,0.,1.,1.,1.,1.,1./), & ! pL
2217  (/0.,0.,0.,0.,0.,0.,0.,1./), & ! pR
2218  (/0.,0.,0.,0.,0.,0.,0./), & ! hEff
2219  'Right column much cooler')
2220 
2221  ! Right column with mixed layer
2222  call find_neutral_surface_positions_continuous(3, &
2223  (/0.,10.,20.,30./), (/22.,18.,14.,10./), (/0.,0.,0.,0./), & ! Left positions, T and S
2224  (/-1.,-1.,-1.,-1./), (/1.,1.,1.,1./), &! Left dRdT and dRdS
2225  (/0.,10.,20.,30./), (/14.,14.,10.,2./), (/0.,0.,0.,0./), & ! Right positions, T and S
2226  (/-1.,-1.,-1.,-1./), (/1.,1.,1.,1./), &! Right dRdT and dRdS
2227  pilro, pirlo, kol, kor, heff)
2228  ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. test_nsp(v, 8, kol, kor, pilro, pirlo, heff, &
2229  (/1,2,3,3,3,3,3,3/), & ! kL
2230  (/1,1,1,1,2,3,3,3/), & ! kR
2231  (/0.,0.,0.,0.,0.,1.,1.,1./), & ! pL
2232  (/0.,0.,0.,0.,0.,0.,0.,1./), & ! pR
2233  (/0.,0.,0.,0.,10.,0.,0./), & ! hEff
2234  'Right column with mixed layer')
2235 
2236  ! Identical columns with mixed layer
2237  call find_neutral_surface_positions_continuous(3, &
2238  (/0.,10.,20.,30./), (/14.,14.,10.,2./), (/0.,0.,0.,0./), & ! Left positions, T and S
2239  (/-1.,-1.,-1.,-1./), (/1.,1.,1.,1./), &! Left dRdT and dRdS
2240  (/0.,10.,20.,30./), (/14.,14.,10.,2./), (/0.,0.,0.,0./), & ! Right positions, T and S
2241  (/-1.,-1.,-1.,-1./), (/1.,1.,1.,1./), &! Right dRdT and dRdS
2242  pilro, pirlo, kol, kor, heff)
2243  ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. test_nsp(v, 8, kol, kor, pilro, pirlo, heff, &
2244  (/1,1,2,2,3,3,3,3/), & ! kL
2245  (/1,1,2,2,3,3,3,3/), & ! kR
2246  (/0.,0.,0.,0.,0.,0.,1.,1./), & ! pL
2247  (/0.,0.,0.,0.,0.,0.,1.,1./), & ! pR
2248  (/0.,10.,0.,10.,0.,10.,0./), & ! hEff
2249  'Identical columns with mixed layer')
2250 
2251  ! Right column with unstable mixed layer
2252  call find_neutral_surface_positions_continuous(3, &
2253  (/0.,10.,20.,30./), (/14.,14.,10.,2./), (/0.,0.,0.,0./), & ! Left positions, T and S
2254  (/-1.,-1.,-1.,-1./), (/1.,1.,1.,1./), &! Left dRdT and dRdS
2255  (/0.,10.,20.,30./), (/10.,14.,12.,4./), (/0.,0.,0.,0./), & ! Right positions, T and S
2256  (/-1.,-1.,-1.,-1./), (/1.,1.,1.,1./), &! Right dRdT and dRdS
2257  pilro, pirlo, kol, kor, heff)
2258  ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. test_nsp(v, 8, kol, kor, pilro, pirlo, heff, &
2259  (/1,2,3,3,3,3,3,3/), & ! kL
2260  (/1,1,1,2,3,3,3,3/), & ! kR
2261  (/0.,0.,0.,0.,0.,0.,.75,1./), & ! pL
2262  (/0.,0.,0.,0.,0.,0.25,1.,1./), & ! pR
2263  (/0.,0.,0.,0.,0.,7.5,0./), & ! hEff
2264  'Right column with unstable mixed layer')
2265 
2266  ! Left column with unstable mixed layer
2267  call find_neutral_surface_positions_continuous(3, &
2268  (/0.,10.,20.,30./), (/10.,14.,12.,4./), (/0.,0.,0.,0./), & ! Left positions, T and S
2269  (/-1.,-1.,-1.,-1./), (/1.,1.,1.,1./), &! Left dRdT and dRdS
2270  (/0.,10.,20.,30./), (/14.,14.,10.,2./), (/0.,0.,0.,0./), & ! Right positions, T and S
2271  (/-1.,-1.,-1.,-1./), (/1.,1.,1.,1./), &! Right dRdT and dRdS
2272  pilro, pirlo, kol, kor, heff)
2273  ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. test_nsp(v, 8, kol, kor, pilro, pirlo, heff, &
2274  (/1,1,1,2,3,3,3,3/), & ! kL
2275  (/1,2,3,3,3,3,3,3/), & ! kR
2276  (/0.,0.,0.,0.,0.,0.25,1.,1./), & ! pL
2277  (/0.,0.,0.,0.,0.,0.,.75,1./), & ! pR
2278  (/0.,0.,0.,0.,0.,7.5,0./), & ! hEff
2279  'Left column with unstable mixed layer')
2280 
2281  ! Two unstable mixed layers
2282  call find_neutral_surface_positions_continuous(3, &
2283  (/0.,10.,20.,30./), (/8.,12.,10.,2./), (/0.,0.,0.,0./), & ! Left positions, T and S
2284  (/-1.,-1.,-1.,-1./), (/1.,1.,1.,1./), &! Left dRdT and dRdS
2285  (/0.,10.,20.,30./), (/10.,14.,12.,4./), (/0.,0.,0.,0./), & ! Right positions, T and S
2286  (/-1.,-1.,-1.,-1./), (/1.,1.,1.,1./), &! Right dRdT and dRdS
2287  pilro, pirlo, kol, kor, heff)
2288  ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. test_nsp(v, 8, kol, kor, pilro, pirlo, heff, &
2289  (/1,1,1,1,2,3,3,3/), & ! kL
2290  (/1,2,3,3,3,3,3,3/), & ! kR
2291  (/0.,0.,0.,0.,0.,0.,0.75,1./), & ! pL
2292  (/0.,0.,0.,0.5,0.5,0.5,1.,1./), & ! pR
2293  (/0.,0.,0.,0.,0.,6.,0./), & ! hEff
2294  'Two unstable mixed layers')
2295 
2296  if (.not. ndiff_unit_tests_continuous) write(*,*) 'Pass'
2297 

References absolute_positions(), find_neutral_surface_positions_continuous(), interface_scalar(), neutral_surface_flux(), test_data1d(), test_fv_diff(), test_fvlsq_slope(), test_ifndp(), and test_nsp().

Referenced by neutral_diffusion_unit_tests().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ ndiff_unit_tests_discontinuous()

logical function mom_neutral_diffusion::ndiff_unit_tests_discontinuous ( logical, intent(in)  verbose)
private
Parameters
[in]verboseIt true, write results to stdout

Definition at line 2301 of file MOM_neutral_diffusion.F90.

2301  logical, intent(in) :: verbose !< It true, write results to stdout
2302  ! Local variables
2303  integer, parameter :: nk = 3
2304  integer, parameter :: ns = nk*4
2305  real, dimension(nk) :: Sl, Sr, Tl, Tr, hl, hr
2306  real, dimension(nk,2) :: TiL, SiL, TiR, SiR
2307  real, dimension(nk,2) :: Pres_l, Pres_r
2308  integer, dimension(ns) :: KoL, KoR
2309  real, dimension(ns) :: PoL, PoR
2310  real, dimension(ns-1) :: hEff, Flx
2311  type(neutral_diffusion_CS) :: CS !< Neutral diffusion control structure
2312  type(EOS_type), pointer :: EOS !< Structure for linear equation of state
2313  type(remapping_CS), pointer :: remap_CS !< Remapping control structure (PLM)
2314  real, dimension(nk,2) :: ppoly_T_l, ppoly_T_r ! Linear reconstruction for T
2315  real, dimension(nk,2) :: ppoly_S_l, ppoly_S_r ! Linear reconstruction for S
2316  real, dimension(nk,2) :: dRdT, dRdS
2317  logical, dimension(nk) :: stable_l, stable_r
2318  integer :: iMethod
2319  integer :: ns_l, ns_r
2320  real :: h_neglect, h_neglect_edge
2321  integer :: k
2322  logical :: v
2323 
2324  v = verbose
2325  ndiff_unit_tests_discontinuous = .false. ! Normally return false
2326  write(*,*) '==== MOM_neutral_diffusion: ndiff_unit_tests_discontinuous ='
2327 !
2328  h_neglect = 1.0e-30 ; h_neglect_edge = 1.0e-10
2329 
2330  ! Unit tests for find_neutral_surface_positions_discontinuous
2331  ! Salinity is 0 for all these tests
2332  allocate(cs%EOS)
2333  call eos_manual_init(cs%EOS, form_of_eos = eos_linear, drho_dt = -1., drho_ds = 0.)
2334  sl(:) = 0. ; sr(:) = 0. ; ; sil(:,:) = 0. ; sir(:,:) = 0.
2335  ppoly_t_l(:,:) = 0.; ppoly_t_r(:,:) = 0.
2336  ppoly_s_l(:,:) = 0.; ppoly_s_r(:,:) = 0.
2337  ! Intialize any control structures needed for unit tests
2338  cs%ref_pres = -1.
2339 
2340  hl = (/10.,10.,10./) ; hr = (/10.,10.,10./)
2341  pres_l(1,1) = 0. ; pres_l(1,2) = hl(1) ; pres_r(1,1) = 0. ; pres_r(1,2) = hr(1)
2342  do k = 2,nk
2343  pres_l(k,1) = pres_l(k-1,2)
2344  pres_l(k,2) = pres_l(k,1) + hl(k)
2345  pres_r(k,1) = pres_r(k-1,2)
2346  pres_r(k,2) = pres_r(k,1) + hr(k)
2347  enddo
2348  cs%delta_rho_form = 'mid_pressure'
2349  cs%neutral_pos_method = 1
2350 
2351  til(1,:) = (/ 22.00, 18.00 /); til(2,:) = (/ 18.00, 14.00 /); til(3,:) = (/ 14.00, 10.00 /);
2352  tir(1,:) = (/ 22.00, 18.00 /); tir(2,:) = (/ 18.00, 14.00 /); tir(3,:) = (/ 14.00, 10.00 /);
2353  call mark_unstable_cells( cs, nk, til, sil, pres_l, stable_l )
2354  call mark_unstable_cells( cs, nk, tir, sir, pres_r, stable_r )
2355  call find_neutral_surface_positions_discontinuous(cs, nk, pres_l, hl, til, sil, ppoly_t_l, ppoly_s_l, stable_l, &
2356  pres_r, hr, tir, sir, ppoly_t_r, ppoly_s_r, stable_r, pol, por, kol, kor, heff)
2357  ndiff_unit_tests_discontinuous = ndiff_unit_tests_discontinuous .or. test_nsp(v, 12, kol, kor, pol, por, heff, &
2358  (/ 1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3 /), & ! KoL
2359  (/ 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 3 /), & ! KoR
2360  (/ 0.00, 0.00, 1.00, 1.00, 0.00, 0.00, 1.00, 1.00, 0.00, 0.00, 1.00, 1.00 /), & ! PoL
2361  (/ 0.00, 0.00, 1.00, 0.00, 0.00, 0.00, 1.00, 0.00, 0.00, 0.00, 1.00, 1.00 /), & ! PoR
2362  (/ 0.00, 10.00, 0.00, 0.00, 0.00, 10.00, 0.00, 0.00, 0.00, 10.00, 0.00 /), & ! hEff
2363  'Identical Columns')
2364 
2365  til(1,:) = (/ 22.00, 18.00 /); til(2,:) = (/ 18.00, 14.00 /); til(3,:) = (/ 14.00, 10.00 /);
2366  tir(1,:) = (/ 20.00, 16.00 /); tir(2,:) = (/ 16.00, 12.00 /); tir(3,:) = (/ 12.00, 8.00 /);
2367  call mark_unstable_cells( cs, nk, til, sil, pres_l, stable_l )
2368  call mark_unstable_cells( cs, nk, tir, sir, pres_r, stable_r )
2369  call find_neutral_surface_positions_discontinuous(cs, nk, pres_l, hl, til, sil, ppoly_t_l, ppoly_s_l, stable_l, &
2370  pres_r, hr, tir, sir, ppoly_t_r, ppoly_s_r, stable_r, pol, por, kol, kor, heff)
2371  ndiff_unit_tests_discontinuous = ndiff_unit_tests_discontinuous .or. test_nsp(v, 12, kol, kor, pol, por, heff, &
2372  (/ 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 3 /), & ! KoL
2373  (/ 1, 1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3 /), & ! KoR
2374  (/ 0.00, 0.50, 1.00, 0.00, 0.50, 0.50, 1.00, 0.00, 0.50, 0.50, 1.00, 1.00 /), & ! PoL
2375  (/ 0.00, 0.00, 0.50, 0.50, 1.00, 0.00, 0.50, 0.50, 1.00, 0.00, 0.50, 1.00 /), & ! PoR
2376  (/ 0.00, 5.00, 0.00, 5.00, 0.00, 5.00, 0.00, 5.00, 0.00, 5.00, 0.00 /), & ! hEff
2377  'Right slightly cooler')
2378 
2379  til(1,:) = (/ 20.00, 16.00 /); til(2,:) = (/ 16.00, 12.00 /); til(3,:) = (/ 12.00, 8.00 /);
2380  tir(1,:) = (/ 22.00, 18.00 /); tir(2,:) = (/ 18.00, 14.00 /); tir(3,:) = (/ 14.00, 10.00 /);
2381  call mark_unstable_cells( cs, nk, til, sil, pres_l, stable_l )
2382  call mark_unstable_cells( cs, nk, tir, sir, pres_r, stable_r )
2383  call find_neutral_surface_positions_discontinuous(cs, nk, pres_l, hl, til, sil, ppoly_t_l, ppoly_s_l, stable_l, &
2384  pres_r, hr, tir, sir, ppoly_t_r, ppoly_s_r, stable_r, pol, por, kol, kor, heff)
2385  ndiff_unit_tests_discontinuous = ndiff_unit_tests_discontinuous .or. test_nsp(v, 12, kol, kor, pol, por, heff, &
2386  (/ 1, 1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3 /), & ! KoL
2387  (/ 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 3 /), & ! KoR
2388  (/ 0.00, 0.00, 0.50, 0.50, 1.00, 0.00, 0.50, 0.50, 1.00, 0.00, 0.50, 1.00 /), & ! PoL
2389  (/ 0.00, 0.50, 1.00, 0.00, 0.50, 0.50, 1.00, 0.00, 0.50, 0.50, 1.00, 1.00 /), & ! PoR
2390  (/ 0.00, 5.00, 0.00, 5.00, 0.00, 5.00, 0.00, 5.00, 0.00, 5.00, 0.00 /), & ! hEff
2391  'Left slightly cooler')
2392 
2393  til(1,:) = (/ 22.00, 20.00 /); til(2,:) = (/ 18.00, 16.00 /); til(3,:) = (/ 14.00, 12.00 /);
2394  tir(1,:) = (/ 32.00, 24.00 /); tir(2,:) = (/ 22.00, 14.00 /); tir(3,:) = (/ 12.00, 4.00 /);
2395  call mark_unstable_cells( cs, nk, til, sil, pres_l, stable_l )
2396  call mark_unstable_cells( cs, nk, tir, sir, pres_r, stable_r )
2397  call find_neutral_surface_positions_discontinuous(cs, nk, pres_l, hl, til, sil, ppoly_t_l, ppoly_s_l, stable_l, &
2398  pres_r, hr, tir, sir, ppoly_t_r, ppoly_s_r, stable_r, pol, por, kol, kor, heff)
2399  ndiff_unit_tests_discontinuous = ndiff_unit_tests_discontinuous .or. test_nsp(v, 12, kol, kor, pol, por, heff, &
2400  (/ 1, 1, 1, 1, 1, 2, 2, 3, 3, 3, 3, 3 /), & ! KoL
2401  (/ 1, 1, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3 /), & ! KoR
2402  (/ 0.00, 0.00, 0.00, 0.00, 1.00, 0.00, 1.00, 0.00, 0.00, 1.00, 1.00, 1.00 /), & ! PoL
2403  (/ 0.00, 1.00, 0.00, 0.00, 0.25, 0.50, 0.75, 1.00, 0.00, 0.00, 0.00, 1.00 /), & ! PoR
2404  (/ 0.00, 0.00, 0.00, 4.00, 0.00, 4.00, 0.00, 0.00, 0.00, 0.00, 0.00 /), & ! hEff
2405  'Right more strongly stratified')
2406 
2407  til(1,:) = (/ 22.00, 18.00 /); til(2,:) = (/ 18.00, 14.00 /); til(3,:) = (/ 14.00, 10.00 /);
2408  tir(1,:) = (/ 14.00, 14.00 /); tir(2,:) = (/ 14.00, 14.00 /); tir(3,:) = (/ 12.00, 8.00 /);
2409  call mark_unstable_cells( cs, nk, til, sil, pres_l, stable_l )
2410  call mark_unstable_cells( cs, nk, tir, sir, pres_r, stable_r )
2411  call find_neutral_surface_positions_discontinuous(cs, nk, pres_l, hl, til, sil, ppoly_t_l, ppoly_s_l, stable_l, &
2412  pres_r, hr, tir, sir, ppoly_t_r, ppoly_s_r, stable_r, pol, por, kol, kor, heff)
2413  ndiff_unit_tests_discontinuous = ndiff_unit_tests_discontinuous .or. test_nsp(v, 12, kol, kor, pol, por, heff, &
2414  (/ 1, 1, 1, 1, 1, 1, 2, 2, 3, 3, 3, 3 /), & ! KoL
2415  (/ 1, 1, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3 /), & ! KoR
2416  (/ 0.00, 0.00, 0.00, 0.00, 0.00, 1.00, 0.00, 1.00, 0.00, 0.50, 1.00, 1.00 /), & ! PoL
2417  (/ 0.00, 1.00, 0.00, 1.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.50, 1.00 /), & ! PoR
2418  (/ 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 5.00, 0.00 /), & ! hEff
2419  'Deep Mixed layer on the right')
2420 
2421  til(1,:) = (/ 14.00, 14.00 /); til(2,:) = (/ 14.00, 12.00 /); til(3,:) = (/ 10.00, 8.00 /);
2422  tir(1,:) = (/ 14.00, 14.00 /); tir(2,:) = (/ 14.00, 14.00 /); tir(3,:) = (/ 14.00, 14.00 /);
2423  call mark_unstable_cells( cs, nk, til, sil, pres_l, stable_l )
2424  call mark_unstable_cells( cs, nk, tir, sir, pres_r, stable_r )
2425  call find_neutral_surface_positions_discontinuous(cs, nk, pres_l, hl, til, sil, ppoly_t_l, ppoly_s_l, stable_l, &
2426  pres_r, hr, tir, sir, ppoly_t_r, ppoly_s_r, stable_r, pol, por, kol, kor, heff)
2427  ndiff_unit_tests_discontinuous = ndiff_unit_tests_discontinuous .or. test_nsp(v, 12, kol, kor, pol, por, heff, &
2428  (/ 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 3 /), & ! KoL
2429  (/ 1, 1, 1, 1, 2, 2, 3, 3, 3, 3, 3, 3 /), & ! KoR
2430  (/ 0.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00 /), & ! PoL
2431  (/ 0.00, 0.00, 0.00, 1.00, 0.00, 1.00, 0.00, 1.00, 1.00, 1.00, 1.00, 1.00 /), & ! PoR
2432  (/ 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00 /), & ! hEff
2433  'Right unstratified column')
2434 
2435  til(1,:) = (/ 14.00, 14.00 /); til(2,:) = (/ 14.00, 12.00 /); til(3,:) = (/ 10.00, 8.00 /);
2436  tir(1,:) = (/ 14.00, 14.00 /); tir(2,:) = (/ 14.00, 14.00 /); tir(3,:) = (/ 12.00, 4.00 /);
2437  call mark_unstable_cells( cs, nk, til, sil, pres_l, stable_l )
2438  call mark_unstable_cells( cs, nk, tir, sir, pres_r, stable_r )
2439  call find_neutral_surface_positions_discontinuous(cs, nk, pres_l, hl, til, sil, ppoly_t_l, ppoly_s_l, stable_l, &
2440  pres_r, hr, tir, sir, ppoly_t_r, ppoly_s_r, stable_r, pol, por, kol, kor, heff)
2441  ndiff_unit_tests_discontinuous = ndiff_unit_tests_discontinuous .or. test_nsp(v, 12, kol, kor, pol, por, heff, &
2442  (/ 1, 1, 1, 1, 1, 1, 2, 2, 2, 3, 3, 3 /), & ! KoL
2443  (/ 1, 1, 1, 1, 2, 2, 3, 3, 3, 3, 3, 3 /), & ! KoR
2444  (/ 0.00, 1.00, 1.00, 1.00, 1.00, 1.00, 0.00, 1.00, 1.00, 0.00, 1.00, 1.00 /), & ! PoL
2445  (/ 0.00, 0.00, 0.00, 1.00, 0.00, 1.00, 0.00, 0.00, 0.00, 0.25, 0.50, 1.00 /), & ! PoR
2446  (/ 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 4.00, 0.00 /), & ! hEff
2447  'Right unstratified column')
2448 
2449  til(1,:) = (/ 14.00, 14.00 /); til(2,:) = (/ 14.00, 10.00 /); til(3,:) = (/ 10.00, 2.00 /);
2450  tir(1,:) = (/ 14.00, 14.00 /); tir(2,:) = (/ 14.00, 10.00 /); tir(3,:) = (/ 10.00, 2.00 /);
2451  call mark_unstable_cells( cs, nk, til, sil, pres_l, stable_l )
2452  call mark_unstable_cells( cs, nk, tir, sir, pres_r, stable_r )
2453  call find_neutral_surface_positions_discontinuous(cs, nk, pres_l, hl, til, sil, ppoly_t_l, ppoly_s_l, stable_l, &
2454  pres_r, hr, tir, sir, ppoly_t_r, ppoly_s_r, stable_r, pol, por, kol, kor, heff)
2455  ndiff_unit_tests_discontinuous = ndiff_unit_tests_discontinuous .or. test_nsp(v, 12, kol, kor, pol, por, heff, &
2456  (/ 1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3 /), & ! KoL
2457  (/ 1, 1, 1, 1, 2, 2, 2, 3, 3, 3, 3, 3 /), & ! KoR
2458  (/ 0.00, 1.00, 1.00, 1.00, 0.00, 0.00, 1.00, 1.00, 0.00, 0.00, 1.00, 1.00 /), & ! PoL
2459  (/ 0.00, 0.00, 0.00, 1.00, 0.00, 0.00, 1.00, 0.00, 0.00, 0.00, 1.00, 1.00 /), & ! PoR
2460  (/ 0.00, 0.00, 0.00, 0.00, 0.00, 10.00, 0.00, 0.00, 0.00, 10.00, 0.00 /), & ! hEff
2461  'Identical columns with mixed layer')
2462 
2463  til(1,:) = (/ 14.00, 12.00 /); til(2,:) = (/ 10.00, 10.00 /); til(3,:) = (/ 8.00, 2.00 /);
2464  tir(1,:) = (/ 14.00, 12.00 /); tir(2,:) = (/ 12.00, 8.00 /); tir(3,:) = (/ 8.00, 2.00 /);
2465  call mark_unstable_cells( cs, nk, til, sil, pres_l, stable_l )
2466  call mark_unstable_cells( cs, nk, tir, sir, pres_r, stable_r )
2467  call find_neutral_surface_positions_discontinuous(cs, nk, pres_l, hl, til, sil, ppoly_t_l, ppoly_s_l, stable_l, &
2468  pres_r, hr, tir, sir, ppoly_t_r, ppoly_s_r, stable_r, pol, por, kol, kor, heff)
2469  ndiff_unit_tests_discontinuous = ndiff_unit_tests_discontinuous .or. test_nsp(v, 12, kol, kor, pol, por, heff, &
2470  (/ 1, 1, 1, 1, 2, 2, 3, 3, 3, 3, 3, 3 /), & ! KoL
2471  (/ 1, 1, 1, 2, 2, 2, 2, 2, 2, 3, 3, 3 /), & ! KoR
2472  (/ 0.00, 0.00, 1.00, 1.00, 0.00, 1.00, 0.00, 0.00, 0.00, 0.00, 1.00, 1.00 /), & ! PoL
2473  (/ 0.00, 0.00, 1.00, 0.00, 0.00, 0.00, 0.00, 1.00, 1.00, 0.00, 1.00, 1.00 /), & ! PoR
2474  (/ 0.00, 10.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 10.00, 0.00 /), & ! hEff
2475  'Left interior unstratified')
2476 
2477  til(1,:) = (/ 12.00, 12.00 /); til(2,:) = (/ 12.00, 10.00 /); til(3,:) = (/ 10.00, 6.00 /);
2478  tir(1,:) = (/ 12.00, 10.00 /); tir(2,:) = (/ 10.00, 12.00 /); tir(3,:) = (/ 8.00, 4.00 /);
2479  call mark_unstable_cells( cs, nk, til, sil, pres_l, stable_l )
2480  call mark_unstable_cells( cs, nk, tir, sir, pres_r, stable_r )
2481  call find_neutral_surface_positions_discontinuous(cs, nk, pres_l, hl, til, sil, ppoly_t_l, ppoly_s_l, stable_l, &
2482  pres_r, hr, tir, sir, ppoly_t_r, ppoly_s_r, stable_r, pol, por, kol, kor, heff)
2483  ndiff_unit_tests_discontinuous = ndiff_unit_tests_discontinuous .or. test_nsp(v, 12, kol, kor, pol, por, heff, &
2484  (/ 1, 1, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3 /), & ! KoL
2485  (/ 1, 1, 1, 1, 1, 1, 2, 2, 3, 3, 3, 3 /), & ! KoR
2486  (/ 0.00, 1.00, 0.00, 0.00, 1.00, 0.00, 0.00, 0.00, 0.00, 0.50, 1.00, 1.00 /), & ! PoL
2487  (/ 0.00, 0.00, 0.00, 0.00, 1.00, 1.00, 0.00, 1.00, 0.00, 0.00, 0.50, 1.00 /), & ! PoR
2488  (/ 0.00, 0.00, 0.00, 10.00, 0.00, 0.00, 0.00, 0.00, 0.00, 5.00, 0.00 /), & ! hEff
2489  'Left mixed layer, Right unstable interior')
2490 
2491  til(1,:) = (/ 14.00, 14.00 /); til(2,:) = (/ 10.00, 10.00 /); til(3,:) = (/ 8.00, 6.00 /);
2492  tir(1,:) = (/ 10.00, 14.00 /); tir(2,:) = (/ 16.00, 16.00 /); tir(3,:) = (/ 12.00, 4.00 /);
2493  call mark_unstable_cells( cs, nk, til, sil, pres_l, stable_l )
2494  call mark_unstable_cells( cs, nk, tir, sir, pres_r, stable_r )
2495  call find_neutral_surface_positions_discontinuous(cs, nk, pres_l, hl, til, sil, ppoly_t_l, ppoly_s_l, stable_l, &
2496  pres_r, hr, tir, sir, ppoly_t_r, ppoly_s_r, stable_r, pol, por, kol, kor, heff)
2497  ndiff_unit_tests_discontinuous = ndiff_unit_tests_discontinuous .or. test_nsp(v, 12, kol, kor, pol, por, heff, &
2498  (/ 1, 1, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3 /), & ! KoL
2499  (/ 1, 1, 1, 1, 1, 1, 2, 2, 3, 3, 3, 3 /), & ! KoR
2500  (/ 0.00, 1.00, 0.00, 1.00, 1.00, 1.00, 1.00, 1.00, 0.00, 0.00, 1.00, 1.00 /), & ! PoL
2501  (/ 0.00, 0.00, 0.00, 0.00, 0.00, 1.00, 0.00, 1.00, 0.00, 0.50, 0.75, 1.00 /), & ! PoR
2502  (/ 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 4.00, 0.00 /), & ! hEff
2503  'Left thick mixed layer, Right unstable mixed')
2504 
2505  til(1,:) = (/ 8.00, 12.00 /); til(2,:) = (/ 12.00, 10.00 /); til(3,:) = (/ 8.00, 4.00 /);
2506  tir(1,:) = (/ 10.00, 14.00 /); tir(2,:) = (/ 14.00, 12.00 /); tir(3,:) = (/ 10.00, 6.00 /);
2507  call mark_unstable_cells( cs, nk, til, sil, pres_l, stable_l )
2508  call mark_unstable_cells( cs, nk, tir, sir, pres_r, stable_r )
2509  call find_neutral_surface_positions_discontinuous(cs, nk, pres_l, hl, til, sil, ppoly_t_l, ppoly_s_l, stable_l, &
2510  pres_r, hr, tir, sir, ppoly_t_r, ppoly_s_r, stable_r, pol, por, kol, kor, heff)
2511  ndiff_unit_tests_discontinuous = ndiff_unit_tests_discontinuous .or. test_nsp(v, 12, kol, kor, pol, por, heff, &
2512  (/ 1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3 /), & ! KoL
2513  (/ 1, 1, 1, 1, 2, 2, 2, 3, 3, 3, 3, 3 /), & ! KoR
2514  (/ 0.00, 1.00, 1.00, 1.00, 0.00, 0.00, 0.00, 1.00, 0.00, 0.00, 0.50, 1.00 /), & ! PoL
2515  (/ 0.00, 0.00, 0.00, 1.00, 0.00, 1.00, 1.00, 0.00, 0.00, 0.50, 1.00, 1.00 /), & ! PoR
2516  (/ 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 5.00, 0.00 /), & ! hEff
2517  'Unstable mixed layers, left cooler')
2518 
2519  call eos_manual_init(cs%EOS, form_of_eos = eos_linear, drho_dt = -1., drho_ds = 2.)
2520  ! Tests for linearized version of searching the layer for neutral surface position
2521  ! EOS linear in T, uniform alpha
2522  cs%max_iter = 10
2523  ! Unit tests require explicit initialization of tolerance
2524  cs%Drho_tol = 0.
2525  cs%x_tol = 0.
2526  ndiff_unit_tests_discontinuous = ndiff_unit_tests_discontinuous .or. (test_rnp(0.5, &
2527  find_neutral_pos_linear(cs, 0., 10., 35., 0., -0.2, 0., &
2528  0., -0.2, 0., 10., -0.2, 0., &
2529  (/12.,-4./), (/34.,0./)), "Temp Uniform Linearized Alpha/Beta"))
2530  ! EOS linear in S, uniform beta
2531  ndiff_unit_tests_discontinuous = ndiff_unit_tests_discontinuous .or. (test_rnp(0.5, &
2532  find_neutral_pos_linear(cs, 0., 10., 35., 0., 0., 0.8, &
2533  0., 0., 0.8, 10., 0., 0.8, &
2534  (/12.,0./), (/34.,2./)), "Salt Uniform Linearized Alpha/Beta"))
2535  ! EOS linear in T/S, uniform alpha/beta
2536  ndiff_unit_tests_discontinuous = ndiff_unit_tests_discontinuous .or. (test_rnp(0.5, &
2537  find_neutral_pos_linear(cs, 0., 10., 35., 0., -0.5, 0.5, &
2538  0., -0.5, 0.5, 10., -0.5, 0.5, &
2539  (/12.,-4./), (/34.,2./)), "Temp/salt Uniform Linearized Alpha/Beta"))
2540  ! EOS linear in T, insensitive to So
2541  ndiff_unit_tests_discontinuous = ndiff_unit_tests_discontinuous .or. (test_rnp(0.5, &
2542  find_neutral_pos_linear(cs, 0., 10., 35., 0., -0.2, 0., &
2543  0., -0.4, 0., 10., -0.6, 0., &
2544  (/12.,-4./), (/34.,0./)), "Temp stratified Linearized Alpha/Beta"))
2545 ! ! EOS linear in S, insensitive to T
2546  ndiff_unit_tests_discontinuous = ndiff_unit_tests_discontinuous .or. (test_rnp(0.5, &
2547  find_neutral_pos_linear(cs, 0., 10., 35., 0., 0., 0.8, &
2548  0., 0., 1.0, 10., 0., 0.5, &
2549  (/12.,0./), (/34.,2./)), "Salt stratified Linearized Alpha/Beta"))
2550  if (.not. ndiff_unit_tests_discontinuous) write(*,*) 'Pass'
2551 

References mom_eos::eos_linear, find_neutral_pos_linear(), find_neutral_surface_positions_discontinuous(), mark_unstable_cells(), test_nsp(), and test_rnp().

Referenced by neutral_diffusion_unit_tests().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ neutral_diffusion()

subroutine, public mom_neutral_diffusion::neutral_diffusion ( type(ocean_grid_type), intent(in)  G,
type(verticalgrid_type), intent(in)  GV,
real, dimension( g %isd: g %ied, g %jsd: g %jed, g %ke), intent(in)  h,
real, dimension( g %isdb: g %iedb, g %jsd: g %jed), intent(in)  Coef_x,
real, dimension( g %isd: g %ied, g %jsdb: g %jedb), intent(in)  Coef_y,
real, intent(in)  dt,
type(tracer_registry_type), pointer  Reg,
type(unit_scale_type), intent(in)  US,
type(neutral_diffusion_cs), pointer  CS 
)

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

Parameters
[in]gOcean grid structure
[in]gvocean vertical grid structure
[in]hLayer thickness [H ~> m or kg m-2]
[in]coef_xdt * Kh * dy / dx at u-points [L2 ~> m2]
[in]coef_ydt * Kh * dx / dy at v-points [L2 ~> m2]
[in]dtTracer time step * I_numitts [T ~> s] (I_numitts in tracer_hordiff)
regTracer registry
[in]usA dimensional unit scaling type
csNeutral diffusion control structure

Definition at line 493 of file MOM_neutral_diffusion.F90.

493  type(ocean_grid_type), intent(in) :: G !< Ocean grid structure
494  type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure
495  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h !< Layer thickness [H ~> m or kg m-2]
496  real, dimension(SZIB_(G),SZJ_(G)), intent(in) :: Coef_x !< dt * Kh * dy / dx at u-points [L2 ~> m2]
497  real, dimension(SZI_(G),SZJB_(G)), intent(in) :: Coef_y !< dt * Kh * dx / dy at v-points [L2 ~> m2]
498  real, intent(in) :: dt !< Tracer time step * I_numitts [T ~> s]
499  !! (I_numitts in tracer_hordiff)
500  type(tracer_registry_type), pointer :: Reg !< Tracer registry
501  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
502  type(neutral_diffusion_CS), pointer :: CS !< Neutral diffusion control structure
503 
504  ! Local variables
505  real, dimension(SZIB_(G),SZJ_(G),CS%nsurf-1) :: uFlx ! Zonal flux of tracer [H conc ~> m conc or conc kg m-2]
506  real, dimension(SZI_(G),SZJB_(G),CS%nsurf-1) :: vFlx ! Meridional flux of tracer
507  ! [H conc ~> m conc or conc kg m-2]
508  real, dimension(SZI_(G),SZJ_(G),G%ke) :: tendency ! tendency array for diagn
509  real, dimension(SZI_(G),SZJ_(G)) :: tendency_2d ! depth integrated content tendency for diagn
510  real, dimension(SZIB_(G),SZJ_(G)) :: trans_x_2d ! depth integrated diffusive tracer x-transport diagn
511  real, dimension(SZI_(G),SZJB_(G)) :: trans_y_2d ! depth integrated diffusive tracer y-transport diagn
512  real, dimension(G%ke) :: dTracer ! change in tracer concentration due to ndiffusion
513 
514  type(tracer_type), pointer :: Tracer => null() ! Pointer to the current tracer
515 
516  integer :: i, j, k, m, ks, nk
517  real :: Idt ! The inverse of the time step [T-1 ~> s-1]
518  real :: h_neglect, h_neglect_edge
519 
520  !### Try replacing both of these with GV%H_subroundoff
521  h_neglect_edge = gv%m_to_H*1.0e-10
522  h_neglect = gv%m_to_H*1.0e-30
523 
524  nk = gv%ke
525 
526  do m = 1,reg%ntr ! Loop over tracer registry
527 
528  tracer => reg%Tr(m)
529 
530  ! for diagnostics
531  if (tracer%id_dfxy_conc > 0 .or. tracer%id_dfxy_cont > 0 .or. tracer%id_dfxy_cont_2d > 0 .or. &
532  tracer%id_dfx_2d > 0 .or. tracer%id_dfy_2d > 0) then
533  idt = 1.0 / dt
534  tendency(:,:,:) = 0.0
535  endif
536 
537  uflx(:,:,:) = 0.
538  vflx(:,:,:) = 0.
539 
540  ! x-flux
541  do j = g%jsc,g%jec ; do i = g%isc-1,g%iec
542  if (g%mask2dCu(i,j)>0.) then
543  call neutral_surface_flux(nk, cs%nsurf, cs%deg, h(i,j,:), h(i+1,j,:), &
544  tracer%t(i,j,:), tracer%t(i+1,j,:), &
545  cs%uPoL(i,j,:), cs%uPoR(i,j,:), &
546  cs%uKoL(i,j,:), cs%uKoR(i,j,:), &
547  cs%uhEff(i,j,:), uflx(i,j,:), &
548  cs%continuous_reconstruction, h_neglect, cs%remap_CS, h_neglect_edge)
549  endif
550  enddo ; enddo
551 
552  ! y-flux
553  do j = g%jsc-1,g%jec ; do i = g%isc,g%iec
554  if (g%mask2dCv(i,j)>0.) then
555  call neutral_surface_flux(nk, cs%nsurf, cs%deg, h(i,j,:), h(i,j+1,:), &
556  tracer%t(i,j,:), tracer%t(i,j+1,:), &
557  cs%vPoL(i,j,:), cs%vPoR(i,j,:), &
558  cs%vKoL(i,j,:), cs%vKoR(i,j,:), &
559  cs%vhEff(i,j,:), vflx(i,j,:), &
560  cs%continuous_reconstruction, h_neglect, cs%remap_CS, h_neglect_edge)
561  endif
562  enddo ; enddo
563 
564  ! Update the tracer concentration from divergence of neutral diffusive flux components
565  do j = g%jsc,g%jec ; do i = g%isc,g%iec
566  if (g%mask2dT(i,j)>0.) then
567 
568  dtracer(:) = 0.
569  do ks = 1,cs%nsurf-1
570  k = cs%uKoL(i,j,ks)
571  dtracer(k) = dtracer(k) + coef_x(i,j) * uflx(i,j,ks)
572  k = cs%uKoR(i-1,j,ks)
573  dtracer(k) = dtracer(k) - coef_x(i-1,j) * uflx(i-1,j,ks)
574  k = cs%vKoL(i,j,ks)
575  dtracer(k) = dtracer(k) + coef_y(i,j) * vflx(i,j,ks)
576  k = cs%vKoR(i,j-1,ks)
577  dtracer(k) = dtracer(k) - coef_y(i,j-1) * vflx(i,j-1,ks)
578  enddo
579  do k = 1, gv%ke
580  tracer%t(i,j,k) = tracer%t(i,j,k) + dtracer(k) * &
581  ( g%IareaT(i,j) / ( h(i,j,k) + gv%H_subroundoff ) )
582  enddo
583 
584  if (tracer%id_dfxy_conc > 0 .or. tracer%id_dfxy_cont > 0 .or. tracer%id_dfxy_cont_2d > 0 ) then
585  do k = 1, gv%ke
586  tendency(i,j,k) = dtracer(k) * g%IareaT(i,j) * idt
587  enddo
588  endif
589 
590  endif
591  enddo ; enddo
592 
593  ! Diagnose vertically summed zonal flux, giving zonal tracer transport from ndiff.
594  ! Note sign corresponds to downgradient flux convention.
595  if (tracer%id_dfx_2d > 0) then
596  do j = g%jsc,g%jec ; do i = g%isc-1,g%iec
597  trans_x_2d(i,j) = 0.
598  if (g%mask2dCu(i,j)>0.) then
599  do ks = 1,cs%nsurf-1
600  trans_x_2d(i,j) = trans_x_2d(i,j) - coef_x(i,j) * uflx(i,j,ks)
601  enddo
602  trans_x_2d(i,j) = trans_x_2d(i,j) * idt
603  endif
604  enddo ; enddo
605  call post_data(tracer%id_dfx_2d, trans_x_2d(:,:), cs%diag)
606  endif
607 
608  ! Diagnose vertically summed merid flux, giving meridional tracer transport from ndiff.
609  ! Note sign corresponds to downgradient flux convention.
610  if (tracer%id_dfy_2d > 0) then
611  do j = g%jsc-1,g%jec ; do i = g%isc,g%iec
612  trans_y_2d(i,j) = 0.
613  if (g%mask2dCv(i,j)>0.) then
614  do ks = 1,cs%nsurf-1
615  trans_y_2d(i,j) = trans_y_2d(i,j) - coef_y(i,j) * vflx(i,j,ks)
616  enddo
617  trans_y_2d(i,j) = trans_y_2d(i,j) * idt
618  endif
619  enddo ; enddo
620  call post_data(tracer%id_dfy_2d, trans_y_2d(:,:), cs%diag)
621  endif
622 
623  ! post tendency of tracer content
624  if (tracer%id_dfxy_cont > 0) then
625  call post_data(tracer%id_dfxy_cont, tendency(:,:,:), cs%diag)
626  endif
627 
628  ! post depth summed tendency for tracer content
629  if (tracer%id_dfxy_cont_2d > 0) then
630  tendency_2d(:,:) = 0.
631  do j = g%jsc,g%jec ; do i = g%isc,g%iec
632  do k = 1, gv%ke
633  tendency_2d(i,j) = tendency_2d(i,j) + tendency(i,j,k)
634  enddo
635  enddo ; enddo
636  call post_data(tracer%id_dfxy_cont_2d, tendency_2d(:,:), cs%diag)
637  endif
638 
639  ! post tendency of tracer concentration; this step must be
640  ! done after posting tracer content tendency, since we alter
641  ! the tendency array.
642  if (tracer%id_dfxy_conc > 0) then
643  do k = 1, gv%ke ; do j = g%jsc,g%jec ; do i = g%isc,g%iec
644  tendency(i,j,k) = tendency(i,j,k) / ( h(i,j,k) + gv%H_subroundoff )
645  enddo ; enddo ; enddo
646  call post_data(tracer%id_dfxy_conc, tendency, cs%diag)
647  endif
648  enddo ! Loop over tracer registry
649 

References neutral_surface_flux().

Referenced by mom_tracer_hor_diff::tracer_hordiff().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ neutral_diffusion_calc_coeffs()

subroutine, public mom_neutral_diffusion::neutral_diffusion_calc_coeffs ( type(ocean_grid_type), intent(in)  G,
type(verticalgrid_type), intent(in)  GV,
type(unit_scale_type), intent(in)  US,
real, dimension( g %isd: g %ied, g %jsd: g %jed, g %ke), intent(in)  h,
real, dimension( g %isd: g %ied, g %jsd: g %jed, g %ke), intent(in)  T,
real, dimension( g %isd: g %ied, g %jsd: g %jed, g %ke), intent(in)  S,
type(neutral_diffusion_cs), pointer  CS 
)

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

Parameters
[in]gOcean grid structure
[in]gvocean vertical grid structure
[in]usA dimensional unit scaling type
[in]hLayer thickness [H ~> m or kg m-2]
[in]tPotential temperature [degC]
[in]sSalinity [ppt]
csNeutral diffusion control structure

Definition at line 264 of file MOM_neutral_diffusion.F90.

264  type(ocean_grid_type), intent(in) :: G !< Ocean grid structure
265  type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure
266  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
267  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h !< Layer thickness [H ~> m or kg m-2]
268  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: T !< Potential temperature [degC]
269  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: S !< Salinity [ppt]
270  type(neutral_diffusion_CS), pointer :: CS !< Neutral diffusion control structure
271 
272  ! Local variables
273  integer :: i, j, k
274  ! Variables used for reconstructions
275  real, dimension(SZK_(G),2) :: ppoly_r_S ! Reconstruction slopes
276  real, dimension(SZI_(G), SZJ_(G)) :: hEff_sum
277  real, dimension(SZI_(G),SZJ_(G)) :: hbl !< bnd. layer depth [m]
278  integer :: iMethod
279  real, dimension(SZI_(G)) :: ref_pres ! Reference pressure used to calculate alpha/beta
280  real, dimension(SZI_(G)) :: rho_tmp ! Routine to calculate drho_dp, returns density which is not used
281  real :: h_neglect, h_neglect_edge
282  integer, dimension(SZI_(G), SZJ_(G)) :: k_top !< Index of the first layer within the boundary
283  real, dimension(SZI_(G), SZJ_(G)) :: zeta_top !< Distance from the top of a layer to the intersection of the
284  !! top extent of the boundary layer (0 at top, 1 at bottom) [nondim]
285  integer, dimension(SZI_(G), SZJ_(G)) :: k_bot !< Index of the last layer within the boundary
286  real, dimension(SZI_(G), SZJ_(G)) :: zeta_bot !< Distance of the lower layer to the boundary layer depth
287  real :: pa_to_H
288 
289  pa_to_h = 1. / gv%H_to_pa
290 
291  k_top(:,:) = 1 ; k_bot(:,:) = 1
292  zeta_top(:,:) = 0. ; zeta_bot(:,:) = 1.
293 
294  ! check if hbl needs to be extracted
295  if (cs%interior_only) then
296  hbl(:,:) = 0.
297  if (ASSOCIATED(cs%KPP_CSp)) call kpp_get_bld(cs%KPP_CSp, hbl, g)
298  if (ASSOCIATED(cs%energetic_PBL_CSp)) call energetic_pbl_get_mld(cs%energetic_PBL_CSp, hbl, g, us)
299  call pass_var(hbl,g%Domain)
300  ! get k-indices and zeta
301  do j=g%jsc-1, g%jec+1 ; do i=g%isc-1,g%iec+1
302  call boundary_k_range(surface, g%ke, h(i,j,:), hbl(i,j), k_top(i,j), zeta_top(i,j), k_bot(i,j), zeta_bot(i,j))
303  enddo; enddo
304  ! TODO: add similar code for BOTTOM boundary layer
305  endif
306 
307  !### Try replacing both of these with GV%H_subroundoff
308  if (gv%Boussinesq) then
309  h_neglect = gv%m_to_H*1.0e-30 ; h_neglect_edge = gv%m_to_H*1.0e-10
310  else
311  h_neglect = gv%kg_m2_to_H*1.0e-30 ; h_neglect_edge = gv%kg_m2_to_H*1.0e-10
312  endif
313 
314  ! If doing along isopycnal diffusion (as opposed to neutral diffusion, set the reference pressure)
315  if (cs%ref_pres>=0.) then
316  ref_pres(:) = cs%ref_pres
317  endif
318 
319  if (cs%continuous_reconstruction) then
320  cs%dRdT(:,:,:) = 0.
321  cs%dRdS(:,:,:) = 0.
322  else
323  cs%T_i(:,:,:,:) = 0.
324  cs%S_i(:,:,:,:) = 0.
325  cs%dRdT_i(:,:,:,:) = 0.
326  cs%dRdS_i(:,:,:,:) = 0.
327  cs%ns(:,:) = 0.
328  cs%stable_cell(:,:,:) = .true.
329  endif
330 
331  ! Calculate pressure at interfaces and layer averaged alpha/beta
332  cs%Pint(:,:,1) = 0.
333  do k=1,g%ke ; do j=g%jsc-1, g%jec+1 ; do i=g%isc-1,g%iec+1
334  cs%Pint(i,j,k+1) = cs%Pint(i,j,k) + h(i,j,k)*gv%H_to_Pa
335  enddo ; enddo ; enddo
336 
337  ! Pressures at the interfaces, this is redundant as P_i(k,1) = P_i(k-1,2) however retain tis
338  ! for now ensure consitency of indexing for diiscontinuous reconstructions
339  if (.not. cs%continuous_reconstruction) then
340  do j=g%jsc-1, g%jec+1 ; do i=g%isc-1,g%iec+1
341  cs%P_i(i,j,1,1) = 0.
342  cs%P_i(i,j,1,2) = h(i,j,1)*gv%H_to_Pa
343  enddo ; enddo
344  do k=2,g%ke ; do j=g%jsc-1, g%jec+1 ; do i=g%isc-1,g%iec+1
345  cs%P_i(i,j,k,1) = cs%P_i(i,j,k-1,2)
346  cs%P_i(i,j,k,2) = cs%P_i(i,j,k-1,2) + h(i,j,k)*gv%H_to_Pa
347  enddo ; enddo ; enddo
348  endif
349 
350  do j = g%jsc-1, g%jec+1
351  ! Interpolate state to interface
352  do i = g%isc-1, g%iec+1
353  if (cs%continuous_reconstruction) then
354  call interface_scalar(g%ke, h(i,j,:), t(i,j,:), cs%Tint(i,j,:), 2, h_neglect)
355  call interface_scalar(g%ke, h(i,j,:), s(i,j,:), cs%Sint(i,j,:), 2, h_neglect)
356  else
357  call build_reconstructions_1d( cs%remap_CS, g%ke, h(i,j,:), t(i,j,:), cs%ppoly_coeffs_T(i,j,:,:), &
358  cs%T_i(i,j,:,:), ppoly_r_s, imethod, h_neglect, h_neglect_edge )
359  call build_reconstructions_1d( cs%remap_CS, g%ke, h(i,j,:), s(i,j,:), cs%ppoly_coeffs_S(i,j,:,:), &
360  cs%S_i(i,j,:,:), ppoly_r_s, imethod, h_neglect, h_neglect_edge )
361  ! In the current ALE formulation, interface values are not exactly at the 0. or 1. of the
362  ! polynomial reconstructions
363  do k=1,g%ke
364  cs%T_i(i,j,k,1) = evaluation_polynomial( cs%ppoly_coeffs_T(i,j,k,:), cs%deg+1, 0. )
365  cs%T_i(i,j,k,2) = evaluation_polynomial( cs%ppoly_coeffs_T(i,j,k,:), cs%deg+1, 1. )
366  cs%S_i(i,j,k,1) = evaluation_polynomial( cs%ppoly_coeffs_S(i,j,k,:), cs%deg+1, 0. )
367  cs%S_i(i,j,k,2) = evaluation_polynomial( cs%ppoly_coeffs_S(i,j,k,:), cs%deg+1, 1. )
368  enddo
369  endif
370  enddo
371 
372  ! Continuous reconstruction
373  if (cs%continuous_reconstruction) then
374  do k = 1, g%ke+1
375  if (cs%ref_pres<0) ref_pres(:) = cs%Pint(:,j,k)
376  call calculate_density_derivs(cs%Tint(:,j,k), cs%Sint(:,j,k), ref_pres, &
377  cs%dRdT(:,j,k), cs%dRdS(:,j,k), g%isc-1, g%iec-g%isc+3, cs%EOS)
378  enddo
379  else ! Discontinuous reconstruction
380  do k = 1, g%ke
381  if (cs%ref_pres<0) ref_pres(:) = cs%Pint(:,j,k)
382  ! Calculate derivatives for the top interface
383  call calculate_density_derivs(cs%T_i(:,j,k,1), cs%S_i(:,j,k,1), ref_pres, &
384  cs%dRdT_i(:,j,k,1), cs%dRdS_i(:,j,k,1), g%isc-1, g%iec-g%isc+3, cs%EOS)
385  if (cs%ref_pres<0) then
386  ref_pres(:) = cs%Pint(:,j,k+1)
387  endif
388  ! Calcualte derivatives at the bottom interface
389  call calculate_density_derivs(cs%T_i(:,j,k,2), cs%S_i(:,j,k,2), ref_pres, &
390  cs%dRdT_i(:,j,k,2), cs%dRdS_i(:,j,k,2), g%isc-1, g%iec-g%isc+3, cs%EOS)
391  enddo
392  endif
393  enddo
394 
395  if (.not. cs%continuous_reconstruction) then
396  do j = g%jsc-1, g%jec+1 ; do i = g%isc-1, g%iec+1
397  call mark_unstable_cells( cs, g%ke, cs%T_i(i,j,:,:), cs%S_i(i,j,:,:), cs%P_i(i,j,:,:), cs%stable_cell(i,j,:) )
398  if (cs%interior_only) then
399  if (.not. cs%stable_cell(i,j,k_bot(i,j))) zeta_bot(i,j) = -1.
400  ! set values in the surface and bottom boundary layer to false.
401  do k = 1, k_bot(i,j)-1
402  cs%stable_cell(i,j,k) = .false.
403  enddo
404  endif
405  enddo ; enddo
406  endif
407 
408  cs%uhEff(:,:,:) = 0.
409  cs%vhEff(:,:,:) = 0.
410  cs%uPoL(:,:,:) = 0.
411  cs%vPoL(:,:,:) = 0.
412  cs%uPoR(:,:,:) = 0.
413  cs%vPoR(:,:,:) = 0.
414  cs%uKoL(:,:,:) = 1
415  cs%vKoL(:,:,:) = 1
416  cs%uKoR(:,:,:) = 1
417  cs%vKoR(:,:,:) = 1
418 
419  ! Neutral surface factors at U points
420  do j = g%jsc, g%jec ; do i = g%isc-1, g%iec
421  if (g%mask2dCu(i,j) > 0.) then
422  if (cs%continuous_reconstruction) then
423  call find_neutral_surface_positions_continuous(g%ke, &
424  cs%Pint(i,j,:), cs%Tint(i,j,:), cs%Sint(i,j,:), cs%dRdT(i,j,:), cs%dRdS(i,j,:), &
425  cs%Pint(i+1,j,:), cs%Tint(i+1,j,:), cs%Sint(i+1,j,:), cs%dRdT(i+1,j,:), cs%dRdS(i+1,j,:), &
426  cs%uPoL(i,j,:), cs%uPoR(i,j,:), cs%uKoL(i,j,:), cs%uKoR(i,j,:), cs%uhEff(i,j,:), &
427  k_bot(i,j), k_bot(i+1,j), 1.-zeta_bot(i,j), 1.-zeta_bot(i+1,j))
428  else
429  call find_neutral_surface_positions_discontinuous(cs, g%ke, &
430  cs%P_i(i,j,:,:), h(i,j,:), cs%T_i(i,j,:,:), cs%S_i(i,j,:,:), cs%ppoly_coeffs_T(i,j,:,:), &
431  cs%ppoly_coeffs_S(i,j,:,:),cs%stable_cell(i,j,:), &
432  cs%P_i(i+1,j,:,:), h(i+1,j,:), cs%T_i(i+1,j,:,:), cs%S_i(i+1,j,:,:), cs%ppoly_coeffs_T(i+1,j,:,:), &
433  cs%ppoly_coeffs_S(i+1,j,:,:), cs%stable_cell(i+1,j,:), &
434  cs%uPoL(i,j,:), cs%uPoR(i,j,:), cs%uKoL(i,j,:), cs%uKoR(i,j,:), cs%uhEff(i,j,:), &
435  hard_fail_heff = cs%hard_fail_heff)
436  endif
437  endif
438  enddo ; enddo
439 
440  ! Neutral surface factors at V points
441  do j = g%jsc-1, g%jec ; do i = g%isc, g%iec
442  if (g%mask2dCv(i,j) > 0.) then
443  if (cs%continuous_reconstruction) then
444  call find_neutral_surface_positions_continuous(g%ke, &
445  cs%Pint(i,j,:), cs%Tint(i,j,:), cs%Sint(i,j,:), cs%dRdT(i,j,:), cs%dRdS(i,j,:), &
446  cs%Pint(i,j+1,:), cs%Tint(i,j+1,:), cs%Sint(i,j+1,:), cs%dRdT(i,j+1,:), cs%dRdS(i,j+1,:), &
447  cs%vPoL(i,j,:), cs%vPoR(i,j,:), cs%vKoL(i,j,:), cs%vKoR(i,j,:), cs%vhEff(i,j,:), &
448  k_bot(i,j), k_bot(i,j+1), 1.-zeta_bot(i,j), 1.-zeta_bot(i,j+1))
449  else
450  call find_neutral_surface_positions_discontinuous(cs, g%ke, &
451  cs%P_i(i,j,:,:), h(i,j,:), cs%T_i(i,j,:,:), cs%S_i(i,j,:,:), cs%ppoly_coeffs_T(i,j,:,:), &
452  cs%ppoly_coeffs_S(i,j,:,:),cs%stable_cell(i,j,:), &
453  cs%P_i(i,j+1,:,:), h(i,j+1,:), cs%T_i(i,j+1,:,:), cs%S_i(i,j+1,:,:), cs%ppoly_coeffs_T(i,j+1,:,:), &
454  cs%ppoly_coeffs_S(i,j+1,:,:), cs%stable_cell(i,j+1,:), &
455  cs%vPoL(i,j,:), cs%vPoR(i,j,:), cs%vKoL(i,j,:), cs%vKoR(i,j,:), cs%vhEff(i,j,:), &
456  hard_fail_heff = cs%hard_fail_heff)
457  endif
458  endif
459  enddo ; enddo
460 
461  ! Continuous reconstructions calculate hEff as the difference between the pressures of the
462  ! neutral surfaces which need to be reconverted to thickness units. The discontinuous version
463  ! calculates hEff from the fraction of the nondimensional fraction of the layer spanned by
464  ! adjacent neutral surfaces.
465  if (cs%continuous_reconstruction) then
466  do k = 1, cs%nsurf-1 ; do j = g%jsc, g%jec ; do i = g%isc-1, g%iec
467  if (g%mask2dCu(i,j) > 0.) cs%uhEff(i,j,k) = cs%uhEff(i,j,k) * pa_to_h
468  enddo ; enddo ; enddo
469  do k = 1, cs%nsurf-1 ; do j = g%jsc-1, g%jec ; do i = g%isc, g%iec
470  if (g%mask2dCv(i,j) > 0.) cs%vhEff(i,j,k) = cs%vhEff(i,j,k) * pa_to_h
471  enddo ; enddo ; enddo
472  endif
473 
474  if (cs%id_uhEff_2d>0) then
475  heff_sum(:,:) = 0.
476  do k = 1,cs%nsurf-1 ; do j=g%jsc,g%jec ; do i=g%isc-1,g%iec
477  heff_sum(i,j) = heff_sum(i,j) + cs%uhEff(i,j,k)
478  enddo ; enddo ; enddo
479  call post_data(cs%id_uhEff_2d, heff_sum, cs%diag)
480  endif
481  if (cs%id_vhEff_2d>0) then
482  heff_sum(:,:) = 0.
483  do k = 1,cs%nsurf-1 ; do j=g%jsc-1,g%jec ; do i=g%isc,g%iec
484  heff_sum(i,j) = heff_sum(i,j) + cs%vhEff(i,j,k)
485  enddo ; enddo ; enddo
486  call post_data(cs%id_vhEff_2d, heff_sum, cs%diag)
487  endif
488 

References 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(), and mom_lateral_boundary_diffusion::surface.

Referenced by mom_tracer_hor_diff::tracer_hordiff().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ neutral_diffusion_end()

subroutine, public mom_neutral_diffusion::neutral_diffusion_end ( type(neutral_diffusion_cs), pointer  CS)

Deallocates neutral_diffusion control structure.

Parameters
csNeutral diffusion control structure

Definition at line 2807 of file MOM_neutral_diffusion.F90.

2807  type(neutral_diffusion_CS), pointer :: CS !< Neutral diffusion control structure
2808 
2809  if (associated(cs)) deallocate(cs)
2810 

◆ neutral_diffusion_init()

logical function, public mom_neutral_diffusion::neutral_diffusion_init ( type(time_type), intent(in), target  Time,
type(ocean_grid_type), intent(in)  G,
type(param_file_type), intent(in)  param_file,
type(diag_ctrl), intent(inout), target  diag,
type(eos_type), intent(in), target  EOS,
type(diabatic_cs), pointer  diabatic_CSp,
type(neutral_diffusion_cs), pointer  CS 
)

Read parameters and allocate control structure for neutral_diffusion module.

Parameters
[in]timeTime structure
[in]gGrid structure
[in,out]diagDiagnostics control structure
[in]param_fileParameter file structure
[in]eosEquation of state
diabatic_cspKPP control structure needed to get BLD
csNeutral diffusion control structure

Definition at line 109 of file MOM_neutral_diffusion.F90.

109  type(time_type), target, intent(in) :: Time !< Time structure
110  type(ocean_grid_type), intent(in) :: G !< Grid structure
111  type(diag_ctrl), target, intent(inout) :: diag !< Diagnostics control structure
112  type(param_file_type), intent(in) :: param_file !< Parameter file structure
113  type(EOS_type), target, intent(in) :: EOS !< Equation of state
114  type(diabatic_CS), pointer :: diabatic_CSp!< KPP control structure needed to get BLD
115  type(neutral_diffusion_CS), pointer :: CS !< Neutral diffusion control structure
116 
117  ! Local variables
118  character(len=256) :: mesg ! Message for error messages.
119  character(len=80) :: string ! Temporary strings
120  logical :: boundary_extrap
121 
122  if (associated(cs)) then
123  call mom_error(fatal, "neutral_diffusion_init called with associated control structure.")
124  return
125  endif
126 
127 
128  ! Log this module and master switch for turning it on/off
129  call log_version(param_file, mdl, version, &
130  "This module implements neutral diffusion of tracers")
131  call get_param(param_file, mdl, "USE_NEUTRAL_DIFFUSION", neutral_diffusion_init, &
132  "If true, enables the neutral diffusion module.", &
133  default=.false.)
134 
135  if (.not.neutral_diffusion_init) then
136  return
137  endif
138 
139  allocate(cs)
140  cs%diag => diag
141  cs%EOS => eos
142  ! call openParameterBlock(param_file,'NEUTRAL_DIFF')
143 
144  ! Read all relevant parameters and write them to the model log.
145  call get_param(param_file, mdl, "NDIFF_CONTINUOUS", cs%continuous_reconstruction, &
146  "If true, uses a continuous reconstruction of T and S when "//&
147  "finding neutral surfaces along which diffusion will happen. "//&
148  "If false, a PPM discontinuous reconstruction of T and S "//&
149  "is done which results in a higher order routine but exacts "//&
150  "a higher computational cost.", default=.true.)
151  call get_param(param_file, mdl, "NDIFF_REF_PRES", cs%ref_pres, &
152  "The reference pressure (Pa) used for the derivatives of "//&
153  "the equation of state. If negative (default), local "//&
154  "pressure is used.", &
155  default = -1.)
156  call get_param(param_file, mdl, "NDIFF_INTERIOR_ONLY", cs%interior_only, &
157  "If true, only applies neutral diffusion in the ocean interior."//&
158  "That is, the algorithm will exclude the surface and bottom"//&
159  "boundary layers.",default = .false.)
160 
161  ! Initialize and configure remapping
162  if (cs%continuous_reconstruction .eqv. .false.) then
163  call get_param(param_file, mdl, "NDIFF_BOUNDARY_EXTRAP", boundary_extrap, &
164  "Extrapolate at the top and bottommost cells, otherwise \n"// &
165  "assume boundaries are piecewise constant", &
166  default=.false.)
167  call get_param(param_file, mdl, "NDIFF_REMAPPING_SCHEME", string, &
168  "This sets the reconstruction scheme used "//&
169  "for vertical remapping for all variables. "//&
170  "It can be one of the following schemes: "//&
171  trim(remappingschemesdoc), default=remappingdefaultscheme)
172  call initialize_remapping( cs%remap_CS, string, boundary_extrapolation = boundary_extrap )
173  call extract_member_remapping_cs(cs%remap_CS, degree=cs%deg)
174  call get_param(param_file, mdl, "NEUTRAL_POS_METHOD", cs%neutral_pos_method, &
175  "Method used to find the neutral position \n"// &
176  "1. Delta_rho varies linearly, find 0 crossing \n"// &
177  "2. Alpha and beta vary linearly from top to bottom, \n"// &
178  " Newton's method for neutral position \n"// &
179  "3. Full nonlinear equation of state, use regula falsi \n"// &
180  " for neutral position", default=3)
181  if (cs%neutral_pos_method > 4 .or. cs%neutral_pos_method < 0) then
182  call mom_error(fatal,"Invalid option for NEUTRAL_POS_METHOD")
183  endif
184 
185  call get_param(param_file, mdl, "DELTA_RHO_FORM", cs%delta_rho_form, &
186  "Determine how the difference in density is calculated \n"// &
187  " full : Difference of in-situ densities \n"// &
188  " no_pressure: Calculated from dRdT, dRdS, but no \n"// &
189  " pressure dependence", &
190  default="mid_pressure")
191  if (cs%neutral_pos_method > 1) then
192  call get_param(param_file, mdl, "NDIFF_DRHO_TOL", cs%drho_tol, &
193  "Sets the convergence criterion for finding the neutral\n"// &
194  "position within a layer in kg m-3.", &
195  default=1.e-10)
196  call get_param(param_file, mdl, "NDIFF_X_TOL", cs%x_tol, &
197  "Sets the convergence criterion for a change in nondim\n"// &
198  "position within a layer.", &
199  default=0.)
200  call get_param(param_file, mdl, "NDIFF_MAX_ITER", cs%max_iter, &
201  "The maximum number of iterations to be done before \n"// &
202  "exiting the iterative loop to find the neutral surface", &
203  default=10)
204  endif
205  call get_param(param_file, mdl, "NDIFF_DEBUG", cs%debug, &
206  "Turns on verbose output for discontinuous neutral "//&
207  "diffusion routines.", &
208  default = .false.)
209  call get_param(param_file, mdl, "HARD_FAIL_HEFF", cs%hard_fail_heff, &
210  "Bring down the model if a problem with heff is detected",&
211  default = .true.)
212  endif
213 
214  if (cs%interior_only) then
215  call extract_diabatic_member(diabatic_csp, kpp_csp=cs%KPP_CSp)
216  call extract_diabatic_member(diabatic_csp, energetic_pbl_csp=cs%energetic_PBL_CSp)
217  if ( .not. ASSOCIATED(cs%energetic_PBL_CSp) .and. .not. ASSOCIATED(cs%KPP_CSp) ) then
218  call mom_error(fatal,"NDIFF_INTERIOR_ONLY is true, but no valid boundary layer scheme was found")
219  endif
220  endif
221 
222 ! call get_param(param_file, mdl, "KHTR", CS%KhTr, &
223 ! "The background along-isopycnal tracer diffusivity.", &
224 ! units="m2 s-1", default=0.0)
225 ! call closeParameterBlock(param_file)
226  if (cs%continuous_reconstruction) then
227  cs%nsurf = 2*g%ke+2 ! Continuous reconstruction means that every interface has two connections
228  allocate(cs%dRdT(szi_(g),szj_(g),szk_(g)+1)) ; cs%dRdT(:,:,:) = 0.
229  allocate(cs%dRdS(szi_(g),szj_(g),szk_(g)+1)) ; cs%dRdS(:,:,:) = 0.
230  else
231  cs%nsurf = 4*g%ke ! Discontinuous means that every interface has four connections
232  allocate(cs%T_i(szi_(g),szj_(g),szk_(g),2)) ; cs%T_i(:,:,:,:) = 0.
233  allocate(cs%S_i(szi_(g),szj_(g),szk_(g),2)) ; cs%S_i(:,:,:,:) = 0.
234  allocate(cs%P_i(szi_(g),szj_(g),szk_(g),2)) ; cs%P_i(:,:,:,:) = 0.
235  allocate(cs%dRdT_i(szi_(g),szj_(g),szk_(g),2)) ; cs%dRdT_i(:,:,:,:) = 0.
236  allocate(cs%dRdS_i(szi_(g),szj_(g),szk_(g),2)) ; cs%dRdS_i(:,:,:,:) = 0.
237  allocate(cs%ppoly_coeffs_T(szi_(g),szj_(g),szk_(g),cs%deg+1)) ; cs%ppoly_coeffs_T(:,:,:,:) = 0.
238  allocate(cs%ppoly_coeffs_S(szi_(g),szj_(g),szk_(g),cs%deg+1)) ; cs%ppoly_coeffs_S(:,:,:,:) = 0.
239  allocate(cs%ns(szi_(g),szj_(g))) ; cs%ns(:,:) = 0.
240  endif
241  ! T-points
242  allocate(cs%Tint(szi_(g),szj_(g),szk_(g)+1)) ; cs%Tint(:,:,:) = 0.
243  allocate(cs%Sint(szi_(g),szj_(g),szk_(g)+1)) ; cs%Sint(:,:,:) = 0.
244  allocate(cs%Pint(szi_(g),szj_(g),szk_(g)+1)) ; cs%Pint(:,:,:) = 0.
245  allocate(cs%stable_cell(szi_(g),szj_(g),szk_(g))) ; cs%stable_cell(:,:,:) = .true.
246  ! U-points
247  allocate(cs%uPoL(g%isd:g%ied,g%jsd:g%jed, cs%nsurf)); cs%uPoL(g%isc-1:g%iec,g%jsc:g%jec,:) = 0.
248  allocate(cs%uPoR(g%isd:g%ied,g%jsd:g%jed, cs%nsurf)); cs%uPoR(g%isc-1:g%iec,g%jsc:g%jec,:) = 0.
249  allocate(cs%uKoL(g%isd:g%ied,g%jsd:g%jed, cs%nsurf)); cs%uKoL(g%isc-1:g%iec,g%jsc:g%jec,:) = 0
250  allocate(cs%uKoR(g%isd:g%ied,g%jsd:g%jed, cs%nsurf)); cs%uKoR(g%isc-1:g%iec,g%jsc:g%jec,:) = 0
251  allocate(cs%uHeff(g%isd:g%ied,g%jsd:g%jed,cs%nsurf-1)); cs%uHeff(g%isc-1:g%iec,g%jsc:g%jec,:) = 0
252  ! V-points
253  allocate(cs%vPoL(g%isd:g%ied,g%jsd:g%jed, cs%nsurf)); cs%vPoL(g%isc:g%iec,g%jsc-1:g%jec,:) = 0.
254  allocate(cs%vPoR(g%isd:g%ied,g%jsd:g%jed, cs%nsurf)); cs%vPoR(g%isc:g%iec,g%jsc-1:g%jec,:) = 0.
255  allocate(cs%vKoL(g%isd:g%ied,g%jsd:g%jed, cs%nsurf)); cs%vKoL(g%isc:g%iec,g%jsc-1:g%jec,:) = 0
256  allocate(cs%vKoR(g%isd:g%ied,g%jsd:g%jed, cs%nsurf)); cs%vKoR(g%isc:g%iec,g%jsc-1:g%jec,:) = 0
257  allocate(cs%vHeff(g%isd:g%ied,g%jsd:g%jed,cs%nsurf-1)); cs%vHeff(g%isc:g%iec,g%jsc-1:g%jec,:) = 0
258 

References mom_diabatic_driver::extract_diabatic_member(), mdl, mom_error_handler::mom_error(), mom_remapping::remappingdefaultscheme, and mom_remapping::remappingschemesdoc.

Here is the call graph for this function:

◆ neutral_diffusion_unit_tests()

logical function, public mom_neutral_diffusion::neutral_diffusion_unit_tests ( logical, intent(in)  verbose)

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

Parameters
[in]verboseIf true, write results to stdout

Definition at line 2024 of file MOM_neutral_diffusion.F90.

2024  logical, intent(in) :: verbose !< If true, write results to stdout
2025 
2026  neutral_diffusion_unit_tests = .false. .or. &
2027  ndiff_unit_tests_continuous(verbose) .or. ndiff_unit_tests_discontinuous(verbose)
2028 
2029 

References ndiff_unit_tests_continuous(), and ndiff_unit_tests_discontinuous().

Referenced by mom_unit_tests::unit_tests().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ neutral_surface_flux()

subroutine mom_neutral_diffusion::neutral_surface_flux ( integer, intent(in)  nk,
integer, intent(in)  nsurf,
integer, intent(in)  deg,
real, dimension(nk), intent(in)  hl,
real, dimension(nk), intent(in)  hr,
real, dimension(nk), intent(in)  Tl,
real, dimension(nk), intent(in)  Tr,
real, dimension(nsurf), intent(in)  PiL,
real, dimension(nsurf), intent(in)  PiR,
integer, dimension(nsurf), intent(in)  KoL,
integer, dimension(nsurf), intent(in)  KoR,
real, dimension(nsurf-1), intent(in)  hEff,
real, dimension(nsurf-1), intent(inout)  Flx,
logical, intent(in)  continuous,
real, intent(in)  h_neglect,
type(remapping_cs), intent(in), optional  remap_CS,
real, intent(in), optional  h_neglect_edge 
)
private

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

Parameters
[in]nkNumber of levels
[in]nsurfNumber of neutral surfaces
[in]degDegree of polynomial reconstructions
[in]hlLeft-column layer thickness [Pa]
[in]hrRight-column layer thickness [Pa]
[in]tlLeft-column layer tracer (conc, e.g. degC)
[in]trRight-column layer tracer (conc, e.g. degC)
[in]pilFractional position of neutral surface within layer KoL of left column
[in]pirFractional position of neutral surface within layer KoR of right column
[in]kolIndex of first left interface above neutral surface
[in]korIndex of first right interface above neutral surface
[in]heffEffective thickness between two neutral surfaces [Pa]
[in,out]flxFlux of tracer between pairs of neutral layers (conc H)
[in]continuousTrue if using continuous reconstruction
[in]h_neglectA negligibly small width for the purpose of cell reconstructions in the same units as h0.
[in]remap_csRemapping control structure used to create sublayers
[in]h_neglect_edgeA negligibly small width for the purpose of edge value calculations in the same units as h0.

Definition at line 1820 of file MOM_neutral_diffusion.F90.

1820  integer, intent(in) :: nk !< Number of levels
1821  integer, intent(in) :: nsurf !< Number of neutral surfaces
1822  integer, intent(in) :: deg !< Degree of polynomial reconstructions
1823  real, dimension(nk), intent(in) :: hl !< Left-column layer thickness [Pa]
1824  real, dimension(nk), intent(in) :: hr !< Right-column layer thickness [Pa]
1825  real, dimension(nk), intent(in) :: Tl !< Left-column layer tracer (conc, e.g. degC)
1826  real, dimension(nk), intent(in) :: Tr !< Right-column layer tracer (conc, e.g. degC)
1827  real, dimension(nsurf), intent(in) :: PiL !< Fractional position of neutral surface
1828  !! within layer KoL of left column
1829  real, dimension(nsurf), intent(in) :: PiR !< Fractional position of neutral surface
1830  !! within layer KoR of right column
1831  integer, dimension(nsurf), intent(in) :: KoL !< Index of first left interface above neutral surface
1832  integer, dimension(nsurf), intent(in) :: KoR !< Index of first right interface above neutral surface
1833  real, dimension(nsurf-1), intent(in) :: hEff !< Effective thickness between two neutral surfaces [Pa]
1834  real, dimension(nsurf-1), intent(inout) :: Flx !< Flux of tracer between pairs of neutral layers (conc H)
1835  logical, intent(in) :: continuous !< True if using continuous reconstruction
1836  real, intent(in) :: h_neglect !< A negligibly small width for the
1837  !! purpose of cell reconstructions
1838  !! in the same units as h0.
1839  type(remapping_CS), optional, intent(in) :: remap_CS !< Remapping control structure used
1840  !! to create sublayers
1841  real, optional, intent(in) :: h_neglect_edge !< A negligibly small width
1842  !! for the purpose of edge value calculations
1843  !! in the same units as h0.
1844  ! Local variables
1845  integer :: k_sublayer, klb, klt, krb, krt, k
1846  real :: T_right_top, T_right_bottom, T_right_layer, T_right_sub, T_right_top_int, T_right_bot_int
1847  real :: T_left_top, T_left_bottom, T_left_layer, T_left_sub, T_left_top_int, T_left_bot_int
1848  real :: dT_top, dT_bottom, dT_layer, dT_ave, dT_sublayer, dT_top_int, dT_bot_int
1849  real, dimension(nk+1) :: Til !< Left-column interface tracer (conc, e.g. degC)
1850  real, dimension(nk+1) :: Tir !< Right-column interface tracer (conc, e.g. degC)
1851  real, dimension(nk) :: aL_l !< Left-column left edge value of tracer (conc, e.g. degC)
1852  real, dimension(nk) :: aR_l !< Left-column right edge value of tracer (conc, e.g. degC)
1853  real, dimension(nk) :: aL_r !< Right-column left edge value of tracer (conc, e.g. degC)
1854  real, dimension(nk) :: aR_r !< Right-column right edge value of tracer (conc, e.g. degC)
1855  ! Discontinuous reconstruction
1856  integer :: iMethod
1857  real, dimension(nk,2) :: Tid_l !< Left-column interface tracer (conc, e.g. degC)
1858  real, dimension(nk,2) :: Tid_r !< Right-column interface tracer (conc, e.g. degC)
1859  real, dimension(nk,deg+1) :: ppoly_r_coeffs_l
1860  real, dimension(nk,deg+1) :: ppoly_r_coeffs_r
1861  real, dimension(nk,deg+1) :: ppoly_r_S_l
1862  real, dimension(nk,deg+1) :: ppoly_r_S_r
1863  logical :: down_flux
1864  ! Setup reconstruction edge values
1865  if (continuous) then
1866  call interface_scalar(nk, hl, tl, til, 2, h_neglect)
1867  call interface_scalar(nk, hr, tr, tir, 2, h_neglect)
1868  call ppm_left_right_edge_values(nk, tl, til, al_l, ar_l)
1869  call ppm_left_right_edge_values(nk, tr, tir, al_r, ar_r)
1870  else
1871  ppoly_r_coeffs_l(:,:) = 0.
1872  ppoly_r_coeffs_r(:,:) = 0.
1873  tid_l(:,:) = 0.
1874  tid_r(:,:) = 0.
1875 
1876  call build_reconstructions_1d( remap_cs, nk, hl, tl, ppoly_r_coeffs_l, tid_l, &
1877  ppoly_r_s_l, imethod, h_neglect, h_neglect_edge )
1878  call build_reconstructions_1d( remap_cs, nk, hr, tr, ppoly_r_coeffs_r, tid_r, &
1879  ppoly_r_s_r, imethod, h_neglect, h_neglect_edge )
1880  endif
1881 
1882  do k_sublayer = 1, nsurf-1
1883  if (heff(k_sublayer) == 0.) then
1884  flx(k_sublayer) = 0.
1885  else
1886  if (continuous) then
1887  klb = kol(k_sublayer+1)
1888  t_left_bottom = ( 1. - pil(k_sublayer+1) ) * til(klb) + pil(k_sublayer+1) * til(klb+1)
1889  klt = kol(k_sublayer)
1890  t_left_top = ( 1. - pil(k_sublayer) ) * til(klt) + pil(k_sublayer) * til(klt+1)
1891  t_left_layer = ppm_ave(pil(k_sublayer), pil(k_sublayer+1) + real(klb-klt), &
1892  al_l(klt), ar_l(klt), tl(klt))
1893 
1894  krb = kor(k_sublayer+1)
1895  t_right_bottom = ( 1. - pir(k_sublayer+1) ) * tir(krb) + pir(k_sublayer+1) * tir(krb+1)
1896  krt = kor(k_sublayer)
1897  t_right_top = ( 1. - pir(k_sublayer) ) * tir(krt) + pir(k_sublayer) * tir(krt+1)
1898  t_right_layer = ppm_ave(pir(k_sublayer), pir(k_sublayer+1) + real(krb-krt), &
1899  al_r(krt), ar_r(krt), tr(krt))
1900  dt_top = t_right_top - t_left_top
1901  dt_bottom = t_right_bottom - t_left_bottom
1902  dt_ave = 0.5 * ( dt_top + dt_bottom )
1903  dt_layer = t_right_layer - t_left_layer
1904  if (signum(1.,dt_top) * signum(1.,dt_bottom) <= 0. .or. signum(1.,dt_ave) * signum(1.,dt_layer) <= 0.) then
1905  dt_ave = 0.
1906  else
1907  dt_ave = dt_layer
1908  endif
1909  flx(k_sublayer) = dt_ave * heff(k_sublayer)
1910  else ! Discontinuous reconstruction
1911  ! Calculate tracer values on left and right side of the neutral surface
1912  call neutral_surface_t_eval(nk, nsurf, k_sublayer, kol, pil, tl, tid_l, deg, imethod, &
1913  ppoly_r_coeffs_l, t_left_top, t_left_bottom, t_left_sub, &
1914  t_left_top_int, t_left_bot_int, t_left_layer)
1915  call neutral_surface_t_eval(nk, nsurf, k_sublayer, kor, pir, tr, tid_r, deg, imethod, &
1916  ppoly_r_coeffs_r, t_right_top, t_right_bottom, t_right_sub, &
1917  t_right_top_int, t_right_bot_int, t_right_layer)
1918 
1919  dt_top = t_right_top - t_left_top
1920  dt_bottom = t_right_bottom - t_left_bottom
1921  dt_sublayer = t_right_sub - t_left_sub
1922  dt_top_int = t_right_top_int - t_left_top_int
1923  dt_bot_int = t_right_bot_int - t_left_bot_int
1924  ! Enforcing the below criterion incorrectly zero out fluxes
1925  !dT_layer = T_right_layer - T_left_layer
1926 
1927  down_flux = dt_top <= 0. .and. dt_bottom <= 0. .and. &
1928  dt_sublayer <= 0. .and. dt_top_int <= 0. .and. &
1929  dt_bot_int <= 0.
1930  down_flux = down_flux .or. &
1931  (dt_top >= 0. .and. dt_bottom >= 0. .and. &
1932  dt_sublayer >= 0. .and. dt_top_int >= 0. .and. &
1933  dt_bot_int >= 0.)
1934  if (down_flux) then
1935  flx(k_sublayer) = dt_sublayer * heff(k_sublayer)
1936  else
1937  flx(k_sublayer) = 0.
1938  endif
1939  endif
1940  endif
1941  enddo
1942 

References interface_scalar(), neutral_surface_t_eval(), ppm_ave(), ppm_left_right_edge_values(), and signum().

Referenced by ndiff_unit_tests_continuous(), and neutral_diffusion().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ neutral_surface_t_eval()

subroutine mom_neutral_diffusion::neutral_surface_t_eval ( integer, intent(in)  nk,
integer, intent(in)  ns,
integer, intent(in)  k_sub,
integer, dimension(ns), intent(in)  Ks,
real, dimension(ns), intent(in)  Ps,
real, dimension(nk), intent(in)  T_mean,
real, dimension(nk,2), intent(in)  T_int,
integer, intent(in)  deg,
integer, intent(in)  iMethod,
real, dimension(nk,deg+1), intent(in)  T_poly,
real, intent(out)  T_top,
real, intent(out)  T_bot,
real, intent(out)  T_sub,
real, intent(out)  T_top_int,
real, intent(out)  T_bot_int,
real, intent(out)  T_layer 
)
private

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

Parameters
[in]nkNumber of cell everages
[in]nsNumber of neutral surfaces
[in]k_subIndex of current neutral layer
[in]ksList of the layers associated with each neutral surface
[in]psList of the positions within a layer of each surface
[in]t_meanCell average of tracer
[in]t_intCell interface values of tracer from reconstruction
[in]degDegree of reconstruction polynomial (e.g. 1 is linear)
[in]imethodMethod of integration to use
[in]t_polyCoefficients of polynomial reconstructions
[out]t_topTracer value at top (across discontinuity if necessary)
[out]t_botTracer value at bottom (across discontinuity if necessary)
[out]t_subAverage of the tracer value over the sublayer
[out]t_top_intTracer value at top interface of neutral layer
[out]t_bot_intTracer value at bottom interface of neutral layer
[out]t_layerCell-average that the the reconstruction belongs to

Definition at line 1948 of file MOM_neutral_diffusion.F90.

1948  integer, intent(in ) :: nk !< Number of cell everages
1949  integer, intent(in ) :: ns !< Number of neutral surfaces
1950  integer, intent(in ) :: k_sub !< Index of current neutral layer
1951  integer, dimension(ns), intent(in ) :: Ks !< List of the layers associated with each neutral surface
1952  real, dimension(ns), intent(in ) :: Ps !< List of the positions within a layer of each surface
1953  real, dimension(nk), intent(in ) :: T_mean !< Cell average of tracer
1954  real, dimension(nk,2), intent(in ) :: T_int !< Cell interface values of tracer from reconstruction
1955  integer, intent(in ) :: deg !< Degree of reconstruction polynomial (e.g. 1 is linear)
1956  integer, intent(in ) :: iMethod !< Method of integration to use
1957  real, dimension(nk,deg+1), intent(in ) :: T_poly !< Coefficients of polynomial reconstructions
1958  real, intent( out) :: T_top !< Tracer value at top (across discontinuity if necessary)
1959  real, intent( out) :: T_bot !< Tracer value at bottom (across discontinuity if necessary)
1960  real, intent( out) :: T_sub !< Average of the tracer value over the sublayer
1961  real, intent( out) :: T_top_int !< Tracer value at top interface of neutral layer
1962  real, intent( out) :: T_bot_int !< Tracer value at bottom interface of neutral layer
1963  real, intent( out) :: T_layer !< Cell-average that the the reconstruction belongs to
1964 
1965  integer :: kl, ks_top, ks_bot
1966 
1967  ks_top = k_sub
1968  ks_bot = k_sub + 1
1969  if ( ks(ks_top) /= ks(ks_bot) ) then
1970  call mom_error(fatal, "Neutral surfaces span more than one layer")
1971  endif
1972  kl = ks(k_sub)
1973  ! First if the neutral surfaces spans the entirety of a cell, then do not search across the discontinuity
1974  if ( (ps(ks_top) == 0.) .and. (ps(ks_bot) == 1.)) then
1975  t_top = t_int(kl,1)
1976  t_bot = t_int(kl,2)
1977  else
1978  ! Search across potential discontinuity at top
1979  if ( (kl > 1) .and. (ps(ks_top) == 0.) ) then
1980  t_top = t_int(kl-1,2)
1981  else
1982  t_top = evaluation_polynomial( t_poly(kl,:), deg+1, ps(ks_top) )
1983  endif
1984  ! Search across potential discontinuity at bottom
1985  if ( (kl < nk) .and. (ps(ks_bot) == 1.) ) then
1986  t_bot = t_int(kl+1,1)
1987  else
1988  t_bot = evaluation_polynomial( t_poly(kl,:), deg+1, ps(ks_bot) )
1989  endif
1990  endif
1991  t_sub = average_value_ppoly(nk, t_mean, t_int, t_poly, imethod, kl, ps(ks_top), ps(ks_bot))
1992  t_top_int = evaluation_polynomial( t_poly(kl,:), deg+1, ps(ks_top))
1993  t_bot_int = evaluation_polynomial( t_poly(kl,:), deg+1, ps(ks_bot))
1994  t_layer = t_mean(kl)
1995 

References mom_remapping::average_value_ppoly(), polynomial_functions::evaluation_polynomial(), and mom_error_handler::mom_error().

Referenced by neutral_surface_flux().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ plm_diff()

subroutine mom_neutral_diffusion::plm_diff ( integer, intent(in)  nk,
real, dimension(nk), intent(in)  h,
real, dimension(nk), intent(in)  S,
integer, intent(in)  c_method,
integer, intent(in)  b_method,
real, dimension(nk), intent(inout)  diff 
)
private

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
[in]nkNumber of levels
[in]hLayer thickness [H ~> m or kg m-2]
[in]sLayer salinity (conc, e.g. ppt)
[in]c_methodMethod to use for the centered difference
[in]b_method=1, use PCM in first/last cell, =2 uses linear extrapolation
[in,out]diffScalar difference across layer (conc, e.g. ppt) determined by the following values for c_method:
  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.

Definition at line 770 of file MOM_neutral_diffusion.F90.

770  integer, intent(in) :: nk !< Number of levels
771  real, dimension(nk), intent(in) :: h !< Layer thickness [H ~> m or kg m-2]
772  real, dimension(nk), intent(in) :: S !< Layer salinity (conc, e.g. ppt)
773  integer, intent(in) :: c_method !< Method to use for the centered difference
774  integer, intent(in) :: b_method !< =1, use PCM in first/last cell, =2 uses linear extrapolation
775  real, dimension(nk), intent(inout) :: diff !< Scalar difference across layer (conc, e.g. ppt)
776  !! determined by the following values for c_method:
777  !! 1. Second order finite difference (not recommended)
778  !! 2. Second order finite volume (used in original PPM)
779  !! 3. Finite-volume weighted least squares linear fit
780  !! \todo The use of c_method to choose a scheme is inefficient
781  !! and should eventually be moved up the call tree.
782 
783  ! Local variables
784  integer :: k
785  real :: hkm1, hk, hkp1, Skm1, Sk, Skp1, diff_l, diff_r, diff_c
786 
787  do k = 2, nk-1
788  hkm1 = h(k-1)
789  hk = h(k)
790  hkp1 = h(k+1)
791 
792  if ( ( hkp1 + hk ) * ( hkm1 + hk ) > 0.) then
793  skm1 = s(k-1)
794  sk = s(k)
795  skp1 = s(k+1)
796  if (c_method==1) then
797  ! Simple centered diff (from White)
798  if ( hk + 0.5 * (hkm1 + hkp1) /= 0. ) then
799  diff_c = ( skp1 - skm1 ) * ( hk / ( hk + 0.5 * (hkm1 + hkp1) ) )
800  else
801  diff_c = 0.
802  endif
803  elseif (c_method==2) then
804  ! Second order accurate centered FV slope (from Colella and Woodward, JCP 1984)
805  diff_c = fv_diff(hkm1, hk, hkp1, skm1, sk, skp1)
806  elseif (c_method==3) then
807  ! Second order accurate finite-volume least squares slope
808  diff_c = hk * fvlsq_slope(hkm1, hk, hkp1, skm1, sk, skp1)
809  endif
810  ! Limit centered slope by twice the side differenced slopes
811  diff_l = 2. * ( sk - skm1 )
812  diff_r = 2. * ( skp1 - sk )
813  if ( signum(1., diff_l) * signum(1., diff_r) <= 0. ) then
814  diff(k) = 0. ! PCM for local extrema
815  else
816  diff(k) = sign( min( abs(diff_l), abs(diff_c), abs(diff_r) ), diff_c )
817  endif
818  else
819  diff(k) = 0. ! PCM next to vanished layers
820  endif
821  enddo
822  if (b_method==1) then ! PCM for top and bottom layer
823  diff(1) = 0.
824  diff(nk) = 0.
825  elseif (b_method==2) then ! Linear extrapolation for top and bottom interfaces
826  diff(1) = ( s(2) - s(1) ) * 2. * ( h(1) / ( h(1) + h(2) ) )
827  diff(nk) = s(nk) - s(nk-1) * 2. * ( h(nk) / ( h(nk-1) + h(nk) ) )
828  endif
829 

References fv_diff(), fvlsq_slope(), and signum().

Referenced by interface_scalar().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ ppm_ave()

real function mom_neutral_diffusion::ppm_ave ( real, intent(in)  xL,
real, intent(in)  xR,
real, intent(in)  aL,
real, intent(in)  aR,
real, intent(in)  aMean 
)
private

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

Parameters
[in]xlFraction position of left bound (0,1)
[in]xrFraction position of right bound (0,1)
[in]alLeft edge scalar value, at x=0
[in]arRight edge scalar value, at x=1
[in]ameanAverage scalar value of cell

Definition at line 732 of file MOM_neutral_diffusion.F90.

732  real, intent(in) :: xL !< Fraction position of left bound (0,1)
733  real, intent(in) :: xR !< Fraction position of right bound (0,1)
734  real, intent(in) :: aL !< Left edge scalar value, at x=0
735  real, intent(in) :: aR !< Right edge scalar value, at x=1
736  real, intent(in) :: aMean !< Average scalar value of cell
737 
738  ! Local variables
739  real :: dx, xave, a6, a6o3
740 
741  dx = xr - xl
742  xave = 0.5 * ( xr + xl )
743  a6o3 = 2. * amean - ( al + ar ) ! a6 / 3.
744  a6 = 3. * a6o3
745 
746  if (dx<0.) then
747  stop 'ppm_ave: dx<0 should not happend!'
748  elseif (dx>1.) then
749  stop 'ppm_ave: dx>1 should not happend!'
750  elseif (dx==0.) then
751  ppm_ave = al + ( ar - al ) * xr + a6 * xr * ( 1. - xr )
752  else
753  ppm_ave = ( al + xave * ( ( ar - al ) + a6 ) ) - a6o3 * ( xr**2 + xr * xl + xl**2 )
754  endif

Referenced by neutral_surface_flux().

Here is the caller graph for this function:

◆ ppm_edge()

real function mom_neutral_diffusion::ppm_edge ( real, intent(in)  hkm1,
real, intent(in)  hk,
real, intent(in)  hkp1,
real, intent(in)  hkp2,
real, intent(in)  Ak,
real, intent(in)  Akp1,
real, intent(in)  Pk,
real, intent(in)  Pkp1,
real, intent(in)  h_neglect 
)
private

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
[in]hkm1Width of cell k-1
[in]hkWidth of cell k
[in]hkp1Width of cell k+1
[in]hkp2Width of cell k+2
[in]akAverage scalar value of cell k
[in]akp1Average scalar value of cell k+1
[in]pkPLM slope for cell k
[in]pkp1PLM slope for cell k+1
[in]h_neglectA negligibly small thickness [H ~> m or kg m-2]

Definition at line 692 of file MOM_neutral_diffusion.F90.

692  real, intent(in) :: hkm1 !< Width of cell k-1
693  real, intent(in) :: hk !< Width of cell k
694  real, intent(in) :: hkp1 !< Width of cell k+1
695  real, intent(in) :: hkp2 !< Width of cell k+2
696  real, intent(in) :: Ak !< Average scalar value of cell k
697  real, intent(in) :: Akp1 !< Average scalar value of cell k+1
698  real, intent(in) :: Pk !< PLM slope for cell k
699  real, intent(in) :: Pkp1 !< PLM slope for cell k+1
700  real, intent(in) :: h_neglect !< A negligibly small thickness [H ~> m or kg m-2]
701 
702  ! Local variables
703  real :: R_hk_hkp1, R_2hk_hkp1, R_hk_2hkp1, f1, f2, f3, f4
704 
705  r_hk_hkp1 = hk + hkp1
706  if (r_hk_hkp1 <= 0.) then
707  ppm_edge = 0.5 * ( ak + akp1 )
708  return
709  endif
710  r_hk_hkp1 = 1. / r_hk_hkp1
711  if (hk<hkp1) then
712  ppm_edge = ak + ( hk * r_hk_hkp1 ) * ( akp1 - ak )
713  else
714  ppm_edge = akp1 + ( hkp1 * r_hk_hkp1 ) * ( ak - akp1 )
715  endif
716 
717  r_2hk_hkp1 = 1. / ( ( 2. * hk + hkp1 ) + h_neglect )
718  r_hk_2hkp1 = 1. / ( ( hk + 2. * hkp1 ) + h_neglect )
719  f1 = 1./ ( ( hk + hkp1) + ( hkm1 + hkp2 ) )
720  f2 = 2. * ( hkp1 * hk ) * r_hk_hkp1 * &
721  ( ( hkm1 + hk ) * r_2hk_hkp1 - ( hkp2 + hkp1 ) * r_hk_2hkp1 )
722  f3 = hk * ( hkm1 + hk ) * r_2hk_hkp1
723  f4 = hkp1 * ( hkp1 + hkp2 ) * r_hk_2hkp1
724 
725  ppm_edge = ppm_edge + f1 * ( f2 * ( akp1 - ak ) - ( f3 * pkp1 - f4 * pk ) )
726 

Referenced by interface_scalar().

Here is the caller graph for this function:

◆ ppm_left_right_edge_values()

subroutine mom_neutral_diffusion::ppm_left_right_edge_values ( integer, intent(in)  nk,
real, dimension(nk), intent(in)  Tl,
real, dimension(nk+1), intent(in)  Ti,
real, dimension(nk), intent(inout)  aL,
real, dimension(nk), intent(inout)  aR 
)
private

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

Parameters
[in]nkNumber of levels
[in]tlLayer tracer (conc, e.g. degC)
[in]tiInterface tracer (conc, e.g. degC)
[in,out]alLeft edge value of tracer (conc, e.g. degC)
[in,out]arRight edge value of tracer (conc, e.g. degC)

Definition at line 2000 of file MOM_neutral_diffusion.F90.

2000  integer, intent(in) :: nk !< Number of levels
2001  real, dimension(nk), intent(in) :: Tl !< Layer tracer (conc, e.g. degC)
2002  real, dimension(nk+1), intent(in) :: Ti !< Interface tracer (conc, e.g. degC)
2003  real, dimension(nk), intent(inout) :: aL !< Left edge value of tracer (conc, e.g. degC)
2004  real, dimension(nk), intent(inout) :: aR !< Right edge value of tracer (conc, e.g. degC)
2005 
2006  integer :: k
2007  ! Setup reconstruction edge values
2008  do k = 1, nk
2009  al(k) = ti(k)
2010  ar(k) = ti(k+1)
2011  if ( signum(1., ar(k) - tl(k))*signum(1., tl(k) - al(k)) <= 0.0 ) then
2012  al(k) = tl(k)
2013  ar(k) = tl(k)
2014  elseif ( sign(3., ar(k) - al(k)) * ( (tl(k) - al(k)) + (tl(k) - ar(k))) > abs(ar(k) - al(k)) ) then
2015  al(k) = tl(k) + 2.0 * ( tl(k) - ar(k) )
2016  elseif ( sign(3., ar(k) - al(k)) * ( (tl(k) - al(k)) + (tl(k) - ar(k))) < -abs(ar(k) - al(k)) ) then
2017  ar(k) = tl(k) + 2.0 * ( tl(k) - al(k) )
2018  endif
2019  enddo

References signum().

Referenced by neutral_surface_flux().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ search_other_column()

real function mom_neutral_diffusion::search_other_column ( type(neutral_diffusion_cs), intent(in)  CS,
integer, intent(in)  ksurf,
real, intent(in)  pos_last,
real, intent(in)  T_from,
real, intent(in)  S_from,
real, intent(in)  P_from,
real, intent(in)  T_top,
real, intent(in)  S_top,
real, intent(in)  P_top,
real, intent(in)  T_bot,
real, intent(in)  S_bot,
real, intent(in)  P_bot,
real, dimension(:), intent(in)  T_poly,
real, dimension(:), intent(in)  S_poly 
)
private

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

Parameters
[in]csNeutral diffusion control structure
[in]ksurfCurrent index of neutral surface
[in]pos_lastLast position within the current layer, used as the lower bound in the rootfinding algorithm
[in]t_fromTemperature at the searched from interface
[in]s_fromSalinity at the searched from interface
[in]p_fromPressure at the searched from interface
[in]t_topTemperature at the searched to top interface
[in]s_topSalinity at the searched to top interface
[in]p_topPressure at the searched to top interface
[in]t_botTemperature at the searched to bottom interface
[in]s_botSalinity at the searched to bottom interface
[in]p_botPressure at the searched to bottom interface
[in]t_polyTemperature polynomial reconstruction coefficients
[in]s_polySalinity polynomial reconstruction coefficients

Definition at line 1393 of file MOM_neutral_diffusion.F90.

1393  type(neutral_diffusion_CS), intent(in ) :: CS !< Neutral diffusion control structure
1394  integer, intent(in ) :: ksurf !< Current index of neutral surface
1395  real, intent(in ) :: pos_last !< Last position within the current layer, used as the lower
1396  !! bound in the rootfinding algorithm
1397  real, intent(in ) :: T_from !< Temperature at the searched from interface
1398  real, intent(in ) :: S_from !< Salinity at the searched from interface
1399  real, intent(in ) :: P_from !< Pressure at the searched from interface
1400  real, intent(in ) :: T_top !< Temperature at the searched to top interface
1401  real, intent(in ) :: S_top !< Salinity at the searched to top interface
1402  real, intent(in ) :: P_top !< Pressure at the searched to top interface
1403  real, intent(in ) :: T_bot !< Temperature at the searched to bottom interface
1404  real, intent(in ) :: S_bot !< Salinity at the searched to bottom interface
1405  real, intent(in ) :: P_bot !< Pressure at the searched to bottom interface
1406  real, dimension(:), intent(in ) :: T_poly !< Temperature polynomial reconstruction coefficients
1407  real, dimension(:), intent(in ) :: S_poly !< Salinity polynomial reconstruction coefficients
1408  ! Local variables
1409  real :: dRhotop, dRhobot
1410  real :: dRdT_top, dRdS_top, dRdT_bot, dRdS_bot
1411  real :: dRdT_from, dRdS_from
1412  real :: P_mid
1413 
1414  ! Calculate the differencei in density at the tops or the bottom
1415  if (cs%neutral_pos_method == 1 .or. cs%neutral_pos_method == 3) then
1416  call calc_delta_rho_and_derivs(cs, t_top, s_top, p_top, t_from, s_from, p_from, drhotop)
1417  call calc_delta_rho_and_derivs(cs, t_bot, s_bot, p_bot, t_from, s_from, p_from, drhobot)
1418  elseif (cs%neutral_pos_method == 2) then
1419  call calc_delta_rho_and_derivs(cs, t_top, s_top, p_top, t_from, s_from, p_from, drhotop, &
1420  drdt_top, drds_top, drdt_from, drds_from)
1421  call calc_delta_rho_and_derivs(cs, t_bot, s_bot, p_bot, t_from, s_from, p_from, drhobot, &
1422  drdt_bot, drds_bot, drdt_from, drds_from)
1423  endif
1424 
1425  ! Handle all the special cases EXCEPT if it connects within the layer
1426  if ( (drhotop > 0.) .or. (ksurf == 1) ) then ! First interface or lighter than anything in layer
1427  pos = pos_last
1428  elseif ( drhotop > drhobot ) then ! Unstably stratified
1429  pos = 1.
1430  elseif ( drhotop < 0. .and. drhobot < 0.) then ! Denser than anything in layer
1431  pos = 1.
1432  elseif ( drhotop == 0. .and. drhobot == 0. ) then ! Perfectly unstratified
1433  pos = 1.
1434  elseif ( drhobot == 0. ) then ! Matches perfectly at the Top
1435  pos = 1.
1436  elseif ( drhotop == 0. ) then ! Matches perfectly at the Bottom
1437  pos = pos_last
1438  else ! Neutral surface within layer
1439  pos = -1
1440  endif
1441 
1442  ! Can safely return if position is >= 0 otherwise will need to find the position within the layer
1443  if (pos>=0) return
1444 
1445  if (cs%neutral_pos_method==1) then
1446  pos = interpolate_for_nondim_position( drhotop, p_top, drhobot, p_bot )
1447  ! For the 'Linear' case of finding the neutral position, the fromerence pressure to use is the average
1448  ! of the midpoint of the layer being searched and the interface being searched from
1449  elseif (cs%neutral_pos_method == 2) then
1450  pos = find_neutral_pos_linear( cs, pos_last, t_from, s_from, p_from, drdt_from, drds_from, &
1451  p_top, drdt_top, drds_top, &
1452  p_bot, drdt_bot, drds_bot, t_poly, s_poly )
1453  elseif (cs%neutral_pos_method == 3) then
1454  pos = find_neutral_pos_full( cs, pos_last, t_from, s_from, p_from, p_top, p_bot, t_poly, s_poly)
1455  endif
1456 

References calc_delta_rho_and_derivs(), find_neutral_pos_full(), find_neutral_pos_linear(), and interpolate_for_nondim_position().

Referenced by find_neutral_surface_positions_discontinuous().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ signum()

real function mom_neutral_diffusion::signum ( real, intent(in)  a,
real, intent(in)  x 
)
private

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

Parameters
[in]aThe magnitude argument
[in]xThe sign (or zero) argument

Definition at line 759 of file MOM_neutral_diffusion.F90.

759  real, intent(in) :: a !< The magnitude argument
760  real, intent(in) :: x !< The sign (or zero) argument
761 
762  signum = sign(a,x)
763  if (x==0.) signum = 0.
764 

Referenced by neutral_surface_flux(), plm_diff(), and ppm_left_right_edge_values().

Here is the caller graph for this function:

◆ test_data1d()

logical function mom_neutral_diffusion::test_data1d ( logical, intent(in)  verbose,
integer, intent(in)  nk,
real, dimension(nk), intent(in)  Po,
real, dimension(nk), intent(in)  Ptrue,
character(len=*), intent(in)  title 
)
private

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

Parameters
[in]verboseIf true, write results to stdout
[in]nkNumber of layers
[in]poCalculated answer
[in]ptrueTrue answer
[in]titleTitle for messages

Definition at line 2652 of file MOM_neutral_diffusion.F90.

2652  logical, intent(in) :: verbose !< If true, write results to stdout
2653  integer, intent(in) :: nk !< Number of layers
2654  real, dimension(nk), intent(in) :: Po !< Calculated answer
2655  real, dimension(nk), intent(in) :: Ptrue !< True answer
2656  character(len=*), intent(in) :: title !< Title for messages
2657 
2658  ! Local variables
2659  integer :: k, stdunit
2660 
2661  test_data1d = .false.
2662  do k = 1,nk
2663  if (po(k) /= ptrue(k)) test_data1d = .true.
2664  enddo
2665 
2666  if (test_data1d .or. verbose) then
2667  stdunit = 6
2668  if (test_data1d) stdunit = 0 ! In case of wrong results, write to error stream
2669  write(stdunit,'(a)') title
2670  do k = 1,nk
2671  if (po(k) /= ptrue(k)) then
2672  test_data1d = .true.
2673  write(stdunit,'(a,i2,2(x,a,f20.16),x,a,1pe22.15,x,a)') &
2674  'k=',k,'Po=',po(k),'Ptrue=',ptrue(k),'err=',po(k)-ptrue(k),'WRONG!'
2675  else
2676  if (verbose) &
2677  write(stdunit,'(a,i2,2(x,a,f20.16),x,a,1pe22.15)') &
2678  'k=',k,'Po=',po(k),'Ptrue=',ptrue(k),'err=',po(k)-ptrue(k)
2679  endif
2680  enddo
2681  endif
2682 

Referenced by ndiff_unit_tests_continuous().

Here is the caller graph for this function:

◆ test_data1di()

logical function mom_neutral_diffusion::test_data1di ( logical, intent(in)  verbose,
integer, intent(in)  nk,
integer, dimension(nk), intent(in)  Po,
integer, dimension(nk), intent(in)  Ptrue,
character(len=*), intent(in)  title 
)
private

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

Parameters
[in]verboseIf true, write results to stdout
[in]nkNumber of layers
[in]poCalculated answer
[in]ptrueTrue answer
[in]titleTitle for messages

Definition at line 2687 of file MOM_neutral_diffusion.F90.

2687  logical, intent(in) :: verbose !< If true, write results to stdout
2688  integer, intent(in) :: nk !< Number of layers
2689  integer, dimension(nk), intent(in) :: Po !< Calculated answer
2690  integer, dimension(nk), intent(in) :: Ptrue !< True answer
2691  character(len=*), intent(in) :: title !< Title for messages
2692 
2693  ! Local variables
2694  integer :: k, stdunit
2695 
2696  test_data1di = .false.
2697  do k = 1,nk
2698  if (po(k) /= ptrue(k)) test_data1di = .true.
2699  enddo
2700 
2701  if (test_data1di .or. verbose) then
2702  stdunit = 6
2703  if (test_data1di) stdunit = 0 ! In case of wrong results, write to error stream
2704  write(stdunit,'(a)') title
2705  do k = 1,nk
2706  if (po(k) /= ptrue(k)) then
2707  test_data1di = .true.
2708  write(stdunit,'(a,i2,2(x,a,i5),x,a)') 'k=',k,'Io=',po(k),'Itrue=',ptrue(k),'WRONG!'
2709  else
2710  if (verbose) &
2711  write(stdunit,'(a,i2,2(x,a,i5))') 'k=',k,'Io=',po(k),'Itrue=',ptrue(k)
2712  endif
2713  enddo
2714  endif
2715 

◆ test_fv_diff()

logical function mom_neutral_diffusion::test_fv_diff ( logical, intent(in)  verbose,
real, intent(in)  hkm1,
real, intent(in)  hk,
real, intent(in)  hkp1,
real, intent(in)  Skm1,
real, intent(in)  Sk,
real, intent(in)  Skp1,
real, intent(in)  Ptrue,
character(len=*), intent(in)  title 
)
private

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

Parameters
[in]verboseIf true, write results to stdout
[in]hkm1Left cell width
[in]hkCenter cell width
[in]hkp1Right cell width
[in]skm1Left cell average value
[in]skCenter cell average value
[in]skp1Right cell average value
[in]ptrueTrue answer [Pa]
[in]titleTitle for messages

Definition at line 2556 of file MOM_neutral_diffusion.F90.

2556  logical, intent(in) :: verbose !< If true, write results to stdout
2557  real, intent(in) :: hkm1 !< Left cell width
2558  real, intent(in) :: hk !< Center cell width
2559  real, intent(in) :: hkp1 !< Right cell width
2560  real, intent(in) :: Skm1 !< Left cell average value
2561  real, intent(in) :: Sk !< Center cell average value
2562  real, intent(in) :: Skp1 !< Right cell average value
2563  real, intent(in) :: Ptrue !< True answer [Pa]
2564  character(len=*), intent(in) :: title !< Title for messages
2565 
2566  ! Local variables
2567  integer :: stdunit
2568  real :: Pret
2569 
2570  pret = fv_diff(hkm1, hk, hkp1, skm1, sk, skp1)
2571  test_fv_diff = (pret /= ptrue)
2572 
2573  if (test_fv_diff .or. verbose) then
2574  stdunit = 6
2575  if (test_fv_diff) stdunit = 0 ! In case of wrong results, write to error stream
2576  write(stdunit,'(a)') title
2577  if (test_fv_diff) then
2578  write(stdunit,'(2(x,a,f20.16),x,a)') 'pRet=',pret,'pTrue=',ptrue,'WRONG!'
2579  else
2580  write(stdunit,'(2(x,a,f20.16))') 'pRet=',pret,'pTrue=',ptrue
2581  endif
2582  endif
2583 

References fv_diff().

Referenced by ndiff_unit_tests_continuous().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ test_fvlsq_slope()

logical function mom_neutral_diffusion::test_fvlsq_slope ( logical, intent(in)  verbose,
real, intent(in)  hkm1,
real, intent(in)  hk,
real, intent(in)  hkp1,
real, intent(in)  Skm1,
real, intent(in)  Sk,
real, intent(in)  Skp1,
real, intent(in)  Ptrue,
character(len=*), intent(in)  title 
)
private

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

Parameters
[in]verboseIf true, write results to stdout
[in]hkm1Left cell width
[in]hkCenter cell width
[in]hkp1Right cell width
[in]skm1Left cell average value
[in]skCenter cell average value
[in]skp1Right cell average value
[in]ptrueTrue answer [Pa]
[in]titleTitle for messages

Definition at line 2588 of file MOM_neutral_diffusion.F90.

2588  logical, intent(in) :: verbose !< If true, write results to stdout
2589  real, intent(in) :: hkm1 !< Left cell width
2590  real, intent(in) :: hk !< Center cell width
2591  real, intent(in) :: hkp1 !< Right cell width
2592  real, intent(in) :: Skm1 !< Left cell average value
2593  real, intent(in) :: Sk !< Center cell average value
2594  real, intent(in) :: Skp1 !< Right cell average value
2595  real, intent(in) :: Ptrue !< True answer [Pa]
2596  character(len=*), intent(in) :: title !< Title for messages
2597 
2598  ! Local variables
2599  integer :: stdunit
2600  real :: Pret
2601 
2602  pret = fvlsq_slope(hkm1, hk, hkp1, skm1, sk, skp1)
2603  test_fvlsq_slope = (pret /= ptrue)
2604 
2605  if (test_fvlsq_slope .or. verbose) then
2606  stdunit = 6
2607  if (test_fvlsq_slope) stdunit = 0 ! In case of wrong results, write to error stream
2608  write(stdunit,'(a)') title
2609  if (test_fvlsq_slope) then
2610  write(stdunit,'(2(x,a,f20.16),x,a)') 'pRet=',pret,'pTrue=',ptrue,'WRONG!'
2611  else
2612  write(stdunit,'(2(x,a,f20.16))') 'pRet=',pret,'pTrue=',ptrue
2613  endif
2614  endif
2615 

References fvlsq_slope().

Referenced by ndiff_unit_tests_continuous().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ test_ifndp()

logical function mom_neutral_diffusion::test_ifndp ( logical, intent(in)  verbose,
real, intent(in)  rhoNeg,
real, intent(in)  Pneg,
real, intent(in)  rhoPos,
real, intent(in)  Ppos,
real, intent(in)  Ptrue,
character(len=*), intent(in)  title 
)
private

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

Parameters
[in]verboseIf true, write results to stdout
[in]rhonegLighter density [kg m-3]
[in]pnegInterface position of lighter density [Pa]
[in]rhoposHeavier density [kg m-3]
[in]pposInterface position of heavier density [Pa]
[in]ptrueTrue answer [Pa]
[in]titleTitle for messages

Definition at line 2620 of file MOM_neutral_diffusion.F90.

2620  logical, intent(in) :: verbose !< If true, write results to stdout
2621  real, intent(in) :: rhoNeg !< Lighter density [kg m-3]
2622  real, intent(in) :: Pneg !< Interface position of lighter density [Pa]
2623  real, intent(in) :: rhoPos !< Heavier density [kg m-3]
2624  real, intent(in) :: Ppos !< Interface position of heavier density [Pa]
2625  real, intent(in) :: Ptrue !< True answer [Pa]
2626  character(len=*), intent(in) :: title !< Title for messages
2627 
2628  ! Local variables
2629  integer :: stdunit
2630  real :: Pret
2631 
2632  pret = interpolate_for_nondim_position(rhoneg, pneg, rhopos, ppos)
2633  test_ifndp = (pret /= ptrue)
2634 
2635  if (test_ifndp .or. verbose) then
2636  stdunit = 6
2637  if (test_ifndp) stdunit = 0 ! In case of wrong results, write to error stream
2638  write(stdunit,'(a)') title
2639  if (test_ifndp) then
2640  write(stdunit,'(4(x,a,f20.16),2(x,a,1pe22.15),x,a)') &
2641  'r1=',rhoneg,'p1=',pneg,'r2=',rhopos,'p2=',ppos,'pRet=',pret,'pTrue=',ptrue,'WRONG!'
2642  else
2643  write(stdunit,'(4(x,a,f20.16),2(x,a,1pe22.15))') &
2644  'r1=',rhoneg,'p1=',pneg,'r2=',rhopos,'p2=',ppos,'pRet=',pret,'pTrue=',ptrue
2645  endif
2646  endif
2647 

References interpolate_for_nondim_position().

Referenced by ndiff_unit_tests_continuous().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ test_nsp()

logical function mom_neutral_diffusion::test_nsp ( logical, intent(in)  verbose,
integer, intent(in)  ns,
integer, dimension(ns), intent(in)  KoL,
integer, dimension(ns), intent(in)  KoR,
real, dimension(ns), intent(in)  pL,
real, dimension(ns), intent(in)  pR,
real, dimension(ns-1), intent(in)  hEff,
integer, dimension(ns), intent(in)  KoL0,
integer, dimension(ns), intent(in)  KoR0,
real, dimension(ns), intent(in)  pL0,
real, dimension(ns), intent(in)  pR0,
real, dimension(ns-1), intent(in)  hEff0,
character(len=*), intent(in)  title 
)
private

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

Parameters
[in]verboseIf true, write results to stdout
[in]nsNumber of surfaces
[in]kolIndex of first left interface above neutral surface
[in]korIndex of first right interface above neutral surface
[in]plFractional position of neutral surface within layer KoL of left column
[in]prFractional position of neutral surface within layer KoR of right column
[in]heffEffective thickness between two neutral surfaces [Pa]
[in]kol0Correct value for KoL
[in]kor0Correct value for KoR
[in]pl0Correct value for pL
[in]pr0Correct value for pR
[in]heff0Correct value for hEff
[in]titleTitle for messages

Definition at line 2721 of file MOM_neutral_diffusion.F90.

2721  logical, intent(in) :: verbose !< If true, write results to stdout
2722  integer, intent(in) :: ns !< Number of surfaces
2723  integer, dimension(ns), intent(in) :: KoL !< Index of first left interface above neutral surface
2724  integer, dimension(ns), intent(in) :: KoR !< Index of first right interface above neutral surface
2725  real, dimension(ns), intent(in) :: pL !< Fractional position of neutral surface within layer KoL of left column
2726  real, dimension(ns), intent(in) :: pR !< Fractional position of neutral surface within layer KoR of right column
2727  real, dimension(ns-1), intent(in) :: hEff !< Effective thickness between two neutral surfaces [Pa]
2728  integer, dimension(ns), intent(in) :: KoL0 !< Correct value for KoL
2729  integer, dimension(ns), intent(in) :: KoR0 !< Correct value for KoR
2730  real, dimension(ns), intent(in) :: pL0 !< Correct value for pL
2731  real, dimension(ns), intent(in) :: pR0 !< Correct value for pR
2732  real, dimension(ns-1), intent(in) :: hEff0 !< Correct value for hEff
2733  character(len=*), intent(in) :: title !< Title for messages
2734 
2735  ! Local variables
2736  integer :: k, stdunit
2737  logical :: this_row_failed
2738 
2739  test_nsp = .false.
2740  do k = 1,ns
2741  test_nsp = test_nsp .or. compare_nsp_row(kol(k), kor(k), pl(k), pr(k), kol0(k), kor0(k), pl0(k), pr0(k))
2742  if (k < ns) then
2743  if (heff(k) /= heff0(k)) test_nsp = .true.
2744  endif
2745  enddo
2746 
2747  if (test_nsp .or. verbose) then
2748  stdunit = 6
2749  if (test_nsp) stdunit = 0 ! In case of wrong results, write to error stream
2750  write(stdunit,'(a)') title
2751  do k = 1,ns
2752  this_row_failed = compare_nsp_row(kol(k), kor(k), pl(k), pr(k), kol0(k), kor0(k), pl0(k), pr0(k))
2753  if (this_row_failed) then
2754  write(stdunit,10) k,kol(k),pl(k),kor(k),pr(k),' <-- WRONG!'
2755  write(stdunit,10) k,kol0(k),pl0(k),kor0(k),pr0(k),' <-- should be this'
2756  else
2757  write(stdunit,10) k,kol(k),pl(k),kor(k),pr(k)
2758  endif
2759  if (k < ns) then
2760  if (heff(k) /= heff0(k)) then
2761  write(stdunit,'(i3,8x,"layer hEff =",2(f20.16,a))') k,heff(k)," .neq. ",heff0(k),' <-- WRONG!'
2762  else
2763  write(stdunit,'(i3,8x,"layer hEff =",f20.16)') k,heff(k)
2764  endif
2765  endif
2766  enddo
2767  endif
2768  if (test_nsp) call mom_error(fatal,"test_nsp failed")
2769 
2770 10 format("ks=",i3," kL=",i3," pL=",f20.16," kR=",i3," pR=",f20.16,a)

References compare_nsp_row(), and mom_error_handler::mom_error().

Referenced by ndiff_unit_tests_continuous(), and ndiff_unit_tests_discontinuous().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ test_rnp()

logical function mom_neutral_diffusion::test_rnp ( real, intent(in)  expected_pos,
real, intent(in)  test_pos,
character(len=*), intent(in)  title 
)
private

Compares output position from refine_nondim_position with an expected value.

Parameters
[in]expected_posThe expected position
[in]test_posThe position returned by the code
[in]titleA label for this test

Definition at line 2793 of file MOM_neutral_diffusion.F90.

2793  real, intent(in) :: expected_pos !< The expected position
2794  real, intent(in) :: test_pos !< The position returned by the code
2795  character(len=*), intent(in) :: title !< A label for this test
2796  ! Local variables
2797  integer :: stdunit = 6 ! Output to standard error
2798  test_rnp = abs(expected_pos - test_pos) > 2*epsilon(test_pos)
2799  if (test_rnp) then
2800  write(stdunit,'(A, f20.16, " .neq. ", f20.16, " <-- WRONG")') title, expected_pos, test_pos
2801  else
2802  write(stdunit,'(A, f20.16, " == ", f20.16)') title, expected_pos, test_pos
2803  endif

Referenced by ndiff_unit_tests_discontinuous().

Here is the caller graph for this function:

Variable Documentation

◆ mdl

character(len=40) mom_neutral_diffusion::mdl = "MOM_neutral_diffusion"

module name

Definition at line 103 of file MOM_neutral_diffusion.F90.

103 character(len=40) :: mdl = "MOM_neutral_diffusion" !< module name

Referenced by neutral_diffusion_init().