MOM_lateral_boundary_diffusion.F90

Go to the documentation of this file.

!> 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
 
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. Bulk-layer approach
                                                                  !! 2. Along layer
  integer :: deg                                                  !< Degree of polynomial reconstruction
  integer :: surface_boundary_scheme                              !< Which boundary layer scheme to use
                                                                  !! 1. ePBL; 2. KPP
  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 log_version(param_file, mdl, version, &
       "This module implements lateral diffusion of tracers near boundaries")
  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) then
    return
  endif
 
  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. Bulk layer approach  \n"//&
                 "2. Along layer approach", default=1)
  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 [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(lateral_boundary_diffusion_cs),  intent(in)    :: cs     !< Control structure for this module
 
  ! Local variables
  real, dimension(SZI_(G),SZJ_(G)) :: hbl                           !< bnd. layer depth [m]
  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)
  if (ASSOCIATED(cs%energetic_PBL_CSp)) call energetic_pbl_get_mld(cs%energetic_PBL_CSp, hbl, g, us)
 
  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
 
    do j = g%jsc-1, g%jec+1
      ! Interpolate state to interface
      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
    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_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,:))
          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,:))
          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)
 
    ! Method #2
    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_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,:))
          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,:))
          endif
        enddo
      enddo
    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.
    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                                         [m]
  real                :: hblt              !< Depth of the boundary layer                               [m]
  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         [m]
  real,                intent(in   ) :: hbl      !< Thickness of the boundary layer         [m]
                                                 !! 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
  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 = 0.
      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 = 1.
    htot = 0.
    if (hbl == 0.) return
    if (hbl >= sum(h(:))) then
      k_top = 1
      zeta_top = 0.
      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_method2
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)
 
  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)                  [m]
  real, dimension(nk),       intent(in   )       :: h_R      !< Layer thickness (right)                 [m]
  real,                      intent(in   )       :: hbl_L    !< Thickness of the boundary boundary
                                                                       !! layer (left)                  [m]
  real,                      intent(in   )       :: hbl_R    !< Thickness of the boundary boundary
                                                             !! layer (right)                           [m]
  real,                      intent(in   )       :: area_L   !< Area of the horizontal grid (left)      [m^2]
  real,                      intent(in   )       :: area_R   !< Area of the horizontal grid (right)     [m^2]
  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 U-point [m^2]
  real, dimension(nk),       intent(  out)       :: F_layer  !< Layerwise diffusive flux at U- or V-point [m^3 conc]
 
  ! Local variables
  real, dimension(nk) :: h_means              !< Calculate the layer-wise harmonic means           [m]
  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                [m]
  real                :: inv_heff             !< Inverse of the harmonic mean of layer thicknesses [m^[-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 [m]
  integer :: k, k_bot_min, k_top_max   !< k-indices, min and max for top and bottom, 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]
 
  f_layer(:) = 0.0
  if (hbl_l == 0. .or. hbl_r == 0.) then
    return
  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)
    ! 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
    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
 
  if (boundary == bottom) then
    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_method1
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)
 
  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)                  [m]
  real, dimension(nk),       intent(in   )       :: h_R      !< Layer thickness (right)                 [m]
  real,                      intent(in   )       :: hbl_L    !< Thickness of the boundary boundary
                                                                       !! layer (left)                  [m]
  real,                      intent(in   )       :: hbl_R    !< Thickness of the boundary boundary
                                                             !! layer (left)                            [m]
  real,                      intent(in   )       :: area_L   !< Area of the horizontal grid (left)      [m^2]
  real,                      intent(in   )       :: area_R   !< Area of the horizontal grid (right)     [m^2]
  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 U-point [m^2]
  real,                      intent(  out)       :: F_bulk   !< The bulk mixed layer lateral flux       [m^3 conc]
  real, dimension(nk),       intent(  out)       :: F_layer  !< Layerwise diffusive flux at U-point     [m^3 conc]
  real, optional, dimension(nk), intent(  out)   :: F_limit  !< The amount of flux not applied due to limiter
                                                             !! F_layer(k) - F_max                      [m^3 conc]
  ! Local variables
  real, dimension(nk) :: h_means              !< Calculate the layer-wise harmonic means           [m]
  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                [m]
  real                :: inv_heff             !< Inverse of the harmonic mean of layer thicknesses [m^[-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 [m]
  integer :: k, k_min, k_max                  !< k-indices, min and max for top and bottom, 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   [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 :: limited                          !< True if the flux limiter was applied
  real    :: hfrac, F_bulk_remain
 
  if (hbl_l == 0. .or. hbl_r == 0.) then
    f_bulk = 0.
    f_layer(:) = 0.
    return
  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)
  f_bulk_remain = f_bulk
  ! 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)
 
    ! 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
 
  elseif (boundary == bottom) then
    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. ) 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.) then
        hfrac = h_means(k)*inv_heff
        f_layer(k) = f_bulk * hfrac
        ! 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
          ! Distribute bulk flux onto layers
          if ( ((boundary == surface) .and. (k == k_min)) .or. ((boundary == bottom) .and. (k == nk)) ) then
            f_layer(k) = f_bulk_remain ! GMM, are not using F_bulk_remain for now. Should we keep it?
          endif
          f_bulk_remain = f_bulk_remain - f_layer(k)
 
          ! Apply flux limiter calculated above
          if (f_max >= 0.) then
            limited = f_layer(k) > f_max
            f_layer(k) = min(f_layer(k),f_max)
          else
            limited = f_layer(k) < f_max
            f_layer(k) = max(f_layer(k),f_max)
          endif
 
          ! GMM, again we are not using F_limit. Should we delete it?
          if (PRESENT(f_limit)) then
            if (limited) then
              f_limit(k) = f_layer(k) - f_max
            else
              f_limit(k) = 0.
            endif
          endif
        else
          ! do not apply a flux on this layer
          f_bulk_remain = f_bulk_remain - f_layer(k)
          f_layer(k) = 0.
        endif
      else
        f_layer(k) = 0.
      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 = 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 = test_boundary_k_range(k_top, zeta_top, k_bot, zeta_bot, 1, 0., 2, 0., 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 = test_boundary_k_range(k_top, zeta_top, k_bot, zeta_bot, 2, 0., 2, 1., 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 = test_boundary_k_range(k_top, zeta_top, k_bot, zeta_bot, 1, 0., 2, 1., 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 = 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 = 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 = test_boundary_k_range(k_top, zeta_top, k_bot, zeta_bot, 1, 0., 2, 0., 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 = test_boundary_k_range(k_top, zeta_top, k_bot, zeta_bot, 1, 0.75, 2, 1., 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 = test_boundary_k_range(k_top, zeta_top, k_bot, zeta_bot, 2, 0.25, 2, 1., 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.
  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 = test_layer_fluxes( verbose, nk, test_name, f_layer, (/-5.0,-5.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 = 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) = 0.
  ppoly0_e_l(2,1) = 1.; 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 = 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 = 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 = 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 = 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.
  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 = 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 = 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 = test_layer_fluxes( verbose, nk, test_name, f_layer, (/-2.,-2./) )
 
  ! 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 = test_layer_fluxes( verbose, nk, test_name, f_layer, (/-7.5,0.0/) )
 
  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 = 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 = 6
 
  test_layer_fluxes = .false.
  do k=1,nk
    if ( f_calc(k) /= f_ans(k) ) then
      test_layer_fluxes = .true.
      write(stdunit,*) "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 = 6
 
  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: Bulk layer](@ref section_method1) (default);
!! * [Method #2: Along layer](@ref section_method2);
!!
!! A brief summary of these methods is provided below.
!!
!! \subsection section_method1 Bulk layer approach (Method #1)
!!
!! 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: 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_method2 Along layer approach (Method #2)
!!
!! This is a 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. Special care (layer reconstruction) must be taken at
!! k_min = min(k_botL, k_bot_R). 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.
!!
!! \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