src/parameterizations/vertical/MOM_tidal_mixing.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		!> Interface to vertical tidal mixing schemes including CVMix tidal mixing.
6		module MOM_tidal_mixing
7
8		use MOM_diag_mediator, only : diag_ctrl, time_type, register_diag_field
9		use MOM_diag_mediator, only : safe_alloc_ptr, post_data
10		use MOM_diagnose_Kdwork, only : vbf_CS
11		use MOM_debugging, only : hchksum
12		use MOM_error_handler, only : MOM_error, is_root_pe, FATAL, WARNING, NOTE
13		use MOM_file_parser, only : openParameterBlock, closeParameterBlock
14		use MOM_file_parser, only : get_param, log_param, log_version, param_file_type
15		use MOM_grid, only : ocean_grid_type
16		use MOM_io, only : slasher, MOM_read_data, field_size
17		use MOM_io, only : read_netCDF_data
18		use MOM_internal_tides, only : int_tide_CS, get_lowmode_loss
19		use MOM_remapping, only : remapping_CS, initialize_remapping, remapping_core_h
20		use MOM_string_functions, only : uppercase, lowercase
21		use MOM_unit_scaling, only : unit_scale_type
22		use MOM_variables, only : thermo_var_ptrs, p3d
23		use MOM_verticalGrid, only : verticalGrid_type
24		use CVMix_tidal, only : CVMix_init_tidal, CVMix_compute_Simmons_invariant
25		use CVMix_tidal, only : CVMix_coeffs_tidal, CVMix_tidal_params_type
26		use CVMix_tidal, only : CVMix_compute_Schmittner_invariant, CVMix_compute_SchmittnerCoeff
27		use CVMix_tidal, only : CVMix_coeffs_tidal_schmittner
28		use CVMix_kinds_and_types, only : CVMix_global_params_type
29		use CVMix_put_get, only : CVMix_put
30
31		implicit none ; private
32
33		#include <MOM_memory.h>
34
35		public tidal_mixing_init
36		public setup_tidal_diagnostics
37		public calculate_tidal_mixing
38		public post_tidal_diagnostics
39		public tidal_mixing_h_amp
40		public tidal_mixing_end
41
42		! A note on unit descriptions in comments: MOM6 uses units that can be rescaled for dimensional
43		! consistency testing. These are noted in comments with units like Z, H, L, and T, along with
44		! their mks counterparts with notation like "a velocity [Z T-1 ~> m s-1]". If the units
45		! vary with the Boussinesq approximation, the Boussinesq variant is given first.
46
47		!> Containers for tidal mixing diagnostics
48		type, public :: tidal_mixing_diags ; private
49		real, allocatable :: Kd_itidal(:,:,:) !< internal tide diffusivity at interfaces
50		!! [H Z T-1 ~> m2 s-1 or kg m-1 s-1]
51		real, allocatable :: Fl_itidal(:,:,:) !< vertical flux of tidal turbulent dissipation
52		!! [H Z2 T-3 ~> m3 s-3 or W m-2]
53		real, allocatable :: Kd_Niku(:,:,:) !< lee-wave diffusivity at interfaces [H Z T-1 ~> m2 s-1 or kg m-1 s-1]
54		real, allocatable :: Kd_Niku_work(:,:,:) !< layer integrated work by lee-wave driven mixing [R Z3 T-3 ~> W m-2]
55		real, allocatable :: Kd_Itidal_Work(:,:,:) !< layer integrated work by int tide driven mixing [R Z3 T-3 ~> W m-2]
56		real, allocatable :: Kd_Lowmode_Work(:,:,:) !< layer integrated work by low mode driven mixing [R Z3 T-3 ~> W m-2]
57		real, allocatable :: N2_int(:,:,:) !< Buoyancy frequency squared at interfaces [T-2 ~> s-2]
58		real, allocatable :: vert_dep_3d(:,:,:) !< The 3-d mixing energy deposition vertical fraction [nondim]?
59		real, allocatable :: Schmittner_coeff_3d(:,:,:) !< The coefficient in the Schmittner et al mixing scheme [nondim]
60		real, allocatable :: tidal_qe_md(:,:,:) !< Input tidal energy dissipated locally,
61		!! interpolated to model vertical coordinate [R Z3 T-3 ~> W m-2]
62		real, allocatable :: Kd_lowmode(:,:,:) !< internal tide diffusivity at interfaces
63		!! due to propagating low modes [H Z T-1 ~> m2 s-1 or kg m-1 s-1]
64		real, allocatable :: Fl_lowmode(:,:,:) !< vertical flux of tidal turbulent
65		!! dissipation due to propagating low modes [H Z2 T-3 ~> m3 s-3 or W m-2]
66		real, allocatable :: TKE_itidal_used(:,:) !< internal tide TKE input at ocean bottom [R Z3 T-3 ~> W m-2]
67		real, allocatable :: N2_bot(:,:) !< bottom squared buoyancy frequency [T-2 ~> s-2]
68		real, allocatable :: N2_meanz(:,:) !< vertically averaged buoyancy frequency [T-2 ~> s-2]
69		real, allocatable :: Polzin_decay_scale_scaled(:,:) !< Vertical scale of decay for tidal dissipation [Z ~> m]
70		real, allocatable :: Polzin_decay_scale(:,:) !< Vertical decay scale for tidal dissipation with Polzin [Z ~> m]
71		real, allocatable :: Simmons_coeff_2d(:,:) !< The Simmons et al mixing coefficient [nondim]
72		end type
73
74		!> Control structure with parameters for the tidal mixing module.
75		type, public :: tidal_mixing_cs ; private
76		logical :: debug = .true. !< If true, do more extensive debugging checks. This is hard-coded.
77
78		! Parameters
79		logical :: int_tide_dissipation = .false. !< Internal tide conversion (from barotropic)
80		!! with the schemes of St Laurent et al (2002) & Simmons et al (2004)
81
82		integer :: Int_tide_profile !< A coded integer indicating the vertical profile
83		!! for dissipation of the internal waves. Schemes that are
84		!! currently encoded are St Laurent et al (2002) and Polzin (2009).
85		logical :: Lee_wave_dissipation = .false. !< Enable lee-wave driven mixing, following
86		!! Nikurashin (2010), with a vertical energy
87		!! deposition profile specified by Lee_wave_profile to be
88		!! St Laurent et al (2002) or Simmons et al (2004) scheme
89
90		integer :: Lee_wave_profile !< A coded integer indicating the vertical profile
91		!! for dissipation of the lee waves. Schemes that are
92		!! currently encoded are St Laurent et al (2002) and
93		!! Polzin (2009).
94		real :: Int_tide_decay_scale !< decay scale for internal wave TKE [Z ~> m]
95
96		real :: Mu_itides !< efficiency for conversion of dissipation
97		!! to potential energy [nondim]
98
99		real :: Gamma_itides !< fraction of local dissipation [nondim]
100
101		real :: Gamma_lee !< fraction of local dissipation for lee waves
102		!! (Nikurashin's energy input) [nondim]
103		real :: Decay_scale_factor_lee !< Scaling factor for the decay scale of lee
104		!! wave energy dissipation [nondim]
105
106		real :: min_zbot_itides !< minimum depth for internal tide conversion [Z ~> m].
107		logical :: Lowmode_itidal_dissipation = .false. !< If true, consider mixing due to breaking low
108		!! modes that have been remotely generated using an internal tidal
109		!! dissipation scheme to specify the vertical profile of the energy
110		!! input to drive diapycnal mixing, along the lines of St. Laurent
111		!! et al. (2002) and Simmons et al. (2004).
112
113		real :: Nu_Polzin !< The non-dimensional constant used in Polzin form of
114		!! the vertical scale of decay of tidal dissipation [nondim]
115
116		real :: Nbotref_Polzin !< Reference value for the buoyancy frequency at the
117		!! ocean bottom used in Polzin formulation of the
118		!! vertical scale of decay of tidal dissipation [T-1 ~> s-1]
119		real :: Polzin_decay_scale_factor !< Scaling factor for the decay length scale
120		!! of the tidal dissipation profile in Polzin [nondim]
121		real :: Polzin_decay_scale_max_factor !< The decay length scale of tidal dissipation
122		!! profile in Polzin formulation should not exceed
123		!! Polzin_decay_scale_max_factor * depth of the ocean [nondim].
124		real :: Polzin_min_decay_scale !< minimum decay scale of the tidal dissipation
125		!! profile in Polzin formulation [Z ~> m]
126
127		real :: TKE_itide_max !< maximum internal tide conversion [R Z3 T-3 ~> W m-2]
128		!! available to mix above the BBL
129
130		real :: utide !< constant tidal amplitude [Z T-1 ~> m s-1] if READ_TIDEAMP is false.
131		real :: kappa_itides !< topographic wavenumber and non-dimensional scaling [Z-1 ~> m-1].
132		real :: kappa_h2_factor !< factor for the product of wavenumber * rms sgs height [nondim]
133		character(len=200) :: inputdir !< The directory in which to find input files
134
135		logical :: use_CVMix_tidal = .false. !< true if CVMix is to be used for determining
136		!! diffusivity due to tidal mixing
137
138		real :: min_thickness !< Minimum thickness allowed [Z ~> m]
139
140		! CVMix-specific parameters
141		integer :: CVMix_tidal_scheme = -1 !< 1 for Simmons, 2 for Schmittner
142		type(CVMix_tidal_params_type) :: CVMix_tidal_params !< A CVMix-specific type with parameters for tidal mixing
143		type(CVMix_global_params_type) :: CVMix_glb_params !< CVMix-specific for Prandtl number only
144		real :: tidal_max_coef !< CVMix-specific maximum allowable tidal
145		!! diffusivity. [Z2 T-1 ~> m2 s-1]
146		real :: tidal_diss_lim_tc !< CVMix-specific dissipation limit depth for
147		!! tidal-energy-constituent data [Z ~> m].
148		type(remapping_CS) :: remap_CS !< The control structure for remapping
149		integer :: remap_answer_date !< The vintage of the order of arithmetic and expressions to use
150		!! for remapping. Values below 20190101 recover the remapping
151		!! answers from 2018, while higher values use more robust
152		!! forms of the same remapping expressions.
153		integer :: tidal_answer_date !< The vintage of the order of arithmetic and expressions in the tidal
154		!! mixing calculations. Values below 20190101 recover the answers
155		!! from the end of 2018, while higher values use updated and more robust
156		!! forms of the same expressions.
157
158		type(int_tide_CS), pointer :: int_tide_CSp=> NULL() !< Control structure for a child module
159
160		! Data containers
161		real, allocatable :: TKE_Niku(:,:) !< Lee wave driven Turbulent Kinetic Energy input
162		!! [R Z3 T-3 ~> W m-2]
163		real, allocatable :: TKE_itidal(:,:) !< The internal Turbulent Kinetic Energy input divided by
164		!! the bottom stratification and in non-Boussinesq mode by
165		!! the near-bottom density [R Z4 H-1 T-2 ~> J m-2 or J m kg-1]
166		real, allocatable :: Nb(:,:) !< The near bottom buoyancy frequency [T-1 ~> s-1].
167		real, allocatable :: mask_itidal(:,:) !< A mask of where internal tide energy is input [nondim]
168		real, allocatable :: h2(:,:) !< Squared bottom depth variance [Z2 ~> m2].
169		real, allocatable :: tideamp(:,:) !< RMS tidal amplitude [Z T-1 ~> m s-1]
170		real, allocatable :: h_src(:) !< tidal constituent input layer thickness [m]
171		real, allocatable :: tidal_qe_2d(:,:) !< Tidal energy input times the local dissipation
172		!! fraction, q*E(x,y), with the CVMix implementation
173		!! of Jayne et al tidal mixing [R Z3 T-3 ~> W m-2].
174		!! TODO: make this E(x,y) only
175		real, allocatable :: tidal_qe_3d_in(:,:,:) !< q*E(x,y,z) with the Schmittner parameterization [R Z3 T-3 ~> W m-2]
176
177
178		! Diagnostics
179		type(diag_ctrl), pointer :: diag => NULL() !< structure to regulate diagnostic output timing
180		type(tidal_mixing_diags) :: dd !< Tidal mixing diagnostic arrays
181
182		!>@{ Diagnostic identifiers
183		integer :: id_TKE_itidal = -1
184		integer :: id_TKE_leewave = -1
185		integer :: id_Kd_itidal = -1
186		integer :: id_Kd_Niku = -1
187		integer :: id_Kd_lowmode = -1
188		integer :: id_Kd_Itidal_Work = -1
189		integer :: id_Kd_Niku_Work = -1
190		integer :: id_Kd_Lowmode_Work = -1
191		integer :: id_Nb = -1
192		integer :: id_N2_bot = -1
193		integer :: id_N2_meanz = -1
194		integer :: id_Fl_itidal = -1
195		integer :: id_Fl_lowmode = -1
196		integer :: id_Polzin_decay_scale = -1
197		integer :: id_Polzin_decay_scale_scaled = -1
198		integer :: id_N2_int = -1
199		integer :: id_Simmons_coeff = -1
200		integer :: id_Schmittner_coeff = -1
201		integer :: id_tidal_qe_md = -1
202		integer :: id_vert_dep = -1
203		!>@}
204
205		end type tidal_mixing_cs
206
207		!>@{ Coded parmameters for specifying mixing schemes
208		character*(20), parameter :: STLAURENT_PROFILE_STRING = "STLAURENT_02"
209		character*(20), parameter :: POLZIN_PROFILE_STRING = "POLZIN_09"
210		integer, parameter :: STLAURENT_02 = 1
211		integer, parameter :: POLZIN_09 = 2
212		character*(20), parameter :: SIMMONS_SCHEME_STRING = "SIMMONS"
213		character*(20), parameter :: SCHMITTNER_SCHEME_STRING = "SCHMITTNER"
214		integer, parameter :: SIMMONS = 1
215		integer, parameter :: SCHMITTNER = 2
216		!>@}
217
218		contains
219
220		!> Initializes internal tidal dissipation scheme for diapycnal mixing
221	1	logical function tidal_mixing_init(Time, G, GV, US, param_file, int_tide_CSp, diag, CS)
222		type(time_type), intent(in) :: Time !< The current time.
223		type(ocean_grid_type), intent(in) :: G !< Grid structure.
224		type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure.
225		type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
226		type(param_file_type), intent(in) :: param_file !< Run-time parameter file handle
227		type(int_tide_CS), pointer :: int_tide_CSp !< A pointer to the internal tides control structure
228		type(diag_ctrl), target, intent(inout) :: diag !< Diagnostics control structure.
229		type(tidal_mixing_cs), intent(inout) :: CS !< This module's control structure.
230
231		! Local variables
232		logical :: use_CVMix_tidal
233		logical :: int_tide_dissipation
234		logical :: read_tideamp
235		integer :: default_answer_date ! The default setting for the various ANSWER_DATE flags.
236		character(len=20) :: tmpstr, int_tide_profile_str
237		character(len=20) :: CVMix_tidal_scheme_str, tidal_energy_type
238		character(len=200) :: filename, h2_file, Niku_TKE_input_file ! Input file names
239		character(len=200) :: tideamp_file ! Input file names or paths
240		character(len=80) :: tideamp_var, rough_var, TKE_input_var ! Input file variable names
241		real :: hamp ! The magnitude of the sub-gridscale bottom depth variance [Z ~> m]
242		real :: utide ! The RMS tidal amplitude [Z T-1 ~> m s-1]
243		real :: max_frac_rough ! A limit on the depth variance as a fraction of the total depth [nondim]
244		real :: prandtl_tidal ! Prandtl number used by CVMix tidal mixing schemes to convert vertical
245		! diffusivities into viscosities [nondim]
246		real :: Niku_scale ! local variable for scaling the Nikurashin TKE flux data [nondim]
247		integer :: i, j, is, ie, js, je
248		integer :: isd, ied, jsd, jed
249
250		! This include declares and sets the variable "version".
251		# include "version_variable.h"
252		character(len=40) :: mdl = "MOM_tidal_mixing" !< This module's name.
253
254	1	is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec
255	1	isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed
256
257		! Read parameters
258		! NOTE: These are read twice because logfile output is streamed and we want
259		! to preserve the ordering of module header before parameters.
260		call get_param(param_file, mdl, "USE_CVMix_TIDAL", use_CVMix_tidal, &
261	1	default=.false., do_not_log=.true.)
262		call get_param(param_file, mdl, "INT_TIDE_DISSIPATION", int_tide_dissipation, &
263	1	default=use_CVMix_tidal, do_not_log=.true.)
264		call log_version(param_file, mdl, version, &
265		"Vertical Tidal Mixing Parameterization", &
266	1	all_default=.not.(use_CVMix_tidal .or. int_tide_dissipation))
267
268		call get_param(param_file, mdl, "USE_CVMix_TIDAL", use_CVMix_tidal, &
269		"If true, turns on tidal mixing via CVMix", &
270	1	default=.false.)
271		call get_param(param_file, mdl, "INT_TIDE_DISSIPATION", int_tide_dissipation, &
272		"If true, use an internal tidal dissipation scheme to "//&
273		"drive diapycnal mixing, along the lines of St. Laurent "//&
274	1	"et al. (2002) and Simmons et al. (2004).", default=use_CVMix_tidal)
275
276		! return if tidal mixing is inactive
277	1	tidal_mixing_init = int_tide_dissipation
278	1	if (.not. tidal_mixing_init) return
279
280	0	CS%debug = CS%debug.and.is_root_pe()
281	0	CS%diag => diag
282	0	if (associated(int_tide_CSp)) CS%int_tide_CSp => int_tide_CSp
283	0	CS%use_CVmix_tidal = use_CVmix_tidal
284	0	CS%int_tide_dissipation = int_tide_dissipation
285
286	0	call get_param(param_file, mdl, "INPUTDIR", CS%inputdir, default=".",do_not_log=.true.)
287	0	CS%inputdir = slasher(CS%inputdir)
288
289		call get_param(param_file, mdl, "DEFAULT_ANSWER_DATE", default_answer_date, &
290		"This sets the default value for the various _ANSWER_DATE parameters.", &
291	0	default=99991231)
292		call get_param(param_file, mdl, "TIDAL_MIXING_ANSWER_DATE", CS%tidal_answer_date, &
293		"The vintage of the order of arithmetic and expressions in the tidal mixing "//&
294		"calculations. Values below 20190101 recover the answers from the end of 2018, "//&
295		"while higher values use updated and more robust forms of the same expressions.", &
296	0	default=default_answer_date, do_not_log=.not.GV%Boussinesq)
297	0	if (.not.GV%Boussinesq) CS%tidal_answer_date = max(CS%tidal_answer_date, 20230701)
298
299		call get_param(param_file, mdl, "REMAPPING_ANSWER_DATE", CS%remap_answer_date, &
300		"The vintage of the expressions and order of arithmetic to use for remapping. "//&
301		"Values below 20190101 result in the use of older, less accurate expressions "//&
302		"that were in use at the end of 2018. Higher values result in the use of more "//&
303		"robust and accurate forms of mathematically equivalent expressions.", &
304	0	default=default_answer_date, do_not_log=.not.GV%Boussinesq)
305	0	if (.not.GV%Boussinesq) CS%remap_answer_date = max(CS%remap_answer_date, 20230701)
306
307	0	if (CS%int_tide_dissipation) then
308
309		! Read in CVMix tidal scheme if CVMix tidal mixing is on
310	0	if (CS%use_CVMix_tidal) then
311		call get_param(param_file, mdl, "CVMIX_TIDAL_SCHEME", CVMix_tidal_scheme_str, &
312		"CVMIX_TIDAL_SCHEME selects the CVMix tidal mixing "//&
313		"scheme with INT_TIDE_DISSIPATION. Valid values are:\n"//&
314		"\t SIMMONS - Use the Simmons et al (2004) tidal \n"//&
315		"\t mixing scheme.\n"//&
316		"\t SCHMITTNER - Use the Schmittner et al (2014) tidal \n"//&
317		"\t mixing scheme.", &
318	0	default=SIMMONS_SCHEME_STRING)
319	0	CVMix_tidal_scheme_str = uppercase(CVMix_tidal_scheme_str)
320
321	0	select case (CVMix_tidal_scheme_str)
322	0	case (SIMMONS_SCHEME_STRING) ; CS%CVMix_tidal_scheme = SIMMONS
323	0	case (SCHMITTNER_SCHEME_STRING) ; CS%CVMix_tidal_scheme = SCHMITTNER
324		case default
325		call MOM_error(FATAL, "tidal_mixing_init: Unrecognized setting "// &
326	0	"#define CVMIX_TIDAL_SCHEME "//trim(CVMix_tidal_scheme_str)//" found in input file.")
327		end select
328		endif ! CS%use_CVMix_tidal
329
330		! Read in vertical profile of tidal energy dissipation
331	0	if ( CS%CVMix_tidal_scheme == SCHMITTNER .or. .not. CS%use_CVMix_tidal) then
332		call get_param(param_file, mdl, "INT_TIDE_PROFILE", int_tide_profile_str, &
333		"INT_TIDE_PROFILE selects the vertical profile of energy "//&
334		"dissipation with INT_TIDE_DISSIPATION. Valid values are:\n"//&
335		"\t STLAURENT_02 - Use the St. Laurent et al exponential \n"//&
336		"\t decay profile.\n"//&
337		"\t POLZIN_09 - Use the Polzin WKB-stretched algebraic \n"//&
338		"\t decay profile.", &
339	0	default=STLAURENT_PROFILE_STRING)
340	0	int_tide_profile_str = uppercase(int_tide_profile_str)
341
342	0	select case (int_tide_profile_str)
343	0	case (STLAURENT_PROFILE_STRING) ; CS%int_tide_profile = STLAURENT_02
344	0	case (POLZIN_PROFILE_STRING) ; CS%int_tide_profile = POLZIN_09
345		case default
346		call MOM_error(FATAL, "tidal_mixing_init: Unrecognized setting "// &
347	0	"#define INT_TIDE_PROFILE "//trim(int_tide_profile_str)//" found in input file.")
348		end select
349		endif
350
351	0	elseif (CS%use_CVMix_tidal) then
352		call MOM_error(FATAL, "tidal_mixing_init: Cannot set INT_TIDE_DISSIPATION to False "// &
353	0	"when USE_CVMix_TIDAL is set to True.")
354		endif
355
356		call get_param(param_file, mdl, "LEE_WAVE_DISSIPATION", CS%Lee_wave_dissipation, &
357		"If true, use an lee wave driven dissipation scheme to "//&
358		"drive diapycnal mixing, along the lines of Nikurashin "//&
359		"(2010) and using the St. Laurent et al. (2002) "//&
360	0	"and Simmons et al. (2004) vertical profile", default=.false.)
361	0	if (CS%lee_wave_dissipation) then
362	0	if (CS%use_CVMix_tidal) then
363		call MOM_error(FATAL, "tidal_mixing_init: Lee wave driven dissipation scheme cannot "// &
364	0	"be used when CVMix tidal mixing scheme is active.")
365		endif
366		call get_param(param_file, mdl, "LEE_WAVE_PROFILE", tmpstr, &
367		"LEE_WAVE_PROFILE selects the vertical profile of energy "//&
368		"dissipation with LEE_WAVE_DISSIPATION. Valid values are:\n"//&
369		"\t STLAURENT_02 - Use the St. Laurent et al exponential \n"//&
370		"\t decay profile.\n"//&
371		"\t POLZIN_09 - Use the Polzin WKB-stretched algebraic \n"//&
372		"\t decay profile.", &
373	0	default=STLAURENT_PROFILE_STRING)
374	0	tmpstr = uppercase(tmpstr)
375	0	select case (tmpstr)
376	0	case (STLAURENT_PROFILE_STRING) ; CS%lee_wave_profile = STLAURENT_02
377	0	case (POLZIN_PROFILE_STRING) ; CS%lee_wave_profile = POLZIN_09
378		case default
379		call MOM_error(FATAL, "tidal_mixing_init: Unrecognized setting "// &
380	0	"#define LEE_WAVE_PROFILE "//trim(tmpstr)//" found in input file.")
381		end select
382		endif
383
384		call get_param(param_file, mdl, "INT_TIDE_LOWMODE_DISSIPATION", CS%Lowmode_itidal_dissipation, &
385		"If true, consider mixing due to breaking low modes that "//&
386		"have been remotely generated; as with itidal drag on the "//&
387		"barotropic tide, use an internal tidal dissipation scheme to "//&
388		"drive diapycnal mixing, along the lines of St. Laurent "//&
389	0	"et al. (2002) and Simmons et al. (2004).", default=.false.)
390
391	0	if ((CS%Int_tide_dissipation .and. (CS%int_tide_profile == POLZIN_09)) .or. &
392		(CS%lee_wave_dissipation .and. (CS%lee_wave_profile == POLZIN_09))) then
393	0	if (CS%use_CVMix_tidal) then
394		call MOM_error(FATAL, "tidal_mixing_init: Polzin scheme cannot "// &
395	0	"be used when CVMix tidal mixing scheme is active.")
396		endif
397		call get_param(param_file, mdl, "NU_POLZIN", CS%Nu_Polzin, &
398		"When the Polzin decay profile is used, this is a "//&
399		"non-dimensional constant in the expression for the "//&
400		"vertical scale of decay for the tidal energy dissipation.", &
401	0	units="nondim", default=0.0697)
402		call get_param(param_file, mdl, "NBOTREF_POLZIN", CS%Nbotref_Polzin, &
403		"When the Polzin decay profile is used, this is the "//&
404		"reference value of the buoyancy frequency at the ocean "//&
405		"bottom in the Polzin formulation for the vertical "//&
406		"scale of decay for the tidal energy dissipation.", &
407	0	units="s-1", default=9.61e-4, scale=US%T_to_s)
408		call get_param(param_file, mdl, "POLZIN_DECAY_SCALE_FACTOR", &
409		CS%Polzin_decay_scale_factor, &
410		"When the Polzin decay profile is used, this is a "//&
411		"scale factor for the vertical scale of decay of the tidal "//&
412	0	"energy dissipation.", default=1.0, units="nondim")
413		call get_param(param_file, mdl, "POLZIN_SCALE_MAX_FACTOR", &
414		CS%Polzin_decay_scale_max_factor, &
415		"When the Polzin decay profile is used, this is a factor "//&
416		"to limit the vertical scale of decay of the tidal "//&
417		"energy dissipation to POLZIN_DECAY_SCALE_MAX_FACTOR "//&
418	0	"times the depth of the ocean.", units="nondim", default=1.0)
419		call get_param(param_file, mdl, "POLZIN_MIN_DECAY_SCALE", CS%Polzin_min_decay_scale, &
420		"When the Polzin decay profile is used, this is the "//&
421		"minimum vertical decay scale for the vertical profile\n"//&
422		"of internal tide dissipation with the Polzin (2009) formulation", &
423	0	units="m", default=0.0, scale=US%m_to_Z)
424		endif
425
426	0	if (CS%Int_tide_dissipation .or. CS%Lee_wave_dissipation) then
427		call get_param(param_file, mdl, "INT_TIDE_DECAY_SCALE", CS%Int_tide_decay_scale, &
428		"The decay scale away from the bottom for tidal TKE with "//&
429		"the new coding when INT_TIDE_DISSIPATION is used.", &
430	0	units="m", default=500.0, scale=US%m_to_Z)
431		call get_param(param_file, mdl, "MU_ITIDES", CS%Mu_itides, &
432		"A dimensionless turbulent mixing efficiency used with "//&
433	0	"INT_TIDE_DISSIPATION, often 0.2.", units="nondim", default=0.2)
434		call get_param(param_file, mdl, "GAMMA_ITIDES", CS%Gamma_itides, &
435		"The fraction of the internal tidal energy that is "//&
436		"dissipated locally with INT_TIDE_DISSIPATION. "//&
437		"THIS NAME COULD BE BETTER.", &
438	0	units="nondim", default=0.3333)
439		call get_param(param_file, mdl, "MIN_ZBOT_ITIDES", CS%min_zbot_itides, &
440		"Turn off internal tidal dissipation when the total "//&
441	0	"ocean depth is less than this value.", units="m", default=0.0, scale=US%m_to_Z)
442		endif
443
444	0	if ( (CS%Int_tide_dissipation .or. CS%Lee_wave_dissipation) .and. &
445		.not. CS%use_CVMix_tidal) then
446
447	0	allocate(CS%Nb(isd:ied,jsd:jed), source=0.)
448	0	allocate(CS%h2(isd:ied,jsd:jed), source=0.)
449	0	allocate(CS%TKE_itidal(isd:ied,jsd:jed), source=0.)
450	0	allocate(CS%mask_itidal(isd:ied,jsd:jed), source=1.)
451
452		call get_param(param_file, mdl, "KAPPA_ITIDES", CS%kappa_itides, &
453		"A topographic wavenumber used with INT_TIDE_DISSIPATION. "//&
454		"The default is 2pi/10 km, as in St.Laurent et al. 2002.", &
455	0	units="m-1", default=8.e-4*atan(1.0), scale=US%Z_to_m)
456
457		call get_param(param_file, mdl, "UTIDE", CS%utide, &
458		"The constant tidal amplitude used with INT_TIDE_DISSIPATION.", &
459	0	units="m s-1", default=0.0, scale=US%m_to_Z*US%T_to_s)
460	0	allocate(CS%tideamp(is:ie,js:je), source=CS%utide)
461
462		call get_param(param_file, mdl, "KAPPA_H2_FACTOR", CS%kappa_h2_factor, &
463		"A scaling factor for the roughness amplitude with "//&
464	0	"INT_TIDE_DISSIPATION.", units="nondim", default=1.0)
465		call get_param(param_file, mdl, "TKE_ITIDE_MAX", CS%TKE_itide_max, &
466		"The maximum internal tide energy source available to mix "//&
467		"above the bottom boundary layer with INT_TIDE_DISSIPATION.", &
468	0	units="W m-2", default=1.0e3, scale=US%W_m2_to_RZ3_T3)
469
470		call get_param(param_file, mdl, "READ_TIDEAMP", read_tideamp, &
471		"If true, read a file (given by TIDEAMP_FILE) containing "//&
472	0	"the tidal amplitude with INT_TIDE_DISSIPATION.", default=.false.)
473	0	if (read_tideamp) then
474	0	if (CS%use_CVMix_tidal) then
475		call MOM_error(FATAL, "tidal_mixing_init: Tidal amplitude files are "// &
476	0	"not compatible with CVMix tidal mixing. ")
477		endif
478		call get_param(param_file, mdl, "TIDEAMP_FILE", tideamp_file, &
479		"The path to the file containing the spatially varying "//&
480	0	"tidal amplitudes with INT_TIDE_DISSIPATION.", default="tideamp.nc")
481	0	filename = trim(CS%inputdir) // trim(tideamp_file)
482	0	call log_param(param_file, mdl, "INPUTDIR/TIDEAMP_FILE", filename)
483		call get_param(param_file, mdl, "TIDEAMP_VARNAME", tideamp_var, &
484		"The name of the tidal amplitude variable in the input file.", &
485	0	default="tideamp")
486		! NOTE: There are certain cases where FMS is unable to read this file, so
487		! we use read_netCDF_data in place of MOM_read_data.
488		call read_netCDF_data(filename, tideamp_var, CS%tideamp, G%domain, &
489	0	rescale=US%m_to_Z*US%T_to_s)
490		endif
491
492		call get_param(param_file, mdl, "H2_FILE", h2_file, &
493		"The path to the file containing the sub-grid-scale "//&
494		"topographic roughness amplitude with INT_TIDE_DISSIPATION.", &
495	0	fail_if_missing=(.not.CS%use_CVMix_tidal))
496	0	filename = trim(CS%inputdir) // trim(h2_file)
497	0	call log_param(param_file, mdl, "INPUTDIR/H2_FILE", filename)
498		call get_param(param_file, mdl, "ROUGHNESS_VARNAME", rough_var, &
499		"The name in the input file of the squared sub-grid-scale "//&
500	0	"topographic roughness amplitude variable.", default="h2")
501		! NOTE: There are certain cases where FMS is unable to read this file, so
502		! we use read_netCDF_data in place of MOM_read_data.
503		call read_netCDF_data(filename, rough_var, CS%h2, G%domain, &
504	0	rescale=US%m_to_Z**2)
505
506		call get_param(param_file, mdl, "FRACTIONAL_ROUGHNESS_MAX", max_frac_rough, &
507		"The maximum topographic roughness amplitude as a fraction of the mean depth, "//&
508		"or a negative value for no limitations on roughness.", &
509	0	units="nondim", default=0.1)
510
511	0	do j=js,je ; do i=is,ie
512	0	if (max(G%meanSL(i,j) + G%bathyT(i,j), 0.0) < CS%min_zbot_itides) CS%mask_itidal(i,j) = 0.0
513	0	CS%tideamp(i,j) = CS%tideamp(i,j) * CS%mask_itidal(i,j) * G%mask2dT(i,j)
514
515		! Restrict rms topo to a fraction (often 10 percent) of the column depth.
516	0	if ((CS%tidal_answer_date < 20190101) .and. (max_frac_rough >= 0.0)) then
517	0	hamp = min(max_frac_rough * max(G%meanSL(i,j) + G%bathyT(i,j), 0.0), sqrt(CS%h2(i,j)))
518	0	CS%h2(i,j) = hamp*hamp
519		else
520	0	if (max_frac_rough >= 0.0) &
521	0	CS%h2(i,j) = min((max_frac_rough * max(G%meanSL(i,j) + G%bathyT(i,j), 0.0))**2, CS%h2(i,j))
522		endif
523
524	0	utide = CS%tideamp(i,j)
525		! Compute the fixed part of internal tidal forcing.
526		! The units here are [R Z4 H-1 T-2 ~> J m-2 or m3 s-2] here. (Note that J m-2 = kg s-2.)
527		CS%TKE_itidal(i,j) = 0.5 * CS%kappa_h2_factor * GV%H_to_RZ * &
528	0	CS%kappa_itides * CS%h2(i,j) * utide*utide
529		enddo ; enddo
530
531		endif
532
533	0	if (CS%Lee_wave_dissipation) then
534
535		call get_param(param_file, mdl, "NIKURASHIN_TKE_INPUT_FILE", Niku_TKE_input_file, &
536		"The path to the file containing the TKE input from lee "//&
537		"wave driven mixing. Used with LEE_WAVE_DISSIPATION.", &
538	0	fail_if_missing=.true.)
539		call get_param(param_file, mdl, "NIKURASHIN_SCALE", Niku_scale, &
540		"A non-dimensional factor by which to scale the lee-wave "//&
541		"driven TKE input. Used with LEE_WAVE_DISSIPATION.", &
542	0	units="nondim", default=1.0)
543
544	0	filename = trim(CS%inputdir) // trim(Niku_TKE_input_file)
545	0	call log_param(param_file, mdl, "INPUTDIR/NIKURASHIN_TKE_INPUT_FILE", filename)
546		call get_param(param_file, mdl, "TKE_INPUT_VAR", TKE_input_var, &
547		"The name in the input file of the turbulent kinetic energy input variable.", &
548	0	default="TKE_input")
549	0	allocate(CS%TKE_Niku(is:ie,js:je), source=0.)
550
551		call MOM_read_data(filename, TKE_input_var, CS%TKE_Niku, G%domain, timelevel=1, & ! ??? timelevel -aja
552	0	scale=Niku_scale*US%W_m2_to_RZ3_T3)
553
554		call get_param(param_file, mdl, "GAMMA_NIKURASHIN",CS%Gamma_lee, &
555		"The fraction of the lee wave energy that is dissipated "//&
556	0	"locally with LEE_WAVE_DISSIPATION.", units="nondim", default=0.3333)
557		call get_param(param_file, mdl, "DECAY_SCALE_FACTOR_LEE",CS%Decay_scale_factor_lee, &
558		"Scaling for the vertical decay scale of the local "//&
559	0	"dissipation of lee wave dissipation.", units="nondim", default=1.0)
560		else
561	0	CS%Decay_scale_factor_lee = -9.e99 ! This should never be used if CS%Lee_wave_dissipation = False
562		endif
563
564		! Configure CVMix
565	0	if (CS%use_CVMix_tidal) then
566
567		! Read in CVMix params
568		!call openParameterBlock(param_file,'CVMix_TIDAL')
569		call get_param(param_file, mdl, "TIDAL_MAX_COEF", CS%tidal_max_coef, &
570		"largest acceptable value for tidal diffusivity", &
571	0	units="m^2/s", default=50e-4, scale=US%m2_s_to_Z2_T) ! the default is 50e-4 in CVMix, 100e-4 in POP.
572		call get_param(param_file, mdl, "TIDAL_DISS_LIM_TC", CS%tidal_diss_lim_tc, &
573		"Min allowable depth for dissipation for tidal-energy-constituent data. "//&
574		"No dissipation contribution is applied above TIDAL_DISS_LIM_TC.", &
575	0	units="m", default=0.0, scale=US%m_to_Z)
576		call get_param(param_file, mdl, 'MIN_THICKNESS', CS%min_thickness, &
577	0	units="m", default=0.001, scale=US%m_to_Z, do_not_log=.True.)
578		call get_param(param_file, mdl, "PRANDTL_TIDAL", prandtl_tidal, &
579		"Prandtl number used by CVMix tidal mixing schemes "//&
580		"to convert vertical diffusivities into viscosities.", &
581	0	units="nondim", default=1.0, do_not_log=.true.)
582	0	call CVMix_put(CS%CVMix_glb_params, 'Prandtl', prandtl_tidal)
583
584		call get_param(param_file, mdl, "TIDAL_ENERGY_TYPE",tidal_energy_type, &
585		"The type of input tidal energy flux dataset. Valid values are"//&
586		"\t Jayne\n"//&
587		"\t ER03 \n",&
588	0	fail_if_missing=.true.)
589		! Check whether tidal energy input format and CVMix tidal mixing scheme are consistent
590	0	if ( .not. ( &
591		(uppercase(tidal_energy_type(1:4)) == 'JAYN' .and. CS%CVMix_tidal_scheme == SIMMONS).or. &
592		(uppercase(tidal_energy_type(1:4)) == 'ER03' .and. CS%CVMix_tidal_scheme == SCHMITTNER) ) )then
593		call MOM_error(FATAL, "tidal_mixing_init: Tidal energy file type ("//&
594		trim(tidal_energy_type)//") is incompatible with CVMix tidal "//&
595	0	" mixing scheme: "//trim(CVMix_tidal_scheme_str) )
596		endif
597	0	CVMix_tidal_scheme_str = lowercase(CVMix_tidal_scheme_str)
598
599		! Set up CVMix
600		call CVMix_init_tidal(CVmix_tidal_params_user = CS%CVMix_tidal_params, &
601		mix_scheme = CVMix_tidal_scheme_str, &
602		efficiency = CS%Mu_itides, &
603		vertical_decay_scale = CS%int_tide_decay_scale*US%Z_to_m, &
604		max_coefficient = CS%tidal_max_coef*US%Z2_T_to_m2_s, &
605		local_mixing_frac = CS%Gamma_itides, &
606	0	depth_cutoff = CS%min_zbot_itides*US%Z_to_m)
607
608	0	call read_tidal_energy(G, GV, US, tidal_energy_type, param_file, CS)
609
610		!call closeParameterBlock(param_file)
611
612		endif ! CVMix on
613
614		! Register Diagnostics fields
615
616	0	if (CS%Int_tide_dissipation .or. CS%Lee_wave_dissipation .or. &
617		CS%Lowmode_itidal_dissipation) then
618
619		CS%id_Kd_itidal = register_diag_field('ocean_model','Kd_itides',diag%axesTi,Time, &
620	0	'Internal Tide Driven Diffusivity', 'm2 s-1', conversion=GV%HZ_T_to_m2_s)
621
622	0	if (CS%use_CVMix_tidal) then
623		CS%id_N2_int = register_diag_field('ocean_model','N2_int',diag%axesTi,Time, &
624	0	'Bouyancy frequency squared, at interfaces', 's-2', conversion=US%s_to_T**2)
625		!> TODO: add units
626	0	if (CS%CVMix_tidal_scheme .eq. SIMMONS) then
627		CS%id_Simmons_coeff = register_diag_field('ocean_model','Simmons_coeff',diag%axesT1,Time, &
628	0	'time-invariant portion of the tidal mixing coefficient using the Simmons', '')
629	0	else if (CS%CVMix_tidal_scheme .eq. SCHMITTNER) then
630		CS%id_Schmittner_coeff = register_diag_field('ocean_model','Schmittner_coeff',diag%axesTL,Time, &
631	0	'time-invariant portion of the tidal mixing coefficient using the Schmittner', '')
632		CS%id_tidal_qe_md = register_diag_field('ocean_model','tidal_qe_md',diag%axesTL,Time, &
633		'input tidal energy dissipated locally interpolated to model vertical coordinates', &
634	0	'W m-2', conversion=US%RZ3_T3_to_W_m2)
635		endif
636		CS%id_vert_dep = register_diag_field('ocean_model','vert_dep',diag%axesTi,Time, &
637	0	'vertical deposition function needed for Simmons et al tidal mixing', '')
638		else
639		CS%id_TKE_itidal = register_diag_field('ocean_model','TKE_itidal',diag%axesT1,Time, &
640		'Internal Tide Driven Turbulent Kinetic Energy', &
641	0	'W m-2', conversion=US%RZ3_T3_to_W_m2)
642		CS%id_Nb = register_diag_field('ocean_model','Nb',diag%axesT1,Time, &
643	0	'Bottom Buoyancy Frequency', 's-1', conversion=US%s_to_T)
644
645		CS%id_Kd_lowmode = register_diag_field('ocean_model','Kd_lowmode',diag%axesTi,Time, &
646		'Internal Tide Driven Diffusivity (from propagating low modes)', &
647	0	'm2 s-1', conversion=GV%HZ_T_to_m2_s)
648
649		CS%id_Fl_itidal = register_diag_field('ocean_model','Fl_itides',diag%axesTi,Time, &
650		'Vertical flux of tidal turbulent dissipation', &
651	0	'm3 s-3', conversion=(GV%H_to_mUS%Z_to_m2US%s_to_T**3))
652
653		CS%id_Fl_lowmode = register_diag_field('ocean_model','Fl_lowmode',diag%axesTi,Time, &
654		'Vertical flux of tidal turbulent dissipation (from propagating low modes)', &
655	0	'm3 s-3', conversion=(GV%H_to_mUS%Z_to_m2US%s_to_T**3))
656
657		CS%id_Polzin_decay_scale = register_diag_field('ocean_model','Polzin_decay_scale', diag%axesT1, Time, &
658		'Vertical decay scale for the tidal turbulent dissipation with Polzin scheme', &
659	0	units='m', conversion=US%Z_to_m)
660
661		CS%id_Polzin_decay_scale_scaled = register_diag_field('ocean_model', &
662		'Polzin_decay_scale_scaled', diag%axesT1, Time, &
663		'Vertical decay scale for the tidal turbulent dissipation with Polzin scheme, '// &
664	0	'scaled by N2_bot/N2_meanz', units='m', conversion=US%Z_to_m)
665
666		CS%id_N2_bot = register_diag_field('ocean_model','N2_b',diag%axesT1,Time, &
667	0	'Bottom Buoyancy frequency squared', 's-2', conversion=US%s_to_T**2)
668
669		CS%id_N2_meanz = register_diag_field('ocean_model','N2_meanz', diag%axesT1, Time, &
670	0	'Buoyancy frequency squared averaged over the water column', 's-2', conversion=US%s_to_T**2)
671
672		CS%id_Kd_Itidal_Work = register_diag_field('ocean_model','Kd_Itidal_Work',diag%axesTL,Time, &
673		'Work done by Internal Tide Diapycnal Mixing', &
674	0	'W m-2', conversion=US%RZ3_T3_to_W_m2)
675
676		CS%id_Kd_Niku_Work = register_diag_field('ocean_model','Kd_Nikurashin_Work',diag%axesTL,Time, &
677		'Work done by Nikurashin Lee Wave Drag Scheme', &
678	0	'W m-2', conversion=US%RZ3_T3_to_W_m2)
679
680		CS%id_Kd_Lowmode_Work = register_diag_field('ocean_model','Kd_Lowmode_Work',diag%axesTL,Time, &
681		'Work done by Internal Tide Diapycnal Mixing (low modes)', &
682	0	'W m-2', conversion=US%RZ3_T3_to_W_m2)
683
684	0	if (CS%Lee_wave_dissipation) then
685		CS%id_TKE_leewave = register_diag_field('ocean_model','TKE_leewave',diag%axesT1,Time, &
686		'Lee wave Driven Turbulent Kinetic Energy', &
687	0	'W m-2', conversion=US%RZ3_T3_to_W_m2)
688		CS%id_Kd_Niku = register_diag_field('ocean_model','Kd_Nikurashin',diag%axesTi,Time, &
689	0	'Lee Wave Driven Diffusivity', 'm2 s-1', conversion=GV%HZ_T_to_m2_s)
690		endif
691		endif ! S%use_CVMix_tidal
692		endif
693
694	0	end function tidal_mixing_init
695
696
697		!> Depending on whether or not CVMix is active, calls the associated subroutine to compute internal
698		!! tidal dissipation and to add the effect of internal-tide-driven mixing to the layer or interface
699		!! diffusivities.
700	0	subroutine calculate_tidal_mixing(dz, j, N2_bot, Rho_bot, N2_lay, N2_int, TKE_to_Kd, max_TKE, &
701	0	G, GV, US, CS, Kd_max, Kv, Kd_lay, Kd_int, VBF)
702		type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure
703		type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure
704		type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
705		real, dimension(SZI_(G),SZK_(GV)), intent(in) :: dz !< The vertical distance across layers [Z ~> m]
706		integer, intent(in) :: j !< The j-index to work on
707		real, dimension(SZI_(G)), intent(in) :: N2_bot !< The near-bottom squared buoyancy
708		!! frequency [T-2 ~> s-2].
709		real, dimension(SZI_(G)), intent(in) :: Rho_bot !< The near-bottom in situ density [R ~> kg m-3]
710		real, dimension(SZI_(G),SZK_(GV)), intent(in) :: N2_lay !< The squared buoyancy frequency of the
711		!! layers [T-2 ~> s-2].
712		real, dimension(SZI_(G),SZK_(GV)+1), intent(in) :: N2_int !< The squared buoyancy frequency at the
713		!! interfaces [T-2 ~> s-2].
714		real, dimension(SZI_(G),SZK_(GV)), intent(in) :: TKE_to_Kd !< The conversion rate between the TKE
715		!! dissipated within a layer and the
716		!! diapycnal diffusivity within that layer,
717		!! usually (~Rho_0 / (G_Earth * dRho_lay))
718		!! [T2 Z-1 ~> s2 m-1]
719		real, dimension(SZI_(G),SZK_(GV)), intent(in) :: max_TKE !< The energy required for a layer to
720		!! entrain to its maximum realizable
721		!! thickness [H Z2 T-3 ~> m3 s-3 or W m-2]
722		type(tidal_mixing_cs), intent(inout) :: CS !< The control structure for this module
723		real, intent(in) :: Kd_max !< The maximum increment for diapycnal
724		!! diffusivity due to TKE-based processes,
725		!! [H Z T-1 ~> m2 s-1 or kg m-1 s-1]
726		!! Set this to a negative value to have no limit.
727		real, dimension(:,:,:), pointer :: Kv !< The "slow" vertical viscosity at each interface
728		!! (not layer!) [H Z T-1 ~> m2 s-1 or Pa s]
729		real, dimension(SZI_(G),SZK_(GV)), &
730		optional, intent(inout) :: Kd_lay !< The diapycnal diffusivity in layers
731		!! [H Z T-1 ~> m2 s-1 or kg m-1 s-1]
732		real, dimension(SZI_(G),SZK_(GV)+1), &
733		optional, intent(inout) :: Kd_int !< The diapycnal diffusivity at interfaces
734		!! [H Z T-1 ~> m2 s-1 or kg m-1 s-1]
735		type(vbf_CS), pointer :: VBF !< A diagnostic structure for vertical buoyancy fluxes
736
737	0	if (CS%Int_tide_dissipation .or. CS%Lee_wave_dissipation .or. CS%Lowmode_itidal_dissipation) then
738	0	if (CS%use_CVMix_tidal) then
739	0	call calculate_CVMix_tidal(dz, j, N2_int, G, GV, US, CS, Kv, Kd_lay, Kd_int)
740		else
741		call add_int_tide_diffusivity(dz, j, N2_bot, Rho_bot, N2_lay, TKE_to_Kd, max_TKE, &
742	0	G, GV, US, CS, Kd_max, Kd_lay, Kd_int, VBF)
743		endif
744		endif
745	0	end subroutine calculate_tidal_mixing
746
747
748		!> Calls the CVMix routines to compute tidal dissipation and to add the effect of internal-tide-driven
749		!! mixing to the interface diffusivities.
750	0	subroutine calculate_CVMix_tidal(dz, j, N2_int, G, GV, US, CS, Kv, Kd_lay, Kd_int)
751		type(ocean_grid_type), intent(in) :: G !< Grid structure.
752		type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure
753		type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
754		type(tidal_mixing_cs), intent(inout) :: CS !< This module's control structure.
755		real, dimension(SZI_(G),SZK_(GV)), intent(in) :: dz !< The vertical distance across layers [Z ~> m]
756		integer, intent(in) :: j !< The j-index to work on
757		real, dimension(SZI_(G),SZK_(GV)+1), intent(in) :: N2_int !< The squared buoyancy
758		!! frequency at the interfaces [T-2 ~> s-2].
759		real, dimension(:,:,:), pointer :: Kv !< The "slow" vertical viscosity at each interface
760		!! (not layer!) [H Z T-1 ~> m2 s-1 or Pa s]
761		real, dimension(SZI_(G),SZK_(GV)), &
762		optional, intent(inout) :: Kd_lay!< The diapycnal diffusivity in the layers
763		!! [H Z T-1 ~> m2 s-1 or kg m-1 s-1].
764		real, dimension(SZI_(G),SZK_(GV)+1), &
765		optional, intent(inout) :: Kd_int!< The diapycnal diffusivity at interfaces
766		!! [H Z T-1 ~> m2 s-1 or kg m-1 s-1].
767		! Local variables
768	0	real, dimension(SZK_(GV)+1) :: Kd_tidal ! tidal diffusivity [m2 s-1]
769	0	real, dimension(SZK_(GV)+1) :: Kv_tidal ! tidal viscosity [m2 s-1]
770	0	real, dimension(SZK_(GV)+1) :: vert_dep ! vertical deposition [nondim]
771	0	real, dimension(SZK_(GV)+1) :: iFaceHeight ! Height of interfaces [m]
772	0	real, dimension(SZK_(GV)+1) :: SchmittnerSocn ! A larger value of the Schmittner coefficint to
773		! use in the Southern Ocean [nondim]. If this is smaller
774		! than Schmittner_coeff, that standard value is used.
775	0	real, dimension(SZK_(GV)) :: cellHeight ! Height of cell centers [m]
776	0	real, dimension(SZK_(GV)) :: tidal_qe_md ! Tidal dissipation energy interpolated from 3d input
777		! to model coordinates [R Z3 T-3 ~> W m-2]
778	0	real, dimension(SZK_(GV)+1) :: N2_int_i ! De-scaled interface buoyancy frequency [s-2]
779	0	real, dimension(SZK_(GV)) :: Schmittner_coeff ! A coefficient in the Schmittner et al (2014) mixing
780		! parameterization [nondim]
781	0	real, dimension(SZK_(GV)) :: h_m ! Cell thickness [m]
782	0	real, allocatable, dimension(:,:) :: exp_hab_zetar ! A badly documented array that appears to be
783		! related to the distribution of tidal mixing energy, with unusual array
784		! extents that are not explained, that is set and used by the CVMix
785		! tidal mixing schemes, perhaps in [m3 kg-1]?
786		real :: dh, hcorr ! Limited thicknesses and a cumulative correction [Z ~> m]
787		real :: Simmons_coeff ! A coefficient in the Simmons et al (2004) mixing parameterization [nondim]
788
789		integer :: i, k, is, ie
790		real, parameter :: rho_fw = 1000.0 ! fresh water density [kg m-3]
791		! TODO: when coupled, get this from CESM (SHR_CONST_RHOFW)
792
793	0	is = G%isc ; ie = G%iec
794
795	0	select case (CS%CVMix_tidal_scheme)
796		case (SIMMONS)
797	0	do i=is,ie
798
799	0	if (G%mask2dT(i,j)<1) cycle
800
801	0	iFaceHeight = 0.0 ! BBL is all relative to the surface
802	0	hcorr = 0.0
803		! Compute cell center depth and cell bottom in meters (negative values in the ocean)
804	0	do k=1,GV%ke
805	0	dh = dz(i,k) ! Nominal thickness to use for increment, in the units of heights
806	0	dh = dh + hcorr ! Take away the accumulated error (could temporarily make dh<0)
807	0	hcorr = min( dh - CS%min_thickness, 0. ) ! If inflating then hcorr<0
808	0	dh = max(dh, CS%min_thickness) ! Limited increment dh>=min_thickness
809	0	cellHeight(k) = iFaceHeight(k) - 0.5 * US%Z_to_m*dh
810	0	iFaceHeight(k+1) = iFaceHeight(k) - US%Z_to_m*dh
811		enddo
812
813		call CVMix_compute_Simmons_invariant( nlev = GV%ke, &
814		energy_flux = US%RZ3_T3_to_W_m2*CS%tidal_qe_2d(i,j), &
815		rho = rho_fw, &
816		SimmonsCoeff = Simmons_coeff, &
817		VertDep = vert_dep, &
818		zw = iFaceHeight, &
819		zt = cellHeight, &
820	0	CVMix_tidal_params_user = CS%CVMix_tidal_params)
821
822		! Since we pass tidal_qe_2d=(CS%Gamma_itides)*tidal_energy_flux_2d, and not tidal_energy_flux_2d in
823		! above subroutine call, we divide Simmons_coeff by CS%Gamma_itides as a corrective step:
824		! TODO: (CS%Gamma_itides)*tidal_energy_flux_2d is unnecessary, directly use tidal_energy_flux_2d
825	0	Simmons_coeff = Simmons_coeff / CS%Gamma_itides
826
827
828		! XXX: Temporary de-scaling of N2_int(i,:) into a temporary variable
829	0	do K=1,GV%ke+1
830	0	N2_int_i(K) = US%s_to_T*2 N2_int(i,K)
831		enddo
832
833		call CVMix_coeffs_tidal( Mdiff_out = Kv_tidal, &
834		Tdiff_out = Kd_tidal, &
835		Nsqr = N2_int_i, &
836		OceanDepth = -iFaceHeight(GV%ke+1),&
837		SimmonsCoeff = Simmons_coeff, &
838		vert_dep = vert_dep, &
839		nlev = GV%ke, &
840		max_nlev = GV%ke, &
841		CVMix_params = CS%CVMix_glb_params, &
842	0	CVMix_tidal_params_user = CS%CVMix_tidal_params)
843
844		! Update diffusivity
845	0	if (present(Kd_lay)) then
846	0	do k=1,GV%ke
847	0	Kd_lay(i,k) = Kd_lay(i,k) + 0.5 * GV%m2_s_to_HZ_T * (Kd_tidal(k) + Kd_tidal(k+1))
848		enddo
849		endif
850	0	if (present(Kd_int)) then
851	0	do K=1,GV%ke+1
852	0	Kd_int(i,K) = Kd_int(i,K) + GV%m2_s_to_HZ_T * Kd_tidal(K)
853		enddo
854		endif
855		! Update viscosity with the proper unit conversion.
856	0	if (associated(Kv)) then
857	0	do K=1,GV%ke+1
858	0	Kv(i,j,K) = Kv(i,j,K) + GV%m2_s_to_HZ_T * Kv_tidal(K) ! Rescale from m2 s-1 to H Z T-1.
859		enddo
860		endif
861
862		! diagnostics
863	0	if (allocated(CS%dd%Kd_itidal)) then
864	0	CS%dd%Kd_itidal(i,j,:) = GV%m2_s_to_HZ_T * Kd_tidal(:)
865		endif
866	0	if (allocated(CS%dd%N2_int)) then
867	0	CS%dd%N2_int(i,j,:) = N2_int(i,:)
868		endif
869	0	if (allocated(CS%dd%Simmons_coeff_2d)) then
870	0	CS%dd%Simmons_coeff_2d(i,j) = Simmons_coeff
871		endif
872	0	if (allocated(CS%dd%vert_dep_3d)) then
873	0	CS%dd%vert_dep_3d(i,j,:) = vert_dep(:)
874		endif
875
876		enddo ! i=is,ie
877
878		case (SCHMITTNER)
879
880		! TODO: correct exp_hab_zetar shapes in CVMix_compute_Schmittner_invariant
881		! and CVMix_compute_SchmittnerCoeff low subroutines
882
883	0	allocate(exp_hab_zetar(GV%ke+1,GV%ke+1))
884
885	0	do i=is,ie
886
887	0	if (G%mask2dT(i,j)<1) cycle
888
889	0	iFaceHeight(:) = 0.0 ! BBL is all relative to the surface
890	0	hcorr = 0.0
891		! Compute heights at cell center and interfaces, and rescale layer thicknesses
892	0	do k=1,GV%ke
893	0	h_m(k) = dz(i,k)*US%Z_to_m ! Rescale thicknesses to m for use by CVmix.
894	0	dh = dz(i,k) ! Nominal thickness to use for increment, in the units of heights
895	0	dh = dh + hcorr ! Take away the accumulated error (could temporarily make dh<0)
896	0	hcorr = min( dh - CS%min_thickness, 0. ) ! If inflating then hcorr<0
897	0	dh = max(dh, CS%min_thickness) ! Limited increment dh>=min_thickness
898	0	cellHeight(k) = iFaceHeight(k) - 0.5 * US%Z_to_m*dh
899	0	iFaceHeight(k+1) = iFaceHeight(k) - US%Z_to_m*dh
900		enddo
901
902	0	SchmittnerSocn = 0.0 ! TODO: compute this
903
904		! form the time-invariant part of Schmittner coefficient term
905		call CVMix_compute_Schmittner_invariant(nlev = GV%ke, &
906		VertDep = vert_dep, &
907		efficiency = CS%Mu_itides, &
908		rho = rho_fw, &
909		exp_hab_zetar = exp_hab_zetar, &
910		zw = iFaceHeight, &
911	0	CVmix_tidal_params_user = CS%CVMix_tidal_params)
912		!TODO: in above call, there is no need to pass efficiency, since it gets
913		! passed via CVMix_init_tidal and stored in CVMix_tidal_params. Change
914		! CVMix API to prevent this redundancy.
915
916		! remap from input z coordinate to model coordinate:
917	0	tidal_qe_md(:) = 0.0
918		call remapping_core_h(CS%remap_cs, size(CS%h_src), CS%h_src, CS%tidal_qe_3d_in(i,j,:), &
919	0	GV%ke, h_m, tidal_qe_md)
920
921		! form the Schmittner coefficient that is based on 3D q*E, which is formed from
922		! summing q_i*TidalConstituent_i over the number of constituents.
923		call CVMix_compute_SchmittnerCoeff( nlev = GV%ke, &
924		energy_flux = US%RZ3_T3_to_W_m2*tidal_qe_md(:), &
925		SchmittnerCoeff = Schmittner_coeff, &
926		exp_hab_zetar = exp_hab_zetar, &
927	0	CVmix_tidal_params_user = CS%CVMix_tidal_params)
928
929		! XXX: Temporary de-scaling of N2_int(i,:) into a temporary variable
930	0	do k=1,GV%ke+1
931	0	N2_int_i(k) = US%s_to_T*2 N2_int(i,k)
932		enddo
933
934		call CVMix_coeffs_tidal_schmittner( Mdiff_out = Kv_tidal, &
935		Tdiff_out = Kd_tidal, &
936		Nsqr = N2_int_i, &
937		OceanDepth = -iFaceHeight(GV%ke+1), &
938		nlev = GV%ke, &
939		max_nlev = GV%ke, &
940		SchmittnerCoeff = Schmittner_coeff, &
941		SchmittnerSouthernOcean = SchmittnerSocn, &
942		CVmix_params = CS%CVMix_glb_params, &
943	0	CVmix_tidal_params_user = CS%CVMix_tidal_params)
944
945		! Update diffusivity
946	0	if (present(Kd_lay)) then
947	0	do k=1,GV%ke
948	0	Kd_lay(i,k) = Kd_lay(i,k) + 0.5 * GV%m2_s_to_HZ_T * (Kd_tidal(k) + Kd_tidal(k+1))
949		enddo
950		endif
951	0	if (present(Kd_int)) then
952	0	do K=1,GV%ke+1
953	0	Kd_int(i,K) = Kd_int(i,K) + (GV%m2_s_to_HZ_T * Kd_tidal(K))
954		enddo
955		endif
956
957		! Update viscosity
958	0	if (associated(Kv)) then
959	0	do K=1,GV%ke+1
960	0	Kv(i,j,K) = Kv(i,j,K) + GV%m2_s_to_HZ_T * Kv_tidal(K) ! Rescale from m2 s-1 to H Z T-1.
961		enddo
962		endif
963
964		! diagnostics
965	0	if (allocated(CS%dd%Kd_itidal)) then
966	0	CS%dd%Kd_itidal(i,j,:) = GV%m2_s_to_HZ_T*Kd_tidal(:)
967		endif
968	0	if (allocated(CS%dd%N2_int)) then
969	0	CS%dd%N2_int(i,j,:) = N2_int(i,:)
970		endif
971	0	if (allocated(CS%dd%Schmittner_coeff_3d)) then
972	0	CS%dd%Schmittner_coeff_3d(i,j,:) = Schmittner_coeff(:)
973		endif
974	0	if (allocated(CS%dd%tidal_qe_md)) then
975	0	CS%dd%tidal_qe_md(i,j,:) = tidal_qe_md(:)
976		endif
977	0	if (allocated(CS%dd%vert_dep_3d)) then
978	0	CS%dd%vert_dep_3d(i,j,:) = vert_dep(:)
979		endif
980		enddo ! i=is,ie
981
982	0	deallocate(exp_hab_zetar)
983
984		case default
985		call MOM_error(FATAL, "tidal_mixing_init: Unrecognized setting "// &
986	0	"#define CVMIX_TIDAL_SCHEME found in input file.")
987		end select
988
989	0	end subroutine calculate_CVMix_tidal
990
991
992		!> This subroutine adds the effect of internal-tide-driven mixing to the layer diffusivities.
993		!! The mechanisms considered are (1) local dissipation of internal waves generated by the
994		!! barotropic flow ("itidal"), (2) local dissipation of internal waves generated by the propagating
995		!! low modes (rays) of the internal tide ("lowmode"), and (3) local dissipation of internal lee waves.
996		!! Will eventually need to add diffusivity due to other wave-breaking processes (e.g. Bottom friction,
997		!! Froude-number-depending breaking, PSI, etc.).
998	0	subroutine add_int_tide_diffusivity(dz, j, N2_bot, Rho_bot, N2_lay, TKE_to_Kd, max_TKE, &
999	0	G, GV, US, CS, Kd_max, Kd_lay, Kd_int, VBF)
1000		type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure
1001		type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure
1002		type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
1003		real, dimension(SZI_(G),SZK_(GV)), intent(in) :: dz !< The vertical distance across layers [Z ~> m]
1004		integer, intent(in) :: j !< The j-index to work on
1005		real, dimension(SZI_(G)), intent(in) :: N2_bot !< The near-bottom squared buoyancy frequency
1006		!! frequency [T-2 ~> s-2].
1007		real, dimension(SZI_(G)), intent(in) :: Rho_bot !< The near-bottom in situ density [R ~> kg m-3]
1008		real, dimension(SZI_(G),SZK_(GV)), intent(in) :: N2_lay !< The squared buoyancy frequency of the
1009		!! layers [T-2 ~> s-2].
1010		real, dimension(SZI_(G),SZK_(GV)), intent(in) :: TKE_to_Kd !< The conversion rate between the TKE
1011		!! dissipated within a layer and the
1012		!! diapycnal diffusivity within that layer,
1013		!! usually (~Rho_0 / (G_Earth * dRho_lay))
1014		!! [T2 Z-1 ~> s2 m-1]
1015		real, dimension(SZI_(G),SZK_(GV)), intent(in) :: max_TKE !< The energy required for a layer
1016		!! to entrain to its maximum realizable
1017		!! thickness [H Z2 T-3 ~> m3 s-3 or W m-2]
1018		type(tidal_mixing_cs), intent(inout) :: CS !< The control structure for this module
1019		real, intent(in) :: Kd_max !< The maximum increment for diapycnal
1020		!! diffusivity due to TKE-based processes
1021		!! [H Z T-1 ~> m2 s-1 or kg m-1 s-1].
1022		!! Set this to a negative value to have no limit.
1023		real, dimension(SZI_(G),SZK_(GV)), &
1024		optional, intent(inout) :: Kd_lay !< The diapycnal diffusivity in layers
1025		!! [H Z T-1 ~> m2 s-1 or kg m-1 s-1].
1026		real, dimension(SZI_(G),SZK_(GV)+1), &
1027		optional, intent(inout) :: Kd_int !< The diapycnal diffusivity at interfaces
1028		!! [H Z T-1 ~> m2 s-1 or kg m-1 s-1].
1029		type(vbf_CS), pointer :: VBF !< A diagnostics structure for vertical buoyancy fluxes
1030
1031		! local
1032
1033		real, dimension(SZI_(G)) :: &
1034	0	dztot, & ! Vertical distance between the top and bottom of the ocean [Z ~> m]
1035	0	dztot_WKB, & ! WKB scaled distance from top to bottom [Z ~> m]
1036	0	TKE_itidal_bot, & ! internal tide TKE at ocean bottom [H Z2 T-3 ~> m3 s-3 or W m-2]
1037	0	TKE_Niku_bot, & ! lee-wave TKE at ocean bottom [H Z2 T-3 ~> m3 s-3 or W m-2]
1038	0	TKE_lowmode_bot, & ! internal tide TKE at ocean bottom lost from all remote low modes [H Z2 T-3 ~> m3 s-3 or W m-2]
1039	0	Inv_int, & ! inverse of TKE decay for int tide over the depth of the ocean [nondim]
1040	0	Inv_int_lee, & ! inverse of TKE decay for lee waves over the depth of the ocean [nondim]
1041	0	Inv_int_low, & ! inverse of TKE decay for low modes over the depth of the ocean [nondim]
1042	0	z0_Polzin, & ! TKE decay scale in Polzin formulation [Z ~> m]
1043	0	z0_Polzin_scaled, & ! TKE decay scale in Polzin formulation [Z ~> m].
1044		! multiplied by N2_bot/N2_meanz to be coherent with the WKB scaled z
1045		! z=int(N2/N2_bot) N2_bot/N2_meanz = int(N2/N2_meanz)
1046		! z0_Polzin_scaled = z0_Polzin * N2_bot/N2_meanz
1047	0	N2_meanz, & ! vertically averaged squared buoyancy frequency [T-2 ~> s-2] for WKB scaling
1048	0	TKE_itidal_rem, & ! remaining internal tide TKE (from barotropic source) [H Z2 T-3 ~> m3 s-3 or W m-2]
1049	0	TKE_Niku_rem, & ! remaining lee-wave TKE [H Z2 T-3 ~> m3 s-3 or W m-2]
1050	0	TKE_lowmode_rem, & ! remaining internal tide TKE (from propagating low mode source) [H Z2 T-3 ~> m3 s-3 or W m-2]
1051	0	TKE_frac_top, & ! fraction of bottom TKE that should appear at top of a layer [nondim]
1052	0	TKE_frac_top_lee, & ! fraction of bottom TKE that should appear at top of a layer [nondim]
1053	0	TKE_frac_top_lowmode, &
1054		! fraction of bottom TKE that should appear at top of a layer [nondim]
1055	0	z_from_bot, & ! distance from bottom [Z ~> m]
1056	0	z_from_bot_WKB ! WKB scaled distance from bottom [Z ~> m]
1057
1058		real :: Kd_add ! Diffusivity to add in a layer [H Z T-1 ~> m2 s-1 or kg m-1 s-1]
1059		real :: TKE_itide_lay ! internal tide TKE imparted to a layer (from barotropic) [H Z2 T-3 ~> m3 s-3 or W m-2]
1060		real :: TKE_Niku_lay ! lee-wave TKE imparted to a layer [H Z2 T-3 ~> m3 s-3 or W m-2]
1061		real :: TKE_lowmode_lay ! internal tide TKE imparted to a layer (from low mode) [H Z2 T-3 ~> m3 s-3 or W m-2]
1062		real :: frac_used ! fraction of TKE that can be used in a layer [nondim]
1063		real :: Izeta ! inverse of TKE decay scale [Z-1 ~> m-1]
1064		real :: Izeta_lee ! inverse of TKE decay scale for lee waves [Z-1 ~> m-1]
1065		real :: z0Ps_num ! The numerator of the unlimited z0_Polzin_scaled [Z T-3 ~> m s-3]
1066		real :: z0Ps_denom ! The denominator of the unlimited z0_Polzin_scaled [T-3 ~> s-3].
1067		real :: z0_psl ! temporary variable [Z ~> m]
1068		real :: TKE_lowmode_tot ! TKE from all low modes [R Z3 T-3 ~> W m-2]
1069
1070		logical :: use_Polzin, use_Simmons
1071		integer :: i, k, is, ie, nz
1072
1073	0	is = G%isc ; ie = G%iec ; nz = GV%ke
1074
1075	0	if (.not.(CS%Int_tide_dissipation .or. CS%Lee_wave_dissipation)) return
1076
1077	0	do i=is,ie ; dztot(i) = 0.0 ; Inv_int(i) = 0.0 ; Inv_int_lee(i) = 0.0 ; Inv_int_low(i) = 0.0 ; enddo
1078	0	do k=1,nz ; do i=is,ie
1079	0	dztot(i) = dztot(i) + dz(i,k)
1080		enddo ; enddo
1081
1082		use_Polzin = ((CS%Int_tide_dissipation .and. (CS%int_tide_profile == POLZIN_09)) .or. &
1083		(CS%lee_wave_dissipation .and. (CS%lee_wave_profile == POLZIN_09)) .or. &
1084	0	(CS%Lowmode_itidal_dissipation .and. (CS%int_tide_profile == POLZIN_09)))
1085		use_Simmons = ((CS%Int_tide_dissipation .and. (CS%int_tide_profile == STLAURENT_02)) .or. &
1086		(CS%lee_wave_dissipation .and. (CS%lee_wave_profile == STLAURENT_02)) .or. &
1087	0	(CS%Lowmode_itidal_dissipation .and. (CS%int_tide_profile == STLAURENT_02)))
1088
1089		! Calculate parameters for vertical structure of dissipation
1090		! Simmons:
1091	0	if ( use_Simmons ) then
1092	0	Izeta = 1.0 / max(CS%Int_tide_decay_scale, GV%dz_subroundoff)
1093	0	Izeta_lee = 1.0 / max(CS%Int_tide_decay_scale*CS%Decay_scale_factor_lee, GV%dz_subroundoff)
1094	0	do i=is,ie
1095	0	CS%Nb(i,j) = sqrt(N2_bot(i))
1096	0	if (allocated(CS%dd%N2_bot)) &
1097	0	CS%dd%N2_bot(i,j) = N2_bot(i)
1098	0	if ( CS%Int_tide_dissipation ) then
1099	0	if (Izeta*dztot(i) > 1.0e-14) then ! L'Hospital's version of Adcroft's reciprocal rule.
1100	0	Inv_int(i) = 1.0 / (1.0 - exp(-Izeta*dztot(i)))
1101		endif
1102		endif
1103	0	if ( CS%Lee_wave_dissipation ) then
1104	0	if (Izeta_lee*dztot(i) > 1.0e-14) then ! L'Hospital's version of Adcroft's reciprocal rule.
1105	0	Inv_int_lee(i) = 1.0 / (1.0 - exp(-Izeta_lee*dztot(i)))
1106		endif
1107		endif
1108	0	if ( CS%Lowmode_itidal_dissipation) then
1109	0	if (Izeta*dztot(i) > 1.0e-14) then ! L'Hospital's version of Adcroft's reciprocal rule.
1110	0	Inv_int_low(i) = 1.0 / (1.0 - exp(-Izeta*dztot(i)))
1111		endif
1112		endif
1113	0	z_from_bot(i) = dz(i,nz)
1114		enddo
1115		endif ! Simmons
1116
1117		! Polzin:
1118	0	if ( use_Polzin ) then
1119		! WKB scaling of the vertical coordinate
1120	0	do i=is,ie ; N2_meanz(i) = 0.0 ; enddo
1121	0	do k=1,nz ; do i=is,ie
1122	0	N2_meanz(i) = N2_meanz(i) + N2_lay(i,k) * dz(i,k)
1123		enddo ; enddo
1124	0	do i=is,ie
1125	0	N2_meanz(i) = N2_meanz(i) / (dztot(i) + GV%dz_subroundoff)
1126	0	if (allocated(CS%dd%N2_meanz)) &
1127	0	CS%dd%N2_meanz(i,j) = N2_meanz(i)
1128		enddo
1129
1130		! WKB scaled z(z=H) z at the surface using the modified Polzin WKB scaling
1131	0	do i=is,ie ; dztot_WKB(i) = dztot(i) ; enddo
1132		! do i=is,ie ; dztot_WKB(i) = 0.0 ; enddo
1133		! do k=1,nz ; do i=is,ie
1134		! dztot_WKB(i) = dztot_WKB(i) + dz(i,k) * N2_lay(i,k) / N2_meanz(i)
1135		! enddo ; enddo
1136		! dztot_WKB(i) = dztot(i) ! Nearly equivalent and simpler
1137
1138	0	do i=is,ie
1139	0	CS%Nb(i,j) = sqrt(N2_bot(i))
1140	0	if (CS%tidal_answer_date < 20190101) then
1141	0	if ((CS%tideamp(i,j) > 0.0) .and. &
1142		(CS%kappa_itides*2 CS%h2(i,j) * CS%Nb(i,j)*3 > 1.0e-14US%T_to_s**3) ) then
1143		z0_Polzin(i) = CS%Polzin_decay_scale_factor * CS%Nu_Polzin * &
1144		CS%Nbotref_Polzin*2 CS%tideamp(i,j) / &
1145	0	( CS%kappa_itides*2 CS%h2(i,j) * CS%Nb(i,j)**3 )
1146	0	if (z0_Polzin(i) < CS%Polzin_min_decay_scale) &
1147	0	z0_Polzin(i) = CS%Polzin_min_decay_scale
1148	0	if (N2_meanz(i) > 1.0e-14US%T_to_s*2 ) then
1149	0	z0_Polzin_scaled(i) = z0_Polzin(i)CS%Nb(i,j)*2 / N2_meanz(i)
1150		else
1151	0	z0_Polzin_scaled(i) = CS%Polzin_decay_scale_max_factor * dztot(i)
1152		endif
1153	0	if (z0_Polzin_scaled(i) > (CS%Polzin_decay_scale_max_factor * dztot(i)) ) &
1154	0	z0_Polzin_scaled(i) = CS%Polzin_decay_scale_max_factor * dztot(i)
1155		else
1156	0	z0_Polzin(i) = CS%Polzin_decay_scale_max_factor * dztot(i)
1157	0	z0_Polzin_scaled(i) = CS%Polzin_decay_scale_max_factor * dztot(i)
1158		endif
1159		else
1160	0	z0Ps_num = (CS%Polzin_decay_scale_factor * CS%Nu_Polzin * CS%Nbotref_Polzin*2) CS%tideamp(i,j)
1161	0	z0Ps_denom = ( CS%kappa_itides*2 CS%h2(i,j) * CS%Nb(i,j) * N2_meanz(i) )
1162	0	if ((CS%tideamp(i,j) > 0.0) .and. &
1163		(z0Ps_num < z0Ps_denom * CS%Polzin_decay_scale_max_factor * dztot(i))) then
1164	0	z0_Polzin_scaled(i) = z0Ps_num / z0Ps_denom
1165
1166	0	if (abs(N2_meanz(i) * z0_Polzin_scaled(i)) < &
1167		CS%Nb(i,j)*2 (CS%Polzin_decay_scale_max_factor * dztot(i))) then
1168	0	z0_Polzin(i) = z0_Polzin_scaled(i) * (N2_meanz(i) / CS%Nb(i,j)**2)
1169		else
1170	0	z0_Polzin(i) = CS%Polzin_decay_scale_max_factor * dztot(i)
1171		endif
1172		else
1173	0	z0_Polzin(i) = CS%Polzin_decay_scale_max_factor * dztot(i)
1174	0	z0_Polzin_scaled(i) = CS%Polzin_decay_scale_max_factor * dztot(i)
1175		endif
1176		endif
1177
1178	0	if (allocated(CS%dd%Polzin_decay_scale)) &
1179	0	CS%dd%Polzin_decay_scale(i,j) = z0_Polzin(i)
1180	0	if (allocated(CS%dd%Polzin_decay_scale_scaled)) &
1181	0	CS%dd%Polzin_decay_scale_scaled(i,j) = z0_Polzin_scaled(i)
1182	0	if (allocated(CS%dd%N2_bot)) &
1183	0	CS%dd%N2_bot(i,j) = CS%Nb(i,j)*CS%Nb(i,j)
1184
1185	0	if (CS%tidal_answer_date < 20190101) then
1186		! These expressions use dimensional constants to avoid NaN values.
1187	0	if ( CS%Int_tide_dissipation .and. (CS%int_tide_profile == POLZIN_09) ) then
1188	0	if (dztot_WKB(i) > 1.0e-14*US%m_to_Z) &
1189	0	Inv_int(i) = ( z0_Polzin_scaled(i) / dztot_WKB(i) ) + 1.0
1190		endif
1191	0	if ( CS%lee_wave_dissipation .and. (CS%lee_wave_profile == POLZIN_09) ) then
1192	0	if (dztot_WKB(i) > 1.0e-14*US%m_to_Z) &
1193	0	Inv_int_lee(i) = ( z0_Polzin_scaled(i)*CS%Decay_scale_factor_lee / dztot_WKB(i) ) + 1.0
1194		endif
1195	0	if ( CS%Lowmode_itidal_dissipation .and. (CS%int_tide_profile == POLZIN_09) ) then
1196	0	if (dztot_WKB(i) > 1.0e-14*US%m_to_Z) &
1197	0	Inv_int_low(i) = ( z0_Polzin_scaled(i) / dztot_WKB(i) ) + 1.0
1198		endif
1199		else
1200		! These expressions give values of Inv_int < 10^14 using a variant of Adcroft's reciprocal rule.
1201	0	Inv_int(i) = 0.0 ; Inv_int_lee(i) = 0.0 ; Inv_int_low(i) = 0.0
1202	0	if ( CS%Int_tide_dissipation .and. (CS%int_tide_profile == POLZIN_09) ) then
1203	0	if (z0_Polzin_scaled(i) < 1.0e14 * dztot_WKB(i)) &
1204	0	Inv_int(i) = ( z0_Polzin_scaled(i) / dztot_WKB(i) ) + 1.0
1205		endif
1206	0	if ( CS%lee_wave_dissipation .and. (CS%lee_wave_profile == POLZIN_09) ) then
1207	0	if (z0_Polzin_scaled(i) < 1.0e14 * dztot_WKB(i)) &
1208	0	Inv_int_lee(i) = ( z0_Polzin_scaled(i)*CS%Decay_scale_factor_lee / dztot_WKB(i) ) + 1.0
1209		endif
1210	0	if ( CS%Lowmode_itidal_dissipation .and. (CS%int_tide_profile == POLZIN_09) ) then
1211	0	if (z0_Polzin_scaled(i) < 1.0e14 * dztot_WKB(i)) &
1212	0	Inv_int_low(i) = ( z0_Polzin_scaled(i) / dztot_WKB(i) ) + 1.0
1213		endif
1214		endif
1215
1216	0	z_from_bot(i) = dz(i,nz)
1217		! Use the new formulation for WKB scaling. N2 is referenced to its vertical mean.
1218	0	if (CS%tidal_answer_date < 20190101) then
1219	0	if (N2_meanz(i) > 1.0e-14US%T_to_s*2 ) then
1220	0	z_from_bot_WKB(i) = dz(i,nz) * N2_lay(i,nz) / N2_meanz(i)
1221	0	else ; z_from_bot_WKB(i) = 0 ; endif
1222		else
1223	0	if (dz(i,nz) * N2_lay(i,nz) < N2_meanz(i) * (1.0e14 * dztot_WKB(i))) then
1224	0	z_from_bot_WKB(i) = dz(i,nz) * N2_lay(i,nz) / N2_meanz(i)
1225	0	else ; z_from_bot_WKB(i) = 0 ; endif
1226		endif
1227		enddo
1228		endif ! Polzin
1229
1230		! Calculate/get dissipation values at bottom
1231		! Both Polzin and Simmons:
1232	0	do i=is,ie
1233		! Dissipation of locally trapped internal tide (non-propagating high modes)
1234	0	if (GV%Boussinesq .or. GV%semi_Boussinesq) then
1235	0	TKE_itidal_bot(i) = min(GV%Z_to_HCS%TKE_itidal(i,j)CS%Nb(i,j), CS%TKE_itide_max)
1236		else
1237		TKE_itidal_bot(i) = min(GV%RZ_to_HRho_bot(i) (CS%TKE_itidal(i,j)*CS%Nb(i,j)), &
1238	0	CS%TKE_itide_max)
1239		endif
1240	0	if (allocated(CS%dd%TKE_itidal_used)) &
1241	0	CS%dd%TKE_itidal_used(i,j) = TKE_itidal_bot(i)
1242	0	TKE_itidal_bot(i) = (GV%RZ_to_H * CS%Mu_itides * CS%Gamma_itides) * TKE_itidal_bot(i)
1243		! Dissipation of locally trapped lee waves
1244	0	TKE_Niku_bot(i) = 0.0
1245	0	if (CS%Lee_wave_dissipation) then
1246	0	TKE_Niku_bot(i) = (GV%RZ_to_H * CS%Mu_itides * CS%Gamma_lee) * CS%TKE_Niku(i,j)
1247		endif
1248		! Dissipation of propagating internal tide (baroclinic low modes; rays) (BDM)
1249	0	TKE_lowmode_tot = 0.0
1250	0	TKE_lowmode_bot(i) = 0.0
1251	0	if (CS%Lowmode_itidal_dissipation) then
1252		! get loss rate due to wave drag on low modes (already multiplied by q)
1253	0	call get_lowmode_loss(i,j,G,CS%int_tide_CSp,"WaveDrag",TKE_lowmode_tot)
1254	0	TKE_lowmode_bot(i) = CS%Mu_itides * GV%RZ_to_H * TKE_lowmode_tot
1255		endif
1256		! Vertical energy flux at bottom
1257	0	TKE_itidal_rem(i) = Inv_int(i) * TKE_itidal_bot(i)
1258	0	TKE_Niku_rem(i) = Inv_int_lee(i) * TKE_Niku_bot(i)
1259	0	TKE_lowmode_rem(i) = Inv_int_low(i) * TKE_lowmode_bot(i)
1260
1261	0	if (allocated(CS%dd%Fl_itidal)) &
1262	0	CS%dd%Fl_itidal(i,j,nz) = TKE_itidal_rem(i) !why is this here? BDM
1263		enddo
1264
1265		! Estimate the work that would be done by mixing in each layer.
1266		! Simmons:
1267	0	if ( use_Simmons ) then
1268	0	do k=nz-1,2,-1 ; do i=is,ie
1269	0	if (max_TKE(i,k) <= 0.0) cycle
1270	0	z_from_bot(i) = z_from_bot(i) + dz(i,k)
1271
1272		! Fraction of bottom flux predicted to reach top of this layer
1273	0	TKE_frac_top(i) = Inv_int(i) * exp(-Izeta * z_from_bot(i))
1274	0	TKE_frac_top_lee(i) = Inv_int_lee(i) * exp(-Izeta_lee * z_from_bot(i))
1275	0	TKE_frac_top_lowmode(i) = Inv_int_low(i) * exp(-Izeta * z_from_bot(i))
1276
1277		! Actual influx at bottom of layer minus predicted outflux at top of layer to give
1278		! predicted power expended
1279	0	TKE_itide_lay = TKE_itidal_rem(i) - TKE_itidal_bot(i) * TKE_frac_top(i)
1280	0	TKE_Niku_lay = TKE_Niku_rem(i) - TKE_Niku_bot(i) * TKE_frac_top_lee(i)
1281	0	TKE_lowmode_lay = TKE_lowmode_rem(i) - TKE_lowmode_bot(i)* TKE_frac_top_lowmode(i)
1282
1283		! Actual power expended may be less than predicted if stratification is weak; adjust
1284	0	if (TKE_itide_lay + TKE_Niku_lay + TKE_lowmode_lay > max_TKE(i,k)) then
1285	0	frac_used = (max_TKE(i,k)) / (TKE_itide_lay + TKE_Niku_lay + TKE_lowmode_lay)
1286	0	TKE_itide_lay = frac_used * TKE_itide_lay
1287	0	TKE_Niku_lay = frac_used * TKE_Niku_lay
1288	0	TKE_lowmode_lay = frac_used * TKE_lowmode_lay
1289		endif
1290
1291		! Calculate vertical flux available to bottom of layer above
1292	0	TKE_itidal_rem(i) = TKE_itidal_rem(i) - TKE_itide_lay
1293	0	TKE_Niku_rem(i) = TKE_Niku_rem(i) - TKE_Niku_lay
1294	0	TKE_lowmode_rem(i) = TKE_lowmode_rem(i) - TKE_lowmode_lay
1295
1296		! Convert power to diffusivity
1297	0	Kd_add = TKE_to_Kd(i,k) * (TKE_itide_lay + TKE_Niku_lay + TKE_lowmode_lay)
1298
1299	0	if (Kd_max >= 0.0) Kd_add = min(Kd_add, Kd_max)
1300	0	if (present(Kd_lay)) then
1301	0	Kd_lay(i,k) = Kd_lay(i,k) + Kd_add
1302		endif
1303
1304	0	if (present(Kd_int)) then
1305	0	Kd_int(i,K) = Kd_int(i,K) + 0.5 * Kd_add
1306	0	Kd_int(i,K+1) = Kd_int(i,K+1) + 0.5 * Kd_add
1307		endif
1308
1309		! diagnostics
1310	0	if (allocated(CS%dd%Kd_itidal).or.(associated(VBF%Kd_itides))) then
1311		! If at layers, CS%dd%Kd_itidal is just TKE_to_Kd(i,k) * TKE_itide_lay
1312		! The following sets the interface diagnostics.
1313	0	Kd_add = TKE_to_Kd(i,k) * TKE_itide_lay
1314	0	if (Kd_max >= 0.0) Kd_add = min(Kd_add, Kd_max)
1315	0	if (allocated(CS%dd%Kd_itidal)) then
1316	0	if (k>1) CS%dd%Kd_itidal(i,j,K) = CS%dd%Kd_itidal(i,j,K) + 0.5*Kd_add
1317	0	if (k<nz) CS%dd%Kd_itidal(i,j,K+1) = CS%dd%Kd_itidal(i,j,K+1) + 0.5*Kd_add
1318		endif
1319	0	if (associated(VBF%Kd_itides)) then
1320		!Not to be confused w/ Kd_itidal (this is to be consistent w/ output parameter names)
1321	0	if (k>1) VBF%Kd_itides(i,j,K) = VBF%Kd_itides(i,j,K) + 0.5*Kd_add
1322	0	if (k<nz) VBF%Kd_itides(i,j,K+1) = VBF%Kd_itides(i,j,K+1) + 0.5*Kd_add
1323		endif
1324		endif
1325	0	if (allocated(CS%dd%Kd_Itidal_work)) &
1326	0	CS%dd%Kd_itidal_work(i,j,k) = GV%H_to_RZ * TKE_itide_lay
1327	0	if (allocated(CS%dd%Fl_itidal)) &
1328	0	CS%dd%Fl_itidal(i,j,k) = TKE_itidal_rem(i)
1329
1330	0	if (allocated(CS%dd%Kd_Niku).or.(associated(VBF%Kd_Niku))) then
1331		! If at layers, CS%dd%Kd_Niku(i,j,K) is just TKE_to_Kd(i,k) * TKE_Niku_lay
1332		! The following sets the interface diagnostics.
1333	0	Kd_add = TKE_to_Kd(i,k) * TKE_Niku_lay
1334	0	if (Kd_max >= 0.0) Kd_add = min(Kd_add, Kd_max)
1335	0	if (allocated(CS%dd%Kd_Niku)) then
1336	0	if (k>1) CS%dd%Kd_Niku(i,j,K) = CS%dd%Kd_Niku(i,j,K) + 0.5*Kd_add
1337	0	if (k<nz) CS%dd%Kd_Niku(i,j,K+1) = CS%dd%Kd_Niku(i,j,K+1) + 0.5*Kd_add
1338		endif
1339	0	if (associated(VBF%Kd_Niku)) then
1340	0	if (k>1) VBF%Kd_Niku(i,j,K) = VBF%Kd_Niku(i,j,K) + 0.5*Kd_add
1341	0	if (k<nz) VBF%Kd_Niku(i,j,K+1) = VBF%Kd_Niku(i,j,K+1) + 0.5*Kd_add
1342		endif
1343		endif
1344		! if (associated(CS%dd%Kd_Niku)) CS%dd%Kd_Niku(i,j,K) = TKE_to_Kd(i,k) * TKE_Niku_lay
1345	0	if (allocated(CS%dd%Kd_Niku_work)) &
1346	0	CS%dd%Kd_Niku_work(i,j,k) = GV%H_to_RZ * TKE_Niku_lay
1347
1348	0	if (allocated(CS%dd%Kd_lowmode).or.(associated(VBF%Kd_lowmode))) then
1349		! If at layers, CS%dd%Kd_lowmode is just TKE_to_Kd(i,k) * TKE_lowmode_lay
1350		! The following sets the interface diagnostics.
1351	0	Kd_add = TKE_to_Kd(i,k) * TKE_lowmode_lay
1352	0	if (Kd_max >= 0.0) Kd_add = min(Kd_add, Kd_max)
1353	0	if (allocated(CS%dd%Kd_lowmode)) then
1354	0	if (k>1) CS%dd%Kd_lowmode(i,j,K) = CS%dd%Kd_lowmode(i,j,K) + 0.5*Kd_add
1355	0	if (k<nz) CS%dd%Kd_lowmode(i,j,K+1) = CS%dd%Kd_lowmode(i,j,K+1) + 0.5*Kd_add
1356		endif
1357	0	if (associated(VBF%Kd_lowmode)) then
1358	0	if (k>1) VBF%Kd_lowmode(i,j,K) = VBF%Kd_lowmode(i,j,K) + 0.5*Kd_add
1359	0	if (k<nz) VBF%Kd_lowmode(i,j,K+1) = VBF%Kd_lowmode(i,j,K+1) + 0.5*Kd_add
1360		endif
1361		endif
1362	0	if (allocated(CS%dd%Kd_lowmode_work)) &
1363	0	CS%dd%Kd_lowmode_work(i,j,k) = GV%H_to_RZ * TKE_lowmode_lay
1364	0	if (allocated(CS%dd%Fl_lowmode)) &
1365	0	CS%dd%Fl_lowmode(i,j,k) = TKE_lowmode_rem(i)
1366		enddo ; enddo
1367		endif ! Simmons
1368
1369		! Polzin:
1370	0	if ( use_Polzin ) then
1371	0	do k=nz-1,2,-1 ; do i=is,ie
1372	0	if (max_TKE(i,k) <= 0.0) cycle
1373	0	z_from_bot(i) = z_from_bot(i) + dz(i,k)
1374	0	if (CS%tidal_answer_date < 20190101) then
1375	0	if (N2_meanz(i) > 1.0e-14US%T_to_s*2 ) then
1376	0	z_from_bot_WKB(i) = z_from_bot_WKB(i) + dz(i,k) * N2_lay(i,k) / N2_meanz(i)
1377	0	else ; z_from_bot_WKB(i) = 0 ; endif
1378		else
1379	0	if (dz(i,k) * N2_lay(i,k) < (1.0e14 * dztot_WKB(i)) * N2_meanz(i)) then
1380	0	z_from_bot_WKB(i) = z_from_bot_WKB(i) + dz(i,k) * N2_lay(i,k) / N2_meanz(i)
1381		endif
1382		endif
1383
1384		! Fraction of bottom flux predicted to reach top of this layer
1385		TKE_frac_top(i) = ( Inv_int(i) * z0_Polzin_scaled(i) ) / &
1386	0	( z0_Polzin_scaled(i) + z_from_bot_WKB(i) )
1387	0	z0_psl = z0_Polzin_scaled(i)*CS%Decay_scale_factor_lee
1388	0	TKE_frac_top_lee(i) = (Inv_int_lee(i) * z0_psl) / (z0_psl + z_from_bot_WKB(i))
1389		TKE_frac_top_lowmode(i) = ( Inv_int_low(i) * z0_Polzin_scaled(i) ) / &
1390	0	( z0_Polzin_scaled(i) + z_from_bot_WKB(i) )
1391
1392		! Actual influx at bottom of layer minus predicted outflux at top of layer to give
1393		! predicted power expended
1394	0	TKE_itide_lay = TKE_itidal_rem(i) - TKE_itidal_bot(i) *TKE_frac_top(i)
1395	0	TKE_Niku_lay = TKE_Niku_rem(i) - TKE_Niku_bot(i) * TKE_frac_top_lee(i)
1396	0	TKE_lowmode_lay = TKE_lowmode_rem(i) - TKE_lowmode_bot(i)*TKE_frac_top_lowmode(i)
1397
1398		! Actual power expended may be less than predicted if stratification is weak; adjust
1399	0	if (TKE_itide_lay + TKE_Niku_lay + TKE_lowmode_lay > max_TKE(i,k)) then
1400	0	frac_used = max_TKE(i,k) / (TKE_itide_lay + TKE_Niku_lay + TKE_lowmode_lay)
1401	0	TKE_itide_lay = frac_used * TKE_itide_lay
1402	0	TKE_Niku_lay = frac_used * TKE_Niku_lay
1403	0	TKE_lowmode_lay = frac_used * TKE_lowmode_lay
1404		endif
1405
1406		! Calculate vertical flux available to bottom of layer above
1407	0	TKE_itidal_rem(i) = TKE_itidal_rem(i) - TKE_itide_lay
1408	0	TKE_Niku_rem(i) = TKE_Niku_rem(i) - TKE_Niku_lay
1409	0	TKE_lowmode_rem(i) = TKE_lowmode_rem(i) - TKE_lowmode_lay
1410
1411		! Convert power to diffusivity
1412	0	Kd_add = TKE_to_Kd(i,k) * (TKE_itide_lay + TKE_Niku_lay + TKE_lowmode_lay)
1413
1414	0	if (Kd_max >= 0.0) Kd_add = min(Kd_add, Kd_max)
1415	0	if (present(Kd_lay)) then
1416	0	Kd_lay(i,k) = Kd_lay(i,k) + Kd_add
1417		endif
1418
1419	0	if (present(Kd_int)) then
1420	0	Kd_int(i,K) = Kd_int(i,K) + 0.5 * Kd_add
1421	0	Kd_int(i,K+1) = Kd_int(i,K+1) + 0.5 * Kd_add
1422		endif
1423
1424		! diagnostics
1425	0	if (allocated(CS%dd%Kd_itidal).or.(associated(VBF%Kd_itides))) then
1426		! If at layers, this is just CS%dd%Kd_itidal(i,j,K) = TKE_to_Kd(i,k) * TKE_itide_lay
1427		! The following sets the interface diagnostics.
1428	0	Kd_add = TKE_to_Kd(i,k) * TKE_itide_lay
1429	0	if (Kd_max >= 0.0) Kd_add = min(Kd_add, Kd_max)
1430	0	if (allocated(CS%dd%Kd_itidal)) then
1431	0	if (k>1) CS%dd%Kd_itidal(i,j,K) = CS%dd%Kd_itidal(i,j,K) + 0.5*Kd_add
1432	0	if (k<nz) CS%dd%Kd_itidal(i,j,K+1) = CS%dd%Kd_itidal(i,j,K+1) + 0.5*Kd_add
1433		endif
1434	0	if (associated(VBF%Kd_itides)) then
1435		!Not to be confused w/ Kd_itidal (this is to be consistent w/ output parameter names)
1436	0	if (k>1) VBF%Kd_itides(i,j,K) = VBF%Kd_itides(i,j,K) + 0.5*Kd_add
1437	0	if (k<nz) VBF%Kd_itides(i,j,K+1) = VBF%Kd_itides(i,j,K+1) + 0.5*Kd_add
1438		endif
1439		endif
1440	0	if (allocated(CS%dd%Kd_Itidal_work)) &
1441	0	CS%dd%Kd_itidal_work(i,j,k) = GV%H_to_RZ * TKE_itide_lay
1442	0	if (allocated(CS%dd%Fl_itidal)) CS%dd%Fl_itidal(i,j,k) = TKE_itidal_rem(i)
1443
1444	0	if (allocated(CS%dd%Kd_Niku).or.(associated(VBF%Kd_Niku))) then
1445		! If at layers, this is just CS%dd%Kd_Niku(i,j,K) = TKE_to_Kd(i,k) * TKE_Niku_lay
1446		! The following sets the interface diagnostics.
1447	0	Kd_add = TKE_to_Kd(i,k) * TKE_Niku_lay
1448	0	if (Kd_max >= 0.0) Kd_add = min(Kd_add, Kd_max)
1449	0	if (allocated(CS%dd%Kd_Niku)) then
1450	0	if (k>1) CS%dd%Kd_Niku(i,j,K) = CS%dd%Kd_Niku(i,j,K) + 0.5*Kd_add
1451	0	if (k<nz) CS%dd%Kd_Niku(i,j,K+1) = CS%dd%Kd_Niku(i,j,K+1) + 0.5*Kd_add
1452		endif
1453	0	if (associated(VBF%Kd_Niku)) then
1454	0	if (k>1) VBF%Kd_Niku(i,j,K) = VBF%Kd_Niku(i,j,K) + 0.5*Kd_add
1455	0	if (k<nz) VBF%Kd_Niku(i,j,K+1) = VBF%Kd_Niku(i,j,K+1) + 0.5*Kd_add
1456		endif
1457		endif
1458		! if (associated(CS%dd%Kd_Niku)) CS%dd%Kd_Niku(i,j,K) = TKE_to_Kd(i,k) * TKE_Niku_lay
1459	0	if (allocated(CS%dd%Kd_Niku_work)) CS%dd%Kd_Niku_work(i,j,k) = GV%H_to_RZ * TKE_Niku_lay
1460
1461	0	if (allocated(CS%dd%Kd_lowmode).or.(associated(VBF%Kd_lowmode))) then
1462		! If at layers, CS%dd%Kd_lowmode is just TKE_to_Kd(i,k) * TKE_lowmode_lay
1463		! The following sets the interface diagnostics.
1464	0	Kd_add = TKE_to_Kd(i,k) * TKE_lowmode_lay
1465	0	if (Kd_max >= 0.0) Kd_add = min(Kd_add, Kd_max)
1466	0	if (allocated(CS%dd%Kd_lowmode)) then
1467	0	if (k>1) CS%dd%Kd_lowmode(i,j,K) = CS%dd%Kd_lowmode(i,j,K) + 0.5*Kd_add
1468	0	if (k<nz) CS%dd%Kd_lowmode(i,j,K+1) = CS%dd%Kd_lowmode(i,j,K+1) + 0.5*Kd_add
1469		endif
1470	0	if (associated(VBF%Kd_lowmode)) then
1471	0	if (k>1) VBF%Kd_lowmode(i,j,K) = VBF%Kd_lowmode(i,j,K) + 0.5*Kd_add
1472	0	if (k<nz) VBF%Kd_lowmode(i,j,K+1) = VBF%Kd_lowmode(i,j,K+1) + 0.5*Kd_add
1473		endif
1474		endif
1475	0	if (allocated(CS%dd%Kd_lowmode_work)) &
1476	0	CS%dd%Kd_lowmode_work(i,j,k) = GV%H_to_RZ * TKE_lowmode_lay
1477	0	if (allocated(CS%dd%Fl_lowmode)) CS%dd%Fl_lowmode(i,j,k) = TKE_lowmode_rem(i)
1478
1479		enddo ; enddo
1480		endif ! Polzin
1481
1482	0	end subroutine add_int_tide_diffusivity
1483
1484		!> Sets up diagnostics arrays for tidal mixing.
1485	0	subroutine setup_tidal_diagnostics(G, GV, CS)
1486		type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure
1487		type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure
1488		type(tidal_mixing_cs), intent(inout) :: CS !< The control structure for this module
1489
1490		! local
1491		integer :: isd, ied, jsd, jed, nz
1492
1493	0	isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed ; nz = GV%ke
1494
1495	0	if ((CS%id_Kd_itidal > 0) .or. (CS%id_Kd_Itidal_work > 0)) &
1496	0	allocate(CS%dd%Kd_itidal(isd:ied,jsd:jed,nz+1), source=0.0)
1497	0	if ((CS%id_Kd_lowmode > 0) .or. (CS%id_Kd_lowmode_work > 0)) &
1498	0	allocate(CS%dd%Kd_lowmode(isd:ied,jsd:jed,nz+1), source=0.0)
1499	0	if (CS%id_Fl_itidal > 0) allocate(CS%dd%Fl_itidal(isd:ied,jsd:jed,nz+1), source=0.0)
1500	0	if (CS%id_Fl_lowmode > 0) allocate(CS%dd%Fl_lowmode(isd:ied,jsd:jed,nz+1), source=0.0)
1501	0	if (CS%id_Polzin_decay_scale > 0) allocate(CS%dd%Polzin_decay_scale(isd:ied,jsd:jed), source=0.0)
1502	0	if (CS%id_N2_bot > 0) allocate(CS%dd%N2_bot(isd:ied,jsd:jed), source=0.0)
1503	0	if (CS%id_N2_meanz > 0) allocate(CS%dd%N2_meanz(isd:ied,jsd:jed), source=0.0)
1504	0	if (CS%id_Polzin_decay_scale_scaled > 0) &
1505	0	allocate(CS%dd%Polzin_decay_scale_scaled(isd:ied,jsd:jed), source=0.0)
1506	0	if ((CS%id_Kd_Niku > 0) .or. (CS%id_Kd_Niku_work > 0)) &
1507	0	allocate(CS%dd%Kd_Niku(isd:ied,jsd:jed,nz+1), source=0.0)
1508	0	if (CS%id_Kd_Niku_work > 0) allocate(CS%dd%Kd_Niku_work(isd:ied,jsd:jed,nz), source=0.0)
1509	0	if (CS%id_Kd_Itidal_work > 0) allocate(CS%dd%Kd_Itidal_work(isd:ied,jsd:jed,nz), source=0.0)
1510	0	if (CS%id_Kd_Lowmode_Work > 0) allocate(CS%dd%Kd_Lowmode_Work(isd:ied,jsd:jed,nz), source=0.0)
1511	0	if (CS%id_TKE_itidal > 0) allocate(CS%dd%TKE_Itidal_used(isd:ied,jsd:jed), source=0.)
1512		! additional diags for CVMix
1513	0	if (CS%id_N2_int > 0) allocate(CS%dd%N2_int(isd:ied,jsd:jed,nz+1), source=0.0)
1514	0	if (CS%id_Simmons_coeff > 0) then
1515	0	if (CS%CVMix_tidal_scheme /= SIMMONS) then
1516		call MOM_error(FATAL, "setup_tidal_diagnostics: Simmons_coeff diagnostics is available "//&
1517	0	"only when CVMix_tidal_scheme is Simmons")
1518		endif
1519	0	allocate(CS%dd%Simmons_coeff_2d(isd:ied,jsd:jed), source=0.0)
1520		endif
1521	0	if (CS%id_vert_dep > 0) allocate(CS%dd%vert_dep_3d(isd:ied,jsd:jed,nz+1), source=0.0)
1522	0	if (CS%id_Schmittner_coeff > 0) then
1523	0	if (CS%CVMix_tidal_scheme /= SCHMITTNER) then
1524		call MOM_error(FATAL, "setup_tidal_diagnostics: Schmittner_coeff diagnostics is available "//&
1525	0	"only when CVMix_tidal_scheme is Schmittner.")
1526		endif
1527	0	allocate(CS%dd%Schmittner_coeff_3d(isd:ied,jsd:jed,nz), source=0.0)
1528		endif
1529	0	if (CS%id_tidal_qe_md > 0) then
1530	0	if (CS%CVMix_tidal_scheme /= SCHMITTNER) then
1531		call MOM_error(FATAL, "setup_tidal_diagnostics: tidal_qe_md diagnostics is available "//&
1532	0	"only when CVMix_tidal_scheme is Schmittner.")
1533		endif
1534	0	allocate(CS%dd%tidal_qe_md(isd:ied,jsd:jed,nz), source=0.0)
1535		endif
1536	0	end subroutine setup_tidal_diagnostics
1537
1538		!> This subroutine offers up diagnostics of the tidal mixing.
1539	0	subroutine post_tidal_diagnostics(G, GV, h ,CS)
1540		type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure
1541		type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure.
1542		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
1543		intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2].
1544		type(tidal_mixing_cs), intent(inout) :: CS !< The control structure for this module
1545
1546	0	if (CS%Int_tide_dissipation .or. CS%Lee_wave_dissipation .or. CS%Lowmode_itidal_dissipation) then
1547	0	if (CS%id_TKE_itidal > 0) call post_data(CS%id_TKE_itidal, CS%dd%TKE_itidal_used, CS%diag)
1548	0	if (CS%id_TKE_leewave > 0) call post_data(CS%id_TKE_leewave, CS%TKE_Niku, CS%diag)
1549	0	if (CS%id_Nb > 0) call post_data(CS%id_Nb, CS%Nb, CS%diag)
1550	0	if (CS%id_N2_bot > 0) call post_data(CS%id_N2_bot, CS%dd%N2_bot, CS%diag)
1551	0	if (CS%id_N2_meanz > 0) call post_data(CS%id_N2_meanz,CS%dd%N2_meanz,CS%diag)
1552
1553	0	if (CS%id_Fl_itidal > 0) call post_data(CS%id_Fl_itidal, CS%dd%Fl_itidal, CS%diag)
1554	0	if (CS%id_Kd_itidal > 0) call post_data(CS%id_Kd_itidal, CS%dd%Kd_itidal, CS%diag)
1555	0	if (CS%id_Kd_Niku > 0) call post_data(CS%id_Kd_Niku, CS%dd%Kd_Niku, CS%diag)
1556	0	if (CS%id_Kd_lowmode> 0) call post_data(CS%id_Kd_lowmode, CS%dd%Kd_lowmode, CS%diag)
1557	0	if (CS%id_Fl_lowmode> 0) call post_data(CS%id_Fl_lowmode, CS%dd%Fl_lowmode, CS%diag)
1558
1559	0	if (CS%id_N2_int> 0) call post_data(CS%id_N2_int, CS%dd%N2_int, CS%diag)
1560	0	if (CS%id_vert_dep> 0) call post_data(CS%id_vert_dep, CS%dd%vert_dep_3d, CS%diag)
1561	0	if (CS%id_Simmons_coeff> 0) call post_data(CS%id_Simmons_coeff, CS%dd%Simmons_coeff_2d, CS%diag)
1562	0	if (CS%id_Schmittner_coeff> 0) call post_data(CS%id_Schmittner_coeff, CS%dd%Schmittner_coeff_3d, CS%diag)
1563	0	if (CS%id_tidal_qe_md> 0) call post_data(CS%id_tidal_qe_md, CS%dd%tidal_qe_md, CS%diag)
1564
1565	0	if (CS%id_Kd_Itidal_Work > 0) &
1566	0	call post_data(CS%id_Kd_Itidal_Work, CS%dd%Kd_Itidal_Work, CS%diag)
1567	0	if (CS%id_Kd_Niku_Work > 0) call post_data(CS%id_Kd_Niku_Work, CS%dd%Kd_Niku_Work, CS%diag)
1568	0	if (CS%id_Kd_Lowmode_Work > 0) &
1569	0	call post_data(CS%id_Kd_Lowmode_Work, CS%dd%Kd_Lowmode_Work, CS%diag)
1570
1571	0	if (CS%id_Polzin_decay_scale > 0 ) &
1572	0	call post_data(CS%id_Polzin_decay_scale, CS%dd%Polzin_decay_scale, CS%diag)
1573	0	if (CS%id_Polzin_decay_scale_scaled > 0 ) &
1574	0	call post_data(CS%id_Polzin_decay_scale_scaled, CS%dd%Polzin_decay_scale_scaled, CS%diag)
1575		endif
1576
1577	0	if (allocated(CS%dd%Kd_itidal)) deallocate(CS%dd%Kd_itidal)
1578	0	if (allocated(CS%dd%Kd_lowmode)) deallocate(CS%dd%Kd_lowmode)
1579	0	if (allocated(CS%dd%Fl_itidal)) deallocate(CS%dd%Fl_itidal)
1580	0	if (allocated(CS%dd%Fl_lowmode)) deallocate(CS%dd%Fl_lowmode)
1581	0	if (allocated(CS%dd%Polzin_decay_scale)) deallocate(CS%dd%Polzin_decay_scale)
1582	0	if (allocated(CS%dd%Polzin_decay_scale_scaled)) deallocate(CS%dd%Polzin_decay_scale_scaled)
1583	0	if (allocated(CS%dd%N2_bot)) deallocate(CS%dd%N2_bot)
1584	0	if (allocated(CS%dd%N2_meanz)) deallocate(CS%dd%N2_meanz)
1585	0	if (allocated(CS%dd%Kd_Niku)) deallocate(CS%dd%Kd_Niku)
1586	0	if (allocated(CS%dd%Kd_Niku_work)) deallocate(CS%dd%Kd_Niku_work)
1587	0	if (allocated(CS%dd%Kd_Itidal_Work)) deallocate(CS%dd%Kd_Itidal_Work)
1588	0	if (allocated(CS%dd%Kd_Lowmode_Work)) deallocate(CS%dd%Kd_Lowmode_Work)
1589	0	if (allocated(CS%dd%TKE_itidal_used)) deallocate(CS%dd%TKE_itidal_used)
1590	0	if (allocated(CS%dd%N2_int)) deallocate(CS%dd%N2_int)
1591	0	if (allocated(CS%dd%vert_dep_3d)) deallocate(CS%dd%vert_dep_3d)
1592	0	if (allocated(CS%dd%Simmons_coeff_2d)) deallocate(CS%dd%Simmons_coeff_2d)
1593	0	if (allocated(CS%dd%Schmittner_coeff_3d)) deallocate(CS%dd%Schmittner_coeff_3d)
1594	0	if (allocated(CS%dd%tidal_qe_md)) deallocate(CS%dd%tidal_qe_md)
1595	0	end subroutine post_tidal_diagnostics
1596
1597		!> This subroutine returns a zonal slice of the topographic roughness amplitudes
1598	0	subroutine tidal_mixing_h_amp(h_amp, G, j, CS)
1599		type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure
1600		real, dimension(SZI_(G)), intent(out) :: h_amp !< The topographic roughness amplitude [Z ~> m]
1601		integer, intent(in) :: j !< j-index of the row to work on
1602		type(tidal_mixing_cs), intent(in) :: CS !< The control structure for this module
1603
1604		integer :: i
1605
1606	0	h_amp(:) = 0.0
1607	0	if ( CS%Int_tide_dissipation .and. .not. CS%use_CVMix_tidal ) then
1608	0	do i=G%isc,G%iec
1609	0	h_amp(i) = sqrt(CS%h2(i,j))
1610		enddo
1611		endif
1612
1613	0	end subroutine tidal_mixing_h_amp
1614
1615		! TODO: move this subroutine to MOM_internal_tide_input module (?)
1616		!> This subroutine read tidal energy inputs from a file.
1617	0	subroutine read_tidal_energy(G, GV, US, tidal_energy_type, param_file, CS)
1618		type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure
1619		type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure
1620		type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
1621		character(len=20), intent(in) :: tidal_energy_type !< The type of tidal energy inputs to read
1622		type(param_file_type), intent(in) :: param_file !< Run-time parameter file handle
1623		type(tidal_mixing_cs), intent(inout) :: CS !< The control structure for this module
1624
1625		! local variables
1626		character(len=200) :: tidal_energy_file ! Input file names or paths
1627		character(len=200) :: tidal_input_var ! Input file variable name
1628		character(len=40) :: mdl = "MOM_tidal_mixing" !< This module's name.
1629		integer :: i, j, isd, ied, jsd, jed
1630		real, allocatable, dimension(:,:) :: &
1631	0	tidal_energy_flux_2d ! Input tidal energy flux at T-grid points [R Z3 T-3 ~> W m-2]
1632
1633	0	isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed
1634
1635		call get_param(param_file, mdl, "TIDAL_ENERGY_FILE", tidal_energy_file, &
1636		"The path to the file containing tidal energy dissipation. "//&
1637	0	"Used with CVMix tidal mixing schemes.", fail_if_missing=.true.)
1638	0	tidal_energy_file = trim(CS%inputdir) // trim(tidal_energy_file)
1639
1640	0	select case (uppercase(tidal_energy_type(1:4)))
1641		case ('JAYN') ! Jayne 2009
1642	0	if (.not. allocated(CS%tidal_qe_2d)) allocate(CS%tidal_qe_2d(isd:ied,jsd:jed))
1643	0	allocate(tidal_energy_flux_2d(isd:ied,jsd:jed))
1644		call get_param(param_file, mdl, "TIDAL_DISSIPATION_VAR", tidal_input_var, &
1645		"The name in the input file of the tidal energy source for mixing.", &
1646	0	default="wave_dissipation")
1647	0	call MOM_read_data(tidal_energy_file, tidal_input_var, tidal_energy_flux_2d, G%domain, scale=US%W_m2_to_RZ3_T3)
1648	0	do j=G%jsc,G%jec ; do i=G%isc,G%iec
1649	0	CS%tidal_qe_2d(i,j) = CS%Gamma_itides * tidal_energy_flux_2d(i,j)
1650		enddo ; enddo
1651	0	deallocate(tidal_energy_flux_2d)
1652		case ('ER03') ! Egbert & Ray 2003
1653	0	call read_tidal_constituents(G, GV, US, tidal_energy_file, param_file, CS)
1654		case default
1655	0	call MOM_error(FATAL, "read_tidal_energy: Unknown tidal energy file type.")
1656		end select
1657
1658	0	end subroutine read_tidal_energy
1659
1660		!> This subroutine reads tidal input energy from a file by constituent.
1661	0	subroutine read_tidal_constituents(G, GV, US, tidal_energy_file, param_file, CS)
1662		type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure
1663		type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure
1664		type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
1665		character(len=200), intent(in) :: tidal_energy_file !< The file from which to read tidal energy inputs
1666		type(param_file_type), intent(in) :: param_file !< Run-time parameter file handle
1667		type(tidal_mixing_cs), intent(inout) :: CS !< The control structure for this module
1668
1669		! local variables
1670		real, parameter :: C1_3 = 1.0/3.0 ! A rational constant [nondim]
1671		real, dimension(SZI_(G),SZJ_(G)) :: &
1672	0	tidal_qk1, & ! qk1 coefficient used in Schmittner & Egbert [nondim]
1673	0	tidal_qo1 ! qo1 coefficient used in Schmittner & Egbert [nondim]
1674		real, allocatable, dimension(:) :: &
1675	0	z_t, & ! depth from surface to midpoint of input layer [Z ~> m]
1676	0	z_w ! depth from surface to top of input layer [Z ~> m]
1677		real, allocatable, dimension(:,:,:) :: &
1678	0	tc_m2, & ! input lunar semidiurnal tidal energy flux [R Z3 T-3 ~> W m-2]
1679	0	tc_s2, & ! input solar semidiurnal tidal energy flux [R Z3 T-3 ~> W m-2]
1680	0	tc_k1, & ! input lunar diurnal tidal energy flux [R Z3 T-3 ~> W m-2]
1681	0	tc_o1 ! input lunar diurnal tidal energy flux [R Z3 T-3 ~> W m-2]
1682		integer, dimension(4) :: nz_in
1683		integer :: k, is, ie, js, je, isd, ied, jsd, jed, i, j
1684
1685	0	is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec
1686	0	isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed
1687
1688		! get number of input levels:
1689	0	call field_size(tidal_energy_file, 'z_t', nz_in)
1690
1691		! allocate local variables
1692	0	allocate(z_t(nz_in(1)), z_w(nz_in(1)) )
1693	0	allocate(tc_m2(isd:ied,jsd:jed,nz_in(1)), &
1694	0	tc_s2(isd:ied,jsd:jed,nz_in(1)), &
1695	0	tc_k1(isd:ied,jsd:jed,nz_in(1)), &
1696	0	tc_o1(isd:ied,jsd:jed,nz_in(1)) )
1697
1698		! allocate CS variables associated with 3d tidal energy dissipation
1699	0	if (.not. allocated(CS%tidal_qe_3d_in)) allocate(CS%tidal_qe_3d_in(isd:ied,jsd:jed,nz_in(1)))
1700	0	if (.not. allocated(CS%h_src)) allocate(CS%h_src(nz_in(1)))
1701
1702		! read in tidal constituents
1703	0	call MOM_read_data(tidal_energy_file, 'M2', tc_m2, G%domain, scale=US%W_m2_to_RZ3_T3)
1704	0	call MOM_read_data(tidal_energy_file, 'S2', tc_s2, G%domain, scale=US%W_m2_to_RZ3_T3)
1705	0	call MOM_read_data(tidal_energy_file, 'K1', tc_k1, G%domain, scale=US%W_m2_to_RZ3_T3)
1706	0	call MOM_read_data(tidal_energy_file, 'O1', tc_o1, G%domain, scale=US%W_m2_to_RZ3_T3)
1707		! Note the hard-coded assumption that z_t and z_w in the file are in centimeters.
1708	0	call MOM_read_data(tidal_energy_file, 'z_t', z_t, scale=0.01*US%m_to_Z)
1709	0	call MOM_read_data(tidal_energy_file, 'z_w', z_w, scale=0.01*US%m_to_Z)
1710
1711	0	do j=js,je ; do i=is,ie
1712	0	if (abs(G%geoLatT(i,j)) < 30.0) then
1713	0	tidal_qk1(i,j) = C1_3
1714	0	tidal_qo1(i,j) = C1_3
1715		else
1716	0	tidal_qk1(i,j) = 1.0
1717	0	tidal_qo1(i,j) = 1.0
1718		endif
1719		enddo ; enddo
1720
1721	0	CS%tidal_qe_3d_in(:,:,:) = 0.0
1722	0	do k=1,nz_in(1)
1723		! Store the input cell thickness in m for use with CVmix.
1724	0	CS%h_src(k) = US%Z_to_m(z_t(k)-z_w(k))2.0
1725		! form tidal_qe_3d_in from weighted tidal constituents
1726	0	do j=js,je ; do i=is,ie
1727	0	if ((z_t(k) <= G%bathyT(i,j) + G%Z_ref) .and. (z_w(k) > CS%tidal_diss_lim_tc)) &
1728		CS%tidal_qe_3d_in(i,j,k) = C1_3tc_m2(i,j,k) + C1_3tc_s2(i,j,k) + &
1729	0	tidal_qk1(i,j)tc_k1(i,j,k) + tidal_qo1(i,j)tc_o1(i,j,k)
1730		enddo ; enddo
1731		enddo
1732
1733		! test if qE is positive
1734	0	if (any(CS%tidal_qe_3d_in<0.0)) then
1735	0	call MOM_error(FATAL, "read_tidal_constituents: Negative tidal_qe_3d_in terms.")
1736		endif
1737
1738		!! collapse 3D qE to 2D qE
1739		!CS%tidal_qe_2d(:,:) = 0.0
1740		!do k=1,nz_in(1) ; do j=js,je ; do i=is,ie
1741		! if (z_t(k) <= G%bathyT(i,j) + G%Z_ref) &
1742		! CS%tidal_qe_2d(i,j) = CS%tidal_qe_2d(i,j) + CS%tidal_qe_3d_in(i,j,k)
1743		!enddo ; enddo ; enddo
1744
1745		! initialize input remapping:
1746		call initialize_remapping(CS%remap_cs, remapping_scheme="PLM", &
1747		boundary_extrapolation=.false., check_remapping=CS%debug, &
1748		answer_date=CS%remap_answer_date, &
1749	0	h_neglect=GV%H_subroundoff, h_neglect_edge=GV%H_subroundoff)
1750
1751	0	deallocate(tc_m2)
1752	0	deallocate(tc_s2)
1753	0	deallocate(tc_k1)
1754	0	deallocate(tc_o1)
1755	0	deallocate(z_t)
1756	0	deallocate(z_w)
1757
1758	0	end subroutine read_tidal_constituents
1759
1760		!> Deallocate fields
1761	0	subroutine tidal_mixing_end(CS)
1762		type(tidal_mixing_cs), intent(inout) :: CS !< This module's control structure, which
1763		!! will be deallocated in this routine.
1764
1765		! TODO: deallocate all the dynamically allocated members here ...
1766	0	if (allocated(CS%tidal_qe_2d)) deallocate(CS%tidal_qe_2d)
1767	0	if (allocated(CS%tidal_qe_3d_in)) deallocate(CS%tidal_qe_3d_in)
1768	0	if (allocated(CS%h_src)) deallocate(CS%h_src)
1769	0	end subroutine tidal_mixing_end
1770
1771	0	end module MOM_tidal_mixing