MOM6
MOM_internal_tides.F90
1 !> Subroutines that use the ray-tracing equations to propagate the internal tide energy density.
2 !!
3 !! \author Benjamin Mater & Robert Hallberg, 2015
4 module mom_internal_tides
5 
6 ! This file is part of MOM6. See LICENSE.md for the license.
7 
8 use mom_debugging, only : is_nan
9 use mom_diag_mediator, only : post_data, query_averaging_enabled, diag_axis_init
10 use mom_diag_mediator, only : register_diag_field, diag_ctrl, safe_alloc_ptr
11 use mom_diag_mediator, only : axes_grp, define_axes_group
12 use mom_domains, only : agrid, to_south, to_west, to_all
13 use mom_domains, only : create_group_pass, do_group_pass, pass_var
14 use mom_domains, only : group_pass_type, start_group_pass, complete_group_pass
15 use mom_error_handler, only : mom_error, fatal, warning, mom_mesg, is_root_pe
16 use mom_file_parser, only : read_param, get_param, log_param, log_version, param_file_type
17 use mom_grid, only : ocean_grid_type
18 use mom_io, only : slasher, vardesc, mom_read_data
19 use mom_restart, only : register_restart_field, mom_restart_cs, restart_init, save_restart
20 use mom_spatial_means, only : global_area_mean
21 use mom_time_manager, only : time_type, time_type_to_real, operator(+), operator(/), operator(-)
22 use mom_unit_scaling, only : unit_scale_type
23 use mom_variables, only : surface, thermo_var_ptrs
24 use mom_verticalgrid, only : verticalgrid_type
25 use mom_wave_structure, only: wave_structure_init, wave_structure, wave_structure_cs
26 
27 !use, intrinsic :: IEEE_ARITHMETIC
28 
29 implicit none ; private
30 
31 #include <MOM_memory.h>
32 
33 public propagate_int_tide !, register_int_tide_restarts
34 public internal_tides_init, internal_tides_end
35 public get_lowmode_loss
36 
37 !> This control structure has parameters for the MOM_internal_tides module
38 type, public :: int_tide_cs ; private
39  logical :: do_int_tides !< If true, use the internal tide code.
40  integer :: nfreq = 0 !< The number of internal tide frequency bands
41  integer :: nmode = 1 !< The number of internal tide vertical modes
42  integer :: nangle = 24 !< The number of internal tide angular orientations
43  integer :: energized_angle = -1 !< If positive, only this angular band is energized for debugging purposes
44  logical :: corner_adv !< If true, use a corner advection rather than PPM.
45  logical :: upwind_1st !< If true, use a first-order upwind scheme.
46  logical :: simple_2nd !< If true, use a simple second order (arithmetic mean) interpolation
47  !! of the edge values instead of the higher order interpolation.
48  logical :: vol_cfl !< If true, use the ratio of the open face lengths to the tracer cell
49  !! areas when estimating CFL numbers. Without aggress_adjust,
50  !! the default is false; it is always true with aggress_adjust.
51  logical :: use_ppmang !< If true, use PPM for advection of energy in angular space.
52 
53  real, allocatable, dimension(:,:) :: refl_angle
54  !< local coastline/ridge/shelf angles read from file
55  ! (could be in G control structure)
56  real :: nullangle = -999.9 !< placeholder value in cells with no reflection
57  real, allocatable, dimension(:,:) :: refl_pref
58  !< partial reflection coeff for each "coast cell"
59  ! (could be in G control structure)
60  logical, allocatable, dimension(:,:) :: refl_pref_logical
61  !< true if reflecting cell with partial reflection
62  ! (could be in G control structure)
63  logical, allocatable, dimension(:,:) :: refl_dbl
64  !< identifies reflection cells where double reflection
65  !! is possible (i.e. ridge cells)
66  ! (could be in G control structure)
67  real, allocatable, dimension(:,:,:,:) :: cp
68  !< horizontal phase speed [m s-1]
69  real, allocatable, dimension(:,:,:,:,:) :: tke_leak_loss
70  !< energy lost due to misc background processes [W m-2]
71  real, allocatable, dimension(:,:,:,:,:) :: tke_quad_loss
72  !< energy lost due to quadratic bottom drag [W m-2]
73  real, allocatable, dimension(:,:,:,:,:) :: tke_froude_loss
74  !< energy lost due to wave breaking [W m-2]
75  real, allocatable, dimension(:,:) :: tke_itidal_loss_fixed
76  !< fixed part of the energy lost due to small-scale drag
77  !! [kg Z-2 ~> kg m-2] here; will be multiplied by N and En to get into [W m-2]
78  real, allocatable, dimension(:,:,:,:,:) :: tke_itidal_loss
79  !< energy lost due to small-scale wave drag [W m-2]
80  real, allocatable, dimension(:,:) :: tot_leak_loss !< Energy loss rates due to misc bakground processes,
81  !! summed over angle, frequency and mode [W m-2]
82  real, allocatable, dimension(:,:) :: tot_quad_loss !< Energy loss rates due to quadratic bottom drag,
83  !! summed over angle, frequency and mode [W m-2]
84  real, allocatable, dimension(:,:) :: tot_itidal_loss !< Energy loss rates due to small-scale drag,
85  !! summed over angle, frequency and mode [W m-2]
86  real, allocatable, dimension(:,:) :: tot_froude_loss !< Energy loss rates due to wave breaking,
87  !! summed over angle, frequency and mode [W m-2]
88  real, allocatable, dimension(:,:) :: tot_allprocesses_loss !< Energy loss rates due to all processes,
89  !! summed over angle, frequency and mode [W m-2]
90  real :: q_itides !< fraction of local dissipation [nondim]
91  real :: en_sum !< global sum of energy for use in debugging
92  type(time_type), pointer :: time => null() !< A pointer to the model's clock.
93  character(len=200) :: inputdir !< directory to look for coastline angle file
94  real :: decay_rate !< A constant rate at which internal tide energy is
95  !! lost to the interior ocean internal wave field.
96  real :: cdrag !< The bottom drag coefficient [nondim].
97  logical :: apply_background_drag
98  !< If true, apply a drag due to background processes as a sink.
99  logical :: apply_bottom_drag
100  !< If true, apply a quadratic bottom drag as a sink.
101  logical :: apply_wave_drag
102  !< If true, apply scattering due to small-scale roughness as a sink.
103  logical :: apply_froude_drag
104  !< If true, apply wave breaking as a sink.
105  real, dimension(:,:,:,:,:), pointer :: en => null()
106  !< The internal wave energy density as a function of (i,j,angle,frequency,mode)
107  real, dimension(:,:,:), pointer :: en_restart => null()
108  !< The internal wave energy density as a function of (i,j,angle); temporary for restart
109  real, allocatable, dimension(:) :: frequency !< The frequency of each band [s-1].
110 
111  type(diag_ctrl), pointer :: diag => null() !< A structure that is used to regulate the
112  !! timing of diagnostic output.
113  type(wave_structure_cs), pointer :: wave_structure_csp => null()
114  !< A pointer to the wave_structure module control structure
115 
116  !>@{ Diag handles
117  ! Diag handles relevant to all modes, frequencies, and angles
118  integer :: id_tot_en = -1, id_tke_itidal_input = -1, id_itide_drag = -1
119  integer :: id_refl_pref = -1, id_refl_ang = -1, id_land_mask = -1
120  integer :: id_dx_cv = -1, id_dy_cu = -1
121  ! Diag handles considering: sums over all modes, frequencies, and angles
122  integer :: id_tot_leak_loss = -1, id_tot_quad_loss = -1, id_tot_itidal_loss = -1
123  integer :: id_tot_froude_loss = -1, id_tot_allprocesses_loss = -1
124  ! Diag handles considering: all modes & freqs; summed over angles
125  integer, allocatable, dimension(:,:) :: &
126  id_en_mode, &
127  id_itidal_loss_mode, &
128  id_allprocesses_loss_mode, &
129  id_ub_mode, &
130  id_cp_mode
131  ! Diag handles considering: all modes, freqs, and angles
132  integer, allocatable, dimension(:,:) :: &
133  id_en_ang_mode, &
134  id_itidal_loss_ang_mode
135  !!@}
136 
137 end type int_tide_cs
138 
139 !> A structure with the active energy loop bounds.
140 type :: loop_bounds_type ; private
141  !>@{ The active loop bounds
142  integer :: ish, ieh, jsh, jeh
143  !!@}
144 end type loop_bounds_type
145 
146 contains
147 
148 !> Calls subroutines in this file that are needed to refract, propagate,
149 !! and dissipate energy density of the internal tide.
150 subroutine propagate_int_tide(h, tv, cn, TKE_itidal_input, vel_btTide, Nb, dt, &
151  G, GV, US, CS)
152  type(ocean_grid_type), intent(inout) :: g !< The ocean's grid structure.
153  type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure.
154  type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
155  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
156  intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2]
157  type(thermo_var_ptrs), intent(in) :: tv !< Pointer to thermodynamic variables
158  !! (needed for wave structure).
159  real, dimension(SZI_(G),SZJ_(G)), intent(in) :: tke_itidal_input !< The energy input to the
160  !! internal waves [W m-2].
161  real, dimension(SZI_(G),SZJ_(G)), intent(in) :: vel_bttide !< Barotropic velocity read
162  !! from file [m s-1].
163  real, dimension(SZI_(G),SZJ_(G)), intent(in) :: nb !< Near-bottom buoyancy frequency [s-1].
164  real, intent(in) :: dt !< Length of time over which these fluxes
165  !! will be applied [s].
166  type(int_tide_cs), pointer :: cs !< The control structure returned by a
167  !! previous call to int_tide_init.
168  real, dimension(SZI_(G),SZJ_(G),CS%nMode), &
169  intent(in) :: cn !< The internal wave speeds of each mode [m s-1].
170  ! Local variables
171  real, dimension(SZI_(G),SZJ_(G),2) :: &
172  test
173  real, dimension(SZI_(G),SZJ_(G),CS%nFreq,CS%nMode) :: &
174  tot_en_mode, & ! energy summed over angles only
175  ub, umax ! near-bottom & max horizontal velocity of wave (modal)
176  real, dimension(SZI_(G),SZJB_(G)) :: &
177  flux_heat_y, &
178  flux_prec_y
179  real, dimension(SZI_(G),SZJ_(G)) :: &
180  tot_en, & ! energy summed over angles, modes, frequencies
181  tot_leak_loss, tot_quad_loss, tot_itidal_loss, tot_froude_loss, tot_allprocesses_loss, &
182  ! energy loss rates summed over angle, freq, and mode
183  drag_scale, & ! bottom drag scale, s-1
184  itidal_loss_mode, allprocesses_loss_mode
185  ! energy loss rates for a given mode and frequency (summed over angles)
186  real :: frac_per_sector, f2, i_rho0, i_d_here, freq2, kmag2
187  real :: c_phase, loss_rate, fr2_max
188  real, parameter :: cn_subro = 1e-100 ! to prevent division by zero
189  real :: en_new, en_check ! for debugging
190  real :: en_initial, delta_e_check ! for debugging
191  real :: tke_froude_loss_check, tke_froude_loss_tot ! for debugging
192  character(len=160) :: mesg ! The text of an error message
193  integer :: a, m, fr, i, j, is, ie, js, je, isd, ied, jsd, jed, nangle, nzm
194  integer :: id_g, jd_g ! global (decomp-invar) indices (for debugging)
195  type(group_pass_type), save :: pass_test, pass_en
196 
197  if (.not.associated(cs)) return
198  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
199  isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed ; nangle = cs%NAngle
200  i_rho0 = 1.0 / gv%Rho0
201 
202  ! Set the wave speeds for the modes, using cg(n) ~ cg(1)/n.**********************
203  ! This is wrong, of course, but it works reasonably in some cases.
204  ! Uncomment if wave_speed is not used to calculate the true values (BDM).
205  !do m=1,CS%nMode ; do j=jsd,jed ; do i=isd,ied
206  ! cn(i,j,m) = cn(i,j,1) / real(m)
207  !enddo ; enddo ; enddo
208 
209  ! Add the forcing.***************************************************************
210  if (cs%energized_angle <= 0) then
211  frac_per_sector = 1.0 / real(cs%nAngle * cs%nMode * cs%nFreq)
212  do m=1,cs%nMode ; do fr=1,cs%nFreq ; do a=1,cs%nAngle ; do j=js,je ; do i=is,ie
213  f2 = 0.25*us%s_to_T**2*((g%CoriolisBu(i,j)**2 + g%CoriolisBu(i-1,j-1)**2) + &
214  (g%CoriolisBu(i-1,j)**2 + g%CoriolisBu(i,j-1)**2))
215  if (cs%frequency(fr)**2 > f2) &
216  cs%En(i,j,a,fr,m) = cs%En(i,j,a,fr,m) + &
217  dt*frac_per_sector*(1-cs%q_itides)*tke_itidal_input(i,j)
218  enddo ; enddo ; enddo ; enddo ; enddo
219  elseif (cs%energized_angle <= cs%nAngle) then
220  frac_per_sector = 1.0 / real(cs%nMode * cs%nFreq)
221  a = cs%energized_angle
222  do m=1,cs%nMode ; do fr=1,cs%nFreq ; do j=js,je ; do i=is,ie
223  f2 = 0.25*us%s_to_T**2*((g%CoriolisBu(i,j)**2 + g%CoriolisBu(i-1,j-1)**2) + &
224  (g%CoriolisBu(i-1,j)**2 + g%CoriolisBu(i,j-1)**2))
225  if (cs%frequency(fr)**2 > f2) &
226  cs%En(i,j,a,fr,m) = cs%En(i,j,a,fr,m) + &
227  dt*frac_per_sector**(1-cs%q_itides)*tke_itidal_input(i,j)
228  enddo ; enddo ; enddo ; enddo
229  else
230  call mom_error(warning, "Internal tide energy is being put into a angular "//&
231  "band that does not exist.")
232  endif
233 
234  ! Pass a test vector to check for grid rotation in the halo updates.
235  do j=jsd,jed ; do i=isd,ied ; test(i,j,1) = 1.0 ; test(i,j,2) = 0.0 ; enddo ; enddo
236  do m=1,cs%nMode ; do fr=1,cs%nFreq
237  call create_group_pass(pass_en, cs%En(:,:,:,fr,m), g%domain)
238  enddo ; enddo
239  call create_group_pass(pass_test, test(:,:,1), test(:,:,2), g%domain, stagger=agrid)
240  call start_group_pass(pass_test, g%domain)
241 
242  ! Apply half the refraction.
243  do m=1,cs%nMode ; do fr=1,cs%nFreq
244  call refract(cs%En(:,:,:,fr,m), cn(:,:,m), cs%frequency(fr), 0.5*dt, g, us, cs%nAngle, cs%use_PPMang)
245  enddo ; enddo
246 
247  ! Check for En<0 - for debugging, delete later
248  do m=1,cs%NMode ; do fr=1,cs%Nfreq ; do a=1,cs%nAngle
249  do j=js,je ; do i=is,ie
250  if (cs%En(i,j,a,fr,m)<0.0) then
251  id_g = i + g%idg_offset ; jd_g = j + g%jdg_offset ! for debugging
252  write(mesg,*) 'After first refraction: En<0.0 at ig=', id_g, ', jg=', jd_g, &
253  'En=',cs%En(i,j,a,fr,m)
254  call mom_error(warning, "propagate_int_tide: "//trim(mesg))
255  cs%En(i,j,a,fr,m) = 0.0
256 ! call MOM_error(FATAL, "propagate_int_tide: stopped due to negative energy.")
257  endif
258  enddo ; enddo
259  enddo ; enddo ; enddo
260 
261  call do_group_pass(pass_en, g%domain)
262 
263  call complete_group_pass(pass_test, g%domain)
264 
265  ! Rotate points in the halos as necessary.
266  call correct_halo_rotation(cs%En, test, g, cs%nAngle)
267 
268  ! Propagate the waves.
269  do m=1,cs%NMode ; do fr=1,cs%Nfreq
270  call propagate(cs%En(:,:,:,fr,m), cn(:,:,m), cs%frequency(fr), dt, g, us, cs, cs%NAngle)
271  enddo ; enddo
272 
273  ! Check for En<0 - for debugging, delete later
274  do m=1,cs%NMode ; do fr=1,cs%Nfreq ; do a=1,cs%nAngle
275  do j=js,je ; do i=is,ie
276  if (cs%En(i,j,a,fr,m)<0.0) then
277  id_g = i + g%idg_offset ; jd_g = j + g%jdg_offset
278  if (abs(cs%En(i,j,a,fr,m))>1.0) then ! only print if large
279  write(mesg,*) 'After propagation: En<0.0 at ig=', id_g, ', jg=', jd_g, &
280  'En=', cs%En(i,j,a,fr,m)
281  call mom_error(warning, "propagate_int_tide: "//trim(mesg))
282 ! call MOM_error(FATAL, "propagate_int_tide: stopped due to negative energy.")
283  endif
284  cs%En(i,j,a,fr,m) = 0.0
285  endif
286  enddo ; enddo
287  enddo ; enddo ; enddo
288 
289  ! Apply the other half of the refraction.
290  do m=1,cs%NMode ; do fr=1,cs%Nfreq
291  call refract(cs%En(:,:,:,fr,m), cn(:,:,m), cs%frequency(fr), 0.5*dt, g, us, cs%NAngle, cs%use_PPMang)
292  enddo ; enddo
293 
294  ! Check for En<0 - for debugging, delete later
295  do m=1,cs%NMode ; do fr=1,cs%Nfreq ; do a=1,cs%nAngle
296  do j=js,je ; do i=is,ie
297  if (cs%En(i,j,a,fr,m)<0.0) then
298  id_g = i + g%idg_offset ; jd_g = j + g%jdg_offset ! for debugging
299  write(mesg,*) 'After second refraction: En<0.0 at ig=', id_g, ', jg=', jd_g, &
300  'En=',cs%En(i,j,a,fr,m)
301  call mom_error(warning, "propagate_int_tide: "//trim(mesg))
302  cs%En(i,j,a,fr,m) = 0.0
303 ! call MOM_error(FATAL, "propagate_int_tide: stopped due to negative energy.")
304  endif
305  enddo ; enddo
306  enddo ; enddo ; enddo
307 
308  ! Apply various dissipation mechanisms.
309  if (cs%apply_background_drag .or. cs%apply_bottom_drag &
310  .or. cs%apply_wave_drag .or. cs%apply_Froude_drag &
311  .or. (cs%id_tot_En > 0)) then
312  tot_en(:,:) = 0.0
313  tot_en_mode(:,:,:,:) = 0.0
314  do m=1,cs%NMode ; do fr=1,cs%Nfreq
315  do j=jsd,jed ; do i=isd,ied ; do a=1,cs%nAngle
316  tot_en(i,j) = tot_en(i,j) + cs%En(i,j,a,fr,m)
317  tot_en_mode(i,j,fr,m) = tot_en_mode(i,j,fr,m) + cs%En(i,j,a,fr,m)
318  enddo ; enddo ; enddo
319  enddo ; enddo
320  endif
321 
322  ! Extract the energy for mixing due to misc. processes (background leakage)------
323  if (cs%apply_background_drag) then
324  do m=1,cs%nMode ; do fr=1,cs%nFreq ; do a=1,cs%nAngle ; do j=jsd,jed ; do i=isd,ied
325  ! Calculate loss rate and apply loss over the time step ; apply the same drag timescale
326  ! to each En component (technically not correct; fix later)
327  cs%TKE_leak_loss(i,j,a,fr,m) = cs%En(i,j,a,fr,m) * cs%decay_rate ! loss rate [Wm-2]
328  cs%En(i,j,a,fr,m) = cs%En(i,j,a,fr,m) / (1.0 + dt *cs%decay_rate) ! implicit update
329  enddo ; enddo ; enddo ; enddo ; enddo
330  endif
331  ! Check for En<0 - for debugging, delete later
332  do m=1,cs%NMode ; do fr=1,cs%Nfreq ; do a=1,cs%nAngle
333  do j=js,je ; do i=is,ie
334  if (cs%En(i,j,a,fr,m)<0.0) then
335  id_g = i + g%idg_offset ; jd_g = j + g%jdg_offset ! for debugging
336  write(mesg,*) 'After leak loss: En<0.0 at ig=', id_g, ', jg=', jd_g, &
337  'En=',cs%En(i,j,a,fr,m)
338  call mom_error(warning, "propagate_int_tide: "//trim(mesg), all_print=.true.)
339  cs%En(i,j,a,fr,m) = 0.0
340 ! call MOM_error(FATAL, "propagate_int_tide: stopped due to negative energy.")
341  endif
342  enddo ; enddo
343  enddo ; enddo ; enddo
344 
345  ! Extract the energy for mixing due to bottom drag-------------------------------
346  if (cs%apply_bottom_drag) then
347  do j=jsd,jed ; do i=isd,ied
348  ! Note the 1 m dimensional scale here. Should this be a parameter?
349  i_d_here = 1.0 / (us%Z_to_m*max(g%bathyT(i,j), 1.0*us%m_to_Z))
350  drag_scale(i,j) = cs%cdrag * sqrt(max(0.0, vel_bttide(i,j)**2 + &
351  tot_en(i,j) * i_rho0 * i_d_here)) * i_d_here
352  enddo ; enddo
353  do m=1,cs%nMode ; do fr=1,cs%nFreq ; do a=1,cs%nAngle ; do j=jsd,jed ; do i=isd,ied
354  ! Calculate loss rate and apply loss over the time step ; apply the same drag timescale
355  ! to each En component (technically not correct; fix later)
356  cs%TKE_quad_loss(i,j,a,fr,m) = cs%En(i,j,a,fr,m) * drag_scale(i,j) ! loss rate
357  cs%En(i,j,a,fr,m) = cs%En(i,j,a,fr,m) / (1.0 + dt *drag_scale(i,j)) ! implicit update
358  enddo ; enddo ; enddo ; enddo ; enddo
359  endif
360  ! Check for En<0 - for debugging, delete later
361  do m=1,cs%NMode ; do fr=1,cs%Nfreq ; do a=1,cs%nAngle
362  do j=js,je ; do i=is,ie
363  if (cs%En(i,j,a,fr,m)<0.0) then
364  id_g = i + g%idg_offset ; jd_g = j + g%jdg_offset ! for debugging
365  write(mesg,*) 'After bottom loss: En<0.0 at ig=', id_g, ', jg=', jd_g, &
366  'En=',cs%En(i,j,a,fr,m)
367  call mom_error(warning, "propagate_int_tide: "//trim(mesg), all_print=.true.)
368  cs%En(i,j,a,fr,m) = 0.0
369 ! call MOM_error(FATAL, "propagate_int_tide: stopped due to negative energy.")
370  !stop
371  endif
372  enddo ; enddo
373  enddo ; enddo ; enddo
374 
375  ! Extract the energy for mixing due to scattering (wave-drag)--------------------
376  ! still need to allow a portion of the extracted energy to go to higher modes.
377  ! First, find velocity profiles
378  if (cs%apply_wave_drag .or. cs%apply_Froude_drag) then
379  do m=1,cs%NMode ; do fr=1,cs%Nfreq
380  ! Calculate modal structure for given mode and frequency
381  call wave_structure(h, tv, g, gv, us, cn(:,:,m), m, cs%frequency(fr), &
382  cs%wave_structure_CSp, tot_en_mode(:,:,fr,m), full_halos=.true.)
383  ! Pick out near-bottom and max horizontal baroclinic velocity values at each point
384  do j=jsd,jed ; do i=isd,ied
385  id_g = i + g%idg_offset ; jd_g = j + g%jdg_offset ! for debugging
386  nzm = cs%wave_structure_CSp%num_intfaces(i,j)
387  ub(i,j,fr,m) = cs%wave_structure_CSp%Uavg_profile(i,j,nzm)
388  umax(i,j,fr,m) = maxval(cs%wave_structure_CSp%Uavg_profile(i,j,1:nzm))
389  enddo ; enddo ! i-loop, j-loop
390  enddo ; enddo ! fr-loop, m-loop
391  endif ! apply_wave or _Froude_drag (Ub or Umax needed)
392  ! Finally, apply loss
393  if (cs%apply_wave_drag) then
394  ! Calculate loss rate and apply loss over the time step
395  call itidal_lowmode_loss(g, us, cs, nb, ub, cs%En, cs%TKE_itidal_loss_fixed, &
396  cs%TKE_itidal_loss, dt, full_halos=.false.)
397  endif
398  ! Check for En<0 - for debugging, delete later
399  do m=1,cs%NMode ; do fr=1,cs%Nfreq ; do a=1,cs%nAngle
400  do j=js,je ; do i=is,ie
401  if (cs%En(i,j,a,fr,m)<0.0) then
402  id_g = i + g%idg_offset ; jd_g = j + g%jdg_offset ! for debugging
403  write(mesg,*) 'After wave drag loss: En<0.0 at ig=', id_g, ', jg=', jd_g, &
404  'En=',cs%En(i,j,a,fr,m)
405  call mom_error(warning, "propagate_int_tide: "//trim(mesg), all_print=.true.)
406  cs%En(i,j,a,fr,m) = 0.0
407 ! call MOM_error(FATAL, "propagate_int_tide: stopped due to negative energy.")
408  endif
409  enddo ; enddo
410  enddo ; enddo ; enddo
411 
412  ! Extract the energy for mixing due to wave breaking-----------------------------
413  if (cs%apply_Froude_drag) then
414  ! Pick out maximum baroclinic velocity values; calculate Fr=max(u)/cg
415  do m=1,cs%NMode ; do fr=1,cs%Nfreq
416  freq2 = cs%frequency(fr)**2
417  do j=jsd,jed ; do i=isd,ied
418  id_g = i + g%idg_offset ; jd_g = j + g%jdg_offset ! for debugging
419  ! Calculate horizontal phase velocity magnitudes
420  f2 = 0.25*us%s_to_T**2*((g%CoriolisBu(i,j)**2 + g%CoriolisBu(i-1,j-1)**2) + &
421  (g%CoriolisBu(i-1,j)**2 + g%CoriolisBu(i,j-1)**2))
422  kmag2 = (freq2 - f2) / (cn(i,j,m)**2 + cn_subro**2)
423  c_phase = 0.0
424  if (kmag2 > 0.0) then
425  c_phase = sqrt(freq2/kmag2)
426  nzm = cs%wave_structure_CSp%num_intfaces(i,j)
427  fr2_max = (umax(i,j,fr,m)/c_phase)**2
428  ! Dissipate energy if Fr>1; done here with an arbitrary time scale
429  if (fr2_max > 1.0) then
430  en_initial = sum(cs%En(i,j,:,fr,m)) ! for debugging
431  ! Calculate effective decay rate [s-1] if breaking occurs over a time step
432  loss_rate = (1/fr2_max - 1.0)/dt
433  do a=1,cs%nAngle
434  ! Determine effective dissipation rate (Wm-2)
435  cs%TKE_Froude_loss(i,j,a,fr,m) = cs%En(i,j,a,fr,m) * abs(loss_rate)
436  ! Update energy
437  en_new = cs%En(i,j,a,fr,m)/fr2_max ! for debugging
438  en_check = cs%En(i,j,a,fr,m) - cs%TKE_Froude_loss(i,j,a,fr,m)*dt ! for debugging
439  ! Re-scale (reduce) energy due to breaking
440  cs%En(i,j,a,fr,m) = cs%En(i,j,a,fr,m)/fr2_max
441  ! Check (for debugging only)
442  if (abs(en_new - en_check) > 1e-10) then
443  call mom_error(warning, "MOM_internal_tides: something is wrong with Fr-breaking.", &
444  all_print=.true.)
445  write(mesg,*) "En_new=", en_new , "En_check=", en_check
446  call mom_error(warning, "MOM_internal_tides: "//trim(mesg), all_print=.true.)
447  endif
448  enddo
449  ! Check (for debugging)
450  delta_e_check = en_initial - sum(cs%En(i,j,:,fr,m))
451  tke_froude_loss_check = abs(delta_e_check)/dt
452  tke_froude_loss_tot = sum(cs%TKE_Froude_loss(i,j,:,fr,m))
453  if (abs(tke_froude_loss_check - tke_froude_loss_tot) > 1e-10) then
454  call mom_error(warning, "MOM_internal_tides: something is wrong with Fr energy update.", &
455  all_print=.true.)
456  write(mesg,*) "TKE_Froude_loss_check=", tke_froude_loss_check, &
457  "TKE_Froude_loss_tot=", tke_froude_loss_tot
458  call mom_error(warning, "MOM_internal_tides: "//trim(mesg), all_print=.true.)
459  endif
460  endif ! Fr2>1
461  endif ! Kmag2>0
462  cs%cp(i,j,fr,m) = c_phase
463  enddo ; enddo
464  enddo ; enddo
465  endif
466  ! Check for En<0 - for debugging, delete later
467  do m=1,cs%NMode ; do fr=1,cs%Nfreq ; do a=1,cs%nAngle
468  do j=js,je ; do i=is,ie
469  if (cs%En(i,j,a,fr,m)<0.0) then
470  id_g = i + g%idg_offset ; jd_g = j + g%jdg_offset
471  write(mesg,*) 'After Froude loss: En<0.0 at ig=', id_g, ', jg=', jd_g, &
472  'En=',cs%En(i,j,a,fr,m)
473  call mom_error(warning, "propagate_int_tide: "//trim(mesg), all_print=.true.)
474  cs%En(i,j,a,fr,m) = 0.0
475 ! call MOM_error(FATAL, "propagate_int_tide: stopped due to negative energy.")
476  !stop
477  endif
478  enddo ; enddo
479  enddo ; enddo ; enddo
480 
481  ! Check for energy conservation on computational domain.*************************
482  do m=1,cs%NMode ; do fr=1,cs%Nfreq
483  call sum_en(g,cs,cs%En(:,:,:,fr,m),'prop_int_tide')
484  enddo ; enddo
485 
486  ! Output diagnostics.************************************************************
487  if (query_averaging_enabled(cs%diag)) then
488  ! Output two-dimensional diagnostistics
489  if (cs%id_tot_En > 0) call post_data(cs%id_tot_En, tot_en, cs%diag)
490  if (cs%id_itide_drag > 0) call post_data(cs%id_itide_drag, drag_scale, cs%diag)
491  if (cs%id_TKE_itidal_input > 0) call post_data(cs%id_TKE_itidal_input, &
492  tke_itidal_input, cs%diag)
493 
494  ! Output 2-D energy density (summed over angles) for each freq and mode
495  do m=1,cs%NMode ; do fr=1,cs%Nfreq ; if (cs%id_En_mode(fr,m) > 0) then
496  tot_en(:,:) = 0.0
497  do a=1,cs%nAngle ; do j=js,je ; do i=is,ie
498  tot_en(i,j) = tot_en(i,j) + cs%En(i,j,a,fr,m)
499  enddo ; enddo ; enddo
500  call post_data(cs%id_En_mode(fr,m), tot_en, cs%diag)
501  endif ; enddo ; enddo
502 
503  ! Output 3-D (i,j,a) energy density for each freq and mode
504  do m=1,cs%NMode ; do fr=1,cs%Nfreq ; if (cs%id_En_ang_mode(fr,m) > 0) then
505  call post_data(cs%id_En_ang_mode(fr,m), cs%En(:,:,:,fr,m) , cs%diag)
506  endif ; enddo ; enddo
507 
508  ! Output 2-D energy loss (summed over angles, freq, modes)
509  tot_leak_loss(:,:) = 0.0
510  tot_quad_loss(:,:) = 0.0
511  tot_itidal_loss(:,:) = 0.0
512  tot_froude_loss(:,:) = 0.0
513  tot_allprocesses_loss(:,:) = 0.0
514  do m=1,cs%NMode ; do fr=1,cs%Nfreq ; do a=1,cs%nAngle ; do j=js,je ; do i=is,ie
515  tot_leak_loss(i,j) = tot_leak_loss(i,j) + cs%TKE_leak_loss(i,j,a,fr,m)
516  tot_quad_loss(i,j) = tot_quad_loss(i,j) + cs%TKE_quad_loss(i,j,a,fr,m)
517  tot_itidal_loss(i,j) = tot_itidal_loss(i,j) + cs%TKE_itidal_loss(i,j,a,fr,m)
518  tot_froude_loss(i,j) = tot_froude_loss(i,j) + cs%TKE_Froude_loss(i,j,a,fr,m)
519  enddo ; enddo ; enddo ; enddo ; enddo
520  do j=js,je ; do i=is,ie
521  tot_allprocesses_loss(i,j) = tot_leak_loss(i,j) + tot_quad_loss(i,j) + &
522  tot_itidal_loss(i,j) + tot_froude_loss(i,j)
523  enddo ; enddo
524  cs%tot_leak_loss = tot_leak_loss
525  cs%tot_quad_loss = tot_quad_loss
526  cs%tot_itidal_loss = tot_itidal_loss
527  cs%tot_Froude_loss = tot_froude_loss
528  cs%tot_allprocesses_loss = tot_allprocesses_loss
529  if (cs%id_tot_leak_loss > 0) then
530  call post_data(cs%id_tot_leak_loss, tot_leak_loss, cs%diag)
531  endif
532  if (cs%id_tot_quad_loss > 0) then
533  call post_data(cs%id_tot_quad_loss, tot_quad_loss, cs%diag)
534  endif
535  if (cs%id_tot_itidal_loss > 0) then
536  call post_data(cs%id_tot_itidal_loss, tot_itidal_loss, cs%diag)
537  endif
538  if (cs%id_tot_Froude_loss > 0) then
539  call post_data(cs%id_tot_Froude_loss, tot_froude_loss, cs%diag)
540  endif
541  if (cs%id_tot_allprocesses_loss > 0) then
542  call post_data(cs%id_tot_allprocesses_loss, tot_allprocesses_loss, cs%diag)
543  endif
544 
545  ! Output 2-D energy loss (summed over angles) for each freq and mode
546  do m=1,cs%NMode ; do fr=1,cs%Nfreq
547  if (cs%id_itidal_loss_mode(fr,m) > 0 .or. cs%id_allprocesses_loss_mode(fr,m) > 0) then
548  itidal_loss_mode(:,:) = 0.0 ! wave-drag processes (could do others as well)
549  allprocesses_loss_mode(:,:) = 0.0 ! all processes summed together
550  do a=1,cs%nAngle ; do j=js,je ; do i=is,ie
551  itidal_loss_mode(i,j) = itidal_loss_mode(i,j) + cs%TKE_itidal_loss(i,j,a,fr,m)
552  allprocesses_loss_mode(i,j) = allprocesses_loss_mode(i,j) + &
553  cs%TKE_leak_loss(i,j,a,fr,m) + cs%TKE_quad_loss(i,j,a,fr,m) + &
554  cs%TKE_itidal_loss(i,j,a,fr,m) + cs%TKE_Froude_loss(i,j,a,fr,m)
555  enddo ; enddo ; enddo
556  call post_data(cs%id_itidal_loss_mode(fr,m), itidal_loss_mode, cs%diag)
557  call post_data(cs%id_allprocesses_loss_mode(fr,m), allprocesses_loss_mode, cs%diag)
558  endif ; enddo ; enddo
559 
560  ! Output 3-D (i,j,a) energy loss for each freq and mode
561  do m=1,cs%NMode ; do fr=1,cs%Nfreq ; if (cs%id_itidal_loss_ang_mode(fr,m) > 0) then
562  call post_data(cs%id_itidal_loss_ang_mode(fr,m), cs%TKE_itidal_loss(:,:,:,fr,m) , cs%diag)
563  endif ; enddo ; enddo
564 
565  ! Output 2-D period-averaged horizontal near-bottom mode velocity for each freq and mode
566  do m=1,cs%NMode ; do fr=1,cs%Nfreq ; if (cs%id_Ub_mode(fr,m) > 0) then
567  call post_data(cs%id_Ub_mode(fr,m), ub(:,:,fr,m), cs%diag)
568  endif ; enddo ; enddo
569 
570  ! Output 2-D horizontal phase velocity for each freq and mode
571  do m=1,cs%NMode ; do fr=1,cs%Nfreq ; if (cs%id_cp_mode(fr,m) > 0) then
572  call post_data(cs%id_cp_mode(fr,m), cs%cp(:,:,fr,m), cs%diag)
573  endif ; enddo ; enddo
574 
575  endif
576 
577 end subroutine propagate_int_tide
578 
579 !> Checks for energy conservation on computational domain
580 subroutine sum_en(G, CS, En, label)
581  type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure.
582  type(int_tide_cs), pointer :: CS !< The control structure returned by a
583  !! previous call to int_tide_init.
584  real, dimension(G%isd:G%ied,G%jsd:G%jed,CS%NAngle), &
585  intent(in) :: En !< The energy density of the internal tides [J m-2].
586  character(len=*), intent(in) :: label !< A label to use in error messages
587  ! Local variables
588  integer :: m,fr,a
589  real :: En_sum, tmpForSumming, En_sum_diff, En_sum_pdiff
590  character(len=160) :: mesg ! The text of an error message
591  real :: days
592 
593  en_sum = 0.0
594  tmpforsumming = 0.0
595  do a=1,cs%nAngle
596  tmpforsumming = global_area_mean(en(:,:,a),g)*g%areaT_global
597  en_sum = en_sum + tmpforsumming
598  enddo
599  en_sum_diff = en_sum - cs%En_sum
600  if (cs%En_sum /= 0.0) then
601  en_sum_pdiff= (en_sum_diff/cs%En_sum)*100.0
602  else
603  en_sum_pdiff= 0.0
604  endif
605  cs%En_sum = en_sum
606  !! Print to screen
607  !if (is_root_pe()) then
608  ! days = time_type_to_real(CS%Time) / 86400.0
609  ! write(mesg,*) trim(label)//': days =', days, ', En_sum=', En_sum, &
610  ! ', En_sum_diff=', En_sum_diff, ', Percent change=', En_sum_pdiff, '%'
611  ! call MOM_mesg(mesg)
612  !if (is_root_pe() .and. (abs(En_sum_pdiff) > 1.0)) &
613  ! call MOM_error(FATAL, "Run stopped due to excessive internal tide energy change.")
614  !endif
615 
616 end subroutine sum_en
617 
618 !> Calculates the energy lost from the propagating internal tide due to
619 !! scattering over small-scale roughness along the lines of Jayne & St. Laurent (2001).
620 subroutine itidal_lowmode_loss(G, US, CS, Nb, Ub, En, TKE_loss_fixed, TKE_loss, dt, full_halos)
621  type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure.
622  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
623  type(int_tide_cs), pointer :: CS !< The control structure returned by a
624  !! previous call to int_tide_init.
625  real, dimension(G%isd:G%ied,G%jsd:G%jed), &
626  intent(in) :: Nb !< Near-bottom stratification [s-1].
627  real, dimension(G%isd:G%ied,G%jsd:G%jed,CS%nFreq,CS%nMode), &
628  intent(inout) :: Ub !< RMS (over one period) near-bottom horizontal
629  !! mode velocity [m s-1].
630  real, dimension(G%isd:G%ied,G%jsd:G%jed), &
631  intent(in) :: TKE_loss_fixed !< Fixed part of energy loss [kg Z-2 ~> kg m-2]
632  !! (rho*kappa*h^2).
633  real, dimension(G%isd:G%ied,G%jsd:G%jed,CS%NAngle,CS%nFreq,CS%nMode), &
634  intent(inout) :: En !< Energy density of the internal waves [J m-2].
635  real, dimension(G%isd:G%ied,G%jsd:G%jed,CS%NAngle,CS%nFreq,CS%nMode), &
636  intent(out) :: TKE_loss !< Energy loss rate [W m-2]
637  !! (q*rho*kappa*h^2*N*U^2).
638  real, intent(in) :: dt !< Time increment [s].
639  logical,optional, intent(in) :: full_halos !< If true, do the calculation over the
640  !! entirecomputational domain.
641  ! Local variables
642  integer :: j,i,m,fr,a, is, ie, js, je
643  real :: En_tot ! energy for a given mode, frequency, and point summed over angles
644  real :: TKE_loss_tot ! dissipation for a given mode, frequency, and point summed over angles
645  real :: TKE_sum_check ! temporary for check summing
646  real :: frac_per_sector ! fraction of energy in each wedge
647  real :: q_itides ! fraction of energy actually lost to mixing (remainder, 1-q, is
648  ! assumed to stay in propagating mode for now - BDM)
649  real :: loss_rate ! approximate loss rate for implicit calc [s-1]
650  real, parameter :: En_negl = 1e-30 ! negilibly small number to prevent division by zero
651 
652  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
653 
654  q_itides = cs%q_itides
655 
656  if (present(full_halos)) then ; if (full_halos) then
657  is = g%isd ; ie = g%ied ; js = g%jsd ; je = g%jed
658  endif ; endif
659 
660  do j=js,je ; do i=is,ie ; do m=1,cs%nMode ; do fr=1,cs%nFreq
661 
662  ! Sum energy across angles
663  en_tot = 0.0
664  do a=1,cs%nAngle
665  en_tot = en_tot + en(i,j,a,fr,m)
666  enddo
667 
668  ! Calculate TKE loss rate; units of [W m-2] here.
669  tke_loss_tot = q_itides * us%Z_to_m**2 * tke_loss_fixed(i,j) * nb(i,j) * ub(i,j,fr,m)**2
670 
671  ! Update energy remaining (this is a pseudo implicit calc)
672  ! (E(t+1)-E(t))/dt = -TKE_loss(E(t+1)/E(t)), which goes to zero as E(t+1) goes to zero
673  if (en_tot > 0.0) then
674  do a=1,cs%nAngle
675  frac_per_sector = en(i,j,a,fr,m)/en_tot
676  tke_loss(i,j,a,fr,m) = frac_per_sector*tke_loss_tot ! Wm-2
677  loss_rate = tke_loss(i,j,a,fr,m) / (en(i,j,a,fr,m) + en_negl) ! s-1
678  en(i,j,a,fr,m) = en(i,j,a,fr,m) / (1.0 + dt*loss_rate)
679  enddo
680  else
681  ! no loss if no energy
682  tke_loss(i,j,:,fr,m) = 0.0
683  endif
684 
685  ! Update energy remaining (this is the old explicit calc)
686  !if (En_tot > 0.0) then
687  ! do a=1,CS%nAngle
688  ! frac_per_sector = En(i,j,a,fr,m)/En_tot
689  ! TKE_loss(i,j,a,fr,m) = frac_per_sector*TKE_loss_tot
690  ! if (TKE_loss(i,j,a,fr,m)*dt <= En(i,j,a,fr,m))then
691  ! En(i,j,a,fr,m) = En(i,j,a,fr,m) - TKE_loss(i,j,a,fr,m)*dt
692  ! else
693  ! call MOM_error(WARNING, "itidal_lowmode_loss: energy loss greater than avalable, "// &
694  ! " setting En to zero.", all_print=.true.)
695  ! En(i,j,a,fr,m) = 0.0
696  ! endif
697  ! enddo
698  !else
699  ! ! no loss if no energy
700  ! TKE_loss(i,j,:,fr,m) = 0.0
701  !endif
702 
703  enddo ; enddo ; enddo ; enddo
704 
705 end subroutine itidal_lowmode_loss
706 
707 !> This subroutine extracts the energy lost from the propagating internal which has
708 !> been summed across all angles, frequencies, and modes for a given mechanism and location.
709 !!
710 !> It can be called from another module to get values from this module's (private) CS.
711 subroutine get_lowmode_loss(i,j,G,CS,mechanism,TKE_loss_sum)
712  integer, intent(in) :: i !< The i-index of the value to be reported.
713  integer, intent(in) :: j !< The j-index of the value to be reported.
714  type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure
715  type(int_tide_cs), pointer :: cs !< The control structure returned by a
716  !! previous call to int_tide_init.
717  character(len=*), intent(in) :: mechanism !< The named mechanism of loss to return
718  real, intent(out) :: tke_loss_sum !< Total energy loss rate due to specified
719  !! mechanism [W m-2].
720 
721  if (mechanism == 'LeakDrag') tke_loss_sum = cs%tot_leak_loss(i,j) ! not used for mixing yet
722  if (mechanism == 'QuadDrag') tke_loss_sum = cs%tot_quad_loss(i,j) ! not used for mixing yet
723  if (mechanism == 'WaveDrag') tke_loss_sum = cs%tot_itidal_loss(i,j) ! currently used for mixing
724  if (mechanism == 'Froude') tke_loss_sum = cs%tot_Froude_loss(i,j) ! not used for mixing yet
725 
726 end subroutine get_lowmode_loss
727 
728 !> Implements refraction on the internal waves at a single frequency.
729 subroutine refract(En, cn, freq, dt, G, US, NAngle, use_PPMang)
730  type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure.
731  integer, intent(in) :: NAngle !< The number of wave orientations in the
732  !! discretized wave energy spectrum.
733  real, dimension(G%isd:G%ied,G%jsd:G%jed,NAngle), &
734  intent(inout) :: En !< The internal gravity wave energy density as a
735  !! function of space and angular resolution,
736  !! [J m-2 radian-1].
737  real, dimension(G%isd:G%ied,G%jsd:G%jed), &
738  intent(in) :: cn !< Baroclinic mode speed [m s-1].
739  real, intent(in) :: freq !< Wave frequency [s-1].
740  real, intent(in) :: dt !< Time step [s].
741  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
742  logical, intent(in) :: use_PPMang !< If true, use PPM for advection rather
743  !! than upwind.
744  ! Local variables
745  integer, parameter :: stencil = 2
746  real, dimension(SZI_(G),1-stencil:NAngle+stencil) :: &
747  En2d
748  real, dimension(1-stencil:NAngle+stencil) :: &
749  cos_angle, sin_angle
750  real, dimension(SZI_(G)) :: &
751  Dk_Dt_Kmag, Dl_Dt_Kmag
752  real, dimension(SZI_(G),0:nAngle) :: &
753  Flux_E
754  real, dimension(SZI_(G),SZJ_(G),1-stencil:NAngle+stencil) :: &
755  CFL_ang
756  real :: f2 ! The squared Coriolis parameter [s-2].
757  real :: favg ! The average Coriolis parameter at a point [s-1].
758  real :: df2_dy, df2_dx ! The x- and y- gradients of the squared Coriolis parameter [s-2 m-1].
759  real :: df_dy, df_dx ! The x- and y- gradients of the Coriolis parameter [s-1 m-1].
760  real :: dlnCn_dx ! The x-gradient of the wave speed divided by itself [m-1].
761  real :: dlnCn_dy ! The y-gradient of the wave speed divided by itself [m-1].
762  real :: Angle_size, dt_Angle_size, angle
763  real :: Ifreq, Kmag2, I_Kmag
764  real, parameter :: cn_subRO = 1e-100
765  integer :: is, ie, js, je, asd, aed, na
766  integer :: i, j, a
767 
768  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; na = size(en,3)
769  asd = 1-stencil ; aed = nangle+stencil
770 
771  ifreq = 1.0 / freq
772 
773  angle_size = (8.0*atan(1.0)) / (real(nangle))
774  dt_angle_size = dt / angle_size
775 
776  do a=asd,aed
777  angle = (real(a) - 0.5) * angle_size
778  cos_angle(a) = cos(angle) ; sin_angle(a) = sin(angle)
779  enddo
780 
781  !### There should also be refraction due to cn.grad(grid_orientation).
782  cfl_ang(:,:,:) = 0.0
783  do j=js,je
784  ! Copy En into angle space with halos.
785  do a=1,na ; do i=is,ie
786  en2d(i,a) = en(i,j,a)
787  enddo ; enddo
788  do a=asd,0 ; do i=is,ie
789  en2d(i,a) = en2d(i,a+nangle)
790  en2d(i,nangle+stencil+a) = en2d(i,stencil+a)
791  enddo ; enddo
792 
793  ! Do the refraction.
794  do i=is,ie
795  f2 = 0.25*us%s_to_T**2 * ((g%CoriolisBu(i,j)**2 + g%CoriolisBu(i-1,j-1)**2) + &
796  (g%CoriolisBu(i,j-1)**2 + g%CoriolisBu(i-1,j)**2))
797  favg = 0.25*us%s_to_T*((g%CoriolisBu(i,j) + g%CoriolisBu(i-1,j-1)) + &
798  (g%CoriolisBu(i,j-1) + g%CoriolisBu(i-1,j)))
799  df2_dx = 0.5*us%s_to_T**2 * ((g%CoriolisBu(i,j)**2 + g%CoriolisBu(i,j-1)**2) - &
800  (g%CoriolisBu(i-1,j)**2 + g%CoriolisBu(i-1,j-1)**2)) * &
801  g%IdxT(i,j)
802  df_dx = 0.5*us%s_to_T*((g%CoriolisBu(i,j) + g%CoriolisBu(i,j-1)) - &
803  (g%CoriolisBu(i-1,j) + g%CoriolisBu(i-1,j-1))) * &
804  g%IdxT(i,j)
805  dlncn_dx = 0.5*( g%IdxCu(i,j) * (cn(i+1,j) - cn(i,j)) / &
806  (0.5*(cn(i+1,j) + cn(i,j)) + cn_subro) + &
807  g%IdxCu(i-1,j) * (cn(i,j) - cn(i-1,j)) / &
808  (0.5*(cn(i,j) + cn(i-1,j)) + cn_subro) )
809  df2_dy = 0.5*us%s_to_T**2 * ((g%CoriolisBu(i,j)**2 + g%CoriolisBu(i-1,j)**2) - &
810  (g%CoriolisBu(i,j-1)**2 + g%CoriolisBu(i-1,j-1)**2)) * &
811  g%IdyT(i,j)
812  df_dy = 0.5*us%s_to_T*((g%CoriolisBu(i,j) + g%CoriolisBu(i-1,j)) - &
813  (g%CoriolisBu(i,j-1) + g%CoriolisBu(i-1,j-1))) * &
814  g%IdyT(i,j)
815  dlncn_dy = 0.5*( g%IdyCv(i,j) * (cn(i,j+1) - cn(i,j)) / &
816  (0.5*(cn(i,j+1) + cn(i,j)) + cn_subro) + &
817  g%IdyCv(i,j-1) * (cn(i,j) - cn(i,j-1)) / &
818  (0.5*(cn(i,j) + cn(i,j-1)) + cn_subro) )
819  kmag2 = (freq**2 - f2) / (cn(i,j)**2 + cn_subro**2)
820  if (kmag2 > 0.0) then
821  i_kmag = 1.0 / sqrt(kmag2)
822  dk_dt_kmag(i) = -ifreq * (favg*df_dx + (freq**2 - f2) * dlncn_dx) * i_kmag
823  dl_dt_kmag(i) = -ifreq * (favg*df_dy + (freq**2 - f2) * dlncn_dy) * i_kmag
824  else
825  dk_dt_kmag(i) = 0.0
826  dl_dt_kmag(i) = 0.0
827  endif
828  enddo
829 
830  ! Determine the energy fluxes in angular orientation space.
831  do a=asd,aed ; do i=is,ie
832  cfl_ang(i,j,a) = (cos_angle(a) * dl_dt_kmag(i) - sin_angle(a) * dk_dt_kmag(i)) * &
833  dt_angle_size
834  if (abs(cfl_ang(i,j,a)) > 1.0) then
835  call mom_error(warning, "refract: CFL exceeds 1.", .true.)
836  if (cfl_ang(i,j,a) > 0.0) then ; cfl_ang(i,j,a) = 1.0 ; else ; cfl_ang(i,j,a) = -1.0 ; endif
837  endif
838  enddo ; enddo
839 
840  ! Advect in angular space
841  if (.not.use_ppmang) then
842  ! Use simple upwind
843  do a=0,na ; do i=is,ie
844  if (cfl_ang(i,j,a) > 0.0) then
845  flux_e(i,a) = cfl_ang(i,j,a) * en2d(i,a)
846  else
847  flux_e(i,a) = cfl_ang(i,j,a) * en2d(i,a+1)
848  endif
849  enddo ; enddo
850  else
851  ! Use PPM
852  do i=is,ie
853  call ppm_angular_advect(en2d(i,:),cfl_ang(i,j,:),flux_e(i,:),nangle,dt,stencil)
854  enddo
855  endif
856 
857  ! Update and copy back to En.
858  do a=1,na ; do i=is,ie
859  !if (En2d(i,a)+(Flux_E(i,A-1)-Flux_E(i,A)) < 0.0) then ! for debugging
860  ! call MOM_error(FATAL, "refract: OutFlux>Available")
861  !endif
862  en(i,j,a) = en2d(i,a) + (flux_e(i,a-1) - flux_e(i,a))
863  enddo ; enddo
864  enddo ! j-loop
865 end subroutine refract
866 
867 !> This subroutine calculates the 1-d flux for advection in angular space using a monotonic
868 !! piecewise parabolic scheme. This needs to be called from within i and j spatial loops.
869 subroutine ppm_angular_advect(En2d, CFL_ang, Flux_En, NAngle, dt, halo_ang)
870  integer, intent(in) :: NAngle !< The number of wave orientations in the
871  !! discretized wave energy spectrum.
872  real, intent(in) :: dt !< Time increment [s].
873  integer, intent(in) :: halo_ang !< The halo size in angular space
874  real, dimension(1-halo_ang:NAngle+halo_ang), &
875  intent(in) :: En2d !< The internal gravity wave energy density as a
876  !! function of angular resolution [J m-2 radian-1].
877  real, dimension(1-halo_ang:NAngle+halo_ang), &
878  intent(in) :: CFL_ang !< The CFL number of the energy advection across angles
879  real, dimension(0:NAngle), intent(out) :: Flux_En !< The time integrated internal wave energy flux
880  !! across angles [J m-2 radian-1].
881  ! Local variables
882  real :: flux
883  real :: u_ang
884  real :: Angle_size
885  real :: I_Angle_size
886  real :: I_dt
887  integer :: a
888  real :: aR, aL, dMx, dMn, Ep, Ec, Em, dA, mA, a6
889 
890  i_dt = 1 / dt
891  angle_size = (8.0*atan(1.0)) / (real(nangle))
892  i_angle_size = 1 / angle_size
893  flux_en(:) = 0
894 
895  do a=0,nangle
896  u_ang = cfl_ang(a)*angle_size*i_dt
897  if (u_ang >= 0.0) then
898  ! Implementation of PPM-H3
899  ep = en2d(a+1)*i_angle_size !MEAN ANGULAR ENERGY DENSITY FOR WEDGE (Jm-2/rad)
900  ec = en2d(a) *i_angle_size !MEAN ANGULAR ENERGY DENSITY FOR WEDGE (Jm-2/rad)
901  em = en2d(a-1)*i_angle_size !MEAN ANGULAR ENERGY DENSITY FOR WEDGE (Jm-2/rad)
902  al = ( 5.*ec + ( 2.*em - ep ) )/6. ! H3 estimate
903  al = max( min(ec,em), al) ; al = min( max(ec,em), al) ! Bound
904  ar = ( 5.*ec + ( 2.*ep - em ) )/6. ! H3 estimate
905  ar = max( min(ec,ep), ar) ; ar = min( max(ec,ep), ar) ! Bound
906  da = ar - al ; ma = 0.5*( ar + al )
907  if ((ep-ec)*(ec-em) <= 0.) then
908  al = ec ; ar = ec ! PCM for local extremum
909  elseif ( da*(ec-ma) > (da*da)/6. ) then
910  al = 3.*ec - 2.*ar !?
911  elseif ( da*(ec-ma) < - (da*da)/6. ) then
912  ar = 3.*ec - 2.*al !?
913  endif
914  a6 = 6.*ec - 3. * (ar + al) ! Curvature
915  ! CALCULATE FLUX RATE (Jm-2/s)
916  flux = u_ang*( ar + 0.5 * cfl_ang(a) * ( ( al - ar ) + a6 * ( 1. - 2./3. * cfl_ang(a) ) ) )
917  !flux = u_ang*( aR - 0.5 * CFL_ang(A) * ( ( aR - aL ) - a6 * ( 1. - 2./3. * CFL_ang(A) ) ) )
918  ! CALCULATE AMOUNT FLUXED (Jm-2)
919  flux_en(a) = dt * flux
920  !Flux_En(A) = (dt * I_Angle_size) * flux
921  else
922  ! Implementation of PPM-H3
923  ep = en2d(a+2)*i_angle_size !MEAN ANGULAR ENERGY DENSITY FOR WEDGE (Jm-2/rad)
924  ec = en2d(a+1)*i_angle_size !MEAN ANGULAR ENERGY DENSITY FOR WEDGE (Jm-2/rad)
925  em = en2d(a) *i_angle_size !MEAN ANGULAR ENERGY DENSITY FOR WEDGE (Jm-2/rad)
926  al = ( 5.*ec + ( 2.*em - ep ) )/6. ! H3 estimate
927  al = max( min(ec,em), al) ; al = min( max(ec,em), al) ! Bound
928  ar = ( 5.*ec + ( 2.*ep - em ) )/6. ! H3 estimate
929  ar = max( min(ec,ep), ar) ; ar = min( max(ec,ep), ar) ! Bound
930  da = ar - al ; ma = 0.5*( ar + al )
931  if ((ep-ec)*(ec-em) <= 0.) then
932  al = ec ; ar = ec ! PCM for local extremum
933  elseif ( da*(ec-ma) > (da*da)/6. ) then
934  al = 3.*ec - 2.*ar
935  elseif ( da*(ec-ma) < - (da*da)/6. ) then
936  ar = 3.*ec - 2.*al
937  endif
938  a6 = 6.*ec - 3. * (ar + al) ! Curvature
939  ! CALCULATE FLUX RATE (Jm-2/s)
940  flux = u_ang*( ar + 0.5 * cfl_ang(a) * ( ( al - ar ) + a6 * ( 1. - 2./3. * cfl_ang(a) ) ) )
941  !flux = u_ang*( aL + 0.5 * CFL_ang(A) * ( ( aR - aL ) + a6 * ( 1. - 2./3. * CFL_ang(A) ) ) )
942  ! CALCULATE AMOUNT FLUXED (Jm-2)
943  flux_en(a) = dt * flux
944  !Flux_En(A) = (dt * I_Angle_size) * flux
945  endif
946  enddo
947 end subroutine ppm_angular_advect
948 
949 !> Propagates internal waves at a single frequency.
950 subroutine propagate(En, cn, freq, dt, G, US, CS, NAngle)
951  type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure.
952  integer, intent(in) :: NAngle !< The number of wave orientations in the
953  !! discretized wave energy spectrum.
954  real, dimension(G%isd:G%ied,G%jsd:G%jed,NAngle), &
955  intent(inout) :: En !< The internal gravity wave energy density as a
956  !! function of space and angular resolution,
957  !! [J m-2 radian-1].
958  real, dimension(G%isd:G%ied,G%jsd:G%jed), &
959  intent(in) :: cn !< Baroclinic mode speed [m s-1].
960  real, intent(in) :: freq !< Wave frequency [s-1].
961  real, intent(in) :: dt !< Time step [s].
962  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
963  type(int_tide_cs), pointer :: CS !< The control structure returned by a
964  !! previous call to int_tide_init.
965  ! Local variables
966  real, dimension(G%IsdB:G%IedB,G%JsdB:G%JedB) :: &
967  speed ! The magnitude of the group velocity at the q points for corner adv [m s-1].
968  integer, parameter :: stencil = 2
969  real, dimension(SZIB_(G),SZJ_(G)) :: &
970  speed_x ! The magnitude of the group velocity at the Cu points [m s-1].
971  real, dimension(SZI_(G),SZJB_(G)) :: &
972  speed_y ! The magnitude of the group velocity at the Cv points [m s-1].
973  real, dimension(0:NAngle) :: &
974  cos_angle, sin_angle
975  real, dimension(NAngle) :: &
976  Cgx_av, Cgy_av, dCgx, dCgy
977  real :: f2 ! The squared Coriolis parameter [s-2].
978  real :: Angle_size, I_Angle_size, angle
979  real :: Ifreq, freq2
980  real, parameter :: cn_subRO = 1e-100
981  type(loop_bounds_type) :: LB
982  integer :: is, ie, js, je, asd, aed, na
983  integer :: ish, ieh, jsh, jeh
984  integer :: i, j, a
985 
986  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; na = size(en,3)
987  asd = 1-stencil ; aed = nangle+stencil
988 
989  ifreq = 1.0 / freq
990  freq2 = freq**2
991 
992  ! Define loop bounds: Need extensions on j-loop so propagate_y
993  ! (done after propagate_x) will have updated values in the halo
994  ! for correct PPM reconstruction. Use if no teleporting and
995  ! no pass_var between propagate_x and propagate_y.
996  !jsh = js-3 ; jeh = je+3 ; ish = is ; ieh = ie
997 
998  ! Define loop bounds: Need 1-pt extensions on loops because
999  ! teleporting eats up a halo point. Use if teleporting.
1000  ! Also requires pass_var before propagate_y.
1001  jsh = js-1 ; jeh = je+1 ; ish = is-1 ; ieh = ie+1
1002 
1003  angle_size = (8.0*atan(1.0)) / real(nangle)
1004  i_angle_size = 1.0 / angle_size
1005 
1006  if (cs%corner_adv) then
1007  ! IMPLEMENT CORNER ADVECTION IN HORIZONTAL--------------------
1008  ! FIND AVERAGE GROUP VELOCITY (SPEED) AT CELL CORNERS
1009  ! NOTE: THIS HAS NOT BE ADAPTED FOR REFLECTION YET (BDM)!!
1010  ! Fix indexing here later
1011  speed(:,:) = 0
1012  do j=jsh-1,jeh ; do i=ish-1,ieh
1013  f2 = us%s_to_T**2 * g%CoriolisBu(i,j)**2
1014  speed(i,j) = 0.25*(cn(i,j) + cn(i+1,j) + cn(i+1,j+1) + cn(i,j+1)) * &
1015  sqrt(max(freq2 - f2, 0.0)) * ifreq
1016  enddo ; enddo
1017  do a=1,na
1018  ! Apply the propagation WITH CORNER ADVECTION/FINITE VOLUME APPROACH.
1019  lb%jsh = js ; lb%jeh = je ; lb%ish = is ; lb%ieh = ie
1020  call propagate_corner_spread(en(:,:,a), a, nangle, speed, dt, g, cs, lb)
1021  enddo ! a-loop
1022  else
1023  ! IMPLEMENT PPM ADVECTION IN HORIZONTAL-----------------------
1024  ! These could be in the control structure, as they do not vary.
1025  do a=0,na
1026  ! These are the angles at the cell edges...
1027  angle = (real(a) - 0.5) * angle_size
1028  cos_angle(a) = cos(angle) ; sin_angle(a) = sin(angle)
1029  enddo
1030 
1031  do a=1,na
1032  cgx_av(a) = (sin_angle(a) - sin_angle(a-1)) * i_angle_size
1033  cgy_av(a) = -(cos_angle(a) - cos_angle(a-1)) * i_angle_size
1034  dcgx(a) = sqrt(0.5 + 0.5*(sin_angle(a)*cos_angle(a) - &
1035  sin_angle(a-1)*cos_angle(a-1)) * i_angle_size - &
1036  cgx_av(a)**2)
1037  dcgy(a) = sqrt(0.5 - 0.5*(sin_angle(a)*cos_angle(a) - &
1038  sin_angle(a-1)*cos_angle(a-1)) * i_angle_size - &
1039  cgy_av(a)**2)
1040  enddo
1041 
1042  do j=jsh,jeh ; do i=ish-1,ieh
1043  f2 = 0.5*us%s_to_T**2 * (g%CoriolisBu(i,j)**2 + g%CoriolisBu(i,j-1)**2)
1044  speed_x(i,j) = 0.5*(cn(i,j) + cn(i+1,j)) * g%mask2dCu(i,j) * &
1045  sqrt(max(freq2 - f2, 0.0)) * ifreq
1046  enddo ; enddo
1047  do j=jsh-1,jeh ; do i=ish,ieh
1048  f2 = 0.5*us%s_to_T**2 * (g%CoriolisBu(i,j)**2 + g%CoriolisBu(i-1,j)**2)
1049  speed_y(i,j) = 0.5*(cn(i,j) + cn(i,j+1)) * g%mask2dCv(i,j) * &
1050  sqrt(max(freq2 - f2, 0.0)) * ifreq
1051  enddo ; enddo
1052 
1053  ! Apply propagation in x-direction (reflection included)
1054  lb%jsh = jsh ; lb%jeh = jeh ; lb%ish = ish ; lb%ieh = ieh
1055  call propagate_x(en(:,:,:), speed_x, cgx_av(:), dcgx(:), dt, g, cs%nAngle, cs, lb)
1056 
1057  ! Check for energy conservation on computational domain (for debugging)
1058  !call sum_En(G,CS,En(:,:,:),'post-propagate_x')
1059 
1060  ! Update halos
1061  call pass_var(en(:,:,:),g%domain)
1062 
1063  ! Apply propagation in y-direction (reflection included)
1064  ! LB%jsh = js ; LB%jeh = je ; LB%ish = is ; LB%ieh = ie ! Use if no teleport
1065  lb%jsh = jsh ; lb%jeh = jeh ; lb%ish = ish ; lb%ieh = ieh
1066  call propagate_y(en(:,:,:), speed_y, cgy_av(:), dcgy(:), dt, g, cs%nAngle, cs, lb)
1067 
1068  ! Check for energy conservation on computational domain (for debugging)
1069  !call sum_En(G,CS,En(:,:,:),'post-propagate_y')
1070 
1071  endif
1072 end subroutine propagate
1073 
1074 !> This subroutine does first-order corner advection. It was written with the hopes
1075 !! of smoothing out the garden sprinkler effect, but is too numerically diffusive to
1076 !! be of much use as of yet. It is not yet compatible with reflection schemes (BDM).
1077 subroutine propagate_corner_spread(En, energized_wedge, NAngle, speed, dt, G, CS, LB)
1078  type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure.
1079  real, dimension(G%isd:G%ied,G%jsd:G%jed), &
1080  intent(inout) :: En !< The energy density integrated over an angular
1081  !! band [W m-2], intent in/out.
1082  real, dimension(G%IsdB:G%IedB,G%Jsd:G%Jed), &
1083  intent(in) :: speed !< The magnitude of the group velocity at the cell
1084  !! corner points [m s-1].
1085  integer, intent(in) :: energized_wedge !< Index of current ray direction.
1086  integer, intent(in) :: NAngle !< The number of wave orientations in the
1087  !! discretized wave energy spectrum.
1088  real, intent(in) :: dt !< Time increment [s].
1089  type(int_tide_cs), pointer :: CS !< The control structure returned by a previous
1090  !! call to continuity_PPM_init.
1091  type(loop_bounds_type), intent(in) :: LB !< A structure with the active energy loop bounds.
1092  ! Local variables
1093  integer :: i, j, k, ish, ieh, jsh, jeh, m
1094  real :: TwoPi, Angle_size
1095  real :: energized_angle ! angle through center of current wedge
1096  real :: theta ! angle at edge of wedge
1097  real :: Nsubrays ! number of sub-rays for averaging
1098  ! count includes the two rays that bound the current wedge,
1099  ! i.e. those at -dtheta/2 and +dtheta/2 from energized angle
1100  real :: I_Nsubwedges ! inverse of number of sub-wedges
1101  real :: cos_thetaDT, sin_thetaDT ! cos(theta)*dt, sin(theta)*dt
1102  real :: xNE,xNW,xSW,xSE,yNE,yNW,ySW,ySE ! corner point coordinates of advected fluid parcel
1103  real :: CFL_xNE,CFL_xNW,CFL_xSW,CFL_xSE,CFL_yNE,CFL_yNW,CFL_ySW,CFL_ySE,CFL_max
1104  real :: xN,xS,xE,xW,yN,yS,yE,yW ! intersection point coordinates of parcel edges and grid
1105  real :: xCrn,yCrn ! grid point contained within advected fluid parcel
1106  real :: xg,yg ! grid point of interest
1107  real :: slopeN,slopeW,slopeS,slopeE, bN,bW,bS,bE ! parameters defining parcel sides
1108  real :: aNE,aN,aNW,aW,aSW,aS,aSE,aE,aC ! sub-areas of advected parcel
1109  real :: a_total ! total area of advected parcel
1110  real :: a1,a2,a3,a4 ! areas used in calculating polygon areas (sub-areas) of advected parcel
1111  real, dimension(G%IsdB:G%IedB,G%Jsd:G%Jed) :: x,y ! coordinates of cell corners
1112  real, dimension(G%IsdB:G%IedB,G%Jsd:G%Jed) :: Idx,Idy ! inverse of dx,dy at cell corners
1113  real, dimension(G%IsdB:G%IedB,G%Jsd:G%Jed) :: dx,dy ! dx,dy at cell corners
1114  real, dimension(2) :: E_new ! energy in cell after advection for subray; set size here to
1115  ! define Nsubrays - this should be made an input option later!
1116 
1117  ish = lb%ish ; ieh = lb%ieh ; jsh = lb%jsh ; jeh = lb%jeh
1118  twopi = (8.0*atan(1.0))
1119  nsubrays = real(size(e_new))
1120  i_nsubwedges = 1./(nsubrays - 1)
1121 
1122  angle_size = twopi / real(nangle)
1123  energized_angle = angle_size * real(energized_wedge - 1) ! for a=1 aligned with x-axis
1124  !energized_angle = Angle_size * real(energized_wedge - 1) + 2.0*Angle_size !
1125  !energized_angle = Angle_size * real(energized_wedge - 1) + 0.5*Angle_size !
1126  x = g%geoLonBu
1127  y = g%geoLatBu
1128  idx = g%IdxBu; dx = g%dxBu
1129  idy = g%IdyBu; dy = g%dyBu
1130 
1131  do j=jsh,jeh; do i=ish,ieh
1132  do m=1,int(nsubrays)
1133  theta = energized_angle - 0.5*angle_size + real(m - 1)*angle_size*i_nsubwedges
1134  if (theta < 0.0) then
1135  theta = theta + twopi
1136  elseif (theta > twopi) then
1137  theta = theta - twopi
1138  endif
1139  cos_thetadt = cos(theta)*dt
1140  sin_thetadt = sin(theta)*dt
1141 
1142  ! corner point coordinates of advected fluid parcel ----------
1143  xg = x(i,j); yg = y(i,j)
1144  xne = xg - speed(i,j)*cos_thetadt
1145  yne = yg - speed(i,j)*sin_thetadt
1146  cfl_xne = (xg-xne)*idx(i,j)
1147  cfl_yne = (yg-yne)*idy(i,j)
1148 
1149  xg = x(i-1,j); yg = y(i-1,j)
1150  xnw = xg - speed(i-1,j)*cos_thetadt
1151  ynw = yg - speed(i-1,j)*sin_thetadt
1152  cfl_xnw = (xg-xnw)*idx(i-1,j)
1153  cfl_ynw = (yg-ynw)*idy(i-1,j)
1154 
1155  xg = x(i-1,j-1); yg = y(i-1,j-1)
1156  xsw = xg - speed(i-1,j-1)*cos_thetadt
1157  ysw = yg - speed(i-1,j-1)*sin_thetadt
1158  cfl_xsw = (xg-xsw)*idx(i-1,j-1)
1159  cfl_ysw = (yg-ysw)*idy(i-1,j-1)
1160 
1161  xg = x(i,j-1); yg = y(i,j-1)
1162  xse = xg - speed(i,j-1)*cos_thetadt
1163  yse = yg - speed(i,j-1)*sin_thetadt
1164  cfl_xse = (xg-xse)*idx(i,j-1)
1165  cfl_yse = (yg-yse)*idy(i,j-1)
1166 
1167  cfl_max = max(abs(cfl_xne),abs(cfl_xnw),abs(cfl_xsw), &
1168  abs(cfl_xse),abs(cfl_yne),abs(cfl_ynw), &
1169  abs(cfl_ysw),abs(cfl_yse))
1170  if (cfl_max > 1.0) then
1171  call mom_error(warning, "propagate_corner_spread: CFL exceeds 1.", .true.)
1172  endif
1173 
1174  ! intersection point coordinates of parcel edges and cell edges ---
1175  if (0.0 <= theta .and. theta < 0.25*twopi) then
1176  xn = x(i-1,j-1)
1177  yw = y(i-1,j-1)
1178  elseif (0.25*twopi <= theta .and. theta < 0.5*twopi) then
1179  xn = x(i,j-1)
1180  yw = y(i,j-1)
1181  elseif (0.5*twopi <= theta .and. theta < 0.75*twopi) then
1182  xn = x(i,j)
1183  yw = y(i,j)
1184  elseif (0.75*twopi <= theta .and. theta <= 1.00*twopi) then
1185  xn = x(i-1,j)
1186  yw = y(i-1,j)
1187  endif
1188  xs = xn
1189  ye = yw
1190 
1191  ! north intersection
1192  slopen = (yne - ynw)/(xne - xnw)
1193  bn = -slopen*xne + yne
1194  yn = slopen*xn + bn
1195  ! west intersection
1196  if (xnw == xsw) then
1197  xw = xnw
1198  else
1199  slopew = (ynw - ysw)/(xnw - xsw)
1200  bw = -slopew*xnw + ynw
1201  xw = (yw - bw)/slopew
1202  endif
1203  ! south intersection
1204  slopes = (ysw - yse)/(xsw - xse)
1205  bs = -slopes*xsw + ysw
1206  ys = slopes*xs + bs
1207  ! east intersection
1208  if (xne == xse) then
1209  xe = xne
1210  else
1211  slopee = (yse - yne)/(xse - xne)
1212  be = -slopee*xse + yse
1213  xe = (ye - be)/slopee
1214  endif
1215 
1216  ! areas --------------------------------------------
1217  ane = 0.0; an = 0.0; anw = 0.0; ! initialize areas
1218  aw = 0.0; asw = 0.0; as = 0.0; ! initialize areas
1219  ase = 0.0; ae = 0.0; ac = 0.0; ! initialize areas
1220  if (0.0 <= theta .and. theta < 0.25*twopi) then
1221  xcrn = x(i-1,j-1); ycrn = y(i-1,j-1)
1222  ! west area
1223  a1 = (yn - ycrn)*(0.5*(xn + xcrn))
1224  a2 = (ycrn - yw)*(0.5*(xcrn + xw))
1225  a3 = (yw - ynw)*(0.5*(xw + xnw))
1226  a4 = (ynw - yn)*(0.5*(xnw + xn))
1227  aw = a1 + a2 + a3 + a4
1228  ! southwest area
1229  a1 = (ycrn - ys)*(0.5*(xcrn + xs))
1230  a2 = (ys - ysw)*(0.5*(xs + xsw))
1231  a3 = (ysw - yw)*(0.5*(xsw + xw))
1232  a4 = (yw - ycrn)*(0.5*(xw + xcrn))
1233  asw = a1 + a2 + a3 + a4
1234  ! south area
1235  a1 = (ye - yse)*(0.5*(xe + xse))
1236  a2 = (yse - ys)*(0.5*(xse + xs))
1237  a3 = (ys - ycrn)*(0.5*(xs + xcrn))
1238  a4 = (ycrn - ye)*(0.5*(xcrn + xe))
1239  as = a1 + a2 + a3 + a4
1240  ! area within cell
1241  a1 = (yne - ye)*(0.5*(xne + xe))
1242  a2 = (ye - ycrn)*(0.5*(xe + xcrn))
1243  a3 = (ycrn - yn)*(0.5*(xcrn + xn))
1244  a4 = (yn - yne)*(0.5*(xn + xne))
1245  ac = a1 + a2 + a3 + a4
1246  elseif (0.25*twopi <= theta .and. theta < 0.5*twopi) then
1247  xcrn = x(i,j-1); ycrn = y(i,j-1)
1248  ! south area
1249  a1 = (ycrn - ys)*(0.5*(xcrn + xs))
1250  a2 = (ys - ysw)*(0.5*(xs + xsw))
1251  a3 = (ysw - yw)*(0.5*(xsw + xw))
1252  a4 = (yw - ycrn)*(0.5*(xw + xcrn))
1253  as = a1 + a2 + a3 + a4
1254  ! southeast area
1255  a1 = (ye - yse)*(0.5*(xe + xse))
1256  a2 = (yse - ys)*(0.5*(xse + xs))
1257  a3 = (ys - ycrn)*(0.5*(xs + xcrn))
1258  a4 = (ycrn - ye)*(0.5*(xcrn + xe))
1259  ase = a1 + a2 + a3 + a4
1260  ! east area
1261  a1 = (yne - ye)*(0.5*(xne + xe))
1262  a2 = (ye - ycrn)*(0.5*(xe + xcrn))
1263  a3 = (ycrn - yn)*(0.5*(xcrn + xn))
1264  a4 = (yn - yne)*(0.5*(xn + xne))
1265  ae = a1 + a2 + a3 + a4
1266  ! area within cell
1267  a1 = (yn - ycrn)*(0.5*(xn + xcrn))
1268  a2 = (ycrn - yw)*(0.5*(xcrn + xw))
1269  a3 = (yw - ynw)*(0.5*(xw + xnw))
1270  a4 = (ynw - yn)*(0.5*(xnw + xn))
1271  ac = a1 + a2 + a3 + a4
1272  elseif (0.5*twopi <= theta .and. theta < 0.75*twopi) then
1273  xcrn = x(i,j); ycrn = y(i,j)
1274  ! east area
1275  a1 = (ye - yse)*(0.5*(xe + xse))
1276  a2 = (yse - ys)*(0.5*(xse + xs))
1277  a3 = (ys - ycrn)*(0.5*(xs + xcrn))
1278  a4 = (ycrn - ye)*(0.5*(xcrn + xe))
1279  ae = a1 + a2 + a3 + a4
1280  ! northeast area
1281  a1 = (yne - ye)*(0.5*(xne + xe))
1282  a2 = (ye - ycrn)*(0.5*(xe + xcrn))
1283  a3 = (ycrn - yn)*(0.5*(xcrn + xn))
1284  a4 = (yn - yne)*(0.5*(xn + xne))
1285  ane = a1 + a2 + a3 + a4
1286  ! north area
1287  a1 = (yn - ycrn)*(0.5*(xn + xcrn))
1288  a2 = (ycrn - yw)*(0.5*(xcrn + xw))
1289  a3 = (yw - ynw)*(0.5*(xw + xnw))
1290  a4 = (ynw - yn)*(0.5*(xnw + xn))
1291  an = a1 + a2 + a3 + a4
1292  ! area within cell
1293  a1 = (ycrn - ys)*(0.5*(xcrn + xs))
1294  a2 = (ys - ysw)*(0.5*(xs + xsw))
1295  a3 = (ysw - yw)*(0.5*(xsw + xw))
1296  a4 = (yw - ycrn)*(0.5*(xw + xcrn))
1297  ac = a1 + a2 + a3 + a4
1298  elseif (0.75*twopi <= theta .and. theta <= 1.00*twopi) then
1299  xcrn = x(i-1,j); ycrn = y(i-1,j)
1300  ! north area
1301  a1 = (yne - ye)*(0.5*(xne + xe))
1302  a2 = (ye - ycrn)*(0.5*(xe + xcrn))
1303  a3 = (ycrn - yn)*(0.5*(xcrn + xn))
1304  a4 = (yn - yne)*(0.5*(xn + xne))
1305  an = a1 + a2 + a3 + a4
1306  ! northwest area
1307  a1 = (yn - ycrn)*(0.5*(xn + xcrn))
1308  a2 = (ycrn - yw)*(0.5*(xcrn + xw))
1309  a3 = (yw - ynw)*(0.5*(xw + xnw))
1310  a4 = (ynw - yn)*(0.5*(xnw + xn))
1311  anw = a1 + a2 + a3 + a4
1312  ! west area
1313  a1 = (ycrn - ys)*(0.5*(xcrn + xs))
1314  a2 = (ys - ysw)*(0.5*(xs + xsw))
1315  a3 = (ysw - yw)*(0.5*(xsw + xw))
1316  a4 = (yw - ycrn)*(0.5*(xw + xcrn))
1317  aw = a1 + a2 + a3 + a4
1318  ! area within cell
1319  a1 = (ye - yse)*(0.5*(xe + xse))
1320  a2 = (yse - ys)*(0.5*(xse + xs))
1321  a3 = (ys - ycrn)*(0.5*(xs + xcrn))
1322  a4 = (ycrn - ye)*(0.5*(xcrn + xe))
1323  ac = a1 + a2 + a3 + a4
1324  endif
1325 
1326  ! energy weighting ----------------------------------------
1327  a_total = ane + an + anw + aw + asw + as + ase + ae + ac
1328  e_new(m) = (ane*en(i+1,j+1) + an*en(i,j+1) + anw*en(i-1,j+1) + &
1329  aw*en(i-1,j) + asw*en(i-1,j-1) + as*en(i,j-1) + &
1330  ase*en(i+1,j-1) + ae*en(i+1,j) + ac*en(i,j)) / (dx(i,j)*dy(i,j))
1331  enddo ! m-loop
1332  ! update energy in cell
1333  en(i,j) = sum(e_new)/nsubrays
1334  enddo ; enddo
1335 end subroutine propagate_corner_spread
1336 
1337 !> Propagates the internal wave energy in the logical x-direction.
1338 subroutine propagate_x(En, speed_x, Cgx_av, dCgx, dt, G, Nangle, CS, LB)
1339  type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure.
1340  integer, intent(in) :: NAngle !< The number of wave orientations in the
1341  !! discretized wave energy spectrum.
1342  real, dimension(G%isd:G%ied,G%jsd:G%jed,Nangle), &
1343  intent(inout) :: En !< The energy density integrated over an angular
1344  !! band [J m-2], intent in/out.
1345  real, dimension(G%IsdB:G%IedB,G%jsd:G%jed), &
1346  intent(in) :: speed_x !< The magnitude of the group velocity at the
1347  !! Cu points [m s-1].
1348  real, dimension(Nangle), intent(in) :: Cgx_av !< The average x-projection in each angular band.
1349  real, dimension(Nangle), intent(in) :: dCgx !< The difference in x-projections between the
1350  !! edges of each angular band.
1351  real, intent(in) :: dt !< Time increment [s].
1352  type(int_tide_cs), pointer :: CS !< The control structure returned by a previous call
1353  !! to continuity_PPM_init.
1354  type(loop_bounds_type), intent(in) :: LB !< A structure with the active energy loop bounds.
1355  ! Local variables
1356  real, dimension(SZI_(G),SZJ_(G)) :: &
1357  EnL, EnR ! Left and right face energy densities [J m-2].
1358  real, dimension(SZIB_(G),SZJ_(G)) :: &
1359  flux_x ! The internal wave energy flux [J s-1].
1360  real, dimension(SZIB_(G)) :: &
1361  cg_p, cg_m, flux1, flux2
1362  !real, dimension(SZI_(G),SZJB_(G),Nangle) :: En_m, En_p
1363  real, dimension(SZI_(G),SZJB_(G),Nangle) :: &
1364  Fdt_m, Fdt_p! Left and right energy fluxes [J]
1365  integer :: i, j, k, ish, ieh, jsh, jeh, a
1366 
1367  ish = lb%ish ; ieh = lb%ieh ; jsh = lb%jsh ; jeh = lb%jeh
1368  do a=1,nangle
1369  ! This sets EnL and EnR.
1370  if (cs%upwind_1st) then
1371  do j=jsh,jeh ; do i=ish-1,ieh+1
1372  enl(i,j) = en(i,j,a) ; enr(i,j) = en(i,j,a)
1373  enddo ; enddo
1374  else
1375  call ppm_reconstruction_x(en(:,:,a), enl, enr, g, lb, simple_2nd=cs%simple_2nd)
1376  endif
1377 
1378  do j=jsh,jeh
1379  ! This is done once with single speed (GARDEN SPRINKLER EFFECT POSSIBLE)
1380  do i=ish-1,ieh
1381  cg_p(i) = speed_x(i,j) * (cgx_av(a))
1382  enddo
1383  call zonal_flux_en(cg_p, en(:,j,a), enl(:,j), enr(:,j), flux1, &
1384  dt, g, j, ish, ieh, cs%vol_CFL)
1385  do i=ish-1,ieh ; flux_x(i,j) = flux1(i); enddo
1386  enddo
1387 
1388  do j=jsh,jeh ; do i=ish,ieh
1389  fdt_m(i,j,a) = dt*flux_x(i-1,j) ! left face influx (J)
1390  fdt_p(i,j,a) = -dt*flux_x(i,j) ! right face influx (J)
1391  enddo ; enddo
1392 
1393  ! test with old (take out later)
1394  !do j=LB%jsh,LB%jeh ; do i=LB%ish,LB%ieh
1395  ! En(i,j,a) = En(i,j,a) - dt* G%IareaT(i,j) * (flux_x(I,j) - flux_x(I-1,j))
1396  !enddo ; enddo
1397 
1398  enddo ! a-loop
1399 
1400  ! Only reflect newly arrived energy; existing energy in incident wedge
1401  ! is not reflected and will eventually propagate out of cell.
1402  ! (only reflects if En > 0)
1403  call reflect(fdt_m(:,:,:), nangle, cs, g, lb)
1404  call teleport(fdt_m(:,:,:), nangle, cs, g, lb)
1405  call reflect(fdt_p(:,:,:), nangle, cs, g, lb)
1406  call teleport(fdt_p(:,:,:), nangle, cs, g, lb)
1407 
1408  ! Update reflected energy (Jm-2)
1409  do j=jsh,jeh ; do i=ish,ieh
1410  !do a=1,CS%nAngle
1411  ! if ((En(i,j,a) + G%IareaT(i,j)*(Fdt_m(i,j,a) + Fdt_p(i,j,a))) < 0.0) then ! for debugging
1412  ! call MOM_error(FATAL, "propagate_x: OutFlux>Available")
1413  ! endif
1414  !enddo
1415  en(i,j,:) = en(i,j,:) + g%IareaT(i,j)*(fdt_m(i,j,:) + fdt_p(i,j,:))
1416  enddo ; enddo
1417 
1418 end subroutine propagate_x
1419 
1420 !> Propagates the internal wave energy in the logical y-direction.
1421 subroutine propagate_y(En, speed_y, Cgy_av, dCgy, dt, G, Nangle, CS, LB)
1422  type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure.
1423  integer, intent(in) :: NAngle !< The number of wave orientations in the
1424  !! discretized wave energy spectrum.
1425  real, dimension(G%isd:G%ied,G%jsd:G%jed,Nangle), &
1426  intent(inout) :: En !< The energy density integrated over an angular
1427  !! band [J m-2], intent in/out.
1428  real, dimension(G%isd:G%ied,G%JsdB:G%JedB), &
1429  intent(in) :: speed_y !< The magnitude of the group velocity at the
1430  !! Cv points [m s-1].
1431  real, dimension(Nangle), intent(in) :: Cgy_av !< The average y-projection in each angular band.
1432  real, dimension(Nangle), intent(in) :: dCgy !< The difference in y-projections between the
1433  !! edges of each angular band.
1434  real, intent(in) :: dt !< Time increment [s].
1435  type(int_tide_cs), pointer :: CS !< The control structure returned by a previous call
1436  !! to continuity_PPM_init.
1437  type(loop_bounds_type), intent(in) :: LB !< A structure with the active energy loop bounds.
1438  ! Local variables
1439  real, dimension(SZI_(G),SZJ_(G)) :: &
1440  EnL, EnR ! South and north face energy densities [J m-2].
1441  real, dimension(SZI_(G),SZJB_(G)) :: &
1442  flux_y ! The internal wave energy flux [J s-1].
1443  real, dimension(SZI_(G)) :: &
1444  cg_p, cg_m, flux1, flux2
1445  !real, dimension(SZI_(G),SZJB_(G),Nangle) :: En_m, En_p
1446  real, dimension(SZI_(G),SZJB_(G),Nangle) :: &
1447  Fdt_m, Fdt_p! South and north energy fluxes [J]
1448  character(len=160) :: mesg ! The text of an error message
1449  integer :: i, j, k, ish, ieh, jsh, jeh, a
1450 
1451  ish = lb%ish ; ieh = lb%ieh ; jsh = lb%jsh ; jeh = lb%jeh
1452  do a=1,nangle
1453  ! This sets EnL and EnR.
1454  if (cs%upwind_1st) then
1455  do j=jsh-1,jeh+1 ; do i=ish,ieh
1456  enl(i,j) = en(i,j,a) ; enr(i,j) = en(i,j,a)
1457  enddo ; enddo
1458  else
1459  call ppm_reconstruction_y(en(:,:,a), enl, enr, g, lb, simple_2nd=cs%simple_2nd)
1460  endif
1461 
1462  do j=jsh-1,jeh
1463  ! This is done once with single speed (GARDEN SPRINKLER EFFECT POSSIBLE)
1464  do i=ish,ieh
1465  cg_p(i) = speed_y(i,j) * (cgy_av(a))
1466  enddo
1467  call merid_flux_en(cg_p, en(:,:,a), enl(:,:), enr(:,:), flux1, &
1468  dt, g, j, ish, ieh, cs%vol_CFL)
1469  do i=ish,ieh ; flux_y(i,j) = flux1(i); enddo
1470  enddo
1471 
1472  do j=jsh,jeh ; do i=ish,ieh
1473  fdt_m(i,j,a) = dt*flux_y(i,j-1) ! south face influx (J)
1474  fdt_p(i,j,a) = -dt*flux_y(i,j) ! north face influx (J)
1475  !if ((En(i,j,a) + G%IareaT(i,j)*(Fdt_m(i,j,a) + Fdt_p(i,j,a))) < 0.0)then ! for debugging
1476  ! call MOM_error(WARNING, "propagate_y: OutFlux>Available prior to reflection", .true.)
1477  ! write(mesg,*) "flux_y_south=",flux_y(i,J-1),"flux_y_north=",flux_y(i,J),"En=",En(i,j,a), &
1478  ! "cn_south=", speed_y(i,J-1) * (Cgy_av(a)), "cn_north=", speed_y(i,J) * (Cgy_av(a))
1479  ! call MOM_error(WARNING, mesg, .true.)
1480  !endif
1481  enddo ; enddo
1482 
1483  ! test with old (take out later)
1484  !do j=jsh,jeh ; do i=ish,ieh
1485  ! En(i,j,a) = En(i,j,a) - dt* G%IareaT(i,j) * (flux_y(i,J) - flux_y(i,J-1))
1486  !enddo ; enddo
1487 
1488  enddo ! a-loop
1489 
1490  ! Only reflect newly arrived energy; existing energy in incident wedge
1491  ! is not reflected and will eventually propagate out of cell.
1492  ! (only reflects if En > 0)
1493  call reflect(fdt_m(:,:,:), nangle, cs, g, lb)
1494  call teleport(fdt_m(:,:,:), nangle, cs, g, lb)
1495  call reflect(fdt_p(:,:,:), nangle, cs, g, lb)
1496  call teleport(fdt_p(:,:,:), nangle, cs, g, lb)
1497 
1498  ! Update reflected energy (Jm-2)
1499  do j=jsh,jeh ; do i=ish,ieh
1500  !do a=1,CS%nAngle
1501  ! if ((En(i,j,a) + G%IareaT(i,j)*(Fdt_m(i,j,a) + Fdt_p(i,j,a))) < 0.0)then ! for debugging
1502  ! call MOM_error(FATAL, "propagate_y: OutFlux>Available", .true.)
1503  ! endif
1504  !enddo
1505  en(i,j,:) = en(i,j,:) + g%IareaT(i,j)*(fdt_m(i,j,:) + fdt_p(i,j,:))
1506  enddo ; enddo
1507 
1508 end subroutine propagate_y
1509 
1510 !> Evaluates the zonal mass or volume fluxes in a layer.
1511 subroutine zonal_flux_en(u, h, hL, hR, uh, dt, G, j, ish, ieh, vol_CFL)
1512  type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure.
1513  real, dimension(SZIB_(G)), intent(in) :: u !< The zonal velocity [m s-1].
1514  real, dimension(SZI_(G)), intent(in) :: h !< Energy density used to calculate the fluxes
1515  !! [J m-2].
1516  real, dimension(SZI_(G)), intent(in) :: hL !< Left- Energy densities in the reconstruction
1517  !! [J m-2].
1518  real, dimension(SZI_(G)), intent(in) :: hR !< Right- Energy densities in the reconstruction
1519  !! [J m-2].
1520  real, dimension(SZIB_(G)), intent(inout) :: uh !< The zonal energy transport [J s-1].
1521  real, intent(in) :: dt !< Time increment [s].
1522  integer, intent(in) :: j !< The j-index to work on.
1523  integer, intent(in) :: ish !< The start i-index range to work on.
1524  integer, intent(in) :: ieh !< The end i-index range to work on.
1525  logical, intent(in) :: vol_CFL !< If true, rescale the ratio of face areas to
1526  !! the cell areas when estimating the CFL number.
1527  ! Local variables
1528  real :: CFL ! The CFL number based on the local velocity and grid spacing [nondim].
1529  real :: curv_3 ! A measure of the thickness curvature over a grid length,
1530  ! with the same units as h_in.
1531  integer :: i
1532 
1533  do i=ish-1,ieh
1534  ! Set new values of uh and duhdu.
1535  if (u(i) > 0.0) then
1536  if (vol_cfl) then ; cfl = (u(i) * dt) * (g%dy_Cu(i,j) * g%IareaT(i,j))
1537  else ; cfl = u(i) * dt * g%IdxT(i,j) ; endif
1538  curv_3 = (hl(i) + hr(i)) - 2.0*h(i)
1539  uh(i) = g%dy_Cu(i,j) * u(i) * &
1540  (hr(i) + cfl * (0.5*(hl(i) - hr(i)) + curv_3*(cfl - 1.5)))
1541  elseif (u(i) < 0.0) then
1542  if (vol_cfl) then ; cfl = (-u(i) * dt) * (g%dy_Cu(i,j) * g%IareaT(i+1,j))
1543  else ; cfl = -u(i) * dt * g%IdxT(i+1,j) ; endif
1544  curv_3 = (hl(i+1) + hr(i+1)) - 2.0*h(i+1)
1545  uh(i) = g%dy_Cu(i,j) * u(i) * &
1546  (hl(i+1) + cfl * (0.5*(hr(i+1)-hl(i+1)) + curv_3*(cfl - 1.5)))
1547  else
1548  uh(i) = 0.0
1549  endif
1550  enddo
1551 end subroutine zonal_flux_en
1552 
1553 !> Evaluates the meridional mass or volume fluxes in a layer.
1554 subroutine merid_flux_en(v, h, hL, hR, vh, dt, G, J, ish, ieh, vol_CFL)
1555  type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure.
1556  real, dimension(SZI_(G)), intent(in) :: v !< The meridional velocity [m s-1].
1557  real, dimension(SZI_(G),SZJ_(G)), intent(in) :: h !< Energy density used to calculate the
1558  !! fluxes [J m-2].
1559  real, dimension(SZI_(G),SZJ_(G)), intent(in) :: hL !< Left- Energy densities in the
1560  !! reconstruction [J m-2].
1561  real, dimension(SZI_(G),SZJ_(G)), intent(in) :: hR !< Right- Energy densities in the
1562  !! reconstruction [J m-2].
1563  real, dimension(SZI_(G)), intent(inout) :: vh !< The meridional energy transport [J s-1].
1564  real, intent(in) :: dt !< Time increment [s].
1565  integer, intent(in) :: J !< The j-index to work on.
1566  integer, intent(in) :: ish !< The start i-index range to work on.
1567  integer, intent(in) :: ieh !< The end i-index range to work on.
1568  logical, intent(in) :: vol_CFL !< If true, rescale the ratio of face
1569  !! areas to the cell areas when estimating
1570  !! the CFL number.
1571  ! Local variables
1572  real :: CFL ! The CFL number based on the local velocity and grid spacing [nondim].
1573  real :: curv_3 ! A measure of the thickness curvature over a grid length,
1574  ! with the same units as h_in.
1575  integer :: i
1576 
1577  do i=ish,ieh
1578  if (v(i) > 0.0) then
1579  if (vol_cfl) then ; cfl = (v(i) * dt) * (g%dx_Cv(i,j) * g%IareaT(i,j))
1580  else ; cfl = v(i) * dt * g%IdyT(i,j) ; endif
1581  curv_3 = hl(i,j) + hr(i,j) - 2.0*h(i,j)
1582  vh(i) = g%dx_Cv(i,j) * v(i) * ( hr(i,j) + cfl * &
1583  (0.5*(hl(i,j) - hr(i,j)) + curv_3*(cfl - 1.5)) )
1584  elseif (v(i) < 0.0) then
1585  if (vol_cfl) then ; cfl = (-v(i) * dt) * (g%dx_Cv(i,j) * g%IareaT(i,j+1))
1586  else ; cfl = -v(i) * dt * g%IdyT(i,j+1) ; endif
1587  curv_3 = hl(i,j+1) + hr(i,j+1) - 2.0*h(i,j+1)
1588  vh(i) = g%dx_Cv(i,j) * v(i) * ( hl(i,j+1) + cfl * &
1589  (0.5*(hr(i,j+1)-hl(i,j+1)) + curv_3*(cfl - 1.5)) )
1590  else
1591  vh(i) = 0.0
1592  endif
1593  enddo
1594 end subroutine merid_flux_en
1595 
1596 !> Reflection of the internal waves at a single frequency.
1597 subroutine reflect(En, NAngle, CS, G, LB)
1598  type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure
1599  integer, intent(in) :: NAngle !< The number of wave orientations in the
1600  !! discretized wave energy spectrum.
1601  real, dimension(G%isd:G%ied,G%jsd:G%jed,NAngle), &
1602  intent(inout) :: En !< The internal gravity wave energy density as a
1603  !! function of space and angular resolution
1604  !! [J m-2 radian-1].
1605  type(int_tide_cs), pointer :: CS !< The control structure returned by a
1606  !! previous call to int_tide_init.
1607  type(loop_bounds_type), intent(in) :: LB !< A structure with the active energy loop bounds.
1608  ! Local variables
1609  real, dimension(G%isd:G%ied,G%jsd:G%jed) :: angle_c
1610  ! angle of boudary wrt equator
1611  real, dimension(G%isd:G%ied,G%jsd:G%jed) :: part_refl
1612  ! fraction of wave energy reflected
1613  ! values should collocate with angle_c
1614  logical, dimension(G%isd:G%ied,G%jsd:G%jed) :: ridge
1615  ! tags of cells with double reflection
1616 
1617  real :: TwoPi ! 2*pi
1618  real :: Angle_size ! size of beam wedge (rad)
1619  real :: angle_wall ! angle of coast/ridge/shelf wrt equator
1620  real, dimension(1:NAngle) :: angle_i ! angle of incident ray wrt equator
1621  real :: angle_r ! angle of reflected ray wrt equator
1622  real, dimension(1:Nangle) :: En_reflected
1623  integer :: i, j, a, a_r, na
1624  !integer :: isd, ied, jsd, jed ! start and end local indices on data domain
1625  ! ! (values include halos)
1626  integer :: isc, iec, jsc, jec ! start and end local indices on PE
1627  ! (values exclude halos)
1628  integer :: ish, ieh, jsh, jeh ! start and end local indices on data domain
1629  ! leaving out outdated halo points (march in)
1630  integer :: id_g, jd_g ! global (decomp-invar) indices
1631 
1632  !isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed
1633  isc = g%isc ; iec = g%iec ; jsc = g%jsc ; jec = g%jec
1634  ish = lb%ish ; ieh = lb%ieh ; jsh = lb%jsh ; jeh = lb%jeh
1635 
1636  twopi = 8.0*atan(1.0)
1637  angle_size = twopi / (real(nangle))
1638 
1639  do a=1,nangle
1640  ! These are the angles at the cell centers
1641  ! (should do this elsewhere since doesn't change with time)
1642  angle_i(a) = angle_size * real(a - 1) ! for a=1 aligned with x-axis
1643  enddo
1644 
1645  angle_c = cs%refl_angle
1646  part_refl = cs%refl_pref
1647  ridge = cs%refl_dbl
1648  en_reflected(:) = 0.0
1649 
1650  !do j=jsc-1,jec+1
1651  do j=jsh,jeh
1652  !do i=isc-1,iec+1
1653  do i=ish,ieh
1654  ! jd_g = j + G%jdg_offset ; id_g = i + G%idg_offset
1655  ! redistribute energy in angular space if ray will hit boundary
1656  ! i.e., if energy is in a reflecting cell
1657  if (angle_c(i,j) /= cs%nullangle) then
1658  do a=1,nangle
1659  if (en(i,j,a) > 0.0) then
1660  ! if ray is incident, keep specified boundary angle
1661  if (sin(angle_i(a) - angle_c(i,j)) >= 0.0) then
1662  angle_wall = angle_c(i,j)
1663  ! if ray is not incident but in ridge cell, use complementary angle
1664  elseif (ridge(i,j)) then
1665  angle_wall = angle_c(i,j) + 0.5*twopi
1666  if (angle_wall > twopi) then
1667  angle_wall = angle_wall - twopi*floor(abs(angle_wall)/twopi)
1668  elseif (angle_wall < 0.0) then
1669  angle_wall = angle_wall + twopi*ceiling(abs(angle_wall)/twopi)
1670  endif
1671  ! if ray is not incident and not in a ridge cell, keep specified angle
1672  else
1673  angle_wall = angle_c(i,j)
1674  endif
1675  ! do reflection
1676  if (sin(angle_i(a) - angle_wall) >= 0.0) then
1677  angle_r = 2.0*angle_wall - angle_i(a)
1678  if (angle_r > twopi) then
1679  angle_r = angle_r - twopi*floor(abs(angle_r)/twopi)
1680  elseif (angle_r < 0.0) then
1681  angle_r = angle_r + twopi*ceiling(abs(angle_r)/twopi)
1682  endif
1683  a_r = nint(angle_r/angle_size) + 1
1684  do while (a_r > nangle) ; a_r = a_r - nangle ; enddo
1685  if (a /= a_r) then
1686  en_reflected(a_r) = part_refl(i,j)*en(i,j,a)
1687  en(i,j,a) = (1.0-part_refl(i,j))*en(i,j,a)
1688  endif
1689  endif
1690  endif
1691  enddo ! a-loop
1692  en(i,j,:) = en(i,j,:) + en_reflected(:)
1693  en_reflected(:) = 0.0
1694  endif
1695  enddo ! i-loop
1696  enddo ! j-loop
1697 
1698  ! Check to make sure no energy gets onto land (only run for debugging)
1699  ! do a=1,NAngle ; do j=jsc,jec ; do i=isc,iec
1700  ! if (En(i,j,a) > 0.001 .and. G%mask2dT(i,j) == 0) then
1701  ! jd_g = j + G%jdg_offset ; id_g = i + G%idg_offset
1702  ! write (mesg,*) 'En=', En(i,j,a), 'a=', a, 'ig_g=',id_g, 'jg_g=',jd_g
1703  ! call MOM_error(FATAL, "reflect: Energy detected out of bounds: "//trim(mesg), .true.)
1704  ! endif
1705  ! enddo ; enddo ; enddo
1706 
1707 end subroutine reflect
1708 
1709 !> Moves energy across lines of partial reflection to prevent
1710 !! reflection of energy that is supposed to get across.
1711 subroutine teleport(En, NAngle, CS, G, LB)
1712  type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure.
1713  integer, intent(in) :: NAngle !< The number of wave orientations in the
1714  !! discretized wave energy spectrum.
1715  real, dimension(G%isd:G%ied,G%jsd:G%jed,NAngle), &
1716  intent(inout) :: En !< The internal gravity wave energy density as a
1717  !! function of space and angular resolution
1718  !! [J m-2 radian-1].
1719  type(int_tide_cs), pointer :: CS !< The control structure returned by a
1720  !! previous call to int_tide_init.
1721  type(loop_bounds_type), intent(in) :: LB !< A structure with the active energy loop bounds.
1722  ! Local variables
1723  real, dimension(G%isd:G%ied,G%jsd:G%jed) :: angle_c
1724  ! angle of boudary wrt equator
1725  real, dimension(G%isd:G%ied,G%jsd:G%jed) :: part_refl
1726  ! fraction of wave energy reflected
1727  ! values should collocate with angle_c
1728  logical, dimension(G%isd:G%ied,G%jsd:G%jed) :: pref_cell
1729  ! flag for partial reflection
1730  logical, dimension(G%isd:G%ied,G%jsd:G%jed) :: ridge
1731  ! tags of cells with double reflection
1732  real :: TwoPi ! size of beam wedge (rad)
1733  real :: Angle_size ! size of beam wedge (rad)
1734  real, dimension(1:NAngle) :: angle_i ! angle of incident ray wrt equator
1735  real, dimension(1:NAngle) :: cos_angle, sin_angle
1736  real :: En_tele ! energy to be "teleported"
1737  character(len=160) :: mesg ! The text of an error message
1738  integer :: i, j, a
1739  !integer :: isd, ied, jsd, jed ! start and end local indices on data domain
1740  ! ! (values include halos)
1741  !integer :: isc, iec, jsc, jec ! start and end local indices on PE
1742  ! ! (values exclude halos)
1743  integer :: ish, ieh, jsh, jeh ! start and end local indices on data domain
1744  ! leaving out outdated halo points (march in)
1745  integer :: id_g, jd_g ! global (decomp-invar) indices
1746  integer :: jos, ios ! offsets
1747  real :: cos_normal, sin_normal, angle_wall
1748  ! cos/sin of cross-ridge normal, ridge angle
1749 
1750  !isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed
1751  !isc = G%isc ; iec = G%iec ; jsc = G%jsc ; jec = G%jec
1752  ish = lb%ish ; ieh = lb%ieh ; jsh = lb%jsh ; jeh = lb%jeh
1753 
1754  twopi = 8.0*atan(1.0)
1755  angle_size = twopi / (real(nangle))
1756 
1757  do a=1,nangle
1758  ! These are the angles at the cell centers
1759  ! (should do this elsewhere since doesn't change with time)
1760  angle_i(a) = angle_size * real(a - 1) ! for a=1 aligned with x-axis
1761  cos_angle(a) = cos(angle_i(a)) ; sin_angle(a) = sin(angle_i(a))
1762  enddo
1763 
1764  angle_c = cs%refl_angle
1765  part_refl = cs%refl_pref
1766  pref_cell = cs%refl_pref_logical
1767  ridge = cs%refl_dbl
1768 
1769  do j=jsh,jeh
1770  do i=ish,ieh
1771  id_g = i + g%idg_offset ; jd_g = j + g%jdg_offset
1772  if (pref_cell(i,j)) then
1773  do a=1,nangle
1774  if (en(i,j,a) > 0) then
1775  ! if ray is incident, keep specified boundary angle
1776  if (sin(angle_i(a) - angle_c(i,j)) >= 0.0) then
1777  angle_wall = angle_c(i,j)
1778  ! if ray is not incident but in ridge cell, use complementary angle
1779  elseif (ridge(i,j)) then
1780  angle_wall = angle_c(i,j) + 0.5*twopi
1781  ! if ray is not incident and not in a ridge cell, keep specified angle
1782  else
1783  angle_wall = angle_c(i,j)
1784  endif
1785  ! teleport if incident
1786  if (sin(angle_i(a) - angle_wall) >= 0.0) then
1787  en_tele = en(i,j,a)
1788  cos_normal = cos(angle_wall + 0.25*twopi)
1789  sin_normal = sin(angle_wall + 0.25*twopi)
1790  ! find preferred zonal offset based on shelf/ridge angle
1791  ios = int(sign(1.,cos_normal))
1792  ! find preferred meridional offset based on shelf/ridge angle
1793  jos = int(sign(1.,sin_normal))
1794  ! find receptive ocean cell in direction of offset
1795  if (.not. pref_cell(i+ios,j+jos)) then
1796  en(i,j,a) = en(i,j,a) - en_tele
1797  en(i+ios,j+jos,a) = en(i+ios,j+jos,a) + en_tele
1798  else
1799  write(mesg,*) 'idg=',id_g,'jd_g=',jd_g,'a=',a
1800  call mom_error(fatal, "teleport: no receptive ocean cell at "//trim(mesg), .true.)
1801  endif
1802  endif ! incidence check
1803  endif ! energy check
1804  enddo ! a-loop
1805  endif ! pref check
1806  enddo ! i-loop
1807  enddo ! j-loop
1808 
1809 end subroutine teleport
1810 
1811 !> Rotates points in the halos where required to accommodate
1812 !! changes in grid orientation, such as at the tripolar fold.
1813 subroutine correct_halo_rotation(En, test, G, NAngle)
1814  type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure
1815  real, dimension(:,:,:,:,:), intent(inout) :: En !< The internal gravity wave energy density as a
1816  !! function of space, angular orientation, frequency,
1817  !! and vertical mode [J m-2 radian-1].
1818  real, dimension(SZI_(G),SZJ_(G),2), &
1819  intent(in) :: test !< An x-unit vector that has been passed through
1820  !! the halo updates, to enable the rotation of the
1821  !! wave energies in the halo region to be corrected.
1822  integer, intent(in) :: NAngle !< The number of wave orientations in the
1823  !! discretized wave energy spectrum.
1824  ! Local variables
1825  real, dimension(G%isd:G%ied,NAngle) :: En2d
1826  integer, dimension(G%isd:G%ied) :: a_shift
1827  integer :: i_first, i_last, a_new
1828  integer :: a, i, j, isd, ied, jsd, jed, m, fr
1829  character(len=160) :: mesg ! The text of an error message
1830  isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
1831 
1832  do j=jsd,jed
1833  i_first = ied+1 ; i_last = isd-1
1834  do i=isd,ied
1835  a_shift(i) = 0
1836  if (test(i,j,1) /= 1.0) then
1837  if (i<i_first) i_first = i
1838  if (i>i_last) i_last = i
1839 
1840  if (test(i,j,1) == -1.0) then ; a_shift(i) = nangle/2
1841  elseif (test(i,j,2) == 1.0) then ; a_shift(i) = -nangle/4
1842  elseif (test(i,j,2) == -1.0) then ; a_shift(i) = nangle/4
1843  else
1844  write(mesg,'("Unrecognized rotation test vector ",2ES9.2," at ",F7.2," E, ",&
1845  &F7.2," N; i,j=",2i4)') &
1846  test(i,j,1), test(i,j,2), g%GeoLonT(i,j), g%GeoLatT(i,j), i, j
1847  call mom_error(fatal, mesg)
1848  endif
1849  endif
1850  enddo
1851 
1852  if (i_first <= i_last) then
1853  ! At least one point in this row needs to be rotated.
1854  do m=1,size(en,5) ; do fr=1,size(en,4)
1855  do a=1,nangle ; do i=i_first,i_last ; if (a_shift(i) /= 0) then
1856  a_new = a + a_shift(i)
1857  if (a_new < 1) a_new = a_new + nangle
1858  if (a_new > nangle) a_new = a_new - nangle
1859  en2d(i,a_new) = en(i,j,a,fr,m)
1860  endif ; enddo ; enddo
1861  do a=1,nangle ; do i=i_first,i_last ; if (a_shift(i) /= 0) then
1862  en(i,j,a,fr,m) = en2d(i,a)
1863  endif ; enddo ; enddo
1864  enddo ; enddo
1865  endif
1866  enddo
1867 end subroutine correct_halo_rotation
1868 
1869 !> Calculates left/right edge values for PPM reconstruction in x-direction.
1870 subroutine ppm_reconstruction_x(h_in, h_l, h_r, G, LB, simple_2nd)
1871  type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure.
1872  real, dimension(SZI_(G),SZJ_(G)), intent(in) :: h_in !< Energy density in a sector (2D).
1873  real, dimension(SZI_(G),SZJ_(G)), intent(out) :: h_l !< Left edge value of reconstruction (2D).
1874  real, dimension(SZI_(G),SZJ_(G)), intent(out) :: h_r !< Right edge value of reconstruction (2D).
1875  type(loop_bounds_type), intent(in) :: LB !< A structure with the active loop bounds.
1876  logical, optional, intent(in) :: simple_2nd !< If true, use the arithmetic mean
1877  !! energy densities as default edge values
1878  !! for a simple 2nd order scheme.
1879  ! Local variables
1880  real, dimension(SZI_(G),SZJ_(G)) :: slp ! The slopes.
1881  real, parameter :: oneSixth = 1./6.
1882  real :: h_ip1, h_im1
1883  real :: dMx, dMn
1884  logical :: use_CW84, use_2nd
1885  character(len=256) :: mesg ! The text of an error message
1886  integer :: i, j, isl, iel, jsl, jel, stencil
1887 
1888  use_2nd = .false. ; if (present(simple_2nd)) use_2nd = simple_2nd
1889  isl = lb%ish-1 ; iel = lb%ieh+1 ; jsl = lb%jsh ; jel = lb%jeh
1890 
1891  ! This is the stencil of the reconstruction, not the scheme overall.
1892  stencil = 2 ; if (use_2nd) stencil = 1
1893 
1894  if ((isl-stencil < g%isd) .or. (iel+stencil > g%ied)) then
1895  write(mesg,'("In MOM_internal_tides, PPM_reconstruction_x called with a ", &
1896  & "x-halo that needs to be increased by ",i2,".")') &
1897  stencil + max(g%isd-isl,iel-g%ied)
1898  call mom_error(fatal,mesg)
1899  endif
1900  if ((jsl < g%jsd) .or. (jel > g%jed)) then
1901  write(mesg,'("In MOM_internal_tides, PPM_reconstruction_x called with a ", &
1902  & "y-halo that needs to be increased by ",i2,".")') &
1903  max(g%jsd-jsl,jel-g%jed)
1904  call mom_error(fatal,mesg)
1905  endif
1906 
1907  if (use_2nd) then
1908  do j=jsl,jel ; do i=isl,iel
1909  h_im1 = g%mask2dT(i-1,j) * h_in(i-1,j) + (1.0-g%mask2dT(i-1,j)) * h_in(i,j)
1910  h_ip1 = g%mask2dT(i+1,j) * h_in(i+1,j) + (1.0-g%mask2dT(i+1,j)) * h_in(i,j)
1911  h_l(i,j) = 0.5*( h_im1 + h_in(i,j) )
1912  h_r(i,j) = 0.5*( h_ip1 + h_in(i,j) )
1913  enddo ; enddo
1914  else
1915  do j=jsl,jel ; do i=isl-1,iel+1
1916  if ((g%mask2dT(i-1,j) * g%mask2dT(i,j) * g%mask2dT(i+1,j)) == 0.0) then
1917  slp(i,j) = 0.0
1918  else
1919  ! This uses a simple 2nd order slope.
1920  slp(i,j) = 0.5 * (h_in(i+1,j) - h_in(i-1,j))
1921  ! Monotonic constraint, see Eq. B2 in Lin 1994, MWR (132)
1922  dmx = max(h_in(i+1,j), h_in(i-1,j), h_in(i,j)) - h_in(i,j)
1923  dmn = h_in(i,j) - min(h_in(i+1,j), h_in(i-1,j), h_in(i,j))
1924  slp(i,j) = sign(1.,slp(i,j)) * min(abs(slp(i,j)), 2. * min(dmx, dmn))
1925  ! * (G%mask2dT(i-1,j) * G%mask2dT(i,j) * G%mask2dT(i+1,j))
1926  endif
1927  enddo ; enddo
1928 
1929  do j=jsl,jel ; do i=isl,iel
1930  ! Neighboring values should take into account any boundaries. The 3
1931  ! following sets of expressions are equivalent.
1932  ! h_im1 = h_in(i-1,j,k) ; if (G%mask2dT(i-1,j) < 0.5) h_im1 = h_in(i,j)
1933  ! h_ip1 = h_in(i+1,j,k) ; if (G%mask2dT(i+1,j) < 0.5) h_ip1 = h_in(i,j)
1934  h_im1 = g%mask2dT(i-1,j) * h_in(i-1,j) + (1.0-g%mask2dT(i-1,j)) * h_in(i,j)
1935  h_ip1 = g%mask2dT(i+1,j) * h_in(i+1,j) + (1.0-g%mask2dT(i+1,j)) * h_in(i,j)
1936  ! Left/right values following Eq. B2 in Lin 1994, MWR (132)
1937  h_l(i,j) = 0.5*( h_im1 + h_in(i,j) ) + onesixth*( slp(i-1,j) - slp(i,j) )
1938  h_r(i,j) = 0.5*( h_ip1 + h_in(i,j) ) + onesixth*( slp(i,j) - slp(i+1,j) )
1939  enddo ; enddo
1940  endif
1941 
1942  call ppm_limit_pos(h_in, h_l, h_r, 0.0, g, isl, iel, jsl, jel)
1943 end subroutine ppm_reconstruction_x
1944 
1945 !> Calculates left/right edge valus for PPM reconstruction in y-direction.
1946 subroutine ppm_reconstruction_y(h_in, h_l, h_r, G, LB, simple_2nd)
1947  type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure.
1948  real, dimension(SZI_(G),SZJ_(G)), intent(in) :: h_in !< Energy density in a sector (2D).
1949  real, dimension(SZI_(G),SZJ_(G)), intent(out) :: h_l !< Left edge value of reconstruction (2D).
1950  real, dimension(SZI_(G),SZJ_(G)), intent(out) :: h_r !< Right edge value of reconstruction (2D).
1951  type(loop_bounds_type), intent(in) :: LB !< A structure with the active loop bounds.
1952  logical, optional, intent(in) :: simple_2nd !< If true, use the arithmetic mean
1953  !! energy densities as default edge values
1954  !! for a simple 2nd order scheme.
1955  ! Local variables
1956  real, dimension(SZI_(G),SZJ_(G)) :: slp ! The slopes.
1957  real, parameter :: oneSixth = 1./6.
1958  real :: h_jp1, h_jm1
1959  real :: dMx, dMn
1960  logical :: use_2nd
1961  character(len=256) :: mesg ! The text of an error message
1962  integer :: i, j, isl, iel, jsl, jel, stencil
1963 
1964  use_2nd = .false. ; if (present(simple_2nd)) use_2nd = simple_2nd
1965  isl = lb%ish ; iel = lb%ieh ; jsl = lb%jsh-1 ; jel = lb%jeh+1
1966 
1967  ! This is the stencil of the reconstruction, not the scheme overall.
1968  stencil = 2 ; if (use_2nd) stencil = 1
1969 
1970  if ((isl < g%isd) .or. (iel > g%ied)) then
1971  write(mesg,'("In MOM_internal_tides, PPM_reconstruction_y called with a ", &
1972  & "x-halo that needs to be increased by ",i2,".")') &
1973  max(g%isd-isl,iel-g%ied)
1974  call mom_error(fatal,mesg)
1975  endif
1976  if ((jsl-stencil < g%jsd) .or. (jel+stencil > g%jed)) then
1977  write(mesg,'("In MOM_internal_tides, PPM_reconstruction_y called with a ", &
1978  & "y-halo that needs to be increased by ",i2,".")') &
1979  stencil + max(g%jsd-jsl,jel-g%jed)
1980  call mom_error(fatal,mesg)
1981  endif
1982 
1983  if (use_2nd) then
1984  do j=jsl,jel ; do i=isl,iel
1985  h_jm1 = g%mask2dT(i,j-1) * h_in(i,j-1) + (1.0-g%mask2dT(i,j-1)) * h_in(i,j)
1986  h_jp1 = g%mask2dT(i,j+1) * h_in(i,j+1) + (1.0-g%mask2dT(i,j+1)) * h_in(i,j)
1987  h_l(i,j) = 0.5*( h_jm1 + h_in(i,j) )
1988  h_r(i,j) = 0.5*( h_jp1 + h_in(i,j) )
1989  enddo ; enddo
1990  else
1991  do j=jsl-1,jel+1 ; do i=isl,iel
1992  if ((g%mask2dT(i,j-1) * g%mask2dT(i,j) * g%mask2dT(i,j+1)) == 0.0) then
1993  slp(i,j) = 0.0
1994  else
1995  ! This uses a simple 2nd order slope.
1996  slp(i,j) = 0.5 * (h_in(i,j+1) - h_in(i,j-1))
1997  ! Monotonic constraint, see Eq. B2 in Lin 1994, MWR (132)
1998  dmx = max(h_in(i,j+1), h_in(i,j-1), h_in(i,j)) - h_in(i,j)
1999  dmn = h_in(i,j) - min(h_in(i,j+1), h_in(i,j-1), h_in(i,j))
2000  slp(i,j) = sign(1.,slp(i,j)) * min(abs(slp(i,j)), 2. * min(dmx, dmn))
2001  ! * (G%mask2dT(i,j-1) * G%mask2dT(i,j) * G%mask2dT(i,j+1))
2002  endif
2003  enddo ; enddo
2004 
2005  do j=jsl,jel ; do i=isl,iel
2006  ! Neighboring values should take into account any boundaries. The 3
2007  ! following sets of expressions are equivalent.
2008  h_jm1 = g%mask2dT(i,j-1) * h_in(i,j-1) + (1.0-g%mask2dT(i,j-1)) * h_in(i,j)
2009  h_jp1 = g%mask2dT(i,j+1) * h_in(i,j+1) + (1.0-g%mask2dT(i,j+1)) * h_in(i,j)
2010  ! Left/right values following Eq. B2 in Lin 1994, MWR (132)
2011  h_l(i,j) = 0.5*( h_jm1 + h_in(i,j) ) + onesixth*( slp(i,j-1) - slp(i,j) )
2012  h_r(i,j) = 0.5*( h_jp1 + h_in(i,j) ) + onesixth*( slp(i,j) - slp(i,j+1) )
2013  enddo ; enddo
2014  endif
2015 
2016  call ppm_limit_pos(h_in, h_l, h_r, 0.0, g, isl, iel, jsl, jel)
2017 end subroutine ppm_reconstruction_y
2018 
2019 !> Limits the left/right edge values of the PPM reconstruction
2020 !! to give a reconstruction that is positive-definite. Here this is
2021 !! reinterpreted as giving a constant thickness if the mean thickness is less
2022 !! than h_min, with a minimum of h_min otherwise.
2023 subroutine ppm_limit_pos(h_in, h_L, h_R, h_min, G, iis, iie, jis, jie)
2024  type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure.
2025  real, dimension(SZI_(G),SZJ_(G)), intent(in) :: h_in !< Thickness of layer (2D).
2026  real, dimension(SZI_(G),SZJ_(G)), intent(inout) :: h_L !< Left edge value (2D).
2027  real, dimension(SZI_(G),SZJ_(G)), intent(inout) :: h_R !< Right edge value (2D).
2028  real, intent(in) :: h_min !< The minimum thickness that can be
2029  !! obtained by a concave parabolic fit.
2030  integer, intent(in) :: iis !< Start i-index for computations
2031  integer, intent(in) :: iie !< End i-index for computations
2032  integer, intent(in) :: jis !< Start j-index for computations
2033  integer, intent(in) :: jie !< End j-index for computations
2034  ! Local variables
2035  real :: curv, dh, scale
2036  character(len=256) :: mesg ! The text of an error message
2037  integer :: i,j
2038 
2039  do j=jis,jie ; do i=iis,iie
2040  ! This limiter prevents undershooting minima within the domain with
2041  ! values less than h_min.
2042  curv = 3.0*(h_l(i,j) + h_r(i,j) - 2.0*h_in(i,j))
2043  if (curv > 0.0) then ! Only minima are limited.
2044  dh = h_r(i,j) - h_l(i,j)
2045  if (abs(dh) < curv) then ! The parabola's minimum is within the cell.
2046  if (h_in(i,j) <= h_min) then
2047  h_l(i,j) = h_in(i,j) ; h_r(i,j) = h_in(i,j)
2048  elseif (12.0*curv*(h_in(i,j) - h_min) < (curv**2 + 3.0*dh**2)) then
2049  ! The minimum value is h_in - (curv^2 + 3*dh^2)/(12*curv), and must
2050  ! be limited in this case. 0 < scale < 1.
2051  scale = 12.0*curv*(h_in(i,j) - h_min) / (curv**2 + 3.0*dh**2)
2052  h_l(i,j) = h_in(i,j) + scale*(h_l(i,j) - h_in(i,j))
2053  h_r(i,j) = h_in(i,j) + scale*(h_r(i,j) - h_in(i,j))
2054  endif
2055  endif
2056  endif
2057  enddo ; enddo
2058 end subroutine ppm_limit_pos
2059 
2060 ! subroutine register_int_tide_restarts(G, param_file, CS, restart_CS)
2061 ! type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure
2062 ! type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters
2063 ! type(int_tide_CS), pointer :: CS !< The control structure returned by a
2064 ! !! previous call to int_tide_init.
2065 ! type(MOM_restart_CS), pointer :: restart_CS !< A pointer to the restart control structure.
2066 
2067 ! ! This subroutine is not currently in use!!
2068 
2069 ! ! This subroutine is used to allocate and register any fields in this module
2070 ! ! that should be written to or read from the restart file.
2071 ! logical :: use_int_tides
2072 ! type(vardesc) :: vd
2073 ! integer :: num_freq, num_angle , num_mode, period_1
2074 ! integer :: isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB, a
2075 ! isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed
2076 ! IsdB = G%IsdB ; IedB = G%IedB ; JsdB = G%JsdB ; JedB = G%JedB
2077 
2078 ! if (associated(CS)) then
2079 ! call MOM_error(WARNING, "register_int_tide_restarts called "//&
2080 ! "with an associated control structure.")
2081 ! return
2082 ! endif
2083 
2084 ! use_int_tides = .false.
2085 ! call read_param(param_file, "INTERNAL_TIDES", use_int_tides)
2086 ! if (.not.use_int_tides) return
2087 
2088 ! allocate(CS)
2089 
2090 ! num_angle = 24
2091 ! call read_param(param_file, "INTERNAL_TIDE_ANGLES", num_angle)
2092 ! allocate(CS%En_restart(isd:ied, jsd:jed, num_angle))
2093 ! CS%En_restart(:,:,:) = 0.0
2094 
2095 ! vd = vardesc("En_restart", &
2096 ! "The internal wave energy density as a function of (i,j,angle,frequency,mode)", &
2097 ! 'h','1','1',"J m-2")
2098 ! call register_restart_field(CS%En_restart, vd, .false., restart_CS)
2099 
2100 ! end subroutine register_int_tide_restarts
2101 
2102 !> This subroutine initializes the internal tides module.
2103 subroutine internal_tides_init(Time, G, GV, US, param_file, diag, CS)
2104  type(time_type), target, intent(in) :: time !< The current model time.
2105  type(ocean_grid_type), intent(inout) :: g !< The ocean's grid structure.
2106  type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure.
2107  type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
2108  type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time
2109  !! parameters.
2110  type(diag_ctrl), target, intent(in) :: diag !< A structure that is used to regulate
2111  !! diagnostic output.
2112  type(int_tide_cs),pointer :: cs !< A pointer that is set to point to the control
2113  !! structure for this module.
2114  ! Local variables
2115  real :: angle_size ! size of wedges, rad
2116  real, allocatable :: angles(:) ! orientations of wedge centers, rad
2117  real, allocatable, dimension(:,:) :: h2 ! topographic roughness scale, m^2
2118  real :: kappa_itides, kappa_h2_factor
2119  ! characteristic topographic wave number
2120  ! and a scaling factor
2121  real, allocatable :: ridge_temp(:,:)
2122  ! array for temporary storage of flags
2123  ! of cells with double-reflecting ridges
2124  logical :: use_int_tides, use_temperature
2125  integer :: num_angle, num_freq, num_mode, m, fr, period_1
2126  integer :: isd, ied, jsd, jed, a, id_ang, i, j
2127  type(axes_grp) :: axes_ang
2128  ! This include declares and sets the variable "version".
2129 #include "version_variable.h"
2130  character(len=40) :: mdl = "MOM_internal_tides" ! This module's name.
2131  character(len=16), dimension(8) :: freq_name
2132  character(len=40) :: var_name
2133  character(len=160) :: var_descript
2134  character(len=200) :: filename
2135  character(len=200) :: refl_angle_file, land_mask_file
2136  character(len=200) :: refl_pref_file, refl_dbl_file
2137  character(len=200) :: dy_cu_file, dx_cv_file
2138  character(len=200) :: h2_file
2139 
2140  isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
2141 
2142  if (associated(cs)) then
2143  call mom_error(warning, "internal_tides_init called "//&
2144  "with an associated control structure.")
2145  return
2146  else
2147  allocate(cs)
2148  endif
2149 
2150  use_int_tides = .false.
2151  call read_param(param_file, "INTERNAL_TIDES", use_int_tides)
2152  cs%do_int_tides = use_int_tides
2153  if (.not.use_int_tides) return
2154 
2155  use_temperature = .true.
2156  call read_param(param_file, "ENABLE_THERMODYNAMICS", use_temperature)
2157  if (.not.use_temperature) call mom_error(fatal, &
2158  "register_int_tide_restarts: internal_tides only works with "//&
2159  "ENABLE_THERMODYNAMICS defined.")
2160 
2161  ! Set number of frequencies, angles, and modes to consider
2162  num_freq = 1 ; num_angle = 24 ; num_mode = 1
2163  call read_param(param_file, "INTERNAL_TIDE_FREQS", num_freq)
2164  call read_param(param_file, "INTERNAL_TIDE_ANGLES", num_angle)
2165  call read_param(param_file, "INTERNAL_TIDE_MODES", num_mode)
2166  if (.not.((num_freq > 0) .and. (num_angle > 0) .and. (num_mode > 0))) return
2167  cs%nFreq = num_freq ; cs%nAngle = num_angle ; cs%nMode = num_mode
2168 
2169  ! Allocate energy density array
2170  allocate(cs%En(isd:ied, jsd:jed, num_angle, num_freq, num_mode))
2171  cs%En(:,:,:,:,:) = 0.0
2172 
2173  ! Allocate phase speed array
2174  allocate(cs%cp(isd:ied, jsd:jed, num_freq, num_mode))
2175  cs%cp(:,:,:,:) = 0.0
2176 
2177  ! Allocate and populate frequency array (each a multiple of first for now)
2178  allocate(cs%frequency(num_freq))
2179  call read_param(param_file, "FIRST_MODE_PERIOD", period_1); ! ADDED BDM
2180  do fr=1,num_freq
2181  cs%frequency(fr) = (8.0*atan(1.0) * (real(fr)) / period_1) ! ADDED BDM
2182  enddo
2183 
2184  ! Read all relevant parameters and write them to the model log.
2185 
2186  cs%Time => time ! direct a pointer to the current model time target
2187 
2188  call get_param(param_file, mdl, "INPUTDIR", cs%inputdir, default=".")
2189  cs%inputdir = slasher(cs%inputdir)
2190 
2191  call log_version(param_file, mdl, version, "")
2192 
2193  call get_param(param_file, mdl, "INTERNAL_TIDE_FREQS", num_freq, &
2194  "The number of distinct internal tide frequency bands "//&
2195  "that will be calculated.", default=1)
2196  call get_param(param_file, mdl, "INTERNAL_TIDE_MODES", num_mode, &
2197  "The number of distinct internal tide modes "//&
2198  "that will be calculated.", default=1)
2199  call get_param(param_file, mdl, "INTERNAL_TIDE_ANGLES", num_angle, &
2200  "The number of angular resolution bands for the internal "//&
2201  "tide calculations.", default=24)
2202 
2203  if (use_int_tides) then
2204  if ((num_freq <= 0) .and. (num_mode <= 0) .and. (num_angle <= 0)) then
2205  call mom_error(warning, "Internal tides were enabled, but the number "//&
2206  "of requested frequencies, modes and angles were not all positive.")
2207  return
2208  endif
2209  else
2210  if ((num_freq > 0) .and. (num_mode > 0) .and. (num_angle > 0)) then
2211  call mom_error(warning, "Internal tides were not enabled, even though "//&
2212  "the number of requested frequencies, modes and angles were all "//&
2213  "positive.")
2214  return
2215  endif
2216  endif
2217 
2218  if (cs%NFreq /= num_freq) call mom_error(fatal, "Internal_tides_init: "//&
2219  "Inconsistent number of frequencies.")
2220  if (cs%NAngle /= num_angle) call mom_error(fatal, "Internal_tides_init: "//&
2221  "Inconsistent number of angles.")
2222  if (cs%NMode /= num_mode) call mom_error(fatal, "Internal_tides_init: "//&
2223  "Inconsistent number of modes.")
2224  if (4*(num_angle/4) /= num_angle) call mom_error(fatal, &
2225  "Internal_tides_init: INTERNAL_TIDE_ANGLES must be a multiple of 4.")
2226 
2227  cs%diag => diag
2228 
2229  call get_param(param_file, mdl, "INTERNAL_TIDE_DECAY_RATE", cs%decay_rate, &
2230  "The rate at which internal tide energy is lost to the "//&
2231  "interior ocean internal wave field.", units="s-1", default=0.0)
2232  call get_param(param_file, mdl, "INTERNAL_TIDE_VOLUME_BASED_CFL", cs%vol_CFL, &
2233  "If true, use the ratio of the open face lengths to the "//&
2234  "tracer cell areas when estimating CFL numbers in the "//&
2235  "internal tide code.", default=.false.)
2236  call get_param(param_file, mdl, "INTERNAL_TIDE_CORNER_ADVECT", cs%corner_adv, &
2237  "If true, internal tide ray-tracing advection uses a "//&
2238  "corner-advection scheme rather than PPM.", default=.false.)
2239  call get_param(param_file, mdl, "INTERNAL_TIDE_SIMPLE_2ND_PPM", cs%simple_2nd, &
2240  "If true, CONTINUITY_PPM uses a simple 2nd order "//&
2241  "(arithmetic mean) interpolation of the edge values. "//&
2242  "This may give better PV conservation properties. While "//&
2243  "it formally reduces the accuracy of the continuity "//&
2244  "solver itself in the strongly advective limit, it does "//&
2245  "not reduce the overall order of accuracy of the dynamic "//&
2246  "core.", default=.false.)
2247  call get_param(param_file, mdl, "INTERNAL_TIDE_UPWIND_1ST", cs%upwind_1st, &
2248  "If true, the internal tide ray-tracing advection uses "//&
2249  "1st-order upwind advection. This scheme is highly "//&
2250  "continuity solver. This scheme is highly "//&
2251  "diffusive but may be useful for debugging.", default=.false.)
2252  call get_param(param_file, mdl, "INTERNAL_TIDE_BACKGROUND_DRAG", &
2253  cs%apply_background_drag, "If true, the internal tide "//&
2254  "ray-tracing advection uses a background drag term as a sink.",&
2255  default=.false.)
2256  call get_param(param_file, mdl, "INTERNAL_TIDE_QUAD_DRAG", cs%apply_bottom_drag, &
2257  "If true, the internal tide ray-tracing advection uses "//&
2258  "a quadratic bottom drag term as a sink.", default=.false.)
2259  call get_param(param_file, mdl, "INTERNAL_TIDE_WAVE_DRAG", cs%apply_wave_drag, &
2260  "If true, apply scattering due to small-scale roughness as a sink.", &
2261  default=.false.)
2262  call get_param(param_file, mdl, "INTERNAL_TIDE_FROUDE_DRAG", cs%apply_Froude_drag, &
2263  "If true, apply wave breaking as a sink.", &
2264  default=.false.)
2265  call get_param(param_file, mdl, "CDRAG", cs%cdrag, &
2266  "CDRAG is the drag coefficient relating the magnitude of "//&
2267  "the velocity field to the bottom stress.", units="nondim", &
2268  default=0.003)
2269  call get_param(param_file, mdl, "INTERNAL_TIDE_ENERGIZED_ANGLE", cs%energized_angle, &
2270  "If positive, only one angular band of the internal tides "//&
2271  "gets all of the energy. (This is for debugging.)", default=-1)
2272  call get_param(param_file, mdl, "USE_PPM_ANGULAR", cs%use_PPMang, &
2273  "If true, use PPM for advection of energy in angular space.", &
2274  default=.false.)
2275  call get_param(param_file, mdl, "GAMMA_ITIDES", cs%q_itides, &
2276  "The fraction of the internal tidal energy that is "//&
2277  "dissipated locally with INT_TIDE_DISSIPATION. "//&
2278  "THIS NAME COULD BE BETTER.", &
2279  units="nondim", default=0.3333)
2280  call get_param(param_file, mdl, "KAPPA_ITIDES", kappa_itides, &
2281  "A topographic wavenumber used with INT_TIDE_DISSIPATION. "//&
2282  "The default is 2pi/10 km, as in St.Laurent et al. 2002.", &
2283  units="m-1", default=8.e-4*atan(1.0))
2284  call get_param(param_file, mdl, "KAPPA_H2_FACTOR", kappa_h2_factor, &
2285  "A scaling factor for the roughness amplitude with n"//&
2286  "INT_TIDE_DISSIPATION.", units="nondim", default=1.0)
2287 
2288  ! Allocate various arrays needed for loss rates
2289  allocate(h2(isd:ied,jsd:jed)) ; h2(:,:) = 0.0
2290  allocate(cs%TKE_itidal_loss_fixed(isd:ied,jsd:jed))
2291  cs%TKE_itidal_loss_fixed = 0.0
2292  allocate(cs%TKE_leak_loss(isd:ied,jsd:jed,num_angle,num_freq,num_mode))
2293  cs%TKE_leak_loss(:,:,:,:,:) = 0.0
2294  allocate(cs%TKE_quad_loss(isd:ied,jsd:jed,num_angle,num_freq,num_mode))
2295  cs%TKE_quad_loss(:,:,:,:,:) = 0.0
2296  allocate(cs%TKE_itidal_loss(isd:ied,jsd:jed,num_angle,num_freq,num_mode))
2297  cs%TKE_itidal_loss(:,:,:,:,:) = 0.0
2298  allocate(cs%TKE_Froude_loss(isd:ied,jsd:jed,num_angle,num_freq,num_mode))
2299  cs%TKE_Froude_loss(:,:,:,:,:) = 0.0
2300  allocate(cs%tot_leak_loss(isd:ied,jsd:jed)) ; cs%tot_leak_loss(:,:) = 0.0
2301  allocate(cs%tot_quad_loss(isd:ied,jsd:jed) ) ; cs%tot_quad_loss(:,:) = 0.0
2302  allocate(cs%tot_itidal_loss(isd:ied,jsd:jed)) ; cs%tot_itidal_loss(:,:) = 0.0
2303  allocate(cs%tot_Froude_loss(isd:ied,jsd:jed)) ; cs%tot_Froude_loss(:,:) = 0.0
2304 
2305  ! Compute the fixed part of the bottom drag loss from baroclinic modes
2306  call get_param(param_file, mdl, "H2_FILE", h2_file, &
2307  "The path to the file containing the sub-grid-scale "//&
2308  "topographic roughness amplitude with INT_TIDE_DISSIPATION.", &
2309  fail_if_missing=.true.)
2310  filename = trim(cs%inputdir) // trim(h2_file)
2311  call log_param(param_file, mdl, "INPUTDIR/H2_FILE", filename)
2312  call mom_read_data(filename, 'h2', h2, g%domain, timelevel=1, scale=us%m_to_Z)
2313  do j=g%jsc,g%jec ; do i=g%isc,g%iec
2314  ! Restrict rms topo to 10 percent of column depth.
2315  h2(i,j) = min(0.01*(g%bathyT(i,j))**2, h2(i,j))
2316  ! Compute the fixed part; units are [kg m-2] here
2317  ! will be multiplied by N and En to get into [W m-2]
2318  cs%TKE_itidal_loss_fixed(i,j) = 0.5*kappa_h2_factor*gv%Rho0*&
2319  kappa_itides * h2(i,j)
2320  enddo ; enddo
2321 
2322  deallocate(h2)
2323 
2324  ! Read in prescribed coast/ridge/shelf angles from file
2325  call get_param(param_file, mdl, "REFL_ANGLE_FILE", refl_angle_file, &
2326  "The path to the file containing the local angle of "//&
2327  "the coastline/ridge/shelf with respect to the equator.", &
2328  fail_if_missing=.false.)
2329  filename = trim(cs%inputdir) // trim(refl_angle_file)
2330  call log_param(param_file, mdl, "INPUTDIR/REFL_ANGLE_FILE", filename)
2331  allocate(cs%refl_angle(isd:ied,jsd:jed)) ; cs%refl_angle(:,:) = cs%nullangle
2332  call mom_read_data(filename, 'refl_angle', cs%refl_angle, &
2333  g%domain, timelevel=1)
2334  ! replace NANs with null value
2335  do j=g%jsc,g%jec ; do i=g%isc,g%iec
2336  if (is_nan(cs%refl_angle(i,j))) cs%refl_angle(i,j) = cs%nullangle
2337  enddo ; enddo
2338  call pass_var(cs%refl_angle,g%domain)
2339 
2340  ! Read in prescribed partial reflection coefficients from file
2341  call get_param(param_file, mdl, "REFL_PREF_FILE", refl_pref_file, &
2342  "The path to the file containing the reflection coefficients.", &
2343  fail_if_missing=.false.)
2344  filename = trim(cs%inputdir) // trim(refl_pref_file)
2345  call log_param(param_file, mdl, "INPUTDIR/REFL_PREF_FILE", filename)
2346  allocate(cs%refl_pref(isd:ied,jsd:jed)) ; cs%refl_pref(:,:) = 1.0
2347  call mom_read_data(filename, 'refl_pref', cs%refl_pref, g%domain, timelevel=1)
2348  !CS%refl_pref = CS%refl_pref*1 ! adjust partial reflection if desired
2349  call pass_var(cs%refl_pref,g%domain)
2350 
2351  ! Tag reflection cells with partial reflection (done here for speed)
2352  allocate(cs%refl_pref_logical(isd:ied,jsd:jed)) ; cs%refl_pref_logical(:,:) = .false.
2353  do j=jsd,jed
2354  do i=isd,ied
2355  ! flag cells with partial reflection
2356  if (cs%refl_angle(i,j) /= cs%nullangle .and. &
2357  cs%refl_pref(i,j) < 1.0 .and. cs%refl_pref(i,j) > 0.0) then
2358  cs%refl_pref_logical(i,j) = .true.
2359  endif
2360  enddo
2361  enddo
2362 
2363  ! Read in double-reflective (ridge) tags from file
2364  call get_param(param_file, mdl, "REFL_DBL_FILE", refl_dbl_file, &
2365  "The path to the file containing the double-reflective ridge tags.", &
2366  fail_if_missing=.false.)
2367  filename = trim(cs%inputdir) // trim(refl_dbl_file)
2368  call log_param(param_file, mdl, "INPUTDIR/REFL_DBL_FILE", filename)
2369  allocate(ridge_temp(isd:ied,jsd:jed)) ; ridge_temp(:,:) = 0.0
2370  call mom_read_data(filename, 'refl_dbl', ridge_temp, g%domain, timelevel=1)
2371  call pass_var(ridge_temp,g%domain)
2372  allocate(cs%refl_dbl(isd:ied,jsd:jed)) ; cs%refl_dbl(:,:) = .false.
2373  do i=isd,ied; do j=jsd,jed
2374  if (ridge_temp(i,j) == 1) then; cs%refl_dbl(i,j) = .true.
2375  else ; cs%refl_dbl(i,j) = .false. ; endif
2376  enddo ; enddo
2377 
2378  ! Read in prescribed land mask from file (if overwriting -BDM).
2379  ! This should be done in MOM_initialize_topography subroutine
2380  ! defined in MOM_fixed_initialization.F90 (BDM)
2381  !call get_param(param_file, mdl, "LAND_MASK_FILE", land_mask_file, &
2382  ! "The path to the file containing the land mask.", &
2383  ! fail_if_missing=.false.)
2384  !filename = trim(CS%inputdir) // trim(land_mask_file)
2385  !call log_param(param_file, mdl, "INPUTDIR/LAND_MASK_FILE", filename)
2386  !G%mask2dCu(:,:) = 1 ; G%mask2dCv(:,:) = 1 ; G%mask2dT(:,:) = 1
2387  !call MOM_read_data(filename, 'land_mask', G%mask2dCu, G%domain, timelevel=1)
2388  !call MOM_read_data(filename, 'land_mask', G%mask2dCv, G%domain, timelevel=1)
2389  !call MOM_read_data(filename, 'land_mask', G%mask2dT, G%domain, timelevel=1)
2390  !call pass_var(G%mask2dCu,G%domain)
2391  !call pass_var(G%mask2dCv,G%domain)
2392  !call pass_var(G%mask2dT,G%domain)
2393 
2394  ! Read in prescribed partial east face blockages from file (if overwriting -BDM)
2395  !call get_param(param_file, mdl, "dy_Cu_FILE", dy_Cu_file, &
2396  ! "The path to the file containing the east face blockages.", &
2397  ! fail_if_missing=.false.)
2398  !filename = trim(CS%inputdir) // trim(dy_Cu_file)
2399  !call log_param(param_file, mdl, "INPUTDIR/dy_Cu_FILE", filename)
2400  !G%dy_Cu(:,:) = 0.0
2401  !call MOM_read_data(filename, 'dy_Cu', G%dy_Cu, G%domain, timelevel=1)
2402  !call pass_var(G%dy_Cu,G%domain)
2403 
2404  ! Read in prescribed partial north face blockages from file (if overwriting -BDM)
2405  !call get_param(param_file, mdl, "dx_Cv_FILE", dx_Cv_file, &
2406  ! "The path to the file containing the north face blockages.", &
2407  ! fail_if_missing=.false.)
2408  !filename = trim(CS%inputdir) // trim(dx_Cv_file)
2409  !call log_param(param_file, mdl, "INPUTDIR/dx_Cv_FILE", filename)
2410  !G%dx_Cv(:,:) = 0.0
2411  !call MOM_read_data(filename, 'dx_Cv', G%dx_Cv, G%domain, timelevel=1)
2412  !call pass_var(G%dx_Cv,G%domain)
2413 
2414  ! Register maps of reflection parameters
2415  cs%id_refl_ang = register_diag_field('ocean_model', 'refl_angle', diag%axesT1, &
2416  time, 'Local angle of coastline/ridge/shelf with respect to equator', 'rad')
2417  cs%id_refl_pref = register_diag_field('ocean_model', 'refl_pref', diag%axesT1, &
2418  time, 'Partial reflection coefficients', '')
2419  cs%id_dx_Cv = register_diag_field('ocean_model', 'dx_Cv', diag%axesT1, &
2420  time, 'North face unblocked width', 'm') ! used if overriding (BDM)
2421  cs%id_dy_Cu = register_diag_field('ocean_model', 'dy_Cu', diag%axesT1, &
2422  time, 'East face unblocked width', 'm') ! used if overriding (BDM)
2423  cs%id_land_mask = register_diag_field('ocean_model', 'land_mask', diag%axesT1, &
2424  time, 'Land mask', 'logical') ! used if overriding (BDM)
2425  ! Output reflection parameters as diags here (not needed every timestep)
2426  if (cs%id_refl_ang > 0) call post_data(cs%id_refl_ang, cs%refl_angle, cs%diag)
2427  if (cs%id_refl_pref > 0) call post_data(cs%id_refl_pref, cs%refl_pref, cs%diag)
2428  if (cs%id_dx_Cv > 0) call post_data(cs%id_dx_Cv, g%dx_Cv, cs%diag)
2429  if (cs%id_dy_Cu > 0) call post_data(cs%id_dy_Cu, g%dy_Cu, cs%diag)
2430  if (cs%id_land_mask > 0) call post_data(cs%id_land_mask, g%mask2dT, cs%diag)
2431 
2432  ! Register 2-D energy density (summed over angles, freq, modes)
2433  cs%id_tot_En = register_diag_field('ocean_model', 'ITide_tot_En', diag%axesT1, &
2434  time, 'Internal tide total energy density', 'J m-2')
2435  ! Register 2-D drag scale used for quadratic bottom drag
2436  cs%id_itide_drag = register_diag_field('ocean_model', 'ITide_drag', diag%axesT1, &
2437  time, 'Interior and bottom drag internal tide decay timescale', 's-1')
2438  !Register 2-D energy input into internal tides
2439  cs%id_TKE_itidal_input = register_diag_field('ocean_model', 'TKE_itidal_input', diag%axesT1, &
2440  time, 'Conversion from barotropic to baroclinic tide, '//&
2441  'a fraction of which goes into rays', 'W m-2')
2442  ! Register 2-D energy losses (summed over angles, freq, modes)
2443  cs%id_tot_leak_loss = register_diag_field('ocean_model', 'ITide_tot_leak_loss', diag%axesT1, &
2444  time, 'Internal tide energy loss to background drag', 'W m-2')
2445  cs%id_tot_quad_loss = register_diag_field('ocean_model', 'ITide_tot_quad_loss', diag%axesT1, &
2446  time, 'Internal tide energy loss to bottom drag', 'W m-2')
2447  cs%id_tot_itidal_loss = register_diag_field('ocean_model', 'ITide_tot_itidal_loss', diag%axesT1, &
2448  time, 'Internal tide energy loss to wave drag', 'W m-2')
2449  cs%id_tot_Froude_loss = register_diag_field('ocean_model', 'ITide_tot_Froude_loss', diag%axesT1, &
2450  time, 'Internal tide energy loss to wave breaking', 'W m-2')
2451  cs%id_tot_allprocesses_loss = register_diag_field('ocean_model', 'ITide_tot_allprocesses_loss', diag%axesT1, &
2452  time, 'Internal tide energy loss summed over all processes', 'W m-2')
2453 
2454  allocate(cs%id_En_mode(cs%nFreq,cs%nMode)) ; cs%id_En_mode(:,:) = -1
2455  allocate(cs%id_En_ang_mode(cs%nFreq,cs%nMode)) ; cs%id_En_ang_mode(:,:) = -1
2456  allocate(cs%id_itidal_loss_mode(cs%nFreq,cs%nMode)) ; cs%id_itidal_loss_mode(:,:) = -1
2457  allocate(cs%id_allprocesses_loss_mode(cs%nFreq,cs%nMode)) ; cs%id_allprocesses_loss_mode(:,:) = -1
2458  allocate(cs%id_itidal_loss_ang_mode(cs%nFreq,cs%nMode)) ; cs%id_itidal_loss_ang_mode(:,:) = -1
2459  allocate(cs%id_Ub_mode(cs%nFreq,cs%nMode)) ; cs%id_Ub_mode(:,:) = -1
2460  allocate(cs%id_cp_mode(cs%nFreq,cs%nMode)) ; cs%id_cp_mode(:,:) = -1
2461 
2462  allocate(angles(cs%NAngle)) ; angles(:) = 0.0
2463  angle_size = (8.0*atan(1.0)) / (real(num_angle))
2464  do a=1,num_angle ; angles(a) = (real(a) - 1) * angle_size ; enddo
2465 
2466  id_ang = diag_axis_init("angle", angles, "Radians", "N", "Angular Orienation of Fluxes")
2467  call define_axes_group(diag, (/ diag%axesT1%handles(1), diag%axesT1%handles(2), id_ang /), axes_ang)
2468 
2469  do fr=1,cs%nFreq ; write(freq_name(fr), '("freq",i1)') fr ; enddo
2470  do m=1,cs%nMode ; do fr=1,cs%nFreq
2471  ! Register 2-D energy density (summed over angles) for each freq and mode
2472  write(var_name, '("Itide_En_freq",i1,"_mode",i1)') fr, m
2473  write(var_descript, '("Internal tide energy density in frequency ",i1," mode ",i1)') fr, m
2474  cs%id_En_mode(fr,m) = register_diag_field('ocean_model', var_name, &
2475  diag%axesT1, time, var_descript, 'J m-2')
2476  call mom_mesg("Registering "//trim(var_name)//", Described as: "//var_descript, 5)
2477 
2478  ! Register 3-D (i,j,a) energy density for each freq and mode
2479  write(var_name, '("Itide_En_ang_freq",i1,"_mode",i1)') fr, m
2480  write(var_descript, '("Internal tide angular energy density in frequency ",i1," mode ",i1)') fr, m
2481  cs%id_En_ang_mode(fr,m) = register_diag_field('ocean_model', var_name, &
2482  axes_ang, time, var_descript, 'J m-2 band-1')
2483  call mom_mesg("Registering "//trim(var_name)//", Described as: "//var_descript, 5)
2484 
2485  ! Register 2-D energy loss (summed over angles) for each freq and mode
2486  ! wave-drag only
2487  write(var_name, '("Itide_wavedrag_loss_freq",i1,"_mode",i1)') fr, m
2488  write(var_descript, '("Internal tide energy loss due to wave-drag from frequency ",i1," mode ",i1)') fr, m
2489  cs%id_itidal_loss_mode(fr,m) = register_diag_field('ocean_model', var_name, &
2490  diag%axesT1, time, var_descript, 'W m-2')
2491  call mom_mesg("Registering "//trim(var_name)//", Described as: "//var_descript, 5)
2492  ! all loss processes
2493  write(var_name, '("Itide_allprocesses_loss_freq",i1,"_mode",i1)') fr, m
2494  write(var_descript, '("Internal tide energy loss due to all processes from frequency ",i1," mode ",i1)') fr, m
2495  cs%id_allprocesses_loss_mode(fr,m) = register_diag_field('ocean_model', var_name, &
2496  diag%axesT1, time, var_descript, 'W m-2')
2497  call mom_mesg("Registering "//trim(var_name)//", Described as: "//var_descript, 5)
2498 
2499  ! Register 3-D (i,j,a) energy loss for each freq and mode
2500  ! wave-drag only
2501  write(var_name, '("Itide_wavedrag_loss_ang_freq",i1,"_mode",i1)') fr, m
2502  write(var_descript, '("Internal tide energy loss due to wave-drag from frequency ",i1," mode ",i1)') fr, m
2503  cs%id_itidal_loss_ang_mode(fr,m) = register_diag_field('ocean_model', var_name, &
2504  axes_ang, time, var_descript, 'W m-2 band-1')
2505  call mom_mesg("Registering "//trim(var_name)//", Described as: "//var_descript, 5)
2506 
2507  ! Register 2-D period-averaged near-bottom horizonal velocity for each freq and mode
2508  write(var_name, '("Itide_Ub_freq",i1,"_mode",i1)') fr, m
2509  write(var_descript, '("Near-bottom horizonal velocity for frequency ",i1," mode ",i1)') fr, m
2510  cs%id_Ub_mode(fr,m) = register_diag_field('ocean_model', var_name, &
2511  diag%axesT1, time, var_descript, 'm s-1')
2512  call mom_mesg("Registering "//trim(var_name)//", Described as: "//var_descript, 5)
2513 
2514  ! Register 2-D horizonal phase velocity for each freq and mode
2515  write(var_name, '("Itide_cp_freq",i1,"_mode",i1)') fr, m
2516  write(var_descript, '("Horizonal phase velocity for frequency ",i1," mode ",i1)') fr, m
2517  cs%id_cp_mode(fr,m) = register_diag_field('ocean_model', var_name, &
2518  diag%axesT1, time, var_descript, 'm s-1')
2519  call mom_mesg("Registering "//trim(var_name)//", Described as: "//var_descript, 5)
2520 
2521  enddo ; enddo
2522 
2523  ! Initialize wave_structure (not sure if this should be here - BDM)
2524  call wave_structure_init(time, g, param_file, diag, cs%wave_structure_CSp)
2525 
2526 end subroutine internal_tides_init
2527 
2528 !> This subroutine deallocates the memory associated with the internal tides control structure
2529 subroutine internal_tides_end(CS)
2530  type(int_tide_cs), pointer :: cs !< A pointer to the control structure returned by a previous
2531  !! call to internal_tides_init, it will be deallocated here.
2532 
2533  if (associated(cs)) then
2534  if (associated(cs%En)) deallocate(cs%En)
2535  if (allocated(cs%frequency)) deallocate(cs%frequency)
2536  if (allocated(cs%id_En_mode)) deallocate(cs%id_En_mode)
2537  if (allocated(cs%id_Ub_mode)) deallocate(cs%id_Ub_mode)
2538  if (allocated(cs%id_cp_mode)) deallocate(cs%id_cp_mode)
2539  deallocate(cs)
2540  endif
2541  cs => null()
2542 end subroutine internal_tides_end
2543 
2544 end module mom_internal_tides
mom_time_manager
Wraps the FMS time manager functions.
Definition: MOM_time_manager.F90:2
mom_variables::surface
Pointers to various fields which may be used describe the surface state of MOM, and which will be ret...
Definition: MOM_variables.F90:38
mom_spatial_means
Functions and routines to take area, volume, mass-weighted, layerwise, zonal or meridional means.
Definition: MOM_spatial_means.F90:2
mom_verticalgrid
Provides a transparent vertical ocean grid type and supporting routines.
Definition: MOM_verticalGrid.F90:2
mom_file_parser::log_version
An overloaded interface to log version information about modules.
Definition: MOM_file_parser.F90:109
mom_diag_mediator
The subroutines here provide convenient wrappers to the fms diag_manager interfaces with additional d...
Definition: MOM_diag_mediator.F90:3
mom_variables::thermo_var_ptrs
Pointers to an assortment of thermodynamic fields that may be available, including potential temperat...
Definition: MOM_variables.F90:82
mom_diag_mediator::axes_grp
A group of 1D axes that comprise a 1D/2D/3D mesh.
Definition: MOM_diag_mediator.F90:96
mom_domains::pass_var
Do a halo update on an array.
Definition: MOM_domains.F90:49
mom_file_parser::param_file_type
A structure that can be parsed to read and document run-time parameters.
Definition: MOM_file_parser.F90:54
mom_file_parser::get_param
An overloaded interface to read and log the values of various types of parameters.
Definition: MOM_file_parser.F90:102
mom_wave_structure
Vertical structure functions for first baroclinic mode wave speed.
Definition: MOM_wave_structure.F90:2
mom_io
This module contains I/O framework code.
Definition: MOM_io.F90:2
mom_restart::mom_restart_cs
A restart registry and the control structure for restarts.
Definition: MOM_restart.F90:72
mom_unit_scaling::unit_scale_type
Describes various unit conversion factors.
Definition: MOM_unit_scaling.F90:14
mom_internal_tides::int_tide_cs
This control structure has parameters for the MOM_internal_tides module.
Definition: MOM_internal_tides.F90:38
mom_diag_mediator::post_data
Make a diagnostic available for averaging or output.
Definition: MOM_diag_mediator.F90:70
mom_verticalgrid::verticalgrid_type
Describes the vertical ocean grid, including unit conversion factors.
Definition: MOM_verticalGrid.F90:24
mom_restart
The MOM6 facility for reading and writing restart files, and querying what has been read.
Definition: MOM_restart.F90:2
mom_domains
Describes the decomposed MOM domain and has routines for communications across PEs.
Definition: MOM_domains.F90:2
mom_variables
Provides transparent structures with groups of MOM6 variables and supporting routines.
Definition: MOM_variables.F90:2
mom_io::mom_read_data
Read a data field from a file.
Definition: MOM_io.F90:74
mom_wave_structure::wave_structure_cs
The control structure for the MOM_wave_structure module.
Definition: MOM_wave_structure.F90:36
mom_file_parser
The MOM6 facility to parse input files for runtime parameters.
Definition: MOM_file_parser.F90:2
mom_grid
Provides the ocean grid type.
Definition: MOM_grid.F90:2
mom_restart::register_restart_field
Register fields for restarts.
Definition: MOM_restart.F90:107
mom_unit_scaling
Provides a transparent unit rescaling type to facilitate dimensional consistency testing.
Definition: MOM_unit_scaling.F90:2
mom_debugging
Provides checksumming functions for debugging.
Definition: MOM_debugging.F90:7
mom_io::vardesc
Type for describing a variable, typically a tracer.
Definition: MOM_io.F90:53
mom_internal_tides
Subroutines that use the ray-tracing equations to propagate the internal tide energy density.
Definition: MOM_internal_tides.F90:4
mom_file_parser::log_param
An overloaded interface to log the values of various types of parameters.
Definition: MOM_file_parser.F90:96
mom_internal_tides::loop_bounds_type
A structure with the active energy loop bounds.
Definition: MOM_internal_tides.F90:140
mom_domains::create_group_pass
Set up a group of halo updates.
Definition: MOM_domains.F90:79
mom_error_handler
Routines for error handling and I/O management.
Definition: MOM_error_handler.F90:2
mom_grid::ocean_grid_type
Ocean grid type. See mom_grid for details.
Definition: MOM_grid.F90:25
mom_diag_mediator::diag_ctrl
The following data type a list of diagnostic fields an their variants, as well as variables that cont...
Definition: MOM_diag_mediator.F90:239
mom_file_parser::read_param
An overloaded interface to read various types of parameters.
Definition: MOM_file_parser.F90:90