mom_lateral_boundary_diffusion Module Reference

Detailed Description

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.

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 (default);
• Method #2: Bulk layer;

A brief summary of these methods is provided below.

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_{k} = -KHTR \times h_{eff}(k) \times (\phi_R(k) - \phi_L(k)), \]

where h_eff is the 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).

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_{bulk} = -KHTR \times h_{eff} \times (\phi_R - \phi_L), \]

where h_eff is the harmonic mean of the boundary layer depth in the left and right columns (

\[ HBL_L \]

and

\[ HBL_R \]

, respectively).

Step #3: decompose F_bulk onto individual layers:

\[ F_{layer}(k) = F_{bulk} \times h_{frac}(k) , \]

where h_{frac} is

\[ h_{frac}(k) = h_u(k) \times \frac{1}{\sum(h_u)}. \]

h_u is the 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_{max} = -0.2 \times [(V_R(k) \times \phi_R(k)) - (V_L(k) \times \phi_L(k))], \]

where V is the cell volume. Why 0.2? t=0 t=inf 0 .2 0 1 0 .2.2.2 0 .2

Harmonic Mean

The harmonic mean (HM) betwen h1 and h2 is defined as:

\[ HM = \frac{2 \times h1 \times h2}{h1 + h2} \]

Data Types
type	lateral_boundary_diffusion_cs
	Sets parameters for lateral boundary mixing module. More...

Functions/Subroutines
logical function, public	lateral_boundary_diffusion_init (Time, G, param_file, diag, diabatic_CSp, CS)
	Initialization routine that reads runtime parameters and sets up pointers to other control structures that might be needed for lateral boundary diffusion. More...
subroutine, public	lateral_boundary_diffusion (G, GV, US, h, Coef_x, Coef_y, dt, Reg, CS)
	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. More...
real function	bulk_average (boundary, nk, deg, h, hBLT, phi, ppoly0_E, ppoly0_coefs, method, k_top, zeta_top, k_bot, zeta_bot)
real function	harmonic_mean (h1, h2)
	Calculate the harmonic mean of two quantities See Harmonic Mean. More...
subroutine, public	boundary_k_range (boundary, nk, h, hbl, k_top, zeta_top, k_bot, zeta_bot)
	Find the k-index range corresponding to the layers that are within the boundary-layer region. More...
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)
	Calculate the lateral boundary diffusive fluxes using the layer by layer method. See Along layer approach (Method #1). More...
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)
	Apply the lateral boundary diffusive fluxes calculated from a 'bulk model' See Bulk layer approach (Method #2). More...
logical function, public	near_boundary_unit_tests (verbose)
	Unit tests for near-boundary horizontal mixing. More...
logical function	test_layer_fluxes (verbose, nk, test_name, F_calc, F_ans)
	Returns true if output of near-boundary unit tests does not match correct computed values and conditionally writes results to stream. More...
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)
	Return true if output of unit tests for boundary_k_range does not match answers. More...

Variables
integer, parameter, public	surface = -1
	Set a value that corresponds to the surface bopundary.
integer, parameter, public	bottom = 1
	Set a value that corresponds to the bottom boundary.
character(len=40)	mdl = "MOM_lateral_boundary_diffusion"
	Name of this module.

Function/Subroutine Documentation

boundary_k_range()

subroutine, public mom_lateral_boundary_diffusion::boundary_k_range	(	integer, intent(in)	boundary,
		integer, intent(in)	nk,
		real, dimension(nk), intent(in)	h,
		real, intent(in)	hbl,
		integer, intent(out)	k_top,
		real, intent(out)	zeta_top,
		integer, intent(out)	k_bot,
		real, intent(out)	zeta_bot
	)

Find the k-index range corresponding to the layers that are within the boundary-layer region.

Parameters

[in]	boundary	SURFACE or BOTTOM [nondim]
[in]	nk	Number of layers [nondim]
[in]	h	Layer thicknesses of the column [H ~> m or kg m-2]
[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)
[out]	k_top	Index of the first layer within the boundary
[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]
[out]	k_bot	Index of the last layer within the boundary
[out]	zeta_bot	Distance of the lower layer to the boundary layer depth (0 at top, 1 at bottom) [nondim]

Definition at line 375 of file MOM_lateral_boundary_diffusion.F90.

  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

bulk_average()

real function mom_lateral_boundary_diffusion::bulk_average ( integer  boundary, integer  nk, integer  deg, real, dimension(nk)  h, real  hBLT, real, dimension(nk)  phi, real, dimension(nk,2)  ppoly0_E, real, dimension(nk,deg+1)  ppoly0_coefs, integer  method, integer  k_top, real  zeta_top, integer  k_bot, real  zeta_bot  )

private

Parameters

boundary	SURFACE or BOTTOM [nondim]
nk	Number of layers [nondim]
deg	Degree of polynomial [nondim]
h	Layer thicknesses [H ~> m or kg m-2]
hblt	Depth of the boundary layer [H ~> m or kg m-2]
phi	Scalar quantity
ppoly0_e	Edge value of polynomial
ppoly0_coefs	Coefficients of polynomial
method	Remapping scheme to use
k_top	Index of the first layer within the boundary
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]
k_bot	Index of the last layer within the boundary
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]

Definition at line 312 of file MOM_lateral_boundary_diffusion.F90.

  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

fluxes_bulk_method()

subroutine mom_lateral_boundary_diffusion::fluxes_bulk_method ( integer, intent(in)  boundary, integer, intent(in)  nk, integer, intent(in)  deg, real, dimension(nk), intent(in)  h_L, real, dimension(nk), intent(in)  h_R, real, intent(in)  hbl_L, real, intent(in)  hbl_R, real, intent(in)  area_L, real, intent(in)  area_R, real, dimension(nk), intent(in)  phi_L, real, dimension(nk), intent(in)  phi_R, real, dimension(nk,deg+1), intent(in)  ppoly0_coefs_L, real, dimension(nk,deg+1), intent(in)  ppoly0_coefs_R, real, dimension(nk,2), intent(in)  ppoly0_E_L, real, dimension(nk,2), intent(in)  ppoly0_E_R, integer, intent(in)  method, real, intent(in)  khtr_u, real, intent(out)  F_bulk, real, dimension(nk), intent(out)  F_layer, logical, intent(in), optional  F_limit, logical, intent(in), optional  linear_decay  )

private

Apply the lateral boundary diffusive fluxes calculated from a 'bulk model' See Bulk layer approach (Method #2).

Parameters

[in]	boundary	Which boundary layer SURFACE or BOTTOM [nondim]
[in]	nk	Number of layers [nondim]
[in]	deg	order of the polynomial reconstruction [nondim]
[in]	h_l	Layer thickness (left) [H ~> m or kg m-2]
[in]	h_r	Layer thickness (right) [H ~> m or kg m-2]
[in]	hbl_l	Thickness of the boundary boundary layer (left) [H ~> m or kg m-2]
[in]	hbl_r	Thickness of the boundary boundary layer (left) [H ~> m or kg m-2]
[in]	area_l	Area of the horizontal grid (left) [L2 ~> m2]
[in]	area_r	Area of the horizontal grid (right) [L2 ~> m2]
[in]	phi_l	Tracer values (left) [conc]
[in]	phi_r	Tracer values (right) [conc]
[in]	ppoly0_coefs_l	Tracer reconstruction (left) [conc]
[in]	ppoly0_coefs_r	Tracer reconstruction (right) [conc]
[in]	ppoly0_e_l	Polynomial edge values (left) [nondim]
[in]	ppoly0_e_r	Polynomial edge values (right) [nondim]
[in]	method	Method of polynomial integration [nondim]
[in]	khtr_u	Horizontal diffusivities times delta t at a velocity point [L2 ~> m2]
[out]	f_bulk	The bulk mixed layer lateral flux [H L2 conc ~> m3 conc]
[out]	f_layer	Layerwise diffusive flux at U- or V-point [H L2 conc ~> m3 conc]
[in]	f_limit	If True, apply a limiter
[in]	linear_decay	If True, apply a linear transition at the base of the boundary layer

Definition at line 598 of file MOM_lateral_boundary_diffusion.F90.

  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

fluxes_layer_method()

subroutine mom_lateral_boundary_diffusion::fluxes_layer_method ( integer, intent(in)  boundary, integer, intent(in)  nk, integer, intent(in)  deg, real, dimension(nk), intent(in)  h_L, real, dimension(nk), intent(in)  h_R, real, intent(in)  hbl_L, real, intent(in)  hbl_R, real, intent(in)  area_L, real, intent(in)  area_R, real, dimension(nk), intent(in)  phi_L, real, dimension(nk), intent(in)  phi_R, real, dimension(nk,deg+1), intent(in)  ppoly0_coefs_L, real, dimension(nk,deg+1), intent(in)  ppoly0_coefs_R, real, dimension(nk,2), intent(in)  ppoly0_E_L, real, dimension(nk,2), intent(in)  ppoly0_E_R, integer, intent(in)  method, real, intent(in)  khtr_u, real, dimension(nk), intent(out)  F_layer, logical, intent(in), optional  linear_decay  )

private

Calculate the lateral boundary diffusive fluxes using the layer by layer method. See Along layer approach (Method #1).

Parameters

[in]	boundary	Which boundary layer SURFACE or BOTTOM [nondim]
[in]	nk	Number of layers [nondim]
[in]	deg	order of the polynomial reconstruction [nondim]
[in]	h_l	Layer thickness (left) [H ~> m or kg m-2]
[in]	h_r	Layer thickness (right) [H ~> m or kg m-2]
[in]	hbl_l	Thickness of the boundary boundary layer (left) [H ~> m or kg m-2]
[in]	hbl_r	Thickness of the boundary boundary layer (right) [H ~> m or kg m-2]
[in]	area_l	Area of the horizontal grid (left) [L2 ~> m2]
[in]	area_r	Area of the horizontal grid (right) [L2 ~> m2]
[in]	phi_l	Tracer values (left) [conc]
[in]	phi_r	Tracer values (right) [conc]
[in]	ppoly0_coefs_l	Tracer reconstruction (left) [conc]
[in]	ppoly0_coefs_r	Tracer reconstruction (right) [conc]
[in]	ppoly0_e_l	Polynomial edge values (left) [nondim]
[in]	ppoly0_e_r	Polynomial edge values (right) [nondim]
[in]	method	Method of polynomial integration [nondim]
[in]	khtr_u	Horizontal diffusivities times delta t at a velocity point [L2 ~> m2]
[out]	f_layer	Layerwise diffusive flux at U- or V-point [H L2 conc ~> m3 conc]
[in]	linear_decay	If True, apply a linear transition at the base of the boundary layer

Definition at line 444 of file MOM_lateral_boundary_diffusion.F90.

  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

harmonic_mean()

real function mom_lateral_boundary_diffusion::harmonic_mean ( real  h1, real  h2  )

private

Calculate the harmonic mean of two quantities See Harmonic Mean.

Parameters

h1	Scalar quantity
h2	Scalar quantity

Definition at line 364 of file MOM_lateral_boundary_diffusion.F90.

  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

lateral_boundary_diffusion()

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

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.

Parameters

[in,out]	g	Grid type
[in]	gv	ocean vertical grid structure
[in]	us	A dimensional unit scaling type
[in]	h	Layer thickness [H ~> m or kg m-2]
[in]	coef_x	dt * Kh * dy / dx at u-points [L2 ~> m2]
[in]	coef_y	dt * Kh * dx / dy at v-points [L2 ~> m2]
[in]	dt	Tracer time step * I_numitts (I_numitts in tracer_hordiff)
	reg	Tracer registry
[in]	cs	Control structure for this module

Definition at line 138 of file MOM_lateral_boundary_diffusion.F90.

  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

lateral_boundary_diffusion_init()

logical function, public mom_lateral_boundary_diffusion::lateral_boundary_diffusion_init	(	type(time_type), intent(in), target	Time,
		type(ocean_grid_type), intent(in)	G,
		type(param_file_type), intent(in)	param_file,
		type(diag_ctrl), intent(inout), target	diag,
		type(diabatic_cs), pointer	diabatic_CSp,
		type(lateral_boundary_diffusion_cs), pointer	CS
	)

Initialization routine that reads runtime parameters and sets up pointers to other control structures that might be needed for lateral boundary diffusion.

Parameters

[in]	time	Time structure
[in]	g	Grid structure
[in]	param_file	Parameter file structure
[in,out]	diag	Diagnostics control structure
	diabatic_csp	KPP control structure needed to get BLD
	cs	Lateral boundary mixing control structure

Definition at line 68 of file MOM_lateral_boundary_diffusion.F90.

  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)

near_boundary_unit_tests()

logical function, public mom_lateral_boundary_diffusion::near_boundary_unit_tests

(

logical, intent(in)

verbose

)

Unit tests for near-boundary horizontal mixing.

Parameters

[in]

verbose

If true, output additional information for debugging unit tests

Definition at line 807 of file MOM_lateral_boundary_diffusion.F90.

  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/) )

test_boundary_k_range()

logical function mom_lateral_boundary_diffusion::test_boundary_k_range ( integer  k_top, real  zeta_top, integer  k_bot, real  zeta_bot, integer  k_top_ans, real  zeta_top_ans, integer  k_bot_ans, real  zeta_bot_ans, character(len=80)  test_name, logical  verbose  )

private

Return true if output of unit tests for boundary_k_range does not match answers.

Parameters

k_top	Index of cell containing top of boundary
zeta_top	Nondimension position
k_bot	Index of cell containing bottom of boundary
zeta_bot	Nondimension position
k_top_ans	Index of cell containing top of boundary
zeta_top_ans	Nondimension position
k_bot_ans	Index of cell containing bottom of boundary
zeta_bot_ans	Nondimension position
test_name	Name of the unit test
verbose	If true always print output

Definition at line 1133 of file MOM_lateral_boundary_diffusion.F90.

  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)

test_layer_fluxes()

logical function mom_lateral_boundary_diffusion::test_layer_fluxes ( logical, intent(in)  verbose, integer, intent(in)  nk, character(len=80), intent(in)  test_name, real, dimension(nk), intent(in)  F_calc, real, dimension(nk), intent(in)  F_ans  )

private

Returns true if output of near-boundary unit tests does not match correct computed values and conditionally writes results to stream.

Parameters

[in]	verbose	If true, write results to stdout
[in]	test_name	Brief description of the unit test
[in]	nk	Number of layers
[in]	f_calc	Fluxes of the unitless tracer from the algorithm [s^-1]
[in]	f_ans	Fluxes of the unitless tracer calculated by hand [s^-1]

Definition at line 1107 of file MOM_lateral_boundary_diffusion.F90.

  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)