MOM6
MOM_wave_speed.F90
1 !> Routines for calculating baroclinic wave speeds
2 module mom_wave_speed
3 
4 ! This file is part of MOM6. See LICENSE.md for the license.
5 
6 use mom_diag_mediator, only : post_data, query_averaging_enabled, diag_ctrl
7 use mom_error_handler, only : mom_error, fatal, warning
8 use mom_file_parser, only : log_version
9 use mom_grid, only : ocean_grid_type
10 use mom_remapping, only : remapping_cs, initialize_remapping, remapping_core_h
11 use mom_unit_scaling, only : unit_scale_type
12 use mom_variables, only : thermo_var_ptrs
13 use mom_verticalgrid, only : verticalgrid_type
14 use mom_eos, only : calculate_density, calculate_density_derivs
15 
16 implicit none ; private
17 
18 #include <MOM_memory.h>
19 
20 public wave_speed, wave_speeds, wave_speed_init, wave_speed_set_param
21 
22 ! A note on unit descriptions in comments: MOM6 uses units that can be rescaled for dimensional
23 ! consistency testing. These are noted in comments with units like Z, H, L, and T, along with
24 ! their mks counterparts with notation like "a velocity [Z T-1 ~> m s-1]". If the units
25 ! vary with the Boussinesq approximation, the Boussinesq variant is given first.
26 
27 !> Control structure for MOM_wave_speed
28 type, public :: wave_speed_cs ; private
29  logical :: use_ebt_mode = .false. !< If true, calculate the equivalent barotropic wave speed instead
30  !! of the first baroclinic wave speed.
31  !! This parameter controls the default behavior of wave_speed() which
32  !! can be overridden by optional arguments.
33  logical :: better_cg1_est = .false. !< If true, use an improved estimate of the first mode
34  !! internal wave speed.
35  real :: mono_n2_column_fraction = 0. !< The lower fraction of water column over which N2 is limited as
36  !! monotonic for the purposes of calculating the equivalent barotropic
37  !! wave speed. This parameter controls the default behavior of
38  !! wave_speed() which can be overridden by optional arguments.
39  real :: mono_n2_depth = -1. !< The depth below which N2 is limited as monotonic for the purposes of
40  !! calculating the equivalent barotropic wave speed [Z ~> m].
41  !! This parameter controls the default behavior of wave_speed() which
42  !! can be overridden by optional arguments.
43  real :: min_speed2 = 0. !< The minimum mode 1 internal wave speed squared [L2 T-2 ~> m2 s-2]
44  real :: wave_speed_tol = 0.001 !< The fractional tolerance with which to solve for the wave
45  !! speeds [nondim]
46  type(remapping_cs) :: remapping_cs !< Used for vertical remapping when calculating equivalent barotropic
47  !! mode structure.
48  logical :: remap_answers_2018 = .true. !< If true, use the order of arithmetic and expressions that
49  !! recover the remapping answers from 2018. If false, use more
50  !! robust forms of the same remapping expressions.
51  type(diag_ctrl), pointer :: diag !< Diagnostics control structure
52 end type wave_speed_cs
53 
54 contains
55 
56 !> Calculates the wave speed of the first baroclinic mode.
57 subroutine wave_speed(h, tv, G, GV, US, cg1, CS, full_halos, use_ebt_mode, mono_N2_column_fraction, &
58  mono_N2_depth, modal_structure, better_speed_est, min_speed, wave_speed_tol)
59  type(ocean_grid_type), intent(in) :: g !< Ocean grid structure
60  type(verticalgrid_type), intent(in) :: gv !< Vertical grid structure
61  type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
62  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
63  intent(in) :: h !< Layer thickness [H ~> m or kg m-2]
64  type(thermo_var_ptrs), intent(in) :: tv !< Thermodynamic variables
65  real, dimension(SZI_(G),SZJ_(G)), intent(out) :: cg1 !< First mode internal wave speed [L T-1 ~> m s-1]
66  type(wave_speed_cs), pointer :: cs !< Control structure for MOM_wave_speed
67  logical, optional, intent(in) :: full_halos !< If true, do the calculation
68  !! over the entire computational domain.
69  logical, optional, intent(in) :: use_ebt_mode !< If true, use the equivalent
70  !! barotropic mode instead of the first baroclinic mode.
71  real, optional, intent(in) :: mono_n2_column_fraction !< The lower fraction
72  !! of water column over which N2 is limited as monotonic
73  !! for the purposes of calculating vertical modal structure.
74  real, optional, intent(in) :: mono_n2_depth !< A depth below which N2 is limited as
75  !! monotonic for the purposes of calculating vertical
76  !! modal structure [Z ~> m].
77  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
78  optional, intent(out) :: modal_structure !< Normalized model structure [nondim]
79  logical, optional, intent(in) :: better_speed_est !< If true, use a more robust estimate of the first
80  !! mode speed as the starting point for iterations.
81  real, optional, intent(in) :: min_speed !< If present, set a floor in the first mode speed
82  !! below which 0 is returned [L T-1 ~> m s-1].
83  real, optional, intent(in) :: wave_speed_tol !< The fractional tolerance for finding the
84  !! wave speeds [nondim]
85 
86  ! Local variables
87  real, dimension(SZK_(G)+1) :: &
88  drho_dt, & ! Partial derivative of density with temperature [R degC-1 ~> kg m-3 degC-1]
89  drho_ds, & ! Partial derivative of density with salinity [R ppt-1 ~> kg m-3 ppt-1]
90  pres, & ! Interface pressure [R L2 T-2 ~> Pa]
91  t_int, & ! Temperature interpolated to interfaces [degC]
92  s_int, & ! Salinity interpolated to interfaces [ppt]
93  h_top, & ! The distance of each filtered interface from the ocean surface [Z ~> m]
94  h_bot, & ! The distance of each filtered interface from the bottom [Z ~> m]
95  gprime ! The reduced gravity across each interface [L2 Z-1 T-2 ~> m s-2].
96  real, dimension(SZK_(G)) :: &
97  igl, igu ! The inverse of the reduced gravity across an interface times
98  ! the thickness of the layer below (Igl) or above (Igu) it, in [T2 L-2 ~> s2 m-2].
99  real, dimension(SZK_(G),SZI_(G)) :: &
100  hf, & ! Layer thicknesses after very thin layers are combined [Z ~> m]
101  tf, & ! Layer temperatures after very thin layers are combined [degC]
102  sf, & ! Layer salinities after very thin layers are combined [ppt]
103  rf ! Layer densities after very thin layers are combined [R ~> kg m-3]
104  real, dimension(SZK_(G)) :: &
105  hc, & ! A column of layer thicknesses after convective istabilities are removed [Z ~> m]
106  tc, & ! A column of layer temperatures after convective istabilities are removed [degC]
107  sc, & ! A column of layer salinites after convective istabilities are removed [ppt]
108  rc, & ! A column of layer densities after convective istabilities are removed [R ~> kg m-3]
109  hc_h ! Hc(:) rescaled from Z to thickness units [H ~> m or kg m-2]
110  real :: i_htot ! The inverse of the total filtered thicknesses [Z ~> m]
111  real :: det, ddet, detkm1, detkm2, ddetkm1, ddetkm2
112  real :: lam ! The eigenvalue [T2 L-2 ~> s m-1]
113  real :: dlam ! The change in estimates of the eigenvalue [T2 L-2 ~> s m-1]
114  real :: lam0 ! The first guess of the eigenvalue [T2 L-2 ~> s m-1]
115  real :: min_h_frac ! [nondim]
116  real :: z_to_pres ! A conversion factor from thicknesses to pressure [R L2 T-2 Z-1 ~> Pa m-1]
117  real, dimension(SZI_(G)) :: &
118  htot, hmin, & ! Thicknesses [Z ~> m]
119  h_here, & ! A thickness [Z ~> m]
120  hxt_here, & ! A layer integrated temperature [degC Z ~> degC m]
121  hxs_here, & ! A layer integrated salinity [ppt Z ~> ppt m]
122  hxr_here ! A layer integrated density [R Z ~> kg m-2]
123  real :: speed2_tot ! overestimate of the mode-1 speed squared [L2 T-2 ~> m2 s-2]
124  real :: cg1_min2 ! A floor in the squared first mode speed below which 0 is returned [L2 T-2 ~> m2 s-2]
125  real :: i_hnew ! The inverse of a new layer thickness [Z-1 ~> m-1]
126  real :: drxh_sum ! The sum of density differences across interfaces times thicknesses [R Z ~> kg m-2]
127  real :: l2_to_z2 ! A scaling factor squared from units of lateral distances to depths [Z2 L-2 ~> 1].
128  real, pointer, dimension(:,:,:) :: t => null(), s => null()
129  real :: g_rho0 ! G_Earth/Rho0 [L2 T-2 Z-1 R-1 ~> m4 s-2 kg-1].
130  real :: c2_scale ! A scaling factor for wave speeds to help control the growth of the determinant
131  ! and its derivative with lam between rows of the Thomas algorithm solver. The
132  ! exact value should not matter for the final result if it is an even power of 2.
133  real :: tol_hfrac ! Layers that together are smaller than this fraction of
134  ! the total water column can be merged for efficiency.
135  real :: tol_solve ! The fractional tolerance with which to solve for the wave speeds [nondim]
136  real :: tol_merge ! The fractional change in estimated wave speed that is allowed
137  ! when deciding to merge layers in the calculation [nondim]
138  real :: rescale, i_rescale
139  integer :: kf(szi_(g)) ! The number of active layers after filtering.
140  integer, parameter :: max_itt = 10
141  real :: lam_it(max_itt), det_it(max_itt), ddet_it(max_itt)
142  logical :: use_eos ! If true, density is calculated from T & S using an equation of state.
143  logical :: better_est ! If true, use an improved estimate of the first mode internal wave speed.
144  logical :: merge ! If true, merge the current layer with the one above.
145  integer :: kc ! The number of layers in the column after merging
146  integer :: i, j, k, k2, itt, is, ie, js, je, nz
147  real :: hw, sum_hc
148  real :: gp ! A limited local copy of gprime [L2 Z-1 T-2 ~> m s-2]
149  real :: n2min ! A minimum buoyancy frequency [T-2 ~> s-2]
150  logical :: l_use_ebt_mode, calc_modal_structure
151  real :: l_mono_n2_column_fraction, l_mono_n2_depth
152  real :: mode_struct(szk_(g)), ms_min, ms_max, ms_sq
153 
154  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
155 
156  if (.not. associated(cs)) call mom_error(fatal, "MOM_wave_speed: "// &
157  "Module must be initialized before it is used.")
158  if (present(full_halos)) then ; if (full_halos) then
159  is = g%isd ; ie = g%ied ; js = g%jsd ; je = g%jed
160  endif ; endif
161 
162  l2_to_z2 = us%L_to_Z**2
163 
164  l_use_ebt_mode = cs%use_ebt_mode
165  if (present(use_ebt_mode)) l_use_ebt_mode = use_ebt_mode
166  l_mono_n2_column_fraction = cs%mono_N2_column_fraction
167  if (present(mono_n2_column_fraction)) l_mono_n2_column_fraction = mono_n2_column_fraction
168  l_mono_n2_depth = cs%mono_N2_depth
169  if (present(mono_n2_depth)) l_mono_n2_depth = mono_n2_depth
170  calc_modal_structure = l_use_ebt_mode
171  if (present(modal_structure)) calc_modal_structure = .true.
172  if (calc_modal_structure) then
173  do k=1,nz; do j=js,je; do i=is,ie
174  modal_structure(i,j,k) = 0.0
175  enddo ; enddo ; enddo
176  endif
177 
178  s => tv%S ; t => tv%T
179  g_rho0 = gv%g_Earth / gv%Rho0
180  ! Simplifying the following could change answers at roundoff.
181  z_to_pres = gv%Z_to_H * (gv%H_to_RZ * gv%g_Earth)
182  use_eos = associated(tv%eqn_of_state)
183 
184  better_est = cs%better_cg1_est ; if (present(better_speed_est)) better_est = better_speed_est
185 
186  if (better_est) then
187  tol_solve = cs%wave_speed_tol ; if (present(wave_speed_tol)) tol_solve = wave_speed_tol
188  tol_hfrac = 0.1*tol_solve ; tol_merge = tol_solve / real(nz)
189  else
190  tol_solve = 0.001 ; tol_hfrac = 0.0001 ; tol_merge = 0.001
191  endif
192 
193  ! The rescaling below can control the growth of the determinant provided that
194  ! (tol_merge*cg1_min2/c2_scale > I_rescale). For default values, this suggests a stable lower
195  ! bound on min_speed of sqrt(nz/(tol_solve*rescale)) or 3e2/1024**2 = 2.9e-4 m/s for 90 layers.
196  ! The upper bound on the rate of increase in the determinant is g'H/c2_scale < rescale or in the
197  ! worst possible oceanic case of g'H < 0.5*10m/s2*1e4m = 5.e4 m2/s2 < 1024**2*c2_scale, suggesting
198  ! that c2_scale can safely be set to 1/(16*1024**2), which would decrease the stable floor on
199  ! min_speed to ~6.9e-8 m/s for 90 layers or 2.33e-7 m/s for 1000 layers.
200  cg1_min2 = cs%min_speed2 ; if (present(min_speed)) cg1_min2 = min_speed**2
201  rescale = 1024.0**4 ; i_rescale = 1.0/rescale
202  c2_scale = us%m_s_to_L_T**2 / 4096.0**2 ! Other powers of 2 give identical results.
203 
204  min_h_frac = tol_hfrac / real(nz)
205 !$OMP parallel do default(none) shared(is,ie,js,je,nz,h,G,GV,US,min_h_frac,use_EOS,T,S,tv,&
206 !$OMP calc_modal_structure,l_use_ebt_mode,modal_structure, &
207 !$OMP l_mono_N2_column_fraction,l_mono_N2_depth,CS, &
208 !$OMP Z_to_pres,cg1,g_Rho0,rescale,I_rescale,L2_to_Z2, &
209 !$OMP better_est,cg1_min2,tol_merge,tol_solve,c2_scale) &
210 !$OMP private(htot,hmin,kf,H_here,HxT_here,HxS_here,HxR_here, &
211 !$OMP Hf,Tf,Sf,Rf,pres,T_int,S_int,drho_dT,drho_dS, &
212 !$OMP drxh_sum,kc,Hc,Hc_H,tC,sc,I_Hnew,gprime,&
213 !$OMP Rc,speed2_tot,Igl,Igu,lam0,lam,lam_it,dlam, &
214 !$OMP mode_struct,sum_hc,N2min,gp,hw, &
215 !$OMP ms_min,ms_max,ms_sq,H_top,H_bot,I_Htot,merge, &
216 !$OMP det,ddet,detKm1,ddetKm1,detKm2,ddetKm2,det_it,ddet_it)
217  do j=js,je
218  ! First merge very thin layers with the one above (or below if they are
219  ! at the top). This also transposes the row order so that columns can
220  ! be worked upon one at a time.
221  do i=is,ie ; htot(i) = 0.0 ; enddo
222  do k=1,nz ; do i=is,ie ; htot(i) = htot(i) + h(i,j,k)*gv%H_to_Z ; enddo ; enddo
223 
224  do i=is,ie
225  hmin(i) = htot(i)*min_h_frac ; kf(i) = 1 ; h_here(i) = 0.0
226  hxt_here(i) = 0.0 ; hxs_here(i) = 0.0 ; hxr_here(i) = 0.0
227  enddo
228  if (use_eos) then
229  do k=1,nz ; do i=is,ie
230  if ((h_here(i) > hmin(i)) .and. (h(i,j,k)*gv%H_to_Z > hmin(i))) then
231  hf(kf(i),i) = h_here(i)
232  tf(kf(i),i) = hxt_here(i) / h_here(i)
233  sf(kf(i),i) = hxs_here(i) / h_here(i)
234  kf(i) = kf(i) + 1
235 
236  ! Start a new layer
237  h_here(i) = h(i,j,k)*gv%H_to_Z
238  hxt_here(i) = (h(i,j,k)*gv%H_to_Z)*t(i,j,k)
239  hxs_here(i) = (h(i,j,k)*gv%H_to_Z)*s(i,j,k)
240  else
241  h_here(i) = h_here(i) + h(i,j,k)*gv%H_to_Z
242  hxt_here(i) = hxt_here(i) + (h(i,j,k)*gv%H_to_Z)*t(i,j,k)
243  hxs_here(i) = hxs_here(i) + (h(i,j,k)*gv%H_to_Z)*s(i,j,k)
244  endif
245  enddo ; enddo
246  do i=is,ie ; if (h_here(i) > 0.0) then
247  hf(kf(i),i) = h_here(i)
248  tf(kf(i),i) = hxt_here(i) / h_here(i)
249  sf(kf(i),i) = hxs_here(i) / h_here(i)
250  endif ; enddo
251  else
252  do k=1,nz ; do i=is,ie
253  if ((h_here(i) > hmin(i)) .and. (h(i,j,k)*gv%H_to_Z > hmin(i))) then
254  hf(kf(i),i) = h_here(i) ; rf(kf(i),i) = hxr_here(i) / h_here(i)
255  kf(i) = kf(i) + 1
256 
257  ! Start a new layer
258  h_here(i) = h(i,j,k)*gv%H_to_Z
259  hxr_here(i) = (h(i,j,k)*gv%H_to_Z)*gv%Rlay(k)
260  else
261  h_here(i) = h_here(i) + h(i,j,k)*gv%H_to_Z
262  hxr_here(i) = hxr_here(i) + (h(i,j,k)*gv%H_to_Z)*gv%Rlay(k)
263  endif
264  enddo ; enddo
265  do i=is,ie ; if (h_here(i) > 0.0) then
266  hf(kf(i),i) = h_here(i) ; rf(kf(i),i) = hxr_here(i) / h_here(i)
267  endif ; enddo
268  endif
269 
270  ! From this point, we can work on individual columns without causing memory to have page faults.
271  do i=is,ie ; if (g%mask2dT(i,j) > 0.5) then
272  if (use_eos) then
273  pres(1) = 0.0 ; h_top(1) = 0.0
274  do k=2,kf(i)
275  pres(k) = pres(k-1) + z_to_pres*hf(k-1,i)
276  t_int(k) = 0.5*(tf(k,i)+tf(k-1,i))
277  s_int(k) = 0.5*(sf(k,i)+sf(k-1,i))
278  h_top(k) = h_top(k-1) + hf(k-1,i)
279  enddo
280  call calculate_density_derivs(t_int, s_int, pres, drho_dt, drho_ds, &
281  tv%eqn_of_state, (/2,kf(i)/) )
282 
283  ! Sum the reduced gravities to find out how small a density difference is negligibly small.
284  drxh_sum = 0.0
285  if (better_est) then
286  ! This is an estimate that is correct for the non-EBT mode for 2 or 3 layers, or for
287  ! clusters of massless layers at interfaces that can be grouped into 2 or 3 layers.
288  ! For a uniform stratification and a huge number of layers uniformly distributed in
289  ! density, this estimate is too large (as is desired) by a factor of pi^2/6 ~= 1.64.
290  if (h_top(kf(i)) > 0.0) then
291  i_htot = 1.0 / (h_top(kf(i)) + hf(kf(i),i)) ! = 1.0 / (H_top(K) + H_bot(K)) for all K.
292  h_bot(kf(i)+1) = 0.0
293  do k=kf(i),2,-1
294  h_bot(k) = h_bot(k+1) + hf(k,i)
295  drxh_sum = drxh_sum + ((h_top(k) * h_bot(k)) * i_htot) * &
296  max(0.0, drho_dt(k)*(tf(k,i)-tf(k-1,i)) + drho_ds(k)*(sf(k,i)-sf(k-1,i)))
297  enddo
298  endif
299  else
300  ! This estimate is problematic in that it goes like 1/nz for a large number of layers,
301  ! but it is an overestimate (as desired) for a small number of layers, by at a factor
302  ! of (H1+H2)**2/(H1*H2) >= 4 for two thick layers.
303  do k=2,kf(i)
304  drxh_sum = drxh_sum + 0.5*(hf(k-1,i)+hf(k,i)) * &
305  max(0.0, drho_dt(k)*(tf(k,i)-tf(k-1,i)) + drho_ds(k)*(sf(k,i)-sf(k-1,i)))
306  enddo
307  endif
308  else
309  drxh_sum = 0.0
310  if (better_est) then
311  h_top(1) = 0.0
312  do k=2,kf(i) ; h_top(k) = h_top(k-1) + hf(k-1,i) ; enddo
313  if (h_top(kf(i)) > 0.0) then
314  i_htot = 1.0 / (h_top(kf(i)) + hf(kf(i),i)) ! = 1.0 / (H_top(K) + H_bot(K)) for all K.
315  h_bot(kf(i)+1) = 0.0
316  do k=kf(i),2,-1
317  h_bot(k) = h_bot(k+1) + hf(k,i)
318  drxh_sum = drxh_sum + ((h_top(k) * h_bot(k)) * i_htot) * max(0.0,rf(k,i)-rf(k-1,i))
319  enddo
320  endif
321  else
322  do k=2,kf(i)
323  drxh_sum = drxh_sum + 0.5*(hf(k-1,i)+hf(k,i)) * max(0.0,rf(k,i)-rf(k-1,i))
324  enddo
325  endif
326  endif
327 
328  ! Find gprime across each internal interface, taking care of convective instabilities by
329  ! merging layers. If the estimated wave speed is too small, simply return zero.
330  if (g_rho0 * drxh_sum <= cg1_min2) then
331  cg1(i,j) = 0.0
332  if (present(modal_structure)) modal_structure(i,j,:) = 0.
333  else
334  ! Merge layers to eliminate convective instabilities or exceedingly
335  ! small reduced gravities. Merging layers reduces the estimated wave speed by
336  ! (rho(2)-rho(1))*h(1)*h(2) / H_tot.
337  if (use_eos) then
338  kc = 1
339  hc(1) = hf(1,i) ; tc(1) = tf(1,i) ; sc(1) = sf(1,i)
340  do k=2,kf(i)
341  if (better_est) then
342  merge = ((drho_dt(k)*(tf(k,i)-tc(kc)) + drho_ds(k)*(sf(k,i)-sc(kc))) * &
343  ((hc(kc) * hf(k,i))*i_htot) < 2.0 * tol_merge*drxh_sum)
344  else
345  merge = ((drho_dt(k)*(tf(k,i)-tc(kc)) + drho_ds(k)*(sf(k,i)-sc(kc))) * &
346  (hc(kc) + hf(k,i)) < 2.0 * tol_merge*drxh_sum)
347  endif
348  if (merge) then
349  ! Merge this layer with the one above and backtrack.
350  i_hnew = 1.0 / (hc(kc) + hf(k,i))
351  tc(kc) = (hc(kc)*tc(kc) + hf(k,i)*tf(k,i)) * i_hnew
352  sc(kc) = (hc(kc)*sc(kc) + hf(k,i)*sf(k,i)) * i_hnew
353  hc(kc) = (hc(kc) + hf(k,i))
354  ! Backtrack to remove any convective instabilities above... Note
355  ! that the tolerance is a factor of two larger, to avoid limit how
356  ! far back we go.
357  do k2=kc,2,-1
358  if (better_est) then
359  merge = ((drho_dt(k2)*(tc(k2)-tc(k2-1)) + drho_ds(k2)*(sc(k2)-sc(k2-1))) * &
360  ((hc(k2) * hc(k2-1))*i_htot) < tol_merge*drxh_sum)
361  else
362  merge = ((drho_dt(k2)*(tc(k2)-tc(k2-1)) + drho_ds(k2)*(sc(k2)-sc(k2-1))) * &
363  (hc(k2) + hc(k2-1)) < tol_merge*drxh_sum)
364  endif
365  if (merge) then
366  ! Merge the two bottommost layers. At this point kc = k2.
367  i_hnew = 1.0 / (hc(kc) + hc(kc-1))
368  tc(kc-1) = (hc(kc)*tc(kc) + hc(kc-1)*tc(kc-1)) * i_hnew
369  sc(kc-1) = (hc(kc)*sc(kc) + hc(kc-1)*sc(kc-1)) * i_hnew
370  hc(kc-1) = (hc(kc) + hc(kc-1))
371  kc = kc - 1
372  else ; exit ; endif
373  enddo
374  else
375  ! Add a new layer to the column.
376  kc = kc + 1
377  drho_ds(kc) = drho_ds(k) ; drho_dt(kc) = drho_dt(k)
378  tc(kc) = tf(k,i) ; sc(kc) = sf(k,i) ; hc(kc) = hf(k,i)
379  endif
380  enddo
381  ! At this point there are kc layers and the gprimes should be positive.
382  do k=2,kc ! Revisit this if non-Boussinesq.
383  gprime(k) = g_rho0 * (drho_dt(k)*(tc(k)-tc(k-1)) + drho_ds(k)*(sc(k)-sc(k-1)))
384  enddo
385  else ! .not.use_EOS
386  ! Do the same with density directly...
387  kc = 1
388  hc(1) = hf(1,i) ; rc(1) = rf(1,i)
389  do k=2,kf(i)
390  if (better_est) then
391  merge = ((rf(k,i) - rc(kc)) * ((hc(kc) * hf(k,i))*i_htot) < 2.0*tol_merge*drxh_sum)
392  else
393  merge = ((rf(k,i) - rc(kc)) * (hc(kc) + hf(k,i)) < 2.0*tol_merge*drxh_sum)
394  endif
395  if (merge) then
396  ! Merge this layer with the one above and backtrack.
397  rc(kc) = (hc(kc)*rc(kc) + hf(k,i)*rf(k,i)) / (hc(kc) + hf(k,i))
398  hc(kc) = (hc(kc) + hf(k,i))
399  ! Backtrack to remove any convective instabilities above... Note
400  ! that the tolerance is a factor of two larger, to avoid limit how
401  ! far back we go.
402  do k2=kc,2,-1
403  if (better_est) then
404  merge = ((rc(k2)-rc(k2-1)) * ((hc(k2) * hc(k2-1))*i_htot) < tol_merge*drxh_sum)
405  else
406  merge = ((rc(k2)-rc(k2-1)) * (hc(k2)+hc(k2-1)) < tol_merge*drxh_sum)
407  endif
408  if (merge) then
409  ! Merge the two bottommost layers. At this point kc = k2.
410  rc(kc-1) = (hc(kc)*rc(kc) + hc(kc-1)*rc(kc-1)) / (hc(kc) + hc(kc-1))
411  hc(kc-1) = (hc(kc) + hc(kc-1))
412  kc = kc - 1
413  else ; exit ; endif
414  enddo
415  else
416  ! Add a new layer to the column.
417  kc = kc + 1
418  rc(kc) = rf(k,i) ; hc(kc) = hf(k,i)
419  endif
420  enddo
421  ! At this point there are kc layers and the gprimes should be positive.
422  do k=2,kc ! Revisit this if non-Boussinesq.
423  gprime(k) = g_rho0 * (rc(k)-rc(k-1))
424  enddo
425  endif ! use_EOS
426 
427  ! Sum the contributions from all of the interfaces to give an over-estimate
428  ! of the first-mode wave speed. Also populate Igl and Igu which are the
429  ! non-leading diagonals of the tridiagonal matrix.
430  if (kc >= 2) then
431  speed2_tot = 0.0
432  if (better_est) then
433  h_top(1) = 0.0 ; h_bot(kc+1) = 0.0
434  do k=2,kc+1 ; h_top(k) = h_top(k-1) + hc(k-1) ; enddo
435  do k=kc,2,-1 ; h_bot(k) = h_bot(k+1) + hc(k) ; enddo
436  i_htot = 0.0 ; if (h_top(kc+1) > 0.0) i_htot = 1.0 / h_top(kc+1)
437  endif
438 
439  if (l_use_ebt_mode) then
440  igu(1) = 0. ! Neumann condition for pressure modes
441  sum_hc = hc(1)
442  n2min = l2_to_z2*gprime(2)/hc(1)
443  do k=2,kc
444  hw = 0.5*(hc(k-1)+hc(k))
445  gp = gprime(k)
446  if (l_mono_n2_column_fraction>0. .or. l_mono_n2_depth>=0.) then
447  if ( ((g%bathyT(i,j)-sum_hc < l_mono_n2_column_fraction*g%bathyT(i,j)) .or. &
448  ((l_mono_n2_depth >= 0.) .and. (sum_hc > l_mono_n2_depth))) .and. &
449  (l2_to_z2*gp > n2min*hw) ) then
450  ! Filters out regions where N2 increases with depth but only in a lower fraction
451  ! of the water column or below a certain depth.
452  gp = us%Z_to_L**2 * (n2min*hw)
453  else
454  n2min = l2_to_z2 * gp/hw
455  endif
456  endif
457  igu(k) = 1.0/(gp*hc(k))
458  igl(k-1) = 1.0/(gp*hc(k-1))
459  sum_hc = sum_hc + hc(k)
460  if (better_est) then
461  ! Estimate that the ebt_mode is sqrt(2) times the speed of the flat bottom modes.
462  speed2_tot = speed2_tot + 2.0 * gprime(k)*((h_top(k) * h_bot(k)) * i_htot)
463  else ! The ebt_mode wave should be faster than the flat-bottom mode, so 0.707 should be > 1?
464  speed2_tot = speed2_tot + gprime(k)*(hc(k-1)+hc(k))*0.707
465  endif
466  enddo
467  !Igl(kc) = 0. ! Neumann condition for pressure modes
468  igl(kc) = 2.*igu(kc) ! Dirichlet condition for pressure modes
469  else ! .not. l_use_ebt_mode
470  do k=2,kc
471  igl(k) = 1.0/(gprime(k)*hc(k)) ; igu(k) = 1.0/(gprime(k)*hc(k-1))
472  if (better_est) then
473  speed2_tot = speed2_tot + gprime(k)*((h_top(k) * h_bot(k)) * i_htot)
474  else
475  speed2_tot = speed2_tot + gprime(k)*(hc(k-1)+hc(k))
476  endif
477  enddo
478  endif
479 
480  if (calc_modal_structure) then
481  mode_struct(:) = 0.
482  mode_struct(1:kc) = 1. ! Uniform flow, first guess
483  endif
484 
485  ! Under estimate the first eigenvalue (overestimate the speed) to start with.
486  if (calc_modal_structure) then
487  lam0 = 0.5 / speed2_tot ; lam = lam0
488  else
489  lam0 = 1.0 / speed2_tot ; lam = lam0
490  endif
491  ! Find the determinant and its derivative with lam.
492  do itt=1,max_itt
493  lam_it(itt) = lam
494  if (l_use_ebt_mode) then
495  ! This initialization of det,ddet imply Neumann boundary conditions for horizontal
496  ! velocity or pressure modes, so that first 3 rows of the matrix are
497  ! / b(1)-lam igl(1) 0 0 0 ... \
498  ! | igu(2) b(2)-lam igl(2) 0 0 ... |
499  ! | 0 igu(3) b(3)-lam igl(3) 0 ... |
500  ! The last two rows of the pressure equation matrix are
501  ! | ... 0 igu(kc-1) b(kc-1)-lam igl(kc-1) |
502  ! \ ... 0 0 igu(kc) b(kc)-lam /
503  call tridiag_det(igu, igl, 1, kc, lam, det, ddet, row_scale=c2_scale)
504  else
505  ! This initialization of det,ddet imply Dirichlet boundary conditions for vertical
506  ! velocity modes, so that first 3 rows of the matrix are
507  ! / b(2)-lam igl(2) 0 0 0 ... |
508  ! | igu(3) b(3)-lam igl(3) 0 0 ... |
509  ! | 0 igu(4) b(4)-lam igl(4) 0 ... |
510  ! The last three rows of the w equation matrix are
511  ! | ... 0 igu(kc-2) b(kc-2)-lam igl(kc-2) 0 |
512  ! | ... 0 0 igu(kc-1) b(kc-1)-lam igl(kc-1) |
513  ! \ ... 0 0 0 igu(kc) b(kc)-lam /
514  call tridiag_det(igu, igl, 2, kc, lam, det, ddet, row_scale=c2_scale)
515  endif
516  ! Use Newton's method iteration to find a new estimate of lam.
517  det_it(itt) = det ; ddet_it(itt) = ddet
518 
519  if ((ddet >= 0.0) .or. (-det > -0.5*lam*ddet)) then
520  ! lam was not an under-estimate, as intended, so Newton's method
521  ! may not be reliable; lam must be reduced, but not by more
522  ! than half.
523  lam = 0.5 * lam
524  dlam = -lam
525  else ! Newton's method is OK.
526  dlam = - det / ddet
527  lam = lam + dlam
528  endif
529 
530  if (calc_modal_structure) then
531  call tdma6(kc, igu, igl, lam, mode_struct)
532  ms_min = mode_struct(1)
533  ms_max = mode_struct(1)
534  ms_sq = mode_struct(1)**2
535  do k = 2,kc
536  ms_min = min(ms_min, mode_struct(k))
537  ms_max = max(ms_max, mode_struct(k))
538  ms_sq = ms_sq + mode_struct(k)**2
539  enddo
540  if (ms_min<0. .and. ms_max>0.) then ! Any zero crossings => lam is too high
541  lam = 0.5 * ( lam - dlam )
542  dlam = -lam
543  mode_struct(1:kc) = abs(mode_struct(1:kc)) / sqrt( ms_sq )
544  else
545  mode_struct(1:kc) = mode_struct(1:kc) / sqrt( ms_sq )
546  endif
547  endif
548 
549  if (abs(dlam) < tol_solve*lam) exit
550  enddo
551 
552  cg1(i,j) = 0.0
553  if (lam > 0.0) cg1(i,j) = 1.0 / sqrt(lam)
554 
555  if (present(modal_structure)) then
556  if (mode_struct(1)/=0.) then ! Normalize
557  mode_struct(1:kc) = mode_struct(1:kc) / mode_struct(1)
558  else
559  mode_struct(1:kc)=0.
560  endif
561  ! Note that remapping_core_h requires that the same units be used
562  ! for both the source and target grid thicknesses, here [H ~> m or kg m-2].
563  do k = 1,kc
564  hc_h(k) = gv%Z_to_H * hc(k)
565  enddo
566  if (cs%remap_answers_2018) then
567  call remapping_core_h(cs%remapping_CS, kc, hc_h(:), mode_struct, &
568  nz, h(i,j,:), modal_structure(i,j,:), &
569  1.0e-30*gv%m_to_H, 1.0e-10*gv%m_to_H)
570  else
571  call remapping_core_h(cs%remapping_CS, kc, hc_h(:), mode_struct, &
572  nz, h(i,j,:), modal_structure(i,j,:), &
573  gv%H_subroundoff, gv%H_subroundoff)
574  endif
575  endif
576  else
577  cg1(i,j) = 0.0
578  if (present(modal_structure)) modal_structure(i,j,:) = 0.
579  endif
580  endif ! cg1 /= 0.0
581  else
582  cg1(i,j) = 0.0 ! This is a land point.
583  if (present(modal_structure)) modal_structure(i,j,:) = 0.
584  endif ; enddo ! i-loop
585  enddo ! j-loop
586 
587 end subroutine wave_speed
588 
589 !> Solve a non-symmetric tridiagonal problem with the sum of the upper and lower diagnonals minus a
590 !! scalar contribution as the leading diagonal.
591 !! This uses the Thomas algorithm rather than the Hallberg algorithm since the matrix is not symmetric.
592 subroutine tdma6(n, a, c, lam, y)
593  integer, intent(in) :: n !< Number of rows of matrix
594  real, dimension(:), intent(in) :: a !< Lower diagonal [T2 L-2 ~> s2 m-2]
595  real, dimension(:), intent(in) :: c !< Upper diagonal [T2 L-2 ~> s2 m-2]
596  real, intent(in) :: lam !< Scalar subtracted from leading diagonal [T2 L-2 ~> s2 m-2]
597  real, dimension(:), intent(inout) :: y !< RHS on entry, result on exit
598 
599  ! Local variables
600  real :: lambda ! A temporary variable in [T2 L-2 ~> s2 m-2]
601  real :: beta(n) ! A temporary variable in [T2 L-2 ~> s2 m-2]
602  real :: I_beta(n) ! A temporary variable in [L2 T-2 ~> m2 s-2]
603  real :: yy(n) ! A temporary variable with the same units as y on entry.
604  integer :: k, m
605 
606  lambda = lam
607  beta(1) = (a(1)+c(1)) - lambda
608  if (beta(1)==0.) then ! lam was chosen too perfectly
609  ! Change lambda and redo this first row
610  lambda = (1. + 1.e-5) * lambda
611  beta(1) = (a(1)+c(1)) - lambda
612  endif
613  i_beta(1) = 1. / beta(1)
614  yy(1) = y(1)
615  do k = 2, n
616  beta(k) = ( (a(k)+c(k)) - lambda ) - a(k) * c(k-1) * i_beta(k-1)
617  ! Perhaps the following 0 needs to become a tolerance to handle underflow?
618  if (beta(k)==0.) then ! lam was chosen too perfectly
619  ! Change lambda and redo everything up to row k
620  lambda = (1. + 1.e-5) * lambda
621  i_beta(1) = 1. / ( (a(1)+c(1)) - lambda )
622  do m = 2, k
623  i_beta(m) = 1. / ( ( (a(m)+c(m)) - lambda ) - a(m) * c(m-1) * i_beta(m-1) )
624  yy(m) = y(m) + a(m) * yy(m-1) * i_beta(m-1)
625  enddo
626  else
627  i_beta(k) = 1. / beta(k)
628  endif
629  yy(k) = y(k) + a(k) * yy(k-1) * i_beta(k-1)
630  enddo
631  ! The units of y change by a factor of [L2 T-2] in the following lines.
632  y(n) = yy(n) * i_beta(n)
633  do k = n-1, 1, -1
634  y(k) = ( yy(k) + c(k) * y(k+1) ) * i_beta(k)
635  enddo
636 
637 end subroutine tdma6
638 
639 !> Calculates the wave speeds for the first few barolinic modes.
640 subroutine wave_speeds(h, tv, G, GV, US, nmodes, cn, CS, full_halos, better_speed_est, &
641  min_speed, wave_speed_tol)
642  type(ocean_grid_type), intent(in) :: g !< Ocean grid structure
643  type(verticalgrid_type), intent(in) :: gv !< Vertical grid structure
644  type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
645  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h !< Layer thickness [H ~> m or kg m-2]
646  type(thermo_var_ptrs), intent(in) :: tv !< Thermodynamic variables
647  integer, intent(in) :: nmodes !< Number of modes
648  real, dimension(G%isd:G%ied,G%jsd:G%jed,nmodes), intent(out) :: cn !< Waves speeds [L T-1 ~> m s-1]
649  type(wave_speed_cs), optional, pointer :: cs !< Control structure for MOM_wave_speed
650  logical, optional, intent(in) :: full_halos !< If true, do the calculation
651  !! over the entire computational domain.
652  logical, optional, intent(in) :: better_speed_est !< If true, use a more robust estimate of the first
653  !! mode speed as the starting point for iterations.
654  real, optional, intent(in) :: min_speed !< If present, set a floor in the first mode speed
655  !! below which 0 is returned [L T-1 ~> m s-1].
656  real, optional, intent(in) :: wave_speed_tol !< The fractional tolerance for finding the
657  !! wave speeds [nondim]
658 
659  ! Local variables
660  real, dimension(SZK_(G)+1) :: &
661  drho_dt, & ! Partial derivative of density with temperature [R degC-1 ~> kg m-3 degC-1]
662  drho_ds, & ! Partial derivative of density with salinity [R ppt-1 ~> kg m-3 ppt-1]
663  pres, & ! Interface pressure [R L2 T-2 ~> Pa]
664  t_int, & ! Temperature interpolated to interfaces [degC]
665  s_int, & ! Salinity interpolated to interfaces [ppt]
666  h_top, & ! The distance of each filtered interface from the ocean surface [Z ~> m]
667  h_bot, & ! The distance of each filtered interface from the bottom [Z ~> m]
668  gprime ! The reduced gravity across each interface [L2 Z-1 T-2 ~> m s-2].
669  real, dimension(SZK_(G),SZI_(G)) :: &
670  hf, & ! Layer thicknesses after very thin layers are combined [Z ~> m]
671  tf, & ! Layer temperatures after very thin layers are combined [degC]
672  sf, & ! Layer salinities after very thin layers are combined [ppt]
673  rf ! Layer densities after very thin layers are combined [R ~> kg m-3]
674  real, dimension(SZK_(G)) :: &
675  igl, igu, & ! The inverse of the reduced gravity across an interface times
676  ! the thickness of the layer below (Igl) or above (Igu) it, in [T2 L-2 ~> s2 m-2].
677  hc, & ! A column of layer thicknesses after convective istabilities are removed [Z ~> m]
678  tc, & ! A column of layer temperatures after convective istabilities are removed [degC]
679  sc, & ! A column of layer salinites after convective istabilities are removed [ppt]
680  rc ! A column of layer densities after convective istabilities are removed [R ~> kg m-3]
681  real :: i_htot ! The inverse of the total filtered thicknesses [Z ~> m]
682  real :: c1_thresh ! if c1 is below this value, don't bother calculating
683  ! cn values for higher modes [L T-1 ~> m s-1]
684  real :: c2_scale ! A scaling factor for wave speeds to help control the growth of the determinant
685  ! and its derivative with lam between rows of the Thomas algorithm solver. The
686  ! exact value should not matter for the final result if it is an even power of 2.
687  real :: det, ddet ! determinant & its derivative of eigen system
688  real :: lam_1 ! approximate mode-1 eigenvalue [T2 L-2 ~> s2 m-2]
689  real :: lam_n ! approximate mode-n eigenvalue [T2 L-2 ~> s2 m-2]
690  real :: dlam ! The change in estimates of the eigenvalue [T2 L-2 ~> s m-1]
691  real :: lammin ! minimum lam value for root searching range [T2 L-2 ~> s2 m-2]
692  real :: lammax ! maximum lam value for root searching range [T2 L-2 ~> s2 m-2]
693  real :: laminc ! width of moving window for root searching [T2 L-2 ~> s2 m-2]
694  real :: det_l,det_r ! determinant value at left and right of window
695  real :: ddet_l,ddet_r ! derivative of determinant at left and right of window
696  real :: det_sub,ddet_sub! derivative of determinant at subinterval endpoint
697  real :: xl,xr ! lam guesses at left and right of window [T2 L-2 ~> s2 m-2]
698  real :: xl_sub ! lam guess at left of subinterval window [T2 L-2 ~> s2 m-2]
699  real,dimension(nmodes) :: &
700  xbl,xbr ! lam guesses bracketing a zero-crossing (root) [T2 L-2 ~> s2 m-2]
701  integer :: numint ! number of widows (intervals) in root searching range
702  integer :: nrootsfound ! number of extra roots found (not including 1st root)
703  real :: z_to_pres ! A conversion factor from thicknesses to pressure [R L2 T-2 Z-1 ~> Pa m-1]
704  real, dimension(SZI_(G)) :: &
705  htot, hmin, & ! Thicknesses [Z ~> m]
706  h_here, & ! A thickness [Z ~> m]
707  hxt_here, & ! A layer integrated temperature [degC Z ~> degC m]
708  hxs_here, & ! A layer integrated salinity [ppt Z ~> ppt m]
709  hxr_here ! A layer integrated density [R Z ~> kg m-2]
710  real :: speed2_tot ! overestimate of the mode-1 speed squared [L2 T-2 ~> m2 s-2]
711  real :: speed2_min ! minimum mode speed (squared) to consider in root searching [L2 T-2 ~> m2 s-2]
712  real :: cg1_min2 ! A floor in the squared first mode speed below which 0 is returned [L2 T-2 ~> m2 s-2]
713  real, parameter :: reduct_factor = 0.5
714  ! A factor used in setting speed2_min [nondim]
715  real :: i_hnew ! The inverse of a new layer thickness [Z-1 ~> m-1]
716  real :: drxh_sum ! The sum of density differences across interfaces times thicknesses [R Z ~> kg m-2]
717  real, pointer, dimension(:,:,:) :: t => null(), s => null()
718  real :: g_rho0 ! G_Earth/Rho0 [L2 T-2 Z-1 R-1 ~> m4 s-2 kg-1].
719  real :: tol_hfrac ! Layers that together are smaller than this fraction of
720  ! the total water column can be merged for efficiency.
721  real :: min_h_frac ! tol_Hfrac divided by the total number of layers [nondim].
722  real :: tol_solve ! The fractional tolerance with which to solve for the wave speeds [nondim].
723  real :: tol_merge ! The fractional change in estimated wave speed that is allowed
724  ! when deciding to merge layers in the calculation [nondim]
725  integer :: kf(szi_(g)) ! The number of active layers after filtering.
726  integer, parameter :: max_itt = 10
727  logical :: use_eos ! If true, density is calculated from T & S using the equation of state.
728  logical :: better_est ! If true, use an improved estimate of the first mode internal wave speed.
729  logical :: merge ! If true, merge the current layer with the one above.
730  integer :: nsub ! number of subintervals used for root finding
731  integer, parameter :: sub_it_max = 4
732  ! maximum number of times to subdivide interval
733  ! for root finding (# intervals = 2**sub_it_max)
734  logical :: sub_rootfound ! if true, subdivision has located root
735  integer :: kc ! The number of layers in the column after merging
736  integer :: sub, sub_it
737  integer :: i, j, k, k2, itt, is, ie, js, je, nz, row, iint, m, ig, jg
738 
739  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
740 
741  if (present(cs)) then
742  if (.not. associated(cs)) call mom_error(fatal, "MOM_wave_speed: "// &
743  "Module must be initialized before it is used.")
744  endif
745 
746  if (present(full_halos)) then ; if (full_halos) then
747  is = g%isd ; ie = g%ied ; js = g%jsd ; je = g%jed
748  endif ; endif
749 
750  s => tv%S ; t => tv%T
751  g_rho0 = gv%g_Earth / gv%Rho0
752  ! Simplifying the following could change answers at roundoff.
753  z_to_pres = gv%Z_to_H * (gv%H_to_RZ * gv%g_Earth)
754  use_eos = associated(tv%eqn_of_state)
755  c1_thresh = 0.01*us%m_s_to_L_T
756  c2_scale = us%m_s_to_L_T**2 / 4096.0**2 ! Other powers of 2 give identical results.
757 
758  better_est = .false. ; if (present(cs)) better_est = cs%better_cg1_est
759  if (present(better_speed_est)) better_est = better_speed_est
760  if (better_est) then
761  tol_solve = 0.001 ; if (present(cs)) tol_solve = cs%wave_speed_tol
762  if (present(wave_speed_tol)) tol_solve = wave_speed_tol
763  tol_hfrac = 0.1*tol_solve ; tol_merge = tol_solve / real(nz)
764  else
765  tol_solve = 0.001 ; tol_hfrac = 0.0001 ; tol_merge = 0.001
766  endif
767  cg1_min2 = 0.0 ; if (present(cs)) cg1_min2 = cs%min_speed2
768  if (present(min_speed)) cg1_min2 = min_speed**2
769 
770  ! Zero out all wave speeds. Values over land or for columns that are too weakly stratified
771  ! are not changed from this zero value.
772  cn(:,:,:) = 0.0
773 
774  min_h_frac = tol_hfrac / real(nz)
775  !$OMP parallel do default(private) shared(is,ie,js,je,nz,h,G,GV,US,min_h_frac,use_EOS,T,S, &
776  !$OMP Z_to_pres,tv,cn,g_Rho0,nmodes,cg1_min2,better_est, &
777  !$OMP c1_thresh,tol_solve,tol_merge)
778  do j=js,je
779  ! First merge very thin layers with the one above (or below if they are
780  ! at the top). This also transposes the row order so that columns can
781  ! be worked upon one at a time.
782  do i=is,ie ; htot(i) = 0.0 ; enddo
783  do k=1,nz ; do i=is,ie ; htot(i) = htot(i) + h(i,j,k)*gv%H_to_Z ; enddo ; enddo
784 
785  do i=is,ie
786  hmin(i) = htot(i)*min_h_frac ; kf(i) = 1 ; h_here(i) = 0.0
787  hxt_here(i) = 0.0 ; hxs_here(i) = 0.0 ; hxr_here(i) = 0.0
788  enddo
789  if (use_eos) then
790  do k=1,nz ; do i=is,ie
791  if ((h_here(i) > hmin(i)) .and. (h(i,j,k)*gv%H_to_Z > hmin(i))) then
792  hf(kf(i),i) = h_here(i)
793  tf(kf(i),i) = hxt_here(i) / h_here(i)
794  sf(kf(i),i) = hxs_here(i) / h_here(i)
795  kf(i) = kf(i) + 1
796 
797  ! Start a new layer
798  h_here(i) = h(i,j,k)*gv%H_to_Z
799  hxt_here(i) = (h(i,j,k)*gv%H_to_Z)*t(i,j,k)
800  hxs_here(i) = (h(i,j,k)*gv%H_to_Z)*s(i,j,k)
801  else
802  h_here(i) = h_here(i) + h(i,j,k)*gv%H_to_Z
803  hxt_here(i) = hxt_here(i) + (h(i,j,k)*gv%H_to_Z)*t(i,j,k)
804  hxs_here(i) = hxs_here(i) + (h(i,j,k)*gv%H_to_Z)*s(i,j,k)
805  endif
806  enddo ; enddo
807  do i=is,ie ; if (h_here(i) > 0.0) then
808  hf(kf(i),i) = h_here(i)
809  tf(kf(i),i) = hxt_here(i) / h_here(i)
810  sf(kf(i),i) = hxs_here(i) / h_here(i)
811  endif ; enddo
812  else
813  do k=1,nz ; do i=is,ie
814  if ((h_here(i) > hmin(i)) .and. (h(i,j,k)*gv%H_to_Z > hmin(i))) then
815  hf(kf(i),i) = h_here(i) ; rf(kf(i),i) = hxr_here(i) / h_here(i)
816  kf(i) = kf(i) + 1
817 
818  ! Start a new layer
819  h_here(i) = h(i,j,k)*gv%H_to_Z
820  hxr_here(i) = (h(i,j,k)*gv%H_to_Z)*gv%Rlay(k)
821  else
822  h_here(i) = h_here(i) + h(i,j,k)*gv%H_to_Z
823  hxr_here(i) = hxr_here(i) + (h(i,j,k)*gv%H_to_Z)*gv%Rlay(k)
824  endif
825  enddo ; enddo
826  do i=is,ie ; if (h_here(i) > 0.0) then
827  hf(kf(i),i) = h_here(i) ; rf(kf(i),i) = hxr_here(i) / h_here(i)
828  endif ; enddo
829  endif
830 
831  ! From this point, we can work on individual columns without causing memory to have page faults.
832  do i=is,ie
833  if (g%mask2dT(i,j) > 0.5) then
834  if (use_eos) then
835  pres(1) = 0.0 ; h_top(1) = 0.0
836  do k=2,kf(i)
837  pres(k) = pres(k-1) + z_to_pres*hf(k-1,i)
838  t_int(k) = 0.5*(tf(k,i)+tf(k-1,i))
839  s_int(k) = 0.5*(sf(k,i)+sf(k-1,i))
840  h_top(k) = h_top(k-1) + hf(k-1,i)
841  enddo
842  call calculate_density_derivs(t_int, s_int, pres, drho_dt, drho_ds, &
843  tv%eqn_of_state, (/2,kf(i)/) )
844 
845  ! Sum the reduced gravities to find out how small a density difference is negligibly small.
846  drxh_sum = 0.0
847  if (better_est) then
848  ! This is an estimate that is correct for the non-EBT mode for 2 or 3 layers, or for
849  ! clusters of massless layers at interfaces that can be grouped into 2 or 3 layers.
850  ! For a uniform stratification and a huge number of layers uniformly distributed in
851  ! density, this estimate is too large (as is desired) by a factor of pi^2/6 ~= 1.64.
852  if (h_top(kf(i)) > 0.0) then
853  i_htot = 1.0 / (h_top(kf(i)) + hf(kf(i),i)) ! = 1.0 / (H_top(K) + H_bot(K)) for all K.
854  h_bot(kf(i)+1) = 0.0
855  do k=kf(i),2,-1
856  h_bot(k) = h_bot(k+1) + hf(k,i)
857  drxh_sum = drxh_sum + ((h_top(k) * h_bot(k)) * i_htot) * &
858  max(0.0, drho_dt(k)*(tf(k,i)-tf(k-1,i)) + drho_ds(k)*(sf(k,i)-sf(k-1,i)))
859  enddo
860  endif
861  else
862  ! This estimate is problematic in that it goes like 1/nz for a large number of layers,
863  ! but it is an overestimate (as desired) for a small number of layers, by at a factor
864  ! of (H1+H2)**2/(H1*H2) >= 4 for two thick layers.
865  do k=2,kf(i)
866  drxh_sum = drxh_sum + 0.5*(hf(k-1,i)+hf(k,i)) * &
867  max(0.0, drho_dt(k)*(tf(k,i)-tf(k-1,i)) + drho_ds(k)*(sf(k,i)-sf(k-1,i)))
868  enddo
869  endif
870  else
871  drxh_sum = 0.0
872  if (better_est) then
873  h_top(1) = 0.0
874  do k=2,kf(i) ; h_top(k) = h_top(k-1) + hf(k-1,i) ; enddo
875  if (h_top(kf(i)) > 0.0) then
876  i_htot = 1.0 / (h_top(kf(i)) + hf(kf(i),i)) ! = 1.0 / (H_top(K) + H_bot(K)) for all K.
877  h_bot(kf(i)+1) = 0.0
878  do k=kf(i),2,-1
879  h_bot(k) = h_bot(k+1) + hf(k,i)
880  drxh_sum = drxh_sum + ((h_top(k) * h_bot(k)) * i_htot) * max(0.0,rf(k,i)-rf(k-1,i))
881  enddo
882  endif
883  else
884  do k=2,kf(i)
885  drxh_sum = drxh_sum + 0.5*(hf(k-1,i)+hf(k,i)) * max(0.0,rf(k,i)-rf(k-1,i))
886  enddo
887  endif
888  endif
889 
890  ! Find gprime across each internal interface, taking care of convective
891  ! instabilities by merging layers.
892  if (g_rho0 * drxh_sum > cg1_min2) then
893  ! Merge layers to eliminate convective instabilities or exceedingly
894  ! small reduced gravities. Merging layers reduces the estimated wave speed by
895  ! (rho(2)-rho(1))*h(1)*h(2) / H_tot.
896  if (use_eos) then
897  kc = 1
898  hc(1) = hf(1,i) ; tc(1) = tf(1,i) ; sc(1) = sf(1,i)
899  do k=2,kf(i)
900  if (better_est) then
901  merge = ((drho_dt(k)*(tf(k,i)-tc(kc)) + drho_ds(k)*(sf(k,i)-sc(kc))) * &
902  ((hc(kc) * hf(k,i))*i_htot) < 2.0 * tol_merge*drxh_sum)
903  else
904  merge = ((drho_dt(k)*(tf(k,i)-tc(kc)) + drho_ds(k)*(sf(k,i)-sc(kc))) * &
905  (hc(kc) + hf(k,i)) < 2.0 * tol_merge*drxh_sum)
906  endif
907  if (merge) then
908  ! Merge this layer with the one above and backtrack.
909  i_hnew = 1.0 / (hc(kc) + hf(k,i))
910  tc(kc) = (hc(kc)*tc(kc) + hf(k,i)*tf(k,i)) * i_hnew
911  sc(kc) = (hc(kc)*sc(kc) + hf(k,i)*sf(k,i)) * i_hnew
912  hc(kc) = (hc(kc) + hf(k,i))
913  ! Backtrack to remove any convective instabilities above... Note
914  ! that the tolerance is a factor of two larger, to avoid limit how
915  ! far back we go.
916  do k2=kc,2,-1
917  if (better_est) then
918  merge = ((drho_dt(k2)*(tc(k2)-tc(k2-1)) + drho_ds(k2)*(sc(k2)-sc(k2-1))) * &
919  ((hc(k2) * hc(k2-1))*i_htot) < tol_merge*drxh_sum)
920  else
921  merge = ((drho_dt(k2)*(tc(k2)-tc(k2-1)) + drho_ds(k2)*(sc(k2)-sc(k2-1))) * &
922  (hc(k2) + hc(k2-1)) < tol_merge*drxh_sum)
923  endif
924  if (merge) then
925  ! Merge the two bottommost layers. At this point kc = k2.
926  i_hnew = 1.0 / (hc(kc) + hc(kc-1))
927  tc(kc-1) = (hc(kc)*tc(kc) + hc(kc-1)*tc(kc-1)) * i_hnew
928  sc(kc-1) = (hc(kc)*sc(kc) + hc(kc-1)*sc(kc-1)) * i_hnew
929  hc(kc-1) = (hc(kc) + hc(kc-1))
930  kc = kc - 1
931  else ; exit ; endif
932  enddo
933  else
934  ! Add a new layer to the column.
935  kc = kc + 1
936  drho_ds(kc) = drho_ds(k) ; drho_dt(kc) = drho_dt(k)
937  tc(kc) = tf(k,i) ; sc(kc) = sf(k,i) ; hc(kc) = hf(k,i)
938  endif
939  enddo
940  ! At this point there are kc layers and the gprimes should be positive.
941  do k=2,kc ! Revisit this if non-Boussinesq.
942  gprime(k) = g_rho0 * (drho_dt(k)*(tc(k)-tc(k-1)) + drho_ds(k)*(sc(k)-sc(k-1)))
943  enddo
944  else ! .not.use_EOS
945  ! Do the same with density directly...
946  kc = 1
947  hc(1) = hf(1,i) ; rc(1) = rf(1,i)
948  do k=2,kf(i)
949  if (better_est) then
950  merge = ((rf(k,i) - rc(kc)) * ((hc(kc) * hf(k,i))*i_htot) < 2.0 * tol_merge*drxh_sum)
951  else
952  merge = ((rf(k,i) - rc(kc)) * (hc(kc) + hf(k,i)) < 2.0*tol_merge*drxh_sum)
953  endif
954  if (merge) then
955  ! Merge this layer with the one above and backtrack.
956  rc(kc) = (hc(kc)*rc(kc) + hf(k,i)*rf(k,i)) / (hc(kc) + hf(k,i))
957  hc(kc) = (hc(kc) + hf(k,i))
958  ! Backtrack to remove any convective instabilities above... Note
959  ! that the tolerance is a factor of two larger, to avoid limit how
960  ! far back we go.
961  do k2=kc,2,-1
962  if (better_est) then
963  merge = ((rc(k2)-rc(k2-1)) * ((hc(k2) * hc(k2-1))*i_htot) < tol_merge*drxh_sum)
964  else
965  merge = ((rc(k2)-rc(k2-1)) * (hc(k2)+hc(k2-1)) < tol_merge*drxh_sum)
966  endif
967  if (merge) then
968  ! Merge the two bottommost layers. At this point kc = k2.
969  rc(kc-1) = (hc(kc)*rc(kc) + hc(kc-1)*rc(kc-1)) / (hc(kc) + hc(kc-1))
970  hc(kc-1) = (hc(kc) + hc(kc-1))
971  kc = kc - 1
972  else ; exit ; endif
973  enddo
974  else
975  ! Add a new layer to the column.
976  kc = kc + 1
977  rc(kc) = rf(k,i) ; hc(kc) = hf(k,i)
978  endif
979  enddo
980  ! At this point there are kc layers and the gprimes should be positive.
981  do k=2,kc ! Revisit this if non-Boussinesq.
982  gprime(k) = g_rho0 * (rc(k)-rc(k-1))
983  enddo
984  endif ! use_EOS
985 
986  !-----------------NOW FIND WAVE SPEEDS---------------------------------------
987  ! ig = i + G%idg_offset ; jg = j + G%jdg_offset
988  ! Sum the contributions from all of the interfaces to give an over-estimate
989  ! of the first-mode wave speed. Also populate Igl and Igu which are the
990  ! non-leading diagonals of the tridiagonal matrix.
991  if (kc >= 2) then
992  ! initialize speed2_tot
993  speed2_tot = 0.0
994  if (better_est) then
995  h_top(1) = 0.0 ; h_bot(kc+1) = 0.0
996  do k=2,kc+1 ; h_top(k) = h_top(k-1) + hc(k-1) ; enddo
997  do k=kc,2,-1 ; h_bot(k) = h_bot(k+1) + hc(k) ; enddo
998  i_htot = 0.0 ; if (h_top(kc+1) > 0.0) i_htot = 1.0 / h_top(kc+1)
999  endif
1000 
1001  ! Calculate Igu, Igl, depth, and N2 at each interior interface
1002  ! [excludes surface (K=1) and bottom (K=kc+1)]
1003  do k=2,kc
1004  igl(k) = 1.0/(gprime(k)*hc(k)) ; igu(k) = 1.0/(gprime(k)*hc(k-1))
1005  if (better_est) then
1006  speed2_tot = speed2_tot + gprime(k)*((h_top(k) * h_bot(k)) * i_htot)
1007  else
1008  speed2_tot = speed2_tot + gprime(k)*(hc(k-1)+hc(k))
1009  endif
1010  enddo
1011 
1012  ! Under estimate the first eigenvalue (overestimate the speed) to start with.
1013  lam_1 = 1.0 / speed2_tot
1014 
1015  ! Find the first eigen value
1016  do itt=1,max_itt
1017  ! calculate the determinant of (A-lam_1*I)
1018  call tridiag_det(igu, igl, 2, kc, lam_1, det, ddet, row_scale=c2_scale)
1019 
1020  ! If possible, use Newton's method iteration to find a new estimate of lam_1
1021  !det = det_it(itt) ; ddet = ddet_it(itt)
1022  if ((ddet >= 0.0) .or. (-det > -0.5*lam_1*ddet)) then
1023  ! lam_1 was not an under-estimate, as intended, so Newton's method
1024  ! may not be reliable; lam_1 must be reduced, but not by more than half.
1025  lam_1 = 0.5 * lam_1
1026  dlam = -lam_1
1027  else ! Newton's method is OK.
1028  dlam = - det / ddet
1029  lam_1 = lam_1 + dlam
1030  endif
1031 
1032  if (abs(dlam) < tol_solve*lam_1) exit
1033  enddo
1034 
1035  if (lam_1 > 0.0) cn(i,j,1) = 1.0 / sqrt(lam_1)
1036 
1037  ! Find other eigen values if c1 is of significant magnitude, > cn_thresh
1038  nrootsfound = 0 ! number of extra roots found (not including 1st root)
1039  if (nmodes>1 .and. kc>=nmodes+1 .and. cn(i,j,1)>c1_thresh) then
1040  ! Set the the range to look for the other desired eigen values
1041  ! set min value just greater than the 1st root (found above)
1042  lammin = lam_1*(1.0 + tol_solve)
1043  ! set max value based on a low guess at wavespeed for highest mode
1044  speed2_min = (reduct_factor*cn(i,j,1)/real(nmodes))**2
1045  lammax = 1.0 / speed2_min
1046  ! set width of interval (not sure about this - BDM)
1047  laminc = 0.5*lam_1
1048  ! set number of intervals within search range
1049  numint = nint((lammax - lammin)/laminc)
1050 
1051  ! Find intervals containing zero-crossings (roots) of the determinant
1052  ! that are beyond the first root
1053 
1054  ! find det_l of first interval (det at left endpoint)
1055  call tridiag_det(igu, igl, 2, kc, lammin, det_l, ddet_l, row_scale=c2_scale)
1056  ! move interval window looking for zero-crossings************************
1057  do iint=1,numint
1058  xr = lammin + laminc * iint
1059  xl = xr - laminc
1060  call tridiag_det(igu, igl, 2, kc, xr, det_r, ddet_r, row_scale=c2_scale)
1061  if (det_l*det_r < 0.0) then ! if function changes sign
1062  if (det_l*ddet_l < 0.0) then ! if function at left is headed to zero
1063  nrootsfound = nrootsfound + 1
1064  xbl(nrootsfound) = xl
1065  xbr(nrootsfound) = xr
1066  else
1067  ! function changes sign but has a local max/min in interval,
1068  ! try subdividing interval as many times as necessary (or sub_it_max).
1069  ! loop that increases number of subintervals:
1070  !call MOM_error(WARNING, "determinant changes sign"// &
1071  ! "but has a local max/min in interval;"//&
1072  ! " reduce increment in lam.")
1073  ! begin subdivision loop -------------------------------------------
1074  sub_rootfound = .false. ! initialize
1075  do sub_it=1,sub_it_max
1076  nsub = 2**sub_it ! number of subintervals; nsub=2,4,8,...
1077  ! loop over each subinterval:
1078  do sub=1,nsub-1,2 ! only check odds; sub = 1; 1,3; 1,3,5,7;...
1079  xl_sub = xl + laminc/(nsub)*sub
1080  call tridiag_det(igu, igl, 2, kc, xl_sub, det_sub, ddet_sub, &
1081  row_scale=c2_scale)
1082  if (det_sub*det_r < 0.0) then ! if function changes sign
1083  if (det_sub*ddet_sub < 0.0) then ! if function at left is headed to zero
1084  sub_rootfound = .true.
1085  nrootsfound = nrootsfound + 1
1086  xbl(nrootsfound) = xl_sub
1087  xbr(nrootsfound) = xr
1088  exit ! exit sub loop
1089  endif ! headed toward zero
1090  endif ! sign change
1091  enddo ! sub-loop
1092  if (sub_rootfound) exit ! root has been found, exit sub_it loop
1093  ! Otherwise, function changes sign but has a local max/min in one of the
1094  ! sub intervals, try subdividing again unless sub_it_max has been reached.
1095  if (sub_it == sub_it_max) then
1096  call mom_error(warning, "wave_speed: root not found "// &
1097  " after sub_it_max subdivisions of original"// &
1098  " interval.")
1099  endif ! sub_it == sub_it_max
1100  enddo ! sub_it-loop-------------------------------------------------
1101  endif ! det_l*ddet_l < 0.0
1102  endif ! det_l*det_r < 0.0
1103  ! exit iint-loop if all desired roots have been found
1104  if (nrootsfound >= nmodes-1) then
1105  ! exit if all additional roots found
1106  exit
1107  elseif (iint == numint) then
1108  ! oops, lamMax not large enough - could add code to increase (BDM)
1109  ! set unfound modes to zero for now (BDM)
1110  ! cn(i,j,nrootsfound+2:nmodes) = 0.0
1111  else
1112  ! else shift interval and keep looking until nmodes or numint is reached
1113  det_l = det_r
1114  ddet_l = ddet_r
1115  endif
1116  enddo ! iint-loop
1117 
1118  ! Use Newton's method to find the roots within the identified windows
1119  do m=1,nrootsfound ! loop over the root-containing widows (excluding 1st mode)
1120  lam_n = xbl(m) ! first guess is left edge of window
1121  do itt=1,max_itt
1122  ! calculate the determinant of (A-lam_n*I)
1123  call tridiag_det(igu, igl, 2, kc, lam_n, det, ddet, row_scale=c2_scale)
1124  ! Use Newton's method to find a new estimate of lam_n
1125  dlam = - det / ddet
1126  lam_n = lam_n + dlam
1127  if (abs(dlam) < tol_solve*lam_1) exit
1128  enddo ! itt-loop
1129  ! calculate nth mode speed
1130  if (lam_n > 0.0) cn(i,j,m+1) = 1.0 / sqrt(lam_n)
1131  enddo ! n-loop
1132  endif ! if nmodes>1 .and. kc>nmodes .and. c1>c1_thresh
1133  endif ! if more than 2 layers
1134  endif ! if drxh_sum < 0
1135  endif ! if not land
1136  enddo ! i-loop
1137  enddo ! j-loop
1138 
1139 end subroutine wave_speeds
1140 
1141 !> Calculate the determinant of a tridiagonal matrix with diagonals a,b-lam,c and its derivative
1142 !! with lam, where lam is constant across rows. Only the ratio of det to its derivative and their
1143 !! signs are typically used, so internal rescaling by consistent factors are used to avoid
1144 !! over- or underflow.
1145 subroutine tridiag_det(a, c, ks, ke, lam, det, ddet, row_scale)
1146  real, dimension(:), intent(in) :: a !< Lower diagonal of matrix (first entry unused)
1147  real, dimension(:), intent(in) :: c !< Upper diagonal of matrix (last entry unused)
1148  integer, intent(in) :: ks !< Starting index to use in determinant
1149  integer, intent(in) :: ke !< Ending index to use in determinant
1150  real, intent(in) :: lam !< Value subtracted from b
1151  real, intent(out):: det !< Determinant
1152  real, intent(out):: ddet !< Derivative of determinant with lam
1153  real, optional, intent(in) :: row_scale !< A scaling factor of the rows of the
1154  !! matrix to limit the growth of the determinant
1155  ! Local variables
1156  real :: detKm1, detKm2 ! Cumulative value of the determinant for the previous two layers.
1157  real :: ddetKm1, ddetKm2 ! Derivative of the cumulative determinant with lam for the previous two layers.
1158  real, parameter :: rescale = 1024.0**4 ! max value of determinant allowed before rescaling
1159  real :: rscl ! A rescaling factor that is applied succesively to each row.
1160  real :: I_rescale ! inverse of rescale
1161  integer :: k ! row (layer interface) index
1162 
1163  i_rescale = 1.0 / rescale
1164  rscl = 1.0 ; if (present(row_scale)) rscl = row_scale
1165 
1166  detkm1 = 1.0 ; ddetkm1 = 0.0
1167  det = (a(ks)+c(ks)) - lam ; ddet = -1.0
1168  do k=ks+1,ke
1169  ! Shift variables and rescale rows to avoid over- or underflow.
1170  detkm2 = row_scale*detkm1 ; ddetkm2 = row_scale*ddetkm1
1171  detkm1 = row_scale*det ; ddetkm1 = row_scale*ddet
1172 
1173  det = ((a(k)+c(k))-lam)*detkm1 - (a(k)*c(k-1))*detkm2
1174  ddet = ((a(k)+c(k))-lam)*ddetkm1 - (a(k)*c(k-1))*ddetkm2 - detkm1
1175 
1176  ! Rescale det & ddet if det is getting too large or too small.
1177  if (abs(det) > rescale) then
1178  det = i_rescale*det ; detkm1 = i_rescale*detkm1
1179  ddet = i_rescale*ddet ; ddetkm1 = i_rescale*ddetkm1
1180  elseif (abs(det) < i_rescale) then
1181  det = rescale*det ; detkm1 = rescale*detkm1
1182  ddet = rescale*ddet ; ddetkm1 = rescale*ddetkm1
1183  endif
1184  enddo
1185 
1186 end subroutine tridiag_det
1187 
1188 !> Initialize control structure for MOM_wave_speed
1189 subroutine wave_speed_init(CS, use_ebt_mode, mono_N2_column_fraction, mono_N2_depth, remap_answers_2018, &
1190  better_speed_est, min_speed, wave_speed_tol)
1191  type(wave_speed_cs), pointer :: cs !< Control structure for MOM_wave_speed
1192  logical, optional, intent(in) :: use_ebt_mode !< If true, use the equivalent
1193  !! barotropic mode instead of the first baroclinic mode.
1194  real, optional, intent(in) :: mono_n2_column_fraction !< The lower fraction of water column over
1195  !! which N2 is limited as monotonic for the purposes of
1196  !! calculating the vertical modal structure.
1197  real, optional, intent(in) :: mono_n2_depth !< The depth below which N2 is limited
1198  !! as monotonic for the purposes of calculating the
1199  !! vertical modal structure [Z ~> m].
1200  logical, optional, intent(in) :: remap_answers_2018 !< If true, use the order of arithmetic and expressions
1201  !! that recover the remapping answers from 2018. Otherwise
1202  !! use more robust but mathematically equivalent expressions.
1203  logical, optional, intent(in) :: better_speed_est !< If true, use a more robust estimate of the first
1204  !! mode speed as the starting point for iterations.
1205  real, optional, intent(in) :: min_speed !< If present, set a floor in the first mode speed
1206  !! below which 0 is returned [L T-1 ~> m s-1].
1207  real, optional, intent(in) :: wave_speed_tol !< The fractional tolerance for finding the
1208  !! wave speeds [nondim]
1209 
1210  ! This include declares and sets the variable "version".
1211 # include "version_variable.h"
1212  character(len=40) :: mdl = "MOM_wave_speed" ! This module's name.
1213 
1214  if (associated(cs)) then
1215  call mom_error(warning, "wave_speed_init called with an "// &
1216  "associated control structure.")
1217  return
1218  else ; allocate(cs) ; endif
1219 
1220  ! Write all relevant parameters to the model log.
1221  call log_version(mdl, version)
1222 
1223  call wave_speed_set_param(cs, use_ebt_mode=use_ebt_mode, mono_n2_column_fraction=mono_n2_column_fraction, &
1224  better_speed_est=better_speed_est, min_speed=min_speed, wave_speed_tol=wave_speed_tol)
1225 
1226  call initialize_remapping(cs%remapping_CS, 'PLM', boundary_extrapolation=.false., &
1227  answers_2018=cs%remap_answers_2018)
1228 
1229 end subroutine wave_speed_init
1230 
1231 !> Sets internal parameters for MOM_wave_speed
1232 subroutine wave_speed_set_param(CS, use_ebt_mode, mono_N2_column_fraction, mono_N2_depth, remap_answers_2018, &
1233  better_speed_est, min_speed, wave_speed_tol)
1234  type(wave_speed_cs), pointer :: cs !< Control structure for MOM_wave_speed
1235  logical, optional, intent(in) :: use_ebt_mode !< If true, use the equivalent
1236  !! barotropic mode instead of the first baroclinic mode.
1237  real, optional, intent(in) :: mono_n2_column_fraction !< The lower fraction of water column over
1238  !! which N2 is limited as monotonic for the purposes of
1239  !! calculating the vertical modal structure.
1240  real, optional, intent(in) :: mono_n2_depth !< The depth below which N2 is limited
1241  !! as monotonic for the purposes of calculating the
1242  !! vertical modal structure [Z ~> m].
1243  logical, optional, intent(in) :: remap_answers_2018 !< If true, use the order of arithmetic and expressions
1244  !! that recover the remapping answers from 2018. Otherwise
1245  !! use more robust but mathematically equivalent expressions.
1246  logical, optional, intent(in) :: better_speed_est !< If true, use a more robust estimate of the first
1247  !! mode speed as the starting point for iterations.
1248  real, optional, intent(in) :: min_speed !< If present, set a floor in the first mode speed
1249  !! below which 0 is returned [L T-1 ~> m s-1].
1250  real, optional, intent(in) :: wave_speed_tol !< The fractional tolerance for finding the
1251  !! wave speeds [nondim]
1252 
1253  if (.not.associated(cs)) call mom_error(fatal, &
1254  "wave_speed_set_param called with an associated control structure.")
1255 
1256  if (present(use_ebt_mode)) cs%use_ebt_mode = use_ebt_mode
1257  if (present(mono_n2_column_fraction)) cs%mono_N2_column_fraction = mono_n2_column_fraction
1258  if (present(mono_n2_depth)) cs%mono_N2_depth = mono_n2_depth
1259  if (present(remap_answers_2018)) cs%remap_answers_2018 = remap_answers_2018
1260  if (present(better_speed_est)) cs%better_cg1_est = better_speed_est
1261  if (present(min_speed)) cs%min_speed2 = min_speed**2
1262  if (present(wave_speed_tol)) cs%wave_speed_tol = wave_speed_tol
1263 
1264 end subroutine wave_speed_set_param
1265 
1266 !> \namespace mom_wave_speed
1267 
1268 !!
1269 !! Subroutine wave_speed() solves for the first baroclinic mode wave speed. (It could
1270 !! solve for all the wave speeds, but the iterative approach taken here means
1271 !! that this is not particularly efficient.)
1272 !!
1273 !! If `e(k)` is the perturbation interface height, this means solving for the
1274 !! smallest eigenvalue (`lam` = 1/c^2) of the system
1275 !!
1276 !! \verbatim
1277 !! -Igu(k)*e(k-1) + (Igu(k)+Igl(k)-lam)*e(k) - Igl(k)*e(k+1) = 0.0
1278 !! \endverbatim
1279 !!
1280 !! with rigid lid boundary conditions e(1) = e(nz+1) = 0.0 giving
1281 !!
1282 !! \verbatim
1283 !! (Igu(2)+Igl(2)-lam)*e(2) - Igl(2)*e(3) = 0.0
1284 !! -Igu(nz)*e(nz-1) + (Igu(nz)+Igl(nz)-lam)*e(nz) = 0.0
1285 !! \endverbatim
1286 !!
1287 !! Here
1288 !! \verbatim
1289 !! Igl(k) = 1.0/(gprime(K)*h(k)) ; Igu(k) = 1.0/(gprime(K)*h(k-1))
1290 !! \endverbatim
1291 !!
1292 !! Alternately, these same eigenvalues can be found from the second smallest
1293 !! eigenvalue of the Montgomery potential (M(k)) calculation:
1294 !!
1295 !! \verbatim
1296 !! -Igl(k)*M(k-1) + (Igl(k)+Igu(k+1)-lam)*M(k) - Igu(k+1)*M(k+1) = 0.0
1297 !! \endverbatim
1298 !!
1299 !! with rigid lid and flat bottom boundary conditions
1300 !!
1301 !! \verbatim
1302 !! (Igu(2)-lam)*M(1) - Igu(2)*M(2) = 0.0
1303 !! -Igl(nz)*M(nz-1) + (Igl(nz)-lam)*M(nz) = 0.0
1304 !! \endverbatim
1305 !!
1306 !! Note that the barotropic mode has been eliminated from the rigid lid
1307 !! interface height equations, hence the matrix is one row smaller. Without
1308 !! the rigid lid, the top boundary condition is simpler to implement with
1309 !! the M equations.
1310 
1311 end module mom_wave_speed
mom_verticalgrid
Provides a transparent vertical ocean grid type and supporting routines.
Definition: MOM_verticalGrid.F90:2
mom_file_parser::log_version
An overloaded interface to log version information about modules.
Definition: MOM_file_parser.F90:109
mom_diag_mediator
The subroutines here provide convenient wrappers to the fms diag_manager interfaces with additional d...
Definition: MOM_diag_mediator.F90:3
mom_variables::thermo_var_ptrs
Pointers to an assortment of thermodynamic fields that may be available, including potential temperat...
Definition: MOM_variables.F90:80
mom_eos
Provides subroutines for quantities specific to the equation of state.
Definition: MOM_EOS.F90:2
mom_remapping::remapping_cs
Container for remapping parameters.
Definition: MOM_remapping.F90:24
mom_unit_scaling::unit_scale_type
Describes various unit conversion factors.
Definition: MOM_unit_scaling.F90:14
mom_remapping
Provides column-wise vertical remapping functions.
Definition: MOM_remapping.F90:2
mom_diag_mediator::post_data
Make a diagnostic available for averaging or output.
Definition: MOM_diag_mediator.F90:70
mom_verticalgrid::verticalgrid_type
Describes the vertical ocean grid, including unit conversion factors.
Definition: MOM_verticalGrid.F90:24
mom_variables
Provides transparent structures with groups of MOM6 variables and supporting routines.
Definition: MOM_variables.F90:2
mom_wave_speed::wave_speed_cs
Control structure for MOM_wave_speed.
Definition: MOM_wave_speed.F90:28
mom_file_parser
The MOM6 facility to parse input files for runtime parameters.
Definition: MOM_file_parser.F90:2
mom_eos::calculate_density_derivs
Calculate the derivatives of density with temperature and salinity from T, S, and P.
Definition: MOM_EOS.F90:81
mom_grid
Provides the ocean grid type.
Definition: MOM_grid.F90:2
mom_unit_scaling
Provides a transparent unit rescaling type to facilitate dimensional consistency testing.
Definition: MOM_unit_scaling.F90:2
mom_wave_speed
Routines for calculating baroclinic wave speeds.
Definition: MOM_wave_speed.F90:2
mom_error_handler
Routines for error handling and I/O management.
Definition: MOM_error_handler.F90:2
mom_grid::ocean_grid_type
Ocean grid type. See mom_grid for details.
Definition: MOM_grid.F90:26
mom_diag_mediator::diag_ctrl
The following data type a list of diagnostic fields an their variants, as well as variables that cont...
Definition: MOM_diag_mediator.F90:240
mom_eos::calculate_density
Calculates density of sea water from T, S and P.
Definition: MOM_EOS.F90:68