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, US, 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, p_surf)
	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, dRdT_ref, dRdS_ref, dRdT_top, dRdS_top, 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

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]	n	Number of levels
[in]	ns	Number of neutral surfaces
[in]	pint	Position of interfaces [R L2 T-2 ~> Pa] or other units
[in]	karr	Index of interface above position
[in]	nparr	Non-dimensional position within layer Karr(:) [nondim]
[in]	k_surface	k-interface to query

Returns

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

Definition at line 1837 of file MOM_neutral_diffusion.F90.

  integer, intent(in) :: n            !< Number of levels
  integer, intent(in) :: ns           !< Number of neutral surfaces
  real,    intent(in) :: Pint(n+1)    !< Position of interfaces [R L2 T-2 ~> Pa] or other units
  integer, intent(in) :: Karr(ns)     !< Index of interface above position
  real,    intent(in) :: NParr(ns)    !< Non-dimensional position within layer Karr(:) [nondim]
  integer, intent(in) :: k_surface    !< k-interface to query
  real                :: absolute_position !< The absolute position of a location [R L2 T-2 ~> Pa]
                                      !! or other units following Pint
  ! Local variables
  integer :: k

  k = karr(k_surface)
  if (k>n) stop 'absolute_position: k>nk is out of bounds!'
  absolute_position = pint(k) + nparr(k_surface) * ( pint(k+1) - pint(k) )

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]	n	Number of levels
[in]	ns	Number of neutral surfaces
[in]	pint	Position of interface [R L2 T-2 ~> Pa] or other units
[in]	karr	Indexes of interfaces about positions
[in]	nparr	Non-dimensional positions within layers Karr(:)

Returns

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

Definition at line 1856 of file MOM_neutral_diffusion.F90.

  integer, intent(in) :: n         !< Number of levels
  integer, intent(in) :: ns        !< Number of neutral surfaces
  real,    intent(in) :: Pint(n+1) !< Position of interface [R L2 T-2 ~> Pa] or other units
  integer, intent(in) :: Karr(ns)  !< Indexes of interfaces about positions
  real,    intent(in) :: NParr(ns) !< Non-dimensional positions within layers Karr(:)

  real,  dimension(ns) :: absolute_positions !< Absolute positions [R L2 T-2 ~> Pa]
                                   !! or other units following Pint

  ! Local variables
  integer :: k_surface, k

  do k_surface = 1, ns
    absolute_positions(k_surface) = absolute_position(n,ns,pint,karr,nparr,k_surface)
  enddo

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

	cs	Neutral diffusion control structure
[in]	t1	Temperature at point 1 [degC]
[in]	s1	Salinity at point 1 [ppt]
[in]	p1_in	Pressure at point 1 [R L2 T-2 ~> Pa]
[in]	t2	Temperature at point 2 [degC]
[in]	s2	Salinity at point 2 [ppt]
[in]	p2_in	Pressure at point 2 [R L2 T-2 ~> Pa]
[out]	drho	Difference in density between the two points [R ~> kg m-3]
[out]	drdt1_out	drho_dt at point 1 [R degC-1 ~> kg m-3 degC-1]
[out]	drds1_out	drho_ds at point 1 [R ppt-1 ~> kg m-3 ppt-1]
[out]	drdt2_out	drho_dt at point 2 [R degC-1 ~> kg m-3 degC-1]
[out]	drds2_out	drho_ds at point 2 [R ppt-1 ~> kg m-3 ppt-1]

Definition at line 1756 of file MOM_neutral_diffusion.F90.

  type(neutral_diffusion_CS)    :: CS        !< Neutral diffusion control structure
  real,           intent(in   ) :: T1        !< Temperature at point 1 [degC]
  real,           intent(in   ) :: S1        !< Salinity at point 1 [ppt]
  real,           intent(in   ) :: p1_in     !< Pressure at point 1 [R L2 T-2 ~> Pa]
  real,           intent(in   ) :: T2        !< Temperature at point 2 [degC]
  real,           intent(in   ) :: S2        !< Salinity at point 2 [ppt]
  real,           intent(in   ) :: p2_in     !< Pressure at point 2 [R L2 T-2 ~> Pa]
  real,           intent(  out) :: drho      !< Difference in density between the two points [R ~> kg m-3]
  real, optional, intent(  out) :: dRdT1_out !< drho_dt at point 1 [R degC-1 ~> kg m-3 degC-1]
  real, optional, intent(  out) :: dRdS1_out !< drho_ds at point 1 [R ppt-1 ~> kg m-3 ppt-1]
  real, optional, intent(  out) :: dRdT2_out !< drho_dt at point 2 [R degC-1 ~> kg m-3 degC-1]
  real, optional, intent(  out) :: dRdS2_out !< drho_ds at point 2 [R ppt-1 ~> kg m-3 ppt-1]
  ! Local variables
  real :: rho1, rho2   ! Densities [R ~> kg m-3]
  real :: p1, p2, pmid ! Pressures [R L2 T-2 ~> Pa]
  real :: drdt1, drdt2 ! Partial derivatives of density with temperature [R degC-1 ~> kg m-3 degC-1]
  real :: drds1, drds2 ! Partial derivatives of density with salinity [R ppt-1 ~> kg m-3 ppt-1]
  real :: drdp1, drdp2 ! Partial derivatives of density with pressure [T2 L-2 ~> s2 m-2]

  ! Use the same reference pressure or the in-situ pressure
  if (cs%ref_pres > 0.) then
    p1 = cs%ref_pres
    p2 = cs%ref_pres
  else
    p1 = p1_in
    p2 = p2_in
  endif

  ! Use the full linear equation of state to calculate the difference in density (expensive!)
  if     (trim(cs%delta_rho_form) == 'full') then
    pmid = 0.5 * (p1 + p2)
    call calculate_density( t1, s1, pmid, rho1, cs%EOS)
    call calculate_density( t2, s2, pmid, rho2, cs%EOS)
    drho = rho1 - rho2
  ! Use the density derivatives at the average of pressures and the differentces int temperature
  elseif (trim(cs%delta_rho_form) == 'mid_pressure') then
    pmid = 0.5 * (p1 + p2)
    if (cs%ref_pres>=0) pmid = cs%ref_pres
    call calculate_density_derivs(t1, s1, pmid, drdt1, drds1, cs%EOS)
    call calculate_density_derivs(t2, s2, pmid, drdt2, drds2, cs%EOS)
    drho = delta_rho_from_derivs( t1, s1, p1, drdt1, drds1, t2, s2, p2, drdt2, drds2)
  elseif (trim(cs%delta_rho_form) == 'local_pressure') then
    call calculate_density_derivs(t1, s1, p1, drdt1, drds1, cs%EOS)
    call calculate_density_derivs(t2, s2, p2, drdt2, drds2, cs%EOS)
    drho = delta_rho_from_derivs( t1, s1, p1, drdt1, drds1, t2, s2, p2, drdt2, drds2)
  else
    call mom_error(fatal, "delta_rho_form is not recognized")
  endif

  if (PRESENT(drdt1_out)) drdt1_out = drdt1
  if (PRESENT(drds1_out)) drds1_out = drds1
  if (PRESENT(drdt2_out)) drdt2_out = drdt2
  if (PRESENT(drds2_out)) drds2_out = drds2

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]	kol	Index of first left interface above neutral surface
[in]	kor	Index of first right interface above neutral surface
[in]	pl	Fractional position of neutral surface within layer KoL of left column
[in]	pr	Fractional position of neutral surface within layer KoR of right column
[in]	kol0	Correct value for KoL
[in]	kor0	Correct value for KoR
[in]	pl0	Correct value for pL
[in]	pr0	Correct value for pR

Definition at line 2829 of file MOM_neutral_diffusion.F90.

  integer,  intent(in) :: KoL   !< Index of first left interface above neutral surface
  integer,  intent(in) :: KoR   !< Index of first right interface above neutral surface
  real,     intent(in) :: pL    !< Fractional position of neutral surface within layer KoL of left column
  real,     intent(in) :: pR    !< Fractional position of neutral surface within layer KoR of right column
  integer,  intent(in) :: KoL0  !< Correct value for KoL
  integer,  intent(in) :: KoR0  !< Correct value for KoR
  real,     intent(in) :: pL0   !< Correct value for pL
  real,     intent(in) :: pR0   !< Correct value for pR

  compare_nsp_row = .false.
  if (kol /= kol0) compare_nsp_row = .true.
  if (kor /= kor0) compare_nsp_row = .true.
  if (pl /= pl0) compare_nsp_row = .true.
  if (pr /= pr0) compare_nsp_row = .true.

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

t1	Temperature at point 1 [degC]
s1	Salinity at point 1 [ppt]
p1	Pressure at point 1 [R L2 T-2 ~> Pa]
drdt1	The partial derivative of density with temperature at point 1 [R degC-1 ~> kg m-3 degC-1]
drds1	The partial derivative of density with salinity at point 1 [R ppt-1 ~> kg m-3 ppt-1]
t2	Temperature at point 2 [degC]
s2	Salinity at point 2 [ppt]
p2	Pressure at point 2 [R L2 T-2 ~> Pa]
drdt2	The partial derivative of density with temperature at point 2 [R degC-1 ~> kg m-3 degC-1]
drds2	The partial derivative of density with salinity at point 2 [R ppt-1 ~> kg m-3 ppt-1]

Definition at line 1818 of file MOM_neutral_diffusion.F90.

  real :: T1    !< Temperature at point 1 [degC]
  real :: S1    !< Salinity at point 1 [ppt]
  real :: P1    !< Pressure at point 1 [R L2 T-2 ~> Pa]
  real :: dRdT1 !< The partial derivative of density with temperature at point 1 [R degC-1 ~> kg m-3 degC-1]
  real :: dRdS1 !< The partial derivative of density with salinity at point 1 [R ppt-1 ~> kg m-3 ppt-1]
  real :: T2    !< Temperature at point 2 [degC]
  real :: S2    !< Salinity at point 2 [ppt]
  real :: P2    !< Pressure at point 2 [R L2 T-2 ~> Pa]
  real :: dRdT2 !< The partial derivative of density with temperature at point 2 [R degC-1 ~> kg m-3 degC-1]
  real :: dRdS2 !< The partial derivative of density with salinity at point 2 [R ppt-1 ~> kg m-3 ppt-1]
  ! Local variables
  real :: drho  ! The density difference [R ~> kg m-3]

  drho = 0.5 * ( (drdt1+drdt2)*(t1-t2) + (drds1+drds2)*(s1-s2))

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]	cs	Control structure with parameters for this module
[in]	z0	Lower bound of position, also serves as the initial guess [nondim]
[in]	t_ref	Temperature at the searched from interface [degC]
[in]	s_ref	Salinity at the searched from interface [ppt]
[in]	p_ref	Pressure at the searched from interface [R L2 T-2 ~> Pa]
[in]	p_top	Pressure at top of layer being searched [R L2 T-2 ~> Pa]
[in]	p_bot	Pressure at bottom of layer being searched [R L2 T-2 ~> Pa]
[in]	ppoly_t	Coefficients of the polynomial reconstruction of T within the layer to be searched [degC]
[in]	ppoly_s	Coefficients of the polynomial reconstruction of T within the layer to be searched [ppt]

Returns

Position where drho = 0 [nondim]

Definition at line 1662 of file MOM_neutral_diffusion.F90.

  type(neutral_diffusion_CS),intent(in) :: CS        !< Control structure with parameters for this module
  real,                      intent(in) :: z0        !< Lower bound of position, also serves as the
                                                     !! initial guess [nondim]
  real,                      intent(in) :: T_ref     !< Temperature at the searched from interface [degC]
  real,                      intent(in) :: S_ref     !< Salinity at the searched from interface [ppt]
  real,                      intent(in) :: P_ref     !< Pressure at the searched from interface [R L2 T-2 ~> Pa]
  real,                      intent(in) :: P_top     !< Pressure at top of layer being searched [R L2 T-2 ~> Pa]
  real,                      intent(in) :: P_bot     !< Pressure at bottom of layer being searched [R L2 T-2 ~> Pa]
  real, dimension(:),        intent(in) :: ppoly_T   !< Coefficients of the polynomial reconstruction of T within
                                                     !! the layer to be searched [degC]
  real, dimension(:),        intent(in) :: ppoly_S   !< Coefficients of the polynomial reconstruction of T within
                                                     !! the layer to be searched [ppt]
  real                                  :: z         !< Position where drho = 0 [nondim]
  ! Local variables
  integer :: iter
  integer :: nterm

  real :: drho_a, drho_b, drho_c ! Density differences [R ~> kg m-3]
  real :: a, b, c     ! Fractional positions [nondim]
  real :: Ta, Tb, Tc  ! Temperatures [degC]
  real :: Sa, Sb, Sc  ! Salinities [ppt]
  real :: Pa, Pb, Pc  ! Pressures [R L2 T-2 ~> Pa]
  integer :: side

  side = 0
  ! Set the first two evaluation to the endpoints of the interval
  b = z0 ; c = 1
  nterm = SIZE(ppoly_t)

  ! Calculate drho at the minimum bound
  tb = evaluation_polynomial( ppoly_t, nterm, b )
  sb = evaluation_polynomial( ppoly_s, nterm, b )
  pb = p_top*(1.-b) + p_bot*b
  call calc_delta_rho_and_derivs(cs, tb, sb, pb, t_ref, s_ref, p_ref, drho_b)

  ! Calculate drho at the maximum bound
  tc = evaluation_polynomial( ppoly_t, nterm, 1. )
  sc = evaluation_polynomial( ppoly_s, nterm, 1. )
  pc = p_bot
  call calc_delta_rho_and_derivs(cs, tc, sc, pc, t_ref, s_ref, p_ref, drho_c)

  if (drho_b >= 0.) then
    z = z0
    return
  elseif (drho_c == 0.) then
    z = 1.
    return
  endif
  if ( sign(1.,drho_b) == sign(1.,drho_c) ) then
    z = z0
    return
  endif

  do iter = 1, cs%max_iter
    ! Calculate new position and evaluate if we have converged
    a = (drho_b*c - drho_c*b)/(drho_b-drho_c)
    ta = evaluation_polynomial( ppoly_t, nterm, a )
    sa = evaluation_polynomial( ppoly_s, nterm, a )
    pa = p_top*(1.-a) + p_bot*a
    call calc_delta_rho_and_derivs(cs, ta, sa, pa, t_ref, s_ref, p_ref, drho_a)
    if (abs(drho_a) < cs%drho_tol) then
      z = a
      return
    endif

    if (drho_a*drho_c > 0.) then
      if ( abs(a-c)<cs%x_tol) then
        z = a
        return
      endif
      c = a ; drho_c = drho_a;
      if (side == -1) drho_b = 0.5*drho_b
      side = -1
    elseif ( drho_b*drho_a > 0 ) then
      if ( abs(a-b)<cs%x_tol) then
        z = a
        return
      endif
      b = a ; drho_b = drho_a
      if (side == 1) drho_c = 0.5*drho_c
      side = 1
    else
      z = a
      return
    endif
  enddo

  z = a

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)  dRdT_ref, real, intent(in)  dRdS_ref, real, intent(in)  dRdT_top, real, intent(in)  dRdS_top, 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]	cs	Control structure with parameters for this module
[in]	z0	Lower bound of position, also serves as the initial guess [nondim]
[in]	t_ref	Temperature at the searched from interface [degC]
[in]	s_ref	Salinity at the searched from interface [ppt]
[in]	drdt_ref	dRho/dT at the searched from interface [R degC-1 ~> kg m-3 degC-1]
[in]	drds_ref	dRho/dS at the searched from interface [R ppt-1 ~> kg m-3 ppt-1]
[in]	drdt_top	dRho/dT at top of layer being searched [R degC-1 ~> kg m-3 degC-1]
[in]	drds_top	dRho/dS at top of layer being searched [R ppt-1 ~> kg m-3 ppt-1]
[in]	drdt_bot	dRho/dT at bottom of layer being searched [R degC-1 ~> kg m-3 degC-1]
[in]	drds_bot	dRho/dS at bottom of layer being searched [R ppt-1 ~> kg m-3 ppt-1]
[in]	ppoly_t	Coefficients of the polynomial reconstruction of T within the layer to be searched [degC].
[in]	ppoly_s	Coefficients of the polynomial reconstruction of S within the layer to be searched [ppt].

Returns

Position where drho = 0 [nondim]

Definition at line 1542 of file MOM_neutral_diffusion.F90.

  type(neutral_diffusion_CS),intent(in) :: CS        !< Control structure with parameters for this module
  real,                      intent(in) :: z0        !< Lower bound of position, also serves as the
                                                     !! initial guess [nondim]
  real,                      intent(in) :: T_ref     !< Temperature at the searched from interface [degC]
  real,                      intent(in) :: S_ref     !< Salinity at the searched from interface [ppt]
  real,                      intent(in) :: dRdT_ref  !< dRho/dT at the searched from interface
                                                     !! [R degC-1 ~> kg m-3 degC-1]
  real,                      intent(in) :: dRdS_ref  !< dRho/dS at the searched from interface
                                                     !! [R ppt-1 ~> kg m-3 ppt-1]
  real,                      intent(in) :: dRdT_top  !< dRho/dT at top of layer being searched
                                                     !! [R degC-1 ~> kg m-3 degC-1]
  real,                      intent(in) :: dRdS_top  !< dRho/dS at top of layer being searched
                                                     !! [R ppt-1 ~> kg m-3 ppt-1]
  real,                      intent(in) :: dRdT_bot  !< dRho/dT at bottom of layer being searched
                                                     !! [R degC-1 ~> kg m-3 degC-1]
  real,                      intent(in) :: dRdS_bot  !< dRho/dS at bottom of layer being searched
                                                     !! [R ppt-1 ~> kg m-3 ppt-1]
  real, dimension(:),        intent(in) :: ppoly_T   !< Coefficients of the polynomial reconstruction of T within
                                                     !! the layer to be searched [degC].
  real, dimension(:),        intent(in) :: ppoly_S   !< Coefficients of the polynomial reconstruction of S within
                                                     !! the layer to be searched [ppt].
  real                                  :: z         !< Position where drho = 0 [nondim]
  ! Local variables
  real :: dRdT_diff  ! Difference in the partial derivative of density with temperature across the
                     ! layer [R degC-1 ~> kg m-3 degC-1]
  real :: dRdS_diff  ! Difference in the partial derivative of density with salinity across the
                     ! layer [R ppt-1 ~> kg m-3 ppt-1]
  real :: drho, drho_dz ! Density anomaly and its derivative with fracitonal position [R ~> kg m-3]
  real :: dRdT_z     ! Partial derivative of density with temperature at a point [R degC-1 ~> kg m-3 degC-1]
  real :: dRdS_z     ! Partial derivative of density with salinity at a point [R ppt-1 ~> kg m-3 ppt-1]
  real :: T_z, dT_dz ! Temperature at a point and its derivative with fractional position [degC]
  real :: S_z, dS_dz ! Salinity at a point and its derivative with fractional position [ppt]
  real :: drho_min, drho_max ! Bounds on density differences [R ~> kg m-3]
  real :: ztest, zmin, zmax ! Fractional positions in the cell [nondim]
  real :: dz         ! Change in position in the cell [nondim]
  real :: a1, a2     ! Fractional weights of the top and bottom values [nondim]
  integer :: iter
  integer :: nterm

  nterm = SIZE(ppoly_t)

  ! Position independent quantities
  drdt_diff = drdt_bot - drdt_top
  drds_diff = drds_bot - drds_top
  ! Initial starting drho (used for bisection)
  zmin = z0        ! Lower bounding interval
  zmax = 1.        ! Maximum bounding interval (bottom of layer)
  a1 = 1. - zmin
  a2 = zmin
  t_z = evaluation_polynomial( ppoly_t, nterm, zmin )
  s_z = evaluation_polynomial( ppoly_s, nterm, zmin )
  drdt_z = a1*drdt_top + a2*drdt_bot
  drds_z = a1*drds_top + a2*drds_bot
  drho_min = 0.5*((drdt_z+drdt_ref)*(t_z-t_ref) + (drds_z+drds_ref)*(s_z-s_ref))

  t_z = evaluation_polynomial( ppoly_t, nterm, 1. )
  s_z = evaluation_polynomial( ppoly_s, nterm, 1. )
  drho_max = 0.5*((drdt_bot+drdt_ref)*(t_z-t_ref) + (drds_bot+drds_ref)*(s_z-s_ref))

  if (drho_min >= 0.) then
    z = z0
    return
  elseif (drho_max == 0.) then
    z = 1.
    return
  endif
  if ( sign(1.,drho_min) == sign(1.,drho_max) ) then
    call mom_error(fatal, "drho_min is the same sign as dhro_max")
  endif

  z = z0
  ztest = z0
  do iter = 1, cs%max_iter
    ! Calculate quantities at the current nondimensional position
    a1 = 1.-z
    a2 = z
    drdt_z    = a1*drdt_top + a2*drdt_bot
    drds_z    = a1*drds_top + a2*drds_bot
    t_z       = evaluation_polynomial( ppoly_t, nterm, z )
    s_z       = evaluation_polynomial( ppoly_s, nterm, z )
    drho = 0.5*((drdt_z+drdt_ref)*(t_z-t_ref) + (drds_z+drds_ref)*(s_z-s_ref))

    ! Check for convergence
    if (abs(drho) <= cs%drho_tol) exit
    ! Update bisection bracketing intervals
    if (drho < 0. .and. drho > drho_min) then
      drho_min = drho
      zmin = z
    elseif (drho > 0. .and. drho < drho_max) then
      drho_max = drho
      zmax = z
    endif

    ! Calculate a Newton step
    dt_dz = first_derivative_polynomial( ppoly_t, nterm, z )
    ds_dz = first_derivative_polynomial( ppoly_s, nterm, z )
    drho_dz = 0.5*( (drdt_diff*(t_z - t_ref) + (drdt_ref+drdt_z)*dt_dz) + &
                    (drds_diff*(s_z - s_ref) + (drds_ref+drds_z)*ds_dz) )

    ztest = z - drho/drho_dz
    ! Take a bisection if z falls out of [zmin,zmax]
    if (ztest < zmin .or. ztest > zmax) then
      if ( drho < 0. ) then
        ztest = 0.5*(z + zmax)
      else
        ztest = 0.5*(zmin + z)
      endif
    endif

    ! Test to ensure we haven't stalled out
    if ( abs(z-ztest) <= cs%x_tol ) exit
    ! Reset for next iteration
    z = ztest
  enddo

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]	nk	Number of levels
[in]	pl	Left-column interface pressure [R L2 T-2 ~> Pa] or other units
[in]	tl	Left-column interface potential temperature [degC]
[in]	sl	Left-column interface salinity [ppt]
[in]	drdtl	Left-column dRho/dT [R degC-1 ~> kg m-3 degC-1]
[in]	drdsl	Left-column dRho/dS [R ppt-1 ~> kg m-3 ppt-1]
[in]	pr	Right-column interface pressure [R L2 T-2 ~> Pa] or other units
[in]	tr	Right-column interface potential temperature [degC]
[in]	sr	Right-column interface salinity [ppt]
[in]	drdtr	Left-column dRho/dT [R degC-1 ~> kg m-3 degC-1]
[in]	drdsr	Left-column dRho/dS [R ppt-1 ~> kg m-3 ppt-1]
[in,out]	pol	Fractional position of neutral surface within layer KoL of left column
[in,out]	por	Fractional position of neutral surface within layer KoR of right column
[in,out]	kol	Index of first left interface above neutral surface
[in,out]	kor	Index of first right interface above neutral surface
[in,out]	heff	Effective thickness between two neutral surfaces [R L2 T-2 ~> Pa] or other units following Pl and Pr.
[in]	bl_kl	Layer index of the boundary layer (left)
[in]	bl_kr	Layer index of the boundary layer (right)
[in]	bl_zl	Nondimensional position of the boundary layer (left)
[in]	bl_zr	Nondimensional position of the boundary layer (right)

Definition at line 933 of file MOM_neutral_diffusion.F90.

  integer,                    intent(in)    :: nk    !< Number of levels
  real, dimension(nk+1),      intent(in)    :: Pl    !< Left-column interface pressure [R L2 T-2 ~> Pa] or other units
  real, dimension(nk+1),      intent(in)    :: Tl    !< Left-column interface potential temperature [degC]
  real, dimension(nk+1),      intent(in)    :: Sl    !< Left-column interface salinity [ppt]
  real, dimension(nk+1),      intent(in)    :: dRdTl !< Left-column dRho/dT [R degC-1 ~> kg m-3 degC-1]
  real, dimension(nk+1),      intent(in)    :: dRdSl !< Left-column dRho/dS [R ppt-1 ~> kg m-3 ppt-1]
  real, dimension(nk+1),      intent(in)    :: Pr    !< Right-column interface pressure [R L2 T-2 ~> Pa] or other units
  real, dimension(nk+1),      intent(in)    :: Tr    !< Right-column interface potential temperature [degC]
  real, dimension(nk+1),      intent(in)    :: Sr    !< Right-column interface salinity [ppt]
  real, dimension(nk+1),      intent(in)    :: dRdTr !< Left-column dRho/dT [R degC-1 ~> kg m-3 degC-1]
  real, dimension(nk+1),      intent(in)    :: dRdSr !< Left-column dRho/dS [R ppt-1 ~> kg m-3 ppt-1]
  real, dimension(2*nk+2),    intent(inout) :: PoL   !< Fractional position of neutral surface within
                                                     !! layer KoL of left column
  real, dimension(2*nk+2),    intent(inout) :: PoR   !< Fractional position of neutral surface within
                                                     !! layer KoR of right column
  integer, dimension(2*nk+2), intent(inout) :: KoL   !< Index of first left interface above neutral surface
  integer, dimension(2*nk+2), intent(inout) :: KoR   !< Index of first right interface above neutral surface
  real, dimension(2*nk+1),    intent(inout) :: hEff  !< Effective thickness between two neutral surfaces
                                                     !! [R L2 T-2 ~> Pa] or other units following Pl and Pr.
  integer, optional,          intent(in)    :: bl_kl !< Layer index of the boundary layer (left)
  integer, optional,          intent(in)    :: bl_kr !< Layer index of the boundary layer (right)
  real, optional,             intent(in)    :: bl_zl !< Nondimensional position of the boundary layer (left)
  real, optional,             intent(in)    :: bl_zr !< Nondimensional position of the boundary layer (right)

  ! Local variables
  integer :: ns                     ! Number of neutral surfaces
  integer :: k_surface              ! Index of neutral surface
  integer :: kl                     ! Index of left interface
  integer :: kr                     ! Index of right interface
  real    :: dRdT, dRdS             ! dRho/dT [kg m-3 degC-1] and dRho/dS [kg m-3 ppt-1] for the neutral surface
  logical :: searching_left_column  ! True if searching for the position of a right interface in the left column
  logical :: searching_right_column ! True if searching for the position of a left interface in the right column
  logical :: reached_bottom         ! True if one of the bottom-most interfaces has been used as the target
  integer :: krm1, klm1
  real    :: dRho, dRhoTop, dRhoBot ! Potential density differences at various points [R ~> kg m-3]
  real    :: hL, hR                 ! Pressure thicknesses [R L2 T-2 ~> Pa]
  integer :: lastK_left, lastK_right ! Layers used during the last iteration
  real    :: lastP_left, lastP_right ! Fractional positions during the last iteration [nondim]
  logical :: interior_limit

  ns = 2*nk+2

  ! Initialize variables for the search
  kr = 1 ;
  kl = 1 ;
  lastp_right = 0.
  lastp_left = 0.
  lastk_right = 1
  lastk_left  = 1
  reached_bottom = .false.

  ! Check to see if we should limit the diffusion to the interior
  interior_limit = PRESENT(bl_kl) .and. PRESENT(bl_kr) .and. PRESENT(bl_zr) .and. PRESENT(bl_zl)

  ! Loop over each neutral surface, working from top to bottom
  neutral_surfaces: do k_surface = 1, ns
    klm1 = max(kl-1, 1)
    if (klm1>nk) stop 'find_neutral_surface_positions(): klm1 went out of bounds!'
    krm1 = max(kr-1, 1)
    if (krm1>nk) stop 'find_neutral_surface_positions(): krm1 went out of bounds!'

    ! Potential density difference, rho(kr) - rho(kl)
    drho = 0.5 * ( ( drdtr(kr) + drdtl(kl) ) * ( tr(kr) - tl(kl) ) &
                 + ( drdsr(kr) + drdsl(kl) ) * ( sr(kr) - sl(kl) ) )
    ! Which column has the lighter surface for the current indexes, kr and kl
    if (.not. reached_bottom) then
      if (drho < 0.) then
        searching_left_column = .true.
        searching_right_column = .false.
      elseif (drho > 0.) then
        searching_right_column = .true.
        searching_left_column = .false.
      else ! dRho == 0.
        if (kl + kr == 2) then ! Still at surface
          searching_left_column = .true.
          searching_right_column = .false.
        else ! Not the surface so we simply change direction
          searching_left_column = .not.  searching_left_column
          searching_right_column = .not.  searching_right_column
        endif
      endif
    endif

    if (searching_left_column) then
      ! Interpolate for the neutral surface position within the left column, layer klm1
      ! Potential density difference, rho(kl-1) - rho(kr) (should be negative)
      drhotop = 0.5 * ( ( drdtl(klm1) + drdtr(kr) ) * ( tl(klm1) - tr(kr) ) &
                     + ( drdsl(klm1) + drdsr(kr) ) * ( sl(klm1) - sr(kr) ) )
      ! Potential density difference, rho(kl) - rho(kr) (will be positive)
      drhobot = 0.5 * ( ( drdtl(klm1+1) + drdtr(kr) ) * ( tl(klm1+1) - tr(kr) ) &
                      + ( drdsl(klm1+1) + drdsr(kr) ) * ( sl(klm1+1) - sr(kr) ) )

      ! Because we are looking left, the right surface, kr, is lighter than klm1+1 and should be denser than klm1
      ! unless we are still at the top of the left column (kl=1)
      if (drhotop > 0. .or. kr+kl==2) then
        pol(k_surface) = 0. ! The right surface is lighter than anything in layer klm1
      elseif (drhotop >= drhobot) then ! Left layer is unstratified
        pol(k_surface) = 1.
      else
        ! Linearly interpolate for the position between Pl(kl-1) and Pl(kl) where the density difference
        ! between right and left is zero. The Pl here are only used to handle massless layers.
        pol(k_surface) = interpolate_for_nondim_position( drhotop, pl(klm1), drhobot, pl(klm1+1) )
      endif
      if (pol(k_surface)>=1. .and. klm1<nk) then ! >= is really ==, when PoL==1 we point to the bottom of the cell
        klm1 = klm1 + 1
        pol(k_surface) = pol(k_surface) - 1.
      endif
      if (real(klm1-lastk_left)+(pol(k_surface)-lastp_left)<0.) then
        pol(k_surface) = lastp_left
        klm1 = lastk_left
      endif
      kol(k_surface) = klm1
      if (kr <= nk) then
        por(k_surface) = 0.
        kor(k_surface) = kr
      else
        por(k_surface) = 1.
        kor(k_surface) = nk
      endif
      if (kr <= nk) then
        kr = kr + 1
      else
        reached_bottom = .true.
        searching_right_column = .true.
        searching_left_column = .false.
      endif
    elseif (searching_right_column) then
      ! Interpolate for the neutral surface position within the right column, layer krm1
      ! Potential density difference, rho(kr-1) - rho(kl) (should be negative)
      drhotop = 0.5 * ( ( drdtr(krm1) + drdtl(kl) ) * ( tr(krm1) - tl(kl) ) + &
                        ( drdsr(krm1) + drdsl(kl) ) * ( sr(krm1) - sl(kl) ) )
      ! Potential density difference, rho(kr) - rho(kl) (will be positive)
      drhobot = 0.5 * ( ( drdtr(krm1+1) + drdtl(kl) ) * ( tr(krm1+1) - tl(kl) ) + &
                        ( drdsr(krm1+1) + drdsl(kl) ) * ( sr(krm1+1) - sl(kl) ) )

      ! Because we are looking right, the left surface, kl, is lighter than krm1+1 and should be denser than krm1
      ! unless we are still at the top of the right column (kr=1)
      if (drhotop >= 0. .or. kr+kl==2) then
        por(k_surface) = 0. ! The left surface is lighter than anything in layer krm1
      elseif (drhotop >= drhobot) then ! Right layer is unstratified
        por(k_surface) = 1.
      else
        ! Linearly interpolate for the position between Pr(kr-1) and Pr(kr) where the density difference
        ! between right and left is zero. The Pr here are only used to handle massless layers.
        por(k_surface) = interpolate_for_nondim_position( drhotop, pr(krm1), drhobot, pr(krm1+1) )
      endif
      if (por(k_surface)>=1. .and. krm1<nk) then ! >= is really ==, when PoR==1 we point to the bottom of the cell
        krm1 = krm1 + 1
        por(k_surface) = por(k_surface) - 1.
      endif
      if (real(krm1-lastk_right)+(por(k_surface)-lastp_right)<0.) then
        por(k_surface) = lastp_right
        krm1 = lastk_right
      endif
      kor(k_surface) = krm1
      if (kl <= nk) then
        pol(k_surface) = 0.
        kol(k_surface) = kl
      else
        pol(k_surface) = 1.
        kol(k_surface) = nk
      endif
      if (kl <= nk) then
        kl = kl + 1
      else
        reached_bottom = .true.
        searching_right_column = .false.
        searching_left_column = .true.
      endif
    else
      stop 'Else what?'
    endif
    if (interior_limit) then
      if (kol(k_surface)<=bl_kl) then
        kol(k_surface) = bl_kl
        if (pol(k_surface)<bl_zl) then
          pol(k_surface) = bl_zl
        endif
      endif
      if (kor(k_surface)<=bl_kr) then
        kor(k_surface) = bl_kr
        if (por(k_surface)<bl_zr) then
          por(k_surface) = bl_zr
        endif
      endif
    endif

    lastk_left = kol(k_surface) ; lastp_left = pol(k_surface)
    lastk_right = kor(k_surface) ; lastp_right = por(k_surface)
    ! Effective thickness
    ! NOTE: This would be better expressed in terms of the layers thicknesses rather
    ! than as differences of position - AJA
    if (k_surface>1) then
      hl = absolute_position(nk,ns,pl,kol,pol,k_surface) - absolute_position(nk,ns,pl,kol,pol,k_surface-1)
      hr = absolute_position(nk,ns,pr,kor,por,k_surface) - absolute_position(nk,ns,pr,kor,por,k_surface-1)
      if ( hl + hr > 0.) then
        heff(k_surface-1) = 2. * hl * hr / ( hl + hr ) ! Harmonic mean of layer thicknesses
      else
        heff(k_surface-1) = 0.
      endif
    endif

  enddo neutral_surfaces

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]	cs	Neutral diffusion control structure
[in]	nk	Number of levels
[in]	pres_l	Left-column interface pressure [R L2 T-2 ~> Pa]
[in]	hcol_l	Left-column layer thicknesses [H ~> m or kg m-2] or other units
[in]	tl	Left-column top interface potential temperature [degC]
[in]	sl	Left-column top interface salinity [ppt]
[in]	ppoly_t_l	Left-column coefficients of T reconstruction [degC]
[in]	ppoly_s_l	Left-column coefficients of S reconstruction [ppt]
[in]	stable_l	True where the left-column is stable
[in]	pres_r	Right-column interface pressure [R L2 T-2 ~> Pa]
[in]	hcol_r	Left-column layer thicknesses [H ~> m or kg m-2] or other units
[in]	tr	Right-column top interface potential temperature [degC]
[in]	sr	Right-column top interface salinity [ppt]
[in]	ppoly_t_r	Right-column coefficients of T reconstruction [degC]
[in]	ppoly_s_r	Right-column coefficients of S reconstruction [ppt]
[in]	stable_r	True where the right-column is stable
[in,out]	pol	Fractional position of neutral surface within layer KoL of left column [nondim]
[in,out]	por	Fractional position of neutral surface within layer KoR of right column [nondim]
[in,out]	kol	Index of first left interface above neutral surface
[in,out]	kor	Index of first right interface above neutral surface
[in,out]	heff	Effective thickness between two neutral surfaces [H ~> m or kg m-2] or other units taken from hcol_l
[in]	zeta_bot_l	Non-dimensional distance to where the boundary layer intersetcs the cell (left) [nondim]
[in]	zeta_bot_r	Non-dimensional distance to where the boundary layer intersetcs the cell (right) [nondim]
[in]	k_bot_l	k-index for the boundary layer (left) [nondim]
[in]	k_bot_r	k-index for the boundary layer (right) [nondim]
[in]	hard_fail_heff	If true (default) bring down the model if the neutral surfaces ever cross [logical]

Definition at line 1187 of file MOM_neutral_diffusion.F90.

  type(neutral_diffusion_CS),     intent(inout) :: CS        !< Neutral diffusion control structure
  integer,                        intent(in)    :: nk        !< Number of levels
  real, dimension(nk,2),          intent(in)    :: Pres_l    !< Left-column interface pressure [R L2 T-2 ~> Pa]
  real, dimension(nk),            intent(in)    :: hcol_l    !< Left-column layer thicknesses [H ~> m or kg m-2]
                                                             !! or other units
  real, dimension(nk,2),          intent(in)    :: Tl        !< Left-column top interface potential temperature [degC]
  real, dimension(nk,2),          intent(in)    :: Sl        !< Left-column top interface salinity [ppt]
  real, dimension(:,:),           intent(in)    :: ppoly_T_l !< Left-column coefficients of T reconstruction [degC]
  real, dimension(:,:),           intent(in)    :: ppoly_S_l !< Left-column coefficients of S reconstruction [ppt]
  logical, dimension(nk),         intent(in)    :: stable_l  !< True where the left-column is stable
  real, dimension(nk,2),          intent(in)    :: Pres_r    !< Right-column interface pressure [R L2 T-2 ~> Pa]
  real, dimension(nk),            intent(in)    :: hcol_r    !< Left-column layer thicknesses [H ~> m or kg m-2]
                                                             !! or other units
  real, dimension(nk,2),          intent(in)    :: Tr        !< Right-column top interface potential temperature [degC]
  real, dimension(nk,2),          intent(in)    :: Sr        !< Right-column top interface salinity [ppt]
  real, dimension(:,:),           intent(in)    :: ppoly_T_r !< Right-column coefficients of T reconstruction [degC]
  real, dimension(:,:),           intent(in)    :: ppoly_S_r !< Right-column coefficients of S reconstruction [ppt]
  logical, dimension(nk),         intent(in)    :: stable_r  !< True where the right-column is stable
  real, dimension(4*nk),          intent(inout) :: PoL       !< Fractional position of neutral surface within
                                                             !! layer KoL of left column [nondim]
  real, dimension(4*nk),          intent(inout) :: PoR       !< Fractional position of neutral surface within
                                                             !! layer KoR of right column [nondim]
  integer, dimension(4*nk),       intent(inout) :: KoL       !< Index of first left interface above neutral surface
  integer, dimension(4*nk),       intent(inout) :: KoR       !< Index of first right interface above neutral surface
  real, dimension(4*nk-1),        intent(inout) :: hEff      !< Effective thickness between two neutral surfaces
                                                             !! [H ~> m or kg m-2] or other units taken from hcol_l
  real, optional,                 intent(in)    :: zeta_bot_L!< Non-dimensional distance to where the boundary layer
                                                             !! intersetcs the cell (left) [nondim]
  real, optional,                 intent(in)    :: zeta_bot_R!< Non-dimensional distance to where the boundary layer
                                                             !! intersetcs the cell (right) [nondim]

  integer, optional,              intent(in)    :: k_bot_L   !< k-index for the boundary layer (left) [nondim]
  integer, optional,              intent(in)    :: k_bot_R   !< k-index for the boundary layer (right) [nondim]
  logical, optional,              intent(in)    :: hard_fail_heff !< If true (default) bring down the model if the
                                                             !! neutral surfaces ever cross [logical]
  ! Local variables
  integer :: ns                     ! Number of neutral surfaces
  integer :: k_surface              ! Index of neutral surface
  integer :: kl_left, kl_right      ! Index of layers on the left/right
  integer :: ki_left, ki_right      ! Index of interfaces on the left/right
  logical :: searching_left_column  ! True if searching for the position of a right interface in the left column
  logical :: searching_right_column ! True if searching for the position of a left interface in the right column
  logical :: reached_bottom         ! True if one of the bottom-most interfaces has been used as the target
  logical :: search_layer
  logical :: fail_heff              ! Fail if negative thickness are encountered.  By default this
                                    ! is true, but it can take its value from hard_fail_heff.
  real    :: dRho                   ! A density difference between columns [R ~> kg m-3]
  real    :: hL, hR                 ! Left and right layer thicknesses [H ~> m or kg m-2] or units from hcol_l
  real    :: lastP_left, lastP_right ! Previous positions for left and right [nondim]
  integer :: k_init_L, k_init_R      ! Starting indices layers for left and right
  real    :: p_init_L, p_init_R      ! Starting positions for left and right [nondim]
  ! Initialize variables for the search
  ns = 4*nk
  ki_right = 1
  ki_left = 1
  kl_left = 1
  kl_right = 1
  lastp_left = 0.
  lastp_right = 0.
  reached_bottom = .false.
  searching_left_column = .false.
  searching_right_column = .false.

  fail_heff = .true.
  if (PRESENT(hard_fail_heff)) fail_heff = hard_fail_heff

  if (PRESENT(k_bot_l) .and. PRESENT(k_bot_r) .and. PRESENT(zeta_bot_l) .and. PRESENT(zeta_bot_r)) then
    k_init_l = k_bot_l; k_init_r = k_bot_r
    p_init_l = zeta_bot_l; p_init_r = zeta_bot_r
    lastp_left = zeta_bot_l; lastp_right = zeta_bot_r
    kl_left = k_bot_l; kl_right = k_bot_r
  else
    k_init_l = 1  ; k_init_r = 1
    p_init_l = 0. ; p_init_r = 0.
  endif
  ! Loop over each neutral surface, working from top to bottom
  neutral_surfaces: do k_surface = 1, ns

    if (k_surface == ns) then
        pol(k_surface) = 1.
        por(k_surface) = 1.
        kol(k_surface) = nk
        kor(k_surface) = nk
    ! If the layers are unstable, then simply point the surface to the previous location
    elseif (.not. stable_l(kl_left)) then
      if (k_surface > 1) then
        pol(k_surface) = ki_left - 1 ! Top interface is at position = 0., Bottom is at position = 1
        kol(k_surface) = kl_left
        por(k_surface) = por(k_surface-1)
        kor(k_surface) = kor(k_surface-1)
      else
        por(k_surface) = p_init_r
        kor(k_surface) = k_init_r
        pol(k_surface) = p_init_l
        kol(k_surface) = k_init_l
      endif
      call increment_interface(nk, kl_left, ki_left, reached_bottom, searching_left_column, searching_right_column)
      searching_left_column = .true.
      searching_right_column = .false.
    elseif (.not. stable_r(kl_right)) then ! Check the right layer for stability
      if (k_surface > 1) then
        por(k_surface) = ki_right - 1 ! Top interface is at position = 0., Bottom is at position = 1
        kor(k_surface) = kl_right
        pol(k_surface) = pol(k_surface-1)
        kol(k_surface) = kol(k_surface-1)
      else
        por(k_surface) = 0.
        kor(k_surface) = 1
        pol(k_surface) = 0.
        kol(k_surface) = 1
      endif
      call increment_interface(nk, kl_right, ki_right, reached_bottom, searching_right_column, searching_left_column)
      searching_left_column = .false.
      searching_right_column = .true.
    else ! Layers are stable so need to figure out whether we need to search right or left
      ! For convenience, the left column uses the searched "from" interface variables, and the right column
      ! uses the searched 'to'. These will get reset in subsequent calc_delta_rho calls

      call calc_delta_rho_and_derivs(cs,                                                                        &
                                     tr(kl_right, ki_right), sr(kl_right, ki_right), pres_r(kl_right,ki_right), &
                                     tl(kl_left, ki_left),   sl(kl_left, ki_left)  , pres_l(kl_left,ki_left),   &
                                     drho)
      if (cs%debug) write(stdout,'(A,I2,A,E12.4,A,I2,A,I2,A,I2,A,I2)') &
          "k_surface=",k_surface, "  dRho=",cs%R_to_kg_m3*drho, &
          "kl_left=",kl_left, "  ki_left=",ki_left, "  kl_right=",kl_right, "  ki_right=",ki_right
      ! Which column has the lighter surface for the current indexes, kr and kl
      if (.not. reached_bottom) then
        if (drho < 0.) then
          searching_left_column  = .true.
          searching_right_column = .false.
        elseif (drho > 0.) then
          searching_left_column  = .false.
          searching_right_column = .true.
        else ! dRho == 0.
          if (  ( kl_left + kl_right == 2 ) .and. (ki_left + ki_right == 2) ) then ! Still at surface
            searching_left_column  = .true.
            searching_right_column = .false.
          else ! Not the surface so we simply change direction
            searching_left_column = .not. searching_left_column
            searching_right_column = .not. searching_right_column
          endif
        endif
      endif
      if (searching_left_column) then
        ! Position of the right interface is known and all quantities are fixed
        por(k_surface) = ki_right - 1.
        kor(k_surface) = kl_right
        pol(k_surface) = search_other_column(cs, k_surface, lastp_left,                                &
                           tr(kl_right, ki_right), sr(kl_right, ki_right), pres_r(kl_right, ki_right), &
                           tl(kl_left,1),          sl(kl_left,1),          pres_l(kl_left,1),          &
                           tl(kl_left,2),          sl(kl_left,2),          pres_l(kl_left,2),          &
                           ppoly_t_l(kl_left,:), ppoly_s_l(kl_left,:))
        kol(k_surface) = kl_left

        if (cs%debug) then
          write(stdout,'(A,I2)') "Searching left layer ", kl_left
          write(stdout,'(A,I2,X,I2)') "Searching from right: ", kl_right, ki_right
          write(stdout,*) "Temp/Salt Reference: ", tr(kl_right,ki_right), sr(kl_right,ki_right)
          write(stdout,*) "Temp/Salt Top L: ", tl(kl_left,1), sl(kl_left,1)
          write(stdout,*) "Temp/Salt Bot L: ", tl(kl_left,2), sl(kl_left,2)
        endif
        call increment_interface(nk, kl_right, ki_right, reached_bottom, searching_right_column, searching_left_column)
        lastp_left = pol(k_surface)
        ! If the right layer increments, then we need to reset the last position on the right
        if ( kl_right == (kor(k_surface) + 1) ) lastp_right = 0.

      elseif (searching_right_column) then
        ! Position of the right interface is known and all quantities are fixed
        pol(k_surface) = ki_left - 1.
        kol(k_surface) = kl_left
        por(k_surface) = search_other_column(cs, k_surface, lastp_right,                         &
                           tl(kl_left, ki_left), sl(kl_left, ki_left), pres_l(kl_left, ki_left), &
                           tr(kl_right,1),       sr(kl_right,1),       pres_r(kl_right,1),       &
                           tr(kl_right,2),       sr(kl_right,2),       pres_r(kl_right,2),       &
                           ppoly_t_r(kl_right,:), ppoly_s_r(kl_right,:))
        kor(k_surface) = kl_right

        if (cs%debug) then
          write(stdout,'(A,I2)') "Searching right layer ", kl_right
          write(stdout,'(A,I2,X,I2)') "Searching from left: ", kl_left, ki_left
          write(stdout,*) "Temp/Salt Reference: ", tl(kl_left,ki_left), sl(kl_left,ki_left)
          write(stdout,*) "Temp/Salt Top L: ", tr(kl_right,1), sr(kl_right,1)
          write(stdout,*) "Temp/Salt Bot L: ", tr(kl_right,2), sr(kl_right,2)
        endif
        call increment_interface(nk, kl_left, ki_left, reached_bottom, searching_left_column, searching_right_column)
        lastp_right = por(k_surface)
        ! If the right layer increments, then we need to reset the last position on the right
        if ( kl_left == (kol(k_surface) + 1) ) lastp_left = 0.
      else
        stop 'Else what?'
      endif
      if (cs%debug)  write(stdout,'(A,I3,A,ES16.6,A,I2,A,ES16.6)') "KoL:", kol(k_surface), " PoL:", pol(k_surface), &
                     "     KoR:", kor(k_surface), " PoR:", por(k_surface)
    endif
    ! Effective thickness
    if (k_surface>1) then
      if ( kol(k_surface) == kol(k_surface-1) .and. kor(k_surface) == kor(k_surface-1) ) then
        hl = (pol(k_surface) - pol(k_surface-1))*hcol_l(kol(k_surface))
        hr = (por(k_surface) - por(k_surface-1))*hcol_r(kor(k_surface))
        if (hl < 0. .or. hr < 0.) then
          if (fail_heff) then
            call mom_error(fatal,"Negative thicknesses in neutral diffusion")
          else
            if (searching_left_column) then
              pol(k_surface) = pol(k_surface-1)
              kol(k_surface) = kol(k_surface-1)
            elseif (searching_right_column) then
              por(k_surface) = por(k_surface-1)
              kor(k_surface) = kor(k_surface-1)
            endif
          endif
        elseif ( hl + hr == 0. ) then
          heff(k_surface-1) = 0.
        else
          heff(k_surface-1) = 2. * ( (hl * hr) / ( hl + hr ) )! Harmonic mean
          if ( kol(k_surface) /= kol(k_surface-1) ) then
            call mom_error(fatal,"Neutral sublayer spans multiple layers")
          endif
          if ( kor(k_surface) /= kor(k_surface-1) ) then
            call mom_error(fatal,"Neutral sublayer spans multiple layers")
          endif
        endif
      else
        heff(k_surface-1) = 0.
      endif
    endif
  enddo neutral_surfaces

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]	hkm1	Left cell width
[in]	hk	Center cell width
[in]	hkp1	Right cell width
[in]	skm1	Left cell average value
[in]	sk	Center cell average value
[in]	skp1	Right cell average value

Definition at line 876 of file MOM_neutral_diffusion.F90.

  real, intent(in) :: hkm1 !< Left cell width
  real, intent(in) :: hk   !< Center cell width
  real, intent(in) :: hkp1 !< Right cell width
  real, intent(in) :: Skm1 !< Left cell average value
  real, intent(in) :: Sk   !< Center cell average value
  real, intent(in) :: Skp1 !< Right cell average value

  ! Local variables
  real :: h_sum, hp, hm

  h_sum = ( hkm1 + hkp1 ) + hk
  if (h_sum /= 0.) h_sum = 1./ h_sum
  hm =  hkm1 + hk
  if (hm /= 0.) hm = 1./ hm
  hp =  hkp1 + hk
  if (hp /= 0.) hp = 1./ hp
  fv_diff = ( hk * h_sum ) * &
            (   ( 2. * hkm1 + hk ) * hp * ( skp1 - sk ) &
              + ( 2. * hkp1 + hk ) * hm * ( sk - skm1 ) )

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]	hkm1	Left cell width
[in]	hk	Center cell width
[in]	hkp1	Right cell width
[in]	skm1	Left cell average value
[in]	sk	Center cell average value
[in]	skp1	Right cell average value

Definition at line 902 of file MOM_neutral_diffusion.F90.

  real, intent(in) :: hkm1 !< Left cell width
  real, intent(in) :: hk   !< Center cell width
  real, intent(in) :: hkp1 !< Right cell width
  real, intent(in) :: Skm1 !< Left cell average value
  real, intent(in) :: Sk   !< Center cell average value
  real, intent(in) :: Skp1 !< Right cell average value

  ! Local variables
  real :: xkm1, xkp1
  real :: h_sum, hx_sum, hxsq_sum, hxy_sum, hy_sum, det

  xkm1 = -0.5 * ( hk + hkm1 )
  xkp1 = 0.5 * ( hk + hkp1 )
  h_sum = ( hkm1 + hkp1 ) + hk
  hx_sum = hkm1*xkm1 + hkp1*xkp1
  hxsq_sum = hkm1*(xkm1**2) + hkp1*(xkp1**2)
  hxy_sum = hkm1*xkm1*skm1 + hkp1*xkp1*skp1
  hy_sum = ( hkm1*skm1 + hkp1*skp1 ) + hk*sk
  det = h_sum * hxsq_sum - hx_sum**2
  if (det /= 0.) then
    !a = ( hxsq_sum * hy_sum - hx_sum*hxy_sum ) / det ! a would be mean of straight line fit
    fvlsq_slope = ( h_sum * hxy_sum - hx_sum*hy_sum ) / det ! Gradient of straight line fit
  else
    fvlsq_slope = 0. ! Adcroft's reciprocal rule
  endif

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]	nk	Number of vertical levels
[in,out]	kl	Current layer (potentially updated)
[in,out]	ki	Current interface
[in,out]	reached_bottom	Updated when kl == nk and ki == 2
[in,out]	searching_this_column	Updated when kl == nk and ki == 2
[in,out]	searching_other_column	Updated when kl == nk and ki == 2

Definition at line 1507 of file MOM_neutral_diffusion.F90.

  integer, intent(in   )                :: nk                     !< Number of vertical levels
  integer, intent(inout)                :: kl                     !< Current layer (potentially updated)
  integer, intent(inout)                :: ki                     !< Current interface
  logical, intent(inout)                :: reached_bottom         !< Updated when kl == nk and ki == 2
  logical, intent(inout)                :: searching_this_column  !< Updated when kl == nk and ki == 2
  logical, intent(inout)                :: searching_other_column !< Updated when kl == nk and ki == 2
  integer :: k

  reached_bottom = .false.
  if (ki == 2) then ! At the bottom interface
    if ((ki == 2) .and. (kl < nk) ) then ! Not at the bottom so just go to the next layer
      kl = kl+1
      ki = 1
    elseif ((kl == nk) .and. (ki==2)) then
      reached_bottom = .true.
      searching_this_column = .false.
      searching_other_column = .true.
    endif
  elseif (ki==1) then ! At the top interface
    ki = 2 ! Next interface is same layer, but bottom interface
  else
    call mom_error(fatal,"Unanticipated eventuality in increment_interface")
  endif

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]	nk	Number of levels
[in]	h	Layer thickness [H ~> m or kg m-2]
[in]	s	Layer scalar (conc, e.g. ppt)
[in,out]	si	Interface scalar (conc, e.g. ppt)
[in]	i_method	=1 use average of PLM edges =2 use continuous PPM edge interpolation
[in]	h_neglect	A negligibly small thickness [H ~> m or kg m-2]

Definition at line 693 of file MOM_neutral_diffusion.F90.

  integer,               intent(in)    :: nk       !< Number of levels
  real, dimension(nk),   intent(in)    :: h        !< Layer thickness [H ~> m or kg m-2]
  real, dimension(nk),   intent(in)    :: S        !< Layer scalar (conc, e.g. ppt)
  real, dimension(nk+1), intent(inout) :: Si       !< Interface scalar (conc, e.g. ppt)
  integer,               intent(in)    :: i_method !< =1 use average of PLM edges
                                                   !! =2 use continuous PPM edge interpolation
  real,                  intent(in)    :: h_neglect !< A negligibly small thickness [H ~> m or kg m-2]
  ! Local variables
  integer :: k, km2, kp1
  real, dimension(nk) :: diff
  real :: Sb, Sa

  call plm_diff(nk, h, s, 2, 1, diff)
  si(1) = s(1) - 0.5 * diff(1)
  if (i_method==1) then
    do k = 2, nk
      ! Average of the two edge values (will be bounded and,
      ! when slopes are unlimited, notionally second-order accurate)
      sa = s(k-1) + 0.5 * diff(k-1) ! Lower edge value of a PLM reconstruction for layer above
      sb = s(k) - 0.5 * diff(k) ! Upper edge value of a PLM reconstruction for layer below
      si(k) = 0.5 * ( sa + sb )
    enddo
  elseif (i_method==2) then
    do k = 2, nk
      ! PPM quasi-fourth order interpolation for edge values following
      ! equation 1.6 in Colella & Woodward, 1984: JCP 54, 174-201.
      km2 = max(1, k-2)
      kp1 = min(nk, k+1)
      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)
    enddo
  endif
  si(nk+1) = s(nk) + 0.5 * diff(nk)

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]	drhoneg	Negative density difference [R ~> kg m-3]
[in]	pneg	Position of negative density difference [R L2 T-2 ~> Pa] or [nondim]
[in]	drhopos	Positive density difference [R ~> kg m-3]
[in]	ppos	Position of positive density difference [R L2 T-2 ~> Pa] or [nondim]

Definition at line 1143 of file MOM_neutral_diffusion.F90.

  real, intent(in) :: dRhoNeg !< Negative density difference [R ~> kg m-3]
  real, intent(in) :: Pneg    !< Position of negative density difference [R L2 T-2 ~> Pa] or [nondim]
  real, intent(in) :: dRhoPos !< Positive density difference [R ~> kg m-3]
  real, intent(in) :: Ppos    !< Position of positive density difference [R L2 T-2 ~> Pa] or [nondim]

  character(len=120) :: mesg

  if (ppos < pneg) then
    call mom_error(fatal, 'interpolate_for_nondim_position: Houston, we have a problem! Ppos<Pneg')
  elseif (drhoneg>drhopos) then
    write(stderr,*) 'dRhoNeg, Pneg, dRhoPos, Ppos=',drhoneg, pneg, drhopos, ppos
    write(mesg,*) 'dRhoNeg, Pneg, dRhoPos, Ppos=', drhoneg, pneg, drhopos, ppos
    call mom_error(warning, 'interpolate_for_nondim_position: '//trim(mesg))
  elseif (drhoneg>drhopos) then !### Does this duplicated test belong here?
    call mom_error(fatal, 'interpolate_for_nondim_position: Houston, we have a problem! dRhoNeg>dRhoPos')
  endif
  if (ppos<=pneg) then ! Handle vanished or inverted layers
    interpolate_for_nondim_position = 0.5
  elseif ( drhopos - drhoneg > 0. ) then
    interpolate_for_nondim_position = min( 1., max( 0., -drhoneg / ( drhopos - drhoneg ) ) )
  elseif ( drhopos - drhoneg == 0) then
    if (drhoneg>0.) then
      interpolate_for_nondim_position = 0.
    elseif (drhoneg<0.) then
      interpolate_for_nondim_position = 1.
    else ! dRhoPos = dRhoNeg = 0
      interpolate_for_nondim_position = 0.5
    endif
  else ! dRhoPos - dRhoNeg < 0
    interpolate_for_nondim_position = 0.5
  endif
  if ( interpolate_for_nondim_position < 0. ) &
    call mom_error(fatal, 'interpolate_for_nondim_position: Houston, we have a problem! Pint < Pneg')
  if ( interpolate_for_nondim_position > 1. ) &
    call mom_error(fatal, 'interpolate_for_nondim_position: Houston, we have a problem! Pint > Ppos')

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]	cs	Neutral diffusion control structure
[in]	nk	Number of levels in a column
[in]	t	Temperature at interfaces [degC]
[in]	s	Salinity at interfaces [ppt]
[in]	p	Pressure at interfaces [R L2 T-2 ~> Pa]
[out]	stable_cell	True if this cell is unstably stratified

Definition at line 1419 of file MOM_neutral_diffusion.F90.

  type(neutral_diffusion_CS), intent(inout) :: CS      !< Neutral diffusion control structure
  integer,                intent(in)    :: nk          !< Number of levels in a column
  real, dimension(nk,2),  intent(in)    :: T           !< Temperature at interfaces [degC]
  real, dimension(nk,2),  intent(in)    :: S           !< Salinity at interfaces [ppt]
  real, dimension(nk,2),  intent(in)    :: P           !< Pressure at interfaces [R L2 T-2 ~> Pa]
  logical, dimension(nk), intent(  out) :: stable_cell !< True if this cell is unstably stratified

  integer :: k, first_stable, prev_stable
  real :: delta_rho ! A density difference [R ~> kg m-3]

  do k = 1,nk
    call calc_delta_rho_and_derivs( cs, t(k,2), s(k,2), max(p(k,2), cs%ref_pres), &
                                        t(k,1), s(k,1), max(p(k,1), cs%ref_pres), delta_rho )
    stable_cell(k) = (delta_rho > 0.)
  enddo

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]

verbose

If true, write results to stdout

Definition at line 2089 of file MOM_neutral_diffusion.F90.

  logical, intent(in) :: verbose !< If true, write results to stdout
  ! Local variables
  integer, parameter         :: nk = 4
  real, dimension(nk+1)      :: TiL, TiR1, TiR2, TiR4, Tio ! Test interface temperatures
  real, dimension(nk)        :: TL                         ! Test layer temperatures
  real, dimension(nk+1)      :: SiL                        ! Test interface salinities
  real, dimension(nk+1)      :: PiL, PiR4                  ! Test interface positions
  real, dimension(2*nk+2)    :: PiLRo, PiRLo               ! Test positions
  integer, dimension(2*nk+2) :: KoL, KoR                   ! Test indexes
  real, dimension(2*nk+1)    :: hEff                       ! Test positions
  real, dimension(2*nk+1)    :: Flx                        ! Test flux
  integer :: k
  logical :: v
  real :: h_neglect

  h_neglect = 1.0e-30

  v = verbose

  ndiff_unit_tests_continuous = .false. ! Normally return false
  write(stdout,*) '==== MOM_neutral_diffusion: ndiff_unit_tests_continuous ='

  ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. &
    test_fv_diff(v,1.,1.,1., 0.,1.,2., 1., 'FV: Straight line on uniform grid')
  ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. &
    test_fv_diff(v,1.,1.,0., 0.,4.,8., 7., 'FV: Vanished right cell')
  ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. &
    test_fv_diff(v,0.,1.,1., 0.,4.,8., 7., 'FV: Vanished left cell')
  ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. &
    test_fv_diff(v,1.,2.,4., 0.,3.,9., 4., 'FV: Stretched grid')
  ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. &
    test_fv_diff(v,2.,0.,2., 0.,1.,2., 0., 'FV: Vanished middle cell')
  ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. &
    test_fv_diff(v,0.,1.,0., 0.,1.,2., 2., 'FV: Vanished on both sides')
  ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. &
    test_fv_diff(v,1.,0.,0., 0.,1.,2., 0., 'FV: Two vanished cell sides')
  ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. &
    test_fv_diff(v,0.,0.,0., 0.,1.,2., 0., 'FV: All vanished cells')

  ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. &
    test_fvlsq_slope(v,1.,1.,1., 0.,1.,2., 1., 'LSQ: Straight line on uniform grid')
  ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. &
    test_fvlsq_slope(v,1.,1.,0., 0.,1.,2., 1., 'LSQ: Vanished right cell')
  ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. &
    test_fvlsq_slope(v,0.,1.,1., 0.,1.,2., 1., 'LSQ: Vanished left cell')
  ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. &
    test_fvlsq_slope(v,1.,2.,4., 0.,3.,9., 2., 'LSQ: Stretched grid')
  ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. &
    test_fvlsq_slope(v,1.,0.,1., 0.,1.,2., 2., 'LSQ: Vanished middle cell')
  ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. &
    test_fvlsq_slope(v,0.,1.,0., 0.,1.,2., 0., 'LSQ: Vanished on both sides')
  ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. &
    test_fvlsq_slope(v,1.,0.,0., 0.,1.,2., 0., 'LSQ: Two vanished cell sides')
  ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. &
    test_fvlsq_slope(v,0.,0.,0., 0.,1.,2., 0., 'LSQ: All vanished cells')

  call interface_scalar(4, (/10.,10.,10.,10./), (/24.,18.,12.,6./), tio, 1, h_neglect)
  !ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. &
  !  test_data1d(5, Tio, (/27.,21.,15.,9.,3./), 'Linear profile, interface temperatures')
  ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. &
    test_data1d(v,5, tio, (/24.,22.5,15.,7.5,6./), 'Linear profile, linear interface temperatures')
  call interface_scalar(4, (/10.,10.,10.,10./), (/24.,18.,12.,6./), tio, 2, h_neglect)
  ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. &
    test_data1d(v,5, tio, (/24.,22.,15.,8.,6./), 'Linear profile, PPM interface temperatures')

  ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. &
    test_ifndp(v,-1.0, 0.,  1.0, 1.0, 0.5, 'Check mid-point')
  ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. &
    test_ifndp(v, 0.0, 0.,  1.0, 1.0, 0.0, 'Check bottom')
  ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. &
    test_ifndp(v, 0.1, 0.,  1.1, 1.0, 0.0, 'Check below')
  ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. &
    test_ifndp(v,-1.0, 0.,  0.0, 1.0, 1.0, 'Check top')
  ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. &
    test_ifndp(v,-1.0, 0., -0.1, 1.0, 1.0, 'Check above')
  ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. &
    test_ifndp(v,-1.0, 0.,  3.0, 1.0, 0.25, 'Check 1/4')
  ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. &
    test_ifndp(v,-3.0, 0.,  1.0, 1.0, 0.75, 'Check 3/4')
  ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. &
    test_ifndp(v, 1.0, 0.,  1.0, 1.0, 0.0, 'Check dRho=0 below')
  ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. &
    test_ifndp(v,-1.0, 0., -1.0, 1.0, 1.0, 'Check dRho=0 above')
  ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. &
    test_ifndp(v, 0.0, 0.,  0.0, 1.0, 0.5, 'Check dRho=0 mid')
  ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. &
    test_ifndp(v,-2.0, .5,  5.0, 0.5, 0.5, 'Check dP=0')

  ! Identical columns
  call find_neutral_surface_positions_continuous(3, &
             (/0.,10.,20.,30./), (/22.,18.,14.,10./), (/0.,0.,0.,0./), & ! Left positions, T and S
             (/-1.,-1.,-1.,-1./), (/1.,1.,1.,1./), &! Left dRdT and dRdS
             (/0.,10.,20.,30./), (/22.,18.,14.,10./), (/0.,0.,0.,0./), & ! Right positions, T and S
             (/-1.,-1.,-1.,-1./), (/1.,1.,1.,1./), &! Right dRdT and dRdS
             pilro, pirlo, kol, kor, heff)
  ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or.  test_nsp(v, 8, kol, kor, pilro, pirlo, heff, &
                                   (/1,1,2,2,3,3,3,3/), & ! KoL
                                   (/1,1,2,2,3,3,3,3/), & ! KoR
                                   (/0.,0.,0.,0.,0.,0.,1.,1./), & ! pL
                                   (/0.,0.,0.,0.,0.,0.,1.,1./), & ! pR
                                   (/0.,10.,0.,10.,0.,10.,0./), & ! hEff
                                   'Identical columns')
  ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. test_data1d(v, 8, &
                                   absolute_positions(3, 8, (/0.,10.,20.,30./), kol, pilro), &
                                   (/0.,0.,10.,10.,20.,20.,30.,30./), '... left positions')
  ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. test_data1d(v, 8, &
                                   absolute_positions(3, 8, (/0.,10.,20.,30./), kor, pirlo), &
                                   (/0.,0.,10.,10.,20.,20.,30.,30./), '... right positions')
  call neutral_surface_flux(3, 2*3+2, 2, (/10.,10.,10./), (/10.,10.,10./), & ! nk, hL, hR
                               (/20.,16.,12./), (/20.,16.,12./), & ! Tl, Tr
                               pilro, pirlo, kol, kor, heff, flx, .true., h_neglect)
  ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. test_data1d(v, 7, flx, &
              (/0.,0.,0.,0.,0.,0.,0./), 'Identical columns, rho flux (=0)')
  call neutral_surface_flux(3, 2*3+2, 2, (/10.,10.,10./), (/10.,10.,10./), & ! nk, hL, hR
                               (/-1.,-1.,-1./), (/1.,1.,1./), & ! Sl, Sr
                               pilro, pirlo, kol, kor, heff, flx, .true., h_neglect)
  ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. test_data1d(v, 7, flx, &
              (/0.,20.,0.,20.,0.,20.,0./), 'Identical columns, S flux')

  ! Right column slightly cooler than left
  call find_neutral_surface_positions_continuous(3, &
             (/0.,10.,20.,30./), (/22.,18.,14.,10./), (/0.,0.,0.,0./), & ! Left positions, T and S
             (/-1.,-1.,-1.,-1./), (/1.,1.,1.,1./), &! Left dRdT and dRdS
             (/0.,10.,20.,30./), (/20.,16.,12.,8./), (/0.,0.,0.,0./), & ! Right positions, T and S
             (/-1.,-1.,-1.,-1./), (/1.,1.,1.,1./), &! Right dRdT and dRdS
             pilro, pirlo, kol, kor, heff)
  ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or.  test_nsp(v, 8, kol, kor, pilro, pirlo, heff, &
                                   (/1,1,2,2,3,3,3,3/), & ! kL
                                   (/1,1,1,2,2,3,3,3/), & ! kR
                                   (/0.,0.5,0.,0.5,0.,0.5,1.,1./), & ! pL
                                   (/0.,0.,0.5,0.,0.5,0.,0.5,1./), & ! pR
                                   (/0.,5.,5.,5.,5.,5.,0./), & ! hEff
                                   'Right column slightly cooler')
  ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. test_data1d(v, 8, &
                                   absolute_positions(3, 8, (/0.,10.,20.,30./), kol, pilro), &
                                   (/0.,5.,10.,15.,20.,25.,30.,30./), '... left positions')
  ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. test_data1d(v, 8, &
                                   absolute_positions(3, 8, (/0.,10.,20.,30./), kor, pirlo), &
                                   (/0.,0.,5.,10.,15.,20.,25.,30./), '... right positions')

  ! Right column slightly warmer than left
  call find_neutral_surface_positions_continuous(3, &
             (/0.,10.,20.,30./), (/22.,18.,14.,10./), (/0.,0.,0.,0./), & ! Left positions, T and S
             (/-1.,-1.,-1.,-1./), (/1.,1.,1.,1./), &! Left dRdT and dRdS
             (/0.,10.,20.,30./), (/24.,20.,16.,12./), (/0.,0.,0.,0./), & ! Right positions, T and S
             (/-1.,-1.,-1.,-1./), (/1.,1.,1.,1./), &! Right dRdT and dRdS
             pilro, pirlo, kol, kor, heff)
  ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or.  test_nsp(v, 8, kol, kor, pilro, pirlo, heff, &
                                   (/1,1,1,2,2,3,3,3/), & ! kL
                                   (/1,1,2,2,3,3,3,3/), & ! kR
                                   (/0.,0.,0.5,0.,0.5,0.,0.5,1./), & ! pL
                                   (/0.,0.5,0.,0.5,0.,0.5,1.,1./), & ! pR
                                   (/0.,5.,5.,5.,5.,5.,0./), & ! hEff
                                   'Right column slightly warmer')

  ! Right column somewhat cooler than left
  call find_neutral_surface_positions_continuous(3, &
             (/0.,10.,20.,30./), (/22.,18.,14.,10./), (/0.,0.,0.,0./), & ! Left positions, T and S
             (/-1.,-1.,-1.,-1./), (/1.,1.,1.,1./), &! Left dRdT and dRdS
             (/0.,10.,20.,30./), (/16.,12.,8.,4./), (/0.,0.,0.,0./), & ! Right positions, T and S
             (/-1.,-1.,-1.,-1./), (/1.,1.,1.,1./), &! Right dRdT and dRdS
             pilro, pirlo, kol, kor, heff)
  ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or.  test_nsp(v, 8, kol, kor, pilro, pirlo, heff, &
                                   (/1,2,2,3,3,3,3,3/), & ! kL
                                   (/1,1,1,1,2,2,3,3/), & ! kR
                                   (/0.,0.,0.5,0.,0.5,1.,1.,1./), & ! pL
                                   (/0.,0.,0.,0.5,0.,0.5,0.,1./), & ! pR
                                   (/0.,0.,5.,5.,5.,0.,0./), & ! hEff
                                   'Right column somewhat cooler')

  ! Right column much colder than left with no overlap
  call find_neutral_surface_positions_continuous(3, &
             (/0.,10.,20.,30./), (/22.,18.,14.,10./), (/0.,0.,0.,0./), & ! Left positions, T and S
             (/-1.,-1.,-1.,-1./), (/1.,1.,1.,1./), &! Left dRdT and dRdS
             (/0.,10.,20.,30./), (/9.,7.,5.,3./), (/0.,0.,0.,0./), & ! Right positions, T and S
             (/-1.,-1.,-1.,-1./), (/1.,1.,1.,1./), &! Right dRdT and dRdS
             pilro, pirlo, kol, kor, heff)
  ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or.  test_nsp(v, 8, kol, kor, pilro, pirlo, heff, &
                                   (/1,2,3,3,3,3,3,3/), & ! kL
                                   (/1,1,1,1,1,2,3,3/), & ! kR
                                   (/0.,0.,0.,1.,1.,1.,1.,1./), & ! pL
                                   (/0.,0.,0.,0.,0.,0.,0.,1./), & ! pR
                                   (/0.,0.,0.,0.,0.,0.,0./), & ! hEff
                                   'Right column much cooler')

  ! Right column with mixed layer
  call find_neutral_surface_positions_continuous(3, &
             (/0.,10.,20.,30./), (/22.,18.,14.,10./), (/0.,0.,0.,0./), & ! Left positions, T and S
             (/-1.,-1.,-1.,-1./), (/1.,1.,1.,1./), &! Left dRdT and dRdS
             (/0.,10.,20.,30./), (/14.,14.,10.,2./), (/0.,0.,0.,0./), & ! Right positions, T and S
             (/-1.,-1.,-1.,-1./), (/1.,1.,1.,1./), &! Right dRdT and dRdS
             pilro, pirlo, kol, kor, heff)
  ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or.  test_nsp(v, 8, kol, kor, pilro, pirlo, heff, &
                                   (/1,2,3,3,3,3,3,3/), & ! kL
                                   (/1,1,1,1,2,3,3,3/), & ! kR
                                   (/0.,0.,0.,0.,0.,1.,1.,1./), & ! pL
                                   (/0.,0.,0.,0.,0.,0.,0.,1./), & ! pR
                                   (/0.,0.,0.,0.,10.,0.,0./), & ! hEff
                                   'Right column with mixed layer')

  ! Identical columns with mixed layer
  call find_neutral_surface_positions_continuous(3, &
             (/0.,10.,20.,30./), (/14.,14.,10.,2./), (/0.,0.,0.,0./), & ! Left positions, T and S
             (/-1.,-1.,-1.,-1./), (/1.,1.,1.,1./), &! Left dRdT and dRdS
             (/0.,10.,20.,30./), (/14.,14.,10.,2./), (/0.,0.,0.,0./), & ! Right positions, T and S
             (/-1.,-1.,-1.,-1./), (/1.,1.,1.,1./), &! Right dRdT and dRdS
             pilro, pirlo, kol, kor, heff)
  ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or.  test_nsp(v, 8, kol, kor, pilro, pirlo, heff, &
                                   (/1,1,2,2,3,3,3,3/), & ! kL
                                   (/1,1,2,2,3,3,3,3/), & ! kR
                                   (/0.,0.,0.,0.,0.,0.,1.,1./), & ! pL
                                   (/0.,0.,0.,0.,0.,0.,1.,1./), & ! pR
                                   (/0.,10.,0.,10.,0.,10.,0./), & ! hEff
                                   'Identical columns with mixed layer')

  ! Right column with unstable mixed layer
  call find_neutral_surface_positions_continuous(3, &
             (/0.,10.,20.,30./), (/14.,14.,10.,2./), (/0.,0.,0.,0./), & ! Left positions, T and S
             (/-1.,-1.,-1.,-1./), (/1.,1.,1.,1./), &! Left dRdT and dRdS
             (/0.,10.,20.,30./), (/10.,14.,12.,4./), (/0.,0.,0.,0./), & ! Right positions, T and S
             (/-1.,-1.,-1.,-1./), (/1.,1.,1.,1./), &! Right dRdT and dRdS
             pilro, pirlo, kol, kor, heff)
  ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or.  test_nsp(v, 8, kol, kor, pilro, pirlo, heff, &
                                   (/1,2,3,3,3,3,3,3/), & ! kL
                                   (/1,1,1,2,3,3,3,3/), & ! kR
                                   (/0.,0.,0.,0.,0.,0.,.75,1./), & ! pL
                                   (/0.,0.,0.,0.,0.,0.25,1.,1./), & ! pR
                                   (/0.,0.,0.,0.,0.,7.5,0./), & ! hEff
                                   'Right column with unstable mixed layer')

  ! Left column with unstable mixed layer
  call find_neutral_surface_positions_continuous(3, &
             (/0.,10.,20.,30./), (/10.,14.,12.,4./), (/0.,0.,0.,0./), & ! Left positions, T and S
             (/-1.,-1.,-1.,-1./), (/1.,1.,1.,1./), &! Left dRdT and dRdS
             (/0.,10.,20.,30./), (/14.,14.,10.,2./), (/0.,0.,0.,0./), & ! Right positions, T and S
             (/-1.,-1.,-1.,-1./), (/1.,1.,1.,1./), &! Right dRdT and dRdS
             pilro, pirlo, kol, kor, heff)
  ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or.  test_nsp(v, 8, kol, kor, pilro, pirlo, heff, &
                                   (/1,1,1,2,3,3,3,3/), & ! kL
                                   (/1,2,3,3,3,3,3,3/), & ! kR
                                   (/0.,0.,0.,0.,0.,0.25,1.,1./), & ! pL
                                   (/0.,0.,0.,0.,0.,0.,.75,1./), & ! pR
                                   (/0.,0.,0.,0.,0.,7.5,0./), & ! hEff
                                   'Left column with unstable mixed layer')

  ! Two unstable mixed layers
  call find_neutral_surface_positions_continuous(3, &
             (/0.,10.,20.,30./), (/8.,12.,10.,2./), (/0.,0.,0.,0./), & ! Left positions, T and S
             (/-1.,-1.,-1.,-1./), (/1.,1.,1.,1./), &! Left dRdT and dRdS
             (/0.,10.,20.,30./), (/10.,14.,12.,4./), (/0.,0.,0.,0./), & ! Right positions, T and S
             (/-1.,-1.,-1.,-1./), (/1.,1.,1.,1./), &! Right dRdT and dRdS
             pilro, pirlo, kol, kor, heff)
  ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or.  test_nsp(v, 8, kol, kor, pilro, pirlo, heff, &
                                   (/1,1,1,1,2,3,3,3/), & ! kL
                                   (/1,2,3,3,3,3,3,3/), & ! kR
                                   (/0.,0.,0.,0.,0.,0.,0.75,1./), & ! pL
                                   (/0.,0.,0.,0.5,0.5,0.5,1.,1./), & ! pR
                                   (/0.,0.,0.,0.,0.,6.,0./), & ! hEff
                                   'Two unstable mixed layers')

  if (.not. ndiff_unit_tests_continuous) write(stdout,*) 'Pass'

ndiff_unit_tests_discontinuous()

logical function mom_neutral_diffusion::ndiff_unit_tests_discontinuous ( logical, intent(in)  verbose )

private

Parameters

[in]

verbose

It true, write results to stdout

Definition at line 2354 of file MOM_neutral_diffusion.F90.

  logical, intent(in) :: verbose !< It true, write results to stdout
  ! Local variables
  integer, parameter          :: nk = 3
  integer, parameter          :: ns = nk*4
  real, dimension(nk)         :: Sl, Sr, Tl, Tr ! Salinities [ppt] and temperatures [degC]
  real, dimension(nk)         :: hl, hr    ! Thicknesses in pressure units [R L2 T-2 ~> Pa]
  real, dimension(nk,2)       :: TiL, SiL, TiR, SiR ! Cell edge salinities [ppt] and temperatures [degC]
  real, dimension(nk,2)       :: Pres_l, Pres_r ! Interface pressures [R L2 T-2 ~> Pa]
  integer, dimension(ns)      :: KoL, KoR
  real, dimension(ns)         :: PoL, PoR
  real, dimension(ns-1)       :: hEff, Flx
  type(neutral_diffusion_CS)  :: CS        !< Neutral diffusion control structure
  type(EOS_type),     pointer :: EOS       !< Structure for linear equation of state
  type(remapping_CS), pointer :: remap_CS  !< Remapping control structure (PLM)
  real, dimension(nk,2)       :: ppoly_T_l, ppoly_T_r ! Linear reconstruction for T
  real, dimension(nk,2)       :: ppoly_S_l, ppoly_S_r ! Linear reconstruction for S
  real, dimension(nk,2)       :: dRdT      !< Partial derivative of density with temperature at
                                           !! cell edges [R degC-1 ~> kg m-3 degC-1]
  real, dimension(nk,2)       :: dRdS      !< Partial derivative of density with salinity at
                                           !! cell edges [R ppt-1 ~> kg m-3 ppt-1]
  logical, dimension(nk)      :: stable_l, stable_r
  integer                     :: iMethod
  integer                     :: ns_l, ns_r
  integer :: k
  logical :: v

  v = verbose
  ndiff_unit_tests_discontinuous = .false. ! Normally return false
  write(stdout,*) '==== MOM_neutral_diffusion: ndiff_unit_tests_discontinuous ='

  ! Unit tests for find_neutral_surface_positions_discontinuous
  ! Salinity is 0 for all these tests
  allocate(cs%EOS)
  call eos_manual_init(cs%EOS, form_of_eos=eos_linear, drho_dt=-1., drho_ds=0.)
  sl(:) = 0. ; sr(:) = 0. ; ; sil(:,:) = 0. ; sir(:,:) = 0.
  ppoly_t_l(:,:) = 0.; ppoly_t_r(:,:) = 0.
  ppoly_s_l(:,:) = 0.; ppoly_s_r(:,:) = 0.
  ! Intialize any control structures needed for unit tests
  cs%ref_pres = -1.

  hl = (/10.,10.,10./) ; hr = (/10.,10.,10./)
  pres_l(1,1) = 0. ; pres_l(1,2) = hl(1) ; pres_r(1,1) = 0. ; pres_r(1,2) = hr(1)
  do k = 2,nk
    pres_l(k,1) = pres_l(k-1,2)
    pres_l(k,2) = pres_l(k,1) + hl(k)
    pres_r(k,1) = pres_r(k-1,2)
    pres_r(k,2) = pres_r(k,1) + hr(k)
  enddo
  cs%delta_rho_form = 'mid_pressure'
  cs%neutral_pos_method = 1

  til(1,:) = (/ 22.00, 18.00 /); til(2,:) = (/ 18.00, 14.00 /); til(3,:) = (/ 14.00, 10.00 /);
  tir(1,:) = (/ 22.00, 18.00 /); tir(2,:) = (/ 18.00, 14.00 /); tir(3,:) = (/ 14.00, 10.00 /);
  call mark_unstable_cells( cs, nk, til, sil, pres_l, stable_l )
  call mark_unstable_cells( cs, nk, tir, sir, pres_r, stable_r )
  call find_neutral_surface_positions_discontinuous(cs, nk, pres_l, hl, til, sil, ppoly_t_l, ppoly_s_l, stable_l, &
           pres_r, hr, tir, sir, ppoly_t_r, ppoly_s_r, stable_r, pol, por, kol, kor, heff)
  ndiff_unit_tests_discontinuous = ndiff_unit_tests_discontinuous .or.  test_nsp(v, 12, kol, kor, pol, por, heff, &
    (/ 1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3 /),  & ! KoL
    (/ 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 3 /),  & ! KoR
    (/ 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
    (/ 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
    (/ 0.00, 10.00, 0.00, 0.00, 0.00, 10.00, 0.00, 0.00, 0.00, 10.00, 0.00 /),  & ! hEff
    'Identical Columns')

  til(1,:) = (/ 22.00, 18.00 /); til(2,:) = (/ 18.00, 14.00 /); til(3,:) = (/ 14.00, 10.00 /);
  tir(1,:) = (/ 20.00, 16.00 /); tir(2,:) = (/ 16.00, 12.00 /); tir(3,:) = (/ 12.00, 8.00 /);
  call mark_unstable_cells( cs, nk, til, sil, pres_l, stable_l )
  call mark_unstable_cells( cs, nk, tir, sir, pres_r, stable_r )
  call find_neutral_surface_positions_discontinuous(cs, nk, pres_l, hl, til, sil, ppoly_t_l, ppoly_s_l, stable_l, &
           pres_r, hr, tir, sir, ppoly_t_r, ppoly_s_r, stable_r, pol, por, kol, kor, heff)
  ndiff_unit_tests_discontinuous = ndiff_unit_tests_discontinuous .or.  test_nsp(v, 12, kol, kor, pol, por, heff, &
    (/ 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 3 /),  & ! KoL
    (/ 1, 1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3 /),  & ! KoR
    (/ 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
    (/ 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
    (/ 0.00, 5.00, 0.00, 5.00, 0.00, 5.00, 0.00, 5.00, 0.00, 5.00, 0.00 /),  & ! hEff
    'Right slightly cooler')

  til(1,:) = (/ 20.00, 16.00 /); til(2,:) = (/ 16.00, 12.00 /); til(3,:) = (/ 12.00, 8.00 /);
  tir(1,:) = (/ 22.00, 18.00 /); tir(2,:) = (/ 18.00, 14.00 /); tir(3,:) = (/ 14.00, 10.00 /);
  call mark_unstable_cells( cs, nk, til, sil, pres_l, stable_l )
  call mark_unstable_cells( cs, nk, tir, sir, pres_r, stable_r )
  call find_neutral_surface_positions_discontinuous(cs, nk, pres_l, hl, til, sil, ppoly_t_l, ppoly_s_l, stable_l, &
           pres_r, hr, tir, sir, ppoly_t_r, ppoly_s_r, stable_r, pol, por, kol, kor, heff)
  ndiff_unit_tests_discontinuous = ndiff_unit_tests_discontinuous .or.  test_nsp(v, 12, kol, kor, pol, por, heff, &
    (/ 1, 1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3 /),  & ! KoL
    (/ 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 3 /),  & ! KoR
    (/ 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
    (/ 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
    (/ 0.00, 5.00, 0.00, 5.00, 0.00, 5.00, 0.00, 5.00, 0.00, 5.00, 0.00 /),  & ! hEff
    'Left slightly cooler')

  til(1,:) = (/ 22.00, 20.00 /); til(2,:) = (/ 18.00, 16.00 /); til(3,:) = (/ 14.00, 12.00 /);
  tir(1,:) = (/ 32.00, 24.00 /); tir(2,:) = (/ 22.00, 14.00 /); tir(3,:) = (/ 12.00, 4.00 /);
  call mark_unstable_cells( cs, nk, til, sil, pres_l, stable_l )
  call mark_unstable_cells( cs, nk, tir, sir, pres_r, stable_r )
  call find_neutral_surface_positions_discontinuous(cs, nk, pres_l, hl, til, sil, ppoly_t_l, ppoly_s_l, stable_l, &
           pres_r, hr, tir, sir, ppoly_t_r, ppoly_s_r, stable_r, pol, por, kol, kor, heff)
  ndiff_unit_tests_discontinuous = ndiff_unit_tests_discontinuous .or.  test_nsp(v, 12, kol, kor, pol, por, heff, &
    (/ 1, 1, 1, 1, 1, 2, 2, 3, 3, 3, 3, 3 /),  & ! KoL
    (/ 1, 1, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3 /),  & ! KoR
    (/ 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
    (/ 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
    (/ 0.00, 0.00, 0.00, 4.00, 0.00, 4.00, 0.00, 0.00, 0.00, 0.00, 0.00 /),  & ! hEff
    'Right more strongly stratified')

  til(1,:) = (/ 22.00, 18.00 /); til(2,:) = (/ 18.00, 14.00 /); til(3,:) = (/ 14.00, 10.00 /);
  tir(1,:) = (/ 14.00, 14.00 /); tir(2,:) = (/ 14.00, 14.00 /); tir(3,:) = (/ 12.00, 8.00 /);
  call mark_unstable_cells( cs, nk, til, sil, pres_l, stable_l )
  call mark_unstable_cells( cs, nk, tir, sir, pres_r, stable_r )
  call find_neutral_surface_positions_discontinuous(cs, nk, pres_l, hl, til, sil, ppoly_t_l, ppoly_s_l, stable_l, &
           pres_r, hr, tir, sir, ppoly_t_r, ppoly_s_r, stable_r, pol, por, kol, kor, heff)
  ndiff_unit_tests_discontinuous = ndiff_unit_tests_discontinuous .or.  test_nsp(v, 12, kol, kor, pol, por, heff, &
    (/ 1, 1, 1, 1, 1, 1, 2, 2, 3, 3, 3, 3 /),  & ! KoL
    (/ 1, 1, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3 /),  & ! KoR
    (/ 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
    (/ 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
    (/ 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 5.00, 0.00 /),  & ! hEff
    'Deep Mixed layer on the right')

  til(1,:) = (/ 14.00, 14.00 /); til(2,:) = (/ 14.00, 12.00 /); til(3,:) = (/ 10.00, 8.00 /);
  tir(1,:) = (/ 14.00, 14.00 /); tir(2,:) = (/ 14.00, 14.00 /); tir(3,:) = (/ 14.00, 14.00 /);
  call mark_unstable_cells( cs, nk, til, sil, pres_l, stable_l )
  call mark_unstable_cells( cs, nk, tir, sir, pres_r, stable_r )
  call find_neutral_surface_positions_discontinuous(cs, nk, pres_l, hl, til, sil, ppoly_t_l, ppoly_s_l, stable_l, &
           pres_r, hr, tir, sir, ppoly_t_r, ppoly_s_r, stable_r, pol, por, kol, kor, heff)
  ndiff_unit_tests_discontinuous = ndiff_unit_tests_discontinuous .or.  test_nsp(v, 12, kol, kor, pol, por, heff, &
    (/ 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 3 /),  & ! KoL
    (/ 1, 1, 1, 1, 2, 2, 3, 3, 3, 3, 3, 3 /),  & ! KoR
    (/ 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
    (/ 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
    (/ 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00 /),  & ! hEff
    'Right unstratified column')

  til(1,:) = (/ 14.00, 14.00 /); til(2,:) = (/ 14.00, 12.00 /); til(3,:) = (/ 10.00, 8.00 /);
  tir(1,:) = (/ 14.00, 14.00 /); tir(2,:) = (/ 14.00, 14.00 /); tir(3,:) = (/ 12.00, 4.00 /);
  call mark_unstable_cells( cs, nk, til, sil, pres_l, stable_l )
  call mark_unstable_cells( cs, nk, tir, sir, pres_r, stable_r )
  call find_neutral_surface_positions_discontinuous(cs, nk, pres_l, hl, til, sil, ppoly_t_l, ppoly_s_l, stable_l, &
           pres_r, hr, tir, sir, ppoly_t_r, ppoly_s_r, stable_r, pol, por, kol, kor, heff)
  ndiff_unit_tests_discontinuous = ndiff_unit_tests_discontinuous .or.  test_nsp(v, 12, kol, kor, pol, por, heff, &
    (/ 1, 1, 1, 1, 1, 1, 2, 2, 2, 3, 3, 3 /),  & ! KoL
    (/ 1, 1, 1, 1, 2, 2, 3, 3, 3, 3, 3, 3 /),  & ! KoR
    (/ 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
    (/ 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
    (/ 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 4.00, 0.00 /),  & ! hEff
    'Right unstratified column')

  til(1,:) = (/ 14.00, 14.00 /); til(2,:) = (/ 14.00, 10.00 /); til(3,:) = (/ 10.00, 2.00 /);
  tir(1,:) = (/ 14.00, 14.00 /); tir(2,:) = (/ 14.00, 10.00 /); tir(3,:) = (/ 10.00, 2.00 /);
  call mark_unstable_cells( cs, nk, til, sil, pres_l, stable_l )
  call mark_unstable_cells( cs, nk, tir, sir, pres_r, stable_r )
  call find_neutral_surface_positions_discontinuous(cs, nk, pres_l, hl, til, sil, ppoly_t_l, ppoly_s_l, stable_l, &
           pres_r, hr, tir, sir, ppoly_t_r, ppoly_s_r, stable_r, pol, por, kol, kor, heff)
  ndiff_unit_tests_discontinuous = ndiff_unit_tests_discontinuous .or.  test_nsp(v, 12, kol, kor, pol, por, heff, &
    (/ 1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3 /),  & ! KoL
    (/ 1, 1, 1, 1, 2, 2, 2, 3, 3, 3, 3, 3 /),  & ! KoR
    (/ 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
    (/ 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
    (/ 0.00, 0.00, 0.00, 0.00, 0.00, 10.00, 0.00, 0.00, 0.00, 10.00, 0.00 /),  & ! hEff
    'Identical columns with mixed layer')

  til(1,:) = (/ 14.00, 12.00 /); til(2,:) = (/ 10.00, 10.00 /); til(3,:) = (/ 8.00, 2.00 /);
  tir(1,:) = (/ 14.00, 12.00 /); tir(2,:) = (/ 12.00, 8.00 /); tir(3,:) = (/ 8.00, 2.00 /);
  call mark_unstable_cells( cs, nk, til, sil, pres_l, stable_l )
  call mark_unstable_cells( cs, nk, tir, sir, pres_r, stable_r )
  call find_neutral_surface_positions_discontinuous(cs, nk, pres_l, hl, til, sil, ppoly_t_l, ppoly_s_l, stable_l, &
           pres_r, hr, tir, sir, ppoly_t_r, ppoly_s_r, stable_r, pol, por, kol, kor, heff)
  ndiff_unit_tests_discontinuous = ndiff_unit_tests_discontinuous .or.  test_nsp(v, 12, kol, kor, pol, por, heff, &
    (/ 1, 1, 1, 1, 2, 2, 3, 3, 3, 3, 3, 3 /),  & ! KoL
    (/ 1, 1, 1, 2, 2, 2, 2, 2, 2, 3, 3, 3 /),  & ! KoR
    (/ 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
    (/ 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
    (/ 0.00, 10.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 10.00, 0.00 /),  & ! hEff
    'Left interior unstratified')

  til(1,:) = (/ 12.00, 12.00 /); til(2,:) = (/ 12.00, 10.00 /); til(3,:) = (/ 10.00, 6.00 /);
  tir(1,:) = (/ 12.00, 10.00 /); tir(2,:) = (/ 10.00, 12.00 /); tir(3,:) = (/ 8.00, 4.00 /);
  call mark_unstable_cells( cs, nk, til, sil, pres_l, stable_l )
  call mark_unstable_cells( cs, nk, tir, sir, pres_r, stable_r )
  call find_neutral_surface_positions_discontinuous(cs, nk, pres_l, hl, til, sil, ppoly_t_l, ppoly_s_l, stable_l, &
           pres_r, hr, tir, sir, ppoly_t_r, ppoly_s_r, stable_r, pol, por, kol, kor, heff)
  ndiff_unit_tests_discontinuous = ndiff_unit_tests_discontinuous .or.  test_nsp(v, 12, kol, kor, pol, por, heff, &
    (/ 1, 1, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3 /),  & ! KoL
    (/ 1, 1, 1, 1, 1, 1, 2, 2, 3, 3, 3, 3 /),  & ! KoR
    (/ 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
    (/ 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
    (/ 0.00, 0.00, 0.00, 10.00, 0.00, 0.00, 0.00, 0.00, 0.00, 5.00, 0.00 /),  & ! hEff
    'Left mixed layer, Right unstable interior')

  til(1,:) = (/ 14.00, 14.00 /); til(2,:) = (/ 10.00, 10.00 /); til(3,:) = (/ 8.00, 6.00 /);
  tir(1,:) = (/ 10.00, 14.00 /); tir(2,:) = (/ 16.00, 16.00 /); tir(3,:) = (/ 12.00, 4.00 /);
  call mark_unstable_cells( cs, nk, til, sil, pres_l, stable_l )
  call mark_unstable_cells( cs, nk, tir, sir, pres_r, stable_r )
  call find_neutral_surface_positions_discontinuous(cs, nk, pres_l, hl, til, sil, ppoly_t_l, ppoly_s_l, stable_l, &
           pres_r, hr, tir, sir, ppoly_t_r, ppoly_s_r, stable_r, pol, por, kol, kor, heff)
  ndiff_unit_tests_discontinuous = ndiff_unit_tests_discontinuous .or.  test_nsp(v, 12, kol, kor, pol, por, heff, &
    (/ 1, 1, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3 /),  & ! KoL
    (/ 1, 1, 1, 1, 1, 1, 2, 2, 3, 3, 3, 3 /),  & ! KoR
    (/ 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
    (/ 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
    (/ 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 4.00, 0.00 /),  & ! hEff
    'Left thick mixed layer, Right unstable mixed')

  til(1,:) = (/ 8.00, 12.00 /); til(2,:) = (/ 12.00, 10.00 /); til(3,:) = (/ 8.00, 4.00 /);
  tir(1,:) = (/ 10.00, 14.00 /); tir(2,:) = (/ 14.00, 12.00 /); tir(3,:) = (/ 10.00, 6.00 /);
  call mark_unstable_cells( cs, nk, til, sil, pres_l, stable_l )
  call mark_unstable_cells( cs, nk, tir, sir, pres_r, stable_r )
  call find_neutral_surface_positions_discontinuous(cs, nk, pres_l, hl, til, sil, ppoly_t_l, ppoly_s_l, stable_l, &
           pres_r, hr, tir, sir, ppoly_t_r, ppoly_s_r, stable_r, pol, por, kol, kor, heff)
  ndiff_unit_tests_discontinuous = ndiff_unit_tests_discontinuous .or.  test_nsp(v, 12, kol, kor, pol, por, heff, &
    (/ 1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3 /),  & ! KoL
    (/ 1, 1, 1, 1, 2, 2, 2, 3, 3, 3, 3, 3 /),  & ! KoR
    (/ 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
    (/ 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
    (/ 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 5.00, 0.00 /),  & ! hEff
    'Unstable mixed layers, left cooler')

  call eos_manual_init(cs%EOS, form_of_eos = eos_linear, drho_dt = -1., drho_ds = 2.)
  ! Tests for linearized version of searching the layer for neutral surface position
  ! EOS linear in T, uniform alpha
  cs%max_iter = 10
  ! Unit tests require explicit initialization of tolerance
  cs%Drho_tol = 0.
  cs%x_tol = 0.
  ndiff_unit_tests_discontinuous = ndiff_unit_tests_discontinuous .or. (test_rnp(0.5, &
             find_neutral_pos_linear(cs, 0., 10., 35., -0.2, 0., &
                                     -0.2, 0., -0.2, 0.,                     &
                                     (/12.,-4./), (/34.,0./)), "Temp Uniform Linearized Alpha/Beta"))
  ! EOS linear in S, uniform beta
  ndiff_unit_tests_discontinuous = ndiff_unit_tests_discontinuous .or. (test_rnp(0.5, &
             find_neutral_pos_linear(cs, 0., 10., 35., 0., 0.8, &
                                     0., 0.8, 0., 0.8,                &
                                    (/12.,0./), (/34.,2./)), "Salt Uniform Linearized Alpha/Beta"))
  ! EOS linear in T/S, uniform alpha/beta
  ndiff_unit_tests_discontinuous = ndiff_unit_tests_discontinuous .or. (test_rnp(0.5,   &
             find_neutral_pos_linear(cs, 0., 10., 35., -0.5, 0.5,                &
                                     -0.5, 0.5, -0.5, 0.5,  &
                                     (/12.,-4./), (/34.,2./)), "Temp/salt Uniform Linearized Alpha/Beta"))
  ! EOS linear in T, insensitive to So
  ndiff_unit_tests_discontinuous = ndiff_unit_tests_discontinuous .or. (test_rnp(0.5, &
             find_neutral_pos_linear(cs, 0., 10., 35., -0.2, 0., &
                                     -0.4, 0., -0.6, 0.,  &
                                     (/12.,-4./), (/34.,0./)), "Temp stratified Linearized Alpha/Beta"))
  ! EOS linear in S, insensitive to T
  ndiff_unit_tests_discontinuous = ndiff_unit_tests_discontinuous .or. (test_rnp(0.5, &
             find_neutral_pos_linear(cs, 0., 10., 35., 0., 0.8,  &
                                      0., 1.0,  0., 0.5,  &
                                     (/12.,0./), (/34.,2./)), "Salt stratified Linearized Alpha/Beta"))
  if (.not. ndiff_unit_tests_discontinuous) write(stdout,*) 'Pass'

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]	g	Ocean grid structure
[in]	gv	ocean vertical grid structure
[in]	h	Layer thickness [H ~> m or kg m-2]
[in]	coef_x	dt * Kh * dy / dx at u-points [L2 ~> m2]
[in]	coef_y	dt * Kh * dx / dy at v-points [L2 ~> m2]
[in]	dt	Tracer time step * I_numitts [T ~> s] (I_numitts in tracer_hordiff)
	reg	Tracer registry
[in]	us	A dimensional unit scaling type
	cs	Neutral diffusion control structure

Definition at line 530 of file MOM_neutral_diffusion.F90.

  type(ocean_grid_type),                     intent(in)    :: G      !< Ocean grid structure
  type(verticalGrid_type),                   intent(in)    :: GV     !< ocean vertical grid structure
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),  intent(in)    :: h      !< Layer thickness [H ~> m or kg m-2]
  real, dimension(SZIB_(G),SZJ_(G)),         intent(in)    :: Coef_x !< dt * Kh * dy / dx at u-points [L2 ~> m2]
  real, dimension(SZI_(G),SZJB_(G)),         intent(in)    :: Coef_y !< dt * Kh * dx / dy at v-points [L2 ~> m2]
  real,                                      intent(in)    :: dt     !< Tracer time step * I_numitts [T ~> s]
                                                                     !! (I_numitts in tracer_hordiff)
  type(tracer_registry_type),                pointer       :: Reg    !< Tracer registry
  type(unit_scale_type),                     intent(in)    :: US     !< A dimensional unit scaling type
  type(neutral_diffusion_CS),                pointer       :: CS     !< Neutral diffusion control structure

  ! Local variables
  real, dimension(SZIB_(G),SZJ_(G),CS%nsurf-1) :: uFlx        ! Zonal flux of tracer [H conc ~> m conc or conc kg m-2]
  real, dimension(SZI_(G),SZJB_(G),CS%nsurf-1) :: vFlx        ! Meridional flux of tracer
                                                              ! [H conc ~> m conc or conc kg m-2]
  real, dimension(SZI_(G),SZJ_(G),G%ke)        :: tendency    ! tendency array for diagn
  real, dimension(SZI_(G),SZJ_(G))             :: tendency_2d ! depth integrated content tendency for diagn
  real, dimension(SZIB_(G),SZJ_(G))            :: trans_x_2d  ! depth integrated diffusive tracer x-transport diagn
  real, dimension(SZI_(G),SZJB_(G))            :: trans_y_2d  ! depth integrated diffusive tracer y-transport diagn
  real, dimension(G%ke)                        :: dTracer     ! change in tracer concentration due to ndiffusion

  type(tracer_type), pointer                   :: Tracer => null() ! Pointer to the current tracer

  integer :: i, j, k, m, ks, nk
  real :: Idt  ! The inverse of the time step [T-1 ~> s-1]
  real :: h_neglect, h_neglect_edge

  if (.not.cs%remap_answers_2018) then
    h_neglect = gv%H_subroundoff ; h_neglect_edge = gv%H_subroundoff
  else
    h_neglect = gv%m_to_H*1.0e-30 ; h_neglect_edge = gv%m_to_H*1.0e-10
  endif

  nk = gv%ke

  do m = 1,reg%ntr ! Loop over tracer registry

    tracer => reg%Tr(m)

    ! for diagnostics
    if (tracer%id_dfxy_conc > 0 .or. tracer%id_dfxy_cont > 0 .or. tracer%id_dfxy_cont_2d > 0 .or. &
        tracer%id_dfx_2d > 0 .or. tracer%id_dfy_2d > 0) then
      idt = 1.0 / dt
      tendency(:,:,:)  = 0.0
    endif

    uflx(:,:,:) = 0.
    vflx(:,:,:) = 0.

    ! x-flux
    do j = g%jsc,g%jec ; do i = g%isc-1,g%iec
      if (g%mask2dCu(i,j)>0.) then
        call neutral_surface_flux(nk, cs%nsurf, cs%deg, h(i,j,:), h(i+1,j,:),       &
                                  tracer%t(i,j,:), tracer%t(i+1,j,:), &
                                  cs%uPoL(i,j,:), cs%uPoR(i,j,:), &
                                  cs%uKoL(i,j,:), cs%uKoR(i,j,:), &
                                  cs%uhEff(i,j,:), uflx(i,j,:), &
                                  cs%continuous_reconstruction, h_neglect, cs%remap_CS, h_neglect_edge)
      endif
    enddo ; enddo

    ! y-flux
    do j = g%jsc-1,g%jec ; do i = g%isc,g%iec
      if (g%mask2dCv(i,j)>0.) then
        call neutral_surface_flux(nk, cs%nsurf, cs%deg, h(i,j,:), h(i,j+1,:),       &
                                  tracer%t(i,j,:), tracer%t(i,j+1,:), &
                                  cs%vPoL(i,j,:), cs%vPoR(i,j,:), &
                                  cs%vKoL(i,j,:), cs%vKoR(i,j,:), &
                                  cs%vhEff(i,j,:), vflx(i,j,:),   &
                                  cs%continuous_reconstruction, h_neglect, cs%remap_CS, h_neglect_edge)
      endif
    enddo ; enddo

    ! Update the tracer concentration from divergence of neutral diffusive flux components
    do j = g%jsc,g%jec ; do i = g%isc,g%iec
      if (g%mask2dT(i,j)>0.) then

        dtracer(:) = 0.
        do ks = 1,cs%nsurf-1
          k = cs%uKoL(i,j,ks)
          dtracer(k) = dtracer(k) + coef_x(i,j)   * uflx(i,j,ks)
          k = cs%uKoR(i-1,j,ks)
          dtracer(k) = dtracer(k) - coef_x(i-1,j) * uflx(i-1,j,ks)
          k = cs%vKoL(i,j,ks)
          dtracer(k) = dtracer(k) + coef_y(i,j)   * vflx(i,j,ks)
          k = cs%vKoR(i,j-1,ks)
          dtracer(k) = dtracer(k) - coef_y(i,j-1) * vflx(i,j-1,ks)
        enddo
        do k = 1, gv%ke
          tracer%t(i,j,k) = tracer%t(i,j,k) + dtracer(k) * &
                          ( g%IareaT(i,j) / ( h(i,j,k) + gv%H_subroundoff ) )
        enddo

        if (tracer%id_dfxy_conc > 0  .or. tracer%id_dfxy_cont > 0 .or. tracer%id_dfxy_cont_2d > 0 ) then
          do k = 1, gv%ke
            tendency(i,j,k) = dtracer(k) * g%IareaT(i,j) * idt
          enddo
        endif

      endif
    enddo ; enddo

    ! Diagnose vertically summed zonal flux, giving zonal tracer transport from ndiff.
    ! Note sign corresponds to downgradient flux convention.
    if (tracer%id_dfx_2d > 0) then
      do j = g%jsc,g%jec ; do i = g%isc-1,g%iec
        trans_x_2d(i,j) = 0.
        if (g%mask2dCu(i,j)>0.) then
          do ks = 1,cs%nsurf-1
            trans_x_2d(i,j) = trans_x_2d(i,j) - coef_x(i,j) * uflx(i,j,ks)
          enddo
          trans_x_2d(i,j) = trans_x_2d(i,j) * idt
        endif
      enddo ; enddo
      call post_data(tracer%id_dfx_2d, trans_x_2d(:,:), cs%diag)
    endif

    ! Diagnose vertically summed merid flux, giving meridional tracer transport from ndiff.
    ! Note sign corresponds to downgradient flux convention.
    if (tracer%id_dfy_2d > 0) then
      do j = g%jsc-1,g%jec ; do i = g%isc,g%iec
        trans_y_2d(i,j) = 0.
        if (g%mask2dCv(i,j)>0.) then
          do ks = 1,cs%nsurf-1
            trans_y_2d(i,j) = trans_y_2d(i,j) - coef_y(i,j) * vflx(i,j,ks)
          enddo
          trans_y_2d(i,j) = trans_y_2d(i,j) * idt
        endif
      enddo ; enddo
      call post_data(tracer%id_dfy_2d, trans_y_2d(:,:), cs%diag)
    endif

    ! post tendency of tracer content
    if (tracer%id_dfxy_cont > 0) then
      call post_data(tracer%id_dfxy_cont, tendency(:,:,:), cs%diag)
    endif

    ! post depth summed tendency for tracer content
    if (tracer%id_dfxy_cont_2d > 0) then
      tendency_2d(:,:) = 0.
      do j = g%jsc,g%jec ; do i = g%isc,g%iec
        do k = 1, gv%ke
          tendency_2d(i,j) = tendency_2d(i,j) + tendency(i,j,k)
        enddo
      enddo ; enddo
      call post_data(tracer%id_dfxy_cont_2d, tendency_2d(:,:), cs%diag)
    endif

    ! post tendency of tracer concentration; this step must be
    ! done after posting tracer content tendency, since we alter
    ! the tendency array.
    if (tracer%id_dfxy_conc > 0) then
      do k = 1, gv%ke ; do j = g%jsc,g%jec ; do i = g%isc,g%iec
        tendency(i,j,k) =  tendency(i,j,k) / ( h(i,j,k) + gv%H_subroundoff )
      enddo ; enddo ; enddo
      call post_data(tracer%id_dfxy_conc, tendency, cs%diag)
    endif
  enddo ! Loop over tracer registry

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,
		real, dimension( g %isd: g %ied, g %jsd: g %jed), intent(in), optional	p_surf
	)

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

Parameters

[in]	g	Ocean grid structure
[in]	gv	ocean vertical grid structure
[in]	us	A dimensional unit scaling type
[in]	h	Layer thickness [H ~> m or kg m-2]
[in]	t	Potential temperature [degC]
[in]	s	Salinity [ppt]
	cs	Neutral diffusion control structure
[in]	p_surf	Surface pressure to include in pressures used for equation of state calculations [R L2 T-2 ~> Pa]

Definition at line 284 of file MOM_neutral_diffusion.F90.

  type(ocean_grid_type),                    intent(in) :: G   !< Ocean grid structure
  type(verticalGrid_type),                  intent(in) :: GV  !< ocean vertical grid structure
  type(unit_scale_type),                    intent(in) :: US  !< A dimensional unit scaling type
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h   !< Layer thickness [H ~> m or kg m-2]
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: T   !< Potential temperature [degC]
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: S   !< Salinity [ppt]
  type(neutral_diffusion_CS),               pointer    :: CS  !< Neutral diffusion control structure
  real, dimension(SZI_(G),SZJ_(G)), optional, intent(in) :: p_surf !< Surface pressure to include in pressures used
                                                              !! for equation of state calculations [R L2 T-2 ~> Pa]

  ! Local variables
  integer, dimension(2) :: EOSdom ! The i-computational domain for the equation of state
  integer :: i, j, k
  ! Variables used for reconstructions
  real, dimension(SZK_(G),2) :: ppoly_r_S       ! Reconstruction slopes
  real, dimension(SZI_(G), SZJ_(G)) :: hEff_sum ! Summed effective face thicknesses [H ~> m or kg m-2]
  real, dimension(SZI_(G),SZJ_(G))  :: hbl      ! Boundary layer depth [H ~> m or kg m-2]
  integer :: iMethod
  real, dimension(SZI_(G)) :: ref_pres ! Reference pressure used to calculate alpha/beta [R L2 T-2 ~> Pa]
  real, dimension(SZI_(G)) :: rho_tmp  ! Routine to calculate drho_dp, returns density which is not used
  real :: h_neglect, h_neglect_edge    ! Negligible thicknesses [H ~> m or kg m-2]
  integer, dimension(SZI_(G), SZJ_(G)) :: k_top  ! Index of the first layer within the boundary
  real,    dimension(SZI_(G), SZJ_(G)) :: zeta_top ! Distance from the top of a layer to the intersection of the
                                                   ! top extent of the boundary layer (0 at top, 1 at bottom) [nondim]
  integer, dimension(SZI_(G), SZJ_(G)) :: k_bot    ! Index of the last layer within the boundary
  real,    dimension(SZI_(G), SZJ_(G)) :: zeta_bot ! Distance of the lower layer to the boundary layer depth
  real :: pa_to_H                      ! A conversion factor from pressure to H units [H T2 R-1 Z-2 ~> m Pa-1 or s2 m-2]

  pa_to_h = 1. / (gv%H_to_RZ * gv%g_Earth)

  k_top(:,:) = 1     ; k_bot(:,:) = 1
  zeta_top(:,:) = 0. ; zeta_bot(:,:) = 0.

  ! Check if hbl needs to be extracted
  if (cs%interior_only) then
    hbl(:,:) = 0.
    if (ASSOCIATED(cs%KPP_CSp)) call kpp_get_bld(cs%KPP_CSp, hbl, g, us, m_to_bld_units=gv%m_to_H)
    if (ASSOCIATED(cs%energetic_PBL_CSp)) &
      call energetic_pbl_get_mld(cs%energetic_PBL_CSp, hbl, g, us, m_to_mld_units=gv%m_to_H)
    call pass_var(hbl, g%Domain)
    ! get k-indices and zeta
    do j=g%jsc-1, g%jec+1 ; do i=g%isc-1,g%iec+1
      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))
    enddo; enddo
    ! TODO: add similar code for BOTTOM boundary layer
  endif

  if (.not.cs%remap_answers_2018) then
    h_neglect = gv%H_subroundoff ; h_neglect_edge = gv%H_subroundoff
  elseif (gv%Boussinesq) then
    h_neglect = gv%m_to_H*1.0e-30 ; h_neglect_edge = gv%m_to_H*1.0e-10
  else
    h_neglect = gv%kg_m2_to_H*1.0e-30 ; h_neglect_edge = gv%kg_m2_to_H*1.0e-10
  endif

  ! If doing along isopycnal diffusion (as opposed to neutral diffusion, set the reference pressure)
  if (cs%ref_pres>=0.) then
    ref_pres(:) = cs%ref_pres
  endif

  if (cs%continuous_reconstruction) then
    cs%dRdT(:,:,:) = 0.
    cs%dRdS(:,:,:) = 0.
  else
    cs%T_i(:,:,:,:) = 0.
    cs%S_i(:,:,:,:) = 0.
    cs%dRdT_i(:,:,:,:) = 0.
    cs%dRdS_i(:,:,:,:) = 0.
    cs%ns(:,:) = 0.
    cs%stable_cell(:,:,:) = .true.
  endif

  ! Calculate pressure at interfaces and layer averaged alpha/beta
  if (present(p_surf)) then
     do j=g%jsc-1,g%jec+1 ; do i=g%isc-1,g%iec+1
       cs%Pint(i,j,1) = p_surf(i,j)
     enddo ; enddo
  else
    cs%Pint(:,:,1) = 0.
  endif
  do k=1,g%ke ; do j=g%jsc-1,g%jec+1 ; do i=g%isc-1,g%iec+1
    cs%Pint(i,j,k+1) = cs%Pint(i,j,k) + h(i,j,k)*(gv%g_Earth*gv%H_to_RZ)
  enddo ; enddo ; enddo

  ! Pressures at the interfaces, this is redundant as P_i(k,1) = P_i(k-1,2) however retain this
  ! for now to ensure consitency of indexing for diiscontinuous reconstructions
  if (.not. cs%continuous_reconstruction) then
    if (present(p_surf)) then
      do j=g%jsc-1,g%jec+1 ; do i=g%isc-1,g%iec+1
        cs%P_i(i,j,1,1) = p_surf(i,j)
        cs%P_i(i,j,1,2) = p_surf(i,j) + h(i,j,1)*(gv%H_to_RZ*gv%g_Earth)
      enddo ; enddo
    else
      do j=g%jsc-1,g%jec+1 ; do i=g%isc-1,g%iec+1
        cs%P_i(i,j,1,1) = 0.
        cs%P_i(i,j,1,2) = h(i,j,1)*(gv%H_to_RZ*gv%g_Earth)
      enddo ; enddo
    endif
    do k=2,g%ke ; do j=g%jsc-1,g%jec+1 ; do i=g%isc-1,g%iec+1
      cs%P_i(i,j,k,1) = cs%P_i(i,j,k-1,2)
      cs%P_i(i,j,k,2) = cs%P_i(i,j,k-1,2) + h(i,j,k)*(gv%H_to_RZ*gv%g_Earth)
    enddo ; enddo ; enddo
  endif

  eosdom(:) = eos_domain(g%HI, halo=1)
  do j = g%jsc-1, g%jec+1
    ! Interpolate state to interface
    do i = g%isc-1, g%iec+1
      if (cs%continuous_reconstruction) then
        call interface_scalar(g%ke, h(i,j,:), t(i,j,:), cs%Tint(i,j,:), 2, h_neglect)
        call interface_scalar(g%ke, h(i,j,:), s(i,j,:), cs%Sint(i,j,:), 2, h_neglect)
      else
        call build_reconstructions_1d( cs%remap_CS, g%ke, h(i,j,:), t(i,j,:), cs%ppoly_coeffs_T(i,j,:,:), &
                                       cs%T_i(i,j,:,:), ppoly_r_s, imethod, h_neglect, h_neglect_edge )
        call build_reconstructions_1d( cs%remap_CS, g%ke, h(i,j,:), s(i,j,:), cs%ppoly_coeffs_S(i,j,:,:), &
                                       cs%S_i(i,j,:,:), ppoly_r_s, imethod, h_neglect, h_neglect_edge )
        ! In the current ALE formulation, interface values are not exactly at the 0. or 1. of the
        ! polynomial reconstructions
        do k=1,g%ke
           cs%T_i(i,j,k,1) = evaluation_polynomial( cs%ppoly_coeffs_T(i,j,k,:), cs%deg+1, 0. )
           cs%T_i(i,j,k,2) = evaluation_polynomial( cs%ppoly_coeffs_T(i,j,k,:), cs%deg+1, 1. )
           cs%S_i(i,j,k,1) = evaluation_polynomial( cs%ppoly_coeffs_S(i,j,k,:), cs%deg+1, 0. )
           cs%S_i(i,j,k,2) = evaluation_polynomial( cs%ppoly_coeffs_S(i,j,k,:), cs%deg+1, 1. )
        enddo
      endif
    enddo

    ! Continuous reconstruction
    if (cs%continuous_reconstruction) then
      do k = 1, g%ke+1
        if (cs%ref_pres<0) ref_pres(:) = cs%Pint(:,j,k)
        call calculate_density_derivs(cs%Tint(:,j,k), cs%Sint(:,j,k), ref_pres, cs%dRdT(:,j,k), &
                                      cs%dRdS(:,j,k), cs%EOS, eosdom)
      enddo
    else ! Discontinuous reconstruction
      do k = 1, g%ke
        if (cs%ref_pres<0) ref_pres(:) = cs%Pint(:,j,k)
        ! Calculate derivatives for the top interface
        call calculate_density_derivs(cs%T_i(:,j,k,1), cs%S_i(:,j,k,1), ref_pres, cs%dRdT_i(:,j,k,1), &
                                      cs%dRdS_i(:,j,k,1), cs%EOS, eosdom)
        if (cs%ref_pres<0) ref_pres(:) = cs%Pint(:,j,k+1)
        ! Calculate derivatives at the bottom interface
        call calculate_density_derivs(cs%T_i(:,j,k,2), cs%S_i(:,j,k,2), ref_pres, cs%dRdT_i(:,j,k,2), &
                                      cs%dRdS_i(:,j,k,2), cs%EOS, eosdom)
      enddo
    endif
  enddo

  if (.not. cs%continuous_reconstruction) then
    do j = g%jsc-1, g%jec+1 ; do i = g%isc-1, g%iec+1
      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,:) )
      if (cs%interior_only) then
        if (.not. cs%stable_cell(i,j,k_bot(i,j))) zeta_bot(i,j) = -1.
        ! set values in the surface and bottom boundary layer to false.
        do k = 1, k_bot(i,j)-1
          cs%stable_cell(i,j,k) = .false.
        enddo
      endif
    enddo ; enddo
  endif

  cs%uhEff(:,:,:) = 0.
  cs%vhEff(:,:,:) = 0.
  cs%uPoL(:,:,:) = 0.
  cs%vPoL(:,:,:) = 0.
  cs%uPoR(:,:,:) = 0.
  cs%vPoR(:,:,:) = 0.
  cs%uKoL(:,:,:) = 1
  cs%vKoL(:,:,:) = 1
  cs%uKoR(:,:,:) = 1
  cs%vKoR(:,:,:) = 1

  ! Neutral surface factors at U points
  do j = g%jsc, g%jec ; do i = g%isc-1, g%iec
    if (g%mask2dCu(i,j) > 0.) then
      if (cs%continuous_reconstruction) then
        call find_neutral_surface_positions_continuous(g%ke,                                    &
                cs%Pint(i,j,:), cs%Tint(i,j,:), cs%Sint(i,j,:), cs%dRdT(i,j,:), cs%dRdS(i,j,:),            &
                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,:),  &
                cs%uPoL(i,j,:), cs%uPoR(i,j,:), cs%uKoL(i,j,:), cs%uKoR(i,j,:), cs%uhEff(i,j,:),           &
                k_bot(i,j), k_bot(i+1,j), zeta_bot(i,j), zeta_bot(i+1,j))
      else
        call find_neutral_surface_positions_discontinuous(cs, g%ke, &
            cs%P_i(i,j,:,:), h(i,j,:), cs%T_i(i,j,:,:), cs%S_i(i,j,:,:), cs%ppoly_coeffs_T(i,j,:,:),           &
            cs%ppoly_coeffs_S(i,j,:,:),cs%stable_cell(i,j,:),                                                  &
            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,:,:), &
            cs%ppoly_coeffs_S(i+1,j,:,:), cs%stable_cell(i+1,j,:),                                             &
            cs%uPoL(i,j,:), cs%uPoR(i,j,:), cs%uKoL(i,j,:), cs%uKoR(i,j,:), cs%uhEff(i,j,:),                   &
            hard_fail_heff = cs%hard_fail_heff)
      endif
    endif
  enddo ; enddo

  ! Neutral surface factors at V points
  do j = g%jsc-1, g%jec ; do i = g%isc, g%iec
    if (g%mask2dCv(i,j) > 0.) then
      if (cs%continuous_reconstruction) then
        call find_neutral_surface_positions_continuous(g%ke,                                              &
                cs%Pint(i,j,:), cs%Tint(i,j,:), cs%Sint(i,j,:), cs%dRdT(i,j,:), cs%dRdS(i,j,:),           &
                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,:), &
                cs%vPoL(i,j,:), cs%vPoR(i,j,:), cs%vKoL(i,j,:), cs%vKoR(i,j,:), cs%vhEff(i,j,:), &
                k_bot(i,j), k_bot(i,j+1), zeta_bot(i,j), zeta_bot(i,j+1))
      else
        call find_neutral_surface_positions_discontinuous(cs, g%ke, &
            cs%P_i(i,j,:,:), h(i,j,:), cs%T_i(i,j,:,:), cs%S_i(i,j,:,:), cs%ppoly_coeffs_T(i,j,:,:),           &
            cs%ppoly_coeffs_S(i,j,:,:),cs%stable_cell(i,j,:),                                                  &
            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,:,:), &
            cs%ppoly_coeffs_S(i,j+1,:,:), cs%stable_cell(i,j+1,:),                                             &
            cs%vPoL(i,j,:), cs%vPoR(i,j,:), cs%vKoL(i,j,:), cs%vKoR(i,j,:), cs%vhEff(i,j,:),                   &
            hard_fail_heff = cs%hard_fail_heff)
      endif
    endif
  enddo ; enddo

  ! Continuous reconstructions calculate hEff as the difference between the pressures of the
  ! neutral surfaces which need to be reconverted to thickness units. The discontinuous version
  ! calculates hEff from the nondimensional fraction of the layer spanned by adjacent neutral
  ! surfaces, so hEff is already in thickness units.
  if (cs%continuous_reconstruction) then
    do k = 1, cs%nsurf-1 ; do j = g%jsc, g%jec ; do i = g%isc-1, g%iec
      if (g%mask2dCu(i,j) > 0.) cs%uhEff(i,j,k) = cs%uhEff(i,j,k) * pa_to_h
    enddo ; enddo ; enddo
    do k = 1, cs%nsurf-1 ; do j = g%jsc-1, g%jec ; do i = g%isc, g%iec
      if (g%mask2dCv(i,j) > 0.) cs%vhEff(i,j,k) = cs%vhEff(i,j,k) * pa_to_h
    enddo ; enddo ; enddo
  endif

  if (cs%id_uhEff_2d>0) then
    heff_sum(:,:) = 0.
    do k = 1,cs%nsurf-1 ; do j=g%jsc,g%jec ; do i=g%isc-1,g%iec
      heff_sum(i,j) = heff_sum(i,j) + cs%uhEff(i,j,k)
    enddo ; enddo ; enddo
    call post_data(cs%id_uhEff_2d, heff_sum, cs%diag)
  endif
  if (cs%id_vhEff_2d>0) then
    heff_sum(:,:) = 0.
    do k = 1,cs%nsurf-1 ; do j=g%jsc-1,g%jec ; do i=g%isc,g%iec
      heff_sum(i,j) = heff_sum(i,j) + cs%vhEff(i,j,k)
    enddo ; enddo ; enddo
    call post_data(cs%id_vhEff_2d, heff_sum, cs%diag)
  endif

neutral_diffusion_end()

subroutine, public mom_neutral_diffusion::neutral_diffusion_end

(

type(neutral_diffusion_cs), pointer

CS

)

Deallocates neutral_diffusion control structure.

Parameters

cs	Neutral diffusion control structure

Definition at line 2863 of file MOM_neutral_diffusion.F90.

  type(neutral_diffusion_CS), pointer :: CS  !< Neutral diffusion control structure

  if (associated(cs)) deallocate(cs)

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(unit_scale_type), intent(in)	US,
		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]	time	Time structure
[in]	g	Grid structure
[in]	us	A dimensional unit scaling type
[in,out]	diag	Diagnostics control structure
[in]	param_file	Parameter file structure
[in]	eos	Equation of state
	diabatic_csp	KPP control structure needed to get BLD
	cs	Neutral diffusion control structure

Definition at line 117 of file MOM_neutral_diffusion.F90.

  type(time_type), target,    intent(in)    :: Time       !< Time structure
  type(ocean_grid_type),      intent(in)    :: G          !< Grid structure
  type(unit_scale_type),      intent(in)    :: US         !< A dimensional unit scaling type
  type(diag_ctrl), target,    intent(inout) :: diag       !< Diagnostics control structure
  type(param_file_type),      intent(in)    :: param_file !< Parameter file structure
  type(EOS_type),  target,    intent(in)    :: EOS        !< Equation of state
  type(diabatic_CS),          pointer       :: diabatic_CSp!< KPP control structure needed to get BLD
  type(neutral_diffusion_CS), pointer       :: CS         !< Neutral diffusion control structure

  ! Local variables
  character(len=256) :: mesg    ! Message for error messages.
  character(len=80)  :: string  ! Temporary strings
  logical :: default_2018_answers
  logical :: boundary_extrap

  if (associated(cs)) then
    call mom_error(fatal, "neutral_diffusion_init called with associated control structure.")
    return
  endif

  ! Log this module and master switch for turning it on/off
  call get_param(param_file, mdl, "USE_NEUTRAL_DIFFUSION", neutral_diffusion_init, &
                 default=.false., do_not_log=.true.)
  call log_version(param_file, mdl, version, &
           "This module implements neutral diffusion of tracers", &
           all_default=.not.neutral_diffusion_init)
  call get_param(param_file, mdl, "USE_NEUTRAL_DIFFUSION", neutral_diffusion_init, &
                 "If true, enables the neutral diffusion module.", &
                 default=.false.)

  if (.not.neutral_diffusion_init) return

  allocate(cs)
  cs%diag => diag
  cs%EOS => eos
 ! call openParameterBlock(param_file,'NEUTRAL_DIFF')

  ! Read all relevant parameters and write them to the model log.
  call get_param(param_file, mdl, "NDIFF_CONTINUOUS", cs%continuous_reconstruction, &
                 "If true, uses a continuous reconstruction of T and S when "//&
                 "finding neutral surfaces along which diffusion will happen. "//&
                 "If false, a PPM discontinuous reconstruction of T and S "//&
                 "is done which results in a higher order routine but exacts "//&
                 "a higher computational cost.", default=.true.)
  call get_param(param_file, mdl, "NDIFF_REF_PRES", cs%ref_pres,                    &
                 "The reference pressure (Pa) used for the derivatives of "//&
                 "the equation of state. If negative (default), local pressure is used.", &
                 units="Pa", default = -1., scale=us%kg_m3_to_R*us%m_s_to_L_T**2)
  call get_param(param_file, mdl, "NDIFF_INTERIOR_ONLY", cs%interior_only, &
                 "If true, only applies neutral diffusion in the ocean interior."//&
                 "That is, the algorithm will exclude the surface and bottom"//&
                 "boundary layers.", default = .false.)

  ! Initialize and configure remapping
  if ( .not.cs%continuous_reconstruction ) then
    call get_param(param_file, mdl, "NDIFF_BOUNDARY_EXTRAP", boundary_extrap, &
                   "Extrapolate at the top and bottommost cells, otherwise   \n"//  &
                   "assume boundaries are piecewise constant",                      &
                   default=.false.)
    call get_param(param_file, mdl, "NDIFF_REMAPPING_SCHEME", string, &
                   "This sets the reconstruction scheme used "//&
                   "for vertical remapping for all variables. "//&
                   "It can be one of the following schemes: "//&
                   trim(remappingschemesdoc), default=remappingdefaultscheme)
    call get_param(param_file, mdl, "DEFAULT_2018_ANSWERS", default_2018_answers, &
                 "This sets the default value for the various _2018_ANSWERS parameters.", &
                 default=.false.)
    call get_param(param_file, mdl, "REMAPPING_2018_ANSWERS", cs%remap_answers_2018, &
                 "If true, use the order of arithmetic and expressions that recover the "//&
                 "answers from the end of 2018.  Otherwise, use updated and more robust "//&
                 "forms of the same expressions.", default=default_2018_answers)
    call initialize_remapping( cs%remap_CS, string, boundary_extrapolation=boundary_extrap, &
                               answers_2018=cs%remap_answers_2018 )
    call extract_member_remapping_cs(cs%remap_CS, degree=cs%deg)
    call get_param(param_file, mdl, "NEUTRAL_POS_METHOD", cs%neutral_pos_method,   &
                   "Method used to find the neutral position                 \n"// &
                   "1. Delta_rho varies linearly, find 0 crossing            \n"// &
                   "2. Alpha and beta vary linearly from top to bottom,      \n"// &
                   "   Newton's method for neutral position                  \n"// &
                   "3. Full nonlinear equation of state, use regula falsi    \n"// &
                   "   for neutral position", default=3)
    if (cs%neutral_pos_method > 4 .or. cs%neutral_pos_method < 0) then
      call mom_error(fatal,"Invalid option for NEUTRAL_POS_METHOD")
    endif

    call get_param(param_file, mdl, "DELTA_RHO_FORM", cs%delta_rho_form,           &
                   "Determine how the difference in density is calculated    \n"// &
                   "  full       : Difference of in-situ densities           \n"// &
                   "  no_pressure: Calculated from dRdT, dRdS, but no        \n"// &
                   "               pressure dependence",                           &
                   default="mid_pressure")
    if (cs%neutral_pos_method > 1) then
      call get_param(param_file, mdl, "NDIFF_DRHO_TOL", cs%drho_tol,            &
                     "Sets the convergence criterion for finding the neutral\n"// &
                     "position within a layer in kg m-3.",                        &
                     default=1.e-10, scale=us%kg_m3_to_R)
      call get_param(param_file, mdl, "NDIFF_X_TOL", cs%x_tol,            &
                     "Sets the convergence criterion for a change in nondim\n"// &
                     "position within a layer.",                        &
                     default=0.)
      call get_param(param_file, mdl, "NDIFF_MAX_ITER", cs%max_iter,              &
                    "The maximum number of iterations to be done before \n"//     &
                     "exiting the iterative loop to find the neutral surface",    &
                     default=10)
    endif
    call get_param(param_file, mdl, "NDIFF_DEBUG", cs%debug,             &
                   "Turns on verbose output for discontinuous neutral "//&
                   "diffusion routines.", &
                   default = .false.)
    call get_param(param_file, mdl, "HARD_FAIL_HEFF", cs%hard_fail_heff, &
                  "Bring down the model if a problem with heff is detected",&
                   default = .true.)
  endif

  ! Store a rescaling factor for use in diagnostic messages.
  cs%R_to_kg_m3 = us%R_to_kg_m3

  if (cs%interior_only) then
    call extract_diabatic_member(diabatic_csp, kpp_csp=cs%KPP_CSp)
    call extract_diabatic_member(diabatic_csp, energetic_pbl_csp=cs%energetic_PBL_CSp)
    if ( .not. ASSOCIATED(cs%energetic_PBL_CSp) .and. .not. ASSOCIATED(cs%KPP_CSp) ) then
      call mom_error(fatal,"NDIFF_INTERIOR_ONLY is true, but no valid boundary layer scheme was found")
    endif
  endif

! call get_param(param_file, mdl, "KHTR", CS%KhTr, &
!                "The background along-isopycnal tracer diffusivity.", &
!                units="m2 s-1", default=0.0)
!  call closeParameterBlock(param_file)
  if (cs%continuous_reconstruction) then
    cs%nsurf = 2*g%ke+2 ! Continuous reconstruction means that every interface has two connections
    allocate(cs%dRdT(szi_(g),szj_(g),szk_(g)+1)) ; cs%dRdT(:,:,:) = 0.
    allocate(cs%dRdS(szi_(g),szj_(g),szk_(g)+1)) ; cs%dRdS(:,:,:) = 0.
  else
    cs%nsurf = 4*g%ke   ! Discontinuous means that every interface has four connections
    allocate(cs%T_i(szi_(g),szj_(g),szk_(g),2))    ; cs%T_i(:,:,:,:) = 0.
    allocate(cs%S_i(szi_(g),szj_(g),szk_(g),2))    ; cs%S_i(:,:,:,:) = 0.
    allocate(cs%P_i(szi_(g),szj_(g),szk_(g),2))    ; cs%P_i(:,:,:,:) = 0.
    allocate(cs%dRdT_i(szi_(g),szj_(g),szk_(g),2)) ; cs%dRdT_i(:,:,:,:) = 0.
    allocate(cs%dRdS_i(szi_(g),szj_(g),szk_(g),2)) ; cs%dRdS_i(:,:,:,:) = 0.
    allocate(cs%ppoly_coeffs_T(szi_(g),szj_(g),szk_(g),cs%deg+1)) ; cs%ppoly_coeffs_T(:,:,:,:) = 0.
    allocate(cs%ppoly_coeffs_S(szi_(g),szj_(g),szk_(g),cs%deg+1)) ; cs%ppoly_coeffs_S(:,:,:,:) = 0.
    allocate(cs%ns(szi_(g),szj_(g)))    ; cs%ns(:,:) = 0.
  endif
  ! T-points
  allocate(cs%Tint(szi_(g),szj_(g),szk_(g)+1)) ; cs%Tint(:,:,:) = 0.
  allocate(cs%Sint(szi_(g),szj_(g),szk_(g)+1)) ; cs%Sint(:,:,:) = 0.
  allocate(cs%Pint(szi_(g),szj_(g),szk_(g)+1)) ; cs%Pint(:,:,:) = 0.
  allocate(cs%stable_cell(szi_(g),szj_(g),szk_(g))) ; cs%stable_cell(:,:,:) = .true.
  ! U-points
  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.
  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.
  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
  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
  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
  ! V-points
  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.
  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.
  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
  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
  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

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]

verbose

If true, write results to stdout

Definition at line 2080 of file MOM_neutral_diffusion.F90.

  logical, intent(in) :: verbose !< If true, write results to stdout

  neutral_diffusion_unit_tests = .false. .or. &
    ndiff_unit_tests_continuous(verbose) .or. ndiff_unit_tests_discontinuous(verbose)

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]	nk	Number of levels
[in]	nsurf	Number of neutral surfaces
[in]	deg	Degree of polynomial reconstructions
[in]	hl	Left-column layer thickness [H ~> m or kg m-2]
[in]	hr	Right-column layer thickness [H ~> m or kg m-2]
[in]	tl	Left-column layer tracer (conc, e.g. degC)
[in]	tr	Right-column layer tracer (conc, e.g. degC)
[in]	pil	Fractional position of neutral surface within layer KoL of left column
[in]	pir	Fractional position of neutral surface within layer KoR of right column
[in]	kol	Index of first left interface above neutral surface
[in]	kor	Index of first right interface above neutral surface
[in]	heff	Effective thickness between two neutral surfaces [H ~> m or kg m-2]
[in,out]	flx	Flux of tracer between pairs of neutral layers (conc H)
[in]	continuous	True if using continuous reconstruction
[in]	h_neglect	A negligibly small width for the purpose of cell reconstructions [H ~> m or kg m-2]
[in]	remap_cs	Remapping control structure used to create sublayers
[in]	h_neglect_edge	A negligibly small width used for edge value calculations if continuous is false [H ~> m or kg m-2]

Definition at line 1877 of file MOM_neutral_diffusion.F90.

  integer,                      intent(in)    :: nk    !< Number of levels
  integer,                      intent(in)    :: nsurf !< Number of neutral surfaces
  integer,                      intent(in)    :: deg   !< Degree of polynomial reconstructions
  real, dimension(nk),          intent(in)    :: hl    !< Left-column layer thickness [H ~> m or kg m-2]
  real, dimension(nk),          intent(in)    :: hr    !< Right-column layer thickness [H ~> m or kg m-2]
  real, dimension(nk),          intent(in)    :: Tl    !< Left-column layer tracer (conc, e.g. degC)
  real, dimension(nk),          intent(in)    :: Tr    !< Right-column layer tracer (conc, e.g. degC)
  real, dimension(nsurf),       intent(in)    :: PiL   !< Fractional position of neutral surface
                                                       !! within layer KoL of left column
  real, dimension(nsurf),       intent(in)    :: PiR   !< Fractional position of neutral surface
                                                       !! within layer KoR of right column
  integer, dimension(nsurf),    intent(in)    :: KoL   !< Index of first left interface above neutral surface
  integer, dimension(nsurf),    intent(in)    :: KoR   !< Index of first right interface above neutral surface
  real, dimension(nsurf-1),     intent(in)    :: hEff  !< Effective thickness between two neutral
                                                       !! surfaces [H ~> m or kg m-2]
  real, dimension(nsurf-1),     intent(inout) :: Flx   !< Flux of tracer between pairs of neutral layers (conc H)
  logical,                      intent(in)    :: continuous !< True if using continuous reconstruction
  real,                         intent(in)    :: h_neglect !< A negligibly small width for the
                                             !! purpose of cell reconstructions [H ~> m or kg m-2]
  type(remapping_CS), optional, intent(in)    :: remap_CS !< Remapping control structure used
                                             !! to create sublayers
  real,               optional, intent(in)    :: h_neglect_edge !< A negligibly small width used for
                                             !! edge value calculations if continuous is false [H ~> m or kg m-2]
  ! Local variables
  integer :: k_sublayer, klb, klt, krb, krt, k
  real :: T_right_top, T_right_bottom, T_right_layer, T_right_sub, T_right_top_int, T_right_bot_int
  real :: T_left_top, T_left_bottom, T_left_layer, T_left_sub, T_left_top_int, T_left_bot_int
  real :: dT_top, dT_bottom, dT_layer, dT_ave, dT_sublayer, dT_top_int, dT_bot_int
  real, dimension(nk+1) :: Til !< Left-column interface tracer (conc, e.g. degC)
  real, dimension(nk+1) :: Tir !< Right-column interface tracer (conc, e.g. degC)
  real, dimension(nk) :: aL_l !< Left-column left edge value of tracer (conc, e.g. degC)
  real, dimension(nk) :: aR_l !< Left-column right edge value of tracer (conc, e.g. degC)
  real, dimension(nk) :: aL_r !< Right-column left edge value of tracer (conc, e.g. degC)
  real, dimension(nk) :: aR_r !< Right-column right edge value of tracer (conc, e.g. degC)
  ! Discontinuous reconstruction
  integer               :: iMethod
  real, dimension(nk,2) :: Tid_l !< Left-column interface tracer (conc, e.g. degC)
  real, dimension(nk,2) :: Tid_r !< Right-column interface tracer (conc, e.g. degC)
  real, dimension(nk,deg+1) :: ppoly_r_coeffs_l
  real, dimension(nk,deg+1) :: ppoly_r_coeffs_r
  real, dimension(nk,deg+1) :: ppoly_r_S_l
  real, dimension(nk,deg+1) :: ppoly_r_S_r
  logical :: down_flux
  ! Setup reconstruction edge values
  if (continuous) then
    call interface_scalar(nk, hl, tl, til, 2, h_neglect)
    call interface_scalar(nk, hr, tr, tir, 2, h_neglect)
    call ppm_left_right_edge_values(nk, tl, til, al_l, ar_l)
    call ppm_left_right_edge_values(nk, tr, tir, al_r, ar_r)
  else
    ppoly_r_coeffs_l(:,:) = 0.
    ppoly_r_coeffs_r(:,:) = 0.
    tid_l(:,:) = 0.
    tid_r(:,:) = 0.

    call build_reconstructions_1d( remap_cs, nk, hl, tl, ppoly_r_coeffs_l, tid_l, &
                                   ppoly_r_s_l, imethod, h_neglect, h_neglect_edge )
    call build_reconstructions_1d( remap_cs, nk, hr, tr, ppoly_r_coeffs_r, tid_r, &
                                   ppoly_r_s_r, imethod, h_neglect, h_neglect_edge )
  endif

  do k_sublayer = 1, nsurf-1
    if (heff(k_sublayer) == 0.) then
      flx(k_sublayer) = 0.
    else
      if (continuous) then
        klb = kol(k_sublayer+1)
        t_left_bottom = ( 1. - pil(k_sublayer+1) ) * til(klb) + pil(k_sublayer+1) * til(klb+1)
        klt = kol(k_sublayer)
        t_left_top = ( 1. - pil(k_sublayer) ) * til(klt) + pil(k_sublayer) * til(klt+1)
        t_left_layer = ppm_ave(pil(k_sublayer), pil(k_sublayer+1) + real(klb-klt), &
                               al_l(klt), ar_l(klt), tl(klt))

        krb = kor(k_sublayer+1)
        t_right_bottom = ( 1. - pir(k_sublayer+1) ) * tir(krb) + pir(k_sublayer+1) * tir(krb+1)
        krt = kor(k_sublayer)
        t_right_top = ( 1. - pir(k_sublayer) ) * tir(krt) + pir(k_sublayer) * tir(krt+1)
        t_right_layer = ppm_ave(pir(k_sublayer), pir(k_sublayer+1) + real(krb-krt), &
                                al_r(krt), ar_r(krt), tr(krt))
        dt_top = t_right_top - t_left_top
        dt_bottom = t_right_bottom - t_left_bottom
        dt_ave = 0.5 * ( dt_top + dt_bottom )
        dt_layer = t_right_layer - t_left_layer
        if (signum(1.,dt_top) * signum(1.,dt_bottom) <= 0. .or. signum(1.,dt_ave) * signum(1.,dt_layer) <= 0.) then
          dt_ave = 0.
        else
          dt_ave = dt_layer
        endif
        flx(k_sublayer) = dt_ave * heff(k_sublayer)
      else ! Discontinuous reconstruction
        ! Calculate tracer values on left and right side of the neutral surface
        call neutral_surface_t_eval(nk, nsurf, k_sublayer, kol, pil, tl, tid_l, deg, imethod, &
                                    ppoly_r_coeffs_l, t_left_top, t_left_bottom, t_left_sub, &
                                    t_left_top_int, t_left_bot_int, t_left_layer)
        call neutral_surface_t_eval(nk, nsurf, k_sublayer, kor, pir, tr, tid_r, deg, imethod, &
                                    ppoly_r_coeffs_r, t_right_top, t_right_bottom, t_right_sub, &
                                    t_right_top_int, t_right_bot_int, t_right_layer)

        dt_top      = t_right_top     - t_left_top
        dt_bottom   = t_right_bottom  - t_left_bottom
        dt_sublayer = t_right_sub     - t_left_sub
        dt_top_int  = t_right_top_int - t_left_top_int
        dt_bot_int  = t_right_bot_int - t_left_bot_int
        ! Enforcing the below criterion incorrectly zero out fluxes
        !dT_layer = T_right_layer - T_left_layer

        down_flux = dt_top <= 0. .and. dt_bottom <= 0. .and.       &
                    dt_sublayer <= 0. .and. dt_top_int <= 0. .and. &
                    dt_bot_int <= 0.
        down_flux = down_flux .or.                                 &
                    (dt_top >= 0. .and. dt_bottom >= 0. .and.      &
                    dt_sublayer >= 0. .and. dt_top_int >= 0. .and. &
                    dt_bot_int >= 0.)
        if (down_flux) then
          flx(k_sublayer) = dt_sublayer * heff(k_sublayer)
        else
          flx(k_sublayer) = 0.
        endif
      endif
    endif
  enddo

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]	nk	Number of cell everages
[in]	ns	Number of neutral surfaces
[in]	k_sub	Index of current neutral layer
[in]	ks	List of the layers associated with each neutral surface
[in]	ps	List of the positions within a layer of each surface
[in]	t_mean	Cell average of tracer
[in]	t_int	Cell interface values of tracer from reconstruction
[in]	deg	Degree of reconstruction polynomial (e.g. 1 is linear)
[in]	imethod	Method of integration to use
[in]	t_poly	Coefficients of polynomial reconstructions
[out]	t_top	Tracer value at top (across discontinuity if necessary)
[out]	t_bot	Tracer value at bottom (across discontinuity if necessary)
[out]	t_sub	Average of the tracer value over the sublayer
[out]	t_top_int	Tracer value at top interface of neutral layer
[out]	t_bot_int	Tracer value at bottom interface of neutral layer
[out]	t_layer	Cell-average that the the reconstruction belongs to

Definition at line 2004 of file MOM_neutral_diffusion.F90.

  integer,                   intent(in   ) :: nk        !< Number of cell everages
  integer,                   intent(in   ) :: ns        !< Number of neutral surfaces
  integer,                   intent(in   ) :: k_sub     !< Index of current neutral layer
  integer, dimension(ns),    intent(in   ) :: Ks        !< List of the layers associated with each neutral surface
  real, dimension(ns),       intent(in   ) :: Ps        !< List of the positions within a layer of each surface
  real, dimension(nk),       intent(in   ) :: T_mean    !< Cell average of tracer
  real, dimension(nk,2),     intent(in   ) :: T_int     !< Cell interface values of tracer from reconstruction
  integer,                   intent(in   ) :: deg       !< Degree of reconstruction polynomial (e.g. 1 is linear)
  integer,                   intent(in   ) :: iMethod   !< Method of integration to use
  real, dimension(nk,deg+1), intent(in   ) :: T_poly    !< Coefficients of polynomial reconstructions
  real,                      intent(  out) :: T_top     !< Tracer value at top (across discontinuity if necessary)
  real,                      intent(  out) :: T_bot     !< Tracer value at bottom (across discontinuity if necessary)
  real,                      intent(  out) :: T_sub     !< Average of the tracer value over the sublayer
  real,                      intent(  out) :: T_top_int !< Tracer value at top interface of neutral layer
  real,                      intent(  out) :: T_bot_int !< Tracer value at bottom interface of neutral layer
  real,                      intent(  out) :: T_layer   !< Cell-average that the the reconstruction belongs to

  integer :: kl, ks_top, ks_bot

  ks_top = k_sub
  ks_bot = k_sub + 1
  if ( ks(ks_top) /= ks(ks_bot) ) then
    call mom_error(fatal, "Neutral surfaces span more than one layer")
  endif
  kl = ks(k_sub)
  ! First if the neutral surfaces spans the entirety of a cell, then do not search across the discontinuity
  if ( (ps(ks_top) == 0.) .and. (ps(ks_bot) == 1.)) then
    t_top = t_int(kl,1)
    t_bot = t_int(kl,2)
  else
    ! Search across potential discontinuity at top
    if ( (kl > 1) .and. (ps(ks_top) == 0.)  ) then
      t_top = t_int(kl-1,2)
    else
      t_top = evaluation_polynomial( t_poly(kl,:), deg+1, ps(ks_top) )
    endif
    ! Search across potential discontinuity at bottom
    if ( (kl < nk) .and. (ps(ks_bot) == 1.) ) then
      t_bot = t_int(kl+1,1)
    else
      t_bot = evaluation_polynomial( t_poly(kl,:), deg+1, ps(ks_bot) )
    endif
  endif
  t_sub = average_value_ppoly(nk, t_mean, t_int, t_poly, imethod, kl, ps(ks_top), ps(ks_bot))
  t_top_int = evaluation_polynomial( t_poly(kl,:), deg+1, ps(ks_top))
  t_bot_int = evaluation_polynomial( t_poly(kl,:), deg+1, ps(ks_bot))
  t_layer = t_mean(kl)

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]	nk	Number of levels
[in]	h	Layer thickness [H ~> m or kg m-2]
[in]	s	Layer salinity (conc, e.g. ppt)
[in]	c_method	Method to use for the centered difference
[in]	b_method	=1, use PCM in first/last cell, =2 uses linear extrapolation
[in,out]	diff	Scalar difference across layer (conc, e.g. ppt) determined by the following values for c_method: Second order finite difference (not recommended) Second order finite volume (used in original PPM) 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 809 of file MOM_neutral_diffusion.F90.

  integer,             intent(in)    :: nk       !< Number of levels
  real, dimension(nk), intent(in)    :: h        !< Layer thickness [H ~> m or kg m-2]
  real, dimension(nk), intent(in)    :: S        !< Layer salinity (conc, e.g. ppt)
  integer,             intent(in)    :: c_method !< Method to use for the centered difference
  integer,             intent(in)    :: b_method !< =1, use PCM in first/last cell, =2 uses linear extrapolation
  real, dimension(nk), intent(inout) :: diff     !< Scalar 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.

  ! Local variables
  integer :: k
  real :: hkm1, hk, hkp1, Skm1, Sk, Skp1, diff_l, diff_r, diff_c

  do k = 2, nk-1
    hkm1 = h(k-1)
    hk = h(k)
    hkp1 = h(k+1)

    if ( ( hkp1 + hk ) * ( hkm1 + hk ) > 0.) then
      skm1 = s(k-1)
      sk = s(k)
      skp1 = s(k+1)
      if (c_method==1) then
        ! Simple centered diff (from White)
        if ( hk + 0.5 * (hkm1 + hkp1) /= 0. ) then
          diff_c = ( skp1 - skm1 ) * ( hk / ( hk + 0.5 * (hkm1 + hkp1) ) )
        else
          diff_c = 0.
        endif
      elseif (c_method==2) then
        ! Second order accurate centered FV slope (from Colella and Woodward, JCP 1984)
        diff_c = fv_diff(hkm1, hk, hkp1, skm1, sk, skp1)
      elseif (c_method==3) then
        ! Second order accurate finite-volume least squares slope
        diff_c = hk * fvlsq_slope(hkm1, hk, hkp1, skm1, sk, skp1)
      endif
      ! Limit centered slope by twice the side differenced slopes
      diff_l = 2. * ( sk - skm1 )
      diff_r = 2. * ( skp1 - sk )
      if ( signum(1., diff_l) * signum(1., diff_r) <= 0. ) then
        diff(k) = 0. ! PCM for local extrema
      else
        diff(k) = sign( min( abs(diff_l), abs(diff_c), abs(diff_r) ), diff_c )
      endif
    else
      diff(k) = 0. ! PCM next to vanished layers
    endif
  enddo
  if (b_method==1) then ! PCM for top and bottom layer
    diff(1) = 0.
    diff(nk) = 0.
  elseif (b_method==2) then ! Linear extrapolation for top and bottom interfaces
    diff(1) = ( s(2) - s(1) ) * 2. * ( h(1) / ( h(1) + h(2) ) )
    diff(nk) = s(nk) - s(nk-1) * 2. * ( h(nk) / ( h(nk-1) + h(nk) ) )
  endif

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]	xl	Fraction position of left bound (0,1)
[in]	xr	Fraction position of right bound (0,1)
[in]	al	Left edge scalar value, at x=0
[in]	ar	Right edge scalar value, at x=1
[in]	amean	Average scalar value of cell

Definition at line 771 of file MOM_neutral_diffusion.F90.

  real, intent(in) :: xL    !< Fraction position of left bound (0,1)
  real, intent(in) :: xR    !< Fraction position of right bound (0,1)
  real, intent(in) :: aL    !< Left edge scalar value, at x=0
  real, intent(in) :: aR    !< Right edge scalar value, at x=1
  real, intent(in) :: aMean !< Average scalar value of cell

  ! Local variables
  real :: dx, xave, a6, a6o3

  dx = xr - xl
  xave = 0.5 * ( xr + xl )
  a6o3 = 2. * amean - ( al + ar ) ! a6 / 3.
  a6 = 3. * a6o3

  if (dx<0.) then
    stop 'ppm_ave: dx<0 should not happend!'
  elseif (dx>1.) then
    stop 'ppm_ave: dx>1 should not happend!'
  elseif (dx==0.) then
    ppm_ave = al + ( ar - al ) * xr + a6 * xr * ( 1. - xr )
  else
    ppm_ave = ( al + xave * ( ( ar - al ) + a6 ) )  - a6o3 * ( xr**2 + xr * xl + xl**2 )
  endif

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]	hkm1	Width of cell k-1
[in]	hk	Width of cell k
[in]	hkp1	Width of cell k+1
[in]	hkp2	Width of cell k+2
[in]	ak	Average scalar value of cell k
[in]	akp1	Average scalar value of cell k+1
[in]	pk	PLM slope for cell k
[in]	pkp1	PLM slope for cell k+1
[in]	h_neglect	A negligibly small thickness [H ~> m or kg m-2]

Definition at line 731 of file MOM_neutral_diffusion.F90.

  real, intent(in) :: hkm1 !< Width of cell k-1
  real, intent(in) :: hk   !< Width of cell k
  real, intent(in) :: hkp1 !< Width of cell k+1
  real, intent(in) :: hkp2 !< Width of cell k+2
  real, intent(in) :: Ak   !< Average scalar value of cell k
  real, intent(in) :: Akp1 !< Average scalar value of cell k+1
  real, intent(in) :: Pk   !< PLM slope for cell k
  real, intent(in) :: Pkp1 !< PLM slope for cell k+1
  real, intent(in) :: h_neglect !< A negligibly small thickness [H ~> m or kg m-2]

  ! Local variables
  real :: R_hk_hkp1, R_2hk_hkp1, R_hk_2hkp1, f1, f2, f3, f4

  r_hk_hkp1 = hk + hkp1
  if (r_hk_hkp1 <= 0.) then
    ppm_edge = 0.5 * ( ak + akp1 )
    return
  endif
  r_hk_hkp1 = 1. / r_hk_hkp1
  if (hk<hkp1) then
    ppm_edge = ak + ( hk * r_hk_hkp1 ) * ( akp1 - ak )
  else
    ppm_edge = akp1 + ( hkp1 * r_hk_hkp1 ) * ( ak - akp1 )
  endif

  r_2hk_hkp1 = 1. / ( ( 2. * hk + hkp1 ) + h_neglect )
  r_hk_2hkp1 = 1. / ( ( hk + 2. * hkp1 ) + h_neglect )
  f1 = 1./ ( ( hk + hkp1) + ( hkm1 + hkp2 ) )
  f2 = 2. * ( hkp1 * hk ) * r_hk_hkp1 * &
            ( ( hkm1 + hk ) * r_2hk_hkp1  - ( hkp2 + hkp1 ) * r_hk_2hkp1 )
  f3 = hk * ( hkm1 + hk ) * r_2hk_hkp1
  f4 = hkp1 * ( hkp1 + hkp2 ) * r_hk_2hkp1

  ppm_edge = ppm_edge + f1 * ( f2 * ( akp1 - ak ) - ( f3 * pkp1 - f4 * pk ) )

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]	nk	Number of levels
[in]	tl	Layer tracer (conc, e.g. degC)
[in]	ti	Interface tracer (conc, e.g. degC)
[in,out]	al	Left edge value of tracer (conc, e.g. degC)
[in,out]	ar	Right edge value of tracer (conc, e.g. degC)

Definition at line 2056 of file MOM_neutral_diffusion.F90.

  integer,                    intent(in)    :: nk !< Number of levels
  real, dimension(nk),        intent(in)    :: Tl !< Layer tracer (conc, e.g. degC)
  real, dimension(nk+1),      intent(in)    :: Ti !< Interface tracer (conc, e.g. degC)
  real, dimension(nk),        intent(inout) :: aL !< Left edge value of tracer (conc, e.g. degC)
  real, dimension(nk),        intent(inout) :: aR !< Right edge value of tracer (conc, e.g. degC)

  integer :: k
  ! Setup reconstruction edge values
  do k = 1, nk
    al(k) = ti(k)
    ar(k) = ti(k+1)
    if ( signum(1., ar(k) - tl(k))*signum(1., tl(k) - al(k)) <= 0.0 ) then
      al(k) = tl(k)
      ar(k) = tl(k)
    elseif ( sign(3., ar(k) - al(k)) * ( (tl(k) - al(k)) + (tl(k) - ar(k))) > abs(ar(k) - al(k)) ) then
      al(k) = tl(k) + 2.0 * ( tl(k) - ar(k) )
    elseif ( sign(3., ar(k) - al(k)) * ( (tl(k) - al(k)) + (tl(k) - ar(k))) < -abs(ar(k) - al(k)) ) then
      ar(k) = tl(k) + 2.0 * ( tl(k) - al(k) )
    endif
  enddo

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]	cs	Neutral diffusion control structure
[in]	ksurf	Current index of neutral surface
[in]	pos_last	Last position within the current layer, used as the lower bound in the root finding algorithm [nondim]
[in]	t_from	Temperature at the searched from interface [degC]
[in]	s_from	Salinity at the searched from interface [ppt]
[in]	p_from	Pressure at the searched from interface [R L2 T-2 ~> Pa]
[in]	t_top	Temperature at the searched to top interface [degC]
[in]	s_top	Salinity at the searched to top interface [ppt]
[in]	p_top	Pressure at the searched to top interface [R L2 T-2 ~> Pa] interface [R L2 T-2 ~> Pa]
[in]	t_bot	Temperature at the searched to bottom interface [degC]
[in]	s_bot	Salinity at the searched to bottom interface [ppt]
[in]	p_bot	Pressure at the searched to bottom interface [R L2 T-2 ~> Pa]
[in]	t_poly	Temperature polynomial reconstruction coefficients [degC]
[in]	s_poly	Salinity polynomial reconstruction coefficients [ppt]

Definition at line 1439 of file MOM_neutral_diffusion.F90.

  type(neutral_diffusion_CS), intent(in   ) :: CS       !< Neutral diffusion control structure
  integer,                    intent(in   ) :: ksurf    !< Current index of neutral surface
  real,                       intent(in   ) :: pos_last !< Last position within the current layer, used as the lower
                                                        !! bound in the root finding algorithm [nondim]
  real,                       intent(in   ) :: T_from   !< Temperature at the searched from interface [degC]
  real,                       intent(in   ) :: S_from   !< Salinity    at the searched from interface [ppt]
  real,                       intent(in   ) :: P_from   !< Pressure at the searched from interface [R L2 T-2 ~> Pa]
  real,                       intent(in   ) :: T_top    !< Temperature at the searched to top interface [degC]
  real,                       intent(in   ) :: S_top    !< Salinity    at the searched to top interface [ppt]
  real,                       intent(in   ) :: P_top    !< Pressure at the searched to top interface [R L2 T-2 ~> Pa]
                                                        !! interface [R L2 T-2 ~> Pa]
  real,                       intent(in   ) :: T_bot    !< Temperature at the searched to bottom interface [degC]
  real,                       intent(in   ) :: S_bot    !< Salinity    at the searched to bottom interface [ppt]
  real,                       intent(in   ) :: P_bot    !< Pressure at the searched to bottom
                                                        !! interface [R L2 T-2 ~> Pa]
  real, dimension(:),         intent(in   ) :: T_poly   !< Temperature polynomial reconstruction coefficients [degC]
  real, dimension(:),         intent(in   ) :: S_poly   !< Salinity    polynomial reconstruction coefficients [ppt]
  ! Local variables
  real :: dRhotop, dRhobot ! Density differences [R ~> kg m-3]
  real :: dRdT_top, dRdT_bot, dRdT_from ! Partial derivatives of density with temperature [R degC-1 ~> kg m-3 degC-1]
  real :: dRdS_top, dRdS_bot, dRdS_from ! Partial derivatives of density with salinity [R ppt-1 ~> kg m-3 ppt-1]

  ! Calculate the differencei in density at the tops or the bottom
  if (cs%neutral_pos_method == 1 .or. cs%neutral_pos_method == 3) then
    call calc_delta_rho_and_derivs(cs, t_top, s_top, p_top, t_from, s_from, p_from, drhotop)
    call calc_delta_rho_and_derivs(cs, t_bot, s_bot, p_bot, t_from, s_from, p_from, drhobot)
  elseif (cs%neutral_pos_method == 2) then
    call calc_delta_rho_and_derivs(cs, t_top, s_top, p_top, t_from, s_from, p_from, drhotop, &
                                   drdt_top, drds_top, drdt_from, drds_from)
    call calc_delta_rho_and_derivs(cs, t_bot, s_bot, p_bot, t_from, s_from, p_from, drhobot, &
                                   drdt_bot, drds_bot, drdt_from, drds_from)
  endif

  ! Handle all the special cases EXCEPT if it connects within the layer
  if ( (drhotop > 0.) .or. (ksurf == 1) ) then      ! First interface or lighter than anything in layer
    pos = pos_last
  elseif ( drhotop > drhobot ) then                 ! Unstably stratified
    pos = 1.
  elseif ( drhotop < 0. .and. drhobot < 0.) then    ! Denser than anything in layer
    pos = 1.
  elseif ( drhotop == 0. .and. drhobot == 0. ) then ! Perfectly unstratified
    pos = 1.
  elseif ( drhobot == 0. ) then                     ! Matches perfectly at the Top
    pos = 1.
  elseif ( drhotop == 0. ) then                     ! Matches perfectly at the Bottom
    pos = pos_last
  else                                              ! Neutral surface within layer
    pos = -1
  endif

  ! Can safely return if position is >= 0 otherwise will need to find the position within the layer
  if (pos>=0) return

  if (cs%neutral_pos_method==1) then
    pos = interpolate_for_nondim_position( drhotop, p_top, drhobot, p_bot )
  ! For the 'Linear' case of finding the neutral position, the fromerence pressure to use is the average
  ! of the midpoint of the layer being searched and the interface being searched from
  elseif (cs%neutral_pos_method == 2) then
    pos = find_neutral_pos_linear( cs, pos_last, t_from, s_from, drdt_from, drds_from, &
                                   drdt_top, drds_top, drdt_bot, drds_bot, t_poly, s_poly )
  elseif (cs%neutral_pos_method == 3) then
    pos = find_neutral_pos_full( cs, pos_last, t_from, s_from, p_from, p_top, p_bot, t_poly, s_poly)
  endif

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]	a	The magnitude argument
[in]	x	The sign (or zero) argument

Definition at line 798 of file MOM_neutral_diffusion.F90.

  real, intent(in) :: a !< The magnitude argument
  real, intent(in) :: x !< The sign (or zero) argument

  signum = sign(a,x)
  if (x==0.) signum = 0.

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]	verbose	If true, write results to stdout
[in]	nk	Number of layers
[in]	po	Calculated answer
[in]	ptrue	True answer
[in]	title	Title for messages

Definition at line 2706 of file MOM_neutral_diffusion.F90.

  logical,             intent(in) :: verbose !< If true, write results to stdout
  integer,             intent(in) :: nk    !< Number of layers
  real, dimension(nk), intent(in) :: Po    !< Calculated answer
  real, dimension(nk), intent(in) :: Ptrue !< True answer
  character(len=*),    intent(in) :: title !< Title for messages

  ! Local variables
  integer :: k, stdunit

  test_data1d = .false.
  do k = 1,nk
    if (po(k) /= ptrue(k)) test_data1d = .true.
  enddo

  if (test_data1d .or. verbose) then
    stdunit = stdout
    if (test_data1d) stdunit = stderr ! In case of wrong results, write to error stream
    write(stdunit,'(a)') title
    do k = 1,nk
      if (po(k) /= ptrue(k)) then
        test_data1d = .true.
        write(stdunit,'(a,i2,2(x,a,f20.16),x,a,1pe22.15,x,a)') &
              'k=',k,'Po=',po(k),'Ptrue=',ptrue(k),'err=',po(k)-ptrue(k),'WRONG!'
      else
        if (verbose) &
          write(stdunit,'(a,i2,2(x,a,f20.16),x,a,1pe22.15)') &
                'k=',k,'Po=',po(k),'Ptrue=',ptrue(k),'err=',po(k)-ptrue(k)
      endif
    enddo
  endif

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]	verbose	If true, write results to stdout
[in]	nk	Number of layers
[in]	po	Calculated answer
[in]	ptrue	True answer
[in]	title	Title for messages

Definition at line 2741 of file MOM_neutral_diffusion.F90.

  logical,                intent(in) :: verbose !< If true, write results to stdout
  integer,                intent(in) :: nk    !< Number of layers
  integer, dimension(nk), intent(in) :: Po    !< Calculated answer
  integer, dimension(nk), intent(in) :: Ptrue !< True answer
  character(len=*),       intent(in) :: title !< Title for messages

  ! Local variables
  integer :: k, stdunit

  test_data1di = .false.
  do k = 1,nk
    if (po(k) /= ptrue(k)) test_data1di = .true.
  enddo

  if (test_data1di .or. verbose) then
    stdunit = stdout
    if (test_data1di) stdunit = stderr ! In case of wrong results, write to error stream
    write(stdunit,'(a)') title
    do k = 1,nk
      if (po(k) /= ptrue(k)) then
        test_data1di = .true.
        write(stdunit,'(a,i2,2(x,a,i5),x,a)') 'k=',k,'Io=',po(k),'Itrue=',ptrue(k),'WRONG!'
      else
        if (verbose) &
          write(stdunit,'(a,i2,2(x,a,i5))') 'k=',k,'Io=',po(k),'Itrue=',ptrue(k)
      endif
    enddo
  endif

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]	verbose	If true, write results to stdout
[in]	hkm1	Left cell width [nondim]
[in]	hk	Center cell width [nondim]
[in]	hkp1	Right cell width [nondim]
[in]	skm1	Left cell average value
[in]	sk	Center cell average value
[in]	skp1	Right cell average value
[in]	ptrue	True answer [nondim]
[in]	title	Title for messages

Definition at line 2610 of file MOM_neutral_diffusion.F90.

  logical,          intent(in) :: verbose !< If true, write results to stdout
  real,             intent(in) :: hkm1  !< Left cell width [nondim]
  real,             intent(in) :: hk    !< Center cell width [nondim]
  real,             intent(in) :: hkp1  !< Right cell width [nondim]
  real,             intent(in) :: Skm1  !< Left cell average value
  real,             intent(in) :: Sk    !< Center cell average value
  real,             intent(in) :: Skp1  !< Right cell average value
  real,             intent(in) :: Ptrue !< True answer [nondim]
  character(len=*), intent(in) :: title !< Title for messages

  ! Local variables
  integer :: stdunit
  real :: Pret

  pret = fv_diff(hkm1, hk, hkp1, skm1, sk, skp1)
  test_fv_diff = (pret /= ptrue)

  if (test_fv_diff .or. verbose) then
    stdunit = stdout
    if (test_fv_diff) stdunit = stderr ! In case of wrong results, write to error stream
    write(stdunit,'(a)') title
    if (test_fv_diff) then
      write(stdunit,'(2(x,a,f20.16),x,a)') 'pRet=',pret,'pTrue=',ptrue,'WRONG!'
    else
      write(stdunit,'(2(x,a,f20.16))') 'pRet=',pret,'pTrue=',ptrue
    endif
  endif

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]	verbose	If true, write results to stdout
[in]	hkm1	Left cell width
[in]	hk	Center cell width
[in]	hkp1	Right cell width
[in]	skm1	Left cell average value
[in]	sk	Center cell average value
[in]	skp1	Right cell average value
[in]	ptrue	True answer
[in]	title	Title for messages

Definition at line 2642 of file MOM_neutral_diffusion.F90.

  logical,          intent(in) :: verbose !< If true, write results to stdout
  real,             intent(in) :: hkm1  !< Left cell width
  real,             intent(in) :: hk    !< Center cell width
  real,             intent(in) :: hkp1  !< Right cell width
  real,             intent(in) :: Skm1  !< Left cell average value
  real,             intent(in) :: Sk    !< Center cell average value
  real,             intent(in) :: Skp1  !< Right cell average value
  real,             intent(in) :: Ptrue !< True answer
  character(len=*), intent(in) :: title !< Title for messages

  ! Local variables
  integer :: stdunit
  real :: Pret

  pret = fvlsq_slope(hkm1, hk, hkp1, skm1, sk, skp1)
  test_fvlsq_slope = (pret /= ptrue)

  if (test_fvlsq_slope .or. verbose) then
    stdunit = stdout
    if (test_fvlsq_slope) stdunit = stderr ! In case of wrong results, write to error stream
    write(stdunit,'(a)') title
    if (test_fvlsq_slope) then
      write(stdunit,'(2(x,a,f20.16),x,a)') 'pRet=',pret,'pTrue=',ptrue,'WRONG!'
    else
      write(stdunit,'(2(x,a,f20.16))') 'pRet=',pret,'pTrue=',ptrue
    endif
  endif

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]	verbose	If true, write results to stdout
[in]	rhoneg	Lighter density [R ~> kg m-3]
[in]	pneg	Interface position of lighter density [nondim]
[in]	rhopos	Heavier density [R ~> kg m-3]
[in]	ppos	Interface position of heavier density [nondim]
[in]	ptrue	True answer [nondim]
[in]	title	Title for messages

Definition at line 2674 of file MOM_neutral_diffusion.F90.

  logical,          intent(in) :: verbose !< If true, write results to stdout
  real,             intent(in) :: rhoNeg !< Lighter density [R ~> kg m-3]
  real,             intent(in) :: Pneg   !< Interface position of lighter density [nondim]
  real,             intent(in) :: rhoPos !< Heavier density [R ~> kg m-3]
  real,             intent(in) :: Ppos   !< Interface position of heavier density [nondim]
  real,             intent(in) :: Ptrue  !< True answer [nondim]
  character(len=*), intent(in) :: title  !< Title for messages

  ! Local variables
  integer :: stdunit
  real :: Pret

  pret = interpolate_for_nondim_position(rhoneg, pneg, rhopos, ppos)
  test_ifndp = (pret /= ptrue)

  if (test_ifndp .or. verbose) then
    stdunit = stdout
    if (test_ifndp) stdunit = stderr ! In case of wrong results, write to error stream
    write(stdunit,'(a)') title
    if (test_ifndp) then
      write(stdunit,'(4(x,a,f20.16),2(x,a,1pe22.15),x,a)') &
            'r1=',rhoneg,'p1=',pneg,'r2=',rhopos,'p2=',ppos,'pRet=',pret,'pTrue=',ptrue,'WRONG!'
    else
      write(stdunit,'(4(x,a,f20.16),2(x,a,1pe22.15))') &
            'r1=',rhoneg,'p1=',pneg,'r2=',rhopos,'p2=',ppos,'pRet=',pret,'pTrue=',ptrue
    endif
  endif

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]	verbose	If true, write results to stdout
[in]	ns	Number of surfaces
[in]	kol	Index of first left interface above neutral surface
[in]	kor	Index of first right interface above neutral surface
[in]	pl	Fractional position of neutral surface within layer KoL of left column
[in]	pr	Fractional position of neutral surface within layer KoR of right column
[in]	heff	Effective thickness between two neutral surfaces [R L2 T-2 ~> Pa]
[in]	kol0	Correct value for KoL
[in]	kor0	Correct value for KoR
[in]	pl0	Correct value for pL
[in]	pr0	Correct value for pR
[in]	heff0	Correct value for hEff
[in]	title	Title for messages

Definition at line 2775 of file MOM_neutral_diffusion.F90.

  logical,                intent(in) :: verbose !< If true, write results to stdout
  integer,                intent(in) :: ns    !< Number of surfaces
  integer, dimension(ns), intent(in) :: KoL   !< Index of first left interface above neutral surface
  integer, dimension(ns), intent(in) :: KoR   !< Index of first right interface above neutral surface
  real, dimension(ns),    intent(in) :: pL    !< Fractional position of neutral surface within layer KoL of left column
  real, dimension(ns),    intent(in) :: pR    !< Fractional position of neutral surface within layer KoR of right column
  real, dimension(ns-1),  intent(in) :: hEff  !< Effective thickness between two neutral surfaces [R L2 T-2 ~> Pa]
  integer, dimension(ns), intent(in) :: KoL0  !< Correct value for KoL
  integer, dimension(ns), intent(in) :: KoR0  !< Correct value for KoR
  real, dimension(ns),    intent(in) :: pL0   !< Correct value for pL
  real, dimension(ns),    intent(in) :: pR0   !< Correct value for pR
  real, dimension(ns-1),  intent(in) :: hEff0 !< Correct value for hEff
  character(len=*),       intent(in) :: title !< Title for messages

  ! Local variables
  integer :: k, stdunit
  logical :: this_row_failed

  test_nsp = .false.
  do k = 1,ns
    test_nsp = test_nsp .or. compare_nsp_row(kol(k), kor(k), pl(k), pr(k), kol0(k), kor0(k), pl0(k), pr0(k))
    if (k < ns) then
      if (heff(k) /= heff0(k)) test_nsp = .true.
    endif
  enddo

  if (test_nsp .or. verbose) then
    stdunit = stdout
    if (test_nsp) stdunit = stderr ! In case of wrong results, write to error stream
    write(stdunit,'(a)') title
    do k = 1,ns
      this_row_failed = compare_nsp_row(kol(k), kor(k), pl(k), pr(k), kol0(k), kor0(k), pl0(k), pr0(k))
      if (this_row_failed) then
        write(stdunit,10) k,kol(k),pl(k),kor(k),pr(k),' <-- WRONG!'
        write(stdunit,10) k,kol0(k),pl0(k),kor0(k),pr0(k),' <-- should be this'
      else
        write(stdunit,10) k,kol(k),pl(k),kor(k),pr(k)
      endif
      if (k < ns) then
        if (heff(k) /= heff0(k)) then
          write(stdunit,'(i3,8x,"layer hEff =",2(f20.16,a))') k,heff(k)," .neq. ",heff0(k),' <-- WRONG!'
        else
          write(stdunit,'(i3,8x,"layer hEff =",f20.16)') k,heff(k)
        endif
      endif
    enddo
  endif
  if (test_nsp) call mom_error(fatal,"test_nsp failed")

10 format("ks=",i3," kL=",i3," pL=",f20.16," kR=",i3," pR=",f20.16,a)

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_pos	The expected position
[in]	test_pos	The position returned by the code
[in]	title	A label for this test

Definition at line 2847 of file MOM_neutral_diffusion.F90.

  real,             intent(in) :: expected_pos !< The expected position
  real,             intent(in) :: test_pos !< The position returned by the code
  character(len=*), intent(in) :: title    !< A label for this test
  ! Local variables
  integer :: stdunit

  stdunit = stdout ! Output to standard error
  test_rnp = abs(expected_pos - test_pos) > 2*epsilon(test_pos)
  if (test_rnp) then
    write(stdunit,'(A, f20.16, " .neq. ", f20.16, " <-- WRONG")') title, expected_pos, test_pos
  else
    write(stdunit,'(A, f20.16, " ==  ", f20.16)') title, expected_pos, test_pos
  endif