20 use cvmix_tidal,
only : cvmix_init_tidal, cvmix_compute_simmons_invariant
21 use cvmix_tidal,
only : cvmix_coeffs_tidal, cvmix_tidal_params_type
22 use cvmix_tidal,
only : cvmix_compute_schmittner_invariant, cvmix_compute_schmittnercoeff
23 use cvmix_tidal,
only : cvmix_coeffs_tidal_schmittner
24 use cvmix_kinds_and_types,
only : cvmix_global_params_type
25 use cvmix_put_get,
only : cvmix_put
27 implicit none ;
private
29 #include <MOM_memory.h>
44 real,
pointer,
dimension(:,:,:) :: &
48 kd_niku_work => null(),&
49 kd_itidal_work => null(),&
50 kd_lowmode_work => null(),&
52 vert_dep_3d => null(),&
53 schmittner_coeff_3d => null()
54 real,
pointer,
dimension(:,:,:) :: tidal_qe_md => null()
56 real,
pointer,
dimension(:,:,:) :: kd_lowmode => null()
58 real,
pointer,
dimension(:,:,:) :: fl_lowmode => null()
60 real,
pointer,
dimension(:,:) :: &
61 tke_itidal_used => null(),&
64 polzin_decay_scale_scaled => null(),&
65 polzin_decay_scale => null(),&
66 simmons_coeff_2d => null()
73 logical :: debug = .true.
76 logical :: int_tide_dissipation = .false.
79 integer :: int_tide_profile
82 logical :: lee_wave_dissipation = .false.
87 integer :: lee_wave_profile
91 real :: int_tide_decay_scale
100 real :: decay_scale_factor_lee
103 real :: min_zbot_itides
104 logical :: lowmode_itidal_dissipation = .false.
113 real :: nbotref_polzin
116 real :: polzin_decay_scale_factor
118 real :: polzin_decay_scale_max_factor
121 real :: polzin_min_decay_scale
124 real :: tke_itide_max
129 real :: kappa_h2_factor
130 character(len=200) :: inputdir
132 logical :: use_cvmix_tidal = .false.
135 real :: min_thickness
138 integer :: cvmix_tidal_scheme = -1
139 type(cvmix_tidal_params_type) :: cvmix_tidal_params
140 type(cvmix_global_params_type) :: cvmix_glb_params
141 real :: tidal_max_coef
142 real :: tidal_diss_lim_tc
147 real,
pointer,
dimension(:,:) :: tke_niku => null()
149 real,
pointer,
dimension(:,:) :: tke_itidal => null()
151 real,
pointer,
dimension(:,:) :: nb => null()
152 real,
pointer,
dimension(:,:) :: mask_itidal => null()
153 real,
pointer,
dimension(:,:) :: h2 => null()
154 real,
pointer,
dimension(:,:) :: tideamp => null()
155 real,
allocatable,
dimension(:) :: h_src
156 real,
allocatable,
dimension(:,:) :: tidal_qe_2d
160 real,
allocatable,
dimension(:,:,:) :: tidal_qe_3d_in
162 logical :: answers_2018
171 integer :: id_tke_itidal = -1
172 integer :: id_tke_leewave = -1
173 integer :: id_kd_itidal = -1
174 integer :: id_kd_niku = -1
175 integer :: id_kd_lowmode = -1
176 integer :: id_kd_itidal_work = -1
177 integer :: id_kd_niku_work = -1
178 integer :: id_kd_lowmode_work = -1
179 integer :: id_nb = -1
180 integer :: id_n2_bot = -1
181 integer :: id_n2_meanz = -1
182 integer :: id_fl_itidal = -1
183 integer :: id_fl_lowmode = -1
184 integer :: id_polzin_decay_scale = -1
185 integer :: id_polzin_decay_scale_scaled = -1
186 integer :: id_n2_int = -1
187 integer :: id_simmons_coeff = -1
188 integer :: id_schmittner_coeff = -1
189 integer :: id_tidal_qe_md = -1
190 integer :: id_vert_dep = -1
210 type(time_type),
intent(in) :: time
215 type(
diag_ctrl),
target,
intent(inout) :: diag
219 logical :: read_tideamp
220 logical :: default_2018_answers
221 character(len=20) :: tmpstr, int_tide_profile_str
222 character(len=20) :: cvmix_tidal_scheme_str, tidal_energy_type
223 character(len=200) :: filename, h2_file, niku_tke_input_file
224 character(len=200) :: tidal_energy_file, tideamp_file
225 real :: utide, hamp, prandtl_tidal, max_frac_rough
227 integer :: i, j, is, ie, js, je
228 integer :: isd, ied, jsd, jed
230 # include "version_variable.h"
231 character(len=40) :: mdl =
"MOM_tidal_mixing"
233 if (
associated(cs))
then
234 call mom_error(warning,
"tidal_mixing_init called when control structure "// &
235 "is already associated.")
243 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
244 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
250 "Vertical Tidal Mixing Parameterization")
251 call get_param(param_file, mdl,
"USE_CVMix_TIDAL", cs%use_CVMix_tidal, &
252 "If true, turns on tidal mixing via CVMix", &
255 call get_param(param_file, mdl,
"INPUTDIR", cs%inputdir, default=
".",do_not_log=.true.)
256 cs%inputdir = slasher(cs%inputdir)
257 call get_param(param_file, mdl,
"INT_TIDE_DISSIPATION", cs%int_tide_dissipation, &
258 "If true, use an internal tidal dissipation scheme to "//&
259 "drive diapycnal mixing, along the lines of St. Laurent "//&
260 "et al. (2002) and Simmons et al. (2004).", default=cs%use_CVMix_tidal)
266 call get_param(param_file, mdl,
"DEFAULT_2018_ANSWERS", default_2018_answers, &
267 "This sets the default value for the various _2018_ANSWERS parameters.", &
269 call get_param(param_file, mdl,
"TIDAL_MIXING_2018_ANSWERS", cs%answers_2018, &
270 "If true, use the order of arithmetic and expressions that recover the "//&
271 "answers from the end of 2018. Otherwise, use updated and more robust "//&
272 "forms of the same expressions.", default=default_2018_answers)
274 if (cs%int_tide_dissipation)
then
277 if (cs%use_CVMix_tidal)
then
278 call get_param(param_file, mdl,
"CVMIX_TIDAL_SCHEME", cvmix_tidal_scheme_str, &
279 "CVMIX_TIDAL_SCHEME selects the CVMix tidal mixing "//&
280 "scheme with INT_TIDE_DISSIPATION. Valid values are:\n"//&
281 "\t SIMMONS - Use the Simmons et al (2004) tidal \n"//&
282 "\t mixing scheme.\n"//&
283 "\t SCHMITTNER - Use the Schmittner et al (2014) tidal \n"//&
284 "\t mixing scheme.", &
286 cvmix_tidal_scheme_str =
uppercase(cvmix_tidal_scheme_str)
288 select case (cvmix_tidal_scheme_str)
292 call mom_error(fatal,
"tidal_mixing_init: Unrecognized setting "// &
293 "#define CVMIX_TIDAL_SCHEME "//trim(cvmix_tidal_scheme_str)//
" found in input file.")
298 if ( cs%CVMix_tidal_scheme.eq.
schmittner .or. .not. cs%use_CVMix_tidal)
then
299 call get_param(param_file, mdl,
"INT_TIDE_PROFILE", int_tide_profile_str, &
300 "INT_TIDE_PROFILE selects the vertical profile of energy "//&
301 "dissipation with INT_TIDE_DISSIPATION. Valid values are:\n"//&
302 "\t STLAURENT_02 - Use the St. Laurent et al exponential \n"//&
303 "\t decay profile.\n"//&
304 "\t POLZIN_09 - Use the Polzin WKB-stretched algebraic \n"//&
305 "\t decay profile.", &
307 int_tide_profile_str =
uppercase(int_tide_profile_str)
309 select case (int_tide_profile_str)
313 call mom_error(fatal,
"tidal_mixing_init: Unrecognized setting "// &
314 "#define INT_TIDE_PROFILE "//trim(int_tide_profile_str)//
" found in input file.")
318 elseif (cs%use_CVMix_tidal)
then
319 call mom_error(fatal,
"tidal_mixing_init: Cannot set INT_TIDE_DISSIPATION to False "// &
320 "when USE_CVMix_TIDAL is set to True.")
323 call get_param(param_file, mdl,
"LEE_WAVE_DISSIPATION", cs%Lee_wave_dissipation, &
324 "If true, use an lee wave driven dissipation scheme to "//&
325 "drive diapycnal mixing, along the lines of Nikurashin "//&
326 "(2010) and using the St. Laurent et al. (2002) "//&
327 "and Simmons et al. (2004) vertical profile", default=.false.)
328 if (cs%lee_wave_dissipation)
then
329 if (cs%use_CVMix_tidal)
then
330 call mom_error(fatal,
"tidal_mixing_init: Lee wave driven dissipation scheme cannot "// &
331 "be used when CVMix tidal mixing scheme is active.")
333 call get_param(param_file, mdl,
"LEE_WAVE_PROFILE", tmpstr, &
334 "LEE_WAVE_PROFILE selects the vertical profile of energy "//&
335 "dissipation with LEE_WAVE_DISSIPATION. Valid values are:\n"//&
336 "\t STLAURENT_02 - Use the St. Laurent et al exponential \n"//&
337 "\t decay profile.\n"//&
338 "\t POLZIN_09 - Use the Polzin WKB-stretched algebraic \n"//&
339 "\t decay profile.", &
346 call mom_error(fatal,
"tidal_mixing_init: Unrecognized setting "// &
347 "#define LEE_WAVE_PROFILE "//trim(tmpstr)//
" found in input file.")
351 call get_param(param_file, mdl,
"INT_TIDE_LOWMODE_DISSIPATION", cs%Lowmode_itidal_dissipation, &
352 "If true, consider mixing due to breaking low modes that "//&
353 "have been remotely generated; as with itidal drag on the "//&
354 "barotropic tide, use an internal tidal dissipation scheme to "//&
355 "drive diapycnal mixing, along the lines of St. Laurent "//&
356 "et al. (2002) and Simmons et al. (2004).", default=.false.)
358 if ((cs%Int_tide_dissipation .and. (cs%int_tide_profile ==
polzin_09)) .or. &
359 (cs%lee_wave_dissipation .and. (cs%lee_wave_profile ==
polzin_09)))
then
360 if (cs%use_CVMix_tidal)
then
361 call mom_error(fatal,
"tidal_mixing_init: Polzin scheme cannot "// &
362 "be used when CVMix tidal mixing scheme is active.")
364 call get_param(param_file, mdl,
"NU_POLZIN", cs%Nu_Polzin, &
365 "When the Polzin decay profile is used, this is a "//&
366 "non-dimensional constant in the expression for the "//&
367 "vertical scale of decay for the tidal energy dissipation.", &
368 units=
"nondim", default=0.0697)
369 call get_param(param_file, mdl,
"NBOTREF_POLZIN", cs%Nbotref_Polzin, &
370 "When the Polzin decay profile is used, this is the "//&
371 "reference value of the buoyancy frequency at the ocean "//&
372 "bottom in the Polzin formulation for the vertical "//&
373 "scale of decay for the tidal energy dissipation.", &
374 units=
"s-1", default=9.61e-4, scale=us%T_to_s)
375 call get_param(param_file, mdl,
"POLZIN_DECAY_SCALE_FACTOR", &
376 cs%Polzin_decay_scale_factor, &
377 "When the Polzin decay profile is used, this is a "//&
378 "scale factor for the vertical scale of decay of the tidal "//&
379 "energy dissipation.", default=1.0, units=
"nondim")
380 call get_param(param_file, mdl,
"POLZIN_SCALE_MAX_FACTOR", &
381 cs%Polzin_decay_scale_max_factor, &
382 "When the Polzin decay profile is used, this is a factor "//&
383 "to limit the vertical scale of decay of the tidal "//&
384 "energy dissipation to POLZIN_DECAY_SCALE_MAX_FACTOR "//&
385 "times the depth of the ocean.", units=
"nondim", default=1.0)
386 call get_param(param_file, mdl,
"POLZIN_MIN_DECAY_SCALE", cs%Polzin_min_decay_scale, &
387 "When the Polzin decay profile is used, this is the "//&
388 "minimum vertical decay scale for the vertical profile\n"//&
389 "of internal tide dissipation with the Polzin (2009) formulation", &
390 units=
"m", default=0.0, scale=us%m_to_Z)
393 if (cs%Int_tide_dissipation .or. cs%Lee_wave_dissipation)
then
394 call get_param(param_file, mdl,
"INT_TIDE_DECAY_SCALE", cs%Int_tide_decay_scale, &
395 "The decay scale away from the bottom for tidal TKE with "//&
396 "the new coding when INT_TIDE_DISSIPATION is used.", &
398 units=
"m", default=500.0, scale=us%m_to_Z)
399 call get_param(param_file, mdl,
"MU_ITIDES", cs%Mu_itides, &
400 "A dimensionless turbulent mixing efficiency used with "//&
401 "INT_TIDE_DISSIPATION, often 0.2.", units=
"nondim", default=0.2)
402 call get_param(param_file, mdl,
"GAMMA_ITIDES", cs%Gamma_itides, &
403 "The fraction of the internal tidal energy that is "//&
404 "dissipated locally with INT_TIDE_DISSIPATION. "//&
405 "THIS NAME COULD BE BETTER.", &
406 units=
"nondim", default=0.3333)
407 call get_param(param_file, mdl,
"MIN_ZBOT_ITIDES", cs%min_zbot_itides, &
408 "Turn off internal tidal dissipation when the total "//&
409 "ocean depth is less than this value.", units=
"m", default=0.0, scale=us%m_to_Z)
412 if ( (cs%Int_tide_dissipation .or. cs%Lee_wave_dissipation) .and. &
413 .not. cs%use_CVMix_tidal)
then
415 call safe_alloc_ptr(cs%Nb,isd,ied,jsd,jed)
416 call safe_alloc_ptr(cs%h2,isd,ied,jsd,jed)
417 call safe_alloc_ptr(cs%TKE_itidal,isd,ied,jsd,jed)
418 call safe_alloc_ptr(cs%mask_itidal,isd,ied,jsd,jed) ; cs%mask_itidal(:,:) = 1.0
420 call get_param(param_file, mdl,
"KAPPA_ITIDES", cs%kappa_itides, &
421 "A topographic wavenumber used with INT_TIDE_DISSIPATION. "//&
422 "The default is 2pi/10 km, as in St.Laurent et al. 2002.", &
423 units=
"m-1", default=8.e-4*atan(1.0), scale=us%Z_to_m)
425 call get_param(param_file, mdl,
"UTIDE", cs%utide, &
426 "The constant tidal amplitude used with INT_TIDE_DISSIPATION.", &
427 units=
"m s-1", default=0.0, scale=us%m_to_Z*us%T_to_s)
428 call safe_alloc_ptr(cs%tideamp,is,ie,js,je) ; cs%tideamp(:,:) = cs%utide
430 call get_param(param_file, mdl,
"KAPPA_H2_FACTOR", cs%kappa_h2_factor, &
431 "A scaling factor for the roughness amplitude with "//&
432 "INT_TIDE_DISSIPATION.", units=
"nondim", default=1.0)
433 call get_param(param_file, mdl,
"TKE_ITIDE_MAX", cs%TKE_itide_max, &
434 "The maximum internal tide energy source available to mix "//&
435 "above the bottom boundary layer with INT_TIDE_DISSIPATION.", &
436 units=
"W m-2", default=1.0e3, scale=us%kg_m3_to_R*us%m_to_Z**3*us%T_to_s**3)
438 call get_param(param_file, mdl,
"READ_TIDEAMP", read_tideamp, &
439 "If true, read a file (given by TIDEAMP_FILE) containing "//&
440 "the tidal amplitude with INT_TIDE_DISSIPATION.", default=.false.)
441 if (read_tideamp)
then
442 if (cs%use_CVMix_tidal)
then
443 call mom_error(fatal,
"tidal_mixing_init: Tidal amplitude files are "// &
444 "not compatible with CVMix tidal mixing. ")
446 call get_param(param_file, mdl,
"TIDEAMP_FILE", tideamp_file, &
447 "The path to the file containing the spatially varying "//&
448 "tidal amplitudes with INT_TIDE_DISSIPATION.", default=
"tideamp.nc")
449 filename = trim(cs%inputdir) // trim(tideamp_file)
450 call log_param(param_file, mdl,
"INPUTDIR/TIDEAMP_FILE", filename)
451 call mom_read_data(filename,
'tideamp', cs%tideamp, g%domain, timelevel=1, scale=us%m_to_Z*us%T_to_s)
454 call get_param(param_file, mdl,
"H2_FILE", h2_file, &
455 "The path to the file containing the sub-grid-scale "//&
456 "topographic roughness amplitude with INT_TIDE_DISSIPATION.", &
457 fail_if_missing=(.not.cs%use_CVMix_tidal))
458 filename = trim(cs%inputdir) // trim(h2_file)
459 call log_param(param_file, mdl,
"INPUTDIR/H2_FILE", filename)
460 call mom_read_data(filename,
'h2', cs%h2, g%domain, timelevel=1, scale=us%m_to_Z**2)
462 call get_param(param_file, mdl,
"FRACTIONAL_ROUGHNESS_MAX", max_frac_rough, &
463 "The maximum topographic roughness amplitude as a fraction of the mean depth, "//&
464 "or a negative value for no limitations on roughness.", &
465 units=
"nondim", default=0.1)
467 do j=js,je ;
do i=is,ie
468 if (g%bathyT(i,j) < cs%min_zbot_itides) cs%mask_itidal(i,j) = 0.0
469 cs%tideamp(i,j) = cs%tideamp(i,j) * cs%mask_itidal(i,j) * g%mask2dT(i,j)
472 if (cs%answers_2018 .and. (max_frac_rough >= 0.0))
then
473 hamp = min(max_frac_rough*g%bathyT(i,j), sqrt(cs%h2(i,j)))
474 cs%h2(i,j) = hamp*hamp
476 if (max_frac_rough >= 0.0) &
477 cs%h2(i,j) = min((max_frac_rough*g%bathyT(i,j))**2, cs%h2(i,j))
480 utide = cs%tideamp(i,j)
483 cs%TKE_itidal(i,j) = 0.5 * cs%kappa_h2_factor * gv%Rho0 * &
484 cs%kappa_itides * cs%h2(i,j) * utide*utide
489 if (cs%Lee_wave_dissipation)
then
491 call get_param(param_file, mdl,
"NIKURASHIN_TKE_INPUT_FILE",niku_tke_input_file, &
492 "The path to the file containing the TKE input from lee "//&
493 "wave driven mixing. Used with LEE_WAVE_DISSIPATION.", &
494 fail_if_missing=.true.)
495 call get_param(param_file, mdl,
"NIKURASHIN_SCALE",niku_scale, &
496 "A non-dimensional factor by which to scale the lee-wave "//&
497 "driven TKE input. Used with LEE_WAVE_DISSIPATION.", &
498 units=
"nondim", default=1.0)
500 filename = trim(cs%inputdir) // trim(niku_tke_input_file)
501 call log_param(param_file, mdl,
"INPUTDIR/NIKURASHIN_TKE_INPUT_FILE", &
503 call safe_alloc_ptr(cs%TKE_Niku,is,ie,js,je) ; cs%TKE_Niku(:,:) = 0.0
504 call mom_read_data(filename,
'TKE_input', cs%TKE_Niku, g%domain, timelevel=1, &
505 scale=us%kg_m3_to_R*us%m_to_Z**3*us%T_to_s**3)
506 cs%TKE_Niku(:,:) = niku_scale * cs%TKE_Niku(:,:)
508 call get_param(param_file, mdl,
"GAMMA_NIKURASHIN",cs%Gamma_lee, &
509 "The fraction of the lee wave energy that is dissipated "//&
510 "locally with LEE_WAVE_DISSIPATION.", units=
"nondim", &
512 call get_param(param_file, mdl,
"DECAY_SCALE_FACTOR_LEE",cs%Decay_scale_factor_lee, &
513 "Scaling for the vertical decay scaleof the local "//&
514 "dissipation of lee waves dissipation.", units=
"nondim", &
517 cs%Decay_scale_factor_lee = -9.e99
521 if (cs%use_CVMix_tidal)
then
525 call get_param(param_file, mdl,
"TIDAL_MAX_COEF", cs%tidal_max_coef, &
526 "largest acceptable value for tidal diffusivity", &
527 units=
"m^2/s", default=50e-4)
528 call get_param(param_file, mdl,
"TIDAL_DISS_LIM_TC", cs%tidal_diss_lim_tc, &
529 "Min allowable depth for dissipation for tidal-energy-constituent data. "//&
530 "No dissipation contribution is applied above TIDAL_DISS_LIM_TC.", &
531 units=
"m", default=0.0, scale=us%m_to_Z)
532 call get_param(param_file, mdl,
"TIDAL_ENERGY_FILE",tidal_energy_file, &
533 "The path to the file containing tidal energy "//&
534 "dissipation. Used with CVMix tidal mixing schemes.", &
535 fail_if_missing=.true.)
536 call get_param(param_file, mdl,
'MIN_THICKNESS', cs%min_thickness, default=0.001, &
538 call get_param(param_file, mdl,
"PRANDTL_TIDAL", prandtl_tidal, &
539 "Prandtl number used by CVMix tidal mixing schemes "//&
540 "to convert vertical diffusivities into viscosities.", &
541 units=
"nondim", default=1.0, &
543 call cvmix_put(cs%CVMix_glb_params,
'Prandtl',prandtl_tidal)
545 tidal_energy_file = trim(cs%inputdir) // trim(tidal_energy_file)
546 call get_param(param_file, mdl,
"TIDAL_ENERGY_TYPE",tidal_energy_type, &
547 "The type of input tidal energy flux dataset. Valid values are"//&
550 fail_if_missing=.true.)
553 (
uppercase(tidal_energy_type(1:4)).eq.
'JAYN' .and. cs%CVMix_tidal_scheme.eq.
simmons).or. &
554 (
uppercase(tidal_energy_type(1:4)).eq.
'ER03' .and. cs%CVMix_tidal_scheme.eq.
schmittner) ) )
then
555 call mom_error(fatal,
"tidal_mixing_init: Tidal energy file type ("//&
556 trim(tidal_energy_type)//
") is incompatible with CVMix tidal "//&
557 " mixing scheme: "//trim(cvmix_tidal_scheme_str) )
559 cvmix_tidal_scheme_str =
lowercase(cvmix_tidal_scheme_str)
562 call cvmix_init_tidal(cvmix_tidal_params_user = cs%CVMix_tidal_params, &
563 mix_scheme = cvmix_tidal_scheme_str, &
564 efficiency = cs%Mu_itides, &
565 vertical_decay_scale = cs%int_tide_decay_scale*us%Z_to_m, &
566 max_coefficient = cs%tidal_max_coef, &
567 local_mixing_frac = cs%Gamma_itides, &
568 depth_cutoff = cs%min_zbot_itides*us%Z_to_m)
578 if (cs%Int_tide_dissipation .or. cs%Lee_wave_dissipation .or. &
579 cs%Lowmode_itidal_dissipation)
then
581 cs%id_Kd_itidal = register_diag_field(
'ocean_model',
'Kd_itides',diag%axesTi,time, &
582 'Internal Tide Driven Diffusivity',
'm2 s-1', conversion=us%Z2_T_to_m2_s)
584 if (cs%use_CVMix_tidal)
then
585 cs%id_N2_int = register_diag_field(
'ocean_model',
'N2_int',diag%axesTi,time, &
586 'Bouyancy frequency squared, at interfaces',
's-2')
588 cs%id_Simmons_coeff = register_diag_field(
'ocean_model',
'Simmons_coeff',diag%axesT1,time, &
589 'time-invariant portion of the tidal mixing coefficient using the Simmons',
'')
590 cs%id_Schmittner_coeff = register_diag_field(
'ocean_model',
'Schmittner_coeff',diag%axesTL,time, &
591 'time-invariant portion of the tidal mixing coefficient using the Schmittner',
'')
592 cs%id_tidal_qe_md = register_diag_field(
'ocean_model',
'tidal_qe_md',diag%axesTL,time, &
593 'input tidal energy dissipated locally interpolated to model vertical coordinates',
'')
594 cs%id_vert_dep = register_diag_field(
'ocean_model',
'vert_dep',diag%axesTi,time, &
595 'vertical deposition function needed for Simmons et al tidal mixing',
'')
598 cs%id_TKE_itidal = register_diag_field(
'ocean_model',
'TKE_itidal',diag%axesT1,time, &
599 'Internal Tide Driven Turbulent Kinetic Energy', &
600 'W m-2', conversion=(us%R_to_kg_m3*us%Z_to_m**3*us%s_to_T**3))
601 cs%id_Nb = register_diag_field(
'ocean_model',
'Nb',diag%axesT1,time, &
602 'Bottom Buoyancy Frequency',
's-1', conversion=us%s_to_T)
604 cs%id_Kd_lowmode = register_diag_field(
'ocean_model',
'Kd_lowmode',diag%axesTi,time, &
605 'Internal Tide Driven Diffusivity (from propagating low modes)', &
606 'm2 s-1', conversion=us%Z2_T_to_m2_s)
608 cs%id_Fl_itidal = register_diag_field(
'ocean_model',
'Fl_itides',diag%axesTi,time, &
609 'Vertical flux of tidal turbulent dissipation', &
610 'm3 s-3', conversion=(us%Z_to_m**3*us%s_to_T**3))
612 cs%id_Fl_lowmode = register_diag_field(
'ocean_model',
'Fl_lowmode',diag%axesTi,time, &
613 'Vertical flux of tidal turbulent dissipation (from propagating low modes)', &
614 'm3 s-3', conversion=(us%Z_to_m**3*us%s_to_T**3))
616 cs%id_Polzin_decay_scale = register_diag_field(
'ocean_model',
'Polzin_decay_scale',diag%axesT1,time, &
617 'Vertical decay scale for the tidal turbulent dissipation with Polzin scheme', &
618 'm', conversion=us%Z_to_m)
620 cs%id_Polzin_decay_scale_scaled = register_diag_field(
'ocean_model', &
621 'Polzin_decay_scale_scaled', diag%axesT1, time, &
622 'Vertical decay scale for the tidal turbulent dissipation with Polzin scheme, '// &
623 'scaled by N2_bot/N2_meanz',
'm', conversion=us%Z_to_m)
625 cs%id_N2_bot = register_diag_field(
'ocean_model',
'N2_b',diag%axesT1,time, &
626 'Bottom Buoyancy frequency squared',
's-2', conversion=us%s_to_T**2)
628 cs%id_N2_meanz = register_diag_field(
'ocean_model',
'N2_meanz', diag%axesT1, time, &
629 'Buoyancy frequency squared averaged over the water column',
's-2', conversion=us%s_to_T**2)
631 cs%id_Kd_Itidal_Work = register_diag_field(
'ocean_model',
'Kd_Itidal_Work',diag%axesTL,time, &
632 'Work done by Internal Tide Diapycnal Mixing', &
633 'W m-2', conversion=(us%R_to_kg_m3*us%Z_to_m**3*us%s_to_T**3))
635 cs%id_Kd_Niku_Work = register_diag_field(
'ocean_model',
'Kd_Nikurashin_Work',diag%axesTL,time, &
636 'Work done by Nikurashin Lee Wave Drag Scheme', &
637 'W m-2', conversion=(us%R_to_kg_m3*us%Z_to_m**3*us%s_to_T**3))
639 cs%id_Kd_Lowmode_Work = register_diag_field(
'ocean_model',
'Kd_Lowmode_Work',diag%axesTL,time, &
640 'Work done by Internal Tide Diapycnal Mixing (low modes)', &
641 'W m-2', conversion=(us%R_to_kg_m3*us%Z_to_m**3*us%s_to_T**3))
643 if (cs%Lee_wave_dissipation)
then
644 cs%id_TKE_leewave = register_diag_field(
'ocean_model',
'TKE_leewave',diag%axesT1,time, &
645 'Lee wave Driven Turbulent Kinetic Energy', &
646 'W m-2', conversion=(us%R_to_kg_m3*us%Z_to_m**3*us%s_to_T**3))
647 cs%id_Kd_Niku = register_diag_field(
'ocean_model',
'Kd_Nikurashin',diag%axesTi,time, &
648 'Lee Wave Driven Diffusivity',
'm2 s-1', conversion=us%Z2_T_to_m2_s)
660 N2_lay, N2_int, Kd_lay, Kd_int, Kd_max, Kv)
664 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)), &
666 real,
dimension(SZI_(G)),
intent(in) :: n2_bot
668 real,
dimension(SZI_(G),SZK_(G)),
intent(in) :: n2_lay
670 real,
dimension(SZI_(G),SZK_(G)+1),
intent(in) :: n2_int
672 integer,
intent(in) :: j
673 real,
dimension(SZI_(G),SZK_(G)),
intent(in) :: tke_to_kd
678 real,
dimension(SZI_(G),SZK_(G)),
intent(in) :: max_tke
681 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)), &
682 intent(inout) :: kd_lay
683 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)+1), &
684 optional,
intent(inout) :: kd_int
686 real,
intent(in) :: kd_max
690 real,
dimension(:,:,:),
pointer :: kv
693 if (cs%Int_tide_dissipation .or. cs%Lee_wave_dissipation .or. cs%Lowmode_itidal_dissipation)
then
694 if (cs%use_CVMix_tidal)
then
698 g, gv, us, cs, n2_lay, kd_lay, kd_int, kd_max)
707 integer,
intent(in) :: j
712 real,
dimension(SZI_(G),SZK_(G)+1),
intent(in) :: N2_int
714 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)), &
716 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)), &
717 intent(inout) :: Kd_lay
718 real,
dimension(:,:,:),
pointer :: Kv
721 real,
dimension(SZK_(G)+1) :: Kd_tidal
722 real,
dimension(SZK_(G)+1) :: Kv_tidal
723 real,
dimension(SZK_(G)+1) :: vert_dep
724 real,
dimension(SZK_(G)+1) :: iFaceHeight
725 real,
dimension(SZK_(G)+1) :: SchmittnerSocn
726 real,
dimension(SZK_(G)) :: cellHeight
727 real,
dimension(SZK_(G)) :: tidal_qe_md
729 real,
dimension(SZK_(G)+1) :: N2_int_i
730 real,
dimension(SZK_(G)) :: Schmittner_coeff
731 real,
dimension(SZK_(G)) :: h_m
732 real,
allocatable,
dimension(:,:) :: exp_hab_zetar
734 integer :: i, k, is, ie
735 real :: dh, hcorr, Simmons_coeff
736 real,
parameter :: rho_fw = 1000.0
738 real :: h_neglect, h_neglect_edge
741 is = g%isc ; ie = g%iec
744 select case (cs%CVMix_tidal_scheme)
748 if (g%mask2dT(i,j)<1) cycle
754 dh = h(i,j,k) * gv%H_to_m
756 hcorr = min( dh - cs%min_thickness, 0. )
757 dh = max( dh, cs%min_thickness )
758 cellheight(k) = ifaceheight(k) - 0.5 * dh
759 ifaceheight(k+1) = ifaceheight(k) - dh
762 call cvmix_compute_simmons_invariant( nlev = g%ke, &
763 energy_flux = cs%tidal_qe_2d(i,j), &
765 simmonscoeff = simmons_coeff, &
766 vertdep = vert_dep, &
769 cvmix_tidal_params_user = cs%CVMix_tidal_params)
774 simmons_coeff = simmons_coeff / cs%Gamma_itides
779 n2_int_i(k) = us%s_to_T**2 * n2_int(i,k)
782 call cvmix_coeffs_tidal( mdiff_out = kv_tidal, &
783 tdiff_out = kd_tidal, &
785 oceandepth = -ifaceheight(g%ke+1),&
786 simmonscoeff = simmons_coeff, &
787 vert_dep = vert_dep, &
790 cvmix_params = cs%CVMix_glb_params, &
791 cvmix_tidal_params_user = cs%CVMix_tidal_params)
795 kd_lay(i,j,k) = kd_lay(i,j,k) + 0.5 * us%m2_s_to_Z2_T * (kd_tidal(k) + kd_tidal(k+1))
799 if (
associated(kv))
then
801 kv(i,j,k) = kv(i,j,k) + us%m2_s_to_Z2_T * kv_tidal(k)
806 if (
associated(dd%Kd_itidal))
then
807 dd%Kd_itidal(i,j,:) = us%m2_s_to_Z2_T*kd_tidal(:)
809 if (
associated(dd%N2_int))
then
810 dd%N2_int(i,j,:) = n2_int(i,:)
812 if (
associated(dd%Simmons_coeff_2d))
then
813 dd%Simmons_coeff_2d(i,j) = simmons_coeff
815 if (
associated(dd%vert_dep_3d))
then
816 dd%vert_dep_3d(i,j,:) = vert_dep(:)
826 allocate(exp_hab_zetar(g%ke+1,g%ke+1))
827 if (gv%Boussinesq)
then
828 h_neglect = gv%m_to_H*1.0e-30 ; h_neglect_edge = gv%m_to_H*1.0e-10
830 h_neglect = gv%kg_m2_to_H*1.0e-30 ; h_neglect_edge = gv%kg_m2_to_H*1.0e-10
836 if (g%mask2dT(i,j)<1) cycle
841 h_m(k) = h(i,j,k)*gv%H_to_m
844 hcorr = min( dh - cs%min_thickness, 0. )
845 dh = max( dh, cs%min_thickness )
846 cellheight(k) = ifaceheight(k) - 0.5 * dh
847 ifaceheight(k+1) = ifaceheight(k) - dh
853 call cvmix_compute_schmittner_invariant(nlev = g%ke, &
854 vertdep = vert_dep, &
855 efficiency = cs%Mu_itides, &
857 exp_hab_zetar = exp_hab_zetar, &
859 cvmix_tidal_params_user = cs%CVMix_tidal_params)
866 call remapping_core_h(cs%remap_cs,
size(cs%h_src), cs%h_src, cs%tidal_qe_3d_in(i,j,:), &
867 g%ke, h_m, tidal_qe_md)
871 call cvmix_compute_schmittnercoeff( nlev = g%ke, &
872 energy_flux = tidal_qe_md(:), &
874 schmittnercoeff = schmittner_coeff, &
875 exp_hab_zetar = exp_hab_zetar, &
876 cvmix_tidal_params_user = cs%CVMix_tidal_params)
880 n2_int_i(k) = us%s_to_T**2 * n2_int(i,k)
883 call cvmix_coeffs_tidal_schmittner( mdiff_out = kv_tidal, &
884 tdiff_out = kd_tidal, &
886 oceandepth = -ifaceheight(g%ke+1), &
887 vert_dep = vert_dep, &
890 schmittnercoeff = schmittner_coeff, &
891 schmittnersouthernocean = schmittnersocn, &
892 cvmix_params = cs%CVMix_glb_params, &
893 cvmix_tidal_params_user = cs%CVMix_tidal_params)
897 kd_lay(i,j,k) = kd_lay(i,j,k) + 0.5 * us%m2_s_to_Z2_T * (kd_tidal(k) + kd_tidal(k+1))
901 if (
associated(kv))
then
903 kv(i,j,k) = kv(i,j,k) + us%m2_s_to_Z2_T * kv_tidal(k)
908 if (
associated(dd%Kd_itidal))
then
909 dd%Kd_itidal(i,j,:) = us%m2_s_to_Z2_T*kd_tidal(:)
911 if (
associated(dd%N2_int))
then
912 dd%N2_int(i,j,:) = n2_int(i,:)
914 if (
associated(dd%Schmittner_coeff_3d))
then
915 dd%Schmittner_coeff_3d(i,j,:) = schmittner_coeff(:)
917 if (
associated(dd%tidal_qe_md))
then
918 dd%tidal_qe_md(i,j,:) = tidal_qe_md(:)
920 if (
associated(dd%vert_dep_3d))
then
921 dd%vert_dep_3d(i,j,:) = vert_dep(:)
925 deallocate(exp_hab_zetar)
928 call mom_error(fatal,
"tidal_mixing_init: Unrecognized setting "// &
929 "#define CVMIX_TIDAL_SCHEME found in input file.")
942 N2_lay, Kd_lay, Kd_int, Kd_max)
946 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)), &
948 real,
dimension(SZI_(G)),
intent(in) :: N2_bot
950 real,
dimension(SZI_(G),SZK_(G)),
intent(in) :: N2_lay
952 integer,
intent(in) :: j
953 real,
dimension(SZI_(G),SZK_(G)),
intent(in) :: TKE_to_Kd
958 real,
dimension(SZI_(G),SZK_(G)),
intent(in) :: max_TKE
961 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)), &
962 intent(inout) :: Kd_lay
963 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)+1), &
964 optional,
intent(inout) :: Kd_int
966 real,
intent(in) :: Kd_max
973 real,
dimension(SZI_(G)) :: &
994 tke_frac_top_lowmode, &
1001 real :: TKE_itide_lay
1002 real :: TKE_Niku_lay
1003 real :: TKE_lowmode_lay
1010 real :: TKE_lowmode_tot
1012 logical :: use_Polzin, use_Simmons
1013 character(len=160) :: mesg
1014 integer :: i, k, is, ie, nz
1018 is = g%isc ; ie = g%iec ; nz = g%ke
1021 if (.not.(cs%Int_tide_dissipation .or. cs%Lee_wave_dissipation))
return
1023 do i=is,ie ; htot(i) = 0.0 ; inv_int(i) = 0.0 ; inv_int_lee(i) = 0.0 ; inv_int_low(i) = 0.0 ;
enddo
1024 do k=1,nz ;
do i=is,ie
1025 htot(i) = htot(i) + gv%H_to_Z*h(i,j,k)
1028 i_rho0 = 1.0 / (gv%Rho0)
1030 use_polzin = ((cs%Int_tide_dissipation .and. (cs%int_tide_profile ==
polzin_09)) .or. &
1031 (cs%lee_wave_dissipation .and. (cs%lee_wave_profile ==
polzin_09)) .or. &
1032 (cs%Lowmode_itidal_dissipation .and. (cs%int_tide_profile ==
polzin_09)))
1033 use_simmons = ((cs%Int_tide_dissipation .and. (cs%int_tide_profile ==
stlaurent_02)) .or. &
1034 (cs%lee_wave_dissipation .and. (cs%lee_wave_profile ==
stlaurent_02)) .or. &
1035 (cs%Lowmode_itidal_dissipation .and. (cs%int_tide_profile ==
stlaurent_02)))
1039 if ( use_simmons )
then
1040 izeta = 1.0 / max(cs%Int_tide_decay_scale, gv%H_subroundoff*gv%H_to_Z)
1041 izeta_lee = 1.0 / max(cs%Int_tide_decay_scale*cs%Decay_scale_factor_lee, &
1042 gv%H_subroundoff*gv%H_to_Z)
1044 cs%Nb(i,j) = sqrt(n2_bot(i))
1045 if (
associated(dd%N2_bot)) dd%N2_bot(i,j) = n2_bot(i)
1046 if ( cs%Int_tide_dissipation )
then
1047 if (izeta*htot(i) > 1.0e-14)
then
1048 inv_int(i) = 1.0 / (1.0 - exp(-izeta*htot(i)))
1051 if ( cs%Lee_wave_dissipation )
then
1052 if (izeta_lee*htot(i) > 1.0e-14)
then
1053 inv_int_lee(i) = 1.0 / (1.0 - exp(-izeta_lee*htot(i)))
1056 if ( cs%Lowmode_itidal_dissipation)
then
1057 if (izeta*htot(i) > 1.0e-14)
then
1058 inv_int_low(i) = 1.0 / (1.0 - exp(-izeta*htot(i)))
1061 z_from_bot(i) = gv%H_to_Z*h(i,j,nz)
1066 if ( use_polzin )
then
1068 do i=is,ie ; n2_meanz(i) = 0.0 ;
enddo
1069 do k=1,nz ;
do i=is,ie
1070 n2_meanz(i) = n2_meanz(i) + n2_lay(i,k) * gv%H_to_Z * h(i,j,k)
1073 n2_meanz(i) = n2_meanz(i) / (htot(i) + gv%H_subroundoff*gv%H_to_Z)
1074 if (
associated(dd%N2_meanz)) dd%N2_meanz(i,j) = n2_meanz(i)
1078 do i=is,ie ; htot_wkb(i) = htot(i) ;
enddo
1086 cs%Nb(i,j) = sqrt(n2_bot(i))
1087 if (cs%answers_2018)
then
1088 if ((cs%tideamp(i,j) > 0.0) .and. &
1089 (cs%kappa_itides**2 * cs%h2(i,j) * cs%Nb(i,j)**3 > 1.0e-14*us%T_to_s**3) )
then
1090 z0_polzin(i) = cs%Polzin_decay_scale_factor * cs%Nu_Polzin * &
1091 cs%Nbotref_Polzin**2 * cs%tideamp(i,j) / &
1092 ( cs%kappa_itides**2 * cs%h2(i,j) * cs%Nb(i,j)**3 )
1093 if (z0_polzin(i) < cs%Polzin_min_decay_scale) &
1094 z0_polzin(i) = cs%Polzin_min_decay_scale
1095 if (n2_meanz(i) > 1.0e-14*us%T_to_s**2 )
then
1096 z0_polzin_scaled(i) = z0_polzin(i)*cs%Nb(i,j)**2 / n2_meanz(i)
1098 z0_polzin_scaled(i) = cs%Polzin_decay_scale_max_factor * htot(i)
1100 if (z0_polzin_scaled(i) > (cs%Polzin_decay_scale_max_factor * htot(i)) ) &
1101 z0_polzin_scaled(i) = cs%Polzin_decay_scale_max_factor * htot(i)
1103 z0_polzin(i) = cs%Polzin_decay_scale_max_factor * htot(i)
1104 z0_polzin_scaled(i) = cs%Polzin_decay_scale_max_factor * htot(i)
1107 z0ps_num = (cs%Polzin_decay_scale_factor * cs%Nu_Polzin * cs%Nbotref_Polzin**2) * cs%tideamp(i,j)
1108 z0ps_denom = ( cs%kappa_itides**2 * cs%h2(i,j) * cs%Nb(i,j) * n2_meanz(i) )
1109 if ((cs%tideamp(i,j) > 0.0) .and. &
1110 (z0ps_num < z0ps_denom * cs%Polzin_decay_scale_max_factor * htot(i)))
then
1111 z0_polzin_scaled(i) = z0ps_num / z0ps_denom
1113 if (abs(n2_meanz(i) * z0_polzin_scaled(i)) < &
1114 cs%Nb(i,j)**2 * (cs%Polzin_decay_scale_max_factor * htot(i)))
then
1115 z0_polzin(i) = z0_polzin_scaled(i) * (n2_meanz(i) / cs%Nb(i,j)**2)
1117 z0_polzin(i) = cs%Polzin_decay_scale_max_factor * htot(i)
1120 z0_polzin(i) = cs%Polzin_decay_scale_max_factor * htot(i)
1121 z0_polzin_scaled(i) = cs%Polzin_decay_scale_max_factor * htot(i)
1125 if (
associated(dd%Polzin_decay_scale)) &
1126 dd%Polzin_decay_scale(i,j) = z0_polzin(i)
1127 if (
associated(dd%Polzin_decay_scale_scaled)) &
1128 dd%Polzin_decay_scale_scaled(i,j) = z0_polzin_scaled(i)
1129 if (
associated(dd%N2_bot)) dd%N2_bot(i,j) = cs%Nb(i,j)*cs%Nb(i,j)
1131 if (cs%answers_2018)
then
1133 if ( cs%Int_tide_dissipation .and. (cs%int_tide_profile ==
polzin_09) )
then
1134 if (htot_wkb(i) > 1.0e-14*us%m_to_Z) &
1135 inv_int(i) = ( z0_polzin_scaled(i) / htot_wkb(i) ) + 1.0
1137 if ( cs%lee_wave_dissipation .and. (cs%lee_wave_profile ==
polzin_09) )
then
1138 if (htot_wkb(i) > 1.0e-14*us%m_to_Z) &
1139 inv_int_lee(i) = ( z0_polzin_scaled(i)*cs%Decay_scale_factor_lee / htot_wkb(i) ) + 1.0
1141 if ( cs%Lowmode_itidal_dissipation .and. (cs%int_tide_profile ==
polzin_09) )
then
1142 if (htot_wkb(i) > 1.0e-14*us%m_to_Z) &
1143 inv_int_low(i) = ( z0_polzin_scaled(i) / htot_wkb(i) ) + 1.0
1147 inv_int(i) = 0.0 ; inv_int_lee(i) = 0.0 ; inv_int_low(i) = 0.0
1148 if ( cs%Int_tide_dissipation .and. (cs%int_tide_profile ==
polzin_09) )
then
1149 if (z0_polzin_scaled(i) < 1.0e14 * htot_wkb(i)) &
1150 inv_int(i) = ( z0_polzin_scaled(i) / htot_wkb(i) ) + 1.0
1152 if ( cs%lee_wave_dissipation .and. (cs%lee_wave_profile ==
polzin_09) )
then
1153 if (z0_polzin_scaled(i) < 1.0e14 * htot_wkb(i)) &
1154 inv_int_lee(i) = ( z0_polzin_scaled(i)*cs%Decay_scale_factor_lee / htot_wkb(i) ) + 1.0
1156 if ( cs%Lowmode_itidal_dissipation .and. (cs%int_tide_profile ==
polzin_09) )
then
1157 if (z0_polzin_scaled(i) < 1.0e14 * htot_wkb(i)) &
1158 inv_int_low(i) = ( z0_polzin_scaled(i) / htot_wkb(i) ) + 1.0
1162 z_from_bot(i) = gv%H_to_Z*h(i,j,nz)
1164 if (cs%answers_2018)
then
1165 if (n2_meanz(i) > 1.0e-14*us%T_to_s**2 )
then
1166 z_from_bot_wkb(i) = gv%H_to_Z*h(i,j,nz) * n2_lay(i,nz) / n2_meanz(i)
1167 else ; z_from_bot_wkb(i) = 0 ;
endif
1169 if (gv%H_to_Z*h(i,j,nz) * n2_lay(i,nz) < n2_meanz(i) * (1.0e14 * htot_wkb(i)))
then
1170 z_from_bot_wkb(i) = gv%H_to_Z*h(i,j,nz) * n2_lay(i,nz) / n2_meanz(i)
1171 else ; z_from_bot_wkb(i) = 0 ;
endif
1180 tke_itidal_bot(i) = min(cs%TKE_itidal(i,j)*cs%Nb(i,j), cs%TKE_itide_max)
1181 if (
associated(dd%TKE_itidal_used)) &
1182 dd%TKE_itidal_used(i,j) = tke_itidal_bot(i)
1183 tke_itidal_bot(i) = (i_rho0 * cs%Mu_itides * cs%Gamma_itides) * tke_itidal_bot(i)
1185 tke_niku_bot(i) = 0.0
1186 if (cs%Lee_wave_dissipation)
then
1187 tke_niku_bot(i) = (i_rho0 * cs%Mu_itides * cs%Gamma_lee) * cs%TKE_Niku(i,j)
1190 tke_lowmode_tot = 0.0
1191 tke_lowmode_bot(i) = 0.0
1192 if (cs%Lowmode_itidal_dissipation)
then
1197 write (mesg,*)
"========", __file__, __line__
1198 call mom_error(fatal,trim(mesg)//
": this block not supported yet. (aa)")
1200 tke_lowmode_bot(i) = cs%Mu_itides * i_rho0 * tke_lowmode_tot
1203 tke_itidal_rem(i) = inv_int(i) * tke_itidal_bot(i)
1204 tke_niku_rem(i) = inv_int_lee(i) * tke_niku_bot(i)
1205 tke_lowmode_rem(i) = inv_int_low(i) * tke_lowmode_bot(i)
1207 if (
associated(dd%Fl_itidal)) dd%Fl_itidal(i,j,nz) = tke_itidal_rem(i)
1212 if ( use_simmons )
then
1213 do k=nz-1,2,-1 ;
do i=is,ie
1214 if (max_tke(i,k) <= 0.0) cycle
1215 z_from_bot(i) = z_from_bot(i) + gv%H_to_Z*h(i,j,k)
1218 tke_frac_top(i) = inv_int(i) * exp(-izeta * z_from_bot(i))
1219 tke_frac_top_lee(i) = inv_int_lee(i) * exp(-izeta_lee * z_from_bot(i))
1220 tke_frac_top_lowmode(i) = inv_int_low(i) * exp(-izeta * z_from_bot(i))
1224 tke_itide_lay = tke_itidal_rem(i) - tke_itidal_bot(i) * tke_frac_top(i)
1225 tke_niku_lay = tke_niku_rem(i) - tke_niku_bot(i) * tke_frac_top_lee(i)
1226 tke_lowmode_lay = tke_lowmode_rem(i) - tke_lowmode_bot(i)* tke_frac_top_lowmode(i)
1229 if (tke_itide_lay + tke_niku_lay + tke_lowmode_lay > (max_tke(i,k)))
then
1230 frac_used = (max_tke(i,k)) / (tke_itide_lay + tke_niku_lay + tke_lowmode_lay)
1231 tke_itide_lay = frac_used * tke_itide_lay
1232 tke_niku_lay = frac_used * tke_niku_lay
1233 tke_lowmode_lay = frac_used * tke_lowmode_lay
1237 tke_itidal_rem(i) = tke_itidal_rem(i) - tke_itide_lay
1238 tke_niku_rem(i) = tke_niku_rem(i) - tke_niku_lay
1239 tke_lowmode_rem(i) = tke_lowmode_rem(i) - tke_lowmode_lay
1242 kd_add = tke_to_kd(i,k) * (tke_itide_lay + tke_niku_lay + tke_lowmode_lay)
1244 if (kd_max >= 0.0) kd_add = min(kd_add, kd_max)
1245 kd_lay(i,j,k) = kd_lay(i,j,k) + kd_add
1247 if (
present(kd_int))
then
1248 kd_int(i,j,k) = kd_int(i,j,k) + 0.5 * kd_add
1249 kd_int(i,j,k+1) = kd_int(i,j,k+1) + 0.5 * kd_add
1253 if (
associated(dd%Kd_itidal))
then
1256 kd_add = tke_to_kd(i,k) * tke_itide_lay
1257 if (kd_max >= 0.0) kd_add = min(kd_add, kd_max)
1258 if (k>1) dd%Kd_itidal(i,j,k) = dd%Kd_itidal(i,j,k) + 0.5*kd_add
1259 if (k<nz) dd%Kd_itidal(i,j,k+1) = dd%Kd_itidal(i,j,k+1) + 0.5*kd_add
1261 if (
associated(dd%Kd_Itidal_work)) &
1262 dd%Kd_itidal_work(i,j,k) = gv%Rho0 * tke_itide_lay
1263 if (
associated(dd%Fl_itidal)) dd%Fl_itidal(i,j,k) = tke_itidal_rem(i)
1265 if (
associated(dd%Kd_Niku))
then
1268 kd_add = tke_to_kd(i,k) * tke_niku_lay
1269 if (kd_max >= 0.0) kd_add = min(kd_add, kd_max)
1270 if (k>1) dd%Kd_Niku(i,j,k) = dd%Kd_Niku(i,j,k) + 0.5*kd_add
1271 if (k<nz) dd%Kd_Niku(i,j,k+1) = dd%Kd_Niku(i,j,k+1) + 0.5*kd_add
1274 if (
associated(dd%Kd_Niku_work)) &
1275 dd%Kd_Niku_work(i,j,k) = gv%Rho0 * tke_niku_lay
1277 if (
associated(dd%Kd_lowmode))
then
1280 kd_add = tke_to_kd(i,k) * tke_lowmode_lay
1281 if (kd_max >= 0.0) kd_add = min(kd_add, kd_max)
1282 if (k>1) dd%Kd_lowmode(i,j,k) = dd%Kd_lowmode(i,j,k) + 0.5*kd_add
1283 if (k<nz) dd%Kd_lowmode(i,j,k+1) = dd%Kd_lowmode(i,j,k+1) + 0.5*kd_add
1285 if (
associated(dd%Kd_lowmode_work)) &
1286 dd%Kd_lowmode_work(i,j,k) = gv%Rho0 * tke_lowmode_lay
1287 if (
associated(dd%Fl_lowmode)) dd%Fl_lowmode(i,j,k) = tke_lowmode_rem(i)
1293 if ( use_polzin )
then
1294 do k=nz-1,2,-1 ;
do i=is,ie
1295 if (max_tke(i,k) <= 0.0) cycle
1296 z_from_bot(i) = z_from_bot(i) + gv%H_to_Z*h(i,j,k)
1297 if (n2_meanz(i) > 1.0e-14*us%T_to_s**2 )
then
1298 z_from_bot_wkb(i) = z_from_bot_wkb(i) &
1299 + gv%H_to_Z * h(i,j,k) * n2_lay(i,k) / n2_meanz(i)
1300 else ; z_from_bot_wkb(i) = 0 ;
endif
1303 tke_frac_top(i) = ( inv_int(i) * z0_polzin_scaled(i) ) / &
1304 ( z0_polzin_scaled(i) + z_from_bot_wkb(i) )
1305 z0_psl = z0_polzin_scaled(i)*cs%Decay_scale_factor_lee
1306 tke_frac_top_lee(i) = (inv_int_lee(i) * z0_psl) / (z0_psl + z_from_bot_wkb(i))
1307 tke_frac_top_lowmode(i) = ( inv_int_low(i) * z0_polzin_scaled(i) ) / &
1308 ( z0_polzin_scaled(i) + z_from_bot_wkb(i) )
1312 tke_itide_lay = tke_itidal_rem(i) - tke_itidal_bot(i) *tke_frac_top(i)
1313 tke_niku_lay = tke_niku_rem(i) - tke_niku_bot(i) * tke_frac_top_lee(i)
1314 tke_lowmode_lay = tke_lowmode_rem(i) - tke_lowmode_bot(i)*tke_frac_top_lowmode(i)
1317 if (tke_itide_lay + tke_niku_lay + tke_lowmode_lay > (max_tke(i,k)))
then
1318 frac_used = (max_tke(i,k)) / (tke_itide_lay + tke_niku_lay + tke_lowmode_lay)
1319 tke_itide_lay = frac_used * tke_itide_lay
1320 tke_niku_lay = frac_used * tke_niku_lay
1321 tke_lowmode_lay = frac_used * tke_lowmode_lay
1325 tke_itidal_rem(i) = tke_itidal_rem(i) - tke_itide_lay
1326 tke_niku_rem(i) = tke_niku_rem(i) - tke_niku_lay
1327 tke_lowmode_rem(i) = tke_lowmode_rem(i) - tke_lowmode_lay
1330 kd_add = tke_to_kd(i,k) * (tke_itide_lay + tke_niku_lay + tke_lowmode_lay)
1332 if (kd_max >= 0.0) kd_add = min(kd_add, kd_max)
1333 kd_lay(i,j,k) = kd_lay(i,j,k) + kd_add
1335 if (
present(kd_int))
then
1336 kd_int(i,j,k) = kd_int(i,j,k) + 0.5 * kd_add
1337 kd_int(i,j,k+1) = kd_int(i,j,k+1) + 0.5 * kd_add
1341 if (
associated(dd%Kd_itidal))
then
1344 kd_add = tke_to_kd(i,k) * tke_itide_lay
1345 if (kd_max >= 0.0) kd_add = min(kd_add, kd_max)
1346 if (k>1) dd%Kd_itidal(i,j,k) = dd%Kd_itidal(i,j,k) + 0.5*kd_add
1347 if (k<nz) dd%Kd_itidal(i,j,k+1) = dd%Kd_itidal(i,j,k+1) + 0.5*kd_add
1349 if (
associated(dd%Kd_Itidal_work)) &
1350 dd%Kd_itidal_work(i,j,k) = gv%Rho0 * tke_itide_lay
1351 if (
associated(dd%Fl_itidal)) dd%Fl_itidal(i,j,k) = tke_itidal_rem(i)
1353 if (
associated(dd%Kd_Niku))
then
1356 kd_add = tke_to_kd(i,k) * tke_niku_lay
1357 if (kd_max >= 0.0) kd_add = min(kd_add, kd_max)
1358 if (k>1) dd%Kd_Niku(i,j,k) = dd%Kd_Niku(i,j,k) + 0.5*kd_add
1359 if (k<nz) dd%Kd_Niku(i,j,k+1) = dd%Kd_Niku(i,j,k+1) + 0.5*kd_add
1362 if (
associated(dd%Kd_Niku_work)) dd%Kd_Niku_work(i,j,k) = gv%Rho0 * tke_niku_lay
1364 if (
associated(dd%Kd_lowmode))
then
1367 kd_add = tke_to_kd(i,k) * tke_lowmode_lay
1368 if (kd_max >= 0.0) kd_add = min(kd_add, kd_max)
1369 if (k>1) dd%Kd_lowmode(i,j,k) = dd%Kd_lowmode(i,j,k) + 0.5*kd_add
1370 if (k<nz) dd%Kd_lowmode(i,j,k+1) = dd%Kd_lowmode(i,j,k+1) + 0.5*kd_add
1372 if (
associated(dd%Kd_lowmode_work)) &
1373 dd%Kd_lowmode_work(i,j,k) = gv%Rho0 * tke_lowmode_lay
1374 if (
associated(dd%Fl_lowmode)) dd%Fl_lowmode(i,j,k) = tke_lowmode_rem(i)
1387 integer :: isd, ied, jsd, jed, nz
1390 isd = g%isd; ied = g%ied; jsd = g%jsd; jed = g%jed; nz = g%ke
1393 if ((cs%id_Kd_itidal > 0) .or. (cs%id_Kd_Itidal_work > 0))
then
1394 allocate(dd%Kd_itidal(isd:ied,jsd:jed,nz+1)) ; dd%Kd_itidal(:,:,:) = 0.0
1396 if ((cs%id_Kd_lowmode > 0) .or. (cs%id_Kd_lowmode_work > 0))
then
1397 allocate(dd%Kd_lowmode(isd:ied,jsd:jed,nz+1)) ; dd%Kd_lowmode(:,:,:) = 0.0
1399 if ( (cs%id_Fl_itidal > 0) )
then
1400 allocate(dd%Fl_itidal(isd:ied,jsd:jed,nz+1)) ; dd%Fl_itidal(:,:,:) = 0.0
1402 if ( (cs%id_Fl_lowmode > 0) )
then
1403 allocate(dd%Fl_lowmode(isd:ied,jsd:jed,nz+1)) ; dd%Fl_lowmode(:,:,:) = 0.0
1405 if ( (cs%id_Polzin_decay_scale > 0) )
then
1406 allocate(dd%Polzin_decay_scale(isd:ied,jsd:jed))
1407 dd%Polzin_decay_scale(:,:) = 0.0
1409 if ( (cs%id_N2_bot > 0) )
then
1410 allocate(dd%N2_bot(isd:ied,jsd:jed)) ; dd%N2_bot(:,:) = 0.0
1412 if ( (cs%id_N2_meanz > 0) )
then
1413 allocate(dd%N2_meanz(isd:ied,jsd:jed)) ; dd%N2_meanz(:,:) = 0.0
1415 if ( (cs%id_Polzin_decay_scale_scaled > 0) )
then
1416 allocate(dd%Polzin_decay_scale_scaled(isd:ied,jsd:jed))
1417 dd%Polzin_decay_scale_scaled(:,:) = 0.0
1419 if ((cs%id_Kd_Niku > 0) .or. (cs%id_Kd_Niku_work > 0))
then
1420 allocate(dd%Kd_Niku(isd:ied,jsd:jed,nz+1)) ; dd%Kd_Niku(:,:,:) = 0.0
1422 if (cs%id_Kd_Niku_work > 0)
then
1423 allocate(dd%Kd_Niku_work(isd:ied,jsd:jed,nz)) ; dd%Kd_Niku_work(:,:,:) = 0.0
1425 if (cs%id_Kd_Itidal_work > 0)
then
1426 allocate(dd%Kd_Itidal_work(isd:ied,jsd:jed,nz))
1427 dd%Kd_Itidal_work(:,:,:) = 0.0
1429 if (cs%id_Kd_Lowmode_Work > 0)
then
1430 allocate(dd%Kd_Lowmode_Work(isd:ied,jsd:jed,nz))
1431 dd%Kd_Lowmode_Work(:,:,:) = 0.0
1433 if (cs%id_TKE_itidal > 0)
then
1434 allocate(dd%TKE_Itidal_used(isd:ied,jsd:jed)) ; dd%TKE_Itidal_used(:,:) = 0.
1437 if (cs%id_N2_int > 0)
then
1438 allocate(dd%N2_int(isd:ied,jsd:jed,nz+1)) ; dd%N2_int(:,:,:) = 0.0
1440 if (cs%id_Simmons_coeff > 0)
then
1441 if (cs%CVMix_tidal_scheme .ne.
simmons)
then
1442 call mom_error(fatal,
"setup_tidal_diagnostics: Simmons_coeff diagnostics is available "//&
1443 "only when CVMix_tidal_scheme is Simmons")
1445 allocate(dd%Simmons_coeff_2d(isd:ied,jsd:jed)) ; dd%Simmons_coeff_2d(:,:) = 0.0
1447 if (cs%id_vert_dep > 0)
then
1448 allocate(dd%vert_dep_3d(isd:ied,jsd:jed,nz+1)) ; dd%vert_dep_3d(:,:,:) = 0.0
1450 if (cs%id_Schmittner_coeff > 0)
then
1451 if (cs%CVMix_tidal_scheme .ne.
schmittner)
then
1452 call mom_error(fatal,
"setup_tidal_diagnostics: Schmittner_coeff diagnostics is available "//&
1453 "only when CVMix_tidal_scheme is Schmittner.")
1455 allocate(dd%Schmittner_coeff_3d(isd:ied,jsd:jed,nz)) ; dd%Schmittner_coeff_3d(:,:,:) = 0.0
1457 if (cs%id_tidal_qe_md > 0)
then
1458 if (cs%CVMix_tidal_scheme .ne.
schmittner)
then
1459 call mom_error(fatal,
"setup_tidal_diagnostics: tidal_qe_md diagnostics is available "//&
1460 "only when CVMix_tidal_scheme is Schmittner.")
1462 allocate(dd%tidal_qe_md(isd:ied,jsd:jed,nz)) ; dd%tidal_qe_md(:,:,:) = 0.0
1470 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)), &
1479 if (cs%Int_tide_dissipation .or. cs%Lee_wave_dissipation .or. cs%Lowmode_itidal_dissipation)
then
1480 if (cs%id_TKE_itidal > 0)
call post_data(cs%id_TKE_itidal, dd%TKE_itidal_used, cs%diag)
1481 if (cs%id_TKE_leewave > 0)
call post_data(cs%id_TKE_leewave, cs%TKE_Niku, cs%diag)
1482 if (cs%id_Nb > 0)
call post_data(cs%id_Nb, cs%Nb, cs%diag)
1483 if (cs%id_N2_bot > 0)
call post_data(cs%id_N2_bot, dd%N2_bot, cs%diag)
1484 if (cs%id_N2_meanz > 0)
call post_data(cs%id_N2_meanz,dd%N2_meanz,cs%diag)
1486 if (cs%id_Fl_itidal > 0)
call post_data(cs%id_Fl_itidal, dd%Fl_itidal, cs%diag)
1487 if (cs%id_Kd_itidal > 0)
call post_data(cs%id_Kd_itidal, dd%Kd_itidal, cs%diag)
1488 if (cs%id_Kd_Niku > 0)
call post_data(cs%id_Kd_Niku, dd%Kd_Niku, cs%diag)
1489 if (cs%id_Kd_lowmode> 0)
call post_data(cs%id_Kd_lowmode, dd%Kd_lowmode, cs%diag)
1490 if (cs%id_Fl_lowmode> 0)
call post_data(cs%id_Fl_lowmode, dd%Fl_lowmode, cs%diag)
1492 if (cs%id_N2_int> 0)
call post_data(cs%id_N2_int, dd%N2_int, cs%diag)
1493 if (cs%id_vert_dep> 0)
call post_data(cs%id_vert_dep, dd%vert_dep_3d, cs%diag)
1494 if (cs%id_Simmons_coeff> 0)
call post_data(cs%id_Simmons_coeff, dd%Simmons_coeff_2d, cs%diag)
1495 if (cs%id_Schmittner_coeff> 0)
call post_data(cs%id_Schmittner_coeff, dd%Schmittner_coeff_3d, cs%diag)
1496 if (cs%id_tidal_qe_md> 0)
call post_data(cs%id_tidal_qe_md, dd%tidal_qe_md, cs%diag)
1498 if (cs%id_Kd_Itidal_Work > 0) &
1499 call post_data(cs%id_Kd_Itidal_Work, dd%Kd_Itidal_Work, cs%diag)
1500 if (cs%id_Kd_Niku_Work > 0)
call post_data(cs%id_Kd_Niku_Work, dd%Kd_Niku_Work, cs%diag)
1501 if (cs%id_Kd_Lowmode_Work > 0) &
1502 call post_data(cs%id_Kd_Lowmode_Work, dd%Kd_Lowmode_Work, cs%diag)
1504 if (cs%id_Polzin_decay_scale > 0 ) &
1505 call post_data(cs%id_Polzin_decay_scale, dd%Polzin_decay_scale, cs%diag)
1506 if (cs%id_Polzin_decay_scale_scaled > 0 ) &
1507 call post_data(cs%id_Polzin_decay_scale_scaled, dd%Polzin_decay_scale_scaled, cs%diag)
1510 if (
associated(dd%Kd_itidal))
deallocate(dd%Kd_itidal)
1511 if (
associated(dd%Kd_lowmode))
deallocate(dd%Kd_lowmode)
1512 if (
associated(dd%Fl_itidal))
deallocate(dd%Fl_itidal)
1513 if (
associated(dd%Fl_lowmode))
deallocate(dd%Fl_lowmode)
1514 if (
associated(dd%Polzin_decay_scale))
deallocate(dd%Polzin_decay_scale)
1515 if (
associated(dd%Polzin_decay_scale_scaled))
deallocate(dd%Polzin_decay_scale_scaled)
1516 if (
associated(dd%N2_bot))
deallocate(dd%N2_bot)
1517 if (
associated(dd%N2_meanz))
deallocate(dd%N2_meanz)
1518 if (
associated(dd%Kd_Niku))
deallocate(dd%Kd_Niku)
1519 if (
associated(dd%Kd_Niku_work))
deallocate(dd%Kd_Niku_work)
1520 if (
associated(dd%Kd_Itidal_Work))
deallocate(dd%Kd_Itidal_Work)
1521 if (
associated(dd%Kd_Lowmode_Work))
deallocate(dd%Kd_Lowmode_Work)
1522 if (
associated(dd%TKE_itidal_used))
deallocate(dd%TKE_itidal_used)
1523 if (
associated(dd%N2_int))
deallocate(dd%N2_int)
1524 if (
associated(dd%vert_dep_3d))
deallocate(dd%vert_dep_3d)
1525 if (
associated(dd%Simmons_coeff_2d))
deallocate(dd%Simmons_coeff_2d)
1526 if (
associated(dd%Schmittner_coeff_3d))
deallocate(dd%Schmittner_coeff_3d)
1527 if (
associated(dd%tidal_qe_md))
deallocate(dd%tidal_qe_md)
1536 character(len=20),
intent(in) :: tidal_energy_type
1537 character(len=200),
intent(in) :: tidal_energy_file
1540 integer :: i, j, isd, ied, jsd, jed, nz
1541 real,
allocatable,
dimension(:,:) :: tidal_energy_flux_2d
1543 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed ; nz = g%ke
1545 select case (
uppercase(tidal_energy_type(1:4)))
1547 if (.not.
allocated(cs%tidal_qe_2d))
allocate(cs%tidal_qe_2d(isd:ied,jsd:jed))
1548 allocate(tidal_energy_flux_2d(isd:ied,jsd:jed))
1549 call mom_read_data(tidal_energy_file,
'wave_dissipation',tidal_energy_flux_2d, g%domain)
1550 do j=g%jsc,g%jec ;
do i=g%isc,g%iec
1551 cs%tidal_qe_2d(i,j) = cs%Gamma_itides * tidal_energy_flux_2d(i,j)
1553 deallocate(tidal_energy_flux_2d)
1557 call mom_error(fatal,
"read_tidal_energy: Unknown tidal energy file type.")
1566 character(len=200),
intent(in) :: tidal_energy_file
1570 real,
parameter :: C1_3 = 1.0/3.0
1571 real,
dimension(SZI_(G),SZJ_(G)) :: &
1574 real,
allocatable,
dimension(:) :: &
1577 real,
allocatable,
dimension(:,:,:) :: &
1582 integer,
dimension(4) :: nz_in
1583 integer :: k, is, ie, js, je, isd, ied, jsd, jed, i, j
1585 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
1586 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
1589 call field_size(tidal_energy_file,
'z_t', nz_in)
1592 allocate(z_t(nz_in(1)), z_w(nz_in(1)) )
1593 allocate(tc_m2(isd:ied,jsd:jed,nz_in(1)), &
1594 tc_s2(isd:ied,jsd:jed,nz_in(1)), &
1595 tc_k1(isd:ied,jsd:jed,nz_in(1)), &
1596 tc_o1(isd:ied,jsd:jed,nz_in(1)) )
1599 if (.not.
allocated(cs%tidal_qe_3d_in))
allocate(cs%tidal_qe_3d_in(isd:ied,jsd:jed,nz_in(1)))
1600 if (.not.
allocated(cs%h_src))
allocate(cs%h_src(nz_in(1)))
1603 call mom_read_data(tidal_energy_file,
'M2', tc_m2, g%domain)
1604 call mom_read_data(tidal_energy_file,
'S2', tc_s2, g%domain)
1605 call mom_read_data(tidal_energy_file,
'K1', tc_k1, g%domain)
1606 call mom_read_data(tidal_energy_file,
'O1', tc_o1, g%domain)
1608 call mom_read_data(tidal_energy_file,
'z_t', z_t, scale=100.0*us%m_to_Z)
1609 call mom_read_data(tidal_energy_file,
'z_w', z_w, scale=100.0*us%m_to_Z)
1611 do j=js,je ;
do i=is,ie
1612 if (abs(g%geoLatT(i,j)) < 30.0)
then
1613 tidal_qk1(i,j) = c1_3
1614 tidal_qo1(i,j) = c1_3
1616 tidal_qk1(i,j) = 1.0
1617 tidal_qo1(i,j) = 1.0
1621 cs%tidal_qe_3d_in(:,:,:) = 0.0
1624 cs%h_src(k) = us%Z_to_m*(z_t(k)-z_w(k))*2.0
1626 do j=js,je ;
do i=is,ie
1627 if ((z_t(k) <= g%bathyT(i,j)) .and. (z_w(k) > cs%tidal_diss_lim_tc)) &
1628 cs%tidal_qe_3d_in(i,j,k) = c1_3*tc_m2(i,j,k) + c1_3*tc_s2(i,j,k) + &
1629 tidal_qk1(i,j)*tc_k1(i,j,k) + tidal_qo1(i,j)*tc_o1(i,j,k)
1649 if (any(cs%tidal_qe_3d_in<0.0))
then
1650 call mom_error(fatal,
"read_tidal_constituents: Negative tidal_qe_3d_in terms.")
1662 boundary_extrapolation=.false., check_remapping=cs%debug)
1678 if (.not.
associated(cs))
return
1681 if (
allocated(cs%tidal_qe_2d))
deallocate(cs%tidal_qe_2d)
1682 if (
allocated(cs%tidal_qe_3d_in))
deallocate(cs%tidal_qe_3d_in)
1683 if (
allocated(cs%h_src))
deallocate(cs%h_src)