MOM_set_diffusivity.F90

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