MOM_wave_speed.F90

!> Routines for calculating baroclinic wave speeds
module mom_wave_speed

! This file is part of MOM6. See LICENSE.md for the license.

use mom_diag_mediator, only : post_data, query_averaging_enabled, diag_ctrl
use mom_error_handler, only : mom_error, fatal, warning
use mom_file_parser, only : log_version
use mom_grid, only : ocean_grid_type
use mom_remapping, only : remapping_cs, initialize_remapping, remapping_core_h
use mom_unit_scaling, only : unit_scale_type
use mom_variables, only : thermo_var_ptrs
use mom_verticalgrid, only : verticalgrid_type
use mom_eos, only : calculate_density, calculate_density_derivs

implicit none ; private

#include <MOM_memory.h>

public wave_speed, wave_speeds, wave_speed_init, wave_speed_set_param

! A note on unit descriptions in comments: MOM6 uses units that can be rescaled for dimensional
! consistency testing. These are noted in comments with units like Z, H, L, and T, along with
! their mks counterparts with notation like "a velocity [Z T-1 ~> m s-1]".  If the units
! vary with the Boussinesq approximation, the Boussinesq variant is given first.

!> Control structure for MOM_wave_speed
type, public :: wave_speed_cs ; private
  logical :: use_ebt_mode = .false.    !< If true, calculate the equivalent barotropic wave speed instead
                                       !! of the first baroclinic wave speed.
                                       !! This parameter controls the default behavior of wave_speed() which
                                       !! can be overridden by optional arguments.
  logical :: better_cg1_est = .false.  !< If true, use an improved estimate of the first mode
                                       !! internal wave speed.
  real :: mono_n2_column_fraction = 0. !< The lower fraction of water column over which N2 is limited as
                                       !! monotonic for the purposes of calculating the equivalent barotropic
                                       !! wave speed. This parameter controls the default behavior of
                                       !! wave_speed() which can be overridden by optional arguments.
  real :: mono_n2_depth = -1.          !< The depth below which N2 is limited as monotonic for the purposes of
                                       !! calculating the equivalent barotropic wave speed [Z ~> m].
                                       !! This parameter controls the default behavior of wave_speed() which
                                       !! can be overridden by optional arguments.
  real :: min_speed2 = 0.              !< The minimum mode 1 internal wave speed squared [L2 T-2 ~> m2 s-2]
  real :: wave_speed_tol = 0.001       !< The fractional tolerance with which to solve for the wave
                                       !! speeds [nondim]
  type(remapping_cs) :: remapping_cs   !< Used for vertical remapping when calculating equivalent barotropic
                                       !! mode structure.
  logical :: remap_answers_2018 = .true.  !< If true, use the order of arithmetic and expressions that
                                       !! recover the remapping answers from 2018.  If false, use more
                                       !! robust forms of the same remapping expressions.
  type(diag_ctrl), pointer :: diag     !< Diagnostics control structure
end type wave_speed_cs

contains

!> Calculates the wave speed of the first baroclinic mode.
subroutine wave_speed(h, tv, G, GV, US, cg1, CS, full_halos, use_ebt_mode, mono_N2_column_fraction, &
                      mono_N2_depth, modal_structure, better_speed_est, min_speed, wave_speed_tol)
  type(ocean_grid_type),            intent(in)  :: G  !< Ocean grid structure
  type(verticalgrid_type),          intent(in)  :: GV !< Vertical grid structure
  type(unit_scale_type),            intent(in)  :: US !< A dimensional unit scaling type
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
                                    intent(in)  :: h  !< Layer thickness [H ~> m or kg m-2]
  type(thermo_var_ptrs),            intent(in)  :: tv !< Thermodynamic variables
  real, dimension(SZI_(G),SZJ_(G)), intent(out) :: cg1 !< First mode internal wave speed [L T-1 ~> m s-1]
  type(wave_speed_cs),              pointer     :: CS !< Control structure for MOM_wave_speed
  logical,                optional, intent(in)  :: full_halos !< If true, do the calculation
                                          !! over the entire computational domain.
  logical,                optional, intent(in)  :: use_ebt_mode !< If true, use the equivalent
                                          !! barotropic mode instead of the first baroclinic mode.
  real,                   optional, intent(in)  :: mono_N2_column_fraction !< The lower fraction
                                          !! of water column over which N2 is limited as monotonic
                                          !! for the purposes of calculating vertical modal structure.
  real,                   optional, intent(in)  :: mono_N2_depth !< A depth below which N2 is limited as
                                          !! monotonic for the purposes of calculating vertical
                                          !! modal structure [Z ~> m].
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
                          optional, intent(out) :: modal_structure !< Normalized model structure [nondim]
  logical, optional, intent(in) :: better_speed_est !< If true, use a more robust estimate of the first
                                     !! mode speed as the starting point for iterations.
  real,    optional, intent(in) :: min_speed !< If present, set a floor in the first mode speed
                                     !! below which 0 is returned [L T-1 ~> m s-1].
  real,    optional, intent(in) :: wave_speed_tol !< The fractional tolerance for finding the
                                     !! wave speeds [nondim]

  ! Local variables
  real, dimension(SZK_(G)+1) :: &
    dRho_dT, &    ! Partial derivative of density with temperature [R degC-1 ~> kg m-3 degC-1]
    dRho_dS, &    ! Partial derivative of density with salinity [R ppt-1 ~> kg m-3 ppt-1]
    pres, &       ! Interface pressure [R L2 T-2 ~> Pa]
    T_int, &      ! Temperature interpolated to interfaces [degC]
    S_int, &      ! Salinity interpolated to interfaces [ppt]
    H_top, &      ! The distance of each filtered interface from the ocean surface [Z ~> m]
    H_bot, &      ! The distance of each filtered interface from the bottom [Z ~> m]
    gprime        ! The reduced gravity across each interface [L2 Z-1 T-2 ~> m s-2].
  real, dimension(SZK_(G)) :: &
    Igl, Igu      ! The inverse of the reduced gravity across an interface times
                  ! the thickness of the layer below (Igl) or above (Igu) it, in [T2 L-2 ~> s2 m-2].
  real, dimension(SZK_(G),SZI_(G)) :: &
    Hf, &         ! Layer thicknesses after very thin layers are combined [Z ~> m]
    Tf, &         ! Layer temperatures after very thin layers are combined [degC]
    Sf, &         ! Layer salinities after very thin layers are combined [ppt]
    Rf            ! Layer densities after very thin layers are combined [R ~> kg m-3]
  real, dimension(SZK_(G)) :: &
    Hc, &         ! A column of layer thicknesses after convective istabilities are removed [Z ~> m]
    Tc, &         ! A column of layer temperatures after convective istabilities are removed [degC]
    Sc, &         ! A column of layer salinites after convective istabilities are removed [ppt]
    Rc, &         ! A column of layer densities after convective istabilities are removed [R ~> kg m-3]
    Hc_H          ! Hc(:) rescaled from Z to thickness units [H ~> m or kg m-2]
  real :: I_Htot  ! The inverse of the total filtered thicknesses [Z ~> m]
  real :: det, ddet, detKm1, detKm2, ddetKm1, ddetKm2
  real :: lam     ! The eigenvalue [T2 L-2 ~> s m-1]
  real :: dlam    ! The change in estimates of the eigenvalue [T2 L-2 ~> s m-1]
  real :: lam0    ! The first guess of the eigenvalue [T2 L-2 ~> s m-1]
  real :: min_h_frac ! [nondim]
  real :: Z_to_pres  ! A conversion factor from thicknesses to pressure [R L2 T-2 Z-1 ~> Pa m-1]
  real, dimension(SZI_(G)) :: &
    htot, hmin, &  ! Thicknesses [Z ~> m]
    H_here, &      ! A thickness [Z ~> m]
    HxT_here, &    ! A layer integrated temperature [degC Z ~> degC m]
    HxS_here, &    ! A layer integrated salinity [ppt Z ~> ppt m]
    HxR_here       ! A layer integrated density [R Z ~> kg m-2]
  real :: speed2_tot ! overestimate of the mode-1 speed squared [L2 T-2 ~> m2 s-2]
  real :: cg1_min2 ! A floor in the squared first mode speed below which 0 is returned [L2 T-2 ~> m2 s-2]
  real :: I_Hnew   ! The inverse of a new layer thickness [Z-1 ~> m-1]
  real :: drxh_sum ! The sum of density differences across interfaces times thicknesses [R Z ~> kg m-2]
  real :: L2_to_Z2 ! A scaling factor squared from units of lateral distances to depths [Z2 L-2 ~> 1].
  real, pointer, dimension(:,:,:) :: T => null(), s => null()
  real :: g_Rho0   ! G_Earth/Rho0 [L2 T-2 Z-1 R-1 ~> m4 s-2 kg-1].
  real :: c2_scale ! A scaling factor for wave speeds to help control the growth of the determinant
                   ! and its derivative with lam between rows of the Thomas algorithm solver.  The
                   ! exact value should not matter for the final result if it is an even power of 2.
  real :: tol_Hfrac ! Layers that together are smaller than this fraction of
                    ! the total water column can be merged for efficiency.
  real :: tol_solve ! The fractional tolerance with which to solve for the wave speeds [nondim]
  real :: tol_merge ! The fractional change in estimated wave speed that is allowed
                    ! when deciding to merge layers in the calculation [nondim]
  real :: rescale, I_rescale
  integer :: kf(szi_(g)) ! The number of active layers after filtering.
  integer, parameter :: max_itt = 10
  real :: lam_it(max_itt), det_it(max_itt), ddet_it(max_itt)
  logical :: use_EOS    ! If true, density is calculated from T & S using an equation of state.
  logical :: better_est ! If true, use an improved estimate of the first mode internal wave speed.
  logical :: merge      ! If true, merge the current layer with the one above.
  integer :: kc         ! The number of layers in the column after merging
  integer :: i, j, k, k2, itt, is, ie, js, je, nz
  real :: hw, sum_hc
  real :: gp      ! A limited local copy of gprime [L2 Z-1 T-2 ~> m s-2]
  real :: N2min   ! A minimum buoyancy frequency [T-2 ~> s-2]
  logical :: l_use_ebt_mode, calc_modal_structure
  real :: l_mono_N2_column_fraction, l_mono_N2_depth
  real :: mode_struct(szk_(g)), ms_min, ms_max, ms_sq

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

  if (.not. associated(cs)) call mom_error(fatal, "MOM_wave_speed: "// &
           "Module must be initialized before it is used.")
  if (present(full_halos)) then ; if (full_halos) then
    is = g%isd ; ie = g%ied ; js = g%jsd ; je = g%jed
  endif ; endif

  l2_to_z2 = us%L_to_Z**2

  l_use_ebt_mode = cs%use_ebt_mode
  if (present(use_ebt_mode)) l_use_ebt_mode = use_ebt_mode
  l_mono_n2_column_fraction = cs%mono_N2_column_fraction
  if (present(mono_n2_column_fraction)) l_mono_n2_column_fraction = mono_n2_column_fraction
  l_mono_n2_depth = cs%mono_N2_depth
  if (present(mono_n2_depth)) l_mono_n2_depth = mono_n2_depth
  calc_modal_structure = l_use_ebt_mode
  if (present(modal_structure)) calc_modal_structure = .true.
  if (calc_modal_structure) then
    do k=1,nz; do j=js,je; do i=is,ie
      modal_structure(i,j,k) = 0.0
    enddo ; enddo ; enddo
  endif

  s => tv%S ; t => tv%T
  g_rho0 = gv%g_Earth / gv%Rho0
  ! Simplifying the following could change answers at roundoff.
  z_to_pres = gv%Z_to_H * (gv%H_to_RZ * gv%g_Earth)
  use_eos = associated(tv%eqn_of_state)

  better_est = cs%better_cg1_est ; if (present(better_speed_est)) better_est = better_speed_est

  if (better_est) then
    tol_solve = cs%wave_speed_tol ; if (present(wave_speed_tol)) tol_solve = wave_speed_tol
    tol_hfrac  = 0.1*tol_solve ; tol_merge = tol_solve / real(nz)
  else
    tol_solve = 0.001 ; tol_hfrac  = 0.0001 ; tol_merge = 0.001
  endif

  ! The rescaling below can control the growth of the determinant provided that
  ! (tol_merge*cg1_min2/c2_scale > I_rescale).  For default values, this suggests a stable lower
  ! bound on min_speed of sqrt(nz/(tol_solve*rescale)) or 3e2/1024**2 = 2.9e-4 m/s for 90 layers.
  ! The upper bound on the rate of increase in the determinant is g'H/c2_scale < rescale or in the
  ! worst possible oceanic case of g'H < 0.5*10m/s2*1e4m = 5.e4 m2/s2 < 1024**2*c2_scale, suggesting
  ! that c2_scale can safely be set to 1/(16*1024**2), which would decrease the stable floor on
  ! min_speed to ~6.9e-8 m/s for 90 layers or 2.33e-7 m/s for 1000 layers.
  cg1_min2 = cs%min_speed2 ; if (present(min_speed)) cg1_min2 = min_speed**2
  rescale = 1024.0**4 ; i_rescale = 1.0/rescale
  c2_scale = us%m_s_to_L_T**2 / 4096.0**2 ! Other powers of 2 give identical results.

  min_h_frac = tol_hfrac / real(nz)
!$OMP parallel do default(none) shared(is,ie,js,je,nz,h,G,GV,US,min_h_frac,use_EOS,T,S,tv,&
!$OMP                                  calc_modal_structure,l_use_ebt_mode,modal_structure, &
!$OMP                                  l_mono_N2_column_fraction,l_mono_N2_depth,CS,   &
!$OMP                                  Z_to_pres,cg1,g_Rho0,rescale,I_rescale,L2_to_Z2, &
!$OMP                                  better_est,cg1_min2,tol_merge,tol_solve,c2_scale) &
!$OMP                          private(htot,hmin,kf,H_here,HxT_here,HxS_here,HxR_here, &
!$OMP                                  Hf,Tf,Sf,Rf,pres,T_int,S_int,drho_dT,drho_dS,   &
!$OMP                                  drxh_sum,kc,Hc,Hc_H,tC,sc,I_Hnew,gprime,&
!$OMP                                  Rc,speed2_tot,Igl,Igu,lam0,lam,lam_it,dlam, &
!$OMP                                  mode_struct,sum_hc,N2min,gp,hw,                 &
!$OMP                                  ms_min,ms_max,ms_sq,H_top,H_bot,I_Htot,merge,   &
!$OMP                                  det,ddet,detKm1,ddetKm1,detKm2,ddetKm2,det_it,ddet_it)
  do j=js,je
    !   First merge very thin layers with the one above (or below if they are
    ! at the top).  This also transposes the row order so that columns can
    ! be worked upon one at a time.
    do i=is,ie ; htot(i) = 0.0 ; enddo
    do k=1,nz ; do i=is,ie ; htot(i) = htot(i) + h(i,j,k)*gv%H_to_Z ; enddo ; enddo

    do i=is,ie
      hmin(i) = htot(i)*min_h_frac ; kf(i) = 1 ; h_here(i) = 0.0
      hxt_here(i) = 0.0 ; hxs_here(i) = 0.0 ; hxr_here(i) = 0.0
    enddo
    if (use_eos) then
      do k=1,nz ; do i=is,ie
        if ((h_here(i) > hmin(i)) .and. (h(i,j,k)*gv%H_to_Z > hmin(i))) then
          hf(kf(i),i) = h_here(i)
          tf(kf(i),i) = hxt_here(i) / h_here(i)
          sf(kf(i),i) = hxs_here(i) / h_here(i)
          kf(i) = kf(i) + 1

          ! Start a new layer
          h_here(i) = h(i,j,k)*gv%H_to_Z
          hxt_here(i) = (h(i,j,k)*gv%H_to_Z)*t(i,j,k)
          hxs_here(i) = (h(i,j,k)*gv%H_to_Z)*s(i,j,k)
        else
          h_here(i) = h_here(i) + h(i,j,k)*gv%H_to_Z
          hxt_here(i) = hxt_here(i) + (h(i,j,k)*gv%H_to_Z)*t(i,j,k)
          hxs_here(i) = hxs_here(i) + (h(i,j,k)*gv%H_to_Z)*s(i,j,k)
        endif
      enddo ; enddo
      do i=is,ie ; if (h_here(i) > 0.0) then
        hf(kf(i),i) = h_here(i)
        tf(kf(i),i) = hxt_here(i) / h_here(i)
        sf(kf(i),i) = hxs_here(i) / h_here(i)
      endif ; enddo
    else
      do k=1,nz ; do i=is,ie
        if ((h_here(i) > hmin(i)) .and. (h(i,j,k)*gv%H_to_Z > hmin(i))) then
          hf(kf(i),i) = h_here(i) ; rf(kf(i),i) = hxr_here(i) / h_here(i)
          kf(i) = kf(i) + 1

          ! Start a new layer
          h_here(i) = h(i,j,k)*gv%H_to_Z
          hxr_here(i) = (h(i,j,k)*gv%H_to_Z)*gv%Rlay(k)
        else
          h_here(i) = h_here(i) + h(i,j,k)*gv%H_to_Z
          hxr_here(i) = hxr_here(i) + (h(i,j,k)*gv%H_to_Z)*gv%Rlay(k)
        endif
      enddo ; enddo
      do i=is,ie ; if (h_here(i) > 0.0) then
        hf(kf(i),i) = h_here(i) ; rf(kf(i),i) = hxr_here(i) / h_here(i)
      endif ; enddo
    endif

    ! From this point, we can work on individual columns without causing memory to have page faults.
    do i=is,ie ; if (g%mask2dT(i,j) > 0.5) then
      if (use_eos) then
        pres(1) = 0.0 ; h_top(1) = 0.0
        do k=2,kf(i)
          pres(k) = pres(k-1) + z_to_pres*hf(k-1,i)
          t_int(k) = 0.5*(tf(k,i)+tf(k-1,i))
          s_int(k) = 0.5*(sf(k,i)+sf(k-1,i))
          h_top(k) = h_top(k-1) + hf(k-1,i)
        enddo
        call calculate_density_derivs(t_int, s_int, pres, drho_dt, drho_ds, &
                                      tv%eqn_of_state, (/2,kf(i)/) )

        ! Sum the reduced gravities to find out how small a density difference is negligibly small.
        drxh_sum = 0.0
        if (better_est) then
          ! This is an estimate that is correct for the non-EBT mode for 2 or 3 layers, or for
          ! clusters of massless layers at interfaces that can be grouped into 2 or 3 layers.
          ! For a uniform stratification and a huge number of layers uniformly distributed in
          ! density, this estimate is too large (as is desired) by a factor of pi^2/6 ~= 1.64.
          if (h_top(kf(i)) > 0.0) then
            i_htot = 1.0 / (h_top(kf(i)) + hf(kf(i),i))  ! = 1.0 / (H_top(K) + H_bot(K)) for all K.
            h_bot(kf(i)+1) = 0.0
            do k=kf(i),2,-1
              h_bot(k) = h_bot(k+1) + hf(k,i)
              drxh_sum = drxh_sum + ((h_top(k) * h_bot(k)) * i_htot) * &
                  max(0.0, drho_dt(k)*(tf(k,i)-tf(k-1,i)) + drho_ds(k)*(sf(k,i)-sf(k-1,i)))
            enddo
          endif
        else
          ! This estimate is problematic in that it goes like 1/nz for a large number of layers,
          ! but it is an overestimate (as desired) for a small number of layers, by at a factor
          ! of (H1+H2)**2/(H1*H2) >= 4 for two thick layers.
          do k=2,kf(i)
            drxh_sum = drxh_sum + 0.5*(hf(k-1,i)+hf(k,i)) * &
                max(0.0, drho_dt(k)*(tf(k,i)-tf(k-1,i)) + drho_ds(k)*(sf(k,i)-sf(k-1,i)))
          enddo
        endif
      else
        drxh_sum = 0.0
        if (better_est) then
          h_top(1) = 0.0
          do k=2,kf(i) ; h_top(k) = h_top(k-1) + hf(k-1,i) ; enddo
          if (h_top(kf(i)) > 0.0) then
            i_htot = 1.0 / (h_top(kf(i)) + hf(kf(i),i))  ! = 1.0 / (H_top(K) + H_bot(K)) for all K.
            h_bot(kf(i)+1) = 0.0
            do k=kf(i),2,-1
              h_bot(k) = h_bot(k+1) + hf(k,i)
              drxh_sum = drxh_sum + ((h_top(k) * h_bot(k)) * i_htot) * max(0.0,rf(k,i)-rf(k-1,i))
            enddo
          endif
        else
          do k=2,kf(i)
            drxh_sum = drxh_sum + 0.5*(hf(k-1,i)+hf(k,i)) * max(0.0,rf(k,i)-rf(k-1,i))
          enddo
        endif
      endif

      !   Find gprime across each internal interface, taking care of convective instabilities by
      ! merging layers.  If the estimated wave speed is too small, simply return zero.
      if (g_rho0 * drxh_sum <= cg1_min2) then
        cg1(i,j) = 0.0
        if (present(modal_structure)) modal_structure(i,j,:) = 0.
      else
        ! Merge layers to eliminate convective instabilities or exceedingly
        ! small reduced gravities.  Merging layers reduces the estimated wave speed by
        ! (rho(2)-rho(1))*h(1)*h(2) / H_tot.
        if (use_eos) then
          kc = 1
          hc(1) = hf(1,i) ; tc(1) = tf(1,i) ; sc(1) = sf(1,i)
          do k=2,kf(i)
            if (better_est) then
              merge = ((drho_dt(k)*(tf(k,i)-tc(kc)) + drho_ds(k)*(sf(k,i)-sc(kc))) * &
                       ((hc(kc) * hf(k,i))*i_htot) < 2.0 * tol_merge*drxh_sum)
            else
              merge = ((drho_dt(k)*(tf(k,i)-tc(kc)) + drho_ds(k)*(sf(k,i)-sc(kc))) * &
                       (hc(kc) + hf(k,i)) < 2.0 * tol_merge*drxh_sum)
            endif
            if (merge) then
              ! Merge this layer with the one above and backtrack.
              i_hnew = 1.0 / (hc(kc) + hf(k,i))
              tc(kc) = (hc(kc)*tc(kc) + hf(k,i)*tf(k,i)) * i_hnew
              sc(kc) = (hc(kc)*sc(kc) + hf(k,i)*sf(k,i)) * i_hnew
              hc(kc) = (hc(kc) + hf(k,i))
              ! Backtrack to remove any convective instabilities above...  Note
              ! that the tolerance is a factor of two larger, to avoid limit how
              ! far back we go.
              do k2=kc,2,-1
                if (better_est) then
                  merge = ((drho_dt(k2)*(tc(k2)-tc(k2-1)) + drho_ds(k2)*(sc(k2)-sc(k2-1))) * &
                           ((hc(k2) * hc(k2-1))*i_htot) < tol_merge*drxh_sum)
                else
                  merge = ((drho_dt(k2)*(tc(k2)-tc(k2-1)) + drho_ds(k2)*(sc(k2)-sc(k2-1))) * &
                           (hc(k2) + hc(k2-1)) < tol_merge*drxh_sum)
                endif
                if (merge) then
                  ! Merge the two bottommost layers.  At this point kc = k2.
                  i_hnew = 1.0 / (hc(kc) + hc(kc-1))
                  tc(kc-1) = (hc(kc)*tc(kc) + hc(kc-1)*tc(kc-1)) * i_hnew
                  sc(kc-1) = (hc(kc)*sc(kc) + hc(kc-1)*sc(kc-1)) * i_hnew
                  hc(kc-1) = (hc(kc) + hc(kc-1))
                  kc = kc - 1
                else ; exit ; endif
              enddo
            else
              ! Add a new layer to the column.
              kc = kc + 1
              drho_ds(kc) = drho_ds(k) ; drho_dt(kc) = drho_dt(k)
              tc(kc) = tf(k,i) ; sc(kc) = sf(k,i) ; hc(kc) = hf(k,i)
            endif
          enddo
          ! At this point there are kc layers and the gprimes should be positive.
          do k=2,kc ! Revisit this if non-Boussinesq.
            gprime(k) = g_rho0 * (drho_dt(k)*(tc(k)-tc(k-1)) + drho_ds(k)*(sc(k)-sc(k-1)))
          enddo
        else  ! .not.use_EOS
          ! Do the same with density directly...
          kc = 1
          hc(1) = hf(1,i) ; rc(1) = rf(1,i)
          do k=2,kf(i)
            if (better_est) then
              merge = ((rf(k,i) - rc(kc)) * ((hc(kc) * hf(k,i))*i_htot) < 2.0*tol_merge*drxh_sum)
            else
              merge = ((rf(k,i) - rc(kc)) * (hc(kc) + hf(k,i)) < 2.0*tol_merge*drxh_sum)
            endif
            if (merge) then
              ! Merge this layer with the one above and backtrack.
              rc(kc) = (hc(kc)*rc(kc) + hf(k,i)*rf(k,i)) / (hc(kc) + hf(k,i))
              hc(kc) = (hc(kc) + hf(k,i))
              ! Backtrack to remove any convective instabilities above...  Note
              ! that the tolerance is a factor of two larger, to avoid limit how
              ! far back we go.
              do k2=kc,2,-1
                if (better_est) then
                  merge = ((rc(k2)-rc(k2-1)) * ((hc(k2) * hc(k2-1))*i_htot) < tol_merge*drxh_sum)
                else
                  merge = ((rc(k2)-rc(k2-1)) * (hc(k2)+hc(k2-1)) < tol_merge*drxh_sum)
                endif
                if (merge) then
                  ! Merge the two bottommost layers.  At this point kc = k2.
                  rc(kc-1) = (hc(kc)*rc(kc) + hc(kc-1)*rc(kc-1)) / (hc(kc) + hc(kc-1))
                  hc(kc-1) = (hc(kc) + hc(kc-1))
                  kc = kc - 1
                else ; exit ; endif
              enddo
            else
              ! Add a new layer to the column.
              kc = kc + 1
              rc(kc) = rf(k,i) ; hc(kc) = hf(k,i)
            endif
          enddo
          ! At this point there are kc layers and the gprimes should be positive.
          do k=2,kc ! Revisit this if non-Boussinesq.
            gprime(k) = g_rho0 * (rc(k)-rc(k-1))
          enddo
        endif  ! use_EOS

        ! Sum the contributions from all of the interfaces to give an over-estimate
        ! of the first-mode wave speed.  Also populate Igl and Igu which are the
        ! non-leading diagonals of the tridiagonal matrix.
        if (kc >= 2) then
          speed2_tot = 0.0
          if (better_est) then
            h_top(1) = 0.0 ; h_bot(kc+1) = 0.0
            do k=2,kc+1 ; h_top(k) = h_top(k-1) + hc(k-1) ; enddo
            do k=kc,2,-1 ; h_bot(k) = h_bot(k+1) + hc(k) ; enddo
            i_htot = 0.0 ; if (h_top(kc+1) > 0.0) i_htot = 1.0 / h_top(kc+1)
          endif

          if (l_use_ebt_mode) then
            igu(1) = 0. ! Neumann condition for pressure modes
            sum_hc = hc(1)
            n2min = l2_to_z2*gprime(2)/hc(1)
            do k=2,kc
              hw = 0.5*(hc(k-1)+hc(k))
              gp = gprime(k)
              if (l_mono_n2_column_fraction>0. .or. l_mono_n2_depth>=0.) then
                if ( ((g%bathyT(i,j)-sum_hc < l_mono_n2_column_fraction*g%bathyT(i,j)) .or. &
                      ((l_mono_n2_depth >= 0.) .and. (sum_hc > l_mono_n2_depth))) .and. &
                     (l2_to_z2*gp > n2min*hw) ) then
                  ! Filters out regions where N2 increases with depth but only in a lower fraction
                  ! of the water column or below a certain depth.
                  gp = us%Z_to_L**2 * (n2min*hw)
                else
                  n2min = l2_to_z2 * gp/hw
                endif
              endif
              igu(k) = 1.0/(gp*hc(k))
              igl(k-1) = 1.0/(gp*hc(k-1))
              sum_hc = sum_hc + hc(k)
              if (better_est) then
                ! Estimate that the ebt_mode is sqrt(2) times the speed of the flat bottom modes.
                speed2_tot = speed2_tot + 2.0 * gprime(k)*((h_top(k) * h_bot(k)) * i_htot)
              else ! The ebt_mode wave should be faster than the flat-bottom mode, so 0.707 should be > 1?
                speed2_tot = speed2_tot + gprime(k)*(hc(k-1)+hc(k))*0.707
              endif
            enddo
           !Igl(kc) = 0. ! Neumann condition for pressure modes
            igl(kc) = 2.*igu(kc) ! Dirichlet condition for pressure modes
          else ! .not. l_use_ebt_mode
            do k=2,kc
              igl(k) = 1.0/(gprime(k)*hc(k)) ; igu(k) = 1.0/(gprime(k)*hc(k-1))
              if (better_est) then
                speed2_tot = speed2_tot + gprime(k)*((h_top(k) * h_bot(k)) * i_htot)
              else
                speed2_tot = speed2_tot + gprime(k)*(hc(k-1)+hc(k))
              endif
            enddo
          endif

          if (calc_modal_structure) then
            mode_struct(:) = 0.
            mode_struct(1:kc) = 1. ! Uniform flow, first guess
          endif

          ! Under estimate the first eigenvalue (overestimate the speed) to start with.
          if (calc_modal_structure) then
            lam0 = 0.5 / speed2_tot ; lam = lam0
          else
            lam0 = 1.0 / speed2_tot ; lam = lam0
          endif
          ! Find the determinant and its derivative with lam.
          do itt=1,max_itt
            lam_it(itt) = lam
            if (l_use_ebt_mode) then
              ! This initialization of det,ddet imply Neumann boundary conditions for horizontal
              ! velocity or pressure modes, so that first 3 rows of the matrix are
              !    /   b(1)-lam  igl(1)      0        0     0  ...  \
              !    |  igu(2)    b(2)-lam   igl(2)     0     0  ...  |
              !    |    0        igu(3)   b(3)-lam  igl(3)  0  ...  |
              ! The last two rows of the pressure equation matrix are
              !    |    ...  0  igu(kc-1)  b(kc-1)-lam  igl(kc-1)  |
              !    \    ...  0     0        igu(kc)     b(kc)-lam  /
              call tridiag_det(igu, igl, 1, kc, lam, det, ddet, row_scale=c2_scale)
            else
              ! This initialization of det,ddet imply Dirichlet boundary conditions for vertical
              ! velocity modes, so that first 3 rows of the matrix are
              !    /  b(2)-lam  igl(2)      0       0     0  ...  |
              !    |  igu(3)  b(3)-lam   igl(3)     0     0  ...  |
              !    |    0       igu(4)  b(4)-lam  igl(4)  0  ...  |
              ! The last three rows of the w equation matrix are
              !    |    ...   0  igu(kc-2)  b(kc-2)-lam  igl(kc-2)     0       |
              !    |    ...   0     0        igu(kc-1)  b(kc-1)-lam  igl(kc-1) |
              !    \    ...   0     0           0        igu(kc)    b(kc)-lam  /
              call tridiag_det(igu, igl, 2, kc, lam, det, ddet, row_scale=c2_scale)
            endif
            ! Use Newton's method iteration to find a new estimate of lam.
            det_it(itt) = det ; ddet_it(itt) = ddet

            if ((ddet >= 0.0) .or. (-det > -0.5*lam*ddet)) then
              ! lam was not an under-estimate, as intended, so Newton's method
              ! may not be reliable; lam must be reduced, but not by more
              ! than half.
              lam = 0.5 * lam
              dlam = -lam
            else  ! Newton's method is OK.
              dlam = - det / ddet
              lam = lam + dlam
            endif

            if (calc_modal_structure) then
              call tdma6(kc, igu, igl, lam, mode_struct)
              ms_min = mode_struct(1)
              ms_max = mode_struct(1)
              ms_sq = mode_struct(1)**2
              do k = 2,kc
                ms_min = min(ms_min, mode_struct(k))
                ms_max = max(ms_max, mode_struct(k))
                ms_sq = ms_sq + mode_struct(k)**2
              enddo
              if (ms_min<0. .and. ms_max>0.) then ! Any zero crossings => lam is too high
                lam = 0.5 * ( lam - dlam )
                dlam = -lam
                mode_struct(1:kc) = abs(mode_struct(1:kc)) / sqrt( ms_sq )
              else
                mode_struct(1:kc) = mode_struct(1:kc) / sqrt( ms_sq )
              endif
            endif

            if (abs(dlam) < tol_solve*lam) exit
          enddo

          cg1(i,j) = 0.0
          if (lam > 0.0) cg1(i,j) = 1.0 / sqrt(lam)

          if (present(modal_structure)) then
            if (mode_struct(1)/=0.) then ! Normalize
              mode_struct(1:kc) = mode_struct(1:kc) / mode_struct(1)
            else
              mode_struct(1:kc)=0.
            endif
            ! Note that remapping_core_h requires that the same units be used
            ! for both the source and target grid thicknesses, here [H ~> m or kg m-2].
            do k = 1,kc
              hc_h(k) = gv%Z_to_H * hc(k)
            enddo
            if (cs%remap_answers_2018) then
              call remapping_core_h(cs%remapping_CS, kc, hc_h(:), mode_struct, &
                                    nz, h(i,j,:), modal_structure(i,j,:), &
                                    1.0e-30*gv%m_to_H, 1.0e-10*gv%m_to_H)
            else
              call remapping_core_h(cs%remapping_CS, kc, hc_h(:), mode_struct, &
                                    nz, h(i,j,:), modal_structure(i,j,:), &
                                    gv%H_subroundoff, gv%H_subroundoff)
            endif
          endif
        else
          cg1(i,j) = 0.0
          if (present(modal_structure)) modal_structure(i,j,:) = 0.
        endif
      endif ! cg1 /= 0.0
    else
      cg1(i,j) = 0.0 ! This is a land point.
      if (present(modal_structure)) modal_structure(i,j,:) = 0.
    endif ; enddo ! i-loop
  enddo ! j-loop

end subroutine wave_speed

!> Solve a non-symmetric tridiagonal problem with the sum of the upper and lower diagnonals minus a
!! scalar contribution as the leading diagonal.
!! This uses the Thomas algorithm rather than the Hallberg algorithm since the matrix is not symmetric.
subroutine tdma6(n, a, c, lam, y)
  integer,            intent(in)    :: n !< Number of rows of matrix
  real, dimension(:), intent(in)    :: a !< Lower diagonal   [T2 L-2 ~> s2 m-2]
  real, dimension(:), intent(in)    :: c !< Upper diagonal   [T2 L-2 ~> s2 m-2]
  real,               intent(in)    :: lam !< Scalar subtracted from leading diagonal [T2 L-2 ~> s2 m-2]
  real, dimension(:), intent(inout) :: y !< RHS on entry, result on exit

  ! Local variables
  real :: lambda     ! A temporary variable in [T2 L-2 ~> s2 m-2]
  real :: beta(n)    ! A temporary variable in [T2 L-2 ~> s2 m-2]
  real :: I_beta(n)  ! A temporary variable in [L2 T-2 ~> m2 s-2]
  real :: yy(n)      ! A temporary variable with the same units as y on entry.
  integer :: k, m

  lambda = lam
  beta(1) = (a(1)+c(1)) - lambda
  if (beta(1)==0.) then ! lam was chosen too perfectly
    ! Change lambda and redo this first row
    lambda = (1. + 1.e-5) * lambda
    beta(1) = (a(1)+c(1)) - lambda
  endif
  i_beta(1) = 1. / beta(1)
  yy(1) = y(1)
  do k = 2, n
    beta(k) = ( (a(k)+c(k)) - lambda ) - a(k) * c(k-1) * i_beta(k-1)
    ! Perhaps the following 0 needs to become a tolerance to handle underflow?
    if (beta(k)==0.) then ! lam was chosen too perfectly
      ! Change lambda and redo everything up to row k
      lambda = (1. + 1.e-5) * lambda
      i_beta(1) = 1. / ( (a(1)+c(1)) - lambda )
      do m = 2, k
        i_beta(m) = 1. / ( ( (a(m)+c(m)) - lambda ) - a(m) * c(m-1) * i_beta(m-1) )
        yy(m) = y(m) + a(m) * yy(m-1) * i_beta(m-1)
      enddo
    else
      i_beta(k) = 1. / beta(k)
    endif
    yy(k) = y(k) + a(k) * yy(k-1) * i_beta(k-1)
  enddo
  ! The units of y change by a factor of [L2 T-2] in the following lines.
  y(n) = yy(n) * i_beta(n)
  do k = n-1, 1, -1
    y(k) = ( yy(k) + c(k) * y(k+1) ) * i_beta(k)
  enddo

end subroutine tdma6

!> Calculates the wave speeds for the first few barolinic modes.
subroutine wave_speeds(h, tv, G, GV, US, nmodes, cn, CS, full_halos, better_speed_est, &
                       min_speed, wave_speed_tol)
  type(ocean_grid_type),                    intent(in)  :: G !< Ocean grid structure
  type(verticalgrid_type),                  intent(in)  :: GV !< Vertical grid structure
  type(unit_scale_type),                    intent(in)  :: US !< A dimensional unit scaling type
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in)  :: h !< Layer thickness [H ~> m or kg m-2]
  type(thermo_var_ptrs),                    intent(in)  :: tv !< Thermodynamic variables
  integer,                                  intent(in)  :: nmodes !< Number of modes
  real, dimension(G%isd:G%ied,G%jsd:G%jed,nmodes), intent(out) :: cn !< Waves speeds [L T-1 ~> m s-1]
  type(wave_speed_cs), optional,            pointer     :: CS !< Control structure for MOM_wave_speed
  logical,             optional,            intent(in)  :: full_halos !< If true, do the calculation
                                                                      !! over the entire computational domain.
  logical, optional, intent(in) :: better_speed_est !< If true, use a more robust estimate of the first
                                     !! mode speed as the starting point for iterations.
  real,    optional, intent(in) :: min_speed !< If present, set a floor in the first mode speed
                                     !! below which 0 is returned [L T-1 ~> m s-1].
  real,    optional, intent(in) :: wave_speed_tol !< The fractional tolerance for finding the
                                     !! wave speeds [nondim]

  ! Local variables
  real, dimension(SZK_(G)+1) :: &
    dRho_dT, &    ! Partial derivative of density with temperature [R degC-1 ~> kg m-3 degC-1]
    dRho_dS, &    ! Partial derivative of density with salinity [R ppt-1 ~> kg m-3 ppt-1]
    pres, &       ! Interface pressure [R L2 T-2 ~> Pa]
    T_int, &      ! Temperature interpolated to interfaces [degC]
    S_int, &      ! Salinity interpolated to interfaces [ppt]
    H_top, &      ! The distance of each filtered interface from the ocean surface [Z ~> m]
    H_bot, &      ! The distance of each filtered interface from the bottom [Z ~> m]
    gprime        ! The reduced gravity across each interface [L2 Z-1 T-2 ~> m s-2].
  real, dimension(SZK_(G),SZI_(G)) :: &
    Hf, &         ! Layer thicknesses after very thin layers are combined [Z ~> m]
    Tf, &         ! Layer temperatures after very thin layers are combined [degC]
    Sf, &         ! Layer salinities after very thin layers are combined [ppt]
    Rf            ! Layer densities after very thin layers are combined [R ~> kg m-3]
  real, dimension(SZK_(G)) :: &
    Igl, Igu, &   ! The inverse of the reduced gravity across an interface times
                  ! the thickness of the layer below (Igl) or above (Igu) it, in [T2 L-2 ~> s2 m-2].
    hc, &         ! A column of layer thicknesses after convective istabilities are removed [Z ~> m]
    tc, &         ! A column of layer temperatures after convective istabilities are removed [degC]
    sc, &         ! A column of layer salinites after convective istabilities are removed [ppt]
    rc            ! A column of layer densities after convective istabilities are removed [R ~> kg m-3]
  real :: I_Htot  ! The inverse of the total filtered thicknesses [Z ~> m]
  real :: c1_thresh  ! if c1 is below this value, don't bother calculating
                     ! cn values for higher modes [L T-1 ~> m s-1]
  real :: c2_scale ! A scaling factor for wave speeds to help control the growth of the determinant
                   ! and its derivative with lam between rows of the Thomas algorithm solver.  The
                   ! exact value should not matter for the final result if it is an even power of 2.
  real :: det, ddet       ! determinant & its derivative of eigen system
  real :: lam_1           ! approximate mode-1 eigenvalue [T2 L-2 ~> s2 m-2]
  real :: lam_n           ! approximate mode-n eigenvalue [T2 L-2 ~> s2 m-2]
  real :: dlam            ! The change in estimates of the eigenvalue [T2 L-2 ~> s m-1]
  real :: lamMin          ! minimum lam value for root searching range [T2 L-2 ~> s2 m-2]
  real :: lamMax          ! maximum lam value for root searching range [T2 L-2 ~> s2 m-2]
  real :: lamInc          ! width of moving window for root searching [T2 L-2 ~> s2 m-2]
  real :: det_l,det_r     ! determinant value at left and right of window
  real :: ddet_l,ddet_r   ! derivative of determinant at left and right of window
  real :: det_sub,ddet_sub! derivative of determinant at subinterval endpoint
  real :: xl,xr           ! lam guesses at left and right of window [T2 L-2 ~> s2 m-2]
  real :: xl_sub          ! lam guess at left of subinterval window [T2 L-2 ~> s2 m-2]
  real,dimension(nmodes) :: &
          xbl,xbr         ! lam guesses bracketing a zero-crossing (root) [T2 L-2 ~> s2 m-2]
  integer :: numint       ! number of widows (intervals) in root searching range
  integer :: nrootsfound  ! number of extra roots found (not including 1st root)
  real :: Z_to_pres ! A conversion factor from thicknesses to pressure [R L2 T-2 Z-1 ~> Pa m-1]
  real, dimension(SZI_(G)) :: &
    htot, hmin, &  ! Thicknesses [Z ~> m]
    H_here, &      ! A thickness [Z ~> m]
    HxT_here, &    ! A layer integrated temperature [degC Z ~> degC m]
    HxS_here, &    ! A layer integrated salinity [ppt Z ~> ppt m]
    HxR_here       ! A layer integrated density [R Z ~> kg m-2]
  real :: speed2_tot ! overestimate of the mode-1 speed squared [L2 T-2 ~> m2 s-2]
  real :: speed2_min ! minimum mode speed (squared) to consider in root searching [L2 T-2 ~> m2 s-2]
  real :: cg1_min2 ! A floor in the squared first mode speed below which 0 is returned [L2 T-2 ~> m2 s-2]
  real, parameter :: reduct_factor = 0.5
                     ! A factor used in setting speed2_min [nondim]
  real :: I_Hnew   ! The inverse of a new layer thickness [Z-1 ~> m-1]
  real :: drxh_sum ! The sum of density differences across interfaces times thicknesses [R Z ~> kg m-2]
  real, pointer, dimension(:,:,:) :: T => null(), s => null()
  real :: g_Rho0   ! G_Earth/Rho0 [L2 T-2 Z-1 R-1 ~> m4 s-2 kg-1].
  real :: tol_Hfrac  ! Layers that together are smaller than this fraction of
                     ! the total water column can be merged for efficiency.
  real :: min_h_frac ! tol_Hfrac divided by the total number of layers [nondim].
  real :: tol_solve  ! The fractional tolerance with which to solve for the wave speeds [nondim].
  real :: tol_merge  ! The fractional change in estimated wave speed that is allowed
                     ! when deciding to merge layers in the calculation [nondim]
  integer :: kf(szi_(g)) ! The number of active layers after filtering.
  integer, parameter :: max_itt = 10
  logical :: use_EOS    ! If true, density is calculated from T & S using the equation of state.
  logical :: better_est ! If true, use an improved estimate of the first mode internal wave speed.
  logical :: merge      ! If true, merge the current layer with the one above.
  integer :: nsub       ! number of subintervals used for root finding
  integer, parameter :: sub_it_max = 4
                        ! maximum number of times to subdivide interval
                        ! for root finding (# intervals = 2**sub_it_max)
  logical :: sub_rootfound ! if true, subdivision has located root
  integer :: kc         ! The number of layers in the column after merging
  integer :: sub, sub_it
  integer :: i, j, k, k2, itt, is, ie, js, je, nz, row, iint, m, ig, jg

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

  if (present(cs)) then
    if (.not. associated(cs)) call mom_error(fatal, "MOM_wave_speed: "// &
           "Module must be initialized before it is used.")
  endif

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

  s => tv%S ; t => tv%T
  g_rho0 = gv%g_Earth / gv%Rho0
  ! Simplifying the following could change answers at roundoff.
  z_to_pres = gv%Z_to_H * (gv%H_to_RZ * gv%g_Earth)
  use_eos = associated(tv%eqn_of_state)
  c1_thresh = 0.01*us%m_s_to_L_T
  c2_scale = us%m_s_to_L_T**2 / 4096.0**2 ! Other powers of 2 give identical results.

  better_est = .false. ; if (present(cs)) better_est = cs%better_cg1_est
  if (present(better_speed_est)) better_est = better_speed_est
  if (better_est) then
    tol_solve = 0.001 ; if (present(cs)) tol_solve = cs%wave_speed_tol
    if (present(wave_speed_tol)) tol_solve = wave_speed_tol
    tol_hfrac  = 0.1*tol_solve ; tol_merge = tol_solve / real(nz)
  else
    tol_solve = 0.001 ; tol_hfrac  = 0.0001 ; tol_merge = 0.001
  endif
  cg1_min2 = 0.0 ; if (present(cs)) cg1_min2 = cs%min_speed2
  if (present(min_speed)) cg1_min2 = min_speed**2

  ! Zero out all wave speeds.  Values over land or for columns that are too weakly stratified
  ! are not changed from this zero value.
  cn(:,:,:) = 0.0

  min_h_frac = tol_hfrac / real(nz)
  !$OMP parallel do default(private) shared(is,ie,js,je,nz,h,G,GV,US,min_h_frac,use_EOS,T,S, &
  !$OMP                                     Z_to_pres,tv,cn,g_Rho0,nmodes,cg1_min2,better_est, &
  !$OMP                                     c1_thresh,tol_solve,tol_merge)
  do j=js,je
    !   First merge very thin layers with the one above (or below if they are
    ! at the top).  This also transposes the row order so that columns can
    ! be worked upon one at a time.
    do i=is,ie ; htot(i) = 0.0 ; enddo
    do k=1,nz ; do i=is,ie ; htot(i) = htot(i) + h(i,j,k)*gv%H_to_Z ; enddo ; enddo

    do i=is,ie
      hmin(i) = htot(i)*min_h_frac ; kf(i) = 1 ; h_here(i) = 0.0
      hxt_here(i) = 0.0 ; hxs_here(i) = 0.0 ; hxr_here(i) = 0.0
    enddo
    if (use_eos) then
      do k=1,nz ; do i=is,ie
        if ((h_here(i) > hmin(i)) .and. (h(i,j,k)*gv%H_to_Z > hmin(i))) then
          hf(kf(i),i) = h_here(i)
          tf(kf(i),i) = hxt_here(i) / h_here(i)
          sf(kf(i),i) = hxs_here(i) / h_here(i)
          kf(i) = kf(i) + 1

          ! Start a new layer
          h_here(i) = h(i,j,k)*gv%H_to_Z
          hxt_here(i) = (h(i,j,k)*gv%H_to_Z)*t(i,j,k)
          hxs_here(i) = (h(i,j,k)*gv%H_to_Z)*s(i,j,k)
        else
          h_here(i) = h_here(i) + h(i,j,k)*gv%H_to_Z
          hxt_here(i) = hxt_here(i) + (h(i,j,k)*gv%H_to_Z)*t(i,j,k)
          hxs_here(i) = hxs_here(i) + (h(i,j,k)*gv%H_to_Z)*s(i,j,k)
        endif
      enddo ; enddo
      do i=is,ie ; if (h_here(i) > 0.0) then
        hf(kf(i),i) = h_here(i)
        tf(kf(i),i) = hxt_here(i) / h_here(i)
        sf(kf(i),i) = hxs_here(i) / h_here(i)
      endif ; enddo
    else
      do k=1,nz ; do i=is,ie
        if ((h_here(i) > hmin(i)) .and. (h(i,j,k)*gv%H_to_Z > hmin(i))) then
          hf(kf(i),i) = h_here(i) ; rf(kf(i),i) = hxr_here(i) / h_here(i)
          kf(i) = kf(i) + 1

          ! Start a new layer
          h_here(i) = h(i,j,k)*gv%H_to_Z
          hxr_here(i) = (h(i,j,k)*gv%H_to_Z)*gv%Rlay(k)
        else
          h_here(i) = h_here(i) + h(i,j,k)*gv%H_to_Z
          hxr_here(i) = hxr_here(i) + (h(i,j,k)*gv%H_to_Z)*gv%Rlay(k)
        endif
      enddo ; enddo
      do i=is,ie ; if (h_here(i) > 0.0) then
        hf(kf(i),i) = h_here(i) ; rf(kf(i),i) = hxr_here(i) / h_here(i)
      endif ; enddo
    endif

    ! From this point, we can work on individual columns without causing memory to have page faults.
    do i=is,ie
      if (g%mask2dT(i,j) > 0.5) then
        if (use_eos) then
          pres(1) = 0.0 ; h_top(1) = 0.0
          do k=2,kf(i)
            pres(k) = pres(k-1) + z_to_pres*hf(k-1,i)
            t_int(k) = 0.5*(tf(k,i)+tf(k-1,i))
            s_int(k) = 0.5*(sf(k,i)+sf(k-1,i))
            h_top(k) = h_top(k-1) + hf(k-1,i)
          enddo
          call calculate_density_derivs(t_int, s_int, pres, drho_dt, drho_ds, &
                                        tv%eqn_of_state, (/2,kf(i)/) )

          ! Sum the reduced gravities to find out how small a density difference is negligibly small.
          drxh_sum = 0.0
          if (better_est) then
            ! This is an estimate that is correct for the non-EBT mode for 2 or 3 layers, or for
            ! clusters of massless layers at interfaces that can be grouped into 2 or 3 layers.
            ! For a uniform stratification and a huge number of layers uniformly distributed in
            ! density, this estimate is too large (as is desired) by a factor of pi^2/6 ~= 1.64.
            if (h_top(kf(i)) > 0.0) then
              i_htot = 1.0 / (h_top(kf(i)) + hf(kf(i),i))  ! = 1.0 / (H_top(K) + H_bot(K)) for all K.
              h_bot(kf(i)+1) = 0.0
              do k=kf(i),2,-1
                h_bot(k) = h_bot(k+1) + hf(k,i)
                drxh_sum = drxh_sum + ((h_top(k) * h_bot(k)) * i_htot) * &
                    max(0.0, drho_dt(k)*(tf(k,i)-tf(k-1,i)) + drho_ds(k)*(sf(k,i)-sf(k-1,i)))
              enddo
            endif
          else
            ! This estimate is problematic in that it goes like 1/nz for a large number of layers,
            ! but it is an overestimate (as desired) for a small number of layers, by at a factor
            ! of (H1+H2)**2/(H1*H2) >= 4 for two thick layers.
            do k=2,kf(i)
              drxh_sum = drxh_sum + 0.5*(hf(k-1,i)+hf(k,i)) * &
                  max(0.0, drho_dt(k)*(tf(k,i)-tf(k-1,i)) + drho_ds(k)*(sf(k,i)-sf(k-1,i)))
            enddo
          endif
        else
          drxh_sum = 0.0
          if (better_est) then
            h_top(1) = 0.0
            do k=2,kf(i) ; h_top(k) = h_top(k-1) + hf(k-1,i) ; enddo
              if (h_top(kf(i)) > 0.0) then
              i_htot = 1.0 / (h_top(kf(i)) + hf(kf(i),i))  ! = 1.0 / (H_top(K) + H_bot(K)) for all K.
              h_bot(kf(i)+1) = 0.0
              do k=kf(i),2,-1
                h_bot(k) = h_bot(k+1) + hf(k,i)
                drxh_sum = drxh_sum + ((h_top(k) * h_bot(k)) * i_htot) * max(0.0,rf(k,i)-rf(k-1,i))
              enddo
            endif
          else
            do k=2,kf(i)
              drxh_sum = drxh_sum + 0.5*(hf(k-1,i)+hf(k,i)) * max(0.0,rf(k,i)-rf(k-1,i))
            enddo
          endif
        endif

        !   Find gprime across each internal interface, taking care of convective
        ! instabilities by merging layers.
        if (g_rho0 * drxh_sum > cg1_min2) then
          ! Merge layers to eliminate convective instabilities or exceedingly
          ! small reduced gravities.  Merging layers reduces the estimated wave speed by
          ! (rho(2)-rho(1))*h(1)*h(2) / H_tot.
          if (use_eos) then
            kc = 1
            hc(1) = hf(1,i) ; tc(1) = tf(1,i) ; sc(1) = sf(1,i)
            do k=2,kf(i)
              if (better_est) then
                merge = ((drho_dt(k)*(tf(k,i)-tc(kc)) + drho_ds(k)*(sf(k,i)-sc(kc))) * &
                         ((hc(kc) * hf(k,i))*i_htot) < 2.0 * tol_merge*drxh_sum)
              else
                merge = ((drho_dt(k)*(tf(k,i)-tc(kc)) + drho_ds(k)*(sf(k,i)-sc(kc))) * &
                         (hc(kc) + hf(k,i)) < 2.0 * tol_merge*drxh_sum)
              endif
              if (merge) then
                ! Merge this layer with the one above and backtrack.
                i_hnew = 1.0 / (hc(kc) + hf(k,i))
                tc(kc) = (hc(kc)*tc(kc) + hf(k,i)*tf(k,i)) * i_hnew
                sc(kc) = (hc(kc)*sc(kc) + hf(k,i)*sf(k,i)) * i_hnew
                hc(kc) = (hc(kc) + hf(k,i))
                ! Backtrack to remove any convective instabilities above...  Note
                ! that the tolerance is a factor of two larger, to avoid limit how
                ! far back we go.
                do k2=kc,2,-1
                  if (better_est) then
                    merge = ((drho_dt(k2)*(tc(k2)-tc(k2-1)) + drho_ds(k2)*(sc(k2)-sc(k2-1))) * &
                             ((hc(k2) * hc(k2-1))*i_htot) < tol_merge*drxh_sum)
                  else
                    merge = ((drho_dt(k2)*(tc(k2)-tc(k2-1)) + drho_ds(k2)*(sc(k2)-sc(k2-1))) * &
                             (hc(k2) + hc(k2-1)) < tol_merge*drxh_sum)
                  endif
                  if (merge) then
                    ! Merge the two bottommost layers.  At this point kc = k2.
                    i_hnew = 1.0 / (hc(kc) + hc(kc-1))
                    tc(kc-1) = (hc(kc)*tc(kc) + hc(kc-1)*tc(kc-1)) * i_hnew
                    sc(kc-1) = (hc(kc)*sc(kc) + hc(kc-1)*sc(kc-1)) * i_hnew
                    hc(kc-1) = (hc(kc) + hc(kc-1))
                    kc = kc - 1
                  else ; exit ; endif
                enddo
              else
                ! Add a new layer to the column.
                kc = kc + 1
                drho_ds(kc) = drho_ds(k) ; drho_dt(kc) = drho_dt(k)
                tc(kc) = tf(k,i) ; sc(kc) = sf(k,i) ; hc(kc) = hf(k,i)
              endif
            enddo
            ! At this point there are kc layers and the gprimes should be positive.
            do k=2,kc ! Revisit this if non-Boussinesq.
              gprime(k) = g_rho0 * (drho_dt(k)*(tc(k)-tc(k-1)) + drho_ds(k)*(sc(k)-sc(k-1)))
            enddo
          else  ! .not.use_EOS
            ! Do the same with density directly...
            kc = 1
            hc(1) = hf(1,i) ; rc(1) = rf(1,i)
            do k=2,kf(i)
              if (better_est) then
                merge = ((rf(k,i) - rc(kc)) * ((hc(kc) * hf(k,i))*i_htot) < 2.0 * tol_merge*drxh_sum)
              else
                merge = ((rf(k,i) - rc(kc)) * (hc(kc) + hf(k,i)) < 2.0*tol_merge*drxh_sum)
              endif
              if (merge) then
                ! Merge this layer with the one above and backtrack.
                rc(kc) = (hc(kc)*rc(kc) + hf(k,i)*rf(k,i)) / (hc(kc) + hf(k,i))
                hc(kc) = (hc(kc) + hf(k,i))
                ! Backtrack to remove any convective instabilities above...  Note
                ! that the tolerance is a factor of two larger, to avoid limit how
                ! far back we go.
                do k2=kc,2,-1
                  if (better_est) then
                    merge = ((rc(k2)-rc(k2-1)) * ((hc(k2) * hc(k2-1))*i_htot) < tol_merge*drxh_sum)
                  else
                    merge = ((rc(k2)-rc(k2-1)) * (hc(k2)+hc(k2-1)) < tol_merge*drxh_sum)
                  endif
                  if (merge) then
                    ! Merge the two bottommost layers.  At this point kc = k2.
                    rc(kc-1) = (hc(kc)*rc(kc) + hc(kc-1)*rc(kc-1)) / (hc(kc) + hc(kc-1))
                    hc(kc-1) = (hc(kc) + hc(kc-1))
                    kc = kc - 1
                  else ; exit ; endif
                enddo
              else
                ! Add a new layer to the column.
                kc = kc + 1
                rc(kc) = rf(k,i) ; hc(kc) = hf(k,i)
              endif
            enddo
            ! At this point there are kc layers and the gprimes should be positive.
            do k=2,kc ! Revisit this if non-Boussinesq.
              gprime(k) = g_rho0 * (rc(k)-rc(k-1))
            enddo
          endif  ! use_EOS

          !-----------------NOW FIND WAVE SPEEDS---------------------------------------
          ! ig = i + G%idg_offset ; jg = j + G%jdg_offset
          !   Sum the contributions from all of the interfaces to give an over-estimate
          ! of the first-mode wave speed.  Also populate Igl and Igu which are the
          ! non-leading diagonals of the tridiagonal matrix.
          if (kc >= 2) then
            ! initialize speed2_tot
            speed2_tot = 0.0
            if (better_est) then
              h_top(1) = 0.0 ; h_bot(kc+1) = 0.0
              do k=2,kc+1 ; h_top(k) = h_top(k-1) + hc(k-1) ; enddo
              do k=kc,2,-1 ; h_bot(k) = h_bot(k+1) + hc(k) ; enddo
              i_htot = 0.0 ; if (h_top(kc+1) > 0.0) i_htot = 1.0 / h_top(kc+1)
            endif

            ! Calculate Igu, Igl, depth, and N2 at each interior interface
            ! [excludes surface (K=1) and bottom (K=kc+1)]
            do k=2,kc
              igl(k) = 1.0/(gprime(k)*hc(k)) ; igu(k) = 1.0/(gprime(k)*hc(k-1))
              if (better_est) then
                speed2_tot = speed2_tot + gprime(k)*((h_top(k) * h_bot(k)) * i_htot)
              else
                speed2_tot = speed2_tot + gprime(k)*(hc(k-1)+hc(k))
              endif
            enddo

            ! Under estimate the first eigenvalue (overestimate the speed) to start with.
            lam_1 = 1.0 / speed2_tot

            ! Find the first eigen value
            do itt=1,max_itt
              ! calculate the determinant of (A-lam_1*I)
              call tridiag_det(igu, igl, 2, kc, lam_1, det, ddet, row_scale=c2_scale)

              ! If possible, use Newton's method iteration to find a new estimate of lam_1
              !det = det_it(itt) ; ddet = ddet_it(itt)
              if ((ddet >= 0.0) .or. (-det > -0.5*lam_1*ddet)) then
                ! lam_1 was not an under-estimate, as intended, so Newton's method
                ! may not be reliable; lam_1 must be reduced, but not by more than half.
                lam_1 = 0.5 * lam_1
                dlam = -lam_1
              else  ! Newton's method is OK.
                dlam = - det / ddet
                lam_1 = lam_1 + dlam
              endif

              if (abs(dlam) < tol_solve*lam_1) exit
            enddo

            if (lam_1 > 0.0) cn(i,j,1) = 1.0 / sqrt(lam_1)

            ! Find other eigen values if c1 is of significant magnitude, > cn_thresh
            nrootsfound = 0    ! number of extra roots found (not including 1st root)
            if (nmodes>1 .and. kc>=nmodes+1 .and. cn(i,j,1)>c1_thresh) then
              ! Set the the range to look for the other desired eigen values
              ! set min value just greater than the 1st root (found above)
              lammin = lam_1*(1.0 + tol_solve)
              ! set max value based on a low guess at wavespeed for highest mode
              speed2_min = (reduct_factor*cn(i,j,1)/real(nmodes))**2
              lammax = 1.0 / speed2_min
              ! set width of interval (not sure about this - BDM)
              laminc = 0.5*lam_1
              ! set number of intervals within search range
              numint = nint((lammax - lammin)/laminc)

              !   Find intervals containing zero-crossings (roots) of the determinant
              ! that are beyond the first root

              ! find det_l of first interval (det at left endpoint)
              call tridiag_det(igu, igl, 2, kc, lammin, det_l, ddet_l, row_scale=c2_scale)
              ! move interval window looking for zero-crossings************************
              do iint=1,numint
                xr = lammin + laminc * iint
                xl = xr - laminc
                call tridiag_det(igu, igl, 2, kc, xr, det_r, ddet_r, row_scale=c2_scale)
                if (det_l*det_r < 0.0) then  ! if function changes sign
                  if (det_l*ddet_l < 0.0) then ! if function at left is headed to zero
                    nrootsfound = nrootsfound + 1
                    xbl(nrootsfound) = xl
                    xbr(nrootsfound) = xr
                  else
                    !   function changes sign but has a local max/min in interval,
                    ! try subdividing interval as many times as necessary (or sub_it_max).
                    ! loop that increases number of subintervals:
                    !call MOM_error(WARNING, "determinant changes sign"// &
                    !            "but has a local max/min in interval;"//&
                    !            " reduce increment in lam.")
                    ! begin subdivision loop -------------------------------------------
                    sub_rootfound = .false. ! initialize
                    do sub_it=1,sub_it_max
                      nsub = 2**sub_it ! number of subintervals; nsub=2,4,8,...
                      ! loop over each subinterval:
                      do sub=1,nsub-1,2 ! only check odds; sub = 1; 1,3; 1,3,5,7;...
                        xl_sub = xl + laminc/(nsub)*sub
                        call tridiag_det(igu, igl, 2, kc, xl_sub, det_sub, ddet_sub, &
                                         row_scale=c2_scale)
                        if (det_sub*det_r < 0.0) then  ! if function changes sign
                          if (det_sub*ddet_sub < 0.0) then ! if function at left is headed to zero
                            sub_rootfound = .true.
                            nrootsfound = nrootsfound + 1
                            xbl(nrootsfound) = xl_sub
                            xbr(nrootsfound) = xr
                            exit ! exit sub loop
                          endif ! headed toward zero
                        endif ! sign change
                      enddo ! sub-loop
                      if (sub_rootfound) exit ! root has been found, exit sub_it loop
                      !   Otherwise, function changes sign but has a local max/min in one of the
                      ! sub intervals, try subdividing again unless sub_it_max has been reached.
                      if (sub_it == sub_it_max) then
                        call mom_error(warning, "wave_speed: root not found "// &
                                       " after sub_it_max subdivisions of original"// &
                                       " interval.")
                      endif ! sub_it == sub_it_max
                    enddo ! sub_it-loop-------------------------------------------------
                  endif ! det_l*ddet_l < 0.0
                endif ! det_l*det_r < 0.0
                ! exit iint-loop if all desired roots have been found
                if (nrootsfound >= nmodes-1) then
                  ! exit if all additional roots found
                  exit
                elseif (iint == numint) then
                  ! oops, lamMax not large enough - could add code to increase (BDM)
                  ! set unfound modes to zero for now (BDM)
                  !   cn(i,j,nrootsfound+2:nmodes) = 0.0
                else
                  ! else shift interval and keep looking until nmodes or numint is reached
                  det_l = det_r
                  ddet_l = ddet_r
                endif
              enddo ! iint-loop

              ! Use Newton's method to find the roots within the identified windows
              do m=1,nrootsfound ! loop over the root-containing widows (excluding 1st mode)
                lam_n = xbl(m) ! first guess is left edge of window
                do itt=1,max_itt
                  ! calculate the determinant of (A-lam_n*I)
                  call tridiag_det(igu, igl, 2, kc, lam_n, det, ddet, row_scale=c2_scale)
                  ! Use Newton's method to find a new estimate of lam_n
                  dlam = - det / ddet
                  lam_n = lam_n + dlam
                  if (abs(dlam) < tol_solve*lam_1)  exit
                enddo ! itt-loop
                ! calculate nth mode speed
                if (lam_n > 0.0) cn(i,j,m+1) = 1.0 / sqrt(lam_n)
              enddo ! n-loop
            endif ! if nmodes>1 .and. kc>nmodes .and. c1>c1_thresh
          endif ! if more than 2 layers
        endif ! if drxh_sum < 0
      endif ! if not land
    enddo ! i-loop
  enddo ! j-loop

end subroutine wave_speeds

!> Calculate the determinant of a tridiagonal matrix with diagonals a,b-lam,c and its derivative
!! with lam, where lam is constant across rows.  Only the ratio of det to its derivative and their
!! signs are typically used, so internal rescaling by consistent factors are used to avoid
!! over- or underflow.
subroutine tridiag_det(a, c, ks, ke, lam, det, ddet, row_scale)
  real, dimension(:), intent(in) :: a     !< Lower diagonal of matrix (first entry unused)
  real, dimension(:), intent(in) :: c     !< Upper diagonal of matrix (last entry unused)
  integer,            intent(in) :: ks    !< Starting index to use in determinant
  integer,            intent(in) :: ke    !< Ending index to use in determinant
  real,               intent(in) :: lam   !< Value subtracted from b
  real,               intent(out):: det   !< Determinant
  real,               intent(out):: ddet  !< Derivative of determinant with lam
  real,     optional, intent(in) :: row_scale !< A scaling factor of the rows of the
                                      !! matrix to limit the growth of the determinant
  ! Local variables
  real :: detKm1, detKm2   ! Cumulative value of the determinant for the previous two layers.
  real :: ddetKm1, ddetKm2 ! Derivative of the cumulative determinant with lam for the previous two layers.
  real, parameter :: rescale = 1024.0**4 ! max value of determinant allowed before rescaling
  real :: rscl      ! A rescaling factor that is applied succesively to each row.
  real :: I_rescale ! inverse of rescale
  integer :: k      ! row (layer interface) index

  i_rescale = 1.0 / rescale
  rscl = 1.0 ; if (present(row_scale)) rscl = row_scale

  detkm1 = 1.0 ; ddetkm1 = 0.0
  det = (a(ks)+c(ks)) - lam ; ddet = -1.0
  do k=ks+1,ke
    ! Shift variables and rescale rows to avoid over- or underflow.
    detkm2 = row_scale*detkm1 ; ddetkm2 = row_scale*ddetkm1
    detkm1 = row_scale*det    ; ddetkm1 = row_scale*ddet

    det =  ((a(k)+c(k))-lam)*detkm1  - (a(k)*c(k-1))*detkm2
    ddet = ((a(k)+c(k))-lam)*ddetkm1 - (a(k)*c(k-1))*ddetkm2 - detkm1

    ! Rescale det & ddet if det is getting too large or too small.
    if (abs(det) > rescale) then
      det = i_rescale*det ; detkm1 = i_rescale*detkm1
      ddet = i_rescale*ddet ; ddetkm1 = i_rescale*ddetkm1
    elseif (abs(det) < i_rescale) then
      det = rescale*det ; detkm1 = rescale*detkm1
      ddet = rescale*ddet ; ddetkm1 = rescale*ddetkm1
    endif
  enddo

end subroutine tridiag_det

!> Initialize control structure for MOM_wave_speed
subroutine wave_speed_init(CS, use_ebt_mode, mono_N2_column_fraction, mono_N2_depth, remap_answers_2018, &
                           better_speed_est, min_speed, wave_speed_tol)
  type(wave_speed_cs), pointer :: CS !< Control structure for MOM_wave_speed
  logical, optional, intent(in) :: use_ebt_mode  !< If true, use the equivalent
                                     !! barotropic mode instead of the first baroclinic mode.
  real,    optional, intent(in) :: mono_N2_column_fraction !< The lower fraction of water column over
                                     !! which N2 is limited as monotonic for the purposes of
                                     !! calculating the vertical modal structure.
  real,    optional, intent(in) :: mono_N2_depth !< The depth below which N2 is limited
                                     !! as monotonic for the purposes of calculating the
                                     !! vertical modal structure [Z ~> m].
  logical, optional, intent(in) :: remap_answers_2018 !< If true, use the order of arithmetic and expressions
                                     !! that recover the remapping answers from 2018.  Otherwise
                                     !! use more robust but mathematically equivalent expressions.
  logical, optional, intent(in) :: better_speed_est !< If true, use a more robust estimate of the first
                                     !! mode speed as the starting point for iterations.
  real,    optional, intent(in) :: min_speed !< If present, set a floor in the first mode speed
                                     !! below which 0 is returned [L T-1 ~> m s-1].
  real,    optional, intent(in) :: wave_speed_tol !< The fractional tolerance for finding the
                                     !! wave speeds [nondim]

  ! This include declares and sets the variable "version".
# include "version_variable.h"
  character(len=40)  :: mdl = "MOM_wave_speed"  ! This module's name.

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

  ! Write all relevant parameters to the model log.
  call log_version(mdl, version)

  call wave_speed_set_param(cs, use_ebt_mode=use_ebt_mode, mono_n2_column_fraction=mono_n2_column_fraction, &
                            better_speed_est=better_speed_est, min_speed=min_speed, wave_speed_tol=wave_speed_tol)

  call initialize_remapping(cs%remapping_CS, 'PLM', boundary_extrapolation=.false., &
                            answers_2018=cs%remap_answers_2018)

end subroutine wave_speed_init

!> Sets internal parameters for MOM_wave_speed
subroutine wave_speed_set_param(CS, use_ebt_mode, mono_N2_column_fraction, mono_N2_depth, remap_answers_2018, &
                                better_speed_est, min_speed, wave_speed_tol)
  type(wave_speed_cs), pointer  :: CS !< Control structure for MOM_wave_speed
  logical, optional, intent(in) :: use_ebt_mode  !< If true, use the equivalent
                                      !! barotropic mode instead of the first baroclinic mode.
  real,    optional, intent(in) :: mono_N2_column_fraction !< The lower fraction of water column over
                                      !! which N2 is limited as monotonic for the purposes of
                                      !! calculating the vertical modal structure.
  real,    optional, intent(in) :: mono_N2_depth !< The depth below which N2 is limited
                                      !! as monotonic for the purposes of calculating the
                                      !! vertical modal structure [Z ~> m].
  logical, optional, intent(in) :: remap_answers_2018 !< If true, use the order of arithmetic and expressions
                                      !! that recover the remapping answers from 2018.  Otherwise
                                      !! use more robust but mathematically equivalent expressions.
  logical, optional, intent(in) :: better_speed_est !< If true, use a more robust estimate of the first
                                     !! mode speed as the starting point for iterations.
  real,    optional, intent(in) :: min_speed !< If present, set a floor in the first mode speed
                                     !! below which 0 is returned [L T-1 ~> m s-1].
  real,    optional, intent(in) :: wave_speed_tol !< The fractional tolerance for finding the
                                     !! wave speeds [nondim]

  if (.not.associated(cs)) call mom_error(fatal, &
     "wave_speed_set_param called with an associated control structure.")

  if (present(use_ebt_mode)) cs%use_ebt_mode = use_ebt_mode
  if (present(mono_n2_column_fraction)) cs%mono_N2_column_fraction = mono_n2_column_fraction
  if (present(mono_n2_depth)) cs%mono_N2_depth = mono_n2_depth
  if (present(remap_answers_2018)) cs%remap_answers_2018 = remap_answers_2018
  if (present(better_speed_est)) cs%better_cg1_est = better_speed_est
  if (present(min_speed)) cs%min_speed2 = min_speed**2
  if (present(wave_speed_tol)) cs%wave_speed_tol = wave_speed_tol

end subroutine wave_speed_set_param

!> \namespace mom_wave_speed

!!
!! Subroutine wave_speed() solves for the first baroclinic mode wave speed.  (It could
!! solve for all the wave speeds, but the iterative approach taken here means
!! that this is not particularly efficient.)
!!
!! If `e(k)` is the perturbation interface height, this means solving for the
!! smallest eigenvalue (`lam` = 1/c^2) of the system
!!
!! \verbatim
!!   -Igu(k)*e(k-1) + (Igu(k)+Igl(k)-lam)*e(k) - Igl(k)*e(k+1) = 0.0
!! \endverbatim
!!
!! with rigid lid boundary conditions e(1) = e(nz+1) = 0.0 giving
!!
!! \verbatim
!!   (Igu(2)+Igl(2)-lam)*e(2) - Igl(2)*e(3) = 0.0
!!   -Igu(nz)*e(nz-1) + (Igu(nz)+Igl(nz)-lam)*e(nz) = 0.0
!! \endverbatim
!!
!! Here
!! \verbatim
!!   Igl(k) = 1.0/(gprime(K)*h(k)) ; Igu(k) = 1.0/(gprime(K)*h(k-1))
!! \endverbatim
!!
!! Alternately, these same eigenvalues can be found from the second smallest
!! eigenvalue of the Montgomery potential (M(k)) calculation:
!!
!! \verbatim
!!   -Igl(k)*M(k-1) + (Igl(k)+Igu(k+1)-lam)*M(k) - Igu(k+1)*M(k+1) = 0.0
!! \endverbatim
!!
!! with rigid lid and flat bottom boundary conditions
!!
!! \verbatim
!!   (Igu(2)-lam)*M(1) - Igu(2)*M(2) = 0.0
!!   -Igl(nz)*M(nz-1) + (Igl(nz)-lam)*M(nz) = 0.0
!! \endverbatim
!!
!! Note that the barotropic mode has been eliminated from the rigid lid
!! interface height equations, hence the matrix is one row smaller.  Without
!! the rigid lid, the top boundary condition is simpler to implement with
!! the M equations.

end module mom_wave_speed