MOM6
Data Types | Functions/Subroutines
mom_wave_speed Module Reference

Detailed Description

Routines for calculating baroclinic wave speeds.

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

   -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

Here 

   Igl(k) = 1.0/(gprime(k)*h(k)) ; Igu(k) = 1.0/(gprime(k)*h(k-1))

Alternately, these same eigenvalues can be found from the second smallest eigenvalue of the Montgomery potential (M(k)) calculation:

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

with rigid lid and flat bottom boundary conditions

   (Igu(2)-lam)*M(1) - Igu(2)*M(2) = 0.0
   -Igl(nz)*M(nz-1) + (Igl(nz)-lam)*M(nz) = 0.0

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.

Data Types

type  wave_speed_cs
 Control structure for MOM_wave_speed. More...
 

Functions/Subroutines

subroutine, public wave_speed (h, tv, G, GV, US, cg1, CS, full_halos, use_ebt_mode, mono_N2_column_fraction, mono_N2_depth, modal_structure)
 Calculates the wave speed of the first baroclinic mode. More...
 
subroutine tdma6 (n, a, b, c, lam, y)
 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. More...
 
subroutine, public wave_speeds (h, tv, G, GV, US, nmodes, cn, CS, full_halos)
 Calculates the wave speeds for the first few barolinic modes. More...
 
subroutine tridiag_det (a, b, c, nrows, lam, det_out, ddet_out, row_scale)
 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. More...
 
subroutine, public wave_speed_init (CS, use_ebt_mode, mono_N2_column_fraction, mono_N2_depth)
 Initialize control structure for MOM_wave_speed. More...
 
subroutine, public wave_speed_set_param (CS, use_ebt_mode, mono_N2_column_fraction, mono_N2_depth)
 Sets internal parameters for MOM_wave_speed. More...
 

Function/Subroutine Documentation

◆ tdma6()

subroutine mom_wave_speed::tdma6 ( integer, intent(in)  n,
real, dimension(n), intent(in)  a,
real, dimension(n), intent(in)  b,
real, dimension(n), intent(in)  c,
real, intent(in)  lam,
real, dimension(n), intent(inout)  y 
)
private

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.

Parameters
[in]nNumber of rows of matrix
[in]aLower diagonal [T2 L-2 ~> s2 m-2]
[in]bLeading diagonal [T2 L-2 ~> s2 m-2]
[in]cUpper diagonal [T2 L-2 ~> s2 m-2]
[in]lamScalar subtracted from leading diagonal [T2 L-2 ~> s2 m-2]
[in,out]yRHS on entry, result on exit

Definition at line 515 of file MOM_wave_speed.F90.

515  integer, intent(in) :: n !< Number of rows of matrix
516  real, dimension(n), intent(in) :: a !< Lower diagonal [T2 L-2 ~> s2 m-2]
517  real, dimension(n), intent(in) :: b !< Leading diagonal [T2 L-2 ~> s2 m-2]
518  real, dimension(n), intent(in) :: c !< Upper diagonal [T2 L-2 ~> s2 m-2]
519  real, intent(in) :: lam !< Scalar subtracted from leading diagonal [T2 L-2 ~> s2 m-2]
520  real, dimension(n), intent(inout) :: y !< RHS on entry, result on exit
521  ! Local variables
522  integer :: k, l
523  real :: beta(n), lambda ! Temporary variables in [T2 L-2 ~> s2 m-2]
524  real :: I_beta(n) ! Temporary variables in [L2 T-2 ~> m2 s-2]
525  real :: yy(n) ! A temporary variable with the same units as y on entry.
526 
527  lambda = lam
528  beta(1) = b(1) - lambda
529  if (beta(1)==0.) then ! lam was chosen too perfectly
530  ! Change lambda and redo this first row
531  lambda = (1. + 1.e-5) * lambda
532  beta(1) = b(1) - lambda
533  endif
534  i_beta(1) = 1. / beta(1)
535  yy(1) = y(1)
536  do k = 2, n
537  beta(k) = ( b(k) - lambda ) - a(k) * c(k-1) * i_beta(k-1)
538  ! Perhaps the following 0 needs to become a tolerance to handle underflow?
539  if (beta(k)==0.) then ! lam was chosen too perfectly
540  ! Change lambda and redo everything up to row k
541  lambda = (1. + 1.e-5) * lambda
542  i_beta(1) = 1. / ( b(1) - lambda )
543  do l = 2, k
544  i_beta(l) = 1. / ( ( b(l) - lambda ) - a(l) * c(l-1) * i_beta(l-1) )
545  yy(l) = y(l) - a(l) * yy(l-1) * i_beta(l-1)
546  enddo
547  else
548  i_beta(k) = 1. / beta(k)
549  endif
550  yy(k) = y(k) - a(k) * yy(k-1) * i_beta(k-1)
551  enddo
552  ! The units of y change by a factor of [L2 T-2] in the following lines.
553  y(n) = yy(n) * i_beta(n)
554  do k = n-1, 1, -1
555  y(k) = ( yy(k) - c(k) * y(k+1) ) * i_beta(k)
556  enddo

Referenced by wave_speed().

Here is the caller graph for this function:

◆ tridiag_det()

subroutine mom_wave_speed::tridiag_det ( real, dimension(:), intent(in)  a,
real, dimension(:), intent(in)  b,
real, dimension(:), intent(in)  c,
integer, intent(in)  nrows,
real, intent(in)  lam,
real, intent(out)  det_out,
real, intent(out)  ddet_out,
real, intent(in), optional  row_scale 
)
private

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.

Parameters
[in]aLower diagonal of matrix (first entry = 0)
[in]bLeading diagonal of matrix (excluding lam)
[in]cUpper diagonal of matrix (last entry = 0)
[in]nrowsSize of matrix
[in]lamValue subtracted from b
[out]det_outDeterminant
[out]ddet_outDerivative of determinant w.r.t. lam
[in]row_scaleA scaling factor of the rows of the matrix to limit the growth of the determinant

Definition at line 1016 of file MOM_wave_speed.F90.

1016  real, dimension(:), intent(in) :: a !< Lower diagonal of matrix (first entry = 0)
1017  real, dimension(:), intent(in) :: b !< Leading diagonal of matrix (excluding lam)
1018  real, dimension(:), intent(in) :: c !< Upper diagonal of matrix (last entry = 0)
1019  integer, intent(in) :: nrows !< Size of matrix
1020  real, intent(in) :: lam !< Value subtracted from b
1021  real, intent(out):: det_out !< Determinant
1022  real, intent(out):: ddet_out !< Derivative of determinant w.r.t. lam
1023  real, optional, intent(in) :: row_scale !< A scaling factor of the rows of the
1024  !! matrix to limit the growth of the determinant
1025  ! Local variables
1026  real, dimension(nrows) :: det ! value of recursion function
1027  real, dimension(nrows) :: ddet ! value of recursion function for derivative
1028  real, parameter:: rescale = 1024.0**4 ! max value of determinant allowed before rescaling
1029  real :: rscl
1030  real :: I_rescale ! inverse of rescale
1031  integer :: n ! row (layer interface) index
1032 
1033  if (size(b) /= nrows) call mom_error(warning, "Diagonal b must be same length as nrows.")
1034  if (size(a) /= nrows) call mom_error(warning, "Diagonal a must be same length as nrows.")
1035  if (size(c) /= nrows) call mom_error(warning, "Diagonal c must be same length as nrows.")
1036 
1037  i_rescale = 1.0/rescale
1038  rscl = 1.0 ; if (present(row_scale)) rscl = row_scale
1039 
1040  det(1) = 1.0 ; ddet(1) = 0.0
1041  det(2) = b(2)-lam ; ddet(2) = -1.0
1042  do n=3,nrows
1043  det(n) = rscl*(b(n)-lam)*det(n-1) - rscl*(a(n)*c(n-1))*det(n-2)
1044  ddet(n) = rscl*(b(n)-lam)*ddet(n-1) - rscl*(a(n)*c(n-1))*ddet(n-2) - det(n-1)
1045  ! Rescale det & ddet if det is getting too large or too small to avoid overflow or underflow.
1046  if (abs(det(n)) > rescale) then
1047  det(n) = i_rescale*det(n) ; det(n-1) = i_rescale*det(n-1)
1048  ddet(n) = i_rescale*ddet(n) ; ddet(n-1) = i_rescale*ddet(n-1)
1049  elseif (abs(det(n)) < i_rescale) then
1050  det(n) = rescale*det(n) ; det(n-1) = rescale*det(n-1)
1051  ddet(n) = rescale*ddet(n) ; ddet(n-1) = rescale*ddet(n-1)
1052  endif
1053  enddo
1054  det_out = det(nrows)
1055  ddet_out = ddet(nrows) / rscl
1056 

References mom_error_handler::mom_error().

Referenced by wave_speeds().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ wave_speed()

subroutine, public mom_wave_speed::wave_speed ( 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(out)  cg1,
type(wave_speed_cs), pointer  CS,
logical, intent(in), optional  full_halos,
logical, intent(in), optional  use_ebt_mode,
real, intent(in), optional  mono_N2_column_fraction,
real, intent(in), optional  mono_N2_depth,
real, dimension( g %isd: g %ied, g %jsd: g %jed, g %ke), intent(out), optional  modal_structure 
)

Calculates the wave speed of the first baroclinic mode.

Parameters
[in]gOcean grid structure
[in]gvVertical grid structure
[in]usA dimensional unit scaling type
[in]hLayer thickness [H ~> m or kg m-2]
[in]tvThermodynamic variables
[out]cg1First mode internal wave speed [L T-1 ~> m s-1]
csControl structure for MOM_wave_speed
[in]full_halosIf true, do the calculation over the entire computational domain.
[in]use_ebt_modeIf true, use the equivalent barotropic mode instead of the first baroclinic mode.
[in]mono_n2_column_fractionThe lower fraction of water column over which N2 is limited as monotonic for the purposes of calculating vertical modal structure.
[in]mono_n2_depthA depth below which N2 is limited as monotonic for the purposes of calculating vertical modal structure [Z ~> m].
[out]modal_structureNormalized model structure [nondim]

Definition at line 51 of file MOM_wave_speed.F90.

51  type(ocean_grid_type), intent(in) :: G !< Ocean grid structure
52  type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure
53  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
54  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
55  intent(in) :: h !< Layer thickness [H ~> m or kg m-2]
56  type(thermo_var_ptrs), intent(in) :: tv !< Thermodynamic variables
57  real, dimension(SZI_(G),SZJ_(G)), intent(out) :: cg1 !< First mode internal wave speed [L T-1 ~> m s-1]
58  type(wave_speed_CS), pointer :: CS !< Control structure for MOM_wave_speed
59  logical, optional, intent(in) :: full_halos !< If true, do the calculation
60  !! over the entire computational domain.
61  logical, optional, intent(in) :: use_ebt_mode !< If true, use the equivalent
62  !! barotropic mode instead of the first baroclinic mode.
63  real, optional, intent(in) :: mono_N2_column_fraction !< The lower fraction
64  !! of water column over which N2 is limited as monotonic
65  !! for the purposes of calculating vertical modal structure.
66  real, optional, intent(in) :: mono_N2_depth !< A depth below which N2 is limited as
67  !! monotonic for the purposes of calculating vertical
68  !! modal structure [Z ~> m].
69  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
70  optional, intent(out) :: modal_structure !< Normalized model structure [nondim]
71 
72  ! Local variables
73  real, dimension(SZK_(G)+1) :: &
74  dRho_dT, & ! Partial derivative of density with temperature [R degC-1 ~> kg m-3 degC-1]
75  dRho_dS, & ! Partial derivative of density with salinity [R ppt-1 ~> kg m-3 ppt-1]
76  pres, & ! Interface pressure [Pa]
77  T_int, & ! Temperature interpolated to interfaces [degC]
78  S_int, & ! Salinity interpolated to interfaces [ppt]
79  gprime ! The reduced gravity across each interface [L2 Z-1 T-2 ~> m s-2].
80  real, dimension(SZK_(G)) :: &
81  Igl, Igu, Igd ! The inverse of the reduced gravity across an interface times
82  ! the thickness of the layer below (Igl) or above (Igu) it.
83  ! Their sum, Igd, is provided for the tridiagonal solver. [T2 L-2 ~> s2 m-2]
84  real, dimension(SZK_(G),SZI_(G)) :: &
85  Hf, & ! Layer thicknesses after very thin layers are combined [Z ~> m]
86  Tf, & ! Layer temperatures after very thin layers are combined [degC]
87  Sf, & ! Layer salinities after very thin layers are combined [ppt]
88  Rf ! Layer densities after very thin layers are combined [R ~> kg m-3]
89  real, dimension(SZK_(G)) :: &
90  Hc, & ! A column of layer thicknesses after convective istabilities are removed [Z ~> m]
91  Tc, & ! A column of layer temperatures after convective istabilities are removed [degC]
92  Sc, & ! A column of layer salinites after convective istabilities are removed [ppt]
93  Rc, & ! A column of layer densities after convective istabilities are removed [R ~> kg m-3]
94  Hc_H ! Hc(:) rescaled from Z to thickness units [H ~> m or kg m-2]
95  real :: det, ddet, detKm1, detKm2, ddetKm1, ddetKm2
96  real :: lam ! The eigenvalue [T2 L-2 ~> s m-1]
97  real :: dlam ! The change in estimates of the eigenvalue [T2 L-2 ~> s m-1]
98  real :: lam0 ! The first guess of the eigenvalue [T2 L-2 ~> s m-1]
99  real :: min_h_frac ! [nondim]
100  real :: Z_to_Pa ! A conversion factor from thicknesses (in Z) to pressure (in Pa)
101  real, dimension(SZI_(G)) :: &
102  htot, hmin, & ! Thicknesses [Z ~> m]
103  H_here, & ! A thickness [Z ~> m]
104  HxT_here, & ! A layer integrated temperature [degC Z ~> degC m]
105  HxS_here, & ! A layer integrated salinity [ppt Z ~> ppt m]
106  HxR_here ! A layer integrated density [R Z ~> kg m-2]
107  real :: speed2_tot ! overestimate of the mode-1 speed squared [L2 T-2 ~> m2 s-2]
108  real :: I_Hnew ! The inverse of a new layer thickness [Z-1 ~> m-1]
109  real :: drxh_sum ! The sum of density diffrences across interfaces times thicknesses [R Z ~> kg m-2]
110  real :: L2_to_Z2 ! A scaling factor squared from units of lateral distances to depths [Z2 L-2 ~> 1].
111  real, parameter :: tol1 = 0.0001, tol2 = 0.001
112  real, pointer, dimension(:,:,:) :: T => null(), s => null()
113  real :: g_Rho0 ! G_Earth/Rho0 [L2 T-2 Z-1 R-1 ~> m4 s-2 kg-1].
114  real :: c2_scale ! A scaling factor for wave speeds to help control the growth of the determinant
115  ! and its derivative with lam between rows of the Thomas algorithm solver. The
116  ! exact value should not matter for the final result if it is an even power of 2.
117  real :: rescale, I_rescale
118  integer :: kf(SZI_(G))
119  integer, parameter :: max_itt = 10
120  real :: lam_it(max_itt), det_it(max_itt), ddet_it(max_itt)
121  logical :: use_EOS ! If true, density is calculated from T & S using an
122  ! equation of state.
123  integer :: kc
124  integer :: i, j, k, k2, itt, is, ie, js, je, nz
125  real :: hw, sum_hc
126  real :: gp ! A limited local copy of gprime [L2 Z-1 T-2 ~> m s-2]
127  real :: N2min ! A minimum buoyancy frequency [T-2 ~> s-2]
128  logical :: l_use_ebt_mode, calc_modal_structure
129  real :: l_mono_N2_column_fraction, l_mono_N2_depth
130  real :: mode_struct(SZK_(G)), ms_min, ms_max, ms_sq
131 
132  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
133 
134  if (.not. associated(cs)) call mom_error(fatal, "MOM_wave_speed: "// &
135  "Module must be initialized before it is used.")
136  if (present(full_halos)) then ; if (full_halos) then
137  is = g%isd ; ie = g%ied ; js = g%jsd ; je = g%jed
138  endif ; endif
139 
140  l2_to_z2 = us%L_to_Z**2
141 
142  l_use_ebt_mode = cs%use_ebt_mode
143  if (present(use_ebt_mode)) l_use_ebt_mode = use_ebt_mode
144  l_mono_n2_column_fraction = cs%mono_N2_column_fraction
145  if (present(mono_n2_column_fraction)) l_mono_n2_column_fraction = mono_n2_column_fraction
146  l_mono_n2_depth = cs%mono_N2_depth
147  if (present(mono_n2_depth)) l_mono_n2_depth = mono_n2_depth
148  calc_modal_structure = l_use_ebt_mode
149  if (present(modal_structure)) calc_modal_structure = .true.
150  if (calc_modal_structure) then
151  do k=1,nz; do j=js,je; do i=is,ie
152  modal_structure(i,j,k) = 0.0
153  enddo ; enddo ; enddo
154  endif
155 
156  s => tv%S ; t => tv%T
157  g_rho0 = gv%g_Earth / gv%Rho0
158  z_to_pa = gv%Z_to_H * gv%H_to_Pa
159  use_eos = associated(tv%eqn_of_state)
160 
161  rescale = 1024.0**4 ; i_rescale = 1.0/rescale
162  ! The following two lines give identical results:
163  ! c2_scale = 16.0 * US%m_s_to_L_T**2
164  c2_scale = us%m_s_to_L_T**2
165 
166  min_h_frac = tol1 / real(nz)
167 !$OMP parallel do default(none) shared(is,ie,js,je,nz,h,G,GV,US,min_h_frac,use_EOS,T,S,tv,&
168 !$OMP calc_modal_structure,l_use_ebt_mode,modal_structure, &
169 !$OMP l_mono_N2_column_fraction,l_mono_N2_depth,CS, &
170 !$OMP Z_to_Pa,cg1,g_Rho0,rescale,I_rescale,L2_to_Z2,c2_scale) &
171 !$OMP private(htot,hmin,kf,H_here,HxT_here,HxS_here,HxR_here, &
172 !$OMP Hf,Tf,Sf,Rf,pres,T_int,S_int,drho_dT, &
173 !$OMP drho_dS,drxh_sum,kc,Hc,Hc_H,Tc,Sc,I_Hnew,gprime,&
174 !$OMP Rc,speed2_tot,Igl,Igu,Igd,lam0,lam,lam_it,dlam, &
175 !$OMP mode_struct,sum_hc,N2min,gp,hw, &
176 !$OMP ms_min,ms_max,ms_sq, &
177 !$OMP det,ddet,detKm1,ddetKm1,detKm2,ddetKm2,det_it,ddet_it)
178  do j=js,je
179  ! First merge very thin layers with the one above (or below if they are
180  ! at the top). This also transposes the row order so that columns can
181  ! be worked upon one at a time.
182  do i=is,ie ; htot(i) = 0.0 ; enddo
183  do k=1,nz ; do i=is,ie ; htot(i) = htot(i) + h(i,j,k)*gv%H_to_Z ; enddo ; enddo
184 
185  do i=is,ie
186  hmin(i) = htot(i)*min_h_frac ; kf(i) = 1 ; h_here(i) = 0.0
187  hxt_here(i) = 0.0 ; hxs_here(i) = 0.0 ; hxr_here(i) = 0.0
188  enddo
189  if (use_eos) then
190  do k=1,nz ; do i=is,ie
191  if ((h_here(i) > hmin(i)) .and. (h(i,j,k)*gv%H_to_Z > hmin(i))) then
192  hf(kf(i),i) = h_here(i)
193  tf(kf(i),i) = hxt_here(i) / h_here(i)
194  sf(kf(i),i) = hxs_here(i) / h_here(i)
195  kf(i) = kf(i) + 1
196 
197  ! Start a new layer
198  h_here(i) = h(i,j,k)*gv%H_to_Z
199  hxt_here(i) = (h(i,j,k)*gv%H_to_Z)*t(i,j,k)
200  hxs_here(i) = (h(i,j,k)*gv%H_to_Z)*s(i,j,k)
201  else
202  h_here(i) = h_here(i) + h(i,j,k)*gv%H_to_Z
203  hxt_here(i) = hxt_here(i) + (h(i,j,k)*gv%H_to_Z)*t(i,j,k)
204  hxs_here(i) = hxs_here(i) + (h(i,j,k)*gv%H_to_Z)*s(i,j,k)
205  endif
206  enddo ; enddo
207  do i=is,ie ; if (h_here(i) > 0.0) then
208  hf(kf(i),i) = h_here(i)
209  tf(kf(i),i) = hxt_here(i) / h_here(i)
210  sf(kf(i),i) = hxs_here(i) / h_here(i)
211  endif ; enddo
212  else
213  do k=1,nz ; do i=is,ie
214  if ((h_here(i) > hmin(i)) .and. (h(i,j,k)*gv%H_to_Z > hmin(i))) then
215  hf(kf(i),i) = h_here(i) ; rf(kf(i),i) = hxr_here(i) / h_here(i)
216  kf(i) = kf(i) + 1
217 
218  ! Start a new layer
219  h_here(i) = h(i,j,k)*gv%H_to_Z
220  hxr_here(i) = (h(i,j,k)*gv%H_to_Z)*gv%Rlay(k)
221  else
222  h_here(i) = h_here(i) + h(i,j,k)*gv%H_to_Z
223  hxr_here(i) = hxr_here(i) + (h(i,j,k)*gv%H_to_Z)*gv%Rlay(k)
224  endif
225  enddo ; enddo
226  do i=is,ie ; if (h_here(i) > 0.0) then
227  hf(kf(i),i) = h_here(i) ; rf(kf(i),i) = hxr_here(i) / h_here(i)
228  endif ; enddo
229  endif
230 
231  ! From this point, we can work on individual columns without causing memory
232  ! to have page faults.
233  do i=is,ie ; if (g%mask2dT(i,j) > 0.5) then
234  if (use_eos) then
235  pres(1) = 0.0
236  do k=2,kf(i)
237  pres(k) = pres(k-1) + z_to_pa*hf(k-1,i)
238  t_int(k) = 0.5*(tf(k,i)+tf(k-1,i))
239  s_int(k) = 0.5*(sf(k,i)+sf(k-1,i))
240  enddo
241  call calculate_density_derivs(t_int, s_int, pres, drho_dt, drho_ds, 2, &
242  kf(i)-1, tv%eqn_of_state, scale=us%kg_m3_to_R)
243 
244  ! Sum the reduced gravities to find out how small a density difference
245  ! is negligibly small.
246  drxh_sum = 0.0
247  do k=2,kf(i)
248  drxh_sum = drxh_sum + 0.5*(hf(k-1,i)+hf(k,i)) * &
249  max(0.0,drho_dt(k)*(tf(k,i)-tf(k-1,i)) + &
250  drho_ds(k)*(sf(k,i)-sf(k-1,i)))
251  enddo
252  else
253  drxh_sum = 0.0
254  do k=2,kf(i)
255  drxh_sum = drxh_sum + 0.5*(hf(k-1,i)+hf(k,i)) * &
256  max(0.0,rf(k,i)-rf(k-1,i))
257  enddo
258  endif
259 
260  if (calc_modal_structure) then
261  mode_struct(:) = 0.
262  endif
263 
264  ! Find gprime across each internal interface, taking care of convective
265  ! instabilities by merging layers.
266  if (drxh_sum <= 0.0) then
267  cg1(i,j) = 0.0
268  else
269  ! Merge layers to eliminate convective instabilities or exceedingly
270  ! small reduced gravities.
271  if (use_eos) then
272  kc = 1
273  hc(1) = hf(1,i) ; tc(1) = tf(1,i) ; sc(1) = sf(1,i)
274  do k=2,kf(i)
275  if ((drho_dt(k)*(tf(k,i)-tc(kc)) + drho_ds(k)*(sf(k,i)-sc(kc))) * &
276  (hc(kc) + hf(k,i)) < 2.0 * tol2*drxh_sum) then
277  ! Merge this layer with the one above and backtrack.
278  i_hnew = 1.0 / (hc(kc) + hf(k,i))
279  tc(kc) = (hc(kc)*tc(kc) + hf(k,i)*tf(k,i)) * i_hnew
280  sc(kc) = (hc(kc)*sc(kc) + hf(k,i)*sf(k,i)) * i_hnew
281  hc(kc) = (hc(kc) + hf(k,i))
282  ! Backtrack to remove any convective instabilities above... Note
283  ! that the tolerance is a factor of two larger, to avoid limit how
284  ! far back we go.
285  do k2=kc,2,-1
286  if ((drho_dt(k2)*(tc(k2)-tc(k2-1)) + drho_ds(k2)*(sc(k2)-sc(k2-1))) * &
287  (hc(k2) + hc(k2-1)) < tol2*drxh_sum) then
288  ! Merge the two bottommost layers. At this point kc = k2.
289  i_hnew = 1.0 / (hc(kc) + hc(kc-1))
290  tc(kc-1) = (hc(kc)*tc(kc) + hc(kc-1)*tc(kc-1)) * i_hnew
291  sc(kc-1) = (hc(kc)*sc(kc) + hc(kc-1)*sc(kc-1)) * i_hnew
292  hc(kc-1) = (hc(kc) + hc(kc-1))
293  kc = kc - 1
294  else ; exit ; endif
295  enddo
296  else
297  ! Add a new layer to the column.
298  kc = kc + 1
299  drho_ds(kc) = drho_ds(k) ; drho_dt(kc) = drho_dt(k)
300  tc(kc) = tf(k,i) ; sc(kc) = sf(k,i) ; hc(kc) = hf(k,i)
301  endif
302  enddo
303  ! At this point there are kc layers and the gprimes should be positive.
304  do k=2,kc ! Revisit this if non-Boussinesq.
305  gprime(k) = g_rho0 * (drho_dt(k)*(tc(k)-tc(k-1)) + &
306  drho_ds(k)*(sc(k)-sc(k-1)))
307  enddo
308  else ! .not.use_EOS
309  ! Do the same with density directly...
310  kc = 1
311  hc(1) = hf(1,i) ; rc(1) = rf(1,i)
312  do k=2,kf(i)
313  if ((rf(k,i) - rc(kc)) * (hc(kc) + hf(k,i)) < 2.0*tol2*drxh_sum) then
314  ! Merge this layer with the one above and backtrack.
315  rc(kc) = (hc(kc)*rc(kc) + hf(k,i)*rf(k,i)) / (hc(kc) + hf(k,i))
316  hc(kc) = (hc(kc) + hf(k,i))
317  ! Backtrack to remove any convective instabilities above... Note
318  ! that the tolerance is a factor of two larger, to avoid limit how
319  ! far back we go.
320  do k2=kc,2,-1
321  if ((rc(k2)-rc(k2-1)) * (hc(k2)+hc(k2-1)) < tol2*drxh_sum) then
322  ! Merge the two bottommost layers. At this point kc = k2.
323  rc(kc-1) = (hc(kc)*rc(kc) + hc(kc-1)*rc(kc-1)) / (hc(kc) + hc(kc-1))
324  hc(kc-1) = (hc(kc) + hc(kc-1))
325  kc = kc - 1
326  else ; exit ; endif
327  enddo
328  else
329  ! Add a new layer to the column.
330  kc = kc + 1
331  rc(kc) = rf(k,i) ; hc(kc) = hf(k,i)
332  endif
333  enddo
334  ! At this point there are kc layers and the gprimes should be positive.
335  do k=2,kc ! Revisit this if non-Boussinesq.
336  gprime(k) = g_rho0 * (rc(k)-rc(k-1))
337  enddo
338  endif ! use_EOS
339 
340  ! Sum the contributions from all of the interfaces to give an over-estimate
341  ! of the first-mode wave speed. Also populate Igl and Igu which are the
342  ! non-leading diagonals of the tridiagonal matrix.
343  if (kc >= 2) then
344  speed2_tot = 0.0
345  if (l_use_ebt_mode) then
346  igu(1) = 0. ! Neumann condition for pressure modes
347  sum_hc = hc(1)
348  n2min = l2_to_z2*gprime(2)/hc(1)
349  do k=2,kc
350  hw = 0.5*(hc(k-1)+hc(k))
351  gp = gprime(k)
352  if (l_mono_n2_column_fraction>0. .or. l_mono_n2_depth>=0.) then
353  if ( ((g%bathyT(i,j)-sum_hc < l_mono_n2_column_fraction*g%bathyT(i,j)) .or. &
354  ((l_mono_n2_depth >= 0.) .and. (sum_hc > l_mono_n2_depth))) .and. &
355  (l2_to_z2*gp > n2min*hw) ) then
356  ! Filters out regions where N2 increases with depth but only in a lower fraction
357  ! of the water column or below a certain depth.
358  gp = us%Z_to_L**2 * (n2min*hw)
359  else
360  n2min = l2_to_z2 * gp/hw
361  endif
362  endif
363  igu(k) = 1.0/(gp*hc(k))
364  igl(k-1) = 1.0/(gp*hc(k-1))
365  speed2_tot = speed2_tot + gprime(k)*(hc(k-1)+hc(k))*0.707
366  sum_hc = sum_hc + hc(k)
367  enddo
368  !Igl(kc) = 0. ! Neumann condition for pressure modes
369  igl(kc) = 2.*igu(kc) ! Dirichlet condition for pressure modes
370  else ! .not. l_use_ebt_mode
371  do k=2,kc
372  igl(k) = 1.0/(gprime(k)*hc(k)) ; igu(k) = 1.0/(gprime(k)*hc(k-1))
373  speed2_tot = speed2_tot + gprime(k)*(hc(k-1)+hc(k))
374  enddo
375  endif
376 
377  if (calc_modal_structure) then
378  mode_struct(1:kc) = 1. ! Uniform flow, first guess
379  endif
380 
381  ! Overestimate the speed to start with.
382  if (calc_modal_structure) then
383  lam0 = 0.5 / speed2_tot ; lam = lam0
384  else
385  lam0 = 1.0 / speed2_tot ; lam = lam0
386  endif
387  ! Find the determinant and its derivative with lam.
388  do itt=1,max_itt
389  lam_it(itt) = lam
390  if (l_use_ebt_mode) then
391  ! This initialization of det,ddet imply Neumann boundary conditions so that first 3 rows
392  ! of the matrix are
393  ! / b(1)-lam igl(1) 0 0 0 ... \
394  ! | igu(2) b(2)-lam igl(2) 0 0 ... |
395  ! | 0 igu(3) b(3)-lam igl(3) 0 ... |
396  ! which is consistent if the eigenvalue problem is for horizontal velocity or pressure modes.
397  !detKm1 = c2_scale*(Igl(1)-lam) ; ddetKm1 = -1.0*c2_scale
398  !det = (Igu(2)+Igl(2)-lam)*detKm1 - (Igu(2)*Igl(1)) ; ddet = (Igu(2)+Igl(2)-lam)*ddetKm1 - detKm1
399  detkm1 = 1.0 ; ddetkm1 = 0.0
400  det = (igl(1)-lam) ; ddet = -1.0
401  if (kc>1) then
402  ! Shift variables and rescale rows to avoid over- or underflow.
403  detkm2 = c2_scale*detkm1 ; ddetkm2 = c2_scale*ddetkm1
404  detkm1 = c2_scale*det ; ddetkm1 = c2_scale*ddet
405  det = (igu(2)+igl(2)-lam)*detkm1 - (igu(2)*igl(1))*detkm2
406  ddet = (igu(2)+igl(2)-lam)*ddetkm1 - (igu(2)*igl(1))*ddetkm2 - detkm1
407  endif
408  ! The last two rows of the pressure equation matrix are
409  ! | ... 0 igu(kc-1) b(kc-1)-lam igl(kc-1) |
410  ! \ ... 0 0 igu(kc) b(kc)-lam /
411  else
412  ! This initialization of det,ddet imply Dirichlet boundary conditions so that first 3 rows
413  ! of the matrix are
414  ! / b(2)-lam igl(2) 0 0 0 ... |
415  ! | igu(3) b(3)-lam igl(3) 0 0 ... |
416  ! | 0 igu43) b(4)-lam igl(4) 0 ... |
417  ! which is consistent if the eigenvalue problem is for vertical velocity modes.
418  detkm1 = 1.0 ; ddetkm1 = 0.0
419  det = (igu(2) + igl(2) - lam) ; ddet = -1.0
420  ! The last three rows of the w equation matrix are
421  ! | ... 0 igu(kc-1) b(kc-1)-lam igl(kc-1) 0 |
422  ! | ... 0 0 igu(kc-1) b(kc-1)-lam igl(kc-1) |
423  ! \ ... 0 0 0 igu(kc) b(kc)-lam /
424  endif
425  do k=3,kc
426  ! Shift variables and rescale rows to avoid over- or underflow.
427  detkm2 = c2_scale*detkm1 ; ddetkm2 = c2_scale*ddetkm1
428  detkm1 = c2_scale*det ; ddetkm1 = c2_scale*ddet
429 
430  det = (igu(k)+igl(k)-lam)*detkm1 - (igu(k)*igl(k-1))*detkm2
431  ddet = (igu(k)+igl(k)-lam)*ddetkm1 - (igu(k)*igl(k-1))*ddetkm2 - detkm1
432 
433  ! Rescale det & ddet if det is getting too large or too small.
434  if (abs(det) > rescale) then
435  det = i_rescale*det ; detkm1 = i_rescale*detkm1
436  ddet = i_rescale*ddet ; ddetkm1 = i_rescale*ddetkm1
437  elseif (abs(det) < i_rescale) then
438  det = rescale*det ; detkm1 = rescale*detkm1
439  ddet = rescale*ddet ; ddetkm1 = rescale*ddetkm1
440  endif
441  enddo
442  ! Use Newton's method iteration to find a new estimate of lam.
443  det_it(itt) = det ; ddet_it(itt) = ddet
444 
445  if ((ddet >= 0.0) .or. (-det > -0.5*lam*ddet)) then
446  ! lam was not an under-estimate, as intended, so Newton's method
447  ! may not be reliable; lam must be reduced, but not by more
448  ! than half.
449  lam = 0.5 * lam
450  dlam = -lam
451  else ! Newton's method is OK.
452  dlam = - det / ddet
453  lam = lam + dlam
454  endif
455 
456  if (calc_modal_structure) then
457  do k = 1,kc
458  igd(k) = igu(k) + igl(k)
459  enddo
460  call tdma6(kc, -igu, igd, -igl, lam, mode_struct)
461  ms_min = mode_struct(1)
462  ms_max = mode_struct(1)
463  ms_sq = mode_struct(1)**2
464  do k = 2,kc
465  ms_min = min(ms_min, mode_struct(k))
466  ms_max = max(ms_max, mode_struct(k))
467  ms_sq = ms_sq + mode_struct(k)**2
468  enddo
469  if (ms_min<0. .and. ms_max>0.) then ! Any zero crossings => lam is too high
470  lam = 0.5 * ( lam - dlam )
471  dlam = -lam
472  mode_struct(1:kc) = abs(mode_struct(1:kc)) / sqrt( ms_sq )
473  else
474  mode_struct(1:kc) = mode_struct(1:kc) / sqrt( ms_sq )
475  endif
476  endif
477 
478  if (abs(dlam) < tol2*lam) exit
479  enddo
480 
481  cg1(i,j) = 0.0
482  if (lam > 0.0) cg1(i,j) = 1.0 / sqrt(lam)
483 
484  if (present(modal_structure)) then
485  if (mode_struct(1)/=0.) then ! Normalize
486  mode_struct(1:kc) = mode_struct(1:kc) / mode_struct(1)
487  else
488  mode_struct(1:kc)=0.
489  endif
490  ! Note that remapping_core_h requires that the same units be used
491  ! for both the source and target grid thicknesses, here [H ~> m or kg m-2].
492  do k = 1,kc
493  hc_h(k) = gv%Z_to_H * hc(k)
494  enddo
495  call remapping_core_h(cs%remapping_CS, kc, hc_h(:), mode_struct, &
496  nz, h(i,j,:), modal_structure(i,j,:), &
497  1.0e-30*gv%m_to_H, 1.0e-10*gv%m_to_H)
498  endif
499  else
500  cg1(i,j) = 0.0
501  if (present(modal_structure)) modal_structure(i,j,:) = 0.
502  endif
503  endif ! cg1 /= 0.0
504  else
505  cg1(i,j) = 0.0 ! This is a land point.
506  if (present(modal_structure)) modal_structure(i,j,:) = 0.
507  endif ; enddo ! i-loop
508  enddo ! j-loop
509 

References mom_error_handler::mom_error(), mom_remapping::remapping_core_h(), and tdma6().

Referenced by mom_lateral_mixing_coeffs::calc_resoln_function(), and mom_diagnostics::calculate_diagnostic_fields().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ wave_speed_init()

subroutine, public mom_wave_speed::wave_speed_init ( type(wave_speed_cs), pointer  CS,
logical, intent(in), optional  use_ebt_mode,
real, intent(in), optional  mono_N2_column_fraction,
real, intent(in), optional  mono_N2_depth 
)

Initialize control structure for MOM_wave_speed.

Parameters
csControl structure for MOM_wave_speed
[in]use_ebt_modeIf true, use the equivalent barotropic mode instead of the first baroclinic mode.
[in]mono_n2_column_fractionThe lower fraction of water column over which N2 is limited as monotonic for the purposes of calculating the vertical modal structure.
[in]mono_n2_depthThe depth below which N2 is limited as monotonic for the purposes of calculating the vertical modal structure [Z ~> m].

Definition at line 1061 of file MOM_wave_speed.F90.

1061  type(wave_speed_CS), pointer :: CS !< Control structure for MOM_wave_speed
1062  logical, optional, intent(in) :: use_ebt_mode !< If true, use the equivalent
1063  !! barotropic mode instead of the first baroclinic mode.
1064  real, optional, intent(in) :: mono_N2_column_fraction !< The lower fraction of water column over
1065  !! which N2 is limited as monotonic for the purposes of
1066  !! calculating the vertical modal structure.
1067  real, optional, intent(in) :: mono_N2_depth !< The depth below which N2 is limited
1068  !! as monotonic for the purposes of calculating the
1069  !! vertical modal structure [Z ~> m].
1070 ! This include declares and sets the variable "version".
1071 #include "version_variable.h"
1072  character(len=40) :: mdl = "MOM_wave_speed" ! This module's name.
1073 
1074  if (associated(cs)) then
1075  call mom_error(warning, "wave_speed_init called with an "// &
1076  "associated control structure.")
1077  return
1078  else ; allocate(cs) ; endif
1079 
1080  ! Write all relevant parameters to the model log.
1081  call log_version(mdl, version)
1082 
1083  call wave_speed_set_param(cs, use_ebt_mode=use_ebt_mode, mono_n2_column_fraction=mono_n2_column_fraction)
1084 
1085  call initialize_remapping(cs%remapping_CS, 'PLM', boundary_extrapolation=.false.)
1086 

References mom_remapping::initialize_remapping(), mom_error_handler::mom_error(), and wave_speed_set_param().

Referenced by mom_diagnostics::mom_diagnostics_init(), and mom_lateral_mixing_coeffs::varmix_init().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ wave_speed_set_param()

subroutine, public mom_wave_speed::wave_speed_set_param ( type(wave_speed_cs), pointer  CS,
logical, intent(in), optional  use_ebt_mode,
real, intent(in), optional  mono_N2_column_fraction,
real, intent(in), optional  mono_N2_depth 
)

Sets internal parameters for MOM_wave_speed.

Parameters
csControl structure for MOM_wave_speed
[in]use_ebt_modeIf true, use the equivalent barotropic mode instead of the first baroclinic mode.
[in]mono_n2_column_fractionThe lower fraction of water column over which N2 is limited as monotonic for the purposes of calculating the vertical modal structure.
[in]mono_n2_depthThe depth below which N2 is limited as monotonic for the purposes of calculating the vertical modal structure [Z ~> m].

Definition at line 1091 of file MOM_wave_speed.F90.

1091  type(wave_speed_CS), pointer :: CS !< Control structure for MOM_wave_speed
1092  logical, optional, intent(in) :: use_ebt_mode !< If true, use the equivalent
1093  !! barotropic mode instead of the first baroclinic mode.
1094  real, optional, intent(in) :: mono_N2_column_fraction !< The lower fraction of water column over
1095  !! which N2 is limited as monotonic for the purposes of
1096  !! calculating the vertical modal structure.
1097  real, optional, intent(in) :: mono_N2_depth !< The depth below which N2 is limited
1098  !! as monotonic for the purposes of calculating the
1099  !! vertical modal structure [Z ~> m].
1100 
1101  if (.not.associated(cs)) call mom_error(fatal, &
1102  "wave_speed_set_param called with an associated control structure.")
1103 
1104  if (present(use_ebt_mode)) cs%use_ebt_mode = use_ebt_mode
1105  if (present(mono_n2_column_fraction)) cs%mono_N2_column_fraction = mono_n2_column_fraction
1106  if (present(mono_n2_depth)) cs%mono_N2_depth = mono_n2_depth
1107 

References mom_error_handler::mom_error().

Referenced by wave_speed_init().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ wave_speeds()

subroutine, public mom_wave_speed::wave_speeds ( real, dimension(szi_(g),szj_(g),szk_(g)), 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,
integer, intent(in)  nmodes,
real, dimension(g%isd:g%ied,g%jsd:g%jed,nmodes), intent(out)  cn,
type(wave_speed_cs), optional, pointer  CS,
logical, intent(in), optional  full_halos 
)

Calculates the wave speeds for the first few barolinic modes.

Parameters
[in]gOcean grid structure
[in]gvVertical grid structure
[in]usA dimensional unit scaling type
[in]hLayer thickness [H ~> m or kg m-2]
[in]tvThermodynamic variables
[in]nmodesNumber of modes
[out]cnWaves speeds [L T-1 ~> m s-1]
csControl structure for MOM_wave_speed
[in]full_halosIf true, do the calculation over the entire computational domain.

Definition at line 561 of file MOM_wave_speed.F90.

561  type(ocean_grid_type), intent(in) :: G !< Ocean grid structure
562  type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure
563  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
564  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h !< Layer thickness [H ~> m or kg m-2]
565  type(thermo_var_ptrs), intent(in) :: tv !< Thermodynamic variables
566  integer, intent(in) :: nmodes !< Number of modes
567  real, dimension(G%isd:G%ied,G%jsd:G%jed,nmodes), intent(out) :: cn !< Waves speeds [L T-1 ~> m s-1]
568  type(wave_speed_CS), optional, pointer :: CS !< Control structure for MOM_wave_speed
569  logical, optional, intent(in) :: full_halos !< If true, do the calculation
570  !! over the entire computational domain.
571  ! Local variables
572  real, dimension(SZK_(G)+1) :: &
573  dRho_dT, & ! Partial derivative of density with temperature [R degC-1 ~> kg m-3 degC-1]
574  dRho_dS, & ! Partial derivative of density with salinity [R ppt-1 ~> kg m-3 ppt-1]
575  pres, & ! Interface pressure [Pa]
576  T_int, & ! Temperature interpolated to interfaces [degC]
577  S_int, & ! Salinity interpolated to interfaces [ppt]
578  gprime ! The reduced gravity across each interface [L2 Z-1 T-2 ~> m s-2].
579  real, dimension(SZK_(G)) :: &
580  Igl, Igu ! The inverse of the reduced gravity across an interface times
581  ! the thickness of the layer below (Igl) or above (Igu) it, in [T2 L-2 ~> s2 m-2].
582  real, dimension(SZK_(G)-1) :: &
583  a_diag, b_diag, c_diag
584  ! diagonals of tridiagonal matrix; one value for each
585  ! interface (excluding surface and bottom) [T2 L-2 ~> s2 m-2]
586  real, dimension(SZK_(G),SZI_(G)) :: &
587  Hf, & ! Layer thicknesses after very thin layers are combined [Z ~> m]
588  Tf, & ! Layer temperatures after very thin layers are combined [degC]
589  Sf, & ! Layer salinities after very thin layers are combined [ppt]
590  Rf ! Layer densities after very thin layers are combined [R ~> kg m-3]
591  real, dimension(SZK_(G)) :: &
592  Hc, & ! A column of layer thicknesses after convective istabilities are removed [Z ~> m]
593  Tc, & ! A column of layer temperatures after convective istabilities are removed [degC]
594  Sc, & ! A column of layer salinites after convective istabilities are removed [ppt]
595  Rc ! A column of layer densities after convective istabilities are removed [R ~> kg m-3]
596  real :: c1_thresh ! if c1 is below this value, don't bother calculating
597  ! cn values for higher modes [L T-1 ~> m s-1]
598  real :: det, ddet ! determinant & its derivative of eigen system
599  real :: lam_1 ! approximate mode-1 eigenvalue [T2 L-2 ~> s2 m-2]
600  real :: lam_n ! approximate mode-n eigenvalue [T2 L-2 ~> s2 m-2]
601  real :: dlam ! increment in lam for Newton's method [T2 L-2 ~> s2 m-2]
602  real :: lamMin ! minimum lam value for root searching range [T2 L-2 ~> s2 m-2]
603  real :: lamMax ! maximum lam value for root searching range [T2 L-2 ~> s2 m-2]
604  real :: lamInc ! width of moving window for root searching [T2 L-2 ~> s2 m-2]
605  real :: det_l,det_r ! determinant value at left and right of window
606  real :: ddet_l,ddet_r ! derivative of determinant at left and right of window
607  real :: det_sub,ddet_sub! derivative of determinant at subinterval endpoint
608  real :: xl,xr ! lam guesses at left and right of window [T2 L-2 ~> s2 m-2]
609  real :: xl_sub ! lam guess at left of subinterval window [T2 L-2 ~> s2 m-2]
610  real,dimension(nmodes) :: &
611  xbl,xbr ! lam guesses bracketing a zero-crossing (root) [T2 L-2 ~> s2 m-2]
612  integer :: numint ! number of widows (intervals) in root searching range
613  integer :: nrootsfound ! number of extra roots found (not including 1st root)
614  real :: min_h_frac
615  real :: Z_to_Pa ! A conversion factor from thicknesses (in Z) to pressure (in Pa)
616  real, dimension(SZI_(G)) :: &
617  htot, hmin, & ! Thicknesses [Z ~> m]
618  H_here, & ! A thickness [Z ~> m]
619  HxT_here, & ! A layer integrated temperature [degC Z ~> degC m]
620  HxS_here, & ! A layer integrated salinity [ppt Z ~> ppt m]
621  HxR_here ! A layer integrated density [R Z ~> kg m-2]
622  real :: speed2_tot ! overestimate of the mode-1 speed squared [L2 T-2 ~> m2 s-2]
623  real :: speed2_min ! minimum mode speed (squared) to consider in root searching [L2 T-2 ~> m2 s-2]
624  real, parameter :: reduct_factor = 0.5
625  ! factor used in setting speed2_min [nondim]
626  real :: I_Hnew ! The inverse of a new layer thickness [Z-1 ~> m-1]
627  real :: drxh_sum ! The sum of density diffrences across interfaces times thicknesses [R Z ~> kg m-2]
628  real, parameter :: tol1 = 0.0001, tol2 = 0.001
629  real, pointer, dimension(:,:,:) :: T => null(), s => null()
630  real :: g_Rho0 ! G_Earth/Rho0 [L2 T-2 Z-1 R-1 ~> m4 s-2 kg-1].
631  integer :: kf(SZI_(G))
632  integer, parameter :: max_itt = 10
633  logical :: use_EOS ! If true, density is calculated from T & S using the equation of state.
634  real, dimension(SZK_(G)+1) :: z_int
635  ! real, dimension(SZK_(G)+1) :: N2 ! The local squared buoyancy frequency [T-2 ~> s-2]
636  integer :: nsub ! number of subintervals used for root finding
637  integer, parameter :: sub_it_max = 4
638  ! maximum number of times to subdivide interval
639  ! for root finding (# intervals = 2**sub_it_max)
640  logical :: sub_rootfound ! if true, subdivision has located root
641  integer :: kc, nrows
642  integer :: sub, sub_it
643  integer :: i, j, k, k2, itt, is, ie, js, je, nz, row, iint, m, ig, jg
644 
645  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
646 
647  if (present(cs)) then
648  if (.not. associated(cs)) call mom_error(fatal, "MOM_wave_speed: "// &
649  "Module must be initialized before it is used.")
650  endif
651 
652  if (present(full_halos)) then ; if (full_halos) then
653  is = g%isd ; ie = g%ied ; js = g%jsd ; je = g%jed
654  endif ; endif
655 
656  s => tv%S ; t => tv%T
657  g_rho0 = gv%g_Earth / gv%Rho0
658  use_eos = associated(tv%eqn_of_state)
659  z_to_pa = gv%Z_to_H * gv%H_to_Pa
660  c1_thresh = 0.01*us%m_s_to_L_T
661 
662  min_h_frac = tol1 / real(nz)
663  !$OMP parallel do default(private) shared(is,ie,js,je,nz,h,G,GV,US,min_h_frac,use_EOS,T,S, &
664  !$OMP Z_to_Pa,tv,cn,g_Rho0,nmodes)
665  do j=js,je
666  ! First merge very thin layers with the one above (or below if they are
667  ! at the top). This also transposes the row order so that columns can
668  ! be worked upon one at a time.
669  do i=is,ie ; htot(i) = 0.0 ; enddo
670  do k=1,nz ; do i=is,ie ; htot(i) = htot(i) + h(i,j,k)*gv%H_to_Z ; enddo ; enddo
671 
672  do i=is,ie
673  hmin(i) = htot(i)*min_h_frac ; kf(i) = 1 ; h_here(i) = 0.0
674  hxt_here(i) = 0.0 ; hxs_here(i) = 0.0 ; hxr_here(i) = 0.0
675  enddo
676  if (use_eos) then
677  do k=1,nz ; do i=is,ie
678  if ((h_here(i) > hmin(i)) .and. (h(i,j,k)*gv%H_to_Z > hmin(i))) then
679  hf(kf(i),i) = h_here(i)
680  tf(kf(i),i) = hxt_here(i) / h_here(i)
681  sf(kf(i),i) = hxs_here(i) / h_here(i)
682  kf(i) = kf(i) + 1
683 
684  ! Start a new layer
685  h_here(i) = h(i,j,k)*gv%H_to_Z
686  hxt_here(i) = (h(i,j,k)*gv%H_to_Z)*t(i,j,k)
687  hxs_here(i) = (h(i,j,k)*gv%H_to_Z)*s(i,j,k)
688  else
689  h_here(i) = h_here(i) + h(i,j,k)*gv%H_to_Z
690  hxt_here(i) = hxt_here(i) + (h(i,j,k)*gv%H_to_Z)*t(i,j,k)
691  hxs_here(i) = hxs_here(i) + (h(i,j,k)*gv%H_to_Z)*s(i,j,k)
692  endif
693  enddo ; enddo
694  do i=is,ie ; if (h_here(i) > 0.0) then
695  hf(kf(i),i) = h_here(i)
696  tf(kf(i),i) = hxt_here(i) / h_here(i)
697  sf(kf(i),i) = hxs_here(i) / h_here(i)
698  endif ; enddo
699  else
700  do k=1,nz ; do i=is,ie
701  if ((h_here(i) > hmin(i)) .and. (h(i,j,k)*gv%H_to_Z > hmin(i))) then
702  hf(kf(i),i) = h_here(i) ; rf(kf(i),i) = hxr_here(i) / h_here(i)
703  kf(i) = kf(i) + 1
704 
705  ! Start a new layer
706  h_here(i) = h(i,j,k)*gv%H_to_Z
707  hxr_here(i) = (h(i,j,k)*gv%H_to_Z)*gv%Rlay(k)
708  else
709  h_here(i) = h_here(i) + h(i,j,k)*gv%H_to_Z
710  hxr_here(i) = hxr_here(i) + (h(i,j,k)*gv%H_to_Z)*gv%Rlay(k)
711  endif
712  enddo ; enddo
713  do i=is,ie ; if (h_here(i) > 0.0) then
714  hf(kf(i),i) = h_here(i) ; rf(kf(i),i) = hxr_here(i) / h_here(i)
715  endif ; enddo
716  endif
717 
718  ! From this point, we can work on individual columns without causing memory
719  ! to have page faults.
720  do i=is,ie
721  if (g%mask2dT(i,j) > 0.5) then
722  if (use_eos) then
723  pres(1) = 0.0
724  do k=2,kf(i)
725  pres(k) = pres(k-1) + z_to_pa*hf(k-1,i)
726  t_int(k) = 0.5*(tf(k,i)+tf(k-1,i))
727  s_int(k) = 0.5*(sf(k,i)+sf(k-1,i))
728  enddo
729  call calculate_density_derivs(t_int, s_int, pres, drho_dt, drho_ds, 2, &
730  kf(i)-1, tv%eqn_of_state, scale=us%kg_m3_to_R)
731 
732  ! Sum the reduced gravities to find out how small a density difference
733  ! is negligibly small.
734  drxh_sum = 0.0
735  do k=2,kf(i)
736  drxh_sum = drxh_sum + 0.5*(hf(k-1,i)+hf(k,i)) * &
737  max(0.0,drho_dt(k)*(tf(k,i)-tf(k-1,i)) + &
738  drho_ds(k)*(sf(k,i)-sf(k-1,i)))
739  enddo
740  else
741  drxh_sum = 0.0
742  do k=2,kf(i)
743  drxh_sum = drxh_sum + 0.5*(hf(k-1,i)+hf(k,i)) * &
744  max(0.0,rf(k,i)-rf(k-1,i))
745  enddo
746  endif
747  ! Find gprime across each internal interface, taking care of convective
748  ! instabilities by merging layers.
749  if (drxh_sum <= 0.0) then
750  cn(i,j,:) = 0.0
751  else
752  ! Merge layers to eliminate convective instabilities or exceedingly
753  ! small reduced gravities.
754  if (use_eos) then
755  kc = 1
756  hc(1) = hf(1,i) ; tc(1) = tf(1,i) ; sc(1) = sf(1,i)
757  do k=2,kf(i)
758  if ((drho_dt(k)*(tf(k,i)-tc(kc)) + drho_ds(k)*(sf(k,i)-sc(kc))) * &
759  (hc(kc) + hf(k,i)) < 2.0 * tol2*drxh_sum) then
760  ! Merge this layer with the one above and backtrack.
761  i_hnew = 1.0 / (hc(kc) + hf(k,i))
762  tc(kc) = (hc(kc)*tc(kc) + hf(k,i)*tf(k,i)) * i_hnew
763  sc(kc) = (hc(kc)*sc(kc) + hf(k,i)*sf(k,i)) * i_hnew
764  hc(kc) = (hc(kc) + hf(k,i))
765  ! Backtrack to remove any convective instabilities above... Note
766  ! that the tolerance is a factor of two larger, to avoid limit how
767  ! far back we go.
768  do k2=kc,2,-1
769  if ((drho_dt(k2)*(tc(k2)-tc(k2-1)) + drho_ds(k2)*(sc(k2)-sc(k2-1))) * &
770  (hc(k2) + hc(k2-1)) < tol2*drxh_sum) then
771  ! Merge the two bottommost layers. At this point kc = k2.
772  i_hnew = 1.0 / (hc(kc) + hc(kc-1))
773  tc(kc-1) = (hc(kc)*tc(kc) + hc(kc-1)*tc(kc-1)) * i_hnew
774  sc(kc-1) = (hc(kc)*sc(kc) + hc(kc-1)*sc(kc-1)) * i_hnew
775  hc(kc-1) = (hc(kc) + hc(kc-1))
776  kc = kc - 1
777  else ; exit ; endif
778  enddo
779  else
780  ! Add a new layer to the column.
781  kc = kc + 1
782  drho_ds(kc) = drho_ds(k) ; drho_dt(kc) = drho_dt(k)
783  tc(kc) = tf(k,i) ; sc(kc) = sf(k,i) ; hc(kc) = hf(k,i)
784  endif
785  enddo
786  ! At this point there are kc layers and the gprimes should be positive.
787  do k=2,kc ! Revisit this if non-Boussinesq.
788  gprime(k) = g_rho0 * (drho_dt(k)*(tc(k)-tc(k-1)) + &
789  drho_ds(k)*(sc(k)-sc(k-1)))
790  enddo
791  else ! .not.use_EOS
792  ! Do the same with density directly...
793  kc = 1
794  hc(1) = hf(1,i) ; rc(1) = rf(1,i)
795  do k=2,kf(i)
796  if ((rf(k,i) - rc(kc)) * (hc(kc) + hf(k,i)) < 2.0*tol2*drxh_sum) then
797  ! Merge this layer with the one above and backtrack.
798  rc(kc) = (hc(kc)*rc(kc) + hf(k,i)*rf(k,i)) / (hc(kc) + hf(k,i))
799  hc(kc) = (hc(kc) + hf(k,i))
800  ! Backtrack to remove any convective instabilities above... Note
801  ! that the tolerance is a factor of two larger, to avoid limit how
802  ! far back we go.
803  do k2=kc,2,-1
804  if ((rc(k2)-rc(k2-1)) * (hc(k2)+hc(k2-1)) < tol2*drxh_sum) then
805  ! Merge the two bottommost layers. At this point kc = k2.
806  rc(kc-1) = (hc(kc)*rc(kc) + hc(kc-1)*rc(kc-1)) / (hc(kc) + hc(kc-1))
807  hc(kc-1) = (hc(kc) + hc(kc-1))
808  kc = kc - 1
809  else ; exit ; endif
810  enddo
811  else
812  ! Add a new layer to the column.
813  kc = kc + 1
814  rc(kc) = rf(k,i) ; hc(kc) = hf(k,i)
815  endif
816  enddo
817  ! At this point there are kc layers and the gprimes should be positive.
818  do k=2,kc ! Revisit this if non-Boussinesq.
819  gprime(k) = g_rho0 * (rc(k)-rc(k-1))
820  enddo
821  endif ! use_EOS
822 
823  !-----------------NOW FIND WAVE SPEEDS---------------------------------------
824  ig = i + g%idg_offset ; jg = j + g%jdg_offset
825  ! Sum the contributions from all of the interfaces to give an over-estimate
826  ! of the first-mode wave speed.
827  if (kc >= 2) then
828  ! Set depth at surface
829  z_int(1) = 0.0
830  ! initialize speed2_tot
831  speed2_tot = 0.0
832  ! Calculate Igu, Igl, depth, and N2 at each interior interface
833  ! [excludes surface (K=1) and bottom (K=kc+1)]
834  do k=2,kc
835  igl(k) = 1.0/(gprime(k)*hc(k)) ; igu(k) = 1.0/(gprime(k)*hc(k-1))
836  z_int(k) = z_int(k-1) + hc(k-1)
837  ! N2(K) = US%L_to_Z**2*gprime(K)/(0.5*(Hc(k)+Hc(k-1)))
838  speed2_tot = speed2_tot + gprime(k)*(hc(k-1)+hc(k))
839  enddo
840  ! Set stratification for surface and bottom (setting equal to nearest interface for now)
841  ! N2(1) = N2(2) ; N2(kc+1) = N2(kc)
842  ! Calculate depth at bottom
843  z_int(kc+1) = z_int(kc)+hc(kc)
844  ! check that thicknesses sum to total depth
845  if (abs(z_int(kc+1)-htot(i)) > 1.e-12*htot(i)) then
846  call mom_error(fatal, "wave_structure: mismatch in total depths")
847  endif
848 
849  ! Define the diagonals of the tridiagonal matrix
850  ! First, populate interior rows
851  do k=3,kc-1
852  row = k-1 ! indexing for TD matrix rows
853  a_diag(row) = -igu(k)
854  b_diag(row) = igu(k)+igl(k)
855  c_diag(row) = -igl(k)
856  enddo
857  ! Populate top row of tridiagonal matrix
858  k=2 ; row = k-1
859  a_diag(row) = 0.0
860  b_diag(row) = igu(k)+igl(k)
861  c_diag(row) = -igl(k)
862  ! Populate bottom row of tridiagonal matrix
863  k=kc ; row = k-1
864  a_diag(row) = -igu(k)
865  b_diag(row) = igu(k)+igl(k)
866  c_diag(row) = 0.0
867  ! Total number of rows in the matrix = number of interior interfaces
868  nrows = kc-1
869 
870  ! Under estimate the first eigenvalue to start with.
871  lam_1 = 1.0 / speed2_tot
872 
873  ! Find the first eigen value
874  do itt=1,max_itt
875  ! calculate the determinant of (A-lam_1*I)
876  call tridiag_det(a_diag(1:nrows),b_diag(1:nrows),c_diag(1:nrows), &
877  nrows,lam_1,det,ddet, row_scale=us%m_s_to_L_T**2)
878  ! Use Newton's method iteration to find a new estimate of lam_1
879  !det = det_it(itt) ; ddet = ddet_it(itt)
880  if ((ddet >= 0.0) .or. (-det > -0.5*lam_1*ddet)) then
881  ! lam_1 was not an under-estimate, as intended, so Newton's method
882  ! may not be reliable; lam_1 must be reduced, but not by more
883  ! than half.
884  lam_1 = 0.5 * lam_1
885  else ! Newton's method is OK.
886  dlam = - det / ddet
887  lam_1 = lam_1 + dlam
888  if (abs(dlam) < tol2*lam_1) then
889  ! calculate 1st mode speed
890  if (lam_1 > 0.0) cn(i,j,1) = 1.0 / sqrt(lam_1)
891  exit
892  endif
893  endif
894  enddo
895 
896  ! Find other eigen values if c1 is of significant magnitude, > cn_thresh
897  nrootsfound = 0 ! number of extra roots found (not including 1st root)
898  if (nmodes>1 .and. kc>=nmodes+1 .and. cn(i,j,1)>c1_thresh) then
899  ! Set the the range to look for the other desired eigen values
900  ! set min value just greater than the 1st root (found above)
901  lammin = lam_1*(1.0 + tol2)
902  ! set max value based on a low guess at wavespeed for highest mode
903  speed2_min = (reduct_factor*cn(i,j,1)/real(nmodes))**2
904  lammax = 1.0 / speed2_min
905  ! set width of interval (not sure about this - BDM)
906  laminc = 0.5*lam_1
907  ! set number of intervals within search range
908  numint = nint((lammax - lammin)/laminc)
909 
910  ! Find intervals containing zero-crossings (roots) of the determinant
911  ! that are beyond the first root
912 
913  ! find det_l of first interval (det at left endpoint)
914  call tridiag_det(a_diag(1:nrows),b_diag(1:nrows),c_diag(1:nrows), &
915  nrows,lammin,det_l,ddet_l, row_scale=us%m_s_to_L_T**2)
916  ! move interval window looking for zero-crossings************************
917  do iint=1,numint
918  xr = lammin + laminc * iint
919  xl = xr - laminc
920  call tridiag_det(a_diag(1:nrows),b_diag(1:nrows),c_diag(1:nrows), &
921  nrows,xr,det_r,ddet_r, row_scale=us%m_s_to_L_T**2)
922  if (det_l*det_r < 0.0) then ! if function changes sign
923  if (det_l*ddet_l < 0.0) then ! if function at left is headed to zero
924  nrootsfound = nrootsfound + 1
925  xbl(nrootsfound) = xl
926  xbr(nrootsfound) = xr
927  else
928  ! function changes sign but has a local max/min in interval,
929  ! try subdividing interval as many times as necessary (or sub_it_max).
930  ! loop that increases number of subintervals:
931  !call MOM_error(WARNING, "determinant changes sign"// &
932  ! "but has a local max/min in interval;"//&
933  ! " reduce increment in lam.")
934  ! begin subdivision loop -------------------------------------------
935  sub_rootfound = .false. ! initialize
936  do sub_it=1,sub_it_max
937  nsub = 2**sub_it ! number of subintervals; nsub=2,4,8,...
938  ! loop over each subinterval:
939  do sub=1,nsub-1,2 ! only check odds; sub = 1; 1,3; 1,3,5,7;...
940  xl_sub = xl + laminc/(nsub)*sub
941  call tridiag_det(a_diag(1:nrows),b_diag(1:nrows),c_diag(1:nrows), &
942  nrows,xl_sub,det_sub,ddet_sub, row_scale=us%m_s_to_L_T**2)
943  if (det_sub*det_r < 0.0) then ! if function changes sign
944  if (det_sub*ddet_sub < 0.0) then ! if function at left is headed to zero
945  sub_rootfound = .true.
946  nrootsfound = nrootsfound + 1
947  xbl(nrootsfound) = xl_sub
948  xbr(nrootsfound) = xr
949  exit ! exit sub loop
950  endif ! headed toward zero
951  endif ! sign change
952  enddo ! sub-loop
953  if (sub_rootfound) exit ! root has been found, exit sub_it loop
954  ! Otherwise, function changes sign but has a local max/min in one of the
955  ! sub intervals, try subdividing again unless sub_it_max has been reached.
956  if (sub_it == sub_it_max) then
957  call mom_error(warning, "wave_speed: root not found "// &
958  " after sub_it_max subdivisions of original"// &
959  " interval.")
960  endif ! sub_it == sub_it_max
961  enddo ! sub_it-loop-------------------------------------------------
962  endif ! det_l*ddet_l < 0.0
963  endif ! det_l*det_r < 0.0
964  ! exit iint-loop if all desired roots have been found
965  if (nrootsfound >= nmodes-1) then
966  ! exit if all additional roots found
967  exit
968  elseif (iint == numint) then
969  ! oops, lamMax not large enough - could add code to increase (BDM)
970  ! set unfound modes to zero for now (BDM)
971  cn(i,j,nrootsfound+2:nmodes) = 0.0
972  else
973  ! else shift interval and keep looking until nmodes or numint is reached
974  det_l = det_r
975  ddet_l = ddet_r
976  endif
977  enddo ! iint-loop
978 
979  ! Use Newton's method to find the roots within the identified windows
980  do m=1,nrootsfound ! loop over the root-containing widows (excluding 1st mode)
981  lam_n = xbl(m) ! first guess is left edge of window
982  do itt=1,max_itt
983  ! calculate the determinant of (A-lam_n*I)
984  call tridiag_det(a_diag(1:nrows),b_diag(1:nrows),c_diag(1:nrows), &
985  nrows,lam_n,det,ddet, row_scale=us%m_s_to_L_T**2)
986  ! Use Newton's method to find a new estimate of lam_n
987  dlam = - det / ddet
988  lam_n = lam_n + dlam
989  if (abs(dlam) < tol2*lam_1) then
990  ! calculate nth mode speed
991  if (lam_n > 0.0) cn(i,j,m+1) = 1.0 / sqrt(lam_n)
992  exit
993  endif ! within tol
994  enddo ! itt-loop
995  enddo ! n-loop
996  else
997  cn(i,j,2:nmodes) = 0.0 ! else too small to worry about
998  endif ! if nmodes>1 .and. kc>nmodes .and. c1>c1_thresh
999  else
1000  cn(i,j,:) = 0.0
1001  endif ! if more than 2 layers
1002  endif ! if drxh_sum < 0
1003  else
1004  cn(i,j,:) = 0.0 ! This is a land point.
1005  endif ! if not land
1006  enddo ! i-loop
1007  enddo ! j-loop
1008 

References mom_error_handler::mom_error(), and tridiag_det().

Referenced by mom_diabatic_driver::diabatic().

Here is the call graph for this function:
Here is the caller graph for this function: