MOM_lateral_boundary_diffusion.F90

!> Calculates and applies diffusive fluxes as a parameterization of lateral mixing (non-neutral) by
!! mesoscale eddies near the top and bottom (to be implemented) boundary layers of the ocean.

module mom_lateral_boundary_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_domains,               only : pass_var
use mom_diag_mediator,         only : diag_ctrl, time_type
use mom_diag_mediator,         only : post_data, register_diag_field
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_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_unit_scaling,          only : unit_scale_type
use mom_verticalgrid,          only : verticalgrid_type
use mom_cvmix_kpp,             only : kpp_get_bld, kpp_cs
use mom_energetic_pbl,         only : energetic_pbl_get_mld, energetic_pbl_cs
use mom_diabatic_driver,       only : diabatic_cs, extract_diabatic_member
use iso_fortran_env,           only : stdout=>output_unit, stderr=>error_unit

implicit none ; private

public near_boundary_unit_tests, lateral_boundary_diffusion, lateral_boundary_diffusion_init
public boundary_k_range

! Private parameters to avoid doing string comparisons for bottom or top boundary layer
integer, public, parameter :: surface = -1 !< Set a value that corresponds to the surface bopundary
integer, public, parameter :: bottom  = 1  !< Set a value that corresponds to the bottom boundary
#include <MOM_memory.h>

!> Sets parameters for lateral boundary mixing module.
type, public :: lateral_boundary_diffusion_cs ; private
  integer :: method                          !< Determine which of the three methods calculate
                                             !! and apply near boundary layer fluxes
                                             !! 1. Along layer
                                             !! 2. Bulk-layer approach (not recommended)
  integer :: deg                             !< Degree of polynomial reconstruction
  integer :: surface_boundary_scheme         !< Which boundary layer scheme to use
                                             !! 1. ePBL; 2. KPP
  logical :: limiter                         !< Controls wether a flux limiter is applied.
                                             !! Only valid when method = 2.
  logical :: linear                          !< If True, apply a linear transition at the base/top of the boundary.
                                             !! The flux will be fully applied at k=k_min and zero at k=k_max.

  type(remapping_cs)              :: remap_cs                     !< Control structure to hold remapping configuration
  type(kpp_cs),           pointer :: kpp_csp => null()            !< KPP control structure needed to get BLD
  type(energetic_pbl_cs), pointer :: energetic_pbl_csp => null()  !< ePBL control structure needed to get BLD
  type(diag_ctrl), pointer :: diag => null()                      !< A structure that is used to
                                                                  !! regulate the timing of diagnostic output.
end type lateral_boundary_diffusion_cs

! This include declares and sets the variable "version".
#include "version_variable.h"
character(len=40) :: mdl = "MOM_lateral_boundary_diffusion"       !< Name of this module

contains

!> Initialization routine that reads runtime parameters and sets up pointers to other control structures that might be
!! needed for lateral boundary diffusion.
logical function lateral_boundary_diffusion_init(Time, G, param_file, diag, diabatic_CSp, CS)
  type(time_type), target,          intent(in)    :: time          !< Time structure
  type(ocean_grid_type),            intent(in)    :: g             !< Grid structure
  type(param_file_type),            intent(in)    :: param_file    !< Parameter file structure
  type(diag_ctrl), target,          intent(inout) :: diag          !< Diagnostics control structure
  type(diabatic_cs),                pointer       :: diabatic_csp  !< KPP control structure needed to get BLD
  type(lateral_boundary_diffusion_cs), pointer    :: cs            !< Lateral boundary mixing control structure

  ! local variables
  character(len=80)  :: string  ! Temporary strings
  logical :: boundary_extrap

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

  ! Log this module and master switch for turning it on/off
  call get_param(param_file, mdl, "USE_LATERAL_BOUNDARY_DIFFUSION", lateral_boundary_diffusion_init, &
                 default=.false., do_not_log=.true.)
  call log_version(param_file, mdl, version, &
           "This module implements lateral diffusion of tracers near boundaries", &
           all_default=.not.lateral_boundary_diffusion_init)
  call get_param(param_file, mdl, "USE_LATERAL_BOUNDARY_DIFFUSION", lateral_boundary_diffusion_init, &
                 "If true, enables the lateral boundary tracer's diffusion module.", &
                 default=.false.)

  if (.not. lateral_boundary_diffusion_init) return

  allocate(cs)
  cs%diag => diag
  call extract_diabatic_member(diabatic_csp, kpp_csp=cs%KPP_CSp)
  call extract_diabatic_member(diabatic_csp, energetic_pbl_csp=cs%energetic_PBL_CSp)

  cs%surface_boundary_scheme = -1
  if ( .not. ASSOCIATED(cs%energetic_PBL_CSp) .and. .not. ASSOCIATED(cs%KPP_CSp) ) then
    call mom_error(fatal,"Lateral boundary diffusion is true, but no valid boundary layer scheme was found")
  endif

  ! Read all relevant parameters and write them to the model log.
  call get_param(param_file, mdl, "LATERAL_BOUNDARY_METHOD", cs%method, &
                 "Determine how to apply boundary lateral diffusion of tracers: \n"//&
                 "1. Along layer approach \n"//&
                 "2. Bulk layer approach (this option is not recommended)", default=1)
  if (cs%method == 2) then
    call get_param(param_file, mdl, "APPLY_LIMITER", cs%limiter, &
                   "If True, apply a flux limiter in the LBD. This is only available \n"//&
                   "when LATERAL_BOUNDARY_METHOD=2.", default=.false.)
  endif
  call get_param(param_file, mdl, "LBD_LINEAR_TRANSITION", cs%linear, &
                 "If True, apply a linear transition at the base/top of the boundary. \n"//&
                 "The flux will be fully applied at k=k_min and zero at k=k_max.", default=.false.)
  call get_param(param_file, mdl, "LBD_BOUNDARY_EXTRAP", boundary_extrap, &
                 "Use boundary extrapolation in LBD code", &
                 default=.false.)
  call get_param(param_file, mdl, "LBD_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)

end function lateral_boundary_diffusion_init

!> Driver routine for calculating lateral diffusive fluxes near the top and bottom boundaries.
!! Two different methods are available:
!! Method 1: lower order representation, calculate fluxes from bulk layer integrated quantities.
!! Method 2: more straight forward, diffusion is applied layer by layer using only information
!! from neighboring cells.
subroutine lateral_boundary_diffusion(G, GV, US, h, Coef_x, Coef_y, dt, Reg, CS)
  type(ocean_grid_type),                intent(inout) :: g   !< Grid type
  type(verticalgrid_type),              intent(in)    :: gv  !< ocean vertical grid structure
  type(unit_scale_type),                intent(in)    :: us  !< A dimensional unit scaling type
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
                                        intent(in)    :: h      !< Layer thickness [H ~> m or kg m-2]
  real, dimension(SZIB_(G),SZJ_(G)),    intent(in)    :: coef_x !< dt * Kh * dy / dx at u-points [L2 ~> m2]
  real, dimension(SZI_(G),SZJB_(G)),    intent(in)    :: coef_y !< dt * Kh * dx / dy at v-points [L2 ~> m2]
  real,                                 intent(in)    :: dt     !< Tracer time step * I_numitts
                                                                !! (I_numitts in tracer_hordiff)
  type(tracer_registry_type),           pointer       :: reg    !< Tracer registry
  type(lateral_boundary_diffusion_cs),  intent(in)    :: cs     !< Control structure for this module

  ! Local variables
  real, dimension(SZI_(G),SZJ_(G)) :: hbl                           !< bnd. layer depth [H ~> m or kg m-2]
  real, dimension(SZI_(G),SZJ_(G),SZK_(G),CS%deg+1) :: ppoly0_coefs !< Coefficients of polynomial
  real, dimension(SZI_(G),SZJ_(G),SZK_(G),2)        :: ppoly0_e     !< Edge values from reconstructions
  real, dimension(SZK_(G),CS%deg+1)                 :: ppoly_s      !< Slopes from reconstruction (placeholder)
  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)) :: uflx        !< Zonal flux of tracer [conc m^3]
  real, dimension(SZIB_(G),SZJ_(G))         :: uflx_bulk   !< Total calculated bulk-layer u-flux for the tracer
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)) :: vflx        !< Meridional flux of tracer [conc m^3]
  real, dimension(SZI_(G),SZJB_(G))         :: vflx_bulk   !< Total calculated bulk-layer v-flux for the tracer
  real, dimension(SZIB_(G),SZJ_(G))         :: uwork_2d    !< Layer summed u-flux transport
  real, dimension(SZI_(G),SZJB_(G))         :: vwork_2d    !< Layer summed v-flux transport
  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
  type(tracer_type), pointer                :: tracer => null() !< Pointer to the current tracer
  integer :: remap_method !< Reconstruction method
  integer :: i,j,k,m      !< indices to loop over
  real    :: idt          !< inverse of the time step [s-1]

  idt = 1./dt
  hbl(:,:) = 0.
  if (ASSOCIATED(cs%KPP_CSp)) call kpp_get_bld(cs%KPP_CSp, hbl, g, us, m_to_bld_units=gv%m_to_H)
  if (ASSOCIATED(cs%energetic_PBL_CSp)) &
    call energetic_pbl_get_mld(cs%energetic_PBL_CSp, hbl, g, us, m_to_mld_units=gv%m_to_H)

  call pass_var(hbl,g%Domain)
  do m = 1,reg%ntr
    tracer => reg%tr(m)

    ! for diagnostics
    if (tracer%id_lbdxy_conc > 0 .or. tracer%id_lbdxy_cont > 0 .or. tracer%id_lbdxy_cont_2d > 0) then
      tendency(:,:,:) = 0.0
    endif

    ! Interpolate state to interface
    do j=g%jsc-1,g%jec+1 ; do i=g%isc-1,g%iec+1
      call build_reconstructions_1d( cs%remap_CS, g%ke, h(i,j,:), tracer%t(i,j,:), ppoly0_coefs(i,j,:,:), &
                                     ppoly0_e(i,j,:,:), ppoly_s, remap_method, gv%H_subroundoff, gv%H_subroundoff)
    enddo ; enddo

    ! Diffusive fluxes in the i-direction
    uflx(:,:,:) = 0.
    vflx(:,:,:) = 0.
    uflx_bulk(:,:) = 0.
    vflx_bulk(:,:) = 0.

    ! Method #1 (layer by layer)
    if (cs%method == 1) then
      do j=g%jsc,g%jec
        do i=g%isc-1,g%iec
          if (g%mask2dCu(i,j)>0.) then
            call fluxes_layer_method(surface, gv%ke, cs%deg, h(i,j,:), h(i+1,j,:), hbl(i,j), hbl(i+1,j),  &
              g%areaT(i,j), g%areaT(i+1,j), tracer%t(i,j,:), tracer%t(i+1,j,:), ppoly0_coefs(i,j,:,:),    &
              ppoly0_coefs(i+1,j,:,:), ppoly0_e(i,j,:,:), ppoly0_e(i+1,j,:,:), remap_method, coef_x(i,j), &
              uflx(i,j,:), cs%linear)
          endif
        enddo
      enddo
      do j=g%jsc-1,g%jec
        do i=g%isc,g%iec
          if (g%mask2dCv(i,j)>0.) then
            call fluxes_layer_method(surface, gv%ke, cs%deg, h(i,j,:), h(i,j+1,:), hbl(i,j), hbl(i,j+1),  &
              g%areaT(i,j), g%areaT(i,j+1), tracer%t(i,j,:), tracer%t(i,j+1,:), ppoly0_coefs(i,j,:,:),    &
              ppoly0_coefs(i,j+1,:,:), ppoly0_e(i,j,:,:), ppoly0_e(i,j+1,:,:), remap_method, coef_y(i,j), &
              vflx(i,j,:), cs%linear)
          endif
        enddo
      enddo

    ! Method #2 (bulk approach)
    elseif (cs%method == 2) then
      do j=g%jsc,g%jec
        do i=g%isc-1,g%iec
          if (g%mask2dCu(i,j)>0.) then
            call fluxes_bulk_method(surface, gv%ke, cs%deg, h(i,j,:), h(i+1,j,:), hbl(i,j), hbl(i+1,j), &
              g%areaT(i,j), g%areaT(i+1,j), tracer%t(i,j,:), tracer%t(i+1,j,:),                         &
              ppoly0_coefs(i,j,:,:), ppoly0_coefs(i+1,j,:,:), ppoly0_e(i,j,:,:),                        &
              ppoly0_e(i+1,j,:,:), remap_method, coef_x(i,j), uflx_bulk(i,j), uflx(i,j,:), cs%limiter,  &
              cs%linear)
          endif
        enddo
      enddo
      do j=g%jsc-1,g%jec
        do i=g%isc,g%iec
          if (g%mask2dCv(i,j)>0.) then
            call fluxes_bulk_method(surface, gv%ke, cs%deg, h(i,j,:), h(i,j+1,:), hbl(i,j), hbl(i,j+1), &
              g%areaT(i,j), g%areaT(i,j+1), tracer%t(i,j,:), tracer%t(i,j+1,:),                         &
              ppoly0_coefs(i,j,:,:), ppoly0_coefs(i,j+1,:,:), ppoly0_e(i,j,:,:),                        &
              ppoly0_e(i,j+1,:,:), remap_method, coef_y(i,j), vflx_bulk(i,j), vflx(i,j,:), cs%limiter,  &
              cs%linear)
          endif
        enddo
      enddo
      ! Post tracer bulk diags
      if (tracer%id_lbd_bulk_dfx>0) call post_data(tracer%id_lbd_bulk_dfx, uflx_bulk*idt, cs%diag)
      if (tracer%id_lbd_bulk_dfy>0) call post_data(tracer%id_lbd_bulk_dfy, vflx_bulk*idt, cs%diag)
    endif

    ! Update the tracer fluxes
    do k=1,gv%ke ; do j=g%jsc,g%jec ; do i=g%isc,g%iec
      if (g%mask2dT(i,j)>0.) then
        tracer%t(i,j,k) = tracer%t(i,j,k) + (( (uflx(i-1,j,k)-uflx(i,j,k)) ) + ( (vflx(i,j-1,k)-vflx(i,j,k) ) ))* &
                          (g%IareaT(i,j)/( h(i,j,k) + gv%H_subroundoff))

        if (tracer%id_lbdxy_conc > 0  .or. tracer%id_lbdxy_cont > 0 .or. tracer%id_lbdxy_cont_2d > 0 ) then
          tendency(i,j,k) = (( (uflx(i-1,j,k)-uflx(i,j,k)) ) + ( (vflx(i,j-1,k)-vflx(i,j,k) ) )) * g%IareaT(i,j) * idt
        endif

      endif
    enddo ; enddo ; enddo

    ! Post the tracer diagnostics
    if (tracer%id_lbd_dfx>0)      call post_data(tracer%id_lbd_dfx, uflx*idt, cs%diag)
    if (tracer%id_lbd_dfy>0)      call post_data(tracer%id_lbd_dfy, vflx*idt, cs%diag)
    if (tracer%id_lbd_dfx_2d>0) then
      uwork_2d(:,:) = 0.
      do k=1,gv%ke ; do j=g%jsc,g%jec ; do i=g%isc-1,g%iec
        uwork_2d(i,j) = uwork_2d(i,j) + (uflx(i,j,k) * idt)
      enddo ; enddo ; enddo
      call post_data(tracer%id_lbd_dfx_2d, uwork_2d, cs%diag)
    endif

    if (tracer%id_lbd_dfy_2d>0) then
      vwork_2d(:,:) = 0.
      do k=1,gv%ke ; do j=g%jsc-1,g%jec ; do i=g%isc,g%iec
        vwork_2d(i,j) = vwork_2d(i,j) + (vflx(i,j,k) * idt)
      enddo ; enddo ; enddo
      call post_data(tracer%id_lbd_dfy_2d, vwork_2d, cs%diag)
    endif

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

    ! post depth summed tendency for tracer content
    if (tracer%id_lbdxy_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_lbdxy_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 and its units.
    if (tracer%id_lbdxy_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_lbdxy_conc, tendency, cs%diag)
    endif

  enddo

end subroutine lateral_boundary_diffusion

!< Calculate bulk layer value of a scalar quantity as the thickness weighted average
real function bulk_average(boundary, nk, deg, h, hBLT, phi, ppoly0_E, ppoly0_coefs, method, k_top, zeta_top, k_bot, &
                           zeta_bot)
  integer             :: boundary          !< SURFACE or BOTTOM                                     [nondim]
  integer             :: nk                !< Number of layers                                      [nondim]
  integer             :: deg               !< Degree of polynomial                                  [nondim]
  real, dimension(nk) :: h                 !< Layer thicknesses                                     [H ~> m or kg m-2]
  real                :: hblt              !< Depth of the boundary layer                           [H ~> m or kg m-2]
  real, dimension(nk) :: phi               !< Scalar quantity
  real, dimension(nk,2)     :: ppoly0_e    !< Edge value of polynomial
  real, dimension(nk,deg+1) :: ppoly0_coefs!< Coefficients of polynomial
  integer                   :: method      !< Remapping scheme to use

  integer             :: k_top             !< Index of the first layer within the boundary
  real                :: zeta_top          !< Fraction of the layer encompassed by the bottom boundary layer
                                           !! (0 if none, 1. if all). For the surface, this is always 0. because
                                           !! integration starts at the surface                     [nondim]
  integer             :: k_bot             !< Index of the last layer within the boundary
  real                :: zeta_bot          !< Fraction of the layer encompassed by the surface boundary layer
                                           !! (0 if none, 1. if all). For the bottom boundary layer, this is always 1.
                                           !! because integration starts at the bottom              [nondim]
  ! Local variables
  real    :: htot !< Running sum of the thicknesses (top to bottom)
  integer :: k    !< k indice


  htot = 0.
  bulk_average = 0.
  if (hblt == 0.) return
  if (boundary == surface) then
    htot = (h(k_bot) * zeta_bot)
    bulk_average = average_value_ppoly( nk, phi, ppoly0_e, ppoly0_coefs, method, k_bot, 0., zeta_bot) * htot
    do k = k_bot-1,1,-1
      bulk_average = bulk_average + phi(k)*h(k)
      htot = htot + h(k)
    enddo
  elseif (boundary == bottom) then
    htot = (h(k_top) * zeta_top)
    ! (note 1-zeta_top because zeta_top is the fraction of the layer)
    bulk_average = average_value_ppoly( nk, phi, ppoly0_e, ppoly0_coefs, method, k_top, (1.-zeta_top), 1.) * htot
    do k = k_top+1,nk
      bulk_average = bulk_average + phi(k)*h(k)
      htot = htot + h(k)
    enddo
  else
    call mom_error(fatal, "bulk_average: a valid boundary type must be provided.")
  endif

  bulk_average = bulk_average / hblt

end function bulk_average

!> Calculate the harmonic mean of two quantities
!! See \ref section_harmonic_mean.
real function harmonic_mean(h1,h2)
  real :: h1 !< Scalar quantity
  real :: h2 !< Scalar quantity
  if (h1 + h2 == 0.) then
    harmonic_mean = 0.
  else
    harmonic_mean = 2.*(h1*h2)/(h1+h2)
  endif
end function harmonic_mean

!> Find the k-index range corresponding to the layers that are within the boundary-layer region
subroutine boundary_k_range(boundary, nk, h, hbl, k_top, zeta_top, k_bot, zeta_bot)
  integer,             intent(in   ) :: boundary !< SURFACE or BOTTOM                       [nondim]
  integer,             intent(in   ) :: nk       !< Number of layers                        [nondim]
  real, dimension(nk), intent(in   ) :: h        !< Layer thicknesses of the column         [H ~> m or kg m-2]
  real,                intent(in   ) :: hbl      !< Thickness of the boundary layer         [H ~> m or kg m-2]
                                                 !! If surface, with respect to zbl_ref = 0.
                                                 !! If bottom, with respect to zbl_ref = SUM(h)
  integer,             intent(  out) :: k_top    !< Index of the first layer within the boundary
  real,                intent(  out) :: zeta_top !< Distance from the top of a layer to the intersection of the
                                                 !! top extent of the boundary layer (0 at top, 1 at bottom)  [nondim]
  integer,             intent(  out) :: k_bot    !< Index of the last layer within the boundary
  real,                intent(  out) :: zeta_bot !< Distance of the lower layer to the boundary layer depth
                                                 !! (0 at top, 1 at bottom)  [nondim]
  ! Local variables
  real :: htot ! Summed thickness [H ~> m or kg m-2]
  integer :: k
  ! Surface boundary layer
  if ( boundary == surface ) then
    k_top = 1
    zeta_top = 0.
    htot = 0.
    k_bot = 1
    zeta_bot = 0.
    if (hbl == 0.) return
    if (hbl >= sum(h(:))) then
      k_bot = nk
      zeta_bot = 1.
      return
    endif
    do k=1,nk
      htot = htot + h(k)
      if ( htot >= hbl) then
        k_bot = k
        zeta_bot = 1 - (htot - hbl)/h(k)
        return
      endif
    enddo
  ! Bottom boundary layer
  elseif ( boundary == bottom ) then
    k_top = nk
    zeta_top = 1.
    k_bot = nk
    zeta_bot = 0.
    htot = 0.
    if (hbl == 0.) return
    if (hbl >= sum(h(:))) then
      k_top = 1
      zeta_top = 1.
      return
    endif
    do k=nk,1,-1
      htot = htot + h(k)
      if (htot >= hbl) then
        k_top = k
        zeta_top = 1 - (htot - hbl)/h(k)
        return
      endif
    enddo
  else
    call mom_error(fatal,"Houston, we've had a problem in boundary_k_range")
  endif

end subroutine boundary_k_range


!> Calculate the lateral boundary diffusive fluxes using the layer by layer method.
!! See \ref section_method1
subroutine fluxes_layer_method(boundary, nk, deg, h_L, h_R, hbl_L, hbl_R, area_L, area_R, phi_L, phi_R, &
                              ppoly0_coefs_L, ppoly0_coefs_R, ppoly0_E_L, ppoly0_E_R, method, khtr_u, &
                              F_layer, linear_decay)

  integer,                   intent(in   )       :: boundary !< Which boundary layer SURFACE or BOTTOM  [nondim]
  integer,                   intent(in   )       :: nk       !< Number of layers                        [nondim]
  integer,                   intent(in   )       :: deg      !< order of the polynomial reconstruction  [nondim]
  real, dimension(nk),       intent(in   )       :: h_L      !< Layer thickness (left)              [H ~> m or kg m-2]
  real, dimension(nk),       intent(in   )       :: h_R      !< Layer thickness (right)             [H ~> m or kg m-2]
  real,                      intent(in   )       :: hbl_L    !< Thickness of the boundary boundary
                                                                       !! layer (left)              [H ~> m or kg m-2]
  real,                      intent(in   )       :: hbl_R    !< Thickness of the boundary boundary
                                                             !! layer (right)                       [H ~> m or kg m-2]
  real,                      intent(in   )       :: area_L   !< Area of the horizontal grid (left)  [L2 ~> m2]
  real,                      intent(in   )       :: area_R   !< Area of the horizontal grid (right) [L2 ~> m2]
  real, dimension(nk),       intent(in   )       :: phi_L    !< Tracer values (left)                [conc]
  real, dimension(nk),       intent(in   )       :: phi_R    !< Tracer values (right)               [conc]
  real, dimension(nk,deg+1), intent(in   )       :: ppoly0_coefs_L !< Tracer reconstruction (left)  [conc]
  real, dimension(nk,deg+1), intent(in   )       :: ppoly0_coefs_R !< Tracer reconstruction (right) [conc]
  real, dimension(nk,2),     intent(in   )       :: ppoly0_E_L !< Polynomial edge values (left)     [nondim]
  real, dimension(nk,2),     intent(in   )       :: ppoly0_E_R !< Polynomial edge values (right)    [nondim]
  integer,                   intent(in   )       :: method   !< Method of polynomial integration    [nondim]
  real,                      intent(in   )       :: khtr_u   !< Horizontal diffusivities times delta t
                                                             !! at a velocity point [L2 ~> m2]
  real, dimension(nk),       intent(  out)       :: F_layer  !< Layerwise diffusive flux at U- or V-point
                                                             !! [H L2 conc ~> m3 conc]
  logical, optional,         intent(in   )       :: linear_decay !< If True, apply a linear transition at the base of
                                                             !! the boundary layer
  ! Local variables
  real, dimension(nk) :: h_means              !< Calculate the layer-wise harmonic means           [H ~> m or kg m-2]
  real                :: khtr_avg             !< Thickness-weighted diffusivity at the u-point     [m^2 s^-1]
                                              !! This is just to remind developers that khtr_avg should be
                                              !! computed once khtr is 3D.
  real                :: heff                 !< Harmonic mean of layer thicknesses           [H ~> m or kg m-2]
  real                :: inv_heff             !< Inverse of the harmonic mean of layer thicknesses
                                              !!  [H-1 ~> m-1 or m2 kg-1]
  real                :: phi_L_avg, phi_R_avg !< Bulk, thickness-weighted tracer averages (left and right column)
                                              !!                                                   [conc m^-3 ]
  real    :: htot                      !< Total column thickness [H ~> m or kg m-2]
  real    :: heff_tot                  !< Total effective column thickness in the transition layer [m]
  integer :: k, k_bot_min, k_top_max   !< k-indices, min and max for bottom and top, respectively
  integer :: k_bot_max, k_top_min      !< k-indices, max and min for bottom and top, respectively
  integer :: k_bot_diff, k_top_diff    !< different between left and right k-indices for bottom and top, respectively
  integer :: k_top_L, k_bot_L          !< k-indices left
  integer :: k_top_R, k_bot_R          !< k-indices right
  real    :: zeta_top_L, zeta_top_R    !< distance from the top of a layer to the boundary
                                       !! layer depth                                     [nondim]
  real    :: zeta_bot_L, zeta_bot_R    !< distance from the bottom of a layer to the boundary
                                       !!layer depth                                      [nondim]
  real    :: h_work_L, h_work_R  !< dummy variables
  real    :: hbl_min             !< minimum BLD (left and right)                          [m]
  real    :: wgt                 !< weight to be used in the linear transition to the interior [nondim]
  real    :: a                   !< coefficient to be used in the linear transition to the interior [nondim]
  logical :: linear              !< True if apply a linear transition

  f_layer(:) = 0.0
  if (hbl_l == 0. .or. hbl_r == 0.) then
    return
  endif

  linear = .false.
  if (PRESENT(linear_decay)) then
    linear = linear_decay
  endif

  ! Calculate vertical indices containing the boundary layer
  call boundary_k_range(boundary, nk, h_l, hbl_l, k_top_l, zeta_top_l, k_bot_l, zeta_bot_l)
  call boundary_k_range(boundary, nk, h_r, hbl_r, k_top_r, zeta_top_r, k_bot_r, zeta_bot_r)

  if (boundary == surface) then
    k_bot_min = min(k_bot_l, k_bot_r)
    k_bot_max = max(k_bot_l, k_bot_r)
    k_bot_diff = (k_bot_max - k_bot_min)

    ! make sure left and right k indices span same range
    if (k_bot_min .ne. k_bot_l) then
      k_bot_l = k_bot_min
      zeta_bot_l = 1.0
    endif
    if (k_bot_min .ne. k_bot_r) then
      k_bot_r= k_bot_min
      zeta_bot_r = 1.0
    endif

    h_work_l = (h_l(k_bot_l) * zeta_bot_l)
    h_work_r = (h_r(k_bot_r) * zeta_bot_r)

    phi_l_avg = average_value_ppoly( nk, phi_l, ppoly0_e_l, ppoly0_coefs_l, method, k_bot_l, 0., zeta_bot_l)
    phi_r_avg = average_value_ppoly( nk, phi_r, ppoly0_e_r, ppoly0_coefs_r, method, k_bot_r, 0., zeta_bot_r)
    heff = harmonic_mean(h_work_l, h_work_r)
    ! tracer flux where the minimum BLD intersets layer
    ! GMM, khtr_avg should be computed once khtr is 3D
    if ((linear) .and. (k_bot_diff .gt. 1)) then
      ! apply linear decay at the base of hbl
      do k = k_bot_min,1,-1
        heff = harmonic_mean(h_l(k), h_r(k))
        f_layer(k) = -(heff * khtr_u) * (phi_r(k) - phi_l(k))
      enddo
      ! heff_total
      heff_tot = 0.0
      do k = k_bot_min+1,k_bot_max, 1
        heff_tot = heff_tot + harmonic_mean(h_l(k), h_r(k))
      enddo

      a = -1.0/heff_tot
      heff_tot = 0.0
      do k = k_bot_min+1,k_bot_max, 1
        heff = harmonic_mean(h_l(k), h_r(k))
        wgt = (a*(heff_tot + (heff * 0.5))) + 1.0
        f_layer(k) = -(heff * khtr_u) * (phi_r(k) - phi_l(k)) * wgt
        heff_tot = heff_tot + heff
      enddo
    else
      f_layer(k_bot_min) = -(heff * khtr_u) * (phi_r_avg - phi_l_avg)
      do k = k_bot_min-1,1,-1
        heff = harmonic_mean(h_l(k), h_r(k))
        f_layer(k) = -(heff * khtr_u) * (phi_r(k) - phi_l(k))
      enddo
    endif
  endif

  if (boundary == bottom) then
    ! TODO: GMM add option to apply linear decay
    k_top_max = max(k_top_l, k_top_r)
    ! make sure left and right k indices span same range
    if (k_top_max .ne. k_top_l) then
      k_top_l = k_top_max
      zeta_top_l = 1.0
    endif
    if (k_top_max .ne. k_top_r) then
      k_top_r= k_top_max
      zeta_top_r = 1.0
    endif

    h_work_l = (h_l(k_top_l) * zeta_top_l)
    h_work_r = (h_r(k_top_r) * zeta_top_r)

    phi_l_avg = average_value_ppoly( nk, phi_l, ppoly0_e_l, ppoly0_coefs_l, method, k_top_l, 1.0-zeta_top_l, 1.0)
    phi_r_avg = average_value_ppoly( nk, phi_r, ppoly0_e_r, ppoly0_coefs_r, method, k_top_r, 1.0-zeta_top_r, 1.0)
    heff = harmonic_mean(h_work_l, h_work_r)

    ! tracer flux where the minimum BLD intersets layer
    f_layer(k_top_max) = (-heff * khtr_u) * (phi_r_avg - phi_l_avg)

    do k = k_top_max+1,nk
      heff = harmonic_mean(h_l(k), h_r(k))
      f_layer(k) = -(heff * khtr_u) * (phi_r(k) - phi_l(k))
    enddo
  endif

end subroutine fluxes_layer_method

!> Apply the lateral boundary diffusive fluxes calculated from a 'bulk model'
!! See \ref section_method2
subroutine fluxes_bulk_method(boundary, nk, deg, h_L, h_R, hbl_L, hbl_R, area_L, area_R, phi_L, phi_R, ppoly0_coefs_L, &
                              ppoly0_coefs_R, ppoly0_E_L, ppoly0_E_R, method, khtr_u, F_bulk, F_layer, F_limit, &
                              linear_decay)

  integer,                   intent(in   )       :: boundary !< Which boundary layer SURFACE or BOTTOM  [nondim]
  integer,                   intent(in   )       :: nk       !< Number of layers                        [nondim]
  integer,                   intent(in   )       :: deg      !< order of the polynomial reconstruction  [nondim]
  real, dimension(nk),       intent(in   )       :: h_L      !< Layer thickness (left)              [H ~> m or kg m-2]
  real, dimension(nk),       intent(in   )       :: h_R      !< Layer thickness (right)             [H ~> m or kg m-2]
  real,                      intent(in   )       :: hbl_L    !< Thickness of the boundary boundary
                                                                       !! layer (left)              [H ~> m or kg m-2]
  real,                      intent(in   )       :: hbl_R    !< Thickness of the boundary boundary
                                                             !! layer (left)                        [H ~> m or kg m-2]
  real,                      intent(in   )       :: area_L   !< Area of the horizontal grid (left)  [L2 ~> m2]
  real,                      intent(in   )       :: area_R   !< Area of the horizontal grid (right) [L2 ~> m2]
  real, dimension(nk),       intent(in   )       :: phi_L    !< Tracer values (left)                [conc]
  real, dimension(nk),       intent(in   )       :: phi_R    !< Tracer values (right)               [conc]
  real, dimension(nk,deg+1), intent(in   )       :: ppoly0_coefs_L !< Tracer reconstruction (left)  [conc]
  real, dimension(nk,deg+1), intent(in   )       :: ppoly0_coefs_R !< Tracer reconstruction (right) [conc]
  real, dimension(nk,2),     intent(in   )       :: ppoly0_E_L !< Polynomial edge values (left)     [nondim]
  real, dimension(nk,2),     intent(in   )       :: ppoly0_E_R !< Polynomial edge values (right)    [nondim]
  integer,                   intent(in   )       :: method   !< Method of polynomial integration    [nondim]
  real,                      intent(in   )       :: khtr_u   !< Horizontal diffusivities times delta t
                                                             !! at a velocity point [L2 ~> m2]
  real,                      intent(  out)       :: F_bulk   !< The bulk mixed layer lateral flux
                                                             !! [H L2 conc ~> m3 conc]
  real, dimension(nk),       intent(  out)       :: F_layer  !< Layerwise diffusive flux at U- or V-point
                                                             !! [H L2 conc ~> m3 conc]
  logical, optional,         intent(in   )       :: F_limit  !< If True, apply a limiter
  logical, optional,         intent(in   )       :: linear_decay !< If True, apply a linear transition at the base of
                                                             !! the boundary layer

  ! Local variables
  real, dimension(nk) :: h_means              !< Calculate the layer-wise harmonic means           [H ~> m or kg m-2]
  real                :: khtr_avg             !< Thickness-weighted diffusivity at the u-point     [m^2 s^-1]
                                              !! This is just to remind developers that khtr_avg should be
                                              !! computed once khtr is 3D.
  real                :: heff                 !< Harmonic mean of layer thicknesses                [H ~> m or kg m-2]
  real                :: heff_tot             !< Total effective column thickness in the transition layer [m]
  real                :: inv_heff             !< Inverse of the harmonic mean of layer thicknesses
                                              !! [H-1 ~> m-1 or m2 kg-1]
  real                :: phi_L_avg, phi_R_avg !< Bulk, thickness-weighted tracer averages (left and right column)
                                              !!                                                   [conc m^-3 ]
  real    :: htot                             ! Total column thickness [H ~> m or kg m-2]
  integer :: k, k_min, k_max                  !< k-indices, min and max for top and bottom, respectively
  integer :: k_diff                           !< difference between k_max and k_min
  integer :: k_top_L, k_bot_L                 !< k-indices left
  integer :: k_top_R, k_bot_R                 !< k-indices right
  real    :: zeta_top_L, zeta_top_R           !< distance from the top of a layer to the
                                              !! boundary layer   [nondim]
  real    :: zeta_bot_L, zeta_bot_R           !< distance from the bottom of a layer to the
                                              !! boundary layer   [nondim]
  real    :: h_work_L, h_work_R               !< dummy variables
  real    :: F_max                            !< The maximum amount of flux that can leave a
                                              !! cell  [m^3 conc]
  logical :: limiter                          !< True if flux limiter should be applied
  logical :: linear                           !< True if apply a linear transition
  real    :: hfrac                            !< Layer fraction wrt sum of all layers [nondim]
  real    :: dphi                             !< tracer gradient                      [conc m^-3]
  real    :: wgt                              !< weight to be used in the linear transition to the
                                              !! interior [nondim]
  real    :: a                                !< coefficient to be used in the linear transition to the
                                              !! interior [nondim]

  f_bulk = 0.
  f_layer(:) = 0.
  if (hbl_l == 0. .or. hbl_r == 0.) then
    return
  endif

  limiter = .false.
  if (PRESENT(f_limit)) then
    limiter = f_limit
  endif
  linear = .false.
  if (PRESENT(linear_decay)) then
    linear = linear_decay
  endif

  ! Calculate vertical indices containing the boundary layer
  call boundary_k_range(boundary, nk, h_l, hbl_l, k_top_l, zeta_top_l, k_bot_l, zeta_bot_l)
  call boundary_k_range(boundary, nk, h_r, hbl_r, k_top_r, zeta_top_r, k_bot_r, zeta_bot_r)

  ! Calculate bulk averages of various quantities
  phi_l_avg  = bulk_average(boundary, nk, deg, h_l, hbl_l, phi_l, ppoly0_e_l, ppoly0_coefs_l, method, k_top_l, &
                            zeta_top_l, k_bot_l, zeta_bot_l)
  phi_r_avg  = bulk_average(boundary, nk, deg, h_r, hbl_r, phi_r, ppoly0_e_r, ppoly0_coefs_r, method, k_top_r, &
                            zeta_top_r, k_bot_r, zeta_bot_r)
  ! Calculate the 'bulk' diffusive flux from the bulk averaged quantities
  ! GMM, khtr_avg should be computed once khtr is 3D
  heff = harmonic_mean(hbl_l, hbl_r)
  f_bulk = -(khtr_u * heff) * (phi_r_avg - phi_l_avg)
  ! Calculate the layerwise sum of the vertical effective thickness. This is different than the heff calculated
  ! above, but is used as a way to decompose the fluxes onto the individual layers
  h_means(:) = 0.
  if (boundary == surface) then
    k_min = min(k_bot_l, k_bot_r)
    k_max = max(k_bot_l, k_bot_r)
    k_diff = (k_max - k_min)
    if ((linear) .and. (k_diff .gt. 1)) then
      do k=1,k_min
        h_means(k) = harmonic_mean(h_l(k),h_r(k))
      enddo
      ! heff_total
      heff_tot = 0.0
      do k = k_min+1,k_max, 1
        heff_tot = heff_tot + harmonic_mean(h_l(k), h_r(k))
      enddo

      a = -1.0/heff_tot
      heff_tot = 0.0
      ! fluxes will decay linearly at base of hbl
      do k = k_min+1,k_max, 1
        heff = harmonic_mean(h_l(k), h_r(k))
        wgt = (a*(heff_tot + (heff * 0.5))) + 1.0
        h_means(k) = harmonic_mean(h_l(k), h_r(k)) * wgt
        heff_tot = heff_tot + heff
      enddo
    else
      ! left hand side
      if (k_bot_l == k_min) then
        h_work_l = h_l(k_min) * zeta_bot_l
      else
        h_work_l = h_l(k_min)
      endif

      ! right hand side
      if (k_bot_r == k_min) then
        h_work_r = h_r(k_min) * zeta_bot_r
      else
        h_work_r = h_r(k_min)
      endif

      h_means(k_min) = harmonic_mean(h_work_l,h_work_r)

      do k=1,k_min-1
        h_means(k) = harmonic_mean(h_l(k),h_r(k))
      enddo
    endif

  elseif (boundary == bottom) then
    !TODO, GMM linear decay is not implemented here
    k_max = max(k_top_l, k_top_r)
    ! left hand side
    if (k_top_l == k_max) then
      h_work_l = h_l(k_max) * zeta_top_l
    else
      h_work_l = h_l(k_max)
    endif

    ! right hand side
    if (k_top_r == k_max) then
      h_work_r = h_r(k_max) * zeta_top_r
    else
      h_work_r = h_r(k_max)
    endif

    h_means(k_max) = harmonic_mean(h_work_l,h_work_r)

    do k=nk,k_max+1,-1
      h_means(k) = harmonic_mean(h_l(k),h_r(k))
    enddo
  endif

  if ( sum(h_means) == 0. .or. f_bulk == 0.) then
    return
  ! Decompose the bulk flux onto the individual layers
  else
    ! Initialize remaining thickness
    inv_heff = 1./sum(h_means)
    do k=1,nk
      if ((h_means(k) > 0.) .and. (phi_l(k) /= phi_r(k))) then
        hfrac = h_means(k)*inv_heff
        f_layer(k) = f_bulk * hfrac

        if (limiter) then
          ! limit the flux to 0.2 of the tracer *gradient*
          ! Why 0.2?
          !  t=0         t=inf
          !   0           .2
          ! 0 1 0       .2.2.2
          !   0           .2
          !
          f_max = -0.2 * ((area_r*(phi_r(k)*h_r(k)))-(area_l*(phi_l(k)*h_r(k))))

          ! check if bulk flux (or F_layer) and F_max have same direction
          if ( sign(1.,f_bulk) == sign(1., f_max)) then
            ! Apply flux limiter calculated above
            if (f_max >= 0.) then
              f_layer(k) = min(f_layer(k),f_max)
            else
              f_layer(k) = max(f_layer(k),f_max)
            endif
          else
            ! do not apply a flux on this layer
            f_layer(k) = 0.
          endif
        else
          dphi = -(phi_r(k) - phi_l(k))
          if (.not. sign(1.,f_bulk) == sign(1., dphi)) then
            ! upgradient, do not apply a flux on this layer
            f_layer(k) = 0.
          endif
        endif ! limited
      endif
    enddo
  endif

end subroutine fluxes_bulk_method

!> Unit tests for near-boundary horizontal mixing
logical function near_boundary_unit_tests( verbose )
  logical,               intent(in) :: verbose !< If true, output additional information for debugging unit tests

  ! Local variables
  integer, parameter    :: nk = 2               ! Number of layers
  integer, parameter    :: deg = 1              ! Degree of reconstruction (linear here)
  integer, parameter    :: method = 1           ! Method used for integrating polynomials
  real, dimension(nk)   :: phi_l, phi_r         ! Tracer values (left and right column)             [ nondim m^-3 ]
  real, dimension(nk)   :: phi_l_avg, phi_r_avg ! Bulk, thickness-weighted tracer averages (left and right column)
  real, dimension(nk,deg+1) :: phi_pp_l, phi_pp_r   ! Coefficients for the linear pseudo-reconstructions
                                                !                                                   [ nondim m^-3 ]

  real, dimension(nk,2) :: ppoly0_e_l, ppoly0_e_r! Polynomial edge values (left and right)          [concentration]
  real, dimension(nk)   :: h_l, h_r             ! Layer thickness (left and right)                  [m]
  real                  :: khtr_u               ! Horizontal diffusivities at U-point               [m^2 s^-1]
  real                  :: hbl_l, hbl_r         ! Depth of the boundary layer (left and right)      [m]
  real                  :: f_bulk               ! Total diffusive flux across the U point           [nondim s^-1]
  real, dimension(nk)   :: f_layer              ! Diffusive flux within each layer at U-point       [nondim s^-1]
  real                  :: h_u, hblt_u          ! Thickness at the u-point                          [m]
  real                  :: khtr_avg             ! Thickness-weighted diffusivity at the u-point     [m^2 s^-1]
  real                  :: heff                 ! Harmonic mean of layer thicknesses                [m]
  real                  :: inv_heff             ! Inverse of the harmonic mean of layer thicknesses [m^[-1]
  character(len=120)    :: test_name            ! Title of the unit test
  integer               :: k_top                ! Index of cell containing top of boundary
  real                  :: zeta_top             ! Nondimension position
  integer               :: k_bot                ! Index of cell containing bottom of boundary
  real                  :: zeta_bot             ! Nondimension position
  real                  :: area_l,area_r        ! Area of grid cell [m^2]
  area_l = 1.; area_r = 1. ! Set to unity for all unit tests

  near_boundary_unit_tests = .false.

  ! Unit tests for boundary_k_range
  test_name = 'Surface boundary spans the entire top cell'
  h_l = (/5.,5./)
  call boundary_k_range(surface, nk, h_l, 5., k_top, zeta_top, k_bot, zeta_bot)
  near_boundary_unit_tests = near_boundary_unit_tests .or. &
                             test_boundary_k_range(k_top, zeta_top, k_bot, zeta_bot, 1, 0., 1, 1., test_name, verbose)

  test_name = 'Surface boundary spans the entire column'
  h_l = (/5.,5./)
  call boundary_k_range(surface, nk, h_l, 10., k_top, zeta_top, k_bot, zeta_bot)
  near_boundary_unit_tests = near_boundary_unit_tests .or. &
                             test_boundary_k_range(k_top, zeta_top, k_bot, zeta_bot, 1, 0., 2, 1., test_name, verbose)

  test_name = 'Bottom boundary spans the entire bottom cell'
  h_l = (/5.,5./)
  call boundary_k_range(bottom, nk, h_l, 5., k_top, zeta_top, k_bot, zeta_bot)
  near_boundary_unit_tests = near_boundary_unit_tests .or. &
                             test_boundary_k_range(k_top, zeta_top, k_bot, zeta_bot, 2, 1., 2, 0., test_name, verbose)

  test_name = 'Bottom boundary spans the entire column'
  h_l = (/5.,5./)
  call boundary_k_range(bottom, nk, h_l, 10., k_top, zeta_top, k_bot, zeta_bot)
  near_boundary_unit_tests = near_boundary_unit_tests .or. &
                             test_boundary_k_range(k_top, zeta_top, k_bot, zeta_bot, 1, 1., 2, 0., test_name, verbose)

  test_name = 'Surface boundary intersects second layer'
  h_l = (/10.,10./)
  call boundary_k_range(surface, nk, h_l, 17.5, k_top, zeta_top, k_bot, zeta_bot)
  near_boundary_unit_tests = near_boundary_unit_tests .or. &
                             test_boundary_k_range(k_top, zeta_top, k_bot, zeta_bot, 1, 0., 2, 0.75, test_name, verbose)

  test_name = 'Surface boundary intersects first layer'
  h_l = (/10.,10./)
  call boundary_k_range(surface, nk, h_l, 2.5, k_top, zeta_top, k_bot, zeta_bot)
  near_boundary_unit_tests = near_boundary_unit_tests .or. &
                             test_boundary_k_range(k_top, zeta_top, k_bot, zeta_bot, 1, 0., 1, 0.25, test_name, verbose)

  test_name = 'Surface boundary is deeper than column thickness'
  h_l = (/10.,10./)
  call boundary_k_range(surface, nk, h_l, 21.0, k_top, zeta_top, k_bot, zeta_bot)
  near_boundary_unit_tests = near_boundary_unit_tests .or. &
                             test_boundary_k_range(k_top, zeta_top, k_bot, zeta_bot, 1, 0., 2, 1., test_name, verbose)

  test_name = 'Bottom boundary intersects first layer'
  h_l = (/10.,10./)
  call boundary_k_range(bottom, nk, h_l, 17.5, k_top, zeta_top, k_bot, zeta_bot)
  near_boundary_unit_tests = near_boundary_unit_tests .or. &
                             test_boundary_k_range(k_top, zeta_top, k_bot, zeta_bot, 1, 0.75, 2, 0., test_name, verbose)

  test_name = 'Bottom boundary intersects second layer'
  h_l = (/10.,10./)
  call boundary_k_range(bottom, nk, h_l, 2.5, k_top, zeta_top, k_bot, zeta_bot)
  near_boundary_unit_tests = near_boundary_unit_tests .or. &
                             test_boundary_k_range(k_top, zeta_top, k_bot, zeta_bot, 2, 0.25, 2, 0., test_name, verbose)

  ! All cases in this section have hbl which are equal to the column thicknesses
  test_name = 'Equal hbl and same layer thicknesses (gradient from right to left)'
  hbl_l = 10; hbl_r = 10
  h_l = (/5.,5./) ; h_r = (/5.,5./)
  phi_l = (/0.,0./) ; phi_r = (/1.,1./)
  phi_pp_l(1,1) = 0.; phi_pp_l(1,2) = 0.
  phi_pp_l(2,1) = 0.; phi_pp_l(2,2) = 0.
  phi_pp_r(1,1) = 1.; phi_pp_r(1,2) = 0.
  phi_pp_r(2,1) = 1.; phi_pp_r(2,2) = 0.
  ppoly0_e_l(1,1) = 0.; ppoly0_e_l(1,2) = 0.
  ppoly0_e_l(2,1) = 0.; ppoly0_e_l(2,2) = 0.
  ppoly0_e_r(1,1) = 1.; ppoly0_e_r(1,2) = 1.
  ppoly0_e_r(2,1) = 1.; ppoly0_e_r(2,2) = 1.
  khtr_u = 1.
  ! Without limiter
  call fluxes_bulk_method(surface, nk, deg, h_l, h_r, hbl_l, hbl_r, area_l, area_r, phi_l, phi_r, phi_pp_l, phi_pp_r, &
                                    ppoly0_e_l, ppoly0_e_r, method, khtr_u, f_bulk, f_layer)
  near_boundary_unit_tests = near_boundary_unit_tests .or. &
                             test_layer_fluxes( verbose, nk, test_name, f_layer, (/-5.0,-5.0/) )

  ! same as above, but with limiter
  call fluxes_bulk_method(surface, nk, deg, h_l, h_r, hbl_l, hbl_r, area_l, area_r, phi_l, phi_r, phi_pp_l, phi_pp_r, &
                                    ppoly0_e_l, ppoly0_e_r, method, khtr_u, f_bulk, f_layer, .true.)
  near_boundary_unit_tests = near_boundary_unit_tests .or. &
                             test_layer_fluxes( verbose, nk, test_name, f_layer, (/-1.0,-1.0/) )

  test_name = 'Equal hbl and same layer thicknesses (gradient from left to right)'
  hbl_l = 10.; hbl_r = 10.
  h_l = (/5.,5./) ; h_r = (/5.,5./)
  phi_l = (/1.,1./) ; phi_r = (/0.,0./)
  phi_pp_l(1,1) = 1.; phi_pp_l(1,2) = 0.
  phi_pp_l(2,1) = 1.; phi_pp_l(2,2) = 0.
  phi_pp_r(1,1) = 0.; phi_pp_r(1,2) = 0.
  phi_pp_r(2,1) = 0.; phi_pp_r(2,2) = 0.
  ppoly0_e_l(1,1) = 1.; ppoly0_e_l(1,2) = 1.
  ppoly0_e_l(2,1) = 1.; ppoly0_e_l(2,2) = 1.
  ppoly0_e_r(1,1) = 0.; ppoly0_e_r(1,2) = 0.
  ppoly0_e_r(2,1) = 0.; ppoly0_e_r(2,2) = 0.
  khtr_u = 1.
  call fluxes_bulk_method(surface, nk, deg, h_l, h_r, hbl_l, hbl_r, area_l, area_r, phi_l, phi_r, phi_pp_l, phi_pp_r,&
                                    ppoly0_e_l, ppoly0_e_r, method, khtr_u, f_bulk, f_layer)
  near_boundary_unit_tests = near_boundary_unit_tests .or. &
                             test_layer_fluxes( verbose, nk, test_name, f_layer, (/5.0,5.0/) )

  test_name = 'Equal hbl and same layer thicknesses (no gradient)'
  hbl_l = 10; hbl_r = 10
  h_l = (/5.,5./) ; h_r = (/5.,5./)
  phi_l = (/1.,1./) ; phi_r = (/1.,1./)
  phi_pp_l(1,1) = 1.; phi_pp_l(1,2) = 0.
  phi_pp_l(2,1) = 1.; phi_pp_l(2,2) = 0.
  phi_pp_r(1,1) = 1.; phi_pp_r(1,2) = 0.
  phi_pp_r(2,1) = 1.; phi_pp_r(2,2) = 0.
  ppoly0_e_l(1,1) = 1.; ppoly0_e_l(1,2) = 1.
  ppoly0_e_l(2,1) = 1.; ppoly0_e_l(2,2) = 1.
  ppoly0_e_r(1,1) = 1.; ppoly0_e_r(1,2) = 1.
  ppoly0_e_r(2,1) = 1.; ppoly0_e_r(2,2) = 1.
  khtr_u = 1.
  call fluxes_bulk_method(surface, nk, deg, h_l, h_r, hbl_l, hbl_r, area_l, area_r, phi_l, phi_r, phi_pp_l, phi_pp_r,&
                                    ppoly0_e_l, ppoly0_e_r, method, khtr_u, f_bulk, f_layer)
  near_boundary_unit_tests = near_boundary_unit_tests .or. &
                             test_layer_fluxes( verbose, nk, test_name, f_layer, (/0.0,0.0/) )

  test_name = 'Equal hbl and different layer thicknesses (gradient right to left)'
  hbl_l = 16.; hbl_r = 16.
  h_l = (/10.,6./) ; h_r = (/6.,10./)
  phi_l = (/0.,0./) ; phi_r = (/1.,1./)
  phi_pp_l(1,1) = 0.; phi_pp_l(1,2) = 0.
  phi_pp_l(2,1) = 0.; phi_pp_l(2,2) = 0.
  phi_pp_r(1,1) = 1.; phi_pp_r(1,2) = 0.
  phi_pp_r(2,1) = 1.; phi_pp_r(2,2) = 0.
  ppoly0_e_l(1,1) = 0.; ppoly0_e_l(1,2) = 0.
  ppoly0_e_l(2,1) = 0.; ppoly0_e_l(2,2) = 0.
  ppoly0_e_r(1,1) = 1.; ppoly0_e_r(1,2) = 1.
  ppoly0_e_r(2,1) = 1.; ppoly0_e_r(2,2) = 1.
  khtr_u = 1.
  call fluxes_bulk_method(surface, nk, deg, h_l, h_r, hbl_l, hbl_r, area_l, area_r, phi_l, phi_r, phi_pp_l, phi_pp_r,&
                                    ppoly0_e_l, ppoly0_e_r, method, khtr_u, f_bulk, f_layer)
  near_boundary_unit_tests = near_boundary_unit_tests .or. &
                             test_layer_fluxes( verbose, nk, test_name, f_layer, (/-8.0,-8.0/) )

  test_name = 'Equal hbl and same layer thicknesses (diagonal tracer values)'
  hbl_l = 10.; hbl_r = 10.
  h_l = (/5.,5./) ; h_r = (/5.,5./)
  phi_l = (/1.,0./) ; phi_r = (/0.,1./)
  phi_pp_l(1,1) = 1.; phi_pp_l(1,2) = 0.
  phi_pp_l(2,1) = 0.; phi_pp_l(2,2) = 0.
  phi_pp_r(1,1) = 0.; phi_pp_r(1,2) = 0.
  phi_pp_r(2,1) = 1.; phi_pp_r(2,2) = 0.
  ppoly0_e_l(1,1) = 1.; ppoly0_e_l(1,2) = 1.
  ppoly0_e_l(2,1) = 0.; ppoly0_e_l(2,2) = 0.
  ppoly0_e_r(1,1) = 0.; ppoly0_e_r(1,2) = 0.
  ppoly0_e_r(2,1) = 1.; ppoly0_e_r(2,2) = 1.
  khtr_u = 1.
  call fluxes_bulk_method(surface, nk, deg, h_l, h_r, hbl_l, hbl_r, area_l, area_r, phi_l, phi_r, phi_pp_l, phi_pp_r,&
                                    ppoly0_e_l, ppoly0_e_r, method, khtr_u, f_bulk, f_layer)
  near_boundary_unit_tests = near_boundary_unit_tests .or. &
                             test_layer_fluxes( verbose, nk, test_name, f_layer, (/0.0,0.0/) )

  test_name = 'Different hbl and different layer thicknesses (gradient from right to left)'
  hbl_l = 12; hbl_r = 20
  h_l = (/6.,6./) ; h_r = (/10.,10./)
  phi_l = (/0.,0./) ; phi_r = (/1.,1./)
  phi_pp_l(1,1) = 0.; phi_pp_l(1,2) = 0.
  phi_pp_l(2,1) = 0.; phi_pp_l(2,2) = 0.
  phi_pp_r(1,1) = 1.; phi_pp_r(1,2) = 0.
  phi_pp_r(2,1) = 1.; phi_pp_r(2,2) = 0.
  ppoly0_e_l(1,1) = 0.; ppoly0_e_l(1,2) = 0.
  ppoly0_e_l(2,1) = 0.; ppoly0_e_l(2,2) = 0.
  ppoly0_e_r(1,1) = 1.; ppoly0_e_r(1,2) = 1.
  ppoly0_e_r(2,1) = 1.; ppoly0_e_r(2,2) = 1.
  khtr_u = 1.
  call fluxes_bulk_method(surface, nk, deg, h_l, h_r, hbl_l, hbl_r, area_l, area_r, phi_l, phi_r, phi_pp_l, phi_pp_r,&
                                    ppoly0_e_l, ppoly0_e_r, method, khtr_u, f_bulk, f_layer)
  near_boundary_unit_tests = near_boundary_unit_tests .or. &
                             test_layer_fluxes( verbose, nk, test_name, f_layer, (/-7.5,-7.5/) )

  ! Cases where hbl < column thickness (polynomial coefficients specified for pseudo-linear reconstruction)

  test_name = 'hbl < column thickness, hbl same, constant concentration each column'
  hbl_l = 2; hbl_r = 2
  h_l = (/1.,2./) ; h_r = (/1.,2./)
  phi_l = (/0.,0./) ; phi_r = (/1.,1./)
  phi_pp_l(1,1) = 0.; phi_pp_l(1,2) = 0.
  phi_pp_l(2,1) = 0.; phi_pp_l(2,2) = 0.
  phi_pp_r(1,1) = 1.; phi_pp_r(1,2) = 0.
  phi_pp_r(2,1) = 1.; phi_pp_r(2,2) = 0.
  ppoly0_e_l(1,1) = 0.; ppoly0_e_l(1,2) = 0.
  ppoly0_e_l(2,1) = 0.; ppoly0_e_l(2,2) = 0.
  ppoly0_e_r(1,1) = 1.; ppoly0_e_r(1,2) = 1.
  ppoly0_e_r(2,1) = 1.; ppoly0_e_r(2,2) = 1.
  khtr_u = 1.
  call fluxes_bulk_method(surface, nk, deg, h_l, h_r, hbl_l, hbl_r, area_l, area_r, phi_l, phi_r, phi_pp_l, phi_pp_r,&
                                    ppoly0_e_l, ppoly0_e_r, method, khtr_u, f_bulk, f_layer)
  near_boundary_unit_tests = near_boundary_unit_tests .or. &
                             test_layer_fluxes( verbose, nk, test_name, f_layer, (/-1.,-1./) )

  test_name = 'hbl < column thickness, hbl same, linear profile right'
  hbl_l = 2; hbl_r = 2
  h_l = (/1.,2./) ; h_r = (/1.,2./)
  phi_l = (/0.,0./) ; phi_r = (/0.5,2./)
  phi_pp_l(1,1) = 0.; phi_pp_l(1,2) = 0.
  phi_pp_l(2,1) = 0.; phi_pp_l(2,2) = 0.
  phi_pp_r(1,1) = 0.; phi_pp_r(1,2) = 1.
  phi_pp_r(2,1) = 1.; phi_pp_r(2,2) = 2.
  khtr_u = 1.
  ppoly0_e_l(1,1) = 0.; ppoly0_e_l(1,2) = 0.
  ppoly0_e_l(2,1) = 0.; ppoly0_e_l(2,2) = 0.
  ppoly0_e_r(1,1) = 0.; ppoly0_e_r(1,2) = 1.
  ppoly0_e_r(2,1) = 1.; ppoly0_e_r(2,2) = 3.
  call fluxes_bulk_method(surface, nk, deg, h_l, h_r, hbl_l, hbl_r, area_l, area_r, phi_l, phi_r, phi_pp_l, phi_pp_r,&
                                    ppoly0_e_l, ppoly0_e_r, method, khtr_u, f_bulk, f_layer)
  near_boundary_unit_tests = near_boundary_unit_tests .or. &
                             test_layer_fluxes( verbose, nk, test_name, f_layer, (/-1.,-1./) )

  test_name = 'hbl < column thickness, hbl same, linear profile right, khtr=2'
  hbl_l = 2; hbl_r = 2
  h_l = (/1.,2./) ; h_r = (/1.,2./)
  phi_l = (/0.,0./) ; phi_r = (/0.5,2./)
  phi_pp_l(1,1) = 0.; phi_pp_l(1,2) = 0.
  phi_pp_l(2,1) = 0.; phi_pp_l(2,2) = 0.
  phi_pp_r(1,1) = 0.; phi_pp_r(1,2) = 1.
  phi_pp_r(2,1) = 1.; phi_pp_r(2,2) = 2.
  khtr_u = 2.
  ppoly0_e_l(1,1) = 0.; ppoly0_e_l(1,2) = 0.
  ppoly0_e_l(2,1) = 0.; ppoly0_e_l(2,2) = 0.
  ppoly0_e_r(1,1) = 0.; ppoly0_e_r(1,2) = 1.
  ppoly0_e_r(2,1) = 1.; ppoly0_e_r(2,2) = 3.
  call fluxes_layer_method(surface, nk, deg, h_l, h_r, hbl_l, hbl_r, area_l, area_r, phi_l, phi_r, phi_pp_l, &
                                    phi_pp_r, ppoly0_e_l, ppoly0_e_r, method, khtr_u, f_layer)
  near_boundary_unit_tests = near_boundary_unit_tests .or. &
                             test_layer_fluxes( verbose, nk, test_name, f_layer, (/-1.,-3./) )

  ! unit tests for layer by layer method
  test_name = 'Different hbl and different column thicknesses (gradient from right to left)'
  hbl_l = 12; hbl_r = 20
  h_l = (/6.,6./) ; h_r = (/10.,10./)
  phi_l = (/0.,0./) ; phi_r = (/1.,1./)
  phi_pp_l(1,1) = 0.; phi_pp_l(1,2) = 0.
  phi_pp_l(2,1) = 0.; phi_pp_l(2,2) = 0.
  phi_pp_r(1,1) = 1.; phi_pp_r(1,2) = 0.
  phi_pp_r(2,1) = 1.; phi_pp_r(2,2) = 0.
  ppoly0_e_l(1,1) = 0.; ppoly0_e_l(1,2) = 0.
  ppoly0_e_l(2,1) = 0.; ppoly0_e_l(2,2) = 0.
  ppoly0_e_r(1,1) = 1.; ppoly0_e_r(1,2) = 1.
  ppoly0_e_r(2,1) = 1.; ppoly0_e_r(2,2) = 1.
  khtr_u = 1.
  call fluxes_layer_method(surface, nk, deg, h_l, h_r, hbl_l, hbl_r, area_l, area_r, phi_l, phi_r, phi_pp_l, &
                                    phi_pp_r, ppoly0_e_l, ppoly0_e_r, method, khtr_u, f_layer)
  near_boundary_unit_tests = near_boundary_unit_tests .or. &
                             test_layer_fluxes( verbose, nk, test_name, f_layer, (/-7.5,-7.5/) )

  test_name = 'Different hbl and different column thicknesses (linear profile right)'

  hbl_l = 15; hbl_r = 6
  h_l = (/10.,10./) ; h_r = (/12.,10./)
  phi_l = (/0.,0./) ; phi_r = (/1.,3./)
  phi_pp_l(1,1) = 0.; phi_pp_l(1,2) = 0.
  phi_pp_l(2,1) = 0.; phi_pp_l(2,2) = 0.
  phi_pp_r(1,1) = 0.; phi_pp_r(1,2) = 2.
  phi_pp_r(2,1) = 2.; phi_pp_r(2,2) = 2.
  ppoly0_e_l(1,1) = 0.; ppoly0_e_l(1,2) = 0.
  ppoly0_e_l(2,1) = 0.; ppoly0_e_l(2,2) = 0.
  ppoly0_e_r(1,1) = 0.; ppoly0_e_r(1,2) = 2.
  ppoly0_e_r(2,1) = 2.; ppoly0_e_r(2,2) = 4.
  khtr_u = 1.
  call fluxes_layer_method(surface, nk, deg, h_l, h_r, hbl_l, hbl_r, area_l, area_r, phi_l, phi_r, phi_pp_l, &
                                    phi_pp_r, ppoly0_e_l, ppoly0_e_r, method, khtr_u, f_layer)
  near_boundary_unit_tests = near_boundary_unit_tests .or. &
                             test_layer_fluxes( verbose, nk, test_name, f_layer, (/-3.75,0.0/) )
end function near_boundary_unit_tests

!> Returns true if output of near-boundary unit tests does not match correct computed values
!! and conditionally writes results to stream
logical function test_layer_fluxes(verbose, nk, test_name, F_calc, F_ans)
  logical,                    intent(in) :: verbose   !< If true, write results to stdout
  character(len=80),          intent(in) :: test_name !< Brief description of the unit test
  integer,                    intent(in) :: nk        !< Number of layers
  real, dimension(nk),        intent(in) :: f_calc    !< Fluxes of the unitless tracer from the algorithm [s^-1]
  real, dimension(nk),        intent(in) :: f_ans     !< Fluxes of the unitless tracer calculated by hand [s^-1]
  ! Local variables
  integer :: k
  integer, parameter :: stdunit = stdout

  test_layer_fluxes = .false.
  do k=1,nk
    if ( f_calc(k) /= f_ans(k) ) then
      test_layer_fluxes = .true.
      write(stdunit,*) "MOM_lateral_boundary_diffusion, UNIT TEST FAILED: ", test_name
      write(stdunit,10) k, f_calc(k), f_ans(k)
    elseif (verbose) then
      write(stdunit,10) k, f_calc(k), f_ans(k)
    endif
  enddo

10 format("Layer=",i3," F_calc=",f20.16," F_ans",f20.16)
end function test_layer_fluxes

!> Return true if output of unit tests for boundary_k_range does not match answers
logical function test_boundary_k_range(k_top, zeta_top, k_bot, zeta_bot, k_top_ans, zeta_top_ans,&
                                       k_bot_ans, zeta_bot_ans, test_name, verbose)
  integer :: k_top               !< Index of cell containing top of boundary
  real    :: zeta_top            !< Nondimension position
  integer :: k_bot               !< Index of cell containing bottom of boundary
  real    :: zeta_bot            !< Nondimension position
  integer :: k_top_ans           !< Index of cell containing top of boundary
  real    :: zeta_top_ans        !< Nondimension position
  integer :: k_bot_ans           !< Index of cell containing bottom of boundary
  real    :: zeta_bot_ans        !< Nondimension position
  character(len=80) :: test_name !< Name of the unit test
  logical :: verbose             !< If true always print output

  integer, parameter :: stdunit = stdout

  test_boundary_k_range = k_top .ne. k_top_ans
  test_boundary_k_range = test_boundary_k_range .or. (zeta_top .ne. zeta_top_ans)
  test_boundary_k_range = test_boundary_k_range .or. (k_bot .ne. k_bot_ans)
  test_boundary_k_range = test_boundary_k_range .or. (zeta_bot .ne. zeta_bot_ans)

  if (test_boundary_k_range) write(stdunit,*) "UNIT TEST FAILED: ", test_name
  if (test_boundary_k_range .or. verbose) then
    write(stdunit,20) "k_top", k_top, "k_top_ans", k_top_ans
    write(stdunit,20) "k_bot", k_bot, "k_bot_ans", k_bot_ans
    write(stdunit,30) "zeta_top", zeta_top, "zeta_top_ans", zeta_top_ans
    write(stdunit,30) "zeta_bot", zeta_bot, "zeta_bot_ans", zeta_bot_ans
  endif

  20 format(a,"=",i3,x,a,"=",i3)
  30 format(a,"=",f20.16,x,a,"=",f20.16)


end function test_boundary_k_range

!> \namespace mom_lateral_boundary_diffusion
!!
!! \section section_LBD The Lateral Boundary Diffusion (LBD) framework
!!
!! The LBD framework accounts for the effects of diabatic mesoscale fluxes
!! within surface and bottom boundary layers. Unlike the equivalent adiabatic
!! fluxes, which is applied along neutral density surfaces, LBD is purely
!! horizontal.
!!
!! The bottom boundary layer fluxes remain to be implemented, although most
!! of the steps needed to do so have already been added and tested.
!!
!! Boundary lateral diffusion can be applied using one of the three methods:
!!
!! * [Method #1: Along layer](@ref section_method2) (default);
!! * [Method #2: Bulk layer](@ref section_method1);
!!
!! A brief summary of these methods is provided below.
!!
!! \subsection section_method1 Along layer approach (Method #1)
!!
!! This is the recommended and more straight forward method where diffusion is
!! applied layer by layer using only information from neighboring cells.
!!
!! Step #1: compute vertical indices containing boundary layer (boundary_k_range).
!! For the TOP boundary layer, these are:
!!
!! k_top, k_bot, zeta_top, zeta_bot
!!
!! Step #2: calculate the diffusive flux at each layer:
!!
!! \f[ F_{k} = -KHTR \times h_{eff}(k) \times (\phi_R(k) - \phi_L(k)),  \f]
!! where h_eff is the [harmonic mean](@ref section_harmonic_mean) of the layer thickness
!! in the left and right columns. This method does not require a limiter since KHTR
!! is already limted based on a diffusive CFL condition prior to the call of this
!! module.
!!
!! Step #3: option to linearly decay the flux from k_bot_min to k_bot_max:
!!
!! If LBD_LINEAR_TRANSITION = True and k_bot_diff > 1, the diffusive flux will decay
!! linearly between the top interface of the layer containing the minimum boundary
!! layer depth (k_bot_min) and the lower interface of the layer containing the
!! maximum layer depth (k_bot_max).
!!
!! \subsection section_method2 Bulk layer approach (Method #2)
!!
!! Apply the lateral boundary diffusive fluxes calculated from a 'bulk model'.This
!! is a lower order representation (Kraus-Turner like approach) which assumes that
!! eddies are acting along well mixed layers (i.e., eddies do not know care about
!! vertical tracer gradients within the boundary layer).
!!
!! Step #1: compute vertical indices containing boundary layer (boundary_k_range).
!! For the TOP boundary layer, these are:
!!
!! k_top, k_bot, zeta_top, zeta_bot
!!
!! Step #2: compute bulk averages (thickness weighted) tracer averages (phi_L and phi_R),
!! then calculate the bulk diffusive flux (F_{bulk}):
!!
!! \f[ F_{bulk} = -KHTR \times h_{eff} \times (\phi_R - \phi_L),  \f]
!! where h_eff is the [harmonic mean](@ref section_harmonic_mean) of the boundary layer depth
!! in the left and right columns (\f[ HBL_L \f] and \f[ HBL_R \f], respectively).
!!
!! Step #3: decompose F_bulk onto individual layers:
!!
!! \f[ F_{layer}(k) = F_{bulk} \times h_{frac}(k) ,  \f]
!!
!! where h_{frac} is
!!
!! \f[ h_{frac}(k) = h_u(k) \times \frac{1}{\sum(h_u)}.  \f]
!!
!! h_u is the [harmonic mean](@ref section_harmonic_mean) of thicknesses at each layer.
!! Special care (layer reconstruction) must be taken at k_min = min(k_botL, k_bot_R).
!!
!! Step #4: option to linearly decay the flux from k_bot_min to k_bot_max:
!!
!! If LBD_LINEAR_TRANSITION = True and k_bot_diff > 1, the diffusive flux will decay
!! linearly between the top interface of the layer containing the minimum boundary
!! layer depth (k_bot_min) and the lower interface of the layer containing the
!! maximum layer depth (k_bot_max).
!!
!! Step #5: limit the tracer flux so that 1) only down-gradient fluxes are applied,
!! and 2) the flux cannot be larger than F_max, which is defined using the tracer
!! gradient:
!!
!! \f[ F_{max} = -0.2 \times [(V_R(k) \times \phi_R(k)) - (V_L(k) \times \phi_L(k))],  \f]
!! where V is the cell volume. Why 0.2?
!!          t=0         t=inf
!!           0           .2
!!         0 1 0       .2.2.2
!!           0           .2
!!
!! \subsection section_harmonic_mean Harmonic Mean
!!
!! The harmonic mean (HM) betwen h1 and h2 is defined as:
!!
!! \f[ HM = \frac{2 \times h1 \times h2}{h1 + h2} \f]
!!
end module mom_lateral_boundary_diffusion