MOM_internal_tide_input.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!> Calculates energy input to the internal tides
7
8use mom_cpu_clock, only : cpu_clock_id, cpu_clock_begin, cpu_clock_end
9use mom_cpu_clock, only : clock_module_driver, clock_module, clock_routine
12use mom_diag_mediator, only : safe_alloc_ptr, post_data, register_diag_field
13use mom_debugging, only : hchksum
14use mom_error_handler, only : mom_error, is_root_pe, fatal, warning, note
15use mom_file_parser, only : get_param, log_param, log_version, param_file_type
16use mom_file_parser, only : read_param
17use mom_forcing_type, only : forcing
18use mom_grid, only : ocean_grid_type
19use mom_io, only : slasher, vardesc, mom_read_data
20use mom_interface_heights, only : thickness_to_dz, find_rho_bottom
23use mom_time_manager, only : time_type, set_time, operator(+), operator(<=)
27use mom_eos, only : calculate_density_derivs, eos_domain
28
29implicit none ; private
30
31#include <MOM_memory.h>
32
35
36! A note on unit descriptions in comments: MOM6 uses units that can be rescaled for dimensional
37! consistency testing. These are noted in comments with units like Z, H, L, and T, along with
38! their mks counterparts with notation like "a velocity [Z T-1 ~> m s-1]". If the units
39! vary with the Boussinesq approximation, the Boussinesq variant is given first.
40
41!> This control structure holds parameters that regulate internal tide energy inputs.
42type, public :: int_tide_input_cs ; private
43 logical :: initialized = .false. !< True if this control structure has been initialized.
44 logical :: debug !< If true, write verbose checksums for debugging.
45 type(diag_ctrl), pointer :: diag => null() !< A structure that is used to
46 !! regulate the timing of diagnostic output.
47 real :: tke_itide_maxi !< Maximum Internal tide conversion
48 !! available to mix above the BBL [H Z2 T-3 ~> m3 s-3 or W m-2]
49 real :: kappa_fill !< Vertical diffusivity used to interpolate sensible values
50 !! of T & S into thin layers [H Z T-1 ~> m2 s-1 or kg m-1 s-1]
51
52 real, allocatable, dimension(:,:,:) :: tke_itidal_coef
53 !< The time-invariant field that enters the TKE_itidal input calculation noting that the
54 !! stratification and perhaps density are time-varying [R Z4 H-1 T-2 ~> J m-2 or J m kg-1].
55 real, allocatable, dimension(:,:,:) :: &
56 tke_itidal_input, & !< The internal tide TKE input at the bottom of the ocean [H Z2 T-3 ~> m3 s-3 or W m-2].
57 tideamp !< The amplitude of the tidal velocities [L T-1 ~> m s-1].
58
59 character(len=200) :: inputdir !< The directory for input files.
60
61 logical :: int_tide_source_test !< If true, apply an arbitrary generation site
62 !! for internal tide testing
63 type(time_type) :: time_max_source !< A time for use in testing internal tides
64 real :: int_tide_source_x !< X Location of generation site
65 !! for internal tide for testing [degrees_E] or [km]
66 real :: int_tide_source_y !< Y Location of generation site
67 !! for internal tide for testing [degrees_N] or [km]
68 integer :: int_tide_source_i !< I Location of generation site
69 integer :: int_tide_source_j !< J Location of generation site
70 logical :: int_tide_use_glob_ij !< Use global indices for generation site
71 integer :: nfreq = 0 !< The number of internal tide frequency bands
72
73
74 !>@{ Diagnostic IDs
75 integer, allocatable, dimension(:) :: id_tke_itidal_itide
76 integer :: id_nb = -1, id_n2_bot = -1
77 !>@}
78end type int_tide_input_cs
79
80!> This type is used to exchange fields related to the internal tides.
81type, public :: int_tide_input_type
82 real, allocatable, dimension(:,:) :: &
83 h2, & !< The squared topographic roughness height [Z2 ~> m2].
84 nb, & !< The bottom stratification [T-1 ~> s-1].
85 rho_bot !< The bottom density or the Boussinesq reference density [R ~> kg m-3].
87
88contains
89
90!> Sets the model-state dependent internal tide energy sources.
91subroutine set_int_tide_input(u, v, h, tv, fluxes, itide, dt, G, GV, US, CS)
92 type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure
93 type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure
94 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
95 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), intent(in) :: u !< The zonal velocity [L T-1 ~> m s-1]
96 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), intent(in) :: v !< The meridional velocity [L T-1 ~> m s-1]
97 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2]
98 type(thermo_var_ptrs), intent(in) :: tv !< A structure containing pointers to the
99 !! thermodynamic fields
100 type(forcing), intent(in) :: fluxes !< A structure of thermodynamic surface fluxes
101 type(int_tide_input_type), intent(inout) :: itide !< A structure containing fields related
102 !! to the internal tide sources.
103 real, intent(in) :: dt !< The time increment [T ~> s].
104 type(int_tide_input_cs), pointer :: cs !< This module's control structure.
105
106 ! Local variables
107 real, dimension(SZI_(G),SZJ_(G)) :: &
108 n2_bot ! The bottom squared buoyancy frequency [T-2 ~> s-2].
109 real, dimension(SZI_(G),SZJ_(G)) :: &
110 rho_bot, & ! The average near-bottom density or the Boussinesq reference density [R ~> kg m-3].
111 h_bot ! Bottom boundary layer thickness [H ~> m or kg m-2].
112 integer, dimension(SZI_(G),SZJ_(G)) :: k_bot ! Bottom boundary layer top layer index.
113
114 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: &
115 t_f, s_f ! The temperature and salinity in [C ~> degC] and [S ~> ppt] with the values in
116 ! the massless layers filled vertically by diffusion.
117 logical :: use_eos ! If true, density is calculated from T & S using an
118 ! equation of state.
119 logical :: avg_enabled ! for testing internal tides (BDM)
120 type(time_type) :: time_end !< For use in testing internal tides (BDM)
121 real :: hz2_t3_to_w_m2 ! unit conversion factor for TKE from internal units
122 ! to mks [T3 kg H-1 Z-2 s-3 ~> kg m-3 or 1]
123 real :: w_m2_to_hz2_t3 ! unit conversion factor for TKE from mks to internal
124 ! units [H Z2 s3 T-3 kg-1 ~> m3 kg-1 or 1]
125
126 integer :: i, j, is, ie, js, je, nz, isd, ied, jsd, jed
127 integer :: i_global, j_global
128 integer :: fr
129
130 hz2_t3_to_w_m2 = gv%H_to_kg_m2*(us%Z_to_m**2)*(us%s_to_T**3)
131 w_m2_to_hz2_t3 = gv%kg_m2_to_H*(us%m_to_Z**2)*(us%T_to_s**3)
132
133 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
134 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
135
136 if (.not.associated(cs)) call mom_error(fatal,"set_diffusivity: "//&
137 "Module must be initialized before it is used.")
138
139 if (.not.cs%initialized) call mom_error(fatal,"set_diffusivity: "//&
140 "Module must be initialized before it is used.")
141
142 use_eos = associated(tv%eqn_of_state)
143
144 ! Smooth the properties through massless layers.
145 if (use_eos) then
146 call vert_fill_ts(h, tv%T, tv%S, cs%kappa_fill*dt, t_f, s_f, g, gv, us, larger_h_denom=.true.)
147 endif
148
149 call find_n2_bottom(g, gv, us, tv, fluxes, h, t_f, s_f, itide%h2, n2_bot, rho_bot, h_bot, k_bot)
150
151 avg_enabled = query_averaging_enabled(cs%diag, time_end=time_end)
152
153 if (gv%Boussinesq .or. gv%semi_Boussinesq) then
154 !$OMP parallel do default(shared)
155 do fr=1,cs%nFreq ; do j=js,je ; do i=is,ie
156 itide%Nb(i,j) = g%mask2dT(i,j) * sqrt(n2_bot(i,j))
157 cs%TKE_itidal_input(i,j,fr) = min(gv%RZ_to_H*gv%Z_to_H*cs%TKE_itidal_coef(i,j,fr)*itide%Nb(i,j), &
158 cs%TKE_itide_maxi)
159 enddo ; enddo ; enddo
160 else
161 !$OMP parallel do default(shared)
162 do fr=1,cs%nFreq ; do j=js,je ; do i=is,ie
163 itide%Nb(i,j) = g%mask2dT(i,j) * sqrt(n2_bot(i,j))
164 itide%Rho_bot(i,j) = g%mask2dT(i,j) * rho_bot(i,j)
165 cs%TKE_itidal_input(i,j,fr) = min((gv%RZ_to_H*gv%RZ_to_H*rho_bot(i,j))*cs%TKE_itidal_coef(i,j,fr)*itide%Nb(i,j), &
166 cs%TKE_itide_maxi)
167 enddo ; enddo ; enddo
168 endif
169
170 if (cs%int_tide_source_test) then
171 cs%TKE_itidal_input(:,:,:) = 0.0
172 if (time_end <= cs%time_max_source) then
173 if (cs%int_tide_use_glob_ij) then
174 do fr=1,cs%nFreq ; do j=js,je ; do i=is,ie
175 i_global = i + g%idg_offset
176 j_global = j + g%jdg_offset
177 if ((i_global == cs%int_tide_source_i) .and. (j_global == cs%int_tide_source_j)) then
178 cs%TKE_itidal_input(i,j,fr) = 1.0*w_m2_to_hz2_t3
179 endif
180 enddo ; enddo ; enddo
181 else
182 do fr=1,cs%nFreq ; do j=js,je ; do i=is,ie
183 ! Input an arbitrary energy point source.id_
184 if (((g%geoLonCu(i-1,j)-cs%int_tide_source_x) * (g%geoLonBu(i,j)-cs%int_tide_source_x) <= 0.0) .and. &
185 ((g%geoLatCv(i,j-1)-cs%int_tide_source_y) * (g%geoLatCv(i,j)-cs%int_tide_source_y) <= 0.0)) then
186 cs%TKE_itidal_input(i,j,fr) = 1.0*w_m2_to_hz2_t3
187 endif
188 enddo ; enddo ; enddo
189 endif
190 endif
191 endif
192
193 if (cs%debug) then
194 call hchksum(n2_bot, "N2_bot", g%HI, haloshift=0, unscale=us%s_to_T**2)
195 call hchksum(cs%TKE_itidal_input,"TKE_itidal_input", g%HI, haloshift=0, &
196 unscale=hz2_t3_to_w_m2)
197 endif
198
199 call enable_averages(dt, time_end, cs%diag)
200
201 do fr=1,cs%nFreq
202 if (cs%id_TKE_itidal_itide(fr) > 0) call post_data(cs%id_TKE_itidal_itide(fr), &
203 cs%TKE_itidal_input(isd:ied,jsd:jed,fr), cs%diag)
204 enddo
205 if (cs%id_Nb > 0) call post_data(cs%id_Nb, itide%Nb, cs%diag)
206 if (cs%id_N2_bot > 0 ) call post_data(cs%id_N2_bot, n2_bot, cs%diag)
207
208 call disable_averaging(cs%diag)
209
210end subroutine set_int_tide_input
211
212!> Estimates the near-bottom buoyancy frequency (N^2).
213subroutine find_n2_bottom(G, GV, US, tv, fluxes, h, T_f, S_f, h2, N2_bot, Rho_bot, h_bot, k_bot)
214 type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure
215 type(verticalgrid_type), intent(in) :: GV !< The ocean's vertical grid structure
216 type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
217 type(thermo_var_ptrs), intent(in) :: tv !< A structure containing pointers to the
218 !! thermodynamic fields
219 type(forcing), intent(in) :: fluxes !< A structure of thermodynamic surface fluxes
220 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2]
221 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: T_f !< Temperature after vertical filtering to
222 !! smooth out the values in thin layers [C ~> degC].
223 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: S_f !< Salinity after vertical filtering to
224 !! smooth out the values in thin layers [S ~> ppt].
225 real, dimension(SZI_(G),SZJ_(G)), intent(in) :: h2 !< Bottom topographic roughness [Z2 ~> m2].
226 real, dimension(SZI_(G),SZJ_(G)), intent(out) :: N2_bot !< The squared buoyancy frequency at the
227 !! ocean bottom [T-2 ~> s-2].
228 real, dimension(SZI_(G),SZJ_(G)), intent(out) :: Rho_bot !< The average density near the ocean
229 !! bottom [R ~> kg m-3]
230 real, dimension(SZI_(G),SZJ_(G)), intent(out) :: h_bot !< Bottom boundary layer thickness [H ~> m or kg m-2]
231 integer, dimension(SZI_(G),SZJ_(G)), intent(out) :: k_bot !< Bottom boundary layer top layer index
232
233 ! Local variables
234 real, dimension(SZI_(G),SZK_(GV)+1) :: &
235 pres, & ! The pressure at each interface [R L2 T-2 ~> Pa].
236 dRho_int ! The unfiltered density differences across interfaces [R ~> kg m-3].
237 real, dimension(SZI_(G),SZK_(GV)) :: dz ! Layer thicknesses in depth units [Z ~> m]
238 real, dimension(SZI_(G)) :: &
239 Temp_int, & ! The temperature at each interface [C ~> degC]
240 Salin_int, & ! The salinity at each interface [S ~> ppt]
241 drho_bot, & ! The density difference at the bottom of a layer [R ~> kg m-3]
242 h_amp, & ! The amplitude of topographic roughness [Z ~> m].
243 hb, & ! The thickness of the water column below the midpoint of a layer [H ~> m or kg m-2]
244 z_from_bot, & ! The distance of a layer center from the bottom [Z ~> m]
245 dRho_dT, & ! The partial derivative of density with temperature [R C-1 ~> kg m-3 degC-1]
246 dRho_dS ! The partial derivative of density with salinity [R S-1 ~> kg m-3 ppt-1].
247
248 real :: dz_int ! The vertical extent of water associated with an interface [Z ~> m]
249 real :: G_Rho0 ! The gravitational acceleration, sometimes divided by the Boussinesq
250 ! density [H T-2 R-1 ~> m4 s-2 kg-1 or m s-2].
251 logical :: do_i(SZI_(G)), do_any
252 integer, dimension(2) :: EOSdom ! The i-computational domain for the equation of state
253 integer :: i, j, k, is, ie, js, je, nz
254
255 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
256 g_rho0 = gv%g_Earth_Z_T2 / gv%H_to_RZ
257 eosdom(:) = eos_domain(g%HI)
258
259 ! Find the (limited) density jump across each interface.
260 do i=is,ie
261 drho_int(i,1) = 0.0 ; drho_int(i,nz+1) = 0.0
262 enddo
263
264 !$OMP parallel do default(none) shared(is,ie,js,je,nz,tv,fluxes,G,GV,US,h,T_f,S_f, &
265 !$OMP h2,N2_bot,Rho_bot,h_bot,k_bot,G_Rho0,EOSdom) &
266 !$OMP private(pres,Temp_Int,Salin_Int,dRho_dT,dRho_dS, &
267 !$OMP dz,hb,dRho_bot,z_from_bot,do_i,h_amp,do_any,dz_int) &
268 !$OMP firstprivate(dRho_int)
269 do j=js,je
270
271 ! Find the vertical distances across layers.
272 call thickness_to_dz(h, tv, dz, j, g, gv)
273
274 if (associated(tv%eqn_of_state)) then
275 if (associated(fluxes%p_surf)) then
276 do i=is,ie ; pres(i,1) = fluxes%p_surf(i,j) ; enddo
277 else
278 do i=is,ie ; pres(i,1) = 0.0 ; enddo
279 endif
280 do k=2,nz
281 do i=is,ie
282 pres(i,k) = pres(i,k-1) + (gv%g_Earth*gv%H_to_RZ)*h(i,j,k-1)
283 temp_int(i) = 0.5 * (t_f(i,j,k) + t_f(i,j,k-1))
284 salin_int(i) = 0.5 * (s_f(i,j,k) + s_f(i,j,k-1))
285 enddo
286 call calculate_density_derivs(temp_int, salin_int, pres(:,k), drho_dt(:), drho_ds(:), &
287 tv%eqn_of_state, eosdom)
288 do i=is,ie
289 drho_int(i,k) = max(drho_dt(i)*(t_f(i,j,k) - t_f(i,j,k-1)) + &
290 drho_ds(i)*(s_f(i,j,k) - s_f(i,j,k-1)), 0.0)
291 enddo
292 enddo
293 else
294 do k=2,nz ; do i=is,ie
295 drho_int(i,k) = (gv%Rlay(k) - gv%Rlay(k-1))
296 enddo ; enddo
297 endif
298
299 ! Find the bottom boundary layer stratification.
300 do i=is,ie
301 hb(i) = 0.0 ; drho_bot(i) = 0.0
302 z_from_bot(i) = 0.5*dz(i,nz)
303 do_i(i) = (g%mask2dT(i,j) > 0.0)
304 h_amp(i) = sqrt(h2(i,j))
305 enddo
306
307 do k=nz,2,-1
308 do_any = .false.
309 do i=is,ie ; if (do_i(i)) then
310 dz_int = 0.5*(dz(i,k) + dz(i,k-1))
311 z_from_bot(i) = z_from_bot(i) + dz_int ! middle of the layer above
312
313 hb(i) = hb(i) + 0.5*(h(i,j,k) + h(i,j,k-1))
314 drho_bot(i) = drho_bot(i) + drho_int(i,k)
315
316 if (z_from_bot(i) > h_amp(i)) then
317 if (k>2) then
318 ! Always include at least one full layer.
319 hb(i) = hb(i) + 0.5*(h(i,j,k-1) + h(i,j,k-2))
320 drho_bot(i) = drho_bot(i) + drho_int(i,k-1)
321 endif
322 do_i(i) = .false.
323 else
324 do_any = .true.
325 endif
326 endif ; enddo
327 if (.not.do_any) exit
328 enddo
329
330 do i=is,ie
331 if (hb(i) > 0.0) then
332 n2_bot(i,j) = (g_rho0 * drho_bot(i)) / hb(i)
333 else ; n2_bot(i,j) = 0.0 ; endif
334 enddo
335
336 if (gv%Boussinesq .or. gv%semi_Boussinesq) then
337 do i=is,ie
338 rho_bot(i,j) = gv%Rho0
339 enddo
340 else
341 ! Average the density over the envelope of the topography.
342 call find_rho_bottom(g, gv, us, tv, h, dz, pres, h_amp, j, rho_bot(:,j), h_bot(:,j), k_bot(:,j))
343 endif
344 enddo
345
346end subroutine find_n2_bottom
347
348!> Returns TKE_itidal_input
349subroutine get_input_tke(G, TKE_itidal_input, nFreq, CS)
350 type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure (in).
351 integer, intent(in) :: nfreq !< number of frequencies
352 real, dimension(SZI_(G),SZJ_(G),nFreq), &
353 intent(out) :: tke_itidal_input !< The energy input to the internal waves
354 !! [H Z2 T-3 ~> m3 s-3 or W m-2].
355 type(int_tide_input_cs), target :: cs !< A pointer that is set to point to the control
356 !! structure for the internal tide input module.
357 integer :: i,j,fr
358
359 do fr=1,nfreq ; do j=g%jsd,g%jed ; do i=g%isd,g%ied
360 tke_itidal_input(i,j,fr) = cs%TKE_itidal_input(i,j,fr)
361 enddo ; enddo ; enddo
362
363end subroutine get_input_tke
364
365!> Returns barotropic tidal velocities
366subroutine get_barotropic_tidal_vel(G, vel_btTide, nFreq, CS)
367 type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure (in).
368 integer, intent(in) :: nfreq !< number of frequencies
369 real, dimension(SZI_(G),SZJ_(G),nFreq), &
370 intent(out) :: vel_bttide !< Barotropic velocity read from file [L T-1 ~> m s-1].
371 type(int_tide_input_cs), target :: cs !< A pointer that is set to point to the control
372 !! structure for the internal tide input module.
373 integer :: i,j,fr
374
375 do fr=1,nfreq ; do j=g%jsd,g%jed ; do i=g%isd,g%ied
376 vel_bttide(i,j,fr) = cs%tideamp(i,j,fr)
377 enddo ; enddo ; enddo
378
379end subroutine get_barotropic_tidal_vel
380
381!> Initializes the data related to the internal tide input module
382subroutine int_tide_input_init(Time, G, GV, US, param_file, diag, CS, itide)
383 type(time_type), intent(in) :: time !< The current model time
384 type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure
385 type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure
386 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
387 type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters
388 type(diag_ctrl), target, intent(inout) :: diag !< structure used to regulate diagnostic output.
389 type(int_tide_input_cs), pointer :: cs !< This module's control structure, which is initialized here.
390 type(int_tide_input_type), pointer :: itide !< A structure containing fields related
391 !! to the internal tide sources.
392 ! Local variables
393 logical :: read_tideamp
394 ! This include declares and sets the variable "version".
395# include "version_variable.h"
396 character(len=40) :: mdl = "MOM_int_tide_input" ! This module's name.
397 character(len=200) :: filename, tideamp_file, h2_file ! Input file names or paths
398 character(len=80) :: tideamp_var, rough_var ! Input file variable names
399 character(len=80) :: var_name
400 character(len=200) :: var_descript
401 character(len=200) :: tidefile_varnames
402
403 real :: mask_itidal ! A multiplicative land mask, 0 or 1 [nondim]
404 real :: max_frac_rough ! The fraction relating the maximum topographic roughness
405 ! to the mean depth [nondim]
406 real :: utide ! constant tidal amplitude [L T-1 ~> m s-1] to be used if
407 ! tidal amplitude file is not present.
408 real :: kappa_h2_factor ! factor for the product of wavenumber * rms sgs height [nondim].
409 real :: kappa_itides ! topographic wavenumber and non-dimensional scaling [L-1 ~> m-1]
410 real :: min_zbot_itides ! Minimum ocean depth for internal tide conversion [Z ~> m].
411 real :: hz2_t3_to_w_m2 ! unit conversion factor for TKE from internal units
412 ! to mks [T3 kg H-1 Z-2 s-3 ~> kg m-3 or 1]
413 real :: w_m2_to_hz2_t3 ! unit conversion factor for TKE from mks to internal
414 ! units [H Z2 s3 T-3 kg-1 ~> m3 kg-1 or 1]
415 integer :: tlen_days !< Time interval from start for adding wave source
416 !! for testing internal tides (BDM)
417 integer :: i, j, is, ie, js, je, isd, ied, jsd, jed
418 integer :: num_freq, fr
419
420 if (associated(cs)) then
421 call mom_error(warning, "int_tide_input_init called with an associated "// &
422 "control structure.")
423 return
424 endif
425 if (associated(itide)) then
426 call mom_error(warning, "int_tide_input_init called with an associated "// &
427 "internal tide input type.")
428 return
429 endif
430 allocate(cs)
431 allocate(itide)
432
433 cs%initialized = .true.
434
435 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
436 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
437
438 cs%diag => diag
439
440 hz2_t3_to_w_m2 = gv%H_to_kg_m2*(us%Z_to_m**2)*(us%s_to_T**3)
441 w_m2_to_hz2_t3 = gv%kg_m2_to_H*(us%m_to_Z**2)*(us%T_to_s**3)
442
443 ! Read all relevant parameters and write them to the model log.
444 call log_version(param_file, mdl, version, "")
445
446 call get_param(param_file, mdl, "INPUTDIR", cs%inputdir, default=".")
447 cs%inputdir = slasher(cs%inputdir)
448
449 call get_param(param_file, mdl, "DEBUG", cs%debug, default=.false., do_not_log=.true.)
450
451 call get_param(param_file, mdl, "MIN_ZBOT_ITIDES", min_zbot_itides, &
452 "Turn off internal tidal dissipation when the total "//&
453 "ocean depth is less than this value.", units="m", default=0.0, scale=us%m_to_Z)
454 call get_param(param_file, mdl, "KD_SMOOTH", cs%kappa_fill, &
455 "A diapycnal diffusivity that is used to interpolate "//&
456 "more sensible values of T & S into thin layers.", &
457 units="m2 s-1", default=1.0e-6, scale=gv%m2_s_to_HZ_T)
458
459 call get_param(param_file, mdl, "UTIDE", utide, &
460 "The constant tidal amplitude used with INT_TIDE_DISSIPATION.", &
461 units="m s-1", default=0.0, scale=us%m_s_to_L_T)
462
463 call read_param(param_file, "INTERNAL_TIDE_FREQS", num_freq)
464 cs%nFreq= num_freq
465
466 allocate(itide%Nb(isd:ied,jsd:jed), source=0.0)
467 allocate(itide%Rho_bot(isd:ied,jsd:jed), source=0.0)
468 allocate(itide%h2(isd:ied,jsd:jed), source=0.0)
469 allocate(cs%TKE_itidal_input(isd:ied,jsd:jed,num_freq), source=0.0)
470 allocate(cs%tideamp(isd:ied,jsd:jed,num_freq), source=utide)
471 allocate(cs%TKE_itidal_coef(isd:ied,jsd:jed, num_freq), source=0.0)
472
473 call get_param(param_file, mdl, "KAPPA_ITIDES", kappa_itides, &
474 "A topographic wavenumber used with INT_TIDE_DISSIPATION. "//&
475 "The default is 2pi/10 km, as in St.Laurent et al. 2002.", &
476 units="m-1", default=8.e-4*atan(1.0), scale=us%L_to_m)
477
478 call get_param(param_file, mdl, "KAPPA_H2_FACTOR", kappa_h2_factor, &
479 "A scaling factor for the roughness amplitude with "//&
480 "INT_TIDE_DISSIPATION.", units="nondim", default=1.0)
481 call get_param(param_file, mdl, "TKE_ITIDE_MAX", cs%TKE_itide_maxi, &
482 "The maximum internal tide energy source available to mix "//&
483 "above the bottom boundary layer with INT_TIDE_DISSIPATION.", &
484 units="W m-2", default=1.0e3, scale=w_m2_to_hz2_t3)
485
486 call get_param(param_file, mdl, "READ_TIDEAMP", read_tideamp, &
487 "If true, read a file (given by TIDEAMP_FILE) containing "//&
488 "the tidal amplitude with INT_TIDE_DISSIPATION.", default=.false.)
489 if (read_tideamp) then
490 call get_param(param_file, mdl, "TIDEAMP_FILE", tideamp_file, &
491 "The path to the file containing the spatially varying "//&
492 "tidal amplitudes with INT_TIDE_DISSIPATION.", default="tideamp.nc")
493 filename = trim(cs%inputdir) // trim(tideamp_file)
494 call log_param(param_file, mdl, "INPUTDIR/TIDEAMP_FILE", filename)
495
496 call read_param(param_file, "INTTIDE_AMP_VARNAMES", tidefile_varnames)
497 do fr=1,num_freq
498 tideamp_var = extractword(tidefile_varnames,fr)
499 call mom_read_data(filename, tideamp_var, cs%tideamp(:,:,fr), g%domain, scale=us%m_s_to_L_T)
500 enddo
501
502 endif
503
504 call get_param(param_file, mdl, "H2_FILE", h2_file, &
505 "The path to the file containing the sub-grid-scale "//&
506 "topographic roughness amplitude with INT_TIDE_DISSIPATION.", &
507 fail_if_missing=.true.)
508 filename = trim(cs%inputdir) // trim(h2_file)
509 call log_param(param_file, mdl, "INPUTDIR/H2_FILE", filename)
510 call get_param(param_file, mdl, "ROUGHNESS_VARNAME", rough_var, &
511 "The name in the input file of the squared sub-grid-scale "//&
512 "topographic roughness amplitude variable.", default="h2")
513 call mom_read_data(filename, rough_var, itide%h2, g%domain, scale=us%m_to_Z**2)
514
515 call get_param(param_file, mdl, "FRACTIONAL_ROUGHNESS_MAX", max_frac_rough, &
516 "The maximum topographic roughness amplitude as a fraction of the mean depth, "//&
517 "or a negative value for no limitations on roughness.", &
518 units="nondim", default=0.1)
519
520 ! The following parameters are used in testing the internal tide code.
521 call get_param(param_file, mdl, "INTERNAL_TIDE_SOURCE_TEST", cs%int_tide_source_test, &
522 "If true, apply an arbitrary generation site for internal tide testing", &
523 default=.false.)
524 if (cs%int_tide_source_test)then
525 call get_param(param_file, mdl, "INTERNAL_TIDE_USE_GLOB_IJ", cs%int_tide_use_glob_ij, &
526 "Use global IJ for internal tide generation source test", default=.false.)
527 call get_param(param_file, mdl, "INTERNAL_TIDE_SOURCE_X", cs%int_tide_source_x, &
528 "X Location of generation site for internal tide", &
529 units=g%x_ax_unit_short, default=1.0, do_not_log=cs%int_tide_use_glob_ij)
530 call get_param(param_file, mdl, "INTERNAL_TIDE_SOURCE_Y", cs%int_tide_source_y, &
531 "Y Location of generation site for internal tide", &
532 units=g%y_ax_unit_short, default=1.0, do_not_log=cs%int_tide_use_glob_ij)
533 call get_param(param_file, mdl, "INTERNAL_TIDE_SOURCE_I", cs%int_tide_source_i, &
534 "I Location of generation site for internal tide", default=0, &
535 do_not_log=.not.cs%int_tide_use_glob_ij)
536 call get_param(param_file, mdl, "INTERNAL_TIDE_SOURCE_J", cs%int_tide_source_j, &
537 "J Location of generation site for internal tide", default=0, &
538 do_not_log=.not.cs%int_tide_use_glob_ij)
539 call get_param(param_file, mdl, "INTERNAL_TIDE_SOURCE_TLEN_DAYS", tlen_days, &
540 "Time interval from start of experiment for adding wave source", &
541 units="days", default=0)
542 cs%time_max_source = time + set_time(0, days=tlen_days)
543
544 if ((cs%int_tide_use_glob_ij) .and. ((cs%int_tide_source_x /= 1.) .or. (cs%int_tide_source_y /= 1.))) then
545 call mom_error(fatal, "MOM_internal_tide_input: "//&
546 "Internal tide source set to use (i,j) indices hence (x,y) geographical coords are meaningless.")
547 endif
548 if ((.not.cs%int_tide_use_glob_ij) .and. ((cs%int_tide_source_i /= 0) .or. (cs%int_tide_source_j /= 0))) then
549 call mom_error(fatal, "MOM_internal_tide_input: "//&
550 "Internal tide source set to use (x,y) geographical coords hence (i,j) indices are meaningless.")
551 endif
552 endif
553
554 do fr=1,num_freq ; do j=js,je ; do i=is,ie
555 mask_itidal = 1.0
556 if (g%meanSL(i,j) + g%bathyT(i,j) < min_zbot_itides) mask_itidal = 0.0
557
558 cs%tideamp(i,j,fr) = cs%tideamp(i,j,fr) * mask_itidal * g%mask2dT(i,j)
559
560 ! Restrict rms topo to a fraction (often 10 percent) of the column depth.
561 if (max_frac_rough >= 0.0) &
562 itide%h2(i,j) = min((max_frac_rough * max(g%meanSL(i,j) + g%bathyT(i,j), 0.0))**2, itide%h2(i,j))
563
564 ! Compute the fixed part of internal tidal forcing; units are [R Z4 H-1 T-2 ~> J m-2 or J m kg-1] here.
565 cs%TKE_itidal_coef(i,j,fr) = 0.5*us%L_to_Z*kappa_h2_factor * gv%H_to_RZ * &
566 kappa_itides * itide%h2(i,j) * cs%tideamp(i,j,fr)**2
567 enddo ; enddo ; enddo
568
569
570 allocate( cs%id_TKE_itidal_itide(num_freq), source=-1)
571
572 do fr=1,num_freq
573 write(var_name, '("TKE_itidal_itide_freq",i1)') fr
574 write(var_descript, '("Internal Tide Driven Turbulent Kinetic Energy in frequency ",i1)') fr
575
576 cs%id_TKE_itidal_itide(fr) = register_diag_field('ocean_model',var_name,diag%axesT1,time, &
577 var_descript, 'W m-2', conversion=hz2_t3_to_w_m2)
578 enddo
579
580 cs%id_Nb = register_diag_field('ocean_model','Nb_itide',diag%axesT1,time, &
581 'Bottom Buoyancy Frequency', 's-1', conversion=us%s_to_T)
582
583 cs%id_N2_bot = register_diag_field('ocean_model','N2_b_itide',diag%axesT1,time, &
584 'Bottom Buoyancy frequency squared', 's-2', conversion=us%s_to_T**2)
585
586end subroutine int_tide_input_init
587
588!> Deallocates any memory related to the internal tide input module.
589subroutine int_tide_input_end(CS)
590 type(int_tide_input_cs), pointer :: cs !< This module's control structure, which is deallocated here.
591
592 if (associated(cs)) deallocate(cs)
593
594end subroutine int_tide_input_end
595
596end module mom_int_tide_input