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>
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 [R Z3 T-3 ~> W m-2].
66 h2, & !< The squared topographic roughness height [Z2 ~> m2].
67 tideamp, & !< The amplitude of the tidal velocities [Z T-1 ~> 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
100 integer :: i, j, k, is, ie, js, je, nz, isd, ied, jsd, jed
102 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
103 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
105 if (.not.
associated(cs))
call mom_error(fatal,
"set_diffusivity: "//&
106 "Module must be initialized before it is used.")
108 use_eos =
associated(tv%eqn_of_state)
112 call vert_fill_ts(h, tv%T, tv%S, cs%kappa_fill*dt, t_f, s_f, g, gv, larger_h_denom=.true.)
115 call find_n2_bottom(h, tv, t_f, s_f, itide%h2, fluxes, g, gv, us, n2_bot)
118 do j=js,je ;
do i=is,ie
119 itide%Nb(i,j) = g%mask2dT(i,j) * sqrt(n2_bot(i,j))
120 itide%TKE_itidal_input(i,j) = min(cs%TKE_itidal_coef(i,j)*itide%Nb(i,j), cs%TKE_itide_max)
123 if (cs%int_tide_source_test)
then
124 itide%TKE_itidal_input(:,:) = 0.0
125 avg_enabled = query_averaging_enabled(cs%diag, time_end=time_end)
126 if (time_end <= cs%time_max_source)
then
127 do j=js,je ;
do i=is,ie
129 if (((g%geoLonCu(i-1,j)-cs%int_tide_source_x) * (g%geoLonBu(i,j)-cs%int_tide_source_x) <= 0.0) .and. &
130 ((g%geoLatCv(i,j-1)-cs%int_tide_source_y) * (g%geoLatCv(i,j)-cs%int_tide_source_y) <= 0.0))
then
131 itide%TKE_itidal_input(i,j) = 1.0*us%kg_m3_to_R*us%m_to_Z**3*us%T_to_s**3
138 call hchksum(n2_bot,
"N2_bot",g%HI,haloshift=0, scale=us%s_to_T**2)
139 call hchksum(itide%TKE_itidal_input,
"TKE_itidal_input",g%HI,haloshift=0, &
140 scale=us%R_to_kg_m3*us%Z_to_m**3*us%s_to_T**3)
143 if (cs%id_TKE_itidal > 0)
call post_data(cs%id_TKE_itidal, itide%TKE_itidal_input, cs%diag)
144 if (cs%id_Nb > 0)
call post_data(cs%id_Nb, itide%Nb, cs%diag)
145 if (cs%id_N2_bot > 0 )
call post_data(cs%id_N2_bot, n2_bot, cs%diag)
150 subroutine find_n2_bottom(h, tv, T_f, S_f, h2, fluxes, G, GV, US, N2_bot)
154 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)),
intent(in) :: h
157 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)),
intent(in) :: T_f
159 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)),
intent(in) :: S_f
161 real,
dimension(SZI_(G),SZJ_(G)),
intent(in) :: h2
162 type(
forcing),
intent(in) :: fluxes
164 real,
dimension(SZI_(G),SZJ_(G)),
intent(out) :: N2_bot
167 real,
dimension(SZI_(G),SZK_(G)+1) :: &
169 real,
dimension(SZI_(G)) :: &
183 logical :: do_i(SZI_(G)), do_any
184 integer :: i, j, k, is, ie, js, je, nz
186 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
187 g_rho0 = (us%L_to_Z**2*gv%g_Earth) / gv%Rho0
191 drho_int(i,1) = 0.0 ; drho_int(i,nz+1) = 0.0
200 if (
associated(tv%eqn_of_state))
then
201 if (
associated(fluxes%p_surf))
then
202 do i=is,ie ; pres(i) = fluxes%p_surf(i,j) ;
enddo
204 do i=is,ie ; pres(i) = 0.0 ;
enddo
208 pres(i) = pres(i) + gv%H_to_Pa*h(i,j,k-1)
209 temp_int(i) = 0.5 * (t_f(i,j,k) + t_f(i,j,k-1))
210 salin_int(i) = 0.5 * (s_f(i,j,k) + s_f(i,j,k-1))
213 drho_dt(:), drho_ds(:), is, ie-is+1, tv%eqn_of_state, scale=us%kg_m3_to_R)
215 drho_int(i,k) = max(drho_dt(i)*(t_f(i,j,k) - t_f(i,j,k-1)) + &
216 drho_ds(i)*(s_f(i,j,k) - s_f(i,j,k-1)), 0.0)
220 do k=2,nz ;
do i=is,ie
221 drho_int(i,k) = (gv%Rlay(k) - gv%Rlay(k-1))
227 hb(i) = 0.0 ; drho_bot(i) = 0.0
228 z_from_bot(i) = 0.5*gv%H_to_Z*h(i,j,nz)
229 do_i(i) = (g%mask2dT(i,j) > 0.5)
230 h_amp(i) = sqrt(h2(i,j))
235 do i=is,ie ;
if (do_i(i))
then
236 dz_int = 0.5*gv%H_to_Z*(h(i,j,k) + h(i,j,k-1))
237 z_from_bot(i) = z_from_bot(i) + dz_int
239 hb(i) = hb(i) + dz_int
240 drho_bot(i) = drho_bot(i) + drho_int(i,k)
242 if (z_from_bot(i) > h_amp(i))
then
245 hb(i) = hb(i) + 0.5*gv%H_to_Z*(h(i,j,k-1) + h(i,j,k-2))
246 drho_bot(i) = drho_bot(i) + drho_int(i,k-1)
253 if (.not.do_any)
exit
257 if (hb(i) > 0.0)
then
258 n2_bot(i,j) = (g_rho0 * drho_bot(i)) / hb(i)
259 else ; n2_bot(i,j) = 0.0 ;
endif
267 type(time_type),
intent(in) :: time
272 type(
diag_ctrl),
target,
intent(inout) :: diag
278 logical :: read_tideamp
280 # include "version_variable.h"
281 character(len=40) :: mdl =
"MOM_int_tide_input"
282 character(len=20) :: tmpstr
283 character(len=200) :: filename, tideamp_file, h2_file
286 real :: max_frac_rough
290 real :: kappa_h2_factor
292 real :: min_zbot_itides
295 integer :: i, j, is, ie, js, je, isd, ied, jsd, jed
297 if (
associated(cs))
then
298 call mom_error(warning,
"int_tide_input_init called with an associated "// &
299 "control structure.")
302 if (
associated(itide))
then
303 call mom_error(warning,
"int_tide_input_init called with an associated "// &
304 "internal tide input type.")
310 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
311 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
318 call get_param(param_file, mdl,
"INPUTDIR", cs%inputdir, default=
".")
319 cs%inputdir = slasher(cs%inputdir)
321 call get_param(param_file, mdl,
"DEBUG", cs%debug, default=.false., do_not_log=.true.)
323 call get_param(param_file, mdl,
"MIN_ZBOT_ITIDES", min_zbot_itides, &
324 "Turn off internal tidal dissipation when the total "//&
325 "ocean depth is less than this value.", units=
"m", default=0.0, scale=us%m_to_Z)
326 call get_param(param_file, mdl,
"KD_SMOOTH", cs%kappa_fill, &
327 "A diapycnal diffusivity that is used to interpolate "//&
328 "more sensible values of T & S into thin layers.", &
329 default=1.0e-6, scale=us%m2_s_to_Z2_T)
331 call get_param(param_file, mdl,
"UTIDE", utide, &
332 "The constant tidal amplitude used with INT_TIDE_DISSIPATION.", &
333 units=
"m s-1", default=0.0, scale=us%m_s_to_L_T)
335 allocate(itide%Nb(isd:ied,jsd:jed)) ; itide%Nb(:,:) = 0.0
336 allocate(itide%h2(isd:ied,jsd:jed)) ; itide%h2(:,:) = 0.0
337 allocate(itide%TKE_itidal_input(isd:ied,jsd:jed)) ; itide%TKE_itidal_input(:,:) = 0.0
338 allocate(itide%tideamp(isd:ied,jsd:jed)) ; itide%tideamp(:,:) = utide
339 allocate(cs%TKE_itidal_coef(isd:ied,jsd:jed)) ; cs%TKE_itidal_coef(:,:) = 0.0
341 call get_param(param_file, mdl,
"KAPPA_ITIDES", kappa_itides, &
342 "A topographic wavenumber used with INT_TIDE_DISSIPATION. "//&
343 "The default is 2pi/10 km, as in St.Laurent et al. 2002.", &
344 units=
"m-1", default=8.e-4*atan(1.0), scale=us%L_to_m)
346 call get_param(param_file, mdl,
"KAPPA_H2_FACTOR", kappa_h2_factor, &
347 "A scaling factor for the roughness amplitude with n"//&
348 "INT_TIDE_DISSIPATION.", units=
"nondim", default=1.0)
349 call get_param(param_file, mdl,
"TKE_ITIDE_MAX", cs%TKE_itide_max, &
350 "The maximum internal tide energy source available to mix "//&
351 "above the bottom boundary layer with INT_TIDE_DISSIPATION.", &
352 units=
"W m-2", default=1.0e3, scale=us%kg_m3_to_R*us%m_to_Z**3*us%T_to_s**3)
354 call get_param(param_file, mdl,
"READ_TIDEAMP", read_tideamp, &
355 "If true, read a file (given by TIDEAMP_FILE) containing "//&
356 "the tidal amplitude with INT_TIDE_DISSIPATION.", default=.false.)
357 if (read_tideamp)
then
358 call get_param(param_file, mdl,
"TIDEAMP_FILE", tideamp_file, &
359 "The path to the file containing the spatially varying "//&
360 "tidal amplitudes with INT_TIDE_DISSIPATION.", default=
"tideamp.nc")
361 filename = trim(cs%inputdir) // trim(tideamp_file)
362 call log_param(param_file, mdl,
"INPUTDIR/TIDEAMP_FILE", filename)
363 call mom_read_data(filename,
'tideamp', itide%tideamp, g%domain, timelevel=1, scale=us%m_s_to_L_T)
366 call get_param(param_file, mdl,
"H2_FILE", h2_file, &
367 "The path to the file containing the sub-grid-scale "//&
368 "topographic roughness amplitude with INT_TIDE_DISSIPATION.", &
369 fail_if_missing=.true.)
370 filename = trim(cs%inputdir) // trim(h2_file)
371 call log_param(param_file, mdl,
"INPUTDIR/H2_FILE", filename)
372 call mom_read_data(filename,
'h2', itide%h2, g%domain, timelevel=1, scale=us%m_to_Z**2)
374 call get_param(param_file, mdl,
"FRACTIONAL_ROUGHNESS_MAX", max_frac_rough, &
375 "The maximum topographic roughness amplitude as a fraction of the mean depth, "//&
376 "or a negative value for no limitations on roughness.", &
377 units=
"nondim", default=0.1)
380 call get_param(param_file, mdl,
"INTERNAL_TIDE_SOURCE_TEST", cs%int_tide_source_test, &
381 "If true, apply an arbitrary generation site for internal tide testing", &
383 if (cs%int_tide_source_test)
then
384 call get_param(param_file, mdl,
"INTERNAL_TIDE_SOURCE_X", cs%int_tide_source_x, &
385 "X Location of generation site for internal tide", default=1.)
386 call get_param(param_file, mdl,
"INTERNAL_TIDE_SOURCE_Y", cs%int_tide_source_y, &
387 "Y Location of generation site for internal tide", default=1.)
388 call get_param(param_file, mdl,
"INTERNAL_TIDE_SOURCE_TLEN_DAYS", tlen_days, &
389 "Time interval from start of experiment for adding wave source", &
390 units=
"days", default=0)
391 cs%time_max_source = time + set_time(0, days=tlen_days)
394 do j=js,je ;
do i=is,ie
396 if (g%bathyT(i,j) < min_zbot_itides) mask_itidal = 0.0
398 itide%tideamp(i,j) = itide%tideamp(i,j) * mask_itidal * g%mask2dT(i,j)
401 if (max_frac_rough >= 0.0) &
402 itide%h2(i,j) = min((max_frac_rough*g%bathyT(i,j))**2, itide%h2(i,j))
405 cs%TKE_itidal_coef(i,j) = 0.5*us%L_to_Z*kappa_h2_factor*gv%Rho0*&
406 kappa_itides * itide%h2(i,j) * itide%tideamp(i,j)**2
410 cs%id_TKE_itidal =
register_diag_field(
'ocean_model',
'TKE_itidal_itide',diag%axesT1,time, &
411 'Internal Tide Driven Turbulent Kinetic Energy', &
412 'W m-2', conversion=us%R_to_kg_m3*us%Z_to_m**3*us%s_to_T**3)
415 'Bottom Buoyancy Frequency',
's-1', conversion=us%s_to_T)
418 'Bottom Buoyancy frequency squared',
's-2', conversion=us%s_to_T**2)
426 if (
associated(cs))
deallocate(cs)