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	0	subroutine Calculate_kappa_shear(u_in, v_in, h, tv, p_surf, kappa_io, tke_io, &
136	0	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	0	diag_N2_init, & ! Diagnostic of N2 as provided to this routine [T-2 ~> s-2]
170	0	diag_S2_init, & ! Diagnostic of S2 as provided to this routine [T-2 ~> s-2]
171	0	diag_N2_mean, & ! Diagnostic of N2 averaged over the timestep applied [T-2 ~> s-2]
172	0	diag_S2_mean ! Diagnostic of S2 averaged over the timestep applied [T-2 ~> s-2]
173		real, dimension(SZI_(G),SZK_(GV)) :: &
174	0	h_2d, & ! A 2-D version of h [H ~> m or kg m-2].
175	0	dz_2d, & ! Vertical distance between interface heights [Z ~> m].
176	0	u_2d, v_2d, & ! 2-D versions of u_in and v_in, converted to [L T-1 ~> m s-1].
177	0	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	0	kappa_2d, & ! 2-D version of kappa_io [H Z T-1 ~> m2 s-1 or Pa s]
180	0	tke_2d ! 2-D version tke_io [Z2 T-2 ~> m2 s-2].
181		real, dimension(SZK_(GV)) :: &
182	0	Idz, & ! The inverse of the thickness of the merged layers [H-1 ~> m2 kg-1].
183	0	h_lay, & ! The layer thickness [H ~> m or kg m-2]
184	0	dz_lay, & ! The geometric layer thickness in height units [Z ~> m]
185	0	u0xdz, & ! The initial zonal velocity times dz [H L T-1 ~> m2 s-1 or kg m-1 s-1]
186	0	v0xdz, & ! The initial meridional velocity times dz [H L T-1 ~> m2 s-1 or kg m-1 s-1]
187	0	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	0	S0xdz ! The initial salinity times dz [S H ~> ppt m or ppt kg m-2].
190		real, dimension(SZK_(GV)+1) :: &
191	0	kappa, & ! The shear-driven diapycnal diffusivity at an interface [H Z T-1 ~> m2 s-1 or Pa s]
192	0	tke, & ! The Turbulent Kinetic Energy per unit mass at an interface [Z2 T-2 ~> m2 s-2].
193	0	kappa_avg, & ! The time-weighted average of kappa [H Z T-1 ~> m2 s-1 or Pa s]
194	0	tke_avg, & ! The time-weighted average of TKE [Z2 T-2 ~> m2 s-2].
195	0	N2_init, & ! N2 as provided to this routine [T-2 ~> s-2].
196	0	S2_init, & ! S2 as provided to this routine [T-2 ~> s-2].
197	0	N2_mean, & ! The time-weighted average of N2 [T-2 ~> s-2].
198	0	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	0	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	0	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	0	is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke
217
218	0	use_temperature = associated(tv%T)
219
220	0	k0dt = dt*CS%kappa_0
221	0	dz_massless = 0.1sqrt((US%Z_to_mGV%m_to_H)*k0dt)
222
223	0	if ((CS%id_N2_init>0) .or. CS%debug) diag_N2_init(:,:,:) = 0.0
224	0	if ((CS%id_S2_init>0) .or. CS%debug) diag_S2_init(:,:,:) = 0.0
225	0	if (CS%id_N2_mean>0) diag_N2_mean(:,:,:) = 0.0
226	0	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	0	do j=js,je
232
233		! Convert layer thicknesses into geometric thickness in height units.
234	0	call thickness_to_dz(h, tv, dz_2d, j, G, GV)
235
236	0	do k=1,nz ; do i=is,ie
237	0	h_2d(i,k) = h(i,j,k)
238	0	u_2d(i,k) = u_in(i,j,k) ; v_2d(i,k) = v_in(i,j,k)
239		enddo ; enddo
240	0	if (use_temperature) then ; do k=1,nz ; do i=is,ie
241	0	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	0	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	0	if (CS%eliminate_massless) then
255	0	nzc = 1
256	0	do k=1,nz
257		! Zero out the thicknesses of all layers, even if they are unused.
258	0	h_lay(k) = 0.0 ; dz_lay(k) = 0.0 ; u0xdz(k) = 0.0 ; v0xdz(k) = 0.0
259	0	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	0	if ((k>CS%nkml) .and. (h_lay(nzc) > 0.0) .and. &
264	0	(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	0	kc(k) = nzc
271	0	h_lay(nzc) = h_lay(nzc) + h_2d(i,k)
272	0	dz_lay(nzc) = dz_lay(nzc) + dz_2d(i,k)
273	0	u0xdz(nzc) = u0xdz(nzc) + u_2d(i,k)*h_2d(i,k)
274	0	v0xdz(nzc) = v0xdz(nzc) + v_2d(i,k)*h_2d(i,k)
275	0	if (use_temperature) then
276	0	T0xdz(nzc) = T0xdz(nzc) + T_2d(i,k)*h_2d(i,k)
277	0	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	0	kc(nz+1) = nzc+1
284
285		! Set up Idz as the inverse of layer thicknesses.
286	0	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	0	kf(1) = 0.0 ; dz_in_lay = h_2d(i,1)
291	0	do k=2,nz
292	0	if (kc(k) > kc(k-1)) then
293	0	kf(k) = 0.0 ; dz_in_lay = h_2d(i,k)
294		else
295	0	kf(k) = dz_in_lay*Idz(kc(k)) ; dz_in_lay = dz_in_lay + h_2d(i,k)
296		endif
297		enddo
298	0	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	0	(G%Coriolis2Bu(I,J-1) + G%Coriolis2Bu(I-1,J)))
319	0	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	0	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	0	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	0	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	0	do K=1,nz+1
357	0	if (kf(K) == 0.0) then
358	0	kappa_2d(i,K) = kappa_avg(kc(K))
359	0	tke_2d(i,K) = tke_avg(kc(K))
360		else
361	0	kappa_2d(i,K) = (1.0-kf(K)) * kappa_avg(kc(K)) + kf(K) * kappa_avg(kc(K)+1)
362	0	tke_2d(i,K) = (1.0-kf(K)) * tke_avg(kc(K)) + kf(K) * tke_avg(kc(K)+1)
363		endif
364		enddo
365	0	do K=1,nz+1
366	0	if (kf(K) == 0.0) then
367	0	if (CS%id_N2_mean>0) diag_N2_mean(i,j,K) = N2_mean(kc(K))
368	0	if (CS%id_S2_mean>0) diag_S2_mean(i,j,K) = S2_mean(kc(K))
369	0	if ((CS%id_N2_init>0) .or. CS%debug) diag_N2_init(i,j,K) = N2_init(kc(K))
370	0	if ((CS%id_S2_init>0) .or. CS%debug) diag_S2_init(i,j,K) = S2_init(kc(K))
371		else
372	0	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	0	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	0	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	0	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	0	do K=1,nz+1
386	0	kappa_2d(i,K) = 0.0 ; tke_2d(i,K) = 0.0
387		enddo
388		endif ; enddo ! i-loop
389
390	0	do K=1,nz+1 ; do i=is,ie
391	0	kappa_io(i,j,K) = G%mask2dT(i,j) * kappa_2d(i,K)
392	0	tke_io(i,j,K) = G%mask2dT(i,j) * tke_2d(i,K)
393	0	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	0	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	0	if (CS%id_Kd_shear > 0) call post_data(CS%id_Kd_shear, kappa_io, CS%diag)
407	0	if (CS%id_TKE > 0) call post_data(CS%id_TKE, tke_io, CS%diag)
408	0	if (CS%id_N2_init > 0) call post_data(CS%id_N2_init, diag_N2_init, CS%diag)
409	0	if (CS%id_S2_init > 0) call post_data(CS%id_S2_init, diag_S2_init, CS%diag)
410	0	if (CS%id_N2_mean > 0) call post_data(CS%id_N2_mean, diag_N2_mean, CS%diag)
411	0	if (CS%id_S2_mean > 0) call post_data(CS%id_S2_mean, diag_S2_mean, CS%diag)
412
413	0	end subroutine Calculate_kappa_shear
414
415
416		!> Subroutine for calculating shear-driven diffusivity and TKE in corner columns
417	12	subroutine Calc_kappa_shear_vertex(u_in, v_in, h, T_in, S_in, tv, p_surf, kappa_io, tke_io, &
418	12	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	24	diag_N2_init, & ! Diagnostic of N2 as provided to this routine [T-2 ~> s-2]
456	24	diag_S2_init, & ! Diagnostic of S2 as provided to this routine [T-2 ~> s-2]
457	24	diag_N2_mean, & ! Diagnostic of N2 averaged over the timestep applied [T-2 ~> s-2]
458	24	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	24	dz_3d ! Vertical distance between interface heights [Z ~> m].
462		real, dimension(SZIB_(G),SZJB_(G),SZK_(GV)+1) :: &
463	24	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	24	h_vert ! Thicknesses interpolated to vertices [H ~> m or kg m-2]
466		real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)) :: &
467	24	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	24	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	24	h_2d, & ! A 2-D version of h interpolated to vertices [H ~> m or kg m-2].
472	24	dz_2d, & ! Vertical distance between interface heights [Z ~> m].
473	24	u_2d, v_2d, & ! 2-D versions of u_in and v_in, converted to [L T-1 ~> m s-1].
474	24	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	24	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	24	tke_2d ! 2-D version tke_io [Z2 T-2 ~> m2 s-2].
479		real, dimension(SZK_(GV)) :: &
480	24	Idz, & ! The inverse of the thickness of the merged layers [H-1 ~> m2 kg-1].
481	24	h_lay, & ! The layer thickness [H ~> m or kg m-2]
482	24	dz_lay, & ! The geometric layer thickness in height units [Z ~> m]
483	24	u0xdz, & ! The initial zonal velocity times dz [L H T-1 ~> m2 s-1 or kg m-1 s-1].
484	24	v0xdz, & ! The initial meridional velocity times dz [H L T-1 ~> m2 s-1 or kg m-1 s-1]
485	24	T0xdz, & ! The initial temperature times dz [C H ~> degC m or degC kg m-2]
486	24	S0xdz ! The initial salinity times dz [S H ~> ppt m or ppt kg m-2]
487		real, dimension(SZK_(GV)+1) :: &
488	24	kappa, & ! The shear-driven diapycnal diffusivity at an interface [H Z T-1 ~> m2 s-1 or Pa s]
489	24	tke, & ! The Turbulent Kinetic Energy per unit mass at an interface [Z2 T-2 ~> m2 s-2].
490	24	kappa_avg, & ! The time-weighted average of kappa [H Z T-1 ~> m2 s-1 or Pa s]
491	24	tke_avg, & ! The time-weighted average of TKE [Z2 T-2 ~> m2 s-2].
492	24	N2_init, & ! N2 as provided to this routine [T-2 ~> s-2].
493	24	S2_init, & ! S2 as provided to this routine [T-2 ~> s-2].
494	24	N2_mean, & ! The time-weighted average of N2 [T-2 ~> s-2].
495	24	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	24	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	24	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	12	isB = G%isc-1 ; ieB = G%iecB ; jsB = G%jsc-1 ; jeB = G%jecB ; nz = GV%ke
524
525	12	if ((CS%id_N2_init>0) .or. CS%debug) diag_N2_init(:,:,:) = 0.0
526	12	if ((CS%id_S2_init>0) .or. CS%debug) diag_S2_init(:,:,:) = 0.0
527	12	if (CS%id_N2_mean>0) diag_N2_mean(:,:,:) = 0.0
528	12	if (CS%id_S2_mean>0) diag_S2_mean(:,:,:) = 0.0
529	8181564	kappa_vertex(:,:,:) = 0.0
530
531	12	use_temperature = associated(tv%T)
532
533	12	k0dt = dt*CS%kappa_0
534	12	dz_massless = 0.1sqrt((US%Z_to_mGV%m_to_H)*k0dt)
535	12	I_Prandtl = 0.0 ; if (CS%Prandtl_turb > 0.0) I_Prandtl = 1.0 / CS%Prandtl_turb
536	12	H_tiny = 0.5 * GV%H_subroundoff
537
538		! Convert layer thicknesses into geometric thickness in height units.
539	12	call thickness_to_dz(h, tv, dz_3d, G, GV, US, halo_size=1)
540
541	12	if (CS%vertex_shear_OBC_bug) then
542		!$OMP parallel do default(shared)
543	912	do k=1,nz
544	6808500	do j=JsB,JeB+1 ; do I=IsB,IeB
545	6807600	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	6753612	do J=JsB,JeB ; do i=IsB,IeB+1
548	6752700	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	744	do J=JsB,JeB
574
575		! Interpolate the various quantities to the corners, using masks.
576	6698532	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	6642900	( (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	6642900	( (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	6642900	GV%H_subroundoff)
585	6642900	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	6642900	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	6642900	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	6642900	(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	6697800	(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	732	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	89304	do I=IsB,IeB ; if ((G%mask2dCu(I,j) + G%mask2dCu(I,j+1)) + &
615	732	(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	61788	if (CS%eliminate_massless) then
620	61788	nzc = 1
621	4695888	do k=1,nz
622		! Zero out the thicknesses of all layers, even if they are unused.
623	4634100	h_lay(k) = 0.0 ; dz_lay(k) = 0.0 ; u0xdz(k) = 0.0 ; v0xdz(k) = 0.0
624	4634100	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	4634100	if ((k>CS%nkml) .and. (h_lay(nzc) > 0.0) .and. &
629	3222012	(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	4634100	kc(k) = nzc
636	4634100	h_lay(nzc) = h_lay(nzc) + h_2d(I,k)
637	4634100	dz_lay(nzc) = dz_lay(nzc) + dz_2d(I,k)
638	4634100	u0xdz(nzc) = u0xdz(nzc) + u_2d(I,k)*h_2d(I,k)
639	4634100	v0xdz(nzc) = v0xdz(nzc) + v_2d(I,k)*h_2d(I,k)
640	4695888	if (use_temperature) then
641	4634100	T0xdz(nzc) = T0xdz(nzc) + T_2d(I,k)*h_2d(I,k)
642	4634100	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	61788	kc(nz+1) = nzc+1
649
650		! Set up Idz as the inverse of layer thicknesses.
651	3345588	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	61788	kf(1) = 0.0 ; dz_in_lay = h_2d(I,1)
656	4634100	do k=2,nz
657	4634100	if (kc(k) > kc(k-1)) then
658	3222012	kf(k) = 0.0 ; dz_in_lay = h_2d(I,k)
659		else
660	1350300	kf(k) = dz_in_lay*Idz(kc(k)) ; dz_in_lay = dz_in_lay + h_2d(I,k)
661		endif
662		enddo
663	61788	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	61788	f2 = G%Coriolis2Bu(I,J)
684	61788	surface_pres = 0.0
685	61788	if (associated(p_surf)) then
686	61788	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	61788	(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	3407376	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	61788	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	61788	if (nz == nzc) then
709	97944	do K=1,nz+1
710	96672	kappa_2d(I,K) = kappa_avg(K)
711	97944	if (CS%all_layer_TKE_bug) then
712	0	tke_2d(I,K) = tke(K)
713		else
714	96672	tke_2d(I,K) = tke_avg(K)
715		endif
716		enddo
717	1272	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	1272	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	1272	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	1272	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	4659732	do K=1,nz+1
731	4659732	if (kf(K) == 0.0) then
732	3248916	kappa_2d(I,K) = kappa_avg(kc(K))
733	3248916	tke_2d(I,K) = tke_avg(kc(K))
734		else
735	1350300	kappa_2d(I,K) = (1.0-kf(K)) * kappa_avg(kc(K)) + kf(K) * kappa_avg(kc(K)+1)
736	1350300	tke_2d(I,K) = (1.0-kf(K)) * tke_avg(kc(K)) + kf(K) * tke_avg(kc(K)+1)
737		endif
738		enddo
739	4659732	do K=1,nz+1
740	4659732	if (kf(K) == 0.0) then
741	3248916	if (CS%id_N2_mean>0) diag_N2_mean(I,J,K) = N2_mean(kc(K))
742	3248916	if (CS%id_S2_mean>0) diag_S2_mean(I,J,K) = S2_mean(kc(K))
743	3248916	if ((CS%id_N2_init>0) .or. CS%debug) diag_N2_init(I,J,K) = N2_init(kc(K))
744	3248916	if ((CS%id_S2_init>0) .or. CS%debug) diag_S2_init(I,J,K) = S2_init(kc(K))
745		else
746	1350300	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	1350300	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	1350300	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	1350300	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	2062368	do K=1,nz+1
760	2062368	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	732	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	744	if (CS%VS_viscosity_bug) then
771	6787836	do K=1,nz+1 ; do I=IsB,IeB
772	6731472	kappa_vertex(I,J,K) = kappa_2d(I,K)
773	6731472	tke_io(I,J,K) = G%mask2dBu(I,J) * tke_2d(I,K)
774	6787104	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	87132	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	86400	kappa_io(i,j,1) = 0.0
791	87120	kappa_io(i,j,nz+1) = 0.0
792		enddo ; enddo
793	12	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	12	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	12	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	6447780	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	6446880	(kappa_vertex(I-1,J,K) + kappa_vertex(I,J-1,K)))
845		enddo ; enddo ; enddo
846		endif
847
848	12	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	12	if (CS%id_Kd_shear > 0) call post_data(CS%id_Kd_shear, kappa_io, CS%diag)
856	12	if (CS%id_TKE > 0) call post_data(CS%id_TKE, tke_io, CS%diag)
857	12	if (CS%id_Kd_vertex > 0) call post_data(CS%id_Kd_vertex, kappa_vertex, CS%diag)
858	12	if (CS%id_N2_init > 0) call post_data(CS%id_N2_init, diag_N2_init, CS%diag)
859	12	if (CS%id_S2_init > 0) call post_data(CS%id_S2_init, diag_S2_init, CS%diag)
860	12	if (CS%id_N2_mean > 0) call post_data(CS%id_N2_mean, diag_N2_mean, CS%diag)
861	12	if (CS%id_S2_mean > 0) call post_data(CS%id_S2_mean, diag_S2_mean, CS%diag)
862
863	12	end subroutine Calc_kappa_shear_vertex
864
865
866		!> This subroutine calculates shear-driven diffusivity and TKE in a single column
867	61788	subroutine kappa_shear_column(kappa, tke, dt, nzc, f2, surface_pres, hlay, dz_lay, &
868	61788	u0xdz, v0xdz, T0xdz, S0xdz, kappa_avg, tke_avg, N2_init, S2_init, &
869	61788	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	123576	u, & ! The zonal velocity after a timestep of mixing [L T-1 ~> m s-1].
915	123576	v, & ! The meridional velocity after a timestep of mixing [L T-1 ~> m s-1].
916	123576	Idz, & ! The inverse of the distance between TKE points [Z-1 ~> m-1].
917	123576	T, & ! The potential temperature after a timestep of mixing [C ~> degC].
918	123576	Sal, & ! The salinity after a timestep of mixing [S ~> ppt].
919	123576	u_test, v_test, & ! Temporary velocities [L T-1 ~> m s-1].
920	123576	T_test, S_test ! Temporary temperatures [C ~> degC] and salinities [S ~> ppt].
921
922		real, dimension(nzc+1) :: &
923	123576	N2, & ! The squared buoyancy frequency at an interface [T-2 ~> s-2].
924	123576	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	123576	dz_Int, & ! The vertical distance with the space surrounding an interface,
927		! as used in calculating kappa and TKE [Z ~> m]
928	123576	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	123576	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	123576	S2, & ! The squared shear at an interface [T-2 ~> s-2].
935	123576	a1, & ! a1 is the coupling between adjacent interfaces in the TKE,
936		! velocity, and density equations [H ~> m or kg m-2]
937	123576	c1, & ! c1 is used in the tridiagonal (and similar) solvers [nondim].
938	123576	k_src, & ! The shear-dependent source term in the kappa equation [T-1 ~> s-1]
939	123576	kappa_src, & ! The shear-dependent source term in the kappa equation [T-1 ~> s-1]
940	123576	kappa_out, & ! The kappa that results from the kappa equation [H Z T-1 ~> m2 s-1 or Pa s]
941	123576	kappa_mid, & ! The average of the initial and predictor estimates of kappa [H Z T-1 ~> m2 s-1 or Pa s]
942	123576	tke_pred, & ! The value of TKE from a predictor step [Z2 T-2 ~> m2 s-2].
943	123576	kappa_pred, & ! The value of kappa from a predictor step [H Z T-1 ~> m2 s-1 or Pa s]
944	123576	pressure, & ! The pressure at an interface [R L2 T-2 ~> Pa].
945	123576	T_int, & ! The temperature interpolated to an interface [C ~> degC].
946	123576	Sal_int, & ! The salinity interpolated to an interface [S ~> ppt].
947	123576	dbuoy_dT, & ! The partial derivative of buoyancy with changes in temperature [Z T-2 C-1 ~> m s-2 degC-1]
948	123576	dbuoy_dS, & ! The partial derivative of buoyancy with changes in salinity [Z T-2 S-1 ~> m s-2 ppt-1]
949	123576	dSpV_dT, & ! The partial derivative of specific volume with changes in temperature [R-1 C-1 ~> m3 kg-1 degC-1]
950	123576	dSpV_dS, & ! The partial derivative of specific volume with changes in salinity [R-1 S-1 ~> m3 kg-1 ppt-1]
951	123576	rho_int, & ! The in situ density interpolated to an interface [R ~> kg m-3]
952	123576	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	123576	K_Q, & ! Diffusivity divided by TKE [H T Z-1 ~> s or kg s m-3]
955	123576	K_Q_tmp, & ! A temporary copy of diffusivity divided by TKE [H T Z-1 ~> s or kg s m-3]
956	123576	local_src_avg, & ! The time-integral of the local source [nondim]
957	123576	tol_min, & ! Minimum tolerated ksrc for the corrector step [T-1 ~> s-1].
958	123576	tol_max, & ! Maximum tolerated ksrc for the corrector step [T-1 ~> s-1].
959	123576	tol_chg, & ! The tolerated kappa change integrated over a timestep [nondim].
960	123576	dist_from_top, & ! The distance from the top surface [Z ~> m].
961	123576	h_from_top, & ! The total thickness above an interface [H ~> m or kg m-2]
962	123576	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	61788	Ri_crit = CS%Rino_crit
1009	61788	gR0 = GV%H_to_RZ * GV%g_Earth
1010	61788	g_R0 = GV%g_Earth_Z_T2 / GV%Rho0
1011	61788	k0dt = dt*CS%kappa_0
1012
1013	61788	I_lz_rescale_sqr = 1.0 ; if (CS%lz_rescale > 0) I_lz_rescale_sqr = 1/(CS%lz_rescale*CS%lz_rescale)
1014
1015	61788	tol_dksrc = CS%kappa_src_max_chg
1016	61788	if (tol_dksrc == 10.0) then
1017		! This is equivalent to the expression below, but avoids changes at roundoff for the default value.
1018	61788	tol_dksrc_low = 0.95
1019		else
1020	0	tol_dksrc_low = (tol_dksrc - 0.5)/tol_dksrc
1021		endif
1022	61788	tol2 = 2.0*CS%kappa_tol_err
1023	61788	dt_refinements = 5 ! Selected so that 1/2^dt_refinements < 1-tol_dksrc_low
1024	61788	use_temperature = .false. ; if (associated(tv%T)) use_temperature = .true.
1025
1026
1027		! Set up Idz as the inverse of layer thicknesses.
1028	3345588	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	61788	I_dz_int(1) = 2.0 / dz_lay(1)
1032	61788	dist_from_top(1) = 0.0 ; h_from_top(1) = 0.0
1033	3283800	do K=2,nzc
1034	3222012	I_dz_int(K) = 2.0 / (dz_lay(k-1) + dz_lay(k))
1035	3222012	dist_from_top(K) = dist_from_top(K-1) + dz_lay(k-1)
1036	3283800	h_from_top(K) = h_from_top(K-1) + hlay(k-1)
1037		enddo
1038	61788	I_dz_int(nzc+1) = 2.0 / dz_lay(nzc)
1039
1040		! Find the inverse of the squared distances from the boundaries.
1041	61788	dist_from_bot = 0.0 ; h_from_bot = 0.0
1042	3283800	do K=nzc,2,-1
1043	3222012	dist_from_bot = dist_from_bot + dz_lay(k)
1044	3222012	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	3222012	((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	3283800	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	61788	if (nzc > 1) then
1055	61788	a1(2) = k0dt*I_dz_int(2)
1056	61788	b1 = 1.0 / (hlay(1) + a1(2))
1057	61788	u(1) = b1 * u0xdz(1) ; v(1) = b1 * v0xdz(1)
1058	61788	T(1) = b1 * T0xdz(1) ; Sal(1) = b1 * S0xdz(1)
1059	61788	c1(2) = a1(2) * b1 ; d1 = hlay(1) * b1 ! = 1 - c1
1060	3222012	do k=2,nzc-1
1061	3160224	bd1 = hlay(k) + d1*a1(k)
1062	3160224	a1(k+1) = k0dt*I_dz_int(k+1)
1063	3160224	b1 = 1.0 / (bd1 + a1(k+1))
1064	3160224	u(k) = b1 * (u0xdz(k) + a1(k)*u(k-1))
1065	3160224	v(k) = b1 * (v0xdz(k) + a1(k)*v(k-1))
1066	3160224	T(k) = b1 * (T0xdz(k) + a1(k)*T(k-1))
1067	3160224	Sal(k) = b1 * (S0xdz(k) + a1(k)*Sal(k-1))
1068	3222012	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	61788	b1 = 1.0 / ((hlay(nzc) + d1a1(nzc)) + k0dtI_dz_int(nzc+1))
1074	61788	u(nzc) = b1 * (u0xdz(nzc) + a1(nzc)*u(nzc-1))
1075	61788	v(nzc) = b1 * (v0xdz(nzc) + a1(nzc)*v(nzc-1))
1076		! For insulating boundary conditions
1077	61788	b1 = 1.0 / (hlay(nzc) + d1*a1(nzc))
1078	61788	T(nzc) = b1 * (T0xdz(nzc) + a1(nzc)*T(nzc-1))
1079	61788	Sal(nzc) = b1 * (S0xdz(nzc) + a1(nzc)*Sal(nzc-1))
1080	3283800	do k=nzc-1,1,-1
1081	3222012	u(k) = u(k) + c1(k+1)u(k+1) ; v(k) = v(k) + c1(k+1)v(k+1)
1082	3283800	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	0	b1 = 1.0 / (hlay(1) + k0dt*I_dz_int(2))
1087	0	u(1) = b1 * u0xdz(1) ; v(1) = b1 * v0xdz(1)
1088	0	b1 = 1.0 / hlay(1)
1089	0	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	61788	h_Int(1) = 0.0 ; h_Int(2) = hlay(1)
1099	61788	dz_Int(1) = 0.0 ; dz_Int(2) = dz_lay(1)
1100	3222012	do K=2,nzc-1
1101	3160224	Norm = 1.0 / (hlay(k)(hlay(k-1)+hlay(k+1)) + 2.0hlay(k-1)*hlay(k+1))
1102	3160224	wt_a = ((hlay(k)+hlay(k+1)) * hlay(k-1)) * Norm
1103	3160224	wt_b = ((hlay(k-1)+hlay(k)) * hlay(k+1)) * Norm
1104	3160224	h_Int(K) = h_Int(K) + hlay(k) * wt_a
1105	3160224	h_Int(K+1) = hlay(k) * wt_b
1106	3160224	dz_Int(K) = dz_Int(K) + dz_lay(k) * wt_a
1107	3222012	dz_Int(K+1) = dz_lay(k) * wt_b
1108		enddo
1109	61788	h_Int(nzc) = h_Int(nzc) + hlay(nzc) ; h_Int(nzc+1) = 0.0
1110	61788	dz_Int(nzc) = dz_Int(nzc) + dz_lay(nzc) ; dz_Int(nzc+1) = 0.0
1111
1112	61788	if (GV%Boussinesq) then
1113	3407376	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	61788	if (use_temperature) then
1129	61788	pressure(1) = surface_pres
1130	3283800	do K=2,nzc
1131	3222012	pressure(K) = pressure(K-1) + gR0*hlay(k-1)
1132	3222012	T_int(K) = 0.5*(T(k-1) + T(k))
1133	3283800	Sal_int(K) = 0.5*(Sal(k-1) + Sal(k))
1134		enddo
1135	61788	if (GV%Boussinesq .or. GV%semi_Boussinesq) then
1136		call calculate_density_derivs(T_int, Sal_int, pressure, dbuoy_dT, dbuoy_dS, &
1137	185364	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	61788	CS%vel_underflow, u, v, T, Sal, N2, S2, GV, US)
1170	3407376	do K=1,nzc+1
1171	3345588	N2_init(K) = N2(K)
1172	3407376	S2_init(K) = S2(K)
1173		enddo
1174
1175
1176		! ----------------------------------------------------
1177		! Iterate
1178		! ----------------------------------------------------
1179	61788	dt_rem = dt
1180	3407376	do K=1,nzc+1
1181	3345588	K_Q(K) = 0.0
1182	3345588	kappa_avg(K) = 0.0 ; tke_avg(K) = 0.0 ; N2_mean(K) = 0.0 ; S2_mean(K) = 0.0
1183	3345588	local_src_avg(K) = 0.0
1184		! Use the grid spacings to scale errors in the source.
1185	3345588	if ( h_Int(K) > 0.0 ) &
1186	3283800	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	229530	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	229530	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	229530	ks_kappa = GV%ke+1 ; ke_kappa = 0
1206	392813	do K=2,nzc ; if (kappa_out(K) > 0.0) then
1207	226250	ks_kappa = K ; exit
1208		endif ; enddo
1209	10228708	do k=nzc,ks_kappa,-1 ; if (kappa_out(K) > 0.0) then
1210	226250	ke_kappa = K ; exit
1211		endif ; enddo
1212	229530	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	229530	if ((ke_kappa < ks_kappa) .or. (itt==CS%max_KS_it)) then
1219	3280	dt_now = dt_rem
1220		else
1221		! Limit dt_now so that \|k_src(k)-kappa_src(k)\| < tol * local_src(k)
1222	226250	dt_test = dt_rem
1223	12317409	do K=2,nzc
1224	12091159	tol_max(K) = kappa_src(K) + tol_dksrc * local_src(K)
1225	12091159	tol_min(K) = kappa_src(K) - tol_dksrc_low * local_src(K)
1226	12317409	tol_chg(K) = tol2 * local_src_avg(K)
1227		enddo
1228
1229	932373	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	869644	T_test, S_test, N2, S2, GV, US, ks_int=ks_kappa, ke_int=ke_kappa)
1240	869644	valid_dt = .true.
1241	869644	Idtt = 1.0 / dt_test
1242	2602594	do K=max(ks_kappa-1,2),min(ke_kappa+1,nzc)
1243	2602594	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	2047528	((Ri_critS2(K) - N2(K)) / (Ri_critS2(K) + CS%FRi_curvature*N2(K)))
1246	2047528	if (CS%restrictive_tolerance_check) then
1247	2047528	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	704136	valid_dt = .false. ; exit
1250		endif
1251		else
1252	0	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	391545	if (0.0 < min(tol_min(K), kappa_src(K) - Idtt*tol_chg(K))) then
1259	1987	valid_dt = .false. ; k_src(K) = 0.0 ; exit
1260		endif
1261		endif
1262		enddo
1263
1264	869644	if (valid_dt) exit
1265	768852	dt_test = 0.5*dt_test
1266		enddo
1267	226250	if ((dt_test < dt_rem) .and. valid_dt) then
1268	116131	dt_inc = 0.5*dt_test
1269	696786	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	580655	N2, S2, GV, US, ks_int=ks_kappa, ke_int=ke_kappa)
1273	580655	valid_dt = .true.
1274	580655	Idtt = 1.0 / (dt_test+dt_inc)
1275	6484804	do K=max(ks_kappa-1,2),min(ke_kappa+1,nzc)
1276	6484804	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	4476258	((Ri_critS2(K) - N2(K)) / (Ri_critS2(K) + CS%FRi_curvature*N2(K)))
1279	4476258	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	202	valid_dt = .false. ; exit
1282		endif
1283		else
1284	1430519	if (0.0 < min(tol_min(K), kappa_src(K) - Idtt*tol_chg(K))) then
1285	2426	valid_dt = .false. ; k_src(K) = 0.0 ; exit
1286		endif
1287		endif
1288		enddo
1289
1290	580655	if (valid_dt) dt_test = dt_test + dt_inc
1291	696786	dt_inc = 0.5*dt_inc
1292		enddo
1293		else
1294	110119	dt_inc = 0.0
1295		endif
1296
1297	226250	dt_now = min(dt_test*(1.0+CS%kappa_tol_err)+dt_inc, dt_rem)
1298	12317409	do K=2,nzc
1299	12317409	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	229530	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	3280	dt_wt = dt_rem / dt ; dt_rem = 0.0
1310	173122	do K=1,nzc+1
1311	169842	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	173122	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	226250	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	12769909	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	226250	nzc, CS, GV, US, K_Q_tmp, tke_pred, kappa_pred)
1328		! call cpu_clock_end(id_clock_KQ)
1329
1330	226250	ks_kappa = GV%ke+1 ; ke_kappa = 0
1331	12769909	do K=1,nzc+1
1332	12543659	kappa_mid(K) = 0.5*(kappa_out(K) + kappa_pred(K))
1333	12543659	if ((kappa_mid(K) > 0.0) .and. (K<ks_kappa)) ks_kappa = K
1334	12769909	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	226250	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	226250	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	226250	dt_wt = dt_now / dt ; dt_rem = dt_rem - dt_now
1350	12769909	do K=1,nzc+1
1351	12543659	kappa_mid(K) = 0.5*(kappa_out(K) + kappa_pred(K))
1352	12543659	kappa_avg(K) = kappa_avg(K) + kappa_mid(K)*dt_wt
1353	12543659	tke_avg(K) = tke_avg(K) + dt_wt0.5(tke_pred(K) + tke(K))
1354	12543659	N2_mean(K) = N2_mean(K) + dt_wt*N2(K)
1355	12543659	S2_mean(K) = S2_mean(K) + dt_wt*S2(K)
1356	12769909	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	229530	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	167742	GV, US)
1368		! call cpu_clock_end(id_clock_project)
1369		endif
1370
1371	229530	if (dt_rem <= 0.0) exit
1372
1373		enddo ! end itt loop
1374
1375	61788	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	2132329	subroutine calculate_projected_state(kappa, u0, v0, T0, S0, dt, nz, dz, I_dz_int, dbuoy_dT, dbuoy_dS, &
1381	2132329	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	2132329	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	2132329	ks = 1 ; ke = nz
1421	2132329	if (present(ks_int)) ks = max(ks_int-1,1)
1422	2132329	if (present(ke_int)) ke = min(ke_int,nz)
1423
1424	2132329	if (ks > ke) return
1425
1426	2132329	if (dt > 0.0) then
1427	2070541	a_b = dt(kappa(ks+1)I_dz_int(ks+1))
1428	2070541	b1 = 1.0 / (dz(ks) + a_b)
1429	2070541	c1(ks+1) = a_b * b1 ; d1 = dz(ks) * b1 ! = 1 - c1
1430
1431	2070541	u(ks) = (b1 * dz(ks))u0(ks) ; v(ks) = (b1 dz(ks))*v0(ks)
1432	2070541	T(ks) = (b1 * dz(ks))T0(ks) ; Sal(ks) = (b1 dz(ks))*S0(ks)
1433	26577899	do K=ks+1,ke-1
1434	24507358	a_a = a_b
1435	24507358	a_b = dt(kappa(K+1)I_dz_int(K+1))
1436	24507358	bd1 = dz(k) + d1*a_a
1437	24507358	b1 = 1.0 / (bd1 + a_b)
1438	24507358	c1(K+1) = a_b * b1 ; d1 = bd1 * b1 ! d1 = 1 - c1
1439
1440	24507358	u(k) = b1 * (dz(k)u0(k) + a_au(k-1))
1441	24507358	v(k) = b1 * (dz(k)v0(k) + a_av(k-1))
1442	24507358	T(k) = b1 * (dz(k)T0(k) + a_aT(k-1))
1443	26577899	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	2070541	a_a = a_b
1450	2070541	b1 = 1.0 / (dz(ke) + d1*a_a)
1451	2070541	T(ke) = b1 * (dz(ke)T0(ke) + a_aT(ke-1))
1452	2070541	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	2070541	if (ke == nz) then
1458	173287	a_b = dt(kappa(nz+1)I_dz_int(nz+1))
1459	173287	b1nz_0 = 1.0 / ((dz(nz) + d1*a_a) + a_b)
1460		else
1461	1897254	b1nz_0 = b1
1462		endif
1463	2070541	u(ke) = b1nz_0 * (dz(ke)u0(ke) + a_au(ke-1))
1464	2070541	v(ke) = b1nz_0 * (dz(ke)v0(ke) + a_av(ke-1))
1465	2070541	if (abs(u(ke)) < vel_under) u(ke) = 0.0
1466	2070541	if (abs(v(ke)) < vel_under) v(ke) = 0.0
1467
1468	28648440	do k=ke-1,ks,-1
1469	26577899	u(k) = u(k) + c1(k+1)*u(k+1)
1470	26577899	v(k) = v(k) + c1(k+1)*v(k+1)
1471	26577899	if (abs(u(k)) < vel_under) u(k) = 0.0
1472	26577899	if (abs(v(k)) < vel_under) v(k) = 0.0
1473	26577899	T(k) = T(k) + c1(k+1)*T(k+1)
1474	28648440	Sal(k) = Sal(k) + c1(k+1)*Sal(k+1)
1475		enddo
1476		else ! dt <= 0.0
1477	3345588	do k=1,nz
1478	3283800	u(k) = u0(k) ; v(k) = v0(k) ; T(k) = T0(k) ; Sal(k) = S0(k)
1479	3283800	if (abs(u(k)) < vel_under) u(k) = 0.0
1480	3345588	if (abs(v(k)) < vel_under) v(k) = 0.0
1481		enddo
1482		endif
1483
1484		! Store the squared shear at interfaces
1485	2132329	S2(1) = 0.0 ; S2(nz+1) = 0.0
1486	2132329	if (ks > 1) &
1487	3	S2(ks) = (((u(ks)-u0(ks-1))2) + ((v(ks)-v0(ks-1))2)) * (US%L_to_ZI_dz_int(ks))*2
1488	31932240	do K=ks+1,ke
1489	31932240	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	2132329	if (ke<nz) &
1492	1897254	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	2132329	N2(1) = 0.0 ; N2(nz+1) = 0.0
1496	2132329	if (ks > 1) &
1497		N2(ks) = max(0.0, I_dz_int(ks) * &
1498	3	(dbuoy_dT(ks) * (T0(ks-1)-T(ks)) + dbuoy_dS(ks) * (S0(ks-1)-Sal(ks))))
1499	31932240	do K=ks+1,ke
1500		N2(K) = max(0.0, I_dz_int(K) * &
1501	31932240	(dbuoy_dT(K) * (T(k-1)-T(k)) + dbuoy_dS(K) * (Sal(k-1)-Sal(k))))
1502		enddo
1503	2132329	if (ke<nz) &
1504		N2(ke+1) = max(0.0, I_dz_int(ke+1) * &
1505	1897254	(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	1364060	subroutine find_kappa_tke(N2, S2, kappa_in, Idz, h_Int, dz_Int, dz_h_Int, I_L2_bdry, f2, &
1511	459060	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	678750	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	1364060	dQdz ! Half the partial derivative of TKE with depth [Z T-2 ~> m s-2].
1548		real, dimension(nz+1) :: &
1549	1364060	dK, & ! The change in kappa [H Z T-1 ~> m2 s-1 or Pa s].
1550	1364060	dQ, & ! The change in TKE [Z2 T-2 ~> m2 s-2].
1551	1364060	cQ, cK, & ! cQ and cK are the upward influences in the tridiagonal and
1552		! hexadiagonal solvers for the TKE and kappa equations [nondim].
1553	1364060	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	1364060	TKE_decay, & ! The local TKE decay rate [T-1 ~> s-1].
1555	1364060	k_src, & ! The source term in the kappa equation [T-1 ~> s-1].
1556	1364060	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	1364060	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	1364060	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	1364060	I_Ld2_debug, & ! A separate version of I_Ld2 for debugging [H-1 Z-1 ~> m-2 or m kg-1].
1628	1364060	kappa_prev, & ! The value of kappa at the start of the current iteration [H Z T-1 ~> m2 s-1 or Pa s]
1629	682030	TKE_prev ! The value of TKE at the start of the current iteration [Z2 T-2 ~> m2 s-2].
1630
1631	682030	c_N2 = CS%C_N2 ; c_S2 = CS%C_S2
1632	682030	q0 = CS%TKE_bg ; kappa0 = CS%kappa_0
1633	682030	TKE_min = max(CS%TKE_bg, 1.0E-20US%m_to_Z2US%T_to_s**2)
1634	682030	Ri_crit = CS%Rino_crit
1635	682030	Ilambda2 = 1.0 / CS%lambda**2
1636	682030	kappa_trunc = CS%kappa_trunc
1637	682030	do_Newton = .false. ; abort_Newton = .false.
1638	682030	tol_err = CS%kappa_tol_err
1639	682030	Newton_err = 0.2 ! This initial value may be automatically reduced later.
1640
1641	682030	ks_kappa = 2 ; ke_kappa = nz ; ks_kappa_prev = 2 ; ke_kappa_prev = nz
1642
1643	682030	ke_src = 0 ; ks_src = nz+1
1644	37118789	do K=2,nz
1645	37118789	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	5317291	((Ri_critS2(K) - N2(K)) / (Ri_critS2(K) + CS%FRi_curvature*N2(K)))
1651	5317291	ke_src = K
1652	5317291	if (ks_src > k) ks_src = K
1653		else
1654	31119468	k_src(K) = 0.0
1655		endif
1656		enddo
1657
1658		! If there is no source anywhere, return kappa(K) = 0.
1659	682030	if (ks_src > ke_src) then
1660	173122	do K=1,nz+1
1661	173122	kappa(K) = 0.0 ; K_Q(K) = 0.0 ; tke(K) = TKE_min
1662		enddo
1663	173122	if (present(kappa_src)) then ; do K=1,nz+1 ; kappa_src(K) = 0.0 ; enddo ; endif
1664	173122	if (present(local_src)) then ; do K=1,nz+1 ; local_src(K) = 0.0 ; enddo ; endif
1665	3280	return
1666		endif
1667
1668	38309727	do K=1,nz+1
1669	37630977	kappa(K) = kappa_in(K)
1670		! TKE_decay(K) = c_nsqrt(N2(K)) + c_ssqrt(S2(K)) ! The expression in JHL.
1671	37630977	TKE_decay(K) = sqrt(c_n2N2(K) + c_s2S2(K))
1672	38309727	if ((kappa(K) > 0.0) .and. (K_Q(K) > 0.0)) then
1673	5721200	TKE(K) = kappa(K) / K_Q(K) ! Perhaps take the max with TKE_min
1674		else
1675	31909777	TKE(K) = TKE_min
1676		endif
1677		enddo
1678		! Apply boundary conditions to kappa.
1679	678750	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	678750	eden2 = kappa0 * Idz(nz)
1685	678750	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	678750	e1(nz+1) = 0.0 ; ome = 1.0
1691		endif
1692	36952227	do k=nz,2,-1
1693	36273477	eden1 = h_Int(K)TKE_decay(K) + ome eden2
1694	36273477	eden2 = kappa0 * Idz(k-1)
1695	36273477	I_eden = 1.0 / (eden2 + eden1)
1696	36952227	e1(K) = eden2 * I_eden ; ome = eden1 * I_eden ! = 1-e1
1697		enddo
1698	678750	e1(1) = 0.0
1699
1700
1701		! Iterate here to convergence to within some tolerance of order tol_err.
1702	1184175	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	1184175	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	1183867	ke_tke = max(ke_kappa,ke_kappa_prev)+1
1717		! aQ is the coupling between adjacent interfaces [Z T-1 ~> m s-1].
1718	55529444	do k=1,min(ke_tke,nz)
1719	55529444	aQ(k) = (0.5(kappa(K)+kappa(K+1)) + kappa0) Idz(k)
1720		enddo
1721	1183867	dQ(1) = -TKE(1)
1722	1183867	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	1183867	tke(1) = q0 ; cQ(2) = 0.0 ; cQcomp = 1.0
1730		endif
1731	54098304	do K=2,ke_tke-1
1732	52914437	dQ(K) = -TKE(K)
1733	52914437	tke_src = dz_h_Int(K)(kappa(K) + kappa0)S2(K) + q0*TKE_decay(K)
1734	52914437	bQd1 = h_Int(K)(TKE_decay(K) + dz_h_Int(K)N2(K)K_Q(K)) + cQcompaQ(k-1)
1735	52914437	bQ = 1.0 / (bQd1 + aQ(k))
1736	52914437	tke(K) = bQ * (h_Int(K)tke_src + aQ(k-1)tke(K-1))
1737	54098304	cQ(K+1) = aQ(k) * bQ ; cQcomp = bQd1 * bQ ! = 1 - cQ
1738		enddo
1739	1183867	if ((ke_tke == nz+1) .and. .not.(tke_noflux_bottom_BC)) then
1740	936594	tke(nz+1) = TKE_min
1741	936594	dQ(nz+1) = 0.0
1742		else
1743	247273	k = ke_tke
1744	247273	tke_src = dz_h_Int(K)kappa0S2(K) + q0*TKE_decay(K) ! Uses that kappa(ke_tke) = 0
1745	247273	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	247273	bQ = 1.0 / ((h_Int(K)TKE_decay(K) + cQcompaQ(k-1)) + aQ(k))
1752	247273	cQ(K+1) = aQ(k) * bQ
1753		! Account for all changes deeper in the water column.
1754	247273	dQ(K) = -TKE(K)
1755		tke(K) = max((bQ * (h_Int(K)tke_src + aQ(k-1)tke(K-1)) + &
1756	247273	cQ(K+1)(tke(K+1) - e1(K+1)tke(K))) / (1.0 - cQ(K+1)*e1(K+1)), TKE_min)
1757	247273	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	1432005	do K=ke_tke+1,nz+1
1762	1432005	dQ(K) = e1(K)*dQ(K-1)
1763	1432005	tke(K) = max(tke(K) + dQ(K), TKE_min)
1764	1432005	if (abs(dQ(K)) < roundoff*tke(K)) exit
1765		enddo
1766	2305531	do K2=K+1,nz
1767	2244327	if (dQ(K2) == 0.0) exit
1768	2119462	dQ(K2) = 0.0
1769		enddo
1770		endif
1771		endif
1772	55282171	do K=ke_tke-1,1,-1
1773	54098304	tke(K) = max(tke(K) + cQ(K+1)*tke(K+1), TKE_min)
1774	55282171	dQ(K) = tke(K) + dQ(K)
1775		enddo
1776
1777		! ----------------------------------------------------
1778		! Calculate kappa, here defined at interfaces.
1779		! ----------------------------------------------------
1780
1781	1183867	ke_kappa_prev = ke_kappa ; ks_kappa_prev = ks_kappa
1782
1783	1183867	dK(1) = 0.0 ! kappa takes boundary values of 0.
1784	1183867	cK(2) = 0.0 ; cKcomp = 1.0
1785	37457344	if (itt == 1) then ; do K=2,nz
1786	36952227	I_Ld2(K) = dz_h_Int(K)(N2(K)Ilambda2 + f2) / tke(K) + I_L2_bdry(K)
1787		enddo ; endif
1788	12540233	do K=2,nz
1789	12535207	dK(K) = -kappa(K)
1790	12535207	if (itt>1) &
1791	5245872	I_Ld2(K) = dz_h_Int(K)(N2(K)Ilambda2 + f2) / tke(K) + I_L2_bdry(K)
1792	12535207	bKd1 = h_Int(K)I_Ld2(K) + cKcompIdz(k-1)
1793	12535207	bK = 1.0 / (bKd1 + Idz(k))
1794
1795	12535207	kappa(K) = bK * (Idz(k-1)kappa(K-1) + h_Int(K) K_src(K))
1796	12535207	cK(K+1) = Idz(k) * bK ; cKcomp = bKd1 * bK ! = 1 - cK(K+1)
1797
1798		! Neglect values that are smaller than kappa_trunc.
1799	12540233	if (kappa(K) < cKcomp*kappa_trunc) then
1800	1311814	kappa(K) = 0.0
1801	1311814	if (K > ke_src) then ; ke_kappa = k-1 ; K_Q(K) = 0.0 ; exit ; endif
1802	11223393	elseif (kappa(K) < 2.0cKcompkappa_trunc) then
1803	135373	kappa(K) = 2.0 * (kappa(K) - cKcomp*kappa_trunc)
1804		endif
1805		enddo
1806	1183867	K_Q(ke_kappa) = kappa(ke_kappa) / tke(ke_kappa)
1807	1183867	dK(ke_kappa) = dK(ke_kappa) + kappa(ke_kappa)
1808	30343845	do K=ke_kappa+2,ke_kappa_prev
1809	30343845	dK(K) = -kappa(K) ; kappa(K) = 0.0 ; K_Q(K) = 0.0
1810		enddo
1811	11356366	do K=ke_kappa-1,2,-1
1812	10172499	kappa(K) = kappa(K) + cK(K+1)*kappa(K+1)
1813		! Neglect values that are smaller than kappa_trunc.
1814	10172499	if (kappa(K) <= kappa_trunc) then
1815	128995	kappa(K) = 0.0
1816	128995	if (K < ks_src) then ; ks_kappa = k+1 ; K_Q(K) = 0.0 ; exit ; endif
1817	10043504	elseif (kappa(K) < 2.0*kappa_trunc) then
1818	16833	kappa(K) = 2.0 * (kappa(K) - kappa_trunc)
1819		endif
1820
1821	10172499	dK(K) = dK(K) + kappa(K)
1822	11356366	K_Q(K) = kappa(K) / tke(K)
1823		enddo
1824	1183867	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	308	ks_kappa_prev = ks_kappa ; ke_kappa_prev = ke_kappa ; ke_kappa = nz
1830	308	ks_kappa = 2
1831	308	dK(1) = 0.0 ; cK(2) = 0.0 ; cKcomp = 1.0 ; dKdQ(1) = 0.0
1832	308	aQ(1) = (0.5(kappa(1)+kappa(2))+kappa0) Idz(1)
1833	308	dQdz(1) = 0.5(TKE(1) - TKE(2))Idz(1)
1834	308	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	308	dQ(1) = 0.0 ; cQ(2) = 0.0 ; cQcomp = 1.0 ; dQmdK(2) = 0.0
1845		endif
1846	2942	do K=2,nz
1847	2636	I_Q = 1.0 / TKE(K)
1848	2636	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	2636	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	2636	decay_term_k = -Idz(k-1)dQmdK(K)dKdQ(K-1) + h_Int(K)*I_Ld2(K)
1856	2636	if (decay_term_k < 0.0) then ; abort_Newton = .true. ; exit ; endif
1857	2636	bK = 1.0 / (Idz(k) + Idz(k-1)*cKcomp + decay_term_k)
1858
1859	2636	cK(K+1) = bK * Idz(k)
1860	2636	cKcomp = bK * (Idz(k-1)*cKcomp + decay_term_k) ! = 1-cK(K+1)
1861	2636	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	2636	dz_Int(K)(N2(K)Ilambda2 + f2) * I_Q*2 kappa(K) )
1867		endif
1868	2636	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	2636	if (dK(K) <= cKcomp*(kappa_trunc - kappa(K))) then
1872	2	dK(K) = -cKcomp*kappa(K)
1873		! if (K > ke_src) then ; ke_kappa = k-1 ; K_Q(K) = 0.0 ; exit ; endif
1874	2634	elseif (dK(K) < cKcomp(2.0kappa_trunc - kappa(K))) then
1875	2	dK(K) = 2.0 * dK(K) - cKcomp(2.0kappa_trunc - kappa(K))
1876		endif
1877
1878		! Solve for dQ(K)...
1879	2636	aQ(k) = (0.5(kappa(K)+kappa(K+1))+kappa0) Idz(k)
1880	2636	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	2636	(aQ(k) * (TKE(K) - TKE(K+1)) - aQ(k-1) * (TKE(K-1) - TKE(K)))
1884	2636	v1 = aQ(k-1) + dQdz(k-1)*dKdQ(K-1)
1885		v2 = (v1dQmdK(K) + dQdz(k-1)cK(K)) + &
1886	2636	((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	2636	decay_term_Q = h_Int(K)TKE_decay(K) - dQdz(k-1)dKdQ(K-1)cQ(K) - v2dKdQ(K)
1891	2636	if (decay_term_Q < 0.0) then ; abort_Newton = .true. ; exit ; endif
1892	2636	bQ = 1.0 / (aQ(k) + (cQcomp*aQ(k-1) + decay_term_Q))
1893
1894	2636	cQ(K+1) = aQ(k) * bQ
1895	2636	cQcomp = (cQcompaQ(k-1) + decay_term_Q) bQ
1896	2636	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	2636	(v2 * dK(K) + tke_src)), cQcomp(-0.5TKE(K)))
1901
1902		! Check whether the next layer will be affected by any nonzero kappas.
1903	2636	if ((itt > 1) .and. (K > ke_src) .and. (dK(K) == 0.0) .and. &
1904	306	((kappa(K) + kappa(K+1)) == 0.0)) then
1905		! Could also do .and. (bQabs(tke_src) < roundoffTKE(K)) then
1906	2	ke_kappa = k-1 ; exit
1907		endif
1908		enddo
1909	308	if ((ke_kappa == nz) .and. (.not. abort_Newton)) then
1910	306	dK(nz+1) = 0.0 ; dKdQ(nz+1) = 0.0
1911	306	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	306	dQ(nz+1) = 0.0
1928		endif
1929	2	elseif (.not. abort_Newton) then
1930		! Alter the first-guess determination of dQ(K).
1931	2	dQ(ke_kappa+1) = dQ(ke_kappa+1) / (1.0 - cQ(ke_kappa+2)*e1(ke_kappa+2))
1932	2	TKE(ke_kappa+1) = max(TKE(ke_kappa+1) + dQ(ke_kappa+1), TKE_min)
1933	2	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	2	dK(K) = 0.0
1942		! Ensure that TKE+dQ will not drop below 0.5*TKE.
1943	2	dQ(K) = max(e1(K)dQ(K-1),-0.5TKE(K))
1944	2	TKE(K) = max(TKE(K) + dQ(K), TKE_min)
1945	2	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	308	if (.not. abort_Newton) then
1950	2942	do K=ke_kappa,2,-1
1951		! Ensure that TKE+dQ will not drop below 0.5*TKE.
1952	2634	dQ(K) = max(dQ(K) + (cQ(K+1)dQ(K+1) + dQmdK(K+1) dK(K+1)), -0.5*TKE(K))
1953	2634	TKE(K) = max(TKE(K) + dQ(K), TKE_min)
1954	2634	dK(K) = dK(K) + (cK(K+1)dK(K+1) + dKdQ(K) dQ(K))
1955		! Truncate away negligibly small values of kappa.
1956	2942	if (dK(K) <= kappa_trunc - kappa(K)) then
1957	2	dK(K) = -kappa(K)
1958	2	kappa(K) = 0.0
1959	2	if ((K < ks_src) .and. (K+1 > ks_kappa)) ks_kappa = K+1
1960	2632	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	2632	kappa(K) = kappa(K) + dK(K)
1966	2632	if (K<=ks_kappa) ks_kappa = 2
1967		endif
1968		enddo
1969	308	dQ(1) = max(dQ(1) + cQ(2)dQ(2) + dQmdK(2) dK(2), TKE_min - TKE(1))
1970	308	TKE(1) = max(TKE(1) + dQ(1), TKE_min)
1971	308	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	1184175	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	679058	Newton_test = Newton_err ; if (do_Newton) Newton_test = 2.0*Newton_err
2008	679058	was_Newton = do_Newton
2009	679058	within_tolerance = .true. ; do_Newton = .true.
2010	23871007	do K=min(ks_kappa,ks_kappa_prev),max(ke_kappa,ke_kappa_prev)
2011	23448939	kappa_mean = kappa0 + (kappa(K) - 0.5*dK(K))
2012	23448939	if (abs(dK(K)) > Newton_test * kappa_mean) then
2013	74562	if (do_Newton) abort_Newton = .true.
2014	74562	within_tolerance = .false. ; do_Newton = .false. ; exit
2015	23374377	elseif (abs(dK(K)) > tol_err * kappa_mean) then
2016	211879	within_tolerance = .false. ; if (.not.do_Newton) exit
2017		endif
2018	23793023	if (abs(dQ(K)) > Newton_test(tke(K) - 0.5dQ(K))) then
2019	13103837	if (do_Newton) abort_Newton = .true.
2020	13103837	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	505117	within_tolerance = .true.
2026	4283199	do K=min(ks_kappa,ks_kappa_prev),max(ke_kappa,ke_kappa_prev)
2027	4283199	if (abs(dK(K)) > tol_err * (kappa0 + (kappa(K) - 0.5*dK(K)))) then
2028	248127	within_tolerance = .false. ; exit
2029		endif
2030		enddo
2031		endif
2032
2033	1184175	if (abort_Newton) then
2034	677308	do_Newton = .false. ; abort_Newton = .false.
2035		! We went to Newton too quickly last time, so restrict the tolerance.
2036	677308	Newton_err = 0.5*Newton_err
2037	677308	ke_kappa_prev = nz
2038	36938390	do K=2,nz ; K_Q(K) = kappa(K) / max(TKE(K), TKE_min) ; enddo
2039		endif
2040
2041	1184175	if (within_tolerance) exit
2042
2043		enddo
2044
2045	678750	if (do_Newton) then ! K_Q needs to be calculated.
2046	2884	do K=1,ks_kappa-1 ; K_Q(K) = 0.0 ; enddo
2047	13835	do K=ks_kappa,ke_kappa ; K_Q(K) = kappa(K) / TKE(K) ; enddo
2048	2886	do K=ke_kappa+1,nz+1 ; K_Q(K) = 0.0 ; enddo
2049		endif
2050
2051	678750	if (present(local_src)) then
2052	226250	local_src(1) = 0.0 ; local_src(nz+1) = 0.0
2053	12317409	do K=2,nz
2054	12091159	diffusive_src = Idz(k-1)(kappa(K-1)-kappa(K)) + Idz(k)(kappa(K+1)-kappa(K))
2055	12091159	chg_by_k0 = kappa0 * ((Idz(k-1)+Idz(k)) / h_Int(K) + I_Ld2(K))
2056	12317409	if (diffusive_src <= 0.0) then
2057	11293448	local_src(K) = K_src(K) + chg_by_k0
2058		else
2059	797711	local_src(K) = (K_src(K) + chg_by_k0) + diffusive_src / h_Int(K)
2060		endif
2061		enddo
2062		endif
2063	678750	if (present(kappa_src)) then
2064	226250	kappa_src(1) = 0.0 ; kappa_src(nz+1) = 0.0
2065	12317409	do K=2,nz
2066	12317409	kappa_src(K) = K_src(K)
2067		enddo
2068		endif
2069
2070	911560	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	0	"kappa-shear equations.", default=.true., do_not_log=just_read)
2291	0	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	1	'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	1	'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	1	'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	1	'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	1	'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	1	'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	0	'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	0	'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	0	'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	0	'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	0	'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