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

A brief summary of these methods is provided below.

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_{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: 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

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_{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. 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.

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)
	Calculate the lateral boundary diffusive fluxes using the layer by layer method. See Along layer approach (Method #2). 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)
	Apply the lateral boundary diffusive fluxes calculated from a 'bulk model' See Bulk layer approach (Method #1). 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. More...
integer, parameter, public	bottom = 1
	Set a value that corresponds to the bottom boundary. More...
character(len=40)	mdl = "MOM_lateral_boundary_diffusion"
	Name of this module. More...

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 [m]
[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)
[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 356 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         [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
 

References bottom, mom_error_handler::mom_error(), and surface.

Referenced by fluxes_bulk_method(), fluxes_layer_method(), near_boundary_unit_tests(), and mom_neutral_diffusion::neutral_diffusion_calc_coeffs().

Here is the call graph for this function:

Here is the caller graph for this function:

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(:,:)  ppoly0_E, real, dimension(:,:)  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 [m]
hblt	Depth of the boundary layer [m]
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 293 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                                         [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
 

References mom_remapping::average_value_ppoly(), bottom, mom_error_handler::mom_error(), and surface.

Referenced by fluxes_bulk_method().

Here is the call graph for this function:

Here is the caller graph for this function:

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, real, dimension(nk), intent(out), optional  F_limit  )

private

Apply the lateral boundary diffusive fluxes calculated from a 'bulk model' See Bulk 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) [m]
[in]	h_r	Layer thickness (right) [m]
[in]	hbl_l	Thickness of the boundary boundary layer (left) [m]
[in]	hbl_r	Thickness of the boundary boundary layer (left) [m]
[in]	area_l	Area of the horizontal grid (left) [m^2]
[in]	area_r	Area of the horizontal grid (right) [m^2]
[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 U-point [m^2]
[out]	f_bulk	The bulk mixed layer lateral flux [m^3 conc]
[out]	f_layer	Layerwise diffusive flux at U-point [m^3 conc]
[out]	f_limit	The amount of flux not applied due to limiter F_layer(k) - F_max [m^3 conc]

Definition at line 537 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)                  [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
 

References bottom, boundary_k_range(), bulk_average(), harmonic_mean(), and surface.

Referenced by lateral_boundary_diffusion(), and near_boundary_unit_tests().

Here is the call graph for this function:

Here is the caller graph for this function:

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  )

private

Calculate the lateral boundary diffusive fluxes using the layer by layer method. See Along 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) [m]
[in]	h_r	Layer thickness (right) [m]
[in]	hbl_l	Thickness of the boundary boundary layer (left) [m]
[in]	hbl_r	Thickness of the boundary boundary layer (right) [m]
[in]	area_l	Area of the horizontal grid (left) [m^2]
[in]	area_r	Area of the horizontal grid (right) [m^2]
[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 U-point [m^2]
[out]	f_layer	Layerwise diffusive flux at U- or V-point [m^3 conc]

Definition at line 424 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)                  [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
 

References mom_remapping::average_value_ppoly(), bottom, boundary_k_range(), harmonic_mean(), and surface.

Referenced by lateral_boundary_diffusion(), and near_boundary_unit_tests().

Here is the call graph for this function:

Here is the caller graph for this function:

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 345 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

Referenced by fluxes_bulk_method(), and fluxes_layer_method().

Here is the caller graph for this function:

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 [m2]
[in]	coef_y	dt * Kh * dx / dy at v-points [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 123 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 [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
 

References mom_energetic_pbl::energetic_pbl_get_mld(), fluxes_bulk_method(), fluxes_layer_method(), mom_cvmix_kpp::kpp_get_bld(), and surface.

Referenced by mom_tracer_hor_diff::tracer_hordiff().

Here is the call graph for this function:

Here is the caller graph for this function:

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

References mom_diabatic_driver::extract_diabatic_member(), mdl, mom_error_handler::mom_error(), mom_remapping::remappingdefaultscheme, and mom_remapping::remappingschemesdoc.

Here is the call graph for this function:

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

References bottom, boundary_k_range(), fluxes_bulk_method(), fluxes_layer_method(), surface, test_boundary_k_range(), and test_layer_fluxes().

Referenced by mom_unit_tests::unit_tests().

Here is the call graph for this function:

Here is the caller graph for this function:

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 1009 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 = 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)
 
 

Referenced by near_boundary_unit_tests().

Here is the caller graph for this function:

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 983 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 = 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)

Referenced by near_boundary_unit_tests().

Here is the caller graph for this function:

Variable Documentation

bottom

integer, parameter, public mom_lateral_boundary_diffusion::bottom = 1

Set a value that corresponds to the bottom boundary.

Definition at line 34 of file MOM_lateral_boundary_diffusion.F90.

integer, public, parameter :: BOTTOM  = 1  !< Set a value that corresponds to the bottom boundary

Referenced by boundary_k_range(), bulk_average(), fluxes_bulk_method(), fluxes_layer_method(), and near_boundary_unit_tests().

mdl

character(len=40) mom_lateral_boundary_diffusion::mdl = "MOM_lateral_boundary_diffusion"

Name of this module.

Definition at line 55 of file MOM_lateral_boundary_diffusion.F90.

character(len=40) :: mdl = "MOM_lateral_boundary_diffusion"       !< Name of this module

Referenced by lateral_boundary_diffusion_init().

surface

integer, parameter, public mom_lateral_boundary_diffusion::surface = -1

Set a value that corresponds to the surface bopundary.

Definition at line 33 of file MOM_lateral_boundary_diffusion.F90.

integer, public, parameter :: SURFACE = -1 !< Set a value that corresponds to the surface bopundary

Referenced by boundary_k_range(), bulk_average(), fluxes_bulk_method(), fluxes_layer_method(), lateral_boundary_diffusion(), near_boundary_unit_tests(), and mom_neutral_diffusion::neutral_diffusion_calc_coeffs().