MOM_neutral_diffusion.F90

!> A column-wise toolbox for implementing neutral diffusion
module mom_neutral_diffusion
 
! This file is part of MOM6. See LICENSE.md for the license.
 
use mom_cpu_clock,             only : cpu_clock_id, cpu_clock_begin, cpu_clock_end
use mom_cpu_clock,             only : clock_module, clock_routine
use mom_diag_mediator,         only : diag_ctrl, time_type
use mom_diag_mediator,         only : post_data, register_diag_field
use mom_eos,                   only : eos_type, eos_manual_init, calculate_compress, calculate_density_derivs
use mom_eos,                   only : calculate_density_second_derivs
use mom_eos,                   only : extract_member_eos, eos_linear, eos_teos10, eos_wright
use mom_error_handler,         only : mom_error, fatal, warning, mom_mesg, is_root_pe
use mom_file_parser,           only : get_param, log_version, param_file_type
use mom_file_parser,           only : openparameterblock, closeparameterblock
use mom_grid,                  only : ocean_grid_type
use mom_neutral_diffusion_aux, only : ndiff_aux_cs_type, set_ndiff_aux_params
use mom_neutral_diffusion_aux, only : mark_unstable_cells, increment_interface, calc_drho, drho_at_pos
use mom_neutral_diffusion_aux, only : search_other_column, interpolate_for_nondim_position, refine_nondim_position
use mom_neutral_diffusion_aux, only : check_neutral_positions
use mom_remapping,             only : remapping_cs, initialize_remapping
use mom_remapping,             only : extract_member_remapping_cs, build_reconstructions_1d
use mom_remapping,             only : average_value_ppoly, remappingschemesdoc, remappingdefaultscheme
use mom_tracer_registry,       only : tracer_registry_type, tracer_type
use mom_verticalgrid,          only : verticalgrid_type
use polynomial_functions,      only : evaluation_polynomial, first_derivative_polynomial
use ppm_functions,             only : ppm_reconstruction, ppm_boundary_extrapolation
use regrid_edge_values,        only : edge_values_implicit_h4
 
implicit none ; private
 
#include <MOM_memory.h>
 
public neutral_diffusion, neutral_diffusion_init, neutral_diffusion_end
public neutral_diffusion_calc_coeffs
public neutral_diffusion_unit_tests
 
!> The control structure for the MOM_neutral_diffusion module
type, public :: neutral_diffusion_cs ; private
  integer :: nkp1   !< Number of interfaces for a column = nk + 1
  integer :: nsurf  !< Number of neutral surfaces
  integer :: deg = 2 !< Degree of polynomial used for reconstructions
  logical :: continuous_reconstruction = .true. !< True if using continuous PPM reconstruction at interfaces
  logical :: refine_position = .false. !< If true, iterate to refine the corresponding positions
                                       !! in neighboring columns
  logical :: debug = .false. !< If true, write verbose debugging messages
  integer :: max_iter !< Maximum number of iterations if refine_position is defined
  real :: tolerance   !< Convergence criterion representing difference from true neutrality
  real :: ref_pres    !< Reference pressure, negative if using locally referenced neutral density
 
  ! Positions of neutral surfaces in both the u, v directions
  real,    allocatable, dimension(:,:,:) :: upol  !< Non-dimensional position with left layer uKoL-1, u-point
  real,    allocatable, dimension(:,:,:) :: upor  !< Non-dimensional position with right layer uKoR-1, u-point
  integer, allocatable, dimension(:,:,:) :: ukol  !< Index of left interface corresponding to neutral surface,
                                                  !! at a u-point
  integer, allocatable, dimension(:,:,:) :: ukor  !< Index of right interface corresponding to neutral surface,
                                                  !! at a u-point
  real,    allocatable, dimension(:,:,:) :: uheff !< Effective thickness at u-point [H ~> m or kg m-2]
  real,    allocatable, dimension(:,:,:) :: vpol  !< Non-dimensional position with left layer uKoL-1, v-point
  real,    allocatable, dimension(:,:,:) :: vpor  !< Non-dimensional position with right layer uKoR-1, v-point
  integer, allocatable, dimension(:,:,:) :: vkol  !< Index of left interface corresponding to neutral surface,
                                                  !! at a v-point
  integer, allocatable, dimension(:,:,:) :: vkor  !< Index of right interface corresponding to neutral surface,
                                                  !! at a v-point
  real,    allocatable, dimension(:,:,:) :: vheff !< Effective thickness at v-point [H ~> m or kg m-2]
  ! Coefficients of polynomial reconstructions for temperature and salinity
  real,    allocatable, dimension(:,:,:,:) :: ppoly_coeffs_t !< Polynomial coefficients for temperature
  real,    allocatable, dimension(:,:,:,:) :: ppoly_coeffs_s !< Polynomial coefficients for salinity
  ! Variables needed for continuous reconstructions
  real,    allocatable, dimension(:,:,:) :: drdt !< dRho/dT [kg m-3 degC-1] at interfaces
  real,    allocatable, dimension(:,:,:) :: drds !< dRho/dS [kg m-3 ppt-1] at interfaces
  real,    allocatable, dimension(:,:,:) :: tint !< Interface T [degC]
  real,    allocatable, dimension(:,:,:) :: sint !< Interface S [ppt]
  real,    allocatable, dimension(:,:,:) :: pint !< Interface pressure [Pa]
  ! Variables needed for discontinuous reconstructions
  real,    allocatable, dimension(:,:,:,:) :: t_i    !< Top edge reconstruction of temperature [degC]
  real,    allocatable, dimension(:,:,:,:) :: s_i    !< Top edge reconstruction of salinity [ppt]
  real,    allocatable, dimension(:,:,:,:) :: drdt_i !< dRho/dT [kg m-3 degC-1] at top edge
  real,    allocatable, dimension(:,:,:,:) :: drds_i !< dRho/dS [kg m-3 ppt-1] at top edge
  integer, allocatable, dimension(:,:)     :: ns     !< Number of interfacs in a column
  logical, allocatable, dimension(:,:,:) :: stable_cell !< True if the cell is stably stratified wrt to the next cell
  type(diag_ctrl), pointer :: diag => null() !< A structure that is used to
                                             !! regulate the timing of diagnostic output.
  integer :: id_uheff_2d = -1 !< Diagnostic IDs
  integer :: id_vheff_2d = -1 !< Diagnostic IDs
 
  real    :: c_p !< heat capacity of seawater (J kg-1 K-1)
  type(eos_type), pointer :: eos !< Equation of state parameters
  type(remapping_cs)       :: remap_cs      !< Remapping control structure used to create sublayers
  type(ndiff_aux_cs_type), pointer :: ndiff_aux_cs !< Store parameters for iteratively finding neutral surface
end type neutral_diffusion_cs
 
! This include declares and sets the variable "version".
#include "version_variable.h"
character(len=40)  :: mdl = "MOM_neutral_diffusion" !< module name
 
contains
 
!> Read parameters and allocate control structure for neutral_diffusion module.
logical function neutral_diffusion_init(Time, G, param_file, diag, EOS, CS)
  type(time_type), target,    intent(in)    :: time       !< Time structure
  type(ocean_grid_type),      intent(in)    :: g          !< Grid structure
  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(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 :: boundary_extrap
  ! For refine_pos
  integer :: max_iter
  real :: drho_tol, xtol, ref_pres
 
  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 log_version(param_file, mdl, version, &
       "This module implements neutral diffusion of tracers")
  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) then
    return
  endif
 
  allocate(cs)
  allocate(cs%ndiff_aux_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.", &
                 default = -1.)
  ! Initialize and configure remapping
  if (cs%continuous_reconstruction .eqv. .false.) then
    call get_param(param_file, mdl, "NDIFF_BOUNDARY_EXTRAP", boundary_extrap, &
                   "Uses a rootfinding approach to find the position of a "//&
                   "neutral surface within a layer taking into account the "//&
                   "nonlinearity of the equation of state and the "//&
                   "polynomial reconstructions of T/S.",                         &
                   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 initialize_remapping( cs%remap_CS, string, boundary_extrapolation = boundary_extrap )
    call extract_member_remapping_cs(cs%remap_CS, degree=cs%deg)
    call get_param(param_file, mdl, "NDIFF_REFINE_POSITION", cs%refine_position, &
                   "Uses a rootfinding approach to find the position of a "//&
                   "neutral surface within a layer taking into account the "//&
                   "nonlinearity of the equation of state and the "//&
                   "polynomial reconstructions of T/S.",                         &
                   default=.false.)
    if (cs%refine_position) then
      call get_param(param_file, mdl, "NDIFF_DRHO_TOL", drho_tol,            &
                     "Sets the convergence criterion for finding the neutral "//&
                     "position within a layer in kg m-3.",                        &
                     default=1.e-10)
      call get_param(param_file, mdl, "NDIFF_X_TOL", xtol,            &
                     "Sets the convergence criterion for a change in nondim "//&
                     "position within a layer.",                        &
                     default=0.)
      call get_param(param_file, mdl, "NDIFF_MAX_ITER", max_iter,              &
                    "The maximum number of iterations to be done before "//&
                     "exiting the iterative loop to find the neutral surface",    &
                     default=10)
      call set_ndiff_aux_params(cs%ndiff_aux_CS, max_iter = max_iter, drho_tol = drho_tol, xtol = xtol)
    endif
    call get_param(param_file, mdl, "NDIFF_DEBUG", cs%debug,             &
                   "Turns on verbose output for discontinuous neutral "//&
                   "diffusion routines.", &
                   default = .false.)
    call set_ndiff_aux_params(cs%ndiff_aux_CS, deg=cs%deg, ref_pres = cs%ref_pres, eos = eos, debug = cs%debug)
  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%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
 
end function neutral_diffusion_init
 
!> Calculate remapping factors for u/v columns used to map adjoining columns to
!! a shared coordinate space.
subroutine neutral_diffusion_calc_coeffs(G, GV, h, T, S, CS)
  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(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
 
  ! Local variables
  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
  integer :: imethod
  real, dimension(SZI_(G)) :: ref_pres ! Reference pressure used to calculate alpha/beta
  real :: h_neglect, h_neglect_edge
  real :: pa_to_h
 
  pa_to_h = 1. / gv%H_to_pa
 
  !### Try replacing both of these with GV%H_subroundoff
  if (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
  cs%Pint(:,:,1) = 0.
  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%H_to_Pa
  enddo ; enddo ; enddo
 
  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 )
      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), g%isc-1, g%iec-g%isc+3, cs%EOS)
      enddo
    else ! Discontinuous reconstruction
      do k = 1, g%ke
        if (cs%ref_pres<0) ref_pres(:) = cs%Pint(:,j,k)
        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), g%isc-1, g%iec-g%isc+3, cs%EOS)
        if (cs%ref_pres<0) then
          ref_pres(:) = cs%Pint(:,j,k+1)
        endif
        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), g%isc-1, g%iec-g%isc+3, cs%EOS)
      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( g%ke, cs%dRdT_i(i,j,:,:), cs%dRdS_i(i,j,:,:), cs%T_i(i,j,:,:), cs%S_i(i,j,:,:), &
                                cs%stable_cell(i,j,:), cs%ns(i,j) )
    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,:) )
      else
        call find_neutral_surface_positions_discontinuous(cs, g%ke, cs%ns(i,j)+cs%ns(i+1,j), &
            cs%Pint(i,j,:), h(i,j,:), cs%T_i(i,j,:,:), cs%S_i(i,j,:,:),                      &
            cs%dRdT_i(i,j,:,:), cs%dRdS_i(i,j,:,:), cs%stable_cell(i,j,:),                   &
            cs%Pint(i+1,j,:), h(i+1,j,:), cs%T_i(i+1,j,:,:), cs%S_i(i+1,j,:,:),              &
            cs%dRdT_i(i+1,j,:,:), cs%dRdS_i(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,:), &
            cs%ppoly_coeffs_T(i,j,:,:), cs%ppoly_coeffs_S(i,j,:,:),                          &
            cs%ppoly_coeffs_T(i+1,j,:,:), cs%ppoly_coeffs_S(i+1,j,:,:))
      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,:) )
      else
        call find_neutral_surface_positions_discontinuous(cs, g%ke, cs%ns(i,j)+cs%ns(i,j+1), &
            cs%Pint(i,j,:), h(i,j,:), cs%T_i(i,j,:,:), cs%S_i(i,j,:,:),                      &
            cs%dRdT_i(i,j,:,:), cs%dRdS_i(i,j,:,:), cs%stable_cell(i,j,:),                   &
            cs%Pint(i,j+1,:), h(i,j+1,:), cs%T_i(i,j+1,:,:), cs%S_i(i,j+1,:,:),              &
            cs%dRdT_i(i,j+1,:,:), cs%dRdS_i(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,:), &
            cs%ppoly_coeffs_T(i,j,:,:), cs%ppoly_coeffs_S(i,j,:,:),                          &
            cs%ppoly_coeffs_T(i,j+1,:,:), cs%ppoly_coeffs_S(i,j+1,:,:))
 
      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 fraction of the nondimensional fraction of the layer occupied by
  ! the... (Please finish this thought. -RWH)
  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
 
end subroutine neutral_diffusion_calc_coeffs
 
!> Update tracer concentration due to neutral diffusion; layer thickness unchanged by this update.
subroutine neutral_diffusion(G, GV, h, Coef_x, Coef_y, dt, Reg, CS)
  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 [m2]
  real, dimension(SZI_(G),SZJB_(G)),         intent(in)    :: coef_y !< dt * Kh * dx / dy at v-points [m2]
  real,                                      intent(in)    :: dt     !< Tracer time step * I_numitts
                                                                     !! (I_numitts in tracer_hordiff)
  type(tracer_registry_type),                pointer       :: reg    !< Tracer registry
  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
  real :: h_neglect, h_neglect_edge
 
  !### Try replacing both of these with GV%H_subroundoff
  h_neglect_edge = gv%m_to_H*1.0e-10
  h_neglect = gv%m_to_H*1.0e-30
 
  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
 
end subroutine neutral_diffusion
 
!> Returns interface scalar, Si, for a column of layer values, S.
subroutine interface_scalar(nk, h, S, Si, i_method, h_neglect)
  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)
 
end subroutine interface_scalar
 
!> Returns the PPM quasi-fourth order edge value at k+1/2 following
!! equation 1.6 in Colella & Woodward, 1984: JCP 54, 174-201.
real function ppm_edge(hkm1, hk, hkp1, hkp2,  Ak, Akp1, Pk, Pkp1, h_neglect)
  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 ) )
 
end function ppm_edge
 
!> Returns the average of a PPM reconstruction between two
!! fractional positions.
real function ppm_ave(xL, xR, aL, aR, aMean)
  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
end function ppm_ave
 
!> A true signum function that returns either -abs(a), when x<0; or abs(a) when x>0; or 0 when x=0.
real function signum(a,x)
  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.
 
end function signum
 
!> 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.
subroutine plm_diff(nk, h, S, c_method, b_method, diff)
  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
 
end subroutine plm_diff
 
!> 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.
real function fv_diff(hkm1, hk, hkp1, Skm1, Sk, Skp1)
  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 ) )
end function fv_diff
 
 
!> 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).
real function fvlsq_slope(hkm1, hk, hkp1, Skm1, Sk, Skp1)
  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
end function fvlsq_slope
 
 
!> Returns positions within left/right columns of combined interfaces using continuous reconstructions of T/S
subroutine find_neutral_surface_positions_continuous(nk, Pl, Tl, Sl, dRdTl, dRdSl, Pr, Tr, Sr, &
                                                     dRdTr, dRdSr, PoL, PoR, KoL, KoR, hEff)
  integer,                    intent(in)    :: nk    !< Number of levels
  real, dimension(nk+1),      intent(in)    :: Pl    !< Left-column interface pressure [Pa]
  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 [kg m-3 degC-1]
  real, dimension(nk+1),      intent(in)    :: dRdSl !< Left-column dRho/dS [kg m-3 ppt-1]
  real, dimension(nk+1),      intent(in)    :: Pr    !< Right-column interface pressure [Pa]
  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 [kg m-3 degC-1]
  real, dimension(nk+1),      intent(in)    :: dRdSr !< Left-column dRho/dS [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 [Pa]
 
  ! 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 and dRho/dS 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, hL, hR
  integer :: lastK_left, lastK_right
  real    :: lastP_left, lastP_right
 
  ns = 2*nk+2
  ! Initialize variables for the search
  kr = 1 ; lastk_right = 1 ; lastp_right = 0.
  kl = 1 ; lastk_left = 1 ; lastp_left = 0.
  reached_bottom = .false.
 
  ! 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.
        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.
        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
 
    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
 
end subroutine find_neutral_surface_positions_continuous
 
!> Higher order version of find_neutral_surface_positions. Returns positions within left/right columns
!! of combined interfaces using intracell reconstructions of T/S
subroutine find_neutral_surface_positions_discontinuous(CS, nk, ns, Pres_l, hcol_l, Tl, Sl, &
                dRdT_l, dRdS_l, stable_l, Pres_r, hcol_r, Tr, Sr, dRdT_r, dRdS_r, stable_r, &
                PoL, PoR, KoL, KoR, hEff, ppoly_T_l, ppoly_S_l, ppoly_T_r, ppoly_S_r)
  type(neutral_diffusion_cs), intent(inout) :: CS  !< Neutral diffusion control structure
  integer,                    intent(in)    :: nk        !< Number of levels
  integer,                    intent(in)    :: ns        !< Number of neutral surfaces
  real, dimension(nk+1),      intent(in)    :: Pres_l    !< Left-column interface pressure [Pa]
  real, dimension(nk),        intent(in)    :: hcol_l    !< Left-column layer thicknesses
  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(nk,2),      intent(in)    :: dRdT_l    !< Left-column, top interface dRho/dT [kg m-3 degC-1]
  real, dimension(nk,2),      intent(in)    :: dRdS_l    !< Left-column, top interface dRho/dS [kg m-3 ppt-1]
  logical, dimension(nk),     intent(in)    :: stable_l  !< Left-column, top interface is stable
  real, dimension(nk+1),      intent(in)    :: Pres_r    !< Right-column interface pressure [Pa]
  real, dimension(nk),        intent(in)    :: hcol_r    !< Left-column layer thicknesses
  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(nk,2),      intent(in)    :: dRdT_r    !< Right-column, top interface dRho/dT [kg m-3 degC-1]
  real, dimension(nk,2),      intent(in)    :: dRdS_r    !< Right-column, top interface dRho/dS [kg m-3 ppt-1]
  logical, dimension(nk),     intent(in)    :: stable_r  !< Right-column, top interface is stable
  real, dimension(4*nk),      intent(inout) :: PoL       !< Fractional position of neutral surface within
                                                         !! layer KoL of left column
  real, dimension(4*nk),      intent(inout) :: PoR       !< Fractional position of neutral surface within
                                                         !! layer KoR of right column
  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 [Pa]
  real, dimension(nk,CS%deg+1), &
                    optional, intent(in)    :: ppoly_T_l !< Left-column coefficients of T reconstruction
  real, dimension(nk,CS%deg+1), &
                    optional, intent(in)    :: ppoly_S_l !< Left-column coefficients of S reconstruction
  real, dimension(nk,CS%deg+1), &
                    optional, intent(in)    :: ppoly_T_r !< Right-column coefficients of T reconstruction
  real, dimension(nk,CS%deg+1), &
                    optional, intent(in)    :: ppoly_S_r !< Right-column coefficients of S reconstruction
 
  ! Local variables
  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
  integer :: k, kl_left_0, kl_right_0
  real    :: dRho, dRhoTop, dRhoBot, hL, hR
  integer :: lastK_left, lastK_right
  real    :: lastP_left, lastP_right
  real    :: min_bound
  real    :: T_other, S_other, P_other, dRdT_other, dRdS_other
  logical, dimension(nk) :: top_connected_l, top_connected_r
  logical, dimension(nk) :: bot_connected_l, bot_connected_r
 
  top_connected_l(:) = .false. ; top_connected_r(:) = .false.
  bot_connected_l(:) = .false. ; bot_connected_r(:) = .false.
 
  ! Check to make sure that polynomial reconstructions were passed if refine_pos defined)
  if (cs%refine_position) then
    if (.not. ( present(ppoly_t_l) .and. present(ppoly_s_l) .and. &
                present(ppoly_t_r) .and. present(ppoly_s_r) ) ) &
        call mom_error(fatal, "fine_neutral_surface_positions_discontinuous: refine_pos is requested, but " //&
                              "polynomial coefficients not available for T and S")
  endif
  do k = 1,nk
    if (stable_l(k)) then
      kl_left = k
      kl_left_0 = k
      exit
    endif
  enddo
  do k = 1,nk
    if (stable_r(k)) then
      kl_right = k
      kl_right_0 = k
      exit
    endif
  enddo
 
  ! Initialize variables for the search
  ki_right = 1 ; lastk_right = 1 ; lastp_right = -1.
  ki_left = 1  ; lastk_left = 1  ; lastp_left = -1.
 
  reached_bottom = .false.
  searching_left_column = .false.
  searching_right_column = .false.
 
  ! Loop over each neutral surface, working from top to bottom
  neutral_surfaces: do k_surface = 1, ns
    ! Potential density difference, rho(kr) - rho(kl)
    drho = 0.5 * ( ( drdt_r(kl_right,ki_right) + drdt_l(kl_left,ki_left) ) * &
                   ( tr(kl_right,ki_right) - tl(kl_left,ki_left) ) &
                 + ( drds_r(kl_right,ki_right) + drds_l(kl_left,ki_left) ) * &
                   ( sr(kl_right,ki_right) - sl(kl_left,ki_left) ) )
    if (cs%debug)  write(*,'(A,I2,A,E12.4,A,I2,A,I2,A,I2,A,I2)') "k_surface=",k_surface,"  dRho=",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_right_column = .true.
        searching_left_column = .false.
      else ! dRho == 0.
        if (  ( kl_left == kl_left_0) .and. ( kl_right == kl_right_0 ) .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
      ! delta_rho is referenced to the right interface T, S, and P
      if (cs%ref_pres>=0.) then
        p_other = cs%ref_pres
      else
        if (ki_right == 1) p_other = pres_r(kl_right)
        if (ki_right == 2) p_other = pres_r(kl_right+1)
      endif
      t_other = tr(kl_right, ki_right)
      s_other = sr(kl_right, ki_right)
      drdt_other = drdt_r(kl_right, ki_right)
      drds_other = drds_r(kl_right, ki_right)
      if (cs%refine_position .and. (lastk_left == kl_left)) then
        call drho_at_pos(cs%ndiff_aux_CS, t_other, s_other, drdt_other, drds_other, pres_l(kl_left),       &
                         pres_l(kl_left+1), ppoly_t_l(kl_left,:), ppoly_s_l(kl_left,:), lastp_left, drhotop)
      else
        drhotop = calc_drho(tl(kl_left,1), sl(kl_left,1), drdt_l(kl_left,1), drds_l(kl_left,1), t_other, s_other, &
                            drdt_other, drds_other)
      endif
      ! Potential density difference, rho(kl) - rho(kl_right,ki_right) (will be positive)
      drhobot = calc_drho(tl(kl_left,2), sl(kl_left,2), drdt_l(kl_left,2), drds_l(kl_left,2), &
                          t_other, s_other, drdt_other, drds_other)
      if (cs%debug) then
        write(*,'(A,I2,A,E12.4,A,E12.4,A,E12.4)') "Searching left layer ", kl_left, &
                                                  " dRhoTop=", drhotop, "  dRhoBot=", drhobot
        write(*,'(A,I2,X,I2)') "Searching from right: ", kl_right, ki_right
        write(*,*) "Temp/Salt Reference: ", t_other, s_other
        write(*,*) "Temp/Salt Top L: ", tl(kl_left,1), sl(kl_left,1)
        write(*,*) "Temp/Salt Bot L: ", tl(kl_left,2), sl(kl_left,2)
      endif
 
      ! Set the position within the starting column
      por(k_surface) = real(ki_right-1)
      kor(k_surface) = kl_right
 
      ! Set position within the searched column
      call search_other_column(drhotop, drhobot, pres_l(kl_left), pres_l(kl_left+1), &
                               lastp_left, lastk_left, kl_left, kl_left_0, ki_left, &
                               top_connected_l, bot_connected_l, pol(k_surface), kol(k_surface), search_layer)
 
      if ( cs%refine_position .and. search_layer ) then
        min_bound = 0.
        if (k_surface > 1) then
          if ( kol(k_surface) == kol(k_surface-1) ) min_bound = pol(k_surface-1)
        endif
        pol(k_surface) = refine_nondim_position( cs%ndiff_aux_CS, t_other, s_other, drdt_other, drds_other, &
                            pres_l(kl_left), pres_l(kl_left+1), ppoly_t_l(kl_left,:), ppoly_s_l(kl_left,:), &
                            drhotop, drhobot, min_bound )
      endif
      if (pol(k_surface) == 0.) top_connected_l(kol(k_surface)) = .true.
      if (pol(k_surface) == 1.) bot_connected_l(kol(k_surface)) = .true.
      call increment_interface(nk, kl_right, ki_right, stable_r, reached_bottom, &
                               searching_right_column, searching_left_column)
 
    elseif (searching_right_column) then
      if (cs%ref_pres>=0.) then
        p_other = cs%ref_pres
      else
        if (ki_left == 1) p_other = pres_l(kl_left)
        if (ki_left == 2) p_other = pres_l(kl_left+1)
      endif
      t_other = tl(kl_left, ki_left)
      s_other = sl(kl_left, ki_left)
      drdt_other = drdt_l(kl_left, ki_left)
      drds_other = drds_l(kl_left, ki_left)
      ! Interpolate for the neutral surface position within the right column, layer krm1
      ! Potential density difference, rho(kr-1) - rho(kl) (should be negative)
 
      if (cs%refine_position .and. (lastk_right == kl_right)) then
        call drho_at_pos(cs%ndiff_aux_CS, t_other, s_other, drdt_other, drds_other, pres_r(kl_right),          &
                         pres_l(kl_right+1), ppoly_t_r(kl_right,:), ppoly_s_r(kl_right,:), lastp_right, drhotop)
      else
        drhotop = calc_drho(tr(kl_right,1), sr(kl_right,1), drdt_r(kl_right,1), drds_r(kl_right,1), &
                    t_other, s_other, drdt_other, drds_other)
      endif
      drhobot = calc_drho(tr(kl_right,2), sr(kl_right,2), drdt_r(kl_right,2), drds_r(kl_right,2), &
                  t_other, s_other, drdt_other, drds_other)
      if (cs%debug) then
        write(*,'(A,I2,A,E12.4,A,E12.4,A,E12.4)') "Searching right layer ", kl_right, &
                                                  "  dRhoTop=", drhotop, "  dRhoBot=", drhobot
        write(*,'(A,I2,X,I2)') "Searching from left: ", kl_left, ki_left
        write(*,*) "Temp/Salt Reference: ", t_other, s_other
        write(*,*) "Temp/Salt Top R: ", tr(kl_right,1), sr(kl_right,1)
        write(*,*) "Temp/Salt Bot R: ", tr(kl_right,2), sr(kl_right,2)
      endif
      ! Set the position within the starting column
      pol(k_surface) = real(ki_left-1)
      kol(k_surface) = kl_left
 
      ! Set position within the searched column
      call search_other_column(drhotop, drhobot, pres_r(kl_right), pres_r(kl_right+1), lastp_right, lastk_right,  &
                               kl_right, kl_right_0, ki_right, top_connected_r, bot_connected_r, por(k_surface),  &
                               kor(k_surface), search_layer)
      if ( cs%refine_position .and. search_layer) then
        min_bound = 0.
        if (k_surface > 1) then
          if ( kor(k_surface) == kor(k_surface-1) )  min_bound = por(k_surface-1)
        endif
        por(k_surface) = refine_nondim_position(cs%ndiff_aux_CS, t_other, s_other, drdt_other, drds_other,      &
                            pres_r(kl_right), pres_r(kl_right+1), ppoly_t_r(kl_right,:), ppoly_s_r(kl_right,:), &
                            drhotop, drhobot, min_bound )
      endif
      if (por(k_surface) == 0.) top_connected_r(kor(k_surface)) = .true.
      if (por(k_surface) == 1.) bot_connected_r(kor(k_surface)) = .true.
      call increment_interface(nk, kl_left, ki_left, stable_l, reached_bottom, &
                               searching_left_column, searching_right_column)
 
    else
      stop 'Else what?'
    endif
    lastk_left = kol(k_surface)  ; lastp_left = pol(k_surface)
    lastk_right = kor(k_surface) ; lastp_right = por(k_surface)
 
    if (cs%debug)  write(*,'(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)
    ! Effective thickness
    if (k_surface>1) then
      ! This is useful as a check to make sure that positions are monotonically increasing
      hl = absolute_position(nk,ns,pres_l,kol,pol,k_surface) - absolute_position(nk,ns,pres_l,kol,pol,k_surface-1)
      hr = absolute_position(nk,ns,pres_r,kor,por,k_surface) - absolute_position(nk,ns,pres_r,kor,por,k_surface-1)
      ! In the case of a layer being unstably stratified, may get a negative thickness. Set the previous position
      ! to the current location
      if ( hl<0. .or. hr<0. ) then
        heff(k_surface-1) = 0.
        call mom_error(warning, "hL or hR is negative")
      elseif ( hl > 0. .and. hr > 0.) 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))
        heff(k_surface-1) = 2. * ( (hl * hr) / ( hl + hr ) )! Harmonic mean
      else
        heff(k_surface-1) = 0.
      endif
    endif
  enddo neutral_surfaces
  if (cs%debug) then
    write (*,*) "==========Start Neutral Surfaces=========="
    do k = 1,ns-1
      if (heff(k)>0.) then
      kl_left = kol(k)
      kl_right = kor(k)
      write (*,'(A,I3,X,ES16.6,X,I3,X,ES16.6)') "Top surface KoL, PoL, KoR, PoR: ", kl_left, pol(k), kl_right, por(k)
      call check_neutral_positions(cs%ndiff_aux_CS, pres_l(kl_left), pres_l(kl_left+1), pres_r(kl_right), &
                                   pres_r(kl_right+1), pol(k), por(k), ppoly_t_l(kl_left,:), ppoly_t_r(kl_right,:), &
                                   ppoly_s_l(kl_left,:), ppoly_s_r(kl_right,:))
      kl_left = kol(k+1)
      kl_right = kor(k+1)
      write (*,'(A,I3,X,ES16.6,X,I3,X,ES16.6)') "Bot surface KoL, PoL, KoR, PoR: ", kl_left, pol(k+1), kl_right, por(k)
      call check_neutral_positions(cs%ndiff_aux_CS, pres_l(kl_left), pres_l(kl_left+1), pres_r(kl_right), &
                                   pres_r(kl_right+1), pol(k), por(k), ppoly_t_l(kl_left,:), ppoly_t_r(kl_right,:), &
                                   ppoly_s_l(kl_left,:), ppoly_s_r(kl_right,:))
      endif
    enddo
    write(*,'(A,E16.6)') "Total thickness of sublayers: ", sum(heff)
    write(*,*) "==========End Neutral Surfaces=========="
  endif
 
end subroutine find_neutral_surface_positions_discontinuous
 
!> Converts non-dimensional position within a layer to absolute position (for debugging)
real function absolute_position(n,ns,Pint,Karr,NParr,k_surface)
  integer, intent(in) :: n            !< Number of levels
  integer, intent(in) :: ns           !< Number of neutral surfaces
  real,    intent(in) :: pint(n+1)    !< Position of interfaces [Pa]
  integer, intent(in) :: karr(ns)     !< Index of interface above position
  real,    intent(in) :: nparr(ns)    !< Non-dimensional position within layer Karr(:)
  integer, intent(in) :: k_surface    !< k-interface to query
  ! 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) )
 
end function absolute_position
 
!> Converts non-dimensional positions within layers to absolute positions (for debugging)
function absolute_positions(n,ns,Pint,Karr,NParr)
  integer, intent(in) :: n            !< Number of levels
  integer, intent(in) :: ns           !< Number of neutral surfaces
  real,    intent(in) :: pint(n+1)    !< Position of interface [Pa]
  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 [Pa]
 
  ! 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
 
end function absolute_positions
 
!> Returns a single column of neutral diffusion fluxes of a tracer.
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)
  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 [Pa]
  real, dimension(nk),          intent(in)    :: hr    !< Right-column layer thickness [Pa]
  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 [Pa]
  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
                                             !! in the same units as h0.
  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
                                             !! for the purpose of edge value calculations
                                             !! in the same units as h0.
  ! 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
 
end subroutine neutral_surface_flux
 
!> Evaluate various parts of the reconstructions to calculate gradient-based flux limter
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)
  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)
 
end subroutine neutral_surface_t_eval
 
!> Discontinuous PPM reconstructions of the left/right edge values within a cell
subroutine ppm_left_right_edge_values(nk, Tl, Ti, aL, aR)
  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
end subroutine ppm_left_right_edge_values
 
!> Returns true if unit tests of neutral_diffusion functions fail. Otherwise returns false.
logical function neutral_diffusion_unit_tests(verbose)
  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)
 
 
end function neutral_diffusion_unit_tests
 
!> Returns true if unit tests of neutral_diffusion functions fail. Otherwise returns false.
logical function ndiff_unit_tests_continuous(verbose)
  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_edge
 
  h_neglect_edge = 1.0e-10 ; h_neglect = 1.0e-30
 
  v = verbose
 
  ndiff_unit_tests_continuous = .false. ! Normally return false
  write(*,*) '==== 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, h_neglect_edge=h_neglect_edge)
  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, h_neglect_edge=h_neglect_edge)
  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
                                   'Indentical 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(*,*) 'Pass'
 
end function ndiff_unit_tests_continuous
 
logical function ndiff_unit_tests_discontinuous(verbose)
  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, hl, hr
  real, dimension(nk,2)       :: til, sil, tir, sir
  real, dimension(nk+1)       :: pres_l, pres_r
  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)       :: poly_t_l, poly_t_r, poly_s,  poly_slope    ! Linear reconstruction for T
  real, dimension(nk,2)       :: drdt, drds
  logical, dimension(nk)      :: stable_l, stable_r
  integer                     :: imethod
  integer                     :: ns_l, ns_r
  real                        :: h_neglect, h_neglect_edge
  integer :: k
  logical :: v
 
  v = verbose
  ndiff_unit_tests_discontinuous = .false. ! Normally return false
  write(*,*) '==== MOM_neutral_diffusion: ndiff_unit_tests_discontinuous ='
 
  h_neglect = 1.0e-30 ; h_neglect_edge = 1.0e-10
 
  ! Unit tests for find_neutral_surface_positions_discontinuous
  ! Salinity is 0 for all these tests
  sl(:) = 0. ; sr(:) = 0. ; poly_s(:,:) = 0. ; sil(:,:) = 0. ; sir(:,:) = 0.
  drdt(:,:) = -1. ; drds(:,:) = 0.
 
  ! Intialize any control structures needed for unit tests
  cs%refine_position = .false.
  cs%ref_pres = -1.
  allocate(remap_cs)
  call initialize_remapping( remap_cs, "PLM", boundary_extrapolation = .true. )
 
  hl = (/10.,10.,10./) ; hr = (/10.,10.,10./) ; pres_l(1) = 0. ; pres_r(1) = 0.
  do k = 1,nk ; pres_l(k+1) = pres_l(k) + hl(k) ; pres_r(k+1) = pres_r(k) + hr(k) ; enddo
  ! Identical columns
  tl = (/20.,16.,12./) ; tr = (/20.,16.,12./)
  call build_reconstructions_1d( remap_cs, nk, hl, tl, poly_t_l, til, poly_slope, imethod, h_neglect, h_neglect_edge )
  call build_reconstructions_1d( remap_cs, nk, hr, tr, poly_t_r, tir, poly_slope, imethod, h_neglect, h_neglect_edge )
  call mark_unstable_cells( nk, drdt, drds, til, sil, stable_l, ns_l )
  call mark_unstable_cells( nk, drdt, drds, tir, sir, stable_r, ns_r )
  call find_neutral_surface_positions_discontinuous(cs, nk, ns_l+ns_r, pres_l, hl, til, sil, drdt, drds, stable_l, &
            pres_r, hr, tir, sir, drdt, drds, 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,2,3,3,3,3/), & ! KoR
                                   (/0.,0.,1.,1.,0.,0.,1.,1.,0.,0.,1.,1./), & ! pL
                                   (/0.,0.,1.,1.,0.,0.,1.,1.,0.,0.,1.,1./), & ! pR
                                   (/0.,10.,0.,0.,0.,10.,0.,0.,0.,10.,0.,0./), & ! hEff
                                   'Identical columns')
  tl = (/20.,16.,12./) ; tr = (/18.,14.,10./)
  call build_reconstructions_1d( remap_cs, nk, hl, tl, poly_t_l, til, poly_slope, imethod, h_neglect, h_neglect_edge )
  call build_reconstructions_1d( remap_cs, nk, hr, tr, poly_t_r, tir, poly_slope, imethod, h_neglect, h_neglect_edge )
  call mark_unstable_cells( nk, drdt, drds, til, sil, stable_l, ns_l )
  call mark_unstable_cells( nk, drdt, drds, tir, sir, stable_r, ns_r )
  call find_neutral_surface_positions_discontinuous(cs, nk, ns_l+ns_r, pres_l, hl, til, sil, drdt, drds, stable_l, &
            pres_r, hr, tir, sir, drdt, drds, 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.0, 0.5, 1.0, 0.0, 0.5, 0.5, 1.0, 0.0, 0.5, 0.5, 1.0, 1.0/), & ! pL
                                   (/0.0, 0.0, 0.5, 0.5, 1.0, 0.0, 0.5, 0.5, 1.0, 0.0, 0.5, 1.0/), & ! pR
                                   (/0.0, 5.0, 0.0, 5.0, 0.0, 5.0, 0.0, 5.0, 0.0, 5.0, 0.0/), & ! hEff
                                   'Right column slightly cooler')
  tl = (/18.,14.,10./) ; tr = (/20.,16.,12./)
  call build_reconstructions_1d( remap_cs, nk, hl, tl, poly_t_l, til, poly_slope, imethod, h_neglect, h_neglect_edge )
  call build_reconstructions_1d( remap_cs, nk, hr, tr, poly_t_r, tir, poly_slope, imethod, h_neglect, h_neglect_edge )
  call mark_unstable_cells( nk, drdt, drds, til, sil, stable_l, ns_l )
  call mark_unstable_cells( nk, drdt, drds, tir, sir, stable_r, ns_r )
  call find_neutral_surface_positions_discontinuous(cs, nk, ns_l+ns_r, pres_l, hl, til, sil, drdt, drds, stable_l, &
            pres_r, hr, tir, sir, drdt, drds, 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.0, 0.0, 0.5, 0.5, 1.0, 0.0, 0.5, 0.5, 1.0, 0.0, 0.5, 1.0/), & ! pL
                                   (/0.0, 0.5, 1.0, 0.0, 0.5, 0.5, 1.0, 0.0, 0.5, 0.5, 1.0, 1.0/), & ! pR
                                   (/0.0, 5.0, 0.0, 5.0, 0.0, 5.0, 0.0, 5.0, 0.0, 5.0, 0.0/), & ! hEff
                                   'Left column slightly cooler')
  tl = (/20.,16.,12./) ; tr = (/14.,10.,6./)
  call build_reconstructions_1d( remap_cs, nk, hl, tl, poly_t_l, til, poly_slope, imethod, h_neglect, h_neglect_edge )
  call build_reconstructions_1d( remap_cs, nk, hr, tr, poly_t_r, tir, poly_slope, imethod, h_neglect, h_neglect_edge )
  call mark_unstable_cells( nk, drdt, drds, til, sil, stable_l, ns_l )
  call mark_unstable_cells( nk, drdt, drds, tir, sir, stable_r, ns_r )
  call find_neutral_surface_positions_discontinuous(cs, nk, ns_l+ns_r, pres_l, hl, til, sil, drdt, drds, stable_l, &
            pres_r, hr, tir, sir, drdt, drds, 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,1,2,2,2,3,3/), & ! KoR
                                   (/0.0, 1.0, 0.0, 0.5, 1.0, 0.0, 0.5, 0.5, 1.0, 1.0, 1.0, 1.0/), & ! pL
                                   (/0.0, 0.0, 0.0, 0.0, 0.5, 0.5, 1.0, 0.0, 0.5, 1.0, 0.0, 1.0/), & ! pR
                                   (/0.0, 0.0, 0.0, 5.0, 0.0, 5.0, 0.0, 5.0, 0.0, 0.0, 0.0/), & ! hEff
                                   'Right column somewhat cooler')
  tl = (/20.,16.,12./) ; tr = (/8.,6.,4./)
  call build_reconstructions_1d( remap_cs, nk, hl, tl, poly_t_l, til, poly_slope, imethod, h_neglect, h_neglect_edge )
  call build_reconstructions_1d( remap_cs, nk, hr, tr, poly_t_r, tir, poly_slope, imethod, h_neglect, h_neglect_edge )
  call mark_unstable_cells( nk, drdt, drds, til, sil, stable_l, ns_l )
  call mark_unstable_cells( nk, drdt, drds, tir, sir, stable_r, ns_r )
  call find_neutral_surface_positions_discontinuous(cs, nk, ns_l+ns_r, pres_l, hl, til, sil, drdt, drds, stable_l, &
            pres_r, hr, tir, sir, drdt, drds, 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,3,3,3,3,3,3,3,3/),  & ! KoL
                                   (/1,1,1,1,1,1,1,1,2,2,3,3/), & ! KoR
                                   (/0.0, 1.0, 0.0, 1.0, 0.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0/), & ! pL
                                   (/0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0/), & ! pR
                                   (/0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0/), & ! hEff
                                   'Right column much cooler')
  tl = (/14.,14.,10./) ; tr = (/14.,14.,10./)
  call build_reconstructions_1d( remap_cs, nk, hl, tl, poly_t_l, til, poly_slope, imethod, h_neglect, h_neglect_edge )
  call build_reconstructions_1d( remap_cs, nk, hr, tr, poly_t_r, tir, poly_slope, imethod, h_neglect, h_neglect_edge )
  call mark_unstable_cells( nk, drdt, drds, til, sil, stable_l, ns_l )
  call mark_unstable_cells( nk, drdt, drds, tir, sir, stable_r, ns_r )
  call find_neutral_surface_positions_discontinuous(cs, nk, ns_l+ns_r, pres_l, hl, til, sil, drdt, drds, stable_l, &
            pres_r, hr, tir, sir, drdt, drds, 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,2,3,3,3,3/), & ! KoR
                                   (/0.,0.,1.,1.,0.,0.,1.,1.,0.,0.,1.,1./), & ! pL
                                   (/0.,0.,1.,1.,0.,0.,1.,1.,0.,0.,1.,1./), & ! pR
                                   (/0.,10.,0.,0.,0.,10.,0.,0.,0.,10.,0.,0./), & ! hEff
                                   'Identical columns with mixed layer')
  tl = (/20.,16.,12./) ; tr = (/14.,14.,10./)
  call build_reconstructions_1d( remap_cs, nk, hl, tl, poly_t_l, til, poly_slope, imethod, h_neglect, h_neglect_edge )
  call build_reconstructions_1d( remap_cs, nk, hr, tr, poly_t_r, tir, poly_slope, imethod, h_neglect, h_neglect_edge )
  call mark_unstable_cells( nk, drdt, drds, til, sil, stable_l, ns_l )
  call mark_unstable_cells( nk, drdt, drds, tir, sir, stable_r, ns_r )
  call find_neutral_surface_positions_discontinuous(cs, nk, ns_l+ns_r, pres_l, hl, til, sil, drdt, drds, stable_l, &
            pres_r, hr, tir, sir, drdt, drds, 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,2,3,3,3/), & ! KoR
                                   (/0.0, 1.0, 0.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.5, 1.0, 1.0/), & ! pL
                                   (/0.0, 0.0, 0.0, 0.0, 1.0, 1.0, 0.0, 0.0, 1.0, 0.0, 0.5, 1.0/), & ! pR
                                   (/0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 5.0, 0.0/), & ! hEff
                                   'Right column with mixed layer')
  tl = (/14.,14.,6./) ; tr = (/12.,16.,8./)
  call build_reconstructions_1d( remap_cs, nk, hl, tl, poly_t_l, til, poly_slope, imethod, h_neglect, h_neglect_edge )
  call build_reconstructions_1d( remap_cs, nk, hr, tr, poly_t_r, tir, poly_slope, imethod, h_neglect, h_neglect_edge )
  call mark_unstable_cells( nk, drdt, drds, til, sil, stable_l, ns_l )
  call mark_unstable_cells( nk, drdt, drds, tir, sir, stable_r, ns_r )
  call find_neutral_surface_positions_discontinuous(cs, nk, ns_l+ns_r, pres_l, hl, til, sil, drdt, drds, stable_l, &
            pres_r, hr, tir, sir, drdt, drds, stable_r, pol, por, kol, kor, heff)
  ndiff_unit_tests_discontinuous = ndiff_unit_tests_discontinuous .or.  test_nsp(v, 10, kol, kor, pol, por, heff, &
                                   (/1,1,1,1,2,2,2,3,3,3/), & ! KoL
                                   (/2,2,2,3,3,3,3,3,3,3/), & ! KoR
                                   (/0.0, 0.0, 0.0, 1.0, 0.0, 1.0, 1.0, 0.0, .75, 1.0/), & ! pL
                                   (/0.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0, .25, 1.0, 1.0/), & ! pR
                                   (/0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 7.5, 0.0/), & ! hEff
                                   'Left mixed layer, right unstable mixed layer')
 
  tl = (/10.,11.,6./) ; tr = (/12.,13.,8./)
  til(:,1) = (/8.,12.,10./) ; til(:,2) = (/12.,10.,2./)
  tir(:,1) = (/10.,14.,12./) ; tir(:,2) = (/14.,12.,4./)
  call mark_unstable_cells( nk, drdt, drds, til, sil, stable_l, ns_l )
  call mark_unstable_cells( nk, drdt, drds, tir, sir, stable_r, ns_r )
  call find_neutral_surface_positions_discontinuous(cs, nk, ns_l+ns_r, pres_l, hl, til, sil, drdt, drds, stable_l, &
            pres_r, hr, tir, sir, drdt, drds, stable_r, pol, por, kol, kor, heff)
  ndiff_unit_tests_discontinuous = ndiff_unit_tests_discontinuous .or.  test_nsp(v, 8, kol, kor, pol, por, heff, &
                                   (/2,2,2,2,2,3,3,3/), & ! KoL
                                   (/2,2,2,3,3,3,3,3/), & ! KoR
                                   (/0.0, 0.0, 0.0, 0.0, 1.0, 0.0, .75, 1.0/), & ! pL
                                   (/0.0, 1.0, 1.0, 0.0, .25, .25, 1.0, 1.0/), & ! pR
                                   (/0.0, 0.0, 0.0, 4.0, 0.0, 7.5, 0.0/), & ! hEff
                                   'Two unstable mixed layers')
  deallocate(remap_cs)
 
  allocate(eos)
  call eos_manual_init(eos, form_of_eos = eos_linear, drho_dt = -1., drho_ds = 2.)
  ! Unit tests for refine_nondim_position
  allocate(cs%ndiff_aux_CS)
  call set_ndiff_aux_params(cs%ndiff_aux_CS, deg = 1, max_iter = 10, drho_tol = 0., xtol = 0., eos = eos)
  ! Tests using Newton's method
  ndiff_unit_tests_discontinuous = ndiff_unit_tests_discontinuous .or. (test_rnp(0.5,refine_nondim_position( &
            cs%ndiff_aux_CS, 20., 35., -1., 2., 0., 1., (/21., -2./), (/35., 0./), -1., 1., 0.), &
            "Temperature stratified (Newton) "))
  ndiff_unit_tests_discontinuous = ndiff_unit_tests_discontinuous .or. (test_rnp(0.5,refine_nondim_position( &
            cs%ndiff_aux_CS, 20., 35., -1., 2., 0., 1., (/20., 0./), (/34., 2./),  -2., 2., 0.), &
            "Salinity stratified    (Newton) "))
  ndiff_unit_tests_discontinuous = ndiff_unit_tests_discontinuous .or. (test_rnp(0.5,refine_nondim_position( &
            cs%ndiff_aux_CS, 20., 35., -1., 2., 0., 1., (/21., -2./), (/34., 2./), -1., 1., 0.), &
            "Temp/Salt stratified   (Newton) "))
  call set_ndiff_aux_params(cs%ndiff_aux_CS, force_brent = .true.)
  ! Tests using Brent's method
  ndiff_unit_tests_discontinuous = ndiff_unit_tests_discontinuous .or. (test_rnp(0.5,refine_nondim_position( &
            cs%ndiff_aux_CS, 20., 35., -1., 2., 0., 1., (/21., -2./), (/35., 0./), -1., 1., 0.), &
            "Temperature stratified (Brent)  "))
  ndiff_unit_tests_discontinuous = ndiff_unit_tests_discontinuous .or. (test_rnp(0.5,refine_nondim_position( &
            cs%ndiff_aux_CS, 20., 35., -1., 2., 0., 1., (/20., 0./), (/34., 2./),  -2., 2., 0.), &
            "Salinity stratified    (Brent)  "))
  ndiff_unit_tests_discontinuous = ndiff_unit_tests_discontinuous .or. (test_rnp(0.5,refine_nondim_position( &
            cs%ndiff_aux_CS, 20., 35., -1., 2., 0., 1., (/21., -2./), (/34., 2./), -1., 1., 0.), &
            "Temp/Salt stratified   (Brent)  "))
  deallocate(eos)
 
  if (.not. ndiff_unit_tests_discontinuous) write(*,*) 'Pass'
 
end function ndiff_unit_tests_discontinuous
 
!> Returns true if a test of fv_diff() fails, and conditionally writes results to stream
logical function test_fv_diff(verbose, hkm1, hk, hkp1, Skm1, Sk, Skp1, Ptrue, title)
  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 [Pa]
  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 = 6
    if (test_fv_diff) stdunit = 0 ! 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
 
end function test_fv_diff
 
!> Returns true if a test of fvlsq_slope() fails, and conditionally writes results to stream
logical function test_fvlsq_slope(verbose, hkm1, hk, hkp1, Skm1, Sk, Skp1, Ptrue, title)
  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 [Pa]
  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 = 6
    if (test_fvlsq_slope) stdunit = 0 ! 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
 
end function test_fvlsq_slope
 
!> Returns true if a test of interpolate_for_nondim_position() fails, and conditionally writes results to stream
logical function test_ifndp(verbose, rhoNeg, Pneg, rhoPos, Ppos, Ptrue, title)
  logical,          intent(in) :: verbose !< If true, write results to stdout
  real,             intent(in) :: rhoneg !< Lighter density [kg m-3]
  real,             intent(in) :: pneg   !< Interface position of lighter density [Pa]
  real,             intent(in) :: rhopos !< Heavier density [kg m-3]
  real,             intent(in) :: ppos   !< Interface position of heavier density [Pa]
  real,             intent(in) :: ptrue  !< True answer [Pa]
  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 = 6
    if (test_ifndp) stdunit = 0 ! 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
 
end function test_ifndp
 
!> Returns true if comparison of Po and Ptrue fails, and conditionally writes results to stream
logical function test_data1d(verbose, nk, Po, Ptrue, title)
  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 = 6
    if (test_data1d) stdunit = 0 ! 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
 
end function test_data1d
 
!> Returns true if comparison of Po and Ptrue fails, and conditionally writes results to stream
logical function test_data1di(verbose, nk, Po, Ptrue, title)
  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 = 6
    if (test_data1di) stdunit = 0 ! 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
 
end function test_data1di
 
!> Returns true if output of find_neutral_surface_positions() does not match correct values,
!! and conditionally writes results to stream
logical function test_nsp(verbose, ns, KoL, KoR, pL, pR, hEff, KoL0, KoR0, pL0, pR0, hEff0, title)
  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 [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 = 6
    if (test_nsp) stdunit = 0 ! 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)
end function test_nsp
 
!> Compares a single row, k, of output from find_neutral_surface_positions()
logical function compare_nsp_row(KoL, KoR, pL, pR, KoL0, KoR0, pL0, pR0)
  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.
end function compare_nsp_row
 
!> Compares output position from refine_nondim_position with an expected value
logical function test_rnp(expected_pos, test_pos, title)
  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 = 6 ! Output to standard error
  test_rnp = expected_pos /= 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
end function test_rnp
!> Deallocates neutral_diffusion control structure
subroutine neutral_diffusion_end(CS)
  type(neutral_diffusion_cs), pointer :: cs  !< Neutral diffusion control structure
 
  if (associated(cs)) deallocate(cs)
 
end subroutine neutral_diffusion_end
 
end module mom_neutral_diffusion