MOM_internal_tides.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!> Subroutines that use the ray-tracing equations to propagate the internal tide energy density.
6!!
7!! \author Benjamin Mater & Robert Hallberg, 2015
9
10use mom_checksums, only : hchksum
11use mom_debugging, only : is_nan
12use mom_diag_mediator, only : post_data, query_averaging_enabled, diag_axis_init
14use mom_diag_mediator, only : register_diag_field, diag_ctrl, safe_alloc_ptr
16use mom_domains, only : agrid, to_south, to_west, to_all, cgrid_ne
17use mom_domains, only : create_group_pass, do_group_pass, pass_var, pass_vector
18use mom_domains, only : group_pass_type, start_group_pass, complete_group_pass
19use mom_error_handler, only : mom_error, fatal, warning, mom_mesg, is_root_pe
20use mom_file_parser, only : read_param, get_param, log_param, log_version, param_file_type
22use mom_grid, only : ocean_grid_type
24use mom_io, only : slasher, mom_read_data, file_exists, axis_info
25use mom_io, only : set_axis_info, get_axis_info, stdout
26use mom_restart, only : register_restart_field, mom_restart_cs, restart_init, save_restart
30use mom_time_manager, only : time_type, time_type_to_real, operator(+), operator(/), operator(-)
35use mpp_domains_mod, only : north_face => north, east_face => east
36
37implicit none ; private
38
39#include <MOM_memory.h>
40
44
45!> This control structure has parameters for the MOM_internal_tides module
46type, public :: int_tide_cs ; private
47 logical :: initialized = .false. !< True if this control structure has been initialized.
48 logical :: do_int_tides !< If true, use the internal tide code.
49 integer :: nfreq = 0 !< The number of internal tide frequency bands
50 integer :: nmode = 1 !< The number of internal tide vertical modes
51 integer :: nangle = 24 !< The number of internal tide angular orientations
52 integer :: energized_angle = -1 !< If positive, only this angular band is energized for debugging purposes
53 real :: dt_itides !< The timestep for internal tides ray-tracing [T ~> s]
54 real :: uniform_test_cg !< Uniform group velocity of internal tide
55 !! for testing internal tides [L T-1 ~> m s-1]
56 logical :: corner_adv !< If true, use a corner advection rather than PPM.
57 logical :: upwind_1st !< If true, use a first-order upwind scheme.
58 logical :: simple_2nd !< If true, use a simple second order (arithmetic mean) interpolation
59 !! of the edge values instead of the higher order interpolation.
60 logical :: vol_cfl !< If true, use the ratio of the open face lengths to the tracer cell
61 !! areas when estimating CFL numbers. Without aggress_adjust,
62 !! the default is false; it is always true with aggress_adjust.
63 logical :: use_ppmang !< If true, use PPM for advection of energy in angular space.
64 logical :: update_kd !< If true, the scheme will modify the diffusivities seen by the dynamics
65 logical :: apply_refraction !< If false, skip refraction (for debugging)
66 logical :: apply_propagation !< If False, do not propagate energy (for debugging)
67 logical :: turn_critical_lat !< If True, rays change direction at critical latitude instead
68 !! of being trapped
69 logical :: reflect_critical_lat !< If True, rays reflect at the critical latitude instead
70 !! of turning parallel to it
71 logical :: debug !< If true, use debugging prints
72 logical :: init_forcing_only !< if True, add TKE forcing only at first step (for debugging)
73 logical :: force_posit_en !< if True, remove subroundoff negative values (needs enhancement)
74 logical :: add_tke_forcing = .true. !< Whether to add forcing, used by init_forcing_only
75
76 real, allocatable, dimension(:,:) :: fraction_tidal_input
77 !< how the energy from one tidal component is distributed
78 !! over the various vertical modes, 2d in frequency and mode [nondim]
79 real, allocatable, dimension(:,:) :: refl_angle
80 !< local coastline/ridge/shelf angles read from file [rad]
81 ! (could be in G control structure)
82 real :: nullangle = -999.9 !< placeholder value in cells with no reflection [rad]
83 real, allocatable, dimension(:,:) :: refl_pref
84 !< partial reflection coeff for each "coast cell" [nondim]
85 ! (could be in G control structure)
86 logical, allocatable, dimension(:,:) :: refl_pref_logical
87 !< true if reflecting cell with partial reflection
88 ! (could be in G control structure)
89 logical, allocatable, dimension(:,:) :: refl_dbl
90 !< identifies reflection cells where double reflection
91 !! is possible (i.e. ridge cells)
92 ! (could be in G control structure)
93 real, allocatable, dimension(:,:) :: trans
94 !< partial transmission coeff for each "coast cell" [nondim]
95 real, allocatable, dimension(:,:) :: residual
96 !< residual of reflection and transmission coeff for each "coast cell" [nondim]
97 real, allocatable, dimension(:,:,:,:) :: cp
98 !< horizontal phase speed [L T-1 ~> m s-1]
99 real, allocatable, dimension(:,:,:,:,:) :: tke_leak_loss
100 !< energy lost due to misc background processes [H Z2 T-3 ~> m3 s-3 or W m-2]
101 real, allocatable, dimension(:,:,:,:,:) :: tke_quad_loss
102 !< energy lost due to quadratic bottom drag [H Z2 T-3 ~> m3 s-3 or W m-2]
103 real, allocatable, dimension(:,:,:,:,:) :: tke_froude_loss
104 !< energy lost due to wave breaking [H Z2 T-3 ~> m3 s-3 or W m-2]
105 real, allocatable, dimension(:,:) :: tke_itidal_loss_fixed
106 !< Fixed part of the energy lost due to small-scale drag [H Z2 L-2 ~> kg m-2] here.
107 !! This will be multiplied by N and the squared near-bottom velocity (and by
108 !! the near-bottom density in non-Boussinesq mode) to get the energy losses
109 !! in [R Z4 H-1 L-2 ~> kg m-2 or m]
110 real, allocatable, dimension(:,:,:,:,:) :: tke_itidal_loss
111 !< energy lost due to small-scale wave drag [H Z2 T-3 ~> m3 s-3 or W m-2]
112 real, allocatable, dimension(:,:,:,:,:) :: tke_residual_loss
113 !< internal tide energy loss due to the residual at slopes [H Z2 T-3 ~> m3 s-3 or W m-2]
114 real, allocatable, dimension(:,:,:,:,:) :: tke_slope_loss
115 !< internal tide energy loss due to the residual at slopes [H Z2 T-3 ~> m3 s-3 or W m-2]
116 real, allocatable, dimension(:,:) :: tke_input_glo_dt
117 !< The integrated energy input to the internal waves [H Z2 L2 T-2 ~> m5 s-2 or J]
118 real, allocatable, dimension(:,:) :: tke_leak_loss_glo_dt
119 !< Integrated energy lost due to misc background processes [H Z2 L2 T-2 ~> m5 s-2 or J]
120 real, allocatable, dimension(:,:) :: tke_quad_loss_glo_dt
121 !< Integrated energy lost due to quadratic bottom drag [H Z2 L2 T-2 ~> m5 s-2 or J]
122 real, allocatable, dimension(:,:) :: tke_froude_loss_glo_dt
123 !< Integrated energy lost due to wave breaking [H Z2 L2 T-2 ~> m5 s-2 or J]
124 real, allocatable, dimension(:,:) :: tke_itidal_loss_glo_dt
125 !< energy lost due to small-scale wave drag [H Z2 T-2 ~> m3 s-2 or J m-2]
126 real, allocatable, dimension(:,:) :: tke_residual_loss_glo_dt
127 !< internal tide energy loss due to the residual at slopes [H Z2 L2 T-2 ~> m5 s-2 or J]
128 real, allocatable, dimension(:,:) :: error_mode
129 !< internal tide energy budget error for each mode [H Z2 L2 T-2 ~> m5 s-2 or J]
130 real, allocatable, dimension(:,:) :: tot_leak_loss !< Energy loss rates due to misc background processes,
131 !! summed over angle, frequency and mode [H Z2 T-3 ~> m3 s-3 or W m-2]
132 real, allocatable, dimension(:,:) :: tot_quad_loss !< Energy loss rates due to quadratic bottom drag,
133 !! summed over angle, frequency and mode [H Z2 T-3 ~> m3 s-3 or W m-2]
134 real, allocatable, dimension(:,:) :: tot_itidal_loss !< Energy loss rates due to small-scale drag,
135 !! summed over angle, frequency and mode [H Z2 T-3 ~> m3 s-3 or W m-2]
136 real, allocatable, dimension(:,:) :: tot_froude_loss !< Energy loss rates due to wave breaking,
137 !! summed over angle, frequency and mode [H Z2 T-3 ~> m3 s-3 or W m-2]
138 real, allocatable, dimension(:,:) :: tot_residual_loss !< Energy loss rates due to residual on slopes,
139 !! summed over angle, frequency and mode [H Z2 T-3 ~> m3 s-3 or W m-2]
140 real, allocatable, dimension(:,:) :: tot_allprocesses_loss !< Energy loss rates due to all processes,
141 !! summed over angle, frequency and mode [H Z2 T-3 ~> m3 s-3 or W m-2]
142 real, allocatable, dimension(:,:,:,:) :: w_struct !< Vertical structure of vertical velocity (normalized)
143 !! for each frequency and each mode [nondim]
144 real, allocatable, dimension(:,:,:,:) :: u_struct !< Vertical structure of horizontal velocity (normalized and
145 !! divided by layer thicknesses) for each frequency and each mode [Z-1 ~> m-1]
146 real, allocatable, dimension(:,:,:) :: u_struct_max !< Maximum of u_struct,
147 !! for each mode [Z-1 ~> m-1]
148 real, allocatable, dimension(:,:,:) :: u_struct_bot !< Bottom value of u_struct,
149 !! for each mode [Z-1 ~> m-1]
150 real, allocatable, dimension(:,:,:) :: int_w2 !< Vertical integral of w_struct squared,
151 !! for each mode [H ~> m or kg m-2]
152 real, allocatable, dimension(:,:,:) :: int_u2 !< Vertical integral of u_struct squared,
153 !! for each mode [H Z-2 ~> m-1 or kg m-4]
154 real, allocatable, dimension(:,:,:) :: int_n2w2 !< Depth-integrated Brunt Vaissalla freqency times
155 !! vertical profile squared, for each mode [H T-2 ~> m s-2 or kg m-2 s-2]
156 real :: q_itides !< fraction of local dissipation [nondim]
157 real :: mixing_effic !< mixing efficiency [nondim]
158 real :: en_sum !< global sum of energy for use in debugging, in MKS units [m5 s-2 or J]
159 real :: en_underflow !< A minuscule amount of energy [H Z2 T-2 ~> m3 s-2 or J m-2]
160 integer :: en_restart_power !< A power factor of 2 by which to multiply the energy in restart [nondim]
161 type(time_type), pointer :: time => null() !< A pointer to the model's clock.
162 type(group_pass_type) :: pass_en !< Pass 5d array Energy as a group of 3d arrays
163 character(len=200) :: inputdir !< directory to look for coastline angle file
164 integer :: itides_adv_limiter !< The type of limiter to use for the energy advection scheme
165 real, allocatable, dimension(:,:,:,:) :: decay_rate_2d !< rate at which internal tide energy is
166 !! lost to the interior ocean internal wave field
167 !! as a function of longitude, latitude, frequency
168 !! and vertical mode [T-1 ~> s-1].
169 real :: cdrag !< The bottom drag coefficient [nondim].
170 real :: drag_min_depth !< The minimum total ocean thickness that will be used in the denominator
171 !! of the quadratic drag terms for internal tides when
172 !! INTERNAL_TIDE_QUAD_DRAG is true [H ~> m or kg m-2]
173 real :: gamma_osborn !< Mixing efficiency from Osborn 1980 [nondim]
174 real :: kd_min !< The minimum diapycnal diffusivity. [L2 T-1 ~> m2 s-1]
175 real :: max_tke_to_kd !< Maximum allowed value for TKE_to_kd [H Z2 T-3 ~> m3 s-3 or W m-2]
176 real :: min_thick_layer_kd !< minimum layer thickness allowed to use with TKE_to_kd [H ~> m or kg m-2]
177 logical :: apply_background_drag
178 !< If true, apply a drag due to background processes as a sink.
179 logical :: apply_bottom_drag
180 !< If true, apply a quadratic bottom drag as a sink.
181 logical :: apply_wave_drag
182 !< If true, apply scattering due to small-scale roughness as a sink.
183 logical :: apply_froude_drag
184 !< If true, apply wave breaking as a sink.
185 real :: en_check_tol !< An energy density tolerance for flagging points with small negative
186 !! internal tide energy [H Z2 T-2 ~> m3 s-2 or J m-2]
187 logical :: apply_residual_drag
188 !< If true, apply sink from residual term of reflection/transmission.
189 logical :: use_2d_decay_rate
190 !< If true, use a spatially varying decay rate for each harmonic.
191 real, allocatable :: en(:,:,:,:,:)
192 !< The internal wave energy density as a function of (i,j,angle,frequency,mode)
193 !! integrated within an angular and frequency band [H Z2 T-2 ~> m3 s-2 or J m-2]
194 real, allocatable :: en_ini_glo(:,:)
195 !< The internal wave energy density as a function of (frequency,mode) spatially
196 !! integrated within an angular and frequency band [H Z2 L2 T-2 ~> m5 s-2 or J]
197 !! only at the start of the routine (for diags)
198 real, allocatable :: en_end_glo(:,:)
199 !< The internal wave energy density as a function of (frequency,mode) spatially
200 !! integrated within an angular and frequency band [H Z2 L2 T-2 ~> m5 s-2 or J]
201 !! only at the end of the routine (for diags)
202 real, allocatable :: en_restart_mode1(:,:,:,:)
203 !< The internal wave energy density as a function of (i,j,angle,freq)
204 !! for mode 1 [H Z2 T-2 ~> m3 s-2 or J m-2]
205 real, allocatable :: en_restart_mode2(:,:,:,:)
206 !< The internal wave energy density as a function of (i,j,angle,freq)
207 !! for mode 2 [H Z2 T-2 ~> m3 s-2 or J m-2]
208 real, allocatable :: en_restart_mode3(:,:,:,:)
209 !< The internal wave energy density as a function of (i,j,angle,freq)
210 !! for mode 3 [H Z2 T-2 ~> m3 s-2 or J m-2]
211 real, allocatable :: en_restart_mode4(:,:,:,:)
212 !< The internal wave energy density as a function of (i,j,angle,freq)
213 !! for mode 4 [H Z2 T-2 ~> m3 s-2 or J m-2]
214 real, allocatable :: en_restart_mode5(:,:,:,:)
215 !< The internal wave energy density as a function of (i,j,angle,freq)
216 !! for mode 5 [H Z2 T-2 ~> m3 s-2 or J m-2]
217
218 real, allocatable, dimension(:) :: frequency !< The frequency of each band [T-1 ~> s-1].
219 real :: int_tide_decay_scale !< vertical decay scale for St Laurent profile [Z ~> m]
220 real :: int_tide_decay_scale_slope !< vertical decay scale for St Laurent profile on slopes [Z ~> m]
221
222 type(wave_speed_cs) :: wave_speed !< Wave speed control structure
223 type(diag_ctrl), pointer :: diag => null() !< A structure that is used to regulate the
224 !! timing of diagnostic output.
225
226 !>@{ Diag handles
227 ! Diag handles relevant to all modes, frequencies, and angles
228 integer :: id_cg1 = -1 ! diagnostic handle for mode-1 speed
229 integer, allocatable, dimension(:) :: id_cn ! diagnostic handle for all mode speeds
230 integer :: id_tot_en = -1
231 integer :: id_refl_pref = -1, id_refl_ang = -1, id_land_mask = -1
232 integer :: id_trans = -1, id_residual = -1
233 integer :: id_dx_cv = -1, id_dy_cu = -1
234 ! Diag handles considering: sums over all modes, frequencies, and angles
235 integer :: id_tot_leak_loss = -1, id_tot_quad_loss = -1, id_tot_itidal_loss = -1
236 integer :: id_tot_froude_loss = -1, id_tot_residual_loss = -1, id_tot_allprocesses_loss = -1
237 ! Diag handles considering: all modes & frequencies; summed over angles
238 integer, allocatable, dimension(:,:) :: &
239 id_en_mode, &
240 id_itidal_loss_mode, &
241 id_leak_loss_mode, &
242 id_quad_loss_mode, &
243 id_froude_loss_mode, &
244 id_residual_loss_mode, &
245 id_allprocesses_loss_mode, &
246 id_itide_drag, &
247 id_ub_mode, &
248 id_cp_mode
249 ! Diag handles considering: all modes, frequencies, and angles
250 integer, allocatable, dimension(:,:) :: &
251 id_en_ang_mode, &
252 id_itidal_loss_ang_mode
253 integer, allocatable, dimension(:) :: &
254 id_tke_itidal_input, &
255 id_ustruct_mode, &
256 id_wstruct_mode, &
257 id_int_w2_mode, &
258 id_int_u2_mode, &
259 id_int_n2w2_mode
260 !>@}
261
262end type int_tide_cs
263
264!> A structure with the active energy loop bounds.
265type :: loop_bounds_type ; private
266 !>@{ The active loop bounds
267 integer :: ish, ieh, jsh, jeh
268 !>@}
269end type loop_bounds_type
270
271!>@{ Enumeration values for numerical schemes
272integer, parameter :: limiter_adv_minmod = 1
273integer, parameter :: limiter_adv_positive = 2
274character*(20), parameter :: limiter_adv_minmod_string = "MINMOD"
275character*(20), parameter :: limiter_adv_positive_string = "POSITIVE"
276!>@}
277
278contains
279
280!> Calls subroutines in this file that are needed to refract, propagate,
281!! and dissipate energy density of the internal tide.
282subroutine propagate_int_tide(h, tv, Nb, Rho_bot, dt, G, GV, US, inttide_input_CSp, CS)
283 type(ocean_grid_type), intent(inout) :: g !< The ocean's grid structure.
284 type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure.
285 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
286 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
287 intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2]
288 type(thermo_var_ptrs), intent(in) :: tv !< Pointer to thermodynamic variables
289 !! (needed for wave structure).
290 real, dimension(SZI_(G),SZJ_(G)), intent(inout) :: nb !< Near-bottom buoyancy frequency [T-1 ~> s-1].
291 !! In some cases the input values are used, but in
292 !! others this is set along with the wave speeds.
293 real, dimension(SZI_(G),SZJ_(G)), intent(in) :: rho_bot !< Near-bottom density or the Boussinesq
294 !! reference density [R ~> kg m-3].
295 real, intent(in) :: dt !< Length of time over which to advance
296 !! the internal tides [T ~> s].
297 type(int_tide_input_cs), intent(in) :: inttide_input_csp !< Internal tide input control structure
298 type(int_tide_cs), intent(inout) :: cs !< Internal tide control structure
299
300 ! Local variables
301 real, dimension(SZI_(G),SZJ_(G),CS%nFreq) :: &
302 tke_itidal_input, & !< The energy input to the internal waves [H Z2 T-3 ~> m3 s-3 or W m-2].
303 vel_bttide !< Barotropic velocity read from file [L T-1 ~> m s-1].
304
305 real, dimension(SZI_(G),SZJ_(G),2) :: &
306 test ! A test unit vector used to determine grid rotation in halos [nondim]
307 real, dimension(SZI_(G),SZJ_(G),CS%nMode) :: &
308 cn ! baroclinic internal gravity wave speeds for each mode [L T-1 ~> m s-1]
309 real, dimension(SZI_(G),SZJ_(G),CS%nFreq,CS%nMode) :: &
310 tot_en_mode, & ! energy summed over angles only [H Z2 T-2 ~> m3 s-2 or J m-2]
311 ub, & ! near-bottom horizontal velocity of wave (modal) [L T-1 ~> m s-1]
312 umax ! Maximum horizontal velocity of wave (modal) [L T-1 ~> m s-1]
313 real, dimension(SZI_(G),SZJ_(G),CS%nFreq,CS%nMode) :: &
314 drag_scale ! bottom drag scale [T-1 ~> s-1]
315 real, dimension(SZI_(G),SZJ_(G)) :: &
316 tot_vel_bttide2, & ! [L2 T-2 ~> m2 s-2]
317 tot_en, & ! energy summed over angles, modes, frequencies [H Z2 T-2 ~> m3 s-2 or J m-2]
318 tot_leak_loss, tot_quad_loss, tot_itidal_loss, tot_froude_loss, tot_residual_loss, tot_allprocesses_loss, &
319 ! energy loss rates summed over angle, freq, and mode [H Z2 T-3 ~> m3 s-3 or W m-2]
320 htot, & ! The vertical sum of the layer thicknesses [H ~> m or kg m-2]
321 itidal_loss_mode, & ! Energy lost due to small-scale wave drag, summed over angles [H Z2 T-3 ~> m3 s-3 or W m-2]
322 leak_loss_mode, &
323 quad_loss_mode, &
324 froude_loss_mode, &
325 residual_loss_mode, &
326 allprocesses_loss_mode ! Total energy loss rates for a given mode and frequency (summed over
327 ! all angles) [H Z2 T-3 ~> m3 s-3 or W m-2]
328 real :: frac_per_sector ! The inverse of the number of angular, modal and frequency bins [nondim]
329 real :: f2 ! The squared Coriolis parameter interpolated to a tracer point [T-2 ~> s-2]
330 real :: kmag2 ! A squared horizontal wavenumber [L-2 ~> m-2]
331 real :: i_d_here ! The inverse of the local water column thickness [H-1 ~> m-1 or m2 kg-1]
332 real :: i_mass ! The inverse of the local water mass [R-1 Z-1 ~> m2 kg-1]
333 real :: i_dt ! The inverse of the timestep [T-1 ~> s-1]
334 real :: dt_sub ! The effective timestep use to subcycle the propagation [T ~> s]
335 real :: en_restart_factor ! A multiplicative factor of the form 2**En_restart_power [nondim]
336 real :: i_en_restart_factor ! The inverse of the restart mult factor [nondim]
337 real :: freq2 ! The frequency squared [T-2 ~> s-2]
338 real :: pe_term ! total potential energy of profile [R Z ~> kg m-2]
339 real :: ke_term ! total kinetic energy of profile [R Z ~> kg m-2]
340 real :: u_mag ! rescaled magnitude of horizontal profile [L Z T-1 ~> m2 s-1]
341 real :: w0 ! rescaled magnitude of vertical profile [Z T-1 ~> m s-1]
342 real :: c_phase ! The phase speed [L T-1 ~> m s-1]
343 ! real :: loss_rate ! An energy loss rate [T-1 ~> s-1]
344 real :: fr2_max ! The column maximum internal wave Froude number squared [nondim]
345 real :: cn_subro ! A tiny wave speed to prevent division by zero [L T-1 ~> m s-1]
346 real :: en_subro ! A tiny energy to prevent division by zero [H Z2 T-2 ~> m3 s-2 or J m-2]
347 real :: en_a, en_b ! Energies for time stepping [H Z2 T-2 ~> m3 s-2 or J m-2]
348 real :: en_sumtmp ! Energies for debugging [H Z2 L2 T-2 ~> m5 s-2 or J]
349 real :: hz2_t2_to_j_m2 ! unit conversion factor for Energy from internal units
350 ! to mks [T2 kg H-1 Z-2 s-2 ~> kg m-3 or 1]
351 real :: j_m2_to_hz2_t2 ! unit conversion factor for Energy from mks to internal
352 ! units [H Z2 s2 T-2 kg-1 ~> m3 kg-1 or 1]
353 character(len=160) :: mesg ! The text of an error message
354 integer :: en_halo_ij_stencil ! The halo size needed for energy advection
355 integer :: a, m, fr, i, j, k, is, ie, js, je, isd, ied, jsd, jed, nangle, nc
356 integer :: id_g, jd_g ! global (decomp-invar) indices (for debugging)
357 integer :: subcycles ! number of subcycles for the propagation
358 type(group_pass_type), save :: pass_test
359 type(time_type) :: time_end
360 logical:: avg_enabled
361
362 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
363 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed ; nangle = cs%NAngle
364
365 hz2_t2_to_j_m2 = gv%H_to_kg_m2*(us%Z_to_m**2)*(us%s_to_T**2)
366 j_m2_to_hz2_t2 = gv%kg_m2_to_H*(us%m_to_Z**2)*(us%T_to_s**2)
367
368 cn_subro = 1e-30*us%m_s_to_L_T
369 en_subro = 1e-30*j_m2_to_hz2_t2
370
371 i_dt = 1.0 / dt
372 en_restart_factor = 2**cs%En_restart_power
373 i_en_restart_factor = 1.0 / en_restart_factor
374
375 if (cs%dt_itides <= 0.) then
376 subcycles = 1
377 else
378 subcycles = ceiling(dt/cs%dt_itides - 0.0001)
379 endif
380 dt_sub = dt / subcycles
381
382 ! initialize local arrays
383 tke_itidal_input(:,:,:) = 0.
384 vel_bttide(:,:,:) = 0.
385 tot_vel_bttide2(:,:) = 0.
386 drag_scale(:,:,:,:) = 0.
387 ub(:,:,:,:) = 0.
388 umax(:,:,:,:) = 0.
389
390 cn(:,:,:) = 0.
391
392 ! Rebuild energy density array from multiple restarts
393 do fr=1,cs%nFreq ; do a=1,cs%nAngle ; do j=jsd,jed ; do i=isd,ied
394 cs%En(i,j,a,fr,1) = cs%En_restart_mode1(i,j,a,fr) * i_en_restart_factor
395 enddo ; enddo ; enddo ; enddo
396
397 if (cs%nMode >= 2) then
398 do fr=1,cs%nFreq ; do a=1,cs%nAngle ; do j=jsd,jed ; do i=isd,ied
399 cs%En(i,j,a,fr,2) = cs%En_restart_mode2(i,j,a,fr) * i_en_restart_factor
400 enddo ; enddo ; enddo ; enddo
401 endif
402
403 if (cs%nMode >= 3) then
404 do fr=1,cs%nFreq ; do a=1,cs%nAngle ; do j=jsd,jed ; do i=isd,ied
405 cs%En(i,j,a,fr,3) = cs%En_restart_mode3(i,j,a,fr) * i_en_restart_factor
406 enddo ; enddo ; enddo ; enddo
407 endif
408
409 if (cs%nMode >= 4) then
410 do fr=1,cs%nFreq ; do a=1,cs%nAngle ; do j=jsd,jed ; do i=isd,ied
411 cs%En(i,j,a,fr,4) = cs%En_restart_mode4(i,j,a,fr) * i_en_restart_factor
412 enddo ; enddo ; enddo ; enddo
413 endif
414
415 if (cs%nMode >= 5) then
416 do fr=1,cs%nFreq ; do a=1,cs%nAngle ; do j=jsd,jed ; do i=isd,ied
417 cs%En(i,j,a,fr,5) = cs%En_restart_mode5(i,j,a,fr) * i_en_restart_factor
418 enddo ; enddo ; enddo ; enddo
419 endif
420
421 if (cs%debug) then
422 ! save initial energy for online budget
423 do m=1,cs%nMode ; do fr=1,cs%nFreq
424 en_sumtmp = 0.
425 do a=1,cs%nAngle
426 en_sumtmp = en_sumtmp + global_area_integral(cs%En(:,:,a,fr,m), g, tmp_scale=hz2_t2_to_j_m2)
427 enddo
428 cs%En_ini_glo(fr,m) = en_sumtmp
429 enddo ; enddo
430 endif
431
432 ! Set properties related to the internal tides, such as the wave speeds, storing some
433 ! of them in the control structure for this module.
434 if (cs%uniform_test_cg > 0.0) then
435 do m=1,cs%nMode ; cn(:,:,m) = cs%uniform_test_cg ; enddo
436 else
437 call wave_speeds(h, tv, g, gv, us, cs%nMode, cn, cs%wave_speed, &
438 cs%w_struct, cs%u_struct, cs%u_struct_max, cs%u_struct_bot, &
439 nb, cs%int_w2, cs%int_U2, cs%int_N2w2, halo_size=2)
440 ! The value of halo_size above would have to be larger if there were
441 ! not a halo update between the calls to propagate_x and propagate_y.
442 ! It can be 1 point smaller if teleport is not used.
443 endif
444
445 call pass_var(cn,g%Domain)
446
447 if (cs%debug) then
448 call hchksum(cn(:,:,1), "CN mode 1", g%HI, haloshift=0, unscale=us%L_to_m*us%s_to_T)
449 call hchksum(cs%w_struct(:,:,:,1), "Wstruct mode 1", g%HI, haloshift=0)
450 call hchksum(cs%u_struct(:,:,:,1), "Ustruct mode 1", g%HI, haloshift=0, unscale=us%m_to_Z)
451 call hchksum(cs%u_struct_bot(:,:,1), "Ustruct_bot mode 1", g%HI, haloshift=0, unscale=us%m_to_Z)
452 call hchksum(cs%u_struct_max(:,:,1), "Ustruct_max mode 1", g%HI, haloshift=0, unscale=us%m_to_Z)
453 call hchksum(cs%int_w2(:,:,1), "int_w2", g%HI, haloshift=0, unscale=gv%H_to_MKS)
454 call hchksum(cs%int_U2(:,:,1), "int_U2", g%HI, haloshift=0, unscale=gv%H_to_mks*us%m_to_Z**2)
455 call hchksum(cs%int_N2w2(:,:,1), "int_N2w2", g%HI, haloshift=0, unscale=gv%H_to_mks*us%s_to_T**2)
456 endif
457
458 ! Set the wave speeds for the modes, using cg(n) ~ cg(1)/n.**********************
459 ! This is wrong, of course, but it works reasonably in some cases.
460 ! Uncomment if wave_speed is not used to calculate the true values (BDM).
461 !do m=1,CS%nMode ; do j=js-2,je+2 ; do i=is-2,ie+2
462 ! cn(i,j,m) = cn(i,j,1) / real(m)
463 !enddo ; enddo ; enddo
464
465 ! Add the forcing.***************************************************************
466
467 if (cs%add_tke_forcing) then
468
469 call get_input_tke(g, tke_itidal_input, cs%nFreq, inttide_input_csp)
470
471 if (cs%debug) then
472 call hchksum(tke_itidal_input(:,:,1), "TKE_itidal_input", g%HI, haloshift=0, &
473 unscale=gv%H_to_mks*(us%Z_to_m**2)*(us%s_to_T)**3)
474 call hchksum(cs%En(:,:,:,1,1), "EnergyIntTides bf input", g%HI, haloshift=0, unscale=hz2_t2_to_j_m2)
475 endif
476
477 if (cs%energized_angle <= 0) then
478 frac_per_sector = 1.0 / real(cs%nAngle)
479 do m=1,cs%nMode ; do fr=1,cs%nFreq ; do a=1,cs%nAngle ; do j=js,je ; do i=is,ie
480 f2 = 0.25*((g%Coriolis2Bu(i,j) + g%Coriolis2Bu(i-1,j-1)) + &
481 (g%Coriolis2Bu(i-1,j) + g%Coriolis2Bu(i,j-1)))
482 if (cs%frequency(fr)**2 > f2) then
483 cs%En(i,j,a,fr,m) = cs%En(i,j,a,fr,m) + (dt*frac_per_sector*(1.0-cs%q_itides) * &
484 cs%fraction_tidal_input(fr,m) * tke_itidal_input(i,j,fr))
485 else
486 ! zero out input TKE value to get correct diagnostics
487 tke_itidal_input(i,j,fr) = 0.
488 endif
489 enddo ; enddo ; enddo ; enddo ; enddo
490 elseif (cs%energized_angle <= cs%nAngle) then
491 frac_per_sector = 1.0
492 a = cs%energized_angle
493 do m=1,cs%nMode ; do fr=1,cs%nFreq ; do j=js,je ; do i=is,ie
494 f2 = 0.25*((g%Coriolis2Bu(i,j) + g%Coriolis2Bu(i-1,j-1)) + &
495 (g%Coriolis2Bu(i-1,j) + g%Coriolis2Bu(i,j-1)))
496 if (cs%frequency(fr)**2 > f2) then
497 cs%En(i,j,a,fr,m) = cs%En(i,j,a,fr,m) + (dt*frac_per_sector*(1.0-cs%q_itides) * &
498 cs%fraction_tidal_input(fr,m) * tke_itidal_input(i,j,fr))
499 else
500 ! zero out input TKE value to get correct diagnostics
501 tke_itidal_input(i,j,fr) = 0.
502 endif
503 enddo ; enddo ; enddo ; enddo
504 else
505 call mom_error(warning, "Internal tide energy is being put into a angular "//&
506 "band that does not exist.")
507 endif
508 endif ! add tke forcing
509
510 if (cs%init_forcing_only) cs%add_tke_forcing=.false.
511
512 if (cs%debug) then
513 call hchksum(cs%En(:,:,:,1,1), "EnergyIntTides af input", g%HI, haloshift=0, unscale=hz2_t2_to_j_m2)
514 ! save forcing for online budget
515 do m=1,cs%nMode ; do fr=1,cs%nFreq
516 en_sumtmp = 0.
517 do a=1,cs%nAngle
518 en_sumtmp = en_sumtmp + global_area_integral(dt*frac_per_sector*(1.0-cs%q_itides)* &
519 cs%fraction_tidal_input(fr,m)*tke_itidal_input(:,:,fr), &
520 g, tmp_scale=hz2_t2_to_j_m2)
521 enddo
522 cs%TKE_input_glo_dt(fr,m) = en_sumtmp
523 enddo ; enddo
524 endif
525
526 ! Pass a test vector to check for grid rotation in the halo updates.
527 do j=jsd,jed ; do i=isd,ied ; test(i,j,1) = 0.0 ; test(i,j,2) = 1.0 ; enddo ; enddo
528 call create_group_pass(pass_test, test(:,:,1), test(:,:,2), g%domain, stagger=agrid)
529 call start_group_pass(pass_test, g%domain)
530
531 if (cs%debug) then
532 call hchksum(cs%En(:,:,:,1,1), "EnergyIntTides af halo", g%HI, haloshift=0, unscale=hz2_t2_to_j_m2)
533 do m=1,cs%nMode ; do fr=1,cs%Nfreq
534 call sum_en(g, gv, us, cs, cs%En(:,:,:,fr,m), 'prop_int_tide: after forcing')
535 if (is_root_pe()) write(stdout,'(A,E18.10)') 'prop_int_tide: after forcing', cs%En_sum
536 enddo ; enddo
537 endif
538
539 call complete_group_pass(pass_test, g%domain)
540
541 ! TKE_slope_loss need to be accumulated but since it is
542 ! passed as inout and accumulated within propagate_x/propagate_y
543 ! it does not need temp array for accumulation
544 cs%TKE_slope_loss(:,:,:,:,:) = 0.
545
546 ! Start subcycling
547 do nc=1,subcycles
548
549 ! Apply half the refraction.
550 if (cs%apply_refraction) then
551 do m=1,cs%nMode ; do fr=1,cs%nFreq
552 call refract(cs%En(:,:,:,fr,m), cn(:,:,m), cs%frequency(fr), 0.5*dt_sub, &
553 g, us, cs%nAngle, cs%use_PPMang)
554 enddo ; enddo
555 endif
556 ! A this point, CS%En is only valid on the computational domain.
557
558 if (cs%force_posit_En) then
559 do m=1,cs%nMode ; do fr=1,cs%Nfreq ; do a=1,cs%nAngle
560 do j=jsd,jed ; do i=isd,ied
561 if (cs%En(i,j,a,fr,m)<0.0) then
562 cs%En(i,j,a,fr,m) = 0.0
563 endif
564 enddo ; enddo
565 enddo ; enddo ; enddo
566 endif
567
568 if (cs%debug) then
569 call hchksum(cs%En(:,:,:,1,1), "EnergyIntTides af refr", g%HI, haloshift=0, unscale=hz2_t2_to_j_m2)
570 do m=1,cs%nMode ; do fr=1,cs%Nfreq
571 call sum_en(g, gv, us, cs, cs%En(:,:,:,fr,m), 'prop_int_tide: after 1/2 refraction')
572 if (is_root_pe()) write(stdout,'(A,E18.10)') 'prop_int_tide: after 1/2 refraction', cs%En_sum
573 enddo ; enddo
574 ! Check for En<0 - for debugging, delete later
575 do m=1,cs%nMode ; do fr=1,cs%Nfreq ; do a=1,cs%nAngle
576 do j=js,je ; do i=is,ie
577 if (cs%En(i,j,a,fr,m)<0.0) then
578 id_g = i + g%idg_offset ; jd_g = j + g%jdg_offset ! for debugging
579 write(mesg,*) 'After first refraction: En<0.0 at ig=', id_g, ', jg=', jd_g, &
580 'En=', hz2_t2_to_j_m2*cs%En(i,j,a,fr,m)
581 call mom_error(warning, "propagate_int_tide: "//trim(mesg))
582 !call MOM_error(FATAL, "propagate_int_tide: stopped due to negative energy.")
583 endif
584 enddo ; enddo
585 enddo ; enddo ; enddo
586 endif
587
588 ! Set the halo size to work on, using similar logic to that used in propagate. This may need
589 ! to be adjusted depending on the advection scheme and whether teleport is used.
590 if (cs%upwind_1st) then ; en_halo_ij_stencil = 2
591 else ; en_halo_ij_stencil = 3 ; endif
592
593 ! Rotate points in the halos as necessary.
594 call do_group_pass(cs%pass_En, g%domain)
595 call correct_halo_rotation(cs%En, test, g, cs%nAngle, halo=en_halo_ij_stencil)
596
597 if (cs%debug) then
598 call hchksum(cs%En(:,:,:,1,1), "EnergyIntTides af halo R", g%HI, haloshift=0, unscale=hz2_t2_to_j_m2)
599 do m=1,cs%nMode ; do fr=1,cs%Nfreq
600 call sum_en(g, gv, us, cs, cs%En(:,:,:,fr,m), 'prop_int_tide: after correct halo rotation')
601 if (is_root_pe()) write(stdout,'(A,E18.10)') 'prop_int_tide: after correct halo rotation', cs%En_sum
602 enddo ; enddo
603 endif
604
605 ! Propagate the waves.
606 do m=1,cs%nMode ; do fr=1,cs%Nfreq
607
608 if (cs%apply_propagation) then
609 call propagate(cs%En(:,:,:,fr,m), cn(:,:,m), cs%frequency(fr), dt_sub, &
610 g, gv, us, cs, cs%NAngle, test(:,:,:), en_halo_ij_stencil, cs%TKE_slope_loss(:,:,:,fr,m))
611 endif
612 enddo ; enddo
613
614 ! Rotate points in the halos as necessary.
615 call do_group_pass(cs%pass_En, g%domain)
616 call correct_halo_rotation(cs%En, test, g, cs%nAngle, halo=en_halo_ij_stencil)
617
618
619 if (cs%force_posit_En) then
620 do m=1,cs%nMode ; do fr=1,cs%Nfreq ; do a=1,cs%nAngle
621 do j=jsd,jed ; do i=isd,ied
622 if (cs%En(i,j,a,fr,m)<0.0) then
623 cs%En(i,j,a,fr,m) = 0.0
624 endif
625 enddo ; enddo
626 enddo ; enddo ; enddo
627 endif
628
629 if (cs%debug) then
630 call hchksum(cs%En(:,:,:,1,1), "EnergyIntTides af prop", g%HI, haloshift=0, unscale=hz2_t2_to_j_m2)
631 do m=1,cs%nMode ; do fr=1,cs%Nfreq
632 call sum_en(g, gv, us, cs, cs%En(:,:,:,fr,m), 'prop_int_tide: after propagate')
633 if (is_root_pe()) write(stdout,'(A,E18.10)') 'prop_int_tide: after propagate', cs%En_sum
634 enddo ; enddo
635 ! Check for En<0 - for debugging, delete later
636 do m=1,cs%nMode ; do fr=1,cs%Nfreq ; do a=1,cs%nAngle
637 do j=js,je ; do i=is,ie
638 if (cs%En(i,j,a,fr,m)<0.0) then
639 id_g = i + g%idg_offset ; jd_g = j + g%jdg_offset
640 if (abs(cs%En(i,j,a,fr,m))>cs%En_check_tol) then ! only print if large
641 write(mesg,*) 'After propagation: En<0.0 at ig=', id_g, ', jg=', jd_g, &
642 'En=', hz2_t2_to_j_m2*cs%En(i,j,a,fr,m)
643 call mom_error(warning, "propagate_int_tide: "//trim(mesg))
644 ! RD propagate produces very little negative energy (diff 2 large numbers), needs fix
645 !call MOM_error(FATAL, "propagate_int_tide: stopped due to negative energy.")
646 endif
647 endif
648 enddo ; enddo
649 enddo ; enddo ; enddo
650 endif
651
652 if (cs%apply_refraction) then
653 ! Apply the other half of the refraction.
654 do m=1,cs%nMode ; do fr=1,cs%Nfreq
655 call refract(cs%En(:,:,:,fr,m), cn(:,:,m), cs%frequency(fr), 0.5*dt_sub, &
656 g, us, cs%NAngle, cs%use_PPMang)
657 enddo ; enddo
658 ! A this point, CS%En is only valid on the computational domain.
659 endif
660
661 if (cs%force_posit_En) then
662 do m=1,cs%nMode ; do fr=1,cs%Nfreq ; do a=1,cs%nAngle
663 do j=jsd,jed ; do i=isd,ied
664 if (cs%En(i,j,a,fr,m)<0.0) then
665 cs%En(i,j,a,fr,m) = 0.0
666 endif
667 enddo ; enddo
668 enddo ; enddo ; enddo
669 endif
670
671 if (cs%debug) then
672 call hchksum(cs%En(:,:,:,1,1), "EnergyIntTides af refr2", g%HI, haloshift=0, unscale=hz2_t2_to_j_m2)
673 do m=1,cs%nMode ; do fr=1,cs%Nfreq
674 call sum_en(g, gv, us, cs, cs%En(:,:,:,fr,m), 'prop_int_tide: after 2/2 refraction')
675 if (is_root_pe()) write(stdout,'(A,E18.10)') 'prop_int_tide: after 2/2 refraction', cs%En_sum
676 enddo ; enddo
677 ! Check for En<0 - for debugging, delete later
678 do m=1,cs%nMode ; do fr=1,cs%Nfreq ; do a=1,cs%nAngle
679 do j=js,je ; do i=is,ie
680 if (cs%En(i,j,a,fr,m)<0.0) then
681 id_g = i + g%idg_offset ; jd_g = j + g%jdg_offset ! for debugging
682 write(mesg,*) 'After second refraction: En<0.0 at ig=', id_g, ', jg=', jd_g, &
683 'En=', hz2_t2_to_j_m2*cs%En(i,j,a,fr,m)
684 call mom_error(warning, "propagate_int_tide: "//trim(mesg))
685 !call MOM_error(FATAL, "propagate_int_tide: stopped due to negative energy.")
686 endif
687 enddo ; enddo
688 enddo ; enddo ; enddo
689 endif
690
691 call do_group_pass(cs%pass_En, g%domain)
692 call correct_halo_rotation(cs%En, test, g, cs%nAngle, halo=en_halo_ij_stencil)
693
694 enddo ! end subcycling
695
696 ! Apply various dissipation mechanisms.
697 if (cs%apply_background_drag .or. cs%apply_bottom_drag &
698 .or. cs%apply_wave_drag .or. cs%apply_Froude_drag &
699 .or. (cs%id_tot_En > 0)) then
700 tot_en(:,:) = 0.0
701 tot_en_mode(:,:,:,:) = 0.0
702 do m=1,cs%nMode ; do fr=1,cs%Nfreq
703 do j=js,je ; do i=is,ie ; do a=1,cs%nAngle
704 tot_en(i,j) = tot_en(i,j) + cs%En(i,j,a,fr,m)
705 tot_en_mode(i,j,fr,m) = tot_en_mode(i,j,fr,m) + cs%En(i,j,a,fr,m)
706 enddo ; enddo ; enddo
707 enddo ; enddo
708 endif
709
710 ! Extract the energy for mixing due to misc. processes (background leakage)------
711 if (cs%apply_background_drag) then
712 do m=1,cs%nMode ; do fr=1,cs%nFreq ; do a=1,cs%nAngle ; do j=js,je ; do i=is,ie
713 ! Calculate loss rate and apply loss over the time step ; apply the same drag timescale
714 ! to each En component (technically not correct; fix later)
715 en_b = cs%En(i,j,a,fr,m) ! save previous value
716 en_a = cs%En(i,j,a,fr,m) / (1.0 + (dt * cs%decay_rate_2d(i,j,fr,m))) ! implicit update
717 cs%TKE_leak_loss(i,j,a,fr,m) = (en_b - en_a) * i_dt ! compute exact loss rate [H Z2 T-3 ~> m3 s-3 or W m-2]
718 cs%En(i,j,a,fr,m) = en_a ! update value
719 enddo ; enddo ; enddo ; enddo ; enddo
720 endif
721
722 if (cs%force_posit_En) then
723 do m=1,cs%nMode ; do fr=1,cs%Nfreq ; do a=1,cs%nAngle
724 do j=jsd,jed ; do i=isd,ied
725 if (cs%En(i,j,a,fr,m)<0.0) then
726 cs%En(i,j,a,fr,m) = 0.0
727 endif
728 enddo ; enddo
729 enddo ; enddo ; enddo
730 endif
731
732 if (cs%debug) then
733 call hchksum(cs%En(:,:,:,1,1), "EnergyIntTides after leak", g%HI, haloshift=0, unscale=hz2_t2_to_j_m2)
734 do m=1,cs%nMode ; do fr=1,cs%Nfreq
735 call sum_en(g, gv, us, cs, cs%En(:,:,:,fr,m), 'prop_int_tide: after background drag')
736 if (is_root_pe()) write(stdout,'(A,E18.10)') 'prop_int_tide: after background drag', cs%En_sum
737 call sum_en(g, gv, us, cs, cs%TKE_leak_loss(:,:,:,fr,m) * dt, 'prop_int_tide: loss after background drag')
738 if (is_root_pe()) write(stdout,'(A,E18.10)') 'prop_int_tide: loss after background drag', cs%En_sum
739 enddo ; enddo
740 ! Check for En<0 - for debugging, delete later
741 do m=1,cs%nMode ; do fr=1,cs%Nfreq ; do a=1,cs%nAngle
742 do j=js,je ; do i=is,ie
743 if (cs%En(i,j,a,fr,m)<0.0) then
744 id_g = i + g%idg_offset ; jd_g = j + g%jdg_offset ! for debugging
745 write(mesg,*) 'After leak loss: En<0.0 at ig=', id_g, ', jg=', jd_g, &
746 'En=', hz2_t2_to_j_m2*cs%En(i,j,a,fr,m)
747 call mom_error(warning, "propagate_int_tide: "//trim(mesg), all_print=.true.)
748 !call MOM_error(FATAL, "propagate_int_tide: stopped due to negative energy.")
749 endif
750 enddo ; enddo
751 enddo ; enddo ; enddo
752 ! save loss term for online budget
753 do m=1,cs%nMode ; do fr=1,cs%nFreq
754 en_sumtmp = 0.
755 do a=1,cs%nAngle
756 en_sumtmp = en_sumtmp + global_area_integral(cs%TKE_leak_loss(:,:,a,fr,m)*dt, g, &
757 tmp_scale=hz2_t2_to_j_m2)
758 enddo
759 cs%TKE_leak_loss_glo_dt(fr,m) = en_sumtmp
760 enddo ; enddo
761 endif
762
763 ! Extract the energy for mixing due to bottom drag-------------------------------
764 if (cs%apply_bottom_drag) then
765 do j=jsd,jed ; do i=isd,ied ; htot(i,j) = 0.0 ; enddo ; enddo
766
767 call get_barotropic_tidal_vel(g, vel_bttide, cs%nFreq, inttide_input_csp)
768
769 do fr=1,cs%Nfreq ; do j=jsd,jed ; do i=isd,ied
770 tot_vel_bttide2(i,j) = tot_vel_bttide2(i,j) + (vel_bttide(i,j,fr)**2)
771 enddo ; enddo ; enddo
772
773 do k=1,gv%ke ; do j=jsd,jed ; do i=isd,ied
774 htot(i,j) = htot(i,j) + h(i,j,k)
775 enddo ; enddo ; enddo
776 if (gv%Boussinesq) then
777 ! This is mathematically equivalent to the form in the option below, but they differ at roundoff.
778 do m=1,cs%NMode ; do fr=1,cs%Nfreq ; do j=jsd,jed ; do i=isd,ied
779 i_d_here = 1.0 / (max(htot(i,j), cs%drag_min_depth))
780 drag_scale(i,j,fr,m) = cs%cdrag * sqrt(max(0.0, us%L_to_Z**2*tot_vel_bttide2(i,j) + &
781 (tot_en_mode(i,j,fr,m) * i_d_here))) * gv%Z_to_H*i_d_here
782 enddo ; enddo ; enddo ; enddo
783 else
784 do m=1,cs%NMode ; do fr=1,cs%Nfreq ; do j=jsd,jed ; do i=isd,ied
785 i_d_here = 1.0 / (max(htot(i,j), cs%drag_min_depth))
786 i_mass = gv%RZ_to_H * i_d_here
787 drag_scale(i,j,fr,m) = (cs%cdrag * (rho_bot(i,j)*i_mass)) * &
788 sqrt(max(0.0, us%L_to_Z**2*tot_vel_bttide2(i,j) + &
789 (tot_en_mode(i,j,fr,m) * i_d_here)))
790 enddo ; enddo ; enddo ; enddo
791 endif
792
793 if (cs%debug) call hchksum(drag_scale(:,:,1,1), "dragscale", g%HI, haloshift=0, unscale=us%s_to_T)
794 if (cs%debug) call hchksum(tot_vel_bttide2(:,:), "tot_vel_btTide2", g%HI, haloshift=0, unscale=us%L_T_to_m_s**2)
795
796 do m=1,cs%nMode ; do fr=1,cs%nFreq ; do a=1,cs%nAngle ; do j=js,je ; do i=is,ie
797 ! Calculate loss rate and apply loss over the time step ; apply the same drag timescale
798 ! to each En component (technically not correct; fix later)
799 en_b = cs%En(i,j,a,fr,m)
800 en_a = cs%En(i,j,a,fr,m) / (1.0 + (dt * drag_scale(i,j,fr,m))) ! implicit update
801 cs%TKE_quad_loss(i,j,a,fr,m) = (en_b - en_a) * i_dt
802 cs%En(i,j,a,fr,m) = en_a
803 enddo ; enddo ; enddo ; enddo ; enddo
804 endif
805
806 if (cs%force_posit_En) then
807 do m=1,cs%nMode ; do fr=1,cs%Nfreq ; do a=1,cs%nAngle
808 do j=jsd,jed ; do i=isd,ied
809 if (cs%En(i,j,a,fr,m)<0.0) then
810 cs%En(i,j,a,fr,m) = 0.0
811 endif
812 enddo ; enddo
813 enddo ; enddo ; enddo
814 endif
815
816 if (cs%debug) then
817 call hchksum(cs%En(:,:,:,1,1), "EnergyIntTides after quad", g%HI, haloshift=0, unscale=hz2_t2_to_j_m2)
818 ! save loss term for online budget
819 do m=1,cs%nMode ; do fr=1,cs%nFreq
820 en_sumtmp = 0.
821 do a=1,cs%nAngle
822 en_sumtmp = en_sumtmp + global_area_integral(cs%TKE_quad_loss(:,:,a,fr,m)*dt, g, &
823 tmp_scale=hz2_t2_to_j_m2)
824 enddo
825 cs%TKE_quad_loss_glo_dt(fr,m) = en_sumtmp
826 enddo ; enddo
827 ! Check for En<0 - for debugging, delete later
828 do m=1,cs%nMode ; do fr=1,cs%Nfreq ; do a=1,cs%nAngle
829 do j=js,je ; do i=is,ie
830 if (cs%En(i,j,a,fr,m)<0.0) then
831 id_g = i + g%idg_offset ; jd_g = j + g%jdg_offset ! for debugging
832 write(mesg,*) 'After bottom loss: En<0.0 at ig=', id_g, ', jg=', jd_g, &
833 'En=', hz2_t2_to_j_m2*cs%En(i,j,a,fr,m)
834 call mom_error(warning, "propagate_int_tide: "//trim(mesg), all_print=.true.)
835 !call MOM_error(FATAL, "propagate_int_tide: stopped due to negative energy.")
836 endif
837 enddo ; enddo
838 enddo ; enddo ; enddo
839 endif
840
841 ! Extract the energy for mixing due to scattering (wave-drag)--------------------
842 ! still need to allow a portion of the extracted energy to go to higher modes.
843 ! First, find velocity profiles
844 if (cs%apply_wave_drag .or. cs%apply_Froude_drag) then
845 do m=1,cs%nMode ; do fr=1,cs%Nfreq
846
847 ! compute near-bottom and max horizontal baroclinic velocity values at each point
848 do j=js,je ; do i=is,ie
849 id_g = i + g%idg_offset ; jd_g = j + g%jdg_offset ! for debugging
850
851 ! Calculate wavenumber magnitude
852 freq2 = cs%frequency(fr)**2
853
854 f2 = (0.25*(g%CoriolisBu(i,j) + g%CoriolisBu(max(i-1,1),max(j-1,1)) + &
855 g%CoriolisBu(i,max(j-1,1)) + g%CoriolisBu(max(i-1,1),j)))**2
856 kmag2 = (freq2 - f2) / ((cn(i,j,m)**2) + (cn_subro**2))
857
858
859 ! Back-calculate amplitude from energy equation
860 if ( (g%mask2dT(i,j) > 0.5) .and. (freq2*kmag2 > 0.0)) then
861 ! Units here are [R Z ~> kg m-2]
862 ke_term = 0.25*gv%H_to_RZ*( (((freq2 + f2) / (freq2*kmag2))*us%L_to_Z**2*cs%int_U2(i,j,m)) + &
863 cs%int_w2(i,j,m) )
864 pe_term = 0.25*gv%H_to_RZ*( cs%int_N2w2(i,j,m) / freq2 )
865
866 if (ke_term + pe_term > 0.0) then
867 w0 = sqrt( gv%H_to_RZ * tot_en_mode(i,j,fr,m) / (ke_term + pe_term) )
868 else
869 !call MOM_error(WARNING, "MOM internal tides: KE + PE <= 0.0; setting to W0 to 0.0")
870 w0 = 0.0
871 endif
872
873 u_mag = w0 * sqrt((freq2 + f2) / (2.0*freq2*kmag2))
874 ! scaled maximum tidal velocity
875 umax(i,j,fr,m) = abs(u_mag * cs%u_struct_max(i,j,m))
876 ! scaled bottom tidal velocity
877 ub(i,j,fr,m) = abs(u_mag * cs%u_struct_bot(i,j,m))
878 else
879 umax(i,j,fr,m) = 0.
880 ub(i,j,fr,m) = 0.
881 endif
882
883 enddo ; enddo ! i-loop, j-loop
884 enddo ; enddo ! fr-loop, m-loop
885 endif ! apply_wave or _Froude_drag (Ub or Umax needed)
886 ! Finally, apply loss
887 if (cs%apply_wave_drag) then
888 ! Calculate loss rate and apply loss over the time step
889 call itidal_lowmode_loss(g, gv, us, cs, nb, rho_bot, ub, cs%En, cs%TKE_itidal_loss_fixed, &
890 cs%TKE_itidal_loss, dt, halo_size=0)
891 endif
892
893 if (cs%force_posit_En) then
894 do m=1,cs%nMode ; do fr=1,cs%Nfreq ; do a=1,cs%nAngle
895 do j=jsd,jed ; do i=isd,ied
896 if (cs%En(i,j,a,fr,m)<0.0) then
897 cs%En(i,j,a,fr,m) = 0.0
898 endif
899 enddo ; enddo
900 enddo ; enddo ; enddo
901 endif
902
903 if (cs%debug) then
904 call hchksum(cs%En(:,:,:,1,1), "EnergyIntTides after wave", g%HI, haloshift=0, unscale=hz2_t2_to_j_m2)
905 do m=1,cs%nMode ; do fr=1,cs%Nfreq
906 call sum_en(g, gv, us, cs, cs%En(:,:,:,fr,m), 'prop_int_tide: before Froude drag')
907 if (is_root_pe()) write(stdout,'(A,E18.10)') 'prop_int_tide: before Froude drag', cs%En_sum
908 enddo ; enddo
909 ! save loss term for online budget, may want to add a debug flag later
910 do m=1,cs%nMode ; do fr=1,cs%nFreq
911 en_sumtmp = 0.
912 do a=1,cs%nAngle
913 en_sumtmp = en_sumtmp + global_area_integral(cs%TKE_itidal_loss(:,:,a,fr,m)*dt, g, &
914 tmp_scale=hz2_t2_to_j_m2)
915 enddo
916 cs%TKE_itidal_loss_glo_dt(fr,m) = en_sumtmp
917 enddo ; enddo
918 ! Check for En<0 - for debugging, delete later
919 do m=1,cs%nMode ; do fr=1,cs%Nfreq ; do a=1,cs%nAngle
920 do j=js,je ; do i=is,ie
921 if (cs%En(i,j,a,fr,m)<0.0) then
922 id_g = i + g%idg_offset ; jd_g = j + g%jdg_offset ! for debugging
923 write(mesg,*) 'After wave drag loss: En<0.0 at ig=', id_g, ', jg=', jd_g, &
924 'En=', hz2_t2_to_j_m2*cs%En(i,j,a,fr,m)
925 call mom_error(warning, "propagate_int_tide: "//trim(mesg), all_print=.true.)
926 !call MOM_error(FATAL, "propagate_int_tide: stopped due to negative energy.")
927 endif
928 enddo ; enddo
929 enddo ; enddo ; enddo
930 endif
931
932 ! Extract the energy for mixing due to wave breaking-----------------------------
933 if (cs%apply_Froude_drag) then
934 ! Pick out maximum baroclinic velocity values; calculate Fr=max(u)/cg
935 do m=1,cs%nMode ; do fr=1,cs%Nfreq
936 freq2 = cs%frequency(fr)**2
937 do j=js,je ; do i=is,ie
938 id_g = i + g%idg_offset ; jd_g = j + g%jdg_offset ! for debugging
939 ! Calculate horizontal phase velocity magnitudes
940 f2 = 0.25*((g%Coriolis2Bu(i,j) + g%Coriolis2Bu(i-1,j-1)) + &
941 (g%Coriolis2Bu(i-1,j) + g%Coriolis2Bu(i,j-1)))
942 kmag2 = (freq2 - f2) / ((cn(i,j,m)**2) + (cn_subro**2))
943 c_phase = 0.0
944 cs%TKE_Froude_loss(i,j,:,fr,m) = 0. ! init for all angles
945 if (kmag2 > 0.0) then
946 c_phase = sqrt(freq2/kmag2)
947 fr2_max = (umax(i,j,fr,m) / c_phase)**2
948 ! Dissipate energy if Fr>1; done here with an arbitrary time scale
949 if (fr2_max > 1.0) then
950 ! Calculate effective decay rate [T-1 ~> s-1] if breaking occurs over a time step
951 !loss_rate = (1.0 - Fr2_max) / (Fr2_max * dt)
952 do a=1,cs%nAngle
953 ! Determine effective dissipation rate (Wm-2)
954 !CS%TKE_Froude_loss(i,j,a,fr,m) = CS%En(i,j,a,fr,m) * abs(loss_rate)
955 en_b = cs%En(i,j,a,fr,m)
956 en_a = cs%En(i,j,a,fr,m)/fr2_max
957 cs%TKE_Froude_loss(i,j,a,fr,m) = (en_b - en_a) * i_dt
958 cs%En(i,j,a,fr,m) = en_a
959 enddo
960 endif ! Fr2>1
961 endif ! Kmag2>0
962 cs%cp(i,j,fr,m) = c_phase
963 enddo ; enddo
964 enddo ; enddo
965 endif
966
967 if (cs%force_posit_En) then
968 do m=1,cs%nMode ; do fr=1,cs%Nfreq ; do a=1,cs%nAngle
969 do j=jsd,jed ; do i=isd,ied
970 if (cs%En(i,j,a,fr,m)<0.0) then
971 cs%En(i,j,a,fr,m) = 0.0
972 endif
973 enddo ; enddo
974 enddo ; enddo ; enddo
975 endif
976
977 if (cs%debug) then
978 call hchksum(cs%En(:,:,:,1,1), "EnergyIntTides after froude", g%HI, haloshift=0, unscale=hz2_t2_to_j_m2)
979 do m=1,cs%nMode ; do fr=1,cs%Nfreq
980 call sum_en(g, gv, us, cs, cs%En(:,:,:,fr,m), 'prop_int_tide: after Froude drag')
981 if (is_root_pe()) write(stdout,'(A,E18.10)') 'prop_int_tide: after Froude drag', cs%En_sum
982 call sum_en(g, gv, us, cs, cs%TKE_Froude_loss(:,:,:,fr,m) * dt, 'prop_int_tide: loss after Froude drag')
983 if (is_root_pe()) write(stdout,'(A,E18.10)') 'prop_int_tide: loss after Froude drag', cs%En_sum
984 enddo ; enddo
985 ! save loss term for online budget, may want to add a debug flag later
986 do m=1,cs%nMode ; do fr=1,cs%nFreq
987 en_sumtmp = 0.
988 do a=1,cs%nAngle
989 en_sumtmp = en_sumtmp + global_area_integral(cs%TKE_Froude_loss(:,:,a,fr,m)*dt, g, &
990 tmp_scale=hz2_t2_to_j_m2)
991 enddo
992 cs%TKE_Froude_loss_glo_dt(fr,m) = en_sumtmp
993 enddo ; enddo
994 ! Check for En<0 - for debugging, delete later
995 do m=1,cs%nMode ; do fr=1,cs%Nfreq ; do a=1,cs%nAngle
996 do j=js,je ; do i=is,ie
997 if (cs%En(i,j,a,fr,m)<0.0) then
998 id_g = i + g%idg_offset ; jd_g = j + g%jdg_offset
999 write(mesg,*) 'After Froude loss: En<0.0 at ig=', id_g, ', jg=', jd_g, &
1000 'En=', hz2_t2_to_j_m2*cs%En(i,j,a,fr,m)
1001 call mom_error(warning, "propagate_int_tide: "//trim(mesg), all_print=.true.)
1002 !call MOM_error(FATAL, "propagate_int_tide: stopped due to negative energy.")
1003 endif
1004 enddo ; enddo
1005 enddo ; enddo ; enddo
1006 endif
1007
1008 ! loss from residual of reflection/transmission coefficients
1009 if (cs%apply_residual_drag) then
1010
1011 do m=1,cs%nMode ; do fr=1,cs%nFreq ; do a=1,cs%nAngle ; do j=js,je ; do i=is,ie
1012 ! implicit form is rewritten to minimize number of divisions
1013 !CS%En(i,j,a,fr,m) = CS%En(i,j,a,fr,m) / (1.0 + dt * CS%TKE_residual_loss(i,j,a,fr,m) / &
1014 ! (CS%En(i,j,a,fr,m) + en_subRO))
1015 ! only compute when partial reflection is present not to create noise elsewhere
1016 if (cs%refl_pref_logical(i,j)) then
1017 en_b = cs%En(i,j,a,fr,m)
1018 en_a = (cs%En(i,j,a,fr,m) * (cs%En(i,j,a,fr,m) + en_subro)) / &
1019 ((cs%En(i,j,a,fr,m) + en_subro) + (dt * cs%TKE_slope_loss(i,j,a,fr,m)))
1020 cs%TKE_residual_loss(i,j,a,fr,m) = (en_b - en_a) * i_dt
1021 cs%En(i,j,a,fr,m) = en_a
1022 endif
1023 enddo ; enddo ; enddo ; enddo ; enddo
1024
1025 else
1026 ! zero out the residual loss term so it does not count towards diagnostics
1027 do m=1,cs%nMode ; do fr=1,cs%nFreq ; do a=1,cs%nAngle ; do j=js,je ; do i=is,ie
1028 cs%TKE_residual_loss(i,j,a,fr,m) = 0.
1029 enddo ; enddo ; enddo ; enddo ; enddo
1030 endif
1031
1032 if (cs%debug) then
1033 call hchksum(cs%En(:,:,:,1,1), "EnergyIntTides after slope", g%HI, haloshift=0, unscale=hz2_t2_to_j_m2)
1034 do m=1,cs%nMode ; do fr=1,cs%Nfreq
1035 call sum_en(g, gv, us, cs, cs%En(:,:,:,fr,m), 'prop_int_tide')
1036 enddo ; enddo
1037 ! save loss term for online budget
1038 do m=1,cs%nMode ; do fr=1,cs%nFreq
1039 en_sumtmp = 0.
1040 do a=1,cs%nAngle
1041 en_sumtmp = en_sumtmp + global_area_integral(cs%TKE_residual_loss(:,:,a,fr,m)*dt, g, &
1042 tmp_scale=hz2_t2_to_j_m2)
1043 enddo
1044 cs%TKE_residual_loss_glo_dt(fr,m) = en_sumtmp
1045 enddo ; enddo
1046 endif
1047
1048 !---- energy budget ----
1049 if (cs%debug) then
1050 ! save final energy for online budget
1051 do m=1,cs%nMode ; do fr=1,cs%nFreq
1052 en_sumtmp = 0.
1053 do a=1,cs%nAngle
1054 en_sumtmp = en_sumtmp + global_area_integral(cs%En(:,:,a,fr,m), g, tmp_scale=hz2_t2_to_j_m2)
1055 enddo
1056 cs%En_end_glo(fr,m) = en_sumtmp
1057 enddo ; enddo
1058
1059 do m=1,cs%nMode ; do fr=1,cs%nFreq
1060 cs%error_mode(fr,m) = cs%En_ini_glo(fr,m) + cs%TKE_input_glo_dt(fr,m) - cs%TKE_leak_loss_glo_dt(fr,m) - &
1061 cs%TKE_quad_loss_glo_dt(fr,m) - cs%TKE_itidal_loss_glo_dt(fr,m) - &
1062 cs%TKE_Froude_loss_glo_dt(fr,m) - cs%TKE_residual_loss_glo_dt(fr,m) - &
1063 cs%En_end_glo(fr,m)
1064 if (is_root_pe()) write(stdout,'(A,F18.10)') &
1065 "error in Energy budget", us%L_to_m**2*hz2_t2_to_j_m2*cs%error_mode(fr,m)
1066 enddo ; enddo
1067 endif
1068
1069 ! Output diagnostics.************************************************************
1070 avg_enabled = query_averaging_enabled(cs%diag, time_end=time_end)
1071 call enable_averages(dt, time_end, cs%diag)
1072
1073 if (query_averaging_enabled(cs%diag)) then
1074 ! Output internal wave modal wave speeds
1075 if (cs%id_cg1 > 0) call post_data(cs%id_cg1, cn(:,:,1),cs%diag)
1076 do m=1,cs%nMode ; if (cs%id_cn(m) > 0) call post_data(cs%id_cn(m), cn(:,:,m), cs%diag) ; enddo
1077
1078 ! Output two-dimensional diagnostics
1079 if (cs%id_tot_En > 0) call post_data(cs%id_tot_En, tot_en, cs%diag)
1080 do fr=1,cs%nFreq
1081 if (cs%id_TKE_itidal_input(fr) > 0) call post_data(cs%id_TKE_itidal_input(fr), &
1082 tke_itidal_input(:,:,fr), cs%diag)
1083 enddo
1084
1085 do m=1,cs%nMode ; do fr=1,cs%nFreq
1086 if (cs%id_itide_drag(fr,m) > 0) call post_data(cs%id_itide_drag(fr,m), drag_scale(:,:,fr,m), cs%diag)
1087 enddo ; enddo
1088
1089 ! Output 2-D energy density (summed over angles) for each frequency and mode
1090 do m=1,cs%nMode ; do fr=1,cs%Nfreq ; if (cs%id_En_mode(fr,m) > 0) then
1091 tot_en(:,:) = 0.0
1092 do a=1,cs%nAngle ; do j=js,je ; do i=is,ie
1093 tot_en(i,j) = tot_en(i,j) + cs%En(i,j,a,fr,m)
1094 enddo ; enddo ; enddo
1095 call post_data(cs%id_En_mode(fr,m), tot_en, cs%diag)
1096 endif ; enddo ; enddo
1097
1098 ! split energy array into multiple restarts
1099 do fr=1,cs%Nfreq ; do a=1,cs%nAngle ; do j=jsd,jed ; do i=isd,ied
1100 cs%En_restart_mode1(i,j,a,fr) = cs%En(i,j,a,fr,1) * en_restart_factor
1101 enddo ; enddo ; enddo ; enddo
1102
1103 if (cs%nMode >= 2) then
1104 do fr=1,cs%Nfreq ; do a=1,cs%nAngle ; do j=jsd,jed ; do i=isd,ied
1105 cs%En_restart_mode2(i,j,a,fr) = cs%En(i,j,a,fr,2) * en_restart_factor
1106 enddo ; enddo ; enddo ; enddo
1107 endif
1108
1109 if (cs%nMode >= 3) then
1110 do fr=1,cs%Nfreq ; do a=1,cs%nAngle ; do j=jsd,jed ; do i=isd,ied
1111 cs%En_restart_mode3(i,j,a,fr) = cs%En(i,j,a,fr,3) * en_restart_factor
1112 enddo ; enddo ; enddo ; enddo
1113 endif
1114
1115 if (cs%nMode >= 4) then
1116 do fr=1,cs%Nfreq ; do a=1,cs%nAngle ; do j=jsd,jed ; do i=isd,ied
1117 cs%En_restart_mode4(i,j,a,fr) = cs%En(i,j,a,fr,4) * en_restart_factor
1118 enddo ; enddo ; enddo ; enddo
1119 endif
1120
1121 if (cs%nMode >= 5) then
1122 do fr=1,cs%Nfreq ; do a=1,cs%nAngle ; do j=jsd,jed ; do i=isd,ied
1123 cs%En_restart_mode5(i,j,a,fr) = cs%En(i,j,a,fr,5) * en_restart_factor
1124 enddo ; enddo ; enddo ; enddo
1125 endif
1126
1127 ! Output 3-D (i,j,a) energy density for each frequency and mode
1128 do m=1,cs%nMode ; do fr=1,cs%Nfreq ; if (cs%id_En_ang_mode(fr,m) > 0) then
1129 call post_data(cs%id_En_ang_mode(fr,m), cs%En(:,:,:,fr,m) , cs%diag)
1130 endif ; enddo ; enddo
1131
1132 ! Output 2-D energy loss (summed over angles, freq, modes)
1133 tot_leak_loss(:,:) = 0.0
1134 tot_quad_loss(:,:) = 0.0
1135 tot_itidal_loss(:,:) = 0.0
1136 tot_froude_loss(:,:) = 0.0
1137 tot_residual_loss(:,:) = 0.0
1138 tot_allprocesses_loss(:,:) = 0.0
1139 do m=1,cs%nMode ; do fr=1,cs%Nfreq ; do a=1,cs%nAngle ; do j=js,je ; do i=is,ie
1140 tot_leak_loss(i,j) = tot_leak_loss(i,j) + cs%TKE_leak_loss(i,j,a,fr,m)
1141 tot_quad_loss(i,j) = tot_quad_loss(i,j) + cs%TKE_quad_loss(i,j,a,fr,m)
1142 tot_itidal_loss(i,j) = tot_itidal_loss(i,j) + cs%TKE_itidal_loss(i,j,a,fr,m)
1143 tot_froude_loss(i,j) = tot_froude_loss(i,j) + cs%TKE_Froude_loss(i,j,a,fr,m)
1144 tot_residual_loss(i,j) = tot_residual_loss(i,j) + cs%TKE_residual_loss(i,j,a,fr,m)
1145 enddo ; enddo ; enddo ; enddo ; enddo
1146 do j=js,je ; do i=is,ie
1147 tot_allprocesses_loss(i,j) = ((((tot_leak_loss(i,j) + tot_quad_loss(i,j)) + &
1148 tot_itidal_loss(i,j)) + tot_froude_loss(i,j)) + &
1149 tot_residual_loss(i,j))
1150 enddo ; enddo
1151 cs%tot_leak_loss = tot_leak_loss
1152 cs%tot_quad_loss = tot_quad_loss
1153 cs%tot_itidal_loss = tot_itidal_loss
1154 cs%tot_Froude_loss = tot_froude_loss
1155 cs%tot_residual_loss = tot_residual_loss
1156 cs%tot_allprocesses_loss = tot_allprocesses_loss
1157 if (cs%id_tot_leak_loss > 0) then
1158 call post_data(cs%id_tot_leak_loss, tot_leak_loss, cs%diag)
1159 endif
1160 if (cs%id_tot_quad_loss > 0) then
1161 call post_data(cs%id_tot_quad_loss, tot_quad_loss, cs%diag)
1162 endif
1163 if (cs%id_tot_itidal_loss > 0) then
1164 call post_data(cs%id_tot_itidal_loss, tot_itidal_loss, cs%diag)
1165 endif
1166 if (cs%id_tot_Froude_loss > 0) then
1167 call post_data(cs%id_tot_Froude_loss, tot_froude_loss, cs%diag)
1168 endif
1169 if (cs%id_tot_residual_loss > 0) then
1170 call post_data(cs%id_tot_residual_loss, tot_residual_loss, cs%diag)
1171 endif
1172 if (cs%id_tot_allprocesses_loss > 0) then
1173 call post_data(cs%id_tot_allprocesses_loss, tot_allprocesses_loss, cs%diag)
1174 endif
1175
1176 ! Output 2-D energy loss (summed over angles) for each frequency and mode
1177 do m=1,cs%nMode ; do fr=1,cs%Nfreq
1178 if (cs%id_itidal_loss_mode(fr,m) > 0 .or. cs%id_allprocesses_loss_mode(fr,m) > 0) then
1179 itidal_loss_mode(:,:) = 0.0 ! wave-drag processes (could do others as well)
1180 leak_loss_mode(:,:) = 0.0
1181 quad_loss_mode(:,:) = 0.0
1182 froude_loss_mode(:,:) = 0.0
1183 residual_loss_mode(:,:) = 0.0
1184 allprocesses_loss_mode(:,:) = 0.0 ! all processes summed together
1185 do a=1,cs%nAngle ; do j=js,je ; do i=is,ie
1186 itidal_loss_mode(i,j) = itidal_loss_mode(i,j) + cs%TKE_itidal_loss(i,j,a,fr,m)
1187 leak_loss_mode(i,j) = leak_loss_mode(i,j) + cs%TKE_leak_loss(i,j,a,fr,m)
1188 quad_loss_mode(i,j) = quad_loss_mode(i,j) + cs%TKE_quad_loss(i,j,a,fr,m)
1189 froude_loss_mode(i,j) = froude_loss_mode(i,j) + cs%TKE_Froude_loss(i,j,a,fr,m)
1190 residual_loss_mode(i,j) = residual_loss_mode(i,j) + cs%TKE_residual_loss(i,j,a,fr,m)
1191 allprocesses_loss_mode(i,j) = allprocesses_loss_mode(i,j) + &
1192 ((((cs%TKE_leak_loss(i,j,a,fr,m) + cs%TKE_quad_loss(i,j,a,fr,m)) + &
1193 cs%TKE_itidal_loss(i,j,a,fr,m)) + cs%TKE_Froude_loss(i,j,a,fr,m)) + &
1194 cs%TKE_residual_loss(i,j,a,fr,m))
1195 enddo ; enddo ; enddo
1196 call post_data(cs%id_itidal_loss_mode(fr,m), itidal_loss_mode, cs%diag)
1197 call post_data(cs%id_leak_loss_mode(fr,m), leak_loss_mode, cs%diag)
1198 call post_data(cs%id_quad_loss_mode(fr,m), quad_loss_mode, cs%diag)
1199 call post_data(cs%id_Froude_loss_mode(fr,m), froude_loss_mode, cs%diag)
1200 call post_data(cs%id_residual_loss_mode(fr,m), residual_loss_mode, cs%diag)
1201 call post_data(cs%id_allprocesses_loss_mode(fr,m), allprocesses_loss_mode, cs%diag)
1202 endif ; enddo ; enddo
1203
1204 ! Output 3-D (i,j,a) energy loss for each frequency and mode
1205 do m=1,cs%nMode ; do fr=1,cs%Nfreq ; if (cs%id_itidal_loss_ang_mode(fr,m) > 0) then
1206 call post_data(cs%id_itidal_loss_ang_mode(fr,m), cs%TKE_itidal_loss(:,:,:,fr,m) , cs%diag)
1207 endif ; enddo ; enddo
1208
1209 ! Output 2-D period-averaged horizontal near-bottom mode velocity for each frequency and mode
1210 do m=1,cs%nMode ; do fr=1,cs%Nfreq ; if (cs%id_Ub_mode(fr,m) > 0) then
1211 call post_data(cs%id_Ub_mode(fr,m), ub(:,:,fr,m), cs%diag)
1212 endif ; enddo ; enddo
1213
1214 do m=1,cs%nMode ; if (cs%id_Ustruct_mode(m) > 0) then
1215 call post_data(cs%id_Ustruct_mode(m), cs%u_struct(:,:,:,m), cs%diag)
1216 endif ; enddo
1217
1218 do m=1,cs%nMode ; if (cs%id_Wstruct_mode(m) > 0) then
1219 call post_data(cs%id_Wstruct_mode(m), cs%w_struct(:,:,:,m), cs%diag)
1220 endif ; enddo
1221
1222 do m=1,cs%nMode ; if (cs%id_int_w2_mode(m) > 0) then
1223 call post_data(cs%id_int_w2_mode(m), cs%int_w2(:,:,m), cs%diag)
1224 endif ; enddo
1225
1226 do m=1,cs%nMode ; if (cs%id_int_U2_mode(m) > 0) then
1227 call post_data(cs%id_int_U2_mode(m), cs%int_U2(:,:,m), cs%diag)
1228 endif ; enddo
1229
1230 do m=1,cs%nMode ; if (cs%id_int_N2w2_mode(m) > 0) then
1231 call post_data(cs%id_int_N2w2_mode(m), cs%int_N2w2(:,:,m), cs%diag)
1232 endif ; enddo
1233
1234 ! Output 2-D horizontal phase velocity for each frequency and mode
1235 do m=1,cs%nMode ; do fr=1,cs%Nfreq ; if (cs%id_cp_mode(fr,m) > 0) then
1236 call post_data(cs%id_cp_mode(fr,m), cs%cp(:,:,fr,m), cs%diag)
1237 endif ; enddo ; enddo
1238
1239 endif
1240
1241 call disable_averaging(cs%diag)
1242
1243end subroutine propagate_int_tide
1244
1245!> Checks for energy conservation on computational domain
1246subroutine sum_en(G, GV, US, CS, En, label)
1247 type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure.
1248 type(verticalgrid_type),intent(in) :: GV !< The ocean's vertical grid structure.
1249 type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
1250 type(int_tide_cs), intent(inout) :: CS !< Internal tide control structure
1251 real, dimension(G%isd:G%ied,G%jsd:G%jed,CS%NAngle), &
1252 intent(in) :: En !< The energy density of the internal tides [H Z2 T-2 ~> m3 s-2 or J m-2].
1253 character(len=*), intent(in) :: label !< A label to use in error messages
1254 ! Local variables
1255 real :: En_sum ! The total energy in MKS units for potential output [m5 s-2 or J]
1256 integer :: a
1257 ! real :: En_sum_diff ! Change in energy from the expected value [m5 s-2 or J]
1258 ! real :: En_sum_pdiff ! Percentage change in energy from the expected value [nondim]
1259 ! character(len=160) :: mesg ! The text of an error message
1260 ! real :: days ! The time in days for use in output messages [days]
1261
1262 en_sum = 0.0
1263 do a=1,cs%nAngle
1264 en_sum = en_sum + global_area_integral(en(:,:,a), g, unscale=gv%H_to_mks*(us%Z_to_m**2)*(us%s_to_T)**2)
1265 enddo
1266 cs%En_sum = en_sum
1267 !En_sum_diff = En_sum - CS%En_sum
1268 !if (CS%En_sum /= 0.0) then
1269 ! En_sum_pdiff= (En_sum_diff/CS%En_sum)*100.0
1270 !else
1271 ! En_sum_pdiff= 0.0
1272 !endif
1273 !! Print to screen
1274 !if (is_root_pe()) then
1275 ! days = time_type_to_real(CS%Time) / 86400.0
1276 ! write(mesg,*) trim(label)//': days =', days, ', En_sum=', En_sum, &
1277 ! ', En_sum_diff=', En_sum_diff, ', Percent change=', En_sum_pdiff, '%'
1278 ! call MOM_mesg(mesg)
1279 !if (is_root_pe() .and. (abs(En_sum_pdiff) > 1.0)) &
1280 ! call MOM_error(FATAL, "Run stopped due to excessive internal tide energy change.")
1281 !endif
1282
1283end subroutine sum_en
1284
1285!> Calculates the energy lost from the propagating internal tide due to
1286!! scattering over small-scale roughness along the lines of Jayne & St. Laurent (2001).
1287subroutine itidal_lowmode_loss(G, GV, US, CS, Nb, Rho_bot, Ub, En, TKE_loss_fixed, TKE_loss, dt, halo_size)
1288 type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure.
1289 type(verticalgrid_type), intent(in) :: GV !< The ocean's vertical grid structure.
1290 type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
1291 type(int_tide_cs), intent(in) :: CS !< Internal tide control structure
1292 real, dimension(G%isd:G%ied,G%jsd:G%jed), &
1293 intent(in) :: Nb !< Near-bottom stratification [T-1 ~> s-1].
1294 real, dimension(G%isd:G%ied,G%jsd:G%jed), &
1295 intent(in) :: Rho_bot !< Near-bottom density [R ~> kg m-3].
1296 real, dimension(G%isd:G%ied,G%jsd:G%jed,CS%nFreq,CS%nMode), &
1297 intent(inout) :: Ub !< RMS (over one period) near-bottom horizontal
1298 !! mode velocity [L T-1 ~> m s-1].
1299 real, dimension(G%isd:G%ied,G%jsd:G%jed), &
1300 intent(in) :: TKE_loss_fixed !< Fixed part of energy loss [R Z4 H-1 L-2 ~> kg m-2 or m]
1301 !! (rho*kappa*h^2) or (kappa*h^2).
1302 real, dimension(G%isd:G%ied,G%jsd:G%jed,CS%NAngle,CS%nFreq,CS%nMode), &
1303 intent(inout) :: En !< Energy density of the internal waves [H Z2 T-2 ~> m3 s-2 or J m-2].
1304 real, dimension(G%isd:G%ied,G%jsd:G%jed,CS%NAngle,CS%nFreq,CS%nMode), &
1305 intent(out) :: TKE_loss !< Energy loss rate [H Z2 T-3 ~> m3 s-3 or W m-2]
1306 !! (q*rho*kappa*h^2*N*U^2).
1307 real, intent(in) :: dt !< Time increment [T ~> s].
1308 integer, optional, intent(in) :: halo_size !< The halo size over which to do the calculations
1309 ! Local variables
1310 integer :: j, i, m, fr, a, is, ie, js, je, halo
1311 real :: En_tot ! energy for a given mode, frequency
1312 ! and point summed over angles [H Z2 T-2 ~> m3 s-2 or J m-2]
1313 real :: TKE_loss_tot ! dissipation for a given mode, frequency
1314 ! and point summed over angles [H Z2 T-3 ~> m3 s-3 or W m-2]
1315 real :: frac_per_sector ! fraction of energy in each wedge [nondim]
1316 real :: q_itides ! fraction of energy actually lost to mixing (remainder, 1-q, is
1317 ! assumed to stay in propagating mode for now - BDM) [nondim]
1318 real :: loss_rate ! approximate loss rate for implicit calc [T-1 ~> s-1]
1319 real :: En_negl ! negligibly small number to prevent division by zero [H Z2 T-2 ~> m3 s-2 or J m-2]
1320 real :: En_a, En_b ! energy before and after timestep [H Z2 T-2 ~> m3 s-2 or J m-2]
1321 real :: I_dt ! The inverse of the timestep [T-1 ~> s-1]
1322 real :: J_m2_to_HZ2_T2 ! unit conversion factor for Energy from mks to internal
1323 ! units [H Z2 s2 T-2 kg-1 ~> m3 kg-1 or 1]
1324
1325 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
1326
1327 j_m2_to_hz2_t2 = gv%m_to_H*(us%m_to_Z**2)*(us%T_to_s**2)
1328
1329 i_dt = 1.0 / dt
1330 q_itides = cs%q_itides
1331 en_negl = 1e-30*j_m2_to_hz2_t2
1332
1333 if (present(halo_size)) then
1334 halo = halo_size
1335 is = g%isc - halo ; ie = g%iec + halo ; js = g%jsc - halo ; je = g%jec + halo
1336 endif
1337
1338 if (cs%debug) then
1339 call hchksum(tke_loss_fixed, "TKE loss fixed", g%HI, haloshift=0, &
1340 unscale=us%RZ_to_kg_m2*(us%Z_to_m**3)*gv%m_to_H*(us%m_to_L**2))
1341 call hchksum(nb(:,:), "Nbottom", g%HI, haloshift=0, unscale=us%s_to_T)
1342 call hchksum(ub(:,:,1,1), "Ubottom", g%HI, haloshift=0, unscale=us%L_to_m*us%s_to_T)
1343 endif
1344
1345 do j=js,je ; do i=is,ie ; do m=1,cs%nMode ; do fr=1,cs%nFreq
1346
1347 ! Sum energy across angles
1348 en_tot = 0.0
1349 do a=1,cs%nAngle
1350 en_tot = en_tot + en(i,j,a,fr,m)
1351 enddo
1352
1353 ! Calculate TKE loss rate; units of [H Z2 T-3 ~> m3 s-3 or W m-2] here.
1354 if (gv%Boussinesq .or. gv%semi_Boussinesq) then
1355 tke_loss_tot = q_itides * gv%RZ_to_H*gv%Z_to_H*tke_loss_fixed(i,j)*nb(i,j)*ub(i,j,fr,m)**2
1356 else
1357 tke_loss_tot = q_itides * (gv%RZ_to_H*gv%RZ_to_H*rho_bot(i,j))*tke_loss_fixed(i,j)*nb(i,j)*ub(i,j,fr,m)**2
1358 endif
1359
1360 ! Update energy remaining (this is a pseudo implicit calc)
1361 ! (E(t+1)-E(t))/dt = -TKE_loss(E(t+1)/E(t)), which goes to zero as E(t+1) goes to zero
1362 if (en_tot > 0.0) then
1363 do a=1,cs%nAngle
1364 frac_per_sector = en(i,j,a,fr,m)/en_tot
1365 tke_loss(i,j,a,fr,m) = frac_per_sector*tke_loss_tot ! [H Z2 T-3 ~> m3 s-3 or W m-2]
1366 loss_rate = tke_loss(i,j,a,fr,m) / (en(i,j,a,fr,m) + en_negl) ! [T-1 ~> s-1]
1367 en_b = en(i,j,a,fr,m)
1368 en_a = en(i,j,a,fr,m) / (1.0 + (dt*loss_rate))
1369 tke_loss(i,j,a,fr,m) = (en_b - en_a) * i_dt ! overwrite with exact value
1370 en(i,j,a,fr,m) = en_a
1371 enddo
1372 else
1373 ! no loss if no energy
1374 do a=1,cs%nAngle
1375 tke_loss(i,j,a,fr,m) = 0.0
1376 enddo
1377 endif
1378
1379 enddo ; enddo ; enddo ; enddo
1380
1381end subroutine itidal_lowmode_loss
1382
1383!> This subroutine extracts the energy lost from the propagating internal which has
1384!> been summed across all angles, frequencies, and modes for a given mechanism and location.
1385!!
1386!> It can be called from another module to get values from this module's (private) CS.
1387subroutine get_lowmode_loss(i,j,G,CS,mechanism,TKE_loss_sum)
1388 integer, intent(in) :: i !< The i-index of the value to be reported.
1389 integer, intent(in) :: j !< The j-index of the value to be reported.
1390 type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure
1391 type(int_tide_cs), intent(in) :: cs !< Internal tide control structure
1392 character(len=*), intent(in) :: mechanism !< The named mechanism of loss to return
1393 real, intent(out) :: tke_loss_sum !< Total energy loss rate due to specified
1394 !! mechanism [H Z2 T-3 ~> m3 s-3 or W m-2].
1395
1396 if (mechanism == 'LeakDrag') tke_loss_sum = cs%tot_leak_loss(i,j)
1397 if (mechanism == 'QuadDrag') tke_loss_sum = cs%tot_quad_loss(i,j)
1398 if (mechanism == 'WaveDrag') tke_loss_sum = cs%tot_itidal_loss(i,j)
1399 if (mechanism == 'Froude') tke_loss_sum = cs%tot_Froude_loss(i,j)
1400 if (mechanism == 'SlopeDrag') tke_loss_sum = cs%tot_residual_loss(i,j)
1401
1402end subroutine get_lowmode_loss
1403
1404
1405!> Returns the values of diffusivity corresponding to various mechanisms
1406subroutine get_lowmode_diffusivity(G, GV, h, tv, US, h_bot, k_bot, j, N2_lay, N2_int, TKE_to_Kd, Kd_max, CS, &
1407 Kd_leak, Kd_quad, Kd_itidal, Kd_Froude, Kd_slope, &
1408 Kd_lay, Kd_int, profile_leak, profile_quad, profile_itidal, &
1409 profile_Froude, profile_slope)
1410
1411 type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure
1412 type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure
1413 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
1414 intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2]
1415 type(thermo_var_ptrs), intent(in) :: tv !< Structure containing pointers to any available
1416 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
1417 real, dimension(SZI_(G)), intent(in) :: h_bot !< Bottom boundary layer thickness [H ~> m or kg m-2]
1418 integer, dimension(SZI_(G)), intent(in) :: k_bot !< Bottom boundary layer top layer index
1419 integer, intent(in) :: j !< The j-index to work on
1420 real, dimension(SZI_(G),SZK_(GV)), intent(in) :: n2_lay !< The squared buoyancy frequency of the
1421 !! layers [T-2 ~> s-2].
1422 real, dimension(SZI_(G),SZK_(GV)+1), intent(in) :: n2_int !< The squared buoyancy frequency of the
1423 !! interfaces [T-2 ~> s-2].
1424 real, dimension(SZI_(G),SZK_(GV)), intent(in) :: tke_to_kd !< The conversion rate between the TKE
1425 !! dissipated within a layer and the
1426 !! diapycnal diffusivity within that layer,
1427 !! usually (~Rho_0 / (G_Earth * dRho_lay))
1428 !! [T2 Z-1 ~> s2 m-1]
1429 real, intent(in) :: kd_max !< The maximum increment for diapycnal
1430 !! diffusivity due to TKE-based processes
1431 !! [H Z T-1 ~> m2 s-1 or kg m-1 s-1].
1432 !! Set this to a negative value to have no limit.
1433 !! [H Z T-1 ~> m2 s-1 or kg m-1 s-1].
1434 type(int_tide_cs), intent(in) :: cs !< The control structure for this module
1435
1436 real, dimension(SZI_(G),SZK_(GV)+1), intent(out) :: kd_leak !< Diffusivity due to background drag
1437 !! [H Z T-1 ~> m2 s-1 or kg m-1 s-1].
1438 real, dimension(SZI_(G),SZK_(GV)+1), intent(out) :: kd_quad !< Diffusivity due to bottom drag
1439 !! [H Z T-1 ~> m2 s-1 or kg m-1 s-1].
1440 real, dimension(SZI_(G),SZK_(GV)+1), intent(out) :: kd_itidal !< Diffusivity due to wave drag
1441 !! [H Z T-1 ~> m2 s-1 or kg m-1 s-1].
1442 real, dimension(SZI_(G),SZK_(GV)+1), intent(out) :: kd_froude !< Diffusivity due to high Froude breaking
1443 !! [H Z T-1 ~> m2 s-1 or kg m-1 s-1].
1444 real, dimension(SZI_(G),SZK_(GV)+1), intent(out) :: kd_slope !< Diffusivity due to critical slopes
1445 !! [H Z T-1 ~> m2 s-1 or kg m-1 s-1].
1446 real, dimension(SZI_(G),SZK_(GV)), intent(inout) :: kd_lay !< The diapycnal diffusivity in layers
1447 !! [H Z T-1 ~> m2 s-1 or kg m-1 s-1].
1448 real, dimension(SZI_(G),SZK_(GV)+1), intent(inout) :: kd_int !< The diapycnal diffusivity at interfaces
1449 !! [H Z T-1 ~> m2 s-1 or kg m-1 s-1].
1450 real, dimension(SZI_(G), SZK_(GV)), intent(out) :: profile_leak !< Normalized profile for background drag
1451 !! [H-1 ~> m-1 or m2 kg-1]
1452 real, dimension(SZI_(G), SZK_(GV)), intent(out) :: profile_quad !< Normalized profile for bottom drag
1453 !! [H-1 ~> m-1 or m2 kg-1]
1454 real, dimension(SZI_(G), SZK_(GV)), intent(out) :: profile_itidal !< Normalized profile for wave drag
1455 !! [H-1 ~> m-1 or m2 kg-1]
1456 real, dimension(SZI_(G), SZK_(GV)), intent(out) :: profile_froude !< Normalized profile for Froude drag
1457 !! [H-1 ~> m-1 or m2 kg-1]
1458 real, dimension(SZI_(G), SZK_(GV)), intent(out) :: profile_slope !< Normalized profile for critical slopes
1459 !! [H-1 ~> m-1 or m2 kg-1]
1460
1461 ! local variables
1462 real :: tke_loss ! temp variable to pass value of internal tides TKE loss [H Z2 T-3 ~> m3 s-3 or W m-2]
1463 real :: renorm_n ! renormalization for N profile [H T-1 ~> m s-1 or kg m-2 s-1]
1464 real :: renorm_n2 ! renormalization for N2 profile [H T-2 ~> m s-2 or kg m-2 s-2]
1465 real :: tmp_stlau ! tmp var for renormalization for StLaurent profile [nondim]
1466 real :: tmp_stlau_slope ! tmp var for renormalization for StLaurent profile [nondim]
1467 real :: renorm_stlau ! renormalization for StLaurent profile [nondim]
1468 real :: renorm_stlau_slope! renormalization for StLaurent profile [nondim]
1469 real :: htot ! total depth of water column [H ~> m or kg m-2]
1470 real :: htmp ! local value of thickness in layers [H ~> m or kg m-2]
1471 real :: h_d ! expomential decay length scale [H ~> m or kg m-2]
1472 real :: h_s ! expomential decay length scale on the slope [H ~> m or kg m-2]
1473 real :: i_h_d ! inverse of expomential decay length scale [H-1 ~> m-1 or m2 kg-1]
1474 real :: i_h_s ! inverse of expomential decay length scale on the slope [H-1 ~> m-1 or m2 kg-1]
1475 real :: tke_to_kd_lim ! limited version of TKE_to_Kd [T2 Z-1 ~> s2 m-1]
1476
1477 ! vertical profiles have units Z-1 for conversion to Kd to be dim correct (see eq 2 of St Laurent GRL 2002)
1478 real, dimension(SZK_(GV)) :: profile_n ! vertical profile varying with N [H-1 ~> m-1 or m2 kg-1]
1479 real, dimension(SZK_(GV)) :: profile_n2 ! vertical profile varying with N2 [H-1 ~> m-1 or m2 kg-1]
1480 real, dimension(SZK_(GV)) :: profile_stlaurent ! vertical profile according to St Laurent 2002
1481 ! [H-1 ~> m-1 or m2 kg-1]
1482 real, dimension(SZK_(GV)) :: profile_stlaurent_slope ! vertical profile according to St Laurent 2002
1483 ! [H-1 ~> m-1 or m2 kg-1]
1484 real, dimension(SZK_(GV)) :: profile_bbl ! vertical profile Heavyside BBL [H-1 ~> m-1 or m2 kg-1]
1485 real, dimension(SZK_(GV)) :: kd_leak_lay ! Diffusivity due to background drag [H Z T-1 ~> m2 s-1 or kg m-1 s-1]
1486 real, dimension(SZK_(GV)) :: kd_quad_lay ! Diffusivity due to bottom drag [H Z T-1 ~> m2 s-1 or kg m-1 s-1]
1487 real, dimension(SZK_(GV)) :: kd_itidal_lay ! Diffusivity due to wave drag [H Z T-1 ~> m2 s-1 or kg m-1 s-1]
1488 real, dimension(SZK_(GV)) :: kd_froude_lay ! Diffusivity due to high Froude breaking [H Z T-1 ~> m2 s-1 or kg m-1 s-1]
1489 real, dimension(SZK_(GV)) :: kd_slope_lay ! Diffusivity due to critical slopes [H Z T-1 ~> m2 s-1 or kg m-1 s-1]
1490
1491 real :: hmin ! A minimum allowable thickness [H ~> m or kg m-2]
1492 real :: h_rmn ! Remaining thickness in k-loop [H ~> m or kg m-2]
1493 real :: frac ! A fraction of thicknesses [nondim]
1494 real :: i_h_bot ! inverse of Bottom boundary layer thickness [H-1 ~> m-1 or m2 kg-1]
1495
1496 real :: verif_n, & ! profile verification [nondim]
1497 verif_n2, & ! profile verification [nondim]
1498 verif_bbl, & ! profile verification [nondim]
1499 verif_stl1,& ! profile verification [nondim]
1500 verif_stl2,& ! profile verification [nondim]
1501 threshold_renorm_n2,& ! Maximum allowable error on N2 profile [H T-2 ~> m s-2 or kg m-2 s-2]
1502 threshold_renorm_n, & ! Maximum allowable error on N profile [H T-1 ~> m s-1 or kg m-2 s-1]
1503 threshold_verif ! Maximum allowable error on verification [nondim]
1504
1505 logical :: non_bous ! fully Non-Boussinesq
1506 integer :: i, k, is, ie, nz
1507
1508 is=g%isc ; ie=g%iec ; nz=gv%ke
1509
1510 non_bous = .not.(gv%Boussinesq .or. gv%semi_Boussinesq)
1511
1512 h_d = cs%Int_tide_decay_scale
1513 h_s = cs%Int_tide_decay_scale_slope
1514 i_h_d = 1 / h_d
1515 i_h_s = 1 / h_s
1516
1517 hmin = 1.0e-6*gv%m_to_H
1518 threshold_renorm_n2 = 1.0e-13 * gv%m_to_H * us%T_to_s**2
1519 threshold_renorm_n = 1.0e-13 * gv%m_to_H * us%T_to_s
1520 threshold_verif = 1.0e-13
1521
1522 ! init output arrays
1523 profile_leak(:,:) = 0.0
1524 profile_quad(:,:) = 0.0
1525 profile_slope(:,:) = 0.0
1526 profile_itidal(:,:) = 0.0
1527 profile_froude(:,:) = 0.0
1528
1529 kd_leak_lay(:) = 0.0
1530 kd_quad_lay(:) = 0.0
1531 kd_itidal_lay(:) = 0.0
1532 kd_froude_lay(:) = 0.0
1533 kd_slope_lay(:) = 0.0
1534
1535 kd_leak(:,:) = 0.0
1536 kd_quad(:,:) = 0.0
1537 kd_itidal(:,:) = 0.0
1538 kd_froude(:,:) = 0.0
1539 kd_slope(:,:) = 0.0
1540
1541 do i=is,ie
1542
1543 ! create vertical profiles for diffusivities in layers
1544 renorm_n = 0.0
1545 renorm_n2 = 0.0
1546 renorm_stlau = 0.0
1547 renorm_stlau_slope = 0.0
1548 tmp_stlau = 0.0
1549 tmp_stlau_slope = 0.0
1550 htot = 0.0
1551 htmp = 0.0
1552 i_h_bot = 1.0 / h_bot(i)
1553
1554 do k=1,nz
1555 ! N-profile
1556 if (n2_lay(i,k) < 0.) call mom_error(warning, "negative buoyancy freq")
1557 renorm_n = renorm_n + (sqrt(max(n2_lay(i,k), 0.)) * h(i,j,k))
1558 ! N2-profile
1559 renorm_n2 = renorm_n2 + (max(n2_lay(i,k), 0.) * h(i,j,k))
1560 ! total depth
1561 htot = htot + h(i,j,k)
1562 enddo
1563
1564 profile_n2(:) = 0.0
1565 profile_n(:) = 0.0
1566 profile_bbl(:) = 0.0
1567 profile_stlaurent(:) = 0.0
1568 profile_stlaurent_slope(:) = 0.0
1569
1570 ! BBL-profile
1571 h_rmn = h_bot(i)
1572 do k=nz,1,-1
1573 if (g%mask2dT(i,j) > 0.0) then
1574 profile_bbl(k) = 0.0
1575 if (h(i,j,k) <= h_rmn) then
1576 profile_bbl(k) = 1.0 * i_h_bot
1577 h_rmn = h_rmn - h(i,j,k)
1578 else
1579 if (h_rmn > 0.0) then
1580 frac = h_rmn / h(i,j,k)
1581 profile_bbl(k) = frac * i_h_bot
1582 h_rmn = h_rmn - frac*h(i,j,k)
1583 endif
1584 endif
1585 endif
1586 enddo
1587
1588 do k=1,nz
1589 if (g%mask2dT(i,j) > 0.0) then
1590 ! N - profile
1591 if (renorm_n > threshold_renorm_n) then
1592 profile_n(k) = sqrt(max(n2_lay(i,k), 0.)) / renorm_n
1593 else
1594 profile_n(k) = 1 / htot
1595 endif
1596
1597 ! N2 - profile
1598 if (renorm_n2 > threshold_renorm_n2) then
1599 profile_n2(k) = max(n2_lay(i,k), 0.) / renorm_n2
1600 else
1601 profile_n2(k) = 1 / htot
1602 endif
1603
1604 ! slope intensified (St Laurent GRL 2002) - profile
1605 ! in paper, z is defined positive upwards, range 0 to -H
1606 ! here depth positive downwards
1607 ! profiles are almost normalized but differ from a few percent
1608 ! so we add a second renormalization factor
1609
1610 ! add first half of layer: get to the layer center
1611 htmp = htmp + 0.5*h(i,j,k)
1612
1613 profile_stlaurent(k) = exp(-i_h_d*(htot-htmp)) / &
1614 (h_d*(1 - exp(-i_h_d*htot)))
1615
1616 profile_stlaurent_slope(k) = exp(-i_h_s*(htot-htmp)) / &
1617 (h_s*(1 - exp(-i_h_s*htot)))
1618
1619 tmp_stlau = tmp_stlau + (profile_stlaurent(k) * h(i,j,k))
1620 tmp_stlau_slope = tmp_stlau_slope + (profile_stlaurent_slope(k) * h(i,j,k))
1621
1622 ! add second half of layer: get to the next interface
1623 htmp = htmp + 0.5*h(i,j,k)
1624 endif
1625 enddo
1626
1627 if (g%mask2dT(i,j) > 0.0) then
1628 ! allow for difference less than verification threshold
1629 renorm_stlau = 1.0
1630 renorm_stlau_slope = 1.0
1631 if (abs(tmp_stlau -1.0) > threshold_verif) renorm_stlau = 1.0 / tmp_stlau
1632 if (abs(tmp_stlau_slope -1.0) > threshold_verif) renorm_stlau_slope = 1.0 / tmp_stlau_slope
1633
1634 do k=1,nz
1635 profile_stlaurent(k) = profile_stlaurent(k) * renorm_stlau
1636 profile_stlaurent_slope(k) = profile_stlaurent_slope(k) * renorm_stlau_slope
1637 enddo
1638 endif
1639
1640 ! verif integrals
1641 if (cs%debug) then
1642 if (g%mask2dT(i,j) > 0.0) then
1643 verif_n = 0.0
1644 verif_n2 = 0.0
1645 verif_bbl = 0.0
1646 verif_stl1 = 0.0
1647 verif_stl2 = 0.0
1648 do k=1,nz
1649 verif_n = verif_n + (profile_n(k) * h(i,j,k))
1650 verif_n2 = verif_n2 + (profile_n2(k) * h(i,j,k))
1651 verif_bbl = verif_bbl + (profile_bbl(k) * h(i,j,k))
1652 verif_stl1 = verif_stl1 + (profile_stlaurent(k) * h(i,j,k))
1653 verif_stl2 = verif_stl2 + (profile_stlaurent_slope(k) * h(i,j,k))
1654 enddo
1655
1656 if (abs(verif_n -1.0) > threshold_verif) then
1657 write(stdout,'(I0,", ",I0,F18.10)') i, j, verif_n
1658 call mom_error(fatal, "mismatch integral for N profile")
1659 endif
1660 if (abs(verif_n2 -1.0) > threshold_verif) then
1661 write(stdout,'(I0,", ",I0,F18.10)') i, j, verif_n2
1662 call mom_error(fatal, "mismatch integral for N2 profile")
1663 endif
1664 if (abs(verif_bbl -1.0) > threshold_verif) then
1665 write(stdout,'(I0,", ",I0,F18.10)') i, j, verif_bbl
1666 call mom_error(fatal, "mismatch integral for bbl profile")
1667 endif
1668 if (abs(verif_stl1 -1.0) > threshold_verif) then
1669 write(stdout,'(I0,", ",I0,F18.10)') i, j, verif_stl1
1670 call mom_error(fatal, "mismatch integral for stl1 profile")
1671 endif
1672 if (abs(verif_stl2 -1.0) > threshold_verif) then
1673 write(stdout,'(I0,", ",I0,F18.10)') i, j, verif_stl2
1674 call mom_error(fatal, "mismatch integral for stl2 profile")
1675 endif
1676
1677 endif
1678 endif
1679
1680 ! note on units: TKE_to_Kd = 1 / ((g/rho0) * drho) Z-1 T2
1681 ! mult by dz gives -1/N2 in T2
1682
1683 ! get TKE loss value and compute diffusivities in layers
1684 if (cs%apply_background_drag) then
1685 call get_lowmode_loss(i, j, g, cs, "LeakDrag", tke_loss)
1686 ! insert logic to switch between profiles here
1687 ! if trim(CS%leak_profile) == "N2" then
1688 profile_leak(i,:) = profile_n2(:)
1689 ! elseif trim(CS%leak_profile) == "N" then
1690 ! profile_leak(:) = profile_N(:)
1691 ! something else
1692 ! endif
1693 kd_leak_lay(:) = 0.
1694 do k=1,nz
1695 ! layer diffusivity for processus
1696 if (h(i,j,k) >= cs%min_thick_layer_Kd) then
1697 tke_to_kd_lim = min(tke_to_kd(i,k), cs%max_TKE_to_Kd)
1698 kd_leak_lay(k) = cs%mixing_effic * tke_loss * tke_to_kd_lim * profile_leak(i,k) * h(i,j,k)
1699 else
1700 kd_leak_lay(k) = 0.
1701 endif
1702 ! add to total Kd in layer
1703 if (cs%update_Kd) kd_lay(i,k) = kd_lay(i,k) + min(kd_leak_lay(k), kd_max)
1704 enddo
1705 endif
1706
1707 if (cs%apply_Froude_drag) then
1708 call get_lowmode_loss(i, j, g, cs, "Froude", tke_loss)
1709 ! insert logic to switch between profiles here
1710 ! if trim(CS%Froude_profile) == "N" then
1711 profile_froude(i,:) = profile_n(:)
1712 ! elseif trim(CS%Froude_profile) == "N2" then
1713 ! profile_Froude(:) = profile_N2(:)
1714 ! something else
1715 ! endif
1716 do k=1,nz
1717 ! layer diffusivity for processus
1718 if (h(i,j,k) >= cs%min_thick_layer_Kd) then
1719 tke_to_kd_lim = min(tke_to_kd(i,k), cs%max_TKE_to_Kd)
1720 kd_froude_lay(k) = cs%mixing_effic * tke_loss * tke_to_kd_lim * profile_froude(i,k) * h(i,j,k)
1721 else
1722 kd_froude_lay(k) = 0.
1723 endif
1724 ! add to total Kd in layer
1725 if (cs%update_Kd) kd_lay(i,k) = kd_lay(i,k) + min(kd_froude_lay(k), kd_max)
1726 enddo
1727 endif
1728
1729 if (cs%apply_wave_drag) then
1730 call get_lowmode_loss(i, j, g, cs, "WaveDrag", tke_loss)
1731 ! insert logic to switch between profiles here
1732 ! if trim(CS%wave_profile) == "StLaurent" then
1733 profile_itidal(i,:) = profile_stlaurent(:)
1734 ! elseif trim(CS%Froude_profile) == "N2" then
1735 ! profile_itidal(:) = profile_N2(:)
1736 ! something else
1737 ! endif
1738 do k=1,nz
1739 ! layer diffusivity for processus
1740 if (h(i,j,k) >= cs%min_thick_layer_Kd) then
1741 tke_to_kd_lim = min(tke_to_kd(i,k), cs%max_TKE_to_Kd)
1742 kd_itidal_lay(k) = cs%mixing_effic * tke_loss * tke_to_kd_lim * profile_itidal(i,k) * h(i,j,k)
1743 else
1744 kd_itidal_lay(k) = 0.
1745 endif
1746 ! add to total Kd in layer
1747 if (cs%update_Kd) kd_lay(i,k) = kd_lay(i,k) + min(kd_itidal_lay(k), kd_max)
1748 enddo
1749 endif
1750
1751 if (cs%apply_residual_drag) then
1752 call get_lowmode_loss(i, j, g, cs, "SlopeDrag", tke_loss)
1753 ! insert logic to switch between profiles here
1754 ! if trim(CS%wave_profile) == "StLaurent" then
1755 profile_slope(i,:) = profile_stlaurent_slope(:)
1756 ! elseif trim(CS%Froude_profile) == "N2" then
1757 ! profile_itidal(:) = profile_N2(:)
1758 ! something else
1759 ! endif
1760 do k=1,nz
1761 ! layer diffusivity for processus
1762 if (h(i,j,k) >= cs%min_thick_layer_Kd) then
1763 tke_to_kd_lim = min(tke_to_kd(i,k), cs%max_TKE_to_Kd)
1764 kd_slope_lay(k) = cs%mixing_effic * tke_loss * tke_to_kd_lim * profile_slope(i,k) * h(i,j,k)
1765 else
1766 kd_slope_lay(k) = 0.
1767 endif
1768 ! add to total Kd in layer
1769 if (cs%update_Kd) kd_lay(i,k) = kd_lay(i,k) + min(kd_slope_lay(k), kd_max)
1770 enddo
1771 endif
1772
1773 if (cs%apply_bottom_drag) then
1774 call get_lowmode_loss(i, j, g, cs, "QuadDrag", tke_loss)
1775 ! insert logic to switch between profiles here
1776 ! if trim(CS%bottom_profile) == "BBL" then
1777 profile_quad(i,:) = profile_bbl(:)
1778 ! elseif trim(CS%bottom_profile) == "N2" then
1779 ! profile_quad(:) = profile_N2(:)
1780 ! something else
1781 ! endif
1782 do k=1,nz
1783 ! layer diffusivity for processus
1784 if (h(i,j,k) >= cs%min_thick_layer_Kd) then
1785 tke_to_kd_lim = min(tke_to_kd(i,k), cs%max_TKE_to_Kd)
1786 kd_quad_lay(k) = cs%mixing_effic * tke_loss * tke_to_kd_lim * profile_quad(i,k) * h(i,j,k)
1787 else
1788 kd_quad_lay(k) = 0.
1789 endif
1790 ! add to total Kd in layer
1791 if (cs%update_Kd) kd_lay(i,k) = kd_lay(i,k) + min(kd_quad_lay(k), kd_max)
1792 enddo
1793 endif
1794
1795 ! interpolate Kd_[] to interfaces and add to Kd_int
1796 if (cs%apply_background_drag) then
1797 do k=1,nz+1
1798 if (k>1) kd_leak(i,k) = 0.5*kd_leak_lay(k-1)
1799 if (k<nz+1) kd_leak(i,k) = kd_leak(i,k) + 0.5*kd_leak_lay(k)
1800 ! add to Kd_int
1801 if (cs%update_Kd) kd_int(i,k) = kd_int(i,k) + min(kd_leak(i,k), kd_max)
1802 enddo
1803 endif
1804
1805 if (cs%apply_wave_drag) then
1806 do k=1,nz+1
1807 if (k>1) kd_itidal(i,k) = 0.5*kd_itidal_lay(k-1)
1808 if (k<nz+1) kd_itidal(i,k) = kd_itidal(i,k) + 0.5*kd_itidal_lay(k)
1809 ! add to Kd_int
1810 if (cs%update_Kd) kd_int(i,k) = kd_int(i,k) + min(kd_itidal(i,k), kd_max)
1811 enddo
1812 endif
1813
1814 if (cs%apply_Froude_drag) then
1815 do k=1,nz+1
1816 if (k>1) kd_froude(i,k) = 0.5*kd_froude_lay(k-1)
1817 if (k<nz+1) kd_froude(i,k) = kd_froude(i,k) + 0.5*kd_froude_lay(k)
1818 ! add to Kd_int
1819 if (cs%update_Kd) kd_int(i,k) = kd_int(i,k) + min(kd_froude(i,k), kd_max)
1820 enddo
1821 endif
1822
1823 if (cs%apply_residual_drag) then
1824 do k=1,nz+1
1825 if (k>1) kd_slope(i,k) = 0.5*kd_slope_lay(k-1)
1826 if (k<nz+1) kd_slope(i,k) = kd_slope(i,k) + 0.5*kd_slope_lay(k)
1827 ! add to Kd_int
1828 if (cs%update_Kd) kd_int(i,k) = kd_int(i,k) + min(kd_slope(i,k), kd_max)
1829 enddo
1830 endif
1831
1832 if (cs%apply_bottom_drag) then
1833 do k=1,nz+1
1834 if (k>1) kd_quad(i,k) = 0.5*kd_quad_lay(k-1)
1835 if (k<nz+1) kd_quad(i,k) = kd_quad(i,k) + 0.5*kd_quad_lay(k)
1836 ! add to Kd_int
1837 if (cs%update_Kd) kd_int(i,k) = kd_int(i,k) + min(kd_quad(i,k), kd_max)
1838 enddo
1839 endif
1840 enddo ! i-loop
1841
1842end subroutine get_lowmode_diffusivity
1843
1844!> Implements refraction on the internal waves at a single frequency.
1845subroutine refract(En, cn, freq, dt, G, US, NAngle, use_PPMang)
1846 type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure.
1847 integer, intent(in) :: NAngle !< The number of wave orientations in the
1848 !! discretized wave energy spectrum.
1849 real, dimension(G%isd:G%ied,G%jsd:G%jed,NAngle), &
1850 intent(inout) :: En !< The internal gravity wave energy density as a
1851 !! function of space and angular resolution,
1852 !! [H Z2 T-2 ~> m3 s-2 or J m-2].
1853 real, dimension(G%isd:G%ied,G%jsd:G%jed), &
1854 intent(in) :: cn !< Baroclinic mode speed [L T-1 ~> m s-1].
1855 real, intent(in) :: freq !< Wave frequency [T-1 ~> s-1].
1856 real, intent(in) :: dt !< Time step [T ~> s].
1857 type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
1858 logical, intent(in) :: use_PPMang !< If true, use PPM for advection rather
1859 !! than upwind.
1860 ! Local variables
1861 integer, parameter :: stencil = 2
1862 real, dimension(SZI_(G),1-stencil:NAngle+stencil) :: &
1863 En2d ! The internal gravity wave energy density in zonal slices [H Z2 T-2 ~> m3 s-2 or J m-2]
1864 real, dimension(1-stencil:NAngle+stencil) :: &
1865 cos_angle, sin_angle ! The cosine and sine of each angle [nondim]
1866 real, dimension(SZI_(G)) :: &
1867 Dk_Dt_Kmag, Dl_Dt_Kmag ! Rates of angular refraction [T-1 ~> s-1]
1868 real, dimension(SZI_(G),0:nAngle) :: &
1869 Flux_E ! The flux of energy between successive angular wedges
1870 ! within a timestep [H Z2 T-2 ~> m3 s-2 or J m-2]
1871 real, dimension(SZI_(G),SZJ_(G),1-stencil:NAngle+stencil) :: &
1872 CFL_ang ! The CFL number of angular refraction [nondim]
1873 real, dimension(G%IsdB:G%IedB,G%jsd:G%jed) :: cn_u !< Internal wave group velocity at U-point [L T-1 ~> m s-1]
1874 real, dimension(G%isd:G%ied,G%JsdB:G%JedB) :: cn_v !< Internal wave group velocity at V-point [L T-1 ~> m s-1]
1875 real, dimension(G%isd:G%ied,G%jsd:G%jed) :: cnmask !< Local mask for group velocity [nondim]
1876 real :: f2 ! The squared Coriolis parameter [T-2 ~> s-2].
1877 real :: favg ! The average Coriolis parameter at a point [T-1 ~> s-1].
1878 real :: df_dy, df_dx ! The x- and y- gradients of the Coriolis parameter [T-1 L-1 ~> s-1 m-1].
1879 real :: dlnCn_dx ! The x-gradient of the wave speed divided by itself [L-1 ~> m-1].
1880 real :: dlnCn_dy ! The y-gradient of the wave speed divided by itself [L-1 ~> m-1].
1881 real :: Angle_size ! The size of each wedge of angles [rad]
1882 real :: dt_Angle_size ! The time step divided by the angle size [T rad-1 ~> s rad-1]
1883 real :: angle ! The central angle of each wedge [rad]
1884 real :: Ifreq ! The inverse of the wave frequency [T ~> s]
1885 real :: Kmag2 ! A squared horizontal wavenumber [L-2 ~> m-2]
1886 real :: I_Kmag ! The inverse of the magnitude of the horizontal wavenumber [L ~> m]
1887 real :: cn_subRO ! A tiny wave speed to prevent division by zero [L T-1 ~> m s-1]
1888 integer :: is, ie, js, je, asd, aed, na
1889 integer :: i, j, a
1890 real :: wgt1, wgt2 ! Weights in an average, both of which may be 0 [nondim]
1891
1892 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; na = size(en,3)
1893 asd = 1-stencil ; aed = nangle+stencil
1894
1895 cnmask(:,:) = merge(0., 1., cn(:,:) == 0.)
1896
1897 do j=js,je ; do i=is-1,ie
1898 ! wgt = 0 if local cn == 0, wgt = 0.5 if both contiguous values != 0
1899 ! and wgt = 1 if neighbour cn == 0
1900 wgt1 = cnmask(i,j) - (0.5 * cnmask(i,j) * cnmask(i+1,j))
1901 wgt2 = cnmask(i+1,j) - (0.5 * cnmask(i,j) * cnmask(i+1,j))
1902 cn_u(i,j) = (wgt1*cn(i,j)) + (wgt2*cn(i+1,j))
1903 enddo ; enddo
1904
1905 do j=js-1,je ; do i=is,ie
1906 wgt1 = cnmask(i,j) - (0.5 * cnmask(i,j) * cnmask(i,j+1))
1907 wgt2 = cnmask(i,j+1) - (0.5 * cnmask(i,j) * cnmask(i,j+1))
1908 cn_v(i,j) = (wgt1*cn(i,j)) + (wgt2*cn(i,j+1))
1909 enddo ; enddo
1910
1911 ifreq = 1.0 / freq
1912 cn_subro = 1e-30*us%m_s_to_L_T
1913 angle_size = (8.0*atan(1.0)) / (real(nangle))
1914 dt_angle_size = dt / angle_size
1915
1916 do a=asd,aed
1917 angle = (real(a) - 0.5) * angle_size
1918 cos_angle(a) = cos(angle) ; sin_angle(a) = sin(angle)
1919 enddo
1920
1921 !### There should also be refraction due to cn.grad(grid_orientation).
1922 cfl_ang(:,:,:) = 0.0
1923 do j=js,je
1924 ! Copy En into angle space with halos.
1925 do a=1,na ; do i=is,ie
1926 en2d(i,a) = en(i,j,a)
1927 enddo ; enddo
1928 do a=asd,0 ; do i=is,ie
1929 en2d(i,a) = en2d(i,a+nangle)
1930 en2d(i,nangle+stencil+a) = en2d(i,stencil+a)
1931 enddo ; enddo
1932
1933 ! Do the refraction.
1934 do i=is,ie
1935 f2 = 0.25* ((g%Coriolis2Bu(i,j) + g%Coriolis2Bu(i-1,j-1)) + &
1936 (g%Coriolis2Bu(i,j-1) + g%Coriolis2Bu(i-1,j)))
1937 favg = 0.25*((g%CoriolisBu(i,j) + g%CoriolisBu(i-1,j-1)) + &
1938 (g%CoriolisBu(i,j-1) + g%CoriolisBu(i-1,j)))
1939 df_dx = 0.5*g%IdxT(i,j)*((g%CoriolisBu(i,j) - g%CoriolisBu(i-1,j-1)) + &
1940 (g%CoriolisBu(i,j-1) - g%CoriolisBu(i-1,j)))
1941 df_dy = 0.5*g%IdyT(i,j)*((g%CoriolisBu(i,j) - g%CoriolisBu(i-1,j-1)) + &
1942 (g%CoriolisBu(i-1,j) - g%CoriolisBu(i,j-1)))
1943
1944 dlncn_dx = g%IdxT(i,j) * (cn_u(i,j) - cn_u(i-1,j)) / (0.5 * (cn_u(i,j) + cn_u(i-1,j)) + cn_subro)
1945 dlncn_dy = g%IdyT(i,j) * (cn_v(i,j) - cn_v(i,j-1)) / (0.5 * (cn_v(i,j) + cn_v(i,j-1)) + cn_subro)
1946
1947 kmag2 = (freq**2 - f2) / (cn(i,j)**2 + cn_subro**2)
1948 if (kmag2 > 0.0) then
1949 i_kmag = 1.0 / sqrt(kmag2)
1950 dk_dt_kmag(i) = -ifreq * (favg*df_dx + (freq**2 - f2) * dlncn_dx) * i_kmag
1951 dl_dt_kmag(i) = -ifreq * (favg*df_dy + (freq**2 - f2) * dlncn_dy) * i_kmag
1952 else
1953 dk_dt_kmag(i) = 0.0
1954 dl_dt_kmag(i) = 0.0
1955 endif
1956 enddo
1957
1958 ! Determine the energy fluxes in angular orientation space.
1959 do a=asd,aed ; do i=is,ie
1960 cfl_ang(i,j,a) = ((cos_angle(a) * dl_dt_kmag(i)) - (sin_angle(a) * dk_dt_kmag(i))) * dt_angle_size
1961 if (abs(cfl_ang(i,j,a)) > 1.0) then
1962 call mom_error(warning, "refract: CFL exceeds 1.", .true.)
1963 if (cfl_ang(i,j,a) > 1.0) then ; cfl_ang(i,j,a) = 1.0 ; else ; cfl_ang(i,j,a) = -1.0 ; endif
1964 endif
1965 enddo ; enddo
1966
1967 ! Advect in angular space
1968 if (.not.use_ppmang) then
1969 ! Use simple upwind
1970 do a=0,na ; do i=is,ie
1971 if (cfl_ang(i,j,a) > 0.0) then
1972 flux_e(i,a) = cfl_ang(i,j,a) * en2d(i,a)
1973 else
1974 flux_e(i,a) = cfl_ang(i,j,a) * en2d(i,a+1)
1975 endif
1976 enddo ; enddo
1977 else
1978 ! Use PPM
1979 do i=is,ie
1980 call ppm_angular_advect(en2d(i,:),cfl_ang(i,j,:),flux_e(i,:),nangle,dt,stencil)
1981 enddo
1982 endif
1983
1984 ! Update and copy back to En.
1985 do a=1,na ; do i=is,ie
1986 !if (En2d(i,a)+(Flux_E(i,A-1)-Flux_E(i,A)) < 0.0) then ! for debugging
1987 ! call MOM_error(FATAL, "refract: OutFlux>Available")
1988 !endif
1989 en(i,j,a) = en2d(i,a) + (flux_e(i,a-1) - flux_e(i,a))
1990 enddo ; enddo
1991 enddo ! j-loop
1992end subroutine refract
1993
1994!> This subroutine calculates the 1-d flux for advection in angular space using a monotonic
1995!! piecewise parabolic scheme. This needs to be called from within i and j spatial loops.
1996subroutine ppm_angular_advect(En2d, CFL_ang, Flux_En, NAngle, dt, halo_ang)
1997 integer, intent(in) :: NAngle !< The number of wave orientations in the
1998 !! discretized wave energy spectrum [nondim]
1999 real, intent(in) :: dt !< Time increment [T ~> s].
2000 integer, intent(in) :: halo_ang !< The halo size in angular space
2001 real, dimension(1-halo_ang:NAngle+halo_ang), &
2002 intent(in) :: En2d !< The internal gravity wave energy density as a
2003 !! function of angular resolution [H Z2 T-2 ~> m3 s-2 or J m-2].
2004 real, dimension(1-halo_ang:NAngle+halo_ang), &
2005 intent(in) :: CFL_ang !< The CFL number of the energy advection across angles [nondim]
2006 real, dimension(0:NAngle), intent(out) :: Flux_En !< The time integrated internal wave energy flux
2007 !! across angles [H Z2 T-2 ~> m3 s-2 or J m-2].
2008 ! Local variables
2009 real :: flux ! The internal wave energy flux across angles [H Z2 T-3 ~> m3 s-3 or W m-2].
2010 real :: u_ang ! Angular propagation speed [Rad T-1 ~> Rad s-1]
2011 real :: Angle_size ! The size of each orientation wedge in radians [Rad]
2012 real :: I_Angle_size ! The inverse of the orientation wedges [Rad-1]
2013 real :: I_dt ! The inverse of the timestep [T-1 ~> s-1]
2014 real :: aR, aL ! Left and right edge estimates of energy density [H Z2 T-2 rad-1 ~> m3 s-2 rad-1 or J m-2 rad-1]
2015 real :: Ep, Ec, Em ! Mean angular energy density for three successive wedges in angular
2016 ! orientation [H Z2 T-2 rad-1 ~> m3 s-2 rad-1 or J m-2 rad-1]
2017 real :: dA, curv_3 ! Difference and curvature of energy density [H Z2 T-2 rad-1 ~> m3 s-2 rad-1 or J m-2 rad-1]
2018 real, parameter :: oneSixth = 1.0/6.0 ! One sixth [nondim]
2019 integer :: a
2020
2021 i_dt = 1.0 / dt
2022 angle_size = (8.0*atan(1.0)) / (real(nangle))
2023 i_angle_size = 1 / angle_size
2024 flux_en(:) = 0
2025
2026 do a=0,nangle
2027 u_ang = cfl_ang(a)*angle_size*i_dt
2028 if (u_ang >= 0.0) then
2029 ! Implementation of PPM-H3
2030 ! Convert wedge-integrated energy density into angular energy densities for three successive
2031 ! wedges around the source wedge for this flux [H Z2 T-2 rad-1 ~> m3 s-2 rad-1 or J m-2 rad-1].
2032 ep = en2d(a+1)*i_angle_size
2033 ec = en2d(a) *i_angle_size
2034 em = en2d(a-1)*i_angle_size
2035 ! Calculate and bound edge values of energy density.
2036 al = ( 5.*ec + ( 2.*em - ep ) ) * onesixth ! H3 estimate
2037 al = max( min(ec,em), al) ; al = min( max(ec,em), al) ! Bound
2038 ar = ( 5.*ec + ( 2.*ep - em ) ) * onesixth ! H3 estimate
2039 ar = max( min(ec,ep), ar) ; ar = min( max(ec,ep), ar) ! Bound
2040 da = ar - al
2041 if ((ep-ec)*(ec-em) <= 0.) then
2042 al = ec ; ar = ec ! use PCM for local extremum
2043 elseif ( 3.0*da*(2.*ec - (ar + al)) > (da*da) ) then
2044 al = 3.*ec - 2.*ar ! Flatten the profile to move the extremum to the left edge
2045 elseif ( 3.0*da*(2.*ec - (ar + al)) < - (da*da) ) then
2046 ar = 3.*ec - 2.*al ! Flatten the profile to move the extremum to the right edge
2047 endif
2048 curv_3 = (ar + al) - 2.0*ec ! Curvature
2049 ! Calculate angular flux rate [H Z2 T-3 ~> m3 s-3 or W m-2]
2050 flux = u_ang*( ar + cfl_ang(a) * ( 0.5*(al - ar) + curv_3 * (cfl_ang(a) - 1.5) ) )
2051 ! Calculate amount of energy fluxed between wedges [H Z2 T-2 ~> m3 s-2 or J m-2]
2052 flux_en(a) = dt * flux
2053 !Flux_En(A) = (dt * I_Angle_size) * flux
2054 else
2055 ! Implementation of PPM-H3
2056 ! Convert wedge-integrated energy density into angular energy densities for three successive
2057 ! wedges around the source wedge for this flux [H Z2 T-2 rad-1 ~> m3 s-2 rad-1 or J m-2 rad-1].
2058 ep = en2d(a+2)*i_angle_size
2059 ec = en2d(a+1)*i_angle_size
2060 em = en2d(a) *i_angle_size
2061 ! Calculate and bound edge values of energy density.
2062 al = ( 5.*ec + ( 2.*em - ep ) ) * onesixth ! H3 estimate
2063 al = max( min(ec,em), al) ; al = min( max(ec,em), al) ! Bound
2064 ar = ( 5.*ec + ( 2.*ep - em ) ) * onesixth ! H3 estimate
2065 ar = max( min(ec,ep), ar) ; ar = min( max(ec,ep), ar) ! Bound
2066 da = ar - al
2067 if ((ep-ec)*(ec-em) <= 0.) then
2068 al = ec ; ar = ec ! use PCM for local extremum
2069 elseif ( 3.0*da*(2.*ec - (ar + al)) > (da*da) ) then
2070 al = 3.*ec - 2.*ar ! Flatten the profile to move the extremum to the left edge
2071 elseif ( 3.0*da*(2.*ec - (ar + al)) < - (da*da) ) then
2072 ar = 3.*ec - 2.*al ! Flatten the profile to move the extremum to the right edge
2073 endif
2074 curv_3 = (ar + al) - 2.0*ec ! Curvature
2075 ! Calculate angular flux rate [H Z2 T-3 ~> m3 s-3 or W m-2]
2076 ! Note that CFL_ang is negative here, so it looks odd compared with equivalent expressions.
2077 flux = u_ang*( al - cfl_ang(a) * ( 0.5*(ar - al) + curv_3 * (-cfl_ang(a) - 1.5) ) )
2078 ! Calculate amount of energy fluxed between wedges [H Z2 T-2 ~> m3 s-2 or J m-2]
2079 flux_en(a) = dt * flux
2080 !Flux_En(A) = (dt * I_Angle_size) * flux
2081 endif
2082 enddo
2083end subroutine ppm_angular_advect
2084
2085!> Propagates internal waves at a single frequency.
2086subroutine propagate(En, cn, freq, dt, G, GV, US, CS, NAngle, test, halo_size, residual_loss)
2087 type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure.
2088 type(verticalgrid_type), intent(in) :: GV !< The ocean's vertical grid structure.
2089 integer, intent(in) :: NAngle !< The number of wave orientations in the
2090 !! discretized wave energy spectrum.
2091 real, dimension(G%isd:G%ied,G%jsd:G%jed,NAngle), &
2092 intent(inout) :: En !< The internal gravity wave energy density as a
2093 !! function of space and angular resolution,
2094 !! [H Z2 T-2 ~> m3 s-2 or J m-2].
2095 real, dimension(G%isd:G%ied,G%jsd:G%jed), &
2096 intent(in) :: cn !< Baroclinic mode speed [L T-1 ~> m s-1].
2097 real, intent(in) :: freq !< Wave frequency [T-1 ~> s-1].
2098 real, intent(in) :: dt !< Time step [T ~> s].
2099 type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
2100 real, dimension(G%isd:G%ied,G%jsd:G%jed,2), intent(in) :: test !< test rotation vector
2101 type(int_tide_cs), intent(inout) :: CS !< Internal tide control structure
2102 integer, intent(in) :: halo_size !< halo size for correct rotation
2103 real, dimension(G%isd:G%ied,G%jsd:G%jed,NAngle), &
2104 intent(inout) :: residual_loss !< internal tide energy loss due
2105 !! to the residual at slopes [H Z2 T-3 ~> m3 s-3 or W m-2].
2106 ! Local variables
2107 integer, parameter :: stencil = 2
2108 real, dimension(SZIB_(G),SZJ_(G)) :: &
2109 speed_x ! The magnitude of the group velocity at the Cu points [L T-1 ~> m s-1].
2110 real, dimension(SZI_(G),SZJB_(G)) :: &
2111 speed_y ! The magnitude of the group velocity at the Cv points [L T-1 ~> m s-1].
2112 real, dimension(0:NAngle) :: &
2113 cos_angle, sin_angle ! The cosine and sine of each angle [nondim]
2114 real, dimension(NAngle) :: &
2115 Cgx_av, & ! The average projection of the wedge into the x-direction [nondim]
2116 Cgy_av, & ! The average projection of the wedge into the y-direction [nondim]
2117 dCgx, & ! The difference in x-projections between the edges of each angular band [nondim].
2118 dCgy ! The difference in y-projections between the edges of each angular band [nondim].
2119 real :: f2 ! The squared Coriolis parameter [T-2 ~> s-2].
2120 real :: Angle_size ! The size of each wedge of angles [rad]
2121 real :: I_Angle_size ! The inverse of the size of each wedge of angles [rad-1]
2122 real :: angle ! The central angle of each wedge [rad]
2123 real :: Ifreq ! The inverse of the frequency [T ~> s]
2124 real :: freq2 ! The frequency squared [T-2 ~> s-2]
2125 type(loop_bounds_type) :: LB
2126 logical :: x_first
2127 integer :: is, ie, js, je, asd, aed, na
2128 integer :: ish, ieh, jsh, jeh
2129 integer :: i, j, a, fr, m
2130
2131 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; na = size(en,3)
2132 asd = 1-stencil ; aed = nangle+stencil
2133
2134 ifreq = 1.0 / freq
2135 freq2 = freq**2
2136
2137 ! Define loop bounds: Need extensions on j-loop so propagate_y
2138 ! (done after propagate_x) will have updated values in the halo
2139 ! for correct PPM reconstruction. Use if no teleporting and
2140 ! no pass_var between propagate_x and propagate_y.
2141 !jsh = js-3 ; jeh = je+3 ; ish = is ; ieh = ie
2142
2143 ! Define loop bounds: Need 1-pt extensions on loops because
2144 ! teleporting eats up a halo point. Use if teleporting.
2145 ! Also requires pass_var before propagate_y.
2146 jsh = js-1 ; jeh = je+1 ; ish = is-1 ; ieh = ie+1
2147
2148 angle_size = (8.0*atan(1.0)) / real(nangle)
2149 i_angle_size = 1.0 / angle_size
2150
2151 x_first = .true. ! x_first = (MOD(G%first_direction,2) == 0)
2152
2153 if (cs%debug) then
2154 do m=1,cs%nMode ; do fr=1,cs%Nfreq
2155 call sum_en(g, gv, us, cs, cs%En(:,:,:,fr,m), 'propagate: top of routine')
2156 if (is_root_pe()) write(stdout,'(A,E18.10)') 'propagate: top of routine', cs%En_sum
2157 enddo ; enddo
2158 endif
2159
2160 ! IMPLEMENT PPM ADVECTION IN HORIZONTAL-----------------------
2161 ! These could be in the control structure, as they do not vary.
2162 do a=0,na
2163 ! These are the angles at the cell edges...
2164 angle = (real(a) - 0.5) * angle_size
2165 cos_angle(a) = cos(angle) ; sin_angle(a) = sin(angle)
2166 enddo
2167
2168 do a=1,na
2169 cgx_av(a) = (sin_angle(a) - sin_angle(a-1)) * i_angle_size
2170 cgy_av(a) = -(cos_angle(a) - cos_angle(a-1)) * i_angle_size
2171 dcgx(a) = sqrt(0.5 + 0.5*(sin_angle(a)*cos_angle(a) - &
2172 sin_angle(a-1)*cos_angle(a-1)) * i_angle_size - &
2173 cgx_av(a)**2)
2174 dcgy(a) = sqrt(0.5 - 0.5*(sin_angle(a)*cos_angle(a) - &
2175 sin_angle(a-1)*cos_angle(a-1)) * i_angle_size - &
2176 cgy_av(a)**2)
2177 enddo
2178
2179 speed_x(:,:) = 0.
2180 do j=jsh,jeh ; do i=ish-1,ieh
2181 f2 = 0.5 * (g%Coriolis2Bu(i,j) + g%Coriolis2Bu(i,j-1))
2182 speed_x(i,j) = 0.5*(cn(i,j) + cn(i+1,j)) * g%mask2dCu(i,j) * &
2183 sqrt(max(freq2 - f2, 0.0)) * ifreq
2184 enddo ; enddo
2185
2186 speed_y(:,:) = 0.
2187 do j=jsh-1,jeh ; do i=ish,ieh
2188 f2 = 0.5 * (g%Coriolis2Bu(i,j) + g%Coriolis2Bu(i-1,j))
2189 speed_y(i,j) = 0.5*(cn(i,j) + cn(i,j+1)) * g%mask2dCv(i,j) * &
2190 sqrt(max(freq2 - f2, 0.0)) * ifreq
2191 enddo ; enddo
2192
2193 call pass_var(speed_x, g%Domain, position=east_face)
2194 call pass_var(speed_y, g%Domain, position=north_face)
2195
2196 call pass_var(en, g%domain)
2197
2198 ! Apply propagation in the first direction (reflection included)
2199 lb%jsh = jsh ; lb%jeh = jeh ; lb%ish = ish ; lb%ieh = ieh
2200 if (x_first) then
2201 call propagate_x(en, speed_x, cgx_av, dcgx, dt, g, us, cs%nAngle, cs, lb, residual_loss, freq2)
2202 else
2203 call propagate_y(en, speed_y, cgy_av, dcgy, dt, g, us, cs%nAngle, cs, lb, residual_loss, freq2)
2204 endif
2205
2206 ! fix underflows
2207 do a=1,na ; do j=jsh,jeh ; do i=ish,ieh
2208 if (abs(en(i,j,a)) < cs%En_underflow) en(i,j,a) = 0.0
2209 enddo ; enddo ; enddo
2210
2211 if (cs%debug) then
2212 do m=1,cs%nMode ; do fr=1,cs%Nfreq
2213 call sum_en(g, gv, us, cs, cs%En(:,:,:,fr,m), 'propagate: after propagate_x')
2214 if (is_root_pe()) write(stdout,'(A,E18.10)') 'propagate: after propagate_x', cs%En_sum
2215 enddo ; enddo
2216 endif
2217
2218 ! Update halos
2219 call pass_var(en, g%domain)
2220 call correct_halo_rotation_2d(en, test, g, nangle, halo=halo_size)
2221
2222 if (cs%debug) then
2223 do m=1,cs%nMode ; do fr=1,cs%Nfreq
2224 call sum_en(g, gv, us, cs, cs%En(:,:,:,fr,m), 'propagate: after halo update')
2225 if (is_root_pe()) write(stdout,'(A,E18.10)') 'propagate: after halo update', cs%En_sum
2226 enddo ; enddo
2227 endif
2228 ! Apply propagation in the second direction (reflection included)
2229 ! LB%jsh = js ; LB%jeh = je ; LB%ish = is ; LB%ieh = ie ! Use if no teleport
2230 lb%jsh = jsh ; lb%jeh = jeh ; lb%ish = ish ; lb%ieh = ieh
2231 if (x_first) then
2232 call propagate_y(en, speed_y, cgy_av, dcgy, dt, g, us, cs%nAngle, cs, lb, residual_loss, freq2)
2233 else
2234 call propagate_x(en, speed_x, cgx_av, dcgx, dt, g, us, cs%nAngle, cs, lb, residual_loss, freq2)
2235 endif
2236
2237 ! fix underflows
2238 do a=1,na ; do j=jsh,jeh ; do i=ish,ieh
2239 if (abs(en(i,j,a)) < cs%En_underflow) en(i,j,a) = 0.0
2240 enddo ; enddo ; enddo
2241
2242 call pass_var(en, g%domain)
2243 call correct_halo_rotation_2d(en, test, g, nangle, halo=halo_size)
2244
2245 if (cs%debug) then
2246 do m=1,cs%nMode ; do fr=1,cs%Nfreq
2247 call sum_en(g, gv, us, cs, cs%En(:,:,:,fr,m), 'propagate: bottom of routine')
2248 if (is_root_pe()) write(stdout,'(A,E18.10)') 'propagate: bottom of routine', cs%En_sum
2249 enddo ; enddo
2250 endif
2251
2252end subroutine propagate
2253
2254
2255!> Propagates the internal wave energy in the logical x-direction.
2256subroutine propagate_x(En, speed_x, Cgx_av, dCgx, dt, G, US, Nangle, CS, LB, residual_loss, freq2)
2257 type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure.
2258 integer, intent(in) :: NAngle !< The number of wave orientations in the
2259 !! discretized wave energy spectrum.
2260 real, dimension(G%isd:G%ied,G%jsd:G%jed,Nangle), &
2261 intent(inout) :: En !< The energy density integrated over an angular
2262 !! band [H Z2 T-2 ~> m3 s-2 or J m-2].
2263 real, dimension(G%IsdB:G%IedB,G%jsd:G%jed), &
2264 intent(in) :: speed_x !< The magnitude of the group velocity at the
2265 !! Cu points [L T-1 ~> m s-1].
2266 real, dimension(Nangle), intent(in) :: Cgx_av !< The average x-projection in each angular band [nondim]
2267 real, dimension(Nangle), intent(in) :: dCgx !< The difference in x-projections between the
2268 !! edges of each angular band [nondim].
2269 real, intent(in) :: dt !< Time increment [T ~> s].
2270 type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
2271 type(int_tide_cs), intent(in) :: CS !< Internal tide control structure
2272 type(loop_bounds_type), intent(in) :: LB !< A structure with the active energy loop bounds.
2273 real, dimension(G%isd:G%ied,G%jsd:G%jed,Nangle), &
2274 intent(inout) :: residual_loss !< internal tide energy loss due
2275 !! to the residual at slopes [H Z2 T-3 ~> m3 s-3 or W m-2].
2276 real, intent(in) :: freq2 !< The square of internal tides frequency [T-2 ~> s-2].
2277
2278 ! Local variables
2279 real, dimension(SZI_(G),SZJ_(G)) :: &
2280 EnL, EnR ! Left and right face energy densities [H Z2 T-2 ~> m3 s-2 or J m-2].
2281 real, dimension(SZIB_(G),SZJ_(G)) :: &
2282 flux_x ! The internal wave energy flux [H Z2 L2 T-3 ~> m5 s-3 or J s-1].
2283 real, dimension(SZIB_(G)) :: &
2284 cg_p, & ! The x-direction group velocity [L T-1 ~> m s-1]
2285 flux1 ! A 1-d copy of the x-direction internal wave energy flux [H Z2 L2 T-3 ~> m5 s-3 or J s-1].
2286 real, dimension(G%isd:G%ied,G%jsd:G%jed,Nangle) :: &
2287 Fdt_m, Fdt_p! Left and right energy fluxes [H Z2 L2 T-2 ~> m5 s-2 or J]
2288 integer :: i, j, ish, ieh, jsh, jeh, a
2289
2290 ish = lb%ish ; ieh = lb%ieh ; jsh = lb%jsh ; jeh = lb%jeh
2291 do a=1,nangle
2292 ! This sets EnL and EnR.
2293 if (cs%upwind_1st) then
2294 do j=jsh,jeh ; do i=ish-1,ieh+1
2295 enl(i,j) = en(i,j,a) ; enr(i,j) = en(i,j,a)
2296 enddo ; enddo
2297 else
2298 call ppm_reconstruction_x(en(:,:,a), enl, enr, g, lb, &
2299 simple_2nd=cs%simple_2nd, adv_limiter=cs%itides_adv_limiter)
2300 endif
2301
2302 do j=jsh,jeh
2303 ! This is done once with single speed (GARDEN SPRINKLER EFFECT POSSIBLE)
2304 do i=ish-1,ieh
2305 cg_p(i) = speed_x(i,j) * (cgx_av(a))
2306 enddo
2307 call zonal_flux_en(cg_p, en(:,j,a), enl(:,j), enr(:,j), flux1, &
2308 dt, g, us, j, ish, ieh, cs%vol_CFL)
2309 do i=ish-1,ieh ; flux_x(i,j) = flux1(i) ; enddo
2310 enddo
2311
2312 do j=jsh,jeh ; do i=ish,ieh
2313 fdt_m(i,j,a) = dt*flux_x(i-1,j) ! left face influx [H Z2 L2 T-2 ~> m5 s-2 or J]
2314 fdt_p(i,j,a) = -dt*flux_x(i,j) ! right face influx [H Z2 L2 T-2 ~> m5 s-2 or J]
2315
2316 ! only compute residual loss on partial reflection cells, remove numerical noise elsewhere
2317 if (cs%refl_pref_logical(i,j)) then
2318 residual_loss(i,j,a) = residual_loss(i,j,a) + &
2319 ((abs(flux_x(i-1,j)) * cs%residual(i,j) * g%IareaT(i,j)) + &
2320 (abs(flux_x(i,j)) * cs%residual(i,j) * g%IareaT(i,j)))
2321 endif
2322 enddo ; enddo
2323
2324 enddo ! a-loop
2325
2326 ! Only reflect newly arrived energy; existing energy in incident wedge is not reflected
2327 ! and will eventually propagate out of cell. (This code only reflects if En > 0.)
2328 call reflect(fdt_m, nangle, cs, g, lb)
2329 !call teleport(Fdt_m, Nangle, CS, G, LB)
2330 call reflect(fdt_p, nangle, cs, g, lb)
2331 !call teleport(Fdt_p, Nangle, CS, G, LB)
2332
2333 ! Update reflected energy [H Z2 T-2 ~> m3 s-2 or J m-2]
2334 do a=1,nangle ; do j=jsh,jeh ; do i=ish,ieh
2335 en(i,j,a) = en(i,j,a) + (g%IareaT(i,j)*(fdt_m(i,j,a) + fdt_p(i,j,a)))
2336 enddo ; enddo ; enddo
2337
2338 ! existing energy at turning latitude should reflect away
2339 if (cs%turn_critical_lat ) then
2340 call turning_latitude(en, nangle, freq2, cs, g, lb)
2341 endif
2342
2343end subroutine propagate_x
2344
2345!> Propagates the internal wave energy in the logical y-direction.
2346subroutine propagate_y(En, speed_y, Cgy_av, dCgy, dt, G, US, Nangle, CS, LB, residual_loss, freq2)
2347 type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure.
2348 integer, intent(in) :: NAngle !< The number of wave orientations in the
2349 !! discretized wave energy spectrum.
2350 real, dimension(G%isd:G%ied,G%jsd:G%jed,Nangle), &
2351 intent(inout) :: En !< The energy density integrated over an angular
2352 !! band [H Z2 T-2 ~> m3 s-2 or J m-2].
2353 real, dimension(G%isd:G%ied,G%JsdB:G%JedB), &
2354 intent(in) :: speed_y !< The magnitude of the group velocity at the
2355 !! Cv points [L T-1 ~> m s-1].
2356 real, dimension(Nangle), intent(in) :: Cgy_av !< The average y-projection in each angular band [nondim]
2357 real, dimension(Nangle), intent(in) :: dCgy !< The difference in y-projections between the
2358 !! edges of each angular band [nondim]
2359 real, intent(in) :: dt !< Time increment [T ~> s].
2360 type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
2361 type(int_tide_cs), intent(in) :: CS !< Internal tide control structure
2362 type(loop_bounds_type), intent(in) :: LB !< A structure with the active energy loop bounds.
2363 real, dimension(G%isd:G%ied,G%jsd:G%jed,Nangle), &
2364 intent(inout) :: residual_loss !< internal tide energy loss due
2365 !! to the residual at slopes [H Z2 T-3 ~> m3 s-3 or W m-2].
2366 real, intent(in) :: freq2 !< The square of internal tides frequency [T-2 ~> s-2].
2367
2368 ! Local variables
2369 real, dimension(SZI_(G),SZJ_(G)) :: &
2370 EnL, EnR ! South and north face energy densities [H Z2 T-2 ~> m3 s-2 or J m-2].
2371 real, dimension(SZI_(G),SZJB_(G)) :: &
2372 flux_y ! The internal wave energy flux [H Z2 L2 T-3 ~> m5 s-3 or J s-1].
2373 real, dimension(SZI_(G)) :: &
2374 cg_p, & ! The y-direction group velocity [L T-1 ~> m s-1]
2375 flux1 ! A 1-d copy of the y-direction internal wave energy flux [H Z2 L2 T-3 ~> m5 s-3 or J s-1].
2376 real, dimension(G%isd:G%ied,G%jsd:G%jed,Nangle) :: &
2377 Fdt_m, Fdt_p! South and north energy fluxes [H Z2 L2 T-2 ~> m5 s-2 or J]
2378 integer :: i, j, ish, ieh, jsh, jeh, a
2379
2380 ish = lb%ish ; ieh = lb%ieh ; jsh = lb%jsh ; jeh = lb%jeh
2381 do a=1,nangle
2382 ! This sets EnL and EnR.
2383 if (cs%upwind_1st) then
2384 do j=jsh-1,jeh+1 ; do i=ish,ieh
2385 enl(i,j) = en(i,j,a) ; enr(i,j) = en(i,j,a)
2386 enddo ; enddo
2387 else
2388 call ppm_reconstruction_y(en(:,:,a), enl, enr, g, lb, &
2389 simple_2nd=cs%simple_2nd, adv_limiter=cs%itides_adv_limiter)
2390 endif
2391
2392 do j=jsh-1,jeh
2393 ! This is done once with single speed (GARDEN SPRINKLER EFFECT POSSIBLE)
2394 do i=ish,ieh
2395 cg_p(i) = speed_y(i,j) * (cgy_av(a))
2396 enddo
2397 call merid_flux_en(cg_p, en(:,:,a), enl(:,:), enr(:,:), flux1, &
2398 dt, g, us, j, ish, ieh, cs%vol_CFL)
2399 do i=ish,ieh ; flux_y(i,j) = flux1(i) ; enddo
2400 enddo
2401
2402 do j=jsh,jeh ; do i=ish,ieh
2403 fdt_m(i,j,a) = dt*flux_y(i,j-1) ! south face influx [H Z2 L2 T-2 ~> m5 s-2 or J]
2404 fdt_p(i,j,a) = -dt*flux_y(i,j) ! north face influx [H Z2 L2 T-2 ~> m5 s-2 or J]
2405
2406 ! only compute residual loss on partial reflection cells, remove numerical noise elsewhere
2407 if (cs%refl_pref_logical(i,j)) then
2408 residual_loss(i,j,a) = residual_loss(i,j,a) + &
2409 ((abs(flux_y(i,j-1)) * cs%residual(i,j) * g%IareaT(i,j)) + &
2410 (abs(flux_y(i,j)) * cs%residual(i,j) * g%IareaT(i,j)))
2411 endif
2412 enddo ; enddo
2413
2414 enddo ! a-loop
2415
2416 ! Only reflect newly arrived energy; existing energy in incident wedge is not reflected
2417 ! and will eventually propagate out of cell. (This code only reflects if En > 0.)
2418 call reflect(fdt_m, nangle, cs, g, lb)
2419 !call teleport(Fdt_m, Nangle, CS, G, LB)
2420 call reflect(fdt_p, nangle, cs, g, lb)
2421 !call teleport(Fdt_p, Nangle, CS, G, LB)
2422
2423 ! Update reflected energy [H Z2 T-2 ~> m3 s-2 or J m-2]
2424 do a=1,nangle ; do j=jsh,jeh ; do i=ish,ieh
2425 en(i,j,a) = en(i,j,a) + g%IareaT(i,j)*(fdt_m(i,j,a) + fdt_p(i,j,a))
2426 enddo ; enddo ; enddo
2427
2428 ! existing energy at turning latitude should reflect away
2429 if (cs%turn_critical_lat ) then
2430 call turning_latitude(en, nangle, freq2, cs, g, lb)
2431 endif
2432
2433end subroutine propagate_y
2434
2435!> Evaluates the zonal mass or volume fluxes in a layer.
2436subroutine zonal_flux_en(u, h, hL, hR, uh, dt, G, US, j, ish, ieh, vol_CFL)
2437 type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure.
2438 real, dimension(SZIB_(G)), intent(in) :: u !< The zonal velocity [L T-1 ~> m s-1].
2439 real, dimension(SZI_(G)), intent(in) :: h !< Energy density used to calculate the fluxes
2440 !! [H Z2 T-2 ~> m3 s-2 or J m-2].
2441 real, dimension(SZI_(G)), intent(in) :: hL !< Left- Energy densities in the reconstruction
2442 !! [H Z2 T-2 ~> m3 s-2 or J m-2].
2443 real, dimension(SZI_(G)), intent(in) :: hR !< Right- Energy densities in the reconstruction
2444 !! [H Z2 T-2 ~> m3 s-2 or J m-2].
2445 real, dimension(SZIB_(G)), intent(out) :: uh !< The zonal energy transport [H Z2 L2 T-3 ~> m5 s-3 or J s-1].
2446 real, intent(in) :: dt !< Time increment [T ~> s].
2447 type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
2448 integer, intent(in) :: j !< The j-index to work on.
2449 integer, intent(in) :: ish !< The start i-index range to work on.
2450 integer, intent(in) :: ieh !< The end i-index range to work on.
2451 logical, intent(in) :: vol_CFL !< If true, rescale the ratio of face areas to
2452 !! the cell areas when estimating the CFL number.
2453 ! Local variables
2454 real :: CFL ! The CFL number based on the local velocity and grid spacing [nondim].
2455 real :: curv_3 ! A measure of the energy density curvature over a grid length [H Z2 T-2 ~> m3 s-2 or J m-2]
2456 integer :: i
2457
2458 do i=ish-1,ieh
2459 ! Set new values of uh and duhdu.
2460 if (u(i) > 0.0) then
2461 if (vol_cfl) then ; cfl = (u(i) * dt) * (g%dy_Cu(i,j) * g%IareaT(i,j))
2462 else ; cfl = u(i) * dt * g%IdxT(i,j) ; endif
2463 curv_3 = (hl(i) + hr(i)) - 2.0*h(i)
2464 uh(i) = g%dy_Cu(i,j) * u(i) * &
2465 (hr(i) + cfl * (0.5*(hl(i) - hr(i)) + curv_3*(cfl - 1.5)))
2466 elseif (u(i) < 0.0) then
2467 if (vol_cfl) then ; cfl = (-u(i) * dt) * (g%dy_Cu(i,j) * g%IareaT(i+1,j))
2468 else ; cfl = -u(i) * dt * g%IdxT(i+1,j) ; endif
2469 curv_3 = (hl(i+1) + hr(i+1)) - 2.0*h(i+1)
2470 uh(i) = g%dy_Cu(i,j) * u(i) * &
2471 (hl(i+1) + cfl * (0.5*(hr(i+1)-hl(i+1)) + curv_3*(cfl - 1.5)))
2472 else
2473 uh(i) = 0.0
2474 endif
2475 enddo
2476end subroutine zonal_flux_en
2477
2478!> Evaluates the meridional mass or volume fluxes in a layer.
2479subroutine merid_flux_en(v, h, hL, hR, vh, dt, G, US, J, ish, ieh, vol_CFL)
2480 type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure.
2481 real, dimension(SZI_(G)), intent(in) :: v !< The meridional velocity [L T-1 ~> m s-1].
2482 real, dimension(SZI_(G),SZJ_(G)), intent(in) :: h !< Energy density used to calculate the
2483 !! fluxes [H Z2 T-2 ~> m3 s-2 or J m-2].
2484 real, dimension(SZI_(G),SZJ_(G)), intent(in) :: hL !< Left- Energy densities in the
2485 !! reconstruction [H Z2 T-2 ~> m3 s-2 or J m-2].
2486 real, dimension(SZI_(G),SZJ_(G)), intent(in) :: hR !< Right- Energy densities in the
2487 !! reconstruction [H Z2 T-2 ~> m3 s-2 or J m-2].
2488 real, dimension(SZI_(G)), intent(out) :: vh !< The meridional energy transport
2489 !! [H Z2 L2 T-3 ~> m5 s-3 or J s-1].
2490 real, intent(in) :: dt !< Time increment [T ~> s].
2491 type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
2492 integer, intent(in) :: J !< The j-index to work on.
2493 integer, intent(in) :: ish !< The start i-index range to work on.
2494 integer, intent(in) :: ieh !< The end i-index range to work on.
2495 logical, intent(in) :: vol_CFL !< If true, rescale the ratio of face
2496 !! areas to the cell areas when estimating
2497 !! the CFL number.
2498 ! Local variables
2499 real :: CFL ! The CFL number based on the local velocity and grid spacing [nondim].
2500 real :: curv_3 ! A measure of the energy density curvature over a grid length [H Z2 T-2 ~> m3 s-2 or J m-2]
2501 integer :: i
2502
2503 do i=ish,ieh
2504 if (v(i) > 0.0) then
2505 if (vol_cfl) then ; cfl = (v(i) * dt) * (g%dx_Cv(i,j) * g%IareaT(i,j))
2506 else ; cfl = v(i) * dt * g%IdyT(i,j) ; endif
2507 curv_3 = (hl(i,j) + hr(i,j)) - 2.0*h(i,j)
2508 vh(i) = g%dx_Cv(i,j) * v(i) * ( hr(i,j) + cfl * &
2509 (0.5*(hl(i,j) - hr(i,j)) + curv_3*(cfl - 1.5)) )
2510 elseif (v(i) < 0.0) then
2511 if (vol_cfl) then ; cfl = (-v(i) * dt) * (g%dx_Cv(i,j) * g%IareaT(i,j+1))
2512 else ; cfl = -v(i) * dt * g%IdyT(i,j+1) ; endif
2513 curv_3 = (hl(i,j+1) + hr(i,j+1)) - 2.0*h(i,j+1)
2514 vh(i) = g%dx_Cv(i,j) * v(i) * ( hl(i,j+1) + cfl * &
2515 (0.5*(hr(i,j+1)-hl(i,j+1)) + curv_3*(cfl - 1.5)) )
2516 else
2517 vh(i) = 0.0
2518 endif
2519 enddo
2520end subroutine merid_flux_en
2521
2522!> Reflection of the internal waves at a single frequency.
2523subroutine reflect(En, NAngle, CS, G, LB)
2524 type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure
2525 integer, intent(in) :: NAngle !< The number of wave orientations in the
2526 !! discretized wave energy spectrum.
2527 real, dimension(G%isd:G%ied,G%jsd:G%jed,NAngle), &
2528 intent(inout) :: En !< The internal gravity wave energy density as a
2529 !! function of space and angular resolution
2530 !! [H Z2 T-2 ~> m3 s-2 or J m-2].
2531 type(int_tide_cs), intent(in) :: CS !< Internal tide control structure
2532 type(loop_bounds_type), intent(in) :: LB !< A structure with the active energy loop bounds.
2533
2534 ! Local variables
2535 real, dimension(G%isd:G%ied,G%jsd:G%jed) :: angle_c
2536 ! angle of boundary wrt equator [rad]
2537 real, dimension(G%isd:G%ied,G%jsd:G%jed) :: part_refl
2538 ! fraction of wave energy reflected
2539 ! values should collocate with angle_c [nondim]
2540 logical, dimension(G%isd:G%ied,G%jsd:G%jed) :: ridge
2541 ! tags of cells with double reflection
2542 real, dimension(1:Nangle) :: En_reflected ! Energy reflected [H Z2 T-2 ~> m3 s-2 or J m-2].
2543
2544 real :: TwoPi ! 2*pi = 6.2831853... [nondim]
2545 real :: Angle_size ! size of beam wedge [rad]
2546 real :: I_Angle_size ! inverse of size of beam wedge [rad-1]
2547 integer :: angle_wall ! angle-bin of coast/ridge/shelf wrt equator
2548 integer :: angle_wall0 ! angle-bin of coast/ridge/shelf wrt equator
2549 integer :: angle_r ! angle-bin of reflected ray wrt equator
2550 integer :: angle_r0 ! angle-bin of reflected ray wrt equator
2551 integer :: angle_to_wall ! angle-bin relative to wall
2552 integer :: a, a0 ! loop index for angles
2553 integer :: i, j
2554 integer :: Nangle_d2 ! Nangle / 2
2555 integer :: isc, iec, jsc, jec ! start and end local indices on PE
2556 ! (values exclude halos)
2557 integer :: ish, ieh, jsh, jeh ! start and end local indices on data domain
2558 ! leaving out outdated halo points (march in)
2559
2560 isc = g%isc ; iec = g%iec ; jsc = g%jsc ; jec = g%jec
2561 ish = lb%ish ; ieh = lb%ieh ; jsh = lb%jsh ; jeh = lb%jeh
2562
2563 twopi = 8.0*atan(1.0)
2564 angle_size = twopi / (real(nangle))
2565 i_angle_size = 1.0 / angle_size
2566 nangle_d2 = (nangle / 2)
2567
2568 ! init local arrays
2569 angle_c(:,:) = cs%nullangle
2570 part_refl(:,:) = 0.
2571 ridge(:,:) = .false.
2572
2573 do j=jsh,jeh ; do i=ish,ieh
2574 if (cs%refl_angle(i,j) /= cs%nullangle) then
2575 angle_c(i,j) = mod(cs%refl_angle(i,j) + twopi, twopi)
2576 endif
2577 part_refl(i,j) = cs%refl_pref(i,j)
2578 ridge(i,j) = cs%refl_dbl(i,j)
2579 enddo ; enddo
2580 en_reflected(:) = 0.0
2581
2582 do j=jsh,jeh ; do i=ish,ieh
2583 ! redistribute energy in angular space if ray will hit boundary
2584 ! i.e., if energy is in a reflecting cell
2585 if (angle_c(i,j) /= cs%nullangle) then
2586 ! refection angle is given in rad, convert to the discrete angle
2587 angle_wall = nint(angle_c(i,j)*i_angle_size) + 1
2588 do a=1,nangle ; if (en(i,j,a) > 0.0) then
2589 ! reindex to 0 -> Nangle-1 for trig
2590 a0 = a - 1
2591 angle_wall0 = angle_wall - 1
2592 ! compute relative angle from wall and use cyclic properties
2593 ! to ensure it is bounded by 0 -> Nangle-1
2594 angle_to_wall = mod((a0 - angle_wall0) + nangle, nangle)
2595
2596 if (ridge(i,j)) then
2597 ! if ray is not incident but in ridge cell, use complementary angle
2598 if ((nangle_d2 < angle_to_wall) .and. (angle_to_wall < nangle)) then
2599 angle_wall0 = mod(angle_wall0 + (nangle_d2 + nangle), nangle)
2600 endif
2601 endif
2602
2603 ! do reflection
2604 if ((0 < angle_to_wall) .and. (angle_to_wall < nangle_d2)) then
2605 angle_r0 = mod(2*angle_wall0 - a0 + nangle, nangle)
2606 angle_r = angle_r0 + 1 !re-index to 1 -> Nangle
2607 if (a /= angle_r) then
2608 en_reflected(angle_r) = part_refl(i,j)*en(i,j,a)
2609 en(i,j,a) = (1.0-part_refl(i,j))*en(i,j,a)
2610 endif
2611 endif
2612 endif ; enddo ! a-loop
2613 do a=1,nangle
2614 en(i,j,a) = en(i,j,a) + en_reflected(a)
2615 en_reflected(a) = 0.0 ! reset values
2616 enddo ! a-loop
2617 endif
2618 enddo ; enddo ! i- and j-loops
2619
2620 ! Check to make sure no energy gets onto land (only run for debugging)
2621 ! do a=1,NAngle ; do j=jsc,jec ; do i=isc,iec
2622 ! if (En(i,j,a) > 0.001 .and. G%mask2dT(i,j) == 0) then
2623 ! write (mesg,*) 'En=', HZ2_T2_to_J_m2*En(i,j,a), 'a=', a, 'ig_g=',i+G%idg_offset, 'jg_g=',j+G%jdg_offset
2624 ! call MOM_error(FATAL, "reflect: Energy detected out of bounds: "//trim(mesg), .true.)
2625 ! endif
2626 ! enddo ; enddo ; enddo
2627
2628end subroutine reflect
2629
2630subroutine turning_latitude(En, NAngle, freq2, CS, G, LB)
2631 type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure
2632 integer, intent(in) :: NAngle !< The number of wave orientations in the
2633 !! discretized wave energy spectrum.
2634 real, dimension(G%isd:G%ied,G%jsd:G%jed,NAngle), &
2635 intent(inout) :: En !< The internal gravity wave energy density as a
2636 !! function of space and angular resolution
2637 !! [H Z2 T-2 ~> m3 s-2 or J m-2].
2638 type(int_tide_cs), intent(in) :: CS !< Internal tide control structure
2639 type(loop_bounds_type), intent(in) :: LB !< A structure with the active energy loop bounds.
2640 real, intent(in) :: freq2 !< The square of the internal tide frequency [T-2 ~> s-2]
2641
2642 ! Local variables
2643 real, dimension(G%isd:G%ied,G%jsd:G%jed) :: angle_c
2644 ! angle of boundary wrt equator [rad]
2645 real, dimension(1:Nangle) :: En_reflected ! Energy reflected [H Z2 T-2 ~> m3 s-2 or J m-2].
2646
2647 real :: TwoPi ! 2*pi = 6.2831853... [nondim]
2648 real :: Angle_size ! size of beam wedge [rad]
2649 real :: I_Angle_size ! inverse of size of beam wedge [rad-1]
2650 real :: f2
2651
2652 integer :: angle_wall ! angle-bin of coast/ridge/shelf wrt equator
2653 integer :: angle_wall0 ! angle-bin of coast/ridge/shelf wrt equator
2654 integer :: angle_r ! angle-bin of reflected ray wrt equator
2655 integer :: angle_r0 ! angle-bin of reflected ray wrt equator
2656 integer :: angle_to_wall ! angle-bin relative to wall
2657 integer :: a, a0 ! loop index for angles
2658 integer :: i, j
2659 integer :: Nangle_d2 ! Nangle / 2
2660 integer :: Nangle_d4p1 ! Nangle / 4 + 1
2661 integer :: Nangle_3d4p1 ! 3*Nangle / 4 + 1
2662 integer :: isc, iec, jsc, jec ! start and end local indices on PE
2663 ! (values exclude halos)
2664 integer :: ish, ieh, jsh, jeh ! start and end local indices on data domain
2665 ! leaving out outdated halo points (march in)
2666
2667 isc = g%isc ; iec = g%iec ; jsc = g%jsc ; jec = g%jec
2668 ish = lb%ish ; ieh = lb%ieh ; jsh = lb%jsh ; jeh = lb%jeh
2669
2670 twopi = 8.0*atan(1.0)
2671 angle_size = twopi / (real(nangle))
2672 i_angle_size = 1.0 / angle_size
2673 nangle_d2 = (nangle / 2)
2674 nangle_d4p1 = (nangle / 4) + 1
2675 nangle_3d4p1 = (3 * nangle / 4) + 1
2676
2677
2678 ! init local arrays
2679 angle_c(:,:) = cs%nullangle
2680 angle_wall = 0
2681 angle_wall0 =0
2682 angle_r = 0
2683 angle_r0 = 0
2684 angle_to_wall = 0
2685
2686 do j=jsh,jeh ; do i=ish,ieh
2687 ! init
2688 angle_wall = 0
2689 angle_wall0 = 0
2690 angle_r = 0
2691 angle_r0 = 0
2692 angle_to_wall = 0
2693
2694 f2 = max(abs(g%Coriolis2Bu(i-1,j)), abs(g%Coriolis2Bu(i,j)), &
2695 abs(g%Coriolis2Bu(i-1,j-1)), abs(g%Coriolis2Bu(i,j-1)))
2696
2697 if (g%CoriolisBu(i,j) < 0. ) then
2698 if (f2 - freq2 >= 0.) then
2699 angle_c(i,j) = 0.5 * twopi
2700 endif
2701 else
2702 if (f2 - freq2 >= 0.) then
2703 angle_c(i,j) = 0.
2704 endif
2705 endif
2706 enddo ; enddo
2707
2708 en_reflected(:) = 0.0
2709
2710 do j=jsh,jeh ; do i=ish,ieh
2711 ! init
2712 angle_wall = 0
2713 angle_wall0 = 0
2714 angle_r = 0
2715 angle_r0 = 0
2716 angle_to_wall = 0
2717
2718 if (angle_c(i,j) /= cs%nullangle) then
2719 ! refection angle is given in rad, convert to the discrete angle
2720 angle_wall = nint(angle_c(i,j)*i_angle_size) + 1
2721 do a=1,nangle ; if (en(i,j,a) > 0.0) then
2722
2723 if (.not. cs%reflect_critical_lat) then
2724
2725 ! turn parallel to critical lat
2726 if ((a > nangle_d4p1) .and. (a < nangle_3d4p1)) then
2727 angle_r0 = nangle_d2
2728 else
2729 angle_r0 = 0
2730 endif
2731 angle_r = angle_r0 + 1 !re-index to 1 -> Nangle
2732
2733 if (a /= angle_r) then
2734 en_reflected(angle_r) = en(i,j,a)
2735 en(i,j,a) = 0.
2736 endif
2737
2738 else
2739
2740 ! reindex to 0 -> Nangle-1 for trig
2741 a0 = a - 1
2742 angle_wall0 = angle_wall - 1
2743 ! compute relative angle from wall and use cyclic properties
2744 ! to ensure it is bounded by 0 -> Nangle-1
2745 angle_to_wall = mod((a0 - angle_wall0) + nangle, nangle)
2746
2747 ! do reflection
2748 if ((0 < angle_to_wall) .and. (angle_to_wall < nangle_d2)) then
2749 angle_r0 = mod(2*angle_wall0 - a0 + nangle, nangle)
2750 angle_r = angle_r0 + 1 !re-index to 1 -> Nangle
2751
2752 if (a /= angle_r) then
2753 en_reflected(angle_r) = en(i,j,a)
2754 en(i,j,a) = 0.
2755 endif
2756 endif
2757 endif
2758 endif ; enddo ! a-loop
2759
2760 do a=1,nangle
2761 en(i,j,a) = en(i,j,a) + en_reflected(a)
2762 en_reflected(a) = 0.0 ! reset values
2763 enddo ! a-loop
2764 endif
2765 enddo ; enddo ! i- and j-loops
2766
2767end subroutine turning_latitude
2768
2769!> Moves energy across lines of partial reflection to prevent
2770!! reflection of energy that is supposed to get across.
2771subroutine teleport(En, NAngle, CS, G, LB)
2772 type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure.
2773 integer, intent(in) :: NAngle !< The number of wave orientations in the
2774 !! discretized wave energy spectrum.
2775 real, dimension(G%isd:G%ied,G%jsd:G%jed,NAngle), &
2776 intent(inout) :: En !< The internal gravity wave energy density as a
2777 !! function of space and angular resolution
2778 !! [H Z2 T-2 ~> m3 s-2 or J m-2].
2779 type(int_tide_cs), intent(in) :: CS !< Internal tide control structure
2780 type(loop_bounds_type), intent(in) :: LB !< A structure with the active energy loop bounds.
2781 ! Local variables
2782 real, dimension(G%isd:G%ied,G%jsd:G%jed) :: angle_c
2783 ! angle of boundary wrt equator [rad]
2784 real, dimension(G%isd:G%ied,G%jsd:G%jed) :: part_refl
2785 ! fraction of wave energy reflected
2786 ! values should collocate with angle_c [nondim]
2787 logical, dimension(G%isd:G%ied,G%jsd:G%jed) :: pref_cell
2788 ! flag for partial reflection
2789 logical, dimension(G%isd:G%ied,G%jsd:G%jed) :: ridge
2790 ! tags of cells with double reflection
2791 real :: TwoPi ! 2*pi = 6.2831853... [nondim]
2792 real :: Angle_size ! size of beam wedge [rad]
2793 real, dimension(1:NAngle) :: angle_i ! angle of incident ray wrt equator [rad]
2794 real, dimension(1:NAngle) :: cos_angle ! Cosine of the beam angle relative to eastward [nondim]
2795 real, dimension(1:NAngle) :: sin_angle ! Sine of the beam angle relative to eastward [nondim]
2796 real :: En_tele ! energy to be "teleported" [H Z2 T-2 ~> m3 s-2 or J m-2]
2797 character(len=160) :: mesg ! The text of an error message
2798 integer :: i, j, a
2799 integer :: ish, ieh, jsh, jeh ! start and end local indices on data domain
2800 ! leaving out outdated halo points (march in)
2801 integer :: id_g, jd_g ! global (decomposition-invariant) indices
2802 integer :: jos, ios ! offsets
2803 real :: cos_normal, sin_normal ! cos/sin of cross-ridge normal direction [nondim]
2804 real :: angle_wall ! The coastline angle or the complementary angle [radians]
2805
2806 ish = lb%ish ; ieh = lb%ieh ; jsh = lb%jsh ; jeh = lb%jeh
2807
2808 twopi = 8.0*atan(1.0)
2809 angle_size = twopi / (real(nangle))
2810
2811 do a=1,nangle
2812 ! These are the angles at the cell centers
2813 ! (should do this elsewhere since doesn't change with time)
2814 angle_i(a) = angle_size * real(a - 1) ! for a=1 aligned with x-axis
2815 cos_angle(a) = cos(angle_i(a)) ; sin_angle(a) = sin(angle_i(a))
2816 enddo
2817
2818 angle_c = cs%refl_angle
2819 part_refl = cs%refl_pref
2820 pref_cell = cs%refl_pref_logical
2821 ridge = cs%refl_dbl
2822
2823 do j=jsh,jeh
2824 do i=ish,ieh
2825 id_g = i + g%idg_offset ; jd_g = j + g%jdg_offset
2826 if (pref_cell(i,j)) then
2827 do a=1,nangle
2828 if (en(i,j,a) > 0) then
2829 ! if ray is incident, keep specified boundary angle
2830 if (sin(angle_i(a) - angle_c(i,j)) >= 0.0) then
2831 angle_wall = angle_c(i,j)
2832 ! if ray is not incident but in ridge cell, use complementary angle
2833 elseif (ridge(i,j)) then
2834 angle_wall = angle_c(i,j) + 0.5*twopi
2835 ! if ray is not incident and not in a ridge cell, keep specified angle
2836 else
2837 angle_wall = angle_c(i,j)
2838 endif
2839 ! teleport if incident
2840 if (sin(angle_i(a) - angle_wall) >= 0.0) then
2841 en_tele = en(i,j,a)
2842 cos_normal = cos(angle_wall + 0.25*twopi)
2843 sin_normal = sin(angle_wall + 0.25*twopi)
2844 ! find preferred zonal offset based on shelf/ridge angle
2845 ios = int(sign(1.,cos_normal))
2846 ! find preferred meridional offset based on shelf/ridge angle
2847 jos = int(sign(1.,sin_normal))
2848 ! find receptive ocean cell in direction of offset
2849 if (.not. pref_cell(i+ios,j+jos)) then
2850 en(i,j,a) = en(i,j,a) - en_tele
2851 en(i+ios,j+jos,a) = en(i+ios,j+jos,a) + en_tele
2852 else
2853 write(mesg,*) 'idg=',id_g,'jd_g=',jd_g,'a=',a
2854 call mom_error(fatal, "teleport: no receptive ocean cell at "//trim(mesg), .true.)
2855 endif
2856 endif ! incidence check
2857 endif ! energy check
2858 enddo ! a-loop
2859 endif ! pref check
2860 enddo ! i-loop
2861 enddo ! j-loop
2862
2863end subroutine teleport
2864
2865!> Rotates points in the halos where required to accommodate
2866!! changes in grid orientation, such as at the tripolar fold.
2867subroutine correct_halo_rotation(En, test, G, NAngle, halo)
2868 type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure
2869 real, dimension(:,:,:,:,:), intent(inout) :: En !< The internal gravity wave energy density as a
2870 !! function of space, angular orientation, frequency,
2871 !! and vertical mode [H Z2 T-2 ~> m3 s-2 or J m-2].
2872 real, dimension(SZI_(G),SZJ_(G),2), &
2873 intent(in) :: test !< An x-unit vector that has been passed through
2874 !! the halo updates, to enable the rotation of the
2875 !! wave energies in the halo region to be corrected [nondim].
2876 integer, intent(in) :: NAngle !< The number of wave orientations in the
2877 !! discretized wave energy spectrum.
2878 integer, intent(in) :: halo !< The halo size over which to do the calculations
2879 ! Local variables
2880 real, dimension(G%isd:G%ied,NAngle) :: En2d ! A zonal row of the internal gravity wave energy density
2881 ! in a frequency band and mode [H Z2 T-2 ~> m3 s-2 or J m-2].
2882 integer, dimension(G%isd:G%ied) :: a_shift
2883 integer :: i_first, i_last, a_new
2884 integer :: a, i, j, ish, ieh, jsh, jeh, m, fr
2885 character(len=160) :: mesg ! The text of an error message
2886 ish = g%isc-halo ; ieh = g%iec+halo ; jsh = g%jsc-halo ; jeh = g%jec+halo
2887
2888 do j=jsh,jeh
2889 i_first = ieh+1 ; i_last = ish-1
2890 do i=ish,ieh
2891 a_shift(i) = 0
2892 if (test(i,j,2) < 0.5) then
2893 if (i<i_first) i_first = i
2894 if (i>i_last) i_last = i
2895 if (test(i,j,2) < -0.5) then ; a_shift(i) = 0.5*nangle
2896 elseif (test(i,j,1) > 0.5) then ; a_shift(i) = -0.25*nangle
2897 elseif (test(i,j,1) < -0.5) then ; a_shift(i) = 0.25*nangle
2898 else
2899 write(mesg,'("Unrecognized rotation test vector ",2ES9.2," at ",F7.2," E, ",&
2900 &F7.2," N; i,j=",2i4)') &
2901 test(i,j,1), test(i,j,2), g%GeoLonT(i,j), g%GeoLatT(i,j), i, j
2902 call mom_error(fatal, mesg)
2903 endif
2904 endif
2905 enddo
2906
2907 if (i_first <= i_last) then
2908 ! At least one point in this row needs to be rotated.
2909 do m=1,size(en,5) ; do fr=1,size(en,4)
2910 do a=1,nangle ; do i=i_first,i_last ; if (a_shift(i) /= 0) then
2911 a_new = a + a_shift(i)
2912 if (a_new < 1) a_new = a_new + nangle
2913 if (a_new > nangle) a_new = a_new - nangle
2914 en2d(i,a_new) = en(i,j,a,fr,m)
2915 endif ; enddo ; enddo
2916 do a=1,nangle ; do i=i_first,i_last ; if (a_shift(i) /= 0) then
2917 en(i,j,a,fr,m) = en2d(i,a)
2918 endif ; enddo ; enddo
2919 enddo ; enddo
2920 endif
2921 enddo
2922end subroutine correct_halo_rotation
2923
2924
2925!> Rotates points in the halos where required to accommodate
2926!! changes in grid orientation, such as at the tripolar fold.
2927subroutine correct_halo_rotation_2d(En, test, G, NAngle, halo)
2928 type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure
2929 real, dimension(:,:,:), intent(inout) :: En !< The internal gravity wave energy density as a
2930 !! function of space, angular orientation, frequency,
2931 !! and vertical mode [H Z2 T-2 ~> m3 s-2 or J m-2].
2932 real, dimension(SZI_(G),SZJ_(G),2), &
2933 intent(in) :: test !< An x-unit vector that has been passed through
2934 !! the halo updates, to enable the rotation of the
2935 !! wave energies in the halo region to be corrected [nondim].
2936 integer, intent(in) :: NAngle !< The number of wave orientations in the
2937 !! discretized wave energy spectrum.
2938 integer, intent(in) :: halo !< The halo size over which to do the calculations
2939 ! Local variables
2940 real, dimension(G%isd:G%ied,NAngle) :: En2d ! A zonal row of the internal gravity wave energy density
2941 ! in a frequency band and mode [H Z2 T-2 ~> m3 s-2 or J m-2].
2942 integer, dimension(G%isd:G%ied) :: a_shift
2943 integer :: i_first, i_last, a_new
2944 integer :: a, i, j, ish, ieh, jsh, jeh
2945 integer :: id_g, jd_g
2946 character(len=160) :: mesg ! The text of an error message
2947 ish = g%isc-halo ; ieh = g%iec+halo ; jsh = g%jsc-halo ; jeh = g%jec+halo
2948
2949 ! top rows
2950 do j=jsh,jeh
2951 !do j= G%jec+1,jeh
2952 i_first = ieh+1 ; i_last = ish-1 ! init
2953 do i=ish,ieh
2954 id_g = i + g%idg_offset ; jd_g = j + g%jdg_offset ! for debugging
2955
2956 a_shift(i) = 0
2957 if (test(i,j,2) < 0.5) then
2958 if (i<i_first) i_first = i
2959 if (i>i_last) i_last = i
2960 if (test(i,j,2) < -0.5) then ; a_shift(i) = 0.5*nangle
2961 elseif (test(i,j,1) > 0.5) then ; a_shift(i) = -0.25*nangle
2962 elseif (test(i,j,1) < -0.5) then ; a_shift(i) = 0.25*nangle
2963 else
2964 write(mesg,'("Unrecognized rotation test vector ",2ES9.2," at ",F7.2," E, ",&
2965 &F7.2," N; i,j=",2i4)') &
2966 test(i,j,1), test(i,j,2), g%GeoLonT(i,j), g%GeoLatT(i,j), i, j
2967 call mom_error(fatal, mesg)
2968 endif
2969 endif
2970 enddo
2971
2972 if (i_first <= i_last) then
2973 ! At least one point in this row needs to be rotated.
2974 do a=1,nangle ; do i=i_first,i_last ; if (a_shift(i) /= 0) then
2975 a_new = a + a_shift(i)
2976 if (a_new < 1) a_new = a_new + nangle
2977 if (a_new > nangle) a_new = a_new - nangle
2978 en2d(i,a_new) = en(i,j,a)
2979 endif ; enddo ; enddo
2980 do a=1,nangle ; do i=i_first,i_last ; if (a_shift(i) /= 0) then
2981 en(i,j,a) = en2d(i,a)
2982 endif ; enddo ; enddo
2983 endif
2984 enddo
2985end subroutine correct_halo_rotation_2d
2986
2987
2988!> Calculates left/right edge values for PPM reconstruction in x-direction.
2989subroutine ppm_reconstruction_x(h_in, h_l, h_r, G, LB, simple_2nd, adv_limiter)
2990 type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure.
2991 real, dimension(SZI_(G),SZJ_(G)), intent(in) :: h_in !< Energy density in a sector (2D)
2992 !! [H Z2 T-2 ~> m3 s-2 or J m-2]
2993 real, dimension(SZI_(G),SZJ_(G)), intent(out) :: h_l !< Left edge value of reconstruction (2D)
2994 !! [H Z2 T-2 ~> m3 s-2 or J m-2]
2995 real, dimension(SZI_(G),SZJ_(G)), intent(out) :: h_r !< Right edge value of reconstruction (2D)
2996 !! [H Z2 T-2 ~> m3 s-2 or J m-2]
2997 type(loop_bounds_type), intent(in) :: LB !< A structure with the active loop bounds.
2998 logical, intent(in) :: simple_2nd !< If true, use the arithmetic mean
2999 !! energy densities as default edge values
3000 !! for a simple 2nd order scheme.
3001 integer, intent(in) :: adv_limiter !< The type of limiter used
3002
3003 ! Local variables
3004 real, dimension(SZI_(G),SZJ_(G)) :: slp ! The slope in energy density times the cell width
3005 ! [H Z2 T-2 ~> m3 s-2 or J m-2]
3006 real, parameter :: oneSixth = 1./6. ! One sixth [nondim]
3007 real :: h_ip1, h_im1 ! The energy densities at adjacent points [H Z2 T-2 ~> m3 s-2 or J m-2]
3008 real :: dMx, dMn ! The maximum and minimum of values of energy density at adjacent points
3009 ! relative to the center point [H Z2 T-2 ~> m3 s-2 or J m-2]
3010 character(len=256) :: mesg ! The text of an error message
3011 integer :: i, j, isl, iel, jsl, jel, stencil
3012
3013 isl = lb%ish-1 ; iel = lb%ieh+1 ; jsl = lb%jsh ; jel = lb%jeh
3014
3015 ! This is the stencil of the reconstruction, not the scheme overall.
3016 stencil = 2 ; if (simple_2nd) stencil = 1
3017
3018 if ((isl-stencil < g%isd) .or. (iel+stencil > g%ied)) then
3019 write(mesg,'("In MOM_internal_tides, PPM_reconstruction_x called with a ", &
3020 & "x-halo that needs to be increased by ",I0,".")') &
3021 stencil + max(g%isd-isl,iel-g%ied)
3022 call mom_error(fatal,mesg)
3023 endif
3024 if ((jsl < g%jsd) .or. (jel > g%jed)) then
3025 write(mesg,'("In MOM_internal_tides, PPM_reconstruction_x called with a ", &
3026 & "y-halo that needs to be increased by ",I0,".")') &
3027 max(g%jsd-jsl,jel-g%jed)
3028 call mom_error(fatal,mesg)
3029 endif
3030
3031 if (simple_2nd) then
3032 do j=jsl,jel ; do i=isl,iel
3033 h_im1 = g%mask2dT(i-1,j) * h_in(i-1,j) + (1.0-g%mask2dT(i-1,j)) * h_in(i,j)
3034 h_ip1 = g%mask2dT(i+1,j) * h_in(i+1,j) + (1.0-g%mask2dT(i+1,j)) * h_in(i,j)
3035 h_l(i,j) = 0.5*( h_im1 + h_in(i,j) )
3036 h_r(i,j) = 0.5*( h_ip1 + h_in(i,j) )
3037 enddo ; enddo
3038 else
3039 do j=jsl,jel ; do i=isl-1,iel+1
3040 if ((g%mask2dT(i-1,j) * g%mask2dT(i,j) * g%mask2dT(i+1,j)) == 0.0) then
3041 slp(i,j) = 0.0
3042 else
3043 ! This uses a simple 2nd order slope.
3044 slp(i,j) = 0.5 * (h_in(i+1,j) - h_in(i-1,j))
3045 ! Monotonic constraint, see Eq. B2 in Lin 1994, MWR (132)
3046 dmx = max(h_in(i+1,j), h_in(i-1,j), h_in(i,j)) - h_in(i,j)
3047 dmn = h_in(i,j) - min(h_in(i+1,j), h_in(i-1,j), h_in(i,j))
3048 slp(i,j) = sign(1.,slp(i,j)) * min(abs(slp(i,j)), 2. * min(dmx, dmn))
3049 ! * (G%mask2dT(i-1,j) * G%mask2dT(i,j) * G%mask2dT(i+1,j))
3050 endif
3051 enddo ; enddo
3052
3053 do j=jsl,jel ; do i=isl,iel
3054 ! Neighboring values should take into account any boundaries. The 3
3055 ! following sets of expressions are equivalent.
3056 ! h_im1 = h_in(i-1,j,k) ; if (G%mask2dT(i-1,j) < 0.5) h_im1 = h_in(i,j)
3057 ! h_ip1 = h_in(i+1,j,k) ; if (G%mask2dT(i+1,j) < 0.5) h_ip1 = h_in(i,j)
3058 h_im1 = g%mask2dT(i-1,j) * h_in(i-1,j) + (1.0-g%mask2dT(i-1,j)) * h_in(i,j)
3059 h_ip1 = g%mask2dT(i+1,j) * h_in(i+1,j) + (1.0-g%mask2dT(i+1,j)) * h_in(i,j)
3060 ! Left/right values following Eq. B2 in Lin 1994, MWR (132)
3061 h_l(i,j) = 0.5*( h_im1 + h_in(i,j) ) + onesixth*( slp(i-1,j) - slp(i,j) )
3062 h_r(i,j) = 0.5*( h_ip1 + h_in(i,j) ) + onesixth*( slp(i,j) - slp(i+1,j) )
3063 enddo ; enddo
3064 endif
3065
3066 select case(adv_limiter)
3067 case (limiter_adv_positive)
3068 call ppm_limit_pos(h_in, h_l, h_r, 0.0, g, isl, iel, jsl, jel)
3069 case (limiter_adv_minmod)
3070 call minmod_limiter(h_in, h_l, h_r, g, isl, iel, jsl, jel)
3071 end select
3072
3073end subroutine ppm_reconstruction_x
3074
3075!> Calculates left/right edge valus for PPM reconstruction in y-direction.
3076subroutine ppm_reconstruction_y(h_in, h_l, h_r, G, LB, simple_2nd, adv_limiter)
3077 type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure.
3078 real, dimension(SZI_(G),SZJ_(G)), intent(in) :: h_in !< Energy density in a sector (2D)
3079 !! [H Z2 T-2 ~> m3 s-2 or J m-2]
3080 real, dimension(SZI_(G),SZJ_(G)), intent(out) :: h_l !< Left edge value of reconstruction (2D)
3081 !! [H Z2 T-2 ~> m3 s-2 or J m-2]
3082 real, dimension(SZI_(G),SZJ_(G)), intent(out) :: h_r !< Right edge value of reconstruction (2D)
3083 !! [H Z2 T-2 ~> m3 s-2 or J m-2]
3084 type(loop_bounds_type), intent(in) :: LB !< A structure with the active loop bounds.
3085 logical, intent(in) :: simple_2nd !< If true, use the arithmetic mean
3086 !! energy densities as default edge values
3087 !! for a simple 2nd order scheme.
3088 integer, intent(in) :: adv_limiter !< The type of limiter used
3089
3090 ! Local variables
3091 real, dimension(SZI_(G),SZJ_(G)) :: slp ! The slope in energy density times the cell width
3092 ! [H Z2 T-2 ~> m3 s-2 or J m-2]
3093 real, parameter :: oneSixth = 1./6. ! One sixth [nondim]
3094 real :: h_jp1, h_jm1 ! The energy densities at adjacent points [H Z2 T-2 ~> m3 s-2 or J m-2]
3095 real :: dMx, dMn ! The maximum and minimum of values of energy density at adjacent points
3096 ! relative to the center point [H Z2 T-2 ~> m3 s-2 or J m-2]
3097 character(len=256) :: mesg ! The text of an error message
3098 integer :: i, j, isl, iel, jsl, jel, stencil
3099
3100 isl = lb%ish ; iel = lb%ieh ; jsl = lb%jsh-1 ; jel = lb%jeh+1
3101
3102 ! This is the stencil of the reconstruction, not the scheme overall.
3103 stencil = 2 ; if (simple_2nd) stencil = 1
3104
3105 if ((isl < g%isd) .or. (iel > g%ied)) then
3106 write(mesg,'("In MOM_internal_tides, PPM_reconstruction_y called with a ", &
3107 & "x-halo that needs to be increased by ",I0,".")') &
3108 max(g%isd-isl,iel-g%ied)
3109 call mom_error(fatal,mesg)
3110 endif
3111 if ((jsl-stencil < g%jsd) .or. (jel+stencil > g%jed)) then
3112 write(mesg,'("In MOM_internal_tides, PPM_reconstruction_y called with a ", &
3113 & "y-halo that needs to be increased by ",I0,".")') &
3114 stencil + max(g%jsd-jsl,jel-g%jed)
3115 call mom_error(fatal,mesg)
3116 endif
3117
3118 if (simple_2nd) then
3119 do j=jsl,jel ; do i=isl,iel
3120 h_jm1 = g%mask2dT(i,j-1) * h_in(i,j-1) + (1.0-g%mask2dT(i,j-1)) * h_in(i,j)
3121 h_jp1 = g%mask2dT(i,j+1) * h_in(i,j+1) + (1.0-g%mask2dT(i,j+1)) * h_in(i,j)
3122 h_l(i,j) = 0.5*( h_jm1 + h_in(i,j) )
3123 h_r(i,j) = 0.5*( h_jp1 + h_in(i,j) )
3124 enddo ; enddo
3125 else
3126 do j=jsl-1,jel+1 ; do i=isl,iel
3127 if ((g%mask2dT(i,j-1) * g%mask2dT(i,j) * g%mask2dT(i,j+1)) == 0.0) then
3128 slp(i,j) = 0.0
3129 else
3130 ! This uses a simple 2nd order slope.
3131 slp(i,j) = 0.5 * (h_in(i,j+1) - h_in(i,j-1))
3132 ! Monotonic constraint, see Eq. B2 in Lin 1994, MWR (132)
3133 dmx = max(h_in(i,j+1), h_in(i,j-1), h_in(i,j)) - h_in(i,j)
3134 dmn = h_in(i,j) - min(h_in(i,j+1), h_in(i,j-1), h_in(i,j))
3135 slp(i,j) = sign(1.,slp(i,j)) * min(abs(slp(i,j)), 2. * min(dmx, dmn))
3136 ! * (G%mask2dT(i,j-1) * G%mask2dT(i,j) * G%mask2dT(i,j+1))
3137 endif
3138 enddo ; enddo
3139
3140 do j=jsl,jel ; do i=isl,iel
3141 ! Neighboring values should take into account any boundaries. The 3
3142 ! following sets of expressions are equivalent.
3143 h_jm1 = g%mask2dT(i,j-1) * h_in(i,j-1) + (1.0-g%mask2dT(i,j-1)) * h_in(i,j)
3144 h_jp1 = g%mask2dT(i,j+1) * h_in(i,j+1) + (1.0-g%mask2dT(i,j+1)) * h_in(i,j)
3145 ! Left/right values following Eq. B2 in Lin 1994, MWR (132)
3146 h_l(i,j) = 0.5*( h_jm1 + h_in(i,j) ) + onesixth*( slp(i,j-1) - slp(i,j) )
3147 h_r(i,j) = 0.5*( h_jp1 + h_in(i,j) ) + onesixth*( slp(i,j) - slp(i,j+1) )
3148 enddo ; enddo
3149 endif
3150
3151 select case(adv_limiter)
3152 case (limiter_adv_positive)
3153 call ppm_limit_pos(h_in, h_l, h_r, 0.0, g, isl, iel, jsl, jel)
3154 case (limiter_adv_minmod)
3155 call minmod_limiter(h_in, h_l, h_r, g, isl, iel, jsl, jel)
3156 end select
3157
3158end subroutine ppm_reconstruction_y
3159
3160!> Limits the left/right edge values of the PPM reconstruction
3161!! to give a reconstruction that is positive-definite. Here this is
3162!! reinterpreted as giving a constant value if the mean value is less
3163!! than h_min, with a minimum of h_min otherwise.
3164subroutine ppm_limit_pos(h_in, h_L, h_R, h_min, G, iis, iie, jis, jie)
3165 type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure.
3166 real, dimension(SZI_(G),SZJ_(G)), intent(in) :: h_in !< Energy density in each sector (2D)
3167 !! [H Z2 T-2 ~> m3 s-2 or J m-2]
3168 real, dimension(SZI_(G),SZJ_(G)), intent(inout) :: h_L !< Left edge value of reconstruction
3169 !! [H Z2 T-2 ~> m3 s-2 or J m-2]
3170 real, dimension(SZI_(G),SZJ_(G)), intent(inout) :: h_R !< Right edge value of reconstruction
3171 !! [H Z2 T-2 ~> m3 s-2 or J m-2]
3172 real, intent(in) :: h_min !< The minimum value that can be
3173 !! obtained by a concave parabolic fit
3174 !! [H Z2 T-2 ~> m3 s-2 or J m-2]
3175 integer, intent(in) :: iis !< Start i-index for computations
3176 integer, intent(in) :: iie !< End i-index for computations
3177 integer, intent(in) :: jis !< Start j-index for computations
3178 integer, intent(in) :: jie !< End j-index for computations
3179 ! Local variables
3180 real :: curv ! The cell-area normalized curvature [H Z2 T-2 ~> m3 s-2 or J m-2]
3181 real :: dh ! The difference between the edge values [H Z2 T-2 ~> m3 s-2 or J m-2]
3182 real :: scale ! A rescaling factor used to give a minimum cell value of at least h_min [nondim]
3183 integer :: i, j
3184
3185 do j=jis,jie ; do i=iis,iie
3186 ! This limiter prevents undershooting minima within the domain with
3187 ! values less than h_min.
3188 curv = 3.0*((h_l(i,j) + h_r(i,j)) - 2.0*h_in(i,j))
3189 if (curv > 0.0) then ! Only minima are limited.
3190 dh = h_r(i,j) - h_l(i,j)
3191 if (abs(dh) < curv) then ! The parabola's minimum is within the cell.
3192 if (h_in(i,j) <= h_min) then
3193 h_l(i,j) = h_in(i,j) ; h_r(i,j) = h_in(i,j)
3194 elseif (12.0*curv*(h_in(i,j) - h_min) < (curv**2 + 3.0*dh**2)) then
3195 ! The minimum value is h_in - (curv^2 + 3*dh^2)/(12*curv), and must
3196 ! be limited in this case. 0 < scale < 1.
3197 scale = 12.0*curv*(h_in(i,j) - h_min) / (curv**2 + 3.0*dh**2)
3198 h_l(i,j) = h_in(i,j) + scale*(h_l(i,j) - h_in(i,j))
3199 h_r(i,j) = h_in(i,j) + scale*(h_r(i,j) - h_in(i,j))
3200 endif
3201 endif
3202 endif
3203 enddo ; enddo
3204end subroutine ppm_limit_pos
3205
3206!> Limits the left/right edge values using the simple minmod limiter
3207!! written in a way that avoids branching in favor of intrinsics
3208subroutine minmod_limiter(h_in, h_L, h_R, G, iis, iie, jis, jie)
3209 type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure.
3210 real, dimension(SZI_(G),SZJ_(G)), intent(in) :: h_in !< Energy density in each sector (2D)
3211 !! [H Z2 T-2 ~> m3 s-2 or J m-2]
3212 real, dimension(SZI_(G),SZJ_(G)), intent(inout) :: h_L !< Left edge value of reconstruction
3213 !! [H Z2 T-2 ~> m3 s-2 or J m-2]
3214 real, dimension(SZI_(G),SZJ_(G)), intent(inout) :: h_R !< Right edge value of reconstruction
3215 !! [H Z2 T-2 ~> m3 s-2 or J m-2]
3216 integer, intent(in) :: iis !< Start i-index for computations
3217 integer, intent(in) :: iie !< End i-index for computations
3218 integer, intent(in) :: jis !< Start j-index for computations
3219 integer, intent(in) :: jie !< End j-index for computations
3220 ! Local variables
3221 real :: sign_h_L, sign_h_R, sign_h_in ! the signs of the edge and center values
3222 real :: sign_h_L_in, sign_h_R_in ! products of signs, detect crossing the zero line
3223 integer :: i, j
3224
3225 do j=jis,jie ; do i=iis,iie
3226
3227 sign_h_l = sign(1.0d0, h_l(i,j))
3228 sign_h_r = sign(1.0d0, h_r(i,j))
3229 sign_h_in = sign(1.0d0, h_in(i,j))
3230
3231 sign_h_l_in = sign_h_l * sign_h_in
3232 sign_h_r_in = sign_h_r * sign_h_in
3233
3234 ! if opposite signs, goes to zero else take the min of edge and centers values
3235 h_l(i,j) = (0.5 * (sign_h_l_in + 1.0)) * (sign_h_l * min(abs(h_l(i,j)), abs(h_in(i,j))))
3236 h_r(i,j) = (0.5 * (sign_h_r_in + 1.0)) * (sign_h_r * min(abs(h_r(i,j)), abs(h_in(i,j))))
3237
3238 enddo ; enddo
3239
3240end subroutine minmod_limiter
3241
3242subroutine register_int_tide_restarts(G, GV, US, param_file, CS, restart_CS)
3243 type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure
3244 type(verticalgrid_type),intent(in):: gv !< The ocean's vertical grid structure.
3245 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
3246 type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters
3247 type(int_tide_cs), pointer :: cs !< Internal tide control structure
3248 type(mom_restart_cs), pointer :: restart_cs !< MOM restart control structure
3249
3250 ! This subroutine is used to allocate and register any fields in this module
3251 ! that should be written to or read from the restart file.
3252 logical :: non_bous ! If true, this run is fully non-Boussinesq
3253 logical :: boussinesq ! If true, this run is fully Boussinesq
3254 logical :: semi_boussinesq ! If true, this run is partially non-Boussinesq
3255 integer :: num_freq, num_angle, num_mode
3256 integer :: isd, ied, jsd, jed, i, j, a, fr, m
3257 character(64) :: units
3258
3259 type(axis_info) :: axes_inttides(2)
3260 real, dimension(:), allocatable :: angles, freqs ! Lables for angles and frequencies [nondim]
3261 real :: hz2_t2_to_j_m2 ! unit conversion factor for Energy from internal to mks [H Z2 T-2 ~> m3 s-2 or J m-2]
3262
3263 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
3264
3265 hz2_t2_to_j_m2 = gv%H_to_MKS*(us%Z_to_m**2)*(us%s_to_T**2)
3266
3267 if (associated(cs)) then
3268 call mom_error(warning, "register_int_tide_restarts called "//&
3269 "with an associated control structure.")
3270 return
3271 endif
3272
3273 allocate(cs)
3274
3275 ! write extra axes
3276 call get_param(param_file, "MOM", "INTERNAL_TIDE_ANGLES", num_angle, default=24)
3277 call get_param(param_file, "MOM", "INTERNAL_TIDE_FREQS", num_freq, default=1)
3278 call get_param(param_file, "MOM", "INTERNAL_TIDE_MODES", num_mode, default=1)
3279
3280 ! define restart units depemding on Boussinesq
3281 call get_param(param_file, "MOM", "BOUSSINESQ", boussinesq, &
3282 "If true, make the Boussinesq approximation.", default=.true., do_not_log=.true.)
3283 call get_param(param_file, "MOM", "SEMI_BOUSSINESQ", semi_boussinesq, &
3284 "If true, do non-Boussinesq pressure force calculations and use mass-based "//&
3285 "thicknesses, but use RHO_0 to convert layer thicknesses into certain "//&
3286 "height changes. This only applies if BOUSSINESQ is false.", &
3287 default=.true., do_not_log=.true.)
3288 non_bous = .not.(boussinesq .or. semi_boussinesq)
3289
3290 units = "J m-2"
3291 if (boussinesq) units = "m3 s-2"
3292
3293 allocate (angles(num_angle))
3294 allocate (freqs(num_freq))
3295
3296 do a=1,num_angle ; angles(a) = a ; enddo
3297 do fr=1,num_freq ; freqs(fr) = fr ; enddo
3298
3299 call set_axis_info(axes_inttides(1), "angle", "", "angle direction", num_angle, angles, "N", 1)
3300 call set_axis_info(axes_inttides(2), "freq", "", "wave frequency", num_freq, freqs, "N", 1)
3301
3302 ! full energy array
3303 allocate(cs%En(isd:ied, jsd:jed, num_angle, num_freq, num_mode), source=0.0)
3304
3305 do m=1,num_mode ; do fr=1,num_freq
3306 call create_group_pass(cs%pass_En, cs%En(:,:,:,fr,m), g%Domain)
3307 enddo ; enddo
3308
3309 ! restart strategy: support for 5d restart is not yet available so we split into
3310 ! 4d restarts. Vertical modes >= 6 are dissipated locally and do not propagate
3311 ! so we only allow for 5 vertical modes and each has its own variable
3312
3313 ! allocate restart arrays
3314 allocate(cs%En_restart_mode1(isd:ied, jsd:jed, num_angle, num_freq), source=0.0)
3315 if (num_mode >= 2) allocate(cs%En_restart_mode2(isd:ied, jsd:jed, num_angle, num_freq), source=0.0)
3316 if (num_mode >= 3) allocate(cs%En_restart_mode3(isd:ied, jsd:jed, num_angle, num_freq), source=0.0)
3317 if (num_mode >= 4) allocate(cs%En_restart_mode4(isd:ied, jsd:jed, num_angle, num_freq), source=0.0)
3318 if (num_mode >= 5) allocate(cs%En_restart_mode5(isd:ied, jsd:jed, num_angle, num_freq), source=0.0)
3319
3320 ! register all 4d restarts and copy into full Energy array when restarting from previous state
3321 call register_restart_field(cs%En_restart_mode1(:,:,:,:), "IW_energy_mode1", .false., restart_cs, &
3322 longname="The internal wave energy density f(i,j,angle,freq) for mode 1", &
3323 units=units, conversion=hz2_t2_to_j_m2, z_grid='1', t_grid="s", &
3324 extra_axes=axes_inttides)
3325
3326 do fr=1,num_freq ; do a=1,num_angle ; do j=jsd,jed ; do i=isd,ied
3327 cs%En(i,j,a,fr,1) = cs%En_restart_mode1(i,j,a,fr)
3328 enddo ; enddo ; enddo ; enddo
3329
3330 if (num_mode >= 2) then
3331 call register_restart_field(cs%En_restart_mode2(:,:,:,:), "IW_energy_mode2", .false., restart_cs, &
3332 longname="The internal wave energy density f(i,j,angle,freq) for mode 2", &
3333 units=units, conversion=hz2_t2_to_j_m2, z_grid='1', t_grid="s", &
3334 extra_axes=axes_inttides)
3335
3336 do fr=1,num_freq ; do a=1,num_angle ; do j=jsd,jed ; do i=isd,ied
3337 cs%En(i,j,a,fr,2) = cs%En_restart_mode2(i,j,a,fr)
3338 enddo ; enddo ; enddo ; enddo
3339
3340 endif
3341
3342 if (num_mode >= 3) then
3343 call register_restart_field(cs%En_restart_mode3(:,:,:,:), "IW_energy_mode3", .false., restart_cs, &
3344 longname="The internal wave energy density f(i,j,angle,freq) for mode 3", &
3345 units=units, conversion=hz2_t2_to_j_m2, z_grid='1', t_grid="s", &
3346 extra_axes=axes_inttides)
3347
3348 do fr=1,num_freq ; do a=1,num_angle ; do j=jsd,jed ; do i=isd,ied
3349 cs%En(i,j,a,fr,3) = cs%En_restart_mode3(i,j,a,fr)
3350 enddo ; enddo ; enddo ; enddo
3351
3352 endif
3353
3354 if (num_mode >= 4) then
3355 call register_restart_field(cs%En_restart_mode4(:,:,:,:), "IW_energy_mode4", .false., restart_cs, &
3356 longname="The internal wave energy density f(i,j,angle,freq) for mode 4", &
3357 units=units, conversion=hz2_t2_to_j_m2, z_grid='1', t_grid="s", &
3358 extra_axes=axes_inttides)
3359
3360 do fr=1,num_freq ; do a=1,num_angle ; do j=jsd,jed ; do i=isd,ied
3361 cs%En(i,j,a,fr,4) = cs%En_restart_mode4(i,j,a,fr)
3362 enddo ; enddo ; enddo ; enddo
3363
3364 endif
3365
3366 if (num_mode >= 5) then
3367 call register_restart_field(cs%En_restart_mode5(:,:,:,:), "IW_energy_mode5", .false., restart_cs, &
3368 longname="The internal wave energy density f(i,j,angle,freq) for mode 5", &
3369 units=units, conversion=hz2_t2_to_j_m2, z_grid='1', t_grid="s", &
3370 extra_axes=axes_inttides)
3371
3372 do fr=1,num_freq ; do a=1,num_angle ; do j=jsd,jed ; do i=isd,ied
3373 cs%En(i,j,a,fr,5) = cs%En_restart_mode5(i,j,a,fr)
3374 enddo ; enddo ; enddo ; enddo
3375
3376 endif
3377
3378end subroutine register_int_tide_restarts
3379
3380!> This subroutine initializes the internal tides module.
3381subroutine internal_tides_init(Time, G, GV, US, param_file, diag, CS)
3382 type(time_type), target, intent(in) :: time !< The current model time.
3383 type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure.
3384 type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure.
3385 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
3386 type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time
3387 !! parameters.
3388 type(diag_ctrl), target, intent(in) :: diag !< A structure that is used to regulate
3389 !! diagnostic output.
3390 type(int_tide_cs), pointer :: cs !< Internal tide control structure
3391
3392 ! Local variables
3393 real :: angle_size ! size of wedges [rad]
3394 real, allocatable :: angles(:) ! orientations of wedge centers [rad]
3395 real, dimension(:,:), allocatable :: h2 ! topographic roughness scale squared [Z2 ~> m2]
3396 real :: kappa_itides ! characteristic topographic wave number [L-1 ~> m-1]
3397 real, dimension(:,:), allocatable :: ridge_temp ! array for temporary storage of flags
3398 ! of cells with double-reflecting ridges [nondim]
3399 real, dimension(:,:), allocatable :: tmp_decay ! a temp array to store decay rates [T-1 ~> s-1]
3400 real :: decay_rate ! A constant rate at which internal tide energy is
3401 ! lost to the interior ocean internal wave field [T-1 ~> s-1].
3402 logical :: use_int_tides, use_temperature
3403 logical :: om4_remap_via_sub_cells ! Use the OM4-era ramap_via_sub_cells for calculating the EBT structure
3404 real :: igw_c1_thresh ! A threshold first mode internal wave speed below which all higher
3405 ! mode speeds are not calculated but simply assigned a speed of 0 [L T-1 ~> m s-1].
3406 real :: kappa_h2_factor ! A roughness scaling factor [nondim]
3407 real :: rms_roughness_frac ! The maximum RMS topographic roughness as a fraction of the
3408 ! nominal ocean depth, or a negative value for no limit [nondim]
3409 real :: period ! A tidal period read from namelist [T ~> s]
3410 real :: hz2_t2_to_j_m2 ! unit conversion factor for Energy from internal units
3411 ! to mks [T2 kg H-1 Z-2 s-2 ~> kg m-3 or 1]
3412 real :: hz2_t3_to_w_m2 ! unit conversion factor for TKE from internal units
3413 ! to mks [T3 kg H-1 Z-2 s-3 ~> kg m-3 or 1]
3414 real :: j_m2_to_hz2_t2 ! unit conversion factor for Energy from mks to internal
3415 ! units [H Z2 s2 T-2 kg-1 ~> m3 kg-1 or 1]
3416 integer :: num_angle, num_freq, num_mode, m, fr
3417 integer :: isd, ied, jsd, jed, a, id_ang, i, j, nz
3418 type(axes_grp) :: axes_ang
3419 ! This include declares and sets the variable "version".
3420# include "version_variable.h"
3421 character(len=40) :: mdl = "MOM_internal_tides" ! This module's name.
3422 character(len=16), dimension(8) :: freq_name
3423 character(len=40) :: var_name
3424 character(len=160) :: var_descript
3425 character(len=200) :: filename
3426 character(len=200) :: refl_angle_file
3427 character(len=200) :: refl_pref_file, refl_dbl_file, trans_file
3428 character(len=200) :: h2_file, decay_file
3429 character(len=80) :: rough_var ! Input file variable names
3430 character(len=80) :: tmpstr
3431
3432 character(len=240), dimension(:), allocatable :: energy_fractions
3433 character(len=240) :: periods
3434
3435 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
3436 nz = gv%ke
3437
3438 hz2_t2_to_j_m2 = gv%H_to_kg_m2*(us%Z_to_m**2)*(us%s_to_T**2)
3439 hz2_t3_to_w_m2 = gv%H_to_kg_m2*(us%Z_to_m**2)*(us%s_to_T**3)
3440 j_m2_to_hz2_t2 = gv%kg_m2_to_H*(us%m_to_Z**2)*(us%T_to_s**2)
3441
3442 cs%initialized = .true.
3443
3444 use_int_tides = .false.
3445 call read_param(param_file, "INTERNAL_TIDES", use_int_tides)
3446 cs%do_int_tides = use_int_tides
3447 if (.not.use_int_tides) return
3448
3449 use_temperature = .true.
3450 call read_param(param_file, "ENABLE_THERMODYNAMICS", use_temperature)
3451 if (.not.use_temperature) call mom_error(fatal, &
3452 "internal_tides_init: internal_tides only works with ENABLE_THERMODYNAMICS defined.")
3453
3454 ! Set number of frequencies, angles, and modes to consider
3455 num_freq = 1 ; num_angle = 24 ; num_mode = 1
3456 call read_param(param_file, "INTERNAL_TIDE_FREQS", num_freq)
3457 call read_param(param_file, "INTERNAL_TIDE_ANGLES", num_angle)
3458 call read_param(param_file, "INTERNAL_TIDE_MODES", num_mode)
3459 if (.not.((num_freq > 0) .and. (num_angle > 0) .and. (num_mode > 0))) return
3460 cs%nFreq = num_freq ; cs%nAngle = num_angle ; cs%nMode = num_mode
3461
3462 allocate(energy_fractions(num_freq))
3463 allocate(cs%fraction_tidal_input(num_freq,num_mode))
3464
3465 call read_param(param_file, "ENERGY_FRACTION_PER_MODE", energy_fractions)
3466
3467 do fr=1,num_freq ; do m=1,num_mode
3468 cs%fraction_tidal_input(fr,m) = extract_real(energy_fractions(fr), " ,", m, 0.)
3469 enddo ; enddo
3470
3471 ! Allocate phase speed array
3472 allocate(cs%cp(isd:ied, jsd:jed, num_freq, num_mode), source=0.0)
3473
3474 ! Allocate and populate frequency array (each a multiple of first for now)
3475 allocate(cs%frequency(num_freq))
3476
3477
3478 call get_param(param_file, mdl, "DEBUG", cs%debug, &
3479 "If true, write out verbose debugging data.", &
3480 default=.false., debuggingparam=.true.)
3481
3482 ! The periods of the tidal constituents for internal tides raytracing
3483 call read_param(param_file, "TIDAL_PERIODS", periods)
3484
3485 do fr=1,num_freq
3486 period = us%s_to_T*extract_real(periods, " ,", fr, 0.)
3487 if (period == 0.) call mom_error(fatal, "MOM_internal_tides: invalid tidal period")
3488 cs%frequency(fr) = 8.0*atan(1.0)/period
3489 enddo
3490
3491 ! Read all relevant parameters and write them to the model log.
3492
3493 cs%Time => time ! direct a pointer to the current model time target
3494
3495 call get_param(param_file, mdl, "INPUTDIR", cs%inputdir, default=".")
3496 cs%inputdir = slasher(cs%inputdir)
3497
3498 call log_version(param_file, mdl, version, "")
3499
3500 call get_param(param_file, mdl, "INTERNAL_TIDE_FREQS", num_freq, &
3501 "The number of distinct internal tide frequency bands "//&
3502 "that will be calculated.", default=1)
3503 call get_param(param_file, mdl, "INTERNAL_TIDE_MODES", num_mode, &
3504 "The number of distinct internal tide modes "//&
3505 "that will be calculated.", default=1)
3506 call get_param(param_file, mdl, "INTERNAL_TIDE_ANGLES", num_angle, &
3507 "The number of angular resolution bands for the internal "//&
3508 "tide calculations.", default=24)
3509 call get_param(param_file, mdl, "DT_ITIDES", cs%dt_itides, &
3510 "The timestep for internal tides ray-tracing scheme. "//&
3511 "If set to -1 (default), it uses the same value as DT_THERM.", &
3512 units="s", default=-1., scale=us%s_to_T)
3513
3514 if (use_int_tides) then
3515 if ((num_freq <= 0) .and. (num_mode <= 0) .and. (num_angle <= 0)) then
3516 call mom_error(warning, "Internal tides were enabled, but the number "//&
3517 "of requested frequencies, modes and angles were not all positive.")
3518 return
3519 endif
3520 else
3521 if ((num_freq > 0) .and. (num_mode > 0) .and. (num_angle > 0)) then
3522 call mom_error(warning, "Internal tides were not enabled, even though "//&
3523 "the number of requested frequencies, modes and angles were all "//&
3524 "positive.")
3525 return
3526 endif
3527 endif
3528
3529 if (cs%NFreq /= num_freq) call mom_error(fatal, "Internal_tides_init: "//&
3530 "Inconsistent number of frequencies.")
3531 if (cs%NAngle /= num_angle) call mom_error(fatal, "Internal_tides_init: "//&
3532 "Inconsistent number of angles.")
3533 if (cs%nMode /= num_mode) call mom_error(fatal, "Internal_tides_init: "//&
3534 "Inconsistent number of modes.")
3535 if (4*(num_angle/4) /= num_angle) call mom_error(fatal, &
3536 "Internal_tides_init: INTERNAL_TIDE_ANGLES must be a multiple of 4.")
3537
3538 cs%diag => diag
3539
3540 call get_param(param_file, mdl, "INTERNAL_TIDES_UPDATE_KD", cs%update_Kd, &
3541 "If true, internal tides ray tracing changes Kd for dynamics.", &
3542 default=.false.)
3543 call get_param(param_file, mdl, "INTERNAL_TIDES_REFRACTION", cs%apply_refraction, &
3544 "If true, internal tides ray tracing does refraction.", &
3545 default=.true.)
3546 call get_param(param_file, mdl, "INTERNAL_TIDES_PROPAGATION", cs%apply_propagation, &
3547 "If true, internal tides ray tracing does propagate.", &
3548 default=.true.)
3549 call get_param(param_file, mdl, "INTERNAL_TIDES_ONLY_INIT_FORCING", cs%init_forcing_only, &
3550 "If true, internal tides ray tracing only applies forcing at first step (debugging).", &
3551 default=.false.)
3552 call get_param(param_file, mdl, "TURN_CRITICAL_LAT", cs%turn_critical_lat, &
3553 "If true, internal tides rays turn at the critical latitude.", &
3554 default=.true.)
3555 call get_param(param_file, mdl, "REFLECT_CRITICAL_LAT", cs%reflect_critical_lat, &
3556 "If true, internal tides rays reflect at the critical latitude. "//&
3557 "If false, rays turn parallel to the critical latitude", &
3558 default=.true.)
3559 call get_param(param_file, mdl, "INTERNAL_TIDES_FORCE_POS_EN", cs%force_posit_En, &
3560 "If true, force energy to be positive by removing subroundoff negative values.", &
3561 default=.true.)
3562 call get_param(param_file, mdl, "KD_MIN", cs%Kd_min, &
3563 "The minimum diapycnal diffusivity.", &
3564 units="m2 s-1", default=2e-6, scale=gv%m2_s_to_HZ_T)
3565 call get_param(param_file, mdl, "MINTHICK_TKE_TO_KD", cs%min_thick_layer_Kd, &
3566 "The minimum thickness allowed with TKE_to_Kd.", &
3567 units="m", default=1e-6, scale=gv%m_to_H)
3568 call get_param(param_file, mdl, "ITIDES_MIXING_EFFIC", cs%mixing_effic, &
3569 "Mixing efficiency for internal tides raytracing", &
3570 units="nondim", default=0.2)
3571 call get_param(param_file, mdl, "MAX_TKE_TO_KD", cs%max_TKE_to_Kd, &
3572 "Limiter for TKE_to_Kd.", &
3573 units="s2 m-1", default=1e9, scale=us%Z_to_m*us%s_to_T**2)
3574 call get_param(param_file, mdl, "INTERNAL_TIDE_DECAY_RATE", decay_rate, &
3575 "The rate at which internal tide energy is lost to the "//&
3576 "interior ocean internal wave field.", &
3577 units="s-1", default=0.0, scale=us%T_to_s)
3578 call get_param(param_file, mdl, "USE_2D_INTERNAL_TIDE_DECAY_RATE", cs%use_2d_decay_rate, &
3579 "If true, use a spatially varying decay rate for leakage loss in the "// &
3580 "internal tide code.", default=.false.)
3581 call get_param(param_file, mdl, "INTERNAL_TIDE_VOLUME_BASED_CFL", cs%vol_CFL, &
3582 "If true, use the ratio of the open face lengths to the "//&
3583 "tracer cell areas when estimating CFL numbers in the "//&
3584 "internal tide code.", default=.false.)
3585 call get_param(param_file, mdl, "INTERNAL_TIDE_SIMPLE_2ND_PPM", cs%simple_2nd, &
3586 "If true, CONTINUITY_PPM uses a simple 2nd order "//&
3587 "(arithmetic mean) interpolation of the edge values. "//&
3588 "This may give better PV conservation properties. While "//&
3589 "it formally reduces the accuracy of the continuity "//&
3590 "solver itself in the strongly advective limit, it does "//&
3591 "not reduce the overall order of accuracy of the dynamic "//&
3592 "core.", default=.false.)
3593 call get_param(param_file, mdl, "INTERNAL_TIDE_UPWIND_1ST", cs%upwind_1st, &
3594 "If true, the internal tide ray-tracing advection uses "//&
3595 "1st-order upwind advection. This scheme is highly "//&
3596 "continuity solver. This scheme is highly "//&
3597 "diffusive but may be useful for debugging.", default=.false.)
3598 call get_param(param_file, mdl, "INTERNAL_TIDE_ADV_LIMITER", tmpstr, &
3599 "Choose the limiter scheme used for the internal tide advection scheme, "//&
3600 "available schemes are: \n"//&
3601 "\t POSITIVE - a positive definite scheme similar to the continuity solver. \n"//&
3602 "\t MINMOD - the simplest limiter.", default=limiter_adv_minmod_string)
3603
3604 tmpstr = uppercase(tmpstr)
3605 select case (tmpstr)
3606 case (limiter_adv_positive_string)
3607 cs%itides_adv_limiter = limiter_adv_positive
3608 case (limiter_adv_minmod_string)
3609 cs%itides_adv_limiter = limiter_adv_minmod
3610 case default
3611 call mom_mesg('internal_tide_init: Advection limiter ="'//trim(tmpstr)//'"', 0)
3612 call mom_error(fatal, "internal_tide_init: Unrecognized setting "// &
3613 "#define INTERNAL_TIDE_ADV_LIMITER "//trim(tmpstr)//" found in input file.")
3614 end select
3615
3616 call get_param(param_file, mdl, "INTERNAL_TIDE_BACKGROUND_DRAG", cs%apply_background_drag, &
3617 "If true, the internal tide ray-tracing advection uses a background drag "//&
3618 "term as a sink.", default=.false.)
3619 call get_param(param_file, mdl, "INTERNAL_TIDE_QUAD_DRAG", cs%apply_bottom_drag, &
3620 "If true, the internal tide ray-tracing advection uses "//&
3621 "a quadratic bottom drag term as a sink.", default=.false.)
3622 call get_param(param_file, mdl, "INTERNAL_TIDE_WAVE_DRAG", cs%apply_wave_drag, &
3623 "If true, apply scattering due to small-scale roughness as a sink.", &
3624 default=.false.)
3625 call get_param(param_file, mdl, "INTERNAL_TIDE_RESIDUAL_DRAG", cs%apply_residual_drag, &
3626 "If true, apply drag due to critical slopes", &
3627 default=.false.)
3628 call get_param(param_file, mdl, "INTERNAL_TIDE_DRAG_MIN_DEPTH", cs%drag_min_depth, &
3629 "The minimum total ocean thickness that will be used in the denominator "//&
3630 "of the quadratic drag terms for internal tides.", &
3631 units="m", default=1.0, scale=gv%m_to_H, do_not_log=.not.cs%apply_bottom_drag)
3632 cs%drag_min_depth = max(cs%drag_min_depth, gv%H_subroundoff)
3633 call get_param(param_file, mdl, "INTERNAL_TIDE_FROUDE_DRAG", cs%apply_Froude_drag, &
3634 "If true, apply wave breaking as a sink.", &
3635 default=.false.)
3636 call get_param(param_file, mdl, "EN_CHECK_TOLERANCE", cs%En_check_tol, &
3637 "An energy density tolerance for flagging points with small negative "//&
3638 "internal tide energy.", &
3639 units="J m-2", default=1.0, scale=j_m2_to_hz2_t2, &
3640 do_not_log=.not.cs%apply_Froude_drag)
3641 call get_param(param_file, mdl, "EN_UNDERFLOW", cs%En_underflow, &
3642 "A small energy density below which Energy is set to zero.", &
3643 units="J m-2", default=1.0e-100, scale=j_m2_to_hz2_t2)
3644 call get_param(param_file, mdl, "EN_RESTART_POWER", cs%En_restart_power, &
3645 "A power factor to save larger values x 2**(power) in restart files.", &
3646 units="nondim", default=0)
3647 call get_param(param_file, mdl, "CDRAG", cs%cdrag, &
3648 "CDRAG is the drag coefficient relating the magnitude of "//&
3649 "the velocity field to the bottom stress.", &
3650 units="nondim", default=0.003)
3651 call get_param(param_file, mdl, "INTERNAL_WAVE_CG1_THRESH", igw_c1_thresh, &
3652 "A minimal value of the first mode internal wave speed below which all higher "//&
3653 "mode speeds are not calculated but are simply reported as 0. This must be "//&
3654 "non-negative for the wave_speeds routine to be used.", &
3655 units="m s-1", default=0.01, scale=us%m_s_to_L_T)
3656 call get_param(param_file, mdl, "REMAPPING_USE_OM4_SUBCELLS", om4_remap_via_sub_cells, &
3657 do_not_log=.true., default=.true.)
3658 call get_param(param_file, mdl, "INTWAVE_REMAPPING_USE_OM4_SUBCELLS", om4_remap_via_sub_cells, &
3659 "If true, use the OM4 remapping-via-subcells algorithm for calculating EBT structure. "//&
3660 "See REMAPPING_USE_OM4_SUBCELLS for details. "//&
3661 "We recommend setting this option to false.", default=om4_remap_via_sub_cells)
3662 call get_param(param_file, mdl, "UNIFORM_TEST_CG", cs%uniform_test_cg, &
3663 "If positive, a uniform group velocity of internal tide for test case", &
3664 default=-1., units="m s-1", scale=us%m_s_to_L_T)
3665 call get_param(param_file, mdl, "INTERNAL_TIDE_ENERGIZED_ANGLE", cs%energized_angle, &
3666 "If positive, only one angular band of the internal tides "//&
3667 "gets all of the energy. (This is for debugging.)", default=-1)
3668 call get_param(param_file, mdl, "USE_PPM_ANGULAR", cs%use_PPMang, &
3669 "If true, use PPM for advection of energy in angular space.", &
3670 default=.false.)
3671 call get_param(param_file, mdl, "GAMMA_ITIDES", cs%q_itides, &
3672 "The fraction of the internal tidal energy that is "//&
3673 "dissipated locally with INT_TIDE_DISSIPATION. "//&
3674 "THIS NAME COULD BE BETTER.", &
3675 units="nondim", default=0.3333)
3676 call get_param(param_file, mdl, "KAPPA_ITIDES", kappa_itides, &
3677 "A topographic wavenumber used with INT_TIDE_DISSIPATION. "//&
3678 "The default is 2pi/10 km, as in St.Laurent et al. 2002.", &
3679 units="m-1", default=8.e-4*atan(1.0), scale=us%L_to_m)
3680 call get_param(param_file, mdl, "KAPPA_H2_FACTOR", kappa_h2_factor, &
3681 "A scaling factor for the roughness amplitude with "//&
3682 "INT_TIDE_DISSIPATION.", units="nondim", default=1.0)
3683 call get_param(param_file, mdl, "GAMMA_OSBORN", cs%gamma_osborn, &
3684 "The mixing efficiency for internan tides from Osborn 1980 ", &
3685 units="nondim", default=0.2)
3686 call get_param(param_file, mdl, "INT_TIDE_DECAY_SCALE", cs%Int_tide_decay_scale, &
3687 "The decay scale away from the bottom for tidal TKE with "//&
3688 "the new coding when INT_TIDE_DISSIPATION is used.", &
3689 units="m", default=500.0, scale=gv%m_to_H)
3690 call get_param(param_file, mdl, "INT_TIDE_DECAY_SCALE_SLOPES", cs%Int_tide_decay_scale_slope, &
3691 "The slope decay scale away from the bottom for tidal TKE with "//&
3692 "the new coding when INT_TIDE_DISSIPATION is used.", &
3693 units="m", default=100.0, scale=gv%m_to_H)
3694
3695 ! Allocate various arrays needed for loss rates
3696 allocate(h2(isd:ied,jsd:jed), source=0.0)
3697 allocate(cs%TKE_itidal_loss_fixed(isd:ied,jsd:jed), source=0.0)
3698 allocate(cs%TKE_leak_loss(isd:ied,jsd:jed,num_angle,num_freq,num_mode), source=0.0)
3699 allocate(cs%TKE_quad_loss(isd:ied,jsd:jed,num_angle,num_freq,num_mode), source=0.0)
3700 allocate(cs%TKE_itidal_loss(isd:ied,jsd:jed,num_angle,num_freq,num_mode), source=0.0)
3701 allocate(cs%TKE_Froude_loss(isd:ied,jsd:jed,num_angle,num_freq,num_mode), source=0.0)
3702 allocate(cs%TKE_residual_loss(isd:ied,jsd:jed,num_angle,num_freq,num_mode), source=0.0)
3703 allocate(cs%TKE_slope_loss(isd:ied,jsd:jed,num_angle,num_freq,num_mode), source=0.0)
3704 allocate(cs%tot_leak_loss(isd:ied,jsd:jed), source=0.0)
3705 allocate(cs%tot_quad_loss(isd:ied,jsd:jed), source=0.0)
3706 allocate(cs%tot_itidal_loss(isd:ied,jsd:jed), source=0.0)
3707 allocate(cs%tot_Froude_loss(isd:ied,jsd:jed), source=0.0)
3708 allocate(cs%tot_residual_loss(isd:ied,jsd:jed), source=0.0)
3709 allocate(cs%u_struct_bot(isd:ied,jsd:jed,num_mode), source=0.0)
3710 allocate(cs%u_struct_max(isd:ied,jsd:jed,num_mode), source=0.0)
3711 allocate(cs%int_w2(isd:ied,jsd:jed,num_mode), source=0.0)
3712 allocate(cs%int_U2(isd:ied,jsd:jed,num_mode), source=0.0)
3713 allocate(cs%int_N2w2(isd:ied,jsd:jed,num_mode), source=0.0)
3714 allocate(cs%w_struct(isd:ied,jsd:jed,1:nz+1,num_mode), source=0.0)
3715 allocate(cs%u_struct(isd:ied,jsd:jed,1:nz,num_mode), source=0.0)
3716 allocate(cs%error_mode(num_freq,num_mode), source=0.0)
3717 allocate(cs%En_ini_glo(num_freq,num_mode), source=0.0)
3718 allocate(cs%En_end_glo(num_freq,num_mode), source=0.0)
3719 allocate(cs%TKE_leak_loss_glo_dt(num_freq,num_mode), source=0.0)
3720 allocate(cs%TKE_quad_loss_glo_dt(num_freq,num_mode), source=0.0)
3721 allocate(cs%TKE_Froude_loss_glo_dt(num_freq,num_mode), source=0.0)
3722 allocate(cs%TKE_itidal_loss_glo_dt(num_freq,num_mode), source=0.0)
3723 allocate(cs%TKE_residual_loss_glo_dt(num_freq,num_mode), source=0.0)
3724 allocate(cs%TKE_input_glo_dt(num_freq,num_mode), source=0.0)
3725 allocate(cs%decay_rate_2d(isd:ied,jsd:jed,num_freq,num_mode), source=0.0)
3726 allocate(tmp_decay(isd:ied,jsd:jed), source=0.0)
3727
3728 if (cs%use_2d_decay_rate) then
3729 call get_param(param_file, mdl, "ITIDES_DECAY_FILE", decay_file, &
3730 "The path to the file containing the decay rates "//&
3731 "for internal tides with USE_2D_INTERNAL_TIDE_DECAY_RATE.", &
3732 fail_if_missing=.true.)
3733 do m=1,num_mode ; do fr=1,num_freq
3734 ! read 2d field for each harmonic
3735 filename = trim(cs%inputdir) // trim(decay_file)
3736 write(var_name, '("decay_rate_freq",i1,"_mode",i1)') fr, m
3737 call mom_read_data(filename, var_name, tmp_decay, g%domain, scale=us%T_to_s)
3738 do j=g%jsc,g%jec ; do i=g%isc,g%iec
3739 cs%decay_rate_2d(i,j,fr,m) = tmp_decay(i,j)
3740 enddo ; enddo
3741 enddo ; enddo
3742 else
3743 do m=1,num_mode ; do fr=1,num_freq ; do j=g%jsc,g%jec ; do i=g%isc,g%iec
3744 cs%decay_rate_2d(i,j,fr,m) = decay_rate
3745 enddo ; enddo ; enddo ; enddo
3746 endif
3747
3748 do m=1,num_mode
3749 call pass_var(cs%decay_rate_2d(:,:,:,m), g%domain)
3750 enddo
3751
3752 ! Compute the fixed part of the bottom drag loss from baroclinic modes
3753 call get_param(param_file, mdl, "H2_FILE", h2_file, &
3754 "The path to the file containing the sub-grid-scale "//&
3755 "topographic roughness amplitude with INT_TIDE_DISSIPATION.", &
3756 fail_if_missing=.true.)
3757 filename = trim(cs%inputdir) // trim(h2_file)
3758 call log_param(param_file, mdl, "INPUTDIR/H2_FILE", filename)
3759 call get_param(param_file, mdl, "ROUGHNESS_VARNAME", rough_var, &
3760 "The name in the input file of the squared sub-grid-scale "//&
3761 "topographic roughness amplitude variable.", default="h2")
3762 call get_param(param_file, mdl, "INTERNAL_TIDE_ROUGHNESS_FRAC", rms_roughness_frac, &
3763 "The maximum RMS topographic roughness as a fraction of the nominal ocean depth, "//&
3764 "or a negative value for no limit.", units="nondim", default=0.1)
3765
3766 call mom_read_data(filename, rough_var, h2, g%domain, scale=us%m_to_Z**2)
3767 do j=g%jsc,g%jec ; do i=g%isc,g%iec
3768 ! Restrict RMS topographic roughness to a fraction (10 percent by default) of the column depth.
3769 if (rms_roughness_frac >= 0.0) then
3770 h2(i,j) = max(min((rms_roughness_frac * max(g%meanSL(i,j) + g%bathyT(i,j), 0.0))**2, h2(i,j)), 0.0)
3771 else
3772 h2(i,j) = max(h2(i,j), 0.0)
3773 endif
3774 ! Compute the fixed part; units are [R Z4 H-1 L-2 ~> kg m-2 or m] here
3775 ! will be multiplied by N and the squared near-bottom velocity (and by the
3776 ! near-bottom density in non-Boussinesq mode) to get into [H Z2 T-3 ~> m3 s-3 or W m-2]
3777 cs%TKE_itidal_loss_fixed(i,j) = 0.5*kappa_h2_factor* gv%H_to_RZ * us%L_to_Z*kappa_itides * h2(i,j)
3778 enddo ; enddo
3779
3780 deallocate(h2)
3781
3782 ! Read in prescribed coast/ridge/shelf angles from file
3783 call get_param(param_file, mdl, "REFL_ANGLE_FILE", refl_angle_file, &
3784 "The path to the file containing the local angle of "//&
3785 "the coastline/ridge/shelf with respect to the equator.", &
3786 fail_if_missing=.false., default='')
3787 filename = trim(cs%inputdir) // trim(refl_angle_file)
3788 allocate(cs%refl_angle(isd:ied,jsd:jed), source=cs%nullangle)
3789 if (file_exists(filename, g%domain)) then
3790 call log_param(param_file, mdl, "INPUTDIR/REFL_ANGLE_FILE", filename)
3791 call mom_read_data(filename, 'refl_angle', cs%refl_angle, g%domain)
3792 else
3793 if (trim(refl_angle_file) /= '' ) call mom_error(fatal, &
3794 "REFL_ANGLE_FILE: "//trim(filename)//" not found")
3795 endif
3796 ! replace NaNs with null value
3797 do j=g%jsc,g%jec ; do i=g%isc,g%iec
3798 if (is_nan(cs%refl_angle(i,j))) cs%refl_angle(i,j) = cs%nullangle
3799 enddo ; enddo
3800 call pass_var(cs%refl_angle, g%domain)
3801
3802 ! Read in prescribed partial reflection coefficients from file
3803 call get_param(param_file, mdl, "REFL_PREF_FILE", refl_pref_file, &
3804 "The path to the file containing the reflection coefficients.", &
3805 fail_if_missing=.false., default='')
3806 filename = trim(cs%inputdir) // trim(refl_pref_file)
3807 allocate(cs%refl_pref(isd:ied,jsd:jed), source=1.0)
3808 if (file_exists(filename, g%domain)) then
3809 call log_param(param_file, mdl, "INPUTDIR/REFL_PREF_FILE", filename)
3810 call mom_read_data(filename, 'refl_pref', cs%refl_pref, g%domain)
3811 else
3812 if (trim(refl_pref_file) /= '' ) call mom_error(fatal, &
3813 "REFL_PREF_FILE: "//trim(filename)//" not found")
3814 endif
3815 !CS%refl_pref = CS%refl_pref*1 ! adjust partial reflection if desired
3816 call pass_var(cs%refl_pref, g%domain)
3817
3818 ! Tag reflection cells with partial reflection (done here for speed)
3819 allocate(cs%refl_pref_logical(isd:ied,jsd:jed), source=.false.)
3820 do j=jsd,jed ; do i=isd,ied
3821 ! flag cells with partial reflection
3822 if ((cs%refl_angle(i,j) /= cs%nullangle) .and. &
3823 (cs%refl_pref(i,j) < 1.0) .and. (cs%refl_pref(i,j) > 0.0)) then
3824 cs%refl_pref_logical(i,j) = .true.
3825 endif
3826 enddo ; enddo
3827
3828 ! Read in double-reflective (ridge) tags from file
3829 call get_param(param_file, mdl, "REFL_DBL_FILE", refl_dbl_file, &
3830 "The path to the file containing the double-reflective ridge tags.", &
3831 fail_if_missing=.false., default='')
3832 filename = trim(cs%inputdir) // trim(refl_dbl_file)
3833 allocate(ridge_temp(isd:ied,jsd:jed), source=0.0)
3834 if (file_exists(filename, g%domain)) then
3835 call log_param(param_file, mdl, "INPUTDIR/REFL_DBL_FILE", filename)
3836 call mom_read_data(filename, 'refl_dbl', ridge_temp, g%domain)
3837 else
3838 if (trim(refl_dbl_file) /= '' ) call mom_error(fatal, &
3839 "REFL_DBL_FILE: "//trim(filename)//" not found")
3840 endif
3841 call pass_var(ridge_temp, g%domain)
3842 allocate(cs%refl_dbl(isd:ied,jsd:jed), source=.false.)
3843 do j=jsd,jed ; do i=isd,ied
3844 cs%refl_dbl(i,j) = (ridge_temp(i,j) == 1)
3845 enddo ; enddo
3846
3847 ! Read in the transmission coefficient and infer the residual
3848 call get_param(param_file, mdl, "TRANS_FILE", trans_file, &
3849 "The path to the file containing the transmission coefficent for internal tides.", &
3850 fail_if_missing=.false., default='')
3851 filename = trim(cs%inputdir) // trim(trans_file)
3852 allocate(cs%trans(isd:ied,jsd:jed), source=0.0)
3853 if (file_exists(filename, g%domain)) then
3854 call log_param(param_file, mdl, "INPUTDIR/TRANS_FILE", filename)
3855 call mom_read_data(filename, 'trans', cs%trans, g%domain)
3856 else
3857 if (trim(trans_file) /= '' ) call mom_error(fatal, &
3858 "TRANS_FILE: "//trim(filename)//" not found")
3859 endif
3860
3861 call pass_var(cs%trans, g%domain)
3862
3863 ! residual
3864 allocate(cs%residual(isd:ied,jsd:jed), source=0.0)
3865 if (cs%apply_residual_drag) then
3866 do j=g%jsc,g%jec ; do i=g%isc,g%iec
3867 if (cs%refl_pref_logical(i,j)) then
3868 cs%residual(i,j) = 1. - (cs%refl_pref(i,j) - cs%trans(i,j))
3869 endif
3870 enddo ; enddo
3871 call pass_var(cs%residual, g%domain)
3872 else
3873 ! report residual of transmission/reflection onto reflection
3874 ! this ensure energy budget is conserved
3875 do j=g%jsc,g%jec ; do i=g%isc,g%iec
3876 if (cs%refl_pref_logical(i,j)) then
3877 cs%refl_pref(i,j) = 1. - cs%trans(i,j)
3878 endif
3879 enddo ; enddo
3880 call pass_var(cs%refl_pref, g%domain)
3881 endif
3882
3883 cs%id_cg1 = register_diag_field('ocean_model', 'cn1', diag%axesT1, &
3884 time, 'First baroclinic mode (eigen) speed', 'm s-1', conversion=us%L_T_to_m_s)
3885 allocate(cs%id_cn(cs%nMode), source=-1)
3886 do m=1,cs%nMode
3887 write(var_name, '("cn_mode",i1)') m
3888 write(var_descript, '("Baroclinic (eigen) speed of mode ",i1)') m
3889 cs%id_cn(m) = register_diag_field('ocean_model',var_name, diag%axesT1, &
3890 time, var_descript, 'm s-1', conversion=us%L_T_to_m_s)
3891 call mom_mesg("Registering "//trim(var_name)//", Described as: "//var_descript, 5)
3892 enddo
3893
3894 ! Register maps of reflection parameters
3895 cs%id_refl_ang = register_diag_field('ocean_model', 'refl_angle', diag%axesT1, &
3896 time, 'Local angle of coastline/ridge/shelf with respect to equator', 'rad')
3897 cs%id_refl_pref = register_diag_field('ocean_model', 'refl_pref', diag%axesT1, &
3898 time, 'Partial reflection coefficients', '')
3899 cs%id_trans = register_diag_field('ocean_model', 'trans', diag%axesT1, &
3900 time, 'Partial transmission coefficients', '')
3901 cs%id_residual = register_diag_field('ocean_model', 'residual', diag%axesT1, &
3902 time, 'Residual of reflection and transmission coefficients', '')
3903 cs%id_dx_Cv = register_diag_field('ocean_model', 'dx_Cv', diag%axesT1, &
3904 time, 'North face unblocked width', 'm', conversion=us%L_to_m)
3905 cs%id_dy_Cu = register_diag_field('ocean_model', 'dy_Cu', diag%axesT1, &
3906 time, 'East face unblocked width', 'm', conversion=us%L_to_m)
3907 cs%id_land_mask = register_diag_field('ocean_model', 'land_mask', diag%axesT1, &
3908 time, 'Land mask', 'nondim')
3909 ! Output reflection parameters as diagnostics here (not needed every timestep)
3910 if (cs%id_refl_ang > 0) call post_data(cs%id_refl_ang, cs%refl_angle, cs%diag)
3911 if (cs%id_refl_pref > 0) call post_data(cs%id_refl_pref, cs%refl_pref, cs%diag)
3912 if (cs%id_trans > 0) call post_data(cs%id_trans, cs%trans, cs%diag)
3913 if (cs%id_residual > 0) call post_data(cs%id_residual, cs%residual, cs%diag)
3914 if (cs%id_dx_Cv > 0) call post_data(cs%id_dx_Cv, g%dx_Cv, cs%diag)
3915 if (cs%id_dy_Cu > 0) call post_data(cs%id_dy_Cu, g%dy_Cu, cs%diag)
3916 if (cs%id_land_mask > 0) call post_data(cs%id_land_mask, g%mask2dT, cs%diag)
3917
3918 ! Register 2-D energy density (summed over angles, freq, modes)
3919 cs%id_tot_En = register_diag_field('ocean_model', 'ITide_tot_En', diag%axesT1, &
3920 time, 'Internal tide total energy density', &
3921 'J m-2', conversion=hz2_t2_to_j_m2)
3922
3923 allocate(cs%id_itide_drag(cs%nFreq, cs%nMode), source=-1)
3924 allocate(cs%id_TKE_itidal_input(cs%nFreq), source=-1)
3925 do fr=1,cs%nFreq
3926 ! Register 2-D energy input into internal tides for each frequency
3927 write(var_name, '("TKE_itidal_input_freq",i1)') fr
3928 write(var_descript, '("a fraction of which goes into rays in frequency ",i1)') fr
3929
3930 cs%id_TKE_itidal_input(fr) = register_diag_field('ocean_model', var_name, diag%axesT1, &
3931 time, 'Conversion from barotropic to baroclinic tide, '//&
3932 var_descript, 'W m-2', conversion=hz2_t3_to_w_m2)
3933 enddo
3934 ! Register 2-D energy losses (summed over angles, freq, modes)
3935 cs%id_tot_leak_loss = register_diag_field('ocean_model', 'ITide_tot_leak_loss', diag%axesT1, &
3936 time, 'Internal tide energy loss to background drag', &
3937 'W m-2', conversion=hz2_t3_to_w_m2)
3938 cs%id_tot_quad_loss = register_diag_field('ocean_model', 'ITide_tot_quad_loss', diag%axesT1, &
3939 time, 'Internal tide energy loss to bottom drag', &
3940 'W m-2', conversion=hz2_t3_to_w_m2)
3941 cs%id_tot_itidal_loss = register_diag_field('ocean_model', 'ITide_tot_itidal_loss', diag%axesT1, &
3942 time, 'Internal tide energy loss to wave drag', &
3943 'W m-2', conversion=hz2_t3_to_w_m2)
3944 cs%id_tot_Froude_loss = register_diag_field('ocean_model', 'ITide_tot_Froude_loss', diag%axesT1, &
3945 time, 'Internal tide energy loss to wave breaking', &
3946 'W m-2', conversion=hz2_t3_to_w_m2)
3947 cs%id_tot_residual_loss = register_diag_field('ocean_model', 'ITide_tot_residual_loss', diag%axesT1, &
3948 time, 'Internal tide energy loss to residual on slopes', &
3949 'W m-2', conversion=hz2_t3_to_w_m2)
3950 cs%id_tot_allprocesses_loss = register_diag_field('ocean_model', 'ITide_tot_allprocesses_loss', diag%axesT1, &
3951 time, 'Internal tide energy loss summed over all processes', &
3952 'W m-2', conversion=hz2_t3_to_w_m2)
3953
3954 allocate(cs%id_En_mode(cs%nFreq,cs%nMode), source=-1)
3955 allocate(cs%id_En_ang_mode(cs%nFreq,cs%nMode), source=-1)
3956 allocate(cs%id_itidal_loss_mode(cs%nFreq,cs%nMode), source=-1)
3957 allocate(cs%id_leak_loss_mode(cs%nFreq,cs%nMode), source=-1)
3958 allocate(cs%id_quad_loss_mode(cs%nFreq,cs%nMode), source=-1)
3959 allocate(cs%id_Froude_loss_mode(cs%nFreq,cs%nMode), source=-1)
3960 allocate(cs%id_residual_loss_mode(cs%nFreq,cs%nMode), source=-1)
3961 allocate(cs%id_allprocesses_loss_mode(cs%nFreq,cs%nMode), source=-1)
3962 allocate(cs%id_itidal_loss_ang_mode(cs%nFreq,cs%nMode), source=-1)
3963 allocate(cs%id_Ub_mode(cs%nFreq,cs%nMode), source=-1)
3964 allocate(cs%id_Ustruct_mode(cs%nMode), source=-1)
3965 allocate(cs%id_Wstruct_mode(cs%nMode), source=-1)
3966 allocate(cs%id_int_w2_mode(cs%nMode), source=-1)
3967 allocate(cs%id_int_U2_mode(cs%nMode), source=-1)
3968 allocate(cs%id_int_N2w2_mode(cs%nMode), source=-1)
3969 allocate(cs%id_cp_mode(cs%nFreq,cs%nMode), source=-1)
3970
3971 allocate(angles(cs%NAngle), source=0.0)
3972 angle_size = (8.0*atan(1.0)) / (real(num_angle))
3973 do a=1,num_angle ; angles(a) = (real(a) - 1) * angle_size ; enddo
3974
3975 id_ang = diag_axis_init("angle", angles, "Radians", "N", "Angular Orientation of Fluxes")
3976 call define_axes_group(diag, (/ diag%axesT1%handles(1), diag%axesT1%handles(2), id_ang /), &
3977 axes_ang, is_h_point=.true.)
3978 do fr=1,cs%nFreq ; write(freq_name(fr), '("freq",i1)') fr ; enddo
3979 do m=1,cs%nMode ; do fr=1,cs%nFreq
3980 ! Register 2-D energy density (summed over angles) for each frequency and mode
3981 write(var_name, '("Itide_En_freq",i1,"_mode",i1)') fr, m
3982 write(var_descript, '("Internal tide energy density in frequency ",i1," mode ",i1)') fr, m
3983 cs%id_En_mode(fr,m) = register_diag_field('ocean_model', var_name, &
3984 diag%axesT1, time, var_descript, 'J m-2', conversion=hz2_t2_to_j_m2)
3985 call mom_mesg("Registering "//trim(var_name)//", Described as: "//var_descript, 5)
3986
3987 ! Register 3-D (i,j,a) energy density for each frequency and mode
3988 write(var_name, '("Itide_En_ang_freq",i1,"_mode",i1)') fr, m
3989 write(var_descript, '("Internal tide angular energy density in frequency ",i1," mode ",i1)') fr, m
3990 cs%id_En_ang_mode(fr,m) = register_diag_field('ocean_model', var_name, &
3991 axes_ang, time, var_descript, 'J m-2 band-1', conversion=hz2_t2_to_j_m2)
3992 call mom_mesg("Registering "//trim(var_name)//", Described as: "//var_descript, 5)
3993
3994 ! Register 2-D energy loss (summed over angles) for each frequency and mode
3995 ! wave-drag only
3996 write(var_name, '("Itide_wavedrag_loss_freq",i1,"_mode",i1)') fr, m
3997 write(var_descript, '("Internal tide energy loss due to wave-drag from frequency ",i1," mode ",i1)') fr, m
3998 cs%id_itidal_loss_mode(fr,m) = register_diag_field('ocean_model', var_name, &
3999 diag%axesT1, time, var_descript, 'W m-2', conversion=hz2_t3_to_w_m2)
4000 call mom_mesg("Registering "//trim(var_name)//", Described as: "//var_descript, 5)
4001 ! Leakage loss
4002 write(var_name, '("Itide_leak_loss_freq",i1,"_mode",i1)') fr, m
4003 write(var_descript, '("Internal tide energy loss due to leakage from frequency ",i1," mode ",i1)') fr, m
4004 cs%id_leak_loss_mode(fr,m) = register_diag_field('ocean_model', var_name, &
4005 diag%axesT1, time, var_descript, 'W m-2', conversion=hz2_t3_to_w_m2)
4006 call mom_mesg("Registering "//trim(var_name)//", Described as: "//var_descript, 5)
4007 ! Quad loss
4008 write(var_name, '("Itide_quad_loss_freq",i1,"_mode",i1)') fr, m
4009 write(var_descript, '("Internal tide energy quad loss from frequency ",i1," mode ",i1)') fr, m
4010 cs%id_quad_loss_mode(fr,m) = register_diag_field('ocean_model', var_name, &
4011 diag%axesT1, time, var_descript, 'W m-2', conversion=hz2_t3_to_w_m2)
4012 call mom_mesg("Registering "//trim(var_name)//", Described as: "//var_descript, 5)
4013 ! Froude loss
4014 write(var_name, '("Itide_froude_loss_freq",i1,"_mode",i1)') fr, m
4015 write(var_descript, '("Internal tide energy Froude loss from frequency ",i1," mode ",i1)') fr, m
4016 cs%id_froude_loss_mode(fr,m) = register_diag_field('ocean_model', var_name, &
4017 diag%axesT1, time, var_descript, 'W m-2', conversion=hz2_t3_to_w_m2)
4018 call mom_mesg("Registering "//trim(var_name)//", Described as: "//var_descript, 5)
4019 ! residual losses
4020 write(var_name, '("Itide_residual_loss_freq",i1,"_mode",i1)') fr, m
4021 write(var_descript, '("Internal tide energy residual loss from frequency ",i1," mode ",i1)') fr, m
4022 cs%id_residual_loss_mode(fr,m) = register_diag_field('ocean_model', var_name, &
4023 diag%axesT1, time, var_descript, 'W m-2', conversion=hz2_t3_to_w_m2)
4024 call mom_mesg("Registering "//trim(var_name)//", Described as: "//var_descript, 5)
4025 ! all loss processes
4026 write(var_name, '("Itide_allprocesses_loss_freq",i1,"_mode",i1)') fr, m
4027 write(var_descript, '("Internal tide energy loss due to all processes from frequency ",i1," mode ",i1)') fr, m
4028 cs%id_allprocesses_loss_mode(fr,m) = register_diag_field('ocean_model', var_name, &
4029 diag%axesT1, time, var_descript, 'W m-2', conversion=hz2_t3_to_w_m2)
4030 call mom_mesg("Registering "//trim(var_name)//", Described as: "//var_descript, 5)
4031
4032 ! Register 3-D (i,j,a) energy loss for each frequency and mode
4033 ! wave-drag only
4034 write(var_name, '("Itide_wavedrag_loss_ang_freq",i1,"_mode",i1)') fr, m
4035 write(var_descript, '("Internal tide energy loss due to wave-drag from frequency ",i1," mode ",i1)') fr, m
4036 cs%id_itidal_loss_ang_mode(fr,m) = register_diag_field('ocean_model', var_name, &
4037 axes_ang, time, var_descript, 'W m-2 band-1', conversion=hz2_t3_to_w_m2)
4038 call mom_mesg("Registering "//trim(var_name)//", Described as: "//var_descript, 5)
4039
4040 ! Register 2-D period-averaged near-bottom horizontal velocity for each frequency and mode
4041 write(var_name, '("Itide_Ub_freq",i1,"_mode",i1)') fr, m
4042 write(var_descript, '("Near-bottom horizontal velocity for frequency ",i1," mode ",i1)') fr, m
4043 cs%id_Ub_mode(fr,m) = register_diag_field('ocean_model', var_name, &
4044 diag%axesT1, time, var_descript, 'm s-1', conversion=us%L_T_to_m_s)
4045 call mom_mesg("Registering "//trim(var_name)//", Described as: "//var_descript, 5)
4046
4047 ! Register 2-D horizontal phase velocity for each frequency and mode
4048 write(var_name, '("Itide_cp_freq",i1,"_mode",i1)') fr, m
4049 write(var_descript, '("Horizontal phase velocity for frequency ",i1," mode ",i1)') fr, m
4050 cs%id_cp_mode(fr,m) = register_diag_field('ocean_model', var_name, &
4051 diag%axesT1, time, var_descript, 'm s-1', conversion=us%L_T_to_m_s)
4052 call mom_mesg("Registering "//trim(var_name)//", Described as: "//var_descript, 5)
4053
4054 ! Register 2-D drag scale used for quadratic bottom drag for each frequency and mode
4055 write(var_name, '("ITide_drag_freq",i1,"_mode",i1)') fr, m
4056 write(var_descript, '("Interior and bottom drag int tide decay timescale in frequency ",i1, " mode ",i1)') fr, m
4057
4058 cs%id_itide_drag(fr,m) = register_diag_field('ocean_model', var_name, diag%axesT1, time, &
4059 's-1', conversion=us%s_to_T)
4060 enddo ; enddo
4061
4062
4063 do m=1,cs%nMode
4064
4065 ! Register 3-D internal tide horizontal velocity profile for each mode
4066 write(var_name, '("Itide_Ustruct","_mode",i1)') m
4067 write(var_descript, '("horizontal velocity profile for mode ",i1)') m
4068 cs%id_Ustruct_mode(m) = register_diag_field('ocean_model', var_name, &
4069 diag%axesTl, time, var_descript, 'm-1', conversion=us%m_to_L)
4070 call mom_mesg("Registering "//trim(var_name)//", Described as: "//var_descript, 5)
4071
4072 ! Register 3-D internal tide vertical velocity profile for each mode
4073 write(var_name, '("Itide_Wstruct","_mode",i1)') m
4074 write(var_descript, '("vertical velocity profile for mode ",i1)') m
4075 cs%id_Wstruct_mode(m) = register_diag_field('ocean_model', var_name, &
4076 diag%axesTi, time, var_descript, 'nondim')
4077 call mom_mesg("Registering "//trim(var_name)//", Described as: "//var_descript, 5)
4078
4079 write(var_name, '("Itide_int_w2","_mode",i1)') m
4080 write(var_descript, '("integral of w2 for mode ",i1)') m
4081 cs%id_int_w2_mode(m) = register_diag_field('ocean_model', var_name, &
4082 diag%axesT1, time, var_descript, 'm', conversion=gv%H_to_m)
4083 call mom_mesg("Registering "//trim(var_name)//", Described as: "//var_descript, 5)
4084
4085 write(var_name, '("Itide_int_U2","_mode",i1)') m
4086 write(var_descript, '("integral of U2 for mode ",i1)') m
4087 cs%id_int_U2_mode(m) = register_diag_field('ocean_model', var_name, &
4088 diag%axesT1, time, var_descript, 'm-1', conversion=us%m_to_Z*gv%H_to_Z)
4089 call mom_mesg("Registering "//trim(var_name)//", Described as: "//var_descript, 5)
4090
4091 write(var_name, '("Itide_int_N2w2","_mode",i1)') m
4092 write(var_descript, '("integral of N2w2 for mode ",i1)') m
4093 cs%id_int_N2w2_mode(m) = register_diag_field('ocean_model', var_name, &
4094 diag%axesT1, time, var_descript, 'm s-2', conversion=gv%H_to_m*us%s_to_T**2)
4095 call mom_mesg("Registering "//trim(var_name)//", Described as: "//var_descript, 5)
4096
4097 enddo
4098
4099 ! Initialize the module that calculates the wave speeds.
4100 call wave_speed_init(cs%wave_speed, gv, c1_thresh=igw_c1_thresh, &
4101 om4_remap_via_sub_cells=om4_remap_via_sub_cells)
4102
4103end subroutine internal_tides_init
4104
4105!> This subroutine deallocates the memory associated with the internal tides control structure
4106subroutine internal_tides_end(CS)
4107 type(int_tide_cs), intent(inout) :: cs !< Internal tide control structure
4108
4109 if (allocated(cs%En)) deallocate(cs%En)
4110 if (allocated(cs%frequency)) deallocate(cs%frequency)
4111 if (allocated(cs%id_En_mode)) deallocate(cs%id_En_mode)
4112 if (allocated(cs%id_Ub_mode)) deallocate(cs%id_Ub_mode)
4113 if (allocated(cs%id_cp_mode)) deallocate(cs%id_cp_mode)
4114 if (allocated(cs%id_Ustruct_mode)) deallocate(cs%id_Ustruct_mode)
4115 if (allocated(cs%id_Wstruct_mode)) deallocate(cs%id_Wstruct_mode)
4116 if (allocated(cs%id_int_w2_mode)) deallocate(cs%id_int_w2_mode)
4117 if (allocated(cs%id_int_U2_mode)) deallocate(cs%id_int_U2_mode)
4118 if (allocated(cs%id_int_N2w2_mode)) deallocate(cs%id_int_N2w2_mode)
4119
4120end subroutine internal_tides_end
4121
4122end module mom_internal_tides