src/parameterizations/vertical/MOM_kappa_shear.F90

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

1		! This file is part of MOM6, the Modular Ocean Model version 6.
2		! See the LICENSE file for licensing information.
3		! SPDX-License-Identifier: Apache-2.0
4
5		!> Shear-dependent mixing following Jackson et al. 2008.
6		module MOM_kappa_shear
7
8		use MOM_cpu_clock, only : cpu_clock_id, cpu_clock_begin, cpu_clock_end
9		use MOM_cpu_clock, only : CLOCK_MODULE_DRIVER, CLOCK_MODULE, CLOCK_ROUTINE
10		use MOM_diag_mediator, only : post_data, register_diag_field, safe_alloc_ptr
11		use MOM_diag_mediator, only : diag_ctrl, time_type
12		use MOM_debugging, only : hchksum, Bchksum
13		use MOM_error_handler, only : MOM_error, is_root_pe, FATAL, WARNING, NOTE
14		use MOM_file_parser, only : get_param, log_version, param_file_type
15		use MOM_grid, only : ocean_grid_type
16		use MOM_interface_heights, only : thickness_to_dz
17		use MOM_unit_scaling, only : unit_scale_type
18		use MOM_variables, only : thermo_var_ptrs
19		use MOM_verticalGrid, only : verticalGrid_type
20		use MOM_EOS, only : calculate_density_derivs
21		use MOM_EOS, only : calculate_density, calculate_specific_vol_derivs
22
23		implicit none ; private
24
25		#include <MOM_memory.h>
26
27		public Calculate_kappa_shear, Calc_kappa_shear_vertex, kappa_shear_init
28		public kappa_shear_is_used, kappa_shear_at_vertex
29
30		! A note on unit descriptions in comments: MOM6 uses units that can be rescaled for dimensional
31		! consistency testing. These are noted in comments with units like Z, H, L, and T, along with
32		! their mks counterparts with notation like "a velocity [Z T-1 ~> m s-1]". If the units
33		! vary with the Boussinesq approximation, the Boussinesq variant is given first.
34
35		!> This control structure holds the parameters that regulate shear mixing
36		type, public :: Kappa_shear_CS ; private
37		real :: RiNo_crit !< The critical shear Richardson number for
38		!! shear-entrainment [nondim]. The theoretical value is 0.25.
39		!! The values found by Jackson et al. are 0.25-0.35.
40		real :: Shearmix_rate !< A nondimensional rate scale for shear-driven
41		!! entrainment [nondim]. The value given by Jackson et al.
42		!! is 0.085-0.089.
43		real :: FRi_curvature !< A constant giving the curvature of the function
44		!! of the Richardson number that relates shear to
45		!! sources in the kappa equation [nondim].
46		!! The values found by Jackson et al. are -0.97 - -0.89.
47		real :: C_N !< The coefficient for the decay of TKE due to
48		!! stratification (i.e. proportional to N*tke) [nondim].
49		!! The values found by Jackson et al. are 0.24-0.28.
50		real :: C_S !< The coefficient for the decay of TKE due to
51		!! shear (i.e. proportional to \|S\|*tke) [nondim].
52		!! The values found by Jackson et al. are 0.14-0.12.
53		real :: lambda !< The coefficient for the buoyancy length scale
54		!! in the kappa equation [nondim].
55		!! The values found by Jackson et al. are 0.82-0.81.
56		real :: lambda2_N_S !< The square of the ratio of the coefficients of
57		!! the buoyancy and shear scales in the diffusivity
58		!! equation, 0 to eliminate the shear scale [nondim].
59		real :: lz_rescale !< A coefficient to rescale the distance to the nearest
60		!! solid boundary. This adjustment is to account for
61		!! regions where 3 dimensional turbulence prevents the
62		!! growth of shear instabilities [nondim].
63		real :: TKE_bg !< The background level of TKE [Z2 T-2 ~> m2 s-2].
64		real :: kappa_0 !< The background diapycnal diffusivity [H Z T-1 ~> m2 s-1 or Pa s]
65		real :: kappa_seed !< A moderately large seed value of diapycnal diffusivity that
66		!! is used as a starting turbulent diffusivity in the iterations
67		!! to finding an energetically constrained solution for the
68		!! shear-driven diffusivity [H Z T-1 ~> m2 s-1 or Pa s]
69		real :: kappa_trunc !< Diffusivities smaller than this are rounded to 0 [H Z T-1 ~> m2 s-1 or Pa s]
70		real :: kappa_tol_err !< The fractional error in kappa that is tolerated [nondim].
71		real :: Prandtl_turb !< Prandtl number used to convert Kd_shear into viscosity [nondim].
72		integer :: nkml !< The number of layers in the mixed layer, as
73		!! treated in this routine. If the pieces of the
74		!! mixed layer are not to be treated collectively,
75		!! nkml is set to 1.
76		integer :: max_RiNo_it !< The maximum number of iterations that may be used
77		!! to estimate the instantaneous shear-driven mixing.
78		integer :: max_KS_it !< The maximum number of iterations that may be used
79		!! to estimate the time-averaged diffusivity.
80		logical :: dKdQ_iteration_bug !< If true. use an older, dimensionally inconsistent estimate of
81		!! the derivative of diffusivity with energy in the Newton's method
82		!! iteration. The bug causes under-corrections when dz > 1m.
83		logical :: KS_at_vertex !< If true, do the calculations of the shear-driven mixing
84		!! at the cell vertices (i.e., the vorticity points).
85		logical :: eliminate_massless !< If true, massless layers are merged with neighboring
86		!! massive layers in this calculation.
87		! I can think of no good reason why this should be false. - RWH
88		real :: vel_underflow !< Velocity components smaller than vel_underflow
89		!! are set to 0 [L T-1 ~> m s-1].
90		real :: kappa_src_max_chg !< The maximum permitted increase in the kappa source within an
91		!! iteration relative to the local source [nondim]. This must be
92		!! greater than 1. The lower limit for the permitted fractional
93		!! decrease is (1 - 0.5/kappa_src_max_chg). These limits could
94		!! perhaps be made dynamic with an improved iterative solver.
95		real :: VS_GeoMean_Kdmin !< A minimum diffusivity for computing the horizontal averages
96		!! when using the geometric mean with VERTEX_SHEAR=True. The model
97		!! is sensitive to this value, which is a drawback of using the
98		!! geometric average as currently implemented.
99		logical :: psurf_bug !< If true, do a simple average of the cell surface pressures to get a
100		!! surface pressure at the corner if VERTEX_SHEAR=True. Otherwise mask
101		!! out any land points in the average.
102		logical :: all_layer_TKE_bug !< If true, report back the latest estimate of TKE instead of the
103		!! time average TKE when there is mass in all layers. Otherwise always
104		!! report the time-averaged TKE, as is currently done when there
105		!! are some massless layers.
106		logical :: VS_viscosity_bug !< If true, use a bug in the calculation of the viscosity that sets
107		!! it to zero for all vertices that are on a coastline.
108		logical :: vertex_shear_OBC_bug !< If false, use extra masking when interpolating thicknesses to velocity
109		!! points for setting up the shear velocities at vertices to avoid using
110		!! external thicknesses at open boundaries. When OBCs are not in use,
111		!! this parameter does not change answers, but true is more efficient.
112		logical :: VS_GeometricMean !< If true use geometric averaging for Kd from vertices to tracer points
113		logical :: VS_ThicknessMean !< If true use thickness weighting when averaging Kd from vertices to
114		!! tracer points
115		logical :: restrictive_tolerance_check !< If false, uses the less restrictive tolerance check to
116		!! determine if a timestep is acceptable for the KS_it outer iteration
117		!! loop, as the code was originally written. True uses the more
118		!! restrictive check.
119		! logical :: layer_stagger = .false. ! If true, do the calculations centered at
120		! layers, rather than the interfaces.
121		logical :: debug = .false. !< If true, write verbose debugging messages.
122		type(diag_ctrl), pointer :: diag => NULL() !< A structure that is used to
123		!! regulate the timing of diagnostic output.
124		!>@{ Diagnostic IDs
125		integer :: id_Kd_shear = -1, id_TKE = -1, id_Kd_vertex = -1, &
126		id_S2_init = -1, id_N2_init = -1, id_S2_mean = -1, id_N2_mean = -1
127		!>@}
128		end type Kappa_shear_CS
129
130		! integer :: id_clock_project, id_clock_KQ, id_clock_avg, id_clock_setup
131
132		contains
133
134		!> Subroutine for calculating shear-driven diffusivity and TKE in tracer columns
135	1	subroutine Calculate_kappa_shear(u_in, v_in, h, tv, p_surf, kappa_io, tke_io, &
136	1	kv_io, dt, G, GV, US, CS)
137		type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure.
138		type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure.
139		type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
140		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
141		intent(in) :: u_in !< Initial zonal velocity [L T-1 ~> m s-1].
142		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
143		intent(in) :: v_in !< Initial meridional velocity [L T-1 ~> m s-1].
144		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
145		intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2].
146		type(thermo_var_ptrs), intent(in) :: tv !< A structure containing pointers to any
147		!! available thermodynamic fields. Absent fields
148		!! have NULL ptrs.
149		real, dimension(:,:), pointer :: p_surf !< The pressure at the ocean surface [R L2 T-2 ~> Pa] (or NULL).
150		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), &
151		intent(inout) :: kappa_io !< The diapycnal diffusivity at each interface
152		!! [H Z T-1 ~> m2 s-1 or kg m-1 s-1]. Initially this
153		!! is the value from the previous timestep, which may
154		!! accelerate the iteration toward convergence.
155		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), &
156		intent(out) :: tke_io !< The turbulent kinetic energy per unit mass at
157		!! each interface (not layer!) [Z2 T-2 ~> m2 s-2].
158		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), &
159		intent(inout) :: kv_io !< The vertical viscosity at each interface
160		!! (not layer!) [H Z T-1 ~> m2 s-1 or Pa s]. This discards any
161		!! previous value (i.e. it is intent out) and
162		!! simply sets Kv = Prandtl * Kd_shear
163		real, intent(in) :: dt !< Time increment [T ~> s].
164		type(Kappa_shear_CS), pointer :: CS !< The control structure returned by a previous
165		!! call to kappa_shear_init.
166
167		! Local variables
168		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1) :: &
169	2	diag_N2_init, & ! Diagnostic of N2 as provided to this routine [T-2 ~> s-2]
170	2	diag_S2_init, & ! Diagnostic of S2 as provided to this routine [T-2 ~> s-2]
171	2	diag_N2_mean, & ! Diagnostic of N2 averaged over the timestep applied [T-2 ~> s-2]
172	2	diag_S2_mean ! Diagnostic of S2 averaged over the timestep applied [T-2 ~> s-2]
173		real, dimension(SZI_(G),SZK_(GV)) :: &
174	2	h_2d, & ! A 2-D version of h [H ~> m or kg m-2].
175	2	dz_2d, & ! Vertical distance between interface heights [Z ~> m].
176	2	u_2d, v_2d, & ! 2-D versions of u_in and v_in, converted to [L T-1 ~> m s-1].
177	2	T_2d, S_2d, rho_2d ! 2-D versions of T [C ~> degC], S [S ~> ppt], and rho [R ~> kg m-3].
178		real, dimension(SZI_(G),SZK_(GV)+1) :: &
179	2	kappa_2d, & ! 2-D version of kappa_io [H Z T-1 ~> m2 s-1 or Pa s]
180	2	tke_2d ! 2-D version tke_io [Z2 T-2 ~> m2 s-2].
181		real, dimension(SZK_(GV)) :: &
182	2	Idz, & ! The inverse of the thickness of the merged layers [H-1 ~> m2 kg-1].
183	2	h_lay, & ! The layer thickness [H ~> m or kg m-2]
184	2	dz_lay, & ! The geometric layer thickness in height units [Z ~> m]
185	2	u0xdz, & ! The initial zonal velocity times dz [H L T-1 ~> m2 s-1 or kg m-1 s-1]
186	2	v0xdz, & ! The initial meridional velocity times dz [H L T-1 ~> m2 s-1 or kg m-1 s-1]
187	2	T0xdz, & ! The initial temperature times thickness [C H ~> degC m or degC kg m-2] or if
188		! temperature is not a state variable, the density times thickness [R H ~> kg m-2 or kg2 m-5]
189	2	S0xdz ! The initial salinity times dz [S H ~> ppt m or ppt kg m-2].
190		real, dimension(SZK_(GV)+1) :: &
191	2	kappa, & ! The shear-driven diapycnal diffusivity at an interface [H Z T-1 ~> m2 s-1 or Pa s]
192	2	tke, & ! The Turbulent Kinetic Energy per unit mass at an interface [Z2 T-2 ~> m2 s-2].
193	2	kappa_avg, & ! The time-weighted average of kappa [H Z T-1 ~> m2 s-1 or Pa s]
194	2	tke_avg, & ! The time-weighted average of TKE [Z2 T-2 ~> m2 s-2].
195	2	N2_init, & ! N2 as provided to this routine [T-2 ~> s-2].
196	2	S2_init, & ! S2 as provided to this routine [T-2 ~> s-2].
197	2	N2_mean, & ! The time-weighted average of N2 [T-2 ~> s-2].
198	2	S2_mean ! The time-weighted average of S2 [T-2 ~> s-2].
199
200		real :: f2 ! The squared Coriolis parameter of each column [T-2 ~> s-2].
201		real :: surface_pres ! The top surface pressure [R L2 T-2 ~> Pa].
202
203		real :: dz_in_lay ! The running sum of the thickness in a layer [H ~> m or kg m-2]
204		real :: k0dt ! The background diffusivity times the timestep [H Z ~> m2 or kg m-1]
205		real :: dz_massless ! A layer thickness that is considered massless [H ~> m or kg m-2]
206		logical :: use_temperature ! If true, temperature and salinity have been
207		! allocated and are being used as state variables.
208
209	2	integer, dimension(SZK_(GV)+1) :: kc ! The index map between the original
210		! interfaces and the interfaces with massless layers
211		! merged into nearby massive layers.
212	2	real, dimension(SZK_(GV)+1) :: kf ! The fractional weight of interface kc+1 for
213		! interpolating back to the original index space [nondim].
214		integer :: is, ie, js, je, i, j, k, nz, nzc
215
216	1	is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke
217
218	1	use_temperature = associated(tv%T)
219
220	1	k0dt = dt*CS%kappa_0
221	1	dz_massless = 0.1sqrt((US%Z_to_mGV%m_to_H)*k0dt)
222
223	1	if ((CS%id_N2_init>0) .or. CS%debug) diag_N2_init(:,:,:) = 0.0
224	1	if ((CS%id_S2_init>0) .or. CS%debug) diag_S2_init(:,:,:) = 0.0
225	1	if (CS%id_N2_mean>0) diag_N2_mean(:,:,:) = 0.0
226	1	if (CS%id_S2_mean>0) diag_S2_mean(:,:,:) = 0.0
227
228		!$OMP parallel do default(private) shared(js,je,is,ie,nz,h,u_in,v_in,use_temperature,tv,G,GV,US, &
229		!$OMP CS,kappa_io,dz_massless,k0dt,p_surf,dt,tke_io,kv_io, &
230		!$OMP diag_N2_init,diag_S2_init,diag_N2_mean,diag_S2_mean)
231	181	do j=js,je
232
233		! Convert layer thicknesses into geometric thickness in height units.
234	180	call thickness_to_dz(h, tv, dz_2d, j, G, GV)
235
236	1429740	do k=1,nz ; do i=is,ie
237	1425600	h_2d(i,k) = h(i,j,k)
238	1429560	u_2d(i,k) = u_in(i,j,k) ; v_2d(i,k) = v_in(i,j,k)
239		enddo ; enddo
240	1429740	if (use_temperature) then ; do k=1,nz ; do i=is,ie
241	1429560	T_2d(i,k) = tv%T(i,j,k) ; S_2d(i,k) = tv%S(i,j,k)
242	0	enddo ; enddo ; else ; do k=1,nz ; do i=is,ie
243	0	rho_2d(i,k) = GV%Rlay(k) ! Could be tv%Rho(i,j,k) ?
244		enddo ; enddo ; endif
245
246		!---------------------------------------
247		! Work on each column.
248		!---------------------------------------
249	64980	do i=is,ie ; if (G%mask2dT(i,j) > 0.0) then
250		! call cpu_clock_begin(id_clock_setup)
251
252		! Store a transposed version of the initial arrays.
253		! Any elimination of massless layers would occur here.
254	46250	if (CS%eliminate_massless) then
255	46250	nzc = 1
256	1063750	do k=1,nz
257		! Zero out the thicknesses of all layers, even if they are unused.
258	1017500	h_lay(k) = 0.0 ; dz_lay(k) = 0.0 ; u0xdz(k) = 0.0 ; v0xdz(k) = 0.0
259	1017500	T0xdz(k) = 0.0 ; S0xdz(k) = 0.0
260
261		! Add a new layer if this one has mass.
262		! if ((h_lay(nzc) > 0.0) .and. (h_2d(i,k) > dz_massless)) nzc = nzc+1
263	1017500	if ((k>CS%nkml) .and. (h_lay(nzc) > 0.0) .and. &
264	345393	(h_2d(i,k) > dz_massless)) nzc = nzc+1
265
266		! Only merge clusters of massless layers.
267		! if ((h_lay(nzc) > dz_massless) .or. &
268		! ((h_lay(nzc) > 0.0) .and. (h_2d(i,k) > dz_massless))) nzc = nzc+1
269
270	1017500	kc(k) = nzc
271	1017500	h_lay(nzc) = h_lay(nzc) + h_2d(i,k)
272	1017500	dz_lay(nzc) = dz_lay(nzc) + dz_2d(i,k)
273	1017500	u0xdz(nzc) = u0xdz(nzc) + u_2d(i,k)*h_2d(i,k)
274	1017500	v0xdz(nzc) = v0xdz(nzc) + v_2d(i,k)*h_2d(i,k)
275	1063750	if (use_temperature) then
276	1017500	T0xdz(nzc) = T0xdz(nzc) + T_2d(i,k)*h_2d(i,k)
277	1017500	S0xdz(nzc) = S0xdz(nzc) + S_2d(i,k)*h_2d(i,k)
278		else
279	0	T0xdz(nzc) = T0xdz(nzc) + rho_2d(i,k)*h_2d(i,k)
280	0	S0xdz(nzc) = S0xdz(nzc) + rho_2d(i,k)*h_2d(i,k)
281		endif
282		enddo
283	46250	kc(nz+1) = nzc+1
284
285		! Set up Idz as the inverse of layer thicknesses.
286	437893	do k=1,nzc ; Idz(k) = 1.0 / h_lay(k) ; enddo
287
288		! Now determine kf, the fractional weight of interface kc when
289		! interpolating between interfaces kc and kc+1.
290	46250	kf(1) = 0.0 ; dz_in_lay = h_2d(i,1)
291	1017500	do k=2,nz
292	1017500	if (kc(k) > kc(k-1)) then
293	345393	kf(k) = 0.0 ; dz_in_lay = h_2d(i,k)
294		else
295	625857	kf(k) = dz_in_lay*Idz(kc(k)) ; dz_in_lay = dz_in_lay + h_2d(i,k)
296		endif
297		enddo
298	46250	kf(nz+1) = 0.0
299		else
300	0	do k=1,nz
301	0	h_lay(k) = h_2d(i,k)
302	0	dz_lay(k) = dz_2d(i,k)
303	0	u0xdz(k) = u_2d(i,k)h_lay(k) ; v0xdz(k) = v_2d(i,k)h_lay(k)
304		enddo
305	0	if (use_temperature) then
306	0	do k=1,nz
307	0	T0xdz(k) = T_2d(i,k)h_lay(k) ; S0xdz(k) = S_2d(i,k)h_lay(k)
308		enddo
309		else
310	0	do k=1,nz
311	0	T0xdz(k) = rho_2d(i,k)h_lay(k) ; S0xdz(k) = rho_2d(i,k)h_lay(k)
312		enddo
313		endif
314	0	nzc = nz
315	0	do k=1,nzc+1 ; kc(k) = k ; kf(k) = 0.0 ; enddo
316		endif
317		f2 = 0.25 * ((G%Coriolis2Bu(I,J) + G%Coriolis2Bu(I-1,J-1)) + &
318	46250	(G%Coriolis2Bu(I,J-1) + G%Coriolis2Bu(I-1,J)))
319	46250	surface_pres = 0.0 ; if (associated(p_surf)) surface_pres = p_surf(i,j)
320
321		! ---------------------------------------------------- I_Ld2_1d, dz_Int_1d
322
323		! Set the initial guess for kappa, here defined at interfaces.
324		! ----------------------------------------------------
325	484143	do K=1,nzc+1 ; kappa(K) = CS%kappa_seed ; enddo
326
327		call kappa_shear_column(kappa, tke, dt, nzc, f2, surface_pres, &
328		h_lay, dz_lay, u0xdz, v0xdz, T0xdz, S0xdz, kappa_avg, &
329		tke_avg, N2_init, S2_init, N2_mean, S2_mean, &
330	46250	tv, CS, GV, US)
331
332		! call cpu_clock_begin(id_clock_setup)
333		! Extrapolate from the vertically reduced grid back to the original layers.
334	46250	if (nz == nzc) then
335	0	do K=1,nz+1
336	0	kappa_2d(i,K) = kappa_avg(K)
337	0	if (CS%all_layer_TKE_bug) then
338	0	tke_2d(i,K) = tke(K)
339		else
340	0	tke_2d(i,K) = tke_avg(K)
341		endif
342		enddo
343	0	if (CS%id_N2_mean>0) then ; do K=1,nz+1
344	0	diag_N2_mean(i,j,K) = N2_mean(K)
345		enddo ; endif
346	0	if (CS%id_S2_mean>0) then ; do K=1,nz+1
347	0	diag_S2_mean(i,j,K) = S2_mean(K)
348		enddo ; endif
349	0	if ((CS%id_N2_init>0) .or. CS%debug) then ; do K=1,nz+1
350	0	diag_N2_init(i,j,K) = N2_init(K)
351		enddo ; endif
352	0	if ((CS%id_S2_init>0) .or. CS%debug) then ; do K=1,nz+1
353	0	diag_S2_init(i,j,K) = S2_init(K)
354		enddo ; endif
355		else
356	1110000	do K=1,nz+1
357	1110000	if (kf(K) == 0.0) then
358	437893	kappa_2d(i,K) = kappa_avg(kc(K))
359	437893	tke_2d(i,K) = tke_avg(kc(K))
360		else
361	625857	kappa_2d(i,K) = (1.0-kf(K)) * kappa_avg(kc(K)) + kf(K) * kappa_avg(kc(K)+1)
362	625857	tke_2d(i,K) = (1.0-kf(K)) * tke_avg(kc(K)) + kf(K) * tke_avg(kc(K)+1)
363		endif
364		enddo
365	1110000	do K=1,nz+1
366	1110000	if (kf(K) == 0.0) then
367	437893	if (CS%id_N2_mean>0) diag_N2_mean(i,j,K) = N2_mean(kc(K))
368	437893	if (CS%id_S2_mean>0) diag_S2_mean(i,j,K) = S2_mean(kc(K))
369	437893	if ((CS%id_N2_init>0) .or. CS%debug) diag_N2_init(i,j,K) = N2_init(kc(K))
370	437893	if ((CS%id_S2_init>0) .or. CS%debug) diag_S2_init(i,j,K) = S2_init(kc(K))
371		else
372	625857	if (CS%id_N2_mean>0) &
373	0	diag_N2_mean(i,j,K) = (1.0-kf(K)) * N2_mean(kc(K)) + kf(K) * N2_mean(kc(K)+1)
374	625857	if (CS%id_S2_mean>0) &
375	0	diag_S2_mean(i,j,K) = (1.0-kf(K)) * S2_mean(kc(K)) + kf(K) * S2_mean(kc(K)+1)
376	625857	if ((CS%id_N2_init>0) .or. CS%debug) &
377	0	diag_N2_init(i,j,K) = (1.0-kf(K)) * N2_init(kc(K)) + kf(K) * N2_init(kc(K)+1)
378	625857	if ((CS%id_S2_init>0) .or. CS%debug) &
379	0	diag_S2_init(i,j,K) = (1.0-kf(K)) * S2_init(kc(K)) + kf(K) * S2_init(kc(K)+1)
380		endif
381		enddo
382		endif
383		! call cpu_clock_end(id_clock_setup)
384		else ! Land points, still inside the i-loop.
385	445200	do K=1,nz+1
386	445200	kappa_2d(i,K) = 0.0 ; tke_2d(i,K) = 0.0
387		enddo
388		endif ; enddo ! i-loop
389
390	1494721	do K=1,nz+1 ; do i=is,ie
391	1490400	kappa_io(i,j,K) = G%mask2dT(i,j) * kappa_2d(i,K)
392	1490400	tke_io(i,j,K) = G%mask2dT(i,j) * tke_2d(i,K)
393	1494540	kv_io(i,j,K) = ( G%mask2dT(i,j) * kappa_2d(i,K) ) * CS%Prandtl_turb
394
395		enddo ; enddo
396
397		enddo ! end of j-loop
398
399	1	if (CS%debug) then
400	0	call hchksum(diag_N2_init, "kappa_shear N2_init", G%HI, unscale=US%s_to_T**2)
401	0	call hchksum(diag_S2_init, "kappa_shear S2_init", G%HI, unscale=US%s_to_T**2)
402	0	call hchksum(kappa_io, "kappa", G%HI, unscale=GV%HZ_T_to_m2_s)
403	0	call hchksum(tke_io, "tke", G%HI, unscale=US%Z_to_m*2US%s_to_T**2)
404		endif
405
406	1	if (CS%id_Kd_shear > 0) call post_data(CS%id_Kd_shear, kappa_io, CS%diag)
407	1	if (CS%id_TKE > 0) call post_data(CS%id_TKE, tke_io, CS%diag)
408	1	if (CS%id_N2_init > 0) call post_data(CS%id_N2_init, diag_N2_init, CS%diag)
409	1	if (CS%id_S2_init > 0) call post_data(CS%id_S2_init, diag_S2_init, CS%diag)
410	1	if (CS%id_N2_mean > 0) call post_data(CS%id_N2_mean, diag_N2_mean, CS%diag)
411	1	if (CS%id_S2_mean > 0) call post_data(CS%id_S2_mean, diag_S2_mean, CS%diag)
412
413	1	end subroutine Calculate_kappa_shear
414
415
416		!> Subroutine for calculating shear-driven diffusivity and TKE in corner columns
417	0	subroutine Calc_kappa_shear_vertex(u_in, v_in, h, T_in, S_in, tv, p_surf, kappa_io, tke_io, &
418	0	kv_io, dt, G, GV, US, CS)
419		type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure.
420		type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure.
421		type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
422		real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &
423		intent(in) :: u_in !< Initial zonal velocity [L T-1 ~> m s-1].
424		real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &
425		intent(in) :: v_in !< Initial meridional velocity [L T-1 ~> m s-1].
426		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
427		intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2].
428		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
429		intent(in) :: T_in !< Layer potential temperatures [C ~> degC]
430		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
431		intent(in) :: S_in !< Layer salinities [S ~> ppt]
432		type(thermo_var_ptrs), intent(in) :: tv !< A structure containing pointers to any
433		!! available thermodynamic fields. Absent fields
434		!! have NULL ptrs.
435		real, dimension(:,:), pointer :: p_surf !< The pressure at the ocean surface [R L2 T-2 ~> Pa]
436		!! (or NULL).
437		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), &
438		intent(out) :: kappa_io !< The diapycnal diffusivity at each interface
439		!! (not layer!) [H Z T-1 ~> m2 s-1 or kg m-1 s-1].
440		real, dimension(SZIB_(G),SZJB_(G),SZK_(GV)+1), &
441		intent(out) :: tke_io !< The turbulent kinetic energy per unit mass at
442		!! each interface (not layer!) [Z2 T-2 ~> m2 s-2].
443		real, dimension(SZIB_(G),SZJB_(G),SZK_(GV)+1), &
444		intent(inout) :: kv_io !< The vertical viscosity at each interface
445		!! [H Z T-1 ~> m2 s-1 or Pa s].
446		!! The previous value is used to initialize kappa
447		!! in the vertex columns as Kappa = Kv/Prandtl
448		!! to accelerate the iteration toward convergence.
449		real, intent(in) :: dt !< Time increment [T ~> s].
450		type(Kappa_shear_CS), pointer :: CS !< The control structure returned by a previous
451		!! call to kappa_shear_init.
452
453		! Local variables
454		real, dimension(SZIB_(G),SZJB_(G),SZK_(GV)+1) :: &
455	0	diag_N2_init, & ! Diagnostic of N2 as provided to this routine [T-2 ~> s-2]
456	0	diag_S2_init, & ! Diagnostic of S2 as provided to this routine [T-2 ~> s-2]
457	0	diag_N2_mean, & ! Diagnostic of N2 averaged over the timestep applied [T-2 ~> s-2]
458	0	diag_S2_mean ! Diagnostic of S2 averaged over the timestep applied [T-2 ~> s-2]
459
460		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: &
461	0	dz_3d ! Vertical distance between interface heights [Z ~> m].
462		real, dimension(SZIB_(G),SZJB_(G),SZK_(GV)+1) :: &
463	0	kappa_vertex ! Diffusivity at interfaces and vertices [H Z T-1 ~> m2 s-1 or Pa s]
464		real, dimension(SZIB_(G),SZJB_(G),SZK_(GV)) :: &
465	0	h_vert ! Thicknesses interpolated to vertices [H ~> m or kg m-2]
466		real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)) :: &
467	0	h_at_u ! A mask-weighted thickness interpolated to u-points [H ~> m or kg m-2]
468		real, dimension(SZI_(G),SZJB_(G),SZK_(GV)) :: &
469	0	h_at_v ! A mask-weighted thickness interpolated to v-points [H ~> m or kg m-2]
470		real, dimension(SZIB_(G),SZK_(GV)) :: &
471	0	h_2d, & ! A 2-D version of h interpolated to vertices [H ~> m or kg m-2].
472	0	dz_2d, & ! Vertical distance between interface heights [Z ~> m].
473	0	u_2d, v_2d, & ! 2-D versions of u_in and v_in, converted to [L T-1 ~> m s-1].
474	0	T_2d, S_2d, rho_2d ! 2-D versions of T [C ~> degC], S [S ~> ppt], and rho [R ~> kg m-3].
475		real, dimension(SZIB_(G),SZK_(GV)+1) :: &
476	0	kappa_2d ! 2-D slice of kappa_vert [H Z T-1 ~> m2 s-1 or Pa s]
477		real, dimension(SZIB_(G),SZK_(GV)+1) :: &
478	0	tke_2d ! 2-D version tke_io [Z2 T-2 ~> m2 s-2].
479		real, dimension(SZK_(GV)) :: &
480	0	Idz, & ! The inverse of the thickness of the merged layers [H-1 ~> m2 kg-1].
481	0	h_lay, & ! The layer thickness [H ~> m or kg m-2]
482	0	dz_lay, & ! The geometric layer thickness in height units [Z ~> m]
483	0	u0xdz, & ! The initial zonal velocity times dz [L H T-1 ~> m2 s-1 or kg m-1 s-1].
484	0	v0xdz, & ! The initial meridional velocity times dz [H L T-1 ~> m2 s-1 or kg m-1 s-1]
485	0	T0xdz, & ! The initial temperature times dz [C H ~> degC m or degC kg m-2]
486	0	S0xdz ! The initial salinity times dz [S H ~> ppt m or ppt kg m-2]
487		real, dimension(SZK_(GV)+1) :: &
488	0	kappa, & ! The shear-driven diapycnal diffusivity at an interface [H Z T-1 ~> m2 s-1 or Pa s]
489	0	tke, & ! The Turbulent Kinetic Energy per unit mass at an interface [Z2 T-2 ~> m2 s-2].
490	0	kappa_avg, & ! The time-weighted average of kappa [H Z T-1 ~> m2 s-1 or Pa s]
491	0	tke_avg, & ! The time-weighted average of TKE [Z2 T-2 ~> m2 s-2].
492	0	N2_init, & ! N2 as provided to this routine [T-2 ~> s-2].
493	0	S2_init, & ! S2 as provided to this routine [T-2 ~> s-2].
494	0	N2_mean, & ! The time-weighted average of N2 [T-2 ~> s-2].
495	0	S2_mean ! The time-weighted average of S2 [T-2 ~> s-2].
496
497		real :: f2 ! The squared Coriolis parameter of each column [T-2 ~> s-2].
498		real :: surface_pres ! The top surface pressure [R L2 T-2 ~> Pa].
499
500		real :: dz_in_lay ! The running sum of the thickness in a layer [H ~> m or kg m-2]
501		real :: k0dt ! The background diffusivity times the timestep [H Z ~> m2 or kg m-1]
502		real :: dz_massless ! A layer thickness that is considered massless [H ~> m or kg m-2]
503		real :: I_hwt ! The inverse of the sum of the adjacent masked thickness weights [H-1 ~> m-1 or m2 kg-1]
504		real :: I_htot ! The inverse of the sum of the thicknesses at adjacent vertices [H-1 ~> m-1 or m2 kg-1]
505		real :: I_Prandtl ! The inverse of the turbulent Prandtl number [nondim].
506		logical :: use_temperature ! If true, temperature and salinity have been
507		! allocated and are being used as state variables.
508
509	0	integer, dimension(SZK_(GV)+1) :: kc ! The index map between the original
510		! interfaces and the interfaces with massless layers
511		! merged into nearby massive layers.
512	0	real, dimension(SZK_(GV)+1) :: kf ! The fractional weight of interface kc+1 for
513		! interpolating back to the original index space [nondim].
514		real :: h_SW, h_SE, h_NW, h_NE ! Thicknesses at adjacent vertices [H ~> m or kg m-2]
515		real :: mks_to_HZ_T ! A factor used to restore dimensional scaling after the geometric mean
516		! diffusivity is taken using thickness weighted powers [H Z s m-2 T-1 ~> 1]
517		! or [H Z m s kg-1 T-1 ~> 1]
518		real :: H_tiny ! A sub-roundoff thickness to use in the denominator when calculating
519		! thickness-weighted averages [H ~> m or kg m-2]
520		integer :: IsB, IeB, JsB, JeB, i, j, k, nz, nzc
521
522		! Diagnostics that should be deleted?
523	0	isB = G%isc-1 ; ieB = G%iecB ; jsB = G%jsc-1 ; jeB = G%jecB ; nz = GV%ke
524
525	0	if ((CS%id_N2_init>0) .or. CS%debug) diag_N2_init(:,:,:) = 0.0
526	0	if ((CS%id_S2_init>0) .or. CS%debug) diag_S2_init(:,:,:) = 0.0
527	0	if (CS%id_N2_mean>0) diag_N2_mean(:,:,:) = 0.0
528	0	if (CS%id_S2_mean>0) diag_S2_mean(:,:,:) = 0.0
529	0	kappa_vertex(:,:,:) = 0.0
530
531	0	use_temperature = associated(tv%T)
532
533	0	k0dt = dt*CS%kappa_0
534	0	dz_massless = 0.1sqrt((US%Z_to_mGV%m_to_H)*k0dt)
535	0	I_Prandtl = 0.0 ; if (CS%Prandtl_turb > 0.0) I_Prandtl = 1.0 / CS%Prandtl_turb
536	0	H_tiny = 0.5 * GV%H_subroundoff
537
538		! Convert layer thicknesses into geometric thickness in height units.
539	0	call thickness_to_dz(h, tv, dz_3d, G, GV, US, halo_size=1)
540
541	0	if (CS%vertex_shear_OBC_bug) then
542		!$OMP parallel do default(shared)
543	0	do k=1,nz
544	0	do j=JsB,JeB+1 ; do I=IsB,IeB
545	0	h_at_u(I,j,k) = G%mask2dCu(I,j) * (h(i,j,k) + h(i+1,j,k)) * 0.5
546		enddo ; enddo
547	0	do J=JsB,JeB ; do i=IsB,IeB+1
548	0	h_at_v(i,J,k) = G%mask2dCv(i,J) * (h(i,j,k) + h(i,j+1,k)) * 0.5
549		enddo ; enddo
550		enddo
551		else
552		! Because G%mask2dCu(I,j) is zero if either G%mask2dT(i,j) or G%mask2dT(i+1,j) except at OBC
553		! faces, the following form give equivalent answers to those above unless OBCs are in use,
554		! although the former is clearly less complicated and costly.
555		!$OMP parallel do default(shared)
556	0	do k=1,nz
557	0	do j=JsB,JeB+1 ; do I=IsB,IeB
558		h_at_u(I,j,k) = G%mask2dCu(I,j) * (G%mask2dT(i,j) * h(i,j,k) + G%mask2dT(i+1,j) * h(i+1,j,k)) / &
559	0	(G%mask2dT(i,j) + G%mask2dT(i+1,j) + 1.0e-36)
560		enddo ; enddo
561	0	do J=JsB,JeB ; do i=IsB,IeB+1
562		h_at_v(i,J,k) = G%mask2dCv(i,J) * (G%mask2dT(i,j) * h(i,j,k) + G%mask2dT(i,j+1) * h(i,j+1,k)) / &
563	0	(G%mask2dT(i,j) + G%mask2dT(i,j+1) + 1.0e-36)
564		enddo ; enddo
565		enddo
566		endif
567
568
569		!$OMP parallel do default(private) shared(jsB,jeB,isB,ieB,nz,h,u_in,v_in,T_in,S_in,h_at_u,h_at_v,dz_3d,H_tiny, &
570		!$OMP use_temperature,tv,G,GV,US,CS,kappa_io, &
571		!$OMP dz_massless,k0dt,p_surf,dt,tke_io,kv_io,kappa_vertex,h_vert,I_Prandtl, &
572		!$OMP diag_N2_init,diag_S2_init,diag_N2_mean,diag_S2_mean)
573	0	do J=JsB,JeB
574
575		! Interpolate the various quantities to the corners, using masks.
576	0	do k=1,nz ; do I=IsB,IeB
577		u_2d(I,k) = ( (u_in(I,j,k) * h_at_u(I,j,k)) + (u_in(I,j+1,k) * h_at_u(I,j+1,k)) ) / &
578	0	( (h_at_u(I,j,k) + h_at_u(I,j+1,k)) + H_tiny )
579		v_2d(I,k) = ( (v_in(i,J,k) * h_at_v(i,J,k)) + (v_in(i+1,J,k) * h_at_v(i+1,J,k)) ) / &
580	0	( (h_at_v(i,J,k) + h_at_v(i+1,J,k)) + H_tiny )
581
582		I_hwt = 1.0 / (((G%mask2dT(i,j) * h(i,j,k) + G%mask2dT(i+1,j+1) * h(i+1,j+1,k)) + &
583		(G%mask2dT(i+1,j) * h(i+1,j,k) + G%mask2dT(i,j+1) * h(i,j+1,k))) + &
584	0	GV%H_subroundoff)
585	0	if (use_temperature) then
586		T_2d(I,k) = ( (G%mask2dT(i,j) * (h(i,j,k) * T_in(i,j,k)) + &
587		G%mask2dT(i+1,j+1) * (h(i+1,j+1,k) * T_in(i+1,j+1,k))) + &
588		(G%mask2dT(i+1,j) * (h(i+1,j,k) * T_in(i+1,j,k)) + &
589	0	G%mask2dT(i,j+1) * (h(i,j+1,k) * T_in(i,j+1,k))) ) * I_hwt
590		S_2d(I,k) = ( (G%mask2dT(i,j) * (h(i,j,k) * S_in(i,j,k)) + &
591		G%mask2dT(i+1,j+1) * (h(i+1,j+1,k) * S_in(i+1,j+1,k))) + &
592		(G%mask2dT(i+1,j) * (h(i+1,j,k) * S_in(i+1,j,k)) + &
593	0	G%mask2dT(i,j+1) * (h(i,j+1,k) * S_in(i,j+1,k))) ) * I_hwt
594		endif
595		h_2d(I,k) = ((G%mask2dT(i,j) * h(i,j,k) + G%mask2dT(i+1,j+1) * h(i+1,j+1,k)) + &
596		(G%mask2dT(i+1,j) * h(i+1,j,k) + G%mask2dT(i,j+1) * h(i,j+1,k)) ) / &
597		((G%mask2dT(i,j) + G%mask2dT(i+1,j+1)) + &
598	0	(G%mask2dT(i+1,j) + G%mask2dT(i,j+1)) + 1.0e-36 )
599		dz_2d(I,k) = ((G%mask2dT(i,j) * dz_3d(i,j,k) + G%mask2dT(i+1,j+1) * dz_3d(i+1,j+1,k)) + &
600		(G%mask2dT(i+1,j) * dz_3d(i+1,j,k) + G%mask2dT(i,j+1) * dz_3d(i,j+1,k)) ) / &
601		((G%mask2dT(i,j) + G%mask2dT(i+1,j+1)) + &
602	0	(G%mask2dT(i+1,j) + G%mask2dT(i,j+1)) + 1.0e-36 )
603		! h_2d(I,k) = 0.25*((h(i,j,k) + h(i+1,j+1,k)) + (h(i+1,j,k) + h(i,j+1,k)))
604		! h_2d(I,k) = (((h(i,j,k)2) + (h(i+1,j+1,k)2)) + &
605		! ((h(i+1,j,k)2) + (h(i,j+1,k)2))) * I_hwt
606		enddo ; enddo
607	0	if (.not.use_temperature) then ; do k=1,nz ; do I=IsB,IeB
608	0	rho_2d(I,k) = GV%Rlay(k)
609		enddo ; enddo ; endif
610
611		!---------------------------------------
612		! Work on each column.
613		!---------------------------------------
614	0	do I=IsB,IeB ; if ((G%mask2dCu(I,j) + G%mask2dCu(I,j+1)) + &
615	0	(G%mask2dCv(i,J) + G%mask2dCv(i+1,J)) > 0.0) then
616		! call cpu_clock_begin(Id_clock_setup)
617		! Store a transposed version of the initial arrays.
618		! Any elimination of massless layers would occur here.
619	0	if (CS%eliminate_massless) then
620	0	nzc = 1
621	0	do k=1,nz
622		! Zero out the thicknesses of all layers, even if they are unused.
623	0	h_lay(k) = 0.0 ; dz_lay(k) = 0.0 ; u0xdz(k) = 0.0 ; v0xdz(k) = 0.0
624	0	T0xdz(k) = 0.0 ; S0xdz(k) = 0.0
625
626		! Add a new layer if this one has mass.
627		! if ((h_lay(nzc) > 0.0) .and. (h_2d(I,k) > dz_massless)) nzc = nzc+1
628	0	if ((k>CS%nkml) .and. (h_lay(nzc) > 0.0) .and. &
629	0	(h_2d(I,k) > dz_massless)) nzc = nzc+1
630
631		! Only merge clusters of massless layers.
632		! if ((h_lay(nzc) > dz_massless) .or. &
633		! ((h_lay(nzc) > 0.0) .and. (h_2d(I,k) > dz_massless))) nzc = nzc+1
634
635	0	kc(k) = nzc
636	0	h_lay(nzc) = h_lay(nzc) + h_2d(I,k)
637	0	dz_lay(nzc) = dz_lay(nzc) + dz_2d(I,k)
638	0	u0xdz(nzc) = u0xdz(nzc) + u_2d(I,k)*h_2d(I,k)
639	0	v0xdz(nzc) = v0xdz(nzc) + v_2d(I,k)*h_2d(I,k)
640	0	if (use_temperature) then
641	0	T0xdz(nzc) = T0xdz(nzc) + T_2d(I,k)*h_2d(I,k)
642	0	S0xdz(nzc) = S0xdz(nzc) + S_2d(I,k)*h_2d(I,k)
643		else
644	0	T0xdz(nzc) = T0xdz(nzc) + rho_2d(I,k)*h_2d(I,k)
645	0	S0xdz(nzc) = S0xdz(nzc) + rho_2d(I,k)*h_2d(I,k)
646		endif
647		enddo
648	0	kc(nz+1) = nzc+1
649
650		! Set up Idz as the inverse of layer thicknesses.
651	0	do k=1,nzc ; Idz(k) = 1.0 / h_lay(k) ; enddo
652
653		! Now determine kf, the fractional weight of interface kc when
654		! interpolating between interfaces kc and kc+1.
655	0	kf(1) = 0.0 ; dz_in_lay = h_2d(I,1)
656	0	do k=2,nz
657	0	if (kc(k) > kc(k-1)) then
658	0	kf(k) = 0.0 ; dz_in_lay = h_2d(I,k)
659		else
660	0	kf(k) = dz_in_lay*Idz(kc(k)) ; dz_in_lay = dz_in_lay + h_2d(I,k)
661		endif
662		enddo
663	0	kf(nz+1) = 0.0
664		else
665	0	do k=1,nz
666	0	h_lay(k) = h_2d(I,k)
667	0	dz_lay(k) = dz_2d(I,k)
668	0	u0xdz(k) = u_2d(I,k)h_lay(k) ; v0xdz(k) = v_2d(I,k)h_lay(k)
669		enddo
670	0	if (use_temperature) then
671	0	do k=1,nz
672	0	T0xdz(k) = T_2d(I,k)h_lay(k) ; S0xdz(k) = S_2d(I,k)h_lay(k)
673		enddo
674		else
675	0	do k=1,nz
676	0	T0xdz(k) = rho_2d(I,k)h_lay(k) ; S0xdz(k) = rho_2d(I,k)h_lay(k)
677		enddo
678		endif
679	0	nzc = nz
680	0	do k=1,nzc+1 ; kc(k) = k ; kf(k) = 0.0 ; enddo
681		endif
682
683	0	f2 = G%Coriolis2Bu(I,J)
684	0	surface_pres = 0.0
685	0	if (associated(p_surf)) then
686	0	if (CS%psurf_bug) then
687		! This is wrong because it is averaging values from land in some places.
688		surface_pres = 0.25 * ((p_surf(i,j) + p_surf(i+1,j+1)) + &
689	0	(p_surf(i+1,j) + p_surf(i,j+1)))
690		else
691		surface_pres = ((G%mask2dT(i,j) * p_surf(i,j) + G%mask2dT(i+1,j+1) * p_surf(i+1,j+1)) + &
692		(G%mask2dT(i+1,j) * p_surf(i+1,j) + G%mask2dT(i,j+1) * p_surf(i,j+1)) ) / &
693		((G%mask2dT(i,j) + G%mask2dT(i+1,j+1)) + &
694	0	(G%mask2dT(i+1,j) + G%mask2dT(i,j+1)) + 1.0e-36 )
695		endif
696		endif
697
698		! ----------------------------------------------------
699		! Set the initial guess for kappa, here defined at interfaces.
700		! ----------------------------------------------------
701	0	do K=1,nzc+1 ; kappa(K) = CS%kappa_seed ; enddo
702
703		call kappa_shear_column(kappa, tke, dt, nzc, f2, surface_pres, &
704		h_lay, dz_lay, u0xdz, v0xdz, T0xdz, S0xdz, kappa_avg, &
705	0	tke_avg, N2_init, S2_init, N2_mean, S2_mean, tv, CS, GV, US)
706		! call cpu_clock_begin(Id_clock_setup)
707		! Extrapolate from the vertically reduced grid back to the original layers.
708	0	if (nz == nzc) then
709	0	do K=1,nz+1
710	0	kappa_2d(I,K) = kappa_avg(K)
711	0	if (CS%all_layer_TKE_bug) then
712	0	tke_2d(I,K) = tke(K)
713		else
714	0	tke_2d(I,K) = tke_avg(K)
715		endif
716		enddo
717	0	if (CS%id_N2_mean>0) then ; do K=1,nz+1
718	0	diag_N2_mean(I,J,K) = N2_mean(K)
719		enddo ; endif
720	0	if (CS%id_S2_mean>0) then ; do K=1,nz+1
721	0	diag_S2_mean(I,J,K) = S2_mean(K)
722		enddo ; endif
723	0	if ((CS%id_N2_init>0) .or. CS%debug) then ; do K=1,nz+1
724	0	diag_N2_init(I,J,K) = N2_init(K)
725		enddo ; endif
726	0	if ((CS%id_S2_init>0) .or. CS%debug) then ; do K=1,nz+1
727	0	diag_S2_init(I,J,K) = S2_init(K)
728		enddo ; endif
729		else
730	0	do K=1,nz+1
731	0	if (kf(K) == 0.0) then
732	0	kappa_2d(I,K) = kappa_avg(kc(K))
733	0	tke_2d(I,K) = tke_avg(kc(K))
734		else
735	0	kappa_2d(I,K) = (1.0-kf(K)) * kappa_avg(kc(K)) + kf(K) * kappa_avg(kc(K)+1)
736	0	tke_2d(I,K) = (1.0-kf(K)) * tke_avg(kc(K)) + kf(K) * tke_avg(kc(K)+1)
737		endif
738		enddo
739	0	do K=1,nz+1
740	0	if (kf(K) == 0.0) then
741	0	if (CS%id_N2_mean>0) diag_N2_mean(I,J,K) = N2_mean(kc(K))
742	0	if (CS%id_S2_mean>0) diag_S2_mean(I,J,K) = S2_mean(kc(K))
743	0	if ((CS%id_N2_init>0) .or. CS%debug) diag_N2_init(I,J,K) = N2_init(kc(K))
744	0	if ((CS%id_S2_init>0) .or. CS%debug) diag_S2_init(I,J,K) = S2_init(kc(K))
745		else
746	0	if (CS%id_N2_mean>0) &
747	0	diag_N2_mean(I,J,K) = (1.0-kf(K)) * N2_mean(kc(K)) + kf(K) * N2_mean(kc(K)+1)
748	0	if (CS%id_S2_mean>0) &
749	0	diag_S2_mean(I,J,K) = (1.0-kf(K)) * S2_mean(kc(K)) + kf(K) * S2_mean(kc(K)+1)
750	0	if ((CS%id_N2_init>0) .or. CS%debug) &
751	0	diag_N2_init(I,J,K) = (1.0-kf(K)) * N2_init(kc(K)) + kf(K) * N2_init(kc(K)+1)
752	0	if ((CS%id_S2_init>0) .or. CS%debug) &
753	0	diag_S2_init(I,J,K) = (1.0-kf(K)) * S2_init(kc(K)) + kf(K) * S2_init(kc(K)+1)
754		endif
755		enddo
756		endif
757		! call cpu_clock_end(Id_clock_setup)
758		else ! Land points, still inside the i-loop.
759	0	do K=1,nz+1
760	0	kappa_2d(I,K) = 0.0 ; tke_2d(I,K) = 0.0
761		enddo
762		endif ; enddo ! i-loop
763
764		! Store the 2-d slices back in the 3-d arrays for restarts or interpolation back to tracer points.
765	0	if (CS%VS_ThicknessMean) then
766	0	do K=1,nz+1 ; do I=IsB,IeB
767	0	h_vert(I,J,k) = h_2d(I,k)
768		enddo ; enddo
769		endif
770	0	if (CS%VS_viscosity_bug) then
771	0	do K=1,nz+1 ; do I=IsB,IeB
772	0	kappa_vertex(I,J,K) = kappa_2d(I,K)
773	0	tke_io(I,J,K) = G%mask2dBu(I,J) * tke_2d(I,K)
774	0	kv_io(I,J,K) = ( G%mask2dBu(I,J) * kappa_vertex(I,J,K) ) * CS%Prandtl_turb
775		enddo ; enddo
776		else
777	0	do K=1,nz+1 ; do I=IsB,IeB
778	0	kappa_vertex(I,J,K) = kappa_2d(I,K)
779	0	tke_io(I,J,K) = tke_2d(I,K)
780	0	kv_io(I,J,K) = kappa_vertex(I,J,K) * CS%Prandtl_turb
781		enddo ; enddo
782		endif
783		enddo ! end of J-loop
784
785		! Set the diffusivities in tracer columns from the values at vertices.
786
787		!$OMP parallel do default(private) shared(G,kappa_io)
788	0	do j=G%jsc,G%jec ; do i=G%isc,G%iec
789		! The turbulent length scales (and hence turbulent diffusivity) should always go to 0 at the top and bottom.
790	0	kappa_io(i,j,1) = 0.0
791	0	kappa_io(i,j,nz+1) = 0.0
792		enddo ; enddo
793	0	if (CS%VS_ThicknessMean .and. CS%VS_GeometricMean) then
794		! This conversion factor is required to allow for arbitrary fractional powers of the diffusivities.
795	0	mks_to_HZ_T = 1.0 / GV%HZ_T_to_MKS
796		!$OMP parallel do default(private) shared(nz,G,GV,CS,kappa_io,kappa_vertex,h_vert)
797	0	do K=2,nz ; do j=G%jsc,G%jec ; do i=G%isc,G%iec
798	0	h_SW = 0.5 * (h_vert(I-1,J-1,k) + h_vert(I-1,J-1,k-1))
799	0	h_NE = 0.5 * (h_vert(I,J,k) + h_vert(I,J,k-1))
800	0	h_NW = 0.5 * (h_vert(I-1,J,k) + h_vert(I-1,J,k-1))
801	0	h_SE = 0.5 * (h_vert(I,J-1,k) + h_vert(I,J-1,k-1))
802	0	if ((h_SW + h_NE) + (h_NW + h_SE) > 0.0) then
803		! The geometric mean is zero if any component is zero, hence the need to use a floor
804		! on the value of kappa_trunc in regions on boundaries of shear zones.
805	0	I_htot = 1.0 / ((h_SW + h_NE) + (h_NW + h_SE))
806		kappa_io(i,j,K) = G%mask2dT(i,j) * mks_to_HZ_T * &
807		( ((GV%HZ_T_to_MKS * max(kappa_vertex(I-1,J-1,K), CS%VS_GeoMean_Kdmin))*(h_SWI_htot) * &
808		(GV%HZ_T_to_MKS * max(kappa_vertex(I,J,K), CS%VS_GeoMean_Kdmin))*(h_NEI_htot)) * &
809		((GV%HZ_T_to_MKS * max(kappa_vertex(I-1,J,K), CS%VS_GeoMean_Kdmin))*(h_NWI_htot) * &
810	0	(GV%HZ_T_to_MKS * max(kappa_vertex(I,J-1,K), CS%VS_GeoMean_Kdmin))*(h_SEI_htot)) )
811		else
812		! If all points have zero thickness, the thickness-weighted geometric mean is undefined, so use
813		! the non-thickness weighted geometric mean instead.
814		kappa_io(i,j,K) = G%mask2dT(i,j) * sqrt(sqrt( &
815		(max(kappa_vertex(I-1,J-1,K),CS%VS_GeoMean_Kdmin) * max(kappa_vertex(I,J,K),CS%VS_GeoMean_Kdmin)) * &
816	0	(max(kappa_vertex(I-1,J,K),CS%VS_GeoMean_Kdmin) * max(kappa_vertex(I,J-1,K),CS%VS_GeoMean_Kdmin)) ))
817		endif
818		enddo ; enddo ; enddo
819	0	elseif (CS%VS_ThicknessMean) then ! Use thickness-weighted arithmetic mean diffusivities.
820		!$OMP parallel do default(private) shared(nz,G,GV,CS,kappa_io,kappa_vertex,h_vert)
821	0	do K=2,nz ; do j=G%jsc,G%jec ; do i=G%isc,G%iec
822	0	h_SW = 0.5 * (h_vert(I-1,J-1,k) + h_vert(I-1,J-1,k-1))
823	0	h_NE = 0.5 * (h_vert(I,J,k) + h_vert(I,J,k-1))
824	0	h_NW = 0.5 * (h_vert(I-1,J,k) + h_vert(I-1,J,k-1))
825	0	h_SE = 0.5 * (h_vert(I,J-1,k) + h_vert(I,J-1,k-1))
826		! The following expression is a thickness weighted arithmetic mean at tracer points:
827	0	I_htot = 1.0 / (((h_SW + h_NE) + (h_NW + h_SE)) + GV%H_subroundoff)
828		kappa_io(i,j,K) = G%mask2dT(i,j) * &
829		(((kappa_vertex(I-1,J-1,K)h_SW) + (kappa_vertex(I,J,K)h_NE)) + &
830	0	((kappa_vertex(I-1,J,K)h_NW) + (kappa_vertex(I,J-1,K)h_SE))) * I_htot
831		enddo ; enddo ; enddo
832	0	elseif (CS%VS_GeometricMean) then ! The geometic mean diffusivities are not thickness weighted.
833		!$OMP parallel do default(private) shared(nz,G,CS,kappa_io,kappa_vertex)
834	0	do K=2,nz ; do j=G%jsc,G%jec ; do i=G%isc,G%iec
835		kappa_io(i,j,K) = G%mask2dT(i,j) * sqrt(sqrt( &
836		(max(kappa_vertex(I-1,J-1,K),CS%VS_GeoMean_Kdmin) * max(kappa_vertex(I,J,K),CS%VS_GeoMean_Kdmin)) * &
837	0	(max(kappa_vertex(I-1,J,K),CS%VS_GeoMean_Kdmin) * max(kappa_vertex(I,J-1,K),CS%VS_GeoMean_Kdmin)) ))
838		enddo ; enddo ; enddo
839		else ! Use a non-thickness weighted arithmetic mean.
840		!$OMP parallel do default(private) shared(nz,G,CS,kappa_io,kappa_vertex)
841	0	do K=2,nz ; do j=G%jsc,G%jec ; do i=G%isc,G%iec
842		kappa_io(i,j,K) = G%mask2dT(i,j) * 0.25 * &
843		((kappa_vertex(I-1,J-1,K) + kappa_vertex(I,J,K)) +&
844	0	(kappa_vertex(I-1,J,K) + kappa_vertex(I,J-1,K)))
845		enddo ; enddo ; enddo
846		endif
847
848	0	if (CS%debug) then
849	0	call Bchksum(diag_N2_init, "shear_vertex N2_init", G%HI, unscale=US%s_to_T**2)
850	0	call Bchksum(diag_S2_init, "shear_vertex S2_init", G%HI, unscale=US%s_to_T**2)
851	0	call hchksum(kappa_io, "kappa", G%HI, unscale=GV%HZ_T_to_m2_s)
852	0	call Bchksum(tke_io, "tke", G%HI, unscale=US%Z_to_m*2US%s_to_T**2)
853		endif
854
855	0	if (CS%id_Kd_shear > 0) call post_data(CS%id_Kd_shear, kappa_io, CS%diag)
856	0	if (CS%id_TKE > 0) call post_data(CS%id_TKE, tke_io, CS%diag)
857	0	if (CS%id_Kd_vertex > 0) call post_data(CS%id_Kd_vertex, kappa_vertex, CS%diag)
858	0	if (CS%id_N2_init > 0) call post_data(CS%id_N2_init, diag_N2_init, CS%diag)
859	0	if (CS%id_S2_init > 0) call post_data(CS%id_S2_init, diag_S2_init, CS%diag)
860	0	if (CS%id_N2_mean > 0) call post_data(CS%id_N2_mean, diag_N2_mean, CS%diag)
861	0	if (CS%id_S2_mean > 0) call post_data(CS%id_S2_mean, diag_S2_mean, CS%diag)
862
863	0	end subroutine Calc_kappa_shear_vertex
864
865
866		!> This subroutine calculates shear-driven diffusivity and TKE in a single column
867	46250	subroutine kappa_shear_column(kappa, tke, dt, nzc, f2, surface_pres, hlay, dz_lay, &
868	46250	u0xdz, v0xdz, T0xdz, S0xdz, kappa_avg, tke_avg, N2_init, S2_init, &
869	46250	N2_mean, S2_mean, tv, CS, GV, US )
870		type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure.
871		real, dimension(SZK_(GV)+1), &
872		intent(inout) :: kappa !< The time-weighted average of kappa [H Z T-1 ~> m2 s-1 or Pa s]
873		real, dimension(SZK_(GV)+1), &
874		intent(out) :: tke !< The Turbulent Kinetic Energy per unit mass at
875		!! an interface [Z2 T-2 ~> m2 s-2].
876		integer, intent(in) :: nzc !< The number of active layers in the column.
877		real, intent(in) :: f2 !< The square of the Coriolis parameter [T-2 ~> s-2].
878		real, intent(in) :: surface_pres !< The surface pressure [R L2 T-2 ~> Pa].
879		real, dimension(SZK_(GV)), &
880		intent(in) :: hlay !< The layer thickness [H ~> m or kg m-2]
881		real, dimension(SZK_(GV)), &
882		intent(in) :: dz_lay !< The geometric layer thickness in height units [Z ~> m]
883		real, dimension(SZK_(GV)), &
884		intent(in) :: u0xdz !< The initial zonal velocity times hlay [H L T-1 ~> m2 s-1 or kg m-1 s-1]
885		real, dimension(SZK_(GV)), &
886		intent(in) :: v0xdz !< The initial meridional velocity times the
887		!! layer thickness [H L T-1 ~> m2 s-1 or kg m-1 s-1]
888		real, dimension(SZK_(GV)), &
889		intent(in) :: T0xdz !< The initial temperature times hlay [C H ~> degC m or degC kg m-2]
890		real, dimension(SZK_(GV)), &
891		intent(in) :: S0xdz !< The initial salinity times hlay [S H ~> ppt m or ppt kg m-2]
892		real, dimension(SZK_(GV)+1), &
893		intent(out) :: kappa_avg !< The time-weighted average of kappa [H Z T-1 ~> m2 s-1 or Pa s]
894		real, dimension(SZK_(GV)+1), &
895		intent(out) :: tke_avg !< The time-weighted average of TKE [Z2 T-2 ~> m2 s-2].
896		real, dimension(SZK_(GV)+1), &
897		intent(out) :: N2_mean !< The time-weighted average of N2 [Z2 T-2 ~> m2 s-2].
898		real, dimension(SZK_(GV)+1), &
899		intent(out) :: S2_mean !< The time-weighted average of S2 [Z2 T-2 ~> m2 s-2].
900		real, dimension(SZK_(GV)+1), &
901		intent(out) :: N2_init !< The initial value of N2 [Z2 T-2 ~> m2 s-2].
902		real, dimension(SZK_(GV)+1), &
903		intent(out) :: S2_init !< The initial value of S2 [Z2 T-2 ~> m2 s-2].
904		real, intent(in) :: dt !< Time increment [T ~> s].
905		type(thermo_var_ptrs), intent(in) :: tv !< A structure containing pointers to any
906		!! available thermodynamic fields. Absent fields
907		!! have NULL ptrs.
908		type(Kappa_shear_CS), pointer :: CS !< The control structure returned by a previous
909		!! call to kappa_shear_init.
910		type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
911
912		! Local variables
913		real, dimension(nzc) :: &
914	92500	u, & ! The zonal velocity after a timestep of mixing [L T-1 ~> m s-1].
915	92500	v, & ! The meridional velocity after a timestep of mixing [L T-1 ~> m s-1].
916	92500	Idz, & ! The inverse of the distance between TKE points [Z-1 ~> m-1].
917	92500	T, & ! The potential temperature after a timestep of mixing [C ~> degC].
918	92500	Sal, & ! The salinity after a timestep of mixing [S ~> ppt].
919	92500	u_test, v_test, & ! Temporary velocities [L T-1 ~> m s-1].
920	92500	T_test, S_test ! Temporary temperatures [C ~> degC] and salinities [S ~> ppt].
921
922		real, dimension(nzc+1) :: &
923	92500	N2, & ! The squared buoyancy frequency at an interface [T-2 ~> s-2].
924	92500	h_Int, & ! The extent of a finite-volume space surrounding an interface,
925		! as used in calculating kappa and TKE [H ~> m or kg m-2]
926	92500	dz_Int, & ! The vertical distance with the space surrounding an interface,
927		! as used in calculating kappa and TKE [Z ~> m]
928	92500	dz_h_Int, & ! The ratio of the vertical distances to the thickness around an
929		! interface [Z H-1 ~> nondim or m3 kg-1]. In non-Boussinesq mode
930		! this is the specific volume, otherwise it is a scaling factor.
931	92500	I_dz_int, & ! The inverse of the distance between velocity & density points
932		! above and below an interface [Z-1 ~> m-1]. This is used to
933		! calculate N2, shear and fluxes.
934	92500	S2, & ! The squared shear at an interface [T-2 ~> s-2].
935	92500	a1, & ! a1 is the coupling between adjacent interfaces in the TKE,
936		! velocity, and density equations [H ~> m or kg m-2]
937	92500	c1, & ! c1 is used in the tridiagonal (and similar) solvers [nondim].
938	92500	k_src, & ! The shear-dependent source term in the kappa equation [T-1 ~> s-1]
939	92500	kappa_src, & ! The shear-dependent source term in the kappa equation [T-1 ~> s-1]
940	92500	kappa_out, & ! The kappa that results from the kappa equation [H Z T-1 ~> m2 s-1 or Pa s]
941	92500	kappa_mid, & ! The average of the initial and predictor estimates of kappa [H Z T-1 ~> m2 s-1 or Pa s]
942	92500	tke_pred, & ! The value of TKE from a predictor step [Z2 T-2 ~> m2 s-2].
943	92500	kappa_pred, & ! The value of kappa from a predictor step [H Z T-1 ~> m2 s-1 or Pa s]
944	92500	pressure, & ! The pressure at an interface [R L2 T-2 ~> Pa].
945	92500	T_int, & ! The temperature interpolated to an interface [C ~> degC].
946	92500	Sal_int, & ! The salinity interpolated to an interface [S ~> ppt].
947	92500	dbuoy_dT, & ! The partial derivative of buoyancy with changes in temperature [Z T-2 C-1 ~> m s-2 degC-1]
948	92500	dbuoy_dS, & ! The partial derivative of buoyancy with changes in salinity [Z T-2 S-1 ~> m s-2 ppt-1]
949	92500	dSpV_dT, & ! The partial derivative of specific volume with changes in temperature [R-1 C-1 ~> m3 kg-1 degC-1]
950	92500	dSpV_dS, & ! The partial derivative of specific volume with changes in salinity [R-1 S-1 ~> m3 kg-1 ppt-1]
951	92500	rho_int, & ! The in situ density interpolated to an interface [R ~> kg m-3]
952	92500	I_L2_bdry, & ! The inverse of the square of twice the harmonic mean
953		! distance to the top and bottom boundaries [H-1 Z-1 ~> m-2 or m kg-1].
954	92500	K_Q, & ! Diffusivity divided by TKE [H T Z-1 ~> s or kg s m-3]
955	92500	K_Q_tmp, & ! A temporary copy of diffusivity divided by TKE [H T Z-1 ~> s or kg s m-3]
956	92500	local_src_avg, & ! The time-integral of the local source [nondim]
957	92500	tol_min, & ! Minimum tolerated ksrc for the corrector step [T-1 ~> s-1].
958	92500	tol_max, & ! Maximum tolerated ksrc for the corrector step [T-1 ~> s-1].
959	92500	tol_chg, & ! The tolerated kappa change integrated over a timestep [nondim].
960	92500	dist_from_top, & ! The distance from the top surface [Z ~> m].
961	92500	h_from_top, & ! The total thickness above an interface [H ~> m or kg m-2]
962	92500	local_src ! The sum of all sources of kappa, including kappa_src and
963		! sources from the elliptic term [T-1 ~> s-1]
964
965		real :: dist_from_bot ! The distance from the bottom surface [Z ~> m].
966		real :: h_from_bot ! The total thickness below and interface [H ~> m or kg m-2]
967		real :: b1 ! The inverse of the pivot in the tridiagonal equations [H-1 ~> m-1 or m2 kg-1].
968		real :: bd1 ! A term in the denominator of b1 [H ~> m or kg m-2].
969		real :: d1 ! 1 - c1 in the tridiagonal equations [nondim]
970		real :: gR0 ! A conversion factor from H to pressure, Rho_0 times g in Boussinesq
971		! mode, or just g when non-Boussinesq [R L2 T-2 H-1 ~> kg m-2 s-2 or m s-2].
972		real :: g_R0 ! g_R0 is a rescaled version of g/Rho [Z R-1 T-2 ~> m4 kg-1 s-2].
973		real :: Norm ! A factor that normalizes two weights to 1 [H-2 ~> m-2 or m4 kg-2].
974		real :: tol_dksrc ! Tolerance for the change in the kappa source within an iteration
975		! relative to the local source [nondim]. This must be greater than 1.
976		real :: tol2 ! The tolerance for the change in the kappa source within an iteration
977		! relative to the average local source over previous iterations [nondim].
978		real :: tol_dksrc_low ! The tolerance for the fractional decrease in ksrc
979		! within an iteration [nondim]. 0 < tol_dksrc_low < 1.
980		real :: Ri_crit ! The critical shear Richardson number for shear-
981		! driven mixing [nondim]. The theoretical value is 0.25.
982		real :: dt_rem ! The remaining time to advance the solution [T ~> s].
983		real :: dt_now ! The time step used in the current iteration [T ~> s].
984		real :: dt_wt ! The fractional weight of the current iteration [nondim].
985		real :: dt_test ! A time-step that is being tested for whether it
986		! gives acceptably small changes in k_src [T ~> s].
987		real :: Idtt ! Idtt = 1 / dt_test [T-1 ~> s-1].
988		real :: dt_inc ! An increment to dt_test that is being tested [T ~> s].
989		real :: wt_a ! The fraction of a layer thickness identified with the interface
990		! above a layer [nondim]
991		real :: wt_b ! The fraction of a layer thickness identified with the interface
992		! below a layer [nondim]
993		real :: k0dt ! The background diffusivity times the timestep [H Z ~> m2 or kg m-1].
994		real :: I_lz_rescale_sqr ! The inverse of a rescaling factor for L2_bdry (Lz) squared [nondim].
995		logical :: valid_dt ! If true, all levels so far exhibit acceptably small changes in k_src.
996		logical :: use_temperature ! If true, temperature and salinity have been
997		! allocated and are being used as state variables.
998		integer :: ks_kappa, ke_kappa ! The k-range with nonzero kappas.
999		integer :: dt_refinements ! The number of 2-fold refinements that will be used
1000		! to estimate the maximum permitted time step. I.e.,
1001		! the resolution is 1/2^dt_refinements.
1002		integer :: k, itt, itt_dt
1003
1004		! This calculation of N2 is for debugging only.
1005		! real, dimension(SZK_(GV)+1) :: &
1006		! N2_debug, & ! A version of N2 for debugging [T-2 ~> s-2]
1007
1008	46250	Ri_crit = CS%Rino_crit
1009	46250	gR0 = GV%H_to_RZ * GV%g_Earth
1010	46250	g_R0 = GV%g_Earth_Z_T2 / GV%Rho0
1011	46250	k0dt = dt*CS%kappa_0
1012
1013	46250	I_lz_rescale_sqr = 1.0 ; if (CS%lz_rescale > 0) I_lz_rescale_sqr = 1/(CS%lz_rescale*CS%lz_rescale)
1014
1015	46250	tol_dksrc = CS%kappa_src_max_chg
1016	46250	if (tol_dksrc == 10.0) then
1017		! This is equivalent to the expression below, but avoids changes at roundoff for the default value.
1018	46250	tol_dksrc_low = 0.95
1019		else
1020	0	tol_dksrc_low = (tol_dksrc - 0.5)/tol_dksrc
1021		endif
1022	46250	tol2 = 2.0*CS%kappa_tol_err
1023	46250	dt_refinements = 5 ! Selected so that 1/2^dt_refinements < 1-tol_dksrc_low
1024	46250	use_temperature = .false. ; if (associated(tv%T)) use_temperature = .true.
1025
1026
1027		! Set up Idz as the inverse of layer thicknesses.
1028	437893	do k=1,nzc ; Idz(k) = 1.0 / dz_lay(k) ; enddo
1029		! Set up I_dz_int as the inverse of the distance between
1030		! adjacent layer centers.
1031	46250	I_dz_int(1) = 2.0 / dz_lay(1)
1032	46250	dist_from_top(1) = 0.0 ; h_from_top(1) = 0.0
1033	391643	do K=2,nzc
1034	345393	I_dz_int(K) = 2.0 / (dz_lay(k-1) + dz_lay(k))
1035	345393	dist_from_top(K) = dist_from_top(K-1) + dz_lay(k-1)
1036	391643	h_from_top(K) = h_from_top(K-1) + hlay(k-1)
1037		enddo
1038	46250	I_dz_int(nzc+1) = 2.0 / dz_lay(nzc)
1039
1040		! Find the inverse of the squared distances from the boundaries.
1041	46250	dist_from_bot = 0.0 ; h_from_bot = 0.0
1042	391643	do K=nzc,2,-1
1043	345393	dist_from_bot = dist_from_bot + dz_lay(k)
1044	345393	h_from_bot = h_from_bot + hlay(k)
1045		! Find the inverse of the squared distances from the boundaries,
1046		I_L2_bdry(K) = ((dist_from_top(K) + dist_from_bot) * (h_from_top(K) + h_from_bot)) / &
1047	345393	((dist_from_top(K) * dist_from_bot) * (h_from_top(K) * h_from_bot))
1048		! reduce the distance by a factor of "lz_rescale"
1049	391643	I_L2_bdry(K) = I_lz_rescale_sqr*I_L2_bdry(K)
1050		enddo
1051
1052		! Determine the velocities and thicknesses after eliminating massless
1053		! layers and applying a time-step of background diffusion.
1054	46250	if (nzc > 1) then
1055	45972	a1(2) = k0dt*I_dz_int(2)
1056	45972	b1 = 1.0 / (hlay(1) + a1(2))
1057	45972	u(1) = b1 * u0xdz(1) ; v(1) = b1 * v0xdz(1)
1058	45972	T(1) = b1 * T0xdz(1) ; Sal(1) = b1 * S0xdz(1)
1059	45972	c1(2) = a1(2) * b1 ; d1 = hlay(1) * b1 ! = 1 - c1
1060	345393	do k=2,nzc-1
1061	299421	bd1 = hlay(k) + d1*a1(k)
1062	299421	a1(k+1) = k0dt*I_dz_int(k+1)
1063	299421	b1 = 1.0 / (bd1 + a1(k+1))
1064	299421	u(k) = b1 * (u0xdz(k) + a1(k)*u(k-1))
1065	299421	v(k) = b1 * (v0xdz(k) + a1(k)*v(k-1))
1066	299421	T(k) = b1 * (T0xdz(k) + a1(k)*T(k-1))
1067	299421	Sal(k) = b1 * (S0xdz(k) + a1(k)*Sal(k-1))
1068	345393	c1(k+1) = a1(k+1) * b1 ; d1 = bd1 * b1 ! d1 = 1 - c1
1069		enddo
1070		! rho or T and S have insulating boundary conditions, u & v use no-slip
1071		! bottom boundary conditions (if kappa0 > 0).
1072		! For no-slip bottom boundary conditions
1073	45972	b1 = 1.0 / ((hlay(nzc) + d1a1(nzc)) + k0dtI_dz_int(nzc+1))
1074	45972	u(nzc) = b1 * (u0xdz(nzc) + a1(nzc)*u(nzc-1))
1075	45972	v(nzc) = b1 * (v0xdz(nzc) + a1(nzc)*v(nzc-1))
1076		! For insulating boundary conditions
1077	45972	b1 = 1.0 / (hlay(nzc) + d1*a1(nzc))
1078	45972	T(nzc) = b1 * (T0xdz(nzc) + a1(nzc)*T(nzc-1))
1079	45972	Sal(nzc) = b1 * (S0xdz(nzc) + a1(nzc)*Sal(nzc-1))
1080	391365	do k=nzc-1,1,-1
1081	345393	u(k) = u(k) + c1(k+1)u(k+1) ; v(k) = v(k) + c1(k+1)v(k+1)
1082	391365	T(k) = T(k) + c1(k+1)T(k+1) ; Sal(k) = Sal(k) + c1(k+1)Sal(k+1)
1083		enddo
1084		else
1085		! This is correct, but probably unnecessary.
1086	278	b1 = 1.0 / (hlay(1) + k0dt*I_dz_int(2))
1087	278	u(1) = b1 * u0xdz(1) ; v(1) = b1 * v0xdz(1)
1088	278	b1 = 1.0 / hlay(1)
1089	278	T(1) = b1 * T0xdz(1) ; Sal(1) = b1 * S0xdz(1)
1090		endif
1091
1092		! This uses half the harmonic mean of thicknesses to provide two estimates
1093		! of the boundary between cells, and the inverse of the harmonic mean to
1094		! weight the two estimates. The net effect is that interfaces around thin
1095		! layers have thin cells, and the total thickness adds up properly.
1096		! The top- and bottom- interfaces have zero thickness, consistent with
1097		! adding additional zero thickness layers.
1098	46250	h_Int(1) = 0.0 ; h_Int(2) = hlay(1)
1099	46250	dz_Int(1) = 0.0 ; dz_Int(2) = dz_lay(1)
1100	345671	do K=2,nzc-1
1101	299421	Norm = 1.0 / (hlay(k)(hlay(k-1)+hlay(k+1)) + 2.0hlay(k-1)*hlay(k+1))
1102	299421	wt_a = ((hlay(k)+hlay(k+1)) * hlay(k-1)) * Norm
1103	299421	wt_b = ((hlay(k-1)+hlay(k)) * hlay(k+1)) * Norm
1104	299421	h_Int(K) = h_Int(K) + hlay(k) * wt_a
1105	299421	h_Int(K+1) = hlay(k) * wt_b
1106	299421	dz_Int(K) = dz_Int(K) + dz_lay(k) * wt_a
1107	345671	dz_Int(K+1) = dz_lay(k) * wt_b
1108		enddo
1109	46250	h_Int(nzc) = h_Int(nzc) + hlay(nzc) ; h_Int(nzc+1) = 0.0
1110	46250	dz_Int(nzc) = dz_Int(nzc) + dz_lay(nzc) ; dz_Int(nzc+1) = 0.0
1111
1112	46250	if (GV%Boussinesq) then
1113	484143	do K=1,nzc+1 ; dz_h_Int(K) = GV%H_to_Z ; enddo
1114		else
1115		! Find an effective average specific volume around an interface.
1116	0	dz_h_Int(1:nzc+1) = 0.0
1117	0	if (hlay(1) > 0.0) dz_h_Int(1) = dz_lay(1) / hlay(1)
1118	0	do K=2,nzc+1
1119	0	if (h_Int(K) > 0.0) then
1120	0	dz_h_Int(K) = dz_Int(K) / h_Int(K)
1121		else
1122	0	dz_h_Int(K) = dz_h_Int(K-1)
1123		endif
1124		enddo
1125		endif
1126
1127		! Calculate thermodynamic coefficients and an initial estimate of N2.
1128	46250	if (use_temperature) then
1129	46250	pressure(1) = surface_pres
1130	391643	do K=2,nzc
1131	345393	pressure(K) = pressure(K-1) + gR0*hlay(k-1)
1132	345393	T_int(K) = 0.5*(T(k-1) + T(k))
1133	391643	Sal_int(K) = 0.5*(Sal(k-1) + Sal(k))
1134		enddo
1135	46250	if (GV%Boussinesq .or. GV%semi_Boussinesq) then
1136		call calculate_density_derivs(T_int, Sal_int, pressure, dbuoy_dT, dbuoy_dS, &
1137	138750	tv%eqn_of_state, (/2,nzc/), scale=-g_R0 )
1138		else
1139		! These should perhaps be combined into a single call to calculate the thermal expansion
1140		! and haline contraction coefficients?
1141		call calculate_specific_vol_derivs(T_int, Sal_int, pressure, dSpV_dT, dSpV_dS, &
1142	0	tv%eqn_of_state, (/2,nzc/) )
1143	0	call calculate_density(T_int, Sal_int, pressure, rho_int, tv%eqn_of_state, (/2,nzc/) )
1144	0	do K=2,nzc
1145	0	dbuoy_dT(K) = GV%g_Earth_Z_T2 * (rho_int(K) * dSpV_dT(K))
1146	0	dbuoy_dS(K) = GV%g_Earth_Z_T2 * (rho_int(K) * dSpV_dS(K))
1147		enddo
1148		endif
1149	0	elseif (GV%Boussinesq .or. GV%semi_Boussinesq) then
1150	0	do K=1,nzc+1 ; dbuoy_dT(K) = -g_R0 ; dbuoy_dS(K) = 0.0 ; enddo
1151		else
1152	0	do K=1,nzc+1 ; dbuoy_dS(K) = 0.0 ; enddo
1153	0	dbuoy_dT(1) = -GV%g_Earth_Z_T2 / GV%Rlay(1)
1154	0	do K=2,nzc
1155	0	dbuoy_dT(K) = -GV%g_Earth_Z_T2 / (0.5*(GV%Rlay(k-1) + GV%Rlay(k)))
1156		enddo
1157	0	dbuoy_dT(nzc+1) = -GV%g_Earth_Z_T2 / GV%Rlay(nzc)
1158		endif
1159
1160		! N2_debug(1) = 0.0 ; N2_debug(nzc+1) = 0.0
1161		! do K=2,nzc
1162		! N2_debug(K) = max((dbuoy_dT(K) * (T0xdz(k-1)Idz(k-1) - T0xdz(k)Idz(k)) + &
1163		! dbuoy_dS(K) * (S0xdz(k-1)Idz(k-1) - S0xdz(k)Idz(k))) * &
1164		! I_dz_int(K), 0.0)
1165		! enddo
1166
1167		! This call just calculates N2 and S2.
1168		call calculate_projected_state(kappa, u, v, T, Sal, 0.0, nzc, hlay, I_dz_int, dbuoy_dT, dbuoy_dS, &
1169	46250	CS%vel_underflow, u, v, T, Sal, N2, S2, GV, US)
1170	484143	do K=1,nzc+1
1171	437893	N2_init(K) = N2(K)
1172	484143	S2_init(K) = S2(K)
1173		enddo
1174
1175
1176		! ----------------------------------------------------
1177		! Iterate
1178		! ----------------------------------------------------
1179	46250	dt_rem = dt
1180	484143	do K=1,nzc+1
1181	437893	K_Q(K) = 0.0
1182	437893	kappa_avg(K) = 0.0 ; tke_avg(K) = 0.0 ; N2_mean(K) = 0.0 ; S2_mean(K) = 0.0
1183	437893	local_src_avg(K) = 0.0
1184		! Use the grid spacings to scale errors in the source.
1185	437893	if ( h_Int(K) > 0.0 ) &
1186	391921	local_src_avg(K) = 0.1 * k0dt * I_dz_int(K) / h_Int(K)
1187		enddo
1188
1189		! call cpu_clock_end(id_clock_setup)
1190
1191		! do itt=1,CS%max_RiNo_it
1192	46250	do itt=1,CS%max_KS_it
1193
1194		! ----------------------------------------------------
1195		! Calculate new values of u, v, rho, N^2 and S.
1196		! ----------------------------------------------------
1197
1198		! call cpu_clock_begin(id_clock_KQ)
1199		call find_kappa_tke(N2, S2, kappa, Idz, h_Int, dz_Int, dz_h_Int, I_L2_bdry, f2, &
1200	46250	nzc, CS, GV, US, K_Q, tke, kappa_out, kappa_src, local_src)
1201		! call cpu_clock_end(id_clock_KQ)
1202
1203		! call cpu_clock_begin(id_clock_avg)
1204		! Determine the range of non-zero values of kappa_out.
1205	46250	ks_kappa = GV%ke+1 ; ke_kappa = 0
1206	373558	do K=2,nzc ; if (kappa_out(K) > 0.0) then
1207	3584	ks_kappa = K ; exit
1208		endif ; enddo
1209	60409	do k=nzc,ks_kappa,-1 ; if (kappa_out(K) > 0.0) then
1210	3584	ke_kappa = K ; exit
1211		endif ; enddo
1212	46250	if (ke_kappa == nzc) kappa_out(nzc+1) = 0.0
1213		! call cpu_clock_end(id_clock_avg)
1214
1215		! Determine how long to use this value of kappa (dt_now).
1216
1217		! call cpu_clock_begin(id_clock_project)
1218	46250	if ((ke_kappa < ks_kappa) .or. (itt==CS%max_KS_it)) then
1219	42666	dt_now = dt_rem
1220		else
1221		! Limit dt_now so that \|k_src(k)-kappa_src(k)\| < tol * local_src(k)
1222	3584	dt_test = dt_rem
1223	21669	do K=2,nzc
1224	18085	tol_max(K) = kappa_src(K) + tol_dksrc * local_src(K)
1225	18085	tol_min(K) = kappa_src(K) - tol_dksrc_low * local_src(K)
1226	21669	tol_chg(K) = tol2 * local_src_avg(K)
1227		enddo
1228
1229	3584	do itt_dt=1,(CS%max_KS_it+1-itt)/2
1230		! The maximum number of times that the time-step is halved in
1231		! seeking an acceptable timestep is reduced with each iteration,
1232		! so that as the maximum number of iterations is approached, the
1233		! whole remaining timestep is used. Typically, an acceptable
1234		! timestep is found long before the minimum is reached, so the
1235		! value of max_KS_it may be unimportant, especially if it is large
1236		! enough.
1237		call calculate_projected_state(kappa_out, u, v, T, Sal, 0.5*dt_test, nzc, hlay, I_dz_int, &
1238		dbuoy_dT, dbuoy_dS, CS%vel_underflow, u_test, v_test, &
1239	3584	T_test, S_test, N2, S2, GV, US, ks_int=ks_kappa, ke_int=ke_kappa)
1240	3584	valid_dt = .true.
1241	3584	Idtt = 1.0 / dt_test
1242	11094	do K=max(ks_kappa-1,2),min(ke_kappa+1,nzc)
1243	11094	if (N2(K) < Ri_crit * S2(K)) then ! Equivalent to Ri < Ri_crit.
1244		K_src(K) = (2.0 * CS%Shearmix_rate * sqrt(S2(K))) * &
1245	3584	((Ri_critS2(K) - N2(K)) / (Ri_critS2(K) + CS%FRi_curvature*N2(K)))
1246	3584	if (CS%restrictive_tolerance_check) then
1247	0	if ((K_src(K) > min(tol_max(K), kappa_src(K) + Idtt*tol_chg(K))) .or. &
1248		(K_src(K) < max(tol_min(K), kappa_src(K) - Idtt*tol_chg(K)))) then
1249	0	valid_dt = .false. ; exit
1250		endif
1251		else
1252	3584	if ((K_src(K) > max(tol_max(K), kappa_src(K) + Idtt*tol_chg(K))) .or. &
1253		(K_src(K) < min(tol_min(K), kappa_src(K) - Idtt*tol_chg(K)))) then
1254	0	valid_dt = .false. ; exit
1255		endif
1256		endif
1257		else
1258	3926	if (0.0 < min(tol_min(K), kappa_src(K) - Idtt*tol_chg(K))) then
1259	0	valid_dt = .false. ; k_src(K) = 0.0 ; exit
1260		endif
1261		endif
1262		enddo
1263
1264	3584	if (valid_dt) exit
1265	0	dt_test = 0.5*dt_test
1266		enddo
1267	3584	if ((dt_test < dt_rem) .and. valid_dt) then
1268	0	dt_inc = 0.5*dt_test
1269	0	do itt_dt=1,dt_refinements
1270		call calculate_projected_state(kappa_out, u, v, T, Sal, 0.5*(dt_test+dt_inc), nzc, hlay, &
1271		I_dz_int, dbuoy_dT, dbuoy_dS, CS%vel_underflow, u_test, v_test, T_test, S_test, &
1272	0	N2, S2, GV, US, ks_int=ks_kappa, ke_int=ke_kappa)
1273	0	valid_dt = .true.
1274	0	Idtt = 1.0 / (dt_test+dt_inc)
1275	0	do K=max(ks_kappa-1,2),min(ke_kappa+1,nzc)
1276	0	if (N2(K) < Ri_crit * S2(K)) then ! Equivalent to Ri < Ri_crit.
1277		K_src(K) = (2.0 * CS%Shearmix_rate * sqrt(S2(K))) * &
1278	0	((Ri_critS2(K) - N2(K)) / (Ri_critS2(K) + CS%FRi_curvature*N2(K)))
1279	0	if ((K_src(K) > max(tol_max(K), kappa_src(K) + Idtt*tol_chg(K))) .or. &
1280		(K_src(K) < min(tol_min(K), kappa_src(K) - Idtt*tol_chg(K)))) then
1281	0	valid_dt = .false. ; exit
1282		endif
1283		else
1284	0	if (0.0 < min(tol_min(K), kappa_src(K) - Idtt*tol_chg(K))) then
1285	0	valid_dt = .false. ; k_src(K) = 0.0 ; exit
1286		endif
1287		endif
1288		enddo
1289
1290	0	if (valid_dt) dt_test = dt_test + dt_inc
1291	0	dt_inc = 0.5*dt_inc
1292		enddo
1293		else
1294	3584	dt_inc = 0.0
1295		endif
1296
1297	3584	dt_now = min(dt_test*(1.0+CS%kappa_tol_err)+dt_inc, dt_rem)
1298	21669	do K=2,nzc
1299	21669	local_src_avg(K) = local_src_avg(K) + dt_now * local_src(K)
1300		enddo
1301		endif ! Are all the values of kappa_out 0?
1302		! call cpu_clock_end(id_clock_project)
1303
1304		! The state has already been projected forward. Now find new values of kappa.
1305
1306	46250	if (ke_kappa < ks_kappa) then
1307		! There is no mixing now, and will not be again.
1308		! call cpu_clock_begin(id_clock_avg)
1309	42666	dt_wt = dt_rem / dt ; dt_rem = 0.0
1310	455306	do K=1,nzc+1
1311	412640	kappa_mid(K) = 0.0
1312		! This would be here but does nothing.
1313		! kappa_avg(K) = kappa_avg(K) + kappa_mid(K)*dt_wt
1314	455306	tke_avg(K) = tke_avg(K) + dt_wt*tke(K)
1315		enddo
1316		! call cpu_clock_end(id_clock_avg)
1317		else
1318		! call cpu_clock_begin(id_clock_project)
1319		call calculate_projected_state(kappa_out, u, v, T, Sal, dt_now, nzc, hlay, I_dz_int, &
1320		dbuoy_dT, dbuoy_dS, CS%vel_underflow, u_test, v_test, &
1321	3584	T_test, S_test, N2, S2, GV, US, ks_int=ks_kappa, ke_int=ke_kappa)
1322		! call cpu_clock_end(id_clock_project)
1323
1324		! call cpu_clock_begin(id_clock_KQ)
1325	28837	do K=1,nzc+1 ; K_Q_tmp(K) = K_Q(K) ; enddo
1326		call find_kappa_tke(N2, S2, kappa_out, Idz, h_Int, dz_Int, dz_h_Int, I_L2_bdry, f2, &
1327	3584	nzc, CS, GV, US, K_Q_tmp, tke_pred, kappa_pred)
1328		! call cpu_clock_end(id_clock_KQ)
1329
1330	3584	ks_kappa = GV%ke+1 ; ke_kappa = 0
1331	28837	do K=1,nzc+1
1332	25253	kappa_mid(K) = 0.5*(kappa_out(K) + kappa_pred(K))
1333	25253	if ((kappa_mid(K) > 0.0) .and. (K<ks_kappa)) ks_kappa = K
1334	28837	if (kappa_mid(K) > 0.0) ke_kappa = K
1335		enddo
1336
1337		! call cpu_clock_begin(id_clock_project)
1338		call calculate_projected_state(kappa_mid, u, v, T, Sal, dt_now, nzc, hlay, I_dz_int, &
1339		dbuoy_dT, dbuoy_dS, CS%vel_underflow, u_test, v_test, &
1340	3584	T_test, S_test, N2, S2, GV, US, ks_int=ks_kappa, ke_int=ke_kappa)
1341		! call cpu_clock_end(id_clock_project)
1342
1343		! call cpu_clock_begin(id_clock_KQ)
1344		call find_kappa_tke(N2, S2, kappa_out, Idz, h_Int, dz_Int, dz_h_Int, I_L2_bdry, f2, &
1345	3584	nzc, CS, GV, US, K_Q, tke_pred, kappa_pred)
1346		! call cpu_clock_end(id_clock_KQ)
1347
1348		! call cpu_clock_begin(id_clock_avg)
1349	3584	dt_wt = dt_now / dt ; dt_rem = dt_rem - dt_now
1350	28837	do K=1,nzc+1
1351	25253	kappa_mid(K) = 0.5*(kappa_out(K) + kappa_pred(K))
1352	25253	kappa_avg(K) = kappa_avg(K) + kappa_mid(K)*dt_wt
1353	25253	tke_avg(K) = tke_avg(K) + dt_wt0.5(tke_pred(K) + tke(K))
1354	25253	N2_mean(K) = N2_mean(K) + dt_wt*N2(K)
1355	25253	S2_mean(K) = S2_mean(K) + dt_wt*S2(K)
1356	28837	kappa(K) = kappa_pred(K) ! First guess for the next iteration.
1357		enddo
1358
1359		! call cpu_clock_end(id_clock_avg)
1360		endif
1361
1362	46250	if (dt_rem > 0.0) then
1363		! Update the values of u, v, T, Sal, N2, and S2 for the next iteration.
1364		! call cpu_clock_begin(id_clock_project)
1365		call calculate_projected_state(kappa_mid, u, v, T, Sal, dt_now, nzc, hlay, I_dz_int, &
1366		dbuoy_dT, dbuoy_dS, CS%vel_underflow, u, v, T, Sal, N2, S2, &
1367	0	GV, US)
1368		! call cpu_clock_end(id_clock_project)
1369		endif
1370
1371	46250	if (dt_rem <= 0.0) exit
1372
1373		enddo ! end itt loop
1374
1375	46250	end subroutine kappa_shear_column
1376
1377		!> This subroutine calculates the velocities, temperature and salinity that
1378		!! the water column will have after mixing for dt with diffusivities kappa. It
1379		!! may also calculate the projected buoyancy frequency and shear.
1380	57002	subroutine calculate_projected_state(kappa, u0, v0, T0, S0, dt, nz, dz, I_dz_int, dbuoy_dT, dbuoy_dS, &
1381	57002	vel_under, u, v, T, Sal, N2, S2, GV, US, ks_int, ke_int)
1382		integer, intent(in) :: nz !< The number of layers (after eliminating massless
1383		!! layers?).
1384		real, dimension(nz+1), intent(in) :: kappa !< The diapycnal diffusivity at interfaces,
1385		!! [H Z T-1 ~> m2 s-1 or Pa s].
1386		real, dimension(nz), intent(in) :: u0 !< The initial zonal velocity [L T-1 ~> m s-1].
1387		real, dimension(nz), intent(in) :: v0 !< The initial meridional velocity [L T-1 ~> m s-1].
1388		real, dimension(nz), intent(in) :: T0 !< The initial temperature [C ~> degC].
1389		real, dimension(nz), intent(in) :: S0 !< The initial salinity [S ~> ppt].
1390		real, intent(in) :: dt !< The time step [T ~> s].
1391		real, dimension(nz), intent(in) :: dz !< The layer thicknesses [H ~> m or kg m-2]
1392		real, dimension(nz+1), intent(in) :: I_dz_int !< The inverse of the distance between successive
1393		!! layer centers [Z-1 ~> m-1].
1394		real, dimension(nz+1), intent(in) :: dbuoy_dT !< The partial derivative of buoyancy with
1395		!! temperature [Z T-2 C-1 ~> m s-2 degC-1].
1396		real, dimension(nz+1), intent(in) :: dbuoy_dS !< The partial derivative of buoyancy with
1397		!! salinity [Z T-2 S-1 ~> m s-2 ppt-1].
1398		real, intent(in) :: vel_under !< Any velocities that are smaller in magnitude
1399		!! than this value are set to 0 [L T-1 ~> m s-1].
1400		real, dimension(nz), intent(inout) :: u !< The zonal velocity after dt [L T-1 ~> m s-1].
1401		real, dimension(nz), intent(inout) :: v !< The meridional velocity after dt [L T-1 ~> m s-1].
1402		real, dimension(nz), intent(inout) :: T !< The temperature after dt [C ~> degC].
1403		real, dimension(nz), intent(inout) :: Sal !< The salinity after dt [S ~> ppt].
1404		real, dimension(nz+1), intent(inout) :: N2 !< The buoyancy frequency squared at interfaces [T-2 ~> s-2].
1405		real, dimension(nz+1), intent(inout) :: S2 !< The squared shear at interfaces [T-2 ~> s-2].
1406		type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure.
1407		type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
1408		integer, optional, intent(in) :: ks_int !< The topmost k-index with a non-zero diffusivity.
1409		integer, optional, intent(in) :: ke_int !< The bottommost k-index with a non-zero
1410		!! diffusivity.
1411
1412		! Local variables
1413	57002	real, dimension(nz+1) :: c1 ! A tridiagonal variable [nondim]
1414		real :: a_a, a_b ! Tridiagonal coupling coefficients [H ~> m or kg m-2]
1415		real :: b1, b1nz_0 ! Tridiagonal variables [H-1 ~> m-1 or m2 kg-1]
1416		real :: bd1 ! A term in the denominator of b1 [H ~> m or kg m-2]
1417		real :: d1 ! A tridiagonal variable [nondim]
1418		integer :: k, ks, ke
1419
1420	57002	ks = 1 ; ke = nz
1421	57002	if (present(ks_int)) ks = max(ks_int-1,1)
1422	57002	if (present(ke_int)) ke = min(ke_int,nz)
1423
1424	57002	if (ks > ke) return
1425
1426	57002	if (dt > 0.0) then
1427	10752	a_b = dt(kappa(ks+1)I_dz_int(ks+1))
1428	10752	b1 = 1.0 / (dz(ks) + a_b)
1429	10752	c1(ks+1) = a_b * b1 ; d1 = dz(ks) * b1 ! = 1 - c1
1430
1431	10752	u(ks) = (b1 * dz(ks))u0(ks) ; v(ks) = (b1 dz(ks))*v0(ks)
1432	10752	T(ks) = (b1 * dz(ks))T0(ks) ; Sal(ks) = (b1 dz(ks))*S0(ks)
1433	11778	do K=ks+1,ke-1
1434	1026	a_a = a_b
1435	1026	a_b = dt(kappa(K+1)I_dz_int(K+1))
1436	1026	bd1 = dz(k) + d1*a_a
1437	1026	b1 = 1.0 / (bd1 + a_b)
1438	1026	c1(K+1) = a_b * b1 ; d1 = bd1 * b1 ! d1 = 1 - c1
1439
1440	1026	u(k) = b1 * (dz(k)u0(k) + a_au(k-1))
1441	1026	v(k) = b1 * (dz(k)v0(k) + a_av(k-1))
1442	1026	T(k) = b1 * (dz(k)T0(k) + a_aT(k-1))
1443	11778	Sal(k) = b1 * (dz(k)S0(k) + a_aSal(k-1))
1444		enddo
1445		! T and S have insulating boundary conditions, u & v use no-slip
1446		! bottom boundary conditions at the solid bottom.
1447
1448		! For insulating boundary conditions or mixing simply stopping, use...
1449	10752	a_a = a_b
1450	10752	b1 = 1.0 / (dz(ke) + d1*a_a)
1451	10752	T(ke) = b1 * (dz(ke)T0(ke) + a_aT(ke-1))
1452	10752	Sal(ke) = b1 * (dz(ke)S0(ke) + a_aSal(ke-1))
1453
1454		! There is no distinction between the effective boundary conditions for
1455		! tracers and velocities if the mixing is separated from the bottom, but if
1456		! the mixing goes all the way to the bottom, use no-slip BCs for velocities.
1457	10752	if (ke == nz) then
1458	0	a_b = dt(kappa(nz+1)I_dz_int(nz+1))
1459	0	b1nz_0 = 1.0 / ((dz(nz) + d1*a_a) + a_b)
1460		else
1461	10752	b1nz_0 = b1
1462		endif
1463	10752	u(ke) = b1nz_0 * (dz(ke)u0(ke) + a_au(ke-1))
1464	10752	v(ke) = b1nz_0 * (dz(ke)v0(ke) + a_av(ke-1))
1465	10752	if (abs(u(ke)) < vel_under) u(ke) = 0.0
1466	10752	if (abs(v(ke)) < vel_under) v(ke) = 0.0
1467
1468	22530	do k=ke-1,ks,-1
1469	11778	u(k) = u(k) + c1(k+1)*u(k+1)
1470	11778	v(k) = v(k) + c1(k+1)*v(k+1)
1471	11778	if (abs(u(k)) < vel_under) u(k) = 0.0
1472	11778	if (abs(v(k)) < vel_under) v(k) = 0.0
1473	11778	T(k) = T(k) + c1(k+1)*T(k+1)
1474	22530	Sal(k) = Sal(k) + c1(k+1)*Sal(k+1)
1475		enddo
1476		else ! dt <= 0.0
1477	437893	do k=1,nz
1478	391643	u(k) = u0(k) ; v(k) = v0(k) ; T(k) = T0(k) ; Sal(k) = S0(k)
1479	391643	if (abs(u(k)) < vel_under) u(k) = 0.0
1480	437893	if (abs(v(k)) < vel_under) v(k) = 0.0
1481		enddo
1482		endif
1483
1484		! Store the squared shear at interfaces
1485	57002	S2(1) = 0.0 ; S2(nz+1) = 0.0
1486	57002	if (ks > 1) &
1487	0	S2(ks) = (((u(ks)-u0(ks-1))2) + ((v(ks)-v0(ks-1))2)) * (US%L_to_ZI_dz_int(ks))*2
1488	414173	do K=ks+1,ke
1489	414173	S2(K) = (((u(k)-u(k-1))2) + ((v(k)-v(k-1))2)) * (US%L_to_ZI_dz_int(K))*2
1490		enddo
1491	57002	if (ke<nz) &
1492	10752	S2(ke+1) = (((u0(ke+1)-u(ke))2) + ((v0(ke+1)-v(ke))2)) * (US%L_to_ZI_dz_int(ke+1))*2
1493
1494		! Store the buoyancy frequency at interfaces
1495	57002	N2(1) = 0.0 ; N2(nz+1) = 0.0
1496	57002	if (ks > 1) &
1497		N2(ks) = max(0.0, I_dz_int(ks) * &
1498	0	(dbuoy_dT(ks) * (T0(ks-1)-T(ks)) + dbuoy_dS(ks) * (S0(ks-1)-Sal(ks))))
1499	414173	do K=ks+1,ke
1500		N2(K) = max(0.0, I_dz_int(K) * &
1501	414173	(dbuoy_dT(K) * (T(k-1)-T(k)) + dbuoy_dS(K) * (Sal(k-1)-Sal(k))))
1502		enddo
1503	57002	if (ke<nz) &
1504		N2(ke+1) = max(0.0, I_dz_int(ke+1) * &
1505	10752	(dbuoy_dT(ke+1) * (T(ke)-T0(ke+1)) + dbuoy_dS(ke+1) * (Sal(ke)-S0(ke+1))))
1506
1507		end subroutine calculate_projected_state
1508
1509		!> This subroutine calculates new, consistent estimates of TKE and kappa.
1510	106836	subroutine find_kappa_tke(N2, S2, kappa_in, Idz, h_Int, dz_Int, dz_h_Int, I_L2_bdry, f2, &
1511	92500	nz, CS, GV, US, K_Q, tke, kappa, kappa_src, local_src)
1512		integer, intent(in) :: nz !< The number of layers to work on.
1513		real, dimension(nz+1), intent(in) :: N2 !< The buoyancy frequency squared at interfaces [T-2 ~> s-2].
1514		real, dimension(nz+1), intent(in) :: S2 !< The squared shear at interfaces [T-2 ~> s-2].
1515		real, dimension(nz+1), intent(in) :: kappa_in !< The initial guess at the diffusivity
1516		!! [H Z T-1 ~> m2 s-1 or Pa s]
1517		real, dimension(nz+1), intent(in) :: h_Int !< The thicknesses associated with interfaces
1518		!! [H ~> m or kg m-2]
1519		real, dimension(nz+1), intent(in) :: dz_Int !< The vertical distances around interfaces [Z ~> m]
1520		real, dimension(nz+1), intent(in) :: dz_h_Int !< The ratio of the vertical distances to the
1521		!! thickness around an interface [Z H-1 ~> nondim or m3 kg-1].
1522		!! In non-Boussinesq mode this is the specific volume.
1523		real, dimension(nz+1), intent(in) :: I_L2_bdry !< The inverse of the squared distance to
1524		!! boundaries [H-1 Z-1 ~> m-2 or m kg-1].
1525		real, dimension(nz), intent(in) :: Idz !< The inverse grid spacing of layers [Z-1 ~> m-1].
1526		real, intent(in) :: f2 !< The squared Coriolis parameter [T-2 ~> s-2].
1527		type(Kappa_shear_CS), pointer :: CS !< A pointer to this module's control structure.
1528		type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure.
1529		type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
1530		real, dimension(nz+1), intent(inout) :: K_Q !< The shear-driven diapycnal diffusivity divided by
1531		!! the turbulent kinetic energy per unit mass at
1532		!! interfaces [H T Z-1 ~> s or kg s m-3].
1533		real, dimension(nz+1), intent(out) :: tke !< The turbulent kinetic energy per unit mass at
1534		!! interfaces [Z2 T-2 ~> m2 s-2].
1535		real, dimension(nz+1), intent(out) :: kappa !< The diapycnal diffusivity at interfaces
1536		!! [H Z T-1 ~> m2 s-1 or Pa s]
1537		real, dimension(nz+1), optional, &
1538		intent(out) :: kappa_src !< The source term for kappa [T-1 ~> s-1]
1539		real, dimension(nz+1), optional, &
1540		intent(out) :: local_src !< The sum of all local sources for kappa
1541		!! [T-1 ~> s-1]
1542		! This subroutine calculates new, consistent estimates of TKE and kappa.
1543
1544		! Local variables
1545		real, dimension(nz) :: &
1546	10752	aQ, & ! aQ is the coupling between adjacent interfaces in the TKE equations [H T-1 ~> m s-1 or kg m-2 s-1]
1547	106836	dQdz ! Half the partial derivative of TKE with depth [Z T-2 ~> m s-2].
1548		real, dimension(nz+1) :: &
1549	106836	dK, & ! The change in kappa [H Z T-1 ~> m2 s-1 or Pa s].
1550	106836	dQ, & ! The change in TKE [Z2 T-2 ~> m2 s-2].
1551	106836	cQ, cK, & ! cQ and cK are the upward influences in the tridiagonal and
1552		! hexadiagonal solvers for the TKE and kappa equations [nondim].
1553	106836	I_Ld2, & ! 1/Ld^2, where Ld is the effective decay length scale for kappa [H-1 Z-1 ~> m-2 or m kg-1]
1554	106836	TKE_decay, & ! The local TKE decay rate [T-1 ~> s-1].
1555	106836	k_src, & ! The source term in the kappa equation [T-1 ~> s-1].
1556	106836	dQmdK, & ! With Newton's method the change in dQ(k-1) due to dK(k) [Z T H-1 ~> s or m3 s kg-1]
1557	106836	dKdQ, & ! With Newton's method the change in dK(k) due to dQ(k) [H Z-1 T-1 ~> s-1 or kg m-3 s-1]
1558	106836	e1 ! The fractional change in a layer TKE due to a change in the
1559		! TKE of the layer above when all the kappas below are 0 [nondim].
1560		! e1 is nondimensional, and 0 < e1 < 1.
1561		real :: tke_src ! The net source of TKE due to mixing against the shear and stratification
1562		! [Z2 T-3 ~> m2 s-3] or [H Z T-3 ~> m2 s-3 or kg m-1 s-3].
1563		! (For convenience, a term involving the non-dissipation of q0 is also included here.)
1564		real :: bQ ! The inverse of the pivot in the tridiagonal equations [T H-1 ~> s m-1 or m2 s kg-1]
1565		real :: bK ! The inverse of the pivot in the tridiagonal equations [Z-1 ~> m-1].
1566		real :: bQd1 ! A term in the denominator of bQ [H T-1 ~> m s-1 or kg m-2 s-1]
1567		real :: bKd1 ! A term in the denominator of bK [Z ~> m].
1568		real :: cQcomp, cKcomp ! 1 - cQ or 1 - cK in the tridiagonal equations [nondim].
1569		real :: c_s2 ! The coefficient for the decay of TKE due to
1570		! shear (i.e. proportional to \|S\|*tke) [nondim].
1571		real :: c_n2 ! The coefficient for the decay of TKE due to
1572		! stratification (i.e. proportional to N*tke) [nondim].
1573		real :: Ri_crit ! The critical shear Richardson number for shear-
1574		! driven mixing [nondim]. The theoretical value is 0.25.
1575		real :: q0 ! The background level of TKE [Z2 T-2 ~> m2 s-2].
1576		real :: Ilambda2 ! 1.0 / CS%lambda**2 [nondim]
1577		real :: TKE_min ! The minimum value of shear-driven TKE that can be
1578		! solved for [Z2 T-2 ~> m2 s-2].
1579		real :: kappa0 ! The background diapycnal diffusivity [H Z T-1 ~> m2 s-1 or Pa s]
1580		real :: kappa_trunc ! Diffusivities smaller than this are rounded to 0 [H Z T-1 ~> m2 s-1 or Pa s]
1581
1582		real :: eden1, eden2 ! Variables used in calculating e1 [H Z-2 ~> m-1 or kg m-4]
1583		real :: I_eden ! The inverse of the denominator in e1 [Z2 H-1 ~> m or m4 kg-1]
1584		real :: ome ! Variables used in calculating e1 [nondim]
1585		real :: diffusive_src ! The diffusive source in the kappa equation [H T-1 ~> m s-1 or kg m-2 s-1]
1586		real :: chg_by_k0 ! The value of k_src that leads to an increase of
1587		! kappa_0 if only the diffusive term is a sink [T-1 ~> s-1]
1588		real :: h_dz_here ! The ratio of the thicknesses to the vertical distances around an interface
1589		! [H Z-1 ~> nondim or kg m-3]. In non-Boussinesq mode this is the density.
1590
1591		real :: kappa_mean ! A mean value of kappa [H Z T-1 ~> m2 s-1 or Pa s]
1592		real :: Newton_test ! The value of relative error that will cause the next
1593		! iteration to use Newton's method [nondim].
1594		! Temporary variables used in the Newton's method iterations.
1595		real :: decay_term_k ! The decay term in the diffusivity equation [Z-1 ~> m-1]
1596		real :: decay_term_Q ! The decay term in the TKE equation - proportional to [H Z-1 T-1 ~> s-1 or kg m-3 s-1]
1597		real :: I_Q ! The inverse of TKE [T2 Z-2 ~> s2 m-2]
1598		real :: kap_src ! A source term in the kappa equation [H T-1 ~> m s-1 or kg m-2 s-1]
1599		real :: v1 ! A temporary variable proportional to [H Z-1 T-1 ~> s-1 or kg m-3 s-1]
1600		real :: v2 ! A temporary variable in [Z T-2 ~> m s-2]
1601		real :: tol_err ! The tolerance for max_err that determines when to
1602		! stop iterating [nondim].
1603		real :: Newton_err ! The tolerance for max_err that determines when to
1604		! start using Newton's method [nondim]. Empirically, an initial
1605		! value of about 0.2 seems to be most efficient.
1606		real, parameter :: roundoff = 1.0e-16 ! A negligible fractional change in TKE [nondim].
1607		! This could be larger but performance gains are small.
1608
1609		logical :: tke_noflux_bottom_BC = .false. ! Specify the boundary conditions
1610		logical :: tke_noflux_top_BC = .false. ! that are applied to the TKE equations.
1611		logical :: do_Newton ! If .true., use Newton's method for the next iteration.
1612		logical :: abort_Newton ! If .true., an Newton's method has encountered a 0
1613		! pivot, and should not have been used.
1614		logical :: was_Newton ! The value of do_Newton before checking convergence.
1615		logical :: within_tolerance ! If .true., all points are within tolerance to
1616		! enable this subroutine to return.
1617		integer :: ks_src, ke_src ! The range indices that have nonzero k_src.
1618		integer :: ks_kappa, ke_kappa, ke_tke ! The ranges of k-indices that are or
1619		integer :: ks_kappa_prev, ke_kappa_prev ! were being worked on.
1620		integer :: itt, k, k2
1621
1622		! These variables are used only for debugging.
1623		logical, parameter :: debug_soln = .false.
1624		real :: K_err_lin ! The imbalance in the K equation [H T-1 ~> m s-1 or kg m-2 s-1]
1625		real :: Q_err_lin ! The imbalance in the Q equation [H Z T-3 ~> m2 s-3 or kg m-1 s-3]
1626		real, dimension(nz+1) :: &
1627	106836	I_Ld2_debug, & ! A separate version of I_Ld2 for debugging [H-1 Z-1 ~> m-2 or m kg-1].
1628	106836	kappa_prev, & ! The value of kappa at the start of the current iteration [H Z T-1 ~> m2 s-1 or Pa s]
1629	53418	TKE_prev ! The value of TKE at the start of the current iteration [Z2 T-2 ~> m2 s-2].
1630
1631	53418	c_N2 = CS%C_N2 ; c_S2 = CS%C_S2
1632	53418	q0 = CS%TKE_bg ; kappa0 = CS%kappa_0
1633	53418	TKE_min = max(CS%TKE_bg, 1.0E-20US%m_to_Z2US%T_to_s**2)
1634	53418	Ri_crit = CS%Rino_crit
1635	53418	Ilambda2 = 1.0 / CS%lambda**2
1636	53418	kappa_trunc = CS%kappa_trunc
1637	53418	do_Newton = .false. ; abort_Newton = .false.
1638	53418	tol_err = CS%kappa_tol_err
1639	53418	Newton_err = 0.2 ! This initial value may be automatically reduced later.
1640
1641	53418	ks_kappa = 2 ; ke_kappa = nz ; ks_kappa_prev = 2 ; ke_kappa_prev = nz
1642
1643	53418	ke_src = 0 ; ks_src = nz+1
1644	434981	do K=2,nz
1645	434981	if (N2(K) < Ri_crit * S2(K)) then ! Equivalent to Ri < Ri_crit.
1646		! Ri = N2(K) / S2(K)
1647		! k_src(K) = (2.0 * CS%Shearmix_rate * sqrt(S2(K))) * &
1648		! ((Ri_crit - Ri) / (Ri_crit + CS%FRi_curvature*Ri))
1649		K_src(K) = (2.0 * CS%Shearmix_rate * sqrt(S2(K))) * &
1650	10752	((Ri_critS2(K) - N2(K)) / (Ri_critS2(K) + CS%FRi_curvature*N2(K)))
1651	10752	ke_src = K
1652	10752	if (ks_src > k) ks_src = K
1653		else
1654	370811	k_src(K) = 0.0
1655		endif
1656		enddo
1657
1658		! If there is no source anywhere, return kappa(K) = 0.
1659	53418	if (ks_src > ke_src) then
1660	455306	do K=1,nz+1
1661	455306	kappa(K) = 0.0 ; K_Q(K) = 0.0 ; tke(K) = TKE_min
1662		enddo
1663	455306	if (present(kappa_src)) then ; do K=1,nz+1 ; kappa_src(K) = 0.0 ; enddo ; endif
1664	455306	if (present(local_src)) then ; do K=1,nz+1 ; local_src(K) = 0.0 ; enddo ; endif
1665	42666	return
1666		endif
1667
1668	86511	do K=1,nz+1
1669	75759	kappa(K) = kappa_in(K)
1670		! TKE_decay(K) = c_nsqrt(N2(K)) + c_ssqrt(S2(K)) ! The expression in JHL.
1671	75759	TKE_decay(K) = sqrt(c_n2N2(K) + c_s2S2(K))
1672	86511	if ((kappa(K) > 0.0) .and. (K_Q(K) > 0.0)) then
1673	7852	TKE(K) = kappa(K) / K_Q(K) ! Perhaps take the max with TKE_min
1674		else
1675	67907	TKE(K) = TKE_min
1676		endif
1677		enddo
1678		! Apply boundary conditions to kappa.
1679	10752	kappa(1) = 0.0 ; kappa(nz+1) = 0.0
1680
1681		! Calculate the term (e1) that allows changes in TKE to be calculated quickly
1682		! below the deepest nonzero value of kappa. If kappa = 0, below interface
1683		! k-1, the final changes in TKE are related by dQ(K+1) = e1(K+1)*dQ(K).
1684	10752	eden2 = kappa0 * Idz(nz)
1685	10752	if (tke_noflux_bottom_BC) then
1686	0	eden1 = h_Int(nz+1)*TKE_decay(nz+1)
1687	0	I_eden = 1.0 / (eden2 + eden1)
1688	0	e1(nz+1) = eden2 * I_eden ; ome = eden1 * I_eden
1689		else
1690	10752	e1(nz+1) = 0.0 ; ome = 1.0
1691		endif
1692	65007	do k=nz,2,-1
1693	54255	eden1 = h_Int(K)TKE_decay(K) + ome eden2
1694	54255	eden2 = kappa0 * Idz(k-1)
1695	54255	I_eden = 1.0 / (eden2 + eden1)
1696	65007	e1(K) = eden2 * I_eden ; ome = eden1 * I_eden ! = 1-e1
1697		enddo
1698	10752	e1(1) = 0.0
1699
1700
1701		! Iterate here to convergence to within some tolerance of order tol_err.
1702	35562	do itt=1,CS%max_RiNo_it
1703
1704		! ----------------------------------------------------
1705		! Calculate TKE
1706		! ----------------------------------------------------
1707
1708		if (debug_soln) then ; do K=1,nz+1 ; kappa_prev(K) = kappa(K) ; TKE_prev(K) = TKE(K) ; enddo ; endif
1709
1710	35562	if (.not.do_Newton) then
1711		! Use separate steps of the TKE and kappa equations, that are
1712		! explicit in the nonlinear source terms, implicit in a linearized
1713		! version of the nonlinear sink terms, and implicit in the linear
1714		! terms.
1715
1716	35562	ke_tke = max(ke_kappa,ke_kappa_prev)+1
1717		! aQ is the coupling between adjacent interfaces [Z T-1 ~> m s-1].
1718	199129	do k=1,min(ke_tke,nz)
1719	199129	aQ(k) = (0.5(kappa(K)+kappa(K+1)) + kappa0) Idz(k)
1720		enddo
1721	35562	dQ(1) = -TKE(1)
1722	35562	if (tke_noflux_top_BC) then
1723	0	tke_src = dz_h_Int(1)kappa0S2(1) + q0 * TKE_decay(1) ! Uses that kappa(1) = 0
1724	0	bQd1 = h_Int(1) * TKE_decay(1)
1725	0	bQ = 1.0 / (bQd1 + aQ(1))
1726	0	tke(1) = bQ * (h_Int(1)*tke_src)
1727	0	cQ(2) = aQ(1) * bQ ; cQcomp = bQd1 * bQ ! = 1 - cQ
1728		else
1729	35562	tke(1) = q0 ; cQ(2) = 0.0 ; cQcomp = 1.0
1730		endif
1731	144193	do K=2,ke_tke-1
1732	108631	dQ(K) = -TKE(K)
1733	108631	tke_src = dz_h_Int(K)(kappa(K) + kappa0)S2(K) + q0*TKE_decay(K)
1734	108631	bQd1 = h_Int(K)(TKE_decay(K) + dz_h_Int(K)N2(K)K_Q(K)) + cQcompaQ(k-1)
1735	108631	bQ = 1.0 / (bQd1 + aQ(k))
1736	108631	tke(K) = bQ * (h_Int(K)tke_src + aQ(k-1)tke(K-1))
1737	144193	cQ(K+1) = aQ(k) * bQ ; cQcomp = bQd1 * bQ ! = 1 - cQ
1738		enddo
1739	35562	if ((ke_tke == nz+1) .and. .not.(tke_noflux_bottom_BC)) then
1740	16188	tke(nz+1) = TKE_min
1741	16188	dQ(nz+1) = 0.0
1742		else
1743	19374	k = ke_tke
1744	19374	tke_src = dz_h_Int(K)kappa0S2(K) + q0*TKE_decay(K) ! Uses that kappa(ke_tke) = 0
1745	19374	if (K == nz+1) then
1746	0	dQ(K) = -TKE(K)
1747	0	bQ = 1.0 / (h_Int(K)TKE_decay(K) + cQcompaQ(k-1))
1748	0	tke(K) = max(TKE_min, bQ * (h_Int(K)tke_src + aQ(k-1)tke(K-1)))
1749	0	dQ(K) = tke(K) + dQ(K)
1750		else
1751	19374	bQ = 1.0 / ((h_Int(K)TKE_decay(K) + cQcompaQ(k-1)) + aQ(k))
1752	19374	cQ(K+1) = aQ(k) * bQ
1753		! Account for all changes deeper in the water column.
1754	19374	dQ(K) = -TKE(K)
1755		tke(K) = max((bQ * (h_Int(K)tke_src + aQ(k-1)tke(K-1)) + &
1756	19374	cQ(K+1)(tke(K+1) - e1(K+1)tke(K))) / (1.0 - cQ(K+1)*e1(K+1)), TKE_min)
1757	19374	dQ(K) = tke(K) + dQ(K)
1758
1759		! Adjust TKE deeper in the water column in case ke_tke increases.
1760		! This might not be strictly necessary?
1761	56300	do K=ke_tke+1,nz+1
1762	56300	dQ(K) = e1(K)*dQ(K-1)
1763	56300	tke(K) = max(tke(K) + dQ(K), TKE_min)
1764	56300	if (abs(dQ(K)) < roundoff*tke(K)) exit
1765		enddo
1766	20260	do K2=K+1,nz
1767	1676	if (dQ(K2) == 0.0) exit
1768	19470	dQ(K2) = 0.0
1769		enddo
1770		endif
1771		endif
1772	179755	do K=ke_tke-1,1,-1
1773	144193	tke(K) = max(tke(K) + cQ(K+1)*tke(K+1), TKE_min)
1774	179755	dQ(K) = tke(K) + dQ(K)
1775		enddo
1776
1777		! ----------------------------------------------------
1778		! Calculate kappa, here defined at interfaces.
1779		! ----------------------------------------------------
1780
1781	35562	ke_kappa_prev = ke_kappa ; ks_kappa_prev = ks_kappa
1782
1783	35562	dK(1) = 0.0 ! kappa takes boundary values of 0.
1784	35562	cK(2) = 0.0 ; cKcomp = 1.0
1785	89817	if (itt == 1) then ; do K=2,nz
1786	65007	I_Ld2(K) = dz_h_Int(K)(N2(K)Ilambda2 + f2) / tke(K) + I_L2_bdry(K)
1787		enddo ; endif
1788	77432	do K=2,nz
1789	77430	dK(K) = -kappa(K)
1790	77430	if (itt>1) &
1791	52020	I_Ld2(K) = dz_h_Int(K)(N2(K)Ilambda2 + f2) / tke(K) + I_L2_bdry(K)
1792	77430	bKd1 = h_Int(K)I_Ld2(K) + cKcompIdz(k-1)
1793	77430	bK = 1.0 / (bKd1 + Idz(k))
1794
1795	77430	kappa(K) = bK * (Idz(k-1)kappa(K-1) + h_Int(K) K_src(K))
1796	77430	cK(K+1) = Idz(k) * bK ; cKcomp = bKd1 * bK ! = 1 - cK(K+1)
1797
1798		! Neglect values that are smaller than kappa_trunc.
1799	77432	if (kappa(K) < cKcomp*kappa_trunc) then
1800	35560	kappa(K) = 0.0
1801	35560	if (K > ke_src) then ; ke_kappa = k-1 ; K_Q(K) = 0.0 ; exit ; endif
1802	41870	elseif (kappa(K) < 2.0cKcompkappa_trunc) then
1803	360	kappa(K) = 2.0 * (kappa(K) - cKcomp*kappa_trunc)
1804		endif
1805		enddo
1806	35562	K_Q(ke_kappa) = kappa(ke_kappa) / tke(ke_kappa)
1807	35562	dK(ke_kappa) = dK(ke_kappa) + kappa(ke_kappa)
1808	64407	do K=ke_kappa+2,ke_kappa_prev
1809	64407	dK(K) = -kappa(K) ; kappa(K) = 0.0 ; K_Q(K) = 0.0
1810		enddo
1811	41870	do K=ke_kappa-1,2,-1
1812	6308	kappa(K) = kappa(K) + cK(K+1)*kappa(K+1)
1813		! Neglect values that are smaller than kappa_trunc.
1814	6308	if (kappa(K) <= kappa_trunc) then
1815	0	kappa(K) = 0.0
1816	0	if (K < ks_src) then ; ks_kappa = k+1 ; K_Q(K) = 0.0 ; exit ; endif
1817	6308	elseif (kappa(K) < 2.0*kappa_trunc) then
1818	0	kappa(K) = 2.0 * (kappa(K) - kappa_trunc)
1819		endif
1820
1821	6308	dK(K) = dK(K) + kappa(K)
1822	41870	K_Q(K) = kappa(K) / tke(K)
1823		enddo
1824	35562	do K=ks_kappa_prev,ks_kappa-2 ; kappa(K) = 0.0 ; K_Q(K) = 0.0 ; enddo
1825
1826		else ! do_Newton is .true.
1827		! Once the solutions are close enough, use a Newton's method solver of the
1828		! whole system to accelerate convergence.
1829	0	ks_kappa_prev = ks_kappa ; ke_kappa_prev = ke_kappa ; ke_kappa = nz
1830	0	ks_kappa = 2
1831	0	dK(1) = 0.0 ; cK(2) = 0.0 ; cKcomp = 1.0 ; dKdQ(1) = 0.0
1832	0	aQ(1) = (0.5(kappa(1)+kappa(2))+kappa0) Idz(1)
1833	0	dQdz(1) = 0.5(TKE(1) - TKE(2))Idz(1)
1834	0	if (tke_noflux_top_BC) then
1835		tke_src = h_Int(1) * (kappa0dz_h_Int(1)S2(1) - (TKE(1) - q0)*TKE_decay(1)) - &
1836	0	aQ(1) * (TKE(1) - TKE(2))
1837
1838	0	bQ = 1.0 / (aQ(1) + h_Int(1)*TKE_decay(1))
1839	0	cQ(2) = aQ(1) * bQ
1840	0	cQcomp = (h_Int(1)TKE_decay(1)) bQ ! = 1 - cQ(2)
1841	0	dQmdK(2) = -dQdz(1) * bQ
1842	0	dQ(1) = bQ * tke_src
1843		else
1844	0	dQ(1) = 0.0 ; cQ(2) = 0.0 ; cQcomp = 1.0 ; dQmdK(2) = 0.0
1845		endif
1846	0	do K=2,nz
1847	0	I_Q = 1.0 / TKE(K)
1848	0	I_Ld2(K) = (N2(K)Ilambda2 + f2) dz_h_Int(K)*I_Q + I_L2_bdry(K)
1849
1850		kap_src = h_Int(K) * (K_src(K) - I_Ld2(K)*kappa(K)) + &
1851	0	Idz(k-1)(kappa(K-1)-kappa(K)) - Idz(k)(kappa(K)-kappa(K+1))
1852
1853		! Ensure that the pivot is always positive, and that 0 <= cK <= 1.
1854		! Otherwise do not use Newton's method.
1855	0	decay_term_k = -Idz(k-1)dQmdK(K)dKdQ(K-1) + h_Int(K)*I_Ld2(K)
1856	0	if (decay_term_k < 0.0) then ; abort_Newton = .true. ; exit ; endif
1857	0	bK = 1.0 / (Idz(k) + Idz(k-1)*cKcomp + decay_term_k)
1858
1859	0	cK(K+1) = bK * Idz(k)
1860	0	cKcomp = bK * (Idz(k-1)*cKcomp + decay_term_k) ! = 1-cK(K+1)
1861	0	if (CS%dKdQ_iteration_bug) then
1862		dKdQ(K) = bK * (Idz(k-1)dKdQ(K-1)cQ(K) + &
1863	0	US%m_to_Z(N2(K)Ilambda2 + f2) * I_Q*2 kappa(K) )
1864		else
1865		dKdQ(K) = bK * (Idz(k-1)dKdQ(K-1)cQ(K) + &
1866	0	dz_Int(K)(N2(K)Ilambda2 + f2) * I_Q*2 kappa(K) )
1867		endif
1868	0	dK(K) = bK * (kap_src + Idz(k-1)dK(K-1) + Idz(k-1)dKdQ(K-1)*dQ(K-1))
1869
1870		! Truncate away negligibly small values of kappa.
1871	0	if (dK(K) <= cKcomp*(kappa_trunc - kappa(K))) then
1872	0	dK(K) = -cKcomp*kappa(K)
1873		! if (K > ke_src) then ; ke_kappa = k-1 ; K_Q(K) = 0.0 ; exit ; endif
1874	0	elseif (dK(K) < cKcomp(2.0kappa_trunc - kappa(K))) then
1875	0	dK(K) = 2.0 * dK(K) - cKcomp(2.0kappa_trunc - kappa(K))
1876		endif
1877
1878		! Solve for dQ(K)...
1879	0	aQ(k) = (0.5(kappa(K)+kappa(K+1))+kappa0) Idz(k)
1880	0	dQdz(k) = 0.5(TKE(K) - TKE(K+1))Idz(k)
1881		tke_src = h_Int(K) * ((dz_h_Int(K) * ((kappa(K) + kappa0)S2(K) - kappa(k)N2(K))) - &
1882		(TKE(k) - q0)*TKE_decay(k)) - &
1883	0	(aQ(k) * (TKE(K) - TKE(K+1)) - aQ(k-1) * (TKE(K-1) - TKE(K)))
1884	0	v1 = aQ(k-1) + dQdz(k-1)*dKdQ(K-1)
1885		v2 = (v1dQmdK(K) + dQdz(k-1)cK(K)) + &
1886	0	((dQdz(k-1) - dQdz(k)) + dz_Int(K)*(S2(K) - N2(K)))
1887
1888		! Ensure that the pivot is always positive, and that 0 <= cQ <= 1.
1889		! Otherwise do not use Newton's method.
1890	0	decay_term_Q = h_Int(K)TKE_decay(K) - dQdz(k-1)dKdQ(K-1)cQ(K) - v2dKdQ(K)
1891	0	if (decay_term_Q < 0.0) then ; abort_Newton = .true. ; exit ; endif
1892	0	bQ = 1.0 / (aQ(k) + (cQcomp*aQ(k-1) + decay_term_Q))
1893
1894	0	cQ(K+1) = aQ(k) * bQ
1895	0	cQcomp = (cQcompaQ(k-1) + decay_term_Q) bQ
1896	0	dQmdK(K+1) = (v2 * cK(K+1) - dQdz(k)) * bQ
1897
1898		! Ensure that TKE+dQ will not drop below 0.5*TKE.
1899		dQ(K) = max(bQ * ((v1 * dQ(K-1) + dQdz(k-1)*dK(k-1)) + &
1900	0	(v2 * dK(K) + tke_src)), cQcomp(-0.5TKE(K)))
1901
1902		! Check whether the next layer will be affected by any nonzero kappas.
1903	0	if ((itt > 1) .and. (K > ke_src) .and. (dK(K) == 0.0) .and. &
1904	0	((kappa(K) + kappa(K+1)) == 0.0)) then
1905		! Could also do .and. (bQabs(tke_src) < roundoffTKE(K)) then
1906	0	ke_kappa = k-1 ; exit
1907		endif
1908		enddo
1909	0	if ((ke_kappa == nz) .and. (.not. abort_Newton)) then
1910	0	dK(nz+1) = 0.0 ; dKdQ(nz+1) = 0.0
1911	0	if (tke_noflux_bottom_BC) then
1912	0	K = nz+1
1913		tke_src = h_Int(K) * (kappa0dz_h_Int(K)S2(K) - (TKE(K) - q0)*TKE_decay(K)) + &
1914	0	aQ(k-1) * (TKE(K-1) - TKE(K))
1915
1916	0	v1 = aQ(k-1) + dQdz(k-1)*dKdQ(K-1)
1917	0	decay_term_Q = max(0.0, h_Int(K)TKE_decay(K) - dQdz(k-1)dKdQ(K-1)*cQ(K))
1918	0	if (decay_term_Q < 0.0) then
1919	0	abort_Newton = .true.
1920		else
1921	0	bQ = 1.0 / (aQ(k) + (cQcomp*aQ(k-1) + decay_term_Q))
1922		! Ensure that TKE+dQ will not drop below 0.5*TKE.
1923	0	dQ(K) = max(bQ * ((v1 * dQ(K-1) + dQdz(k-1)dK(K-1)) + tke_src), -0.5TKE(K))
1924	0	TKE(K) = max(TKE(K) + dQ(K), TKE_min)
1925		endif
1926		else
1927	0	dQ(nz+1) = 0.0
1928		endif
1929	0	elseif (.not. abort_Newton) then
1930		! Alter the first-guess determination of dQ(K).
1931	0	dQ(ke_kappa+1) = dQ(ke_kappa+1) / (1.0 - cQ(ke_kappa+2)*e1(ke_kappa+2))
1932	0	TKE(ke_kappa+1) = max(TKE(ke_kappa+1) + dQ(ke_kappa+1), TKE_min)
1933	0	do k=ke_kappa+2,nz+1
1934		if (debug_soln .and. (K < nz+1)) then
1935		! Ignore this source?
1936		aQ(k) = (0.5(kappa(K)+kappa(K+1))+kappa0) Idz(k)
1937		! tke_src_norm = ((kappa0dz_Int(K)S2(K) - h_Int(K)(TKE(K)-q0)TKE_decay(K)) - &
1938		! (aQ(k) * (TKE(K) - TKE(K+1)) - aQ(k-1) * (TKE(K-1) - TKE(K))) ) / &
1939		! (aQ(k) + (aQ(k-1) + h_Int(K)*TKE_decay(K)))
1940		endif
1941	0	dK(K) = 0.0
1942		! Ensure that TKE+dQ will not drop below 0.5*TKE.
1943	0	dQ(K) = max(e1(K)dQ(K-1),-0.5TKE(K))
1944	0	TKE(K) = max(TKE(K) + dQ(K), TKE_min)
1945	0	if (abs(dQ(K)) < roundoff*TKE(K)) exit
1946		enddo
1947		if (debug_soln) then ; do K2=K+1,nz+1 ; dQ(K2) = 0.0 ; dK(K2) = 0.0 ; enddo ; endif
1948		endif
1949	0	if (.not. abort_Newton) then
1950	0	do K=ke_kappa,2,-1
1951		! Ensure that TKE+dQ will not drop below 0.5*TKE.
1952	0	dQ(K) = max(dQ(K) + (cQ(K+1)dQ(K+1) + dQmdK(K+1) dK(K+1)), -0.5*TKE(K))
1953	0	TKE(K) = max(TKE(K) + dQ(K), TKE_min)
1954	0	dK(K) = dK(K) + (cK(K+1)dK(K+1) + dKdQ(K) dQ(K))
1955		! Truncate away negligibly small values of kappa.
1956	0	if (dK(K) <= kappa_trunc - kappa(K)) then
1957	0	dK(K) = -kappa(K)
1958	0	kappa(K) = 0.0
1959	0	if ((K < ks_src) .and. (K+1 > ks_kappa)) ks_kappa = K+1
1960	0	elseif (dK(K) < 2.0*kappa_trunc - kappa(K)) then
1961	0	dK(K) = 2.0dK(K) - (2.0kappa_trunc - kappa(K))
1962	0	kappa(K) = max(kappa(K) + dK(K), 0.0) ! The max is for paranoia.
1963	0	if (K<=ks_kappa) ks_kappa = 2
1964		else
1965	0	kappa(K) = kappa(K) + dK(K)
1966	0	if (K<=ks_kappa) ks_kappa = 2
1967		endif
1968		enddo
1969	0	dQ(1) = max(dQ(1) + cQ(2)dQ(2) + dQmdK(2) dK(2), TKE_min - TKE(1))
1970	0	TKE(1) = max(TKE(1) + dQ(1), TKE_min)
1971	0	dK(1) = 0.0
1972		endif
1973
1974		! Check these solutions for consistency.
1975		! The unit conversions here have not been carefully tested.
1976		if (debug_soln) then ; do K=2,nz
1977		! In these equations, K_err_lin and Q_err_lin should be at round-off levels
1978		! compared with the dominant terms, perhaps, h_IntI_Ld2kappa and
1979		! h_IntTKE_decayTKE. The exception is where, either 1) the decay term has been
1980		! been increased to ensure a positive pivot, or 2) negative TKEs have been
1981		! truncated, or 3) small or negative kappas have been rounded toward 0.
1982		I_Q = 1.0 / TKE(K)
1983		I_Ld2_debug(K) = (N2(K)Ilambda2 + f2) dz_h_Int(K)*I_Q + I_L2_bdry(K)
1984
1985		kap_src = h_Int(K) * (k_src(K) - I_Ld2(K)*kappa_prev(K)) + &
1986		(Idz(k-1)*(kappa_prev(k-1)-kappa_prev(k)) - &
1987		Idz(k)*(kappa_prev(k)-kappa_prev(k+1)))
1988		K_err_lin = -Idz(k-1)(dK(K-1)-dK(K)) + Idz(k)(dK(K)-dK(K+1)) + &
1989		h_Int(K)I_Ld2_debug(K)dK(K) - kap_src - &
1990		dz_Int(K)(N2(K)Ilambda2 + f2)I_Q2kappa_prev(K) * dQ(K)
1991
1992		h_dz_here = 0.0 ; if (abs(dz_h_Int(K)) > 0.0) h_dz_here = 1.0 / dz_h_Int(K)
1993		tke_src = h_Int(K) * ((kappa_prev(K) + kappa0)*S2(K) - &
1994		kappa_prev(K)N2(K) - (TKE_prev(K) - q0)h_dz_here*TKE_decay(K)) - &
1995		(aQ(k) * (TKE_prev(K) - TKE_prev(K+1)) - aQ(k-1) * (TKE_prev(K-1) - TKE_prev(K)))
1996		Q_err_lin = tke_src + (aQ(k-1) * (dQ(K-1)-dQ(K)) - aQ(k) * (dQ(k)-dQ(k+1))) - &
1997		0.5(TKE_prev(K)-TKE_prev(K+1))Idz(k) * (dK(K) + dK(K+1)) - &
1998		0.5(TKE_prev(K)-TKE_prev(K-1))Idz(k-1)* (dK(K-1) + dK(K)) + &
1999		dz_Int(K) * (dK(K) * (S2(K) - N2(K)) - dQ(K)*TKE_decay(K))
2000		enddo ; endif
2001
2002		endif ! End of the Newton's method solver.
2003
2004		! Test kappa for convergence...
2005	35562	if ((tol_err < Newton_err) .and. (.not.abort_Newton)) then
2006		! A lower tolerance is used to switch to Newton's method than to switch back.
2007	10752	Newton_test = Newton_err ; if (do_Newton) Newton_test = 2.0*Newton_err
2008	10752	was_Newton = do_Newton
2009	10752	within_tolerance = .true. ; do_Newton = .true.
2010	36765	do K=min(ks_kappa,ks_kappa_prev),max(ke_kappa,ke_kappa_prev)
2011	31447	kappa_mean = kappa0 + (kappa(K) - 0.5*dK(K))
2012	31447	if (abs(dK(K)) > Newton_test * kappa_mean) then
2013	3584	if (do_Newton) abort_Newton = .true.
2014	3584	within_tolerance = .false. ; do_Newton = .false. ; exit
2015	27863	elseif (abs(dK(K)) > tol_err * kappa_mean) then
2016	1850	within_tolerance = .false. ; if (.not.do_Newton) exit
2017		endif
2018	33181	if (abs(dQ(K)) > Newton_test(tke(K) - 0.5dQ(K))) then
2019	20043	if (do_Newton) abort_Newton = .true.
2020	20043	do_Newton = .false. ; if (.not.within_tolerance) exit
2021		endif
2022		enddo
2023
2024		else ! Newton's method will not be used again, so no need to check.
2025	24810	within_tolerance = .true.
2026	31600	do K=min(ks_kappa,ks_kappa_prev),max(ke_kappa,ke_kappa_prev)
2027	31600	if (abs(dK(K)) > tol_err * (kappa0 + (kappa(K) - 0.5*dK(K)))) then
2028	19376	within_tolerance = .false. ; exit
2029		endif
2030		enddo
2031		endif
2032
2033	35562	if (abort_Newton) then
2034	10752	do_Newton = .false. ; abort_Newton = .false.
2035		! We went to Newton too quickly last time, so restrict the tolerance.
2036	10752	Newton_err = 0.5*Newton_err
2037	10752	ke_kappa_prev = nz
2038	65007	do K=2,nz ; K_Q(K) = kappa(K) / max(TKE(K), TKE_min) ; enddo
2039		endif
2040
2041	35562	if (within_tolerance) exit
2042
2043		enddo
2044
2045	10752	if (do_Newton) then ! K_Q needs to be calculated.
2046	0	do K=1,ks_kappa-1 ; K_Q(K) = 0.0 ; enddo
2047	0	do K=ks_kappa,ke_kappa ; K_Q(K) = kappa(K) / TKE(K) ; enddo
2048	0	do K=ke_kappa+1,nz+1 ; K_Q(K) = 0.0 ; enddo
2049		endif
2050
2051	10752	if (present(local_src)) then
2052	3584	local_src(1) = 0.0 ; local_src(nz+1) = 0.0
2053	21669	do K=2,nz
2054	18085	diffusive_src = Idz(k-1)(kappa(K-1)-kappa(K)) + Idz(k)(kappa(K+1)-kappa(K))
2055	18085	chg_by_k0 = kappa0 * ((Idz(k-1)+Idz(k)) / h_Int(K) + I_Ld2(K))
2056	21669	if (diffusive_src <= 0.0) then
2057	14159	local_src(K) = K_src(K) + chg_by_k0
2058		else
2059	3926	local_src(K) = (K_src(K) + chg_by_k0) + diffusive_src / h_Int(K)
2060		endif
2061		enddo
2062		endif
2063	10752	if (present(kappa_src)) then
2064	3584	kappa_src(1) = 0.0 ; kappa_src(nz+1) = 0.0
2065	21669	do K=2,nz
2066	21669	kappa_src(K) = K_src(K)
2067		enddo
2068		endif
2069
2070	99668	end subroutine find_kappa_tke
2071
2072		!> This subroutine initializes the parameters that regulate shear-driven mixing
2073	1	function kappa_shear_init(Time, G, GV, US, param_file, diag, CS)
2074		type(time_type), intent(in) :: Time !< The current model time.
2075		type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure.
2076		type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure.
2077		type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
2078		type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time
2079		!! parameters.
2080		type(diag_ctrl), target, intent(inout) :: diag !< A structure that is used to regulate diagnostic
2081		!! output.
2082		type(Kappa_shear_CS), pointer :: CS !< A pointer that is set to point to the control
2083		!! structure for this module
2084		logical :: kappa_shear_init !< True if module is to be used, False otherwise
2085
2086		! Local variables
2087		real :: KD_normal ! The KD of the main model, read here only as a parameter
2088		! for setting the default of KD_SMOOTH [Z2 T-1 ~> m2 s-1]
2089		real :: kappa_0_default ! The default value for KD_KAPPA_SHEAR_0 [Z2 T-1 ~> m2 s-1]
2090		logical :: merge_mixedlayer
2091		integer :: number_of_OBC_segments
2092		logical :: debug_shear
2093		logical :: enable_bugs ! If true, the defaults for recently added bug-fix flags are set to
2094		! recreate the bugs, or if false bugs are only used if actively selected.
2095		logical :: just_read ! If true, this module is not used, so only read the parameters.
2096		! This include declares and sets the variable "version".
2097		# include "version_variable.h"
2098		character(len=40) :: mdl = "MOM_kappa_shear" ! This module's name.
2099
2100	1	if (associated(CS)) then
2101		call MOM_error(WARNING, "kappa_shear_init called with an associated "// &
2102	0	"control structure.")
2103	0	return
2104		endif
2105	1	allocate(CS)
2106
2107		! The Jackson-Hallberg-Legg shear mixing parameterization uses the following
2108		! 6 nondimensional coefficients. That paper gives 3 best fit parameter sets.
2109		! Ri_Crit Rate FRi_Curv K_buoy TKE_N TKE_Shear
2110		! p1: 0.25 0.089 -0.97 0.82 0.24 0.14
2111		! p2: 0.30 0.085 -0.94 0.86 0.26 0.13
2112		! p3: 0.35 0.088 -0.89 0.81 0.28 0.12
2113		! Future research will reveal how these should be modified to take
2114		! subgridscale inhomogeneity into account.
2115
2116		! Set default, read and log parameters
2117	1	call get_param(param_file, mdl, "USE_JACKSON_PARAM", kappa_shear_init, default=.false., do_not_log=.true.)
2118		call log_version(param_file, mdl, version, &
2119		"Parameterization of shear-driven turbulence following Jackson, Hallberg and Legg, JPO 2008", &
2120	1	log_to_all=.true., debugging=kappa_shear_init, all_default=.not.kappa_shear_init)
2121		call get_param(param_file, mdl, "USE_JACKSON_PARAM", kappa_shear_init, &
2122		"If true, use the Jackson-Hallberg-Legg (JPO 2008) "//&
2123	1	"shear mixing parameterization.", default=.false.)
2124	1	just_read = .not.kappa_shear_init
2125		call get_param(param_file, mdl, "VERTEX_SHEAR", CS%KS_at_vertex, &
2126		"If true, do the calculations of the shear-driven mixing "//&
2127		"at the cell vertices (i.e., the vorticity points).", &
2128	1	default=.false., do_not_log=just_read)
2129		call get_param(param_file, mdl, "ENABLE_BUGS_BY_DEFAULT", enable_bugs, &
2130	1	default=.true., do_not_log=.true.) ! This is logged from MOM.F90.
2131		call get_param(param_file, mdl, "VERTEX_SHEAR_VISCOSITY_BUG", CS%VS_viscosity_bug, &
2132		"If true, use a bug in vertex shear that zeros out viscosities at "//&
2133		"vertices on coastlines.", &
2134	1	default=enable_bugs, do_not_log=just_read.or.(.not.CS%KS_at_vertex))
2135		call get_param(param_file, mdl, "OBC_NUMBER_OF_SEGMENTS", number_of_OBC_segments, &
2136	1	default=0, do_not_log=.true.)
2137		call get_param(param_file, mdl, "VERTEX_SHEAR_OBC_BUG", CS%vertex_shear_OBC_bug, &
2138		"If false, use extra masking when interpolating thicknesses to velocity "//&
2139		"points for setting up the shear velocities at vertices to avoid using "//&
2140		"external thicknesses at open boundaries. When OBCs are not in use, "//&
2141		"this parameter does not change answers, but true is more efficient.", &
2142		default=enable_bugs, &
2143	1	do_not_log=just_read.or.(.not.CS%KS_at_vertex).or.(number_of_OBC_segments<=0))
2144		! Use OBC settings to set the default for VERTEX_SHEAR_OBC_BUG?
2145		call get_param(param_file, mdl, "VERTEX_SHEAR_GEOMETRIC_MEAN", CS%VS_GeometricMean, &
2146		"If true, use a geometric mean for moving diffusivity from "//&
2147		"vertices to tracer points. False uses algebraic mean.", &
2148	1	default=.false., do_not_log=just_read.or.(.not.CS%KS_at_vertex))
2149		call get_param(param_file, mdl, "VERTEX_SHEAR_THICKNESS_MEAN", CS%VS_ThicknessMean, &
2150		"If true, apply thickness weighting to horizontal averagings of diffusivity "//&
2151		"to tracer points in the kappa shear solver.", &
2152	1	default=.false.)
2153	1	if (CS%VS_GeometricMean) then
2154		call get_param(param_file, mdl, "VERTEX_SHEAR_GEOMETRIC_MEAN_KDMIN", &
2155		CS%VS_GeoMean_Kdmin, "If using the geometric mean in vertex shear, "//&
2156		"use this minimum value for Kd. This is an ad-hoc parameter, the "//&
2157		"diffusivities on the edge of shear regions are sensitive to the choice.",&
2158	0	units="m2 s-1", default=0.0, scale=GV%m2_s_to_HZ_T, do_not_log=just_read)
2159		endif
2160		call get_param(param_file, mdl, "RINO_CRIT", CS%RiNo_crit, &
2161		"The critical Richardson number for shear mixing.", &
2162	1	units="nondim", default=0.25, do_not_log=just_read)
2163		call get_param(param_file, mdl, "SHEARMIX_RATE", CS%Shearmix_rate, &
2164		"A nondimensional rate scale for shear-driven entrainment. "//&
2165		"Jackson et al find values in the range of 0.085-0.089.", &
2166	1	units="nondim", default=0.089, do_not_log=just_read)
2167		call get_param(param_file, mdl, "MAX_RINO_IT", CS%max_RiNo_it, &
2168		"The maximum number of iterations that may be used to "//&
2169		"estimate the Richardson number driven mixing.", &
2170	1	units="nondim", default=50, do_not_log=just_read)
2171		call get_param(param_file, mdl, "KD", KD_normal, &
2172	1	units="m2 s-1", scale=US%m2_s_to_Z2_T, default=0.0, do_not_log=.true.)
2173	1	kappa_0_default = max(Kd_normal, 1.0e-7*US%m2_s_to_Z2_T)
2174		call get_param(param_file, mdl, "KD_KAPPA_SHEAR_0", CS%kappa_0, &
2175		"The background diffusivity that is used to smooth the "//&
2176		"density and shear profiles before solving for the "//&
2177		"diffusivities. The default is the greater of KD and 1e-7 m2 s-1.", &
2178		units="m2 s-1", default=kappa_0_default*US%Z2_T_to_m2_s, scale=GV%m2_s_to_HZ_T, &
2179	1	do_not_log=just_read)
2180		call get_param(param_file, mdl, "KD_SEED_KAPPA_SHEAR", CS%kappa_seed, &
2181		"A moderately large seed value of diapycnal diffusivity that is used as a "//&
2182		"starting turbulent diffusivity in the iterations to find an energetically "//&
2183		"constrained solution for the shear-driven diffusivity.", &
2184	1	units="m2 s-1", default=1.0, scale=GV%m2_s_to_HZ_T)
2185		call get_param(param_file, mdl, "KD_TRUNC_KAPPA_SHEAR", CS%kappa_trunc, &
2186		"The value of shear-driven diffusivity that is considered negligible "//&
2187		"and is rounded down to 0. The default is 1% of KD_KAPPA_SHEAR_0.", &
2188		units="m2 s-1", default=0.01CS%kappa_0GV%HZ_T_to_m2_s, scale=GV%m2_s_to_HZ_T, &
2189	1	do_not_log=just_read)
2190		call get_param(param_file, mdl, "FRI_CURVATURE", CS%FRi_curvature, &
2191		"The nondimensional curvature of the function of the "//&
2192		"Richardson number in the kappa source term in the "//&
2193	1	"Jackson et al. scheme.", units="nondim", default=-0.97, do_not_log=just_read)
2194		call get_param(param_file, mdl, "TKE_N_DECAY_CONST", CS%C_N, &
2195		"The coefficient for the decay of TKE due to "//&
2196		"stratification (i.e. proportional to N*tke). "//&
2197		"The values found by Jackson et al. are 0.24-0.28.", &
2198	1	units="nondim", default=0.24, do_not_log=just_read)
2199		! call get_param(param_file, mdl, "LAYER_KAPPA_STAGGER", CS%layer_stagger, &
2200		! default=.false., do_not_log=just_read)
2201		call get_param(param_file, mdl, "TKE_SHEAR_DECAY_CONST", CS%C_S, &
2202		"The coefficient for the decay of TKE due to shear (i.e. "//&
2203		"proportional to \|S\|*tke). The values found by Jackson "//&
2204	1	"et al. are 0.14-0.12.", units="nondim", default=0.14, do_not_log=just_read)
2205		call get_param(param_file, mdl, "KAPPA_BUOY_SCALE_COEF", CS%lambda, &
2206		"The coefficient for the buoyancy length scale in the "//&
2207		"kappa equation. The values found by Jackson et al. are "//&
2208	1	"in the range of 0.81-0.86.", units="nondim", default=0.82, do_not_log=just_read)
2209		call get_param(param_file, mdl, "KAPPA_N_OVER_S_SCALE_COEF2", CS%lambda2_N_S, &
2210		"The square of the ratio of the coefficients of the "//&
2211		"buoyancy and shear scales in the diffusivity equation, "//&
2212		"Set this to 0 (the default) to eliminate the shear scale. "//&
2213		"This is only used if USE_JACKSON_PARAM is true.", &
2214	1	units="nondim", default=0.0, do_not_log=just_read)
2215		call get_param(param_file, mdl, "LZ_RESCALE", CS%lz_rescale, &
2216		"A coefficient to rescale the distance to the nearest solid boundary. "//&
2217		"This adjustment is to account for regions where 3 dimensional turbulence "//&
2218		"prevents the growth of shear instabilities [nondim].", &
2219	1	units="nondim", default=1.0)
2220		call get_param(param_file, mdl, "KAPPA_SHEAR_TOL_ERR", CS%kappa_tol_err, &
2221		"The fractional error in kappa that is tolerated. "//&
2222		"Iteration stops when changes between subsequent "//&
2223		"iterations are smaller than this everywhere in a "//&
2224		"column. The peak diffusivities usually converge most "//&
2225		"rapidly, and have much smaller errors than this.", &
2226	1	units="nondim", default=0.1, do_not_log=just_read)
2227		call get_param(param_file, mdl, "TKE_BACKGROUND", CS%TKE_bg, &
2228		"A background level of TKE used in the first iteration "//&
2229		"of the kappa equation. TKE_BACKGROUND could be 0.", &
2230	1	units="m2 s-2", default=0.0, scale=US%m_to_Z*2US%T_to_s**2)
2231		call get_param(param_file, mdl, "KAPPA_SHEAR_ELIM_MASSLESS", CS%eliminate_massless, &
2232		"If true, massless layers are merged with neighboring "//&
2233		"massive layers in this calculation. The default is "//&
2234		"true and I can think of no good reason why it should "//&
2235		"be false. This is only used if USE_JACKSON_PARAM is true.", &
2236	1	default=.true., do_not_log=just_read)
2237		call get_param(param_file, mdl, "MAX_KAPPA_SHEAR_IT", CS%max_KS_it, &
2238		"The maximum number of iterations that may be used to "//&
2239		"estimate the time-averaged diffusivity.", &
2240	1	default=13, do_not_log=just_read)
2241		call get_param(param_file, mdl, "PRANDTL_TURB", CS%Prandtl_turb, &
2242		"The turbulent Prandtl number applied to shear instability.", &
2243	1	units="nondim", default=1.0, do_not_log=just_read)
2244		call get_param(param_file, mdl, "VEL_UNDERFLOW", CS%vel_underflow, &
2245		"A negligibly small velocity magnitude below which velocity components are set "//&
2246		"to 0. A reasonable value might be 1e-30 m/s, which is less than an "//&
2247		"Angstrom divided by the age of the universe.", &
2248	1	units="m s-1", default=0.0, scale=US%m_s_to_L_T, do_not_log=just_read)
2249		call get_param(param_file, mdl, "KAPPA_SHEAR_MAX_KAP_SRC_CHG", CS%kappa_src_max_chg, &
2250		"The maximum permitted increase in the kappa source within an iteration relative "//&
2251		"to the local source; this must be greater than 1. The lower limit for the "//&
2252		"permitted fractional decrease is (1 - 0.5/kappa_src_max_chg). These limits "//&
2253		"could perhaps be made dynamic with an improved iterative solver.", &
2254	1	default=10.0, units="nondim", do_not_log=just_read)
2255
2256		call get_param(param_file, mdl, "DEBUG", CS%debug, &
2257		"If true, write out verbose debugging data.", &
2258	1	default=.false., debuggingParam=.true., do_not_log=just_read)
2259		call get_param(param_file, mdl, "DEBUG_KAPPA_SHEAR", debug_shear, &
2260		"If true, write debugging data for the kappa-shear code.", &
2261	1	default=.false., debuggingParam=.true., do_not_log=.true.)
2262	1	if (debug_shear) CS%debug = .true.
2263		call get_param(param_file, mdl, "KAPPA_SHEAR_VERTEX_PSURF_BUG", CS%psurf_bug, &
2264		"If true, do a simple average of the cell surface pressures to get a pressure "//&
2265		"at the corner if VERTEX_SHEAR=True. Otherwise mask out any land points in "//&
2266	1	"the average.", default=.false., do_not_log=(just_read .or. (.not.CS%KS_at_vertex)))
2267
2268		call get_param(param_file, mdl, "KAPPA_SHEAR_ITER_BUG", CS%dKdQ_iteration_bug, &
2269		"If true, use an older, dimensionally inconsistent estimate of the "//&
2270		"derivative of diffusivity with energy in the Newton's method iteration. "//&
2271	1	"The bug causes undercorrections when dz > 1 m.", default=.false., do_not_log=just_read)
2272		call get_param(param_file, mdl, "KAPPA_SHEAR_ALL_LAYER_TKE_BUG", CS%all_layer_TKE_bug, &
2273		"If true, report back the latest estimate of TKE instead of the time average "//&
2274		"TKE when there is mass in all layers. Otherwise always report the time "//&
2275		"averaged TKE, as is currently done when there are some massless layers.", &
2276	1	default=.false., do_not_log=just_read)
2277		call get_param(param_file, mdl, "USE_RESTRICTIVE_TOLERANCE_CHECK", CS%restrictive_tolerance_check, &
2278		"If true, uses the more restrictive tolerance check to determine if a timestep "//&
2279		"is acceptable for the KS_it outer iteration loop. False uses the original less "//&
2280	1	"restrictive check.", default=.false., do_not_log=just_read)
2281		! id_clock_KQ = cpu_clock_id('Ocean KS kappa_shear', grain=CLOCK_ROUTINE)
2282		! id_clock_avg = cpu_clock_id('Ocean KS avg', grain=CLOCK_ROUTINE)
2283		! id_clock_project = cpu_clock_id('Ocean KS project', grain=CLOCK_ROUTINE)
2284		! id_clock_setup = cpu_clock_id('Ocean KS setup', grain=CLOCK_ROUTINE)
2285
2286	1	CS%nkml = 1
2287	1	if (GV%nkml>0) then
2288		call get_param(param_file, mdl, "KAPPA_SHEAR_MERGE_ML",merge_mixedlayer, &
2289		"If true, combine the mixed layers together before solving the "//&
2290	1	"kappa-shear equations.", default=.true., do_not_log=just_read)
2291	1	if (merge_mixedlayer) CS%nkml = GV%nkml
2292		endif
2293
2294		! Forego remainder of initialization if not using this scheme
2295	1	if (.not. kappa_shear_init) return
2296
2297	1	CS%diag => diag
2298
2299		CS%id_Kd_shear = register_diag_field('ocean_model','Kd_shear', diag%axesTi, Time, &
2300	1	'Shear-driven Diapycnal Diffusivity at horizontal tracer points', 'm2 s-1', conversion=GV%HZ_T_to_m2_s)
2301	1	if (CS%KS_at_vertex) then
2302		CS%id_TKE = register_diag_field('ocean_model','TKE_shear', diag%axesBi, Time, &
2303	0	'Shear-driven Turbulent Kinetic Energy at horizontal vertices', 'm2 s-2', conversion=US%Z_to_m*2US%s_to_T**2)
2304		CS%id_Kd_vertex = register_diag_field('ocean_model','Kd_shear_vertex', diag%axesBi, Time, &
2305	0	'Shear-driven Diapycnal Diffusivity at horizontal vertices', 'm2 s-1', conversion=GV%HZ_T_to_m2_s)
2306		CS%id_S2_init = register_diag_field('ocean_model','S2_shear_in', diag%axesBi, Time, &
2307	0	'Interface shear squared at horizontal vertices, as input to kappa-shear', 's-2', conversion=US%s_to_T**2)
2308		CS%id_N2_init = register_diag_field('ocean_model','N2_shear_in', diag%axesBi, Time, &
2309	0	'Interface stratification at horizontal vertices, as input to kappa-shear', 's-2', conversion=US%s_to_T**2)
2310		CS%id_S2_mean = register_diag_field('ocean_model','S2_shear_mean', diag%axesBi, Time, &
2311		'Interface shear squared at horizontal vertices, averaged over timestep in kappa-shear', &
2312	0	's-2', conversion=US%s_to_T**2)
2313		CS%id_N2_mean = register_diag_field('ocean_model','N2_shear_mean', diag%axesBi, Time, &
2314		'Interface stratification at horizontal vertices, averaged over timestep in kappa-shear', &
2315	0	's-2', conversion=US%s_to_T**2)
2316		else
2317		CS%id_TKE = register_diag_field('ocean_model','TKE_shear', diag%axesTi, Time, &
2318		'Shear-driven Turbulent Kinetic Energy at horizontal tracer points', &
2319	1	'm2 s-2', conversion=US%Z_to_m*2US%s_to_T**2)
2320		CS%id_S2_init = register_diag_field('ocean_model','S2_shear_in', diag%axesTi, Time, &
2321	1	'Interface shear squared at horizontal tracer points, as input to kappa-shear', 's-2', conversion=US%s_to_T**2)
2322		CS%id_N2_init = register_diag_field('ocean_model','N2_shear_in', diag%axesTi, Time, &
2323		'Interface stratification at horizontal tracer points, as input to kappa-shear', &
2324	1	's-2', conversion=US%s_to_T**2)
2325		CS%id_S2_mean = register_diag_field('ocean_model','S2_shear_mean', diag%axesTi, Time, &
2326		'Interface shear squared at horizontal tracer points, averaged over timestep in kappa-shear', &
2327	1	's-2', conversion=US%s_to_T**2)
2328		CS%id_N2_mean = register_diag_field('ocean_model','N2_shear_mean', diag%axesTi, Time, &
2329		'Interface stratification at horizontal tracer points, averaged ove timestep in kappa-shear', &
2330	1	's-2', conversion=US%s_to_T**2)
2331		endif
2332
2333	1	end function kappa_shear_init
2334
2335		!> This function indicates to other modules whether the Jackson et al shear mixing
2336		!! parameterization will be used without needing to duplicate the log entry.
2337	4	logical function kappa_shear_is_used(param_file)
2338		type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters
2339
2340		! Local variables
2341		character(len=40) :: mdl = "MOM_kappa_shear" ! This module's name.
2342		! This function reads the parameter "USE_JACKSON_PARAM" and returns its value.
2343
2344		call get_param(param_file, mdl, "USE_JACKSON_PARAM", kappa_shear_is_used, &
2345	4	default=.false., do_not_log=.true.)
2346	4	end function kappa_shear_is_used
2347
2348		!> This function indicates to other modules whether the Jackson et al shear mixing parameterization
2349		!! will be used at the vertices without needing to duplicate the log entry. It returns false if
2350		!! the Jackson et al scheme is not used or if it is used via calculations at the tracer points.
2351	4	logical function kappa_shear_at_vertex(param_file)
2352		type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters
2353
2354		! Local variables
2355		character(len=40) :: mdl = "MOM_kappa_shear" ! This module's name.
2356		logical :: do_kappa_shear
2357		! This function returns true only if the parameters "USE_JACKSON_PARAM" and "VERTEX_SHEAR" are both true.
2358
2359	4	kappa_shear_at_vertex = .false.
2360
2361		call get_param(param_file, mdl, "USE_JACKSON_PARAM", do_kappa_shear, &
2362	4	default=.false., do_not_log=.true.)
2363	4	if (do_Kappa_Shear) &
2364		call get_param(param_file, mdl, "VERTEX_SHEAR", kappa_shear_at_vertex, &
2365		"If true, do the calculations of the shear-driven mixing "//&
2366		"at the cell vertices (i.e., the vorticity points).", &
2367	4	default=.false., do_not_log=.true.)
2368
2369	4	end function kappa_shear_at_vertex
2370
2371		!> \namespace mom_kappa_shear
2372		!!
2373		!! By Laura Jackson and Robert Hallberg, 2006-2008
2374		!!
2375		!! This file contains the subroutines that determine the diapycnal
2376		!! diffusivity driven by resolved shears, as specified by the
2377		!! parameterizations described in Jackson and Hallberg (JPO, 2008).
2378		!!
2379		!! The technique by which the 6 equations (for kappa, TKE, u, v, T,
2380		!! and S) are solved simultaneously has been dramatically revised
2381		!! from the previous version. The previous version was not converging
2382		!! in some cases, especially near the surface mixed layer, while the
2383		!! revised version does. The revised version solves for kappa and
2384		!! TKE with shear and stratification fixed, then marches the density
2385		!! and velocities forward with an adaptive (and aggressive) time step
2386		!! in a predictor-corrector-corrector emulation of a trapezoidal
2387		!! scheme. Run-time-settable parameters determine the tolerance to
2388		!! which the kappa and TKE equations are solved and the minimum time
2389		!! step that can be taken.
2390
2391	0	end module MOM_kappa_shear