MOM6
MOM_lateral_boundary_diffusion.F90
1 !> Calculates and applies diffusive fluxes as a parameterization of lateral mixing (non-neutral) by
2 !! mesoscale eddies near the top and bottom (to be implemented) boundary layers of the ocean.
3 
4 module mom_lateral_boundary_diffusion
5 
6 ! This file is part of MOM6. See LICENSE.md for the license.
7 
8 use mom_cpu_clock, only : cpu_clock_id, cpu_clock_begin, cpu_clock_end
9 use mom_cpu_clock, only : clock_module, clock_routine
10 use mom_domains, only : pass_var
11 use mom_diag_mediator, only : diag_ctrl, time_type
12 use mom_diag_mediator, only : post_data, register_diag_field
13 use mom_error_handler, only : mom_error, fatal, warning, mom_mesg, is_root_pe
14 use mom_file_parser, only : get_param, log_version, param_file_type
15 use mom_file_parser, only : openparameterblock, closeparameterblock
16 use mom_grid, only : ocean_grid_type
17 use mom_remapping, only : remapping_cs, initialize_remapping
18 use mom_remapping, only : extract_member_remapping_cs, build_reconstructions_1d
19 use mom_remapping, only : average_value_ppoly, remappingschemesdoc, remappingdefaultscheme
20 use mom_tracer_registry, only : tracer_registry_type, tracer_type
21 use mom_unit_scaling, only : unit_scale_type
22 use mom_verticalgrid, only : verticalgrid_type
23 use mom_cvmix_kpp, only : kpp_get_bld, kpp_cs
24 use mom_energetic_pbl, only : energetic_pbl_get_mld, energetic_pbl_cs
25 use mom_diabatic_driver, only : diabatic_cs, extract_diabatic_member
26 use iso_fortran_env, only : stdout=>output_unit, stderr=>error_unit
27 
28 implicit none ; private
29 
30 public near_boundary_unit_tests, lateral_boundary_diffusion, lateral_boundary_diffusion_init
31 public boundary_k_range
32 
33 ! Private parameters to avoid doing string comparisons for bottom or top boundary layer
34 integer, public, parameter :: surface = -1 !< Set a value that corresponds to the surface bopundary
35 integer, public, parameter :: bottom = 1 !< Set a value that corresponds to the bottom boundary
36 #include <MOM_memory.h>
37 
38 !> Sets parameters for lateral boundary mixing module.
39 type, public :: lateral_boundary_diffusion_cs ; private
40  integer :: method !< Determine which of the three methods calculate
41  !! and apply near boundary layer fluxes
42  !! 1. Along layer
43  !! 2. Bulk-layer approach (not recommended)
44  integer :: deg !< Degree of polynomial reconstruction
45  integer :: surface_boundary_scheme !< Which boundary layer scheme to use
46  !! 1. ePBL; 2. KPP
47  logical :: limiter !< Controls wether a flux limiter is applied.
48  !! Only valid when method = 2.
49  logical :: linear !< If True, apply a linear transition at the base/top of the boundary.
50  !! The flux will be fully applied at k=k_min and zero at k=k_max.
51 
52  type(remapping_cs) :: remap_cs !< Control structure to hold remapping configuration
53  type(kpp_cs), pointer :: kpp_csp => null() !< KPP control structure needed to get BLD
54  type(energetic_pbl_cs), pointer :: energetic_pbl_csp => null() !< ePBL control structure needed to get BLD
55  type(diag_ctrl), pointer :: diag => null() !< A structure that is used to
56  !! regulate the timing of diagnostic output.
57 end type lateral_boundary_diffusion_cs
58 
59 ! This include declares and sets the variable "version".
60 #include "version_variable.h"
61 character(len=40) :: mdl = "MOM_lateral_boundary_diffusion" !< Name of this module
62 
63 contains
64 
65 !> Initialization routine that reads runtime parameters and sets up pointers to other control structures that might be
66 !! needed for lateral boundary diffusion.
67 logical function lateral_boundary_diffusion_init(Time, G, param_file, diag, diabatic_CSp, CS)
68  type(time_type), target, intent(in) :: time !< Time structure
69  type(ocean_grid_type), intent(in) :: g !< Grid structure
70  type(param_file_type), intent(in) :: param_file !< Parameter file structure
71  type(diag_ctrl), target, intent(inout) :: diag !< Diagnostics control structure
72  type(diabatic_cs), pointer :: diabatic_csp !< KPP control structure needed to get BLD
73  type(lateral_boundary_diffusion_cs), pointer :: cs !< Lateral boundary mixing control structure
74 
75  ! local variables
76  character(len=80) :: string ! Temporary strings
77  logical :: boundary_extrap
78 
79  if (ASSOCIATED(cs)) then
80  call mom_error(fatal, "lateral_boundary_diffusion_init called with associated control structure.")
81  return
82  endif
83 
84  ! Log this module and master switch for turning it on/off
85  call get_param(param_file, mdl, "USE_LATERAL_BOUNDARY_DIFFUSION", lateral_boundary_diffusion_init, &
86  default=.false., do_not_log=.true.)
87  call log_version(param_file, mdl, version, &
88  "This module implements lateral diffusion of tracers near boundaries", &
89  all_default=.not.lateral_boundary_diffusion_init)
90  call get_param(param_file, mdl, "USE_LATERAL_BOUNDARY_DIFFUSION", lateral_boundary_diffusion_init, &
91  "If true, enables the lateral boundary tracer's diffusion module.", &
92  default=.false.)
93 
94  if (.not. lateral_boundary_diffusion_init) return
95 
96  allocate(cs)
97  cs%diag => diag
98  call extract_diabatic_member(diabatic_csp, kpp_csp=cs%KPP_CSp)
99  call extract_diabatic_member(diabatic_csp, energetic_pbl_csp=cs%energetic_PBL_CSp)
100 
101  cs%surface_boundary_scheme = -1
102  if ( .not. ASSOCIATED(cs%energetic_PBL_CSp) .and. .not. ASSOCIATED(cs%KPP_CSp) ) then
103  call mom_error(fatal,"Lateral boundary diffusion is true, but no valid boundary layer scheme was found")
104  endif
105 
106  ! Read all relevant parameters and write them to the model log.
107  call get_param(param_file, mdl, "LATERAL_BOUNDARY_METHOD", cs%method, &
108  "Determine how to apply boundary lateral diffusion of tracers: \n"//&
109  "1. Along layer approach \n"//&
110  "2. Bulk layer approach (this option is not recommended)", default=1)
111  if (cs%method == 2) then
112  call get_param(param_file, mdl, "APPLY_LIMITER", cs%limiter, &
113  "If True, apply a flux limiter in the LBD. This is only available \n"//&
114  "when LATERAL_BOUNDARY_METHOD=2.", default=.false.)
115  endif
116  call get_param(param_file, mdl, "LBD_LINEAR_TRANSITION", cs%linear, &
117  "If True, apply a linear transition at the base/top of the boundary. \n"//&
118  "The flux will be fully applied at k=k_min and zero at k=k_max.", default=.false.)
119  call get_param(param_file, mdl, "LBD_BOUNDARY_EXTRAP", boundary_extrap, &
120  "Use boundary extrapolation in LBD code", &
121  default=.false.)
122  call get_param(param_file, mdl, "LBD_REMAPPING_SCHEME", string, &
123  "This sets the reconstruction scheme used "//&
124  "for vertical remapping for all variables. "//&
125  "It can be one of the following schemes: "//&
126  trim(remappingschemesdoc), default=remappingdefaultscheme)
127  call initialize_remapping( cs%remap_CS, string, boundary_extrapolation = boundary_extrap )
128  call extract_member_remapping_cs(cs%remap_CS, degree=cs%deg)
129 
130 end function lateral_boundary_diffusion_init
131 
132 !> Driver routine for calculating lateral diffusive fluxes near the top and bottom boundaries.
133 !! Two different methods are available:
134 !! Method 1: lower order representation, calculate fluxes from bulk layer integrated quantities.
135 !! Method 2: more straight forward, diffusion is applied layer by layer using only information
136 !! from neighboring cells.
137 subroutine lateral_boundary_diffusion(G, GV, US, h, Coef_x, Coef_y, dt, Reg, CS)
138  type(ocean_grid_type), intent(inout) :: g !< Grid type
139  type(verticalgrid_type), intent(in) :: gv !< ocean vertical grid structure
140  type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
141  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
142  intent(in) :: h !< Layer thickness [H ~> m or kg m-2]
143  real, dimension(SZIB_(G),SZJ_(G)), intent(in) :: coef_x !< dt * Kh * dy / dx at u-points [L2 ~> m2]
144  real, dimension(SZI_(G),SZJB_(G)), intent(in) :: coef_y !< dt * Kh * dx / dy at v-points [L2 ~> m2]
145  real, intent(in) :: dt !< Tracer time step * I_numitts
146  !! (I_numitts in tracer_hordiff)
147  type(tracer_registry_type), pointer :: reg !< Tracer registry
148  type(lateral_boundary_diffusion_cs), intent(in) :: cs !< Control structure for this module
149 
150  ! Local variables
151  real, dimension(SZI_(G),SZJ_(G)) :: hbl !< bnd. layer depth [H ~> m or kg m-2]
152  real, dimension(SZI_(G),SZJ_(G),SZK_(G),CS%deg+1) :: ppoly0_coefs !< Coefficients of polynomial
153  real, dimension(SZI_(G),SZJ_(G),SZK_(G),2) :: ppoly0_e !< Edge values from reconstructions
154  real, dimension(SZK_(G),CS%deg+1) :: ppoly_s !< Slopes from reconstruction (placeholder)
155  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)) :: uflx !< Zonal flux of tracer [conc m^3]
156  real, dimension(SZIB_(G),SZJ_(G)) :: uflx_bulk !< Total calculated bulk-layer u-flux for the tracer
157  real, dimension(SZI_(G),SZJB_(G),SZK_(G)) :: vflx !< Meridional flux of tracer [conc m^3]
158  real, dimension(SZI_(G),SZJB_(G)) :: vflx_bulk !< Total calculated bulk-layer v-flux for the tracer
159  real, dimension(SZIB_(G),SZJ_(G)) :: uwork_2d !< Layer summed u-flux transport
160  real, dimension(SZI_(G),SZJB_(G)) :: vwork_2d !< Layer summed v-flux transport
161  real, dimension(SZI_(G),SZJ_(G),G%ke) :: tendency !< tendency array for diagn
162  real, dimension(SZI_(G),SZJ_(G)) :: tendency_2d !< depth integrated content tendency for diagn
163  type(tracer_type), pointer :: tracer => null() !< Pointer to the current tracer
164  integer :: remap_method !< Reconstruction method
165  integer :: i,j,k,m !< indices to loop over
166  real :: idt !< inverse of the time step [s-1]
167 
168  idt = 1./dt
169  hbl(:,:) = 0.
170  if (ASSOCIATED(cs%KPP_CSp)) call kpp_get_bld(cs%KPP_CSp, hbl, g, us, m_to_bld_units=gv%m_to_H)
171  if (ASSOCIATED(cs%energetic_PBL_CSp)) &
172  call energetic_pbl_get_mld(cs%energetic_PBL_CSp, hbl, g, us, m_to_mld_units=gv%m_to_H)
173 
174  call pass_var(hbl,g%Domain)
175  do m = 1,reg%ntr
176  tracer => reg%tr(m)
177 
178  ! for diagnostics
179  if (tracer%id_lbdxy_conc > 0 .or. tracer%id_lbdxy_cont > 0 .or. tracer%id_lbdxy_cont_2d > 0) then
180  tendency(:,:,:) = 0.0
181  endif
182 
183  ! Interpolate state to interface
184  do j=g%jsc-1,g%jec+1 ; do i=g%isc-1,g%iec+1
185  call build_reconstructions_1d( cs%remap_CS, g%ke, h(i,j,:), tracer%t(i,j,:), ppoly0_coefs(i,j,:,:), &
186  ppoly0_e(i,j,:,:), ppoly_s, remap_method, gv%H_subroundoff, gv%H_subroundoff)
187  enddo ; enddo
188 
189  ! Diffusive fluxes in the i-direction
190  uflx(:,:,:) = 0.
191  vflx(:,:,:) = 0.
192  uflx_bulk(:,:) = 0.
193  vflx_bulk(:,:) = 0.
194 
195  ! Method #1 (layer by layer)
196  if (cs%method == 1) then
197  do j=g%jsc,g%jec
198  do i=g%isc-1,g%iec
199  if (g%mask2dCu(i,j)>0.) then
200  call fluxes_layer_method(surface, gv%ke, cs%deg, h(i,j,:), h(i+1,j,:), hbl(i,j), hbl(i+1,j), &
201  g%areaT(i,j), g%areaT(i+1,j), tracer%t(i,j,:), tracer%t(i+1,j,:), ppoly0_coefs(i,j,:,:), &
202  ppoly0_coefs(i+1,j,:,:), ppoly0_e(i,j,:,:), ppoly0_e(i+1,j,:,:), remap_method, coef_x(i,j), &
203  uflx(i,j,:), cs%linear)
204  endif
205  enddo
206  enddo
207  do j=g%jsc-1,g%jec
208  do i=g%isc,g%iec
209  if (g%mask2dCv(i,j)>0.) then
210  call fluxes_layer_method(surface, gv%ke, cs%deg, h(i,j,:), h(i,j+1,:), hbl(i,j), hbl(i,j+1), &
211  g%areaT(i,j), g%areaT(i,j+1), tracer%t(i,j,:), tracer%t(i,j+1,:), ppoly0_coefs(i,j,:,:), &
212  ppoly0_coefs(i,j+1,:,:), ppoly0_e(i,j,:,:), ppoly0_e(i,j+1,:,:), remap_method, coef_y(i,j), &
213  vflx(i,j,:), cs%linear)
214  endif
215  enddo
216  enddo
217 
218  ! Method #2 (bulk approach)
219  elseif (cs%method == 2) then
220  do j=g%jsc,g%jec
221  do i=g%isc-1,g%iec
222  if (g%mask2dCu(i,j)>0.) then
223  call fluxes_bulk_method(surface, gv%ke, cs%deg, h(i,j,:), h(i+1,j,:), hbl(i,j), hbl(i+1,j), &
224  g%areaT(i,j), g%areaT(i+1,j), tracer%t(i,j,:), tracer%t(i+1,j,:), &
225  ppoly0_coefs(i,j,:,:), ppoly0_coefs(i+1,j,:,:), ppoly0_e(i,j,:,:), &
226  ppoly0_e(i+1,j,:,:), remap_method, coef_x(i,j), uflx_bulk(i,j), uflx(i,j,:), cs%limiter, &
227  cs%linear)
228  endif
229  enddo
230  enddo
231  do j=g%jsc-1,g%jec
232  do i=g%isc,g%iec
233  if (g%mask2dCv(i,j)>0.) then
234  call fluxes_bulk_method(surface, gv%ke, cs%deg, h(i,j,:), h(i,j+1,:), hbl(i,j), hbl(i,j+1), &
235  g%areaT(i,j), g%areaT(i,j+1), tracer%t(i,j,:), tracer%t(i,j+1,:), &
236  ppoly0_coefs(i,j,:,:), ppoly0_coefs(i,j+1,:,:), ppoly0_e(i,j,:,:), &
237  ppoly0_e(i,j+1,:,:), remap_method, coef_y(i,j), vflx_bulk(i,j), vflx(i,j,:), cs%limiter, &
238  cs%linear)
239  endif
240  enddo
241  enddo
242  ! Post tracer bulk diags
243  if (tracer%id_lbd_bulk_dfx>0) call post_data(tracer%id_lbd_bulk_dfx, uflx_bulk*idt, cs%diag)
244  if (tracer%id_lbd_bulk_dfy>0) call post_data(tracer%id_lbd_bulk_dfy, vflx_bulk*idt, cs%diag)
245  endif
246 
247  ! Update the tracer fluxes
248  do k=1,gv%ke ; do j=g%jsc,g%jec ; do i=g%isc,g%iec
249  if (g%mask2dT(i,j)>0.) then
250  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) ) ))* &
251  (g%IareaT(i,j)/( h(i,j,k) + gv%H_subroundoff))
252 
253  if (tracer%id_lbdxy_conc > 0 .or. tracer%id_lbdxy_cont > 0 .or. tracer%id_lbdxy_cont_2d > 0 ) then
254  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
255  endif
256 
257  endif
258  enddo ; enddo ; enddo
259 
260  ! Post the tracer diagnostics
261  if (tracer%id_lbd_dfx>0) call post_data(tracer%id_lbd_dfx, uflx*idt, cs%diag)
262  if (tracer%id_lbd_dfy>0) call post_data(tracer%id_lbd_dfy, vflx*idt, cs%diag)
263  if (tracer%id_lbd_dfx_2d>0) then
264  uwork_2d(:,:) = 0.
265  do k=1,gv%ke ; do j=g%jsc,g%jec ; do i=g%isc-1,g%iec
266  uwork_2d(i,j) = uwork_2d(i,j) + (uflx(i,j,k) * idt)
267  enddo ; enddo ; enddo
268  call post_data(tracer%id_lbd_dfx_2d, uwork_2d, cs%diag)
269  endif
270 
271  if (tracer%id_lbd_dfy_2d>0) then
272  vwork_2d(:,:) = 0.
273  do k=1,gv%ke ; do j=g%jsc-1,g%jec ; do i=g%isc,g%iec
274  vwork_2d(i,j) = vwork_2d(i,j) + (vflx(i,j,k) * idt)
275  enddo ; enddo ; enddo
276  call post_data(tracer%id_lbd_dfy_2d, vwork_2d, cs%diag)
277  endif
278 
279  ! post tendency of tracer content
280  if (tracer%id_lbdxy_cont > 0) then
281  call post_data(tracer%id_lbdxy_cont, tendency, cs%diag)
282  endif
283 
284  ! post depth summed tendency for tracer content
285  if (tracer%id_lbdxy_cont_2d > 0) then
286  tendency_2d(:,:) = 0.
287  do j=g%jsc,g%jec ; do i=g%isc,g%iec
288  do k=1,gv%ke
289  tendency_2d(i,j) = tendency_2d(i,j) + tendency(i,j,k)
290  enddo
291  enddo ; enddo
292  call post_data(tracer%id_lbdxy_cont_2d, tendency_2d, cs%diag)
293  endif
294 
295  ! post tendency of tracer concentration; this step must be
296  ! done after posting tracer content tendency, since we alter
297  ! the tendency array and its units.
298  if (tracer%id_lbdxy_conc > 0) then
299  do k=1,gv%ke ; do j=g%jsc,g%jec ; do i=g%isc,g%iec
300  tendency(i,j,k) = tendency(i,j,k) / ( h(i,j,k) + gv%H_subroundoff )
301  enddo ; enddo ; enddo
302  call post_data(tracer%id_lbdxy_conc, tendency, cs%diag)
303  endif
304 
305  enddo
306 
307 end subroutine lateral_boundary_diffusion
308 
309 !< Calculate bulk layer value of a scalar quantity as the thickness weighted average
310 real function bulk_average(boundary, nk, deg, h, hBLT, phi, ppoly0_E, ppoly0_coefs, method, k_top, zeta_top, k_bot, &
311  zeta_bot)
312  integer :: boundary !< SURFACE or BOTTOM [nondim]
313  integer :: nk !< Number of layers [nondim]
314  integer :: deg !< Degree of polynomial [nondim]
315  real, dimension(nk) :: h !< Layer thicknesses [H ~> m or kg m-2]
316  real :: hblt !< Depth of the boundary layer [H ~> m or kg m-2]
317  real, dimension(nk) :: phi !< Scalar quantity
318  real, dimension(nk,2) :: ppoly0_e !< Edge value of polynomial
319  real, dimension(nk,deg+1) :: ppoly0_coefs!< Coefficients of polynomial
320  integer :: method !< Remapping scheme to use
321 
322  integer :: k_top !< Index of the first layer within the boundary
323  real :: zeta_top !< Fraction of the layer encompassed by the bottom boundary layer
324  !! (0 if none, 1. if all). For the surface, this is always 0. because
325  !! integration starts at the surface [nondim]
326  integer :: k_bot !< Index of the last layer within the boundary
327  real :: zeta_bot !< Fraction of the layer encompassed by the surface boundary layer
328  !! (0 if none, 1. if all). For the bottom boundary layer, this is always 1.
329  !! because integration starts at the bottom [nondim]
330  ! Local variables
331  real :: htot !< Running sum of the thicknesses (top to bottom)
332  integer :: k !< k indice
333 
334 
335  htot = 0.
336  bulk_average = 0.
337  if (hblt == 0.) return
338  if (boundary == surface) then
339  htot = (h(k_bot) * zeta_bot)
340  bulk_average = average_value_ppoly( nk, phi, ppoly0_e, ppoly0_coefs, method, k_bot, 0., zeta_bot) * htot
341  do k = k_bot-1,1,-1
342  bulk_average = bulk_average + phi(k)*h(k)
343  htot = htot + h(k)
344  enddo
345  elseif (boundary == bottom) then
346  htot = (h(k_top) * zeta_top)
347  ! (note 1-zeta_top because zeta_top is the fraction of the layer)
348  bulk_average = average_value_ppoly( nk, phi, ppoly0_e, ppoly0_coefs, method, k_top, (1.-zeta_top), 1.) * htot
349  do k = k_top+1,nk
350  bulk_average = bulk_average + phi(k)*h(k)
351  htot = htot + h(k)
352  enddo
353  else
354  call mom_error(fatal, "bulk_average: a valid boundary type must be provided.")
355  endif
356 
357  bulk_average = bulk_average / hblt
358 
359 end function bulk_average
360 
361 !> Calculate the harmonic mean of two quantities
362 !! See \ref section_harmonic_mean.
363 real function harmonic_mean(h1,h2)
364  real :: h1 !< Scalar quantity
365  real :: h2 !< Scalar quantity
366  if (h1 + h2 == 0.) then
367  harmonic_mean = 0.
368  else
369  harmonic_mean = 2.*(h1*h2)/(h1+h2)
370  endif
371 end function harmonic_mean
372 
373 !> Find the k-index range corresponding to the layers that are within the boundary-layer region
374 subroutine boundary_k_range(boundary, nk, h, hbl, k_top, zeta_top, k_bot, zeta_bot)
375  integer, intent(in ) :: boundary !< SURFACE or BOTTOM [nondim]
376  integer, intent(in ) :: nk !< Number of layers [nondim]
377  real, dimension(nk), intent(in ) :: h !< Layer thicknesses of the column [H ~> m or kg m-2]
378  real, intent(in ) :: hbl !< Thickness of the boundary layer [H ~> m or kg m-2]
379  !! If surface, with respect to zbl_ref = 0.
380  !! If bottom, with respect to zbl_ref = SUM(h)
381  integer, intent( out) :: k_top !< Index of the first layer within the boundary
382  real, intent( out) :: zeta_top !< Distance from the top of a layer to the intersection of the
383  !! top extent of the boundary layer (0 at top, 1 at bottom) [nondim]
384  integer, intent( out) :: k_bot !< Index of the last layer within the boundary
385  real, intent( out) :: zeta_bot !< Distance of the lower layer to the boundary layer depth
386  !! (0 at top, 1 at bottom) [nondim]
387  ! Local variables
388  real :: htot ! Summed thickness [H ~> m or kg m-2]
389  integer :: k
390  ! Surface boundary layer
391  if ( boundary == surface ) then
392  k_top = 1
393  zeta_top = 0.
394  htot = 0.
395  k_bot = 1
396  zeta_bot = 0.
397  if (hbl == 0.) return
398  if (hbl >= sum(h(:))) then
399  k_bot = nk
400  zeta_bot = 1.
401  return
402  endif
403  do k=1,nk
404  htot = htot + h(k)
405  if ( htot >= hbl) then
406  k_bot = k
407  zeta_bot = 1 - (htot - hbl)/h(k)
408  return
409  endif
410  enddo
411  ! Bottom boundary layer
412  elseif ( boundary == bottom ) then
413  k_top = nk
414  zeta_top = 1.
415  k_bot = nk
416  zeta_bot = 0.
417  htot = 0.
418  if (hbl == 0.) return
419  if (hbl >= sum(h(:))) then
420  k_top = 1
421  zeta_top = 1.
422  return
423  endif
424  do k=nk,1,-1
425  htot = htot + h(k)
426  if (htot >= hbl) then
427  k_top = k
428  zeta_top = 1 - (htot - hbl)/h(k)
429  return
430  endif
431  enddo
432  else
433  call mom_error(fatal,"Houston, we've had a problem in boundary_k_range")
434  endif
435 
436 end subroutine boundary_k_range
437 
438 
439 !> Calculate the lateral boundary diffusive fluxes using the layer by layer method.
440 !! See \ref section_method1
441 subroutine fluxes_layer_method(boundary, nk, deg, h_L, h_R, hbl_L, hbl_R, area_L, area_R, phi_L, phi_R, &
442  ppoly0_coefs_L, ppoly0_coefs_R, ppoly0_E_L, ppoly0_E_R, method, khtr_u, &
443  F_layer, linear_decay)
444 
445  integer, intent(in ) :: boundary !< Which boundary layer SURFACE or BOTTOM [nondim]
446  integer, intent(in ) :: nk !< Number of layers [nondim]
447  integer, intent(in ) :: deg !< order of the polynomial reconstruction [nondim]
448  real, dimension(nk), intent(in ) :: h_L !< Layer thickness (left) [H ~> m or kg m-2]
449  real, dimension(nk), intent(in ) :: h_R !< Layer thickness (right) [H ~> m or kg m-2]
450  real, intent(in ) :: hbl_L !< Thickness of the boundary boundary
451  !! layer (left) [H ~> m or kg m-2]
452  real, intent(in ) :: hbl_R !< Thickness of the boundary boundary
453  !! layer (right) [H ~> m or kg m-2]
454  real, intent(in ) :: area_L !< Area of the horizontal grid (left) [L2 ~> m2]
455  real, intent(in ) :: area_R !< Area of the horizontal grid (right) [L2 ~> m2]
456  real, dimension(nk), intent(in ) :: phi_L !< Tracer values (left) [conc]
457  real, dimension(nk), intent(in ) :: phi_R !< Tracer values (right) [conc]
458  real, dimension(nk,deg+1), intent(in ) :: ppoly0_coefs_L !< Tracer reconstruction (left) [conc]
459  real, dimension(nk,deg+1), intent(in ) :: ppoly0_coefs_R !< Tracer reconstruction (right) [conc]
460  real, dimension(nk,2), intent(in ) :: ppoly0_E_L !< Polynomial edge values (left) [nondim]
461  real, dimension(nk,2), intent(in ) :: ppoly0_E_R !< Polynomial edge values (right) [nondim]
462  integer, intent(in ) :: method !< Method of polynomial integration [nondim]
463  real, intent(in ) :: khtr_u !< Horizontal diffusivities times delta t
464  !! at a velocity point [L2 ~> m2]
465  real, dimension(nk), intent( out) :: F_layer !< Layerwise diffusive flux at U- or V-point
466  !! [H L2 conc ~> m3 conc]
467  logical, optional, intent(in ) :: linear_decay !< If True, apply a linear transition at the base of
468  !! the boundary layer
469  ! Local variables
470  real, dimension(nk) :: h_means !< Calculate the layer-wise harmonic means [H ~> m or kg m-2]
471  real :: khtr_avg !< Thickness-weighted diffusivity at the u-point [m^2 s^-1]
472  !! This is just to remind developers that khtr_avg should be
473  !! computed once khtr is 3D.
474  real :: heff !< Harmonic mean of layer thicknesses [H ~> m or kg m-2]
475  real :: inv_heff !< Inverse of the harmonic mean of layer thicknesses
476  !! [H-1 ~> m-1 or m2 kg-1]
477  real :: phi_L_avg, phi_R_avg !< Bulk, thickness-weighted tracer averages (left and right column)
478  !! [conc m^-3 ]
479  real :: htot !< Total column thickness [H ~> m or kg m-2]
480  real :: heff_tot !< Total effective column thickness in the transition layer [m]
481  integer :: k, k_bot_min, k_top_max !< k-indices, min and max for bottom and top, respectively
482  integer :: k_bot_max, k_top_min !< k-indices, max and min for bottom and top, respectively
483  integer :: k_bot_diff, k_top_diff !< different between left and right k-indices for bottom and top, respectively
484  integer :: k_top_L, k_bot_L !< k-indices left
485  integer :: k_top_R, k_bot_R !< k-indices right
486  real :: zeta_top_L, zeta_top_R !< distance from the top of a layer to the boundary
487  !! layer depth [nondim]
488  real :: zeta_bot_L, zeta_bot_R !< distance from the bottom of a layer to the boundary
489  !!layer depth [nondim]
490  real :: h_work_L, h_work_R !< dummy variables
491  real :: hbl_min !< minimum BLD (left and right) [m]
492  real :: wgt !< weight to be used in the linear transition to the interior [nondim]
493  real :: a !< coefficient to be used in the linear transition to the interior [nondim]
494  logical :: linear !< True if apply a linear transition
495 
496  f_layer(:) = 0.0
497  if (hbl_l == 0. .or. hbl_r == 0.) then
498  return
499  endif
500 
501  linear = .false.
502  if (PRESENT(linear_decay)) then
503  linear = linear_decay
504  endif
505 
506  ! Calculate vertical indices containing the boundary layer
507  call boundary_k_range(boundary, nk, h_l, hbl_l, k_top_l, zeta_top_l, k_bot_l, zeta_bot_l)
508  call boundary_k_range(boundary, nk, h_r, hbl_r, k_top_r, zeta_top_r, k_bot_r, zeta_bot_r)
509 
510  if (boundary == surface) then
511  k_bot_min = min(k_bot_l, k_bot_r)
512  k_bot_max = max(k_bot_l, k_bot_r)
513  k_bot_diff = (k_bot_max - k_bot_min)
514 
515  ! make sure left and right k indices span same range
516  if (k_bot_min .ne. k_bot_l) then
517  k_bot_l = k_bot_min
518  zeta_bot_l = 1.0
519  endif
520  if (k_bot_min .ne. k_bot_r) then
521  k_bot_r= k_bot_min
522  zeta_bot_r = 1.0
523  endif
524 
525  h_work_l = (h_l(k_bot_l) * zeta_bot_l)
526  h_work_r = (h_r(k_bot_r) * zeta_bot_r)
527 
528  phi_l_avg = average_value_ppoly( nk, phi_l, ppoly0_e_l, ppoly0_coefs_l, method, k_bot_l, 0., zeta_bot_l)
529  phi_r_avg = average_value_ppoly( nk, phi_r, ppoly0_e_r, ppoly0_coefs_r, method, k_bot_r, 0., zeta_bot_r)
530  heff = harmonic_mean(h_work_l, h_work_r)
531  ! tracer flux where the minimum BLD intersets layer
532  ! GMM, khtr_avg should be computed once khtr is 3D
533  if ((linear) .and. (k_bot_diff .gt. 1)) then
534  ! apply linear decay at the base of hbl
535  do k = k_bot_min,1,-1
536  heff = harmonic_mean(h_l(k), h_r(k))
537  f_layer(k) = -(heff * khtr_u) * (phi_r(k) - phi_l(k))
538  enddo
539  ! heff_total
540  heff_tot = 0.0
541  do k = k_bot_min+1,k_bot_max, 1
542  heff_tot = heff_tot + harmonic_mean(h_l(k), h_r(k))
543  enddo
544 
545  a = -1.0/heff_tot
546  heff_tot = 0.0
547  do k = k_bot_min+1,k_bot_max, 1
548  heff = harmonic_mean(h_l(k), h_r(k))
549  wgt = (a*(heff_tot + (heff * 0.5))) + 1.0
550  f_layer(k) = -(heff * khtr_u) * (phi_r(k) - phi_l(k)) * wgt
551  heff_tot = heff_tot + heff
552  enddo
553  else
554  f_layer(k_bot_min) = -(heff * khtr_u) * (phi_r_avg - phi_l_avg)
555  do k = k_bot_min-1,1,-1
556  heff = harmonic_mean(h_l(k), h_r(k))
557  f_layer(k) = -(heff * khtr_u) * (phi_r(k) - phi_l(k))
558  enddo
559  endif
560  endif
561 
562  if (boundary == bottom) then
563  ! TODO: GMM add option to apply linear decay
564  k_top_max = max(k_top_l, k_top_r)
565  ! make sure left and right k indices span same range
566  if (k_top_max .ne. k_top_l) then
567  k_top_l = k_top_max
568  zeta_top_l = 1.0
569  endif
570  if (k_top_max .ne. k_top_r) then
571  k_top_r= k_top_max
572  zeta_top_r = 1.0
573  endif
574 
575  h_work_l = (h_l(k_top_l) * zeta_top_l)
576  h_work_r = (h_r(k_top_r) * zeta_top_r)
577 
578  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)
579  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)
580  heff = harmonic_mean(h_work_l, h_work_r)
581 
582  ! tracer flux where the minimum BLD intersets layer
583  f_layer(k_top_max) = (-heff * khtr_u) * (phi_r_avg - phi_l_avg)
584 
585  do k = k_top_max+1,nk
586  heff = harmonic_mean(h_l(k), h_r(k))
587  f_layer(k) = -(heff * khtr_u) * (phi_r(k) - phi_l(k))
588  enddo
589  endif
590 
591 end subroutine fluxes_layer_method
592 
593 !> Apply the lateral boundary diffusive fluxes calculated from a 'bulk model'
594 !! See \ref section_method2
595 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, &
596  ppoly0_coefs_R, ppoly0_E_L, ppoly0_E_R, method, khtr_u, F_bulk, F_layer, F_limit, &
597  linear_decay)
598 
599  integer, intent(in ) :: boundary !< Which boundary layer SURFACE or BOTTOM [nondim]
600  integer, intent(in ) :: nk !< Number of layers [nondim]
601  integer, intent(in ) :: deg !< order of the polynomial reconstruction [nondim]
602  real, dimension(nk), intent(in ) :: h_L !< Layer thickness (left) [H ~> m or kg m-2]
603  real, dimension(nk), intent(in ) :: h_R !< Layer thickness (right) [H ~> m or kg m-2]
604  real, intent(in ) :: hbl_L !< Thickness of the boundary boundary
605  !! layer (left) [H ~> m or kg m-2]
606  real, intent(in ) :: hbl_R !< Thickness of the boundary boundary
607  !! layer (left) [H ~> m or kg m-2]
608  real, intent(in ) :: area_L !< Area of the horizontal grid (left) [L2 ~> m2]
609  real, intent(in ) :: area_R !< Area of the horizontal grid (right) [L2 ~> m2]
610  real, dimension(nk), intent(in ) :: phi_L !< Tracer values (left) [conc]
611  real, dimension(nk), intent(in ) :: phi_R !< Tracer values (right) [conc]
612  real, dimension(nk,deg+1), intent(in ) :: ppoly0_coefs_L !< Tracer reconstruction (left) [conc]
613  real, dimension(nk,deg+1), intent(in ) :: ppoly0_coefs_R !< Tracer reconstruction (right) [conc]
614  real, dimension(nk,2), intent(in ) :: ppoly0_E_L !< Polynomial edge values (left) [nondim]
615  real, dimension(nk,2), intent(in ) :: ppoly0_E_R !< Polynomial edge values (right) [nondim]
616  integer, intent(in ) :: method !< Method of polynomial integration [nondim]
617  real, intent(in ) :: khtr_u !< Horizontal diffusivities times delta t
618  !! at a velocity point [L2 ~> m2]
619  real, intent( out) :: F_bulk !< The bulk mixed layer lateral flux
620  !! [H L2 conc ~> m3 conc]
621  real, dimension(nk), intent( out) :: F_layer !< Layerwise diffusive flux at U- or V-point
622  !! [H L2 conc ~> m3 conc]
623  logical, optional, intent(in ) :: F_limit !< If True, apply a limiter
624  logical, optional, intent(in ) :: linear_decay !< If True, apply a linear transition at the base of
625  !! the boundary layer
626 
627  ! Local variables
628  real, dimension(nk) :: h_means !< Calculate the layer-wise harmonic means [H ~> m or kg m-2]
629  real :: khtr_avg !< Thickness-weighted diffusivity at the u-point [m^2 s^-1]
630  !! This is just to remind developers that khtr_avg should be
631  !! computed once khtr is 3D.
632  real :: heff !< Harmonic mean of layer thicknesses [H ~> m or kg m-2]
633  real :: heff_tot !< Total effective column thickness in the transition layer [m]
634  real :: inv_heff !< Inverse of the harmonic mean of layer thicknesses
635  !! [H-1 ~> m-1 or m2 kg-1]
636  real :: phi_L_avg, phi_R_avg !< Bulk, thickness-weighted tracer averages (left and right column)
637  !! [conc m^-3 ]
638  real :: htot ! Total column thickness [H ~> m or kg m-2]
639  integer :: k, k_min, k_max !< k-indices, min and max for top and bottom, respectively
640  integer :: k_diff !< difference between k_max and k_min
641  integer :: k_top_L, k_bot_L !< k-indices left
642  integer :: k_top_R, k_bot_R !< k-indices right
643  real :: zeta_top_L, zeta_top_R !< distance from the top of a layer to the
644  !! boundary layer [nondim]
645  real :: zeta_bot_L, zeta_bot_R !< distance from the bottom of a layer to the
646  !! boundary layer [nondim]
647  real :: h_work_L, h_work_R !< dummy variables
648  real :: F_max !< The maximum amount of flux that can leave a
649  !! cell [m^3 conc]
650  logical :: limiter !< True if flux limiter should be applied
651  logical :: linear !< True if apply a linear transition
652  real :: hfrac !< Layer fraction wrt sum of all layers [nondim]
653  real :: dphi !< tracer gradient [conc m^-3]
654  real :: wgt !< weight to be used in the linear transition to the
655  !! interior [nondim]
656  real :: a !< coefficient to be used in the linear transition to the
657  !! interior [nondim]
658 
659  f_bulk = 0.
660  f_layer(:) = 0.
661  if (hbl_l == 0. .or. hbl_r == 0.) then
662  return
663  endif
664 
665  limiter = .false.
666  if (PRESENT(f_limit)) then
667  limiter = f_limit
668  endif
669  linear = .false.
670  if (PRESENT(linear_decay)) then
671  linear = linear_decay
672  endif
673 
674  ! Calculate vertical indices containing the boundary layer
675  call boundary_k_range(boundary, nk, h_l, hbl_l, k_top_l, zeta_top_l, k_bot_l, zeta_bot_l)
676  call boundary_k_range(boundary, nk, h_r, hbl_r, k_top_r, zeta_top_r, k_bot_r, zeta_bot_r)
677 
678  ! Calculate bulk averages of various quantities
679  phi_l_avg = bulk_average(boundary, nk, deg, h_l, hbl_l, phi_l, ppoly0_e_l, ppoly0_coefs_l, method, k_top_l, &
680  zeta_top_l, k_bot_l, zeta_bot_l)
681  phi_r_avg = bulk_average(boundary, nk, deg, h_r, hbl_r, phi_r, ppoly0_e_r, ppoly0_coefs_r, method, k_top_r, &
682  zeta_top_r, k_bot_r, zeta_bot_r)
683  ! Calculate the 'bulk' diffusive flux from the bulk averaged quantities
684  ! GMM, khtr_avg should be computed once khtr is 3D
685  heff = harmonic_mean(hbl_l, hbl_r)
686  f_bulk = -(khtr_u * heff) * (phi_r_avg - phi_l_avg)
687  ! Calculate the layerwise sum of the vertical effective thickness. This is different than the heff calculated
688  ! above, but is used as a way to decompose the fluxes onto the individual layers
689  h_means(:) = 0.
690  if (boundary == surface) then
691  k_min = min(k_bot_l, k_bot_r)
692  k_max = max(k_bot_l, k_bot_r)
693  k_diff = (k_max - k_min)
694  if ((linear) .and. (k_diff .gt. 1)) then
695  do k=1,k_min
696  h_means(k) = harmonic_mean(h_l(k),h_r(k))
697  enddo
698  ! heff_total
699  heff_tot = 0.0
700  do k = k_min+1,k_max, 1
701  heff_tot = heff_tot + harmonic_mean(h_l(k), h_r(k))
702  enddo
703 
704  a = -1.0/heff_tot
705  heff_tot = 0.0
706  ! fluxes will decay linearly at base of hbl
707  do k = k_min+1,k_max, 1
708  heff = harmonic_mean(h_l(k), h_r(k))
709  wgt = (a*(heff_tot + (heff * 0.5))) + 1.0
710  h_means(k) = harmonic_mean(h_l(k), h_r(k)) * wgt
711  heff_tot = heff_tot + heff
712  enddo
713  else
714  ! left hand side
715  if (k_bot_l == k_min) then
716  h_work_l = h_l(k_min) * zeta_bot_l
717  else
718  h_work_l = h_l(k_min)
719  endif
720 
721  ! right hand side
722  if (k_bot_r == k_min) then
723  h_work_r = h_r(k_min) * zeta_bot_r
724  else
725  h_work_r = h_r(k_min)
726  endif
727 
728  h_means(k_min) = harmonic_mean(h_work_l,h_work_r)
729 
730  do k=1,k_min-1
731  h_means(k) = harmonic_mean(h_l(k),h_r(k))
732  enddo
733  endif
734 
735  elseif (boundary == bottom) then
736  !TODO, GMM linear decay is not implemented here
737  k_max = max(k_top_l, k_top_r)
738  ! left hand side
739  if (k_top_l == k_max) then
740  h_work_l = h_l(k_max) * zeta_top_l
741  else
742  h_work_l = h_l(k_max)
743  endif
744 
745  ! right hand side
746  if (k_top_r == k_max) then
747  h_work_r = h_r(k_max) * zeta_top_r
748  else
749  h_work_r = h_r(k_max)
750  endif
751 
752  h_means(k_max) = harmonic_mean(h_work_l,h_work_r)
753 
754  do k=nk,k_max+1,-1
755  h_means(k) = harmonic_mean(h_l(k),h_r(k))
756  enddo
757  endif
758 
759  if ( sum(h_means) == 0. .or. f_bulk == 0.) then
760  return
761  ! Decompose the bulk flux onto the individual layers
762  else
763  ! Initialize remaining thickness
764  inv_heff = 1./sum(h_means)
765  do k=1,nk
766  if ((h_means(k) > 0.) .and. (phi_l(k) /= phi_r(k))) then
767  hfrac = h_means(k)*inv_heff
768  f_layer(k) = f_bulk * hfrac
769 
770  if (limiter) then
771  ! limit the flux to 0.2 of the tracer *gradient*
772  ! Why 0.2?
773  ! t=0 t=inf
774  ! 0 .2
775  ! 0 1 0 .2.2.2
776  ! 0 .2
777  !
778  f_max = -0.2 * ((area_r*(phi_r(k)*h_r(k)))-(area_l*(phi_l(k)*h_r(k))))
779 
780  ! check if bulk flux (or F_layer) and F_max have same direction
781  if ( sign(1.,f_bulk) == sign(1., f_max)) then
782  ! Apply flux limiter calculated above
783  if (f_max >= 0.) then
784  f_layer(k) = min(f_layer(k),f_max)
785  else
786  f_layer(k) = max(f_layer(k),f_max)
787  endif
788  else
789  ! do not apply a flux on this layer
790  f_layer(k) = 0.
791  endif
792  else
793  dphi = -(phi_r(k) - phi_l(k))
794  if (.not. sign(1.,f_bulk) == sign(1., dphi)) then
795  ! upgradient, do not apply a flux on this layer
796  f_layer(k) = 0.
797  endif
798  endif ! limited
799  endif
800  enddo
801  endif
802 
803 end subroutine fluxes_bulk_method
804 
805 !> Unit tests for near-boundary horizontal mixing
806 logical function near_boundary_unit_tests( verbose )
807  logical, intent(in) :: verbose !< If true, output additional information for debugging unit tests
808 
809  ! Local variables
810  integer, parameter :: nk = 2 ! Number of layers
811  integer, parameter :: deg = 1 ! Degree of reconstruction (linear here)
812  integer, parameter :: method = 1 ! Method used for integrating polynomials
813  real, dimension(nk) :: phi_l, phi_r ! Tracer values (left and right column) [ nondim m^-3 ]
814  real, dimension(nk) :: phi_l_avg, phi_r_avg ! Bulk, thickness-weighted tracer averages (left and right column)
815  real, dimension(nk,deg+1) :: phi_pp_l, phi_pp_r ! Coefficients for the linear pseudo-reconstructions
816  ! [ nondim m^-3 ]
817 
818  real, dimension(nk,2) :: ppoly0_e_l, ppoly0_e_r! Polynomial edge values (left and right) [concentration]
819  real, dimension(nk) :: h_l, h_r ! Layer thickness (left and right) [m]
820  real :: khtr_u ! Horizontal diffusivities at U-point [m^2 s^-1]
821  real :: hbl_l, hbl_r ! Depth of the boundary layer (left and right) [m]
822  real :: f_bulk ! Total diffusive flux across the U point [nondim s^-1]
823  real, dimension(nk) :: f_layer ! Diffusive flux within each layer at U-point [nondim s^-1]
824  real :: h_u, hblt_u ! Thickness at the u-point [m]
825  real :: khtr_avg ! Thickness-weighted diffusivity at the u-point [m^2 s^-1]
826  real :: heff ! Harmonic mean of layer thicknesses [m]
827  real :: inv_heff ! Inverse of the harmonic mean of layer thicknesses [m^[-1]
828  character(len=120) :: test_name ! Title of the unit test
829  integer :: k_top ! Index of cell containing top of boundary
830  real :: zeta_top ! Nondimension position
831  integer :: k_bot ! Index of cell containing bottom of boundary
832  real :: zeta_bot ! Nondimension position
833  real :: area_l,area_r ! Area of grid cell [m^2]
834  area_l = 1.; area_r = 1. ! Set to unity for all unit tests
835 
836  near_boundary_unit_tests = .false.
837 
838  ! Unit tests for boundary_k_range
839  test_name = 'Surface boundary spans the entire top cell'
840  h_l = (/5.,5./)
841  call boundary_k_range(surface, nk, h_l, 5., k_top, zeta_top, k_bot, zeta_bot)
842  near_boundary_unit_tests = near_boundary_unit_tests .or. &
843  test_boundary_k_range(k_top, zeta_top, k_bot, zeta_bot, 1, 0., 1, 1., test_name, verbose)
844 
845  test_name = 'Surface boundary spans the entire column'
846  h_l = (/5.,5./)
847  call boundary_k_range(surface, nk, h_l, 10., k_top, zeta_top, k_bot, zeta_bot)
848  near_boundary_unit_tests = near_boundary_unit_tests .or. &
849  test_boundary_k_range(k_top, zeta_top, k_bot, zeta_bot, 1, 0., 2, 1., test_name, verbose)
850 
851  test_name = 'Bottom boundary spans the entire bottom cell'
852  h_l = (/5.,5./)
853  call boundary_k_range(bottom, nk, h_l, 5., k_top, zeta_top, k_bot, zeta_bot)
854  near_boundary_unit_tests = near_boundary_unit_tests .or. &
855  test_boundary_k_range(k_top, zeta_top, k_bot, zeta_bot, 2, 1., 2, 0., test_name, verbose)
856 
857  test_name = 'Bottom boundary spans the entire column'
858  h_l = (/5.,5./)
859  call boundary_k_range(bottom, nk, h_l, 10., k_top, zeta_top, k_bot, zeta_bot)
860  near_boundary_unit_tests = near_boundary_unit_tests .or. &
861  test_boundary_k_range(k_top, zeta_top, k_bot, zeta_bot, 1, 1., 2, 0., test_name, verbose)
862 
863  test_name = 'Surface boundary intersects second layer'
864  h_l = (/10.,10./)
865  call boundary_k_range(surface, nk, h_l, 17.5, k_top, zeta_top, k_bot, zeta_bot)
866  near_boundary_unit_tests = near_boundary_unit_tests .or. &
867  test_boundary_k_range(k_top, zeta_top, k_bot, zeta_bot, 1, 0., 2, 0.75, test_name, verbose)
868 
869  test_name = 'Surface boundary intersects first layer'
870  h_l = (/10.,10./)
871  call boundary_k_range(surface, nk, h_l, 2.5, k_top, zeta_top, k_bot, zeta_bot)
872  near_boundary_unit_tests = near_boundary_unit_tests .or. &
873  test_boundary_k_range(k_top, zeta_top, k_bot, zeta_bot, 1, 0., 1, 0.25, test_name, verbose)
874 
875  test_name = 'Surface boundary is deeper than column thickness'
876  h_l = (/10.,10./)
877  call boundary_k_range(surface, nk, h_l, 21.0, k_top, zeta_top, k_bot, zeta_bot)
878  near_boundary_unit_tests = near_boundary_unit_tests .or. &
879  test_boundary_k_range(k_top, zeta_top, k_bot, zeta_bot, 1, 0., 2, 1., test_name, verbose)
880 
881  test_name = 'Bottom boundary intersects first layer'
882  h_l = (/10.,10./)
883  call boundary_k_range(bottom, nk, h_l, 17.5, k_top, zeta_top, k_bot, zeta_bot)
884  near_boundary_unit_tests = near_boundary_unit_tests .or. &
885  test_boundary_k_range(k_top, zeta_top, k_bot, zeta_bot, 1, 0.75, 2, 0., test_name, verbose)
886 
887  test_name = 'Bottom boundary intersects second layer'
888  h_l = (/10.,10./)
889  call boundary_k_range(bottom, nk, h_l, 2.5, k_top, zeta_top, k_bot, zeta_bot)
890  near_boundary_unit_tests = near_boundary_unit_tests .or. &
891  test_boundary_k_range(k_top, zeta_top, k_bot, zeta_bot, 2, 0.25, 2, 0., test_name, verbose)
892 
893  ! All cases in this section have hbl which are equal to the column thicknesses
894  test_name = 'Equal hbl and same layer thicknesses (gradient from right to left)'
895  hbl_l = 10; hbl_r = 10
896  h_l = (/5.,5./) ; h_r = (/5.,5./)
897  phi_l = (/0.,0./) ; phi_r = (/1.,1./)
898  phi_pp_l(1,1) = 0.; phi_pp_l(1,2) = 0.
899  phi_pp_l(2,1) = 0.; phi_pp_l(2,2) = 0.
900  phi_pp_r(1,1) = 1.; phi_pp_r(1,2) = 0.
901  phi_pp_r(2,1) = 1.; phi_pp_r(2,2) = 0.
902  ppoly0_e_l(1,1) = 0.; ppoly0_e_l(1,2) = 0.
903  ppoly0_e_l(2,1) = 0.; ppoly0_e_l(2,2) = 0.
904  ppoly0_e_r(1,1) = 1.; ppoly0_e_r(1,2) = 1.
905  ppoly0_e_r(2,1) = 1.; ppoly0_e_r(2,2) = 1.
906  khtr_u = 1.
907  ! Without limiter
908  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, &
909  ppoly0_e_l, ppoly0_e_r, method, khtr_u, f_bulk, f_layer)
910  near_boundary_unit_tests = near_boundary_unit_tests .or. &
911  test_layer_fluxes( verbose, nk, test_name, f_layer, (/-5.0,-5.0/) )
912 
913  ! same as above, but with limiter
914  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, &
915  ppoly0_e_l, ppoly0_e_r, method, khtr_u, f_bulk, f_layer, .true.)
916  near_boundary_unit_tests = near_boundary_unit_tests .or. &
917  test_layer_fluxes( verbose, nk, test_name, f_layer, (/-1.0,-1.0/) )
918 
919  test_name = 'Equal hbl and same layer thicknesses (gradient from left to right)'
920  hbl_l = 10.; hbl_r = 10.
921  h_l = (/5.,5./) ; h_r = (/5.,5./)
922  phi_l = (/1.,1./) ; phi_r = (/0.,0./)
923  phi_pp_l(1,1) = 1.; phi_pp_l(1,2) = 0.
924  phi_pp_l(2,1) = 1.; phi_pp_l(2,2) = 0.
925  phi_pp_r(1,1) = 0.; phi_pp_r(1,2) = 0.
926  phi_pp_r(2,1) = 0.; phi_pp_r(2,2) = 0.
927  ppoly0_e_l(1,1) = 1.; ppoly0_e_l(1,2) = 1.
928  ppoly0_e_l(2,1) = 1.; ppoly0_e_l(2,2) = 1.
929  ppoly0_e_r(1,1) = 0.; ppoly0_e_r(1,2) = 0.
930  ppoly0_e_r(2,1) = 0.; ppoly0_e_r(2,2) = 0.
931  khtr_u = 1.
932  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,&
933  ppoly0_e_l, ppoly0_e_r, method, khtr_u, f_bulk, f_layer)
934  near_boundary_unit_tests = near_boundary_unit_tests .or. &
935  test_layer_fluxes( verbose, nk, test_name, f_layer, (/5.0,5.0/) )
936 
937  test_name = 'Equal hbl and same layer thicknesses (no gradient)'
938  hbl_l = 10; hbl_r = 10
939  h_l = (/5.,5./) ; h_r = (/5.,5./)
940  phi_l = (/1.,1./) ; phi_r = (/1.,1./)
941  phi_pp_l(1,1) = 1.; phi_pp_l(1,2) = 0.
942  phi_pp_l(2,1) = 1.; phi_pp_l(2,2) = 0.
943  phi_pp_r(1,1) = 1.; phi_pp_r(1,2) = 0.
944  phi_pp_r(2,1) = 1.; phi_pp_r(2,2) = 0.
945  ppoly0_e_l(1,1) = 1.; ppoly0_e_l(1,2) = 1.
946  ppoly0_e_l(2,1) = 1.; ppoly0_e_l(2,2) = 1.
947  ppoly0_e_r(1,1) = 1.; ppoly0_e_r(1,2) = 1.
948  ppoly0_e_r(2,1) = 1.; ppoly0_e_r(2,2) = 1.
949  khtr_u = 1.
950  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,&
951  ppoly0_e_l, ppoly0_e_r, method, khtr_u, f_bulk, f_layer)
952  near_boundary_unit_tests = near_boundary_unit_tests .or. &
953  test_layer_fluxes( verbose, nk, test_name, f_layer, (/0.0,0.0/) )
954 
955  test_name = 'Equal hbl and different layer thicknesses (gradient right to left)'
956  hbl_l = 16.; hbl_r = 16.
957  h_l = (/10.,6./) ; h_r = (/6.,10./)
958  phi_l = (/0.,0./) ; phi_r = (/1.,1./)
959  phi_pp_l(1,1) = 0.; phi_pp_l(1,2) = 0.
960  phi_pp_l(2,1) = 0.; phi_pp_l(2,2) = 0.
961  phi_pp_r(1,1) = 1.; phi_pp_r(1,2) = 0.
962  phi_pp_r(2,1) = 1.; phi_pp_r(2,2) = 0.
963  ppoly0_e_l(1,1) = 0.; ppoly0_e_l(1,2) = 0.
964  ppoly0_e_l(2,1) = 0.; ppoly0_e_l(2,2) = 0.
965  ppoly0_e_r(1,1) = 1.; ppoly0_e_r(1,2) = 1.
966  ppoly0_e_r(2,1) = 1.; ppoly0_e_r(2,2) = 1.
967  khtr_u = 1.
968  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,&
969  ppoly0_e_l, ppoly0_e_r, method, khtr_u, f_bulk, f_layer)
970  near_boundary_unit_tests = near_boundary_unit_tests .or. &
971  test_layer_fluxes( verbose, nk, test_name, f_layer, (/-8.0,-8.0/) )
972 
973  test_name = 'Equal hbl and same layer thicknesses (diagonal tracer values)'
974  hbl_l = 10.; hbl_r = 10.
975  h_l = (/5.,5./) ; h_r = (/5.,5./)
976  phi_l = (/1.,0./) ; phi_r = (/0.,1./)
977  phi_pp_l(1,1) = 1.; phi_pp_l(1,2) = 0.
978  phi_pp_l(2,1) = 0.; phi_pp_l(2,2) = 0.
979  phi_pp_r(1,1) = 0.; phi_pp_r(1,2) = 0.
980  phi_pp_r(2,1) = 1.; phi_pp_r(2,2) = 0.
981  ppoly0_e_l(1,1) = 1.; ppoly0_e_l(1,2) = 1.
982  ppoly0_e_l(2,1) = 0.; ppoly0_e_l(2,2) = 0.
983  ppoly0_e_r(1,1) = 0.; ppoly0_e_r(1,2) = 0.
984  ppoly0_e_r(2,1) = 1.; ppoly0_e_r(2,2) = 1.
985  khtr_u = 1.
986  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,&
987  ppoly0_e_l, ppoly0_e_r, method, khtr_u, f_bulk, f_layer)
988  near_boundary_unit_tests = near_boundary_unit_tests .or. &
989  test_layer_fluxes( verbose, nk, test_name, f_layer, (/0.0,0.0/) )
990 
991  test_name = 'Different hbl and different layer thicknesses (gradient from right to left)'
992  hbl_l = 12; hbl_r = 20
993  h_l = (/6.,6./) ; h_r = (/10.,10./)
994  phi_l = (/0.,0./) ; phi_r = (/1.,1./)
995  phi_pp_l(1,1) = 0.; phi_pp_l(1,2) = 0.
996  phi_pp_l(2,1) = 0.; phi_pp_l(2,2) = 0.
997  phi_pp_r(1,1) = 1.; phi_pp_r(1,2) = 0.
998  phi_pp_r(2,1) = 1.; phi_pp_r(2,2) = 0.
999  ppoly0_e_l(1,1) = 0.; ppoly0_e_l(1,2) = 0.
1000  ppoly0_e_l(2,1) = 0.; ppoly0_e_l(2,2) = 0.
1001  ppoly0_e_r(1,1) = 1.; ppoly0_e_r(1,2) = 1.
1002  ppoly0_e_r(2,1) = 1.; ppoly0_e_r(2,2) = 1.
1003  khtr_u = 1.
1004  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,&
1005  ppoly0_e_l, ppoly0_e_r, method, khtr_u, f_bulk, f_layer)
1006  near_boundary_unit_tests = near_boundary_unit_tests .or. &
1007  test_layer_fluxes( verbose, nk, test_name, f_layer, (/-7.5,-7.5/) )
1008 
1009  ! Cases where hbl < column thickness (polynomial coefficients specified for pseudo-linear reconstruction)
1010 
1011  test_name = 'hbl < column thickness, hbl same, constant concentration each column'
1012  hbl_l = 2; hbl_r = 2
1013  h_l = (/1.,2./) ; h_r = (/1.,2./)
1014  phi_l = (/0.,0./) ; phi_r = (/1.,1./)
1015  phi_pp_l(1,1) = 0.; phi_pp_l(1,2) = 0.
1016  phi_pp_l(2,1) = 0.; phi_pp_l(2,2) = 0.
1017  phi_pp_r(1,1) = 1.; phi_pp_r(1,2) = 0.
1018  phi_pp_r(2,1) = 1.; phi_pp_r(2,2) = 0.
1019  ppoly0_e_l(1,1) = 0.; ppoly0_e_l(1,2) = 0.
1020  ppoly0_e_l(2,1) = 0.; ppoly0_e_l(2,2) = 0.
1021  ppoly0_e_r(1,1) = 1.; ppoly0_e_r(1,2) = 1.
1022  ppoly0_e_r(2,1) = 1.; ppoly0_e_r(2,2) = 1.
1023  khtr_u = 1.
1024  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,&
1025  ppoly0_e_l, ppoly0_e_r, method, khtr_u, f_bulk, f_layer)
1026  near_boundary_unit_tests = near_boundary_unit_tests .or. &
1027  test_layer_fluxes( verbose, nk, test_name, f_layer, (/-1.,-1./) )
1028 
1029  test_name = 'hbl < column thickness, hbl same, linear profile right'
1030  hbl_l = 2; hbl_r = 2
1031  h_l = (/1.,2./) ; h_r = (/1.,2./)
1032  phi_l = (/0.,0./) ; phi_r = (/0.5,2./)
1033  phi_pp_l(1,1) = 0.; phi_pp_l(1,2) = 0.
1034  phi_pp_l(2,1) = 0.; phi_pp_l(2,2) = 0.
1035  phi_pp_r(1,1) = 0.; phi_pp_r(1,2) = 1.
1036  phi_pp_r(2,1) = 1.; phi_pp_r(2,2) = 2.
1037  khtr_u = 1.
1038  ppoly0_e_l(1,1) = 0.; ppoly0_e_l(1,2) = 0.
1039  ppoly0_e_l(2,1) = 0.; ppoly0_e_l(2,2) = 0.
1040  ppoly0_e_r(1,1) = 0.; ppoly0_e_r(1,2) = 1.
1041  ppoly0_e_r(2,1) = 1.; ppoly0_e_r(2,2) = 3.
1042  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,&
1043  ppoly0_e_l, ppoly0_e_r, method, khtr_u, f_bulk, f_layer)
1044  near_boundary_unit_tests = near_boundary_unit_tests .or. &
1045  test_layer_fluxes( verbose, nk, test_name, f_layer, (/-1.,-1./) )
1046 
1047  test_name = 'hbl < column thickness, hbl same, linear profile right, khtr=2'
1048  hbl_l = 2; hbl_r = 2
1049  h_l = (/1.,2./) ; h_r = (/1.,2./)
1050  phi_l = (/0.,0./) ; phi_r = (/0.5,2./)
1051  phi_pp_l(1,1) = 0.; phi_pp_l(1,2) = 0.
1052  phi_pp_l(2,1) = 0.; phi_pp_l(2,2) = 0.
1053  phi_pp_r(1,1) = 0.; phi_pp_r(1,2) = 1.
1054  phi_pp_r(2,1) = 1.; phi_pp_r(2,2) = 2.
1055  khtr_u = 2.
1056  ppoly0_e_l(1,1) = 0.; ppoly0_e_l(1,2) = 0.
1057  ppoly0_e_l(2,1) = 0.; ppoly0_e_l(2,2) = 0.
1058  ppoly0_e_r(1,1) = 0.; ppoly0_e_r(1,2) = 1.
1059  ppoly0_e_r(2,1) = 1.; ppoly0_e_r(2,2) = 3.
1060  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, &
1061  phi_pp_r, ppoly0_e_l, ppoly0_e_r, method, khtr_u, f_layer)
1062  near_boundary_unit_tests = near_boundary_unit_tests .or. &
1063  test_layer_fluxes( verbose, nk, test_name, f_layer, (/-1.,-3./) )
1064 
1065  ! unit tests for layer by layer method
1066  test_name = 'Different hbl and different column thicknesses (gradient from right to left)'
1067  hbl_l = 12; hbl_r = 20
1068  h_l = (/6.,6./) ; h_r = (/10.,10./)
1069  phi_l = (/0.,0./) ; phi_r = (/1.,1./)
1070  phi_pp_l(1,1) = 0.; phi_pp_l(1,2) = 0.
1071  phi_pp_l(2,1) = 0.; phi_pp_l(2,2) = 0.
1072  phi_pp_r(1,1) = 1.; phi_pp_r(1,2) = 0.
1073  phi_pp_r(2,1) = 1.; phi_pp_r(2,2) = 0.
1074  ppoly0_e_l(1,1) = 0.; ppoly0_e_l(1,2) = 0.
1075  ppoly0_e_l(2,1) = 0.; ppoly0_e_l(2,2) = 0.
1076  ppoly0_e_r(1,1) = 1.; ppoly0_e_r(1,2) = 1.
1077  ppoly0_e_r(2,1) = 1.; ppoly0_e_r(2,2) = 1.
1078  khtr_u = 1.
1079  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, &
1080  phi_pp_r, ppoly0_e_l, ppoly0_e_r, method, khtr_u, f_layer)
1081  near_boundary_unit_tests = near_boundary_unit_tests .or. &
1082  test_layer_fluxes( verbose, nk, test_name, f_layer, (/-7.5,-7.5/) )
1083 
1084  test_name = 'Different hbl and different column thicknesses (linear profile right)'
1085 
1086  hbl_l = 15; hbl_r = 6
1087  h_l = (/10.,10./) ; h_r = (/12.,10./)
1088  phi_l = (/0.,0./) ; phi_r = (/1.,3./)
1089  phi_pp_l(1,1) = 0.; phi_pp_l(1,2) = 0.
1090  phi_pp_l(2,1) = 0.; phi_pp_l(2,2) = 0.
1091  phi_pp_r(1,1) = 0.; phi_pp_r(1,2) = 2.
1092  phi_pp_r(2,1) = 2.; phi_pp_r(2,2) = 2.
1093  ppoly0_e_l(1,1) = 0.; ppoly0_e_l(1,2) = 0.
1094  ppoly0_e_l(2,1) = 0.; ppoly0_e_l(2,2) = 0.
1095  ppoly0_e_r(1,1) = 0.; ppoly0_e_r(1,2) = 2.
1096  ppoly0_e_r(2,1) = 2.; ppoly0_e_r(2,2) = 4.
1097  khtr_u = 1.
1098  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, &
1099  phi_pp_r, ppoly0_e_l, ppoly0_e_r, method, khtr_u, f_layer)
1100  near_boundary_unit_tests = near_boundary_unit_tests .or. &
1101  test_layer_fluxes( verbose, nk, test_name, f_layer, (/-3.75,0.0/) )
1102 end function near_boundary_unit_tests
1103 
1104 !> Returns true if output of near-boundary unit tests does not match correct computed values
1105 !! and conditionally writes results to stream
1106 logical function test_layer_fluxes(verbose, nk, test_name, F_calc, F_ans)
1107  logical, intent(in) :: verbose !< If true, write results to stdout
1108  character(len=80), intent(in) :: test_name !< Brief description of the unit test
1109  integer, intent(in) :: nk !< Number of layers
1110  real, dimension(nk), intent(in) :: f_calc !< Fluxes of the unitless tracer from the algorithm [s^-1]
1111  real, dimension(nk), intent(in) :: f_ans !< Fluxes of the unitless tracer calculated by hand [s^-1]
1112  ! Local variables
1113  integer :: k
1114  integer, parameter :: stdunit = stdout
1115 
1116  test_layer_fluxes = .false.
1117  do k=1,nk
1118  if ( f_calc(k) /= f_ans(k) ) then
1119  test_layer_fluxes = .true.
1120  write(stdunit,*) "MOM_lateral_boundary_diffusion, UNIT TEST FAILED: ", test_name
1121  write(stdunit,10) k, f_calc(k), f_ans(k)
1122  elseif (verbose) then
1123  write(stdunit,10) k, f_calc(k), f_ans(k)
1124  endif
1125  enddo
1126 
1127 10 format("Layer=",i3," F_calc=",f20.16," F_ans",f20.16)
1128 end function test_layer_fluxes
1129 
1130 !> Return true if output of unit tests for boundary_k_range does not match answers
1131 logical function test_boundary_k_range(k_top, zeta_top, k_bot, zeta_bot, k_top_ans, zeta_top_ans,&
1132  k_bot_ans, zeta_bot_ans, test_name, verbose)
1133  integer :: k_top !< Index of cell containing top of boundary
1134  real :: zeta_top !< Nondimension position
1135  integer :: k_bot !< Index of cell containing bottom of boundary
1136  real :: zeta_bot !< Nondimension position
1137  integer :: k_top_ans !< Index of cell containing top of boundary
1138  real :: zeta_top_ans !< Nondimension position
1139  integer :: k_bot_ans !< Index of cell containing bottom of boundary
1140  real :: zeta_bot_ans !< Nondimension position
1141  character(len=80) :: test_name !< Name of the unit test
1142  logical :: verbose !< If true always print output
1143 
1144  integer, parameter :: stdunit = stdout
1145 
1146  test_boundary_k_range = k_top .ne. k_top_ans
1147  test_boundary_k_range = test_boundary_k_range .or. (zeta_top .ne. zeta_top_ans)
1148  test_boundary_k_range = test_boundary_k_range .or. (k_bot .ne. k_bot_ans)
1149  test_boundary_k_range = test_boundary_k_range .or. (zeta_bot .ne. zeta_bot_ans)
1150 
1151  if (test_boundary_k_range) write(stdunit,*) "UNIT TEST FAILED: ", test_name
1152  if (test_boundary_k_range .or. verbose) then
1153  write(stdunit,20) "k_top", k_top, "k_top_ans", k_top_ans
1154  write(stdunit,20) "k_bot", k_bot, "k_bot_ans", k_bot_ans
1155  write(stdunit,30) "zeta_top", zeta_top, "zeta_top_ans", zeta_top_ans
1156  write(stdunit,30) "zeta_bot", zeta_bot, "zeta_bot_ans", zeta_bot_ans
1157  endif
1158 
1159  20 format(a,"=",i3,x,a,"=",i3)
1160  30 format(a,"=",f20.16,x,a,"=",f20.16)
1161 
1162 
1163 end function test_boundary_k_range
1164 
1165 !> \namespace mom_lateral_boundary_diffusion
1166 !!
1167 !! \section section_LBD The Lateral Boundary Diffusion (LBD) framework
1168 !!
1169 !! The LBD framework accounts for the effects of diabatic mesoscale fluxes
1170 !! within surface and bottom boundary layers. Unlike the equivalent adiabatic
1171 !! fluxes, which is applied along neutral density surfaces, LBD is purely
1172 !! horizontal.
1173 !!
1174 !! The bottom boundary layer fluxes remain to be implemented, although most
1175 !! of the steps needed to do so have already been added and tested.
1176 !!
1177 !! Boundary lateral diffusion can be applied using one of the three methods:
1178 !!
1179 !! * [Method #1: Along layer](@ref section_method2) (default);
1180 !! * [Method #2: Bulk layer](@ref section_method1);
1181 !!
1182 !! A brief summary of these methods is provided below.
1183 !!
1184 !! \subsection section_method1 Along layer approach (Method #1)
1185 !!
1186 !! This is the recommended and more straight forward method where diffusion is
1187 !! applied layer by layer using only information from neighboring cells.
1188 !!
1189 !! Step #1: compute vertical indices containing boundary layer (boundary_k_range).
1190 !! For the TOP boundary layer, these are:
1191 !!
1192 !! k_top, k_bot, zeta_top, zeta_bot
1193 !!
1194 !! Step #2: calculate the diffusive flux at each layer:
1195 !!
1196 !! \f[ F_{k} = -KHTR \times h_{eff}(k) \times (\phi_R(k) - \phi_L(k)), \f]
1197 !! where h_eff is the [harmonic mean](@ref section_harmonic_mean) of the layer thickness
1198 !! in the left and right columns. This method does not require a limiter since KHTR
1199 !! is already limted based on a diffusive CFL condition prior to the call of this
1200 !! module.
1201 !!
1202 !! Step #3: option to linearly decay the flux from k_bot_min to k_bot_max:
1203 !!
1204 !! If LBD_LINEAR_TRANSITION = True and k_bot_diff > 1, the diffusive flux will decay
1205 !! linearly between the top interface of the layer containing the minimum boundary
1206 !! layer depth (k_bot_min) and the lower interface of the layer containing the
1207 !! maximum layer depth (k_bot_max).
1208 !!
1209 !! \subsection section_method2 Bulk layer approach (Method #2)
1210 !!
1211 !! Apply the lateral boundary diffusive fluxes calculated from a 'bulk model'.This
1212 !! is a lower order representation (Kraus-Turner like approach) which assumes that
1213 !! eddies are acting along well mixed layers (i.e., eddies do not know care about
1214 !! vertical tracer gradients within the boundary layer).
1215 !!
1216 !! Step #1: compute vertical indices containing boundary layer (boundary_k_range).
1217 !! For the TOP boundary layer, these are:
1218 !!
1219 !! k_top, k_bot, zeta_top, zeta_bot
1220 !!
1221 !! Step #2: compute bulk averages (thickness weighted) tracer averages (phi_L and phi_R),
1222 !! then calculate the bulk diffusive flux (F_{bulk}):
1223 !!
1224 !! \f[ F_{bulk} = -KHTR \times h_{eff} \times (\phi_R - \phi_L), \f]
1225 !! where h_eff is the [harmonic mean](@ref section_harmonic_mean) of the boundary layer depth
1226 !! in the left and right columns (\f[ HBL_L \f] and \f[ HBL_R \f], respectively).
1227 !!
1228 !! Step #3: decompose F_bulk onto individual layers:
1229 !!
1230 !! \f[ F_{layer}(k) = F_{bulk} \times h_{frac}(k) , \f]
1231 !!
1232 !! where h_{frac} is
1233 !!
1234 !! \f[ h_{frac}(k) = h_u(k) \times \frac{1}{\sum(h_u)}. \f]
1235 !!
1236 !! h_u is the [harmonic mean](@ref section_harmonic_mean) of thicknesses at each layer.
1237 !! Special care (layer reconstruction) must be taken at k_min = min(k_botL, k_bot_R).
1238 !!
1239 !! Step #4: option to linearly decay the flux from k_bot_min to k_bot_max:
1240 !!
1241 !! If LBD_LINEAR_TRANSITION = True and k_bot_diff > 1, the diffusive flux will decay
1242 !! linearly between the top interface of the layer containing the minimum boundary
1243 !! layer depth (k_bot_min) and the lower interface of the layer containing the
1244 !! maximum layer depth (k_bot_max).
1245 !!
1246 !! Step #5: limit the tracer flux so that 1) only down-gradient fluxes are applied,
1247 !! and 2) the flux cannot be larger than F_max, which is defined using the tracer
1248 !! gradient:
1249 !!
1250 !! \f[ F_{max} = -0.2 \times [(V_R(k) \times \phi_R(k)) - (V_L(k) \times \phi_L(k))], \f]
1251 !! where V is the cell volume. Why 0.2?
1252 !! t=0 t=inf
1253 !! 0 .2
1254 !! 0 1 0 .2.2.2
1255 !! 0 .2
1256 !!
1257 !! \subsection section_harmonic_mean Harmonic Mean
1258 !!
1259 !! The harmonic mean (HM) betwen h1 and h2 is defined as:
1260 !!
1261 !! \f[ HM = \frac{2 \times h1 \times h2}{h1 + h2} \f]
1262 !!
1263 end module mom_lateral_boundary_diffusion
mom_lateral_boundary_diffusion::lateral_boundary_diffusion_cs
Sets parameters for lateral boundary mixing module.
Definition: MOM_lateral_boundary_diffusion.F90:39
mom_verticalgrid
Provides a transparent vertical ocean grid type and supporting routines.
Definition: MOM_verticalGrid.F90:2
mom_tracer_registry::tracer_type
The tracer type.
Definition: MOM_tracer_registry.F90:38
mom_file_parser::log_version
An overloaded interface to log version information about modules.
Definition: MOM_file_parser.F90:109
mom_diag_mediator
The subroutines here provide convenient wrappers to the fms diag_manager interfaces with additional d...
Definition: MOM_diag_mediator.F90:3
mom_energetic_pbl
Energetically consistent planetary boundary layer parameterization.
Definition: MOM_energetic_PBL.F90:2
mom_domains::pass_var
Do a halo update on an array.
Definition: MOM_domains.F90:54
mom_remapping::remapping_cs
Container for remapping parameters.
Definition: MOM_remapping.F90:24
mom_tracer_registry
This module contains the tracer_registry_type and the subroutines that handle registration of tracers...
Definition: MOM_tracer_registry.F90:5
mom_file_parser::param_file_type
A structure that can be parsed to read and document run-time parameters.
Definition: MOM_file_parser.F90:54
mom_file_parser::get_param
An overloaded interface to read and log the values of various types of parameters.
Definition: MOM_file_parser.F90:102
mom_unit_scaling::unit_scale_type
Describes various unit conversion factors.
Definition: MOM_unit_scaling.F90:14
mom_remapping
Provides column-wise vertical remapping functions.
Definition: MOM_remapping.F90:2
mom_diag_mediator::post_data
Make a diagnostic available for averaging or output.
Definition: MOM_diag_mediator.F90:70
mom_verticalgrid::verticalgrid_type
Describes the vertical ocean grid, including unit conversion factors.
Definition: MOM_verticalGrid.F90:24
mom_domains
Describes the decomposed MOM domain and has routines for communications across PEs.
Definition: MOM_domains.F90:2
mom_cpu_clock
Wraps the MPP cpu clock functions.
Definition: MOM_cpu_clock.F90:2
mom_file_parser
The MOM6 facility to parse input files for runtime parameters.
Definition: MOM_file_parser.F90:2
mom_tracer_registry::tracer_registry_type
Type to carry basic tracer information.
Definition: MOM_tracer_registry.F90:138
mom_energetic_pbl::energetic_pbl_cs
This control structure holds parameters for the MOM_energetic_PBL module.
Definition: MOM_energetic_PBL.F90:34
mom_lateral_boundary_diffusion
Calculates and applies diffusive fluxes as a parameterization of lateral mixing (non-neutral) by meso...
Definition: MOM_lateral_boundary_diffusion.F90:4
mom_diabatic_driver
This routine drives the diabatic/dianeutral physics for MOM.
Definition: MOM_diabatic_driver.F90:2
mom_grid
Provides the ocean grid type.
Definition: MOM_grid.F90:2
mom_unit_scaling
Provides a transparent unit rescaling type to facilitate dimensional consistency testing.
Definition: MOM_unit_scaling.F90:2
mom_cvmix_kpp
Provides the K-Profile Parameterization (KPP) of Large et al., 1994, via CVMix.
Definition: MOM_CVMix_KPP.F90:2
mom_error_handler
Routines for error handling and I/O management.
Definition: MOM_error_handler.F90:2
mom_cvmix_kpp::kpp_cs
Control structure for containing KPP parameters/data.
Definition: MOM_CVMix_KPP.F90:71
mom_grid::ocean_grid_type
Ocean grid type. See mom_grid for details.
Definition: MOM_grid.F90:26
mom_diabatic_driver::diabatic_cs
Control structure for this module.
Definition: MOM_diabatic_driver.F90:92
mom_diag_mediator::diag_ctrl
The following data type a list of diagnostic fields an their variants, as well as variables that cont...
Definition: MOM_diag_mediator.F90:240