src/parameterizations/vertical/MOM_set_diffusivity.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		!> Calculate vertical diffusivity from all mixing processes
6		module MOM_set_diffusivity
7
8		use MOM_bkgnd_mixing, only : calculate_bkgnd_mixing, bkgnd_mixing_init, bkgnd_mixing_cs
9		use MOM_bkgnd_mixing, only : bkgnd_mixing_end
10		use MOM_cpu_clock, only : cpu_clock_id, cpu_clock_begin, cpu_clock_end
11		use MOM_cpu_clock, only : CLOCK_MODULE_DRIVER, CLOCK_MODULE, CLOCK_ROUTINE
12		use MOM_CVMix_ddiff, only : CVMix_ddiff_init, CVMix_ddiff_end, CVMix_ddiff_cs
13		use MOM_CVMix_ddiff, only : compute_ddiff_coeffs
14		use MOM_CVMix_shear, only : calculate_CVMix_shear, CVMix_shear_init, CVMix_shear_cs
15		use MOM_CVMix_shear, only : CVMix_shear_end
16		use MOM_diag_mediator, only : diag_ctrl, time_type
17		use MOM_diag_mediator, only : post_data, register_diag_field
18		use MOM_diagnose_kdwork, only : vbf_CS
19		use MOM_debugging, only : hchksum, uvchksum, Bchksum, hchksum_pair
20		use MOM_EOS, only : calculate_density, calculate_density_derivs, EOS_domain
21		use MOM_error_handler, only : MOM_error, is_root_pe, FATAL, WARNING, NOTE
22		use MOM_error_handler, only : callTree_showQuery
23		use MOM_error_handler, only : callTree_enter, callTree_leave, callTree_waypoint
24		use MOM_file_parser, only : get_param, log_param, log_version, param_file_type
25		use MOM_forcing_type, only : forcing, optics_type
26		use MOM_full_convection, only : full_convection
27		use MOM_grid, only : ocean_grid_type
28		use MOM_interface_heights, only : thickness_to_dz, find_rho_bottom
29		use MOM_internal_tides, only : int_tide_CS, get_lowmode_loss, get_lowmode_diffusivity
30		use MOM_intrinsic_functions, only : invcosh
31		use MOM_io, only : slasher, MOM_read_data
32		use MOM_isopycnal_slopes, only : vert_fill_TS
33		use MOM_kappa_shear, only : calculate_kappa_shear, kappa_shear_init, Kappa_shear_CS
34		use MOM_kappa_shear, only : calc_kappa_shear_vertex, kappa_shear_at_vertex
35		use MOM_open_boundary, only : ocean_OBC_type, OBC_segment_type, OBC_NONE
36		use MOM_open_boundary, only : OBC_DIRECTION_E, OBC_DIRECTION_W, OBC_DIRECTION_N, OBC_DIRECTION_S
37		use MOM_string_functions, only : uppercase
38		use MOM_tidal_mixing, only : tidal_mixing_CS, calculate_tidal_mixing, tidal_mixing_h_amp
39		use MOM_tidal_mixing, only : setup_tidal_diagnostics, post_tidal_diagnostics
40		use MOM_tidal_mixing, only : tidal_mixing_init, tidal_mixing_end
41		use MOM_unit_scaling, only : unit_scale_type
42		use MOM_variables, only : thermo_var_ptrs, vertvisc_type, p3d
43		use MOM_verticalGrid, only : verticalGrid_type
44		use user_change_diffusivity, only : user_change_diff, user_change_diff_init
45		use user_change_diffusivity, only : user_change_diff_end, user_change_diff_CS
46
47		implicit none ; private
48
49		#include <MOM_memory.h>
50
51		public set_diffusivity
52		public set_BBL_TKE
53		public set_diffusivity_init
54		public set_diffusivity_end
55
56		! A note on unit descriptions in comments: MOM6 uses units that can be rescaled for dimensional
57		! consistency testing. These are noted in comments with units like Z, H, L, and T, along with
58		! their mks counterparts with notation like "a velocity [Z T-1 ~> m s-1]". If the units
59		! vary with the Boussinesq approximation, the Boussinesq variant is given first.
60
61		!> This control structure contains parameters for MOM_set_diffusivity.
62		type, public :: set_diffusivity_CS ; private
63		logical :: initialized = .false. !< True if this control structure has been initialized.
64		logical :: debug !< If true, write verbose checksums for debugging.
65
66		logical :: bulkmixedlayer !< If true, a refined bulk mixed layer is used with
67		!! GV%nk_rho_varies variable density mixed & buffer layers.
68		real :: FluxRi_max !< The flux Richardson number where the stratification is
69		!! large enough that N2 > omega2 [nondim]. The full expression
70		!! for the Flux Richardson number is usually
71		!! FLUX_RI_MAX*N2/(N2+OMEGA2). The default is 0.2.
72		logical :: bottomdraglaw !< If true, the bottom stress is calculated with a
73		!! drag law c_drag\|u\|u.
74		logical :: BBL_mixing_as_max !< If true, take the maximum of the diffusivity
75		!! from the BBL mixing and the other diffusivities.
76		!! Otherwise, diffusivities from the BBL_mixing is added.
77		logical :: use_LOTW_BBL_diffusivity !< If true, use simpler/less precise, BBL diffusivity.
78		logical :: LOTW_BBL_use_omega !< If true, use simpler/less precise, BBL diffusivity.
79		real :: Von_Karm !< The von Karman constant as used in the BBL diffusivity calculation
80		!! [nondim]. See (http://en.wikipedia.org/wiki/Von_Karman_constant)
81		real :: BBL_effic !< Efficiency with which the energy extracted
82		!! by bottom drag drives BBL diffusion in the original BBL scheme, times
83		!! conversion factors between the natural units of mean kinetic energy
84		!! and those those used for TKE [Z2 L-2 ~> nondim].
85		real :: ePBL_BBL_effic !< efficiency with which the energy extracted
86		!! by bottom drag drives BBL diffusion in the ePBL BBL scheme [nondim]
87		logical :: ePBL_BBL_mstar !< logical if the bottom boundary layer uses an mstar x ustar^3 formulation
88		!! needed here to know whether or not to populate the bottom ustar
89		real :: cdrag !< quadratic drag coefficient [nondim]
90		real :: dz_BBL_avg_min !< A minimal distance over which to average to determine the average
91		!! bottom boundary layer density [Z ~> m]
92		real :: IMax_decay !< Inverse of a maximum decay scale for
93		!! bottom-drag driven turbulence [H-1 ~> m-1 or m2 kg-1].
94		real :: Kv !< The interior vertical viscosity [H Z T-1 ~> m2 s-1 or Pa s]
95		real :: Kd !< interior diapycnal diffusivity [H Z T-1 ~> m2 s-1 or kg m-1 s-1]
96		real :: Kd_min !< minimum diapycnal diffusivity [H Z T-1 ~> m2 s-1 or kg m-1 s-1]
97		real :: Kd_max !< maximum increment for diapycnal diffusivity [H Z T-1 ~> m2 s-1 or kg m-1 s-1]
98		!! Set to a negative value to have no limit.
99		real :: Kd_add !< uniform diffusivity added everywhere without
100		!! filtering or scaling [H Z T-1 ~> m2 s-1 or kg m-1 s-1]
101		real :: Kd_smooth !< Vertical diffusivity used to interpolate more
102		!! sensible values of T & S into thin layers [H Z T-1 ~> m2 s-1 or kg m-1 s-1]
103		type(diag_ctrl), pointer :: diag => NULL() !< structure to regulate diagnostic output timing
104
105		logical :: limit_dissipation !< If enabled, dissipation is limited to be larger
106		!! than the following:
107		real :: dissip_min !< Minimum dissipation [R Z2 T-3 ~> W m-3]
108		real :: dissip_N0 !< Coefficient a in minimum dissipation = a+b*N [R Z2 T-3 ~> W m-3]
109		real :: dissip_N1 !< Coefficient b in minimum dissipation = a+b*N [R Z2 T-2 ~> J m-3]
110		real :: dissip_N2 !< Coefficient c in minimum dissipation = c*N2 [R Z2 T-1 ~> J s m-3]
111		real :: dissip_Kd_min !< Minimum Kd [H Z T-1 ~> m2 s-1 or kg m-1 s-1], with dissipation Rho0Kd_minN^2
112
113		real :: omega !< Earth's rotation frequency [T-1 ~> s-1]
114		logical :: ML_radiation !< allow a fraction of TKE available from wind work
115		!! to penetrate below mixed layer base with a vertical
116		!! decay scale determined by the minimum of
117		!! (1) The depth of the mixed layer, or
118		!! (2) An Ekman length scale.
119		!! Energy available to drive mixing below the mixed layer is
120		!! given by E = ML_RAD_COEFFMSTARUSTAR**3. Optionally, if
121		!! ML_rad_TKE_decay is true, this is further reduced by a factor
122		!! of exp(-h_ML*Idecay_len_TkE), where Idecay_len_TKE is
123		!! calculated the same way as in the mixed layer code.
124		!! The diapycnal diffusivity is KD(k) = E/(N2(k)+OMEGA2),
125		!! where N2 is the squared buoyancy frequency [T-2 ~> s-2] and OMEGA2
126		!! is the rotation rate of the earth squared.
127		real :: ML_rad_kd_max !< Maximum diapycnal diffusivity due to turbulence radiated from
128		!! the base of the mixed layer [H Z T-1 ~> m2 s-1 or kg m-1 s-1].
129		real :: ML_rad_efold_coeff !< Coefficient to scale penetration depth [nondim]
130		real :: ML_rad_coeff !< Coefficient which scales MSTAR*USTAR^3 to obtain energy
131		!! available for mixing below mixed layer base [nondim]
132		logical :: ML_rad_bug !< If true use code with a bug that reduces the energy available
133		!! in the transition layer by a factor of the inverse of the energy
134		!! deposition lenthscale (in m).
135		logical :: ML_rad_TKE_decay !< If true, apply same exponential decay
136		!! to ML_rad as applied to the other surface
137		!! sources of TKE in the mixed layer code.
138		real :: ustar_min !< A minimum value of ustar to avoid numerical
139		!! problems [Z T-1 ~> m s-1]. If the value is small enough,
140		!! this parameter should not affect the solution.
141		real :: TKE_decay !< ratio of natural Ekman depth to TKE decay scale [nondim]
142		real :: mstar !< ratio of friction velocity cubed to
143		!! TKE input to the mixed layer [nondim]
144		logical :: ML_use_omega !< If true, use absolute rotation rate instead
145		!! of the vertical component of rotation when
146		!! setting the decay scale for mixed layer turbulence.
147		real :: ML_omega_frac !< When setting the decay scale for turbulence, use
148		!! this fraction [nondim] of the absolute rotation rate blended
149		!! with the local value of f, as f^2 ~= (1-of)f^2 + of4*omega^2.
150		logical :: user_change_diff !< If true, call user-defined code to change diffusivity.
151		logical :: useKappaShear !< If true, use the kappa_shear module to find the
152		!! shear-driven diapycnal diffusivity.
153		logical :: Vertex_Shear !< If true, do the calculations of the shear-driven mixing
154		!! at the cell vertices (i.e., the vorticity points).
155		logical :: use_CVMix_shear !< If true, use one of the CVMix modules to find
156		!! shear-driven diapycnal diffusivity.
157		logical :: double_diffusion !< If true, enable double-diffusive mixing using an old method.
158		logical :: use_CVMix_ddiff !< If true, enable double-diffusive mixing via CVMix.
159		logical :: use_tidal_mixing !< If true, activate tidal mixing diffusivity.
160		logical :: use_int_tides !< If true, use internal tides ray tracing
161		logical :: simple_TKE_to_Kd !< If true, uses a simple estimate of Kd/TKE that
162		!! does not rely on a layer-formulation.
163		real :: Max_Rrho_salt_fingers !< max density ratio for salt fingering [nondim]
164		real :: Max_salt_diff_salt_fingers !< max salt diffusivity for salt fingers [H Z T-1 ~> m2 s-1 or kg m-1 s-1]
165		real :: Kv_molecular !< Molecular viscosity for double diffusive convection [H Z T-1 ~> m2 s-1 or Pa s]
166
167		integer :: answer_date !< The vintage of the order of arithmetic and expressions in this module's
168		!! calculations. Values below 20190101 recover the answers from the
169		!! end of 2018, while higher values use updated and more robust forms
170		!! of the same expressions. Values above 20240630 use more accurate
171		!! expressions for cases where USE_LOTW_BBL_DIFFUSIVITY is true. Values
172		!! above 20250301 use less confusing expressions to set the bottom-drag
173		!! generated diffusivity when USE_LOTW_BBL_DIFFUSIVITY is false.
174		integer :: LOTW_BBL_answer_date !< The vintage of the order of arithmetic and expressions
175		!! in the LOTW_BBL calculations. Values below 20240630 recover the
176		!! original answers, while higher values use more accurate expressions.
177		!! This only applies when USE_LOTW_BBL_DIFFUSIVITY is true.
178		integer :: drag_diff_answer_date !< The vintage of the order of arithmetic in the drag diffusivity
179		!! calculations. Values above 20250301 use less confusing expressions
180		!! to set the bottom-drag generated diffusivity when
181		!! USE_LOTW_BBL_DIFFUSIVITY is false.
182
183		character(len=200) :: inputdir !< The directory in which input files are found
184		type(user_change_diff_CS), pointer :: user_change_diff_CSp => NULL() !< Control structure for a child module
185		type(Kappa_shear_CS), pointer :: kappaShear_CSp => NULL() !< Control structure for a child module
186		type(CVMix_shear_cs), pointer :: CVMix_shear_csp => NULL() !< Control structure for a child module
187		type(CVMix_ddiff_cs), pointer :: CVMix_ddiff_csp => NULL() !< Control structure for a child module
188		type(bkgnd_mixing_cs), pointer :: bkgnd_mixing_csp => NULL() !< Control structure for a child module
189		type(int_tide_CS), pointer :: int_tide_CSp => NULL() !< Control structure for a child module
190		type(tidal_mixing_cs) :: tidal_mixing !< Control structure for a child module
191
192		!>@{ Diagnostic IDs
193		integer :: id_maxTKE = -1, id_TKE_to_Kd = -1, id_Kd_user = -1
194		integer :: id_Kd_layer = -1, id_Kd_BBL = -1, id_N2 = -1
195		integer :: id_Kd_Work = -1, id_KT_extra = -1, id_KS_extra = -1, id_R_rho = -1
196		integer :: id_Kd_bkgnd = -1, id_Kv_bkgnd = -1, id_Kd_leak = -1
197		integer :: id_Kd_quad = -1, id_Kd_itidal = -1, id_Kd_Froude = -1, id_Kd_slope = -1
198		integer :: id_prof_leak = -1, id_prof_quad = -1, id_prof_itidal= -1
199		integer :: id_prof_Froude= -1, id_prof_slope = -1, id_bbl_thick = -1, id_kbbl = -1
200		integer :: id_Kd_Work_added = -1
201		!>@}
202
203		end type set_diffusivity_CS
204
205		!> This structure has memory for used in calculating diagnostics of diffusivity
206		type diffusivity_diags
207		real, pointer, dimension(:,:,:) :: &
208		N2_3d => NULL(), & !< squared buoyancy frequency at interfaces [T-2 ~> s-2]
209		Kd_user => NULL(), & !< user-added diffusivity at interfaces [H Z T-1 ~> m2 s-1 or kg m-1 s-1]
210		Kd_BBL => NULL(), & !< BBL diffusivity at interfaces [H Z T-1 ~> m2 s-1 or kg m-1 s-1]
211		Kd_Work => NULL(), & !< layer integrated work by diapycnal mixing [R Z3 T-3 ~> W m-2]
212		Kd_Work_added => NULL(), & !< layer integrated work by added mixing [R Z3 T-3 ~> W m-2]
213		maxTKE => NULL(), & !< energy required to entrain to h_max [H Z2 T-3 ~> m3 s-3 or W m-2]
214		Kd_bkgnd => NULL(), & !< Background diffusivity at interfaces [H Z T-1 ~> m2 s-1 or kg m-1 s-1]
215		Kv_bkgnd => NULL(), & !< Viscosity from background diffusivity at interfaces [H Z T-1 ~> m2 s-1 or Pa s]
216		KT_extra => NULL(), & !< Double diffusion diffusivity for temperature [H Z T-1 ~> m2 s-1 or kg m-1 s-1]
217		KS_extra => NULL(), & !< Double diffusion diffusivity for salinity [H Z T-1 ~> m2 s-1 or kg m-1 s-1]
218		drho_rat => NULL(), & !< The density difference ratio used in double diffusion [nondim].
219		Kd_leak => NULL(), & !< internal tides leakage diffusivity [H Z T-1 ~> m2 s-1 or kg m-1 s-1]
220		Kd_quad => NULL(), & !< internal tides bottom drag diffusivity [H Z T-1 ~> m2 s-1 or kg m-1 s-1]
221		Kd_itidal => NULL(), & !< internal tides wave drag diffusivity [H Z T-1 ~> m2 s-1 or kg m-1 s-1]
222		Kd_Froude => NULL(), & !< internal tides high Froude diffusivity [H Z T-1 ~> m2 s-1 or kg m-1 s-1]
223		Kd_slope => NULL(), & !< internal tides critical slopes diffusivity [H Z T-1 ~> m2 s-1 or kg m-1 s-1]
224		prof_leak => NULL(), & !< vertical profile for leakage [H-1 ~> m-1 or m2 kg-1]
225		prof_quad => NULL(), & !< vertical profile for bottom drag [H-1 ~> m-1 or m2 kg-1]
226		prof_itidal => NULL(), & !< vertical profile for wave drag [H-1 ~> m-1 or m2 kg-1]
227		prof_Froude => NULL(), & !< vertical profile for Froude drag [H-1 ~> m-1 or m2 kg-1]
228		prof_slope => NULL() !< vertical profile for critical slopes [H-1 ~> m-1 or m2 kg-1]
229		real, pointer, dimension(:,:) :: bbl_thick => NULL(), & !< bottom boundary layer thickness [H ~> m or kg m-2]
230		kbbl => NULL() !< top of bottom boundary layer
231
232		real, pointer, dimension(:,:,:) :: TKE_to_Kd => NULL()
233		!< conversion rate (~1.0 / (G_Earth + dRho_lay)) between TKE
234		!! dissipated within a layer and Kd in that layer [T2 Z-1 ~> s2 m-1]
235
236		end type diffusivity_diags
237
238		!>@{ CPU time clocks
239		integer :: id_clock_kappaShear, id_clock_CVMix_ddiff
240		!>@}
241
242		contains
243
244	24	subroutine set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, optics, visc, dt, Kd_int, &
245	0	G, GV, US, CS, VBF, Kd_lay, Kd_extra_T, Kd_extra_S)
246		type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure.
247		type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure.
248		type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
249		real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &
250		intent(in) :: u !< The zonal velocity [L T-1 ~> m s-1].
251		real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &
252		intent(in) :: v !< The meridional velocity [L T-1 ~> m s-1].
253		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
254		intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2].
255		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
256		intent(in) :: u_h !< Zonal velocity interpolated to h points [L T-1 ~> m s-1].
257		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
258		intent(in) :: v_h !< Meridional velocity interpolated to h points [L T-1 ~> m s-1].
259		type(thermo_var_ptrs), intent(inout) :: tv !< Structure with pointers to thermodynamic
260		!! fields. Out is for tv%TempxPmE.
261		type(forcing), intent(in) :: fluxes !< A structure of thermodynamic surface fluxes
262		type(optics_type), pointer :: optics !< A structure describing the optical
263		!! properties of the ocean.
264		type(vertvisc_type), intent(inout) :: visc !< Structure containing vertical viscosities, bottom
265		!! boundary layer properties and related fields.
266		real, intent(in) :: dt !< Time increment [T ~> s].
267		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), &
268		intent(out) :: Kd_int !< Diapycnal diffusivity at each interface
269		!! [H Z T-1 ~> m2 s-1 or kg m-1 s-1].
270		type(set_diffusivity_CS), pointer :: CS !< Module control structure.
271		type(vbf_CS), pointer :: VBF !< A diagnostic control structure for vertical buoyancy fluxes
272		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
273		optional, intent(out) :: Kd_lay !< Diapycnal diffusivity of each layer
274		!! [H Z T-1 ~> m2 s-1 or kg m-1 s-1].
275
276		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), &
277		optional, intent(out) :: Kd_extra_T !< The extra diffusivity at interfaces of
278		!! temperature due to double diffusion relative
279		!! to the diffusivity of density
280		!! [H Z T-1 ~> m2 s-1 or kg m-1 s-1]
281		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), &
282		optional, intent(out) :: Kd_extra_S !< The extra diffusivity at interfaces of
283		!! salinity due to double diffusion relative
284		!! to the diffusivity of density
285		!! [H Z T-1 ~> m2 s-1 or kg m-1 s-1]
286
287		! local variables
288	24	real :: N2_bot(SZI_(G)) ! Bottom squared buoyancy frequency [T-2 ~> s-2]
289	24	real :: rho_bot(SZI_(G)) ! In situ near-bottom density [T-2 ~> s-2]
290	24	real :: h_bot(SZI_(G)) ! Bottom boundary layer thickness [H ~> m or kg m-2]
291	24	integer :: k_bot(SZI_(G)) ! Bottom boundary layer thickness top layer index
292
293		type(diffusivity_diags) :: dd ! structure with arrays of available diags
294
295
296		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: &
297	24	T_f, S_f ! Temperature and salinity [C ~> degC] and [S ~> ppt] with properties in massless layers
298		! filled vertically by diffusion or the properties after full convective adjustment.
299
300		real, dimension(SZI_(G),SZK_(GV)) :: &
301	24	N2_lay, & !< Squared buoyancy frequency associated with layers [T-2 ~> s-2]
302	24	Kd_lay_2d, & !< The layer diffusivities [H Z T-1 ~> m2 s-1 or kg m-1 s-1]
303	24	dz, & !< Height change across layers [Z ~> m]
304	24	maxTKE, & !< Energy required to entrain to h_max [H Z2 T-3 ~> m3 s-3 or W m-2]
305	24	prof_leak_2d, & !< vertical profile for leakage [Z-1 ~> m-1]
306	24	prof_quad_2d, & !< vertical profile for bottom drag [Z-1 ~> m-1]
307	24	prof_itidal_2d, & !< vertical profile for wave drag [Z-1 ~> m-1]
308	24	prof_Froude_2d, & !< vertical profile for Froude drag [Z-1 ~> m-1]
309	24	prof_slope_2d, & !< vertical profile for critical slopes [Z-1 ~> m-1]
310	24	TKE_to_Kd !< Conversion rate (~1.0 / (G_Earth + dRho_lay)) between
311		!< TKE dissipated within a layer and Kd in that layer [T2 Z-1 ~> s2 m-1]
312
313		real, dimension(SZI_(G),SZK_(GV)+1) :: &
314	24	N2_int, & !< squared buoyancy frequency associated at interfaces [T-2 ~> s-2]
315	24	Kd_int_2d, & !< The interface diffusivities [H Z T-1 ~> m2 s-1 or kg m-1 s-1]
316	24	Kv_bkgnd, & !< The background diffusion related interface viscosities [H Z T-1 ~> m2 s-1 or Pa s]
317	24	Kd_leak_2d, & !< internal tides leakage diffusivity [H Z T-1 ~> m2 s-1 or kg m-1 s-1]
318	24	Kd_quad_2d, & !< internal tides bottom drag diffusivity [H Z T-1 ~> m2 s-1 or kg m-1 s-1]
319	24	Kd_itidal_2d, & !< internal tides wave drag diffusivity [H Z T-1 ~> m2 s-1 or kg m-1 s-1]
320	24	Kd_Froude_2d, & !< internal tides high Froude diffusivity [H Z T-1 ~> m2 s-1 or kg m-1 s-1]
321	24	Kd_slope_2d, & !< internal tides critical slopes diffusivity [H Z T-1 ~> m2 s-1 or kg m-1 s-1]
322	24	dRho_int, & !< Locally referenced potential density difference across interfaces [R ~> kg m-3]
323	24	KT_extra, & !< Double diffusion diffusivity of temperature [H Z T-1 ~> m2 s-1 or kg m-1 s-1]
324	24	KS_extra !< Double diffusion diffusivity of salinity [H Z T-1 ~> m2 s-1 or kg m-1 s-1]
325
326		real :: dissip ! local variable for dissipation calculations [Z2 R T-3 ~> W m-3]
327		real :: Omega2 ! squared absolute rotation rate [T-2 ~> s-2]
328
329		logical :: use_EOS ! If true, compute density from T/S using equation of state.
330		logical :: TKE_to_Kd_used ! If true, TKE_to_Kd and maxTKE need to be calculated.
331	24	integer :: kb(SZI_(G)) ! The index of the lightest layer denser than the
332		! buffer layer, or -1 without a bulk mixed layer.
333		logical :: showCallTree ! If true, show the call tree.
334
335		integer :: i, j, k, is, ie, js, je, nz, isd, ied, jsd, jed
336
337		real :: kappa_dt_fill ! diffusivity times a timestep used to fill massless layers [H Z ~> m2 or kg m-1]
338
339	12	is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke
340	12	isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed
341	12	showCallTree = callTree_showQuery()
342	0	if (showCallTree) call callTree_enter("set_diffusivity(), MOM_set_diffusivity.F90")
343
344	12	if (.not.associated(CS)) call MOM_error(FATAL,"set_diffusivity: "//&
345	0	"Module must be initialized before it is used.")
346
347	12	if (.not.CS%initialized) call MOM_error(FATAL,"set_diffusivity: "//&
348	0	"Module must be initialized before it is used.")
349
350	12	if (CS%answer_date < 20190101) then
351		! These hard-coded dimensional parameters are being replaced.
352	0	kappa_dt_fill = 1.e-3GV%m2_s_to_HZ_T 7200.*US%s_to_T
353		else
354	12	kappa_dt_fill = CS%Kd_smooth * dt
355		endif
356	12	Omega2 = CS%omega * CS%omega
357
358	12	use_EOS = associated(tv%eqn_of_state)
359
360	12	if ((CS%use_CVMix_ddiff .or. CS%double_diffusion) .and. &
361		.not.(present(Kd_extra_T) .and. present(Kd_extra_S))) &
362		call MOM_error(FATAL, "set_diffusivity: both Kd_extra_T and Kd_extra_S must be present "//&
363	0	"when USE_CVMIX_DDIFF or DOUBLE_DIFFUSION are true.")
364
365		TKE_to_Kd_used = (CS%use_tidal_mixing .or. CS%ML_radiation .or. &
366		CS%use_int_tides .or. &
367	12	(CS%bottomdraglaw .and. .not.CS%use_LOTW_BBL_diffusivity))
368
369		! Set Kd_lay, Kd_int and Kv_slow to constant values, mostly to fill the halos.
370	12	if (present(Kd_lay)) Kd_lay(:,:,:) = CS%Kd
371	8000988	Kd_int(:,:,:) = CS%Kd
372	12	if (present(Kd_extra_T)) Kd_extra_T(:,:,:) = 0.0
373	12	if (present(Kd_extra_S)) Kd_extra_S(:,:,:) = 0.0
374	12	if (associated(visc%Kv_slow)) visc%Kv_slow(:,:,:) = CS%Kv
375
376		! Set up arrays for diagnostics.
377
378	12	if (CS%id_N2 > 0) allocate(dd%N2_3d(isd:ied,jsd:jed,nz+1), source=0.0)
379	12	if (CS%id_Kd_user > 0) allocate(dd%Kd_user(isd:ied,jsd:jed,nz+1), source=0.0)
380	12	if (CS%id_Kd_Work > 0) allocate(dd%Kd_Work(isd:ied,jsd:jed,nz), source=0.0)
381	12	if (CS%id_Kd_Work_added > 0) allocate(dd%Kd_Work_added(isd:ied,jsd:jed,nz), source=0.0)
382	12	if (CS%id_maxTKE > 0) allocate(dd%maxTKE(isd:ied,jsd:jed,nz), source=0.0)
383	12	if (CS%id_TKE_to_Kd > 0) allocate(dd%TKE_to_Kd(isd:ied,jsd:jed,nz), source=0.0)
384	12	if ((CS%double_diffusion) .and. (CS%id_KT_extra > 0)) &
385	0	allocate(dd%KT_extra(isd:ied,jsd:jed,nz+1), source=0.0)
386	12	if ((CS%double_diffusion) .and. (CS%id_KS_extra > 0)) &
387	0	allocate(dd%KS_extra(isd:ied,jsd:jed,nz+1), source=0.0)
388	12	if (CS%id_R_rho > 0) allocate(dd%drho_rat(isd:ied,jsd:jed,nz+1), source=0.0)
389	12	if (CS%id_Kd_BBL > 0 .or. associated(VBF%Kd_BBL)) &
390	0	allocate(dd%Kd_BBL(isd:ied,jsd:jed,nz+1), source=0.0)
391
392	12	if (CS%id_Kd_bkgnd > 0) allocate(dd%Kd_bkgnd(isd:ied,jsd:jed,nz+1), source=0.)
393	12	if (CS%id_Kv_bkgnd > 0) allocate(dd%Kv_bkgnd(isd:ied,jsd:jed,nz+1), source=0.)
394
395	12	if (CS%id_kbbl > 0) allocate(dd%kbbl(isd:ied,jsd:jed), source=0.)
396	12	if (CS%id_bbl_thick > 0) allocate(dd%bbl_thick(isd:ied,jsd:jed), source=0.)
397	12	if (CS%id_Kd_leak > 0) allocate(dd%Kd_leak(isd:ied,jsd:jed,nz+1), source=0.)
398	12	if (CS%id_Kd_quad > 0) allocate(dd%Kd_quad(isd:ied,jsd:jed,nz+1), source=0.)
399	12	if (CS%id_Kd_itidal > 0) allocate(dd%Kd_itidal(isd:ied,jsd:jed,nz+1), source=0.)
400	12	if (CS%id_Kd_Froude > 0) allocate(dd%Kd_Froude(isd:ied,jsd:jed,nz+1), source=0.)
401	12	if (CS%id_Kd_slope > 0) allocate(dd%Kd_slope(isd:ied,jsd:jed,nz+1), source=0.)
402
403	12	if (CS%id_prof_leak > 0) allocate(dd%prof_leak(isd:ied,jsd:jed,nz), source=0.)
404	12	if (CS%id_prof_quad > 0) allocate(dd%prof_quad(isd:ied,jsd:jed,nz), source=0.)
405	12	if (CS%id_prof_itidal > 0) allocate(dd%prof_itidal(isd:ied,jsd:jed,nz), source=0.)
406	12	if (CS%id_prof_Froude > 0) allocate(dd%prof_Froude(isd:ied,jsd:jed,nz), source=0.)
407	12	if (CS%id_prof_slope > 0) allocate(dd%prof_slope(isd:ied,jsd:jed,nz), source=0.)
408
409		! set up arrays for tidal mixing diagnostics
410	12	if (CS%use_tidal_mixing) &
411	0	call setup_tidal_diagnostics(G, GV, CS%tidal_mixing)
412
413	12	if (CS%useKappaShear) then
414	12	if (CS%debug) then
415	0	call hchksum_pair("before calc_KS [uv]_h", u_h, v_h, G%HI, unscale=US%L_T_to_m_s)
416		endif
417	12	call cpu_clock_begin(id_clock_kappaShear)
418	12	if (CS%Vertex_shear) then
419		call full_convection(G, GV, US, h, tv, T_f, S_f, fluxes%p_surf, &
420	12	kappa_dt_fill, halo=1)
421
422		call calc_kappa_shear_vertex(u, v, h, T_f, S_f, tv, fluxes%p_surf, visc%Kd_shear, &
423	12	visc%TKE_turb, visc%Kv_shear_Bu, dt, G, GV, US, CS%kappaShear_CSp)
424	8000988	if (associated(visc%Kv_shear)) visc%Kv_shear(:,:,:) = 0.0 ! needed for other parameterizations
425	12	if (CS%debug) then
426	0	call hchksum(visc%Kd_shear, "after calc_KS_vert visc%Kd_shear", G%HI, unscale=GV%HZ_T_to_m2_s)
427	0	call Bchksum(visc%Kv_shear_Bu, "after calc_KS_vert visc%Kv_shear_Bu", G%HI, unscale=GV%HZ_T_to_m2_s)
428	0	call Bchksum(visc%TKE_turb, "after calc_KS_vert visc%TKE_turb", G%HI, unscale=US%Z_to_m*2US%s_to_T**2)
429		endif
430		else
431		! Changes: visc%Kd_shear ; Sets: visc%Kv_shear and visc%TKE_turb
432		call calculate_kappa_shear(u_h, v_h, h, tv, fluxes%p_surf, visc%Kd_shear, visc%TKE_turb, &
433	0	visc%Kv_shear, dt, G, GV, US, CS%kappaShear_CSp)
434	0	if (CS%debug) then
435	0	call hchksum(visc%Kd_shear, "after calc_KS visc%Kd_shear", G%HI, unscale=GV%HZ_T_to_m2_s)
436	0	call hchksum(visc%Kv_shear, "after calc_KS visc%Kv_shear", G%HI, unscale=GV%HZ_T_to_m2_s)
437	0	call hchksum(visc%TKE_turb, "after calc_KS visc%TKE_turb", G%HI, unscale=US%Z_to_m*2US%s_to_T**2)
438		endif
439		endif
440	12	call cpu_clock_end(id_clock_kappaShear)
441	12	if (showCallTree) call callTree_waypoint("done with calculate_kappa_shear (set_diffusivity)")
442	12	if (associated(VBF%Kd_KS)) then ; do K=1,nz+1 ; do i=is,ie ; do j=js,je
443	0	VBF%Kd_KS(i,j,K) = visc%Kd_shear(i,j,K)
444		enddo ; enddo ; enddo ; endif
445	0	elseif (CS%use_CVMix_shear) then
446		!NOTE{BGR}: this needs to be cleaned up. It works in 1D case, but has not been tested outside.
447	0	call calculate_CVMix_shear(u_h, v_h, h, tv, visc%Kd_shear, visc%Kv_shear, G, GV, US, CS%CVMix_shear_CSp)
448	0	if (CS%debug) then
449	0	call hchksum(visc%Kd_shear, "after CVMix_shear visc%Kd_shear", G%HI, unscale=GV%HZ_T_to_m2_s)
450	0	call hchksum(visc%Kv_shear, "after CVMix_shear visc%Kv_shear", G%HI, unscale=GV%HZ_T_to_m2_s)
451		endif
452	0	elseif (associated(visc%Kv_shear)) then
453	0	visc%Kv_shear(:,:,:) = 0.0 ! needed if calculate_kappa_shear is not enabled
454		endif
455
456		! Smooth the properties through massless layers.
457	12	if (use_EOS) then
458	12	if (CS%debug) then
459	0	call hchksum(tv%T, "before vert_fill_TS tv%T", G%HI, unscale=US%C_to_degC)
460	0	call hchksum(tv%S, "before vert_fill_TS tv%S", G%HI, unscale=US%S_to_ppt)
461	0	call hchksum(h, "before vert_fill_TS h",G%HI, unscale=GV%H_to_m)
462		endif
463		!$omp target enter data map(to: tv%T, tv%S, h)
464		!$omp target enter data map(alloc: T_f, S_f)
465	12	call vert_fill_TS(h, tv%T, tv%S, kappa_dt_fill, T_f, S_f, G, GV, US, larger_h_denom=.true.)
466		!$omp target exit data map(from: T_f, S_f)
467		!$omp target exit data map(release: tv%T, tv%S, h)
468	12	if (CS%debug) then
469	0	call hchksum(tv%T, "after vert_fill_TS tv%T", G%HI, unscale=US%C_to_degC)
470	0	call hchksum(tv%S, "after vert_fill_TS tv%S", G%HI, unscale=US%S_to_ppt)
471	0	call hchksum(h, "after vert_fill_TS h",G%HI, unscale=GV%H_to_m)
472		endif
473		endif
474
475		! Calculate the diffusivities, Kd_lay and Kd_int, for each layer and interface. This would
476		! be an appropriate place to add a depth-dependent parameterization or another explicit
477		! parameterization of Kd.
478
479		!$OMP parallel do default(shared) private(dRho_int,N2_lay,Kd_lay_2d,Kd_int_2d,Kv_bkgnd,N2_int,dz, &
480		!$OMP N2_bot,rho_bot,h_bot,k_bot,KT_extra,KS_extra,TKE_to_Kd,maxTKE,dissip,kb) &
481		!$OMP if(.not. CS%use_CVMix_ddiff)
482	732	do j=js,je
483
484		! Set up variables related to the stratification.
485	720	call find_N2(h, tv, T_f, S_f, fluxes, j, G, GV, US, CS, dRho_int, N2_lay, N2_int, N2_bot, rho_bot, h_bot, k_bot)
486
487	720	if (associated(dd%N2_3d)) then
488	0	do K=1,nz+1 ; do i=is,ie ; dd%N2_3d(i,j,K) = N2_int(i,K) ; enddo ; enddo
489		endif
490
491		! Add background mixing
492	720	call calculate_bkgnd_mixing(h, tv, N2_lay, Kd_lay_2d, Kd_int_2d, Kv_bkgnd, j, G, GV, US, CS%bkgnd_mixing_csp)
493		! Update Kv and 3-d diffusivity diagnostics.
494	720	if (associated(visc%Kv_slow)) then ; do K=1,nz+1 ; do i=is,ie
495	0	visc%Kv_slow(i,j,K) = visc%Kv_slow(i,j,K) + Kv_bkgnd(i,K)
496		enddo ; enddo ; endif
497	720	if (CS%id_Kv_bkgnd > 0) then ; do K=1,nz+1 ; do i=is,ie
498	0	dd%Kv_bkgnd(i,j,K) = Kv_bkgnd(i,K)
499		enddo ; enddo ; endif
500	720	if (CS%id_Kd_bkgnd > 0) then ; do K=1,nz+1 ; do i=is,ie
501	0	dd%Kd_bkgnd(i,j,K) = Kd_int_2d(i,K)
502		enddo ; enddo ; endif
503	720	if (associated(VBF%Kd_bkgnd)) then ; do K=1,nz+1 ; do i=is,ie
504	0	VBF%Kd_bkgnd(i,j,K) = Kd_int_2d(i,K)
505		enddo ; enddo ; endif
506
507		! Double-diffusion (old method)
508	720	if (CS%double_diffusion) then
509	0	call double_diffusion(tv, h, T_f, S_f, j, G, GV, US, CS, KT_extra, KS_extra)
510		! One of Kd_extra_T and Kd_extra_S is always 0. Kd_extra_S is positive for salt fingering.
511		! Kd_extra_T is positive for double diffusive convection.
512	0	do K=2,nz ; do i=is,ie
513	0	if (KS_extra(i,K) > KT_extra(i,K)) then ! salt fingering
514	0	Kd_lay_2d(i,k-1) = Kd_lay_2d(i,k-1) + 0.5 * KT_extra(i,K)
515	0	Kd_lay_2d(i,k) = Kd_lay_2d(i,k) + 0.5 * KT_extra(i,K)
516	0	Kd_extra_S(i,j,K) = KS_extra(i,K) - KT_extra(i,K)
517	0	Kd_extra_T(i,j,K) = 0.0
518	0	elseif (KT_extra(i,K) > 0.0) then ! double-diffusive convection
519	0	Kd_lay_2d(i,k-1) = Kd_lay_2d(i,k-1) + 0.5 * KS_extra(i,K)
520	0	Kd_lay_2d(i,k) = Kd_lay_2d(i,k) + 0.5 * KS_extra(i,K)
521	0	Kd_extra_T(i,j,K) = KT_extra(i,K) - KS_extra(i,K)
522	0	Kd_extra_S(i,j,K) = 0.0
523		else ! There is no double diffusion at this interface.
524	0	Kd_extra_T(i,j,K) = 0.0
525	0	Kd_extra_S(i,j,K) = 0.0
526		endif
527		enddo ; enddo
528	0	if (associated(dd%KT_extra)) then ; do K=1,nz+1 ; do i=is,ie
529	0	dd%KT_extra(i,j,K) = KT_extra(i,K)
530		enddo ; enddo ; endif
531
532	0	if (associated(dd%KS_extra)) then ; do K=1,nz+1 ; do i=is,ie
533	0	dd%KS_extra(i,j,K) = KS_extra(i,K)
534		enddo ; enddo ; endif
535
536	0	if (associated(VBF%Kd_ddiff_T)) then ; do K=1,nz+1 ; do i=is,ie
537	0	VBF%Kd_ddiff_T(i,j,K) = KT_extra(i,K)
538		enddo ; enddo ; endif
539	0	if (associated(VBF%Kd_ddiff_S)) then ; do K=1,nz+1 ; do i=is,ie
540	0	VBF%Kd_ddiff_S(i,j,K) = KS_extra(i,K)
541		enddo ; enddo ; endif
542		endif
543
544		! Apply double diffusion via CVMix
545		! GMM, we need to pass HBL to compute_ddiff_coeffs, but it is not yet available.
546	720	if (CS%use_CVMix_ddiff) then
547	0	call cpu_clock_begin(id_clock_CVMix_ddiff)
548	0	if (associated(dd%drho_rat)) then
549		call compute_ddiff_coeffs(h, tv, G, GV, US, j, Kd_extra_T, Kd_extra_S, &
550	0	CS%CVMix_ddiff_csp, dd%drho_rat)
551		else
552	0	call compute_ddiff_coeffs(h, tv, G, GV, US, j, Kd_extra_T, Kd_extra_S, CS%CVMix_ddiff_csp)
553		endif
554	0	if (associated(VBF%Kd_ddiff_T)) then ; do K=1,nz+1 ; do i=is,ie
555	0	VBF%Kd_ddiff_T(i,j,K) = KT_extra(i,K)
556		enddo ; enddo ; endif
557	0	if (associated(VBF%Kd_ddiff_S)) then ; do K=1,nz+1 ; do i=is,ie
558	0	VBF%Kd_ddiff_S(i,j,K) = KS_extra(i,K)
559		enddo ; enddo ; endif
560	0	call cpu_clock_end(id_clock_CVMix_ddiff)
561		endif
562
563		! Calculate conversion ratios from TKE to layer diffusivities.
564	720	if (TKE_to_Kd_used) then
565	0	call find_TKE_to_Kd(h, tv, dRho_int, N2_lay, j, dt, G, GV, US, CS, TKE_to_Kd, maxTKE, kb)
566	0	if (associated(dd%maxTKE)) then ; do k=1,nz ; do i=is,ie
567	0	dd%maxTKE(i,j,k) = maxTKE(i,k)
568		enddo ; enddo ; endif
569	0	if (associated(dd%TKE_to_Kd)) then ; do k=1,nz ; do i=is,ie
570	0	dd%TKE_to_Kd(i,j,k) = TKE_to_Kd(i,k)
571		enddo ; enddo ; endif
572		endif
573
574		! Add the input turbulent diffusivity.
575	720	if (CS%useKappaShear .or. CS%use_CVMix_shear) then
576	6447600	do K=2,nz ; do i=is,ie
577	6446880	Kd_int_2d(i,K) = visc%Kd_shear(i,j,K) + 0.5 * (Kd_lay_2d(i,k-1) + Kd_lay_2d(i,k))
578		enddo ; enddo
579	87120	do i=is,ie
580	86400	Kd_int_2d(i,1) = visc%Kd_shear(i,j,1) ! This isn't actually used. It could be 0.
581	87120	Kd_int_2d(i,nz+1) = 0.0
582		enddo
583	6534720	do k=1,nz ; do i=is,ie
584	6534000	Kd_lay_2d(i,k) = Kd_lay_2d(i,k) + 0.5 * (visc%Kd_shear(i,j,K) + visc%Kd_shear(i,j,K+1))
585		enddo ; enddo
586		else
587	0	do i=is,ie
588	0	Kd_int_2d(i,1) = Kd_lay_2d(i,1) ; Kd_int_2d(i,nz+1) = 0.0
589		enddo
590	0	do K=2,nz ; do i=is,ie
591	0	Kd_int_2d(i,K) = 0.5 * (Kd_lay_2d(i,k-1) + Kd_lay_2d(i,k))
592		enddo ; enddo
593		endif
594
595	720	if (CS%ML_radiation .or. CS%use_tidal_mixing .or. associated(dd%Kd_Work)) then
596	0	call thickness_to_dz(h, tv, dz, j, G, GV)
597		endif
598
599		! Add the ML_Rad diffusivity.
600	720	if (CS%ML_radiation) then
601	0	call add_MLrad_diffusivity(dz, fluxes, tv, j, Kd_int_2d, G, GV, US, CS, TKE_to_Kd, Kd_lay_2d)
602		endif
603
604		! Add the Nikurashin and / or tidal bottom-driven mixing
605	720	if (CS%use_tidal_mixing) &
606		call calculate_tidal_mixing(dz, j, N2_bot, rho_bot, N2_lay, N2_int, TKE_to_Kd, &
607		maxTKE, G, GV, US, CS%tidal_mixing, &
608	0	CS%Kd_max, visc%Kv_slow, Kd_lay_2d, Kd_int_2d, VBF)
609
610		! Add diffusivity from internal tides ray tracing
611	720	if (CS%use_int_tides) then
612
613		call get_lowmode_diffusivity(G, GV, h, tv, US, h_bot, k_bot, j, N2_lay, N2_int, TKE_to_Kd, CS%Kd_max, &
614		CS%int_tide_CSp, Kd_leak_2d, Kd_quad_2d, Kd_itidal_2d, Kd_Froude_2d, Kd_slope_2d, &
615		Kd_lay_2d, Kd_int_2d, prof_leak_2d, prof_quad_2d, prof_itidal_2d, prof_froude_2d, &
616	0	prof_slope_2d)
617
618	0	if (CS%id_kbbl > 0) then ; do i=is,ie
619	0	dd%kbbl(i,j) = k_bot(i)
620		enddo ; endif
621	0	if (CS%id_bbl_thick > 0) then ; do i=is,ie
622	0	dd%bbl_thick(i,j) = h_bot(i)
623		enddo ; endif
624	0	if (CS%id_Kd_leak > 0) then ; do K=1,nz+1 ; do i=is,ie
625	0	dd%Kd_leak(i,j,K) = Kd_leak_2d(i,K)
626		enddo ; enddo ; endif
627	0	if (CS%id_Kd_quad > 0) then ; do K=1,nz+1 ; do i=is,ie
628	0	dd%Kd_quad(i,j,K) = Kd_quad_2d(i,K)
629		enddo ; enddo ; endif
630	0	if (CS%id_Kd_itidal > 0) then ; do K=1,nz+1 ; do i=is,ie
631	0	dd%Kd_itidal(i,j,K) = Kd_itidal_2d(i,K)
632		enddo ; enddo ; endif
633	0	if (CS%id_Kd_Froude > 0) then ; do K=1,nz+1 ; do i=is,ie
634	0	dd%Kd_Froude(i,j,K) = Kd_Froude_2d(i,K)
635		enddo ; enddo ; endif
636	0	if (CS%id_Kd_slope > 0) then ; do K=1,nz+1 ; do i=is,ie
637	0	dd%Kd_slope(i,j,K) = Kd_slope_2d(i,K)
638		enddo ; enddo ; endif
639	0	if (associated (VBF%Kd_leak)) then ; do K=1,nz+1 ; do i=is,ie
640	0	VBF%Kd_leak(i,j,K) = min(Kd_leak_2d(i,K), CS%Kd_max)
641		enddo ; enddo ; endif
642	0	if (associated (VBF%Kd_quad)) then ; do K=1,nz+1 ; do i=is,ie
643	0	VBF%Kd_quad(i,j,K) = min(Kd_quad_2d(i,K), CS%Kd_max)
644		enddo ; enddo ; endif
645	0	if (associated (VBF%Kd_itidal)) then ; do K=1,nz+1 ; do i=is,ie
646	0	VBF%Kd_itidal(i,j,K) = min(Kd_itidal_2d(i,K), CS%Kd_max)
647		enddo ; enddo ; endif
648	0	if (associated (VBF%Kd_Froude)) then ; do K=1,nz+1 ; do i=is,ie
649	0	VBF%Kd_Froude(i,j,K) = min(Kd_Froude_2d(i,K), CS%Kd_max)
650		enddo ; enddo ; endif
651	0	if (associated (VBF%Kd_slope)) then ; do K=1,nz+1 ; do i=is,ie
652	0	VBF%Kd_slope(i,j,K) = min(Kd_slope_2d(i,K), CS%Kd_max)
653		enddo ; enddo ; endif
654
655	0	if (CS%id_prof_leak > 0) then ; do k=1,nz ; do i=is,ie
656	0	dd%prof_leak(i,j,k) = prof_leak_2d(i,k)
657		enddo ; enddo ; endif
658	0	if (CS%id_prof_quad > 0) then ; do k=1,nz ; do i=is,ie
659	0	dd%prof_quad(i,j,k) = prof_quad_2d(i,k)
660		enddo ; enddo ; endif
661	0	if (CS%id_prof_itidal > 0) then ; do k=1,nz ; do i=is,ie
662	0	dd%prof_itidal(i,j,k) = prof_itidal_2d(i,k)
663		enddo ; enddo ; endif
664	0	if (CS%id_prof_Froude > 0) then ; do k=1,nz ; do i=is,ie
665	0	dd%prof_Froude(i,j,k) = prof_Froude_2d(i,k)
666		enddo ; enddo ; endif
667	0	if (CS%id_prof_slope > 0) then ; do k=1,nz ; do i=is,ie
668	0	dd%prof_slope(i,j,k) = prof_slope_2d(i,k)
669		enddo ; enddo ; endif
670		endif
671
672		! This adds the diffusion sustained by the energy extracted from the flow by the bottom drag.
673	720	if (CS%bottomdraglaw .and. (CS%BBL_effic > 0.0)) then
674	720	if (CS%use_LOTW_BBL_diffusivity) then
675		call add_LOTW_BBL_diffusivity(h, u, v, tv, fluxes, visc, j, N2_int, Rho_bot, Kd_int_2d, &
676	720	G, GV, US, CS, dd%Kd_BBL, Kd_lay_2d)
677		else
678		call add_drag_diffusivity(h, u, v, tv, fluxes, visc, j, TKE_to_Kd, &
679	0	maxTKE, kb, rho_bot, G, GV, US, CS, Kd_lay_2d, Kd_int_2d, dd%Kd_BBL)
680		endif
681	720	if (associated(VBF%Kd_BBL)) then ; do K=1,nz+1 ; do i=is,ie
682	0	VBF%Kd_BBL(i,j,K) = dd%Kd_BBL(i,j,K)
683		enddo ; enddo ; endif
684		endif
685
686	720	if (CS%limit_dissipation) then
687		! This calculates the dissipation ONLY from Kd calculated in this routine
688		! dissip has units of W/m3 (= kg/m3 * m2/s * 1/s2)
689		! 1) a global constant,
690		! 2) a dissipation proportional to N (aka Gargett) and
691		! 3) dissipation corresponding to a (nearly) constant diffusivity.
692	0	do K=2,nz ; do i=is,ie
693		dissip = max( CS%dissip_min, & ! Const. floor on dissip.
694		CS%dissip_N0 + CS%dissip_N1 * sqrt(N2_int(i,K)), & ! Floor aka Gargett
695	0	CS%dissip_N2 * N2_int(i,K)) ! Floor of Kd_min*rho0/F_Ri
696		Kd_int_2d(i,K) = max(Kd_int_2d(i,K) , & ! Apply floor to Kd
697	0	dissip * (CS%FluxRi_max / (GV%H_to_RZ * (N2_int(i,K) + Omega2))))
698		enddo ; enddo
699		endif
700
701		! Optionally add a uniform diffusivity at the interfaces.
702	720	if (CS%Kd_add > 0.0) then
703	0	do K=1,nz+1 ; do i=is,ie
704	0	Kd_int_2d(i,K) = Kd_int_2d(i,K) + CS%Kd_add
705		enddo ; enddo
706	0	VBF%Kd_add = CS%Kd_add
707		endif
708
709		! Copy the 2-d slices into the 3-d array that is exported.
710	6621840	do K=1,nz+1 ; do i=is,ie
711	6621120	Kd_int(i,j,K) = Kd_int_2d(i,K)
712		enddo ; enddo
713
714	720	if (CS%limit_dissipation) then
715		! This calculates the layer dissipation ONLY from Kd calculated in this routine
716		! dissip has units of W/m3 (= kg/m3 * m2/s * 1/s2)
717		! 1) a global constant,
718		! 2) a dissipation proportional to N (aka Gargett) and
719		! 3) dissipation corresponding to a (nearly) constant diffusivity.
720	0	do k=2,nz-1 ; do i=is,ie
721		dissip = max( CS%dissip_min, & ! Const. floor on dissip.
722		CS%dissip_N0 + CS%dissip_N1 * sqrt(N2_lay(i,k)), & ! Floor aka Gargett
723	0	CS%dissip_N2 * N2_lay(i,k)) ! Floor of Kd_min*rho0/F_Ri
724		Kd_lay_2d(i,k) = max(Kd_lay_2d(i,k) , & ! Apply floor to Kd
725	0	dissip * (CS%FluxRi_max / (GV%H_to_RZ * (N2_lay(i,k) + Omega2))))
726		enddo ; enddo
727		endif
728
729	720	if (associated(dd%Kd_Work)) then
730	0	do k=1,nz ; do i=is,ie
731	0	dd%Kd_Work(i,j,k) = GV%H_to_RZ * Kd_lay_2d(i,k) * N2_lay(i,k) * dz(i,k) ! Watt m-2 = kg s-3
732		enddo ; enddo
733		endif
734
735		! Optionally add a uniform diffusivity to the layers.
736	720	if ((CS%Kd_add > 0.0) .and. (present(Kd_lay))) then
737	0	do k=1,nz ; do i=is,ie
738	0	Kd_lay_2d(i,k) = Kd_lay_2d(i,k) + CS%Kd_add
739		enddo ; enddo
740		endif
741
742	720	if (associated(dd%Kd_Work_added)) then
743	0	do k=1,nz ; do i=is,ie
744	0	dd%Kd_Work_added(i,j,k) = GV%H_to_RZ * CS%Kd_add * N2_lay(i,k) * dz(i,k) ! Watt m-2 = kg s-3
745		enddo ; enddo
746		endif
747
748		! Copy the 2-d slices into the 3-d array that is exported; this was done above for Kd_int.
749	732	if (present(Kd_lay)) then ; do k=1,nz ; do i=is,ie
750	0	Kd_lay(i,j,k) = Kd_lay_2d(i,k)
751		enddo ; enddo ; endif
752		enddo ! j-loop
753
754	12	if (CS%user_change_diff) then
755		call user_change_diff(h, tv, G, GV, US, CS%user_change_diff_CSp, Kd_lay, Kd_int, &
756	0	T_f, S_f, dd%Kd_user)
757		endif
758
759	12	if (CS%debug) then
760	0	if (present(Kd_lay)) call hchksum(Kd_lay, "Kd_lay", G%HI, haloshift=0, unscale=GV%HZ_T_to_m2_s)
761
762	0	if (CS%useKappaShear) call hchksum(visc%Kd_shear, "Turbulent Kd", G%HI, haloshift=0, unscale=GV%HZ_T_to_m2_s)
763
764	0	if (CS%use_CVMix_ddiff) then
765	0	call hchksum(Kd_extra_T, "MOM_set_diffusivity: Kd_extra_T", G%HI, haloshift=0, unscale=GV%HZ_T_to_m2_s)
766	0	call hchksum(Kd_extra_S, "MOM_set_diffusivity: Kd_extra_S", G%HI, haloshift=0, unscale=GV%HZ_T_to_m2_s)
767		endif
768
769	0	if (allocated(visc%kv_bbl_u) .and. allocated(visc%kv_bbl_v)) then
770		call uvchksum("BBL Kv_bbl_[uv]", visc%kv_bbl_u, visc%kv_bbl_v, G%HI, &
771		haloshift=0, symmetric=.true., unscale=GV%HZ_T_to_m2_s, &
772	0	scalar_pair=.true.)
773		endif
774
775	0	if (allocated(visc%bbl_thick_u) .and. allocated(visc%bbl_thick_v)) then
776		call uvchksum("BBL bbl_thick_[uv]", visc%bbl_thick_u, visc%bbl_thick_v, &
777		G%HI, haloshift=0, symmetric=.true., unscale=US%Z_to_m, &
778	0	scalar_pair=.true.)
779		endif
780
781	0	if (allocated(visc%Ray_u) .and. allocated(visc%Ray_v)) then
782		call uvchksum("Ray_[uv]", visc%Ray_u, visc%Ray_v, G%HI, 0, &
783	0	symmetric=.true., unscale=GV%H_to_m*US%s_to_T, scalar_pair=.true.)
784		endif
785
786		endif
787
788	12	if (CS%debug) then
789	0	if (CS%id_prof_leak > 0) call hchksum(dd%prof_leak, "leakage_profile", G%HI, haloshift=0, unscale=GV%m_to_H)
790	0	if (CS%id_prof_slope > 0) call hchksum(dd%prof_slope, "slope_profile", G%HI, haloshift=0, unscale=GV%m_to_H)
791	0	if (CS%id_prof_Froude > 0) call hchksum(dd%prof_Froude, "Froude_profile", G%HI, haloshift=0, unscale=GV%m_to_H)
792	0	if (CS%id_prof_quad > 0) call hchksum(dd%prof_quad, "quad_profile", G%HI, haloshift=0, unscale=GV%m_to_H)
793	0	if (CS%id_prof_itidal > 0) call hchksum(dd%prof_itidal, "itidal_profile", G%HI, haloshift=0, unscale=GV%m_to_H)
794	0	if (CS%id_TKE_to_Kd > 0) call hchksum(dd%TKE_to_Kd, "TKE_to_Kd", G%HI, haloshift=0, unscale=US%m_to_ZUS%T_to_s*2)
795	0	if (CS%id_Kd_leak > 0) call hchksum(dd%Kd_leak, "Kd_leak", G%HI, haloshift=0, unscale=GV%HZ_T_to_m2_s)
796	0	if (CS%id_Kd_quad > 0) call hchksum(dd%Kd_quad, "Kd_quad", G%HI, haloshift=0, unscale=GV%HZ_T_to_m2_s)
797	0	if (CS%id_Kd_itidal > 0) call hchksum(dd%Kd_itidal, "Kd_itidal", G%HI, haloshift=0, unscale=GV%HZ_T_to_m2_s)
798	0	if (CS%id_Kd_Froude > 0) call hchksum(dd%Kd_Froude, "Kd_Froude", G%HI, haloshift=0, unscale=GV%HZ_T_to_m2_s)
799	0	if (CS%id_Kd_slope > 0) call hchksum(dd%Kd_slope, "Kd_slope", G%HI, haloshift=0, unscale=GV%HZ_T_to_m2_s)
800		endif
801
802		! post diagnostics
803	12	if (present(Kd_lay) .and. (CS%id_Kd_layer > 0)) call post_data(CS%id_Kd_layer, Kd_lay, CS%diag)
804
805		! background mixing
806	12	if (CS%id_Kd_bkgnd > 0) call post_data(CS%id_Kd_bkgnd, dd%Kd_bkgnd, CS%diag)
807	12	if (CS%id_Kv_bkgnd > 0) call post_data(CS%id_Kv_bkgnd, dd%Kv_bkgnd, CS%diag)
808
809	12	if (CS%id_kbbl > 0) call post_data(CS%id_kbbl, dd%kbbl, CS%diag)
810	12	if (CS%id_bbl_thick > 0) call post_data(CS%id_bbl_thick, dd%bbl_thick, CS%diag)
811	12	if (CS%id_Kd_leak > 0) call post_data(CS%id_Kd_leak, dd%Kd_leak, CS%diag)
812	12	if (CS%id_Kd_slope > 0) call post_data(CS%id_Kd_slope, dd%Kd_slope, CS%diag)
813	12	if (CS%id_Kd_Froude > 0) call post_data(CS%id_Kd_Froude, dd%Kd_Froude, CS%diag)
814	12	if (CS%id_Kd_quad > 0) call post_data(CS%id_Kd_quad, dd%Kd_quad, CS%diag)
815	12	if (CS%id_Kd_itidal > 0) call post_data(CS%id_Kd_itidal, dd%Kd_itidal, CS%diag)
816
817	12	if (CS%id_prof_leak > 0) call post_data(CS%id_prof_leak, dd%prof_leak, CS%diag)
818	12	if (CS%id_prof_slope > 0) call post_data(CS%id_prof_slope, dd%prof_slope, CS%diag)
819	12	if (CS%id_prof_Froude > 0) call post_data(CS%id_prof_Froude, dd%prof_Froude, CS%diag)
820	12	if (CS%id_prof_quad > 0) call post_data(CS%id_prof_quad, dd%prof_quad, CS%diag)
821	12	if (CS%id_prof_itidal > 0) call post_data(CS%id_prof_itidal, dd%prof_itidal, CS%diag)
822
823		! tidal mixing
824	12	if (CS%use_tidal_mixing) &
825	0	call post_tidal_diagnostics(G, GV, h, CS%tidal_mixing)
826
827	12	if (CS%id_N2 > 0) call post_data(CS%id_N2, dd%N2_3d, CS%diag)
828	12	if (CS%id_Kd_Work > 0) call post_data(CS%id_Kd_Work, dd%Kd_Work, CS%diag)
829	12	if (CS%id_Kd_Work_added > 0) call post_data(CS%id_Kd_Work_added, dd%Kd_Work_added, CS%diag)
830	12	if (CS%id_maxTKE > 0) call post_data(CS%id_maxTKE, dd%maxTKE, CS%diag)
831	12	if (CS%id_TKE_to_Kd > 0) call post_data(CS%id_TKE_to_Kd, dd%TKE_to_Kd, CS%diag)
832
833	12	if (CS%id_Kd_user > 0) call post_data(CS%id_Kd_user, dd%Kd_user, CS%diag)
834
835		! double diffusive mixing
836	12	if (CS%double_diffusion) then
837	0	if (CS%id_KT_extra > 0) call post_data(CS%id_KT_extra, dd%KT_extra, CS%diag)
838	0	if (CS%id_KS_extra > 0) call post_data(CS%id_KS_extra, dd%KS_extra, CS%diag)
839	12	elseif (CS%use_CVMix_ddiff) then
840	0	if (CS%id_KT_extra > 0) call post_data(CS%id_KT_extra, Kd_extra_T, CS%diag)
841	0	if (CS%id_KS_extra > 0) call post_data(CS%id_KS_extra, Kd_extra_S, CS%diag)
842	0	if (CS%id_R_rho > 0) call post_data(CS%id_R_rho, dd%drho_rat, CS%diag)
843		endif
844	12	if (CS%id_Kd_BBL > 0) call post_data(CS%id_Kd_BBL, dd%Kd_BBL, CS%diag)
845
846	12	if (associated(dd%N2_3d)) deallocate(dd%N2_3d)
847	12	if (associated(dd%Kd_Work)) deallocate(dd%Kd_Work)
848	12	if (associated(dd%Kd_Work_added)) deallocate(dd%Kd_Work_added)
849	12	if (associated(dd%Kd_user)) deallocate(dd%Kd_user)
850	12	if (associated(dd%maxTKE)) deallocate(dd%maxTKE)
851	12	if (associated(dd%TKE_to_Kd)) deallocate(dd%TKE_to_Kd)
852	12	if (associated(dd%KT_extra)) deallocate(dd%KT_extra)
853	12	if (associated(dd%KS_extra)) deallocate(dd%KS_extra)
854	12	if (associated(dd%drho_rat)) deallocate(dd%drho_rat)
855	12	if (associated(dd%Kd_BBL)) deallocate(dd%Kd_BBL)
856	12	if (associated(dd%Kd_bkgnd)) deallocate(dd%Kd_bkgnd)
857	12	if (associated(dd%Kv_bkgnd)) deallocate(dd%Kv_bkgnd)
858
859	12	if (associated(dd%Kd_leak)) deallocate(dd%Kd_leak)
860	12	if (associated(dd%Kd_quad)) deallocate(dd%Kd_quad)
861	12	if (associated(dd%Kd_itidal)) deallocate(dd%Kd_itidal)
862	12	if (associated(dd%Kd_Froude)) deallocate(dd%Kd_Froude)
863	12	if (associated(dd%Kd_slope)) deallocate(dd%Kd_slope)
864	12	if (associated(dd%prof_leak)) deallocate(dd%prof_leak)
865	12	if (associated(dd%prof_quad)) deallocate(dd%prof_quad)
866	12	if (associated(dd%prof_itidal)) deallocate(dd%prof_itidal)
867	12	if (associated(dd%prof_Froude)) deallocate(dd%prof_Froude)
868	12	if (associated(dd%prof_slope)) deallocate(dd%prof_slope)
869
870	12	if (showCallTree) call callTree_leave("set_diffusivity()")
871
872	12	end subroutine set_diffusivity
873
874		!> Convert turbulent kinetic energy to diffusivity
875	0	subroutine find_TKE_to_Kd(h, tv, dRho_int, N2_lay, j, dt, G, GV, US, CS, &
876	0	TKE_to_Kd, maxTKE, kb)
877		type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure
878		type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure
879		type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
880		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
881		intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2]
882		type(thermo_var_ptrs), intent(in) :: tv !< Structure containing pointers to any available
883		!! thermodynamic fields.
884		real, dimension(SZI_(G),SZK_(GV)+1), intent(in) :: dRho_int !< Change in locally referenced potential density
885		!! across each interface [R ~> kg m-3].
886		real, dimension(SZI_(G),SZK_(GV)), intent(in) :: N2_lay !< The squared buoyancy frequency of the
887		!! layers [T-2 ~> s-2].
888		integer, intent(in) :: j !< j-index of row to work on
889		real, intent(in) :: dt !< Time increment [T ~> s].
890		type(set_diffusivity_CS), pointer :: CS !< Diffusivity control structure
891		real, dimension(SZI_(G),SZK_(GV)), intent(out) :: TKE_to_Kd !< The conversion rate between the
892		!! TKE dissipated within a layer and the
893		!! diapycnal diffusivity within that layer,
894		!! usually (~Rho_0 / (G_Earth * dRho_lay))
895		!! [T2 Z-1 ~> s2 m-1]
896		real, dimension(SZI_(G),SZK_(GV)), intent(out) :: maxTKE !< The energy required to for a layer to entrain to its
897		!! maximum realizable thickness [H Z2 T-3 ~> m3 s-3 or W m-2]
898		integer, dimension(SZI_(G)), intent(out) :: kb !< Index of lightest layer denser than the buffer
899		!! layer, or -1 without a bulk mixed layer.
900		! Local variables
901		real, dimension(SZI_(G),SZK_(GV)) :: &
902	0	ds_dsp1, & ! coordinate variable (sigma-2) difference across an
903		! interface divided by the difference across the interface
904		! below it [nondim]
905	0	dsp1_ds, & ! inverse coordinate variable (sigma-2) difference
906		! across an interface times the difference across the
907		! interface above it [nondim]
908	0	rho_0, & ! Layer potential densities relative to surface pressure [R ~> kg m-3]
909	0	dz, & ! Height change across layers [Z ~> m]
910	0	maxEnt ! maxEnt is the maximum value of entrainment from below (with
911		! compensating entrainment from above to keep the layer
912		! density from changing) that will not deplete all of the
913		! layers above or below a layer within a timestep [H ~> m or kg m-2].
914		real, dimension(SZI_(G)) :: &
915	0	htot, & ! total thickness above or below a layer, or the
916		! integrated thickness in the BBL [H ~> m or kg m-2].
917	0	mFkb, & ! total thickness in the mixed and buffer layers times ds_dsp1 [H ~> m or kg m-2]
918	0	p_ref, & ! array of tv%P_Ref pressures [R L2 T-2 ~> Pa]
919	0	Rcv_kmb, & ! coordinate density in the lowest buffer layer [R ~> kg m-3]
920	0	p_0 ! An array of 0 pressures [R L2 T-2 ~> Pa]
921
922		real :: dh_max ! maximum amount of entrainment a layer could undergo before
923		! entraining all fluid in the layers above or below [H ~> m or kg m-2]
924		real :: dRho_lay ! density change across a layer [R ~> kg m-3]
925		real :: Omega2 ! rotation rate squared [T-2 ~> s-2]
926		real :: grav ! Gravitational acceleration [Z T-2 ~> m s-2]
927		real :: G_Rho0 ! Gravitational acceleration divided by Boussinesq reference density
928		! [Z R-1 T-2 ~> m4 s-2 kg-1]
929		real :: G_IRho0 ! Alternate calculation of G_Rho0 with thickness rescaling factors
930		! [Z2 T-2 R-1 H-1 ~> m4 s-2 kg-1 or m7 kg-2 s-2]
931		real :: I_dt ! 1/dt [T-1 ~> s-1]
932		real :: dz_neglect ! A negligibly small height change [Z ~> m]
933		real :: hN2pO2 ! h (N^2 + Omega^2), in [Z T-2 ~> m s-2].
934	0	logical :: do_i(SZI_(G))
935
936		integer, dimension(2) :: EOSdom ! The i-computational domain for the equation of state
937		integer :: i, k, is, ie, nz, i_rem, kmb, kb_min
938	0	is = G%isc ; ie = G%iec ; nz = GV%ke
939
940	0	I_dt = 1.0 / dt
941	0	Omega2 = CS%omega**2
942	0	dz_neglect = GV%dZ_subroundoff
943	0	grav = GV%g_Earth_Z_T2
944	0	G_Rho0 = grav / GV%Rho0
945	0	if (CS%answer_date < 20190101) then
946	0	G_IRho0 = grav * GV%H_to_Z*2 GV%RZ_to_H
947		else
948	0	G_IRho0 = GV%H_to_Z*G_Rho0
949		endif
950
951		! Find the vertical distances across layers.
952	0	call thickness_to_dz(h, tv, dz, j, G, GV)
953
954		! Simple but coordinate-independent estimate of Kd/TKE
955	0	if (CS%simple_TKE_to_Kd) then
956	0	do k=1,nz ; do i=is,ie
957	0	hN2pO2 = dz(i,k) * (N2_lay(i,k) + Omega2) ! Units of Z T-2.
958	0	if (hN2pO2 > 0.) then
959	0	TKE_to_Kd(i,k) = 1.0 / hN2pO2 ! Units of T2 H-1.
960	0	else ; TKE_to_Kd(i,k) = 0. ; endif
961		! The maximum TKE conversion we allow is really a statement
962		! about the upper diffusivity we allow. Kd_max must be set.
963	0	maxTKE(i,k) = hN2pO2 * CS%Kd_max ! Units of H Z2 T-3.
964		enddo ; enddo
965	0	kb(is:ie) = -1 ! kb should not be used by any code in non-layered mode -AJA
966	0	return
967		endif
968
969		! Determine kb - the index of the shallowest active interior layer.
970	0	if (CS%bulkmixedlayer) then
971	0	kmb = GV%nk_rho_varies
972	0	do i=is,ie ; p_0(i) = 0.0 ; p_ref(i) = tv%P_Ref ; enddo
973	0	EOSdom(:) = EOS_domain(G%HI)
974	0	do k=1,nz
975	0	call calculate_density(tv%T(:,j,k), tv%S(:,j,k), p_0, rho_0(:,k), tv%eqn_of_state, EOSdom)
976		enddo
977	0	call calculate_density(tv%T(:,j,kmb), tv%S(:,j,kmb), p_ref, Rcv_kmb, tv%eqn_of_state, EOSdom)
978
979	0	kb_min = kmb+1
980	0	do i=is,ie
981		! Determine the next denser layer than the buffer layer in the
982		! coordinate density (sigma-2).
983	0	do k=kmb+1,nz-1 ; if (Rcv_kmb(i) <= GV%Rlay(k)) exit ; enddo
984	0	kb(i) = k
985
986		! Backtrack, in case there are massive layers above that are stable
987		! in sigma-0.
988	0	do k=kb(i)-1,kmb+1,-1
989	0	if (rho_0(i,kmb) > rho_0(i,k)) exit
990	0	if (h(i,j,k)>2.0*GV%Angstrom_H) kb(i) = k
991		enddo
992		enddo
993
994	0	call set_density_ratios(h, tv, kb, G, GV, US, CS, j, ds_dsp1, rho_0)
995		else ! not bulkmixedlayer
996	0	kb_min = 2 ; kmb = 0
997	0	do i=is,ie ; kb(i) = 1 ; enddo
998	0	call set_density_ratios(h, tv, kb, G, GV, US, CS, j, ds_dsp1)
999		endif
1000
1001		! Determine maxEnt - the maximum permitted entrainment from below by each
1002		! interior layer.
1003	0	do k=2,nz-1 ; do i=is,ie
1004	0	dsp1_ds(i,k) = 1.0 / ds_dsp1(i,k)
1005		enddo ; enddo
1006	0	do i=is,ie ; dsp1_ds(i,nz) = 0.0 ; enddo
1007
1008	0	if (CS%bulkmixedlayer) then
1009	0	kmb = GV%nk_rho_varies
1010	0	do i=is,ie
1011	0	htot(i) = h(i,j,kmb)
1012	0	mFkb(i) = 0.0
1013	0	if (kb(i) < nz) mFkb(i) = ds_dsp1(i,kb(i)) * (h(i,j,kmb) - GV%Angstrom_H)
1014		enddo
1015	0	do k=1,kmb-1 ; do i=is,ie
1016	0	htot(i) = htot(i) + h(i,j,k)
1017	0	mFkb(i) = mFkb(i) + ds_dsp1(i,k+1)*(h(i,j,k) - GV%Angstrom_H)
1018		enddo ; enddo
1019		else
1020	0	do i=is,i
1021	0	maxEnt(i,1) = 0.0 ; htot(i) = h(i,j,1) - GV%Angstrom_H
1022		enddo
1023		endif
1024	0	do k=kb_min,nz-1 ; do i=is,ie
1025	0	if (k == kb(i)) then
1026	0	maxEnt(i,kb(i)) = mFkb(i)
1027	0	elseif (k > kb(i)) then
1028	0	if (CS%answer_date < 20190101) then
1029	0	maxEnt(i,k) = (1.0/dsp1_ds(i,k))*(maxEnt(i,k-1) + htot(i))
1030		else
1031	0	maxEnt(i,k) = ds_dsp1(i,k)*(maxEnt(i,k-1) + htot(i))
1032		endif
1033	0	htot(i) = htot(i) + (h(i,j,k) - GV%Angstrom_H)
1034		endif
1035		enddo ; enddo
1036
1037	0	do i=is,ie
1038	0	htot(i) = h(i,j,nz) - GV%Angstrom_H ; maxEnt(i,nz) = 0.0
1039	0	do_i(i) = (G%mask2dT(i,j) > 0.0)
1040		enddo
1041	0	do k=nz-1,kb_min,-1
1042	0	i_rem = 0
1043	0	do i=is,ie ; if (do_i(i)) then
1044	0	if (k<kb(i)) then ; do_i(i) = .false. ; cycle ; endif
1045	0	i_rem = i_rem + 1 ! Count the i-rows that are still being worked on.
1046	0	maxEnt(i,k) = MIN(maxEnt(i,k), dsp1_ds(i,k+1)*maxEnt(i,k+1) + htot(i))
1047	0	htot(i) = htot(i) + (h(i,j,k) - GV%Angstrom_H)
1048		endif ; enddo
1049	0	if (i_rem == 0) exit
1050		enddo ! k-loop
1051
1052		! Now set maxTKE and TKE_to_Kd.
1053	0	do i=is,ie
1054	0	maxTKE(i,1) = 0.0 ; TKE_to_Kd(i,1) = 0.0
1055	0	maxTKE(i,nz) = 0.0 ; TKE_to_Kd(i,nz) = 0.0
1056		enddo
1057	0	do k=2,kmb ; do i=is,ie
1058	0	maxTKE(i,k) = 0.0
1059	0	TKE_to_Kd(i,k) = 1.0 / ((N2_lay(i,k) + Omega2) * (dz(i,k) + dz_neglect))
1060		enddo ; enddo
1061	0	do k=kmb+1,kb_min-1 ; do i=is,ie
1062		! These are the properties in the deeper mixed and buffer layers, and
1063		! should perhaps be revisited.
1064	0	maxTKE(i,k) = 0.0 ; TKE_to_Kd(i,k) = 0.0
1065		enddo ; enddo
1066	0	do k=kb_min,nz-1 ; do i=is,ie
1067	0	if (k<kb(i)) then
1068	0	maxTKE(i,k) = 0.0
1069	0	TKE_to_Kd(i,k) = 0.0
1070		else
1071		! maxTKE is found by determining the kappa that gives maxEnt.
1072		! kappa_max = I_dt * dRho_int(i,K+1) * maxEnt(i,k) * &
1073		! G_IRho0(h(i,j,k) + dh_max) / (G_Rho0dRho_lay)
1074		! maxTKE(i,k) = GV%g_Earth_Z_T2 * dRho_lay * kappa_max
1075		! dRho_int should already be non-negative, so the max is redundant?
1076	0	dh_max = maxEnt(i,k) * (1.0 + dsp1_ds(i,k))
1077	0	dRho_lay = 0.5 * max(dRho_int(i,K) + dRho_int(i,K+1), 0.0)
1078
1079		! TKE_to_Kd should be rho_InSitu / G_Earth * (delta rho_InSitu)
1080		! The omega^2 term in TKE_to_Kd is due to a rescaling of the efficiency of turbulent
1081		! mixing by a factor of N^2 / (N^2 + Omega^2), as proposed by Melet et al., 2013?
1082	0	if (allocated(tv%SpV_avg)) then
1083		maxTKE(i,k) = I_dt * ((GV%H_to_RZ * grav * tv%SpV_avg(i,j,k)*2) &
1084		(0.5max(dRho_int(i,K+1) + dsp1_ds(i,k)dRho_int(i,K), 0.0))) * &
1085	0	((h(i,j,k) + dh_max) * maxEnt(i,k))
1086	0	TKE_to_Kd(i,k) = 1.0 / (grav * tv%SpV_avg(i,j,k) * dRho_lay + CS%omega*2 (dz(i,k) + dz_neglect))
1087		else
1088		maxTKE(i,k) = I_dt * (G_IRho0 * &
1089		(0.5max(dRho_int(i,K+1) + dsp1_ds(i,k)dRho_int(i,K), 0.0))) * &
1090	0	((h(i,j,k) + dh_max) * maxEnt(i,k))
1091	0	TKE_to_Kd(i,k) = 1.0 / (G_Rho0 * dRho_lay + CS%omega*2 (dz(i,k) + dz_neglect))
1092		endif
1093		endif
1094		enddo ; enddo
1095
1096		end subroutine find_TKE_to_Kd
1097
1098		!> Calculate Brunt-Vaisala frequency, N^2.
1099	720	subroutine find_N2(h, tv, T_f, S_f, fluxes, j, G, GV, US, CS, dRho_int, &
1100	1440	N2_lay, N2_int, N2_bot, Rho_bot, h_bot, k_bot)
1101		type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure
1102		type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure
1103		type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
1104		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
1105		intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2]
1106		type(thermo_var_ptrs), intent(in) :: tv !< Structure containing pointers to any available
1107		!! thermodynamic fields.
1108		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
1109		intent(in) :: T_f !< layer temperature with the values in massless layers
1110		!! filled vertically by diffusion [C ~> degC].
1111		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
1112		intent(in) :: S_f !< Layer salinities with values in massless
1113		!! layers filled vertically by diffusion [S ~> ppt].
1114		type(forcing), intent(in) :: fluxes !< A structure of thermodynamic surface fluxes
1115		integer, intent(in) :: j !< j-index of row to work on
1116		type(set_diffusivity_CS), pointer :: CS !< Diffusivity control structure
1117		real, dimension(SZI_(G),SZK_(GV)+1), &
1118		intent(out) :: dRho_int !< Change in locally referenced potential density
1119		!! across each interface [R ~> kg m-3].
1120		real, dimension(SZI_(G),SZK_(GV)+1), &
1121		intent(out) :: N2_int !< The squared buoyancy frequency at the interfaces [T-2 ~> s-2].
1122		real, dimension(SZI_(G),SZK_(GV)), &
1123		intent(out) :: N2_lay !< The squared buoyancy frequency of the layers [T-2 ~> s-2].
1124		real, dimension(SZI_(G)), intent(out) :: N2_bot !< The near-bottom squared buoyancy frequency [T-2 ~> s-2].
1125		real, dimension(SZI_(G)), intent(out) :: Rho_bot !< Near-bottom density [R ~> kg m-3].
1126		real, dimension(SZI_(G)), optional, intent(out) :: h_bot !< Bottom boundary layer thickness [H ~> m or kg m-2].
1127		integer, dimension(SZI_(G)), optional, intent(out) :: k_bot !< Bottom boundary layer top layer index.
1128
1129		! Local variables
1130		real, dimension(SZI_(G),SZK_(GV)+1) :: &
1131	1440	pres, & ! pressure at each interface [R L2 T-2 ~> Pa]
1132	1440	dRho_int_unfilt, & ! unfiltered density differences across interfaces [R ~> kg m-3]
1133	1440	dRho_dT, & ! partial derivative of density wrt temp [R C-1 ~> kg m-3 degC-1]
1134	1440	dRho_dS ! partial derivative of density wrt saln [R S-1 ~> kg m-3 ppt-1]
1135		real, dimension(SZI_(G),SZK_(GV)) :: &
1136	1440	dz ! Height change across layers [Z ~> m]
1137		real, dimension(SZI_(G)) :: &
1138	1440	Temp_int, & ! temperature at each interface [C ~> degC]
1139	1440	Salin_int, & ! salinity at each interface [S ~> ppt]
1140	1440	drho_bot, & ! A density difference [R ~> kg m-3]
1141	1440	h_amp, & ! The topographic roughness amplitude [Z ~> m].
1142	1440	dz_BBL_avg, & ! The distance over which to average to find the near-bottom density [Z ~> m]
1143	1440	hb, & ! The thickness of the bottom layer [H ~> m or kg m-2]
1144	1440	z_from_bot ! The height above the bottom [Z ~> m]
1145
1146		real :: dz_int ! Vertical distance associated with an interface [Z ~> m]
1147		real :: G_Rho0 ! Gravitational acceleration, perhaps divided by Boussinesq reference density,
1148		! times some unit conversion factors [H T-2 R-1 ~> m4 s-2 kg-1 or m s-2].
1149		real :: H_neglect ! A negligibly small thickness [H ~> m or kg m-2]
1150
1151	1440	logical :: do_i(SZI_(G)), do_any
1152		integer, dimension(2) :: EOSdom ! The i-computational domain for the equation of state
1153		integer :: i, k, is, ie, nz
1154
1155	720	is = G%isc ; ie = G%iec ; nz = GV%ke
1156	720	G_Rho0 = GV%g_Earth_Z_T2 / GV%H_to_RZ
1157	720	H_neglect = GV%H_subroundoff
1158
1159		! Find the (limited) density jump across each interface.
1160	87120	do i=is,ie
1161	86400	dRho_int(i,1) = 0.0 ; dRho_int(i,nz+1) = 0.0
1162	87120	dRho_int_unfilt(i,1) = 0.0 ; dRho_int_unfilt(i,nz+1) = 0.0
1163		enddo
1164	720	if (associated(tv%eqn_of_state)) then
1165	720	if (associated(fluxes%p_surf)) then
1166	87120	do i=is,ie ; pres(i,1) = fluxes%p_surf(i,j) ; enddo
1167		else
1168	0	do i=is,ie ; pres(i,1) = 0.0 ; enddo
1169		endif
1170	720	EOSdom(:) = EOS_domain(G%HI)
1171	54000	do K=2,nz
1172	6446880	do i=is,ie
1173	6393600	pres(i,K) = pres(i,K-1) + (GV%g_EarthGV%H_to_RZ)h(i,j,k-1)
1174	6393600	Temp_Int(i) = 0.5 * (T_f(i,j,k) + T_f(i,j,k-1))
1175	6446880	Salin_Int(i) = 0.5 * (S_f(i,j,k) + S_f(i,j,k-1))
1176		enddo
1177		call calculate_density_derivs(Temp_int, Salin_int, pres(:,K), dRho_dT(:,K), dRho_dS(:,K), &
1178	53280	tv%eqn_of_state, EOSdom)
1179	6447600	do i=is,ie
1180		dRho_int(i,K) = max(dRho_dT(i,K)*(T_f(i,j,k) - T_f(i,j,k-1)) + &
1181	6393600	dRho_dS(i,K)*(S_f(i,j,k) - S_f(i,j,k-1)), 0.0)
1182		dRho_int_unfilt(i,K) = max(dRho_dT(i,K)*(tv%T(i,j,k) - tv%T(i,j,k-1)) + &
1183	6446880	dRho_dS(i,K)*(tv%S(i,j,k) - tv%S(i,j,k-1)), 0.0)
1184		enddo
1185		enddo
1186		else
1187	0	do K=2,nz ; do i=is,ie
1188	0	dRho_int(i,K) = GV%Rlay(k) - GV%Rlay(k-1)
1189		enddo ; enddo
1190		endif
1191
1192		! Find the vertical distances across layers.
1193	720	call thickness_to_dz(h, tv, dz, j, G, GV)
1194
1195		! Set the buoyancy frequencies.
1196	6534720	do k=1,nz ; do i=is,ie
1197		N2_lay(i,k) = G_Rho0 * 0.5*(dRho_int(i,K) + dRho_int(i,K+1)) / &
1198	6534000	(h(i,j,k) + H_neglect)
1199		enddo ; enddo
1200	87120	do i=is,ie ; N2_int(i,1) = 0.0 ; N2_int(i,nz+1) = 0.0 ; enddo
1201	6447600	do K=2,nz ; do i=is,ie
1202		N2_int(i,K) = G_Rho0 * dRho_int(i,K) / &
1203	6446880	(0.5*(h(i,j,k-1) + h(i,j,k) + H_neglect))
1204		enddo ; enddo
1205
1206		! Find the bottom boundary layer stratification, and use this in the deepest layers.
1207	87120	do i=is,ie
1208	86400	hb(i) = 0.0 ; dRho_bot(i) = 0.0 ; h_amp(i) = 0.0
1209	86400	z_from_bot(i) = 0.5*dz(i,nz)
1210	87120	do_i(i) = (G%mask2dT(i,j) > 0.0)
1211		enddo
1212	720	if (CS%use_tidal_mixing) call tidal_mixing_h_amp(h_amp, G, j, CS%tidal_mixing)
1213
1214	720	do k=nz,2,-1
1215	720	do_any = .false.
1216	87120	do i=is,ie ; if (do_i(i)) then
1217	60168	dz_int = 0.5*(dz(i,k) + dz(i,k-1))
1218	60168	z_from_bot(i) = z_from_bot(i) + dz_int ! middle of the layer above
1219
1220	60168	hb(i) = hb(i) + 0.5*(h(i,j,k) + h(i,j,k-1))
1221	60168	drho_bot(i) = drho_bot(i) + dRho_int(i,K)
1222
1223	60168	if (z_from_bot(i) > h_amp(i)) then
1224	60168	if (k>2) then
1225		! Always include at least one full layer.
1226	60168	hb(i) = hb(i) + 0.5*(h(i,j,k-1) + h(i,j,k-2))
1227	60168	drho_bot(i) = drho_bot(i) + dRho_int(i,K-1)
1228		endif
1229	60168	do_i(i) = .false.
1230		else
1231	0	do_any = .true.
1232		endif
1233		endif ; enddo
1234	720	if (.not.do_any) exit
1235		enddo
1236
1237	87120	do i=is,ie
1238	86400	if (hb(i) > 0.0) then
1239	60168	N2_bot(i) = (G_Rho0 * drho_bot(i)) / hb(i)
1240	26232	else ; N2_bot(i) = 0.0 ; endif
1241	86400	z_from_bot(i) = 0.5*dz(i,nz)
1242	87120	do_i(i) = (G%mask2dT(i,j) > 0.0)
1243		enddo
1244
1245	720	do k=nz,2,-1
1246	720	do_any = .false.
1247	87120	do i=is,ie ; if (do_i(i)) then
1248	60168	dz_int = 0.5*(dz(i,k) + dz(i,k-1))
1249	60168	z_from_bot(i) = z_from_bot(i) + dz_int ! middle of the layer above
1250
1251	60168	N2_int(i,K) = N2_bot(i)
1252	60168	if (k>2) N2_lay(i,k-1) = N2_bot(i)
1253
1254	60168	if (z_from_bot(i) > h_amp(i)) then
1255	60168	if (k>2) N2_int(i,K-1) = N2_bot(i)
1256	60168	do_i(i) = .false.
1257		else
1258	0	do_any = .true.
1259		endif
1260		endif ; enddo
1261	720	if (.not.do_any) exit
1262		enddo
1263
1264	720	if (associated(tv%eqn_of_state)) then
1265	6621840	do K=1,nz+1 ; do i=is,ie
1266	6621120	dRho_int(i,K) = dRho_int_unfilt(i,K)
1267		enddo ; enddo
1268		endif
1269
1270		! Average over the larger of the envelope of the topography or a minimal distance.
1271	87120	do i=is,ie ; dz_BBL_avg(i) = max(h_amp(i), CS%dz_BBL_avg_min) ; enddo
1272	720	call find_rho_bottom(G, GV, US, tv, h, dz, pres, dz_BBL_avg, j, Rho_bot, h_bot, k_bot)
1273
1274	1440	end subroutine find_N2
1275
1276		!> This subroutine sets the additional diffusivities of temperature and
1277		!! salinity due to double diffusion, using the same functional form as is
1278		!! used in MOM4.1, and taken from the appendix of Danabasoglu et al. (2006), which updates
1279		!! what was in Large et al. (1994). All the coefficients here should probably
1280		!! be made run-time variables rather than hard-coded constants.
1281		!!
1282		!! \todo Find reference for NCAR tech note above.
1283	0	subroutine double_diffusion(tv, h, T_f, S_f, j, G, GV, US, CS, Kd_T_dd, Kd_S_dd)
1284		type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure.
1285		type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure.
1286		type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
1287		type(thermo_var_ptrs), intent(in) :: tv !< Structure containing pointers to any available
1288		!! thermodynamic fields; absent fields have NULL ptrs.
1289		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
1290		intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2].
1291		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
1292		intent(in) :: T_f !< layer temperatures with the values in massless layers
1293		!! filled vertically by diffusion [C ~> degC].
1294		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
1295		intent(in) :: S_f !< Layer salinities with values in massless
1296		!! layers filled vertically by diffusion [S ~> ppt].
1297		integer, intent(in) :: j !< Meridional index upon which to work.
1298		type(set_diffusivity_CS), pointer :: CS !< Module control structure.
1299		real, dimension(SZI_(G),SZK_(GV)+1), &
1300		intent(out) :: Kd_T_dd !< Interface double diffusion diapycnal
1301		!! diffusivity for temp [H Z T-1 ~> m2 s-1 or kg m-1 s-1]
1302		real, dimension(SZI_(G),SZK_(GV)+1), &
1303		intent(out) :: Kd_S_dd !< Interface double diffusion diapycnal
1304		!! diffusivity for saln [H Z T-1 ~> m2 s-1 or kg m-1 s-1]
1305
1306		real, dimension(SZI_(G)) :: &
1307	0	dRho_dT, & ! partial derivatives of density with respect to temperature [R C-1 ~> kg m-3 degC-1]
1308	0	dRho_dS, & ! partial derivatives of density with respect to salinity [R S-1 ~> kg m-3 ppt-1]
1309	0	pres, & ! pressure at each interface [R L2 T-2 ~> Pa]
1310	0	Temp_int, & ! temperature at interfaces [C ~> degC]
1311	0	Salin_int ! Salinity at interfaces [S ~> ppt]
1312
1313		real :: alpha_dT ! density difference between layers due to temp diffs [R ~> kg m-3]
1314		real :: beta_dS ! density difference between layers due to saln diffs [R ~> kg m-3]
1315
1316		real :: Rrho ! vertical density ratio [nondim]
1317		real :: diff_dd ! factor for double-diffusion [nondim]
1318		real :: Kd_dd ! The dominant double diffusive diffusivity [H Z T-1 ~> m2 s-1 or kg m-1 s-1]
1319		real :: prandtl ! flux ratio for diffusive convection regime [nondim]
1320
1321		integer, dimension(2) :: EOSdom ! The i-computational domain for the equation of state
1322		integer :: i, k, is, ie, nz
1323	0	is = G%isc ; ie = G%iec ; nz = GV%ke
1324
1325	0	if (associated(tv%eqn_of_state)) then
1326	0	do i=is,ie
1327	0	pres(i) = 0.0 ; Kd_T_dd(i,1) = 0.0 ; Kd_S_dd(i,1) = 0.0
1328	0	Kd_T_dd(i,nz+1) = 0.0 ; Kd_S_dd(i,nz+1) = 0.0
1329		enddo
1330	0	if (associated(tv%p_surf)) then ; do i=is,ie ; pres(i) = tv%p_surf(i,j) ; enddo ; endif
1331	0	EOSdom(:) = EOS_domain(G%HI)
1332	0	do K=2,nz
1333	0	do i=is,ie
1334	0	pres(i) = pres(i) + (GV%g_EarthGV%H_to_RZ)h(i,j,k-1)
1335	0	Temp_Int(i) = 0.5 * (T_f(i,j,k-1) + T_f(i,j,k))
1336	0	Salin_Int(i) = 0.5 * (S_f(i,j,k-1) + S_f(i,j,k))
1337		enddo
1338		call calculate_density_derivs(Temp_int, Salin_int, pres, dRho_dT, dRho_dS, &
1339	0	tv%eqn_of_state, EOSdom)
1340
1341	0	do i=is,ie
1342	0	alpha_dT = -1.0dRho_dT(i) (T_f(i,j,k-1) - T_f(i,j,k))
1343	0	beta_dS = dRho_dS(i) * (S_f(i,j,k-1) - S_f(i,j,k))
1344
1345	0	if ((alpha_dT > beta_dS) .and. (beta_dS > 0.0)) then ! salt finger case
1346	0	Rrho = min(alpha_dT / beta_dS, CS%Max_Rrho_salt_fingers)
1347	0	diff_dd = 1.0 - ((RRho-1.0)/(CS%Max_Rrho_salt_fingers-1.0))
1348	0	Kd_dd = CS%Max_salt_diff_salt_fingers * diff_dddiff_dddiff_dd
1349	0	Kd_T_dd(i,K) = 0.7 * Kd_dd
1350	0	Kd_S_dd(i,K) = Kd_dd
1351	0	elseif ((alpha_dT < 0.) .and. (beta_dS < 0.) .and. (alpha_dT > beta_dS)) then ! diffusive convection
1352	0	Rrho = alpha_dT / beta_dS
1353	0	Kd_dd = CS%Kv_molecular * 0.909 * exp(4.6 * exp(-0.54 * (1/Rrho - 1)))
1354	0	prandtl = 0.15*Rrho
1355	0	if (Rrho > 0.5) prandtl = (1.85-0.85/Rrho)*Rrho
1356	0	Kd_T_dd(i,K) = Kd_dd
1357	0	Kd_S_dd(i,K) = prandtl * Kd_dd
1358		else
1359	0	Kd_T_dd(i,K) = 0.0 ; Kd_S_dd(i,K) = 0.0
1360		endif
1361		enddo
1362		enddo
1363		endif
1364
1365	0	end subroutine double_diffusion
1366
1367		!> This routine adds diffusion sustained by flow energy extracted by bottom drag.
1368	0	subroutine add_drag_diffusivity(h, u, v, tv, fluxes, visc, j, TKE_to_Kd, maxTKE, &
1369	0	kb, rho_bot, G, GV, US, CS, Kd_lay, Kd_int, Kd_BBL)
1370		type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure
1371		type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure
1372		type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
1373		real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &
1374		intent(in) :: u !< The zonal velocity [L T-1 ~> m s-1]
1375		real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &
1376		intent(in) :: v !< The meridional velocity [L T-1 ~> m s-1]
1377		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
1378		intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2]
1379		type(thermo_var_ptrs), intent(in) :: tv !< Structure containing pointers to any available
1380		!! thermodynamic fields.
1381		type(forcing), intent(in) :: fluxes !< A structure of thermodynamic surface fluxes
1382		type(vertvisc_type), intent(in) :: visc !< Structure containing vertical viscosities, bottom
1383		!! boundary layer properties and related fields
1384		integer, intent(in) :: j !< j-index of row to work on
1385		real, dimension(SZI_(G),SZK_(GV)), intent(in) :: TKE_to_Kd !< The conversion rate between the TKE
1386		!! TKE dissipated within a layer and the
1387		!! diapycnal diffusivity within that layer,
1388		!! usually (~Rho_0 / (G_Earth * dRho_lay))
1389		!! [T2 Z-1 ~> s2 m-1]
1390		real, dimension(SZI_(G),SZK_(GV)), intent(in) :: maxTKE !< The energy required to for a layer to entrain to its
1391		!! maximum-realizable thickness [H Z2 T-3 ~> m3 s-3 or W m-2]
1392		integer, dimension(SZI_(G)), intent(in) :: kb !< Index of lightest layer denser than the buffer
1393		!! layer, or -1 without a bulk mixed layer
1394		real, dimension(SZI_(G)), intent(in) :: rho_bot !< In situ density averaged over a near-bottom
1395		!! region [R ~> kg m-3]
1396		type(set_diffusivity_CS), pointer :: CS !< Diffusivity control structure
1397		real, dimension(SZI_(G),SZK_(GV)), intent(inout) :: Kd_lay !< The diapycnal diffusivity in layers,
1398		!! [H Z T-1 ~> m2 s-1 or kg m-1 s-1]
1399		real, dimension(SZI_(G),SZK_(GV)+1), intent(inout) :: Kd_int !< The diapycnal diffusivity at interfaces,
1400		!! [H Z T-1 ~> m2 s-1 or kg m-1 s-1]
1401		real, dimension(:,:,:), pointer :: Kd_BBL !< Interface BBL diffusivity
1402		!! [H Z T-1 ~> m2 s-1 or kg m-1 s-1]
1403
1404		! This routine adds diffusion sustained by flow energy extracted by bottom drag.
1405
1406		real, dimension(SZK_(GV)+1) :: &
1407	0	Rint ! coordinate density of an interface [R ~> kg m-3]
1408		real, dimension(SZI_(G)) :: &
1409	0	htot, & ! total thickness above or below a layer, or the
1410		! integrated thickness in the BBL [H ~> m or kg m-2].
1411	0	rho_htot, & ! running integral with depth of density [R H ~> kg m-2 or kg2 m-5]
1412	0	gh_sum_top, & ! BBL value of g'h that can be supported by
1413		! the local ustar, times R0_g [R H ~> kg m-2 or kg2 m-5]
1414	0	Rho_top, & ! density at top of the BBL [R ~> kg m-3]
1415	0	TKE, & ! turbulent kinetic energy available to drive
1416		! bottom-boundary layer mixing in a layer [H Z2 T-3 ~> m3 s-3 or W m-2]
1417	0	I2decay ! inverse of twice the TKE decay scale [H-1 ~> m-1 or m2 kg-1].
1418
1419		real :: TKE_to_layer ! TKE used to drive mixing in a layer [H Z2 T-3 ~> m3 s-3 or W m-2]
1420		real :: TKE_Ray ! TKE from layer Rayleigh drag used to drive mixing in layer [H Z2 T-3 ~> m3 s-3 or W m-2]
1421		real :: TKE_here ! TKE that goes into mixing in this layer [H Z2 T-3 ~> m3 s-3 or W m-2]
1422		real :: dRl, dRbot ! temporaries holding density differences [R ~> kg m-3]
1423		real :: cdrag_sqrt ! square root of the drag coefficient [nondim]
1424		real :: ustar_h ! Ustar at a thickness point rescaled into thickness
1425		! flux units [H T-1 ~> m s-1 or kg m-2 s-1].
1426		real :: absf ! average absolute Coriolis parameter around a thickness point [T-1 ~> s-1]
1427		real :: R0_g ! Rho0 / G_Earth [R T2 H-1 ~> kg s2 m-4 or s2 m-1]
1428		real :: delta_Kd ! increment to Kd from the bottom boundary layer mixing [H Z T-1 ~> m2 s-1 or kg m-1 s-1]
1429		logical :: Rayleigh_drag ! Set to true if Rayleigh drag velocities
1430		! defined in visc, on the assumption that this
1431		! extracted energy also drives diapycnal mixing.
1432
1433	0	logical :: domore, do_i(SZI_(G))
1434		logical :: do_diag_Kd_BBL
1435
1436		integer :: i, k, is, ie, nz, i_rem, kb_min
1437	0	is = G%isc ; ie = G%iec ; nz = GV%ke
1438
1439	0	do_diag_Kd_BBL = associated(Kd_BBL)
1440
1441	0	if (.not.(CS%bottomdraglaw .and. (CS%BBL_effic > 0.0))) return
1442
1443	0	cdrag_sqrt = sqrt(CS%cdrag)
1444	0	TKE_Ray = 0.0 ; Rayleigh_drag = .false.
1445	0	if (allocated(visc%Ray_u) .and. allocated(visc%Ray_v)) Rayleigh_drag = .true.
1446
1447	0	R0_g = GV%H_to_RZ / GV%g_Earth_Z_T2
1448
1449	0	do K=2,nz ; Rint(K) = 0.5*(GV%Rlay(k-1)+GV%Rlay(k)) ; enddo
1450
1451	0	kb_min = max(GV%nk_rho_varies+1,2)
1452
1453		! The turbulence decay scale is 0.5*ustar/f from K&E & MOM_vertvisc.F90
1454		! Any turbulence that makes it into the mixed layers is assumed
1455		! to be relatively small and is discarded.
1456	0	do i=is,ie
1457	0	ustar_h = visc%ustar_BBL(i,j)
1458	0	if (associated(fluxes%ustar_tidal)) then
1459	0	if (allocated(tv%SpV_avg)) then
1460	0	ustar_h = ustar_h + GV%RZ_to_Hrho_bot(i) fluxes%ustar_tidal(i,j)
1461		else
1462	0	ustar_h = ustar_h + GV%Z_to_H * fluxes%ustar_tidal(i,j)
1463		endif
1464		endif
1465		absf = 0.25 * ((abs(G%CoriolisBu(I-1,J-1)) + abs(G%CoriolisBu(I,J))) + &
1466	0	(abs(G%CoriolisBu(I-1,J)) + abs(G%CoriolisBu(I,J-1))))
1467	0	if ((ustar_h > 0.0) .and. (absf > 0.5CS%IMax_decayustar_h)) then
1468	0	I2decay(i) = absf / ustar_h
1469		else
1470		! The maximum decay scale should be something of order 200 m.
1471		! If ustar_h = 0, this is land so this value doesn't matter.
1472	0	I2decay(i) = 0.5*CS%IMax_decay
1473		endif
1474	0	if (CS%drag_diff_answer_date <= 20250301) then
1475	0	TKE(i) = ((CS%BBL_effic * cdrag_sqrt) * exp(-I2decay(i)h(i,j,nz)) ) visc%BBL_meanKE_loss_sqrtCd(i,j)
1476		else
1477	0	TKE(i) = (CS%BBL_effic * exp(-I2decay(i)h(i,j,nz)) ) visc%BBL_meanKE_loss(i,j)
1478		endif
1479
1480	0	if (associated(fluxes%BBL_tidal_dis)) &
1481		TKE(i) = TKE(i) + fluxes%BBL_tidal_dis(i,j) * GV%RZ_to_H * &
1482	0	(CS%BBL_effic * exp(-I2decay(i)*h(i,j,nz)))
1483
1484		! Distribute the work over a BBL of depth 20^2 ustar^2 / g' following
1485		! Killworth & Edwards (1999) and Zilitikevich & Mironov (1996).
1486		! Rho_top is determined by finding the density where
1487		! integral(bottom, Z) (rho(z') - rho(Z)) dz' = rho_0 400 ustar^2 / g
1488
1489	0	gh_sum_top(i) = R0_g * 400.0 * ustar_h**2
1490
1491	0	do_i(i) = (G%mask2dT(i,j) > 0.0)
1492	0	htot(i) = h(i,j,nz)
1493	0	rho_htot(i) = GV%Rlay(nz)*(h(i,j,nz))
1494	0	Rho_top(i) = GV%Rlay(1)
1495	0	if (CS%bulkmixedlayer .and. do_i(i)) Rho_top(i) = GV%Rlay(kb(i)-1)
1496		enddo
1497
1498	0	do k=nz-1,2,-1 ; domore = .false.
1499	0	do i=is,ie ; if (do_i(i)) then
1500	0	htot(i) = htot(i) + h(i,j,k)
1501	0	rho_htot(i) = rho_htot(i) + GV%Rlay(k)*(h(i,j,k))
1502	0	if (htot(i)*GV%Rlay(k-1) <= (rho_htot(i) - gh_sum_top(i))) then
1503		! The top of the mixing is in the interface atop the current layer.
1504	0	Rho_top(i) = (rho_htot(i) - gh_sum_top(i)) / htot(i)
1505	0	do_i(i) = .false.
1506	0	elseif (k <= kb(i)) then ; do_i(i) = .false.
1507	0	else ; domore = .true. ; endif
1508		endif ; enddo
1509	0	if (.not.domore) exit
1510		enddo ! k-loop
1511
1512	0	do i=is,ie ; do_i(i) = (G%mask2dT(i,j) > 0.0) ; enddo
1513	0	do k=nz-1,kb_min,-1
1514	0	i_rem = 0
1515	0	do i=is,ie ; if (do_i(i)) then
1516	0	if (k<kb(i)) then ; do_i(i) = .false. ; cycle ; endif
1517	0	i_rem = i_rem + 1 ! Count the i-rows that are still being worked on.
1518		! Apply vertical decay of the turbulent energy. This energy is
1519		! simply lost.
1520	0	TKE(i) = TKE(i) * exp(-I2decay(i) * (h(i,j,k) + h(i,j,k+1)))
1521
1522	0	if (maxTKE(i,k) <= 0.0) cycle
1523
1524		! This is an analytic integral where diffusivity is a quadratic function of
1525		! rho that goes asymptotically to 0 at Rho_top (vaguely following KPP?).
1526	0	if (TKE(i) > 0.0) then
1527	0	if (Rint(K) <= Rho_top(i)) then
1528	0	TKE_to_layer = TKE(i)
1529		else
1530	0	dRl = Rint(K+1) - Rint(K) ; dRbot = Rint(K+1) - Rho_top(i)
1531		TKE_to_layer = TKE(i) * dRl * &
1532	0	(3.0dRbot(Rint(K) - Rho_top(i)) + dRl2) / (dRbot3)
1533		endif
1534	0	else ; TKE_to_layer = 0.0 ; endif
1535
1536		! TKE_Ray has been initialized to 0 above.
1537	0	if (Rayleigh_drag) TKE_Ray = 0.5CS%BBL_effic G%IareaT(i,j) * &
1538		(((G%areaCu(I-1,j) * visc%Ray_u(I-1,j,k) * u(I-1,j,k)**2) + &
1539		(G%areaCu(I,j) * visc%Ray_u(I,j,k) * u(I,j,k)**2)) + &
1540		((G%areaCv(i,J-1) * visc%Ray_v(i,J-1,k) * v(i,J-1,k)**2) + &
1541	0	(G%areaCv(i,J) * visc%Ray_v(i,J,k) * v(i,J,k)**2)))
1542
1543	0	if (TKE_to_layer + TKE_Ray > 0.0) then
1544	0	if (CS%BBL_mixing_as_max) then
1545	0	if (TKE_to_layer + TKE_Ray > maxTKE(i,k)) &
1546	0	TKE_to_layer = maxTKE(i,k) - TKE_Ray
1547
1548	0	TKE(i) = TKE(i) - TKE_to_layer
1549
1550	0	if (Kd_lay(i,k) < (TKE_to_layer + TKE_Ray) * TKE_to_Kd(i,k)) then
1551	0	delta_Kd = (TKE_to_layer + TKE_Ray) * TKE_to_Kd(i,k) - Kd_lay(i,k)
1552	0	if ((CS%Kd_max >= 0.0) .and. (delta_Kd > CS%Kd_max)) then
1553	0	delta_Kd = CS%Kd_max
1554	0	Kd_lay(i,k) = Kd_lay(i,k) + delta_Kd
1555		else
1556	0	Kd_lay(i,k) = (TKE_to_layer + TKE_Ray) * TKE_to_Kd(i,k)
1557		endif
1558	0	Kd_int(i,K) = Kd_int(i,K) + 0.5 * delta_Kd
1559	0	Kd_int(i,K+1) = Kd_int(i,K+1) + 0.5 * delta_Kd
1560	0	if (do_diag_Kd_BBL) then
1561	0	Kd_BBL(i,j,K) = Kd_BBL(i,j,K) + 0.5 * delta_Kd
1562	0	Kd_BBL(i,j,K+1) = Kd_BBL(i,j,K+1) + 0.5 * delta_Kd
1563		endif
1564		endif
1565		else
1566	0	if (Kd_lay(i,k) >= maxTKE(i,k) * TKE_to_Kd(i,k)) then
1567	0	TKE_here = 0.0
1568	0	TKE(i) = TKE(i) + TKE_Ray
1569	0	elseif (Kd_lay(i,k) + (TKE_to_layer + TKE_Ray) * TKE_to_Kd(i,k) > &
1570		maxTKE(i,k) * TKE_to_Kd(i,k)) then
1571	0	TKE_here = ((TKE_to_layer + TKE_Ray) + Kd_lay(i,k) / TKE_to_Kd(i,k)) - maxTKE(i,k)
1572	0	TKE(i) = (TKE(i) - TKE_here) + TKE_Ray
1573		else
1574	0	TKE_here = TKE_to_layer + TKE_Ray
1575	0	TKE(i) = TKE(i) - TKE_to_layer
1576		endif
1577	0	if (TKE(i) < 0.0) TKE(i) = 0.0 ! This should be unnecessary?
1578
1579	0	if (TKE_here > 0.0) then
1580	0	delta_Kd = TKE_here * TKE_to_Kd(i,k)
1581	0	if (CS%Kd_max >= 0.0) delta_Kd = min(delta_Kd, CS%Kd_max)
1582	0	Kd_lay(i,k) = Kd_lay(i,k) + delta_Kd
1583	0	Kd_int(i,K) = Kd_int(i,K) + 0.5 * delta_Kd
1584	0	Kd_int(i,K+1) = Kd_int(i,K+1) + 0.5 * delta_Kd
1585	0	if (do_diag_Kd_BBL) then
1586	0	Kd_BBL(i,j,K) = Kd_BBL(i,j,K) + 0.5 * delta_Kd
1587	0	Kd_BBL(i,j,K+1) = Kd_BBL(i,j,K+1) + 0.5 * delta_Kd
1588		endif
1589		endif
1590		endif
1591		endif
1592
1593		! This may be risky - in the case that there are exactly zero
1594		! velocities at 4 neighboring points, but nonzero velocities
1595		! above the iterations would stop too soon. I don't see how this
1596		! could happen in practice. RWH
1597	0	if ((TKE(i)<= 0.0) .and. (TKE_Ray == 0.0)) then
1598	0	do_i(i) = .false. ; i_rem = i_rem - 1
1599		endif
1600
1601		endif ; enddo
1602	0	if (i_rem == 0) exit
1603		enddo ! k-loop
1604
1605		end subroutine add_drag_diffusivity
1606
1607		!> Calculates a BBL diffusivity use a Prandtl number 1 diffusivity with a law of the
1608		!! wall turbulent viscosity, up to a BBL height where the energy used for mixing has
1609		!! consumed the mechanical TKE input.
1610	720	subroutine add_LOTW_BBL_diffusivity(h, u, v, tv, fluxes, visc, j, N2_int, Rho_bot, Kd_int, &
1611	720	G, GV, US, CS, Kd_BBL, Kd_lay)
1612		type(ocean_grid_type), intent(in) :: G !< Grid structure
1613		type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure
1614		type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
1615		real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &
1616		intent(in) :: u !< u component of flow [L T-1 ~> m s-1]
1617		real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &
1618		intent(in) :: v !< v component of flow [L T-1 ~> m s-1]
1619		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
1620		intent(in) :: h !< Layer thickness [H ~> m or kg m-2]
1621		type(thermo_var_ptrs), intent(in) :: tv !< Structure containing pointers to any available
1622		!! thermodynamic fields.
1623		type(forcing), intent(in) :: fluxes !< Surface fluxes structure
1624		type(vertvisc_type), intent(in) :: visc !< Structure containing vertical viscosities, bottom
1625		!! boundary layer properties and related fields.
1626		integer, intent(in) :: j !< j-index of row to work on
1627		real, dimension(SZI_(G),SZK_(GV)+1), &
1628		intent(in) :: N2_int !< Square of Brunt-Vaisala at interfaces [T-2 ~> s-2]
1629		real, dimension(SZI_(G)), intent(in) :: rho_bot !< In situ density averaged over a near-bottom
1630		!! region [R ~> kg m-3]
1631		real, dimension(SZI_(G),SZK_(GV)+1), &
1632		intent(inout) :: Kd_int !< Interface net diffusivity [H Z T-1 ~> m2 s-1 or kg m-1 s-1]
1633		type(set_diffusivity_CS), pointer :: CS !< Diffusivity control structure
1634		real, dimension(:,:,:), pointer :: Kd_BBL !< Interface BBL diffusivity [H Z T-1 ~> m2 s-1 or kg m-1 s-1]
1635		real, dimension(SZI_(G),SZK_(GV)), &
1636		optional, intent(inout) :: Kd_lay !< Layer net diffusivity [H Z T-1 ~> m2 s-1 or kg m-1 s-1]
1637
1638		! Local variables
1639	720	real :: dz(SZI_(G),SZK_(GV)) ! Height change across layers [Z ~> m]
1640	720	real :: dz_above(SZK_(GV)+1) ! Distance from each interface to the surface [Z ~> m]
1641		real :: TKE_column ! net TKE input into the column [H Z2 T-3 ~> m3 s-3 or W m-2]
1642		real :: BBL_meanKE_dis ! Sum of tidal and mean kinetic energy dissipation in the bottom boundary layer, which
1643		! can act as a source of TKE [H L2 T-3 ~> m3 s-3 or W m-2]
1644		real :: TKE_remaining ! remaining TKE available for mixing in this layer and above [H Z2 T-3 ~> m3 s-3 or W m-2]
1645		real :: TKE_consumed ! TKE used for mixing in this layer [H Z2 T-3 ~> m3 s-3 or W m-2]
1646		real :: TKE_Kd_wall ! TKE associated with unlimited law of the wall mixing [H Z2 T-3 ~> m3 s-3 or W m-2]
1647		real :: cdrag_sqrt ! square root of the drag coefficient [nondim]
1648		real :: ustar ! value of ustar at a thickness point [H T-1 ~> m s-1 or kg m-2 s-1].
1649		real :: ustar2 ! The square of ustar [H2 T-2 ~> m2 s-2 or kg2 m-4 s-2]
1650		real :: absf ! average absolute value of Coriolis parameter around a thickness point [T-1 ~> s-1]
1651		real :: dz_int ! Distance between the center of the layers around an interface [Z ~> m]
1652		real :: z_bot ! Distance to interface K from bottom [Z ~> m]
1653		real :: h_bot ! Total thickness between interface K and the bottom [H ~> m or kg m-2]
1654		real :: D_minus_z ! Distance between interface k and the surface [Z ~> m]
1655		real :: total_depth ! Total distance between the seafloor and the sea surface [Z ~> m]
1656		real :: Idecay ! Inverse of decay scale used for "Joule heating" loss of TKE with
1657		! height [H-1 ~> m-1 or m2 kg-1].
1658		real :: Kd_wall ! Law of the wall diffusivity [H Z T-1 ~> m2 s-1 or kg m-1 s-1]
1659		real :: Kd_lower ! diffusivity for lower interface [H Z T-1 ~> m2 s-1 or kg m-1 s-1]
1660		real :: ustar_D ! The extent of the water column times u* [H Z T-1 ~> m2 s-1 or Pa s].
1661		real :: N2_min ! Minimum value of N2 to use in calculation of TKE_Kd_wall [T-2 ~> s-2]
1662		logical :: Rayleigh_drag ! Set to true if there are Rayleigh drag velocities defined in visc, on
1663		! the assumption that this extracted energy also drives diapycnal mixing.
1664		integer :: i, k
1665		logical :: do_diag_Kd_BBL
1666
1667	720	if (.not.(CS%bottomdraglaw .and. (CS%BBL_effic > 0.0))) return
1668	720	do_diag_Kd_BBL = associated(Kd_BBL)
1669
1670	720	N2_min = 0.
1671	720	if (CS%LOTW_BBL_use_omega) N2_min = CS%omega**2
1672
1673		! Determine whether to add Rayleigh drag contribution to TKE
1674	720	Rayleigh_drag = .false.
1675	720	if (allocated(visc%Ray_u) .and. allocated(visc%Ray_v)) Rayleigh_drag = .true.
1676	720	cdrag_sqrt = sqrt(CS%cdrag)
1677
1678		! Find the vertical distances across layers.
1679	720	call thickness_to_dz(h, tv, dz, j, G, GV)
1680
1681	87120	do i=G%isc,G%iec ! Developed in single-column mode
1682
1683		! Column-wise parameters.
1684		absf = 0.25 * ((abs(G%CoriolisBu(I-1,J-1)) + abs(G%CoriolisBu(I,J))) + &
1685	86400	(abs(G%CoriolisBu(I-1,J)) + abs(G%CoriolisBu(I,J-1)))) ! Non-zero on equator!
1686
1687		! u* at the bottom [H T-1 ~> m s-1 or kg m-2 s-1].
1688	86400	ustar = visc%ustar_BBL(i,j)
1689	86400	ustar2 = ustar**2
1690		! In add_drag_diffusivity(), fluxes%ustar_tidal is also added in. There is no
1691		! double-counting because the logic surrounding the calls to add_drag_diffusivity()
1692		! and add_LOTW_BBL_diffusivity() only calls one of the two routines.
1693	86400	if (associated(fluxes%ustar_tidal)) then
1694	0	if (allocated(tv%SpV_avg)) then
1695	0	ustar = ustar + GV%RZ_to_Hrho_bot(i) fluxes%ustar_tidal(i,j)
1696		else
1697	0	ustar = ustar + GV%Z_to_H * fluxes%ustar_tidal(i,j)
1698		endif
1699		endif
1700
1701		! The maximum decay scale should be something of order 200 m. We use the smaller of u*/f and
1702		! (IMax_decay)^-1 as the decay scale. If ustar = 0, this is land so this value doesn't matter.
1703	86400	Idecay = CS%IMax_decay
1704	86400	if ((ustar > 0.0) .and. (absf > CS%IMax_decay * ustar)) Idecay = absf / ustar
1705
1706		! Energy input at the bottom [H Z2 T-3 ~> m3 s-3 or W m-2].
1707		! (Note that visc%BBL_meanKE_loss is in [H L2 T-3 ~> m3 s-3 or W m-2], set in set_BBL_TKE().)
1708	86400	BBL_meanKE_dis = visc%BBL_meanKE_loss(i,j)
1709		! Add in tidal dissipation energy at the bottom [H Z2 T-3 ~> m3 s-3 or W m-2].
1710		! Note that BBL_tidal_dis is in [R Z L2 T-3 ~> W m-2].
1711	86400	if (associated(fluxes%BBL_tidal_dis)) &
1712	0	BBL_meanKE_dis = BBL_meanKE_dis + fluxes%BBL_tidal_dis(i,j) * GV%RZ_to_H
1713	86400	TKE_column = CS%BBL_effic * BBL_meanKE_dis ! Only use a fraction of the mechanical dissipation for mixing.
1714
1715	86400	TKE_remaining = TKE_column
1716	86400	if (CS%LOTW_BBL_answer_date > 20240630) then
1717	86400	dz_above(1) = GV%dz_subroundoff ! This could perhaps be 0 instead.
1718	6566400	do K=2,GV%ke+1
1719	6566400	dz_above(K) = dz_above(K-1) + dz(i,k-1)
1720		enddo
1721	86400	total_depth = dz_above(GV%ke+1)
1722		else
1723	0	total_depth = ( sum(dz(i,:)) + GV%dz_subroundoff ) ! Total column thickness [Z ~> m].
1724		endif
1725	86400	ustar_D = ustar * total_depth
1726	86400	h_bot = 0.
1727	86400	z_bot = 0.
1728	86400	Kd_lower = 0. ! Diffusivity on bottom boundary.
1729
1730		! Work upwards from the bottom, accumulating work used until it exceeds the available TKE input
1731		! at the bottom.
1732	6480720	do K=GV%ke,2,-1
1733	6393600	dz_int = 0.5 * (dz(i,k-1) + dz(i,k))
1734
1735		! Add in additional energy input from bottom-drag against slopes (sides)
1736	6393600	if (Rayleigh_drag) TKE_remaining = TKE_remaining + &
1737		0.5CS%BBL_effic G%IareaT(i,j) * &
1738		(((G%areaCu(I-1,j) * visc%Ray_u(I-1,j,k) * u(I-1,j,k)**2) + &
1739		(G%areaCu(I,j) * visc%Ray_u(I,j,k) * u(I,j,k)**2)) + &
1740		((G%areaCv(i,J-1) * visc%Ray_v(i,J-1,k) * v(i,J-1,k)**2) + &
1741	6393600	(G%areaCv(i,J) * visc%Ray_v(i,J,k) * v(i,J,k)**2)))
1742
1743		! Exponentially decay TKE across the thickness of the layer.
1744		! This is energy loss in addition to work done as mixing, apparently to Joule heating.
1745	6393600	TKE_remaining = exp(-Idecayh(i,j,k)) TKE_remaining
1746
1747	6393600	z_bot = z_bot + dz(i,k) ! Distance between upper interface of layer and the bottom [Z ~> m].
1748	6393600	h_bot = h_bot + h(i,j,k) ! Thickness between upper interface of layer and the bottom [H ~> m or kg m-2].
1749	6393600	if (CS%LOTW_BBL_answer_date > 20240630) then
1750	6393600	D_minus_z = dz_above(K)
1751		else
1752	0	D_minus_z = max(total_depth - z_bot, 0.) ! Distance from the interface to the surface [Z ~> m].
1753		endif
1754
1755		! Diffusivity using law of the wall, limited by rotation, at height z [H Z T-1 ~> m2 s-1 or kg m-1 s-1].
1756		! This calculation is at the upper interface of the layer
1757	6393600	if ( ustar_D + absf * ( h_bot * D_minus_z ) == 0.) then
1758	0	Kd_wall = 0.
1759		else
1760		Kd_wall = ((CS%von_karm * ustar2) * (z_bot * D_minus_z)) &
1761	6393600	/ (ustar_D + absf * (h_bot * D_minus_z))
1762		endif
1763
1764		! TKE associated with Kd_wall [H Z2 T-3 ~> m3 s-3 or W m-2].
1765		! This calculation is for the volume spanning the interface.
1766	6393600	TKE_Kd_wall = Kd_wall * dz_int * max(N2_int(i,K), N2_min)
1767
1768		! Now bound Kd such that the associated TKE is no greater than available TKE for mixing.
1769	6393600	if (TKE_Kd_wall > 0.) then
1770	4452432	TKE_consumed = min(TKE_Kd_wall, TKE_remaining)
1771	4452432	Kd_wall = (TKE_consumed / TKE_Kd_wall) * Kd_wall ! Scale Kd so that only TKE_consumed is used.
1772		else
1773		! Either N2=0 or dh = 0.
1774	1941168	if (TKE_remaining > 0.) then
1775	0	Kd_wall = CS%Kd_max
1776		else
1777	1941168	Kd_wall = 0.
1778		endif
1779	1941168	TKE_consumed = 0.
1780		endif
1781
1782		! Now use up the appropriate about of TKE associated with the diffusivity chosen
1783	6393600	TKE_remaining = TKE_remaining - TKE_consumed ! Note this will be non-negative
1784
1785		! Add this BBL diffusivity to the model net diffusivity.
1786	6393600	Kd_int(i,K) = Kd_int(i,K) + Kd_wall
1787	6393600	if (present(Kd_lay)) Kd_lay(i,k) = Kd_lay(i,k) + 0.5 * (Kd_wall + Kd_lower)
1788	6393600	Kd_lower = Kd_wall ! Store for next layer up.
1789	6480000	if (do_diag_Kd_BBL) Kd_BBL(i,j,K) = Kd_wall
1790		enddo ! k
1791		enddo ! i
1792
1793	720	end subroutine add_LOTW_BBL_diffusivity
1794
1795		!> This routine adds effects of mixed layer radiation to the layer diffusivities.
1796	0	subroutine add_MLrad_diffusivity(dz, fluxes, tv, j, Kd_int, G, GV, US, CS, TKE_to_Kd, Kd_lay)
1797		type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure
1798		type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure
1799		type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
1800		real, dimension(SZI_(G),SZK_(GV)), intent(in) :: dz !< Height change across layers [Z ~> m]
1801		type(forcing), intent(in) :: fluxes !< Surface fluxes structure
1802		type(thermo_var_ptrs), intent(in) :: tv !< Structure containing pointers to any available
1803		!! thermodynamic fields.
1804		integer, intent(in) :: j !< The j-index to work on
1805		real, dimension(SZI_(G),SZK_(GV)+1), intent(inout) :: Kd_int !< The diapycnal diffusivity at interfaces
1806		!! [H Z T-1 ~> m2 s-1 or kg m-1 s-1]
1807		type(set_diffusivity_CS), pointer :: CS !< Diffusivity control structure
1808		real, dimension(SZI_(G),SZK_(GV)), intent(in) :: TKE_to_Kd !< The conversion rate between the TKE
1809		!! TKE dissipated within a layer and the
1810		!! diapycnal diffusivity witin that layer,
1811		!! usually (~Rho_0 / (G_Earth * dRho_lay))
1812		!! [T2 Z-1 ~> s2 m-1]
1813		real, dimension(SZI_(G),SZK_(GV)), &
1814		optional, intent(inout) :: Kd_lay !< The diapycnal diffusivity in layers
1815		!! [H Z T-1 ~> m2 s-1 or kg m-1 s-1].
1816
1817		! This routine adds effects of mixed layer radiation to the layer diffusivities.
1818
1819	0	real, dimension(SZI_(G)) :: h_ml ! Mixed layer thickness [Z ~> m]
1820	0	real, dimension(SZI_(G)) :: TKE_ml_flux ! Mixed layer TKE flux [H Z2 T-3 ~> m3 s-3 or W m-2]
1821	0	real, dimension(SZI_(G)) :: I_decay ! A decay rate [Z-1 ~> m-1].
1822	0	real, dimension(SZI_(G)) :: Kd_mlr_ml ! Diffusivities associated with mixed layer radiation
1823		! [H Z T-1 ~> m2 s-1 or kg m-1 s-1]
1824
1825		real :: f_sq ! The square of the local Coriolis parameter or a related variable [T-2 ~> s-2].
1826		real :: h_ml_sq ! The square of the mixed layer thickness [Z2 ~> m2]
1827		real :: u_star_H ! ustar converted to thickness based units [H T-1 ~> m s-1 or kg m-2 s-1]
1828		real :: ustar_sq ! ustar squared [Z2 T-2 ~> m2 s-2]
1829		real :: Kd_mlr ! A diffusivity associated with mixed layer turbulence radiation
1830		! [H Z T-1 ~> m2 s-1 or kg m-1 s-1]
1831		real :: I_rho ! The inverse of the Boussinesq reference density [R-1 ~> m3 kg-1]
1832		real :: C1_6 ! 1/6 [nondim]
1833		real :: Omega2 ! rotation rate squared [T-2 ~> s-2].
1834		real :: z1 ! layer thickness times I_decay [nondim]
1835		real :: I_decay_len2_TKE ! Squared inverse decay lengthscale for TKE from the bulk mixed
1836		! layer code [Z-2 ~> m-2]
1837		real :: dz_neglect ! A negligibly small height change [Z ~> m]
1838
1839	0	logical :: do_any, do_i(SZI_(G))
1840		integer :: i, k, is, ie, nz, kml
1841	0	is = G%isc ; ie = G%iec ; nz = GV%ke
1842
1843	0	Omega2 = CS%omega**2
1844	0	C1_6 = 1.0 / 6.0
1845	0	kml = GV%nkml
1846	0	dz_neglect = GV%dz_subroundoff
1847	0	I_rho = GV%H_to_Z * GV%RZ_to_H ! == 1.0 / GV%Rho0 ! This is not used when fully non-Boussinesq.
1848
1849	0	if (.not.CS%ML_radiation) return
1850
1851	0	do i=is,ie ; h_ml(i) = 0.0 ; do_i(i) = (G%mask2dT(i,j) > 0.0) ; enddo
1852	0	do k=1,kml ; do i=is,ie ; h_ml(i) = h_ml(i) + dz(i,k) ; enddo ; enddo
1853
1854	0	do i=is,ie ; if (do_i(i)) then
1855	0	if (CS%ML_omega_frac >= 1.0) then
1856	0	f_sq = 4.0 * Omega2
1857		else
1858		f_sq = 0.25 * ((G%Coriolis2Bu(I,J) + G%Coriolis2Bu(I-1,J-1)) + &
1859	0	(G%Coriolis2Bu(I,J-1) + G%Coriolis2Bu(I-1,J)))
1860	0	if (CS%ML_omega_frac > 0.0) &
1861	0	f_sq = CS%ML_omega_frac * 4.0 * Omega2 + (1.0 - CS%ML_omega_frac) * f_sq
1862		endif
1863
1864		! Determine the energy flux out of the mixed layer and its vertical decay scale.
1865	0	if (associated(fluxes%ustar) .and. (GV%Boussinesq .or. .not.associated(fluxes%tau_mag))) then
1866	0	ustar_sq = max(fluxes%ustar(i,j), CS%ustar_min)**2
1867	0	u_star_H = GV%Z_to_H * fluxes%ustar(i,j)
1868	0	elseif (allocated(tv%SpV_avg)) then
1869	0	ustar_sq = max(fluxes%tau_mag(i,j) * tv%SpV_avg(i,j,1), CS%ustar_min**2)
1870	0	u_star_H = GV%RZ_to_H * sqrt(fluxes%tau_mag(i,j) / tv%SpV_avg(i,j,1))
1871		else ! This semi-Boussinesq form is mathematically equivalent to the Boussinesq version above.
1872		! Differs at roundoff: ustar_sq = max(fluxes%tau_mag(i,j) * I_rho, CS%ustar_min**2)
1873	0	ustar_sq = max((sqrt(fluxes%tau_mag(i,j) * I_rho))2, CS%ustar_min2)
1874	0	u_star_H = GV%RZ_to_H * sqrt(fluxes%tau_mag(i,j) * GV%Rho0)
1875		endif
1876	0	TKE_ml_flux(i) = (CS%mstar * CS%ML_rad_coeff) * (ustar_sq * u_star_H)
1877	0	I_decay_len2_TKE = CS%TKE_decay*2 (f_sq / ustar_sq)
1878
1879	0	if (CS%ML_rad_TKE_decay) &
1880	0	TKE_ml_flux(i) = TKE_ml_flux(i) * exp(-h_ml(i) * sqrt(I_decay_len2_TKE))
1881
1882		! Calculate the inverse decay scale
1883	0	h_ml_sq = (CS%ML_rad_efold_coeff * (h_ml(i)+dz_neglect))**2
1884	0	I_decay(i) = sqrt((I_decay_len2_TKE * h_ml_sq + 1.0) / h_ml_sq)
1885
1886		! Average the dissipation layer kml+1, using
1887		! a more accurate Taylor series approximations for very thin layers.
1888	0	z1 = dz(i,kml+1) * I_decay(i)
1889	0	if (z1 > 1e-5) then
1890	0	Kd_mlr = TKE_ml_flux(i) * TKE_to_Kd(i,kml+1) * (1.0 - exp(-z1))
1891		else
1892	0	Kd_mlr = TKE_ml_flux(i) * TKE_to_Kd(i,kml+1) * (z1 * (1.0 - z1 * (0.5 - C1_6 * z1)))
1893		endif
1894	0	Kd_mlr_ml(i) = min(Kd_mlr, CS%ML_rad_kd_max)
1895	0	TKE_ml_flux(i) = TKE_ml_flux(i) * exp(-z1)
1896		endif ; enddo
1897
1898	0	if (present(Kd_lay)) then
1899	0	do k=1,kml+1 ; do i=is,ie ; if (do_i(i)) then
1900	0	Kd_lay(i,k) = Kd_lay(i,k) + Kd_mlr_ml(i)
1901		endif ; enddo ; enddo
1902		endif
1903	0	do K=2,kml+1 ; do i=is,ie ; if (do_i(i)) then
1904	0	Kd_int(i,K) = Kd_int(i,K) + Kd_mlr_ml(i)
1905		endif ; enddo ; enddo
1906	0	if (kml<=nz-1) then ; do i=is,ie ; if (do_i(i)) then
1907	0	Kd_int(i,Kml+2) = Kd_int(i,Kml+2) + 0.5 * Kd_mlr_ml(i)
1908		endif ; enddo ; endif
1909
1910	0	do k=kml+2,nz-1
1911	0	do_any = .false.
1912	0	do i=is,ie ; if (do_i(i)) then
1913	0	z1 = dz(i,k)*I_decay(i)
1914	0	if (CS%ML_Rad_bug) then
1915		! These expressions are dimensionally inconsistent. -RWH
1916		! This is supposed to be the integrated energy deposited in the layer,
1917		! not the average over the layer as in these expressions.
1918	0	if (z1 > 1e-5) then
1919		Kd_mlr = (TKE_ml_flux(i) * TKE_to_Kd(i,k)) * & ! Units of H Z T-1
1920	0	US%m_to_Z * ((1.0 - exp(-z1)) / dz(i,k)) ! Units of m-1
1921		else
1922		Kd_mlr = (TKE_ml_flux(i) * TKE_to_Kd(i,k)) * & ! Units of H Z T-1
1923	0	US%m_to_Z * (I_decay(i) * (1.0 - z1 * (0.5 - C1_6*z1))) ! Units of m-1
1924		endif
1925		else
1926	0	if (z1 > 1e-5) then
1927	0	Kd_mlr = (TKE_ml_flux(i) * TKE_to_Kd(i,k)) * (1.0 - exp(-z1))
1928		else
1929	0	Kd_mlr = (TKE_ml_flux(i) * TKE_to_Kd(i,k)) * (z1 * (1.0 - z1 * (0.5 - C1_6*z1)))
1930		endif
1931		endif
1932	0	Kd_mlr = min(Kd_mlr, CS%ML_rad_kd_max)
1933	0	if (present(Kd_lay)) then
1934	0	Kd_lay(i,k) = Kd_lay(i,k) + Kd_mlr
1935		endif
1936	0	Kd_int(i,K) = Kd_int(i,K) + 0.5 * Kd_mlr
1937	0	Kd_int(i,K+1) = Kd_int(i,K+1) + 0.5 * Kd_mlr
1938
1939	0	TKE_ml_flux(i) = TKE_ml_flux(i) * exp(-z1)
1940	0	if (TKE_ml_flux(i) * I_decay(i) < 0.1 * CS%Kd_min * Omega2) then
1941	0	do_i(i) = .false.
1942	0	else ; do_any = .true. ; endif
1943		endif ; enddo
1944	0	if (.not.do_any) exit
1945		enddo
1946
1947	0	end subroutine add_MLrad_diffusivity
1948
1949		!> This subroutine calculates several properties related to bottom
1950		!! boundary layer turbulence.
1951	12	subroutine set_BBL_TKE(u, v, h, tv, fluxes, visc, G, GV, US, CS, OBC)
1952		type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure
1953		type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure
1954		type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
1955		real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &
1956		intent(in) :: u !< The zonal velocity [L T-1 ~> m s-1]
1957		real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &
1958		intent(in) :: v !< The meridional velocity [L T-1 ~> m s-1]
1959		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
1960		intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2]
1961		type(thermo_var_ptrs), intent(in) :: tv !< Structure with pointers to thermodynamic fields
1962		type(forcing), intent(in) :: fluxes !< A structure of thermodynamic surface fluxes
1963		type(vertvisc_type), intent(inout) :: visc !< Structure containing vertical viscosities, bottom
1964		!! boundary layer properties and related fields.
1965		type(set_diffusivity_CS), pointer :: CS !< Diffusivity control structure
1966		type(ocean_OBC_type), pointer :: OBC !< Open boundaries control structure.
1967
1968		! This subroutine calculates several properties related to bottom
1969		! boundary layer turbulence.
1970
1971		real, dimension(SZI_(G)) :: &
1972	24	htot ! Running sum of the depth in the BBL [Z ~> m].
1973
1974		real, dimension(SZIB_(G)) :: &
1975	24	uhtot, & ! running integral of u in the BBL [Z L T-1 ~> m2 s-1]
1976	24	ustar, & ! bottom boundary layer piston velocity [H T-1 ~> m s-1 or kg m-2 s-1].
1977	24	u2_bbl ! square of the mean zonal velocity in the BBL [L2 T-2 ~> m2 s-2]
1978
1979	24	real :: vhtot(SZI_(G)) ! running integral of v in the BBL [Z L T-1 ~> m2 s-1]
1980
1981		real, dimension(SZI_(G),SZJB_(G)) :: &
1982	24	vstar, & ! ustar at at v-points [H T-1 ~> m s-1 or kg m-2 s-1].
1983	24	v2_bbl ! square of average meridional velocity in BBL [L2 T-2 ~> m2 s-2]
1984		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: &
1985	24	dz ! The vertical distance between interfaces around a layer [Z ~> m]
1986
1987		real :: cdrag_sqrt ! Square root of the drag coefficient [nondim]
1988		real :: hvel ! thickness at velocity points [Z ~> m]
1989
1990	12	logical :: domore, do_i(SZI_(G))
1991		integer :: i, j, k, is, ie, js, je, nz
1992		integer :: l_seg
1993		logical :: local_open_u_BC, local_open_v_BC
1994		logical :: has_obc
1995
1996	12	local_open_u_BC = .false.
1997	12	local_open_v_BC = .false.
1998	12	if (associated(OBC)) then
1999	0	local_open_u_BC = OBC%open_u_BCs_exist_globally
2000	0	local_open_v_BC = OBC%open_v_BCs_exist_globally
2001		endif
2002
2003	12	is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke
2004
2005	12	if (.not.associated(CS)) call MOM_error(FATAL,"set_BBL_TKE: "//&
2006	0	"Module must be initialized before it is used.")
2007
2008	12	if (.not.CS%initialized) call MOM_error(FATAL,"set_BBL_TKE: "//&
2009	0	"Module must be initialized before it is used.")
2010
2011	12	if (.not.CS%bottomdraglaw .or. (CS%BBL_effic<=0.0 .and. CS%ePBL_BBL_effic<=0.0 .and. &
2012		(.not.CS%ePBL_BBL_mstar))) then
2013	0	if (allocated(visc%ustar_BBL)) then
2014	0	do j=js,je ; do i=is,ie ; visc%ustar_BBL(i,j) = 0.0 ; enddo ; enddo
2015		endif
2016	0	if (allocated(visc%BBL_meanKE_loss)) then
2017	0	do j=js,je ; do i=is,ie ; visc%BBL_meanKE_loss(i,j) = 0.0 ; enddo ; enddo
2018		endif
2019	0	if (allocated(visc%BBL_meanKE_loss_sqrtCd)) then
2020	0	do j=js,je ; do i=is,ie ; visc%BBL_meanKE_loss_sqrtCd(i,j) = 0.0 ; enddo ; enddo
2021		endif
2022	0	return
2023		endif
2024
2025	12	cdrag_sqrt = sqrt(CS%cdrag)
2026
2027		! Find the vertical distances across layers.
2028	12	call thickness_to_dz(h, tv, dz, G, GV, US, halo_size=1)
2029
2030		!$OMP parallel default(shared) private(do_i,vhtot,htot,domore,hvel,uhtot,ustar,u2_bbl)
2031		!$OMP do
2032	744	do J=js-1,je
2033		! Determine ustar and the square magnitude of the velocity in the bottom boundary layer.
2034		! Together these give the TKE source and vertical decay scale.
2035	88572	do i=is,ie
2036	88572	do_i(i) = .false. ; vstar(i,J) = 0.0 ; vhtot(i) = 0.0 ; htot(i) = 0.0
2037		enddo
2038	732	if (allocated(visc%Kv_bbl_v)) then
2039	88572	do i=is,ie ; if ((G%mask2dCv(i,J) > 0.0) .and. (cdrag_sqrt*visc%bbl_thick_v(i,J) > 0.0)) then
2040	58728	do_i(i) = .true.
2041	58728	vstar(i,J) = visc%Kv_bbl_v(i,J) / (cdrag_sqrt*visc%bbl_thick_v(i,J))
2042		endif ; enddo
2043		endif
2044		!### What about terms from visc%Ray?
2045
2046	38028	do k=nz,1,-1
2047	38028	domore = .false.
2048	4601388	do i=is,ie ; if (do_i(i)) then
2049		! Determine if grid point is an OBC
2050	1323384	has_obc = .false.
2051	1323384	if (local_open_v_BC) then
2052	0	l_seg = abs(OBC%segnum_v(i,J))
2053	0	if (l_seg /= 0) has_obc = OBC%segment(l_seg)%open
2054		endif
2055
2056		! Compute h based on OBC state
2057	1323384	if (has_obc) then
2058	0	if (OBC%segnum_v(i,J) > 0) then ! OBC_DIRECTION_N
2059	0	hvel = dz(i,j,k)
2060		else
2061	0	hvel = dz(i,j+1,k)
2062		endif
2063		else
2064	1323384	hvel = 0.5*(dz(i,j,k) + dz(i,j+1,k))
2065		endif
2066
2067	1323384	if ((htot(i) + hvel) >= visc%bbl_thick_v(i,J)) then
2068	58728	vhtot(i) = vhtot(i) + (visc%bbl_thick_v(i,J) - htot(i))*v(i,J,k)
2069	58728	htot(i) = visc%bbl_thick_v(i,J)
2070	58728	do_i(i) = .false.
2071		else
2072	1264656	vhtot(i) = vhtot(i) + hvel*v(i,J,k)
2073	1264656	htot(i) = htot(i) + hvel
2074	1264656	domore = .true.
2075		endif
2076		endif ; enddo
2077	38028	if (.not.domore) exit
2078		enddo
2079	88584	do i=is,ie ; if ((G%mask2dCv(i,J) > 0.0) .and. (htot(i) > 0.0)) then
2080	58728	v2_bbl(i,J) = (vhtot(i)vhtot(i)) / (htot(i)htot(i))
2081		else
2082	29112	v2_bbl(i,J) = 0.0
2083		endif ; enddo
2084		enddo
2085		!$OMP do
2086	732	do j=js,je
2087	87840	do I=is-1,ie
2088	87840	do_i(I) = .false. ; ustar(I) = 0.0 ; uhtot(I) = 0.0 ; htot(I) = 0.0
2089		enddo
2090	720	if (allocated(visc%bbl_thick_u)) then
2091	87840	do I=is-1,ie ; if ((G%mask2dCu(I,j) > 0.0) .and. (cdrag_sqrt*visc%bbl_thick_u(I,j) > 0.0)) then
2092	59376	do_i(I) = .true.
2093	59376	ustar(I) = visc%Kv_bbl_u(I,j) / (cdrag_sqrt*visc%bbl_thick_u(I,j))
2094		endif ; enddo
2095		endif
2096
2097	36540	do k=nz,1,-1 ; domore = .false.
2098	4457880	do I=is-1,ie ; if (do_i(I)) then
2099		! Determine if grid point is an OBC
2100	1338876	has_obc = .false.
2101	1338876	if (local_open_u_BC) then
2102	0	l_seg = abs(OBC%segnum_u(I,j))
2103	0	if (l_seg /= 0) has_obc = OBC%segment(l_seg)%open
2104		endif
2105
2106		! Compute h based on OBC state
2107	1338876	if (has_obc) then
2108	0	if (OBC%segnum_u(I,j) > 0) then ! OBC_DIRECTION_E
2109	0	hvel = dz(i,j,k)
2110		else ! OBC_DIRECTION_W
2111	0	hvel = dz(i+1,j,k)
2112		endif
2113		else
2114	1338876	hvel = 0.5*(dz(i,j,k) + dz(i+1,j,k))
2115		endif
2116
2117	1338876	if ((htot(I) + hvel) >= visc%bbl_thick_u(I,j)) then
2118	59376	uhtot(I) = uhtot(I) + (visc%bbl_thick_u(I,j) - htot(I))*u(I,j,k)
2119	59376	htot(I) = visc%bbl_thick_u(I,j)
2120	59376	do_i(I) = .false.
2121		else
2122	1279500	uhtot(I) = uhtot(I) + hvel*u(I,j,k)
2123	1279500	htot(I) = htot(I) + hvel
2124	1279500	domore = .true.
2125		endif
2126		endif ; enddo
2127	36540	if (.not.domore) exit
2128		enddo
2129	87840	do I=is-1,ie ; if ((G%mask2dCu(I,j) > 0.0) .and. (htot(i) > 0.0)) then
2130	59376	u2_bbl(I) = (uhtot(I)uhtot(I)) / (htot(I)htot(I))
2131		else
2132	27744	u2_bbl(I) = 0.0
2133		endif ; enddo
2134
2135	87132	do i=is,ie
2136		visc%ustar_BBL(i,j) = sqrt(0.5G%IareaT(i,j) &
2137		(((G%areaCu(I-1,j)(ustar(I-1)ustar(I-1))) + &
2138		(G%areaCu(I,j)(ustar(I)ustar(I)))) + &
2139		((G%areaCv(i,J-1)(vstar(i,J-1)vstar(i,J-1))) + &
2140	86400	(G%areaCv(i,J)(vstar(i,J)vstar(i,J)))) ) )
2141		visc%BBL_meanKE_loss(i,j) = cdrag_sqrt * &
2142		((((G%areaCu(I-1,j)(ustar(I-1)u2_bbl(I-1))) + &
2143		(G%areaCu(I,j) * (ustar(I)*u2_bbl(I)))) + &
2144		((G%areaCv(i,J-1)(vstar(i,J-1)v2_bbl(i,J-1))) + &
2145	86400	(G%areaCv(i,J) * (vstar(i,J)v2_bbl(i,J)))) )G%IareaT(i,j))
2146		! The following line could be omitted if SET_DIFF_ANSWER_DATE > 20250301 and EPBL_BBL_EFFIC_BUG is false.
2147		visc%BBL_meanKE_loss_sqrtCd(i,j) = &
2148		((((G%areaCu(I-1,j)(ustar(I-1)u2_bbl(I-1))) + &
2149		(G%areaCu(I,j) * (ustar(I)*u2_bbl(I)))) + &
2150		((G%areaCv(i,J-1)(vstar(i,J-1)v2_bbl(i,J-1))) + &
2151	87120	(G%areaCv(i,J) * (vstar(i,J)v2_bbl(i,J)))) )G%IareaT(i,j))
2152		enddo
2153		enddo
2154		!$OMP end parallel
2155
2156		end subroutine set_BBL_TKE
2157
2158	0	subroutine set_density_ratios(h, tv, kb, G, GV, US, CS, j, ds_dsp1, rho_0)
2159		type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure.
2160		type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure.
2161		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
2162		intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2].
2163		type(thermo_var_ptrs), intent(in) :: tv !< Structure containing pointers to any
2164		!! available thermodynamic fields; absent
2165		!! fields have NULL ptrs.
2166		integer, dimension(SZI_(G)), intent(in) :: kb !< Index of lightest layer denser than the buffer
2167		!! layer, or -1 without a bulk mixed layer.
2168		type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
2169		type(set_diffusivity_CS), pointer :: CS !< Control structure returned by previous
2170		!! call to diabatic_entrain_init.
2171		integer, intent(in) :: j !< Meridional index upon which to work.
2172		real, dimension(SZI_(G),SZK_(GV)), intent(out) :: ds_dsp1 !< Coordinate variable (sigma-2)
2173		!! difference across an interface divided by
2174		!! the difference across the interface below
2175		!! it [nondim]
2176		real, dimension(SZI_(G),SZK_(GV)), &
2177		optional, intent(in) :: rho_0 !< Layer potential densities relative to
2178		!! surface press [R ~> kg m-3].
2179
2180		! Local variables
2181		real :: g_R0 ! g_R0 is a rescaled version of g/Rho [L2 Z-1 R-1 T-2 ~> m4 kg-1 s-2]
2182		real :: eps, tmp ! nondimensional temporary variables [nondim]
2183	0	real :: a(SZK_(GV)), a_0(SZK_(GV)) ! nondimensional temporary variables [nondim]
2184	0	real :: p_ref(SZI_(G)) ! an array of tv%P_Ref pressures [R L2 T-2 ~> Pa]
2185	0	real :: Rcv(SZI_(G),SZK_(GV)) ! coordinate density in the mixed and buffer layers [R ~> kg m-3]
2186		real :: I_Drho ! The inverse of the coordinate density difference between
2187		! layers [R-1 ~> m3 kg-1]
2188
2189		integer, dimension(2) :: EOSdom ! The i-computational domain for the equation of state
2190		integer :: i, k, k3, is, ie, nz, kmb
2191	0	is = G%isc ; ie = G%iec ; nz = GV%ke
2192
2193	0	do k=2,nz-1
2194	0	if (GV%g_prime(k+1) /= 0.0) then
2195	0	if (GV%Boussinesq .or. GV%Semi_Boussinesq) then
2196	0	do i=is,ie
2197	0	ds_dsp1(i,k) = GV%g_prime(k) / GV%g_prime(k+1)
2198		enddo
2199		else ! Use a mathematically equivalent form that avoids any dependency on RHO_0.
2200	0	do i=is,ie
2201	0	ds_dsp1(i,k) = (GV%Rlay(k) - GV%Rlay(k-1)) / (GV%Rlay(k+1) - GV%Rlay(k))
2202		enddo
2203		endif
2204		else
2205	0	do i=is,ie
2206	0	ds_dsp1(i,k) = 1.
2207		enddo
2208		endif
2209		enddo
2210
2211	0	if (CS%bulkmixedlayer) then
2212	0	g_R0 = GV%g_Earth / (GV%Rho0)
2213	0	kmb = GV%nk_rho_varies
2214	0	eps = 0.1
2215	0	do i=is,ie ; p_ref(i) = tv%P_Ref ; enddo
2216	0	EOSdom(:) = EOS_domain(G%HI)
2217	0	do k=1,kmb
2218	0	call calculate_density(tv%T(:,j,k), tv%S(:,j,k), p_ref, Rcv(:,k), tv%eqn_of_state, EOSdom)
2219		enddo
2220	0	do i=is,ie
2221	0	if (kb(i) <= nz-1) then
2222		! Set up appropriately limited ratios of the reduced gravities of the
2223		! interfaces above and below the buffer layer and the next denser layer.
2224	0	k = kb(i)
2225
2226	0	if (GV%Boussinesq .or. GV%Semi_Boussinesq) then
2227	0	I_Drho = g_R0 / GV%g_prime(k+1)
2228		else
2229	0	I_Drho = 1.0 / (GV%Rlay(k+1) - GV%Rlay(k))
2230		endif
2231		! The indexing convention for a is appropriate for the interfaces.
2232	0	do k3=1,kmb
2233	0	a(k3+1) = (GV%Rlay(k) - Rcv(i,k3)) * I_Drho
2234		enddo
2235	0	if ((present(rho_0)) .and. (a(kmb+1) < 2.0epsds_dsp1(i,k))) then
2236		! If the buffer layer nearly matches the density of the layer below in the
2237		! coordinate variable (sigma-2), use the sigma-0-based density ratio if it is
2238		! greater (and stable).
2239	0	if ((rho_0(i,k) > rho_0(i,kmb)) .and. &
2240		(rho_0(i,k+1) > rho_0(i,k))) then
2241	0	I_Drho = 1.0 / (rho_0(i,k+1)-rho_0(i,k))
2242	0	a_0(kmb+1) = min((rho_0(i,k)-rho_0(i,kmb)) * I_Drho, ds_dsp1(i,k))
2243	0	if (a_0(kmb+1) > a(kmb+1)) then
2244	0	do k3=2,kmb
2245	0	a_0(k3) = a_0(kmb+1) + (rho_0(i,kmb)-rho_0(i,k3-1)) * I_Drho
2246		enddo
2247	0	if (a(kmb+1) <= eps*ds_dsp1(i,k)) then
2248	0	do k3=2,kmb+1 ; a(k3) = a_0(k3) ; enddo
2249		else
2250		! Alternative... tmp = 0.5(1.0 - cos(PI(a(K2+1)/(eps*ds_dsp1(i,k)) - 1.0)) )
2251	0	tmp = a(kmb+1)/(eps*ds_dsp1(i,k)) - 1.0
2252	0	do k3=2,kmb+1 ; a(k3) = tmpa(k3) + (1.0-tmp)a_0(k3) ; enddo
2253		endif
2254		endif
2255		endif
2256		endif
2257
2258	0	ds_dsp1(i,k) = MAX(a(kmb+1),1e-5)
2259
2260	0	do k3=2,kmb
2261		! ds_dsp1(i,k3) = MAX(a(k3),1e-5)
2262		! Deliberately treat convective instabilities of the upper mixed
2263		! and buffer layers with respect to the deepest buffer layer as
2264		! though they don't exist. They will be eliminated by the upcoming
2265		! call to the mixedlayer code anyway.
2266		! The indexing convention is appropriate for the interfaces.
2267	0	ds_dsp1(i,k3) = MAX(a(k3),ds_dsp1(i,k))
2268		enddo
2269		endif ! (kb(i) <= nz-1)
2270		enddo ! I-loop.
2271		endif ! bulkmixedlayer
2272
2273	0	end subroutine set_density_ratios
2274
2275	1	subroutine set_diffusivity_init(Time, G, GV, US, param_file, diag, CS, int_tide_CSp, halo_TS, &
2276		double_diffuse, physical_OBL_scheme)
2277		type(time_type), intent(in) :: Time !< The current model time
2278		type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure.
2279		type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure.
2280		type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
2281		type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time
2282		!! parameters.
2283		type(diag_ctrl), target, intent(inout) :: diag !< A structure used to regulate diagnostic output.
2284		type(set_diffusivity_CS), pointer :: CS !< pointer set to point to the module control
2285		!! structure.
2286		type(int_tide_CS), pointer :: int_tide_CSp !< Internal tide control structure
2287		integer, intent(out) :: halo_TS !< The halo size of tracer points that must be
2288		!! valid for the calculations in set_diffusivity.
2289		logical, intent(out) :: double_diffuse !< This indicates whether some version
2290		!! of double diffusion is being used.
2291		logical, intent(in) :: physical_OBL_scheme !< If true, a physically based
2292		!! parameterization (like KPP or ePBL or a bulk mixed
2293		!! layer) is used outside of set_diffusivity to
2294		!! specify the mixing that occurs in the ocean's
2295		!! surface boundary layer.
2296
2297		! Local variables
2298		real :: decay_length ! The maximum decay scale for the BBL diffusion [H ~> m or kg m-2]
2299		logical :: ML_use_omega
2300		integer :: default_answer_date ! The default setting for the various ANSWER_DATE flags.
2301		! This include declares and sets the variable "version".
2302		# include "version_variable.h"
2303		character(len=40) :: mdl = "MOM_set_diffusivity" ! This module's name.
2304		real :: vonKar ! The von Karman constant as used for mixed layer viscosity [nondim]
2305		real :: Kd_z ! The background diapycnal diffusivity in [Z2 T-1 ~> m2 s-1] for use
2306		! in setting the default for other diffusivities.
2307		real :: omega_frac_dflt ! The default value for the fraction of the absolute rotation rate
2308		! that is used in place of the absolute value of the local Coriolis
2309		! parameter in the denominator of some expressions [nondim]
2310		logical :: Bryan_Lewis_diffusivity ! If true, the background diapycnal diffusivity uses
2311		! the Bryan-Lewis (1979) style tanh profile.
2312		logical :: use_regridding ! If true, use the ALE algorithm rather than layered
2313		! isopycnal or stacked shallow water mode.
2314		logical :: TKE_to_Kd_used ! If true, TKE_to_Kd and maxTKE need to be calculated.
2315		integer :: is, ie, js, je
2316		integer :: isd, ied, jsd, jed
2317
2318	1	if (associated(CS)) then
2319		call MOM_error(WARNING, "diabatic_entrain_init called with an associated "// &
2320	0	"control structure.")
2321	0	return
2322		endif
2323	1	allocate(CS)
2324
2325	1	CS%initialized = .true.
2326
2327	1	is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec
2328	1	isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed
2329
2330	1	CS%diag => diag
2331	1	CS%int_tide_CSp => int_tide_CSp
2332
2333		! These default values always need to be set.
2334	1	CS%BBL_mixing_as_max = .true.
2335	1	CS%cdrag = 0.003 ; CS%BBL_effic = 0.0
2336	1	CS%bulkmixedlayer = (GV%nkml > 0)
2337
2338		! Read all relevant parameters and write them to the model log.
2339	1	call log_version(param_file, mdl, version, "")
2340
2341	1	call get_param(param_file, mdl, "INPUTDIR", CS%inputdir, default=".")
2342	1	CS%inputdir = slasher(CS%inputdir)
2343		call get_param(param_file, mdl, "FLUX_RI_MAX", CS%FluxRi_max, &
2344		"The flux Richardson number where the stratification is "//&
2345		"large enough that N2 > omega2. The full expression for "//&
2346		"the Flux Richardson number is usually "//&
2347	1	"FLUX_RI_MAX*N2/(N2+OMEGA2).", units="nondim", default=0.2)
2348		call get_param(param_file, mdl, "OMEGA", CS%omega, &
2349	1	"The rotation rate of the earth.", units="s-1", default=7.2921e-5, scale=US%T_to_s)
2350
2351		call get_param(param_file, mdl, "DEFAULT_ANSWER_DATE", default_answer_date, &
2352		"This sets the default value for the various _ANSWER_DATE parameters.", &
2353	1	default=99991231)
2354		call get_param(param_file, mdl, "SET_DIFF_ANSWER_DATE", CS%answer_date, &
2355		"The vintage of the order of arithmetic and expressions in the set diffusivity "//&
2356		"calculations. Values below 20190101 recover the answers from the end of 2018, "//&
2357		"while higher values use updated and more robust forms of the same expressions. "//&
2358		"Values above 20250301 also use less confusing expressions to set the bottom-drag "//&
2359		"generated diffusivity when USE_LOTW_BBL_DIFFUSIVITY is false.", &
2360	1	default=default_answer_date, do_not_log=.not.GV%Boussinesq)
2361	1	if (.not.GV%Boussinesq) CS%answer_date = max(CS%answer_date, 20230701)
2362
2363		! CS%use_tidal_mixing is set to True if an internal tidal dissipation scheme is to be used.
2364		CS%use_tidal_mixing = tidal_mixing_init(Time, G, GV, US, param_file, &
2365	1	CS%int_tide_CSp, diag, CS%tidal_mixing)
2366
2367		call get_param(param_file, mdl, "ML_RADIATION", CS%ML_radiation, &
2368		"If true, allow a fraction of TKE available from wind "//&
2369		"work to penetrate below the base of the mixed layer "//&
2370		"with a vertical decay scale determined by the minimum "//&
2371		"of: (1) The depth of the mixed layer, (2) an Ekman "//&
2372	1	"length scale.", default=.false.)
2373	1	if (CS%ML_radiation) then
2374		! This give a minimum decay scale that is typically much less than Angstrom.
2375	0	CS%ustar_min = 2e-4 * CS%omega * (GV%Angstrom_Z + GV%dZ_subroundoff)
2376
2377		call get_param(param_file, mdl, "ML_RAD_EFOLD_COEFF", CS%ML_rad_efold_coeff, &
2378		"A coefficient that is used to scale the penetration "//&
2379		"depth for turbulence below the base of the mixed layer. "//&
2380	0	"This is only used if ML_RADIATION is true.", units="nondim", default=0.2)
2381		call get_param(param_file, mdl, "ML_RAD_BUG", CS%ML_rad_bug, &
2382		"If true use code with a bug that reduces the energy available "//&
2383		"in the transition layer by a factor of the inverse of the energy "//&
2384	0	"deposition lenthscale (in m).", default=.false.)
2385		call get_param(param_file, mdl, "ML_RAD_KD_MAX", CS%ML_rad_kd_max, &
2386		"The maximum diapycnal diffusivity due to turbulence "//&
2387		"radiated from the base of the mixed layer. "//&
2388		"This is only used if ML_RADIATION is true.", &
2389	0	units="m2 s-1", default=1.0e-3, scale=GV%m2_s_to_HZ_T)
2390		call get_param(param_file, mdl, "ML_RAD_COEFF", CS%ML_rad_coeff, &
2391		"The coefficient which scales MSTAR*USTAR^3 to obtain "//&
2392		"the energy available for mixing below the base of the "//&
2393		"mixed layer. This is only used if ML_RADIATION is true.", &
2394	0	units="nondim", default=0.2)
2395		call get_param(param_file, mdl, "ML_RAD_APPLY_TKE_DECAY", CS%ML_rad_TKE_decay, &
2396		"If true, apply the same exponential decay to ML_rad as "//&
2397		"is applied to the other surface sources of TKE in the "//&
2398	0	"mixed layer code. This is only used if ML_RADIATION is true.", default=.true.)
2399		call get_param(param_file, mdl, "MSTAR", CS%mstar, &
2400		"The ratio of the friction velocity cubed to the TKE "//&
2401	0	"input to the mixed layer.", units="nondim", default=1.2)
2402		call get_param(param_file, mdl, "TKE_DECAY", CS%TKE_decay, &
2403		"The ratio of the natural Ekman depth to the TKE decay scale.", &
2404	0	units="nondim", default=2.5)
2405		call get_param(param_file, mdl, "ML_USE_OMEGA", ML_use_omega, &
2406		"If true, use the absolute rotation rate instead of the "//&
2407		"vertical component of rotation when setting the decay "//&
2408	0	"scale for turbulence.", default=.false., do_not_log=.true.)
2409	0	omega_frac_dflt = 0.0
2410	0	if (ML_use_omega) then
2411	0	call MOM_error(WARNING, "ML_USE_OMEGA is depricated; use ML_OMEGA_FRAC=1.0 instead.")
2412	0	omega_frac_dflt = 1.0
2413		endif
2414		call get_param(param_file, mdl, "ML_OMEGA_FRAC", CS%ML_omega_frac, &
2415		"When setting the decay scale for turbulence, use this "//&
2416		"fraction of the absolute rotation rate blended with the "//&
2417		"local value of f, as sqrt((1-of)f^2 + of4*omega^2).", &
2418	0	units="nondim", default=omega_frac_dflt)
2419		endif
2420
2421		call get_param(param_file, mdl, "BOTTOMDRAGLAW", CS%bottomdraglaw, &
2422		"If true, the bottom stress is calculated with a drag "//&
2423		"law of the form c_drag\|u\|u. The velocity magnitude "//&
2424		"may be an assumed value or it may be based on the actual "//&
2425	1	"velocity in the bottommost HBBL, depending on LINEAR_DRAG.", default=.true.)
2426	1	if (CS%bottomdraglaw) then
2427		call get_param(param_file, mdl, "CDRAG", CS%cdrag, &
2428		"The drag coefficient relating the magnitude of the "//&
2429		"velocity field to the bottom stress. CDRAG is only used "//&
2430	1	"if BOTTOMDRAGLAW is true.", units="nondim", default=0.003)
2431		call get_param(param_file, mdl, "BBL_EFFIC", CS%BBL_effic, &
2432		"The efficiency with which the energy extracted by bottom drag drives BBL "//&
2433		"diffusion. This is only used if BOTTOMDRAGLAW is true.", &
2434	1	units="nondim", default=0.20, scale=US%L_to_Z**2)
2435		call get_param(param_file, mdl, "EPBL_BBL_EFFIC", CS%ePBL_BBL_effic, &
2436	1	units="nondim", default=0.0, do_not_log=.true.)
2437		call get_param(param_file, mdl, "EPBL_BBL_USE_MSTAR", CS%ePBL_BBL_mstar, &
2438	1	default=.false., do_not_log=.true.)
2439		call get_param(param_file, mdl, "BBL_MIXING_MAX_DECAY", decay_length, &
2440		"The maximum decay scale for the BBL diffusion, or 0 to allow the mixing "//&
2441		"to penetrate as far as stratification and rotation permit. The default "//&
2442		"for now is 200 m. This is only used if BOTTOMDRAGLAW is true.", &
2443	1	units="m", default=200.0, scale=GV%m_to_H)
2444
2445	1	CS%IMax_decay = 0.0
2446	1	if (decay_length > 0.0) CS%IMax_decay = 1.0/decay_length
2447		call get_param(param_file, mdl, "BBL_MIXING_AS_MAX", CS%BBL_mixing_as_max, &
2448		"If true, take the maximum of the diffusivity from the "//&
2449		"BBL mixing and the other diffusivities. Otherwise, "//&
2450		"diffusivity from the BBL_mixing is simply added.", &
2451	1	default=.true.)
2452		call get_param(param_file, mdl, "USE_LOTW_BBL_DIFFUSIVITY", CS%use_LOTW_BBL_diffusivity, &
2453		"If true, uses a simple, imprecise but non-coordinate dependent, model "//&
2454		"of BBL mixing diffusivity based on Law of the Wall. Otherwise, uses "//&
2455	1	"the original BBL scheme.", default=.false.)
2456	1	if (CS%use_LOTW_BBL_diffusivity) then
2457		call get_param(param_file, mdl, "LOTW_BBL_USE_OMEGA", CS%LOTW_BBL_use_omega, &
2458		"If true, use the maximum of Omega and N for the TKE to diffusion "//&
2459	1	"calculation. Otherwise, N is N.", default=.true.)
2460		call get_param(param_file, mdl, 'VON_KARMAN_CONST', vonKar, &
2461		'The value the von Karman constant as used for mixed layer viscosity.', &
2462	1	units='nondim', default=0.41)
2463		call get_param(param_file, mdl, 'VON_KARMAN_BBL', CS%von_Karm, &
2464		'The value the von Karman constant as used in calculating the BBL diffusivity.', &
2465	1	units='nondim', default=vonKar)
2466		endif
2467		else
2468	0	CS%use_LOTW_BBL_diffusivity = .false. ! This parameterization depends on a u* from viscous BBL
2469		endif
2470		call get_param(param_file, mdl, "LOTW_BBL_ANSWER_DATE", CS%LOTW_BBL_answer_date, &
2471		"The vintage of the order of arithmetic and expressions in the LOTW_BBL "//&
2472		"calculations. Values below 20240630 recover the original answers, while "//&
2473		"higher values use more accurate expressions. This only applies when "//&
2474		"USE_LOTW_BBL_DIFFUSIVITY is true.", &
2475	1	default=default_answer_date, do_not_log=.not.CS%use_LOTW_BBL_diffusivity)
2476		call get_param(param_file, mdl, "DRAG_DIFFUSIVITY_ANSWER_DATE", CS%drag_diff_answer_date, &
2477		"The vintage of the order of arithmetic in the drag diffusivity calculations. "//&
2478		"Values above 20250301 use less confusing expressions to set the bottom-drag "//&
2479		"generated diffusivity when USE_LOTW_BBL_DIFFUSIVITY is false. ", &
2480	1	default=CS%answer_date, do_not_log=CS%use_LOTW_BBL_diffusivity.or.(CS%BBL_effic<=0.0))
2481
2482		CS%id_Kd_BBL = register_diag_field('ocean_model', 'Kd_BBL', diag%axesTi, Time, &
2483	1	'Bottom Boundary Layer Diffusivity', 'm2 s-1', conversion=GV%HZ_T_to_m2_s)
2484
2485		call get_param(param_file, mdl, "DZ_BBL_AVG_MIN", CS%dz_BBL_avg_min, &
2486		"A minimal distance over which to average to determine the average bottom "//&
2487	1	"boundary layer density.", units="m", default=0.0, scale=US%m_to_Z)
2488
2489		TKE_to_Kd_used = (CS%use_tidal_mixing .or. CS%ML_radiation .or. &
2490	1	(CS%bottomdraglaw .and. .not.CS%use_LOTW_BBL_diffusivity))
2491		call get_param(param_file, mdl, "SIMPLE_TKE_TO_KD", CS%simple_TKE_to_Kd, &
2492		"If true, uses a simple estimate of Kd/TKE that will "//&
2493		"work for arbitrary vertical coordinates. If false, "//&
2494		"calculates Kd/TKE and bounds based on exact energetics "//&
2495		"for an isopycnal layer-formulation.", &
2496	1	default=.false., do_not_log=.not.TKE_to_Kd_used)
2497
2498		call get_param(param_file, mdl, "INTERNAL_TIDES", CS%use_int_tides, &
2499		"If true, use the code that advances a separate set of "//&
2500	1	"equations for the internal tide energy density.", default=.false.)
2501
2502		! set parameters related to the background mixing
2503	1	call bkgnd_mixing_init(Time, G, GV, US, param_file, CS%diag, CS%bkgnd_mixing_csp, physical_OBL_scheme)
2504
2505		call get_param(param_file, mdl, "KV", CS%Kv, &
2506		"The background kinematic viscosity in the interior. "//&
2507		"The molecular value, ~1e-6 m2 s-1, may be used.", &
2508	1	units="m2 s-1", scale=GV%m2_s_to_HZ_T, fail_if_missing=.true.)
2509
2510		call get_param(param_file, mdl, "KD", Kd_z, &
2511		"The background diapycnal diffusivity of density in the "//&
2512		"interior. Zero or the molecular value, ~1e-7 m2 s-1, "//&
2513	1	"may be used.", default=0.0, units="m2 s-1", scale=US%m2_s_to_Z2_T)
2514	1	CS%Kd = (GV%m2_s_to_HZ_TUS%Z2_T_to_m2_s) Kd_z
2515		call get_param(param_file, mdl, "KD_MIN", CS%Kd_min, &
2516		"The minimum diapycnal diffusivity.", &
2517	1	units="m2 s-1", default=0.01Kd_zUS%Z2_T_to_m2_s, scale=GV%m2_s_to_HZ_T)
2518		call get_param(param_file, mdl, "KD_MAX", CS%Kd_max, &
2519		"The maximum permitted increment for the diapycnal "//&
2520		"diffusivity from TKE-based parameterizations, or a negative "//&
2521	1	"value for no limit.", units="m2 s-1", default=-1.0, scale=GV%m2_s_to_HZ_T)
2522	1	if (CS%simple_TKE_to_Kd) then
2523	0	if (CS%Kd_max<=0.) call MOM_error(FATAL, &
2524	0	"set_diffusivity_init: To use SIMPLE_TKE_TO_KD, KD_MAX must be set to >0.")
2525		call get_param(param_file, mdl, "USE_REGRIDDING", use_regridding, &
2526	0	do_not_log=.true., default=.false.)
2527	0	if (use_regridding) call MOM_error(WARNING, &
2528	0	"set_diffusivity_init: SIMPLE_TKE_TO_KD can not be used reliably with USE_REGRIDDING.")
2529		endif
2530
2531		call get_param(param_file, mdl, "KD_ADD", CS%Kd_add, &
2532		"A uniform diapycnal diffusivity that is added "//&
2533		"everywhere without any filtering or scaling.", &
2534	1	units="m2 s-1", default=0.0, scale=GV%m2_s_to_HZ_T)
2535	1	if (CS%use_LOTW_BBL_diffusivity .and. CS%Kd_max<=0.) call MOM_error(FATAL, &
2536		"set_diffusivity_init: KD_MAX must be set (positive) when "// &
2537	0	"USE_LOTW_BBL_DIFFUSIVITY=True.")
2538		call get_param(param_file, mdl, "KD_SMOOTH", CS%Kd_smooth, &
2539		"A diapycnal diffusivity that is used to interpolate "//&
2540		"more sensible values of T & S into thin layers.", &
2541	1	units="m2 s-1", default=1.0e-6, scale=GV%m2_s_to_HZ_T)
2542
2543		call get_param(param_file, mdl, "DEBUG", CS%debug, &
2544		"If true, write out verbose debugging data.", &
2545	1	default=.false., debuggingParam=.true.)
2546
2547		call get_param(param_file, mdl, "USER_CHANGE_DIFFUSIVITY", CS%user_change_diff, &
2548	1	"If true, call user-defined code to change the diffusivity.", default=.false.)
2549
2550		call get_param(param_file, mdl, "DISSIPATION_MIN", CS%dissip_min, &
2551		"The minimum dissipation by which to determine a lower "//&
2552		"bound of Kd (a floor).", &
2553	1	units="W m-3", default=0.0, scale=US%W_m2_to_RZ3_T3*US%Z_to_m)
2554		call get_param(param_file, mdl, "DISSIPATION_N0", CS%dissip_N0, &
2555		"The intercept when N=0 of the N-dependent expression "//&
2556		"used to set a minimum dissipation by which to determine "//&
2557		"a lower bound of Kd (a floor): A in eps_min = A + B*N.", &
2558	1	units="W m-3", default=0.0, scale=US%W_m2_to_RZ3_T3*US%Z_to_m)
2559		call get_param(param_file, mdl, "DISSIPATION_N1", CS%dissip_N1, &
2560		"The coefficient multiplying N, following Gargett, used to "//&
2561		"set a minimum dissipation by which to determine a lower "//&
2562		"bound of Kd (a floor): B in eps_min = A + B*N", &
2563	1	units="J m-3", default=0.0, scale=US%W_m2_to_RZ3_T3US%Z_to_mUS%s_to_T)
2564		call get_param(param_file, mdl, "DISSIPATION_KD_MIN", CS%dissip_Kd_min, &
2565		"The minimum vertical diffusivity applied as a floor.", &
2566	1	units="m2 s-1", default=0.0, scale=GV%m2_s_to_HZ_T)
2567
2568		CS%limit_dissipation = (CS%dissip_min>0.) .or. (CS%dissip_N1>0.) .or. &
2569	1	(CS%dissip_N0>0.) .or. (CS%dissip_Kd_min>0.)
2570	1	CS%dissip_N2 = 0.0
2571	1	if (CS%FluxRi_max > 0.0) &
2572	1	CS%dissip_N2 = CS%dissip_Kd_min * GV%H_to_RZ / CS%FluxRi_max
2573
2574		CS%id_Kd_bkgnd = register_diag_field('ocean_model', 'Kd_bkgnd', diag%axesTi, Time, &
2575	1	'Background diffusivity added by MOM_bkgnd_mixing module', 'm2/s', conversion=GV%HZ_T_to_m2_s)
2576		CS%id_Kv_bkgnd = register_diag_field('ocean_model', 'Kv_bkgnd', diag%axesTi, Time, &
2577	1	'Background viscosity added by MOM_bkgnd_mixing module', 'm2/s', conversion=GV%HZ_T_to_m2_s)
2578
2579	1	if (CS%use_int_tides) then
2580		CS%id_kbbl = register_diag_field('ocean_model', 'kbbl', diag%axesT1, Time, &
2581	0	'BBL index at h points', 'nondim')
2582		CS%id_bbl_thick = register_diag_field('ocean_model', 'bbl_thick', diag%axesT1, Time, &
2583	0	'BBL thickness at h points', 'm', conversion=US%Z_to_m)
2584		CS%id_Kd_leak = register_diag_field('ocean_model', 'Kd_leak', diag%axesTi, Time, &
2585	0	'internal tides leakage viscosity added by MOM_internal tides module', 'm2/s', conversion=GV%HZ_T_to_m2_s)
2586		CS%id_Kd_Froude = register_diag_field('ocean_model', 'Kd_Froude', diag%axesTi, Time, &
2587	0	'internal tides Froude viscosity added by MOM_internal tides module', 'm2/s', conversion=GV%HZ_T_to_m2_s)
2588		CS%id_Kd_itidal = register_diag_field('ocean_model', 'Kd_itidal', diag%axesTi, Time, &
2589	0	'internal tides wave drag viscosity added by MOM_internal tides module', 'm2/s', conversion=GV%HZ_T_to_m2_s)
2590		CS%id_Kd_quad = register_diag_field('ocean_model', 'Kd_quad', diag%axesTi, Time, &
2591	0	'internal tides bottom viscosity added by MOM_internal tides module', 'm2/s', conversion=GV%HZ_T_to_m2_s)
2592		CS%id_Kd_slope = register_diag_field('ocean_model', 'Kd_slope', diag%axesTi, Time, &
2593	0	'internal tides slope viscosity added by MOM_internal tides module', 'm2/s', conversion=GV%HZ_T_to_m2_s)
2594		CS%id_prof_leak = register_diag_field('ocean_model', 'prof_leak', diag%axesTl, Time, &
2595	0	'internal tides leakage profile added by MOM_internal tides module', 'm-1', conversion=GV%m_to_H)
2596		CS%id_prof_Froude = register_diag_field('ocean_model', 'prof_Froude', diag%axesTl, Time, &
2597	0	'internal tides Froude profile added by MOM_internal tides module', 'm-1', conversion=GV%m_to_H)
2598		CS%id_prof_itidal = register_diag_field('ocean_model', 'prof_itidal', diag%axesTl, Time, &
2599	0	'internal tides wave drag profile added by MOM_internal tides module', 'm-1', conversion=GV%m_to_H)
2600		CS%id_prof_quad = register_diag_field('ocean_model', 'prof_quad', diag%axesTl, Time, &
2601	0	'internal tides bottom profile added by MOM_internal tides module', 'm-1', conversion=GV%m_to_H)
2602		CS%id_prof_slope = register_diag_field('ocean_model', 'prof_slope', diag%axesTl, Time, &
2603	0	'internal tides slope profile added by MOM_internal tides module', 'm-1', conversion=GV%m_to_H)
2604		endif
2605
2606		CS%id_Kd_layer = register_diag_field('ocean_model', 'Kd_layer', diag%axesTL, Time, &
2607	1	'Diapycnal diffusivity of layers (as set)', 'm2 s-1', conversion=GV%HZ_T_to_m2_s)
2608
2609	1	if (CS%use_tidal_mixing) then
2610		CS%id_Kd_Work = register_diag_field('ocean_model', 'Kd_Work', diag%axesTL, Time, &
2611	0	'Work done by Diapycnal Mixing', 'W m-2', conversion=US%RZ3_T3_to_W_m2)
2612		CS%id_Kd_Work_added = register_diag_field('ocean_model', 'Kd_Work_added', diag%axesTL, Time, &
2613	0	'Work done by additional mixing Kd_add', 'W m-2', conversion=US%RZ3_T3_to_W_m2)
2614		CS%id_maxTKE = register_diag_field('ocean_model', 'maxTKE', diag%axesTL, Time, &
2615	0	'Maximum layer TKE', 'm3 s-3', conversion=(GV%H_to_mUS%Z_to_m2US%s_to_T**3))
2616		CS%id_TKE_to_Kd = register_diag_field('ocean_model', 'TKE_to_Kd', diag%axesTL, Time, &
2617	0	'Convert TKE to Kd', 's2 m', conversion=GV%HZ_T_to_m2_s(GV%m_to_HUS%m_to_Z*2US%T_to_s**3))
2618		CS%id_N2 = register_diag_field('ocean_model', 'N2', diag%axesTi, Time, &
2619		'Buoyancy frequency squared', 's-2', conversion=US%s_to_T**2, cmor_field_name='obvfsq', &
2620		cmor_long_name='Square of seawater buoyancy frequency', &
2621	0	cmor_standard_name='square_of_brunt_vaisala_frequency_in_sea_water')
2622		endif
2623
2624	1	if (CS%user_change_diff) &
2625		CS%id_Kd_user = register_diag_field('ocean_model', 'Kd_user', diag%axesTi, Time, &
2626	0	'User-specified Extra Diffusivity', 'm2 s-1', conversion=GV%HZ_T_to_m2_s)
2627
2628		call get_param(param_file, mdl, "DOUBLE_DIFFUSION", CS%double_diffusion, &
2629		"If true, increase diffusivities for temperature or salinity based on the "//&
2630		"double-diffusive parameterization described in Large et al. (1994).", &
2631	1	default=.false.)
2632
2633	1	if (CS%double_diffusion) then
2634		call get_param(param_file, mdl, "MAX_RRHO_SALT_FINGERS", CS%Max_Rrho_salt_fingers, &
2635		"Maximum density ratio for salt fingering regime.", &
2636	0	default=1.9, units="nondim")
2637		call get_param(param_file, mdl, "MAX_SALT_DIFF_SALT_FINGERS", CS%Max_salt_diff_salt_fingers, &
2638		"Maximum salt diffusivity for salt fingering regime.", &
2639	0	default=1.e-4, units="m2 s-1", scale=GV%m2_s_to_HZ_T)
2640		call get_param(param_file, mdl, "KV_MOLECULAR", CS%Kv_molecular, &
2641		"Molecular viscosity for calculation of fluxes under double-diffusive "//&
2642	0	"convection.", default=1.5e-6, units="m2 s-1", scale=GV%m2_s_to_HZ_T)
2643		! The default molecular viscosity follows the CCSM4.0 and MOM4p1 defaults.
2644		endif ! old double-diffusion
2645
2646	1	if (CS%user_change_diff) then
2647	0	call user_change_diff_init(Time, G, GV, US, param_file, diag, CS%user_change_diff_CSp)
2648		endif
2649
2650		call get_param(param_file, mdl, "BRYAN_LEWIS_DIFFUSIVITY", Bryan_Lewis_diffusivity, &
2651		"If true, use a Bryan & Lewis (JGR 1979) like tanh "//&
2652		"profile of background diapycnal diffusivity with depth. "//&
2653	1	"This is done via CVMix.", default=.false., do_not_log=.true.)
2654	1	if (CS%use_tidal_mixing .and. Bryan_Lewis_diffusivity) &
2655		call MOM_error(FATAL,"MOM_Set_Diffusivity: "// &
2656	0	"Bryan-Lewis and internal tidal dissipation are both enabled. Choose one.")
2657
2658	1	CS%useKappaShear = kappa_shear_init(Time, G, GV, US, param_file, CS%diag, CS%kappaShear_CSp)
2659	1	CS%Vertex_Shear = kappa_shear_at_vertex(param_file)
2660
2661	1	if (CS%useKappaShear) &
2662	1	id_clock_kappaShear = cpu_clock_id('(Ocean kappa_shear)', grain=CLOCK_MODULE)
2663
2664		! CVMix shear-driven mixing
2665	1	CS%use_CVMix_shear = CVMix_shear_init(Time, G, GV, US, param_file, CS%diag, CS%CVMix_shear_csp)
2666
2667		! CVMix double diffusion mixing
2668	1	CS%use_CVMix_ddiff = CVMix_ddiff_init(Time, G, GV, US, param_file, CS%diag, CS%CVMix_ddiff_csp)
2669	1	if (CS%use_CVMix_ddiff) &
2670	0	id_clock_CVMix_ddiff = cpu_clock_id('(Double diffusion via CVMix)', grain=CLOCK_MODULE)
2671
2672	1	if (CS%double_diffusion .and. CS%use_CVMix_ddiff) then
2673		call MOM_error(FATAL, 'set_diffusivity_init: '// &
2674		'Multiple double-diffusion options selected (DOUBLE_DIFFUSION and '//&
2675	0	'USE_CVMIX_DDIFF), please disable all but one option to proceed.')
2676		endif
2677
2678	1	if (CS%double_diffusion .or. CS%use_CVMix_ddiff) then
2679		CS%id_KT_extra = register_diag_field('ocean_model', 'KT_extra', diag%axesTi, Time, &
2680	0	'Double-diffusive diffusivity for temperature', 'm2 s-1', conversion=GV%HZ_T_to_m2_s)
2681		CS%id_KS_extra = register_diag_field('ocean_model', 'KS_extra', diag%axesTi, Time, &
2682	0	'Double-diffusive diffusivity for salinity', 'm2 s-1', conversion=GV%HZ_T_to_m2_s)
2683		endif
2684	1	if (CS%use_CVMix_ddiff) then
2685		CS%id_R_rho = register_diag_field('ocean_model', 'R_rho', diag%axesTi, Time, &
2686	0	'Double-diffusion density ratio', 'nondim')
2687		endif
2688
2689	1	halo_TS = 0
2690	1	if (CS%Vertex_Shear) halo_TS = 1
2691
2692	1	double_diffuse = (CS%double_diffusion .or. CS%use_CVMix_ddiff)
2693
2694		end subroutine set_diffusivity_init
2695
2696		!> Clear pointers and deallocate memory
2697	1	subroutine set_diffusivity_end(CS)
2698		type(set_diffusivity_CS), intent(inout) :: CS !< Control structure for this module
2699
2700	1	call bkgnd_mixing_end(CS%bkgnd_mixing_csp)
2701
2702	1	if (CS%use_tidal_mixing) &
2703	0	call tidal_mixing_end(CS%tidal_mixing)
2704
2705	1	if (CS%user_change_diff) call user_change_diff_end(CS%user_change_diff_CSp)
2706
2707	1	if (associated(CS%CVMix_ddiff_CSp)) deallocate(CS%CVMix_ddiff_CSp)
2708
2709	1	if (CS%use_CVMix_shear) then
2710	0	call CVMix_shear_end(CS%CVMix_shear_CSp)
2711	0	deallocate(CS%CVMix_shear_CSp)
2712		endif
2713
2714		! NOTE: CS%kappaShear_CSp is always allocated, even if unused
2715	1	deallocate(CS%kappaShear_CSp)
2716	1	end subroutine set_diffusivity_end
2717
2718	0	end module MOM_set_diffusivity