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, 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, 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)
 Calculate the determinant of a tridiagonal matrix with diagonals a,b-lam,c where lam is constant across rows. 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
[in]bLeading diagonal
[in]cUpper diagonal
[in]lamScalar subtracted from leading diagonal
[in,out]yRHS on entry, result on exit

Definition at line 464 of file MOM_wave_speed.F90.

Referenced by wave_speed().

464  integer, intent(in) :: n !< Number of rows of matrix
465  real, dimension(n), intent(in) :: a !< Lower diagonal
466  real, dimension(n), intent(in) :: b !< Leading diagonal
467  real, dimension(n), intent(in) :: c !< Upper diagonal
468  real, intent(in) :: lam !< Scalar subtracted from leading diagonal
469  real, dimension(n), intent(inout) :: y !< RHS on entry, result on exit
470  ! Local variables
471  integer :: k, l
472  real :: beta(n), yy(n), lambda
473  lambda = lam
474  beta(1) = b(1) - lambda
475  if (beta(1)==0.) then ! lam was chosen too perfectly
476  ! Change lambda and redo this first row
477  lambda = (1. + 1.e-5) * lambda
478  beta(1) = b(1) - lambda
479  endif
480  beta(1) = 1. / beta(1)
481  yy(1) = y(1)
482  do k = 2, n
483  beta(k) = ( b(k) - lambda ) - a(k) * c(k-1) * beta(k-1)
484  if (beta(k)==0.) then ! lam was chosen too perfectly
485  ! Change lambda and redo everything up to row k
486  lambda = (1. + 1.e-5) * lambda
487  beta(1) = 1. / ( b(1) - lambda )
488  do l = 2, k
489  beta(l) = 1. / ( ( b(l) - lambda ) - a(l) * c(l-1) * beta(l-1) )
490  yy(l) = y(l) - a(l) * yy(l-1) * beta(l-1)
491  enddo
492  else
493  beta(k) = 1. / beta(k)
494  endif
495  yy(k) = y(k) - a(k) * yy(k-1) * beta(k-1)
496  enddo
497  y(n) = yy(n) * beta(n)
498  do k = n-1, 1, -1
499  y(k) = ( yy(k) - c(k) * y(k+1) ) * beta(k)
500  enddo
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 
)
private

Calculate the determinant of a tridiagonal matrix with diagonals a,b-lam,c where lam is constant across rows.

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

Definition at line 1047 of file MOM_wave_speed.F90.

Referenced by wave_speeds().

1047  real, dimension(:), intent(in) :: a !< Lower diagonal of matrix (first entry = 0)
1048  real, dimension(:), intent(in) :: b !< Leading diagonal of matrix (excluding lam)
1049  real, dimension(:), intent(in) :: c !< Upper diagonal of matrix (last entry = 0)
1050  integer, intent(in) :: nrows !< Size of matrix
1051  real, intent(in) :: lam !< Value subtracted from b
1052  real, intent(out):: det_out !< Determinant
1053  real, intent(out):: ddet_out !< Derivative of determinant w.r.t. lam
1054  ! Local variables
1055  real, dimension(nrows) :: det ! value of recursion function
1056  real, dimension(nrows) :: ddet ! value of recursion function for derivative
1057  real, parameter:: rescale = 1024.0**4 ! max value of determinant allowed before rescaling
1058  real :: i_rescale ! inverse of rescale
1059  integer :: n ! row (layer interface) index
1060 
1061  if (size(b) .ne. nrows) call mom_error(warning, "Diagonal b must be same length as nrows.")
1062  if (size(a) .ne. nrows) call mom_error(warning, "Diagonal a must be same length as nrows.")
1063  if (size(c) .ne. nrows) call mom_error(warning, "Diagonal c must be same length as nrows.")
1064 
1065  i_rescale = 1.0/rescale
1066 
1067  det(1) = 1.0 ; ddet(1) = 0.0
1068  det(2) = b(2)-lam ; ddet(2) = -1.0
1069  do n=3,nrows
1070  det(n) = (b(n)-lam)*det(n-1) - (a(n)*c(n-1))*det(n-2)
1071  ddet(n) = (b(n)-lam)*ddet(n-1) - (a(n)*c(n-1))*ddet(n-2) - det(n-1)
1072  ! Rescale det & ddet if det is getting too large.
1073  if (abs(det(n)) > rescale) then
1074  det(n) = i_rescale*det(n) ; det(n-1) = i_rescale*det(n-1)
1075  ddet(n) = i_rescale*ddet(n) ; ddet(n-1) = i_rescale*ddet(n-1)
1076  endif
1077  enddo
1078  det_out = det(nrows)
1079  ddet_out = ddet(nrows)
1080 
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,
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]hLayer thickness (m or kg/m2)
[in]tvThermodynamic variables
[out]cg1First mode internal wave speed (m/s)
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.
[out]modal_structureNormalized model structure (non-dim)

Definition at line 45 of file MOM_wave_speed.F90.

References mom_eos::calculate_density_derivs(), 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().

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

Definition at line 1085 of file MOM_wave_speed.F90.

References wave_speed_set_param().

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

1085  type(wave_speed_cs), pointer :: cs !< Control structure for MOM_wave_speed
1086  logical, optional, intent(in) :: use_ebt_mode !< If true, use the equivalent
1087  !! barotropic mode instead of the first baroclinic mode.
1088  real, optional, intent(in) :: mono_n2_column_fraction !< The lower fraction of water column over which N2 is limited
1089  !! as monotonic for the purposes of calculating vertical modal structure.
1090  real, optional, intent(in) :: mono_n2_depth !< The depth below which N2 is limited
1091  !! as monotonic for the purposes of calculating vertical modal structure.
1092 ! This include declares and sets the variable "version".
1093 #include "version_variable.h"
1094  character(len=40) :: mdl = "MOM_wave_speed" ! This module's name.
1095 
1096  if (associated(cs)) then
1097  call mom_error(warning, "wave_speed_init called with an "// &
1098  "associated control structure.")
1099  return
1100  else ; allocate(cs) ; endif
1101 
1102  ! Write all relevant parameters to the model log.
1103  call log_version(mdl, version)
1104 
1105  call wave_speed_set_param(cs, use_ebt_mode=use_ebt_mode, mono_n2_column_fraction=mono_n2_column_fraction)
1106 
1107  call initialize_remapping(cs%remapping_CS, 'PLM', boundary_extrapolation=.false.)
1108 
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 vertical modal structure.
[in]mono_n2_depthThe depth below which N2 is limited as monotonic for the purposes of calculating vertical modal structure.

Definition at line 1113 of file MOM_wave_speed.F90.

Referenced by wave_speed_init().

1113  type(wave_speed_cs), pointer :: cs !< Control structure for MOM_wave_speed
1114  logical, optional, intent(in) :: use_ebt_mode !< If true, use the equivalent
1115  !! barotropic mode instead of the first baroclinic mode.
1116  real, optional, intent(in) :: mono_n2_column_fraction !< The lower fraction of water column over which N2 is limited
1117  !! as monotonic for the purposes of calculating vertical modal structure.
1118  real, optional, intent(in) :: mono_n2_depth !< The depth below which N2 is limited
1119  !! as monotonic for the purposes of calculating vertical modal structure.
1120 
1121  if (.not.associated(cs)) call mom_error(fatal, &
1122  "wave_speed_set_param called with an associated control structure.")
1123 
1124  if (present(use_ebt_mode)) cs%use_ebt_mode=use_ebt_mode
1125  if (present(mono_n2_column_fraction)) cs%mono_N2_column_fraction=mono_n2_column_fraction
1126  if (present(mono_n2_depth)) cs%mono_N2_depth=mono_n2_depth
1127 
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,
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]hLayer thickness (m or kg/m2)
[in]tvThermodynamic variables
[in]nmodesNumber of modes
[out]cnWaves speeds (m/s)
csControl structure for MOM_wave_speed
[in]full_halosIf true, do the calculation over the entire computational domain.

Definition at line 505 of file MOM_wave_speed.F90.

References tridiag_det().

Referenced by mom_diabatic_driver::diabatic().

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