MOM_tidal_mixing.F90

1! This file is part of MOM6, the Modular Ocean Model version 6.
2! See the LICENSE file for licensing information.
3! SPDX-License-Identifier: Apache-2.0
4
5!> Interface to vertical tidal mixing schemes including CVMix tidal mixing.
7
9use mom_diag_mediator, only : safe_alloc_ptr, post_data
11use mom_debugging, only : hchksum
12use mom_error_handler, only : mom_error, is_root_pe, fatal, warning, note
14use mom_file_parser, only : get_param, log_param, log_version, param_file_type
15use mom_grid, only : ocean_grid_type
16use mom_io, only : slasher, mom_read_data, field_size
17use mom_io, only : read_netcdf_data
24use cvmix_tidal, only : cvmix_init_tidal, cvmix_compute_simmons_invariant
25use cvmix_tidal, only : cvmix_coeffs_tidal, cvmix_tidal_params_type
26use cvmix_tidal, only : cvmix_compute_schmittner_invariant, cvmix_compute_schmittnercoeff
27use cvmix_tidal, only : cvmix_coeffs_tidal_schmittner
28use cvmix_kinds_and_types, only : cvmix_global_params_type
29use cvmix_put_get, only : cvmix_put
30
31implicit none ; private
32
33#include <MOM_memory.h>
34
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
48type, 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]
72end type
73
74!> Control structure with parameters for the tidal mixing module.
75type, 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
205end type tidal_mixing_cs
206
207!>@{ Coded parmameters for specifying mixing schemes
208character*(20), parameter :: stlaurent_profile_string = "STLAURENT_02"
209character*(20), parameter :: polzin_profile_string = "POLZIN_09"
210integer, parameter :: stlaurent_02 = 1
211integer, parameter :: polzin_09 = 2
212character*(20), parameter :: simmons_scheme_string = "SIMMONS"
213character*(20), parameter :: schmittner_scheme_string = "SCHMITTNER"
214integer, parameter :: simmons = 1
215integer, parameter :: schmittner = 2
216!>@}
217
218contains
219
220!> Initializes internal tidal dissipation scheme for diapycnal mixing
221logical 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 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
255 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 default=.false., do_not_log=.true.)
262 call get_param(param_file, mdl, "INT_TIDE_DISSIPATION", int_tide_dissipation, &
263 default=use_cvmix_tidal, do_not_log=.true.)
264 call log_version(param_file, mdl, version, &
265 "Vertical Tidal Mixing Parameterization", &
266 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 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 "et al. (2002) and Simmons et al. (2004).", default=use_cvmix_tidal)
275
276 ! return if tidal mixing is inactive
277 tidal_mixing_init = int_tide_dissipation
278 if (.not. tidal_mixing_init) return
279
280 cs%debug = cs%debug.and.is_root_pe()
281 cs%diag => diag
282 if (associated(int_tide_csp)) cs%int_tide_CSp => int_tide_csp
283 cs%use_CVmix_tidal = use_cvmix_tidal
284 cs%int_tide_dissipation = int_tide_dissipation
285
286 call get_param(param_file, mdl, "INPUTDIR", cs%inputdir, default=".",do_not_log=.true.)
287 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 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 default=default_answer_date, do_not_log=.not.gv%Boussinesq)
297 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 default=default_answer_date, do_not_log=.not.gv%Boussinesq)
305 if (.not.gv%Boussinesq) cs%remap_answer_date = max(cs%remap_answer_date, 20230701)
306
307 if (cs%int_tide_dissipation) then
308
309 ! Read in CVMix tidal scheme if CVMix tidal mixing is on
310 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 default=simmons_scheme_string)
319 cvmix_tidal_scheme_str = uppercase(cvmix_tidal_scheme_str)
320
321 select case (cvmix_tidal_scheme_str)
322 case (simmons_scheme_string) ; cs%CVMix_tidal_scheme = simmons
323 case (schmittner_scheme_string) ; cs%CVMix_tidal_scheme = schmittner
324 case default
325 call mom_error(fatal, "tidal_mixing_init: Unrecognized setting "// &
326 "#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 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 default=stlaurent_profile_string)
340 int_tide_profile_str = uppercase(int_tide_profile_str)
341
342 select case (int_tide_profile_str)
343 case (stlaurent_profile_string) ; cs%int_tide_profile = stlaurent_02
344 case (polzin_profile_string) ; cs%int_tide_profile = polzin_09
345 case default
346 call mom_error(fatal, "tidal_mixing_init: Unrecognized setting "// &
347 "#define INT_TIDE_PROFILE "//trim(int_tide_profile_str)//" found in input file.")
348 end select
349 endif
350
351 elseif (cs%use_CVMix_tidal) then
352 call mom_error(fatal, "tidal_mixing_init: Cannot set INT_TIDE_DISSIPATION to False "// &
353 "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 "and Simmons et al. (2004) vertical profile", default=.false.)
361 if (cs%lee_wave_dissipation) then
362 if (cs%use_CVMix_tidal) then
363 call mom_error(fatal, "tidal_mixing_init: Lee wave driven dissipation scheme cannot "// &
364 "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 default=stlaurent_profile_string)
374 tmpstr = uppercase(tmpstr)
375 select case (tmpstr)
376 case (stlaurent_profile_string) ; cs%lee_wave_profile = stlaurent_02
377 case (polzin_profile_string) ; cs%lee_wave_profile = polzin_09
378 case default
379 call mom_error(fatal, "tidal_mixing_init: Unrecognized setting "// &
380 "#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 "et al. (2002) and Simmons et al. (2004).", default=.false.)
390
391 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 if (cs%use_CVMix_tidal) then
394 call mom_error(fatal, "tidal_mixing_init: Polzin scheme cannot "// &
395 "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 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 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 "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 "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 units="m", default=0.0, scale=us%m_to_Z)
424 endif
425
426 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 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 "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 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 "ocean depth is less than this value.", units="m", default=0.0, scale=us%m_to_Z)
442 endif
443
444 if ( (cs%Int_tide_dissipation .or. cs%Lee_wave_dissipation) .and. &
445 .not. cs%use_CVMix_tidal) then
446
447 allocate(cs%Nb(isd:ied,jsd:jed), source=0.)
448 allocate(cs%h2(isd:ied,jsd:jed), source=0.)
449 allocate(cs%TKE_itidal(isd:ied,jsd:jed), source=0.)
450 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 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 units="m s-1", default=0.0, scale=us%m_to_Z*us%T_to_s)
460 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 "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 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 "the tidal amplitude with INT_TIDE_DISSIPATION.", default=.false.)
473 if (read_tideamp) then
474 if (cs%use_CVMix_tidal) then
475 call mom_error(fatal, "tidal_mixing_init: Tidal amplitude files are "// &
476 "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 "tidal amplitudes with INT_TIDE_DISSIPATION.", default="tideamp.nc")
481 filename = trim(cs%inputdir) // trim(tideamp_file)
482 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 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 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 fail_if_missing=(.not.cs%use_CVMix_tidal))
496 filename = trim(cs%inputdir) // trim(h2_file)
497 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 "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 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 units="nondim", default=0.1)
510
511 do j=js,je ; do i=is,ie
512 if (max(g%meanSL(i,j) + g%bathyT(i,j), 0.0) < cs%min_zbot_itides) cs%mask_itidal(i,j) = 0.0
513 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 if ((cs%tidal_answer_date < 20190101) .and. (max_frac_rough >= 0.0)) then
517 hamp = min(max_frac_rough * max(g%meanSL(i,j) + g%bathyT(i,j), 0.0), sqrt(cs%h2(i,j)))
518 cs%h2(i,j) = hamp*hamp
519 else
520 if (max_frac_rough >= 0.0) &
521 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 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 cs%kappa_itides * cs%h2(i,j) * utide*utide
529 enddo ; enddo
530
531 endif
532
533 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 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 units="nondim", default=1.0)
543
544 filename = trim(cs%inputdir) // trim(niku_tke_input_file)
545 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 default="TKE_input")
549 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 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 "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 "dissipation of lee wave dissipation.", units="nondim", default=1.0)
560 else
561 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 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 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 units="m", default=0.0, scale=us%m_to_Z)
576 call get_param(param_file, mdl, 'MIN_THICKNESS', cs%min_thickness, &
577 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 units="nondim", default=1.0, do_not_log=.true.)
582 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 fail_if_missing=.true.)
589 ! Check whether tidal energy input format and CVMix tidal mixing scheme are consistent
590 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 " mixing scheme: "//trim(cvmix_tidal_scheme_str) )
596 endif
597 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 depth_cutoff = cs%min_zbot_itides*us%Z_to_m)
607
608 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 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 'Internal Tide Driven Diffusivity', 'm2 s-1', conversion=gv%HZ_T_to_m2_s)
621
622 if (cs%use_CVMix_tidal) then
623 cs%id_N2_int = register_diag_field('ocean_model','N2_int',diag%axesTi,time, &
624 'Bouyancy frequency squared, at interfaces', 's-2', conversion=us%s_to_T**2)
625 !> TODO: add units
626 if (cs%CVMix_tidal_scheme .eq. simmons) then
627 cs%id_Simmons_coeff = register_diag_field('ocean_model','Simmons_coeff',diag%axesT1,time, &
628 'time-invariant portion of the tidal mixing coefficient using the Simmons', '')
629 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 '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 '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 '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 'W m-2', conversion=us%RZ3_T3_to_W_m2)
642 cs%id_Nb = register_diag_field('ocean_model','Nb',diag%axesT1,time, &
643 '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 '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 'm3 s-3', conversion=(gv%H_to_m*us%Z_to_m**2*us%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 'm3 s-3', conversion=(gv%H_to_m*us%Z_to_m**2*us%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 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 '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 '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 '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 '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 '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 'W m-2', conversion=us%RZ3_T3_to_W_m2)
683
684 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 '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 '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
694end 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.
700subroutine calculate_tidal_mixing(dz, j, N2_bot, Rho_bot, N2_lay, N2_int, TKE_to_Kd, max_TKE, &
701 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 if (cs%Int_tide_dissipation .or. cs%Lee_wave_dissipation .or. cs%Lowmode_itidal_dissipation) then
738 if (cs%use_CVMix_tidal) then
739 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 g, gv, us, cs, kd_max, kd_lay, kd_int, vbf)
743 endif
744 endif
745end 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.
750subroutine 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 real, dimension(SZK_(GV)+1) :: Kd_tidal ! tidal diffusivity [m2 s-1]
769 real, dimension(SZK_(GV)+1) :: Kv_tidal ! tidal viscosity [m2 s-1]
770 real, dimension(SZK_(GV)+1) :: vert_dep ! vertical deposition [nondim]
771 real, dimension(SZK_(GV)+1) :: iFaceHeight ! Height of interfaces [m]
772 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 real, dimension(SZK_(GV)) :: cellHeight ! Height of cell centers [m]
776 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 real, dimension(SZK_(GV)+1) :: N2_int_i ! De-scaled interface buoyancy frequency [s-2]
779 real, dimension(SZK_(GV)) :: Schmittner_coeff ! A coefficient in the Schmittner et al (2014) mixing
780 ! parameterization [nondim]
781 real, dimension(SZK_(GV)) :: h_m ! Cell thickness [m]
782 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 is = g%isc ; ie = g%iec
794
795 select case (cs%CVMix_tidal_scheme)
796 case (simmons)
797 do i=is,ie
798
799 if (g%mask2dT(i,j)<1) cycle
800
801 ifaceheight = 0.0 ! BBL is all relative to the surface
802 hcorr = 0.0
803 ! Compute cell center depth and cell bottom in meters (negative values in the ocean)
804 do k=1,gv%ke
805 dh = dz(i,k) ! Nominal thickness to use for increment, in the units of heights
806 dh = dh + hcorr ! Take away the accumulated error (could temporarily make dh<0)
807 hcorr = min( dh - cs%min_thickness, 0. ) ! If inflating then hcorr<0
808 dh = max(dh, cs%min_thickness) ! Limited increment dh>=min_thickness
809 cellheight(k) = ifaceheight(k) - 0.5 * us%Z_to_m*dh
810 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 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 simmons_coeff = simmons_coeff / cs%Gamma_itides
826
827
828 ! XXX: Temporary de-scaling of N2_int(i,:) into a temporary variable
829 do k=1,gv%ke+1
830 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 cvmix_tidal_params_user = cs%CVMix_tidal_params)
843
844 ! Update diffusivity
845 if (present(kd_lay)) then
846 do k=1,gv%ke
847 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 if (present(kd_int)) then
851 do k=1,gv%ke+1
852 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 if (associated(kv)) then
857 do k=1,gv%ke+1
858 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 if (allocated(cs%dd%Kd_itidal)) then
864 cs%dd%Kd_itidal(i,j,:) = gv%m2_s_to_HZ_T * kd_tidal(:)
865 endif
866 if (allocated(cs%dd%N2_int)) then
867 cs%dd%N2_int(i,j,:) = n2_int(i,:)
868 endif
869 if (allocated(cs%dd%Simmons_coeff_2d)) then
870 cs%dd%Simmons_coeff_2d(i,j) = simmons_coeff
871 endif
872 if (allocated(cs%dd%vert_dep_3d)) then
873 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 allocate(exp_hab_zetar(gv%ke+1,gv%ke+1))
884
885 do i=is,ie
886
887 if (g%mask2dT(i,j)<1) cycle
888
889 ifaceheight(:) = 0.0 ! BBL is all relative to the surface
890 hcorr = 0.0
891 ! Compute heights at cell center and interfaces, and rescale layer thicknesses
892 do k=1,gv%ke
893 h_m(k) = dz(i,k)*us%Z_to_m ! Rescale thicknesses to m for use by CVmix.
894 dh = dz(i,k) ! Nominal thickness to use for increment, in the units of heights
895 dh = dh + hcorr ! Take away the accumulated error (could temporarily make dh<0)
896 hcorr = min( dh - cs%min_thickness, 0. ) ! If inflating then hcorr<0
897 dh = max(dh, cs%min_thickness) ! Limited increment dh>=min_thickness
898 cellheight(k) = ifaceheight(k) - 0.5 * us%Z_to_m*dh
899 ifaceheight(k+1) = ifaceheight(k) - us%Z_to_m*dh
900 enddo
901
902 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 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 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 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 cvmix_tidal_params_user = cs%CVMix_tidal_params)
928
929 ! XXX: Temporary de-scaling of N2_int(i,:) into a temporary variable
930 do k=1,gv%ke+1
931 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 cvmix_tidal_params_user = cs%CVMix_tidal_params)
944
945 ! Update diffusivity
946 if (present(kd_lay)) then
947 do k=1,gv%ke
948 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 if (present(kd_int)) then
952 do k=1,gv%ke+1
953 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 if (associated(kv)) then
959 do k=1,gv%ke+1
960 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 if (allocated(cs%dd%Kd_itidal)) then
966 cs%dd%Kd_itidal(i,j,:) = gv%m2_s_to_HZ_T*kd_tidal(:)
967 endif
968 if (allocated(cs%dd%N2_int)) then
969 cs%dd%N2_int(i,j,:) = n2_int(i,:)
970 endif
971 if (allocated(cs%dd%Schmittner_coeff_3d)) then
972 cs%dd%Schmittner_coeff_3d(i,j,:) = schmittner_coeff(:)
973 endif
974 if (allocated(cs%dd%tidal_qe_md)) then
975 cs%dd%tidal_qe_md(i,j,:) = tidal_qe_md(:)
976 endif
977 if (allocated(cs%dd%vert_dep_3d)) then
978 cs%dd%vert_dep_3d(i,j,:) = vert_dep(:)
979 endif
980 enddo ! i=is,ie
981
982 deallocate(exp_hab_zetar)
983
984 case default
985 call mom_error(fatal, "tidal_mixing_init: Unrecognized setting "// &
986 "#define CVMIX_TIDAL_SCHEME found in input file.")
987 end select
988
989end 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.).
998subroutine add_int_tide_diffusivity(dz, j, N2_bot, Rho_bot, N2_lay, TKE_to_Kd, max_TKE, &
999 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 dztot, & ! Vertical distance between the top and bottom of the ocean [Z ~> m]
1035 dztot_WKB, & ! WKB scaled distance from top to bottom [Z ~> m]
1036 TKE_itidal_bot, & ! internal tide TKE at ocean bottom [H Z2 T-3 ~> m3 s-3 or W m-2]
1037 TKE_Niku_bot, & ! lee-wave TKE at ocean bottom [H Z2 T-3 ~> m3 s-3 or W m-2]
1038 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 Inv_int, & ! inverse of TKE decay for int tide over the depth of the ocean [nondim]
1040 Inv_int_lee, & ! inverse of TKE decay for lee waves over the depth of the ocean [nondim]
1041 Inv_int_low, & ! inverse of TKE decay for low modes over the depth of the ocean [nondim]
1042 z0_Polzin, & ! TKE decay scale in Polzin formulation [Z ~> m]
1043 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 n2_meanz, & ! vertically averaged squared buoyancy frequency [T-2 ~> s-2] for WKB scaling
1048 tke_itidal_rem, & ! remaining internal tide TKE (from barotropic source) [H Z2 T-3 ~> m3 s-3 or W m-2]
1049 tke_niku_rem, & ! remaining lee-wave TKE [H Z2 T-3 ~> m3 s-3 or W m-2]
1050 tke_lowmode_rem, & ! remaining internal tide TKE (from propagating low mode source) [H Z2 T-3 ~> m3 s-3 or W m-2]
1051 tke_frac_top, & ! fraction of bottom TKE that should appear at top of a layer [nondim]
1052 tke_frac_top_lee, & ! fraction of bottom TKE that should appear at top of a layer [nondim]
1053 tke_frac_top_lowmode, &
1054 ! fraction of bottom TKE that should appear at top of a layer [nondim]
1055 z_from_bot, & ! distance from bottom [Z ~> m]
1056 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 is = g%isc ; ie = g%iec ; nz = gv%ke
1074
1075 if (.not.(cs%Int_tide_dissipation .or. cs%Lee_wave_dissipation)) return
1076
1077 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 do k=1,nz ; do i=is,ie
1079 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 (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 (cs%Lowmode_itidal_dissipation .and. (cs%int_tide_profile == stlaurent_02)))
1088
1089 ! Calculate parameters for vertical structure of dissipation
1090 ! Simmons:
1091 if ( use_simmons ) then
1092 izeta = 1.0 / max(cs%Int_tide_decay_scale, gv%dz_subroundoff)
1093 izeta_lee = 1.0 / max(cs%Int_tide_decay_scale*cs%Decay_scale_factor_lee, gv%dz_subroundoff)
1094 do i=is,ie
1095 cs%Nb(i,j) = sqrt(n2_bot(i))
1096 if (allocated(cs%dd%N2_bot)) &
1097 cs%dd%N2_bot(i,j) = n2_bot(i)
1098 if ( cs%Int_tide_dissipation ) then
1099 if (izeta*dztot(i) > 1.0e-14) then ! L'Hospital's version of Adcroft's reciprocal rule.
1100 inv_int(i) = 1.0 / (1.0 - exp(-izeta*dztot(i)))
1101 endif
1102 endif
1103 if ( cs%Lee_wave_dissipation ) then
1104 if (izeta_lee*dztot(i) > 1.0e-14) then ! L'Hospital's version of Adcroft's reciprocal rule.
1105 inv_int_lee(i) = 1.0 / (1.0 - exp(-izeta_lee*dztot(i)))
1106 endif
1107 endif
1108 if ( cs%Lowmode_itidal_dissipation) then
1109 if (izeta*dztot(i) > 1.0e-14) then ! L'Hospital's version of Adcroft's reciprocal rule.
1110 inv_int_low(i) = 1.0 / (1.0 - exp(-izeta*dztot(i)))
1111 endif
1112 endif
1113 z_from_bot(i) = dz(i,nz)
1114 enddo
1115 endif ! Simmons
1116
1117 ! Polzin:
1118 if ( use_polzin ) then
1119 ! WKB scaling of the vertical coordinate
1120 do i=is,ie ; n2_meanz(i) = 0.0 ; enddo
1121 do k=1,nz ; do i=is,ie
1122 n2_meanz(i) = n2_meanz(i) + n2_lay(i,k) * dz(i,k)
1123 enddo ; enddo
1124 do i=is,ie
1125 n2_meanz(i) = n2_meanz(i) / (dztot(i) + gv%dz_subroundoff)
1126 if (allocated(cs%dd%N2_meanz)) &
1127 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 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 do i=is,ie
1139 cs%Nb(i,j) = sqrt(n2_bot(i))
1140 if (cs%tidal_answer_date < 20190101) then
1141 if ((cs%tideamp(i,j) > 0.0) .and. &
1142 (cs%kappa_itides**2 * cs%h2(i,j) * cs%Nb(i,j)**3 > 1.0e-14*us%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 ( cs%kappa_itides**2 * cs%h2(i,j) * cs%Nb(i,j)**3 )
1146 if (z0_polzin(i) < cs%Polzin_min_decay_scale) &
1147 z0_polzin(i) = cs%Polzin_min_decay_scale
1148 if (n2_meanz(i) > 1.0e-14*us%T_to_s**2 ) then
1149 z0_polzin_scaled(i) = z0_polzin(i)*cs%Nb(i,j)**2 / n2_meanz(i)
1150 else
1151 z0_polzin_scaled(i) = cs%Polzin_decay_scale_max_factor * dztot(i)
1152 endif
1153 if (z0_polzin_scaled(i) > (cs%Polzin_decay_scale_max_factor * dztot(i)) ) &
1154 z0_polzin_scaled(i) = cs%Polzin_decay_scale_max_factor * dztot(i)
1155 else
1156 z0_polzin(i) = cs%Polzin_decay_scale_max_factor * dztot(i)
1157 z0_polzin_scaled(i) = cs%Polzin_decay_scale_max_factor * dztot(i)
1158 endif
1159 else
1160 z0ps_num = (cs%Polzin_decay_scale_factor * cs%Nu_Polzin * cs%Nbotref_Polzin**2) * cs%tideamp(i,j)
1161 z0ps_denom = ( cs%kappa_itides**2 * cs%h2(i,j) * cs%Nb(i,j) * n2_meanz(i) )
1162 if ((cs%tideamp(i,j) > 0.0) .and. &
1163 (z0ps_num < z0ps_denom * cs%Polzin_decay_scale_max_factor * dztot(i))) then
1164 z0_polzin_scaled(i) = z0ps_num / z0ps_denom
1165
1166 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 z0_polzin(i) = z0_polzin_scaled(i) * (n2_meanz(i) / cs%Nb(i,j)**2)
1169 else
1170 z0_polzin(i) = cs%Polzin_decay_scale_max_factor * dztot(i)
1171 endif
1172 else
1173 z0_polzin(i) = cs%Polzin_decay_scale_max_factor * dztot(i)
1174 z0_polzin_scaled(i) = cs%Polzin_decay_scale_max_factor * dztot(i)
1175 endif
1176 endif
1177
1178 if (allocated(cs%dd%Polzin_decay_scale)) &
1179 cs%dd%Polzin_decay_scale(i,j) = z0_polzin(i)
1180 if (allocated(cs%dd%Polzin_decay_scale_scaled)) &
1181 cs%dd%Polzin_decay_scale_scaled(i,j) = z0_polzin_scaled(i)
1182 if (allocated(cs%dd%N2_bot)) &
1183 cs%dd%N2_bot(i,j) = cs%Nb(i,j)*cs%Nb(i,j)
1184
1185 if (cs%tidal_answer_date < 20190101) then
1186 ! These expressions use dimensional constants to avoid NaN values.
1187 if ( cs%Int_tide_dissipation .and. (cs%int_tide_profile == polzin_09) ) then
1188 if (dztot_wkb(i) > 1.0e-14*us%m_to_Z) &
1189 inv_int(i) = ( z0_polzin_scaled(i) / dztot_wkb(i) ) + 1.0
1190 endif
1191 if ( cs%lee_wave_dissipation .and. (cs%lee_wave_profile == polzin_09) ) then
1192 if (dztot_wkb(i) > 1.0e-14*us%m_to_Z) &
1193 inv_int_lee(i) = ( z0_polzin_scaled(i)*cs%Decay_scale_factor_lee / dztot_wkb(i) ) + 1.0
1194 endif
1195 if ( cs%Lowmode_itidal_dissipation .and. (cs%int_tide_profile == polzin_09) ) then
1196 if (dztot_wkb(i) > 1.0e-14*us%m_to_Z) &
1197 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 inv_int(i) = 0.0 ; inv_int_lee(i) = 0.0 ; inv_int_low(i) = 0.0
1202 if ( cs%Int_tide_dissipation .and. (cs%int_tide_profile == polzin_09) ) then
1203 if (z0_polzin_scaled(i) < 1.0e14 * dztot_wkb(i)) &
1204 inv_int(i) = ( z0_polzin_scaled(i) / dztot_wkb(i) ) + 1.0
1205 endif
1206 if ( cs%lee_wave_dissipation .and. (cs%lee_wave_profile == polzin_09) ) then
1207 if (z0_polzin_scaled(i) < 1.0e14 * dztot_wkb(i)) &
1208 inv_int_lee(i) = ( z0_polzin_scaled(i)*cs%Decay_scale_factor_lee / dztot_wkb(i) ) + 1.0
1209 endif
1210 if ( cs%Lowmode_itidal_dissipation .and. (cs%int_tide_profile == polzin_09) ) then
1211 if (z0_polzin_scaled(i) < 1.0e14 * dztot_wkb(i)) &
1212 inv_int_low(i) = ( z0_polzin_scaled(i) / dztot_wkb(i) ) + 1.0
1213 endif
1214 endif
1215
1216 z_from_bot(i) = dz(i,nz)
1217 ! Use the new formulation for WKB scaling. N2 is referenced to its vertical mean.
1218 if (cs%tidal_answer_date < 20190101) then
1219 if (n2_meanz(i) > 1.0e-14*us%T_to_s**2 ) then
1220 z_from_bot_wkb(i) = dz(i,nz) * n2_lay(i,nz) / n2_meanz(i)
1221 else ; z_from_bot_wkb(i) = 0 ; endif
1222 else
1223 if (dz(i,nz) * n2_lay(i,nz) < n2_meanz(i) * (1.0e14 * dztot_wkb(i))) then
1224 z_from_bot_wkb(i) = dz(i,nz) * n2_lay(i,nz) / n2_meanz(i)
1225 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 do i=is,ie
1233 ! Dissipation of locally trapped internal tide (non-propagating high modes)
1234 if (gv%Boussinesq .or. gv%semi_Boussinesq) then
1235 tke_itidal_bot(i) = min(gv%Z_to_H*cs%TKE_itidal(i,j)*cs%Nb(i,j), cs%TKE_itide_max)
1236 else
1237 tke_itidal_bot(i) = min(gv%RZ_to_H*rho_bot(i) * (cs%TKE_itidal(i,j)*cs%Nb(i,j)), &
1238 cs%TKE_itide_max)
1239 endif
1240 if (allocated(cs%dd%TKE_itidal_used)) &
1241 cs%dd%TKE_itidal_used(i,j) = tke_itidal_bot(i)
1242 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 tke_niku_bot(i) = 0.0
1245 if (cs%Lee_wave_dissipation) then
1246 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 tke_lowmode_tot = 0.0
1250 tke_lowmode_bot(i) = 0.0
1251 if (cs%Lowmode_itidal_dissipation) then
1252 ! get loss rate due to wave drag on low modes (already multiplied by q)
1253 call get_lowmode_loss(i,j,g,cs%int_tide_CSp,"WaveDrag",tke_lowmode_tot)
1254 tke_lowmode_bot(i) = cs%Mu_itides * gv%RZ_to_H * tke_lowmode_tot
1255 endif
1256 ! Vertical energy flux at bottom
1257 tke_itidal_rem(i) = inv_int(i) * tke_itidal_bot(i)
1258 tke_niku_rem(i) = inv_int_lee(i) * tke_niku_bot(i)
1259 tke_lowmode_rem(i) = inv_int_low(i) * tke_lowmode_bot(i)
1260
1261 if (allocated(cs%dd%Fl_itidal)) &
1262 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 if ( use_simmons ) then
1268 do k=nz-1,2,-1 ; do i=is,ie
1269 if (max_tke(i,k) <= 0.0) cycle
1270 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 tke_frac_top(i) = inv_int(i) * exp(-izeta * z_from_bot(i))
1274 tke_frac_top_lee(i) = inv_int_lee(i) * exp(-izeta_lee * z_from_bot(i))
1275 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 tke_itide_lay = tke_itidal_rem(i) - tke_itidal_bot(i) * tke_frac_top(i)
1280 tke_niku_lay = tke_niku_rem(i) - tke_niku_bot(i) * tke_frac_top_lee(i)
1281 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 if (tke_itide_lay + tke_niku_lay + tke_lowmode_lay > max_tke(i,k)) then
1285 frac_used = (max_tke(i,k)) / (tke_itide_lay + tke_niku_lay + tke_lowmode_lay)
1286 tke_itide_lay = frac_used * tke_itide_lay
1287 tke_niku_lay = frac_used * tke_niku_lay
1288 tke_lowmode_lay = frac_used * tke_lowmode_lay
1289 endif
1290
1291 ! Calculate vertical flux available to bottom of layer above
1292 tke_itidal_rem(i) = tke_itidal_rem(i) - tke_itide_lay
1293 tke_niku_rem(i) = tke_niku_rem(i) - tke_niku_lay
1294 tke_lowmode_rem(i) = tke_lowmode_rem(i) - tke_lowmode_lay
1295
1296 ! Convert power to diffusivity
1297 kd_add = tke_to_kd(i,k) * (tke_itide_lay + tke_niku_lay + tke_lowmode_lay)
1298
1299 if (kd_max >= 0.0) kd_add = min(kd_add, kd_max)
1300 if (present(kd_lay)) then
1301 kd_lay(i,k) = kd_lay(i,k) + kd_add
1302 endif
1303
1304 if (present(kd_int)) then
1305 kd_int(i,k) = kd_int(i,k) + 0.5 * kd_add
1306 kd_int(i,k+1) = kd_int(i,k+1) + 0.5 * kd_add
1307 endif
1308
1309 ! diagnostics
1310 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 kd_add = tke_to_kd(i,k) * tke_itide_lay
1314 if (kd_max >= 0.0) kd_add = min(kd_add, kd_max)
1315 if (allocated(cs%dd%Kd_itidal)) then
1316 if (k>1) cs%dd%Kd_itidal(i,j,k) = cs%dd%Kd_itidal(i,j,k) + 0.5*kd_add
1317 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 if (associated(vbf%Kd_itides)) then
1320 !Not to be confused w/ Kd_itidal (this is to be consistent w/ output parameter names)
1321 if (k>1) vbf%Kd_itides(i,j,k) = vbf%Kd_itides(i,j,k) + 0.5*kd_add
1322 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 if (allocated(cs%dd%Kd_Itidal_work)) &
1326 cs%dd%Kd_itidal_work(i,j,k) = gv%H_to_RZ * tke_itide_lay
1327 if (allocated(cs%dd%Fl_itidal)) &
1328 cs%dd%Fl_itidal(i,j,k) = tke_itidal_rem(i)
1329
1330 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 kd_add = tke_to_kd(i,k) * tke_niku_lay
1334 if (kd_max >= 0.0) kd_add = min(kd_add, kd_max)
1335 if (allocated(cs%dd%Kd_Niku)) then
1336 if (k>1) cs%dd%Kd_Niku(i,j,k) = cs%dd%Kd_Niku(i,j,k) + 0.5*kd_add
1337 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 if (associated(vbf%Kd_Niku)) then
1340 if (k>1) vbf%Kd_Niku(i,j,k) = vbf%Kd_Niku(i,j,k) + 0.5*kd_add
1341 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 if (allocated(cs%dd%Kd_Niku_work)) &
1346 cs%dd%Kd_Niku_work(i,j,k) = gv%H_to_RZ * tke_niku_lay
1347
1348 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 kd_add = tke_to_kd(i,k) * tke_lowmode_lay
1352 if (kd_max >= 0.0) kd_add = min(kd_add, kd_max)
1353 if (allocated(cs%dd%Kd_lowmode)) then
1354 if (k>1) cs%dd%Kd_lowmode(i,j,k) = cs%dd%Kd_lowmode(i,j,k) + 0.5*kd_add
1355 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 if (associated(vbf%Kd_lowmode)) then
1358 if (k>1) vbf%Kd_lowmode(i,j,k) = vbf%Kd_lowmode(i,j,k) + 0.5*kd_add
1359 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 if (allocated(cs%dd%Kd_lowmode_work)) &
1363 cs%dd%Kd_lowmode_work(i,j,k) = gv%H_to_RZ * tke_lowmode_lay
1364 if (allocated(cs%dd%Fl_lowmode)) &
1365 cs%dd%Fl_lowmode(i,j,k) = tke_lowmode_rem(i)
1366 enddo ; enddo
1367 endif ! Simmons
1368
1369 ! Polzin:
1370 if ( use_polzin ) then
1371 do k=nz-1,2,-1 ; do i=is,ie
1372 if (max_tke(i,k) <= 0.0) cycle
1373 z_from_bot(i) = z_from_bot(i) + dz(i,k)
1374 if (cs%tidal_answer_date < 20190101) then
1375 if (n2_meanz(i) > 1.0e-14*us%T_to_s**2 ) then
1376 z_from_bot_wkb(i) = z_from_bot_wkb(i) + dz(i,k) * n2_lay(i,k) / n2_meanz(i)
1377 else ; z_from_bot_wkb(i) = 0 ; endif
1378 else
1379 if (dz(i,k) * n2_lay(i,k) < (1.0e14 * dztot_wkb(i)) * n2_meanz(i)) then
1380 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 ( z0_polzin_scaled(i) + z_from_bot_wkb(i) )
1387 z0_psl = z0_polzin_scaled(i)*cs%Decay_scale_factor_lee
1388 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 ( 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 tke_itide_lay = tke_itidal_rem(i) - tke_itidal_bot(i) *tke_frac_top(i)
1395 tke_niku_lay = tke_niku_rem(i) - tke_niku_bot(i) * tke_frac_top_lee(i)
1396 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 if (tke_itide_lay + tke_niku_lay + tke_lowmode_lay > max_tke(i,k)) then
1400 frac_used = max_tke(i,k) / (tke_itide_lay + tke_niku_lay + tke_lowmode_lay)
1401 tke_itide_lay = frac_used * tke_itide_lay
1402 tke_niku_lay = frac_used * tke_niku_lay
1403 tke_lowmode_lay = frac_used * tke_lowmode_lay
1404 endif
1405
1406 ! Calculate vertical flux available to bottom of layer above
1407 tke_itidal_rem(i) = tke_itidal_rem(i) - tke_itide_lay
1408 tke_niku_rem(i) = tke_niku_rem(i) - tke_niku_lay
1409 tke_lowmode_rem(i) = tke_lowmode_rem(i) - tke_lowmode_lay
1410
1411 ! Convert power to diffusivity
1412 kd_add = tke_to_kd(i,k) * (tke_itide_lay + tke_niku_lay + tke_lowmode_lay)
1413
1414 if (kd_max >= 0.0) kd_add = min(kd_add, kd_max)
1415 if (present(kd_lay)) then
1416 kd_lay(i,k) = kd_lay(i,k) + kd_add
1417 endif
1418
1419 if (present(kd_int)) then
1420 kd_int(i,k) = kd_int(i,k) + 0.5 * kd_add
1421 kd_int(i,k+1) = kd_int(i,k+1) + 0.5 * kd_add
1422 endif
1423
1424 ! diagnostics
1425 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 kd_add = tke_to_kd(i,k) * tke_itide_lay
1429 if (kd_max >= 0.0) kd_add = min(kd_add, kd_max)
1430 if (allocated(cs%dd%Kd_itidal)) then
1431 if (k>1) cs%dd%Kd_itidal(i,j,k) = cs%dd%Kd_itidal(i,j,k) + 0.5*kd_add
1432 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 if (associated(vbf%Kd_itides)) then
1435 !Not to be confused w/ Kd_itidal (this is to be consistent w/ output parameter names)
1436 if (k>1) vbf%Kd_itides(i,j,k) = vbf%Kd_itides(i,j,k) + 0.5*kd_add
1437 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 if (allocated(cs%dd%Kd_Itidal_work)) &
1441 cs%dd%Kd_itidal_work(i,j,k) = gv%H_to_RZ * tke_itide_lay
1442 if (allocated(cs%dd%Fl_itidal)) cs%dd%Fl_itidal(i,j,k) = tke_itidal_rem(i)
1443
1444 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 kd_add = tke_to_kd(i,k) * tke_niku_lay
1448 if (kd_max >= 0.0) kd_add = min(kd_add, kd_max)
1449 if (allocated(cs%dd%Kd_Niku)) then
1450 if (k>1) cs%dd%Kd_Niku(i,j,k) = cs%dd%Kd_Niku(i,j,k) + 0.5*kd_add
1451 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 if (associated(vbf%Kd_Niku)) then
1454 if (k>1) vbf%Kd_Niku(i,j,k) = vbf%Kd_Niku(i,j,k) + 0.5*kd_add
1455 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 if (allocated(cs%dd%Kd_Niku_work)) cs%dd%Kd_Niku_work(i,j,k) = gv%H_to_RZ * tke_niku_lay
1460
1461 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 kd_add = tke_to_kd(i,k) * tke_lowmode_lay
1465 if (kd_max >= 0.0) kd_add = min(kd_add, kd_max)
1466 if (allocated(cs%dd%Kd_lowmode)) then
1467 if (k>1) cs%dd%Kd_lowmode(i,j,k) = cs%dd%Kd_lowmode(i,j,k) + 0.5*kd_add
1468 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 if (associated(vbf%Kd_lowmode)) then
1471 if (k>1) vbf%Kd_lowmode(i,j,k) = vbf%Kd_lowmode(i,j,k) + 0.5*kd_add
1472 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 if (allocated(cs%dd%Kd_lowmode_work)) &
1476 cs%dd%Kd_lowmode_work(i,j,k) = gv%H_to_RZ * tke_lowmode_lay
1477 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
1482end subroutine add_int_tide_diffusivity
1483
1484!> Sets up diagnostics arrays for tidal mixing.
1485subroutine 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 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed ; nz = gv%ke
1494
1495 if ((cs%id_Kd_itidal > 0) .or. (cs%id_Kd_Itidal_work > 0)) &
1496 allocate(cs%dd%Kd_itidal(isd:ied,jsd:jed,nz+1), source=0.0)
1497 if ((cs%id_Kd_lowmode > 0) .or. (cs%id_Kd_lowmode_work > 0)) &
1498 allocate(cs%dd%Kd_lowmode(isd:ied,jsd:jed,nz+1), source=0.0)
1499 if (cs%id_Fl_itidal > 0) allocate(cs%dd%Fl_itidal(isd:ied,jsd:jed,nz+1), source=0.0)
1500 if (cs%id_Fl_lowmode > 0) allocate(cs%dd%Fl_lowmode(isd:ied,jsd:jed,nz+1), source=0.0)
1501 if (cs%id_Polzin_decay_scale > 0) allocate(cs%dd%Polzin_decay_scale(isd:ied,jsd:jed), source=0.0)
1502 if (cs%id_N2_bot > 0) allocate(cs%dd%N2_bot(isd:ied,jsd:jed), source=0.0)
1503 if (cs%id_N2_meanz > 0) allocate(cs%dd%N2_meanz(isd:ied,jsd:jed), source=0.0)
1504 if (cs%id_Polzin_decay_scale_scaled > 0) &
1505 allocate(cs%dd%Polzin_decay_scale_scaled(isd:ied,jsd:jed), source=0.0)
1506 if ((cs%id_Kd_Niku > 0) .or. (cs%id_Kd_Niku_work > 0)) &
1507 allocate(cs%dd%Kd_Niku(isd:ied,jsd:jed,nz+1), source=0.0)
1508 if (cs%id_Kd_Niku_work > 0) allocate(cs%dd%Kd_Niku_work(isd:ied,jsd:jed,nz), source=0.0)
1509 if (cs%id_Kd_Itidal_work > 0) allocate(cs%dd%Kd_Itidal_work(isd:ied,jsd:jed,nz), source=0.0)
1510 if (cs%id_Kd_Lowmode_Work > 0) allocate(cs%dd%Kd_Lowmode_Work(isd:ied,jsd:jed,nz), source=0.0)
1511 if (cs%id_TKE_itidal > 0) allocate(cs%dd%TKE_Itidal_used(isd:ied,jsd:jed), source=0.)
1512 ! additional diags for CVMix
1513 if (cs%id_N2_int > 0) allocate(cs%dd%N2_int(isd:ied,jsd:jed,nz+1), source=0.0)
1514 if (cs%id_Simmons_coeff > 0) then
1515 if (cs%CVMix_tidal_scheme /= simmons) then
1516 call mom_error(fatal, "setup_tidal_diagnostics: Simmons_coeff diagnostics is available "//&
1517 "only when CVMix_tidal_scheme is Simmons")
1518 endif
1519 allocate(cs%dd%Simmons_coeff_2d(isd:ied,jsd:jed), source=0.0)
1520 endif
1521 if (cs%id_vert_dep > 0) allocate(cs%dd%vert_dep_3d(isd:ied,jsd:jed,nz+1), source=0.0)
1522 if (cs%id_Schmittner_coeff > 0) then
1523 if (cs%CVMix_tidal_scheme /= schmittner) then
1524 call mom_error(fatal, "setup_tidal_diagnostics: Schmittner_coeff diagnostics is available "//&
1525 "only when CVMix_tidal_scheme is Schmittner.")
1526 endif
1527 allocate(cs%dd%Schmittner_coeff_3d(isd:ied,jsd:jed,nz), source=0.0)
1528 endif
1529 if (cs%id_tidal_qe_md > 0) then
1530 if (cs%CVMix_tidal_scheme /= schmittner) then
1531 call mom_error(fatal, "setup_tidal_diagnostics: tidal_qe_md diagnostics is available "//&
1532 "only when CVMix_tidal_scheme is Schmittner.")
1533 endif
1534 allocate(cs%dd%tidal_qe_md(isd:ied,jsd:jed,nz), source=0.0)
1535 endif
1536end subroutine setup_tidal_diagnostics
1537
1538!> This subroutine offers up diagnostics of the tidal mixing.
1539subroutine 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 if (cs%Int_tide_dissipation .or. cs%Lee_wave_dissipation .or. cs%Lowmode_itidal_dissipation) then
1547 if (cs%id_TKE_itidal > 0) call post_data(cs%id_TKE_itidal, cs%dd%TKE_itidal_used, cs%diag)
1548 if (cs%id_TKE_leewave > 0) call post_data(cs%id_TKE_leewave, cs%TKE_Niku, cs%diag)
1549 if (cs%id_Nb > 0) call post_data(cs%id_Nb, cs%Nb, cs%diag)
1550 if (cs%id_N2_bot > 0) call post_data(cs%id_N2_bot, cs%dd%N2_bot, cs%diag)
1551 if (cs%id_N2_meanz > 0) call post_data(cs%id_N2_meanz,cs%dd%N2_meanz,cs%diag)
1552
1553 if (cs%id_Fl_itidal > 0) call post_data(cs%id_Fl_itidal, cs%dd%Fl_itidal, cs%diag)
1554 if (cs%id_Kd_itidal > 0) call post_data(cs%id_Kd_itidal, cs%dd%Kd_itidal, cs%diag)
1555 if (cs%id_Kd_Niku > 0) call post_data(cs%id_Kd_Niku, cs%dd%Kd_Niku, cs%diag)
1556 if (cs%id_Kd_lowmode> 0) call post_data(cs%id_Kd_lowmode, cs%dd%Kd_lowmode, cs%diag)
1557 if (cs%id_Fl_lowmode> 0) call post_data(cs%id_Fl_lowmode, cs%dd%Fl_lowmode, cs%diag)
1558
1559 if (cs%id_N2_int> 0) call post_data(cs%id_N2_int, cs%dd%N2_int, cs%diag)
1560 if (cs%id_vert_dep> 0) call post_data(cs%id_vert_dep, cs%dd%vert_dep_3d, cs%diag)
1561 if (cs%id_Simmons_coeff> 0) call post_data(cs%id_Simmons_coeff, cs%dd%Simmons_coeff_2d, cs%diag)
1562 if (cs%id_Schmittner_coeff> 0) call post_data(cs%id_Schmittner_coeff, cs%dd%Schmittner_coeff_3d, cs%diag)
1563 if (cs%id_tidal_qe_md> 0) call post_data(cs%id_tidal_qe_md, cs%dd%tidal_qe_md, cs%diag)
1564
1565 if (cs%id_Kd_Itidal_Work > 0) &
1566 call post_data(cs%id_Kd_Itidal_Work, cs%dd%Kd_Itidal_Work, cs%diag)
1567 if (cs%id_Kd_Niku_Work > 0) call post_data(cs%id_Kd_Niku_Work, cs%dd%Kd_Niku_Work, cs%diag)
1568 if (cs%id_Kd_Lowmode_Work > 0) &
1569 call post_data(cs%id_Kd_Lowmode_Work, cs%dd%Kd_Lowmode_Work, cs%diag)
1570
1571 if (cs%id_Polzin_decay_scale > 0 ) &
1572 call post_data(cs%id_Polzin_decay_scale, cs%dd%Polzin_decay_scale, cs%diag)
1573 if (cs%id_Polzin_decay_scale_scaled > 0 ) &
1574 call post_data(cs%id_Polzin_decay_scale_scaled, cs%dd%Polzin_decay_scale_scaled, cs%diag)
1575 endif
1576
1577 if (allocated(cs%dd%Kd_itidal)) deallocate(cs%dd%Kd_itidal)
1578 if (allocated(cs%dd%Kd_lowmode)) deallocate(cs%dd%Kd_lowmode)
1579 if (allocated(cs%dd%Fl_itidal)) deallocate(cs%dd%Fl_itidal)
1580 if (allocated(cs%dd%Fl_lowmode)) deallocate(cs%dd%Fl_lowmode)
1581 if (allocated(cs%dd%Polzin_decay_scale)) deallocate(cs%dd%Polzin_decay_scale)
1582 if (allocated(cs%dd%Polzin_decay_scale_scaled)) deallocate(cs%dd%Polzin_decay_scale_scaled)
1583 if (allocated(cs%dd%N2_bot)) deallocate(cs%dd%N2_bot)
1584 if (allocated(cs%dd%N2_meanz)) deallocate(cs%dd%N2_meanz)
1585 if (allocated(cs%dd%Kd_Niku)) deallocate(cs%dd%Kd_Niku)
1586 if (allocated(cs%dd%Kd_Niku_work)) deallocate(cs%dd%Kd_Niku_work)
1587 if (allocated(cs%dd%Kd_Itidal_Work)) deallocate(cs%dd%Kd_Itidal_Work)
1588 if (allocated(cs%dd%Kd_Lowmode_Work)) deallocate(cs%dd%Kd_Lowmode_Work)
1589 if (allocated(cs%dd%TKE_itidal_used)) deallocate(cs%dd%TKE_itidal_used)
1590 if (allocated(cs%dd%N2_int)) deallocate(cs%dd%N2_int)
1591 if (allocated(cs%dd%vert_dep_3d)) deallocate(cs%dd%vert_dep_3d)
1592 if (allocated(cs%dd%Simmons_coeff_2d)) deallocate(cs%dd%Simmons_coeff_2d)
1593 if (allocated(cs%dd%Schmittner_coeff_3d)) deallocate(cs%dd%Schmittner_coeff_3d)
1594 if (allocated(cs%dd%tidal_qe_md)) deallocate(cs%dd%tidal_qe_md)
1595end subroutine post_tidal_diagnostics
1596
1597!> This subroutine returns a zonal slice of the topographic roughness amplitudes
1598subroutine 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 h_amp(:) = 0.0
1607 if ( cs%Int_tide_dissipation .and. .not. cs%use_CVMix_tidal ) then
1608 do i=g%isc,g%iec
1609 h_amp(i) = sqrt(cs%h2(i,j))
1610 enddo
1611 endif
1612
1613end 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.
1617subroutine 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 tidal_energy_flux_2d ! Input tidal energy flux at T-grid points [R Z3 T-3 ~> W m-2]
1632
1633 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 "Used with CVMix tidal mixing schemes.", fail_if_missing=.true.)
1638 tidal_energy_file = trim(cs%inputdir) // trim(tidal_energy_file)
1639
1640 select case (uppercase(tidal_energy_type(1:4)))
1641 case ('JAYN') ! Jayne 2009
1642 if (.not. allocated(cs%tidal_qe_2d)) allocate(cs%tidal_qe_2d(isd:ied,jsd:jed))
1643 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 default="wave_dissipation")
1647 call mom_read_data(tidal_energy_file, tidal_input_var, tidal_energy_flux_2d, g%domain, scale=us%W_m2_to_RZ3_T3)
1648 do j=g%jsc,g%jec ; do i=g%isc,g%iec
1649 cs%tidal_qe_2d(i,j) = cs%Gamma_itides * tidal_energy_flux_2d(i,j)
1650 enddo ; enddo
1651 deallocate(tidal_energy_flux_2d)
1652 case ('ER03') ! Egbert & Ray 2003
1653 call read_tidal_constituents(g, gv, us, tidal_energy_file, param_file, cs)
1654 case default
1655 call mom_error(fatal, "read_tidal_energy: Unknown tidal energy file type.")
1656 end select
1657
1658end subroutine read_tidal_energy
1659
1660!> This subroutine reads tidal input energy from a file by constituent.
1661subroutine 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 tidal_qk1, & ! qk1 coefficient used in Schmittner & Egbert [nondim]
1673 tidal_qo1 ! qo1 coefficient used in Schmittner & Egbert [nondim]
1674 real, allocatable, dimension(:) :: &
1675 z_t, & ! depth from surface to midpoint of input layer [Z ~> m]
1676 z_w ! depth from surface to top of input layer [Z ~> m]
1677 real, allocatable, dimension(:,:,:) :: &
1678 tc_m2, & ! input lunar semidiurnal tidal energy flux [R Z3 T-3 ~> W m-2]
1679 tc_s2, & ! input solar semidiurnal tidal energy flux [R Z3 T-3 ~> W m-2]
1680 tc_k1, & ! input lunar diurnal tidal energy flux [R Z3 T-3 ~> W m-2]
1681 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 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
1686 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
1687
1688 ! get number of input levels:
1689 call field_size(tidal_energy_file, 'z_t', nz_in)
1690
1691 ! allocate local variables
1692 allocate(z_t(nz_in(1)), z_w(nz_in(1)) )
1693 allocate(tc_m2(isd:ied,jsd:jed,nz_in(1)), &
1694 tc_s2(isd:ied,jsd:jed,nz_in(1)), &
1695 tc_k1(isd:ied,jsd:jed,nz_in(1)), &
1696 tc_o1(isd:ied,jsd:jed,nz_in(1)) )
1697
1698 ! allocate CS variables associated with 3d tidal energy dissipation
1699 if (.not. allocated(cs%tidal_qe_3d_in)) allocate(cs%tidal_qe_3d_in(isd:ied,jsd:jed,nz_in(1)))
1700 if (.not. allocated(cs%h_src)) allocate(cs%h_src(nz_in(1)))
1701
1702 ! read in tidal constituents
1703 call mom_read_data(tidal_energy_file, 'M2', tc_m2, g%domain, scale=us%W_m2_to_RZ3_T3)
1704 call mom_read_data(tidal_energy_file, 'S2', tc_s2, g%domain, scale=us%W_m2_to_RZ3_T3)
1705 call mom_read_data(tidal_energy_file, 'K1', tc_k1, g%domain, scale=us%W_m2_to_RZ3_T3)
1706 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 call mom_read_data(tidal_energy_file, 'z_t', z_t, scale=0.01*us%m_to_Z)
1709 call mom_read_data(tidal_energy_file, 'z_w', z_w, scale=0.01*us%m_to_Z)
1710
1711 do j=js,je ; do i=is,ie
1712 if (abs(g%geoLatT(i,j)) < 30.0) then
1713 tidal_qk1(i,j) = c1_3
1714 tidal_qo1(i,j) = c1_3
1715 else
1716 tidal_qk1(i,j) = 1.0
1717 tidal_qo1(i,j) = 1.0
1718 endif
1719 enddo ; enddo
1720
1721 cs%tidal_qe_3d_in(:,:,:) = 0.0
1722 do k=1,nz_in(1)
1723 ! Store the input cell thickness in m for use with CVmix.
1724 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 do j=js,je ; do i=is,ie
1727 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_3*tc_m2(i,j,k) + c1_3*tc_s2(i,j,k) + &
1729 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 if (any(cs%tidal_qe_3d_in<0.0)) then
1735 call mom_error(fatal, "read_tidal_constituents: Negative tidal_qe_3d_in terms.")
1736 endif
1737
1738 !! collapse 3D q*E to 2D q*E
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 h_neglect=gv%H_subroundoff, h_neglect_edge=gv%H_subroundoff)
1750
1751 deallocate(tc_m2)
1752 deallocate(tc_s2)
1753 deallocate(tc_k1)
1754 deallocate(tc_o1)
1755 deallocate(z_t)
1756 deallocate(z_w)
1757
1758end subroutine read_tidal_constituents
1759
1760!> Deallocate fields
1761subroutine 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 if (allocated(cs%tidal_qe_2d)) deallocate(cs%tidal_qe_2d)
1767 if (allocated(cs%tidal_qe_3d_in)) deallocate(cs%tidal_qe_3d_in)
1768 if (allocated(cs%h_src)) deallocate(cs%h_src)
1769end subroutine tidal_mixing_end
1770
1771end module mom_tidal_mixing