MOM_internal_tides.F90

Go to the documentation of this file.

module mom_internal_tides
!***********************************************************************
!*                   GNU General Public License                        *
!* This file is a part of MOM.                                         *
!*                                                                     *
!* MOM is free software; you can redistribute it and/or modify it and  *
!* are expected to follow the terms of the GNU General Public License  *
!* as published by the Free Software Foundation; either version 2 of   *
!* the License, or (at your option) any later version.                 *
!*                                                                     *
!* MOM is distributed in the hope that it will be useful, but WITHOUT  *
!* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY  *
!* or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public    *
!* License for more details.                                           *
!*                                                                     *
!* For the full text of the GNU General Public License,                *
!* write to: Free Software Foundation, Inc.,                           *
!*           675 Mass Ave, Cambridge, MA 02139, USA.                   *
!* or see:   http://www.gnu.org/licenses/gpl.html                      *
!***********************************************************************
!
!********+*********+*********+*********+*********+*********+*********+**
!*                                                                     *
!*  By Benjamin Mater & Robert Hallberg, 2015                          *
!*                                                                     *
!*    This program contains the subroutines that use the ray-tracing   *
!*  equations to propagate the internal tide energy density.           *
!*                                                                     *
!*  Macros written all in capital letters are defined in MOM_memory.h. *
!*                                                                     *
!*     A small fragment of the grid is shown below:                    *
!*                                                                     *
!*    j+1  x ^ x ^ x   At x:  q                                        *
!*    j+1  > o > o >   At ^:  v, tauy                                  *
!*    j    x ^ x ^ x   At >:  u, taux                                  *
!*    j    > o > o >   At o:  h, fluxes.                               *
!*    j-1  x ^ x ^ x                                                   *
!*        i-1  i  i+1  At x & ^:                                       *
!*           i  i+1    At > & o:                                       *
!*                                                                     *
!*  The boundaries always run through q grid points (x).               *
!*                                                                     *
!********+*********+*********+*********+*********+*********+*********+**

use mom_debugging,     only : is_nan
use mom_diag_mediator, only : post_data, query_averaging_enabled, diag_axis_init
use mom_diag_mediator, only : register_diag_field, diag_ctrl, safe_alloc_ptr
use mom_diag_mediator, only : axes_grp, define_axes_group
use mom_domains, only       : agrid, to_south, to_west, to_all
use mom_domains, only       : create_group_pass, do_group_pass, pass_var
use mom_domains, only       : group_pass_type, start_group_pass, complete_group_pass
use mom_error_handler, only : mom_error, fatal, warning, mom_mesg, is_root_pe
use mom_file_parser, only   : read_param, get_param, log_param, log_version, param_file_type
use mom_grid, only          : ocean_grid_type
use mom_io, only            : slasher, vardesc
use mom_restart, only       : register_restart_field, mom_restart_cs, restart_init, save_restart
use mom_spatial_means, only : global_area_mean
use mom_time_manager, only  : time_type, operator(+), operator(/), operator(-)
use mom_time_manager, only  : get_time, get_date, set_time, set_date
use mom_time_manager, only  : time_type_to_real
use mom_variables, only     : surface, thermo_var_ptrs
use mom_verticalgrid, only  : verticalgrid_type
use fms_mod, only           : read_data
use mom_wave_structure, only: wave_structure_init, wave_structure, wave_structure_cs

!   Forcing is a structure containing pointers to the forcing fields
! which may be used to drive MOM.  All fluxes are positive downward.
!   Surface is a structure containing pointers to various fields that
! may be used describe the surface state of MOM.

!use, intrinsic :: IEEE_ARITHMETIC

implicit none ; private

#include <MOM_memory.h>

public propagate_int_tide !, register_int_tide_restarts
public internal_tides_init, internal_tides_end
public get_lowmode_loss

type, public :: int_tide_cs ; private
  logical :: do_int_tides    ! If true, use the internal tide code.
  integer :: nfreq = 0
  integer :: nmode = 1
  integer :: nangle = 24
  integer :: energized_angle = -1
  logical :: corner_adv      ! If true, use a corner advection rather than PPM.
  logical :: upwind_1st      ! If true, use a first-order upwind scheme.
  logical :: simple_2nd      ! If true, use a simple second order (arithmetic
                             ! mean) interpolation of the edge values instead
                             ! of the higher order interpolation.
  logical :: vol_cfl         ! If true, use the ratio of the open face lengths
                             ! to the tracer cell areas when estimating CFL
                             ! numbers.  Without aggress_adjust, the default is
                             ! false; it is always true with.
  logical :: use_ppmang      ! If true, use PPM for advection of energy in
                             ! angular space.

  real, allocatable, dimension(:,:) :: refl_angle
                        ! local coastline/ridge/shelf angles read from file
                        ! (could be in G control structure)
  real :: nullangle = -999.9 ! placeholder value in cell with no reflection
  real, allocatable, dimension(:,:) :: refl_pref
                        ! partial reflection coeff for each ``coast cell"
                        ! (could be in G control structure)
  logical, allocatable, dimension(:,:) :: refl_pref_logical
                        ! true if reflecting cell with partial reflection
                        ! (could be in G control structure)
  logical, allocatable, dimension(:,:) :: refl_dbl
                        ! identifies reflection cells where double reflection
                        ! is possible (i.e. ridge cells)
                        ! (could be in G control structure)
  real, allocatable, dimension(:,:,:,:) :: cp
                        ! horizontal phase speed [m s-1]
  real, allocatable, dimension(:,:,:,:,:) :: tke_leak_loss
                        ! energy lost due to misc background processes [W m-2]
  real, allocatable, dimension(:,:,:,:,:) :: tke_quad_loss
                        ! energy lost due to quadratic bottom drag [W m-2]
  real, allocatable, dimension(:,:,:,:,:) :: tke_froude_loss
                        ! energy lost due to wave breaking [W m-2]
  real, allocatable, dimension(:,:) :: tke_itidal_loss_fixed
                        ! fixed part of the energy lost due to small-scale drag
                        ! [kg m-2] here; will be multiplied by N and En to get
                        ! into [W m-2]
  real, allocatable, dimension(:,:,:,:,:) :: tke_itidal_loss
                        ! energy lost due to small-scale wave drag [W m-2]
  real, allocatable, dimension(:,:) :: tot_leak_loss, tot_quad_loss, &
                                       tot_itidal_loss, tot_froude_loss, tot_allprocesses_loss
                        ! energy loss rates summed over angle, freq, and mode
  real :: q_itides      ! fraction of local dissipation (nondimensional)
  real :: en_sum        ! global sum of energy for use in debugging
  type(time_type),pointer    :: time
                        ! The current model time
  character(len=200) :: inputdir
                        ! directory to look for coastline angle file
  real :: decay_rate    ! A constant rate at which internal tide energy is
                        ! lost to the interior ocean internal wave field.
  real :: cdrag         ! The bottom drag coefficient for MEKE (non-dim).
  logical :: apply_background_drag
                        ! If true, apply a drag due to background processes as a sink.
  logical :: apply_bottom_drag
                        ! If true, apply a quadratic bottom drag as a sink.
  logical :: apply_wave_drag
                        ! If true, apply scattering due to small-scale
                        ! roughness as a sink.
  logical :: apply_froude_drag
                        ! If true, apply wave breaking as a sink.
  real, dimension(:,:,:,:,:), pointer :: &
    en                  ! The internal wave energy density as a function of
                        ! (i,j,angle,frequency,mode)
  real, dimension(:,:,:), pointer :: &
    en_restart          ! The internal wave energy density as a function of
                        ! (i,j,angle); temporary for restart
  real, allocatable, dimension(:) :: &
    frequency           ! The frequency of each band.

  real    :: int_tide_source_x ! delete later
                               ! X Location of generation site
                               ! for internal tide for testing
  real    :: int_tide_source_y ! delete later
                               ! Y Location of generation site
                               ! for internal tide for testing

  type(diag_ctrl), pointer :: diag ! A structure that is used to regulate the
                        ! timing of diagnostic output.
  type(wave_structure_cs),  pointer :: wave_structure_csp => null()

  ! Diag handles relevant to all modes, frequencies, and angles
  integer :: id_itide_drag = -1
  integer :: id_refl_pref = -1, id_refl_ang = -1, id_land_mask = -1
  integer :: id_dx_cv = -1, id_dy_cu = -1
  integer :: id_tke_itidal_input = -1
  ! Diag handles considering: sums over all modes, frequencies, and angles
  integer :: id_tot_en = -1, &
             id_tot_leak_loss = -1, &
             id_tot_quad_loss = -1, &
             id_tot_itidal_loss = -1, &
             id_tot_froude_loss = -1, &
             id_tot_allprocesses_loss = -1
  ! Diag handles considering: all modes & freqs; summed over angles
  integer, allocatable, dimension(:,:) :: &
             id_en_mode, &
             id_itidal_loss_mode, &
             id_allprocesses_loss_mode, &
             id_ub_mode, &
             id_cp_mode
  ! Diag handles considering: all modes, freqs, and angles
  integer, allocatable, dimension(:,:) :: &
             id_en_ang_mode, &
             id_itidal_loss_ang_mode

end type int_tide_cs

type :: loop_bounds_type ; private
  integer :: ish, ieh, jsh, jeh
end type loop_bounds_type

contains

!> This subroutine calls other subroutines in this file that are needed to
!! refract, propagate, and dissipate energy density of the internal tide.
subroutine propagate_int_tide(h, tv, cn, TKE_itidal_input, vel_btTide, Nb, dt, &
                              G, GV, CS)
  type(ocean_grid_type),            intent(inout) :: G  !< The ocean's grid structure.
  type(verticalgrid_type),          intent(in)    :: GV !< The ocean's vertical grid structure.
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
                                    intent(in)    :: h  !< Layer thicknesses, in H
                                                        !! (usually m or kg m-2).
  type(thermo_var_ptrs),            intent(in)    :: tv !< Pointer to thermodynamic variables
                                                        !! (needed for wave structure).
  real, dimension(SZI_(G),SZJ_(G)), intent(in)    :: TKE_itidal_input !< The energy input to the
                                                        !! internal waves, in W m-2.
  real, dimension(SZI_(G),SZJ_(G)), intent(in)    :: vel_btTide       !< Barotropic velocity read
                                                        !! from file, in m s-1.
  real, dimension(SZI_(G),SZJ_(G)), intent(in)    :: Nb !< Near-bottom buoyancy frequency, in s-1.
  real,                             intent(in)    :: dt !< Length of time over which these fluxes
                                                        !! will be applied, in s.
  type(int_tide_cs),                pointer       :: CS !< A pointer to the control structure
                                                        !! returned by a previous call to
                                                        !! int_tide_init.
  real, dimension(SZI_(G),SZJ_(G),CS%nMode), &
                                    intent(in)    :: cn

  ! This subroutine calls other subroutines in this file that are needed to
  ! refract, propagate, and dissipate energy density of the internal tide.
  !
  ! Arguments:
  !  (in)      h -  Layer thickness, in m or kg m-2  (needed for wave structure).
  !  (in)      tv - Pointer to thermodynamic variables (needed for wave structure).
  !  (in)      cn - Internal gravity wave speeds of modes, in m s-1.
  !  (in)      TKE_itidal_input - The energy input to the internal waves, in W m-2.
  !  (in)      vel_btTide - Barotropic velocity read from file, in m s-1
  !  (in)      Nb - Near-bottom buoyancy frequency, in s-1
  !  (in)      dt - Length of time over which these fluxes will be applied, in s.
  !  (in)      G - The ocean's grid structure.
  !  (in)      GV - The ocean's vertical grid structure.
  !  (in)      CS - A pointer to the control structure returned by a previous
  !                 call to int_tide_init.
  real, dimension(SZI_(G),SZJ_(G),2) :: &
    test
  real, dimension(SZI_(G),SZJ_(G),CS%nFreq,CS%nMode) :: &
    tot_En_mode, & ! energy summed over angles only
    Ub, Umax       ! near-bottom & max horizontal velocity of wave (modal)
  real, dimension(SZI_(G),SZJB_(G)) :: &
    flux_heat_y, &
    flux_prec_y
  real, dimension(SZI_(G),SZJ_(G)) :: &
    tot_En, &      ! energy summed over angles, modes, frequencies
    tot_leak_loss, tot_quad_loss, tot_itidal_loss, tot_Froude_loss, tot_allprocesses_loss, &
                   ! energy loss rates summed over angle, freq, and mode
    drag_scale, &  ! bottom drag scale, s-1
    itidal_loss_mode, allprocesses_loss_mode
                   ! energy loss rates for a given mode and frequency (summed over angles)
  real :: frac_per_sector, f2, I_rho0, I_D_here, freq2, Kmag2
  real :: c_phase, loss_rate, Fr2_max
  real, parameter :: cn_subRO = 1e-100               ! to prevent division by zero
  real :: En_new, En_check                           ! for debugging
  real :: En_initial, Delta_E_check                  ! for debugging
  real :: TKE_Froude_loss_check, TKE_Froude_loss_tot ! for debugging
  integer :: a, m, fr, i, j, is, ie, js, je, isd, ied, jsd, jed, nAngle, nzm
  integer :: id_g, jd_g         ! global (decomp-invar) indices (for debugging)
  type(group_pass_type), save :: pass_test, pass_En

  if (.not.associated(cs)) return
  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
  isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed ; nangle = cs%NAngle
  i_rho0 = 1.0 / gv%Rho0

  ! Set the wave speeds for the modes, using cg(n) ~ cg(1)/n.**********************
  ! This is wrong, of course, but it works reasonably in some cases.
  ! Uncomment if wave_speed is not used to calculate the true values (BDM).
  !do m=1,CS%nMode ; do j=jsd,jed ; do i=isd,ied
  !  cn(i,j,m) = cn(i,j,1) / real(m)
  !enddo ; enddo ; enddo

  ! Add the forcing.***************************************************************
  if (cs%energized_angle <= 0) then
    frac_per_sector = 1.0 / real(cs%nangle * cs%nmode * cs%nfreq)
    do m=1,cs%nMode ; do fr=1,cs%nFreq ; do a=1,cs%nAngle ; do j=js,je ; do i=is,ie
      f2 = 0.25*((g%CoriolisBu(i,j)**2 + g%CoriolisBu(i-1,j-1)**2) + &
                 (g%CoriolisBu(i,j-1)**2 + g%CoriolisBu(i-1,j)**2))
      if (cs%frequency(fr)**2 > f2) &
        cs%En(i,j,a,fr,m) = cs%En(i,j,a,fr,m) + &
                            dt*frac_per_sector*(1-cs%q_itides)*tke_itidal_input(i,j)
    enddo ; enddo ; enddo ; enddo ; enddo
  elseif (cs%energized_angle <= cs%nAngle) then
    frac_per_sector = 1.0 / real(cs%nmode * cs%nfreq)
    a = cs%energized_angle
    do m=1,cs%nMode ; do fr=1,cs%nFreq ; do j=js,je ; do i=is,ie
      f2 = 0.25*((g%CoriolisBu(i,j)**2 + g%CoriolisBu(i-1,j-1)**2) + &
                 (g%CoriolisBu(i,j-1)**2 + g%CoriolisBu(i-1,j)**2))
      if (cs%frequency(fr)**2 > f2) &
        cs%En(i,j,a,fr,m) = cs%En(i,j,a,fr,m) + &
                            dt*frac_per_sector**(1-cs%q_itides)*tke_itidal_input(i,j)
    enddo ; enddo ; enddo ; enddo
  else
    call mom_error(warning, "Internal tide energy is being put into a angular "//&
                            "band that does not exist.")
  endif

  ! Pass a test vector to check for grid rotation in the halo updates.
  do j=jsd,jed ; do i=isd,ied ; test(i,j,1) = 1.0 ; test(i,j,2) = 0.0 ; enddo ; enddo
  do m=1,cs%nMode ; do fr=1,cs%nFreq
    call create_group_pass(pass_en, cs%En(:,:,:,fr,m), g%domain)
  enddo; enddo
  call create_group_pass(pass_test, test(:,:,1), test(:,:,2), g%domain, stagger=agrid)
  call start_group_pass(pass_test, g%domain)

  ! Apply half the refraction.*****************************************************
  do m=1,cs%nMode ; do fr=1,cs%nFreq
    call refract(cs%En(:,:,:,fr,m), cn(:,:,m), cs%frequency(fr), 0.5*dt, g, cs%nAngle, cs%use_PPMang)
  enddo ; enddo

  ! Check for En<0 - for debugging, delete later
  do m=1,cs%NMode ; do fr=1,cs%Nfreq ; do a=1,cs%nAngle
    do j=js,je ; do i=is,ie
      if (cs%En(i,j,a,fr,m)<0.0) then
        id_g = i + g%idg_offset ; jd_g = j + g%jdg_offset ! for debugging
        print *, 'After first refraction: En<0.0 at ig=', id_g, ', jg=', jd_g
        print *, 'En=',cs%En(i,j,a,fr,m)
        print *, 'Setting En to zero'; cs%En(i,j,a,fr,m) = 0.0
        !stop
      endif
    enddo ; enddo
  enddo ; enddo ; enddo

  call do_group_pass(pass_en, g%domain)

  call complete_group_pass(pass_test, g%domain)

  ! Rotate points in the halos as necessary.
  call correct_halo_rotation(cs%En, test, g, cs%nAngle)

  ! Propagate the waves.***********************************************************
  do m=1,cs%NMode ; do fr=1,cs%Nfreq
    call propagate(cs%En(:,:,:,fr,m), cn(:,:,m), cs%frequency(fr), dt, g, cs, cs%NAngle)
  enddo ; enddo

  ! Check for En<0 - for debugging, delete later
  do m=1,cs%NMode ; do fr=1,cs%Nfreq ; do a=1,cs%nAngle
    do j=js,je ; do i=is,ie
      if (cs%En(i,j,a,fr,m)<0.0) then
        id_g = i + g%idg_offset ; jd_g = j + g%jdg_offset
        cs%En(i,j,a,fr,m) = 0.0
        if(abs(cs%En(i,j,a,fr,m))>1.0)then! only print if large
          print *, 'After propagation: En<0.0 at ig=', id_g, ', jg=', jd_g
          print *, 'En=',cs%En(i,j,a,fr,m)
          print *, 'Setting En to zero'
          !stop
        endif
      endif
    enddo ; enddo
  enddo ; enddo ; enddo

  !! Test if energy has passed coast for debugging only; delete later
  !do j=js,je
  !  do i=is,ie
  !    id_g = i + G%idg_offset ; jd_g = j + G%jdg_offset
  !    if (id_g == 106 .and. jd_g == 55 ) then
        !print *, 'After propagation:'
        !print *, 'En_O  =', CS%En(i,j,:,1,1), 'refl_angle=', CS%refl_angle(i,j)
        !print *, 'En_W =', CS%En(i-1,j,:,1,1), 'refl_angle=', CS%refl_angle(i-1,j)
        !print *, 'En_NW =', CS%En(i-1,j+1,:,1,1), 'refl_angle=', CS%refl_angle(i-1,j+1)
        !print *, 'En_N =', CS%En(i,j+1,:,1,1), 'refl_angle=', CS%refl_angle(i,j+1)
        !print *, 'En_NE =', CS%En(i+1,j+1,:,1,1), 'refl_angle=', CS%refl_angle(i+1,j+1)
        !print *, 'En_E =', CS%En(i+1,j,:,1,1), 'refl_angle=', CS%refl_angle(i+1,j)
 !     endif
 !   enddo
 ! enddo

  ! Apply the other half of the refraction.****************************************
  do m=1,cs%NMode ; do fr=1,cs%Nfreq
    call refract(cs%En(:,:,:,fr,m), cn(:,:,m), cs%frequency(fr), 0.5*dt, g, cs%NAngle, cs%use_PPMang)
  enddo ; enddo

  ! Check for En<0 - for debugging, delete later
  do m=1,cs%NMode ; do fr=1,cs%Nfreq ; do a=1,cs%nAngle
    do j=js,je ; do i=is,ie
      if (cs%En(i,j,a,fr,m)<0.0) then
        id_g = i + g%idg_offset ; jd_g = j + g%jdg_offset ! for debugging
        print *, 'After second refraction: En<0.0 at ig=', id_g, ', jg=', jd_g
        !stop
      endif
    enddo ; enddo
  enddo ; enddo ; enddo

  ! Apply various dissipation mechanisms.******************************************
  if (cs%apply_background_drag .or. cs%apply_bottom_drag &
      .or. cs%apply_wave_drag .or. cs%apply_Froude_drag &
      .or. (cs%id_tot_En > 0)) then
    tot_en(:,:) = 0.0
    tot_en_mode(:,:,:,:) = 0.0
    do m=1,cs%NMode ; do fr=1,cs%Nfreq
      do j=jsd,jed ; do i=isd,ied ; do a=1,cs%nAngle
        tot_en(i,j) = tot_en(i,j) + cs%En(i,j,a,fr,m)
        tot_en_mode(i,j,fr,m) = tot_en_mode(i,j,fr,m) + cs%En(i,j,a,fr,m)
      enddo ; enddo ; enddo
    enddo ; enddo
  endif

  ! Extract the energy for mixing due to misc. processes (background leakage)------
  if (cs%apply_background_drag) then
    do m=1,cs%nMode ; do fr=1,cs%nFreq ; do a=1,cs%nAngle ; do j=jsd,jed ; do i=isd,ied
    ! Calculate loss rate and apply loss over the time step ; apply the same drag timescale
    ! to each En component (technically not correct; fix later)
    cs%TKE_leak_loss(i,j,a,fr,m)  = cs%En(i,j,a,fr,m) * cs%decay_rate ! loss rate [Wm-2]
    cs%En(i,j,a,fr,m) = cs%En(i,j,a,fr,m) / (1.0 + dt *cs%decay_rate) ! implicit update
    enddo ; enddo ; enddo ; enddo ; enddo
  endif
  ! Check for En<0 - for debugging, delete later
  do m=1,cs%NMode ; do fr=1,cs%Nfreq ; do a=1,cs%nAngle
    do j=js,je ; do i=is,ie
      if (cs%En(i,j,a,fr,m)<0.0) then
        id_g = i + g%idg_offset ; jd_g = j + g%jdg_offset ! for debugging
        print *, 'After leak loss: En<0.0 at ig=', id_g, ', jg=', jd_g
        !stop
      endif
    enddo ; enddo
  enddo ; enddo ; enddo

  ! Extract the energy for mixing due to bottom drag-------------------------------
  if (cs%apply_bottom_drag) then
    do j=jsd,jed ; do i=isd,ied
      i_d_here = 1.0 / max(g%bathyT(i,j), 1.0)
      drag_scale(i,j) = cs%cdrag * sqrt(max(0.0, vel_bttide(i,j)**2 + &
                        tot_en(i,j) * i_rho0 * i_d_here)) * i_d_here
    enddo ; enddo
    do m=1,cs%nMode ; do fr=1,cs%nFreq ; do a=1,cs%nAngle ; do j=jsd,jed ; do i=isd,ied
      ! Calculate loss rate and apply loss over the time step ; apply the same drag timescale
      ! to each En component (technically not correct; fix later)
      cs%TKE_quad_loss(i,j,a,fr,m)  = cs%En(i,j,a,fr,m) * drag_scale(i,j) ! loss rate
      cs%En(i,j,a,fr,m) = cs%En(i,j,a,fr,m) / (1.0 + dt *drag_scale(i,j)) ! implicit update
    enddo ; enddo ; enddo ; enddo ; enddo
  endif
  ! Check for En<0 - for debugging, delete later
  do m=1,cs%NMode ; do fr=1,cs%Nfreq ; do a=1,cs%nAngle
    do j=js,je ; do i=is,ie
      if (cs%En(i,j,a,fr,m)<0.0) then
        id_g = i + g%idg_offset ; jd_g = j + g%jdg_offset ! for debugging
        print *, 'After bottom loss: En<0.0 at ig=', id_g, ', jg=', jd_g
        !stop
      endif
    enddo ; enddo
  enddo ; enddo ; enddo

  ! Extract the energy for mixing due to scattering (wave-drag)--------------------
  ! still need to allow a portion of the extracted energy to go to higher modes.
  ! First, find velocity profiles
  if (cs%apply_wave_drag .or. cs%apply_Froude_drag) then
    do m=1,cs%NMode ; do fr=1,cs%Nfreq
      ! Calculate modal structure for given mode and frequency
      call wave_structure(h, tv, g, gv, cn(:,:,m), m, cs%frequency(fr), &
                          cs%wave_structure_CSp, tot_en_mode(:,:,fr,m), full_halos=.true.)
      ! Pick out near-bottom and max horizontal baroclinic velocity values at each point
      do j=jsd,jed ; do i=isd,ied
        id_g = i + g%idg_offset ; jd_g = j + g%jdg_offset ! for debugging
        nzm = cs%wave_structure_CSp%num_intfaces(i,j)
        ub(i,j,fr,m) = cs%wave_structure_CSp%Uavg_profile(i,j,nzm)
        umax(i,j,fr,m) = maxval(cs%wave_structure_CSp%Uavg_profile(i,j,1:nzm))
        !! for debugging print profile, etc. Delete later
        !if(id_g .eq. 260 .and. &
        !   jd_g .eq. 50 .and. &
        !   tot_En_mode(i,j,1,1)>500.0) then
        !  print *, 'Profiles for mode ',m,' and frequency ',fr
        !  print *, 'id_g=', id_g, 'jd_g=', jd_g
        !  print *, 'c',m,'=',   cn(i,j,m)
        !  print *, 'nzm=',  nzm
        !  print *, 'z=',    CS%wave_structure_CSp%z_depths(i,j,1:nzm)
        !  print *, 'N2=',   CS%wave_structure_CSp%N2(i,j,1:nzm)
        !  print *, 'Ub=',   Ub(i,j,fr,m)
        !  print *, 'Umax=', Umax(i,j,fr,m)
        !  print *, 'Upro=', CS%wave_structure_CSp%Uavg_profile(i,j,1:nzm)
        !  print *, 'Wpro=', CS%wave_structure_CSp%W_profile(i,j,1:nzm)
        !  print *, 'En',m,'=',   tot_En_mode(i,j,fr,m)
        !  if (m==3) stop
        !endif ! for debug - delete later
      enddo ; enddo ! i-loop, j-loop
    enddo ; enddo ! fr-loop, m-loop
  endif ! apply_wave or _Froude_drag (Ub or Umax needed)
  ! Finally, apply loss
  if (cs%apply_wave_drag) then
    ! Calculate loss rate and apply loss over the time step
    call itidal_lowmode_loss(g, cs, nb, ub, cs%En, cs%TKE_itidal_loss_fixed, &
                             cs%TKE_itidal_loss, dt, full_halos=.false.)
  endif
  ! Check for En<0 - for debugging, delete later
  do m=1,cs%NMode ; do fr=1,cs%Nfreq ; do a=1,cs%nAngle
    do j=js,je ; do i=is,ie
      if (cs%En(i,j,a,fr,m)<0.0) then
        id_g = i + g%idg_offset ; jd_g = j + g%jdg_offset ! for debugging
        print *, 'After wave drag loss: En<0.0 at ig=', id_g, ', jg=', jd_g
        !stop
      endif
    enddo ; enddo
  enddo ; enddo ; enddo

  ! Extract the energy for mixing due to wave breaking-----------------------------
  if (cs%apply_Froude_drag) then
    ! Pick out maximum baroclinic velocity values; calculate Fr=max(u)/cg
    do m=1,cs%NMode ; do fr=1,cs%Nfreq
      freq2 = cs%frequency(fr)**2
      do j=jsd,jed ; do i=isd,ied
        id_g = i + g%idg_offset ; jd_g = j + g%jdg_offset ! for debugging
        ! Calculate horizontal phase velocity magnitudes
        f2 = 0.25*(g%CoriolisBu(i,j)**2 + g%CoriolisBu(i-1,j)**2 + &
                 g%CoriolisBu(i,j-1)**2 + g%CoriolisBu(i-1,j-1)**2 )
        kmag2 = (freq2 - f2) / (cn(i,j,m)**2 + cn_subro**2)
        c_phase = 0.0
        if (kmag2 > 0.0) then
          c_phase = sqrt(freq2/kmag2)
          nzm = cs%wave_structure_CSp%num_intfaces(i,j)
          fr2_max = (umax(i,j,fr,m)/c_phase)**2
          ! Dissipate energy if Fr>1; done here with an arbitrary time scale
          if (fr2_max > 1.0) then
            en_initial = sum(cs%En(i,j,:,fr,m)) ! for debugging
            ! Calculate effective decay rate (s-1) if breaking occurs over a time step
            loss_rate = (1/fr2_max - 1.0)/dt
            do a=1,cs%nAngle
              ! Determine effective dissipation rate (Wm-2)
              cs%TKE_Froude_loss(i,j,a,fr,m) = cs%En(i,j,a,fr,m) * abs(loss_rate)
              ! Update energy
              en_new = cs%En(i,j,a,fr,m)/fr2_max ! for debugging
              en_check = cs%En(i,j,a,fr,m) - cs%TKE_Froude_loss(i,j,a,fr,m)*dt ! for debugging
              ! Re-scale (reduce) energy due to breaking
              cs%En(i,j,a,fr,m) = cs%En(i,j,a,fr,m)/fr2_max
              ! Check (for debugging only)
             if (abs(en_new - en_check) > 1e-10) then
               call mom_error(warning, "MOM_internal_tides: something's wrong with Fr-breaking.")
               print *, "En_new=", en_new
               print *, "En_check=", en_check
              endif
            enddo
            ! Check (for debugging)
            delta_e_check = en_initial - sum(cs%En(i,j,:,fr,m))
            tke_froude_loss_check = abs(delta_e_check)/dt
            tke_froude_loss_tot = sum(cs%TKE_Froude_loss(i,j,:,fr,m))
            if (abs(tke_froude_loss_check - tke_froude_loss_tot) > 1e-10) then
              call mom_error(warning, "MOM_internal_tides: something's wrong with Fr energy update.")
              print *, "TKE_Froude_loss_check=", tke_froude_loss_check
              print *, "TKE_Froude_loss_tot=", tke_froude_loss_tot
            endif
          endif ! Fr2>1
        endif ! Kmag2>0
        cs%cp(i,j,fr,m) = c_phase
      enddo ; enddo
    enddo ; enddo
  endif
  ! Check for En<0 - for debugging, delete later
  do m=1,cs%NMode ; do fr=1,cs%Nfreq ; do a=1,cs%nAngle
    do j=js,je ; do i=is,ie
      if (cs%En(i,j,a,fr,m)<0.0) then
        id_g = i + g%idg_offset ; jd_g = j + g%jdg_offset
        print *, 'After Froude loss: En<0.0 at ig=', id_g, ', jg=', jd_g
        !stop
      endif
    enddo ; enddo
  enddo ; enddo ; enddo

  ! Check for energy conservation on computational domain.*************************
  do m=1,cs%NMode ; do fr=1,cs%Nfreq
    !print *, 'sum_En: mode(',m,'), freq(',fr,'):'
    call sum_en(g,cs,cs%En(:,:,:,fr,m),'prop_int_tide')
  enddo ; enddo

  ! Output diagnostics.************************************************************
  if (query_averaging_enabled(cs%diag)) then
    ! Output two-dimensional diagnostistics
    if (cs%id_tot_En > 0)     call post_data(cs%id_tot_En, tot_en, cs%diag)
    if (cs%id_itide_drag > 0) call post_data(cs%id_itide_drag, drag_scale, cs%diag)
    if (cs%id_TKE_itidal_input > 0) call post_data(cs%id_TKE_itidal_input, &
                                                   tke_itidal_input, cs%diag)

    ! Output 2-D energy density (summed over angles) for each freq and mode
    do m=1,cs%NMode ; do fr=1,cs%Nfreq ; if (cs%id_En_mode(fr,m) > 0) then
      tot_en(:,:) = 0.0
      do a=1,cs%nAngle ; do j=js,je ; do i=is,ie
        tot_en(i,j) = tot_en(i,j) + cs%En(i,j,a,fr,m)
      enddo ; enddo ; enddo
      call post_data(cs%id_En_mode(fr,m), tot_en, cs%diag)
    endif ; enddo ; enddo

    ! Output 3-D (i,j,a) energy density for each freq and mode
    do m=1,cs%NMode ; do fr=1,cs%Nfreq ; if (cs%id_En_ang_mode(fr,m) > 0) then
      call post_data(cs%id_En_ang_mode(fr,m), cs%En(:,:,:,fr,m) , cs%diag)
    endif ; enddo ; enddo

    ! Output 2-D energy loss (summed over angles, freq, modes)
    tot_leak_loss(:,:)   = 0.0
    tot_quad_loss(:,:)   = 0.0
    tot_itidal_loss(:,:) = 0.0
    tot_froude_loss(:,:) = 0.0
    tot_allprocesses_loss(:,:) = 0.0
    do m=1,cs%NMode ; do fr=1,cs%Nfreq ; do a=1,cs%nAngle ; do j=js,je ; do i=is,ie
      tot_leak_loss(i,j)   = tot_leak_loss(i,j)   + cs%TKE_leak_loss(i,j,a,fr,m)
      tot_quad_loss(i,j)   = tot_quad_loss(i,j)   + cs%TKE_quad_loss(i,j,a,fr,m)
      tot_itidal_loss(i,j) = tot_itidal_loss(i,j) + cs%TKE_itidal_loss(i,j,a,fr,m)
      tot_froude_loss(i,j) = tot_froude_loss(i,j) + cs%TKE_Froude_loss(i,j,a,fr,m)
    enddo ; enddo ; enddo ; enddo ; enddo
    do j=js,je ; do i=is,ie
      tot_allprocesses_loss(i,j) = tot_leak_loss(i,j) + tot_quad_loss(i,j) + &
                          tot_itidal_loss(i,j) + tot_froude_loss(i,j)
    enddo ; enddo
    cs%tot_leak_loss         = tot_leak_loss
    cs%tot_quad_loss         = tot_quad_loss
    cs%tot_itidal_loss       = tot_itidal_loss
    cs%tot_Froude_loss       = tot_froude_loss
    cs%tot_allprocesses_loss = tot_allprocesses_loss
    if (cs%id_tot_leak_loss > 0) then
      call post_data(cs%id_tot_leak_loss, tot_leak_loss, cs%diag)
    endif
    if (cs%id_tot_quad_loss > 0) then
      call post_data(cs%id_tot_quad_loss, tot_quad_loss, cs%diag)
    endif
    if (cs%id_tot_itidal_loss > 0) then
      call post_data(cs%id_tot_itidal_loss, tot_itidal_loss, cs%diag)
    endif
    if (cs%id_tot_Froude_loss > 0) then
      call post_data(cs%id_tot_Froude_loss, tot_froude_loss, cs%diag)
    endif
    if (cs%id_tot_allprocesses_loss > 0) then
      call post_data(cs%id_tot_allprocesses_loss, tot_allprocesses_loss, cs%diag)
    endif

    ! Output 2-D energy loss (summed over angles) for each freq and mode
    do m=1,cs%NMode ; do fr=1,cs%Nfreq
    if (cs%id_itidal_loss_mode(fr,m) > 0 .or. cs%id_allprocesses_loss_mode(fr,m) > 0) then
      itidal_loss_mode(:,:)       = 0.0 ! wave-drag processes (could do others as well)
      allprocesses_loss_mode(:,:) = 0.0 ! all processes summed together
      do a=1,cs%nAngle ; do j=js,je ; do i=is,ie
        itidal_loss_mode(i,j)       = itidal_loss_mode(i,j) + cs%TKE_itidal_loss(i,j,a,fr,m)
        allprocesses_loss_mode(i,j) = allprocesses_loss_mode(i,j) + &
                                 cs%TKE_leak_loss(i,j,a,fr,m) + cs%TKE_quad_loss(i,j,a,fr,m) + &
                                 cs%TKE_itidal_loss(i,j,a,fr,m) + cs%TKE_Froude_loss(i,j,a,fr,m)
      enddo ; enddo ; enddo
      call post_data(cs%id_itidal_loss_mode(fr,m), itidal_loss_mode, cs%diag)
      call post_data(cs%id_allprocesses_loss_mode(fr,m), allprocesses_loss_mode, cs%diag)
    endif ; enddo ; enddo

    ! Output 3-D (i,j,a) energy loss for each freq and mode
    do m=1,cs%NMode ; do fr=1,cs%Nfreq ; if (cs%id_itidal_loss_ang_mode(fr,m) > 0) then
      call post_data(cs%id_itidal_loss_ang_mode(fr,m), cs%TKE_itidal_loss(:,:,:,fr,m) , cs%diag)
    endif ; enddo ; enddo

    ! Output 2-D period-averaged horizontal near-bottom mode velocity for each freq and mode
    do m=1,cs%NMode ; do fr=1,cs%Nfreq ; if (cs%id_Ub_mode(fr,m) > 0) then
      call post_data(cs%id_Ub_mode(fr,m), ub(:,:,fr,m), cs%diag)
    endif ; enddo ; enddo

    ! Output 2-D horizontal phase velocity for each freq and mode
    do m=1,cs%NMode ; do fr=1,cs%Nfreq ; if (cs%id_cp_mode(fr,m) > 0) then
      call post_data(cs%id_cp_mode(fr,m), cs%cp(:,:,fr,m), cs%diag)
    endif ; enddo ; enddo

  endif

end subroutine propagate_int_tide

!> This subroutine checks for energy conservation on computational domain
subroutine sum_en(G, CS, En, label)
  type(ocean_grid_type),  intent(in)    :: G    !< The ocean's grid structure.
  type(int_tide_cs), pointer            :: CS
  real, dimension(G%isd:G%ied,G%jsd:G%jed,CS%NAngle), intent(in) :: En
  character(len=*), intent(in) :: label

  ! This subroutine checks for energy conservation on computational domain
  integer :: m,fr,a
  real :: En_sum, tmpForSumming, En_sum_diff, En_sum_pdiff
  integer :: seconds
  real :: Isecs_per_day = 1.0 / 86400.0
  real :: days

  call get_time(cs%Time, seconds)
  days = real(seconds) * Isecs_per_day

  en_sum = 0.0;
  tmpforsumming = 0.0
  do a=1,cs%nAngle
    tmpforsumming = global_area_mean(en(:,:,a),g)*g%areaT_global
    en_sum = en_sum + tmpforsumming
  enddo
  en_sum_diff = en_sum - cs%En_sum
  if (cs%En_sum .ne. 0.0) then
    en_sum_pdiff= (en_sum_diff/cs%En_sum)*100.0
  else
    en_sum_pdiff= 0.0;
  endif
  cs%En_sum = en_sum
  !! Print to screen
  !if (is_root_pe()) then
  !  print *, label,':','days =', days
  !  print *, 'En_sum=', En_sum
  !  print *, 'En_sum_diff=',    En_sum_diff
  !  print *, 'Percent change=', En_sum_pdiff, '%'
  !  !if (abs(En_sum_pdiff) > 1.0) then ; stop ; endif
  !endif

end subroutine sum_en

!> This subroutine calculates the energy lost from the propagating internal tide due to
!! scattering over small-scale roughness along the lines of Jayne & St. Laurent (2001).
subroutine itidal_lowmode_loss(G, CS, Nb, Ub, En, TKE_loss_fixed, TKE_loss, dt, full_halos)
  type(ocean_grid_type),     intent(in)    :: G  !< The ocean's grid structure.
  type(int_tide_cs),         pointer       :: CS
  real, dimension(G%isd:G%ied,G%jsd:G%jed), &
                             intent(in)    :: Nb !< Near-bottom stratification, in s-1.
  real, dimension(G%isd:G%ied,G%jsd:G%jed,CS%nFreq,CS%nMode), &
                             intent(inout) :: Ub !< Rms (over one period) near-bottom horizontal
                                                 !! mode velocity , in m s-1.
  real, dimension(G%isd:G%ied,G%jsd:G%jed), &
                             intent(in) :: TKE_loss_fixed !< Fixed part of energy loss,
                                                 !! in kg m-2 (rho*kappa*h^2).
  real, dimension(G%isd:G%ied,G%jsd:G%jed,CS%NAngle,CS%nFreq,CS%nMode), &
                             intent(inout) :: En !< Energy density of the internal waves, in J m-2.
  real, dimension(G%isd:G%ied,G%jsd:G%jed,CS%NAngle,CS%nFreq,CS%nMode), &
                             intent(out)   :: TKE_loss    !< Energy loss rate, in W m-2
                                                 !! (q*rho*kappa*h^2*N*U^2).
  real,                      intent(in)    :: dt !< Time increment, in s.
  logical,optional,          intent(in)    :: full_halos  !< If true, do the calculation over the
                                                 !! entirecomputational domain.

  ! This subroutine calculates the energy lost from the propagating internal tide due to
  ! scattering over small-scale roughness along the lines of Jayne & St. Laurent (2001).
  !
  ! Arguments:
  !  (in)      Nb - near-bottom stratification, in s-1.
  !  (in)      Ub - rms (over one period) near-bottom horizontal mode velocity , in m s-1.
  !  (inout)   En - energy density of the internal waves, in J m-2.
  !  (in)      TKE_loss_fixed - fixed part of energy loss, in kg m-2 (rho*kappa*h^2)
  !  (out)     TKE_loss - energy loss rate, in W m-2 (q*rho*kappa*h^2*N*U^2)
  !  (in)      dt - time increment, in s
  !  (in,opt)  full_halos - If true, do the calculation over the entire
  !                         computational domain.

  integer :: j,i,m,fr,a, is, ie, js, je
  real    :: En_tot          ! energy for a given mode, frequency, and point summed over angles
  real    :: TKE_loss_tot    ! dissipation for a given mode, frequency, and point summed over angles
  real    :: TKE_sum_check   ! temporary for check summing
  real    :: frac_per_sector ! fraction of energy in each wedge
  real    :: q_itides        ! fraction of energy actually lost to mixing (remainder, 1-q, is
                             ! assumed to stay in propagating mode for now - BDM)
  real    :: loss_rate       ! approximate loss rate for implicit calc, s-1
  real, parameter :: En_negl = 1e-30 ! negilibly small number to prevent division by zero

  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec

  q_itides = cs%q_itides

  if (present(full_halos)) then ; if (full_halos) then
    is = g%isd ; ie = g%ied ; js = g%jsd ; je = g%jed
  endif ; endif

  do j=js,je ; do i=is,ie ; do m=1,cs%nMode ; do fr=1,cs%nFreq

    ! Sum energy across angles
    en_tot = 0.0
    do a=1,cs%nAngle
      en_tot = en_tot + en(i,j,a,fr,m)
    enddo

    ! Calculate TKE loss rate; units of [W m-2] here.
    tke_loss_tot = q_itides * tke_loss_fixed(i,j) * nb(i,j) * ub(i,j,fr,m)**2

    ! Update energy remaining (this is a pseudo implicit calc)
    ! (E(t+1)-E(t))/dt = -TKE_loss(E(t+1)/E(t)), which goes to zero as E(t+1) goes to zero
    if (en_tot > 0.0) then
      do a=1,cs%nAngle
        frac_per_sector = en(i,j,a,fr,m)/en_tot
        tke_loss(i,j,a,fr,m) = frac_per_sector*tke_loss_tot           ! Wm-2
        loss_rate = tke_loss(i,j,a,fr,m) / (en(i,j,a,fr,m) + en_negl) ! s-1
        en(i,j,a,fr,m) = en(i,j,a,fr,m) / (1.0 + dt*loss_rate)
      enddo
    else
      ! no loss if no energy
      tke_loss(i,j,:,fr,m) = 0.0
    endif

    ! Update energy remaining (this is the old explicit calc)
    !if (En_tot > 0.0) then
    !  do a=1,CS%nAngle
    !    frac_per_sector = En(i,j,a,fr,m)/En_tot
    !    TKE_loss(i,j,a,fr,m) = frac_per_sector*TKE_loss_tot
    !    if(TKE_loss(i,j,a,fr,m)*dt <= En(i,j,a,fr,m))then
    !      En(i,j,a,fr,m) = En(i,j,a,fr,m) - TKE_loss(i,j,a,fr,m)*dt
    !    else
    !      call MOM_error(WARNING, "itidal_lowmode_loss: energy loss greater than avalable, "// &
    !                        " setting En to zero.")
    !      En(i,j,a,fr,m) = 0.0
    !    endif
    !  enddo
    !else
    !  ! no loss if no energy
    !  TKE_loss(i,j,:,fr,m) = 0.0
    !endif

  enddo ; enddo ; enddo ; enddo

end subroutine itidal_lowmode_loss

!> This subroutine extracts the energy lost from the propagating internal which has
!> been summed across all angles, frequencies, and modes for a given mechanism and location.
!> It can be called from another module to get values from this module's (private) CS.
subroutine get_lowmode_loss(i,j,G,CS,mechanism,TKE_loss_sum)
  integer,                intent(in)  :: i,j
  type(ocean_grid_type),  intent(in)  :: G            !< The ocean's grid structure
  type(int_tide_cs),      pointer     :: CS
  character(len=*),       intent(in)  :: mechanism
  real,                   intent(out) :: TKE_loss_sum !< Total energy loss rate due to specified
                                                      !! mechanism, in W m-2.

  ! This subroutine extracts the energy lost from the propagating internal which has
  ! been summed across all angles, frequencies, and modes for a given mechanism and location.
  ! It can be called from another module to get values from this module's (private) CS.
  !
  ! Arguments:
  !  (out)      TKE_loss_sum - total energy loss rate due to specified mechanism, in W m-2.

  if(mechanism == 'LeakDrag') tke_loss_sum = cs%tot_leak_loss(i,j)   ! not used for mixing yet
  if(mechanism == 'QuadDrag') tke_loss_sum = cs%tot_quad_loss(i,j)   ! not used for mixing yet
  if(mechanism == 'WaveDrag') tke_loss_sum = cs%tot_itidal_loss(i,j) ! currently used for mixing
  if(mechanism == 'Froude')   tke_loss_sum = cs%tot_Froude_loss(i,j) ! not used for mixing yet

end subroutine get_lowmode_loss

!> This subroutine does refraction on the internal waves at a single frequency.
subroutine refract(En, cn, freq, dt, G, NAngle, use_PPMang)
  type(ocean_grid_type),  intent(in)    :: G    !< The ocean's grid structure.
  integer,                intent(in)    :: NAngle
  real, dimension(G%isd:G%ied,G%jsd:G%jed,NAngle), &
                          intent(inout) :: En   !< The internal gravity wave energy density as a
                                                !! function of space and angular resolution,
                                                !! in J m-2 radian-1.
  real, dimension(G%isd:G%ied,G%jsd:G%jed),        &
                          intent(in)    :: cn   !< Baroclinic mode speed, in m s-1.
  real,                   intent(in)    :: freq !< Wave frequency, in s-1.
  real,                   intent(in)    :: dt   !< Time step, in s.
  logical,                intent(in)    :: use_PPMang !< If true, use PPM for advection rather
                                                !! than upwind.
  !  This subroutine does refraction on the internal waves at a single frequency.

  ! Arguments:
  ! (inout) En - the internal gravity wave energy density as a function of space
  !              and angular resolution, in J m-2 radian-1.
  ! (in)    cn - baroclinic mode speed, in m s-1
  ! (in)    freq - wave frequency, in s-1
  ! (in)    dt - time step, in s
  ! (in)    use_PPMang - if true, use PPM for advection rather than upwind

  integer, parameter :: stencil = 2
  real, dimension(SZI_(G),1-stencil:NAngle+stencil) :: &
    En2d
  real, dimension(1-stencil:NAngle+stencil) :: &
    cos_angle, sin_angle
  real, dimension(SZI_(G)) :: &
    Dk_Dt_Kmag, Dl_Dt_Kmag
  real, dimension(SZI_(G),0:nAngle) :: &
    Flux_E
  real, dimension(SZI_(G),SZJ_(G),1-stencil:NAngle+stencil) :: &
    CFL_ang
  real :: f2              ! The squared Coriolis parameter, in s-2.
  real :: favg            ! The average Coriolis parameter at a point, in s-1.
  real :: df2_dy, df2_dx  ! The x- and y- gradients of the squared Coriolis parameter, in s-2 m-1.
  real :: df_dy, df_dx    ! The x- and y- gradients of the Coriolis parameter, in s-1 m-1.
  real :: dlnCn_dx        ! The x-gradient of the wave speed divided by itself in m-1.
  real :: dlnCn_dy        ! The y-gradient of the wave speed divided by itself in m-1.
  real :: Angle_size, dt_Angle_size, angle
  real :: Ifreq, Kmag2, I_Kmag
  real, parameter :: cn_subRO = 1e-100
  integer :: is, ie, js, je, asd, aed, na
  integer :: i, j, a

  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; na = size(en,3)
  asd = 1-stencil ; aed = nangle+stencil

  ifreq = 1.0 / freq

  angle_size = (8.0*atan(1.0)) / (real(nangle))
  dt_angle_size = dt / angle_size

  do a=asd,aed
    angle = (real(A) - 0.5) * angle_size
    cos_angle(a) = cos(angle) ; sin_angle(a) = sin(angle)
  enddo

  !### There should also be refraction due to cn.grad(grid_orientation).
  cfl_ang(:,:,:) = 0.0;
  do j=js,je
  ! Copy En into angle space with halos.
    do a=1,na ; do i=is,ie
      en2d(i,a) = en(i,j,a)
    enddo ; enddo
    do a=asd,0 ; do i=is,ie
      en2d(i,a) = en2d(i,a+nangle)
      en2d(i,nangle+stencil+a) = en2d(i,stencil+a)
    enddo ; enddo

  ! Do the refraction.
    do i=is,ie
      f2 = 0.25*((g%CoriolisBu(i,j)**2 + g%CoriolisBu(i-1,j-1)**2) + &
                 (g%CoriolisBu(i,j-1)**2 + g%CoriolisBu(i-1,j)**2))
      favg = 0.25*((g%CoriolisBu(i,j) + g%CoriolisBu(i-1,j-1)) + &
                 (g%CoriolisBu(i,j-1) + g%CoriolisBu(i-1,j)))
      df2_dx = 0.5*((g%CoriolisBu(i,j)**2 + g%CoriolisBu(i,j-1)**2) - &
                    (g%CoriolisBu(i-1,j)**2 + g%CoriolisBu(i-1,j-1)**2)) * &
               g%IdxT(i,j)
      df_dx = 0.5*((g%CoriolisBu(i,j) + g%CoriolisBu(i,j-1)) - &
                    (g%CoriolisBu(i-1,j) + g%CoriolisBu(i-1,j-1))) * &
               g%IdxT(i,j)
      dlncn_dx = 0.5*( g%IdxCu(i,j) * (cn(i+1,j) - cn(i,j)) / &
                       (0.5*(cn(i+1,j) + cn(i,j)) + cn_subro) + &
                       g%IdxCu(i-1,j) * (cn(i,j) - cn(i-1,j)) / &
                       (0.5*(cn(i,j) + cn(i-1,j)) + cn_subro) )
      df2_dy = 0.5*((g%CoriolisBu(i,j)**2 + g%CoriolisBu(i-1,j)**2) - &
                    (g%CoriolisBu(i,j-1)**2 + g%CoriolisBu(i-1,j-1)**2)) * &
               g%IdyT(i,j)
      df_dy = 0.5*((g%CoriolisBu(i,j) + g%CoriolisBu(i-1,j)) - &
                    (g%CoriolisBu(i,j-1) + g%CoriolisBu(i-1,j-1))) * &
               g%IdyT(i,j)
      dlncn_dy = 0.5*( g%IdyCv(i,j) * (cn(i,j+1) - cn(i,j)) / &
                       (0.5*(cn(i,j+1) + cn(i,j)) + cn_subro) + &
                       g%IdyCv(i,j-1) * (cn(i,j) - cn(i,j-1)) / &
                       (0.5*(cn(i,j) + cn(i,j-1)) + cn_subro) )
      kmag2 = (freq**2 - f2) / (cn(i,j)**2 + cn_subro**2)
      if (kmag2 > 0.0) then
        i_kmag = 1.0 / sqrt(kmag2)
        dk_dt_kmag(i) = -ifreq * (favg*df_dx + (freq**2 - f2) * dlncn_dx) * i_kmag
        dl_dt_kmag(i) = -ifreq * (favg*df_dy + (freq**2 - f2) * dlncn_dy) * i_kmag
      else
        dk_dt_kmag(i) = 0.0
        dl_dt_kmag(i) = 0.0
      endif
    enddo

    ! Determine the energy fluxes in angular orientation space.
    do a=asd,aed ; do i=is,ie
      cfl_ang(i,j,a) = (cos_angle(a) * dl_dt_kmag(i) - sin_angle(a) * dk_dt_kmag(i)) * &
                dt_angle_size
      if (abs(cfl_ang(i,j,a)) > 1.0) then
        call mom_error(warning, "refract: CFL exceeds 1.", .true.)
        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
      endif
    enddo; enddo

    ! Advect in angular space
    if(.not.use_ppmang) then
      ! Use simple upwind
      do  a=0,na ; do i=is,ie
        if (cfl_ang(i,j,a) > 0.0) then
          flux_e(i,a) = cfl_ang(i,j,a) * en2d(i,a)
        else
          flux_e(i,a) = cfl_ang(i,j,a) * en2d(i,a+1)
        endif
      enddo; enddo
    else
      ! Use PPM
      do i=is,ie
        call ppm_angular_advect(en2d(i,:),cfl_ang(i,j,:),flux_e(i,:),nangle,dt,stencil)
      enddo
    endif

  ! Update and copy back to En.
    do a=1,na ; do i=is,ie
      !if(En2d(i,a)+(Flux_E(i,A-1)-Flux_E(i,A)) < 0.0)then ! for debugging
      !  print *,"refract: OutFlux>Available" ; !stop
      !endif
      en(i,j,a) = en2d(i,a) + (flux_e(i,a-1) - flux_e(i,a))
    enddo ; enddo
  enddo ! j-loop
end subroutine refract

!> This subroutine calculates the 1-d flux for advection in angular space
!! using a monotonic piecewise parabolic scheme. Should be within i and j spatial
!! loops.
subroutine ppm_angular_advect(En2d, CFL_ang, Flux_En, NAngle, dt, halo_ang)
  integer,                     intent(in)    :: NAngle
  real,                        intent(in)    :: dt
  integer,                     intent(in)    :: halo_ang
  real, dimension(1-halo_ang:NAngle+halo_ang),   &
                               intent(in)    :: En2d
  real, dimension(1-halo_ang:NAngle+halo_ang),   &
                               intent(in)    :: CFL_ang
  real, dimension(0:NAngle),   intent(out)   :: Flux_En

  !   This subroutine calculates the 1-d flux for advection in angular space
  ! using a monotonic piecewise parabolic scheme. Should be within i and j spatial
  ! loops
  real :: flux
  real :: u_ang
  real :: Angle_size
  real :: I_Angle_size
  real :: I_dt
  integer :: a
  real :: aR, aL, dMx, dMn, Ep, Ec, Em, dA, mA, a6

  i_dt = 1 / dt
  angle_size = (8.0*atan(1.0)) / (real(nangle))
  i_angle_size = 1 / angle_size
  flux_en(:) = 0

  do a=0,nangle
    u_ang = cfl_ang(a)*angle_size*i_dt
    if (u_ang >= 0.0) then
      ! Implementation of PPM-H3
      ep = en2d(a+1)*i_angle_size !MEAN ANGULAR ENERGY DENSITY FOR WEDGE (Jm-2/rad)
      ec = en2d(a)  *i_angle_size !MEAN ANGULAR ENERGY DENSITY FOR WEDGE (Jm-2/rad)
      em = en2d(a-1)*i_angle_size !MEAN ANGULAR ENERGY DENSITY FOR WEDGE (Jm-2/rad)
      al = ( 5.*ec + ( 2.*em - ep ) )/6. ! H3 estimate
      al = max( min(ec,em), al) ; al = min( max(ec,em), al) ! Bound
      ar = ( 5.*ec + ( 2.*ep - em ) )/6. ! H3 estimate
      ar = max( min(ec,ep), ar) ; ar = min( max(ec,ep), ar) ! Bound
      da = ar - al ; ma = 0.5*( ar + al )
      if ((ep-ec)*(ec-em) <= 0.) then
        al = ec ; ar = ec ! PCM for local extremum
      elseif ( da*(ec-ma) > (da*da)/6. ) then
        al = 3.*ec - 2.*ar !?
      elseif ( da*(ec-ma) < - (da*da)/6. ) then
        ar = 3.*ec - 2.*al !?
      endif
      a6 = 6.*ec - 3. * (ar + al) ! Curvature
      ! CALCULATE FLUX RATE (Jm-2/s)
      flux = u_ang*( ar + 0.5 * cfl_ang(a) * ( ( al - ar ) + a6 * ( 1. - 2./3. * cfl_ang(a) ) ) )
      !flux = u_ang*( aR - 0.5 * CFL_ang(A) * ( ( aR - aL ) - a6 * ( 1. - 2./3. * CFL_ang(A) ) ) )
      ! CALCULATE AMOUNT FLUXED (Jm-2)
      flux_en(a) = dt * flux
      !Flux_En(A) = (dt * I_Angle_size) * flux
    else
      ! Implementation of PPM-H3
      ep = en2d(a+2)*i_angle_size !MEAN ANGULAR ENERGY DENSITY FOR WEDGE (Jm-2/rad)
      ec = en2d(a+1)*i_angle_size !MEAN ANGULAR ENERGY DENSITY FOR WEDGE (Jm-2/rad)
      em = en2d(a)  *i_angle_size !MEAN ANGULAR ENERGY DENSITY FOR WEDGE (Jm-2/rad)
      al = ( 5.*ec + ( 2.*em - ep ) )/6. ! H3 estimate
      al = max( min(ec,em), al) ; al = min( max(ec,em), al) ! Bound
      ar = ( 5.*ec + ( 2.*ep - em ) )/6. ! H3 estimate
      ar = max( min(ec,ep), ar) ; ar = min( max(ec,ep), ar) ! Bound
      da = ar - al ; ma = 0.5*( ar + al )
      if ((ep-ec)*(ec-em) <= 0.) then
        al = ec ; ar = ec ! PCM for local extremum
      elseif ( da*(ec-ma) > (da*da)/6. ) then
        al = 3.*ec - 2.*ar
      elseif ( da*(ec-ma) < - (da*da)/6. ) then
        ar = 3.*ec - 2.*al
      endif
      a6 = 6.*ec - 3. * (ar + al) ! Curvature
      ! CALCULATE FLUX RATE (Jm-2/s)
      flux = u_ang*( ar + 0.5 * cfl_ang(a) * ( ( al - ar ) + a6 * ( 1. - 2./3. * cfl_ang(a) ) ) )
      !flux = u_ang*( aL + 0.5 * CFL_ang(A) * ( ( aR - aL ) + a6 * ( 1. - 2./3. * CFL_ang(A) ) ) )
      ! CALCULATE AMOUNT FLUXED (Jm-2)
      flux_en(a) = dt * flux
      !Flux_En(A) = (dt * I_Angle_size) * flux
    endif
  enddo
end subroutine ppm_angular_advect

!> This subroutine does refraction on the internal waves at a single frequency.
subroutine propagate(En, cn, freq, dt, G, CS, NAngle)
  type(ocean_grid_type),  intent(inout) :: G    !< The ocean's grid structure.
  integer,                intent(in)    :: NAngle
  real, dimension(G%isd:G%ied,G%jsd:G%jed,NAngle), &
                          intent(inout) :: En   !< The internal gravity wave energy density as a
                                                !! function of space and angular resolution,
                                                !! in J m-2 radian-1.
  real, dimension(G%isd:G%ied,G%jsd:G%jed),        &
                          intent(in)    :: cn   !< Baroclinic mode speed, in m s-1.
  real,                   intent(in)    :: freq !< Wave frequency, in s-1.
  real,                   intent(in)    :: dt   !< Time step, in s.
  type(int_tide_cs),      pointer       :: CS
  !  This subroutine does refraction on the internal waves at a single frequency.

  ! Arguments:
  ! (inout) En - the internal gravity wave energy density as a function of space
  !              and angular resolution, in J m-2 radian-1.
  ! (in)    cn - baroclinic mode speed, in m s-1
  ! (in)    freq - wave frequency, in s-1
  ! (in)    dt - time step, in s

  real, dimension(G%IsdB:G%IedB,G%JsdB:G%JedB) :: &
    speed  ! The magnitude of the group velocity at the q points for corner adv, in m s-1.
  integer, parameter :: stencil = 2
  real, dimension(SZIB_(G),SZJ_(G)) :: &
    speed_x  ! The magnitude of the group velocity at the Cu points, in m s-1.
  real, dimension(SZI_(G),SZJB_(G)) :: &
    speed_y  ! The magnitude of the group velocity at the Cv points, in m s-1.
  real, dimension(0:NAngle) :: &
    cos_angle, sin_angle
  real, dimension(NAngle) :: &
    Cgx_av, Cgy_av, dCgx, dCgy
  real :: f2   ! The squared Coriolis parameter, in s-2.
  real :: Angle_size, I_Angle_size, angle
  real :: Ifreq, freq2
  real, parameter :: cn_subRO = 1e-100
  type(loop_bounds_type) :: LB
  integer :: is, ie, js, je, asd, aed, na
  integer :: ish, ieh, jsh, jeh
  integer :: i, j, a

  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; na = size(en,3)
  asd = 1-stencil ; aed = nangle+stencil

  ifreq = 1.0 / freq
  freq2 = freq**2

  ! Define loop bounds: Need extensions on j-loop so propagate_y
  ! (done after propagate_x) will have updated values in the halo
  ! for correct PPM reconstruction. Use if no teleporting and
  ! no pass_var between propagate_x and propagate_y.
  !jsh = js-3 ; jeh = je+3 ; ish = is ; ieh = ie

  ! Define loop bounds: Need 1-pt extensions on loops because
  ! teleporting eats up a halo point. Use if teleporting.
  ! Also requires pass_var before propagate_y.
  jsh = js-1 ; jeh = je+1 ; ish = is-1 ; ieh = ie+1

  angle_size = (8.0*atan(1.0)) / real(nangle)
  i_angle_size = 1.0 / angle_size

  if (cs%corner_adv) then
    ! IMPLEMENT CORNER ADVECTION IN HORIZONTAL--------------------
    ! FIND AVERAGE GROUP VELOCITY (SPEED) AT CELL CORNERS;
    ! NOTE: THIS HAS NOT BE ADAPTED FOR REFLECTION YET (BDM)!!
    ! Fix indexing here later
    speed(:,:) = 0;
    do j=jsh-1,jeh ; do i=ish-1,ieh
      f2 = g%CoriolisBu(i,j)**2
      speed(i,j) = 0.25*(cn(i,j) + cn(i+1,j) + cn(i+1,j+1) + cn(i,j+1)) * &
                     sqrt(max(freq2 - f2, 0.0)) * ifreq
    enddo ; enddo
    do a=1,na
      ! Apply the propagation WITH CORNER ADVECTION/FINITE VOLUME APPROACH.
      lb%jsh = js ; lb%jeh = je ; lb%ish = is ; lb%ieh = ie
      call propagate_corner_spread(en(:,:,a), a, nangle, speed, dt, g, cs, lb)
    end do ! a-loop
  else
    ! IMPLEMENT PPM ADVECTION IN HORIZONTAL-----------------------
    ! These could be in the control structure, as they do not vary.
    do a=0,na
      ! These are the angles at the cell edges...
      angle = (real(A) - 0.5) * angle_size
      cos_angle(a) = cos(angle) ; sin_angle(a) = sin(angle)
    enddo

    do a=1,na
      cgx_av(a) = (sin_angle(a) - sin_angle(a-1)) * i_angle_size
      cgy_av(a) = -(cos_angle(a) - cos_angle(a-1)) * i_angle_size
      dcgx(a) = sqrt(0.5 + 0.5*(sin_angle(a)*cos_angle(a) - &
                                sin_angle(a-1)*cos_angle(a-1)) * i_angle_size - &
                     cgx_av(a)**2)
      dcgy(a) = sqrt(0.5 - 0.5*(sin_angle(a)*cos_angle(a) - &
                                sin_angle(a-1)*cos_angle(a-1)) * i_angle_size - &
                     cgy_av(a)**2)
    enddo

    do j=jsh,jeh ; do i=ish-1,ieh
      f2 = 0.5*(g%CoriolisBu(i,j)**2 + g%CoriolisBu(i,j-1)**2)
      speed_x(i,j) = 0.5*(cn(i,j) + cn(i+1,j)) * g%mask2dCu(i,j) * &
                     sqrt(max(freq2 - f2, 0.0)) * ifreq
    enddo ; enddo
    do j=jsh-1,jeh ; do i=ish,ieh
      f2 = 0.5*(g%CoriolisBu(i,j)**2 + g%CoriolisBu(i-1,j)**2)
      speed_y(i,j) = 0.5*(cn(i,j) + cn(i,j+1)) * g%mask2dCv(i,j) * &
                     sqrt(max(freq2 - f2, 0.0)) * ifreq
    enddo ; enddo

    ! Apply propagation in x-direction (reflection included)
    lb%jsh = jsh ; lb%jeh = jeh ; lb%ish = ish ; lb%ieh = ieh
    call propagate_x(en(:,:,:), speed_x, cgx_av(:), dcgx(:), dt, g, cs%nAngle, cs, lb)

    ! Check for energy conservation on computational domain (for debugging)
    !call sum_En(G,CS,En(:,:,:),'post-propagate_x')

    ! Update halos
    call pass_var(en(:,:,:),g%domain)

    ! Apply propagation in y-direction (reflection included)
    ! LB%jsh = js ; LB%jeh = je ; LB%ish = is ; LB%ieh = ie ! Use if no teleport
    lb%jsh = jsh ; lb%jeh = jeh ; lb%ish = ish ; lb%ieh = ieh
    call propagate_y(en(:,:,:), speed_y, cgy_av(:), dcgy(:), dt, g, cs%nAngle, cs, lb)

    ! Check for energy conservation on computational domain (for debugging)
    !call sum_En(G,CS,En(:,:,:),'post-propagate_y')

  endif
end subroutine propagate

!> This subroutine does first-order corner advection. It was written with the hopes
!! of smoothing out the garden sprinkler effect, but is too numerically diffusive to
!! be of much use as of yet. It is not yet compatible with reflection schemes (BDM).
subroutine propagate_corner_spread(En, energized_wedge, NAngle, speed, dt, G, CS, LB)
  type(ocean_grid_type),  intent(in)    :: G     !< The ocean's grid structure.
  real, dimension(G%isd:G%ied,G%jsd:G%jed),   &
                          intent(inout) :: En    !< The energy density integrated over an angular
                                                 !! band, in W m-2, intent in/out.
  real, dimension(G%IsdB:G%IedB,G%Jsd:G%Jed), &
                          intent(in)    :: speed !< The magnitude of the group velocity at the cell
                                                 !! corner points, in m s-1.
  integer,                intent(in)    :: energized_wedge !< Index of current ray direction.
  integer,                intent(in)    :: NAngle
  real,                   intent(in)    :: dt    !< Time increment in s.
  type(int_tide_cs),      pointer       :: CS    !< The control structure returned by a previous
                                                 !! call to continuity_PPM_init.
  type(loop_bounds_type), intent(in)    :: LB    !< A structure with the active energy loop bounds.

  !   This subroutine does first-order corner advection. It was written with the hopes
  ! of smoothing out the garden sprinkler effect, but is too numerically diffusive to
  ! be of much use as of yet. It is not yet compatible with reflection schemes (BDM).

  ! Arguments: En - The energy density integrated over an angular band, in W m-2,
  !                 intent in/out.
  !  (in)      energized_wedge - index of current ray direction
  !  (in)      speed - The magnitude of the group velocity at the cell corner
  !                      points, in m s-1.
  !  (in)      dt - Time increment in s.
  !  (in)      G - The ocean's grid structure.
  !  (in)      CS - The control structure returned by a previous call to
  !                 continuity_PPM_init.
  !  (in)      LB - A structure with the active energy loop bounds.

  integer :: i, j, k, ish, ieh, jsh, jeh, m
  real :: TwoPi, Angle_size
  real :: energized_angle ! angle through center of current wedge
  real :: theta ! angle at edge of wedge
  real :: Nsubrays     ! number of sub-rays for averaging;
                       ! count includes the two rays that bound the current wedge,
                       ! i.e. those at -dtheta/2 and +dtheta/2 from energized angle
  real :: I_Nsubwedges ! inverse of number of sub-wedges
  real :: cos_thetaDT, sin_thetaDT ! cos(theta)*dt, sin(theta)*dt
  real :: xNE,xNW,xSW,xSE,yNE,yNW,ySW,ySE ! corner point coordinates of advected fluid parcel
  real :: CFL_xNE,CFL_xNW,CFL_xSW,CFL_xSE,CFL_yNE,CFL_yNW,CFL_ySW,CFL_ySE,CFL_max
  real :: xN,xS,xE,xW,yN,yS,yE,yW ! intersection point coordinates of parcel edges and grid
  real :: xCrn,yCrn ! grid point contained within advected fluid parcel
  real :: xg,yg ! grid point of interest
  real :: slopeN,slopeW,slopeS,slopeE, bN,bW,bS,bE ! parameters defining parcel sides
  real :: aNE,aN,aNW,aW,aSW,aS,aSE,aE,aC ! sub-areas of advected parcel
  real :: a_total ! total area of advected parcel
  real :: a1,a2,a3,a4 ! areas used in calculating polygon areas (sub-areas) of advected parcel
  real, dimension(G%IsdB:G%IedB,G%Jsd:G%Jed) :: x,y ! coordinates of cell corners
  real, dimension(G%IsdB:G%IedB,G%Jsd:G%Jed) :: Idx,Idy ! inverse of dx,dy at cell corners
  real, dimension(G%IsdB:G%IedB,G%Jsd:G%Jed) :: dx,dy ! dx,dy at cell corners
  real, dimension(2) :: E_new ! energy in cell after advection for subray; set size here to
                              ! define Nsubrays - this should be made an input option later!

  ish = lb%ish ; ieh = lb%ieh ; jsh = lb%jsh ; jeh = lb%jeh
  twopi = (8.0*atan(1.0))
  nsubrays = real(size(e_new))
  i_nsubwedges = 1./(nsubrays - 1)

  angle_size = twopi / real(nangle)
  energized_angle = angle_size * real(energized_wedge - 1) ! for a=1 aligned with x-axis
  !energized_angle = Angle_size * real(energized_wedge - 1) + 2.0*Angle_size !
  !energized_angle = Angle_size * real(energized_wedge - 1) + 0.5*Angle_size !
  x = g%geoLonBu
  y = g%geoLatBu
  idx = g%IdxBu; dx = g%dxBu
  idy = g%IdyBu; dy = g%dyBu

  do j=jsh,jeh; do i=ish,ieh
    do m=1,int(nsubrays)
      theta = energized_angle - 0.5*angle_size + real(m - 1)*Angle_size*I_Nsubwedges
      if (theta < 0.0) then
        theta = theta + twopi
      elseif (theta > twopi) then
        theta = theta - twopi
      endif
      cos_thetadt = cos(theta)*dt
      sin_thetadt = sin(theta)*dt

      ! corner point coordinates of advected fluid parcel ----------
      xg = x(i,j); yg = y(i,j)
      xne = xg - speed(i,j)*cos_thetadt
      yne = yg - speed(i,j)*sin_thetadt
      cfl_xne = (xg-xne)*idx(i,j)
      cfl_yne = (yg-yne)*idy(i,j)

      xg = x(i-1,j); yg = y(i-1,j)
      xnw = xg - speed(i-1,j)*cos_thetadt
      ynw = yg - speed(i-1,j)*sin_thetadt
      cfl_xnw = (xg-xnw)*idx(i-1,j)
      cfl_ynw = (yg-ynw)*idy(i-1,j)

      xg = x(i-1,j-1); yg = y(i-1,j-1)
      xsw = xg - speed(i-1,j-1)*cos_thetadt
      ysw = yg - speed(i-1,j-1)*sin_thetadt
      cfl_xsw = (xg-xsw)*idx(i-1,j-1)
      cfl_ysw = (yg-ysw)*idy(i-1,j-1)

      xg = x(i,j-1); yg = y(i,j-1)
      xse = xg - speed(i,j-1)*cos_thetadt
      yse = yg - speed(i,j-1)*sin_thetadt
      cfl_xse = (xg-xse)*idx(i,j-1)
      cfl_yse = (yg-yse)*idy(i,j-1)

      cfl_max = max(abs(cfl_xne),abs(cfl_xnw),abs(cfl_xsw), &
                    abs(cfl_xse),abs(cfl_yne),abs(cfl_ynw), &
                    abs(cfl_ysw),abs(cfl_yse))
      if (cfl_max > 1.0) then
        call mom_error(warning, "propagate_corner_spread: CFL exceeds 1.", .true.)
      endif

      ! intersection point coordinates of parcel edges and cell edges ---
      if (0.0 <= theta .and. theta < 0.25*twopi) then
          xn = x(i-1,j-1)
          yw = y(i-1,j-1)
      elseif (0.25*twopi <= theta .and. theta < 0.5*twopi) then
          xn = x(i,j-1)
          yw = y(i,j-1)
      elseif (0.5*twopi <= theta .and. theta < 0.75*twopi) then
          xn = x(i,j)
          yw = y(i,j)
      elseif (0.75*twopi <= theta .and. theta <= 1.00*twopi) then
          xn = x(i-1,j)
          yw = y(i-1,j)
      endif
      xs = xn
      ye = yw

      ! north intersection
      slopen = (yne - ynw)/(xne - xnw)
      bn = -slopen*xne + yne
      yn = slopen*xn + bn
      ! west intersection
      if (xnw == xsw) then
        xw = xnw
      else
        slopew = (ynw - ysw)/(xnw - xsw)
        bw = -slopew*xnw + ynw
        xw = (yw - bw)/slopew
      endif
      ! south intersection
      slopes = (ysw - yse)/(xsw - xse)
      bs = -slopes*xsw + ysw
      ys = slopes*xs + bs
      ! east intersection
      if (xne == xse) then
        xe = xne
      else
        slopee = (yse - yne)/(xse - xne)
        be = -slopee*xse + yse
        xe = (ye - be)/slopee
      endif

      ! areas --------------------------------------------
      ane = 0.0; an = 0.0; anw = 0.0; ! initialize areas
      aw = 0.0; asw = 0.0; as = 0.0; ! initialize areas
      ase = 0.0; ae = 0.0; ac = 0.0; ! initialize areas
      if (0.0 <= theta .and. theta < 0.25*twopi) then
          xcrn = x(i-1,j-1); ycrn = y(i-1,j-1);
          ! west area
          a1 = (yn - ycrn)*(0.5*(xn + xcrn))
          a2 = (ycrn - yw)*(0.5*(xcrn + xw))
          a3 = (yw - ynw)*(0.5*(xw + xnw))
          a4 = (ynw - yn)*(0.5*(xnw + xn))
          aw = a1 + a2 + a3 + a4
          ! southwest area
          a1 = (ycrn - ys)*(0.5*(xcrn + xs))
          a2 = (ys - ysw)*(0.5*(xs + xsw))
          a3 = (ysw - yw)*(0.5*(xsw + xw))
          a4 = (yw - ycrn)*(0.5*(xw + xcrn))
          asw = a1 + a2 + a3 + a4
          ! south area
          a1 = (ye - yse)*(0.5*(xe + xse))
          a2 = (yse - ys)*(0.5*(xse + xs))
          a3 = (ys - ycrn)*(0.5*(xs + xcrn))
          a4 = (ycrn - ye)*(0.5*(xcrn + xe))
          as = a1 + a2 + a3 + a4
          ! area within cell
          a1 = (yne - ye)*(0.5*(xne + xe))
          a2 = (ye - ycrn)*(0.5*(xe + xcrn))
          a3 = (ycrn - yn)*(0.5*(xcrn + xn))
          a4 = (yn - yne)*(0.5*(xn + xne))
          ac = a1 + a2 + a3 + a4
      elseif (0.25*twopi <= theta .and. theta < 0.5*twopi) then
          xcrn = x(i,j-1); ycrn = y(i,j-1);
          ! south area
          a1 = (ycrn - ys)*(0.5*(xcrn + xs))
          a2 = (ys - ysw)*(0.5*(xs + xsw))
          a3 = (ysw - yw)*(0.5*(xsw + xw))
          a4 = (yw - ycrn)*(0.5*(xw + xcrn))
          as = a1 + a2 + a3 + a4
          ! southeast area
          a1 = (ye - yse)*(0.5*(xe + xse))
          a2 = (yse - ys)*(0.5*(xse + xs))
          a3 = (ys - ycrn)*(0.5*(xs + xcrn))
          a4 = (ycrn - ye)*(0.5*(xcrn + xe))
          ase = a1 + a2 + a3 + a4
          ! east area
          a1 = (yne - ye)*(0.5*(xne + xe))
          a2 = (ye - ycrn)*(0.5*(xe + xcrn))
          a3 = (ycrn - yn)*(0.5*(xcrn + xn))
          a4 = (yn - yne)*(0.5*(xn + xne))
          ae = a1 + a2 + a3 + a4
          ! area within cell
          a1 = (yn - ycrn)*(0.5*(xn + xcrn))
          a2 = (ycrn - yw)*(0.5*(xcrn + xw))
          a3 = (yw - ynw)*(0.5*(xw + xnw))
          a4 = (ynw - yn)*(0.5*(xnw + xn))
          ac = a1 + a2 + a3 + a4
      elseif (0.5*twopi <= theta .and. theta < 0.75*twopi) then
          xcrn = x(i,j); ycrn = y(i,j);
          ! east area
          a1 = (ye - yse)*(0.5*(xe + xse))
          a2 = (yse - ys)*(0.5*(xse + xs))
          a3 = (ys - ycrn)*(0.5*(xs + xcrn))
          a4 = (ycrn - ye)*(0.5*(xcrn + xe))
          ae = a1 + a2 + a3 + a4
          ! northeast area
          a1 = (yne - ye)*(0.5*(xne + xe))
          a2 = (ye - ycrn)*(0.5*(xe + xcrn))
          a3 = (ycrn - yn)*(0.5*(xcrn + xn))
          a4 = (yn - yne)*(0.5*(xn + xne))
          ane = a1 + a2 + a3 + a4
          ! north area
          a1 = (yn - ycrn)*(0.5*(xn + xcrn))
          a2 = (ycrn - yw)*(0.5*(xcrn + xw))
          a3 = (yw - ynw)*(0.5*(xw + xnw))
          a4 = (ynw - yn)*(0.5*(xnw + xn))
          an = a1 + a2 + a3 + a4
          ! area within cell
          a1 = (ycrn - ys)*(0.5*(xcrn + xs))
          a2 = (ys - ysw)*(0.5*(xs + xsw))
          a3 = (ysw - yw)*(0.5*(xsw + xw))
          a4 = (yw - ycrn)*(0.5*(xw + xcrn))
          ac = a1 + a2 + a3 + a4
      elseif (0.75*twopi <= theta .and. theta <= 1.00*twopi) then
          xcrn = x(i-1,j); ycrn = y(i-1,j);
          ! north area
          a1 = (yne - ye)*(0.5*(xne + xe))
          a2 = (ye - ycrn)*(0.5*(xe + xcrn))
          a3 = (ycrn - yn)*(0.5*(xcrn + xn))
          a4 = (yn - yne)*(0.5*(xn + xne))
          an = a1 + a2 + a3 + a4
          ! northwest area
          a1 = (yn - ycrn)*(0.5*(xn + xcrn))
          a2 = (ycrn - yw)*(0.5*(xcrn + xw))
          a3 = (yw - ynw)*(0.5*(xw + xnw))
          a4 = (ynw - yn)*(0.5*(xnw + xn))
          anw = a1 + a2 + a3 + a4;
          ! west area
          a1 = (ycrn - ys)*(0.5*(xcrn + xs))
          a2 = (ys - ysw)*(0.5*(xs + xsw))
          a3 = (ysw - yw)*(0.5*(xsw + xw))
          a4 = (yw - ycrn)*(0.5*(xw + xcrn))
          aw = a1 + a2 + a3 + a4
          ! area within cell
          a1 = (ye - yse)*(0.5*(xe + xse))
          a2 = (yse - ys)*(0.5*(xse + xs))
          a3 = (ys - ycrn)*(0.5*(xs + xcrn))
          a4 = (ycrn - ye)*(0.5*(xcrn + xe))
          ac = a1 + a2 + a3 + a4
      endif

      ! energy weighting ----------------------------------------
      a_total = ane + an + anw + aw + asw + as + ase + ae + ac
      e_new(m) = (ane*en(i+1,j+1) + an*en(i,j+1) + anw*en(i-1,j+1) + &
                  aw*en(i-1,j) + asw*en(i-1,j-1) + as*en(i,j-1) + &
                  ase*en(i+1,j-1) + ae*en(i+1,j) + ac*en(i,j)) / (dx(i,j)*dy(i,j))
    enddo ! m-loop
    ! update energy in cell
    en(i,j) = sum(e_new)/nsubrays
  enddo; enddo
end subroutine propagate_corner_spread

! #@# This subroutine needs a doxygen description
subroutine propagate_x(En, speed_x, Cgx_av, dCgx, dt, G, Nangle, CS, LB)
  type(ocean_grid_type),   intent(in)    :: G  !< The ocean's grid structure.
  integer,                 intent(in)    :: NAngle
  real, dimension(G%isd:G%ied,G%jsd:G%jed,Nangle),   &
                           intent(inout) :: En !< The energy density integrated over an angular
                                               !! band, in J m-2, intent in/out.
  real, dimension(G%IsdB:G%IedB,G%jsd:G%jed),        &
                           intent(in)    :: speed_x !< The magnitude of the group velocity at the
                                               !! Cu points, in m s-1.
  real, dimension(Nangle), intent(in)    :: Cgx_av, dCgx
  real,                    intent(in)    :: dt !< Time increment in s.
  type(int_tide_cs),       pointer       :: CS !< The control structure returned by a previous call
                                               !! to continuity_PPM_init.
  type(loop_bounds_type),  intent(in)    :: LB !< A structure with the active energy loop bounds.

  ! Arguments: En - The energy density integrated over an angular band, in J m-2,
  !                 intent in/out.
  !  (in)      speed_x - The magnitude of the group velocity at the Cu
  !                      points, in m s-1.
  !  (in)      dt - Time increment in s.
  !  (in)      G - The ocean's grid structure.
  !  (in)      CS - The control structure returned by a previous call to
  !                 continuity_PPM_init.
  !  (in)      LB - A structure with the active energy loop bounds.
  real, dimension(SZI_(G),SZJ_(G)) :: &
    EnL, EnR    ! Left and right face energy densities, in J m-2.
  real, dimension(SZIB_(G),SZJ_(G)) :: &
    flux_x      ! The internal wave energy flux, in J s-1.
  real, dimension(SZIB_(G)) :: &
    cg_p, cg_m, flux1, flux2
  !real, dimension(SZI_(G),SZJB_(G),Nangle) :: En_m, En_p
  real, dimension(SZI_(G),SZJB_(G),Nangle) :: &
    Fdt_m, Fdt_p! Left and right energy fluxes, in J
  integer :: i, j, k, ish, ieh, jsh, jeh, a

  ish = lb%ish ; ieh = lb%ieh ; jsh = lb%jsh ; jeh = lb%jeh
  do a=1,nangle
  ! This sets EnL and EnR.
    if (cs%upwind_1st) then
      do j=jsh,jeh ; do i=ish-1,ieh+1
        enl(i,j) = en(i,j,a) ; enr(i,j) = en(i,j,a)
      enddo ; enddo
    else
      call ppm_reconstruction_x(en(:,:,a), enl, enr, g, lb, simple_2nd=cs%simple_2nd)
    endif

    do j=jsh,jeh
      ! This is done once with single speed (GARDEN SPRINKLER EFFECT POSSIBLE)
      do i=ish-1,ieh
        cg_p(i) = speed_x(i,j) * (cgx_av(a))
      enddo
      call zonal_flux_en(cg_p, en(:,j,a), enl(:,j), enr(:,j), flux1, &
                         dt, g, j, ish, ieh, cs%vol_CFL)
      do i=ish-1,ieh ; flux_x(i,j) = flux1(i); enddo
    enddo

    do j=jsh,jeh ; do i=ish,ieh
      fdt_m(i,j,a) = dt*flux_x(i-1,j) ! left face influx  (J)
      fdt_p(i,j,a) = -dt*flux_x(i,j)  ! right face influx (J)
    enddo ; enddo

    ! test with old (take out later)
    !do j=LB%jsh,LB%jeh ; do i=LB%ish,LB%ieh
    !  En(i,j,a) = En(i,j,a) - dt* G%IareaT(i,j) * (flux_x(I,j) - flux_x(I-1,j))
    !enddo ; enddo

  enddo ! a-loop

  ! Only reflect newly arrived energy; existing energy in incident wedge
  ! is not reflected and will eventually propagate out of cell.
  ! (only reflects if En > 0)
  call reflect(fdt_m(:,:,:), nangle, cs, g, lb)
  call teleport(fdt_m(:,:,:), nangle, cs, g, lb)
  call reflect(fdt_p(:,:,:), nangle, cs, g, lb)
  call teleport(fdt_p(:,:,:), nangle, cs, g, lb)

  ! Update reflected energy (Jm-2)
  do j=jsh,jeh ; do i=ish,ieh
    !do a=1,CS%nAngle
    !  if((En(i,j,a) + G%IareaT(i,j)*(Fdt_m(i,j,a) + Fdt_p(i,j,a))) < 0.0)then ! for debugging
    !    print *,"propagate_x: OutFlux>Available" ; !stop
    !  endif
    !enddo
    en(i,j,:) = en(i,j,:) + g%IareaT(i,j)*(fdt_m(i,j,:) + fdt_p(i,j,:))
  enddo ; enddo

end subroutine propagate_x

! #@# This subroutine needs a doxygen description.
subroutine propagate_y(En, speed_y, Cgy_av, dCgy, dt, G, Nangle, CS, LB)
  type(ocean_grid_type),   intent(in)    :: G  !< The ocean's grid structure.
  integer,                 intent(in)    :: NAngle
  real, dimension(G%isd:G%ied,G%jsd:G%jed,Nangle), &
                           intent(inout) :: En !< The energy density integrated over an angular
                                               !! band, in J m-2, intent in/out.
  real, dimension(G%isd:G%ied,G%JsdB:G%JedB),      &
                           intent(in)    :: speed_y !< The magnitude of the group velocity at the
                                               !! Cv points, in m s-1.
  real, dimension(Nangle), intent(in)    :: Cgy_av, dCgy
  real,                    intent(in)    :: dt !< Time increment in s.
  type(int_tide_cs),       pointer       :: CS !< The control structure returned by a previous call
                                               !! to continuity_PPM_init.
  type(loop_bounds_type),  intent(in)    :: LB !< A structure with the active energy loop bounds.

  ! Arguments: En - The energy density integrated over an angular band, in J m-2,
  !                 intent in/out.
  !  (in)      speed_y - The magnitude of the group velocity at the Cv
  !                      points, in m s-1.
  !  (in)      dt - Time increment in s.
  !  (in)      G - The ocean's grid structure.
  !  (in)      CS - The control structure returned by a previous call to
  !                 continuity_PPM_init.
  !  (in)      LB - A structure with the active energy loop bounds.
  real, dimension(SZI_(G),SZJ_(G)) :: &
    EnL, EnR    ! South and north face energy densities, in J m-2.
  real, dimension(SZI_(G),SZJB_(G)) :: &
    flux_y      ! The internal wave energy flux, in J s-1.
  real, dimension(SZI_(G)) :: &
    cg_p, cg_m, flux1, flux2
  !real, dimension(SZI_(G),SZJB_(G),Nangle) :: En_m, En_p
  real, dimension(SZI_(G),SZJB_(G),Nangle) :: &
    Fdt_m, Fdt_p! South and north energy fluxes, in J
  integer :: i, j, k, ish, ieh, jsh, jeh, a

  ish = lb%ish ; ieh = lb%ieh ; jsh = lb%jsh ; jeh = lb%jeh
  do a=1,nangle
    ! This sets EnL and EnR.
    if (cs%upwind_1st) then
      do j=jsh-1,jeh+1 ; do i=ish,ieh
        enl(i,j) = en(i,j,a) ; enr(i,j) = en(i,j,a)
      enddo ; enddo
    else
      call ppm_reconstruction_y(en(:,:,a), enl, enr, g, lb, simple_2nd=cs%simple_2nd)
    endif

    do j=jsh-1,jeh
      !   This is done once with single speed (GARDEN SPRINKLER EFFECT POSSIBLE)
      do i=ish,ieh
        cg_p(i) = speed_y(i,j) * (cgy_av(a))
      enddo
      call merid_flux_en(cg_p, en(:,:,a), enl(:,:), enr(:,:), flux1, &
                         dt, g, j, ish, ieh, cs%vol_CFL)
      do i=ish,ieh ; flux_y(i,j) = flux1(i); enddo
    enddo

    do j=jsh,jeh ; do i=ish,ieh
      fdt_m(i,j,a) = dt*flux_y(i,j-1) ! south face influx (J)
      fdt_p(i,j,a) = -dt*flux_y(i,j)  ! north face influx (J)
      !if((En(i,j,a) + G%IareaT(i,j)*(Fdt_m(i,j,a) + Fdt_p(i,j,a))) < 0.0)then ! for debugging
      !  print *,"propagate_y: OutFlux>Available prior to reflection" ; !stop
      !  print *,"flux_y_south=",flux_y(i,J-1)
      !  print *,"flux_y_north=",flux_y(i,J)
      !  print *,"En=",En(i,j,a)
      !  print *,"cn_south=", speed_y(i,J-1) * (Cgy_av(a))
      !  print *,"cn_north=", speed_y(i,J) * (Cgy_av(a))
      !endif
    enddo ; enddo

    ! test with old (take out later)
    !do j=jsh,jeh ; do i=ish,ieh
    ! En(i,j,a) = En(i,j,a) - dt* G%IareaT(i,j) * (flux_y(i,J) - flux_y(i,J-1))
    !enddo ; enddo

  enddo ! a-loop

  ! Only reflect newly arrived energy; existing energy in incident wedge
  ! is not reflected and will eventually propagate out of cell.
  ! (only reflects if En > 0)
  call reflect(fdt_m(:,:,:), nangle, cs, g, lb)
  call teleport(fdt_m(:,:,:), nangle, cs, g, lb)
  call reflect(fdt_p(:,:,:), nangle, cs, g, lb)
  call teleport(fdt_p(:,:,:), nangle, cs, g, lb)

  ! Update reflected energy (Jm-2)
  do j=jsh,jeh ; do i=ish,ieh
    !do a=1,CS%nAngle
    !  if((En(i,j,a) + G%IareaT(i,j)*(Fdt_m(i,j,a) + Fdt_p(i,j,a))) < 0.0)then ! for debugging
    !    print *,"propagate_y: OutFlux>Available" ; !stop
    !  endif
    !enddo
    en(i,j,:) = en(i,j,:) + g%IareaT(i,j)*(fdt_m(i,j,:) + fdt_p(i,j,:))
  enddo ; enddo

end subroutine propagate_y

!> This subroutines evaluates the zonal mass or volume fluxes in a layer.
subroutine zonal_flux_en(u, h, hL, hR, uh, dt, G, j, ish, ieh, vol_CFL)
  type(ocean_grid_type),     intent(in)    :: G  !< The ocean's grid structure.
  real, dimension(SZIB_(G)), intent(in)    :: u  !< The zonal velocity, in m s-1.
  real, dimension(SZI_(G)),  intent(in)    :: h  !< Energy density used to calculate the fluxes,
                                                 !! in J m-2.
  real, dimension(SZI_(G)),  intent(in)    :: hL !< Left- Energy densities in the reconstruction,
                                                 !! in J m-2.
  real, dimension(SZI_(G)),  intent(in)    :: hR !< Right- Energy densities in the reconstruction,
                                                 !! in J m-2.
  real, dimension(SZIB_(G)), intent(inout) :: uh !< The zonal energy transport,
                                                 !! in J s-1.
  real,                      intent(in)    :: dt !< Time increment in s.
  integer,                   intent(in)    :: j, ish, ieh !< The index range to work on.
  logical,                   intent(in)    :: vol_CFL !< If true, rescale the ratio of face areas to
                                                 !! the cell areas when estimating the CFL number.

  !   This subroutines evaluates the zonal mass or volume fluxes in a layer.
  !
  ! Arguments: u - Zonal velocity, in m s-1.
  !  (in)      h - Energy density used to calculate the fluxes, in J m-2.
  !  (in)      hL, hR - Left- and right- Energy densities in the reconstruction, in J m-2.
  !  (out)     uh - The zonal energy transport, in J s-1.
  !  (in)      dt - Time increment in s.
  !  (in)      G - The ocean's grid structure.
  !  (in)      j, ish, ieh - The index range to work on.
  !  (in)      vol_CFL - If true, rescale the ratio of face areas to the cell
  !                      areas when estimating the CFL number.
  real :: CFL  ! The CFL number based on the local velocity and grid spacing, ND.
  real :: curv_3 ! A measure of the thickness curvature over a grid length,
                 ! with the same units as h_in.
  integer :: i

  do i=ish-1,ieh
    ! Set new values of uh and duhdu.
    if (u(i) > 0.0) then
      if (vol_cfl) then ; cfl = (u(i) * dt) * (g%dy_Cu(i,j) * g%IareaT(i,j))
      else ; cfl = u(i) * dt * g%IdxT(i,j) ; endif
      curv_3 = (hl(i) + hr(i)) - 2.0*h(i)
      uh(i) = g%dy_Cu(i,j) * u(i) * &
          (hr(i) + cfl * (0.5*(hl(i) - hr(i)) + curv_3*(cfl - 1.5)))
    elseif (u(i) < 0.0) then
      if (vol_cfl) then ; cfl = (-u(i) * dt) * (g%dy_Cu(i,j) * g%IareaT(i+1,j))
      else ; cfl = -u(i) * dt * g%IdxT(i+1,j) ; endif
      curv_3 = (hl(i+1) + hr(i+1)) - 2.0*h(i+1)
      uh(i) = g%dy_Cu(i,j) * u(i) * &
          (hl(i+1) + cfl * (0.5*(hr(i+1)-hl(i+1)) + curv_3*(cfl - 1.5)))
    else
      uh(i) = 0.0
    endif
  enddo
end subroutine zonal_flux_en

!> This subroutines evaluates the meridional mass or volume fluxes in a layer.
subroutine merid_flux_en(v, h, hL, hR, vh, dt, G, J, ish, ieh, vol_CFL)
  type(ocean_grid_type),            intent(in)    :: G  !< The ocean's grid structure.
  real, dimension(SZI_(G)),         intent(in)    :: v  !< The meridional velocity, in m s-1.
  real, dimension(SZI_(G),SZJ_(G)), intent(in)    :: h  !< Energy density used to calculate the
                                                        !! fluxes, in J m-2.
  real, dimension(SZI_(G),SZJ_(G)), intent(in)    :: hL !< Left- Energy densities in the
                                                        !! reconstruction, in J m-2.
  real, dimension(SZI_(G),SZJ_(G)), intent(in)    :: hR !< Right- Energy densities in the
                                                        !! reconstruction, in J m-2.
  real, dimension(SZI_(G)),         intent(inout) :: vh !< The meridional energy transport,
                                                        !! in J s-1.
  real,                             intent(in)    :: dt !< Time increment in s.
  integer,                          intent(in)    :: J, ish, ieh !< The index range to work on.
  logical,                          intent(in)    :: vol_CFL !< If true, rescale the ratio of face
                                                        !! areas to the cell areas when estimating
                                                        !! the CFL number.
  !   This subroutines evaluates the meridional mass or volume fluxes in a layer.
  !
  ! Arguments: v - Meridional velocity, in m s-1.
  !  (in)      h - Energy density used to calculate the fluxes, in J m-2.
  !  (in)      hL, hR - Left- and right- Energy densities in the reconstruction, in J m-2.
  !  (out)     vh - The meridional energy transport, in J s-1.
  !  (in)      dt - Time increment in s.
  !  (in)      G - The ocean's grid structure.
  !  (in)      J, ish, ieh - The index range to work on.
  !  (in)      vol_CFL - If true, rescale the ratio of face areas to the cell
  !                       areas when estimating the CFL number.
  real :: CFL ! The CFL number based on the local velocity and grid spacing, ND.
  real :: curv_3 ! A measure of the thickness curvature over a grid length,
                 ! with the same units as h_in.
  integer :: i

  do i=ish,ieh
    if (v(i) > 0.0) then
      if (vol_cfl) then ; cfl = (v(i) * dt) * (g%dx_Cv(i,j) * g%IareaT(i,j))
      else ; cfl = v(i) * dt * g%IdyT(i,j) ; endif
      curv_3 = hl(i,j) + hr(i,j) - 2.0*h(i,j)
      vh(i) = g%dx_Cv(i,j) * v(i) * ( hr(i,j) + cfl * &
          (0.5*(hl(i,j) - hr(i,j)) + curv_3*(cfl - 1.5)) )
    elseif (v(i) < 0.0) then
      if (vol_cfl) then ; cfl = (-v(i) * dt) * (g%dx_Cv(i,j) * g%IareaT(i,j+1))
      else ; cfl = -v(i) * dt * g%IdyT(i,j+1) ; endif
      curv_3 = hl(i,j+1) + hr(i,j+1) - 2.0*h(i,j+1)
      vh(i) = g%dx_Cv(i,j) * v(i) * ( hl(i,j+1) + cfl * &
          (0.5*(hr(i,j+1)-hl(i,j+1)) + curv_3*(cfl - 1.5)) )
    else
      vh(i) = 0.0
    endif
  enddo
end subroutine merid_flux_en

!> This subroutine does reflection of the internal waves at a single frequency.
subroutine reflect(En, NAngle, CS, G, LB)
  type(ocean_grid_type),  intent(in)    :: G    !< The ocean's grid structure
  integer,                intent(in)    :: NAngle
  real, dimension(G%isd:G%ied,G%jsd:G%jed,NAngle), &
                          intent(inout) :: En
  type(int_tide_cs),      pointer       :: CS
  type(loop_bounds_type), intent(in)    :: LB

  !  This subroutine does reflection of the internal waves at a single frequency.

  real, dimension(G%isd:G%ied,G%jsd:G%jed) :: angle_c
                                           ! angle of boudary wrt equator
  real, dimension(G%isd:G%ied,G%jsd:G%jed) :: part_refl
                                           ! fraction of wave energy reflected
                                           ! values should collocate with angle_c
  logical, dimension(G%isd:G%ied,G%jsd:G%jed) :: ridge
                                           ! tags of cells with double reflection

  real    :: TwoPi                         ! 2*pi
  real    :: Angle_size                    ! size of beam wedge (rad)
  real    :: angle_wall                    ! angle of coast/ridge/shelf wrt equator
  real, dimension(1:NAngle) :: angle_i     ! angle of incident ray wrt equator
  real    :: angle_r                       ! angle of reflected ray wrt equator
  real, dimension(1:Nangle) :: En_reflected
  integer :: i, j, a, a_r, na
  !integer :: isd, ied, jsd, jed   ! start and end local indices on data domain
  !                                ! (values include halos)
  integer :: isc, iec, jsc, jec   ! start and end local indices on PE
                                  ! (values exclude halos)
  integer :: ish, ieh, jsh, jeh   ! start and end local indices on data domain
                                  ! leaving out outdated halo points (march in)
  integer :: id_g, jd_g           ! global (decomp-invar) indices

  !isd = G%isd  ; ied = G%ied  ; jsd = G%jsd  ; jed = G%jed
  isc = g%isc  ; iec = g%iec  ; jsc = g%jsc  ; jec = g%jec
  ish = lb%ish ; ieh = lb%ieh ; jsh = lb%jsh ; jeh = lb%jeh

  twopi = 8.0*atan(1.0);
  angle_size = twopi / (real(nangle))

  do a=1,nangle
    ! These are the angles at the cell centers
    ! (should do this elsewhere since doesn't change with time)
    angle_i(a) = angle_size * real(a - 1) ! for a=1 aligned with x-axis
  enddo

  angle_c   = cs%refl_angle
  part_refl = cs%refl_pref
  ridge     = cs%refl_dbl
  en_reflected(:) = 0.0

  !do j=jsc-1,jec+1
  do j=jsh,jeh
    jd_g = j + g%jdg_offset
    !do i=isc-1,iec+1
    do i=ish,ieh
      id_g = i + g%idg_offset
      ! redistribute energy in angular space if ray will hit boundary
      ! i.e., if energy is in a reflecting cell
      if (angle_c(i,j) .ne. cs%nullangle) then
        do a=1,nangle
          if (en(i,j,a) > 0.0) then
            ! if ray is incident, keep specified boundary angle
            if (sin(angle_i(a) - angle_c(i,j)) >= 0.0) then
              angle_wall = angle_c(i,j)
            ! if ray is not incident but in ridge cell, use complementary angle
            elseif (ridge(i,j)) then
              angle_wall = angle_c(i,j) + 0.5*twopi
              if (angle_wall > twopi) then
                angle_wall = angle_wall - twopi*floor(abs(angle_wall)/twopi)
              elseif (angle_wall < 0.0) then
                angle_wall = angle_wall + twopi*ceiling(abs(angle_wall)/twopi)
              endif
            ! if ray is not incident and not in a ridge cell, keep specified angle
            else
              angle_wall = angle_c(i,j)
            endif
            ! do reflection
            if (sin(angle_i(a) - angle_wall) >= 0.0) then
              angle_r = 2.0*angle_wall - angle_i(a)
              if (angle_r > twopi) then
                angle_r = angle_r - twopi*floor(abs(angle_r)/twopi)
              elseif (angle_r < 0.0) then
                angle_r = angle_r + twopi*ceiling(abs(angle_r)/twopi)
              endif
              a_r = nint(angle_r/angle_size) + 1
              do while (a_r > nangle) ; a_r = a_r - nangle ; enddo
              if (a .ne. a_r) then
                en_reflected(a_r) = part_refl(i,j)*en(i,j,a)
                en(i,j,a)   = (1.0-part_refl(i,j))*en(i,j,a)
              endif
            endif
          endif
        enddo ! a-loop
        en(i,j,:) = en(i,j,:) + en_reflected(:)
        en_reflected(:) = 0.0
      endif
    enddo ! i-loop
  enddo ! j-loop

  ! Check to make sure no energy gets onto land (only run for debugging)
  !do j=jsc,jec
  !  jd_g = j + G%jdg_offset
  !  do i=isc,iec
  !    id_g = i + G%idg_offset
  !    do a=1,NAngle
  !      if (En(i,j,a) > 0.001 .and. G%mask2dT(i,j) == 0) then
  !        print *, 'En=', En(i,j,a), 'a=', a, 'ig_g=',id_g, 'jg_g=',jd_g
  !        !stop 'Energy detected out of bounds!'
  !      endif
  !    enddo ! a-loop
  !  enddo ! i-loop
  !enddo ! j-loop

end subroutine reflect

!> This subroutine moves energy across lines of partial reflection to prevent
!! reflection of energy that is supposed to get across.
subroutine teleport(En, NAngle, CS, G, LB)
  type(ocean_grid_type),  intent(in)    :: G    !< The ocean's grid structure.
  integer,                intent(in)    :: NAngle
  real, dimension(G%isd:G%ied,G%jsd:G%jed,NAngle), &
                          intent(inout) :: En
  type(int_tide_cs),      pointer       :: CS
  type(loop_bounds_type), intent(in)    :: LB

  ! This subroutine moves energy across lines of partial reflection to prevent
  ! reflection of energy that is supposed to get across.

  real, dimension(G%isd:G%ied,G%jsd:G%jed)    :: angle_c
                                              ! angle of boudary wrt equator
  real, dimension(G%isd:G%ied,G%jsd:G%jed)    :: part_refl
                                              ! fraction of wave energy reflected
                                              ! values should collocate with angle_c
  logical, dimension(G%isd:G%ied,G%jsd:G%jed) :: pref_cell
                                              ! flag for partial reflection
  logical, dimension(G%isd:G%ied,G%jsd:G%jed) :: ridge
                                           ! tags of cells with double reflection
  real                        :: TwoPi      ! size of beam wedge (rad)
  real                        :: Angle_size ! size of beam wedge (rad)
  real, dimension(1:NAngle)   :: angle_i    ! angle of incident ray wrt equator
  real, dimension(1:NAngle)   :: cos_angle, sin_angle
  real                        :: En_tele    ! energy to be "teleported"
  integer :: i, j, a
  !integer :: isd, ied, jsd, jed    ! start and end local indices on data domain
  !                                 ! (values include halos)
  !integer :: isc, iec, jsc, jec    ! start and end local indices on PE
  !                                 ! (values exclude halos)
  integer :: ish, ieh, jsh, jeh     ! start and end local indices on data domain
                                    ! leaving out outdated halo points (march in)
  integer :: id_g, jd_g             ! global (decomp-invar) indices
  integer :: jos, ios               ! offsets
  real    :: cos_normal, sin_normal, angle_wall
                                    ! cos/sin of cross-ridge normal, ridge angle

  !isd = G%isd  ; ied = G%ied  ; jsd = G%jsd  ; jed = G%jed
  !isc = G%isc  ; iec = G%iec  ; jsc = G%jsc  ; jec = G%jec
  ish = lb%ish ; ieh = lb%ieh ; jsh = lb%jsh ; jeh = lb%jeh

  twopi = 8.0*atan(1.0)
  angle_size = twopi / (real(nangle))

  do a=1,nangle
    ! These are the angles at the cell centers
    ! (should do this elsewhere since doesn't change with time)
    angle_i(a) = angle_size * real(a - 1) ! for a=1 aligned with x-axis
    cos_angle(a) = cos(angle_i(a)) ; sin_angle(a) = sin(angle_i(a))
  enddo

  angle_c   = cs%refl_angle
  part_refl = cs%refl_pref
  pref_cell = cs%refl_pref_logical
  ridge     = cs%refl_dbl

  do j=jsh,jeh
    do i=ish,ieh
      id_g = i + g%idg_offset ; jd_g = j + g%jdg_offset
      if (pref_cell(i,j)) then
        do a=1,nangle
          if (en(i,j,a) > 0) then
            ! if ray is incident, keep specified boundary angle
            if (sin(angle_i(a) - angle_c(i,j)) >= 0.0) then
              angle_wall = angle_c(i,j)
            ! if ray is not incident but in ridge cell, use complementary angle
            elseif (ridge(i,j)) then
              angle_wall = angle_c(i,j) + 0.5*twopi
            ! if ray is not incident and not in a ridge cell, keep specified angle
            else
              angle_wall = angle_c(i,j)
            endif
            ! teleport if incident
            if (sin(angle_i(a) - angle_wall) >= 0.0) then
              en_tele = en(i,j,a)
              cos_normal = cos(angle_wall + 0.25*twopi)
              sin_normal = sin(angle_wall + 0.25*twopi)
              ! find preferred zonal offset based on shelf/ridge angle
              ios = int(sign(1.,cos_normal))
              ! find preferred meridional offset based on shelf/ridge angle
              jos = int(sign(1.,sin_normal))
              ! find receptive ocean cell in direction of offset
              if (.not. pref_cell(i+ios,j+jos)) then
                en(i,j,a) = en(i,j,a) - en_tele
                en(i+ios,j+jos,a) = en(i+ios,j+jos,a) + en_tele
              else
                call mom_error(warning, "teleport: no receptive ocean cell", .true.)
                print *, 'idg=',id_g,'jd_g=',jd_g,'a=',a
                stop
              endif
            endif ! incidence check
          endif ! energy check
        enddo ! a-loop
      endif ! pref check
    enddo ! i-loop
  enddo ! j-loop

end subroutine teleport

!> This subroutine rotates points in the halos where required to accomodate
!! changes in grid orientation, such as at the tripolar fold.
subroutine correct_halo_rotation(En, test, G, NAngle)
  type(ocean_grid_type),              intent(in)    :: G    !< The ocean's grid structure
  real, dimension(:,:,:,:,:),         intent(inout) :: En
  real, dimension(SZI_(G),SZJ_(G),2), intent(in)    :: test
  integer,                            intent(in)    :: NAngle
  !   This subroutine rotates points in the halos where required to accomodate
  ! changes in grid orientation, such as at the tripolar fold.

  real, dimension(G%isd:G%ied,NAngle) :: En2d
  integer, dimension(G%isd:G%ied) :: a_shift
  integer :: i_first, i_last, a_new
  integer :: a, i, j, isd, ied, jsd, jed, m, fr
  character(len=80) :: mesg
  isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed

  do j=jsd,jed
    i_first = ied+1 ; i_last = isd-1
    do i=isd,ied
      a_shift(i) = 0
      if (test(i,j,1) /= 1.0) then
        if (i<i_first) i_first = i
        if (i>i_last) i_last = i

        if (test(i,j,1) == -1.0) then ; a_shift(i) = nangle/2
        elseif (test(i,j,2) == 1.0) then ; a_shift(i) = -nangle/4
        elseif (test(i,j,2) == -1.0) then ; a_shift(i) = nangle/4
        else
          write(mesg,'("Unrecognized rotation test vector ",2ES9.2," at ",F7.2," E, ",&
                       &F7.2," N; i,j=",2i4)') &
                test(i,j,1), test(i,j,2), g%GeoLonT(i,j), g%GeoLatT(i,j), i, j
          call mom_error(fatal, mesg)
        endif
      endif
    enddo

    if (i_first <= i_last) then
      ! At least one point in this row needs to be rotated.
      do m=1,size(en,5) ; do fr=1,size(en,4)
        do a=1,nangle ; do i=i_first,i_last ; if (a_shift(i) /= 0) then
          a_new = a + a_shift(i)
          if (a_new < 1) a_new = a_new + nangle
          if (a_new > nangle) a_new = a_new - nangle
          en2d(i,a_new) = en(i,j,a,fr,m)
        endif ; enddo ; enddo
        do a=1,nangle ; do i=i_first,i_last ; if (a_shift(i) /= 0) then
          en(i,j,a,fr,m) = en2d(i,a)
        endif ; enddo ; enddo
      enddo ; enddo
    endif
  enddo
end subroutine correct_halo_rotation

!> This subroutine calculates left/right edge values for PPM reconstruction.
subroutine ppm_reconstruction_x(h_in, h_l, h_r, G, LB, simple_2nd)
  type(ocean_grid_type),            intent(in)  :: G    !< The ocean's grid structure.
  real, dimension(SZI_(G),SZJ_(G)), intent(in)  :: h_in !< Energy density in a sector (2D).
  real, dimension(SZI_(G),SZJ_(G)), intent(out) :: h_l  !< Left edge value of reconstruction (2D).
  real, dimension(SZI_(G),SZJ_(G)), intent(out) :: h_r  !< Right edge value of reconstruction (2D).
  type(loop_bounds_type),           intent(in)  :: LB   !< A structure with the active loop bounds.
  logical, optional,                intent(in)  :: simple_2nd !< If true, use the arithmetic mean
                                                        !! energy densities as default edge values
                                                        !! for a simple 2nd order scheme.

  ! This subroutine calculates left/right edge values for PPM reconstruction.
  ! Arguments: h_in    - Energy density in a sector (2D)
  !  (out)     h_l,h_r - left/right edge value of reconstruction (2D)
  !  (in)      G - The ocean's grid structure.
  !  (in)      LB - A structure with the active loop bounds.
  !  (in, opt) simple_2nd - If true, use the arithmetic mean energy densities as
  !                         default edge values for a simple 2nd order scheme.

  ! Local variables with useful mnemonic names.
  real, dimension(SZI_(G),SZJ_(G))  :: slp ! The slopes.
  real, parameter :: oneSixth = 1./6.
  real :: h_ip1, h_im1
  real :: dMx, dMn
  logical :: use_CW84, use_2nd
  character(len=256) :: mesg
  integer :: i, j, isl, iel, jsl, jel, stencil

  use_2nd = .false. ; if (present(simple_2nd)) use_2nd = simple_2nd
  isl = lb%ish-1 ; iel = lb%ieh+1 ; jsl = lb%jsh ; jel = lb%jeh

  ! This is the stencil of the reconstruction, not the scheme overall.
  stencil = 2 ; if (use_2nd) stencil = 1

  if ((isl-stencil < g%isd) .or. (iel+stencil > g%ied)) then
    write(mesg,'("In MOM_internal_tides, PPM_reconstruction_x called with a ", &
               & "x-halo that needs to be increased by ",i2,".")') &
               stencil + max(g%isd-isl,iel-g%ied)
    call mom_error(fatal,mesg)
  endif
  if ((jsl < g%jsd) .or. (jel > g%jed)) then
    write(mesg,'("In MOM_internal_tides, PPM_reconstruction_x called with a ", &
               & "y-halo that needs to be increased by ",i2,".")') &
               max(g%jsd-jsl,jel-g%jed)
    call mom_error(fatal,mesg)
  endif

  if (use_2nd) then
    do j=jsl,jel ; do i=isl,iel
      h_im1 = g%mask2dT(i-1,j) * h_in(i-1,j) + (1.0-g%mask2dT(i-1,j)) * h_in(i,j)
      h_ip1 = g%mask2dT(i+1,j) * h_in(i+1,j) + (1.0-g%mask2dT(i+1,j)) * h_in(i,j)
      h_l(i,j) = 0.5*( h_im1 + h_in(i,j) )
      h_r(i,j) = 0.5*( h_ip1 + h_in(i,j) )
    enddo ; enddo
  else
    do j=jsl,jel ; do i=isl-1,iel+1
      if ((g%mask2dT(i-1,j) * g%mask2dT(i,j) * g%mask2dT(i+1,j)) == 0.0) then
        slp(i,j) = 0.0
      else
        ! This uses a simple 2nd order slope.
        slp(i,j) = 0.5 * (h_in(i+1,j) - h_in(i-1,j))
        ! Monotonic constraint, see Eq. B2 in Lin 1994, MWR (132)
        dmx = max(h_in(i+1,j), h_in(i-1,j), h_in(i,j)) - h_in(i,j)
        dmn = h_in(i,j) - min(h_in(i+1,j), h_in(i-1,j), h_in(i,j))
        slp(i,j) = sign(1.,slp(i,j)) * min(abs(slp(i,j)), 2. * min(dmx, dmn))
                ! * (G%mask2dT(i-1,j) * G%mask2dT(i,j) * G%mask2dT(i+1,j))
      endif
    enddo; enddo

    do j=jsl,jel ; do i=isl,iel
      ! Neighboring values should take into account any boundaries.  The 3
      ! following sets of expressions are equivalent.
    ! h_im1 = h_in(i-1,j,k) ; if (G%mask2dT(i-1,j) < 0.5) h_im1 = h_in(i,j)
    ! h_ip1 = h_in(i+1,j,k) ; if (G%mask2dT(i+1,j) < 0.5) h_ip1 = h_in(i,j)
      h_im1 = g%mask2dT(i-1,j) * h_in(i-1,j) + (1.0-g%mask2dT(i-1,j)) * h_in(i,j)
      h_ip1 = g%mask2dT(i+1,j) * h_in(i+1,j) + (1.0-g%mask2dT(i+1,j)) * h_in(i,j)
      ! Left/right values following Eq. B2 in Lin 1994, MWR (132)
      h_l(i,j) = 0.5*( h_im1 + h_in(i,j) ) + onesixth*( slp(i-1,j) - slp(i,j) )
      h_r(i,j) = 0.5*( h_ip1 + h_in(i,j) ) + onesixth*( slp(i,j) - slp(i+1,j) )
    enddo; enddo
  endif

  call ppm_limit_pos(h_in, h_l, h_r, 0.0, g, isl, iel, jsl, jel)
end subroutine ppm_reconstruction_x

!> This subroutine calculates left/right edge valus for PPM reconstruction.
subroutine ppm_reconstruction_y(h_in, h_l, h_r, G, LB, simple_2nd)
  type(ocean_grid_type),            intent(in)  :: G    !< The ocean's grid structure.
  real, dimension(SZI_(G),SZJ_(G)), intent(in)  :: h_in !< Energy density in a sector (2D).
  real, dimension(SZI_(G),SZJ_(G)), intent(out) :: h_l  !< Left edge value of reconstruction (2D).
  real, dimension(SZI_(G),SZJ_(G)), intent(out) :: h_r  !< Right edge value of reconstruction (2D).
  type(loop_bounds_type),           intent(in)  :: LB   !< A structure with the active loop bounds.
  logical, optional,                intent(in)  :: simple_2nd !< If true, use the arithmetic mean
                                                        !! energy densities as default edge values
                                                        !! for a simple 2nd order scheme.

  ! This subroutine calculates left/right edge valus for PPM reconstruction.
  ! Arguments: h_in    - Energy density in a sector (2D)
  !  (out)     h_l,h_r - left/right edge value of reconstruction (2D)
  !  (in)      G - The ocean's grid structure.
  !  (in)      LB - A structure with the active loop bounds.
  !  (in, opt) simple_2nd - If true, use the arithmetic mean energy densities as
  !                         default edge values for a simple 2nd order scheme.

  ! Local variables with useful mnemonic names.
  real, dimension(SZI_(G),SZJ_(G))  :: slp ! The slopes.
  real, parameter :: oneSixth = 1./6.
  real :: h_jp1, h_jm1
  real :: dMx, dMn
  logical :: use_2nd
  character(len=256) :: mesg
  integer :: i, j, isl, iel, jsl, jel, stencil

  use_2nd = .false. ; if (present(simple_2nd)) use_2nd = simple_2nd
  isl = lb%ish ; iel = lb%ieh ; jsl = lb%jsh-1 ; jel = lb%jeh+1

  ! This is the stencil of the reconstruction, not the scheme overall.
  stencil = 2 ; if (use_2nd) stencil = 1

  if ((isl < g%isd) .or. (iel > g%ied)) then
    write(mesg,'("In MOM_internal_tides, PPM_reconstruction_y called with a ", &
               & "x-halo that needs to be increased by ",i2,".")') &
               max(g%isd-isl,iel-g%ied)
    call mom_error(fatal,mesg)
  endif
  if ((jsl-stencil < g%jsd) .or. (jel+stencil > g%jed)) then
    write(mesg,'("In MOM_internal_tides, PPM_reconstruction_y called with a ", &
                 & "y-halo that needs to be increased by ",i2,".")') &
                 stencil + max(g%jsd-jsl,jel-g%jed)
    call mom_error(fatal,mesg)
  endif

  if (use_2nd) then
    do j=jsl,jel ; do i=isl,iel
      h_jm1 = g%mask2dT(i,j-1) * h_in(i,j-1) + (1.0-g%mask2dT(i,j-1)) * h_in(i,j)
      h_jp1 = g%mask2dT(i,j+1) * h_in(i,j+1) + (1.0-g%mask2dT(i,j+1)) * h_in(i,j)
      h_l(i,j) = 0.5*( h_jm1 + h_in(i,j) )
      h_r(i,j) = 0.5*( h_jp1 + h_in(i,j) )
    enddo ; enddo
  else
    do j=jsl-1,jel+1 ; do i=isl,iel
      if ((g%mask2dT(i,j-1) * g%mask2dT(i,j) * g%mask2dT(i,j+1)) == 0.0) then
        slp(i,j) = 0.0
      else
        ! This uses a simple 2nd order slope.
        slp(i,j) = 0.5 * (h_in(i,j+1) - h_in(i,j-1))
        ! Monotonic constraint, see Eq. B2 in Lin 1994, MWR (132)
        dmx = max(h_in(i,j+1), h_in(i,j-1), h_in(i,j)) - h_in(i,j)
        dmn = h_in(i,j) - min(h_in(i,j+1), h_in(i,j-1), h_in(i,j))
        slp(i,j) = sign(1.,slp(i,j)) * min(abs(slp(i,j)), 2. * min(dmx, dmn))
                ! * (G%mask2dT(i,j-1) * G%mask2dT(i,j) * G%mask2dT(i,j+1))
      endif
    enddo ; enddo

    do j=jsl,jel ; do i=isl,iel
      ! Neighboring values should take into account any boundaries.  The 3
      ! following sets of expressions are equivalent.
      h_jm1 = g%mask2dT(i,j-1) * h_in(i,j-1) + (1.0-g%mask2dT(i,j-1)) * h_in(i,j)
      h_jp1 = g%mask2dT(i,j+1) * h_in(i,j+1) + (1.0-g%mask2dT(i,j+1)) * h_in(i,j)
      ! Left/right values following Eq. B2 in Lin 1994, MWR (132)
      h_l(i,j) = 0.5*( h_jm1 + h_in(i,j) ) + onesixth*( slp(i,j-1) - slp(i,j) )
      h_r(i,j) = 0.5*( h_jp1 + h_in(i,j) ) + onesixth*( slp(i,j) - slp(i,j+1) )
    enddo ; enddo
  endif

  call ppm_limit_pos(h_in, h_l, h_r, 0.0, g, isl, iel, jsl, jel)
end subroutine ppm_reconstruction_y

!> This subroutine limits the left/right edge values of the PPM reconstruction
!! to give a reconstruction that is positive-definite.  Here this is
!! reinterpreted as giving a constant thickness if the mean thickness is less
!! than h_min, with a minimum of h_min otherwise.
subroutine ppm_limit_pos(h_in, h_L, h_R, h_min, G, iis, iie, jis, jie)
  type(ocean_grid_type),            intent(in)     :: G     !< The ocean's grid structure.
  real, dimension(SZI_(G),SZJ_(G)), intent(in)     :: h_in  !< Thickness of layer (2D).
  real, dimension(SZI_(G),SZJ_(G)), intent(inout)  :: h_L   !< Left edge value (2D).
  real, dimension(SZI_(G),SZJ_(G)), intent(inout)  :: h_R   !< Right edge value (2D).
  real,                             intent(in)     :: h_min !< The minimum thickness that can be
                                                            !! obtained by a concave parabolic fit.
  integer,                          intent(in)     :: iis, iie, jis, jie !< Index range for
                                                            !! computation.

  ! This subroutine limits the left/right edge values of the PPM reconstruction
  ! to give a reconstruction that is positive-definite.  Here this is
  ! reinterpreted as giving a constant thickness if the mean thickness is less
  ! than h_min, with a minimum of h_min otherwise.
  ! Arguments: h_in    - thickness of layer (2D)
  !  (inout)   h_L     - left edge value (2D)
  !  (inout)   h_R     - right edge value (2D)
  !  (in)      h_min   - The minimum thickness that can be obtained by a
  !                      concave parabolic fit.
  !  (in)      iis, iie, jis, jie - Index range for computation.
  !  (in)      G - The ocean's grid structure.

  ! Local variables
  real    :: curv, dh, scale
  character(len=256) :: mesg
  integer :: i,j

  do j=jis,jie ; do i=iis,iie
    ! This limiter prevents undershooting minima within the domain with
    ! values less than h_min.
    curv = 3.0*(h_l(i,j) + h_r(i,j) - 2.0*h_in(i,j))
    if (curv > 0.0) then ! Only minima are limited.
      dh = h_r(i,j) - h_l(i,j)
      if (abs(dh) < curv) then ! The parabola's minimum is within the cell.
        if (h_in(i,j) <= h_min) then
          h_l(i,j) = h_in(i,j) ; h_r(i,j) = h_in(i,j)
        elseif (12.0*curv*(h_in(i,j) - h_min) < (curv**2 + 3.0*dh**2)) then
          ! The minimum value is h_in - (curv^2 + 3*dh^2)/(12*curv), and must
          ! be limited in this case.  0 < scale < 1.
          scale = 12.0*curv*(h_in(i,j) - h_min) / (curv**2 + 3.0*dh**2)
          h_l(i,j) = h_in(i,j) + scale*(h_l(i,j) - h_in(i,j))
          h_r(i,j) = h_in(i,j) + scale*(h_r(i,j) - h_in(i,j))
        endif
      endif
    endif
  enddo ; enddo
end subroutine ppm_limit_pos

! subroutine register_int_tide_restarts(G, param_file, CS, restart_CS)
!   type(ocean_grid_type), intent(inout) :: G    !< The ocean's grid structure
!   type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters
!   type(int_tide_CS),     pointer :: CS
!   type(MOM_restart_CS),  pointer :: restart_CS

!   ! This subroutine is not currently in use!!

!   ! Arguments: G - The ocean's grid structure.
!   !  (in)      param_file - A structure indicating the open file to parse for
!   !                         model parameter values.
!   !  (in/out)  CS - A pointer that is set to point to the control structure
!   !                 for this module.
!   !  (in)      restart_CS - A pointer to the restart control structure.
!   ! This subroutine is used to allocate and register any fields in this module
!   ! that should be written to or read from the restart file.
!   logical :: use_int_tides
!   type(vardesc) :: vd
!   integer :: num_freq, num_angle , num_mode, period_1
!   integer :: isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB, a
!   isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed
!   IsdB = G%IsdB ; IedB = G%IedB ; JsdB = G%JsdB ; JedB = G%JedB

!   if (associated(CS)) then
!     call MOM_error(WARNING, "register_int_tide_restarts called "//&
!                              "with an associated control structure.")
!     return
!   endif

!   use_int_tides = .false.
!   call read_param(param_file, "INTERNAL_TIDES", use_int_tides)
!   if (.not.use_int_tides) return

!   allocate(CS)

!   num_angle = 24
!   call read_param(param_file, "INTERNAL_TIDE_ANGLES", num_angle)
!   allocate(CS%En_restart(isd:ied, jsd:jed, num_angle))
!   CS%En_restart(:,:,:) = 0.0

!   vd = vardesc("En_restart", &
!     "The internal wave energy density as a function of (i,j,angle,frequency,mode)", &
!     'h','1','1',"J m-2")
!   call register_restart_field(CS%En_restart, vd, .false., restart_CS)

!   !--------------------check----------------------------------------------
!   if (is_root_pe()) then
!     print *,'register_int_tide_restarts: CS and CS%En_restart allocated!'
!     print *,'register_int_tide_restarts: CS%En_restart registered!'
!     print *,'register_int_tide_restarts: done!'
!   endif
!   !-----------------------------------------------------------------------

! end subroutine register_int_tide_restarts

! #@# This subroutine needs a doxygen comment.
subroutine internal_tides_init(Time, G, GV, param_file, diag, CS)
  type(time_type), target,   intent(in)    :: Time !< The current model time.
  type(ocean_grid_type),     intent(inout) :: G    !< The ocean's grid structure.
  type(verticalgrid_type),   intent(in)    :: GV   !< The ocean's vertical grid structure.
  type(param_file_type),     intent(in)    :: param_file !< A structure to parse for run-time
                                                   !! parameters.
  type(diag_ctrl), target,   intent(in)    :: diag !< A structure that is used to regulate
                                                   !! diagnostic output.
  type(int_tide_cs),pointer                :: CS   !< A pointer that is set to point to the control
                                                   !! structure for this module.

  ! Arguments: Time - The current model time.
  !  (in)      G - The ocean's grid structure.
  !  (in)      GV - The ocean's vertical grid structure.
  !  (in)      param_file - A structure indicating the open file to parse for
  !                         model parameter values.
  !  (in)      diag - A structure that is used to regulate diagnostic output.
  !  (in/out)  CS - A pointer that is set to point to the control structure
  !                 for this module
  real                              :: Angle_size ! size of wedges, rad
  real, allocatable                 :: angles(:)  ! orientations of wedge centers, rad
  real, allocatable, dimension(:,:) :: h2         ! topographic roughness scale, m^2
  real                              :: kappa_itides, kappa_h2_factor
                                                  ! characteristic topographic wave number
                                                  ! and a scaling factor
  real, allocatable                 :: ridge_temp(:,:)
                                                  ! array for temporary storage of flags
                                                  ! of cells with double-reflecting ridges
  logical                           :: use_int_tides, use_temperature
  integer                           :: num_angle, num_freq, num_mode, m, fr, period_1
  integer :: isd, ied, jsd, jed, a, id_ang, i, j
  type(axes_grp) :: axes_ang
  ! This include declares and sets the variable "version".
#include "version_variable.h"
  character(len=40)  :: mdl = "MOM_internal_tides" ! This module's name.
  character(len=16), dimension(8) :: freq_name
  character(len=40)  :: var_name
  character(len=160) :: var_descript
  character(len=200) :: filename
  character(len=200) :: refl_angle_file, land_mask_file
  character(len=200) :: refl_pref_file, refl_dbl_file
  character(len=200) :: dy_Cu_file, dx_Cv_file
  character(len=200) :: h2_file

  isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed

  if (associated(cs)) then
    call mom_error(warning, "internal_tides_init called "//&
                             "with an associated control structure.")
    return
  else
    allocate(cs)
  endif

  use_int_tides = .false.
  call read_param(param_file, "INTERNAL_TIDES", use_int_tides)
  cs%do_int_tides = use_int_tides
  if (.not.use_int_tides) return

  use_temperature = .true.
  call read_param(param_file, "ENABLE_THERMODYNAMICS", use_temperature)
  if (.not.use_temperature) call mom_error(fatal, &
    "register_int_tide_restarts: internal_tides only works with "//&
    "ENABLE_THERMODYNAMICS defined.")

  ! Set number of frequencies, angles, and modes to consider
  num_freq = 1 ; num_angle = 24 ; num_mode = 1
  call read_param(param_file, "INTERNAL_TIDE_FREQS", num_freq)
  call read_param(param_file, "INTERNAL_TIDE_ANGLES", num_angle)
  call read_param(param_file, "INTERNAL_TIDE_MODES", num_mode)
  if (.not.((num_freq > 0) .and. (num_angle > 0) .and. (num_mode > 0))) return
  cs%nFreq = num_freq ; cs%nAngle = num_angle ; cs%nMode = num_mode

  ! Allocate energy density array
  allocate(cs%En(isd:ied, jsd:jed, num_angle, num_freq, num_mode))
  cs%En(:,:,:,:,:) = 0.0

  ! Allocate phase speed array
  allocate(cs%cp(isd:ied, jsd:jed, num_freq, num_mode))
  cs%cp(:,:,:,:) = 0.0

  ! Allocate and populate frequency array (each a multiple of first for now)
  allocate(cs%frequency(num_freq))
  call read_param(param_file, "FIRST_MODE_PERIOD", period_1); ! ADDED BDM
  do fr=1,num_freq
    cs%frequency(fr) = (8.0*atan(1.0) * (real(fr)) / period_1) ! ADDED BDM
  enddo

  ! Read all relevant parameters and write them to the model log.

  cs%Time => time ! direct a pointer to the current model time target

  call get_param(param_file, mdl, "INPUTDIR", cs%inputdir, default=".")
                 cs%inputdir = slasher(cs%inputdir)

  call log_version(param_file, mdl, version, "")

  call get_param(param_file, mdl, "INTERNAL_TIDE_FREQS", num_freq, &
                 "The number of distinct internal tide frequency bands \n"//&
                 "that will be calculated.", default=1)
  call get_param(param_file, mdl, "INTERNAL_TIDE_MODES", num_mode, &
                 "The number of distinct internal tide modes \n"//&
                 "that will be calculated.", default=1)
  call get_param(param_file, mdl, "INTERNAL_TIDE_ANGLES", num_angle, &
                 "The number of angular resolution bands for the internal \n"//&
                 "tide calculations.", default=24)

  if (use_int_tides) then
    if ((num_freq <= 0) .and. (num_mode <= 0) .and. (num_angle <= 0)) then
      call mom_error(warning, "Internal tides were enabled, but the number "//&
             "of requested frequencies, modes and angles were not all positive.")
      return
    endif
  else
    if ((num_freq > 0) .and. (num_mode > 0) .and. (num_angle > 0)) then
      call mom_error(warning, "Internal tides were not enabled, even though "//&
             "the number of requested frequencies, modes and angles were all "//&
             "positive.")
      return
    endif
  endif

  if (cs%NFreq /= num_freq) call mom_error(fatal, "Internal_tides_init: "//&
      "Inconsistent number of frequencies.")
  if (cs%NAngle /= num_angle) call mom_error(fatal, "Internal_tides_init: "//&
      "Inconsistent number of angles.")
  if (cs%NMode /= num_mode) call mom_error(fatal, "Internal_tides_init: "//&
      "Inconsistent number of modes.")
  if (4*(num_angle/4) /= num_angle) call mom_error(fatal, &
    "Internal_tides_init: INTERNAL_TIDE_ANGLES must be a multiple of 4.")

  cs%diag => diag

  call get_param(param_file, mdl, "INTERNAL_TIDE_DECAY_RATE", cs%decay_rate, &
                 "The rate at which internal tide energy is lost to the \n"//&
                 "interior ocean internal wave field.", units="s-1", default=0.0)
  call get_param(param_file, mdl, "INTERNAL_TIDE_VOLUME_BASED_CFL", cs%vol_CFL, &
                 "If true, use the ratio of the open face lengths to the \n"//&
                 "tracer cell areas when estimating CFL numbers in the \n"//&
                 "internal tide code.", default=.false.)
  call get_param(param_file, mdl, "INTERNAL_TIDE_CORNER_ADVECT", cs%corner_adv, &
                 "If true, internal tide ray-tracing advection uses a \n"//&
                 " corner-advection scheme rather than PPM.\n", default=.false.)
  call get_param(param_file, mdl, "INTERNAL_TIDE_SIMPLE_2ND_PPM", cs%simple_2nd, &
                 "If true, CONTINUITY_PPM uses a simple 2nd order \n"//&
                 "(arithmetic mean) interpolation of the edge values. \n"//&
                 "This may give better PV conservation propterties. While \n"//&
                 "it formally reduces the accuracy of the continuity \n"//&
                 "solver itself in the strongly advective limit, it does \n"//&
                 "not reduce the overall order of accuracy of the dynamic \n"//&
                 "core.", default=.false.)
  call get_param(param_file, mdl, "INTERNAL_TIDE_UPWIND_1ST", cs%upwind_1st, &
                 "If true, the internal tide ray-tracing advection uses \n"//&
                 "1st-order upwind advection.  This scheme is highly \n"//&
                 "continuity solver.  This scheme is highly \n"//&
                 "diffusive but may be useful for debugging.", default=.false.)
  call get_param(param_file, mdl, "INTERNAL_TIDE_BACKGROUND_DRAG", &
                 cs%apply_background_drag, "If true, the internal tide \n"//&
                 "ray-tracing advection uses a background drag term as a sink.",&
                 default=.false.)
  call get_param(param_file, mdl, "INTERNAL_TIDE_QUAD_DRAG", cs%apply_bottom_drag, &
                 "If true, the internal tide ray-tracing advection uses \n"//&
                 "a quadratic bottom drag term as a sink.", default=.false.)
  call get_param(param_file, mdl, "INTERNAL_TIDE_WAVE_DRAG", cs%apply_wave_drag, &
                 "If true, apply scattering due to small-scale roughness as a sink.", &
                 default=.false.)
  call get_param(param_file, mdl, "INTERNAL_TIDE_FROUDE_DRAG", cs%apply_Froude_drag, &
                 "If true, apply wave breaking as a sink.", &
                 default=.false.)
  call get_param(param_file, mdl, "CDRAG", cs%cdrag, &
                 "CDRAG is the drag coefficient relating the magnitude of \n"//&
                 "the velocity field to the bottom stress.", units="nondim", &
                 default=0.003)
  call get_param(param_file, mdl, "INTERNAL_TIDE_ENERGIZED_ANGLE", cs%energized_angle, &
                 "If positive, only one angular band of the internal tides \n"//&
                 "gets all of the energy.  (This is for debugging.)", default=-1)
  call get_param(param_file, mdl, "USE_PPM_ANGULAR", cs%use_PPMang, &
                 "If true, use PPM for advection of energy in angular  \n"//&
                 "space.", default=.false.)
  call get_param(param_file, mdl, "GAMMA_ITIDES", cs%q_itides, &
                 "The fraction of the internal tidal energy that is \n"//&
                 "dissipated locally with INT_TIDE_DISSIPATION.  \n"//&
                 "THIS NAME COULD BE BETTER.", &
                 units="nondim", default=0.3333)
    call get_param(param_file, mdl, "KAPPA_ITIDES", kappa_itides, &
               "A topographic wavenumber used with INT_TIDE_DISSIPATION. \n"//&
               "The default is 2pi/10 km, as in St.Laurent et al. 2002.", &
               units="m-1", default=8.e-4*atan(1.0))
  call get_param(param_file, mdl, "KAPPA_H2_FACTOR", kappa_h2_factor, &
               "A scaling factor for the roughness amplitude with n"//&
               "INT_TIDE_DISSIPATION.",  units="nondim", default=1.0)

  ! Allocate various arrays needed for loss rates
  allocate(h2(isd:ied,jsd:jed)) ; h2(:,:) = 0.0
  allocate(cs%TKE_itidal_loss_fixed(isd:ied,jsd:jed))
    cs%TKE_itidal_loss_fixed = 0.0
  allocate(cs%TKE_leak_loss(isd:ied,jsd:jed,num_angle,num_freq,num_mode))
    cs%TKE_leak_loss(:,:,:,:,:) = 0.0
  allocate(cs%TKE_quad_loss(isd:ied,jsd:jed,num_angle,num_freq,num_mode))
    cs%TKE_quad_loss(:,:,:,:,:) = 0.0
  allocate(cs%TKE_itidal_loss(isd:ied,jsd:jed,num_angle,num_freq,num_mode))
    cs%TKE_itidal_loss(:,:,:,:,:) = 0.0
  allocate(cs%TKE_Froude_loss(isd:ied,jsd:jed,num_angle,num_freq,num_mode))
    cs%TKE_Froude_loss(:,:,:,:,:) = 0.0
  allocate(cs%tot_leak_loss(isd:ied,jsd:jed))   ; cs%tot_leak_loss(:,:)   = 0.0
  allocate(cs%tot_quad_loss(isd:ied,jsd:jed) )  ; cs%tot_quad_loss(:,:)   = 0.0
  allocate(cs%tot_itidal_loss(isd:ied,jsd:jed)) ; cs%tot_itidal_loss(:,:) = 0.0
  allocate(cs%tot_Froude_loss(isd:ied,jsd:jed)) ; cs%tot_Froude_loss(:,:) = 0.0

  ! Compute the fixed part of the bottom drag loss from baroclinic modes
  call get_param(param_file, mdl, "H2_FILE", h2_file, &
          "The path to the file containing the sub-grid-scale \n"//&
          "topographic roughness amplitude with INT_TIDE_DISSIPATION.", &
          fail_if_missing=.true.)
  filename = trim(cs%inputdir) // trim(h2_file)
  call log_param(param_file, mdl, "INPUTDIR/H2_FILE", filename)
  call read_data(filename, 'h2', h2, domain=g%domain%mpp_domain, timelevel=1)
  do j=g%jsc,g%jec ; do i=g%isc,g%iec
    ! Restrict rms topo to 10 percent of column depth.
    h2(i,j) = min(0.01*g%bathyT(i,j)**2, h2(i,j))
    ! Compute the fixed part; units are [kg m-2] here;
    ! will be multiplied by N and En to get into [W m-2]
    cs%TKE_itidal_loss_fixed(i,j) = 0.5*kappa_h2_factor*gv%Rho0*&
         kappa_itides * h2(i,j)
  enddo; enddo

  ! Read in prescribed coast/ridge/shelf angles from file
  call get_param(param_file, mdl, "REFL_ANGLE_FILE", refl_angle_file, &
               "The path to the file containing the local angle of \n"//&
               "the coastline/ridge/shelf with respect to the equator.", &
               fail_if_missing=.false.)
  filename = trim(cs%inputdir) // trim(refl_angle_file)
  call log_param(param_file, mdl, "INPUTDIR/REFL_ANGLE_FILE", filename)
  allocate(cs%refl_angle(isd:ied,jsd:jed)) ; cs%refl_angle(:,:) = cs%nullangle
  call read_data(filename, 'refl_angle', cs%refl_angle, &
                 domain=g%domain%mpp_domain, timelevel=1)
  ! replace NANs with null value
  do j=g%jsc,g%jec ; do i=g%isc,g%iec
    if(is_nan(cs%refl_angle(i,j))) cs%refl_angle(i,j) = cs%nullangle
  enddo ; enddo
  call pass_var(cs%refl_angle,g%domain)

  ! Read in prescribed partial reflection coefficients from file
  call get_param(param_file, mdl, "REFL_PREF_FILE", refl_pref_file, &
               "The path to the file containing the reflection coefficients.", &
               fail_if_missing=.false.)
  filename = trim(cs%inputdir) // trim(refl_pref_file)
  call log_param(param_file, mdl, "INPUTDIR/REFL_PREF_FILE", filename)
  allocate(cs%refl_pref(isd:ied,jsd:jed)) ; cs%refl_pref(:,:) = 1.0
  call read_data(filename, 'refl_pref', cs%refl_pref, &
                 domain=g%domain%mpp_domain, timelevel=1)
  !CS%refl_pref = CS%refl_pref*1 ! adjust partial reflection if desired
  call pass_var(cs%refl_pref,g%domain)

  ! Tag reflection cells with partial reflection (done here for speed)
  allocate(cs%refl_pref_logical(isd:ied,jsd:jed)) ; cs%refl_pref_logical(:,:) = .false.
  do j=jsd,jed
    do i=isd,ied
      ! flag cells with partial reflection
      if (cs%refl_angle(i,j) .ne. cs%nullangle .and. &
        cs%refl_pref(i,j) < 1.0 .and. cs%refl_pref(i,j) > 0.0) then
        cs%refl_pref_logical(i,j) = .true.
      endif
    enddo
  enddo

  ! Read in double-reflective (ridge) tags from file
  call get_param(param_file, mdl, "REFL_DBL_FILE", refl_dbl_file, &
               "The path to the file containing the double-reflective ridge tags.", &
               fail_if_missing=.false.)
  filename = trim(cs%inputdir) // trim(refl_dbl_file)
  call log_param(param_file, mdl, "INPUTDIR/REFL_DBL_FILE", filename)
  allocate(ridge_temp(isd:ied,jsd:jed)) ; ridge_temp(:,:) = 0.0
  call read_data(filename, 'refl_dbl', ridge_temp, &
                 domain=g%domain%mpp_domain, timelevel=1)
  call pass_var(ridge_temp,g%domain)
  allocate(cs%refl_dbl(isd:ied,jsd:jed)) ; cs%refl_dbl(:,:) = .false.
  do i=isd,ied; do j=jsd,jed
    if (ridge_temp(i,j) == 1) then; cs%refl_dbl(i,j) = .true.
    else ; cs%refl_dbl(i,j) = .false. ; endif
  enddo; enddo

  ! Read in prescribed land mask from file (if overwriting -BDM).
  ! This should be done in MOM_initialize_topography subroutine
  ! defined in MOM_fixed_initialization.F90 (BDM)
  !call get_param(param_file, mdl, "LAND_MASK_FILE", land_mask_file, &
  !             "The path to the file containing the land mask.", &
  !             fail_if_missing=.false.)
  !filename = trim(CS%inputdir) // trim(land_mask_file)
  !call log_param(param_file, mdl, "INPUTDIR/LAND_MASK_FILE", filename)
  !G%mask2dCu(:,:) = 1 ; G%mask2dCv(:,:) = 1 ; G%mask2dT(:,:)  = 1
  !call read_data(filename, 'land_mask', G%mask2dCu, &
  !               domain=G%domain%mpp_domain, timelevel=1)
  !call read_data(filename, 'land_mask', G%mask2dCv, &
  !               domain=G%domain%mpp_domain, timelevel=1)
  !call read_data(filename, 'land_mask', G%mask2dT, &
  !               domain=G%domain%mpp_domain, timelevel=1)
  !call pass_var(G%mask2dCu,G%domain)
  !call pass_var(G%mask2dCv,G%domain)
  !call pass_var(G%mask2dT,G%domain)

  ! Read in prescribed partial east face blockages from file (if overwriting -BDM)
  !call get_param(param_file, mdl, "dy_Cu_FILE", dy_Cu_file, &
  !             "The path to the file containing the east face blockages.", &
  !             fail_if_missing=.false.)
  !filename = trim(CS%inputdir) // trim(dy_Cu_file)
  !call log_param(param_file, mdl, "INPUTDIR/dy_Cu_FILE", filename)
  !G%dy_Cu(:,:) = 0.0
  !call read_data(filename, 'dy_Cu', G%dy_Cu, &
  !               domain=G%domain%mpp_domain, timelevel=1)
  !call pass_var(G%dy_Cu,G%domain)

  ! Read in prescribed partial north face blockages from file (if overwriting -BDM)
  !call get_param(param_file, mdl, "dx_Cv_FILE", dx_Cv_file, &
  !             "The path to the file containing the north face blockages.", &
  !             fail_if_missing=.false.)
  !filename = trim(CS%inputdir) // trim(dx_Cv_file)
  !call log_param(param_file, mdl, "INPUTDIR/dx_Cv_FILE", filename)
  !G%dx_Cv(:,:) = 0.0
  !call read_data(filename, 'dx_Cv', G%dx_Cv, &
  !               domain=G%domain%mpp_domain, timelevel=1)
  !call pass_var(G%dx_Cv,G%domain)

  ! For debugging - delete later
  call get_param(param_file, mdl, "INTERNAL_TIDE_SOURCE_X", cs%int_tide_source_x, &
                "X Location of generation site for internal tide", default=1.)
  call get_param(param_file, mdl, "INTERNAL_TIDE_SOURCE_Y", cs%int_tide_source_y, &
                "Y Location of generation site for internal tide", default=1.)

  ! Register maps of reflection parameters
  cs%id_refl_ang = register_diag_field('ocean_model', 'refl_angle', diag%axesT1, &
                 time, 'Local angle of coastline/ridge/shelf with respect to equator', 'rad')
  cs%id_refl_pref = register_diag_field('ocean_model', 'refl_pref', diag%axesT1, &
                 time, 'Partial reflection coefficients', '')
  cs%id_dx_Cv = register_diag_field('ocean_model', 'dx_Cv', diag%axesT1, &
                 time, 'North face unblocked width', 'm') ! used if overriding (BDM)
  cs%id_dy_Cu = register_diag_field('ocean_model', 'dy_Cu', diag%axesT1, &
                 time, 'East face unblocked width', 'm')  ! used if overriding (BDM)
  cs%id_land_mask = register_diag_field('ocean_model', 'land_mask', diag%axesT1, &
                 time, 'Land mask', 'logical')            ! used if overriding (BDM)
  ! Output reflection parameters as diags here (not needed every timestep)
  if (cs%id_refl_ang > 0)   call post_data(cs%id_refl_ang, cs%refl_angle, cs%diag)
  if (cs%id_refl_pref > 0)  call post_data(cs%id_refl_pref, cs%refl_pref, cs%diag)
  if (cs%id_dx_Cv > 0)      call post_data(cs%id_dx_Cv, g%dx_Cv, cs%diag)
  if (cs%id_dy_Cu > 0)      call post_data(cs%id_dy_Cu, g%dy_Cu, cs%diag)
  if (cs%id_land_mask > 0)  call post_data(cs%id_land_mask, g%mask2dT, cs%diag)

  ! Register 2-D energy density (summed over angles, freq, modes)
  cs%id_tot_En = register_diag_field('ocean_model', 'ITide_tot_En', diag%axesT1, &
                 time, 'Internal tide total energy density', 'J m-2')
  ! Register 2-D drag scale used for quadratic bottom drag
  cs%id_itide_drag = register_diag_field('ocean_model', 'ITide_drag', diag%axesT1, &
                 time, 'Interior and bottom drag internal tide decay timescale', 's-1')
  !Register 2-D energy input into internal tides
  cs%id_TKE_itidal_input = register_diag_field('ocean_model', 'TKE_itidal_input', diag%axesT1, &
                 time, 'Conversion from barotropic to baroclinic tide, \n'//&
                 'a fraction of which goes into rays', 'W m-2')
  ! Register 2-D energy losses (summed over angles, freq, modes)
  cs%id_tot_leak_loss = register_diag_field('ocean_model', 'ITide_tot_leak_loss', diag%axesT1, &
                time, 'Internal tide energy loss to background drag', 'W m-2')
  cs%id_tot_quad_loss = register_diag_field('ocean_model', 'ITide_tot_quad_loss', diag%axesT1, &
                time, 'Internal tide energy loss to bottom drag', 'W m-2')
  cs%id_tot_itidal_loss = register_diag_field('ocean_model', 'ITide_tot_itidal_loss', diag%axesT1, &
                time, 'Internal tide energy loss to wave drag', 'W m-2')
  cs%id_tot_Froude_loss = register_diag_field('ocean_model', 'ITide_tot_Froude_loss', diag%axesT1, &
                time, 'Internal tide energy loss to wave breaking', 'W m-2')
  cs%id_tot_allprocesses_loss = register_diag_field('ocean_model', 'ITide_tot_allprocesses_loss', diag%axesT1, &
                time, 'Internal tide energy loss summed over all processes', 'W m-2')

  allocate(cs%id_En_mode(cs%nFreq,cs%nMode)) ; cs%id_En_mode(:,:) = -1
  allocate(cs%id_En_ang_mode(cs%nFreq,cs%nMode)) ; cs%id_En_ang_mode(:,:) = -1
  allocate(cs%id_itidal_loss_mode(cs%nFreq,cs%nMode)) ; cs%id_itidal_loss_mode(:,:) = -1
  allocate(cs%id_allprocesses_loss_mode(cs%nFreq,cs%nMode)) ; cs%id_allprocesses_loss_mode(:,:) = -1
  allocate(cs%id_itidal_loss_ang_mode(cs%nFreq,cs%nMode)) ; cs%id_itidal_loss_ang_mode(:,:) = -1
  allocate(cs%id_Ub_mode(cs%nFreq,cs%nMode)) ; cs%id_Ub_mode(:,:) = -1
  allocate(cs%id_cp_mode(cs%nFreq,cs%nMode)) ; cs%id_cp_mode(:,:) = -1

  allocate(angles(cs%NAngle)) ; angles(:) = 0.0
  angle_size = (8.0*atan(1.0)) / (real(num_angle))
  do a=1,num_angle ; angles(a) = (real(a) - 1) * angle_size ; enddo

  id_ang = diag_axis_init("angle", angles, "Radians", "N", "Angular Orienation of Fluxes")
  call define_axes_group(diag, (/ diag%axesT1%handles(1), diag%axesT1%handles(2), id_ang /), axes_ang)

  do fr=1,cs%nFreq ; write(freq_name(fr), '("freq",i1)') fr ; enddo
  do m=1,cs%nMode ; do fr=1,cs%nFreq
    ! Register 2-D energy density (summed over angles) for each freq and mode
    write(var_name, '("Itide_En_freq",i1,"_mode",i1)') fr, m
    write(var_descript, '("Internal tide energy density in frequency ",i1," mode ",i1)') fr, m
    cs%id_En_mode(fr,m) = register_diag_field('ocean_model', var_name, &
                 diag%axesT1, time, var_descript, 'J m-2')
    call mom_mesg("Registering "//trim(var_name)//", Described as: "//var_descript, 5)

    ! Register 3-D (i,j,a) energy density for each freq and mode
    write(var_name, '("Itide_En_ang_freq",i1,"_mode",i1)') fr, m
    write(var_descript, '("Internal tide angular energy density in frequency ",i1," mode ",i1)') fr, m
    cs%id_En_ang_mode(fr,m) = register_diag_field('ocean_model', var_name, &
                 axes_ang, time, var_descript, 'J m-2 band-1')
    call mom_mesg("Registering "//trim(var_name)//", Described as: "//var_descript, 5)

    ! Register 2-D energy loss (summed over angles) for each freq and mode
    ! wave-drag only
    write(var_name, '("Itide_wavedrag_loss_freq",i1,"_mode",i1)') fr, m
    write(var_descript, '("Internal tide energy loss due to wave-drag from frequency ",i1," mode ",i1)') fr, m
    cs%id_itidal_loss_mode(fr,m) = register_diag_field('ocean_model', var_name, &
                 diag%axesT1, time, var_descript, 'W m-2')
    call mom_mesg("Registering "//trim(var_name)//", Described as: "//var_descript, 5)
    ! all loss processes
    write(var_name, '("Itide_allprocesses_loss_freq",i1,"_mode",i1)') fr, m
    write(var_descript, '("Internal tide energy loss due to all processes from frequency ",i1," mode ",i1)') fr, m
    cs%id_allprocesses_loss_mode(fr,m) = register_diag_field('ocean_model', var_name, &
                 diag%axesT1, time, var_descript, 'W m-2')
    call mom_mesg("Registering "//trim(var_name)//", Described as: "//var_descript, 5)

    ! Register 3-D (i,j,a) energy loss for each freq and mode
    ! wave-drag only
    write(var_name, '("Itide_wavedrag_loss_ang_freq",i1,"_mode",i1)') fr, m
    write(var_descript, '("Internal tide energy loss due to wave-drag from frequency ",i1," mode ",i1)') fr, m
    cs%id_itidal_loss_ang_mode(fr,m) = register_diag_field('ocean_model', var_name, &
                 axes_ang, time, var_descript, 'W m-2 band-1')
    call mom_mesg("Registering "//trim(var_name)//", Described as: "//var_descript, 5)

    ! Register 2-D period-averaged near-bottom horizonal velocity for each freq and mode
    write(var_name, '("Itide_Ub_freq",i1,"_mode",i1)') fr, m
    write(var_descript, '("Near-bottom horizonal velocity for frequency ",i1," mode ",i1)') fr, m
    cs%id_Ub_mode(fr,m) = register_diag_field('ocean_model', var_name, &
                 diag%axesT1, time, var_descript, 'm s-1')
    call mom_mesg("Registering "//trim(var_name)//", Described as: "//var_descript, 5)

    ! Register 2-D horizonal phase velocity for each freq and mode
    write(var_name, '("Itide_cp_freq",i1,"_mode",i1)') fr, m
    write(var_descript, '("Horizonal phase velocity for frequency ",i1," mode ",i1)') fr, m
    cs%id_cp_mode(fr,m) = register_diag_field('ocean_model', var_name, &
                 diag%axesT1, time, var_descript, 'm s-1')
    call mom_mesg("Registering "//trim(var_name)//", Described as: "//var_descript, 5)

  enddo ; enddo

  ! Initialize wave_structure (not sure if this should be here - BDM)
  call wave_structure_init(time, g, param_file, diag, cs%wave_structure_CSp)

end subroutine internal_tides_init


subroutine internal_tides_end(CS)
  type(int_tide_cs),            pointer :: CS
  ! Arguments:  CS - A pointer to the control structure returned by a previous
  !                  call to internal_tides_init, it will be deallocated here.

  if (associated(cs)) then
    if (associated(cs%En)) deallocate(cs%En)
    if (allocated(cs%frequency)) deallocate(cs%frequency)
    if (allocated(cs%id_En_mode)) deallocate(cs%id_En_mode)
    if (allocated(cs%id_Ub_mode)) deallocate(cs%id_Ub_mode)
    if (allocated(cs%id_cp_mode)) deallocate(cs%id_cp_mode)
    deallocate(cs)
  endif
  cs => null()
end subroutine internal_tides_end


end module mom_internal_tides