← back to index

src/tracer/MOM_hor_bnd_diffusion.F90

portedportable, not yet portedexecuted, not portableexecutable, not hit by this run

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
8module MOM_hor_bnd_diffusion
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
19use MOM_remapping, only : remapping_CS, initialize_remapping, reintegrate_column
20use MOM_remapping, only : extract_member_remapping_CS, remapping_core_h
21use MOM_remapping, only : remappingSchemesDoc, remappingDefaultScheme
22use MOM_spatial_means, only : global_mass_integral
23use MOM_tracer_registry, only : tracer_registry_type, tracer_type
24use MOM_unit_scaling, only : unit_scale_type
25use MOM_variables, only : vertvisc_type
26use MOM_verticalGrid, only : verticalGrid_type
27use MOM_CVMix_KPP, only : KPP_get_BLD, KPP_CS
28use MOM_energetic_PBL, only : energetic_PBL_get_MLD, energetic_PBL_CS
29use MOM_diabatic_driver, only : diabatic_CS, extract_diabatic_member
30use MOM_io, only : stdout, stderr
31
32implicit none ; private
33
34public near_boundary_unit_tests, hor_bnd_diffusion, hor_bnd_diffusion_init
35public boundary_k_range, hor_bnd_diffusion_end
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.
811logical 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
971 if (ASSOCIATED(CS)) then
980 call MOM_error(FATAL, "hor_bnd_diffusion_init called with associated control structure.")
990 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, &
1041 default=.false., do_not_log=.true.)
105 call log_version(param_file, mdl, version, &
106 "This module implements horizontal diffusion of tracers near boundaries", &
1071 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.", &
1101 default=.false.)
1111 if (.not. hor_bnd_diffusion_init) return
112
1130 allocate(CS)
1140 CS%diag => diag
1150 CS%H_subroundoff = GV%H_subroundoff
1160 call extract_diabatic_member(diabatic_CSp, KPP_CSp=CS%KPP_CSp)
1170 call extract_diabatic_member(diabatic_CSp, energetic_PBL_CSp=CS%energetic_PBL_CSp)
118
119 ! max. number of vertical layers
1200 CS%hbd_nk = 2 + (GV%ke*2)
121 ! allocate the hbd grids and k_max
1220 allocate(CS%hbd_grd_u(SZIB_(G),SZJ_(G),CS%hbd_nk), source=0.0)
1230 allocate(CS%hbd_grd_v(SZI_(G),SZJB_(G),CS%hbd_nk), source=0.0)
1240 allocate(CS%hbd_u_kmax(SZIB_(G),SZJ_(G)), source=0)
1250 allocate(CS%hbd_v_kmax(SZI_(G),SZJB_(G)), source=0)
126
1270 CS%surface_boundary_scheme = -1
1280 if ( .not. ASSOCIATED(CS%energetic_PBL_CSp) .and. .not. ASSOCIATED(CS%KPP_CSp) ) then
1290 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"//&
1350 "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, &
1370 "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, &
1390 "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", &
1420 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: "//&
1470 trim(remappingSchemesDoc), default=remappingDefaultScheme)
148 call get_param(param_file, mdl, "REMAPPING_USE_OM4_SUBCELLS", om4_remap_via_sub_cells, &
1490 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. "//&
1540 "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., &
1600 h_neglect=CS%H_subroundoff, h_neglect_edge=CS%H_subroundoff)
1610 call extract_member_remapping_CS(CS%remap_CS, degree=CS%deg)
162 call get_param(param_file, mdl, "DEBUG", debug, &
1630 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.", &
1660 default=debug, debuggingParam=.true.)
167
1680 id_clock_hbd = cpu_clock_id('(Ocean HBD)', grain=CLOCK_MODULE)
169
1700end 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
1780subroutine 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
1930 real, dimension(SZI_(G),SZJ_(G)) :: hbl !< Boundary layer depth [H ~> m or kg m-2]
1940 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)) :: uFlx !< Zonal flux of tracer [conc H L2 ~> conc m3 or conc kg]
1950 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)) :: vFlx !< Meridional flux of tracer
196 !! [conc H L2 ~> conc m3 or conc kg]
1970 real, dimension(SZIB_(G),SZJ_(G)) :: uwork_2d !< Layer summed u-flux transport
198 !! [conc H L2 ~> conc m3 or conc kg]
1990 real, dimension(SZI_(G),SZJB_(G)) :: vwork_2d !< Layer summed v-flux transport
200 !! [conc H L2 ~> conc m3 or conc kg]
2010 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].
2070 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]
2120 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
2200 call cpu_clock_begin(id_clock_hbd)
2210 Idt = 1./dt
222
2230 if (associated(visc%h_ML)) then
2240 hbl(:,:) = visc%h_ML(:,:)
225 else
2260 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.
2290 call pass_var(hbl, G%Domain, halo=1)
230
231 ! build HBD grid
2320 call hbd_grid(SURFACE, G, GV, hbl, h, CS)
233
2340 do m = 1,Reg%ntr
235 ! current tracer
2360 tracer => Reg%tr(m)
237
2380 if (CS%debug) then
2390 call hchksum(tracer%t, "before HBD "//tracer%name, G%HI, unscale=tracer%conc_scale)
240 endif
241
242 ! for diagnostics
2430 if (tracer%id_hbdxy_conc > 0 .or. tracer%id_hbdxy_cont > 0 .or. tracer%id_hbdxy_cont_2d > 0 .or. CS%debug) then
2440 tendency(:,:,:) = 0.0
2450 tracer_old(:,:,:) = tracer%t(:,:,:)
246 endif
247
248 ! Diffusive fluxes in the i- and j-direction
2490 uFlx(:,:,:) = 0.
2500 vFlx(:,:,:) = 0.
251
252 ! HBD layer by layer
2530 do j=G%jsc,G%jec
2540 do i=G%isc-1,G%iec
2550 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), &
2590 CS%hbd_grd_u(I,j,:), CS)
260 endif
261 enddo
262 enddo
2630 do J=G%jsc-1,G%jec
2640 do i=G%isc,G%iec
2650 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), &
2690 CS%hbd_grd_v(i,J,:), CS)
270 endif
271 enddo
272 enddo
273
274 ! Update the tracer fluxes
2750 do k=1,GV%ke ; do j=G%jsc,G%jec ; do i=G%isc,G%iec
2760 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) ) ))* &
2780 G%IareaT(i,j) / ( h(i,j,k) + GV%H_subroundoff )
279
2800 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))) * &
2820 G%IareaT(i,j) * Idt
283 endif
284 endif
285 enddo ; enddo ; enddo
286
287 ! Do user controlled underflow of the tracer concentrations.
2880 if (tracer%conc_underflow > 0.0) then
2890 do k=1,GV%ke ; do j=G%jsc,G%jec ; do i=G%isc,G%iec
2900 if (abs(tracer%t(i,j,k)) < tracer%conc_underflow) tracer%t(i,j,k) = 0.0
291 enddo ; enddo ; enddo
292 endif
293
2940 if (CS%debug) then
2950 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
2970 tracer_int_prev = global_mass_integral(h, G, GV, tracer_old, scale=tracer%conc_scale)
2980 tracer_int_end = global_mass_integral(h, G, GV, tracer%t, scale=tracer%conc_scale)
2990 write(mesg,*) 'Total '//tracer%name//' before/after HBD:', tracer_int_prev, tracer_int_end
3000 call MOM_mesg(mesg)
301 endif
302
303 ! Post the tracer diagnostics
3040 if (tracer%id_hbd_dfx>0) call post_data(tracer%id_hbd_dfx, uFlx(:,:,:)*Idt, CS%diag)
3050 if (tracer%id_hbd_dfy>0) call post_data(tracer%id_hbd_dfy, vFlx(:,:,:)*Idt, CS%diag)
3060 if (tracer%id_hbd_dfx_2d>0) then
3070 uwork_2d(:,:) = 0.
3080 do k=1,GV%ke ; do j=G%jsc,G%jec ; do I=G%isc-1,G%iec
3090 uwork_2d(I,j) = uwork_2d(I,j) + (uFlx(I,j,k) * Idt)
310 enddo ; enddo ; enddo
3110 call post_data(tracer%id_hbd_dfx_2d, uwork_2d, CS%diag)
312 endif
313
3140 if (tracer%id_hbd_dfy_2d>0) then
3150 vwork_2d(:,:) = 0.
3160 do k=1,GV%ke ; do J=G%jsc-1,G%jec ; do i=G%isc,G%iec
3170 vwork_2d(i,J) = vwork_2d(i,J) + (vFlx(i,J,k) * Idt)
318 enddo ; enddo ; enddo
3190 call post_data(tracer%id_hbd_dfy_2d, vwork_2d, CS%diag)
320 endif
321
322 ! post tendency of tracer content
3230 if (tracer%id_hbdxy_cont > 0) then
3240 call post_data(tracer%id_hbdxy_cont, tendency, CS%diag)
325 endif
326
327 ! post depth summed tendency for tracer content
3280 if (tracer%id_hbdxy_cont_2d > 0) then
3290 tendency_2d(:,:) = 0.
3300 do j=G%jsc,G%jec ; do i=G%isc,G%iec
3310 do k=1,GV%ke
3320 tendency_2d(i,j) = tendency_2d(i,j) + tendency(i,j,k)
333 enddo
334 enddo ; enddo
3350 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.
3410 if (tracer%id_hbdxy_conc > 0) then
3420 do k=1,GV%ke ; do j=G%jsc,G%jec ; do i=G%isc,G%iec
3430 tendency(i,j,k) = tendency(i,j,k) / ( h(i,j,k) + CS%H_subroundoff )
344 enddo ; enddo ; enddo
3450 call post_data(tracer%id_hbdxy_conc, tendency, CS%diag)
346 endif
347
348 enddo
349
3500 call cpu_clock_end(id_clock_hbd)
351
3520end subroutine hor_bnd_diffusion
353
354!> Build the HBD grid where tracers will be remapped to.
3550subroutine 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
3660 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
3700 CS%hbd_grd_u(:,:,:) = 0.0
3710 CS%hbd_grd_v(:,:,:) = 0.0
3720 CS%hbd_u_kmax(:,:) = 0
3730 CS%hbd_v_kmax(:,:) = 0
374
3750 do j=G%jsc,G%jec
3760 do I=G%isc-1,G%iec
3770 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), &
3790 CS%H_subroundoff, dz_top)
3800 nk = SIZE(dz_top)
3810 if (nk > CS%hbd_nk) then
3820 write(*,*)'nk, CS%hbd_nk', nk, CS%hbd_nk
3830 call MOM_error(FATAL,"Houston, we've had a problem in hbd_grid, u-points (nk cannot be > CS%hbd_nk)")
384 endif
385
3860 CS%hbd_u_kmax(I,j) = nk
387
388 ! set the HBD grid to dz_top
3890 do k=1,nk
3900 CS%hbd_grd_u(I,j,k) = dz_top(k)
391 enddo
3920 deallocate(dz_top)
393 endif
394 enddo
395 enddo
396
3970 do J=G%jsc-1,G%jec
3980 do i=G%isc,G%iec
3990 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), &
4010 CS%H_subroundoff, dz_top)
402
4030 nk = SIZE(dz_top)
4040 if (nk > CS%hbd_nk) then
4050 write(*,*)'nk, CS%hbd_nk', nk, CS%hbd_nk
4060 call MOM_error(FATAL,"Houston, we've had a problem in hbd_grid, v-points (nk cannot be > CS%hbd_nk)")
407 endif
408
4090 CS%hbd_v_kmax(i,J) = nk
410
411 ! set the HBD grid to dz_top
4120 do k=1,nk
4130 CS%hbd_grd_v(i,J,k) = dz_top(k)
414 enddo
4150 deallocate(dz_top)
416 endif
417 enddo
418 enddo
419
4200end subroutine hbd_grid
421
422!> Calculate the harmonic mean of two quantities [arbitrary]
423!! See \ref section_harmonic_mean.
4240real function harmonic_mean(h1,h2)
425 real :: h1 !< Scalar quantity [arbitrary]
426 real :: h2 !< Scalar quantity [arbitrary]
4270 if (h1 + h2 == 0.) then
4280 harmonic_mean = 0.
429 else
4300 harmonic_mean = 2.*(h1*h2)/(h1+h2)
431 endif
4320end function harmonic_mean
433
434!> Returns the location of the minimum value in a 1D array
435!! between indices s and e.
4360integer 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
4460 minimum = x(s) ! assume the first is the min
4470 location = s ! record its position
4480 do i = s+1, e ! start with next elements
4490 if (x(i) < minimum) then ! if x(i) less than the min?
4500 minimum = x(i) ! Yes, a new minimum found
4510 location = i ! record its position
452 endif
453 enddo
4540 find_minimum = location ! return the position
4550end function find_minimum
456
457!> Swaps the values of its two formal arguments.
4580subroutine 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
4650 tmp = a
4660 a = b
4670 b = tmp
4680end subroutine swap
469
470!> Receives a 1D array x and sorts it into ascending order.
4710subroutine 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
4780 do i = 1, n-1
4790 location = find_minimum(x, i, n) ! find min from this to last
4800 call swap(x(i), x(location)) ! swap this and the minimum
481 enddo
4820end subroutine sort
483
484!> Returns the unique values in a 1D array.
4850subroutine 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
4920 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
4970 limit = .false.
4980 if (present(val_max)) then
4990 limit = .true.
5000 if (val_max > MAXVAL(val)) then
5010 if (is_root_pe()) write(*,*)'val_max, MAXVAL(val)',val_max, MAXVAL(val)
5020 call MOM_error(FATAL,"Houston, we've had a problem in unique (val_max cannot be > MAXVAL(val))")
503 endif
504 endif
505
5060 tmp(:) = 0.
5070 min_val = MINVAL(val)-1
5080 max_val = MAXVAL(val)
5090 i = 0
5100 do while (min_val<max_val)
5110 i = i+1
5120 min_val = MINVAL(val, mask=val>min_val)
5130 tmp(i) = min_val
514 enddo
5150 ii = i
5160 if (limit) then
5170 do j=1,ii
5180 if (tmp(j) <= val_max) i = j
519 enddo
520 endif
5210 allocate(val_unique(i), source=tmp(1:i))
5220end 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.
5290subroutine 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
5420 real, dimension(nk+1) :: eta_L, eta_R!< Interfaces in the left and right columns [H ~> m or kg m-2]
5430 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]
5450 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
5520 n = (2*nk)+3
5530 allocate(eta_all(n))
554 ! compute and merge interfaces
5550 eta_L(:) = 0.0 ; eta_R(:) = 0.0 ; eta_all(:) = 0.0
5560 kk = 0
5570 do k=2,nk+1
5580 eta_L(k) = eta_L(k-1) + h_L(k-1)
5590 eta_R(k) = eta_R(k-1) + h_R(k-1)
5600 kk = kk + 2
5610 eta_all(kk) = eta_L(k)
5620 eta_all(kk+1) = eta_R(k)
563 enddo
564
565 ! add hbl_L and hbl_R into eta_all
5660 eta_all(kk+2) = hbl_L
5670 eta_all(kk+3) = hbl_R
568
569 ! find maximum depth
5700 min_depth = MIN(MAXVAL(eta_L), MAXVAL(eta_R))
5710 max_bld = MAX(hbl_L, hbl_R)
5720 max_depth = MIN(min_depth, max_bld)
573
574 ! sort eta_all
5750 call sort(eta_all, n)
576 ! remove duplicates from eta_all and sets maximum depth
5770 call unique(eta_all, n, eta_unique, max_depth)
578
5790 nk1 = SIZE(eta_unique)
5800 allocate(h(nk1-1))
5810 do k=1,nk1-1
5820 h(k) = (eta_unique(k+1) - eta_unique(k)) + H_subroundoff
583 enddo
5840end 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.
5880subroutine 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 !
6060 F_max = -0.2 * ((area_R*(phi_R*h_R))-(area_L*(phi_L*h_L)))
607
6080 if ( SIGN(1.,F_layer) == SIGN(1., F_max)) then
609 ! Apply flux limiter calculated above
6100 if (F_max >= 0.) then
6110 F_layer = MIN(F_layer,F_max)
612 else
6130 F_layer = MAX(F_layer,F_max)
614 endif
615 else
6160 F_layer = 0.0
617 endif
6180end subroutine flux_limiter
619
620!> Find the k-index range corresponding to the layers that are within the boundary-layer region
6210subroutine 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
6390 if ( boundary == SURFACE ) then
6400 k_top = 1
6410 zeta_top = 0.
6420 htot = 0.
6430 k_bot = 1
6440 zeta_bot = 0.
6450 if (hbl == 0.) return
6460 if (hbl >= SUM(h(:))) then
6470 k_bot = nk
6480 zeta_bot = 1.
6490 return
650 endif
6510 do k=1,nk
6520 htot = htot + h(k)
6530 if ( htot >= hbl) then
6540 k_bot = k
6550 zeta_bot = 1 - (htot - hbl)/h(k)
6560 return
657 endif
658 enddo
659
660 ! Bottom boundary layer
6610 elseif ( boundary == BOTTOM ) then
6620 k_top = nk
6630 zeta_top = 1.
6640 k_bot = nk
6650 zeta_bot = 0.
6660 htot = 0.
6670 if (hbl == 0.) return
6680 if (hbl >= SUM(h(:))) then
6690 k_top = 1
6700 zeta_top = 1.
6710 return
672 endif
6730 do k=nk,1,-1
6740 htot = htot + h(k)
6750 if (htot >= hbl) then
6760 k_top = k
6770 zeta_top = 1 - (htot - hbl)/h(k)
6780 return
679 endif
680 enddo
681 else
6820 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
6890subroutine fluxes_layer_method(boundary, ke, hbl_L, hbl_R, h_L, h_R, phi_L, phi_R, &
6900 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
7130 real, allocatable :: phi_L_z(:) !< Tracer values in the ztop grid (left) [conc]
7140 real, allocatable :: phi_R_z(:) !< Tracer values in the ztop grid (right) [conc]
7150 real, allocatable :: F_layer_z(:) !< Diffusive flux at U/V-point in the ztop grid [H L2 conc ~> m3 conc]
7160 real, allocatable :: khtr_ul_z(:) !< khtr_u at layer centers in the ztop grid [H L2 conc ~> m3 conc]
7170 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]
7190 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
7360 F_layer(:) = 0.0
7370 khtr_ul(:) = 0.0
7380 if (hbl_L == 0. .or. hbl_R == 0.) then
7390 return
740 endif
741
742 ! allocate arrays
7430 allocate(phi_L_z(nk), source=0.0)
7440 allocate(phi_R_z(nk), source=0.0)
7450 allocate(F_layer_z(nk), source=0.0)
7460 allocate(khtr_ul_z(nk), source=0.0)
747
748 ! remap tracer to dz_top
7490 call remapping_core_h(CS%remap_cs, ke, h_L(:), phi_L(:), nk, dz_top(:), phi_L_z(:))
7500 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
7530 do k = 1,ke
7540 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
7570 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
7610 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
7640 call boundary_k_range(boundary, nk, dz_top, hbl_L, k_top_L, zeta_top_L, k_bot_L, zeta_bot_L)
7650 call boundary_k_range(boundary, nk, dz_top, hbl_R, k_top_R, zeta_top_R, k_bot_R, zeta_bot_R)
766
7670 if (boundary == SURFACE) then
7680 k_bot_min = MIN(k_bot_L, k_bot_R)
7690 k_bot_max = MAX(k_bot_L, k_bot_R)
7700 k_bot_diff = (k_bot_max - k_bot_min)
771
772 ! tracer flux where the minimum BLD intersects layer
7730 if ((CS%linear) .and. (k_bot_diff > 1)) then
774 ! apply linear decay at the base of hbl
7750 do k = k_bot_min,1,-1
7760 F_layer_z(k) = -(dz_top(k) * khtr_ul_z(k)) * (phi_R_z(k) - phi_L_z(k))
7770 if (CS%limiter_remap) call flux_limiter(F_layer_z(k), area_L, area_R, phi_L_z(k), &
7780 phi_R_z(k), dz_top(k), dz_top(k))
779 enddo
7800 htot = 0.0
7810 do k = k_bot_min+1,k_bot_max, 1
7820 htot = htot + dz_top(k)
783 enddo
784
7850 a = -1.0/htot
7860 htot = 0.
7870 do k = k_bot_min+1,k_bot_max, 1
7880 wgt = (a*(htot + (dz_top(k) * 0.5))) + 1.0
7890 F_layer_z(k) = -(dz_top(k) * khtr_ul_z(k)) * (phi_R_z(k) - phi_L_z(k)) * wgt
7900 htot = htot + dz_top(k)
7910 if (CS%limiter_remap) call flux_limiter(F_layer_z(k), area_L, area_R, phi_L_z(k), &
7920 phi_R_z(k), dz_top(k), dz_top(k))
793 enddo
794 else
7950 do k = k_bot_min,1,-1
7960 F_layer_z(k) = -(dz_top(k) * khtr_ul_z(k)) * (phi_R_z(k) - phi_L_z(k))
7970 if (CS%limiter_remap) call flux_limiter(F_layer_z(k), area_L, area_R, phi_L_z(k), &
7980 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)
8060 call reintegrate_column(nk, dz_top(:), F_layer_z(:), ke, h_vel(:), F_layer(:))
807
808 ! used to avoid fluxes below hbl
8090 if (CS%linear) then
8100 htot_max = MAX(hbl_L, hbl_R)
811 else
8120 htot_max = MIN(hbl_L, hbl_R)
813 endif
814
8150 tmp1 = 0.0 ; tmp2 = 0.0
8160 do k = 1,ke
817 ! apply flux_limiter
8180 if (CS%limiter .and. F_layer(k) /= 0.) then
8190 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
8230 if (MAX(tmp1+(h_L(k)*0.5), tmp2+(h_R(k)*0.5)) > htot_max) then
8240 F_layer(k) = 0.
825 endif
826
8270 tmp1 = tmp1 + h_L(k)
8280 tmp2 = tmp2 + h_R(k)
829 enddo
830
831 ! deallocated arrays
8320 deallocate(phi_L_z)
8330 deallocate(phi_R_z)
8340 deallocate(F_layer_z)
8350 deallocate(khtr_ul_z)
836
8370end subroutine fluxes_layer_method
838
839!> Unit tests for near-boundary horizontal mixing
8400logical 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]
8460 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
8590 allocate(CS)
860 ! fill required fields in CS
8610 CS%linear=.false.
8620 CS%H_subroundoff = 1.0E-20
8630 CS%debug=.false.
8640 CS%limiter=.false.
8650 CS%limiter_remap=.false.
8660 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., &
8700 h_neglect=CS%H_subroundoff, h_neglect_edge=CS%H_subroundoff)
8710 call extract_member_remapping_CS(CS%remap_CS, degree=CS%deg)
8720 allocate(CS%hbd_grd_u(1,1,CS%hbd_nk), source=0.0)
8730 allocate(CS%hbd_u_kmax(1,1), source=0)
8740 near_boundary_unit_tests = .false.
8750 write(stdout,*) '==== MOM_hor_bnd_diffusion ======================='
876
877 ! Unit tests for boundary_k_range
8780 test_name = 'Surface boundary spans the entire top cell'
8790 h_L = (/5.,5./)
8800 call boundary_k_range(SURFACE, nk, h_L, 5., k_top, zeta_top, k_bot, zeta_bot)
881 near_boundary_unit_tests = near_boundary_unit_tests .or. &
8820 test_boundary_k_range(k_top, zeta_top, k_bot, zeta_bot, 1, 0., 1, 1., test_name, verbose)
883
8840 test_name = 'Surface boundary spans the entire column'
8850 h_L = (/5.,5./)
8860 call boundary_k_range(SURFACE, nk, h_L, 10., k_top, zeta_top, k_bot, zeta_bot)
887 near_boundary_unit_tests = near_boundary_unit_tests .or. &
8880 test_boundary_k_range(k_top, zeta_top, k_bot, zeta_bot, 1, 0., 2, 1., test_name, verbose)
889
8900 test_name = 'Bottom boundary spans the entire bottom cell'
8910 h_L = (/5.,5./)
8920 call boundary_k_range(BOTTOM, nk, h_L, 5., k_top, zeta_top, k_bot, zeta_bot)
893 near_boundary_unit_tests = near_boundary_unit_tests .or. &
8940 test_boundary_k_range(k_top, zeta_top, k_bot, zeta_bot, 2, 1., 2, 0., test_name, verbose)
895
8960 test_name = 'Bottom boundary spans the entire column'
8970 h_L = (/5.,5./)
8980 call boundary_k_range(BOTTOM, nk, h_L, 10., k_top, zeta_top, k_bot, zeta_bot)
899 near_boundary_unit_tests = near_boundary_unit_tests .or. &
9000 test_boundary_k_range(k_top, zeta_top, k_bot, zeta_bot, 1, 1., 2, 0., test_name, verbose)
901
9020 test_name = 'Surface boundary intersects second layer'
9030 h_L = (/10.,10./)
9040 call boundary_k_range(SURFACE, nk, h_L, 17.5, k_top, zeta_top, k_bot, zeta_bot)
905 near_boundary_unit_tests = near_boundary_unit_tests .or. &
9060 test_boundary_k_range(k_top, zeta_top, k_bot, zeta_bot, 1, 0., 2, 0.75, test_name, verbose)
907
9080 test_name = 'Surface boundary intersects first layer'
9090 h_L = (/10.,10./)
9100 call boundary_k_range(SURFACE, nk, h_L, 2.5, k_top, zeta_top, k_bot, zeta_bot)
911 near_boundary_unit_tests = near_boundary_unit_tests .or. &
9120 test_boundary_k_range(k_top, zeta_top, k_bot, zeta_bot, 1, 0., 1, 0.25, test_name, verbose)
913
9140 test_name = 'Surface boundary is deeper than column thickness'
9150 h_L = (/10.,10./)
9160 call boundary_k_range(SURFACE, nk, h_L, 21.0, k_top, zeta_top, k_bot, zeta_bot)
917 near_boundary_unit_tests = near_boundary_unit_tests .or. &
9180 test_boundary_k_range(k_top, zeta_top, k_bot, zeta_bot, 1, 0., 2, 1., test_name, verbose)
919
9200 test_name = 'Bottom boundary intersects first layer'
9210 h_L = (/10.,10./)
9220 call boundary_k_range(BOTTOM, nk, h_L, 17.5, k_top, zeta_top, k_bot, zeta_bot)
923 near_boundary_unit_tests = near_boundary_unit_tests .or. &
9240 test_boundary_k_range(k_top, zeta_top, k_bot, zeta_bot, 1, 0.75, 2, 0., test_name, verbose)
925
9260 test_name = 'Bottom boundary intersects second layer'
9270 h_L = (/10.,10./)
9280 call boundary_k_range(BOTTOM, nk, h_L, 2.5, k_top, zeta_top, k_bot, zeta_bot)
929 near_boundary_unit_tests = near_boundary_unit_tests .or. &
9300 test_boundary_k_range(k_top, zeta_top, k_bot, zeta_bot, 2, 0.25, 2, 0., test_name, verbose)
931
9320 if (.not. near_boundary_unit_tests) write(stdout,*) 'Passed boundary_k_range'
933
934 ! unit tests for sorting array and finding unique values
9350 test_name = 'Sorting array'
9360 eta1 = (/1., 0., 0.1/)
9370 call sort(eta1, nk+1)
938 near_boundary_unit_tests = near_boundary_unit_tests .or. &
9390 test_layer_fluxes( verbose, nk+1, test_name, eta1, (/0., 0.1, 1./) )
940
9410 test_name = 'Unique values'
9420 call unique((/0., 1., 1., 2./), nk+2, h1)
943 near_boundary_unit_tests = near_boundary_unit_tests .or. &
9440 test_layer_fluxes( verbose, nk+1, test_name, h1, (/0., 1., 2./) )
9450 deallocate(h1)
946
9470 test_name = 'Unique values with maximum depth'
9480 call unique((/0., 1., 1., 2., 3./), nk+3, h1, 2.)
949 near_boundary_unit_tests = near_boundary_unit_tests .or. &
9500 test_layer_fluxes( verbose, nk+1, test_name, h1, (/0., 1., 2./) )
9510 deallocate(h1)
952
9530 if (.not. near_boundary_unit_tests) write(stdout,*) 'Passed sort and unique'
954
955 ! unit tests for merge_interfaces
9560 test_name = 'h_L = h_R and BLD_L = BLD_R'
9570 call merge_interfaces(nk, (/1., 2./), (/1., 2./), 1.5, 1.5, CS%H_subroundoff, h1)
958 near_boundary_unit_tests = near_boundary_unit_tests .or. &
9590 test_layer_fluxes( verbose, nk, test_name, h1, (/1., 0.5/) )
9600 deallocate(h1)
961
9620 test_name = 'h_L = h_R and BLD_L /= BLD_R'
9630 call merge_interfaces(nk, (/1., 2./), (/1., 2./), 0.5, 1.5, CS%H_subroundoff, h1)
964 near_boundary_unit_tests = near_boundary_unit_tests .or. &
9650 test_layer_fluxes( verbose, nk+1, test_name, h1, (/0.5, 0.5, 0.5/) )
9660 deallocate(h1)
967
9680 test_name = 'h_L /= h_R and BLD_L = BLD_R'
9690 call merge_interfaces(nk, (/1., 3./), (/2., 2./), 1.5, 1.5, CS%H_subroundoff, h1)
970 near_boundary_unit_tests = near_boundary_unit_tests .or. &
9710 test_layer_fluxes( verbose, nk, test_name, h1, (/1., 0.5/) )
9720 deallocate(h1)
973
9740 test_name = 'h_L /= h_R and BLD_L /= BLD_R'
9750 call merge_interfaces(nk, (/1., 3./), (/2., 2./), 0.5, 1.5, CS%H_subroundoff, h1)
976 near_boundary_unit_tests = near_boundary_unit_tests .or. &
9770 test_layer_fluxes( verbose, nk+1, test_name, h1, (/0.5, 0.5, 0.5/) )
9780 deallocate(h1)
979
9800 test_name = 'Left deeper than right, h_L /= h_R and BLD_L /= BLD_R'
9810 call merge_interfaces(nk, (/2., 3./), (/2., 2./), 1.0, 2.0, CS%H_subroundoff, h1)
982 near_boundary_unit_tests = near_boundary_unit_tests .or. &
9830 test_layer_fluxes( verbose, nk, test_name, h1, (/1., 1./) )
9840 deallocate(h1)
985
9860 test_name = 'Left has zero thickness, h_L /= h_R and BLD_L = BLD_R'
9870 call merge_interfaces(nk, (/4., 0./), (/2., 2./), 2.0, 2.0, CS%H_subroundoff, h1)
988 near_boundary_unit_tests = near_boundary_unit_tests .or. &
9890 test_layer_fluxes( verbose, nk-1, test_name, h1, (/2./) )
9900 deallocate(h1)
991
9920 test_name = 'Left has zero thickness, h_L /= h_R and BLD_L /= BLD_R'
9930 call merge_interfaces(nk, (/4., 0./), (/2., 2./), 1.0, 2.0, CS%H_subroundoff, h1)
994 near_boundary_unit_tests = near_boundary_unit_tests .or. &
9950 test_layer_fluxes( verbose, nk, test_name, h1, (/1., 1./) )
9960 deallocate(h1)
997
9980 test_name = 'Right has zero thickness, h_L /= h_R and BLD_L = BLD_R'
9990 call merge_interfaces(nk, (/2., 2./), (/0., 4./), 2.0, 2.0, CS%H_subroundoff, h1)
1000 near_boundary_unit_tests = near_boundary_unit_tests .or. &
10010 test_layer_fluxes( verbose, nk-1, test_name, h1, (/2./) )
10020 deallocate(h1)
1003
10040 test_name = 'Right has zero thickness, h_L /= h_R and BLD_L /= BLD_R'
10050 call merge_interfaces(nk, (/2., 2./), (/0., 4./), 1.0, 2.0, CS%H_subroundoff, h1)
1006 near_boundary_unit_tests = near_boundary_unit_tests .or. &
10070 test_layer_fluxes( verbose, nk, test_name, h1, (/1., 1./) )
10080 deallocate(h1)
1009
10100 test_name = 'Right deeper than left, h_L /= h_R and BLD_L = BLD_R'
10110 call merge_interfaces(nk+1, (/2., 2., 0./), (/2., 2., 1./), 4., 4., CS%H_subroundoff, h1)
1012 near_boundary_unit_tests = near_boundary_unit_tests .or. &
10130 test_layer_fluxes( verbose, nk, test_name, h1, (/2., 2./) )
10140 deallocate(h1)
1015
10160 test_name = 'Right and left small values at bottom, h_L /= h_R and BLD_L = BLD_R'
10170 call merge_interfaces(nk+2, (/2., 2., 1., 1./), (/1., 1., .5, .5/), 3., 3., CS%H_subroundoff, h1)
1018 near_boundary_unit_tests = near_boundary_unit_tests .or. &
10190 test_layer_fluxes( verbose, nk+2, test_name, h1, (/1., 1., .5, .5/) )
10200 deallocate(h1)
1021
10220 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
10250 test_name = 'Equal hbl and same layer thicknesses (gradient from right to left)'
10260 hbl_L = 2. ; hbl_R = 2.
10270 h_L = (/2.,2./) ; h_R = (/2.,2./)
10280 phi_L = (/0.,0./) ; phi_R = (/1.,1./)
10290 khtr_u = (/1.,1.,1./)
10300 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, &
10320 khtr_u, F_layer, 1., 1., CS%hbd_u_kmax(1,1), CS%hbd_grd_u(1,1,:), CS)
1033 near_boundary_unit_tests = near_boundary_unit_tests .or. &
10340 test_layer_fluxes( verbose, nk, test_name, F_layer, (/-2.0,0.0/) )
1035
10360 test_name = 'Equal hbl and same layer thicknesses (gradient from left to right)'
10370 hbl_L = 2. ; hbl_R = 2.
10380 h_L = (/2.,2./) ; h_R = (/2.,2./)
10390 phi_L = (/2.,1./) ; phi_R = (/1.,1./)
10400 khtr_u = (/0.5,0.5,0.5/)
10410 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, &
10430 khtr_u, F_layer, 1., 1., CS%hbd_u_kmax(1,1), CS%hbd_grd_u(1,1,:), CS)
1044 near_boundary_unit_tests = near_boundary_unit_tests .or. &
10450 test_layer_fluxes( verbose, nk, test_name, F_layer, (/1.0,0.0/) )
1046
10470 test_name = 'hbl < column thickness, hbl same, linear profile right, khtr=2'
10480 hbl_L = 2 ; hbl_R = 2
10490 h_L = (/1.,2./) ; h_R = (/1.,2./)
10500 phi_L = (/0.,0./) ; phi_R = (/0.5,2./)
10510 khtr_u = (/2.,2.,2./)
10520 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, &
10540 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.
1056 near_boundary_unit_tests = near_boundary_unit_tests .or. &
10570 test_layer_fluxes( verbose, nk, test_name, F_layer, (/-1.0,-4.0/) )
1058
10590 test_name = 'Different hbl and different column thicknesses (zero gradient)'
10600 hbl_L = 12 ; hbl_R = 20
10610 h_L = (/6.,6./) ; h_R = (/10.,10./)
10620 phi_L = (/1.,1./) ; phi_R = (/1.,1./)
10630 khtr_u = (/1.,1.,1./)
10640 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, &
10660 khtr_u, F_layer, 1., 1., CS%hbd_u_kmax(1,1), CS%hbd_grd_u(1,1,:), CS)
1067 near_boundary_unit_tests = near_boundary_unit_tests .or. &
10680 test_layer_fluxes( verbose, nk, test_name, F_layer, (/0.,0./) )
1069
10700 test_name = 'Different hbl and different column thicknesses (gradient from left to right)'
1071
10720 hbl_L = 15 ; hbl_R = 10.
10730 h_L = (/10.,5./) ; h_R = (/10.,0./)
10740 phi_L = (/1.,1./) ; phi_R = (/0.,0./)
10750 khtr_u = (/1.,1.,1./)
10760 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, &
10780 khtr_u, F_layer, 1., 1., CS%hbd_u_kmax(1,1), CS%hbd_grd_u(1,1,:), CS)
1079 near_boundary_unit_tests = near_boundary_unit_tests .or. &
10800 test_layer_fluxes( verbose, nk, test_name, F_layer, (/10.,0.0/) )
1081
10820 if (.not. near_boundary_unit_tests) write(stdout,*) 'Passed fluxes_layer_method'
1083
10840end 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
10880logical 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
10970 test_layer_fluxes = .false.
10980 do k=1,nk
10990 if ( F_calc(k) /= F_ans(k) ) then
11000 test_layer_fluxes = .true.
11010 write(stdout,*) "MOM_hor_bnd_diffusion, UNIT TEST FAILED: ", test_name
11020 write(stdout,10) k, F_calc(k), F_ans(k)
11030 elseif (verbose) then
11040 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)
11090end function test_layer_fluxes
1110
1111!> Return true if output of unit tests for boundary_k_range does not match answers
11120logical 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
11250 test_boundary_k_range = k_top /= k_top_ans
11260 test_boundary_k_range = test_boundary_k_range .or. (zeta_top /= zeta_top_ans)
11270 test_boundary_k_range = test_boundary_k_range .or. (k_bot /= k_bot_ans)
11280 test_boundary_k_range = test_boundary_k_range .or. (zeta_bot /= zeta_bot_ans)
1129
11300 if (test_boundary_k_range) write(stdout,*) "UNIT TEST FAILED: ", test_name
11310 if (test_boundary_k_range .or. verbose) then
11320 write(stdout,20) "k_top", k_top, "k_top_ans", k_top_ans
11330 write(stdout,20) "k_bot", k_bot, "k_bot_ans", k_bot_ans
11340 write(stdout,30) "zeta_top", zeta_top, "zeta_top_ans", zeta_top_ans
11350 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
11420end function test_boundary_k_range
1143
1144!> Same as hbd_grid, but only used in the unit tests.
11450subroutine 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
11540 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
11580 CS%hbd_grd_u(1,1,:) = 0.0
11590 CS%hbd_u_kmax(1,1) = 0
1160
11610 call merge_interfaces(2, h_L, h_R, hbl_L, hbl_R, CS%H_subroundoff, dz_top)
11620 nk = SIZE(dz_top)
11630 if (nk > CS%hbd_nk) then
11640 write(*,*)'nk, CS%hbd_nk', nk, CS%hbd_nk
11650 call MOM_error(FATAL,"Houston, we've had a problem in hbd_grid_test, (nk cannot be > CS%hbd_nk)")
1166 endif
1167
11680 CS%hbd_u_kmax(1,1) = nk
1169
1170 ! set the HBD grid to dz_top
11710 do k=1,nk
11720 CS%hbd_grd_u(1,1,k) = dz_top(k)
1173 enddo
11740 deallocate(dz_top)
1175
11760end subroutine hbd_grid_test
1177
1178!> Deallocates hor_bnd_diffusion control structure
11791subroutine hor_bnd_diffusion_end(CS)
1180 type(hbd_CS), pointer :: CS !< Horizontal boundary diffusion control structure
1181
11821 if (associated(CS)) deallocate(CS)
1183
11841end 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!!
12510end module MOM_hor_bnd_diffusion