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