mom_wave_structure Module Reference

Detailed Description

Vertical structure functions for first baroclinic mode wave speed.

Data Types
type	wave_structure_cs
	The control structure for the MOM_wave_structure module. More...

Functions/Subroutines
subroutine, public	wave_structure (h, tv, G, GV, US, cn, ModeNum, freq, CS, En, full_halos)
	This subroutine determines the internal wave velocity structure for any mode. More...
subroutine	tridiag_solver (a, b, c, h, y, method, x)
	Solves a tri-diagonal system Ax=y using either the standard Thomas algorithm (TDMA_T) or its more stable variant that invokes the "Hallberg substitution" (TDMA_H). More...
subroutine, public	wave_structure_init (Time, G, param_file, diag, CS)
	Allocate memory associated with the wave structure module and read parameters. More...

Function/Subroutine Documentation

tridiag_solver()

subroutine mom_wave_structure::tridiag_solver ( real, dimension(:), intent(in)  a, real, dimension(:), intent(in)  b, real, dimension(:), intent(in)  c, real, dimension(:), intent(in)  h, real, dimension(:), intent(in)  y, character(len=*), intent(in)  method, real, dimension(:), intent(out)  x  )

private

Solves a tri-diagonal system Ax=y using either the standard Thomas algorithm (TDMA_T) or its more stable variant that invokes the "Hallberg substitution" (TDMA_H).

Parameters

[in]	a	lower diagonal with first entry equal to zero.
[in]	b	middle diagonal.
[in]	c	upper diagonal with last entry equal to zero.
[in]	h	vector of values that have already been added to b; used for systems of the form (e.g. average layer thickness in vertical diffusion case): [ -alpha(k-1/2) ] * e(k-1) + [ alpha(k-1/2) + alpha(k+1/2) + h(k) ] * e(k) + [ -alpha(k+1/2) ] * e(k+1) = y(k) where a(k)=[-alpha(k-1/2)], b(k)=[alpha(k-1/2)+alpha(k+1/2) + h(k)], and c(k)=[-alpha(k+1/2)]. Only used with TDMA_H method.
[in]	y	vector of known values on right hand side.
[in]	method	A string describing the algorithm to use
[out]	x	vector of unknown values to solve for.

Definition at line 543 of file MOM_wave_structure.F90.

  real, dimension(:), intent(in)  :: a !< lower diagonal with first entry equal to zero.
  real, dimension(:), intent(in)  :: b !< middle diagonal.
  real, dimension(:), intent(in)  :: c !< upper diagonal with last entry equal to zero.
  real, dimension(:), intent(in)  :: h !< vector of values that have already been added to b; used
           !! for systems of the form (e.g. average layer thickness in vertical diffusion case):
           !! [ -alpha(k-1/2) ]                       * e(k-1) +
           !! [  alpha(k-1/2) + alpha(k+1/2) + h(k) ] * e(k)   +
           !! [ -alpha(k+1/2) ]                       * e(k+1) = y(k)
           !! where a(k)=[-alpha(k-1/2)], b(k)=[alpha(k-1/2)+alpha(k+1/2) + h(k)],
           !! and c(k)=[-alpha(k+1/2)]. Only used with TDMA_H method.
  real, dimension(:), intent(in)  :: y !< vector of known values on right hand side.
  character(len=*),   intent(in)  :: method !< A string describing the algorithm to use
  real, dimension(:), intent(out) :: x !< vector of unknown values to solve for.
  ! Local variables
  integer :: nrow                         ! number of rows in A matrix
!  real, allocatable, dimension(:,:) :: A_check ! for solution checking
!  real, allocatable, dimension(:)   :: y_check ! for solution checking
  real, allocatable, dimension(:) :: c_prime, y_prime, q, alpha
                                          ! intermediate values for solvers
  real    :: Q_prime, beta                ! intermediate values for solver
  integer :: k                            ! row (e.g. interface) index
  integer :: i,j
 
  nrow = size(y)
  allocate(c_prime(nrow))
  allocate(y_prime(nrow))
  allocate(q(nrow))
  allocate(alpha(nrow))
!  allocate(A_check(nrow,nrow))
!  allocate(y_check(nrow))
 
  if (method == 'TDMA_T') then
    ! Standard Thomas algoritim (4th variant).
    ! Note: Requires A to be non-singular for accuracy/stability
    c_prime(:) = 0.0       ; y_prime(:) = 0.0
    c_prime(1) = c(1)/b(1) ; y_prime(1) = y(1)/b(1)
 
    ! Forward sweep
    do k=2,nrow-1
      c_prime(k) = c(k)/(b(k)-a(k)*c_prime(k-1))
    enddo
    !print *, 'c_prime=', c_prime(1:nrow)
    do k=2,nrow
      y_prime(k) = (y(k)-a(k)*y_prime(k-1))/(b(k)-a(k)*c_prime(k-1))
    enddo
    !print *, 'y_prime=', y_prime(1:nrow)
    x(nrow) = y_prime(nrow)
 
    ! Backward sweep
    do k=nrow-1,1,-1
      x(k) = y_prime(k)-c_prime(k)*x(k+1)
    enddo
    !print *, 'x=',x(1:nrow)
 
    ! Check results - delete later
    !do j=1,nrow ; do i=1,nrow
    !  if (i==j)then ;       A_check(i,j) = b(i)
    !  elseif (i==j+1)then ; A_check(i,j) = a(i)
    !  elseif (i==j-1)then ; A_check(i,j) = c(i)
    !  endif
    !enddo ; enddo
    !print *, 'A(2,1),A(2,2),A(1,2)=', A_check(2,1), A_check(2,2), A_check(1,2)
    !y_check = matmul(A_check,x)
    !if (all(y_check /= y))then
    !  print *, "tridiag_solver: Uh oh, something's not right!"
    !  print *, "y=", y
    !  print *, "y_check=", y_check
    !endif
 
  elseif (method == 'TDMA_H') then
    ! Thomas algoritim (4th variant) w/ Hallberg substitution.
    ! For a layered system where k is at interfaces, alpha{k+1/2} refers to
    ! some property (e.g. inverse thickness for mode-structure problem) of the
    ! layer below and alpha{k-1/2} refers to the layer above.
    ! Here, alpha(k)=alpha{k+1/2} and alpha(k-1)=alpha{k-1/2}.
    ! Strictly speaking, this formulation requires A to be a non-singular,
    ! symmetric, diagonally dominant matrix, with h>0.
    ! Need to add a check for these conditions.
    do k=1,nrow-1
      if (abs(a(k+1)-c(k)) > 1.e-10*(abs(a(k+1))+abs(c(k)))) then
        call mom_error(warning, "tridiag_solver: matrix not symmetric; need symmetry when invoking TDMA_H")
      endif
    enddo
    alpha = -c
    ! Alpha of the bottom-most layer is not necessarily zero. Therefore,
    ! back out the value from the provided b(nrow and h(nrow) values
    alpha(nrow)   = b(nrow)-h(nrow)-alpha(nrow-1)
    ! Prime other variables
    beta       = 1/b(1)
    y_prime(:) = 0.0       ; q(:) = 0.0
    y_prime(1) = beta*y(1) ; q(1) = beta*alpha(1)
    q_prime    = 1-q(1)
 
    ! Forward sweep
    do k=2,nrow-1
      beta = 1/(h(k)+alpha(k-1)*q_prime+alpha(k))
      if (isnan(beta))then ; print *, "Tridiag_solver: beta is NAN" ; endif
      q(k) = beta*alpha(k)
      y_prime(k) = beta*(y(k)+alpha(k-1)*y_prime(k-1))
      q_prime = beta*(h(k)+alpha(k-1)*q_prime)
    enddo
    if ((h(nrow)+alpha(nrow-1)*q_prime+alpha(nrow)) == 0.0)then
      call mom_error(fatal, "Tridiag_solver: this system is not stable.") ! ; overriding beta(nrow)
      ! This has hard-coded dimensions: beta = 1/(1e-15) ! place holder for unstable systems - delete later
    else
      beta = 1/(h(nrow)+alpha(nrow-1)*q_prime+alpha(nrow))
    endif
    y_prime(nrow) = beta*(y(nrow)+alpha(nrow-1)*y_prime(nrow-1))
    x(nrow) = y_prime(nrow)
    ! Backward sweep
    do k=nrow-1,1,-1
      x(k) = y_prime(k)+q(k)*x(k+1)
    enddo
    !print *, 'yprime=',y_prime(1:nrow)
    !print *, 'x=',x(1:nrow)
  endif
 
  deallocate(c_prime,y_prime,q,alpha)
! deallocate(A_check,y_check)
 

wave_structure()

subroutine, public mom_wave_structure::wave_structure	(	real, dimension( g %isd: g %ied, g %jsd: g %jed, g %ke), intent(in)	h,
		type(thermo_var_ptrs), intent(in)	tv,
		type(ocean_grid_type), intent(in)	G,
		type(verticalgrid_type), intent(in)	GV,
		type(unit_scale_type), intent(in)	US,
		real, dimension( g %isd: g %ied, g %jsd: g %jed), intent(in)	cn,
		integer, intent(in)	ModeNum,
		real, intent(in)	freq,
		type(wave_structure_cs), pointer	CS,
		real, dimension( g %isd: g %ied, g %jsd: g %jed), intent(in), optional	En,
		logical, intent(in), optional	full_halos
	)

This subroutine determines the internal wave velocity structure for any mode.

This subroutine solves for the eigen vector [vertical structure, e(k)] associated with the first baroclinic mode speed [i.e., smallest eigen value (lam = 1/c^2)] of the system d2e/dz2 = -(N2/cn2)e, or (A-lam*I)e = 0, where A = -(1/N2)(d2/dz2), lam = 1/c^2, and I is the identity matrix. 2nd order discretization in the vertical lets this system be represented as

-Igu(k)*e(k-1) + (Igu(k)+Igl(k)-lam)*e(k) - Igl(k)*e(k+1) = 0.0

with rigid lid boundary conditions e(1) = e(nz+1) = 0.0 giving

(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

where, upon noting N2 = reduced gravity/layer thickness, we get Igl(k) = 1.0/(gprime(k)*H(k)) ; Igu(k) = 1.0/(gprime(k)*H(k-1))

The eigen value for this system is approximated using "wave_speed." This subroutine uses these eigen values (mode speeds) to estimate the corresponding eigen vectors (velocity structure) using the "inverse iteration with shift" method. The algorithm is

Pick a starting vector reasonably close to mode structure and with unit magnitude, b_guess For n=1,2,3,... Solve (A-lam*I)e = e_guess for e Set e_guess=e/|e| and repeat, with each iteration refining the estimate of e

Parameters

[in]	g	The ocean's grid structure.
[in]	gv	The ocean's vertical grid structure.
[in]	us	A dimensional unit scaling type
[in]	h	Layer thicknesses [H ~> m or kg m-2]
[in]	tv	A structure pointing to various thermodynamic variables.
[in]	cn	The (non-rotational) mode internal gravity wave speed [m s-1].
[in]	modenum	Mode number
[in]	freq	Intrinsic wave frequency [s-1].
	cs	The control structure returned by a previous call to wave_structure_init.
[in]	en	Internal wave energy density [J m-2].
[in]	full_halos	If true, do the calculation over the entire computational domain.

Definition at line 92 of file MOM_wave_structure.F90.

  type(ocean_grid_type),                    intent(in)  :: G  !< The ocean's grid structure.
  type(verticalGrid_type),                  intent(in)  :: GV !< The ocean's 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 thicknesses [H ~> m or kg m-2]
  type(thermo_var_ptrs),                    intent(in)  :: tv !< A structure pointing to various
                                                              !! thermodynamic variables.
  real, dimension(SZI_(G),SZJ_(G)),         intent(in)  :: cn !< The (non-rotational) mode internal
                                                              !! gravity wave speed [m s-1].
  integer,                                  intent(in)  :: ModeNum !< Mode number
  real,                                     intent(in)  :: freq !< Intrinsic wave frequency [s-1].
  type(wave_structure_CS),                  pointer     :: CS !< The control structure returned by a
                                                              !! previous call to wave_structure_init.
  real, dimension(SZI_(G),SZJ_(G)), &
                                  optional, intent(in)  :: En !< Internal wave energy density [J m-2].
  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, dRho_dS, &
    pres, T_int, S_int, &
    gprime        ! The reduced gravity across each interface [m2 Z-1 s-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 [s2 m-2].
  real, dimension(SZK_(G),SZI_(G)) :: &
    Hf, Tf, Sf, Rf
  real, dimension(SZK_(G)) :: &
    Hc, Tc, Sc, Rc, &
    det, ddet
  real, dimension(SZI_(G),SZJ_(G)) :: &
    htot
  real :: lam
  real :: min_h_frac
  real :: H_to_pres
  real, dimension(SZI_(G)) :: &
    hmin, &  ! Thicknesses [Z ~> m].
    H_here, HxT_here, HxS_here, HxR_here
  real :: speed2_tot
  real :: I_Hnew, drxh_sum
  real, parameter :: tol1  = 0.0001, tol2 = 0.001
  real, pointer, dimension(:,:,:) :: T => null(), s => null()
  real :: g_Rho0  ! G_Earth/Rho0 in m5 Z-1 s-2 kg-1.
  real :: rescale, I_rescale
  integer :: kf(SZI_(G))
  integer, parameter :: max_itt = 1 ! number of times to iterate in solving for eigenvector
  real, parameter    :: cg_subRO = 1e-100 ! a very small number
  real, parameter    :: a_int = 0.5 ! value of normalized integral: \int(w_strct^2)dz = a_int
  real               :: I_a_int     ! inverse of a_int
  real               :: f2          ! squared Coriolis frequency
  real               :: Kmag2       ! magnitude of horizontal wave number squared
  logical            :: use_EOS     ! If true, density is calculated from T & S using an
                                    ! equation of state.
  real, dimension(SZK_(G)+1) :: w_strct, u_strct, W_profile, Uavg_profile, z_int, N2
                                        ! local representations of variables in CS; note,
                                        ! not all rows will be filled if layers get merged!
  real, dimension(SZK_(G)+1) :: w_strct2, u_strct2
                                        ! squared values
  real, dimension(SZK_(G))   :: dz      ! thicknesses of merged layers (same as Hc I hope)
  real, dimension(SZK_(G)+1) :: dWdz_profile ! profile of dW/dz
  real                       :: w2avg   ! average of squared vertical velocity structure funtion
  real                       :: int_dwdz2, int_w2, int_N2w2, KE_term, PE_term, W0
                                        ! terms in vertically averaged energy equation
  real                       :: gp_unscaled ! A version of gprime rescaled to [m s-2].
  real, dimension(SZK_(G)-1) :: lam_z   ! product of eigen value and gprime(k); one value for each
                                        ! interface (excluding surface and bottom)
  real, dimension(SZK_(G)-1) :: a_diag, b_diag, c_diag
                                        ! diagonals of tridiagonal matrix; one value for each
                                        ! interface (excluding surface and bottom)
  real, dimension(SZK_(G)-1) :: e_guess ! guess at eigen vector with unit amplitde (for TDMA)
  real, dimension(SZK_(G)-1) :: e_itt   ! improved guess at eigen vector (from TDMA)
  real    :: Pi
  integer :: kc
  integer :: i, j, k, k2, itt, is, ie, js, je, nz, nzm, row, ig, jg, ig_stop, jg_stop
 
  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
  i_a_int = 1/a_int
 
  !if (present(CS)) then
    if (.not. associated(cs)) call mom_error(fatal, "MOM_wave_structure: "// &
           "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
 
  pi = (4.0*atan(1.0))
 
  s => tv%S ; t => tv%T
  g_rho0 = us%L_T_to_m_s**2 * gv%g_Earth /gv%Rho0
  use_eos = associated(tv%eqn_of_state)
 
  h_to_pres = gv%Z_to_H*gv%H_to_Pa
  rescale = 1024.0**4 ; i_rescale = 1.0/rescale
 
  min_h_frac = tol1 / real(nz)
 
  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,j) = 0.0 ; enddo
    do k=1,nz ; do i=is,ie ; htot(i,j) = htot(i,j) + h(i,j,k)*gv%H_to_Z ; enddo ; enddo
 
    do i=is,ie
      hmin(i) = htot(i,j)*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 ! use_EOS?
 
    ! From this point, we can work on individual columns without causing memory
    ! to have page faults.
    do i=is,ie ; if (cn(i,j)>0.0)then
      !----for debugging, remove later----
      ig = i + g%idg_offset ; jg = j + g%jdg_offset
      !if (ig == CS%int_tide_source_x .and. jg == CS%int_tide_source_y) then
      !-----------------------------------
      if (g%mask2dT(i,j) > 0.5) then
 
        lam = 1/(cn(i,j)**2)
 
        ! Calculate drxh_sum
        if (use_eos) then
          pres(1) = 0.0
          do k=2,kf(i)
            pres(k) = pres(k-1) + h_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))
          enddo
          call calculate_density_derivs(t_int, s_int, pres, drho_dt, drho_ds, 2, &
                                        kf(i)-1, tv%eqn_of_state)
 
          ! 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 ! use_EOS?
 
        !   Find gprime across each internal interface, taking care of convective
        ! instabilities by merging layers.
        if (drxh_sum >= 0.0) then
          ! 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 STRUCTURE-------------------------------------
          ! Construct and solve tridiagonal system for the interior interfaces
          ! Note that kc   = number of layers,
          !           kc+1 = nzm = number of interfaces,
          !           kc-1 = number of interior interfaces (excluding surface and bottom)
          ! Also, note that "K" refers to an interface, while "k" refers to the layer below.
          ! Need at least 3 layers (2 internal interfaces) to generate a matrix, also
          ! need number of layers to be greater than the mode number
          if (kc >= modenum + 1) then
            ! Set depth at surface
            z_int(1) = 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%m_to_Z**2*gprime(k)/(0.5*(hc(k)+hc(k-1)))
            enddo
            ! Set stratification for surface and bottom (setting equal to nearest interface for now)
            n2(1) = n2(2) ; n2(kc+1) = n2(kc)
            ! Calcualte 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,j)) > 1.e-14*htot(i,j)) then
              call mom_error(fatal, "wave_structure: mismatch in total depths")
            endif
 
            ! Note that many of the calcluation from here on revert to using vertical
            ! distances in m, not Z.
 
            ! Populate interior rows of tridiagonal matrix; must multiply through by
            ! gprime to get tridiagonal matrix to the symmetrical form:
            ! [-1/H(k-1)]e(k-1) + [1/H(k-1)+1/H(k)-lam_z]e(k) + [-1/H(k)]e(k+1) = 0,
            ! where lam_z = lam*gprime is now a function of depth.
            ! Frist, populate interior rows
            do k=3,kc-1
              row = k-1 ! indexing for TD matrix rows
              gp_unscaled = us%m_to_Z*gprime(k)
              lam_z(row) = lam*gp_unscaled
              a_diag(row) = gp_unscaled*(-igu(k))
              b_diag(row) = gp_unscaled*(igu(k)+igl(k)) - lam_z(row)
              c_diag(row) = gp_unscaled*(-igl(k))
              if (isnan(lam_z(row)))then  ; print *, "Wave_structure: lam_z(row) is NAN" ; endif
              if (isnan(a_diag(row)))then ; print *, "Wave_structure: a(k) is NAN" ; endif
              if (isnan(b_diag(row)))then ; print *, "Wave_structure: b(k) is NAN" ; endif
              if (isnan(c_diag(row)))then ; print *, "Wave_structure: c(k) is NAN" ; endif
            enddo
            ! Populate top row of tridiagonal matrix
            k=2 ; row = k-1 ;
            gp_unscaled = us%m_to_Z*gprime(k)
            lam_z(row) = lam*gp_unscaled
            a_diag(row) = 0.0
            b_diag(row) = gp_unscaled*(igu(k)+igl(k)) - lam_z(row)
            c_diag(row) = gp_unscaled*(-igl(k))
            ! Populate bottom row of tridiagonal matrix
            k=kc ; row = k-1
            gp_unscaled = us%m_to_Z*gprime(k)
            lam_z(row) = lam*gp_unscaled
            a_diag(row) = gp_unscaled*(-igu(k))
            b_diag(row) = gp_unscaled*(igu(k)+igl(k)) - lam_z(row)
            c_diag(row) = 0.0
 
            ! Guess a vector shape to start with (excludes surface and bottom)
            e_guess(1:kc-1) = sin((z_int(2:kc)/htot(i,j)) *pi)
            e_guess(1:kc-1) = e_guess(1:kc-1)/sqrt(sum(e_guess(1:kc-1)**2))
 
            ! Perform inverse iteration with tri-diag solver
            do itt=1,max_itt
              call tridiag_solver(a_diag(1:kc-1),b_diag(1:kc-1),c_diag(1:kc-1), &
                                  -lam_z(1:kc-1),e_guess(1:kc-1),"TDMA_H",e_itt)
              e_guess(1:kc-1) = e_itt(1:kc-1)/sqrt(sum(e_itt(1:kc-1)**2))
            enddo ! itt-loop
            w_strct(2:kc) = e_guess(1:kc-1)
            w_strct(1)    = 0.0 ! rigid lid at surface
            w_strct(kc+1) = 0.0 ! zero-flux at bottom
 
            ! Check to see if solver worked
            ig_stop = 0 ; jg_stop = 0
            if (isnan(sum(w_strct(1:kc+1))))then
              print *, "Wave_structure: w_strct has a NAN at ig=", ig, ", jg=", jg
              if (i<g%isc .or. i>g%iec .or. j<g%jsc .or. j>g%jec)then
                print *, "This is occuring at a halo point."
              endif
              ig_stop = ig ; jg_stop = jg
            endif
 
            ! Normalize vertical structure function of w such that
            ! \int(w_strct)^2dz = a_int (a_int could be any value, e.g., 0.5)
            nzm = kc+1 ! number of layer interfaces after merging
                       !(including surface and bottom)
            w2avg = 0.0
            do k=1,nzm-1
              dz(k) = us%Z_to_m*hc(k)
              w2avg = w2avg + 0.5*(w_strct(k)**2+w_strct(k+1)**2)*dz(k)
            enddo
            !### Some mathematical cancellations could occur in the next two lines.
            w2avg = w2avg / htot(i,j)
            w_strct = w_strct / sqrt(htot(i,j)*w2avg*i_a_int)
 
            ! Calculate vertical structure function of u (i.e. dw/dz)
            do k=2,nzm-1
              u_strct(k) = 0.5*((w_strct(k-1) - w_strct(k)  )/dz(k-1) + &
                           (w_strct(k)   - w_strct(k+1))/dz(k))
            enddo
            u_strct(1)   = (w_strct(1)   -  w_strct(2) )/dz(1)
            u_strct(nzm)  = (w_strct(nzm-1)-  w_strct(nzm))/dz(nzm-1)
 
            ! Calculate wavenumber magnitude
            f2 = us%s_to_T**2 * g%CoriolisBu(i,j)**2
            !f2 = 0.25*US%s_to_T**2 *((G%CoriolisBu(I,J)**2 + G%CoriolisBu(I-1,J-1)**2) + &
            !    (G%CoriolisBu(I,J-1)**2 + G%CoriolisBu(I-1,J)**2))
            kmag2 = (freq**2 - f2) / (cn(i,j)**2 + cg_subro**2)
 
            ! Calculate terms in vertically integrated energy equation
            int_dwdz2 = 0.0 ; int_w2 = 0.0 ; int_n2w2 = 0.0
            u_strct2 = u_strct(1:nzm)**2
            w_strct2 = w_strct(1:nzm)**2
            ! vertical integration with Trapezoidal rule
            do k=1,nzm-1
              int_dwdz2 = int_dwdz2 + 0.5*(u_strct2(k)+u_strct2(k+1))*dz(k)
              int_w2    = int_w2    + 0.5*(w_strct2(k)+w_strct2(k+1))*dz(k)
              int_n2w2  = int_n2w2  + 0.5*(w_strct2(k)*n2(k)+w_strct2(k+1)*n2(k+1))*dz(k)
            enddo
 
            ! Back-calculate amplitude from energy equation
            if (kmag2 > 0.0) then
              ke_term = 0.25*gv%Rho0*( (1+f2/freq**2)/kmag2*int_dwdz2 + int_w2 )
              pe_term = 0.25*gv%Rho0*( int_n2w2/freq**2 )
              if (en(i,j) >= 0.0) then
                w0 = sqrt( en(i,j)/(ke_term + pe_term) )
              else
                call mom_error(warning, "wave_structure: En < 0.0; setting to W0 to 0.0")
                print *, "En(i,j)=", en(i,j), " at ig=", ig, ", jg=", jg
                w0 = 0.0
              endif
              ! Calculate actual vertical velocity profile and derivative
              w_profile    = w0*w_strct
              dwdz_profile = w0*u_strct
              ! Calculate average magnitude of actual horizontal velocity over a period
              uavg_profile = abs(dwdz_profile) * sqrt((1+f2/freq**2)/(2.0*kmag2))
            else
              w_profile    = 0.0
              dwdz_profile = 0.0
              uavg_profile = 0.0
            endif
 
            ! Store values in control structure
            cs%w_strct(i,j,1:nzm)     = w_strct(:)
            cs%u_strct(i,j,1:nzm)     = u_strct(:)
            cs%W_profile(i,j,1:nzm)   = w_profile(:)
            cs%Uavg_profile(i,j,1:nzm)= uavg_profile(:)
            cs%z_depths(i,j,1:nzm)    = us%Z_to_m*z_int(:)
            cs%N2(i,j,1:nzm)          = n2(:)
            cs%num_intfaces(i,j)      = nzm
          else
            ! If not enough layers, default to zero
            nzm = kc+1
            cs%w_strct(i,j,1:nzm)     = 0.0
            cs%u_strct(i,j,1:nzm)     = 0.0
            cs%W_profile(i,j,1:nzm)   = 0.0
            cs%Uavg_profile(i,j,1:nzm)= 0.0
            cs%z_depths(i,j,1:nzm)    = 0.0 ! could use actual values
            cs%N2(i,j,1:nzm)          = 0.0 ! could use with actual values
            cs%num_intfaces(i,j)       = nzm
          endif  ! kc >= 3 and kc > ModeNum + 1?
        endif ! drxh_sum >= 0?
      !else     ! if at test point - delete later
      !  return ! if at test point - delete later
      !endif    ! if at test point - delete later
      endif ! mask2dT > 0.5?
    else
      ! if cn=0.0, default to zero
      nzm                       = nz+1! could use actual values
      cs%w_strct(i,j,1:nzm)     = 0.0
      cs%u_strct(i,j,1:nzm)     = 0.0
      cs%W_profile(i,j,1:nzm)   = 0.0
      cs%Uavg_profile(i,j,1:nzm)= 0.0
      cs%z_depths(i,j,1:nzm)    = 0.0 ! could use actual values
      cs%N2(i,j,1:nzm)          = 0.0 ! could use with actual values
      cs%num_intfaces(i,j)       = nzm
    endif ; enddo ! if cn>0.0? ; i-loop
  enddo ! j-loop
 

wave_structure_init()

subroutine, public mom_wave_structure::wave_structure_init	(	type(time_type), intent(in)	Time,
		type(ocean_grid_type), intent(in)	G,
		type(param_file_type), intent(in)	param_file,
		type(diag_ctrl), intent(in), target	diag,
		type(wave_structure_cs), pointer	CS
	)

Allocate memory associated with the wave structure module and read parameters.

Parameters

[in]	time	The current model time.
[in]	g	The ocean's grid structure.
[in]	param_file	A structure to parse for run-time parameters.
[in]	diag	A structure that is used to regulate diagnostic output.
	cs	A pointer that is set to point to the control structure for this module.

Definition at line 667 of file MOM_wave_structure.F90.

  type(time_type),         intent(in) :: Time !< The current model time.
  type(ocean_grid_type),   intent(in) :: G    !< The ocean's grid structure.
  type(param_file_type),   intent(in) :: param_file !< A structure to parse for run-time
                                              !! parameters.
  type(diag_ctrl), target, intent(in) :: diag !< A structure that is used to regulate
                                              !! diagnostic output.
  type(wave_structure_CS), pointer    :: CS   !< A pointer that is set to point to the
                                              !! control structure for this module.
! This include declares and sets the variable "version".
#include "version_variable.h"
  character(len=40)  :: mdl = "MOM_wave_structure"  ! This module's name.
  integer :: isd, ied, jsd, jed, nz
 
  isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed ; nz = g%ke
 
  if (associated(cs)) then
    call mom_error(warning, "wave_structure_init called with an "// &
                            "associated control structure.")
    return
  else ; allocate(cs) ; endif
 
  call get_param(param_file, mdl, "INTERNAL_TIDE_SOURCE_X", cs%int_tide_source_x, &
                 "X Location of generation site for internal tide", default=1.)
  call get_param(param_file, mdl, "INTERNAL_TIDE_SOURCE_Y", cs%int_tide_source_y, &
                 "Y Location of generation site for internal tide", default=1.)
 
  cs%diag => diag
 
  ! Allocate memory for variable in control structure; note,
  ! not all rows will be filled if layers get merged!
  allocate(cs%w_strct(isd:ied,jsd:jed,nz+1))
  allocate(cs%u_strct(isd:ied,jsd:jed,nz+1))
  allocate(cs%W_profile(isd:ied,jsd:jed,nz+1))
  allocate(cs%Uavg_profile(isd:ied,jsd:jed,nz+1))
  allocate(cs%z_depths(isd:ied,jsd:jed,nz+1))
  allocate(cs%N2(isd:ied,jsd:jed,nz+1))
  allocate(cs%num_intfaces(isd:ied,jsd:jed))
 
  ! Write all relevant parameters to the model log.
  call log_version(param_file, mdl, version, "")