MOM_hor_bnd_diffusion.F90

1! This file is part of MOM6, the Modular Ocean Model version 6.
2! See the LICENSE file for licensing information.
3! SPDX-License-Identifier: Apache-2.0
4
5!> Calculates and applies diffusive fluxes as a parameterization of horizontal mixing (non-neutral) by
6!! mesoscale eddies near the top and bottom (to be implemented) boundary layers of the ocean.
7
9
10use mom_cpu_clock, only : cpu_clock_id, cpu_clock_begin, cpu_clock_end
11use mom_cpu_clock, only : clock_module
12use mom_checksums, only : hchksum
13use mom_domains, only : pass_var
14use mom_diag_mediator, only : diag_ctrl, time_type
15use mom_diag_mediator, only : post_data, register_diag_field
16use mom_error_handler, only : mom_error, mom_mesg, fatal, is_root_pe
17use mom_file_parser, only : get_param, log_version, param_file_type
18use mom_grid, only : ocean_grid_type
21use mom_remapping, only : remappingschemesdoc, remappingdefaultscheme
23use mom_tracer_registry, only : tracer_registry_type, tracer_type
30use mom_io, only : stdout, stderr
31
32implicit none ; private
33
36
37! Private parameters to avoid doing string comparisons for bottom or top boundary layer
38integer, public, parameter :: surface = -1 !< Set a value that corresponds to the surface boundary
39integer, public, parameter :: bottom = 1 !< Set a value that corresponds to the bottom boundary
40#include <MOM_memory.h>
41
42!> Sets parameters for horizontal boundary mixing module.
43type, public :: hbd_cs ; private
44 logical :: debug !< If true, write verbose checksums for debugging.
45 integer :: deg !< Degree of polynomial reconstruction.
46 integer :: hbd_nk !< Maximum number of levels in the HBD grid [nondim]
47 integer :: surface_boundary_scheme !< Which boundary layer scheme to use
48 !! 1. ePBL; 2. KPP
49 logical :: limiter !< Controls whether a flux limiter is applied in the
50 !! native grid (default is true).
51 logical :: limiter_remap !< Controls whether a flux limiter is applied in the
52 !! remapped grid (default is false).
53 logical :: linear !< If True, apply a linear transition at the base/top of the boundary.
54 !! The flux will be fully applied at k=k_min and zero at k=k_max.
55 real :: h_subroundoff !< A thickness that is so small that it can be added to a thickness of
56 !! Angstrom or larger without changing it at the bit level [H ~> m or kg m-2].
57 !! If Angstrom is 0 or exceedingly small, this is negligible compared to 1e-17 m.
58 ! HBD dynamic grids
59 real, allocatable, dimension(:,:,:) :: hbd_grd_u !< HBD thicknesses at t-points adjacent to
60 !! u-points [H ~> m or kg m-2]
61 real, allocatable, dimension(:,:,:) :: hbd_grd_v !< HBD thicknesses at t-points adjacent to
62 !! v-points (left and right) [H ~> m or kg m-2]
63 integer, allocatable, dimension(:,:) :: hbd_u_kmax !< Maximum vertical index in hbd_grd_u [nondim]
64 integer, allocatable, dimension(:,:) :: hbd_v_kmax !< Maximum vertical index in hbd_grd_v [nondim]
65 type(remapping_cs) :: remap_cs !< Control structure to hold remapping configuration.
66 type(kpp_cs), pointer :: kpp_csp => null() !< KPP control structure needed to get BLD.
67 type(energetic_pbl_cs), pointer :: energetic_pbl_csp => null() !< ePBL control structure needed to get BLD.
68 type(diag_ctrl), pointer :: diag => null() !< A structure that is used to
69 !! regulate the timing of diagnostic output.
70end type hbd_cs
71
72! This include declares and sets the variable "version".
73#include "version_variable.h"
74character(len=40) :: mdl = "MOM_hor_bnd_diffusion" !< Name of this module
75integer :: id_clock_hbd !< CPU clock for hbd
76
77contains
78
79!> Initialization routine that reads runtime parameters and sets up pointers to other control structures that might be
80!! needed for horizontal boundary diffusion.
81logical function hor_bnd_diffusion_init(Time, G, GV, US, param_file, diag, diabatic_CSp, CS)
82 type(time_type), target, intent(in) :: time !< Time structure
83 type(ocean_grid_type), intent(in) :: g !< Grid structure
84 type(verticalgrid_type), intent(in) :: gv !< ocean vertical grid structure
85 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
86 type(param_file_type), intent(in) :: param_file !< Parameter file structure
87 type(diag_ctrl), target, intent(inout) :: diag !< Diagnostics control structure
88 type(diabatic_cs), pointer :: diabatic_csp !< KPP control structure needed to get BLD
89 type(hbd_cs), pointer :: cs !< Horizontal boundary mixing control structure
90
91 ! local variables
92 character(len=80) :: string ! Temporary strings
93 logical :: boundary_extrap ! controls if boundary extrapolation is used in the HBD code
94 logical :: om4_remap_via_sub_cells ! Use the OM4-era ramap_via_sub_cells for HBD
95 logical :: debug !< If true, write verbose checksums for debugging purposes
96
97 if (ASSOCIATED(cs)) then
98 call mom_error(fatal, "hor_bnd_diffusion_init called with associated control structure.")
99 return
100 endif
101
102 ! Log this module and master switch for turning it on/off
103 call get_param(param_file, mdl, "USE_HORIZONTAL_BOUNDARY_DIFFUSION", hor_bnd_diffusion_init, &
104 default=.false., do_not_log=.true.)
105 call log_version(param_file, mdl, version, &
106 "This module implements horizontal diffusion of tracers near boundaries", &
107 all_default=.not.hor_bnd_diffusion_init)
108 call get_param(param_file, mdl, "USE_HORIZONTAL_BOUNDARY_DIFFUSION", hor_bnd_diffusion_init, &
109 "If true, enables the horizontal boundary tracer's diffusion module.", &
110 default=.false.)
111 if (.not. hor_bnd_diffusion_init) return
112
113 allocate(cs)
114 cs%diag => diag
115 cs%H_subroundoff = gv%H_subroundoff
116 call extract_diabatic_member(diabatic_csp, kpp_csp=cs%KPP_CSp)
117 call extract_diabatic_member(diabatic_csp, energetic_pbl_csp=cs%energetic_PBL_CSp)
118
119 ! max. number of vertical layers
120 cs%hbd_nk = 2 + (gv%ke*2)
121 ! allocate the hbd grids and k_max
122 allocate(cs%hbd_grd_u(szib_(g),szj_(g),cs%hbd_nk), source=0.0)
123 allocate(cs%hbd_grd_v(szi_(g),szjb_(g),cs%hbd_nk), source=0.0)
124 allocate(cs%hbd_u_kmax(szib_(g),szj_(g)), source=0)
125 allocate(cs%hbd_v_kmax(szi_(g),szjb_(g)), source=0)
126
127 cs%surface_boundary_scheme = -1
128 if ( .not. ASSOCIATED(cs%energetic_PBL_CSp) .and. .not. ASSOCIATED(cs%KPP_CSp) ) then
129 call mom_error(fatal,"Horizontal boundary diffusion is true, but no valid boundary layer scheme was found")
130 endif
131
132 ! Read all relevant parameters and write them to the model log.
133 call get_param(param_file, mdl, "HBD_LINEAR_TRANSITION", cs%linear, &
134 "If True, apply a linear transition at the base/top of the boundary. \n"//&
135 "The flux will be fully applied at k=k_min and zero at k=k_max.", default=.false.)
136 call get_param(param_file, mdl, "APPLY_LIMITER", cs%limiter, &
137 "If True, apply a flux limiter in the native grid.", default=.true.)
138 call get_param(param_file, mdl, "APPLY_LIMITER_REMAP", cs%limiter_remap, &
139 "If True, apply a flux limiter in the remapped grid.", default=.false.)
140 call get_param(param_file, mdl, "HBD_BOUNDARY_EXTRAP", boundary_extrap, &
141 "Use boundary extrapolation in HBD code", &
142 default=.false.)
143 call get_param(param_file, mdl, "HBD_REMAPPING_SCHEME", string, &
144 "This sets the reconstruction scheme used "//&
145 "for vertical remapping for all variables. "//&
146 "It can be one of the following schemes: "//&
147 trim(remappingschemesdoc), default=remappingdefaultscheme)
148 call get_param(param_file, mdl, "REMAPPING_USE_OM4_SUBCELLS", om4_remap_via_sub_cells, &
149 do_not_log=.true., default=.true.)
150
151 call get_param(param_file, mdl, "HBD_REMAPPING_USE_OM4_SUBCELLS", om4_remap_via_sub_cells, &
152 "If true, use the OM4 remapping-via-subcells algorithm for horizontal boundary diffusion. "//&
153 "See REMAPPING_USE_OM4_SUBCELLS for details. "//&
154 "We recommend setting this option to false.", default=om4_remap_via_sub_cells)
155
156 ! GMM, TODO: add HBD params to control optional arguments in initialize_remapping.
157 call initialize_remapping( cs%remap_CS, string, boundary_extrapolation=boundary_extrap, &
158 om4_remap_via_sub_cells=om4_remap_via_sub_cells, &
159 check_reconstruction=.false., check_remapping=.false., &
160 h_neglect=cs%H_subroundoff, h_neglect_edge=cs%H_subroundoff)
161 call extract_member_remapping_cs(cs%remap_CS, degree=cs%deg)
162 call get_param(param_file, mdl, "DEBUG", debug, &
163 default=.false., debuggingparam=.true., do_not_log=.true.)
164 call get_param(param_file, mdl, "HBD_DEBUG", cs%debug, &
165 "If true, write out verbose debugging data in the HBD module.", &
166 default=debug, debuggingparam=.true.)
167
168 id_clock_hbd = cpu_clock_id('(Ocean HBD)', grain=clock_module)
169
170end function hor_bnd_diffusion_init
171
172!> Driver routine for calculating horizontal diffusive fluxes near the top and bottom boundaries.
173!! Diffusion is applied using only information from neighboring cells, as follows:
174!! 1) remap tracer to a z* grid (HBD grid)
175!! 2) calculate diffusive tracer fluxes (F) in the HBD grid using a layer by layer approach
176!! 3) remap fluxes to the native grid
177!! 4) update tracer by adding the divergence of F
178subroutine hor_bnd_diffusion(G, GV, US, h, Coef_x, Coef_y, dt, Reg, visc, CS)
179 type(ocean_grid_type), intent(inout) :: g !< Grid type
180 type(verticalgrid_type), intent(in) :: gv !< ocean vertical grid structure
181 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
182 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thickness [H ~> m or kg m-2]
183 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)+1), intent(in) :: coef_x !< dt * Kh * dy / dx at u-points [L2 ~> m2]
184 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)+1), intent(in) :: coef_y !< dt * Kh * dx / dy at v-points [L2 ~> m2]
185 real, intent(in) :: dt !< Tracer time step * I_numitts
186 !! (I_numitts in tracer_hordiff) [T ~> s]
187 type(tracer_registry_type), pointer :: reg !< Tracer registry
188 type(vertvisc_type), intent(in) :: visc !< Structure with vertical viscosities,
189 !! boundary layer properties and related fields
190 type(hbd_cs), pointer :: cs !< Control structure for this module
191
192 ! Local variables
193 real, dimension(SZI_(G),SZJ_(G)) :: hbl !< Boundary layer depth [H ~> m or kg m-2]
194 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)) :: uflx !< Zonal flux of tracer [conc H L2 ~> conc m3 or conc kg]
195 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)) :: vflx !< Meridional flux of tracer
196 !! [conc H L2 ~> conc m3 or conc kg]
197 real, dimension(SZIB_(G),SZJ_(G)) :: uwork_2d !< Layer summed u-flux transport
198 !! [conc H L2 ~> conc m3 or conc kg]
199 real, dimension(SZI_(G),SZJB_(G)) :: vwork_2d !< Layer summed v-flux transport
200 !! [conc H L2 ~> conc m3 or conc kg]
201 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: tendency !< tendency array for diagnostics at first in
202 !! [H conc T-1 ~> m conc s-1 or kg m-2 conc s-1],
203 !! then converted to [conc T-1 ~> conc s-1].
204 ! For temperature these units are
205 ! [C H T-1 ~> degC m s-1 or degC kg m-2 s-1] and
206 ! then [C T-1 ~> degC s-1].
207 real, dimension(SZI_(G),SZJ_(G)) :: tendency_2d !< depth integrated content tendency for diagnostics in
208 !! [H conc T-1 ~> m conc s-1 or kg m-2 conc s-1].
209 !! For temperature these units are
210 !! [C H T-1 ~> degC m s-1 or degC kg m-2 s-1].
211 type(tracer_type), pointer :: tracer => null() !< Pointer to the current tracer [conc]
212 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: tracer_old !< local copy of the initial tracer concentration,
213 !! only used to compute tendencies [conc].
214 real :: tracer_int_prev !< Globally integrated tracer before HBD is applied, in mks units [conc kg]
215 real :: tracer_int_end !< Integrated tracer after HBD is applied, in mks units [conc kg]
216 real :: idt !< inverse of the time step [T-1 ~> s-1]
217 character(len=256) :: mesg !< Message for error messages.
218 integer :: i, j, k, m !< indices to loop over
219
220 call cpu_clock_begin(id_clock_hbd)
221 idt = 1./dt
222
223 if (associated(visc%h_ML)) then
224 hbl(:,:) = visc%h_ML(:,:)
225 else
226 call mom_error(fatal, "hor_bnd_diffusion requires that visc%h_ML is associated.")
227 endif
228 ! This halo update is probably not necessary because visc%h_ML has valid halo data.
229 call pass_var(hbl, g%Domain, halo=1)
230
231 ! build HBD grid
232 call hbd_grid(surface, g, gv, hbl, h, cs)
233
234 do m = 1,reg%ntr
235 ! current tracer
236 tracer => reg%tr(m)
237
238 if (cs%debug) then
239 call hchksum(tracer%t, "before HBD "//tracer%name, g%HI, unscale=tracer%conc_scale)
240 endif
241
242 ! for diagnostics
243 if (tracer%id_hbdxy_conc > 0 .or. tracer%id_hbdxy_cont > 0 .or. tracer%id_hbdxy_cont_2d > 0 .or. cs%debug) then
244 tendency(:,:,:) = 0.0
245 tracer_old(:,:,:) = tracer%t(:,:,:)
246 endif
247
248 ! Diffusive fluxes in the i- and j-direction
249 uflx(:,:,:) = 0.
250 vflx(:,:,:) = 0.
251
252 ! HBD layer by layer
253 do j=g%jsc,g%jec
254 do i=g%isc-1,g%iec
255 if (g%mask2dCu(i,j)>0.) then
256 call fluxes_layer_method(surface, gv%ke, hbl(i,j), hbl(i+1,j), &
257 h(i,j,:), h(i+1,j,:), tracer%t(i,j,:), tracer%t(i+1,j,:), &
258 coef_x(i,j,:), uflx(i,j,:), g%areaT(i,j), g%areaT(i+1,j), cs%hbd_u_kmax(i,j), &
259 cs%hbd_grd_u(i,j,:), cs)
260 endif
261 enddo
262 enddo
263 do j=g%jsc-1,g%jec
264 do i=g%isc,g%iec
265 if (g%mask2dCv(i,j)>0.) then
266 call fluxes_layer_method(surface, gv%ke, hbl(i,j), hbl(i,j+1), &
267 h(i,j,:), h(i,j+1,:), tracer%t(i,j,:), tracer%t(i,j+1,:), &
268 coef_y(i,j,:), vflx(i,j,:), g%areaT(i,j), g%areaT(i,j+1), cs%hbd_v_kmax(i,j), &
269 cs%hbd_grd_v(i,j,:), cs)
270 endif
271 enddo
272 enddo
273
274 ! Update the tracer fluxes
275 do k=1,gv%ke ; do j=g%jsc,g%jec ; do i=g%isc,g%iec
276 if (g%mask2dT(i,j)>0.) then
277 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) ) ))* &
278 g%IareaT(i,j) / ( h(i,j,k) + gv%H_subroundoff )
279
280 if (tracer%id_hbdxy_conc > 0 .or. tracer%id_hbdxy_cont > 0 .or. tracer%id_hbdxy_cont_2d > 0 ) then
281 tendency(i,j,k) = ((uflx(i-1,j,k)-uflx(i,j,k)) + (vflx(i,j-1,k)-vflx(i,j,k))) * &
282 g%IareaT(i,j) * idt
283 endif
284 endif
285 enddo ; enddo ; enddo
286
287 ! Do user controlled underflow of the tracer concentrations.
288 if (tracer%conc_underflow > 0.0) then
289 do k=1,gv%ke ; do j=g%jsc,g%jec ; do i=g%isc,g%iec
290 if (abs(tracer%t(i,j,k)) < tracer%conc_underflow) tracer%t(i,j,k) = 0.0
291 enddo ; enddo ; enddo
292 endif
293
294 if (cs%debug) then
295 call hchksum(tracer%t, "after HBD "//tracer%name, g%HI, unscale=tracer%conc_scale)
296 ! tracer (native grid) integrated tracer amounts before and after HBD
297 tracer_int_prev = global_mass_integral(h, g, gv, tracer_old, scale=tracer%conc_scale)
298 tracer_int_end = global_mass_integral(h, g, gv, tracer%t, scale=tracer%conc_scale)
299 write(mesg,*) 'Total '//tracer%name//' before/after HBD:', tracer_int_prev, tracer_int_end
300 call mom_mesg(mesg)
301 endif
302
303 ! Post the tracer diagnostics
304 if (tracer%id_hbd_dfx>0) call post_data(tracer%id_hbd_dfx, uflx(:,:,:)*idt, cs%diag)
305 if (tracer%id_hbd_dfy>0) call post_data(tracer%id_hbd_dfy, vflx(:,:,:)*idt, cs%diag)
306 if (tracer%id_hbd_dfx_2d>0) then
307 uwork_2d(:,:) = 0.
308 do k=1,gv%ke ; do j=g%jsc,g%jec ; do i=g%isc-1,g%iec
309 uwork_2d(i,j) = uwork_2d(i,j) + (uflx(i,j,k) * idt)
310 enddo ; enddo ; enddo
311 call post_data(tracer%id_hbd_dfx_2d, uwork_2d, cs%diag)
312 endif
313
314 if (tracer%id_hbd_dfy_2d>0) then
315 vwork_2d(:,:) = 0.
316 do k=1,gv%ke ; do j=g%jsc-1,g%jec ; do i=g%isc,g%iec
317 vwork_2d(i,j) = vwork_2d(i,j) + (vflx(i,j,k) * idt)
318 enddo ; enddo ; enddo
319 call post_data(tracer%id_hbd_dfy_2d, vwork_2d, cs%diag)
320 endif
321
322 ! post tendency of tracer content
323 if (tracer%id_hbdxy_cont > 0) then
324 call post_data(tracer%id_hbdxy_cont, tendency, cs%diag)
325 endif
326
327 ! post depth summed tendency for tracer content
328 if (tracer%id_hbdxy_cont_2d > 0) then
329 tendency_2d(:,:) = 0.
330 do j=g%jsc,g%jec ; do i=g%isc,g%iec
331 do k=1,gv%ke
332 tendency_2d(i,j) = tendency_2d(i,j) + tendency(i,j,k)
333 enddo
334 enddo ; enddo
335 call post_data(tracer%id_hbdxy_cont_2d, tendency_2d, cs%diag)
336 endif
337
338 ! post tendency of tracer concentration; this step must be
339 ! done after posting tracer content tendency, since we alter
340 ! the tendency array and its units.
341 if (tracer%id_hbdxy_conc > 0) then
342 do k=1,gv%ke ; do j=g%jsc,g%jec ; do i=g%isc,g%iec
343 tendency(i,j,k) = tendency(i,j,k) / ( h(i,j,k) + cs%H_subroundoff )
344 enddo ; enddo ; enddo
345 call post_data(tracer%id_hbdxy_conc, tendency, cs%diag)
346 endif
347
348 enddo
349
350 call cpu_clock_end(id_clock_hbd)
351
352end subroutine hor_bnd_diffusion
353
354!> Build the HBD grid where tracers will be remapped to.
355subroutine hbd_grid(boundary, G, GV, hbl, h, CS)
356 integer, intent(in ) :: boundary !< Which boundary layer SURFACE or BOTTOM [nondim]
357 type(ocean_grid_type), intent(inout) :: G !< Grid type
358 type(verticalgrid_type), intent(in) :: GV !< ocean vertical grid structure
359 real, dimension(SZI_(G),SZJ_(G)), &
360 intent(in) :: hbl !< Boundary layer depth [H ~> m or kg m-2]
361 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
362 intent(in) :: h !< Layer thickness in the native grid [H ~> m or kg m-2]
363 type(hbd_cs), pointer :: CS !< Horizontal diffusion control structure
364
365 ! Local variables
366 real, allocatable :: dz_top(:) !< temporary HBD grid given by merge_interfaces [H ~> m or kg m-2]
367 integer :: nk, i, j, k !< number of layers in the HBD grid, and integers used in do-loops
368
369 ! reset arrays
370 cs%hbd_grd_u(:,:,:) = 0.0
371 cs%hbd_grd_v(:,:,:) = 0.0
372 cs%hbd_u_kmax(:,:) = 0
373 cs%hbd_v_kmax(:,:) = 0
374
375 do j=g%jsc,g%jec
376 do i=g%isc-1,g%iec
377 if (g%mask2dCu(i,j)>0.) then
378 call merge_interfaces(gv%ke, h(i,j,:), h(i+1,j,:), hbl(i,j), hbl(i+1,j), &
379 cs%H_subroundoff, dz_top)
380 nk = SIZE(dz_top)
381 if (nk > cs%hbd_nk) then
382 write(*,*)'nk, CS%hbd_nk', nk, cs%hbd_nk
383 call mom_error(fatal,"Houston, we've had a problem in hbd_grid, u-points (nk cannot be > CS%hbd_nk)")
384 endif
385
386 cs%hbd_u_kmax(i,j) = nk
387
388 ! set the HBD grid to dz_top
389 do k=1,nk
390 cs%hbd_grd_u(i,j,k) = dz_top(k)
391 enddo
392 deallocate(dz_top)
393 endif
394 enddo
395 enddo
396
397 do j=g%jsc-1,g%jec
398 do i=g%isc,g%iec
399 if (g%mask2dCv(i,j)>0.) then
400 call merge_interfaces(gv%ke, h(i,j,:), h(i,j+1,:), hbl(i,j), hbl(i,j+1), &
401 cs%H_subroundoff, dz_top)
402
403 nk = SIZE(dz_top)
404 if (nk > cs%hbd_nk) then
405 write(*,*)'nk, CS%hbd_nk', nk, cs%hbd_nk
406 call mom_error(fatal,"Houston, we've had a problem in hbd_grid, v-points (nk cannot be > CS%hbd_nk)")
407 endif
408
409 cs%hbd_v_kmax(i,j) = nk
410
411 ! set the HBD grid to dz_top
412 do k=1,nk
413 cs%hbd_grd_v(i,j,k) = dz_top(k)
414 enddo
415 deallocate(dz_top)
416 endif
417 enddo
418 enddo
419
420end subroutine hbd_grid
421
422!> Calculate the harmonic mean of two quantities [arbitrary]
423!! See \ref section_harmonic_mean.
424real function harmonic_mean(h1,h2)
425 real :: h1 !< Scalar quantity [arbitrary]
426 real :: h2 !< Scalar quantity [arbitrary]
427 if (h1 + h2 == 0.) then
428 harmonic_mean = 0.
429 else
430 harmonic_mean = 2.*(h1*h2)/(h1+h2)
431 endif
432end function harmonic_mean
433
434!> Returns the location of the minimum value in a 1D array
435!! between indices s and e.
436integer function find_minimum(x, s, e)
437 integer, intent(in) :: s !< start index
438 integer, intent(in) :: e !< end index
439 real, dimension(e), intent(in) :: x !< 1D array to be checked [arbitrary]
440
441 ! local variables
442 real :: minimum ! Minimum value in the same units as x [arbitrary]
443 integer :: location
444 integer :: i
445
446 minimum = x(s) ! assume the first is the min
447 location = s ! record its position
448 do i = s+1, e ! start with next elements
449 if (x(i) < minimum) then ! if x(i) less than the min?
450 minimum = x(i) ! Yes, a new minimum found
451 location = i ! record its position
452 endif
453 enddo
454 find_minimum = location ! return the position
455end function find_minimum
456
457!> Swaps the values of its two formal arguments.
458subroutine swap(a, b)
459 real, intent(inout) :: a !< First value to be swapped [arbitrary]
460 real, intent(inout) :: b !< Second value to be swapped [arbitrary]
461
462 ! local variables
463 real :: tmp ! A temporary copy of a [arbitrary]
464
465 tmp = a
466 a = b
467 b = tmp
468end subroutine swap
469
470!> Receives a 1D array x and sorts it into ascending order.
471subroutine sort(x, n)
472 integer, intent(in ) :: n !< Number of points in the array
473 real, dimension(n), intent(inout) :: x !< 1D array to be sorted [arbitrary]
474
475 ! local variables
476 integer :: i, location
477
478 do i = 1, n-1
479 location = find_minimum(x, i, n) ! find min from this to last
480 call swap(x(i), x(location)) ! swap this and the minimum
481 enddo
482end subroutine sort
483
484!> Returns the unique values in a 1D array.
485subroutine unique(val, n, val_unique, val_max)
486 integer, intent(in ) :: n !< Number of points in the array.
487 real, dimension(n), intent(in ) :: val !< 1D array to be checked [arbitrary]
488 real, dimension(:), allocatable, intent(inout) :: val_unique !< Returned 1D array with unique values [arbitrary]
489 real, optional, intent(in ) :: val_max !< sets the maximum value in val_unique to
490 !! this value [arbitrary]
491 ! local variables
492 real, dimension(n) :: tmp ! The list of unique values [arbitrary]
493 integer :: i, j, ii
494 real :: min_val, max_val ! The minimum and maximum values in the list [arbitrary]
495 logical :: limit
496
497 limit = .false.
498 if (present(val_max)) then
499 limit = .true.
500 if (val_max > maxval(val)) then
501 if (is_root_pe()) write(*,*)'val_max, MAXVAL(val)',val_max, maxval(val)
502 call mom_error(fatal,"Houston, we've had a problem in unique (val_max cannot be > MAXVAL(val))")
503 endif
504 endif
505
506 tmp(:) = 0.
507 min_val = minval(val)-1
508 max_val = maxval(val)
509 i = 0
510 do while (min_val<max_val)
511 i = i+1
512 min_val = minval(val, mask=val>min_val)
513 tmp(i) = min_val
514 enddo
515 ii = i
516 if (limit) then
517 do j=1,ii
518 if (tmp(j) <= val_max) i = j
519 enddo
520 endif
521 allocate(val_unique(i), source=tmp(1:i))
522end subroutine unique
523
524
525!> Given layer thicknesses (and corresponding interfaces) and BLDs in two adjacent columns,
526!! return a set of 1-d layer thicknesses whose interfaces cover all interfaces in the left
527!! and right columns plus the two BLDs. This can be used to accurately remap tracer tendencies
528!! in both columns.
529subroutine merge_interfaces(nk, h_L, h_R, hbl_L, hbl_R, H_subroundoff, h)
530 integer, intent(in ) :: nk !< Number of layers [nondim]
531 real, dimension(nk), intent(in ) :: h_L !< Layer thicknesses in the left column [H ~> m or kg m-2]
532 real, dimension(nk), intent(in ) :: h_R !< Layer thicknesses in the right column [H ~> m or kg m-2]
533 real, intent(in ) :: hbl_L !< Thickness of the boundary layer in the left column
534 !! [H ~> m or kg m-2]
535 real, intent(in ) :: hbl_R !< Thickness of the boundary layer in the right column
536 !! [H ~> m or kg m-2]
537 real, intent(in ) :: H_subroundoff !< GV%H_subroundoff [H ~> m or kg m-2]
538 real, dimension(:), allocatable, intent(inout) :: h !< Combined thicknesses [H ~> m or kg m-2]
539
540 ! Local variables
541 integer :: n !< Number of layers in eta_all
542 real, dimension(nk+1) :: eta_L, eta_R!< Interfaces in the left and right columns [H ~> m or kg m-2]
543 real, dimension(:), allocatable :: eta_all !< Combined list of interfaces in the left and right columns
544 !! plus hbl_L and hbl_R [H ~> m or kg m-2]
545 real, dimension(:), allocatable :: eta_unique !< Combined list of unique interfaces (eta_L, eta_R), possibly
546 !! hbl_L and hbl_R [H ~> m or kg m-2]
547 real :: min_depth !< Minimum depth [H ~> m or kg m-2]
548 real :: max_depth !< Maximum depth [H ~> m or kg m-2]
549 real :: max_bld !< Deepest BLD [H ~> m or kg m-2]
550 integer :: k, kk, nk1 !< loop indices (k and kk) and array size (nk1)
551
552 n = (2*nk)+3
553 allocate(eta_all(n))
554 ! compute and merge interfaces
555 eta_l(:) = 0.0 ; eta_r(:) = 0.0 ; eta_all(:) = 0.0
556 kk = 0
557 do k=2,nk+1
558 eta_l(k) = eta_l(k-1) + h_l(k-1)
559 eta_r(k) = eta_r(k-1) + h_r(k-1)
560 kk = kk + 2
561 eta_all(kk) = eta_l(k)
562 eta_all(kk+1) = eta_r(k)
563 enddo
564
565 ! add hbl_L and hbl_R into eta_all
566 eta_all(kk+2) = hbl_l
567 eta_all(kk+3) = hbl_r
568
569 ! find maximum depth
570 min_depth = min(maxval(eta_l), maxval(eta_r))
571 max_bld = max(hbl_l, hbl_r)
572 max_depth = min(min_depth, max_bld)
573
574 ! sort eta_all
575 call sort(eta_all, n)
576 ! remove duplicates from eta_all and sets maximum depth
577 call unique(eta_all, n, eta_unique, max_depth)
578
579 nk1 = SIZE(eta_unique)
580 allocate(h(nk1-1))
581 do k=1,nk1-1
582 h(k) = (eta_unique(k+1) - eta_unique(k)) + h_subroundoff
583 enddo
584end subroutine merge_interfaces
585
586!> Calculates the maximum flux that can leave a cell and uses that to apply a
587!! limiter to F_layer.
588subroutine flux_limiter(F_layer, area_L, area_R, phi_L, phi_R, h_L, h_R)
589 real, intent(inout) :: F_layer !< Tracer flux to be checked [H L2 conc ~> m3 conc]
590 real, intent(in) :: area_L !< Area of left cell [L2 ~> m2]
591 real, intent(in) :: area_R !< Area of right cell [L2 ~> m2]
592 real, intent(in) :: h_L !< Thickness of left cell [H ~> m or kg m-2]
593 real, intent(in) :: h_R !< Thickness of right cell [H ~> m or kg m-2]
594 real, intent(in) :: phi_L !< Tracer concentration in the left cell [conc]
595 real, intent(in) :: phi_R !< Tracer concentration in the right cell [conc]
596
597 ! local variables
598 real :: F_max !< maximum flux allowed [conc H L2 ~> conc m3 or conc kg]
599 ! limit the flux to 0.2 of the tracer *gradient*
600 ! Why 0.2?
601 ! t=0 t=inf
602 ! 0 .2
603 ! 0 1 0 .2.2.2
604 ! 0 .2
605 !
606 f_max = -0.2 * ((area_r*(phi_r*h_r))-(area_l*(phi_l*h_l)))
607
608 if ( sign(1.,f_layer) == sign(1., f_max)) then
609 ! Apply flux limiter calculated above
610 if (f_max >= 0.) then
611 f_layer = min(f_layer,f_max)
612 else
613 f_layer = max(f_layer,f_max)
614 endif
615 else
616 f_layer = 0.0
617 endif
618end subroutine flux_limiter
619
620!> Find the k-index range corresponding to the layers that are within the boundary-layer region
621subroutine boundary_k_range(boundary, nk, h, hbl, k_top, zeta_top, k_bot, zeta_bot)
622 integer, intent(in ) :: boundary !< SURFACE or BOTTOM [nondim]
623 integer, intent(in ) :: nk !< Number of layers [nondim]
624 real, dimension(nk), intent(in ) :: h !< Layer thicknesses of the column [H ~> m or kg m-2]
625 real, intent(in ) :: hbl !< Thickness of the boundary layer [H ~> m or kg m-2]
626 !! If surface, with respect to zbl_ref = 0.
627 !! If bottom, with respect to zbl_ref = SUM(h)
628 integer, intent( out) :: k_top !< Index of the first layer within the boundary
629 real, intent( out) :: zeta_top !< Distance from the top of a layer to the intersection of the
630 !! top extent of the boundary layer (0 at top, 1 at bottom) [nondim]
631 integer, intent( out) :: k_bot !< Index of the last layer within the boundary
632 real, intent( out) :: zeta_bot !< Distance of the lower layer to the boundary layer depth
633 !! (0 at top, 1 at bottom) [nondim]
634 ! Local variables
635 real :: htot ! Summed thickness [H ~> m or kg m-2]
636 integer :: k
637
638 ! Surface boundary layer
639 if ( boundary == surface ) then
640 k_top = 1
641 zeta_top = 0.
642 htot = 0.
643 k_bot = 1
644 zeta_bot = 0.
645 if (hbl == 0.) return
646 if (hbl >= sum(h(:))) then
647 k_bot = nk
648 zeta_bot = 1.
649 return
650 endif
651 do k=1,nk
652 htot = htot + h(k)
653 if ( htot >= hbl) then
654 k_bot = k
655 zeta_bot = 1 - (htot - hbl)/h(k)
656 return
657 endif
658 enddo
659
660 ! Bottom boundary layer
661 elseif ( boundary == bottom ) then
662 k_top = nk
663 zeta_top = 1.
664 k_bot = nk
665 zeta_bot = 0.
666 htot = 0.
667 if (hbl == 0.) return
668 if (hbl >= sum(h(:))) then
669 k_top = 1
670 zeta_top = 1.
671 return
672 endif
673 do k=nk,1,-1
674 htot = htot + h(k)
675 if (htot >= hbl) then
676 k_top = k
677 zeta_top = 1 - (htot - hbl)/h(k)
678 return
679 endif
680 enddo
681 else
682 call mom_error(fatal,"Houston, we've had a problem in boundary_k_range")
683 endif
684
685end subroutine boundary_k_range
686
687!> Calculate the horizontal boundary diffusive fluxes using the layer by layer method.
688!! See \ref section_method
689subroutine fluxes_layer_method(boundary, ke, hbl_L, hbl_R, h_L, h_R, phi_L, phi_R, &
690 khtr_u, F_layer, area_L, area_R, nk, dz_top, CS)
691
692 integer, intent(in ) :: boundary !< Which boundary layer SURFACE or BOTTOM [nondim]
693 integer, intent(in ) :: ke !< Number of layers in the native grid [nondim]
694 real, intent(in ) :: hbl_L !< Thickness of the boundary boundary
695 !! layer (left) [H ~> m or kg m-2]
696 real, intent(in ) :: hbl_R !< Thickness of the boundary boundary
697 !! layer (right) [H ~> m or kg m-2]
698 real, dimension(ke), intent(in ) :: h_L !< Thicknesses in the native grid (left) [H ~> m or kg m-2]
699 real, dimension(ke), intent(in ) :: h_R !< Thicknesses in the native grid (right) [H ~> m or kg m-2]
700 real, dimension(ke), intent(in ) :: phi_L !< Tracer values in the native grid (left) [conc]
701 real, dimension(ke), intent(in ) :: phi_R !< Tracer values in the native grid (right) [conc]
702 real, dimension(ke+1),intent(in ) :: khtr_u !< Horizontal diffusivities times the time step
703 !! at a velocity point and vertical interfaces [L2 ~> m2]
704 real, dimension(ke), intent( out) :: F_layer !< Layerwise diffusive flux at U- or V-point
705 !! in the native grid [H L2 conc ~> m3 conc]
706 real, intent(in ) :: area_L !< Area of the horizontal grid (left) [L2 ~> m2]
707 real, intent(in ) :: area_R !< Area of the horizontal grid (right) [L2 ~> m2]
708 integer, intent(in ) :: nk !< Number of layers in the HBD grid [nondim]
709 real, dimension(nk), intent(in ) :: dz_top !< The HBD z grid [H ~> m or kg m-2]
710 type(hbd_cs), pointer :: CS !< Horizontal diffusion control structure
711
712 ! Local variables
713 real, allocatable :: phi_L_z(:) !< Tracer values in the ztop grid (left) [conc]
714 real, allocatable :: phi_R_z(:) !< Tracer values in the ztop grid (right) [conc]
715 real, allocatable :: F_layer_z(:) !< Diffusive flux at U/V-point in the ztop grid [H L2 conc ~> m3 conc]
716 real, allocatable :: khtr_ul_z(:) !< khtr_u at layer centers in the ztop grid [H L2 conc ~> m3 conc]
717 real, dimension(ke) :: h_vel !< Thicknesses at u- and v-points in the native grid
718 !! The harmonic mean is used to avoid zero values [H ~> m or kg m-2]
719 real, dimension(ke) :: khtr_ul !< khtr_u at the vertical layer of the native grid [L2 ~> m2]
720 real :: htot !< Total column thickness [H ~> m or kg m-2]
721 integer :: k !< Index used in the vertical direction
722 integer :: k_bot_min !< Minimum k-index for the bottom
723 integer :: k_bot_max !< Maximum k-index for the bottom
724 integer :: k_bot_diff !< Difference between bottom left and right k-indices
725 integer :: k_top_L, k_bot_L !< k-indices left native grid
726 integer :: k_top_R, k_bot_R !< k-indices right native grid
727 real :: zeta_top_L, zeta_top_R !< distance from the top of a layer to the boundary
728 !! layer depth in the native grid [nondim]
729 real :: zeta_bot_L, zeta_bot_R !< distance from the bottom of a layer to the boundary
730 !! layer depth in the native grid [nondim]
731 real :: wgt !< weight to be used in the linear transition to the interior [nondim]
732 real :: a !< coefficient used in the linear transition to the interior [nondim]
733 real :: tmp1, tmp2 !< dummy variables [H ~> m or kg m-2]
734 real :: htot_max !< depth below which no fluxes should be applied [H ~> m or kg m-2]
735
736 f_layer(:) = 0.0
737 khtr_ul(:) = 0.0
738 if (hbl_l == 0. .or. hbl_r == 0.) then
739 return
740 endif
741
742 ! allocate arrays
743 allocate(phi_l_z(nk), source=0.0)
744 allocate(phi_r_z(nk), source=0.0)
745 allocate(f_layer_z(nk), source=0.0)
746 allocate(khtr_ul_z(nk), source=0.0)
747
748 ! remap tracer to dz_top
749 call remapping_core_h(cs%remap_cs, ke, h_l(:), phi_l(:), nk, dz_top(:), phi_l_z(:))
750 call remapping_core_h(cs%remap_cs, ke, h_r(:), phi_r(:), nk, dz_top(:), phi_r_z(:))
751
752 ! thicknesses at velocity points & khtr_u at layer centers
753 do k = 1,ke
754 h_vel(k) = harmonic_mean(h_l(k), h_r(k))
755 ! GMM, writing 0.5 * (A(k) + A(k+1)) as A(k) + 0.5 * (A(k+1) - A(k)) to recover
756 ! answers with depth-independent khtr
757 khtr_ul(k) = khtr_u(k) + 0.5 * (khtr_u(k+1) - khtr_u(k))
758 enddo
759
760 ! remap khtr_ul to khtr_ul_z
761 call remapping_core_h(cs%remap_cs, ke, h_vel(:), khtr_ul(:), nk, dz_top(:), khtr_ul_z(:))
762
763 ! Calculate vertical indices containing the boundary layer in dz_top
764 call boundary_k_range(boundary, nk, dz_top, hbl_l, k_top_l, zeta_top_l, k_bot_l, zeta_bot_l)
765 call boundary_k_range(boundary, nk, dz_top, hbl_r, k_top_r, zeta_top_r, k_bot_r, zeta_bot_r)
766
767 if (boundary == surface) then
768 k_bot_min = min(k_bot_l, k_bot_r)
769 k_bot_max = max(k_bot_l, k_bot_r)
770 k_bot_diff = (k_bot_max - k_bot_min)
771
772 ! tracer flux where the minimum BLD intersects layer
773 if ((cs%linear) .and. (k_bot_diff > 1)) then
774 ! apply linear decay at the base of hbl
775 do k = k_bot_min,1,-1
776 f_layer_z(k) = -(dz_top(k) * khtr_ul_z(k)) * (phi_r_z(k) - phi_l_z(k))
777 if (cs%limiter_remap) call flux_limiter(f_layer_z(k), area_l, area_r, phi_l_z(k), &
778 phi_r_z(k), dz_top(k), dz_top(k))
779 enddo
780 htot = 0.0
781 do k = k_bot_min+1,k_bot_max, 1
782 htot = htot + dz_top(k)
783 enddo
784
785 a = -1.0/htot
786 htot = 0.
787 do k = k_bot_min+1,k_bot_max, 1
788 wgt = (a*(htot + (dz_top(k) * 0.5))) + 1.0
789 f_layer_z(k) = -(dz_top(k) * khtr_ul_z(k)) * (phi_r_z(k) - phi_l_z(k)) * wgt
790 htot = htot + dz_top(k)
791 if (cs%limiter_remap) call flux_limiter(f_layer_z(k), area_l, area_r, phi_l_z(k), &
792 phi_r_z(k), dz_top(k), dz_top(k))
793 enddo
794 else
795 do k = k_bot_min,1,-1
796 f_layer_z(k) = -(dz_top(k) * khtr_ul_z(k)) * (phi_r_z(k) - phi_l_z(k))
797 if (cs%limiter_remap) call flux_limiter(f_layer_z(k), area_l, area_r, phi_l_z(k), &
798 phi_r_z(k), dz_top(k), dz_top(k))
799 enddo
800 endif
801 endif
802
803 !GMM, TODO: boundary == BOTTOM
804
805 ! remap flux to h_vel (native grid)
806 call reintegrate_column(nk, dz_top(:), f_layer_z(:), ke, h_vel(:), f_layer(:))
807
808 ! used to avoid fluxes below hbl
809 if (cs%linear) then
810 htot_max = max(hbl_l, hbl_r)
811 else
812 htot_max = min(hbl_l, hbl_r)
813 endif
814
815 tmp1 = 0.0 ; tmp2 = 0.0
816 do k = 1,ke
817 ! apply flux_limiter
818 if (cs%limiter .and. f_layer(k) /= 0.) then
819 call flux_limiter(f_layer(k), area_l, area_r, phi_l(k), phi_r(k), h_l(k), h_r(k))
820 endif
821
822 ! if tracer point is below htot_max, set flux to zero
823 if (max(tmp1+(h_l(k)*0.5), tmp2+(h_r(k)*0.5)) > htot_max) then
824 f_layer(k) = 0.
825 endif
826
827 tmp1 = tmp1 + h_l(k)
828 tmp2 = tmp2 + h_r(k)
829 enddo
830
831 ! deallocated arrays
832 deallocate(phi_l_z)
833 deallocate(phi_r_z)
834 deallocate(f_layer_z)
835 deallocate(khtr_ul_z)
836
837end subroutine fluxes_layer_method
838
839!> Unit tests for near-boundary horizontal mixing
840logical function near_boundary_unit_tests( verbose )
841 logical, intent(in) :: verbose !< If true, output additional information for debugging unit tests
842
843 ! Local variables
844 integer, parameter :: nk = 2 ! Number of layers
845 real, dimension(nk+1) :: eta1 ! Updated interfaces with one extra value [m]
846 real, dimension(:), allocatable :: h1 ! Updated list of layer thicknesses or other field [m] or [arbitrary]
847 real, dimension(nk) :: phi_l, phi_r ! Tracer values (left and right column) [conc]
848 real, dimension(nk) :: h_l, h_r ! Layer thickness (left and right) [m]
849 real, dimension(nk+1) :: khtr_u ! Horizontal diffusivities at U-point and interfaces[m2 s-1]
850 real :: hbl_l, hbl_r ! Depth of the boundary layer (left and right) [m]
851 real, dimension(nk) :: f_layer ! Diffusive flux within each layer at U-point [conc m3 s-1]
852 character(len=120) :: test_name ! Title of the unit test
853 integer :: k_top ! Index of cell containing top of boundary
854 real :: zeta_top ! Fractional position in the cell of the top [nondim]
855 integer :: k_bot ! Index of cell containing bottom of boundary
856 real :: zeta_bot ! Fractional position in the cell of the bottom [nondim]
857 type(hbd_cs), pointer :: cs
858
859 allocate(cs)
860 ! fill required fields in CS
861 cs%linear=.false.
862 cs%H_subroundoff = 1.0e-20
863 cs%debug=.false.
864 cs%limiter=.false.
865 cs%limiter_remap=.false.
866 cs%hbd_nk = 2 + (2*2)
867 call initialize_remapping( cs%remap_CS, 'PLM', boundary_extrapolation=.true., &
868 om4_remap_via_sub_cells=.true., & ! ### see fail below when using fixed remapping alg.
869 check_reconstruction=.true., check_remapping=.true., &
870 h_neglect=cs%H_subroundoff, h_neglect_edge=cs%H_subroundoff)
871 call extract_member_remapping_cs(cs%remap_CS, degree=cs%deg)
872 allocate(cs%hbd_grd_u(1,1,cs%hbd_nk), source=0.0)
873 allocate(cs%hbd_u_kmax(1,1), source=0)
875 write(stdout,*) '==== MOM_hor_bnd_diffusion ======================='
876
877 ! Unit tests for boundary_k_range
878 test_name = 'Surface boundary spans the entire top cell'
879 h_l = (/5.,5./)
880 call boundary_k_range(surface, nk, h_l, 5., k_top, zeta_top, k_bot, zeta_bot)
882 test_boundary_k_range(k_top, zeta_top, k_bot, zeta_bot, 1, 0., 1, 1., test_name, verbose)
883
884 test_name = 'Surface boundary spans the entire column'
885 h_l = (/5.,5./)
886 call boundary_k_range(surface, nk, h_l, 10., k_top, zeta_top, k_bot, zeta_bot)
888 test_boundary_k_range(k_top, zeta_top, k_bot, zeta_bot, 1, 0., 2, 1., test_name, verbose)
889
890 test_name = 'Bottom boundary spans the entire bottom cell'
891 h_l = (/5.,5./)
892 call boundary_k_range(bottom, nk, h_l, 5., k_top, zeta_top, k_bot, zeta_bot)
894 test_boundary_k_range(k_top, zeta_top, k_bot, zeta_bot, 2, 1., 2, 0., test_name, verbose)
895
896 test_name = 'Bottom boundary spans the entire column'
897 h_l = (/5.,5./)
898 call boundary_k_range(bottom, nk, h_l, 10., k_top, zeta_top, k_bot, zeta_bot)
900 test_boundary_k_range(k_top, zeta_top, k_bot, zeta_bot, 1, 1., 2, 0., test_name, verbose)
901
902 test_name = 'Surface boundary intersects second layer'
903 h_l = (/10.,10./)
904 call boundary_k_range(surface, nk, h_l, 17.5, k_top, zeta_top, k_bot, zeta_bot)
906 test_boundary_k_range(k_top, zeta_top, k_bot, zeta_bot, 1, 0., 2, 0.75, test_name, verbose)
907
908 test_name = 'Surface boundary intersects first layer'
909 h_l = (/10.,10./)
910 call boundary_k_range(surface, nk, h_l, 2.5, k_top, zeta_top, k_bot, zeta_bot)
912 test_boundary_k_range(k_top, zeta_top, k_bot, zeta_bot, 1, 0., 1, 0.25, test_name, verbose)
913
914 test_name = 'Surface boundary is deeper than column thickness'
915 h_l = (/10.,10./)
916 call boundary_k_range(surface, nk, h_l, 21.0, k_top, zeta_top, k_bot, zeta_bot)
918 test_boundary_k_range(k_top, zeta_top, k_bot, zeta_bot, 1, 0., 2, 1., test_name, verbose)
919
920 test_name = 'Bottom boundary intersects first layer'
921 h_l = (/10.,10./)
922 call boundary_k_range(bottom, nk, h_l, 17.5, k_top, zeta_top, k_bot, zeta_bot)
924 test_boundary_k_range(k_top, zeta_top, k_bot, zeta_bot, 1, 0.75, 2, 0., test_name, verbose)
925
926 test_name = 'Bottom boundary intersects second layer'
927 h_l = (/10.,10./)
928 call boundary_k_range(bottom, nk, h_l, 2.5, k_top, zeta_top, k_bot, zeta_bot)
930 test_boundary_k_range(k_top, zeta_top, k_bot, zeta_bot, 2, 0.25, 2, 0., test_name, verbose)
931
932 if (.not. near_boundary_unit_tests) write(stdout,*) 'Passed boundary_k_range'
933
934 ! unit tests for sorting array and finding unique values
935 test_name = 'Sorting array'
936 eta1 = (/1., 0., 0.1/)
937 call sort(eta1, nk+1)
939 test_layer_fluxes( verbose, nk+1, test_name, eta1, (/0., 0.1, 1./) )
940
941 test_name = 'Unique values'
942 call unique((/0., 1., 1., 2./), nk+2, h1)
944 test_layer_fluxes( verbose, nk+1, test_name, h1, (/0., 1., 2./) )
945 deallocate(h1)
946
947 test_name = 'Unique values with maximum depth'
948 call unique((/0., 1., 1., 2., 3./), nk+3, h1, 2.)
950 test_layer_fluxes( verbose, nk+1, test_name, h1, (/0., 1., 2./) )
951 deallocate(h1)
952
953 if (.not. near_boundary_unit_tests) write(stdout,*) 'Passed sort and unique'
954
955 ! unit tests for merge_interfaces
956 test_name = 'h_L = h_R and BLD_L = BLD_R'
957 call merge_interfaces(nk, (/1., 2./), (/1., 2./), 1.5, 1.5, cs%H_subroundoff, h1)
959 test_layer_fluxes( verbose, nk, test_name, h1, (/1., 0.5/) )
960 deallocate(h1)
961
962 test_name = 'h_L = h_R and BLD_L /= BLD_R'
963 call merge_interfaces(nk, (/1., 2./), (/1., 2./), 0.5, 1.5, cs%H_subroundoff, h1)
965 test_layer_fluxes( verbose, nk+1, test_name, h1, (/0.5, 0.5, 0.5/) )
966 deallocate(h1)
967
968 test_name = 'h_L /= h_R and BLD_L = BLD_R'
969 call merge_interfaces(nk, (/1., 3./), (/2., 2./), 1.5, 1.5, cs%H_subroundoff, h1)
971 test_layer_fluxes( verbose, nk, test_name, h1, (/1., 0.5/) )
972 deallocate(h1)
973
974 test_name = 'h_L /= h_R and BLD_L /= BLD_R'
975 call merge_interfaces(nk, (/1., 3./), (/2., 2./), 0.5, 1.5, cs%H_subroundoff, h1)
977 test_layer_fluxes( verbose, nk+1, test_name, h1, (/0.5, 0.5, 0.5/) )
978 deallocate(h1)
979
980 test_name = 'Left deeper than right, h_L /= h_R and BLD_L /= BLD_R'
981 call merge_interfaces(nk, (/2., 3./), (/2., 2./), 1.0, 2.0, cs%H_subroundoff, h1)
983 test_layer_fluxes( verbose, nk, test_name, h1, (/1., 1./) )
984 deallocate(h1)
985
986 test_name = 'Left has zero thickness, h_L /= h_R and BLD_L = BLD_R'
987 call merge_interfaces(nk, (/4., 0./), (/2., 2./), 2.0, 2.0, cs%H_subroundoff, h1)
989 test_layer_fluxes( verbose, nk-1, test_name, h1, (/2./) )
990 deallocate(h1)
991
992 test_name = 'Left has zero thickness, h_L /= h_R and BLD_L /= BLD_R'
993 call merge_interfaces(nk, (/4., 0./), (/2., 2./), 1.0, 2.0, cs%H_subroundoff, h1)
995 test_layer_fluxes( verbose, nk, test_name, h1, (/1., 1./) )
996 deallocate(h1)
997
998 test_name = 'Right has zero thickness, h_L /= h_R and BLD_L = BLD_R'
999 call merge_interfaces(nk, (/2., 2./), (/0., 4./), 2.0, 2.0, cs%H_subroundoff, h1)
1001 test_layer_fluxes( verbose, nk-1, test_name, h1, (/2./) )
1002 deallocate(h1)
1003
1004 test_name = 'Right has zero thickness, h_L /= h_R and BLD_L /= BLD_R'
1005 call merge_interfaces(nk, (/2., 2./), (/0., 4./), 1.0, 2.0, cs%H_subroundoff, h1)
1007 test_layer_fluxes( verbose, nk, test_name, h1, (/1., 1./) )
1008 deallocate(h1)
1009
1010 test_name = 'Right deeper than left, h_L /= h_R and BLD_L = BLD_R'
1011 call merge_interfaces(nk+1, (/2., 2., 0./), (/2., 2., 1./), 4., 4., cs%H_subroundoff, h1)
1013 test_layer_fluxes( verbose, nk, test_name, h1, (/2., 2./) )
1014 deallocate(h1)
1015
1016 test_name = 'Right and left small values at bottom, h_L /= h_R and BLD_L = BLD_R'
1017 call merge_interfaces(nk+2, (/2., 2., 1., 1./), (/1., 1., .5, .5/), 3., 3., cs%H_subroundoff, h1)
1019 test_layer_fluxes( verbose, nk+2, test_name, h1, (/1., 1., .5, .5/) )
1020 deallocate(h1)
1021
1022 if (.not. near_boundary_unit_tests) write(stdout,*) 'Passed merge interfaces'
1023
1024 ! All cases in this section have hbl which are equal to the column thicknesses
1025 test_name = 'Equal hbl and same layer thicknesses (gradient from right to left)'
1026 hbl_l = 2. ; hbl_r = 2.
1027 h_l = (/2.,2./) ; h_r = (/2.,2./)
1028 phi_l = (/0.,0./) ; phi_r = (/1.,1./)
1029 khtr_u = (/1.,1.,1./)
1030 call hbd_grid_test(surface, hbl_l, hbl_r, h_l, h_r, cs)
1031 call fluxes_layer_method(surface, nk, hbl_l, hbl_r, h_l, h_r, phi_l, phi_r, &
1032 khtr_u, f_layer, 1., 1., cs%hbd_u_kmax(1,1), cs%hbd_grd_u(1,1,:), cs)
1034 test_layer_fluxes( verbose, nk, test_name, f_layer, (/-2.0,0.0/) )
1035
1036 test_name = 'Equal hbl and same layer thicknesses (gradient from left to right)'
1037 hbl_l = 2. ; hbl_r = 2.
1038 h_l = (/2.,2./) ; h_r = (/2.,2./)
1039 phi_l = (/2.,1./) ; phi_r = (/1.,1./)
1040 khtr_u = (/0.5,0.5,0.5/)
1041 call hbd_grid_test(surface, hbl_l, hbl_r, h_l, h_r, cs)
1042 call fluxes_layer_method(surface, nk, hbl_l, hbl_r, h_l, h_r, phi_l, phi_r, &
1043 khtr_u, f_layer, 1., 1., cs%hbd_u_kmax(1,1), cs%hbd_grd_u(1,1,:), cs)
1045 test_layer_fluxes( verbose, nk, test_name, f_layer, (/1.0,0.0/) )
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 khtr_u = (/2.,2.,2./)
1052 call hbd_grid_test(surface, hbl_l, hbl_r, h_l, h_r, cs)
1053 call fluxes_layer_method(surface, nk, hbl_l, hbl_r, h_l, h_r, phi_l, phi_r, &
1054 khtr_u, f_layer, 1., 1., cs%hbd_u_kmax(1,1), cs%hbd_grd_u(1,1,:), cs)
1055 ! ### This test fails when om4_remap_via_sub_cells=.false.
1057 test_layer_fluxes( verbose, nk, test_name, f_layer, (/-1.0,-4.0/) )
1058
1059 test_name = 'Different hbl and different column thicknesses (zero gradient)'
1060 hbl_l = 12 ; hbl_r = 20
1061 h_l = (/6.,6./) ; h_r = (/10.,10./)
1062 phi_l = (/1.,1./) ; phi_r = (/1.,1./)
1063 khtr_u = (/1.,1.,1./)
1064 call hbd_grid_test(surface, hbl_l, hbl_r, h_l, h_r, cs)
1065 call fluxes_layer_method(surface, nk, hbl_l, hbl_r, h_l, h_r, phi_l, phi_r, &
1066 khtr_u, f_layer, 1., 1., cs%hbd_u_kmax(1,1), cs%hbd_grd_u(1,1,:), cs)
1068 test_layer_fluxes( verbose, nk, test_name, f_layer, (/0.,0./) )
1069
1070 test_name = 'Different hbl and different column thicknesses (gradient from left to right)'
1071
1072 hbl_l = 15 ; hbl_r = 10.
1073 h_l = (/10.,5./) ; h_r = (/10.,0./)
1074 phi_l = (/1.,1./) ; phi_r = (/0.,0./)
1075 khtr_u = (/1.,1.,1./)
1076 call hbd_grid_test(surface, hbl_l, hbl_r, h_l, h_r, cs)
1077 call fluxes_layer_method(surface, nk, hbl_l, hbl_r, h_l, h_r, phi_l, phi_r, &
1078 khtr_u, f_layer, 1., 1., cs%hbd_u_kmax(1,1), cs%hbd_grd_u(1,1,:), cs)
1080 test_layer_fluxes( verbose, nk, test_name, f_layer, (/10.,0.0/) )
1081
1082 if (.not. near_boundary_unit_tests) write(stdout,*) 'Passed fluxes_layer_method'
1083
1084end function near_boundary_unit_tests
1085
1086!> Returns true if output of near-boundary unit tests does not match correct computed values
1087!! and conditionally writes results to stream
1088logical function test_layer_fluxes(verbose, nk, test_name, F_calc, F_ans)
1089 logical, intent(in) :: verbose !< If true, write results to stdout
1090 character(len=80), intent(in) :: test_name !< Brief description of the unit test
1091 integer, intent(in) :: nk !< Number of layers
1092 real, dimension(nk), intent(in) :: f_calc !< Fluxes or other quantity from the algorithm [arbitrary]
1093 real, dimension(nk), intent(in) :: f_ans !< Expected value calculated by hand [arbitrary]
1094 ! Local variables
1095 integer :: k
1096
1097 test_layer_fluxes = .false.
1098 do k=1,nk
1099 if ( f_calc(k) /= f_ans(k) ) then
1100 test_layer_fluxes = .true.
1101 write(stdout,*) "MOM_hor_bnd_diffusion, UNIT TEST FAILED: ", test_name
1102 write(stdout,10) k, f_calc(k), f_ans(k)
1103 elseif (verbose) then
1104 write(stdout,10) k, f_calc(k), f_ans(k)
1105 endif
1106 enddo
1107
110810 format("Layer=",i3," F_calc=",f20.16," F_ans",f20.16)
1109end function test_layer_fluxes
1110
1111!> Return true if output of unit tests for boundary_k_range does not match answers
1112logical function test_boundary_k_range(k_top, zeta_top, k_bot, zeta_bot, k_top_ans, zeta_top_ans,&
1113 k_bot_ans, zeta_bot_ans, test_name, verbose)
1114 integer :: k_top !< Index of cell containing top of boundary
1115 real :: zeta_top !< Fractional position in the cell of the top boundary [nondim]
1116 integer :: k_bot !< Index of cell containing bottom of boundary
1117 real :: zeta_bot !< Fractional position in the cell of the bottom boundary [nondim]
1118 integer :: k_top_ans !< Expected index of cell containing top of boundary
1119 real :: zeta_top_ans !< Expected fractional position of the top boundary [nondim]
1120 integer :: k_bot_ans !< Expected index of cell containing bottom of boundary
1121 real :: zeta_bot_ans !< Expected fractional position of the bottom boundary [nondim]
1122 character(len=80) :: test_name !< Name of the unit test
1123 logical :: verbose !< If true always print output
1124
1125 test_boundary_k_range = k_top /= k_top_ans
1126 test_boundary_k_range = test_boundary_k_range .or. (zeta_top /= zeta_top_ans)
1127 test_boundary_k_range = test_boundary_k_range .or. (k_bot /= k_bot_ans)
1128 test_boundary_k_range = test_boundary_k_range .or. (zeta_bot /= zeta_bot_ans)
1129
1130 if (test_boundary_k_range) write(stdout,*) "UNIT TEST FAILED: ", test_name
1131 if (test_boundary_k_range .or. verbose) then
1132 write(stdout,20) "k_top", k_top, "k_top_ans", k_top_ans
1133 write(stdout,20) "k_bot", k_bot, "k_bot_ans", k_bot_ans
1134 write(stdout,30) "zeta_top", zeta_top, "zeta_top_ans", zeta_top_ans
1135 write(stdout,30) "zeta_bot", zeta_bot, "zeta_bot_ans", zeta_bot_ans
1136 endif
1137
1138 20 format(a,"=",i3,1x,a,"=",i3)
1139 30 format(a,"=",f20.16,1x,a,"=",f20.16)
1140
1141
1142end function test_boundary_k_range
1143
1144!> Same as hbd_grid, but only used in the unit tests.
1145subroutine hbd_grid_test(boundary, hbl_L, hbl_R, h_L, h_R, CS)
1146 integer, intent(in) :: boundary !< Which boundary layer SURFACE or BOTTOM [nondim]
1147 real, intent(in) :: hbl_L !< Boundary layer depth, left [H ~> m or kg m-2]
1148 real, intent(in) :: hbl_R !< Boundary layer depth, right [H ~> m or kg m-2]
1149 real, dimension(2), intent(in) :: h_L !< Layer thickness in the native grid, left [H ~> m or kg m-2]
1150 real, dimension(2), intent(in) :: h_R !< Layer thickness in the native grid, right [H ~> m or kg m-2]
1151 type(hbd_cs), pointer :: CS !< Horizontal diffusion control structure
1152
1153 ! Local variables
1154 real, allocatable :: dz_top(:) !< temporary HBD grid given by merge_interfaces [H ~> m or kg m-2]
1155 integer :: nk, k !< number of layers in the HBD grid, and integers used in do-loops
1156
1157 ! reset arrays
1158 cs%hbd_grd_u(1,1,:) = 0.0
1159 cs%hbd_u_kmax(1,1) = 0
1160
1161 call merge_interfaces(2, h_l, h_r, hbl_l, hbl_r, cs%H_subroundoff, dz_top)
1162 nk = SIZE(dz_top)
1163 if (nk > cs%hbd_nk) then
1164 write(*,*)'nk, CS%hbd_nk', nk, cs%hbd_nk
1165 call mom_error(fatal,"Houston, we've had a problem in hbd_grid_test, (nk cannot be > CS%hbd_nk)")
1166 endif
1167
1168 cs%hbd_u_kmax(1,1) = nk
1169
1170 ! set the HBD grid to dz_top
1171 do k=1,nk
1172 cs%hbd_grd_u(1,1,k) = dz_top(k)
1173 enddo
1174 deallocate(dz_top)
1175
1176end subroutine hbd_grid_test
1177
1178!> Deallocates hor_bnd_diffusion control structure
1179subroutine hor_bnd_diffusion_end(CS)
1180 type(hbd_cs), pointer :: cs !< Horizontal boundary diffusion control structure
1181
1182 if (associated(cs)) deallocate(cs)
1183
1184end subroutine hor_bnd_diffusion_end
1185
1186!> \namespace mom_hor_bnd_diffusion
1187!!
1188!! \section section_HBD The Horizontal Boundary Diffusion (HBD) framework
1189!!
1190!! The HBD framework accounts for the effects of diabatic mesoscale fluxes
1191!! within surface and bottom boundary layers. Unlike the equivalent adiabatic
1192!! fluxes, which is applied along neutral density surfaces, HBD is purely
1193!! horizontal. To assure that diffusive fluxes are strictly horizontal
1194!! regardless of the vertical coordinate system, this method relies on
1195!! regridding/remapping techniques.
1196!!
1197!! The bottom boundary layer fluxes remain to be implemented, although some
1198!! of the steps needed to do so have already been added and tested.
1199!!
1200!! Horizontal boundary diffusion is applied as follows:
1201!!
1202!! 1) remap tracer to a z* grid (HBD grid)
1203!! 2) calculate diffusive tracer fluxes (F) in the HBD grid using a layer by layer approach (@ref section_method)
1204!! 3) remap fluxes to the native grid
1205!! 4) update tracer by adding the divergence of F
1206!!
1207!! \subsection section_method Along layer approach
1208!!
1209!! Here diffusion is applied layer by layer using only information from neighboring cells.
1210!!
1211!! Step #1: define vertical grid using interfaces and surface boundary layers from left and right
1212!! columns (see merge_interfaces).
1213!!
1214!! Step #2: compute vertical indices containing boundary layer (boundary_k_range).
1215!! For the TOP boundary layer, these are:
1216!!
1217!! k_top, k_bot, zeta_top, zeta_bot
1218!!
1219!! Step #2: calculate the diffusive flux at each layer:
1220!!
1221!! \f[ F_{k} = -KHTR \times h_{eff}(k) \times (\phi_R(k) - \phi_L(k)), \f]
1222!! where h_eff is the [harmonic mean](@ref section_harmonic_mean) of the layer thickness
1223!! in the left and right columns.
1224!!
1225!! Step #3: option to linearly decay the flux from k_bot_min to k_bot_max:
1226!!
1227!! If HBD_LINEAR_TRANSITION = True and k_bot_diff > 1, the diffusive flux will decay
1228!! linearly between the top interface of the layer containing the minimum boundary
1229!! layer depth (k_bot_min) and the lower interface of the layer containing the
1230!! maximum layer depth (k_bot_max).
1231!!
1232!! Step #4: remap the fluxes back to the native grid. This is done at velocity points, whose vertical grid
1233!! is determined using [harmonic mean](@ref section_harmonic_mean). To assure monotonicity,
1234!! tracer fluxes are limited so that 1) only down-gradient fluxes are applied,
1235!! and 2) the flux cannot be larger than F_max, which is defined using the tracer
1236!! gradient:
1237!!
1238!! \f[ F_{max} = -0.2 \times [(V_R(k) \times \phi_R(k)) - (V_L(k) \times \phi_L(k))], \f]
1239!! where V is the cell volume. Why 0.2?
1240!! t=0 t=inf
1241!! 0 .2
1242!! 0 1 0 .2.2.2
1243!! 0 .2
1244!!
1245!! \subsection section_harmonic_mean Harmonic Mean
1246!!
1247!! The harmonic mean (HM) between h1 and h2 is defined as:
1248!!
1249!! \f[ HM = \frac{2 \times h1 \times h2}{h1 + h2} \f]
1250!!
1251end module mom_hor_bnd_diffusion