6 use mom_cpu_clock,
only : cpu_clock_id, cpu_clock_begin, cpu_clock_end
7 use mom_cpu_clock,
only : clock_module_driver, clock_module, clock_routine
23 implicit none ;
private
25 #include <MOM_memory.h>
27 public set_int_tide_input, int_tide_input_init, int_tide_input_end
44 real,
allocatable,
dimension(:,:) :: tke_itidal_coef
46 character(len=200) :: inputdir
48 logical :: int_tide_source_test
50 type(time_type) :: time_max_source
51 real :: int_tide_source_x
53 real :: int_tide_source_y
58 integer :: id_tke_itidal = -1, id_nb = -1, id_n2_bot = -1
64 real,
allocatable,
dimension(:,:) :: &
65 tke_itidal_input, & !< The internal tide TKE input at the bottom of the ocean [W m-2].
66 h2, & !< The squared topographic roughness height [Z2 ~> m2].
67 tideamp, & !< The amplitude of the tidal velocities [m s-1].
74 subroutine set_int_tide_input(u, v, h, tv, fluxes, itide, dt, G, GV, US, CS)
78 real,
dimension(SZIB_(G),SZJ_(G),SZK_(G)),
intent(in) :: u
79 real,
dimension(SZI_(G),SZJB_(G),SZK_(G)),
intent(in) :: v
80 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)),
intent(in) :: h
83 type(
forcing),
intent(in) :: fluxes
86 real,
intent(in) :: dt
89 real,
dimension(SZI_(G),SZJ_(G)) :: &
92 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)) :: &
97 logical :: avg_enabled
98 type(time_type) :: time_end
101 integer :: i, j, k, is, ie, js, je, nz
102 integer :: isd, ied, jsd, jed
105 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
106 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
108 if (.not.
associated(cs))
call mom_error(fatal,
"set_diffusivity: "//&
109 "Module must be initialized before it is used.")
111 use_eos =
associated(tv%eqn_of_state)
115 call vert_fill_ts(h, tv%T, tv%S, cs%kappa_fill*dt*us%s_to_T, t_f, s_f, g, gv, larger_h_denom=.true.)
118 call find_n2_bottom(h, tv, t_f, s_f, itide%h2, fluxes, g, gv, us, n2_bot)
121 do j=js,je ;
do i=is,ie
122 itide%Nb(i,j) = g%mask2dT(i,j) * us%s_to_T*sqrt(n2_bot(i,j))
123 itide%TKE_itidal_input(i,j) = min(cs%TKE_itidal_coef(i,j)*itide%Nb(i,j), cs%TKE_itide_max)
126 if (cs%int_tide_source_test)
then
127 itide%TKE_itidal_input(:,:) = 0.0
128 avg_enabled = query_averaging_enabled(cs%diag, time_end=time_end)
129 if (time_end <= cs%time_max_source)
then
130 do j=js,je ;
do i=is,ie
132 if (((g%geoLonCu(i-1,j)-cs%int_tide_source_x) * (g%geoLonBu(i,j)-cs%int_tide_source_x) <= 0.0) .and. &
133 ((g%geoLatCv(i,j-1)-cs%int_tide_source_y) * (g%geoLatCv(i,j)-cs%int_tide_source_y) <= 0.0))
then
134 itide%TKE_itidal_input(i,j) = 1.0
141 call hchksum(n2_bot,
"N2_bot",g%HI,haloshift=0, scale=us%s_to_T**2)
142 call hchksum(itide%TKE_itidal_input,
"TKE_itidal_input",g%HI,haloshift=0)
145 if (cs%id_TKE_itidal > 0)
call post_data(cs%id_TKE_itidal, itide%TKE_itidal_input, cs%diag)
146 if (cs%id_Nb > 0)
call post_data(cs%id_Nb, itide%Nb, cs%diag)
147 if (cs%id_N2_bot > 0 )
call post_data(cs%id_N2_bot, n2_bot, cs%diag)
149 end subroutine set_int_tide_input
152 subroutine find_n2_bottom(h, tv, T_f, S_f, h2, fluxes, G, GV, US, N2_bot)
156 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)),
intent(in) :: h
159 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)),
intent(in) :: T_f
161 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)),
intent(in) :: S_f
163 real,
dimension(SZI_(G),SZJ_(G)),
intent(in) :: h2
164 type(
forcing),
intent(in) :: fluxes
166 real,
dimension(SZI_(G),SZJ_(G)),
intent(out) :: N2_bot
169 real,
dimension(SZI_(G),SZK_(G)+1) :: &
171 real,
dimension(SZI_(G)) :: &
185 logical :: do_i(SZI_(G)), do_any
186 integer :: i, j, k, is, ie, js, je, nz
187 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
188 g_rho0 = (us%L_to_Z**2*gv%g_Earth) / gv%Rho0
192 drho_int(i,1) = 0.0 ; drho_int(i,nz+1) = 0.0
201 if (
associated(tv%eqn_of_state))
then
202 if (
associated(fluxes%p_surf))
then
203 do i=is,ie ; pres(i) = fluxes%p_surf(i,j) ;
enddo
205 do i=is,ie ; pres(i) = 0.0 ;
enddo
209 pres(i) = pres(i) + gv%H_to_Pa*h(i,j,k-1)
210 temp_int(i) = 0.5 * (t_f(i,j,k) + t_f(i,j,k-1))
211 salin_int(i) = 0.5 * (s_f(i,j,k) + s_f(i,j,k-1))
214 drho_dt(:), drho_ds(:), is, ie-is+1, tv%eqn_of_state)
216 drho_int(i,k) = max(drho_dt(i)*(t_f(i,j,k) - t_f(i,j,k-1)) + &
217 drho_ds(i)*(s_f(i,j,k) - s_f(i,j,k-1)), 0.0)
221 do k=2,nz ;
do i=is,ie
222 drho_int(i,k) = gv%Rlay(k) - gv%Rlay(k-1)
228 hb(i) = 0.0 ; drho_bot(i) = 0.0
229 z_from_bot(i) = 0.5*gv%H_to_Z*h(i,j,nz)
230 do_i(i) = (g%mask2dT(i,j) > 0.5)
231 h_amp(i) = sqrt(h2(i,j))
236 do i=is,ie ;
if (do_i(i))
then
237 dz_int = 0.5*gv%H_to_Z*(h(i,j,k) + h(i,j,k-1))
238 z_from_bot(i) = z_from_bot(i) + dz_int
240 hb(i) = hb(i) + dz_int
241 drho_bot(i) = drho_bot(i) + drho_int(i,k)
243 if (z_from_bot(i) > h_amp(i))
then
246 hb(i) = hb(i) + 0.5*gv%H_to_Z*(h(i,j,k-1) + h(i,j,k-2))
247 drho_bot(i) = drho_bot(i) + drho_int(i,k-1)
254 if (.not.do_any)
exit
258 if (hb(i) > 0.0)
then
259 n2_bot(i,j) = (g_rho0 * drho_bot(i)) / hb(i)
260 else ; n2_bot(i,j) = 0.0 ;
endif
264 end subroutine find_n2_bottom
267 subroutine int_tide_input_init(Time, G, GV, US, param_file, diag, CS, itide)
268 type(time_type),
intent(in) :: time
273 type(
diag_ctrl),
target,
intent(inout) :: diag
279 logical :: read_tideamp
281 #include "version_variable.h"
282 character(len=40) :: mdl =
"MOM_int_tide_input"
283 character(len=20) :: tmpstr
284 character(len=200) :: filename, tideamp_file, h2_file
287 real :: max_frac_rough
291 real :: kappa_h2_factor
293 real :: min_zbot_itides
296 integer :: i, j, is, ie, js, je, isd, ied, jsd, jed
298 if (
associated(cs))
then
299 call mom_error(warning,
"int_tide_input_init called with an associated "// &
300 "control structure.")
303 if (
associated(itide))
then
304 call mom_error(warning,
"int_tide_input_init called with an associated "// &
305 "internal tide input type.")
311 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
312 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
319 call get_param(param_file, mdl,
"INPUTDIR", cs%inputdir, default=
".")
320 cs%inputdir = slasher(cs%inputdir)
322 call get_param(param_file, mdl,
"DEBUG", cs%debug, default=.false., do_not_log=.true.)
324 call get_param(param_file, mdl,
"MIN_ZBOT_ITIDES", min_zbot_itides, &
325 "Turn off internal tidal dissipation when the total "//&
326 "ocean depth is less than this value.", units=
"m", default=0.0, scale=us%m_to_Z)
327 call get_param(param_file, mdl,
"KD_SMOOTH", cs%kappa_fill, &
328 "A diapycnal diffusivity that is used to interpolate "//&
329 "more sensible values of T & S into thin layers.", &
330 default=1.0e-6, scale=us%m2_s_to_Z2_T)
332 call get_param(param_file, mdl,
"UTIDE", utide, &
333 "The constant tidal amplitude used with INT_TIDE_DISSIPATION.", &
334 units=
"m s-1", default=0.0)
336 allocate(itide%Nb(isd:ied,jsd:jed)) ; itide%Nb(:,:) = 0.0
337 allocate(itide%h2(isd:ied,jsd:jed)) ; itide%h2(:,:) = 0.0
338 allocate(itide%TKE_itidal_input(isd:ied,jsd:jed)) ; itide%TKE_itidal_input(:,:) = 0.0
339 allocate(itide%tideamp(isd:ied,jsd:jed)) ; itide%tideamp(:,:) = utide
340 allocate(cs%TKE_itidal_coef(isd:ied,jsd:jed)) ; cs%TKE_itidal_coef(:,:) = 0.0
342 call get_param(param_file, mdl,
"KAPPA_ITIDES", kappa_itides, &
343 "A topographic wavenumber used with INT_TIDE_DISSIPATION. "//&
344 "The default is 2pi/10 km, as in St.Laurent et al. 2002.", &
345 units=
"m-1", default=8.e-4*atan(1.0))
347 call get_param(param_file, mdl,
"KAPPA_H2_FACTOR", kappa_h2_factor, &
348 "A scaling factor for the roughness amplitude with n"//&
349 "INT_TIDE_DISSIPATION.", units=
"nondim", default=1.0)
350 call get_param(param_file, mdl,
"TKE_ITIDE_MAX", cs%TKE_itide_max, &
351 "The maximum internal tide energy source available to mix "//&
352 "above the bottom boundary layer with INT_TIDE_DISSIPATION.", &
353 units=
"W m-2", default=1.0e3)
355 call get_param(param_file, mdl,
"READ_TIDEAMP", read_tideamp, &
356 "If true, read a file (given by TIDEAMP_FILE) containing "//&
357 "the tidal amplitude with INT_TIDE_DISSIPATION.", default=.false.)
358 if (read_tideamp)
then
359 call get_param(param_file, mdl,
"TIDEAMP_FILE", tideamp_file, &
360 "The path to the file containing the spatially varying "//&
361 "tidal amplitudes with INT_TIDE_DISSIPATION.", default=
"tideamp.nc")
362 filename = trim(cs%inputdir) // trim(tideamp_file)
363 call log_param(param_file, mdl,
"INPUTDIR/TIDEAMP_FILE", filename)
364 call mom_read_data(filename,
'tideamp', itide%tideamp, g%domain, timelevel=1)
367 call get_param(param_file, mdl,
"H2_FILE", h2_file, &
368 "The path to the file containing the sub-grid-scale "//&
369 "topographic roughness amplitude with INT_TIDE_DISSIPATION.", &
370 fail_if_missing=.true.)
371 filename = trim(cs%inputdir) // trim(h2_file)
372 call log_param(param_file, mdl,
"INPUTDIR/H2_FILE", filename)
373 call mom_read_data(filename,
'h2', itide%h2, g%domain, timelevel=1, scale=us%m_to_Z**2)
375 call get_param(param_file, mdl,
"FRACTIONAL_ROUGHNESS_MAX", max_frac_rough, &
376 "The maximum topographic roughness amplitude as a fraction of the mean depth, "//&
377 "or a negative value for no limitations on roughness.", &
378 units=
"nondim", default=0.1)
381 call get_param(param_file, mdl,
"INTERNAL_TIDE_SOURCE_TEST", cs%int_tide_source_test, &
382 "If true, apply an arbitrary generation site for internal tide testing", &
384 if (cs%int_tide_source_test)
then
385 call get_param(param_file, mdl,
"INTERNAL_TIDE_SOURCE_X", cs%int_tide_source_x, &
386 "X Location of generation site for internal tide", default=1.)
387 call get_param(param_file, mdl,
"INTERNAL_TIDE_SOURCE_Y", cs%int_tide_source_y, &
388 "Y Location of generation site for internal tide", default=1.)
389 call get_param(param_file, mdl,
"INTERNAL_TIDE_SOURCE_TLEN_DAYS", tlen_days, &
390 "Time interval from start of experiment for adding wave source", &
391 units=
"days", default=0)
392 cs%time_max_source = time + set_time(0, days=tlen_days)
395 do j=js,je ;
do i=is,ie
397 if (g%bathyT(i,j) < min_zbot_itides) mask_itidal = 0.0
399 itide%tideamp(i,j) = itide%tideamp(i,j) * mask_itidal * g%mask2dT(i,j)
402 if (max_frac_rough >= 0.0) &
403 itide%h2(i,j) = min((max_frac_rough*g%bathyT(i,j))**2, itide%h2(i,j))
406 cs%TKE_itidal_coef(i,j) = 0.5*kappa_h2_factor*gv%Rho0*&
407 kappa_itides * us%Z_to_m**2*itide%h2(i,j) * itide%tideamp(i,j)**2
411 cs%id_TKE_itidal = register_diag_field(
'ocean_model',
'TKE_itidal_itide',diag%axesT1,time, &
412 'Internal Tide Driven Turbulent Kinetic Energy',
'W m-2')
414 cs%id_Nb = register_diag_field(
'ocean_model',
'Nb_itide',diag%axesT1,time, &
415 'Bottom Buoyancy Frequency',
's-1')
417 cs%id_N2_bot = register_diag_field(
'ocean_model',
'N2_b_itide',diag%axesT1,time, &
418 'Bottom Buoyancy frequency squared',
's-2', conversion=us%s_to_T**2)
420 end subroutine int_tide_input_init
423 subroutine int_tide_input_end(CS)
426 if (
associated(cs))
deallocate(cs)
428 end subroutine int_tide_input_end