MOM_wave_speed.F90

Go to the documentation of this file.

!> 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.
  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.
  type(remapping_cs) :: remapping_cs   !< Used for vertical remapping when calculating equivalent barotropic
                                       !! mode structure.
  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)
  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]
 
  ! 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 [Pa]
    t_int, &      ! Temperature interpolated to interfaces [degC]
    s_int, &      ! Salinity interpolated to interfaces [ppt]
    gprime        ! The reduced gravity across each interface [L2 Z-1 T-2 ~> m s-2].
  real, dimension(SZK_(G)) :: &
    igl, igu, igd ! The inverse of the reduced gravity across an interface times
                  ! the thickness of the layer below (Igl) or above (Igu) it.
                  ! Their sum, Igd, is provided for the tridiagonal solver.  [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 :: 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_pa  ! A conversion factor from thicknesses (in Z) to pressure (in Pa)
  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 :: i_hnew   ! The inverse of a new layer thickness [Z-1 ~> m-1]
  real :: drxh_sum ! The sum of density diffrences 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, parameter :: tol1  = 0.0001, tol2 = 0.001
  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 :: rescale, i_rescale
  integer :: kf(szi_(g))
  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.
  integer :: kc
  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
  z_to_pa = gv%Z_to_H * gv%H_to_Pa
  use_eos = associated(tv%eqn_of_state)
 
  rescale = 1024.0**4 ; i_rescale = 1.0/rescale
  ! The following two lines give identical results:
  ! c2_scale = 16.0 * US%m_s_to_L_T**2
  c2_scale = us%m_s_to_L_T**2
 
  min_h_frac = tol1 / 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_Pa,cg1,g_Rho0,rescale,I_rescale,L2_to_Z2,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,           &
!$OMP                                  drho_dS,drxh_sum,kc,Hc,Hc_H,Tc,Sc,I_Hnew,gprime,&
!$OMP                                  Rc,speed2_tot,Igl,Igu,Igd,lam0,lam,lam_it,dlam, &
!$OMP                                  mode_struct,sum_hc,N2min,gp,hw,                 &
!$OMP                                  ms_min,ms_max,ms_sq,                            &
!$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
        do k=2,kf(i)
          pres(k) = pres(k-1) + z_to_pa*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))
        enddo
        call calculate_density_derivs(t_int, s_int, pres, drho_dt, drho_ds, 2, &
                                      kf(i)-1, tv%eqn_of_state, scale=us%kg_m3_to_R)
 
        ! Sum the reduced gravities to find out how small a density difference
        ! is negligibly small.
        drxh_sum = 0.0
        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
      else
        drxh_sum = 0.0
        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
 
      if (calc_modal_structure) then
        mode_struct(:) = 0.
      endif
 
  !   Find gprime across each internal interface, taking care of convective
  ! instabilities by merging layers.
      if (drxh_sum <= 0.0) then
        cg1(i,j) = 0.0
      else
        ! Merge layers to eliminate convective instabilities or exceedingly
        ! small reduced gravities.
        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 ((drho_dt(k)*(tf(k,i)-tc(kc)) + drho_ds(k)*(sf(k,i)-sc(kc))) * &
                (hc(kc) + hf(k,i)) < 2.0 * tol2*drxh_sum) 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 ((drho_dt(k2)*(tc(k2)-tc(k2-1)) + drho_ds(k2)*(sc(k2)-sc(k2-1))) * &
                    (hc(k2) + hc(k2-1)) < tol2*drxh_sum) 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 ((rf(k,i) - rc(kc)) * (hc(kc) + hf(k,i)) < 2.0*tol2*drxh_sum) 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 ((rc(k2)-rc(k2-1)) * (hc(k2)+hc(k2-1)) < tol2*drxh_sum) 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 (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))
              speed2_tot = speed2_tot + gprime(k)*(hc(k-1)+hc(k))*0.707
              sum_hc = sum_hc + hc(k)
            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))
              speed2_tot = speed2_tot + gprime(k)*(hc(k-1)+hc(k))
            enddo
          endif
 
          if (calc_modal_structure) then
            mode_struct(1:kc) = 1. ! Uniform flow, first guess
          endif
 
          ! 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 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  ...  |
              ! which is consistent if the eigenvalue problem is for horizontal velocity or pressure modes.
             !detKm1 = c2_scale*(Igl(1)-lam) ; ddetKm1 = -1.0*c2_scale
             !det = (Igu(2)+Igl(2)-lam)*detKm1 - (Igu(2)*Igl(1)) ; ddet = (Igu(2)+Igl(2)-lam)*ddetKm1 - detKm1
              detkm1 = 1.0 ; ddetkm1 = 0.0
              det = (igl(1)-lam) ; ddet = -1.0
              if (kc>1) then
                ! Shift variables and rescale rows to avoid over- or underflow.
                detkm2 = c2_scale*detkm1 ; ddetkm2 = c2_scale*ddetkm1
                detkm1 = c2_scale*det    ; ddetkm1 = c2_scale*ddet
                det = (igu(2)+igl(2)-lam)*detkm1 - (igu(2)*igl(1))*detkm2
                ddet = (igu(2)+igl(2)-lam)*ddetkm1 - (igu(2)*igl(1))*ddetkm2 - detkm1
              endif
              ! 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  /
            else
              ! This initialization of det,ddet imply Dirichlet boundary conditions 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       igu43)  b(4)-lam  igl(4)  0  ...  |
              ! which is consistent if the eigenvalue problem is for vertical velocity modes.
              detkm1 = 1.0 ; ddetkm1 = 0.0
              det = (igu(2) + igl(2) - lam) ; ddet = -1.0
              ! The last three rows of the w equation matrix are
              !    |    ...   0  igu(kc-1)  b(kc-1)-lam  igl(kc-1)     0       |
              !    |    ...   0     0        igu(kc-1)  b(kc-1)-lam  igl(kc-1) |
              !    \    ...   0     0           0        igu(kc)    b(kc)-lam  /
            endif
            do k=3,kc
              ! Shift variables and rescale rows to avoid over- or underflow.
              detkm2 = c2_scale*detkm1 ; ddetkm2 = c2_scale*ddetkm1
              detkm1 = c2_scale*det    ; ddetkm1 = c2_scale*ddet
 
              det = (igu(k)+igl(k)-lam)*detkm1 - (igu(k)*igl(k-1))*detkm2
              ddet = (igu(k)+igl(k)-lam)*ddetkm1 - (igu(k)*igl(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
            ! 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
              do k = 1,kc
                igd(k) = igu(k) + igl(k)
              enddo
              call tdma6(kc, -igu, igd, -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) < tol2*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
            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)
          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 a scalar contribution to the leading diagonal.
!! This uses the Thomas algorithm rather than the Hallberg algorithm since the matrix is not symmetric.
subroutine tdma6(n, a, b, c, lam, y)
  integer,            intent(in)    :: n !< Number of rows of matrix
  real, dimension(n), intent(in)    :: a !< Lower diagonal   [T2 L-2 ~> s2 m-2]
  real, dimension(n), intent(in)    :: b !< Leading diagonal [T2 L-2 ~> s2 m-2]
  real, dimension(n), 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(n), intent(inout) :: y !< RHS on entry, result on exit
  ! Local variables
  integer :: k, l
  real :: beta(n), lambda  ! Temporary variables in [T2 L-2 ~> s2 m-2]
  real :: I_beta(n)        ! Temporary variables in [L2 T-2 ~> m2 s-2]
  real :: yy(n)            ! A temporary variable with the same units as y on entry.
 
  lambda = lam
  beta(1) = b(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) = b(1) - lambda
  endif
  i_beta(1) = 1. / beta(1)
  yy(1) = y(1)
  do k = 2, n
    beta(k) = ( b(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. / ( b(1) - lambda )
      do l = 2, k
        i_beta(l) = 1. / ( ( b(l) - lambda ) - a(l) * c(l-1) * i_beta(l-1) )
        yy(l) = y(l) - a(l) * yy(l-1) * i_beta(l-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)
  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.
  ! 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 [Pa]
    t_int, &      ! Temperature interpolated to interfaces [degC]
    s_int, &      ! Salinity interpolated to interfaces [ppt]
    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)-1) :: &
    a_diag, b_diag, c_diag
                  ! diagonals of tridiagonal matrix; one value for each
                  ! interface (excluding surface and bottom) [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]
  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 :: 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            ! increment in lam for Newton's method [T2 L-2 ~> s2 m-2]
  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 :: min_h_frac
  real :: z_to_pa  ! A conversion factor from thicknesses (in Z) to pressure (in Pa)
  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, parameter :: reduct_factor = 0.5
                     ! 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 diffrences across interfaces times thicknesses [R Z ~> kg m-2]
  real, parameter :: tol1  = 0.0001, tol2 = 0.001
  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].
  integer :: kf(szi_(g))
  integer, parameter :: max_itt = 10
  logical :: use_eos    ! If true, density is calculated from T & S using the equation of state.
  real, dimension(SZK_(G)+1) :: z_int
  ! real, dimension(SZK_(G)+1) :: N2  ! The local squared buoyancy frequency [T-2 ~> s-2]
  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, nrows
  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
  use_eos = associated(tv%eqn_of_state)
  z_to_pa = gv%Z_to_H * gv%H_to_Pa
  c1_thresh = 0.01*us%m_s_to_L_T
 
  min_h_frac = tol1 / 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_Pa,tv,cn,g_Rho0,nmodes)
  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
          do k=2,kf(i)
            pres(k) = pres(k-1) + z_to_pa*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))
          enddo
          call calculate_density_derivs(t_int, s_int, pres, drho_dt, drho_ds, 2, &
                                        kf(i)-1, tv%eqn_of_state, scale=us%kg_m3_to_R)
 
          ! Sum the reduced gravities to find out how small a density difference
          ! is negligibly small.
          drxh_sum = 0.0
          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
        else
          drxh_sum = 0.0
          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
    !   Find gprime across each internal interface, taking care of convective
    ! instabilities by merging layers.
        if (drxh_sum <= 0.0) then
          cn(i,j,:) = 0.0
        else
          ! Merge layers to eliminate convective instabilities or exceedingly
          ! small reduced gravities.
          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 ((drho_dt(k)*(tf(k,i)-tc(kc)) + drho_ds(k)*(sf(k,i)-sc(kc))) * &
                  (hc(kc) + hf(k,i)) < 2.0 * tol2*drxh_sum) 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 ((drho_dt(k2)*(tc(k2)-tc(k2-1)) + drho_ds(k2)*(sc(k2)-sc(k2-1))) * &
                      (hc(k2) + hc(k2-1)) < tol2*drxh_sum) 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 ((rf(k,i) - rc(kc)) * (hc(kc) + hf(k,i)) < 2.0*tol2*drxh_sum) 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 ((rc(k2)-rc(k2-1)) * (hc(k2)+hc(k2-1)) < tol2*drxh_sum) 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.
          if (kc >= 2) then
            ! Set depth at surface
            z_int(1) = 0.0
            ! initialize speed2_tot
            speed2_tot = 0.0
            ! 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))
              z_int(k) = z_int(k-1) + hc(k-1)
              ! N2(K) = US%L_to_Z**2*gprime(K)/(0.5*(Hc(k)+Hc(k-1)))
              speed2_tot = speed2_tot + gprime(k)*(hc(k-1)+hc(k))
            enddo
            ! Set stratification for surface and bottom (setting equal to nearest interface for now)
            ! N2(1) = N2(2) ; N2(kc+1) = N2(kc)
            ! Calculate depth at bottom
            z_int(kc+1) = z_int(kc)+hc(kc)
            ! check that thicknesses sum to total depth
            if (abs(z_int(kc+1)-htot(i)) > 1.e-12*htot(i)) then
              call mom_error(fatal, "wave_structure: mismatch in total depths")
            endif
 
            ! Define the diagonals of the tridiagonal matrix
            ! First, populate interior rows
            do k=3,kc-1
              row = k-1 ! indexing for TD matrix rows
              a_diag(row) = -igu(k)
              b_diag(row) = igu(k)+igl(k)
              c_diag(row) = -igl(k)
            enddo
            ! Populate top row of tridiagonal matrix
            k=2 ; row = k-1
            a_diag(row) = 0.0
            b_diag(row) = igu(k)+igl(k)
            c_diag(row) = -igl(k)
            ! Populate bottom row of tridiagonal matrix
            k=kc ; row = k-1
            a_diag(row) = -igu(k)
            b_diag(row) = igu(k)+igl(k)
            c_diag(row) = 0.0
            ! Total number of rows in the matrix = number of interior interfaces
            nrows = kc-1
 
            ! Under estimate the first eigenvalue 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(a_diag(1:nrows),b_diag(1:nrows),c_diag(1:nrows), &
                                      nrows,lam_1,det,ddet, row_scale=us%m_s_to_L_T**2)
              ! 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
              else  ! Newton's method is OK.
                dlam = - det / ddet
                lam_1 = lam_1 + dlam
                if (abs(dlam) < tol2*lam_1) then
                  ! calculate 1st mode speed
                  if (lam_1 > 0.0) cn(i,j,1) = 1.0 / sqrt(lam_1)
                  exit
                endif
              endif
            enddo
 
            ! 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 + tol2)
              ! 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(a_diag(1:nrows),b_diag(1:nrows),c_diag(1:nrows), &
                               nrows,lammin,det_l,ddet_l, row_scale=us%m_s_to_L_T**2)
              ! move interval window looking for zero-crossings************************
              do iint=1,numint
                xr = lammin + laminc * iint
                xl = xr - laminc
                call tridiag_det(a_diag(1:nrows),b_diag(1:nrows),c_diag(1:nrows), &
                                 nrows,xr,det_r,ddet_r, row_scale=us%m_s_to_L_T**2)
                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(a_diag(1:nrows),b_diag(1:nrows),c_diag(1:nrows), &
                                 nrows,xl_sub,det_sub,ddet_sub, row_scale=us%m_s_to_L_T**2)
                        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(a_diag(1:nrows),b_diag(1:nrows),c_diag(1:nrows), &
                                   nrows,lam_n,det,ddet, row_scale=us%m_s_to_L_T**2)
                  ! Use Newton's method to find a new estimate of lam_n
                  dlam = - det / ddet
                  lam_n = lam_n + dlam
                  if (abs(dlam) < tol2*lam_1) then
                    ! calculate nth mode speed
                    if (lam_n > 0.0) cn(i,j,m+1) = 1.0 / sqrt(lam_n)
                    exit
                  endif ! within tol
                enddo ! itt-loop
              enddo ! n-loop
            else
              cn(i,j,2:nmodes) = 0.0 ! else too small to worry about
            endif ! if nmodes>1 .and. kc>nmodes .and. c1>c1_thresh
          else
            cn(i,j,:) = 0.0
          endif ! if more than 2 layers
        endif ! if drxh_sum < 0
      else
        cn(i,j,:) = 0.0 ! This is a land point.
      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, b, c, nrows, lam, det_out, ddet_out, row_scale)
  real, dimension(:), intent(in) :: a !< Lower diagonal of matrix (first entry = 0)
  real, dimension(:), intent(in) :: b !< Leading diagonal of matrix (excluding lam)
  real, dimension(:), intent(in) :: c !< Upper diagonal of matrix (last entry = 0)
  integer,            intent(in) :: nrows !< Size of matrix
  real,               intent(in) :: lam !< Value subtracted from b
  real,               intent(out):: det_out !< Determinant
  real,               intent(out):: ddet_out !< Derivative of determinant w.r.t. 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, dimension(nrows) :: det ! value of recursion function
  real, dimension(nrows) :: ddet ! value of recursion function for derivative
  real, parameter:: rescale = 1024.0**4 ! max value of determinant allowed before rescaling
  real :: rscl
  real :: I_rescale ! inverse of rescale
  integer :: n      ! row (layer interface) index
 
  if (size(b) /= nrows) call mom_error(warning, "Diagonal b must be same length as nrows.")
  if (size(a) /= nrows) call mom_error(warning, "Diagonal a must be same length as nrows.")
  if (size(c) /= nrows) call mom_error(warning, "Diagonal c must be same length as nrows.")
 
  i_rescale = 1.0/rescale
  rscl = 1.0 ; if (present(row_scale)) rscl = row_scale
 
  det(1) = 1.0      ; ddet(1) = 0.0
  det(2) = b(2)-lam ; ddet(2) = -1.0
  do n=3,nrows
    det(n)  = rscl*(b(n)-lam)*det(n-1)  - rscl*(a(n)*c(n-1))*det(n-2)
    ddet(n) = rscl*(b(n)-lam)*ddet(n-1) - rscl*(a(n)*c(n-1))*ddet(n-2) - det(n-1)
    ! Rescale det & ddet if det is getting too large or too small to avoid overflow or underflow.
    if (abs(det(n)) > rescale) then
      det(n)  = i_rescale*det(n)  ; det(n-1)  = i_rescale*det(n-1)
      ddet(n) = i_rescale*ddet(n) ; ddet(n-1) = i_rescale*ddet(n-1)
    elseif (abs(det(n)) < i_rescale) then
      det(n)  = rescale*det(n)  ; det(n-1)  = rescale*det(n-1)
      ddet(n) = rescale*ddet(n) ; ddet(n-1) = rescale*ddet(n-1)
    endif
  enddo
  det_out = det(nrows)
  ddet_out = ddet(nrows) / rscl
 
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)
  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].
! 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)
 
  call initialize_remapping(cs%remapping_CS, 'PLM', boundary_extrapolation=.false.)
 
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)
  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].
 
  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
 
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